李丁山 ,屈文璋 ,許 誠 ,欒曉東 ,孫 溥 ,姚 斌 ,冉啟順,潘 笑,卓賢軍,郭 銳,閆建峰

(1.中國艦船研究院,北京,100192;2.中國自然資源航空物探遙感中心,北京,100083;3.中國科學(xué)院地質(zhì)與地球物理研究所,北京,100029)

0 引言

目前,在應(yīng)用上較為成熟的水下目標(biāo)探測手段是水聲探測和被動(dòng)磁異常探測。水聲相對(duì)于其他信號(hào)在海水中具有衰減慢、分辨率較高等特點(diǎn),吊放聲吶和聲吶浮標(biāo)是典型的水聲探測裝備。被動(dòng)磁異常探測技術(shù)具有連續(xù)工作、反應(yīng)迅速及定位精度高等優(yōu)點(diǎn),可作為聲吶探測發(fā)現(xiàn)目標(biāo)后的輔助探測手段,用于目標(biāo)確定及精確定位[1]。目前已有多種型號(hào)磁探儀(如氦光泵磁探儀、銫光泵磁探儀、核磁共振式磁探儀及超導(dǎo)磁探儀等)投入使用。但是,隨著水下目標(biāo)聲、磁隱身技術(shù)的快速發(fā)展,傳統(tǒng)聲探測和被動(dòng)磁探測效能大幅降低。

基于上述問題,聚焦于利用航空瞬變電磁法(airborne transient electromagnetic,ATEM)對(duì)水下目標(biāo)進(jìn)行主動(dòng)電磁探測的方法研究。ATEM 具備主動(dòng)探測能力,對(duì)高導(dǎo)體的反應(yīng)敏感、探測速度快,可實(shí)現(xiàn)大面積快速勘探。在地球物理勘探領(lǐng)域,該方法通過分析具有不同電性的地下礦藏和圍巖之間所感應(yīng)出的二次場信號(hào),可得到地下電性結(jié)構(gòu)分布信息,從而實(shí)現(xiàn)地下資源探測[2]。同理,由于海水的均勻?qū)щ娦?而水下高導(dǎo)目標(biāo)體相對(duì)海水呈現(xiàn)明顯低阻[3],這為利用航空瞬變電磁法探測水中高導(dǎo)體提供了良好的物性基礎(chǔ)。

文中擬從ATEM 基本原理出發(fā),分別基于一維和三維模型進(jìn)行正演計(jì)算,研究不同仿真環(huán)境下水中高導(dǎo)體的垂直磁場響應(yīng)特征,并結(jié)合不同信噪比數(shù)據(jù)模擬實(shí)際干擾進(jìn)行反演計(jì)算,探索該方法在水下目標(biāo)探測領(lǐng)域的可行性。

1 ATEM 基本原理

ATEM 是一種電磁感應(yīng)原理為基礎(chǔ)的主動(dòng)源電磁探測方法[4]。通常是將電磁探測設(shè)備(包含發(fā)射線圈、接收線圈及發(fā)射機(jī)等)搭載在飛行平臺(tái)上利用發(fā)射線圈(磁偶源)向海水中發(fā)射不同波形的電磁波(如三角波、半正弦波等),在一次場關(guān)斷期間,線圈接收由水中高導(dǎo)體感應(yīng)產(chǎn)生的二次渦流場,該二次場含有豐富的高導(dǎo)體信息,通過分析接收線圈內(nèi)感應(yīng)電動(dòng)勢的變化,進(jìn)而識(shí)別水中高導(dǎo)體的位置和尺寸特征?；驹砣鐖D1 所示,紅色長方體代表高導(dǎo)體,虛線代表電磁場。

2 一維仿真計(jì)算

為探索航空瞬變電磁法在水下目標(biāo)探測領(lǐng)域的可行性,從基本原理層面出發(fā),首先構(gòu)建海水—高導(dǎo)電層—海水—海底巖石的4 層層狀模型,再采用數(shù)值濾波方法進(jìn)行數(shù)值計(jì)算,最終研究該方法能否對(duì)水下高導(dǎo)電層進(jìn)行辨別。

2.1 模型構(gòu)建

圖1 中高導(dǎo)體存在于海水之中、海底巖石之上,將此場景簡化成在一維層狀模型中某一深度存在高導(dǎo)電層,據(jù)此構(gòu)建層狀模型。直角坐標(biāo)系中,層狀大地在x和y方向是均勻無限的,以海平面為z=0,向下為正。其中: σ1,σ2,···,σn分別代表每層的電導(dǎo)率;d1,d2,···,dn分別代表每層的厚度,采用中心回線裝置,圓形發(fā)射線圈半徑為r,接收線圈設(shè)置在回線的中心,距離海面高度為h(h取值為正),如圖2 所示。

2.2 正演理論

對(duì)于圓形發(fā)射中心回線源而言,該裝置下接收線圈中心處垂直磁場響應(yīng)Hz是關(guān)于頻率 ω的函數(shù),可用下式表達(dá)[5]

式中:J1(·)為1 階貝塞爾函數(shù);I為發(fā)射電流;a為發(fā)射線圈半徑;λ為積分變量;rTE為反射系數(shù),和模型介質(zhì)參數(shù)有關(guān)。

對(duì)上式頻率域作傅里葉逆變換可得到時(shí)間域的垂直磁場響應(yīng)

式中:t為時(shí)間;Δ為采樣間隔;Wn為濾波系數(shù),n表示濾波系數(shù)個(gè)數(shù)。

實(shí)測數(shù)據(jù)感應(yīng)電動(dòng)勢與磁場分量存在以下關(guān)系

式中:SRx為接收線圈的面積;H為磁場強(qiáng)度分量。

2.3 垂直磁場響應(yīng)特征分析

基于一維層狀模型,采用47 點(diǎn)1 階漢克爾系數(shù)和300 點(diǎn)正弦數(shù)值濾波系數(shù)進(jìn)行數(shù)值計(jì)算[6],設(shè)計(jì)了2 種仿真場景,分別將高導(dǎo)電層放在不同深度,收發(fā)線圈放置在不同高度,計(jì)算不同場景下接收線圈中心處垂直磁場響應(yīng)值dBz/dt。

1) 不同目標(biāo)層深度時(shí)垂直磁場響應(yīng)

根據(jù)目標(biāo)層的高導(dǎo)電屬性,將其電阻率設(shè)為0.000 05 Ω·m,目標(biāo)層厚1 m,嵌入海水(深度300 m)之中,其中心距離海面分別設(shè)為5、10、25、50、75 和100 m;以中心距離海面5 m 為例,該層狀模型參數(shù)見表1。

發(fā)射和接收線圈在海面15 m 高處,圓形發(fā)射線圈半徑15 m,發(fā)射電流110 A 時(shí),得到該目標(biāo)層在不同深度的垂直磁場響應(yīng)衰減曲線,見圖3。從圖3 中可以看到,隨著高導(dǎo)電層深度的增加,所觀測到瞬變電磁響應(yīng)曲線越接近于背景場,當(dāng)深度為100 m 時(shí),其響應(yīng)曲線與背景曲線基本重合,意味著此時(shí)并不能有效探測到高導(dǎo)層。這主要是因?yàn)殡S著目標(biāo)層與收發(fā)線圈之間距離的增大,目標(biāo)體受發(fā)射線圈一次場的影響在減小,由此生成的二次場也在迅速減小。但整體可以看出,該方法對(duì)一定深度內(nèi)的高導(dǎo)電層具有明顯的分辨能力。

2) 不同收發(fā)高度時(shí)垂直磁場響應(yīng)

以目標(biāo)層在75 m 水深處為模型,分別將收發(fā)線圈設(shè)定在水面以上10、5 和0.5 m 處進(jìn)行信號(hào)發(fā)射與接收,發(fā)射磁矩(發(fā)射電流大小與線圈面積的乘積)不變,得到不同收發(fā)高度時(shí)的垂直磁場響應(yīng)曲線,如圖4 所示。

圖4 不同收發(fā)高度時(shí)垂直磁場響應(yīng)對(duì)比圖Fig.4 Comparison of vertical magnetic field responses at different transmitting and receiving heights

圖4 中可以看出,3 種收發(fā)高度下有無高導(dǎo)層得到的垂直磁場響應(yīng)曲線在早期時(shí)兩者基本重合,但在晚期4.9 ms 左右時(shí)兩者會(huì)有明顯的分叉走勢,說明可以明顯分辨海水中高導(dǎo)層。值得考慮的是,在晚期4.9 ms 之后,隨著收發(fā)高度增大,其二次場值逐漸減小,實(shí)際情況下,二次場可能會(huì)淹沒在背景噪聲之中。因此實(shí)際應(yīng)用中需要設(shè)法增強(qiáng)二次場值,以提高信噪比。

3 三維仿真計(jì)算

文中從一維層狀模型出發(fā),對(duì)海水中有無高導(dǎo)電層進(jìn)行仿真研究,結(jié)論表明該方法一定范圍內(nèi)能夠分辨出高導(dǎo)電層,但并未考慮三維目標(biāo)體尺寸、規(guī)模等因素,故從三維模型正演出發(fā)對(duì)水下三維目標(biāo)體垂直磁場響應(yīng)特征進(jìn)行深入研究。

3.1 三維模型正演方法

海水中的三維脈沖響應(yīng)可以通過在空間中數(shù)值求解麥克斯韋方程組得到。對(duì)比頻率-時(shí)間轉(zhuǎn)換算法,利用時(shí)域有限體積法直接在時(shí)域上求解麥克斯韋方程組[7],得到海水的脈沖響應(yīng),該方法可以在早期階段達(dá)到較高的仿真精度。對(duì)于激勵(lì)源的選取,考慮到實(shí)際的脈沖發(fā)射系統(tǒng)有一定的帶寬,因此將理想的脈沖源換作高斯源[8]。對(duì)于空間網(wǎng)格的剖分,采用典型Yee 方法[9],采用不等間距的六面體將邊緣電場和面磁通量離散化(見圖5),電場分量E在邊緣(藍(lán)色箭頭)上定義,磁通分量B在面(橙色箭頭)上定義。文中給出了電偶源下電場離散表達(dá)式,且每個(gè)時(shí)間步長的電場值利用迭代法求解,即

圖5 三維空間離散化Fig.5 Spatial discretization in three-dimensional space

式中:e為離散電場;Mfμ和Meσ分別為磁導(dǎo)率和電導(dǎo)率矩陣;C為離散旋度矩陣。

3.2 正演計(jì)算

在三維仿真計(jì)算中,利用圖1 中圓形中心回線形式進(jìn)行信號(hào)收發(fā),發(fā)射線圈半徑r=15 m,發(fā)射電流分別選用1 A 和10 A;發(fā)射/接收線圈離海面距離h=10 m;高導(dǎo)體長75 m,寬6 m,高8 m;電阻率0.000 05 Ω·m,高導(dǎo)體中心離海面距離d=10 m,高導(dǎo)體軸線和線圈的直徑在同一平面內(nèi)。據(jù)此,計(jì)算了高導(dǎo)體中心正上方存在高導(dǎo)體和不存在高導(dǎo)體時(shí)的dBz/dt響應(yīng),如圖6 所示。

圖6 不同電流下垂直磁場響應(yīng)對(duì)比圖Fig.6 Comparison of vertical magnetic field responses under different currents

圖6 中可以看出: 電流1 A 和10 A 情況下,有無目標(biāo)體得到dBz/dt衰減曲線在4 ms 左右均有明顯的分叉現(xiàn)象,說明該方法對(duì)高導(dǎo)體有明顯的異常響應(yīng)特征;同時(shí),發(fā)現(xiàn)當(dāng)發(fā)射電流變大,其二次場信號(hào)值也變大,信噪比也隨之提高。

4 反演成像目標(biāo)檢測

通過前面不同算例的正演計(jì)算,證明了該方法對(duì)高導(dǎo)體具有明顯的異常響應(yīng)特征,這為后續(xù)反演計(jì)算奠定了堅(jiān)實(shí)的基礎(chǔ)。反演成像是通過觀測數(shù)據(jù)進(jìn)行反演擬合得到介質(zhì)模型,然而滿足擬合閾值的模型眾多,需要在反演過程中增加正則化約束條件,以減少反演的多解性,實(shí)現(xiàn)目標(biāo)介質(zhì)屬性(尺寸、位置及電阻率)的有效檢測[10]。OCCAM反演算法具有對(duì)初始模型依賴程度小,運(yùn)算穩(wěn)定收斂等優(yōu)點(diǎn),能夠?qū)ふ以跇O小目標(biāo)體下符合數(shù)據(jù)的光滑模型[11],因此文中將該算法用來反演計(jì)算進(jìn)行目標(biāo)檢測。

4.1 OCCAM 反演基本原理

OCCAM 算法是一種基于光滑模型約束的算法[12],在數(shù)據(jù)擬合差和模型光滑度之間尋求最佳權(quán)衡,使得目標(biāo)函數(shù)值最小,可以表示為

式中:U(m)為目標(biāo)函數(shù)值;R(m)為模型約束項(xiàng),也稱模型粗糙度;λ為拉格朗日因子;‖dobs-F(m)‖為數(shù)據(jù)擬合項(xiàng),其中dobs為觀測數(shù)據(jù),F(m)為正演數(shù)據(jù);W為數(shù)據(jù)歸一化的對(duì)角矩陣;X*為‖Wdobs-WF(m)‖2項(xiàng)的期望值。

使?U/?mi+1=0即可,通過求解方程組得到模型的一般解為

反演的計(jì)算過程就是每一次將式(9)得到的模型參數(shù)mi+1代入目標(biāo)函數(shù)式(6)中計(jì)算,不停迭代尋找出滿足閾值條件的最佳模型參數(shù)為止[13]。該算法快速精準(zhǔn)成像的關(guān)鍵在于拉格朗日因子的合理選取、初始模型結(jié)構(gòu)的設(shè)定以及雅克比矩陣的快速求解。

4.2 噪聲分析

為了研究噪聲對(duì)目標(biāo)檢測效果的影響,以一維層狀模型為例,假設(shè)正演計(jì)算得到時(shí)間采樣窗口下的dBz/dt值所受噪聲服從高斯分布,且不同時(shí)間窗口該值互不相關(guān),則可以認(rèn)為dBz/dt值所受的噪聲服從高斯白噪聲[14]。因此,在正演仿真數(shù)據(jù)中疊加不同信噪比的高斯白噪聲,討論其對(duì)目標(biāo)檢測效果的影響。

其中,發(fā)射條件為圓形發(fā)射線圈半徑30 m,收發(fā)線圈距離海面10 m,發(fā)射電流大小100 A。模型總共3 層,第1 層和第3 層均為海水層,第2 層為高導(dǎo)電目標(biāo)層(0.005 Ω·m),高導(dǎo)薄層深度及厚度均為5 m,第3 層的∞代表半空間。觀測數(shù)據(jù)信噪比分別為10、20、30 和40 dB,不同信噪比的數(shù)據(jù)如圖7 所示,反演結(jié)果如圖8 所示。

圖7 不同信噪比數(shù)據(jù)質(zhì)量對(duì)比圖Fig.7 Data quality comparison of different SNRs

圖8 不同信噪比數(shù)據(jù)反演結(jié)果圖Fig.8 Inversion results of different SNR data

從圖7 中明顯看出,信噪比越低,正演數(shù)據(jù)所受干擾越嚴(yán)重,觀測數(shù)據(jù)質(zhì)量越低。圖8 中紅色線段表示正演模型,藍(lán)色線段代表不同信噪比下的反演結(jié)果,可以看出: 信噪比為40 dB 時(shí),可準(zhǔn)確反演出高導(dǎo)薄層的深度和厚度值;信噪比降為30 dB時(shí),反演效果變差,但是對(duì)深度值能給出較好的約束;當(dāng)信噪比減小為20 dB 和10 dB 后,雖然能對(duì)高導(dǎo)薄層有異常響應(yīng),但無法準(zhǔn)確判斷出深度及厚度信息。因此,認(rèn)為當(dāng)信噪比大于等于10 dB時(shí),利用OCCAM 反演算法對(duì)水中有無高導(dǎo)體能準(zhǔn)確判斷;信噪比越高,反演計(jì)算判斷效果越好;當(dāng)達(dá)到40 dB 時(shí),該方法能夠?qū)δ繕?biāo)層所處深度及厚度進(jìn)行有效檢測。

4.3 目標(biāo)體導(dǎo)電性對(duì)反演效果影響分析

在40 dB 信噪比情況下,將高導(dǎo)薄層的電阻率值分別設(shè)為4.2 節(jié)中高導(dǎo)薄層電阻率(0.005 Ω·m)的1/10(0.000 5 Ω·m)和1/100(0.000 05 Ω·m),其他參數(shù)相同,對(duì)應(yīng)的反演電阻率值分別繪制見圖9。可以看出: 在逐漸降低高導(dǎo)薄層的電阻率后,OCCAM算法計(jì)算得到的反演電阻率值能對(duì)目標(biāo)層所處深度進(jìn)行有效檢測,但厚度信息卻無法辨別。同時(shí)看出,當(dāng)目標(biāo)層電導(dǎo)率值逐漸變大時(shí),得到的水下深部反演電阻率值逐漸接近高導(dǎo)目標(biāo)層的電阻率值,這主要是由于高導(dǎo)電層對(duì)深部海水電性信息的屏蔽作用。

圖9 40 dB 下目標(biāo)體不同電阻率值反演結(jié)果圖Fig.9 Inversion results of different resistivity values of target at 40 dB

5 結(jié)論

通過一維計(jì)算可知: 1) 在一定深度內(nèi)、一定收發(fā)高度下,該方法對(duì)水中高導(dǎo)電層具有明顯的分辨能力;2) 考慮實(shí)際噪聲情況下,應(yīng)該選取更優(yōu)的裝置參數(shù)(如收發(fā)線圈高度更低),獲得更大的二次場值,保證更佳的分辨能力。通過三維模型正演計(jì)算可知: 有無目標(biāo)體得到dBz/dt衰減曲線有明顯的分叉現(xiàn)象,說明該方法對(duì)高導(dǎo)體有明顯的異常響應(yīng)特征。

利用不同信噪比合成觀測數(shù)據(jù)進(jìn)行反演計(jì)算,可以得到: 1) 當(dāng)信噪比大于等于10 dB 時(shí),OCCAM反演算法能對(duì)水中有無高導(dǎo)體進(jìn)行準(zhǔn)確判斷;2)信噪比越高,反演計(jì)算效果越好;3) 當(dāng)達(dá)到40 dB時(shí),該算法能夠?qū)δ繕?biāo)層所處深度及厚度進(jìn)行有效檢測;4) 高導(dǎo)薄層的電阻率逐漸減小后,OCCAM算法仍能對(duì)目標(biāo)體所處深度進(jìn)行有效檢測,但厚度信息卻無法辨別。

綜上,采用航空瞬變電磁法對(duì)水下高導(dǎo)體探測在理論上是可行的,在未來可考慮作為水下目標(biāo)探測的補(bǔ)充手段。但是該方法到應(yīng)用還需要剔除水下目標(biāo)體的激電效應(yīng),否則會(huì)導(dǎo)致在實(shí)測數(shù)據(jù)中經(jīng)常出現(xiàn)符號(hào)反轉(zhuǎn)現(xiàn)象;同時(shí)還需考慮裝置形式所存在的運(yùn)動(dòng)噪聲該如何克服等核心問題;對(duì)于反演成像問題,當(dāng)數(shù)據(jù)信噪比得到提高后,可考慮通過并行計(jì)算快速獲取目標(biāo)模型。