2.5維復電阻率反演并行算法設計

陳 實， 吳文鸝， 顧觀文， 梁 萌

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

?

2.5維復電阻率反演并行算法設計

陳 實， 吳文鸝， 顧觀文， 梁 萌

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

近年來2.5維復電阻率反演取得了一些理論和應用方面的成果，但計算過程復雜、計算量巨大，特別是正演計算過程中計算量大導致了反演速度較慢。這里在2.5維復電阻率反演計算基礎上，將反演算法中的正演過程的排列和頻率循環分割成多個規模較小的程序段，并發執行正演過程中的每個排列和頻率循環的求解計算過程，并根據并行計算的節點數目或者計算核心數目，動態調整每個節點或者計算核心的計算量，從而提高了算法的計算效率和適應性。

復電阻率; 并行計算; 2.5維反演

0 前言

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

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

趙廣茂[1]對2.5 維復電阻率反演理論過程做了初步推導。范翠松[2]對該反演過程進行了進一步說明。

為提高反演穩定性，對參數進行歸一化處理即令

x=(ln(ρ0))，ln(m)，ln(τ),ln(c))，并以此構建單一排列的反演目標函數：

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

(1)

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

(2)

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

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

(3)

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

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

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

3)調用正演計算模塊進行正演計算并生成形成反演矩陣。

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

5)調用正演計算模塊進行計算誤差。

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

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

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

表1 計算模塊時間統計表

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

2 正演過程并行算法實現

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

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

圖1 2.5 維復電阻率反演計算流程圖

圖2 并行計算過程分解

2.2 并行過程程序設計

由于整體程序采用Fortran實現，所以上述算法中涉及的并行過程可以采用支持Fortran的并行開發包MPI或者OpenMP實現[8]。其中MPI適合用于多進程并行計算[9]，主要用于將計算程序分解到網絡中的多個計算節點上，實現并行計算或者將程序分解成多個計算進程，每個計算進程雖然工作在同一個計算機上，但是每個計算進程有獨立的數據與存儲空間，計算進程相互獨立，無法實現數據直接訪問，計算進程間數據傳輸依靠MPI_Send和MPI_Recieve函數通過網絡數據傳輸實現[10]。OpenMP主要用于實現單進程多線程并行運算，其并行過程采用線程實現，多個并行過程數據共享采用欄柵或者臨界區實現，比較適合程序結構較復雜，中間共享數據結構較多的程序實現并行計算[11]。針對2.5 維復電阻率反演算法程序，采用OpenMP實現并行計算。

3 反演并行算法算例分析

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

圖3 異常體位置

反演每次迭代運行時間對比見表2。

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

保持該算例10個觀測點、6個工作頻率，分別取1個排列、2個排列、3個排列、4個排列、5個排列和6個排列重新構造算例，進行串行和并行運算時間對比，對比結果見圖5。

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

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

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

圖4 反演斷面圖

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

圖6 并行計算時間和計算核心數關系圖

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

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

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

表3 不同平臺運行時間對比

4 結論

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

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

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

[2] 范翠松,李桐林,嚴加永.2.5維復電阻率反演及其應用試驗[J].地球物理學報,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].北京:地質出版社,2005. LI J M. Geoelectricity and electrical prospecting[M].Beijing: Geological Press,2005.(In Chinese)

[5] 安學慶.求解大型稀疏線性方程組算法研究[D].鄭州:鄭州大學,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].地震地質,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.并行程序設計 [M] . 陸鑫達,譯.北京:機械工業出版社, 2005. BARRY WILKINSON MICHAEL ALLEN. Design of parallel programing[M].Lu Xinda, Beijing: Mechanical Industry Press, 2005.(In Chinese)

[9] 蔣英,雷永梅.基于MPI的幾種算法的并行編程通用算法[J].計算機工程與應用.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環境下MPI并行編程與算法實現研究[J].航空計算技術.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]潘小敏,皮維超,盛新慶.基于共享內存的高效OpenMP并行多層快速多極子算法[J].北京理工大學學報.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項目(2014AA06A610)；物化探所所長基金項目(AS2013J07)

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

1001-1749(2015)04-0416-06

P 631.3

A

10.3969/j.issn.1001-1749.2015.04.02