擴展的近似解析離散化方法及彈性波方程數值模擬

2017-11-01 23:56:46段焱文安一凡王笑叢桂志先

石油地球物理勘探 2017年5期

關鍵詞：方法

汪勇段焱文安一凡王笑叢桂志先

(①長江大學地球物理與石油資源學院，湖北武漢 430100； ②油氣資源與勘探技術教育部重點實驗室(長江大學)，湖北武漢 430100； ③東方地球物理公司國際勘探事業部，河北涿州 072751)

擴展的近似解析離散化方法及彈性波方程數值模擬

汪勇*①②段焱文①②安一凡①②王笑叢③桂志先①②

汪勇，段焱文，安一凡，王笑叢，桂志先.擴展的近似解析離散化方法及彈性波方程數值模擬.石油地球物理勘探，2017,52(5):928-940,955.

地震波場數值模擬是研究地震波理論、偏移成像和地震反演等工作的基礎，提高數值模擬的精度具有重要意義。在前人的研究基礎上，提出了一種擴展的近似解析離散化數值模擬方法，并從理論上對該方法的精度、數值頻散、穩定性和計算效率進行了分析。該方法在時間差分上比擴展前提高了一階精度、誤差最大降低了88%。與其他近似解析離散化類方法一樣，具有算子半徑小和適應粗網格步長的優勢，在最小主波長內僅需使用5.9個網格點。與四階Lax-Wendroff修正格式和交錯網格有限差分格式的頻散曲線和模擬結果對比，驗證了該方法能更好地壓制數值頻散。二維各向同性均勻介質和水平層狀介質模型數值模擬的地震彈性波場特征清晰準確，說明了該方法的實用性。

擴展最優近似解析離散化彈性波方程數值模擬穩定性條件數值頻散

1 引言

地震波場數值模擬是由已知的巖石結構和物理參數建立地震地質模型，依據不同的算法模擬地震波在地下介質中的傳播過程，計算在地面或地下各觀測點地震記錄的一種方法。目前，常用的數值模擬方法主要有射線追蹤法[1]和波動方程法，其中波動方程法有偽譜法[2，3]、有限元法[4]、邊界元法[5]、譜元法[6]和有限差分法。20世紀 60 年代，Alterman等[7]首先將有限差分方法用于地震波場的數值模擬，該方法現已成為地震波場數值模擬中應用最為廣泛的一種方法，差分格式也由早期的低階差分發展到高階差分，由常規網格發展到交錯網格、旋轉交錯網格、可變網格等[8-16]。近似解析離散化方法(Nearly-Analytic Discrete Method，簡稱NADM)是20世紀90年代發展起來的一種有限差分數值模擬方法[17，18]，其基本思想是在時間上采用泰勒公式展開，在空間上利用截斷的泰勒公式構造高階插值函數來逼近空間偏導數，與其他有限差分法的根本區別是它利用空間節點的位移和梯度共同逼近變量的空間偏導數，從而提高計算精度，在大網格下可有效地壓制數值頻散，提高計算效率。楊頂輝及其研究團隊對該方法進行了深入的研究[19-25]。

本文在前人研究的基礎上，對優化的近似解析離散化方法進行了擴展，即將位移場的泰勒展開由5項擴展為6項，將時間差分精度提高到4階，進而提高了數值模擬精度。

2 擴展的近似解析離散方法原理

近似解析離散化類數值模擬方法的基本思想是：在時間上采用泰勒公式展開，在空間上利用截斷的泰勒公式構造高階插值函數來逼近空間偏導數，從而得到原函數的一個最優近似。該方法的特點是在求解偏微分方程時包含質點位移和速度的各階空間偏導數，應用了原函數、各階偏導數之間的相互聯系，能有效減少離散過程中地震信息的丟失、提高數值計算的精度。以二維聲波方程為例討論該類方法數值模擬的基本思想。

二維聲波波動方程為

(1)

設滿足聲波方程的解為U，定義

(2)

(3)

式中m、k分別表示S對x和t的偏導數的階次。

利用截斷的泰勒公式表示n+1時刻的U，可以得到

(4)

式中：i和j分別為x和z方向空間網格點序號；n為時間序號。P及其時間的偏導數可以由波動方程式(1)轉化為U和S對空間的高階偏導，即

(5)

(6)

(7)

根據式(3)，如果式(6)中S對空間的二階偏導數可用U表示的一階向后時間差分近似，即

(8)

則稱為原始的近似解析離散化方法。

(9)

(10)

文獻[24]指出，由于采用式(8)表示的一階向后差分格式求取S，降低了整體數值模擬精度，導致了NADM方法精度僅為O(Δt2+Δx4+Δz4)。

為了提高時間層推進精度，對S同樣用泰勒公式展開，則可以得到速度場S的時間層推進差分格式

(11)

式(4)～式(7)以及式(11)構成了改進的近似解析離散化方法(Improving Nearly-Analytic Discrete Method，簡稱INADM)的差分格式。該方法與原方法相比，在增加了存儲量和計算量的情況下，提高了數模擬值時間精度。

此外，如果利用截斷的泰勒公式表示n-1時刻的U，可以得到

(12)

將式(4)和式(12)相加，可以得到最優近似解析離散化方法(Optimum Nearly-Analytic Discrete Method，簡稱ONADM)的差分格式

(13)

ONADM方法的核心是消除了差分格式式(4)中的速度項，避免了原方法中降低時間精度的因素，從而提高了整體差分精度。文獻[21]指出ONADM的精度為O(Δt2+Δx4+Δz4)。與INADM方法相比，該方法由于不涉及速度項，所以在計算效率和存儲空間方面得到了提高，但其缺點是不能進行帶速度項的波動方程數值模擬，如黏彈介質波動方程模擬等。

INADM和ONADM方法有各自的優缺點，ONADM的精度比INADM高，但INADM具有更廣泛的適用性，能夠對黏彈介質等復雜介質波動方程進行數值模擬。所以本文選擇對INADM方法進行修正或擴展，提高INADM中位移和速度項的差分精度是本文的研究目的。

(14)

因為按照INADM方法計算時，原本就需要計算速度場S，所以式(14)中的速度高階空間導數(≤5次)是容易直接求取的，這就為提高泰勒展開精度提供了基礎，進而U和S的時間層推進公式可以擴展為

(15)

(16)

式(15)和式(16)即為擴展的近似解析離散化方法(Expand Nearly-Analytic Discrete Method，簡稱ENADM)的差分公式。式中P的各階時間導數分別用式(5)～式(7)和式(14)計算，其中涉及到的空間高階偏導數則與其他NAD類方法一樣，通過質點的位移場和速度場及其梯度逼近計算。在上述差分公式推導中，雖然介質的地震波速度v是常數，但在數值模擬中可以隨空間位置而改變，即上述方法依然可以用于層狀介質和非均勻介質等復雜介質的數值模擬研究。

3 ENADM方法分析

3.1 精度分析

NAD類有限差分法中的空間導數是利用待求空間網格點周圍點的位移及梯度共同逼近， Yang等[19]指出，它們的空間精度為四階，所以ONADM、INADM和ENADM的空間精度相同，在后面的頻散關系討論中，也可以看出其頻散曲線基本一致。ONADM采用式(13)進行位移的時間層推進，由于消掉了速度項，所以它的時間精度為四階。INADM的位移和速度分別采用式(4)和式(11)進行時間層推進，位移的時間精度為四階，而速度為三階，所以整體時間精度應為三階。ENADM在INADM基礎上進行了擴展，其位移和速度按式(15)和式(16)分別進行計算，所以位移和速度的時間精度分別為五階和四階，整體的時間精度為四階，比INADM提高一階時間精度，所以ENADM的差分精度為O(Δt4+Δx4+Δz4)。

應用二維平面諧波初值問題，比較四階Lax-Wendroff修正格式(簡稱LWC)、最優近似解析離散化方法(ONADM)、改進的近似解析離散化方法(INADM)和擴展的近似解析離散化方法(ENADM)數值模擬精度。二維平面波初值問題可以表示為

(17)

式中：θ0是初始時刻波陣面法線方向(即傳播方向)與x軸的夾角；f0是平面簡諧波的頻率。其準確解析解為

(18)

LWC差分格式[26]為

(19)

二維波場數值模擬中，假設平面波頻率f0=20Hz,θ0=π/4。設置均勻介質模型，波速為4000m/s，模型長度和深度均為2000m，縱、橫向網格間距相同。在不同空間步長和時間步長條件下，計算前述四種數值模擬結果的相對誤差。相對誤差定義為

(20)

表1 四種方法在不同情況下的最大相對誤差(%)比較

圖1 四種差分方法的相對誤差曲線(a)Δx=10m，Δt=0.5ms； (b)Δx=15m，Δt=0.5ms； (c)Δx=20m，Δt=1.0ms； (d)Δx=25m，Δt=2.0ms

從表1和圖1可以看出，四種數值模擬方法相對誤差均隨時間步長或空間步長的增大而增加。從誤差統計來看，ENADM具有最高的數值精度，INADM和ONADM次之，而LWC的數值精度最低、誤差增長最快。當空間步長為25m，時間步長為2ms時，ENADM最大相對誤差僅為0.45%，而擴展前的INADM卻有3.74%，其相對誤差最大降低了88%，能更好地用于大尺度模型的長時間數值模擬。

3.2 頻散分析

數值模擬中由于求解的誤差而導致數值波速與頻率有關，即產生所謂的數值頻散現象，這不是地震波特征真實的反映，所以頻散關系分析是判斷一種數值模擬方法優劣的重要依據，其結果決定了該方法能否或在何種條件下適用于地震波正反演計算。選取LWC、ONADM、INADM、ENADM和交錯網格有限差分(Staggered-grid Finite Difference Method，簡稱SFDM)五種方法進行二維頻散關系的比較和分析，以說明ENADM方法在壓制數值頻散方面的優越性。

ENADM方法頻散關系(推導過程見附錄A)為

φxsinφx+4(cosφz-1)+φzsinφz]-

3[φzsinφz+2(cosφz-1)]-2(cosφxcosφz-

cosφx-cosφz+1)}

(21)

確定φ后，代入上述頻散關系，可以解得對應的ω。定義數值波速與真實速度的比值為

(22)

在理想情況下，如果不存在數值頻散則速度比γ恒等于1。γ偏離1越大，說明該方法的數值頻散越嚴重，反之則說明該方法能更好地壓制數值頻散。

同理可以得到其他四種差分方法的頻散關系。INADM頻散關系為

φxsinφx+4(cosφz-1)+φzsinφz]-

(23)

ONADM頻散關系為

exp(-iω)=2-exp(iω)+α2[4cosφx+Δxkxsinφx+

4cosφz+Δxkzsinφz-8]+

(24)

LWC頻散關系為

cosω-1

(25)

SFDM頻散關系[27]為

(26)

圖2為五種方法在不同的α和θ時的頻散關系曲線。

圖2 五種差分方法在不同參數時的頻散曲線

取φ∈[0,π]作為橫坐標，它是波數與空間步長的乘積。當空間步長確定后，φ的增加表示頻率隨之增加。另一方面，單位波長內采樣點數M=2π/φ，所以橫坐標也可以看作M由∞逐漸減小至2。從圖中的頻散曲線可以看出：①隨著空間采樣點數的減小，5種方法的頻散現象逐步加劇，而NAD類方法的數值頻散比LWC和SFDM方法要小，其頻散曲線更趨近于1，說明了這類方法在壓制數值頻散方面的更具優勢；②假設數值速度在理論速度的99.9%以內表示不存在頻散，則NAD類、SFDM和LWC對應的最小的φ分別為0.337、0.274和0.131，所以它們在最小主波長內需要使用的網格點數分別為5.9、7.3和15.3個；③由于三種NAD類方法在求取空間偏導數時所用差分格式和精度相同，所以它們的頻散曲線基本重合。在θ=π/4,α=0.45時，ENADM方法的速度比最接近1，說明該方法在大Courant數和對角線方向傳播時，它能更好地壓制頻散；④不同的θ，即波在不同的傳播方向上，速度比不同，說明存在一定的數值各向異性，α越大越明顯。當α=0.45時，ENADM、INADM、ONADM、LWC和SFDM在三個方向上的速度比最大相差分別為8.2%、8.4%、8.3%、20.7%和9.2%，說明NAD類方法對數值各向異性也有很好的壓制，且ENADM效果最好。

圖3 二維聲波均勻模型不同方法800 ms時刻模擬的波場快照

圖4比較了三種方法在(x=300m，z=1500m)處接收到的波形記錄。圖4a中的黑色實線為ENADM數值模擬結果，紅色虛線為精確解析解，二者基本重合，模擬結果非常理想。與同是四階精度的LWC和SFDM相比，ENADM的質點振動曲線沒有LWC和SFDM中的“拖尾”現象，在粗網格數值模擬時具有更高的計算精度。需要說明的是，數值模擬中三種方法均采用了相同的時間和空間網格，但SFDM采用一階速度—應力方程，模擬結果是應力分量地震記錄； ENADM和LWC采用二階位移方程，模擬結果是位移地震記錄，導致了三種方法所得地震記錄振幅在數量級上存在一定的差別。為了比較，圖4中的三個震動曲線經過了振幅均一化處理。

圖4 在(x=300m，z=1500m)處三種方法波形記錄

3.3 穩定性分析

穩定性條件是有限差分數值模擬中一個非常重要的問題，如果不滿足穩定性條件，則前面各層的舍入誤差將會影響到后面各層的值，后面的誤差的影響越來越大，會使有限差分的結果完全錯誤。本文采用Fourier方法[28]對ENADM的差分格式進行穩定性分析。記

(27)

(28)

式中vmax為波的最大傳播速度。

3.4 計算效率分析

對上面的二維均勻模型進行聲波數值模擬，采樣時間為1s。LWC和三種NAD類方法使用的空間網格間距Δx和時間網格間距Δt的選擇滿足無數值頻散和穩定性的要求。四種方法的計算效率、使用的數組大小、個數如表2所示。

表2 不同方法計算效率比較對比

從表2可看出，由于ENADM算法提高了時間精度，導致計算時所需要的內存和時間增加，與擴展前的INADM相比，其占用內存增加了28.6%，計算時間增加了28.4%，即ENADM方法在增加存儲和計算量的條件下提高了數值模擬的精度。同時，因為計算了速度分量，所以很容易進行含速度項的波動方程數值模擬研究，在實際數值模擬中，可以根據介質類型和精度要求，選擇合適的數值模擬方法。

4 模型試算

采用ENADM方法對二維各向同性完全彈性介質的彈性波方程進行數值模擬，波動方程為

(29)

式中：vP和vS分別為縱、橫波波速；ux和uz分別為彈性波在x和z方向的位移，對其進行泰勒公式展開可得

(30)

4.1 均勻介質模型

均勻模型介質為泊松體，尺寸為2000m×2000m，縱波速度為4000m/s，橫波速度為2310m/s。由于最小速度為2310m/s，按頻散分析的結果，每個波長內需采樣5.9個點，可取的最大空間步長為9.8m，為確保精度，設定縱、橫向空間間隔均為8m。由于最大速度為4000m/s，按穩定性條件，設置時間步長為0.8ms。在模型中心處激發40Hz的Ricker子波等能量震源。圖5為220ms時刻的波場快照，波場非常清晰，沒有出現數值頻散現象。圖5中外層為縱波波前，內層為橫波波前。由圖可以看出，采用等能量震源時：①圖5a顯示的位移水平分量上，當x為1000m時，由于縱波質點垂向振動，橫波質點橫向振動，所以縱波位移為0，橫波最大，且縱波波前相位以x=1000m左右對稱。當z=1000m時，則縱波最大，而橫波為0，且橫波波前相位以z=1000m上下對稱。②圖5b顯示的位移垂直分量上，縱波和橫波波前以x=1000m和z=1000m為對稱軸，相位相反。

圖5 均勻介質模型彈性波方程ENADM數值模擬波場快照

4.2 水平層狀介質模型

設置兩層水平層狀介質模型每層厚度均為750m，模型長度為1500m。第一層縱、橫波速度分別為3800m/s和2600m/s；第二層縱、橫波速度分別為4500m/s和2900m/s。在地面中間處(x=750m，z=0m)激發20Hz雷克子波震源，震源加載到x和z方向位移上。縱、橫向空間步長均為15m，時間步長為1ms，采用PML人工邊界條件[29，30]。圖6為地面接收到的地震記錄，長度為1s，記錄的縱橫波波場清晰，沒有數值頻散。圖中ZP和ZS分別表示直達縱波和直達橫波，PP表示反射縱波，SS表示反射橫波，PS和SP表示在反射界面形成的轉換波，它們達到檢波點的時間相同，相互干涉疊加。

總體來看，雙層模型中地面模擬接收到的地震記錄非常清晰，沒有數值頻散和不穩定數值結果，地震記錄中的直達波、反射波和轉換波顯示清楚，說明ENADM算法可以有效地模擬彈性波在多層各向同性介質模型中的傳播。

圖6 水平層狀介質模型彈性波方程ENADM數值模擬地面地震記錄

5 結論與建議

本文在前人的研究基礎上，提出并推導了二維各向同性介質彈性波方程的擴展近似解析離散化數值模擬方法。通過理論分析和模型試算，進行了該方法的模擬精度、頻散關系分析、穩定性研究及計算效率分析。研究結果表明，擴展的近似解析離散化數值模擬方法具有模擬精度高、適用粗網格計算的優勢，在大尺度模型波場模擬時具有明顯的優勢，在研究深部地球物理和石油地球物理勘探等方面有廣泛的應用前景。本文僅將該方法應用于二維各向同性介質彈性波模擬，但該方法能進一步推廣到三維、各向異性、黏彈及雙相介質等復雜介質的模擬，也可以采用同樣的思路將其他NAD類方法(如加權近似解析離散)進行擴展，為研究地震波在復雜介質中的傳播規律提供了一種新的方法和思路。

[1] Cerveny V.Seismic rays and ray intensities in inhomogeneous anisotropic media.Geophyscial Journal International，1972，29(1)：1-33.

[2] Kosloff D D，Baysal E.Forward modeling by a Fourier method.Geophysics，1982，47(10)：1402-1412.

[3] 唐懷谷，何兵壽.一階聲波方程時間四階精度差分格式的偽譜法求解.石油地球物理勘探，2017，52(1)：71-80. Tang Huaigu，He Bingshou.Pseudo spectrum method of first-order acoustic wave equation finite-difference schemes with fourth-order time difference accuracy.OGP，2017，52(1)：71-80.

[4] Drake L A.Rayleigh waves at a continental boundary by the finite element method.Bulletin of the Seismological Society of America，1972，62(5)：1259-1268.

[5] Bouchon M，Coutant O.Calculation of synthetic seismograms in a laterally varying medium by the bounda-ry element-discrete wavenumber method.Bulletin of the Seismological Society of America，1994，84(6)：1869-1881.

[6] Komatissch D，Barnes C，Tromp J.Simulation of anisotropic wave propagation based upon a spectral element method.Geophysics，2000，65(4)：1251-1260.

[7] Alterman Z，Karal F C.Propagation of seismic wave in layered media by finite difference methods.Bulletin of the Seismological Society of America，1968，58(1)：367-398.

[8] Kelly K R，Ward R W，Treitel S et al.Synthetic seismograms：a finite-difference approach.Geophysics，1976，41(1)：2-27.

[9] Virieux J.P-SV wave propagation in heterogeneous media：Velocity-stress finite-difference method.Geo-physics，1986，51(4)：889-901.

[10] Graves R W.Simulation seismic wave propagation in 3D elastic media using staggered-grid finite differences.Bulletin of the Seismological Society of America，1996，86(4)：1091-1106.

[11] 董良國，馬在田，曹景忠等.一階彈性波方程交錯網格高階差分解法.地球物理學報，2000，43(3)：411-419. Dong Liangguo，Ma Zaitian，Cao Jingzhong et al.The staggered-grid high-order difference method of first-order elastic equation.Chinese Journal of Geophysics，2000，43(3)：411-419.

[12] 董良國，馬在田，曹景忠.一階彈性波方程交錯網格高階差分解法穩定性研究.地球物理學報，2000，43(6)：856-864. Dong Liangguo，Ma Zaitian，Cao Jingzhong.The stability study of the staggered-grid high-order difference method of first-order elastic equation.Chinese Journal of Geophysics，2000，43(6)：856-864.

[13] 董良國.彈性波數值模擬中的吸收邊界條件.石油地球物理勘探，1999，34(1)：46-56. Dong Liangguo.Absorptive boundary condition in elastic-wave numerical modeling.OGP，1999，34(1)：46-56.

[14] 王書強，楊頂輝，楊寬德.彈性波方程的緊致差分方法.清華大學學報(自然科學版)，2002，42(8)：1128-1131. Wang Shuqiang，Yang Dinghui，Yang Kuande.Com-pact finite difference scheme for elastic equations.Journal of Tsinghua University(Science and Techno-logy)，2002，42(8)：1128-1131.

[15] Saenger E H，Gold N，Shapiro S A.Modeling the propa-gation of elastic waves using a modified finite-difference grid.Wave Motion，2000，31(1)：77-92.

[16] 徐文才，楊國權，李振春等.橫向各向同性介質擬聲波一階速度—應力方程.石油地球物理勘探，2016，51(1)：87-96. Xu Wencai，Yang Guoquan，Li Zhenchun et al.First order velocity-stress equation in TI media.OGP，2016，51(1)：87-96.

[17] Kondoh Y，Hosaka Y，Ishii K.Kerrel optimum nearly analytical discrimination algorithm applied to parabo-lic and hyperbolic equations.Computers & Mathema-tics with Applications，1994，27(3)：59-90.

[18] 楊頂輝，滕吉文，張中杰.三分量地震波場的近似解析離散模擬技術.地球物理學報，1996，39(增刊)：283-291. Yang Dinghui，Teng Jiwen，Zhang Zhongjie.Nearly-analytical discretization modeling technique of 3-component seismic wave-fields.Chinese Journal of Geophysics，1996，39(S)：283-291.

[19] Yang D H，Teng J W，Zhang Z J et al.A nearly-analytic discrete method for acoustic and elastic wave equation.Bulletin of the Seismological Society of America，2003，93(2)：882-890.

[20] Yang D H，Lu M，Wu R S et al.An optimal nearly-ana-lytic discrete method for 2D acoustic and elastic wave equations.Bulletin of the Seismological Society of America，2004，94(5)：1982-1991.

[21] Yang D H，Song G J，Lu M.Optimally accurate nearly analytic discrete scheme for wave-field simulation in 3D anisotropic media.Bulletin of the Seismological Society of America，2007，97(5)：1557-1569.

[22] Yang D H，Song G J，Chen S et al.An improved nearly analytical discrete method：an efficient tool to simulate the seismic response of 2-D porous structures.Journal of Geophysics & Engineering，2007，4(1)：40-52.

[23] 盧明.改進的近似解析離散化方法及彈性波波場模擬[學位論文].北京：清華大學數學科學系，2004.

[24] 宋國杰.三維彈性波方程的改進近似解析離散化方法及波場模擬[學位論文].北京：清華大學數學科學系，2008.

[25] Tong P，Yang D H，Hua B L.High accuracy wave simulation-revised derivation，numerical analysis and testing of a nearly analytic integration discrete method for solving acoustic equation.International Journal of Solid and Structures，2011，48(1)：56-70.

[26] Lukacova M M，Warnecke G.Lax-Wendroff type se-cond order evolution Galerkin methods for multidimensional hyperbolic systems.East-West Journal of Numerical Mathematics，2000，8(2)：127-152.

[27] 董良國，李培明.地震波傳播數值模擬中的頻散問題.天然氣工業，2004，24(6)：53-56. Dong Liangguo，Li Peiming.Dispersive problem in seismic wave propagation numerical modeling.Natural Gas Industry，2004，24(6)：53-56.

[28] 張文生.科學計算中的偏微分方程有限差分法.北京：高等教育出版社，2006.

[29] 高正輝，孫建國，孫章慶等.基于完全匹配層構建新方法的2.5維聲波近似方程數值模擬.石油地球物理勘探，2016，51(6)：1128-1133. Gao Zhenghui，Sun Jianguo，Sun Zhangqing et al.2.5D acoustic approximate equation numerical simulation with a new construction method of the perfectly matched layer.OGP，2016，51(6)：1128-1133.

[30] 汪勇，段焱文，王婷等.優化近似解析離散化方法的二維彈性波波場分離模擬.石油地球物理勘探，2017，52(3)：458-467. Wang Yong，Duan Yanwen，Wang Ting et al.Numerical simulation of elastic wave separation in 2D isotropic medium with the optimal nearly-analytic discretization.OGP，2017，52(3)：458-467.

附錄AENADM頻散關系

將式(5)～式(7)和式(4)代入式(15),可以得到ENADM的差分格式為