基于HHT的艦船地震波信號實時檢測算法

孫慶鵬，周海文

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

0 引言

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

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

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

1 艦船地震波信號HHT分析

1.1 HHT信號分析原理

HHT分析是由經驗模態分解（EMD）和Hilbert變換組成。EMD方法能夠自適應地將復雜信號分解成一系列具有不同時間尺度的IMF分量和1個殘余分量，從高頻到低頻的順序依次分解，如下式[16]：

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

對EMD分解得到的IMFi分量進行Hilbert變換，可以在不造成信息缺失的前提下，得到一個復IMF分量信號，即

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

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

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

邊際譜提供了對每個頻率的幅值測量，表達了整個時間段內幅值的累積。

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

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

Hilbert能量譜提供了對每個頻率的能量大小，表示的是每段頻率再整個時間長度內所累計的能量。

1.2 實測艦船地震信號HHT分析

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

圖1 試驗場景圖Fig.1 Experimental environment

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

圖2 實測貨船地震波數據Fig.2 Measured seismic wave data of cargo ship

圖3為信號經EMD分解后得到的13個IMF分量和1個殘余項（Res）。從圖3中可以看出，隨著EMD分解的進行，本征固有模態函數IMF的頻率從高到低，并且幅值越來越小。最后1項為分解殘余項，該部分表征信號總體變化趨勢，可以通過消除該項進行零點校準和溫度校正。

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

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

圖4 信號HHT時頻圖Fig.4 Time-frequency diagram of HHT signal

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

圖5 信號功率譜Fig.5 Signal power spectrum

圖6 信號邊際譜Fig.6 Signal marginal spectrum

圖7 信號瞬時能量譜Fig.7 Signal instantaneous energy spectrum

圖8 信號Hilbert能量譜Fig.8 Signal Hilbert energy spectrum

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

2 基于Hilbert能量譜的實時檢測算法

2.1 構造檢測特征量

由前面的實測艦船地震波面HHT分析可知，艦船地震波信號的能量集中于10～30 Hz之間，因此，可將信號的10～30 Hz頻段內的能量和作為檢測特征量，對艦船目標進行檢測。檢測特征量計算公式如下：

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

2.2 實時檢測算法

為實現對艦船地震波信號的實時檢測，基于上述構造的特征量，利用滑動窗進行數據截取和檢測特征量更新。檢測時，設置檢測閾值，當特征量連續多次大于檢測閾值時，判斷有信號出現。算法的原理框如圖9所示。

圖9 地震波信號實時檢測算法Fig.9 Detection algorithm for seismic wave signal

算法的具體步驟如下：

1）信號預處理，利用帶通濾波器對信號進行濾波，帶通濾波器的截止頻率為10 Hz和30 Hz。

2）Hilbert能量譜計算。利用滑動窗W截取一段數據，對該段數據求解Hilbert能量譜：

式中，H(w,t)為信號的HHT時頻譜。

3）特征提取。提取w1～w2頻段內的能量，按照式（8）計算i時刻特征量。

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

5）目標檢測。從初始時刻開始，若從 i時刻起連續 3次特征量大于閾值，即當iT＞Ui-1、Ti+1＞Ui、且Ti+2＞Ui+1時，則判定存在目標信號。

3 算例

3.1 實測算例

根據上述的地震波信號實時檢測算法，對2種貨船實測的地震波數據進行試驗，該數據是在前面的同樣試驗條件下獲取的。

檢測試驗時按照算法步驟進行，其中數據滑動窗的窗長N為1 024、步長L為460，信號的采樣率Fs為460 Hz，即每秒更新1次特征量，閾值調整窗口K為30，閾值調整參數μ為1.5。

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

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

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

從上述結果可知，本文所提出的檢測算法能夠有效檢測地震波信號，且其檢測時刻點都在目標中心經過正橫之前，具有較好的實時性和一定的檢測距離。

3.2 仿真算例

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

第1節中截取的Z方向上的數據長度為21 s，目標中心通過時刻大約在 10 s處。現假設環境背景噪聲的長度為180 s，將目標信號每隔10 s疊加到其中，可疊加6段目標信號。環境背景可用高斯白噪聲代替，其大小根據設定的信噪比SNR來計算，信噪比SNR的計算公式為10log（信號時域峰-峰值/環境背景噪聲峰-峰值）。

仿真時，將信噪比設置為-6 dB，如圖12為疊加前后的地震波信號，從圖中可以看出，目標信號完全被淹沒在背景噪聲中。下面對仿真信號進行試驗（檢測參數按照3.1節中設置）。

圖12 仿真信號Fig.12 Simulated signal

利用本文提出的檢測算法進行檢測，檢測結果如圖13所示。圖中紅色實體“□”為檢測點，從圖中可以看出，在目標出現處特征量都有明顯的變化，且算法能夠檢測到目標信號。

圖13 仿真信號檢測結果Fig.13 Detection results of simulated signal

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

表1 仿真信號檢測時間Table 1 Detection time of simulated signal

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

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

4 結束語

本文針對艦船地震波信號的非平穩特點，利用非平穩處理方法HHT對信號進行分析，得到了艦船地震波信號的 Hilbert能量譜，并基于信號的Hilbert能量譜特點，提出基于滑動小波能量譜的艦船地震波檢測算法，重點討論了算法的檢測效果及信噪比對檢測算法性能的影響。基于實測和仿真數據試驗表明：本文所提出的地震波信號實時檢測算法能夠有效地實時對目標信號進行檢測；在正常環境中檢測距離大約為60 m，在SNR為-6 dB時檢測距離大約在10 m左右，信噪比越小檢測距離越近。從不同信噪比仿真信號檢測的結果中可以得到，在信噪比為-15 dB時，檢測概率達到96%。

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