磁電阻率法二維有限元數值模擬

尤 淼，張金會，解建建

（安徽省勘查技術院，安徽 合肥 230031)

磁電阻率法（Magnetometric resistivity method,MMR）是通過測量兩點間由人工直流電場激發的磁場的勘探方法。主要研究人工傳導電流和電化學作用引起的極化電流在空間各點產生測磁場變化規律，達到找礦目的。數值計算方面，鄧靖武[1]（2005）使用三維有限差分法研究了磁電法正演理論；楊從濤[2]（2010）使用二維有限單元法進行了磁電法的數值計算。

為了解磁電阻率法的特點和異常體的響應曲線特征，使用基于矩形網格的有限單元法進行二維磁電法數值模擬并進行模型計算。

1 基本原理

穩定電流場滿足的Maxwell方程組為：

存在本構方程：

σ為電導率，μ為真空磁導率。

由Maxwell方程組和本構方程可以看出，要獲得磁場強度H，需首先計算電場強度E。由于電場E是矢量場，引入標量電位u，令電場E表示為標量電位u的梯度形式：

使用有限單元法計算點電流源在研究區域引起的標量電位u，在利用上式，可計算電場強度E，進而求出穩定電流場激發出的磁場強度H。

2 點源二維地電斷面的正演問題

假設在地表A處，存在一電流源I，電流密度矢量為 。研究區域Ω為二維地電斷面，與構造走向垂直，邊界為Γ，如圖1。

圖1 二維電場示意圖

滿足如下關系：

根據奧斯特-高斯公式，有：

δ函數滿足積分關系：

因此有如下關系：

因為電位u與電流密度矢量的關系為：

σ為介質的電導率，整理上式，可得電位滿足的微分方程：

考慮地面Γs上的電流，因空氣中電阻率可視為無窮大，所以流向空氣中的電流密度為0，即電流密度在地面上的法向分量為0，因此電位的法向導數為0：

在無窮遠邊界Γ∞上，可假設研究區域的不均勻性對邊界上的電位沒有影響，其表達式為：

其中，c是比例系數，r"是點源到該邊界點的距離。有限元求解區域需足夠大，以滿足此邊界條件。

因研究區域為三維場源、二維構造，對于此類邊值問題，需要進行傅里葉變換，變換到二維波數域中進行求解，經過傅里葉變換后的二維邊值問題為：

其中，σ為電導率，k為波數域中的參數，U是波數域中的電位，I是供電電流，δ(A)是供電點的位置，為波數域中電位對地表邊界的外法向。K0為第二類零階修正貝塞爾函數，K1為第二類一階修正貝塞爾函數，cos(r,n)為邊界上A點到該邊界點的矢徑r與該點外法向n之間夾角的余弦。點源二維電場的內邊界條件是自然邊界條件，在泛函極值過程中將自動滿足。

使用加權余量法，得到點源二維電場的變分問題為：

3 網格剖分

研究區域采用如下圖所示的矩形網格剖分，有限單元法采用矩形單元，雙線性插值，在區域中心進行網格加密處理，邊界處網格稀疏化。

圖2 矩形網格剖分示意圖

4 總體集成

將區域積分分解為各矩形單元上的積分，根據疊加原理，擴展成全體節點組成的線性方程組：

KU=P

有限單元法計算中，剛度矩陣K采用變帶寬存儲，以節約計算機內存空間。使用Cholesky分解法求解上式的方程組，獲得波數域k中各節點的電位值U。

5 傅里葉反變換

對于每個節點，由不同波數k計算出一組U，再使用傅立葉反變換計算出三維空間的電位u。

由于u(x,y,-z)=u(x,y,z)，故采用余弦形式，數值積分方法：

采用最優化方法，選取5點波數進行傅里葉反變換[3]。

6 計算電流密度

由有限單元法計算獲得各節點電位U后，根據下式計算各節點電流密度j：

在直角坐標系下展開，得到：

7 計算磁場

已知穩定電流磁場滿足如下關系式：

將上式在直角坐標系下展開，得到電流密度矢量和磁場矢量的關系式：

將所有節點電流密度和磁場關系式關系整理，得到如下方程組：

其中，jx1表示1號節點在x方向的電流密度分量值，上式中，節點總數為M×N。求解上述方程組，即可得到二維研究區域內，各節點的磁場分量。

8 模型驗證

網格剖分在水平方向為中間區域網格數40，水平間距1m；左右區域各12個稀疏網格，成等比例剖分。剖分區域水平方向總長度為200m。垂直方向25個網格，淺部加密剖分，深部稀疏剖分。最大深度為800m。

8.1 模型一

模型一為均勻半空間模型，電阻率為10Ω·m，水平偶極源AB，偶極中心深度為6m，偶極距4m。地表各節點磁場Hy分量如下圖所示。

圖3 模型一地表各節點磁場Hy分量

8.2 模型二

模型二為均勻半空間模型，電阻率為10Ω·m，垂直偶極源AB，偶極中心深度為6m，偶極距10m。地表各節點磁場Hy分量如下圖所示。

圖4 模型二地表各節點磁場Hy分量

8.3 模型三

模型三為二維模型，地下背景電阻率為10Ω·m，水平偶極源AB，偶極中心深度為6m，偶極距4m。在地表測線中心向下深度30m處，存在一個規模4m×4m的高阻異常，電阻率為100Ω·m。

圖5 模型三剖面圖

計算獲得的地表各節點磁場分量如下圖所示。

圖6 模型三地表節點磁場Hy分量曲線圖

9 結論

本文使用二維矩形網格剖分的有限單元法，計算了幾個典型地電模型的磁電法響應。因采用矩形網格剖分線性插值，數值解的結果曲線不夠平滑，需使用三角形等更精細的剖分方式或高階插值方式。二維計算中，磁場所求解方程組為欠定，使用最小二乘法求解，會造成精度的損失，需考慮使用三維點源電位進行計算。