0 引言

在三维地震勘探中，受地表障碍物、观测设备、施工条件等因素的限制，采集到的数据往往出现不规则、缺道、坏道等现象。地震数据的这种不规则性会严重影响后续的多次波压制^[1]、偏移成像^[2]、四维地震监测^[3]等环节。为了消除地震数据不规则的影响，一般采用地震数据规则化方法重建缺失数据。但是，由于数据规则化本质上是一个欠定的反演问题，因此一直以来都是地震数据处理领域十分具有挑战性的研究课题。

目前，地震数据规则化方法主要分为三类：第一类是在信号变换域引入稀疏约束恢复地震数据，如Radon变换^[4]、傅里叶变换^[5-8]、小波变换^[9]、曲波变换^[10-11]等；第二类是基于褶积算子的预测滤波器方法，根据已知数据提取局部地震属性(如局部倾角)，再扩展该信息以恢复缺失数据^[12-18]；第三类是模型驱动方法，利用已知的地下模型信息，采用积分延拓算子重建缺失数据，常见的延拓算子包括炮检距延拓^[19]、炮点延拓^[20]、方位角延拓^[21]等。其中前两类方法属于数据驱动的规则化方法，

以往基于模型驱动的数据规则化方法虽然利用了已知的地下模型信息，但是处理精度往往过度依赖于模型的准确性。当地下模型存在一定误差时，规则化结果不理想。本文采用叠前时间偏移和反偏移算子作为积分延拓算子，在成像空间保留了一个维度的冗余信息(如地表炮检距)，使映射到成像空间的过程中不丢失模型的误差信息；再通过引入更加稳定的预条件算子，可以在成像空间更好地提取和应用反射波同相轴特征，从而改善数据规则化质量。最后通过理论模型和实际资料处理验证本文方法能否够有效改善数据规则化效果，改进最终成像质量。

1 方法原理 1.1 最小二乘偏移理论

基于一阶Born近似理论，反射地震数据d可以通过地下反射系数模型m线性映射得到

$ \mathit{\boldsymbol{d}} = \mathit{\boldsymbol{Lm}} $

(1)

式中L为线性映射算子，也称为Born正演算子，可以是叠前/叠后反时间偏移算子，也可以是叠前/叠后反深度偏移算子。该算子与速度模型有关，由于规则化过程中速度模型固定，因此式(1)忽略了速度这个参数。最小二乘偏移最小化目标函数为

$ J\left( \mathit{\boldsymbol{m}} \right) = \frac{1}{2}{\left\| {\mathit{\boldsymbol{Lm}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right\|^2} + \frac{\mu }{2}{\left\| {\mathit{\boldsymbol{Dm}}} \right\|^2} $

(2)

式中：d_obs为观测数据；D为正则化算子；μ为阻尼因子。式(2)的显式解为

$ \mathit{\boldsymbol{m}} = {({\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}} + \mu {\mathit{\boldsymbol{D}}^{\rm{T}}}\mathit{\boldsymbol{D}})^{ - 1}}{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{d}}_{{\rm{obs}}}} $

(3)

定义一个新的向量p，满足

$ \mathit{\boldsymbol{m}} = \mathit{\boldsymbol{Pp}} $

(4)

则式(2)变为

$ J\left( \mathit{\boldsymbol{p}} \right) = \frac{1}{2}{\left\| {\mathit{\boldsymbol{LPp}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right\|^2} + \frac{\mu }{2}{\left\| \mathit{\boldsymbol{p}} \right\|^2} $

(5)

式中P=D^-1为预条件算子。上式的解为

$ \mathit{\boldsymbol{p}} = {({\mathit{\boldsymbol{P}}^{\rm{T}}}{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{LP}} + \mu \mathit{\boldsymbol{I}})^{ - 1}}{\mathit{\boldsymbol{P}}^{\rm{T}}}{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{d}}_{{\rm{obs}}}} $

(6)

式中I为单位矩阵。将式(6)代入式(4)，可得

$ \mathit{\boldsymbol{m}} = \mathit{\boldsymbol{P}}{({\mathit{\boldsymbol{P}}^{\rm{T}}}{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{LP}} + \mu \mathit{\boldsymbol{I}})^{ - 1}}{\mathit{\boldsymbol{P}}^{\rm{T}}}{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{d}}_{{\rm{obs}}}} $

(7)

Clearbout^[22]指出，合理的预条件算子比对应的正则化算子可以更快地加速迭代反演的收敛速度。

由于Born正演(或反偏移)算子L一般无法用矩阵显式表达，而且模型和数据的规模很大，因此一般采用共轭梯度法近似求解式(7)，得到分辨率更高、振幅保真度更好的反射系数模型。在相同的理论框架下，选择不同的反偏移算子，可以构建不同的最小二乘偏移算法，如最小二乘叠前时间偏移^[23-24]、最小二乘逆时偏移^[25-26]等。

1.2 地震数据规则化

当观测数据不规则时，式(5)可以修改为

$ J\left( \mathit{\boldsymbol{p}} \right) = \frac{1}{2}{\left\| {\mathit{\boldsymbol{KLPp}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right\|^2} + \frac{\mu }{2}{\left\| \mathit{\boldsymbol{p}} \right\|^2} $

(8)

式中K为Mask算子，该算子的作用是将缺失地震道置0，非缺失地震道保持不变。通过共轭梯度法最小化目标函数式(8)得到最优解p_opt，再应用预条件算子和反偏移算子，则可以重建缺失地震数据

$ {\mathit{\boldsymbol{d}}_{{\rm{reg}}}} = \mathit{\boldsymbol{LP}}{\mathit{\boldsymbol{p}}_{{\rm{opt}}}} $

(9)

重建后的地震数据可以表示为

$ {\mathit{\boldsymbol{d}}_{{\rm{reg}}}} = \mathit{\boldsymbol{LP}}{({\mathit{\boldsymbol{P}}^{\rm{T}}}{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{KLP}} + \mu \mathit{\boldsymbol{I}})^{ - 1}}{\mathit{\boldsymbol{P}}^{\rm{T}}}{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{K}}^{\rm{T}}}{\mathit{\boldsymbol{d}}_{{\rm{obs}}}} $

(10)

分三种情况分析式(10)中各项对数据规则化的影响。

第一种情况，忽略Hessian矩阵求逆过程(式(10)中的求逆项)，即直接采用偏移+反偏移进行数据重建，式(10)变为

$ {\mathit{\boldsymbol{d}}_{{\rm{reg}}}} = \mathit{\boldsymbol{L}}{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{K}}^{\rm{T}}}{\mathit{\boldsymbol{d}}_{{\rm{obs}}}} $

(11)

这种做法虽然也能够起到数据插值的作用，但是振幅信息不可靠。通过设计一个匹配滤波器可以在一定程度上解决这个问题，但是匹配滤波器的估计本身也是一个欠定的反问题，存在多解性。

第二种情况，忽略预条件项，即直接采用最小二乘偏移+反偏移进行数据重建，式(10)变为

$ {\mathit{\boldsymbol{d}}_{{\rm{reg}}}} = \mathit{\boldsymbol{L}}{({\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{KL}})^{ - 1}}{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{K}}^{\rm{T}}}{\mathit{\boldsymbol{d}}_{{\rm{obs}}}} $

(12)

这种做法可以解决偏移+反偏移中存在的振幅匹配问题，但是无法消除数据缺失带来的成像噪声，该噪声反投影到数据域，会降低插值数据的信噪比。

第三种情况，完整地求解式(10)，即采用预条件最小二乘偏移+反偏移进行数据重建，这种做法既考虑了振幅不匹配的问题，又可以压制数据不规则带来的噪声。

实际上，预条件算子的选取对数据规则化效果的影响很大。本文采用因果积分算子作为预条件算子，该算子的作用是使振幅随炮检距平滑变化。当速度比较准确时，共成像点道集拉平，满足振幅平滑变化的条件；当速度误差较大时，相邻道之间的振幅差异也不会太大，因此也基本能满足条件；当速度误差很大时，需要进一步采用更加合适的预条件算子，以保留原始数据中的剩余旅行时信息。

1.3 叠前时间偏移/反偏移

综合考虑算子的灵活性、计算效率以及预处理阶段对规则化的需求，本文采用Kirchhoff叠前时间反偏移作为线性正演算子L，该算子可以表示为

$ \mathit{\boldsymbol{d}}\left( {y, H, t} \right) = \mathop \sum \limits_{x, h} w\left( {x, y, H, \tau } \right)\mathit{\boldsymbol{m}}\left( {x, h, \tau } \right) $

(13)

式中：m(x，h，τ)为时间偏移共炮检距道集；d(y，H，t)为叠前共炮检距数据；x和h分别为成像域中心点坐标和炮检距；τ为成像域双程旅行时；y和H分别为数据域中心点坐标和炮检距，且H∈[h-Δh, h+Δh], Δh为炮检距间隔；w为加权函数，合适的加权函数可以加快反演的收敛速度。在各向同性介质和直射线假设下，双程旅行时满足双平方根方程

$ t = \sqrt {\frac{{{\tau ^2}}}{4} + \frac{{{{\left( {y - x - h} \right)}^2}}}{{{v^2}}}} + \sqrt {\frac{{{\tau ^2}}}{4} + \frac{{{{\left( {y - x + h} \right)}^2}}}{{{v^2}}}} $

(14)

式中v为均方根速度。

正演算子式(13)的共轭转置L^T为Kirchhoff叠前时间偏移，其表达式为

$ \mathit{\boldsymbol{m}}\left( {x, h, \tau } \right) = \mathop \sum \limits_{y, H} w\left( {x, y, H, \tau } \right)\mathit{\boldsymbol{d}}\left( {y, H, t} \right) $

(15)

1.4 速度的影响

叠前时间偏移共成像点道集属于扩展模型的一种，该模型除了和物理位置有关，还允许模型随炮检距变化。因此，不论背景速度是否正确，使用该模型均可以模拟出和观测数据相同的反射波旅行时信息，都可以在不发生周波跳跃的前提下正确地匹配波形信息^[27]，即不论速度正确与否，随着迭代次数的增加，目标函数式(2)中的第一项(数据拟合项)总是可以逐渐趋近于零。但是，如果速度模型不准确，式(2)中的正则化算子需要适应速度不准确带来的旅行匹配误差。

2 数值试验 2.1 Marmousi模型测试

图 1a为深度域Marmousi层速度模型，图 1b为Dix公式转换后的时间域均方根速度模型。图 1c是公开的波动方程有限差分正演共炮点道集^[28]经过炮检距分选后的(时间，中心点，炮检距)域的叠前数据体。图 1d为图 1c数据随机剔除50%地震道后的叠前数据体。

使用均方根速度场(图 1b)和不同的叠前数据体进行叠前时间偏移和最小二乘叠前时间偏移试验。首先对比炮检距为-0.325km的共炮检距成像剖面。图 2a和图 2b为使用完整的地震数据(图 1c)分别进行常规叠前时间偏移和最小二乘叠前时间偏移得到的共炮检距成像剖面，可以看出：尽管Marmousi模型构造复杂，速度横向变化剧烈，但是时间偏移仍然可以大致准确地对地下构造形态进行成像；经过最小二乘反演之后的成像结果分辨率得到显著提高。图 2c、图 2d和图 2e为使用不完整数据(图 1d)分别进行常规叠前时间偏移、最小二乘叠前时间偏移和预条件最小二乘叠前时间偏移得到的共炮检距成像剖面，可以看出，数据缺失给成像结果带来许多噪声，降低了信噪比，经过最小二乘反演，虽然可以提高分辨率，但是并不能消除数据缺失的影响。通过加入合理的先验信息(比如本文所采用的共成像点道集同相轴沿炮检距平滑变化)，既可以提高成像分辨率，又可以压制成像噪声。

图 3a和图 3b为采用完整数据分别进行常规叠前时间偏移和最小二乘叠前时间偏移的共成像点道集，可以看出：尽管模型复杂，速度横向变化剧烈，但是采用Dix公式转化而来的均方根速度场(图 1b)，共成像点道集基本可以被拉平，该信息可以作为先验信息约束反演成像；另外，经过最小二乘反演，成像分辨率有一定的提升，而且浅层大炮检距的拉伸效应也在一定程度上得到压制。图 3c、图 3d、图 3e为采用不完整地震数据分别进行常规叠前时间偏移、最小二乘叠前时间偏移和预条件最小二乘叠前时间偏移的共成像点道集，可以看出，和共炮检距成像剖面类似，数据缺失在成像道集中也引入噪声，通过最小二乘反演能够提高分辨率，但是无法压制这类噪声，利用共成像点道集同相轴沿炮检距平滑变化这个先验信息，再进行预条件最小二乘叠前时间偏移，不仅能够提高分辨率，而且能够有效压制数据缺失带来的噪声。

对上述成像结果应用叠前时间反偏移即可以进行数据规则化。图 4为5.300km位置处的CMP道集规则化结果，图 5为规则化结果与原始道集误差。由图 4和图 5可以看出，直接采用常规偏移方法也可以对数据进行插值，但是，由于没有考虑Hessian矩阵的影响，插值结果的振幅信息不准确，需要考虑振幅匹配问题，而这个问题也是一个欠定问题，存在多解性。通过最小二乘反演可以解决振幅匹配问题，但是在数据缺失严重的区域，没有先验信息做约束，容易将成像结果中的噪声反投影到数据空间，使插值结果信噪比较低。预条件最小二乘偏移规则化可以同时解决振幅不匹配和信噪比的问题，使得插值结果振幅属性更加可靠、信噪比更高。

需要注意的是，对于本文方法，如果均方根速度比较准确，可以采用沿着炮检距方向的平滑算子作为预条件算子，但是如果均方根速度不准确，则需要采用更加复杂的预条件算子，比如双曲Radon变换、非稳态平面波构建滤波器^[8]等。尽管如此，通过偏移算子的作用，成像空间依然会比数据空间更简单，因此有助于简化预条件算子的设计，这也是模型驱动类数据规则化方法的意义所在。

2.2 实际资料处理

为了进一步验证本文方法的有效性，对海上实际不规则地震数据进行规则化处理。图 6a为海上实际观测的叠前数据体，受拖缆漂移等因素的影响，观测数据中存在明显的地震道缺失；图 6b为叠加速度分析得到的速度场，本文使用该速度场进行时间域成像，并将其作为模型先验信息进行数据规则化处理。

使用图 6b所示的速度场，采用常规叠前时间偏移、最小二乘叠前时间偏移、预条件最小二乘叠前时间偏移进行成像。图 7是炮检距为-0.425km的三种方法共炮检距成像剖面，可以看出，受数据不规则的影响，常规偏移结果中存在明显的成像噪声。经过最小二乘反演，成像分辨率得到显著提高，但是信噪比依然较低。通过在共成像点道集引入平滑约束，进行预条件最小二乘反演之后，不仅可以提高成像分辨率，而且可以有效压制成像噪声。

图 8为常规叠前时间偏移、最小二乘叠前时间偏移、预条件最小二乘叠前时间偏移后得到的位于2.360km处的炮检距域共成像点道集，可以看出，共成像点道集基本拉平，说明图 6b所示的速度场比较准确。另外，也能够明显的看出，受数据缺失的影响，常规叠前时间偏移成像道集的信噪比较低(图 8a)。通过最小二乘反演进行成像，虽然可以提高分辨率，但是无法压制数据缺失带来的噪声(图 8b)。进一步在成像道集引入同相轴沿炮检距方向平滑变化这一先验信息，可以显著提高道集的信噪比，同时也能够保持反演带来的分辨率(图 8c)。

最后，应用叠前时间反偏移对缺失数据进行重建。图 9a为2.360km处的CMP道集，该道集存在明显的数据缺失。图 9b~图 9d分别为采用常规叠前时间偏移、最小二乘叠前时间偏移和预条件最小二乘叠前时间偏移的规则化结果。从结果可以看出，对常规叠前时间偏移结果进行反偏移，虽然也可以起到数据插值的作用，但是由于没有考虑Hessian矩阵的影响，插值出来的振幅信息不可靠。经过最小二乘反演，虽然可以解决振幅信息不匹配的问题，但是受不规则数据引入噪声的影响，插值结果的信噪比较低。通过预条件化技术既可以恢复振幅信息，又能够提高规则化后数据的信噪比。

3 结论与讨论

本文在最小二乘反演的理论框架下，采用最小二乘叠前时间偏移进行数据规则化，利用合理的预条件算子改善数据规则化的处理效果。从理论模型和实际资料处理结果可以得到如下几点结论及认识：

(1) 叠前时间偏移共成像点道集是一种扩展模型，可以缓解模型驱动类数据规则化方法对模型的过度依赖；

(2) 最小二乘叠前时间偏移可以提高成像分辨率，改善振幅保真度，是本文规则化方法的中间结果，可以用作质量监控和后续的成像结果分析；

(3) 成像空间比数据空间的信息更加简单，更有利于设计预条件算子；

(4) 合理的预条件算子可以改善规则化的效果。当然，在复杂地表、复杂地质条件下，也可以采用更加复杂的成像算子，以简化成像空间预条件算子的设计，从而进一步改善规则化效果，但是相应地计算成本会增加，对速度模型的依赖性也会变强。