周 翔 李長(zhǎng)偉,2 萬(wàn)文武 李忠乾3

（1.桂林理工大學(xué)地球科學(xué)學(xué)院；2.廣西隱伏金屬礦產(chǎn)勘查重點(diǎn)實(shí)驗(yàn)室）

大地電磁法是一種天然源電磁勘探方法，探測(cè)深度大，被廣泛應(yīng)用于區(qū)域構(gòu)造調(diào)查、石油勘探、地殼上地幔研究、深部找礦等領(lǐng)域［1-4］。隨著現(xiàn)代技術(shù)的發(fā)展，數(shù)值模擬成為了大地電磁正演的有效手段，常見(jiàn)的數(shù)值模擬方法包括積分方程法、有限差分法、有限體積法等。積分方程法具有較高的求解精度，但只適用于簡(jiǎn)單模型，不適合用于復(fù)雜模型［5］。有限差分法實(shí)現(xiàn)簡(jiǎn)單，但是只能采用結(jié)構(gòu)化網(wǎng)格，同樣不能用于復(fù)雜模型［6］；有限體積法是基于有限差分和有限元發(fā)展起來(lái)的，兼顧了有限差分和有限元的優(yōu)點(diǎn)，表達(dá)形式簡(jiǎn)單，適合處理帶任意地形的復(fù)雜模型，但是難以實(shí)現(xiàn)高精度的模擬［7-8］；有限單元法網(wǎng)格剖分靈活，可以處理復(fù)雜模型［9-12］。

譜元法是基于有限元和譜方法發(fā)展起來(lái)的一種數(shù)值模擬方法，其優(yōu)點(diǎn)在于一方面具有有限元處理復(fù)雜地形的能力，另一方面具有譜方法的高精度特性。

譜元法最早被Patera［13］提出，并應(yīng)用于流體力學(xué)當(dāng)中，同樣的模型和網(wǎng)格下，相對(duì)其他幾種方法計(jì)算量低，并且具有較高的計(jì)算精度。因?yàn)樽V元法的諸多優(yōu)點(diǎn)，國(guó)內(nèi)外學(xué)者對(duì)其展開(kāi)了一系列的研究。如Lim等［14］利用譜元法實(shí)現(xiàn)了三維腔體中電阻率正演成像。Lee［15］采用了混合階高斯—洛巴托—勒讓德（Gauss-Lobatto-Legendre，簡(jiǎn)稱(chēng)GLL）多項(xiàng)式的譜元法用于求解矢量電磁場(chǎng)波動(dòng)方程，實(shí)現(xiàn)了3D譜元法的數(shù)值模擬研究。黃鑫［16］采用基于規(guī)則六面體網(wǎng)格譜元法進(jìn)行頻率域航空電磁正演。Zhou等［17］將插值節(jié)點(diǎn)和積分節(jié)點(diǎn)選在同一GLL點(diǎn)上，解決了有限元在低頻時(shí)收斂速度慢和難以迭代求解線(xiàn)性方程組的問(wèn)題。陳漢波等［18］基于異常場(chǎng)的方法，結(jié)合混合階矢量基函數(shù)，利用不完全LU分解的IDR（s）迭代算法求解線(xiàn)性方程組，實(shí)現(xiàn)了任意各向異性介質(zhì)的海洋可控源三維譜元法正演。Pasquetti等［19］將三角譜元法和四邊形網(wǎng)格譜元法進(jìn)行了對(duì)比研究，積分節(jié)點(diǎn)采用最優(yōu)化點(diǎn)集Fekete節(jié)點(diǎn)，計(jì)算結(jié)果表明，三角網(wǎng)格譜元法的條件數(shù)隨階數(shù)的增長(zhǎng)快于四邊形網(wǎng)格譜元法。張帥胤等［20］采用一種三角形到四邊形的映射，對(duì)復(fù)雜區(qū)域上橢圓形方程的混合邊界問(wèn)題，建立三角單元的Legendre譜元法，應(yīng)用于若干不規(guī)則區(qū)域問(wèn)題的計(jì)算，通過(guò)數(shù)值算例驗(yàn)證該方法的有效性。方小姣等［21］基于譜元法實(shí)現(xiàn)了大地電磁二維數(shù)值模擬。Li等［22］精確推導(dǎo)了關(guān)于譜元的二次多項(xiàng)式。趙廷剛等［23］推導(dǎo)了擬譜格式雪比切夫—勒讓德（Chebyshev-Legendre）擬譜格式，并進(jìn)行了數(shù)值實(shí)驗(yàn)，取得良好的效果。張榮欣等［24］用雪比切夫多項(xiàng)式求解了均勻流場(chǎng)中的聲波問(wèn)題。朱昌允等［25］采用Chebyshev譜元方法結(jié)合并行計(jì)算，求解三維區(qū)域的Helmholtz方程問(wèn)題。劉玲［26］研究了基于高斯—洛巴托—雪比切夫（Gauss-Lobatto-Chebyshev，簡(jiǎn)稱(chēng)GLC）多項(xiàng)式譜元法的頻率域三維電磁正演模擬，利用GLC正交多項(xiàng)式構(gòu)建譜元矢量基函數(shù)，結(jié)合多項(xiàng)式性質(zhì)，推導(dǎo)單元積分解析表達(dá)式，實(shí)現(xiàn)高精度計(jì)算。崔海洋［27］實(shí)現(xiàn)了基于OpenMP大地電磁二維有限元正演模擬。韓思旭等［28］實(shí)現(xiàn)了基于CUDA并行計(jì)算的大地電磁二維有限元數(shù)值模擬，加速比達(dá)20多倍。劉慶等［29］采用了基于CPU（OpenMP）和GPU（CUDA）異構(gòu)并行處理的方式，實(shí)現(xiàn)了基于GPU并行的大地電磁二維正演。汪茂等［30］實(shí)現(xiàn)了基于大地電磁二維反演的MPI和CUDA并行算法的研究，加速比能達(dá)到2.15～3.09倍。杜偉［31］實(shí)現(xiàn)了基于GPU加速的直流電法的正演與大地電磁反演。謝悅等［32］研究了濃度對(duì)流擴(kuò)散方程的并行計(jì)算和MATLAB的高效實(shí)現(xiàn)方法。蘇輝等［33］研究了基于MATLAB的有限元方法的GPU加速。

本項(xiàng)目對(duì)大地電磁二維并行譜元法正演進(jìn)行了研究。從電磁場(chǎng)滿(mǎn)足的2階偏微分方程出發(fā)得到赫姆霍茲方程，采用高階基函數(shù)離散赫姆霍茲方程，用bicgstab求解線(xiàn)性方程組后得到電磁場(chǎng)和視電阻率的分布。為進(jìn)一步提高計(jì)算速度，采用OpenMP編程模式實(shí)現(xiàn)了多個(gè)頻點(diǎn)的并行計(jì)算。通過(guò)國(guó)際模型COMMEMI2D-1的計(jì)算對(duì)比，表明譜元法相較于有限元，能夠采用較少網(wǎng)格，到達(dá)同等的精度，有效地減少了計(jì)算時(shí)間；對(duì)常見(jiàn)的地電模型采用并行計(jì)算，實(shí)現(xiàn)頻點(diǎn)之間的并行，大大減少了正演計(jì)算的時(shí)間。

1 控制方程

大地電磁場(chǎng)滿(mǎn)足如下方程：

式中，E為電場(chǎng)強(qiáng)度；D為電位移矢量，本構(gòu)關(guān)系為D=εE；H為磁場(chǎng)強(qiáng)度；B為磁感應(yīng)強(qiáng)度，本構(gòu)關(guān)系為B=μH；J為電流密度，本構(gòu)關(guān)系為J=σE；ρ代表電流密度，其中ε，μ，σ分別表示介電常數(shù)、磁導(dǎo)率、電導(dǎo)率。

由式（1）和B=μH可得

采用時(shí)變因子e-iωt，公式（2）可寫(xiě)成

式中，ω為角頻率。式（5）可寫(xiě)成

再對(duì)式（7）求旋度，可得到赫姆霍茲方程

采用法拉第定律H=(1/iωμ?×E)求解磁場(chǎng)，并根據(jù)公式計(jì)算視電阻率

式中，Ex和Hy為電磁場(chǎng)的分量。

2 譜元法實(shí)現(xiàn)

為了用譜元法求解方程（8），將計(jì)算區(qū)域離散成四邊形單元。為了方便數(shù)值積分，通過(guò)坐標(biāo)映射，將離散單元轉(zhuǎn)換成標(biāo)準(zhǔn)參考單元，在參考單元內(nèi)利用GLL正交多項(xiàng)式進(jìn)行數(shù)值逼近，GLL正交多項(xiàng)式定義在［-1，1］區(qū)間。參考單元與物理單元的轉(zhuǎn)化關(guān)系見(jiàn)圖1。

物理坐標(biāo)(x,y)與參考坐標(biāo)(ξ,η)之間的映射關(guān)系為

式中，mx=xe+1-xe，my=ye+1-ye為物理單元2個(gè)方向的長(zhǎng)。

二重積分的變量替換公式為

其中，J為雅可比矩陣，

經(jīng)過(guò)物理單元與參考單元的轉(zhuǎn)化，將單元積分的計(jì)算轉(zhuǎn)化到(ξ,η)坐標(biāo)。

采用增加插值基函數(shù)階數(shù)，或者增加剖分單元個(gè)數(shù)來(lái)提高計(jì)算精度。在譜元法中，當(dāng)插值正交多項(xiàng)式階數(shù)選為N=1，與傳統(tǒng)的線(xiàn)性方法一致，當(dāng)插值基函數(shù)階數(shù)選為N=2，與二次插值的有限元方法一致。理論上來(lái)講，階數(shù)越高，效果越逼近，但同時(shí)占用內(nèi)存越大，計(jì)算機(jī)消耗越多。通常4到8階最好，為使效果最優(yōu)，內(nèi)存最小的情況下，本研究選擇了4階Legendre插值。

一維參考域的N階GLL插值多項(xiàng)式的表達(dá)式：

二維等參單元基函數(shù)由一階正交多項(xiàng)式構(gòu)成，參考單元下的基函數(shù)表達(dá)式如下。

其中，?rs(ξ,η)是等參單元的基函數(shù)Nξ，Nη表示在參考單元ξ，η方向基函數(shù)的階數(shù)，電場(chǎng)E(ξ,η)可寫(xiě)成

其中，Na=Nξ(Nη+1)+Nη(Nξ+1)是參考域節(jié)點(diǎn)數(shù)是參考域單元節(jié)點(diǎn)上的電場(chǎng)分量。

將微分方程（5）轉(zhuǎn)化成積分問(wèn)題，利用伽遼金殘差法處理電場(chǎng)矢量方程，令單元內(nèi)加權(quán)殘差為0，即

其中，vi為加權(quán)函數(shù)。

通過(guò)格林第一定律和狄利克雷邊界條件n×E|?L=0，?L是區(qū)域的邊界，式（18）可改寫(xiě)成

選擇權(quán)函數(shù)為插值基函數(shù)??，得到單元內(nèi)的方程

對(duì)式（20）可改寫(xiě)成矩陣形式

式中，Se為剛度矩陣；Oe為質(zhì)量矩陣；k為單元總數(shù)。

大地電磁正演需要計(jì)算多個(gè)頻點(diǎn)上的電磁響應(yīng)，以得到地下由淺至深的電阻率變化情況。每個(gè)頻點(diǎn)都需要進(jìn)行譜元求解，為了提高計(jì)算效率，考慮了基于頻點(diǎn)的并行計(jì)算策略。將計(jì)算任務(wù)按頻點(diǎn)進(jìn)行劃分，利用OpenMP并行機(jī)制為每一個(gè)頻點(diǎn)分配一個(gè)section，等待所有頻點(diǎn)任務(wù)結(jié)束，將結(jié)果匯總輸出。其中，并行加速比=單機(jī)的串行執(zhí)行時(shí)間/N個(gè)進(jìn)程的并行計(jì)算時(shí)間，并行效率=并行加速比/并行的進(jìn)程個(gè)數(shù)。并行計(jì)算流程如圖2所示。

3 數(shù)值計(jì)算結(jié)果與分析

基于上述分析，采用MATLAB平臺(tái)編寫(xiě)了二維大地電磁并行譜元法的程序。

3.1 算例1

算例1采用了國(guó)際模型中COMMEMI2D-1模型對(duì)本算法進(jìn)行驗(yàn)證，模型示意如圖3所示，中間有1個(gè)低阻的異常體，低阻異常體模型位于X方向-0.5～0.5 km、Y方向0.25～2.25 km，設(shè)電阻率ρ2=0.5Ω·m，背景電阻率ρ1=100Ω·m。計(jì)算區(qū)域設(shè)置為X方向從-10～10 km，Y方向?yàn)?～4 km。采用譜元法計(jì)算得到的視電阻曲線(xiàn)如圖4所示，可以看出，計(jì)算結(jié)果與IE方法的計(jì)算結(jié)果相吻合，最大誤差不超過(guò)5%。

在本算例中，將本方法與有限元進(jìn)行對(duì)比，如表1所示。譜元采用20乘20的網(wǎng)格，有限單元法采用60×80的網(wǎng)格，二者到達(dá)了相同的精度。在同等精度下，SEM相較于FEM剖分更少的網(wǎng)格、用了更少的計(jì)算時(shí)間。

?

3.2 算例2

算例2為地塹模型（圖5），地塹電阻率ρ2=100Ω·m，在地塹上面有1層覆蓋層，電阻率ρ1=10Ω·m。設(shè)置空氣層電阻率為1010Ω·m，凹陷部分為1個(gè)倒梯形，梯形上邊長(zhǎng)為4 km，梯形下邊長(zhǎng)為2 km，梯形高為2 km。在頻率范圍0.01～1 Hz內(nèi)進(jìn)行正演計(jì)算，繪制了X方向-3～3 km等值線(xiàn)圖。地塹視電阻率等值線(xiàn)見(jiàn)圖6。從視電阻率圖可以看出，地塹輪廓較為明顯，很好地反映了輪廓的邊界，其低阻異常的范圍能夠反映模型的基本特征，這與正演所設(shè)定的地電模型完全吻合。地塹的結(jié)果進(jìn)一步驗(yàn)證了本算法的正確性。

表2為多個(gè)頻率情況下地塹模型串并行計(jì)算時(shí)間對(duì)比。與串行相比，并行更慢的原因是程序需要花費(fèi)時(shí)間在并行管理，單個(gè)線(xiàn)程沒(méi)有體現(xiàn)并行的意義。2個(gè)線(xiàn)程的加速比為1.659 8倍，當(dāng)采用6個(gè)線(xiàn)程的時(shí)候，加速比可以達(dá)到2.561 1倍。圖7是OpenMP并行效率，可以看出，隨著線(xiàn)程數(shù)的增加，加速比增加，計(jì)算時(shí)間減少，但是加速比與線(xiàn)程數(shù)不是同比例增加（圖8），即并行效率降低了，這是因?yàn)殡S著線(xiàn)程數(shù)的增加，線(xiàn)程之間的通訊量和并行管理開(kāi)銷(xiāo)增加。

?

4 結(jié)論

將譜元法應(yīng)用于大地電磁二維正演，開(kāi)發(fā)了算法程序，通過(guò)數(shù)值模擬進(jìn)行測(cè)試。國(guó)際模型COMMEMI2D-1驗(yàn)證了算法的有效性和準(zhǔn)確性，與有限元對(duì)比，譜元法在同等精度下只需剖分更少的網(wǎng)格，用更少的計(jì)算時(shí)間；利用地電模型測(cè)試OpenMP并行效果，基于頻點(diǎn)的并行策略能大大提高正演速度，可以為反演的快速實(shí)現(xiàn)奠定基礎(chǔ)。