陳 實(shí)， 吳文鸝， 顧觀文， 梁 萌

(中國地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所，廊坊 065000)

?

2.5維復(fù)電阻率反演并行算法設(shè)計

陳 實(shí)， 吳文鸝， 顧觀文， 梁 萌

(中國地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所，廊坊 065000)

近年來2.5維復(fù)電阻率反演取得了一些理論和應(yīng)用方面的成果，但計算過程復(fù)雜、計算量巨大，特別是正演計算過程中計算量大導(dǎo)致了反演速度較慢。這里在2.5維復(fù)電阻率反演計算基礎(chǔ)上，將反演算法中的正演過程的排列和頻率循環(huán)分割成多個規(guī)模較小的程序段，并發(fā)執(zhí)行正演過程中的每個排列和頻率循環(huán)的求解計算過程，并根據(jù)并行計算的節(jié)點(diǎn)數(shù)目或者計算核心數(shù)目，動態(tài)調(diào)整每個節(jié)點(diǎn)或者計算核心的計算量，從而提高了算法的計算效率和適應(yīng)性。

復(fù)電阻率; 并行計算; 2.5維反演

0 前言

在2.5 維復(fù)電阻率反演過程中，正演模擬計算過程比較復(fù)雜、數(shù)據(jù)處理與計算量大限制了該反演算法的實(shí)用性，因此在原有計算精度的基礎(chǔ)上提高計算速度，已經(jīng)成為該領(lǐng)域研究的重點(diǎn)。提高程序運(yùn)算速度主要采取兩種方法：①從算法入手改進(jìn)算法的設(shè)計,提高算法的計算效率(如采用更有效求解復(fù)線性方程組算法等)；②從計算設(shè)備入手，采用多核或者計算機(jī)集群并行計算來實(shí)現(xiàn)為程序加速。這里采用多核計算機(jī)，在保持2.5 維復(fù)電阻率反演算法主體結(jié)構(gòu)的基礎(chǔ)上，將反演過程中的正演過程按照排列和頻率循環(huán)分割形成新的計算線程，并將新的計算線程分解到不同的計算核心上，以實(shí)現(xiàn)程序的并行計算，減少了算法的運(yùn)行時間，提高了算法的運(yùn)算速度，程序結(jié)構(gòu)相對簡單，可以在較少修改的情況下部署到其他計算平臺上。

1 2.5維復(fù)電阻率反演原理與并行計算分析

趙廣茂[1]對2.5 維復(fù)電阻率反演理論過程做了初步推導(dǎo)。范翠松[2]對該反演過程進(jìn)行了進(jìn)一步說明。

為提高反演穩(wěn)定性，對參數(shù)進(jìn)行歸一化處理即令

x=(ln(ρ0))，ln(m)，ln(τ),ln(c))，并以此構(gòu)建單一排列的反演目標(biāo)函數(shù)：

φ=‖Uaa-UaA(x)‖+ ‖Uφφ-Uφφ(x)‖+μλ‖R‖

(1)

式中：a和φ分別為觀測電場的振幅向量與相位向量；A(x)和φ(x)分別為電場振幅與相位的正演函數(shù)；Ua和Uφ為歸一化對角陣；λ是阻尼因子；μ為縮放系數(shù)；R為模型的二階光滑度算子矩陣。 將式(1)在xk鄰域線性展開并極小化得到實(shí)系數(shù)反演方程組：

(2)

其中：ΔdA和Δdφ分別為電場振幅和相位的絕對擬合差向量；LA和Lφ分別為電場振幅與相位的靈敏度矩陣；Δxk為模型修正向量,通過求解方程(2)即可得到xk+1，進(jìn)入下次迭代過程。當(dāng)利用N個排列進(jìn)行反演時，由于公式(2)為線性過程，通過求解左端靈敏度相關(guān)項(xiàng)，并作線性累加求和，并將右端項(xiàng)線性累加求和整理得到新方程組(3)。

(P+μλRTR)·Δxk=B

(3)

其中：P為靈敏度矩陣求和；B為右端項(xiàng)擬合差向量求和。通過式(3)得到線性方程組并求解該線性方程組即可實(shí)現(xiàn)多個排列的繁衍計算。根據(jù)上述理論分析可以得到2.5維復(fù)電阻率，其反演計算流程如圖1所示。

1)讀取程序使用的數(shù)據(jù)和模型參數(shù)。

2)計算反演光滑度算子矩陣。

3)調(diào)用正演計算模塊進(jìn)行正演計算并生成形成反演矩陣。

4)調(diào)用反演模塊求解反演線性方程組[5]，并更新模型參數(shù)。

5)調(diào)用正演計算模塊進(jìn)行計算誤差。

6)判斷誤差是否收斂。誤差不收斂則返回步驟4)，誤差收斂繼續(xù)判斷誤差是否小于停止運(yùn)算指定誤差，誤差大于停止運(yùn)算指定誤差則返回步驟3)進(jìn)行下一次迭代。誤差小于或等于停止運(yùn)算指定誤差則停止運(yùn)算輸出反演結(jié)果。

實(shí)驗(yàn)中運(yùn)用的算例觀測方式為偶極-偶極，模型共進(jìn)行了6個排列的數(shù)據(jù)采集，每個排列10個觀測點(diǎn)，發(fā)射偶極長度為100 m，接收偶極長度為100 m。數(shù)據(jù)采集的工作頻率分別為0.02 Hz、0.1 Hz、1 Hz、8 Hz、32 Hz、128 Hz。模型計算采用有限元法，地電模型剖分網(wǎng)格見圖2，水平方向剖分為55塊，深度方向地面以上剖分為10層，地面以下剖分為18層。

根據(jù)上述剖分過程和計算流程，對原有串行算法各個中間過程進(jìn)行時間統(tǒng)計，其中讀取數(shù)據(jù)和模型參數(shù)用時為0.52 s，計算反演光滑度矩陣用時為42.2 s。迭代一次計算過程中，正演計算模型變化量并形成反演矩陣的過程計算時間為7 537.13 s，求解反演矩陣并更新模型過程為57.32 s，正演誤差計算過程共用時間1 758.96 s。

表1 計算模塊時間統(tǒng)計表

根據(jù)上述時間統(tǒng)計圖正演計算過程占整個迭代過程計算時間的79%，反演過程中正演計算合成過程占整個迭代計算時間的18.4%，反演求解線性方程組過程只占迭代計算時間的1%左右。如果能將正演計算和反演中正演合成過程實(shí)現(xiàn)并行計算就可以大幅度提高整個迭代計算過程的運(yùn)算速度。

2 正演過程并行算法實(shí)現(xiàn)

2.1 并行算法計算量分配方案

正演過程每個排列和頻率的計算過程相對獨(dú)立并且結(jié)構(gòu)相似，適合進(jìn)行并行計算，在并行計算結(jié)束時，將每一個計算單元計算出的地表場合成反演線性方程組的左端項(xiàng)和右端項(xiàng)。反演算法正演過程并行分解粒度可以用排列級、頻率級[6]、波數(shù)級和觀測點(diǎn)級。根據(jù)算法實(shí)際實(shí)現(xiàn)過程中所選擇的實(shí)驗(yàn)平臺，最終確定并行分解粒度為規(guī)模適中的頻率級，即正演計算過程中有m個排列、每個排列有n個頻率，每個頻率存在t個計算波數(shù)，實(shí)驗(yàn)平臺存在u個計算核心，將程序分解為m*n個計算線程，其中每個計算核心分配m*n/u個計算線程[7]，其流程描述見圖2。

圖1 2.5 維復(fù)電阻率反演計算流程圖

圖2 并行計算過程分解

2.2 并行過程程序設(shè)計

由于整體程序采用Fortran實(shí)現(xiàn)，所以上述算法中涉及的并行過程可以采用支持Fortran的并行開發(fā)包MPI或者OpenMP實(shí)現(xiàn)[8]。其中MPI適合用于多進(jìn)程并行計算[9]，主要用于將計算程序分解到網(wǎng)絡(luò)中的多個計算節(jié)點(diǎn)上，實(shí)現(xiàn)并行計算或者將程序分解成多個計算進(jìn)程，每個計算進(jìn)程雖然工作在同一個計算機(jī)上，但是每個計算進(jìn)程有獨(dú)立的數(shù)據(jù)與存儲空間，計算進(jìn)程相互獨(dú)立，無法實(shí)現(xiàn)數(shù)據(jù)直接訪問，計算進(jìn)程間數(shù)據(jù)傳輸依靠MPI_Send和MPI_Recieve函數(shù)通過網(wǎng)絡(luò)數(shù)據(jù)傳輸實(shí)現(xiàn)[10]。OpenMP主要用于實(shí)現(xiàn)單進(jìn)程多線程并行運(yùn)算，其并行過程采用線程實(shí)現(xiàn)，多個并行過程數(shù)據(jù)共享采用欄柵或者臨界區(qū)實(shí)現(xiàn)，比較適合程序結(jié)構(gòu)較復(fù)雜，中間共享數(shù)據(jù)結(jié)構(gòu)較多的程序?qū)崿F(xiàn)并行計算[11]。針對2.5 維復(fù)電阻率反演算法程序，采用OpenMP實(shí)現(xiàn)并行計算。

3 反演并行算法算例分析

實(shí)驗(yàn)所選工作計算機(jī)為DELL雙核12線程核工作站，工作頻率為3.46 GHz，內(nèi)存為24 G。算例模型為6個排列，每個排列10個觀測點(diǎn)，6個工作頻率，地電模型水平方向剖分為55塊，深度方向剖分為28層。異常體位置見圖3，其中異常體的中心為(600,-200),長度和寬度分別為200 m。

圖3 異常體位置

反演每次迭代運(yùn)行時間對比見表2。

并行反演迭代20次運(yùn)行結(jié)果如圖4，其中異常體中心在(600,-200),結(jié)果與算例模型相符，與串行程序結(jié)果相同。

保持該算例10個觀測點(diǎn)、6個工作頻率，分別取1個排列、2個排列、3個排列、4個排列、5個排列和6個排列重新構(gòu)造算例，進(jìn)行串行和并行運(yùn)算時間對比，對比結(jié)果見圖5。

由圖5可以看出，串行計算過程計算時間隨計算量增加而線性增加，并行計算過程每個排列被按照頻率數(shù)分解為6個計算線程，計算核心有12個。排列數(shù)為1時計算線程共6個，分配在6個計算核心上，其余6個計算核心閑置，并行運(yùn)算時間為串行運(yùn)算時間的20.7 %；排列數(shù)為2時，計算線程為12個，分配在12個計算線程上，12個計算核心全部參加運(yùn)算，運(yùn)算時間比1個排列稍有增加，并行運(yùn)算時間為串行運(yùn)算時間的11.3 %；排列數(shù)大于等于3時，前12個計算線程分別在12個計算核心上進(jìn)行計算，待某個計算核心計算完成后，在剩余計算線程選擇一個調(diào)度到空閑的計算核心進(jìn)行計算，其中排列數(shù)為3時，并行運(yùn)算時間為串行運(yùn)算時間的13.0 %，排列數(shù)為4時，并行運(yùn)算時間為串行運(yùn)算時間的10.2 %，排列數(shù)為5時，并行運(yùn)算時間為串行運(yùn)算時間的12.4 %，排列數(shù)為6時，并行運(yùn)算時間為串行運(yùn)算時間的10.8 %。

由此可以看出，并行計算時間也隨計算量緩慢增加，其斜率與總計算線程數(shù)和計算核心數(shù)目的比值成正比，假設(shè)為每個計算核心的計算時間Ti，并行計算時間曲線的斜率和Ti的最大值MAX(Ti)成正比。

表2 每次迭代并行和串行時間對比

圖4 反演斷面圖

圖5 不同排列串行與并行時間對比

圖6 并行計算時間和計算核心數(shù)關(guān)系圖

保持該算例6個排列、10個觀測點(diǎn)、6個工作頻率，將程序移植HP9000 Superdome服務(wù)器上，該機(jī)型體系結(jié)構(gòu)為IA-64,CPU工作頻率為1 598 MHz，計算核心分別取1～64個，其運(yùn)行時間見圖6。

并行計算時間隨著計算核心數(shù)增加而迅速下降，由于算例中計算量被分為36組共36個計算線程，當(dāng)計算核心增加到36以上時，前36個計算核心參與計算，后面的計算核心并未參與工作，計算時間與36個計算核心時基本相同。

算例實(shí)驗(yàn)使用了兩種運(yùn)行環(huán)境，其中Dell T7500工作站(2核6線程)具有12個計算核心全部參與運(yùn)算，CPU工作頻率為3 460 MHz，HP9000 Superdome(64核)具有64個運(yùn)算核心，CPU工作頻率為1 598 MHz，其中參與運(yùn)算的工作核心數(shù)為36，其運(yùn)行時間對比見表3。其中Dell T7500工作站為當(dāng)前最新的計算工作站，其指令結(jié)構(gòu)和編譯器優(yōu)化程度較高，所以其整體加速效果優(yōu)于HP9000 Superdome服務(wù)器。

表3 不同平臺運(yùn)行時間對比

4 結(jié)論

1)反演并行計算程序在保持原有精度的基礎(chǔ)上，使原有串行程序能夠在多核計算機(jī)上進(jìn)行并行運(yùn)算，充分利用計算平臺的計算資源，提高了程序的運(yùn)算速度。

2)在劃分并行任務(wù)的顆粒度時，在充分使用計算核心的前提下，每次分配給計算核心的任務(wù)量盡量大。因?yàn)樵诓⑿腥蝿?wù)開始時需要對每個計算核心分配資源，每次分配資源耗時相對較多，例如開辟多維數(shù)組，如果計算任務(wù)相對較小，則本次分配資源占用的比例相對較高，程序的運(yùn)行效率就會比較低。

[1] 趙廣茂.帶地形的復(fù)電阻率2.5維電磁場正反演研究[D].長春：吉林大學(xué),2009. ZHAO G M. Research of complex resistivity 2.5D electromagnetic forword and inversion with topography[D]. Changchun:Jilin University,2009.(In Chinese)

[2] 范翠松,李桐林,嚴(yán)加永.2.5維復(fù)電阻率反演及其應(yīng)用試驗(yàn)[J].地球物理學(xué)報,2012,55(12):4045-4046. FAN C S,LI T L,YAN J Y. Research and application experiment on 2.5D SIP inversion[J].Chinese Journal of Geophysics,2012,55(12):4045-4046.(In Chinese)

[3] 王家映.地球物理反演理論[M].北京:高等教育出版社,2002. WANG J Y. Inversion theory of physical geography[M].Beijing:Higher education press,2002.(In Chinese)

[4] 李金銘.地電場與電法勘探[M].北京:地質(zhì)出版社,2005. LI J M. Geoelectricity and electrical prospecting[M].Beijing: Geological Press,2005.(In Chinese)

[5] 安學(xué)慶.求解大型稀疏線性方程組算法研究[D].鄭州:鄭州大學(xué),2006. AN X Q. The solution of large-scale sparse linear system of equations algorithm research[D]. Zhengzhou: Zhengzhou University,2006.(In Chinese)

[6] 李焱,胡祥云,吳桂桔,等.基于MPI的二維大地電磁正演的并行計算[J].地震地質(zhì),2010,32(3): 392-401. LI Y,HU X Y, WU G J, et al. Parallel computation of 2-D magnetotelluric forward modeling based on MPI[J].Seismology and Geology,2010,32(3):392-401.(In Chinese)

[7] 顧觀文, 梁萌, 吳文鸝.基于并行處理的CSAMT擬二維反演解釋[J].物探與化探.2010，34(3):400-402. GU G W,LIANG M ,WU W L.The quasi two dimensional inversion and interpretation of CSAMT based on parallel processing.Geophysical and Geochemical Exploration,2010,34(3):400-402.(In Chinese)

[8] Barry Wilkinson Michael Allen.并行程序設(shè)計 [M] . 陸鑫達(dá),譯.北京:機(jī)械工業(yè)出版社, 2005. BARRY WILKINSON MICHAEL ALLEN. Design of parallel programing[M].Lu Xinda, Beijing: Mechanical Industry Press, 2005.(In Chinese)

[9] 蔣英,雷永梅.基于MPI的幾種算法的并行編程通用算法[J].計算機(jī)工程與應(yīng)用.2003,(03):140-141. JIANG Y,LEI Y M. Parallel programming universal algorithm of several algorithms based on MPI[J]. Computer Engineering and Applications.2003,(03):140-141.(In Chinese)

[10]曾志峰.Linux環(huán)境下MPI并行編程與算法實(shí)現(xiàn)研究[J].航空計算技術(shù).2004,32(2):62-64. ZENG ZH F. Under the environment of Linux MPI parallel programming and algorithm implementation research[J].Aeronautical Computer Technique.2004,32(2):62-64.(In Chinese)

[11]潘小敏,皮維超,盛新慶.基于共享內(nèi)存的高效OpenMP并行多層快速多極子算法[J].北京理工大學(xué)學(xué)報.2012.2(32):165-166. PAN X M,PI W CH,SHENG X Q. Efficient parellelization of multilevel fast multipole algrithm based on OpenMP[J].Transactions of Beijing Institute of Technology.2012.2(32):165-166.(In Chinese)

Design of 2.5 dimensional complex resistivity inversion parallel computing scheme

CHEN Shi,WU Wen-li,GU Guan-wen,LIANG Meng

(Institute of Geophysical and Geochemical Exploration, CAGS, Langfang 065000, China)

Some theoretical and application aspects of results have been made for 2.5 dimensional complex resistivity inversion, but complex mathematical calcultionsin forward processes slow down the inversion speed. Based on the 2.5 dimensional complex resistivity inversion，the forward processes of the inversion algorithm will be divided into much smaller forward cycles accordingtoits arrangement and frequency.the amount of calculations will be assigned dynamically to each node or the calculation core depending on the number of nodes or the calculation of the number of parallel computing core, therefore the adaptability and efficiency of the algorithm will be greatly improved in this way.

complex; resistivity; parallel algorithmparallel algorithm; 2.5 dimensional inversion

2014-12-29 改回日期：2015-04-24

國家863項(xiàng)目(2014AA06A610)；物化探所所長基金項(xiàng)目(AS2013J07)

陳實(shí)(1982-)，男，助理工程師，現(xiàn)主要從事軟件工程以及并行計算與三維可視化軟件開發(fā)工作，E-mail:chen884488shi@126.com。

1001-1749(2015)04-0416-06

P 631.3

A

10.3969/j.issn.1001-1749.2015.04.02