各向异性介质Low-rank有限差分法纯qP波叠前平面波最小二乘逆时偏移

各向异性介质Low-rank有限差分法纯qP波叠前平面波最
小二乘逆时偏移
黄金强;李振春
【摘要】拟声波最小二乘逆时偏移是一种极具潜力的地震波成像工具,但该方法遭受各向异性拟声波近似的限制,TTI介质正演模拟不稳定、反偏移记录中遭受伪横波二次扰动及数值频散假象,另外拟声波最小二乘逆时偏移还面临计算效率低、收敛速度慢、对速度等模型参数依赖性高等问题.为了克服各向异性拟声波最小二乘逆时偏移的缺陷,在反演框架下,本文借助Low-rank有限差分算法首次提出并实现了TTI介质纯qP波线性正演模拟及纯qP波最小二乘逆时偏移;为了进一步提升反演成像效率,同时改善反演成像方法对模型参数误差的依赖性及对地震数据噪声的适应性,通过引入叠前平面波优化策略,发展了TTI介质纯qP波叠前平面波最小二乘逆时偏移成像方法.在编程实现方法的基础上,通过开展模型成像测试,展示了本方法的优势和潜力:一方面加快了反演成像效率,另一方面也提升了方法的抗噪性,同时还降低了方法对模型参数的依赖性.
【期刊名称】《地球物理学报》
【年(卷),期】2019(062)008
【总页数】24页(P3106-3129)
【关键词】各向异性介质;Low-rank有限差分;纯qP波;最小二乘逆时偏移;叠前平面波最小二乘逆时偏移
【作者】黄金强;李振春
【作者单位】贵州大学资源与环境工程学院,贵阳550025;中国石油大学(华东)地球科学与技术学院,山东青岛,266580
【正文语种】中文
【中图分类】P631
0 引言
最小二乘反演成像方法早已成为油气勘探领域的一种高分辨率探测工具，但由于基于各向同性假设的声波最小二乘偏移方法忽略了介质的各向异性，可能导致强各向异性地层及其下覆构造出现成像误差拉大、分辨率降低、偏移噪声严重、甚至迭代发散等问题，严重损坏反演结果的质量(李振春等，2017).此外，在“声学近似”理论下，各向异性介质拟声波方程偏移及反偏移算子始终存在纵横波耦合及不稳定现象，因此，常规的各向异性拟声波最小二乘偏移也面临诸多缺陷.鉴于此，要开展各向异性储层油气勘探与开发，关键在于实现针对纯qP波信息的快速且精确成像.
Tarantola(1984)建立了比较完备的最小二乘反演基础理论，其基本思想是利用模拟数据与观测数据的最佳匹配程度来设计目标泛函，进而达到更新并改善成像结果直至最佳成像的目的，这种从偏移成像过渡到反演成像的理论突破引起了国内外专家学者的深入探讨.由于当时计算机技术的限制，早期的最小二乘偏移实现主要依托于Kirchhoff积分算子(Nemeth et al.，1999；Duquet et al.，2000)，因其计算成本低，被广泛应用于商业软件中；随后，伴随着相移法等单程波算子的不断运用，出现了实现方式灵活多样的最小二乘单程波偏移(Kuehl and Sacchi，2002；Clapp et al.，2005；Kaplan et al.，2010；Huang and Schuster，2012)；近年来，Dai和Schuster(2009)率先提出并实现了最小二乘逆时偏移及多震源优化
策略，在双程波方程的约束下实现反演成像，理论严谨、成像优势显著，迅速成为工业界与学术界探讨的前沿热点(Tang，2009；Dai et al.，2011，2012；Wong et al.，2010，2011；Dong et al.，2012).为了适用于不同类型的地球介质，国
内外学者先后提出了弹性波最小二乘逆时偏移(Gu et al.，2017；Feng and Schuster，2017；Stanton and Sacchi，2017；Duan et al.，2017；Chen
and Sacchi，2017)、黏介质最小二乘逆时偏移(李振春等，2014；Dutta and Schuster，2014；Guitton et al.，2007；Sun et al.，2016；Guo and McMechan，2018)与各向异性介质最小二乘逆时偏移(李振春等，2017)等反演
成像方法；紧接着，针对反演成像中的不足，又相继发展了不依赖于地震子波(李
庆洋等，2017)、缓解速度模型依赖性(Dai and Schuster，2013)、利用多次波信息进行反演成像(刘伊克等，2018；刘学建和刘伊克，2016；Li et al.，2017；Wong et al.，2015)等的一系列优化措施，用以弥补最小二乘逆时偏移的不足，
并进一步改善反演成像质量；在迭代算法方面，先后实现了最速下降法、共轭梯度法、抛物拟合法及L-BFGS方法等迭代优化方案(Nocedal，1980；Hou and Symes，2015；Vigh and Starr，2008)；在目标泛函构建上，通过修改成本函数设计了基于l2模、基于l1模、基于混合l1/l2模与基于互相关等方式的目标泛函
构建方式(Gu et al.， 2017；李庆洋等，2016；敖瑞德等，2015).但调研发现，
各向异性介质最小二乘逆时偏移成像方法研究较少，且依然存在如下两方面的不足：首先，拟声波最小二乘逆时偏移仍然存在计算效率低、收敛速度慢、对速度等模型参数依赖性高等问题.各向异性介质最小二乘偏移只在速度模型与各向异性参数都
较为准确的情况下，才能保证反演结果有效地收敛，若模型参数存在较大的系统误差，则很可能出现迭代收敛慢甚至发散(李振春等，2017).从已出版的文献中可见，针对上述问题已有大量的解决方案，Dai和Schuster(2013)提出采用平面波最小
二乘逆时偏移来提升计算速率，并发展了叠前平面波最小二乘逆时偏移进一步改善
方法对速度误差的依赖性；李振春等(2017)将平面波编码策略应用到VTI介质最
小二乘逆时偏移成像中，实现了基于平面波加速的VTI介质最小二乘逆时偏移，
在一定程度上加快了反演成像效率，但该方法对模型参数误差的敏感性依然较高；Liu等(2013)通过引入空移成像条件，提出了扩展最小二乘逆时偏移，有效地降低了最小二乘逆时偏移对速度模型的依赖性；在此基础上，Hou和Symes(2016)进一步采用加权共轭梯度法来加快收敛速度，缓解了收敛速度慢的问题.为了改善各
向异性介质最小二乘逆时偏移对各向异性参数的依赖性，本文欲将叠前平面波优化策略应用于各向异性介质最小二乘逆时偏移中.
其次，常规的拟声波方程反偏移记录存在伪横波二次扰动，并伴有低速伪横波引入的数值频散假象，当该方法拓展到TTI介质时，在对称轴倾角急剧变化的强各向异性区域极易出现模拟不稳定现象，导致地震波异常振幅掩盖真实的有效波信息，这是各向异性声学近似理论的固有缺陷，在最小二乘反演成像中表现更为剧烈(黄金
强等，2017).从传统的拟声波方程出发，根据扰动理论与VTI介质Born近似假设，所推导的线性化正演算子也即反偏移算子遭受伪横波影响，在各向异性区域将出现伪横波二次扰动，不利于迭代的进行.为了准确模拟复杂介质中纯qP波的传播特征，Zhan等(2012，2013)提出在声学近似的基础上进一步简化相速度公式，进而导
出完全不含伪横波干扰的TTI介质纯qP波控制方程，但是该类方程在时空域表现为拟微分形式，求解困难，且难以运用Born近似求取纯qP波线性化正演算子；
伪横波在声学近似下之所以表现为菱形干扰，是因为令沿对称轴方向的横波速度为零时，横波群速度曲线正好为菱形形态，伪横波产生的本质原因在于对qP波相速度公式中的根式项进行了平方处理，进而引入了增根(Chu et al.，2013；Sheng and Zhou，2014；黄建平等，2016；程玖兵等，2013).基于此，Song等(2016)提出采用函数拟合的方法来近似相速度的根式项，进而消除伪横波干扰，从而实现纯qP波数值模拟，但该方法在某些相角区间内误差较大，且难以扩展到TTI介质
中.
为了克服各向异性声学近似理论的固有缺陷，黄金强和李振春(2017)将Fomel等(2010)提出的Low-rank分解策略运用于各向异性介质正演模拟中，由于Low-rank分解无需对频散关系进行近似，因而保证了地震波运动学特征的准确性且消除了伪横波假象及由此产生的一系列问题，随后将该方法成功应用于逆时偏移成像中，成像效果显著，但计算成本较高；为了继承有限差分的效率优势，Song等(2011，2013)将Low-rank分解与有限差分结合，提出了效率更高的Low-rank 有限差分法，并推广到TTI介质高精度正演模拟中，准确且高效地模拟出了TTI介质纯qP波信息，编程方式简单灵活，有很好的推广价值.对此，本文选用Low-rank有限差分法来设计各向异性介质纯qP波线性化正演算子，进而避免传统拟声波最小二乘逆时偏移的缺陷.
为了消除伪横波干扰等对反演成像的影响，本文在前人研究基础上，将Low-rank 有限差分算法应用到各向异性TTI介质正、反向波场延拓中，首次推导出TTI介质纯qP波线性化正演算子，进而在反演框架下实现TTI介质纯qP波最小二乘逆时偏移；为了提升反演成像的计算效率，同时改善本方法对模型参数误差的依赖性，及对地震数据噪声的适应性，通过引入叠前平面波优化策略，最终形成一种优化的各向异性纯qP波叠前平面波最小二乘逆时偏移成像方法；在编程实现方法的基础上，通过开展模型成像测试，展示了本方法的优势和潜力：一方面加快了反演成像的效率，另一方面提升了方法的抗噪性，同时还降低了方法对模型参数的依赖性.
1 方法原理
1.1 Low-rank有限差分法波场延拓
准确高效的纯qP波正、反向延拓算子是实现各向异性介质正演模拟、逆时偏移、最小二乘逆时偏移以及多参数全波形反演的第一步.近年来，基于Low-rank有限差分法的波场延拓算法因其计算效率高、无伪横波干扰等优势，正被广泛应用于各
向异性介质地震波正演模拟、逆时偏移处理实践中(Song and Alkhalifah，2013；黄金强和李振春，2017；黄金强等，2018).根据前人研究，记在空间位置为x、t
时刻的地震波场为p(x,t)，在相同空间位置、上一时刻(t-Δt)的地震波场为p(x,t-
Δt)，则在下一时刻的地震波场p(x,t+Δt)满足如下关系：
p(x,t+Δt)
+p(xR,t)]+p(x,t)+δ(x-xs)f(t),
(1)
其中，f(t)代表震源子波函数，xs为震源坐标位置，δ(x)为狄拉克函数，矩阵G代表空间x处的Low-rank有限差分权系数，与模型参数有关，并构成一一对应关系，具体的计算步骤可参见(黄金强等，2018)；为三维整数矢量，代表Low-rank有限差分加权所需的第l个网格点对应位置坐标，L代表差分模板长度(即所用的邻近网格点数)，决定差分阶数.
分析公式(1)可见：①该差分系数与Taylor差分系数不同，从频散关系来看，方刚等(2014)证实它是一种全局优化的差分系数，而Taylor差分系数仅对低波数成分
有效，只是局部优化的差分系数，因此该差分系数在频散压制方面具有一定的优势；
②从公式特点来看，通过对邻近两个时刻的地震波场进行加权求和即可实现地震波场的前向延拓，无需计算复杂的空间偏导数算子及混合偏导数算子，易于计算机编程；③同时，该种延拓方式具有Low-rank分解无伪横波干扰等优势：能够准确计算出纯qP波信息、适应较大时间空间步长、适应高频震源，还兼具有限差分法快速高效的优点，计算效率较高，节省成本；④对于不同类型的地球介质拥有完全相同的延拓模式和计算量，因为没有波动方程的概念，所以不会因为介质复杂度增加而增加计算量，方法及程序都易于直接拓展到其他介质中，唯一不同之处在于差分系数G的值不尽相同，需预先计算并存储.
如果说函数是数集到数集的映射，泛函是空间到数集的映射，那么算子可以理解为
是空间到空间的映射，即正演算子可定义为模型空间到数据空间的映射，当然实现该类映射的地球物理数值计算方法有很多，如双程波类有限差分正演、伪谱法正演、谱元法正演、单程波类正演、射线类正演等，这里的Low-rank有限差分法也是一种双程波类正演算子，但是没有显式的波动方程，而是一个离散的半解析的递推公式.因此，类似于矩阵表示方法，公式(1)也可以表示为矩阵算子的形式：
d=Am,
(2)
其中，m表示模型空间，这里代表网格化的速度及各向异性参数；A表示由Low-rank有限差分法所设计的双程波正演算子，它由模型参数决定；d表示数据空间，也即通过数值模拟得到的地震数据.
1.2 Low-rank有限差分法线性正演与伴随
Low-rank有限差分是消除各向异性介质中存在伪横波干扰，实现纯qP波场延拓
的有效工具.若已知沿着对称轴方向的qP波背景相慢度sp0及背景各向异性参数，并选定差分系数模板，可以确定Low-rank有限差分系数矩阵G0，同时给定震源坐标xs及震源函数f(t)，可依据公式(1)构建背景波场或入射波场p0(x,t+Δt)：
+p 0(x R,t)]+p 0(x,t)+δ(x-x s)f(t).
(3)
反偏移是利用成像结果或者说近似反射系数通过两次正演模拟来获取线性正演记录，这也是最小二乘逆时偏移成像的核心步骤.由于基于拟声波方程的反偏移算子存在
伪横波噪声，造成两次正演所得到的反偏移记录面临伪横波二次扰动、严重的数值频散和不稳定现象.为此，下文将借助Low-rank有限差分工具，推导TTI介质纯qP波线性化正演模拟算子，即纯qP波反偏移算子，这是本文的创新所在.定义慢
度扰动这里sp表示地下真实相慢度，Δsp表示扰动慢度.根据地震波摄动理论与
Born近似理论，背景慢度sp0与背景波场p0存在线性关系，并采用背景波场p0代替总场p，可导出基于Low-rank有限差分法的线性化地震散射波场延拓公式：
(4)
这里，ps为扰动波场，也即地震散射波场.公式(4)清楚地反映了反射系数模型与扰动波场之间的定量关系.实际上，计算散射波场仍需进行两步正演模拟，第一步是
基于公式(3)通过背景模型参数与震源函数来计算背景波场p0，当然，Low-rank
有限差分系数矩阵根据背景模型参数及选定的差分模板确定；第二步是基于公式(4)采用相同的差分系数矩阵，通过加载入射波场的一阶偏导数与反射系数的乘积作为二次震源激发来获得散射波场ps.为了便于表示，上述两步正演可用矩阵算子表示为：
dcal=Lm,
(5)
其中，m表示反射系数模型，L表示线性化正演算子，dcal表示散射正演计算所
得的地震数据.
同理，假设在炮点xs处激发，在检波点坐标xg处接收的观测炮记录为d(xg,t;xs)，则通过伴随状态法可以重构伴随波场p*，基于Low-rank有限差分法的伴随方程
可表示为如下形式：
(6)
其中，p*(x,t)表示t时刻的检波点波场或伴随波场，也即反传波场.观察可见：在时间上，重构反传波场是逆时传播，而计算震源波场是正向传播，重构检波点波场
p*与计算震源波场p0所采用的差分系数矩阵G0完全相同.
通过应用震源波场的二阶偏导数与检波点波场的零延迟互相关，可得炮点坐标为
xs所对应的单炮成像结果：
(7)
式中，m(x;xs)表示炮点坐标为xs所对应的单炮成像结果.为了简化上述成像公式，用矩阵算子符号代替公式(6)和公式(7)，进而将其改写为：
mmig=LTdcal,
(8)
其中，与线性正演算子L类似，LT表示偏移成像算子，mmig表示单炮成像结果.值得一提的是，本文给出的线性正演算子、偏移成像算子既源于波动方程线性化之后的结果，但又不同于波动方程的线性化形式，因为Low-rank有限差分法地震波场延拓公式并不是显式的解析解形式，只是伪解析解，所以只能通过数值试验来验证其共轭特征.幸运的是，有限差分与Low-rank有限差分等数值计算方法并不影
响算子的性质特征，因此上述两算子完全互为共轭，这也是进行最小二乘逆时偏移的基本要求.
1.3 Lrfd-LSRTM与Lrfd-Pre-PLSRTM
对于经典的炮域最小二乘逆时偏移，在理论上，反射系数模型m与观测系统类型
无关，对于Ns炮地震数据集，其目标泛函可表示为：
这里，m表示待求反射系数模型.每次迭代的最终成像结果是通过叠加所有的独立单炮成像结果所得.从公式可以看出，在编程实现迭代更新的过程中，①第一步是对观测数据做一次双程波RTM获得初始反射系数模型；随后，②输入初始模型与各向异性参数模型，利用公式(5)计算出预测炮数据，并与观测数据做差获得数据残差；③根据公式(8)，对炮数据残差做偏移进而计算出当前炮的反射系数模型修正量；④在迭代中，对每一炮数据是逐炮分别进行计算，最后叠加所有的修正量获得总的模型更新量，最终更新模型得到叠加后的成像结果，循环往复直至残差数据与预定阈值相吻合为止.从步骤来看，这与标准的声波最小二乘逆时偏移(LSRTM)基本一致，唯一的不同点在于采纳的线性正演算子为基于Low-rank有限差分的适用于各向异性介质的线性正演算子，偏移成像算子也是基于Low-rank有限差分的适用于各向异性介质的偏移成像算子，因此本文将该方法简称为Lrfd-LSRTM.
在实现Lrfd-LSRTM成像中，需慎重选择如下关键参数：①首先是差分模板L，L 越大，波场延拓精度越高，但计算效率变低，因此在选择中应权衡精度和效率；②其次是Low-rank分解中秩的选择，本文采用预设误差与最大秩相结合的办法，即秩的选择优先满足Low-rank分解中的误差值小于预设误差，但秩的容许值不能超过给定的最大秩；③再就是迭代次数，在理论上，稳健的最小二乘偏移方法随着迭代次数的增加，成像精度逐步升高，但对成像质量的改善力度将明显减缓，本文测试中设置最大迭代次数为30次.炮域LSRTM虽然缓解了诸如低频噪声强、深部能量弱、振幅不均衡等系列问题，但是当速度及各向异性参数模型的误差较大时，来自不同炮的单炮成像结果在相同成像位置处并不一定完全相同，这就很可能导致叠加后的成像结果出现成像模糊、甚至收敛变慢或发散，即炮域LSRTM对速度误差敏感性较高.此外，炮域LSRTM对数据噪声的适应能力较差，且计算成本高、内存消耗大.
为了进一步提升反演成像效率，同时改善反演成像方法对模型参数误差的依赖性及对地震数据噪声的适应性，本文将叠前平面波编码策略引入到各向异性介质Low-rank有限差分法最小二乘逆时偏移成像中.其中，平面波编码提升了反演成像效率、减小了计算机内存消耗，而叠前平面波编程实现方式进一步改善了算法的稳健性，在迭代中自动优化方法的抗噪性并降低方法对模型误差的敏感性.所谓叠前平面波
编码，主要是优化了平面波最小二乘逆时偏移的实现方式，其第一步仍然是对炮域地震数据按照线性时移的方式合成平面波域地震数据，在平面波域或射线参数p 域，合成地震记录可表示为
(10)
d(xg,t;xs)表示炮域的地震数据；d(xg,t;p)表示平面波域地震数据，即平面波记录；δ(t-p·xs)是时移函数，*表示卷积运算，xs是每一炮震源坐标，xg是接收点坐标，时移量为p·xs，它随震源坐标线性变化.因此，平面波编码过程可以表述为：首先
将炮域的地震数据进行逐炮时移，然后进行线性叠加，最终将多炮地震数据压缩成1个平面波记录.当然选择不同的射线参数，其时移量就不同，也就合成了不同的
平面波道集.
第二步是对合成的Np(≪Ns)个平面波道集进行逐个偏移成像，获得Np个独立的
成像结果.具体方式为采用相同的编码方式合成对应的平面波震源，随后激发产生
平面波震源波场，将合成的平面波道集作为边界条件反传获得平面波反传波场，最后将两者进行零延迟互相关即可获得不同平面波道集对应的成像结果
mi(i=1,…,Np).
与平面波偏移不同的是，需对上述Np个平面波成像结果进行分别存储，而不是叠加成一个剖面.Low-rank有限差分法叠前平面波LSRTM的处理流程与Lrfd-LSRTM类似，只是在每次迭代中，观测记录为平面波记录，需按照公式(10)对激
发震源进行平面波编码处理.此外，每个平面波道集所对应的成像结果都不进行叠
加成像处理，而是对其进行分别存储，且在迭代中独自地进行更新，相互之间互不影响，即在进行下一次迭代时，输入的反射系数模型为叠加之前对应的前一次成像结果.这种做法的优势在于，有利于减弱LSRTM对偏移速度与各向异性参数模型
误差的依赖性、抑制数据噪声、且方便提取平面波域共成像点道集用于质量监控等；缺点是略微增加了额外内存，内存大小与平面波道集的个数成正比，将上述成像方法简记为Lrfd-Pre-PLSRTM.
对于规则观测系统，LSRTM、PLSRTM、Pre-PLSRTM方法的计算时间及I/O成
本如表1所示，表中可见：①PLSRTM与Pre-PLSRTM的计算时间小于LSRTM，主要是由于Np≪Ns；②Pre-PLSRTM的I/O成本明显高于LSRTM与PLSRTM，这是因为各个叠前平面波成像结果需分别存储、独立更新.值得一提的是，最速下
降法LSRTM需要两次基于散射波方程的线性正演，其中一次是求取反偏移数据，用来匹配观测数据，另一次是对梯度进行线性正演，用于计算更新步长，加上一次基于伴随方程的逆时偏移，用以计算梯度更新量.
表1 LSRTM、PLSRTM、Pre-PLSRTM的计算时间及I/O成本对比Table 1 Comparison of computation time and I/O costs of LSRTM, PLSRTM and
Pre-PLSRTM计算时间主要变量的I/O成本(单位:×4 字节
数)LSRTMNiter×Ns×(2×Tmod+Trtm)(ng×nt)+(3×nx×nz)+(nt×nx×nz)+(6×nx×nz)炮数据+3×波场变量+正传波场+ (2×成像结果+2×梯度+2×更新
量)PLSRTMNiter×Np×(2×Tmod+Trtm)(ng×nt)+(3×nx×nz)+(nt×nx×nz)+(6
×nx×nz)炮数据+3×波场变量+正传波场+ (2×成像结果+2×梯度+2×更新
量)Pre-
PLSRTMNiter×Np×(2×Tmod+Trtm)(ng×nt)+(3×nx×nz)+(nt×nx×nz)+(6×n
x×nz) ×Np炮数据+3×波场变量+正传波场+ Np×(2×成像结果+2×梯度+2×更
新量)说明:各叠前平面波成像结果需分别存储、独立更新
注：Niter—迭代次数，Ns—地震总炮数，Np—合成平面波个数(≪Ns)，Tmod—单炮线性正演时间，Trtm—单炮逆时偏移时间，nx—模型横向网格大小，nz—模型纵向网格大小，nt—总的时间点数，ng—单炮接受道数.
2 模型测试
本文首先选用简单Salt模型进行Low-rank有限差分法纯qP波线性正演模拟、纯qP波最小二乘逆时偏移成像以及纯qP波叠前平面波最小二乘逆时偏移成像测试，并与拟声波线性正演模拟及相应的反演成像结果进行对比来验证本方法的正确性；随后，将本方法应用于含有不同强度噪声的地震数据成像试验中，进一步测试本方法的抗噪性；在验证了方法正确性及对噪声地震数据有效性的基础上，开展复杂MarII模型成像试处理，以测验方法对复杂地质模型的正、反演适应性，同时，测试了本方法对速度及各向异性参数模型误差的敏感性，突出本方法的优势所在；最后，将本方法应用于实际资料试处理中.下列测试统一采用Cerjan衰减边界条件处理人工边界反射，在成像处理中，使用震源归一化互相关成像条件来补偿中深部成像能量.
2.1 简单Salt模型验证
2.1.1 模型参数
图1a—d分别展示了Salt模型对应的真实纵波速度、偏移纵波速度以及各向异性参数.真实纵波速度用来合成观测记录，偏移纵波速度是通过对真实纵波速度的平
滑获得，用以表征背景速度.根据定义，通过背景速度与真实速度能够求出反射系
数模型(如图1e所示)，这是反演成像试图通过迭代来逼近的理论值，也是成像试
验中用以比较成像质量优劣的重要参考.需要说明的是，在常规的有限差分法波场
计算中，通常需要对各向异性参数进行平滑处理，以提高算法的稳定性，但对于Low-rank有限差分算法，无需平滑各向异性参数，这是因为Low-rank分解更偏向于处理线性相关的传播矩阵，平滑处理只会增加传播矩阵的秩.。