基于MPI的OBC三維多波多分量地震觀測正演模擬并行算法實(shí)現(xiàn)

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

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

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

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

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

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

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

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

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

(1a)

(1b)

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

(2)

(3)

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

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

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

(4)

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

(5a)

(5b)

(5c)

(5d)

2 基于MPI的三維并行計(jì)算的整體設(shè)計(jì)

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

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

2.1 三維并行相鄰進(jìn)程的交換設(shè)計(jì)

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

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

圖2 三維相鄰進(jìn)程交換示意圖解

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

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

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

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

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

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

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

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

6) 同步驟4)一樣，進(jìn)行應(yīng)力場的計(jì)算。

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

2.3 并行效率分析

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

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

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

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

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

(6)

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

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

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

表3 不同計(jì)算進(jìn)程的并行效率

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

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

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

參 考 文 獻(xiàn)

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

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

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

[3] 張樹林，李緒宣，姜立紅.海上多波多分量地震技術(shù)新進(jìn)展與發(fā)展方向[J].物探化探計(jì)算技術(shù)，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] 董良國，馬在田，曹景忠，等.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法[J].地球物理學(xué)報(bào)，2000,43(3)：411-419

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

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

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

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

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

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

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

[13] 都志輝.高性能計(jì)算并行編程技術(shù)—MPI并行程序設(shè)計(jì)[M].北京：清華大學(xué)出版社,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] 王德利,雍運(yùn)動(dòng),韓立國，等.三維粘彈介質(zhì)地震波場有限差分并行模擬[J].西北地震學(xué)報(bào),2007,29(1):30-34

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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