引言

地震波在地下传播过程中，随着传播距离的增加，由于球面扩散和介质的吸收作用，会出现振幅衰减、视频率逐渐下降、高频相对低频衰减更快、主频下移、有效频带变窄的纵横向变化特征，且易使同条测线的单炮之间、同一单炮的道与道之间的能量存在明显差异，由此得到的成像剖面中深层能量弱、同相轴模糊、高频缺少、成像分辨率降低、难以识别精细构造和薄储层等。球面扩散效应会导致地震波振幅呈$ a/t $规律衰减，其中，a为球面扩散系数，$ t $为反射波旅行时间，均匀介质时$ a = 1/v $，层状介质时$ a = {v_1}/v_{{\rm{rms}}}^2 $，连续介质时$ a = {v_0}/v_{{\rm{rms}}}^2 $，$ v $为介质速度，$ v_{{\rm{rms}}}^{} $为均方根速度，$ v_1^{} $为层状介质第1层的速度，$ v_0^{} $为连续介质初始速度。大地吸收主要是由地震波在传播过程中地下介质颗粒间的内摩擦引起，介质的多相性、黏弹性等均会造成地震波振幅的复杂衰减，难以用一种机制进行描述（王云专等，1998；李胜强等，2020）。因此在实际地震资料处理中，需对地震波传播过程中由于球面扩散和大地吸收作用导致的振幅和频率衰减特征进行分析，并采用相应技术进行补偿，消除与反射界面反射系数无关的因素导致的振幅变化，使补偿后的振幅仅与反射系数有关，从而达到利用补偿后的振幅相对变化反映反射界面反射系数变化的目的（田媛媛等，2020）。

球面扩散和大地吸收补偿方法较多。在时频域方面，Newman（1973）首先较系统地分析和推导了均匀介质、层状介质和连续介质模型的球面扩散补偿公式。吕牛顿（1986）推导了波速随深度线性增加的连续介质中水平地质界面反射波的球面扩散补偿因子。高军等（1996）基于小波变换理论，提出了时频域球面扩散和吸收补偿的处理方法。王云专等（1998）通过分时窗、分频带处理的方式拟合得到球面扩散及频率吸收衰减曲线，并进行分频补偿。凌云等（2005）在原有时频域球面扩散与吸收衰减补偿方法的基础上，进一步提出了空间域补偿方法，通过对某油田三维地震数据的处理、解释和开发钻井的验证，表明该补偿方法可在消除近地表对地震属性空间影响的同时，满足相对保持储层振幅信息的处理要求。羊屋三维处理、解释一体化方法研究组（2002）采用时频域球面扩散和吸收补偿加地表一致性反褶积获得了较好的振幅和频率补偿结果，并实现了相对保持振幅的要求，研究发现了1块具有第二类AVO特性、常规成像条件下无法识别的地质异常体。刘财等（2005）研究了层状倾斜均匀介质地层中地震波球面扩散特征，并给出了相应球面扩散补偿公式。李鲲鹏等（2000）针对常规方法存在的不足，在地层吸收补偿中使用了小波包分解理论。刘喜武等（2006）使用了广义Ｓ变换，并与短时傅里叶变换和连续小波变换的应用效果进行了比较。张固澜等（2010）使用了改进的广义Ｓ变换，该变换具有兼顾频率分辨率和时间分辨率的特点。孙佳林等（2012）使用了曲波变换，对地震数据进行分频分角度补偿分析。杨学亭等（2014）提出了基于连续小波变换的时频域地震波能量衰减补偿方法，提高了地震资料的分辨率。杨存等（2013）、苏世龙等（2015）分别对时频域地震波能量衰减补偿方法进行了应用研究。张龙等（2015）将时频域球面扩散与吸收补偿技术应用于水陆交互区三维地震资料优化处理中，补偿后的资料接近陆地资料。田媛媛等（2020）将改进的地震波几何扩散补偿方法在表层黄土塬介质条件下进行应用，取得了较好的效果。

针对地震波衰减，地球物理学家们还提出了其他的补偿方法，如Hale（1982）提出了Hale算法，但该方法存在相位不相容性问题；Varela等（1993）对Hale算法进行了改进，解决了其不相容性问题；Hargreaves等（1987）提出了相移反Q滤波法；赵建勋等（1992）研究了串联反Q滤波方法，在实际地震数据处理中使用变Q，提高了结果的准确性；赵圣亮等（1994）提出了频域高频能量补偿方法；Wang（2002）提出了基于分层常Q模型的稳定反Q滤波方法，之后Wang（2006）又基于Q值随旅行时间变化的模型对该方法进行了改进；Ferber（2005）利用滤波组方法进行地震波能量的补偿，其中Q值可随时间任意变换；王小杰等（2019）将基于小波分频的Q值补偿技术应用于东海中深层油气勘探，有效补偿了东海陆架盆地中深层地震资料的振幅和相位，提高了中深层地震资料的质量；李胜强等（2020）总结了当前高分辨率深反射地震探测采集处理中的主要振幅补偿技术。

本文基于已有研究成果，系统梳理了时频域球面扩散和吸收补偿方法，研发了在时间、频率、炮检距和炮域实现的最小二乘地震波4D球面扩散与吸收衰减补偿方法。

1.   地震波衰减的物理机制分析

地震波在层状介质中传播时，总的衰减模式为（Newman，1973）：

式中，$ t $为传播时间；$ f $为频率；$ d(f,t) $为$ t $时刻的振幅，$ d(f,0) $为入射波振幅；$ c $为吸收系数，量纲为$ {m^{ - 1}} $，与介质的固有弹性参数有关，且与频率成正比$ c(f) $；$ v_1^{} $为层状介质第1层的速度；$v_{{\rm{rms}}}^{}$为均方根速度；$ {p_j} $为第$ j $个界面的反射系数，j为反射界面顺序号。

由式（1）可知，地震波在地下传播过程中，其衰减与传播距离和频率等因素有关，振幅总体上呈指数规律衰减，该衰减逐点时变，且频率不同，衰减特征不同。

研究认为，球面扩散和大地吸收主要引起地震波能量在时域、频域、炮检距域和炮域4个域内的衰减变化。因此，只有在上述4个域内对地震波衰减规律进行四维统计分析，才能有效合理地补偿球面扩散和大地吸收导致的地震波能量衰减，从而有效提高地震勘探的成像质量和分辨率。

2.   最小二乘4D球面扩散与吸收补偿方法原理

2.1   基于正弦函数的分频滤波

根据实际地震数据特点、频率分解和波形重构精度要求，选用正弦函数构建分频滤波函数，如图1所示，该函数可有效地将地震数据分解到一系列窄频带，同时具有较小的吉布斯效应和频带间的频率泄漏，具有较高的重构精度，可满足实际地震数据的处理要求，这是本方法提高地震波衰减补偿效果和成像分辨率的关键技术之一。

图 1  正弦分频滤波函数示意

Figure 1.  Scheme of cosine frequency division filter function

下载: 全尺寸图片 幻灯片

基于正弦函数的分频滤波表达式为：

式中，$ d(t) $为输入的地震数据；$ \omega (f,t) $为正弦分频滤波函数；$ D(f,t) $为分频后的数据。

采用自由分频策略，即可根据实际地震资料的频谱和分频扫描结果，自由灵活的设计分频策略，将地震数据分解为一系列的频带，各频带尺度可不同。该策略有利于准确定位需要相对加强或压制的频带，尽可能恢复有效信号能量。

2.2   最小二乘高阶e指数衰减曲线拟合与振幅补偿

2.2.1   最小二乘高阶e指数衰减曲线拟合方法

用于拟合地震波衰减曲线的函数应包含球面扩散与大地吸收衰减特性，基于式（1），采用高阶e指数函数作为球面扩散与大地吸收衰减的拟合函数。由于地震波不同频率成分的吸收衰减特点不同，故在不同频段分别进行拟合（凌云等，2005）。

分频策略下，仅考虑地震道时域振幅衰减时，采用的振幅衰减高阶e指数拟合函数为：

式中，$ D(f,t) $为衰减后的振幅，$ D(f,0) $为初始振幅；$ f $为频率；$ t $为旅行时间；$ c^k(f) $为复杂衰减曲线的不同阶数的拟合系数，其中$ k $为阶数。

对每个频带而言，式（3）e指数多项式整体可分为2部分，其零阶和一阶项可表示常规球面扩散，其高阶项可表示更复杂的球面扩散和大地吸收变化，该函数可更科学合理的描述各频带的球面扩散和大地吸收特征，以此为基础可实现随时间和频率变化的球面扩散和大地吸收衰减补偿，可较好的解决常规球面扩散和大地吸收补偿方法仅能补偿振幅随时间的衰减问题，及常规采用区域速度和常数吸收系数导致的工区内不同点存在剩余球面扩散和大地吸收问题。

在正弦分频策略下，针对每个频带地震数据，将地震道划分为$ N $个计算时窗，每个时窗内有M个样点，则各时窗的均方根振幅为：

式中，$ D(f,i) $为第$ i $个时窗的均方根振幅，$ i = 1,2,3, \cdots ,N $；$ d(f,{t_j}) $为该时窗内第$ j $个数据样点的振幅，$ j = 1,2,3, \cdots ,M $。

第$ i $个时窗的归一化均方根振幅为：

$ D_{}'(f,i) $构成的曲线可认为是实际地震波衰减曲线，由式（1）可知，其总体上呈e指数规律衰减。

设用于拟合第$ f $频带的振幅吸收衰减曲线的高阶e指数函数为：

采用最小二乘法通过${{p}}(f,t)$拟合实际地震波衰减曲线，构建目标函数为：

式中，$ {t_N} $表示各时窗中间位置的时间。

通常$ N \gg k $，在最小二乘意义下通过迭代求解，可估算$ c_{}^0(f) $，$ c_{}^1(f) $，···，$ c_{}^k(f) $等拟合系数，由此可拟合复杂地震波衰减曲线。

2.2.2   分频模型炮集数据的吸收衰减函数拟合与振幅补偿

在地表地震地质条件、地下结构相对稳定的地区，可认为区域内的大地吸收衰减特征是宏观稳定的，则大地吸收衰减的空间高频变化可认为是由近地表变化和地震波激发接收条件及耦合变化等引起的，可通过模型炮数据和每炮数据间统计振幅补偿曲线的空间高频差异对近地表等因素引起的空间能量变化进行补偿。模型炮数据选用激发和接收能量较好、干扰小、信噪比高、质量好、有代表性的炮集数据。首先拟合各频带模型炮集数据的吸收衰减函数：

对所有输入炮集的相应频带数据进行振幅补偿：

式中，$ D_{}'(f,t) $为振幅补偿后的结果。

2.2.3   分频输入炮集数据的高频吸收衰减函数拟合与振幅补偿

以$ D_{}'(f,t) $为新的输入炮集数据，使用式（6）拟合各频带输入炮集数据的高频吸收衰减函数${{p}}(f,t)$，对相应的频带进行补偿：

在上述补偿完成后，进一步对输入炮集的各频带进行统计，求取随炮检距变化的补偿因子，消除因炮检距不同而引起的振幅吸收变化，从而实现各频带输入地震数据的时间、频率、炮检距和炮域4D球面扩散与大地吸收衰减补偿。

2.3   波形数据重构

对于每个炮集，将补偿后的各频带数据进行波形重构，形成最终成果数据，即：

式中，$ d'(t) $为最终补偿输出的成果数据。

经上述处理后完成了4D球面扩散和大地吸收补偿。

相关研究表明，在地表地震地质条件和地下结构相对稳定的地区，球面扩散和大地吸收衰减特征具有空间宏观稳定性，而近地表等因素引起了激发接收能量的高频变化；球面扩散和大地吸收衰减引起的振幅变化是时间、频率、炮检距和炮位置的函数，且其影响程度远大于地下储层等反射系数的空间变化。因此，在陆上开展地震勘探须对球面扩散和大地吸收因素引起的振幅衰减进行有效补偿，对于岩性勘探而言，补偿方法还需满足相对保幅的处理要求（凌云等，2005）。本方法在时间、频率、炮检距和炮域内统计地震波振幅衰减特征曲线，并进行自适应针对性补偿，原本相对低频衰减更快的高频信号可得到更多补偿，同时具有减弱部分低频干扰的效果，从而实现补偿球面扩散和大地吸收因素引起的振幅在时间、频率和空间上变化的目的。

但该方法与其他方法一样，仅能在时间、频率、炮检距和炮域内补偿振幅衰减，无法补偿近地表因素引起的地震波波形空间变化，这个问题可通过高分辨率宽带Ricker子波多道反褶积等地震资料处理技术予以解决。

3.   实际处理效果分析

本文以西南某县城活动断层调查项目中某地震测线单炮数据为例进行分析。目标区地处四川省西南部，横断山系高山峡谷区南段，属高山、中山深切割地貌，总体呈南低北高地势，山脉呈南北向的陡峭状。目标区内最低海拔1 750 m。Sp1480炮振幅补偿处理测试效果如图2所示。对于原始单炮，由于近地表存在较厚的松散砾石层，地震波能量衰减较严重，总体上呈浅层能量急剧衰减、中深层能量衰减减缓的特征，且炮集上存在面波等强线性干扰，干扰波发生频散，中深层信息无法识别。对于采用多频带拟合补偿方式，经地震波4D球面扩散与吸收补偿后的单炮，通过细分频带，分别进行振幅衰减曲线拟合计算和振幅补偿，中高频段相对低频段获得更多补偿，在视觉效果上出现了大量连续、细的实同相轴。对于采用全频带统一拟合补偿，经地震波4D球面扩散与吸收补偿后的单炮，不分频带统一按新方法进行补偿，同相轴连续性变好，整体波组特征保护较好，也出现了部分实同相轴，但新出现的同相轴数量和分辨率低于采用多频带拟合补偿方式得到的单炮。

图 2  某实际数据处理效果对比

Figure 2.  Comparison of actual data processing effect

下载: 全尺寸图片 幻灯片

经新方法处理后，地震炮集记录的面貌得到了较大改善，尤其是中深层的中高频成分得到了较好补偿，且低频噪声受到一定压制，总体上浅中深地层的能量较均衡，证明了新方法的有效性和低频压噪性。

对应于图2数据的频率-振幅谱如图3所示，浅层频率-振幅谱分析时窗为：道号79～120，时间216～410 ms；深层频率-振幅谱分析时窗为：道号79～120，时间700～1 000 ms。由图3可知，对于浅层和深层，本方法均有效补偿了中高频成分，其能量显著增强，使主频得到了提高，频带得到展宽，从而提高了地震资料的质量，有利于后期地震资料的处理和解释。

图 3  实际处理数据的频率-振幅谱

Figure 3.  Frequency-amplitude spectrums of actual data processed

下载: 全尺寸图片 幻灯片

对应于图2数据的振幅变化曲线如图4所示，分析时窗为：道号72～75，时间132～1 100 ms。由图4可知，本方法的振幅补偿曲线更理想，各时间振幅曲线较稳定，其中商业软件振幅补偿后单炮振幅曲线波动变化大，且有部分时间补偿不足、部分时间补偿过大同时存在的问题，这说明新方法更好的拟合了实际能量衰减趋势。

图 4  实际处理数据的振幅变化曲线

Figure 4.  Amplitude curves of actual data processed

下载: 全尺寸图片 幻灯片

振幅补偿前、后的叠加剖面如图5所示，为便于对比，其他处理参数不变。对于原始叠加剖面，浅层能量较强，随着时间的增加，能量快速衰减，500 ms以下能量较弱，无法识别有效信息。对于经某商业软件振幅补偿后的叠加剖面，道集上残留的噪声严重影响叠加剖面质量，浅层大量有效信号被噪声污染，无法成像，如150 ms附近的主要反射同相轴被部分掩盖，无法有效识别和连续追踪，更无法识别主要反射层间的薄层反射。

图 5  振幅补偿前、后的叠加剖面

Figure 5.  Stack sections before/after amplitude compensation

下载: 全尺寸图片 幻灯片

由图5可知，新方法得到的叠加剖面整体面貌改善较大，浅中深层能量较均衡，有效中高频信息得到明显加强，低频噪声和强能量干扰得到一定程度压制，信噪比明显提高，复合波明显分开，地层接触关系清晰，有效同相轴清晰连续可追踪，层间薄层弱反射信号得到加强，分辨率得到大幅度提高，波组特征明显，相对保幅性较好。

对应于图5数据的频率-振幅谱如图6所示，浅层频率-振幅谱分析时窗为：CDP号338～700，时间121～320 ms；深层频率-振幅谱分析时窗为：CDP号338～700，时间700～1 000 ms。由图6可知，本方法处理后，对于浅层和深层，中高频成分均得到了较好补偿，低频噪声被相对压制，主频上移，频带被大幅展宽，效果较明显，验证了新方法的合理性。

图 6  振幅补偿前、后叠加剖面的频率-振幅谱

Figure 6.  Frequency-amplitude spectrums of stack sections before/after amplitude compensation

下载: 全尺寸图片 幻灯片

4.   结论

（1）实现了时间、频率、炮检距和炮域内地震波4D球面扩散和吸收补偿，解决了常规振幅补偿仅能补偿振幅随距离衰减而不能补偿振幅随频率衰减的问题。

（2）通过最小二乘高阶e指数曲线拟合求取各频带能量吸收衰减曲线，可更精确地对球面扩散和大地吸收进行补偿，并较好地保持了振幅的相对关系，且能够自适应求取补偿因子，减小了人为因素造成的影响，可适应更多地区的不同能量衰减模式。

（3）本方法采用了正弦函数分频，可实现对频率段的较细划分，且有较高的重构精度、较低的漏频和吉布斯效应，提高了振幅补偿的精确性和稳定性。

（4）实际处理结果表明，本方法可较好地补偿由于球面扩散和大地吸收因素造成的中深层中高频成分的损失，使主频上移，频带展宽，且能相对压制低频噪声，从而有效提高了地震资料的分辨率，为后续高分辨率处理和反演创造了良好条件。本方法有望在浅层地震勘探、发震构造调查、薄层油气或煤炭资源勘探中发挥作用。

（5）实际处理中应注意，当地震记录上存在较强的面波等噪声干扰时，会影响能量吸收衰减曲线求取的准确性，应先进行压制面波等去噪处理，然后再应用本方法进行处理。