基于二維彈性波方程的高階NAD-SPRK算法

陳 麗，朱興文，張朝元*

（1.大理大學工程學院，云南大理 671003；2.大理大學數學與計算機學院，云南大理 671003）

開發精度高、效率高的模擬地震波計算方法是目前涉及地震波動方程開展反演的一個重要方向。過去的計算格式〔1-5〕在進行地震波模擬時會造成嚴重的數值頻散現象而達不到當前規模大的地震波模擬效度，楊頂輝教授團隊于2003年首次在地震波正演模擬中引入近似解析離散化（nearly analytic discrete，NAD）算子〔6〕，目前已獲得系列NAD算子的數值算法〔7-10〕，這些算法的模擬頻散效果均很不錯，然而這些算法僅為四階的空間精度。

為提高地震波計算模擬效度，本文針對具有哈密爾頓系統的二維彈性波方程，結合離散空間高階偏導數的八階NAD算子和離散時間導數的二階辛分部Runge-Kutta算法〔4，11〕，得到了八階NADSPRK算法。針對該方法，從理論和數值計算兩方面研究了其穩定性條件、數值頻散和計算效率。結果表明：同四階NSPRK算法〔11〕、八階Lax-Wendroff（LWC）算法〔12〕和八階交錯網格（SG）算法〔12〕相比，八階NAD-SPRK算法壓制數值頻散的能力顯著優于傳統數值計算方法，且具有最小的數值誤差和最高的計算效率。

1 方法推導

設二維彈性波方程為：

其中，j=1，3，fi和ui分別為i方向的力源分量和位移分量，ρ=ρ（x，z）為介質密度，σij為應力張量。

由應力與張力之間的物理關系，方程（1）式的向量方程為：

記vi=?ui/?t（i=1，2，3），則v=（v1，v2，v3）T，方程（2）式為：

其中D稱為偏微分算子矩陣，且

O為三階零方陣，I為三階單位矩陣。

此時，表達式（4）擁有哈密爾頓系統模式，故采納哈密爾頓系統模式的求解流程對二維彈性波方程表達式（1）的對應方程（4）式進行求解。具體求解如下：

1）由于涉及位移u及粒子速度v的有關二階、三階等偏導數存在于高階偏微分算子矩陣D中，因此首先需要對方程（4）式右邊的D進行空間高階偏導數離散化。逼近位移u和粒子速度v的二階和三階偏導數的計算公式詳見文獻〔12-13〕。

2）對D進行離散化處理后，表達式（4）成為關于時間的一個半離散化哈密爾頓系統模式。再采納辛分部的二階Runge-Kutta算法〔4，11〕對半離散化哈密爾頓系統表達式（4）進行求解，計算公式為：

其中，△t表示時間步長，V1表示計算的中間變量。

為加快計算速度及減少存儲空間，簡化（5）式得到計算公式：

其中D2=D·D。

以上二維彈性波方程表達式（1）的對應求解表達式（5）和式（6）記為八階NAD-SPRK算法。

為了揭示八階NAD-SPRK算法在進行地震波場數值模擬時的數值頻散和計算效率情況，下面從理論和數值計算兩方面分析該方法的穩定性條件、數值頻散關系和計算效率，并同四階NSPRK算法、八階LWC算法和八階SG算法進行比較。

2 穩定性分析

在利用有限差分方法進行波場數值模擬時，為保證計算的穩定性，應該選取恰當的模擬參數，就是在選定空間增量h及波速c后，時間增量△t的選擇也要達到相應的穩定條件。本文使用Fourier方法〔13〕推導出二維聲波情形下八階NAD-SPRK算法的穩定性，其穩定結果為：

其中，h=△x=△z為空間步長，α=c△t/h為庫朗數，αmax為最大庫朗數。

對于彈性波情況，推導出彈性波方程數值方法的精確或解析穩定性條件通常是很復雜的〔5〕。當應

用八階NAD-SPRK算法求解二維彈性波方程時，其穩定性條件由系數估計技術（Vichnevetsky，1979）不能直接確定，而是由聲波情形下最大波速時的穩定性條件來近似代替。因此，估計時間網格尺寸應滿足以下穩定性條件：

式中tmax為使八階NAD-SPRK算法在二維彈性波情況下保持穩定的最大時間步長，cmax為最大縱波速度。

3 數值頻散

下面通過模擬波場來研究八階NAD-SPRK算法的數值頻散情況。為此，本文挑選二維彈性波介質具有均勻各向同性特點的方程如下：

此處，u1，u3分別表示位移在方向x和方向z的分量，ρ為介質密度，λ，μ為Lame參數，f1和f3為震源函數在方向x和方向z的分量。

在數值實驗中，選擇模型的參數為λ＝6.0 GPa、μ＝24 GPa、ρ＝1.5 g/cm3。網格點數為301×301，空間網格步長和時間步長分別為△x＝△z=40 m和△t＝1×10-3s。震源位于計算區域中心，且Ricker子波為f0的震源函數，且f1=f（t）、f3＝0。接收器位于R（5.25 km，4.8 km）。其中，f（t）表達式為：

圖1分別為由八階NAD-SPRK算法、四階NSPRK算法、八階LWC算法和八階SG算法通過模擬計算給出在Ricker子波頻率為f0＝25 Hz下從時間t＝0 s到t＝0.68 s在接收點R處的波形記錄。可以看到圖1（a）、圖1（b）和圖1（d）的波形幾乎相同，而由八階LWC算法模擬產生的波形記錄存在一定的數值頻散，見圖1（c）。

圖1 4種算法模擬產生在均勻各向同性介質和頻率為f0=25 Hz的Ricker子波下的波形記錄

圖2分別給出了由4種方法模擬產生在Ricker子波頻率為f0=40 Hz下從時間t=0 s到t=0.68 s在接收器R處的波形記錄。從圖2（c）和圖2（d）可了解到，經八階LWC算法和八階SG算法模擬計算出的波形記錄出現嚴重的數值頻散現象。而經八階NAD-SPRK算法和四階NSPRK算法計算得到的圖2（a）和圖2（b）的波形記錄清晰，且在較高頻率下數值頻散不明顯。尤其是由八階NAD-SPRK算法計算得到的圖2（a）在如此高的頻率情況下幾乎沒有可見的數值頻散。結果表明，八階NAD-SPRK算法在抑制數值頻散方面比八階LWC算法和八階SG算法更有效。

圖2 4種算法模擬產生在均勻各向同性介質和頻率為f0=40 Hz的Ricker子波下的波形記錄

4 計算效率

此部分，將呈現八階NAD-SPRK算法數值計算的模擬結果，也展示該方法數值計算的模擬效率。此處，同樣選擇介質具有均勻各向同性特點的二維彈性波表達式（9）。

此實驗中，模型的參數為λ=4.704 GPa、μ=8.4 GPa和ρ=2.1 g/cm3。網格點數為301×301，空間網格步長和時間步長分別為△x=△z=40 m和△t=0.004 s。震源位于計算區域中心，震源函數為f0=12 Hz的Ricker子波，且f1=f（t）和f3=0。

圖3分別為粗網格（△x=△z=40 m）下由八階NAD-SPRK算法和八階LWC算法在t=1.6 s時刻產生的水平位移方向和垂直位移方向的波場瞬時快照。從圖3可知用這兩種算法計算模擬彈性波給出的波前快照大致上是相同的。雖然在相同網格點數下，八階NAD-SPRK算法的計算成本比八階LWC算法高，這是由于八階NAD-SPRK算法模擬計算時需要同時計算位移、粒子速度及其梯度等變量。但是，八階NAD-SPRK算法模擬生成的波場快照（圖3（a）（c））在空間步長為40 m時數值頻散較小，而八階LWC算法數值頻散較嚴重（圖3（b）（d））。結果表明，八階NAD-SPRK算法能有效應用于粗網格大尺度情形下的數值模擬。

為了有效消除數值頻散，在與圖3有相同的庫朗數條件下，由八階LWC算法在細網格△x=△z=20 m和對應601×601的網格點數條件下在t=1.6 s時刻計算模擬出的波場快照見圖4，此時沒有可見的數值頻散。在達到相同消除數值頻散條件下，對于同一計算域，在△x=△z=40 m粗網格下，八階NAD-SPRK算法的網格點數僅為301×301。可見，八階NAD-SPRK算法模擬計算時所需的內存需求量約為八階LWC算法的25.08%。

對比圖3（a）（c）和圖4（a）（b），本文所提出的方法在相同庫朗數下，可以提供與八階LWC算法在細網格上相同的精度，但它們的計算成本是不同的。八階NAD-SPRK算法模擬得到圖3（a）（c）所消耗的CPU時間約185 s，而八階LWC算法模擬得到圖4（a）（b）所消耗的CPU時間約293 s。這揭示在沒有數值頻散可見的狀況下，八階NAD-SPRK算法的計算效率約為八階LWC算法在細網格下的1.6倍。以上實驗均在Intel（R）Core（TM）2Duo CPU和內存為2.33 G的計算機環境下進行。

圖3 粗網格△x=△z＝40 m下由八階NAD-SPRK算法和八階LWC算法生成的在t＝1.6 s時刻的波場快照

圖4 八階LWC算法在細網格△x=△z＝20 m下t＝1.6 s時刻的波場快照

為了進一步驗證八階NAD-SPRK算法的有效性，比較了八階NAD-SPRK算法和八階LWC算法在模擬二維彈性波方程的數值解。Dablain〔1〕的研究結果表明，在細網格條件下八階LWC算法提供的數值結果近似等于解析解，并與偽譜方法模擬的結果相等。

八階NAD-SPRK算法在空間網格步長和時間步長分別為△x=△z＝40 m和△t=0.004 s下模擬得到的結果與八階LWC算法在空間和時間步長分別為△x=△z=20 m和△t=0.002 s下模擬產生的結果基本相等。圖5為由八階NAD-SPRK算法在粗網格（△x=△z=40 m）上和八階LWC算法在細網格（△x=△z=20 m）上模擬得到的彈性波波形記錄，其中圖5（a）和圖5（b）分別為在接收器R（6.6 km，7.6 km）處的水平位移方向和垂直位移方向的波形記錄。結果表明，八階NAD-SPRK算法針對均勻各向異性情形在粗網格條件下可以得到相等的模擬結果，并且計算量和計算機內存都大大地減少。

圖5 細網格（△x=△z＝20 m，△t＝0.002 s）下的八階LWC算法和粗網格（△x=△z＝40 m，△t＝0.004 s）下的八階NAD-SPRK算法在接收器R（6.6 km，7.6 km）處位移水平方向對比的波形記錄

5 結論

為使地震波數值計算的正演方法的精度更高和速度更快，本文在楊頂輝教師團隊的研究基礎上，針對擁有哈密爾頓特征的二維彈性波方程，結合離散空間高階偏導數的八階NAD算子和離散時間導數的二階辛分部Runge-Kutta算法，推導出了一種八階NAD-SPRK算法。結合理論和數值兩角度進行了該方法的穩定性分析、數值頻散探索和計算效率比較等工作。數值結果表明：同四階NSPRK算法、八階LWC算法和八階SG算法相比，八階NAD-SPRK算法在壓制數值頻散上明顯好于傳統的數值模擬方法，計算量和計算機內存也都大大地減少。在達到相同消除數值頻散條件下，八階NADSPRK算法模擬計算時所需的內存需求量約為八階LWC算法的25.08%，計算效率約為八階LWC算法的1.6倍。可見，八階NAD-SPRK算法在全波形反演和波動方程的逆時偏移等問題上有望得到廣泛推廣和有效應用。