一种稳定的TTI介质逆时偏移方法及其应用.doc

1、一种稳定的TTI介质逆时偏移方法及其应用Yu Zhang, Houzhu Zhang, Guanquan Zhang摘要：非均匀倾斜横向各向同性(丁门)介质中的正演和逆时偏移存在严重的不稳定 问题，尤其是在构造倾斜剧烈变化的区域。基于原始各向异性弹性波方程，我们推 导出一个稳定的了门介质声波方程。特别是，我们推导出了一个垂直横向各向同性(VTI)介质波动方程组，这个方程组与它的弹性对应部分等效。然后引入了旋转史 标系中的向共桅微分算子，用来使TT】介质声波方程稳定。与传统的公式相比，这 个新的方程组不会增加数值复杂度，并且不论是用伪谱法还是舟阶精度有限差分法 进行求解，都能得到一个稳定的解。最

2、后的几个实例表明，我们的方法能够得到稳 定商质量的TTI介质逆时偏移成像结果。1引言逆时偏移(RTM)已经成为了进行复杂构造成像时常用的方法(Baysal等， 1983； McMcchan, 1983； Whitmore, 1983)。它提供了 一种能够有效解决速度横向 剧烈变化的成像方法，并且成像不受倾角的限制。宽方位角采集数据的引入提供 了更好的地下构造照明和更多的方位角信息(Etgen, 2019),在很大意义上这能够 辅助确定各向异性参数，而这些各向异性参数对于高陡各向异性覆盅层下的成像 是很有意义的。覆是在倾斜岩丘侧翼上的泥岩或页岩层可看做是TD各向异性问 题(Alkhalifah和

3、Bednar, 2000)。这种地下的组合关系在墨西哥湾的深水区是很常 见的，尤其是在盐丘隆起形成的小盆地内。对这种各向异性倾斜对称轴的忽略不 仅仅会导致成像模糊、盐体侧翼移位等问题，还会使盐丘底部和盐下构造成像交 形、失意、扭曲等(Huang等，2019)(,各向同性介质和叮】介质的逆时偏移已经得到广泛应用，但是三维非均匀T。 介质中逆时偏移的应用还不成熟。这种不成熟的原因在于推导非垂直对称轴介质 系统的数值公式所遇到的困难以及该公式引起的不稳定问题。相对于各向同性和 VT】的情况，TT1介质逆时偏移的实施也会带来更高的计算成本。Alkhalifah (2000)推导出用于横向各向同性(T1

4、)介质中对压缩波(P波)进 行正演模拟的伪声波方程。这个方程由考虑Thomsen参数，和5的P波频散关系 推导而来，并且假设沿对称轴方向的剪切波(S波)速度为零(Thomsen, 1986； Alkahalifah, 1998)。它包含波场在时间和空间上的四阶偏导数，这带来了一定的 数值计算复杂度。基于同一频散关系，Zhou等人(2019a)得到了由两个方程组 成的VT1介质摘合的二阶偏微分方程组，对于数值求解来说，这个方程组更便捷。 之后，Zhou等人(2019b)将他们的方程组推广到 TO 介质；Lcasge等人(2019) 展示了一些基于该方程组的简单模型三维数值算例。在实际应用中，对复

5、杂TT】介质已经发现了由Zhou等推导的 TO 介质伪声 波方程(2019b)或者它们的等价变分形式(Fowler等，2019)引起的不稳定性问 题(Fletcher等,2019; Zhang和Zhang, 2019) (,为了得到一个稳定的TT1介质逆 时偏移算法，我们需要从基本的弹性波动方程组入手来寻求使TTI介质声波方程 能够包含更多精确物理量的方法。比如，Fletcher等(2019)提出了加入非零S波速 度项的方法来克服这个问题。本文提出了一个简单的与之不同的解决方案，这个 方案首次被报道于此前的一次会议论文集(Zhang和Zhang, 2019)。我们重新推 导了各向异性弹性波方程

6、，得到了一个等价的TH介质声波方程，它在保持波物 理能量的同时也保留了波的运动学特征。基于一个与它的弹性对应部分等效的 叮】介质方程组，我们引入了旋转坐标系下的肯共匏的差分算子来使TT】声波方 程稳定。考虑波传播过程中波场能量的保持，我们证明了新方程组是稳定的。之 后我们用高阶精度有限差分法来求解这个新的TTI声波方程以避免数值频散。数 值试算的例子表明我们的方法提供了一个稳定的TT1介质逆时偏移方法，满足了 生产需求。我们也展示了墨西哥湾一个宽方位角数据的例子来表明TT1逆时偏移 是如何帮助在复杂地质区域进行盐下成像的。2稳定的声波方程传统的各向异性声波方程是从下面的运动学关系中推导出来的(

7、Tsvankin, 1996):«（1 + 2£sin,8（£ - 5）sin2 3 cos? y/(1)式中表示p波相速度，力表示相速度方向相对各向异性对称轴的夹角，入表 示沿对称轴方向的速度，%三4。)。和5是Thomsen无量纲各向异性参数 (Thomsen, 1986)o在式中，我们假设了沿对称轴的S波速度是零来简化这个关系。根据Alkhalifali (2000)的工作，Zhou等(2019b)将TD介质中声波传播的物理过 程进行公式化，得到一个耦合的二阶偏微分方程组。三维情况下的方程组可写作 (Fletcher 等,2019)：L W =(I + 2&

8、quot;仆 + J' + ?)+/,% 01(2)J堤=2( 一加4+/加 + )式中产是拟压力波场，是引入的一个辅助波场来减少数值计算。”表、H五 和“或是与对称轴对齐的旋转坐标系下的二阶基分算子。Hxt =cos2 0cos2 0r + sin2 0cos2 0- + sin2 0-rdr6)广(3)Ql+ sin 2cos2 0-sin0sin2e- cossin 20dxdydydzdzdx=sin，J + cos。J - sin2。(4)。广d.y-oxay(5)= cos- sirr 0_7 + sin- 0siir 0- + cos 6- 。厂 &z2) d2d

9、2d2+ sin 20sirr 0bsin 0sin 26+ cossin 26dxdydydzdzdx式中。是倾角，。是各向异性方位角。旋转坐标系定义为、XV 9r'COSCOS一 sin。,cossin 0sin 0cos6COS0sin 8sin。(6)下面的关系式可以用来简化计算这里我们注意到，与P波运动学关系（1）相比，方程组（2）有两个本征值:(8)皿)= L + £Sy/± JQ + 2-sin犷/一8()一6)s2印cos?犷片一 5 一式中的减号引入了对称轴方向之外的非零S波能量，在正演和偏移中我们把这些 能量视为噪音。由于这些s波的存在，我们需要

10、确保£2 5来避免式（8）中出现 不符合物理规律的S波速度（Aki和Richards, 1980）。即使在£之5的限制条件下，对于在对称轴方向存在强速度变化的非均匀TO 介质，方程（2）和它的变体依然会引起计算的不稳定（Zhang和Zhan3 2019）。 为了解决不稳定问题，我们从VT1介质弹性波方程组开始探究，方程中各向异性 对称轴是垂直的。通过设置4 =。, Hooke定律可以写作（Thomsen, 1986）:1 + 2£=01 1 + 2£J1 + 251 + 2£ 、/1 + 26 口 '1 + 2s Jl + 2s e22

11、i I"式中是密度，叫和%分别是应力和应变张量矩阵的对角元素。它们通过下面两 式与位移向量（、,、,七）联系起来：(10)(H)对式（9）中的三个方程分别对y和z求一阶微分，应用关系式（10）和（11）,可以得到：6% _ 1 6 吊尸=万曲 ©2“ 1 q< =dr p dy d%.1 dJ .dr p dz如果我们定义：(12)加g噜+ (1 + 2曙+ 423誓) 0*1 + 2£炭 + (1 + 2£格 + 八 + 2。鲁) 由生+却(13)(14)偏移中我们通常假设密度是常量。在这个假设下，式(14)可以写为:(15)据我们所知，式(15

12、)最初由Duvcncck等人(2019)发表，其中和，.分别是水 平和垂直应力分量。事实上，从式(9)和定义13中，我们能够看出P = b“=b22 和 r = %。在式(2)的VTI介质情况下，如果我们定义)=广+ 7,则式(2)可以写为：1» =，1 + 2£ 0苏+ 石忑才011 + 26 1J n°平a2 r)(16)只有当5在空间上为一定值时，式(15)和(16)是等价的。在这种情况下，我 们有 = = )+" r = J1 + 2 Jr(17)式(15)和(16)的不同在于等式右侧的第一个矩阵。在式(15)中，矩阵l + 2f J1 + 26

13、 )3 + 26(18)&dy dz )在物理约束条件£>6(19)下，是自共规且非负定解的。可以证明，方程(12)和(15)保持了波传播的物理能量，并且保证了波传 播过程中的稳定性。由于方程(12)和(15)与弹性波方程(9)是等价的，所以 这并不奇怪；它们本就应该如式(9) 一样保持波传播的能量。另一方面，由于式 (16)是仅仅由运动学关系推导出来的，它的动力学特征值得怀疑，因此可能会 导致不稳定问题。所以，在RTM计算过程中，我们应该选择方程(12)或(15) 来进行波场计算，而不要使用方程(16)。另外，VTI介质声波方程(15)比它的 弹性等效方程(12)更简

14、洁，这使得实际应用中可以大量节省计算时间，所以生 产中式(15)更有吸引力。下一步，我们通过将坐标系旋转到TTI介质各向异性对称轴的方向，将式(15) 由叮】推广到TO 的情况，得到1 d2 pl + 2e、J1 + 25(20)为了达到稳定性要求，坐标旋转后的差分算子G。，G6和应该如二阶微分算子(或，J和三)一样，是自共扼且非负定解的。例如，我们可以将这三个算 次 守才子写为:右穴=(5),(5)q=(oj(2)在式(21)中我们令D- =coscos + sincos-sin dxdy dz(21)(22)(OJ是算子心的转置。对于式(20),如果我们定义能量为E(，)=UT快长(前一图

15、苏)其中ap + brbp - ar% =，、 、(24)能够证明，在求解TT1声波方程(20)的过程中，能量E(f)总是守恒的，这就是 说乡刖=0(25)dt因此，对于TTI介质声波正演和偏移，方程(20)是稳定的。3应用方程(20)可以通过下面的时间方向上4阶有限差分公式进行数值求解(Etgcn, 1986):(26)式中加是离散时间步长，/是密度矩阵，并且有(27)Dispersion o( Finite Difference Methods图1不同中心有限差分算法逼近一阶微分的频散曲线(八k)=k )c黑 色曲线代表理论上的频散关系；蓝色曲线是40阶传统有限差分算法 的频散关系；红色曲

16、线是本文中用来计算实例的有限差分算法。空间差分算子矩阵。可以通过伪谱法(Kosbff和Baysal, 1982)或者显式高 阶有限差分方法(Etgcn, 1986)计算。由于式(26)中的差分算子都是1阶的， 所以通常选用交错网格(Viricux, 1986； Sacn空r等，2(X)0)有限差分方法。然而, 我们发现这样的方法需要将模型参数内插到不同的计算网格上，这样就会引起数据不准确。另外，为了得到高的精度，交错网格方法通常假设在三个空间方向上一致采样，这降低了它们应用的灵活性，增加了计算成本。因此，我们决定在所有三个空间方向上都采用中心差分来计算一阶导数。伪谱法和高阶显式有限差分方法是已

17、发展的用来求解TTI介质声波方程的两 种方法。由于伪谱法需要对地震数据进行多次三维傅里叶变换，这导致了巨大的计算量，所以在生产中为了提高效率，我们更倾向于使用有限差分法。然而，使用中心差分来逼近一 阶微分的方法会在高 波数区产生数值频 散。由于Nuquist波数 处奇点的存在，传统 中心有限差分的收敛 性是很慢的。如图1 所示，即使使用40阶 的传统有限差分来逼 近一阶微分，在 80%Nyquist波数以上 的区域仍然出现了强 数值频散(蓝色曲 线)。为了提高生产中 TT1逆时偏移的精度，我们发展了极高阶的有限差分算法，对于任意空间方向，它都能够使得微分计算的精度达到92%的Nyquist波数

18、（红色曲线）。这是我们下一步用来计算例子 所用的算法。在我们的TT1有限差分算法中，对&我们仍然需要严格应用Courani - Friedrichs-Lewy （CFL）条件来避免数值模拟中的不稳定问题（Courant等，1967）。 另外，/的选取需要保证每个波长内至少有5个采样点，这样可以有效减少时间 方向上的数值频散（Zhan&等，2019）o4数值算例第一个例子是三维脉冲响应测试。使用的模型是为墨西哥湾实际数据项目而 建立的各向异性模型。测试中，所用的网格尺寸为= &v = 1230,偏移最高频 率为15Hz。偏移中，我们选择&： = 58以及&

19、 = 6.69ms来确保数值稳定性和抑 制频散。使用传统TH声波方程（2）进行的逆时偏移导致了 一个不稳定的成像结 果（图2a）o从这次早先的测试中，我们认识到克服这种不稳定性是将TT1介质 逆时偏移投入到实际生产过程中的一个很重要的步骤。我们用方程（20）实施了 逆时偏移之后，得到了一个稳定的成像结果（图2b）。这个例子表明，我们提出 的TT1逆时偏移方程能够解决计算中的不稳定问题。图2冲击响应测试。（a）基于公式（2）的TTI逆时偏移引起数值不稳定；（b）采用公式（15）的TTI RTM 提供了一个解决数值不稳定问题的方案。ff5 ,号一二一.三WW-图3 （a）为某小盆地内各向同性偏移得

20、到的共成像点道集°为了便于观察，来自互相垂直测线的两个 常生会期息云本周一常佳的西州 u里缶立线会嗝（力伽出的物幅 六的她推）.出局Ea95001900028500b) 1000,CZ意，'(xXi -3000 -9600(19000128500,9000 1r7000ifsMiMOdrsooDistance (m)96001d000-28500Disunce (m)图4 （a, b）为Garden Banks地区各向同性模型的两种方位角（90°和0° ）的逆时偏移结果；（c, d）为 该地区TTI模型的两种方位角逆时偏移结果。第二个实例是墨西哥湾Gard

21、en Banks区域的一个将两条互相垂直的地震测线 结合起来的双方位角项目。该项目的目标是利用双方位角数据提供的良好的地下 地质结构照明来改进复杂地质区域的盐下成像。该项目最初采用各向同性速度更 新进行处理。经过双方位角层析成像的多次迭代之后，我们用控制波束偏移（CBM） 方法（Tin&和Waiig, 2019）对比了从两种方位角数据得到的共成像点位置处的 道集。具有负偏移距的90°方位角道集显示在左侧，在此处可见形成了 “超道集” 或“蝴蝶状道集”（Dewey等，2019）；正偏移距的0°方位角道集显示在右侧。在 陡倾角盆地区域，每条测线的共成像点道集（ClGs）

22、都不是平的。由图3a可见， 在共成像点位置处，0°方位角测线数据表现出过校正剩余时差，这需要更高的速 度来拉平道集。90°方位角测线表现出欠校正剩余时差，需要降低速度来拉平道 集。对各向同性偏移来说，这样的过校正或欠校正问题很难找到合适的解决方案。 我们对该项目应用VTI介质成像工作流程，依然出现了类似的现象。在具有陡倾 角盐体侧翼的小盆地内，速度随方位角的变化可以用TD各向异性来解释。因此, 处理小组决定将TT1各向异性应用到速度模型的更新中，用以拉平每个方位角的 共成像点道集。图3b显示了应用 TO 层析成像后得到的共成像点道集。与图3a 的各向同性层析得到的道集相比，

23、TO层析成像改善了道集的平滑度，消除了 0° 和90°方位角之间测线的不一致的剩余时差。图4是各向同性和TH逆时偏移成 像结果之间的对比。对于所有两个方位角，我们提出的TD逆时偏移算法能够得 到稳定的成像结果，并且改善了盐下区域的成像聚焦，同时改变了沉积层位同相 轴的位置。采用TF1逆时偏移的盐下成像结果表现出了更简单的构造，表明了门模型具有更多的地质意义。臼2019年以来，TTI逆时偏移已经广泛应用于墨西哥湾宽方位角(WAZ)项 目。图5展示了墨西哥湾*alkcr Ridge地区宽方位角数据的例子。我们在具有倾斜 侧翼的小盆地区域测试了上述TO逆时偏移方法。为了更好地对比

24、两种方法，我 们分别应用各向同性和各向异性层析成像技术来建立各向同性和TD介质速度模 型，与之对应，分别应用各向同性逆时偏移和TO逆时偏移方法对数据进行偏移 处理。从图中可以明显看出，与各向同性逆时偏移相比，TTI逆时偏移能够产生 一个更好的成像结果，尤其是在沿着陡倾角盐体的侧翼和小盆地下侧的盐下区域。 这个例子表明，通过应用宽方位角数据中的方位角信息，TO逆时偏移能提升传 统各向同性逆时偏移的地震成像精度。Distance (fn)图5 (a)为墨西哥湾论Iker Ridge地区各向同性模型宽方位角数据的逆时偏移结果；(b)为该地区TTI模 型宽方位角数据逆时偏移结果o5讨论传统TTI介质声

25、波方程是由频散关系推导出来的，如式(Do仅当模型参数 是常量的时候，方程里的微分才是有效合理的；直接将这些方程应用到一般TH 介质会引起不稳定问题。这告诉我们各向异性是一种弹性现象，它的声波近似不 能仅仅是一个频散关系，而需要更多的物理关系约束。我们通过引入式(20)中 的自共施算子矩阵，确保波传播过程中能量守恒，避免了正演和偏移中的异常振 幅，在一定程度上解决了上述问题。文献中也可以看到一些其他的解决方案。比如，FHchcr等人(2019)提出的 加入非零SV波速度项的方法，有(28)仅=4+2f 1/京+ ”芹o Yp 正祈一 11 + 26认0上网 0 Vp-r制0其中，%是SV波沿对称

26、轴方向的速度，»和，是拟压力波场。根据Duvcncck和 Bakker (2019)的研究，方程(28)能够改善稳定性,但是不能保证一定稳定。Liu等人(2019)提出如下的公式：l +% 0 x|P(29)西I a 。 "J 这个公式看上去类似于式(20),但是由于等号右侧的微分算子矩阵不是自共施的, 所以它并不能保证在一般TT】介质中的稳定性要求。我们注意到Duvcncck和Bakker (2019)最近提出了稳定的具有向共辄旋转差 分算子的TT1声波方程组，可写为了布*44 A；2丫1 J He A 耳 Bj 入 r )(30)二闺岛+与&) 加。工1-(用内 2 + &网),2=6瓦-心附32(用内 3 + &i&3)_L_£lf /?