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

李祎昕， 譚捍東

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

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

李祎昕1,2， 譚捍東1

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

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

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

0 引言

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

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

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

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

1.1 聲波波動方程

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

(1)

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

1.2 有限差分算法

應用有限差分法求解波動方程，在構造差分方程時采用混合網格(Themixed-gridstencil)9點差分格式，該模板聯合了經典笛卡爾坐標系中的交錯網格模板(staggered-gridstencil)和45°旋轉網格模板(45°rotated-gridstencil)[4]。

經過推導，可以得到線性代數方程組(推導過程見附錄A)：

AP=S

(2)

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

1.3 吸收邊界

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

(3)

(4)

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

1.4 稀疏矩陣方程求解

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

阻抗矩陣A是非對稱、非正定的，并且是一個高度稀疏的矩陣。調用MUMPS采用了直接法中的LU因式分解法，LU因式分解法是一種十分適合多炮正演的方法，可以在多炮正演中大大降低運算量。因為每一個震源右端項S都可以重復使用LU因式分解結果，求解過程都是簡單的線性解方程。

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

(5)

MUMPS軟件包可以在串行環境中運用，也可以在并行環境中運用，在并行環境中運用時支持MPI接口。可以通過MPI調用多個進程實現解方程的并行加速，我們是在并行環境中實現了解方程計算。

2 二維頻率域正演并行算法設計

2.1 設計思路和難點

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

在MPI并行計算中，進程與進程之間的通訊是比較費時的，頻繁的通訊在路由器帶寬限制下極易造成網絡阻塞，導致大部分時間耗費在進程間相互通訊的等待上，不能發揮并行計算的優勢。因此提高算法計算效率的關鍵是盡量減少進程間的通訊。

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

2.2 進程組和通訊域

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

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

圖2 MPI進程組和通訊器關系示意圖Fig.2 Relationships between MPI process group and MPI communicator

2.3 多通訊域算法設計

MPI多通訊域算法實現流程如下：

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

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

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

4)將通訊域賦給MUMPS實現解方程并行。

5)收集正演計算結果。

部分關鍵程序段如下：

MPI初始化時定義新的通訊域，例如將4個進程分成兩個通訊域，通訊域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的參數把頻率分組。

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

nw !正演頻率數

Nfreqgroup !頻率組數(與進程分組對應)

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，運行4個MPI進程P0、P1、P2、P3。將4個進程分成兩組，P0和P1為通訊域1，P2和P3為通訊域2，將4個頻率分成兩組放入兩個通訊域內，在通訊域1內，應用頻率循環求解f1、f2，求解單個頻率的過程中調用MUMPS函數庫(進程P0、P1)求解。同時，在通訊域2中調用MUMPS函數庫(進程P2、P3)求解頻率f3、f4。最后再把正演計算結果收集起來(圖3)。

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

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

3 模型算例及效率檢測

建立如圖4所示的速度模型，水平方向的網格數nx=801,垂直方向的網格數nz=801，網格間距h=40 m，介質速度v1=4 000 m/s，v2=6 000 m/s。觀測系統炮點在地下5 000 m處激發，炮間距為1 000 m，接收點在地下4 500 m處，道間距160 m。選取由低到高8個頻率進行正演。(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集群下配置相同的兩個節點進行測試。兩個節點的配置為：CPU:Intel(R) Core(TM)i7950@ 3.07 GHz 4核8線程，內存8 G。

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

3.1 單炮測試

統計了單炮下8個進程下開1、2、4、8個通訊域

正演8個頻率(1個通訊域相當于頻率循環8次，8個通訊域相當于8個頻率并行計算)需要的時間，并且計算了相對加速比和效率(表1)。由表1可以看出，隨著通訊域個數的增加，相對加速比和效率都提高。這是因為劃分了大的通訊粒度之后減少了通訊損耗的時間。

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

3.2 多炮測試

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

3.3 正演結果繪圖

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

表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 結論

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

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

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

4)在設計通訊域個數時需要考慮內存的開銷，因為每個通訊域內都需要進行獨立的并行解方程計算。也就是說N個通訊域在計算過程中需要同時存儲N個稀疏矩陣，這就增加了N倍內存的開銷。

圖5 正演結果對比圖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點格式及數值實現[J].地球物理學報,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]曹健, 陳景波, 曹書紅 .頻率域波動方程正演基于多重網格預條件的迭代算法研究[J].地球物理學報,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]. 中國海洋大學學報(自然科學版),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].地球物理學進展,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]張明財，熊章強，張大洲. 基于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

國家自然科學基金項目(41374078)；國土資源部地質調查項目(12120113086100，12120113101300)

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

譚捍東(1966-)，男，教授，主要從事電法勘探理論及應用研究，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