一種利用Nystr?m離散與FFT快速褶積的散射地震波并行計算方法

徐楊楊， 孫建國 ， 商耀達

吉林大學地球探測科學與技術學院， 長春 130026

0 引言

類似于光學散射和電磁散射問題的數值模擬正演方法，常規的聲波散射數值模擬方法也分為兩大類：基于微分方程理論的有限差分法(Finite Difference)、有限元法(Finite Element)等和基于積分方程理論的邊界元法、矩量法、Nystr?m法等(Fu et al.,1997;Fu and Wu,2000; Hsiao et al.,2000；劉立彬等，2020；邱長凱等，2020；徐楊楊等，2021).微分方程類的數值算法具有容易實現、計算速度快等優點，但是網格剖分需要在全空間上進行，對于較小的異常體也會生成規模非常大的離散矩陣，并且容易受到數值頻散的影響.而積分方程類的數值算法具有半解析性質、只在散射異常體上進行網格離散即可.因此，對小型散射異常體的散射波場數值模擬，積分方程法生成的離散矩陣規模小，降低了計算和存儲的要求，其計算速度和精度要明顯優于基于微分方程的數值算法.對地震散射問題中的大尺度復雜地質條件而言，由于積分方程法形成的線性方程組系數矩陣的密實性，并且離散矩陣元素的計算過程復雜，計算耗時、存儲量大，相比于計算電磁學，基于積分方程的數值方法在地球物理中的應用受到了非常大的限制(王德智，2015).

針對積分方程法處理散射問題時所存在的問題，國內外學者在積分方程快速算法方面作了大量研究工作.由于求解積分方程數值法生成的大型線性方程組非常耗時，有學者提出采用近似解法來避免解方程組.Zhdanov和Fang(1996)、Zhdanov等(2000)提出了一種擬線性的近似解法，認為散射異常場是背景場與反射系數的乘積，反射系數是最小二乘問題的解，該方法在電磁散射正反演中取得較好的應用效果.Sun(2005)根據計算電磁學的擬解析線性近似技術，將電反射率張量的概念引入直流電場積分方程數值法中，給出了擬線性級數的表達形式，提供了一種新的直流電場的正反演方法.孫建國(2006)借鑒電磁散射問題中的擬線性近似技術，通過假設散射波場與背景波場之間的線性關系，建立了聲波散射的擬解析近似表達式.擬解析近似是半解析解，只需計算數值積分，該方法可以有效的避開了復雜介質情況下大型線性方程組的求解.

殷長春等(2018)利用插值技術將準線性近似算法應用于三維情況下航空電磁數據的數值模擬計算，在稀疏網格上計算散射場和點反射系數，然后進行線性插值，最后根據擬線性近似技術來求解電磁異常體的散射波場.另一類方法利用格林函數的Toeplitz性質，通過快速傅氏變換(Fast Fourier Transform，簡稱FFT)來加速線性方程組的迭代類求解方法中的矩陣向量乘積，稱為積分方程快速傅里葉變換法(Integral Equation Fast Fourier Transform，簡稱IE-FFT)(Peters and Volakis,1988; Volakis and Barkeshli,1991)，達到降低存儲和計算復雜的目的，如共軛梯度快速傅里葉變換(Conjugate Gradient Fast Fourier Transform，簡稱CG-FFT).Osnabrugge 等(2016)在背景波數中引入虛部分量，得到了帶阻尼因子的背景格林函數.通過加入預條件因子進一步加速了迭代的收斂性，同時利用FFT實現了循環卷積的快速算法.該方法在光學散射數值模擬中表現出了精度高、計算效率快的優點.

此外，傳統的求解積分方程的矩量法計算復雜，離散后的線性方程組系數矩陣每個元素都需要計算內積積分，計算量大，效率低.雖然快速多極子方法(Fast Multiple Method，簡稱FMM)和多層快速多極子方法(Multievel Fast Multiple Algorithm，簡稱MLFMA)(Rokhlin,1990; ?akir,2006; ?akr,2009；崔曉娜和劉洪，2020)能夠有效地降低存儲和計算復雜度，但對積分核依賴性強，且實現過程較為復雜.

除了上述方法之外，在計算數學中還有一種解積分方程的方法，稱為Nystr?m法.該方法利用適當的數值求積公式來近似積分算子(Nyst?m離散)，可以避免生成系數矩陣元素時的大量積分計算，相較于傳統的矩量法計算成本較低.如果積分方程的核函數只由Green函數組成，則基于數值積分公式的Nystr?m法所生成系數矩陣只對積分核進行簡單計算并且具有Toeplitz性質，計算復雜度低且易實現并行化.對于奇異性的處理，可以采用階數較高的插值基函數、奇異性減法等技術均可處理奇異性(Canino et al.，1998).

然而，在地震散射問題中出現的L-S方程的積分核除了Green函數之外還包含有速度擾動項.因此，為了使得積分方程的核函數只包含Green函數，首先將L-S方程改寫為一個與原積分方程等價的方程(等價L-S方程)，并對該方程進行逐點歸一化，進而得到在離散后其系數矩陣具有Toeplitz性質的歸一化等價L-S方程.然后，對歸一化等價L-S方程進行Nystr?m離散，從而得到其系數矩陣具有Toeplitz性質的線性方程組，極大地降低了存儲量.同時，為了提高空間褶積的計算效率，將歸一化等價L-S方程變換到波數域，在其傅里葉譜空間中實現矩陣和向量之間的乘積，進而實現對線性方程組的快速求解.通過簡單的分析可知，上述兩個措施可將內存和計算復雜度由O(N×N)分別降為O(N)和O(NlogN)(N=NxNz，Nx、Nz分別為x、z方向上的網格剖分數量).此外，因數值模擬在頻率域進行，各頻率下的波場值計算是獨立的.因此，為進一步提高計算效率，設計了基于MPI+OpenMP的進程級和線程級混合并行的實現模式，在有效降低其內存消耗的基礎上最大限度地提高積分方程法散射波場數值模擬的計算效率.

1 等價與歸一化等價L-S方程

利用擾動理論，可將描述散射問題的Helmholtz方程轉化成為Lippmann-Schwinger(L-S)積分方程：

U(r)=Ub(r)+Us(r)

(1)

由于格林函數G(r,r′)中的Hankel函數在r=r′點包含In(k0|r-r′|)項而奇異，考慮Cauchy主值意義下的收斂性，(1)式右端的積分項在點r∈Ω上一致連續(Fu et al.,1997)，即L-S方程具有弱奇異性.應用Fredholm擇一定理，可知L-S方程存在唯一解.

(2)

(3)

這是一個積分核為Green函數的積分方程.因此，如果首先對方程(2)進行某種方式的數值離散，然后再進行形如方程(3)的逐點歸一化處理，則可得到其系數矩陣具有Toeplitz性質的線性方程組.

2 等價L-S方程的Nystr?m離散

令Pn是插值算子，{rj,j=1,2,…,n}?Ω是插值基點，{lj(r′),j=j=1,2,…,n}是插值基函數，滿足lj(r′i)=δij.于是：

(4)

定義弱奇異積分算子：

(5)

則其近似算子滿足：

(6)

其中權函數：

(7)

(8)

(9)

定義常元插值算子：

(10)

其中，Mj為剖分單元Vj的中心.此時，需要計算弱奇異積分：

(11)

從而得到權系數.由于當r∈Vj且r=r′時G(r,r′)奇異，所以必須要對(11)進行奇異性消除處理，具體見劉寧和孫建國(2007)或張金會和孫建國(2009) .

3 線性方程組的逐次超松弛迭代

Lippmann-Schwinger積分方程的常元Nystr?m近似生成的系數矩陣為稠密矩陣，并且每一個矩陣元都必須計算格林函數在剖分元上的弱奇異積分.當散射體規模較大時，網格剖分數量也很大，在二維情況下，系數矩陣由N2個元素組成(N=NxNz，Nx、Nz分別為x、z方向上的網格剖分數量).采用高斯消元法或LU分解等直接法求解時需要存儲整個矩陣及其逆矩陣，內存需求為O(N2)并且矩陣求逆計算復雜度為O(N3).而迭代法無需存儲逆矩陣，其復雜度下降為O(N2).常用的迭代算法有Gauss-Seisel迭代、超松弛迭代(SOR)、共軛梯度法、廣義最小殘差法以及其他Krylov子空間法(Kleinman et al., 1990a,b; Kleinman and Van Den Berg,1991; Yu and Fu, 2014).本文采用逐次超松弛迭代法(Kleinman and Van Den Berg, 1991)來求解L-S方程離散生成的線性系統.

(12)

(13)

式中，βn為復逐次超松弛因子.根據Kleinman和Van Den Berg(1991)的超松弛迭代在光學散射中的應用經驗，通過在迭代過程中最小化殘差‖rn‖ 來得到逐次復松弛因子βn，即取βn使得 ‖rn‖ 在Hilbert空間中取極小：

(14)

4 矩陣向量乘積的FFT實現

(15a)

(15b)

其中：

(15c)

對于原L-S積分方程，核函數為α(r′)G(r,r′).由于速度擾動函數在每個面元上的值不同，導致帶速度擾動權的積分核經數值離散后的系數矩陣不具有Toeplitz性質.因此，有必要對L-S方程進行改寫，使得積分核僅包含Green函數，保證離散后的系數矩陣具有Toeplitz性質.

任給一n階循環矩陣C=circ(c0,c1,…,cn)可以被傅里葉矩陣F對角化(Davis，2013)，即:

(16)

(17)

對于一個n階對稱Toeolitz矩陣T，可以通過將它嵌到2n階循環矩陣C中，來計算矩陣向量乘積，即:

(18)

從而有:

(19)

根據式(19)可以得到矩陣向量乘積Tx.

(20)

其中:

0≤i≤Nx-1，

(21)

5 頻率域散射波場正演的并行化

MPI并行標準庫通過消息傳遞來實現程序的并行化，是一種進程級的并行工具.它具有效率高、可移植性強等優勢，是目前實際并行計算中普遍采用的并行編程環境.OpenMP作為一種在內存共享基礎上實現多個線程間的并行計算編程接口，通過簡單的高級語言指令即可實現多線程并行化計算(Quinn，2004).Nystr?m法求解Lippmann-Schwinger積分方程是在頻率域中實現的，各頻率片之間是獨立計算的，可以通過MPI并行處理技術實現頻率波場計算并行化.此外，當地下散射異常體規模較大時，剖分網格數量將會很大，求解單頻率下的線性方程組時通常含有大量稠密矩陣向量乘法運算相關的for循環，在單個節點上使用OpenMP實現for循環的并行化可進一步提高效率.基于MPI+OpenMP的混合并行模式可以減少MPI的通信開銷，充分利用資源，得到更高的加速比，在有效降低其內存消耗的基礎上最大限度地提高積分方程法散射波場數值模擬的計算效率.

5.1 基于MPI的頻率并行

基于MPI并行實現頻率片并行化的算法設計過程中，為獲得單炮所有頻率下的波場值，可以采用主從式并行模式，主進程進行消息發送與接收以及其他輸入輸出任務，從進程負責所有頻率的波場值計算.具體并行實現為主進程進行各從進程頻率數和偏移量的計算并將消息發送給從進程，當從進程計算完相應頻率數的波場值后，主進程接收各從進程發送的相應頻率數、偏移量和波場值，最終得到所有接收道和所有頻率波場值生成的二維數組.基于MPI并行的散射波場數值模擬實現步驟如下：

(1)調用MPI_Init函數,并生成通信域，參數初始化；

(2)根據集群中的節點數量，調用MPI_Comm_size函數來確定MPI_COMM_WORLD中的總進程數量(一個節點執行一個進程)，并調用MPI_Comm_rank確定進程編號；

(3)主進程根據設定的總進程數計算各從進程需要計算的頻率數和頻率偏移量并發送給相應的從進程；

(4)各從進程接收所需計算的頻率數和頻率偏移量，并計算相應頻率數下的所有道的波場值，并將計算結果發送回主進程；

(5)主進程接收從進程發送的頻率數、頻率偏移量和計算結果，對所有頻率下的波場值進行逆傅氏變換得到時域單炮記錄；

(6)調用MPI_Finalize，結束并行環境.

5.2 基于OpenMP的節點內for循環并行

L-S方程的Nystr?m近似迭代解法實現過程中涉及以for循環形式出現的稠密矩陣向量乘法運算，而線程級并行的OpenMP并不需要進行復雜的線程創建、同步、銷毀等工作.基于OpenMP的節點內循環并行化實現步驟如下：

(1)按照單個計算節點上的CPU核的數目，調用omp_set_num_threads設定線程數；

(2)在無數據相關性的for循環外部調用#pragma omp parallel for語句，將該循環并行化；

(3)調用private(〈variable list〉)和shared(〈variable list〉)設置并行域內的私有變量和共享變量.

6 數值試驗

下面通過數值試驗來測試基于Nystr?m法的IE-FFT并行算法在解決地震散射波場正演模擬中的有效性以及計算效率.首先設計了球形單體散射模型，通過與直接解法得到的數值結果進行對比，分析IE-FFT算法在存儲和計算效率上的優勢.然后用更為復雜的SEG/EAGE的鹽丘模型對快速算法進行進一步測試.數值試驗環境：Linux操作系統，包含11個節點，每個節點2個CPU，每個CPU包含8個核(CPU型號：Intel(R) Xeon(R) CPU E5-2620 0 @ 2.00 GHz )，使用基于MPI+OpenMP的并行算法對散射模型進行單炮數值模擬.

6.1 球體散射模型

球狀散射體速度模型為均勻的背景介質中含有一半徑200 m的高速球狀散射體(如圖1).模型大小為3000 m×3000 m，球散射體速度為3000 m×3000 m，背景速度2000 m·s-1，震源為主頻15 Hz的雷克子波，震源位置rs=(1500 m,10 m)，檢波器置于z=10 m的界面上，空間采樣間隔為dx=dz=5 m，時間采樣間隔為4 ms，記錄時間2.4 s.

圖1 球體散射模型Fig.1 Sphere scattering model

基于MPI+OpenMP，開辟11個進程進行基于MPI的頻率并行計算，每個進程開6個線程進行基于OpenMP的for循環并行計算，分別采用直接法和IE-FFT迭代算法計算球體散射模型正演結果.對于單球體散射速度模型，僅需要在散射異常體上進行離散.兩種方法得到的單頻波場(f=15 Hz)和單炮記錄結果如圖2和圖3，IE-FFT迭代快速算法與直接法求解得出的數值結果具有較好的一致性.為了更細致的分析，分別從兩種算法得到的單炮記錄中隨機抽取2道進行對比(如圖4)，結果顯示出，IE-FFT算法得到的數值結果在同相軸能量上與直接求解法結果有極小差別，量級在誤差允許范圍之內，兩種方法同相軸相位對應一致.

圖2 球體散射模型直接法與IE-FFT迭代法15Hz單頻波場數值結果(a) 直接法單頻波場實部； (b) 直接法單頻波場虛部； (c) IE-FFT迭代法單頻波場實部； (d) IE-FFT迭代法單頻波場虛部.Fig.2 Numerical results of 15Hz single-frequency wave field using the direct method and IE-FFT iterative method on the sphere scattering model(a) The real part of the single-frequency wave field using direct method； (b) The imaginary part of the single-frequency wave field using direct method； (c) The real part of the single-frequency wave field using the IE-FFT iterative method； (d) The imaginary part of the single-frequency wave field using the IE-FFT iterative method.

圖3 球體散射模型直接法(a) 與IE-FFT迭代法; (b) 數值結果.Fig.3 Numerical results using the direct method (a) and IE-FFT iterative method (b) on the sphere scattering model

圖4 球體散射模型直接法(實線)與IE-FFT迭代法(點劃線)單道對比(a) 第40道單道對比； (b) 第150道彈道對比.Fig.4 Comparison between the direct method (solid line) and the IE-FFT iterative method (dotted line) for the sphere scattering model(a) Comparison of trace 40； (b) Comparison of trace 150.

球體散射模型直接法與IE-FFT迭代法的CPU時間與內存消耗對比見表1.數值實現過程中，直接法耗費CPU時間為4624.27 s，每個進程內存消耗為515 M，IE-FFT法迭代50次耗費CPU計算時間為396.37 s，每個進程內存消耗為78 M.IE-FFT迭代算法的數值結果可以較好地吻合直接法結果，而計算速度是直接法的11倍，并且應用格林函數的Toeplitz性質降低了系數矩陣的內存需求，內存消耗僅約為直接法的1/7.

6.2 SEG/EAGE鹽丘模型

當地下介質中包含較大尺度的復雜散射體且速度擾動更強時，積分方程法應用于此類介質模型時需要在全空間中進行離散計算.而直接法因其離散后的矩陣為NxNz階方陣，普通的微機難以獲取相應的棧內存來存儲系數矩陣且計算效率不滿足實際應用需求.為了驗證IE-FFT迭代算法在復雜強散射介質中的適應性，選取SEG/GAGE鹽丘速度模型.模型參數如圖5所示，模型大小為7800 m×2090 m，背景速度1500 m·s-1，震源子波為主頻15 Hz的雷克子波，源位置rs=(3800 m,10 m)，檢波器置于z=10 m的界面上，迭代法空間采樣間隔為dx=dz=10 m，時間采樣間隔為4 ms，記錄時間4 s.

圖5 SEG/EAGE鹽丘模型Fig.5 SEG/EAGE salt dome model

IE-FFT快速迭代算法得到的單頻波場(f=15 Hz)和單炮地震記錄如圖6和圖7所示.從單炮地震記錄上可以清晰分辨鹽丘頂端以及巖體右側橫向起伏變化劇烈處的散射同相軸，對強散射體的識別能力較強，并且模型底部測試層的同相軸清晰、能量強，并未受到高速異常體的屏蔽影響.在直接法無法求解大型線性方程組的情況下，IE-FFT快速迭代法迭代100次耗費CPU計算時間和每個進程內存消耗為為9466.54 s和321 M(見表1).

圖6 IE-FFT迭代算法15HZ單頻波場數值結果(a) 單頻波場實部； (b) 單頻波場虛部.Fig.6 Numerical results of 15Hz single-frequency wave field using the IE-FFT iterative method on the SEG/EAGE salt dome model(a) The real part of the single-frequency wave field； (b) The imaginary part of the single-frequency wave field.

圖7 IE-FFT迭代算法散射波場數值結果Fig.7 Numerical result using the IE-FFT iterative method on the SEG/EAGE salt dome model

表1 直接法與IE-FFT迭代法CPU時間與內存消耗Table 1 CPU time and memory consumption of the direct method and the IE-FFT iterative method

7 結論

傳統上，對L-S方程的數值求解一般采用矩量法.此時，相應線性方程組的系數矩陣元素都是以內積積分的形式出現的.雖然可以利用數值積分法計算這些內積積分且可以獲得較高的精度，但是其計算過程復雜、對計算機資源要求高且效率低.與此相反，基于數值積分公式的Nystr?m離散可使系數矩陣只包含對積分核(Green函數)的積分計算.因此，從計算效率方面考慮，對于大尺度地質條件下的地震散射波場計算問題，Nystr?m法會更加適合.事實上，通過對L-S方程的變形處理，可使得最終形成的線性方程組的系數矩陣具有Toeplitz性質，進而可將存儲量由對N階矩陣直接進行存儲的O(N×N)降為O(N).此外，利用FFT和IFFT計算離散空間褶積，可以加速線性方程組求解迭代計算過程中的矩矢乘積，使其計算復雜度由O(N×N)降為O(NlogN).再有，根據各單頻波場在計算上相互獨立的特點，采用了基于MPI+OpenMP的進程級和線程級相結合的并行化計算方案.數值結果表明，相較于傳統的積分方程數值算法，利用Nystr?m離散改寫后的等價L-S方程且同時采用FFT計算離散褶積和并行計算的數值求解方案可以得到與直接求解法相當的數值結果，在保證計算精度的同時，內存和CPU耗時則顯著降低，有效地克服了傳統數值算法在求解L-S方程時系數矩陣存儲和線性方程組求解效率低的問題.