2019, Vol. 62
2. 海洋国家实验室海洋矿产资源评价与探测技术功能实验室, 山东青岛 266071;
3. 贵州大学资源与环境工程学院, 贵阳 550025;
4. 中国石化中原物探研究院, 河南濮阳 457001
, HUANG JianPing1,2
, LI ZhenChun1,2, HUANG JinQiang3, LI QingYang4, ZHU Feng1,2, LIU MengLi1,2
2. Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao Shandong 266071, China;
3. Guizhou University Resource and Environmental Engineering College, Guiyang Guizhou 550025, China;
4. Sinopec Geophysical Prospecting Institute of Zhongyuan, Puyang Henan 457001, China
大量的生产实例表明，实际沉积地震地层大多会呈现地震各向异性特征，至少是弱各向异性(Tsvankin and Thomsen, 1994).在地震勘探中，备受学者关注的各向异性介质类型主要有：(1)具有垂直对称轴的横向各向同性(VTI，Transversely Isotropy with a Vertical axis of symmetry)介质；(2)具有水平对称轴的横向各向同性(HTI，Transversely Isotropy with a Horizontal axis of symmetry)介质；(3)具有倾斜对称轴的横向各向同性(TTI，Transversely Isotropy with a Tilted axis of symmetry)介质；(4)具有三个相互正交对称面的正交各向异性(OA，Orthorhombic Anisotropy)介质(Tatham and McCormack, 1991).在常规地震资料处理中，一般假设地下介质为各向同性，对于地下大尺度构造，这种近似的影响是可以忽略的，然而，随着油气勘探的深入，研究重心逐渐转向小尺度目标，如果忽略介质的各向异性会导致速度提取不准确，特别是造成目标区域成像不准确，进而影响储层预测、油藏描述、复杂构造成像以及岩性成像等等.为此，随着对地震资料分辨率以及成像质量要求的提高，地下介质的各向异性在地震资料处理中的影响已经不容忽视(吴国忱，2006).
近年来，许多学者围绕各向异性拟声波正演模拟和偏移成像技术展开了大量的研究工作.各向异性介质弹性波波动方程可以准确描述地震波在地下介质中的传播规律，然而在各向异性介质中，纵横波在传播过程中相互耦合，很难进行分离成像.鉴于此，Alkhalifah(1998)提出TI介质声学近似，即将沿对称轴方向上的横波速度设为零，以达到简化波动方程的目的.声学近似这一简化手段，促进了TI介质qP波简化方程构建以及逆时偏移技术的快速发展.根据qP波方程建立机制，可划分为如下两类：(1)完全解耦的纯qP波方程.针对解耦这一问题，Klié和Toro(2001)较早开始研究，并推导了一个近似的VTI介质声波方程；之后，Du等(2005)和Zhang等(2005)基于弱各向异性假设条件和平方根近似，推导出时间-波数域qP波解耦方程；Pestana等(2011)和Zhan等(2011)基于快速展开法实现了TI介质中纯qP波方程的解耦.近来，Xu和Zhou(2014a, 2014b)采用微分算子分解策略，成功模拟出纯qP波.尽管此类解耦方程在运动学上是等价的，但其动力学特性难以保证，且其波场的数量值大小也迥然不同.(2)qP波-qSV波耦合的方程.Alkhalifah(2000)基于声学近似假设，推导出了VTI介质四阶波动方程，由于其包含高阶项，Zhou等(2006a)通过引入辅助函数对其进行降阶处理，推导出了VTI介质的二阶耦合的偏微分方程.但对于TTI介质(Zhou et al., 2006b)，在倾角剧烈变化的区域，会存在严重的不稳定现象.针对不稳定这一问题，虽然一些学者提出通过引入非零横波项进行修正以及将不稳定区域设为椭圆各向异性等方法进行处理(Fletcher et al., 2009)，但却无法克服方程本身的不稳定性.Duveneck等(2008)、Zhang和Zhang(2009)、Zhang等(2011)从弹性波的基本理论(胡克定律和运动方程)出发，推导出一套相对稳定的VTI拟声波波动方程.这些研究进程成为各向异性偏移成像技术快速发展的重要组成部分.但从总体上看，利用反演策略来处理各向异性介质中的偏移成像问题仍属少见.
相比于射线类成像方法(如Krichhoff偏移、高斯束偏移)和单程波成像方法，逆时偏移(Reverse Time Migration，RTM)成像在成像精度和对复杂构造的适应力上具有明显优势，成为目前地震成像的重要手段.然而常规逆时偏移，往往利用正向传播的震源波场和反向传播的检波波场做互相关进行成像，其偏移算子仅仅是正演算子的共轭转置，而不是它的逆(Claerbout，1992).由于采集孔径、复杂地下构造以及有限波场带宽的影响，常规偏移方法只能对地下构造模糊成像，通常只能提供较准确的地下构造信息，无法满足岩性油气藏勘探开发的需求.最小二乘偏移(Least-Squares Migration，LSM)采用模型匹配数据的思想，将成像视为最小二乘意义下的反演问题，通过共轭梯度法(或最速下降法、高斯牛顿法等)迭代使误差函数达到最小，得到分辨率更高、振幅保真性更好的成像结果.最小二乘反演思想是由Tarantola(1984)提出的，随后，围绕最小二乘偏移技术，国内外学者进行了大量的研究工作.最小二乘偏移思想早期主要被应用于射线理论(Nemeth et al., 1999；Dai et al., 2011；黄建平等，2013；刘玉金等，2013)和单程波理论(Küehl and Sacchi, 2003；Clapp et al., 2005；周华敏等，2014；黄建平等，2016a；Zhu et al., 2017)；随后基于双程波理论的最小二乘逆时偏移(Least-Squares Reverse Time Migration，LSRTM)方法逐渐发展起来(Dai et al., 2012；王彦飞，2013；黄建平等，2014；李振春等，2014；黄建平等，2015；李庆洋等，2016a；李庆洋等，2016b；李振春等，2017).
当前LSRTM理论算法大多是基于二阶常密度标量声波方程(LSRTM based on Second-order，S-LSRTM)建立的.而地下实际地层是变密度的，基于波形匹配策略的常规LSRTM算法无法考虑密度变化对振幅的影响，进而无法获得相对保真的成像结果，失去了LSRTM的重要优势.一阶速度-应力方程可以方便地处理变密度介质，为此，李庆洋等(2016c)发展了一套基于一阶速度-应力方程的各向同性LSRTM(LSRTM based on First-order, F-LSRTM)算法，并取得了较好地应用效果.然而，实际地下介质不仅是变密度的，而且大多存在各向异性，对于各向异性地区，如果仅仅采用各向同性LSRTM算法，将会导致模拟数据与实际观测数据不匹配，进而可能会出现偏移噪声严重、成像位置难以准确识别、迭代不收敛等问题，影响反演成像质量.VTI介质一阶速度-应力方程能够很好地适应地下各向异性变密度介质，更加符合实际地质情况.因此，本文首先实现了基于一阶速度-应力方程的VTI介质LSRTM成像方法，通过理论公式推导，最终形成一套完整的基于VTI介质一阶速度-应力方程的LSRTM理论算法.最后通过模型测试验证了本文方法的正确性和优越性.
1 方法原理 1.1 一阶速度-应力方程VTI介质拟声波方程表征了地震波在各向异性介质中的传播规律，是后续正演模拟、偏移成像以及波形反演的研究基础.常用的二维VTI介质拟声波一阶速度-应力方程为(Duveneck et al., 2008)：
|
(1) |
其中，p，q分别表示应力的水平分量和垂直分量，u，w分别表示质点速度的水平分量和垂直分量，ρ为密度，ε，δ为表征各向异性特征的Thomsen参数，s为介质慢度，fx，fz分别是x，z方向的加载震源.
为了简化，可将公式(1)写成矩阵向量形式：
|
(2) |
其中，
|
表示正演模拟算子，
根据场的扰动理论，不失一般性，慢度场的平方s2可以表示为背景慢度场的平方s02与扰动慢度场的平方Δs2之和，即
|
(3) |
同样地，各向异性参数场也可以分解为背景场与扰动场的叠加(Dutta and Schuster, 2014)，即：
|
(4) |
其中，ε0，δ0表示背景场，Δε，Δδ表示扰动场.为了便于计算，文中参考Dutta和Schuster(2014)在多参数条件下提出的扰动简化方式，假设各向异性参数扰动较小，可以忽略不计，即，
|
(5) |
另外，地震波场满足叠加原理，在介质存在扰动时，地震波总波场U可以看作是背景介质产生的波场U0和扰动介质产生的波场Us的叠加，即
|
(6) |
背景波场与原始波场均满足原始波动方程(1)，将公式(3)、(4)、(5)、(6)代入公式(1)，并与背景波场满足的波动方程作差后应用Born近似，用背景波场U0代替地震波总波场U，即可推导出扰动波场Us的表达式(具体推导见附录A)：
|
(7) |
利用矩阵向量可以表示为：
|
(8) |
其中，扰动波场
从数学的角度来看，偏移算子是正演模拟算子的伴随(Plessix，2010)，故而本文利用伴随状态法求解偏移算子.由伴随状态法可知：
|
(9) |
这里，L*表示L的伴随矩阵，U*=(u* w* p* q*)表示伴随波场.
将L代入公式(9)进行计算求解得：
|
(10) |
公式(10)即为伴随算子.
基于反演思想的最小二乘逆时偏移(LSRTM)旨在获取地下最优的介质模型，使得正演模拟数据(即反偏移记录)与观测记录残差的模最小，是一个最小范数问题，定义目标泛函如下：
|
(11) |
其中，|| ||2表示L2-范数，〈，〉表示标量乘，Uobs=(uobs wobs pobs qobs)T表示地震波场的观测记录，即代表地下波场的真实传播过程.
为了便于求解公式(11)，这里引入目标泛函的梯度，即，
|
(12) |
对公式(8)左右两端的m同时求导得，
|
(13) |
这里，
将公式(13)代入公式(12)，并结合公式(9)得：
|
(14) |
其中，*表示伴随，Us－Uobs表示残差向量.这里定义伴随波场：
|
(15) |
进一步得到，
|
(16) |
将公式(10)代入公式(16)，整理可得伴随方程如下：
|
(17) |
观察公式(1)和(17)不难发现，基于VTI介质一阶速度-应力方程推导出的伴随方程，与原始波动方程之间有较大差异，而基于二阶标量声波方程得到的伴随方程和原始波动方程是相同的(李庆洋等，2016c)，这是一阶方程与二阶方程LSRTM算法的本质区别.
另外，结合公式(14)、(15)，梯度公式可简化为
|
(18) |
将上式展开即可得出梯度更新公式：
|
(19) |
成像结果m(k)可以通过最速下降法(Steepest Descent Method，SDM)或者共轭梯度法(Conjugate Gradient Method，CGM)进行迭代求解，由于CGM克服了SDM收敛速度慢的特点，是求解大型非线性问题最有效的算法之一，因此本文采用CGM进行迭代计算，具体过程如下：
|
(20) |
式中，α为计算步长；β为共轭梯度修正因子；z为共轭梯度；k为迭代次数；P为预条件算子，引入适当的约束条件可以使得反演问题更加稳定且加快收敛速度，本文采用震源波场的照明补偿作为预条件算子.
基于上述推导，可以归纳总结出VTI-LSRTM算法的具体实现流程，如图 1所示.
|
图 1 VTI-LSRTM算法实现流程 Fig. 1 Flow chart of LSRTM for VTI media |
为了验证本文方法的正确性及适应性，本文分别用气块模型、层状模型以及复杂Marmousi模型对算法进行测试.其中气块模型主要用来测试算法的正确性、偏移噪声压制能力及其在各向异性介质成像时的优越性；利用层状模型来定性分析密度的引入对成像分辨率及保幅性的影响；利用Marmousi模型来测试成像算法的稳定性及对复杂模型的适应性，重点考查非均匀变密度情况的适应性；最后，讨论了不同初始速度模型以及各向异性Thomsen参数存在误差时，成像算法的稳定性.
2.1 气块模型首先，本文以一个气块模型为例，验证算法的正确性.气块模型所对应的各参数场如图 2所示，其中，偏移速度场是由真实速度场经平滑后得到，进而求出慢度扰动.此外，密度场是由Gardner公式计算得出，计算公式为ρ=0.31ν0.25.模型参数如下：网格大小为401×201，纵、横向网格间距均为10 m；计算参数设置如下：观测系统为全速度场接收，采用主频为25 Hz的雷克子波作为震源，共激发201炮，炮间隔为20 m，每炮401道，道间隔为10 m，记录的时间长度为1.6 s，时间采样间隔为0.8 ms，采样点数为2001.观测记录是由线性正演算子经高阶有限差分方法模拟得到，边界条件采用改进的PML边界(黄建平等，2016b).
|
图 2 气块模型 (a)真实纵波速度；(b)偏移纵波速度；(c) ε模型；(d) δ模型. Fig. 2 Parameter fields of gas block model (a) True P-velocity model; (b) Migration velocity model; (c) ε model; (d) δ model. |
图 3(a，c)分别给出了各向同性声波RTM和VTI-RTM成像结果；图 3(b，d)分别为各向同性声波F-LSRTM和VTI-F-LSRTM成像结果.从RTM成像结果(图 3(a，c))可以看出，浅部地层含较强的偏移噪声和震源效应，信噪比较低，且深层倾斜同相轴能量相对较弱，虽然其基本能刻画地下层位信息，但振幅均衡性较差.而F-LSRTM成像结果(图 3(b，d))的成像质量有明显的改善，基本压制了浅层噪声和震源效应，且深部能量得到有效补偿，振幅更加均衡，成像分辨率明显提高.
|
图 3 成像结果 (a)各向同性声波RTM；(b)各向同性声波LSRTM迭代30次；(c) VTI-RTM；(d) VTI-LSRTM迭代30次. Fig. 3 Image result (a) Isotropic acoustic RTM; (b) Isotropic acoustic LSRTM with 30 iterations; (c) VTI-RTM; (d) VTI-LSRTM with 30 iterations. |
从各向同性声波偏移结果(图 3(a, b))可以看出，在梯形状各向异性地层中，包含有较强的伪横波干扰(如图中实线箭头所指处)，此外，倾斜界面成像位置有明显上倾的趋势(如图矩形框中所示)，相对真实位置有较大偏离，这是由于地层速度各向异性的存在，地震波沿其他方向传播的速度大于垂向速度，导致地震波走时相对减小所致，且由于倾斜界面明显偏移不足，倾斜界面拉长，进而在深部地层底部出现了明显的交叉(如图中虚线箭头所示)等偏移假象，容易对后续解释工作造成误导.而在各向异性偏移结果(图 3(c，d))中，消除了伪横波干扰、倾斜地层上倾、偏移假象等现象，成像位置更为准确，成像质量明显提高.
值得注意的是，各向同性声波RTM中出现的伪横波干扰、深部交叉现象等现象，并未随着F-LSRTM的迭代得到有效地压制和改善，反而变得更严重，这是由于地层的各向异性影响了地震波的传播规律，使得在各向同性假设下的模拟数据与实际观测记录不匹配，以致于在迭代过程中，真实有效构造信息难以准确充分识别，成像假象被当成有效信息处理，进而影响成像质量；而VTI-F-LSRTM考虑了介质各向异性的影响，使得模拟数据与实际观测数据匹配程度较好，随着迭代次数的增加，噪声得到压制，振幅更加均衡，成像质量更理想.
图 4为各向同性F-LSRTM(实线)和VTI-F-LSRTM(虚线)的残差收敛曲线.由图可知，反演成像过程较为稳定，随着迭代次数的增加，两种F-LSRTM算法的残差均快速减小，并逐渐收敛到一个相对稳定的值，而VTI-F-LSRTM收敛效果明显优于各向同性F-LSRTM，收敛到一个更低的阈值，这是由于VTI-F-LSRTM方法考虑了介质的各向异性，使得模拟数据与观测数据匹配较好.
|
图 4 数据残差归一化收敛曲线 Fig. 4 The normalized convergence curves of data residual |
为了进一步测试密度对成像结果的影响，本文设计了简单层状模型进行测试，模型对应的速度、密度参数设置如图 5a所示，密度设置基本满足Gardner公式，慢度扰动可看作理想的成像结果，如图 5b所示.详细参数如下，网格大小为150×150，纵、横向网格间距均为10 m；采用主频为25 Hz的雷克子波作为震源，一共激发75炮，炮间隔为20 m，每炮150道，道间隔为10 m，记录的时间长度为1.2 s，时间采样间隔为1.0 ms.炮记录是由线性正演算子经高阶有限差分方法模拟得到.
|
图 5 平层模型 (a)参数模型；(b)慢度扰动. Fig. 5 Layer model (a) Parameter model; (b) Slowness perturbation. |
为了对比分析密度对成像值的影响，消除其他因素(规则网格、交错网格等等)的干扰，这里采用一阶方程分别进行常密度(Constant Density, CD)CD-LSRTM以及变密度(Variable Density, VD)VD-LSRTM，成像结果如图 6所示.二者成像结果界面位置差异较小，表明密度对构造成像的影响较小.进一步考查成像结果对成像振幅的影响，如图 7为提取的偏移距700 m处的单道记录对比图，黑色实线为慢度扰动，红色实线为变密度结果，绿色虚线为常密度结果.从图 7的单道成像结果对比和局部放大图可以看出：变密度结果相对于常密度方程与真实慢度扰动有更好的一致性，在成像过程中考虑地下介质的密度变化，能更好地实现岩性成像.
|
图 6 成像结果对比(迭代10次) (a) VD-LSRTM；(b) CD-LSRTM. Fig. 6 Comparison of imaging results (a) VD-LSRTM; (b) CD-LSRTM. |
|
图 7 成像结果和真实慢度扰动单道对比 Fig. 7 Single trace contrast between imaging result and real slowness perturbation at the distance of 700 m |
进一步，考查新方法对复杂模型的适应性，Marmousi模型作为典型的标准复杂速度模型，含有强横向变速和不整合界面等复杂精细构造，常用于检测偏移成像方法对复杂构造的适应能力，可验证算法在复杂地质构造条件下的有效性和适应性.图 8a为Marmousi模型真实纵波速度，偏移速度场是由真实纵波速度经平滑后得到，如图 8b所示，借鉴Alkhalifah(1997)的做法，构建了Thomsen参数模型，如图 8(c, d)所示.密度场是由Gardner公式计算得出，计算公式为ρ=0.31ν0.25.各参数设置如下：模型网格大小为369×188，纵、横向网格间距均为12.5 m，观测系统为全速度场接收，采用主频为20 Hz的雷克子波作为震源，一共激发185炮，炮间隔为25 m，道间隔为12.5 m，时间采样间隔为0.8 ms，采样点数为5001.合成数据Ⅰ由Born正演算子经高阶有限差分方法模拟得到，合成数据Ⅱ由常规高阶有限差分正演模拟得到.需特别指出的是，记录数据考虑了密度变化的影响.
|
图 8 Marmousi模型 (a)真实纵波速度；(b)偏移速度；(c) ε模型；(d) δ模型. Fig. 8 Parameter fields of Marmousi model (a) True P-velocity model; (b) Migration velocity model; (c) ε model; (d) δ model. |
慢度扰动可视为理想的成像结果，如图 9a所示.图 9(b—f)是利用合成数据Ⅰ得到的成像结果(RTMⅠ)，图 9(g，h)是利用合成数据Ⅱ得到的成像结果(RTMⅡ).图 9b是各向同性F-LSRTMⅠ成像结果，分析可知，由于没有考虑介质各向异性的影响，致使同相轴紊乱甚至错位，绕射波不收敛，断层边界及深部构造难以分辨，且包含有大量噪声，成像质量较差；图 9c为VTI-F-RTMⅠ成像结果，剖面内包含大量的低频噪声及震源效应，地下真实构造被淹没，难以准确辨认；图 9(d, e)分别为VTI-F-LSRTMⅠ迭代15次和VTI-S-LSRTMⅠ迭代30次的成像结果，可以看出，二者的成像效果均明显优于图 9b，反射同相轴正确归位，成像质量均显著提高，与图 9c相比，低频噪声及震源效应得到了有效压制，地下真实构造得以准确识别，但进一步观察VTI-S-LSRTMⅠ结果不难发现，成像剖面内含有轻微噪声(如图中箭头所指处)，成像质量整体上略逊于VTI-F-LSRTMⅠ结果(图 9(d, f))，分辨率较低，究其原因可能是因为：(1)密度变化带来的影响；(2)S-LSRTMⅠ采用规则网格算法，而F-LSRTMⅠ采用交错网格算法，相比规则网格算法，交错网格频散压制效果更好、成像精度更高.图 9g为VTI-F-RTMⅡ成像结果，图 9h为VTI-F-LSRTMⅡ迭代15次的成像结果.与图 9g相比，图 9h成像振幅的均衡性有明显的增强，偏移噪声得到有效压制，陡倾断层刻画清晰，成像分辨率明显提高，成像质量显著改善.但仔细观察不难发现，VTI-F-LSRTMⅡ在经过多次迭代后，放大了一些成像噪声，这是采用最小二乘偏移方法所不可避免的问题(周华敏等，2014)，如何有效解决这一问题将是本文后续努力的方向.总的说来，本文所提VTI-F-LSRTM算法，能够有效压制偏移噪声，对小尺度构造刻画更好，能量补偿效果显著，成像分辨率明显提高.
|
图 9 慢度扰动及成像结果 (a)慢度扰动；(b)各向同性声波LSRTMⅠ迭代30次；(c) VTI-F-RTMⅠ；(d) VTI-F-LSRTMⅠ迭代15次；(e) VTI-S-LSRTMⅠ迭代30次；(f) VTI-F-LSRTMⅠ迭代30次；(g) VTI-F-RTMⅡ；(h) VTI-F-LSRTMⅡ迭代15次. Fig. 9 Image result (a) Slowness perturbation; (b) Isotropic acoustic LSRTMⅠ with 30 iterations; (c) VTI-F-RTMⅠ; (d) VTI-F-LSRTMⅠ with 15 iterations; (e) VTI-S-LSRTMⅠ with 30 iterations; (f) VTI-F-LSRTMⅠ with 30 iterations; (g) VTI-F-RTMⅡ; (h) VTI-F-LSRTMⅡ with 15 iterations. |
图 10给出了归一化后的数据残差收敛曲线(基于合成数据Ⅰ)，由图可知，反演成像过程较为稳定，随着迭代次数的增加，几种LSRTM算法的残差均收敛到一个相对稳定的值.而ISO-F-LSRTM(实线)忽略了各向异性的影响，致使模拟数据与观测数据差异较大，因此数据残差只收敛到68%左右，而前文提到的简单气块模型，采用ISO-F-LSRTM算法时，数据残差可以收敛至20%左右，由此可以看出，地下各向异性介质的复杂程度对反偏移数据有很大的影响，进而影响数据残差的收敛性.通过对比得到如下认识：(1)当研究对象为复杂各向异性介质时，尤其西部探区采集数据成像时，不能忽略各向异性的影响；(2)VTI-S-LSRTM (带圈实线)忽略了密度的影响，使模拟数据与观测数据之间也有一定的差异，但这种差异要比各向异性的影响小的多，此外，两种方程之间差分网格的差异也具有一定的影响，因而数据残差仅收敛至28%左右；(3)利用VTI-F-LSRTM(虚线)算法模拟出的数据与观测数据更为匹配，数据残差能够快速收敛到一个更低的水平，进而得到更好的结果，更能体现出本文算法的优势.
|
图 10 数据残差归一化收敛曲线 Fig. 10 Normalized convergence curve of data residuals |
为了测试本文方法对地下复杂各向异性介质中各个参数场(即速度、ε和δ)的敏感性.下面采取控制变量法，利用复杂Marmousi模型分别对速度、ε和δ进行敏感性测试.
2.4.1 速度敏感性测试在该项测试中，假设各向异性参数ε和δ正确，通过引入速度误差来分析速度对偏移成像的影响.偏移速度相对误差大小设计如下：3.16%、6.58%、9.45%.
成像结果如图 11所示，从图中可以看出，当速度误差为3.16%时，成像剖面中即包含大量的干扰噪声，同相轴成像位置不准确，振幅不均衡，成像剖面紊乱，成像质量较差，并且，随着速度相对误差的增大，成像位置偏差亦随之变大(如矩形框所示)，成像效果也变得更差，地下真实构造难以准确识别定位(如椭圆框所示).由此可见，速度误差对偏移结果有很大的影响.
|
图 11 含不同速度误差的成像结果(均迭代30次) (a)真实速度；(b)误差为3.16%；(c)误差为6.58%；(d)误差为9.45%. Fig. 11 Imaging results with different velocity errors (30 iterations) (a) True; (b) Error is 3.16%; (c) Error is 6.58%; (d) Error is 9.45%. |
带有不同速度误差的残差收敛曲线如图 12所示，经分析可知，速度误差严重影响反偏移数据，使其与观测数据不匹配，收敛效果很差，随着速度相对误差的增大，收敛效果随之变差，且当速度误差增长至9.45%时，误差曲线基本不收敛.
|
图 12 含不同速度误差的残差收敛曲线 Fig. 12 Residual convergence curve with different velocity errors |
在该项测试中，假设速度和各向异性参数δ正确，通过引入参数ε误差来分析ε对偏移成像的影响.参数ε相对误差大小设计如下：9.45%、19.75%、29.98%.
各成像结果如图 13所示，从图中可以看出，当有ε误差存在时，成像剖面内存在少量偏移噪声，但从整体上来看，成像质量依然很高，然而，随着ε相对误差的增大，偏移噪声逐渐增多，绕射波不收敛，同相轴成像位置相对不准确，虽然大尺度构造依然能够大致辨别，但地下复杂小尺度构造难以准确识别，整个成像剖面相对紊乱，成像效果逐渐变差(如椭圆框所示).
|
图 13 含不同ε误差的成像结果(均迭代30次) (a)真实ε；(b)误差为9.45%；(c)误差为19.75%；(d)误差为29.98%. Fig. 13 Imaging results with different ε errors (30 iterations) (a) True; (b) Error is 9.45%; (c) Error is 19.75%; (d) Error is 29.98%. |
残差收敛曲线如图 14所示，可以看出，与速度误差相比，参数ε误差对反偏移结果影响较小，残差收敛效果一般，但随着ε相对误差的增大，收敛效果亦随之变差.
|
图 14 含不同ε误差的残差收敛曲线 Fig. 14 Residual convergence curve with different ε errors |
在该项测试中，假设速度和各向异性参数ε正确，通过引入参数δ误差来分析δ对偏移成像的影响.参数δ相对误差大小设计如下：9.45%、28.35%、50.19%.
成像结果如图 15所示，从图中可以看出，当有δ误差存在时，基本不影响成像结果，成像剖面比较清晰，绕射波基本收敛，同相轴成像位置基本准确，并且当δ相对误差增大至50.19%时，成像剖面依然清晰，成像效果较好.
|
图 15 含不同δ误差的成像结果(均迭代30次) (a)真实δ；(b)误差为9.45%；(c)误差为28.35%；(d)误差为50.19%. Fig. 15 Imaging results with different δ errors (30 iterations) (a) True; (b) Error is 9.45%; (c) Error is 28.35%; (d) Error is 50.19%. |
残差收敛曲线如图 16所示，当δ相对误差增大至50.19%时，残差收敛曲线仍收敛至5%以下.经分析可知，参数δ误差对反偏移数据基本没有影响，反偏移数据与观测记录基本匹配，收敛效果较好，而且，随着δ相对误差的增大，收敛效果基本不变.
|
图 16 含不同δ误差的残差收敛曲线 Fig. 16 Residual convergence curve with different δ errors |
总的来说，速度误差对成像结果的影响最大，其次是参数ε，而参数δ对成像结果的影响最小，几乎可以忽略.所以，在实际生产资料处理过程中，速度是偏移成像的核心，一点微小的误差即会造成地层成像位置不准确，成像质量变差.此外，介质的各向异性也不容无视，对于复杂Marmousi模型来说，当参数ε相对误差为29.98%时，成像结果中出现较多偏移噪声，成像剖面相对紊乱，绕射波不收敛，精细构造难以识别，且数据残差收敛效果较差.
3 结论与认识本文综合考虑地下介质的各向异性以及介质的非均匀性，在各向同性LSRTM的基础上，发展了基于一阶速度-应力方程的非均匀变密度VTI介质LSRTM算法.通过数值模拟分析可以得出如下几点认识：
(1) 与各向同性LSRTM算法相比，VTI-LSRTM考虑了介质的各向异性，消除了介质各向异性引起的成像误差，使模拟数据与观测数据更加匹配，进而得到更好的成像结果.
(2) 现有的LSRTM理论算法大多是以二阶常密度标量声波方程为基础，忽略了密度变化的影响，无法有效处理地下变密度介质.而本文LSRTM算法是基于一阶方程展开的，便于使用交错网格，比规则网格算法拥有更高的精度，压制频散效果更好，且本文算法易于处理地下非均匀变密度介质.随着迭代次数的递增，数据残差曲线(图 10)能够快速稳定收敛，成像质量亦随之提高，体现了本文算法一定的有效性.
(3) 通过对速度、各向异性参数误差的敏感性分析可知，速度误差对成像结果的影响最大，其次是参数ε，而参数δ对成像结果的影响最小，这对于后续多参数的波形反演、各向异性参数提取等具有一定的借鉴作用.
本文基于一阶速度-应力方程的VTI介质最小二乘偏移方法，能够较好的实现反射波归位和绕射波收敛，且能考虑地下介质的密度变化，但仍然存在一些需要完善和改进的地方：(1)从参数敏感性分析测试可以看出，LSRTM算法对速度和各向异性参数的依赖性比较强，这必将制约其实际发展，可通过研究反演策略以提高速度和各向异性参数反演的准确性，同时修改目标泛函以减弱方法对参数的敏感性；(2)考虑到LSRTM算法计算量较大，可引入GPU加速、编码策略等进一步提高计算效率；3)由于预条件算子和正则化算子等能够提高算法稳定性、加快收敛速度，可进行针对性研究以进一步优化算法.此外，本文仅针对VTI介质拟声波方程进行研究，需进一步研究更为复杂的弹性波，亦或是将其推广至TTI介质等等以增强其实用性.
附录A在该附录中，给出了本文扰动波场(7)式的推导过程：
首先，将公式(3)、(4)、(5)、(6)代入公式(1)得，
|
(A1) |
另外，背景波场满足如下波动方程：
|
(A2) |
将(A2)代入(A1)并结合假设条件公式(5)得，
|
(A3) |
然后应用Born近似条件，即用背景波场p0、q0分别代替总波场p0+ps、q0+qs，并进行简单移项整理得，
|
(A4) |
这里，利用模型参数m代替(－Δs2).公式(A4)表明，扰动波场与慢度扰动呈线性关系.通过观察分析可知，扰动波场与背景波场具有相同的波动方程，仅在震源项有所差别，由此看出，扰动波场是由背景波场关于时间的一阶导数与慢度扰动的乘积作为震源激发产生的.
Alkhalifah T. 1997. An anisotropic Marmousi model. Stanford Exploration Project, SEP-95: 265-282.
|
Alkhalifah T. 1998. Acoustic approximations for processing in transversely isotropic media. Geophysics, 63(2): 623-631. DOI:10.1190/1.1444361 |
Alkhalifah T. 2000. An acoustic wave equation for anisotropic media. Geophysics, 65(4): 1239-1250. DOI:10.1190/1.1444815 |
Claerbout J F. 1992. Earth Soundings Analysis-Processing Versus Inversion. Boston: Blackwell Scientific Publications.
|
Clapp M L, Clapp R G, Biondi B L. 2005. Regularized least-squares inversion for 3-D subsalt imaging.//75th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1814-1817.
|
Dai W, Wang X, Schuster G T. 2011. Least-squares migration of multisource data with a deblurring filter. Geophysics, 76(5). |
Dai W, Fowler P, Schuster G T. 2012. Multi-source least-squares reverse time migration. Geophysical Prospecting, 60(4): 681-695. DOI:10.1111/gpr.2012.60.issue-4 |
Du X, Bancroft J C, Lines L R. 2005. Reverse-time migration for tilted TI media.//75th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1930-1934.
|
Dutta G, Schuster G T. 2014. Attenuation compensation for least-squares reverse time migration using the viscoacoustic-wave equation. Geophysics, 79(6): S251-S262. DOI:10.1190/geo2013-0414.1 |
Duveneck E P, Milcik P, Bakker P M, et al. 2008. Acoustic VTI wave equations and their application for anisotropic reverse-time migration: 78th Annual International Meeting, SEG, Expanded Abstracts, 2186-2189.
|
Fletcher R P, Du X, Fowler P J. 2009. Reverse time migration in tilted transversely isotropic (TTI) media. Geophysics, 74(6): WCA179-WCA187. DOI:10.1190/1.3269902 |
Huang J P, Li Z C, Kong X, et al. 2013. A study of least-squares migration imaging method for fractured-type carbonate reservoir. Chinese Journal of Geophysics (in Chinese), 56(5): 1716-1725. DOI:10.6038/cjg20130529 |
Huang J P, Cao X L, Li Z C, et al. 2014. Least square reverse time migration in high resolution imaging of near surface. Oil Geophysical Prospecting (in Chinese), 49(1): 107-112. |
Huang J P, Li C, Li Q Y, et al. 2015. Least-squares reverse time migration with static plane-wave encoding. Chinese Journal of Geophysics (in Chinese), 58(6): 2046-2056. DOI:10.6038/cjg20150619 |
Huang J P, Li C, Li Q Y. 2016a. Least-squares Migration Imaging Theory and Method. Beijing: Science Press.
|
Huang J P, Yang Y, Li Z C, et al. 2016b. Lebedev grid finite-difference modeling for irregular free-surface and stability analysis based on M-PML boundary condition. Journal of China University of Petroleum (in Chinese), 40(4): 47-56. |
Küehl H, Sscchi M D. 2003. Least-squares wave-equation migration for AVP/AVA inversion. Geophysics, 68(1): 262-273. DOI:10.1190/1.1543212 |
Klié H, Toro W. 2001. A new acoustic wave equation for modeling in anisotropic media.//71st Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1171-1174.
|
Li Q Y, Huang J P, Li Z C, et al. 2016a. Optimized multi-source least-squares reverse time migration. Oil Geophysical Prospecting (in Chinese), 51(2): 334-341. |
Li Q Y, Huang J P, Li Z C, et al. 2016b. Mean-residual normalized cross-correlation least-squares reverse time migration and its application. Chinese Journal of Geophysics (in Chinese), 59(8): 3006-3015. DOI:10.6038/cjg20160823 |
Li Q Y, Huang J P, Li Z C, et al. 2016c. Multi-source least-squares reverse time migration based on first-order velocity-stress wave equation. Chinese Journal of Geophysics (in Chinese), 59(12): 4666-4676. DOI:10.6038/cjg20161226 |
Li Z C, Guo Z B, Tian K. 2014. Least-squares reverse time migration in visco-acoustic medium. Chinese Journal of Geophysics (in Chinese), 57(1): 214-228. DOI:10.6038/cjg20140118 |
Li Z C, Huang J Q, Huang J P, et al. 2017. Fast least-squares reverse time migration based on plane-wave encoding for VTI media. Chinese Journal of Geophysics (in Chinese), 60(1): 240-257. DOI:10.6038/cjg20170120 |
Liu Y J, Li Z C, Wu D, et al. 2013. The research on local slope constrained least-squares migration. Chinese Journal of Geophysics (in Chinese), 56(3): 1003-1011. DOI:10.6038/cjg20130328 |
Nemeth T, Wu C J, Schuster G T. 1999. Least-squares migration of incomplete reflection data. Geophysics, 64(1): 208-221. DOI:10.1190/1.1444517 |
Pestana R C, Ursin B, Stoffa P L. 2011. Separate P-and SV-wave equations for VTI media.//81st Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 163-167.
|
Plessix R E. 2010. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications. Geophysical Journal of the Royal Astronomical Society, 167(2): 495-503. |
Tarantola A. 1984. Linearized inversion of seismic reflection data. Geophysical Prospecting, 32(6): 998-1015. DOI:10.1111/gpr.1984.32.issue-6 |
Tatham R H, McCormack M D. 1991. Multicomponent Seismology in Petroleum Exploration. Society of Exploration Geophysicist. |
Tsvankin I, Thomsen L. 1994. Nonhyperbolic reflection moveout in anisotropic media. Geophysics, 59(8): 1290-1304. DOI:10.1190/1.1443686 |
Wang Y F. 2013. Comparison of interferometric migration and preconditioned regularizing least squares migration inversion methods in seismic imaging. Chinese Journal of Geophysics (in Chinese), 56(1): 230-238. DOI:10.6038/cjg20130123 |
Wu G C. 2006. Seismic Wave Propagation and Imaging in Anisotropic Media (in Chinese). Dongying: China University of Petroleum Press.
|
Xu S, Zhou H B. 2014a. Pure quasi-P wave calculation in anisotropic media.//84th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 3887-3891.
|
Xu S, Zhou H B. 2014b. Accurate simulations of pure quasi-P-waves in complex anisotropic media. Geophysics, 79(6): T341-T348. DOI:10.1190/geo2014-0242.1 |
Zhan G, Pestana R C, Stoffa P L. 2011. An acoustic wave equation for pure P wave in 2D TTI media.//81st Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 168-173.
|
Zhang L B, Rector Ⅲ J W, Hoversten G M. 2005. Finite-difference modelling of wave propagation in acoustic tilted TI media. Geophysical Prospecting, 53(6): 843-852. DOI:10.1111/gpr.2005.53.issue-6 |
Zhang Y, Zhang H Z. 2009. A stable TTI reverse time migration and its implementation.//79th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
|
Zhang Y, Zhang H, Zhang G Q. 2011. A stable TTI reverse time migration and its implementation. Geophysics, 76(3): WA3-WA11. DOI:10.1190/1.3554411 |
Zhou H, Zhang G, Bloor R. 2006a. An anisotropic acoustic wave equation for VTI media.//68th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
|
Zhou H B, Bloor R, Zhang G Q. 2006b. An anisotropic acoustic wave equation for modeling and migration in 2D TTI media.//76th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 194-198.
|
Zhou H M, Chen S C, Ren H R, et al. 2014. One-way wave equation least-squares migration based on illumination compensation. Chinese Journal of Geophysics (in Chinese), 44(1): 369-374. DOI:10.6038/cjg20140823 |
Zhu F, Huang J P, Li Z C, et al. 2017. PLSFFD pre-stack depth migration for VTI medium.//87th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
|
黄建平, 李振春, 孔雪, 等. 2013. 碳酸盐岩裂缝型储层最小二乘偏移成像方法研究. 地球物理学报, 56(5): 1716-1725. DOI:10.6038/cjg20130529 |
黄建平, 曹晓莉, 李振春, 等. 2014. 最小二乘逆时偏移在近地表高精度成像中的应用. 石油地球物理勘探, 49(1): 107-112. |
黄建平, 李闯, 李庆洋, 等. 2015. 一种基于平面波静态编码的最小二乘逆时偏移方法. 地球物理学报, 58(6): 2046-2056. DOI:10.6038/cjg20150619 |
黄建平, 李闯, 李庆洋. 2016a. 最小二乘偏移成像理论及方法. 北京: 科学出版社.
|
黄建平, 杨宇, 李振春, 等. 2016b. 基于M-PML边界的Lebedev网格起伏地表正演模拟方法及稳定性分析. 中国石油大学学报:自然科学版, 40(4): 47-56. |
李庆洋, 黄建平, 李振春, 等. 2016a. 优化的多震源最小二乘逆时偏移. 石油地球物理勘探, 51(2): 334-341. |
李庆洋, 黄建平, 李振春, 等. 2016b. 去均值归一化互相关最小二乘逆时偏移及其应用. 地球物理学报, 59(8): 3006-3015. DOI:10.6038/cjg20160823 |
李庆洋, 黄建平, 李振春, 等. 2016c. 基于一阶速度-应力方程的多震源最小二乘逆时偏移. 地球物理学报, 59(12): 4666-4676. DOI:10.6038/cjg20161226 |
李振春, 郭振波, 田坤. 2014. 黏声介质最小平方逆时偏移. 地球物理学报, 57(1): 214-228. DOI:10.6038/cjg20140118 |
李振春, 黄金强, 黄建平, 等. 2017. 基于平面波加速的VTI介质最小二乘逆时偏移. 地球物理学报, 60(1): 240-257. DOI:10.6038/cjg20170120 |
刘玉金, 李振春, 吴丹, 等. 2013. 局部倾角约束最小二乘偏移方法研究. 地球物理学报, 56(3): 1003-1011. DOI:10.6038/cjg20130328 |
王彦飞. 2013. 地震波干涉偏移及预条件正则化最小二乘偏移成像方法对比. 地球物理学报, 56(1): 230-238. DOI:10.6038/cjg20130123 |
吴国忱. 2006. 各向异性介质地震波传播与成像. 东营: 中国石油大学出版社.
|
周华敏, 陈生昌, 任浩然, 等. 2014. 基于照明补偿的单程波最小二乘偏移. 地球物理学报, 57(8): 2644-2655. DOI:10.6038/cjg20140823 |