基于CUDA架構(gòu)并行算法的帶地形AMT二維反演實(shí)現(xiàn)與應(yīng)用

韓思旭， 陳小斌， 陳衛(wèi)營， 羅 強(qiáng)， 宋婉婷

(1.廣東省地球物理探礦大隊(duì)， 廣州 510800; 2.廣東省地質(zhì)物探工程勘察院， 廣州 510800; 3.中國地震局地殼應(yīng)力研究所， 北京100085; 4.中國科學(xué)院地質(zhì)與地球物理研究所， 中國科學(xué)院礦產(chǎn)資源研究重點(diǎn)實(shí)驗(yàn)室， 北京 100029)

音頻大地電磁法(audio-frequency magnetotelluric method，AMT)通過觀測和處理天然電磁場信號獲取地下一定范圍內(nèi)的電性結(jié)構(gòu)分布，具有頻帶寬(0.1 Hz～100 kHz)、探測深度大、對二維構(gòu)造反應(yīng)明顯、施工輕便等優(yōu)勢，已成為礦產(chǎn)勘探、油氣調(diào)查、工程勘查等領(lǐng)域的重要地球物理探測工具[1]。反演是AMT資料處理過程中的關(guān)鍵步驟。目前，二維反演是AMT數(shù)據(jù)處理實(shí)際應(yīng)用中最為有效和普遍的手段，如何提高二維反演精度與速度是保障探測結(jié)果解釋可靠性和提升數(shù)據(jù)處理效率的關(guān)鍵問題[2]。

AMT數(shù)據(jù)的二維反演和一維反演在反演原理和反演方法上并無本質(zhì)區(qū)別, 只是二維反演處理的不再只是單個(gè)測點(diǎn)上的數(shù)據(jù), 而是包含整條測線(垂直于構(gòu)造走向) 上所有測點(diǎn)的數(shù)據(jù)，并且還需同時(shí)考慮介質(zhì)電性在縱向和橫向兩個(gè)方向上的變化，因此二維反演需要擬合的數(shù)據(jù)量大大增加[3]。此外，相較于一維反演，二維反演處理的觀測參數(shù)更多, 如TE(transverse-electric)和TM(transverse-magnetic)兩個(gè)極化方向的電阻率和相位，或者由垂直磁場計(jì)算而來的傾子Tipper等[4]，這就使得模型數(shù)據(jù)相對參數(shù)的偏導(dǎo)數(shù)矩陣(雅可比矩陣)的階數(shù)較大, 造成計(jì)算量大幅上升。因此在處理大網(wǎng)格、多頻率模型時(shí)，二維反演需要較長的計(jì)算時(shí)間。二維反演的主要計(jì)算量來自于二維正演計(jì)算及偏導(dǎo)數(shù)矩陣的計(jì)算。為此，眾多學(xué)者提出了多種快速二維反演方法，從二維正演或雅可比矩陣計(jì)算入手，發(fā)展了降低計(jì)算量、提升計(jì)算速度的方法[5-8]，極大地推動(dòng)了AMT二維反演的實(shí)用化。

隨著高性能計(jì)算機(jī)特別是并行計(jì)算技術(shù)的蓬勃發(fā)展，利用并行計(jì)算能力成為實(shí)現(xiàn)高精度、快速度AMT數(shù)據(jù)反演的一種可行手段[9-10]。最先發(fā)展起來的并行計(jì)算方式主要基于多核CPU(central processing unit)進(jìn)行，其中代表性的是基于MPI(message passing interface)環(huán)境的并行計(jì)算，大量的研究成果顯示出并行計(jì)算在電磁正反演中的良好加速效果[11-15]。但是，受CPU核數(shù)限制，這種并行計(jì)算的實(shí)現(xiàn)往往需要借助超級計(jì)算機(jī)或者計(jì)算機(jī)集群，造成計(jì)算成本較高[16]，同時(shí)CPU程序邏輯的復(fù)雜度也限定了程序執(zhí)行的指令并行性。而GPU(graphics processing unit)的最大特點(diǎn)是它擁有超多計(jì)算核心，且每個(gè)核心都可以模擬一個(gè)CPU的計(jì)算功能。目前GPU已經(jīng)發(fā)展為多線程多核處理器，具有強(qiáng)大的并行計(jì)算能力，近年來將GPU并行計(jì)算能力與實(shí)際的科學(xué)工程計(jì)算需求相結(jié)合已表現(xiàn)出巨大的潛力[17-19]。尤其是2006年NVIDIA公司為實(shí)現(xiàn)GPU通用計(jì)算提出編程統(tǒng)一計(jì)算設(shè)備架構(gòu)(compute unified device architecture，CUDA)后，使得在科學(xué)計(jì)算中一些過去必須由高性能計(jì)算集群或超級計(jì)算機(jī)處理的計(jì)算任務(wù)可被轉(zhuǎn)換到單個(gè)小型計(jì)算機(jī)上進(jìn)行，大大降低了硬件成本和運(yùn)行成本。已有研究成果表明，基于CUDA架構(gòu)的并行算法可在AMT正反演加速中發(fā)揮出色的作用[20-21]。

為此，在CUDA 架構(gòu)下開展帶地形AMT數(shù)據(jù)的二維反演并行算法研究，旨在利用GPU強(qiáng)大的計(jì)算能力及并行計(jì)算技術(shù)實(shí)現(xiàn)高精度、快速度的AMT數(shù)據(jù)二維反演，并為后續(xù)開展并行三維反演建立基礎(chǔ)。

1 AMT二維正演與反演

1.1 AMT二維有限元正演模擬

對于一般二維各向同性介質(zhì)，若取坐標(biāo)系x軸方向?yàn)榈刭|(zhì)體走向(即電性沿此方向不變)，z軸垂直向下，則與之對應(yīng)的二維MT邊界值問題的變分問題[22]為

(1)

TE和TM兩種模式的邊界條件基本一致，僅在上邊界AB位置的選擇上有所不同。對于TE模式，上邊界位于空氣中，并離地表足夠遠(yuǎn)，對于TM模式上邊界則位于地表。

為適應(yīng)地形變化，本研究使用雙線性任意四邊形單元對區(qū)域Ω進(jìn)行剖分。將求解區(qū)域剖分為多個(gè)四邊形單元，在離散的每個(gè)單元上進(jìn)行線性化處理，最后合成總體系數(shù)矩陣為

(2)

式(2)中：K=K1e+K2e+K3e為有限元總體系數(shù)矩陣。令式(2)的變分值為零，得線性方程組為

Ku=0

(3)

代入頂邊界條件u|AB=1.0，當(dāng)求解TE模式時(shí)AB位于空氣中頂邊界，當(dāng)求解TM模式時(shí)AB位于地表。由于方程中的矩陣K是稀疏對稱正定的，為了節(jié)約內(nèi)存空間按變帶寬存儲，本文中采用稀疏矩陣快速回代的Cholesky分解法求解方程。

1.2 ARIA二維反演

自適應(yīng)正則化反演方法(adaptive regularized inversion algorithm ,ARIA)是陳小斌等[23]提出的一種反演方法，相較于傳統(tǒng)的Occam算法，ARIA采用簡單有效的正則化因子求取算法，其收斂時(shí)間更快，擬合效果基本相同。其總目標(biāo)函數(shù)可歸結(jié)為

Φ=Φd+λΦm→min

(4)

式(4)中：Φd=[d-A(m)]TWd[d-A(m)]為數(shù)據(jù)目標(biāo)函數(shù)；Φm=mTCxm+mTCzm為模型目標(biāo)函數(shù)；λ為數(shù)據(jù)目標(biāo)函數(shù)和模型粗糙度的正則化調(diào)節(jié)因子，即

Φ=[d-A(m)]TWd[d-A(m)]+

λ[mTCxm+mTCzm]

(5)

根據(jù)目標(biāo)函數(shù)極小化原則，可得反演矩陣方程為

WdJΔm+λ(Cx+Cz)Δm=WdΔd-

λ(Cx+Cz)m

(6)

式(6)中：J為正演數(shù)據(jù)對模型的偏導(dǎo)數(shù)矩陣。

這樣每次迭代模型的更新為mi+1=mi+Δm。

為保證反演過程中數(shù)據(jù)目標(biāo)函數(shù)值不受數(shù)據(jù)誤差的影響，對數(shù)據(jù)方差進(jìn)行了規(guī)范化處理，使得正則化因子的取值與數(shù)據(jù)誤差的大小基本無關(guān)，即

(7)

(8)

反演中每一步計(jì)算都需要極小化正則化因子，此過程需要重復(fù)的正則化求解反演方程組，較為耗時(shí)，為了減少計(jì)算量本文根據(jù)數(shù)據(jù)目標(biāo)函數(shù)和模型目標(biāo)函數(shù)的關(guān)系進(jìn)行自適應(yīng)正則化因子調(diào)節(jié)[23]，關(guān)系式為

(9)

2 CUDA架構(gòu)下AMT二維正反演并行計(jì)算實(shí)現(xiàn)策略

AMT正反演中涉及多個(gè)頻率、多個(gè)測點(diǎn)和大量網(wǎng)格計(jì)算單元。傳統(tǒng)串行算法是將每個(gè)頻點(diǎn)、每個(gè)測點(diǎn)、每個(gè)網(wǎng)格單元依次計(jì)算，當(dāng)測點(diǎn)、頻點(diǎn)、網(wǎng)格單元較多時(shí)勢必涉及巨大計(jì)算量，造成計(jì)算速度慢。如果將多個(gè)測點(diǎn)、頻點(diǎn)、網(wǎng)格單元的計(jì)算分成多個(gè)模塊的線程來進(jìn)行獨(dú)立執(zhí)行，這樣就可以提高計(jì)算速度，從而縮短計(jì)算時(shí)間，這便是并行計(jì)算的主要思想。

在AMT二維反演計(jì)算中，靈敏度的計(jì)算、反演矩陣方程組的正則化求解、模型響應(yīng)部分為主要的耗時(shí)計(jì)算部分，因此是需要并行的部分。其中AMT靈敏度的計(jì)算和模型的響應(yīng)部分都需要對多個(gè)頻率進(jìn)行單獨(dú)計(jì)算，且各頻率對應(yīng)的電磁場值間是相互獨(dú)立的；反演矩陣方程組的正則化求解過程中需要對矩陣A(M×N)進(jìn)行正則化得到矩陣B(M×N)，且這個(gè)正則化的過程中矩陣B每個(gè)元素的計(jì)算是相互獨(dú)立的，根據(jù)這些特點(diǎn)可以將程序按頻率或者矩陣的行數(shù)和列數(shù)進(jìn)行分粒度，將每個(gè)頻點(diǎn)或矩陣的每數(shù)元素分配到各個(gè)進(jìn)程同時(shí)進(jìn)行計(jì)算，并行執(zhí)行。

二維反演計(jì)算中模型靈敏度矩陣的計(jì)算、反演方程的正則化處理以及求解部分主要的耗時(shí)計(jì)算部分，因此同樣的將其定為需要并行的部分。在AMT二維反演靈敏度矩陣計(jì)算中，各個(gè)測點(diǎn)及各個(gè)模型所涉及的計(jì)算都是相互獨(dú)立的；反演方程的正則計(jì)算中，正則化后的每個(gè)元素計(jì)算都是相互獨(dú)立的；正則化后的矩陣方程進(jìn)行不完全三角分解的預(yù)條件共軛梯度法(incomplete cholesky factorization preconditioned conjugate gradient，ICCG)求解迭代時(shí)將涉及到大量的矩陣、向量的乘積計(jì)算，也是可以采取并行計(jì)算加速的。因此在二維反演靈敏度矩陣計(jì)算中將程序按測點(diǎn)、縱向網(wǎng)格分塊數(shù)量、橫向網(wǎng)格分塊數(shù)量進(jìn)行粒度分布，將各個(gè)測點(diǎn)針對各個(gè)有限單元網(wǎng)格的偏導(dǎo)數(shù)計(jì)算進(jìn)程同時(shí)進(jìn)行，并行執(zhí)行；在反演方程的正則計(jì)算中，將程序按模型數(shù)量分成多個(gè)線程塊，在同一個(gè)線程塊中又將其分成多個(gè)線程；在正則化后的矩陣方程進(jìn)行ICCG求解迭代時(shí)每一步的矩陣與向量的乘積計(jì)算使用并行化處理。

本文中程序?qū)崿F(xiàn)的過程使用的是主從并行模式，即將程序的執(zhí)行分主函數(shù)進(jìn)程和子函數(shù)進(jìn)程，主進(jìn)函數(shù)程控制全局的數(shù)據(jù)結(jié)構(gòu)，主要負(fù)責(zé)子進(jìn)程任務(wù)數(shù)量的劃分、派發(fā)、結(jié)果回收及輸出結(jié)果；子函數(shù)進(jìn)程主要負(fù)責(zé)執(zhí)行由主進(jìn)程分發(fā)指定的任務(wù)并把結(jié)果回傳給主函數(shù)進(jìn)程。其中正演計(jì)算部分具體的實(shí)現(xiàn)過程是：主進(jìn)程讀入網(wǎng)格剖分?jǐn)?shù)據(jù)及網(wǎng)格單元電阻率值，并將其傳送給子進(jìn)程，按頻率劃分給子進(jìn)程各自計(jì)算模型的響應(yīng)，并將計(jì)算結(jié)果傳給主進(jìn)程，主進(jìn)程整合所有計(jì)算的結(jié)果，并將其輸出；反演靈敏度計(jì)算部分具體的實(shí)現(xiàn)過程是：主進(jìn)程讀入正演所需的數(shù)據(jù)及模型數(shù)據(jù)，并將其傳送給子進(jìn)程，按測點(diǎn)數(shù)進(jìn)行劃分子進(jìn)程各自計(jì)算測點(diǎn)位置的電磁場值關(guān)于模型偏導(dǎo)數(shù)，以及按測點(diǎn)數(shù)和有限元模型單元數(shù)進(jìn)行劃分子進(jìn)程各自計(jì)算測點(diǎn)響應(yīng)關(guān)于模型單元的偏導(dǎo)數(shù)，并將計(jì)算結(jié)果傳給主進(jìn)程，主進(jìn)程整合所有計(jì)算的結(jié)果，并將其輸出。圖1所示為二維正演并行計(jì)算簡化流程圖，圖2所示為二維反演模型靈敏度計(jì)算并行計(jì)算簡化流程圖。在并行計(jì)算中，CUDA并行計(jì)算需要CPU來對GPU分配計(jì)算線程，由GPU來執(zhí)行計(jì)算，最終再將計(jì)算結(jié)果返回CPU。GPU作為設(shè)備來完成設(shè)備(Device)代碼命令部分，這部分代碼中CPU分配的進(jìn)程與GPU的計(jì)算眾核對應(yīng)位置關(guān)系是隨機(jī)的，并且獨(dú)立完成，CPU作為主機(jī)(Host)對GPU下達(dá)線程分配命令，并等待GPU計(jì)算返回結(jié)果。

圖1 AMT二維正演CPU+GPU并行計(jì)算流程圖

圖2 基于CPU+GPU的偏導(dǎo)數(shù)矩陣并行計(jì)算流程圖

3 程序準(zhǔn)確性驗(yàn)證

3.1 一維層狀模型

為了驗(yàn)證程序算法的正確性，首先設(shè)計(jì)了一個(gè)8層一維模型，各層電阻率ρ和厚度h分別為ρ=[10、200、20、100、2、500、50、10] Ω·m，h=[100、500、200、1 000、1 000、5 200、10 000] m。正演計(jì)算頻率為55個(gè)，最高頻率為100 000 Hz，最低頻率為0.000 01 Hz，每個(gè)數(shù)量級按6個(gè)頻點(diǎn)對數(shù)等間隔取值。將正演結(jié)果(電阻率和相位)與一維解析解進(jìn)行對比, 如圖3所示。結(jié)果表明，該二維正演程序所獲結(jié)果與一維解析結(jié)果基本完全擬合，驗(yàn)證了該算法及程序是可行的且是準(zhǔn)確的。

圖3 8層模型二維正演計(jì)算結(jié)果與解析結(jié)果對比

然后利用本文所提出的反演方法對上述二維正演結(jié)果進(jìn)行一維、二維反演，結(jié)果如圖4所示。由反演結(jié)果可以看出，一維、二維反演結(jié)果基本一致，都能較好地將數(shù)據(jù)恢復(fù)到真實(shí)模型附近。說明一維反演算法是可靠的，基于CUDA并行計(jì)算的二維反演也是正確的。

圖4 8層模型1D和2D反演結(jié)果對比圖

3.2 二維模型

接下來，通過三個(gè)復(fù)雜程度依次遞增的二維模型驗(yàn)證本程序的準(zhǔn)確性，結(jié)果如圖5～圖7所示?？梢钥闯?，本文算法對于單獨(dú)TE反演時(shí)對低阻體的形態(tài)恢復(fù)較好，單獨(dú)TM反演時(shí)對高阻體的形態(tài)恢復(fù)較好，這與AMT中TE、TM兩種模式對高低阻體反應(yīng)靈敏度特性有關(guān)。但總體來說無論是對于不帶地形的簡單二維模型(圖5)，還是帶地形的二維簡單(圖6)和較復(fù)雜模型(圖7)，本文算法都能很好地復(fù)原模型的真實(shí)結(jié)構(gòu)，特別是當(dāng)采用TE+TM模式反演時(shí)，反演結(jié)果對模型中的高、低阻異常的形態(tài)和位置都做出了準(zhǔn)確地刻畫。由此可說明，該程序在處理二維模型的正反演中具有較高的精度。

圖5 簡單二維模型反演結(jié)果

圖6 帶地形的簡單二維模型反演結(jié)果

圖7 帶地形的復(fù)雜二維模型反演結(jié)果

4 并行加速及壓縮效果對比

為了說明CDUA并行計(jì)算方法在正演計(jì)算速度上的優(yōu)越性，與基于CPU串行計(jì)算的有限單元法計(jì)算程序進(jìn)行了對比，同時(shí)測試了5種不同網(wǎng)格剖分在GPU硬件環(huán)境下的執(zhí)行效率(硬件環(huán)境為Intel? core(TM) i7-9750 CPU 2.60GHz及NvidiaGeForce GTX1660 Ti)，測試結(jié)果如表1、表2所示。從表中可以看出無論在何種測點(diǎn)、頻率、網(wǎng)格配置下，GPU計(jì)算速度都要比CPU的計(jì)算速度快，且計(jì)算網(wǎng)格越大加速比越高，加度比計(jì)算公式為

表1 靈敏度矩陣串行與并行計(jì)算效率對比

表2 矩陣正則化非壓縮串行與壓縮并行計(jì)算效率對比

(10)

由此可以說明本文所提的并行反演方法可行，在同樣的測點(diǎn)、頻點(diǎn)、模型網(wǎng)格情況下可以更快地得到反演結(jié)果，在實(shí)際應(yīng)用中就能提高效率，從而可以開展更多測點(diǎn)的長剖面測量工作。

5 應(yīng)用實(shí)例

上述數(shù)值模擬結(jié)果表明，采用本文提出的并行計(jì)算方案，可以極大地提升AMT數(shù)據(jù)的二維反演效率。最后，以云南某金屬礦的AMT探測數(shù)據(jù)反演為例，說明本文提出的二維反演方法在實(shí)際應(yīng)用中的效果。如圖8所示，測量區(qū)位于云南永平某銅鈷多金屬礦區(qū)，區(qū)內(nèi)發(fā)育的褶皺主要為向斜，局部發(fā)育一些小型褶皺。向斜緊靠F4礦化帶西側(cè)展布，軸向30°，軸部地層為花開左組下段(J2h1)，兩翼為漾江組(J1y)及麥初箐組上段(T3m2)，西翼被F4斷層所切，出現(xiàn)J2h1與T3m2接觸，并出現(xiàn)強(qiáng)烈蝕變及銅礦化。礦段內(nèi)礦體主要賦存于向斜西翼三疊系與侏羅系過渡部位，主要受F5號帶控制，分布于礦段中部，為一組平行斷裂帶。斷裂帶產(chǎn)狀：260°～330°∠40°～65°，傾角呈現(xiàn)南陡北緩的特征，北側(cè)向西轉(zhuǎn)折，走向?yàn)楸蔽?，切過侏羅紀(jì)和三疊紀(jì)地層。斷裂帶的成分主要為蝕變碎裂巖，銅礦化普遍發(fā)育，主要銅礦物有砷黝銅礦、黃銅礦、斑銅礦、輝銅礦、孔雀石、銅藍(lán)等。

1為侏羅系中統(tǒng);2為侏羅系下統(tǒng);3為三疊系;4為三疊系;5為出露工業(yè)礦體;6為出露低品位礦體;7為地質(zhì)界線;8為實(shí)測斷層;9為推測斷層;10為背斜;11為向斜

為獲得目前正在開采的F5號帶所控制礦體南沿以及F2、F3、F4號帶的深部形態(tài)，對礦山南沿的礦產(chǎn)資源潛力做出初步評價(jià)，在測區(qū)實(shí)施了AMT測量。共布置2條AMT測線，點(diǎn)距20 m，觀測頻率為1.0～32 000 Hz，每數(shù)量級10頻點(diǎn)對數(shù)等間隔分布，TM標(biāo)量測量，使用儀器為重慶地質(zhì)儀器廠生產(chǎn)的CLEM陣列電磁法系統(tǒng)。

以L01線數(shù)據(jù)為例進(jìn)行并行二維反演，并對反演結(jié)果進(jìn)行地質(zhì)資料解釋與驗(yàn)證。L01線長度為2 400 m，共包含121個(gè)TM標(biāo)量測量點(diǎn)。反演初始模型設(shè)置橫向129塊，每塊20 m，頂層厚度10 m，往深部按1.1倍逐步增加，共34層。設(shè)置好參數(shù)之后，程序根據(jù)具體測點(diǎn)的實(shí)際坐標(biāo)分布適當(dāng)?shù)刈詣?dòng)調(diào)節(jié)橫向?qū)挾?，使得測點(diǎn)盡量落在模型塊的橫向中心位置；縱向模型位置分布根據(jù)地形起伏自動(dòng)貼合，一方面適應(yīng)地形起伏，另一方面使得模型往深部逐步變平緩。初始模型電阻率設(shè)置為100 Ω·m，設(shè)置好反演參數(shù)之后進(jìn)行自動(dòng)迭代反演。在Win10 64位，i7-9750H CPU 2.60 GHz，NVIDIA GeForce GTX1660 GPU配置上運(yùn)行，迭代15次，擬合殘差3.37%，用時(shí)約30 min，結(jié)果如圖9所示。

由圖9可以看出，受斷裂帶影響測區(qū)地層的電性結(jié)構(gòu)較為錯(cuò)亂，淺層地表電性特征與地質(zhì)揭露的實(shí)際資料基本相符，在F5號帶附近實(shí)施的鉆孔所揭示的中淺部地層結(jié)構(gòu)也能與本結(jié)果相對應(yīng)。測線中部和南部出現(xiàn)的兩處大面積的低阻異常區(qū)域是F3和F4號斷裂帶的反映，它們在水平位置上能夠與實(shí)際情況相吻合，并顯示出向深部繼續(xù)延伸的趨勢。該二維反演結(jié)果很好地刻畫了地下電性結(jié)構(gòu)的形態(tài)，表明本文提出的并行二維反演算法具有較高的反演精度，且具有可觀的計(jì)算效率，在實(shí)測數(shù)據(jù)的反演處理中可以獲得很好的效果。

圖9 L01線二維反演高程-電阻率斷面圖

6 結(jié)論

正演、靈敏度矩陣計(jì)算、反演方程組正則化是影響MT數(shù)據(jù)二維反演效率的主要部分，而它們具有很好的并行性，因此利用GPU強(qiáng)大的計(jì)算能力及并行計(jì)算技術(shù)可極大地提高AMT二維反演程序的運(yùn)行效率。在CUDA 架構(gòu)下實(shí)現(xiàn)了基于有限元和自適應(yīng)正則化反演算法的帶地形AMT數(shù)據(jù)的二維反演并行計(jì)算，并通過理論模型驗(yàn)證了方法的有效性和準(zhǔn)確性。數(shù)值模擬結(jié)果表明，該并行算法可顯著提升靈敏度矩陣計(jì)算和反演方程組正則化求解的效率，最高加速比可達(dá)10倍以上。通過對數(shù)值模擬結(jié)果以及野外實(shí)測數(shù)據(jù)的處理表明，本文所提出的AMT二維并行反演方法在保證高精度的前提下具有很高的計(jì)算效率，可為今后研究及實(shí)際生產(chǎn)應(yīng)用提供一種有力工具。