周竹生，王 紅

(1. 中南大學 有色金屬成礦預測教育部重點實驗室，長沙 410083；2. 中南大學 地球科學與信息物理學院，長沙 410083)

改進BISQ模型地震波場數(shù)值模擬中的邊界處理

周竹生1,2，王 紅1,2

(1. 中南大學 有色金屬成礦預測教育部重點實驗室，長沙 410083；2. 中南大學 地球科學與信息物理學院，長沙 410083)

從改進BISQ模型雙相介質(zhì)所對應的一階速度—應力運動方程出發(fā)，構(gòu)建2×2N階交錯網(wǎng)格有限差分模擬算法，為了盡可能地減小或消除數(shù)值模擬中由人工邊界引起的虛假反射，建立完全匹配層(PML)吸收邊界的2×2N階交錯網(wǎng)格有限差分算法。詳細地討論了PML吸收邊界條件的構(gòu)建及其有限差分算法的實現(xiàn)。通過MATLAB編程進行波場模擬，將加入PML吸收邊界、常規(guī)指數(shù)衰減吸收邊界及未加吸收邊界的3種數(shù)值模擬結(jié)果進行對比，論證PML吸收邊界能十分有效地吸收邊界反射。

改進BISQ模型；雙相介質(zhì)；有限差分；交錯網(wǎng)格；完全匹配層(PML)

相對單相介質(zhì)而言，雙相介質(zhì)能夠更好地描述實際地層結(jié)構(gòu)和性質(zhì)，因此，雙相介質(zhì)的數(shù)值模擬更有利于實際介質(zhì)中的地震波傳播規(guī)律的認識和研究[1-4]。在利用計算機進行數(shù)值模擬計算過程中，由于計算區(qū)域有限，在波動方程數(shù)值模擬中常出現(xiàn)的問題就是來自邊界的人工反射，這種反射是由縮小計算區(qū)域而引起的。避免邊界反射的最簡單的方法是增大計算區(qū)域，但這種方法會延遲邊界反射和超出模擬的最大時間，顯然增加了計算量[5]。因此，出現(xiàn)了其它避免邊界反射的方法，例如吸收邊界條件。

好的吸收邊界條件不僅要吸收效果好，而且邊界應能盡量貼近主要計算區(qū)域，這樣可以節(jié)約計算機存儲空間和計算時間。計算吸收邊界的方法有很多種。最早的吸收邊界條件是由Lysmer和Kuhlemeyer提出的[6]，常被稱為黏性邊界條件，它的缺點是吸收人工反射波的效果較差。使用最為普遍的是CLAYTON和ENGQUIST[7]，ENGQUIST 和 MAJDA[8]提出的二維吸收邊界條件，他們用擬微分算子的工具將波動算子作分解，在人工邊界上保留向外傳播的部分，然后用Pade逼近，得到不同精度的邊界微分方程。這種邊界條件能完全吸收垂直入射時的反射波，且二階以上的Clayton-Engquist條件能較好吸收入射角在一定范圍內(nèi)的反射波，但當入射角接近π/2時，吸收效果很不好。雖然旁軸近似法[9]實現(xiàn)起來比較簡單且需要的網(wǎng)格數(shù)較少，能完全吸收垂直入射時的反射波，同時能較好吸收入射角在一定范圍內(nèi)的反射波，但它只在一定的角度和頻率范圍內(nèi)是有效吸收的。

針對上述不足，本文作者擬采用完全匹配層(PML)吸收邊界法來進行邊界處理，并通過計算實例，將其與傳統(tǒng)的指數(shù)衰減法進行對比，據(jù)此說明 PML方法的有效性。

1 數(shù)值模擬中的交錯網(wǎng)格高階有限差分

在改進BISQ模型的雙相介質(zhì)地震波場數(shù)值模擬中采用的一階速度—應力方程[10]為

稱為Biot流系數(shù)；Ks為固體體積模量；流體體積模量Kf可由流體聲速來確定；α為有效應力的孔隙彈性系數(shù)；vx、vz、Vx和 Vz分別表示固體和流體的位移分量；xxσ、zzσ分別是固體x、z方向的正應力；xzσ是固體的切應力；sρ和fρ分別為固體和流體密度；φ為孔隙度；η為流體黏滯系數(shù)；aρ為固-流耦合附加密度；μ和λ為彈性常數(shù)；k為滲透率；P為流體壓力。

對方程(1)進行高階差分近似，所構(gòu)建的空間 2N階、時間2階精度的有限差分格式如方程(2)所示，其他方程的差分格式與此類似。其中i、j、k表示x、z、t的離散序號；U、W、P、Q和S分別表示固相vx、vz、σxx、σzz和 σxz；u、w、F 和 R 分別表示流相 Vx、Vz、S 和 P。

2 邊界條件

2.1 指數(shù)衰減邊界條件

指數(shù)衰減的原理是對網(wǎng)格周圍的耗散采用質(zhì)點的速度和應力值同乘上一個小于1的函數(shù)來衰減。例如，采用CERJAN等[11]建議的吸收邊界，則吸收衰減函數(shù)定為

式中：I為給定的吸收邊界帶寬度的節(jié)點數(shù)；i為吸收邊界內(nèi)的節(jié)點號；α為衰減系數(shù)，其值的選定與I的大小密切相關，且對吸收效果的影響很大。這種吸收衰減方法的特點是對每個方向的反射波都采取相同的處理方式，實施方便簡單，計算所占內(nèi)存小。缺點是給定的吸收介質(zhì)的邊界要寬，不允許靠近邊界的區(qū)域出現(xiàn)介質(zhì)的不連續(xù)性，計算精度不高，吸收效果較一般。

2.2 PML邊界條件

完全匹配層(PML)的概念由Berenger于1994年首先提出[12]，它是目前消除邊界反射的理想方法之一。最初它被運用于時域電磁場有限差分邊界問題的解決中，后來人們不斷將其推廣到其他領域，并在聲波與彈性波的地震波數(shù)值模擬中得到成功應用[13-16]。PML吸收邊界條件假設存在一種各向異性有耗介質(zhì)(PML層)，地震波進入PML層后迅速衰減，通過適當選取參數(shù)，能夠使以不同入射角進入邊界層的任意頻率和任意極化的地震波保持其特征阻抗和相速度均固定，從而使地震波能夠完全匹配自由空間，理論上不產(chǎn)生任何反射[17-19]。

下面運用分裂式方法來實現(xiàn)完全匹配層(PML)思想，并給出相應的差分格式。

圖1所示為PML吸收邊界的區(qū)域劃分示意圖，其中，區(qū)域A1、A2、A3、A4分別為左邊界、右邊界、頂邊界和底邊界，區(qū)域 B1、B2、B3和 B4分別代表四個角區(qū)域，區(qū)域C是模擬工作的主要計算區(qū)域。

圖1 PML吸收邊界區(qū)域劃分示意圖Fig. 1 PML absorbing boundary zoning diagram

PML介質(zhì)中存在8個分量，在進行邊界條件處理時，需要對每個方向上的分量沿界面垂直和平行兩個方向進行分裂，為了消除人工分界面上的反射，在界面的垂直方向上需引入吸收衰減系數(shù)。

在二維空間條件下，對一階速度—應力彈性波動方程運用時域變量分裂法，對 PML介質(zhì)中的每個分量進行波場變量分離，并對于每個波場變量分裂如下：

展開為

式中：上標x和z代表相應的空間導數(shù)。根據(jù)上式對方程組(3)進行分裂，得到帶有衰減因子的 PML吸收邊界系統(tǒng)方程組的差分格式如下：

式中：Vmax代表PML介質(zhì)層的速度(最大縱波速度)；α= 10-6代表理論反射系數(shù)；L代表匹配層寬度；xi為匹配層區(qū)域與內(nèi)部區(qū)域界面的橫向距離；zi為匹配層區(qū)域與內(nèi)部區(qū)域界面的縱向距離；系數(shù)a=0.25，b=0.75。

數(shù)值模擬計算時，衰減系數(shù)的取值根據(jù)其所在的區(qū)域變化而改變，在不同區(qū)域其取值不同。根據(jù)圖1，在左邊界A1和右邊界A2區(qū)域：＞0，=0；在上邊界A3和下邊界A4區(qū)域：=0，＞0；在四個角B1、B2、B3、B4區(qū)域：＞0，＞0。

3 波場模擬

為檢驗模擬算法的正確性，本文作者建立了一個各向同性雙相介質(zhì)模型，模型網(wǎng)格點數(shù)為200×200，網(wǎng)格間距為5 m，時間間隔為0.1 ms，震源位于模型中心，采用主頻為30 Hz的雷克子波，模型的物性參數(shù)設置如下：固體參數(shù)有拉梅系數(shù)12.72 GPa，剪切模量 6.84 GPa，固體體積模量 56.81 GPa，密度 2 738 kg/m3，耦合附加密度83 kg/m3，縱波速度3 000 m/s，橫波速度2 000 m/s；流體參數(shù)有體積模量1.46 GPa，密度454 kg/m3，孔隙度0.237 8，粘滯系數(shù)1×10-6Pa·s，滲透率10 md。采用精度為O(Δt2+Δx16)的交錯網(wǎng)格有限差分法，對固體和流體的水平分量和垂直分量分別進行計算，數(shù)值模擬的結(jié)果如圖2所示。

由圖2可以看出，在雙相介質(zhì)中存在第二類縱波，從內(nèi)到外依次為慢縱波(p)，橫波(S)和快縱波(P)，其衰減規(guī)律明顯不同于第一類縱波，第一類縱波在固相和流相都引起較強的振動，第二類縱波引起較弱的固相振動，較強的流相振動，這種關系取決于雙相介質(zhì)的性質(zhì)。其中第一類縱波的相位在固相和流相部分是相同的，第二類縱波的相位在固相和流相部分則相反。

4 邊界處理效果對比

根據(jù)以上模型，對固體和流體的水平分量和垂直分量分別進行計算，并采用不同方法對人工邊界進行處理，所得數(shù)值模擬結(jié)果如下(見圖3～5)：

圖3所示為對人工邊界未進行任何處理時所得到的150 ms時刻的波場快照。由圖3可見，在邊界處，產(chǎn)生了明顯的邊界反射波，必將嚴重干擾計算波場的記錄特征。

圖4所示為對人工邊界采用了指數(shù)衰減處理時所得到的150 ms時刻的波場快照，此時，邊界反射波明顯減弱，有效波場得到突出，但邊界反射很難吸收干凈，對計算波場仍存影響。

圖5所示為對人工邊界采用了PML吸收邊界條件處理后所得到的 150 ms時刻的波場快照。從圖 5可清楚看出，加了 PML吸收邊界以后，由于該邊界條件很好地削弱了來自人工邊界的反射，使得模擬波場的特征變得非常清晰。

圖2 100 ms時刻的波場快照(時間2階，空間16階)Fig. 2 Snapshots of 100 ms (two-order time and sixteen-order space): (a) Solid phase horizontal components; (b) Solid phase vertical components; (c) Flow phase horizontal components; (d) Flow phase vertical components

圖3 未加邊界條件時150 ms時刻的波場快照(時間2階，空間16階)Fig. 3 Snapshots of 150 ms (two-order time and sixteen-order space) with non-absorbing boundary: (a) Solid phase horizontal components; (b) Solid phase vertical components; (c) Flow phase horizontal components; (d) Flow phase vertical components

圖4 加指數(shù)衰減邊界時150 ms時刻的波場快照(時間2階，空間16階)Fig. 4 Snapshots of 150 ms (two-order time and sixteen-order space) with exponential boundary: (a) Solid phase horizontal components; (b) Solid phase vertical components; (c) Flow phase horizontal components; (d) Flow phase vertical components

圖5 加PML邊界條件時150 ms時刻的波場快照Fig. 5 Snapshots of 150 ms (two-order time and sixteen-order space) with PML: (a) Solid phase horizontal components; (b) Solid phase vertical components; (c) Flow phase horizontal components; (d) Flow phase vertical components

5 結(jié)論

1) 將交錯網(wǎng)格有限差分法運用于改進BISQ模型的雙相介質(zhì)，給出了其相應的一階速度—應力波動方程在不同時間和空間精度條件下的交錯網(wǎng)格差分算法格式，實現(xiàn)了雙相介質(zhì)的地震波場數(shù)值模擬。模擬波場符合雙相介質(zhì)的理論特征，表明所建立的算法格式是正確、可行的。

2) 針對波場數(shù)值模擬過程中的邊界處理問題，建立了完全匹配層(PML)邊界的2×2N階交錯網(wǎng)格高階有限差分格式，通過與未加吸收邊界及加常規(guī)指數(shù)衰減吸收邊界所得到的模擬波場進行對比，說明 PML邊界條件能十分有效地吸收邊界反射，驗證了 PML吸收邊界條件的高效、優(yōu)越和可行性。

REFERENCES

[1] BIOT M A. Theory of propagation of elastic waves in a fluidsaturated porous solid (Part I): Low-frequency range [J]. The Journal of the Acoustical Society of America, 1956, 28(2): 168-178.

[2] BIOT M A. Theory of propagation of elastic waves in a fluid-saturated porous solid (Part II): Higher frequency range [J].The Journal of the Acoustical Society of America, 1956, 28(2):179-191.

[3] DVORKIN J, NUR A. Dynamic poroelasticity: A unified model with the squirt and the Biot mechanisms [J]. Geophysics, 1993,58(4): 524-533.

[4] DIALLO M S, APPEL E. Acoustic wave propagation in saturated porous media: reformulation of the Biot/Squirt flow theory [J]. Journal of Applied Geophysics, 2000, 44(4):313-325.

[5] LIU Y, SEN M K. A hybrid scheme for absorbing edge reflections in numerical modeling of wave propagation [J].Geophysics, 2010, 75(2): A1-A6.

[6] DAI N, VAFIDIS A, KANASEWICH E. Composite absorbing boundaries for the numerical simulation of seismic waves [J].Bulletin of the Seismological Society of America, 1994, 84(1):185-191.

[7] CLAYTON R, ENGQUIST B. Absorbing boundary conditions for acoustic and elastic wave equations [J]. Bulletin of the Seismological Society of America, 1977, 67(6): 1529-1540.

[8] ENGQUIST B, MAJDA A. Radiation boundary conditions for acoustic and elastic wave calculations [J]. Communications on Pure and Applied Mathematics, 1979, 32(3): 313-357.

[9] 夏 凡, 董良國, 馬在田. 三維彈性波數(shù)值模擬中的吸收邊界條件[J]. 地球物理學報, 2004, 47(1): 132-136.XIA Fan, DONG Liang-guo, MA Zai-tian. Absorbing boundary conditions in 3D elastic-wave numerical modeling [J]. Chinese Journal of Geophysics, 2004, 47(1): 132-136.

[10] 王 俊. 基于改進 BISQ模型的雙相介質(zhì)波場數(shù)值模擬方法研究[D]. 山東: 中國石油大學(華東), 2009: 41-45.WANG Jun. Method study of wavefield numerical simulation in two-phase medium based on the reformulated BISQ mechanism[D]. Shandong: China University of Petroleum (East China),2009: 41-45.

[11] CERJAN C, KOSLOFF D, KOSLOFF R, RESHEF M. A nonreflecting boundary condition for discrete acoustic and elastic wave equations [J]. Geophysics, 1985, 50(4): 705-708.

[12] BERENGER J P. A perfectly matched layer for the absorption of electromagnetic waves [J]. Journal of Computational Physics,1994, 114(2): 185-200.

[13] CHEW W C, LIU Q H. Perfectly matched layers for elastodynamics: a new absorbing boundary condition [J]. Journal of Computational Acoustics, 1996, 4(4): 341-359.

[14] LIU Q H, TAO J P. The perfectly matched layer for acoustic waves in absorptive media [J]. Journal of the Acoustical Society of America, 1997, 102(4): 2072-2082.

[15] ZENG Y Q, HE J Q, LIU Q H. The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media [J]. Geophysics, 2001, 66(4): 1258-1266.

[16] COLLINO F, TSOGKA C. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media [J]. Geophysics, 2001, 66(1):294-307.

[17] 王永剛, 邢文軍, 謝萬學, 朱兆林. 完全匹配層吸收邊界條件的研究[J]. 中國石油大學學報: 自然科學版, 2007, 31(1):19-24.WANG Yong-gang, XING Wen-jun, XIE Wan-xue, ZHU Zhao-lin. Study of absorbing boundary condition by perfectly matched layer [J]. Journal of China University of Petroleum:Edition of Natural Science, 2007, 31(1): 19-24.

[18] 趙海波, 王秀明, 王 東, 陳 浩. 完全匹配層吸收邊界在孔隙介質(zhì)彈性波模擬中的應用[J]. 地球物理學報, 2007, 50(2):581-591.ZHAO Hai-bo, WANG Xiu-ming, WANG Dong, CHEN Hao.Applications of the boundary absorption using a perfectly matched layer for elastic wave simulation in poroelastic media[J].Chinese Journal of Geophysics, 2007, 50(2): 581-591.

[19] 熊章強, 張大洲, 秦 臻, 周文斌. 瑞雷波數(shù)值模擬中的邊界條件處理及模擬實例分析[J]. 中南大學學報:自然科學版,2008, 39(4): 824-830.XIONG Zhang-qiang, ZHANG Da-zhou, QIN Zheng, ZHOU Wen-bin. Boundary conditions and case analysis of numerical modeling of Rayleigh wave [J]. Journal of Central South University: Science and Technology, 2008, 39(4): 824-830.

[20] 李 寧. 完美匹配層理論及其在地震波場數(shù)值模擬中的應用[D]. 哈爾濱: 中國地震局工程力學研究所, 2006: 46-49.LI Ning. The theory and the application of perfectly matched layer in the seismic wave simulation[D]. Harbin: Institute of Engineering Mechanics, China Earthquake Administration, 2006:46-49.

Boundary treatment for numerical simulation of seismic waves based on reformulated BISQ model

ZHOU Zhu-sheng1,2, WANG Hong1,2
(1. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education,Central South University, Changsha 410083, China;2. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)

Based on first-order velocity-stress equations of the reformulated BISQ model in double-phase media, the simulation algorithm of finite difference of 2-order time and 2N-order space staggered-grid was built. Meanwhile, in order to minimum effects caused by the artificial boundary in the numerical simulation, the construction of the perfectly matched layer (PML) absorbing boundary condition and the realization of the finite-difference algorithm were discussed in detail. The arithmetic was realized with MATLAB. Compared with the conventional decaying exponential absorbing boundary and the non-absorbing boundary, the PML boundary can more effectively attenuate reflections, which is supported by the wave field modeling.

reformulated BISQ model; double-phase medium; finite-difference method; staggered-grid; perfectly matched layer (PML)

P631.4

A

1004-0609(2012)03-0976-09

國家自然科學基金資助項目(41074085，40804027)；湖南省自然基金重點資助項目(09JJ3084)

2011-12-01；

2012-01-04

周竹生，教授，博士；電話：13707316924；E-mail: 672390136@qq.com

(編輯 何學鋒)