田 坤，黃建平，李振春，李 娜，孔 雪，李慶洋

1.中國(guó)石油大學(xué)地球科學(xué)與技術(shù)學(xué)院，山東 青島 266580

2.中石油東方地球物理公司研究院，河北 涿州 072750

0 引言

由于計(jì)算能力的快速發(fā)展，在過(guò)去幾十年里對(duì)地震波在復(fù)雜介質(zhì)中傳播數(shù)值模擬方法的研究也越來(lái)越廣泛和深入，其中應(yīng)用最廣泛的是有限差分法［1］。在有限差分算法中，計(jì)算模型不可能無(wú)限大，也就是模擬區(qū)域都是有界的，這樣就需要在人工邊界處對(duì)波進(jìn)行吸收衰減。吸收邊界條件有許多種，其中比較有效的是Bérenger［2］在電磁波模擬中提出的完全匹配層（perfectly matched layer，PML）吸收邊界條件。其后Chew等［3］引入復(fù)數(shù)伸展坐標(biāo)系對(duì)PML吸收邊界條件進(jìn)行公式化，指出PML的本質(zhì)是在復(fù)數(shù)伸展坐標(biāo)系中進(jìn)行坐標(biāo)變換。Rappaport［4］證明了PML介質(zhì)等價(jià)于在吸收邊界區(qū)域引入各向異性介質(zhì)。近年來(lái)，PML吸收邊界條件也被應(yīng)用到聲波和彈性波的有限差分正演模擬中［5－7］。大量研究［3，8］表明，PML 吸收邊界條件比旁軸近似吸收邊界條件［9］、Higdon吸收邊界條件［10］、廖氏吸收邊界條件［11］和指數(shù)衰減吸收邊界條件［12］等具有更優(yōu)越的吸收性能。此外，Teixeira等［13］和Chew等［14］在柱坐標(biāo)系和球坐標(biāo)系中實(shí)現(xiàn)了PML吸收邊界條件。在國(guó)內(nèi)，研究者對(duì)PML吸收邊界條件也做了一些討論［15－17］。

但是，傳統(tǒng)的PML也存在一定的缺陷，離散后的PML反射系數(shù)不嚴(yán)格為0，尤其在大入射角的情況下更為嚴(yán)重；在這種情況下，相當(dāng)一部分能量以反射波的形式被送回到主研究區(qū)域［18－20］。而大入射角的情況是普遍存在的，比如薄片區(qū)域、震源位置接近研究區(qū)域邊緣、大偏移距接收等等。傳統(tǒng)PML對(duì)窄區(qū)域的面波吸收也有一定的問(wèn)題。為了克服這些問(wèn)題，Kuzuoglu等［21］提出了復(fù)頻移拉伸函數(shù)以改進(jìn)坐標(biāo)變換關(guān)系，在復(fù)頻移拉伸函數(shù)中引入了額外的2個(gè)衰減參數(shù)，從而改善PML的吸收性能，稱之為復(fù)頻移PML（CFS－PML）。采用復(fù)頻移拉伸函數(shù)后，PML算法不易用傳統(tǒng)的分裂變量的形式實(shí)現(xiàn)，而不分裂變量計(jì)算卷積又大大增加計(jì)算成本。Drossaert等［22］提出了基于遞歸積分的不分裂復(fù)頻移PML邊界條件，實(shí)際上就是采用復(fù)頻移拉伸函數(shù)，通過(guò)引入2個(gè)新的輔助張量，用遞歸積分的方法對(duì)其進(jìn)行實(shí)現(xiàn)，從而避免直接計(jì)算卷積，減少計(jì)算量。張顯文等［23］推導(dǎo)了基于遞歸積分的不分裂復(fù)頻移PML彈性波方程交錯(cuò)網(wǎng)格高階差分法，并對(duì)其進(jìn)行了實(shí)現(xiàn)。Qin Zhen等［24］利用輔助差分方程技術(shù)實(shí)現(xiàn)了復(fù)頻移PML的計(jì)算。無(wú)論是遞歸積分還是輔助差分方程，都可以在不分裂的情形下避免直接計(jì)算卷積，減少計(jì)算量；但是這2種方法推導(dǎo)過(guò)程復(fù)雜繁瑣，不易實(shí)現(xiàn)。筆者推導(dǎo)了一種遞推卷積方法，在交錯(cuò)網(wǎng)格情形下利用遞推的形式計(jì)算卷積，推導(dǎo)過(guò)程直觀易懂，易于編程，應(yīng)用方便簡(jiǎn)單，計(jì)算成本沒(méi)有太大變化。采用這種方法來(lái)求解彈性波方程，與常規(guī)PML的結(jié)果進(jìn)行比較，可知，利用遞推卷積計(jì)算的復(fù)頻移PML能夠有效地改善困難情況下的吸收效果，證明了本文方法的正確性和有效性。

1 基本原理

PML技術(shù)本質(zhì)上是將波動(dòng)方程在PML層內(nèi)進(jìn)行復(fù)坐標(biāo)變換［3］。對(duì)于變換后的坐標(biāo)，方程及其解的形式不變，但對(duì)于原坐標(biāo)，解是衰減的。傳統(tǒng)的PML中頻率域坐標(biāo)變換關(guān)系（以x方向?yàn)槔?/p>

式中：sx是復(fù)拉伸函數(shù)；x是原坐標(biāo)是變換后的坐標(biāo)；ω是圓頻率是衰減系數(shù)。復(fù)頻移拉伸函數(shù)則是對(duì)式（3）進(jìn)行了擴(kuò)展，使其形式更一般化［19］：

其中：αx≥0；Κx≥1?？梢钥闯觯瑥?fù)拉伸函數(shù)是復(fù)頻移拉伸函數(shù)在αx＝0、Κx＝1情況下的特例。

將坐標(biāo)變換關(guān)系變換回時(shí)域，用sx（t）表示1/sx的傅里葉反變換，可以得到

可以看出，坐標(biāo)變換在時(shí)域是卷積關(guān)系，重要的是如何避免直接計(jì)算卷積。

對(duì)式（4）取倒數(shù)并進(jìn)行傅里葉變換可得

式中：δ（x）為沖激函數(shù)；H（t）為單位階躍函數(shù)。

這樣，時(shí)域坐標(biāo)變換關(guān)系就分為2部分：其中第一部分容易計(jì)算，只要計(jì)算空間偏導(dǎo)除以系數(shù)Κ就可以了；最主要的是計(jì)算第二部分的卷積。

在離散情況下，將nΔt時(shí)刻的卷積寫為［20］

由于采用交錯(cuò)網(wǎng)格，則

其中，

因?yàn)槭剑?1）是簡(jiǎn)單的指數(shù)形式，所以可以把式（10）寫成遞推的形式：

其 中：為的簡(jiǎn)記形式；為的簡(jiǎn)記形式。

這樣，在物理區(qū)域進(jìn)行正常計(jì)算，在PML內(nèi)通過(guò)下面的關(guān)系進(jìn)行坐標(biāo)變換后的計(jì)算就可以了：

其中，Ψx可以通過(guò)式（12）的遞推得到。

值得注意的是，式（12）與（13）中的場(chǎng)變量在時(shí)間上不滿足交錯(cuò)布局性。筆者采用如下策略：

1）首先進(jìn)行坐標(biāo)變換：

將式（12）代入式（14）可得

2）利用式（12）更新。

這就是復(fù)頻移PML的遞推卷積實(shí)現(xiàn)的基本原理。其實(shí)就是直接將卷積替換，在離散交錯(cuò)網(wǎng)格情形下通過(guò)遞推計(jì)算，直觀易懂，方便簡(jiǎn)單。它可以很容易地推廣到y(tǒng)，z方向。4個(gè)角上的情形與傳統(tǒng)的PML類似，所有方向同時(shí)考慮就可以了。

2 正演模擬計(jì)算與對(duì)比

2.1 模擬無(wú)限大區(qū)域時(shí)復(fù)頻移PML和常規(guī)PML的吸收性能

為了驗(yàn)證本方法的正確性與有效性，對(duì)各向同性彈性介質(zhì)進(jìn)行了試算，并與傳統(tǒng)的PML的計(jì)算結(jié)果進(jìn)行了對(duì)比。所有算例均采用交錯(cuò)網(wǎng)格高階有限差分法進(jìn)行計(jì)算，時(shí)間二階，空間六階。圖1是層狀介質(zhì)模型，采樣點(diǎn)為301×301，網(wǎng)格間距5m×5 m，上層縱橫波速度分別為3 000m/s、1 600m/s，下層縱橫波速度分別為3 300m/s、1 900m/s，密度均為2 800kg/m3，震源網(wǎng)格點(diǎn)坐標(biāo)為（150，5），為與z軸成逆時(shí)針135°夾角的集中力震源，道間距5 m，采樣率0.5ms，采樣時(shí)間1.5s。圖2是震源主頻為20Hz的z分量炮記錄對(duì)比。如圖中箭頭所示，可以清楚地看到，常規(guī)PML的炮記錄在遠(yuǎn)偏移距處有比較強(qiáng)的虛假反射，尤其是P－P反射波，其兩翼已經(jīng)完全畸變。這是由于在遠(yuǎn)偏移距處波入射到匹配層的入射角過(guò)大造成的。而復(fù)頻移PML的炮記錄則要干凈很多，說(shuō)明復(fù)頻移PML對(duì)掠射情況下波的吸收效果比常規(guī)PML要好得多。為了更直觀地觀察，又對(duì)0.3s時(shí)的z分量波場(chǎng)快照進(jìn)行了對(duì)比，如圖3a、b所示。通過(guò)對(duì)比可以清楚地看到，黑框范圍內(nèi)，上界面附近遠(yuǎn)偏移距處，常規(guī)PML的波場(chǎng)快照有明顯的反射能量，而復(fù)頻移PML的波場(chǎng)快照則幾乎看不到這一點(diǎn)。為了更清楚地顯示，對(duì)上述波場(chǎng)快照的黑框部分進(jìn)行放大，如圖3c、d所示。通過(guò)二者的對(duì)比可以更清楚地看到上述現(xiàn)象（如圖3中黑色箭頭所示）。這進(jìn)一步說(shuō)明了復(fù)頻移PML較常規(guī)PML能夠有效地改善掠射情況下的吸收效果。

圖1 層狀介質(zhì)模型Fig.1 Model of layer medium

為了更進(jìn)一步說(shuō)明這個(gè)問(wèn)題，選擇同樣的速度模型，改變震源主頻后進(jìn)行了同樣的計(jì)算對(duì)比。圖4是震源主頻為30Hz的z分量炮記錄對(duì)比結(jié)果?？梢钥闯觯撼Ｒ?guī)PML的炮記錄在如圖中黑色箭頭所示的大偏移距處有比較強(qiáng)烈的假反射能量，而且隨偏移距增加假反射能量越強(qiáng)烈；這是由于偏移距越大，波入射到匹配層的入射角越大，掠射情況越嚴(yán)重。而復(fù)頻移PML比常規(guī)PML的吸收效果要更好一些。為了更好地說(shuō)明，同樣對(duì)0.3s時(shí)的z分量波場(chǎng)快照進(jìn)行了對(duì)比，如圖5所示。通過(guò)黑色方框和箭頭所示區(qū)域的對(duì)比，也可以驗(yàn)證在掠射情況下復(fù)頻移PML比常規(guī)PML吸收效果好。為驗(yàn)證上述結(jié)論，又抽取了3個(gè)單道在常規(guī)PML和復(fù)頻移PML情況下0～1.2s和0.3～1.2s的波形進(jìn)行了對(duì)比，如圖6所示，圖中3個(gè)單道與炮點(diǎn)縱坐標(biāo)相同，而且隨著橫坐標(biāo)數(shù)值的增大，與炮點(diǎn)距離增加，偏移距也就越大，相應(yīng)地對(duì)匹配層的入射角也越大。從圖6可以看到：在單道1常規(guī)PML和復(fù)頻移PML的效果基本一樣；但是在單道2，常規(guī)PML和復(fù)頻移PML的波形開(kāi)始出現(xiàn)明顯的區(qū)別，常規(guī)PML的波形產(chǎn)生了一些不應(yīng)該有的虛假震蕩；在單道3二者的差別更加明顯，常規(guī)PML的虛假震蕩也更大。這也說(shuō)明了在遠(yuǎn)偏移距處即大入射角的情況下復(fù)頻移PML比常規(guī)PML吸收效果好；而且隨著偏移距增大，入射角增加，效果對(duì)比更加明顯，常規(guī)PML中的虛假反射影響更加嚴(yán)重。

以上結(jié)論可與文獻(xiàn)［22－24］中的結(jié)果進(jìn)行類比，結(jié)論相似；尤其是與文獻(xiàn)［24］，結(jié)果基本一致。這是因?yàn)檫@幾種方法在坐標(biāo)變換中都是采用相同的復(fù)頻移拉伸函數(shù)，所不同的只是實(shí)現(xiàn)方法。相較而言，本文方法更直觀易懂，易于編程；而且本文模擬震源接近研究區(qū)域邊緣的情況而非長(zhǎng)條形介質(zhì)，更符合大多數(shù)的實(shí)際情況。

圖2 主頻為20Hz的炮記錄對(duì)比Fig.2 Comparison of seismogram when dominant frequency is 20Hz

圖3 主頻為20Hz的0.3s波場(chǎng)快照及局部放大圖Fig.3 Snapshot and its partial enlarged detail of 0.3swhen dominant frequency is 20Hz

圖4 主頻為30Hz的炮記錄對(duì)比Fig.4 Comparison of seismogram when dominant frequency is 30Hz

圖5 主頻為30Hz的0.3s波場(chǎng)快照及局部放大圖Fig.5 Snapshot and its partial enlarged detail of 0.3swhen dominant frequency is 30Hz

為考察掠射情況下常規(guī)PML的虛假反射對(duì)復(fù)雜介質(zhì)波場(chǎng)尤其是繞射等弱能量波場(chǎng)的影響以及復(fù)頻移PML在這種情況下的吸收性能，對(duì)含散射體的兩層介質(zhì)模型進(jìn)行了計(jì)算對(duì)比。圖7是含散射體介質(zhì)模型，在網(wǎng)格點(diǎn)坐標(biāo)（225，15）即近地表處有一散射體，散射體縱橫波速度分別為2 700、1 300m/s，其他參數(shù)與前相同。圖8是震源主頻為30Hz的z分量炮記錄對(duì)比。如圖中黑色箭頭所示，可以清楚地看到常規(guī)PML的炮記錄在兩翼有很強(qiáng)的虛假反射能量，嚴(yán)重影響炮記錄的質(zhì)量；尤其是繞射波，由于能量較弱，相對(duì)來(lái)說(shuō)對(duì)其的影響更嚴(yán)重，如直達(dá)縱波的繞射橫波右翼，畸變就比較嚴(yán)重，而直達(dá)橫波的繞射縱波右翼，則幾乎被掩蓋。復(fù)頻移PML的炮記錄相對(duì)來(lái)說(shuō)就好很多，各個(gè)波形都很清楚，這也說(shuō)明了復(fù)頻移PML較常規(guī)PML的效果優(yōu)勢(shì)。圖9是0.3s時(shí)的z分量波場(chǎng)快照對(duì)比，也可以看出常規(guī)PML在掠射情況下產(chǎn)生的比較強(qiáng)的虛假反射對(duì)繞射波等弱能量波場(chǎng)有相當(dāng)?shù)挠绊?，這會(huì)對(duì)實(shí)際應(yīng)用中的成像、反演等其他處理產(chǎn)生不利影響，而復(fù)頻移PML則會(huì)改善這一點(diǎn)。

2.2 復(fù)頻移PML與常規(guī)PML的計(jì)算量與存儲(chǔ)量

二維情況下，在存儲(chǔ)量方面，對(duì)于分裂的常規(guī)PML來(lái)說(shuō)，如果不存儲(chǔ)總波場(chǎng)，那么需要存儲(chǔ)數(shù)組的最大個(gè)數(shù)是10。但是，每次循環(huán)都需要由分裂的波場(chǎng)分量相加計(jì)算總波場(chǎng)。如果存儲(chǔ)總波場(chǎng)，則所需要存儲(chǔ)的數(shù)組最大個(gè)數(shù)為15。對(duì)于不分裂的復(fù)頻移PML來(lái)說(shuō)，所需存儲(chǔ)的數(shù)組最大個(gè)數(shù)為13，與常規(guī)PML相比沒(méi)有太大變化。在計(jì)算量方面，不存儲(chǔ)總波場(chǎng)的常規(guī)PML計(jì)算時(shí)間為295s，復(fù)頻移PML計(jì)算時(shí)間為213s。硬件條件相同，CPU為英特爾E8400，主頻3GHz。可見(jiàn)，復(fù)頻移PML的計(jì)算量由于不分裂而有所減少。

圖6 不同單道的復(fù)頻移PML和常規(guī)PML波形對(duì)比Fig.6 Wave contrast figure of different receivers between CFS－PML and classical PML

2.3 窄區(qū)域自由表面條件下PML和復(fù)頻移PML對(duì)面波的吸收性能

由于越來(lái)越高的工程研究和實(shí)際應(yīng)用要求，對(duì)自由表面條件下近地表地震波傳播的模擬是重要和基本的。瑞雷波相速度與頻率的關(guān)系已經(jīng)被廣泛地用于淺層橫波速度估計(jì)中［25］。因此，產(chǎn)生包含準(zhǔn)確瑞雷波信息的合成地震記錄是任何近地表地震模擬的首要目標(biāo)。然而，在含傳統(tǒng)PML的數(shù)值模擬中，如果模擬區(qū)域頂部自由表面距離底部的PML比較近，由于低頻瑞雷波的傳播方向與底部PML的走向基本平行，瑞雷波在傳播過(guò)程中會(huì)與底部PML相互作用，降低吸收效果，甚至產(chǎn)生不穩(wěn)定性［26］。復(fù)頻移PML同樣會(huì)改善這一點(diǎn)。下面以帶自由表面的各向同性彈性均勻介質(zhì)為例進(jìn)行說(shuō)明。模型大小為15km×1km，密度為2 700kg/m3，縱橫波速度分別為3 200m/s、1 870m/s，震源采用主頻為2.5Hz的雷克子波，置于靠近左邊界距自由地表0.15km的地方。圖10是不同時(shí)刻的z分量波場(chǎng)快照對(duì)比。2s時(shí)面波正要和體波分開(kāi)，可以看出兩者結(jié)果幾乎沒(méi)有差別，但是圖10a左邊有一些輕微的干擾噪音。這是由于震源太靠近左邊界，大量能量突然進(jìn)入匹配層，傳統(tǒng)PML吸收不完全所致，這也說(shuō)明了復(fù)頻移PML比傳統(tǒng)PML吸收效率要高。6s時(shí)面波已在介質(zhì)中傳播相當(dāng)一段距離并與體波完全分離，可以看出圖10c中面波已經(jīng)有所畸變并在其后面產(chǎn)生了一些虛假波形；這是由于面波與底部PML相互作用，降低吸收效果，并影響到物理區(qū)域的緣故。而圖10d中復(fù)頻移PML的吸收效果就很好，瑞雷面波的特征得到了保持。

圖7 含散射體介質(zhì)模型Fig.7 Layer media model with scatterer

3 結(jié)論與討論

1）常規(guī)PML在離散后反射系數(shù)不為0，尤其在大入射角、大偏移距的情況下會(huì)產(chǎn)生很強(qiáng)的假反射，降低吸收效果，對(duì)真實(shí)波場(chǎng)產(chǎn)生比較嚴(yán)重的影響，這種影響隨入射角、偏移距的增加而增大。

2）基于遞推卷積的復(fù)頻移PML技術(shù)能夠增強(qiáng)掠射情況的波場(chǎng)衰減，降低反射系數(shù)，有效改善匹配層的吸收效果。

3）對(duì)于含散射體的復(fù)雜介質(zhì)模型，常規(guī)PML也會(huì)產(chǎn)生假反射，并更加嚴(yán)重地影響散射場(chǎng)等弱能量波場(chǎng)，而且會(huì)對(duì)實(shí)際應(yīng)用中的偏移、反演等處理產(chǎn)生不利影響；而基于遞推卷積的復(fù)頻移PML技術(shù)同樣可以改善吸收效果，降低不利影響。

4）對(duì)于窄區(qū)域自由表面條件下面波的模擬，傳統(tǒng)PML也會(huì)對(duì)波場(chǎng)產(chǎn)生比較嚴(yán)重的影響，而基于遞推卷積的復(fù)頻移PML同樣會(huì)改善這一點(diǎn)。

5）基于遞推卷積的復(fù)頻移PML技術(shù)采用不分裂變量的遞推卷積方法進(jìn)行計(jì)算，實(shí)現(xiàn)時(shí)通過(guò)遞推來(lái)計(jì)算卷積，使其存儲(chǔ)量和計(jì)算量與常規(guī)PML相比都沒(méi)有太大變化，不會(huì)增加計(jì)算成本。

圖8 含散射體模型的z分量炮記錄對(duì)比Fig.8 Comparison of z component seismogram of the model with scatterer

圖9 含散射體模型的z分量0.3s波場(chǎng)快照及局部放大圖對(duì)比Fig.9 Comparison of z component snapshot and its partial enlarged detail of 0.3s

本文基于遞推卷積的復(fù)頻移PML技術(shù)能夠有效改善傳統(tǒng)PML在困難情況下的吸收效果，而且不增加計(jì)算成本，這對(duì)很多如薄片區(qū)域、震源接近研究區(qū)域邊緣、大偏移距接收的體波模擬和窄區(qū)域的面波模擬等情況有很好的應(yīng)用前景。對(duì)于偏移成像等處理中的吸收邊界條件，這也是一種較好的選擇。但是復(fù)頻移PML沒(méi)有改變常規(guī)PML的基本思想，所以它還是存在一些問(wèn)題，比如對(duì)一些各向異性介質(zhì)有固有的不穩(wěn)定性、離散后匹配層反射系數(shù)不嚴(yán)格為0等不足還是不能很好地克服。筆者在匹配層外圍采用的是Dirichlet邊界條件，后續(xù)研究中可以將其替換為旁軸吸收邊界條件以進(jìn)一步提高離散后的吸收效果；另外對(duì)于二階位移波動(dòng)方程的復(fù)頻移PML吸收邊界條件也有待進(jìn)一步研究?；谶f推卷積的復(fù)頻移PML技術(shù)以及后續(xù)相關(guān)方法的進(jìn)一步研究，將有利于西部碳酸鹽巖探區(qū)復(fù)雜近地表速度模型下的正演模擬，尤其是對(duì)碳酸鹽巖探區(qū)近地表散射波機(jī)理認(rèn)識(shí)的研究，為將來(lái)碳酸鹽巖探區(qū)高精度勘探開(kāi)發(fā)服務(wù)。

圖10 帶自由地表模型的z分量波場(chǎng)快照對(duì)比Fig.10 Comparison of z component snapshots with free surface

（References）：

［1］Virieux J.P－SV Wave Propagation in Heterogeneous Media：Velocity－Stress Finite－Difference Method［J］.Geophysics，1986，51（4）：889－901.

［2］Bérenger J P.A Perfectly Matched Layer for the Absorption of Electromagnetic Waves［J］.Journal of Computational Physics，1994，114（2）：185－200.

［3］Chew W C，Weedon W H.A 3－D Perfectly Matched Medium from Modified Maxwell’s Equations with Stretched Coordinates［J］.Microwave and Optical Technology Letters，1994，7（13）：599－604.

［4］Rappaport C M.Perfectly Matched Absorbing Boundary Conditions Based on Anisotropic Lossy Mapping of Space［J］.IEEE Microwave Guided Wave Lett，1995，5（3）：90－92.

［5］Liu Q H，Tao J P.The Perfectly Matched Layer for Acoustic Waves in Absorptive Media［J］.Journal of the Acoustical Society of America，1997，102（4）：2072－2082.

［6］Wang T，Oristaglio M L.3－D Simulation of GPR Surveys Over Pipes in Dispersive Soils［J］.Geophysics，2000，65（5）：1560－1568.

［7］Zeng Y Q，He J Q，Liu Q H.The Application of the Perfectly Matched Layer in Numerical Modeling of Wave Propagation in Poroelastic Media［J］.Geophysics，2001，66（4）：1258－1266.

［8］Chen Y H，Chew W C，Oristaglio M L.Application of Perfectly Matched Layers to the Transient Modeling of Subsurface EM Problems［J］.Geophysics，1997，62（6）：1730－1736.

［9］Engquist B，Majda A.Absorbing Boundary Conditions for the Numerical Simulation of Waves［J］.Mathematics of Computation，1977，31：629－651.

［10］Higdon R L.Absorbing Boundary Condition for E－lastic Waves［J］.Geophysics，1991，56（2）：231－241.

［11］Liao Z P，Wong H L，Yang B P，et al.A Transmitting Boundary for Transient Wave Analysis［J］.Scientia Sinica，1984，27（10）：1063－1076.

［12］Marfurt K J.Accuracy of Finite－Difference and Finite－Element Modeling of the Scalar and Elastic Wave Equations［J］.Geophysics，1984，49（5）：533－549.

［13］Teixeira F L，Chew W C.On Causality and Dynamic Stability of Perfectly Matched Layers for FDTD Simulations［J］.IEEE Transactions on Microwave Theory and Techniques，1999，47（6）：775－785.

［14］Chew W C，Liu Q H.Perfectly Matched Layers for Elastodynamics：A New Absorbing Boundary Condition［J］.Journal of Computational Acoustics，1996，4（3）：341－359.

［15］劉四新，曾昭發(fā)，徐波.地質(zhì)雷達(dá)數(shù)值模擬中有損耗介質(zhì)吸收邊界條件的實(shí)現(xiàn)［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2005，35（3）：378－381.Liu Sixin，Zeng Zhaofa，Xu Bo.Realization of Absorbing Boundary Condition with Lossy Media for Ground Penetrating Radar Simulation［J］.Journal of Jilin University：Earth Science Edition，2005，35（3）：378－381.

［16］殷文.基于頻率域高階有限差分法的正演模擬及并行算法［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2008，38（1）：144－151.Yin Wen.Forward Modeling and Parallel Algorithm Based on High－Order Finite－Difference Method in Frequency Domain［J］.Journal of Jilin University：Earth Science Edition，2008，38（1）：144－151.

［17］韓利，韓立國(guó)，李翔，等.二階聲波方程頻域PML邊界條件及頻域變網(wǎng)格步長(zhǎng)并行計(jì)算［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2011，41（4）：1226－1232.Han Li，Han Liguo，Li Xiang，et al.PML Boundary Conditions for Second－Order Acoustic Wave Equations and Variable Grid Parallel Computation in Frequency－Domain Modeling［J］.Journal of Jilin University：Earth Science Edition，2011，41（4）：1226－1232.

［18］Festa G，Vilotte J P.The Newmark Scheme as Velocity－Stress Time－Staggering：An Efficient PML Implementation for Spectral－Element Simulations of Elastrodynamics［J］.Geophysical Journal International，2005，161：789－812.

［19］Roden J A，Gedney S D.Convolution PML（CPML）：An Efficient FDTD Implementation of the CFS－PML for Arbitrary Media［J］.Microwave and Optical Technology Letters，2000，27（5）：334－339.

［20］Komatitsch D，Martin R.An Unsplit Convolutional Perfectly Matched Layer Improved at Grazing Inci－dence for the Seismic Wave Equation［J］.Geophysics，2007，72（5）：SM155－SM167.

［21］Kuzuoglu M，Mittra R.Frequency Dependence of the Constitutive Parameters of Causal Perfectly Matched Anisotropic Absorbers［J］.IEEE Microwave and Guided Wave Letters，1996，6（12）：447－449.

［22］Drossaert H F，Giannopoulos A.A Nonsplit Complex Frequency－Shifted PML Based on Recursive Integration for FDTD Modeling of Elastic Waves［J］.Geophysics，2007，72（2）：T9－T17.

［23］張顯文，韓立國(guó)，黃玲，等.基于遞歸積分的復(fù)頻移PML彈性波方程交錯(cuò)網(wǎng)格高階差分法［J］.地球物理學(xué)報(bào)，2009，52（7）：1800－1807.Zhang Xianwen，Han Liguo，Huang Ling，et al.A Staggered－Grid High－Order Difference Method of Complex Frequency－Shifted PML Based on Recursive Integration for Elastic Wave Equation［J］.Chinese Journal of Geophysics，2009，52（7）：1800－1807.

［24］Qin Z，Lu M H，Zheng X D，et al.The Implementation of an Improved NPML Absorbing Boundary Condition in Elastic Wave Modeling［J］.Applied Geophysics，2009，6（2）：113－121.

［25］Socco L V，F(xiàn)oti S，Boiero D.Surface－Wave Analysis for Building Near－Surface Velocity Models－Established Approaches and New Perspectives［J］.Geophysics，2010，75（5）：A83－A102.

［26］Festa G，Delavaud E，Vilotte J P.Interaction Between Surface Waves and Absorbing Boundaries for Wave Propagation in Geological Basins：2DNumerical Simulations［J］.Geophysical Research Letters，2005，32：L20306.