0 引言
由于地层的非完全弹性性质，地震波在传播过程中会有一定的能量衰减和频散现象，对地震资料的分辨率有较大影响。反Q滤波(Q为品质因子)是常用的地层吸收衰减补偿方法，它最早由Hale[1]根据Futterman[2]的数学模型提出，但该方法计算量太大。Bickel 等[3]提出的褶积反Q滤波方法也存在计算量大的问题。Hargreaves等[4]提出了类Stolt偏移的相位校正方法，进一步减小了运算量。裴江云和何樵登[5]给出了可应用于低信噪比条件的最佳意义反Q滤波。Wang[6]对Hargreaves方法进行了改进，提出了一种更高效更稳定的常Q层反Q滤波方法。刘财等[7]研究发现了品质因子与吸收系数的关系,并由此提出频域衰减补偿方法。Wang[8]通过在补偿算子中引入稳定因子提出了一种反Q滤波稳定算法。
品质因子Q的准确估计是反Q滤波成功与否的关键，直接影响了补偿效果的好坏。Gladwin等[9]提出上升时间原理法。Jannsen等[10]提出子波模拟法。王辉等[11]基于上升时间原理，提出时间域相邻道地震波衰减成像方法，消除了震源能量对Q值估计的影响。然而由于实际地震资料的信噪比较低，导致时域方法估计Q值精度不高。频谱比值法是一种频域估计Q值的方法，它利用频率和振幅谱比值对数关系的斜率值来计算Q值[12, 13]。频谱比值法在微噪声的情况下效果更好[14, 15]。
与频域分析相比，时频谱分解可以更细致地刻画信号的局部频谱特征。因为衰减因子与地震信号的时间和频率均有关，所以可以通过时频变换把地震信号转换到时频域进行衰减补偿。较早的时频变换方法是短时傅里叶变换(STFT)。Gabor等[16]提出能够分析非平稳信号的STFT。后来Morlet等[17]提出了连续小波变换(CWT)。该方法可以自适应调整窗口宽度，是一种多分辨率分析方法。但是小波变换中小波基的构造有许多方法，容易产生多解性；而且小波系数无法提供波形频率部分的精确估计值，尤其在高频段。Stockwell等[18]提出了基于Morlet小波变换的S变换方法，它兼备短时傅里叶变换和小波变换的优点，但S变换的分辨率比较低。Liu等[19]提出了局部时频变换的方法，用傅里叶级数匹配目标信号并利用整形正则化最小二乘条件求解方程。该方法能够调节选取谱分解的频率范围和频率采样间隔，因此比S变换更加灵活并具有与CWT相近的分辨率。由于该变换在最小二乘意义下可逆，因此适合于处理非平稳信号的反Q补偿。
在常见的地层Q模型中，Kolsky衰减模型因为其简单方便被广泛应用[20]。本文将地下介质分为不同Q值的若干层，利用局部时频变换把信号转为时频谱后，通过频谱比值法求取各层Q值，最后利用基于Kolsky衰减模型的大地滤波算子在时频域实现信号的衰减补偿。
1 理论基础 1.1 局部时频变换基本原理局部时频变换的基本原理是利用正则化条件约束最小二乘问题以解决数学欠定问题。它是一种能够表征非平稳地震信号时变频率特征的方法。该方法的关键思想是使输入数据和其所有谐波分量间的误差最小化，通过控制时间分辨率来进行正则化约束，实现非平稳回归分析。
傅里叶级数定义为关于正弦和余弦函数总和的函数。假设一个物理可实现信号f(x)的周期为L，则其傅里叶级数为
其中：n=0,±1,±2,±3,±4,…。傅里叶级数的复数表达式为
其中：An是傅里叶系数；Ψn(x)=ei(2πnx/L)。时变系数An(x)的绝对值为信号的时频谱，方程(2)提供了时频转换的反运算。该时频变换的频率范围和频率采样间隔可以通过Nyquist频率来决定或者由使用者自行分配。
非平稳回归分析允许系数An随x变化，可以在最小二乘准则下解如下目标函数以求得系数An(x)：
因为该最小二乘问题(式(3))的未知数要多于限制条件，所以该最小化问题是欠定的。解决办法是添加额外的正则化条件，以限制估计系数可能的变化。经典的约束条件有Tikhonov正则化，该正则化方法能够使目标函数满足局部光滑。整形正则化[21]是另一种可选的约束条件，该方法具有比经典Tikhonov正则化更少的条件数、更快的收敛速度和更少的调节参数，并且参数选取比较灵活，因此具有更好的适用性。局部时频变换即是利用整形正则化来求解欠定问题的时频转换方法。
如图 1所示，局部时频变换能够准确地在时频谱中将两个具有不同起始频率和调频系数的线性信号的合成信号分解开来，直观反映了信号的时频谱特征。
 |
| a. chip信号; b. chip信号的局部时频谱。图 1 chip信号及其时频谱 Fig. 1 Chip signal and its time-frequency spectrum |
| |
对于复杂的非平稳信号，局部时频变换依然有着很好的处理效果。图 2为对非平稳信号做局部时频分解和重构,并计算重构后信号与原信号的误差。通过对比误差和原信号的大小，我们可以发现误差与原信号的比值不足百分之一，误差几乎可以忽略不计。因此，局部时频变换具有精度较高的重构能力。
 |
| a. 非平稳信号；b. 非平稳信号局部时频谱；c. 重构信号；d. 误差分析。 图 2 非平稳信号时频谱及其信号重构 Fig. 2 Time-frequency spectrum and reconstruction of the non-stationary signal |
| |
首先假设衰减因子λ和频率f遵循严格线性关系：
式中：v是相速度；Q是品质因子。此时定义相速度如下： 式中，vr和fr分别表示参考相速度和参考频率。参考频率一般选为频谱中部的频率值。 随后Wang[22]研究认为fr影响频散校正的准确度，将fr设为一个较大值能有效恢复频散效应。此时波数k变为如下复数形式：
将公式(6)代入公式(7)，复波数k(f)变为
在频率域求解平面波得到频率域内原信号S0(f)与衰减信号S(f)的方程式为
式中,ζ表示波前面沿射线路径从震源到检波器距离。将式(8)代入式(9)并忽略
(因为Q值一般为20~200或更大)，可以得到:
式中，t是地震波传播时间。定义大地滤波算子G：
公式(11)右侧的e指数实部表示振幅衰减部分指数，虚部表示频散和相位变化部分指数。把G替换到式(10)中并对所有的频率进行叠加，得到时间域内的衰减地震道为
式(12)便是基于Kolsky模型的时间域衰减计算公式。根据衰减模型得到补偿因子为
所以时频域地震波衰减补偿公式为
1.3 频谱比值法求Q值
地震信号在均匀吸收介质中的振幅谱如下：
其中：z是地震信号传播距离；A(t)是与频率有关的变量；S(f,t0)是初始频谱。当旅行时为t1和t2时,我们可以得到： (17),(16)两式相除得:
令d(f)等于式(18)的右边部分,则式(18)可以表示为
对式(19)求取对数得 式中，C是取决于反射系数的常量。式(20)所表示的直线斜率为 此时我们得到Q值为
在Q值的计算过程中，需要对谱比法求得的数组求斜率，本次研究采取最小二乘线性拟合方法。
2 理论模型测试本文采用Ricker子波与反射系数序列褶积后的复波模型作为理论模型进行测试。Richer子波主频为50 Hz，时间采样间隔为0.001 s，最大延续时间为2 s。
根据实际地层衰减规律，首先建立变Q地层模型。选取不同时深段三层品质因子的地层模型：地震信号传播时间在0~1 s段，Q值为120；在1~2 s段，Q值为160；在2~3 s段，Q值为200。在该地层模型的基础上建立原始复波信号(图 3a)。利用Kolsky衰减模型(式(12))对原始复波信号进行衰减，衰减后的信号如图 3b所示。复波信号理论时频谱如图 4a所示，分别使用S变换和局部时频变换得到的时频谱如图 4b、c所示。可以看出，局部时频变换的分辨率和准确度要明显优于S变换。抽取不同时深的频谱值(图 5a、b)，将其做商后取对数，使用最小二乘线性拟合得到近似线性的数据，如图 5c所示；把直线斜率代入式(22)求得0~1 s段Q值，其他段Q值同法可求。
 |
| a.原始复波信号；b.衰减后复波信号。 图 3 复波信号衰减前后对比图 Fig. 3 Comparison of complex wavelet signal and its attenuation |
| |
 |
| a.理论时频谱；b.S变换时频谱；c.局部时频谱。 图 4 复波信号时频谱 Fig. 4 Time-frequency spectrum of complex wavelet signal |
| |
 |
| a.0.2 s处频谱值；b.0.5 s处频谱值；c.最小二乘线性拟合后的对数域频谱比值。 图 5 频谱比值法示意图 Fig. 5 Spectrum ratio method |
| |
得到地层Q值后，我们分别对衰减后的复波信号做S变换衰减补偿和局部时频变换衰减补偿，结果如图 6a、b所示。可以发现，两种变换都使衰减信号的振幅得到补偿，畸变的波形得到恢复。比较两者补偿后的误差(图 6c、d)可以看出，局部时频变换的误差远小于S变换。
 |
| a.S变换补偿后的复波信号；b.局部时频变换补偿后的复波信号；c.S变换补偿信号与原信号差值；d.局部时频变换补偿信号与原信号差值。图 6 补偿后的复波模型及其补偿误差 Fig. 6 Compensational complex wavelet signal and the residual error |
| |
本文使用基于局部时频变换的衰减补偿方法进行某地区实际二维叠前资料处理。图 7a是原始地震剖面，共201道，时间长度为1.5 s，采样间隔为2 ms。根据频谱比值法计算的数据发现，将该地震资料划分为浅层和深层两段不同Q值地层进行补偿比较合适。对地震单道数据计算后，将该地震资料分为0~0.7 s、0.7~1.5 s两层不同Q值的地层。图 7b是利用本文方法补偿后地震剖面图。为了更好地观察补偿前后的差异，对地震剖面的第20~40道，0.2~0.4 s和0.8~1.0 s分别做补偿前后的频谱归一化对比，如图 7c、d所示。由图 7发现，补偿后地震波的同向轴变得更加清晰，地层吸收衰减得到补偿，其中尤以深层明显。实际资料处理结果表明，基于局部时频变换的衰减补偿方法能够有效补偿地震波衰减，提高地震资料的分辨率。
 |
| a.地震剖面；b.补偿后地震剖面；c.浅层振幅谱归一化对比；d.深层振幅谱归一化对比。 图 7 某地区二维叠前地震资料及其局部时频变换域补偿后结果 Fig. 7 Compensation of seicsmic data using local time-frequency transform |
| |
局部时频变换域地震波吸收衰减补偿方法在理论模型测试和实际地震资料处理中都取得了较好的效果，从中我们得到了以下两点结论：
1)局部时频变换经理论模型测试证明其能很好地表征信号的时频特点，具有较高的重构精度，而且参数选取比较灵活，在最小二乘意义下可逆，因此非常适合处理非平稳信号的反Q滤波。
2)本文方法对中深层地震波吸收衰减补偿效果更好，有效补偿了中深部地层的能量衰减，提高了地震资料的分辨率。
| [1] | Hale D. An Inverse Q-Filter[R]. San Francisco:Stanford University, 1981:231-243. |
| [2] | Futterman W I. Dispersive Body Waves[J]. Geophy-sics Res, 1962, 67(4):5279-5291. |
| [3] | Bickel S H, Natarajan R R. Plane-Wave Q Deconvolution[J]. Geophysics, 1985, 50(9):1426-1439. |
| [4] | Hargreaves N D, Calvert A J. Inverse Q Filtering by Fourier Transform[J]. Geophysics, 1991, 56(4):519-527. |
| [5] | 裴江云, 何樵登. 基于Kjartansson模型的反Q滤波[J].地球物理学进展, 1994, 9(1):90-100. Pei Jiangyun, He Qiaodeng. Inverse Q Filtering According to Kjartansson Model[J]. Progress in Geophysics,1994, 9(1):90-100. |
| [6] | Wang Yanghua. A Stable and Efficient Approach of Inverse Q Filtering[J]. Geophysics, 2002, 67(2):657-663. |
| [7] | 刘财, 刘洋, 王典,等.一种频域吸收衰减补偿方法[J]. 石油物探, 2005, 44(2):116-118. Liu Cai, Liu Yang, Wang Dian, et al. A Method to Compensate Strata Absorption and Attenuation in Frequency Domain[J].Geophysical Prospecting for Petroleum, 2005, 44(2):116-118. |
| [8] | Wang Y H. Inverse Q-Filter for Seismic Resolution Enhancement[J]. Geophysics, 2006, 71(3):51-60. |
| [9] | Gladwin M T, Stacey F D. Anelastic Degradation of Acoustic Pulses in Rocks[J]. Phys Earth Plan Int, 1974, 8(4):332-336. |
| [10] | Jannsen D, Voss J, Theilen F. Comparison of Methods to Determine Q in Shallow Marine Sediments from Vertical Reflection Seismograms[J]. Geophysical Prospecting, 1985, 23(4), 479-497. |
| [11] | 王辉, 常旭, 刘尹克.时间域相邻道地震波衰减成像研究[J].地球物理学报,2001,44(3):396-403. Wang Hui, Chang Xu, Liu Yinke. Seismic Neighbo-ring Traces Attenuation Tomography in Time Domain[J]. Chinese Journal of Geophysics, 2001,44(3):396-403. |
| [12] | Bath M. Spectral Analysis in Geophysics[M]. New York:Elsevier,1974. |
| [13] | 冯晅, 张瑾, 刘财,等. 基于改进的Kolsky模型波场延拓公式的纵波Q值、横波Q值估计[J]. 吉林大学学报(地球科学版), 2014, 44(1):359-368. Feng Xuan, Zhang Jin, Liu Cai, et al. Estimation of P-and S-Wave Quality Factors Based on the Formula of the Wave-Field Continuation in Modified Kolsky Model[J]. Journal of Jilin University (Earth Science Edition), 2014, 44(1):359-368. |
| [14] | Tonn R. The Determination of the Seismic Quality Factor Q from VSP Data:A Comparison of Different Computational Methods[J]. Geophysical Prospecting, 1991, 39(1):1-27. |
| [15] | Tonn R. Comparison of Seven Methods for the Computation of Q[J]. Physics of the Earth and Planetary Interiors,1989,55(3):259-268. |
| [16] | Gabor D. Theory of Communication[J]. Journal of the Institution, 1946, 93(26):429-441. |
| [17] | Morlet J, Arens G, Fourgeau E, et al. Wave Propagation and Sampling Theory:Part I:Complex Signal and Scattering in Multilayered Media[J]. Geophy-sics, 1982, 47(2):203-221. |
| [18] | Stockwell R G, Mansinha L, Lowe R P. Localization of the Complex Spectrum:The S Transform[J]. IEEE Transactions on Signal Processing, 1996, 44(4):998-1001. |
| [19] | Liu Yang,Fomel Sergey. Seismic Data Analysis Using Local Time-Frequency Decomposition[J]. Geophysical Prospecting,2013,61(3):516-525. |
| [20] | Kolsky H. The Propagation of Stress Pulses in Viscoelastic Solids[J]. Philosophical Magazine,1956,1(8):693-710. |
| [21] | Fomel S. Shaping Regularization in Geophysical-Estimation Problems[J]. Geophysics, 2007, 72(2):R29-R36. |
| [22] | Wang Y H. Seismic Inverse Q Filtering[M]. New York:Wiley-Black Well,2008:1-248. |