交錯網(wǎng)格高階差分算法的改進(jìn)

白亞東，馮志民，楊慶節(jié)

（1.中國地質(zhì)大學(xué) 地球物理與信息技術(shù)學(xué)院，北京 100083；2.寧夏地球物理地球化學(xué)勘查院，銀川 750004；3.吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，長春 130026）

交錯網(wǎng)格高階差分算法的改進(jìn)

白亞東1，2，馮志民2，楊慶節(jié)3

（1.中國地質(zhì)大學(xué) 地球物理與信息技術(shù)學(xué)院，北京 100083；2.寧夏地球物理地球化學(xué)勘查院，銀川 750004；3.吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，長春 130026）

利用Taylor公式展開導(dǎo)出了交錯網(wǎng)格中可導(dǎo)函數(shù)任意次導(dǎo)數(shù)的任意偶數(shù)階精度的差分近似式及相應(yīng)的差分系數(shù)，從而改進(jìn)了常規(guī)高精度交錯網(wǎng)格有限差分算法中對各個空間導(dǎo)數(shù)采用不一致精度的問題。采用推導(dǎo)的交錯網(wǎng)格高階差分算法對一階彈性波動方程進(jìn)行數(shù)值模擬，得到了精確的數(shù)值模擬結(jié)果，證明了推導(dǎo)的交錯網(wǎng)格高階差分算法的正確性。

交錯網(wǎng)格；有限差分；數(shù)值模擬

0 引言

數(shù)值模擬問題出現(xiàn)在多個學(xué)科分支中，地震波場的數(shù)值模擬就是理解地震波在復(fù)雜介質(zhì)中傳播的不可或缺的工具。有限差分和偽譜法是兩種使用較為廣泛的地震波場數(shù)值模擬方法，雖然理論上講偽譜法是高階有限差分法中空間差分的階數(shù)達(dá)到無窮時的極限情況［1］，但是有限差分由于自身的易實(shí)現(xiàn)性，仍然得到了許多科研人員的重視。

有限差分算法不僅應(yīng)用于地震波場模擬上，地震波的偏移成像和反演等也多有應(yīng)用。為了提高有限差分算法的計算精度、穩(wěn)定性和適用性，該方法出現(xiàn)了多種變種方法：變網(wǎng)格有限差分［2－3］、不連續(xù)網(wǎng)格有限差分［4］、不規(guī)則網(wǎng)格有限差分［5］、交錯網(wǎng)格有限差分［6－7］、旋轉(zhuǎn)交錯網(wǎng)格有限差分［8－9］、可變時間步長有限差分［10］、自適應(yīng)可變空間步長網(wǎng)格有限差分以及隱式有限差分［11－12］。

交錯網(wǎng)格是由 Madariage［13］最早提出，Virieux［6］首先將其使用到一階速度－應(yīng)力方程中。Virieux發(fā)展的交錯網(wǎng)格的精度為O（Δt2＋Δx2），與常規(guī)網(wǎng)格相比既沒有增加計算量又沒有增加存儲空間，局部精度卻提高了4倍，收斂速度也有所提高。Levander［14］又將交錯網(wǎng)格有限差分的精度提高到O（Δt2＋Δx4）。由于交錯網(wǎng)格有限差分與傳統(tǒng)有限差分比較，具有更高的精度和更好的穩(wěn)定性，常被應(yīng)用于近地表復(fù)雜地形的模擬，粘彈性波場模擬［15－16］也已經(jīng)使用了交錯網(wǎng)格有限差分算法。

董良國等［17］發(fā)展了更高精度的交錯網(wǎng)格有限方法，精度達(dá)到了O（Δt4＋Δx2 N），為了使用較少的時間層，不增加計算存儲空間，將速度（應(yīng)力）對時間的三次導(dǎo)數(shù)準(zhǔn)確地轉(zhuǎn)化為應(yīng)力（速度）對空間的導(dǎo)數(shù)，把交錯網(wǎng)格和高階差分法成功地結(jié)合在一起，這種高階交錯網(wǎng)格有限差分解法，能夠適用于求解二維或者三維任意介質(zhì)情況下由速度－應(yīng)力表示的彈性波方程；李振春等［18］提出一種新的基于高階交錯網(wǎng)格技術(shù)，通過改變網(wǎng)格空間步長實(shí)現(xiàn)局部網(wǎng)格加密，彌補(bǔ)了常規(guī)網(wǎng)格中因采樣不足導(dǎo)致的頻散現(xiàn)象；李萌［19］在董良國的工作基礎(chǔ)上改進(jìn)了差分格式中三次導(dǎo)數(shù)的差分近似方法；劉洋等［20］通過引入Claerbout［21］的高階分母相，推導(dǎo)出一種新的隱式交錯網(wǎng)格有限差分公式，該公式對一階空間導(dǎo)數(shù)的（2 N＋2）階精度的差分，相當(dāng)于常規(guī)交錯網(wǎng)格有限差分公式的4 N階精度的差分，且在提高精度的同時，并不增加計算成本；劉洋等［11］還將Claerbout的思想與正規(guī)網(wǎng)格有限差分方法相結(jié)合，同樣得到了很好的效果。

作者著眼于常規(guī)交錯網(wǎng)格算法中對空間任意導(dǎo)數(shù)的高階差分，以Taylor公式展開為基礎(chǔ)，推導(dǎo)出了可導(dǎo)函數(shù)任意次導(dǎo)數(shù)的任意偶數(shù)階精度的差分近似式以及差分系數(shù)，從而改進(jìn)了O（Δt4＋Δx2 N）一階速度－應(yīng)力方程組差分格式和差分系數(shù).這里推導(dǎo)的差分格式與常規(guī)的高階交錯網(wǎng)格差分格式相似，但當(dāng)變量對空間求三次導(dǎo)數(shù)時，給出了更精確的差分系數(shù)，而差分系數(shù)的精確度正是決定交錯網(wǎng)格有限差分算法精度的關(guān)鍵。比較本方法和常規(guī)方法在相同條件和相同精度下的正演波形，結(jié)果顯示本方法的波形與解析解的誤差更小，穩(wěn)定性更好。

1 改進(jìn)的交錯網(wǎng)格差分格式

在二維TI介質(zhì)中，假定外力為零，一階速度－應(yīng)力彈性波方程為：

其中：vx、vz為分別為速度在x、z方向上的分量；τxx、τzz、τxz為應(yīng)力分量；ρ為密度；cij為介質(zhì)的彈性常數(shù)。式中的速度和應(yīng)力呈現(xiàn)一種耦合遞推的關(guān)系，適用于交錯網(wǎng)格差分算法。

1.1 時間上的差分近似

式中：Δt為時間步長，令M＝2，式（6）就是常規(guī)的時間4階精度差分近似。為了減少計算內(nèi)存，利用速度和應(yīng)力的耦合關(guān)系，得到式（1）―式（5）的時間4階精度差分近似，以式（1）為例：

1.2 空間上的差分近似

為了改進(jìn)精度為“O（Δt4＋Δx2 N）”的常規(guī)差分格式，這里對式（7）中各變量對空間二次、三次導(dǎo)數(shù)的差分格式進(jìn)行改進(jìn)，求解函數(shù)任意次（奇數(shù)次和偶數(shù)次）導(dǎo)數(shù)的任意偶數(shù)階精度的差分系數(shù)。

首先推導(dǎo)離散化函數(shù)值在半網(wǎng)格點(diǎn)上時，差分近似式的差分系數(shù)如圖1所示，f（x）的n次導(dǎo)數(shù)的差分使用的離散化函數(shù)值在半網(wǎng)格點(diǎn)上。

圖1 交錯的網(wǎng)格點(diǎn)Fig.1 Form of staggered points

式（8）寫成矩陣的形式有：

其中

簡單記為：X＝B·D，所以F＝D－1·B－1·S，若設(shè)F＝C·S，則有

式中C就是差分系數(shù)矩陣。由于D是初等矩陣，所以

則式（10）可寫為式（11）。

方便起見，這時只觀察式（11）中第一個式子，兩邊取轉(zhuǎn)置有，

因此有

即

通過求解式（12）可以得到一次導(dǎo)數(shù)差分格式的2 N階精度的差分系數(shù)，同理式（11）中其他的N－1個式子都可以得到類似式（12）的方程，即

式中C（N－i＋1）in表示函數(shù)f（x）的2i－1次導(dǎo)數(shù)的2（N －i＋1）階精度的差分系數(shù)（i，n＝1，2，…，N）。該式包含常規(guī)求一次導(dǎo)數(shù)差分近似式系數(shù)的矩陣方程，由于N不能小于i，所以此時的差分精度為2（N －i＋1）階，i＝1時也成立。因此有差分近似式：

然后推導(dǎo)離散化函數(shù)值在整網(wǎng)格點(diǎn)上時，差分近似式的差分系數(shù)如圖2所示，f（x）的n次導(dǎo)數(shù)的差分使用的離散化函數(shù)值在整網(wǎng)格點(diǎn)上.

令函數(shù)f（x）在x＝x0±nΔx處進(jìn)行Taylor展開，接下來的推導(dǎo)過程相同，不在詳述細(xì)節(jié)。

式（17）、式（18）分別是圖2所示情況下函數(shù)f（x）的任意奇數(shù)次和任意偶數(shù)次導(dǎo)數(shù)的差分近似式，其差分精度與交錯網(wǎng)格上的一致。同時給出了差分系數(shù)的求解矩陣方程。

1.3 改進(jìn)后的差分格式

在確定了差分系數(shù)之后，將對應(yīng)的變量對空間一次、二次和三次導(dǎo)數(shù)差分近似式代入到一階速度－應(yīng)力方程式（1）―式（5）中，經(jīng)過整理，得到其精度為O（Δt4＋Δx2 N）的改進(jìn)后準(zhǔn)確的差分格式，以式（1）為例：

圖2 非交錯的網(wǎng)格點(diǎn)Fig.2 Form of non－staggered points

表1 常用差分精度的穩(wěn)定性條件Tab.1 Sbility conditions of common difference accuracy

2 穩(wěn)定性條件

通過傅里葉分析方法，在二維TI介質(zhì)XOZ面內(nèi)，空間網(wǎng)格間距、時間采樣間隔和波速要滿足式（20）。

由式（20）計算得到幾種常用差分精度的穩(wěn)定性條件，如表1所示，其中。

3 精度對比以及模型試算

3.1 精度對比

為了驗證改進(jìn)后的交錯網(wǎng)格有限差分法的精度和正確性，將改進(jìn)前后的正演結(jié)果進(jìn)行比較。

各向同性均勻介質(zhì)如圖3所示，物性參數(shù)縱波速度為3 000m／s，橫波速度為2 000m／s，介質(zhì)密度為2 000kg／m3。震源坐標(biāo)為（1 000m，1 000 m），采用主頻為25Hz的Richer子波激發(fā)。時間采樣間隔為0.001s，模型網(wǎng)格大小為200m×200 m，水平和垂直空間步長都為10m，模型自由表面采用PML吸收邊界。

圖3 各向同性均勻介質(zhì)模型Fig.3 Isotropic medium model

圖4 250ms時模型波場快照Fig.4 Snapshots of Fig3at 250ms

分別使用不同精度的常規(guī)交錯網(wǎng)格有限差分公式和本文推導(dǎo)的交錯網(wǎng)格有限差分公式對一階速度－應(yīng)力波動方程（方程（1）－方程（5））進(jìn)行模擬，模擬結(jié)果圖4中展示。

圖4分別是精度為O（Δt4＋Δx4）和O（Δt4＋Δx10）常規(guī)交錯網(wǎng)格有限差分算法改進(jìn)前后的250 ms時的波場快照，可以看出，改進(jìn)前、后的高精度交錯網(wǎng)格有限差分算法在相同精度下，改進(jìn)后的模擬結(jié)果的精度有較大的提高，波場頻散情況明顯改善。為了更加清晰地對比改進(jìn)前后的精度，將其精度為O（Δt4＋Δx10）正演波形同時與解析解波形進(jìn)行比較。圖5給出的是模型中坐標(biāo)為（500m，500m）質(zhì)點(diǎn)處的振動記錄數(shù)值解與解析解的比較，由圖5可見，改進(jìn)后的數(shù)值解波形與解析解更加一致，說明改進(jìn)后的高階交錯網(wǎng)格有限差分算法精度更高。

圖5 算法改進(jìn)前后的質(zhì)點(diǎn)振動記錄與解析解的對比Fig.5 Comparison of particle vibration recorded before and after the improved algorithm and analytic solutions

3.2 雙層模型試算

為了研究波場在層狀介質(zhì)中的傳播規(guī)律，建立層狀地質(zhì)模型，圖6是分層均勻介質(zhì)模型，每層的物性參數(shù)如圖6所示。網(wǎng)格大小為250m×250m，空間步長為10m，兩個介質(zhì)分界面分別在500m深和1 000m深處。震源坐標(biāo)為（200m，1 250m），計算時采用改進(jìn)后的時間4階空間10階精度的差分格式，用主頻為25Hz的Richer子波激發(fā)。模型四周均采用PML吸收邊界。

圖6 分層介質(zhì)模型Fig.6 Layered medium model

圖7分別是精度為O（Δt4＋Δx8）的改進(jìn)后算法的400ms時的波場快照，可以看出，完善后的高精度交錯網(wǎng)格有限差分算法，能精確地模擬彈性波在地下層狀介質(zhì)中的傳播特性，各種反射和透射波都可以觀察到。在自由表面也沒有看到明顯的反射，說明將PML吸收邊界應(yīng)用到改進(jìn)后的交錯網(wǎng)格算法中，對彈性波也有很好的吸收效果。

圖7 400ms時模型波場快照Fig.7 Snapshots of Fig6at 400ms

4 結(jié)論

交錯網(wǎng)格有限差分以其高效、精確、實(shí)用等優(yōu)點(diǎn)常被用于波動方程模擬。作者通過Taylor公式展開法推導(dǎo)了函數(shù)任意次導(dǎo)數(shù)的任意階精度的差分格式及其差分系數(shù)，改進(jìn)了傳統(tǒng)的高精度交錯網(wǎng)格有限差分公式。通過精度對比及模型試算的結(jié)果表明，在計算量增加不大情況下，改進(jìn)后的高精度交錯網(wǎng)格有限差分公式比常規(guī)傳統(tǒng)的公式精度更高。

［1］ 唐小平，白超英，劉寬厚.偽譜法分離波動方程彈性波模擬［J］.石油地球物理勘探，2012，47（1）：19－26.

TANG X P，BAI CH Y，LIU K H.Elastic wavefield simulation using separated equations through pseudospectral method［J］.OGP，2012，47（1）：19－26.（In Chinese）

［2］ WANG Y，XU J，GERARD T S.Viscoelastic wave simulation in basins by a variable－grid finite－difference method［J］.Bull Seismol Soc Am，2001，91（6）：1741 －1749.

［3］ NARAYANJ P，KUMAR S.A fourth order accurate SH－wave staggered grid finite－difference algorithm with variable grid size and VGR－stress imaging technique［J］.Pure applGeophys，2008，165：271－294.

［4］ SHIN A，HIROYUKI F.3Dfinite－difference method using discontinuous grids［J］.Bull SeismolSoc Am，1999，89：918－930.

［5］ OPR?AL I，ZAHRADNíK J.Elastic finite－difference method for irregular grids［J］.Geophysics，1999，64：240－250.

［6］ VIRIEUX J.SH－wave propagation in heterogeneous media：velocity－stress finite－difference method：Geophysics［J］.1984，49（11），1933－1957.

［7］ THOMAS B，ERIK H S.Accuracy of heterogeneous staggered－grid finite－difference modeling of rayleigh waves［J］.Geophysics，2006，71（4）：T109－T115.

［8］ ERIK H S，NORBERT G，SERGE A S.Modeling the propagation of elastic waves using a modified finitedifference grid［J］.Wave Motion，2000，31：77－92.

［9］ REESHIDEV B，MRINAL K S.Finite－difference modelling of S－wave splitting in anisotropic media［J］.Geophysical Prospecting，2008，56：293－312.

［10］TESSMER E.Seismic finite－difference modeling with spatially varying time steps［J］.Geophysics，2000，65 （4）：1290－1293.

［11］LIU Y，MRINAL K S.A practical implicit finitedifference method：examples from seismic modeling ［J］.Journal of Geophysics and Engineering，2009（6）：231－249.

［12］LIU Y，MRINAL K S.Finite－difference modeling with adaptive variable－length spatial operators［J］.Geophysics，2011，76（4）：T79－T89.

［13］MADARIAGA R.Dynamics of an expanding circular fault［J］.Bull Seismol Soc Am，1976，66，639－666.

［14］LEVANDER A R.Fourth－order finite－difference PSV seismograms［J］.Geophysics，1988，53（11）：1425 －1436.

［15］OPERTO S，VIRIEUX J，AMESTOY P，et al.3Dfinite－difference frequency－domain modeling of visco－acoustic wave propagation using a massively parallel direct solver：Afeasibility study［J］.Geophysics，2007，72：SM195－SM211.

［16］LIAO J P，WANG H Z，AND MA Z T.2－D elastic wave modeling with frequency－space 25－point finitedifference operators［J］.Applied Geophysics，2009，6 （3）：259－266.

［17］董良國，馬在田，曹景中.一階彈性波方程交錯網(wǎng)格高階差分解法穩(wěn)定性研究［J］.地球物理學(xué)報，2000，43 （3）：411－419.

DONG L G，MA Z T，CAO J ZH.The Staggered－Grid High－Order Difference Method of One－Order E－lastic Equation［J］.Chinese Journalof Geophysics，2000，43（3）：411－419.（In Chinese）

［18］李振春，張慧，張華.一階彈性波方程的變網(wǎng)格高階有限差分?jǐn)?shù)值模擬［J］.石油地球物理勘探，2008，43（6）：711－716.

LI ZH CH，ZHANG H，ZHANG H.Variable－grid high－order finite－difference numeric simulation of first －order elastic wave equation［J］.OGP，2008，43（6）：711－716.（In Chinese）

［19］李萌.改進(jìn)的一階彈性波方程交錯網(wǎng)格高階差分解法［J］.地球物理學(xué)進(jìn)展，2009，24（3）：1065－1068.

LI M.Improvement to the staggered－grid high－order difference method of one－order elastic wave equation ［J］.Progress in Geophys，2009，24（3）：1065－1068.（In Chinese）

［20］LIU Y，MRINAL K S.An implicit staggered－grid finite－difference method for seismic modeling［J］.Geophysical Journal International，2009，179：459－474.

［21］CLAERBOUT J F.The craft of wavefield extrapolation，in Imaging theEarth's Interior［M］.Oxford：Blackwell Scientific Publications，1985.

Improvement of staggered grid finite difference algorithm

BAI Ya－dong1，2，F(xiàn)ENG Zhi－min2，YANG Qing－jie3
（1.School of Geophysics and Information Technology，China University of Geosciences，Beijing 100083，China；2.Ningxia Academy of Geophysical and Geochemical Exploration，Ningxia 750004，China；3.College of GeoExploration Science and Technology，Jilin University，Changchun 130026，China）

A finite difference scheme with any－order accuracy and corresponding coefficients of any number of derivatives of functions which have any number of derivatives by Taylor formula is deduced in this paper.It is helpful to improve the staggered grid finite difference algorithm which has inconsistent accuracy in difference space derivatives.We make a simulation of seismic response with conventional FD scheme and the new FD scheme respectively.As a result，the new FD scheme is more stabilized and more accurate.

staggered grid；finite difference；numerical simulation

P 631.4

：A

10.3969／j.issn.1001－1749.2015.06.07

1001－1749（2015）06－0709－07

2014－11－26改回日期：2015－03－24

寧夏回族自治區(qū)地勘基金項目（2014－地勘基金07）

白亞東（1975－），男，高級工程師，長期從事地球物理勘探工作，E－mail：byd1975＠163.com。