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

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

0 引言

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

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

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

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

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

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

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

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

2 單位分解積分

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

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

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

圖1 有限覆蓋示意圖

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

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

其中?k(x)為計算點對節(jié)點k的形函數。

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

3.1 節(jié)點局部積分域構造

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

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

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

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

式中：S為積分域的預估面積；nS為積分域內預估節(jié)點數。二維情況下，ds可使用兩個坐標方向表示

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

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

圖2 中節(jié)點k的坐標為xk=(xk，zk)，則矩形節(jié)點局部積分域的四個頂點坐標為

如果該矩形節(jié)點局部積分域超出了全局計算域，則應將矩形節(jié)點局部積分域劃在計算域內，超出計算域Ω的部分應去掉。當域Ω的邊界為矩形時，邊界上矩形節(jié)點局部積分域頂點坐標為

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

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

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

第i（i=1，2，…，ndx×ndz）個子域的四個頂點坐標分別為

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

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

3.2 局部弱式無單元法正演

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

應用單位分解積分法對原積分進行離散化后，將計算域Ω（包含邊界?！蓿┥系娜钟蚍e分轉化成節(jié)點局部積分域Ωk（包含邊界）上的積分之和，即將全局弱式轉化為局部弱式，完成了直流電阻率LWFEFM 的主要積分部分（圖4）。為方便起見，本文采用圓形支持域［50-52］，即給定積分點（即高斯點）支持域內節(jié)點數nq后，按距離選取積分點最近的nq個節(jié)點作為該積分點的支持域。

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

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

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

其中?i(i=1，2，…，nq)為高斯點xg支持域內nq個節(jié)點對應的形函數。對于單個節(jié)點xk和單個高斯點xg，有

上式為積分點子系數矩陣，其矩陣元素為

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

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

當無單元法形函數具有Kronecker delta 函數性質時，上式右端項為

當無單元法形函數不具有Kronecker delta 函數性質時，式（33）的右端項為

如圖5 所示，當節(jié)點局部積分域相無不重疊且恰好覆蓋計算域時，節(jié)點局部積分域內任意一點對應該節(jié)點的單位分解函數為1，這種情況即為基于全局弱式積分的數值模擬法，例如EFGM 和FEM。在EFGM 中節(jié)點局部積分域可視為背景積分單元，在FEM 中節(jié)點局部積分域可視為基于節(jié)點拓撲關系剖分的網格單元。相比于基于全局弱式的EFGM，LWF-EFM 在求解K(1)、K(2)兩項積分式時不需要在全局域內剖分背景積分單元，是真正意義上的無單元法。

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

4 模型算例

4.1 層狀模型

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

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

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

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

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

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

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

4.2 二維矩形異常體模型

建立一個均勻半空間模型，介質電阻率ρ1=100 Ω·m；在均勻半空間中有一個電阻率ρ2=1000 Ω·m的二維矩形異常體，異常體寬10 m（x方向），高8 m（z方向），中心點坐標（x，z）為（0，15 m）［52］。計算域Ω水平方向寬200 m（x：-100～100 m），垂直方向高100 m（z：0～100 m）。

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

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

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

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

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

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

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

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

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

4.3 復雜地電模型

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

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

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

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

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

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

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

為了模擬地下復雜形態(tài)異常體，在地形起伏模型基礎上加入具有復雜形態(tài)的兩個近直立異常體（圖16中藍色和紅色區(qū)域）。其中，紅色區(qū)域為高阻異常體，位于地形起伏模型山脊頂點的下方，電阻率為ρ2；藍色區(qū)域為低阻異常體，位于地形起伏模型山谷低點的下方，電阻率為ρ3。

圖16 復雜地電模型節(jié)點示意圖

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

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

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

5 結論

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