周 麗，顧漢明，成景旺，劉春成，劉志斌，楊小春

(1.中國地質(zhì)大學(武漢)地球物理與空間信息學院，湖北武漢430074;2.中國地質(zhì)大學(武漢)地球內(nèi)部多尺度成像湖北省重點實驗室，湖北武漢430074;3.長江大學地球物理與石油資源學院，湖北武漢430100;4.中海油研究總院，北京100027)

海上多波多分量地震勘探采用海底OBC采集方式，將檢波器安置在海底接收。該檢波器內(nèi)部包含了壓力檢波器(水聽器)和3個相互垂直的速度檢波器(X分量、Y分量和Z分量)，因此也將海上多波多分量技術稱為海上四分量(M4C)地震勘探技術[1]。這種采集方式不僅能夠接收傳統(tǒng)的縱波，而且能夠接收橫波以及轉(zhuǎn)換波，是真正意義上的全波勘探，可以攜帶更加豐富的海底信息。目前為止M4C技術已解決了海上勘探中的多個難題[2-5]。

海上作業(yè)實驗成本高，耗時費力，所得出的觀測系統(tǒng)參數(shù)不具有完備性。而基于三維波動方程的并行數(shù)值模擬不僅能夠快速模擬復雜介質(zhì)中波的傳播特征，而且能模擬各種不同的觀測系統(tǒng)參數(shù)下的炮集記錄，基于模擬記錄的偏移成像處理效果優(yōu)選觀測系統(tǒng)參數(shù)。三維交錯網(wǎng)格有限差分由于其占內(nèi)存小、模擬精度高、易于并行等優(yōu)點已被很多學者用于復雜介質(zhì)的地震正演模擬中。Virieux[6-7]使用了時間二階精度和空間二階精度的交錯網(wǎng)格有限差分格式，對一階應力-速度方程進行了差分離散，模擬了P-SV波和P-SH波的傳播;Levander[8]在空間和時間都為二階精度的基礎上，將空間上的有限差分格式改為四階，而時間仍保持二階精度，將模擬精度進一步提高;董良國等[9]首先進行了一階彈性波方程的交錯網(wǎng)格高階有限差分模擬，推導了高階有限差分的格式，給出差分格式的穩(wěn)定性條件;裴正林[10-11]將交錯網(wǎng)格有限差分引入到了三維中，進行了三維各向同性介質(zhì)和各向異性介質(zhì)中彈性波的模擬。

三維地震波波動方程數(shù)值模擬對計算機的內(nèi)存要求很大，目前單版PC機不能滿足較大規(guī)模的正演計算，尤其是適應海上寬方位OBC地震采集的正演模擬。大型的并行計算機、對稱多處理機、分布式共享存儲處理機、PC-Cluster集群以及圖形處理器(Graphic Processing Unit,GPU)計算等的出現(xiàn)，都為地震數(shù)值模擬提供了一個很好的計算平臺。其中PC-Cluster集群搭建方便，運算性能好，相對來說較廉價，已經(jīng)成為了高性能的主流產(chǎn)品[12]。常用的并行模式主要有MPI和并行虛擬機(Parallel Virtual Machine,PVM)，其中MPI是目前廣泛使用的并行編程工具，它具有移植性好、功能強大、高效實用等多種特點[13]。MPI是一個庫，共有上百個函數(shù)調(diào)用接口，目前廣泛使用的為OPENMPI和MPICH，本文的并行計算是在MPICH2平臺上實現(xiàn)的。

并行計算為大規(guī)模的三維地震數(shù)值模擬注入了新的活力。Furumura等[14]采用了偽譜法與有限差分混合并行法對1999年的臺灣天然地震進行了模擬;Bohlen[15]采用交錯網(wǎng)格有限差分法進行了粘彈介質(zhì)的波動方程模擬;Micikevicius[16]使用CUDA(Compute Unified Device Architecture)平臺實現(xiàn)了三維有限差分的GPU計算;王德利等[17]采用了Bohlen的方法進行三維粘彈介質(zhì)的地震波場并行模擬;王月英[18]采用有限元法進行了聲波波動方程的并行計算;何兵壽等[19]從二維彈性波動方程出發(fā)，給出了適用求解的計算空間劃分方法及通信方案;魏星等[20]和秦艷芳等[21]結(jié)合偽譜法和有限差分的各自優(yōu)點，在X和Y方向使用偽譜法計算，在Z方向使用交錯網(wǎng)格計算，并在Z方向上進行了并行設計;公續(xù)飛等[22]描述了三維彈性波有限差分法的GPU實現(xiàn)方法;段玉婷等[23]對精細積分法的實現(xiàn)過程進行GPU并行，實現(xiàn)了三維地震正演模擬;張明財?shù)萚24]基于MPI進行了三維瑞雷波的有限差分并行模擬。

我們針對海底OBC采集特點，采用基于MPI的交錯網(wǎng)格有限差分法進行了海上OBC三維地震多波多分量地震觀測的正演模擬，并行設計中將計算任務沿X，Y，Z3個方向分配成多個子區(qū)域，每個子區(qū)域采用一個進程進行計算，通過簡單的3層模型驗證該并行方案的效率，最后進行了海上某靶區(qū)OBC多波多分量地震觀測的正演并行數(shù)值模擬。

1 三維彈性波方程的有限差分正演算法

在三維各向同性介質(zhì)中，利用交錯網(wǎng)格有限差分法進行正演模擬時采用的波動方程形式為一階應力-速度方程，其可以表示為

(1a)

(1b)

式中：vx，vy和vz表示速度分量;σij(i,j=x,y,z)表示應力分量;λ和μ為拉梅常數(shù)。在交錯網(wǎng)格有限差分中，變量的導數(shù)是在相應的變量網(wǎng)格半節(jié)點上計算的。為了保證在空間上具有偶數(shù)階精度，彈性參數(shù)以及波場分量的空間分布如表1和圖1所示。其中正應力σxx，σyy和σzz置于節(jié)點(i,j,k)上;剪應力σxy，σxz和σyz分別置于節(jié)點(i+1/2,j+1/2,k)，(i+1/2,j,k+1/2)和(i,j+1/2,k+1/2)上;速度分量vx，vy和vz分別置于節(jié)點(i+1/2,j,k)，(i,j+1/2,k)和(i,j,k+1/2)上。同時保證在時間上具有偶數(shù)階精度，應力分量在半時刻t+Δt/2取值，速度分量在整時刻t取值。波動方程中一階空間導數(shù)采用公式(2)的離散格式[25]:

(2)

(3)

表1 彈性波波場分量和彈性參數(shù)的空間位置

圖1 交錯網(wǎng)格應力分量與速度分量的空間位置

PML完全匹配層法主要衰減沿著模型邊界法向方向傳播的地震波，不衰減平行于模型邊界傳播的地震波。完全匹配層法最主要的是衰減因子的求取，本文衰減因子d(x)的計算如下：

(4)

式中：vPmax為最大縱波速度;δ和R分別表示PML層的寬度和理論反射系數(shù)。根據(jù)PML的方程分解思路，設v=v‖+v⊥，σ=σ‖+σ⊥，且dx，dy和dz分別為X,Y,Z方向上的阻尼因子。以vx變量為例(其它變量類推)，三維各向同性介質(zhì)中最佳匹配吸收邊界條件可以表示為

(5a)

(5b)

(5c)

(5d)

2 基于MPI的三維并行計算的整體設計

三維彈性波正演模擬中，速度分量和應力分量一共有9個變量，再加上密度ρ和拉梅常數(shù)λ，μ，則在正演過程中至少要申請12個三維數(shù)組，假設每個變量采用雙精度浮點數(shù)占8個字節(jié)且三維模型空間網(wǎng)格數(shù)為Nx,Ny,Nz，則計算過程中一共需要申請的內(nèi)存大約為Nx×Ny×Nz×12×8個字節(jié)?？梢娙S彈性波正演所需內(nèi)存很大，尤其是在模型較大的時候，單個PC機根本無法進行計算，必須通過并行將所需內(nèi)存劃分到多個子空間中分別存儲，才能實現(xiàn)三維波動方程的求解。

交錯網(wǎng)格有限差分格式中，無論是速度分量還是應力分量，任何網(wǎng)格點處波場值的計算只需要周圍有限個點即可，因此計算過程具有良好的局部性，非常有利于并行計算。計算過程中，可將整個模型分成多個子區(qū)域進行計算，所有子區(qū)域并行執(zhí)行，各個子區(qū)域只需要交換子區(qū)域邊界處的波場值便能達到整個模型的并行計算。

2.1 三維并行相鄰進程的交換設計

假設正演模型網(wǎng)格數(shù)為N=NxNyNz，其中Nx，Ny，Nz分別為X，Y和Z方向的網(wǎng)格點數(shù)。并行計算中，假設X，Y，Z方向分別采用Px，Py，Pz個進程計算，則總的計算進程數(shù)為Nproc=PxPyPz，每個子空間計算網(wǎng)格數(shù)為Npe=N/Nproc。計算所需內(nèi)存被分配到Nproc個子空間上，每個子空間具有自己獨立的內(nèi)存空間。由此可見，除了由于通信需要，每個子空間邊界需要增加的網(wǎng)格點數(shù)以外，該并行方案在計算過程中不重復占用空間，程序并沒有因為并行增加額外的巨大的內(nèi)存空間。

若有限差分過程中空間導數(shù)采用2N階高階精度，則根據(jù)交錯網(wǎng)格有限差分格式，每個子空間在X，Y，Z3個方向均需要增加N個網(wǎng)格點作為緩沖過渡層，用來與相鄰子進程間進行通信與交換，保存從相鄰子空間交換得來的數(shù)據(jù)。圖2為三維并行過程中相鄰子空間的數(shù)據(jù)通信示意圖，每個子空間均需要和周圍的上、下、左、右、前、后6個子空間進行交換。圖2中淺灰色部分表示發(fā)送源地址，深灰色部分表示接收地址。并行計算過程中，將所有計算進程按X，Z，Y的方向順序依次進行標識排序，則每個進程對應一個標識myid。同時設置一個三維進程坐標系，其中X，Y，Z方向的取值范圍分別為：0～Px-1，0～Py-1和0～Pz-1，因此每個進程對應一個三維坐標，那么X坐標為0的子模塊左邊界、X坐標為Px-1的子模塊右邊界、Y坐標為0的子模塊前邊界、Y坐標為Py-1的子模塊后邊界和Z坐標為Pz-1的子模塊下邊界要用吸收邊界條件，而Z坐標為0的子模塊上邊界采用吸收邊界條件或自由邊界條件。同時這些邊界處的衰減因子d(x)按照公式(4)計算，而每個進程的其它邊界則令d(x)=0。

圖2 三維相鄰進程交換示意圖解

現(xiàn)假設某一進程標識號為myid，則與該進程交換的周圍6個進程的標識號分別為myid-1(左)、myid+1(右)、myid-Px(上)、myid+Px(下)、myid-PxPz(前)和myid+PxPz(后)。有的進程(當X坐標為0或Px-1，Y坐標為0或Py-1，Z坐標為0或Pz-1)的某些相鄰進程不存在時，則引入虛擬進程MPI_PROC_NULL。虛擬進程是不存在的假想進程，在MPI中的主要作用是充當真實進程通信的目或源。當一個真實進程向一個虛擬進程發(fā)送數(shù)據(jù)或從一個虛擬進程接收數(shù)據(jù)時，該真實進程會立即正確返回，如同執(zhí)行了一個空操作。

2.2 基于MPI的海上OBC地震觀測并行模擬計算流程

結(jié)合交錯網(wǎng)格有限差分格式以及海上OBC觀測系統(tǒng)的特點，設計了并行計算流程。

1) 首先，調(diào)用MPI初始化程序，啟動并行計算環(huán)境，并獲取進程序號以及總的計算進程個數(shù)Nproc。

2) 主進程(0進程)讀取正演相關參數(shù)、震源文件以及接收點文件。其中正演參數(shù)包括時間步長、空間步長、震源主頻、單炮記錄時間長度、各個方向進程個數(shù)(Px，Py，Pz)。判斷正演穩(wěn)定性條件是否滿足，以及總的計算進程PxPyPz是否等于Nproc，如果不滿足這些條件，則重新修改參數(shù)。

3) 各個進程讀取相應的模型縱、橫波速度值和密度值，并根據(jù)讀取的縱、橫波速度和密度求取拉梅參數(shù)等。根據(jù)橫波速度值(橫波速度由0值變成非零值的分界處)求取每個接收點對應的Z方向海底位置，用來保存海底OBC的波場值;同時計算該進程的三維坐標，根據(jù)坐標值判斷每個進程的哪些邊界需要進行邊界條件的處理以及使用哪種邊界條件(自由邊界條件還是PML邊界條件)。

4) 時間循環(huán)開始，判斷震源屬于哪一個進程，則在該進程對應的網(wǎng)格點上加入震源函數(shù)。

5) 各個進程速度場并行計算，同時根據(jù)步驟3)的判斷，進行邊界的處理。將需要交換的邊界處的速度場保存至交換緩沖區(qū)，然后和周圍的6個進程進行交換，以便進行下一個時刻的應力場的計算。交換過程中，首先確定需要交換的周圍6個進程的標識號，然后調(diào)用標準非阻塞通信進行進程之間的數(shù)據(jù)交換與通信。

6) 同步驟4)一樣，進行應力場的計算。

7) 判斷各個進程中各個接收點所在位置，進行波場信息規(guī)約，輸出指定時刻的快照，保存該時刻接收點的波場值。判斷時間循環(huán)是否結(jié)束，若沒有結(jié)束，跳出步驟4)繼續(xù)計算。

2.3 并行效率分析

本次并行測試所使用集群的硬件環(huán)境配置為：①管理節(jié)點一個，計算節(jié)點30個;②節(jié)點CPU為Intel(R)Quad Core E5520 Xeon(R) CPU，2.26GHz×2，8核;③內(nèi)存為16GB(4×4GB)，1066MHz;④磁盤采用15T磁盤陣列。

通過一個簡單的三維層狀模型(圖3)來進行并行設計的效率性分析。該模型大小為3000m×2400m×1800m，網(wǎng)格數(shù)為Nx=600，Ny=480，Nz=360，空間步長dx=dy=dz=5m，時間步長dt=0.5ms，記錄長度選取2s，正演過程中采用空間4階、時間2階計算精度，震源采用30Hz的雷克子波，且激發(fā)點位于(1500m，1200m，150m)處，模型周邊均使用吸收邊界條件，模型參數(shù)見表2。為了方便觀察地震波在各層之內(nèi)的傳播特征，分別在X方向1500m，Y方向1700m，Z方向900m截取快照剖面，圖3b到圖3d 分別為1000ms時刻3個分量的波場快照。快照圖中，各個界面波場特征明顯，可以看到第一層海水中只存在縱波沒有橫波，而且在vx分量的波場快照中，X方向1500m(經(jīng)過震源)剖面上沒有波場傳播，這也驗證了各向同性介質(zhì)中縱波源激發(fā)不會產(chǎn)生SH波的傳播。圖3中①代表第一界面產(chǎn)生的反射縱波;②代表第一界面的折射縱波;③代表直達波;④代表第一界面的透射橫波;⑤代表第一界面的透射縱波;⑥代表第一界面的透射縱波在第二界面產(chǎn)生的透射縱波;⑦代表第一界面的透射縱波在第二界面產(chǎn)生的反射橫波;⑧代表第一界面的透射縱波在第二界面產(chǎn)生的透射橫波;⑨代表第一界面的透射縱波在第二界面產(chǎn)生的反射縱波;代表第一界面的透射縱波在第二界面產(chǎn)生的反射縱波，又入射到第一界面產(chǎn)生反射縱波;代表第一界面的透射橫波在第二界面產(chǎn)生的反射縱波;代表第一界面的透射橫波在第二界面產(chǎn)生的反射橫波;代表第一界面的透射橫波在第二界面產(chǎn)生的透射橫波;代表第一界面的透射橫波在第二界面產(chǎn)生的透射縱波。

圖3 1000ms時刻的波場快照a 三維層狀模型; b vx分量; c vy分量; d vz分量

vP/(m·s-1)vS/(m·s-1)h/m第1層15000750第2層200011541350第3層30001732

并行計算過程中，分析不同進程數(shù)以及不同進程分配方案的并行效率，通常是由加速比S來定義：

(6)

其中，WS為計算過程中的串行所用時間;WP為計算過程中的并行所用時間，包括并行計算時間和通信時間;P為計算所用的總進程個數(shù)。表3給出了不同計算進程的并行效率。

根據(jù)交錯網(wǎng)格有限差分格式，在計算速度分量的過程中，X方向需要交換的應力分量為σxx，σxy和σxz;Y方向需要交換的應力分量為σyy，σxy和σyz;Z方向需要交換的應力分量為σzz，σxz和σyz。在計算應力分量的過程中，X，Y，Z方向需要交換的速度分量均為vx，vy，vz。因此每個方向與相鄰節(jié)點都需要交換3個變量，同時結(jié)合交錯網(wǎng)格有限差分的空間差分精度N，則該模型正演過程中各個試驗單個計算節(jié)點與周圍相鄰的節(jié)點需要交換的最多數(shù)據(jù)量為(單位為3N)：當PxPyPz=2×1×1，單個進程數(shù)據(jù)交換量為480×360=172800;當PxPyPz=2×2×1，單個進程數(shù)據(jù)交換量為300×360+240×360=194400;當PxPyPz=2×2×2，單個進程數(shù)據(jù)交換量為300×240+300×180+240×180=169200;當PxPyPz=5×3×2，單個進程數(shù)據(jù)交換量為120×160+(120×180+160×180)×2=120000。

各個不同試驗的計算結(jié)果如表3所示，從表3可以看出：①隨著計算進程的增加，每個進程分配的計算網(wǎng)格點數(shù)減少，單個進程占用內(nèi)存減少。②隨著計算進程的增加，正演計算時間明顯減少，同時加速比隨著進程數(shù)的增加而增加。但是由于串行部分的存在以及交換進程數(shù)的增加帶來通信時間上的增加，故加速比增加，速度逐漸變小。③單個方向上并行，滿足每個進程的交換進程最少(最多為2個);3個方向上并行，滿足每個進程的交換信息數(shù)最少。

表3 不同計算進程的并行效率

3 海上某靶區(qū)OBC三維多波多分量地震觀測的數(shù)值模擬試驗

實際海上某靶區(qū)三維地質(zhì)模型如圖4所示，大小為4.5km×10.0km×6.6km。采集所用觀測系統(tǒng)如圖5所示，一共有8條電纜接收，每條電纜長度為7200m，有240道接收，道間距為30m，纜間距為400m，選取中間的一個震源點來激發(fā)。正演過程中空間網(wǎng)格大小為10m×10m×10m，時間步長為0.9ms，記錄長度為6s，為了提高正演精度，采用空間8階有限差分算子進行正演模擬。將正演任務分為60個進程并行計算，其中X，Y，Z方向分別采用的計算進程個數(shù)為Px=3，Py=5，Pz=4，正演所用時間為19.8h。圖6到圖9 為正演得到的多波多分量(4C)記錄。圖中直達波、折射波、轉(zhuǎn)換橫波以及深層海底反射波均能看到。由于接收點位置與震源位置有一定的距離，因此折射波先于直達波到達，尤其是在離震源位置較遠的電纜上這種情況更加明顯。除了獲得了海上4C分量模擬地震數(shù)據(jù)之外，同時還基于散度和旋度對彈性波場進行分離，得到了分離后的縱波和轉(zhuǎn)換橫波記錄，如圖10到圖11所示。可以看出，OBC采集不僅能夠接收縱波，還能接收到轉(zhuǎn)換橫波，分離后的縱波分量和水聽器分量記錄大致相同，說明了波場分離的正確性。

圖4 海上某靶區(qū)三維地質(zhì)模型

圖5 海上某靶區(qū)OBC地震采集觀測系統(tǒng)(紅色為震源)

圖6 海上某靶區(qū)OBC多波多分量(4C)模擬地震記錄vx分量

圖7 海上某靶區(qū)OBC多波多分量(4C)模擬地震記錄vy分量

圖8 海上某靶區(qū)OBC多波多分量(4C)模擬地震記錄vz分量

圖9 海上某靶區(qū)OBC多波多分量(4C)模擬地震記錄水聽器分量

圖10 海上某靶區(qū)OBC模擬資料波場分離后的縱波分量

圖11 海上某靶區(qū)OBC模擬資料波場分離后的橫波分量

4 結(jié)束語

交錯網(wǎng)格有限差分格式具有良好的局部性，非常有利于并行計算。基于MPI的彈性波有限差分正演模擬，將海量三維計算所需內(nèi)存分配到不同的進程上，減少了各個進程的內(nèi)存需求，而且每個進程計算三維模型的一部分區(qū)域，各個區(qū)域并行計算。從X,Y,Z 3個方向?qū)⒉⑿腥蝿談澐?，選擇更加靈活，能夠合理利用多個計算節(jié)點，提高了正演模擬效率?；贛PI的海上多波多分量地震觀測的正演模擬可以為研究海上寬方位采集、觀測系統(tǒng)論證以及高分辨率成像提供基礎數(shù)據(jù)。

參 考 文 獻

[1]CaldwellJ.Marinemulticomponentseismology[J].TheLeadingEdge,1999,18(11)：1274-1282

[2] 張樹林，夏斌，何家雄.海上多波多分量地震采集技術的應用[J].天然氣地球科學，2005，16(1):103-107

ZhangSL,XiaB,HeJX.Theoffshoremulti-wavesandmulti-componentsseismicacquisitiontechnique:acasestudyofYinggehaiBasin[J].NaturalGasGeoscience,2005,16(1):103-107

[3] 張樹林，李緒宣，姜立紅.海上多波多分量地震技術新進展與發(fā)展方向[J].物探化探計算技術，2000，22(2):97-107

ZhangSL,LiXX,JiangLH.ImprovementanddevelopmentofChinaoffshoremultiwaveandmulticomponentseismictechnique[J].ComputingTechniquesforGeophysicalandGeochemicalExploration,2000,22(2):97-107

[4] 劉海波，全海燕，陳浩林,等.海上多波多分量地震采集綜述[J].中國石油勘探，2007,12(3):52-57

LiuHB,QuanHY,ChenHL,etal.Overviewofmarinemulti-waveandmulti-componentseismicacquisition[J].ChinaPetroleumExploration,2007,12(3):52-57

[5]HeHY,ZhouHZ,FuDD.Multicomponentseismic-explorationinYinggehaiBasin[J].TheLeadingEdge,2002,21(9):914-920

[6]VirieuxJ.SHwavepropagationinheterogenousmedia:velocity-stressfinite-differencemethod[J].Geophysics,1984,49(11):1933-1957

[7]VirieuxJ.P-SVwavepropagationinheterogenousmedia:velocity-stressfinite-differencemethod[J].Geophysics,1986,51(4):889-901

[8]LevanderAR.Fourth-orderfinite-differenceP-SVseismograms[J].Geophysics,1988,53(11):1425-1436

[9] 董良國，馬在田，曹景忠，等.一階彈性波方程交錯網(wǎng)格高階差分解法[J].地球物理學報，2000,43(3)：411-419

DongLG,MaZT,CaoJZ,etal.Astaggered-gridhigh-orderdifferencemethodofone-orderelasticwaveequation[J].ChineseJournalofGeophysics,2000,43(3):411-419

[10] 裴正林.三維各向同性介質(zhì)彈性波方程交錯網(wǎng)格高階有限差分模擬[J].石油物探，2005,44(4)：308-315

PeiZL.Numericalsimulationofelasticwaveequationin3-Disotropicmediawithstaggered-gridhigh-orderdifferencemethod[J].GeophysicalProspectingforPetroleum,2005,44(4):308-315

[11] 裴正林.三維各向異性介質(zhì)中彈性波方程交錯網(wǎng)格高階有限差分法數(shù)值模擬[J].石油大學學報(自然科學版)，2004,28(5)：23-29

PeiZL.Three-dimensionalnumericalsimulationofelasticwavepropagationin3-Danisotropicmediawithstaggered-gridhigh-orderdifferencemethod[J].JournaloftheUniversityofPetroleum,China,2004,28(5):23-29

[12] 張軍華,臧勝濤,單聯(lián)瑜，等.高性能計算的發(fā)展現(xiàn)狀及趨勢[J].石油地球物理勘探,2010,45(6):918-925

ZhangJH,ZangST,ShanLY,etal.Developmentstatusandtrendsforhighperformancecomputing[J].OilGeophysicalProspectingforPetroleum,2010,45(6):918-925

[13] 都志輝.高性能計算并行編程技術—MPI并行程序設計[M].北京：清華大學出版社,2001：13-15

YuZH.Parallelprogramtechniqueofhighperformancecomputation—MPIparallelprogramdesign[M].Beijing:TsinghuaUniversityPress,2001：13-15

[14]FurumuraT,KoketsuK,WenKL.ParallelPSM/FDMhybridsimulationofgroundmotionsfromthe1999Chi-chi,Taiwan,earthquake[J].PureandAppliedGeophysics,2002,159(9):2133-2146

[15]BohlenT.Parallel3-Dviscoelasticfinitedifferenceseismicmodeling[J].Computers&Geosciences,2002,28(8):887-899

[16]MicikeviciusP.3DfinitedifferencecomputationonGPUsusingCUDA[J].ExpandedAbstractsof2ndWorkshoponGeneralPurposeProcessingonGraphicsProcessingUnits,2009,79-84

[17] 王德利,雍運動,韓立國，等.三維粘彈介質(zhì)地震波場有限差分并行模擬[J].西北地震學報,2007,29(1):30-34

WangDL,YongYD,HanLG,etal.Parallelsimulationoffinitedifferenceforseismicwavefiledmodelingin3-Dviscoelasticmedia[J].NorthwesternSeismologicalJournal,2007,29(1):30-34

[18] 王月英.基于MPI的三維波動方程有限元法并行正演模擬[J].石油物探,2009,48(3):221-243

WangYY.3Dwave-equationfinite-elementforwardmodelingbasedonmessagepassinginterfaceparallelcomputing[J].GeophysicalProspectingforPetroleum,2009,48(3):221-243

[19] 何兵壽,張會星,韓令賀.彈性波方程正演的粗粒度并行算法[J].地球物理學進展,2010,25(2):650-656

HeBS,ZhangHX,HanLH.Forwardmodelingofelasticwaveequationwithcoarse-grainedalgorithem[J].ProgressinGeophysic,2010,25(2):650-656

[20] 魏星,王彥賓,陳曉非.模擬地震波場的偽譜和高階有限差分混合方法[J].地震學報,2010,32(4):392-400

WeiX,WangYB,ChenXF.HybridPSM/FDMmethodforseismicwavefieldsimulation[J].ActaSeismologicaSinica,2010,32(4):392-400

[21] 秦艷芳,王彥賓.地震波傳播的三維偽譜和高階有限差分混合方法并行模擬[J].地震學報,2012,34(2):147-156

QinYF,WangYB.Three-dimensionalparallelhybridPSM/FDMsimulationofseismicwavepropagation[J].ActaSeismologicaSinica,2012,34(2):147-156

[22] 公續(xù)飛,杜啟振.三維彈性波有限差分數(shù)值模擬及其GPU并行實現(xiàn)[C]∥中國地球物理學會第二十七屆年會論文集.長沙：中國地球物理學會，2011：485

GongXF,DuQZ.Finite-differencesimulationof3DelasticwavefieldanditsimplementationonGPU[C]∥ChineseGeophysicalSociety27thAnnualMeetingProceeding.Changsha:ChineseGeophysicalSociety,2011:485

[23] 段玉婷，李靖宇，胡天悅.基于GPU的三維精細積分法正演模擬[C]∥中國地球物理學會第二十七屆年會論文集.長沙：中國地球物理學會，2011:484

DuanYT,LiJY,HuTY.Applicationofthe3DpreciseintegrationmethodforseismicmodelingbasedonGPU[C]∥ChineseGeophysicalSociety27thAnnualMeetingProceeding.Changsha:ChineseGeophysicalSociety,2011:484

[24] 張明財,熊章強,張大洲.基于MPI的三維瑞雷面波有限差分并行模擬[J].石油物探,2013,52(4):354-362

ZhangMC,XiongZQ,ZhangDZ.3DfinitedifferenceparallelsimulationofRayleighwavebasedonmessagepassinginterface[J].GeophysicalProspectingforPetroleum,2005,44(4):308-315

[25] 牟永光,裴正林.三維復雜介質(zhì)地震數(shù)值模擬[M].北京：石油工業(yè)出版社,2005:16-18

MouYG,PeiZL.Seismicnumericalmodelingfor3-Dcomplexmedia[M].Beijing:PetroleumIndustryPress,2005:16-18