張明財，熊章強，張大洲

（1.中國水電工程顧問集團公司西北勘測設(shè)計研究院，陜西西安710043；2.中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南長沙410083）

瑞雷面波憑借其衰減小、信噪比高、抗干擾能力強、分辨率高以及在層狀介質(zhì)中所具有的頻散特性，在地震勘探中有著廣泛的應(yīng)用，其正演模擬也成為國內(nèi)外學(xué)者的研究熱點。Hestholm［1］通過模擬對比二維、三維空間中的瑞雷面波發(fā)現(xiàn)，二維介質(zhì)難以模擬實際的面波傳播特性。董良國等［2－3］對一階彈性波方程交錯網(wǎng)格高階差分解法進行了研究，給出了差分方程的穩(wěn)定性條件。Igel等［4］給出了三維球坐標系下的交錯網(wǎng)格有限差分公式，并對地震波場進行了數(shù)值模擬。Peter等［5］利用四階有限差分交錯網(wǎng)格對三維介質(zhì)中的瑞雷面波進行了數(shù)值模擬，研究了鏡像法處理自由地表時瑞雷面波的精度與最小波長內(nèi)網(wǎng)格點數(shù)的關(guān)系。熊章強［6］通過正演模擬，對復(fù)雜介質(zhì)下的瑞雷面波傳播特征進行研究，利用高階交錯網(wǎng)格有限差分對垂直界面、斜界面、垂直低速帶等復(fù)雜模型進行了正演模擬計算，并提取相應(yīng)的頻散曲線進行了分析。周竹生等［7］利用交錯網(wǎng)格有限差分對彈性介質(zhì)中的瑞雷面波進行正演模擬研究，在人工邊界上采用了變系數(shù)吸收邊界條件，避免了角點繞射，在層狀介質(zhì)（軟弱夾層、連續(xù)介質(zhì)）模型中取得了較好的應(yīng)用效果；裴正林［8］給出了三維各向同性介質(zhì)中一階應(yīng)力－速度彈性波方程交錯網(wǎng)格任意偶數(shù)階精度有限差分格式，推導(dǎo)出三維各向同性彈性介質(zhì)完全匹配層吸收邊界條件公式和相應(yīng)的交錯網(wǎng)格高階有限差分格式。

自由地表的處理是瑞雷面波數(shù)值模擬的重要部分，許多學(xué)者經(jīng)過研究提出的方法主要有如下幾種：以真空形式處理應(yīng)力的自由邊界法［9］、應(yīng)力鏡像法［10］、橫向各向同性介質(zhì)代替自由邊界法［11］和聲學(xué)－彈性介質(zhì)邊界代替自由邊界法［12］。由于聲學(xué)－彈性介質(zhì)邊界代替自由邊界法符合物理規(guī)律，易于編程實現(xiàn)，應(yīng)用效果良好，我們利用該方法處理自由地表。

基于單PC機的串行算法勢必會受到當前計算機數(shù)據(jù)存儲量小、運算速度慢的限制，采用并行計算技術(shù)增加數(shù)據(jù)存儲量、減少正演模擬時間是一條非常有效的途徑。王月英［13］采用基于 MPI的并行算法對三維波動方程進行了有限單元正演模擬，并驗證了該算法的可行性和有效性；Minkoff［14］研究了三維速度－應(yīng)力有限差分地震波數(shù)值模擬方法，通過數(shù)值模擬提高了運算規(guī)模和效率；何兵壽等［15］研究了利用有限差分法并行求解該方程的基本思路與方法，極大地提高了數(shù)值求解彈性波方程的計算效率。

MPI是一種比較著名的應(yīng)用于并行環(huán)境的消息傳遞標準，其目的是為基于消息傳遞的并行程序設(shè)計提供一個高效、可擴展、統(tǒng)一的編程環(huán)境。MPI由消息傳遞庫函數(shù)具體實現(xiàn)，目前應(yīng)用比較廣泛 的 有 MPICH 和 LAMMPI［16］。我 們 選 用MPICH 1.2.7版本。

1 三維彈性波動方程的有限差分數(shù)值模擬算法

在三維各向同性彈性介質(zhì)中，一階速度－應(yīng)力方程可以表述為［1］

式中：i，j，k∈｛x，y，z｝；δij為 Kronecker記號；ρ表示介質(zhì)密度；fi，vi分別表示體力和速度的i分量；σij表示應(yīng)力分量；λ，μ為拉梅常數(shù)。

在交錯網(wǎng)格有限差分計算中，為了保證在空間上具有偶數(shù)階精度，將正應(yīng)力 （σxx，σyy，σzz）和模型參數(shù) （ρ， λ，μ）置于節(jié)點 （i，j，k），切 應(yīng) 力σxy置于節(jié)點（i＋1／2，j＋1／2，k），切應(yīng)力σxz置于節(jié)點（i＋1／2，j，k＋1／2），切應(yīng)力σyz置于節(jié)點（i，j＋1／2，k＋1／2），速度分量u，v，w 分別置于 （i＋1／2，j，k），（i，j＋1／2，k），（i，j，k＋1／2），如圖1所示。

圖1 交錯網(wǎng)格中速度分量和應(yīng)力分量的空間位置關(guān)系

為了使交錯網(wǎng)格有限差分在時間上具有2 M階精度，在t＋Δt／2時刻計算速度分量u，v，w，而在t時刻計算應(yīng)力分量σxx，σxy，σxz，σyy，σyz和σzz。一階速度－應(yīng)力彈性波動方程（1）可用交錯網(wǎng)格有限差分求解。f（x，t）對于時間的導(dǎo)數(shù)?f（x，t）／?t可用中心差分表示，對于空間的導(dǎo)數(shù)?f（x，t）／?x可用向前、向后差分表示，即［2－3］

其中，Δt表示時間步長，Δh表示空間離散網(wǎng)格大小。（2）式即為時間域求導(dǎo)的有限差分離散格式，（3）式和（4）式為空間域求導(dǎo)的有限差分離散格式。將（2）式至（4）式代入一階速度 －應(yīng)力方程（1）中，即可求解出彈性波動方程。

2 自由邊界條件

在有限差分計算中，自由邊界條件直接影響著瑞雷面波數(shù)值模擬的正確性。表1羅列了幾種處理自由地表邊界條件的方法在三維情況下的具體實施過程。我們選用聲學(xué)－彈性介質(zhì)邊界代替自由邊界法，在三維介質(zhì)中將自由邊界網(wǎng)格點的介質(zhì)參數(shù)按表1進行特殊處理，自由邊界以下網(wǎng)格點的介質(zhì)參數(shù)不變，速度分量和應(yīng)力分量按照正常的有限差分離散格式進行計算。

表1 自由地表處理方法

3 并行設(shè)計

利用交錯網(wǎng)格有限差分進行瑞雷面波數(shù)值模擬的特點如下：①在時間域中，速度分量（n＋1）／2時刻的波場值由其（n－1）／2時刻的波場值和應(yīng)力分量n時刻的波場值計算得到；應(yīng)力分量n＋1時刻的波場值由其n－1時刻的波場值和速度分量（n－1）／2時刻的波場值計算得到；②在空間域中，速度和應(yīng)力分量在任意網(wǎng)格點處的波場值只由各自周圍有限個網(wǎng)格點上的波場值計算得到。據(jù)此可知，運算任務(wù)具有良好的局部性和并行度，可將其分解為若干個子任務(wù)，這些子任務(wù)可以并行地執(zhí)行，以期達到并行求解波動方程的目的。

3.1 模型劃分

設(shè)總的模型規(guī)模為N＝Nx×Ny×Nz，其中Nx，Ny，Nz分別為x，y，z方向的網(wǎng)格點數(shù)。將模型沿x，y，z方向分別分解為Px，Py，Pz個小區(qū)間，即可將總的模型規(guī)模分解到NP＝Px×Py×Pz個子進程中，每個子進程中的網(wǎng)格規(guī)模為N／NP。各個子進程均具有自己獨立的內(nèi)存空間，數(shù)據(jù)在整個計算空間中不重復(fù)占用內(nèi)存，程序并沒有因為并行運算而額外增加內(nèi)存開銷。

3.2 數(shù)據(jù)通信

圖2 相鄰進程數(shù)據(jù)通信

在空間四階有限差分數(shù)值模擬時，每個子進程中的模型空間均需沿x，y，z方向劃分出2個網(wǎng)格點組成過渡層。子進程在每個時間步長計算出該進程模型空間的波場值（包括應(yīng)力分量和速度分量）后，必須與其相鄰的進程間通過過渡層進行數(shù)據(jù)通信（發(fā)送或接收）。圖2為二維情形數(shù)據(jù)通信示意圖，其中，淺灰色區(qū)域表示數(shù)據(jù)發(fā)送源地址，深灰色區(qū)域表示數(shù)據(jù)接收目的地址，數(shù)據(jù)按箭頭方向傳遞。圖2中進程1的右部與進程2的左部通信，進程1的上部與進程4的下部通信，進程2的上部與進程3的下部通信，進程3的左部與進程4的右部通信。由于進程1的左部、下部，進程2的右部、下部，進程3的右部、上部及進程4的左部、上部沒有相鄰的進程，故不需要設(shè)置過渡層。

3.3 基于MPI的有限差分計算程序流程

基于MPI的有限差分計算程序設(shè)計流程如圖3所示，具體步驟描述如下。

1）初始化MPI獲取通信器及進程總數(shù)和進程序號。

2）每個子進程讀取指定模型空間的模型參數(shù)，包括縱波速度vP，橫波速度vS，密度ρ，主進程讀取檢波器坐標、震源信息。

3）主進程將震源信息、檢波器坐標廣播到子進程中，將各子進程中的模型參數(shù)進行歸約，判斷頻散條件和穩(wěn)定性條件。

4）時間循環(huán)開始，各個子進程判斷震源坐標是否位于本進程計算區(qū)域，如果是，則讀取震源函數(shù)，將震源施加于相應(yīng)的位置。

5）各個子進程計算波場，并根據(jù)模型關(guān)系處理自由邊界和吸收邊界。以圖2所示的進程為例，進程3、進程4的上邊界為自由邊界，進程1、進程4的左邊界，進程2、進程3的右邊界，進程1、進程2的下邊界為吸收邊界。

6）進程之間通過阻塞式標準消息發(fā)送和接受函數(shù)進行數(shù)據(jù)的交換，更新每個進程中的波場信息。

7）對各個進程中的波場信息進行歸約，由主進程輸出指定時刻的波場快照。

8）判斷時間循環(huán)是否結(jié)束，如果沒有結(jié)束，則轉(zhuǎn)至步驟4）再次循環(huán)；如果結(jié)束，則調(diào)用 MPI函數(shù)MPI＿Finalize退出MPI，結(jié)束程序。

3.4 PC機多進程并行效率比較

設(shè)計一模型，大小為100m×100m×100m，網(wǎng)格大小dx＝dy＝dz＝1m，總的網(wǎng)格數(shù)為N＝1×106個，采樣間隔為0.4ms，總采樣時間長度為0.1s。其三維速度－應(yīng)力有限差分數(shù)值模擬過程中共需存儲9個三維數(shù)組（3個速度分量和6個應(yīng)力分量），在單機運行時，需消耗9×N×8字節(jié)的內(nèi)存，若將2臺PC機組成一個PC集群，則每臺PC機僅需消耗9×N×4字節(jié)的內(nèi)存。在雙核PC機上分別開設(shè)1，2，4，8個進程對上述模型進行有限差分數(shù)值模擬，統(tǒng)計出不同進程數(shù)時的串行部分計算時間WS，并行部分計算時間 WP，采用Amdnhl加速性能定律來求取加速比S：

圖3 基于MPI的有限差分計算程序流程

式中：P為總進程數(shù)。表2給出了不同進程數(shù)的串行、并行CPU時間及加速比，其中進程數(shù)為1表示只進行串行計算。由表2可知：

1）隨著進程數(shù)P的增加，每個進程的網(wǎng)格規(guī)模會減小，為待模擬網(wǎng)格規(guī)模的1／P。

2）進行并行計算（即進程數(shù)大于1時）后，串行計算的時間略有增加，并行計算的時間略有減少，但總的運行時間約為只進行串行計算時的1／2，這是因為計算所用的PC機只具備2個運算核心，當進程數(shù)大于2時，每個核心運行多個進程，單核中的兩個進程之間增加的通信任務(wù)會消耗CPU時間。

3）加速比S隨著進程數(shù)P的增加而增加，但與進程數(shù)增加速度的差距越來越大，這是因為①進程數(shù)越多，勢必會增大各進程之間的通信時間；②并行并非嚴格執(zhí)行，某些進程中會執(zhí)行一些必要的輸入輸出，此時其余的進程可能處于等待狀態(tài)。

表2 不同進程數(shù)并行計算加速效率比較

4 瑞雷面波模擬算例

4.1 均勻半空間中的瑞雷面波

均勻各向同性彈性半空間介質(zhì)是一種理想的模型，也是最早證明瑞雷面波存在的假設(shè)，它為研究瑞雷面波的某些傳播特性提供了簡單方便的方法。為了驗證本文有限差分并行計算程序的正確性，設(shè)計均勻各向同性介質(zhì)模型，大小為400m×400m×400m，網(wǎng)格間隔dx＝dy＝dz＝1m，縱波速度vP＝1 000m／s，橫波速度vS＝577m／s，密度ρ＝2 100kg／m3；震源為25Hz的高斯一階導(dǎo)數(shù)函數(shù)，位于坐標（200m，200m，0）處，采樣間隔為0.4ms，總采樣時間長度為400ms，延時40ms。在高性能計算機上申請4個節(jié)點，每個節(jié)點上申請8個CPU，共計32個CPU，每個CPU開設(shè)1個計算進程，1個進程中的網(wǎng)格規(guī)模為100×104個。

圖4為距震源50m（150m，200m，0）和100m（100m，200m，0）處接收的地震信號與解析解對比的結(jié)果，解析解由Lamb問題解析法得到，其中的Green函數(shù)計算部分由Cagniard－de Hoop法解析得到。圖4a和圖4b分別表示距震源50m（150m，200m，0）處的水平和垂直方向速度分量，圖4c和圖4d分別表示距震源100m（100m，200m，0）處的水平和垂直方向速度分量，兩處速度分量都與解析解吻合很好，說明本文有限差分并行計算程序是正確的，用聲學(xué)－彈性介質(zhì)近似自由地表是行之有效的。

圖4 數(shù)值解與解析解對比

圖5為60，160，260ms時vx，vy，vz分量的波場快照，以震源為頂點對模型區(qū)域做切割（即x＝200m，y＝200m，z＝0）處理。分析波場快照可知：瑞雷面波沿自由地表（z＝0）傳播，其能量強于體波，且衰減很慢，特別是從vz分量的波場快照可以發(fā)現(xiàn)，沿深度方向，瑞雷面波的能量大多集中于25m以內(nèi)，即一個波長之內(nèi)，深度方向超過一個波長后，P波逐漸變得清晰，說明瑞雷面波的穿透深度為一個波長。

為了更加清晰地分析波場快照，選取y＝200m所在平面內(nèi)的vz分量260ms時波場快照，如圖6所示，震源位于（x＝200m，z＝0）處。在波場快照中可以很清楚地看到能量很強的瑞雷面波、橫波和縱波。此外還可以看出，在吸收邊界處沒有產(chǎn)生強烈的反射波，說明程序中使用的吸收邊界條件是有效的，達到了預(yù)期的目的。

圖7為各分量的波場記錄。由圖7不難發(fā)現(xiàn)：各分量波場記錄中面波的能量強于直達P波，瑞雷面波的同相軸都為直線，說明瑞雷面波在均勻各向同性介質(zhì)中沒有發(fā)生頻散現(xiàn)象，這與地震勘探理論一致，也進一步證明了本文有限差分計算程序的正確性。

應(yīng)用相移法提取vx分量的頻散曲線，結(jié)果如圖8所示。由圖8可以直觀地看出均勻介質(zhì)中瑞雷面波的頻散曲線為一條直線，說明在均勻介質(zhì)中瑞雷面波不會發(fā)生頻散。

4.2 水平層狀速度遞增模型介質(zhì)中的瑞雷面波

通常情況下地層速度隨著深度的增大而增大，故設(shè)計一個3層速度遞增模型，如表3所示。數(shù)值模擬過程中采用25Hz的高斯一階導(dǎo)數(shù)作為激發(fā)震源，模型大小為200m×200m×100m，網(wǎng)格間距取dx＝dy＝dz＝0.5m，模型網(wǎng)格為400個×400個×200個，開辟32個進程進行并行計算。

表3 速度遞增模型參數(shù)

圖9為水平層狀速度遞增模型vx，vy，vz分量的波場快照。由圖9可以看出，由于各層介質(zhì)的參數(shù)不同，波場快照中出現(xiàn)了很明顯的頻散，大部分能量聚集在模型表面附近。由于分層界面的存在，震源激發(fā)的SV波在界面處發(fā)生反射，此反射波以大于臨界角的角度入射到自由地表會轉(zhuǎn)換為P波，該P波為非均勻平面波；震源激發(fā)的P波以球面波的形式在介質(zhì)中傳播，該球面波可由Sommerfeld積分公式分解為很多非均勻平面波［17］，其中的非均勻SV波和由反射形成的非均勻P波疊合形成頻散的瑞雷面波。隨著時間的推移，在各界面處發(fā)生多次反射，形成復(fù)雜的非均勻P波和SV波，這些波相互疊合形成多種模態(tài)的瑞雷面波，其幅度隨著深度的增加呈指數(shù)遞減。

圖10a為x＝100m平面上記錄到的vz分量單炮波形記錄。為了進一步分析頻散現(xiàn)象，利用相移法提取vz分量波形記錄頻散曲線，如圖10b所示。圖10b中白色方塊為采用快速標量傳遞矩陣法［18］得到的解析解，可見基階模式和高階模式的頻散曲線都吻合得比較好。這進一步說明了瑞雷面波模擬方法的正確性。

圖10 單炮波形記錄及其頻散曲線

5 結(jié)束語

彈性波有限差分數(shù)值模擬具有良好的并行特征，基于MPI的并行計算通過將待模擬區(qū)域劃分為若干個子進程中的子區(qū)域，每個子進程相互通信完成消息傳遞，有效地增加了三維彈性波數(shù)值模擬的網(wǎng)格規(guī)模，提高了計算效率，為進一步研究三維空間中瑞雷面波的傳播特性奠定了基礎(chǔ)。

致謝：感謝中南大學(xué)高性能計算中心提供計算平臺。

［1］Hestholm S O.Finite－difference seismic wave modeling including surface topography［D］.Texas：Univer－sity of Houston，1999

［2］董良國，馬在田，曹景忠，等.一階彈性波方程交錯網(wǎng)格高階差分解法［J］.地球物理學(xué)報，2000，43（3）：411－419 Dong L G，Ma Z T，Cao J Z，et al.A staggered－grid high－order difference method of one－order elastic wave equation ［J］.Chinese Journal Geophysics，2000，43（3）：411－419

［3］董良國，馬在田，曹景忠.一階彈性波方程交錯網(wǎng)格高階差分解法穩(wěn)定性問題［J］.地球物理學(xué)報，2000，43（6）：856－864 Dong L G，Ma Z T，Cao J Z.A study on stability of the staggered－grid high－order difference method of first－order elastic wave equation［J］.Chinese Journal Geophysics，2000，43（6）：856－864

［4］Igel H，Nissen－Meyer T，Jahnke G.Wave propagation in 3Dspherical sections：effects of subduction zones［J］.Physics of the Earth and Planetary Interiors，2002，132（1）：219－234

［5］Peter M，Jozef K，Ladislav H.3Dfourth－order staggered－grid finite－difference schemes：stability and grid dispersion［J］.Bulletin of the Seismological Society of America，2000，90（3）：587－603

［6］熊章強.復(fù)雜介質(zhì)中瑞雷面波的正演模擬及傳播特征研究［D］.武漢：中國地質(zhì)大學(xué)（武漢），2006 Xiong Z Q.Forward modeling and propagation characteristics research of Rayleigh wave in complex media［D］.Wuhan：China University of Geosciences（Wuhan），2006

［7］周竹生，劉喜亮，熊孝雨.彈性介質(zhì)中瑞雷面波有限差分法正演模擬［J］.地球物理學(xué)報，2007，50（2）：567－573 Zhou Z S，Liu X L，Xiong X Y.Finite－difference modeling of Rayleigh surface wave in elastic media［J］.Chinese Journal Geophysics，2007，50（2）：567－573

［8］裴正林.三維各向同性介質(zhì)彈性波方程交錯網(wǎng)格高階有限差分法模擬［J］.石油物探，2005，44（4）：308－315 Pei Z L.Numerical simulation of elastic wave equation in 3－D isotropic media with staggered－grid highorder difference method［J］.Geophysical Prospecting for Petroleum，2005，44（4）：308－315

［9］Zahradník J，Priolo E.Heterogeneous formulations of elastodynamic equations and finite－difference schemes［J］.Geophysical Journal International，1995，120：663－676

［10］Graves R W.Simulating seismic wave propagation in 3Delastic media using staggered－grid finite differences［J］.Bulletin of the Seismological Society of A－merica，1996，86（4）：1091－1106

［11］Mittet R.Free－surface boundary conditions for elastic staggered－grid modeling schemes［J］.Geophysics，2002，67（5）：1616－1623

［12］Xu Y X，Xia J H，Miller R D.Numerical investigation of implementation of air－earth boundary by a－coustic－elastic boundary approach［J］.Geophysics，2007，72（5）：147－153

［13］王月英.基于MPI的的三維波動方程有限元法并行正演模擬［J］.石油物探，2009，48（3）：221－225 Wang Y Y.3Dwave－equation finite－element forward modeling based on message passing interface parallel computing［J］.Geophysical Prospecting for Petroleum，2009，48（3）：221－225

［14］Minkoff S E.Spatial parallelism of a 3Dfinite difference velocity－stress elastic wave propagation code［J］.SIAM Journal on Scientific Computing，2002，24（1）：1－19

［15］何兵壽，張會星，韓令賀.彈性波方程正演的粗粒度并行算法［J］.地球物理學(xué)進展，2010，25（2）：650－656 He B S，Zhang H X，Han L H.Forward modeling of elastic wave equation with coarse－grained parallel algorithm［J］.Progress in Geophysics，2010，25（2）：650－656

［16］張林波，遲學(xué)斌，莫則堯，等.并行計算導(dǎo)論［M］.北京：清華大學(xué)出版社，2006：175－182 Zhang L B，Chi X B，Mo Z Y，et al.Introduction to parallel computing［M］.Beijing：Tsinghua University Press，2006：175－182

［17］巴特.地震學(xué)的數(shù)學(xué)問題［M］.鄭治真，譯.北京：科學(xué)出版社，1976：229－234 Bart M.Mathematical problem of seismology［M］.Zhen Z Z，translator.Beijing：Science Press，1976：229－234

［18］凡友華，肖柏勛，劉家琦.計算層狀介質(zhì)中軸對稱柱面瑞雷面波頻散函數(shù)的δ矩陣法［J］.物探與化探，2001，25（2）：109－116 Fan Y Y，Xiao B X，Liu J Q.Theδmatrix method for the dispersion function computation of axis symmetrical cylindrical Rayleigh wave in multilayered media［J］.Geophysics ＆ Geochemical Exploration，2001，25（2）：109－116