李立元 栗蘋 李國林 張廣為
(北京理工大學 機電動態控制重點實驗室, 北京 100081)
近年來,在各代干擾機、微波武器、大功率信號等一系列的干擾下,無線電引信極有可能因此出現早炸或失效的情況。 太赫茲波的工作頻率避開了高功率微波武器、雷達、干擾機等常用的工作頻段,再加上太赫茲波在大氣中的衰減較大[1-2],主要適用于近程探測,因此太赫茲波應用于引信之上具有天然的優勢[3]。
對地的無線電引信獲得的回波能量取決于地面散射系數,而地面散射系數又取決于地面類型、落角等一系列因素[4-5]。 面對不同類型的地面,回波能量具有極大的差異。 而在很長一部分中國邊境線的高原上,由于空氣的稀薄,同樣對回波能量的大小造成了影響。 綜上所述,獲取太赫茲信號在高原不同地貌下的回波特性,并實現對回波數據快速分類從而確定炸高,對于保證引信的毀傷效能顯得尤為重要。
因此,實地采集了高原各種地貌下的回波數據,本文僅僅對灌木地形下的回波雙譜用4 種積分雙譜進行降維與分類,從而在占用很少資源的情況下實現對引信炸高的識別。
高階譜是通過累積量函數(cum)來定義的,因此也被稱作累積量譜。 假如是一個零均值的k階平穩過程,則將該過程的k階累積量定義為

式中:ω為傅里葉變換到頻域得到的頻率分量;功率譜、雙譜(三階譜)、三譜(四階譜)分別為k=2,3,4 時的特殊情況。
本文所用雙譜(三階累積量譜)為

雙譜是二維函數,可以直接通過卷積神經網絡(CNN)等方式直接進行分類,但由于數據計算量過于龐大,嚴重限制了雙譜在引信上的應用。因此,可以按照一定的規則對雙譜值進行積分,將二維函數轉化為一維函數[6],從而大大減少分類時的難度及計算量,通過雙譜的對稱性,計算量則會進一步減少。 根據選定的積分路徑,將積分雙譜分為以下4 種[7]。
1.2.1 徑向積分雙譜
如圖1 所示,在雙譜平面上,積分路徑取過原點的直線,得到的積分值即為徑向積分雙譜(radially integrated bispectrum, RIB)。

圖1 徑向積分雙譜Fig.1 Radially integrated bispectrum
RIB 具有一個明顯的缺陷:如果斜率的步長選擇過大,那即使積分路徑上的所有ω2均取整,仍然會有點被遺漏,而這些點可能包含重要的信息。 如果斜率的步長選擇過小,路徑上取整后的ω2必然還會有所重復,從而造成一些雙譜值被重復積分,從而導致雙譜積分值特征不明顯。
1.2.2 軸向積分雙譜
如圖2 所示,在雙譜平面上,積分路徑取平行于ω1軸或平行于ω2軸的直線,得到的積分值即為軸向積分雙譜(axially integrated bispectrum,AIB)。AIB 使用了每一個雙譜值并確保雙譜值不會被重復使用,但其丟失了絕大部分相位信息。

圖2 軸向積分雙譜Fig.2 Axially integrated bispectrum
1.2.3 圓周積分雙譜
如圖3 所示,在雙譜平面上,積分路徑取以原點為圓心的圓周,得到的積分值即為圓周積分雙譜(circularly integrated bispectrum, CIB)。

圖3 圓周積分雙譜Fig.3 Circularly integrated bispectrum
同樣的,由于步長的原因,CIB 在積分時必定會遺漏一部分雙譜值。
1.2.4 圍線積分雙譜
如圖4 所示,在雙譜平面上,積分路徑取以原點為中心的正方形,得到的積分值即為圍線積分雙譜(surrounding-line integral bispectrum, SIB)。SIB 同樣使用了每一個雙譜值并確保雙譜值不會被重復使用,同時保留了信號的尺度信息及部分相位信息[8-9]。

圖4 圍線積分雙譜Fig.4 Surrounding-line integral bispectrum
2.1.1 目標亮點模型
如圖5(a)所示,以地面為目標時,高原上地面粗糙,檢測出的回波信號必然是多個目標回波信號的疊加,因此可能導致無法識別出正確的炸高。 圖5(b)為實采數據頻譜,通過圖像可以明顯看出,頻率值具有多個峰值,說明檢測到了不同距離目標的回波信號。

圖5 多目標回波示意圖Fig.5 Schematic diagram of multi-target echo
大部分無線電引信以地面分布的對象作為目標,為了準確地探測地面目標,必須建立一個有效的回波信號模型。
目標回波是目標對于信號源發射的信號做出的一種響應,目標的固有散射特性和發射電磁波都會對回波信號造成影響[10]。 回波特性與入射的初始時間無關,因此可以將目標看做一個線性時不變系統,回波就是對入射波的響應。 可以將單個目標的傳遞函數Hn(ω)定義為
Hn(ω)=An(ω)exp(iωτ+ i?) (4)
式中:An(ω)為幅度反射系數,是關于頻率的函數;τ為時延,是關于目標反射距離的函數,τ=2R/c,R為參考點到信號源的波程,c為光速;?為回波產生的相位差。
因此,函數可以視為一個多亮點疊加的模型,總的傳遞函數表示為

2.1.2 信號仿真
如圖6 所示,通過在Simulink 模型中添加2 個不同的距離,從而得到雙亮點接收信號的數字特征。 模型中默認幅度反射系數為相同值。

圖6 雙亮點回波Simulink 模型Fig.6 Simulink model of double bright spot echo
如圖7 所示,能很明顯看出目標亮點模型的仿真結果與圖5(b)中的實采信號包絡具有極高的相似度,一定程度上模擬出了地面目標的信號特征,之后還可以通過添加地雜波進行進一步研究。

圖7 雙目標亮點模型仿真結果Fig.7 Simulation results of double-target highlight model
如圖8 所示,將前端通過云臺綁定在標志桿上,通過改變標志桿的高度h及發射信號與地面的夾角,從而獲取不同的回波信號。 在本次實驗中,高度h分別取2 m、3 m、4 m、5 m,發射信號與地面的夾角θ分別取20°、45°、65°、85°。

圖8 數據采集平臺Fig.8 Data collection platform
2.3.1 經驗模態分解
經驗模態分解(empirical mode decomposition,EMD)是Hilbert-Huang 變換(HHT)最重要的步驟,HHT 利用EMD,得到信號各階的本征模態函數(intrinsic mode function,IMF)。EMD 過程:首先,對原始信號進行處理,求出各個局部極大值點。 然后,用插值法連接這些點組成上包絡線,同樣找到極小值點連接這些點組成下包絡線,用原始信號減去包絡線均值就得到了第一個residue,如果其不滿足IMF 定義(極值點和過零點數目至多相差一個,且包絡線均值為0),就以該residue 作為原始信號繼續重復上面的操作,直到得到第一階IMF。 最后,用原始信號減去該IMF,得到新的信號,作為初始信號再重復上面的篩選過程,直到所得到的剩余部分為單一信號或其值小于設定值時,分解完成。
2.3.2 數據分解
將回波數據通過EMD 得到數據的各階IMF及殘余分量,對各階本征模態函數進一步提取各種雙譜特征用于驗證分類的準確性。
圖9 和圖10 為初始信號及各階IMF 示意圖,其中采樣頻率為1 MHz,即采樣時間間隔為1 μs。 圖10 中縱坐標無實際意義,僅為將各階本征模態函數直觀地表示在同一張圖中的相對參考位置。

圖9 初始信號Fig.9 Initial signal

圖10 各階IMF 示意圖Fig.10 Schematic diagram of IMF in each stage
為了更直觀地看出IMF 信號瞬時頻率的變化,每兩階提取一個IMF 信號的瞬時頻率。 從圖11可以明顯看出,各階IMF 分量的瞬時頻率隨著階數的增加而降低,因此提取不同距離回波信號的同階IMF 分量雙譜特征,再次進行分類。 半數以上的當多次分類結果為同一結果時,結束對下一階IMF 分量進行分類,否則持續此步驟,此時得到的信號為最終分類結果。 回波信號分類流程如圖12所示。

圖12 回波信號分類流程Fig.12 Flow chart of echo signal classification
2.4.1 灌木數據對比
在不同的高度上,回波信號具有不同的雙譜特征,如圖13 所示,分別為在θ=45°時距地面2,3,4,5 m 處的雙譜分布。
如圖13 所示,以f1、f2這2 個頻率作為坐標軸來表示信號的偏斜程度在頻域上的分解,分別得到在θ=45°時距地面,2,3,4,5 m 處的雙譜分布。 可以明顯地看出,2,3,4,5 m 處的雙譜分布在幅值上具有明顯的區別,但對于工作平臺很小的引信來說,實時性是一個非常重要的指標,將二維數據轉換為一維數據能夠將數據量直接從x2(x為fft 次數)降為x,因此,科學有效地提取出一維雙譜切片特征對于降低運算速度具有極大的幫助。

圖13 θ =45°時不同高度雙譜分布Fig.13 Bispectrum distribution at different heights at θ =45°
2.4.2 四種積分雙譜特征
以θ=45°時4 個不同高度的4 種積分雙譜為例,得到的結果如圖14 所示。

圖14 θ =45°時不同積分雙譜分布Fig.14 Bispectrum distribution by different integral bispectrum at θ =45°
從圖14 中可以看出,將二維的雙譜特征按照4 種方式進行積分并歸一化之后,仍然具有顯著的差異,保證了一維數據的分類處理仍具有很高的可行性。
目前,主要采用的分類方法有概率神經網絡(PNN)[11]、支持向量機(SVM)[12-13]、卷積神經網絡(CNN)[14-15],由于最鄰近節點算法(KNN)計算簡便、計算速度快[16],符合引信實時性的要求,因此本文采用KNN 算法。
最鄰近節點算法也稱為k最近鄰算法,該方法簡單且直觀,如圖15 所示。 每一個樣本x通過計算到其他樣本的距離,用與其最鄰近的k個樣本來進行分類,將樣本x分為k個樣本中數量最多的類別。 由于同樣情況下回波信號的樣本集重疊很多,用在本文非常合適。 該方法無需訓練,適合應用于引信上。

圖15 最鄰近節點算法Fig.15 KNN algorithm

將4 種積分雙譜得到的4 個高度、4 個角度的數據各取30 組,利用KNN 對其進行分類,得到結果如圖16 ~圖19 所示,黑色方框中的數據分別表示在對應數據下的分類正確的數量及所占分類數據總量的比例。

圖16 RIB 分類結果Fig.16 RIB classification results

圖17 AIB 分類結果Fig.17 AIB classification results

圖18 CIB 分類結果Fig.18 CIB classification results
3.2.1 RIB 分類結果分析
對20°、45°、65°、85°這4 個角度120 組分類數據RIB 積分數據進行分類,最終可以得到4 次分類的準確率分別為61.7%、61.7%、69.2%、66.7%,利用RIB 對480 組數據分類的總準確率為64.8%。
3.2.2 AIB 分類結果分析
對20°、45°、65°、85°這4 個角度120 組分類數據AIB 積分數據進行分類,最終可以得到4 次分類的準確率分別為89.2%、91.7%、80.0%、82.5%,利用AIB 對480 組數據分類的總準確率為85.8%。
3.2.3 CIB 分類結果分析
對20°、45°、65°、85°這4 個角度120 組分類數據CIB 積分數據進行分類,最終可以得到4 次分類的準確率分別為84.2%、80.8%、80.0%、80.0%,利用CIB 對480 組數據分類的總準確率為81.3%。
3.2.4 SIB 分類結果分析
對20°、45°、65°、85°這4 個角度120 組分類數據SIB 積分數據進行分類,最終可以得到4 次分類的準確率分別為81. 7%、64. 2%、60. 0%、53.3%,利用SIB 對480 組數據分類的總準確率為64.8%。
顯然,利用KNN 對AIB 和CIB 得到的一維數據分類結果更為準確,準確率可以達到80% 以上,而對RIB 和SIB 得到的一維數據分類結果并不理想,分類的準確率只有65%左右。
通過表1、表2 可以看出,利用經EMD 對不同頻率的數據進行特征提取后,分類結果有了明顯的改進,成功率最多可提升約8.3%,分類總成功率提升3.7%。

表1 分類結果匯總Table 1 Summary of classification results

表2 優化雙譜分類結果Table 2 Optimization of bispectral classification results
本文主要圍繞雙譜分析在信號處理中的應用,對實測高原地區灌木的回波進行了分析,通過4 種不同的積分雙譜方式對雙譜特征進行了降維,進而對一維特征進行了分類,分類結果表明:
1) 利用KNN 算法,得到了對2 m、3 m、4 m、5 m這4 個高度的分類成功率,最高可以達到74.1%。
2) 雖然在理論上SIB 算法利用了所有的點數和相位信息,但在實際數據分類中效果并不好。而AIB 與CIB 兩種算法準確率較高,分別可以達到85.4% ~85.8%和81.5% ~81.3%,其中對于45°的數據,AIB 分類效果更為卓越,可以達到91.7%。
3) 利用EMD 對數據進行分解后,優化分類結果有了明顯的提高,其中AIB 成功率最高可達92.5%。
4) 說明太赫茲波段的引信完全可以通過AIB 和CIB 兩種算法實現快速準確的分類,利用HHT 優化的雙譜數據對于提升分類識別結果的準確度具有一定的作用。