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

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

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

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

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

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

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

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

擴展最優(yōu)近似解析離散化 彈性波方程 數(shù)值模擬 穩(wěn)定性條件 數(shù)值頻散

1 引言

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

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

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

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

二維聲波波動方程為

(1)

設(shè)滿足聲波方程的解為U，定義

(2)

(3)

式中m、k分別表示S對x和t的偏導(dǎo)數(shù)的階次。

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

(4)

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

(5)

(6)

(7)

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

(8)

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

(9)

(10)

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

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

(11)

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

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

(12)

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

(13)

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

INADM和ONADM方法有各自的優(yōu)缺點，ONADM的精度比INADM高，但INADM具有更廣泛的適用性，能夠?qū)︷椊橘|(zhì)等復(fù)雜介質(zhì)波動方程進行數(shù)值模擬。所以本文選擇對INADM方法進行修正或擴展，提高INADM中位移和速度項的差分精度是本文的研究目的。

(14)

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

(15)

(16)

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

3 ENADM方法分析

3.1 精度分析

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

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

(17)

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

(18)

LWC差分格式[26]為

(19)

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

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

3.2 頻散分析

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

ENADM方法頻散關(guān)系(推導(dǎo)過程見附錄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)

確定φ后，代入上述頻散關(guān)系，可以解得對應(yīng)的ω。定義數(shù)值波速與真實速度的比值為

(22)

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

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

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

(23)

ONADM頻散關(guān)系為

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

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

(24)

LWC頻散關(guān)系為

cosω-1

(25)

SFDM頻散關(guān)系[27]為

(26)

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

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

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

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

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

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

3.3 穩(wěn)定性分析

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

(27)

(28)

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

3.4 計算效率分析

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

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

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

4 模型試算

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

(29)

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

(30)

4.1 均勻介質(zhì)模型

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

圖5 均勻介質(zhì)模型彈性波方程ENADM數(shù)值模擬波場快照

4.2 水平層狀介質(zhì)模型

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

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

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

5 結(jié)論與建議

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

[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] 董良國，馬在田，曹景忠等.一階彈性波方程交錯網(wǎng)格高階差分解法.地球物理學(xué)報，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] 董良國，馬在田，曹景忠.一階彈性波方程交錯網(wǎng)格高階差分解法穩(wěn)定性研究.地球物理學(xué)報，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] 董良國.彈性波數(shù)值模擬中的吸收邊界條件.石油地球物理勘探，1999，34(1)：46-56. Dong Liangguo.Absorptive boundary condition in elastic-wave numerical modeling.OGP，1999，34(1)：46-56.

[14] 王書強，楊頂輝，楊寬德.彈性波方程的緊致差分方法.清華大學(xué)學(xué)報(自然科學(xué)版)，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] 徐文才，楊國權(quán)，李振春等.橫向各向同性介質(zhì)擬聲波一階速度—應(yīng)力方程.石油地球物理勘探，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] 楊頂輝，滕吉文，張中杰.三分量地震波場的近似解析離散模擬技術(shù).地球物理學(xué)報，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] 盧明.改進的近似解析離散化方法及彈性波波場模擬[學(xué)位論文].北京：清華大學(xué)數(shù)學(xué)科學(xué)系，2004.

[24] 宋國杰.三維彈性波方程的改進近似解析離散化方法及波場模擬[學(xué)位論文].北京：清華大學(xué)數(shù)學(xué)科學(xué)系，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] 董良國，李培明.地震波傳播數(shù)值模擬中的頻散問題.天然氣工業(yè)，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] 張文生.科學(xué)計算中的偏微分方程有限差分法.北京：高等教育出版社，2006.

[29] 高正輝，孫建國，孫章慶等.基于完全匹配層構(gòu)建新方法的2.5維聲波近似方程數(shù)值模擬.石油地球物理勘探，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] 汪勇，段焱文，王婷等.優(yōu)化近似解析離散化方法的二維彈性波波場分離模擬.石油地球物理勘探，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頻散關(guān)系

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

(A-1)

(A-2)

同樣地，可以求得u和s的各階空間偏導(dǎo)數(shù)表達式，有

(A-3)

(A-4)

(A-5)

(A-6)

(A-7)

(A-8)

(A-9)

(A-10)

(A-11)

式(A-1)右邊的第3項可表示為

Δxkxsin(Δxkx)+4cos(Δxkz)+

(A-12)

第4項可表示為

(A-13)

第5項可表示為

(A-14)

第6項可表示為

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

(A-15)

將式(A-12)～式(A-15)代入式(A-1)可以得到ENADM差分格式的頻散關(guān)系為

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

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

2(cosφxcosφz-sinφx-sinφz+1)}

(A-16)

附錄B二維ENADM算法的增長矩陣

ENADM算法的增長矩陣G6×6各元素表達式為

g11= 1+(2α2-α4)(cosθ1+cosθ2-2)+

sin(-θ1+θ2)-4sinθ1-6sinθ2]+

sin(-θ1+θ2)-4sinθ1-6sinθ2]+

6sinθ1]

6sinθ1]

-cosθ1-2]

cos(θ1-θ2)-2cosθ2-2cosθ1+2]

g42=-iαvsinθ1+i2α3vsinθ1

g43=-iαvsinθ2+i2α3vsinθ2

g44= 1+2α2(cosθ1+cosθ2-2)-α4(cosθ1+

2cosθ2-2cosθ1+2]

式中：θ1=Δxkx;θ2=Δxkz。

附錄C二維彈性波方程位移的時間高階導(dǎo)數(shù)

(本文編輯：宜明理)

汪勇 博士，副教授，1979年生；2001年獲湖北大學(xué)計算機軟件專業(yè)工學(xué)學(xué)士學(xué)位；2009年獲中國地質(zhì)大學(xué)(武漢)地球探測與信息技術(shù)專業(yè)工學(xué)碩士學(xué)位；2013年獲中國地質(zhì)大學(xué)(武漢)地球探測與信息技術(shù)專業(yè)工學(xué)博士學(xué)位；現(xiàn)在長江大學(xué)地球物理與石油資源學(xué)院主要從事地震波場理論和油氣儲層預(yù)測方面的研究。

1000-7210(2017)05-0928-13

P631

A

10.13810/j.cnki.issn.1000-7210.2017.05.005

*湖北省武漢市蔡甸區(qū)蔡甸街大學(xué)路111號長江大學(xué)地球物理與石油資源學(xué)院，430100。Email：cdwangyong@yangtzeu.edu.cn

本文于2016年10月18日收到，最終修改稿于2017年8月1日收到。

本項研究受國家“973”計劃項目(2013CB228605)和中國石油科技創(chuàng)新基金項目(2015D-5006-0301)聯(lián)合資助。