低信噪比下多目標調制譜軸頻自動檢測算法

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

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

艦船輻射噪聲[1]調制譜特征是水聲目標識別的重要特征，具有很好的穩健性。對調制譜進行諧波分析可以獲得螺旋槳的軸頻、葉頻和槳葉數等特征[2-4]，其中，軸頻提取對目標識別具有重要意義。

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

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

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

本文將調制譜中的軸頻及其倍頻線譜視為螞蟻尋優的“路徑”，通過改進尋優策略，使其適用于調制譜中線譜地提取。

1.1 蟻群算法描述

蟻群算法的狀態轉移概率公式為

(1)

(2)

式中，dij為節點i到節點j的距離。

螞蟻每完成一次循環后路徑信息素更新為

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

(3)

(4)

式中：τij(t+1)為更新后節點i到節點j的信息素強度；ρ為信息素揮發系數；m為蟻群數量；k為第k只螞蟻；Δτij(t)為節點i到節點j的信息素增加量，其中Δτij(t)的計算方式采用蟻周型，如式(5)所示。

(5)

式中：Q為信息素增加的強度系數；Lk為循環結束后螞蟻經過的路徑長度。

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

采用柵格法[17]將DEMON譜圖劃分成一個個小方格，將螞蟻隨機布放于小方格之中，令其對整個DEMON譜圖進行搜索。其中，在應用算法時存在以下幾個問題。

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

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

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

(3) 在經典蟻群算法中，啟發函數以及信息素的更新是根據螞蟻距目的地的距離來確定的，而在DEMON 譜圖中，由于“目的地”未知，所以無法依據“距離”這一標準來確定。

基于以上幾點問題，提出一種基于蟻群算法的線譜提取算法，并作以下改進。

(1) 對螞蟻的搜尋范圍進行限定，即每只螞蟻只能在時間軸方向進行搜索，而不能橫向搜索，如圖2所示。縱坐標為時間，橫坐標為頻率。

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

(2) 提出一種新的代價函數作為決定啟發函數和信息素更新的標準，如式(6)所示。

(6)

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

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

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

(7)

(8)

(9)

式中：f(Pi)為線譜Pi在i點處的頻率值；g(Pi)為一個標記變量，其定義如式(9)所示，即當g(Pi)小于設定的線譜幅值門限count時，則認為此處出現斷點，否則設為1。

則新的啟發函數為

(10)

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

(3) 隨著迭代次數的增加，路徑上的信息素濃度易積累過高，導致螞蟻的路徑選擇能力在一定程度上降低，極可能增加陷入局部最優解的概率。因此，對信息素增加的強度系數Q進行改進，使其隨迭代次數的增加自適應的調整，如式(11)所示，當迭代次數大于設定的閾值時，Q值隨迭代次數的增加逐漸減小，從而降低算法陷入局部最優的概率。

(11)

式中，Ne為當前迭代次數。

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

算法處理步驟如下：

步驟1對獲取的信號進行解調處理，得到DEMON譜圖，并將其柵格化處理。

步驟2初始化禁忌表，禁忌表用以記錄螞蟻訪問過的節點。

步驟3將M只螞蟻隨機放置于DEMON譜圖起始時刻的各頻率點處，并計算每只螞蟻到下個節點的代價函數P_cos(ξ)。

步驟4由代價函數P_cos(ξ)計算得到啟發因子ηij(t)。

步驟6更新禁忌表，并記錄每只螞蟻走過的路徑。

步驟7當螞蟻遍歷整個DEMON譜圖的時間軸時，一次迭代結束，并重復上述過程。

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

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

2 諧波庫匹配算法

由于螺旋槳結構設計以及環境等因素，所提取的線譜中可能會包含虛假線譜，也可能缺失軸頻及其諧波線譜。并且，在多目標情況下，線譜結構更復雜，軸頻地自動提取也更困難。常用的軸頻提取算法包括最大公約數法、余數門限法[19]、倍頻法等，此類算法在某些多目標情況下效果較差。針對這種情況，提出一種諧波庫匹配算法，自動提取多目標軸頻。

算法處理流程如下：

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

步驟2根據現有船舶螺旋槳結構可知，螺旋槳槳葉數通常最多為7葉。因此，根據已經提取的n條線譜建立一個大小為n×7的諧波庫矩陣H。

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

(12)

式中，Δ為誤差控制量，依據實際情況確定，用以表示由非嚴格諧波情況造成的誤差。

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

(13)

步驟5此時，基頻序列{Bm}中可能由于非嚴格諧波關系存在弱小的誤差，將{Bm}中的每個基頻返回到DEMON譜圖中，分別計算基頻為Bm時的各線譜的能量Pm。

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

3 仿真及試驗驗證

3.1 算法處理流程

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

步驟1對輸入信號S進行絕對值解調，獲得DEMON譜。

步驟2利用文獻[20]中的排序截短算法對DEMON譜進行背景均衡，剔除趨勢項。

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

步驟4利用諧波庫匹配算法實現多目標信號軸頻地自動提取。

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

3.2 仿真驗證

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

3.2.1 背景均衡

解調后的原始的DEMON譜圖，如圖5所示。背景均衡后的DEMON譜圖，如圖6所示。由圖5、圖6可知，此時信噪比較低，線譜能量較弱，部分線譜幾乎不可見，經過背景均衡，剔除趨勢項后，算法的背景噪聲有所降低。

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

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

3.2.2 線譜提取

經過最優路徑算法和本文提出的蟻群算法所提出的線譜，如圖7和圖8所示。由圖7、圖8可知，最優路徑算法可以提取部分較弱的線譜，但同時會遺漏部分線譜，如18 Hz附近的兩條線譜會漏掉一根。本文所提算法效果較好，即使部分肉眼幾乎不可見的線譜也能較好地提取出來。

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

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

不同信噪比下，兩種算法線譜提取性能的比較，如圖9所示。為了比較兩種算法的線譜提取能力，現定義線譜提取率P，如式(14)所示。由圖9可知，最優路徑算法無法提取相鄰較近的兩根線譜，其線譜提取的準確率始終低于本文算法。這是因為本文所提算法是一種全局尋優的算法，當信噪比較低時，算法受離散能量較高的噪點的影響較小，而最優路徑算法受能量較高的噪點的影響較大，從而提取能力隨信噪比的降低下降較為嚴重。但本文算法花費時間較最優路徑算法長，用時15 s，最優路徑算法用時5 s，下一步將通過并行處理等手段提高算法的運算速度。

(14)

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

3.2.3 軸頻自動提取

(1) 最大公約數法

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

表1 最大公約數法所得結果Tab.1 Results obtained by the maximum common divisor method

(2) 余數門限法

該算法首先定義了一個3×n的矩陣，n為設定的線譜的個數，矩陣的第一列用于存放線譜，第二列存放線譜的初始權值ω，第三列存放判決依據F，F最大時所對應的頻率即為所求軸頻，如表2所示。由表2可知，當F最大時，頻率為9 Hz，即所求軸頻為9 Hz，只能求得其中的某一個目標的軸頻，因此該算法在某些多目標情況下是不適用的。

表2 余數門限法所得結果Tab.2 Results obtained by the remainder threshold algorithm

(3) 諧波庫匹配算法

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

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

3.3 海試數據驗證

該信號為多目標信號，經過先驗信息以及人工綜合判定，確定兩目標軸頻分別為1.52 Hz和1.88 Hz，信號采樣頻率為44 100 Hz，頻率分辨率為0.01 Hz，線譜提取參數設置同仿真。

3.3.1 背景均衡

解調后的原始DEMON譜圖，如圖10所示。其中部分線譜能量較弱，幾乎淹沒于背景噪聲中。背景均衡后的DEMON譜圖，如圖11所示。經過背景均衡，算法的背景噪聲有所降低。

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

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

3.3.2 線譜提取

經過最優路徑算法和本文提出的蟻群算法所提出的線譜，如圖12和圖13所示。由圖12、圖13可知，最優路徑算法可以提取部分較弱的線譜，但同時會遺漏部分線譜，例如2.96 Hz，4.44 Hz處的線譜。本文所提算法可以將能量非常低的線譜提取出來，效果較好。

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

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

3.3.3 軸頻自動提取

最大公約數法得到的結果，如表4所示。Sq3最大值為20，所對應的頻率為1.82 Hz，與實際軸頻有所差距，并且只能提取單個軸頻，效果欠佳。

余數門限法所求結果，如表5所示。當F最大時對應的頻率為1.52 Hz，即所求軸頻為1.52 Hz，只能求得單個目標的軸頻，因此在某些多目標情況下是不適用的。

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

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

本文所提算法的結果，如表6所示。由表6可知，算法得到的基頻分為兩類，每一類中的Pm的最大值所對應的基頻分別為1.52 Hz和1.88 Hz，取其中最大的兩個值并與噪聲背景均值比較，最終可得，目標的軸頻分別為1.52 Hz和1.88 Hz，與實際軸頻一致，效果較好。

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

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

4 結 論

本文針對多目標情況，提出了一種軸頻自動提取算法。首先利用提出的蟻群算法提取DEMON譜圖中的線譜，相比于最優路徑算法，可以更好地提取弱線譜。然后利用提出的諧波庫匹配算法實現多目標軸頻地自動提取，解決了以往最大公約數算法和余數門限法只能提取單個軸頻以及在某些多目標情況下提取錯誤軸頻的問題，效果較好。