殷哲琦，張志勇，黃臨平，藍澤鸞，周 乾

(東華理工大學核工程與地球物理學院，江西南昌 330013)

自Coggon(1971)將有限單元法引入地球物理領(lǐng)域起，有限單元法在地球物理電磁場計算問題中得到了廣泛應用;陳樂壽(1981)提出了以矩形單元為主結(jié)合使用部分三角形單元對有限單元法進行改進，提高了對分界面形狀的適應能力;Wannamaker等(1986)利用有限單元法對帶地形的二維大地電磁(MT)進行數(shù)值模擬;Bouchard(1988)對大地電磁進行了二維地形改正的研究，表明起伏地形會產(chǎn)生明顯的假異常，說明在進行MT數(shù)據(jù)解釋時應消除地形的影響;Key等(2006)人利用非結(jié)構(gòu)化自適應有限元網(wǎng)格對二維大地電磁進行研究，研究表明非結(jié)構(gòu)化和自適應算法可以很好地對復雜2D模型進行剖分，自適應算法采用的是一種后驗誤差估計，其目的是使網(wǎng)格剖分能更好的滿足實際模型，通過不斷迭代得到最佳網(wǎng)格剖分，減少計算誤差從而提高計算精度。國內(nèi)學者也對有限元算法進行了研究;徐世浙等(1995)利用有限單元法對電導率連續(xù)分層變化的層狀介質(zhì)的大地電磁進行了正演研究;為提高計算精度，史明娟等(1997)對基于二次插值四邊形網(wǎng)格剖分進行了大地電磁正演研究;柳建新等(2009)采用自適應有限元算法，通過非結(jié)構(gòu)化自動網(wǎng)格剖分技術(shù)，使得網(wǎng)格剖分較好擬合復雜地形的邊界;由于有限單元法進行數(shù)值模擬的網(wǎng)格邊界為截斷邊界，影響計算精度，張繼鋒等(2010)模擬了MT的截斷邊界對計算結(jié)果的影響。

本文利用二次場進行CSEM數(shù)值模擬，在非均勻介質(zhì)地電結(jié)構(gòu)中傳播電磁總場可以分為背景場和二次場的總和，研究工作從麥克斯韋方程組出發(fā)，導出了同時考慮電阻率與磁導率異常的可控源電磁法二次場方程;采用了基于二叉樹結(jié)構(gòu)的三角單元剖分，該剖分方法網(wǎng)格生成速度快、節(jié)點編號與訪問容易，具有較好的擬合地質(zhì)單元與地形的效果;并推導了三角單元雙線性與雙二次插值的單元鋼度矩陣表達式;分別對電阻率異常及磁導率異常模型進行了試算，結(jié)果表明二次場算法無需對場源區(qū)域進行剖分，可有效的減小計算區(qū)間，提高計算效率;另二次場算法邊界處理簡單，有利用提高計算精度，解決了低頻計算不穩(wěn)定問題。

1 基于二次場的有限單元法二維大地電磁場數(shù)值模擬原理

考慮地下電導率與磁導率變化，則二次電磁場滿足下面形式的麥克斯韋方程組(劉小軍，2007):

其中取場值負諧變(e－iωt)，ω為角頻率;Eb、Hb為一次電場與磁場，Es，Hs為二次電場與磁場=σiωε為復電導率，σ為電導率，σ=Δσ+σ0，其中Δσ為異常電導率、σ0背景電導率，ε為介電常數(shù)(取真空中介電常數(shù)ε0=8.8542×10－12F·m－1); μ磁導率、Δμ為異常磁導率，其中=iωμ= iωΔμ(取真空中磁導率μ0=4π×10－7H·m－1)。

二維條件下，取地質(zhì)體走向為y軸，場源為沿y軸向無限延伸的線電流源，則走向方向二次電場滿足偏微分方程(劉云，2010)，

式(2)為線源二維可控源二次場方程。

如視(2)式右端項為“源項”，則可寫成以下形式，

構(gòu)造(3)式的泛函

微分方程(3)的解與求解泛函(4)的極值問題等價。

對研究區(qū)域剖分，泛函式(4)可表示為離散形式(馬為，2008)

對(5)變分式，得到有限單元計算方程組

求解上式即可得到節(jié)點的場值。

2 網(wǎng)格剖分

鑒于各種地表的地球物理觀測方法，一般觀測數(shù)據(jù)對淺部分辨率高、向深部逐漸降低，基于這樣的認識，設計了基于二叉樹的收縮網(wǎng)格剖分算法(Li，2008)。算法原理，假設由模型底部向地表分辨率變高，采用四邊形剖分，則其具有二叉樹式的結(jié)構(gòu);剖分形成的四邊形可進一步三角化，實現(xiàn)三角剖分;利用對三角單元邊的訪問可實現(xiàn)二次插值;當然有二次插值結(jié)點的三角形還可進一步離散為四個三角形實現(xiàn)二次剖分，典型的剖分如圖1 (張志勇，2013)。

3 鋼度矩陣計算

收縮網(wǎng)格剖分的最后單元為一次和二次插值三角形，利用自然坐標可推導得單元上的鋼度矩陣(Liu，1990)。三角單元節(jié)點逆時針編號i，j，m，如果有二次插值節(jié)點則為i，j，m，r，p，q，其中r為邊i，j中點、p為邊j，m中點、q為邊m，i中點。如圖2所示，下文中a，b，c的定義來自徐世浙(1994)。

3.1 泛函式第一項單元積分

當選擇雙線性插值，K1e矩陣元素為

當選擇雙二次插值，K1e矩陣元素為

圖1 基于二叉樹的剖分結(jié)構(gòu)Fig.1 The split structure based on the binary tree

圖2 三角單元結(jié)構(gòu)Fig.2 Triangle unit structure

3.2 泛函式(4)第二項單元積分

當選擇雙線性插值，則有

當選擇雙二次插值，則有

式(7)～(8)中Δ為三角形單元的面積。

3.3 泛函式源項計算

nz表示外法線方向與z軸夾角余弦。

當選擇雙線性插值，K5e三角單元閉合曲面積分，其中邊上的積分矩陣元素為

當選擇雙二次插值

其中l(wèi)為計算邊的長度。

將K5e中nx換為nz即可得到K6e。

4 算例分析

在均勻半空間內(nèi)設置三個截面為長方形(寬400 m，高500 m)的二度體，二度體上頂埋深100 m，間距1 600 m，模型設置如圖3。在地表采用進行可控源電磁測量，觀測頻率為(1～214 Hz)。

圖3 計算模型Fig.3 Model of test

算例一:電阻率模型，取異常體磁導率、介電常數(shù)等于圍巖，圍巖電阻率ρ0=100 Ω·m，M1異常體ρ1=50 Ω·m，M2異常體ρ2=100 Ω·m，M3異常體ρ3=200 Ω·m。

計算結(jié)果如圖4，在只考慮電阻率影響條件下，M1形成低阻圈閉異常，最低電阻率等值線70 Ω·m;M3形成高阻圈閉異常，最高電阻率等值線130 Ω·m;過渡帶與近區(qū)數(shù)據(jù)計算穩(wěn)定。

圖4 算例一計算結(jié)果Fig.4 Contours of apparent resistivity and phase for first test

算例二:磁導率與電阻率模型，取異常體介電常數(shù)等于圍巖，圍巖電阻率ρ0=100 Ω·m，磁導率為μ0，M1異常體ρ1=50 Ω·m，μ1=2μ0，M2異常體ρ2=100 Ω·m，μ2=2μ0，M3異常體ρ3=200 Ω ·m，μ3=2μ0。

計算結(jié)果如圖5，同時考慮電阻率、磁導率影響條件下，三個異常體均形成明顯異常。M1形成低阻圈閉異常，最低電阻率等值線80 Ω·m;M2形成高阻圈閉異常，最高電阻率等值線120 Ω·m;M3形成高阻圈閉異常，最高電阻率等值線130 Ω·m;過渡帶與近區(qū)數(shù)據(jù)計算穩(wěn)定。磁導率對視電阻率明顯的影響，使視電阻率增大。

5 結(jié)論

本文引入二次場算法實現(xiàn)了二維CSEM正演模擬。研究了同時考慮電阻率與磁導率的可控源電磁法二次場方程，采用基于二叉樹結(jié)構(gòu)的收縮網(wǎng)格剖分，推導了三角單元雙線性與雙二次插值單元鋼度矩陣，實現(xiàn)了有限單元模擬。實踐表明，二次場算法無需對場源區(qū)域剖分，可有效的減小剖分區(qū)間，提高計算效率;與總場算法相比邊界處理簡單;解決可控源低頻計算不穩(wěn)定問題、提高了計算精度。對電阻率異常模型及同時考慮電阻率與磁導率異常的試算，試算表明磁導率對可控源視電阻率產(chǎn)生明顯影響，使電阻率變大。

圖5 算例二計算結(jié)果Fig.5 Contours of apparent resistivity and phase for second test

陳樂壽.1981.有限元法在大地電磁場正演計算中的應用及改進［J］.地球科學，(2):241-260.

陳小斌，張翔，胡文寶.2000.有限元直接迭代算法在MT二維正演計算中的應用［J］.石油地球物理勘探，35(4):487-496.

柳建新，蔣鵬飛，童孝忠，等.2009.不完全 LU分解預處理的BICGSTAB算法在大地電磁二維正演模擬中的應用［J］.中南大學學報(自然科學版)，40(2):484-491.

劉小軍，王家林，于鵬.2007.基于二次場的二維大地電磁有限元法數(shù)值模擬［J］.同濟大學學報(自然科學版)，35(08):1113-1117.

劉云，王緒本.2010.大地電磁二維自適應地形有限元正演模擬［J］.地震地質(zhì)，32(03):382-391.

馬為，陳小斌，趙國澤.2008.大地電磁測深二維正演中輔助場的新算法［J］.地震地質(zhì)，30(02):525-533.

史明娟，徐世浙，劉斌.1997.大地電磁二次函數(shù)插值的有限元法正演模擬［J］.地球物理學報，40(3):421-430.

徐世浙.1994.地球物理中的有限單元法［M］.北京:科學出版社: 229-247.

徐世浙，于濤，李予國，等.1995.電導率分塊連續(xù)變化的二維MT有限元模擬(Ⅰ)［J］.高校地質(zhì)學報，1(2):65-73.

張繼鋒，湯井田，王燁，等.2010.基于預處理共軛梯度的大地電磁快速正演［J］.中南大學學報:自然科學版，41(5):1877-1882.

張志勇，劉慶成.2013.基于收縮二叉樹結(jié)構(gòu)網(wǎng)格剖分的大地電磁二維有限單元法正演研究.石油地球物理勘探，48(3):482-487.

Bouchard M C A K.1988.Two-dimensional terrain correction in magnetotelluric surveys［J］.Geophysics，53(6):854-862.

Coggon J.H.1971.Electromagnetic and electrical modeling by the finite element method［J］.Geophysics，36(1):132-135.

Davis T.A.2006.Direct methods for sparse linear systems［M］.SIAM，99-133.

Key K，Weiss C.2006.Adaptive finite-element modeling using unstructured grids:The 2D magnetotelluric example［J］.Geophysics，71 (6):G291-G299.

Li Y，Pek J.2008.Adaptive finite element modelling of two-dimensional magnetotelluric fields in general anisotropic media［J］.Geophysical Journal International，175:942-954.

Liu J.1990.The Role of Elimination Trees in Sparse Factorization［J］.SIAM Journal on Matrix Analysis and Applications，11(1):134-172.

Wannamaker P.E，Stodt J.A，Rijo L.1986.Two-dimensional topographic responses in magnetotelluric modeled using finite elements［J］.Geophysics，51(11):2131-2144.

Wannamaker P E，Stodt J A，Rijo L.1987.A stable finite element solution for two-dimensional magnetotelluric modelling［J］.Geophysical Journal Research，88(01):277-296.