段沛然, 谷丙洛, 李振春, 李青陽

1 中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院, 青島 266580

2 深層油氣重點(diǎn)實(shí)驗(yàn)室, 青島 266580

3 中國石油長慶油田分公司勘探開發(fā)研究院, 西安 710018

0 引言

目前,油氣地震勘探的重點(diǎn)已由老油田、老區(qū)塊的簡單構(gòu)造逐漸向山地、丘陵等復(fù)雜地表區(qū)域轉(zhuǎn)移(劉紅偉等,2011;岳玉波等,2012).這些區(qū)域勘探面臨地表起伏大、近地表橫向速度變化劇烈等難題,給地震資料偏移成像帶來巨大的挑戰(zhàn)(裴正林,2004;黃建平等,2014;楊宏偉等,2022).對于復(fù)雜起伏地表地震資料的偏移成像方法主要有兩類：一是在偏移前采用高程基準(zhǔn)面靜校正,將地震資料校正到浮動或者固定基準(zhǔn)面上,然后再實(shí)施偏移成像(Berryhill,1979;Yilmaz and Lucas,1986;李振春等,2007);二是基于起伏地表直接進(jìn)行偏移成像(Wiggins,1984;Beasley and Lynn,1992;Margrave and Yao,2000;何英等,2002;Shragge and Sava,2005;程玖兵等,2006;Liu et al.,2011).當(dāng)基巖出露,地表高程變化大,近地表橫向速度變化劇烈時,靜校正往往會扭曲波場,降低地震成像的質(zhì)量.從起伏地表直接進(jìn)行偏移成像的方法,其算法設(shè)計隱含靜校正過程,具有自然適應(yīng)復(fù)雜地表條件的能力,一直是是業(yè)內(nèi)研究的熱點(diǎn)和目標(biāo).

由起伏地表直接偏移的成像方法包括射線類偏移、單程波方程偏移和雙程波方程偏移.與射線類(Beylkin,1985;Gray and Bleistein,2009)和單程波方程(Claerbout,1971;Mulder and Plessix,2004)偏移方法相比,基于雙程波方程的逆時偏移不受傾角和偏移孔徑限制,避免了上下行波分離,可對多種波場(一次反射波、多次波、棱柱波等)準(zhǔn)確成像,因此在復(fù)雜地下構(gòu)造成像方面具有巨大的理論優(yōu)勢(Whitmore,1983;Baysal et al.,1983;Gu et al.,2015).波場延拓作為逆時偏移算法的核心環(huán)節(jié),其關(guān)鍵在于高精度的地震波數(shù)值模擬方法.目前,起伏地表地震波數(shù)值模擬方法研究主要分為網(wǎng)格劃分和自由地表邊界條件兩個部分.其中,網(wǎng)格劃分可分為規(guī)則網(wǎng)格、曲變網(wǎng)格、非結(jié)構(gòu)網(wǎng)格以及無網(wǎng)格.基于規(guī)則網(wǎng)格的有限差分方法(Finite Difference Method,FDM)(Kelly et al.,1976;Marfurt,1984;Virieux,1986)具有計算效率高、易于實(shí)現(xiàn)、適用于任意復(fù)雜介質(zhì)等優(yōu)點(diǎn),是目前應(yīng)用最為廣泛的地震波數(shù)值模擬方法(Ren et al.,2017;吳國忱等,2018; 吳建魯和吳國忱,2018; 段沛然等,2020).為了解決FDM中的起伏地表問題,Jih等(1988)利用具有不同角度的線段擬合地表起伏,采用單邊近似法處理不規(guī)則自由表面條件實(shí)現(xiàn)起伏地表地震波數(shù)值模擬.該方法需采用精細(xì)網(wǎng)格剖分實(shí)現(xiàn)線段貼合起伏地表,因此計算效率較低.Hayashi等(2001)在FDM中引入局部變網(wǎng)格剖分模擬起伏地表自由表面的瑞雷波.變網(wǎng)格由于依賴于網(wǎng)格形變程度,僅能減少雜散衍射,無法完全消除矩形網(wǎng)格階梯近似所引起的衍射(孫成禹等,2008; 李振春等,2014).Tessmer等(1992)將變網(wǎng)格思想引入交錯網(wǎng)格FDM中,但其僅有二階精度(Hestholm and Ruud,2002),難以應(yīng)用于實(shí)際.Zhang和Chen(2006)根據(jù)實(shí)際地表起伏劃分曲線網(wǎng)格,實(shí)現(xiàn)起伏地表貼體網(wǎng)格FDM.de la Puente等(2014)采用基于完全交錯網(wǎng)格的垂向變換網(wǎng)格和模擬地震波建模的網(wǎng)格變形策略,實(shí)現(xiàn)了起伏地表條件下的彈性波數(shù)值模擬.Shragge 和Tapley(2017)采用廣義結(jié)構(gòu)網(wǎng)格法求解張量聲波方程.Chen 等(2022)將浸入邊界法與平面波透射邊界相結(jié)合,提出一種適用于起伏地表的浸入吸收邊界,該方法基于常規(guī)的矩形網(wǎng)格剖分,簡單高效.貼體網(wǎng)格是一類曲變網(wǎng)格,適用于同位網(wǎng)格和交錯網(wǎng)格,使用范圍廣且較易生成.但是貼體網(wǎng)格要求地形起伏的函數(shù)處處可導(dǎo),因此不適用于地表起伏劇烈的區(qū)域(如崎嶇地表、近垂直斷崖等地質(zhì)條件).非結(jié)構(gòu)網(wǎng)格地震波數(shù)值模擬的實(shí)現(xiàn)方法主要有限元法(K?ser and Iske,2005;董興鵬和楊頂輝,2017)、邊界元法(符力耘和牟永光,1994; 何彥鋒等,2013)、譜元法(Komatitsch et al.,2000;Cristini and Komatitsch,2012)、間斷伽遼金法(Cockburn et al.,2000;K?ser and Dumbser,2006)、格子法(徐義和張劍鋒,2008;孫輝等,2019)等.這些方法均基于非結(jié)構(gòu)網(wǎng)格或非規(guī)則網(wǎng)格算子的全空間離散,能夠很好的貼合起伏地表,但是這些方法大多采用全局積分方程弱形式,計算量大,長時間模擬存在不穩(wěn)定問題.

無網(wǎng)格法是近些年來興起的一類數(shù)值計算新方法,它是有限元等傳統(tǒng)方法的重要補(bǔ)充和發(fā)展.基于無網(wǎng)格的徑向基函數(shù)有限差分方法(Radius-Basis-Function Finite Difference Method,RBF-FDM)能夠有效避免虛假反射,具有節(jié)點(diǎn)剖分隨模型速度變化,復(fù)雜地質(zhì)體剖分靈活等優(yōu)點(diǎn)(劉立彬等,2020).RBF是RBF-FDM中的一類徑向?qū)ΨQ函數(shù),用于計算差分權(quán)重,并求解不同類型的偏微分方程(Tolstykh and Shirobokov,2003; Fornberg et al.,2013; Fornberg and Flyer,2015;Dehghan and Mohammadi,2019).在地震領(lǐng)域,Martin 等(2013, 2015)首次將RBF-FDM應(yīng)用于聲波方程,實(shí)現(xiàn)了二維非均勻介質(zhì)的數(shù)值模擬.Takekawa等(2014a,b)研究了RBF-FDM中的自由表面條件,并將哈密頓粒子法引入RBF-FDM中,進(jìn)行了起伏地表聲波數(shù)值模擬.Takekawa和Mikada(2018)在RBF-FDM中提出了一種適用于無網(wǎng)格的自由表面邊界條件,并將其應(yīng)用于頻域彈性波數(shù)值模擬.Li等(2017)進(jìn)一步將含時-空域優(yōu)化算子的RBF-FDM應(yīng)用于聲波數(shù)值模擬,同時采用混合吸收邊界條件壓制邊界反射.Takekawa等(2016,2018)通過增加相鄰節(jié)點(diǎn)數(shù)目來提高數(shù)值模擬精度,將多參數(shù)泰勒展開應(yīng)用于非均勻介質(zhì)RBF-FDM.Liu等(2019)利用RBF-FDM進(jìn)行頻率域聲波方程數(shù)值模擬,采用完全匹配層(Perfectly Matched Layers,PML)邊界條件來壓制邊界反射.Duan等(2021)提出了一種基于頻散關(guān)系和穩(wěn)定性條件約束節(jié)點(diǎn)間隔的自適應(yīng)節(jié)點(diǎn)剖分方法,實(shí)現(xiàn)了靈活穩(wěn)定的RBF-FDM聲波數(shù)值模擬.Wu等(2021,2022)實(shí)現(xiàn)了基于RBF-FDM的聲波逆時偏移和全波形反演.

除了起伏地表地震波數(shù)值模擬外,國內(nèi)外已有多位學(xué)者對起伏地表標(biāo)量波場逆時偏移成像方法進(jìn)行了研究(張曉波等,2022;孫志廣等,2023;姚宗惠等,2023).徐義(2008)提出了一種適用于起伏地表的格子法疊前逆時度偏移方法.Zhang 等(2010)提出了一種基于地表起伏非反射遞歸算法的疊后逆時偏移算法.劉紅偉等(2011)給出了一種實(shí)現(xiàn)起伏地表疊前逆時偏移的方法,該方法通過簡化自由邊界條件,結(jié)合GPU加速算法實(shí)現(xiàn)了高效穩(wěn)定的起伏地表疊前逆時偏移.Lan等(2014)基于貼體網(wǎng)格地形“平化”策略較好地擬合平滑地表起伏,實(shí)現(xiàn)了貼體網(wǎng)格的起伏地表逆時偏移.Liu等(2017)提出了一種適用于起伏地表的非結(jié)構(gòu)網(wǎng)格法疊前逆時偏移.然而,實(shí)際地震波是包含縱橫波的矢量場,將縱橫波的耦合特性綜合應(yīng)用于基于矢量波場的彈性波逆時偏移,彌補(bǔ)了標(biāo)量波成像在巖性識別和構(gòu)造解釋等方面存在多解性的缺陷.曲英銘等(2015)提出一種基于曲坐標(biāo)系的分層坐標(biāo)變換法,實(shí)現(xiàn)起伏地表條件下的彈性波疊前逆時偏移.Qu等(2021)將此方法推廣到具有起伏地表的棱柱波彈性波逆時偏移方法中.Naveed和Chen(2019)利用曲網(wǎng)格進(jìn)行波場延拓,利用貼體網(wǎng)格處理起伏地表實(shí)現(xiàn)了彈性波逆時偏移.Zhong等(2021)基于Zhang和McMechan(2010)提出的波場分離方法,利用地形起伏相關(guān)的濾波器實(shí)現(xiàn)了起伏地表彈性波逆時偏移.上述起伏地表彈性波逆時偏移方法均以結(jié)構(gòu)網(wǎng)格剖分為基礎(chǔ),計算量較大.面對當(dāng)前復(fù)雜的勘探目標(biāo),進(jìn)一步研究起伏地表彈性波逆時偏移方法仍具有重要意義.

綜上所述,基于RBF-FDM的地震波數(shù)值模擬方法研究主要局限于水平地表聲學(xué)介質(zhì),而針對起伏地表彈性介質(zhì)的RBF-FDM數(shù)值模擬的理論和方法研究尚不完善.另外,針對起伏地表彈性波逆時偏移,目前主要采用曲坐標(biāo)系下的貼體網(wǎng)格等方法進(jìn)行波場延拓,盡管曲網(wǎng)格能較好的貼合地表起伏,但對于起伏地表的自由邊界條件適用性不足.本文在聲波RBF-FDM數(shù)值模擬基礎(chǔ)上,提出一種優(yōu)化的起伏地表無網(wǎng)格彈性波RTM方法.通過定量分析基于RBF-FDM的彈性波數(shù)值模擬的頻散關(guān)系與穩(wěn)定性條件,本文發(fā)展了高精度的QR徑向基函數(shù)有限差分方法(QR-RBF-FDM),同時提出一種優(yōu)化的起伏地表自適應(yīng)節(jié)點(diǎn)剖分技術(shù),形成了新的基于無網(wǎng)格的起伏地表彈性波數(shù)值模擬方法.由于起伏地表無網(wǎng)格邊界條件加載困難,本文提出了適用于無網(wǎng)格的精確自由邊界條件和彈性波混合吸收邊界條件.此外,本文將此無網(wǎng)格RBF-FD應(yīng)用于精確的縱橫波場矢量分解公式(Gu et al.,2019),實(shí)現(xiàn)了起伏地表彈性波逆時偏移成像.通過對高斯山丘模型,起伏凹陷模型和起伏地表Marmousi-2模型的數(shù)值試算驗(yàn)證了本文提出方法在起伏地表彈性波數(shù)值模擬與逆時偏移成像中的可行性與有效性.

1 徑向基函數(shù)有限差分方法彈性波數(shù)值模擬

1.1 徑向基函數(shù)彈性波有限差分方法

各向同性介質(zhì)中,彈性波二階位移方程可表示為

(1)

其中,u(x,z,t)和w(x,z,t)表示在空間x、z方向上的位移波場分量,ρ(x,z)表示介質(zhì)密度,λ(x,z)和μ(x,z)為拉梅常數(shù).在起伏地表條件下,傳統(tǒng)的FDM很難滿足要求.基于無網(wǎng)格節(jié)點(diǎn)離散模型的RBF-FDM能夠有效避免虛假反射,具有剖分節(jié)點(diǎn)隨模型參數(shù)變化,起伏界面處理靈活等優(yōu)勢(劉立彬等,2020).本文將RBF-FDM引入彈性波波動方程數(shù)值模擬中,以便獲得高精度的彈性波地震波場.

RBF-FDM是基于一組局部無序的離散節(jié)點(diǎn),根據(jù)與節(jié)點(diǎn)距離相關(guān)的徑向基函數(shù)構(gòu)造差分算子矩陣實(shí)現(xiàn)偏微分方程.在任一時刻,若給定Ω區(qū)域內(nèi)存在波場U(x)及無序離散節(jié)點(diǎn)x={x1,x2,…,xn},節(jié)點(diǎn)位置及參數(shù)由點(diǎn)生成方法可知,波場U(x)可由徑向基函數(shù)φ(r)近似表示為

(2)

式中,r節(jié)點(diǎn)間距(中心點(diǎn)與其最近鄰近節(jié)點(diǎn)的距離),rj=‖x-xj‖表示鄰近節(jié)點(diǎn)xj距離中心節(jié)點(diǎn)x的距離,‖·‖表示歐幾里得范數(shù);N為中心節(jié)點(diǎn)x差分模板所需的鄰近節(jié)點(diǎn)個數(shù),根據(jù)鄰近當(dāng)前任意點(diǎn)的距離選取差分模板的范圍,選擇N個最近點(diǎn),記為差分模板,差分模板中的點(diǎn)按距離遠(yuǎn)近排序;{c1,c2,…,cN}為待定的RBF-FD系數(shù).根據(jù)無網(wǎng)格RBF-FD基本原理,對于任意線性偏微分算子L,中心點(diǎn)節(jié)點(diǎn)x0處波場U(x)|x=x0的任意空間偏導(dǎo)數(shù)可以表示為

(3)

(4)

其中,徑向基函數(shù)系數(shù)矩陣的條件數(shù)是影響RBF-FD差分系數(shù)求解的重要因素,為保證差分方程的穩(wěn)定,條件數(shù)的取值范圍應(yīng)在106～108(Flyer et al.,2014;Martin et al.,2015;Li et al.,2017).實(shí)際中,該條件數(shù)由當(dāng)前差分模板的節(jié)點(diǎn)分布,差分格式,節(jié)點(diǎn)間隔和徑向基函數(shù)類型及形變參數(shù)等相關(guān).根據(jù)條件數(shù)經(jīng)驗(yàn)取值范圍(106～108),在類均勻分布(如圖1)的模型中得到了一組節(jié)點(diǎn)數(shù)與節(jié)點(diǎn)間隔和形變參數(shù)的最佳取值范圍(見表2,參見Li et al., 2017).

表1 徑向基函數(shù)有限差分方法中常用的徑向基函數(shù)Table 1 Commonly used radial basis function in RBF-FDM

采用RBF-FDM求解彈性波方程,首先通過時間二階差分來代替式(1)中的時間微分算子如下：

表2 類均勻分布情況下RBF-FD最優(yōu)差分節(jié)點(diǎn)數(shù)、節(jié)點(diǎn)間隔和形變參數(shù)的選取Table 2 Optimal ranges of parameters for RBF-FD with different numbers of neighbors in quasi-uniform distribution

(5)

其中,上標(biāo)l為當(dāng)前時刻,τ為時間步長.采用RBF-FDM差分算子對式(1)中u分量的空間導(dǎo)數(shù)進(jìn)行離散：

(6)

(7)

本文采用RBF-FD中常用的K-D算法(KNN-Search)尋點(diǎn),能夠不占用計算量找到計算域中任一中心節(jié)點(diǎn)差分模板中的N個鄰近差分節(jié)點(diǎn).

1.2 數(shù)值頻散分析與穩(wěn)定性條件

在平面波條件下假設(shè)下,彈性波場可以表示為

(8)

(9)

(10)

(11)

對式(11)中兩式相乘得：

+β2(A+B)]=0,

(12)

根據(jù)式(12)可得：

(13)

(14)

將式(10)中A,B代入式(11)化簡可得RBF-FD頻散關(guān)系為

(15)

其中：

(16)

其中,θ表示相角.

穩(wěn)定性條件也是衡量微分方程數(shù)值求解算法有效性的重要指標(biāo).本文根據(jù)傅里葉譜分析法(Von Neumann,1950)推導(dǎo)RBF-FD情況下的穩(wěn)定性條件.此處,將式(1)表示為矩陣形式：

(17)

為了方便起見,進(jìn)一步將式(17)簡寫為

(18)

式中U(t)=[u(t),w(t)]T為多分量波場矢量,Q為二階矩陣.

應(yīng)用RBF-FD的思想求解彈性波方程,對式(18)中的時間導(dǎo)數(shù)進(jìn)行二階有限差分離散可得：

+O(Δt2M),

(19)

將式(18)代入式(19)得：

(20)

設(shè)G(kx,kz)是Q(kx,kz)的傅里葉變換對,對式(20)進(jìn)行空間傅里葉變換可得：

U(t+Δt,kx,kz)-2U(t,kx,kz)+U(t-Δt,kx,kz)

(21)

將U(t,kx,kz)的系數(shù)矩陣記為

(22)

則式(7)改寫為

Un+1=2Un-Un-1+AUn,

(23)

式中n表示當(dāng)前時刻,根據(jù)傅里葉分析方法,令：

(24)

聯(lián)立式(23)、(24)得：

(25)

式中,E為單位矩陣,B為增長矩陣.

根據(jù)Von Neumann穩(wěn)定性條件可知,當(dāng)增長矩陣特征值的絕對值小于或等于1時,遞歸是穩(wěn)定的,即時間穩(wěn)定性條件是增長矩陣B的譜半徑不大于1.由式(25)可知,B的特征值β應(yīng)滿足：

|βE-B|=|β2E-βA-2E|=0,

(26)

由式(26)可知,β的取值與矩陣A有關(guān),因此先求取的A特征值γ：

(27)

根據(jù)式(22),若G的特征值為η,求取A特征值γ需要首先求取G的特征值η.根據(jù)矩陣Q的表達(dá)形式,η應(yīng)滿足：

(28)

(29)

根據(jù)行列式規(guī)則式(28)可化簡為

+Dkzz)2=0,

(30)

利用2次多項式的求解方法計算得到η的2個特征值為

(31)

由式(27)可知,矩陣A特征值γ為

(32)

由于A為2階方陣,且有2個互異的特征值,因此矩陣A可對角化,有：

A=PΛP-1,

(33)

式中P為可逆矩陣,Λ為A的特征值構(gòu)成的對角陣,具體表示為

Λ=diag{γ1,γ2},

(34)

由式(27)可得：

|β2E-βΛ-2E|=0,

(35)

由式(35)可得B的特征值β和A的特征值γ的關(guān)系式為

β2-βγ-2=0,

(36)

即：

(37)

(38)

若系數(shù)矩陣A的特征值γ的絕對值小于或等于1,則滿足穩(wěn)定性條件.對比|γ1|和|γ2|表達(dá)式可知,RBF-FD差分格式需滿足的穩(wěn)定性條件是：

(39)

2 QR徑向基函數(shù)有限差分方法

在常規(guī)RBF-FDM中,形變參數(shù)ε趨近于0時,徑向基函數(shù)趨于平坦而擬合精度越高.然而平坦的徑向基函數(shù)會導(dǎo)致差分矩陣的病態(tài),影響波場延拓的穩(wěn)定性.Flyer等(2012)研究了含泰勒多項式的RBF-FDM,提高了數(shù)值模擬的穩(wěn)定性,但該方法依賴模型參數(shù)的取值,在處理復(fù)雜構(gòu)造模型時,構(gòu)造泰勒多項式所增加的參數(shù)將會導(dǎo)致差分矩陣的病態(tài),使得波場傳播不穩(wěn)定.為了提高模擬精度且不落入病態(tài)問題,Fornberg等(2013)提出一種無需形變參數(shù)的高斯QR基函數(shù),將其成功用于插值問題.本文將高斯QR基函數(shù)的思想引入RBF-FD中,在笛卡爾坐標(biāo)系下構(gòu)造差分矩陣,將高斯徑向基函數(shù)轉(zhuǎn)換至極坐標(biāo)系,并展開為切比雪夫多項式的形式,利用QR分解法求解由切比雪夫多項式構(gòu)成的差分矩陣,有效解決了徑向基函數(shù)的形變參數(shù)取值和矩陣病態(tài)問題,在保證精度的同時增強(qiáng)波場傳播的穩(wěn)定性.

QR徑向基函數(shù)法利用切比雪夫多項式T(x)對高斯徑向基函數(shù)φ(r)=e-i(εr)2級數(shù)展開(Fornberg et al., 2011),并將其變換至極坐標(biāo)系(r,θ)得到：

(40)

其中切比雪夫多項式表示為

(41)

其中,j≥0,0≤m≤(j-p)/2;當(dāng)p=0時,j為偶數(shù),當(dāng)p=1時,j為奇數(shù);公式(40)中尺度因子dj,m為

(42)

當(dāng)形變參數(shù)ε適當(dāng)取值時,尺度因子dj,m的取值范圍為[0,1];公式(40)中參數(shù)cj,m和sj,m為

(43)

其中參數(shù)b0=1,bm=2,m>0;參數(shù)t0=1/2,tj=1,j>0;F2為超越方程組,其中參數(shù)αj,m=(j-2m+p+1)/2,βj,m=[j-2m+1,(j-2m+p+1)/2];由此,公式(40)可改寫為矩陣形式：

(44)

其中,矩陣C由參數(shù)cj,m和sj,m構(gòu)成,三角陣D由尺度因子dj,m構(gòu)成.

利用QR分解對矩陣C進(jìn)行分解代入公式(44)得：

φ(x)=Q·R·D·T(x),

(45)

將公式(45)代入公式(40)得：

(46)

(47)

求解公式(47)可得到優(yōu)化的差分系數(shù).由于切比雪夫多項式的定義域?yàn)閞∈[-1,1],計算差分系數(shù)時需將當(dāng)前點(diǎn)的差分模板預(yù)先縮放到[-1,1]的范圍內(nèi).

圖1 類均勻節(jié)點(diǎn)剖分示意圖

這里采用一個均勻介質(zhì)模型研究QR徑向基函數(shù)有限差分方法(QR-RBF-FDM)的擬合精度、穩(wěn)定性及其頻譜分析.介質(zhì)縱波速度為3000 m·s-1, 橫波速度為1732 m·s-1,模型大小為1000 m×1000 m,節(jié)點(diǎn)間隔為10 m,時間步長為0.5 ms,最大記錄時間為1.0 s.采用20 Hz雷克子波作為縱波震源,位于近(500 m, 20 m)處.圖1給出了類均勻節(jié)點(diǎn)分布圖,紅點(diǎn)和藍(lán)點(diǎn)分別表示任意中心差分點(diǎn)及其對應(yīng)的無網(wǎng)格差分模板包含的N個差分節(jié)點(diǎn).圖2給出不同方法得到的點(diǎn)(500 m, 0 m)處的波形圖及誤差曲線圖.其中,解析解(analytical_solution)由Cagniard-De Hoop法(De Hoop, 1960)得到,FDM表示時間2階空間10階精度的有限差分.由圖2a單道記錄與圖2b振幅誤差曲線圖可見：QR-RBF-FDM 和 FDM的波形和振幅與解析解基本吻合,而RBF-FDM的模擬結(jié)果取決于形變參數(shù)取值.當(dāng)形變參數(shù)為0.01時,RBF-FDM的波形和振幅與解析解基本吻合,而當(dāng)形變參數(shù)為0.1或0.001時,RBF-FDM的模擬結(jié)果與解析解偏差很大.圖2c為單道記錄對應(yīng)的頻譜圖,QR-RBF-FDM能夠有效抑制形變參數(shù)導(dǎo)致的差分矩陣病態(tài)問題,獲得與FDM和解析解相近的頻譜;常規(guī)的RBF-FDM采用不適合的形變參數(shù)將引起波形震蕩,時間空間頻散分別出現(xiàn)在波前的淺部和尾部,同時引起了頻譜畸變,對頻譜影響較大.由此可知,QR-RBF-FDM比常規(guī)的RBF-FDM具有更高的模擬精度,且形變參數(shù)的取值不影響QR-RBF-FDM的結(jié)果.圖3為QR-RBF-FDM、常規(guī)RBF-FDM與FDM的計算效率對比.與FDM相比,QR-RBF-FDM與RBF-FDM的計算成本主要來源于差分系數(shù)計算和差分方程迭代兩個方面.其中,差分方程迭代過程與FDM基本一致,而盡管二者僅需在整個計算過程中計算一次差分系數(shù),但由于每個節(jié)點(diǎn)進(jìn)行差分系數(shù)矩陣求解仍然需要消耗巨大的計算成本,因此QR-RBF-FDM與RBF-FDM在計算成本和計算時間方面為FDM的4倍左右(如圖3).QR-RBF-FDM求解差分系數(shù)過程計算相對常規(guī)的RBF-FDM更加復(fù)雜,因此計算時間較長,但考慮到無需選擇形變參數(shù),因此可避免因形變參數(shù)選擇不當(dāng)而產(chǎn)生的計算冗余.

圖2 均勻模型(500 m, 0 m)處不同方法(a) 單道記錄波形; (b) 振幅誤差曲線; (c) 對應(yīng)的頻譜.

圖3 均勻模型不同方法歸一化計算時間對比圖

3 起伏地表自適應(yīng)節(jié)點(diǎn)剖分方法

節(jié)點(diǎn)剖分是RBF-FDM的第一步,也是連接速度模型與數(shù)值模擬算法的關(guān)鍵.Duan等(2021)提出的自適應(yīng)節(jié)點(diǎn)剖分法利用頻散關(guān)系和穩(wěn)定性條件約束節(jié)點(diǎn)間隔的取值范圍,能夠在矩形計算域內(nèi)生成適應(yīng)模型速度的非均勻節(jié)點(diǎn)分布.但考慮起伏地表的不規(guī)則邊界與地下介質(zhì)的不均勻性,Duan等(2021)的方法將不再適用.因此,在Duan等(2021)的基礎(chǔ)上,本文提出一種自適應(yīng)起伏地表節(jié)點(diǎn)剖分法,邊界節(jié)點(diǎn)貼體起伏地表,并在內(nèi)部節(jié)點(diǎn)剖分過程中,對起伏地表和地下構(gòu)造變換劇烈的區(qū)域依據(jù)速度變化的梯度大小進(jìn)行自適應(yīng)調(diào)節(jié).起伏地表節(jié)點(diǎn)剖分方法(見圖4)采用如下的步驟：

(2)利用自適應(yīng)節(jié)點(diǎn)剖分法進(jìn)行內(nèi)部節(jié)點(diǎn)剖分,剖分方式如圖4bi—iii所示.(i) 設(shè)定最下層邊界節(jié)點(diǎn)為潛在點(diǎn)位.(ii) 根據(jù)速度模型參數(shù)確定各起始點(diǎn)位的節(jié)點(diǎn)間距,設(shè)置以起始點(diǎn)位為中心,半徑為節(jié)點(diǎn)間距的圓域.(iii) 刪除此圓域內(nèi)所有其他的潛在點(diǎn)位,在刪除節(jié)點(diǎn)的圓周范圍內(nèi)放置等距的新點(diǎn)位.由下層潛在點(diǎn)位迭代直至與邊界節(jié)點(diǎn)距離小于節(jié)點(diǎn)間距的一半d/2.記錄內(nèi)部節(jié)點(diǎn)的初始坐標(biāo)(ξ,η)對應(yīng)的速度模型參數(shù),如圖4b(iv)所示.

(3)刪除與邊界節(jié)點(diǎn)間距小于d/2的內(nèi)部節(jié)點(diǎn),如圖4c所示.

(4)根據(jù)自適應(yīng)節(jié)點(diǎn)控制方程(式(51)),將內(nèi)部節(jié)點(diǎn)由初始坐標(biāo)調(diào)整到新的坐標(biāo)(x,z),通常,需要四次以上的迭代獲得滿足精度要求的節(jié)點(diǎn)分布,得到每個節(jié)點(diǎn)的坐標(biāo)及對應(yīng)的速度模型參數(shù),如圖4d所示.

圖4 起伏地表自適應(yīng)節(jié)點(diǎn)剖分過程示意圖

為了計算起伏地表上節(jié)點(diǎn)的法線方向,步驟1中地表節(jié)點(diǎn)的坐標(biāo)變換到極坐標(biāo)系：

(48)

式中(x,z)表示笛卡爾坐標(biāo)系中地表節(jié)點(diǎn)的位置向量,(r,θ)表示極坐標(biāo)系下地表節(jié)點(diǎn)的位置向量,取速度模型左上頂點(diǎn)作為極點(diǎn),則關(guān)于極角θ的梯度為

(49)

由此可知,法向量的分量nx和nz可表示為

(50)

其中sgn(·) 表示常用符號函數(shù),g表示梯度.在步驟4中,引入自適應(yīng)網(wǎng)格理論(Saltzman and Brackbill, 1982)中的控制方程對自適應(yīng)節(jié)點(diǎn)剖分法生成的無網(wǎng)格節(jié)點(diǎn)位置進(jìn)一步調(diào)整,控制方程為

(51)

(52)

J=xξzη-xηzξ,

(53)

其中J表示雅可比多項式,其中一階導(dǎo)數(shù)?x/?ξ,?x/?η,?z/?ξ和?z/?η簡寫為xξ,xη,zξ和zη.

無網(wǎng)格與基于坐標(biāo)變換法的貼體網(wǎng)格和曲網(wǎng)格各有優(yōu)劣.(1)貼體網(wǎng)格和曲網(wǎng)格是起伏地表條件下的地震波模擬中的重要方法,其網(wǎng)格易于生成且在曲坐標(biāo)系下可直接使用水平地表假設(shè)的自由邊界條件.與這兩種網(wǎng)格相比,無網(wǎng)格的優(yōu)勢在于無需坐標(biāo)變換,直接采用彈性波、聲波等方程直接進(jìn)行差分,后期研究中更適用于復(fù)雜的黏性介質(zhì)或各向異性方程.(2)貼體網(wǎng)格和曲網(wǎng)格在生成過程中,僅能通過一些密度函數(shù)進(jìn)行疏密調(diào)節(jié),筆者認(rèn)為此種調(diào)節(jié)疏密通常是不可控的,因?yàn)榫W(wǎng)格數(shù)不變的情況下調(diào)節(jié)一個部分疏密必將使得另一部分密疏,一定程度上將導(dǎo)致數(shù)值頻散以及穩(wěn)定性問題.與這兩種方法相比,無網(wǎng)格以模型速度、密度等參數(shù)為依據(jù),直接生成滿足模型剖分條件的節(jié)點(diǎn)分布,不存在調(diào)節(jié)后一部分點(diǎn)過密一部分點(diǎn)稀疏的現(xiàn)象,因此相對更加具有數(shù)值穩(wěn)定性.(3)貼體網(wǎng)格和曲網(wǎng)格本質(zhì)上仍然是固定的網(wǎng)格間隔,這就造成在深層、中深層進(jìn)行數(shù)值模擬,網(wǎng)格過密造成冗余.通過文獻(xiàn)調(diào)研可知,貼體網(wǎng)格往往需要與變網(wǎng)格相互結(jié)合以提高其計算效率.與之相比,無網(wǎng)格節(jié)點(diǎn)在剖分過程中根據(jù)速度、密度等參數(shù)可自適應(yīng)調(diào)節(jié)節(jié)點(diǎn)間距,因此在深層、中深層等復(fù)雜構(gòu)造中可由細(xì)網(wǎng)格自然過渡到粗網(wǎng)格.

4 與地形相關(guān)的邊界條件

4.1 無網(wǎng)格自由邊界條件

RBF-FDM的優(yōu)勢在于處理起伏自由邊界.在自由表面是固體地球與空氣層之間的物性突變面,常規(guī)自由邊界條件通常假設(shè)垂直自由表面的應(yīng)力為0(Lan and Zhang, 2011),二維水平自由表面(z=0)處的表達(dá)式為

(54)

地形起伏會明顯影響地震波場的傳播,此時水平自由邊界條件的假設(shè)將不適用.為了在RBF-FDM中精確處理起伏地表,需要將無網(wǎng)格節(jié)點(diǎn)自然放置在起伏的曲面地形上,利用起伏地表自適應(yīng)節(jié)點(diǎn)剖分方法獲得起伏界面的法線分量,在此基礎(chǔ)上重新推導(dǎo)起伏界面的精確自由邊界條件.圖5顯示了規(guī)則網(wǎng)格和無網(wǎng)格的剖分示意圖.

圖5 網(wǎng)格剖分示意圖(a) 規(guī)則網(wǎng)格; (b) 無網(wǎng)格節(jié)點(diǎn).紅色曲線表示起伏地表, 藍(lán)色節(jié)點(diǎn)表示內(nèi)部節(jié)點(diǎn),紅色節(jié)點(diǎn)表示地表節(jié)點(diǎn),箭頭表示地表節(jié)點(diǎn)的法向量方向.

因此,本文采用Dirichlet邊界條件推導(dǎo)出無網(wǎng)格節(jié)點(diǎn)的自由曲面邊界條件,將其應(yīng)用于貼合實(shí)際地形的RBF-FDM具有更高的精度.Dirichlet邊界條件如下：

(55)

其中σxx和σzz為正應(yīng)力,σxz為切應(yīng)力.在此基礎(chǔ)上,通過自適應(yīng)節(jié)點(diǎn)生成方法求出了與起伏地表垂直的外法向x分量nx和z分量nz.因此,應(yīng)力分量可以表示為

(56)

將式(56)代入式(55),經(jīng)過化簡可以得到起伏地表自由邊界條件：

(57)

圖6 差分模板示意圖紅點(diǎn)為中心點(diǎn),藍(lán)點(diǎn)為中心點(diǎn)對應(yīng)的差分點(diǎn),橢圓區(qū)域?yàn)檫吔绻?jié)點(diǎn)差分模板,圓形區(qū)域?yàn)閮?nèi)部節(jié)點(diǎn)差分模板.

其中,地表界面上的法向量的分量nx和nz可式(56)計算得到.在規(guī)則網(wǎng)格FDM中,通常需在自由表面上方引入虛擬點(diǎn)以滿足自由表面處邊界條件的差分格式.與之相比,在RBF-FDM中,計算節(jié)點(diǎn)可以放在任意位置進(jìn)行差分,邊界節(jié)點(diǎn)完全貼體起伏地表分布.因此無需特殊處理,在邊界節(jié)點(diǎn)處采用如圖6橢圓區(qū)域所示差分模板,利用QR-RBF-FDM獲得邊界節(jié)點(diǎn)處的空間導(dǎo)數(shù)并運(yùn)用于式(57)實(shí)現(xiàn)無網(wǎng)格的自由邊界條件.其中一階空間導(dǎo)數(shù)可由RBF-FDM近似得到：

(58)

其中x0為地表節(jié)點(diǎn)的坐標(biāo),mj和nj表示由QR-RBF-FDM計算得到的一階差分系數(shù),同理可得垂直分量w的一階空間導(dǎo)數(shù)的表示形式.本文由Dirichlet邊界條件出發(fā),推導(dǎo)了基于RBF-FD的自由邊界條件,相對于與應(yīng)力鏡像法等常規(guī)方法而言,該方法無需設(shè)置虛擬層,在笛卡爾坐標(biāo)系下直接求解邊界條件,處理方式簡單且計算精度高.本文在數(shù)值算例中將對此進(jìn)行驗(yàn)證.

4.2 無網(wǎng)格彈性波混合吸收邊界條件

利用RBF-FDM進(jìn)行起伏地表彈性波逆時偏移時,必須考慮起伏地表情況下的邊界吸收問題.基于Li 等(2017)提出的分段光滑曲線邊界節(jié)點(diǎn)生成方法,本文在現(xiàn)有邊界節(jié)點(diǎn)基礎(chǔ)上,基于彈性波Clayton和Engquist(CE)單程波方程推導(dǎo)了無網(wǎng)格起伏地表條件下的混合吸收邊界條件.起伏地表條件下的彈性波CE單程波方程表示為(任志明,2016)：

(59)

(60)

將式(59)寫為位移分量方程的形式可得：

(61)

通過引入邊界區(qū)域每個節(jié)點(diǎn)對應(yīng)的最近內(nèi)點(diǎn),為單程波方程提供穩(wěn)定的差分：

(62)

其中,位移分量u和w的上標(biāo)-1,0,1表示前一時刻,當(dāng)前時刻與后一時刻,邊界區(qū)域的節(jié)點(diǎn)(表示為下標(biāo)B)與其對應(yīng)的最近內(nèi)點(diǎn)(表示為下標(biāo)I)構(gòu)成的差分項Al和Bl(l=1,2,3,4)為

(63)

(64)

(65)

在此基礎(chǔ)上,在邊界區(qū)域的同時采用雙程波方程與上述單程波方程,即可構(gòu)成起伏地表條件下的無網(wǎng)格混合吸收邊界.混合吸收邊界條件實(shí)現(xiàn)步驟(Liu and Sen, 2010)如下：

(3)對于邊界區(qū)域,將雙程波波場和單程波波場加權(quán)平均得到最終的波場：

(66)

其中,加權(quán)系數(shù)為

aBi=(i-1)/N.

(67)

5 基于徑向基函數(shù)有限差分方法的彈性波逆時偏移

在無網(wǎng)格模型參數(shù)化的基礎(chǔ)上,我們開展了基于RBF-FDM的彈性波逆時偏移成像方法應(yīng)用研究.逆時偏移的理論基礎(chǔ)是Claerbout(1971)提出的時間一致性原理,即反射面存在于地層內(nèi)下行波波至?xí)r間與某一上行波波至?xí)r間相一致之處.因此QR-RBF-FDM彈性波逆時偏移成像可分為三個部分：一是基于QR-RBF-FDM震源波場的正向延拓,二是基于QR-RBF-FDM檢波點(diǎn)波場的逆向延拓,三是應(yīng)用成像條件獲得偏移成像結(jié)果.波場分離是彈性波RTM至關(guān)重要的環(huán)節(jié).Gu等(2019)解決了亥姆霍茲分解導(dǎo)致的波場振幅和相位畸變,建立了精確的波場分解和重新組合方程,并實(shí)現(xiàn)了3D標(biāo)量波彈性波逆時偏移,該方法能夠直接獲得具有正確振幅、相位和物理單位的偏移結(jié)果.本文基于Gu等(2019)的方法,我們推導(dǎo)得到在無網(wǎng)格模型參數(shù)化條件下通過RBF-FDM實(shí)現(xiàn)的精確縱橫波波場分離方程,在正向和逆向波場延拓過程中,利用無網(wǎng)格震源波場和檢波點(diǎn)波場進(jìn)行精確的縱橫波波場分離,并應(yīng)用互相關(guān)成像條件得到成像結(jié)果.

根據(jù)Gu等(2019)的方法推導(dǎo),精確的波場分離方程如下：

(68)

(69)

其中,vp和vs為P-和S-速度.由于無網(wǎng)格法的本質(zhì)是局部無序的離散節(jié)點(diǎn),無法直接求解式(68)、(69)中的積分.考慮到矢量場V僅有子波項與時間有關(guān),因此在正向和反向波場計算過程中,對震源子波做二次積分提取式(68)、(69)中的積分,矢量場中P和S勢則已包含積分,即可去掉式(68)、(69)中的積分.將式(68)代入式(69)可得：

(70)

根據(jù)散度、旋度和梯度的關(guān)系求解式(68)得：

(71)

采用QR-RBF-FDM對式(71)進(jìn)行離散,代入式(6)二階空間算子可得：

(72)

(73)

其中x=(x,y,z)為空間向量;其中sm和rn分別表示源波場的m波模式和接收器波場的n波模式.

6 數(shù)值算例

6.1 高斯起伏模型

本節(jié)通過兩類高斯起伏地表模型(凸型模型和凹型模型)驗(yàn)證在自由邊界條件情況下本文方法的有效性,如圖7a、b所示.模型大小為2 km×1.5 km,其中凸型高斯函數(shù)為

(74)

凹型高斯函數(shù)為

(75)

我們分別對兩類模型采用自適應(yīng)起伏地表節(jié)點(diǎn)剖分法和規(guī)則網(wǎng)格剖分法進(jìn)行離散化.圖8展示了自適應(yīng)起伏地表節(jié)點(diǎn)剖分法離散的結(jié)果,由于本測試選用均勻速遞模型以測試起伏自由地表下QR-RBF-FDM的可靠性和穩(wěn)定性,在速度層上為真空,地表采用自由地表條件,矩形邊界采用無網(wǎng)格彈性波混合吸收條件.模型參數(shù)設(shè)置如下：縱波速度為3000 m·s-1,橫波速度2000 m·s-1,而對應(yīng)的自適應(yīng)起伏地表節(jié)點(diǎn)剖分結(jié)果呈現(xiàn)出節(jié)點(diǎn)間隔相同的特征,其間隔大小為10 m(事實(shí)上當(dāng)速度非均勻條件下,自適應(yīng)起伏地表節(jié)點(diǎn)剖分結(jié)果的節(jié)點(diǎn)間隔會呈現(xiàn)出不同速遞區(qū)域的差異化現(xiàn)象,我們將在后續(xù)模型測試中展示).圖9展示了規(guī)則網(wǎng)格剖分結(jié)果,其正交方向的網(wǎng)格間隔統(tǒng)一為10 m.圖9中紅線部分為高斯起伏界面位置,由于規(guī)則網(wǎng)格剖分無法精確捕捉真實(shí)的界面位置,因此實(shí)際數(shù)值模擬過程中,將選取縱向最接近真實(shí)界面位置的點(diǎn)作為該模型的自由邊界頂點(diǎn).

兩類高斯模型測試過程中,時間間隔統(tǒng)一為0.5 ms,采樣點(diǎn)數(shù)為2001,總時長為1 s.震源選取主頻為20 Hz的雷克子波,震源位置坐標(biāo)為 (0.4 km, 0 km).與常規(guī)FDM相比(如圖10所示),圖10展示了利用QR-RBF-FDM所得的不同時刻凸型高斯起伏地表模型波場快照,圖8d和圖11中對地震相做了詳細(xì)的分析,并在記錄上對波的類型做了相應(yīng)標(biāo)注.對比可知,常規(guī)FDM和本文方法均能夠較好地處理起伏地表,波場快照中能夠清晰地呈現(xiàn)出直達(dá)P波、反射PP波、反射PS波和面波R.但與常規(guī)FDM不同,QR-RBF-FDM在起伏構(gòu)造中沒有產(chǎn)生臺階散射,模擬結(jié)果更準(zhǔn)確,反射波同相軸更清晰.對比二者數(shù)值結(jié)果可得以下結(jié)論：(1)在0.2 s之前,波在水平地表傳播,QR-RBF-FDM和常規(guī)FDM的模擬結(jié)果基本一致;(2)0.2 s至0.6 s之間,波經(jīng)過了起伏地表,常規(guī)FDM模擬的波場中產(chǎn)生大量“波紋狀”的次級散射波,這些是由規(guī)則網(wǎng)格離散所造成的臺階散射效應(yīng).與之不同,QR-RBF-FDM的地表節(jié)點(diǎn)分布完全貼合了起伏地形,因此QR-RBF-FDM模擬的波場不存在上述次級散射,展現(xiàn)出了該方法對于復(fù)雜區(qū)域模擬的優(yōu)勢;(3)0.6 s之后,波穿過起伏地表并在水平地表繼續(xù)傳播,該階段QR-RBF-FDM和常規(guī)FDM的模擬結(jié)果基本相同,該結(jié)論與0.2 s之前的情況一致.圖12和圖13展示了凹型高斯起伏地表模型波場模擬結(jié)果,對比該結(jié)果可以得出與凸型高斯起伏地表模型類似的結(jié)論.進(jìn)一步地,對比圖11和圖13中波場快照可知,凸起和凹陷模型在傳播過程中均會產(chǎn)生直達(dá)波P、反射縱波PP、反射橫波PS和面波R.與凸起模型相比,0.6時地震波圍繞凹陷模型成回轉(zhuǎn)形式向前擴(kuò)散,無論直達(dá)波或是面波在正方向產(chǎn)生了很強(qiáng)的散射.圖14和圖15分別為兩類模型采用QR-RBF-FDM和FDM模擬得到的地震記錄,從中可以很清晰地看出QR-RBF-FDM相較于FDM在消除階梯散射上的優(yōu)勢,且面波同相軸清晰.綜上可得,常規(guī)FDM和QR-RBF-FDM方法有很大的區(qū)別,QR-RBF-FDM得到的波場在起伏界面處非常穩(wěn)定,并且不會產(chǎn)生大量次級散射,能夠準(zhǔn)確地對自由邊界條件下的起伏地表模型進(jìn)行數(shù)值模擬.

圖7 均勻介質(zhì)起伏地表模型(a) 凸型高斯模型; (b) 凹型高斯模型.

圖8 均勻介質(zhì)起伏地表模型自適應(yīng)起伏地表節(jié)點(diǎn)剖分結(jié)果(a) 凸型高斯模型; (b) 凹型高斯模型.

圖10 凸型高斯起伏地表模型常規(guī)FDM波場快照(a) 0.2 s水平分量; (b) 0.4 s水平分量; (c) 0.6 s水平分量; (d) 0.2 s垂直分量; (e) 0.4 s垂直分量; (f) 0.6 s垂直分量.

圖11 凸型高斯起伏地表模型QR-RBF-FDM波場快照(a) 0.2 s水平分量; (b) 0.4 s水平分量; (c) 0.6 s水平分量; (d) 0.2 s垂直分量; (e) 0.4 s垂直分量; (f) 0.6 s垂直分量.

6.2 起伏凹陷模型

本節(jié)利用起伏凹陷模型驗(yàn)證非均勻模型下自適應(yīng)起伏地表節(jié)點(diǎn)剖分法和QR-RBF-FDM的穩(wěn)定性和適應(yīng)性,此時地表及矩形邊界采用無網(wǎng)格彈性波混合吸收條件.該模型見圖16,模型大小為6 km×4 km,地表起伏函數(shù)為6.1節(jié)兩類高斯函數(shù)的組合.圖17展示了自適應(yīng)起伏地表節(jié)點(diǎn)剖分法離散的結(jié)果.與圖17b中常規(guī)節(jié)點(diǎn)剖分方法采用統(tǒng)一尺度步長相比,圖17a基于自適應(yīng)起伏地表節(jié)點(diǎn)剖分方法中節(jié)點(diǎn)間距隨著速度增加而變大,范圍在10 m至15 m之間.基于該模型進(jìn)行數(shù)值模擬,時間步長為0.5 ms,記錄時間長度為3 s,選取主頻為20 Hz的雷克子波作為震源時間函數(shù),檢波點(diǎn)均勻布設(shè)于地表,接收范圍為4 km,相鄰檢波點(diǎn)間距為10 m.

圖18給出了不同時刻QR-RBF-FDM模擬的波場快照,震源位于地表橫向3 km處.對比不同時刻波場可知,在近地表區(qū)域模擬的波場不存臺階散射,波傳播至深層區(qū)域時(即節(jié)點(diǎn)間隔增大區(qū)域)仍然穩(wěn)定且無頻散現(xiàn)象.由圖19的不同炮點(diǎn)位置(地表x方向坐標(biāo)1 km, 3 km 和5 km處)激發(fā)的地震記錄可以觀察到,不論在地表任何位置激發(fā)的震源,QR-RBF-FDM都能很好地處理起伏地表,采用無網(wǎng)格彈性波混合吸收條件后地震記錄中可以清晰地看出直達(dá)波,反射波,轉(zhuǎn)換波的發(fā)育而無面波與散射波,驗(yàn)證了方法的正確性.綜上可知,本文提出的QR-RBF-FDM不僅通過靈活的自適應(yīng)起伏地表節(jié)點(diǎn)剖分法提供一種穩(wěn)定的離散節(jié)點(diǎn)模型,而且與常規(guī)RBF-FDM相比無需選取形變參數(shù),大幅提升了無網(wǎng)格算法的計算效率和計算精度.與常規(guī)曲坐標(biāo)系算法相比,QR-RBF-FDM差分形式更為簡單,更易于推廣.

圖12 凹型高斯起伏地表模型常規(guī)FDM波場快照(a) 0.2 s水平分量; (b) 0.4 s水平分量; (c) 0.6 s水平分量; (d) 0.2 s垂直分量; (e) 0.4 s垂直分量; (f) 0.6 s垂直分量.

圖13 凹型高斯起伏地表模型QR-RBF-FDM波場快照(a) 0.2 s水平分量; (b) 0.4 s水平分量; (c) 0.6 s水平分量; (d) 0.2 s垂直分量; (e) 0.4 s垂直分量; (f) 0.6 s垂直分量.

圖14 凸型高斯起伏地表模型QR-RBF-FDM和FDM地震記錄(a) QR-RBF-FDM水平分量; (b) QR-RBF-FDM垂直分量; (c) FDM水平分量; (d) FDM垂直分量.

圖16 起伏凹陷模型(a) 縱波速度模型; (b) 橫波速度模型.

圖17 起伏凹陷模型節(jié)點(diǎn)剖分結(jié)果(a) 自適應(yīng)起伏地表節(jié)點(diǎn)剖分結(jié)果; (b) 常規(guī)節(jié)點(diǎn)剖分結(jié)果.

圖18 起伏凹陷模型QR-RBF-FDM波場快照(a) 0.5 s水平分量; (b) 1.0 s水平分量; (c) 1.5 s水平分量; (d) 0.5 s垂直分量; (e) 1.0 s垂直分量; (f) 1.5 s垂直分量.

圖19 起伏凹陷模型QR-RBF-FDM地震記錄炮點(diǎn)位于地表1 km：(a)水平分量;(b)垂直分量; 炮點(diǎn)位于地表3 km：(c)水平分量;(d)垂直分量; 炮點(diǎn)位于地表5 km：(e)水平分量;(f)垂直分量.

圖20 起伏凹陷模型多炮成像結(jié)果(a) PP成像; (b) PS成像.

圖21 起伏地表Marmousi-2模型(a) 縱波速度模型; (b) 橫波速度模型.

此處,我們進(jìn)一步地利用上述方法實(shí)現(xiàn)了ERTM,采用60炮主頻20 Hz的雷克子波作為震源時間函數(shù)在地表附近放炮,中間激發(fā)兩端接收,炮間距100 m.最大記錄時長3 s,每炮400道接收,道間距10 m.圖20給出了多炮疊加后的PP和PS剖面,從中可以看出,QR-RBF-RTM可靈活地處理起伏地表,且地下界面位置準(zhǔn)確,同相軸清晰連續(xù),該數(shù)值結(jié)果驗(yàn)證了本文方法的正確性.

6.3 起伏地表的Marmousi模型

本節(jié)利用經(jīng)典Marmousi-2模型驗(yàn)證復(fù)雜模型下自適應(yīng)起伏地表節(jié)點(diǎn)剖分法和QR-RBF-FDM的穩(wěn)定性和適應(yīng)性,此時地表采用無網(wǎng)格彈性波混合吸收條件.模型見圖21,模型大小為13.6 km×2.8 km,采用Foothill模型的起伏地形作為Marmousi的地形起伏.圖22展示了自適應(yīng)起伏地表節(jié)點(diǎn)剖分法離散的結(jié)果.圖22基于自適應(yīng)起伏地表節(jié)點(diǎn)剖分方法中節(jié)點(diǎn)間距隨著速度增加而變大,范圍在5～18 m之間.基于該模型進(jìn)行數(shù)值模擬,時間步長為0.5 ms,記錄時間長度為4.5 s,選取主頻為20 Hz的雷克子波作為震源時間函數(shù),檢波點(diǎn)均勻布設(shè)于地表,接收范圍為10 km,相鄰檢波點(diǎn)間距為10 m.

圖23給出了不同時刻水平和垂直分量的波場快照,震源位于地表橫向3 km處.圖24為縱橫波分離的波場快照,其中左列為縱波波場分量,右列為橫波波場分量.圖23表明本文方法對于地表高程變化劇烈的復(fù)雜地表模型能夠?qū)崿F(xiàn)高精度數(shù)值模擬,圖24表明QR-RBF-FD應(yīng)用于分離方程能夠有效分離縱橫波.由圖25地震記錄可見,在復(fù)雜模型下彈性波場仍然穩(wěn)定,且無矩形網(wǎng)格剖分引起的虛假反射現(xiàn)象.由圖25的不同炮點(diǎn)位置(地表x方向坐標(biāo)3.2 km, 6.5 km和9.8 km處)激發(fā)的地震記錄可以觀察到,不論在地表任何位置激發(fā)的震源,在記錄中均可以清晰地看出,直達(dá)波、反射波和轉(zhuǎn)換波的發(fā)育.由于起伏地表的存在,反射波和轉(zhuǎn)換波的雙曲線特性發(fā)生規(guī)則扭曲,驗(yàn)證了在復(fù)雜構(gòu)造和起伏地表情況下QR-RBF-FDM的適應(yīng)性.

此處,我們進(jìn)一步地利用上述方法實(shí)現(xiàn)了ERTM,采用135炮主頻20 Hz的雷克子波作為震源時間函數(shù)在地表附近放炮,中間激發(fā)兩端接收,炮間距100 m,最大記錄時長4.5 s,每炮1000道接收,道間距10 m.圖26給出了多炮疊加后的PP和PS成像結(jié)果剖面,從中可以看出,成像結(jié)果與偏移速度模型的構(gòu)造界面信息吻合,且地下界面位置準(zhǔn)確,無論是PP成像還是PS額成像,同相軸均清晰連續(xù),具有較高的信噪比,表明基于QR-RBF- ERTM方法對起伏地表下得到復(fù)雜模型具有良好的適應(yīng)性.

7 結(jié)論

本文提出一種自適應(yīng)節(jié)點(diǎn)剖分法實(shí)現(xiàn)起伏地表無網(wǎng)格剖分,在此基礎(chǔ)上提出QRRBF-FDM,推導(dǎo)出彈性波模擬中RBF-FDM的頻散關(guān)系和穩(wěn)定性條件.同時,提出了一種適用于起伏地表的無網(wǎng)格起伏自由地表條件和彈性波混合邊界條件,實(shí)現(xiàn)了基于無網(wǎng)格的起伏地表彈性波數(shù)值模擬.此外,推導(dǎo)出一種適用于無網(wǎng)格的精確矢量波分離方法,實(shí)現(xiàn)了縱橫波分離的無網(wǎng)格彈性波逆時偏移方法.通過對高斯起伏模型,起伏凹陷模型和起伏地表的Marmousi-2模型的數(shù)值模擬和偏移試算,得出到以下幾點(diǎn)認(rèn)識：

(1)自適應(yīng)節(jié)點(diǎn)剖分法不僅能與地表完全匹配,而且節(jié)點(diǎn)分布密度隨地下介質(zhì)的速度變化而變化,在速度大的區(qū)域節(jié)點(diǎn)分布稀疏,速度小的區(qū)域節(jié)點(diǎn)分布密集.同時采用自適應(yīng)節(jié)點(diǎn)控制方程對節(jié)點(diǎn)位置進(jìn)行調(diào)整,保證節(jié)點(diǎn)分布的均勻性.基于該方法,利用無形變參數(shù)擾動QR-RBF-FDM進(jìn)行數(shù)值模擬,既保持了傳統(tǒng)FDM的高精度,又彌補(bǔ)了傳統(tǒng)FDM在自由起伏地表處理中缺乏的靈活性的缺點(diǎn).

(2)在起伏地表條件下,貼體網(wǎng)格與曲網(wǎng)格是模擬地震波的重要方法,其網(wǎng)格易于生成,而且可利用水平地表假設(shè)的自由邊界條件,在曲坐標(biāo)系下直接使用.無網(wǎng)格具有利用彈性波、聲波等方程進(jìn)行直接差分的優(yōu)點(diǎn),無需坐標(biāo)變換,后期研究中更適用于復(fù)雜的黏性介質(zhì)或各向異性方程.需要指出的是,無網(wǎng)格數(shù)值模擬的差分模板仍基于FDM的數(shù)值差分思想,其傳播的穩(wěn)定性依賴于差分模板的均勻性和連續(xù)性,因此在現(xiàn)階段仍然無法突破梯度非連續(xù)模型,這將是后續(xù)研究工作的重點(diǎn)與難點(diǎn).

(3)基于Gu等(2019)的基礎(chǔ)上,本文推導(dǎo)并實(shí)現(xiàn)了無網(wǎng)格模型參數(shù)化條件下的起伏地表彈性波逆時偏移.模型試算結(jié)果表明,該方法對地表附近和復(fù)雜深層成像效果較好,而且震源的干擾效應(yīng)能準(zhǔn)確反演起伏地表形態(tài).

(4)應(yīng)用QR-RBF-FDM,以無網(wǎng)格模型參數(shù)化為基礎(chǔ),實(shí)現(xiàn)逆時偏移,對未來起伏地表復(fù)雜構(gòu)造、復(fù)雜介質(zhì)處理具有廣闊的應(yīng)用前景;同時,QR-RBF-FDM也可直接應(yīng)用于最小二乘偏移、全波形反演.

圖22 起伏地表Marmousi-2模型自適應(yīng)起伏地表節(jié)點(diǎn)剖分結(jié)果

圖23 起伏地表Marmousi-2模型RBF-FDM波場快照(a) 1.0 s水平分量; (b) 1.0 s垂直分量; (c) 2.0 s水平分量; (d) 2.0 s垂直分量; (e) 3.0 s水平分量; (f) 3.0 s垂直分量.

圖24 起伏地表Marmousi-2模型RBF-FDM波場快照(a) 1.0 s縱波分量; (b) 1.0 s橫波分量; (c) 2.0 s縱波分量; (d) 2.0 s橫波分量; (e) 3.0 s縱波分量; (f) 3.0 s橫波分量.

圖25 起伏地表Marmousi-2模型RFB-FDM地震記錄炮點(diǎn)位于地表3.2 km： (a) 水平分量; (b) 垂直分量; 炮點(diǎn)位于地表6.5 km： (c) 水平分量; (d) 垂直分量; 炮點(diǎn)位于地表9.8 km： (e) 水平分量; (f) 垂直分量.

圖26 起伏地表Marmousi-2模型彈性波逆時偏移結(jié)果 (a) PP成像; (b) PS成像.