復(fù)電阻率法二維正演并行計算研究

張志勇 , 譚捍東, 王春陽， 汪 茂， 何中波， 梁 鵬

(1.核工業(yè)北京地質(zhì)研究院中核集團鈾資源勘查與評價技術(shù)重點實驗室，北京 100029；2.中國地質(zhì)大學(xué)(北京) 地球物理與信息技術(shù)學(xué)院，北京 100083;3.中水東北勘測設(shè)計研究有限責(zé)任公司，長春 130021)

0 引言

復(fù)電阻率法(CR)又稱為譜激發(fā)極化法(SIP)是以巖礦石電阻率的頻譜特性差異為基礎(chǔ)，通過觀測大地的視復(fù)電阻率頻譜或激發(fā)極化場的衰減曲線，來尋找極化目標體的一種電法勘探方法[1]。相對其他電法勘探分支方法，復(fù)電阻率法有以下優(yōu)點：①該方法能得到多個電性參數(shù)的信息，結(jié)合這些參數(shù)的對比解釋可提供更豐富的地電信息；②采用輕便的采集裝備，使其在復(fù)雜地形地區(qū)開展勘探工作更有利。目前，復(fù)電阻率法已成功應(yīng)用在固體礦產(chǎn)勘探[2-4]、水文地質(zhì)[5]、油氣資源[6]、監(jiān)測環(huán)境[7]、工程地質(zhì)[8]等方面。

目前計算機技術(shù)發(fā)展迅速，將MPI等并行計算技術(shù)引入到地球物理數(shù)值模擬中是當(dāng)前研究熱點。當(dāng)網(wǎng)格剖分單元較多，多發(fā)射源、多頻率的情況下，完成復(fù)電阻率法二維有限元數(shù)值模擬需要的計算時間較長，因此有必要開發(fā)復(fù)電阻率法二維正演并行算法來提高其計算效率。

1 MPI簡介

MPI(Message Passing Interface)是一個“庫”而非一門語言，其具有C++和FORTRAN90的函數(shù)調(diào)用接口[9]。筆者采用Fortran和MPICH2相結(jié)合的開發(fā)工具編寫并行程序。

MPI程序設(shè)計包含兩種并行模式：①主從模式；②對等模式[10](圖1)。

1)主從模式，也就是在所有進程中定義其中的一個作為主進程，除主進程外的所有進程都當(dāng)作子進程。通過主進程可以進行廣播數(shù)據(jù)或發(fā)送數(shù)據(jù)等一系列操作，同時也可以將分配好的任務(wù)分發(fā)給子進程。等到所有的子進程都完成計算以后，主進程通過收集操作從每個子進程得到任務(wù)結(jié)果，然后再看是否需要下一步的操作。

2)對等模式，就是所實現(xiàn)的MPI程序中各個進程間的計算量和功能大體一致。

圖1 并行模式分類Fig.1 Parallel mode classification(a)對等模式；(b)主從模式

2 復(fù)電阻率法二維正演

2.1 二維數(shù)值模擬

二維地電模型如圖2所示，其中構(gòu)造走向方向為y軸，假設(shè)其電性參數(shù)僅僅在x-z平面是變化的而沿y軸方向不變。

圖2 二維地電模型示意圖Fig.2 2D geoelectric model

假設(shè)時間因子是eiωt，忽略位移電流的影響，二次場滿足的麥克斯韋方程組為[11]

(1)

傅里葉變換

(2)

沿走向方向?qū)κ?1)作傅里葉變換，整理后得到波數(shù)域中滿足的麥克斯韋方程組

(3)

(4)

利用四節(jié)點的等參單元法進行數(shù)值模擬，將式(3)和式(4)分別運用伽里金法、格林公式，且正演使用第一類邊界條件(邊界網(wǎng)格上的電磁場值為零)，可以得到本文二維數(shù)值模擬滿足的有限單元法計算公式[12-16]：

(6)

2.2 Cole-Cole模型

為同時考慮激電效應(yīng)和電磁效應(yīng)，將大地的實際電阻率由Cole-Cole模型定義得到的復(fù)電阻率替換。Pelton等[17]通過對大量巖石、礦石標本的一些研究，總結(jié)出就均勻巖石、礦石而言，它們的復(fù)電阻率頻譜可以用Cole-Cole模型來表示。

(7)

其中:ρ(iω)是復(fù)電阻率；m是極化率；ρ0是零頻電阻率；τ是時間常數(shù)；c是頻率相關(guān)系數(shù)。

2.3 有限單元計算公式

引入Cole-Cole后，通過單元分析和系數(shù)矩陣存儲，最終正演計算形成的矩陣方程可寫成如下形式

(8)

k1(i,j)=

(9)

(10)

(11)

k4(i,j)=

(12)

(13)

(14)

2.4 矩陣的壓縮存儲

正演矩陣方程(8)中；左端項系數(shù)矩陣K是稀疏且對稱的，全部存儲K矩陣需要的內(nèi)存較大。由于矩陣K每行的非零元素個數(shù)有限，充分利用這個特點，可以有效地節(jié)約內(nèi)存空間。結(jié)合有限單元法的特點，單元分析時每個節(jié)點僅與其周圍的幾個節(jié)點參與運算。以筆者采用的4節(jié)點等參單元為例，最終形成的矩陣中每行最多只有18個非零元素，在網(wǎng)格剖分后，可根據(jù)有限元單元分析特點得到非零元素的行列號。采用CSR(按行壓縮存儲)方法，通過三個一維向量就可完成矩陣的存儲，便于使用共軛梯度法來求解方程組。

2.5 求解方程組

正演形成的方程組采用林紹忠[18]介紹的SSOR—PCG(對稱逐步超松弛預(yù)處理共軛梯度法)方法解方程。

解形如Ax=b的方程時，SSOR法的預(yù)處理矩陣M為：

M=(2-ω)-1(D/ω+L)(D/ω)-1(D/ω+L)T

(15)

其中:D為A的對角陣;L為A的嚴格下三角矩陣;ω為松弛因子，本文中將松弛因子取為“1”，初始x為零向量。那么迭代格式可以簡化為

(16)

2.6 正演正確性驗證

圖3 二維棱柱體模型Fig.3 Two-dimensional prism model

3 二維正演并行算法總體設(shè)計

在對開發(fā)的復(fù)電阻率法二維正演串行程序深入分析后，可以發(fā)現(xiàn)正演中，其大部分計算時間都花費在頻率循環(huán)部分求解電磁場值上面，每個發(fā)射源對應(yīng)頻率的方程組的求解等都是相互獨立的且互不影響的，并行性較好，因此我們采用主從模式，基于多發(fā)射源和多頻率來實現(xiàn)正演并行計算，并對其并行效率進行分析。并行的基本思路是：①將MPI環(huán)境初始化MPI_INIT( )；②主進程讀入所需的參數(shù)文件；③主進程廣播數(shù)據(jù)到各子進程MPI_Bcast( )；④各進程分別計算自己的任務(wù)；⑤主進程收集和整理所有進程的結(jié)果MPI_GATHERV()；⑥結(jié)束MPI并行環(huán)境MPI_FINALIZE( )[19-21]。

圖4 二維棱柱體模型計算得到Ex實虛部對比曲線Fig.4 The Ex response comparison between the 2D finite element results and the results from 1DCSEM(a)f=8.0 Hz; (b)f=8.0 Hz; (c)f=16.0 Hz; (d)f=16.0 Hz

3.1 模型

如圖5所示，在背景電阻率為100 Ω·m的地下介質(zhì)中存在一個長200 m厚100 m的異常體，其頂部埋深為160 m，其中ρ0=10 Ω·m，m=0.5，c=0.5，τ=30 s。采用的裝置為偶極-偶極，沿x方向水平放置電偶極子，xz方向網(wǎng)格剖分為121×64。以下正演并行程序均在18個節(jié)點的集群上運算，其計算環(huán)境均為：Linux操作系統(tǒng)，CPU-Intel(R) Xeon(R) CPU E5-2620 v2 @ 2.1GHz，內(nèi)存128G，開發(fā)語言為Fortran，編譯器為ifort，并行環(huán)境為MPICH2。

圖5 低阻高級化模型Fig.5 Low resistivity and high polarization model

各個發(fā)射源對應(yīng)頻率的計算是獨立的，基于發(fā)射源對應(yīng)的頻率來進行并行計算的思路，其流程如圖6所示。

3.1.1 效率分析

采用圖5所示的模型，21個發(fā)射源(其位置從X= -1 000 m到X= 1 000 m之間每隔100 m一個，共21個源)，從X= -1 000 m到X= 1 000 m之間每隔50 m布置一個測點，每個排列除發(fā)射源點外共40個接收點。5個頻率(f= 0.1 Hz, 1.0 Hz, 8.0 Hz, 32.0 Hz, 128.0 Hz)，用不同節(jié)點數(shù)來進行并行計算，將其計算時間和串行程序計算時間統(tǒng)計如表1所示。

圖6 正演并行流程圖Fig.6 Forward parallel flow chart

程序運行模式節(jié)點數(shù)量每個節(jié)點參與計算的進程數(shù)各節(jié)點對應(yīng)進程計算的源對應(yīng)的頻率數(shù)量網(wǎng)格大小時間/s并行加速比并行效率/%串行程序11105121×642207.39無無并行程序2153,52121×641147.471.92496.23135,35,35121×64745.762.96098.74127,26,26,26121×64635.783.47286.85121,21,21,21,21121×64501.224.40488.16118,18,18,17,17,17121×64452.274.88181.37115,15,15,15,15,15,15121×64377.525.84783.5

從表1可以看出，在每個節(jié)點只有一個進程參與計算的前提下，當(dāng)參與并行計算的節(jié)點數(shù)增加時，加速比增大。當(dāng)有7個節(jié)點參與并行計算時，其計算速度約為串行程序計算速度的5.847倍，很大程度地降低了復(fù)電祖率法二維正演所需時間。但是，從表2中并行效率這一列看出，并不是參與并行的節(jié)點數(shù)越多，并行效率越高。當(dāng)參與并行的節(jié)點數(shù)從2增加到3時，并行效率增加，而當(dāng)參與并行的節(jié)點數(shù)從3增加到4(或節(jié)點數(shù)從5增加到6)時，并行效率反而降低，這是由于節(jié)點等待或節(jié)點間的數(shù)據(jù)交換等通信耗費一些時間。

3.1.2 正確性驗證

采用多發(fā)射源和多頻率進行正演計算時，筆者僅展示其中一部分數(shù)據(jù)的對比結(jié)果。如圖7所示，圖7中Tx、Rx分別代表發(fā)射源和接收點的位置，串行計算結(jié)果用紅色圓圈表示，并行計算結(jié)果用藍色點表示，計算結(jié)果完全吻合。

4 結(jié)論

當(dāng)發(fā)射源、頻率較多時，完成復(fù)電阻率法二維正演所需時間較長，我們通過引入MPI并行計算，實現(xiàn)了復(fù)電阻率法二維正演并行算法，極大地提高了正演的計算效率。通過理論模型的試算結(jié)果表明，該并行算法是可靠的，穩(wěn)定的，較好地解決了復(fù)電阻率法二維數(shù)值模擬的計算速度問題，為其反演并行計算打下良好的基礎(chǔ)。對其并行計算時間分析可知，當(dāng)每個節(jié)點的進程數(shù)都為“1”時，隨著節(jié)點數(shù)的增多，其并行加速比逐漸增大，但其并行效率并不與參與計算的節(jié)點數(shù)成正比。

圖7 正演并行結(jié)果與串行結(jié)果對比圖Fig.7 The comparison between forward modeling parallel and serial results(a)Tx=-1 000 m Rx=-300 m; (b)Tx=-1 000 m Rx=-300 m;(c)Tx=-1 000 m Rx=0 m; (d)Tx=-1 000 m Rx=0 m;(e)Tx=-1 000 m Rx=300 m; (f)Tx=-1 000 m Rx=300 m;(g)Tx=-500 m Rx=100 m; (h)Tx=-500 m Rx=100 m

參考文獻：

[1] 楊進. 環(huán)境與工程地球物理[M]. 北京: 地質(zhì)出版社, 2011.

YANG J. Environmental and engineering geophysics[M]. Beijing: Geological Publishing House, 2011. (In Chinese)

[2] 楊振威, 嚴加永, 陳向斌. 頻譜激電法在安徽沙溪斑巖銅礦中的應(yīng)用[J]. 地球物理學(xué)進展, 2013, 28(4): 2014-2023.

YANG Z W, YAN J Y, CHEN X B. The application of spectral induced polarization in Shaxi porphyry copper in Anhui province[J]. Progress in Geophysics, 2013, 28(4): 2014-2023. (In Chinese)

[3] 曹蔚杰, 單明良, 高勇, 等. 復(fù)電阻率法(CR)在銅鉬礦勘查中的應(yīng)用效果[J]. 工程地球物理學(xué)報, 2014, 11(2): 203-207.

CAO W J, SHAN M L, GAO Y ，et al. The application of complex resistivity method(CR) to copper-molybdenum Mine exploration[J]. Chinese Journal of Engineering Geophysics, 2014, 11(2): 203-207. (In Chinese)

[4] 徐自生, 張麗, 唐偉, 等. 復(fù)電阻率法(CR)在內(nèi)蒙古那吉河鉛鋅礦探測中的應(yīng)用[J]. 礦產(chǎn)勘查, 2015(2):165-170.

XU Z S, ZHANG L, TANG W, et al. Application of complex resistivity method (CR) in prospecting lead-zinc deposit in Inner Mongolia[J]. Mineral Exploration, 2015(2): 165-170. (In Chinese)

[5] REVIl A, KARAOULIS M, JOHNSON T, et al. Review: some low-frequency electrical methods for subsurface characterization and monitor in hydrogeology[J]. Hydrogeology Journal, 2012, 20(4): 617-658.

[6] 許傳建, 徐自生, 楊志成, 等. 復(fù)電阻率(CR)法探測油氣藏的應(yīng)用效果[J]. 石油地球物理勘探, 2004, A: 31-35.

XU C J, XU Z S, YANG Z C, et al. The application of complex resistivity method(CR) to oil gas exploration[J]. Oil Geophysical Prospecting, 2004, A: 31-35. (In Chinese)

[7]FLORES OROZCO A, KEMNA A, OBERD?RSTER C, et al. Delineation of subsurface hydrocarbon contamination at a former hydrogenation plant using spectral induced polarization imaging[J]. Journal of Contaminant Hydrology, 2012, 136-137: 131-144.

[8] BREEDE K, KEMNA A, ESSER O, et al. Spectral induced polarization measurements on variably saturated sand-clay mixtures[J]. Near surface geophysics, 2012, 10(6):479-489.

[9] 張武生, 薛巍, 李建江, 等. MPI并行程序設(shè)計實例教程[M]. 北京: 清華大學(xué)出版社, 2009.

ZHANG W S, XUE W, LI J J, et al. Parallel programming with MPI[M]. Beijing: Tsinghua University Press, 2009. (In Chinese)

[10] 都志輝, 李三立, 陳渝, 等. 高性能計算之并行編程技術(shù)——MPI并行程序設(shè)計[M]. 北京: 清華大學(xué)出版社, 2001.

DU Z H, LI S L, CHEN Y, et al. High performance parallel computing technology——MPI parallel programming[M]. Beijing: Tsinghua University Press, 2001. (In Chinese)

[11] NABIGHIAN M. N. Electromagnetic Methods in applied geophysics[M].Theory. Society of Exploration Geophysicists, 1988.

[12] 徐世浙. 地球物理中的有限單元法[M]. 北京: 科學(xué)出版社, 1994.

XU S Z. Finite element method in geophysics[M]. Beijing: Science Press, 1994. (In Chinese)

[13] 張斌. 可控源音頻大地電磁二維正反演研究[D]. 北京: 中國地質(zhì)大學(xué), 2012.

ZHANG B. Research of 2D CSAMT forward and inversion[D]. Beijing: China University of Geosciences(Beijing), 2012. (In Chinese)

[14] 張志勇. 復(fù)電阻率法二維有限元正演研究[D]. 北京: 中國地質(zhì)大學(xué), 2013.

ZHANG Z Y. 2D FEM forward modeling research of complex resistivity method[D]. Beijing: China University of Geosciences(Beijing), 2013. (In Chinese)

[15] 王寒冰. 極低頻電磁法頻率域二維正反演研究[D]. 北京: 中國地質(zhì)大學(xué), 2014.

WANG H B. Research of 2D frequency-domain ELF forward and inversion[D]. Beijing: China University of Geosciences, 2014. (In Chinese)

[16] 桂兵. 可控源音頻大地電磁法張量數(shù)據(jù)二維反演研究[D]. 北京: 中國地質(zhì)大學(xué), 2015.

GUI B. Research of 2D CSAMT tensor data inversion[D]. Beijing: China University of Geosciences, 2015. (In Chinese)

[17] PELTON W H, WARD S H, HALLOF P G, et al. Mineral discrimination and removal of inductive coupling with multifrequency IP[J]. Geophysics, 1978, 43(3): 588-609.

[18] 林紹忠.對稱逐步超松弛預(yù)處理共軛梯度法的改進迭代格式[J].數(shù)值計算與計算機應(yīng)用, 1997(12):266-270.

LIN S Z. Improved iterative format of symmetric successive over relaxation-preconditioned conjugated gradient method[J]. Journal of Numerical Methods and Computer Applications, 1997(12):266-270. (In Chinese)

[19] 譚捍東. 大地電磁法三維正反演問題研究[D]. 北京: 中國地質(zhì)大學(xué), 2000.

TAN H D. A study on 3d MT forward and inversion[D]. Beijing: China University of Geosciences, 2000. (In Chinese)

[20] 胡祥云, 李焱, 楊文采, 等. 大地電磁三維數(shù)據(jù)空間反演并行算法研究[J]. 地球物理學(xué)報, 2012, 55(12): 3969-3978.

HU X Y, LI Y, YANG W C, et al. Three-dimensional magnetotelluric parallel inversion using data space method[J]. Chinese Journal of Geophysics, 2012, 55(12): 3969-3978. (In Chinese)

[21] 林昌洪. 三維大地電磁共軛梯度反演算法研究[D]. 北京: 中國地質(zhì)大學(xué), 2004.

LIN C H. A study on 3d MT non-linear conjugate gradient inversion[D]. Beijing: China University of Geosciences, 2004. (In Chinese)

Astudyonparallelcomputationfortwo-dimensionalcomplexresistivityforwardmodeling