周聪，刘江平，曾祥芝，范承余，曹进

(1.中国地质大学(武汉)地球物理与空间信息学院，湖北 武汉，430074；2.安徽省高速公路控股集团有限公司，安徽 合肥，230051)

在二维地震勘探中，因为受野外施工条件(如障碍物、禁区等)、成本、效益等多方面因素的限制，野外采集道间距往往过大，容易产生空间假频，从而影响地震资料的多道处理。而利用地震道插值能加密空间采样率，防止频散出现和提高信噪比[1]。目前，地震道内插的方法有很多。Spitz[2]提出在f-x域中进行地震道内插，该方法计算量较大且精度较低。此后，国九英等[3-5]在此基础上做了更深入的研究。Wang等[6-7]将小波变换理论引入地震道插值中，并指出小波插值方法是一种能够消除假频并扩展带宽的地震道插值方法。Trad等[8-9]采用预优共轭梯度法对大矩阵求逆，用双曲Radon变换对CMP道集数据进行插值处理，取得较好的插值效果。崔兴福等[10]也利用小波变换在时间和频率域良好的局部化性质，通过一维小波变换重构公式实现地震道内插。王维红等[11]基于带限正反最小平方抛物Radon变换的Levinson递推算法，给出利用抛物Radon变换进行地震道重建外推的基本原理以及叠前地震数据规则化地处理流程。Zhou等[12]提出在Wavelet-Radon域校正假频以及地震数据的插值，该方法可以非常可靠地恢复高频信号，在合成数据和实际数据的处理中都取得很好的效果。Herrmann等[13-14]利用Curvelet变换的稀疏特性，成功将其应用于地震数据的重建。刘国昌等[15]在Curvelet变换的基础上，提出基于POCS算法的Curvelet变换震数据插值方法。在以往的反射地震勘探中，面波往往被认为是一种具有明显高振幅和频率特性的规则干扰波而被压制[16]。而在多波多分量地震勘探中，由于横波静校正的需要，许多学者都利用瑞雷面波的频散信息反演浅层横波速度结构，进而求取静校正量[17-18]。面波的反演主要是基于频散曲线的反演，而对面波数据进行有效的插值可以提高频率速度谱的信噪比，有利于频散曲线的精确提取。在大时空采样地震记录中，由于面波的速度很低，其记录中的高频成分容易出现超周期现象，降低其频率速度谱的信噪比，而且容易出现假的高阶能量。而对面波插值，能够在一定程度上消除假频，提高其频率速度谱的信噪比。为此，本文作者利用面波的线性特征，首先对原始面波数据进行线性动校正(LNMO)处理，然后利用二维小波变换对其插值，最后反线性动校(RLNMO)恢复信号。理论和实际数据插值前后频率速度谱的变化表明：本文方法能有效提高频率速度谱的信噪比和频散曲线的提取精度。

1 小波插值的基本原理

以一维小波变换来说明小波插值的基本原理。根据小波变换的基本原理，信号在分辨率为2-m-1空间的近似部分可以由分辨率为2-m空间的近似部分和细节部分通过小波重构来得到[6]。如图1所示，将原始信号c0看成是分辨率为20的空间的近似部分，如果该空间的细节部分d0可以忽略或估计，那么分辨率为2-1空间的近似信号c-1就可以由原始信号得到。

图1 小波重构示意图Fig.1 Sketch of wavelet reconstruction

反变换可表示为：

在式(4)中，是已知的，而未知。可以通过以下2种方法来确定。

(1)方法1。由式(4)可知，若分辨率为20空间的细节部分可以忽略，即令，则可得：

在很多情况下，是可以忽略的。这时，信号通过插值来增加信息量是一个有效的处理过程。

(2)方法2。由于忽略相当于将信号进行简单的去噪处理，剔除该信号中的高频部分，因此，对于需要保留高频信息的信号来说，方法1可能达不到该要求。在图像插值里面，一般是利用原图像来确定高频信息[19-20]。具体做法是分解原图像得到高频信息，由于不同分辨率下各个对应子图间存在相似性，可以通过一般插值(如3次样条)得到和原始图像同一分辨率下的高频信号。王真理等[7]指出为得到的信息，可对原始数据进行空间方向高通滤波，或采用频率调制方法，提高其频率，并保持波形与原剖面的相似性。

由于利用方法2得到的会包含一些不真实的信息以及其实现的复杂性，因此，本文尝试采用方法1。在二维小波重构中，需要将3个细节部分置0。通过面波的频率速度谱来分析插值的效果，特别是其高频部分频散能量的变化。本文选取的母小波是db4小波。

2 模拟数据试验结果分析

下面对数值模拟的体波和面波记录做小波插值，分别从体波插值结果的误差以及面波插值前后频率速度谱的变化来说明本文插值方法的有效性。

2.1 体波插值

如图2所示地震记录为二维合成记录，道间距2.5 m，道数为480，时间采样率为0.5 ms，采样长度800 ms，雷克子波主频35 Hz，同相轴的倾角分别为1 ms/道和0.5 ms/道，对该记录做3次小波插值试验。将原始数据分别均匀抽道成5，10和20 m，然后插值成2.5 m，并与原始合成记录进行比较。由于小波插值是同时对时间和空间进行插值，而在地震插值中我们主要关心道间插值，所以，二维插值实验中舍弃时间方向上插值的数据，只考虑道间插值。

图2 二维合成地震记录Fig.2 2D synthetic seismic data

小波插值实验结果见图3～5。从图3～5可见：随着道间距和插值倍数的增加，小波插值的误差也越来越大。图3(b)显示的插值结果很好，仅仅是在同相轴相交处稍差，而图5(b)中倾角大的同相轴的插值结果要远差于倾角小的同相轴。所以，对于不考虑细节部分的小波插值方法来说，虽然其能够在一定程度上消除假频，但较严重的假频反过来又会影响插值的效果。

2.2 面波插值

如图6(a)所示为根据3层介质模拟出的面波记录，其道间距为2 m，时间采样率为0.5 ms，雷克子波主频为20 Hz。将其均匀抽道得到道间距4 m的试验信号，如图6(b)所示。利用小波变换直接对该记录插值，得到结果如图7(a)所示。由图7(a)可见：插值效果不理想，插值出来的中高频部分几乎是拷贝其左边道的同相轴信息。若同相轴斜率较大，则其高频部分容易出现假频信息，会影响插值效果。可以通过2种数据处理的方法来处理：一是通过坐标旋转，将数据转动某一角度使得大倾角同相轴的斜率变小，但此时在新坐标系中的纵横向采样率发生变化，难以处理且可行性有待验证；另一种是对于线性信号来说(如面波)，可以先对原始数据做类似线性动校正的处理，然后进行小波插值，最后反线性动校正恢复信号即可。后者即相当于单坐标轴的旋转。

图3 小波插值试验(1倍插值)Fig.3 Wavelet interpolation experiment (1× interpolation)

图4 小波插值试验(2倍插值)Fig.4 Wavelet interpolation experiment (2× interpolation)

图5 小波插值试验(3倍插值)Fig.5 Wavelet interpolation experiment (3× interpolation)

图6 理论面波记录Fig.6 Synthetic rayleigh waves

图7 小波变换对瑞雷面波插值Fig.7 Rayleigh wave interpolation using wavelet transform

下面对图6(b)所示的面波记录进行一定程度的线性动校正处理。由于要往上“拉伸”面波，因此，需要在面波记录的零时刻前端补零，以防止面波被拉出边界。经线性动校处理后的结果如图7(b)所示，再对该记录做小波插值，然后做反线性动校处理恢复信号，如图7(c)所示。将图6(a)，7(a)和7(c)中第36道抽出来进行对比，结果见图7(d)。由图7(d)可见：经过线性动校处理的插值效果要优于直接利用小波插值，特别是在面波记录的后半部，如0.4～0.6 s之间。从整个记录剖面来看，除了记录2端的几道，结合线性动校处理的小波插值能够很好地恢复面波信号。对小波插值方法来说，记录2端道的插值效果较差是其缺陷，本文作者认为在小波插值过程中，默认记录以外的部分为零值，所以，插值出来的边界道能量要弱。该处理方法通过单坐标轴的旋转，可减小线性同相轴的斜率，避开直接对大时空采样地震记录进行小波插值，对于同相轴为近似线性的面波记录的道内插有较好的效果。该处理过程的重点部分是线性动校正量的估算。以面波为例，面波在地震记录上呈“扫帚”状，可以大致画出“扫帚”的上、下界限，即可以估算出面波的最大速度和最小速度；然后可以取平均速度作为参考速度；最后，利用道间距计算出校正量即可。

2.3 面波的频率速度谱的变化

在转换波静校正中，地表横波速度模型的获取尤其关键，而面波速度非常趋近于横波，因此，通过面波反演来获取横波速度是可行的。面波的反演主要是基于频散曲线的反演，而对面波数据进行插值可以提高频率速度谱的信噪比，有利于频散曲线的提取。

下面本文从插值前后面波频率速度谱的变化来分析插值的结果。提取频散曲线可以通过2个线性变换实现[21-22]：首先对原始数据做倾斜叠加(即线性Radon变换)，将数据变为慢度-截距域(p-τ)；其次，沿截距τ方向做一维傅里叶变换得到慢度-频率域(p-f)；最后，将慢度转换成速度即可获得速度-频率域(v-f)。

如图8(a)～(c)所示分别是图6(a)，6(b)和7(b)所示面波通过线性Radon变换得到的频率速度谱。对比图8(a)和8(b)可见：均匀抽道后，在能谱的高频高速区域出现一个较强的假高阶能量。这可能是由于抽道后，面波的高频部分产生超周期现象，其高频信息的同相轴斜率由小变大，因此，其高频信息的速度增大，出现了如图8(b)中箭头A所示的假高阶能量。在箭头B所示的低速区域，也产生干扰信息。对面波插值后，如图8(c)所示，箭头A处的假高阶能量团明显减弱，而箭头B处的噪音干扰也得到改善。从插值前后频率速度谱的变化情况说明本文的插值方法可以对面波的高频信息进行有效插值，并能够提高其频率速度谱的信噪比。

3 实例分析

下面将本文的小波插值方法应用于实际数据。同样，根据记录插值前后面波频率速度谱的变化来验证插值方法的有效性。

3.1 野外地震记录插值

图8 理论面波数据的频率速度谱Fig.8 Dispersive images in f-v domain of synthetic data

如图9(a)所示为野外采集的原始地震记录，时间采样率为2 ms，记录长度为7 s，道间距为20 m，记录道数36道(起始道为零道)。图9(a)中上部黑色部分为体波。图9(b)所示为将原始地震记录均匀抽道后得到的地震记录，道间距为40 m。利用本文的插值方法对该记录进行插值，得到图9(c)所示的插值结果。对比图9(a)和9(c)可以看到：面波插值效果较好，得到的面波记录与原始记录中的面波基本一致，没有带入背景噪声。将图9(a)和9(c)中第16道抽出来对比，得到图9(d)。可以看到面波主体部分的波形匹配较好，只是在面波的尾部如时间轴4.3 s附近重构效果稍差。从图9(a)，(c)和(d)均可以发现体波部分发生较大变化。这是由于在插值前进行线性动校处理，体波同相轴多是非线性的，在“拉平”面波的同时，将体波的同相轴破坏。因此，本文的插值方法无法同时对体波和面波进行有效插值。由于体波速度较高，其假频现象要比面波记录中的少。当需要对体波进行有效插值时，可以直接利用小波变换对体波插值。

图9 实际地震数据Fig.9 Real seismic data

3.2 插值前后频散曲线的变化

如图10(a)～(c)所示分别是图9(a)～(c)通过线性Radon变换得到的频率-速度谱。对比图10(a)和(b)可以看到：抽道前后其频率速度谱发生明显变化，抽道后地震记录的频率速度谱的信噪比降低很多，在高频高速的区域出现了3个假能量团，而且2个能量特别强，这有可能是面波高频率信息超周期的现象造成的；在高频低速的区域也出现了较多的噪音干扰；其二阶的频散能量发生明显弯曲，尾部的能量团也产生断裂，对面波插值后，高频高速区域的3个假能量团均已消失，说明面波中高频部分的假频信息得到改善；高频低速区域的噪音干扰也已消除，二阶的频散能量团较插值前更圆滑，而且其尾部断裂的能量团也得到改善。对比图10(a)和(c)可见：除了三阶频散能量团没有较好地恢复外，其他的与原始记录的频率速度谱基本相同。这可能是由于抽道后该部分信息已经缺失，而插值无法恢复缺失的信息，只能突出已有的有效信息，压制干扰信息。因此，实际资料的频率速度谱的变化也说明本文插值方法能够较好地应用于实际，提高频散-速度谱的信噪比，有利于精确提取面波频散曲线。

图10 实际地震资料的频率速度谱Fig.10 Dispersive images in f-v domain of real seismic data

4 结论

(1)基于小波变换的基本原理，采用db4母小波，实现了对瑞雷面波的时空插值。

(2)根据瑞雷面波低速和线性特点，将线性动校正引入小波插值中，有效实现了对瑞雷面波的插值，改善了瑞雷面波的插值效果。

(3)在大时空采样情况下，对瑞雷面波进行插值，可有效提高频率速度谱的信噪比，有利于瑞雷面波频散曲线的精确提取。

[1]罗丹.f-k域地震道插值方法研究[D].成都: 成都理工大学地球物理学院, 2009: 1-3.LUO Dan.The research of methods to seismic trace interpolation in f-k domain[D].Chengdu: Chengdu University of Technology.College of Geophysics, 2009: 1-3.

[2]Spitz S.Seismic trace interpolation in the f-x domain[J].Geophysics, 1991, 56(6): 785-794.

[3]国九英, 周兴元, 俞寿朋.f-x域等道距道内插[J].石油地球物理勘探, 1996, 31(1): 28-34.GUO Jiu-ying, ZHOU Xing-yuan, YU Shou-peng.Iso-interval trace interpolation in f-x domain[J].Oil Geophysical Prospecting,1996, 31(1): 28-34.

[4]Porsani M J.Seismic trace interpolation using half-step prediction filters[J].Geophysics, 1999, 64(5): 1461-1467.

[5]WANG Yang-hua.Seismic trace interpolation in the f-x-y domain[J].Geophysics, 2002, 67(4): 1232-1239.

[6]WANG Zhen-li, LI Yan-da.Trace interpolation using wavelet transform[C]//64th SEG Annual International Meeting, Los Angles, USA: SEG, 1994, 13: 729-730.

[7]王真理, 李衍达.基于小波的信号部分重构与插值[J].清华大学学报: 自然科学版, 1995, 35(5): 67-72.WANG Zhen-li, LI Yan-da.Signal partial reconstruction and interpolation based on wavelet transform[J].Journal of Tsinghua University: Science and Technology, 1995, 35(5): 67-72.

[8]Trad D O, Ulrych T J, Sacchi M D.Accurate interpolation with high-resolution time-variant Radon transforms[J].Geophysics,2002, 67(2): 644-656.

[9]Trad D O.Interpolation and multiple attenuation with migration operators[J].Geophysics, 2003, 68(6): 2043-2054.

[10]崔兴福, 刘东奇, 张关泉.小波变换实现地震道插值[J].石油地球物理勘探, 2003, 38(增刊): 93-97.CUI Xing-fu, LIU Dong-qi, ZHANG Guan-quan.Seismic trace interpolation using wavelet transform[J].Oil Geophysical Prospecting, 2003, 38(Suppl): 93-97.

[11]王维红, 刘洪.抛物Radon变换法近偏移距波场外推[J].地球物理学进展, 2005, 20(2): 289-293.WANG Wei-hong, LIU Hong.Near offset wavefield extrapolation based on parabolic radon transform[J].Progress in Geophysics, 2005, 20(2): 289-293.

[12]ZHOU Yu, Ferguson J, McMechan G, et al.Wavelet-Radon domain dealiasing and interpolation of seismicdata[J].Geophysics, 2007, 72(2): 41-49.

[13]Herrmann F J, Hennenfent G.Non-parametric seismic data recovery with curvelet frames[J].Geophys J Int, 2008, 148(1):233-248.

[14]Naghizadeh M, Sacchi D.Beyond alias hierarchical scale curvelet interpolation of regularly and irregularly sampled seismic data[J].Geophysics, 2010, 75(6): WB189-WB202.

[15]刘国昌, 陈小宏, 郭志峰, 等.基于Curvelet变换的缺失地震数据插值方法[J].石油地球物理勘探, 2011, 46(2): 237-246.LIU Guo-chang, CHEN Xiao-hong, GUO Zhi-feng, et al.Missing seismic data rebuilding by interpolation based on Curvelet transform[J].Oil Geophysical Prospecting, 2011, 46(2):237-246.

[16]张恒磊, 刘天佑, 李红巧.Curvelet域面波衰减方法研究[J].中南大学学报: 自然科学版, 2011, 42(8): 2372-2378.ZHANG Heng-lei, LIU Tian-you, LI Hong-qiao.Attenuation of surface wave in Curvelet domain[J].Journal of Central South University: Science and Technology, 2011, 42(8): 2372-2378.

[17]孟小红, 郭良辉.利用地震瑞利波速度反演求取P-SV波横波静校正量[J].石油地球物理勘探, 2008, 42(4): 448-453.MENG Xiao-hong, GUO Liang-hui.Using velocity inversion of seismic Rayleigh wave to compute S-wave statics of P-SV wave[J].Oil Geophysical Prospecting, 2007, 42(4): 448-453.

[18]Bansal R, Ross W, Lee S, et al.A novel approach to estimating near-surface S-wave velocity and converted-wave receiver statics[C]//79th SEG Annual International Meeting, Houston,USA: SEG, 2009, 28: 1192-1196.

[19]卢珏.基于小波的图像插值研究[J].武汉理工大学学报, 2003,25(1): 81-83.LU Jue.Research of image interpolation based on wavelet[J].Journal of Wuhan University of Technology, 2003, 25(1): 81-83.

[20]陶洪久, 柳健, 田金文.基于小波变换和插值的超分辨率图像处理算法[J].武汉理工大学学报, 2002, 24(8): 63-66.TAO Hong-jiu, LIU Jian, TIAN Jin-wen.Surper resolution image processing algorithm based on wavelet transform and bilinear interpolation[J].Journal of Wuhan University of Technology, 2002, 24(8): 63-66.

[21]George A M, Mathew J Y.Analysis of dispersive waves by wave field transformation[J].Geophysics, 1981, 46(6): 869-874.

[22]LUO Yin-he, XIA Jiang-hai, Richard D M, et al.Rayleigh-wave dispersive energy imaging using a high-resolution linear radon transform[J].Pure and Applied Geophysics, 2008, 165(5):903-922.