随着地震勘探程度的不断深入, 人们对油气勘探精度的要求越来越高, 地下深部复杂构造成像成为了业界研究的重点, 逆时偏移由于不受角度限制和复杂构造成像较好的特点^[1-3]而得到广泛应用。相对传统逆时偏移, 最小二乘逆时偏移具有更好的振幅均衡性、更高的分辨率和保真度, 对地下复杂构造有更高的成像精度, 但计算量大的问题限制了它的进一步发展。

早在20世纪80年代, TARANTOLA^[4]就提出了最小二乘偏移理论, 但由于计算量远远大于常规偏移, 因此发展缓慢。随着计算机硬件的发展, 最小二乘偏移理论得到进一步丰富。NEMETH等^[5]、KUEHL等^[6]和黄建平等^[7]先后将最小二乘偏移理论应用于单程波方程, 证明了最小二乘偏移可以大幅度提高偏移成像的精度。DAI等^[8-11]提出了最小二乘逆时偏移理论, 并将多震源编码方式引入到最小二乘逆时偏移, 提高了计算效率。屠宁^[12]将压缩感知理论应用于最小二乘逆时偏移, 大大降低了最小二乘逆时偏移的计算成本。目前, 多震源偏移技术被广泛用于提高最小二乘逆时偏移的计算效率, 但由于超道集之间的相关性较高, 导致其成像结果存在较强的串扰噪声, 收敛效果较差。李振春等^[13]、刘玉金等^[14]和FAN等^[15]将合理的约束条件引入最小二乘逆时偏移, 提高了收敛速度, 压制了串扰噪声。DAI等^[9]将去模糊化滤波算子引入多震源最小二乘逆时偏移中, 在加快收敛速度的同时, 压制了串扰噪声。FOMEL^[16]和DUTTA^[17]将Seislet变换算子作为预条件算子, 压制了多震源编码数据直接成像产生的串扰噪声。还有一些学者通过采用平面波编码等方式降低超道集的相关性来压制串扰噪声。李庆洋等^[18]和李闯等^[19]分别研究了基于多震源编码和平面波编码的最小二乘逆时偏移, 指出平面波编码策略进一步提高了最小二乘逆时偏移的计算效率及收敛速度。随后, 李振春等^[20]将平面波最小二乘逆时偏移(PLSRTM)应用于各向异性介质成像, 使其更符合实际地下情况。但由于平面波最小二乘逆时偏移存在一定的偏移假象, 导致收敛效果较差, 影响最终的成像精度。LI等^[21-22]首先将奇异谱分析算子作为正则化算子引入平面波最小二乘逆时偏移, 压制偏移假象, 提高了成像精度; 然后将随机优化与奇异谱分析算子相结合, 进一步提高平面波最小二乘逆时偏移的计算效率和收敛速度。朱峰等^[23]将构造信息引入平面波最小二乘偏移, 提高了成像精度。

本文在实现平面波最小二乘逆时偏移的基础上, 利用基于局部倾角的构造预测方法获取构造信息, 计算构造导向滤波算子^[24-25], 并将其作为正则化算子引入PLSRTM, 压制偏移假象。通过简单模型和HESS模型试算验证了方法的适用性和有效性。

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

对于各向同性介质, LSRTM理论已经很成熟, 波场传播过程通常采用常密度声波方程表达^[8]:

$ \frac{1}{{{\mathit{\boldsymbol{v}}^2}}}\frac{{{\partial ^2}\mathit{\boldsymbol{u}}}}{{\partial {t^2}}} - {\nabla ^2}\mathit{\boldsymbol{u}} = \mathit{\boldsymbol{f}} $

(1)

式中:v为速度场, f为震源项, u为波场。

根据Born近似理论, 分别对速度场和波场进行Born近似可得:

$ \mathit{\boldsymbol{v}} = {\mathit{\boldsymbol{v}}_0} + \delta \mathit{\boldsymbol{v}} $

(2)

$ \mathit{\boldsymbol{u}} = {\mathit{\boldsymbol{u}}_0} + \delta \mathit{\boldsymbol{u}} $

(3)

式中:v₀为背景速度场, δv为扰动速度场, u₀为背景波场, δu为扰动波场。

同时, 假设反射系数为:

$ \mathit{\boldsymbol{m}} = \frac{{2\delta \mathit{\boldsymbol{v}}}}{{{\mathit{\boldsymbol{v}}_0}}} $

(4)

将公式(2)~公式(4)代入公式(1), 可得到基于Born正演的声波方程, 即反偏移方程:

$ \frac{1}{{\mathit{\boldsymbol{v}}_0^2}}\frac{{{\partial ^2}\delta \mathit{\boldsymbol{u}}}}{{\partial {t^2}}} - {\nabla ^2}\delta \mathit{\boldsymbol{u}} = \mathit{\boldsymbol{m}}\frac{1}{{\mathit{\boldsymbol{v}}_0^2}}\frac{{{\partial ^2}{\mathit{\boldsymbol{u}}_0}}}{{\partial {t^2}}} $

(5)

LSRTM可以在数据域和成像域实现。在数据域实现LSRTM的关键步骤为求取目标函数和梯度, 目标函数可以表示为:

$ \mathit{\boldsymbol{J}} = \frac{1}{2}\left\| {\mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}}} \right\|_2^2 $

(6)

梯度可以表示为:

$ \frac{{\partial \mathit{\boldsymbol{J}}}}{{\partial \mathit{\boldsymbol{m}}}} = {\mathit{\boldsymbol{L}}^{\rm{T}}}\left( {\mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}}} \right) $

(7)

其中, L为Born正演算子, L^T为Born正演算子的伴随算子, d为观测到的炮记录数据, m为偏移结果。本文将Born正演算子的伴随算子近似为正演算子的逆, 即偏移算子。

对目标函数的求解采用共轭梯度法, 由于需要多次迭代, 计算量巨大。本文通过引入正则化约束条件来加快收敛速度:

$ \mathit{\boldsymbol{J}} = \frac{1}{2}\left\| {\mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}}} \right\|_2^2 + \frac{1}{2}\lambda \left\| {\mathit{\boldsymbol{zm}}} \right\|_2^2 $

(8)

式中:λ是正则化参数, 作为先验信息约束成像结果; z为用来约束成像结果的正则化算子。

1.2 基于平面波编码的最小二乘逆时偏移 1.2.1 平面波逆时偏移

为了实现平面波逆时偏移(PRTM), 首先将炮记录进行平面波编码, 其相当于一个倾斜叠加的过程。

选用不同的射线参数p, 对原始炮记录进行倾斜叠加^[8]:

$ p = \frac{{{\rm{d}}t}}{{{\rm{d}}x}} = \frac{{\sin \theta }}{{v\left( \theta \right)}} $

(9)

式中:θ为平面波入射的角度, v(θ)为地表沿θ方向上的速度。

经过平面波编码的炮记录d(x, p)可以表示为:

$ \mathit{\boldsymbol{d}}\left( {x,p} \right) = \sum\limits_{{\mathit{\boldsymbol{x}}_{\rm{s}}}} {\mathit{\boldsymbol{d}} * \mathit{\boldsymbol{\delta }}\left( {t - p \cdot \left| {{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right|} \right)} $

(10)

式中, δ(t－p·|x_s|)为时移函数, *为卷积运算, x_s为震源坐标, p·|x_s|为时移量。

平面波逆时偏移是将经过平面波编码的震源波场和以炮记录作为边界反传获得的检波点波场互相关, 并将各个平面波射线参数得到的偏移结果进行叠加:

$ \mathit{\boldsymbol{m}}\left( {x,p} \right) = \sum\limits_{i = 1}^N {\mathit{\boldsymbol{m}}\left( {{p_i}} \right)} $

(11)

式中, m(x, p)为最终PRTM成像结果, N为平面波射线参数的个数, m(p_i)为第i个平面波射线参数获得的成像结果, p_i为第i个平面波射线参数。

1.2.2 平面波最小二乘逆时偏移

将平面波编码技术引入LSRTM即PLSRTM。通过对观测数据和计算数据进行平面波编码, 推导出PLSRTM的目标函数, 即将公式(10)和公式(11)代入公式(8), 得到正则化约束的目标函数:

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{J}} = \frac{1}{2}\left\| {\mathit{\boldsymbol{Lm}}\left( {x,p} \right) - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}\left( {x,p} \right)} \right\|_2^2 + }\\ {\;\;\;\;\;\;\frac{1}{2}{\lambda ^2}\left\| {\mathit{\boldsymbol{zm}}\left( {x,p} \right)\left| {_2^2} \right.} \right.} \end{array} $

(12)

梯度表示为:

$ \frac{{\partial \mathit{\boldsymbol{J}}}}{{\partial \mathit{\boldsymbol{m}}}} = {\mathit{\boldsymbol{L}}^{\rm{T}}}\left[ {\mathit{\boldsymbol{Lm}}\left( {x,p} \right) - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}\left( {x,p} \right)} \right] - {\lambda ^2}{\mathit{\boldsymbol{z}}^{\rm{T}}}\mathit{\boldsymbol{zm}}\left( {x,p} \right) $

(13)

式中, d_obs(x, p)为经过平面波编码的观测数据。

1.3 平面波最小二乘逆时偏移的优化

本文采用基于局部倾角的构造预测方法来获取构造信息^[24-25], 其中局部倾角信息通过平面波分解滤波器(PWD)建立数据的局部平面波模型来求取, 如图 1所示。首先将预测步长作为半径, 利用预测因子由相邻地震道预测得到目标地震道; 然后利用构造信息计算出构造导向滤波算子, 并将其作为正则化算子引入PLSRTM, 即基于构造导向滤波算子的平面波最小二乘逆时偏移(SPLSRTM)。

地震剖面由一系列地震道组成, 即s=[s₁ s₂ … s_N]^T, 其中N为地震道的总数, s_i为第i道数据。用线性算子^[25]表示平面波分解算子:

$ \mathit{\boldsymbol{B}} = \mathit{\boldsymbol{D}}\left( \sigma \right)\mathit{\boldsymbol{s}} $

(14)

式中, B为原始道与预测道的残差, D是非平稳PWD算子。

$ \mathit{\boldsymbol{D}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{I}}&0&0& \cdots &0\\ { - {Q_{1,2}}\left( {{\sigma _1}} \right)}&\mathit{\boldsymbol{I}}&0& \cdots &0\\ 0&{ - {Q_{2,3}}\left( {{\sigma _2}} \right)}&\mathit{\boldsymbol{I}}& \cdots &0\\ \vdots&\vdots&\vdots &{}& \vdots \\ 0&0& \cdots &{ - {Q_{N - 1,N}}\left( {{\sigma _{N - 1}}} \right)}&\mathit{\boldsymbol{I}} \end{array}} \right] $

(15)

式中:I为单位矩阵, σ_i为局部倾角, Q_{i, j}(σ_i)为由第i道到第j道的预测因子。本文采用共轭梯度法来计算预测误差B, 通过预测误差B求取地震道的局部倾角σ, 计算出每一道的预测因子Q_{i, j}(σ_i), 进而得到构造导向滤波算子T=D^-1。

将构造导向滤波算子T作为正则化算子代入新变量a(x, p):

$ \mathit{\boldsymbol{a}}\left( {x,p} \right) = \mathit{\boldsymbol{Tm}}\left( {x,p} \right) $

(16)

将公式(16)代入公式(12), 得到SPLSRTM目标函数:

$ \begin{array}{l} \mathit{\boldsymbol{J}} = \frac{1}{2}\left\| {\mathit{\boldsymbol{LDa}}\left( {x,p} \right) - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}\left( {x,p} \right)} \right\|_2^2 + \\ \;\;\;\;\;\;\frac{1}{2}{\lambda ^2}\left\| {\mathit{\boldsymbol{a}}\left( {x,p} \right)} \right\|_2^2 \end{array} $

(17)

本文采用共轭梯度法求解目标函数(17), 得到最优解a(x, p), 并代入公式(16), 得到最终偏移成像结果m(x, p)。整个过程相当于对每次迭代的梯度利用构造导向滤波器进行滤波, 以压制偏移假象, 从而提高观测数据与计算数据的匹配程度, 达到加快收敛速度的目的。

2 应用实例 2.1 简单模型

选用的简单模型如图 2所示, 横向和纵向网格数分别为1024和231, 横向和纵向间隔分别为10m和8m。选取主频为20Hz的雷克子波作为震源进行Born正演模拟, 炮点数为1024, 每炮1024道, 记录长度为2.5s, 采样间隔1ms。将正演模拟得到的炮记录进行平面波编码, 得到平面波射线参数p在-0.1~0.1ms/m的40个平面波记录。同时, 利用多震源采集方法获得10炮同时激发、每炮1024道接收的10个超道集Born正演模拟记录。

在SPLSRTM中, 每次迭代需要沿着成像倾角的方向进行平滑, 并且成像倾角信息的准确性会影响最终的滤波效果及成像精度。本文对反射系数和第一次迭代梯度进行了倾角拾取对比, 如图 3所示。可以看出, 从第一次迭代梯度中获取的倾角信息具有较高的准确性。在偏移过程中, 可以从每一次迭代梯度中获取局部倾角信息。

分别对10个多震源编码的超道集记录进行了多震源最小二乘逆时偏移(MLSRTM), 对5个和10个平面波记录进行了PLSRTM偏移, 如图 4所示。对5个和10个平面波记录进行了SPLSRTM偏移, 如图 5所示。

结合图 4和图 5可以看出, 在相同计算量的前提下, PLSRTM的偏移结果明显好于MLSRTM。随着平面波射线参数p的增加或迭代次数的增加, 偏移假象在一定程度上得到压制, 但大倾角等复杂构造的PLSRTM偏移结果成像精度较低。在迭代次数和平面波射线参数相同的前提下, SPLSRTM相对PLSRTM对大倾角等复杂构造(红色方框部分)具有更高的成像精度, 同相轴能量更均衡, 噪声得到了更好的压制。

图 6对比了不同情形的收敛曲线, 可以看出:①PLSRTM数据残差相对MLSRTM更接近0, 收敛效果更好; ②随着平面波射线参数的增加, PLSRTM的数据残差减小的速度变快; ③在相同平面波射线参数情况下, SPLSRTM相对PLSRTM的数据残差更接近0, 减小的速度更快。

2.2 HESS模型

为了更好地验证本文方法的适用性, 选取构造更为复杂的HESS模型(图 7)进行了模拟测试。模型横向和纵向网格数分别为1024和453, 横向和纵向采样间隔均为30m。选取主频为20Hz的雷克子波作为震源进行了Born正演模拟, 炮点数为1024, 每炮1024道, 记录长度为4.1s, 采样间隔1ms。将正演模拟得到的炮记录进行平面波编码, 得到平面波射线参数p在-0.1~0.1ms/m的40个平面波记录。

对5个平面波记录分别进行PLSRTM和SPLSRTM偏移, 成像结果和收敛曲线分别如图 8和图 9所示。

对比图 8a和图 8b可见, 随着迭代次数的增加, PLSRTM可以在一定程度上压制串扰偏移假象, 成像结果接近反射系数模型。对比图 8a和图 8c可知, 在相同迭代次数的情况下, SPLSRTM的偏移假象更少, 收敛效果更好, 同相轴能量更为均衡。由图 9可知, 当模型比较复杂的时候, SPLSRTM的收敛效果同样好于PLSRTM, 收敛速度也更快。

为了更好地分析PLSRTM和SPLSRTM两种偏移方法对于复杂模型成像精度的差别, 在图 8b和图 8d中选取红色方框区域进行放大显示, 如图 10所示。对比图 10a和图 10b可以看出, 选取较少的平面波记录进行偏移成像时, PLSRTM偏移结果存在严重的偏移假象; 而SPLSRTM偏移结果同相轴能量相对均衡, 偏移假象得到压制, 成像精度较高。

3 结论

本文将平面波编码的炮记录直接应用于最小二乘逆时偏移以提高其计算效率, 同时针对复杂构造的成像精度及收敛等问题, 将构造导向滤波算子引入平面波最小二乘逆时偏移, 实现了基于构造导向滤波算子的平面波最小二乘逆时偏移。通过简单和复杂构造模型测试, 得到以下几点认识:

1) 在相同计算量的前提下, PLSRTM相对MLSRTM的串扰噪声影响较小, 收敛效果改善明显。

2) 随着平面波射线参数的增加, PLSRTM的收敛效果得到一定改善, 但是计算量增加, 且复杂构造的成像精度较低。

3) 对于复杂构造模型, SPLSRTM可以在保证计算效率的前提下压制偏移假象, 提高成像精度, 改善收敛效果。

需要指出的是, 本文模型试算结果基于理论模型数据, 速度模型已知, 并且在波场传播过程中进行了各向同性假设。在实际应用中, 由于很难提供准确的速度参数场和各向异性参数场, 而最小二乘逆时偏移对速度和各向异性参数很敏感, 因此如何降低最小二乘逆时偏移对于速度和各向异性参数的敏感性是未来研究的重要方向。