基于MPI的三維瑞雷面波有限差分并行模擬

張明財(cái)，熊章強(qiáng)，張大洲

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

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

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

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

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

1 三維彈性波動(dòng)方程的有限差分?jǐn)?shù)值模擬算法

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

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

在交錯(cuò)網(wǎng)格有限差分計(jì)算中，為了保證在空間上具有偶數(shù)階精度，將正應(yīng)力 （σxx，σyy，σzz）和模型參數(shù) （ρ， λ，μ）置于節(jié)點(diǎn) （i，j，k），切 應(yīng) 力σxy置于節(jié)點(diǎn)（i＋1／2，j＋1／2，k），切應(yīng)力σxz置于節(jié)點(diǎn)（i＋1／2，j，k＋1／2），切應(yīng)力σyz置于節(jié)點(diǎn)（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 交錯(cuò)網(wǎng)格中速度分量和應(yīng)力分量的空間位置關(guān)系

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

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

2 自由邊界條件

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

表1 自由地表處理方法

3 并行設(shè)計(jì)

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

3.1 模型劃分

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

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

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

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

3.3 基于MPI的有限差分計(jì)算程序流程

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

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

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

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

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

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

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

7）對(duì)各個(gè)進(jìn)程中的波場信息進(jìn)行歸約，由主進(jìn)程輸出指定時(shí)刻的波場快照。

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

3.4 PC機(jī)多進(jìn)程并行效率比較

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

圖3 基于MPI的有限差分計(jì)算程序流程

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

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

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

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

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

4 瑞雷面波模擬算例

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

均勻各向同性彈性半空間介質(zhì)是一種理想的模型，也是最早證明瑞雷面波存在的假設(shè)，它為研究瑞雷面波的某些傳播特性提供了簡單方便的方法。為了驗(yàn)證本文有限差分并行計(jì)算程序的正確性，設(shè)計(jì)均勻各向同性介質(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ù)，位于坐標(biāo)（200m，200m，0）處，采樣間隔為0.4ms，總采樣時(shí)間長度為400ms，延時(shí)40ms。在高性能計(jì)算機(jī)上申請(qǐng)4個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)上申請(qǐng)8個(gè)CPU，共計(jì)32個(gè)CPU，每個(gè)CPU開設(shè)1個(gè)計(jì)算進(jìn)程，1個(gè)進(jìn)程中的網(wǎng)格規(guī)模為100×104個(gè)。

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)束語

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

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

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

［2］董良國，馬在田，曹景忠，等.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法［J］.地球物理學(xué)報(bào)，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］董良國，馬在田，曹景忠.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法穩(wěn)定性問題［J］.地球物理學(xué)報(bào)，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］熊章強(qiáng).復(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é)報(bào)，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ì)彈性波方程交錯(cuò)網(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的的三維波動(dòng)方程有限元法并行正演模擬［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］何兵壽，張會(huì)星，韓令賀.彈性波方程正演的粗粒度并行算法［J］.地球物理學(xué)進(jìn)展，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é)斌，莫?jiǎng)t堯，等.并行計(jì)算導(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］凡友華，肖柏勛，劉家琦.計(jì)算層狀介質(zhì)中軸對(duì)稱柱面瑞雷面波頻散函數(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