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

張志勇 , 譚捍東, 王春陽(yáng)， 汪 茂， 何中波， 梁 鵬

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

0 引言

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

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

1 MPI簡(jiǎn)介

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

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

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

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

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

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

2.1 二維數(shù)值模擬

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

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

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

(1)

傅里葉變換

(2)

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

(3)

(4)

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

(6)

2.2 Cole-Cole模型

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

(7)

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

2.3 有限單元計(jì)算公式

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

(8)

k1(i,j)=

(9)

(10)

(11)

k4(i,j)=

(12)

(13)

(14)

2.4 矩陣的壓縮存儲(chǔ)

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

2.5 求解方程組

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

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

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

(15)

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

(16)

2.6 正演正確性驗(yàn)證

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

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

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

圖4 二維棱柱體模型計(jì)算得到Ex實(shí)虛部對(duì)比曲線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ì)中存在一個(gè)長(zhǎng)200 m厚100 m的異常體，其頂部埋深為160 m，其中ρ0=10 Ω·m，m=0.5，c=0.5，τ=30 s。采用的裝置為偶極-偶極，沿x方向水平放置電偶極子，xz方向網(wǎng)格剖分為121×64。以下正演并行程序均在18個(gè)節(jié)點(diǎn)的集群上運(yùn)算，其計(jì)算環(huán)境均為：Linux操作系統(tǒng)，CPU-Intel(R) Xeon(R) CPU E5-2620 v2 @ 2.1GHz，內(nèi)存128G，開(kāi)發(fā)語(yǔ)言為Fortran，編譯器為ifort，并行環(huán)境為MPICH2。

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

各個(gè)發(fā)射源對(duì)應(yīng)頻率的計(jì)算是獨(dú)立的，基于發(fā)射源對(duì)應(yīng)的頻率來(lái)進(jìn)行并行計(jì)算的思路，其流程如圖6所示。

3.1.1 效率分析

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

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

程序運(yùn)行模式節(jié)點(diǎn)數(shù)量每個(gè)節(jié)點(diǎn)參與計(jì)算的進(jìn)程數(shù)各節(jié)點(diǎn)對(duì)應(yīng)進(jìn)程計(jì)算的源對(duì)應(yīng)的頻率數(shù)量網(wǎng)格大小時(shí)間/s并行加速比并行效率/%串行程序11105121×642207.39無(wú)無(wú)并行程序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可以看出，在每個(gè)節(jié)點(diǎn)只有一個(gè)進(jìn)程參與計(jì)算的前提下，當(dāng)參與并行計(jì)算的節(jié)點(diǎn)數(shù)增加時(shí)，加速比增大。當(dāng)有7個(gè)節(jié)點(diǎn)參與并行計(jì)算時(shí)，其計(jì)算速度約為串行程序計(jì)算速度的5.847倍，很大程度地降低了復(fù)電祖率法二維正演所需時(shí)間。但是，從表2中并行效率這一列看出，并不是參與并行的節(jié)點(diǎn)數(shù)越多，并行效率越高。當(dāng)參與并行的節(jié)點(diǎn)數(shù)從2增加到3時(shí)，并行效率增加，而當(dāng)參與并行的節(jié)點(diǎn)數(shù)從3增加到4(或節(jié)點(diǎn)數(shù)從5增加到6)時(shí)，并行效率反而降低，這是由于節(jié)點(diǎn)等待或節(jié)點(diǎn)間的數(shù)據(jù)交換等通信耗費(fèi)一些時(shí)間。

3.1.2 正確性驗(yàn)證

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

4 結(jié)論

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

圖7 正演并行結(jié)果與串行結(jié)果對(duì)比圖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

參考文獻(xiàn)：

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

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

[2] 楊振威, 嚴(yán)加永, 陳向斌. 頻譜激電法在安徽沙溪斑巖銅礦中的應(yīng)用[J]. 地球物理學(xué)進(jìn)展, 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é)報(bào), 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)蒙古那吉河鉛鋅礦探測(cè)中的應(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)法探測(cè)油氣藏的應(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è)計(jì)實(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] 都志輝, 李三立, 陳渝, 等. 高性能計(jì)算之并行編程技術(shù)——MPI并行程序設(shè)計(jì)[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]. 北京: 中國(guó)地質(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]. 北京: 中國(guó)地質(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]. 北京: 中國(guó)地質(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]. 北京: 中國(guó)地質(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] 林紹忠.對(duì)稱逐步超松弛預(yù)處理共軛梯度法的改進(jìn)迭代格式[J].數(shù)值計(jì)算與計(jì)算機(jī)應(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] 譚捍東. 大地電磁法三維正反演問(wèn)題研究[D]. 北京: 中國(guó)地質(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é)報(bào), 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]. 北京: 中國(guó)地質(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