基于二次場算法的CSEM二維有限單元法正演

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

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

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

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

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

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

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

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

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

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

構造(3)式的泛函

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

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

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

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

2 網格剖分

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

3 鋼度矩陣計算

收縮網格剖分的最后單元為一次和二次插值三角形，利用自然坐標可推導得單元上的鋼度矩陣(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 基于二叉樹的剖分結構Fig.1 The split structure based on the binary tree

圖2 三角單元結構Fig.2 Triangle unit structure

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

當選擇雙線性插值，則有

當選擇雙二次插值，則有

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

3.3 泛函式源項計算

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

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

當選擇雙二次插值

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

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

4 算例分析

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

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

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

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

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

算例二:磁導率與電阻率模型，取異常體介電常數等于圍巖，圍巖電阻率ρ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。

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

5 結論

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

圖5 算例二計算結果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.基于二次場的二維大地電磁有限元法數值模擬［J］.同濟大學學報(自然科學版)，35(08):1113-1117.

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

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

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

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

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

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

張志勇，劉慶成.2013.基于收縮二叉樹結構網格剖分的大地電磁二維有限單元法正演研究.石油地球物理勘探，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.