孫華慶 王 丹

（海軍潛艇學(xué)院 青島 266000）

1 引言

磁探測(cè)系統(tǒng)是對(duì)潛艇聲納探測(cè)的有效補(bǔ)充，可實(shí)現(xiàn)對(duì)潛艇目標(biāo)的高精度定位。針對(duì)復(fù)雜信號(hào)常見的非平穩(wěn)性、非線性特點(diǎn)，美國(guó)華裔Huang等在1998年提出了具有創(chuàng)新性的信號(hào)處理方法——經(jīng)驗(yàn)?zāi)B(tài)分解方法（Empirical Mode Decomposition，EMD）［1～2］。

EMD算法對(duì)磁異常信號(hào)的處理基本能達(dá)到從強(qiáng)噪聲信號(hào)中提取目標(biāo)信號(hào)的目的，但是對(duì)于本身磁異常信號(hào)強(qiáng)度弱、背景噪聲信號(hào)頻帶較寬的情況提取有一定誤差。主要原因是IMF分解不正常，稱模態(tài)混疊現(xiàn)象［3］，該現(xiàn)象有兩種情況，一種是一個(gè)IMF中包含著差異極大的時(shí)間特征尺度，另外一種是相近的特征時(shí)間尺度分布在不同的IMF中，會(huì)導(dǎo)致IMF波形不能真實(shí)反映該頻率的波形或相鄰的兩個(gè)IMF存在波形混疊現(xiàn)象。這種現(xiàn)象和輸入信號(hào)特征有關(guān)，特別是非線性非平穩(wěn)的脈沖信號(hào)、間歇信號(hào)及含高斯噪聲信號(hào)。

為解決這一模態(tài)混疊現(xiàn)象，Wu等在2009年提出了總體平均經(jīng)驗(yàn)?zāi)B(tài)分解算法（Ensemble Empir?ical Mode Decomposition，EEMD）［4］。Wu利用高斯白噪聲的特點(diǎn)，在含噪信號(hào)中加入高斯白噪聲信號(hào)，對(duì)該合成信號(hào)EMD分解出IMF后，再重復(fù)向原始信號(hào)中加入不同的高斯噪聲進(jìn)行EMD分解，在多次重復(fù)加噪分解后進(jìn)行標(biāo)準(zhǔn)化平均，恢復(fù)原數(shù)據(jù)分布大小。

在EEMD過(guò)程中，對(duì)原始信號(hào)加入高斯噪聲能夠改變信號(hào)中極值點(diǎn)位置及幅值，從而在IMF分解前改變上下包絡(luò)線的形態(tài)，在提取的低階分量中可以把原有信號(hào)噪聲、添加的高斯噪聲和部分信號(hào)中可能存在的間歇信號(hào)分離出來(lái)［5～6］。

對(duì)經(jīng)驗(yàn)?zāi)B(tài)分解算法的分析中，CEEMD相對(duì)于另外兩種更加優(yōu)化，但是在IMF分量重構(gòu)的過(guò)程中，需要考慮如何確定濾除低頻分量的階數(shù)。濾掉的分量多，重構(gòu)信號(hào)平滑性好，沒有噪聲的干擾，但目標(biāo)信號(hào)有遺失；濾掉分量少，信號(hào)保存較完整，但是噪聲干擾太大目標(biāo)信號(hào)特征無(wú)法提取。在分量選擇過(guò)程中，需要對(duì)兩個(gè)指標(biāo)進(jìn)行平衡優(yōu)化，根據(jù)原始信號(hào)的特點(diǎn)及選擇最合適的分解階數(shù)。

2 潛艇磁異常信號(hào)仿真

當(dāng)磁探儀與目標(biāo)位置點(diǎn)距離大于2.5倍的磁性目標(biāo)長(zhǎng)度時(shí)，該目標(biāo)可以視為一個(gè)磁偶極子［7］。距離磁性目標(biāo)r(x,y,z)處的磁場(chǎng)可表示為

式中：m(mx,my,mz)為潛艇的磁矩，μ0為真空磁導(dǎo)率，r= ||r為潛艇到磁探儀的距離［8］。

仿真假設(shè)條件如下，在水深200m的某試驗(yàn)海區(qū)，在地磁坐標(biāo)系中，設(shè)試驗(yàn)海域?yàn)榘霃?00m的圓形，潛艇位于海域中心，速度為4m/s。距離海底50m。設(shè)某潛艇磁矩

圖1 原始目標(biāo)信號(hào)仿真圖

在上述原始目標(biāo)信號(hào)的基礎(chǔ)上添加信噪比-20dB高斯白噪聲，這種情況下，潛艇磁異常信號(hào)基本淹沒在背景噪聲中，無(wú)法從時(shí)域信號(hào)中分辨出信號(hào)的特征甚至不能發(fā)現(xiàn)潛艇信號(hào)存在。

圖2 加噪信號(hào)仿真圖

3 EEMD算法信號(hào)處理

信號(hào)迭代過(guò)程中，計(jì)算前后數(shù)據(jù)的標(biāo)準(zhǔn)偏差（Standard Deviation，SD），若標(biāo)準(zhǔn)偏差小于某一預(yù)設(shè)閾值，則停止。信號(hào)分解表達(dá)式為

圖3是對(duì)加噪信號(hào)進(jìn)行EEMD分解的結(jié)果。從結(jié)果可以看到，共進(jìn)行了11階的分解，得到11個(gè)IMF分量。第一階IMF1和第二階IMF2基本全是高頻的高斯噪聲，信號(hào)信息基本分辨不出；從第三階IMF3開始，磁異常信號(hào)逐漸顯示出來(lái)，特別是到了第五階IMF5，噪聲基本消失，磁異常信號(hào)在該層的分離膠明顯，和原始信號(hào)曲線的變化趨勢(shì)類似；到第十一階IMF11，信號(hào)逐漸失真，向單調(diào)曲線變化，說(shuō)明磁異常信號(hào)在這個(gè)階層的頻率上很少；最后是余項(xiàng)Res，余項(xiàng)基本成為一條單調(diào)曲線，無(wú)法獲取信號(hào)的信息。

由EEMD的IMF分解結(jié)果可知，如在IMF3中，EEMD的極值點(diǎn)分布更均勻，說(shuō)明EEMD算法達(dá)到了均化極值點(diǎn)分布的目的，經(jīng)過(guò)數(shù)次平均，減少添加的噪聲的影響，從而減小了模態(tài)混疊效應(yīng)。

4 最優(yōu)降噪算法推導(dǎo)

青島理工大學(xué)的鄭一教授針對(duì)隨鉆鉆井液脈沖信號(hào)容易受噪聲干擾的特點(diǎn)，利用EEMD算法將信號(hào)分解不同的固有模態(tài)函數(shù)，進(jìn)行降噪濾波，對(duì)脈沖信號(hào)方波整形處理，引入“影響因子”的概念，綜合考慮重構(gòu)信號(hào)的逼近度和相關(guān)度。該方法對(duì)鉆井液脈沖信號(hào)的處理的結(jié)果合理有效［9］。

基于EEMD的最優(yōu)降噪算法，針對(duì)受高斯噪聲污染的潛艇磁異常信號(hào)，提出建立EEMD最優(yōu)降噪算法濾波器，即將高低頻的分量進(jìn)行按頻率分離。使用兩個(gè)指標(biāo)進(jìn)行選擇最佳降噪階數(shù)。

圖3 EEMD分解IMF分量及Rs余項(xiàng)圖

1）相似性指標(biāo)

設(shè)磁異常信號(hào)為，將信號(hào)進(jìn)行CEEMD分解，構(gòu)建低通濾波器。取得n層固有模態(tài)分量IMF及最后一個(gè)剩余分量residual之和，IMF1～I(xiàn)MF3之和……，即

則Ex1,Ex2,Ex3,…,Exn是分離重構(gòu)后的信號(hào)，該信號(hào)中依次頻率降低，含噪比例依次降低，目標(biāo)信號(hào)成分比例依次降低。因此需要找到一個(gè)合適的分離降噪層數(shù)，既能保證信號(hào)的相對(duì)完整性，又能盡可能多地濾掉高頻信號(hào)。

為表征信號(hào)的相似性，用標(biāo)準(zhǔn)差（均方誤差）進(jìn)行衡量。將n個(gè)處理過(guò)的去噪信號(hào)與原始信號(hào)分別作差，得n個(gè)信號(hào)的殘差，記為C1,C2,C3,…,Cn，對(duì)殘差進(jìn)行均方誤差分析?？紤]便于不同信號(hào)特征比較，統(tǒng)一標(biāo)準(zhǔn)，使其殘差標(biāo)準(zhǔn)差指標(biāo)范圍在0～1之間。

殘差標(biāo)準(zhǔn)差矩陣為

則歸一化后為

2）光滑性指標(biāo)

為描述曲線的光滑性指標(biāo)，將任意一段曲線在x0點(diǎn)分兩左右段，設(shè)為P（t），Q（t），在一定步長(zhǎng)的小區(qū)間內(nèi)，P，Q的曲率越接近，說(shuō)明該段的平滑性越好；如果很多一定步長(zhǎng)連續(xù)的小區(qū)間的平滑性都較好，說(shuō)明曲線整體的光滑性好。在x0點(diǎn)左右小區(qū)間段的曲率為

設(shè)區(qū)域取值步長(zhǎng)為l，則

在x0點(diǎn)位置，定義曲線光滑度的指標(biāo)計(jì)算公式：

其中，對(duì)于每一條曲線

由光滑度定義可知，SN越小越接近0，說(shuō)明曲線越光滑，當(dāng)SN=0時(shí)，說(shuō)明在x0點(diǎn)小區(qū)域內(nèi)曲線為直線。為分析整個(gè)曲線的光滑性，設(shè)定某一間隔取x0的值，將整個(gè)曲線用多個(gè)連續(xù)的小區(qū)間表達(dá)出來(lái)，每個(gè)小區(qū)間的光滑度值集中到該曲線對(duì)應(yīng)的矩陣SN中。

對(duì)每條曲線的光滑度矩陣求標(biāo)準(zhǔn)差，光滑度標(biāo)準(zhǔn)差組成的矩陣為

同理，為方便比較，使其光滑度標(biāo)準(zhǔn)差指標(biāo)范圍在0～1之間。

建立目標(biāo)函數(shù)FM，引入權(quán)重參數(shù)α(0＜α＜1)，根據(jù)目標(biāo)信號(hào)特征平衡光滑性和相似性的矛盾，進(jìn)而選擇最佳的降噪層數(shù)，實(shí)現(xiàn)復(fù)雜磁異常信號(hào)的提取。α(0＜α＜0.5)時(shí)，追求降噪信號(hào)曲線更光滑；α(0.5＜α＜1)時(shí)，追求降噪信號(hào)與目標(biāo)真實(shí)信號(hào)相似性。

當(dāng)FM取最小值時(shí)，對(duì)應(yīng)的降噪階數(shù)即最優(yōu)。

5 信號(hào)分析

表1 優(yōu)化算法指標(biāo)表（α=0.5）

當(dāng)α=0.5時(shí)，濾波過(guò)程中信號(hào)的平滑性和相似性同等權(quán)重，從高頻到低頻，隨著信號(hào)重構(gòu)分量減少，合成后的信號(hào)特征與目標(biāo)真實(shí)信號(hào)特征相似度降低，合成信號(hào)中的噪聲含量降低，曲線更加平滑，信號(hào)整體完整性降低。FM參數(shù)先降低后升高，在IMF5時(shí)取最小值，根據(jù)上節(jié)推導(dǎo)公式，目標(biāo)函數(shù)FM取最小值0.3201時(shí)有最佳降噪階數(shù)，即濾掉前5階高頻分量，得到的濾波信號(hào)保留了目標(biāo)磁異常信號(hào)的特征，同時(shí)平滑性較好，沒有噪聲的干擾。

利用EEMD算法對(duì)不同階數(shù)的信號(hào)進(jìn)行重構(gòu)合成，結(jié)果如下：

圖4 不同階數(shù)EEMD信號(hào)累加重構(gòu)圖

圖4（a）～（d）分別代表4階降噪到7階降噪結(jié)果曲線與原始信號(hào)對(duì)比圖。由圖可知4階降噪結(jié)果異常小突起過(guò)多，在進(jìn)行磁異常信號(hào)檢測(cè)時(shí)容易造成誤判；第5階降噪后曲線能夠反映原始目標(biāo)信號(hào)變化趨勢(shì)，兩側(cè)曲線波動(dòng)較小，不會(huì)對(duì)異常提取造成干擾噪聲基本消失，原始目標(biāo)信號(hào)特征顯現(xiàn)出來(lái)，可以進(jìn)行特征提?。坏?階降噪過(guò)于平滑，曲線無(wú)法表達(dá)原始目標(biāo)信號(hào)的特征，目標(biāo)信號(hào)信息丟失過(guò)多，第7階及后重構(gòu)曲線基本為一條直線，不具有借鑒意義。

利用平滑性和相似性兩個(gè)指標(biāo)建立目標(biāo)函數(shù)，對(duì)EEMD算法的降噪階數(shù)進(jìn)行優(yōu)化選擇，仿真結(jié)果表明該最優(yōu)降噪算法原理簡(jiǎn)單，階數(shù)選擇合適，具有一定穩(wěn)定性。

6 結(jié)語(yǔ)

本文仿真分析了潛艇磁異常經(jīng)典模型——磁偶極子模型，在潛艇磁異常信號(hào)中加入高強(qiáng)度噪聲，利用總體平均經(jīng)驗(yàn)?zāi)B(tài)分解算法對(duì)該信號(hào)處理。IMF分解過(guò)程中發(fā)現(xiàn)，噪聲大部分集中在前兩階，對(duì)分解后的各層信號(hào)進(jìn)行重構(gòu)合成，結(jié)果基本能將潛艇磁異常特征表現(xiàn)出來(lái)。成對(duì)的添加高斯噪聲，在抑制模態(tài)混疊的同時(shí)減小了增加噪聲的影響。建立指標(biāo)目標(biāo)函數(shù)FM，綜合考慮相似性和平滑性的基礎(chǔ)上選擇最優(yōu)階數(shù)?；贓EMD的最優(yōu)降噪濾波算法對(duì)仿真數(shù)據(jù)的處理中效果較好。

本文結(jié)果表明，EEMD算法在處理強(qiáng)噪聲干擾的磁異常的信號(hào)時(shí)，能夠較好地提取目標(biāo)信號(hào)特征。本文的結(jié)果對(duì)今后水下潛艇磁異常探測(cè)信號(hào)檢測(cè)算法的發(fā)展有一定借鑒意義。

參考文獻(xiàn)

［1］HUANG N E，SHEN Z，LONG S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlin?ear and non—stationary time series analysis［C］//Proc.Roy.Soc.London，1998（454）：903-995.

［2］HUANG N E，SHEN Z，LONG R S.A new view of nonlin?ear water waves-the Hilbert spectrum，Ann.Rev.Fluid Mech，1999（31）：417-457.

［3］GAI G H.The processing of rotor startup slgnals based on empirical mode decomposition［J］.Mechanical Systems and Signal Processing，2006（20）：225-235.

［4］WU Z H，HUANG N E.Ensemble empirical mode decom?position：A noise assisted data analysis method［J］.Ad?vances in Adaptive Data Analysis，2009（1）：1-41.

［5］YEH J R，SHIEH J S.Complementary ensemble empiri?cal mode decomposition：A noise enhanced data analysis method［J］，Advances in Adaptive Data Analysis，2010，2（2）：135-156.

［6］鄭近德，程軍圣，楊宇.改進(jìn)的EEMD算法及其應(yīng)用研究［J］，振動(dòng)與沖擊，2013，32（21）：21-26.

［7］張朝陽(yáng)，肖昌漢，高俊吉，等.磁性物體磁偶極子模型試用性的試驗(yàn)研究［J］.應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào)，2010，18（5）：862-868.

［8］單志超，曲曉慧，楊日杰，等.潛艇航向?qū)χ鄙龣C(jī)磁異探潛的影響［J］. 火力與指揮控制，2013，38（2）：62-64.

［9］鄭一，孫曉峰，陳健，岳軍.基于集合經(jīng)驗(yàn)?zāi)B(tài)分解的隨鉆脈沖信號(hào)優(yōu)良降噪整形算法［J］.石油勘探與開發(fā)，2012，39（6）：750-753.