0 引言

零相位子波地震资料具有分辨率高、与反射界面对应关系明确以及同相轴横向连续性好等优点^[1-3]，因此确保地震资料的零相位性是有效开展构造解释和储层研究的前提。地震资料的处理建立在非零相位子波条件下才能达到预期效果，例如反褶积处理只有在最小相位子波时才能将子波压缩为脉冲，而实际地震子波是混合相位的，因而地震资料经过反褶积处理后存在剩余相位^[4]。此外，实际反射系数序列也不是严格的白噪，加之计算中反褶积算子有限长等因素的影响，都将导致反褶积结果为非零相位。剩余相位的存在会极大降低地震资料的品质，如何消除剩余相位、实现地震资料零相位化一直是地球物理领域研究的难点与热点。

Levy等^[5]、Longbottom等^[6]、White^[7]以及Van der Baan^[8]通过常相位旋转方式和峰度最大值准则估计子波的剩余相位；赫永峰等^[9]、陈必远等^[10]在窄带滤波的基础上，对不同数据进行常相位校正，然后将校正后的窄带数据叠加得到分频相位校正结果；单联瑜等^[11]针对最大方差模计算效率低和解析法计算精度不高的缺点，提出一种新的相位校正方法，即对常用的相位校正判别准则进行改进，实现精确的相位校正；李强等^[12]以零相位合成记录作为标准道，利用最大相似度准则求取地震数据的剩余相位；张亚南等^[13]通过ARMA模型和振幅谱单位化的方式构造纯相位滤波器以描述地震子波相位残余，在最大方差模准则约束下，采用改进的粒子群算法针对残余的子波相位进行准确寻优；徐刚等^[14]在常相位假设的基础上，以离散信息熵作为量化标准对无井记录的地震数据进行相位扫描，进而确定能够将剩余相位影响降至最低的相位值；高少武等^[15]通过构建纯相位因子特征方程计算相位随频率f变化的纯相位因子，实现剩余相位的非常数校正；刘俊州等^[16]研究了基于偏斜度准则的常相位旋转地震子波相位估计方法的优势，并分析了偏斜度准则与峰度准则的差异性及其适用条件；刘传奇等^[17]在常相位假设下，利用不同主频的Ricker子波在相同相位下旋转的峰值时差与相位旋转量具有线性关系这一规律，实现地震资料剩余相位角的估计。

从以上讨论可以看出，前人关于相位校正的研究主要是基于常相位假设。当地震资料的剩余相位在有效频带内为一常量或变化不大时，常相位校正技术可以有效消除剩余相位的影响。但事实上，由于受采集、处理以及地层滤波等各种复杂因素的影响^[18]，地震子波的剩余相位具有频变的特征，例如地层Q吸收会导致相位随频率f发生变化^[19-21]。此时，对不同频率采用同一个常数相位进行校正很难达到有效消除频变剩余相位影响的预期效果。因此，开展频变剩余相位校正并实现地震资料零相位化的研究十分必要。

本文提出零相位合成记录约束下的匹配滤波相位校正技术，该技术将井旁地震道作为输入、零相位合成记录作为期望输出，利用匹配滤波算法得到二者间的相位匹配因子，然后将相位匹配因子与地震数据褶积得到零相位化的地震数据。该方法不仅考虑了实际地震数据剩余相位的频变特性，还充分利用测井合成记录进行了高精度相位校正。理论模型和实际资料处理验证了本文方法能够有效消除频变剩余相位影响，且抗噪性强，具有较高的实用价值。

1 方法原理

常规零相位合成记录约束下相位校正方法流程^{[12, 22]}是：首先用零相位子波与基于测井数据得到的反射系数序列进行褶积，产生零相位合成记录；然后在预先设定的范围内按照一定的步长对井旁地震道进行相位旋转；最后求出经不同相位旋转下的井旁地震道与零相位合成记录的误差平方和，最小误差平方和对应的相位即为地震数据的剩余相位。该方法没有考虑子波剩余相位的频变特征，对于剩余相位具有频变特征的地震资料而言，很难达到期望的校正效果。

在地震资料处理中，匹配滤波通常用于陆上不同震源数据的匹配，或浅海与陆地过渡区域地震数据匹配^[23-25]。设b(t)是长度为n的输入数据；d(t)是长度为m+n的期望输出数据；a(t)是长度为m的匹配滤波器因子。输入数据b(t)经过a(t)滤波后的实际输出为

$ y(t) = a(t)*b(t) = \sum\limits_{\tau = 0}^m a (\tau )b(t - \tau ) $

(1)

匹配滤波的实际输出与期望输出的总能量误差为

$ \begin{aligned} E & =\sum_{t=0}^{m+n}[y(t)-d(t)]^2 \\ & =\sum_{t=0}^{n+n}\left[\sum_{\tau=0}^m a(\tau) b(t-\tau)-d(t)\right]^2 \end{aligned} $

(2)

为使总能量误差最小，匹配滤波因子a(τ)需满足

$ \frac{{\partial E}}{{\partial a(s)}} = 2\sum\limits_{t = 0}^{m + n} {\left[ {\sum\limits_{\tau = 0}^m a (\tau )b(t - \tau ) - d(t)} \right]} b(t - s) = 0 $

(3)

整理式(3)，可得

$ \begin{aligned} \sum_{\tau=0}^m a(\tau) & \sum_{t=0}^{m+n} b(t-\tau) b(t-s) \\ & =\sum_{t=0}^{m+n} d(t) b(t-s) \quad s=0,1, \cdots, m \end{aligned} $

(4)

令

$ {r_{bb}}(\tau - s) = \sum\limits_{t = 0}^{m + n} b (t - \tau )b(t - s) $

(5)

$ {r_{db}}(s) = \sum\limits_{t = 0}^{m + n} d (t)b(t - s) $

(6)

式中：r_bb(τ－s)是时间延迟为τ－s的b(t)的自相关；r_db(s)是时间延迟为s的b(t)与d(t)的互相关。于是式(4)可以写成

$ \sum\limits_{\tau = 0}^m {{r_{bb}}} (\tau - s)a(\tau ) = {r_{db}}(s)\quad s = 0, 1, \cdots , m $

(7)

式(7)的矩阵表达形式为

$ \left( {\begin{array}{*{20}{c}} {(1 + \lambda ){r_{bb}}(0)}&{{r_{bb}}(1)}& \cdots &{{r_{bb}}(m)}\\ {\begin{array}{*{20}{l}} \begin{array}{l} {r_{bb}}(1)\\ \end{array} \end{array}}&{\begin{array}{*{20}{l}} \begin{array}{l} (1 + \lambda ){r_{bb}}(0)\\ \end{array} \end{array}}&{\begin{array}{*{20}{l}} \cdots \\ \vdots \end{array}}&{{r_{bb}}(m - 1)}\\ {{r_{bb}}(m)}&{{r_{bb}}(m - 1)}& \cdots &{(1 + \lambda ){r_{bb}}(0)} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {a(0)}\\ {a(1)}\\ \vdots \\ {a(m)} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{r_{db}}(0)}\\ {{r_{db}}(1)}\\ \vdots \\ {{r_{db}}(m)} \end{array}} \right) $

(8)

式中：λ为白噪系数，其作用是防止匹配滤波因子产生畸变；左端b(t)的自相关矩阵为一对称矩阵，称为Toeplitz矩阵。求解这种矩阵方程的方法很多，其中Levision递推算法可对矩阵快速求解，从而得到匹配滤波因子a(t)。

如果将井旁地震道作为输入，零相位合成记录作为输出，则基于式(8)得到的匹配滤波因子a(t)与地震数据褶积后，其结果与零相位合成记录总能量误差最小，即最佳逼近零相位。

本文实现地震资料剩余相位校正的处理步骤如下：

(1) 地震资料频谱分析，求出与地震资料振幅谱相当的零相位子波；

(2) 利用测井数据中的声波、密度资料得到反射系数序列，并与上一步求出的零相位子波褶积制作合成记录；

(3) 将井旁地震道作为输入，步骤(2)中的零相位合成记录作为输出，利用式(8)求出匹配滤波因子a(t)；

(4) 将匹配滤波因子a(t)与地震数据褶积，实现地震数据剩余相位校正。

从上述流程可以看出，该方法不受常相位假设的限制，能够在零相位合成记录的约束下，根据地震资料本身的相位特征求出对应的相位匹配因子a(t)。将a(t)与地震数据褶积后的结果可以消除剩余相位的影响，从而改善剖面质量，提高地震资料的分辨率。

2 理论模型试算

利用地层Q吸收频散效应得到相位随频率发生变化的子波^[26]，具体做法是：①基于常Q模型，模拟得到传播不同时间后的地震子波；②通过反Q滤波进行振幅补偿，得到只存在相位变化、无振幅能量衰减的地震子波。主频为35 Hz的Ricker子波在Q=120的地层中传播时间分别为0.6、0.9、1.2、1.5 s时，对应的频变相位子波及其振幅谱、相位谱如图 1所示(为方便对比和相位谱的求取，图中的地震子波按传播时间进行了时移)，通过反Q滤波振幅补偿使频变相位子波的振幅谱与原始Ricker子波的振幅谱保持一致，频变相位子波波形及振幅的变化只与相位有关。为方便后面叙述，依次将0.6、0.9、1.2、1.5 s传播时间的子波记为频变相位子波1、2、3和4(图 1左)。可以发现，与原始零相位子波相比，频变相位子波的波形发生显著变化，主峰值偏离零时刻，左旁瓣能量增强，右旁瓣能量减弱；相位频变现象越明显，波形变化越显著，当相位频变到一定程度时，子波主峰值的能量与左旁瓣的能量基本相当(例如频变相位子波3和4)，此时无法通过子波的主峰值识别地层反射界面。

将零相位Ricker子波作为标准道，对图 1左中的频变相位子波进行常相位校正，求得频变相位子波1~4对应的剩余相位分别为－63°、－85°、－102°和－117°，基于求得的剩余相位对频变相位子波进行常相位校正，得到图 2所示结果。由图可见，经过校正后的频变相位子波在一定程度上降低了剩余相位的影响，例如频变相位子波1的相位变化范围由校正前的0°~－115°变为校正后的63°~－52°，主峰值也基本对应到零时刻，但与零相位子波的波形仍然存在一定差异；另外，校正后频变相位子波右旁瓣的能量强于左旁瓣，且接近主峰值的能量，导致仍然无法通过主峰值识别地层反射界面。因此，当地震子波存在明显频变剩余相位时，常相位校正技术不能完全消除剩余相位的影响，得不到真正的零相位子波。

将图 1左中的频变相位子波作为输入、主频为35 Hz的零相位Ricker子波作为输出，基于本文的匹配滤波算法得到如图 3所示的相位匹配因子。可见，相位匹配因子的相位谱也随频率发生变化，在同一频率处，相位匹配因子与频变相位子波的相位谱数值大小相等、符号相反。由于时间域的褶积等同于频率域的乘法，因此将地震数据与相位匹配因子褶积后，相当于大小相等、符号相反的子波相位谱与匹配因子的相位谱相加，故二者褶积后可以有效消除不同频率处的剩余相位，实现分频相位校正。基于匹配滤波算法得到的相位校正结果如图 4所示，从波形(图 4左)和相位谱(图 4右)可知，以上4种频变相位子波经过本文方法处理后均变为零相位子波，证实了本文方法的有效性。

以上是基于子波开展的研究，为进一步验证本文方法的有效性，采用频变合成记录开展了深入研究。分别将主频为35 Hz的零相位Ricker子波、频变相位子波3与随机反射系数(图 5左)进行褶积，得到图 5右所示的合成记录。与零相位合成记录相比，频变剩余相位造成地震记录上存在明显“下拉”现象，且局部存在波形变化(图 5右箭头指示处)。

以零相位合成记录作为标准道，对频变相位合成记录进行常相位校正，得到图 6所示结果。由图可见，频变相位合成记录经过常相位校正后，基本解决了地震数据“下拉”现象，但局部位置的波形特征、振幅值仍与零相位合成记录存在差异(图 6箭头指示处)。

将图 5右频变相位合成记录作为输入、零相位合成记录作为输出，利用本文方法得到图 7左所示相位匹配因子。由于图 5左反射系数是随机生成的，因此图 5右频变相位记录各道之间互不相关，然而图 7左得到的相位匹配因子却是一致的，这说明相位匹配因子取决于地震记录中子波的相位信息，不受反射系数干扰。将图 5b的频变相位合成记录与图 7a对应的相位匹配因子进行褶积，得到相位校正结果(图 7右)，由图可见，本文方法处理后的频变相位记录与零相位记录完全一致，说明基于零相位合成记录作为期望输出得到的相位匹配因子也能够实现地震数据的分频相位校正。

实际地震数据不可避免会受到随机噪声的干扰，因此处理方法的抗噪性直接关系到其实用性。向图 5b的频变相位记录加入不同水平的随机噪声以测试本文方法的抗噪性。仍然将零相位合成记录作为输出，以含噪频变相位记录作为输入，利用匹配滤波算法得到图 8所示的相位匹配因子。可见，虽然不同信噪比(SNR)频变相位记录对应的相位匹配因子局部存在略微差异，但整体上基本一致。将图 8中的相位匹配因子与其对应的含噪频变相位记录褶积后，得到图 9所示的相位校正结果。与零相位记录对比可知，即使在低信噪比(SNR=5.7 dB)时，仍然能够取得较好效果，表明本文方法抗噪性强。

3 实际资料处理

A1井为渤海M油田的评价井，在东二下段Ⅰ、Ⅲ油组均钻遇油层，为落实这两套油层的地震响应特征，根据目的层段地震资料频谱分析结果，采用主频为20 Hz的Ricker子波制作零相位合成记录进行井震标定。

基于原始地震数据得到的井震标定结果如图 10a所示，零相位合成记录与井旁地震道的相关系数为0.73，对应的时深关系与A1井原始VSP时深关系相比向下“漂移”18 ms。合成记录与原始地震数据的主要波组特征基本一致，但局部细节存在差异，例如东二下段Ⅲ油组顶在零相位合成记录上为强波峰响应，而在原始地震资料上却对应弱波峰响应，二者振幅能量差异非常明显。

以零相位合成记录作为标准道，对原始地震数据进行常相位校正(求得的剩余相角为－55°)。常相位校正后地震数据的井震标定结果如图 10b所示。此时零相位合成记录与井旁地震道的相关系数为0.78，对应的时深关系与A1井原始VSP时深关系相比向下“漂移”14 ms。从标定结果可以看出，原始地震数据经过常相位校正后，东二下段Ⅲ油组顶对应的波峰振幅能量得到加强，与零相位合成记录的吻合程度有一定的提高，但在细节上仍然存在差异。

采用本文方法对原始地震数据进行相位校正，得到图 10c所示的井震标定结果，此时合成记录与井旁地震道的相关系数提高到0.92，对应的时深关系与A1井原始VSP时深关系相比向下“漂移”10 ms。从图 10c可以看出，基于本文方法处理后，目的层东二下段Ⅰ、Ⅲ油组在地震资料与合成记录上的响应特征基本一致，有效解决了井震之间局部细节差异大的问题，同时分辨率也得到一定提高(红色箭头指示处)，为目的层的构造、储层研究提供了可靠的地震资料。

利用反射系数序列与维纳滤波算法相结合的确定性子波提取技术^[27]依次提取原始资料、常相位校正后资料以及本文方法处理后的地震资料目的层段的子波。图 11为基于原始地震资料提取的子波及其振幅谱、相位谱，从子波波形及相位谱可以看出，原始地震资料存在明显的剩余相位，且剩余相位具有明显的频变特征，最大剩余相位可达－90°。原始地震资料经过常相位校正后仍然存在剩余相位(图 12)，在有效频带范围内，剩余相位变化范围为－4°~43°，但剩余相位频变现象明显减弱。基于本文方法处理后地震资料提取的地震子波关于t=0对称(图 13a)，与零相位子波特征基本一致，且不同频率成分的相位基本为零(图 13c)。

利用未参与相位校正处理的A2井进一步验证本文方法的可靠性。图 14a~图 14c分别为插入A1井、A2井零相位合成记录的原始剖面、常相位校正以及本文方法校正处理后的地震剖面，两口井合成记录的子波均是主频为20 Hz的零相位Ricker子波。对比可见，基于本文方法处理后的地震数据(图 14c)与A2井合成记录的匹配关系最好，而在图 14a和图 14b中都存在与A2井合成记录不相匹配的现象(图中黑色箭头指示处)。在振幅谱一致的情况下(图 15)，图 14c地震剖面的分辨率也是最高的。例如，东二下段Ⅰ油组顶真实的地震响应该为波谷(图中红色虚线)，而在图 14a、图 14b中局部显示为复波，经过本文方法处理后，图 14c的分辨率得到显著提高，有效地恢复了东二下段Ⅰ油组顶对应的波谷响应(图中红色箭头指示处)。

4 结束语

针对常相位校正技术存在无法有效消除频变剩余相位的影响这一问题，本文提出零相位合成记录约束下利用匹配滤波算法实现地震资料的分频相位校正，即将井旁地震道作为输入、零相位合成记录作为输出，利用匹配滤波算法求出二者间的相位匹配因子。在同一频率处，相位匹配因子与地震子波的相位谱数值大小相等、符号相反，因此实际地震数据与相位匹配因子褶积后，可以消除各个频率成分的剩余相位，即实现地震资料剩余相位的分频校正。理论模型和实际资料的处理结果验证了本文方法的有效性，且该方法抗噪声干扰能力强，具有较高的实用价值。

在实际应用中，一定要选择可靠性较高的测井数据制作零相位合成记录，以保证相位匹配因子准确可靠。此外，实际地震资料的剩余相位除具有频变特点外，还具有时变、空变的特点，因此本文方法只有针对目的层位地震数据进行处理时才能得到最佳效果。由于常相位是频变相位的一种特例，故本文方法也适应于常相位校正。