三維TI介質(zhì)中多波走時(shí)層析成像

黃國嬌 孫江兵 白超英 錢 衛(wèi)

(①河海大學(xué)地球科學(xué)與工程學(xué)院地質(zhì)科學(xué)與工程系，江蘇南京 211100； ②江蘇省地質(zhì)資料館，江蘇南京 210018； ③長安大學(xué)地質(zhì)工程與測繪學(xué)院地球物理系，陜西西安 710054)

1 引言

地震層析成像技術(shù)經(jīng)過近40年的不斷發(fā)展和完善，已成為研究地球內(nèi)部結(jié)構(gòu)以及各向異性的有效技術(shù)之一。隨著油氣等礦產(chǎn)資源的不斷開發(fā)以及對(duì)地殼和上地幔結(jié)構(gòu)的深入研究，對(duì)地震成像分辨率的要求也越來越高。

傳統(tǒng)的各向異性層析成像是在弱各向異性假設(shè)下，主要利用走時(shí)信息進(jìn)行線性反演。如：Hirahara等[1]利用qP波走時(shí)反演三維各向異性參數(shù)；Jech[2]進(jìn)行了三維 TI介質(zhì)走時(shí)反演；Chapman等[3]提出了一種適用于各向異性介質(zhì)中的線性化走時(shí)層析成像方法，并應(yīng)用于井間層析成像。上述研究以弱各向異性假設(shè)為前提，即所有正、反演理論均基于各向同性介質(zhì)背景下的一階擾動(dòng)而建立，簡化了數(shù)學(xué)計(jì)算的難度，且在反演迭代中不用重復(fù)進(jìn)行射線追蹤計(jì)算(雅克比矩陣元素固定不變)，從而有利于進(jìn)行線性反演重建介質(zhì)屬性 。然而，弱各向異性假設(shè)忽略了介質(zhì)橫向變化，反演結(jié)果可信度不高。

近年來，各向異性層析成像已經(jīng)發(fā)展到利用多種地震信息進(jìn)行非線性反演。如前人在遺傳算法、非線性解析法波動(dòng)方程反演等方面進(jìn)行了大量的理論研究[4-6]；Zheng等[7]利用qP波的慢度和質(zhì)點(diǎn)振動(dòng)測量值確定弱各向異性參數(shù)；阮愛國等[8]推導(dǎo)了弱各向異性介質(zhì)中的地震波速，同時(shí)介紹了反演上地幔各向異性參數(shù)的具體算法；郭飚[9]從弱各向異性條件出發(fā)構(gòu)建地震各向異性P波走時(shí)層析成像算法；傅旦丹等[10]采用坐標(biāo)擾動(dòng)法反演一維正交各向異性介質(zhì)的模型參數(shù)；陳竟一等[11]采用模擬退火方法討論了利用VSP觀測的(水平界面)反射波走時(shí)反演橢圓各向異性速度；劉玉柱等[12]實(shí)施了VTI介質(zhì)中的多參數(shù)聯(lián)合走時(shí)反演；劉瑞合等[13]討論了VTI介質(zhì)中各向異性參數(shù)的反演策略；Wang等[14]在介質(zhì)的對(duì)稱軸與界面垂直的假設(shè)條件下，提出了基于網(wǎng)格射線追蹤算法的TTI介質(zhì)走時(shí)層析成像方法；Wang[15]利用修改的牛頓迭代法求解穩(wěn)定非線性問題，采用最小變化約束策略避免每次迭代時(shí)計(jì)算雅可比矩陣元素，并在井間觀測方式下反演了各向異性介質(zhì)模型參數(shù)；Waheed等[16]提出了針對(duì)VTI介質(zhì)的伴隨狀態(tài)法初至波走時(shí)層析成像法，利用梯度法求解非線性問題，避免了直接計(jì)算偏導(dǎo)數(shù)矩陣。然而，縱觀現(xiàn)有的各向異性反演方法，大部分是研究弱各向異性介質(zhì)或簡單的介質(zhì)模型，如一維正交各向異性介質(zhì)模型、橢圓各向異性介質(zhì)及VTI介質(zhì)。

為實(shí)現(xiàn)一般各向異性介質(zhì)中的非線性反演，Zhou等[17-19]討論了走時(shí)關(guān)于各向異性彈性參數(shù)的偏導(dǎo)數(shù)計(jì)算問題，采用初至波(qP、qSV、qSH)走時(shí)分別反演了各向異性參數(shù)。對(duì)于地震層析成像而言分辨率至關(guān)重要，而僅用初至波資料在一定程度上限制了成像深度和空間分辨率。因此，如何提高成像分辨率是一個(gè)亟待解決的問題。黃光南等[20]利用反射qP波走時(shí)實(shí)現(xiàn)了二維TI介質(zhì)層析成像，認(rèn)為利用多波走時(shí)聯(lián)合有助于改善反演結(jié)果。目前，已經(jīng)實(shí)現(xiàn)了二維TI介質(zhì)多震相旅行時(shí)層析成像[21]。

本文研究了三維TI介質(zhì)中多波走時(shí)層析成像技術(shù)。將各向同性介質(zhì)中的分區(qū)多步不規(guī)則最短路徑算法[22，23]推廣至三維TI介質(zhì)中，實(shí)現(xiàn)多次透射、反射以及轉(zhuǎn)換波的追蹤計(jì)算，采用各向同性地震多震相走時(shí)層析成像技術(shù)[24，25]，結(jié)合有關(guān)群速度、相速度偏導(dǎo)數(shù)計(jì)算[17]，通過聯(lián)合多波走時(shí)同時(shí)反演5個(gè)各向異性彈性參數(shù)。數(shù)值模擬結(jié)果表明，多波走時(shí)資料聯(lián)合反演可有效提高成像的空間分辨率。

2 多波射線追蹤計(jì)算方法簡介

將各向同性介質(zhì)中分區(qū)多步不規(guī)則最短路徑算法推廣至各向異性介質(zhì)中進(jìn)行多次波的射線追蹤計(jì)算，其中在單元節(jié)點(diǎn)上采用Zhou等[26]推導(dǎo)的公式計(jì)算群速度。在分區(qū)多步計(jì)算中的分區(qū)是指在模型參數(shù)化時(shí)，可根據(jù)反射界面的具體分布情況將模型分成若干個(gè)獨(dú)立的計(jì)算區(qū)域(相鄰的區(qū)域由界面連接)。多步計(jì)算則是根據(jù)所要計(jì)算地震震相的種類從震源所在的區(qū)域逐區(qū)(步)進(jìn)行計(jì)算。具體來說，自震源所在的單元開始進(jìn)行掃描，按照一定的步長進(jìn)行波前擴(kuò)展，若當(dāng)前區(qū)域所有網(wǎng)格單元內(nèi)節(jié)點(diǎn)都計(jì)算完畢后，波前就停留在該區(qū)的速度離散界面上。若追蹤下行透射波(或透射轉(zhuǎn)換波)，則從速度界面中走時(shí)最小的點(diǎn)開始(根據(jù)惠更斯原理，視為次級(jí)震源)，繼續(xù)新區(qū)內(nèi)的波前擴(kuò)展和射線追蹤；若追蹤上行反射波(或反射轉(zhuǎn)換波)，則從速度界面中走時(shí)最小的點(diǎn)開始，繼續(xù)原區(qū)內(nèi)的波前擴(kuò)展和射線追蹤；若存在轉(zhuǎn)換波，則調(diào)用相應(yīng)的速度模型。重復(fù)上述步驟則可實(shí)現(xiàn)多震相地震射線追蹤計(jì)算。

三維TI介質(zhì)一般可以用7個(gè)物理量m∈{a11,a13,a33,a44,a66,θ0,φ0}(5個(gè)彈性模量參數(shù)和對(duì)稱軸角度)來表述。為演示該算法的射線追蹤能力，給出了一個(gè)兩層TI介質(zhì)中多波地震射線路徑追蹤結(jié)果(圖1)?？梢?，無論是一次反射波(圖1a)還是多次反射轉(zhuǎn)換波(圖1b)，該算法均能進(jìn)行追蹤計(jì)算。

3 反演算法

一般而言，非線性多波(含反射或折射波)走時(shí)聯(lián)合同時(shí)反演，可歸結(jié)為帶約束的阻尼最小二乘最優(yōu)化問題，其反演問題的二范數(shù)目標(biāo)函數(shù)為

(1)

式(1)的法方程為[27]

(2)

式(2)可用循環(huán)迭代的共軛梯度法進(jìn)行求解[28]。目前的關(guān)鍵問題是如何求解具有偏導(dǎo)性質(zhì)的雅克比矩陣A中的元素。

雅克比矩陣元素是走時(shí)關(guān)于彈性參數(shù)的導(dǎo)數(shù)，采用Zhou等[19]推導(dǎo)的一階走時(shí)擾動(dòng)方程求解。該公式不同于傳統(tǒng)的基于特征向量的走時(shí)擾動(dòng)方程[29，30]，不需要計(jì)算特征向量，能夠避免兩個(gè)準(zhǔn)剪切波的奇點(diǎn)問題。

圖1 兩層各向異性介質(zhì)中的多波地震射線 路徑追蹤結(jié)果

(a)一次透射、反射波射線路徑(黑線為P1P1(上角表示上行波，下角表示下行波)，綠線為P1P2P2P1)；(b)多次透射、反射轉(zhuǎn)換波射線路徑(黑線為P1SV1P1，綠線為P1SV1P2SV2P1)

第一層彈性模量參數(shù)分別為：a11=9.09(km/s)2；a13=2.98(km/s)2；a33=7.54(km/s)2；a44=2.28(km/s)2；a66=3.84(km/s)2。對(duì)稱軸角度θ0=0°、φ0=0°。第二層彈性模量參數(shù)分別為：a11=20.31(km/s)2；a13=9.58(km/s)2；a33=22.29(km/s)2；a44=8.35(km/s)2；a66=11.36(km/s)2。對(duì)稱軸角度θ0=45°、φ0=315°

沿著射線路徑R的走時(shí)一階擾動(dòng)方程可以用相速度的導(dǎo)數(shù)表示為

(3)

式中：m為模型參數(shù)，m∈{a11(x),a13(x),a33(x),a44(x),a66(x),θ0(x),φ0(x)}，x=(x1,x2,x3)為射線路徑坐標(biāo)； ds=|dx|。

(4)

根據(jù)式(4)就可以得到雅克比矩陣元素的計(jì)算公式

(5)

值得注意的是，無論是射線路徑和走時(shí)的追蹤計(jì)算，還是反演方程的建立，或是雅克比矩陣的計(jì)算，文中方法都未假設(shè)介質(zhì)為弱各向異性，并且每次迭代都重新追蹤射線路徑和走時(shí)，并計(jì)算雅克比矩陣元素。因此，該方法適用于TI介質(zhì)。

4 數(shù)值模擬實(shí)例及分析

4.1 斷層模型

對(duì)斷層模型(圖3a)模擬井間地震及VSP成像，炮檢排列如圖2a所示。在反演模擬實(shí)驗(yàn)中，采用分區(qū)多步最短路徑法追蹤計(jì)算斷層模型(圖3a)的6類波型(初至波qP、qSV和qSH及模型底界面的一次反射波qPqP、qSVqSV和qSHqSH)的走時(shí)作為反演的觀測走時(shí)。圖2為射線路徑及走時(shí)示意圖。由圖可見，qP速度最快，qSV和qSH波分離，符合各向異性介質(zhì)中波的傳播特點(diǎn)。將斷層上盤介質(zhì)彈性參數(shù)作為初始模型反演斷層的彈性模量參數(shù)。

圖2 射線路徑及走時(shí)示意圖 (a)初至波qP、qSV、qSH射線路徑； (b)反射波qPqP、qSVqSV、qSHqSH射線路徑； (c)初至波qP、qSV、qSH走時(shí)曲線； (d)反射波qPqP、qSVqSV、qSHqSH走時(shí)曲線 將121個(gè)檢波器均勻放置于地表(間距為5m)，40個(gè)檢波器、40個(gè)炮點(diǎn)均置于4口鉆井中(模型直立的四條棱，每個(gè)鉆孔 放置10個(gè)檢波器，間隔為5m)，且炮點(diǎn)與檢波器位置相同，當(dāng)其中某口井作為激發(fā)井時(shí)，則另外3口井作為接收井

TI介質(zhì)的5個(gè)彈性參數(shù)被分成兩部分，一是與qP、qSV波相關(guān)的4個(gè)彈性參數(shù)(a11、a13、a33和a44)，二是與qSH波相關(guān)的2個(gè)彈性參數(shù)(a44和a66)。當(dāng)采用某一種波反演時(shí)，只能得到與之對(duì)應(yīng)的彈性參數(shù)。圖3為初至qP波反演彈性參數(shù)的剖面結(jié)果。從整體上分析，y=0剖面反演結(jié)果(圖3b)明顯好于y=15m剖面(圖3c)，這主要是因?yàn)閥=0剖面上既有震源、又有接收點(diǎn)，射線平面即為y=0剖面，該剖面上彈性參數(shù)的反演退化為二維反演，射線覆蓋相對(duì)較好。其中a44的反演效果最好，而其他三個(gè)彈性參數(shù)(a11、a13、a33)雖然恢復(fù)了斷層形態(tài)和位置，但是未完全反演出斷層下盤介質(zhì)的彈性參數(shù)，說明qP波對(duì)a44最敏感。

當(dāng)聯(lián)合初至波qP、qSV和qSH走時(shí)反演彈性參數(shù)時(shí)，能夠反演5個(gè)彈性參數(shù)(圖4)。對(duì)比初至qP波反演彈性參數(shù)的剖面結(jié)果(圖3)與初至波qP、qSV和qSH聯(lián)合反演彈性參數(shù)的剖面結(jié)果(圖4)發(fā)現(xiàn)：后者的y=0剖面的5個(gè)彈性參數(shù)反演結(jié)果(圖4b)幾乎與真實(shí)模型(圖4a)完全一致，無論斷層的起伏、形態(tài)和位置，還是斷層上、下盤介質(zhì)的彈性參數(shù)值均得到了很好的重建；后者的y=15m剖面的5個(gè)彈性參數(shù)的反演效果不是很理想(圖4c)，但較前者的y=15m剖面(圖3c)有所改善，依稀可見斷層的形態(tài)和位置。

圖5為初至波qP、qSV、qSH和反射波qPqP、qSVqSV、qSHqSH聯(lián)合反演彈性參數(shù)的剖面結(jié)果。由圖可見，當(dāng)加入反射波后，成像分辨率顯著提高，特別是分辨率較差的剖面(圖4c)的斷層形態(tài)和位置也得到很好的重建(圖5c)。這主要是因?yàn)榉瓷洳╭PqP、qSVqSV、qSHqSH的射線路徑(圖2b)與初至波qP、qSV和qSH的射線路徑(圖2a)存在較大區(qū)別，當(dāng)加入反射波之后，射線對(duì)模型的覆蓋率及角度覆蓋明顯增加，進(jìn)而成像分辨率顯著提高。因此，聯(lián)合多波提高走時(shí)成像分辨率，實(shí)際上是通過增加有效的射線覆蓋和角度覆蓋實(shí)現(xiàn)的。

4.2 高、低速異常體模型

圖6為真實(shí)各向異性介質(zhì)模型，各向異性彈性模量參數(shù)及相應(yīng)的Thomsen參數(shù)如表1所示。模擬井間地震及VSP成像時(shí)的觀測系統(tǒng)同斷層模型(圖2a)。實(shí)驗(yàn)中采用6類波型的走時(shí)信息，即初至波qP、qSV和qSH與模型底界面的一次反射波qPqP、 qSVqSV和qSHqSH。將背景介質(zhì)彈性參數(shù)作為初始模型反演高、低速異常體的彈性模量參數(shù)。

圖3 初至波qP反演彈性參數(shù)的剖面結(jié)果 (a)真實(shí)剖面；(b)y=0剖面；(c)y=15m剖面

斷層模型尺度為50m×50m×50m，斷層上盤介質(zhì)的彈性模量參數(shù)為a11=9.09(km/s)2、a13=2.981(km/s)2、a33=7.535(km/s)2、a44=2.274(km/s)2、a66=3.843(km/s)2,對(duì)稱軸角度θ0=45°,φ0=0°。斷層下盤介質(zhì)的彈性模量參數(shù)為a11=13.84(km/s)2、a13=4.245(km/s)2、a33=11.34(km/s)2、a44=3.345(km/s)2、a66=5.051(km/s)2，對(duì)稱軸角度為θ0=45°、φ0=0°

表1 高、低速異常體彈性參數(shù)及相應(yīng)的Thomsen參數(shù)表

圖7為初至波qP彈性參數(shù)反演結(jié)果。由圖可見，反演的彈性參數(shù)能大致反映高、低速異常體的位置和形態(tài)，且高速異常體對(duì)應(yīng)的彈性參數(shù)的反演效果明顯好于低速異常體，這主要是因?yàn)楦咚佼惓ｓw對(duì)地震射線具有匯集作用； 在反演的4個(gè)彈性參數(shù)中，a44的反演分辨率最高，在10(圖7a)、15m切片(圖7b)和15m剖面(圖7d)上的高速異常體的形態(tài)和幅度恢復(fù)得較好，這說明a44對(duì)qP波群(相)速度起決定作用； 隨著深度的增加，成像分辨率越來越低，如20m切片(圖7c)只能反演異常位置，不能較好地恢復(fù)異常形狀和幅度，不同位置的剖面圖(圖7d、圖7e)也證實(shí)了這種現(xiàn)象，即無論是高速異常體還是低速異常體，均未較好地重建下部的異常輪廓。

圖4 初至波qP、qSV和qSH聯(lián)合反演彈性參數(shù)的剖面結(jié)果 (a)真實(shí)剖面； (b)y=0剖面； (c)y=15m剖面

為了有效提高走時(shí)成像分辨率，進(jìn)行了初至波qP、qSV和qSH聯(lián)合反演成像(圖8)及初至波qP、qSV和qSH與反射波qPqP、qSVqSV、qSH qSH聯(lián)合反演成像(圖9)。分析初至波qP反演彈性參數(shù)的結(jié)果(圖7)和初至波qP、qSV和qSH聯(lián)合反演彈性參數(shù)的結(jié)果(圖8)發(fā)現(xiàn)，后者不僅能反演5個(gè)彈性參數(shù)(a11、a13、a33、a44和a66)，而且反演分辨率明顯提高，很好地重建了高、低速異常的形態(tài)、幅度及深度(只是20m深度切片的結(jié)果差些)，特別是也能大致看出前者中分辨率較差的a33(圖7)的異常形態(tài)和位置。

分析初至波qP、qSV、qSH和反射波qPqP、qSVqSV、qSH qSH聯(lián)合反演彈性參數(shù)的結(jié)果(圖9)可知，反演的5個(gè)彈性參數(shù)均很好地重建了高、低速異常體，所恢復(fù)的高、低速異常體的形狀、位置、幅度接近于真實(shí)模型，異常邊界也清晰可見。因此，聯(lián)合多波走時(shí)資料可以有效提高成像分辨率。

模型尺度為50m×50m×50m，模型中包含兩個(gè)速度異常體(紅色為高速異常體，藍(lán)色為低速異常體)，大小均為10m×10m×10m，深度范圍為10～20m，背景和異常體均為TI介質(zhì)，對(duì)稱軸角度都為(45°,0°)

圖7 初至波qP彈性參數(shù)反演結(jié)果 (a)z=10m切片； (b)z=15m切片； (c)z=20m切片； (d)y=15m剖面； (e)y=35m剖面

圖8 初至波qP、qSV和qSH彈性參數(shù)聯(lián)合反演結(jié)果 (a)z=10m切片； (b)z=15m切片； (c)z=20m切片； (d)y=15m剖面； (e)y=35m剖面

圖9 初至波qP、qSV、qSH與反射波qPqP、qSVqSV、qSHqSH彈性參數(shù)聯(lián)合反演結(jié)果 (a)z=10m切片； (b)z=15m切片； (c)z=20m切片； (d)y=15m剖面； (e)y=35m剖面

5 結(jié)束語

地球內(nèi)部地震各向異性介質(zhì)廣泛存在，前人在研究地球內(nèi)部地震各向異性時(shí)都假設(shè)介質(zhì)為弱各向異性，且多為VTI或HTI介質(zhì)，從而簡化了算法，并根據(jù)擾動(dòng)理論更新模型參數(shù)，即不需要每次迭代都進(jìn)行射線追蹤及雅克比矩陣元素的計(jì)算。但是這種處理忽略了介質(zhì)的橫向不均勻性，結(jié)果往往很難與實(shí)際情況吻合。為了進(jìn)行TI介質(zhì)的多波走時(shí)層析成像，本文利用一階走時(shí)擾動(dòng)方程計(jì)算雅克比矩陣元素，并結(jié)合分區(qū)多步不規(guī)則最短路徑射線追蹤算法和共軛梯度法求解帶約束的阻尼最小二乘問題，形成了一套可用于TI介質(zhì)中多波走時(shí)聯(lián)合反演各向異性彈性參數(shù)的技術(shù)。數(shù)值模擬實(shí)驗(yàn)證明該方法不僅能有效地反演各向異性彈性參數(shù)，而且還能提高成像分辨率，是一種實(shí)用的反演算法。

[1] Hirahara K and Yuzo I.Traveltime inversion for three-dimensional P-wave velocity anisotropy.Journal of Physics of the Earth,1984,32(3):197-218.

[2] Jech J.Three-dimensional inversion problem for inhomogeneous transversely isotropic media.Studia Geophysica et Geodetica,1988,32(2):136-143.

[3] Chapman C H and Pratt G R.Traveltime tomography in anisotropic media—Ⅰ.Theory.Geophysical Journal International,1992,109(1):1-19.

[4] 周輝,何樵登.非線性各向異性波形反演方法.石油地球物理勘探,1995,30(6):725-735. Zhou Hui and He Qiaodeng.Nonlinear waveform inversion in anisotropic medium.OGP,1995,30(6):725-735.

[5] 張文生,何樵登.利用VSP和跨孔走時(shí)反演橫向各向同性介質(zhì)中的彈性常數(shù).石油地球物理勘探,1997,32(3):345-356. Zhang Wensheng,He Qiaodeng.Deriving the elastic constants of transversal isotropic medium from VSP and crosshole travel times.OGP,1997,32(3):345-356.

[6] 張秉銘,張中杰.一種新的地層彈性參數(shù)直接反演方法.地震學(xué)報(bào),2000,22(6):654-660. Zhang Bingming,Zhang Zhongjie.A new way of directly elastic parameters inversion.Acta Seismologica Sinica,2000,22(6):654-660.

[7] Zheng X Y and Psenick I.Local determination of weak anisotropy parameters from qP-wave slowness and particle motion measurements.Pure and Applied Geophysics,2002,159(7-8):1881-1905.

[8] 阮愛國,王椿鏞.上地幔各向異性的反演方法.西北地震學(xué)報(bào),2002,24(2):104-112. Ruan Aiguo,Wang Chunyong.The inversion methods of the upper mantle anisotropy.Northwestern Seismological Journal,2002,24(2):104-112.

[9] 郭飚.非均勻各向異性介質(zhì)的地震P波走時(shí)層析成像研究[學(xué)位論文].北京:中國地震局地震研究所,2010.

[10] 傅旦丹,何樵登,劉一峰等.一種新的全局最優(yōu)化反演計(jì)算方法.石油地球物理勘探,2000,35(4):536-542. Fu Dandan,He Qiaodeng,Liu Yifeng et al.A novel global-optimized inversion method.OGP,2000,35(4):536-542.

[11] 陳竟一,張中杰.利用有偏VSP反射波旅行時(shí)反演橢圓各向異性速度—模擬退火方法.中國地球物理學(xué)會(huì)第二十屆年會(huì)論文集,2004,501.

[12] 劉玉柱,王光銀,董良國等.VTI介質(zhì)多參數(shù)聯(lián)合走時(shí)層析成像方法.地球物理學(xué)報(bào),2014,57(10):3402-3410. Liu Yuzhu,Wang Guangyin,Dong Liangguo et al.Joint inversion of VTI parameters using nonlinear traveltime tomography.Chinese Journal of Geophy-sics,2014,57(10):3402-3410.

[13] 劉瑞合,趙金玉,印興耀等.VTI介質(zhì)各向異性參數(shù)層析反演策略與應(yīng)用.石油地球物理勘探,2017,52(3):484-490. Liu Ruihe,Zhao Jinyu,Yin Xingyao et al.Strategy of anisotropic parameter tomography inversion in VTI medium.OGP,2017,52(3):484-490.

[14] Wang X X,Tsvankin I.Ray-based gridded tomography for tilted transversely isotropic media.Geophysics,2013,78(1):C11-C23.

[15] Wang Y H.Seismic ray tracing in anisotropic media:A modified Newton algorithm for solving highly nonlinear systems.Geophysics,2014,79(1):T1-T7.

[16] Waheed B U,Flagg G and Yarman C E.First-arrival traveltime tomography for anisotropic media using the adjoint-state method.Geophysics,2016,81(4):R147-R155.

[17] Zhou B and Greenhalgh S A.Analytic expressions for the velocity sensitivity to the elastic moduli for the most general anisotropic media.Geophysical Prospecting,2005,53(4):619-641.

[18] Zhou B and Greenhalgh S A.Nonlinear traveltime inversion scheme for crosshole seismic tomography in titled transversly isotropic media.Geophysics,2008,73(4):D17-D33.

[19] Zhou B and Greenhalgh S A.Nonlinear traveltime inversion for 3D seismic tomography in strongly anisotropic media.Geophysical Journal International,2008,172(1):383-394.

[20] 黃光南,Zhou Bing,鄧居智等.各向異性TI介質(zhì)qP反射波走時(shí)層析成像.地球物理學(xué)報(bào),2015,58(6):2035-2045. Huang Guangnan,Zhou Bing,Deng Juzhi et al.Tra-veltime tomography of qP reflection waves in anisotropic TI media.Chinese Journal of Geophysics,2015,58(6):2035-2045.

[21] 黃國嬌,白超英,錢衛(wèi).一般TI介質(zhì)中多震相旅行時(shí)層析成像——以井間成像為例.石油地球物理勘探,2016,51(1):115-126. Huang Guojiao,Bai Chaoying,Qian Wei.Multi-phase traveltime tomography in general anisotropic TI media:an example of crosshole tomography.OGP,2016,51(1)：115-126.

[22] Bai C Y,Tang X P and Zhao R.2D/3D multiple transmitted,converted and reflected arrivals in complex layered media with the modified shortest path method.Geophysical Journal International,2009,179(1):201-214.

[23] Bai C Y,Huang G J and Zhao R.2D/3D irregular shortest-path raytracing for multiple arrivals and its applications.Geophysical Journal International,2010,183(3):1596-1612.

[24] 黃國嬌,白超英.二維復(fù)雜層狀介質(zhì)中地震多波走時(shí)聯(lián)合反演成像.地球物理學(xué)報(bào),2010,53(12):2972-2981. Huang Guojiao,Bai Chaoying.Simultaneous inversion with multiple traveltimes within 2-D complex layered media.Chinese Journal of Geophysics,2010,53(12):2972-2981.

[25] 白超英,黃國嬌,李忠生.三維復(fù)雜層狀介質(zhì)中多震相走時(shí)聯(lián)合反演成像.地球物理學(xué)報(bào),2011,54(1):182-192. Bai Chaoying,Huang Guojiao,Li Zhongsheng.Simultaneous inversion combining multiple-phase travel times within 3D complex layered media.Chinese Journal of Geophysics,2011,54(1):182-192.

[26] Zhou B and Greenhalgh S A.Raypath and traveltime computations for 2D transversely isotropic media with dipping symmetry axes.Exploration Geophysics,2006,37(2):150-159.

[27] Menke W.Geophysical Data Analysis:Discrete Inverse Theory.Academic Press Inc,San Diego,California,1984,40-126.

[28] Zhou B,Greenhalgh S A,Sinadinovski C.Iterative algorithm for the damped minimum norm,least-squares and constrained problem in seismic tomography.Exploration Geophysics,1992,23(3):497-505.

[29] Crampin S.Effective anisotropic elastic constants for wave-propagation through cracked solids.Geophysical Journal of the Royal Astronomical Society,1984,76(1):135-145.