彈性波傳播的低數(shù)值頻散波場模擬方法

張朝元

（大理學(xué)院 數(shù)學(xué)與計算機(jī)學(xué)院，云南 大理671003）

研究地震波在地球介質(zhì)中的傳播規(guī)律是認(rèn)識地球內(nèi)部結(jié)構(gòu)和進(jìn)行油氣資源勘探等研究活動的重要內(nèi)容。波動方程數(shù)值模擬方法是了解地震波傳播規(guī)律的一項(xiàng)基本手段，為地震波的傳播機(jī)理研究提供了重要依據(jù)。彈性波方程是一種常見的地震波方程，被廣泛地應(yīng)用于地震波研究中。準(zhǔn)確正演模擬彈性波傳播是研究地震各向異性和AVO反演的基礎(chǔ)。因此，本文基于彈性波方程，研究地震波傳播的正演數(shù)值模擬方法。

常見的正演數(shù)值模擬方法依據(jù)微分項(xiàng)處理的數(shù)值方式而不同，如射線追蹤法［1］、有限元法［2］、偽譜法［3］、反射率法［4］和有限差分法［5－7］等。利用這些方法求解波動方程數(shù)值解時，會遇到因離散求解波動方程而產(chǎn)生的偽波動現(xiàn)象，即所謂的網(wǎng)格頻散或數(shù)值頻散，尤其是在有較大梯度的波場附近或網(wǎng)格點(diǎn)較粗時，頻散現(xiàn)象更為嚴(yán)重，這樣極大地降低了地震波模擬精度。為此，人們開展了大量的研究工作并從頻散前校正和頻散后校正提出大量壓制數(shù)值頻散的數(shù)值模擬方法［8－16］。頻散前校正方法以高階數(shù)值格式和新型高階微分算子為代表。高階數(shù)值格式［6，7］不一定意味著會降低數(shù)值頻散［9］，且存在邊界處理困難和不易并行等缺點(diǎn)。以近似解析離散化算子［17］（記為NAD）為代表的新型高階微分算子［18－22］包含位移和梯度等信息，因此，NAD算子不僅能有效地壓制數(shù)值頻散且適應(yīng)于大規(guī)模并行計算。頻散后校正方法主要有通量校正傳輸技術(shù)［10］（記為FCT）。Boris和Book等人把用于流體動力學(xué)數(shù)值仿真問題的通量校正傳輸技術(shù)推廣應(yīng)用于地震波數(shù)值模擬上，取得了許多有價值的結(jié)果［8－16］。

為有效地壓制數(shù)值頻散以提高地震波場模擬精度，在前人的基礎(chǔ)上，本文提出了彈性波傳播的八階FNRK方法。該方法是將八階NAD算子［23－25］替代差分算子進(jìn)行頻散前校正而采用FCT技術(shù)［10］作為頻散后校正的一種 Runge－Kutta方法［26］?；陔p層均勻各向同性介質(zhì)和三層橫向各向同性介質(zhì)的地震波數(shù)值實(shí)驗(yàn)表明，同八階 NAD-RK方法、八階Lax-Wendroff（LWC）方法和八階交錯網(wǎng)格（SG）方法相比，本文提出的八階FNRK方法能有效地壓制數(shù)值頻散并最大程度地提高數(shù)值模擬記錄的信噪比，能有效地應(yīng)用于彈性波復(fù)雜介質(zhì)模擬研究。

1 方法原理

下面，我們以二維彈性波方程來闡述本文提出的八階FNRK方法的基本原理。設(shè)二維彈性波傳播方程如下

其中：下標(biāo)j＝1，3；ρ＝ρ（x，z）和σij分別代表介質(zhì)密度和應(yīng)力，ui和fi分別為i方向的位移分量和力源函數(shù)分量。在地震波場模擬實(shí)驗(yàn)中，通常選擇隨時間變化的震源函數(shù)為

1.1 空間網(wǎng)格離散

為了有效壓制因空間離散化所產(chǎn)生的數(shù)值頻散，我們引入近似解析離散算子來替代差分算子進(jìn)行空間偏導(dǎo)數(shù)高階離散。

其中G是二階微分算子。

再令vi＝?ui/?t（i＝1，2，3），V＝（v1，v2，v3）T，則方程（3）式變?yōu)?/p>

又令W＝（U，V）T，則式（4）變?yōu)?/p>

這樣，我們得到了與彈性波傳播方程（1）式等價的方程組（5）式，（5）式中的微分算子L包含了關(guān)于空間的二階和三階偏導(dǎo)數(shù)。在這里，我們選擇利用近似解析離散算子來對空間偏導(dǎo)數(shù)進(jìn)行八階離散，具體的離散公式［23－25］見附錄 A。

1.2 時間離散

為求解具有常微分方程組形式的彈性波傳播方程等價方程（5）式，我們采用三級三階的Runge-Kutta方法［26］來對時間導(dǎo)數(shù)進(jìn)行離散，簡化的計算公式如下

其中Wn＝W（nΔt）。把向量W＝（U，V）T代入（6）式，就有

其中偏微分算子G2＝G·G。

把用于求解彈性波傳播方程（1）式的離散公式（6）式或（7）式簡稱為八階 NAD-RK方法。

1.3 八階FNRK方法

由（6）式或（7）式模擬地震彈性波在介質(zhì)模型中的傳播，能得到數(shù)值模型中所有網(wǎng)格點(diǎn)在任意時刻各方向的位移值和粒子速度值。為在數(shù)值模擬中有效消除數(shù)值頻散的影響以提高計算精度，我們引入通量校正傳輸技術(shù)［10］。采用FCT技術(shù)來校正八階NAD-RK方法模擬彈性波傳播后的數(shù)值頻散，即得到基于彈性波傳播方程的低數(shù)值頻散波場模擬方法，稱為八階FNRK方法。該方法包括下面5個步驟。

（1）給定初始值（其 中

（2）利用離散公式（6）式或（7）式來計算（其中n≥0）。

（3）漫射通量計算。

①計算第n個和n＋1個時間步長的漫射通量為

其中：λ，μ（0≤λ，μ≤1）為漫散通量校正系數(shù)。

②在壓制網(wǎng)格頻散引起的偽波動同時也導(dǎo)致了一定的波場振幅損失，為此利用漫射通量A，B來平滑離散方程（6）式或（7）式的解

（4）反漫射通量計算過程。

①用被修改過的解和八階 NAD-RK 方法模擬得到的解來計算反漫射通量，即

②利用反漫射通量修正平滑解，即可得到消除數(shù)值頻散后的正確解

其中：

這里，m（1，2）表示向量X，Z，Sx，Sz中的第m個元素，max代表取最大值，min代表取最小值，sign代表取數(shù)學(xué)符號，abs表示取絕對值。

（5）接著回到步驟（2），重復(fù)步驟（2）至步驟（4），如此循環(huán)直到滿足結(jié)束條件，即可得到所需的結(jié)果。

2 數(shù)值模擬

為驗(yàn)證本文提出的八階FNRK方法在彈性波模擬中具有低數(shù)值頻散，我們用八階FNRK方法來模擬彈性波在均勻各向同性（記為HI）介質(zhì)和橫向各向同性（記為TI）介質(zhì)中的傳播。

2.1 雙層HI介質(zhì)

地震波在介質(zhì)不連續(xù)面發(fā)生反射、折射，并攜帶著地質(zhì)構(gòu)造信息的上行波被地震臺站接收而形成所謂的地震記錄?；谶@些地震記錄，由地震波反演成像方法推斷出地下結(jié)構(gòu)，從而引導(dǎo)人們對深層內(nèi)部結(jié)構(gòu)進(jìn)行認(rèn)識。因此，準(zhǔn)確模擬地震波在地下介質(zhì)波阻抗不連續(xù)處的反射、折射具有重要意義。而雙層介質(zhì)模型是多層介質(zhì)模型的基礎(chǔ)，為此，我們首先考察八階FNRK方法模擬彈性波在雙層均勻各向同性介質(zhì)中的傳播。

設(shè)雙層HI介質(zhì)模型的計算區(qū)域?yàn)?≤x，z≤9km，層間界面位于深度z0≤4.5km處。上層介質(zhì)中的彈性系數(shù)和密度分別為λ1＝4.8 GPa，μ1＝6.5GPa，ρ1＝3.2g/cm3；下層介質(zhì)中的彈性系數(shù)和密度分別為λ2＝6.0GPa，μ2＝24.0GPa，ρ2＝1.5g/cm3。震源位于模型中心，震源函數(shù)為f0＝25Hz的Ricker子波，且f1（x）＝f3（x）＝f（t）。接收器位置為R（5km，4km）。設(shè)時間步長為Δt＝3ms，空間步長為Δx＝Δz＝30m。

圖1分別為八階FNRK方法（圖1-A）、八階NAD-RK方法（圖1-B）、八階LWC方法（圖1-C）和八階SG方法（圖1-D）模擬得到t＝1.2s時刻在水平位移方向的波場瞬時快照。圖2分別為由八階FNRK 方法（圖2-A）、八階 NAD-RK 方法（圖2-B）、八階 LWC方法（圖2-C）和八階SG 方法（圖2-D）模擬得到在R處從時間t＝0s到t＝1.2s在水平位移方向的波形記錄。從圖1和圖2可看出，由八階FNRK方法模擬得到的波場快照和波形記錄無任何可見的數(shù)值頻散，而同階的NAD-RK方法、LWC方法和SG方法卻有著嚴(yán)重的數(shù)值頻散。這說明八階FNRK方法在粗網(wǎng)格條件下對于彈性波在強(qiáng)間斷層狀介質(zhì)中有著很好的波傳播模擬效果。

2.2 三層TI介質(zhì)

為檢驗(yàn)八階FNRK方法的模擬效果和優(yōu)越性，我們選擇三層TI介質(zhì)模型進(jìn)行波場模擬，同時與八階NAD-RK方法、八階LWC方法和八階SG方法進(jìn)行比較。其中，三層TI介質(zhì)模型中的彈性系數(shù)和介質(zhì)密度如表1所示。設(shè)模型的計算區(qū)域?yàn)?≤x，z≤12km，采用頻率f0＝20Hz的Ricker子波作為震源，且f1＝f3＝f（t），震源位于計算區(qū)域中心，取空間步長為Δx＝Δz＝40m，取時間步長為Δt＝2ms。

圖3給出了粗網(wǎng)格（Δx＝Δz＝40m）情況下在t＝1.2s時刻分別由八階FNRK方法（圖3－A）、八階 NAD-RK 方法（圖3-B）、八階 LWC方法（圖3-C）和八階SG方法（圖3-D）模擬得到在水平位移方向的波場瞬時快照。從圖3可以看出，4種方法模擬得到的數(shù)值波場幾乎一樣，但是圖3-B表明八階NAD-RK方法有輕度的數(shù)值頻散誤差，而圖3-C和圖3-D表明八階LWC方法和八階SG方法均有著嚴(yán)重的數(shù)值頻散現(xiàn)象。而在同樣的計算條件下，本文提出的八階FNRK方法得到了清晰的波場快照（圖3-A），無可見的數(shù)值頻散。可見，在粗網(wǎng)格條件下八階FNRK方法在壓制數(shù)值頻散上明顯優(yōu)越于八階LWC方法和八階SG方法，也好于八階NAD-RK方法。

圖1 雙層介質(zhì)下t＝1.2s時刻水平位移方向的波場瞬時快照Fig.1 Snapshots of wave fields at time t＝1.2sand in the horizontal displacement direction under the two-layer media

圖2 雙層介質(zhì)下在接收點(diǎn)R（5km，4km）處的波形記錄Fig.2 Waveforms at the receiver point R（5km，4km）and in the horizontal displacement direction under the two-layer media

圖3 粗網(wǎng)格（Δx＝Δz＝40m）條件下t＝1.2s時刻水平位移方向的波場瞬時快照Fig.3 Snapshots of wave fields at time t＝1.2sin the horizontal displacement direction under the coarse grid（Δx＝Δz＝40m）

表1 三層介質(zhì)模型中的參數(shù)Table 1 Parameters used in a three-layer model

3 結(jié)論

如何有效地消除波場模擬中的數(shù)值頻散是地震波正演數(shù)值模擬方法研究的重要內(nèi)容之一。基于彈性傳播方程，本文以近似解析離散算子替代差分算子對空間偏導(dǎo)數(shù)進(jìn)行八階離散，以Runge－Kutta方法對時間導(dǎo)數(shù)進(jìn)行三階離散，從而實(shí)現(xiàn)頻散前的校正；而采用通量校正傳輸技術(shù)進(jìn)行頻散后校正，最終得到一種頻散前校正和頻散后校正有效結(jié)合的綜合校正數(shù)值頻散的方法，即八階FNRK方法。為考察本文所提出的八階FNRK方法在壓制數(shù)值頻散上的效果，我們將其應(yīng)用于模擬彈性波在雙層HI介質(zhì)和三層TI介質(zhì)中的傳播。模擬結(jié)果表明，該方法在壓制數(shù)值頻散上明顯優(yōu)越于傳統(tǒng)的高階有限差分方法。同時，該方法對在大尺度復(fù)雜介質(zhì)中的地震波傳播有著很強(qiáng)的模擬適應(yīng)能力。因此，本文提出的八階FNRK方法有望廣泛地用于大尺度復(fù)雜地質(zhì)構(gòu)造介質(zhì)中的地震波場模擬和分析。

［1］Chapman C H，Shearer P M.Ray tracing in azimuthally anisotropic media-Ⅱ：Quasi shear wave coupling［J］.Geophys J，1989，96：65－83.

［2］楊頂輝.雙相各向異性介質(zhì)中彈性波方程的有限元解法及波場模擬［J］.地球物理學(xué)報，2002，45（4）：575－583.Yang D H.Finite element method of the elastic wave equation and wave field simulation in two-phase anisotropic media［J］.Chinese Journal of Geophysics，2002，45（4）：575－583.（In Chinese）

［3］Komatitsch D，Vilotte J P.The spectral element method：An efficient tool to simulate the seismic responses of 2Dand 3Dgeological structures［J］.Bull Seism Soc Am，1998，88：368－392.

［4］Chen X F.A systematic and efficient method of computing normal modes for multi-layered half space［J］.Geophys J Int，1993，115：391－409.

［5］Kelly K R，Wave R W，Treitel S.Synthetic seismograms：A finite-difference approach［J］.Geophysics，1976，41：2－27.

［6］Dablain M A.The application of high-order differencing to scalar wave equation［J］.Geophysics，1986，51：54－66.

［7］董良國，馬在田，曹景忠，等.一階彈性波方程交錯網(wǎng)格高階差分解法［J］.地球物理學(xué)報，2000，43（3）：411－419.Dong L G，Ma Z T，Cao J Z，etal.A staggered-grid high-order difference method of one-order elastic wave equation［J］.Chinese Journal of Geophysics，2000，43（3）：411－419.（In Chinese）

［8］Book D L，Boris J P，Hain K.Flux-corrected transport，II：generalization of the method［J］.J Comput Phys，1975，18（3）：248－283.

［9］Fei T，Larner K.Elimination of numerical dispersion in finite-difference modeling and migration by fluxcorrected transport［J］.Geophysics，1995，60（6）：1830－1842.

［10］楊頂輝，滕吉文.各向異性介質(zhì)中三分量地震記錄的FCT有限差分模擬［J］.石油地球物理勘探，1997，32（2）：181－190.Yang D H，Teng J W.FCT finite difference modeling of three-component seismic records in anisotropic medium［J］.OGP，1997，32（2）：181－190.（In Chinese）

［11］Yang D H，Liu E，Zhang Z J，etal.Finitedifference modeling in two-dimensional anisotropic media using a flux-corrected transport technique［J］.Geophys J Int，2002，148：320－328.

［12］鄭海山，張中杰.橫向各向同性（VTI）介質(zhì)中非線性地震波場模擬［J］.地球物理學(xué)報，2005，48（3）：660－671.Zheng H S，Zhang Z J.Synthetic seismograms of nonlinear seismic waves in isotropic（VTI）media［J］.Chinese Journal of Geophysics，2005，48（3）：660－671.（In Chinese）

［13］吳國忱，王華忠.波場模擬中的數(shù)值頻散分析與校正策略［J］.地球物理學(xué)進(jìn)展，2005，20（1）：58－65.Wu G C，Wang H Z.Analysis of numerical dispersion in wave-field simulation［J］.Progress in Geophysics，2005，20（1）：58－65.（In Chinese）

［14］Zheng H S，Zhang Z J，Liu E.Non-linear seismic wave propagation in anisotropic media using fluxcorrected transport technique［J］.Geophys J Int，2006，165：943－956.

［15］王珺，楊長春，馮英杰.用優(yōu)化通量校正傳輸技術(shù)壓制數(shù)值模擬的頻散［J］.勘探地球物理進(jìn)展，2007，30（4）：252－256.Wang J，Yang C C，F(xiàn)eng Y J.Elimination of numerical dispersion in finite-difference modeling by optimized flux-corrected transport［J］.PEG，2007，30（4）：252－256.（In Chinese）

［16］陳可洋，楊微，劉洪林，等.二階彈性波動方程高精度交錯網(wǎng)格波場分離數(shù)值模擬［J］.物探與化探，2009，33（6）：700－704.Chen K Y，Yang W，Liu H L，etal.Wave field separation numerical modeling of second order elastic wave equation by high-precision staggered-grid finite difference scheme［J］.Geophysical ＆ Geochemical exploration，2009，33（6）：700－704.（In Chinese）

［17］Kondoh Y，Hosaka Y，Ishii K.Kernel optimum nearly-analytical discretization algorithm applied to parabolic and hyperbolic equations［J］.Computers Math Appl，1994，27（3）：59－90.

［18］Yang D H，Teng J W，Zhang Z J，etal.A nearlyanalytic discrete method for acoustic and elastic wave equations in anisotropic media［J］.Bull Seism Soc Am，2003，93（2）：882－890.

［19］Yang D H，Peng J M，Lu M，etal.Optimal nearlyanalytic discrete approximation to the scalar wave equation［J］.Bull Seism Soc Am，2006，96（3）：1114－1130.

［20］Yang D H，Song G J，Chen S，etal.An improved nearly analytical discrete method：an efficient tool to simulate the seismic response of 2-D porous structures［J］.J Geophys Eng，2007，4：40－52.

［21］Yang D H，Wang N，Chen S，etal.An explicit method based on the implicit Runge-Kutta algorithm for solving the wave equations［J］.Bull Seism Soc Am，2009，99（6）：3340－3354.

［22］王磊，楊頂輝，鄧小英.非均勻介質(zhì)中地震波應(yīng)力場的WNAD方法及其數(shù)值模擬［J］.地球物理學(xué)報，2009，52（6）：1526－1535.Wang L，Yang D H，Deng X Y.A WNAD method for seismic stress-field modeling in heterogeneous media［J］.Chinese J Geophys，2009，52（6）：1526－1535.（In Chinese）

［23］Tong P，Yang D H，Hua B L，etal.A high-order stereo-modeling method for solving wave equations［J］.Bulletin of the Seismological Society of America，2013，103（2）：811－833.

［24］Zhang C Y，Li X，Ma X，etal.A Runge-Kutta method with using eighth-order nearly-analytic spatial discretization operator for solving a 2Dacoustic wave equation［J］.Journal of Seismic Exploration，2014，23（3）：279－302.

［25］Zhang C Y，Ma X，Yang L，etal.Symplectic partitioned Runge-Kutta method based on the eighthorder nearly analytic discrete operator and its wavefield simulations［J］.Applied Geophysics，2014，11（1）：89－106.

［26］Qiu J X，Li T G，Khoo B C.Simulations of compressible two-medium flow by Runge-Kutta discontinuous Galerkin methods with the ghost fluid method［J］.Commun Computs Phys，2008，3：479－504.

附錄A：空間高階偏導(dǎo)數(shù)的逼近公式

利用局部插值近似方法，可以得到位移和粒子速度關(guān)于空間的二階和三階偏導(dǎo)數(shù)的逼近公式，具體推導(dǎo)過程見參考文獻(xiàn)［23－25］。這里，只給出了關(guān)于位移的高階偏導(dǎo)數(shù)的計算公式

其中：Δx，Δz分別代表x方向和z方向的空間步長。

把（A-1）～（A-7）式中的u換成v就得到關(guān)于粒子速度的二階和三階偏導(dǎo)數(shù)的逼近公式。