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

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

?

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

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

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

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

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

0 前言

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

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

趙廣茂[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個(gè)排列進(jìn)行反演時(shí)，由于公式(2)為線性過程，通過求解左端靈敏度相關(guān)項(xiàng)，并作線性累加求和，并將右端項(xiàng)線性累加求和整理得到新方程組(3)。

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

(3)

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

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

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

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

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

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

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

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

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

表1 計(jì)算模塊時(shí)間統(tǒng)計(jì)表

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

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

2.1 并行算法計(jì)算量分配方案

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

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

圖2 并行計(jì)算過程分解

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

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

3 反演并行算法算例分析

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

圖3 異常體位置

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

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

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

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

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

表2 每次迭代并行和串行時(shí)間對比

圖4 反演斷面圖

圖5 不同排列串行與并行時(shí)間對比

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

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

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

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

表3 不同平臺(tái)運(yùn)行時(shí)間對比

4 結(jié)論

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

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

[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é)報(bào),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ì)算[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è)計(jì) [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ì)算機(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].航空計(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é)報(bào).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)主要從事軟件工程以及并行計(jì)算與三維可視化軟件開發(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