2017, Vol. 32
反褶积是地震资料处理的一个必须部分 (Robinson and Treitel, 1967; 伊尔马滋, 2006)，生产中常用的是基于Wiener最小平方反褶积技术.根据褶积理论Wiener反褶积是从已知的地震道中估算地层的反射系数，但是震源子波也是未知变量，为了求得反射系数需要两个假设条件：(1) 震源子波必须是最小相位的，这意味着震源子波的振幅 (等价于其自相关) 可以通过解地震道的能量谱的平方根获得，然后对震源子波的振幅的自然对数求Hilbert变换可求得震源子波的相位；(2) 反射系数是统计意义上的白噪的、平稳的 (其振幅谱大致是一常数).常规的Wiener反褶积假设震源子波是平稳的即震源子波在其旅行过程中是一成不变的，这与实际地层反射相差较大，至少有两个显而易见的物理影响没有考虑：地层黏滞性衰减 (Kjartansson, 1979) 和short-path多次波 (Margrave et al., 2011).常数Q理论 (Kjartansson, 1979) 描述地层黏滞性衰减，其数值呈现走时和频率的指数衰减.
地震道的非平稳性的解法大致分为两类 (Margrave et al., 2011)：直接的时变反褶积和反Q滤波法，Clarke (1968)提出一种基于最优的Wiener滤波的时域非平稳的反褶积，Koehler和Taner (1985)提出一种广义的时变反褶积的数学理论.Bickel和Natarajan (1985)采用反Q滤波作为在平面波域的一种反褶积类型，所有的反Q滤波方法都面临一个共同的难题—反滤波的稳定性 (Wang, 2002)，因为Q滤波引起走时和频率的指数衰减，所有反Q滤波肯定导致走时和频率的指数增长并因此产生不稳定.Wang (2002)提出了一种基于波场延拓理论的反Q滤波方法，可同时进行振幅补偿和相位校正，然后Wang (2006)又进一步在Gabor域将实现这种稳定算法，即推广到Q模型随深度或走时连续变化的情形，其计算结果更准确，在此基础上陈增保等 (2014)提出一种在Gabor域实现的带限稳定反Q滤波算法，它适用于一般Q模型，即根据地震道的平均Q值和信噪比差异设计不同的时变带通滤波器.
针对地震子波的时变的非平稳性特征，彭才等 (2007)提出了基于动态褶积模型的子波估计方法，高静怀等 (2009)利用基于反射地震记录变子波模型提取地震子波进而提高地震记录，这些方法大都对地震记录采用分段处理，即首先在时间上进行分段，假设每段内的地震子波是时不变的，进而在每段内提取出一个时不变的地震子波，张漫漫等 (2014)提出一种基于时频域谱模拟的时变子波提取方法，另外郭廷超等 (2015)在传统谱模拟反褶积的基础上，提出了一种基于S变换的时变谱模拟反褶积方法.Zhang等 (2007)针对分层反Q滤波方法，通过在Gabor变换谱上拾取增益控制频率，进而计算增益控制振幅，有效地控制了振幅补偿算子对噪声的放大.
不同于常规的Wiener反褶积，基于Gabor变换的反褶积方法求得一个非平稳的反褶积算子并用于校正震源子波的波形变化和能量衰减，基于Gabor变换反褶积可描述含有地层衰减影响的非平稳性真实野外模型、是常规反褶积的自然延伸，Gabor变换是小时窗的Fourier变换，这些分割的小时窗部分POU (a partition of unity) 的叠加总和是1个单位.对于含有地层衰减影响的非平稳性褶积模型，基于Gabor变换的反褶积可以因子分解为三项：震源子波、衰减函数和反射系数的Gabor变换.
本文利用Margrave等 (2011)提出的基于Gabor变换的反褶积方法，在Gabor域实现对地震波能量衰减的补偿然后在时间域求得地层反射系数，与Wiener反褶积相比较，基于Gabor变换反褶积可补偿中深层的能量衰减并因此拓宽有效频带和提高时间分辨率，求得的地层反射系数更接近真实的地层反射系数，理论模型和实际数据的应用效果都较好.
1 方法原理 1.1 非平稳性反褶积平稳性的地震道褶积模型满足公式为
|
(1) |
在公式 (1) 的褶积模型中，震源子波是不随走时而变化的，这与实际地层的物理影响是不符合的，例如地震波传播过程中遇到的几何扩散、地层黏滞性、能量转换吸收等 (赵建勋和倪克森, 1992; 杨学亭等, 2014)，而基于Gabor变换的反褶积适应实际的地层影响例如地层Q能量吸收、震源子波随走时变化等情况.
首先Gabor变换使用Gaussian窗函数Ωk(t)，公式为
|
(2) |
这里Δτ是Gaussian窗间距，T是一半Gausssian窗宽度.Gabor变换要求∑kΩk(t)=1, ∀t，这样构造的函数叫做POU.
基于Gabor变换反褶积的因子分解中的衰减函数为
|
(3) |
然后结合公式 (1) 的褶积模型，可得到如下公式为
|
(4) |
这就是地震道非平稳性反褶积模型公式 (Margrave et al., 2002)，当Q值→∞，公式 (4) 的积分就变为反射系数r(t) 的Fourier变换r(f) 即常规反褶积模型如同公式 (1).对于有限值Q，公式 (4) 精确的描述地震波旅行过程中的非平稳性的特征.
公式 (4) 可写成矩阵形式为
|
(5) |
W是代表震源子波w(t) 的Toeplitz矩阵，A是代表衰减函数a (τ, t－τ) 的矩阵，R是代表反射系数的向量.WA说明了震源子波在实际地层传播过程中不是一成不变的，而是随走时变化，因此更能代表实际物理过程.
为解得反射系数R，公式 (5) 变为
|
(6) |
公式 (6) 说明基于Gabor变换反褶积就是求取反算子 (WA)－1的过程.公式 (7) 为
|
(7) |
根据公式 (7)，对于某时刻t，地震道的Gabor变换大致是震源子波的Fourier变换、衰减函数和地震反射系数的Gabor变换的乘积.
公式 (7) 有明确的物理意义，将一个很小的时窗应用到地震道上然后做Fourier变换即时频分析 (高静怀和张兵, 2012)，这个时窗越小 (大于子波长度)，它越能描述震源子波在地层旅行过程中的非平稳性特征，它就越接近真实的地层物理响应.
从地震记录里求取地层反射系数的过程就是反褶积，因此在Gabor域满足反射系数为
|
(8) |
μ是一个很小的常数，Amax对应
实际资料的震源子波都是最小相位的，因此创建一个最小相位的地震子波w(图 1) 主频是20 Hz，理论地震道 (图 2a) 是此地震子波与反射系数 (图 3a, 图 4a) 褶积生成，其不包含任何能量衰减例如Q吸收衰减，这个地震道因为深层的高频信息没有被衰减所以能清楚地分辨出在深层 (1 s以下) 的信号，如图 2中的褐色方框所示.
|
图 1 最小相位子波 Figure 1 The minimum phase wavelet |
|
图 2 不同Q值对地震道的影响 Figure 2 The result of a trace with different Q value |
|
图 3 两种反褶积方法求得的反射系数比较 (a) 对图 2中地震道做Wiener反褶积求得的反射系数；(b) 对图 2中地震道做Gabor域反褶积求得的反射系数. Figure 3 The comparison of reflection coefficient with different deconvolution (a) The result of different Q value in Fig.2 with Wiener deconvolution; (b) The result of different Q value in Fig.2 with Gabor deconvolution. |
|
图 4 关联度分析 Figure 4 The correlation analysis (the blue line indicates Wiener deconvolution result, the blue one indicates Gabor deconvolution result) |
为符合实际的地层吸收衰减而测试不同Q值对该地震道的影响，图 2中b、c、d、e、f分别对应的Q值是20、40、60、80、100，对比分析发现对于衰减较弱 (Q>=80) 的地震道2e、2f深层 (1 s以下) 能较能清楚的分辨地层反射信号，因此得出结论是在求地层反射系数即反褶积过程中必须考虑地震波能量的衰减.
对图 2中不同Q值的地震道做Wiener反褶积得到的反射系数 (图 3a中b、c、d、e、f)，其波形特征特别是深层 (1.2 s以下) 与原始反射系数 (图 3a中a) 差别较大，如图 3中的褐色方框所示；同样的对图 2中这些不同Q值的地震道做基于Gabor变换反褶积得到的反射系数 (图 3b中b、c、d、e、f)，其波形特征与原始反射系数 (图 3a或b中a) 基本相似，尤其对于能量衰减不剧烈的地震道 (Q>20)，值得说明的是在深层 (1.2 s以下) Gabor域反褶积求得的反射系数明显好于Wiener反褶积的结果，与原始反射系数基本接近如图 3b中褐色方框所示.
为检验反褶积求得的反射系数与原始反射系数 (图 3a或b中a) 的匹配程度，将这些求得的反射系数与原始反射系数做关联度分析，Wiener反褶积的结果中最好的关联度 (Q=100时) 也只能达到0.47而且还有时延即lag≠0，如图 4中的蓝色曲线所示.
基于Gabor变换反褶积的结果其关联度最好的 (Q=100时) 能达到0.79并且没有时延 (lag=0)，如图 4中的红色曲线所示.这表明基于Gabor变换反褶积考虑到地震子波不是一成不变的并伴有能量衰减的非平稳性特征，因此基于Gabor变换反褶积较好地补偿了中深层能量衰减，并相应地提高了中深层的分辨率 (图 3b中1 s以下).
2.2 实际地震数据应用基于Gabor变换反褶积处理实际地震数据时需要注意的以下两个方面：(1) 在野外资料采集过程中，震源产生的每一炮的震源子波是不同的；(2) 由于测线较长，测线上每一个共中心点 (CMP) 下的地质构造不同并因此导致在每个CMP点的震源子波的衰减系数不同即Q值的估计值不同.针对地震资料上述两个方面的非平稳特征，基于Gabor域反褶积分别采取以下措施：对于第1方面需要做子波反褶积，具体做法是在野外放炮采集时需要记录每一炮的远场子波，并在室内的资料处理过程中每一炮的反褶积用上相对应的远场子波；对于第2方面需要估算在每个CMP点的震源子波的衰减系数即Q值，处理流程如图 5所示.
|
图 5 Gabor域反褶积的处理流程 Figure 5 The workflow of Gabor deconvolution |
对于浅水如渤海湾的地震资料的处理流程中，反褶积是很关键的处理步骤，直接影响最终的偏移成像质量，所以在反褶积这个处理节点上分别试验了不同的反褶积方法：t-x域Wiener反褶积、τ-p域预测反褶积和Gabor域反褶积，并比较其相应的叠加剖面，分别如图 6中a，b，c，d所示.
|
图 6 某航行线反褶积叠加对比 (a) 原始叠加；(b) t-x域Wiener反褶积叠加；(c) 某航行线τ-p域预测反褶积叠加；(d) Gabor域反褶积叠加. Figure 6 The comparison of stack profile with different deconvolution (a) The raw stack; (b) The stack with Wiener deconvolution in T-X domain; (c) The stack with prediction decovolution in τ-p domain; (d) The stack with Gabor deconvolution. |
单从反褶积的结果来看，t-x域Wiener反褶积与τ-p域预测反褶积大致相当，Gabor域反褶积的频谱较好.仅从频谱分析 (图 7) 不能决定选择何种反褶积方法，需要更充分的试验，然后以该试验线为中心上下各扩5 km，用这个小3D数据体分别完成t-x域Wiener反褶积、τ-p域预测反褶积、Gabor域反褶积的整个处理流程即得到最终的叠前时间偏移剖面.
|
图 7 各种反褶积方法频谱分析 Figure 7 The frequency analysis with different deconvolution (the green line indicates raw stack; the wathet one indicates Wiener result; the blue one indicates τ-p result; the red one indicates Gabor result) |
比较图 8中a和b的局部成像 (如绿色箭头所示)，在局部地震同相轴的连续性方面τ-p域预测反褶积不如t-x域Wiener反褶积.下面重点比较t-x域Wiener反褶积和Gabor反褶积的偏移结果，由于偏移后的叠加剖面存在部分随机噪声，经去除后结果比较 (图 9).
|
图 8 反褶积+叠前时间偏移的效果对比 (a) t-x域Wiener反褶积的偏移剖面；(b) τ-p域反褶积的偏移剖面. Figure 8 The comparison of the seismic data with different convolution following the same PSTM (a) The PSTM stack profile with Wiener deconvolution; (b) The PSTM stack profile with prediction deconvolution in τ-p domain. |
|
图 9 反褶积+叠前时间偏移的效果对比 (a) t-x域Wiener反褶积的偏移剖面；(b) Gabor域反褶积的偏移剖面. Figure 9 The comparison of the seismic data with different convolution following the same PSTM (a) The PSTM stack profile with Wiener deconvolution; (b) The PSTM stack profile with Gabor deconvolution. |
经过更充分的试验和对比，分别与t-x域Wiener反褶积、τ-p域反褶积相比较，在提高剖面的时间分辨率和地质波组的丰富程度方面，Gabor域反褶积 (图 9b) 有更好的效果，尤其是图 9中a和b的椭圆形区域内效果较为明显，所以最终选用Gabor域反褶积.
3 结论常用的Wiener反褶积存在某些不合理性，原因在于地层的吸收衰减导致震源子波在其旅行过程中无论波形还是频谱都是变化的、具有非平稳性特征.基于Gabor变换反褶积考虑到震源子波的非平稳性因此更适应实际情况.理论模型的测试表明当地层衰减不太剧烈 (Q>20) 时基于Gabor变换反褶积能很好的求出地层的反射系数包括深层.实际资料处理也发现较之于Wiener反褶积和τ-p域反褶积，基于Gabor变换反褶积的剖面有效频带稍宽、中深层的成像和分辨率都相对较好.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！| [] | Aki K, Richards P G. 2002. Quantitative Seismology[M]. 2nd ed. Sansalito, CA:University Science Books. |
| [] | Bickel S H, Natarajan R R. 1985. Plane-wave Q deconvolution[J]. Geophysics, 50(9): 1426–1439. DOI:10.1190/1.1442011 |
| [] | Chen Z B, Chen X H, Li J Y, et al. 2014. A band-limited and robust inverse Q filtering algorithm[J]. Oil Geophysical Prospecting (in Chinese), 49(1): 68–75. |
| [] | Clarke G K C. 1968. Time-varying deconvolution filters[J]. Geophysics, 33(6): 936–944. DOI:10.1190/1.1439987 |
| [] | Gabor D. 1946. Theory of communication[J]. Journal of the Institution of Electrical Engineers-Part Ⅲ:Radio and Communication Engineering, 93(26): 429–457. DOI:10.1049/ji-3-2.1946.0074 |
| [] | Gao J H, Wang L L, Zhao W. 2009. Enhancing resolution of seismic traces based on the changing wavelet model of the seismogram[J]. Chinese Journal of Geophysics (in Chinese), 52(5): 1289–1300. DOI:10.3969/j.issn.0001-5733.2009.05.018 |
| [] | Gao J H, Zhang B. 2012. Seismic attenuation parameter estimation based on the time-frequency filtering[J]. Oil Geophysical Prospecting (in Chinese), 47(6): 931–936. |
| [] | Guo T C, Cao W J, Tao C J, et al. 2015. Research on time-varying spectral modeling deconvolution method[J]. Geophysical Prospecting for Petroleum (in Chinese), 54(1): 36–42. |
| [] | Kjartansson E. 1979. Constant Q-wave propagation and attenuation[J]. Journal of Geophysical Research, 84(B9): 4737–4748. DOI:10.1029/JB084iB09p04737 |
| [] | Koehler F, Taner M T. 1985. The use of the conjugate-gradient algorithm in the computation of predictive deconvolution operators[J]. Geophysics, 50(12): 2752–2758. DOI:10.1190/1.1441895 |
| [] | Margrave G, Lamoureux M P, Henley D C. 2011. Margrave, Gabor deconvolution:Estimating reflectivity by nonstationary deconvolution of seismic data[J]. Geophysics, 76(3): W15–W30. DOI:10.1190/1.3560167 |
| [] | Margrave G F, Lamoureux M P, Grossman J P, et al. 2002. Gabor deconvolution of seismic data for source waveform and Q correction[C].//72nd Annual International Meeting. SEG Technical Program. Expanded Abstracts, 2190-2193. |
| [] | Peng C, Zhu S J, Huang Z Y, et al. 2007. Dynamic wavelet estimation based on the dynamic convolution model[J]. Geophysical Prospecting for Petroleum (in Chinese), 46(4): 324–328. |
| [] | Robinson E A, Treitel S. 1967. Principles of digital wiener filtering[J]. Geophysical Prospecting, 15(3): 311–332. DOI:10.1111/gpr.1967.15.issue-3 |
| [] | Wang Y H. 2002. A stable and efficient approach of inverse Q filtering[J]. Geophysics, 67(2): 657–663. DOI:10.1190/1.1468627 |
| [] | Wang Y H. 2006. Inverse Q-filter for seismic resolution enhancement[J]. Geophysics, 71(3): V51–V60. DOI:10.1190/1.2192912 |
| [] | Yang X T, Liu C, Liu Y, et al. 2014. The attenuation compensation of seismic wave energy in time-frequency domain based on the continuous wavelet transform[J]. Geophysical Prospecting for Petroleum (in Chinese), 53(5): 523–529. |
| [] | Yilmaz Ö. 2006. Seismic Data Analysis (in Chinese)[M]. Liu H S, Wang K B, Tong S Y, et al., Trans. Beijing:Petroleum Industry Press, 110-160, 180-190. |
| [] | Zhang M M, Dai Y S, Zhang Y N, et al. 2014. Time-variant seismic wavelet estimation method based on spectral modeling in time-frequency domain[J]. Geophysical Prospecting for Petroleum (in Chinese), 53(6): 675–682. |
| [] | Zhang X W, Han L G, Zhang F J, et al. 2007. An inverse Q-filter algorithm based on stable wavefield continuation[J]. Applied Geophysics, 4(4): 263–270. DOI:10.1007/s11770-007-0040-9 |
| [] | Zhao J X, Ni K S. 1992. Cascaded inverse Q filtering and the application[J]. Oil Geophysical Prospecting (in Chinese), 27(6): 722–730. |
| [] | 陈增保, 陈小宏, 李景叶, 等. 2014. 一种带限稳定的反Q滤波算法[J]. 石油地球物理勘探, 49(1): 68–75. |
| [] | 高静怀, 汪玲玲, 赵伟. 2009. 基于反射地震记录变子波模型提高地震记录分辨率[J]. 地球物理学报, 52(5): 1289–1300. DOI:10.3969/j.issn.0001-5733.2009.05.018 |
| [] | 高静怀, 张兵. 2012. 基于时频滤波的吸收衰减参数估算[J]. 石油地球物理勘探, 47(6): 931–936. |
| [] | 郭廷超, 曹文俊, 陶长江, 等. 2015. 时变谱模拟反褶积方法研究[J]. 石油物探, 54(1): 36–42. |
| [] | 彭才, 朱仕军, 黄中玉, 等. 2007. 基于动态褶积模型的动态子波估计[J]. 石油物探, 46(4): 324–328. |
| [] | 杨学亭, 刘财, 刘洋, 等. 2014. 基于连续小波变换的时频域地震波能量衰减补偿[J]. 石油物探, 53(5): 523–529. |
| [] | 渥·伊尔马滋. 2006. 地震资料分析-地震资料处理、反演和解释[M]. 刘怀山, 王克斌, 童思友, 等, 译. 北京: 石油工业出版社, 110-160, 180-190. |
| [] | 张漫漫, 戴永寿, 张亚南, 等. 2014. 基于时频域谱模拟的时变子波估计方法[J]. 石油物探, 53(6): 675–682. |
| [] | 赵建勋, 倪克森. 1992. 串联反Q滤波及其应用[J]. 石油地球物理勘探, 27(6): 722–730. |