射線雅可比的直接數值計算方法及應用

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

(①中國科學院測量與地球物理研究所計算與勘探地球物理研究中心，湖北武漢 430077；②大地測量與地球動力學國家重點實驗室，湖北武漢430077； ③中國科學院大學，北京100049；④東方地球物理公司物探技術研究中心，河北涿州 072751)

0 引言

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

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

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

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

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

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

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

1 射線雅可比的定義與計算

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

1.1 動力學射線追蹤

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

(1)

(2)

式中m和n表示與中心射線垂直的平面坐標。

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

(3)

進而可得式(2)的矩陣V中各個二階偏導數的具體表達式為

(4)

(5)

此時

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

1.2 射線雅可比的直接數值計算

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

(6)

(7)

J=Ω(u)·e

(8)

于是，可根據

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

(9)

計算關于任意參數u的雅可比J(u)(u=T，s)。

根據矢量微分可得

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

(10)

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

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

(11)

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

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

(12)

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

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

(13)

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

圖1 初始角面元計算示意圖

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

2 兩種射線雅可比計算方法的比較

射線雅可比既可為正，也可為負，還可以為零。雅可比為零表示射線族產生交叉，局部波前面匯聚成一點，這樣的點稱為焦散點。對相同的射線，分別應用動力學射線追蹤方法和直接數值計算方法計算射線雅可比，通過得出的雅可比曲線及其射線族展開，比較兩種方法的準確性。

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

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

圖2 速度分布及對應的4條射線(標號為①～④)

圖3 對應圖2中的4條射線兩種方法計算的雅可比曲線

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

圖5 對應圖4中的4條射線兩種方法計算的雅可比曲線

由以上測試結果可知，對于介質速度擾動較小的情況，兩種射線雅可比計算方法獲得的雅可比曲線能很好地符合。當介質速度變化較大時，速度場的二階偏導數不穩定，動力學射線追蹤計算的射線雅可比難以保證準確性；這時基于物理定義的射線雅可比直接數值計算方法無需計算速度場的二階偏導數，求得的雅可比具有更好的穩定性和準確性。

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

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

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

表1 二維介質中兩種射線雅可比計算方法的耗時

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

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

表2 三維介質中兩種射線雅可比計算方法的耗時

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

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

3 逆散射保幅偏移成像

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

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

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

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

(14)

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

(15)

(16)

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

(17)

其中

(18)

式中：v0(s)與v0(r)分別為炮點與檢波點的波速； dt為射線追蹤的時間采樣間隔；Jdt(s)和Jdt(r)分別為關于炮點與檢波點的基準雅可比；ε為雅可比閾值參數(通常設為0.1%～1.0%)；J′為平滑閾值處理后的雅可比。以J′代替J，則式(16)變為

(19)

式(19)能夠克服在焦散點附近過小的J引起的奇異問題。另外，當J較大時，J′→J，保證了非焦散區域計算結果的準確性。

3.2 偏移成像數值試驗

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

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

4 結論

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