射線雅可比的直接數(shù)值計(jì)算方法及應(yīng)用

梁 全 毛偉建 歐陽威 李武群 張建磊

(①中國科學(xué)院測量與地球物理研究所計(jì)算與勘探地球物理研究中心，湖北武漢 430077；②大地測量與地球動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，湖北武漢430077； ③中國科學(xué)院大學(xué)，北京100049；④東方地球物理公司物探技術(shù)研究中心，河北涿州 072751)

0 引言

地震偏移成像技術(shù)是現(xiàn)今地震數(shù)據(jù)處理的核心技術(shù)之一，在油氣資源探測中起重要作用[1]?？碧斤L(fēng)險(xiǎn)和勘探難度的不斷加大使偏移成像的振幅保真受到重視[2-5]，真振幅疊前深度偏移能獲得更準(zhǔn)確的地下介質(zhì)物性參數(shù)信息，提高油氣勘探精度。地震波真振幅偏移成像旨在得出能更準(zhǔn)確地反映地下反射系數(shù)或介質(zhì)參數(shù)擾動(dòng)的成像振幅[6-7]，從而能進(jìn)行高精度儲(chǔ)層反演。

研究表明，基于線性化反演理論的最小二乘偏移是實(shí)現(xiàn)真振幅成像的有效方法[8]，也是今后保幅成像方法的發(fā)展趨勢[9]。但因當(dāng)前計(jì)算機(jī)硬件條件難以承受Hessian矩陣數(shù)量龐大的元素、迭代更新求解反問題的極高計(jì)算成本、Hessian算子難以精確求逆等，極大地限制了最小二乘偏移方法在此領(lǐng)域的推廣應(yīng)用[9]。作為最小二乘偏移成像特例的常規(guī)一次成像(即Hessian矩陣為單位陣)在實(shí)際數(shù)據(jù)成像處理中仍具現(xiàn)實(shí)意義。對于常規(guī)一次成像，真振幅偏移成像的本質(zhì)是波場傳播的準(zhǔn)確模擬。

射線理論真振幅偏移成像的基本思想[1，3，6]是通過加權(quán)繞射疊加抵消波的傳播引起的振幅變化，提高成像結(jié)果的準(zhǔn)確性。射線理論中，基于射線雅可比計(jì)算輸運(yùn)方程解的方法給出了地震波沿射線的格林函數(shù)振幅表達(dá)式[3，6，10]，該格林函數(shù)普遍存在于射線[11-19]及高斯束類偏移中[20-23]。因此，射線雅可比對于積分法真振幅偏移成像意義重大。

通常利用射線追蹤計(jì)算射線雅可比。傳統(tǒng)射線追蹤方法旨在實(shí)現(xiàn)走時(shí)和射線路徑的快速計(jì)算[24-27]，較少涉及雅可比的計(jì)算，主要應(yīng)用于塊狀構(gòu)造模型的快速正演模擬。為適應(yīng)疊前偏移成像的需求，新型網(wǎng)格類射線追蹤方法不斷涌現(xiàn)[28-31]，按其基本原理可分為最短路徑法[32-33]、有限差分法[34-35]和波前構(gòu)建法[28，36]等三種。最短路徑法和有限差分法側(cè)重于快速計(jì)算網(wǎng)格模型中單次波至走時(shí)。在保幅偏移成像處理中，需依賴走時(shí)場在空間上的二階偏導(dǎo)數(shù)信息計(jì)算所涉及的波前幾何擴(kuò)散[37-40]。

Vanelle等[40]從計(jì)算效率和數(shù)據(jù)存儲(chǔ)的角度對基于走時(shí)場計(jì)算波前幾何擴(kuò)散的方法做了系統(tǒng)研究，即從走時(shí)場平方的泰勒展開式出發(fā)，提出了基于雙曲線走時(shí)插值計(jì)算波前幾何擴(kuò)散的方法。對于此雙曲線走時(shí)插值方法，在存在多波至的復(fù)雜區(qū)域，走時(shí)場在空間的一階偏導(dǎo)數(shù)和二階偏導(dǎo)數(shù)均不連續(xù)(泰勒展開的條件不滿足)，通過插值系數(shù)得出的波前幾何擴(kuò)散因子與動(dòng)力學(xué)射線追蹤所得結(jié)果存在較大差異[40]。提高用雙曲線走時(shí)插值方法計(jì)算波前擴(kuò)散的精度需考慮多波至。然而，在用網(wǎng)格差分計(jì)算走時(shí)空間二階偏導(dǎo)數(shù)時(shí)，多波至走時(shí)的利用十分復(fù)雜。基于走時(shí)場計(jì)算波前幾何擴(kuò)散的方法涉及走時(shí)場的空間二階偏導(dǎo)數(shù)的計(jì)算[10]，該二階偏導(dǎo)數(shù)需根據(jù)網(wǎng)格節(jié)點(diǎn)上的走時(shí)進(jìn)行差分?jǐn)?shù)值計(jì)算而獲得。網(wǎng)格節(jié)點(diǎn)上的走時(shí)精度通常不高，差分計(jì)算走時(shí)場的二階偏導(dǎo)數(shù)不很準(zhǔn)確，且難以較好地處理復(fù)雜區(qū)域多波至問題，因此這種依賴走時(shí)場空間二階偏導(dǎo)數(shù)的方法不適合波前幾何擴(kuò)散的高精度計(jì)算[10]。

Popov等[41]最早提出計(jì)算波前幾何擴(kuò)散的動(dòng)力學(xué)射線追蹤系統(tǒng)，并將其廣泛應(yīng)用于波前構(gòu)建法射線追蹤的振幅計(jì)算和各種射線類保幅疊前深度偏移[12，36，41-42]，在高斯束偏移中尤為常見[20-23]。動(dòng)力學(xué)射線追蹤通過求解動(dòng)力學(xué)常微分方程組實(shí)現(xiàn)射線雅可比的計(jì)算[10，41]。該方法基于射線一階擾動(dòng)理論[43]。在速度擾動(dòng)較小的變速介質(zhì)中，通過動(dòng)力學(xué)射線追蹤能快速有效地計(jì)算射線雅可比； 但當(dāng)速度擾動(dòng)較大時(shí)，其空間上二階導(dǎo)數(shù)很不穩(wěn)定，這時(shí)動(dòng)力學(xué)射線追蹤求得的射線雅可比難以準(zhǔn)確反映波前的實(shí)際變化，導(dǎo)致得到偽焦散點(diǎn)，遺漏真焦散點(diǎn)。

為實(shí)現(xiàn)更準(zhǔn)確的射線雅可比計(jì)算以提高GRT逆散射偏移成像[6，7，44-45]精度，本文從射線雅可比定義出發(fā)[10]，提出利用射線管橫截面元直接數(shù)值計(jì)算的方法求取射線雅可比。即通過數(shù)值計(jì)算射線管橫截面元與初始角面元二者的面積比得出射線雅可比，所得結(jié)果具有明確的物理意義，能準(zhǔn)確識別射線上的真實(shí)焦散點(diǎn)。在逆散射保幅偏移成像條件中，權(quán)函數(shù)的分母項(xiàng)存在兩個(gè)射線雅可比[6-7]，射線的焦散會(huì)導(dǎo)致權(quán)函數(shù)奇異。為此，文中繪出一種合理的雅可比平滑閾值方法，避免了焦散引起的奇異問題?；诓煌乃俣饶Ｐ?，對動(dòng)力學(xué)射線追蹤和直接數(shù)值計(jì)算兩種方法計(jì)算的射線雅可比進(jìn)行了測試和比較。結(jié)果表明：在速度擾動(dòng)較小的介質(zhì)中，上述兩種方法得出的射線雅可比一致； 在速度擾動(dòng)較大的介質(zhì)中，動(dòng)力學(xué)射線追蹤難以獲得較準(zhǔn)確的射線雅可比，直接數(shù)值計(jì)算的射線雅可比則能準(zhǔn)確地反映波前實(shí)際變化。最后將兩種射線雅可比計(jì)算方法分別應(yīng)用于鹽丘模型的逆散射保幅偏移成像，直接數(shù)值計(jì)算法獲得的成像剖面更清晰，進(jìn)一步驗(yàn)證了射線雅可比直接數(shù)值計(jì)算法的有效性和適用性。

1 射線雅可比的定義與計(jì)算

射線雅可比是射線的參數(shù)坐標(biāo)與直角坐標(biāo)之間轉(zhuǎn)換的雅可比矩陣行列式，具有與射線管面元相關(guān)的物理定義[10]，通常用J表示射線雅可比。在速度場具有連續(xù)的一階和二階偏導(dǎo)數(shù)的假設(shè)下，Popov等[41]基于擾動(dòng)理論導(dǎo)出了計(jì)算射線雅可比的動(dòng)力學(xué)射線追蹤方程組，該方程組涉及兩個(gè)傳輸矩陣P和Q的求解，其中以Q表示與波前幾何擴(kuò)散有關(guān)的射線雅可比矩陣、 detQ表示射線雅可比。理論上，detQ=J，但速度模型復(fù)雜時(shí)，動(dòng)力學(xué)射線追蹤法計(jì)算的射線雅可比與物理定義的直接數(shù)值計(jì)算方法求得結(jié)果的差異不可忽略。這里討論兩種射線雅可比計(jì)算方法的原理，并分析各自的特點(diǎn)。

1.1 動(dòng)力學(xué)射線追蹤

動(dòng)力學(xué)射線追蹤是計(jì)算射線雅可比、得出波前幾何擴(kuò)散和射線振幅的一種常用途徑，它在射線及高斯束類保幅偏移成像中的應(yīng)用十分普遍。動(dòng)力學(xué)射線追蹤滿足以下微分方程組[41]

(1)

(2)

式中m和n表示與中心射線垂直的平面坐標(biāo)。

設(shè)et、em和en為構(gòu)成射線中心坐標(biāo)系的一組正交基向量，且這三個(gè)基向量依次表示中心射線的單位切向量(傾角和方位角分別為θ和φ)、m坐標(biāo)的單位切向量和n坐標(biāo)的單位切向量，則這組正交基向量可由θ和φ表示為

(3)

進(jìn)而可得式(2)的矩陣V中各個(gè)二階偏導(dǎo)數(shù)的具體表達(dá)式為

(4)

(5)

此時(shí)

Q反映波前的幾何擴(kuò)散情況。式(1)基于射線的一階擾動(dòng)理論，介質(zhì)的速度擾動(dòng)較小時(shí)，根據(jù)該式進(jìn)行動(dòng)力學(xué)射線追蹤能夠較準(zhǔn)確地得出波場的動(dòng)力學(xué)信息； 當(dāng)介質(zhì)的速度場復(fù)雜時(shí)，二階導(dǎo)數(shù)矩陣V變化較大，得出的動(dòng)力學(xué)參數(shù)也會(huì)不穩(wěn)定和不準(zhǔn)確。在保幅疊前深度偏移成像中，為了使計(jì)算的射線振幅更穩(wěn)定，通常會(huì)將速度場的二階導(dǎo)數(shù)進(jìn)行一定的限制，比如直接將其設(shè)置為零，但這樣的處理不利于成像振幅保真。

1.2 射線雅可比的直接數(shù)值計(jì)算

射線雅可比是表示射線上點(diǎn)的射線參數(shù)坐標(biāo)γ1，γ2，u與廣義笛卡爾坐標(biāo)x，y，z之間轉(zhuǎn)換的行列式。γ1和γ2為決定射線初始入射角的射線參數(shù)，u為沿著射線的單調(diào)參數(shù)，u通常為射線長度s或時(shí)間T。從射線參數(shù)坐標(biāo)到笛卡爾坐標(biāo)轉(zhuǎn)換的雅可比行列式為[10]

(6)

(7)

J=Ω(u)·e

(8)

于是，可根據(jù)

J=guJ(u)=J(s)=gTJ(T)

(9)

計(jì)算關(guān)于任意參數(shù)u的雅可比J(u)(u=T，s)。

根據(jù)矢量微分可得

dΩ(u)=Ω(u)dγ1dγ2

(10)

式中： dΩ(u)表示射線管截等u面的矢量微元面； dγ1、dγ2為初始入射角微元面。當(dāng)u=T時(shí)，矢量微元面有著明確的物理意義，這時(shí)dΩ(u)=dΩ(T)表示射線管截波前面的矢量微元面。dΩ(T)在射線雅可比J的計(jì)算中起著關(guān)鍵作用。

聯(lián)合式(7)、式(9)和式(10)，可得

(11)

式中：N是垂直于波前面的單位法向量； dΩ(T)的單位向量為±N，且有dΩ(T)=dΩ(T)·N。

由于e為射線的單位切向量，N為波前面單位法向量，因此e·N=vp/vg，其中vp為相速度，vg為群速度。以有限差分代替式(11)中微分，得到基于射線管面元的射線雅可比直接數(shù)值計(jì)算方法

(12)

式中： ΔS(T)表示T時(shí)刻射線管截波前面的橫截面元； ΔΩ表示射線管初始角面元，該初始角面元是射線管初始切線構(gòu)成的圓錐側(cè)面截單位球面所得面元。

式(12)中初始角面元的數(shù)值計(jì)算可通過計(jì)算一個(gè)側(cè)棱長為一個(gè)單位的正三棱錐的底面面積實(shí)現(xiàn)，圖1是初始角面元面積數(shù)值計(jì)算示意圖。設(shè)中心射線的初始切向量與輔助射線的初始切向量夾角(圖1中紅色向量與黃色向量夾角)為θless(為小量)，由空間幾何知識易得中心射線對應(yīng)的射線管初始角面元面積(圓錐底面上的正三角形面積)為

(13)

基于上式，式(12)中ΔS(T)表示圖1中以黃色的向量為初始切向量的三條輔助的射線在T時(shí)刻的三個(gè)端點(diǎn)構(gòu)成的三角形面積。在二維介質(zhì)情況下，需要計(jì)算兩條輔助的射線路徑，式(12)中的面元ΔS(T)和角面元ΔΩ分別退化為線元和初始角元(兩條輔助射線的初始切線截單位圓所得的弧長)。

圖1 初始角面元計(jì)算示意圖

圓錐頂點(diǎn)代表炮點(diǎn)； 紅色向量(所在直線為圓錐軸)代表中心射線的初始入射切向量； 三個(gè)黃色的單位向量(終點(diǎn)構(gòu)成正三角形)代表三條輔助射線的初始入射切向量式(12)很好地避免了動(dòng)力學(xué)射線追蹤(式(1))計(jì)算波前幾何擴(kuò)散涉及到的速度場二階偏導(dǎo)數(shù)不穩(wěn)定問題，是具有明確物理意義的射線雅可比數(shù)值計(jì)算方法，求得的結(jié)果具有很好的準(zhǔn)確性和穩(wěn)定性。式(12)不僅能夠較好地應(yīng)用于各向同性介質(zhì)背景下射線雅可比的計(jì)算，而且還能夠很好地適應(yīng)各向異性介質(zhì)背景下波前幾何擴(kuò)散和振幅的計(jì)算，可有效避免求解復(fù)雜的各向異性動(dòng)力學(xué)射線追蹤方程組中存在的簡化計(jì)算和不準(zhǔn)確問題。

2 兩種射線雅可比計(jì)算方法的比較

射線雅可比既可為正，也可為負(fù)，還可以為零。雅可比為零表示射線族產(chǎn)生交叉，局部波前面匯聚成一點(diǎn)，這樣的點(diǎn)稱為焦散點(diǎn)。對相同的射線，分別應(yīng)用動(dòng)力學(xué)射線追蹤方法和直接數(shù)值計(jì)算方法計(jì)算射線雅可比，通過得出的雅可比曲線及其射線族展開，比較兩種方法的準(zhǔn)確性。

圖2為設(shè)計(jì)的變化較平緩的速度模型及以該模型背景計(jì)算的4條射線，模型的橫向和深度方向均有251個(gè)采樣點(diǎn)，橫向采樣間隔為40m， 深度方向采樣間隔為20m， 炮點(diǎn)在地表橫向第126個(gè)采樣點(diǎn)上，從左到右四條射線的入射角分別為-27°、 -9°、 9°、 27°。圖3給出了兩種方法求得的4條射線隨時(shí)間變化的雅可比曲線。在圖3的4個(gè)子圖中，藍(lán)色虛線是基于動(dòng)力學(xué)射線追蹤方法計(jì)算的結(jié)果，紅色實(shí)線是基于直接數(shù)值計(jì)算方法計(jì)算的結(jié)果，二者重合。可見對于每條射線，動(dòng)力學(xué)射線追蹤與直接數(shù)值計(jì)算兩種方法得出的射線雅可比相同。

圖4給出了鹽丘速度模型和以該模型為背景計(jì)算的4條射線，模型橫向有649個(gè)采樣點(diǎn)，采樣間隔為24.384m(80ft)，深度方向有300個(gè)采樣點(diǎn)，采樣間隔為12.192m(40ft)。射線追蹤的炮點(diǎn)位于地表橫向第421個(gè)采樣點(diǎn)上，從左到右4條射線的入射角分別為-36°、-13.5°、5.2°、31.5°。圖5給出了應(yīng)用兩種方法求得的4條射線的雅可比曲線。在圖5的4個(gè)子圖中，藍(lán)色的虛線是動(dòng)力學(xué)射線追蹤方法計(jì)算的結(jié)果，紅色的實(shí)線是直接數(shù)值計(jì)算方法計(jì)算的結(jié)果?？梢娚渚€①的兩條雅可比曲線與射線④的兩條雅可比曲線仍具有較好的一致性，射線②和射線③的各兩條雅可比曲線均呈現(xiàn)出較大的偏差。對于射線②，動(dòng)力學(xué)射線追蹤方法計(jì)算的雅可比曲線經(jīng)過了原點(diǎn)以外的零點(diǎn)，直接數(shù)值計(jì)算方法計(jì)算的雅可比曲線沒有經(jīng)過原點(diǎn)以外的零點(diǎn)。為更好地分析射線②兩種方法計(jì)算的雅可比的差異，以這條射線為中心進(jìn)行射線族展開，圖6為兩種方法計(jì)算的雅可比曲線(圖6a)及其射線族(圖6b)。在[-13.7°，-13.3°]內(nèi)以0.04°的采樣間隔給定了11個(gè)入射角計(jì)算得出的射線族(射線族所在的時(shí)間段包括0.75s以后圖6a中藍(lán)色曲線的零點(diǎn)所在時(shí)刻)，從圖6b中的射線族可以看出，這個(gè)射線族中的射線沒有產(chǎn)生交叉，表明射線上實(shí)際并不存在焦散點(diǎn)。因此，圖6a中正確的雅可比曲線不應(yīng)該經(jīng)過原點(diǎn)以外的零點(diǎn)，動(dòng)力學(xué)射線追蹤方法計(jì)算的雅可比得出了偽焦散點(diǎn)，與實(shí)際不符，直接數(shù)值計(jì)算方法計(jì)算的結(jié)果更準(zhǔn)確。與圖6相反，圖7為圖4中射線③兩種方法計(jì)算的雅可比曲線(圖7a)及其射線族展開(圖7b)，射線族中產(chǎn)生的交叉現(xiàn)象是射線上存在焦散點(diǎn)的典型特征。因此，圖7a中正確的雅可比曲線應(yīng)該經(jīng)過原點(diǎn)以外的零點(diǎn)，圖中可見動(dòng)力學(xué)射線追蹤方法計(jì)算的雅可比曲線未能經(jīng)過原點(diǎn)以外的零點(diǎn)，直接數(shù)值計(jì)算方法計(jì)算的雅可比則準(zhǔn)確識別出了焦散點(diǎn)。

圖2 速度分布及對應(yīng)的4條射線(標(biāo)號為①～④)

圖3 對應(yīng)圖2中的4條射線兩種方法計(jì)算的雅可比曲線

圖4 速度分布及4條射線(標(biāo)號為①～④)位置

圖5 對應(yīng)圖4中的4條射線兩種方法計(jì)算的雅可比曲線

由以上測試結(jié)果可知，對于介質(zhì)速度擾動(dòng)較小的情況，兩種射線雅可比計(jì)算方法獲得的雅可比曲線能很好地符合。當(dāng)介質(zhì)速度變化較大時(shí)，速度場的二階偏導(dǎo)數(shù)不穩(wěn)定，動(dòng)力學(xué)射線追蹤計(jì)算的射線雅可比難以保證準(zhǔn)確性；這時(shí)基于物理定義的射線雅可比直接數(shù)值計(jì)算方法無需計(jì)算速度場的二階偏導(dǎo)數(shù)，求得的雅可比具有更好的穩(wěn)定性和準(zhǔn)確性。

圖6 中心射線②的雅可比曲線(a)及其射線族展開(b，射線族上不存在焦散點(diǎn))

圖7 中心射線③的雅可比曲線(a)及其射線族展開(b，射線族上存在焦散點(diǎn))

由于射線雅可比的直接數(shù)值計(jì)算方法需要額外計(jì)算輔助的射線路徑(二維介質(zhì)需要計(jì)算兩條輔助射線，三維介質(zhì)需要計(jì)算三條輔助射線)，因此，計(jì)算效率較動(dòng)力學(xué)射線追蹤稍有降低，但在可接受的范圍內(nèi)。表1給出了圖2的二維介質(zhì)中相同的參數(shù)下兩種射線雅可比計(jì)算方法的耗時(shí)，可見，兩種方法的計(jì)算效率在同一個(gè)數(shù)量級。

表1 二維介質(zhì)中兩種射線雅可比計(jì)算方法的耗時(shí)

圖8為三維介質(zhì)的速度剖面顯示(左)及其中一條射線基于兩種方法計(jì)算的雅可比曲線(右)，表2給出了該圖中相同的參數(shù)下兩種射線雅可比計(jì)算方法的耗時(shí)，可見，三維情況下兩種方法的計(jì)算效率也在同一個(gè)數(shù)量級。

無論是動(dòng)力學(xué)射線追蹤還是直接數(shù)值計(jì)算，求取射線追蹤表的計(jì)算耗時(shí)相對整個(gè)偏移成像的耗時(shí)來說是很小的，兩種雅可比計(jì)算方法在計(jì)算效率上的差異可以忽略。對于基于物理定義的直接數(shù)值計(jì)算方法，控制中心射線初始切向量與輔助射線的初始切向量夾角θless充分小(通常可設(shè)為0.001°～0.005°)，能夠很好地保證中心射線與輔助射線的一致性。在一些較極端的情況下，即使對于很小的θless，輔助射線可能也會(huì)出現(xiàn)較大的偏折，導(dǎo)致射線族失去一致性，這是介質(zhì)復(fù)雜時(shí)運(yùn)動(dòng)學(xué)射線追蹤計(jì)算的射線場存在固有的“陰影區(qū)”引起的，不是直接數(shù)值計(jì)算方法本身的問題；對于這種很極端的情況，射線雅可比的較準(zhǔn)確計(jì)算涉及更加復(fù)雜的射線理論。

表2 三維介質(zhì)中兩種射線雅可比計(jì)算方法的耗時(shí)

基于射線雅可比的直接數(shù)值計(jì)算方法具有明確的物理意義，避免了動(dòng)力學(xué)射線追蹤涉及的不穩(wěn)定的速度場空間二階偏導(dǎo)數(shù)的計(jì)算。直接數(shù)值計(jì)算方法求得的射線雅可比能夠更好地反映波前的變化，準(zhǔn)確識別出真正的焦散點(diǎn)，通過對雅可比進(jìn)行合理的閾值能夠有效避免保幅疊前深度偏移中焦散引起的奇異問題。

圖8 三維速度剖面(a)與其中的一條射線基于兩種方法計(jì)算的射線雅可比曲線(b)

3 逆散射保幅偏移成像

射線雅可比的計(jì)算對于GRT逆散射保幅偏移成像非常重要，該雅可比是其真振幅成像條件的組成部分[6-7，44-45]，復(fù)雜介質(zhì)中射線雅可比計(jì)算的準(zhǔn)確性會(huì)影響成像。將射線雅可比的直接數(shù)值計(jì)算方法應(yīng)用到逆散射保幅偏移成像，對比根據(jù)動(dòng)力學(xué)射線追蹤得到的成像結(jié)果，進(jìn)一步說明雅可比直接數(shù)值計(jì)算方法的有效性與適用性。

3.1 GRT逆散射保幅偏移成像基本原理

根據(jù)散射理論，地下介質(zhì)被分解成背景介質(zhì)和擾動(dòng)介質(zhì)，地表記錄的散射波場是入射波場作用于擾動(dòng)介質(zhì)的產(chǎn)物。在一階Born近似的條件下，用波場振幅的射線理論解代替格林函數(shù)，地表的聲波散射場滿足以下積分方程[6]

δ[t-φ(s,x,r)]

(14)

I(x)=(R*Us)(x)

(15)

(16)

式中：B為權(quán)函數(shù)；J(x,s)與J(x,r)分別為炮點(diǎn)到成像點(diǎn)的射線雅可比與檢波點(diǎn)到成像點(diǎn)的射線雅可比；θmig為成像點(diǎn)的偏移張角；θs為炮點(diǎn)的射線傾角；θr為檢波點(diǎn)的射線傾角；c0(x)為成像點(diǎn)的聲波速度。權(quán)函數(shù)B的分母中存在兩個(gè)射線雅可比，在焦散點(diǎn)處雅可比為零，權(quán)函數(shù)奇異。為避免焦散引起的奇異值問題，給出以下雅可比平滑閾值處理方案

(17)

其中

(18)

式中：v0(s)與v0(r)分別為炮點(diǎn)與檢波點(diǎn)的波速； dt為射線追蹤的時(shí)間采樣間隔；Jdt(s)和Jdt(r)分別為關(guān)于炮點(diǎn)與檢波點(diǎn)的基準(zhǔn)雅可比；ε為雅可比閾值參數(shù)(通常設(shè)為0.1%～1.0%)；J′為平滑閾值處理后的雅可比。以J′代替J，則式(16)變?yōu)?/p>

(19)

式(19)能夠克服在焦散點(diǎn)附近過小的J引起的奇異問題。另外，當(dāng)J較大時(shí)，J′→J，保證了非焦散區(qū)域計(jì)算結(jié)果的準(zhǔn)確性。

3.2 偏移成像數(shù)值試驗(yàn)

分別將基于動(dòng)力學(xué)射線追蹤和直接數(shù)值計(jì)算方法求取的射線雅可比應(yīng)用于SEG/EAGE鹽丘模型(其各種參數(shù)均與圖4對應(yīng)參數(shù)相同)的偏移成像?？傆?jì)模擬了325炮的地震數(shù)據(jù)，炮間距為48.768m(160ft)，每炮2593道接收，道間距為6.096m(20ft)。從基于動(dòng)力學(xué)射線追蹤方法計(jì)算射線雅可比獲得的成像剖面(圖9a)可看出，在鹽丘的內(nèi)部和下方出現(xiàn)了一些離散的隨機(jī)噪聲。將基于直接數(shù)值計(jì)算方法計(jì)算射線雅可比獲得的成像剖面(圖9b)與圖9a對比，可見隨機(jī)噪聲得到了很好的壓制，獲得了更清晰的成像剖面。因此，直接數(shù)值計(jì)算方法在復(fù)雜介質(zhì)中的射線雅可比計(jì)算具有更強(qiáng)適應(yīng)性。

圖9 基于動(dòng)力學(xué)射線追蹤(a)和直接數(shù)值計(jì)算(b)兩種射線雅可比求取方法所得逆散射保幅偏移成像剖面(左)及局部放大(右)

4 結(jié)論

射線雅可比在計(jì)算波前幾何擴(kuò)散因子和射線格林函數(shù)振幅的過程中起到了關(guān)鍵的作用。該雅可比對于基于射線理論的真振幅偏移成像意義重大，是GRT逆散射偏移真振幅成像條件權(quán)函數(shù)的組成部分。在復(fù)雜介質(zhì)中，速度場的二階偏導(dǎo)數(shù)變化較大，動(dòng)力學(xué)射線追蹤計(jì)算的雅可比不準(zhǔn)確，難以識別真正的焦散點(diǎn)，會(huì)導(dǎo)致形成偽焦散點(diǎn)。為更好地實(shí)現(xiàn)GRT逆散射保幅偏移成像，本文研究了更準(zhǔn)確的射線雅可比計(jì)算方法，從其物理定義出發(fā)，利用射線管的橫截面元直接數(shù)值計(jì)算的方法求取射線雅可比。直接數(shù)值計(jì)算方法避免了速度場二階偏導(dǎo)數(shù)的計(jì)算，得出的雅可比更穩(wěn)定，能夠準(zhǔn)確地識別出真正的焦散點(diǎn)。同時(shí)，給出了一種合理的雅可比平滑閾值方法，避免了疊前深度保幅偏移中焦散點(diǎn)處的奇異問題。將兩種計(jì)算射線雅可比的方法分別應(yīng)用到逆散射保幅偏移成像，直接數(shù)值計(jì)算獲得了更清晰的成像剖面，進(jìn)一步驗(yàn)證了射線雅可比直接數(shù)值計(jì)算方法的有效性與適用性。