散射波積分方程的Adomian分解解法

汪燚林，董良國

同濟大學海洋地質國家重點實驗室, 上海 200092

0 引言

在目前的地震勘探中，主流反演框架還是基于合成地震數據與觀測地震數據的最佳匹配，最經典的反演方法仍然是基于最優化理論的梯度導引類方法.這類地震反演框架經常要面對的一個關鍵問題是：已知背景模型和背景地震波場，如何更精確地計算一定模型參數擾動下的攝動波場.

地震波散射場與介質參數擾動之間的關系，通常可以通過Green函數的積分方程表達，最為常見的是Lippmann-Schwinger積分方程與Ricatti積分方程.從反演角度看，這些積分方程是非線性的，通過對積分項中的Green函數及總波場進行某些近似來獲得其近似解，其中，Born和Rytov近似應用最為廣泛.這兩種近似在散射波模擬(Devaney, 1981; Snieder and Lomax, 1996; Spetzler and Snieder, 2001, 2004)、層析成像(Woodward, 1992; Liu et al., 2009; Feng et al., 2020b)以及走時和波形反演(Luo and Schuster, 1991; Tarantola, 1984; Pratt et al., 1998)中被大量應用.

然而，對散射場和參數攝動之間實施這樣簡單的線性近似，直接影響了散射波場的模擬精度，必然會降低反演的精度和分辨率，也直接導致了目前地震波形反演方法對初始模型的強烈依賴這個痼疾.在模擬散射場時，如果能夠根據不同情況和需求，利用精度更高的模擬手段，突破傳統的對于散射波場的一階近似，有助于提高復雜介質中散射場的模擬精度，就有可能克服上述成像反演中存在的問題(Innanen, 2009; Albertin et al., 2013; Wu and Zheng, 2014; Jakobsen and Wu, 2016).

要將散射序列中的高階項應用到后續的反演成像中，需要研究散射序列的求解方法，并從正問題的角度認識散射序列中高階項的物理含義和作用.

在現有文獻中，關于Born散射序列和Rytov散射序列的求解，主要是基于Helmholtz方程或者積分方程，通過相應物理量的級數展開而得到的.對非線性Lippmann-Schwinger積分方程進行迭代求解，就可以將非線性的積分方程表達為顯式求和過程的Born散射序列解(Yura et al., 1983; Tsihrintzis and Devaney, 2000a; Osnabrugge et al., 2016; 符力耘, 2010; 董興朋等, 2016), 這種Born散射序列的求解過程相對簡單，但理論支撐不足.基于波場的Rytov變換表達式，結合非均勻介質Helmholtz方程和復相位攝動函數的級數展開，可以得到Rytov散射序列解(Schmeltzer and Robert, 1967; Yura et al., 1983; Tsihrintzis and Devaney, 2000b; Kim and Tinin, 2009),這一類求解方法思路簡單，但推導求解過程比較繁瑣.另外，Marks(2006)從Rytov變換的極限表達式出發并結合復相位函數的級數展開，推導了一種包括了傳統的Born序列和Rytov序列在內的混合序列，通過控制混合參數可以得到不同的散射序列.然而，該方法的推導過程復雜，且并未直接得到Born散射序列和Rytov散射序列.最近，Feng等(2018, 2020a)基于Ricatti非線性積分方程發展了一種在時間域迭代求解Ricatti積分方程的新方法，且得到了一種形式的Rytov序列解.

由Adomian發展起來的分解算法是一種高效求解線性或者非線性積分方程的方法(Adomian, 1990)，它的基本思想就是將積分方程的解表達成一個無窮級數和的形式，并利用一個多項式迭代求解，從低階至高階獨立的求得級數表達式中的各項，進而得到不同精度的逼近解.該方法很好地保留了原方程的物理性質，自該方法提出以來，很多學者對這一算法進行了應用并證明了其解的收斂性(Adomian, 1990; Cherruault et al., 1992; Cherruault and Adomian, 1993; Wazwaz, 1997).同時，針對該方法在計算量、收斂性以及對于復雜問題的處理等方面存在的一些問題，國內外學者也相繼提出了一些改進的方法(Wazwaz, 1999; Wazwaz and El-Sayed, 2001; Duan and Rach, 2011; Zhang and Liang, 2018; 解烈軍, 2013; 郭銀鳳, 2019).總的來說，在面臨一些非線性問題時，Adomian分解方法被認為是一種適應范圍廣泛、計算過程簡單且收斂速度較快的方法(陳一鳴等, 2013; 牛紅玲等, 2013).

在勘探地震領域，研究如何更好地數值求解散射波積分方程、更精確地計算地震波散射場，不但對于認識和理解地震波傳播過程中的散射機理，更好地研究散射波場和介質擾動之間的關系有重要的理論意義，對于提高反演和成像精度也有重要的現實意義.因此，本文將Adomian分解算法應用至散射波積分方程的求解，實現了在統一理論框架下高效求解Lippmann-Schwinger積分方程和Ricatti積分方程，得到了地震波Born散射序列和Rytov散射序列，為下一步發展非線性反演方法提供了精度更高的散射波場正演手段.

1 散射場的積分表達

在空間位置r處介質速度為c(r)的地下模型中，在rs處激發的地震波場u(r;rs)滿足標量波方程

(1)

u(r;rs)=

式中的g(r;r′)為背景介質下的格林函數.

從反演介質參數的角度看，Lippmann-Schwinger積分方程(2)為非線性Fredholm積分方程.但從波場正演的角度看，由于(2)式中的積分函數與u(r′;rs)之間為線性關系，因此Lippmann-Schwinger方程實際上是第一類線性Fredholm積分方程.

設φ(r;rs)為復相位攝動函數，將總波場和背景波場的關系表達為

(3)

則(1)式的解可以表達為下列Ricatti積分方程的形式(Wu, 2003)：

(4)

上述Ricatti積分方程中的積分函數與φ(r;rs)之間為非線性關系.因此，無論是從波場正演還是從地震反演的角度看，Ricatti方程都是均勻第一類非線性Fredholm積分方程，這一點與Lippmann-Schwinger積分方程不同.

通過求解上述Lippmann-Schwinger積分方程，可以確定散射場為Δu=u-u0;通過求解上述Ricatti積分方程，可以確定散射場為Δu=u0(eφ-1).

2 散射波積分方程的Adomian分解解法

2.1 求解非線性積分方程的Adomian分解方法簡介

Adomian分解方法(Adomian decomposition method)是計算數學中求解Fredholm積分方程的常用方法(Adomian, 1990)，下面對該方法進行簡要介紹.

將非線性Fredholm積分方程

(5)

中的u(x)表達為下列分解序列：

(6)

而非線性項F(u(t))是u(t)的非線性函數，將其表達為下列形式的Adomian多項式

(7)

其中，n=0,1,2,….

將(6)和(7)式代入(5)式，有：

(8)

進而可以得到非線性Fredholm積分方程(5)的序列解：

(9)

2.2 利用Adomian分解方法求解Lippmann-Schwinger積分方程

由于線性積分方程是非線性積分方程的一個特例，下面我們把上述Adomian分解方法應用于求解(2)式的Lippmann-Schwinger方程的散射波序列解.

利用(7)式的Adomian多項式，可得：

(10)

將以上各式代入(9)式，即利用Adomian分解方法，可得：

u0(r;rs)=u0(r;rs),

?

(11)

上式中，n滿足n≥0.這樣，就可以得到Lippmann-Schwinger積分方程的解為

+u2(r;rs)+…,

(12)

上述(12)式實際上就是標量波動方程(1)的Born序列解.其中，u0(r;rs)為背景場，而u1,u2,u3,…分別稱為背景場u0(r;rs)的一階、二階、三階以及更高階的校正項，即總的散射場可以表達為

+u3(r;rs)….

(13)

可以看到，一階校正項u1(r;rs)即為常規的Born近似散射場，而u2(r;rs)、u3(r;rs)、…、un(r;rs)則稱為Born散射序列解中的高階校正項.

2.3 利用Adomian分解方法求解Ricatti積分方程

下面利用Adomian分解方法求解(4)式的Ricatti方程的散射波序列解.

(4)式的Ricatti積分方程又可寫為

(14)

利用(7)式可得：

(15)

將以上各式代入(9)式，即利用Adomian分解方法，可以得到非線性Ricatti積分方程的Rytov序列解：

(16)

由上述方程(16)依次計算得到各階校正項φm(r;rs)(m≥0)后，就可以得到包含Rytov散射序列中M階項的散射場：

(17)

其中，M≥1.

當M=1時,(17)式即為常規的Rytov近似散射場Δu(r;rs)=u0(r;rs){exp[φ0(r;rs)]-1}.而當M>1時，(17)式就是對傳統Rytov近似散射場進一步校正的高階Rytov序列逼近的結果.顯然，計算Rytov散射序列中第n階項要用到Rytov序列中第0項至第n-1項，而計算Born散射序列中的第n階項un(r;rs)只需要用到Born散射序列中的第n-1階項un-1(r;rs)，因為un(r;rs)是當前背景模型下的散射源un-1(r′;rs)ε(r′)產生的散射場.

3 散射序列解的數值計算方法

利用波動方程(1)，對真實速度模型和背景速度模型分別進行兩次地震波模擬，兩次模擬結果之差即為該模型擾動下的真實散射場.

理論上，從Born序列解(11)式可以發現，在當前背景模型c0(r)和介質擾動ε(r)的基礎上，分別以ui(r;rs)(i=0,1,2…)作為背景場進行一次地震波模擬，就可以得到散射場的Born序列解中的校正項ui+1(r;rs)(i=0,1,2…).從Rytov序列解(16)式可以發現，在當前背景模型c0(r)和介質擾動ε(r)的基礎上，以φi(r;rs)(i=0,1,2,…)的不同組合作為復合震源進行一次地震波模擬，就可以得到散射場的Rytov序列解中的校正項φi+1(r;rs)(i=0,1,2,…).

實際計算中，結合Born散射序列表達式(11)和Rytov散射序列表達式(16)，可以建立起兩個散射序列解中的各階校正項之間的對應關系：

(18)

在(18)式中，各等式左側為Born序列中的各階校正項(u1,u2,u3,…)與背景波場(u0)，右側為Rytov序列中的各階校正項(φ0,φ1,φ2,…).在本文中，先計算出Born序列解中的各階校正項ui，再根據(18)式，可以求得Rytov序列解對應的各階校正項φi，再分別根據(13)和(17)式便可求得對應不同階的Born序列逼近和Rytov序列逼近下的散射場.

4 數值算例

顯然，通過Adomian分解解法得到的散射波場的Born序列解和Rytov序列解，其中的一階校正項就是傳統的散射波場的Born近似解和Rytov近似解.在滿足一定的條件下，在一階近似解的基礎上，引入Born序列解和Rytov序列解的高階項，可以更精確地描述地震波散射場.

下面利用三個模型試驗，通過和傳統一階近似結果對比，分析散射序列解中的高階項對一階近似散射場的校正作用.

4.1 模型一

首先選取了一個一維模型(圖1)進行測試，震源和接收點分別位于模型的頂端和底端，因此這里僅分析前向散射波場.模型深度1.5 km，真實速度在1900 m·s-1和2100 m·s-1之間變化(圖1中藍線)，背景模型速度為常數2000 m·s-1(圖1中紅虛線).

圖1 一維速度擾動模型與觀測方式(藍線：真實速度；紅虛線：背景速度)模型頂端星號表示炮點位置，模型底端三角形表示檢波點位置.

各階Born、Rytov序列逼近的前向散射波形和真實散射波形對比結果如圖2、圖3所示.由于多次前向散射能量相對較弱，因此將1～2 s的散射波記錄單獨顯示在圖2b和圖3b中.可以看到，無論是對Born散射序列還是Rytov散射序列，不管是能量比較強的一次前向散射還是能量相對較弱的多次前向散射，在一階近似(即Born近似和Rytov近似)的基礎上，添加各級高階項會使得序列解(包括走時和振幅)更加逼近真解，散射場的描述精度得到提高.

另外，從圖2和圖3還可以看出，一階近似幾乎沒有刻畫多次前向散射的能力.考慮高階校正項后，除了可以對一次前向散射進行校正以外(見圖2a和圖3a)，還可以明顯提高對多次前向散射的刻畫能

圖2 各階Born散射序列逼近解與散射波真解對比

圖3 各階Rytov散射序列逼近解與散射波真解對比

力(見圖2b和圖3b).

4.2 模型二

第二個試驗選取了一個二維高斯球速度模型(圖4)進行測試.半徑為150 m的高速異常體位于模型中心，背景為2000 m·s-1的常速模型，異常的速度擾動最大值為300 m·s-1.點震源位于模型上方(1500,100)m處，在距離高速異常體中心半徑為1300 m的右半圓周上大約每10°設置一個檢波器，共記錄19個不同角度的散射波信號，以便觀測不同散射角度上的散射波場.

圖4 高斯球速度模型與觀測系統

由于不同散射角度散射振幅相差很大，為更清晰地對比各個角度不同階散射波形和真實散射波形的匹配情況，對不同角度的散射波進行f′i(t)=fi(t)/max(abs(fi(t)i=1,2,3))的歸一化.其中，i=1,2,3分別表示真實散射場、一階近似散射場和高階散射序列逼近的散射場，fi(t)和f′i(t)分別是同一散射角度歸一化前后的散射波形，max(abs(fi(t)))為該散射角不同散射波數據絕對值的最大值.

圖5、圖8分別顯示了19道不同角度的不同階Born和Rytov序列逼近的散射場和真實散射場的歸一化后的波形對比結果，圖6、圖9分別對應圖5、圖8中高階Born序列和高階Rytov序列逼近解的各道歸一化因子max(abs(fi(t)i=3))的對數隨角度的變化，它體現了不同散射角度的散射場振幅的相對變化.可以發現，在常速背景模型情況下，這個高斯球速度異常在不同的散射角度上產生的散射波場能量差別很大，小角度的前向散射波振幅和大角度的背向散射波振幅相差近4個數量級.小角度的前向散射波能量強、頻帶寬、高頻成分豐富；隨著散射角的增大，散射波能量逐漸變弱、頻帶變窄并向低頻移動.

圖5 不同散射角歸一化后的真實散射場、一階Born近似和高階Born序列逼近散射場

圖6 圖5中高階Born序列逼近解的各道歸一化因子的對數隨角度的變化

圖7 40°～90°真實散射場與不同階Born逼近的散射場波形對比

圖8 不同散射角歸一化后真實散射場、一階Rytov近似和高階Rytov序列逼近的散射場

圖9 圖8中高階Rytov序列逼近解的各道歸一化因子的對數隨角度的變化

為進一步顯示地震散射序列中各高階項對于一階Born和Rytov近似的校正作用，圖7、圖10分別顯示了散射角約40°～90°的不同階Born和Rytov序列逼近的散射波形和真實散射波形的對比結果.可以看到，高斯球速度擾動模型下，對于大部分散射角而言，在一階近似(即Born近似和Rytov近似)的基礎上(紅線)，添加散射序列的二階項后(藍線)的結果與真實散射場更加逼近，再添加散射序列的三階項后(綠線)精度進一步提高.可見，對該模型而言，在許多散射角度上，逐步添加各級高階項計算得到的散射波場(包括走時和振幅)會逐漸逼近真解，散射波場的模擬精度逐漸提高，說明了考慮散射序列中的高階項是很必要的.特別是在一階近似(即Born和Rytov近似)的基礎上，二階項和三階項的作用尤其明顯.

圖10 40°～90°真實散射場與不同階Rytov序列逼近的散射場波形對比

4.3 模型三

第三個試驗選取了一個相對復雜的Sigsbee速度模型(圖11b)進行測試.背景速度為真實模型的平滑模型(圖11a).點震源位于模型上方(2,0.1)km位置處，這里，我們主要關注前向散射波，因此，在模型下方深度為1.4 km位置處間隔10 m共放置了400個檢波器(圖11b)，用于記錄模型的前向散射序列波場.

圖11 (a)Sigsbee平滑背景與(b)真實速度模型及觀測方式

圖12顯示了Sigsbee模型在當前背景速度下的真實前向散射波場記錄，圖13為不同階Born序列解的前向散射波形記錄.當前背景速度條件下，可以發現，一階Born近似所描述的波現象和真實散射波場記錄基本一致，另外，Born散射序列的二階項、三階項對于散射場的刻畫仍然有一定貢獻，且能量上呈逐漸減弱的趨勢.進一步從圖14中的兩個檢波器接收的散射波記錄可以看到，Born散射序列的高階項是對于一階Born近似的校正，不斷添加Born序列的高階項，可以得到更精確的散射波場模擬結果.

圖12 真實前向散射波場記錄

圖13 (a)一階Born近似前向散射波場(u1)；(b)Born序列二階前向散射波場(u2)；(c)Born序列三階前向散射波場(u3)(均在圖12所示的同一色標范圍內顯示)

圖14 在水平位置1.0 km(a)及3.5 km(b)處檢波器接收到的真實散射場與Born序列逼近波場

圖15a顯示了一階Rytov近似前向散射波記錄,圖15b顯示了二階Rytov逼近和一階Rytov近似的前向散射波記錄的差，圖15c顯示了三階Rytov逼近和二階Rytov逼近的前向散射波記錄的差.當前背景速度條件下，和Born散射序列的前三項類似，Rytov散射序列的前三項體現在散射波記錄上的能量逐漸減弱(圖15a—c).另外，進一步從圖16中兩個檢波器接收到散射波記錄可以看到，當前背景速度下，Rytov散射序列的高階項是對于一階Rytov近似的一種校正，不斷添加Rytov散射序列的高階項，可以得到更精確的散射波場模擬結果.

圖15 (a)一階Rytov近似前向散射波場；(b)二階Rytov逼近與一階Rytov近似前向散射波場之差；(c)三階Rytov逼近與二階Rytov逼近的前向散射波場之差(均在如圖12所示的同一色標范圍內顯示)

圖16 在水平位置1.0 km(a)及3.5 km(b)處檢波器接收到的真實散射場與Rytov序列逼近波場

在模型相對復雜，背景模型較好的情況下，逐步添加散射序列的高階項計算得到的散射波場(包括走時和振幅)會逐漸逼近真解，散射波場模擬的精度逐漸提高，說明了在一定條件下，復雜介質情況下考慮散射序列中的高階項是必要的.

5 討論

本文通過Born及Rytov散射序列來表達散射波場，這必然會涉及序列是否收斂的問題，只有序列級數收斂才有實際物理意義.上面三個數值試驗說明：在速度異常體尺度不太大、速度擾動量不太大、模型本身沒有非常復雜以及背景速度模型不太差的條件下，散射場的這兩種序列解還是穩定收斂的，也就是說，通過引入Born序列解和Rytov序列解的高階項，可以更精確地描述真實的地震波散射場.

我們發現，這兩個散射序列級數的收斂性與介質模型的復雜程度、背景介質模型的精度、介質擾動的幾何尺度和擾動強度均有關，也與觀測的空間位置(如不同散射角)有關.譬如，對散射角較大的背向散射和散射角較小的前向散射，這兩個散射序列級數就表現出不同的收斂性質，可見，地震波場散射序列的收斂性是一個十分復雜的問題.Jakobsen和Wu(2016)用強擾動介質的數值實驗揭示了Born散射序列在某些條件下不收斂的現象，但是，文章也沒有定量的給出Born散射序列的收斂條件，另外，文章沒有涉及本文討論的Rytov散射序列.Born和Rytov散射序列在理論上本身存在一定的關聯性，我們發現，兩者在收斂條件方面也具有很大的一致性.然而，影響散射序列收斂性這一問題的因素眾多，不同介質情況下，很難定量給出散射序列能夠對散射場進行精確描述(散射序列收斂)的條件.總的來說，在介質擾動的幾何尺度不太大、擾動量不太大(或背景模型不太差)，散射波強度不太強的條件下，散射序列還是趨于穩定收斂的.因此，在實際應用中，應當根據具體的物理問題，仔細研究散射序列的收斂性質，明確散射序列的實際物理含義，再加以合理運用.

結合前文中數值算例部分的內容，我們在該部分進一步設計實驗，僅考察速度擾動量的大小這一因素對于散射序列收斂性的影響，對散射序列存在的收斂性問題進行舉例論證，以明確散射序列的收斂具有嚴格的收斂條件.

選取類似模型測試二中的二維高斯球速度模型進行對比測試.其他參數與模型測試二相同，只是將圖4中模型的速度最大擾動由300 m·s-1(15%)提高到1600 m·s-1(80%).在這個較差的初始速度模型條件下，圖17顯示了部分散射角的不同階Born散射序列逼近和真實散射波的記錄對比結果，圖18顯示了部分散射角的不同階Rytov散射序列逼近和真實散射波的記錄對比結果.在模型擾動量增大至80%的情況下，首先可以看到，從小角度散射到中角度散射，散射波的能量(振幅)逐漸降低，其次，在散射角相對較大(50°和60°)、散射波能量相對較弱時，此時的一階近似(Born近似及Rytov近似)相對真解引入了較大的近似誤差，但此時引入相應的高階項仍然實現了對于一階近似的校正.然而，在當前介質強擾動情況下，在散射角較小(10°～40°)、散射波能量較強時，引入散射序列中的高階項使得散射序列解更加偏離了真解，說明這時的散射序列不再穩定收斂，散射序列逼近解也失去了相應的物理含義.

圖17 10°～60°真實散射場與不同階Born序列逼近的散射場波形對比

圖18 10°～60°真實散射場與不同階Rytov序列逼近的散射場波形對比

目前，關于散射序列問題的研究主要集中于對正問題的討論，一些學者嘗試通過對傳統散射序列實施重整化(Wu and Zheng, 2014; Jakobsen and Wu, 2016)、數學變換(Osnabrugge et al., 2016; Eftekhar et al., 2018)或者發展新的散射序列(Feng et al., 2018, 2020a; Jakobsen et al., 2020)等思路，以改善散射序列的相關特性.其中，Jakobsen 和 Wu(2016)也是從正問題的角度，考慮介質參數存在強擾動時，設計數值實驗系統的討論了De Wolf序列在收斂特性(精度)方面相較于Born散射序列的優勢，并指出，這種改進的散射序列表達不僅有助于我們理解強散射介質中的波現象，也是進一步發展非線性反演方法的基礎.

反問題方面，目前的地震成像和反演通常基于的Born近似，相當于一階線性近似，從而產生了線性Fréchet導數，這就是目前占統治地位的基于局部線性化的迭代反演，其中，全波形反演(FWI)就是一個典型代表.然而，目前地震波形反演中出現的幾個致命問題，譬如對初始模型的強烈依賴問題，其根源就在于對介質參數擾動和波場攝動之間的關系實施了簡單的線性近似(Born近似或Rytov近似).實際地震波場與介質參數之間是強烈的非線性關系，如果能夠很好地利用散射序列中的高階項，就有可能規避目前成像和反演中存在的諸多問題.

6 結論

地震反演的精度很大程度上取決于正問題的描述精度，建立更精確的模型擾動和波場擾動之間的關系十分關鍵.目前，基于Lippmann-Schwinger積分方程和Ricatti積分方程近似下建立的散射場的一階近似方法直接影響了后續成像和反演的精度.本文將計算數學領域中的Adomian分解方法應用至求解散射場的Lippmann-Schwinger和Ricatti積分方程中，簡潔、高效地得到了散射場積分方程的序列解.數值試驗結果表明，如果采用經典的一階Born或者Rytov近似，某些情況下難以保證近似散射場描述的精度.在一定的條件下，散射序列解收斂時，與傳統的Born和Rytov近似解相比，引入散射序列中的高階項可以更精確地描述地震波散射場.