周嘉欣,熊 凱,楊 鍇

(同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院,上海200092)

直角坐標(biāo)系下立體層析FRECHET導(dǎo)數(shù)求取修正問題探討

周嘉欣,熊 凱,楊 鍇

(同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院,上海200092)

直角坐標(biāo)系下立體層析反演所需的FRECHET導(dǎo)數(shù)是地表觀測(cè)到的擾動(dòng)量對(duì)模型空間初始擾動(dòng)量的一階導(dǎo)數(shù)。當(dāng)傍軸射線沒有到達(dá)地表時(shí),FRECHET導(dǎo)數(shù)的計(jì)算會(huì)產(chǎn)生誤差,從而影響反演精度,必須對(duì)其進(jìn)行修正。利用中心射線地表出射點(diǎn)坐標(biāo)、傍軸射線最后一點(diǎn)坐標(biāo)以及傍軸射線最后一點(diǎn)與其地表預(yù)期出射點(diǎn)坐標(biāo)差之間存在的三角關(guān)系,換算出中心射線地表出射點(diǎn)與傍軸射線地表預(yù)期出射點(diǎn)的坐標(biāo)差和初始位置擾動(dòng)之間的線性關(guān)系,從而求取正確的FRECHET導(dǎo)數(shù)。上述換算與前人提出的多追蹤一小段射線的計(jì)算公式等價(jià),但是計(jì)算成本更低。數(shù)值算例證實(shí)了修正方法的正確性以及直角坐標(biāo)系下立體層析實(shí)施這種修正的必要性。

立體層析反演;傍軸射線;FRECHET導(dǎo)數(shù);直角坐標(biāo)系

在工業(yè)界實(shí)施偏移速度建模的歷史中,基于射線理論的層析反演扮演著非常重要的角色[1-4]。即便在全波形反演成為學(xué)術(shù)界和工業(yè)界研究熱點(diǎn)的今天,層析反演依然不可或缺。學(xué)者們已經(jīng)證明,在全波形反演的第一階段,即反演低波數(shù)速度成分的階段,全波形反演基本上等價(jià)于層析反演[5-6]。在諸多層析反演方法中,立體層析是一種極具特色的層析反演算法,它重新定義了數(shù)據(jù)空間和模型空間,炮檢點(diǎn)坐標(biāo)和對(duì)應(yīng)的射線參數(shù)被納入其中,實(shí)現(xiàn)了對(duì)速度結(jié)構(gòu)、反射點(diǎn)位置和反射層局部形態(tài)的聯(lián)合反演。目前,立體層析提供的高分辨率速度模型不僅能夠滿足偏移成像的需求,而且能夠?yàn)槿ㄐ畏囱菀约皟?chǔ)層反演提供高質(zhì)量的初始模型[7-13]。

傳統(tǒng)立體層析方法是在射線中心坐標(biāo)系下建立層析矩陣方程的。BILLETTE等[13]和LAMBARE[14]在射線中心坐標(biāo)系下推導(dǎo)了二維立體層析所需的模型空間各分量擾動(dòng)與數(shù)據(jù)空間各分量擾動(dòng)之間的線性關(guān)系,成功實(shí)現(xiàn)了射線中心坐標(biāo)系下的二維立體層析反演。在此基礎(chǔ)上,CHALARD等[15]實(shí)現(xiàn)了射線中心坐標(biāo)系下的三維立體層析,在三維射線中心坐標(biāo)系下使用降階漢密爾頓算子推導(dǎo)出三維立體層析所需的各個(gè)導(dǎo)數(shù)關(guān)系。YANG等[16]指出,三維射線擾動(dòng)理論在射線中心坐標(biāo)系下實(shí)現(xiàn)并非最佳選擇,因?yàn)樵谏渚€中心坐標(biāo)系下構(gòu)造三維傍軸射線追蹤所需的傳播矩陣時(shí),涉及到對(duì)垂直于中心射線的各個(gè)方向求關(guān)于速度二階導(dǎo)數(shù)的計(jì)算,這個(gè)計(jì)算非常耗時(shí)。在三維直角坐標(biāo)系下實(shí)現(xiàn)傍軸射線追蹤則不存在上述問題。

盡管在直角坐標(biāo)系下實(shí)施三維傍軸射線追蹤對(duì)于降低計(jì)算成本有很大優(yōu)勢(shì),但是有一個(gè)關(guān)鍵的環(huán)節(jié)——FRECHET導(dǎo)數(shù)矩陣的計(jì)算存在誤差,需要對(duì)其進(jìn)行修正。在直角坐標(biāo)系下實(shí)施傍軸射線追蹤時(shí),除射線垂直出射地表外,一般情況下中心射線與傍軸射線不可能同時(shí)到達(dá)地表。當(dāng)中心射線到達(dá)地表時(shí),傍軸射線往往還沒有到達(dá)。建立立體層析矩陣應(yīng)該使用傍軸射線到達(dá)地表時(shí)其對(duì)應(yīng)的傳播矩陣中的各個(gè)元素,如果采用傍軸射線沒有到達(dá)地表時(shí)的傳播矩陣元素,那么地表接收到的數(shù)據(jù)空間的擾動(dòng)量與地下模型空間的擾動(dòng)量之間的線性關(guān)系將無法通過射線擾動(dòng)理論正確得到。這意味著FRECHET導(dǎo)數(shù)矩陣的計(jì)算將產(chǎn)生誤差,只有對(duì)這個(gè)誤差實(shí)施修正才能正確實(shí)現(xiàn)直角坐標(biāo)系下的立體層析。

NAG[17]提出多追蹤一小段傍軸射線使得傍軸射線到達(dá)地表的思路,需要計(jì)算出對(duì)應(yīng)于多追蹤這一小段的積分步長(zhǎng),它涉及到地表面法向量與坐標(biāo)差向量的點(diǎn)積以及除法運(yùn)算。如果對(duì)每條射線都實(shí)施上述運(yùn)算,則會(huì)增加不少計(jì)算量。通過分析我們發(fā)現(xiàn),當(dāng)中心射線到達(dá)地表但傍軸射線并沒有到達(dá)地表時(shí),中心射線地表出射點(diǎn)坐標(biāo)、傍軸射線最后一點(diǎn)坐標(biāo)以及傍軸射線最后一點(diǎn)與其地表預(yù)期出射點(diǎn)的坐標(biāo)差之間存在一個(gè)近似的三角關(guān)系。利用這個(gè)三角關(guān)系,在傍軸射線沒有到達(dá)地表的情況下依然可以換算出立體層析所需要的中心射線地表出射點(diǎn)坐標(biāo)與傍軸射線地表預(yù)期出射點(diǎn)的坐標(biāo)差和地下射線出射點(diǎn)初始位置之間的線性關(guān)系,從而求取正確的FRECHET導(dǎo)數(shù)。本文首先從理論上證明了上述三角關(guān)系換算公式與NAG[17]多追蹤一小段射線的計(jì)算公式等價(jià),然后通過數(shù)值算例證實(shí)了該換算公式的正確性。理論數(shù)據(jù)反演結(jié)果表明,該換算公式能使直角坐標(biāo)系下立體層析反演的精度得到明顯提高。

1 直角坐標(biāo)系下立體層析線性關(guān)系

立體層析同時(shí)利用了走時(shí)信息和同相軸的局部斜率信息來估計(jì)背景速度模型,因此其數(shù)據(jù)空間和模型空間與傳統(tǒng)走時(shí)層析方法明顯不同。圖1給出了二維立體層析數(shù)據(jù)空間各分量與模型空間各分量。地表能接收到的數(shù)據(jù)為d=(Sx,Rx,TSR,pSx,pRx),其中Sx,Rx為炮檢點(diǎn)的橫坐標(biāo);TSR為射線對(duì)雙程旅行時(shí);pSx,pRx為炮檢點(diǎn)處射線參數(shù)水平分量。一個(gè)局部同相軸可以唯一地對(duì)應(yīng)于一個(gè)炮點(diǎn)和一個(gè)接收點(diǎn)位置Sx和Rx,疊前道集內(nèi)拾取的斜率pSx和pRx,以及雙程走時(shí)TSR。“立體層析”這一提法來源于其數(shù)據(jù)空間引入的炮檢點(diǎn)處射線參數(shù)水平分量pSx和pRx對(duì)反演產(chǎn)生的強(qiáng)約束[14]。模型空間為m=(x0,z0,θS,θR,v),其中x0,z0為地下反射點(diǎn)坐標(biāo);θS表示地下反射點(diǎn)向炮點(diǎn)一側(cè)出射的射線與法線方向的夾角;θR表示地下反射點(diǎn)向檢波點(diǎn)一側(cè)出射的射線與法線方向的夾角;v表示介質(zhì)速度。

圖1 二維立體層析數(shù)據(jù)空間各分量與模型空間各分量

立體層析將傳統(tǒng)反射層析由炮點(diǎn)到反射點(diǎn)再到檢波點(diǎn)的反射過程變?yōu)橛傻叵路瓷潼c(diǎn)分別向炮點(diǎn)和檢波點(diǎn)發(fā)出射線的兩個(gè)透射過程。數(shù)據(jù)空間殘差Δd和模型空間殘差Δm之間的線性關(guān)系可以表示為:

(1)

(1)式中的系數(shù)矩陣F即為二維立體層析的FRECHET導(dǎo)數(shù)矩陣。其完整形式如下:

(2)

式中:σ為尺度量綱因子,用以解決數(shù)據(jù)量級(jí)差別過大的問題。注意地面炮點(diǎn)一側(cè)觀測(cè)到的坐標(biāo)和射線參數(shù)與檢波點(diǎn)一側(cè)對(duì)應(yīng)的地下反射點(diǎn)處張角和方位角無關(guān)。這一事實(shí)在檢波點(diǎn)處同樣存在,因此有

(3)

另由Fermat原理,在射線兩端及背景速度固定不變時(shí),走時(shí)不變且最小。因此模型空間地下散射角的一階擾動(dòng)對(duì)旅行時(shí)TSR不會(huì)造成影響,有

(4)

根據(jù)(3)～(4)式,可將FRECHET導(dǎo)數(shù)矩陣簡(jiǎn)化為:

(5)

注意(5)式中除了旅行時(shí)TSR關(guān)于模型空間的導(dǎo)數(shù)無需通過射線擾動(dòng)理論計(jì)算外,其余各個(gè)分量關(guān)于模型空間的導(dǎo)數(shù)均需要通過射線擾動(dòng)理論計(jì)算,具體推導(dǎo)過程見下文。

2 二維直角坐標(biāo)系下的漢密爾頓傳播算子及傍軸射線追蹤系統(tǒng)

二維直角坐標(biāo)系下漢密爾頓算子形式為:

(6)

為表達(dá)方便,我們不妨用x1表示x,x3表示z,p1表示px,p3表示pz,那么二維直角坐標(biāo)系下的射線追蹤方程組可表示為:

(7)

由漢密爾頓方法,可得其具體形式為:

(8a)

(8b)

方程(8)由4個(gè)一階微分方程組成。當(dāng)某條參考射線(中心射線)的初始射線場(chǎng)有一個(gè)微小的擾動(dòng)量時(shí),擾動(dòng)的射線會(huì)在中心射線附近傳播,該擾動(dòng)射線稱為傍軸射線。傍軸射線通過計(jì)算其相對(duì)中心射線的相空間擾動(dòng)量得到。對(duì)漢密爾頓系統(tǒng)兩端取全微分,可得到直角坐標(biāo)系下的傍軸射線追蹤系統(tǒng):

(9)

將上述傍軸射線追蹤系統(tǒng)寫為如下矩陣形式:

(10)

其中,δη=(δx1,δx3,δp1,δp3)T,

S=

(10)式即為二維直角坐標(biāo)系下的運(yùn)動(dòng)學(xué)傍軸射線追蹤系統(tǒng),表達(dá)了射線場(chǎng)的初始擾動(dòng)[δ(x1)0,δ(x3)0,δ(p1)0,δ(p3)0]與射線軌跡上任意一點(diǎn)的射線場(chǎng)擾動(dòng)(δx1,δx3,δp1,δp3)之間的線性關(guān)系。二維直角坐標(biāo)系漢密爾頓形式下的速度擾動(dòng)(δv)對(duì)漢密爾頓傳播算子的線性影響可以表達(dá)為:

(11)

其中δH為漢密爾頓擾動(dòng)算子:

(12)

(13)

對(duì)(11)式表示的漢密爾頓算子關(guān)于δζ=(δx1,δx3,δp1,δp3)T全微分,可得:

(14)

A矩陣表示背景場(chǎng)的傍軸系數(shù)矩陣;B矩陣表示速度擾動(dòng)造成的傍軸系數(shù)矩陣,各元素可以通過對(duì)擾動(dòng)的漢密爾頓傳播算子關(guān)于射線參數(shù)和空間坐標(biāo)參數(shù)的導(dǎo)數(shù)求得,具體表示為:

根據(jù)GILBERT等[18],可求得(14)式的傳播矩陣形式解為:

(15)

射線傳播矩陣Q-1(t1,t0)是傍軸射線追蹤的一組基本解,即分別代入(1,0,0,0)、(0,1,0,0)、(0,0,1,0)、(0,0,0,1)這4組初始值計(jì)算的4組解構(gòu)成的4×4矩陣:

因此(15)式可以寫為如下顯式形式:

(16)

其中,

用(16)式可以建立二維直角坐標(biāo)系下立體層析所需數(shù)據(jù)空間分量中除旅行時(shí)TSR外所有分量對(duì)于模型分量的一階線性擾動(dòng)關(guān)系。為符合表達(dá)上的習(xí)慣,以下我們依然用x表示x1,z表示x3,px表示p1,pz表示p3,這時(shí)炮點(diǎn)坐標(biāo)和炮點(diǎn)射線參數(shù)與模型空間各分量的線性關(guān)系可以表達(dá)如下:

(17)

其中,

(18)

3 直角坐標(biāo)系下的FRECHET導(dǎo)數(shù)誤差來源及修正思路

如第2節(jié)所述,(17)式提供了二維直角坐標(biāo)系下立體層析所需的除旅行時(shí)對(duì)模型空間外的所有FRECHET導(dǎo)數(shù)的計(jì)算方法。其中傳播矩陣Q中的元素必須是傍軸射線到達(dá)地表時(shí)的元素。但實(shí)際計(jì)算中存在的問題是當(dāng)中心射線到達(dá)地表時(shí)射線追蹤即刻停止,而這時(shí)傍軸射線往往還沒有到達(dá)地表(如圖2所示),使得Q矩陣元素對(duì)應(yīng)的空間位置在x(t)處,而不是我們希望的x(t′)處。這個(gè)明顯的誤差意味著(5)式所示的FRECHET導(dǎo)數(shù)矩陣中的?Sx/?x0,?Sx/?z0,?pSx/?x0,?pRx/?z0,?Rx/?x0,?Rz/?z0,?pRx/?x0,?pRx/?z0將產(chǎn)生錯(cuò)誤。尤其當(dāng)出射到地表的中心射線出射角比較大時(shí)(即地下反射傾角較大時(shí)),這個(gè)錯(cuò)誤將愈發(fā)嚴(yán)重。

圖2 直角坐標(biāo)系下傍軸射線到達(dá)地表處擾動(dòng)關(guān)系

針對(duì)上述問題,NAG[17]提出將中心射線再多追一小段,使得傍軸射線到達(dá)地表。但是如何計(jì)算多追這一段的積分步長(zhǎng)dtu？這里不妨將x(t′)處的相空間表達(dá)為dy,x(t)處的相空間表達(dá)為δy,則

顯然有:

(19)

這里y0代表中心射線追蹤到x0時(shí)的相空間量。借助(6)式漢密爾頓形式,(19)式可進(jìn)一步表示為:

(20)

(21)

由(21)式可知,只要對(duì)中心射線再追蹤出dtu對(duì)應(yīng)的一小段,就可以得到對(duì)應(yīng)的傍軸射線到達(dá)地表的傳播矩陣及出射坐標(biāo)等一系列信息,再根據(jù)(17)式得到正確的FRECHET導(dǎo)數(shù)值。

上述過程即修正FRECHET導(dǎo)數(shù)計(jì)算誤差的常規(guī)途徑。本文根據(jù)中心射線地表出射點(diǎn)坐標(biāo)、傍軸射線最后一點(diǎn)坐標(biāo)以及傍軸射線最后一點(diǎn)與其地表預(yù)期出射點(diǎn)的坐標(biāo)差之間存在的一個(gè)近似三角關(guān)系,直接對(duì)FRECHET導(dǎo)數(shù)公式進(jìn)行修改,使得解決誤差問題的方案更加直觀簡(jiǎn)潔。

4 無需多追蹤一步的修正算法

如果直接將(21)式和(17)式代入(20)式,可得:

(22)

對(duì)(22)式中x0,z0分別求導(dǎo),得:

(23)

(23)式為最新的地表坐標(biāo)以及射線參數(shù)p關(guān)于地下反射點(diǎn)坐標(biāo)(x0,z0)的導(dǎo)數(shù)形式,相當(dāng)于在沒有多追蹤一段的情形下,得到了正確的FRECHET導(dǎo)數(shù)公式。如果不做此修正,FRECHET導(dǎo)數(shù)直接用(17)式進(jìn)行求取,結(jié)果為:

(24)

將(24)式與(23)式進(jìn)行對(duì)比可以看出,沒有考慮到達(dá)地表時(shí)的這一段dtu確實(shí)會(huì)帶來誤差。如果對(duì)(22)式中的px0,pz0分別求導(dǎo),即可得到數(shù)據(jù)空間關(guān)于模型空間中散射角的導(dǎo)數(shù)公式,這里不再贅述其推導(dǎo)過程。這種修正的實(shí)質(zhì)是,盡管求出了dtu,但未必要多追這一段,將其代入(20)式再對(duì)有關(guān)的模型空間分量求導(dǎo)也可以獲得等價(jià)的結(jié)果。

在正演驗(yàn)證本文提出的修正方法前,這里給出本文方法與直接進(jìn)行射線追蹤的幾何等價(jià)性證明。如圖3所示,從O′點(diǎn)繼續(xù)進(jìn)行射線追蹤至B點(diǎn),此時(shí)對(duì)應(yīng)的傍軸射線恰好達(dá)到地表O點(diǎn)處。

圖3 中心射線與傍軸射線幾何關(guān)系

容易看出x=OO′=AO′+OA,注意AB其實(shí)就是(17)式中的δz,顯然有:

(25)

由地表O′點(diǎn)處的射線參數(shù)信息,可以求得射線出射到地表的出射角α,有:

而OA是傍軸擾動(dòng)量在x方向的投影,也可以套用(17)式求得。應(yīng)用上述各條件,有:

(26)

(26)式與(22)式完全一致,從另一個(gè)角度驗(yàn)證了本文提出的方法在幾何上的正確性。

5 數(shù)值實(shí)驗(yàn)

5.1 傍軸射線追蹤公式精度驗(yàn)證

如圖4a所示,從地下坐標(biāo)x=0,z=8000m處出射一條中心射線至地表。顯然中心射線對(duì)應(yīng)的傍軸射線并未到達(dá)地表,傍軸射線最后一點(diǎn)的位置與地表有較大的距離(圖4b),若將這一點(diǎn)對(duì)應(yīng)的擾動(dòng)量當(dāng)作地表接收到的擾動(dòng)量,必然會(huì)影響反演的精度。

圖4 中心射線及其傍軸射線(a)與射線到達(dá)地表時(shí)放大顯示(b)

傳統(tǒng)方法建議將中心射線再多追蹤一段使得傍軸射線的最后一點(diǎn)也到達(dá)地表,以求取正確的擾動(dòng)量。但根據(jù)(20)式所代表的三角關(guān)系,我們認(rèn)為無須多追蹤一小段射線,通過修正的FRECHET導(dǎo)數(shù)(即(23)式)就可以實(shí)現(xiàn)正確的反演。本節(jié)先通過一個(gè)簡(jiǎn)單的數(shù)值算例證明本文的修正公式能夠預(yù)測(cè)出正確的地表觀測(cè)擾動(dòng)量。

如圖4b所示,C點(diǎn)為中心射線到達(dá)地表時(shí)對(duì)應(yīng)的傍軸射線最后一點(diǎn),設(shè)E點(diǎn)為傍軸射線到達(dá)地表時(shí)的位置,將C點(diǎn)分別投影到x,z軸,得到以C點(diǎn)作為地表接收點(diǎn)擾動(dòng)量與實(shí)際E點(diǎn)對(duì)應(yīng)擾動(dòng)量縱橫坐標(biāo)誤差分別為CD,DE。在此數(shù)值算例中,DE=236m,CD=58m。

圖5 對(duì)中心射線多追蹤一段AB

5.2 利用修正FRECHET導(dǎo)數(shù)的反演精度測(cè)試

圖6為采用三角剖分方式建立的常規(guī)垂直梯度多層穹隆模型?；谏渚€追蹤正演方法獲得模擬炮數(shù)據(jù)共計(jì)101炮,炮間距20m,每炮401道,道間距10m。圖7a顯示了某一炮對(duì)應(yīng)的射線路徑;圖7b顯示了反演所用的初始模型。

圖6 多層穹隆模型

圖8顯示了直角坐標(biāo)系下FRECHET導(dǎo)數(shù)未做修正時(shí)的最終反演結(jié)果,FRECHET導(dǎo)數(shù)利用(24)式計(jì)算得到。由圖8可以明顯看出,模型空間中“傾角條”(即圖8中藍(lán)色線)的反演在3～5反射層的構(gòu)造翼部存在較大誤差,最深部的地層產(chǎn)狀也遠(yuǎn)未達(dá)到理想效果,迭代80次后最終的泛函殘差達(dá)到e=3.4×102。

采取本文修正方法((23)式)對(duì)FRECHET導(dǎo)數(shù)進(jìn)行修改之后,直角坐標(biāo)系下的反演結(jié)果大為改進(jìn)。如圖9所示,對(duì)應(yīng)于6個(gè)反射層的構(gòu)造形態(tài)與正演模型相比非常接近,其泛函殘差在10步之內(nèi)即收斂到102數(shù)量級(jí)。迭代80次后最終泛函殘差1.8×10,相比修正前的殘差降低了一個(gè)數(shù)量級(jí)?？梢姳疚男拚椒▽?duì)于直角坐標(biāo)系下的立體層析而言不僅是正確的,而且是不可或缺的。

圖7 基于射線追蹤正演的理論數(shù)據(jù)及用于反演的初始模型a 真實(shí)速度模型及某一炮的射線路徑; b 初始速度模型

圖8 應(yīng)用修正前FRECHET導(dǎo)數(shù)的直角坐標(biāo)系下第80次迭代反演結(jié)果

圖9 應(yīng)用本文修正方法后直角坐標(biāo)系下第80次迭代反演結(jié)果

圖10對(duì)比了修正FRECHET導(dǎo)數(shù)前后前20次迭代目標(biāo)泛函的下降曲線,可以看出,修正后(藍(lán)色虛線)第20次迭代時(shí)的目標(biāo)誤差泛函下降速度已明顯減小,收斂速度加快。這從另一個(gè)角度說明了本文修正算法的正確性。

圖10 修正FRECHET導(dǎo)數(shù)前后前20次迭代目標(biāo)誤差泛函下降速度對(duì)比

6 結(jié)束語

在直角坐標(biāo)系下求取立體層析所需的FRECHET導(dǎo)數(shù)時(shí),必須將傍軸射線的最后一點(diǎn)擾動(dòng)修正至地表,否則會(huì)導(dǎo)致FRECHET導(dǎo)數(shù)求取產(chǎn)生誤差,影響反演精度。不同于多追蹤一小段中心射線以保證相應(yīng)的傍軸射線到達(dá)地表的傳統(tǒng)修正思路,本文根據(jù)中心射線地表出射點(diǎn)坐標(biāo)、傍軸射線最后一點(diǎn)坐標(biāo)以及傍軸射線最后一點(diǎn)與其地表預(yù)期出射點(diǎn)的坐標(biāo)差之間存在的三角關(guān)系,通過一個(gè)簡(jiǎn)潔的換算獲得了FRECHET導(dǎo)數(shù)修正公式,理論推導(dǎo)和數(shù)值實(shí)驗(yàn)證明了這個(gè)修正公式與前人提議的多追一小段射線的修正公式等價(jià),同時(shí)也證明了這個(gè)修正是保證直角坐標(biāo)系下立體層析反演精度不可或缺的環(huán)節(jié)。

[1] AL-YAHYA K.Velocity analysis by iterative profile migration[J].Geophysics,1989,54(6):718-729

[2] STORK C.Reflection tomography in the postmigrated domain[J].Geophysics,2012,77(5):680-692

[3] 邵榮峰,方伍寶,蔡杰雄,等.高斯束層析偏移速度建模方法及應(yīng)用[J].石油物探,2016,55(1):91-99 SHAO R F,FANG W B,CAI J X,et al.A method of migration velocity analysis based on Gaussian beam tomography and its application[J].Geophysical Prospecting for Petroleum,2016,55(1):91-99

[4] 劉玉柱,楊積忠.有效利用初至信息的偏移距加權(quán)地震層析成像方法[J].石油物探,2014,53(1):99-115 LIU Y Z,YANG J Z.Offset-weighted seismic tomography aimed at using the first arrival efficiently[J].Geophysical Prospecting for Petroleum,2014,53(1):99-115

[5] SYMES W W.Migration velocity analysis and waveform inversion[J].Geophysical Prospecting,2008,56(6):765-790

[6] BIONDO B,ALMOMIN A.Simultaneous inversion of full data bandwidth by tomographic full-waveform inversion[J].Geophysics,2014,79(3):WA129-WA140

[7] VINCENT P,LAMBARE G,OPERTO S,et al.Building starting models for full waveform inversion from wide-aperture data by stereotomography[J].Geophysical Prospecting,2013,61(S1):109-137

[8] 倪瑤,楊鍇,陳寶書.立體層析反演方法理論分析與應(yīng)用測(cè)試[J].石油物探,2013,52(2):430-438 NI Y,YANG K,CHEN B S.Stereotomography inversion method theory and application testing[J].Geophysical Prospecting for Petroleum,2013,52(2):430-436

[9] REN L J,SUN X D,LI Z C.The stereotomography based on eigen-wave attributes[R].Beijing:International Geophysical Conference,2014

[10] 李振偉,楊鍇,倪瑤,等.基于立體層析反演的偏移速度建模應(yīng)用研究[J].石油物探,2014,53(4):444-452 LI Z W,YANG K,NI Y,et al.Migration velocity analysis with stereo-tomography[J].Geophysical Prospecting for Petroleum,2014,53 (4):444-452

[11] 王宇翔,楊鍇,楊小椿,等.基于梯度平方結(jié)構(gòu)張量算法的高密度二維立體層析反演[J].地球物理學(xué)報(bào),2016,59(1):263-276 WANG Y X,YANG K,YANG X C,et al.A high-density stereo-tomography method based on the gradient square structure tensors algorithm[J].Chinese Journal of Geophysics,59(1):263-276

[12] BRILLATZ C,VIGEE L,COLEOU T,et al.Getting closer to the geology contrasts by starting stratigraphic inversion with velocity models from new tomographic methods[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:4763-4767

[13] BILLETTE F,LAMBARE G.Velocity macro-model estimation from seismic reflection data by stereotomography[J].Geophysical Journal International,1998,135(2):671-690

[14] LAMBARE G.Stereotompgraphy[J].Geophysics,2008,73(5):25-34

[15] CHALARD E,PODVIN P,LAMBARE G.Estimation of velocity macro-models by 3D stereotomography[J].Expanded Abstracts of 70thAnnual Internat SEG Mtg,2000:2257-2260

[16] YANG K,XING F Y,LI Z W.3D stereotomography applied to the deep-sea data acquired in the South China Sea,part I:theoretical foundation[J].Expanded Abstracts of 86thAnnual Internat SEG Mtg,2016:4263-4267

[17] NAG S.PP/PS anisotropic stereotomography[J].Journal of Geophysics International,2010,181(1):427-452

[18] GILBERT F,BACKUS G E.Propagator matrices in elastic wave and vibration problems[J].Studia Geophysica et Geodaetica,1966,10(3):271-271

(編輯：戴春秋)

On the correction of the Frechet derivatives of stereotomography in Cartesian coordinate

ZHOU Jiaxin,XIONG Kai,YANG Kai

(SchoolofOceanandEarthScience,TongjiUniversity,Shanghai200092,China)

The Frechet derivative of stereotomography in the Cartesian coordinate system is a first-order derivative of disturbance quantity observed in surface to initial model space.If the endpoint of the paraxial ray has to be corrected to the surface,a significant error will be introduced into the Frechet derivative.The inversion accuracy will be reduced and should be modified.Differing from the conventional method to solve this problem by tracing the central ray a little longer,we use the triangular relationship among the coordinates of the emergence point for the central ray at the surface,the endpoint of the paraxial ray at the surface and the lateral difference between the endpoint of the paraxial ray and the expected emergence point at the surface,to convert a linear relationship between the endpoint coordinate difference of the central ray and paraxial ray in the surface and the initial perturbation of the ray start-point.Then the correct Frechet derivatives can be achieved based on this conversion.The method presented is exactly equivalent to the previous method with a lower computational cost.Numerical examples demonstrate the correctness and the necessity of the presented correction for stereotomography in the Cartesian coordinate system.

stereotomography,linear relationship,Frechet derivatives,Cartesian coordinate

2016-08-22;改回日期：2016-11-15。

周嘉欣(1992—),女,碩士在讀,主要研究方向?yàn)榱Ⅲw層析反演。

楊鍇(1972—),男,教授,博士生導(dǎo)師,主要研究方向?yàn)榈卣鸩ǖ膫鞑ァ⒊上衽c反演。

國(guó)家自然科學(xué)基金面上項(xiàng)目(41574098,41630964)資助。

P631

A

1000-1441(2017)04-0491-09

10.3969/j.issn.1000-1441.2017.04.004

This research is financially supported by the National Natural Science Foundation of China (Grant Nos.41274117,41574098).