基于MPI+OpenMP的時間域航空電磁快速正演算法

任運通，李貅，齊彥福，曹華科

(長安大學 地質(zhì)工程與測繪學院，陜西 西安 710054)

0 引言

時間域航空電磁法采用機載移動平臺，通過線圈發(fā)射大功率一次場電磁波信號，利用接收線圈接收到感應的二次磁場，實現(xiàn)對地下電性結(jié)構(gòu)及能源礦產(chǎn)分布的精細探測。該方法具有分辨率高、探測速度快、通行性極好等諸多優(yōu)點,已被廣泛應用于油氣資源、環(huán)境工程和金屬礦產(chǎn)等領域[1-6]。然而，傳統(tǒng)的電阻率成像和一維反演技術(shù)已經(jīng)無法滿足三維問題解釋的要求，開展三維反演解釋技術(shù)研究是當前的熱點問題，作為反演的基礎，時間域航空電磁三維正演算法受到越來越多的關注。

目前，時間域航空電磁三維有限元正演方法主要分為間接法和直接法兩種。其中的間接法是基于頻率域麥克斯韋方程組，首先計算頻率域響應，再通過時頻轉(zhuǎn)換方法轉(zhuǎn)換到時間域[7-10]。由于間接算法存在穩(wěn)定性差等諸多問題，人們普遍將研究方向轉(zhuǎn)移到了時間域直接算法，直接求解時間域電磁響應。2010 年，Um等[11]采用非結(jié)構(gòu)化網(wǎng)格剖分策略，對電性源的海洋三維瞬變電磁場進行了時間域矢量有限元的三維正演模擬計算。2016年，李賀等[12]人研究了直接時間域矢量有限元瞬變電磁三維正演模擬方法。2017年，齊彥福等人[13-14]基于航空電磁footprint的局部網(wǎng)格策略，利用非結(jié)構(gòu)矢量有限元方法，對復雜介質(zhì)情況下的時間域三維航空電磁響應進行了進一步模擬研究。然而，由于航空電磁法采用同時移動發(fā)射源和測點的方式進行觀測，且采樣密集，每次移動測點均需求解一次正演方程，產(chǎn)生巨大的計算量，傳統(tǒng)的串行計算方法無法滿足未來三維反演的效率要求?？紤]到并行計算可以有效提高程序的運算效率，解決相同數(shù)量的問題所需執(zhí)行的時間更短，開展了基于MPI+OpenMP并行技術(shù)的時間域航空電磁快速正演算法研究。

目前廣為使用的并行加速技術(shù)有兩種，分別是MPI(Message Passing Interface)和OpenMP(Open Multi-Processing)并行技術(shù)。MPI和OpenMP均需要與計算機語言結(jié)合使用，因兩者較好的通信性和可移植性而受到廣泛地應用。在尺度上，MPI多為基于多臺處理機上的跨節(jié)點并行方法，OpenMP主要通過一些指令集對現(xiàn)有的Fortran或C/C++等程序進行擴展。并行技術(shù)在提高計算效率上具有顯著優(yōu)點，將該方法應用到地球物理數(shù)值模擬和數(shù)據(jù)處理中具有良好的前景。早在1997年，Newman和Alumbaugh[15]對三維大地電磁正演方法以及靈敏度矩陣進行了深入研究，成功將并行加速技術(shù)應用到了電磁反演當中，實現(xiàn)了多個節(jié)點的并行加速計算，極大推動了并行計算技術(shù)在地球物理領域的發(fā)展。國內(nèi)對于并行計算技術(shù)應用到地球物理的研究，起步雖然稍晚于國外，發(fā)展卻并不遜色。2005年榮瑩等[16]實現(xiàn)了基于MPI的處理機集群并行計算系統(tǒng)平臺的構(gòu)建，在現(xiàn)有的硬件條件下提高計算力； 2006年，譚捍東[17]等結(jié)合MPI的優(yōu)越性，通過頻點并行的方式成功實現(xiàn)了大地電磁三維正演的并行計算，得到了很高的加速比，驗證了并行算法的穩(wěn)定性以及高效性；2011年，李小康[18]研究了頻率域航空電磁法有限單元二維正演的并行計算；2015年，陳輝[19]研究了基于MPI的航空瞬變電磁一維正反演，提高了航空電磁一維模擬技術(shù)的計算效率。目前，對于并行加速技術(shù)在電磁方法數(shù)值模擬中的應用主要集中在大地電磁法，而航空電磁正演的并行算法研究依然停留在一維正反演和二維的頻率域正演階段。

本文將并行技術(shù)與時間域有限元算法相結(jié)合，基于三維航空電磁局部網(wǎng)格之間相互獨立的特性[20]，通過采用MPI+OpenMP的并行加速策略，在保證正演結(jié)果可靠性的基礎上，極大提高了時間域航空電磁有限元三維正演的計算效率。首先通過與串行有限元法和有限體積法進行對比檢驗本文并行程序的可靠性，然后進行計算效率分析，討論并行加速比與硬件條件的關系，最后將該并行算法應用于復雜起伏地表模型，模擬其航空電磁響應。

1 方法理論

1.1 時間域航空電磁有限元正演理論

時間域麥克斯韋方程組可以表示為[21]

(1)

(2)

(3)

(4)

其中，t為時間，r為位置矢量，e(r,t)、j(r,t)、d(r,t)及b(r,t)分別表示r處在t時刻電場強度、電流密度、電位移矢量以及磁感應強度，q和μ分別表示累積電荷和磁導率。本文假定介質(zhì)的磁導率與自由空間磁導率相同，即μ=μ0=4π×10-7H/m，介電常數(shù)與自由空間介電常數(shù)ε0相同，即ε=ε0=8.85×10-12F/m。電磁場之間具有如下本構(gòu)關系：

j(r,t)=σe(r,t)+js(r,t) ，

(5)

d(r,t)=εe(r,t) ，

(6)

b(r,t)=μh(r,t) ，

(7)

其中：js(r,t)表示外部施加的源電流密度，h(r,t)和σ分別表示磁場強度和電導率。通過消去磁場，可以獲得電場擴散方程

(8)

本文采用非結(jié)構(gòu)矢量有限元方法進行空間離散，四面體單元中任意位置的電場可以表示為：

(9)

(10)

(11)

其中：M、S分別為質(zhì)量和剛度矩陣，J為電流源項，針對每個單元，可由下面三式給出

(12)

(13)

(14)

其中V是單元的體積。由于各個單元相互獨立，因此Mk、Sk和Jk可以采用OpenMP并行技術(shù)加速計算。針對發(fā)射源，將發(fā)射線圈分解為若干段導線，每段近似為一個電偶極子[24]每個電偶極子可以表示為：

js(r,t)=δ(r-rs)vI(t)dl。

(15)

上式中：I(t)表示t時刻的電流強度，dl表示電偶極子長度，v表示電流方向，δ表示脈沖函數(shù)。然后，利用二階后推歐拉公式對有限元方程(11)進行時間域離散：

(3A+2ΔtB)ei+2(t)=A[4ei+1(t)-ei(t)]-2ΔtSi+2，

(16)

其中：Δt為時間步長，ei(t)表示第i時刻的電場值。式(16)亦可簡寫為：

Fe=P，

(17)

其中：F表示大型稀疏系數(shù)矩陣，P是右端項。當發(fā)射階躍電流波形時，在0時刻之前發(fā)射線圈中供恒定電流，在全空間產(chǎn)生穩(wěn)定的磁場，根據(jù)楞次定律可知空間中任意位置的電場均為0，故電場的初始條件為：

e(r,0)=0 。

(18)

采用狄利克雷邊界條件，根據(jù)電磁波在空間中的幾何指數(shù)衰減規(guī)律，可以通過加大計算區(qū)域的擴邊范圍，使其滿足

e|Γ=0 ，

(19)

其中Γ表示計算區(qū)域的外邊界。最后，采用直接求解器Pardiso對線性方程組進行求解，即可獲得全空間的電場值，再利用法拉第電磁感應定律(式(1))計算磁場響應。

1.2 局部網(wǎng)格正演計算技術(shù)

大量的數(shù)值實驗表明三維數(shù)值模擬結(jié)果的精度受到網(wǎng)格質(zhì)量的影響非常大，因此需要對導電大地和發(fā)射源進行精細的網(wǎng)格剖分來獲得高精度正演結(jié)果。傳統(tǒng)的時間域航空電磁三維正演模擬算法采用全局網(wǎng)格加密方式，即通過一套網(wǎng)格對所有測點的電磁響應進行模擬。然而，由于航空電磁法采樣密集且觀測剖面長，導致正演網(wǎng)格單元數(shù)量巨大。如果采用全局網(wǎng)格進行正演模擬，將產(chǎn)生巨大的計算量，嚴重降低計算效率?？紤]到航空電磁系統(tǒng)影響范圍有限(如圖1所示)，本文采用局部網(wǎng)格加密技術(shù)，針對每個測點或者相鄰的幾個觀測點分別設計獨立的網(wǎng)格，并分別在局部網(wǎng)格上進行正演模擬，以此在保證結(jié)果可靠性的同時提高計算效率。

圖1 航空系統(tǒng)影響范圍示意Fig.1 Moving footprint for time-domian airborne EM system

局部網(wǎng)格技術(shù)根據(jù)航空電磁系統(tǒng)的影響范圍設計合理的局部精細網(wǎng)格。圖2展示了全局網(wǎng)格和局部網(wǎng)格，其中圖2a表示全局網(wǎng)格，圖2b、c、d則分別表示針對測線上不同測點所設計的相互獨立的局部網(wǎng)格。從圖中可知全局網(wǎng)格對測線上所有測點、測線下方地空分界面以及異常體部分均進行了相對細致的網(wǎng)格剖分，因此造成網(wǎng)格數(shù)量較大，導致正演模擬的計算量劇增，計算效率降低。而局部網(wǎng)格僅僅對異常體、當前測點位置以及其下方的地空分界面進行局部加密，對單個測點進行單獨求解，因此所使用的網(wǎng)格數(shù)量大大降低，計算量減小。

2 并行加速

2.1 基于OpenMP的并行加速

OpenMP采用分叉—合并(Fork-Join)執(zhí)行模式。一個OpenMP程序開始于一個單一的線程，該進程又稱作主線程，通過并行指令對程序的并行區(qū)間進行定義，在這個區(qū)間中程序塊由多個線程自動分配任務及并行執(zhí)行，線程數(shù)由主線程決定，在一個大型程序中可以嵌套多個OpenMP并行區(qū)間，其中并行區(qū)間由一組$OMP指令進行開啟(!OMP PARALLEL)和關閉(!OMP END PARALLEL)。圖3是OpenMP并行模式的示意圖。

2.2 基于MPI的并行加速

MPI需要用頭文件use mpi進行聲明，然后在主進程(即節(jié)點)中進行初始化，開啟MPI并行環(huán)境，獲取進程個數(shù)及編號，再將多個任務合理地分配到各個進程并傳輸對應的數(shù)據(jù)，然后分別進行計算，最后匯集到主進程并關閉并行環(huán)境，如圖4所示。

圖2 全局網(wǎng)格與局部網(wǎng)格對比 (a為全局網(wǎng)格；b、c、d為局部網(wǎng)格)Fig.2 Global grid vs. local grid ((a) is a Global grid; (b) (c) (d) is a local grid)

圖3 OpenMP并行模式示意Fig.3 OpenMP parallel mode schematic

2.3 并行程序設計及實現(xiàn)

由于航空電磁系統(tǒng)在每個觀測點處的靈敏區(qū)域相對于總計算區(qū)域來說遠遠不及，因此僅僅在觀測點的一定位置范圍內(nèi)對網(wǎng)格進行加密即可滿足計算的精度要求，針對每個觀測點分別單獨設計獨立的網(wǎng)格來提高計算效率。

時間域航空電磁三維正演各個測點計算相互獨立，不存在數(shù)據(jù)依賴關系，具有非常好的并行性，在三維正演的基礎上實現(xiàn)了基于MPI的多測點并行計算，在獨立網(wǎng)格內(nèi)實現(xiàn)了基于OpenMP的單元矩陣并行計算。程序采用主從模式，整個程序的基本結(jié)構(gòu)及運行流程由主進程把控，主要負責讀取模型參數(shù)、數(shù)據(jù)傳輸以及對并行任務的分配，子進程負責接收來自主進程傳遞的消息、對分配的任務進行計算以及輸出本進程的計算結(jié)果。為了充分利用計算資源提高計算效率，主進程在完成任務分配和數(shù)據(jù)傳輸?shù)墓ぷ骱笠矃⑴c到測點的并行計算當中。

具體實現(xiàn)的流程如圖5所示，主進程讀取模型的整套網(wǎng)格參數(shù)，隨后將其傳遞給其他子進程，之后各進程根據(jù)分配的網(wǎng)格參數(shù)計算模型的響應結(jié)果，為了防止數(shù)據(jù)傳輸時發(fā)生通信沖突，采用各進程分別輸出各自的計算結(jié)果。

圖4 MPI并行模式示意Fig.4 MPI parallel mode schematic

圖5 航空瞬變電磁三維正演并行計算流程Fig.5 Airborne transient electromagnetic three-dimensional Forward modeling parallel computing flow chart

3 模型算例

3.1 精度驗證

為檢驗本文所開發(fā)的時間域航空電磁多粒度并行正演方法的正確性，設計了如圖6所示的水平板狀體模型，發(fā)射波形為階躍波，發(fā)射電流強度為435 A，通過全局網(wǎng)格和局部網(wǎng)格兩種剖分方式，對異常體上方y(tǒng)=0 m剖面上-400～400 m范圍內(nèi)的32個測點，分別進行正演模擬，且將正演結(jié)果與有限體積方法的正演結(jié)果進行比對，進行數(shù)值精度驗證。

圖6 水平板狀體模型示意Fig.6 Schematic diagram of the horizontal plate model

使用的集群服務器有16個型號為Intel(R) Xeon(R) CPU E5-2609 v4 @ 1.70 GHz的邏輯CPU，每個CPU擁有2G內(nèi)存和4個線程。采用時間域有限元全局網(wǎng)格法、基于多粒度并行加速策略的時間域有限元局部網(wǎng)格法以及有限體積方法分別計算的電磁響應結(jié)果如圖7所示，可以清楚地看出，上述3種不同方法所計算的結(jié)果吻合得非常好，且均對模型中的地下板狀體有明顯反映。這一結(jié)果有效驗證了本文并行算法的精度。

3.2 并行性能分析

為了分析基于MPI+OpenMP并行加速策略的時間域航空電磁三維正演程序的并行加速效果，采用上述水平塊狀體模型進行測算，根據(jù)控制變量法，控制計算測點數(shù)為32和線程數(shù)為4保持不變，依次對不同數(shù)量的進程進行比較，通過加速比和并行效率兩個參數(shù)來評估本文設計的三維正演并行程序。

令原有的串行程序在集群服務器單進程上的運行時間為Ts，經(jīng)過多粒度并行優(yōu)化后，運行所需時間為Tp，其中P表示開啟的進程數(shù)，則加速比Sp可表示為：

Sp=Ts/Tp，

(20)

并行效率可定義為：

圖7 全局網(wǎng)格串行計算和局部網(wǎng)格并行計算電磁響應結(jié)果對比Fig.7 Comparison of electromagnetic response results between global grid serial computing and local grid parallel computing

ep=Sp/P。

(21)

測試結(jié)果見表1。在計算測點數(shù)一定時，隨著進程數(shù)的增加，計算時間呈非線性下降趨勢，加速比逐漸增大，相反并行效率逐漸降低；當開啟16個進程時，采用多粒度并行加速后執(zhí)行耗費的時間僅僅是串行程序耗時的1/10，同時并行效率達到最低。

表1 基于MPI+OpenMP并行加速策略的行正演計算統(tǒng)計Table 1 Line forward calculation statistics based on MPI+OpenMP parallel acceleration strategy

由圖8可以發(fā)現(xiàn)以下特點：并行化后的正演程序加速比與開啟的進程數(shù)并不是整數(shù)倍關系，而是呈非線性增長，且并行效率整體呈降低趨勢。這是因為在集群服務器上開啟并行環(huán)境后，會占用一定內(nèi)存，且進程間進行數(shù)據(jù)傳遞需要占用時間，而且任務量與開啟的進程數(shù)不匹配時，計算任務少的進程會率先完成計算量，但該進程并不會結(jié)束，而是要等待未完成任務的進程，從而導致時間的損耗。所以開啟P個進程，加速比Sp并不能達到P，且隨著進程數(shù)的增加Sp與P的差值越來越大。

圖8 不同進程數(shù)的加速比(a)和并行效率(b)Fig.8 Acceleration ratio (a) and parallel efficiency (b) for different node numbers

3.3 復雜模型計算

考慮到地形起伏不平對航空電磁正演結(jié)果的影響較為嚴重，采取四面體網(wǎng)格進行加密剖分后的計算量較大，從而導致計算時間較長，因此利用多粒度并行加速技術(shù)進行加速計算。導入地形文件數(shù)據(jù)模擬起伏地表，并設計地下埋藏有塊狀及傾斜板狀良導體(圖9)。設高阻圍巖的電導率為0.01 S/m，塊狀體和板狀體的電導率為1 S/m；傾斜版狀體頂部埋深60 m，垂直深度250 m，塊狀體頂部埋深50 m，邊長200 m。采用中心回線裝置，飛行高度30 m，發(fā)射線框半徑15 m，線圈匝數(shù)為1，采用階躍波激發(fā)。全區(qū)設置11條測線，每條測線布設171個測點，共計1 881個測點，每個網(wǎng)格單元數(shù)約90 000個，正演過程共需進行9 350次矩陣分解，327 250次回代。若采用傳統(tǒng)串行計算時，總耗時約1 968.623 min(32.8 h)，采用并行加速開啟16個進程后總計算時間約194.336 min(3.24 h)，消耗內(nèi)存約5.2 G，效率提高了10.13倍,如表2所示。

表2 計算情況統(tǒng)計Table 2 Calculation statistics table

圖10是主剖面(x=0測線)的多測道圖，可以看出在早期表現(xiàn)出地形的響應，到了晚期，地形影響逐漸變小，表現(xiàn)出異常體的響應。圖11呈現(xiàn)出整個測區(qū)在四個時間點時的航空電磁響應結(jié)果，從圖中可以明顯看出在時間早期，僅存在地形的響應，由于地形凸起，電磁響應呈現(xiàn)出相對低異常，在x=0，y=-500附近，相對低異常達到極大值，由于塊狀異常體處于凸起地形下方，相對于傾斜板埋藏較深，故隨著時間延長，傾斜板的異常響應率先顯示出來，隨后出現(xiàn)塊狀體異常響應，且由于地下埋藏有異常體，導致地形的影響減弱，到晚期時，電磁波穿過異常體，僅剩地形的電磁響應。

圖9 復雜模型計算Fig.9 Calculation of complex models

圖10 復雜模型主剖面(x=0測線)多測道曲線Fig.10 Multiple trajectory map of complex model main section (x=0 line)

4 結(jié)論

利用MPI實現(xiàn)了三維正演模擬中的多測點并行計算和實現(xiàn)了獨立網(wǎng)格內(nèi)單元矩陣計算的OpenMP并行加速，最終實現(xiàn)了基于MPI+OpenMP并行加速策略對時間域航空電磁三維正演方法的并行化，并得出并行化后的正演程序，其加速比隨著開啟的進程 數(shù)增加呈非線性增長。

圖11 復雜地形響應曲面Fig.11 Complex terrain anomaly response surface map

由于“并行開銷”的存在，開啟的進程數(shù)越多，并行效率整體呈現(xiàn)出下降趨勢，即使只開啟一個進程，效率也比串行效率低。地形起伏模型的數(shù)值模擬體現(xiàn)了本文并行優(yōu)化后的正演方法的高效性，正確高效的三維正演為三維反演提供了可能。