基于正交貼體網(wǎng)格的黏彈介質(zhì)地震波模擬

劉志強(qiáng)，黃磊，李鋼柱，牛興國(guó)，張曉萌

（1.內(nèi)蒙古農(nóng)業(yè)大學(xué)水利與土木建筑工程學(xué)院，內(nèi)蒙古呼和浩特 010018；2.內(nèi)蒙古有色地質(zhì)礦業(yè)（集團(tuán)）物探勘查公司，內(nèi)蒙古呼和浩特 010010）

0 引言

受到地層黏滯性的影響，地震波在地下介質(zhì)中傳播時(shí)能量會(huì)衰減，這種衰減在地表附近尤為嚴(yán)重。若用完全彈性模型模擬地下真實(shí)介質(zhì)，結(jié)果與實(shí)際差別較大。研究表明，黏彈模型可以較精確地反映地層的吸收作用，近似代替實(shí)際地下介質(zhì)更為合理。很多學(xué)者對(duì)黏彈介質(zhì)中的地震波場(chǎng)進(jìn)行了數(shù)值模擬。Robertsson等［1］設(shè)計(jì)了二維速度—應(yīng)力交錯(cuò)網(wǎng)格模擬方法，并推廣至三維黏彈模型［2］；Xu等［3-4］實(shí)現(xiàn)了三維起伏地表黏彈模型的數(shù)值模擬，并通過(guò)定義新的復(fù)合記憶變量有效節(jié)省了計(jì)算機(jī)內(nèi)存；?tekl等［5］在頻率域用旋轉(zhuǎn)有限差分算子對(duì)黏彈模型進(jìn)行了精確的數(shù)值模擬。在中國(guó)，宋常瑜等［6］利用有限差分法進(jìn)行了井間數(shù)值黏彈模擬，并分析了地震波在井間的傳播特征；唐啟軍等［7］利用有限差分法模擬了黏彈單斜各向異性介質(zhì)中的地震波場(chǎng)；陳玉璽等［8］研究了三維黏彈介質(zhì)地震波場(chǎng)數(shù)值模擬技術(shù)；陳豪等［9］利用交錯(cuò)網(wǎng)格有限差分法模擬了三維黏彈性介質(zhì)中的地震波場(chǎng)，并分析了品質(zhì)因子對(duì)波場(chǎng)特征的影響。

常用的數(shù)值模擬方法包括有限差分法［10-12］、偽譜法［13-14］、有限元法［15-16］、普元法［17-18］等。其中有限差分法因具有原理簡(jiǎn)單、易編程實(shí)現(xiàn)、計(jì)算效率高、對(duì)內(nèi)存需求相對(duì)較小等特點(diǎn)，在地震波數(shù)值模擬中使用最廣泛。常規(guī)有限差分法采用笛卡爾坐標(biāo)系下的規(guī)則網(wǎng)格，在模擬起伏地表時(shí)必然會(huì)出現(xiàn)階梯狀的邊界，從而引起虛假散射波。為減弱這種虛假散射波，需采用精細(xì)網(wǎng)格，這會(huì)導(dǎo)致存儲(chǔ)量的增加和計(jì)算量的增大。因此，發(fā)展了可變網(wǎng)格和不規(guī)則網(wǎng)格地震波數(shù)值模擬方法。Hayashi等［19］采用空間變網(wǎng)格法加密起伏地表附近的網(wǎng)格，從而盡量消除階梯狀網(wǎng)格帶來(lái)的誤差，同時(shí)能減小計(jì)算量。董良國(guó)［20］和Tarrass 等［21］利用縱向坐標(biāo)變換思路，將起伏邊界變換為水平邊界，然后在變換后的矩形網(wǎng)格中采用有限差分法求解地震波方程。但是他們的方法只適用于地形起伏較緩的情況，在起伏較大時(shí)容易出現(xiàn)不穩(wěn)定現(xiàn)象。Lan 等［22］、Sun等［23］在貼體網(wǎng)格下使用同位網(wǎng)格有限差分法求解波動(dòng)方程，能夠適應(yīng)任意起伏地表，計(jì)算精度較高，但為了滿足自由邊界條件需要做繁雜的坐標(biāo)旋轉(zhuǎn)和插值運(yùn)算。為此，蔣麗麗［24］、丘磊［25］、李慶洋等［26］提出采用正交貼體網(wǎng)格對(duì)起伏地表模型進(jìn)行網(wǎng)格剖分。

在對(duì)現(xiàn)有各種起伏地表地震波數(shù)值模擬方法進(jìn)行研究的基礎(chǔ)上，本文將Zhang 等［27］提出的改進(jìn)RL（Ryskin & Leal）正交貼體網(wǎng)格生成技術(shù)［28］引入起伏地表下黏彈介質(zhì)的網(wǎng)格剖分，采用牽引力鏡像法實(shí)現(xiàn)起伏地表處的自由邊界條件。由于采用了正交貼體網(wǎng)格，可以直接實(shí)施笛卡爾坐標(biāo)系下的自由邊界條件，從而避免了繁雜的坐標(biāo)旋轉(zhuǎn)或插值運(yùn)算。在曲線坐標(biāo)系中，通常情況下變量并非嚴(yán)格定義在半網(wǎng)格點(diǎn)上［29］，因而必須進(jìn)行復(fù)雜的插值運(yùn)算求取沒(méi)有定義在交錯(cuò)位置上的變量值，這不僅會(huì)降低計(jì)算效率，還極大地影響了差分算法的精度。因此，本文采用Boegy等［30］提出的選擇性濾波同位網(wǎng)格有限差分法計(jì)算曲線坐標(biāo)系下的一階速度—應(yīng)力方程。最后通過(guò)數(shù)值算例驗(yàn)證本文方法的精度和適用性。

1 網(wǎng)格剖分方法

正交貼體網(wǎng)格生成實(shí)際上是一個(gè)數(shù)學(xué)變換過(guò)程，即由任意形狀的物理域(x，z)變換到直角四邊形的計(jì)算域(ξ，η)（圖1）。在這個(gè)變換中點(diǎn)與點(diǎn)之間存在一一對(duì)應(yīng)關(guān)系，在數(shù)學(xué)上可以表示為

圖1 網(wǎng)格映射示意圖

RL 正交貼體網(wǎng)格可由（x，z）平面上的一對(duì)Laplace方程生成［28］

式中：x，ξ≡?x/?ξ表示x對(duì)ξ的一階導(dǎo)數(shù)；f是畸變函數(shù)，由ξ和η方向的刻度因子給出，具體形式為

式中：hξ和hη分別表示ξ和η方向的正交刻度因子；gξξ和gηη是包含坐標(biāo)系信息的度量張量分量，計(jì)算公式為

如果f= 1，式（2）就變成保角變換，則f= 1稱為完全光滑條件。對(duì)于地表起伏不大的模型，通過(guò)式（2）基本可以確保內(nèi)部網(wǎng)格具有較優(yōu)秀的正交性和光滑性。但是，當(dāng)?shù)乇砥鸱容^大時(shí)，為了保證邊界上的正交性，往往會(huì)造成內(nèi)部網(wǎng)格點(diǎn)的嚴(yán)重畸變（圖2）。為此，Zhang等［27］引入ξ、η方向的光滑刻度因子(ξ)i和(η)j，并通過(guò)與正交刻度因子加權(quán)平均得到改進(jìn)的畸變函數(shù)

圖2 內(nèi)部網(wǎng)格點(diǎn)嚴(yán)重畸變的網(wǎng)格剖分示意圖

式中：i和j分別為ξ和η方向的網(wǎng)格點(diǎn)序號(hào)；Nξ和Nη是ξ和η方向的網(wǎng)格點(diǎn)總數(shù)；rξ和rη是經(jīng)驗(yàn)參數(shù)，可以取常數(shù)，也可以根據(jù)以下公式計(jì)算

利用改進(jìn)的RL 算法兼顧了網(wǎng)格剖分的正交性和光滑性，避免了起伏地表附近網(wǎng)格點(diǎn)的畸變（圖3）。

圖3 改進(jìn)后的RL 正交貼體網(wǎng)格剖分示意圖

2 數(shù)值模擬方法

2.1 曲線坐標(biāo)系下的黏彈介質(zhì)波動(dòng)方程

利用鏈?zhǔn)椒▌t和笛卡爾坐標(biāo)系下的Kelvin 黏彈介質(zhì)一階速度—應(yīng)力方程，可以推導(dǎo)出曲線坐標(biāo)系下的二維黏彈介質(zhì)一階速度—應(yīng)力波動(dòng)方程。矩陣形式的波動(dòng)方程可表示為

式中：U為由速度和應(yīng)力組成的列矩陣；為U對(duì)時(shí)間的一階微分；A、B、C、D是系數(shù)矩陣。具體有

其中：vx、vz分別表示速度的水平分量和垂直分量；τxx、τzz、τxz為應(yīng)力分量；ρ為密度；λ、μ、λ'和μ'為拉梅系數(shù)，與縱、橫波速度vP、vS及縱、橫波品質(zhì)因子QP、QS的關(guān)系為

式中ω為圓頻率。

2.2 空間離散算法

對(duì)式（11）的空間導(dǎo)數(shù)，本文用優(yōu)化的7 點(diǎn)同位網(wǎng)格中心差分格式［31］求解。該格式中所有變量和模型介質(zhì)參數(shù)都定義在同一個(gè)網(wǎng)格點(diǎn)上（圖4）。對(duì)于一維問(wèn)題的差分格式為

圖4 優(yōu)化的7 點(diǎn)同位網(wǎng)格有限差分示意圖

式中：Δξ為空間步長(zhǎng)；aj為中心差分系數(shù)。

中心差分格式只能用于內(nèi)部網(wǎng)格點(diǎn)上變量偏導(dǎo)數(shù)的離散求解，對(duì)于模型邊界附近的網(wǎng)格點(diǎn)，由于缺少模型外網(wǎng)格點(diǎn)上的變量值，因此必須使用非中心差分格式求解。本文釆用Berland 等［32］提出的7 點(diǎn)非中心差分格式計(jì)算邊界網(wǎng)格點(diǎn)上的偏導(dǎo)數(shù)

式中：P=1或2，Q=5或4；bj為非中心差分系數(shù)。

同位網(wǎng)格差分格式在求解一階速度—應(yīng)力波動(dòng)方程時(shí)會(huì)產(chǎn)生高頻的格點(diǎn)振蕩現(xiàn)象。為了消除格點(diǎn)振蕩帶來(lái)的影響，對(duì)式（18）和式（19）的計(jì)算結(jié)果做選擇性濾波處理

式中：ud為經(jīng)過(guò)濾波后的波場(chǎng)；F為濾波衰減函數(shù)；σ∈[0，1]是濾波衰減因子；dj為濾波算法的權(quán)系數(shù)。當(dāng)P=Q=3時(shí)，對(duì)中心差分格式計(jì)算結(jié)果進(jìn)行濾波；P=2、Q=4 和P=1、Q=5 時(shí)，對(duì)非中心差分格式計(jì)算結(jié)果進(jìn)行濾波。

2.3 邊界條件的實(shí)施

圖5 是自由界面附近垂直方向?qū)?shù)的離散求解示意圖。設(shè)j=n為自由界面。圖中點(diǎn)A、B為內(nèi)部網(wǎng)格點(diǎn)，采用7 點(diǎn)中心差分格式求解。點(diǎn)C、D為邊界網(wǎng)格點(diǎn)，分別采用非中心差分格式求解。在正交貼體網(wǎng)格下，由于地表處網(wǎng)格線上的法向量與地表正交，因此可以直接采用水平界面時(shí)的自由邊界條件。根據(jù)自由邊界條件，地表網(wǎng)格點(diǎn)上的兩個(gè)應(yīng)力分量為零，即

圖5 自由邊界處有限差分格式示意圖

自由地表以上虛擬網(wǎng)格點(diǎn)上的應(yīng)力則可以根據(jù)應(yīng)力鏡像法求得，即

在自由地表處對(duì)式（11）中的應(yīng)力水平分量進(jìn)行求解時(shí)，需要用到速度分量在η方向上的導(dǎo)數(shù)。將式（21）代入式（11）并整理，可得速度自由邊界條件

利用式（23）可求出vx在η方向的導(dǎo)數(shù)，vz在η方向的導(dǎo)數(shù)可以用相同的方法求取。因此，利用式（21）～式（23）就可以完整的實(shí)施自由邊界條件。

其他邊界采用基于輔助微分方程的完美匹配層（ADE-PML）吸收邊界條件［33］消除人工邊界反射。

3 數(shù)值算例分析

3.1 起伏地表均勻介質(zhì)模型

為了測(cè)試本文方法對(duì)起伏地表模型的有效性和精度，設(shè)計(jì)如圖6a 所示的起伏地表均勻介質(zhì)模型，參數(shù)為：vP=3.3 km/s，vS=2.2 km/s，ρ=2.0 g/cm3，QP=25，QS=20。圖6b 為利用本文方法得到的正交貼體網(wǎng)格，其計(jì)算域中兩個(gè)方向的空間網(wǎng)格步長(zhǎng)都為5 m?？梢钥闯?，本文提出的網(wǎng)格剖分法不僅能夠精確描述起伏地表，而且在邊界上具有精確的正交性。主頻為75 Hz的Ricker子波作為震源，位于（500 m，5 m）處，分別采用本文方法和規(guī)則網(wǎng)格有限差分法進(jìn)行模擬。圖7為利用兩種方法模擬的0.15 s時(shí)刻的波場(chǎng)快照，可以看出，本文方法能夠有效壓制邊界處階梯狀網(wǎng)格引起的虛假散射波。圖8為利用本文方法模擬的直達(dá)波在（500 m，600 m）處的單道地震記錄與黏彈解析解及彈性波解析解的對(duì)比，可以看出，本文方法的數(shù)值解與黏彈解析解完全吻合，且與彈性波解相比振幅明顯衰減。

圖8 本文方法模擬的直達(dá)波與黏彈解析解及彈性解析解的對(duì)比

3.2 起伏地表多層介質(zhì)模型

為了測(cè)試本文方法對(duì)起伏地表模型的適用性，建立了起伏地表兩層和三層黏彈介質(zhì)模型（圖9），各層參數(shù)如表1所示。采用與起伏地表均勻模型相同的方法對(duì)起伏地表兩層和三層模型進(jìn)行網(wǎng)格剖分（圖6b，地形相同，網(wǎng)格剖分相同）。

表1 起伏地表多層介質(zhì)模型參數(shù)

圖9 起伏地表兩層（a）和三層（b）黏彈介質(zhì)模型

以主頻為75 Hz 的Ricker 子波作為震源，位于（500 m，5 m），圖10和圖11分別為兩層模型本文方法模擬的0.20和0.25 s時(shí)刻的波場(chǎng)快照。由圖可見(jiàn)，地震波在傳播到山峰和山谷的拐點(diǎn)時(shí)都會(huì)產(chǎn)生強(qiáng)烈的散射波和繞射波，且發(fā)生面波和體波之間的相互轉(zhuǎn)換，散射面波的能量要強(qiáng)于散射體波；當(dāng)?shù)卣鸩▊鞑サ浇缑鏁r(shí)會(huì)產(chǎn)生反射波、透射波和轉(zhuǎn)換波。圖12為三層模型本文方法模擬的兩分量單炮記錄。從圖中可以觀察到起伏地表拐點(diǎn)處產(chǎn)生的強(qiáng)烈散射波和兩個(gè)介質(zhì)分界面上的反射波、轉(zhuǎn)換波。由于地層界面是崎嶇的，因此反射波同相軸發(fā)生了畸變，不再是以炮點(diǎn)為中心的雙曲線。起伏地表兩層和三層黏彈介質(zhì)模型計(jì)算結(jié)果表明本文方法對(duì)復(fù)雜模型同樣有較強(qiáng)的適用性。

圖10 兩層模型本文方法模擬的0.20 s 時(shí)刻的波場(chǎng)快照

圖11 兩層模型本文方法模擬的0.25 s 時(shí)刻的波場(chǎng)快照

圖12 三層模型本文方法模擬的單炮記錄

4 結(jié)論

地震波數(shù)值模擬是研究地層中地震波傳播規(guī)律的有效方法。本文采用改進(jìn)的RL 正交貼體網(wǎng)格生成方法對(duì)起伏地表下的黏彈性介質(zhì)進(jìn)行網(wǎng)格剖分，并利用優(yōu)化的選擇性濾波同位網(wǎng)格有限差分法模擬曲線坐標(biāo)系下的一階速度—應(yīng)力方程。數(shù)值算例表明：

（1）改進(jìn)的RL 算法能夠?qū)ζ鸱乇砟Ｐ瓦M(jìn)行精確的正交貼體網(wǎng)格剖分，從而有效消除階梯狀網(wǎng)格引起的虛假散射波；

（2）網(wǎng)格的正交性使得在實(shí)施自由邊界條件時(shí)無(wú)需做復(fù)雜的坐標(biāo)轉(zhuǎn)換和插值運(yùn)算，從而有效提高了模擬精度；

（3）因?yàn)槎S和三維使用的是相同的差分格式，所以本文算法很容易向三維推廣。