一種基于全局優(yōu)化的交錯(cuò)網(wǎng)格有限差分法

印興耀 劉 博 楊鳳英

（中國(guó)山東青島266580中國(guó)石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院）

引言

地震波數(shù)值模擬是在地震波傳播理論的基礎(chǔ)上，通過數(shù)值計(jì)算來模擬地震波在地下介質(zhì)中的傳播（董良國(guó)，2003），是研究地震波傳播特性與地球介質(zhì)參數(shù)關(guān)系的重要手段.交錯(cuò)網(wǎng)格有限差分技術(shù)已經(jīng)廣泛應(yīng)用于各向異性介質(zhì)（Igel et al，1995；Collino，Tsogka，2001；裴正林，王尚旭，2005；殷文等，2006；何燕，2008）、孔隙介質(zhì)（Dai et al，1995；Zeng et al，2001；孫林潔，2012）、黏彈性介質(zhì)（白曉寅，2008；孫成禹等，2010）的波動(dòng)方程模擬中，但是常規(guī)交錯(cuò)網(wǎng)格在模擬非均勻性較強(qiáng)的介質(zhì)和各向異性介質(zhì)時(shí)需對(duì)介質(zhì)參數(shù)進(jìn)行平均或內(nèi)插（陳浩等，2006），低階插值可能會(huì)導(dǎo)致精度降低，高階插值則會(huì)導(dǎo)致計(jì)算量增大.Gold等（1997）提出了旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分算法，Saenger等進(jìn)一步研究了旋轉(zhuǎn)交錯(cuò)網(wǎng)格的穩(wěn)定性與頻散的關(guān)系（Saenger et al，2000；Saenger，Bohlen，2004）.在旋轉(zhuǎn)交錯(cuò)網(wǎng)格中采用沿網(wǎng)格對(duì)角線差分的方式，可以避免彈性模量的平均與插值，因而更準(zhǔn)確.隨后一些研究人員通過旋轉(zhuǎn)交錯(cuò)網(wǎng)格的完全匹配層（perfectly matched layer，簡(jiǎn)寫為PML）邊界條件的實(shí)現(xiàn)，利用旋轉(zhuǎn)交錯(cuò)網(wǎng)格技術(shù)對(duì)各向異性介質(zhì)（李敏，劉洋，2012）和黏彈介質(zhì)（嚴(yán)紅勇，劉洋，2012）等進(jìn)行了數(shù)值模擬.但由于旋轉(zhuǎn)交錯(cuò)網(wǎng)格中差分是沿對(duì)角線進(jìn)行，步長(zhǎng)是同等條件下常規(guī)交錯(cuò)網(wǎng)格的■2倍（網(wǎng)格是正方形的情況下），若用常規(guī)的差分格式會(huì)導(dǎo)致離散的梯度算子和散度算子在計(jì)算時(shí)誤差增大，因而容易出現(xiàn)數(shù)值頻散.而如果使用高階有限差分則會(huì)因需要更多的網(wǎng)格點(diǎn)而增大計(jì)算量與內(nèi)存.緊致有限差分是一種隱式差分格式，在同等網(wǎng)格數(shù)的條件下，具有比常規(guī)差分格式更高的精度以及更低的數(shù)值頻散，故成為目前關(guān)注的熱點(diǎn).Lele（1992）分析了緊致差分格式的分辨率特性；王書強(qiáng)等（2002）對(duì)彈性波方程的緊致差分格式及其與中心差分的誤差進(jìn)行了研究；Du等（2009）基于交錯(cuò)網(wǎng)格緊致有限差分進(jìn)行了橫向各向同性介質(zhì)的正演模擬；Boersma（2011）利用6階緊致有限差分求解了Navier-Stokes方程.

實(shí)現(xiàn)有限差分時(shí)，差分系數(shù)的確定會(huì)影響數(shù)值模擬的精度.有限差分系數(shù)的求取一般包括泰勒展開和最優(yōu)化算法兩種方法（Liu，2013）.基于泰勒展開得到的差分系數(shù)在波數(shù)較小時(shí)精度高，但隨波數(shù)增大其精度會(huì)降低；基于最優(yōu)化算法的差分系數(shù)求取方法是在給定精度與差分算子長(zhǎng)度的情況下盡可能拓寬波數(shù)范圍（Holberg，1987）.Shan（2009），Kosloff等（2010），Zhou和Zhang（2011）以及Liu（2013）等對(duì)基于最小二乘算法的差分算子優(yōu)化方法進(jìn)行了深入研究；Zhang和Yao（2012）則基于模擬退火算法對(duì)中心差分算子進(jìn)行了優(yōu)化.

本文擬結(jié)合旋轉(zhuǎn)交錯(cuò)網(wǎng)格與緊致有限差分技術(shù)，基于模擬退火全局優(yōu)化算法對(duì)緊致差分算子進(jìn)行優(yōu)化，拓寬數(shù)值模擬的波數(shù)范圍，在此基礎(chǔ)上進(jìn)行頻散分析，最后通過數(shù)值模擬驗(yàn)證該方法的可行性和有效性.

1 旋轉(zhuǎn)交錯(cuò)網(wǎng)格緊致有限差分方法

1.1 一階速度－應(yīng)力方程

在彈性波正演模擬中，除了使用二階方程外，還常常采用一階速度－應(yīng)力彈性波方程，其主要優(yōu)點(diǎn)是無需對(duì)彈性常數(shù)進(jìn)行空間差分（Virieux，1984）.這里根據(jù)運(yùn)動(dòng)平衡方程和本構(gòu)關(guān)系給出二維情況下體力為零時(shí)各向異性介質(zhì)的一階速度－應(yīng)力方程：

式中：v，σ分別代表速度和應(yīng)力，1表示x方向，3表示z方向；cij為介質(zhì)彈性張量矩陣的元素，當(dāng)極端各向異性介質(zhì)時(shí)，該矩陣有21個(gè)獨(dú)立參數(shù)；當(dāng)二維時(shí)，不考慮y分量，故張量矩陣中與y分量有關(guān)參數(shù)為0，只有5個(gè)獨(dú)立參數(shù)（c11，c13＝c31，c15＝c51，c33，c35＝c53，c55）.

1.2 緊致有限差分

在地震波正演模擬中，采用常規(guī)的中心有限差分往往需要較多的網(wǎng)格點(diǎn)才能比較有效地壓制數(shù)值頻散，不利于邊界的處理（王書強(qiáng)等，2002）.緊致有限差分是一種隱式差分格式，只需要較少的網(wǎng)格點(diǎn)即可有效地壓制數(shù)值頻散，且同等階次的精度要比常規(guī)中心有限差分格式高，因此本文采用緊致有限差分格式，彌補(bǔ)旋轉(zhuǎn)交錯(cuò)網(wǎng)格因差分步長(zhǎng)大導(dǎo)致梯度和散度算子在計(jì)算時(shí)誤差增大因而更容易產(chǎn)生數(shù)值頻散的不足.

常規(guī)的2（M－1）點(diǎn)2 M階緊致有限差分格式（Kim，Lee，1996）為

式中，Δx為網(wǎng)格的空間步長(zhǎng)，a和bn分別為差分的權(quán)系數(shù)，可以通過在第i個(gè)網(wǎng)格點(diǎn)的泰勒展開，然后對(duì)比對(duì)應(yīng)系數(shù)來求解

即可得到2（M－1）點(diǎn)2 M階緊致有限差分格式的差分系數(shù)（Liu，Sen，2009），例如當(dāng)M＝5時(shí)，可以得到8點(diǎn)10階的緊致差分格式系數(shù)為a＝0.257 894 74，b1＝0.889 871 16，b2＝0.216 121 16，b3＝－0.004 701 2，b4＝0.000 151 55.

1.3 旋轉(zhuǎn)交錯(cuò)網(wǎng)格及完全匹配層吸收邊界條件

在二維旋轉(zhuǎn)交錯(cuò)網(wǎng)格上，密度和速度的各個(gè)分量定義在相同位置，彈性模量和應(yīng)力的各個(gè)分量定義在另外的相同位置上，如圖1所示.旋轉(zhuǎn)交錯(cuò)網(wǎng)格在計(jì)算時(shí)分為兩步（陳浩等，2006）：首先，計(jì)算沿著對(duì)角線方向即?x，?z方向進(jìn)行差分，得到相關(guān)物理量的對(duì)角線方向的空間一階導(dǎo)數(shù)；然后，通過坐標(biāo)旋轉(zhuǎn)的換算關(guān)系，將得到的兩個(gè)對(duì)角線方向的差分進(jìn)行線性組合從而獲得坐標(biāo)軸方向即x，z方向的一階空間導(dǎo)數(shù)，其換算關(guān)系為（Saenger et al，2000）：

圖1 旋轉(zhuǎn)交錯(cuò)網(wǎng)格及完全匹配層吸收邊界示意圖x，z為坐標(biāo)軸方向，x?，z?為對(duì)角線方向；灰色區(qū)域表示沿著x和z方向均進(jìn)行衰減，白色區(qū)域表示只沿著x方向或者z方向進(jìn)行衰減Fig.1 Schematic diagram of rotated staggeredgrid and perfectly matched layer（PML）absorbing boundary xand zare the coordinate directions，andx?andz? are the diagonal directions.Gray areas indicate that waves are absorbed along both xand z directions and white areas indicate that waves are absorbed only along the xor zdirection

利用波動(dòng)方程進(jìn)行數(shù)值模擬，一個(gè)關(guān)鍵問題就是邊界條件.為了解決邊界反射問題，引入彈性波方程的PML邊界條件.陳浩等（2006）指出，在旋轉(zhuǎn)交錯(cuò)網(wǎng)格中，PML吸收邊界條件的處理方式及吸收效果與常規(guī)交錯(cuò)網(wǎng)格中幾乎是一樣的.為此，引入圖1所示的PML吸收層.PML邊界的基本做法是在研究區(qū)域四周引入PML，波在PML中傳播時(shí)不會(huì)產(chǎn)生反射，并且隨傳播距離按一定規(guī)律衰減.當(dāng)波傳播到PML邊界時(shí)，波場(chǎng)近似為零，也不會(huì)產(chǎn)生反射（王守東，2003）.

以式（1）中的vx為例來說明旋轉(zhuǎn)交錯(cuò)網(wǎng)格緊致有限差分算法PML的實(shí)現(xiàn)方法：

求出.其中

類似地可以推導(dǎo)出其它分量的旋轉(zhuǎn)交錯(cuò)網(wǎng)格緊致有限差分下的PML控制方程.

1.4 緊致有限差分的全局優(yōu)化

根據(jù)平面波理論，令

式中，β＝kΔx／2，0≤β≤π／2.

為了使式（2）左右兩邊的誤差盡可能小，須滿足β－β＊≈0的β取值范圍應(yīng)盡可能地大.

本文中，采用模擬退火全局最優(yōu)算法（Kirkpatrick et al，1983）來優(yōu)化旋轉(zhuǎn)交錯(cuò)網(wǎng)格緊致有限差分算子.為此，建立以下目標(biāo)函數(shù)：

Kim和Lee（1996）指出，當(dāng)θ接近最大值π／2時(shí)，會(huì)出現(xiàn)很多難以控制的誤差，從而導(dǎo)致優(yōu)化效果不佳，因此本文選擇θ＝rπ／2，其中0＜r＜1.表1列出了優(yōu)化的8點(diǎn)10階，10點(diǎn)12階，12點(diǎn)14階以及14點(diǎn)16階的緊致有限差分系數(shù).

表1 優(yōu)化的10—16階緊致有限差分系數(shù)Table 1 Optimized coefficients of 10—16-order compact finite-difference

2 頻散分析

由式（12）可以看出，β接近于β＊的程度表征了優(yōu)化的緊致有限差分算子的精度，故定義

使α≈1的β取值范圍越大，數(shù)值頻散越小.

圖2給出了常規(guī)的基于泰勒展開得到的差分算子與本文提出的全局優(yōu)化算子的精度對(duì)比.可以看出，對(duì)于常規(guī)的基于泰勒展開得到的緊致有限差分算子，隨著空間階次的增大，使α≈1的波數(shù)范圍變得越寬，即數(shù)值模擬的精度隨階次的增大而增大.顯然，本文得到的全局優(yōu)化的緊致有限差分算子在相同最大誤差下使α≈1的波數(shù)范圍比常規(guī)緊致有限差分算子還要大.從圖2中還可以看出，優(yōu)化的8點(diǎn)10階緊致有限差分算子的頻散比常規(guī)14點(diǎn)16階緊致有限差分算子小，即優(yōu)化的10階精度要高于常規(guī)的16階精度，從而可以使用較少的網(wǎng)格點(diǎn)來達(dá)到較高的精度，節(jié)省計(jì)算內(nèi)存.

圖2 10階常規(guī)有限差分，10—16階常規(guī)緊致有限差分及對(duì)應(yīng)的優(yōu)化緊致有限差分的頻散對(duì)比 （k代表波數(shù)，Δx表示空間網(wǎng)格間距）Fig.2 Dispersion comparison by 10-order conventional finite-difference，10—16-order conventional compact finite-difference and the corresponding optimized compact finitedifference.k denotes the wavenumber，andΔxis spatial grid spacing

由圖2中的10階常規(guī)有限差分、10階常規(guī)緊致有限差分與10階優(yōu)化緊致有限差分的頻散曲線可以看出，10階常規(guī)緊致有限差分格式的波數(shù)范圍較10階常規(guī)有限差分要寬，而經(jīng)過優(yōu)化的10階緊致有限差分格式的波數(shù)范圍則更進(jìn)一步拓寬.設(shè)最大允許誤差為εmax＝｜α－1｜0≤β≤βmax，取εmax＝0.2%，分別計(jì)算得到10階常規(guī)有限差分、10階常規(guī)緊致有限差分及10階優(yōu)化緊致有限差分的最大波數(shù)范圍βmax分別為0.847 5，1.037 5，1.458 0，即10階常規(guī)緊致有限差分的波數(shù)范圍比常規(guī)差分格式拓寬了1.22倍，而經(jīng)過全局優(yōu)化后，波數(shù)范圍提升至1.72倍.考慮到正演模擬中最常采用Δx＝Δz的情況，則對(duì)于旋轉(zhuǎn)交錯(cuò)網(wǎng)格，Δr比常規(guī)交錯(cuò)網(wǎng)格差分步長(zhǎng)Δx及Δz增大了2倍，因而可以采用優(yōu)化系數(shù)的緊致有限差分來彌補(bǔ)這一不足，從而有效地壓制數(shù)值頻散.需要指出的是，優(yōu)化系數(shù)的緊致有限差分格式也可以很方便地應(yīng)用于PML邊界條件，而不需要特殊的處理.

3 模型試算

3.1 模型1

為了驗(yàn)證全局優(yōu)化的旋轉(zhuǎn)交錯(cuò)網(wǎng)格緊致有限差分的精度比常規(guī)旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分方法有所提升，首先采用一個(gè)簡(jiǎn)單的各向同性模型進(jìn)行測(cè)試（模型大小為3 000m×3 000m，在測(cè)試中并未添加任何邊界條件），震源采用主頻30Hz的雷克子波，z方向集中力源激發(fā)，網(wǎng)格大小為Δx＝Δz＝10m，時(shí)間采樣為Δt＝1ms.縱橫波速度分別為4 000m／s和2 600m／s，模型密度為2 300kg／m3.圖3給出了3種有限差分條件下400ms時(shí)刻x分量波場(chǎng)快照，均采用旋轉(zhuǎn)交錯(cuò)網(wǎng)格系統(tǒng).圖3a采用常規(guī)10階常規(guī)有限差分格式；圖3b采用8點(diǎn)10階，10點(diǎn)12階，12點(diǎn)14階，14點(diǎn)16階常規(guī)緊致有限差分格式；圖3c采用優(yōu)化的8點(diǎn)10階緊致有限差分格式.可以看出，10階常規(guī)有限差分與10階緊致有限差分對(duì)比，后者頻散要?。欢鴮?duì)于圖3b中的緊致有限差分格式，隨著空間階次的增大，其精度越高，頻散越??；對(duì)比圖3b與圖3c可以看出，優(yōu)化的8點(diǎn)10階緊致有限差分的數(shù)值頻散甚至比常規(guī)14點(diǎn)16階緊致有限差分還要小，這與圖2所示是一致的.圖4給出了從圖3中提取的z＝500m處（即50個(gè)網(wǎng)格點(diǎn)處）的波形曲線，以此來對(duì)比以上3種差分格式的數(shù)值頻散.可以看出，10階常規(guī)有限差分算法在波形外出現(xiàn)了數(shù)值抖動(dòng)，說明數(shù)值頻散比較高，模擬精度相對(duì)較低，而10階緊致有限差分和優(yōu)化的10階緊致差分精度較高，從虛線框中可以進(jìn)一步看出優(yōu)化的10階緊致有限差分算法的精度更高一些.

圖3 10階常規(guī)有限差分（a），10—16階常規(guī)緊致有限差分（b）與優(yōu)化的10階緊致有限差分（c）得到的波場(chǎng)快照Fig.3 Snapshots of wave field by using 10-order conventional finite-difference（a），10－16-order conventional compact finite-difference（b）and optimized 10-order compact finite-difference（c）

3.2 模型2

圖4 10階常規(guī)有限差分，10階常規(guī)緊致有限差分及10階優(yōu)化緊致有限差分波形對(duì)比圖Fig.4 Comparison of waveforms by using 10-order conventional finite-difference（FD）（dark line），10-order conventional compact FD （blue line）and optimized 10-order compact FD （red line）

模型2的幾何構(gòu)造如圖5所示，其中凹陷構(gòu)造上部及下部均為各向同性介質(zhì)，最下層為橫向各向同性（VTI）介質(zhì).模型縱向和橫向長(zhǎng)為3 000m，具體參數(shù)見表2.采用優(yōu)化的旋轉(zhuǎn)交錯(cuò)網(wǎng)格緊致有限差分方法進(jìn)行正演模擬，網(wǎng)格大小為Δx＝Δz＝10m，時(shí)間采樣 為Δt＝1ms.震源采用30Hz雷克子波，在模型的（1 150m，1 500m）處以縱波源形式激發(fā).模擬過程中，采用PML邊界條件（PML層部分同樣采用優(yōu)化的旋轉(zhuǎn)交錯(cuò)網(wǎng)格緊致有限差分算法），厚度層為30個(gè)網(wǎng)格點(diǎn).為體現(xiàn)吸收效果包含了PML層，圖6給出了600ms時(shí)刻x及z分量的波場(chǎng)快照.可以看出，震源激發(fā)的縱波傳播至分界面時(shí)，發(fā)生了反射與透射，可以見到PP波、PS波等，并可在凹陷處看到斷面波，波場(chǎng)復(fù)雜，但波前面清晰，數(shù)值頻散小.從圖7所示的單炮記錄中也可以看出優(yōu)化的旋轉(zhuǎn)交錯(cuò)網(wǎng)格緊致有限差分算法具有較高的模擬精度，各類波均可得到清晰體現(xiàn)，并且能夠證明所加PML邊界條件的吸收效果很好，無明顯人為邊界反射產(chǎn)生.

圖5 凹陷模型示意圖Fig.5 Schematic diagram of the depression model

圖6 600ms時(shí)刻優(yōu)化的旋轉(zhuǎn)交錯(cuò)網(wǎng)格緊致有限差分正演x分量（左）和z分量（右）的波場(chǎng)快照（虛線框外圍代表PML層）Fig.6 Wavefield snapshots of x-（left）and z-component（right）at 600ms using optimized rotated staggered-grid compact finite-difference method（PML layers are outside the dashed box）

表2 凹陷模型參數(shù)Table 2 Parameters of the depression model

圖7 凹陷模型優(yōu)化的旋轉(zhuǎn)交錯(cuò)網(wǎng)格緊致有限差分正演x分量（左）和z分量（右）單炮記錄Fig.7 Records of x-（left）and z-component（right）of the depression model using optimized rotated staggered-grid compact finite-difference method

4 討論與結(jié)論

旋轉(zhuǎn)交錯(cuò)網(wǎng)格在地震波正演模擬中由于其可以避免在模擬中因彈性模量插值或平均而產(chǎn)生的誤差和其對(duì)邊界條件處理的便捷性而被廣泛應(yīng)用.但由于其差分方向是沿網(wǎng)格的對(duì)角線方向，因此相對(duì)于常規(guī)的交錯(cuò)網(wǎng)格而言，其差分步長(zhǎng)較大，梯度算子、散度算子在計(jì)算時(shí)更容易產(chǎn)生誤差，這意味著旋轉(zhuǎn)交錯(cuò)網(wǎng)格更容易產(chǎn)生數(shù)值頻散，從而影響模擬的精度.若采用小網(wǎng)格，雖然可以壓制數(shù)值頻散，卻又大大增加了運(yùn)算量及內(nèi)存.因此本文采用旋轉(zhuǎn)交錯(cuò)網(wǎng)格與緊致有限差分技術(shù)相結(jié)合的方法來提高模擬精度，并在此基礎(chǔ)上，基于模擬退火算法對(duì)緊致有限差分進(jìn)行全局優(yōu)化，進(jìn)一步提高計(jì)算精度，壓制數(shù)值頻散.

數(shù)值頻散分析結(jié)果表明，優(yōu)化的旋轉(zhuǎn)交錯(cuò)網(wǎng)格緊致有限差分算法相對(duì)于普通旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分算法具有更寬的波數(shù)范圍，這意味著優(yōu)化的旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分算法可以在采用更大的空間采樣間隔或者更高的震源主頻時(shí)依然有較高的精度和較低的數(shù)值頻散；且經(jīng)過全局優(yōu)化的10階緊致有限差分算子比常規(guī)的16階緊致有限差分算子具有更高的精度，即可以采用較少的網(wǎng)格點(diǎn)來達(dá)到較高的模擬精度，因此可以節(jié)省計(jì)算內(nèi)存.

數(shù)值模擬實(shí)驗(yàn)進(jìn)一步證實(shí)了本文提出的優(yōu)化的旋轉(zhuǎn)交錯(cuò)網(wǎng)格緊致有限差分算法的正確性與可行性.凹陷模型中，各類波均具有比較清晰的響應(yīng)，揭示了波場(chǎng)的傳播規(guī)律.同時(shí)波場(chǎng)快照與單炮記錄中無明顯的人為邊界反射產(chǎn)生，也證明了模擬采用的PML吸收邊界條件的有效性.

但是本文提出的方法亦存在不足之處.在求解緊致差分格式時(shí)需要解三角矩陣線性方程，較常規(guī)顯式有限差分格式增加額外的計(jì)算負(fù)擔(dān)，且不利于在GPU等平臺(tái)直接并行，例如圖3a和圖3c中的400ms時(shí)長(zhǎng)的計(jì)算時(shí)間分別為46.5s和85.6s，后者計(jì)算耗時(shí)約為前者的1.8倍.在科學(xué)計(jì)算中，效率與計(jì)算精度常常不可兼得，隨著計(jì)算機(jī)科學(xué)計(jì)算能力的提高，可以采用GPU加速的LU分解算法或者使用LAPACK等優(yōu)化程序包來克服求解大型矩陣的計(jì)算效率問題.在實(shí)際應(yīng)用中需權(quán)衡利弊，選擇適合的算法.

衷心感謝審稿專家提出的寶貴意見和建議.

白曉寅.2008.基于地震波衰減理論的地層吸收參數(shù)提取方法研究［D］.東營(yíng)：中國(guó)石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院：24-54.

Bai X Y.2008.The Study on Method of Layer Absorbing Parameters Extraction Based on the Theory of Seismic Wave Attenuation［D］.Dongying：School of Geosciences，China University of Petroleum：24-54（in Chinese）.

陳浩，王秀明，趙海波.2006.旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分及其完全匹配層吸收邊界條件［J］.科學(xué)通報(bào)，51（17）：1985-1994.

Chen H，Wang X M，Zhao H B.2006.Rotated staggered grid finite-difference and the PML boundary［J］.Chinese Science Bulletin，51（17）：1985-1994（in Chinese）.

董良國(guó).2003.地震波數(shù)值模擬與反演中幾個(gè)關(guān)鍵問題研究［D］.上海：同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院：1-2.

Dong L G.2003.Several Key Problems in Seismic Wave Numerical Simulation and Inversion［D］.Shanghai：School of Ocean and Earth Science，Tongji University：1-2（in Chinese）.

何燕.2008.正交各向異性彈性波高階有限差分正演模擬研究［D］.東營(yíng)：中國(guó)石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院：45-80.

He Y.2008.High-Order Finite-Difference Forward Modeling of Elastic-Wave in Orthorhombic Anisotropic Media［D］.Dongying：School of Geosciences，China University of Petroleum：45-80（in Chinese）.

李敏，劉洋.2012.高階旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分方法模擬TTI介質(zhì)中橫波分裂［J］.物探與化探，36（6）：934-940.

Li M，Liu Y.2012.Modeling of the S-wave splitting in TTI media using high-order rotated staggered grid scheme［J］.Geophysical and Geochemical Exploration，36（6）：934-940（in Chinese）.

裴正林，王尚旭.2005.任意傾斜各向異性介質(zhì)中彈性波波場(chǎng)交錯(cuò)網(wǎng)格高階有限差分法模擬［J］.地震學(xué)報(bào)，27（4）：441-451.

Pei Z L，Wang S X.2005.A staggered-grid high-order finite-difference modeling for elastic wave field in arbitrary tilt anisotropic media［J］.Acta Seismologica Sinica，27（4）：441-451（in Chinese）.

孫成禹，肖云飛，印興耀，彭洪超.2010.黏彈介質(zhì)波動(dòng)方程有限差分解的穩(wěn)定性研究［J］.地震學(xué)報(bào)，32（2）：147-156.

Sun C Y，Xiao Y F，Yin X Y，Peng H C.2010.Study on the stability of finite difference solution of visco-elastic wave equations［J］.Acta Seismologica Sinica，32（2）：147-156（in Chinese）.

孫林潔.2012.非均勻孔隙介質(zhì)有限差分正演模擬方法研究［D］.青島：中國(guó)石油大學(xué) （華東）地球科學(xué)與技術(shù)學(xué)院：64-108.

Sun L J.2012.Finite-Difference Modeling Research in Heterogeneous Porous Media［D］.Qingdao：School of Geosciences，China University of Petroleum：64-108（in Chinese）.

王守東.2003.聲波方程完全匹配層吸收邊界［J］.石油地球物理勘探，38（1）：31-34.

Wang S D.2003.Absorbing boundary condition for acoustic wave equation by perfectly matched layer［J］.Oil Geophysi-cal Prospecting，38（1）：31-34（in Chinese）.

王書強(qiáng)，楊頂輝，楊寬德.2002.彈性波方程的緊致差分方法［J］.清華大學(xué)學(xué)報(bào)：自然科學(xué)版，42（8）：1128-1131.

Wang S Q，Yang D H，Yang K D.2002.Compact finite difference scheme for elastic equations［J］.Journal of Tsinghua University：Science and Technology，42（8）：1128-1131（in Chinese）.

嚴(yán)紅勇，劉洋.2012.黏彈TTI介質(zhì)中旋轉(zhuǎn)交錯(cuò)網(wǎng)格高階有限差分?jǐn)?shù)值模擬［J］.地球物理學(xué)報(bào)，55（4）：1354-1365.

Yan H Y，Liu Y.2012.Rotated staggered grid high-order finite-difference numerical modeling for wave propagation in viscoelastic TTI media［J］.Chinese Journal of Geophysics，55（4）：1354-1365（in Chinese）.

殷文，印興耀，吳國(guó)忱，梁鍇.2006.高精度頻率域彈性波方程有限差分方法及波場(chǎng)模擬［J］.地球物理學(xué)報(bào)，49（2）：561-568.

Yin W，Yin X Y，Wu G C，Liang K.2006.The method of finite difference of high precision elastic wave equations in the frequency domain and wave field simulation［J］.Chinese Journal of Geophysics，49（2）：561-568（in Chinese）.

Boersma B J.2011.A 6th order staggered compact finite difference method for the incompressible Navier-Stokes and scalar transport equations［J］.J Comput Phys，230（12）：4940-4954.

Collino F，Tsogka C.2001.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media［J］.Geophysics，66（1）：294-307.

Dai N，Vafidis A，Kanasewich E R.1995.Wave propagation in heterogeneous，porous media：A velocity-stress，finitedifference method［J］.Geophysics，60（2）：327-340.

Du Q Z，Li B，Hou B.2009.Numerical modeling of seismic wavefields in transversely isotropic media with a compact staggered-grid finite difference scheme［J］.Appl Geophys，6（1）：42-49.

Gold N，Shapiro S A，Burr E.1997.Modelling of high contrasts in elastic media using a modified finite difference scheme［C］∥SEG Technical Program Expanded Abstracts.Dallas，Texas：SEG：1850-1853.

Holberg O.1987.Computational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena［J］.Geophys Prospect，35（6）：629-655.

Igel H，Mora P，Riollet B.1995.Anisotropic wave propagation through finite-difference grids［J］.Geophysics，60（4）：1203-1216.

Kim J W，Lee D J.1996.Optimized compact finite difference schemes with maximum resolution［J］.AIAA J，34（5）：887-893.

Kirkpatrick S，Gelatt C D Jr，Vecchi M P.1983.Optimization by simulated annealing［J］.Science，220（4598）：671-680.

Kosloff D，Pestana R C，Tal-Ezer H.2010.Acoustic and elastic numerical wave simulations by recursive spatial derivative operators［J］.Geophysics，75（6）：T167-T174.

Lele S K.1992.Compact finite difference schemes with spectral-like resolution［J］.J Comput Phys，103（1）：16-42.

Liu Y，Sen M K.2009.An implicit staggered-grid finite-difference method for seismic modelling［J］.Geophys J Int，179（1）：459-474.

Liu Y.2013.Globally optimal finite-difference schemes based on least squares［J］.Geophysics，78（4）：T113-T132.

Saenger E H，Bohlen T.2004.Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid［J］.Geophysics，69（2）：583-591.

Saenger E H，Gold N，Shapiro S A.2000.Modeling the propagation of elastic waves using a modified finite-difference grid［J］.Wave Motion，31（1）：77-92.

Shan G.2009.Optimized implicit finite-difference and Fourier finite-difference migration for VTI media［J］.Geophysics，74（6）：WCA189-WCA197.

Virieux J.1984.SH-wave propagation in heterogeneous media：Velocity-stress finite-difference method［J］.Geophysics，49（11）：1933-1957.

Zeng Y Q，He J Q，Liu Q H.2001.The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media［J］.Geophysics，66（4）：1258-1266.

Zhang J H，Yao Z X.2012.Optimized finite-difference operator for broadband seismic wave modeling［J］.Geophysics，78（1）：A13-A18.

Zhou H B，Zhang G Q.2011.Prefactored optimized compact finite-difference schemes for second spatial derivatives［J］.Geophysics，76（5）：WB87-WB95.