基于CPML-RML組合邊界條件粘彈TTI介質(zhì)旋轉(zhuǎn)交錯網(wǎng)格有限差分正演模擬

張奎濤,顧漢明,劉少勇,劉春成,陳寶書,張 立,肖逸飛

(1.中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院,地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室,湖北武漢430074;2.中海油研究總院有限責(zé)任公司,北京100028)

實(shí)際地下介質(zhì)中充滿流體、裂縫與裂隙,會同時表現(xiàn)出各向異性與粘彈性特征,對地震波的振幅和相位造成影響,因此,在進(jìn)行地震波場數(shù)值模擬時,這種粘滯性特征不能忽視。粘彈各向異性現(xiàn)象在實(shí)驗(yàn)室與波場測量中得到了證實(shí)[1]。

對于粘彈各向異性介質(zhì),國內(nèi)外學(xué)者進(jìn)行了大量的研究。CARCIONE[2-3]基于廣義標(biāo)準(zhǔn)線性固體(GSLS)模型,建立了線性粘彈性各向異性本構(gòu)關(guān)系及相應(yīng)的波動方程,為后續(xù)在粘彈性各向異性介質(zhì)中的應(yīng)用進(jìn)行了大量的研究[4-7]。受制于實(shí)際計算和存儲能力,地震數(shù)值模擬通常僅能在有限空間進(jìn)行,因此,需要在計算區(qū)域周圍增加人工邊界條件以降低邊界反射干擾。完全匹配層(PML)吸收邊界條件是一種典型的人工邊界條件,由于其穩(wěn)健和高效的吸收特性而被廣泛使用[8-12]。它由BERENGER[13]于1994年提出,之后許多學(xué)者在其基礎(chǔ)上加以改進(jìn),使其變得更加高效。其中最具代表性的是RODEN等[14]提出的卷積完全匹配層(CPML)吸收邊界條件,由于增加了額外的控制參數(shù)(α),其對低頻反射與掠射波有較好的吸收效果。隨后,許多學(xué)者將CPML吸收邊界條件應(yīng)用到一階速度-應(yīng)力方程的數(shù)值模擬中,并取得了較滿意的效果[15-18]。LI等[19]在CPML吸收邊界條件中引入了輔助記憶變量,避免了卷積操作并將其成功應(yīng)用到二階波動方程數(shù)值模擬中。WANG等[20]使用CPML吸收邊界條件對分?jǐn)?shù)時間導(dǎo)數(shù)的粘聲波方程進(jìn)行了數(shù)值模擬。然而,傳統(tǒng)的PML吸收邊界條件在大炮檢距掠射情況下吸收效果不理想,對于某些介質(zhì)存在不穩(wěn)定性[21],而CPML吸收邊界條件在復(fù)雜介質(zhì)構(gòu)造情況下,對低頻反射的吸收效果不佳,干擾了有效波。另一類人工邊界條件是隨機(jī)介質(zhì)層(RML)邊界條件。RML最早由CLAPP[22]提出,隨后有學(xué)者將其應(yīng)用于逆時偏移中[23-24],方法是在邊界區(qū)域設(shè)置隨機(jī)介質(zhì),使得波場傳播到邊界區(qū)域時發(fā)生不相干散射,削弱邊界反射的能量,達(dá)到減少人工邊界反射干擾的目的。

有許多不同的網(wǎng)格被應(yīng)用于有限差分正演模擬中[25-27]。其中最為經(jīng)典并應(yīng)用最廣泛的是標(biāo)準(zhǔn)交錯網(wǎng)格(stand-staggered-grid,SSG)。SSG由MADARIAGA[28]在研究斷層破碎時所提出,隨后被廣泛應(yīng)用于有限差分?jǐn)?shù)值模擬中[29-31]。交錯網(wǎng)格是將彈性參數(shù)定義在整網(wǎng)格點(diǎn)或半網(wǎng)格點(diǎn)上,對于包含多參數(shù)模型的數(shù)值模擬,在計算某些量的空間導(dǎo)數(shù)時,存在所需模型參數(shù)在相應(yīng)空間節(jié)點(diǎn)沒有定義的問題,需要對模型空間進(jìn)行插值或近似處理。這些模型空間的處理會引入計算誤差,在模型參數(shù)變化較為劇烈區(qū)域甚至?xí)霈F(xiàn)數(shù)值計算不穩(wěn)定現(xiàn)象。為了解決SSG在多參數(shù)模型數(shù)值模擬中存在的問題,SAENGER等[32-33]提出了旋轉(zhuǎn)交錯網(wǎng)格(rotated-staggered-grid,RSG)算法,因其將彈性參數(shù)定義在整網(wǎng)格點(diǎn)或網(wǎng)格中心上,在計算空間導(dǎo)數(shù)時無需進(jìn)行插值處理,從而避免了插值處理帶來的計算誤差,提高了數(shù)值模擬的精度。此外,由于RSG是利用網(wǎng)格對角線上的信息來求解空間導(dǎo)數(shù),因此比SSG更加穩(wěn)定,得到的數(shù)值解更加可靠。由于RSG具有較好的穩(wěn)定性,被眾多學(xué)者應(yīng)用于各種復(fù)雜介質(zhì)的數(shù)值模擬中[34-36]。

本文在前人研究的基礎(chǔ)上,采用基于GSLS模型的粘彈TTI介質(zhì)波動方程進(jìn)行數(shù)值模擬。在邊界處理上,利用隨機(jī)邊界的散射特點(diǎn)進(jìn)一步改善CPML吸收邊界條件的吸收能力,提出將CPML吸收邊界條件與RML邊界條件相結(jié)合的策略,并推導(dǎo)出相應(yīng)的旋轉(zhuǎn)交錯網(wǎng)格有限差分格式,利用均勻介質(zhì)模型及復(fù)雜介質(zhì)Hess模型對組合邊界條件的正確性及有效性進(jìn)行了測試及對比分析。

1 理論與方法

1.1 粘彈TTI介質(zhì)波動方程

根據(jù)CARCIONE[2-3]基于GSLS模型建立的線性粘彈性各向異性理論,可得到粘彈TTI介質(zhì)波動方程。

應(yīng)力-應(yīng)變關(guān)系:

(1)

(2)

其中,K為廣義非壓縮模量,G為廣義剛性模量,cij為彈性系數(shù),eij為記憶變量分量。

記憶變量方程為:

(3)

式中:Qv(v=1,2)分別為縱、橫波品質(zhì)因子;f0為地震子波主頻。

(5)

其中,

速度與應(yīng)力關(guān)系可用動量守恒方程表示:

(6)

(7)

1.2 常規(guī)SPML吸收邊界條件

本文考慮二維情況,在復(fù)數(shù)伸展坐標(biāo)下,(6)式可以寫成如下形式[38]:

(8)

(9)

對(9)式求偏導(dǎo),有:

(10)

式中:sx,sz為復(fù)頻伸展函數(shù)。

(11)

將(10)式與(11)式代入(8)式,并在x,z方向分裂有:

將(12)式變換到時間域便可得到SPML吸收邊界條件[18],即

1.3 CPML吸收邊界條件

以(6)式為例(忽略外力),推導(dǎo)出適用于粘彈TTI介質(zhì)的二維CPML吸收邊界條件。對于CPML吸收邊界條件,引入了兩個額外的衰減參數(shù)κ和α,使得新的復(fù)頻伸展函數(shù)變?yōu)?

(14)

式中:κi,αi為輔助衰減系數(shù),且κi和αi滿足κi≥1,αi≥0,i=x,z。

整理(14)式得到:

(15)

將(10)式與(15)式代入(8)式中的第1個方程(以第1個方程為例,后面類似)有:

(16)

令

(17)

則(16)式可寫成:

(18)

整理(17)式得到:

本次共入選受試者236例，其中試驗(yàn)組脫落9例、對照組脫落6例。最終208例患者進(jìn)入符合方案數(shù)據(jù)集（PPS），221例患者進(jìn)入全分析數(shù)據(jù)集（FAS），231例患者進(jìn)入安全性數(shù)據(jù)集（SS）。全部病例均簽署知情同意書。

(19)

將(19)式變換到時間域,有:

(20)

不失一般性地,對于方程

(21)

其通解為:

(22)

(22)式在時間域離散形式下有:

(23)

即有:

(24)

對比(22)式,(18)式的解為:

(25)

式中:

(26)

i=x,z

則(18)式在時間域的形式為:

(27)

式中:Txx,Txz可采用(25)式迭代求解得到。對比(6)式中vx分量項(xiàng)與(27)式可以發(fā)現(xiàn),只需將?i替換成?i/κi+Tij(i,j=x,z)即可得到相應(yīng)的粘彈TTI介質(zhì)的CPML吸收邊界條件。

1.4 CPML-RML組合邊界條件

一般而言,對于單一隨機(jī)邊界條件,由于隨機(jī)邊界區(qū)域填充的為隨機(jī)介質(zhì),可被認(rèn)為由一個個散射點(diǎn)組成,當(dāng)波場傳播到該區(qū)域時,波場發(fā)生隨機(jī)散射形成非相干能量,從而達(dá)到減少人工邊界反射的目的。然而,由于散射方向是任意的,因此,還是有部分較強(qiáng)的能量會被散射到有效物理計算區(qū)域,從而污染了原有的地震波場;另一方面,在某些大角度入射的情況下,CPML吸收邊界條件對低頻反射及掠射波的吸收能力有限,可采用隨機(jī)介質(zhì)將其能量散射。因此,本文利用隨機(jī)介質(zhì)層(RML)邊界的散射特點(diǎn),提出CPML吸收邊界條件與RML邊界條件相組合的策略,進(jìn)一步提升邊界的吸收效果。

圖1為CPML-RML組合邊界示意圖。其中,隨機(jī)邊界區(qū)域填充的隨機(jī)介質(zhì)可采用隨機(jī)介質(zhì)建模得到。在具體實(shí)施過程中,根據(jù)隨機(jī)介質(zhì)建模理論公式,給定自相關(guān)長度(分散程度)、相關(guān)半徑、方向角度和方差等隨機(jī)介質(zhì)建模參數(shù),然后將計算區(qū)域的彈性參數(shù)取均值作為隨機(jī)邊界區(qū)域彈性參數(shù)的背景值,最后采用隨機(jī)介質(zhì)建模程序得到隨機(jī)彈性參數(shù),并將其賦值在隨機(jī)邊界區(qū)域的每個網(wǎng)格點(diǎn)處,其中隨機(jī)邊界區(qū)域的彈性參數(shù)滿足其均值為背景值,方差為給定的方差。因此,當(dāng)波場傳播到邊界時,首先會進(jìn)入CPML邊界區(qū)域進(jìn)行吸收衰減,然后會進(jìn)入隨機(jī)邊界區(qū)域,能量會被散射,使得最終的人工邊界反射被極大地衰減,從而達(dá)到非常好的吸收效果。

圖1 CPML-RML組合邊界示意圖解

1.5 旋轉(zhuǎn)交錯網(wǎng)格有限差分格式

旋轉(zhuǎn)交錯網(wǎng)格(RSG)的基本思想是將彈性參數(shù)定義在整網(wǎng)格點(diǎn)或網(wǎng)格中心點(diǎn)上,如圖2所示。

圖2 旋轉(zhuǎn)交錯網(wǎng)格示意圖解

在網(wǎng)格點(diǎn)①處存放vi、ρ,在網(wǎng)格點(diǎn)②處存放σij、eij、cij,然后在對角線方向上計算相應(yīng)量的空間導(dǎo)數(shù),此過程中避免了插值操作,使得編程實(shí)現(xiàn)更簡單,正演模擬的精度和穩(wěn)定性更高。其計算公式為:

(28a)

(28b)

式中:i,j,n分別為x,z,t方向樣點(diǎn)序號;u為速度波場或應(yīng)力波場;Lx,Lz分別為x,z方向上的微分算子;Δx,Δz分別為x,z方向上的空間離散步長;N為差分階數(shù)的1/2;m為差分階數(shù)序號;cm為差分系數(shù),可采用Taylor展開推導(dǎo)得到。

二維情況下的粘彈TTI介質(zhì)波動方程CPML吸收邊界條件旋轉(zhuǎn)交錯網(wǎng)格有限差分格式為:

(29a)

(29b)

(30a)

(30b)

(30c)

(31a)

(31b)

(31c)

(32)

式中:u為速度波場或應(yīng)力波場;Δt為時間步長。

2 數(shù)值模擬與分析

2.1 二維均勻介質(zhì)模型

為了對比SPML吸收邊界條件、CPML吸收邊界條件與CPML-RML組合邊界條件的吸收效果,本文設(shè)計了一個大小為2500m×500m的均勻介質(zhì)模型,網(wǎng)格間距為5m×5m,時間步長為1ms,總時間采樣點(diǎn)數(shù)為2000,共2s記錄長度,采用30Hz主頻的雷克子波,空間差分階數(shù)為20階,邊界層厚度為250m(對于組合邊界,CPML和RML邊界厚度各125m),震源的位置為(750m,150m),接收排列范圍為0～2500m,排列深度為150m,接收點(diǎn)間距為10m,共251道接收,其均勻介質(zhì)彈性參數(shù)見表1。

圖3、圖4和圖5分別給出了SPML吸收邊界條件、CPML吸收邊界條件和CPML-RML組合邊界條件的x分量地震記錄及波場快照。波場快照從上到下對應(yīng)的時間依次為0.4s、0.6s和1.0s,并且為了區(qū)分對比,地震記錄及波場快照均做了增益處理,使其振幅處于同一水平范圍,即在制圖過程中,將所有數(shù)據(jù)的振幅統(tǒng)一限制在某一更低(相比于原始數(shù)據(jù)范圍)的數(shù)值范圍內(nèi)(例如將所有波場快照數(shù)據(jù)的振幅限制在-0.2～0.2),這樣既能達(dá)到增益效果,又能實(shí)現(xiàn)振幅對比。從圖3、圖4和圖5中可以看出,在0.6s之前,不論采用哪種邊界條件,其吸收效果均很好,未看見明顯的邊界反射,說明3種吸收邊界條件在波場傳播的前期均有好的吸收效果。但值得注意的是,當(dāng)波場傳播到后期(1.0s),即波場大角度入射到邊界時,SPML吸收邊界條件會出現(xiàn)如紅色箭頭所示的邊界反射,即SPML吸收邊界條件不能很好地將大角度的邊界反射吸收;而CPML吸收邊界條件會出現(xiàn)弱邊界反射,但其效果優(yōu)于SPML吸收邊界條件;CPML-RML組合邊界條件從地震記錄及波場快照上均未出現(xiàn)邊界反射。說明與SPML吸收邊界條件和CPML吸收邊界條件相比,本文提出的CPML-RML組合邊界條件在吸收人工邊界反射方面更加優(yōu)越。

表1 粘彈TTI均勻介質(zhì)彈性參數(shù)

圖3 SPML吸收邊界條件的x分量地震記錄(a)與波場快照(b)

圖4 CPML吸收邊界條件的x分量地震記錄(a)與波場快照(b)

圖5 CPML-RML組合邊界條件的x分量地震記錄(a)與波場快照(b)

圖6對比了3種邊界條件第211道的地震記錄,圖6b為圖6a部分區(qū)域的放大結(jié)果。從圖6可以看出,SPML吸收邊界條件出現(xiàn)了較明顯的邊界反射,CPML吸收邊界條件則出現(xiàn)了弱邊界反射,而CPML-RML組合邊界條件未出現(xiàn)邊界反射。

不同邊界條件的計算效率如表2所示。表2中的計算效率是在Win10系統(tǒng)、i5CPU(2.7GHz)、4G內(nèi)存以及VS2013+IVF2013編程環(huán)境下得出的。從表2中可以明顯地看到,粘彈TTI介質(zhì)情形下,與SPML吸收邊界條件相比,CPML-RML組合邊界的計算效率高出18.2%,這是由于SPML吸收邊界條件需要將方程進(jìn)行分裂而導(dǎo)致的,而CPML吸收邊界條件與CPML-RML組合邊界條件計算效率相當(dāng);同樣地,由于粘彈性TTI介質(zhì)方程比彈性TTI介質(zhì)方程多了記憶變量方程,因此,在相同的邊界條件下,彈性TTI介質(zhì)方程的計算效率比粘彈性TTI介質(zhì)方程高49.1%,但粘彈性TTI介質(zhì)方程更能夠反映地下真實(shí)介質(zhì)情況;另一方面,各向同性介質(zhì)的計算效率比彈性TTI介質(zhì)的計算效率高6.3%。

圖6 3種邊界條件第211道地震記錄對比a 顯示時間范圍為0～2.0s; b 顯示時間范圍為0.865～2.000s

因此,隨著介質(zhì)復(fù)雜程度的提高,計算效率隨之降低,但更加接近地下真實(shí)介質(zhì)的情況。此外,由于CPML吸收邊界條件不需要對方程進(jìn)行分裂,因而其計算效率比SPML吸收邊界條件高,實(shí)現(xiàn)簡單,編程方便。由于CPML-RML組合邊界條件的吸收效果優(yōu)于CPML吸收邊界條件,因此更加適用于粘彈性介質(zhì)方程的數(shù)值模擬。

表2 不同邊界條件的計算效率

2.2 二維復(fù)雜Hess模型

為了說明CPML-RML組合邊界條件對復(fù)雜介質(zhì)的適應(yīng)性及正確性,采用二維Hess粘彈TTI介質(zhì)模型進(jìn)行測試。二維Hess模型大小為3000m×2500m,網(wǎng)格間距為5m×5m,時間步長為0.2ms,總時間采樣點(diǎn)數(shù)為12500,2.5s記錄長度,采用30Hz主頻的雷克子波,空間差分階數(shù)為20階,邊界層厚度為300m(對于組合邊界,CPML和RML邊界厚度各150m),震源位置為(2000m,10m),接收排列范圍為0～3000m,排列深度為0,接收點(diǎn)間距為10m,共301道接收。

圖8給出了彈性TTI介質(zhì)與粘彈TTI介質(zhì)采用組合邊界條件在0.6s與0.8s時刻的波場快照。波場快照上未出現(xiàn)邊界反射,說明本文提出的組合邊界策略在復(fù)雜構(gòu)造介質(zhì)中的正確性及有效性;同時,也可以從波場快照中明顯地看見彈性TTI介質(zhì)與粘彈TTI介質(zhì)在能量及旅行時上的差異,表明粘彈介質(zhì)存在振幅衰減及相位延遲。圖9給出了彈性TTI介質(zhì)與粘彈性TTI介質(zhì)采用組合邊界條件得到的地震記錄。由圖9也可以看出,彈性TTI介質(zhì)與粘彈性TTI介質(zhì)在能量及相位延遲上的差異,同時,在粘彈TTI介質(zhì)地震記錄中,淺部能量強(qiáng)、深部能量弱,與實(shí)際采集得到的地震記錄特征相同,而彈性TTI介質(zhì)未表現(xiàn)出該特征,說明采用粘彈TTI介質(zhì)進(jìn)行數(shù)值模擬更加符合實(shí)際地質(zhì)情況,為研究地下波場特征提供了更多的理論依據(jù)。

圖7 二維Hess粘彈TTI介質(zhì)(θ=0,φ=60°)模型參數(shù)a 縱波速度vP; b 橫波速度vS; c 縱波各向異性參數(shù)ε; d 橫波各向異性參數(shù)γ; e 變異系數(shù)δ; f 密度ρ; g 縱波品質(zhì)因子Q1; h 橫波品質(zhì)因子Q2

圖8 彈性TTI介質(zhì)與粘彈性TTI介質(zhì)采用組合邊界條件不同時刻波場快照a 彈性TTI介質(zhì)0.6s時刻波場快照; b 粘彈TTI介質(zhì)0.6s時刻波場快照; c 彈性TTI介質(zhì)0.8s時刻波場快照; d 粘彈TTI介質(zhì)0.8s時刻波場快照

圖9 彈性TTI介質(zhì)(a)與粘彈性TTI介質(zhì)(b)采用組合邊界條件得到的地震記錄

圖10為粘彈TTI介質(zhì)分別采用SPML吸收邊界條件、CPML吸收邊界條件與CPML-RML組合邊界條件得到的地震記錄。對比圖10a、圖10b和圖10c 可以看出,采用SPML吸收邊界條件得到的地震記錄上出現(xiàn)了邊界反射,采用CPML吸收邊界條件得到的地震記錄上出現(xiàn)了低頻反射,而采用CPML-RML組合邊界條件未出現(xiàn)邊界反射,說明本文提出的組合邊界條件吸收效果好。

圖10 粘彈性TTI介質(zhì)采用不同邊界條件得到的地震記錄a SPML吸收邊界條件; b CPML吸收邊界條件; c CPML-RML組合邊界條件

3 結(jié)束語

本文在前人對人工邊界條件研究的基礎(chǔ)上,提出了CPML吸收邊界條件與RML邊界條件相結(jié)合的策略,并將其應(yīng)用于基于GSLS模型的粘彈TTI介質(zhì)數(shù)值模擬中,同時給出了適用于粘彈TTI介質(zhì)的CPML吸收邊界條件旋轉(zhuǎn)交錯網(wǎng)格有限差分算法,以提高數(shù)值模擬的穩(wěn)定性及數(shù)值解的精確性。采用二維均勻介質(zhì)模型和復(fù)雜Hess模型測試了本文組合邊界條件的適應(yīng)性及正確性,驗(yàn)證了其吸收效果,得到以下結(jié)論:

1) 與SPML吸收邊界條件、CPML吸收邊界條件相比,本文提出的組合邊界條件在不增加計算量的情況下具有更好的吸收效果,能有效減少人工邊界反射、提高數(shù)值模擬的精度,獲得更好的數(shù)值模擬結(jié)果,適合于粘彈TTI介質(zhì)的數(shù)值模擬;

2) 與TTI介質(zhì)相比,粘彈TTI介質(zhì)中的地震波場表現(xiàn)出明顯的振幅衰減和相位延遲,地震記錄中淺層能量強(qiáng)、深層能量弱,較好地表現(xiàn)出了地下真實(shí)介質(zhì)的粘滯性特征。說明本文組合邊界條件適用于復(fù)雜地下介質(zhì)。