有限單元法求解水中二維體問題正演研究

林文東， 張志勇， 劉慶成， 湯洪志， 周 峰

(1.東華理工大學 核工程與地球物理學院，撫州 344000；2.中南大學 地球科學與信息物理學院，長沙 410083；3.東華理工大學 核資源與環境教育部重點實驗室，南昌 330013)

0 引言

潛艇是現代海軍的重要組成部分，近年來，潛艇的發展十分迅速，隱蔽性更強，使反潛變得越來越困難。早期以聲響訊號探測水下的人造物體被應用的最為廣泛，但隨著減聲降噪技術的發展，聲納技術遇到嚴峻挑戰，且由于海洋環境的復雜，使靈敏度受到一定影響，同時聲納探測還有自身的諸如“聲影區”的局限，非聲探測技術將發揮重要的作用[1]。

磁異常探測技術是目前潛艇非聲探測中得到應用最為廣泛的，該技術具有分類能力好、定位精度高、執行時間短、成本低、不受淺海復雜環境影響等優點。但由于海洋背景磁場形成很強的磁場噪聲，因此必須研究海洋背景磁場的特點和規律，以便從復雜的磁場信號中提取艦船或水下目標產生的磁場信號[2]。磁異常探測受自然條件影響較大，如霧、雨、雪、低云和風浪超過5級時，探測距離明顯下降，虛警率也上升。當云層下緣低于100 m或水平視距小于1 000 m時，基本不能使用。同時隨著潛艇消磁材料的發展，磁異常探測的效果也有所下降。

相對于傳統潛艇探測技術，電磁法探測是近幾年才開始起步的，夏立新譯著的文章介紹了瞬變電磁法(TEM)探測掩埋水雷的可行性；孫東平等闡述了復雜電磁環境及其對潛艇隱蔽性的影響；傅金祝[3]給出了利用電磁法探測的初步研究成果。目前，國內相關研究大多數處于理論研究，實用性不足，并且缺少物理模型實驗性的技術研究工作，極大地束縛了反潛裝備的研發。

天然電磁場是由太陽微粒輻射(太陽風)作用下形成的地球磁層和電離層的變化，在地球表面有限區域內可視為平面波，其傳播深度主要依賴于振動頻率或者場的變化，大地電磁測深是通過對地面電磁場的觀測，來研究地下巖礦石電阻率的分布規律的一種物探方法，其理論的關鍵是研究地面電磁場與地下巖礦石的電阻率存在的關系，電磁測深法廣泛應用于地球內部結構、深部地質研究、深部找礦等領域，作者將電磁法探測應用到海洋中，同時引入有限單元法來研究水中二維體在天然電磁場中的異常場，以達到探測目的。有限單元的計算包括四個主要環節：①建立與所研究問題相符的數學物理方程與邊界條件；②對計算模型進行離散化處理；③依據邊值問題所滿足的泛函，通過選擇離散單元上的插值函數，形成鋼度矩陣；④對稀疏線性系統進行求解，得到邊值問題的解。相對于只考慮電阻率的情況，作者在此基礎上加入了磁導率的概念，效果優于只考慮電阻率的情況。

圖1 全球電場、磁場強度平均振幅譜特征圖(Compbell,1967)Fig．1 Typical average amplitude spectrum of geomagnetic variations (Compbell,1967)

對于水下二維體的電磁法探測，有限單元的離散單元采用三角形剖分，三角單元可以更好地擬合復雜地形與地質單元。為了提高擬合效果，Key[4]等采用非結構化的自適應剖分算法進行海洋大地電磁的正演；Li[5]等利用非結構化網格進行了各向異性問題的二維大地電磁正演。非結構化網格剖分算法，有利于對地質單元與地形的擬合，且最大可能地保證剖分單元的質量、減少擴展區域單元數量。通過模型試算驗證了二維大地電磁有限元算法在水下二維體探測中具有較高的精度和適用性，可以為電磁法探測水下潛艇打下堅實的基礎。

1 天然電磁場下水中二維體的有限單元數值模擬基本原理

圖1為全球電場、磁場強度平均振幅譜特征圖，是1967年Compbell的研究成果，從圖1可以看出，大地電磁測深法所觀測的電磁場信號十分微弱，電場振幅最低僅為0.01 mV/km，磁場的振幅最低為0.001 nT，因此要求觀測儀器要有很高的精度，同時有效地識別、抑制干擾噪聲也是至關重要的。

根據麥克斯韋方程，角頻率為ω(時間因子為e-iωt)的定態電磁場方程是[6-7]

▽×E=iωμH

(1)

▽×H=(σ-iωε)E

(2)

其中μ是介質的磁導率；σ是電導率；ε是介電常數。

對于海中二維體，取走向為z軸，x軸與z軸垂直，y軸垂直向上，當平面電磁波以任何角度入射海面時，海水中的電磁波總以平面波形式，幾乎垂直的向下傳播，將式(1)與式(2)按分量展開，可得到兩個獨立方程組，命名為E型波和H型波，也稱TE模式和TM模式[6]。

E型：

(3)

H型：

(4)

從式(3)和式(4)中的第二、三式中解出Hx、Hy和Ex、Ey，并分別帶入第一式，得到偏微分方程[6-8]：

(σ-iωε)Ez=0

(5)

iωμHz=0

(6)

式(5)和式(6)可以統一表示成

▽·(τ▽u)+λu=0

(7)

(8)

u|Γ top= 1

(9)

δF(u)=0

(10)

其中 Ω為研究區域；Γtop為研究上邊界，E型波取空氣邊界、H型波取海面邊界；Γbottom為研究區域底邊界。

根據有限單元剖分結果，泛函式(8)可表示為離散形式[6,9]

(11)

根據變分式子(10)，得到有限單元計算方程組

Ku=0

(12)

將邊界條件式(3b)代入式(5)求解，即可得到節點的場值。

2 網格的三角剖分方法

鑒于各種地表的地球物理觀測方法，一般都是淺部觀測數據的分辨率高，到深部逐漸變差，基于這樣的認識，在研究海中問題時，作者引入網格的三角剖分方法,剖分方法引自文章[10]。

圖2 三角剖分Fig．2 Triangle split

3 水下二維體的電磁剖分單元鋼度矩陣計算

網格剖分的最后單元為一次和二次插值三角形，利用一維、二維自然坐標[6]推導可得到單元上的鋼度矩陣。三角單元節點逆時針編號i、j、m，若有二次插值節點則為i、j、m、r、p、q，其中r為邊i,j為中點；p為邊j,m中點、q為邊m,i中點。下文中a、b、c的定義來自參考文獻[6]。

3.1 泛函式(11)第一項單元積分

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

(13)

k,h∈{i,j,k}

(14)

w,ω,n∈{i,j,k}&n≠w≠ω

(15)

(16)

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

當選擇雙線性插值，則有

還原氧化石墨烯的形貌特征通過掃描電鏡和透射電鏡進行表征。圖3是還原氧化石墨烯的掃描電鏡圖，(a)和(b)圖是正面圖，能夠看到，經遠紅外光波制備的rGO表面是皺縮的，片層雖小但分布均勻。(c)和(d)圖是截面圖，能夠發現，rGO呈絮狀堆疊，片與片之間呈現較好的網狀結構，這是靜電噴霧以及遠紅外光波處理的結果，這種方法制備的石墨烯片層比較均勻，幾乎是單層結構，這大大防止了石墨烯的團聚，而且得到的材料擁有更大的比表面積，對在超級電容器上的應用提供了理論支撐。

(17)

當選擇雙二次插值，則有

(18)

3.3 泛函式(11)第三項單元積分

當選擇雙線性插值，則有

(19)

當選擇雙二次插值，則有

(20)

4 算例分析

計算區域水平方向距離50 m，垂直方向取以2為底的對數,垂向空氣剖分15層，海水剖分35層，模擬計算的空氣、海水和水下二維體的電阻率、介電常數、磁導率如表1所示。計算其在TE、TM模式下的視電阻率。

表1 空氣、海水和水下低阻體的電阻率、介電常數、磁導率

TE、TM計算結果如圖3至圖6所示。

圖3 TM模式，只考慮電阻率的視電阻率(a)與視相位(b)斷面圖Fig．3 TM model, the section map of apparent resistivity (a) and apparent phase (b) when only consider the resistivity

圖4 TM模式，同時考慮電阻率與磁導率的視電阻率(a)與視相位(b)斷面圖Fig．4 TM model, the section map of apparent resistivity (a) and apparent phase(b) when consider the resistivity and magneticpermeability

圖5 TE模式，只考慮電阻率的視電阻率(a)與視相位(b)斷面圖Fig．5 TE model, the section map of apparent resistivity (a) and apparent phase (b) when only consider the resistivity

圖6 TE模式，同時考慮電阻率與磁導率的視電阻率(a)與視相位(b)斷面圖Fig．6 TE model, the section map of apparent resistivity (a) and apparent phase (b) when consider the resistivity and magneticpermeability

圖3與圖4是在TM模式下的模擬計算結果，由計算結果可見，在只考慮電阻率的情況下，在上底30 m處左右形成低阻帶狀異常(圖3)；同時考慮電阻率與磁導率，在相同位置處形成的異常曲線是封閉的，但異常值明顯很大，可以很好地刻畫出異常區域(圖4)。TE模式的計算結果如圖5、圖6所示，圖5中在上底30 m處左右形成高阻異常；在同時考慮電阻率與磁導率的圖6中，在相同位置處形成的高阻異常值更大一些。

5 結論

作者研究了電磁法探測水中二維體的問題，基于大地電磁問題的有限單元正演，通過模型試算，驗證了二維大地電磁有限元算法在水下二維體探測中具有較高的精度和適用性，得出如下結論：

(1)引入二維大地電磁有限單元正演來探測水下二維體,具有較高的精度和適用性。

(2)TM模式下，在低阻時，同時引入磁導率和電阻率的效果比只考慮電阻率探測異常效果要好；在引入磁導率的情況下，可以更好地圈定異常區域，從而確定水下二維體的位置。

(3)TE模式下，在低阻時，同時引入磁導率和電阻率的效果比只考慮電阻率得到的異常值要大，但是也可以刻畫出存在異常的區域。

參考文獻：

[1] 申會堂,張玉敏.海洋中的探測新技術[J].艦船科學技術,2001,1(3):61-64.

[2] 祝傳剛,陳標,查玉龍，等. 基于海洋磁異常特征的目標定位理論和方法[J]. 青島大學學報：自然科學版, 2008, 21(2):40-47.

[3] 傅金祝.探測掩埋水雷和臨時爆炸裝置的磁和電磁方法[J].水雷戰與艦船防護,2011, 19(1):65-69.

[4] KEY K, WEISS C. Adaptive finite-element modeling using unstructured grids: The 2D magnetotelluric example[J].Geophysics, 2006,71(6): 291-299.

[5] LI Y, PEK J. Adaptive finite element modelling of two-dimensional magnetotelluric fields in general anisotropic media[J].Geophys.J.Int., 2008,175: 942-954.

[6] 徐世浙.地球物理中的有限單元法[M].北京:科學出版社, 1994.

[7] MOHAMMAD I H,SRIGUTOMO W,SUTARNO D.The modeling of 2D controlled source audio magnetotelluric(CSAMT) responses using finite element method[J].Journal of Electromagnetic Analysis and Applications,2012(4):293-304.

[8] TONG X,LIU J,XIE W.Three-dimensional forward modeling for magnetotelluric sounding by finite element method[J].J.Cent.South Ualv.Technol,2009,16:136-142.

[9] 石宇，羅天涯，熊彬，等.有限元法求解各向異性介質中大地電磁場邊值問題研究[J].地球科學前沿,2012(2):246-251.

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