慢度平方三角網格模型VTI介質二維QP波各向異性立體層析反演

周嘉欣， 楊鍇， 邵煒棟

1 中煤科工集團常州研究院有限公司， 江蘇常州 213015 2 天地(常州)自動化股份有限公司， 江蘇常州 213015 3 同濟大學海洋地質國家重點實驗室， 上海 200092 4 阿里巴巴集團， 杭州 311121

0 引言

各向異性在地下介質中廣泛存在，在地震波傳播中各向異性會同時影響信號的相位和振幅.對于大偏移距數據，不考慮各向異性會導致速度估計發生很大誤差.這促使近年來學術界和工業界在研究基于地震波的反演、成像以及儲層建模技術時，必須將算法從各向同性拓展到各向異性(Tsvankin et al.，2010).因此面向各向異性參數估計的層析成像技術成為現代地震反演與成像技術中的一個重要組成部分.然而，由于各向異性層析成像方法引入了更多未知參量，使得參數估計的不確定性和復雜度都大大增加，因此長期以來一直是一個非常有挑戰的研究課題.

基于走時的各向異性層析成像已經有大量研究工作(Chapman and Pratt，1992；Zhou et al., 2008; Wang and Tsvankin, 2013；盧明輝等, 2005； 劉玉柱等，2014；黃光南等，2015).20世紀末以來，一種全新的層析成像算法——立體層析，開始為學術界和工業界所知.立體層析將局部相干同相軸在炮集與檢波點道集內的射線參數(或P參數)水平分量，以及炮點坐標和檢波點坐標都納入到反演的數據空間之中，重新定義了層析反演的數據分量和模型分量，使得立體層析成為射線類層析反演方法中唯一一種可以同時反演速度結構、反射點位置與反射層局部形態的方法(Billette and Lambare, 1998).國內外學者在各向同性立體層析的理論和實踐方面已有很多研究工作(Chauris et al., 2002; Billette et al, 2003; Alerini et al., 2007, 2008; 倪瑤等，2013, 2014；任麗娟等，2014；李振偉等，2014).Prieux等(2013)利用立體層析為全波形反演建立了優秀的初始模型.Li等(2015)將塊狀約束加入立體層析FRECHET導數矩陣的建立過程，使得反演結果更具有地質一致性.王宇翔等(2016)將結構張量算法用于P參數水平分量的搜索，實現了高密度的二維立體層析反演.楊鍇等(2016)基于非降階漢密爾頓算子實現了三維直角坐標系下的立體層析反演方法.Xing等(2017)并將其用于實際資料處理.Xiong等(2017)將結構張量用于成像域，將人工拾取與運動學反偏移相結合，制定了更為合理的成像域數據點拾取策略，形成了一種更為實用的兩步法策略.

各向異性立體層析方面的研究則相對較少，Barbosa等(2008)討論了橢圓與非橢圓各向異性介質中PP波立體層析的FRECHET導數求取問題.Nag等(2010)推導了基于PP/PS波數據空間、適用于一般各向異性的立體層析FRECHET導數公式.Tavakoli 等(2017)則將上述工作在伴隨狀態法的理論框架下進行了初步實現.

前人實施各向異性立體層析方法各異，但尚未有人在三角網格模型下予以討論.使用三角網格模型是考慮到它一種非常稀疏的模型描述方式，利用三角網格來描述模型空間將可以大幅壓縮層析矩陣的規模，在提高求解精度的同時壓縮了計算成本.Yang等(2017)，邵煒棟等(2017)基于慢度平方三角剖分模型實現了二維各向同性立體層析.他們采用地質界面約束下的慢度平方模型三角網格剖分描述地質模型(Hale and Cohen, 1991)，能夠自然的在模型描述中引入地質界面，該建模方法能夠以最稀疏的方式描述塊體內的緩變與劇變的波阻抗界面，從而極大的壓縮層析反演的模型空間.考慮到各向異性對地震波運動學特征的影響屬于次級效應，如果想要保證各向異性參數的反演精度，采用盡可能稀疏的模型表達方式不僅是必要的，而且在實踐中更是必須的.慢度平方三角網格模型描述的另一個好處是射線追蹤有解析解，使得計算更為高效，從正演的角度對節省計算成本也做出了貢獻.

在Yang 等(2017)、邵煒棟等(2017)工作的基礎上，本文基于考慮界面的各向異性射線擾動理論(Farra, 1989),將模型空間從各向同性擴展到VTI介質，并逐一求取了反演所需的各個數據空間分量對任一模型空間分量的偏導數，從而建立起慢度平方三角網格剖分下二維VTI介質中的qP波各向異性立體層析所需的FRECHET導數矩陣.本文所使用的各向異性參數化方式為ε以及η=δ-ε，也就是說需要反演的各向異性參數為ε和η.背景介質定義為橢圓各向異性VTI介質，即η=0的情況.由于在慢度平方橢圓各向異性介質中的射線追蹤有解析解，基于射線擾動理論從背景橢圓各向異性介質擾動到非橢圓各向異性介質時能夠獲得半解析解，這使得射線追蹤以及FRECHET導數的建立變得簡潔高效.可以認為是本文方法最顯著的一個優勢.

本文展示的敏感度測試表明，慢度平方三角剖分模型下建立的立體層析FRECHET導數矩陣中的元素具有較好的敏感度，可以用于二維VTI介質QP波各向異性立體層析反演.由于模型空間參數的增加，對所有參數進行同時反演變得相當困難，因此本文設計了相應的分步驟反演工作流程，首先反演除了各向異性參數之外的其他模型空間分量如慢度、界面信息、反射點位置等，在獲得比較可靠的上述初始模型空間信息基礎上，再反演各向異性參數獲得各類模型參數信息，來保證各向異性參數反演可以得到穩定的結果，理論算例證實了上述觀點，為后續的實際應用奠定了理論基礎.

1 各向異性射線理論簡介

為保證文章的自洽性，首先對各向異性介質射線理論做一簡介.Christoffel方程是各向異性射線理論研究的起點.

1.1 Christoffel方程

det(Γjk-Gδjk)=0,

(1)

其中G(x,p)=1，是方程的特征值，對應的方程特征向量為相應波的極化矢量.Γjk是Christoffel 矩陣的元素,Γjk=piplcijkl，密度歸一化后的彈性參數aijkl=cijkl/ρ.其中cijkl是各向異性彈性參數,pi=?T/?xi表示相慢度各分量.T為走時，ρ表示密度.

由于彈性剛度常數的對稱性，可以將這個系統寫成Voigt矩陣形式.彈性剛度矩陣元素帶有兩個下標，即：cijkl=Cm n(m,n=1,2,…,6).其對應關系如下：11→1，22→2，33→3，23和32→4，13和31→5，12和21→6.由此，原本3×3×3×3的彈性剛度系數矩陣cijkl可以用6×6的矩陣Cm n表示.本文將使用前文提及的密度歸一化后的彈性參數aijkl及其對應的6×6的矩陣Am n.

常見的橫向各向同性介質(TI介質)具有柱對稱軸.在這種介質中，方程(1)存在三個獨立的解:一個qP波和兩個不同方向的qS波.由于地球內部大部分沉積巖都表現出TI介質的各向異性，因此以對稱軸垂直時的VTI模型為基礎的各向異性TI介質研究顯得尤為重要.該模型常常被用于描述由周期性的薄互層、巖石內部結構和平行排列的微裂隙引起的各向異性.本文研究以對稱軸沿著z軸方向即垂直于各向同性平面介質(VTI介質).在此類介質中方程(1)中三個特征值可以解析表達為(2)式：

(2)

其中

Farra(1989)據此定義了各向異性介質中的Hamiltonian：

(3)

(4)

根據(3)式所示的各向異性Hamiltonian，可得對應的運動學射線追蹤方程：

(5a)

(5b)

1.2 慢度平方三角網格中VTI橢圓各向異性介質的二維qP波射線追蹤

本節回顧Farra(1989)在背景介質為VTI橢圓各向異性下導出的qP波射線追蹤方程.由于該方程形式簡潔，使得在背景介質求取關于ε的FRECHET導數變得相當容易.

根據Thomsen參數表征有：

(6)

(7)

結合(2)式，可以得到如下關系：

(8)

+4(δ-ε)A33.

(9)

(10)

其中

(11)

在橢圓各向異性介質(δ-ε=0)情況下，顯然有ΔG=0，這時將(10)式代回(3)式可以得到橢圓各向異性介質下的qP波Hamiltonian：

(12)

顯然可以將ΔG視為由η=δ-ε的小擾動造成的變量，將(11)式代入(3)式，得到對應的由橢圓各向異性擾動至非橢圓各向異性的ΔH:

(13)

令H=HqP0+ΔH，即可得到本文使用的VTI介質中的qP波Hamlitonian:

(14)

將(12)式所示的橢圓各向異性背景下qP波Hamiltonian重寫為

(15)

(16a)

(16b)

顯然，方程(16)的解析表達式如下：

+x(τ0),

(17)

走時可由積分變量τ與走時T的關系dT=u2du得到：

(18)

如果僅僅在橢圓各向異性介質下實施立體層析反演，通過(17)、(18)式即可容易求得數據空間中X坐標，走時T，射線參數P關于ε的偏導數：

(19a)

(19b)

(19c)

(19d)

2 慢度三角網格模型下VTI介質二維QP波立體層析的數據空間與模型空間簡介

考慮到文章的自洽性，首先利用圖1展示常規二維各向同性立體層析中數據空間與模型空間的各個分量.對于二維空間的一個射線對而言，假設地表水平，那么地表能觀測到的數據為d=(Sx,Rx,Tsr,

圖1 二維立體層析數據空間各分量與模型空間各分量Fig.1 The model-space components and the data-space components in 2D stereotomography

Psx,Prx),其中Sx,Rx為炮檢點的X坐標；Tsr為這個射線對的總走時；Psx,Prx為炮檢點處的射線參數水平分量.這些數據分量與地下的模型空間m=(x0,z0,θs,θr,V)有關，其中x0,z0為地下反射點坐標,θs,θr代表從炮點與檢波點一側出射的地下散射角,V代表地下介質速度信息.其中共炮點道集和共檢波點道集內的局部相干同相軸的斜率等同于檢波點處和炮點處慢度矢量的水平分量(Billette and Lambaré，1998).

而在慢度平方三角網格下VTI介質中討論二維QP波各向異性立體層析時，數據空間分量不變，其模型空間變為：m=(x0,z0,θs,θr,xn,zn,u2,ε,η)其中xn和zn是構成界面的節點的橫、縱坐標,u2表示慢度平方.ε和η是之前定義過的各向異性參數.數據空間的殘差Δd與模型空間殘差Δm之間的關系由射線擾動理論建立，這里可以簡寫為

Δd=FΔm.

(20)

(20)式代表了一個矩陣方程，其系數矩陣F即為立體層析FRECHET偏導數矩陣.在慢度平方介質下的完整形式可以表達如下：

(21)

在這個矩陣中，顯然有：

(22)

其次，根據費馬原理，走時T在速度結構、激發點與接收點位置都固定時，已經是最小走時.因此可以認為散射角對走時T不存在一階擾動影響.即可認為

(23)

由于在該理論方法下，層內各向異性參數被設為常數(即每個三角形網格中也為常數)，因此射線相慢度p對各向異性參數并沒有一階影響，即

將上述關系代入(21)式，可以得到簡化后的三角網格立體層析矩陣(24)式.(24)式中σ的含義為針對各類數據的尺度加權因子，用于平衡不同類型的數據對于立體層析反演的誤差泛函的貢獻比例，避免由于物理量綱不同造成數據類型量級之間差別過大的問題，參數的大小需要通過做具體的測試得到.

(24)

立體層析反演的理論基礎是射線擾動理論(Farra and Madariaga,1987).射線擾動理論的重要性在于，它在傍軸射線追蹤和漢密爾頓傳播算子的基礎上建立了速度擾動與描述射線傳播的相空間各參數擾動量之間的一階線性關系.三角網格剖分由于在模型描述中引入了界面，因此需要在擾動理論中考慮界面擾動對立體層析數據空間的影響.Farra等(1989)提出的考慮界面擾動的射線擾動理論恰恰是建立三角網格立體層析FRECHET導數矩陣所需要的.Yang等(2017)，邵煒棟等(2017)基于該理論建立了各向同性介質下三角網格立體層析FRECHET導數矩陣.由于各向異性參數的加入，本文不僅僅利用傳統的射線擾動理論，同時也結合解析公式引入了數據空間分量關于各向異性參數的FRECHET導數.第3、4節根據上述認識逐一建立三角網格VTI介質QP波二維立體層析中數據空間分量對模型空間分量之間的一階線性擾動關系.

3 三角網格慢度平方模型下VTI介質二維QP波各向異性立體層析FRECHET導數求取

在上文列出的FRECHET導數矩陣中，前9列對應于各向同性二維三角網格立體層析層析FRECHET導數，其推導方法請參考邵煒棟等(2017)，這里不再重復.本文重點討論與各向異性參數ε和η相關的這部分導數，即(24)式中的最后兩列如何導出.

3.1 基于射線擾動理論推導數據空間關于ε和η的一階擾動關系

以下基于射線擾動理論推導數據空間各分量關于各向異性參數ε和η的一階擾動公式.我們首先定義一個新的Hamiltonian表達式：

(25)

按照射線擾動理論可以寫出

(26)

(27)

公式(26)的另一種更為簡潔的寫法為

y(τ)=y0(τ)+Δy(τ).

(28)

為推導數據空間關于各向異性參數ε、η與初始相空間擾動的一階擾動關系，不妨令初始相空間擾動為Δy0=(Δx0, Δp0)，需要考察的擾動為各向異性參數的擾動Δε和Δη，慢度的擾動Δu.對(28)式進行微分即有

(29)

利用傳播矩陣(Gilbert and Backus，1966)，可以得到式(29)的傳播矩陣解：

(30)

其中右端第二項中ΔB的計算如下：

(31)

公式(30)中的傳播矩陣P0在橢圓各向異性介質中有很簡潔的形式，因為橢圓各向異性背景介質中的擾動射線方程有解析解，其形式為：

Δx(τ)=Δpx0(τ-τ0)+Δx(τ0),

Δz(τ)=Δp0(τ-τ0)+Δz(τ0),

Δpx=Δpx0,

Δpz=Δpz0.

(32a)

將(32a)式整理為與傳播矩陣的乘積形式

(32b)

即可看出其傳播矩陣P0矩陣的形式確實非常簡單.

基于(31)，(32)式展開(30)式即可獲得地表數據空間中觀測坐標X、射線參數P關于各向異性參數的一階擾動關系，如(33)式所示.由于本文將一個地質塊內的各向異性參數設定為常數，在這種情況下地表觀測的P參數對各向異性參數沒有敏感度.

Δz(τ)=Δz(τ0)+Δpx(τ0)+(τ-τ0)Δpz(τ0)

(33a)

Δx(τ)=Δx(τ0)+Δpx(τ0)+(τ-τ0)Δpx(τ0)

(33b)

Δpx(τ)=Δpx0(τ0), Δpz(τ)=Δpz0(τ0).

(33c)

因此很容易寫出地面觀測時，坐標X、射線參數Px關于各向異性參數η的Frechet導數為

(34a)

(34b)

而X、Px關于ε的偏導數，容易看出其形式與(19a)(19c)式完全相同，這里不再重復.

關于走時T、η的一階擾動關系，Farra等(1989)定義了如下走時擾動公式：

(35)

因此很容易寫出T關于η的FRECHET導數為

(36)

從(35)式可以看出，走時ΔT與ε無關.因此ΔT關于ε的偏導數依靠背景橢圓各向異性公式即可求得，與(19d)式相同，這里也不再重復.

至此當射線在一個地質塊內傳播時，立體層析數據空間關于ε和η的偏導數已經全部得到.

3.2 非橢圓各向異性VTI介質中的半解析射線追蹤及上述一階擾動關系的數值驗證

由于非橢圓各向異性介質中的走時計算公式可以理解為背景橢圓各向異性計算公式(18)與擾動量計算公式(35)之和(Farra et al., 1989)，因此有如下走時半解析公式：

(37)

顯然坐標與慢度矢量在非橢圓各向異性介質下的半解析射線追蹤方程亦可以同樣的形式寫出：

(38)

(38)式即為本文采用的慢度平方VTI介質中的各向異性射線追蹤方程.本文方法實現下的射線追蹤方程在層內計算時可以直接得出層位出射點的射線坐標與慢度的解析解，加之以層位之間的連續性條件(Shao et al.,2017)即可實現模型內完整的射線追蹤.

圖2 黑色曲線展示了橢圓各向異性介質下的射線追蹤，灰色曲線為在橢圓各向異性上加以δ的擾動的小擾動計算得到的射線.其中=1/3000/3000 s2·m-2; kz=-1×10-8;ε= 0.3; A33=9.0; A44=5.19; Δδ=-0.2; Sx=500 m; Sz=0.4 m; θ=30°Fig.2 The dark curve shows ray tracing in elliptical anisotropic media, when the gray curve shows the perturbed ray with small perturbation δ in elliptical anisotropic media. Parameters are given as =1/3000/3000 s2·m-2;kz=-1×10-8;ε=0.3; A33=9.0; A44=5.19; Δδ=-0.2;Sx=500 m;Sz=0.4 m; θ=30°

圖3 半解析射線追蹤公式(38)精度測試紅色點線表示由(38)式計算的射線路徑，而綠色實線則是表示常規龍格庫塔方法解得的射線追蹤.測試參數 = 0.11;2ε2=0.2;η=δ-ε2=0.1.Fig.3 Precision test of semi-analytical ray tracing equationThe red dotted line represents the ray calculated by equation (38), when the green solid line represents the ray calculated by Runge-Kutta methods. Parameters are given as = 0.11;ε=0.2; η=δ-ε=0.1.

4 射線穿越多個界面時FRECHET導數推導及其敏感度測試

4.1 射線穿越多個界面時FRECHET導數推導

如上一節所述，公式(19)、(36)表達了當射線在一個地質塊內傳播時，立體層析數據空間關于ε和η的偏導數.然而地下介質一般而言不可能是單個地質塊組成的，也就是說射線一般會穿越多個地質界面達到地表.因此公式(19)、(36)還需要做與通過界面有關的某種修正才可以用于多層三角網格模型.這個問題并非是慢度平方各向異性立體層析所獨有，邵煒棟等(2017)在建立慢度平方各向同性立體層析FRECHET導數時就已經仔細考慮過.這里簡要回顧他們的處理方式，并將其用于本文在射線穿越多個界面時的各向異性參數的FRECHET導數求取.根據邵煒棟等(2017)中的推導，過界面后射線擾動量與過界面前傍軸射線末端的擾動量關系如下：

(39)

(39)式代表了射線過一個界面時，根據矩陣Π和矩陣T即可求得過界面之后的過界面之后的擾動量.其中Π矩陣表示計算中心射線到達地表時對應的傍軸射線的信息，而T矩陣則表示考慮Snell定律和參數連續性后得到的過界面后的傍軸射線信息.在計算下一層的傍軸擾動量時，便可把此次擾動量作為下一層的初始擾動量進行計算.各向異性下這部分的推導與各向同性一致，但具體計算時需要代入式(25)所示的各向異性Hamiltonian，T矩陣和Π矩陣的具體形式請參考邵煒棟等(2017).

由(39)式，可以將射線穿越n個界面時(19，36)式中X關于各向異性參數的偏導數改寫為：

(40a)

(40b)

走時T并不屬于(39)式定義的射線相空間元素，因此在射線穿越多個界面時其偏導數形式并沒有變化.

4.2 與各向異性參數ε和η有關FRECHET導數敏感度測試

在得到三角網格二維QP波各向異性立體層析所需的數據空間各分量與模型空間各分量之間的一階線性擾動關系之后，以下基于三角網格理論模型中的一個射線對信息，來測試不同數據分量與不同模型分量之間的敏感度.圖4a展示了一個三層慢度平方模型，射線穿過2個界面，3個地層，共5個節點.數據空間信息依然通過射線正演獲得，射線的出射點在(1.8 km，2 km)處，在z=0 m處可觀測到數據空間(Sx,Rx,Tsr,Psx,Prx).由于邵煒棟等(2017)已經展示了各向同性下的數據分量與不同模型分量之間的偏導數敏感度，這里不再重復.這里僅展示數據空間分量Psx,Tsr對與各向異性參數相關的模型空間分量ε和η的偏導數敏感度測試.

圖4b，4c，4d，4e展示了數據空間分量炮點坐標Sx和走時T關于各向異性參數ε,η的偏導數測試.在初始參數上加以±20%的擾動，從而計算得到這些曲線.以圖4b為例，偏導數?Sx(ε)/?ε是由傍軸射線追蹤算得，可以清楚看到直線Sx(ε0)+(?Sx(ε)/?ε)Δε與在真實介質中射線追蹤得到的曲線Sx(ε0+Δε)相切,圖4證明了Sx,T和ε,Δη在模型空間中的敏感性是存在的，根據(40)、(19)式計算得到的FRECHET導數結果是正確的.

圖4 敏感度測試模型以及三角網格立體層析中Sx, T關于ε和η的敏感度測試(a) 模型和射線對； (b) Sx/ε; (c) Sx/η; (d) Tsr/ε; (e) Tsr/η.Fig.4 Sensitivity test of Frechet derivatives for some selected data components (Sx, T) and selected model components (ε,η) of triangulated stereotomography(a) Test model and a ray pair; (b) Sx/ε; (c) Sx/η; (d) Tsr/ε; (e) Tsr/η.

為測試多參數聯合反演的可行性，這里采用了模型分辨率矩陣進行驗證.模型分辨率矩陣是分析多參數反演耦合特性的有效工具，其計算公式如(41)所示.其中F為Frechet矩陣，λ可取0.將R歸一化后，R如果接近單位陣I，則說明各個模型空間之間是解耦的.如果對角線外有值，便可以根據橫縱坐標的信息查看具體的模型分量之間的耦合性.

R=[FTF+λI]-1FTF.

(41)

選取一個簡單的兩層模型，在(6,2.99)km處分別以140°、200°出射兩條射線，計算出數據空間關于模型空間各個分量的Frechet矩陣進行上述運算，再套用公式(41)即可得到如圖5b所示的模型分辨率矩陣.注意圖5b從左到右依次為界面的兩個節點的縱坐標；兩個反射點的橫、縱坐標；兩個散射角、一層的背景慢度，兩個各向異性參數，從上到下有同樣的參數順序.可以看出該矩陣呈現出較好的對角化程度.值得關注的是最后兩個各向異性參數顯然在本文方法下可以進行較為穩定的反演，這也證明了最后實現多參數聯合反演的可行性.

圖5 模型分辨率矩陣(a) 用于測試模型分辨率速度模型; (b) 模型分辨率矩陣測試結果.Fig.5 Model resolution matrix(a) Test Vp0 model of model resolution matrix; (b) Result of test of model resolution matrix.

5 三角網格立體層析理論數據反演測試

5.1 三層簡單模型各向異性參數反演測試

由于模型空間加入了各向異性參數，使得直接進行多參數同時反演變得尤為困難.本節首先設計一個簡單三層模型的理論試驗，且只反演各向異性參數.用來檢驗在速度、節點等其他模型空間分量位置都設置為正確值的情況下，各向異性參數的反演精度.本次實驗共計使用了40根射線，并由這些射線來生成本次反演所需的的數據空間分量(Sx,Rx,T).圖6a,6d分別展示了真實的各向異性參數模型ε和η，圖6b,6e為兩者所對應的初始模型，均設為0(即為各向同性)，而圖6c,6f則為反演迭代5輪后得到的最終各向異性參數模型ε和η.

圖6g展示了泛函下降曲線，反演第二輪就已接近真實值，證明了之前推導的關于ε和η的偏導數的正確性.然而，對于任何一種各向異性層析反演方法來說，同時反演所有模型參數都是非常困難的.圖7展示了當界面信息不準時，即便速度為真實值時各向異性參數的反演結果.

圖6 三層模型各向異性參數反演測試(a) 真實ε模型； (b) 真實η模型； (c) 初始ε模型； (d) 初始η模型； (e) 5輪迭代反演后得到的ε模型； (f) 5輪迭代反演后得到的η模型； (g) 泛函下降曲線.Fig.6 Inversion test of three-layer model(a) True ε model； (b) True η model； (c) Initial ε model； (d) Initial η model； (e) ε model after 5th inversion； (f) η model after 5th inversion； (g) Decline of the cost function during the iterations.

根據圖7的結果可以看出，層位信息對反演結果有著非常重要的影響，若將較差的層位信息直接用于與各向異性參數的聯合反演，得到的結果將不符合精度要求.同理，如果VP0不準確時，對各向異性參數的反演精度同樣會產生負面影響.為了保證其穩健性，設計一個多步驟反演工作流程是非常必要的.本文給出了一個實用的二維三角網格下各向異性qP波立體層析的基本工作步驟：

圖7 層位信息不對時反演得到的結果(a) 反演得到的ε模型； (b) 反演得到的η模型.Fig.7 The inversion when other information is incorrect(a) Final ε model； (b) Final η model.

(1)對原始數據做疊前深度偏移(PSDM)；

(2)在初始PSDM結果上拾取層位信息，通過運動學反偏移得到立體層析數據空間；

(3)利用拾取的層位節點信息建立三角網格模型，初始各向異性參數為常數.

(4)基于手頭的數據空間信息，進行常規二維各向同性三角網格立體層析反演，得到除各向異性參數以外的參數的初始模型；

(5)固定其他參數不變，單獨反演慢度平方u2,ε和η來作為各向異性參數的初始模型；

(6)通過步驟(4)和(5)得到的較好的模型，可以進行慢度平方、界面節點位置、反射點位置、散射角和各向異性參數的聯合反演.

由于適用于各向異性介質的運動學反偏移正在研究之中，因此在5.2節的聯合反演測試中，我們并沒有真的使用運動學反偏移來獲得數據空間，我們的數據空間依然是采取從界面出發向上進行射線追蹤到達地表獲得的.也就是說5.2節的反演是直接從第4步開始的，第4步是一個各向同性介質下的三角網格立體層析，具體實現過程可參考邵煒棟等(2017).

5.2 理論數據聯合反演測試

為了驗證5.1節提出的工作流程的有效性，本實驗在如圖8a所示的穹狀速度模型中進行，橫向12 km，縱向6 km，分為7個地質塊體，中間有6個界面.該測試僅使用了180個射線對來生成實驗所需的模型空間.根據5.1節所示工作流程，首先暫時忽略各向異性帶來的影響，利用圖9b所示的Vp0和節點的初始模型，先進行常規二維各向同性三角網格立體層析反演，得到除各向異性參數以外的慢度、節點位置等更為貼近真實模型的初始模型(如圖9c).可以看出相對于圖9b所示的最初始模型，經過幾輪常規各向同性立體層析反演后，圖8c所示模型兩側上翹部分得到了明顯改善，但與真實模型(圖9a)仍有差距.這時使用常規各向同性立體層析反演得到的較好的速度和節點模型(圖9c)，固定模型空間所有其他參數，僅僅進行兩個各向異性參數ε和η的同時反演，就可以得到如圖10c和圖11c所示的ε和η模型的反演中間結果.最后再用圖9c,10c,11c的模型進行速度、節點、各向異性參數的多參數聯合反演，就得到如圖9d、10d、11d所示的最終聯合反演后的速度、節點、各向異性參數模型.從圖12所示的模型殘差上可以看出，最終的反演結果與真實模型已經比較接近.可以證實5.1節提出的工作流程的有效性和之前推導的偏導數的正確性.

圖8 三角網格剖分建立的真實模型Fig.8 True triangulated model

圖9 VP0和節點位置的反演(a) VP0和節點真實模型； (b) VP0和節點初始模型； (c) 控制各向異性參數不變，單獨反演速度和節點； (d) 聯合反演得到最終的VP0和節點.Fig.9 Estimation of VP0 and the interface location(a) Exact model; (b) Very initial model before step (4); (c) Initial model for joint inversion after step (4); (d) Final output model after joint inversion.

圖11 關于η反演結果顯示(a) η真實模型; (b) η初始模型; (c) 用單獨反演的速度和界面信息反演得到的η; (d) 聯合反演得到最終η.Fig.11 Estimation of η(a) Exact model; (b) Very initial model before step (5); (c) Initial model for joint inversion after step (5); (d) Final output model after joint inversion.

圖12 模型殘差顯示(a) 真實VP0-反演最終VP0; (b)真實ε-反演最終ε; (c)真實η-反演最終η.Fig.12 The misfit between the exact model and the final output model(a) Exact VP0-the final output VP0; (b) Exact ε-the final output ε; (c) Exact η-the final output η.

6 討論與結論

慢度平方三角網格剖分模型下實施立體層析具有自然描述界面、同時更新慢度和界面、射線追蹤有解析等優勢.本文在Yang等(2017)、邵煒棟等(2017)實現了慢度平方三角網格二維立體層析反演方法的基礎上，進一步發展了慢度平方三角剖分模型下的VTI介質QP波二維立體層析反演方法.本文所使用的各向異性參數化方式為ε以及η=δ-ε.背景介質定義為橢圓各向異性VTI介質，即η=0的情況.由于在慢度平方橢圓各向異性介質中的射線追蹤有與各向同性介質下形式相近的解析解，使得到非橢圓各向異性介質情形下亦能獲得半解析解，這使得射線追蹤以及FRECHET導數的求取變得簡潔高效.這是本文方法相對于其他方法最為顯著的一個優勢.

在上述理論前提下，我們逐一求取了反演所需的各個數據空間分量對任一模型空間分量的偏導數，從而建立起慢度平方三角網格剖分下二維VTI介質中的QP波各向異性立體層析所需的FRECHET導數矩陣.嚴格的敏感度測試和理論數據的算例一證實了我們推導的FRECHET導數矩陣中的元素具有較好的敏感度，可以用于二維VTI介質QP波各向異性立體層析反演.由于模型空間參數的增加，對所有參數進行同時反演變得相當困難，因此設計了相應的工作流程，首先反演除了各向異性參數之外的其他模型空間分量如慢度、界面信息、反射點位置等，在獲得比較可靠的上述初始模型空間信息基礎上，再反演各向異性參數獲得各類模型參數信息，保證了各向異性參數反演可以得到穩定的結果.理論算例二證實了上述觀點.后續將基于各向異性波動方程正演模擬的理論數據和實際數據進一步發展本方法，使得本文提出的算法真正具有實用價值.

致謝作者感謝自然科學基金面上項目(41874118)對于這項研究的支持.