黃光南， ZHOU Bing， 鄧居智， 李紅星， 李澤林， 張華

1 東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室， 南昌 330013 2 中國地質(zhì)大學(xué)(武漢)地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 430074 3 中國石油大學(xué)(北京)油氣資源與探測國家重點(diǎn)實(shí)驗(yàn)室， 北京 102249 4 Petroleum Geosciences, the Petroleum Institute, P.O. Box 2533, Abu Dhabi, UAE

?

各向異性TI介質(zhì)qP反射波走時(shí)層析成像

黃光南1,2,3， ZHOU Bing4， 鄧居智1， 李紅星1， 李澤林1， 張華1

1 東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室， 南昌 330013 2 中國地質(zhì)大學(xué)(武漢)地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 430074 3 中國石油大學(xué)(北京)油氣資源與探測國家重點(diǎn)實(shí)驗(yàn)室， 北京 102249 4 Petroleum Geosciences, the Petroleum Institute, P.O. Box 2533, Abu Dhabi, UAE

地震走時(shí)層析成像是反演地層各向異性參數(shù)分布的有效方法，但是關(guān)于地震各向異性介質(zhì)走時(shí)層析成像的研究并不多，其技術(shù)遠(yuǎn)遠(yuǎn)沒有達(dá)到成熟的階段.在野外數(shù)據(jù)采集時(shí)，地表反射波觀測方式相對井間和垂直地震剖面觀測方式的成本更低，利用qP反射波走時(shí)反演各向異性參數(shù)具有更加廣泛的實(shí)用價(jià)值.本文實(shí)現(xiàn)的TI介質(zhì)地震走時(shí)層析成像方法結(jié)合了TI介質(zhì)反射波射線追蹤算法、走時(shí)擾動(dòng)方程和非線性共軛梯度算法，它可以對任意強(qiáng)度的TI介質(zhì)模型進(jìn)行反演，文中嘗試?yán)胵P反射波走時(shí)重建TI介質(zhì)模型的參數(shù)圖像.利用qP反射波對層狀介質(zhì)模型和塊狀異常體模型進(jìn)行走時(shí)反演，由于qP波相速度對彈性模量參數(shù)和Thomsen參數(shù)的偏微分不同，所以可以分別反演彈性模量參數(shù)和Thomsen參數(shù).數(shù)值模擬結(jié)果表明：利用qP反射波可以反演出TI介質(zhì)模型的彈性模量參數(shù)與Thomsen參數(shù)，不同模型的走時(shí)迭代反演達(dá)到了較好的收斂效果，與各向同性介質(zhì)走時(shí)反演結(jié)果相比較，各向異性介質(zhì)走時(shí)反演結(jié)果具有較好的識別能力.

各向異性介質(zhì)；各向異性參數(shù)；彈性模量參數(shù)；走時(shí)層析成像

1 引言

室內(nèi)研究與野外觀測表明大多數(shù)巖石和地層存在明顯的各向異性.引起地層各向異性的主要因素有：巖石的結(jié)構(gòu)各向異性，地層方向應(yīng)力引起的各向異性，巖性各向異性，巖石晶體定向排列引起的各向異性，以及巖石裂隙裂縫引起的各向異性(Crampin, 1984; Helbig, 1981; 李芳等, 2012).地震波在各向同性介質(zhì)與各向異性介質(zhì)中的傳播特征不同，在各向異性介質(zhì)有三種體波：qP波，qSV波和qSH波，這三種波各自以自身的傳播速度與相互垂直的極化方向傳播(Zhou and Greenhalgh, 2006).因此，如果地層具有各向異性，那么地震走時(shí)層析成像方法必須適用于這種介質(zhì)模型，否則，地震走時(shí)反演結(jié)果可能會(huì)產(chǎn)生錯(cuò)誤的地質(zhì)解釋.

經(jīng)過40多年的發(fā)展歷程，地震各向同性介質(zhì)走時(shí)層析成像方法已經(jīng)比較成熟，很多學(xué)者利用各種射線追蹤算法、數(shù)學(xué)反演算法和走時(shí)擾動(dòng)方程得到了不同的地震走時(shí)層析成像方法.這里列舉一些在勘探地球物理學(xué)方面有代表性的地震走時(shí)層析成像方法.Zelt等(2006)的有限差分初至波走時(shí)層析成像方法得到了廣泛的應(yīng)用，利用這種初至波走時(shí)層析成像方法在地下水污染的圈定和調(diào)查方面取得了良好的應(yīng)用效果.Tryggvason等(2009)結(jié)合初至波走時(shí)層析和靜校正方法，將其運(yùn)用于核廢料處置的勘探選址研究，所得地震剖面的可靠性與分辨率均得到了較大提高.Zhou(2006)提出了地震多尺度可變形地層走時(shí)層析成像方法，數(shù)值測試結(jié)果表明該方法比基于晶格或者網(wǎng)格節(jié)點(diǎn)的走時(shí)層析成像方法具有較高的成像精度.劉玉柱和楊積忠(2014)利用先驗(yàn)信息在初至波走時(shí)層析成像過程中加入正則化約束條件，反演結(jié)果使原來隱藏的特征變得更加清晰.

關(guān)于地震各向異性介質(zhì)走時(shí)層析成像的研究也在不斷升溫，周輝和何樵登(1995)利用層狀TI介質(zhì)的P波和SH波走時(shí)，結(jié)合雙曲線擬合法和遺傳算法反演得到了五個(gè)彈性模量參數(shù)值.張文生和何樵登(1999)利用跨孔qP波和qSV波進(jìn)行了統(tǒng)計(jì)法走時(shí)層析成像，提高了計(jì)算結(jié)果的穩(wěn)定性和精度.盧明輝等(2005)利用正交各向異性介質(zhì)P波走時(shí)反演得到了單層介質(zhì)的Thomsen參數(shù).劉玉柱等(2014)利用多參數(shù)聯(lián)合走時(shí)層析成像方法對VTI介質(zhì)進(jìn)行走時(shí)反演取得了較好的數(shù)值測試結(jié)果.國外對地震各向異性介質(zhì)走時(shí)層析成像的研究相對較多，Cerveny(2001)是較早從事各向異性走時(shí)層析成像的學(xué)者，他運(yùn)用各向異性介質(zhì)“程函方程”的哈密爾頓形式，推導(dǎo)出了各向異性介質(zhì)走時(shí)擾動(dòng)方程.后續(xù)學(xué)者利用該走時(shí)擾動(dòng)方程和不同的觀測系統(tǒng)對各種各向異性介質(zhì)模型進(jìn)行反演取得了不同的效果.Jech(1988)對該走時(shí)擾動(dòng)方程進(jìn)行了修改，并且實(shí)現(xiàn)了三維各向異性介質(zhì)走時(shí)層析成像方法.Chapman和Pratt(1992)、Pratt和Chapman(1992)提出了各向異性介質(zhì)走時(shí)反演的線性化方法，利用井間觀測方式進(jìn)行了走時(shí)反演.Wang 和Tsvankin(2013)提出了網(wǎng)格化的各向異性介質(zhì)走時(shí)層析成像方法，該方法極大地提高了成像的精度.Wang(2014)對各向異性介質(zhì)彎曲射線法產(chǎn)生的高度非線性問題，提出了修改的牛頓迭代法用于穩(wěn)定非線性問題的解，并利用井間觀測方式反演了各向異性介質(zhì)模型參數(shù).縱觀大多數(shù)各向異性介質(zhì)走時(shí)層析成像方法，它們都是采用了弱各向異性介質(zhì)假設(shè)條件.實(shí)際上，巖石和地層的各向異性強(qiáng)度是不確定的，將它們假設(shè)為弱各向異性介質(zhì)具有不合理性.Zhou等(2008)結(jié)合各向異性介質(zhì)初至波射線追蹤算法、走時(shí)擾動(dòng)方程和非線性共軛梯度，實(shí)現(xiàn)了地震各向異性初至波走時(shí)層析成像方法.這種地震走時(shí)層析成像方法的優(yōu)點(diǎn)是能夠適用于任意強(qiáng)度的各向異性介質(zhì)模型，它的各向異性介質(zhì)走時(shí)擾動(dòng)方程避免了qSV波和qSH波的奇異值問題，非線性共軛梯度算法可以同時(shí)反演出幾個(gè)彈性模量參數(shù)(或者Thomsen參數(shù)).

國內(nèi)外關(guān)于地震各向異性介質(zhì)走時(shí)層析成像的文章不是很多，其技術(shù)遠(yuǎn)遠(yuǎn)沒有達(dá)到成熟的階段，因此仍然有很多的理論和數(shù)值模擬研究值得探索.本文以Zhou等(2008)的各向異性初至波走時(shí)層析成像方法為基礎(chǔ)，利用反射波射線追蹤算法實(shí)現(xiàn)了TI介質(zhì)qP反射波非線性走時(shí)層析成像.由于地表反射波觀測方式相比井間和垂直地震剖面觀測方式的成本更低，因此TI介質(zhì)qP反射波非線性走時(shí)層析成像方法具有更加廣泛的實(shí)用價(jià)值.在數(shù)值模擬部分，利用qP反射波走時(shí)重建層狀TI介質(zhì)模型和塊狀異常體TI介質(zhì)模型的圖像，用于測試各向異性介質(zhì)參數(shù)重建和成像的可能性.通過分析各種模型的彈性模量參數(shù)與Thomsen參數(shù)反演結(jié)果，驗(yàn)證TI介質(zhì)qP反射波走時(shí)層析成像方法的正確性.

2 方法原理

2.1 各向異性介質(zhì)走時(shí)擾動(dòng)方程

Cerveny推導(dǎo)出了各向異性介質(zhì)走時(shí)擾動(dòng)方程的表達(dá)式.很多學(xué)者對該公式做了不同的修改，Zhou推導(dǎo)出了各向異性介質(zhì)走時(shí)擾動(dòng)方程的三種表達(dá)形式，這三種表達(dá)式被稱為一階走時(shí)擾動(dòng)方程的波速偏微分形式，其中最為簡單的是走時(shí)擾動(dòng)方程的相速度偏微分形式(Zhou et al., 2008)

(1)

(2)

各向異性模型參數(shù)化后，可以得到相速度vp(m,n)與群速度vg(m,n)模型.

(3)

2.2 雅可比矩陣表達(dá)式

當(dāng)震源與檢波器的數(shù)量比較大時(shí)，它們之間存在大量的射線路徑，公式(3)可以寫成非線性反演的矩陣形式(Zhou et al., 2008)

(l=1,2,…)

(4)

式中上標(biāo)-g表示廣義逆矩陣，Wd和Wm分別代表數(shù)據(jù)與模型參數(shù)的權(quán)重矩陣，λ是平衡數(shù)據(jù)擬合度與模型粗糙度之間的變量.(4)式可以利用非線性共軛梯度算法進(jìn)行求解(Zhouetal., 1992)，其他變量的表達(dá)式如下:

(5)

(6)

3 數(shù)值模擬

下面利用各向異性介質(zhì)qP反射波走時(shí)層析成像進(jìn)行數(shù)值模擬.數(shù)值模擬包括兩部分：第一部分是利用qP反射波走時(shí)反演TI介質(zhì)模型的彈性模量參數(shù)；第二部分是利用qP反射波走時(shí)反演TI介質(zhì)模型的Thomsen參數(shù).測試模型包括層狀介質(zhì)模型和塊狀異常體模型.Zhou對各向異性介質(zhì)初至波走時(shí)層析成像做過相應(yīng)的研究工作，并且指出當(dāng)?shù)貙訉ΨQ軸傾角不固定時(shí)，算法很難反演得到合理的傾角參數(shù)分布.原因是相速度和正演走時(shí)均對地層對稱軸傾角非常敏感，因此本文假設(shè)模型的地層對稱軸傾角為已知，只反演模型的彈性模量參數(shù)或者Thomsen參數(shù)分布.數(shù)值模擬實(shí)驗(yàn)的觀測走時(shí)是射線追蹤得到的正演走時(shí)，并且添加了一定百分比的高斯噪聲.在非線性走時(shí)反演時(shí)，Wd和Wm是根據(jù)數(shù)據(jù)和模型先驗(yàn)信息確定的權(quán)重矩陣，用以減小反演的多解性.文中假設(shè)先驗(yàn)信息的權(quán)重相同，即數(shù)據(jù)的權(quán)重矩陣和模型的權(quán)重矩陣均為單位矩陣(Wd=Id和Wm=Im).阻尼因子λ是通過多次數(shù)值測試，比較反演的收斂情況而選取的數(shù)值大小.在反演迭代過程中，利用了阿爾法中值濾波器剔除每次迭代結(jié)果的高波數(shù)噪聲.地震走時(shí)層析成像的反演算法是非線性共軛梯度，Zhou等(1992)將其與不同反演算法做過數(shù)值和實(shí)際數(shù)據(jù)試驗(yàn)，證明了它是一種有效而快速的反演算法.

3.1 qP反射波走時(shí)反演彈性模量參數(shù)模型

下面分別對層狀介質(zhì)模型和塊狀異常體模型進(jìn)行TI介質(zhì)反射波非線性走時(shí)反演，求取模型的彈性模量參數(shù).

(1)層狀介質(zhì)模型的彈性模量參數(shù)反演.層狀介質(zhì)模型的長度為800 m，深度為600 m，第一層介質(zhì)的彈性模量參數(shù)為c11=8.81 GPa,c13=2.78 GPa,c33=7.34 GPa,c44=2.28 GPa,c66=3.88 GPa，對稱軸傾角為θ0=0°.第二層介質(zhì)的彈性模量參數(shù)為c11=21.49 GPa,c13=10.85 GPa,c33=19.54 GPa,c44=8.71 GPa,c66=8.45 GPa，對稱軸傾角為θ0=90°.第三層介質(zhì)的彈性模量參數(shù)為c11=13.86 GPa,c13=6.23 GPa,c33=12.36 GPa,c44=4.35 GPa,c66=5.09 GPa，對稱軸傾角為θ0=0°.采用地表激發(fā)和接收的觀測方式，炮點(diǎn)與檢波點(diǎn)位置相同，道間距為25 m，炮點(diǎn)與檢波點(diǎn)均為33個(gè)，模型的網(wǎng)格間距為50 m.將第一層介質(zhì)參數(shù)作為初始模型用以反演第二層和第三層介質(zhì)的彈性模量參數(shù).為了改善反演成像的效果，我們對模型設(shè)置了三個(gè)反射界面，充分利用三個(gè)反射界面的反射波走時(shí)信息.圖1a為層狀彈性模量參數(shù)模型的反射波射線路徑分布圖，射線覆蓋了整個(gè)模型范圍，圖1b是利用相同的各向異性介質(zhì)正演走時(shí)結(jié)合縱波速度走時(shí)層析成像所得速度反演結(jié)果，可以看出速度反演結(jié)果效果不理想，第一層速度值變化較大，第二層下界面位置不正確.圖2為層狀介質(zhì)模型的彈性模量參數(shù)(a—c—e—g)與反演結(jié)果(b—d—f—h)的對比圖，左側(cè)為真實(shí)彈性模量參數(shù)，右側(cè)為反演結(jié)果，從整體上分析，c11,c13,c33和c44均得到了較好的反演效果.由迭代反演收斂曲線(圖9)可知反演得到了較好的收斂水平，表明反演結(jié)果與真實(shí)模型相接近.

(2)塊狀異常體模型的彈性模量參數(shù)反演.塊狀異常體模型的長度為800 m，深度為500 m，背景介質(zhì)的彈性模量參數(shù)為c11=12.0 GPa,c13=5.0 GPa,c33=10.0 GPa,c44=4.0 GPa,c66=6.0 GPa，對稱軸傾角為θ0=0°.塊狀異常體的彈性模量參數(shù)為c11=20.0 GPa,c13=10.0 GPa,c33=15.0 GPa,c44=8.0 GPa,c66=12.0 GPa，對稱軸傾角為θ0=45°.采用地表激發(fā)和接收的觀測方式，炮點(diǎn)與檢波點(diǎn)位置重合，水平道間隔為16 m，炮點(diǎn)與檢波點(diǎn)均為51個(gè)，模型的網(wǎng)格間距為50 m.將背景介質(zhì)參數(shù)作為初始模型用以反演塊狀異常體的彈性模量參數(shù).為了改善反演成像的效果，我們對模型設(shè)置了兩個(gè)反射界面，一個(gè)是塊狀異常體的上邊界，另一個(gè)是模型的底界面.圖3a為塊狀彈性模量參數(shù)模型的反射波射線路徑分布圖，射線覆蓋了塊狀異常體區(qū)域，圖3b是利用相同的各向異性介質(zhì)正演走時(shí)結(jié)合縱波速度走時(shí)層析成像所得速度反演結(jié)果，反演所得異常體的形狀能夠得到識別，但是異常體的數(shù)值變化非常大，整體效果不太理想.圖4為塊狀異常體模型的彈性模量參數(shù)(a—c—e—g)與反演結(jié)果(b—d—f—h)的對比圖，左側(cè)為真實(shí)彈性模量參數(shù)，右側(cè)為反演結(jié)果，c13,c33,c44均得到了較好的反演效果，由于c11的敏感性較弱，它的反演效果相對較差.如果某個(gè)參數(shù)的速度偏導(dǎo)數(shù)相對于其他參數(shù)的速度偏導(dǎo)數(shù)不靈敏，那么利用觀測走時(shí)數(shù)據(jù)來重建該參數(shù)剖面是不太現(xiàn)實(shí)的.國內(nèi)外文獻(xiàn)均是針對敏感性較強(qiáng)的各向異性參數(shù)進(jìn)行反演，而很難把靈敏度弱的各向異性參數(shù)反演出來.由迭代反演收斂曲線(圖9)可知反演的走時(shí)數(shù)據(jù)擬合較好，說明反演結(jié)果與真實(shí)模型相接近.

3.2 qP反射波走時(shí)反演Thomsen參數(shù)模型

下面分別對層狀介質(zhì)模型和塊狀異常體模型進(jìn)行TI介質(zhì)反射波非線性走時(shí)反演，求取模型的Thomsen參數(shù).

(1)層狀介質(zhì)模型的Thomsen參數(shù)反演.層狀介質(zhì)模型的長度為800 m，深度為600 m，第一層介質(zhì)的Thomsen參數(shù)為α0=2.71 km·s-1,β0=1.51 km·s-1,ε=0.10,δ*=-0.13,γ=0.35，對稱軸傾角為θ0=0°；第二層介質(zhì)的Thomsen參數(shù)為α0=4.72 km·s-1,β0=2.89 km·s-1,ε=0.26,δ*=0.17,γ=0.17，對稱軸傾角為θ0=90°；第三層介質(zhì)的Thomsen參數(shù)為α0=3.37 km·s-1,β0=1.83 km·s-1,ε=0.18,δ*=0.02,γ=0.26，對稱軸傾角為θ0=0°.采用地表激發(fā)和接收的觀測方式，炮點(diǎn)與檢波點(diǎn)位置相同，水平道間距為25 m，炮點(diǎn)與檢波點(diǎn)均為33個(gè)，模型的網(wǎng)格間距為50 m.將第一層參數(shù)作為初始模型用以反演第二層和第三層的Thomsen參數(shù).為了改善反演成像的效果，我們對模型設(shè)置了三個(gè)反射界面，充分利用了三個(gè)反射界面的反射波走時(shí)信息.圖5a為層狀Thomsen參數(shù)模型的反射波射線路徑分布圖，射線覆蓋了整個(gè)模型范圍，圖5b是利用相同的各向異性介質(zhì)正演走時(shí)結(jié)合縱波速度走時(shí)層析成像所得速度反演結(jié)果，反演所得結(jié)果的中間層速度值變化比較大，并且中間層上下兩個(gè)邊界位置不易識別.圖6為層狀介質(zhì)模型的Thomsen參數(shù)(a—c—e—g)與反演結(jié)果(b—d—f—h)的對比圖，左側(cè)為真實(shí)Thomsen參數(shù)，右側(cè)為反演結(jié)果，α0,ε和δ*均得到了較好的反演效果，因?yàn)棣?的敏感性較弱，它的反演效果稍差，其界面位置很難分辨清楚.由迭代反演收斂曲線(圖9)可知觀測走時(shí)與正演走時(shí)之間的殘差較小，反演得到的結(jié)果與真實(shí)模型相接近.

圖1 層狀介質(zhì)彈性參數(shù)模型的射線分布與速度反演結(jié)果Fig.1 The raypath distribution and inverted velocity profile of the elastic moduli layered model

圖2 層狀介質(zhì)模型的彈性模量參數(shù)與反演結(jié)果Fig.2 The true elastic moduli parameters and inverted results for the elastic moduli layered model

圖3 塊狀異常體彈性參數(shù)模型的射線分布與速度反演結(jié)果Fig.3 The raypath distribution and inverted velocity profile of the elastic moduli blocky abnormal body model

圖4 塊狀異常體模型的彈性模量參數(shù)與反演結(jié)果Fig.4 The true elastic moduli parameters and inverted results for the elastic moduli blocky abnormal body model

圖5 層狀介質(zhì)Thomsen參數(shù)模型的射線分布與速度反演結(jié)果Fig.5 The raypath distribution and inverted velocity profile of the Thomsen layered model

圖6 層狀介質(zhì)模型的Thomsen參數(shù)與反演結(jié)果Fig.6 The true elastic moduli parameters and inverted results for the Thomsen layered model

(2)塊狀異常體模型的Thomsen參數(shù)反演.塊狀異常體模型的長度為800m，深度為500m，背景介質(zhì)的Thomsen參數(shù)為α0=3.16 km·s-1,β0=2.0 km·s-1,ε=0.25,δ*=0.3,γ=0.25，對稱軸傾角為θ0=0°；塊狀異常體的Thomsen參數(shù)為α0=3.87 km·s-1,β0=2.83 km·s-1,ε=0.18,δ*=1.14,γ=0.25，對稱軸傾角為θ0=45°.采用地表激發(fā)和接收的觀測方式，炮點(diǎn)與檢波點(diǎn)位置重合，水平道間隔為16 m，炮點(diǎn)與檢波點(diǎn)均為51個(gè)，模型的網(wǎng)格間距為50 m.將背景介質(zhì)作為初始模型用以反演塊狀異常體的Thomsen參數(shù).為了改善反演成像的效果，我們對模型設(shè)置了兩個(gè)反射界面，一個(gè)是塊狀異常體的上邊界，另一個(gè)是模型的底界面.圖7a為塊狀異常體Thomsen參數(shù)模型的反射波射線路徑分布圖，射線覆蓋了塊狀異常體區(qū)域，圖7b是利用相同的各向異性介質(zhì)正演走時(shí)結(jié)合縱波速度走時(shí)層析成像所得速度反演結(jié)果，可以看出速度反演結(jié)果效果非常不理想，完全無法識別異常體的形態(tài)和位置.圖8為塊狀異常體模型的Thomsen參數(shù)(a—c—e)與反演結(jié)果(b—d—f)的對比圖，左側(cè)為真實(shí)Thomsen參數(shù)，右側(cè)為反演結(jié)果，參數(shù)α0,β0和δ*均得到了較好的反演效果，由于塊狀異常體的ε值比背景介質(zhì)的ε值更低，反演算法比較難反演出ε參數(shù).由圖9的收斂曲線可知該模型的走時(shí)反演得到了較好的收斂效果，表明反演模型與真實(shí)模型的相似性較好.

圖7 塊狀異常體Thomsen參數(shù)模型的射線分布與速度反演結(jié)果Fig.7 The raypath distribution and inverted velocity profile of the Thomsen blocky abnormal body model

圖8 塊狀異常體模型的Thomsen參數(shù)與反演結(jié)果Fig.8 The true elastic moduli parameters and inverted results for the Thomsen blocky abnormal body model

圖9 四種模型的迭代反演收斂曲線Fig.9 The iterative inversion convergence curve for these four models

4 結(jié)論

各向異性現(xiàn)象在地球介質(zhì)中普遍存在，地震走時(shí)層析成像需要考慮地層的各向異性.本文的各向異性TI介質(zhì)qP反射波走時(shí)層析成像方法適用于任意強(qiáng)度的各向異性介質(zhì)，該方法相比井間和垂直地震剖面走時(shí)層析成像方法具有更加廣泛的應(yīng)用價(jià)值.由于qP波相速度對彈性模量參數(shù)和Thomsen參數(shù)的速度偏微分不同，所以利用qP反射波可以分別反演出彈性模量參數(shù)和Thomsen參數(shù).在數(shù)值模擬部分，利用qP反射波分別反演層狀介質(zhì)模型和塊狀異常體模型，均得到了較好的反演結(jié)果.這些模型的迭代反演收斂曲線均得到了較好的收斂效果，表明反演結(jié)果與真實(shí)模型之間相接近.與各向同性介質(zhì)走時(shí)反演結(jié)果相比較，各向異性介質(zhì)走時(shí)反演結(jié)果具有較好的識別能力.數(shù)值模擬結(jié)果證明了各向異性TI介質(zhì)qP反射波走時(shí)層析成像算法的正確性.由于受地表觀測系統(tǒng)的限制，反射波走時(shí)層析成像相比井間直達(dá)波走時(shí)層析成像的反演效果較差，在固定檢波器數(shù)量的情況下，應(yīng)最大限度地增加qP反射波射線覆蓋范圍以減小走時(shí)反演的多解性和非線性.根據(jù)反演理論，已知的信息量越多反演的不確定性越弱，如果能夠綜合利用qP波、qSV波和qSH波進(jìn)行多波走時(shí)聯(lián)合反演，將有助于改善各向異性參數(shù)的反演效果，因此各向異性TI介質(zhì)多波走時(shí)聯(lián)合反演值得更深入的研究.

Cerveny V. 2001. Seismic Ray Theory. Cambridge: Cambridge University Press.

Chapman C H, Pratt R G. 1992. Traveltime tomography in anisotropic media—I Theory.GeophysicalJournalInternational, 109(1): 1-19.Crampin S. 1984. Effective anisotropic elastic constants for wave-propagation through cracked solids.GeophysicalJournaloftheRoyalAstronomicalSociety, 76(1): 135-145.

Helbig K. 1981. Systematic classification of layered-induced transverse isotropy.GeophysicalProspecting, 29(4): 550-577.Jech J. 1988. Three-dimensional inversion problem for inhomogeneous transversely isotropic media.StudiaGeophysicaetGeodetica, 32(2): 136-143.Kreyszig E. 1993. Advanced Engineering Mathematics. New York: John Wiley & Sons Inc, 790.

Li F, Cao S Y, Yao J. 2012. Calculation of phase and group velocities in an arbitrary anisotropic medium.ChineseJournalofGeophysics(in Chinese), 55(10): 3420-3426.

Liu Y Z, Wang G Y, Dong L G, et al. 2014. Joint inversion of VTI parameters using nonlinear traveltime tomography.ChineseJournalofGeophysics(in Chinese), 57(10): 3402-3410.

Liu Y Z, Yang J Z. 2014. Offset-weighted seismic tomography aimed at using the first arrival efficiently.GeophysicalProspectingforPetroleum(in Chinese), 53(1): 99-115.Lu M H, Tang J H, Yang H Z, et al. 2005. P-wave traveltime analysis and Thomsen parameters inversion in orthorhombic media.ChineseJournalofGeophysics(in Chinese), 48(5): 1167-1171.Pratt R G, Chapman C H. 1992. Traveltime tomography in anisotropic media—II Application.GeophysicalJournalInternational, 109(1): 20-37.Tryggvason A, Schmelzbach C, Juhlin C. 2009. Traveltime tomographic inversion with simultaneous static corrections—well worth the effort.Geophysics, 74(6): WCB25-WCB33.

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

Wang Y H. 2014. Seismic ray tracing in anisotropic media: A modified Newton algorithm for solving highly nonlinear systems.Geophysics, 79(1): T1-T7.Zelt C A, Azaria A, Levander A. 2006. 3D seismic refraction traveltime tomography at a groundwater contamination site.Geophysics, 71(5): H67-H78.

Zhang W S, He Q D. 1999. Statistical tomography in anisotropic media by using crosshole traveltimes of qP wave and qSV wave.ChineseJournalofGeophysics(in Chinese), 42(S1): 195-203.

Zhou B, Greenhalgh S A, Sinadinovski C. 1992. Iterative algorithms for the damped minimum norm, least squares and constraining problem in seismic tomography.ExplorationGeophysics, 23: 497-505.

Zhou B, Greenhalgh S A. 2005. Analytic expressions for the velocity sensitivity to the elastic moduli for the most general anisotropic media.GeophysicalProspecting, 53(4): 619-641.Zhou B, Greenhalgh S A. 2006. Raypath and traveltime computations for 2D transversely isotropic media with dipping symmetry axes.ExplorationGeophysics, 37(2): 150-159.Zhou B, Greenhalgh S A, Green A. 2008. Nonlinear traveltime inversion scheme for crosshole seismic tomography in tilted transversely isotropic media.Geophysics, 73(4): D17-D33.Zhou H, He Q D. 1995. Seismic traveltime inversion for transversely

isotropic media.GeophysicalProspectingforPetroleum(in Chinese), 34(3): 63-68.Zhou H W. 2006. Multiscale deformable-layer tomography.Geophysics, 71(3): R11-R19.

附中文參考文獻(xiàn)

李芳, 曹思遠(yuǎn), 姚健. 2012. 任意各向異性介質(zhì)相(群)速度的計(jì)算. 地球物理學(xué)報(bào), 55(10): 3420-3426, doi: 10.6038/j.issn.0001-5733.2012.10.025.

劉玉柱, 王光銀, 董良國等. 2014. VTI介質(zhì)多參數(shù)聯(lián)合走時(shí)層析成像方法. 地球物理學(xué)報(bào), 57(10): 3402-3410.

劉玉柱, 楊積忠. 2014. 有效利用初至信息的偏移距加權(quán)地震層析成像方法. 石油物探, 53(1): 99-115.

盧明輝, 唐建侯, 楊慧珠等. 2005. 正交各向異性介質(zhì)P波走時(shí)分析及Thomsen參數(shù)反演. 地球物理學(xué)報(bào), 48(5): 1167-1171, doi: 10.3321/j.issn:0001-5733.2005.05.026.

張文生, 何樵登. 1999. 各向異性介質(zhì)中利用跨孔qP波和qSV波走時(shí)的統(tǒng)計(jì)法成像. 地球物理學(xué)報(bào), 42(S1): 195-203.

周輝, 何樵登. 1995. 橫向各向同性介質(zhì)走時(shí)反演. 石油物探, 34(3): 63-68.

(本文編輯 何燕)

Traveltime tomography of qP reflection waves in anisotropic TI media

HUANG Guang-Nan1,2,3， ZHOU Bing4， DENG Ju-Zhi1， LI Hong-Xing1， LI Ze-Lin1， ZHANG Hua1

1FundamentalScienceonRadioactiveGeologyandExplorationTechnologyLaboratory,EastChinaInstituteofTechnology,Nanchang330013,China2HubeiSubsurfaceMulti-scaleImagingKeyLaboratory(SMIL),ChinaUniversityofGeosciences,Wuhan430074,China3StateKeyLaboratoryofPetroleumResourcesandProspecting,ChinaUniversityofPetroleum,Beijing102249,China4PetroleumGeosciences,thePetroleumInstitute,P.O.Box2533,AbuDhabi,UAE

Seismic traveltime tomography is an effective method to invert underground anisotropic parameters. Research about the anisotropic traveltime tomographic method remains insufficient. Most studies on this issue are based on weak anisotropy assumption which was proposed by Leon Thomsen according to previous research. But the degree of anisotropic property of strata and rocks is uncertain in real applications. Therefore, a kind of traveltime tomography by use of qP reflection waves is proposed to invert arbitrary anisotropic TI media.First arrival traveltime tomography is often used for near-surface and crosswell anisotropic parameters inversion. In fact, traveltime tomography of qP reflection waves has a much more extensive practical value, because the surface observation system is more economic and general than the systems of crossholes and vertical seismic profiles. Here, we present a new nonlinear traveltime inversion method for the surface observation system which combines several important features: (1) A robust reflected wave ray tracing method is used for arbitrary anisotropic TI media; (2) The first-order traveltime perturbation equation is not the eigenvector form of Cerveny′s linearized formula, which suffers from a singularity problem for the two quasi-shear waves; (3) An effective computation of the Jacobian matrix is employed for an arbitrary anisotropic TI media; and (4) It adopts a fast, local minimization search style of nonlinear inversion. These features make the traveltime tomography of qP reflection waves can be used to invert anisotropic TI media with arbitrary anisotropic degree.In numerical simulation section, this seismic reflection tomography method was used to invert for a layered model and a blocky abnormal body model, respectively. The elastic moduli parameter and Thomsen parameter of these anisotropic models were obtained, respectively. Firstly, we used the traveltimes of qP reflection waves to invert the elastic moduli parameters of the layered model and blocky abnormal body model. (1) According to the tomograms of the layered model, the velocity image is not reconstructed very well although the ray paths can cover the whole model perfectly. The velocity values of the first and third layers differ from the true velocity value greatly. And the bottom interface of the second layer is not recovered very well. Fortunately, the elastic moduli parameters are all reconstructed very well. We can see these images are all very close to sections of the true elastic moduli parameters. (2) After inversion of the blocky abnormal body model, the velocity image is not reconstructed very well. Because the velocity values of the blocky abnormal body differ from the true velocity values greatly. On the contrary, the elastic moduli parameters are very close to the profiles of true elastic moduli parameters. But thec11parameter is not as good as the other three parameters. Secondly, we used the traveltimes of qP reflection waves to invert the Thomsen parameters of the layered model and blocky abnormal body model. (1) According to the tomograms of the layered model, the velocity image is still not reconstructed very well, because the velocity values of the second layer differ from the true velocity values greatly. The Thomsen parameters are reconstructed very well, except the image ofβ0is slightly not so good. We can see the inverted and true images ofα0,εandδ*are very close to each other. (2) After inversion of the blocky abnormal body model, the velocity image is not reconstructed correctly, because its numerical value and shape differ from the true velocity model greatly. However, the profiles ofα0,β0andδ*are reconstructed successfully, except theεparameter. The images of these three inverted parameters are very close to that of the true parameters. Besides, these four numerical experiments all reach satisfactory convergence levels for different anisotropic models. Hence, the simulation results show that traveltime of qP reflection waves can invert anisotropic parameters correctly.The traveltimes of qP reflection waves have been used to invert anisotropic parameters for the layered model and blocky abnormal body model successfully. The images of the reconstructed results are very close to the images of the true anisotropic parameters. The tomographic method proposed here has superior capability in recovering theoretical models compared with the isotropic traveltime tomographic method. The more information is used, the less uncertainty of the inversion results according to the inverse theory. In other words, the inversion results can be improved greatly if the traveltimes of qP, qSV and qSH waves are used in the traveltime inversion. Therefore, joint traveltime inversion of qP, qSV and qSH waves should worth further research.

Anisotropic media; Anisotropic parameters; Elastic moduli parameters; Traveltime tomography

10.6038/cjg20150618.

東華理工大學(xué)博士科研啟動(dòng)基金(DHBK2013212)，中國地質(zhì)大學(xué)(武漢)地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室基金(SMIL-2015-10)，國家科技重大專項(xiàng)(2011ZX05024-001-02)，國家自然科學(xué)基金(41104073, 41004048, 41364004, 41104074, 41304097)和國家科技支撐計(jì)劃(2011BAB04B03)聯(lián)合資助.

黃光南，男，1983年生，講師，主要從事地震走時(shí)層析成像方法研究.E-mail:bobking2@126.com

10.6038/cjg20150618

P631

2014-06-25，2015-03-16收修定稿

黃光南，ZHOU Bing，鄧居智等.2015. 各向異性TI介質(zhì)qP反射波走時(shí)層析成像.地球物理學(xué)報(bào)，58(6):2035-2045,

Huang G N, Zhou B, Deng J Z，et al. 2015. Traveltime tomography of qP reflection waves in anisotropic TI media.ChineseJ.Geophys. (in Chinese),58(6):2035-2045,doi:10.6038/cjg20150618.