预条件弹性介质最小二乘逆时偏移
刘梦丽,
徐兴荣,
王小卫,
胡书华, ...
岩性油气藏 ![]() ![]() |
![]() ![]() ![]() |
弹性波动理论的发展以及长期的勘探实践均已表明,地震波是以矢量波的形式在地球介质中传播的。建立在声波介质假设上的问题只能主要围绕纵波速度和密度进行,忽略了横波,不符合精细勘探的要求[1-3]。为此针对精细岩性油气藏勘探,发展弹性介质假设下的成像方法迫在眉睫。逆时偏移是当前偏移方法中较为精确的深度偏移方法且发展迅速,已经从叠后发展到叠前,从二维到三维,从声波到弹性波,从各向同性到各向异性[4-7]。弹性波逆时偏移主要分为2个过程,首先构建震源波场和检波波场,利用成像条件得到成像结果。1986年Chang等[8]首次完成了各向同性介质弹性波叠前逆时偏移;随后Chang等[9-10]又深入研究,完成了2D和3D弹性波逆时偏移;董良国等[11]分析不同地形起伏情况下自由边界的情况,模拟出地表起伏情况下弹性波复杂的传播现象,实现了弹性波起伏地表自由边界条件的数值计算;为了深入研究弹性波波场传播规律,陈可洋[12]将波印廷矢量应用于多波多分量各向异性介质弹性波波动方程中,实现了上行波、下行波、左行波和右行波的波场分离。针对弹性波成像的串扰问题,王维红等[13]依据各向同性介质一阶速度-应力方程组,利用Helmholtz分解提取纯纵波和纯横波波场,使用震源归一化的互相关成像条件获得纯波成像,避免了纵横波串扰问题。张伟等[14]利用高阶交错网格有限差分数值方法构建了矢量的纵横波波场,将震源归一化的内积成像条件应用于分离后的纯纵波和横波矢量场,由此得到的转换波成像避免了传统弹性波成像方法中出现的极性反转。成像条件是影响弹性逆时偏移的另一个重要因素。杜启振等[15]基于激发时间成像条件实现了横向各向同性介质中的多波多分量叠前逆时偏移;张晓语等[16]基于能量守恒定理及能量密度,实现了弹性波能量互相关成像,能够压制背向散射;张智等[17]提出稳定的激发振幅成像条件,在震源波场和检波器波场的传播过程中,保存最大能量密度的时刻和相应的波场值,归一化后获得角度依赖的反射系数成像剖面,与归一化成像条件相比,稳定激发振幅成像条件具有更小的计算量。
常规的逆时偏移使用的是反偏移算子的共轭,而不是它的逆,因此造成成像结果模糊化,成像精度难以满足精细油气勘探的精度要求[18-19]。为了弥补常规逆时偏移的缺陷,从最小二乘思想的角度重新审视逆时偏移的过程[20-25]。近年来,对最小二乘逆时偏移(LSRTM)的研究逐渐引起了国内外学者的关注。其中,为了提升效率,减少耗时,Dai等[26-27]引入了相位编码;针对高精度成像这一目标,黄建平等[28]优化了LSRTM算法;LSRTM还可应用于黏声介质,李振春等[29]完成了黏声介质LSRTM算法,改善成像质量。考虑到LSRTM的优势,国内外学者将弹性波逆时偏移(ERTM)与最小二乘偏移(LSM)结合发展出弹性波最小二乘逆时偏移(EL‐ SRTM)来进行复杂构造的保幅成像处理[30-31]。
从弹性波波动方程出发,推导弹性波线性正演以及逆时偏移算子表达式;然后从误差函数出发,推导Hessian矩阵的表达式,以此作为梯度预条件算子,实现一种预条件弹性波最小二乘逆时偏移(PELSRTM)算法;最后,通过SEG/EAGE二维盐丘模型测试P-ELSRTM对复杂模型成像时的有效性及其优势,并对比和分析常规ELSRTM及P-ELSRTM对低信噪比炮数据的适应性,以期证明预处理对反演过程中的改进效果。
1 弹性波最小二乘逆时偏移 1.1 线性化反偏移和偏移算子在非均匀各向同性弹性介质中,二维弹性波一阶速度-应力方程[30-31]可表示为
$ \left\{ {\begin{array}{*{20}{l}} {\rho \frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} t}} - \left( {\frac{{\partial {\kern 1pt} {\sigma _{{\rm{xx}}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} {\sigma _{{\rm{xz}}}}}}{{\partial {\kern 1pt} z}}} \right) = 0}\\ {\rho \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} t}} - \left( {\frac{{\partial {\kern 1pt} {\sigma _{{\rm{zz}}}}}}{{\partial {\kern 1pt} z}} + \frac{{\partial {\kern 1pt} {\sigma _{{\rm{xz}}}}}}{{\partial {\kern 1pt} x}}} \right) = 0}\\ {\frac{{\partial {\kern 1pt} {\sigma _{{\rm{xx}}}}}}{{\partial {\kern 1pt} t}} - \lambda \left( {\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}} \right) - 2\mu \frac{{\partial {u_{\rm{x}}}}}{{\partial x}} = {S_{{\rm{xx}}}}}\\ {\frac{{\partial {\kern 1pt} {\sigma _{{\rm{zz}}}}}}{{\partial {\kern 1pt} t}} - \lambda \left( {\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}} \right) - 2\mu \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}} = {S_{{\rm{zz}}}}}\\ {\frac{{\partial {\kern 1pt} {\sigma _{{\rm{xz}}}}}}{{\partial {\kern 1pt} t}} - \mu \left( {\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} z}} + \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} x}}} \right) = 0} \end{array}} \right. $ | (1) |
式中:Sxx和Szz分别为水平方向和垂直方向上的震源;ux和uz分别为介质中在水平方向和垂直方向上质点振动速度分量;σxx,σxz,σzz分别代表应力分量;ρ为介质质量密度,kg/cm3,常数;λ和μ均为拉梅常数。拉梅常数λ和μ、密度ρ、P波速度vP、S波速度vS之间的关系为:λ = ρ(vP2 - 2 vS2),μ = ρvS2。
与背景介质和扰动介质相对应的是背景波场P =[ux,uz,σxx,σzz,σxz]和扰动波场δu =[Δux,Δuz,Δσxx,Δ σzz,Δσxz]。总波场同样满足波动方程,即
$ \left\{ \begin{array}{l} \begin{array}{*{20}{c}} {\rho \frac{{\partial ({u_{\rm{x}}} + \Delta {u_{\rm{x}}})}}{{\partial {\kern 1pt} t}} - \left[ {\frac{{\partial ({\sigma _{{\rm{xx}}}} + \Delta {\sigma _{{\rm{xx}}}})}}{{\partial {\kern 1pt} x}} + \frac{{\partial ({\sigma _{{\rm{xz}}}} + \Delta {\sigma _{{\rm{xx}}}})}}{{\partial {\kern 1pt} z}}} \right] = 0}\\ {\rho \frac{{\partial ({u_{\rm{z}}} + \Delta {u_{\rm{z}}})}}{{\partial {\kern 1pt} t}} - \left[ {\frac{{\partial ({\sigma _{{\rm{zz}}}} + \Delta {\sigma _{{\rm{zz}}}})}}{{\partial {\kern 1pt} z}} + \frac{{\partial ({\sigma _{{\rm{xz}}}} + \Delta {\sigma _{{\rm{xz}}}})}}{{\partial {\kern 1pt} x}}} \right] = 0} \end{array}\\ \begin{array}{*{20}{l}} {\frac{{\partial ({\sigma _{{\rm{xx}}}} + \Delta {\sigma _{{\rm{xx}}}})}}{{\partial {\kern 1pt} t}} - (\lambda + \delta \lambda )\left[ {\frac{{\partial ({u_{\rm{x}}} + \Delta {u_{\rm{x}}})}}{{\partial {\kern 1pt} x}} + \frac{{\partial ({u_{\rm{z}}} + \Delta {u_{\rm{z}}})}}{{\partial {\kern 1pt} z}}} \right] - 2(\mu + \delta \mu )\frac{{\partial ({u_{\rm{x}}} + \Delta {u_{\rm{x}}})}}{{\partial {\kern 1pt} x}} = {S_{{\rm{xx}}}}}\\ {\frac{{\partial ({\sigma _{{\rm{zz}}}} + \Delta {\sigma _{{\rm{zz}}}})}}{{\partial {\kern 1pt} t}} - (\lambda + \delta \lambda )\left[ {\frac{{\partial ({u_{\rm{x}}} + \Delta {u_{\rm{x}}})}}{{\partial {\kern 1pt} x}} + \frac{{\partial ({u_{\rm{z}}} + \Delta {u_{\rm{z}}})}}{{\partial {\kern 1pt} z}}} \right] - 2(\mu + \delta \mu )\frac{{\partial ({u_{\rm{z}}} + \Delta {u_{\rm{z}}})}}{{\partial {\kern 1pt} z}} = {S_{{\rm{zz}}}}} \end{array}\\ \frac{{\partial ({\sigma _{{\rm{xz}}}} + \delta {\sigma _{{\rm{xz}}}})}}{{\partial t}} - (\mu + \delta \mu )\left[ {\frac{{\partial ({u_{\rm{x}}} + \delta {\kern 1pt} {u_{\rm{x}}})}}{{\partial {\kern 1pt} z}} + \frac{{\partial ({u_{\rm{z}}} + \delta {\kern 1pt} {u_{\rm{z}}})}}{{\partial {\kern 1pt} x}}} \right] = 0 \end{array} \right. $ | (2) |
将式(2)与式(1)相减,并且在模型扰动足够小的假设下,用背景波场P代替总波场Ptotal,可得到基于Born近似的关于扰动波场δu的波动方程,可以表示为
$ \left\{ \begin{array}{l} \begin{array}{*{20}{l}} {\rho \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} t}} - \left( {\frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{xx}}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{xz}}}}}}{{\partial {\kern 1pt} z}}} \right) = 0}\\ {\rho \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} t}} - \left( {\frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{zz}}}}}}{{\partial {\kern 1pt} z}} + \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{xz}}}}}}{{\partial {\kern 1pt} x}}} \right) = 0} \end{array}\\ \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{xx}}}}}}{{\partial {\kern 1pt} t}} - \lambda \left( {\frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}} \right) - 2\mu \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} = (\lambda + 2\mu ){m_{\rm{P}}}\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} + [(\lambda + 2\mu ){m_{\rm{P}}} - 2\mu {m_{\rm{S}}}]\frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}\\ \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{zz}}}}}}{{\partial {\kern 1pt} t}} - \lambda \left( {\frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}} \right) - 2\mu \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}} = (\lambda + 2\mu ){m_{\rm{P}}}\frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}} + [(\lambda + 2\mu ){m_{\rm{P}}} - 2\mu {m_{\rm{S}}}]\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}}\\ \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{xz}}}}}}{{\partial {\kern 1pt} t}} - \mu \left( {\frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} z}} + \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} x}}} \right) = \mu {m_{\rm{S}}}\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} z}} + \mu {m_{\rm{S}}}\frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} x}} \end{array} \right. $ | (3) |
式中:Δux和Δuz均为质点振动速度扰动项;Δσxx,Δσxz,Δσzz均为应力扰动项;P,S波反射系数模型可用速度扰动δvP,δvS和背景速度vP,vS表示为mP (x) = 2 δvP/vP,mS (x) = 2 δvS/vS,反映了介质速度扰动的相对大小。
从式(3)可以看出,要得到质点振动速度扰动项和应力扰动项需要经历2个过程:第1个过程是先通过式(1)得到背景波场P0;第2个过程是利用式(3)进行计算,通过背景慢度模型和新的虚拟震源项再次激发,再传播到地下x位置处形成散射波场。另外,从式(3)的右端表达式可以看出,质点振动速度扰动项和应力扰动项的产生是由背景介质参数与介质速度扰动作为新的震源,而且两者与介质速度扰动项以及应力扰动项呈线性关系。
$ L{\kern 1pt} m = d $ | (4) |
式中:m为反射系数模型的矩阵形式;d表示矩阵形式的观测数据,d =(dx,dz,0,0,0)T;L为基于Born近似下的线性反偏移算子,它建立了介质参数与正演数据之间的联系。
采用时间2阶、空间10阶精度的交错网格有限差分法对式(1)进行差分离散。在进行波场模拟时,由于人工截断边界会导致边界反射,采用完全匹配层(PML)吸收边界,来压制边界反射提高模拟区域的精度与信噪比。
本文借助伴随状态法得到伴随算子及相应的伴随波场控制方程,由伴随状态法可得
$ \langle {L^*}w,\delta u\rangle = \langle w,L{\kern 1pt} \delta {\kern 1pt} u\rangle $ | (5) |
式中:w为伴随波场向量,
则伴随方程可表示为
$ \left\{ \begin{array}{l} \begin{array}{*{20}{c}} {\rho \frac{{\partial {\kern 1pt} {w_{\rm{x}}}}}{{\partial {\kern 1pt} t}} - \frac{{\partial [(\lambda + 2\mu ){{\hat \sigma }_{{\rm{xx}}}}]}}{{\partial {\kern 1pt} x}} - \frac{{\partial (\lambda {\kern 1pt} {{\hat \sigma }_{{\rm{zz}}}})}}{{\partial {\kern 1pt} x}} - \frac{{\partial (\mu {\kern 1pt} {{\hat \sigma }_{{\rm{xz}}}})}}{{\partial {\kern 1pt} z}} = - \delta {d_{\rm{x}}}({x_{\rm{g}}},t;{x_{\rm{s}}})}\\ {\rho \frac{{\partial {\kern 1pt} {w_{\rm{z}}}}}{{\partial {\kern 1pt} t}} - \frac{{\partial [(\lambda + 2\mu ){{\hat \sigma }_{{\rm{zz}}}}]}}{{\partial {\kern 1pt} z}} - \frac{{\partial (\lambda {\kern 1pt} {{\hat \sigma }_{{\rm{xx}}}})}}{{\partial {\kern 1pt} z}} - \frac{{\partial (\mu {\kern 1pt} {{\hat \sigma }_{{\rm{xz}}}})}}{{\partial {\kern 1pt} x}} = - \delta {d_{\rm{z}}}({x_{\rm{g}}},t;{x_{\rm{s}}})} \end{array}\\ \begin{array}{*{20}{l}} {\frac{{\partial {\kern 1pt} {\kern 1pt} {{\hat \sigma }_{{\rm{xx}}}}}}{{\partial {\kern 1pt} t}} - \frac{{\partial {\kern 1pt} {w_{\rm{x}}}}}{{\partial {\kern 1pt} x}} = 0}\\ {\frac{{\partial {\kern 1pt} {\kern 1pt} {{\hat \sigma }_{{\rm{zz}}}}}}{{\partial {\kern 1pt} t}} - \frac{{\partial {\kern 1pt} {w_{\rm{z}}}}}{{\partial {\kern 1pt} z}} = 0}\\ {\frac{{\partial {\kern 1pt} {\kern 1pt} {{\hat \sigma }_{{\rm{xz}}}}}}{{\partial {\kern 1pt} t}} - \frac{{\partial {\kern 1pt} {w_{\rm{x}}}}}{{\partial {\kern 1pt} z}} - \frac{{\partial {\kern 1pt} {w_{\rm{z}}}}}{{\partial {\kern 1pt} x}} = 0} \end{array} \end{array} \right. $ | (6) |
式中:δd为残差向量,
为了得到与记录数据最佳匹配的偏移结果,把最小二乘逆时偏移引入反演中,其误差函数定义为
$ \chi (m) = \frac{1}{2}{\left\| {Lm - {d_{{\rm{obs}}}}} \right\|^2} = \frac{1}{2}{\left\| {\delta d} \right\|^2} $ | (7) |
根据伴随状态方法,求解式(7)可得λ,μ的梯度表达式为
$ \left\{ {\begin{array}{*{20}{l}} {\delta {\kern 1pt} \lambda = \frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} \lambda }} = - \int {\left[ {\left( {\frac{{\partial {u_{\rm{x}}}}}{{\partial x}} + \frac{{\partial {u_{\rm{z}}}}}{{\partial z}}} \right)({{\hat \sigma }_{{\rm{xx}}}} + {{\hat \sigma }_{{\rm{zz}}}})} \right]{\rm{d}}{\kern 1pt} t} }\\ {\delta {\kern 1pt} \mu = \frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} \mu }} = - \int {\left[ {2\frac{{\partial {u_{\rm{x}}}}}{{\partial x}}{{\hat \sigma }_{{\rm{xx}}}} + 2\frac{{\partial {u_{\rm{z}}}}}{{\partial z}}{{\hat \sigma }_{{\rm{zz}}}} + \left( {\frac{{\partial {u_{\rm{x}}}}}{{\partial z}} + \frac{{\partial {u_{\rm{z}}}}}{{\partial x}}} \right){{\hat \sigma }_{{\rm{xz}}}}} \right]{\rm{d}}{\kern 1pt} t} } \end{array}} \right. $ | (8) |
vP和vS的梯度可由链式公式求得,即
$ \left\{ {\begin{array}{*{20}{l}} {\delta {v_{\rm{P}}} = 2{v_{\rm{P}}}\rho \frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} \lambda }}}\\ {\delta {v_{\rm{S}}} = - 4{v_{\rm{S}}}\rho \frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} \lambda }} + 2{v_{\rm{S}}}\rho \frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} \mu }}} \end{array}} \right. $ | (9) |
梯度计算是最小二乘成像的核心问题之一,对梯度做预条件,加速问题的收敛。本文中用于后续迭代更新的梯度为g = (δvP,δvS)。将式(9)代入式(7)可得δvP和δvS的梯度表达式,即
$ \left\{ {\begin{array}{*{20}{l}} {\delta {v_{\rm{P}}} = - 2{v_{\rm{P}}}\rho \int {\left[ {\left( {\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}} \right)({{\hat \sigma }_{{\rm{xx}}}} + {{\hat \sigma }_{{\rm{zz}}}})} \right]{\rm{d}}{\kern 1pt} t} }\\ {\delta {v_{\rm{S}}} = - 2{v_{\rm{S}}}\rho \int {\left[ { - 2\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}}{{\hat \sigma }_{{\rm{xx}}}} - 2\frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}{{\hat \sigma }_{{\rm{zz}}}} + \left( {\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} z}} + \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} x}}} \right){{\hat \sigma }_{{\rm{xz}}}}} \right]{\rm{d}}{\kern 1pt} t} } \end{array}} \right. $ | (10) |
式(7)所对应的解的估计称最小二乘偏移成像的结果,真实反射模型的最小二乘解:
$ m = {({L^*}L)^{ - 1}}{L^*}{d_{{\rm{obs}}}} $ | (11) |
传统的偏移方法是通过反偏移算子的共轭转置实现的,会造成地震信息缺失。从式(9)可看出,最小二乘偏移反演通过Hessian算子L* L对偏移成像进行修正,难点就在于求解一个精确的Hessian并求逆。Hessian算子H是误差函数在m处对模型的二阶导数,即
$ \begin{array}{*{20}{l}} {H = \frac{{{\partial ^2}\chi }}{{\partial {\kern 1pt} {m^2}}} = \frac{{{\partial ^2}\chi }}{{\partial {\kern 1pt} {m^2}}}(d - {d_{{\rm{obs}}}}) + \frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} m}}\frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} m}}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{{\partial {\kern 1pt} {L^*}}}{{\partial {\kern 1pt} m}}(d - {d_{{\rm{obs}}}}) + {L^*}L} \end{array} $ | (12) |
从式(12)可以看出,H由2部分构成,其中,
式(1)用矩阵形式可表示为
$ S{\kern 1pt} w = f $ | (13) |
其中:
$ \begin{array}{l} S = \left[ {\begin{array}{*{20}{c}} {\rho \frac{\partial }{{\partial {\kern 1pt} t}}}&0&{ - \frac{\partial }{{\partial {\kern 1pt} x}}}&0&{ - \frac{\partial }{{\partial {\kern 1pt} z}}}\\ 0&{\rho \frac{\partial }{{\partial {\kern 1pt} t}}}&0&{ - \frac{\partial }{{\partial {\kern 1pt} z}}}&{ - \frac{\partial }{{\partial {\kern 1pt} x}}}\\ { - (\lambda + 2\mu )\frac{\partial }{{\partial {\kern 1pt} x}}}&{ - \lambda \frac{\partial }{{\partial {\kern 1pt} z}}}&{\frac{\partial }{{\partial {\kern 1pt} t}}}&0&0\\ { - \lambda \frac{\partial }{{\partial {\kern 1pt} x}}}&{ - (\lambda + 2\mu )\frac{\partial }{{\partial {\kern 1pt} z}}}&0&{\frac{\partial }{{\partial t}}}&0\\ { - \mu \frac{\partial }{{\partial {\kern 1pt} z}}}&{ - \mu \frac{\partial }{{\partial {\kern 1pt} x}}}&0&0&{\frac{\partial }{{\partial {\kern 1pt} t}}} \end{array}} \right],\\ w = \left[ {\begin{array}{*{20}{l}} {{u_{\rm{x}}}}\\ {{u_{\rm{z}}}}\\ {{\sigma _{{\rm{xx}}}}}\\ {{\sigma _{{\rm{zz}}}}}\\ {{\sigma _{{\rm{xz}}}}} \end{array}} \right],f = \left[ {\begin{array}{*{20}{l}} 0\\ 0\\ {{S_{{\rm{xx}}}}}\\ {{S_{{\rm{zz}}}}}\\ 0 \end{array}} \right]。\end{array} $ |
对式(13)求导可得
$ {\frac{{\partial {\kern 1pt} S}}{{\partial {\kern 1pt} m}}w + S\frac{{\partial {\kern 1pt} w}}{{\partial {\kern 1pt} m}} = 0} $ | (14) |
移项:
$ {\frac{{\partial {\kern 1pt} w}}{{\partial {\kern 1pt} m}} = - {S^{ - 1}}\frac{{\partial {\kern 1pt} S}}{{\partial {\kern 1pt} m}}w} $ | (15) |
通过汉森矩阵的线性项来近似代替全汉森矩阵,即
$ H = {L^*}L = {\left( {\frac{{\partial {\kern 1pt} w}}{{\partial {\kern 1pt} m}}} \right)^*}\left( {\frac{{\partial {\kern 1pt} w}}{{\partial {\kern 1pt} m}}} \right) = {\left( {{S^{ - 1}}\frac{{\partial {\kern 1pt} S}}{{\partial {\kern 1pt} m}}{\kern 1pt} {\kern 1pt} w} \right)^*}\left( {{S^{ - 1}}\frac{{\partial {\kern 1pt} S}}{{\partial {\kern 1pt} m}}w} \right) $ | (16) |
用拉梅常数表示的预条件算子表达式为
$ \left\{ {\begin{array}{*{20}{l}} {{H^\lambda } = \int {\frac{{{{({{\dot \sigma }_{{\rm{xx}}}} + {{\dot \sigma }_{{\rm{zz}}}})}^2}}}{{2{{(\lambda + \mu )}^2}}}{\rm{d}}{\kern 1pt} t} }\\ {{H^\mu } = \int {\left[ {\frac{{\dot \sigma _{{\rm{xz}}}^2}}{{{\mu ^2}}} + \frac{{{{({{\dot \sigma }_{{\rm{xx}}}} + {{\dot \sigma }_{{\rm{zz}}}})}^2}}}{{2{{(\lambda + \mu )}^2}}} + \frac{{{{({{\dot \sigma }_{{\rm{xx}}}} - {{\dot \sigma }_{{\rm{zz}}}})}^2}}}{{2{\mu ^2}}}} \right]{\rm{d}}{\kern 1pt} t} } \end{array}} \right. $ | (17) |
利用共轭梯度法迭代求解式(7)得到模拟数据与观测数据的差别最小的最佳反射系数模型,以P表示预条件算子,梯度g = (δvP, δvS),具体实现流程如下:
(1) 输入初始模型(初始模型赋0值,即第1次迭代的结果为逆时偏移的偏移剖面),观测数据dobs,梯度误差阈值δ。
(2) 计算第i次迭代共轭梯度dki和共轭方向修正因子β
$ \left\{ {\begin{array}{*{20}{l}} {d{\kern 1pt} {k^i} = P{\kern 1pt} {g^i} + \beta {\kern 1pt} d{\kern 1pt} {k^{i - 1}}}\\ {\beta = \frac{{{g^i}P{\kern 1pt} {g^i}}}{{{g^{i - 1}}P{\kern 1pt} {g^{i - 1}}}}} \end{array}} \right. $ | (18) |
(3) 计算步长α并对成像结果进行更新,得到mi + 1
$ \left\{ {\begin{array}{*{20}{l}} {\alpha = \frac{{{{(d{\kern 1pt} {k^i})}^T}{g^i}}}{{{{(L{\kern 1pt} d{\kern 1pt} {k^i})}^T}L{\kern 1pt} d{\kern 1pt} {k^i}}}}\\ {{m^{i + 1}} = {m^i} - \alpha {\kern 1pt} d{\kern 1pt} {k^i}} \end{array}} \right. $ | (19) |
(4) 判断是否满足终止迭代条件,若满足,则退出迭代并保存成像结果;不满足,则回到第(2),(3)步继续迭代直到满足迭代终止条件,最后一次迭代即为最终的成像结果。
2 模型试算 2.1 盐丘模型采用SEG/EAGE提供的盐丘模型来测试常规ELSRTM和P-ELSRTM对模型的成像效果,纵横波真实速度场和反射系数模型如图 1所示。设定纵横波速度比为1.73,密度为常数。模型横向网格点数为645,纵向网格点数为150,采样间隔为10 m。利用该模型主要考察P-ELSRTM成像算法是否能够使盐下构造准确成像以及对振幅的补偿效果:由于盐丘的屏蔽作用,导致盐下构造的振幅较小,另外,盐丘侧翼发育的高陡小构造会引起严重的构造假象。观测系统采用全接收方式,每炮的网格间隔数为10 m,共60炮。震源类型为爆炸震源,震源函数采用主频为25 Hz的雷克子波。
![]() |
下载原图 图 1 SEG/EAGE二维盐丘模型 Fig. 1 SEG-EAGE 2-D salt dome model |
ERTM成像结果如图 2所示,逆时偏移成像中使用全波动方程模拟波场,可以实现对陡倾角构造的成像,但是在ERTM剖面中含有较强的低波数噪声和震源效应[图 2(a)和图 2(b)]。应用拉普拉斯算子去噪之后,相比未作滤波的ERTM结果得到了改善[图 2(c)和图 2(d)],模型浅部的低频噪音得到了压制,但成像能量不均衡。使用预条件算子和滤波之后,相比滤波后的ERTM成像剖面质量得到进一步提高[图 2(e)和图 2(f)],盐下构造清晰,且模型浅部的低频噪音和震源效应也得到了更好地压制。
![]() |
下载原图 图 2 ERTM成像结果对比 Fig. 2 Comparison of ERTM imaging results |
为了讨论预条件ELSRTM的有效性及优点,分别对常规ELSRTM和预条件P-ELSRTM成像结果进行对比(图 3):①未作滤波的常规ELSRTM在迭代次数为5和30次的成像结果中都残留较严重的低波数噪音,影响对盐丘之上沉积层的识别[图 3 (a)-(d)],而应用高通滤波预条件算子后,成像剖面中的低频噪音得到了压制,剖面振幅更加均衡[图 3(e)-(h)],证明了滤波算子的有效性;②滤波后的常规ELSRTM成像结果中盐体边界清晰,盐下成像相对于ERTM结果更为准确,但仍存在深部成像能量不足的现象,并且对顶部低频噪声的压制不彻底[图 3(e)-(h)];③应用照明补偿预条件算子和高通滤波后的结果[图 3(i)-(l)]表明,采用预条件成像方法实现了对复杂模型深部小断层较为清晰的成像,盐下断层及基底均能正确成像,并且成像剖面信噪比更高、能量更加均衡。
![]() |
下载原图 图 3 ELSRTM成像结果对比 Fig. 3 Comparison of ELSRTM imaging results |
分别从常规ELSRTM和P-ELSRTM的成像结果中抽取偏移距为4 900 m的地震单道以对比各方法的保幅性(迭代30次),结果如图 4所示。从图 4可看出:①在模型的浅部区域由于低频噪音和震源效应的存在,常规ELSRTM的振幅曲线会有明显跳动,而P-ELSRTM的振幅曲线抖动幅度较小;②在模型的500 m,830 m,1 430 m,4 500 m(反射界面)处,常规ELSRTM的振幅与真实反射系数之间有较大的差异,而P-ELSRTM的振幅与真实反射系数更加接近。因此,在浅-中深层P-ELSRTM都比常规ELSRTM具有更好的保幅性。
![]() |
下载原图 图 4 ELSRTM与P-ELSRTM成像结果单道振幅对比(第1次及第30次迭代) Fig. 4 Amplitude of single trace of ELSRTM and P-ELSRTM imaging results |
常规ELSRTM和P-ELSRTM的成像结果与真实结果的归一化残差随迭代次数变化的结果如图 5所示。从图 5可以看出:P-ELSRTM的成像反演结果在前60次迭代中都能稳定收敛,最终收敛到22.6%;常规ELSRTM随着迭代次数的增加只能收敛到52%便停止收敛。因此,使用梯度预条件算子能够使复杂模型成像结果的残差收敛更快且达到更低的值。
![]() |
下载原图 图 5 P-ELSRTM与ELSRTM误差随迭代次数的变化图(迭代60次) Fig. 5 Residual convergence curves of P-ELSRTM and ELSRTM |
在这次测试中,CPU型号为Intel(R)Xeon(R) CPU E5-2650 v2@2.60 GHz,P-ELSRTM和常规EL‐ SRTM均串行运行,迭代一次所花费的时间分别为4 667 s和4 642 s。基于不同最优化算法的LSRTM结果的对比如表 1所列。由于计算汉森算子所花费的计算量相对于偏移运算较小,因此表 1中忽略了汉森算子的计算量,仅对比了偏移运算的计算量。P-ELSRTM算法与ELSRTM算法的计算量相近,但前者具有更高的成像质量、更快的收敛速度及更好的稳定性。
![]() |
下载CSV 表 1 ELSRTM成像结果对比 Table 1 Comparison of ELSRTM imaging results |
野外地震资料采集的数据中不可避免地含有环境以及人为因素造成的随机噪声。为了测试PELSRTM对低信噪比数据的适应性,对炮记录添加低信噪比的随机噪声,使其X分量和Z分量信噪比为0.5 dB(图 6)。从ELSRTM和P-ELSRTM的成像结果(图 7)可看出,对于低信噪比数据,ELSRTM和P-ELSRTM结果中都含有类似于随机噪声的偏移假象,但相比常规ELSRTM,P-ELSRTM能够更好地压制低频噪音。从残差收敛曲线可以看出,P-ELSRTM的收敛速度较常规ELSRTM快,并且迭代60次后其残差较小(图 8)。因此,P-ELSRTM在处理含噪数据方面有较大优势,抗噪性更好。
![]() |
下载原图 图 6 低信噪比单炮数据 Fig. 6 Low SNR shot data |
![]() |
下载原图 图 7 低信噪比数据成像结果 Fig. 7 Imaging results of low SNR data |
![]() |
下载原图 图 8 含噪音数据的P-ELSRTM与ELSRTM误差随迭代次数的变化图(迭代60次) Fig. 8 Residual convergence curves of P-ELSRTM and ELSRTM with noise |
(1) P-ELSRTM和ELSRTM能够有效压制浅层低频噪音,消除震源效应,有效改善常规ERTM的成像质量,具有更高的成像分辨率。
(2) P-ELSRTM在中深部具有更好的保幅性和分辨率,有利于深部及盐下构造的识别;P-ELSRTM比常规ELSRTM的残差曲线更快地收敛。在对含噪数据的成像测试中,P-ELSRTM比常规ELSRTM可以收敛到一个更低的残差水平。
(3) 本文使用一阶位移波动方程得到的纵横波是耦合在一起的,因此通过迭代无法完全消除纵横波串扰噪音。针对该问题,进行纵波和横波波场矢量波场分离可进一步提高成像精度。
[1] |
TENG Y C, DAI T F. Finite-element prestack reverse-time migration for elastic waves. Geophysical Prospecting for Petroleum, 1988, 5(1): 1204-1208. |
[2] |
SUN R, WANG A. Scalar reverse-time depth migration of prestack elastic seismic data. Geophysics, 2001, 66(5): 1519. DOI:10.1190/1.1487098 |
[3] |
DUAN Y, SAVA P. Scalar imaging condition for elastic reverse time migration. Geophysics, 2015, 80(4): S127-S136. |
[4] |
刘红伟, 李博, 刘洪, 等. 地震叠前逆时偏移高阶有限差分算法及GPU实现. 地球物理学报, 2010, 53(7): 1725-1733. LIU H W, LI B, LIU H, et al. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation. Chinese Journal of Geophysics, 2010, 53(7): 1725-1733. DOI:10.3969/j.issn.0001-5733.2010.07.024 |
[5] |
徐兴荣, 王西文, 王宇超, 等. 基于波场分离理论的逆时偏移成像条件研究及应用. 地球物理学进展, 2012, 27(5): 2084-2090. XU X R, WANG X W, WANG Y C, et al. Study and application of imaging condition for reverse-time migration based on wave-fields separation. Progress in Geophysics, 2012, 27(5): 2084-2090. |
[6] |
陈可洋. 逆时成像技术在大庆探区复杂构造成像中的应用. 岩性油气藏, 2017, 29(6): 91-100. CHEN K Y. Application of reverse-time migration technology to complex structural imaging in Daqing exploration area. Lithologic Reservoirs, 2017, 29(6): 91-100. DOI:10.3969/j.issn.1673-8926.2017.06.012 |
[7] |
陈可洋. 几种地震观测方式的逆时成像分析. 岩性油气藏, 2013, 25(1): 95-101. CHEN K Y. Reverse-time migration analysis of several seismic observation models. Lithologic Reservoirs, 2013, 25(1): 95-101. DOI:10.3969/j.issn.1673-8926.2013.01.018 |
[8] |
CHANG W F, MCMECHAN G A. Reverse-time migration of offset vertical seismic profiling data using the excitation-time imaging condition. Geophysics, 1986, 51(1): 139-140. |
[9] |
CHANG W F, MCMECHAN G A. Elastic reverse-time migration. Geophysics, 1987, 52(3): 243-256. |
[10] |
CHANG W F, MCMECHAN G A. 3-D elastic prestack, reversetime depth migration. Geophysics, 1994, 59(4): 597-609. |
[11] |
董良国, 郭晓玲, 吴晓丰, 等. 起伏地表弹性波传播有限差分法数值模拟. 天然气工业, 2007, 27(10): 38-41. DONG L G, GUO X L, WU X F, et al. Finite difference numerical simulation for the elastic wave propagation in rugged topography. Natural Gas Industry, 2007, 27(10): 38-41. DOI:10.3321/j.issn:1000-0976.2007.10.010 |
[12] |
陈可洋. 各向异性弹性介质方向行波波场分离正演数值模拟. 岩性油气藏, 2014, 26(5): 91-96. CHEN K Y. Wave field separating numerical simulation of anisotropic elastic medium directional one-way wave. Lithologic Reservoirs, 2014, 26(5): 91-96. DOI:10.3969/j.issn.1673-8926.2014.05.017 |
[13] |
王维红, 张伟, 石颖, 等. 基于波场分离的弹性波逆时偏移. 地球物理学报, 2017, 60(7): 2813-2824. WANG W H, ZHANG W, SHI Y, et al. Elastic reverse time migration based on wavefield separation. Chinese Journal of Geophysics, 2017, 60(7): 2813-2824. |
[14] |
张伟, 石颖. 矢量分离纵横波场的弹性波逆时偏移. 地球物理学进展, 2017, 32(4): 1728-1734. ZHANG W, SHI Y. Elastic reverse time migration based on vector decomposition of P-and S-wavefields. Progress in Geophysics, 2017, 32(4): 1728-1734. |
[15] |
杜启振, 秦童. 横向各向同性介质弹性波多分量叠前逆时偏移. 地球物理学报, 2009, 52(3): 801-807. DU Q Z, QIN T. Multicomponent prestack reverse-time migration of elastic waves in transverse isotropic medium. Chinese Journal of Geophysics, 2009, 52(3): 801-807. |
[16] |
张晓语, 杜启振, 张树奎, 等. 基于一阶弹性波方程的能量互相关成像条件. 地球物理学报, 2019, 62(1): 289-297. ZHANG X Y, DU Q Z, ZHANG S K, et al. The energy crosscorrelation imaging based on first-order elastic wave equations. Chinese Journal of Geophysics, 2019, 62(1): 289-297. |
[17] |
张智, 刘有山, 徐涛, 等. 弹性波逆时偏移中的稳定激发振幅成像条件. 地球物理学报, 2013, 56(10): 3523-3533. ZHANG Z, LIU Y S, XU T, et al. A stable excitation amplitude imaging condition for reverse time migration in elastic wave equation. Chinese Journal of Geophysics, 2013, 56(10): 3523-3533. DOI:10.6038/cjg20131027 |
[18] |
王华忠, 王雄文, 王西文. 地震波反演的基本问题分析. 岩性油气藏, 2012, 24(6): 1-9. WANG H Z, WANG X W, WANG X W. Analysis of the basic problems of seismic wave inversion. Lithologic Reservoirs, 2012, 24(6): 1-9. DOI:10.3969/j.issn.1673-8926.2012.06.002 |
[19] |
任浩然, 王华忠, 黄光辉. 地震波反演成像方法的理论分析与对比. 岩性油气藏, 2012, 24(5): 12-18. REN H R, WANG H Z, HUANG G H. Theoretical analysis and comparison of seismic wave inversion and imaging methods. Lithologic Reservoirs, 2012, 24(5): 12-18. DOI:10.3969/j.issn.1673-8926.2012.05.002 |
[20] |
李庆洋, 黄建平, 李振春. 基于Student's t分布的不依赖子波最小二乘逆时偏移. 地球物理学报, 2017, 60(12): 4790-4800. LI Q Y, HUANG J P, LI Z C. Source-independent least-squares reverse time migration using student's t distribution. Chinese Journal of Geophysics, 2017, 60(12): 4790-4800. DOI:10.6038/cjg20171220 |
[21] |
NEMETH T, WU C, SCHUSTER G T. Least-squares migration of incomplete reflection data. Geophysics, 1999, 64(1): 208-221. |
[22] |
LI C, HUANG J P, LI Z C, et al. Regularized least-squares migration of simultaneous-source seismic data with adaptive singular spectrum analysis. Petroleum Science, 2017, 1(14): 61-74. |
[23] |
李振春, 李闯, 黄建平, 等. 基于先验模型约束的最小二乘逆时偏移方法. 石油地球物理勘探, 2016, 51(4): 738-744. LI Z C, LI C, HUANG J P, et al. Regularized least-squares reverse time migration with prior model. Oil Geophysical Prospecting, 2016, 51(4): 738-744. |
[24] |
李闯, 黄建平, 李振春, 等. 预条件最小二乘逆时偏移方法. 石油地球物理勘探, 2016, 51(3): 513-520. LI C, HUANG J P, LI Z C, et al. Preconditioned least-squares reverse time migration. Oil Geophysical Prospecting, 2016, 51(3): 513-520. |
[25] |
刘梦丽, 黄建平, 李闯, 等. 基于角度滤波成像的最小二乘逆时偏移. 石油地球物理勘探, 2018, 53(3): 469-477. LIU M L, HUANG J P, LI C, et al. Least-squares reverse time migration based on the angular filtering imaging condition. Oil Geophysical Prospecting, 2018, 53(3): 469-477. |
[26] |
DAI W, BOONYASIRIWAT C, SCHUSTER G T. Multi-source least-squares reverse time migration. Geophysical Prospecting, 2012, 60(4): 681-695. DOI:10.1111/j.1365-2478.2012.01092.x |
[27] |
DAI W, SCHUSTER G T. Plane-wave least-squares reversetime migration. Geophysics, 2013, 78(4): S165-S177. |
[28] |
黄建平, 曹晓莉, 李振春, 等. 最小二乘逆时偏移在近地表高精度成像中的应用. 石油地球物理勘探, 2014, 49(1): 107-112. HUANG J P, CAO X L, LI Z C, et al. Least square reverse time migration in high resolution imaging of near surface. Oil Geophysical Prospecting, 2014, 49(1): 107-112. |
[29] |
李振春, 郭振波, 田坤. 黏声介质最小平方逆时偏移. 地球物理学报, 2014, 57(1): 214-228. LI Z C, GUO Z B, TIAN K. Least-squares reverse time migration in visco-acoustic media. Chinese Journal of Geophysics, 2014, 57(1): 214-228. |
[30] |
GU B L, LI Z C, YANG P, et al. Elastic least-squares reverse time migration with hybrid l1/l2 misfit function. Geophysics, 2017, 82(3): S271-S291. |
[31] |
REN Z M, LIU Y, SEN M K. Least-squares reverse time migration in elastic media. Geophysical Journal International, 2017, 208(2): 1103-1125. DOI:10.1093/gji/ggw443 |
[32] |
LI C, GAO J, GAO Z, et al. Periodic plane-wave least-squares reverse time migration for diffractions. Geophysics, 2020, 85(4): S185-S198. |
[33] |
LI C, GAO J, WANG R, et al. Enhancing subsurface scatters using reflection-damped plane-wave least-squares reverse time migration. IEEE Geoence and Remote Sensing Letters, 2020, 17(4): 706-710. |
[34] |
YANG J D, ZHU H J. Viscoacoustic least-squares reverse-time migration using a time-domain complex-valued wave equation. Geophysics, 2019, 84(5): 1-130. |
[35] |
郭旭, 黄建平, 李振春, 等. 基于一阶速度-应力方程的VTI介质最小二乘逆时偏移. 地球物理学报, 2019, 62(6): 2188-2202. GUO X, HUANG J P, LI Z C, et al. Least-squares reverse time migration based on first-order velocity-stress wave equation in VTI media. Chinese Journal of Geophysics, 2019, 62(6): 2188-2202. |