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

(1.東華理工大學(xué) 核工程與地球物理學(xué)院，撫州 344000；2.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083；3.東華理工大學(xué) 核資源與環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，南昌 330013)

0 引言

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

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

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

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

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

對(duì)于水下二維體的電磁法探測(cè)，有限單元的離散單元采用三角形剖分，三角單元可以更好地?cái)M合復(fù)雜地形與地質(zhì)單元。為了提高擬合效果，Key[4]等采用非結(jié)構(gòu)化的自適應(yīng)剖分算法進(jìn)行海洋大地電磁的正演；Li[5]等利用非結(jié)構(gòu)化網(wǎng)格進(jìn)行了各向異性問(wèn)題的二維大地電磁正演。非結(jié)構(gòu)化網(wǎng)格剖分算法，有利于對(duì)地質(zhì)單元與地形的擬合，且最大可能地保證剖分單元的質(zhì)量、減少擴(kuò)展區(qū)域單元數(shù)量。通過(guò)模型試算驗(yàn)證了二維大地電磁有限元算法在水下二維體探測(cè)中具有較高的精度和適用性，可以為電磁法探測(cè)水下潛艇打下堅(jiān)實(shí)的基礎(chǔ)。

1 天然電磁場(chǎng)下水中二維體的有限單元數(shù)值模擬基本原理

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

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

▽×E=iωμH

(1)

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

(2)

其中μ是介質(zhì)的磁導(dǎo)率；σ是電導(dǎo)率；ε是介電常數(shù)。

對(duì)于海中二維體，取走向?yàn)閦軸，x軸與z軸垂直，y軸垂直向上，當(dāng)平面電磁波以任何角度入射海面時(shí)，海水中的電磁波總以平面波形式，幾乎垂直的向下傳播，將式(1)與式(2)按分量展開(kāi)，可得到兩個(gè)獨(dú)立方程組，命名為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)可以統(tǒng)一表示成

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

(7)

(8)

u|Γ top= 1

(9)

δF(u)=0

(10)

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

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

(11)

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

Ku=0

(12)

將邊界條件式(3b)代入式(5)求解，即可得到節(jié)點(diǎn)的場(chǎng)值。

2 網(wǎng)格的三角剖分方法

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

圖2 三角剖分Fig．2 Triangle split

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

網(wǎng)格剖分的最后單元為一次和二次插值三角形，利用一維、二維自然坐標(biāo)[6]推導(dǎo)可得到單元上的鋼度矩陣。三角單元節(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)。下文中a、b、c的定義來(lái)自參考文獻(xiàn)[6]。

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

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

(13)

k,h∈{i,j,k}

(14)

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

(15)

(16)

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

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

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

(17)

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

(18)

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

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

(19)

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

(20)

4 算例分析

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

表1 空氣、海水和水下低阻體的電阻率、介電常數(shù)、磁導(dǎo)率

TE、TM計(jì)算結(jié)果如圖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模式，同時(shí)考慮電阻率與磁導(dǎo)率的視電阻率(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模式，同時(shí)考慮電阻率與磁導(dǎo)率的視電阻率(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模式下的模擬計(jì)算結(jié)果，由計(jì)算結(jié)果可見(jiàn)，在只考慮電阻率的情況下，在上底30 m處左右形成低阻帶狀異常(圖3)；同時(shí)考慮電阻率與磁導(dǎo)率，在相同位置處形成的異常曲線是封閉的，但異常值明顯很大，可以很好地刻畫(huà)出異常區(qū)域(圖4)。TE模式的計(jì)算結(jié)果如圖5、圖6所示，圖5中在上底30 m處左右形成高阻異常；在同時(shí)考慮電阻率與磁導(dǎo)率的圖6中，在相同位置處形成的高阻異常值更大一些。

5 結(jié)論

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

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

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

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

參考文獻(xiàn)：

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

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

[3] 傅金祝.探測(cè)掩埋水雷和臨時(shí)爆炸裝置的磁和電磁方法[J].水雷戰(zhàn)與艦船防護(hù),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].北京:科學(xué)出版社, 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] 石宇，羅天涯，熊彬，等.有限元法求解各向異性介質(zhì)中大地電磁場(chǎng)邊值問(wèn)題研究[J].地球科學(xué)前沿,2012(2):246-251.

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