田雪豐

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

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

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

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

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

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

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

(1)

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

(2)

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

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

(3)

其中Δt為時(shí)間剖分步長(zhǎng)。

2 時(shí)間2M階，空間2N階差分近似

2.1 時(shí)間2M階差分近似

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

(4)

(5)

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

(6)

2.2 空間2N階差分近似

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

(7)

其中Δx為x方向的空間剖分步長(zhǎng)。

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

(8)

其中Δz為z方向的空間剖分步長(zhǎng)。

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

2.3 時(shí)間2M階，空間2N階差分格式

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

(9)

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

(10)

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

下面以時(shí)間4階、空間4階精度差分格式為例對(duì)上述過(guò)程進(jìn)行簡(jiǎn)要說(shuō)明。

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

(11)

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

(12)

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

(13)

(14)

(15)

(16)

(17)

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

(18)

(19)

(20)

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

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

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

(21)

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

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

(22)

則(6)式記為

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

(23)

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

Vn+1/2=Un

(24)

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

Wn=(Un,Vn)T

(25)

可得：

(26)

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

(27)

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

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

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

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

(28)

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

(29)

由此可見(jiàn)，求取A的特征值γ需要先求取G的特征值η。η滿(mǎn)足

|ηE-G|=0

(30)

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

(31)

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

(32)

(33)

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

(34)

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

(35)

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

(36)

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

A=PΛP-1

(37)

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

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

(38)

因此式(28)可寫(xiě)為：

|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個(gè)值中，有4個(gè)為無(wú)實(shí)部的復(fù)數(shù)，另外1個(gè)為0。當(dāng)γ=0時(shí)，β=±1滿(mǎn)足Von Neumann穩(wěn)定性條件。當(dāng)γ≠0時(shí)，令γ=ia，a為實(shí)數(shù)。若γ2+4<0，即a>2或a<-2，則β為沒(méi)有實(shí)部的復(fù)數(shù)

(43)

(44)

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

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

(45)

注意到空間差分系數(shù)Cn[N]在n為奇數(shù)時(shí)為正值，在n為偶數(shù)時(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)定性條件。需要指出的是，在同一個(gè)差分格式中，也可以將不同階數(shù)的時(shí)間導(dǎo)數(shù)用不同精度的空間差分進(jìn)行逼近，例如，在時(shí)間4階差分格式中，將1階時(shí)間導(dǎo)數(shù)利用空間4階差分格式逼近，將3階時(shí)間導(dǎo)數(shù)利用空間2階差分格式逼近，則差分格式的穩(wěn)定性條件為：

≤1 (48)

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

≤1

(49)

4 應(yīng)用實(shí)例

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

設(shè)定一個(gè)各向同性均勻介質(zhì)模型，模型各項(xiàng)參數(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階時(shí)間精度，4階空間精度差分格式進(jìn)行正演模擬。由表1可知，滿(mǎn)足穩(wěn)定性條件的最大時(shí)間步長(zhǎng)為0.001 51s。

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論與建議

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

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

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