基于HHT的艦船地震波信號(hào)實(shí)時(shí)檢測算法

孫慶鵬，周海文

（1.中國人民解放軍92213部隊(duì)，廣東 湛江 524064；2.海軍工程大學(xué) 電器工程學(xué)院，湖北 武漢 430033）

0 引言

隨著反水雷技術(shù)的不斷進(jìn)步，水雷引信中常用的聲、磁、水壓等物理場信號(hào)不斷被抑制，水雷引信對目標(biāo)艦船的檢測識(shí)別變得越來越困難。艦船海底地震波是近些年發(fā)現(xiàn)的新物理場，其是艦船在淺海航行時(shí)引起的振動(dòng)、噪聲和水體擾動(dòng)以壓力波的形式經(jīng)流體介質(zhì)傳播到海底，進(jìn)而引起海底介質(zhì)的振動(dòng)，在海底界面和海底以彈性波的形式向周圍傳播[1]。艦船地震波信號(hào)的成分集中于低頻段，其具有衰減小、傳播距離遠(yuǎn)、震動(dòng)幅度大等特點(diǎn)，可以作為水雷引信的補(bǔ)充[2]。

目前，俄、美兩國已經(jīng)在水雷引信系統(tǒng)中使用艦船地震波場聯(lián)合其他物理場對目標(biāo)進(jìn)行檢測和識(shí)別[3-4]，但出于保密，其相關(guān)技術(shù)并未報(bào)道。國內(nèi)的研究主要集中于解釋艦船海底地震波形成機(jī)理和傳播規(guī)律上[5-8]，著重分析了艦船地震波波動(dòng)成分，并研究了其在不同海底介質(zhì)中的傳播特性。對于艦船地震波信號(hào)的應(yīng)用，國內(nèi)只初步探討了其可用于艦船目標(biāo)檢測和識(shí)別的可能性，例如：文獻(xiàn)[9]初步探討了艦船地震波信號(hào)用于艦船目標(biāo)檢測的可行性；文獻(xiàn)[10]指出艦船地震波信號(hào)能夠進(jìn)行遠(yuǎn)距離的檢測；沈陽理工大學(xué)將艦船等效為固定的點(diǎn)震源，利用艦船地震波中的Scholte波特點(diǎn)對艦船目標(biāo)進(jìn)行探測和識(shí)別[11-12]。

艦船地震波信號(hào)是一種非平穩(wěn)信號(hào)，文獻(xiàn)[13]基于短時(shí)傅里葉變換提出一種實(shí)時(shí)的艦船檢測算法，但短時(shí)傅里葉變換的時(shí)域和頻域分辨率固定，而且其所用檢測閾值是固定的，檢測效果受環(huán)境影響較大；文獻(xiàn)[14]中首次運(yùn)用了小波變換對艦船地震波信號(hào)進(jìn)行分析，但只探討了小波變化在特征提取中的優(yōu)勢，并未研究具體的艦船地震波信號(hào)檢測算法，并且小波變換也是以傅里葉變換為基礎(chǔ)，其小波基的選擇也會(huì)影響信號(hào)分析的結(jié)果。HHT具有簡單高效、自適應(yīng)性強(qiáng)、時(shí)頻分辨率高等優(yōu)點(diǎn)，更適用于處理非線性、非平穩(wěn)信號(hào)[15]。本文將利用HHT方法對實(shí)測艦船地震波信號(hào)進(jìn)行分析，然后提出基于 Hilbert能量譜的艦船地震波信號(hào)實(shí)時(shí)檢測算法。

1 艦船地震波信號(hào)HHT分析

1.1 HHT信號(hào)分析原理

HHT分析是由經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）和Hilbert變換組成。EMD方法能夠自適應(yīng)地將復(fù)雜信號(hào)分解成一系列具有不同時(shí)間尺度的IMF分量和1個(gè)殘余分量，從高頻到低頻的順序依次分解，如下式[16]：

式中：IMFi為EMD分解得到的一些本征模態(tài)分量；rn(t)為殘余分量。

對EMD分解得到的IMFi分量進(jìn)行Hilbert變換，可以在不造成信息缺失的前提下，得到一個(gè)復(fù)IMF分量信號(hào)，即

式中：ai(t)為復(fù)IMF分量的瞬時(shí)幅值；θi(t)為瞬時(shí)相位。根據(jù)瞬時(shí)頻率定義可得瞬時(shí)頻率：

因此，可以得到希爾伯特時(shí)頻譜，也稱Hilbert譜，記為

如果將 H(w,t)對時(shí)間積分，就能得到 Hilbert邊際譜，即

邊際譜提供了對每個(gè)頻率的幅值測量，表達(dá)了整個(gè)時(shí)間段內(nèi)幅值的累積。

將式（4）幅值取平方，對其頻率積分，可得Hilbert瞬時(shí)能量為

Hilbert瞬時(shí)能量提供了信號(hào)能量隨時(shí)間變化的情況。將幅值的平方對時(shí)間積分，便可得Hilbert能量譜：

Hilbert能量譜提供了對每個(gè)頻率的能量大小，表示的是每段頻率再整個(gè)時(shí)間長度內(nèi)所累計(jì)的能量。

1.2 實(shí)測艦船地震信號(hào)HHT分析

根據(jù)上述HHT分析原理，現(xiàn)對一組實(shí)測的貨船地震波信號(hào)進(jìn)行HHT分析。該數(shù)據(jù)是在漢江某處由地震波采集系統(tǒng)（采集加速度矢量）測量獲得，試驗(yàn)時(shí)采樣頻率為 460 Hz，周圍無明顯的工業(yè)噪聲污染，測量目標(biāo)為江上正在行進(jìn)的貨船，航速大約為10 kn，沿水流方向直線行駛。圖1為當(dāng)時(shí)實(shí)測的試驗(yàn)場景。目標(biāo)A/B船的噸位大約為2 000 t，航速分別為10 kn、13 kn，采集系統(tǒng)放在沙地中，與目標(biāo)的正橫距離大約為50 m。

圖1 試驗(yàn)場景圖Fig.1 Experimental environment

文獻(xiàn)[10]中對艦船地震波信號(hào)的時(shí)頻進(jìn)行分析，指出艦船地震波信號(hào)在垂直方向上的檢測距離最遠(yuǎn)，頻域線譜衰減時(shí)間長。圖2為實(shí)測的貨船垂直軸方向的地震波信號(hào)，從圖中可知當(dāng)艦船經(jīng)過測量系統(tǒng)時(shí)，加速度有明顯的幅度變化。同時(shí)，考慮到數(shù)據(jù)較多，在這截取垂直軸 az方向的一部分加速度數(shù)據(jù)進(jìn)行HHT分析，如圖2（a）中的紅色框中的一段數(shù)據(jù)，進(jìn)行HHT分析。

圖2 實(shí)測貨船地震波數(shù)據(jù)Fig.2 Measured seismic wave data of cargo ship

圖3為信號(hào)經(jīng)EMD分解后得到的13個(gè)IMF分量和1個(gè)殘余項(xiàng)（Res）。從圖3中可以看出，隨著EMD分解的進(jìn)行，本征固有模態(tài)函數(shù)IMF的頻率從高到低，并且幅值越來越小。最后1項(xiàng)為分解殘余項(xiàng)，該部分表征信號(hào)總體變化趨勢，可以通過消除該項(xiàng)進(jìn)行零點(diǎn)校準(zhǔn)和溫度校正。

圖3 貨船地震波信號(hào)EMD分解Fig.3 EMD decomposition of cargo ship seismic wave signal

圖4為信號(hào)的HHT信號(hào)時(shí)頻圖，該圖很好地表達(dá)了信號(hào)幅值隨時(shí)間和頻率的變化過程，用HHT時(shí)頻圖可以很好地完成對突變信號(hào)的檢測，如圖4中所示，突變信號(hào)發(fā)生在第20 s，突變信號(hào)的頻率集中在10～30 Hz之間。

圖4 信號(hào)HHT時(shí)頻圖Fig.4 Time-frequency diagram of HHT signal

圖5為信號(hào)利用周期圖法得到的信號(hào)功率譜。圖 6-8分別為根據(jù)公式（5）-（7）得到的地震波信號(hào)的邊際譜、瞬時(shí)能量譜及 Hilbert能量譜。對比圖5和圖6可知，信號(hào)的功率譜往往會(huì)低估低頻部分的幅值；瞬時(shí)能量譜表達(dá)了信號(hào)能量隨時(shí)間變化情況，從圖7中可以明顯看出在20 s左右信號(hào)能量發(fā)生變化，這與信號(hào)（圖2（b））動(dòng)態(tài)變化相一致；Hilbert能量譜反映信號(hào)每個(gè)頻段在整個(gè)時(shí)間長度內(nèi)所累積的總能量，從圖 8中可知，信號(hào)能量集中在10～30 Hz之間，在100～200 Hz頻段內(nèi)存在一些線譜，這是由于信號(hào)中含有目標(biāo)水中傳播的直達(dá)噪聲信號(hào)，對比圖5和圖8，可以看出 Hilbert能量譜能夠更好地反映信號(hào)高頻段的特征。

圖5 信號(hào)功率譜Fig.5 Signal power spectrum

圖6 信號(hào)邊際譜Fig.6 Signal marginal spectrum

圖7 信號(hào)瞬時(shí)能量譜Fig.7 Signal instantaneous energy spectrum

圖8 信號(hào)Hilbert能量譜Fig.8 Signal Hilbert energy spectrum

從功率譜、邊際譜以及信號(hào)的 Hilbert能量譜都能清楚地看出，信號(hào)的能量集中在10～30 Hz之間，該明顯的頻域特征能夠用于艦船目標(biāo)的地震波信號(hào)檢測。

2 基于Hilbert能量譜的實(shí)時(shí)檢測算法

2.1 構(gòu)造檢測特征量

由前面的實(shí)測艦船地震波面HHT分析可知，艦船地震波信號(hào)的能量集中于10～30 Hz之間，因此，可將信號(hào)的10～30 Hz頻段內(nèi)的能量和作為檢測特征量，對艦船目標(biāo)進(jìn)行檢測。檢測特征量計(jì)算公式如下：

式中：Ei(a)為 i時(shí)刻的 Hilbert能量譜；w1和 w2為選取的頻段的下限截止頻率和上限截止頻率，這里可取w1為10 Hz，w2為30 Hz。

2.2 實(shí)時(shí)檢測算法

為實(shí)現(xiàn)對艦船地震波信號(hào)的實(shí)時(shí)檢測，基于上述構(gòu)造的特征量，利用滑動(dòng)窗進(jìn)行數(shù)據(jù)截取和檢測特征量更新。檢測時(shí)，設(shè)置檢測閾值，當(dāng)特征量連續(xù)多次大于檢測閾值時(shí)，判斷有信號(hào)出現(xiàn)。算法的原理框如圖9所示。

圖9 地震波信號(hào)實(shí)時(shí)檢測算法Fig.9 Detection algorithm for seismic wave signal

算法的具體步驟如下：

1）信號(hào)預(yù)處理，利用帶通濾波器對信號(hào)進(jìn)行濾波，帶通濾波器的截止頻率為10 Hz和30 Hz。

2）Hilbert能量譜計(jì)算。利用滑動(dòng)窗W截取一段數(shù)據(jù)，對該段數(shù)據(jù)求解Hilbert能量譜：

式中，H(w,t)為信號(hào)的HHT時(shí)頻譜。

3）特征提取。提取w1～w2頻段內(nèi)的能量，按照式（8）計(jì)算i時(shí)刻特征量。

4）自適應(yīng)閾值確定。從初始時(shí)刻開始，選取本時(shí)刻 i之前的K個(gè)特征量平均值作為基閾值U0(i)，后期可根據(jù)實(shí)際情況利用閾值調(diào)整因子μ（＞1）進(jìn)行調(diào)整，即

5）目標(biāo)檢測。從初始時(shí)刻開始，若從 i時(shí)刻起連續(xù) 3次特征量大于閾值，即當(dāng)iT＞Ui-1、Ti+1＞Ui、且Ti+2＞Ui+1時(shí)，則判定存在目標(biāo)信號(hào)。

3 算例

3.1 實(shí)測算例

根據(jù)上述的地震波信號(hào)實(shí)時(shí)檢測算法，對2種貨船實(shí)測的地震波數(shù)據(jù)進(jìn)行試驗(yàn)，該數(shù)據(jù)是在前面的同樣試驗(yàn)條件下獲取的。

檢測試驗(yàn)時(shí)按照算法步驟進(jìn)行，其中數(shù)據(jù)滑動(dòng)窗的窗長N為1 024、步長L為460，信號(hào)的采樣率Fs為460 Hz，即每秒更新1次特征量，閾值調(diào)整窗口K為30，閾值調(diào)整參數(shù)μ為1.5。

檢測的結(jié)果如圖10-11所示，圖10為A貨船的檢測結(jié)果，圖11為B貨船檢測結(jié)果。從兩張圖中可以看出，當(dāng)目標(biāo)接近測量系統(tǒng)時(shí)，特征量均發(fā)生較為明顯的變化，兩張圖中的紅色實(shí)心“□”為信號(hào)檢測時(shí)刻點(diǎn)，豎直虛線為檢測點(diǎn)的對應(yīng)時(shí)刻，A貨船中心通過是大約在164 s左右通過，B貨船中心大約在177 s左右通過。從圖10-11可知，A貨船在 153 s時(shí)被檢測到，B貨船在 166 s時(shí)被檢測到，目標(biāo)被檢測到均在目標(biāo)中心通過之前完成，且A貨船在經(jīng)過前12 s完成，B貨船在前11 s完成，試驗(yàn)時(shí)貨船的速度大約為10 kn，換算獲得檢測的距離大約在 60 m，該作用距離基本能夠滿足水雷的作戰(zhàn)需求。

圖10 A貨船檢測結(jié)果Fig.10 Detection results of cargo ship A

圖11 B貨船檢測結(jié)果Fig.11 Detection results of cargo ship B

從上述結(jié)果可知，本文所提出的檢測算法能夠有效檢測地震波信號(hào)，且其檢測時(shí)刻點(diǎn)都在目標(biāo)中心經(jīng)過正橫之前，具有較好的實(shí)時(shí)性和一定的檢測距離。

3.2 仿真算例

3.1節(jié)中所用的實(shí)測數(shù)據(jù)信噪比較高，所以檢測效果較為明顯，而現(xiàn)實(shí)應(yīng)用中信噪比往往無法控制，因此還需驗(yàn)證算法在低信噪比條件下的檢測性能。由于實(shí)測時(shí)沒有獲得信噪比低的數(shù)據(jù)，現(xiàn)截取A貨船Z方向上的一段實(shí)測數(shù)據(jù)疊加背景噪聲代替低信噪比數(shù)據(jù)。

第1節(jié)中截取的Z方向上的數(shù)據(jù)長度為21 s，目標(biāo)中心通過時(shí)刻大約在 10 s處?，F(xiàn)假設(shè)環(huán)境背景噪聲的長度為180 s，將目標(biāo)信號(hào)每隔10 s疊加到其中，可疊加6段目標(biāo)信號(hào)。環(huán)境背景可用高斯白噪聲代替，其大小根據(jù)設(shè)定的信噪比SNR來計(jì)算，信噪比SNR的計(jì)算公式為10log（信號(hào)時(shí)域峰-峰值/環(huán)境背景噪聲峰-峰值）。

仿真時(shí)，將信噪比設(shè)置為-6 dB，如圖12為疊加前后的地震波信號(hào)，從圖中可以看出，目標(biāo)信號(hào)完全被淹沒在背景噪聲中。下面對仿真信號(hào)進(jìn)行試驗(yàn)（檢測參數(shù)按照3.1節(jié)中設(shè)置）。

圖12 仿真信號(hào)Fig.12 Simulated signal

利用本文提出的檢測算法進(jìn)行檢測，檢測結(jié)果如圖13所示。圖中紅色實(shí)體“□”為檢測點(diǎn)，從圖中可以看出，在目標(biāo)出現(xiàn)處特征量都有明顯的變化，且算法能夠檢測到目標(biāo)信號(hào)。

圖13 仿真信號(hào)檢測結(jié)果Fig.13 Detection results of simulated signal

表1中列出了檢測點(diǎn)和目標(biāo)中心通過時(shí)刻。從表1中可知，在目標(biāo)中心通過之前就能夠完成對目標(biāo)的檢測，但是檢測時(shí)刻都在目標(biāo)經(jīng)過前2 s左右，換算為距離大約是 10 m左右，即檢測距離約為10 m。對比 3.1節(jié)中的檢測結(jié)果，可以得出 SNR的變化會(huì)影響本算法的檢測距離，SNR越小檢測距離越近。

表1 仿真信號(hào)檢測時(shí)間Table 1 Detection time of simulated signal

為計(jì)算檢測算法在低信噪比情況下的檢測概率，現(xiàn)按表2中信噪比生成仿真數(shù)據(jù)，各信噪比數(shù)據(jù)段中存在100段目標(biāo)信號(hào)，檢測結(jié)果如表2所示。從表中可知，在信噪比為-20 dB時(shí)，信號(hào)的檢測概率大約為86%，虛警概率大約為6%，若信噪比大于-10 dB，檢測概率為100%。因此，由仿真結(jié)果可以判斷，本文提出的檢測地震波信號(hào)檢測算法在較低的信噪比環(huán)境中仍具有較高的檢測概率。

表2 不同信噪比下的檢測結(jié)果Table 2 Detection results of different SNR

4 結(jié)束語

本文針對艦船地震波信號(hào)的非平穩(wěn)特點(diǎn)，利用非平穩(wěn)處理方法HHT對信號(hào)進(jìn)行分析，得到了艦船地震波信號(hào)的 Hilbert能量譜，并基于信號(hào)的Hilbert能量譜特點(diǎn)，提出基于滑動(dòng)小波能量譜的艦船地震波檢測算法，重點(diǎn)討論了算法的檢測效果及信噪比對檢測算法性能的影響?；趯?shí)測和仿真數(shù)據(jù)試驗(yàn)表明：本文所提出的地震波信號(hào)實(shí)時(shí)檢測算法能夠有效地實(shí)時(shí)對目標(biāo)信號(hào)進(jìn)行檢測；在正常環(huán)境中檢測距離大約為60 m，在SNR為-6 dB時(shí)檢測距離大約在10 m左右，信噪比越小檢測距離越近。從不同信噪比仿真信號(hào)檢測的結(jié)果中可以得到，在信噪比為-15 dB時(shí)，檢測概率達(dá)到96%。

本文在驗(yàn)證檢測算法性能和計(jì)算檢測概率時(shí)，利用的是在漢江岸邊所測數(shù)據(jù)，目標(biāo)為普通的貨船，對于應(yīng)用于水雷中的地震波引信應(yīng)該考慮接受的是艦船海底地震波信號(hào)。因此，為進(jìn)一步驗(yàn)證算法，接下來將進(jìn)行海底艦船目標(biāo)檢測試驗(yàn)。