徐劍俠， 張偉， 陳曉非*

1 中國(guó)科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院， 合肥 230026 2 南方科技大學(xué)地球與空間科學(xué)系， 廣東深圳 518055

0 引言

有限差分算法是地震學(xué)中廣泛應(yīng)用的重要算法(Alterman and Karal, 1968; Bayliss et al., 1986; Saenger and Bohlen, 2004)，可以簡(jiǎn)潔高效地模擬二維/三維非均勻模型中地震波波場(chǎng)，尤其是在Virieux(1986)提出了基于應(yīng)力-速度的一階差分格式后得到了廣泛的應(yīng)用.

地形起伏會(huì)明顯影響地震波場(chǎng)的傳播.為了在有限差分中精確處理起伏地表，需要解決網(wǎng)格劃分和自由地表邊界條件引入兩個(gè)問(wèn)題.對(duì)于網(wǎng)格劃分，現(xiàn)有工作已經(jīng)研究比較完善，主要有規(guī)則網(wǎng)格，非結(jié)構(gòu)網(wǎng)格和垂向變換/貼體網(wǎng)格等方法.規(guī)則網(wǎng)格容易生成和使用，最簡(jiǎn)單的方法是在加密網(wǎng)格后采用臺(tái)階狀網(wǎng)格近似描繪起伏地表(Robertsson, 1996; Wang and Zhang, 2004)；非結(jié)構(gòu)網(wǎng)格適用性好，較易生成網(wǎng)格，但是對(duì)精度有一定的影響(K?ser et al., 2001)；垂向變換/貼體網(wǎng)格思路相近，依照實(shí)際地表起伏劃分曲線(xiàn)網(wǎng)格，使用比較方便，適用于同位網(wǎng)格和交錯(cuò)網(wǎng)格(Hestholm and Ruud, 1994, 1998; Zhang and Chen, 2006)，但是對(duì)網(wǎng)格變形范圍具有一定的限制(Tessmer et al., 1992; Hestholm and Ruud, 1994).自由地表邊界條件近年來(lái)具有較大發(fā)展.對(duì)于速度分量法向?qū)?shù)的需求，多通過(guò)牽引力為零邊界條件，轉(zhuǎn)化為速度分量水平導(dǎo)數(shù)計(jì)算.對(duì)于應(yīng)力分量法向?qū)?shù)計(jì)算，主要實(shí)現(xiàn)方法大體有：?jiǎn)芜叢罘址椒▽?shí)施方便，通過(guò)在地表格點(diǎn)上采用單邊差分計(jì)算，可以滿(mǎn)足全局二階空間差分精度，并被成功拓展到復(fù)雜介質(zhì)(Nilsson et al., 2007; Appel? and Petersson, 2009; Lan and Zhang, 2011,2012; 蘭海強(qiáng)等, 2011; Liu et al., 2018)；應(yīng)力鏡像法是地表邊界條件在水平規(guī)則地表下的簡(jiǎn)化表述，因此直接用于起伏地表存在一定近似(例如Robertsson, 1996)，部分工作也考慮了起伏地表法向?qū)恳Φ挠绊?，和牽引力鏡像法物理思路一致(Solano et al., 2016)；牽引力鏡像法考慮了起伏地表法向量變化，通過(guò)牽引力關(guān)于地表反對(duì)稱(chēng)表述自由地表邊界條件，對(duì)于自由地表邊界條件表述更加準(zhǔn)確，在起伏地表邊界滿(mǎn)足四階空間差分精度(Zhang and Chen, 2006; Zhang et al., 2012b)；Taylor多項(xiàng)式方法適用于規(guī)則網(wǎng)格，通過(guò)引入Taylor多項(xiàng)式擬合外推虛擬點(diǎn)內(nèi)波場(chǎng)信息，同時(shí)考慮了自由地表法向影響(Lombard et al., 2008)；在地球物理勘探領(lǐng)域不關(guān)注地表反射，也有工作將吸收層作用于起伏地表(Rao and Wang, 2013, 2018)；還有工作通過(guò)樣條或者插值多項(xiàng)式外推物理量，進(jìn)而處理自由地表邊界條件(Mulder, 2017; Mulder and Huiskes, 2017)等.上述方法已經(jīng)可以在笛卡爾坐標(biāo)系框架下較好地處理起伏地表邊界條件，其中基于貼體網(wǎng)格牽引力鏡像法(Zhang and Chen, 2006)在貼體網(wǎng)格中將牽引力關(guān)于地表反對(duì)稱(chēng)，可以簡(jiǎn)潔清晰地表征自由地表物理邊界條件，同時(shí)具有較高空間差分精度.

當(dāng)研究區(qū)域-全球尺度的地震波傳播問(wèn)題時(shí)，必須要考慮地球曲率的影響.雖然可以通過(guò)網(wǎng)格變形(Appel? and Petersson, 2009)或者展平變換(Li et al., 2014)等方法處理，但是選用極坐標(biāo)系(二維)或球坐標(biāo)系(三維)處理起來(lái)更為簡(jiǎn)潔方便(Igel and Weber, 1995; Toyokuni and Takenaka, 2006; Zhang et al., 2012a; Wang et al., 2013).

目前在極/球坐標(biāo)系中，有限差分法不能較好處理起伏地表邊界條件的，所以我們希望可以將牽引力鏡像法(Zhang and Chen, 2006)推廣到極/球坐標(biāo)系有限差分中.作為第一步工作，本文集中于將牽引力鏡像法推廣到二維極坐標(biāo)系有限差分中處理起伏地表邊界條件.

1 數(shù)學(xué)原理

1.1 規(guī)則地表情況下極坐標(biāo)系方程

(1)

運(yùn)動(dòng)學(xué)方程和廣義胡克定律可以寫(xiě)為

σ=C:ε,

(2)

(3)

1.2 起伏地表情況下極坐標(biāo)系方程

當(dāng)存在起伏地形情況下，離散網(wǎng)格和實(shí)際地形不一致將引起的虛假波形，因此采用貼體網(wǎng)格貼合實(shí)際地形的有限差分(Fornberg, 1988; Zhang and Chen, 2006)具有更高精度.參照前人工作(Zhang et al., 2012a)，以ζ為橫向坐標(biāo)η為徑向坐標(biāo)，引入貼合地形的貼體網(wǎng)格ζ=ζ(r,θ),η=η(r,θ)，將物理空間中不規(guī)則的貼體網(wǎng)格映射到計(jì)算空間中的規(guī)則網(wǎng)格，參見(jiàn)圖 1.應(yīng)用鏈?zhǔn)角髮?dǎo)法則，(1)式中的哈密頓算子和極坐標(biāo)基矢量關(guān)系變?yōu)?4)式形式：

(4)

通過(guò)(4)式中正交基矢量的微分關(guān)系，可以得到坐標(biāo)基矢量在貼體網(wǎng)格坐標(biāo)系中的空間導(dǎo)數(shù)，參見(jiàn) (5)式：

圖1 物理網(wǎng)格-計(jì)算網(wǎng)格示意圖通過(guò)ζ=ζ(r,θ),η=η(r,θ)將不規(guī)則的物理網(wǎng)格(r,θ)映射為規(guī)則的計(jì)算網(wǎng)格(η,ζ).Fig.1 Physical-computational spaceThe operators ζ=ζ(r,θ),η=η(r,θ) mapping from irregular physical grid (r,θ) to regular computational grid (η,ζ).

(5)

通過(guò)r和θ的正交關(guān)系，可以得到貼體網(wǎng)格的空間空間變形系數(shù)矩陣：

(6)

將(4)、(5)、(6)式代入(2)式，經(jīng)過(guò)化簡(jiǎn)可以得到貼體網(wǎng)格下彈性動(dòng)力學(xué)方程組：

(7)

1.3 極坐標(biāo)系下貼體網(wǎng)格中牽引力鏡像法

(8)

將(8)式代入牽引力公式后，自由地表處牽引力為0的條件可以寫(xiě)為

Tη=nη·σ=0,

(9)

在對(duì)(9)式進(jìn)行時(shí)間域求導(dǎo)后，將(7)式中應(yīng)力變量時(shí)間導(dǎo)數(shù)表達(dá)式代入，整理后可以得到速度分量在η方向?qū)?shù)使用速度分量、速度分量在ξ方向?qū)?shù)表達(dá)的具體形式，參見(jiàn)(10)式.由于體貼網(wǎng)格中所有變量在ξ方向的導(dǎo)數(shù)全部可解，所以自由地表處速度分量在η方向?qū)?shù)可以由(10)式計(jì)算.

(10)

由于自由地表邊界條件此時(shí)僅有牽引力關(guān)于地表反對(duì)稱(chēng)的假設(shè)，所以需要把(7)式中所有應(yīng)力分量η方向的導(dǎo)數(shù)通過(guò)牽引力及其鏡像表達(dá).將(9)式中牽引力表達(dá)式對(duì)η求導(dǎo)，通過(guò)使用(5)式和(6)式化簡(jiǎn)整理可得(11)式.至此，(7)式中自由地表處所有應(yīng)力分量η方向?qū)?shù)，可以用組合的形式，由應(yīng)力分量、牽引力分量及牽引力分量η方向?qū)?shù)表示.

(11)

至此，在自由地表處(7)式中所涉及的自由地表邊界條件施加完成.

1.4 其他數(shù)值計(jì)算處理

仿照Z(yǔ)hang和Chen(2006)，本研究所用空間差分模板采用DRP/opt MacCormack格式(Hixon, 1998)，時(shí)間步積分采用4-6階Runge-Kutta積分(Hu et al.,1996).

對(duì)于自由地表下臨近格點(diǎn)，如果差分模板跨過(guò)自由地表，此時(shí)速度/應(yīng)力分量η方向?qū)?shù)需要額外處理.其中應(yīng)力分量η方向?qū)?shù)可以繼續(xù)沿用(11)式，但是速度分量η方向?qū)?shù)無(wú)法采用(10)式求解，此處沿用Zhang和Chen(2006)處理方法，采用緊致差分格式(Hixon and Turkel, 2000)，在保證4階精度情況下計(jì)算相關(guān)格點(diǎn)速度分量η方向?qū)?shù).

本算法中震源引入方法與Zhang和Chen(2006)一致，空間上的分布為3個(gè)格點(diǎn)半寬度的光滑高斯形，二維分布采用兩個(gè)一維函數(shù)乘積構(gòu)成以提高精度(Tornberg and Engquist, 2004).

本算法在非自由地表的吸收邊界處，采用e指數(shù)衰減層處理(Cerjan et al., 1985).

2 算例對(duì)比

我們需要計(jì)算一系列算例和參考結(jié)果對(duì)比，驗(yàn)證本算法正確性.考慮到Zhang和Chen(2006)可以在直角坐標(biāo)系下很好地處理地表起伏，所以我們將其方法結(jié)果作為參照解.通過(guò)選用相同介質(zhì)模型和貼體網(wǎng)格，分別使用直角坐標(biāo)系和極坐標(biāo)系進(jìn)行計(jì)算，在歸一化后對(duì)二者進(jìn)行對(duì)比檢驗(yàn).

對(duì)于DRP/opt MacCormack格式，傳播20到30個(gè)波長(zhǎng)的距離上穩(wěn)定性條件為每波長(zhǎng)內(nèi)網(wǎng)格點(diǎn)6～8個(gè)(Hixon, 1998).本文選取橫波內(nèi)每波長(zhǎng)網(wǎng)格點(diǎn)25，CFL約0.3，以提高計(jì)算精度.

我們選擇模型為均勻介質(zhì)，縱波波速為3.0 km·s-1，橫波波速為1.2 km·s-1，密度為1800 kg·m-3；震源為爆炸源，震源深度為0.5 km，震源時(shí)間函數(shù)為Ricker子波，中心頻率為0.96 Hz，時(shí)間延遲為1.625 s；計(jì)算時(shí)間步長(zhǎng)為0.005 s；空間網(wǎng)格采用貼體網(wǎng)格劃分，在橫向/θ方向網(wǎng)格間距不大于0.05 km.對(duì)地表起伏光滑變化的情況，垂向/R方向采用垂向變換網(wǎng)格(Ruud and Hestholm, 2001)，網(wǎng)格初始間距為0.05 km，根據(jù)地表起伏進(jìn)行一定的拉伸或壓縮.

在圖2—5中，每幅圖的最上方為地表附近網(wǎng)格劃分示意圖，*代表震源，并用S標(biāo)識(shí)震中位置，A、B則為兩臺(tái)站，方便在震相圖中對(duì)應(yīng)識(shí)別，空心正三角形為密集布置臺(tái)站位置，實(shí)心倒三角為稀疏布置臺(tái)站，圖中網(wǎng)格示意圖為計(jì)算網(wǎng)格抽稀后結(jié)果，以清晰顯示.為方便直觀對(duì)比波形，我們將極坐標(biāo)系下波形結(jié)果，全部轉(zhuǎn)換到相應(yīng)直角坐標(biāo)系下顯示.每幅圖的中間為密集臺(tái)站模擬結(jié)果對(duì)比，以StaVx和StaVz分別標(biāo)識(shí)水平方向和垂直方向震動(dòng)，可以用于檢驗(yàn)本方法波形同相軸的一致性.每幅圖的最下方為波形顯示比例放大3倍后的稀疏臺(tái)站的波形對(duì)比，以ZoomVx和ZoomVz分別標(biāo)識(shí)水平方向和垂直方向震動(dòng)，可以精細(xì)檢驗(yàn)本文方法的振幅和相位.在波形對(duì)比圖中，紅色波形為采用Zhang和Chen(2006)方法在直角坐標(biāo)系下計(jì)算的結(jié)果，黑色波形為采用本文結(jié)果在極坐標(biāo)系下的結(jié)果.在波形圖中，不同臺(tái)站波形，按照臺(tái)站沿自由表面到震中距離排列；所有極坐標(biāo)系下計(jì)算，均設(shè)定震源正上方處地表半徑為100 km，計(jì)算網(wǎng)格為扇形，不包括極坐標(biāo)系原點(diǎn).

圖2 極坐標(biāo)下規(guī)則地表波形對(duì)比Fig.2 Waveform comparison in polar coordinate system with regular topography

圖3 極坐標(biāo)系下起伏地表情況波形對(duì)比Fig.3 Waveform comparison in polar coordinate system with irregular topography

圖4 水平地表情況波形對(duì)比Fig.4 Waveform comparison under horizontal surface

圖5 水平地表疊加光滑起伏地表情況波形對(duì)比Fig.5 Waveform comparison under horizontal surface with irregular topography

首先我們對(duì)比規(guī)則地表和僅存在局部起伏地表兩種情況.圖2為本文算法規(guī)則地表下退化結(jié)果的對(duì)比；圖3為在規(guī)則地表上疊加光滑地形起伏時(shí)結(jié)果對(duì)比，所加山谷深度和山峰高度均為2 km.可以看到在這兩種情況下本文算法和參照結(jié)果吻合，因此本方法適用于局部光滑地形起伏情況下自由地表邊界條件處理.

考慮到本文方法在極坐標(biāo)系下使用貼體網(wǎng)格，同樣可以處理直角坐標(biāo)系下水平地表情況，此時(shí)的水平地表在極坐標(biāo)系中可以認(rèn)為是全局光滑變形的起伏地表，對(duì)比結(jié)果參見(jiàn)圖 4.圖 5對(duì)比了在極坐標(biāo)系中全局光滑變形基礎(chǔ)上，疊加小尺度變形的波形對(duì)比結(jié)果.本文方法在這些對(duì)比中和參照結(jié)果吻合，因此本方法適用于全局光滑地形起伏情況下自由地表邊界條件處理，同時(shí)也適用于全局光滑地形起伏和局部地形起伏耦合情況.

上述各算例中，本文方法計(jì)算波形和參考算法結(jié)果高度一致，說(shuō)明本文方法準(zhǔn)確、有效，可以精確處理極坐標(biāo)系有限差分方法中起伏地表處自由邊界條件.

3 總結(jié)與展望

本文基于極坐標(biāo)系貼體網(wǎng)格有限差分，引入貼體網(wǎng)格和牽引力鏡像方法處理起伏地表邊界條件，并和直角坐標(biāo)系下已有方法對(duì)比，算法得到了較全面檢驗(yàn)，證明本文方法正確、有效.

本文針對(duì)二維極坐標(biāo)系，利用鏈?zhǔn)角髮?dǎo)法則展開(kāi)哈密頓算子到貼體網(wǎng)格坐標(biāo)系中，并利用矢量正交性質(zhì)化簡(jiǎn)方程.這些方法全部可以方便地拓展到三維球坐標(biāo)系，可以預(yù)期拓展后的方法對(duì)更加精確模擬實(shí)際地震波場(chǎng)在地球中傳播具有重要意義.

本文研究關(guān)注點(diǎn)為起伏地表邊界條件，模擬區(qū)域均為扇形區(qū)域，由于極坐標(biāo)系固有的周期性，在施加循環(huán)邊界條件后，本工作模擬區(qū)域可以方便地從扇形區(qū)域擴(kuò)展到完整圓環(huán)區(qū)域；現(xiàn)階段本文方法未包含極坐標(biāo)系原點(diǎn)，因此無(wú)法模擬穿過(guò)原點(diǎn)及附近區(qū)域的相關(guān)震相，例如天然地震中的PKIKP震相，也需要進(jìn)一步完善，Wang等(2001)采用的球心處理方法，可以為本文方法完善提供重要參考.

致謝感謝南方科技大學(xué)張振國(guó)老師對(duì)本文完成提供的參考意見(jiàn)和幫助.