基于航空瞬變電磁法的水下高導體探測方法

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

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

0 引言

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

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

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

1 ATEM 基本原理

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

2 一維仿真計算

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

2.1 模型構建

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

2.2 正演理論

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

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

對上式頻率域作傅里葉逆變換可得到時間域的垂直磁場響應

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

實測數(shù)據(jù)感應電動勢與磁場分量存在以下關系

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

2.3 垂直磁場響應特征分析

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

1) 不同目標層深度時垂直磁場響應

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

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

2) 不同收發(fā)高度時垂直磁場響應

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

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

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

3 三維仿真計算

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

3.1 三維模型正演方法

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

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

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

3.2 正演計算

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

圖6 不同電流下垂直磁場響應對比圖Fig.6 Comparison of vertical magnetic field responses under different currents

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

4 反演成像目標檢測

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

4.1 OCCAM 反演基本原理

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

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

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

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

4.2 噪聲分析

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

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

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

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

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

4.3 目標體導電性對反演效果影響分析

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

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

5 結論

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

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

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