馬 凱， 陳 喆， 王易川， 程玉勝

(海軍潛艇學(xué)院，山東 青島 266000)

艦船輻射噪聲[1]調(diào)制譜特征是水聲目標(biāo)識(shí)別的重要特征，具有很好的穩(wěn)健性。對(duì)調(diào)制譜進(jìn)行諧波分析可以獲得螺旋槳的軸頻、葉頻和槳葉數(shù)等特征[2-4]，其中，軸頻提取對(duì)目標(biāo)識(shí)別具有重要意義。

許多學(xué)者對(duì)軸頻地提取進(jìn)行了深入研究[5-6]。文獻(xiàn)[7]模擬聲納兵處理DEMON(detection of modulation spectrum on noise)譜的方法，利用雙門(mén)限提高了水聲目標(biāo)輻射噪聲調(diào)制譜檢測(cè)能力，并通過(guò)提取到的線譜間隔估計(jì)軸頻，算法在單目標(biāo)情況下效果較好。文獻(xiàn)[8]通過(guò)最大公約數(shù)法提取目標(biāo)軸頻，但最大公約數(shù)法在某些多目標(biāo)情況下不適用。Di Martino等[9]提出一種基于代價(jià)函數(shù)的線譜提取算法，用于低頻線譜地提取。該算法通過(guò)設(shè)定一代價(jià)函數(shù)，使符合線譜特征的點(diǎn)可以提取出來(lái)，對(duì)于弱線譜地提取效果較好，但該算法只能提取單根線譜。李山等[10]通過(guò)設(shè)定頻率滑動(dòng)窗實(shí)現(xiàn)了多根線譜地提取，但算法并不是全局最優(yōu)的，當(dāng)信噪比較低時(shí)，算法效果并不理想。針對(duì)上述情況，本文提出一種基于蟻群算法的調(diào)制譜線譜提取算法。蟻群算法是由Dorigo等[11]提出的一種比較成熟的全局優(yōu)化算法，被廣泛應(yīng)用于機(jī)器人路徑規(guī)劃、任務(wù)規(guī)劃問(wèn)題等[12-15]。該算法是一種模擬螞蟻覓食行為的智能群算法，具有正反饋、強(qiáng)魯棒性、并行性、較好的適應(yīng)性等優(yōu)點(diǎn)[16]。

1 基于蟻群算法的線譜提取算法

蟻群算法的原理源于模擬螞蟻種群的覓食行為，螞蟻在覓食時(shí)會(huì)在搜尋路徑上留下一種名為信息素的分泌物，經(jīng)過(guò)此路徑的螞蟻越多，信息素濃度就越高。而信息素會(huì)對(duì)后面的螞蟻產(chǎn)生一種正向的吸引作用，即信息素濃度越大，則后面的螞蟻選擇此條路線的概率越大。同時(shí)，部分信息素會(huì)揮發(fā)，有利于蟻群進(jìn)行其他路線的搜索，提高全局搜索能力，避免陷入局部最優(yōu)。

本文將調(diào)制譜中的軸頻及其倍頻線譜視為螞蟻尋優(yōu)的“路徑”，通過(guò)改進(jìn)尋優(yōu)策略，使其適用于調(diào)制譜中線譜地提取。

1.1 蟻群算法描述

蟻群算法的狀態(tài)轉(zhuǎn)移概率公式為

(1)

(2)

式中，dij為節(jié)點(diǎn)i到節(jié)點(diǎn)j的距離。

螞蟻每完成一次循環(huán)后路徑信息素更新為

τij(t+1)=(1-ρ)τij(t)+Δτij(t)

(3)

(4)

式中：τij(t+1)為更新后節(jié)點(diǎn)i到節(jié)點(diǎn)j的信息素強(qiáng)度；ρ為信息素?fù)]發(fā)系數(shù)；m為蟻群數(shù)量；k為第k只螞蟻；Δτij(t)為節(jié)點(diǎn)i到節(jié)點(diǎn)j的信息素增加量，其中Δτij(t)的計(jì)算方式采用蟻周型，如式(5)所示。

(5)

式中：Q為信息素增加的強(qiáng)度系數(shù)；Lk為循環(huán)結(jié)束后螞蟻經(jīng)過(guò)的路徑長(zhǎng)度。

1.2 基于蟻群算法的線譜提取算法

采用柵格法[17]將DEMON譜圖劃分成一個(gè)個(gè)小方格，將螞蟻隨機(jī)布放于小方格之中，令其對(duì)整個(gè)DEMON譜圖進(jìn)行搜索。其中，在應(yīng)用算法時(shí)存在以下幾個(gè)問(wèn)題。

(1) 在經(jīng)典的蟻群算法中，螞蟻可以隨機(jī)選擇周?chē)?個(gè)方向進(jìn)行搜尋，如圖1所示。但在DEMON譜圖中，每條線譜在每一時(shí)刻的頻率是單一的，即螞蟻只能隨著時(shí)間增加或減少的方向進(jìn)行搜索，而不能橫向搜索。

圖1 螞蟻可選搜尋方向圖Fig.1 Available search directions for ants

(2) 在經(jīng)典蟻群算法中，螞蟻的出發(fā)點(diǎn)和目的地一般是固定的，而在DEMON譜圖中，螞蟻的出發(fā)點(diǎn)則位于起始時(shí)刻任意頻率點(diǎn)處，目的地位于終止時(shí)刻任意頻率點(diǎn)處，位置不固定。

(3) 在經(jīng)典蟻群算法中，啟發(fā)函數(shù)以及信息素的更新是根據(jù)螞蟻距目的地的距離來(lái)確定的，而在DEMON 譜圖中，由于“目的地”未知，所以無(wú)法依據(jù)“距離”這一標(biāo)準(zhǔn)來(lái)確定。

基于以上幾點(diǎn)問(wèn)題，提出一種基于蟻群算法的線譜提取算法，并作以下改進(jìn)。

(1) 對(duì)螞蟻的搜尋范圍進(jìn)行限定，即每只螞蟻只能在時(shí)間軸方向進(jìn)行搜索，而不能橫向搜索，如圖2所示?？v坐標(biāo)為時(shí)間，橫坐標(biāo)為頻率。

圖2 螞蟻在DEMON譜圖中的可選搜尋方向圖Fig.2 Available search directions for ants in lofargram

(2) 提出一種新的代價(jià)函數(shù)作為決定啟發(fā)函數(shù)和信息素更新的標(biāo)準(zhǔn)，如式(6)所示。

(6)

式中：ξ為線譜路徑；d(i)為線譜的頻率連續(xù)性特征；G(i)為線譜的軌跡連續(xù)性特征；A(i)為線譜的強(qiáng)度特征，即線譜路徑上的各個(gè)頻率點(diǎn)上的幅值；eps為算法精度，用以防止出現(xiàn)分母為0的情況。則由式(5)可知，當(dāng)線譜的能量越強(qiáng)，即A(i)越大，頻率連續(xù)性越好，即d(i)越小，軌跡連續(xù)性越好，即G(i)越大，則線譜代價(jià)函數(shù)越小，即P_cos(ξ)越小。

d(i)，G(i)，A(i)的定義式分別為

d(i)=|f(Pi)-f(Pi+1)|+eps

(7)

(8)

(9)

式中：f(Pi)為線譜Pi在i點(diǎn)處的頻率值；g(Pi)為一個(gè)標(biāo)記變量，其定義如式(9)所示，即當(dāng)g(Pi)小于設(shè)定的線譜幅值門(mén)限count時(shí)，則認(rèn)為此處出現(xiàn)斷點(diǎn)，否則設(shè)為1。

則新的啟發(fā)函數(shù)為

(10)

由式(10)可知，P_cos(ξ)越小，即搜尋路徑越接近線譜時(shí)，啟發(fā)函數(shù)值越大，符合啟發(fā)函數(shù)的設(shè)計(jì)原則。

(3) 隨著迭代次數(shù)的增加，路徑上的信息素濃度易積累過(guò)高，導(dǎo)致螞蟻的路徑選擇能力在一定程度上降低，極可能增加陷入局部最優(yōu)解的概率。因此，對(duì)信息素增加的強(qiáng)度系數(shù)Q進(jìn)行改進(jìn)，使其隨迭代次數(shù)的增加自適應(yīng)的調(diào)整，如式(11)所示，當(dāng)?shù)螖?shù)大于設(shè)定的閾值時(shí)，Q值隨迭代次數(shù)的增加逐漸減小，從而降低算法陷入局部最優(yōu)的概率。

(11)

式中，Ne為當(dāng)前迭代次數(shù)。

算法處理流程，如圖3所示。其中C為設(shè)定的迭代次數(shù)。

算法處理步驟如下：

步驟1對(duì)獲取的信號(hào)進(jìn)行解調(diào)處理，得到DEMON譜圖，并將其柵格化處理。

步驟2初始化禁忌表，禁忌表用以記錄螞蟻訪問(wèn)過(guò)的節(jié)點(diǎn)。

步驟3將M只螞蟻隨機(jī)放置于DEMON譜圖起始時(shí)刻的各頻率點(diǎn)處，并計(jì)算每只螞蟻到下個(gè)節(jié)點(diǎn)的代價(jià)函數(shù)P_cos(ξ)。

步驟4由代價(jià)函數(shù)P_cos(ξ)計(jì)算得到啟發(fā)因子ηij(t)。

步驟6更新禁忌表，并記錄每只螞蟻?zhàn)哌^(guò)的路徑。

步驟7當(dāng)螞蟻遍歷整個(gè)DEMON譜圖的時(shí)間軸時(shí)，一次迭代結(jié)束，并重復(fù)上述過(guò)程。

步驟8記錄代價(jià)函數(shù)最小的n條路徑，并進(jìn)行二次閾值判斷，最終得到所需要的線譜。

圖3 算法處理流程圖Fig.3 Flowchart of algorithm processing

2 諧波庫(kù)匹配算法

由于螺旋槳結(jié)構(gòu)設(shè)計(jì)以及環(huán)境等因素，所提取的線譜中可能會(huì)包含虛假線譜，也可能缺失軸頻及其諧波線譜。并且，在多目標(biāo)情況下，線譜結(jié)構(gòu)更復(fù)雜，軸頻地自動(dòng)提取也更困難。常用的軸頻提取算法包括最大公約數(shù)法、余數(shù)門(mén)限法[19]、倍頻法等，此類(lèi)算法在某些多目標(biāo)情況下效果較差。針對(duì)這種情況，提出一種諧波庫(kù)匹配算法，自動(dòng)提取多目標(biāo)軸頻。

算法處理流程如下：

步驟1將提取到的n條線譜從小到大排列，得到線譜序列{fm}，m=1,2,…,M。

步驟2根據(jù)現(xiàn)有船舶螺旋槳結(jié)構(gòu)可知，螺旋槳槳葉數(shù)通常最多為7葉。因此，根據(jù)已經(jīng)提取的n條線譜建立一個(gè)大小為n×7的諧波庫(kù)矩陣H。

步驟3將每根線譜fm分別假設(shè)為軸頻的1～7次諧波，其中1次諧波為軸頻。定義一品質(zhì)因數(shù)Sk，用以記錄fm為k次諧波時(shí)，序列{fm}中符合式(12)條件的線譜個(gè)數(shù)。

(12)

式中，Δ為誤差控制量，依據(jù)實(shí)際情況確定，用以表示由非嚴(yán)格諧波情況造成的誤差。

步驟4將線譜fm的品質(zhì)因數(shù)Sk存入諧波庫(kù)H(m,k)中，統(tǒng)計(jì)每條線譜fm中品質(zhì)因數(shù)最大時(shí)的諧波次數(shù)kmax，并按照式(13)計(jì)算出線譜fm的基頻Bm。

(13)

步驟5此時(shí)，基頻序列{Bm}中可能由于非嚴(yán)格諧波關(guān)系存在弱小的誤差，將{Bm}中的每個(gè)基頻返回到DEMON譜圖中，分別計(jì)算基頻為Bm時(shí)的各線譜的能量Pm。

步驟6將頻率之差小于δHz的基頻歸為一類(lèi)，分別計(jì)算每一類(lèi)中Pm最大時(shí)的基頻Bm，并比較這幾類(lèi)中Pm值的大小，取最大的兩個(gè)P1和P2,其中P1≥P2，若P2大于噪聲背景均值，則基頻為B1和B2，否則，認(rèn)為目標(biāo)為單目標(biāo)，基頻為B1。

3 仿真及試驗(yàn)驗(yàn)證

3.1 算法處理流程

本文算法處理流程，如圖4所示。

步驟1對(duì)輸入信號(hào)S進(jìn)行絕對(duì)值解調(diào)，獲得DEMON譜。

步驟2利用文獻(xiàn)[20]中的排序截短算法對(duì)DEMON譜進(jìn)行背景均衡，剔除趨勢(shì)項(xiàng)。

步驟3通過(guò)提出的蟻群算法提取DEMON譜中的線譜。

步驟4利用諧波庫(kù)匹配算法實(shí)現(xiàn)多目標(biāo)信號(hào)軸頻地自動(dòng)提取。

圖4 本文算法處理流程圖Fig.4 Algorithm processing flow chart of this paper

3.2 仿真驗(yàn)證

仿真中假設(shè)有兩個(gè)目標(biāo)，軸頻分別為6 Hz，9 Hz，5葉槳，信噪比設(shè)為-17 dB，采樣頻率為5 000 Hz，其中在15 Hz處添加了一條干擾線譜，并且為了模擬實(shí)際信號(hào)中的非嚴(yán)格周期現(xiàn)象，分別在軸頻倍頻處加入一弱小的擾動(dòng)。仿真參數(shù)設(shè)置：α=7，β=7，Rho=0.1，Q=100，m=400，L=100，所有參數(shù)均經(jīng)過(guò)蒙特卡洛仿真試驗(yàn)確定，綜合性能最優(yōu)。

3.2.1 背景均衡

解調(diào)后的原始的DEMON譜圖，如圖5所示。背景均衡后的DEMON譜圖，如圖6所示。由圖5、圖6可知，此時(shí)信噪比較低，線譜能量較弱，部分線譜幾乎不可見(jiàn)，經(jīng)過(guò)背景均衡，剔除趨勢(shì)項(xiàng)后，算法的背景噪聲有所降低。

圖5 原始DEMON譜圖Fig.5 Original DEMON

圖6 背景均衡后的DEMON譜圖Fig.6 DEMON after background equalization

3.2.2 線譜提取

經(jīng)過(guò)最優(yōu)路徑算法和本文提出的蟻群算法所提出的線譜，如圖7和圖8所示。由圖7、圖8可知，最優(yōu)路徑算法可以提取部分較弱的線譜，但同時(shí)會(huì)遺漏部分線譜，如18 Hz附近的兩條線譜會(huì)漏掉一根。本文所提算法效果較好，即使部分肉眼幾乎不可見(jiàn)的線譜也能較好地提取出來(lái)。

圖7 最優(yōu)路徑算法所提線譜Fig.7 Line spectra proposed by optimal path algorithm

圖8 本文算法所提線譜Fig.8 Line spectra proposed by the algorithm in this paper

不同信噪比下，兩種算法線譜提取性能的比較，如圖9所示。為了比較兩種算法的線譜提取能力，現(xiàn)定義線譜提取率P，如式(14)所示。由圖9可知，最優(yōu)路徑算法無(wú)法提取相鄰較近的兩根線譜，其線譜提取的準(zhǔn)確率始終低于本文算法。這是因?yàn)楸疚乃崴惴ㄊ且环N全局尋優(yōu)的算法，當(dāng)信噪比較低時(shí)，算法受離散能量較高的噪點(diǎn)的影響較小，而最優(yōu)路徑算法受能量較高的噪點(diǎn)的影響較大，從而提取能力隨信噪比的降低下降較為嚴(yán)重。但本文算法花費(fèi)時(shí)間較最優(yōu)路徑算法長(zhǎng)，用時(shí)15 s，最優(yōu)路徑算法用時(shí)5 s，下一步將通過(guò)并行處理等手段提高算法的運(yùn)算速度。

(14)

圖9 兩種算法的線譜提取性能對(duì)比圖Fig.9 Comparison of the extraction performance of the two algorithms

3.2.3 軸頻自動(dòng)提取

(1) 最大公約數(shù)法

該算法首先計(jì)算了所提線譜中任意兩條線譜之差，并定義了3種品質(zhì)因數(shù)：Sq1，Sq2和Sq3。其中：Sq1為上述差值中相同數(shù)的個(gè)數(shù)；Sq2為上述差值中能被所提線譜頻率除盡的數(shù)的個(gè)數(shù)；Sq3為Sq1和Sq2的乘積，Sq3最大時(shí)所對(duì)應(yīng)的頻率即為軸頻。Sq3最大值為77，所對(duì)應(yīng)的頻率為3 Hz，很明顯不是所求軸頻，所以該算法在某些多目標(biāo)情況下不適用，如表1所示。

表1 最大公約數(shù)法所得結(jié)果Tab.1 Results obtained by the maximum common divisor method

(2) 余數(shù)門(mén)限法

該算法首先定義了一個(gè)3×n的矩陣，n為設(shè)定的線譜的個(gè)數(shù)，矩陣的第一列用于存放線譜，第二列存放線譜的初始權(quán)值ω，第三列存放判決依據(jù)F，F(xiàn)最大時(shí)所對(duì)應(yīng)的頻率即為所求軸頻，如表2所示。由表2可知，當(dāng)F最大時(shí)，頻率為9 Hz，即所求軸頻為9 Hz，只能求得其中的某一個(gè)目標(biāo)的軸頻，因此該算法在某些多目標(biāo)情況下是不適用的。

表2 余數(shù)門(mén)限法所得結(jié)果Tab.2 Results obtained by the remainder threshold algorithm

(3) 諧波庫(kù)匹配算法

本文所提算法經(jīng)式(12)處理后的基頻，如表3所示。由表3可知，算法得到的基頻分為兩類(lèi)，每一類(lèi)中的Pm的最大值所對(duì)應(yīng)的基頻分別為6 Hz和9 Hz，與噪聲背景均值比較最終可得目標(biāo)的軸頻分別為6 Hz和9 Hz，與仿真設(shè)定的軸頻一致，效果較好。

表3 本文算法所得結(jié)果Tab.3 The results obtained by the algorithm in this paper

3.3 海試數(shù)據(jù)驗(yàn)證

該信號(hào)為多目標(biāo)信號(hào)，經(jīng)過(guò)先驗(yàn)信息以及人工綜合判定，確定兩目標(biāo)軸頻分別為1.52 Hz和1.88 Hz，信號(hào)采樣頻率為44 100 Hz，頻率分辨率為0.01 Hz，線譜提取參數(shù)設(shè)置同仿真。

3.3.1 背景均衡

解調(diào)后的原始DEMON譜圖，如圖10所示。其中部分線譜能量較弱，幾乎淹沒(méi)于背景噪聲中。背景均衡后的DEMON譜圖，如圖11所示。經(jīng)過(guò)背景均衡，算法的背景噪聲有所降低。

圖10 原始DEMON譜圖(海試數(shù)據(jù))Fig.10 Original DEMON(sea trial data)

圖11 背景均衡后的DEMON譜圖(海試數(shù)據(jù))Fig.11 DEMON after background equalization(sea trial data)

3.3.2 線譜提取

經(jīng)過(guò)最優(yōu)路徑算法和本文提出的蟻群算法所提出的線譜，如圖12和圖13所示。由圖12、圖13可知，最優(yōu)路徑算法可以提取部分較弱的線譜，但同時(shí)會(huì)遺漏部分線譜，例如2.96 Hz，4.44 Hz處的線譜。本文所提算法可以將能量非常低的線譜提取出來(lái)，效果較好。

圖12 最優(yōu)路徑算法所提線譜(海試數(shù)據(jù))Fig.12 Line spectra proposed by optimal path algorithm(sea trial data)

圖13 本文算法所提線譜(海試數(shù)據(jù))Fig.13 Line spectra proposed by the algorithm in this paper(sea trial data)

3.3.3 軸頻自動(dòng)提取

最大公約數(shù)法得到的結(jié)果，如表4所示。Sq3最大值為20，所對(duì)應(yīng)的頻率為1.82 Hz，與實(shí)際軸頻有所差距，并且只能提取單個(gè)軸頻，效果欠佳。

余數(shù)門(mén)限法所求結(jié)果，如表5所示。當(dāng)F最大時(shí)對(duì)應(yīng)的頻率為1.52 Hz，即所求軸頻為1.52 Hz，只能求得單個(gè)目標(biāo)的軸頻，因此在某些多目標(biāo)情況下是不適用的。

表4 最大公約數(shù)法所得結(jié)果(海試數(shù)據(jù))Tab.4 Results obtained by the maximum common divisor method(sea trial data)

表5 余數(shù)門(mén)限法所得結(jié)果(海試數(shù)據(jù))Tab.5 Results obtained by the remainder threshold algorithm(sea trial data)

本文所提算法的結(jié)果，如表6所示。由表6可知，算法得到的基頻分為兩類(lèi)，每一類(lèi)中的Pm的最大值所對(duì)應(yīng)的基頻分別為1.52 Hz和1.88 Hz，取其中最大的兩個(gè)值并與噪聲背景均值比較，最終可得，目標(biāo)的軸頻分別為1.52 Hz和1.88 Hz，與實(shí)際軸頻一致，效果較好。

表6 本文算法所得結(jié)果(海試數(shù)據(jù))Tab.6 The results obtained by the algorithm in this paper (sea trial data)

綜上分析可知，本文所提出的基于蟻群算法的線譜提取算法相較于最優(yōu)路徑算法效果更好，可以提出能量很弱的線譜；在多目標(biāo)條件下，提出的基于諧波庫(kù)匹配算法相比于過(guò)去的最大公約數(shù)法以及余數(shù)門(mén)限法效果更好。最大公約數(shù)法及余數(shù)門(mén)限法在某些多目標(biāo)情況下所提出的軸頻可能不是真實(shí)軸頻，且只能提取單個(gè)軸頻，本文算法則能夠解決此問(wèn)題。

4 結(jié) 論

本文針對(duì)多目標(biāo)情況，提出了一種軸頻自動(dòng)提取算法。首先利用提出的蟻群算法提取DEMON譜圖中的線譜，相比于最優(yōu)路徑算法，可以更好地提取弱線譜。然后利用提出的諧波庫(kù)匹配算法實(shí)現(xiàn)多目標(biāo)軸頻地自動(dòng)提取，解決了以往最大公約數(shù)算法和余數(shù)門(mén)限法只能提取單個(gè)軸頻以及在某些多目標(biāo)情況下提取錯(cuò)誤軸頻的問(wèn)題，效果較好。