三維廣域電磁法自適應(yīng)有限元正演

周峰，徐錦通，張志勇，藍(lán)澤鸞，陳輝，湯洪志

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

0 引言

廣域電磁法是由何繼善院士及其團(tuán)隊(duì)潛心研究多年提出的一種頻率域有源電磁方法［1-4］。相比可控源音頻大地電磁法（CSAMT），廣域電磁法僅需單分量電磁場即可定義視參數(shù)，而無需正交電場和磁場，野外工作效率較高［5-8］。另外，由于CSAMT 沿用了卡尼亞視電阻率計(jì)算公式，需要利用正交電場和磁場計(jì)算阻抗，導(dǎo)致過渡區(qū)和近區(qū)的數(shù)據(jù)不滿足平面波的假設(shè)，因而這部分?jǐn)?shù)據(jù)不能參與資料解譯工作，這在一定程度上影響了CSAMT 勘探效果。相反，廣域電磁法摒棄了遠(yuǎn)區(qū)的限制，擴(kuò)展了頻率域可控場源觀測區(qū)范圍，具有勘探深度大、觀測范圍廣、工作效率高等優(yōu)點(diǎn)，對地下地質(zhì)結(jié)構(gòu)的探測效果明顯。近年來，廣域電磁法已經(jīng)在油氣、頁巖氣等礦產(chǎn)資源、地下水資源以及環(huán)境工程等領(lǐng)域得到了廣泛應(yīng)用，取得了較好的應(yīng)用效果［1，9-10］。

現(xiàn)階段廣域電磁法資料解譯工作主要集中在一維、二維［6，11］，迫切需要開展三維廣域電磁法資料解譯技術(shù)研究。三維正演模擬是資料解譯的重要環(huán)節(jié)，同時(shí)又是認(rèn)識(shí)復(fù)雜地質(zhì)結(jié)構(gòu)響應(yīng)特征的必要途徑。因此，研究高效率、高精度的三維廣域電磁法正演方法成為目前研究的熱點(diǎn)。

三維數(shù)值模擬方法主要包含積分方程法［12-16］、有限差分法［17-19］、有限體積法［20-22］、譜元法［23-24］以及有限單元法［25-29］。每一種數(shù)值模擬方法都具有不同的應(yīng)用特點(diǎn)，如：積分方程（IE）法是一種純二次場算法，只需要對異常體進(jìn)行剖分，具有較高的計(jì)算精度，但該方法形成的系統(tǒng)矩陣是密實(shí)矩陣，求解效率低，需要開展相應(yīng)的加速算法（快速FFT）［14，30］；有限差分法數(shù)值模擬方法理論簡單，網(wǎng)格拓?fù)潢P(guān)系簡潔［17］，但刻畫復(fù)雜地形及模型較困難；有限單元法［31-32］是目前應(yīng)用最廣的數(shù)值模擬方法之一，隨著非結(jié)構(gòu)網(wǎng)格離散技術(shù)的發(fā)展，對起伏地形和復(fù)雜地質(zhì)結(jié)構(gòu)模型的剖分更加精細(xì)［33］。隨著自適應(yīng)技術(shù)的發(fā)展，基于后驗(yàn)誤差估計(jì)策略的自適應(yīng)有限元技術(shù)得到飛速發(fā)展［34-35］。自適應(yīng)加密技術(shù)主要采用兩類策略：一類基于全局誤差估計(jì)的自適應(yīng)加密策略；另一類是面向目標(biāo)的自適應(yīng)加密策略。前者主要關(guān)注整體數(shù)值誤差；后者重點(diǎn)考慮的則是觀測點(diǎn)的數(shù)值誤差，對非重點(diǎn)區(qū)域不做過多的考慮，避免了不必要的網(wǎng)格離散。當(dāng)前，自適應(yīng)有限元正演技術(shù)已在頻率域電磁法數(shù)值模擬中得到較廣泛的應(yīng)用［36-38］。但是，關(guān)于廣域電磁法的自適應(yīng)正演算法鮮有報(bào)道。此外，基于有限體積法的電磁法數(shù)值模擬近年來得到了更多關(guān)注，尤其是與非結(jié)構(gòu)網(wǎng)格相結(jié)合的應(yīng)用發(fā)展較快。譜元法是有限元和譜方法相結(jié)合而發(fā)展起來的一種數(shù)值模擬方法，兼顧了譜方法的高精度和有限元的靈活性［23-24］，近年來在電磁法數(shù)值模擬方面得到了廣泛關(guān)注。目前，該方法網(wǎng)格離散大多采用結(jié)構(gòu)化網(wǎng)格，而關(guān)于非結(jié)構(gòu)化的譜方法數(shù)值模擬技術(shù)報(bào)道較少［39-40］。

綜上所述，為了更好地模擬復(fù)雜地電模型及提高廣域電磁法正演精度，本文采用自適應(yīng)有限元技術(shù)并結(jié)合非結(jié)構(gòu)四面體單元實(shí)現(xiàn)三維廣域電磁法的高精度、高效率正演計(jì)算。首先，從廣域電磁法滿足的微分方程出發(fā)，推導(dǎo)出三維廣域電磁法滿足的電場方程，并采用矢量有限元對電場方程進(jìn)行離散，形成線性方程組；其次，開展以電流密度法向連續(xù)后驗(yàn)誤差估計(jì)策略為依據(jù)的自適應(yīng)有限元算法進(jìn)行高精度電磁場分量計(jì)算，并采用一種穩(wěn)健的Brent-Dekker算法計(jì)算廣域視電阻率；最后，設(shè)計(jì)地電模型驗(yàn)證自適應(yīng)算法的正確性，并對廣域視電阻率響應(yīng)特征進(jìn)行分析。

1 方法理論

1.1 廣域電磁法邊值問題

假定負(fù)諧時(shí)間因子為e-iωt，電場E和磁場H滿足Maxwell方程

式中：σ表示電導(dǎo)率；μ表示磁導(dǎo)率；σ-iωε為復(fù)電導(dǎo)率，其中ω為角頻率，ε為介電常數(shù)；Js為外部電流密度。對式（1）兩邊取旋度，可得雙旋度結(jié)構(gòu)的電場方程

以式（5）為控制方程，利用Galerkin 有限元推導(dǎo)有限元方程。設(shè)余量

令r在每個(gè)計(jì)算單元Ωe內(nèi)滿足

式中Ne表示矢量基函數(shù)。將式（6）代入式（7），經(jīng)矢量恒等變換，式（7）可改寫為

采用非結(jié)構(gòu)四面體單元離散整個(gè)求解區(qū)域，引入矢量基函數(shù)N，則四面體單元內(nèi)任意一點(diǎn)的電場可由各條棱邊的電場和矢量基函數(shù)表示

式中：K=C+M為剛度矩陣，其中；E為未知電場；為外部場源積分項(xiàng)。本文將場源置于四面體單元的各邊［41］，假定源的中點(diǎn)坐標(biāo)為(x0，y0，z0)，外部電流密度Js可表示為

式中：I為電流密度；δ為Delta函數(shù)。為了對長導(dǎo)線源進(jìn)行積分，可將有限長源看作多個(gè)偶極源的疊加。此外，本文采用擴(kuò)邊技術(shù)加載截?cái)噙吔鐥l件，即E|Γ=0，這里Γ為無窮遠(yuǎn)邊界。

1.2 基于后驗(yàn)誤差估計(jì)的自適應(yīng)正演策略

為保證有限元數(shù)值解的精度，需人為對網(wǎng)格單元進(jìn)行優(yōu)化，這一過程耗時(shí)且費(fèi)力。近年來，基于后驗(yàn)誤差估計(jì)的自適應(yīng)有限元技術(shù)逐步應(yīng)用于地球物理數(shù)值模擬，提高了正演計(jì)算的精度。為了實(shí)現(xiàn)三維多場源廣域電磁法自適應(yīng)正演模擬，需要構(gòu)建三維廣域電磁法問題對應(yīng)的對偶問題

式中：T表示測點(diǎn)所涉及四面體單元個(gè)數(shù)；Vk為第k個(gè)單元的體積；I0=(1，1，1)為三分量的單位源；W表示電場對應(yīng)的對偶場；κ表示權(quán)重系數(shù)，取值為觀測點(diǎn)到場源中心之間的距離的立方，其作用是平衡場源對每個(gè)觀測點(diǎn)的影響［33］。式（13）與式（5）具有相同的結(jié)構(gòu)形式，采用Galerkin矢量有限元對式（13）進(jìn)行離散，可得到相同的剛度矩陣K，即對偶問題形成的線性方程組為

后驗(yàn)誤差分析是自適應(yīng)網(wǎng)格加密的核心引擎［33-35］。為此，本文采用以電流密度法向連續(xù)為后驗(yàn)估計(jì)的網(wǎng)格自適應(yīng)細(xì)化方案［33，36-37］，單元的后驗(yàn)誤差e可根據(jù)下式求得

式中：?Δi表示四面體的第i個(gè)面；J=σ·E表示電流密度，J1、J2分別表示共面?Δi的面電流密度；n表示面法向。因此，本文利用式（15）評(píng)價(jià)單元的數(shù)值誤差、標(biāo)記需要細(xì)化的四面體單元。同時(shí)，為了提高測點(diǎn)計(jì)算精度，對對偶問題線性方程（式（14）） 進(jìn)行求解，對偶問題所對應(yīng)的單元誤差w可根據(jù)下式求得

式中：W1、W2分別表示共面?Δi的對偶場；σ1、σ2分別表示相鄰四面體單元的電導(dǎo)率。然后，采用加權(quán)思想，最終獲得加權(quán)誤差估計(jì)為

通過上述過程可以標(biāo)記每個(gè)單元的加權(quán)誤差，然后按照一定比例對四面體進(jìn)行細(xì)化處理，從而得到新的網(wǎng)格單元；然后，重復(fù)上述過程，直到有限元解的精度滿足要求或達(dá)到最大網(wǎng)格細(xì)化迭代次數(shù)為止；最后，通過得到的電場或磁場分量計(jì)算廣域視電阻率。需要指出，可控源電磁法在場源位置不滿足電流密度法向連續(xù)的條件，因而包含場源的單位誤差會(huì)很大。為此，本文利用局部加密技術(shù)細(xì)化場源，在計(jì)算單元誤差時(shí)將源所在單元剔除，以降低源的影響，保證細(xì)化過程朝預(yù)設(shè)目標(biāo)進(jìn)行。

1.3 廣域視電阻率的計(jì)算

廣域電磁法的特點(diǎn)是采用電磁場單分量進(jìn)行計(jì)算，利用均勻半空間條件下電磁場分量關(guān)于電阻率的非線性方程定義廣域電磁視電阻率，克服了卡尼亞電阻率定義中對遠(yuǎn)區(qū)的要求，拓展了人工場源的觀測范圍，簡化了野外觀測裝備，可顯著提高工作效率。為了闡述廣域視電阻率的優(yōu)勢，以x方向布設(shè)的電偶源為例，均勻半空間的視電阻率為

式中：r表示收發(fā)距；k表示波數(shù)；I表示電流強(qiáng)度；Ex是沿觀測方向x兩點(diǎn)M、N 間的電場；為M、N 間的電位差，其中表示觀測電極M、N 之間的距離；，其中dL表示偶極源的長度；FE-Ex(ikr)=1-3 sin2?+e-ikr(1+ikr)，其中?表示測點(diǎn)和偶極源中心連線與偶極源之間的夾角，此項(xiàng)涉及電阻率、頻率、場源與觀測點(diǎn)等相關(guān)參數(shù)。為提高廣域電磁視電阻率計(jì)算過程的穩(wěn)定性，本文將Brent-Dekker求根算法［42-43］引入廣域視電阻率的計(jì)算，并與二分法和直接迭代法［44］計(jì)算結(jié)果進(jìn)行對比。

假設(shè)一個(gè)電阻率為100 Ω·m 的均勻半空間，在原點(diǎn)設(shè)置一個(gè)沿x方向的場源，觀測點(diǎn)位于(0，5000 m，0)。文獻(xiàn)［43］對基于Brent-Dekker 算法的求解過程做了簡單描述，具體過程如下。

首先，假設(shè)在區(qū)間(a0，b0)內(nèi)求解函數(shù)f(x)，若f(a0)與f(b0)符號(hào)相反，則說明在區(qū)間(a0，b0)存在根，否則無根。利用該算法進(jìn)行迭代求解會(huì)涉及三個(gè)參數(shù)：ak-1，bk-1，bk-2，其中ak-1是第k-1 次迭代所得區(qū)域的下邊界，bk-1表示第k-1次迭代所得到的迭代點(diǎn)，bk-2是上一次迭代的迭代點(diǎn)或預(yù)估值。反復(fù)迭代，直到f(ak-1)與f(bk-1)符號(hào)相反，以保障(ak-1，bk-1)內(nèi)存在解，并滿足|f(ak-1)|＜|f(bk-1)|。

然后，對每一步迭代進(jìn)行逆二次插值

如果f(bk-1) 與f(bk-2)相等，則使用割線法

式中bk-2取值為ak-1。若要保證sk被接納，需要滿足兩個(gè)條件：①sk在或內(nèi)；②若上次迭代采用二分法計(jì)算，則必須滿足，若上次迭代采用逆二次插值法或割線法，則必須滿足，否則，不能保證求解過程的穩(wěn)定性。表1所示為均勻半空間模型不同算法廣域視電阻率計(jì)算結(jié)果。

表1 均勻半空間模型不同算法廣域視電阻率計(jì)算結(jié)果

從表1 可知，Brent-Dekker 算法相較于直接迭代法，在全區(qū)迭代次數(shù)相對穩(wěn)定、相同精度條件下，在過渡帶（8 Hz）和近區(qū)（1 Hz）的迭代次數(shù)明顯少于二分法。因此，為穩(wěn)定求解三維模型視電阻率，本文采用Brent-Dekker算法。

2 算例分析

2.1 層狀模型

為了驗(yàn)證算法的正確性，設(shè)計(jì)圖1 所示的三層層狀介質(zhì)模型。場源沿x方向布設(shè)，長度為10 m，中心坐標(biāo)為（0，0，0），發(fā)射電流為1 A，設(shè)置14個(gè)發(fā)射頻率，范圍為0.1～4096 Hz。沿x軸-6000～6000 m 范圍布設(shè)了一條觀測剖面，測點(diǎn)間距為200 m。求解區(qū)域?yàn)?［-40 km，40 km］3。該模型的卡尼亞視電阻率和廣域視電阻率剖面見圖2，0.1 Hz和8 Hz所對應(yīng)的視電阻率和相對誤差error曲線見圖3和圖4。

圖1 層狀模型示意圖

圖2 層狀模型卡尼亞視電阻率（ 上）和廣域視電阻率（下）剖面

圖3 層狀模型0.1 Hz 時(shí)1-th、3-th 以及13-th 網(wǎng)格下的卡尼亞（上）和廣域電磁（下）視電阻率（a）和相對誤差（b）曲線

圖4 層狀模型8 Hz 時(shí)1-th、3-th 以及13-th 網(wǎng)格下的卡尼亞（上）和廣域電磁（下）視電阻率（a）和相對誤差（b）曲線

圖3 和圖4 分別展示了不同層次（1-th，3-th，13-th）網(wǎng)格下得到的卡尼亞和廣域電磁視電阻率曲線及誤差分布。由圖可見，隨著網(wǎng)格密度的增加，計(jì)算精度明顯提升，相對誤差在4%以內(nèi)。13-th網(wǎng)格條件下得到視電阻率與解析解吻合最好，誤差相對較小。還可以看出，卡尼亞視電阻率相對廣域電磁視電阻率受非平面波的影響更大。圖5為以電流密度法向連續(xù)的自適應(yīng)有限技術(shù)細(xì)化三種層次網(wǎng)格（1-th、3-th、13-th）后的加權(quán)誤差β，可見隨著網(wǎng)格密度增加，誤差整體出現(xiàn)降低趨勢，符合預(yù)期。

圖5 層狀模型不同網(wǎng)格剖分下的加權(quán)誤差β

2.2 梯形山模型

為了分析廣域視電阻率受地形影響的特征，設(shè)計(jì)圖6 所示的梯形山模型。地下半空間的電阻率為100 Ω·m，空氣部分的電阻率為1×108Ω·m；求解區(qū)域?yàn)?［-10 km，10 km］3； 場源沿y方向布設(shè)，場源中心坐標(biāo)為（5 km，0，0），長度為100 m，發(fā)射電流為1 A。在地面沿y方向布設(shè)了一條觀測剖面，觀測范圍為［-1500 m，1500 m］。分別計(jì)算0.1、2、4 Hz的電場、磁場、卡尼亞視電阻率以及廣域視電阻率（圖7和圖8）。

圖6 梯形山地形模型示意圖

圖7 梯形山地形模型電場分量Ey（上）和磁場分量Hx（下）的實(shí)部（a）和虛部（b）曲線

圖8 梯形山地形模型卡尼亞和廣域視電阻率曲線

圖7 為不同頻率下電場分量Ey和磁場分量Hx的實(shí)部和虛部，可見地形會(huì)引起電、磁場的變化，在地形起伏的拐點(diǎn)位置畸變最明顯。該模型的正地形相當(dāng)于一個(gè)低阻異常體，對電流有吸引作用，因而地形正上方的曲線向下彎曲。

圖8 是不同頻率卡尼亞視電阻率和廣域視電阻率曲線，可見隨著頻率增高，卡尼亞視電阻率曲線出現(xiàn)抬升特征，數(shù)據(jù)進(jìn)入了近區(qū)，無法真實(shí)客觀地反映地下介質(zhì)電阻率信息。然而，整個(gè)頻率段的廣域視電阻率值受場源影響較小，一定程度上能夠避免卡尼亞視電阻率因場源影響而被抬升的情況。頻率為0.1 Hz時(shí)此現(xiàn)象最明顯，幾乎在每一個(gè)觀測點(diǎn)上這兩種視電阻率的差值均達(dá)到一個(gè)數(shù)量級(jí)。與卡尼亞視電阻率曲線相比，廣域視電阻率值更接近真實(shí)情況，說明了后者具有更大勘探深度，具備近區(qū)勘探能力。

2.3 球狀異常體模型

設(shè)計(jì)圖9 所示的球狀異常體模型。計(jì)算區(qū)域?yàn)閇-10 km，10 km]3。沿x方向布設(shè)線性源，源中心坐標(biāo)為(0，-5 km，0)，長度為1 m，發(fā)射電流為1 A。在地面設(shè)置一個(gè)正方形觀測區(qū)域：x=[-1000 m，1000 m]，y= [-1000 m，1000 m]。分別計(jì)算1、2、8、16 Hz 的電場分量Ex、磁場分量Hy、卡尼亞視電阻率以及Ex的廣域視電阻率，結(jié)果見圖10～圖13。

圖9 球狀異常體模型示意圖

圖10 球狀異常體模型1 Hz（a）和2 Hz（b）電場Ex（上）和磁場Hy（下）分布

圖11 球狀異常體模型1 Hz（a）和2 Hz（b）卡尼亞視電阻率（上）和廣域視電阻率（下）分布

圖12 球狀異常體模型8 Hz（a）和16 Hz（b）電場分量Ex（上）和磁場分量Hy（下）分布

圖13 球狀異常體模型8 Hz（a）和16 Hz（b）卡尼亞視電阻率（上）和廣域視電阻率（下）分布

由圖10和圖11可見，越靠近場源，電磁場越強(qiáng)，符合物理規(guī)律。磁場對地下異常體的敏感度弱于電場，因而電場Ex受低阻異常體的影響出現(xiàn)等值線彎曲現(xiàn)象，而磁場Hy則基本不受影響，即：低阻異常體能夠吸引電流，而對磁場則沒有影響。圖中兩種視電阻率等值線形態(tài)區(qū)別較大，主要是受場源的影響，卡尼亞視電阻率在靠近場源區(qū)域已經(jīng)進(jìn)入近區(qū)或過渡區(qū)，因而電阻率明顯高于背景電阻率，其電阻率信息不能反映真實(shí)的地下電阻率情況；廣域視電阻率則完全不同，在靠近場源的區(qū)域不會(huì)過早進(jìn)入近區(qū)或過渡區(qū)，因而其視電阻率信息能夠較真實(shí)地反映地下介質(zhì)的地電信息。因此，廣域電磁法能在一定程度上避免遠(yuǎn)區(qū)的限制，具有更大的勘探深度。

圖12展示了球狀異常體模型在8、16 Hz時(shí)的電場分量Ex、磁場分量Hy分布，圖13是對應(yīng)的卡尼亞視電阻率和廣域視電阻率平面等值線圖。從圖12可以看出，電場和磁場平面分布特征與圖7類似，不同之處在于隨著頻率增高，電場幅值增加，而磁場幅值減小。圖13中的視電阻率等值線分布特征與圖8類似，不同的是卡尼亞視電阻率受非平面波的影響變?nèi)?，僅較少的測點(diǎn)進(jìn)入了近區(qū)，大部分測點(diǎn)的視電阻率接近真實(shí)值。因此，廣域視電阻率能夠很好地反映低阻異常體的存在，在異常體正上方呈現(xiàn)低阻特征。

總之，在同等收發(fā)距的前提下，廣域電磁法摒棄了“遠(yuǎn)區(qū)”的限制，不受非平面波的影響，與CSAMT相比具有更大的勘探范圍和深度，得到視電阻率能夠較真實(shí)地反映地下介質(zhì)的電性分布，適用性更強(qiáng)，準(zhǔn)確度更高。

3 結(jié)論

本文采用基于電流密度法向連續(xù)為后驗(yàn)策略的自適應(yīng)有限元算法，實(shí)現(xiàn)了三維廣域電磁法正演計(jì)算。通過模型測試得出以下幾點(diǎn)認(rèn)識(shí)：

（1）針對廣域視電阻率的計(jì)算，采用了一種穩(wěn)健的Brent-Dekker 求根算法，該算法結(jié)合二分法、割線法等算法的優(yōu)勢，使求根過程更穩(wěn)健，迭代次數(shù)相對適中，確保了計(jì)算結(jié)果的正確性。

（2）將電流密度法向連續(xù)的后驗(yàn)誤差估計(jì)策略的自適應(yīng)有限元算法引入三維廣域電磁法正演模擬，提高了廣域視電阻率的計(jì)算精度。通過與解析解對比，驗(yàn)證了算法的正確性。該算法測試發(fā)現(xiàn)，場源附近單元受場源影響較大，會(huì)造成場源附近單元的過度細(xì)化，影響了場源積分的求解效率。

（3）對球狀異常體和梯形山地形等典型模型進(jìn)行正演計(jì)算，分析了廣域電磁法的探測深度與范圍。在同等情況下，CSAMT 因受到非平面波的影響，視電阻率曲線發(fā)生抬升現(xiàn)象，不能很好地反映真實(shí)的地下介質(zhì)信息，而廣域視電阻率摒棄了近區(qū)的限制，得到的視電阻率更真實(shí)可靠，拓展了有源電磁法的觀測范圍，同時(shí)增加了勘探深度。模型計(jì)算結(jié)果表明，地形對廣域電磁法和CSAMT 都有影響，尤其在地形拐點(diǎn)，畸變尤為明顯。因此，在進(jìn)行資料解釋時(shí)需要考慮地形的影響，盡量降低地形影響帶來的解釋誤差。