尤 淼，張金會，解建建

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

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

為了解磁電阻率法的特點和異常體的響應(yīng)曲線特征，使用基于矩形網(wǎng)格的有限單元法進(jìn)行二維磁電法數(shù)值模擬并進(jìn)行模型計算。

1 基本原理

穩(wěn)定電流場滿足的Maxwell方程組為：

存在本構(gòu)方程：

σ為電導(dǎo)率，μ為真空磁導(dǎo)率。

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

使用有限單元法計算點電流源在研究區(qū)域引起的標(biāo)量電位u，在利用上式，可計算電場強度E，進(jìn)而求出穩(wěn)定電流場激發(fā)出的磁場強度H。

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

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

圖1 二維電場示意圖

滿足如下關(guān)系：

根據(jù)奧斯特-高斯公式，有：

δ函數(shù)滿足積分關(guān)系：

因此有如下關(guān)系：

因為電位u與電流密度矢量的關(guān)系為：

σ為介質(zhì)的電導(dǎo)率，整理上式，可得電位滿足的微分方程：

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

在無窮遠(yuǎn)邊界Γ∞上，可假設(shè)研究區(qū)域的不均勻性對邊界上的電位沒有影響，其表達(dá)式為：

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

因研究區(qū)域為三維場源、二維構(gòu)造，對于此類邊值問題，需要進(jìn)行傅里葉變換，變換到二維波數(shù)域中進(jìn)行求解，經(jīng)過傅里葉變換后的二維邊值問題為：

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

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

3 網(wǎng)格剖分

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

圖2 矩形網(wǎng)格剖分示意圖

4 總體集成

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

KU=P

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

5 傅里葉反變換

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

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

采用最優(yōu)化方法，選取5點波數(shù)進(jìn)行傅里葉反變換[3]。

6 計算電流密度

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

在直角坐標(biāo)系下展開，得到：

7 計算磁場

已知穩(wěn)定電流磁場滿足如下關(guān)系式：

將上式在直角坐標(biāo)系下展開，得到電流密度矢量和磁場矢量的關(guān)系式：

將所有節(jié)點電流密度和磁場關(guān)系式關(guān)系整理，得到如下方程組：

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

8 模型驗證

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

8.1 模型一

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

圖3 模型一地表各節(jié)點磁場Hy分量

8.2 模型二

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

圖4 模型二地表各節(jié)點磁場Hy分量

8.3 模型三

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

圖5 模型三剖面圖

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

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

9 結(jié)論

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