基于自適應(yīng)網(wǎng)格的仿真型有限差分地震波數(shù)值模擬

劉志強(qiáng)， 孫建國， 孫輝， 劉明忱， 高正輝， 石秀林

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026

?

基于自適應(yīng)網(wǎng)格的仿真型有限差分地震波數(shù)值模擬

劉志強(qiáng)， 孫建國*， 孫輝， 劉明忱， 高正輝， 石秀林

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026

在復(fù)雜山地和復(fù)雜海底條件下，地表和海底的劇烈起伏對地震波數(shù)值模擬提出了更高的要求.常規(guī)有限差分法采用矩形網(wǎng)格對模型進(jìn)行網(wǎng)格剖分，由于矩形網(wǎng)格自身的限制，起伏地表或起伏海底只能由一系列階梯狀折線代替，從而引起人為虛假繞射波.此外，在模擬液-固界面的反射波時，如果界面與網(wǎng)格線不一致，則需要更密的網(wǎng)格才能得到精確的結(jié)果.為了解決上述問題，本文將自適應(yīng)網(wǎng)格生成技術(shù)引入到起伏海底速度模型的網(wǎng)格剖分中，采用高階仿真型有限差分法(MFD)對曲線坐標(biāo)下的聲波方程波進(jìn)行了數(shù)值模擬.利用自適應(yīng)網(wǎng)格生成技術(shù)對速度模型進(jìn)行網(wǎng)格剖分不僅可以準(zhǔn)確地描述模型邊界，而且可以有效消除虛假繞射波.高階仿真型有限差分法可以有效壓制頻散提高計(jì)算精度.模型試算結(jié)果表明，本文方法對復(fù)雜海底模型具有很好的適應(yīng)性.

虛假繞射； 數(shù)值頻散； 自適應(yīng)網(wǎng)格； 高階仿真型有限差分法

1 引言

地震波數(shù)值模擬中，首先要做的就是對模型進(jìn)行網(wǎng)格剖分，網(wǎng)格質(zhì)量對方程求解的精度和效率都起著至關(guān)重要的作用.傳統(tǒng)意義下好的網(wǎng)格性質(zhì)，如大小均勻、長寬適當(dāng)、相互正交等，有時并不能獲得精確的數(shù)值解.人們通過研究證實(shí)，依據(jù)物理問題中參數(shù)的變化來構(gòu)造網(wǎng)格，往往能夠獲得與物理問題更相適應(yīng)、更高精度的數(shù)值解.用變分法來構(gòu)造網(wǎng)格的思想最早由Barfield(1970)提出.Brackbill和Saltzman(1982)基于Barfield(1970)的思想，給出了描述網(wǎng)格疏密性、正交性和光滑性的泛函.Antonios(1988)對Brackbill和Saltzman(1982)的方法進(jìn)行了改進(jìn)，提出了一種帶方向控制的泛函來控制網(wǎng)格走向的方法.在國內(nèi)也有很多人對該方法進(jìn)行了相關(guān)研究.康紅文等(2002)對自適應(yīng)網(wǎng)格中的權(quán)函數(shù)進(jìn)行了研究.于明(2004)在前人的基礎(chǔ)上改進(jìn)了光滑性控制的變分表達(dá)式，使其與疏密程度和正交性控制方程具有一致的量綱，并且添加了網(wǎng)格正規(guī)性的控制方程.

滿足地質(zhì)模型要求的曲線網(wǎng)格生成后，就需要在這些離散點(diǎn)上計(jì)算波場.有限差分法簡單靈活，是地球物理中最常用的數(shù)值計(jì)算方法.有限差分法地震波數(shù)值模擬由Alterman和Karal(1968)提出.Boore(1972)和Kelly等(1976)將該方法推廣至非均勻介質(zhì).但是他們的有限差分法存在計(jì)算精度低，容易導(dǎo)致計(jì)算不穩(wěn)定等問題.Madariaga(1976)首次提出了在交錯網(wǎng)格中求解一階速度-應(yīng)力彈性波方程的算法來提高有限差分法的精度和穩(wěn)定性.Virieux(1984)將該方法推廣到非均勻介質(zhì)地震波數(shù)值模擬中.Crase(1990)給出了基于交錯網(wǎng)格系統(tǒng)的高階差分格式.傳統(tǒng)有限差分法采用矩形網(wǎng)格對模型進(jìn)行網(wǎng)格剖分，由于矩形網(wǎng)格自身的限制，起伏的海底界面只能由一系列的階梯狀折線代替，從而造成數(shù)值模擬過程中的人為虛假繞射波.此外，在模擬液-固構(gòu)造下的地震波場時，如果界面與網(wǎng)格線不一致，則在界面處需要更細(xì)的網(wǎng)格才能得到精確的結(jié)果(van Vossen et al.，2002).因此，對于傳統(tǒng)的規(guī)則網(wǎng)格有限差分法需要對整個計(jì)算區(qū)域采用較密的網(wǎng)格，而這不僅增加了內(nèi)存需求而且降低了計(jì)算效率.為了解決上述問題，Moczo(1989)提出了可變網(wǎng)格的思想；Jastram和Behle(1992)將變網(wǎng)格的思想運(yùn)用到聲波方程數(shù)值模擬中；Jastram和Tessmer(1994)將該思想運(yùn)用到彈性波數(shù)值模擬中，提出了縱向網(wǎng)格步長逐漸變化的算法；朱生旺等(2007)將該算法的空間差分精度進(jìn)一步提高，使該算法的模擬精度和靈活性更高.但由于它們依舊是基于規(guī)則網(wǎng)格實(shí)現(xiàn)的，并不能從根本上消除階梯狀網(wǎng)格引起的虛假繞射波.為此人們發(fā)展了不規(guī)則網(wǎng)格有限差分法.Opr?al和Zahradniík(1999)給出了非規(guī)則介質(zhì)中的二階波動方程的矩形不規(guī)則網(wǎng)格差分方法.Pitarka(1999)提出了各項(xiàng)同性介質(zhì)中矩形不規(guī)則網(wǎng)格的有限差分方法.張劍鋒(1998)通過坐標(biāo)映射得到了非規(guī)則網(wǎng)格差分方法.但是他們的方法計(jì)算精度較低.

為了有效消除人為虛假繞射波，本文將基于變分原理的自適應(yīng)網(wǎng)格生成技術(shù)引入到起伏海底速度模型的網(wǎng)格剖分中，采用曲線坐標(biāo)系下的高階仿真型有限差分法對聲波方程進(jìn)行了數(shù)值模擬.該方法生成的網(wǎng)格具有以下特點(diǎn):(1) 可以更加準(zhǔn)確地描述模型界面; (2) 可以盡量保證網(wǎng)格的光滑性和正交性; (3) 網(wǎng)格剖分更加合理.高階仿真型有限差分法可以有效壓制數(shù)值頻散、提高計(jì)算精度.

2 方法原理

2.1 自適應(yīng)網(wǎng)格生成

考慮任意形狀的物理域(x,z)通過坐標(biāo)變換可以轉(zhuǎn)換到直角四邊形的計(jì)算域(ξ,η)(如圖1)，那么存在坐標(biāo)系變換的Jacobi式J=xξzη-xηzξ，并且成立坐標(biāo)變換

ξx=zη/J,ξz=-xη/J,ηx=-zξ/J,ηz=xξ/J.

(1)

圖1 網(wǎng)格映射示意圖Fig.1 The sketch map of grid mapping

根據(jù)微分幾何理論，J表示兩個坐標(biāo)系中面積微元之比，即dxdz=|J|dξdη.|J|的大小能反映網(wǎng)格疏密程度，可以用J和一個用來調(diào)節(jié)網(wǎng)格疏密的權(quán)函數(shù)w(x,z)的乘積來調(diào)節(jié)整個物理域網(wǎng)格的疏密性(Brackbill，1982)：

(2)

式中w(x,z)為加權(quán)函數(shù)，D為整個模型區(qū)域.

定義表示映射后坐標(biāo)光滑性的泛函：

(3)

定義表示網(wǎng)格正交性的泛函：

(4)

通常，生成自適應(yīng)網(wǎng)格時需要同時考慮所有性質(zhì).故加權(quán)組合上述泛函取極小值時的Euler方程展開式，可得(見附錄A)：

(5)

式中

ai=λvavi+λsasi+λoaoi，

bi=λvbvi+λsbsi+λoboi，

ci=λvcvi+λscsi+λocoi， i=1,2,3

λv,λs,λo為正的加權(quán)系數(shù)(可以根據(jù)問題的需要來選擇數(shù)值大小).

式(5)即為最終得到的自適應(yīng)網(wǎng)格生成控制方程，可以采用Gauss-Seidel法來求解.權(quán)函數(shù)w(x,z)可以根據(jù)具體問題選取,本文以速度模型的梯度來構(gòu)造權(quán)函數(shù).可令權(quán)函數(shù)為

(6)

式中ε為用來調(diào)整權(quán)函數(shù)的大小的系數(shù)，ux是速度的橫向梯度，uz是速度的縱向梯度.

2.2 曲線坐標(biāo)系下的高階仿真型有限差分

設(shè)物理平面(x,z)和計(jì)算平面(ξ,η)存在如下的一一對應(yīng)關(guān)系:

(7)

根據(jù)鏈?zhǔn)椒▌t，變量u的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)可表示成如下形式：

(8)

利用這些導(dǎo)數(shù)可以求得梯度(GRID)、散度(DIV)和拉普拉斯(LAP)算子的差分近似，但是該方法不僅方程復(fù)雜而且很難保證解的物理性質(zhì).因此我們用映射法和支撐算子法來建立不規(guī)則坐標(biāo)系下的DIV和GRID的仿真型有限差分近似.

定義離散函數(shù)空間Hξ和Hη來表示向量A=(AX,AZ)的水平分量和垂直分量；定義離散空間HC表示標(biāo)量U.它們在笛卡爾坐標(biāo)系下的分布如圖2所示.

圖2 離散函數(shù)的空間分布Fig.2 The spatial distribution of the discrete function

為了獲得二維曲線坐標(biāo)系下DIV和GRAD差分近似，我們需要一維差分算子和映射算子.?/? ξ的差分算子可分為如下兩種類型:

Dξ∶Hξ→HC;

則四階差分算子的具體形式如下(Castillo et al,2000)：

(9)

式中Aξ表示向量A的水平分量，U表示標(biāo)量.

具體形式如下：

(10)

式中Aη表示向量A的垂直分量，U表示標(biāo)量.

定義映射算子：

具體形式如下：

(11)

(12)

(13)

式中AX,AZ是向量A的水平分量和垂直分量，

).

利用支撐算子法得到GRID的離散近似(見附錄B).

(14)

(15)

式中GXi,j+1/2,GZi+1/2,j為梯度的水平分量和垂直分量.

3 曲線坐標(biāo)系下的聲波方程數(shù)值求解

非均勻各項(xiàng)同性介質(zhì)中的一階聲波方程為

(16)

其中，p為壓強(qiáng)，vx，vz為質(zhì)點(diǎn)速度，ρ是密度，vp是縱波速度.它們在笛卡爾坐標(biāo)系下的分布如圖3.

圖3 一階聲波方程的變量分布Fig.3 The variable distribution of first-order acoustic equation

利用四階仿真型有限差分求解方程(16)中的空間導(dǎo)數(shù)得

(17)

把方程(16)改寫成如下形式：

(18)

式中U=(p,vx,vz)，F(xiàn)是關(guān)于波場變量U的函數(shù).

令Un=U(nΔt)，則可用p步顯式Runge-Kutta算法對方程(19)進(jìn)行迭代求解(Bogey，2004)：

(19)

其中，αs是Runge-Kutta算法的系數(shù)，Δt是時間步長.

4 模型試算

下面給出五個模型算例，用以驗(yàn)證本文方法的精度、有效性和適用性.第一個是均勻半空間速度模型，通過與解析解對比，用來驗(yàn)證本文方法的精度.第二和第三個是起伏海底兩層速度模型，通過與規(guī)則網(wǎng)格有限差法對比，來驗(yàn)證方法的有效性.第四個和第五個是復(fù)雜海底速度模型，用來驗(yàn)證本文方法的適應(yīng)性.

4.1 均勻速度模型

首先利用本文方法對均勻半空間速度模型進(jìn)行模擬.模型速度為2000 m·s-1，密度為1000 kg·m-3.數(shù)值解所用空間網(wǎng)格大小為10 m，時間采樣間隔為1 ms,震源是主頻為30 Hz的雷克子波，位于(20 m,20 m)處，檢波點(diǎn)位于(120 m,10 m).圖4a為相應(yīng)的波形對比，圖中Analysis曲線為Cagniard-De Hoop算法得到的解析解，MFD曲線為采用本文方法得到的數(shù)值解.由圖可見,本文方法得到的直達(dá)波與解析解基本吻合.圖4b為直達(dá)波振幅譜對比.由圖可見，本文方法得到的直達(dá)波振幅譜和解析解的振幅譜基本一致.

4.2 海底背斜速度模型

圖5a為海底單背斜速度模型，模型大小為2000 m×2000 m,第一層速度為1500 m·s-1，第二層速度為2000 m·s-1，密度為1000 m·s-1.計(jì)算域中的空間網(wǎng)格尺寸為10 m，時間采樣間隔為1 ms，最大記錄時間為1.6 s,震源采用主頻為30 Hz的雷克子波，位于(1000 m,20 m)處，檢波點(diǎn)深度為10 m.圖5b為利用本文方法得到的自適應(yīng)網(wǎng)格示意圖，由圖可以看出，利用本文方法對海底背斜模型進(jìn)行自適應(yīng)網(wǎng)格剖分不僅可以更準(zhǔn)確描述海底界面,而且在界面附近有很好的密集性、光滑性和正交性.圖6為波場快照對比、圖7為模擬地震記錄對比，由圖可以看出，相對與規(guī)則網(wǎng)格有限差分法，本文方法可以有效消除階梯狀網(wǎng)格引起的人為虛假繞射波.圖8為不同方法得到的反射波最大振幅與解析解的誤差對比；實(shí)線表示利用本文方法得到的反射波最大振幅誤差，虛線表示規(guī)則網(wǎng)格有限差分法得到的反射波最大振幅誤差.由圖可以看出，相對與規(guī)則網(wǎng)格有限差分法，本文方法在模擬起伏海底反射波時具有更高的計(jì)算精度.

圖4 均勻半空間速度模型中解析解和MFD得到的數(shù)值解對比(a) 波形對比; (b) 振幅譜對比.Fig.4 The contract of analytical and MFD solution in the homogeneous half-space velocity model(a) The contract of waveform; (b) The contract of amplitude spectrum.

圖5 海底單背斜模型(a)和自適應(yīng)網(wǎng)格示意圖(b)Fig.5 The single anticline model under the sea (a) and the sketch map of adaptive grid (b)

圖6 0.8 s時的波場快照對比(a) 本文方法； (b) 規(guī)則網(wǎng)格有限差分法.Fig.6 The contract of wavefield snapshot at 0.8 s(a) Our method； (b) Regular grid finite-difference method.

4.3 起伏海底速度模型

圖9a為起伏海底速度模型，模型大小為6000 m×4000 m，第一層速度為1500 m·s-1，第二層速度為2000 m·s-1，密度為1000 m·s-1.計(jì)算域中的空間網(wǎng)格尺寸為10 m，時間采樣間隔為1 ms，最大記錄時間為6 s,震源是主頻為30 Hz的雷克子波，位于(1000 m,20 m)處，檢波點(diǎn)深度為10 m.圖9b為利用本文方法得到的自適應(yīng)網(wǎng)格示意圖，由圖可以看出，本文提出的網(wǎng)格剖分法對劇烈的起伏海底同樣有很好的適應(yīng)性.圖10為波場快照對比，圖11為模擬地震記錄對比.由圖可見，本文方法相對與規(guī)則網(wǎng)格有限差分法不僅可以消除虛假繞射波，而且可以有效壓制數(shù)值頻散.

圖7 模擬地震記錄對比(a) 本文方法； (b) 規(guī)則網(wǎng)格有限差分法.Fig.7 The contract of seismic record(a) Our method； (b) Regular grid finite-difference method.

圖8 本文方法和規(guī)則網(wǎng)格有限差分法得到的反射波最大振幅誤差Fig.8 The errors of the maximum amplitude of the reflected waves by the proposed method and the MFD method

4.4 復(fù)雜海底速度模型

圖12為起伏海底多層模型，模型大小為6000 m×4500 m，速度分布如圖所示.計(jì)算域中的空間網(wǎng)格尺寸為10 m，時間采樣間隔為1 ms，最大記錄時間為6 s,震源是主頻為30 Hz的雷克子波，位于(3000 m,20 m)處，檢波點(diǎn)深度為10 m.圖13為利用本文方法得到的自適應(yīng)網(wǎng)格示意圖，由圖可見，利用本文提出的自適應(yīng)網(wǎng)格剖分法對復(fù)雜海底模型進(jìn)行網(wǎng)格剖分不僅可以更好地描述每一層界面，而且在每個界面附近網(wǎng)格明顯密集，從而提高了網(wǎng)格的分辨率.圖14為利用本文方法得到的模擬地震記錄和2 s時的波場快照.由圖可見，本文方法對復(fù)雜海底模型同樣有很好的適應(yīng)性.

圖9 起伏海底模型(a)和自適應(yīng)網(wǎng)格示意圖(b)Fig.9 Rough seabed model (a) and the sketch map of adaptive grid (b)

圖10 1.2 s的波場快照對比(a) 本文方法； (b) 規(guī)則網(wǎng)格有限差分法.Fig.10 The contract of wavefield snapshot at 1.2 s(a) Our method； (b) Regular grid finite-difference method.

圖11 模擬地震記錄對比(a) 本文方法； (b) 規(guī)則網(wǎng)格有限差分法.Fig.11 The contract of seismic record(a) The seismic record of our method； (b) Regular grid finite-difference method.

圖12 復(fù)雜海底模型Fig.12 The complex sea bottom model

圖13 復(fù)雜海底模型自適應(yīng)網(wǎng)格示意圖Fig.13 The sketch map of adaptive grid of complex seabed model

圖14 本文方法得到的復(fù)雜海底模型模擬地震記錄(a)和2 s時的波場快照(b)Fig.14 The seismic record (a) and snapshot at 2 s (b) of complex seabed velocity model

圖15 Marmousi速度模型Fig.15 Marmsousi velocity model

圖16 Marmousi模型自適應(yīng)網(wǎng)格示意圖Fig.16 The sketch map of adaptive grid of Marmousi

圖17 Marmousi模型單炮模擬地震記錄(a)和0.8 s時的波場快照(a)Fig.17 Single-shot seismic record (a) and snapshot at 2 s (b) of Marmousi model

4.5 Marmousi模型

圖15為Marmousi速度模型，模型大小為10000 m×3000 m.計(jì)算域中的空間網(wǎng)格尺寸為10 m，時間采樣間隔為1 ms，最大記錄時間為3 s，震源是主頻為30 Hz的雷克子波.圖16為利用本文方法得到的自適應(yīng)網(wǎng)格示意圖.由于數(shù)據(jù)量較大網(wǎng)格細(xì)節(jié)不能清晰展現(xiàn),但從整體上可以看出網(wǎng)格線在模型梯度大的區(qū)域顏色較深呈密集分布,而在梯度小的區(qū)域顏色較淺網(wǎng)格線相對稀疏.圖17a為利用本文方法得到的單炮模擬地震記錄，炮點(diǎn)位置為(3000 m,20 m),道間距為20 m,一共是250道.圖17b為0.8 s時刻的波場快照.由圖可以看出,利用本文方法對Marmousi模型進(jìn)行數(shù)值模擬不僅可以消除階梯狀網(wǎng)格引起的虛假繞射波，而且可以有效壓制數(shù)值頻散提高計(jì)算精度，從而驗(yàn)證了本文方法對復(fù)雜海底模型的強(qiáng)適應(yīng)性.

5 結(jié)論

為了有效消除規(guī)則網(wǎng)格有限差分法引起的虛假繞射波，本文將自適應(yīng)網(wǎng)格生成技術(shù)引入到地震波數(shù)值模擬的網(wǎng)格剖分中，采用仿真型有限差分法對曲線坐標(biāo)系下的聲波方程進(jìn)行了數(shù)值模擬.通過理論分析與模型試算可以得到以下結(jié)論：(1)利用自適應(yīng)網(wǎng)格對起伏海底進(jìn)行網(wǎng)格剖分相對與規(guī)則網(wǎng)格，不僅可以更準(zhǔn)確描述界面起伏，而且在界面附近可以加密網(wǎng)格.從而在不增加計(jì)算量的前提下消除階梯狀網(wǎng)格引起的虛假繞射波.(2)高階仿真型有限差分法相對于規(guī)則網(wǎng)格有限差分法可以有效降低數(shù)值頻散、提高計(jì)算精度.(3)相對于規(guī)則網(wǎng)格有限差分法，本文方法在相對較大的網(wǎng)格間距下同樣可以有效消除階梯網(wǎng)格引起的繞射波，從而可以有效減少內(nèi)存，提高計(jì)算效率.

附錄A 自適應(yīng)網(wǎng)格生成方法

定義表示映射后坐標(biāo)疏密性的泛函(Brackbill,1982)

(A1)

式中(x,z)為物理平面上的網(wǎng)格點(diǎn)，w(x,z)為加權(quán)函數(shù)，J為雅克比行列式，D為整個模型區(qū)域.

對式(A1)變量代換，該泛函變?yōu)?/p>

(A2)

泛函(A2)對應(yīng)變分問題的Euler方程為

(A3)

(A4)式中(x,z)是物理平面上的點(diǎn)，(ξ,η)是計(jì)算平面上的點(diǎn)

(A5)

定義一個表示映射后坐標(biāo)光滑性的泛函

(A6)

對式(A6)進(jìn)行變量代換，該泛函變成

(A7)

如果網(wǎng)格光滑性好，則式(A7)取極小值.故式(A7)的Euler方程為

(A8)

將式(A8)展開化簡得

bs1xξ ξ+bs2xξ η+bs3xη η+as1zξ ξ+as2zξ η+as3zη η=0,

as1xξ ξ+as2xξ η+as3xη η+cs1zξ ξ+cs2zξ η+cs3zη η=0,

(A9)

式中

as3=xξxη-zξzη, bs1=J+2xξzη,

cs3=J+2xξzη.

(A10)

在考慮網(wǎng)格的正交性，定義表示坐標(biāo)正交程度的泛函

(A11)

對式(A11)變量代換，該泛函變成

(A12)

式(A12)變分問題的Euler方程為

(A13)

將式(A13)展開并化成標(biāo)準(zhǔn)形式方程

bo1xξ ξ+bo2xξ η+bo3xη η+ao1zξ ξ+ao2zξ η+ao3zη η=0,

ao1xξ ξ+ao2xξ η+ao3xη η+co1zξ ξ+co2zξ η+co3zη η=0,

(A14)

式中

(A15)

通常，生成自適應(yīng)網(wǎng)格時需要同時考慮所有性質(zhì).設(shè)表示網(wǎng)格疏密性，光滑性和正交性的正的加權(quán)系數(shù)λv,λs,λo，做代數(shù)運(yùn)算λv×(A4)+λs×(A9)+λo×(A14)得

(A16)

式中

ai=λvavi+λsasi+λoaoi，

bi=λvbvi+λsbsi+λoboi，

ci=λvcvi+λscsi+λocoi， i=1,2,3.

附錄B 支撐算子法

在支撐算子法中散度和梯度滿足如下的積分關(guān)系：∫Vudiv w dV+∫V(w,grad u)dV=∮u(w,n)dS,

(B1)

式(B1)的左邊積分也可以寫成如下內(nèi)積形式：

(f,g)H=∫VfgdV,(a,b)H=∫V(a,b)dV,

(B2)

如果方程在邊界上的積分等于零則(B1)式變?yōu)?/p>

(u,div w)H+(grad u,w)H=0,

(B3)

因此微分算子div和grad互為負(fù)共軛：

grad=-div*.

(B4)

根據(jù)鏈?zhǔn)椒▌t向量A=(AX,AZ)在(ξ,η)坐標(biāo)系下的散度為

(B5)

式中

{(AX zη)ξ-(AX zξ)η}/J,

(B6)

把(B5)式帶到(B1)式的第一項(xiàng)中得

+[(AZ xξ)η-(AZ xη)ξ]u}dξdη,

(B7)

把(B6)式的第二項(xiàng)寫成如下形式：

∫V(A,grad u)dV=∫V(AX·GX+AZ·GZ)Jdξdη,

(B8)

利用分步積分法求(B7)并與(B8)比較得

(B9)

Alterman Z, Karal F C. 1968.Propagation of elastic waves in layered media by finite difference methods.BulletinoftheSeismologicalSocietyofAmerica, 58(1): 367-398.

Antonios E G, Antonie J E.1988. Directional control in grid generation.JournalofComputationalPhysics, 74(2): 422-439.Barfield W D. 1970. An optimal mesh generator for Lagrangian hydrodynamic calculations in two space dimensions.JournalofComputationalPhysics, 6(3): 417-429.

Bogey C, Bailly C. 2004. A family of low dispersive and low dissipative explicit schemes for flow and noise computations.JournalofComputationalPhysics, 194(1): 194-214.

Boore D M. 1972. Finite difference methods for seismic wave propagation in heterogeneous materials.MethodsinComputationalPhysics, 11: 1-37.

Brackbill J U, Saltzman J S. 1982. Adaptive zoning for singular problems in two dimensions.JournalofComputationalPhysics, 46(3): 342-368. Castillo J E, Hyman J M, Shashkov M J, et al. 2000. High-order mimetic finite difference methods on nonuniform grids.∥Proceedings of the Third International Conference on Spectral and High Order Methods. Houston, Texas: Special Issue of Houston Journal of Mathematics.

Crase E. 1990. High-order (space and time) finite-difference modeling of the elastic wave equation.∥ 60th SEG Annual Meeting.San Francisco, California: Society of Exploration Geophysicists, 987-991.

Jastram C, Behle A. 1992. Acoustic modelling on a grid of vertically varying spacing.GeophysicalProspecting, 40(2): 157-169.

Jastram C, Tessmer E. 1994. Elastic modelling on a grid with vertically varying spacing.GeophysicalProspecting, 42(4): 357-370.

Kang H W, Gu X Q, Liu C J. 2002. Study on weigh functions of adaptive grid technique.ActaMechanicaSinica(in Chinese), 34(5): 790-795. Kelly K R, Ward R W, Treitel S, et al. 1976. Synthetic seismograms: A finite-difference approach.Geophysics, 41(1): 2-27. Madariaga R. 1976. Dynamics of an expanding circular fault.BulletinoftheSeismologicalSocietyofAmerica, 66(3): 639-666. Moczo P. 1989.Finite-difference technique for SH-waves in 2-D media using irregular grids-application to the seismic response problem.GeophysicalJournalInternational, 99(2): 321-329.

Opr?al I, Zahradník J. 1999.Elastic finite-difference method for irregular grids.Geophysics, 64(1): 240-250.

Pitarka A. 1999. 3D elastic finite-difference modeling of seismic motion using staggered grids with nonuniform spacing.BulletinoftheSeismologicalSocietyofAmerica, 89(1): 54-68.

van Vossen R, Robertsson J O A, Chapman C H. 2002. Finite-difference modeling of wave propagation in a fluid-solid configuration.Geophysics, 67(2): 618-624.Virieux J.1984. SH-wave propagation in heterogeneous media: velocity-stress finite-difference method.Geophysics, 49(11): 1933-1942.

Yu M. 2004. Generation of 2-D adaptive structured grids by variational methods.ChineseJournalofComputationalPhysics(in Chinese), 21(1): 27-34.

Zhang J F. 1998. Non-orthogonal grid finite-difference method for numerical simulation of elastic wave propagation.ActaGeophysicaSinica(in Chinese), 41(S): 357-366. Zhu S W, Qu S L, Wei X C, etal. 2007. Numeric simulation by grid-various finite-difference elastic wave equation.OilGeophysicalProspecting(in Chinese), 42(6): 634-639.

附中文參考文獻(xiàn)

康紅文, 谷湘潛, 柳崇健. 2002. 自適應(yīng)網(wǎng)格技術(shù)中的權(quán)函數(shù)問題研究. 力學(xué)學(xué)報(bào), 34(5): 790-795.

于明. 2004. 二維自適應(yīng)結(jié)構(gòu)網(wǎng)格的變分生成方法. 計(jì)算物理, 21(1): 27-34.

張劍鋒. 1998. 彈性波數(shù)值模擬的非規(guī)則網(wǎng)格差分法. 地球物理學(xué)報(bào), 41(S): 357-366.

朱生旺, 曲壽利, 魏修成等. 2007.變網(wǎng)格有限差分彈性波方程數(shù)值模擬方法. 石油地球物理勘探, 42(6): 634-639.

(本文編輯 胡素芳)

Mimetic finite-difference numerical simulation of seismic wave based on the adaptive grid

LIU Zhi-Qiang, SUN Jian-Guo*, SUN Hui, LIU Ming-Chen, GAO Zheng-Hui, SHI Xiu-Lin

CollegeofGeoExplorationScienceandTechnology,JilinUniversity,Changchun130026,China

The drastic fluctuation of the earth′s surface and the sea bottom put forward a higher requirement to the numerical simulation of seismic waves. In the conventional finite difference methods, the rectangular mesh is used to divide the model. However, the rugged sea bottom and earth surface can only be represented by staircase lines, which will causes purious diffractions. Furthermore, when simulating the waves reflected at liquid-solid interfaces, denser meshes are needed to obtain sufficient accurate results for the interface does not coincide with the grid line.Thus, for the conventional finite difference method, denser meshes should be adopted in the whole study region, which will cause the increase of the computational cost and the memory.To solve the above problems, we first introduce the adaptive mesh generation technique to divide the velocity model with rugged sea bottom. Then, we employ high-order mimetic finite difference(MFD) method to simulate the acoustic equation under the curvilinear coordinate system. The adaptive mesh generation technique can not only represent the boundary of the model precisely but also suppress the spurious diffractions effectively. The high-order mimetic finite difference method can suppress the dispersion effectively and improve the computational precision. The numerical result show that the proposed method has a good adaptability to the complex sea bottom model.Keywords Spurious diffractions; Numerical dispersion; Adaptive grid; Mimetic finite-difference

10.6038/cjg20161225.

國家自然科學(xué)基金(41274120，41404085，41504084)資助.

劉志強(qiáng),男,1987年生,吉林大學(xué)地球探測與信息技術(shù)專業(yè)在讀博士，主要從事地震波數(shù)值模擬研究. E-mail: 490681597@qq.com

10.6038/cjg20161225

P631

2016-03-31，2016-09-20收修定稿

劉志強(qiáng)， 孫建國， 孫輝等. 2016. 基于自適應(yīng)網(wǎng)格的仿真型有限差分地震波數(shù)值模擬. 地球物理學(xué)報(bào)，59(12):4654-4665,

Liu Z Q, Sun J G, Sun H, et al. 2016. Mimetic finite-difference numerical simulation of seismic wave based on the adaptive grid.ChineseJ.Geophys. (in Chinese),59(12):4654-4665,doi:10.6038/cjg20161225.

*通訊作者 孫建國，男，1956年生，教授，主要從事波動理論與成像技術(shù)、地震資料處理方法與解釋技術(shù)等方面的教學(xué)和研究工作．

E-mail: sun_jg@jlu.edu.cn