廣義Rytov近似有限頻初至層析

馮 波 WU Ru-Shan 羅 飛 許榮偉 王華忠

(①同濟大學海洋與地球科學學院，上海 200092；②同濟大學海洋高等研究院，上海 200092；③Modeling and Imaging Laboratory，University of California，Santa Cruz，CA 95060，USA)

0 引言

地震波旅行時常用于反演地下速度結構。在射線層析中，由于引入了高頻近似假設，旅行時只與連接炮檢對的幾何射線上的速度擾動有關。對于有限頻帶的地震波，若速度非均勻體的尺度小于菲涅耳體的寬度時，地震波散射起主導作用。由此發展了繞射層析理論[1-2]，可以較好地處理地震波的一階繞射效應。基于Born近似， Luo等[3]給出了地震波旅行時Fréchet導數(旅行時敏感度核函數)，可以定量描述速度擾動對地震波旅行時擾動的影響。通過波動方程線性化近似(Born近似或Rytov近似)，Woodward[4]提出了“波路徑”的概念，作為無限高頻假設下的射線路徑向有限頻地震波傳播的推廣。借助伴隨狀態理論，可以實現(波形、旅行時、振幅等)Fréchet導數的高效計算[5-11]。

通常旅行時Fréchet導數由Born近似導出[12-18]，但Born近似成立條件較為苛刻。相比于Born近似， Rytov近似可以更好地描述由于速度擾動引起的前向散射波的相位擾動[19]，因而更適用于旅行時反演。基于Rytov近似也可以構造相應的旅行時敏感度核函數[20-27]。

然而，Born近似與Rytov近似均屬于弱散射近似，要求速度擾動遠小于背景速度[28-29]。為克服弱散射近似的制約， Feng等[30]提出了可以適用于強擾動模型的廣義Rytov近似理論(Generalized Rytov approximation，GRA)。對于前向散射及小角度傳播，GRA可以較為準確地預測地震波的相位擾動，因而更適用于旅行時反演。本文將GRA理論應用于傳統的有限頻層析，構造一種適用范圍更廣的旅行時敏感度核函數并用于初至旅行時層析。模型數據測試結果表明，本文提出的廣義Rytov近似旅行時層析可以建立高精度的近地表速度模型，且收斂速度較快。

1 理論

1.1 廣義Rytov近似有限頻旅行時敏感度核函數

頻率域標量聲波方程可以表示為

(2+k2)u0=0

(1)

式中：k(x)=ω/v0(x)，為背景速度場v0(x)中的地震波數；u0(x，ω)為背景波場；2為Laplace算子。

記擾動后的速度場為v(x)，波場(總場)為u(x，ω)，總場與背景場的關系為

u(x,ω)=u0(x,ω)exp[ψ(x,ω)]

(2)

式中ψ(x，ω)為復相位。

在前向散射(散射角定義如圖1所示)及小角度傳播假設下，復相位可以用廣義Rytov近似[31]表示為

圖1 散射角定義：地震波入射方向與散射方向的夾角前向散射對應小散射角，背向散射對應大散射角

(3)

式(3)引起的相位延遲可以表示為

(4)

有限頻帶地震信號的旅行時擾動可以表示為單頻諧波的相位延遲的加權疊加[23-26，32]，有

(5)

結合式(3)～式(5)，可以建立帶限信號旅行時擾動與慢度擾動的線性關系

(6)

式中：xs和xr分別為源和檢波點坐標；KGRA(x;xr,xs)為廣義Rytov近似旅行時敏感度核函數(Generalized Rytov approximation based traveltime sensitivity kernel， 簡稱為GRA核函數)，可表示為

(7)

顯然，顯式計算式(7)需要計算頻率域格林函數[25]，因而計算量較大，尤其是對于三維問題。為了避免顯式計算核函數，本文采用馮波等[26]給出的方法，通過引入如下輔助函數

(8)

利用Parseval定理，將式(7)中對角頻率的積分轉化為對時間積分，則

KGRA(x;xr,xs)

G0(x,ω;xr)]dω}

(9)

式中：T表示地震記錄時長；g0表示時間域格林函數；“*”表示褶積(與上標“*”表示復共軛不同)。

若定義旅行時伴隨場為

(10)

式(9)可以重寫為

KGRA(x;xr,xs)=

(11)

1.2 旅行時層析反問題求解

基于旅行時殘差L2范數的誤差泛函可表示為

(12)

式中：m=(m1,m2,…,mNv)T∈M，為模型空間M中的向量，本文為慢度模型(采用網格參數化方式)，其中Nv為網格點數；Δt=(Δt1,Δt2,…,ΔtNd)T∈D，為數據空間D中的向量(旅行時殘差)，其中Nd為炮檢對數。

利用Gauss-Newton反演算法，式(12)可以迭代求解

(13)

式中：k為迭代次數；α為步長，可以通過線性搜索方法計算；J為目標函數梯度；為線性Hessian矩陣(忽略旅行時對模型參數的二階導數)，其中K為KGRA的矩陣形式；P為光滑算子。

根據式(12)，泛函梯度可表示為

(14)

基于廣義Rytov近似，?Δt/?m即為式(11)給出的核函數。將式(11)代入式(14)，則泛函數梯度可表示為

J(mk)(x)=(KTΔtk)(x)

(15)

考慮到顯式計算Hessian矩陣需要巨大的計算量及存儲量，且求逆困難，因此，將對Hessian矩陣的求逆轉化為求解如下法方程的近似解

Ha(mk)Δmk=J(mk)

(16)

本文用共軛梯度(Conjugate gradient，CG)法求解式(16)獲得模型更新方向，并采用隱式矩陣向量乘方法[26]直接計算Hessian矩陣與向量的乘，因而無需顯式計算及存儲Hessian矩陣。具體計算公式詳見附錄A。

2 數值試驗

2.1 高斯擾動模型

圖2為一個含有高斯異常體的二維10km×5km速度模型，網格間距為10m。其背景速度v0=3000m/s，異常體的速度擾動可表示為

圖2 含有高斯擾動(ε= 50%)的速度模型及觀測系統紅色倒三角為檢波點位置

(17)

首先，用時—空域高階有限差分求解聲波方程，用真實模型(包含高斯異常體的速度模型)和背景模型(即均勻背景速度模型)進行正演(震源是主頻為15Hz的Ricker子波)，得到的地震記錄分別作為觀測數據和模擬數據。然后，用互相關(Cross-correlation，簡記為CC)函數測量觀測數據和模擬數據的時差，作為真實時差的測量值(即觀測時差)。隨后，分別用射線(Ray)核函數(速度網格內的射線長度)、Rytov近似旅行時敏感度核函數[26]及本文提出的GRA旅行時敏感度核函數預測旅行時擾動(即核函數與模型擾動做內積)。圖3為震源在地表(5.0km，0)和(0，0)處、201個檢波器均勻分布在z=5.0km處(圖2中紅色倒三角所示)旅行時擾動對比。由圖3可以看出：①射線和GRA核函數的預測時差完全重合，但與Rytov近似預測結果有較大偏差。②對于小角度前向散射(圖3a中對應為：震源位于地表(5.0km，0)，檢波器橫坐標在(5.0km，5.0km)附近)，GRA核函數的預測時差與觀測時差基本一致。證明了GRA核函數在前向小角度散射時具有較好的線性特征，克服了傳統波動方程線性化近似中的速度小擾動假設的制約。這與GRA理論預測相符(GRA理論指出，對于小角度前向散射，GRA可以較為精確地預測地震波的相位擾動)。③隨著散射角增大，GRA核函數的預測時差逐漸偏離觀測時差，并逐漸與Rytov近似預測結果趨于一致。根據GRA理論推導可知[30]，Rytov近似為GRA的一階近似，本質上都隱含了小角傳播假設。因此，隨著散射角增大，GRA的成立條件(即前向小角散射)不再成立，導致預測時差偏離觀測值。對于圖3b中的震源位置(左上角(0，0)處)，前向小角度散射對應檢波器橫坐標在10km附近，三類核函數預測結果與圖3a一致。

圖3 震源在地表不同位置三種旅行時敏感度核函數的預測時差與測量時差對比(a)(xs，zs)=(5.0km，0)；(b)(xs，zs)=(0，0)

上述測試表明：GRA核函數精度與散射角密切相關。對于前向小角度散射，GRA可以較為準確地預測地震波的旅行時擾動。隨著散射角不斷增大，GRA核函數精度逐漸降低，但仍好于Rytov近似。雖然對于大角度散射，GRA核函數的預測結果并不準確，甚至與觀測時差偏離較遠，但并不會導致層析反問題求解失敗。原因在于：層析反問題通過最優化算法迭代求解時差目標函數，而非直接對線性方程組直接求逆。即使正問題的線性化近似不嚴格成立，只要反問題的凸性較好，依然可以采用最優化算法求解目標泛函。

為了驗證GRA初至層析方法的有效性、盡可能消除觀測系統對反演結果的影響，設計了一個理想觀測系統如圖4a所示，100個檢波器均勻分布高斯異常體模型的四周，圍繞模型四周激發100次，震源是主頻為15Hz的Ricker子波(對應波長λ0=0.2km，與高斯異常體尺度參數關系為a= 5λ0)。

圖4 高斯模型觀測系統(a)及GRA旅行時反演結果(b)黑色正三角形表示震源位置；紅色倒三角形表示檢波點位置

反演采用的初始模型為均勻背景速度場(v0=3000m/s)，并用互相關方法測量觀測數據和模擬數據的時移量作為初至旅行時殘差。反演采用SEIS-COPE數值優化軟件包(SEISCOPE optimization toolbox)中的截斷牛頓算法(對于旅行時層析，截斷牛頓算法退化為高斯牛頓算法)實現最優化，其中，梯度和Hessian矩陣與向量乘由附錄A中給出的隱式矩陣與向量乘方法計算。在迭代反演中，采用2次CG內迭代求解式(16)，獲得模型更新方向，并用線性搜索方法計算迭代步長。反演收斂準則為判斷相對旅行時誤差(J(mk)/J(m0))是否小于預設的閾值(σ=1.0×10-4)。經過10次迭代(圖5為反演收斂曲線)，反演結果如圖4b所示，與真實速度模型在視覺上幾乎沒有差別。為定量對比反演結果，對反演的擾動模型(反演結果減去初始模型)與真實的高斯異常體進行橫向和縱向抽線對比(圖6)，反演的速度擾動與真實的高斯擾動基本吻合，證明了對于完備的觀測系統及簡單高斯速度模型，GRA初至層析可以得到高精度的反演結果，證實了方法的有效性。

圖5 GRA初至層析的收斂曲線

圖6 高斯模型層析(藍色)與真實(紅色)速度擾動曲線(a)x=2.5km；(b)z=2.5km

2.2 Marmousi-2模型

應用Marmousi-2 模型驗證GRA初至層析對復雜速度模型及不完備的觀測系統的有效性。為模擬陸上地震采集，切除了原始速度模型中上覆的水體部分，并在模型右側擴邊3km(圖7a)。x和z方向網格數目分別為2001和301，網格間距均為10m。

圖7 Marmousi-2模型(a)及GRA初至層析的速度更新量(b)

設計了一個單邊觀測系統：震源和檢波器均放置在地表，起始炮點位于x=0，炮間距為100m，共171炮；最小和最大炮檢距分別為400m和5000m，道間距為20m(每炮共231道)。最大炮檢距較大，因此可以用初至波旅行時反演近地表速度。采用有限差分求解波動方程正演地震信號，震源是主頻為15Hz的Ricker子波。由于觀測地震記錄較為復雜，采用自動方法[31]拾取初至，作為觀測記錄的初至旅行時。

初始模型采用隨深度線性遞增的速度場(z=0，v=1000m/s；z=3.5km，v=3000m/s)。為降低計算量，采用20m×20m網格剖分進行速度反演，并用相同的初至拾取方法拾取模擬數據的旅行時，進而計算初至時差。反演實現框架與高斯模型一致(相對旅行時誤差閾值σ=1.0×10-2)。經過5次迭代反演收斂，得到的速度更新量如圖7b所示，淺層的高速異常得到了較好的恢復。其收斂曲線如圖8所示，可以看出，高斯—牛頓算法收斂速度較快。

圖8 GRA初至層析的收斂曲線

為定量對比反演結果，抽取z=0.1、0.2、0.4、0.6km處速度橫向變化曲線，如圖9所示。可以看出，在近地表(z=0.1、0.2km)，反演結果(藍色曲線)逼近真實的高速隆起，反演精度較高(圖9a和圖9b)。隨著深度增加，反演精度逐漸降低(圖9c和圖9d)，但反演結果仍然與真實速度模型的整體趨勢相吻合。

圖9 不同深度處初始(綠色)、真實(紅色)、GRA層析(藍色)速度橫向變化曲線對比(a)z=0.1km；(b)z=0.2km；(c)z=0.4km；(d)z=0.6km

Marmousi-2模型實驗結果表明，利用中、小炮檢距的初至波旅行時信息，GRA初至層析可以建立高精度的近地表速度模型。

3 討論

旅行時(差)的精確測量是所有旅行時反演類方法所共有的問題。本文的重點為討論如何將廣義Rytov近似方法與有限頻層析理論的結合，是筆者過去研究工作[26]的理論推廣和應用，因而沒有討論如何準確測量初至波時差。若震源子波已知且模擬數據與觀測數據的初至波形基本一致，此時測量初至波時差相對簡單。反之，若震源子波未知或由于地下介質復雜導致初至波形畸變嚴重，此時需要采用與震源無關的反演類算法消除子波波形畸變對旅行時測量的影響，如雙差旅行時(double-difference traveltime)[33-34]測量等。

4 結論

本文提出了一種基于廣義Rytov近似的有限頻旅行時敏感度核函數，克服了傳統有限頻層析方法中弱速度擾動假設的制約。對于前向散射及小角度傳播，廣義Rytov近似旅行時敏感度核函數可以較為準確地預測地震波的旅行時擾動，從而可以提高反演精度并加速反演收斂。在計算方面，本文采用了隱式矩陣向量乘方法直接計算Hessian矩陣與向量的乘，無需顯式計算及存儲Hessian矩陣，進而可以提高高斯—牛頓迭代反演算法的效率。復雜的Marmousi-2模型測試結果表明，GRA初至層析可以建立高精度的近地表速度模型。

附錄A 隱式矩陣向量乘

Hessian矩陣與向量積可表示為

Hap=KTKp=KTh

(A-1)

式中：p=p(x)為模型空間中的向量；h=h(xr,xs)為數據空間中的向量。因此需要計算兩次矩陣向量積。

首先， ?(xr,xs)，Kp可以表示為

[G0(x,ω;xr)q*(xr,ω;xs)]dω}dx

(A-2)

式中up(xr,ω;xs)滿足如下Born近似

(A-3)

利用Parseval定理，式(A-2)可改寫為

(A-4)

式中

uq(xr,t;xs)=2πRe[q*(xr,t;xs)]

=2πRe[FT-1q(xr,ω;xs)]*

(A-5)

up(xr,t;xs)=FT-1up(xr,ω;xs)

(A-6)

同時，式(A-3)在時—空域滿足如下方程組

(A-7)

式中f為震源函數。

其次，根據式(15)，對于空間任意一點x，KTh可以表示為

(A-8)

式中

(A-9)

為單炮旅行時伴隨波場。