水力壓裂微地震震源定位與震源機(jī)制聯(lián)合反演研究

李　宏,楊心超,朱海波,張　偉,曲壽利

(1.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;2.南方科技大學(xué),廣東深圳518055)

水力壓裂改造技術(shù)已成為頁巖氣和致密砂巖氣等非常規(guī)低滲透性儲層有效的增產(chǎn)技術(shù),對壓裂過程中產(chǎn)生的微地震活動進(jìn)行監(jiān)測,是評價壓裂改造效果的一種重要手段[1-2]。微地震監(jiān)測主要有井中監(jiān)測與地面監(jiān)測兩種方式。井中監(jiān)測方式采集的數(shù)據(jù)質(zhì)量較高,但受檢波器數(shù)量與觀測范圍限制,對震源參數(shù)(震源位置與機(jī)制)反演約束較差,而且如果在壓裂井附近缺少觀測井的情況下則無法進(jìn)行監(jiān)測[3];地面監(jiān)測方式觀測范圍更大,覆蓋次數(shù)更高,因而可以避免此類問題的出現(xiàn),但其主要問題是受地表環(huán)境因素影響,資料品質(zhì)較差、信噪比較低[4-7]。

傳統(tǒng)的地面微地震監(jiān)測資料的處理通常分為兩個步驟:先進(jìn)行震源定位,然后利用定位信息進(jìn)行震源機(jī)制反演[8]。因?yàn)榈孛嫖⒌卣鸨O(jiān)測資料信噪比較低,初至拾取較為困難[9-10],所以格點(diǎn)搜索類方法成為震源定位反演的一種常用方法[11]。該方法對震源位置空間進(jìn)行格點(diǎn)離散并搜索,目標(biāo)函數(shù)采用疊加函數(shù)的形式,疊加能量最大的格點(diǎn)即為震源位置。實(shí)際資料顯示微地震震源機(jī)制有很強(qiáng)的雙力偶震源成分,其振幅極性在地表呈四象限分布[12]。在進(jìn)行震源定位時,目標(biāo)函數(shù)若采用直接疊加的方式,振幅極性的反轉(zhuǎn)會造成微震疊加過程中能量的抵消。采用絕對值疊加函數(shù)可避免這個問題,但噪聲水平會同步放大,并不能很好地突出微震事件的疊加能量。疊加函數(shù)的選擇會直接影響震源定位的精度,在定位偏差較大的情況下進(jìn)行震源機(jī)制反演,會將定位誤差帶入到震源機(jī)制反演中,造成震源機(jī)制反演誤差較大,這對后續(xù)的微地震監(jiān)測資料解釋造成困擾。

目前主要有兩種方法解決振幅極性反轉(zhuǎn)帶來的問題。一種是基于互相關(guān)的地震干涉偏移類方法[13]。ZHEBEL等[14-15]通過對相鄰道做互相關(guān)處理來消除振幅極性反轉(zhuǎn)的影響,從而實(shí)現(xiàn)了震源定位與震源機(jī)制聯(lián)合反演;LI等[16]基于該方法提出了加權(quán)彈性波干涉成像方法,并應(yīng)用于多波多分量記錄的微地震定位。另外一種方法是通過震源定位與震源機(jī)制聯(lián)合反演,利用震源機(jī)制來校正微地震記錄的極性。LIANG等[17]采用格點(diǎn)掃描疊加方法,對離散震源位置與震源機(jī)制組成的全空間進(jìn)行離散掃描,并用震源機(jī)制校正記錄極性。ANIKIEV等[18]在每個搜索格點(diǎn)進(jìn)行震源機(jī)制反演來獲得震源機(jī)制,并用此震源機(jī)制來計算極性校正后的疊加函數(shù)。

本文介紹了一種格點(diǎn)掃描的震源定位與震源機(jī)制聯(lián)合反演的迭代搜索算法。首先根據(jù)常規(guī)震源定位方法獲得震源初始位置,采用P波初動震源機(jī)制反演獲得震源機(jī)制,用該震源機(jī)制校正微地震記錄并進(jìn)行定位,迭代此過程以獲得較準(zhǔn)確的震源位置,然后采用P波輻射花樣約束的震源機(jī)制反演方法,以獲得較準(zhǔn)確的震源機(jī)制。最后給出了模型測試結(jié)果,驗(yàn)證該方法的有效性和優(yōu)越性;并將該方法應(yīng)用于中國四川盆地某壓裂井的地面微地震監(jiān)測實(shí)際資料的處理。

1　方法原理

傳統(tǒng)的地面微地震監(jiān)測處理流程是先進(jìn)行震源定位再進(jìn)行震源機(jī)制反演,這種方式會造成誤差的傳遞;而對由震源位置與震源機(jī)制組成的震源參數(shù)全空間進(jìn)行格點(diǎn)搜索的聯(lián)合反演方式,由于搜索樣本數(shù)目太大,對計算能力有較高的要求[17]。

因此本文提出了一種震源位置與震源機(jī)制迭代搜索方法,技術(shù)流程如圖1所示。采用這種方法既可以避免震源參數(shù)全空間同時離散造成的搜索格點(diǎn)急劇增長,也可以疊加目標(biāo)函數(shù)振幅校正使其更具有實(shí)際物理意義,并且根據(jù)迭代過程中對震源機(jī)制精度需求的不同,本文采用了兩種震源機(jī)制反演方法,進(jìn)一步提高反演速度。

我們首先采用常規(guī)的定位算法對微地震事件進(jìn)行初始定位,獲得震源初始位置后,進(jìn)行震源定位和震源機(jī)制迭代反演,具體如下:首先進(jìn)行P波初動極性反演獲得震源機(jī)制,然后用此震源機(jī)制來校正微地震記錄,并采用極性校正疊加函數(shù)重新進(jìn)行定位,最后迭代計算直至震源位置收斂,此時可以獲得較為準(zhǔn)確的震源位置與初步的震源機(jī)制結(jié)果。

由于地面觀測系統(tǒng)覆蓋區(qū)域有限,可能無法對某一震源機(jī)制節(jié)面位置進(jìn)行采樣,因而P波初動極性震源機(jī)制反演存在多解性問題[19]。但即使得不到準(zhǔn)確的震源機(jī)制解,也可以提供正確的檢波點(diǎn)初動極性,故不會影響定位結(jié)果。在獲得準(zhǔn)確的定位結(jié)果后,我們可進(jìn)一步采用更精確震源機(jī)制反演方法來反演震源機(jī)制,本文采用的是輻射花樣約束的震源機(jī)制反演方法。

圖1　迭代搜索算法流程

在采用掃描疊加搜索方法進(jìn)行震源定位時,疊加目標(biāo)函數(shù)一般有以下3種方式:

(1)

(2)

(3)

式中:ui為第i道的記錄;τi,x為震源在x點(diǎn)到第i檢波點(diǎn)的旅行時;si為第i道P波初動極性;f1(t)為疊加函數(shù)的絕對值;f2(t)為絕對值疊加函數(shù);f3(t)為極性校正疊加函數(shù)。

假設(shè)微地震震源機(jī)制主要成分為雙力偶型震源,則其振幅在地表呈四象限分布。若采用公式(1)直接疊加方式,則導(dǎo)致正負(fù)區(qū)域振幅能量相互抵消。甚至當(dāng)震源傾角為近90°,觀測系統(tǒng)在震源上方對稱分布時,即使事件震級很大,但因f1(t)值為零,故無法對其進(jìn)行準(zhǔn)確定位。如若目標(biāo)函數(shù)采用公式(2)絕對值疊加方式,雖然可以避免出現(xiàn)振幅抵消的情況,但該函數(shù)會同時提高噪聲水平,不利于目標(biāo)函數(shù)極值拾取。本文在迭代反演過程中采用的是公式(3)極性校正疊加:利用震源機(jī)制得到的極性來校正微地震記錄,然后進(jìn)行疊加。這種方法可以避免采用公式(1)與公式(2)帶來的問題,以獲得較好的疊加效果[17-18]。

在整個反演流程中用到了兩種震源機(jī)制反演方法。首先在震源定位和震源機(jī)制迭代反演過程中采用了P波初動極性震源機(jī)制反演算法[19],這種算法只利用了P波初動極性信息,故計算速度較快,而且可以提供一個大致準(zhǔn)確的震源機(jī)制,以滿足定位時振幅校正的需要。該算法對震源坐標(biāo)系在地平坐標(biāo)系中的微震進(jìn)行全空間掃描,計算矛盾極性比ψ:

(4)

式中:NC為矛盾極性數(shù),NT為極性總數(shù)。當(dāng)矛盾極

性比ψ最小時,P波初動極性與觀測極性擬合最佳,即為求取的一組震源機(jī)制節(jié)面。

獲得準(zhǔn)確震源定位后,為了避免P波初動極性震源機(jī)制反演的多解性,獲得更為精確的震源機(jī)制,本文進(jìn)一步采用P波輻射花樣約束下的震源機(jī)制反演算法[20]。該算法通過對震源機(jī)制全空間掃描,求取與實(shí)際P波輻射花樣吻合度最好的一組震源機(jī)制節(jié)面,采用的反演目標(biāo)函數(shù)是實(shí)際觀測和理論計算之間的相關(guān)系數(shù)R(吻合度)[6]:

(5)

2　模型測試

為了驗(yàn)證本文方法的可靠性,我們利用理論模型數(shù)據(jù)進(jìn)行了測試。由于受上、下地層與圍壓因素影響,目的層埋深超過1000m時,水力壓力裂縫呈高角度裂縫[21],所以設(shè)計微地震震源模型采用的裂縫參數(shù):方位、傾角、滑動角分別為20°,90°,40°,位于(0,0,1500m)處。速度模型和觀測系統(tǒng)如圖2所示,在深度為1000m處有一高速層;檢波陣列呈米字狀分布,共320道檢波器,檢波距為50m,最大偏移距為2000m,最大偏移距與震源深度比約為1.3,與實(shí)際水力壓裂地面微地震監(jiān)測觀測系統(tǒng)分布相似。本文采用CHEN等[22]的廣義反射-透射系數(shù)正演模擬算法生成理論地震記錄,如圖3所示。

圖2　微地震震源模型的速度模型(a)與觀測系統(tǒng)(b)

從圖3中拾取初至P波同相軸并截取初至P波波形,通過與理論子波互相關(guān),計算各道初至P波的初動極性和振幅,利用本文提出的聯(lián)合反演方法反演震源位置與震源機(jī)制參數(shù)。其定位結(jié)果為(0,0,1500m),震源機(jī)制反演結(jié)果為(21°,90°,40°),如表1中所示的信噪比為∞時的結(jié)果。而采用傳統(tǒng)反演方法,得到的位置為(0,0,1508m),震源機(jī)制為(21°,90°,40°),如表2 中所示的信噪比為∞時的結(jié)果?？梢钥闯?在沒有噪聲的情況下,兩種方法都得到較準(zhǔn)確的結(jié)果,但聯(lián)合反演方法與傳統(tǒng)方法相比,垂直方向上定位精度更高。圖4給出了模型理論P(yáng)波振幅與聯(lián)合反演得到的振幅對比結(jié)果,圖中圓形表示理論振幅,方形表示反演結(jié)果,圖標(biāo)大小表示振幅相對大小,紅色和藍(lán)色分別表示正負(fù)極性。可以看出,理論振幅與反演振幅無論極性還是振幅相對大小分布都非常吻合,驗(yàn)證了本文提出方法的有效性。

圖3　微地震震源模型(20°,90°,40°)的正演模擬微地震記錄(地面垂直分量)

圖4　微地震震源模型理論P(yáng)波振幅與反演結(jié)果振幅對比

為了驗(yàn)證本文方法的抗噪性,我們從實(shí)際野外地面微地震監(jiān)測記錄中截取一段背景噪聲(圖5),并與圖3所示的理論地震記錄合成不同噪聲水平的地震記錄,并對此記錄進(jìn)行反演測試。表1給出了本文提出的聯(lián)合反演方法的反演結(jié)果,表2給出了傳統(tǒng)反演方法的反演結(jié)果。圖6對比了兩種方法定位結(jié)果的誤差,其中藍(lán)線為聯(lián)合反演的定位誤差,紅線為傳統(tǒng)方法的定位誤差,可以看出,相比于傳統(tǒng)反演方法,本文的聯(lián)合反演方法定位精度更高,且抗噪性能更好。對比表1與表2可以發(fā)現(xiàn),聯(lián)合反演方法在垂直方向上具有更高的精度,說明該方法可以提高垂直方向上的反演約束。圖7對比了不同信噪比條件下利用兩種反演方法進(jìn)行震源機(jī)制反演的結(jié)果誤差,圖7a為傳統(tǒng)反演方法震源機(jī)制反演誤差隨信噪比變化結(jié)果,圖7b為聯(lián)合反演方法震源機(jī)制反演誤差隨信噪比變化結(jié)果?？梢钥吹?無論反演精度還是穩(wěn)定性,聯(lián)合反演方法都要比傳統(tǒng)方法好。

傳統(tǒng)反演方法將震源定位與震源機(jī)制反演分開進(jìn)行,并且在震源定位時采用絕對值疊加函數(shù),使得幅值正負(fù)分布的噪聲變成幅值為正的噪聲。隨著噪聲水平的增加目標(biāo)函數(shù)信噪比逐步降低,計算目標(biāo)函數(shù)極值誤差逐步增加,所以其定位誤差也逐步增加。在定位誤差增大的情況下,再進(jìn)行震源機(jī)制反演,必然會將定位誤差引入,導(dǎo)致震源機(jī)制反演誤差增大。而在聯(lián)合反演方法中,定位時采用的是極性校正疊加函數(shù),而不是絕對值疊加,這樣疊加時噪聲仍呈正負(fù)隨機(jī)分布,通過疊加可以有效壓制噪聲,使得定位精度的提高,并有助于震源機(jī)制反演精度的提高。

另外本文提出的聯(lián)合反演方法迭代收斂非常迅速,一般2到3個循環(huán)即可收斂。原因在于迭代過程中只要反演獲得的P波初動極性與觀測極性擬合較好,其定位結(jié)果會迅速收斂。P波初動極性震源機(jī)制反演的主要目的是為震源定位提供一個準(zhǔn)確的極性校正函數(shù),而該反演算法也可以快速完成。所以聯(lián)合反演方法所用的運(yùn)算時間大概是傳統(tǒng)反演方法2～3倍,而與全空間離散掃描方法相比,反演速度提高了10倍以上。

圖5　實(shí)際野外地面微地震監(jiān)測記錄中截取的噪聲

信噪比定位結(jié)果/mxyz誤差/m震源機(jī)制反演結(jié)果/(°)方位角傾角滑動角吻合度∞001500021.090.040.00.9964100014928.021.090.040.00.996680014928.021.090.040.00.99596412149215.021.788.140.00.992744815049.924.186.041.80.98702-44152020.827.684.537.70.9576

表2　理論模型傳統(tǒng)反演方法結(jié)果

圖6　不同信噪比條件下傳統(tǒng)反演方法與聯(lián)合反演方法定位誤差對比

圖7　不同信噪比條件下傳統(tǒng)反演方法(a)與聯(lián)合反演方法(b)震源機(jī)制反演結(jié)果誤差對比

3　實(shí)際資料應(yīng)用

采用本文提出的聯(lián)合反演方法對四川盆地某壓裂井地面微地震實(shí)際監(jiān)測資料進(jìn)行震源定位與震源機(jī)制聯(lián)合反演研究。利用聲波測井曲線建立速度模型如圖8a 所示,壓裂監(jiān)測所采用的觀測系統(tǒng)如圖8b所示,射孔段深度約為3294～3369m。共10條測線,呈放射狀分布,相鄰測線方位角相差約36°,每條測線有125個檢波點(diǎn),道間距為25m,最小偏移距為300m,最大偏移距為3400m。本次監(jiān)測共識別了1073個微地震事件,從中選擇109個較強(qiáng)微地震事件進(jìn)行處理分析。

圖9給出了一個典型微地震事件記錄,表3給出了對該微地震事件分別采用傳統(tǒng)反演方法和本文提出的反演方法的反演結(jié)果?？梢钥闯?利用聯(lián)合反演方法得到的震源機(jī)制的反演結(jié)果吻合度更高,可以認(rèn)為反演結(jié)果更為可靠。圖10對比了聯(lián)合反演結(jié)果振幅與實(shí)際振幅,可以看出極性、振幅大小以及分布吻合得很好。圖11顯示了所處理的109個微地震事件的震源機(jī)制吻合度,平均值為0.9452,與傳統(tǒng)方法得到的震源機(jī)制吻合度平均值0.9380[6]相比有所提升,可以認(rèn)為本次反演結(jié)果精度整體有所提升。

圖12對比了聯(lián)合反演方法與傳統(tǒng)反演方法震源定位結(jié)果,結(jié)合圖13顯示的由聯(lián)合反演方法得到的震源位置與震源機(jī)制結(jié)果,可以看出聯(lián)合反演方法得到的震源節(jié)面連續(xù)性更好,水力壓裂產(chǎn)生的縫網(wǎng)發(fā)育方向以北偏東方向?yàn)橹?可以認(rèn)為其更能準(zhǔn)確反映裂縫的延展走向。圖14給出的震源機(jī)制分布顯示,壓裂裂縫面的延展方向也是北偏東方向(10°～20°),兩者較為一致。

圖8　四川盆地某壓裂井微震速度模型(a)與觀測系統(tǒng)(b)

圖9　實(shí)際資料中某一典型微地震記錄

反演方法震源位置/m震源機(jī)制/(°)xyz方位角傾角滑動角吻合度傳統(tǒng)反演36963860323312.975.085.90.9618聯(lián)合反演36963820321311.472.086.80.9738

圖10　某一典型微震事件聯(lián)合反演振幅與實(shí)際振幅對比

圖11　聯(lián)合反演中P波振幅約束震源機(jī)制反演目標(biāo)函數(shù)值

圖12　聯(lián)合反演方法震源定位結(jié)果與傳統(tǒng)反演方法震源定位結(jié)果對比

圖13　聯(lián)合反演方法得到的震源位置與震源機(jī)制結(jié)果顯示

圖14　聯(lián)合反演方法得到的震源機(jī)制(方位角、傾角、滑動角)分布a 高角度震源機(jī)制下的方位角; b 高角度震源機(jī)制下的傾角; c 高角度震源機(jī)制下的滑動角; d 低角度震源機(jī)制下的方位角; e 低角度震源機(jī)制下的傾角;f 低角度震源機(jī)制下的滑動角

4　結(jié)論

本文介紹了一種基于格點(diǎn)掃描方式的震源定位與震源機(jī)制聯(lián)合反演的迭代搜索算法。與傳統(tǒng)震源參數(shù)反演方法相比,該方法采用定位結(jié)果與震源機(jī)制結(jié)果互相約束,可以同時提高震源位置與震源機(jī)制的反演精度;并且采用迭代搜索方式,避免了同時搜索震源位置與震源機(jī)制空間造成的計算量過大問題,提高了反演速度。

如果在進(jìn)行格點(diǎn)搜索時考慮采用多尺度搜索策略,可以進(jìn)一步降低搜索格點(diǎn)數(shù)目,提高計算效率。另外若將所有格點(diǎn)旅行時預(yù)先計算好并存儲起來,可避免迭代過程中重復(fù)計算,搜索效率將進(jìn)一步提升。

在獲得震源位置與震源機(jī)制信息后,可以進(jìn)一步結(jié)合地質(zhì)資料、地層應(yīng)力背景以及巖心信息來選擇具體裂縫破裂面,用于裂縫空間幾何展布和縫網(wǎng)演化過程等壓裂解釋工作。

[1]MAXWELL S.Microseismic:Growth born from success[J].The Leading Edge,2010,29(3):338-343

[2]朱海波,楊心超,廖如剛,等.基于微地震裂縫參數(shù)反演的解釋與應(yīng)用研究[J].石油物探,2017,56(1):150-157

ZHU H B,YANG X C,LIAO R G,et al.Microseismic fracture interpretation and application based on parameters inversion[J].Geophysical Prospecting for Petroleum,2017,56(1):150-157

[3]MAXWELL S C,RUTLEDGE J,JONES R,et al.Petroleum reservoir characterization using downhole microseismic monitoring[J].Geophysics,2010,75(5):A75-A129

[4]DUNCAN P,EISNER L.Reservoir characterization using surface microseismic monitoring[J].Geophysics,2010,75(5):A139-A146

[5]宋維琪,喻志超,楊勤勇,等.低信噪比微地震事件初至拾取方法研究[J].石油物探,2013,52(6):596-601

SONG W Q,YU Z C,YANG Q Y,et al.The first arrival picking method of microseismic events with low S/N[J].Geophysical Prospecting for Petroleum,2013,52(6):596-601

[6]楊心超,朱海波,李宏,等.基于P波輻射花樣的壓裂微地震震源機(jī)制反演方法研究及應(yīng)用[J].石油物探,2016,55(5):640-648

YANG X C,ZHU H B,LI H,et al.Microseismic focal mechanism inversion based on P-wave radiation pattern and its application[J].Geophysical Prospecting for Petroleum,2016,55(5):640-648

[7]程磊磊,姜宇東,崔樹果,等.基于強(qiáng)事件約束的微地震剩余靜校正量估算方法[J].石油物探,2015,54(6):690-698

CHENG L L,JIANG Y D,CUI S G,et al.The estimation of microseismic residual statics with the constraint of strong-event smoothness[J].Geophysical Prospecting for Petroleum,2015,54(6):690-698

[8]朱海波,楊心超,王瑜,等.水力壓裂微地震監(jiān)測的震源機(jī)制反演方法應(yīng)用研究[J].石油物探,2014,53(5):556-561

Zhu H B,Yang X C,Wang Y,et el.The application of microseismic source mechanism inversion in hydraulic fracturing monitoring[J].Geophysical Prospecting for Petroleum,2014,53(5):556-561

[9]譚玉陽,何川,曹耐.基于多道相似系數(shù)的微地震事件自動識別[J].石油物探,2015,54(2):126-132

TAN Y Y,HE C,CAO N.Automatic microseismic event detection based on multi-channel semblance coefficient[J].Geophysical Prospecting for Petroleum,2015,54(2):126-132

[10]盛冠群,李振春,王維波,等.基于小波分解與高階統(tǒng)計量的微地震初至拾取方法研究[J].石油物探,2015,54(4):388-395

SHENG G Q,LI Z C,WANG W B,et al.A new automatic detection method of microseismic events based on wavelet decomposition and high-order statistics[J].Geophysical Prospecting for Petroleum,2015,54(4):388-395

[11]KAO H,SHAN S J.Rapid identification of earthquake rupture plane using Source-Scanning Algorithm[J].Geophysical Journal International,2007,168(3):1011-1020

[12]BAIG A,URBANCIC T.Microseismic moment tensors:a path to understanding frac growth[J].The Leading Edge,2010,29(3):320-324

[13]SCHUSTER G T,Yu J,Sheng J,et al.Interferometric/daylight seismic imaging[J].Geophysical Journal International,2004,157(2):838-852

[14]ZHEBEL O,EISNER L.Simultanous microseismic event localization and source mechanism determination[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1033-1037

[15]ZHEBEL O,EISNER L.Simultaneous microseismic event localization and source mechanism determination[J].Geophysics,2015,80(1):KS1-KS9

[16]LI L,CHEN H,WANG X M.Weighted elastic wave interferometric imaging of microseismic source location[J].Applied Geophysics,2015,12(2):221-234

[17]LIANG C,YU Y,YANG Y,et al.Joint inversion of source location and focal mechanism of microseismicity[J].Geophysics,2016,81(2):KS103-KS111

[18]ANIKIEV D,VALENTA J,STANK F,et al.Joint location and source mechanism inversion of microseismic events:benchmarking on seismicity induced by hydraulic fracturing[J].Geophysical Journal International,2014,198(1):249-258

[19]許忠淮,閻明,趙仲和.由多個小地震推斷的華北地區(qū)構(gòu)造應(yīng)力場的方向[J].地震學(xué)報,1983,5(3):15-26

XU Z H,YAN M,ZHAO Z H.Evaluation of the direction of tectonic stress in north China from recorded data of a large number of small earthquakes[J].Acta Seismologica Sinica,1983,5(3):15-26

[20]楊心超,朱海波,崔樹果,等.P波初動震源機(jī)制解在水力壓裂微地震監(jiān)測中的應(yīng)用[J].石油物探,2015,54(1):43-50

YANG X C,ZHU H B,CUI S G,et al.Application of P-wave first-motion focal mechanism solutions in microseismic monitoring for hydraulic fracturing[J].Geophysical Prospecting for Petroleum,2015,54(1):43-50

[21]FISHER M K,WARPINSKI N R.Hydraulic-fracture-height growth:real data[J].SPE Production & Operations,2012,27(1):8-19

[22]CHEN X.Seismogram synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method.I.Theory of two-dimensional SH case[J].Bulletin of the Seismological Society of America,1990,80(6A):1696-1724