李青陽(yáng) 吳國(guó)忱* 段沛然

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

0 引言

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

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

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

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

1 彈性波波動(dòng)方程

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

(1)

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

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

(2)

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

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

(3)

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

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

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

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

(4)

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

(5)

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

2.2 交錯(cuò)網(wǎng)格法

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

(6)

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

(7)

式中vx、vz為彈性波場(chǎng)x、z方向速度分量。

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

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

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

(8)

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

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

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

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

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

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

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

圖2 交錯(cuò)網(wǎng)格與準(zhǔn)規(guī)則網(wǎng)格剖分對(duì)比

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

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

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

(17)

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

(18)

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

(19a)

(19b)

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

(20)

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

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

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

(21)

(22)

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

(23)

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

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

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

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

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

圖5 縱波源波場(chǎng)水平分量(a)與垂直分量(b)對(duì)比

圖6 橫波源波場(chǎng)水平分量(a)與垂直分量(b)對(duì)比

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

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

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

5.1 層狀模型

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

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

圖8 層狀模型

(24)

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

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

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

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

5.2 Marmousi-2模型

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

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

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

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

圖11 圖10局部放大顯示

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

圖12 Mamoursi-2模型

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

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

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

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

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

6 結(jié)論

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

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

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

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

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