準規(guī)則網(wǎng)格高階有限差分法非均質(zhì)彈性波波場模擬

李青陽 吳國忱* 段沛然

(①中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580；②海洋國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室，山東青島 266071)

0 引言

地震波場數(shù)值模擬是已知地下介質(zhì)結(jié)構(gòu)和參數(shù)，基于射線理論或波動方程理論模擬地震波在地下傳播的一種技術(shù)。主要有有限元法、有限差分法和偽譜法等。近年來，學(xué)者們逐漸提出了一些新的數(shù)值模擬方法，如傅里葉積分法[1-2]、一步外推法[3-4]和低秩近似法[5-7]等。

有限差分法因其具有計算效率高、內(nèi)存占用小、方便實現(xiàn)等優(yōu)點，被廣泛應(yīng)用于地震波波動方程數(shù)值模擬[8-10]。Alterman等[11]最先提出層狀介質(zhì)二階彈性波方程的離散數(shù)值解，其實質(zhì)就是均勻介質(zhì)的數(shù)值解。Boore[12]提出非均勻介質(zhì)二階彈性波有限差分法，Kelly等[8]對這一方法進行了改進。Madariaga[13]提出一階速度—應(yīng)力彈性波方程交錯網(wǎng)格的有限差分法，Virieux[14-15]將交錯網(wǎng)格法用于P-SV和SH波速度—應(yīng)力方程的數(shù)值模擬。Levander[16]在Madariaga方法的基礎(chǔ)上提出四階精度交錯網(wǎng)格有限差分法，提高了數(shù)值模擬的穩(wěn)定性。Graves[17]將交錯網(wǎng)格有限差分應(yīng)用于三維彈性波介質(zhì)。Moczo等[18]推導(dǎo)了縱橫波的穩(wěn)定性條件，并指出階數(shù)越小頻散越嚴重。Kristek等[19]提出了三維空間四階速度—應(yīng)力交錯網(wǎng)格有限差分算法，有效壓制了數(shù)值頻散，提高了數(shù)值模擬精度。國內(nèi)許多學(xué)者對于有限差分法進行了大量研究，如：劉洋等[20]提出任意偶數(shù)階精度差分格式；董良國等[21]基于交錯網(wǎng)格提出時間四階、空間高階有限差分法，實現(xiàn)一階彈性波方程高精度數(shù)值求解；殷文等[22]在頻率域?qū)崿F(xiàn)了有限差分法，在提高模擬精度的同時很好地壓制了頻散；楊慶節(jié)等[23]在董良國等[21]研究的基礎(chǔ)上優(yōu)化了差分格式，提高了算法的穩(wěn)定性和精度；杜啟振等[24]聯(lián)合高階有限差分法和偽譜法共同求解聲波方程，針對復(fù)雜介質(zhì)的數(shù)值模擬有較好的穩(wěn)定性。

目前應(yīng)用最廣泛的是Virieux[15]、Levander[16]的交錯網(wǎng)格有限差分法，不僅可以對流體—固體耦合介質(zhì)進行建模，且能夠適應(yīng)于高泊松比介質(zhì)(包括流體)，同時將密度的空間分布以及變化考慮在內(nèi)，對于非均勻彈性介質(zhì)下地震波傳播建模是一個非常不錯的選擇。然而交錯網(wǎng)格有限差分法面臨內(nèi)存需求量大、計算效率低的問題[25-27]。為此，Di Bartolo等[28]提出等效交錯網(wǎng)格法，實現(xiàn)了二階變密度聲波方程的數(shù)值模擬，該方法精度與交錯網(wǎng)格的正演模擬精度一致，從數(shù)學(xué)上也能夠證明二者之間的等價性，且在存儲上的優(yōu)勢非常明顯，二維情況下內(nèi)存占用量僅相當(dāng)于交錯網(wǎng)格法的三分之一。

受Di Bartolo等[28]研究的啟發(fā)，本文應(yīng)用準規(guī)則網(wǎng)格(Quasi-regular grid，QRG)正演算法模擬非均勻彈性介質(zhì)地震波傳播。首先給出各向同性非均勻彈性介質(zhì)下通常使用的一階速度—應(yīng)力彈性波方程和二階位移彈性波方程，在特殊假設(shè)條件下分析二者與經(jīng)典二階均勻彈性波位移方程的聯(lián)系；然后回顧了QRG與交錯網(wǎng)格剖分算法，在此基礎(chǔ)上提出彈性波QRG高階有限差分算法，分析不同算法的內(nèi)存占用情況，以體現(xiàn)QRG法的優(yōu)越性；從數(shù)學(xué)角度證明了QRG法和交錯網(wǎng)格法在模擬彈性波傳播過程中的等價性，并給出不同震源的加載方式以及邊界條件和穩(wěn)定性條件； 最后采用簡單層狀模型和復(fù)雜Marmousi-2模型進行數(shù)值模擬，通過與交錯網(wǎng)格法和規(guī)則網(wǎng)格法的對比，驗證QRG法對非均勻介質(zhì)數(shù)值模擬的準確性及適用性。

1 彈性波波動方程

二維彈性波傳統(tǒng)二階位移—應(yīng)力方程為

(1)

式中：u和w分別為x、z方向的位移波場；λ和μ為彈性介質(zhì)中的拉梅常數(shù)；ρ為介質(zhì)密度；τxx和τzz表示正應(yīng)力；τxz表示切應(yīng)力；fx、fz和fxz是震源項。

將式(1)中應(yīng)力分量代入位移分量中，就得到各向同性非均勻二階位移彈性波方程

(2)

式中sx和sz為純位移方程震源項。需要注意的是，在離散情況下，加載震源方式不同會導(dǎo)致二階方程與二階位移—應(yīng)力方程模擬結(jié)果不同。

在彈性介質(zhì)均勻假設(shè)條件下，式(2)退化為經(jīng)典的二階彈性波位移方程

(3)

2 不同網(wǎng)格剖分的正演算法

2.1 規(guī)則網(wǎng)格法

規(guī)則網(wǎng)格法是地震波傳播模擬最常用的有限差分法，具體思想是對均勻假設(shè)下的彈性波方程(式(3))中的連續(xù)偏導(dǎo)數(shù)項應(yīng)用中心差分近似得到離散的近似表達式。因為均勻假設(shè)不涉及彈性參數(shù)和密度的空間變化，顯然該方法僅在假定模型參數(shù)變化忽略不計或只考慮運動學(xué)的情況下滿足精度要求，是一種利用規(guī)則網(wǎng)格剖分得到一個穩(wěn)定的解析近似解的方法。在規(guī)則網(wǎng)格法中位移被限制在空間步長為Δx、Δz和時間間隔為Δt的每個離散點上，然后通過泰勒級數(shù)展開，并用已知的中心差分算子來近似彈性波方程的時空連續(xù)偏導(dǎo)項。

以位移水平分量為例，時間二階導(dǎo)數(shù)差分近似為

(4)

對于空間導(dǎo)數(shù)，二階有限差分近似為

(5)

式中假設(shè)Δx=Δz=h。從式(5)可知，在同一方向的二階導(dǎo)數(shù)融合了半網(wǎng)格的思想，空間步長是h，而混合求導(dǎo)項的空間步長必須是2h。

2.2 交錯網(wǎng)格法

Virieux[15]的交錯網(wǎng)格法是通過建立半網(wǎng)格進行中心差分運算，實現(xiàn)對非均勻彈性波方程進行數(shù)值求解。例如空間二階精度差分格式可表示為

(6)

類似地，對于時間的中心差分運算也采用半網(wǎng)格思想，具體的遞推流程如圖1所示。經(jīng)典交錯網(wǎng)格提出后最先是應(yīng)用于彈性波一階速度—應(yīng)力方程

(7)

式中vx、vz為彈性波場x、z方向速度分量。

圖1 彈性介質(zhì)交錯網(wǎng)格模擬示意圖

2.3 準規(guī)則網(wǎng)格法

受Di Bartolo等[28]等效交錯網(wǎng)格法啟發(fā)，本文提出QRG高階有限差分法，利用交錯網(wǎng)格對位移分量剖分，實現(xiàn)非均勻介質(zhì)彈性波場數(shù)值模擬。該方法是建立在二階位移方程之上，實際上該方程與一階速度—應(yīng)力方程是等價的，等價的基礎(chǔ)是位移與速度滿足

(8)

將式(8)代入一階速度—應(yīng)力方程(式(7))，則

轉(zhuǎn)換為二階位移—應(yīng)力方程(式(1))。雖然二階位移—應(yīng)力方程的時間離散并沒有落在半網(wǎng)格上，即無半時間網(wǎng)格，但是空間離散仍然采用交錯網(wǎng)格法剖分。

圖2為QRG法與交錯網(wǎng)格法網(wǎng)格剖分對比，可以看出，在二維情況下，交錯網(wǎng)格法需要存儲應(yīng)力(τxx、τzz、τxz)、位移(u、w)共5個變量場，而QRG法只需存儲2個位移場(u、w)；在三維情況下，交錯網(wǎng)格法需要存儲應(yīng)力(τxx、τyy、τzz、τxy、τyz、τxz)、位移(u、v、w)共9個變量場，而QRG法只需存儲3個位移場(u、v、w)，因此QRG法內(nèi)存占用更少。

QRG法剖分下的任意階差分格式為

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

式中：al和am是差分系數(shù)(表1)； 2N為差分階次；b為空間導(dǎo)數(shù)項對應(yīng)的拉梅常數(shù)。

式(9)～式(16)中位移分量的下標是嚴格按照圖2中準規(guī)則網(wǎng)格剖分方法給出的；式(9)～式(12)對應(yīng)非均勻彈性波位移水平分量，其下標與正三角形一一對應(yīng)；式(13)～式(16)對應(yīng)位移垂直分量，其下標對應(yīng)于倒三角形。

圖2 交錯網(wǎng)格與準規(guī)則網(wǎng)格剖分對比

表1 任意階差分系數(shù)

3 QRG與交錯網(wǎng)格法等價性證明

以彈性波方程的水平分量為例，式(2a)忽略震源項可表示為

(17)

為方便書寫，選用空間二階精度近似，則式(17)QRG差分格式為

(18)

二階位移—應(yīng)力方程(式(1))交錯網(wǎng)格差分格式為

(19a)

(19b)

將式(19)中的應(yīng)力項代入水平位移項中，整理可得

(20)

對式(20)應(yīng)用加法交換律即得到與QRG法差分格式(式(18))完全相同的結(jié)果，上述結(jié)果在空間任意階差分格式、推廣至三維的差分格式均成立。垂直分量的推導(dǎo)過程與之類似，不再贅述。

4 震源、穩(wěn)定性和邊界條件

QRG法震源的加載需要考慮其網(wǎng)格剖分的特殊性。二階位移—應(yīng)力方程整理為二階方程的過程中，對二階位移—應(yīng)力方程所加載的震源進行了一次空間求導(dǎo)運算，以位移水平分量為例

(21)

(22)

通常正演模擬加載的是點震源，在空間坐標可以表達為f×δ(i，k)，則震源的空間導(dǎo)數(shù)為

(23)

由上式可知，對f作中心差分求導(dǎo)為零值，因此在數(shù)值模擬中，二階純位移方程無法得到與位移—應(yīng)力在應(yīng)力上加載震源模擬的結(jié)果。鑒于此，選擇將震源加載在位移分量上，其物理意義可以這樣描述：對空間一點加力，使該點位移產(chǎn)生響應(yīng)并等效于力源，位移響應(yīng)即為新的震源加載在該點，這樣就避免震源對空間求導(dǎo)。因此，規(guī)則網(wǎng)格法、QRG法和交錯網(wǎng)格法同時加載位移震源模擬結(jié)果是相同的。

對于經(jīng)典規(guī)則網(wǎng)格法的位移加載震源，通常有四點法和八點法，而QRG法因為網(wǎng)格剖分不同，所以震源加載的方式也會有所不同。以四點法為例，圖3和圖4是兩種方法縱波震源、橫波震源的加載方式對比示意圖，規(guī)則網(wǎng)格法震源中虛線表示水平和垂直方向加載位移的合成，由于QRG法的位移是相互交錯排列，因此加載方式與規(guī)則網(wǎng)格不同。

圖3 規(guī)則網(wǎng)格法(左)與準規(guī)則網(wǎng)格法(右)縱波震源加載示意圖

圖5和圖6分別是QRG法不同的震源波場。因為QRG法的震源加載在位移上，所以縱波/橫波震源加載過程中，會因為網(wǎng)格的離散導(dǎo)致殘留的橫波/縱波存在，該問題在規(guī)則網(wǎng)格法震源加載中同樣存在。對于這種情況可以采用二階位移—應(yīng)力方程和二階純位移方程相互耦合的方法[31]在應(yīng)力上加載震源，得到干凈的縱波/橫波波場，如圖7所示。

圖4 規(guī)則網(wǎng)格法(左)與準規(guī)則網(wǎng)格法(右)橫波震源加載示意圖

圖5 縱波源波場水平分量(a)與垂直分量(b)對比

圖6 橫波源波場水平分量(a)與垂直分量(b)對比

之前證明了QRG法與交錯網(wǎng)格法的等價性，因此QRG法的穩(wěn)定性條件與交錯網(wǎng)格法一致[21]。對于有限模型空間造成的人工反射問題，QRG法與規(guī)則網(wǎng)格法加載方式類似，通常會采用吸收邊界條件來模擬無限介質(zhì)，常見的方法有波動方程分解法[32]、旁軸近似法[33]、阻尼衰減法[34]、完全匹配層吸收邊界[35]，本文采用完全匹配層吸收邊界條件。

圖7 二階位移—應(yīng)力方程和二階

5 數(shù)值模擬結(jié)果及分析

5.1 層狀模型

首先通過簡單模型驗證QRG法的正演精度。層狀模型(圖8)尺寸為1500m×1000m，其中第二層與第三層縱、橫波速度相同，但密度相差很大；震源是主頻為20Hz雷克子波，位于(750m，0m)處，空間步長統(tǒng)一為5m，時間采樣間隔為0.4ms。

圖9為層狀介質(zhì)模型交錯網(wǎng)格法、QRG法和規(guī)則網(wǎng)格法模擬的0.4s波場快照，可以看出，交錯網(wǎng)格法和QRG法對于密度異常界面能夠產(chǎn)生反射、透射波，而規(guī)則網(wǎng)格法模擬忽略了密度參數(shù)信息(式(3))，不能有效識別密度界面，表明交錯網(wǎng)格法和QRG法有更高的模型參數(shù)利用率。圖10為抽取接收點位于x=500m和x=1000m的水平分量地震記錄，可見：三種方法模擬記錄中直達波(0.2s)與一次反射波(0.4～1.6s)的旅行時信息基本一致，說明模型參數(shù)的空間變化對于波的旅行時影響微乎其微；QRG法和交錯網(wǎng)格法相位信息基本吻合，振幅略微有些許差異，是由計算過程中誤差累積造成的；二者與規(guī)則網(wǎng)格對比差異很大，主要體現(xiàn)在反射波的相位與振幅信息，是介質(zhì)的非均勻性造成的。對式(2)進行簡單數(shù)學(xué)變換可得

圖8 層狀模型

(24)

圖9 層狀介質(zhì)模型不同方法模擬的0.4s波場快照

對比上式與式(3)可知，非均勻介質(zhì)彈性波方程位移分量比均勻介質(zhì)下方程多出了關(guān)于拉梅常數(shù)的導(dǎo)數(shù)項Au和Aw，直達波是從炮點沿第一層介質(zhì)直接傳播至檢波點，該過程滿足介質(zhì)均勻假設(shè)，即拉梅常數(shù)導(dǎo)數(shù)為0，因此QRG法的直達波模擬結(jié)果與規(guī)則網(wǎng)格法結(jié)果一致。在層狀模型中，介質(zhì)的非均勻性主要體現(xiàn)在在彈性間斷面上，而在層間介質(zhì)仍然滿足均勻假設(shè)，因此Au和Aw僅對界面處的反射、透射波的振幅和相位產(chǎn)生影響，而不改變旅行時信息。

圖11為圖10的局部放大顯示(0.3～0.8s)，可見0.35～0.45s處第一界面的交錯網(wǎng)格法和QRG法的反射波旅行時、相位和振幅基本一致，在0.5～0.6s處第二界面的反射波旅行時和相位基本吻合、振幅有所差異。而規(guī)則網(wǎng)格法除了旅行時信息與前二者吻合，相位和振幅均偏小。在0.6～0.7s處交錯網(wǎng)格法和QRG法可見純密度界面的反射波信息，而在規(guī)則網(wǎng)格法記錄中觀察不到。

綜上所述，相比于規(guī)則網(wǎng)格法，QRG法能夠較準確地模擬非均勻彈性介質(zhì)下的地震波傳播，其旅行時、相位及振幅信息均與交錯網(wǎng)格法模擬結(jié)果基本一致，在兼具模擬精度和節(jié)省內(nèi)存的考慮下，QRG法在振幅上的微弱差異完全可以接受。

5.2 Marmousi-2模型

采用復(fù)雜Marmousi-2模型(圖12)驗證本文方法對復(fù)雜模型的適用性和穩(wěn)定性。模型尺寸為 6800m×1400m，空間步長為5m；震源采用主頻為20Hz的雷克子波，時間采樣間隔為0.4ms，震源位于(3400m，0m)處；檢波器均勻布置在模型表面，間隔為5m。

圖13為Marmousi-2模型采用QRG法和交錯網(wǎng)格法模擬的0.6s時刻水平分量波場快照，圖14為垂直分量波場快照，兩種方法的時空差分精度均為時間二階和空間十階。波場信息完整且無明顯頻散，證明QRG法對復(fù)雜模型數(shù)值模擬穩(wěn)定性良好，從波場快照可見本文QRG法對復(fù)雜模型模擬結(jié)果非常穩(wěn)定。

圖15為兩種方法模擬的水平分量道集，圖16為兩種方法模擬的垂直分量道集，從兩種方法模擬記錄結(jié)果可見，QRG與交錯網(wǎng)格法的模擬精度基本一致。圖17為兩種方法模擬結(jié)果的單道水平及垂直分量對比，可見對復(fù)雜的Marmousi-2模型，兩種方法模擬的地震信號旅行時、相位和波形信息基本一致。

圖10 層狀模型不同方法模擬的水平分量地震記錄對比

圖11 圖10局部放大顯示

綜上所述，本文提出QRG法高階有限差分法在非均勻彈性介質(zhì)復(fù)雜模型中的展示出良好的穩(wěn)定性，同時通過與交錯網(wǎng)格法結(jié)果的對比，體現(xiàn)了QRG法所具有的優(yōu)勢。

圖12 Mamoursi-2模型

圖13 Mamoursi-2兩種方法模擬的0.6s水平分量波場快照

圖14 Mamoursi-2兩種方法模擬的0.6s垂直分量波場快照

圖15 Mamoursi-2兩種方法模擬水平分量的道集

圖16 Mamoursi-2兩種方法模擬垂直分量的道集

圖17 Marmousi-2兩種方法模擬的水平(a)及垂直(b)分量記錄對比

6 結(jié)論

針對非均勻彈性介質(zhì)中地震波傳播正演模擬問題，本文提出一種新的QRG方法，對彈性波方程不同的位移分量采用交錯的思想排列，并對位移分量做新的中心差分運算。通過算法分析和模型測試得到以下認識：

(1)本文提出的QRG法與交錯網(wǎng)格法具有相同的模擬精度，并且優(yōu)于傳統(tǒng)的規(guī)則法，同時具有良好的穩(wěn)定性；

(2)QRG法在內(nèi)存方面具有很大的優(yōu)勢，對同一個二維模型分別用QRG法和交錯網(wǎng)格法做一次正演模擬，QRG法的內(nèi)存使用量比交錯法減少60%，如推廣到三維，則減少66.7%；

(3)QRG法加載縱/橫波震源會產(chǎn)生與規(guī)則網(wǎng)格法一樣微弱的橫/縱波殘留，在一定誤差允許范圍內(nèi)是可以忽略不計的。如果要得到干凈的波場，可以通過交錯網(wǎng)格法與QRG法相互耦合的方式解決這一問題。

QRG法計算量較大，如何引入一些優(yōu)化系數(shù)方法降低計算復(fù)雜度，是今后的研究方向； 另外，如何將QRG法應(yīng)用到逆時偏移、最小二乘偏移以及全波形反演以提高地震偏移成像與反演的精度，也是今后的研究方向。