基于聲波散射理論的零井源距VSP波場數值模擬方法

顏 紅 芹

(婁底職業技術學院，湖南婁底 417000)

0 引言

基于VSP數據的波形反演[1]的實質是，在一定條件約束下，通過多次迭代，最小化模型參數的正演模擬合成數據與觀測數據的殘差[2-4]，因而波動方程正演模擬的數值計算效率和計算精度是制約全波形反演應用于實際的重要因素[5]。針對高精度的VSP正演數值模擬問題，國內外學者展開了一系列研究。鄒延延等[6]比較完整地總結了VSP正演方法，包括射線追蹤法、高斯射線束法[7]、反射率法[8]和有限差分法。其中，有限差分法[9-10]由于計算效率高而成為VSP波形反演中最常用的正演模擬方法，但該方法受限于頻散條件及穩定性條件，需要嚴格的子波主頻和網格比匹配。另外，人工邊界[11]吸收不完全也會引入多余的干擾。

基于體積分方程的正演模擬方法則很好地克服了以上有限差分算法存在的問題，不受正演子波頻率的限制，且正演過程不會引入人工邊界導致的干擾。Fokkema等[12]介紹了聲波互易性原理在地震勘探中的應用，在均勻介質假設下，把二維一階應力—速度聲波方程組轉化為邊界積分和體積積分形式，將該正演問題轉化為求解第二類Fredholm積分方程，至此拉開了研究波動方程積分數值解法的序幕。Lam等[13]將WKBJ近似方法引入背景介質格林函數的求解中，利用Neumann級數迭代方法求解一維聲波方程的近似解，但是數值計算方法忽略了地下介質和波場之間嚴格的非線性特征，并且計算效率較低。Yang等[14]提出了一種共軛梯度—傅里葉變換(CG-FFT)法求解三維彈性波動方程散射問題，該方法在高頻情況下仍具有較高的計算效率和計算精度，不足之處在于該方法選定均勻介質作為背景介質，模型較復雜時，格林函數不能滿足計算精度要求。Haffinger等[15]在求解聲波方程散射問題時采用了一種介于Neumann級數迭代法和共軛梯度迭代方法之間的一種數值求解方法，計算過程比較復雜。

針對體積分聲波方程，在前人的研究基礎上，本文提出了一種提高積分方程計算效率的數值計算方法，求解零井源距VSP波場。首先將積分方程轉化為矩陣相乘的形式，給出了在常背景介質假設下的格林函數以及正演算子的矩陣表達格式，并給出了均勻變化的背景介質假設下的格林函數表達式；通過預條件最小二乘方法實現了波場的直接快速精確求解；最后，以階梯狀模型和根據實際聲波測井數據建立的模型驗證本文數值計算方法的可行性和有效性。

1 理論及方法

本文從三個方面闡述基于散射理論的聲波方程正演問題的快速、穩定數值求解方法。首先，詳細介紹了由散射理論將聲波方程微分方程轉換為Lippman-Schwinger積分方程的過程以及求解該積分方程的迭代計算方法；然后，將該第二類Fredholm積分方程通過簡化的正演算子轉化為大型線性方程組的求解問題；最后，介紹了求解該問題的預條件最小二乘數值計算方法。

1.1 聲波方程散射問題

在常密度假設下，頻率—空間域聲波在地下各向同性介質中傳播滿足Helmholtz方程

(1)

式中：k(x)=ω/c(x)為波數，其中c(x)為介質的真實縱波速度場；u(x，xs,ω)為頻率—空間域的聲波全波場；xs(xs，ys，zs)為源點坐標；f(ω)為震源子波，ω=2πf為角頻率。

根據聲波散射理論，全波場u(x，xs，ω)可描述為在背景介質中傳播的入射波場u0(x，xs，ω)和由速度擾動產生的散射波場us(x，xs，ω)的疊加，即

u(x，xs，ω)=u0(x，xs，ω)+us(x，xs，ω)

(2)

入射波場u0滿足

(3)

式中k0(x)=ω/c0(x)，其中c0(x)為背景速度，一般通過對真實縱波速度場c(x)做低通濾波獲得。引入背景介質對應的格林函數G(x，xs，ω)，可得入射波場的解析解為

u0(x，xs，ω)=f(ω)G(x，xs，ω)

(4)

引入對比函數

(5)

根據式(1)～式(3)，散射波場us(x，xs，ω)滿足

(6)

在Sommerfeld輻射邊界條件假設下，利用互易定理，式(6)所示的聲波散射波場的微分方程等價于如下的積分方程

(7)

上式也稱數據方程，描述了散射場與對比函數和全波場之間的非線性關系。式(7)兩邊同時加上入射波場u0(x，xs，ω)，則全波場可表示為

u(x，xs，ω)=u0(x，xs，ω)+

(8)

式(8)即為求解Helmholtz方程空間全波場的Lippman-Schwinger積分方程，對應于數據方程，被稱為目標方程，描述了全波場與對比函數之間的非線性關系，被廣泛應用于積分方程對比源反演[16-17]。從物理意義上看，通過散射理論，該積分方程將聲波全波場分解為在背景介質中傳播的入射波場以及將真實介質中異于背景介質的散射體看作二次震源所產生的散射波場之和；從數學形式上看，該積分方程是典型的第二類Fredholm積分方程，格林函數和對比函數的乘積構成了積分核，通常采用迭代方法[18]求解。

1.2 數值迭代求解方法

簡單起見，將式(8)寫為

(9)

(1)給定迭代次數，迭代誤差限；

(2)選定入射波場u0為全波場u的零級近似；

1.3 最小二乘數值解法及對角預條件算子

為實現基于聲波體積分方程的零井源距VSP波場快速數值模擬方法，將式(9)寫為矩陣形式

Au=g

(10)

式中

(11)

I為單位矩陣；A為線性方程組的系數矩陣；

g=u0

(12)

假設震源點坐標為zs，根據留數定理[21]可得常速度背景介質假設下一維頻率域格林函數G的解析解

(13)

將空間網格離散為z1，z2，…，zn，可得格林函數的離散格式

(14)

(15)

另外，數值積分具有可以對網格進行非等間距剖分的特點，因此可以根據速度的復雜情況更加合理地進行網格剖分，因此本文方法能更好地適應較復雜的速度模型。

常速度背景介質假設往往會導致波場旅行時的畸變。在均勻變化的速度背景介質假設下，基于WKBJ近似，格林函數G解析解為

G(z,zs,ω)

(16)

類似于常速度背景介質假設，可以得到均勻變化的速度背景介質下的離散形式格林函數G以及正演算子L。

至此，將式(8)的聲波方程散射問題轉化為如式(10)的大型線性方程組的求解問題。

由于離散的格林函數矩陣G是對稱正定的，因此，考慮到系數矩陣A的對稱性和正定性，引入穩定高效的最小二乘數值解法求解式(10)，實現全波場的直接求解。最小二乘法通過最小化誤差的平方和，尋找解的最佳函數匹配。引入目標函數

(17)

u=A-1g或u=(ATA)-1ATg

(18)

由上式可知，進行聲波方程正演時，在每個頻率點處，基于數值積分的方法只需做一次矩陣的求逆運算。相比于Neumann級數迭代法，上述求解方法比較簡潔，數值測試表明本文方法計算效率較高。

另外，當頻率較高時，格林函數矩陣G的條件數往往較大，會造成數值計算的不穩定[23]。引入對角矩陣D=diag(A)作為預條件算子，求解DAu=Dg，而不是直接求解Au=g。則式(10)的最小二乘解為

u=(DA)-1Dg

(19)

2 數值分析

2.1 階梯狀模型

為了驗證本文方法的有效性、可行性以及解的精確性，設計了一個如圖1所示的階梯狀速度模型。模型的深度采樣間隔為5m，深度采樣點數為301，正演采用主頻為35Hz的零相位Ricker子波。記錄長度為2.5s，時間采樣間隔為2ms。在頻率域，由背景速度計算格林函數G以及正演算子L。對0～250Hz的頻率范圍采用0.4Hz的頻率間隔，分別利用上述預條件最小二乘方法和Neumann級數迭代方法求解，進行626次正演，得到頻率域的全波場u，對其進行傅里葉反變換得到時間域的VSP全波場記錄如圖2a所示(兩種方法的計算結果相差無幾，僅展示一個)。圖中可見清晰的直達波和反射波，其后可見明顯的多次波，說明該方法可以有效地模擬包括直達波、一次反射波以及多次波在內的完整的全波場信息。在同等軟硬件配置條件下，表1列出了本文的數值積分方法、Neumann級數迭代方法計算該階梯狀模型波場時的運算時間，可見本文的數值積分法在計算聲波積分方程時具有相對較高的計算效率。

圖1 層狀速度模型

為了說明本文方法計算精度的優勢，與時間2階、空間2階差分法進行對比。采用10層完全匹配層吸收邊界和相同的模擬參數，時間2階、空間2階交錯網格有差分法正演的零井源距VSP全波場記錄如圖2b所示。與圖2a相比，二者的旅行時信息和波形信息基本一致。但是，聲波方程交錯網格有限差分方法需要滿足數值計算的穩定性條件

(20)

式中：vmax為最大速度；am為2M階空間差分算子系數[24]。

在相同的網格參數情況下，本文方法的計算結果不受數值頻散以及人工邊界反射的干擾，而有限差分方法正演受到選取的正演子波主頻和網格參數的限制，計算的波場中存在嚴重的數值頻散假象

圖2 層狀模型不同方法正演的零井源距VSP記錄(a)本文方法； (b)時間2階、空間2階有限差分法； (c)時間2階、空間8階有限差分法

(圖2b中矩形黃框所示)。另外，由于PML吸收邊界不能完全吸收到達邊界的波場，因此波場中存在人工邊界反射(圖2b中矩形紅框所示); 而吸收效果更好的邊界條件往往會導致計算量的增加[11]。對于圖1速度模型，保持空間采樣間隔和時間總記錄長度不變，減小時間采樣間隔為0.5ms，采用300層完全匹配層吸收邊界，時間2階、空間8階交錯網格有限差分法正演結果如圖2c所示，可以較好地壓制數值頻散和邊界反射。從VSP記錄0.25～1.00s、0～250m范圍的放大顯示(圖2中小圖)可見，圖2b中存在邊界反射干擾和數值頻散假象，而圖2a中則不見這些干擾，圖2c中的干擾得到了有效壓制。

抽取圖2切除直達波后的第一個接收點的記錄，如圖3所示，有5個(4個正相和1個負相)明顯的反射波形。可以看出：三種方法正演結果的波峰到達時刻一致，但時間2階、空間2階有限差分法(紅色實線)受到空間采樣間隔和子波頻率的限制，在每一個反射界面處，由于數值頻散產生能量較強的續至波形；時間2階、空間8階有限差分方法(藍色虛線)的結果能夠有效壓制頻散現象。作為聲波微分方程的解析解，對式(9)進行數值求解比用差分法對微分方程直接求解精度更高，但當空間采樣點數相對較多時，數值計算相對耗時，可以根據不同的需求選用不同的數值求解方法。

表1 不同方法的計算時間對比

圖3 層狀模型不同方法正演結果的單道對比

黑色實線、紅色實線和藍色虛線分別為本文方法、聲波方程時間2階空間2階有限差分法和時間2階空間8階有限差分法正演結果不同頻率下格林函數矩陣G的條件數如圖4a所示，可見當正演頻率較高時，格林函數矩陣G的條件數會存在一些異常大的值以及振蕩。圖4b展示了矩陣A(紅色虛線)和DA(藍色實線)條件數的對數及二者之差的對數(黑色實線)隨頻率的變化情況。由圖可知，在150Hz范圍以內，隨著頻率的增加，條件數逐漸增大；在頻率范圍150～200Hz條件數變化比較劇烈；頻率超過200Hz時，條件數逐漸增加，同時，二者條件數之差也有相似的規律。值得注意的是，頻率較高時，矩陣DA比矩陣A的條件數大約減小了10%，亦即在理論上，引入預條件算子D小幅度提高了計算的穩定性。另外，在每個頻率下，矩陣DA的條件數都大于0，根據條件數的定義

K(A)‖A‖‖A-1‖

(21)

可知矩陣DA是滿秩的。由矩陣論可知，應用D=Diag(A)作為預條件算子后，可以較好地保證解的一致性。

圖4 不同頻率下三個矩陣的條件數a)格林函數矩陣G； (b)系數矩陣A(紅色虛線)和系數矩陣DA(藍色實線)及其差(黑色實線)

2.2 聲波測井數據模型

實際聲波測井資料(圖5a)采樣間隔為0.125m，井位處的地層傾角接近0°。應用Backus平均方法[25-26]獲得間隔為4、10m重采樣速度曲線(圖5b和圖5c)。

為了測試不同深度采樣間隔對地震響應的影響，利用本文提出的預條件最小二乘方法求解兩種采樣間隔模型對應的波場。設置地震記錄的時間長度為0.5s，采樣間隔為2ms；正演采用主頻為35Hz的零相位Ricker子波。在頻率域，對0～250Hz的頻率范圍采用2Hz的頻率間隔，進行126次正演，得到頻率域的全波場，對其進行傅里葉反變換，得到時間域的零井源距VSP波場記錄(圖6a和6b)。如果采用如上的網格參數，則有限差分法計算不穩定。減小時間采樣間隔到0.5ms，記錄長度為0.5s，采用50層的完全匹配層吸收邊界，時間2階、空間8階有限差分法模擬的零井源距VSP波場記錄如圖6c和6d所示。由圖6可見，如采用適當的有限差分網格參數以及適當層數的PML吸收邊界，兩種方法計算的波場基本一致，都很好地模擬了聲波在實際介質中的傳播，但較小深度采樣間隔的速度模型對應的波場細節更加豐富。表2展示了不同方法實現數值模擬波場時所需要的計算時間，可見，當空間采樣點數比較小時，本文的數值積分法和有限差分法的計算效率非常接近。同時，本文的數值積分法的理論計算效率與網格參數的關系還有待進一步深入地研究。

圖5 實際聲波測井(a)及其4m(b)和10m(c)間隔重采樣曲線

圖6 實際測井數據速度模型兩種方法計算的零井源距VSP記錄(a)采樣間隔為4m，本文方法； (b)采樣間隔為10m，本文方法； (c)采樣間隔為4m，有限差分法； (d)采樣間隔為10m，有限差分法

表2 不同空間采用間隔兩種方法模擬用時

切除直達波之后，從圖6a與圖6b中分別抽取第一個接收點的記錄并歸一化，其波形、頻譜分別如圖7所示。由圖可見，當空間采樣間隔為10m時，在地震有效頻段內存在嚴重的陷頻現象，并且缺少有效的中高頻段信息；當空間采樣間隔為4m時，在有效地震頻段內，該頻譜在中高頻率段的信息更加完整和豐富。因此，對于相同的測井速度，重采樣間隔較小的速度模型對應的地震記錄具有更強的分辨能力。

圖7 空間采樣間隔為10m(黑線)和4m(紅線)本文方法模擬記錄的單道波形(a)及其頻譜(b)對比

3 結論

基于聲波散射理論，本文提出了一種精確求解聲波方程得到零井源距VSP波場記錄的數值計算方法。將聲波微分方程推導為Lippman-Schwinger積分方程，再基于矩形積分公式，將其改寫為矩陣相乘的形式，將Fredholm第二類積分方程的聲波散射正問題轉化為凸函數的最優化問題，借助預條件最小二乘方法實現了聲波波場的穩定求解。

本文方法將Sommerfeld輻射邊界條件引入積分公式推導中，因此不會在零井源距VSP波場中引入邊界反射。同時不受子波頻率和空間采樣間隔的影響而產生數值頻散，為后續進行高精度零井源距VSP波形反演奠定了較好的基礎。

通過階梯狀速度模型和實際聲波測井數據建立的速度模型，驗證了方法的可行性和有效性。模型正演結果表明，本文提出的方法克服了有限差分方法求解聲波微分方程時受限于頻散條件和穩定性條件，以及人工邊界吸收不完全的固有缺陷，具有明顯的優勢。

另外，其他可以增強計算穩定性或者提高解的精度以及計算效率的預條件算子值得進一步深入研究。將本文的方法推廣到二維聲波方程和彈性波方程，模擬非零井源距VSP波場，也是下一步的研究方向。