實現(xiàn)復頻移完全匹配層吸收邊界條件的遞推卷積方法

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

1.中國石油大學地球科學與技術學院，山東 青島 266580

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

0 引言

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

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

1 基本原理

PML技術本質上是將波動方程在PML層內進行復坐標變換［3］。對于變換后的坐標，方程及其解的形式不變，但對于原坐標，解是衰減的。傳統(tǒng)的PML中頻率域坐標變換關系（以x方向為例）為

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

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

將坐標變換關系變換回時域，用sx（t）表示1/sx的傅里葉反變換，可以得到

可以看出，坐標變換在時域是卷積關系，重要的是如何避免直接計算卷積。

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

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

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

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

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

其中，

因為式（11）是簡單的指數(shù)形式，所以可以把式（10）寫成遞推的形式：

其 中：為的簡記形式；為的簡記形式。

這樣，在物理區(qū)域進行正常計算，在PML內通過下面的關系進行坐標變換后的計算就可以了：

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

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

1）首先進行坐標變換：

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

2）利用式（12）更新。

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

2 正演模擬計算與對比

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

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

圖1 層狀介質模型Fig.1 Model of layer medium

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

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

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

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

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

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

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

2.2 復頻移PML與常規(guī)PML的計算量與存儲量

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

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

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

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

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

3 結論與討論

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

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

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

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

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

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

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

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

圖10 帶自由地表模型的z分量波場快照對比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ā)，徐波.地質雷達數(shù)值模擬中有損耗介質吸收邊界條件的實現(xiàn)［J］.吉林大學學報：地球科學版，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］.吉林大學學報：地球科學版，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］韓利，韓立國，李翔，等.二階聲波方程頻域PML邊界條件及頻域變網(wǎng)格步長并行計算［J］.吉林大學學報：地球科學版，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］張顯文，韓立國，黃玲，等.基于遞歸積分的復頻移PML彈性波方程交錯網(wǎng)格高階差分法［J］.地球物理學報，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.