劉彬， 榮棉水， 孫偉家， 黃建平， 巴振寧

1 中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院， 青島 266580 2 北京工業(yè)大學(xué)城市建設(shè)學(xué)部， 北京 100124 3 中國科學(xué)院地質(zhì)與地球物理研究所, 地球與行星物理重點(diǎn)實(shí)驗(yàn)室, 北京 100029 4 天津大學(xué)土木工程系， 天津 300372

0 引言

地球深部圈層及沉積盆地是一種多層分區(qū)非均勻介質(zhì)系統(tǒng)，其中大尺度不規(guī)則邊界的反射和透射決定了地震波的傳播特征，而小尺度層內(nèi)隨機(jī)非均勻介質(zhì)引起波的散射和衰減.近一個多世紀(jì)以來，針對這種復(fù)雜介質(zhì)系統(tǒng)地震波傳播的理論和數(shù)值研究得到了極大的發(fā)展，有限差分法、有限元法和譜元法等普適的全數(shù)值模擬技術(shù)得到了廣泛應(yīng)用.然而，這些數(shù)值模擬方法均需對研究區(qū)域進(jìn)行空間離散，有限尺寸的空間離散對不規(guī)則邊界和體非均勻介質(zhì)的近似程度有限，影響了地震模擬的精度，而半解析的數(shù)值模擬方法處理這些問題具有巨大的優(yōu)勢，能根據(jù)介質(zhì)非均質(zhì)程度和計算波長來調(diào)整模擬精度.本文擬綜合幾種半解析數(shù)值模擬方法，利用其優(yōu)勢，形成多層分區(qū)非均勻介質(zhì)地震波傳播模擬的最佳解決方案.

大多數(shù)半解析的地震數(shù)值模擬方法的建立源于邊界積分方程，其優(yōu)勢在于精確刻畫不規(guī)則邊界和解析實(shí)施地層界面邊界條件，實(shí)現(xiàn)對不規(guī)則界面反射、透射和散射的精確模擬.這些半解析數(shù)值模擬方法主要分為三類：早期的Aki-Larner方法(Aki and Larner, 1970; Bard, 1982)、Bouchon-Campillo 離散波數(shù)方法(Bouchon, 1982, 1985; Campillo and Bouchon, 1985; Kawase, 1988; Bouchon et al., 1989)和直接(Dravinski, 1983; 符力耘和牟永光, 1994; Fu et al., 2002)與間接(Sánchez-Sesma and Esquivel, 1979; Sánchez-Sesma and Campillo, 1991; Ba and Liang, 2010; Ba and Xiao, 2016)的邊界元法.基于含有邊界積分與體積分的廣義Lippmann-Schwinger積分方程，這些半解析數(shù)值方法進(jìn)一步拓展至模擬含有不規(guī)則邊界和體非均勻性的層狀介質(zhì)地震波傳播(Fu, 2002; Fu and Bouchon, 2004; 管西竹等，2011).然而，基于積分方程的全波形數(shù)值解數(shù)據(jù)存儲和計算耗時巨大，不適用于超大模型模擬和高頻記錄合成.為此，發(fā)展了基于波場連接技術(shù)的積分方程數(shù)值模擬方法(Fu and Wu, 2001; Ge et al., 2005)和各種積分方程逼近求解地震模擬方法(Hu et al., 2009; Yu et al., 2010; Yu and Fu, 2012, 2013, 2014). 這些逼近方法以損失多次散射的模擬精度為代價來有效提高計算效率，適用于介質(zhì)速度橫向變化不大的地震模擬.通過積分方程的頻率-波數(shù)域表征，上述半解析積分方程數(shù)值模擬方法可轉(zhuǎn)換為高效率的R/T遞推傳播矩陣法.該發(fā)展方向的主要貢獻(xiàn)源于Kennett (1974)極大發(fā)展了早期的Thomson-Haskell矩陣方法(Gilbert and MacDonald, 1960)，形成水平層狀介質(zhì)地震模擬的經(jīng)典方法(Kennett, 1983; Luco and Apsel, 1983; Maupin and Kennett, 1987)，在地震學(xué)領(lǐng)域得到廣泛應(yīng)用.該方法進(jìn)一步拓展至模擬不規(guī)則邊界地層，形成不規(guī)則邊界均勻地層廣義R/T遞推傳播矩陣法(Chen, 1990, 1995, 1996; Cao et al., 2004; Ge and Chen, 2007, 2008)，其優(yōu)勢在于精確刻畫不規(guī)則地層界面，精確模擬界面的反射、透射和散射效應(yīng).最近，Liu等(2020)在該方法基礎(chǔ)上，利用廣義Lippmann-Schwinger積分方程的頻率-波數(shù)域表征(Fu, 2006)，結(jié)合R/T遞推傳播矩陣算法，采用薄板離散技術(shù)，建立了非均質(zhì)地層薄板化全局廣義R/T遞推傳播矩陣法.其優(yōu)勢在于：可根據(jù)地震模擬波長來調(diào)整薄板厚度或根據(jù)速度橫向變化程度采用高階薄板近似，準(zhǔn)確模擬地層內(nèi)部隨機(jī)非均質(zhì)體散射引起的地震波衰減效應(yīng).但由于采用薄板離散，該方法弱化了對不規(guī)則地層界面的精細(xì)刻畫，對界面反射和透射的模擬精度有待研究.

本文擬吸收上述兩種廣義R/T遞推傳播矩陣方法各自的優(yōu)勢，形成一種非均勻介質(zhì)、不規(guī)則邊界的全局廣義R/T遞推傳播矩陣混合方法，采用該混合方法求解二維情況下邊界不規(guī)則、層內(nèi)非均質(zhì)的復(fù)雜模型中SH地震波場，特別是用于起伏近地表復(fù)雜介質(zhì)地震波場的模擬.重點(diǎn)研究與波長相關(guān)的縱/橫向離散精度以及與非均質(zhì)強(qiáng)度相關(guān)的模擬精度，并開展了與邊界元法的比較研究.

1 兩種全局廣義R/T遞推傳播矩陣地震模擬方法

均勻介質(zhì)的不規(guī)則邊界全局廣義R/T遞推傳播矩陣(GGRTM)法(Chen,1990,1995,1996)將積分方程的位移項(xiàng)和應(yīng)力項(xiàng)在邊界上展開，使對位移項(xiàng)的求解轉(zhuǎn)換為對展開系數(shù)的求解，再通過引入全局廣義反射/透射系數(shù)矩陣并建立該矩陣與界面相關(guān)參數(shù)矩陣的聯(lián)系來實(shí)現(xiàn)波場中任意位置位移的求解.非均勻介質(zhì)全局廣義R/T遞推傳播矩陣法是不規(guī)則邊界均勻介質(zhì)廣義R/T遞推傳播矩陣法在非均勻介質(zhì)情況下的推廣，其基本原理為：首先將整個計算模型進(jìn)行薄層離散，之后將對含不規(guī)則邊界的薄層邊界積分與對非均勻介質(zhì)薄層的體積分相結(jié)合，利用積分方程的離散矩陣同位移和應(yīng)力及上下行波的關(guān)系獲得修正反射/透射矩陣，最后將修正反射/透射矩陣轉(zhuǎn)換為全局廣義反射/透射矩陣并結(jié)合震源層的上下行波實(shí)現(xiàn)波場求解.本文方法的推導(dǎo)過程詳見附錄.

非均勻介質(zhì)全局廣義R/T遞推傳播矩陣方法基于廣義Lipmann-Schwinger積分方程.第一步，先將模型劃分成數(shù)個厚度很小的水平薄層.對于積分方程中的邊界積分項(xiàng)直接使用虛擬邊界上的物理參數(shù)進(jìn)行計算；而對于體積分部分，若模型為一般的沉積地層，在地震勘探頻帶范圍內(nèi)采用一階Born近似，利用每個薄層上下界面的邊界積分來近似替代體積分，實(shí)現(xiàn)體積分項(xiàng)的降次，同時也可以滿足模擬精度要求.而對于橫向速度變化大的薄板(含鹽丘)可采用二階或高階Born近似.第二步，將邊界積分項(xiàng)和體積分項(xiàng)進(jìn)行離散并組裝成矩陣方程組.第三步，將通過矩陣方程組和相鄰薄層間的位移和應(yīng)力關(guān)系，將矩陣方程組轉(zhuǎn)換成全局修正反射/透射矩陣，并進(jìn)一步轉(zhuǎn)換至全局廣義反射/透射矩陣.第四步，利用震源層的總上下行波振幅和全局廣義反射/透射矩陣，在層間使用迭代關(guān)系進(jìn)行任意位置的波場求解.

以上兩種全局廣義R/T遞推傳播矩陣地震模擬方法均屬于半解析的地震數(shù)值模擬方法，其原理和計算流程相似.不規(guī)則界面全局廣義R/T遞推傳播矩陣方法的優(yōu)勢在于能準(zhǔn)確描述模型中的不規(guī)則界面，實(shí)現(xiàn)了對界面的反射、透射及散射效應(yīng)的精確模擬，但其并未考慮內(nèi)部介質(zhì)的非均勻性.而非均勻介質(zhì)全局廣義R/T遞推傳播矩陣法的優(yōu)勢恰恰在于能更好地描述地層內(nèi)部隨機(jī)非均質(zhì)體散射對地震波衰減的影響，但該方法在實(shí)施過程中采用薄層近似來處理含不規(guī)則界面的地層，弱化了對不規(guī)則地層界面的精細(xì)刻畫，因此根據(jù)計算模型的特點(diǎn)采用以上兩種方法混合求解有利于充分發(fā)揮兩種方法的優(yōu)勢.在實(shí)際情況中，地層中不規(guī)則界面與非均勻介質(zhì)往往是同時存在的，當(dāng)需研究的模型主要由不規(guī)則界面地層構(gòu)成且含有隨機(jī)非均勻介質(zhì)層(如巖性強(qiáng)非均質(zhì)的沖積山谷)時，對該非均質(zhì)層可采用基于薄層近似的非均勻介質(zhì)全局廣義R/T遞推傳播矩陣法，其他地層則采用傳統(tǒng)的不規(guī)則界面全局廣義R/T遞推傳播矩陣方法求解，從而通過混合實(shí)施兩種方法能實(shí)現(xiàn)對同一模型中不規(guī)則界面與介質(zhì)非均勻性的準(zhǔn)確模擬.

2 模擬精度分析

在本節(jié)中，我們對多種工況進(jìn)行了正演計算，比較了不同剖分精度以及內(nèi)部介質(zhì)不同速度擾動量情況下的模擬結(jié)果.圖1給出了一個已知各地層背景速度的二維層狀模型，該模型第一、二層分界面為傾斜界面，模型從上至下各層速度分別為2700 m·s-1、3000 m·s-1、3600 m·s-1和4000 m·s-1，地震子波采用主頻為15 Hz的雷克子波，震源位于模型橫向距離1500 m，深度800 m處，101個檢波器以30 m的間隔均勻分布在地表.為了討論本文方法在不同剖分精度下的計算情況，我們首先對圖1所示的模型根據(jù)單波長內(nèi)不同單元數(shù)量來進(jìn)行離散.

圖1 一個含傾斜界面的二維層狀模型Fig.1 A 2-D layered model with an oblique interface

圖2 均勻情況下不同橫向剖分精度時第51道振幅曲線比較(自上至下分別為10、5、3、2單元每波長)Fig.2 51st trace amplitude curve comparison between different horizontal discrete precision in homogeneous situation (from top to bottom: 10,5,3,2 elements per wavelength)

在本節(jié)的研究中，我們分別考慮了橫向和縱向的離散精度.首先考慮均勻介質(zhì)情形下的橫向離散精度.對于圖1所示模型，各地層中最低波速為2700 m·s-1，當(dāng)震源主頻為15 Hz時相應(yīng)的最短波長為180 m.一般情況下，單位波長內(nèi)單元剖分?jǐn)?shù)量為3個即可保證計算精度(Campillo,1987; Fu, 2002)，即剖分精度為60 m時可滿足條件.為揭示不同離散精度下計算結(jié)果的差異，我們將單個波長內(nèi)單元剖分的數(shù)量分別取為10、5、3、2，并分別計算這四種情況下的地震波場.為方便對比，這里抽取了位于模型正中的第51道的振幅曲線進(jìn)行比較，如圖2所示.我們首先計算了單個波長內(nèi)單元剖分?jǐn)?shù)量為25的結(jié)果，在圖2中以黑色實(shí)線表示，進(jìn)而給出了單波長內(nèi)單元剖分?jǐn)?shù)量為10、5、3、2四種情況下同一道的振幅與其進(jìn)行比較，在圖2中以離散點(diǎn)表示.

由圖2可以看出，黑色實(shí)線所示的單波長內(nèi)單元剖分?jǐn)?shù)量為25的振幅曲線最為精確，橫向剖分精度為每波長10、5、3個單元時同一道的振幅與該計算結(jié)果非常吻合，進(jìn)一步說明單波長內(nèi)單元剖分?jǐn)?shù)為3時就已有足夠的計算精度.而當(dāng)剖分精度為每波長2個單元時，同一道振幅與該計算結(jié)果有較為明顯的差異，說明在此單元剖分?jǐn)?shù)情況下已難以保證計算精度.

圖3 一個含傾斜界面和非均質(zhì)的二維層狀模型Fig.3 A 2-D layered model with an oblique interface and heterogeneous layer

圖4 非均勻情況下不同縱向剖分精度時第51道振幅曲線比較(a) 第二層介質(zhì)速度擾動為5% ; (b) 第二層介質(zhì)速度擾動為25% .(a)、(b)中子圖自上至下分別代表單位波長內(nèi)薄層數(shù)量分別為10、5、3、2.Fig.4 51st amplitude curve comparison between different vertical discrete precision in heterogeneous situation(a) 5% perturbation of velocity for the second layer, (b) 25% perturbation of velocity for the second layer. Subgraphs from top to bottom in (a) and (b) stand for situations of 10,5,3,2 slabs per wavelength.

為了研究非均勻介質(zhì)情況下本文方法的模擬精度，我們對圖1所示模型中的第二層介質(zhì)加入兩種不同程度的波速擾動，使其體現(xiàn)不同程度的介質(zhì)非均勻性，如圖3所示.研究中分別考慮了介質(zhì)弱擾動和強(qiáng)擾動兩種情形，弱擾動對應(yīng)所示模型中第二層介質(zhì)5%的速度擾動，對于背景速度為3000 m·s-1的第二層介質(zhì)，加入5%的速度擾動意味著每個體積元的速度分布在2850 m·s-1至3150 m·s-1.類似地，強(qiáng)擾動對應(yīng)該層25%的速度擾動，即每個體積元的速度分布在2250 m·s-1至3750 m·s-1.我們分別計算了兩種情況下不同薄層離散厚度(即縱向離散精度)下第51道的單道記錄結(jié)果，而兩種情況下橫向離散精度在滿足采樣定理的前提下保持相同，同為單個波長下6個體積元.與均勻介質(zhì)精度研究類似，首先計算單個波長內(nèi)薄層剖分?jǐn)?shù)量為25的結(jié)果，在圖4中以黑色線表示，單波長內(nèi)薄層剖分?jǐn)?shù)為10、5、3、2四種情況在圖4中以離散點(diǎn)表示.需要特別說明的是，在本節(jié)中為了方便比較，我們將所有情況下的振幅曲線都進(jìn)行了歸一化處理(即振幅分布從0至1).圖4中給出了弱擾動和強(qiáng)擾動兩種情況下不同單波長內(nèi)薄層剖分?jǐn)?shù)量時的中間道振幅曲線的比較.

由圖4可知，當(dāng)以黑色實(shí)線所示的單波長內(nèi)薄層剖分?jǐn)?shù)量為25的振幅曲線作為比較對象時，在內(nèi)部介質(zhì)弱擾動情況下，縱向剖分精度滿足10、5、3個薄層每波長時，同一道的振幅與該計算結(jié)果仍非常吻合，而當(dāng)剖分精度僅滿足每波長2個薄層時，同一道振幅與該計算結(jié)果存在著一定的差異，說明介質(zhì)弱擾動情形接近于均勻介質(zhì)的計算結(jié)果.在介質(zhì)強(qiáng)擾動情況下，縱向剖分精度為5、3個薄層每波長時同一道的振幅與該計算結(jié)果已出現(xiàn)輕微的偏離，而當(dāng)剖分精度僅滿足每波長2個薄層時，同一道振幅與該計算結(jié)果差異非常大，表明介質(zhì)的強(qiáng)擾動對保證計算精度所需的最低單元剖分?jǐn)?shù)已產(chǎn)生明顯的影響，對于含有強(qiáng)擾動的非均勻介質(zhì)，保證計算精度需采用更大的單波長內(nèi)的薄層剖分?jǐn)?shù).

3 混合法實(shí)施

本節(jié)中我們針對一個包含著不規(guī)則邊界的層狀介質(zhì)，混合使用不規(guī)則界面全局廣義R/T遞推傳播矩陣方法和非均勻介質(zhì)全局廣義R/T遞推傳播矩陣方法進(jìn)行地震模擬，并與邊界元法(Fu, 2002)的計算結(jié)果進(jìn)行比較.邊界元方法在均質(zhì)情況下進(jìn)行地震波場計算時需要對整個模型的所有地層進(jìn)行邊界的剖分，分別計算邊界單元之間的相互影響，即格林函數(shù)，組裝成參數(shù)矩陣后統(tǒng)一進(jìn)行求解.而在非均質(zhì)情況進(jìn)行計算還需要對內(nèi)部介質(zhì)進(jìn)行剖分，再計算體積元與邊界單元以及體積元與體積元之間的相互影響之后統(tǒng)一組裝成矩陣，再進(jìn)行矩陣方程組的求解.

圖5 一個含不規(guī)則邊界的層狀均勻介質(zhì)Fig.5 A homogeneous layered media with irregular boundary

首先考慮內(nèi)部介質(zhì)均勻的情況，如圖5中的模型所示，該模型第一、二層界面為不規(guī)則界面，模型從上至下各層的速度分別為2700 m·s-1、3000 m·s-1、3600 m·s-1和4000 m·s-1，地震子波采用主頻為15 Hz的雷克子波，震源位于模型橫向距離1500 m，深度800 m處，101個檢波器以30 m的間隔均勻分布在地表.

圖6為對以上模型分別使用混合法和邊界元法得到的單炮記錄的比較圖，圖中黑線為混合法的計算結(jié)果，藍(lán)線為邊界元法的計算結(jié)果.從中可以看出兩種方法的計算結(jié)果基本一致.我們從兩種方法的計算結(jié)果中抽取第51道和第71道的記錄進(jìn)行波形和相位的比較，如圖7所示.從圖中可以看出，采用混合法和邊界元法的計算結(jié)果幾乎完全一致，表明了混合方法在處理非規(guī)則界面均勻介質(zhì)問題時的準(zhǔn)確性.

圖6 層狀均勻介質(zhì)混合法與邊界元法的單炮記錄比較(黑色：混合法，藍(lán)色：邊界元法)Fig.6 Single shot record comparison between hybrid method and boundary element method in layered homogeneous media (black: hybrid method, blue: boundary element method)

圖7 層狀均勻介質(zhì)混合法與邊界元法的第51道和第71道記錄振幅比較(上：第51道，下：第71道，黑線：混合法，點(diǎn)線：邊界元法)Fig.7 51st and 71st trace amplitude comparison between hybrid method and boundary element method in layered homogeneous media (top: 51st trace, bottom: 71st trace, black line: hybrid method, dot line: boundary element method)

圖8 一個含不規(guī)則的邊界的層狀非均勻介質(zhì)Fig.8 A heterogeneous layered media with irregular boundary

圖9 5%速度擾動下層狀非均勻介質(zhì)混合法與邊界元法的單炮記錄比較(黑色：混合法，藍(lán)色：邊界元法)Fig.9 Single shot record comparison between hybrid method and boundary element method in layered heterogeneous media with 5% velocity perturbation (black: hybrid method, blue: boundary element method)

圖10 25%速度擾動下層狀非均勻介質(zhì)混合法與邊界元法的單炮記錄比較(黑色：混合法，藍(lán)色：邊界元法)Fig.10 Single shot record comparison between hybrid method and boundary element method in layered heterogeneous media with 25% velocity perturbation (black: hybrid method, blue: boundary element method)

圖11 層狀非均勻介質(zhì)混合法與邊界元法的第51道記錄振幅比較(上：5%速度擾動，下：25%速度擾動，黑線：混合法，點(diǎn)線：邊界元法)Fig.11 51st trace amplitude comparison between hybrid method and boundary element method in layered homogeneous media (top: 5% velocity perturbation, bottom: 25% velocity perturbation, black line: hybrid method, dot line: boundary element method)

進(jìn)一步考慮介質(zhì)非均勻的情況，我們在圖5所示模型中震源位置的上方增加了一個200 m厚的水平非均質(zhì)層，如圖8所示.非均質(zhì)層分布在深度500 m至700 m，其背景速度仍為3600 m·s-1，在實(shí)際內(nèi)部介質(zhì)剖分中，首先按照非均質(zhì)層的體積元個數(shù)生成一個元素為-1到1的隨機(jī)數(shù)矩陣，乘以速度擾動(如：5%)，并利用此矩陣和背景速度生成實(shí)際速度場，得到每個體積元的實(shí)際速度.與第2節(jié)中的數(shù)值算例類似，我們在這里也分別考慮該非均質(zhì)層為弱擾動或強(qiáng)擾動兩種情況，在該層中分別加入5%和25%的速度擾動，并和邊界元法的計算結(jié)果進(jìn)行比較，如圖9和圖10所示，同樣針對弱、強(qiáng)兩種擾動情況，我們從中抽取了第51道的記錄進(jìn)行振幅比較，結(jié)果如圖11所示.從圖9的地震記錄比較圖中可以看出，在弱擾動情況下，混合法計算出的波場仍能和邊界元法保持較好的一致；而在圖10中所示強(qiáng)擾動情況下，混合法的計算結(jié)果已經(jīng)和邊界元法的計算結(jié)果有一定的差異，這是由于在強(qiáng)擾動情況下橫向速度變化變大，而混合法中的非均勻介質(zhì)全局廣義R/T遞推傳播矩陣方法基于一階Born近似，對于體積分的計算精度不夠所致.而從圖11中的同一道地震記錄比較可以看出，弱擾動的情況下，單道記錄的波形也基本相似，而對于強(qiáng)擾動情況，反射波的相位和振幅已經(jīng)和邊界元的計算結(jié)果發(fā)生了偏差，而多次波部分也存在相位和振幅差異.

4 結(jié)論

本文將常規(guī)的不規(guī)則邊界R/T遞推傳播矩陣方法進(jìn)行了拓展，建立了一種用于求解二維SH波情況下含不規(guī)則邊界、多層非均勻介質(zhì)波傳播問題的全局廣義R/T遞推傳播矩陣混合方法.這一方法將不規(guī)則邊界、多層非均勻介質(zhì)內(nèi)傳播的地震波分解為入射波、邊界散射波和體散射波三個部分，總波場采用廣義Lippmann-Schwinger積分方程進(jìn)行精確疊加求解.該方法的關(guān)鍵在于將不規(guī)則邊界與介質(zhì)的非均勻性進(jìn)行了解耦，對半空間內(nèi)不規(guī)則邊界產(chǎn)生的散射波，通過將半空間劃分為一系列水平非均質(zhì)薄層，計算含不規(guī)則邊界的薄層的傳播算子來求解；對于考慮介質(zhì)非均勻性的一般的沉積地層，直接采用邊界積分替代體積分，即對每個薄層進(jìn)行一階Born近似并獲得傳播算子求解，而對于速度橫向變化大的薄板(含鹽丘)可采用二階或高階Born近似來提高模擬精度.需要說明的是，對于介質(zhì)速度空間平緩變化的情況，我們沿用了前人的方法(Fu，2002；Fu and Bouchon, 2004)，將非均勻介質(zhì)波動方程簡化為均勻背景介質(zhì)波動方程，該簡化可直接將半空間均勻介質(zhì)的格林函數(shù)作為波動方程的解.采用本文發(fā)展的混合方法，模型中任一位置的地震波場均可通過震源項(xiàng)和整個半空間的廣義傳播矩陣經(jīng)過迭代實(shí)現(xiàn)求解.本文重點(diǎn)研究了非均質(zhì)地層薄板化R/T傳播矩陣法的模擬精度，包括與波長相關(guān)的縱/橫向離散精度以及與非均質(zhì)強(qiáng)度相關(guān)的模擬精度，并與邊界元的計算結(jié)果進(jìn)行了比較.對同時含有不規(guī)則邊界地層和隨機(jī)非均質(zhì)地層的復(fù)雜模型采用混合方法進(jìn)行地震波傳播數(shù)值模擬，結(jié)果表明本文提出的混合法是一種解決復(fù)雜模型高精度地震模擬的有效方法.

作為一種可以替代有限差分或有限元的高效數(shù)值計算方法，本文提出的全局廣義R/T遞推傳播矩陣混合法在計算速度和精度上都有很大優(yōu)勢.該方法根據(jù)采樣定理進(jìn)行薄板離散，其極限精度取決于薄板離散的個數(shù)，隨著薄板(層)個數(shù)的增加，計算量也會相應(yīng)增大.然而該方法提供了一種靈活的選擇，即在廣域均質(zhì)區(qū)域可以通過增大薄板厚度以提高計算速度，在非均質(zhì)區(qū)域則可通過減小薄板厚度以提高計算精度.這樣的變剖分精度的方法尤其適用于含非均勻?qū)拥膶訝罱橘|(zhì)模型數(shù)值求解，在計算精度、效率上相比有限元、有限差分方法有明顯的優(yōu)勢，在未來地震學(xué)相關(guān)研究中用于處理不規(guī)則邊界和非均勻介質(zhì)波場計算有廣泛的應(yīng)用前景.需要指出的是本文只考慮SH波情形，將本文方法用于P-SV波波場求解尚待進(jìn)一步的研究工作.

附錄

本文方法的推導(dǎo)過程分為如下幾個步驟：第一步，通過整個模型的物理性質(zhì)，獲得廣義Lippmann-Schwinger積分方程表達(dá)式，進(jìn)而轉(zhuǎn)換到頻率域下并進(jìn)行橫向離散與薄層離散；第二步，使用Waterman(1975)的波場展開方法將積分方程組轉(zhuǎn)換成關(guān)于位移和應(yīng)力的矩陣方程組；第三步，將矩陣方程組中的參數(shù)矩陣轉(zhuǎn)換成修正R/T矩陣進(jìn)而轉(zhuǎn)換成全局廣義R/T矩陣，從震源層出發(fā)，利用震源層的上下行波和全局廣義R/T矩陣實(shí)現(xiàn)波場中任一點(diǎn)的波場求解.以上三步的詳細(xì)推導(dǎo)過程如下：

A 頻率域下廣義Lippmann-Schwinger積分方程的推導(dǎo)及離散

圖A1中為一散射介質(zhì)模型，假設(shè)震源S(ω)位于r0處，則此散射介質(zhì)中的地震響應(yīng)滿足以下標(biāo)量Helmholtz波動方程：

(A1)

圖A1 二維層狀非均勻介質(zhì)半空間模型示意圖Fig.A1 A 2-D heterogeneous layered media in half space

其中ω為角頻率，k=ω/v(r)為波數(shù).

根據(jù)(A1)式，在r處的地震響應(yīng)u(r)可表示為入射波場u(0)(r)、邊界散射u(1)(r)和體積散射u(2)(r)三部分，其表達(dá)式分別為：

(A2)

(A3)

(A4)

綜合式(A2)至(A4)，可以得到廣義Lippmann-Schwinger積分方程，即：

(A5)

其中C(r)為與邊界幾何特征相關(guān)的參數(shù).

為了求解以上積分方程，我們使用離散波數(shù)法對公式(A5)中的邊界積分項(xiàng)和體積分項(xiàng)進(jìn)行離散.假設(shè)介質(zhì)具有局部非均勻性，介質(zhì)按長度L在水平方向上按周期性分布n次，則任意層內(nèi)的格林函數(shù)的離散波數(shù)形式為：

(A6)

(A7)

式中相應(yīng)參數(shù)的表達(dá)式為：

(A8)

(A9)

(A10)

特別地，在求解式(A10)中的體積分項(xiàng)時采用梯形法則來進(jìn)行簡化，將非均質(zhì)層劃分成多個薄層并采用上下邊界的邊界積分來近似代替體積分，相當(dāng)于對每一個薄層做一次一階Born近似，近似后得到的表達(dá)式為：

(A11)

圖A2 對不規(guī)則邊界與介質(zhì)非均勻性進(jìn)行區(qū)分的分層模型Fig.A2 A layered model which separate the irregular boundary and interior heterogeneity

至此，我們就可以針對圖A1中的原始模型進(jìn)行薄層剖分，并對薄層類型進(jìn)行分類，如圖A2所示.

舉例來說，對圖A2所示的分層模型，若由Γ1與Γ2以及兩側(cè)的人工邊界圍成的薄層中包含著不規(guī)則邊界，則使用不規(guī)則界面全局廣義R/T遞推傳播矩陣方法來進(jìn)行計算；若由Γ2與Γ3以及兩側(cè)的人工邊界圍成的薄層中包含著非均勻介質(zhì)，則使用非均勻介質(zhì)全局廣義R/T遞推傳播矩陣方法進(jìn)行計算.其中，含不規(guī)則邊界的層中的廣義Lippmann-Schwinger積分方程僅含有a(j)(x,z)項(xiàng)，而含非均質(zhì)薄層的方程則同時包含a(j)(x,z)和b(j)(x,z)項(xiàng).我們將每個薄層對應(yīng)的方程組分成兩類，分別使用適合該薄層類型的方法進(jìn)行混合求解.

B 矩陣方程組的轉(zhuǎn)化

根據(jù)式(A7)中的表達(dá)式，提取其中的位移項(xiàng)α與應(yīng)力項(xiàng)β，將與位移項(xiàng)和應(yīng)力項(xiàng)相乘的參數(shù)歸為三類：(1)與位移項(xiàng)相乘的邊界參數(shù)Q；(2)與應(yīng)力項(xiàng)相乘的邊界參數(shù)P；(3)與位移項(xiàng)相乘的內(nèi)部介質(zhì)參數(shù)M.根據(jù)Waterman(1975)的波場展開，并分別考慮被觀測點(diǎn)相對于薄層的相對位置，式(A7)可以分為如下兩種情況：

① 當(dāng)z<ξ(j-1)且δsj=0時：

(B1)

② 當(dāng)z>ξ(j)且δsj=0時：

(B2)

(B3)

(B4)

(B5)

(B6)

(B7)

(B8)

(B9)

(B10)

(B11)

(B12)

(B13)

(B14)

(B15)

(B16)

(B17)

(B18)

(B19)

(B20)

其中，Δξ(j)(x′)為表述界面起伏的參數(shù).需要注意的是，僅在處理不規(guī)則界面部分時才需要計算公式(B17)至(B20)部分，在處理非均質(zhì)部分時，由于界面起伏為0，公式(B3)至(B14)中的積分項(xiàng)的值都為1，僅需要計算積分項(xiàng)前的系數(shù).

C 全局修正反射/透射矩陣及全局廣義反射/透射矩陣的推導(dǎo)

(C1)

同時，層內(nèi)的上下行波和修正反射/透射矩陣滿足以下條件：

(C2)

(C3)

①當(dāng)薄層位于震源層之上時，有：

(C4)

②當(dāng)薄層位于震源層之下時，有：

(C5)

至此，模型中任意一點(diǎn)的波場可以通過從震源層出發(fā)的上下行波與全局廣義R/T矩陣通過迭代關(guān)系獲得.