殷哲琦，張志勇，黃臨平，藍(lán)澤鸞，周 乾

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

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

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

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

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

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

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

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

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

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

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

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

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

求解上式即可得到節(jié)點(diǎn)的場(chǎng)值。

2 網(wǎng)格剖分

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

3 鋼度矩陣計(jì)算

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

3.1 泛函式第一項(xiàng)單元積分

當(dāng)選擇雙線性插值，K1e矩陣元素為

當(dāng)選擇雙二次插值，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)第二項(xiàng)單元積分

當(dāng)選擇雙線性插值，則有

當(dāng)選擇雙二次插值，則有

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

3.3 泛函式源項(xiàng)計(jì)算

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

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

當(dāng)選擇雙二次插值

其中l(wèi)為計(jì)算邊的長(zhǎng)度。

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

4 算例分析

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

圖3 計(jì)算模型Fig.3 Model of test

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

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

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

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

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

5 結(jié)論

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

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

陳樂(lè)壽.1981.有限元法在大地電磁場(chǎng)正演計(jì)算中的應(yīng)用及改進(jìn)［J］.地球科學(xué)，(2):241-260.

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

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

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

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

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

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

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

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

張繼鋒，湯井田，王燁，等.2010.基于預(yù)處理共軛梯度的大地電磁快速正演［J］.中南大學(xué)學(xué)報(bào):自然科學(xué)版，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.