田雪豐

(中國煤炭地質(zhì)總局地球物理勘探研究院，河北 涿州 072750)

基于有限差分法的彈性波模擬或成像處理[1-7]，受差分格式穩(wěn)定性條件的限制，每種差分格式的時間步長和空間步長的比值(簡稱時空步長之比)都被限制在一定范圍內(nèi)。為了精細(xì)地對復(fù)雜地質(zhì)模型的地震響應(yīng)進(jìn)行數(shù)值模擬，要求空間網(wǎng)格步長足夠小。因此，受限于差分格式的穩(wěn)定性要求，必須選取小的時間步長。時間步長越小，則計算的時間步數(shù)越多，計算效率越低。

基于交錯網(wǎng)格的一階彈性波方程數(shù)值求解技術(shù)[1-2,4,6,8]相比于二階彈性波方程，由于具有頻散小，收斂速度快的優(yōu)點，在彈性波的模擬和偏移中得到了廣泛應(yīng)用[8-13]。穩(wěn)定性條件是交錯網(wǎng)格差分算法的重要研究內(nèi)容[2,6,13]，Virieux[14]首先給出了三維情況下各向同性介質(zhì)中一階彈性波方程的交錯網(wǎng)格的時間2階、空間2階差分格式的穩(wěn)定性條件。Levander[15]在Virieux的基礎(chǔ)上發(fā)展了一階彈性波交錯網(wǎng)格的差分格式，提出交錯網(wǎng)格的空間差分格式可以為任意精度，并給出了時間2階精度、空間4階精度的差分格式及其穩(wěn)定性條件。此后，高階空間精度的差分格式廣泛地應(yīng)用于一階彈性波交錯網(wǎng)格。

董良國[8]通過將速度(應(yīng)力)對時間的導(dǎo)數(shù)轉(zhuǎn)化為應(yīng)力(速度)對空間的導(dǎo)數(shù)，使時間高階導(dǎo)數(shù)的計算在一個時間層內(nèi)完成，在不增加計算所需內(nèi)存的前提下，得到了任意階時間精度和空間精度的差分格式。

本文通過標(biāo)準(zhǔn)傅立葉分析得到了這種任意階時間精度的差分格式的穩(wěn)定性條件，并提出一種將不同階數(shù)的時間導(dǎo)數(shù)轉(zhuǎn)換為不同精度的空間差分的時間高階差分方法。本文的穩(wěn)定性分析結(jié)果指出：在一階彈性波交錯網(wǎng)格中，時間高階差分格式比時間2階差分格式的穩(wěn)定性條件更寬松，因而可以提高時空步長之比，減少數(shù)值模擬的時間步數(shù)。此外，本文指出，由于空間差分引起的數(shù)值頻散是由空間上的網(wǎng)格離散造成的，而時間差分引起的數(shù)值頻散是由時間上的網(wǎng)格離散造成的，兩者來源不同。因此時間精度的提高不能減小空間差分引起的數(shù)值頻散，壓制頻散時要根據(jù)頻散的來源進(jìn)行壓制。

1 一階彈性波速度應(yīng)力方程

在不考慮體力情況下，各向同性介質(zhì)一階彈性波速度-應(yīng)力方程可寫作[4,8]：

(1)

將(1)式表示為矩陣形式，有：

(2)

上式中Q為式(1)中等號右邊的五階方陣，

U(t)=[vx(t)，vz(t)，τxx(t+Δt/2)，τzz(t+Δt/2)，τxz(t+Δt/2)]T

(3)

其中Δt為時間剖分步長。

2 時間2M階，空間2N階差分近似

2.1 時間2M階差分近似

應(yīng)用交錯網(wǎng)格的思想[8,13]求解一階速度-應(yīng)力方程，根據(jù)一元函數(shù)的Taylor展式，得到速度和應(yīng)力的時間2M階差分近似，以vx和τxx分量為例，有：

(4)

(5)

由此可以將(4)(5)所代表的等式聯(lián)立為矩陣形式：

(6)

2.2 空間2N階差分近似

根據(jù)交錯網(wǎng)格思想，1階空間導(dǎo)數(shù)的2N階精度差分近似表達(dá)式如下：

(7)

其中Δx為x方向的空間剖分步長。

同理可得N方向上的1階導(dǎo)數(shù)的2N階精度差分近似為：

(8)

其中Δz為z方向的空間剖分步長。

(7)、(8)兩式中的系數(shù)Cn(N)可通過將不同節(jié)點上的Taylor展式聯(lián)立，利用待定系數(shù)法求解得到。

2.3 時間2M階，空間2N階差分格式

根據(jù)(2)式可得：

(9)

將上式代入(6)式中，得到：

(10)

將Q中的空間差分算子x和Dz用(7)、(8)兩式所示的差分格式代替，則得到時間2M階，空間2N階差分近似表達(dá)式。

下面以時間4階、空間4階精度差分格式為例對上述過程進(jìn)行簡要說明。

根據(jù)矩陣乘法，可得：

(11)

令(10)式中的時間差分精度為4階，即M為2，得到：

(12)

將(11)式代入(12)式中，以方程形式表示，得到下列各式：

(13)

(14)

(15)

(16)

(17)

將(13)—(17)中的空間微分用空間4階精度差分格式代替，得到時間4階、空間4階的一階彈性波方程差分格式?？臻g1次導(dǎo)數(shù)用(7)(8)兩式所示的差分格式進(jìn)行代替，結(jié)合微分算子的乘法運算方法，可得：

(18)

(19)

(20)

3 一階彈性波方程交錯網(wǎng)格二維空間穩(wěn)定性分析

3.1 一階彈性波速度-應(yīng)力方程Fourier分析

對(10)式進(jìn)行空間傅立葉變換，得到：

(21)

其中G(kx,kz)是Q(x,z)的傅立葉變換對。

將U(t,kx,kz)的系數(shù)矩陣記為A：

(22)

則(6)式記為

Un+1/2=Un-1/2+A·Un

(23)

式中n為時間網(wǎng)格序號，根據(jù)傅立葉分析方法，令：

Vn+1/2=Un

(24)

聯(lián)立(22)和(23)，并令：

Wn=(Un,Vn)T

(25)

可得：

(26)

(26)式中的E為單位矩陣，將式中的增長矩陣記為：

(27)

根據(jù)Von Neumann穩(wěn)定性條件可知，(10)式穩(wěn)定的條件是增長矩陣B的譜半徑不大于1，即B的模長最大的特征值的模不大于1。

3.2 差分格式譜半徑與系數(shù)矩陣A的關(guān)系

由(27)式可知，B的特征值β滿足如下關(guān)系

|βE-B|=|β2E-βA-E|=0

(28)

可見β的取值與矩陣A有關(guān)，因此先求取A的特征值γ。根據(jù)(22)式，若G的特征值為η，則：

(29)

由此可見，求取A的特征值γ需要先求取G的特征值η。η滿足

|ηE-G|=0

(30)

依據(jù)矩陣Q的表示形式，式(30)可寫為：

(31)

式中Dkx，Dkz為空間微分算子Dx和Dz的傅立葉變換對。由(4)(5)可得：

(32)

(33)

依據(jù)行列式計算規(guī)則，(31)式簡化為：

(34)

根據(jù)阿貝爾定理，對于一般的5次方程不存在用根式表達(dá)的求根公式，由(34)式可以看出，η存在一個單重根0。在不考慮0根的前提下，(34)式中可約去一個η，使η的最高階數(shù)降為4，因此，可利用4次多項式的求解方法計算得到η的5個特征值為：

(35)

由(29)式所示的γ和η的關(guān)系，得到矩陣A的5個特征值為

(36)

由于A為5階方陣，且有5個互異的特征值，因此矩陣A可對角化，有：

A=PΛP-1

(37)

式中P為可逆矩陣，Λ為A的特征值所構(gòu)成的對角陣，即：

Λ=diag(γ1，γ2，γ3，γ4，γ5)

(38)

因此式(28)可寫為：

|P(β2E-βΛ-E)P-1|=0

(39)

由于P為可逆矩陣，可以證明：

|(β2E-βΛ-E|=0

(40)

因此，B的特征值β和A的特征值γ有如下關(guān)系：

β2-βγ-1=0

(41)

解得：

(42)

由(32)、(33)、(36)式可知，γ的5個值中，有4個為無實部的復(fù)數(shù)，另外1個為0。當(dāng)γ=0時，β=±1滿足Von Neumann穩(wěn)定性條件。當(dāng)γ≠0時，令γ=ia，a為實數(shù)。若γ2+4<0，即a>2或a<-2，則β為沒有實部的復(fù)數(shù)

(43)

(44)

因此，當(dāng)|a|≤2時，式(10)式滿足Von Neumann穩(wěn)定性條件。至此，得到如下結(jié)論：當(dāng)(36)式所示的系數(shù)矩陣A的最大模長的特征值的模不大于2時，(10)式所示的差分格式穩(wěn)定。

3.3 差分格式的穩(wěn)定性條件

(45)

注意到空間差分系數(shù)Cn[N]在n為奇數(shù)時為正值，在n為偶數(shù)時為負(fù)值，因此可將(-1)n+1Cn[N]記為|Cn[N]|。由于彈性常數(shù)λ和μ均大于0，因此|γ1|>|γ3|，由此可知差分格式穩(wěn)定的條件是：

≤1

(46)

≤1

(47)

表1 各階差分精度的穩(wěn)定性條件

表1是根據(jù)(47)式得到的部分差分格式的穩(wěn)定性條件。需要指出的是，在同一個差分格式中，也可以將不同階數(shù)的時間導(dǎo)數(shù)用不同精度的空間差分進(jìn)行逼近，例如，在時間4階差分格式中，將1階時間導(dǎo)數(shù)利用空間4階差分格式逼近，將3階時間導(dǎo)數(shù)利用空間2階差分格式逼近，則差分格式的穩(wěn)定性條件為：

≤1 (48)

在時間6階精度差分格式中，將1階時間導(dǎo)數(shù)利用空間2階差分逼近，將3階時間導(dǎo)數(shù)和5階時間導(dǎo)數(shù)利用空間4階差分逼近，則差分格式的穩(wěn)定性條件為：

≤1

(49)

4 應(yīng)用實例

4.1 時間高階精度和時間2階精度差分格式正演模擬

設(shè)定一個各向同性均勻介質(zhì)模型，模型各項參數(shù)為：密度ρ=2 000kg/m3，縱波速度υp=2 000m/s，橫波速度υs=1 300m/s。設(shè)定震源表達(dá)式為：

{1-2[πf(t-100)Δt]2}e-[πf(t-100)]2

(50)

其中f為震源子波主頻，(t-100)意味著震源延遲(100Δt)ms激發(fā)。

令f=30Hz,Δh=5m，選取2階時間精度，4階空間精度差分格式進(jìn)行正演模擬。由表1可知，滿足穩(wěn)定性條件的最大時間步長為0.001 51s。

令Δt=0.001 5s，計算空間401×401個網(wǎng)格點，震源位置(201,201)，第250個采樣時刻(375ms)的vx分量波前快照如圖1(a)；令Δt=0.001 6s，記錄第250個采樣時刻(400ms)的vx分量波前快照如圖1(b)?？梢钥闯?，使用2階空間精度差分格式進(jìn)行正演， 當(dāng)Δt>Δtmax， 數(shù)值解變得不穩(wěn)定； 只要Δt<Δtmax，

圖1 不同差分階數(shù)和時間步長的波場快照Figure 1 Wave field snapshot of different difference orders and time step lengths

數(shù)值解穩(wěn)定且能準(zhǔn)確描述波前。

仍然令f=30Hz，Δh=5m。選取時間4階、空間4階精度差分格式進(jìn)行正演模擬。由表1可知，滿足穩(wěn)定性條件的最大時間步長為0.004 31s。

令Δt=0.004 4s，計算空間801×801個網(wǎng)格點，震源位置(401，401)，記錄第250個采樣時刻，即1 100ms的波前快照，如圖1(c)，可以看到，當(dāng)Δt>Δtmax時，數(shù)值解變得不穩(wěn)定。當(dāng)Δt=0.004 3s，記錄第250個采樣時刻，即1 075ms的波前快照。如圖1(d)，Δt=0.004 3s時，Δt<Δtmax，雖然波前得到了成像，但出現(xiàn)了許多不應(yīng)有的“雜波”(箭頭所指處)。這些“雜波”是由于時間步長過大所產(chǎn)生的由時間差分引起的頻散。

仍然令f=30Hz,Δh=5m。選取4階時間精度、4階空間精度差分格式進(jìn)行正演模擬。令Δt=0.003 2s，計算空間401×401個網(wǎng)格點，震源位置(201，201)，記錄t=230Δt，即736ms的波前快照，如圖1(e)。

雖然為了避免數(shù)值頻散，實際應(yīng)用中時間步長不能取到穩(wěn)定性條件條件允許的最大值。但是比起時間2階、空間4階精度差分格式穩(wěn)定性條件允許的時間步長最大值Δtmax=0.001 5s，時間4階、空間4階精度的時間步長Δt已經(jīng)可以取到其兩倍以上。

4.2 提高時間差分精度對空間差分引起的數(shù)值頻散的作用

以圖1中的模型和初始條件為例，令子波主頻f=60Hz，Δh=5m。選取2階時間精度、4階空間精度差分格式進(jìn)行正演模擬。時間剖分步長為0.001 3s，計算空間401×401個網(wǎng)格點，震源位置(201,201)，t=350Δt(455ms)的波前快照如圖2a，在子波頻率變高，波長變短，網(wǎng)格步長不變的情況下，一個波長內(nèi)的節(jié)點數(shù)減少，因此由空間采樣率不足造成了較為嚴(yán)重的頻散。

將時間差分精度提高到4階，其他條件不變，得到如圖2b的正演結(jié)果。二者的比較表明，時間差分精度的提高并沒有降低空間差分引起的數(shù)值頻散。相反，由空間差分引起的數(shù)值頻散變得更加嚴(yán)重。

保持時間差分精度為4階，將空間差分精度提高至8階，得到如圖2c的正演結(jié)果。隨著空間差分精度的提高，空間差分引起的數(shù)值頻散顯著降低。

從理論上看，由時間差分引起的數(shù)值頻散是時間采樣不足引起的，而空間差分?jǐn)?shù)值頻散是由空間采樣不足引起的[17]，因此數(shù)值頻散需要根據(jù)頻散的來源進(jìn)行壓制，不能通過提高時間(空間)差分精度來壓制空間(時間)差分引起的頻散。

圖2 不同差分階數(shù)條件下455ms時刻的波場快照Figure 2 Wave field snapshot at time 455ms underdifferent difference orders

5 結(jié)論與建議

本文在利用一階速度-應(yīng)力方程進(jìn)行彈性波正演模擬的過程中，通過將速度(應(yīng)力)對時間的導(dǎo)數(shù)轉(zhuǎn)化為應(yīng)力(速度)對空間的導(dǎo)數(shù)，使得高階時間高階導(dǎo)數(shù)的計算可以只通過一個時間層上的空間差分進(jìn)行逼近，從而降低了計算所需的內(nèi)存。同時，本文對Un+1/2=Un-1/2+A·Un類型的差分格式進(jìn)行了穩(wěn)定性分析，得到了一階彈性波方程高階時間差分格式相比一階彈性波方程2階時間差分格式的穩(wěn)定性條件寬松的結(jié)論。

時間高階差分格式允許在數(shù)值模擬過程中使用更大的時間步長，從而減少數(shù)值模擬所需的時間層數(shù)。但目前該格式存在時間層內(nèi)計算量大的問題。由于不同階數(shù)的時間導(dǎo)數(shù)可以用不同精度的空間差分進(jìn)行逼近，因此在日后的研究中，可以利用不同精度的空間差分格式進(jìn)行組合，進(jìn)一步放寬穩(wěn)定性條件的約束。

空間差分引起的數(shù)值頻散是由空間上的網(wǎng)格離散造成的，而時間差分引起的數(shù)值頻散是由時間上的網(wǎng)格離散造成的，兩者來源不同。因此壓制頻散時要根據(jù)頻散的來源進(jìn)行壓制，不能寄希望于通過提高一方的差分精度去壓制另一方的頻散。當(dāng)時間步長過大時，由時間差分引起的數(shù)值頻散會變得嚴(yán)重。