基于無單元局部弱式法的2.5維直流電阻率正演

麻昌英，閆玲玲，姚振岸*，柳建新，趙文學(xué)，周聰

（1.江西省防震減災(zāi)與工程地質(zhì)災(zāi)害探測工程研究中心（東華理工大學(xué)），江西南昌 330013； 2.中南大學(xué)有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點(diǎn)實(shí)驗(yàn)室，湖南長沙 410083； 3.東華理工大學(xué)地球物理與測控技術(shù)學(xué)院，江西南昌 330013）

0 引言

目前直流電阻率法正演模擬主要采用積分方程法、有限差分法和有限單元法（Finite element method，F(xiàn)EM）。相比積分方程法和有限差分法，F(xiàn)EM 的優(yōu)勢是適應(yīng)于任意地形和任意復(fù)雜的地質(zhì)目標(biāo)的模擬，具有適應(yīng)性強(qiáng)、精度高、計(jì)算效率較高等特點(diǎn)，是目前直流電阻率法和電磁法正演模擬中主要的數(shù)值模擬方法［1-10］。近年來，在常規(guī)FEM 基礎(chǔ)上，基于非結(jié)構(gòu)化單元和自適應(yīng)加密策略的自適應(yīng)FEM 得到了快速發(fā)展［11-24］。自適應(yīng)FEM 利用圖形幾何學(xué)相關(guān)單元剖分器對任意分布節(jié)點(diǎn)進(jìn)行網(wǎng)格剖分。非結(jié)構(gòu)化單元適用于起伏地形和地下形態(tài)復(fù)雜目標(biāo)體的正演模擬，可在局部區(qū)域自適應(yīng)地加密單元，模擬過程更靈活，具有較高的模擬精度?；诜墙Y(jié)構(gòu)化單元的FEM 可很好地模擬復(fù)雜起伏地形和地下復(fù)雜目標(biāo)體，但不規(guī)則單元剖分不可避免。由于使用預(yù)定義的、基于節(jié)點(diǎn)連接信息的單元進(jìn)行模擬，不可避免地需要進(jìn)行網(wǎng)格單元剖分［25-28］。隨著計(jì)算機(jī)硬件的快速發(fā)展，以及野外實(shí)際生產(chǎn)中待開發(fā)區(qū)域地質(zhì)條件日趨復(fù)雜，對資料解釋及勘探精細(xì)化的要求不斷提高，尤其對復(fù)雜地電模型的模擬靈活性和模擬精度需求更高。發(fā)展創(chuàng)新型的精細(xì)化、高靈活性和更高模擬精度的地球物理勘探數(shù)據(jù)解譯系統(tǒng)，已成為地球物理勘探領(lǐng)域的研究熱點(diǎn)。

相比傳統(tǒng)數(shù)值模擬方法，基于全局弱式的無單元Galerkin 法（Element-Free Galerkin Method，EFGM）無需進(jìn)行基于節(jié)點(diǎn)連接信息的單元剖分［29］，可方便地根據(jù)實(shí)際需求任意布置和加密節(jié)點(diǎn)，對任意復(fù)雜地電模型及局部加密等問題具有很強(qiáng)的適應(yīng)性，是一種高靈活性、高精度的正演模擬方法［30-33］。在地球物理領(lǐng)域，無單元Galerkin 法在地震波場［34-40］、探地雷達(dá)［27，41-42］、大地電磁［43-47］及重力勘探［48-49］等正演模擬研究中取得了良好的研究成果。在直流電阻率法正演模擬中，麻昌英等［50-51］、Ma 等［52］分別利用徑向基點(diǎn)插值（RPIM）形函數(shù)和移動(dòng)最小二乘（MLS）形函數(shù)開展了2.5維直流電阻率EFGM 正演，實(shí)現(xiàn)了任意復(fù)雜形體、地形起伏、局部節(jié)點(diǎn)加密的直流電阻率法的正演模擬。研究成果表明，EFGM 能有效地解決傳統(tǒng)數(shù)值模擬方法受單元約束的局限性，與傳統(tǒng)數(shù)值模擬方法相輔相成，具有較高的靈活性、適應(yīng)性和模擬精度，尤其在任意復(fù)雜地電模型的正演中效果很好，在地球物理領(lǐng)域具有良好的發(fā)展前景。

EFGM 是一種基于全局弱式的無單元法，需要借助于全局域剖分背景單元進(jìn)行積分計(jì)算，對單元仍有一定的依賴性，一定程度上損失了無單元的屬性。為進(jìn)一步減小對單元的依賴性，避免使用全局域背景單元積分，基于局部弱式的無單元法得以快速發(fā)展。Atluri 等［53］在彈性靜力學(xué)局部邊界積分方程的基礎(chǔ)上，提出了無網(wǎng)格局部Petrov-Galerkin 法（Mesh-Free Local Petrov-Galerkin Method，MLPGM）。該方法利用加權(quán)余量法，將全局域積分轉(zhuǎn)化為節(jié)點(diǎn)局部域積分。MLPGM 應(yīng)用于大地電磁場的二維正演模擬［28，54］，展現(xiàn)了較強(qiáng)的靈活性和適應(yīng)性及較高的計(jì)算精度。然而，與EFGM 引入權(quán)函數(shù)不同，MLPGM 采用加權(quán)余量法進(jìn)行方程離散，使用了更多的模擬參數(shù)；同時(shí)，離散方程按節(jié)點(diǎn)編號進(jìn)行組裝，系數(shù)矩陣稀疏但不對稱。Melenk 等［55］和Babu?ka 等［56］提出了單位分解有限單元法（Partition of Unity Finite Element Method，PUFEM）和單位分解法（Partition of Unity Method，PUM），并進(jìn)行了嚴(yán)格的數(shù)學(xué)證明和論述；Carpinteri等［57］將PUM 積分應(yīng)用于EFGM，將全局域積分轉(zhuǎn)化為節(jié)點(diǎn)局部域積分，實(shí)現(xiàn)了基于局部弱式的無單元法。隨后，該方法在各個(gè)領(lǐng)域得到深入研究和廣泛的應(yīng)用［58-62］。基于局部弱式的無單元法與EFGM 類似，計(jì)算過程基本一致，均繼承了FEM 的一些優(yōu)點(diǎn)，形成稀疏對稱的系數(shù)矩陣；但基于局部弱式的無單元法利用了單位分解積分將全局域積分轉(zhuǎn)化為節(jié)點(diǎn)局部域積分，擺脫了背景單元的限制，是一種真正的無單元法，具有更高的靈活性和更強(qiáng)的適應(yīng)性。

本文基于戴前偉等［41］、Feng 等［42］、麻昌英等［50］、Ma 等［52］提出的探地雷達(dá)和直流電阻率EFGM 正演方法，將單位分解積分應(yīng)用于直流電阻率EFGM 正演，提出直流電阻率無單元局部弱式法（Local Weak Form Element-Free Method，LWF-EFM），實(shí)現(xiàn)了高靈活性、高精度的2.5 維直流電阻率正演模擬。將上述方法分別應(yīng)用于層狀模型及二維異常體模型的正演模擬，驗(yàn)證了本文算法的有效性，結(jié)果表明本文方法相比于FEM 和無單元Galerkin 法具有更強(qiáng)的靈活性和更高的模擬精度，可為高精度直流電阻率資料解釋提供技術(shù)支持。

1 直流電阻率法2.5 維邊值問題的全局弱式形式

當(dāng)采用第三類邊界條件時(shí)，2.5 維直流電阻率法波數(shù)域總電位U(x，λ)的邊值問題及變分問題可分別表示為［63］

式中：σ表示電導(dǎo)率；λ表示波數(shù)；I0表示電流；δ(A)表示狄拉克函數(shù)，其中A是場源點(diǎn)；K0、K1分別表示第二類零階、一階修正貝塞爾函數(shù)；Ω表示問題域（全局域）；Γs表示地表邊界，?！逓榻?cái)噙吔?；rA為場源點(diǎn)A到截?cái)噙吔绲木嚯x；n為邊界上的單位外法向；cos(rA，n)表示n與矢徑rA構(gòu)成的夾角余弦。式（2）經(jīng)EFGM 推導(dǎo)可獲得基于全局弱式的三項(xiàng)積分式K(1)、K(2)和F［50-52］

式中ΦT=[?1，?2，…，?n]為形函數(shù)向量，其中n表示支持域Ω中的場節(jié)點(diǎn)個(gè)數(shù)。不同正演方法采用不同的形函數(shù)。在無單元法中可采用徑向基點(diǎn)插值法（Radial Point Interpolation Method，RPIM）形函數(shù)［50-51］和移動(dòng)最小二乘（Moving Least Squares，MLS）形函數(shù)［52］；FEM 則采用基于單元構(gòu)造的形函數(shù)，這些形函數(shù)均具有單位分解性質(zhì)。由于RPIM 形函數(shù)具有Kronecker delta 函數(shù)性質(zhì)、數(shù)值穩(wěn)定性較好、適應(yīng)性強(qiáng)等優(yōu)點(diǎn)，本文采用RPIM 形函數(shù)進(jìn)行正演計(jì)算。采用單位分解積分法可將式（3）～式（5）中的全局域積分轉(zhuǎn)化到局部積分域求解。

2 單位分解積分

覆蓋和單位分解是單位分解積分［25，56］的基礎(chǔ)。設(shè)Ω是定義在Rd(d=1，2，3)上的有界開域，表示Ω的閉包，QN表示Ω中N個(gè)點(diǎn)的集合

基于QN定義有限開覆蓋，表示N個(gè)中點(diǎn)在xk(k=1，2，…，N)的圓域（或矩形域），滿足

則稱ON是Ω的一個(gè)有限覆蓋。式（7）中rk為點(diǎn)xk(k=1，2，…，N)對應(yīng)的覆蓋Ωk的半徑。如圖1 所示，在有界區(qū)域Ω內(nèi)由有限個(gè)節(jié)點(diǎn)xk(k=1，2，…，N)對應(yīng)的構(gòu)成了Ω上的一個(gè)有限覆蓋。

圖1 有限覆蓋示意圖

由于式（12）具有單位分解的性質(zhì)，因此稱這種積分方式為單位分解積分（Partition of Unity Quadrature，PUQ）。式（12）是單位分解積分的基礎(chǔ)，通過該式可將有界區(qū)域Ω（包含邊界）上的積分轉(zhuǎn)化成節(jié)點(diǎn)子域Ωk上的積分之和。若將單位分解積分應(yīng)用到式（3）～式（5），則積分式將用到兩個(gè)函數(shù)集：形函數(shù)Φ(x)和單位分解函數(shù)Ψ(x)。單位分解函數(shù)和形函數(shù)在本質(zhì)上是兩個(gè)不同的函數(shù)系，兩者相互獨(dú)立，任意具有單位分解性質(zhì)的函數(shù)，均滿足式（9）～式（11），因此都可作為單位分解函數(shù)。由單位分解積分定理可知，單位分解函數(shù)不會(huì)對積分計(jì)算精度產(chǎn)生影響。由于形函數(shù)（如RPIM、MLS、FEM 形函數(shù)等）均具有單位分解性質(zhì)，在求解域Ω內(nèi)任意一點(diǎn)的形函數(shù)之和為1，即

即形函數(shù)Φ(x)滿足式（9）～式（11），因此本文將正演方法中構(gòu)造的形函數(shù)Φ(x)作為單位分解函數(shù)，即

其中?k(x)為計(jì)算點(diǎn)對節(jié)點(diǎn)k的形函數(shù)。

3 直流電阻率法2.5 維LWF-EFM正演

3.1 節(jié)點(diǎn)局部積分域構(gòu)造

在LWF-EFM 中，首先需要為每一個(gè)節(jié)點(diǎn)構(gòu)造節(jié)點(diǎn)局部積分域，目的是將全局域積分轉(zhuǎn)化為節(jié)點(diǎn)局部域積分；然后，計(jì)算域Ω中的集合QN={x1，x2，…，xN}，xk∈Ω，以節(jié)點(diǎn)xk(k=1，2，…，N)為中心構(gòu)造計(jì)算域Ω的一個(gè)有限覆蓋，覆蓋Ωk的形式可以是圓形、橢圓形、矩形等。由式（12）可知，通過單位分解積分，可將全局域Ω積分轉(zhuǎn)化為節(jié)點(diǎn)局部域Ωk(k=1，2，…，N)的積分之和，其中節(jié)點(diǎn)局部積分域即為節(jié)點(diǎn)的覆蓋，因此覆蓋構(gòu)造即為節(jié)點(diǎn)局部積分域構(gòu)造。為了方便布置高斯點(diǎn)進(jìn)行積分計(jì)算，本文采用矩形節(jié)點(diǎn)局部積分域。對于矩形節(jié)點(diǎn)局部積分域，在二維情況下通常采用兩個(gè)方向的尺寸參數(shù)確定范圍。下面以直流電阻率數(shù)據(jù)為例，說明LWF-EFM 矩形節(jié)點(diǎn)局部積分域的確定方法。

對于一個(gè)節(jié)點(diǎn)，其節(jié)點(diǎn)局部積分域的尺寸ds為

式中：αs為積分域的無量綱尺寸，一般選1.0～3.0，可通過數(shù)值試驗(yàn)選取；p為節(jié)點(diǎn)局部積分域內(nèi)節(jié)點(diǎn)間的平均距離，當(dāng)節(jié)點(diǎn)均勻分布時(shí)，p取相鄰兩節(jié)點(diǎn)間的距離，節(jié)點(diǎn)非均勻分布時(shí)，p取節(jié)點(diǎn)積分域內(nèi)結(jié)點(diǎn)平均距離。一維情況下，定義平均節(jié)點(diǎn)距離為

式中：Ds為ds的預(yù)估值；nDs為積分域內(nèi)預(yù)估節(jié)點(diǎn)數(shù)。二維情況下，節(jié)點(diǎn)平均距離定義為

式中：S為積分域的預(yù)估面積；nS為積分域內(nèi)預(yù)估節(jié)點(diǎn)數(shù)。二維情況下，ds可使用兩個(gè)坐標(biāo)方向表示

式中：ds=[dsxdsz]；αsp=[αsxpxαszpz]。矩形節(jié)點(diǎn)局部積分域如圖2所示。

圖2 矩形節(jié)點(diǎn)局部積分域示意圖

圖2 中節(jié)點(diǎn)k的坐標(biāo)為xk=(xk，zk)，則矩形節(jié)點(diǎn)局部積分域的四個(gè)頂點(diǎn)坐標(biāo)為

如果該矩形節(jié)點(diǎn)局部積分域超出了全局計(jì)算域，則應(yīng)將矩形節(jié)點(diǎn)局部積分域劃在計(jì)算域內(nèi)，超出計(jì)算域Ω的部分應(yīng)去掉。當(dāng)域Ω的邊界為矩形時(shí)，邊界上矩形節(jié)點(diǎn)局部積分域頂點(diǎn)坐標(biāo)為

其中：x1，x2＜xmin；x3，x4＞xmax；z1，z4＜zmin；z2，z3＞zmax。

為保證數(shù)值積分的計(jì)算精度，通常將節(jié)點(diǎn)局部積分域進(jìn)一步劃分成ndx×ndz個(gè)子域，然后對每個(gè)子域單獨(dú)布置高斯積分點(diǎn)，計(jì)算積分并求和。例如：當(dāng)ndx=ndz=2 時(shí)，采用4 個(gè)子域，節(jié)點(diǎn)局部積分域細(xì)分方案如圖3所示。

圖3 ndx=ndz=2 時(shí)節(jié)點(diǎn)局部積分域再劃分子域示意圖

第i（i=1，2，…，ndx×ndz）個(gè)子域的四個(gè)頂點(diǎn)坐標(biāo)分別為

式中：dnx=(x4-x1)/ndx；dnz=(z2-z1)/ndz。

由于采用節(jié)點(diǎn)局部積分域計(jì)算，在地形起伏區(qū)域內(nèi)地表上的節(jié)點(diǎn)積分域與地形相交，靠近地表的節(jié)點(diǎn)局部積分域也可能與地形相交，對地形以外的區(qū)域不需要進(jìn)行積分計(jì)算，特別是地形起伏形態(tài)復(fù)雜時(shí)，需要對該部分節(jié)點(diǎn)局部積分域進(jìn)行特殊處理，使得計(jì)算區(qū)域更加貼合地形。為了更好地模擬計(jì)算域Ω的地表邊界，當(dāng)積分域落在計(jì)算域Ω地表邊界之外時(shí)，可省略Ω之外的積分域，僅對Ω內(nèi)和邊界進(jìn)行積分。根據(jù)上述節(jié)點(diǎn)局部積分域構(gòu)造方法，可知節(jié)點(diǎn)局部積分域的大小與節(jié)點(diǎn)疏密程度成正相關(guān)，即在節(jié)點(diǎn)密集分布區(qū)域?qū)?jié)點(diǎn)局部積分域自適應(yīng)縮小尺寸，反之則放大尺寸，可自動(dòng)獲得與節(jié)點(diǎn)分布相適應(yīng)的一組節(jié)點(diǎn)局部積分域。因此，節(jié)點(diǎn)局部積分域越小，越有利于模擬地形形態(tài)，因此對節(jié)點(diǎn)局部積分域進(jìn)一步細(xì)分有利于LWF-EFM 對地形起伏模型的精確模擬。

3.2 局部弱式無單元法正演

將單位分解積分式（12）分別代入式（3）和式（4），得到

應(yīng)用單位分解積分法對原積分進(jìn)行離散化后，將計(jì)算域Ω（包含邊界Γ∞）上的全局域積分轉(zhuǎn)化成節(jié)點(diǎn)局部積分域Ωk（包含邊界）上的積分之和，即將全局弱式轉(zhuǎn)化為局部弱式，完成了直流電阻率LWFEFM 的主要積分部分（圖4）。為方便起見，本文采用圓形支持域［50-52］，即給定積分點(diǎn)（即高斯點(diǎn)）支持域內(nèi)節(jié)點(diǎn)數(shù)nq后，按距離選取積分點(diǎn)最近的nq個(gè)節(jié)點(diǎn)作為該積分點(diǎn)的支持域。

圖4 基于單位分解積分的無單元法求解示意圖

假設(shè)節(jié)點(diǎn)局部積分域Ωk中采用了Ng個(gè)高斯點(diǎn)xg=(xg，zg)，g=1，2，…，Ng進(jìn)行積分計(jì)算，高斯點(diǎn)xg對應(yīng)的權(quán)重和雅可比值分別為wg和Jg，高斯點(diǎn)xg的支持域Ωq內(nèi)包含nq個(gè)節(jié)點(diǎn)，則有

展開式（23）～式（25），得到

其中?i(i=1，2，…，nq)為高斯點(diǎn)xg支持域內(nèi)nq個(gè)節(jié)點(diǎn)對應(yīng)的形函數(shù)。對于單個(gè)節(jié)點(diǎn)xk和單個(gè)高斯點(diǎn)xg，有

上式為積分點(diǎn)子系數(shù)矩陣，其矩陣元素為

式中i，j=1，2，…，nq，表示支持域Ωq中的節(jié)點(diǎn)局部編號。

對于式（5），若采用N個(gè)節(jié)點(diǎn)對計(jì)算域Ω進(jìn)行離散，可得

當(dāng)無單元法形函數(shù)具有Kronecker delta 函數(shù)性質(zhì)時(shí)，上式右端項(xiàng)為

當(dāng)無單元法形函數(shù)不具有Kronecker delta 函數(shù)性質(zhì)時(shí)，式（33）的右端項(xiàng)為

如圖5 所示，當(dāng)節(jié)點(diǎn)局部積分域相無不重疊且恰好覆蓋計(jì)算域時(shí)，節(jié)點(diǎn)局部積分域內(nèi)任意一點(diǎn)對應(yīng)該節(jié)點(diǎn)的單位分解函數(shù)為1，這種情況即為基于全局弱式積分的數(shù)值模擬法，例如EFGM 和FEM。在EFGM 中節(jié)點(diǎn)局部積分域可視為背景積分單元，在FEM 中節(jié)點(diǎn)局部積分域可視為基于節(jié)點(diǎn)拓?fù)潢P(guān)系剖分的網(wǎng)格單元。相比于基于全局弱式的EFGM，LWF-EFM 在求解K(1)、K(2)兩項(xiàng)積分式時(shí)不需要在全局域內(nèi)剖分背景積分單元，是真正意義上的無單元法。

圖5 節(jié)點(diǎn)局部積分域相互不重疊且恰好覆蓋計(jì)算域示意圖

4 模型算例

4.1 層狀模型

采用具有解析解的三層電阻率模型（圖6）驗(yàn)證直流電阻率LWF-EFM 算法的有效性。建立水平方向?qū)?80 m（x：-190～190 m），垂直方向深108 m（z：0～108 m）的矩形計(jì)算域。在地表埋設(shè)52 根電極作為供電和觀測電極，在地表x=0處進(jìn)行對稱四極裝置測深觀測。節(jié)點(diǎn)分布如圖7 所示，對計(jì)算域采用不規(guī)則節(jié)點(diǎn)進(jìn)行離散，靠近電極和地層界面區(qū)域節(jié)點(diǎn)密度較大，其他區(qū)域則采用稀疏的節(jié)點(diǎn)分布，總節(jié)點(diǎn)數(shù)為14280。分別采用LWF-EFM、EFGM 和FEM 對22組不同收發(fā)距（，其中表示電極A、B 間的距離）觀測數(shù)據(jù)進(jìn)行正演模擬。對于LWF-EFM，節(jié)點(diǎn)局部積分域以節(jié)點(diǎn)為中心分為4 個(gè)子積分域，即ndx=ndz=2，在子積分域中采用16 個(gè)高斯點(diǎn)，即G=16，支持域內(nèi)使用4個(gè)節(jié)點(diǎn)。對于EFGM 模擬，背景單元內(nèi)采用16個(gè)高斯點(diǎn)，即G=16，支持域內(nèi)使用4個(gè)節(jié)點(diǎn)。

圖6 層狀模型及對稱四極裝置觀測示意圖

圖7 層狀模型節(jié)點(diǎn)分布示意圖

圖8 和圖9 分別是ρ2為150、50 Ω·m 時(shí)的視電阻率和相對誤差曲線?？梢钥闯觯簾o論中間層為高阻或低阻，LWF-EFM 的模擬結(jié)果與解析解均吻合較好，相對誤差較小。表1 為模擬結(jié)果平均相對誤差，可見LWF-EFM 模擬結(jié)果平均相對誤差與EFGM 相差不大，均小于1%，低于FEM 平均相對誤差，表明了LWF-EFM 的有效性。層狀模型不同方法模擬結(jié)果分析表明，在計(jì)算域較小時(shí)，F(xiàn)EM 受到邊界影響較大，EFGM 和LWF-EFM 可獲得高精度的模擬結(jié)果。采用合適的模擬參數(shù)時(shí)，相比于FEM，LWF-EFM 可獲得更高的模擬精度，與EFGM 模擬精度相當(dāng)。

表1 層狀模型不同方法視電阻率平均相對誤差

圖8 層狀模型ρ2=150 Ω·m 時(shí)不同方法模擬視電阻率（左）和相對誤差（右）曲線

圖9 層狀模型ρ2=50 Ω·m 時(shí)不同方法模擬視電阻率（左）和相對誤差（右）曲線

4.2 二維矩形異常體模型

建立一個(gè)均勻半空間模型，介質(zhì)電阻率ρ1=100 Ω·m；在均勻半空間中有一個(gè)電阻率ρ2=1000 Ω·m的二維矩形異常體，異常體寬10 m（x方向），高8 m（z方向），中心點(diǎn)坐標(biāo)（x，z）為（0，15 m）［52］。計(jì)算域Ω水平方向?qū)?00 m（x：-100～100 m），垂直方向高100 m（z：0～100 m）。

分別采用三種類型節(jié)點(diǎn)分布離散計(jì)算域：第1種為均勻節(jié)點(diǎn)分布，節(jié)點(diǎn)間距為1 m，節(jié)點(diǎn)總數(shù)為20301；第2 種節(jié)點(diǎn)分布是在第1 種節(jié)點(diǎn)分布基礎(chǔ)上，對靠近場源的淺層區(qū)域和異常體區(qū)域進(jìn)行節(jié)點(diǎn)加密，節(jié)點(diǎn)總數(shù)為25150（圖10a）；第3 種節(jié)點(diǎn)分布是在第2 種節(jié)點(diǎn)分布基礎(chǔ)上，對計(jì)算域外圍區(qū)域的節(jié)點(diǎn)進(jìn)行抽稀，節(jié)點(diǎn)總數(shù)為12878（圖10b）。在地表水平方向-58～58 m 范圍內(nèi)以2 m 的間隔均勻布置59 根電極，進(jìn)行直流電阻率法溫納裝置［64］觀測。下文若不加以說明，均采用該電極布置和觀測裝置。為減小邊界影響，采用FEM 進(jìn)行正演模擬時(shí)，首先在第2 種節(jié)點(diǎn)分布基礎(chǔ)上對計(jì)算域進(jìn)行擴(kuò)邊，擴(kuò)邊后的范圍為2000 m×1000 m；然后，分別使用第1種節(jié)點(diǎn)分布和第3 種節(jié)點(diǎn)分布，采用EFGM 和LWFEFM 進(jìn)行正演模擬。EFGM 模擬采用1 m×1 m 的背景積分單元。這兩種方法計(jì)算過程中，每個(gè)積分單元（或節(jié)點(diǎn)子局部積分域）均采用16 個(gè)高斯積分點(diǎn)。

圖10 二維模型模擬采用的剖分節(jié)點(diǎn)分布示意圖［52］

圖11為基于第2種節(jié)點(diǎn)分布、采用FEM 計(jì)算的電阻率擬斷面。由于第2種節(jié)點(diǎn)分布相比于其他兩種節(jié)點(diǎn)分布節(jié)點(diǎn)數(shù)最多，且在靠近場源和異常體周圍區(qū)域進(jìn)行了節(jié)點(diǎn)加密，因此該節(jié)點(diǎn)分布條件下的計(jì)算結(jié)果更精確，故將圖11所示結(jié)果作為標(biāo)準(zhǔn)結(jié)果。圖12a和圖12b 為采用EFGM 獲得的模擬結(jié)果，圖12c 和圖12d 為EFGM 模擬結(jié)果與參考標(biāo)準(zhǔn)結(jié)果的相對誤差。對比圖12c和圖12d可知，使用第3種節(jié)點(diǎn)分布的相對誤差較小，大部分區(qū)域相對誤差接近0，而使用第1 種節(jié)點(diǎn)分布時(shí)，在淺層和異常體周圍區(qū)域相對誤差明顯相對較大。

圖11 二維模型基于第2 種節(jié)點(diǎn)的FEM 視電阻率擬斷面

圖12 二維模型EFGM 視電阻率擬斷面（上）及相對誤差剖面（下）

圖13為采用LWF-EFM 獲得的模擬結(jié)果。對比圖13a、圖13b與圖11可知，LWF-EFM 的模擬結(jié)果與參考結(jié)果基本一致，表明LWF-EFM 模擬結(jié)果是可靠的。對比圖13c和圖13d 可知，與EFGM 類似，LWFEFM 使用第3 種節(jié)點(diǎn)分布得到的模擬結(jié)果相對誤差較小，大部分區(qū)域相對誤差接近0，而采用第1 種節(jié)點(diǎn)分布時(shí)，在淺層和異常體周圍區(qū)域相對誤差明顯相對較大。

圖13 二維模型LWF-EFM 視電阻率擬斷面（上）及相對誤差剖面（下）

表2為第1種和第2種節(jié)點(diǎn)分布的節(jié)點(diǎn)數(shù)、模擬結(jié)果最大相對誤差、平均相對誤差和計(jì)算時(shí)間統(tǒng)計(jì)?？芍捎玫?種節(jié)點(diǎn)分布時(shí)最大相對誤差和平均相對誤差均大于第3種節(jié)點(diǎn)分布；同時(shí)，由于第1 種節(jié)點(diǎn)分布使用了更多的節(jié)點(diǎn)，計(jì)算成本顯著更高。對比相對誤差統(tǒng)計(jì)結(jié)果可知，與EFGM 一樣，LWF-EFM 也可采用不規(guī)則的節(jié)點(diǎn)分布，在場源和異常體周圍區(qū)域加密節(jié)點(diǎn)，在計(jì)算域外圍區(qū)域使用稀疏的節(jié)點(diǎn)分布，可在使用更少節(jié)點(diǎn)的條件下獲得與均勻節(jié)點(diǎn)分布相同或更高的模擬精度；同時(shí)，使用更少的節(jié)點(diǎn)也可節(jié)省計(jì)算成本。相比EFGM，LWF-EFM 計(jì)算效率較低，主要原因是不同節(jié)點(diǎn)分布的節(jié)點(diǎn)局部積分域之間通常存在很多重合區(qū)域，實(shí)際計(jì)算積分面積大于計(jì)算域Ω，導(dǎo)致計(jì)算量增加。但由于采用節(jié)點(diǎn)局部積分域進(jìn)行積分，僅需要設(shè)定節(jié)點(diǎn)局部積分域參數(shù)，不再需要背景單元剖分，相對需要背景單元的EFGM 對單元的依賴性更小，是真正的無單元法，正演過程更方便。同時(shí)，在計(jì)算域內(nèi)任意區(qū)域進(jìn)行節(jié)點(diǎn)加密或者稀釋均十分方便，可根據(jù)節(jié)點(diǎn)分布情況自動(dòng)構(gòu)造節(jié)點(diǎn)局部積分域，無需其他處理，具有更高的靈活性。

表2 二維模型LWF-EFM 不同節(jié)點(diǎn)分布的節(jié)點(diǎn)數(shù)、模擬相對誤差和計(jì)算時(shí)間

4.3 復(fù)雜地電模型

建立圖14a 所示水平地形起伏模型，地下空間介質(zhì)電阻率ρ1=100 Ω·m。計(jì)算域Ω覆蓋x方向300 m（-150～150 m），垂直方向（z）最大深度為115 m，山脊最高點(diǎn)高度為7.6 m，位于x=-14 m 處，山谷最低處位于x=10 m，高度為3.4 m［51］。

圖14 地形起伏模型節(jié)點(diǎn)示意圖

基于圖14a 所示節(jié)點(diǎn)分布，采用不規(guī)則節(jié)點(diǎn)對場源附近的淺層區(qū)域和地形起伏區(qū)域進(jìn)行節(jié)點(diǎn)加密，在計(jì)算域的外圍區(qū)域使用稀疏節(jié)點(diǎn)分布（圖14b）。分別采用LWF-EFM、EFGM 和FEM 對模型進(jìn)行正演，模擬直流電阻率法的視電阻率。采用EFGM 模擬時(shí)，在地形起伏區(qū)域采用任意四邊形背景單元模擬地形。為減小邊界影響，采用FEM 模擬時(shí)，在上述計(jì)算域的基礎(chǔ)上進(jìn)行擴(kuò)邊，擴(kuò)邊后水平方向?qū)?000 m，垂直方向高1000 m。

圖15 為分別采用地形起伏模型未加密節(jié)點(diǎn)分布（圖14a）和地表加密節(jié)點(diǎn)分布（圖14b）時(shí)，EFGM、FEM 和LWF-EFM 模擬結(jié)果。由圖可知，無論采用未加密節(jié)點(diǎn)還是加密節(jié)點(diǎn)，EFGM 和FEM 的計(jì)算結(jié)果基本一致。其中，EFGM 通過采用任意四邊形背景單元，很好地模擬了地形起伏。對比圖15c（左）與圖15a（左）、圖15b（左）可知，在局部區(qū)域LWF-EFM與EFGM 和FEM 的模擬結(jié)果有較明顯的差異，表明在地形起伏區(qū)域，若采用較稀疏的節(jié)點(diǎn)分布，LWFEFM 對地形起伏模型的模擬精度較低。其主要原因是LWF-EFM 是基于節(jié)點(diǎn)局部積分域進(jìn)行積分，在地形起伏區(qū)域節(jié)點(diǎn)局部積分域不可避免地會(huì)與地形相交，因而不能準(zhǔn)確地模擬地形，導(dǎo)致模擬結(jié)果出現(xiàn)偏差。在地形起伏區(qū)域，LWF-EFM 采用節(jié)點(diǎn)局部積分域進(jìn)行積分，理論上當(dāng)節(jié)點(diǎn)局部積分域與地形相交時(shí)，需要將位于地形以外的積分域區(qū)域刪除。相比于EFGM 和基于非結(jié)構(gòu)化三角形的FEM，LWF-EFM處理復(fù)雜地形與積分域的難度和工作量會(huì)很大，因此需要提高LWF-EFM 對地形的模擬精度，

圖15 地形起伏模型采用未加密節(jié)點(diǎn)（左）和加密節(jié)點(diǎn)（右）時(shí)不同方法正演視電阻率擬斷面

首先，由于LWF-EFM 通常根據(jù)節(jié)點(diǎn)分布構(gòu)造節(jié)點(diǎn)局部積分域，本文對地表起伏段進(jìn)行合理的節(jié)點(diǎn)加密，使得靠近地形起伏區(qū)域的節(jié)點(diǎn)局部積分域使用較小尺寸的剖分單元，可更精確地描述地形起伏形態(tài)，提高對地形的模擬精度；同時(shí)，在地表起伏區(qū)域和源附近加密節(jié)點(diǎn)，也有利于提高模擬精度。其次，當(dāng)高斯點(diǎn)位于地形以外時(shí)，本文采用對高斯點(diǎn)不做計(jì)算的處理策略。雖然這可能出現(xiàn)同一子域內(nèi)部分高斯點(diǎn)進(jìn)行了計(jì)算而另一部分高斯點(diǎn)不做計(jì)算的情況，造成一定的積分誤差，但這樣處理方法簡單方便，同時(shí)可保證積分區(qū)域盡可能地貼合地形。

圖15c（右）為LWF-EFM 經(jīng)上述兩種方法處理后的模擬結(jié)果。對比圖15c（左）與圖15c（右）可知，在地形起伏區(qū)域進(jìn)行加密節(jié)點(diǎn)可明顯改善LWFEFM 對地形起伏模型的模擬效果。對比圖15c（右）與圖15a（右）、圖15b（右）可知，若采用加密節(jié)點(diǎn)，三種算法對地形起伏模型的正演模擬結(jié)果基本一致，表明采用上述兩種方法可有效提高LWF-EFM 對地形起伏的模擬精度。此外，由于LWF-EFM 采用節(jié)點(diǎn)局部積分域積分，相對于EFGM 需要考慮背景單元分布與節(jié)點(diǎn)分布相適應(yīng)的問題［52］，LWF-EFM在節(jié)點(diǎn)加密方面更方便，設(shè)計(jì)任意節(jié)點(diǎn)分布具有更高的靈活性和適應(yīng)性。

為了模擬地下復(fù)雜形態(tài)異常體，在地形起伏模型基礎(chǔ)上加入具有復(fù)雜形態(tài)的兩個(gè)近直立異常體（圖16中藍(lán)色和紅色區(qū)域）。其中，紅色區(qū)域?yàn)楦咦璁惓ｓw，位于地形起伏模型山脊頂點(diǎn)的下方，電阻率為ρ2；藍(lán)色區(qū)域?yàn)榈妥璁惓ｓw，位于地形起伏模型山谷低點(diǎn)的下方，電阻率為ρ3。

圖16 復(fù)雜地電模型節(jié)點(diǎn)示意圖

使用本文方法分別對ρ2=200 Ω·m、ρ3=20 Ω·m和ρ2=5000 Ω·m、ρ3=1 Ω·m 兩種情況下的復(fù)雜地電模型進(jìn)行正演模擬，結(jié)果如圖17 所示。分析圖17 可知，山脊地形可以引起明顯的低阻異常，山谷地形可以引起明顯的高阻異常。對比分析圖17a和圖17b可知，在山脊和山谷下方分別加入高阻和低阻異常體后，若異常體的電阻率與圍巖的電阻率差異不大時(shí)，高阻和低阻異常體引起的異?；颈坏匦我鸬漠惓ｄ螞]；隨著異常體與圍巖之間電阻率差異的增加，異常體引起的異常逐步顯現(xiàn)，但未能完全湮沒地形的影響。

圖17 復(fù)雜地電模型LWF-EFM 模擬視電阻率擬斷面圖

因此，對復(fù)雜地電模型正演結(jié)果綜合分析表明，直流電阻率法受地形影響較大，地形之下的異常體（如直立異常體）產(chǎn)生的視電阻率異常很可能被地形引起的異常掩蓋。

5 結(jié)論

本文提出一種靈活、高精度的直流電阻率無單元局部弱式正演方法。該方法將單位分解積分應(yīng)用于基于全局弱式積分的直流電阻率EFGM 正演，將全局域積分轉(zhuǎn)化為節(jié)點(diǎn)局部積分域積分。對不同地電模型應(yīng)用三種方法（EFGM，F(xiàn)EM，LWF-EFM）進(jìn)行數(shù)值模擬，結(jié)果驗(yàn)證了本文算法應(yīng)用于直流電阻率正演模擬的正確性和有效性。與EFGM 相比，本文算法無需在全局域內(nèi)進(jìn)行背景積分單元剖分；與FEM 相比，本文算法不需要進(jìn)行基于節(jié)點(diǎn)拓?fù)潢P(guān)系的單元剖分，是一種更靈活性、適應(yīng)性較強(qiáng)的高精度無單元正演模擬算法。模擬結(jié)果表明，由于本文算法采用的不同節(jié)點(diǎn)局部積分域之間通常存在很多重合區(qū)域，導(dǎo)致實(shí)際積分面積比計(jì)算域Ω大，計(jì)算量也隨之增加。地形起伏模型模擬結(jié)果分析表明，相比于EFGM 和FEM，本文算法在地形起伏區(qū)域節(jié)點(diǎn)局部積分域不可避免地會(huì)與地形相交，不利于地形起伏模型的精確模擬，因而本文對地形起伏區(qū)域進(jìn)行節(jié)點(diǎn)加密，對地形以外的高斯點(diǎn)不做計(jì)算，明顯提高了地形起伏模型的模擬精度。