磁異常信號奇異值分解的隨機共振檢測方法

李啟飛，吳 芳，韓蕾蕾，范趙鵬，李沛宗

(1.海軍航空大學(xué)，山東 煙臺 264001；2.解放軍91550部隊，遼寧 大連 116000；3.解放軍91001部隊，北京 100000)

0 引言

磁異常探測作為水下目標(biāo)探測的有效手段，已被廣泛應(yīng)用于軍事領(lǐng)域。軍事領(lǐng)域主要使用航空器作為載體，搭載磁探測平臺，對所在空域進(jìn)行磁場探測。當(dāng)水下存在目標(biāo)時，目標(biāo)會造成附近磁場的異常波動。對該異常波動信號進(jìn)行處理，從而實現(xiàn)對水下目標(biāo)的檢測，是本文的研究內(nèi)容。磁異常探測領(lǐng)域研究主要分為兩個部分：一是水下目標(biāo)磁場模型的建立，主要方法包括邊界積分法[1-2]、磁偶極子法[3]、有限元法[4]、積分方程法[1, 5]等；二是磁異常信號的處理，側(cè)重于對水下目標(biāo)磁場信號的檢測，目前主流方法包括信號的時域特征分析、頻譜分析、時頻特性分析[6]、基于標(biāo)準(zhǔn)正交基函數(shù)[7](Orthonormalized Basis Function，OBF)分解、熵濾波器檢測[8]等。

在實際檢測過程中，傳感器采樣信號過低的信噪比，使得目標(biāo)信號淹沒在背景噪聲中，從時域上進(jìn)行分析，很難直接對微弱的磁異常信號進(jìn)行檢測。OBF算法只有在目標(biāo)信號和噪聲合適的情況下才能夠使用，熵濾波器對于低信噪比檢測能力較差。為了提高對水下目標(biāo)磁異常信號的檢測能力，本文針對上述問題，提出了磁異常信號奇異值分解的隨機共振檢測方法。

1 探測原理及結(jié)果

1.1 水下目標(biāo)磁場仿真

通過等效磁源對潛艇磁場的分布進(jìn)行模擬的方法叫做等效磁源法，主要的等效磁體模型有單個磁偶極子模型、旋轉(zhuǎn)橢球體模型和磁偶極子陣列模型等方法。文獻(xiàn)[9]中指出，在2.5倍物體長度以上的空間，完全可以把磁性物體視做磁偶極子處理，滿足一般工程應(yīng)用的需求。潛艇一般與反潛機距離百米以上，磁偶極子模型作為等效磁源對潛艇磁場進(jìn)行仿真，是完全可行的。

在航空磁探反潛過程中，潛艇作為鐵磁性目標(biāo)，其磁場H可以表示為：

(1)

(2)

通過計算，可以得到：

(3)

式(3)中，H為等效磁偶極子磁場強度，Hx,Hy,Hz是磁場強度的三分量，H=ixHx+iyHy+izHz，ix、iy、iz為坐標(biāo)系中的單位向量。

1.2 奇異值分解原理

在工程應(yīng)用中，由于噪聲的嚴(yán)重干擾，目標(biāo)信號往往淹沒在隨機的噪聲中，并且噪聲與目標(biāo)信號的頻域相互重疊，導(dǎo)致傳統(tǒng)的算法難以對低信噪比的信號進(jìn)行檢測。奇異值分解在信號處理領(lǐng)域有著獨特的作用，該方法具有理想的去相關(guān)特性[10]，能夠在強噪聲背景下提取有用信號和特征信息[11-12]。

奇異值分解是一種非線性濾波，該方法將時間序列從矩陣角度出發(fā)，將時間序列矩陣分解到一些列奇異值對應(yīng)的時頻子空間中，從而對信號進(jìn)行了提取。

步驟1 將輸入信號X構(gòu)造成Hankel矩陣A，對矩陣A進(jìn)行奇異值分解，并提取信號Xi。

將時間序列X=[x(1),x(2),…,x(N)]構(gòu)造成Hankel矩陣A：

(4)

其中，N=m+n-1。

將Hankel矩陣A進(jìn)行奇異值分解，并表達(dá)成以下形式：

(5)

式(5)中，Σ=diag(σ1,σ2,…,σn)，σi為矩陣A的奇異值，即Σ為奇異值σi降序排序。矩陣A就可以構(gòu)造為一系列奇異值σi與子矩陣Ai乘積之和，并且子矩陣Ai的形式如下：

(6)

取矩陣Ai的第一行全部元素Xi(t),t=1,2,…，n和最后一列的后m-1個元素Xi(t),t=n+1,…,m+n-1且N=m+n-1，將其組合成新的時間序列Xi：

Xi=[xi(1),xi(2),…,xi(N)]

(7)

奇異值分解所得子信號Xi為原始信號X的線性變換，所以當(dāng)輸入信號為不存在目標(biāo)磁異常信號時，輸出信號Xi仍是高斯白噪聲。

步驟2 對Xi進(jìn)行FFT，并找出在特征頻段能量最大的分解信號F[Xi]。

奇異值分解的目的就是尋找在水下信號特征頻段處能量最大的奇異值分解信號F[Xi]，該子信號Xi即為最優(yōu)子信號。取分解階次n，不同階次n所對應(yīng)奇異值σi的大小，σi的分布呈現(xiàn)明顯的階梯形，表明噪聲與信號分離情況較好[11]。如式(8)所示，對各奇異值σi下的分解信號Xi在特征頻段的能量Ei進(jìn)行計算，積分區(qū)間為[0,1]是因為目標(biāo)磁異常信號的能量主要集中在0～1 Hz。

(8)

子信號Xi特定頻段能量Ei最大時的子信號，即為最優(yōu)子信號Xi。

1.3 隨機共振原理

隨機共振是1981年Benzi研究古冰川氣象問題時提出的概念，并在生物、電子、信號處理、故障檢測[13]方面得到了發(fā)展[14-15]。隨機共振是指在非線性系統(tǒng)的作用下，噪聲會在一定程度上加強系統(tǒng)的輸出響應(yīng)。

經(jīng)典的隨機共振模型可用如下的非線性郎之萬方程[16](Langevin Equation)表示：

(9)

其中，V(x)為勢函數(shù)，表示如下：

(10)

式(10)中，a,b為隨機共振系統(tǒng)的結(jié)構(gòu)參數(shù)，x為系統(tǒng)的輸出，S+Γ隨機共振系統(tǒng)的輸入信號，S+Γ為水下目標(biāo)磁場信號，Γ為零均值高斯白噪聲。

但當(dāng)磁異常信號出現(xiàn)時，隨機共振系統(tǒng)的輸出會出現(xiàn)波動，從而使輸出x的標(biāo)準(zhǔn)差(The Standard Deviation，STD)出現(xiàn)異常。因此，可以先依次在隨機共振系統(tǒng)的輸出信號x中截取一段長度為N時間序列，隨后用其輸出x的標(biāo)準(zhǔn)差作為檢測統(tǒng)計量Y。

對于非線性郎之萬方程，可以通過龍格庫塔法[17](Runge-Kutta Methods)進(jìn)行求解：

(11)

式(11)中，τ是采樣間隔時間，u(n)是S+Γ奇異值分解后得到的時間序列，xi(n)是隨機共振系統(tǒng)的輸出。

檢驗統(tǒng)計量Y表示式如下：

(12)

2 奇異值分解-隨機共振(SVD-RS)算法

如圖1所示，本文對信號的處理一共有兩個模塊：一是奇異值分解(SVD)模塊；二是隨機共振(RS)模塊。信號通過SVD模塊，提取特征頻段能量最大信號，隨后將其通過隨機共振模塊，并用統(tǒng)計量Y作為SVD-RS檢測器的檢驗統(tǒng)計量。

圖1 檢測流程圖Fig.1 Test flow chart of SVD-RS system

根據(jù)潛艇磁場產(chǎn)生機理，將潛艇磁矩分為縱向、橫向、垂向3個固定磁矩和縱向、橫向、垂向3個感應(yīng)磁矩共6個磁矩，分別記為mil，mit，miv，mpl，mpt，mpv。在一次探測過程中，一般不考慮地磁場變化，那么miv也可以認(rèn)為不變，因此考慮垂向磁矩mv時不分開考慮。

在仿真情況下，水下目標(biāo)的磁矩各分量取mil=5×104A·m2，mit=0.27mil，mpl=1.5mil，mpt=0.18mil，mv=0.85mil，磁矩信息來源于文獻(xiàn)[18]。默認(rèn)飛機與目標(biāo)的高度差為300 m，反潛機航速默認(rèn)為100 m/s，磁力計采樣頻率設(shè)為100 Hz。隨機共振系統(tǒng)的參數(shù)a=1,b=70，檢驗統(tǒng)計量Y的序列窗口長度為N=300。

通過磁偶極子模型仿真水下目標(biāo)信號，并疊加噪聲，圖2是輸入信噪比SNRin=-3.29 dB時，輸入總場信號和目標(biāo)磁異常信號的時域圖和頻域圖。其中輸入信噪比SNRin定義為在磁異常信號出現(xiàn)的時間段(4～6 s)的輸入信號的信噪比。

圖2 時域圖和頻域圖Fig.2 The time-domain and frequency-domain diagrams

由圖3(a)可知，一階奇異值σ1在信號X中所占權(quán)重最大，且由圖3(b)可知，一階分解信號X在特征頻段上的能量E1最強，故將提取X1作為輸入信號X在奇異值分解系統(tǒng)的輸出信號。圖3(c)是信號X1的時域圖，輸入信號X在奇異值分解后所提取的信號X1在4 s處出現(xiàn)波谷，在6 s處出現(xiàn)波峰，基本還原了水下目標(biāo)磁場信號S的時域特征。圖3(d)為X1的頻域圖，奇異值分解明顯對信號非特征頻段的能量進(jìn)行了有效的抑制，從而有效地提取了目標(biāo)信號。

將最優(yōu)子信號Xi與磁異常信號進(jìn)行互相關(guān)運算，結(jié)果如圖4所示。當(dāng)時延為0 s時，互相關(guān)幅值最大，在其他時延情況下，互相關(guān)函數(shù)的值明顯低于0 s的值，有效驗證了奇異值分解對信號波形提取能力。

如圖5所示，對輸入信號和加性白噪聲進(jìn)行奇異值分解，取一階信號,其輸出信噪比SNRout為1.07 dB?？梢钥闯?，奇異值分解對信噪比有明顯的增益效果，且奇異值分解的方法能夠有效地從特征頻段中提取信號，從而恢復(fù)波形。

圖3 奇異值分解的結(jié)果Fig.3 The result of singular value decomposition

圖4 互相關(guān)函數(shù)Fig.4 Cross correlation function

圖5 最優(yōu)子信號X1與背景磁場Fig.5 Optimal subsignals X1 and background magnetic field

一般來說，我們選擇系統(tǒng)參數(shù)a=b=1，系統(tǒng)輸出初值xi(1)=0。但為了系統(tǒng)更好的收斂效率，通過對參數(shù)的多次調(diào)整，取a=1,b=70。圖6(a)為奇異值分解子信號Xi，同時也是隨機共振系統(tǒng)的輸入信號S+Γ。圖6(b)是隨機共振系統(tǒng)的輸出x的波形圖，可以看出，在目標(biāo)信號出現(xiàn)的時間段4～6 s附近，輸出x存在明顯波動。圖6(c)是檢驗統(tǒng)計量Y歸一化后的波形圖，截取窗長度N=300，以標(biāo)準(zhǔn)差作為檢驗統(tǒng)計量可以較好體現(xiàn)出是否出現(xiàn)目標(biāo)。

圖6 隨機共振系統(tǒng)的響應(yīng)分析Fig.6 Response analysis of stochastic resonance systems

3 仿真實驗

為了驗證本算法的性能,分別對不同噪聲強度D下磁異常信號的檢測進(jìn)行的仿真分析，并對無目標(biāo)情況下的算法的輸出進(jìn)行了研究。圖7對奇異值分解前后的輸入信噪比SNRin、輸出信噪比SNRout進(jìn)行了不同條件下的仿真，其中磁異常信號強度不變，噪聲強度D=1,2，…，5?？梢钥闯?，在低輸入信噪比SNRin條件下，奇異值分解仍然能夠較好地提高輸出信噪比SNRout。并且當(dāng)輸入信號中不存在目標(biāo)磁異常信號時，即輸入信號僅為強度為D的高斯白噪聲信號。

圖7 信噪比變化圖Fig.7 Variation of SNR

圖8為在不同噪聲強度D下，經(jīng)過奇異值分解的最優(yōu)子信號X1與磁異常信號進(jìn)行互相關(guān)運算的結(jié)果?？梢钥闯?，互相關(guān)函數(shù)在時延為0 s達(dá)到峰值，即最優(yōu)子信號X1基本恢復(fù)了信號時域特征。無目標(biāo)時，輸入為高斯白噪聲，經(jīng)過線性的奇異值分解，輸出仍是高斯白噪聲，并且與原輸入高斯信號的互相關(guān)函數(shù)值為在0 s處的輸出幅度較低。

圖8 互相關(guān)函數(shù)圖Fig.8 Cross correlation function

圖9為隨機共振系統(tǒng)的檢驗統(tǒng)計量Y歸一化后在不同時間的值，其中窗長度為N，隨機共振系統(tǒng)結(jié)構(gòu)參數(shù)a=1,b=70?？梢钥闯?，在不同噪聲強度下，X1通過隨機共振系統(tǒng)后的檢驗統(tǒng)計量能夠較好的檢測出信號。在6 s處，存在信號的同時，Y也達(dá)到了最大值；當(dāng)不存在信號時，檢驗統(tǒng)計量Y僅為0.3，明顯可以檢測出不存在信號。

圖9 檢驗統(tǒng)計量Fig.9 Response analysis of stochastic resonance systems

4 結(jié)論

本文提出了磁異常信號奇異值分解的隨機共振檢測方法。該方法將奇異值分解子信號通過隨機共振系統(tǒng)，從強噪聲中較好地恢復(fù)了水下目標(biāo)信號的波形，并建立了有效的檢測統(tǒng)計量。仿真實驗結(jié)果表明，輸入信號在信噪比僅為-12 dB時，經(jīng)過奇異值分解，仍能較好地恢復(fù)波形，并極大地提高了信噪；隨后通過隨機共振系統(tǒng)后，亦能較好地檢測信號。