艦船輻射噪聲廣義多尺度數學形態學特征提取與應用研究

郭 政， 趙 梅， 胡長青， 倪俊帥

(1. 中國科學院 聲學研究所東海研究站，上海 201815； 2. 中國科學院大學，北京 100049)

當前水聲目標識別技術的發展過程中，目標艦船識別一直是重點研究方向。艦船輻射噪聲中包含了目標艦船的航行狀態、自體振動及船體與水體相互作用等物理規律，因此有效提取表征這種物理規律的艦船輻射噪聲特征參數成為了值得研究的問題。

艦船輻射噪聲信號是一種非線性、非平穩的時變信號，因而線性方法在艦船輻射噪聲特征提取中存在一定局限性[1-2]。近年來許多學者討論了包括Lyapunov指數、關聯維數、分形特征、樣本熵等在內的非線性特征在艦船輻射噪聲特征提取中的應用，取得了一定成果。

章新華等[3]研究了艦船輻射噪聲的混沌特性，嘗試提取了艦船輻射噪聲的Lyapunov指數和分形維數。陳向東等[4]以相空間重構理論為基礎討論了艦船輻射噪聲的非線性，驗證了非線性方法應用于艦船輻射噪聲特征提取的可行性。劉朝暉等[5]提取了水聲目標分形特征向量并應用于水聲目標識別中。焦義民等[6]在構造了非線性頻譜字典的基礎上提出非線性譜特征，較好地反映了艦船輻射噪聲的低頻特征。李余興等[7]在集合經驗模態分解(ensemble empirical mode decomposition，EEMD)基礎上提取了艦船輻射噪聲的樣本熵。陳哲等提取了艦船輻射噪聲的多尺度排列熵，在多個尺度上描述了艦船輻射噪聲的復雜度特征。通過不同方法提取的艦船輻射噪聲非線性特征參數從各個方面反映了艦船輻射噪聲非線性、非平穩的特性，在不同艦船目標之間有良好的區分度。但多數非線性特征計算仍較為復雜，因此有必要進行方法改進，在保證提取的特征參數穩定性的基礎上簡化運算，采用更簡潔的方法提取艦船輻射噪聲的非線性特征，以減少計算所需時長。

數學形態學方法是一種非線性信號處理方法，其核心思想為利用各種形態特征不同的結構元素提取目標信號或進行信號降噪，能在分離背景噪聲的同時將成分復雜信號分解為若干具有物理意義的部分[8-9]，以實現目標信號提取及后續處理；且數學形態學本質上僅涉及加、減、取極值等運算，計算簡單，易于實現，因而在圖像處理、機械故障診斷以及模式識別等領域已有廣泛應用[10-13]。艦船輻射噪聲信號作為一種非線性的一維時變信號，包含產生機理相異的構成成分，與數學形態學方法有較高契合度。本文將物理意義明晰且計算簡單數學形態學方法應用于艦船輻射噪聲的非線性特征提取，與熵的概念結合并加以改進，給出新的非線性特征計算方法。

本文在討論數學形態學基本思想及數學形態譜計算方法的基礎上，以腐蝕運算替代開運算，簡化了多尺度數學形態譜計算過程，并借鑒多尺度熵(multi-scale entropy，MSE)的計算方法提出了一種廣義多尺度數學形態腐蝕譜熵(generalized multiscale pattern erosion spectrum entropy，GMPESE)計算方法，進而應用于實測艦船輻射噪聲，進行目標識別驗證。目標識別結果證明了該方法的可行性。

1 數學形態譜基礎理論

1.1 一維數學形態學運算[14]

一維形態學的基本運算包括腐蝕、膨脹、開運算和閉運算四種。假設f(n)和g(m)兩離散信號的定義域分別為：F= (0,1,…,N-1) 和G=(0,1,…,M-1)，且M≤N。定義g(m)為結構元素，則f(n)關于結構元素g(m)的腐蝕運算定義為

(fΘg)(n)=min[f(n+m)-g(m)]

(1)

膨脹運算定義為

(f⊕g)(n)=max[f(n-m)+g(n)]

(2)

在此基礎上，f(n)關于g(m)的開運算定義為

f°g=(fΘg⊕g)(n)

(3)

閉運算定義為

f·g=(f⊕gΘg)(n)

(4)

在一維形態學運算中，結構元素g起到一種輔助濾波的作用，即通過與原信號進行簡單的加減和求極值運算對一維信號進行濾波和信息提取。對于單位結構元素g，可定義λg為尺度λ下的多尺度結構元素。多尺度運算結構元素可以表示為

(5)

式中，共進行λ次腐蝕或膨脹運算。

1.2 多尺度數學形態譜

若f(n)為一維時間序列函數，g(m)為實值函數，則f(n)的數學形態譜可由式(6)計算

(6)

式中：A(f)為一維函數f在定義域內的面積； 當λ≥0時，SPS為開運算形態譜；當λ<0時，SPS為閉運算形態譜，λ的正負區間分別為信號本身的結構信息與背景信息，而兩者具有一致性，因此通常對形態譜正區間即開運算形態譜進行研究。與頻譜、功率譜等譜函數能夠描述信號在不同頻率上的成分類似，數學形態譜也可以從多個結構元素尺度上對信號進行描述。對于一維離散信號，給定信號一維離散化定義區間λ∈[0,+∞)，式(6)可簡化如式(7)所示

SPS(f,λ,g)=A[f°λg-f°(λ+1)g],λ≥0

(7)

2 數學形態腐蝕譜及廣義多尺度數學形態腐蝕譜熵特征計算方法

2.1 改進的數學形態譜——數學形態腐蝕譜

基于形態譜的特征提方法取在故障診斷等領域已有較多應用，且表現良好。但傳統形態譜方法應用于艦船輻射噪聲特征提取時，對不同的目標艦船區分度并不理想。對傳統形態譜進行改進，以腐蝕運算代替式(6)中的開運算，定義數學形態腐蝕譜(pattern erosion spectrum，PES)為

(8)

與前述過程相似，式(8)的一維離散簡化為

SPES(f,λ,g)=A[fΘλg-fΘ (λ+1)g],λ≥0

(9)

2.2 廣義多尺度數學形態腐蝕譜熵計算

自Shannon提出信息熵的概念描述信息的不確定性以來，以信息熵為基礎逐步發展出了K-S熵、近似熵、模糊熵、樣本熵和排列熵等反映信號復雜度和不確定性的熵值計算方法。

根據事件λ出現的概率q(λ)，可以計算信息熵

(10)

由于數學形態譜可在若干個結構元素尺度上描述信號，能較好的契合“某一事件出現的概率”這一概念，因此可以引入信息熵，應用數學形態學方法計算信號熵值EPESE

(11)

式中，EPESE為數學形態腐蝕譜熵(pattern erosion spectrum entropy，PESE)，q(λ)=SPS(f,λ,g)/∑SPS(f,λ,g)。

改進后的數學形態腐蝕譜熵與傳統數學形態譜熵在尺度上特性相似，即可以從多個結構元素尺度上分析評估信號，但僅在單一時間尺度上描述信號的非線性特性，未能從更全面的時間尺度分析信號的固有特性。這一點與樣本熵的局限性類似。為此，借鑒Costa等[15-16]提出的MSE計算方法，在多個時間尺度上衡量信號的非線性與復雜性特性，并將所考量的多時間尺度推廣到多時間尺度及多結構元素尺度，給出一種廣義多尺度數學形態腐蝕譜熵的計算步驟：

步驟1選取待處理艦船輻射噪聲信號x(i)，其粗粒化時間序列y(τ)可表示為

(12)

此時，粗粒化過程相當于按矩陣Y進行相空間重構

(13)

式中，M=N-(m-1)τ，m為嵌入維數，τ為時間尺度因子。對應時間尺度因子τ=1,2,…,τmax計算得到的時間序列長度為floor(N/τ)，floor為向下取整。因此當τ=1時，時間序列為原時間序列。

步驟2根據式(11)計算各粗粒化時間序列的形態學腐蝕譜熵值EPESE，并將其組合為粗粒化時間尺度因子的函數，得到廣義多尺度數學形態腐蝕譜熵

EGMPESE(x,λ,τ)=EPESE(y(τ),λ,τ)

(14)

3 艦船輻射噪聲試驗及實測數據概況

為驗證GMPESE特征提取方法在不同信道和環境噪聲條件下進行目標艦船輻射噪聲特征提取可行性，本文選取中國科學院聲學研究所東海研究站于2018年6月千島湖試驗及2020年1月東海試驗中采集的艦船輻射數據進行處理。

兩次試驗中設備布放與環境情況均有所不同。

(1) 千島湖試驗中試驗地點水深為64 m。試驗以潛標方式進行水聽器布放，水聽器工作帶寬為1 Hz～12 kHz，采樣頻率為32 kHz。兩個自容式水聽器及兩個TD分別對應固定在潛標上，入水后位于9.3 m與31.4 m深處，如圖1(a)所示。由TD數據估計得到千島湖試驗水域聲速剖面，如圖1(b)所示，為典型千島湖夏季聲速剖面。

圖1 千島湖試驗設備布放及試驗水域聲速剖面Fig.1 Equipment deployment and sound velocity profile of Qiandao Lake experiment

(2) 東海試驗中試驗地點水深為38 m。試驗采用鋼質硬連接方式布置了工作帶寬為20 Hz～20 kHz，采樣頻率為48 kHz的八元水聽器水平陣列，水聽器接收陣列固定在浮標上，如圖2(a)所示。水平陣兩端安裝有溫深儀，入水后位于8 m深處。此外，由于當日風速較大，水流較急，試驗選擇以關閉主機并拋錨的試驗船牽引浮標的方式，保證水聽器陣列不產生大范圍漂移。試驗海域聲速剖面由CTD測得，如圖2(b)所示，為典型東海淺海的冬季聲速剖面。

圖2 東海試驗設備布放及試驗水域聲速剖面Fig.2 Equipment deployment and sound velocity profile of East China Sea experiment

兩次試驗于不同季節進行。由圖2可知，其聲速剖面分別為正、負梯度，涵蓋了典型的冬、夏兩個季節的淺海(湖)聲信道情況。

選取千島湖試驗中獲取的Q1～Q4船的艦船輻射噪聲數據及東海試驗中獲取的D1～D5船的艦船輻射噪聲數據，用于后續分析及艦船目標識別驗證。Q1～Q4及D1～D5艦船輻射噪聲時域波形，如圖3所示。

圖3 實測艦船輻射噪聲時域波形Fig.3 Time-domain waveform of measured ship radiated noise

下文中，首先以東海試驗數據為例進行GMPESE相關參數及影響因素的分析，然后綜合考慮千島湖試驗數據和東海試驗數據，分析對比不同目標艦船之間GMPESE特征及MSE特征的區分度，并進行目標識別驗證。

4 參數選取影響分析

在數學形態學方法應用中，結構元素的種類與其參數選取是一個重要問題。目前常見的結構元素包括扁平結構元素、三角結構元素、半圓結構元素、正弦結構元素等。為保證信號形狀在運算過程中相對穩定，通常選取與原信號形狀接近的結構元素進行數學形態學運算。

對同一信號分別應用扁平結構元素、三角結構元素及半圓結構元素的數學形態學運算，結果如圖4所示。由圖4(c)可知，因半圓結構元素與信號形狀差異較大，進行多次數學形態學運算后難以保持信號原有形狀；由圖4(a)和圖4(b)可知，因三角結構元素及扁平結構元素與信號形狀差異較小，進行多次數學形態學運算后較好地保持了信號原有形狀，兩者相比半圓結構元素更適合于一維信號數學形態學運算，且應用扁平結構元素進行數學形態學運算時信號幅值不會產生過大變化。

選取東海試驗中不同艦船輻射噪聲信號D1～D5，分別應用扁平結構元素、三角結構元素和半圓結構元素計算20維GMPESE，結果如圖5所示。

圖5(a)在20維時間尺度上對艦船輻射噪聲信號D1～D5有很好的區分度，圖5(b)和圖5(c)則在20維時間尺度上對艦船輻射噪聲信號D1～D5區分度不佳。這一點同樣證明了前文的論述，即相對半圓結構元素和三角結構元素，扁平結構元素因其對信號形狀和幅值影響較小，信號失真度較低，因此應用扁平結構元素計算GMPESE特征區分度較好。

----為原信號； 為λ=10的膨脹與腐蝕運算； ……為λ=20的膨脹與腐蝕運算。圖4 三種不同結構元素示意Fig.4 Three different kinds of structural elements

圖5 應用不同結構元素的GMPESEFig.5 GMPESE of different kinds of structural elements

由式(14)可知，結構元素尺度因子λ和時間尺度因子τ取值均可影響GMPESE特征計算結果。而嚴格意義上來說，時間尺度因子并非結構元素的相關參數，但GMPESE特征的計算方法借鑒自MSE特征，在計算GMPESE特征的粗粒化過程中，時間尺度因子的值決定了粗粒化計算后時間序列的時長，進而影響該時間尺度下GMPESE特征計算結果，這與結構元素尺度因子對PES的影響機制相近。另一方面，時間尺度因子與結構元素尺度因子表現形式類似，因此將兩者放在同一部分進行討論。圖6中：圖6(a)為在計算GMPESE特征過程中，結構元素尺度因子λ取值對多尺度形態譜熵值的影響；圖6(b)為在由PES計算GMPESE特征過程中，時間尺度因子τ的取值對GMPESE最終值的影響。明顯可見當因子取值λ>20，τ>20時，PES和GMPESE變化幅度減小，甚至趨近平穩。在這種PES和GMPESE變化趨近平穩的情況下，該尺度因子取值對應的特征參數值對區分不同艦船目標貢獻度很小，同時會增加計算復雜度，因此在計算GMPESE時取λ=20，τ=20。該結論與文獻[17-18]一致。

圖6 兩種尺度因子對特征參數計算影響Fig.6 The influence of structural element scale factors and time scale factors on the calculation of characteristic parameters

選定結構元素并確定結構元素尺度因子及時間尺度因子后，下一步應進行結構元素長度或高度參數選取。結構元素的參數選取并無明確準則，已有研究多依據經驗參數進行取值。由于本文選用扁平結構元素進行數學形態學運算，本文以扁平結構元素參數的尋優為例，通過非經驗方法進行選取結構元素參數的嘗試，并給出一種基于貓群算法(cat swarm optimization，CSO)的尋優方法。CSO是一種由粒子群算法(particle swarm optimization，PSO)發展而來的尋優算法，本文借助CSO能跳出局部最優達到全局最優的特點對扁平結構元素的長度參數進行尋優。

此處參數尋優目的是尋找令目標識別準確率更高的結構元素長度參數，因此選取目標識別準確率作為適應度函數。參數優化的基本流程如圖7所示。

圖7 結構元素參數優化流程Fig.7 Structural element parameter optimization process

與PSO相似，CSO種群規模與迭代次數越大，越容易跳出局部最優值，尋優效果越好，但同時也需要更多時間。為平衡尋優效果和耗費時長，設置種群規模為20，以迭代50次為結束條件，在1～50內尋找合適的結構元素長度。經過迭代運算后，得出最優長度為2，即g=[0,0]為最優扁平結構元素。

當結構元素g的長度取1～15時，提取的GMPESE特征用于艦船目標識別的準確率，如圖8所示。由圖8可知，當結構元素長度為2時目標識別準確率達到最大值，與尋優結果一致。

圖8 目標識別準確率隨結構元素長度的變化Fig.8 The trend of target recognition accuracy with the length of structural elements

本章以東海試驗數據為例討論了結構元素及其他參數選取對GMPESE特征的影響，得出如下結論：

(1) 扁平結構元素對原信號幅值及形狀影響較小，更適用于艦船輻射噪聲數學形態學運算。

(2) 結構元素尺度因子及時間尺度因子分別取λ=20，τ=20，結構元素取g=[0,0]更適合艦船輻射噪聲GMPESE特征計算。

5 基于廣義多尺度數學形態腐蝕譜熵的艦船目標識別

5.1 特征及區分度分析

GMPESE特征借鑒了MSE特征計算方法，與MSE的基本思想有相通之處，因此首先對不同目標艦船輻射噪聲的GMPESE及MSE特征區分度進行比對分析。分別取410段及900段時長為1 s的千島湖試驗及東海試驗中不同目標艦船輻射噪聲，結構元素尺度因子及時間尺度因子分別取λ=20，τ=20，結構元素取g=[0,0]，計算其MSE及GMPESE，用SVM-RFE方法選取其前三階特征參數，并進行對比，不同艦船輻射噪聲在特征空間中的分布，如圖9和圖10所示。

圖9 不同艦船輻射噪聲在前3維GMPESE與MSE特征空間中分布(東海)Fig.9 Different ship radiated noise’s distribution in the first three-dimensional GMPESE and MSE feature space (East China Sea)

圖9(a)、圖9(b)分別為東海試驗中不同艦船輻射噪聲的前3維MSE特征及GMPESE特征分布。由圖9可知，在僅取前3維情況下，海試試驗中5種目標艦船的GMPESE與MSE特征均有不錯的區分度。圖10(a)、圖10(b)分別為千島湖試驗中不同艦船輻射噪聲的前3維MSE特征及GMPESE特征分布。由圖10可知：湖試試驗中4種目標艦船的MSE特征在特征空間中有所重疊，但大致有一定區分度；4種目標艦船的GMPESE特征在特征空間有所重疊，但區分度好于MSE特征。考慮到東海試驗和千島湖試驗中不同目標艦船之間差異程度并不相同，故而同樣的特征提取方法應用于不同艦船輻射噪聲數據得出的特征區分度不同。此外，經驗證在去除前3維特征，僅使用排序靠后的若干維特征的情況下，仍可在一定程度上區分目標艦船。綜上所述，可以作出推斷：僅靠前3維特征可能并不能很好地識別不同艦船目標，提取多維特征有其必要性。

圖10 不同艦船輻射噪聲在前3維GMPESE與MSE特征空間中分布(千島湖)Fig.10 Different ship radiated noise’s distribution in the first three-dimensional GMPESE and MSE feature space (Qiandao Lake)

5.2 艦船目標識別

在實際獲取艦船輻射噪聲信號時，很多時候難以保證水聽器記錄的信號都是有效的艦船輻射噪聲信號，這也就意味著有效信號常是片段化的，從短時長的信號中提取的特征參數是否有效也是判斷特征參數性能是否優秀的驗證指標之一。因此，驗證多尺度形態學方法對信號時長的依賴程度有其必要性。

分別取410段及900段時長0.1 s～1.0 s的千島湖試驗及東海試驗中不同目標艦船輻射噪聲，結構元素尺度因子及時間尺度因子分別取λ=20，τ=20，結構元素取g=[0,0]，計算MSE及GMPESE特征參數，記錄其特征參數計算的運算耗時，并應用支持向量機進行識別，對比其識別準確率及運算耗時隨信號時長的變化，如圖11所示。

圖11 MSE、GMPESE方法的目標識別準確率及運算耗時隨信號時長的變化Fig.11 The target recognition accuracy and calculation time of MSE and GMPESE’s tendency with the signal duration

由圖11(a)可知，MSE特征提取算法運算耗時與信號時長呈指數關系，即隨著信號時長的增加，計算信號MSE特征的運算耗時及其增幅均明顯增加，且在信號時長達到0.5 s時MSE特征才能穩定反映艦船輻射噪聲的非線性特性，使艦船目標識別率達到92%的可接受水平；由圖11(b)可知，GMPESE特征提取算法運算耗時可以近似看作信號時長的線性函數，在信號時長達到0.3 s時GMPESE特征即可穩定反映艦船輻射噪聲的非線性特性，使艦船目標識別率達到94.67%的可接受水平；由圖11(c)可知，MSE特征提取算法運算耗時與信號時長呈指數關系，計算信號MSE特征的運算耗時及其增幅均明顯增加，且在信號時長達到0.3 s時MSE特征才能穩定反映艦船輻射噪聲的非線性特性，使艦船目標識別率達到82.52%的可接受水平；由圖11(d)可知，GMPESE特征提取算法運算耗時可以近似看作信號時長的線性函數，在信號時長達到0.2 s時GMPESE特征即可穩定反映艦船輻射噪聲的非線性特性，使艦船目標識別率達到86.41%的可接受水平。綜合對比圖11(a)～圖11(d)可知，GMPESE方法運算耗時明顯遠低于MSE，兩者相差兩個數量級；得到可用于艦船目標識別的穩定GMPESE特征所需的艦船輻射噪聲信號時長，較MSE特征也更短。

圖11(c)、圖11(d)與圖10結果相對照，同樣證明了僅靠前3維特征并不能很好地區分不同艦船目標這一推斷的合理性，因此實際應用中取τ=20，計算多維GMPESE是有必要的。

6 結 論

本文在數學形態學方法的基礎上，提出一種廣義多尺度數學形態腐蝕譜熵特征提取方法，并在此基礎上對實測艦船輻射噪聲進行了特征提取及目標識別，結果表明：

(1) 數學形態學方法能較好反映艦船輻射噪聲的非線性特征，描述艦船輻射噪聲復雜度的廣義多尺度數學形態腐蝕譜熵特征可以作為不同目標艦船輻射噪聲的有效識別手段。

(2) GMPESE特征僅涉及加減和取極值運算，因此計算效率高，運算耗時更少，適合實際工程應用。

(3) GMPESE特征計算所需信號時長更短(千島湖環境0.2 s，東海環境0.3 s)，穩定性更好，適合實際水聲環境下艦船目標識別。