周嘉欣， 楊鍇， 邵煒棟

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

0 引言

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

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

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

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

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

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

1 各向異性射線理論簡介

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

1.1 Christoffel方程

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

(1)

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

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

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

(2)

其中

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

(3)

(4)

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

(5a)

(5b)

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

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

根據(jù)Thomsen參數(shù)表征有：

(6)

(7)

結(jié)合(2)式，可以得到如下關(guān)系：

(8)

+4(δ-ε)A33.

(9)

(10)

其中

(11)

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

(12)

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

(13)

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

(14)

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

(15)

(16a)

(16b)

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

+x(τ0),

(17)

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

(18)

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

(19a)

(19b)

(19c)

(19d)

2 慢度三角網(wǎng)格模型下VTI介質(zhì)二維QP波立體層析的數(shù)據(jù)空間與模型空間簡介

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

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

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

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

Δd=FΔm.

(20)

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

(21)

在這個矩陣中，顯然有：

(22)

其次，根據(jù)費馬原理，走時T在速度結(jié)構(gòu)、激發(fā)點與接收點位置都固定時，已經(jīng)是最小走時.因此可以認(rèn)為散射角對走時T不存在一階擾動影響.即可認(rèn)為

(23)

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

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

(24)

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

3 三角網(wǎng)格慢度平方模型下VTI介質(zhì)二維QP波各向異性立體層析FRECHET導(dǎo)數(shù)求取

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

3.1 基于射線擾動理論推導(dǎo)數(shù)據(jù)空間關(guān)于ε和η的一階擾動關(guān)系

以下基于射線擾動理論推導(dǎo)數(shù)據(jù)空間各分量關(guān)于各向異性參數(shù)ε和η的一階擾動公式.我們首先定義一個新的Hamiltonian表達(dá)式：

(25)

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

(26)

(27)

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

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

(28)

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

(29)

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

(30)

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

(31)

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

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

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

Δpx=Δpx0,

Δpz=Δpz0.

(32a)

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

(32b)

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

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

Δ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)

因此很容易寫出地面觀測時，坐標(biāo)X、射線參數(shù)Px關(guān)于各向異性參數(shù)η的Frechet導(dǎo)數(shù)為

(34a)

(34b)

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

關(guān)于走時T、η的一階擾動關(guān)系，F(xiàn)arra等(1989)定義了如下走時擾動公式：

(35)

因此很容易寫出T關(guān)于η的FRECHET導(dǎo)數(shù)為

(36)

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

至此當(dāng)射線在一個地質(zhì)塊內(nèi)傳播時，立體層析數(shù)據(jù)空間關(guān)于ε和η的偏導(dǎo)數(shù)已經(jīng)全部得到.

3.2 非橢圓各向異性VTI介質(zhì)中的半解析射線追蹤及上述一階擾動關(guān)系的數(shù)值驗證

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

(37)

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

(38)

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

圖2 黑色曲線展示了橢圓各向異性介質(zhì)下的射線追蹤，灰色曲線為在橢圓各向異性上加以δ的擾動的小擾動計算得到的射線.其中=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)式計算的射線路徑，而綠色實線則是表示常規(guī)龍格庫塔方法解得的射線追蹤.測試參數(shù) = 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導(dǎo)數(shù)推導(dǎo)及其敏感度測試

4.1 射線穿越多個界面時FRECHET導(dǎo)數(shù)推導(dǎo)

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

(39)

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

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

(40a)

(40b)

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

4.2 與各向異性參數(shù)ε和η有關(guān)FRECHET導(dǎo)數(shù)敏感度測試

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

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

圖4 敏感度測試模型以及三角網(wǎng)格立體層析中Sx, T關(guān)于ε和η的敏感度測試(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/η.

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

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

(41)

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

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

5 三角網(wǎng)格立體層析理論數(shù)據(jù)反演測試

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

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

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

圖6 三層模型各向異性參數(shù)反演測試(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.

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

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

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

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

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

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

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

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

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

5.2 理論數(shù)據(jù)聯(lián)合反演測試

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

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

圖9 VP0和節(jié)點位置的反演(a) VP0和節(jié)點真實模型； (b) VP0和節(jié)點初始模型； (c) 控制各向異性參數(shù)不變，單獨反演速度和節(jié)點； (d) 聯(lián)合反演得到最終的VP0和節(jié)點.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 關(guān)于η反演結(jié)果顯示(a) η真實模型; (b) η初始模型; (c) 用單獨反演的速度和界面信息反演得到的η; (d) 聯(lián)合反演得到最終η.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 討論與結(jié)論

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

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

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