李祎昕， 譚捍東

(1.中國地質(zhì)大學(xué) 地球物理與信息技術(shù)學(xué)院，北京 100083；2.遼寧省物測勘查院 ，沈陽 110121)

基于MPI多通訊域的地震波頻率域二維正演并行算法研究

李祎昕1,2， 譚捍東1

(1.中國地質(zhì)大學(xué) 地球物理與信息技術(shù)學(xué)院，北京 100083；2.遼寧省物測勘查院 ，沈陽 110121)

地震波場頻率域正演是頻率域全波形反演的基礎(chǔ)。針對反演計算量巨大的問題，利用頻率域二維波形正演算法中頻率和炮點(diǎn)計算的獨(dú)立性，開發(fā)出粗細(xì)粒度結(jié)合的MPI并行算法：粗粒度為頻率并行，細(xì)粒度為單個頻率解方程并行。實(shí)現(xiàn)方法是將正演頻率分組放入不同的MPI通訊域內(nèi)，每個通訊域內(nèi)單個頻率求解方程過程，采用基于MPI的MUMPS(多波前大規(guī)模并行稀疏直接解法器)軟件包并行加速。模型測試結(jié)果表明：MPI多通訊域并行算法計算結(jié)果正確，計算效率顯著提高，加速效果穩(wěn)定。

地震正演; 并行算法; MPI通訊域

0 引言

全波形方法在反演過程中可以使用地下波場的全部信息，是一種能夠揭示復(fù)雜地質(zhì)背景下構(gòu)造與儲層物性的方法。Tarantola[1]提出基于最小二乘法的二維時間域全波形反演方法，為全波形反演方法奠定基礎(chǔ)；Pratt等[2-3]將全波形反演推廣到頻率域，實(shí)現(xiàn)了高斯牛頓和全牛頓方法的非均勻速度模型反演，取得良好效果。國內(nèi)、外學(xué)者對全波形理論的深入研究表明，相比于時間域，頻率域全波形反演能更好地應(yīng)用于反演的層析成像，頻率域反演只需要反演有限個頻率點(diǎn)，可以提高計算效率。對于粘彈性介質(zhì)，頻率域正演模擬比時間域正演模擬更容易實(shí)現(xiàn)。頻率域正演模擬時,其衰減系數(shù)可以是頻率的函數(shù)，實(shí)現(xiàn)起來方便。另外，在數(shù)值計算過程中，頻率域全波形反演在多炮情況下計算更加高效，這是由于頻率域正演過程中直接法求解大型稀疏方程組時，不同炮(即方程組中的不同右端項(xiàng))使用同一LU分解因子，而時間域不同炮的波場必須分別求解。目前，頻率域正演數(shù)值模擬常用的方法有有限元法和有限差分法，與有限元法相比，有限差分方法計算時間短，結(jié)果穩(wěn)定，比較容易實(shí)現(xiàn)。但是與其他地震數(shù)值模擬方法相比，全波形方法數(shù)據(jù)量和計算量都很大，目前還難以得到廣泛應(yīng)用[4-14]。

隨著計算機(jī)技術(shù)地發(fā)展，以及對勘探復(fù)雜度要求地提高，并行計算逐漸在地球物理中得到應(yīng)用。并行計算是由運(yùn)行在多個部件上的小任務(wù)合作來求解一個規(guī)模很大的計算問題的一種方法。簡單的來說就是“分而治之”。與串行相比并行可以降低單個問題的求解時間，增加問題求解的規(guī)模，提高問題求解的精度。并行程序的開發(fā)需要特殊的編程環(huán)境，目前常用的并行編程環(huán)境有MPI、OpenMP、MPV、Cuda等。其中MPI(Message Passing Interface)是并行程序的開發(fā)中最常用的編程工具，它的功能十分強(qiáng)大，可以在多種平臺上應(yīng)用，而且計算速度快。

在全波形反演過程中，需要多次正演來擬合實(shí)際波場。正演模擬的精確程度和計算效率，直接影響了最終反演的效果，但是常規(guī)的全波形正演模擬是非常耗時的。因此，開發(fā)高效的地震波正演并行算法十分有必要。目前，地震波正演并行算法研究多在時間域內(nèi)進(jìn)行[15-18]，采用的并行方式多為將計算區(qū)域劃分為多個計算任務(wù)并賦給不同的進(jìn)程執(zhí)行，這種方式降低了內(nèi)存消耗，但各個進(jìn)程需頻繁通信來交換計算區(qū)域邊界上的計算結(jié)果而降低了并行計算效率。針對這個問題，這里開發(fā)了基于MPI的頻率和解方程兩重并行的頻率域二維波形正演并行計算方法。

1 頻率域二維波形正演理論

1.1 聲波波動方程

二維頻率域聲波方程在各向同性介質(zhì)下的波場傳播：

(1)

其中：ρ(x,z)代表密度；κ(x,z)為復(fù)體積模量；ω為頻率；p(x,z,ω)為壓力場；s(x,z,ω)為震源。衰減的影響可通過體積模量表達(dá)式中的復(fù)P波速度實(shí)現(xiàn)。

1.2 有限差分算法

應(yīng)用有限差分法求解波動方程，在構(gòu)造差分方程時采用混合網(wǎng)格(Themixed-gridstencil)9點(diǎn)差分格式，該模板聯(lián)合了經(jīng)典笛卡爾坐標(biāo)系中的交錯網(wǎng)格模板(staggered-gridstencil)和45°旋轉(zhuǎn)網(wǎng)格模板(45°rotated-gridstencil)[4]。

經(jīng)過推導(dǎo)，可以得到線性代數(shù)方程組(推導(dǎo)過程見附錄A)：

AP=S

(2)

其中：A為復(fù)數(shù)阻抗矩陣，包含了頻率、介質(zhì)速度、密度等信息；P為某一個頻率下的壓力場向量；S為某一個頻率下的震源向量。P和S在水平方向和垂直方向的網(wǎng)格維數(shù)分別為Nx和Nz。

1.3 吸收邊界

聲波波動方程模擬的是無限介質(zhì)，正演的模型是有限的，對邊界的截斷會造成界面反射，影響正演模擬精度。為了消除界面反射，采用PML(Perfectly Matched Layers)最佳匹配層吸收邊界。

(3)

(4)

圖1 最佳匹配層吸收邊界示意圖Fig.1 The sketch map of perfectly matched layers

1.4 稀疏矩陣方程求解

在求解稀疏矩陣的過程中調(diào)用了MUMPS(Multifrontal Massively Parallel Sparse Direct Solver)軟件包，MUMPS中文名是多波前大規(guī)模并行稀疏直接解法器，是基于直接法中的多波前方法開發(fā)出來的。多波前法是波前法的一種改進(jìn)方法，在分解策略上將大型稀疏矩陣轉(zhuǎn)化為很多小的稠密矩陣分解，減小了運(yùn)算量，能夠有效地求解大規(guī)模稀疏矩陣。使用時可以調(diào)用fortran和c的接口。

阻抗矩陣A是非對稱、非正定的，并且是一個高度稀疏的矩陣。調(diào)用MUMPS采用了直接法中的LU因式分解法，LU因式分解法是一種十分適合多炮正演的方法，可以在多炮正演中大大降低運(yùn)算量。因?yàn)槊恳粋€震源右端項(xiàng)S都可以重復(fù)使用LU因式分解結(jié)果，求解過程都是簡單的線性解方程。

LU[P1P2…Pn]=[S1S2…Sn]

(5)

MUMPS軟件包可以在串行環(huán)境中運(yùn)用，也可以在并行環(huán)境中運(yùn)用，在并行環(huán)境中運(yùn)用時支持MPI接口?？梢酝ㄟ^MPI調(diào)用多個進(jìn)程實(shí)現(xiàn)解方程的并行加速，我們是在并行環(huán)境中實(shí)現(xiàn)了解方程計算。

2 二維頻率域正演并行算法設(shè)計

2.1 設(shè)計思路和難點(diǎn)

從二維頻率域波形正演的原理可知，頻率域波形正演需要求解由低到高一系列頻率的波動方程的解，每個頻率是獨(dú)立計算的，不存在累積誤差。另外，在正演單個頻率時，多炮點(diǎn)求解方程組的過程中也可以調(diào)用基于MPI的并行解方程軟件求解。算法的設(shè)計意圖是將并行算法應(yīng)用到這兩種計算中，實(shí)現(xiàn)頻率和解方程的兩重并行。

在MPI并行計算中，進(jìn)程與進(jìn)程之間的通訊是比較費(fèi)時的，頻繁的通訊在路由器帶寬限制下極易造成網(wǎng)絡(luò)阻塞，導(dǎo)致大部分時間耗費(fèi)在進(jìn)程間相互通訊的等待上，不能發(fā)揮并行計算的優(yōu)勢。因此提高算法計算效率的關(guān)鍵是盡量減少進(jìn)程間的通訊。

在并行算法實(shí)現(xiàn)過程中遇到了一個難點(diǎn)：單個頻率求解方程過程中，解方程并行軟件會自動調(diào)用全部的進(jìn)程。所以求解第一個頻率用到的進(jìn)程不能同時求解第二個頻率，與我們的設(shè)計思路矛盾。因此，需要開發(fā)出一種算法對MPI的進(jìn)程進(jìn)行重新組織，以便各個頻率之間可以并行計算。

2.2 進(jìn)程組和通訊域

為了闡明MPI多通訊域算法首先要引入兩個概念，①M(fèi)PI進(jìn)程組；②通訊域。進(jìn)程組(process group)表示MPI的進(jìn)程的集合。通訊器(communicator)表示進(jìn)程之間的通訊環(huán)境。通訊器分為2類：①域內(nèi)通信器；(intracommunicator)②域間通信器(intercommunicator)。域內(nèi)通訊器表示同一個進(jìn)程組內(nèi)部之間的通訊；域間通信器表示不同進(jìn)程組之間的通訊。MPI在初始化時會默認(rèn)創(chuàng)建一個通訊器叫MPI_COMM_WORLD，所有的進(jìn)程都是在這個通訊器的內(nèi)部進(jìn)行通訊的。

多通訊域概念就是從這里來的，如圖2所示，將已有的進(jìn)程分成兩組group1和group2，基于這兩組創(chuàng)建兩個通訊器comm1和comm2，這樣單個通訊器內(nèi)部就可以實(shí)現(xiàn)不與另外通訊器發(fā)生關(guān)系的域內(nèi)通訊了。在程序之間通訊比較復(fù)雜時，將通訊域劃分為多個，可以更加方便、有序的管理進(jìn)程。

圖2 MPI進(jìn)程組和通訊器關(guān)系示意圖Fig.2 Relationships between MPI process group and MPI communicator

2.3 多通訊域算法設(shè)計

MPI多通訊域算法實(shí)現(xiàn)流程如下：

1)MPI初始化時定義新的通訊域。

2)修改輸入文件的參數(shù)把頻率分組。

3)將分好的頻率組放入通訊域中。

4)將通訊域賦給MUMPS實(shí)現(xiàn)解方程并行。

5)收集正演計算結(jié)果。

部分關(guān)鍵程序段如下：

MPI初始化時定義新的通訊域，例如將4個進(jìn)程分成兩個通訊域，通訊域1(P0、P1)通訊域2(P2、P3)。

####################

CALL MPI_INIT(infompi)! Define Host and Slave Processors CALL MPI_COMM_RANK(MPI_COMM_WORLD,mype,infompi)CALL MPI_COMM_SIZE(MPI_COMM_WORLD,nproc,infompi)

column=int(mype/2)

key=mype-column*2

CALL MPI_Comm_split(MPI_Comm_world,column,key,Comm_column)

CALL MPI_Comm_rank(Comm_column,mype_line,infompi)

CALL MPI_Comm_size(Comm_column,nproc_line,infompi)

CALL MPI_Barrier(MPI_COMM_WORLD,infompi)

###########################

修改輸入文件fwt.par的參數(shù)把頻率分組。

##################################

nw !正演頻率數(shù)

Nfreqgroup !頻率組數(shù)(與進(jìn)程分組對應(yīng))

iwg(nfreqgroup) !把頻率分段，每組的頻率劃分

###################################

將分好的頻率組放入通訊域中

###################################

select case (column)

case(0)

ifreqgroup=1

write(*,*) mype,mype_line,nproc_line

case(1)

ifreqgroup=2

write(*,*) mype,mype_line,nproc_line

end select

####################################

將通訊域賦給MUMPS

####################################

! INITIALIZE MUMPS AND ALLOCATE ARRAYS IRN, JCN AND A FOR MUMPSIF ( mype_line .eq. 0 ) THEN

WRITE(*,*) "Initialize MUMPS"

ENDIF

! Parallel sessionid%COMM = Comm_column

! Unsymmetric matrix

id%SYM = 0

! Host working

id%PAR = 1

! Initilization of MUMPS

id%JOB = -1

! Initialize MUMPS

CALL cmumps( id )

######################################

例如正演4個頻率f1、f2、f3、f4，運(yùn)行4個MPI進(jìn)程P0、P1、P2、P3。將4個進(jìn)程分成兩組，P0和P1為通訊域1，P2和P3為通訊域2，將4個頻率分成兩組放入兩個通訊域內(nèi)，在通訊域1內(nèi)，應(yīng)用頻率循環(huán)求解f1、f2，求解單個頻率的過程中調(diào)用MUMPS函數(shù)庫(進(jìn)程P0、P1)求解。同時，在通訊域2中調(diào)用MUMPS函數(shù)庫(進(jìn)程P2、P3)求解頻率f3、f4。最后再把正演計算結(jié)果收集起來(圖3)。

圖3 正演流程示意圖Fig.3 Flow chart of forward modeling

多通訊域代碼可以通過輸入文件靈活修改通訊域個數(shù)和頻率的分配，在進(jìn)程數(shù)和頻率數(shù)增加時，可以分配更多的通訊域進(jìn)行計算。通訊域個數(shù)越多，每個通訊域正演的頻率越少，進(jìn)程之間的通訊減少，效率提高。

3 模型算例及效率檢測

建立如圖4所示的速度模型，水平方向的網(wǎng)格數(shù)nx=801,垂直方向的網(wǎng)格數(shù)nz=801，網(wǎng)格間距h=40 m，介質(zhì)速度v1=4 000 m/s，v2=6 000 m/s。觀測系統(tǒng)炮點(diǎn)在地下5 000 m處激發(fā)，炮間距為1 000 m，接收點(diǎn)在地下4 500 m處，道間距160 m。選取由低到高8個頻率進(jìn)行正演。(f1=0.393 700 8 Hz，f2=0.787 401 6，f3=1.181 102，f4=1.574 803，f5=1.968 504，f6=2.362 20，f7=2.755 90，f8=3.149 60)

采用MPI集群下配置相同的兩個節(jié)點(diǎn)進(jìn)行測試。兩個節(jié)點(diǎn)的配置為：CPU:Intel(R) Core(TM)i7950@ 3.07 GHz 4核8線程，內(nèi)存8 G。

圖4 速度模型fvFig.4 Velocity model fv

3.1 單炮測試

統(tǒng)計了單炮下8個進(jìn)程下開1、2、4、8個通訊域

正演8個頻率(1個通訊域相當(dāng)于頻率循環(huán)8次，8個通訊域相當(dāng)于8個頻率并行計算)需要的時間，并且計算了相對加速比和效率(表1)。由表1可以看出，隨著通訊域個數(shù)的增加，相對加速比和效率都提高。這是因?yàn)閯澐至舜蟮耐ㄓ嵙６戎鬁p少了通訊損耗的時間。

表1 多機(jī)測試中不同通訊域下相對加速比和效率Tab.1 Speedup and efficiency of different MPI communicators in two computers test

3.2 多炮測試

改變模型接收文件的炮數(shù)為1炮、5炮、10炮、20炮，統(tǒng)計模型(8進(jìn)程8個頻率)在1、2、4、8通訊域下的計算時間(表2)。由表2可以看出，在多炮計算中多通訊域的代碼加速效果穩(wěn)定。

3.3 正演結(jié)果繪圖

為了證明多通訊域代碼的正確性，繪制了非并行代碼和并行的兩通訊域代碼在f1=0.393 700 8 Hz和f2=0.787 401 6 Hz時頻率域正演波場快照fmap1和fmap2。圖5中可以看出，fmap1和fmap2。波形相同，說明多通訊域代碼的計算結(jié)果是可靠的。

表2 不同通訊域下1炮、5炮、10炮、20炮的模型計算時間Tab.2 The compute time in different MPI communicators of 1-shot，5-shots，10-shots，20-shots models

4 結(jié)論

1)對比并行和非并行算法的頻率域二維正演模擬算法，兩者計算結(jié)果一致，但并行算法利用了多計算核的優(yōu)勢，極大地提高了計算效率。

2)二維頻率域波形正演并行算法在單機(jī)測試，多機(jī)測試，多炮計算中表現(xiàn)都很穩(wěn)定，加速比顯著提高，這是由于多通訊域劃分的并行方式改變了計算粒度，進(jìn)程之間的通訊開銷降低了。

3)多通訊域代碼可以根據(jù)硬件條件和頻率個數(shù)進(jìn)行通訊域的劃分，具有計算的靈活性。

4)在設(shè)計通訊域個數(shù)時需要考慮內(nèi)存的開銷，因?yàn)槊總€通訊域內(nèi)都需要進(jìn)行獨(dú)立的并行解方程計算。也就是說N個通訊域在計算過程中需要同時存儲N個稀疏矩陣，這就增加了N倍內(nèi)存的開銷。

圖5 正演結(jié)果對比圖Fig.5 The comparison of forward modeling results(a)非并行程序頻率域波場快照fmap1；(b)非并行程序頻率域波場快照fmap2； (c)兩通訊域頻率域波場快照fmap1；(d)兩通訊域頻率域波場快照fmap2

[1] TARANTOLA A. Inversion of seismic reflection data in the acoustic approximation[J]. Geophysics, 1984, 49(8): 1259-1266.

[2] PRATT R. G, SHIN, C，HICKS,G J. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J]. Geophys. J. Int., 1998,133:341-362.

[3] PRATT R.G.. Seismic waveform inversion in the frequency domain, part 1:theory and verification in a physical scale model[J]. Geophysics, 1999，64:888-901.

[4] HUSTEDT B., OPERTO S.， VIRIEUX J. Mixed-grid and staggered-grid finite-difference methods for frequency-domain acoustic wave modelling[J]. Geophysical Journal International,2004, 157:1269-1296.

[5] RAVAUT C.,OPERTO S.,IMPROTA L.,et al.Herrero and P.Dell'Aversana Multiscale imaging of complex structures from multifold wide-aperture seismic data by ferquency-domain full-waveform tomography:application to a thrust belt[J].Geophysical Journal International, 2004，159:1032-1056.

[6] JO C. H., SHIN C， SUH, J. H. An optimal 9-point, finite-difference, frequency-space, 2D scalar wave extrapolator[J]. Geophysics,1996，61(2): 529-537.

[7] MARFURT, K. Accuracy of finite-difference and finite-elements modeling inscalar and elastic wave equation[J]. Geophysics, 1984，49: 533-549.

[8] SAENGER E. H., GOLD N. SHAPIRO A. Modeling the propagation of elastic waves using a modified finite-difference grid[J]. Wave motion, 2000，31:77-92.

[9] SHIN C., YOON K., MARFURT K.J. Efficient calculation of a partial derivative wavefield using reciprocity for seismic imaging and inversion[J]. Geophysics, 2001，66(6): 1856-1863.

[10]SIRGUE L.PRATT R. G..Efficient waveform inversion and imaging: a strategy for selecting temporal frequencies[J]. Geophysics, 2004. 69: 231-248.

[11]OPERTO S, RAVAUT C, IMPROTA L,et al.Quantitative imaging of complex structures from dense wide-aperture seismic data by multiscale traveltime and waveform inversions: a case study[J]. Geophysical Prospecting, 2004.52: 625-651.

[12]STEKL I，PRATT R.G.Accurate viscoelastic modeling by frequency-domain finite differences using rotated operators[J].Geophysics,1998，63(5):1779-1794.

[13]曹書紅, 陳景波 .聲波方程頻率域高精度正演的17點(diǎn)格式及數(shù)值實(shí)現(xiàn)[J].地球物理學(xué)報,2012(10): 3440-3449. CAO S H, CHEN J B .A 17-point scheme and its numerical implementation for high-accuracy modeling of frequency-domain acoustic equation[J].Chinese Journal Geophysics,2012(10): 3440-3449.(In Chinese)

[14]曹健, 陳景波, 曹書紅 .頻率域波動方程正演基于多重網(wǎng)格預(yù)條件的迭代算法研究[J].地球物理學(xué)報,2015,58(3): 1002-1012. CAO J, CHEN J B, CAO S H .Studies on iterative algorithms for modeling of frequency-domain wave equation based on multi-grid precondition[J].Chinese Journal Geophysics,2015,58(3): 1002-1012.(In Chinese).

[15]宋鵬,解闖,李金山,等，基于MPI+OpenMP的三維聲波方程正演模擬[J]. 中國海洋大學(xué)學(xué)報(自然科學(xué)版),2015(09):97-102+129. SONG P, XIE C, LI J S，et al .The 3D modeling of acoustic wave equation based on MPI+OpenMp[J].Periodical of Ocean University of China, 2015(09):97-102+129. (In Chinese)

[16]王月英.基于MPI的三維波動方程有限元法并行正演模擬[J]. 石油物探,2009,48(3):221-225. WANG Y Y.3D wave-equation finite-element forward modeling based on message passing interface parallel computing[J]. Geophysical Prospecting for Petroleum,2009,48(3):221-225.(In Chinese)

[17]何兵壽, 張會星, 韓令賀 .彈性波方程正演的粗粒度并行算法[J].地球物理學(xué)進(jìn)展,2010，25(2):650-656. HE B S, ZHANG H X, HAN L H.Forward modelling of elastic wave equation with coarse-grained parallel algorithm[J].Progress in Geophysics,2010，25(2): 650-656.(In Chinese)

[18]張明財，熊章強(qiáng)，張大洲. 基于MPI的三維瑞雷面波有限差分并行模擬[J]. 石油物探, 2013, 52(4): 354-362. ZHANG M C, XIONG Z Q, ZHANG D Z. 3D finite difference parallel simulation of Rayleigh wave based on message passing interface[J]. Geophysical prospecting for petroleum, 2013, 52(4): 354-362. (In Chinese)

Seismic frequency-domain 2D forward modeling based on parallel algorithms

LI Yixin1,2， TAN Handong1

(1.School of Geophysics and Information Technology China University of Geosciences，Beijing 100083,China;2.Geophysical prospecting exploration institute of Liaoning province，Shenyang 110121,China)

Frequency-domain forward is the basis of frequency-domain full waveform inversion. In order to solve the computationally intensive in inversion, we apply a combination of coarse-fined grained MPI parallel algorithm based on frequencies and multi-source simulations computing independently in frequency-domain forward modeling. We use coarse grained in frequency parallel and fine grained in solve equations parallel. To realization the method, we put different frequency group of forward modeling in different MPI communicators to make the frequencies computing independently. In each communicator, we use MUMPS (Multi-frontal Massively Parallel Sparse Direct Solver) based on MPI to speedup the method of solving equations. The MPI parallel algorithm can provide reference to other geophysical methods of the multi-source and multi-frequency forward computational efficiency.

seismic forward modeling; parallel algorithm; MPI communicators

2016-01-12 改回日期：2016-02-23

國家自然科學(xué)基金項(xiàng)目(41374078)；國土資源部地質(zhì)調(diào)查項(xiàng)目(12120113086100，12120113101300)

李祎昕(1988-),女，碩士，主要從事地球物理地震算法研究，E-mail：875298227@qq.com。

譚捍東(1966-)，男，教授，主要從事電法勘探理論及應(yīng)用研究，E-mail：thd@cugb.edu.cn。

1001-1749(2017)01-0052-07

P 631.4

A

10.3969/j.issn.1001-1749.2017.01.08