基于Marmousi模型的CFS-NPML分?jǐn)?shù)階粘聲波正演模擬

廖萬平， 王興建， 李卿武， 王崇名， 張俊杰

(1.成都理工大學(xué) 地球物理學(xué)院，成都 610059;2.油氣藏地質(zhì)及開發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室(成都理工大學(xué)),成都 610059)

0 引言

地震波場(chǎng)正演模擬是預(yù)先構(gòu)建數(shù)學(xué)模型或?qū)嶓w模型，然后通過研究地震波在模型中的傳播規(guī)律與實(shí)際地震勘探資料進(jìn)行對(duì)比，通過不斷修正模型盡可能使模擬數(shù)據(jù)達(dá)到與探測(cè)資料接近，以至于能夠更精確更高效地推斷地下地層情況。但是在數(shù)值模擬常會(huì)遇到一個(gè)問題，在有界區(qū)域內(nèi)模擬聲波的傳播，如果強(qiáng)行加截?cái)噙吔纾瑫?huì)導(dǎo)致虛假反射現(xiàn)象[1]，因此為了在有限區(qū)域內(nèi)達(dá)到聲波正演模擬的真實(shí)效果，需要一種反射非常小甚至不產(chǎn)生反射的人工吸收邊界條件[2]。

目前已有許多方法，用以消除邊界反射。Cerjan等[3]通過在計(jì)算機(jī)邊界引入損耗介質(zhì)來衰減向外傳播的波，雖然取得了一定吸收效果，但是其對(duì)大角度邊界反射吸收方面卻不盡人意；Clayton[4]吸收邊界在旁軸近似理論的基礎(chǔ)上導(dǎo)出，但是其在特定的入射角和頻率范圍內(nèi)，才具有較好的吸收效果；Liao Z P[5]基于牛頓后向差分公式推廣了多次透射理論，應(yīng)用平面波原理及二次插值方法的到，其實(shí)現(xiàn)簡(jiǎn)單，且精度能夠很好的滿足大多數(shù)工程計(jì)算，然而其吸收邊界條件的不穩(wěn)定性限制了其應(yīng)用；Bérenger[6]提出了一種相對(duì)于傳統(tǒng)吸收邊界更高效的PML邊界條件，其能夠吸收來自各個(gè)方向、各種頻率的波，而不產(chǎn)生任何反射，但隨著對(duì)邊界的研究和應(yīng)用的不斷深入，學(xué)者們發(fā)現(xiàn)當(dāng)入射波傳播方向于PML區(qū)域和內(nèi)部區(qū)域接近平行時(shí)會(huì)導(dǎo)致其無法對(duì)高角度入射波/掠射波以及瞬逝波完全吸收。隨著簡(jiǎn)單高效的CFS-NPML技術(shù)提出，其不僅解決了分裂場(chǎng)引起的數(shù)值誤差，還提高了對(duì)高角度入射波/掠射波和瞬逝波的吸收效果，使得邊界條件達(dá)到滿足實(shí)際生產(chǎn)的要求[7]。CFS-NPML吸收邊界的優(yōu)越性其已被廣泛應(yīng)用于電磁波、地震波等波的傳播數(shù)值模擬問題中，尤其是在地震波正演模擬領(lǐng)域，均獲得了豐碩的研究成果。但是到目前為止鮮有基于CFS-NPML吸收邊界在分?jǐn)?shù)階粘聲波方面的應(yīng)用。

這里以粘聲波動(dòng)方程的速度-應(yīng)力方程為基礎(chǔ)，完成了對(duì)C-E吸收邊界、廖氏吸收邊界、PML吸收邊界、CFS-NPML吸收邊界在均勻介質(zhì)下的對(duì)比，然后更進(jìn)一步將PML吸收邊界和CFS-NPML吸收邊界應(yīng)用到分?jǐn)?shù)階粘聲[9]條件下，并采用較為復(fù)雜的Marmousi模型進(jìn)行數(shù)值模擬。結(jié)果表明與傳統(tǒng)方法PML邊界條件相比，CFS-NPML吸收邊界具有更好的魯棒性、可行性以及吸收效果，并且能夠直觀高效地反映聲波在粘性介質(zhì)中的傳播衰減規(guī)律。

1 方法與原理

1.1 分?jǐn)?shù)階粘聲波動(dòng)方程

分?jǐn)?shù)階粘聲波動(dòng)方程描述的是聲波在粘性介質(zhì)中的傳播規(guī)律由Kjartansson[9]提出。聲波在二維均勻介質(zhì)中傳播時(shí)的質(zhì)點(diǎn)運(yùn)動(dòng)平衡方程為式(1)。

(1)

式中：p為聲壓；ux和uz分別為質(zhì)點(diǎn)在x、z軸方向的位移分量。

根據(jù)Caputo[10]分?jǐn)?shù)階導(dǎo)數(shù)定義,線性黏彈性介質(zhì)中本構(gòu)關(guān)系的褶積形式可轉(zhuǎn)化為乘積形式，即

(2)

結(jié)合式(1)、式(2)和應(yīng)變-位移關(guān)系(幾何方程)可得式(3)。

(3)

根據(jù)式(1)、式(2)可推導(dǎo)出常Q衰減介質(zhì)的聲波速度-應(yīng)力方程為式(4)。

(4)

式中：vx、vz表示質(zhì)點(diǎn)在x、z軸方向的振動(dòng)速度；D1-2γ為時(shí)間1-2γ的階導(dǎo)數(shù)。

聲波相速度頻散關(guān)系和衰減系數(shù)分別為：

(5)

(6)

式中：cp為相速度；α為衰減系數(shù)。可以看出，兩者均與頻率有關(guān)：前者表明常Q衰減介質(zhì)中波的傳播速度是頻率的函數(shù)，即存在頻散現(xiàn)象；后者說明介質(zhì)對(duì)地震波能量吸收有頻率選擇性。

1.2 完全匹配層(PML)

完全吸收邊界條件思想就是在需要計(jì)算場(chǎng)值區(qū)域內(nèi)采用速度-應(yīng)力方程，在計(jì)算區(qū)域外使用一定厚度的截?cái)噙吔纾?dāng)波場(chǎng)運(yùn)動(dòng)到邊界時(shí)，不會(huì)發(fā)生反射，而是穿過邊界，在邊界區(qū)域內(nèi)場(chǎng)值按指數(shù)方式衰減，達(dá)到消弱邊界反射的效果。在PML吸收層區(qū)域內(nèi)，對(duì)頻率空間域的x和z方向偏導(dǎo)分別作如下替換：

(7)

式中：ω為角頻率；dx和dz分別表示為和方向的阻尼因子。

在PML邊界條件應(yīng)用時(shí)選擇合適的衰減系數(shù)極為重要，選擇不當(dāng)會(huì)導(dǎo)致嚴(yán)重的邊界反射現(xiàn)象。因此為達(dá)到降低反射現(xiàn)象的最佳效果，筆者采用Collino等[12]提出的一種PML邊界與內(nèi)部計(jì)算區(qū)域截?cái)噙吔缇嚯x為指數(shù)關(guān)系的一種衰減函數(shù)：

(8)

式中：L為PML邊界層的厚度；r為離PML層內(nèi)邊界的距離；Vmax為最大縱波度；Δγ為空間步長(zhǎng)；R為理論反射系數(shù)。

1.3 非分裂復(fù)頻移完全匹配層(CFS-NPML)

PML邊界條件經(jīng)過一般性推導(dǎo)將時(shí)間域變換導(dǎo)頻率域，通過引入復(fù)伸展坐標(biāo)。在頻率域中復(fù)伸展坐標(biāo)可以表示為式(9)與式(10)。

(9)

(10)

對(duì)式(9)、式(10)式引入復(fù)頻移伸展(CFS)函數(shù)。其形式為式(11)。

(11)

式中：γ和α為復(fù)頻移伸展函數(shù)中尺度因子和頻移因子，滿足γ≥1和α≥0。當(dāng)這兩參數(shù)同時(shí)取等號(hào)時(shí)就變?yōu)槌Ｒ?guī)PML邊界條件。參數(shù)γ主要對(duì)廣角入射波/掠射波以及瞬逝波進(jìn)行吸收。α主要對(duì)波的低頻成分進(jìn)行吸收，避免衰減系數(shù)出現(xiàn)奇異值。本文參數(shù)表達(dá)式如下

dθ(θ)=dmax(r/L)Pd

(12)

γθ(θ)=1+(Gmax-1)(r/L)Pk

(13)

αθ(θ)=Amax[1-(r/L)Pα]

(14)

(15)

式中：L為邊界網(wǎng)格厚度；Pd、Pk、Pα為多項(xiàng)式衰減因子階數(shù)，一般是在[1,4]之間，這里分別取2、3、1；Δγ為空間步長(zhǎng)；vmax為最大縱波速度；r為層內(nèi)點(diǎn)距邊界與非邊界的距離；Gmax通常介于1～20，這里取1；Amax=1.0*pi*f0，f0為震源主頻；R為邊界吸收系數(shù)，這里取0.000 1%；

1.4 吸收邊界厚度優(yōu)選分析

從圖1給出的關(guān)系可以看出衰減項(xiàng)dθ(θ)、最大衰減項(xiàng)dmax的大小與吸收邊界厚度成反比。由于邊界的存在導(dǎo)致地震波傳到邊界時(shí)會(huì)產(chǎn)生反射波，而反射強(qiáng)度的大小取決于衰減項(xiàng)和最大衰減項(xiàng)的值[13]。當(dāng)CFS-NPML層數(shù)減小時(shí)，衰減項(xiàng)、最大衰減項(xiàng)增大，導(dǎo)致數(shù)值反射強(qiáng)度增強(qiáng)，從而影響CFS-NPML吸收效果收衰減效果。由圖1可以看出，當(dāng)邊界層厚度在20層到30層時(shí)曲線變化越來越小，越來越趨于平緩，綜合邊界條件的實(shí)用性和經(jīng)濟(jì)性考慮，即達(dá)到好的吸收衰減特征，同時(shí)減小數(shù)值反射波，CFS-NPML吸收邊界選取20層～30層較為理想。

圖1 吸收層厚度與衰減項(xiàng)、最大衰減項(xiàng)關(guān)系Fig.1 Relation between absorption layer thickness and attenuation term and maximum attenuation term(a)吸收層厚度與衰減項(xiàng)關(guān)系;(b)吸收層厚度與最大衰減項(xiàng)關(guān)系

2 數(shù)值算例

2.1 均勻介質(zhì)模型

為了對(duì)比驗(yàn)證本文的CFS-NPML邊界條件吸收效果和穩(wěn)定性，圖2、圖3分別建立了C-E吸收邊界、廖氏吸收邊界、PML吸收邊界、CFS-NPML吸收邊界在均勻介質(zhì)條件下的數(shù)值模擬和第200道波場(chǎng)記錄。模型離散網(wǎng)格大小為400×300，空間網(wǎng)格采樣步長(zhǎng)X=Z=5 m，時(shí)間采樣步長(zhǎng)為0.000 1 s，震源位于(200,150)，震源采用30 Hz的雷克子波。由圖2、圖3可以明顯看出C-E吸收邊界能吸收大部分到達(dá)邊界的聲波，但還是存在小部分虛假反射。廖氏吸收邊界能夠較好吸收到達(dá)邊界的聲波，但仍然存在虛假反射現(xiàn)象。由于圖3(c)和圖3(d)采用更高效的PML吸收邊界條件，因此其穩(wěn)定性能和吸收效果均好于傳統(tǒng)C-E吸收邊界、廖氏吸收邊界均未出現(xiàn)明顯虛假反射現(xiàn)象。與PML邊界條件相比CFS-NPML吸收邊界在PML邊界基礎(chǔ)上加入了能夠吸收廣角入射波/掠射波以及瞬逝波和低頻成分的吸收參數(shù)γ和α,能夠有效吸收和壓制反射波。從結(jié)果可以看出，經(jīng)典PML邊界到CFS-NPML邊界吸收效果和正演模擬精度在不斷提升。

圖2 不同吸收邊界在175 ms時(shí)刻波場(chǎng)快照Fig.2 Snapshot of the wave field at 175 ms at different absorption boundaries(a)C-E吸收邊界;(b)廖氏吸收邊界;(c)PML吸收邊界;(d)CFS-NPML吸收邊界

圖3 不同吸收邊界在175 ms時(shí)刻第200道記錄Fig.3 Different absorption boundaries were recorded at the 200th track at 175 ms(a)C-E吸收邊界波場(chǎng)記錄;(b)廖氏吸收邊界波場(chǎng)記錄;(c)PML吸收邊界波場(chǎng)記錄;(d)CFS-NPML吸收邊界波場(chǎng)記錄

2.2 Marmousi正演模型

為了進(jìn)一步對(duì)比PML吸收邊界和CFS-NPML吸收邊界的吸收效果，同時(shí)驗(yàn)證CFS-NPML邊界吸收條件在分?jǐn)?shù)階粘聲波正演模擬中的魯棒性、可行性，筆者采用具有復(fù)雜地質(zhì)構(gòu)造和較強(qiáng)橫向速度變化，且被廣泛使用的標(biāo)準(zhǔn)Marmousi模型。參數(shù)模型見圖4，模型是一個(gè)計(jì)算大小約為6.6 km×2.4 km的狹長(zhǎng)區(qū)域，離散網(wǎng)格點(diǎn)為692×263，網(wǎng)格間距為dx=dx=10 m,震源采用零相位的雷克子波主頻30 Hz，可使分辨率達(dá)到極限，采樣點(diǎn)時(shí)間間隔為1 ms，震源位于(250,1)處，結(jié)合吸收效果及計(jì)算成本問題，作者選擇邊界上20層作為吸收區(qū)域，分別采用PML吸收邊界和CFS-NPML吸收邊界進(jìn)行分?jǐn)?shù)階粘聲波數(shù)值模擬，得到的地震記錄如圖5所示。

圖4 Marmousi的衰減和速度模型Fig.4 Marmousi's attenuation and velocity model(a)衰減模型;(b)速度模型

圖5 PML邊界和CFS-NPML邊界下的Marmousi地震記錄Fig.5 Marmousi seismic records under PML boundary and CFS-NPML boundary(a)PML邊界;(b)CFS-NPML邊界

根據(jù)圖5可以看出，在粘滯性介質(zhì)中地震波場(chǎng)能量由淺到深，逐漸變?nèi)酰顚臃瓷湔穹鶞p弱,分辨率逐漸降低。在深度約3 km分界線以下PML吸收邊界條件下的地震記錄消失，而CFS-NPML吸收邊界條件下的地震記錄仍能分辨。但當(dāng)?shù)竭_(dá)4 km深度時(shí)，傳統(tǒng)PML邊界條件下的數(shù)值模擬出現(xiàn)不穩(wěn)定現(xiàn)象，而在CFS-NPML邊界條件下的地震記錄未出現(xiàn)任何明顯變化，由此可看出CFS-NPML吸收邊界能夠有效地壓制虛假反射，阻斷了反射波淹沒地震記錄提高了信噪比。圖6為在不同邊界條件下放大的振幅和能量變化圖，可以明顯看出，使用CFS-NPML邊界條件進(jìn)行正演模擬的振幅大小略強(qiáng)于PML邊界正演和能量變化相比與PML邊界條件穩(wěn)定，從而證明CFS-NPML邊界穩(wěn)定性和吸收效果上強(qiáng)于PML邊界條件。

圖6 不同邊界振幅變化放大圖和能量衰減變化放大圖Fig.6 Amplification diagram of amplitude changes at different boundaries variation diagram of energy attenuation(a)振幅變化放大圖;(b)能量衰減變化放大圖

為進(jìn)一步說明CFS-NPML吸收邊界的吸收效果及穩(wěn)定性、可行性，給出了PML和CFS-PML邊界條件下的不同時(shí)刻波場(chǎng)快照。對(duì)比圖7(a)和圖7(b)可以看出當(dāng)t=800 ms時(shí)，實(shí)現(xiàn)的兩種吸收邊界均能很好地反應(yīng)聲波在粘滯性介質(zhì)中的傳播規(guī)律，并且在邊界處，沒有看到明顯的虛假反射現(xiàn)象，可以驗(yàn)證兩種吸收邊界均能夠有效吸收到達(dá)邊界的聲波。通過圖8可以看出，當(dāng)t=2 800 ms時(shí)，PML與CFS-NPML邊界吸收條件在邊界條件內(nèi)發(fā)生微妙變化，PML邊界處的波場(chǎng)開始模糊看不清，而CFS-NPML邊界處的能夠保持其穩(wěn)定的波場(chǎng)記錄。根據(jù)圖9可以明顯看出，當(dāng)t=4 800 ms時(shí)，CFS-NPML邊界條件內(nèi)的波場(chǎng)值頻率大于PML邊界內(nèi)的頻率，且保持其粘性衰減特性，這與前面主要對(duì)波的低頻成分進(jìn)行吸收一致。

圖7 t=800 ms處 PML邊界和CFS-NPML邊界的波場(chǎng)快照Fig.7 t=800 ms snapshot of the wave field of the T-PML boundary and the CFS-NPML boundary(a)PML邊界;(b)CFS-NPML邊界

圖8 t=2 800 ms處 PML邊界和CFS-NPML邊界的波場(chǎng)快照Fig.8 t=2 800 ms snapshot of the wave field of the T-PML boundary and the CFS-NPML boundary(a)PML邊界;(b)CFS-NPML邊界

圖9 t=4 800 ms處 PML邊界和CFS-NPML邊界的波場(chǎng)快照Fig.9 t=4 800 ms snapshot of the wave field of the T-PML boundary and the CFS-NPML boundary(a)PML邊界;(b)CFS-NPML邊界

復(fù)頻移PML雖然在計(jì)算時(shí)加入了更多的輔助變量，但是與經(jīng)典PML的分裂式相比，采用了非分裂形式，在求解運(yùn)算時(shí)減少了大量的卷積運(yùn)算。從表1也可以看出，PML在黏聲正演模擬中所需時(shí)間比CFS-NPML所需時(shí)間多，所占內(nèi)存多，可以看出在同樣條件下CFS-NPML在計(jì)算效率和存儲(chǔ)空間上的優(yōu)勢(shì)。

表1 不同邊界條件黏聲正演模擬計(jì)算時(shí)間和存儲(chǔ)量

3 結(jié)論

筆者通過對(duì)四種吸收邊界條件對(duì)虛假反射的吸收效果進(jìn)行對(duì)比，然后更進(jìn)一步地進(jìn)行對(duì)基于Marmousi模型的CFS-NPML分?jǐn)?shù)階粘聲正演模擬研究，將PML、CFS-NPML邊界條件應(yīng)用到接近實(shí)際地層條件的分?jǐn)?shù)階粘聲正演模擬中，并采用Marmousi模型來測(cè)試邊界條件在復(fù)雜地層下的穩(wěn)定性、可行性。通過對(duì)PML、CFS-NPML兩種邊界條件在數(shù)值模擬中的效果對(duì)比可得到如下結(jié)論：

與C-E吸收邊界、廖氏吸收邊界相比PML和CFS-NPML邊界條件能夠明顯地壓制波場(chǎng)的虛假反射，并保持其穩(wěn)定性；與傳統(tǒng)PML邊界條件相比，CFS-NPML邊界能夠在分?jǐn)?shù)階粘聲正演模擬中保持其魯棒性、可行性，并且能夠有效吸收高角度入射波/掠射波以及瞬逝波，避免發(fā)生虛假反射現(xiàn)象而淹沒地震記錄，提高了數(shù)值模擬信噪比；在CFS-NPML邊界條件內(nèi)的分?jǐn)?shù)階粘聲正演模擬能夠很好地反映聲波因地下介質(zhì)的粘滯性，而產(chǎn)生對(duì)波的吸收和衰減等特性。有助于我們更好地認(rèn)識(shí)地震波在粘性介質(zhì)中的傳播規(guī)律，對(duì)補(bǔ)償因地下介質(zhì)粘滯性因素影響而損失的地震波能量具有深遠(yuǎn)意義。