周毅恒, 楊 軍, 夏賽強, 呂明久
(空軍預警學院, 湖北 武漢 430019)
直升機飛行過程中,其旋翼高速轉動,會使雷達回波信號產生額外的頻率調制,稱為微多普勒現象,此時回波信號中會包含直升機的旋翼長度、旋轉頻率等信息,如果能夠較為精確地提取出這些物理特征,則可實現對機型進行穩健識別,對國土安全具有重要意義。
對于旋翼類目標的微動參數估計問題,目前主要有以下幾種方法:一是基于變換域的微動參數估計方法,在時頻域通過逆Radon變換、Hough變換等方法將邊緣檢測問題轉換為峰值提取問題,此類方法估計精度較高,魯棒性強,但是運算量較大;二是基于圖像域的估計方法,通過對時頻圖進行骨架提取的操作獲得目標微動曲線,進而實現對微動目標的參數估計,該方法易受噪聲影響,隨著信噪比(signal to noise ratio,SNR)下降,估計精度將會越來越差。實測數據表明,閃爍現象在旋翼回波中普遍存在,上述方法在該條件下均難以有效估計出目標的微動參數,主要問題在于當旋翼目標微動回波出現閃爍時,回波信號幅度項增加了辛格(sinc(·))函數調制,時頻域中會出現余弦包絡、零頻帶和閃爍帶,雖然余弦包絡包含目標的微動信息,但是易受其余微多普勒特征的影響,檢測難度很大。
近年來,深度學習在圖像識別、信號處理領域得到了廣泛運用,為解決上述問題提供了一種新的思路。本文在研究閃爍現象下旋翼回波時頻圖特點的基礎上,提出一種有效的旋翼微動參數估計方法,針對零頻帶、閃爍帶以及大量無規律分布的噪點影響余弦包絡檢測的問題,將其都等效為時頻圖像中的噪聲,利用去噪卷積神經網絡(denosing convolutional neural network,DnCNN)架構分別訓練去噪、去閃爍網絡,消除噪聲、零頻帶和閃爍帶的影響。針對傳統逆Radon變換采用遍歷法進行微動參數搜索,計算量較大的問題,采用黃金分割法對其搜索過程進行改進,提高運算速度,進而完成對旋翼目標微動參數的估計。
不失一般性,假設雷達旋翼回波已經經過平動補償,建立以雷達所在位置為原點的OXYZ空間直角坐標系,如圖1所示。該直升機擁有一個旋翼,多個葉片,以轉動角頻率正向旋轉。各個旋翼葉片長度均為,雷達到直升機旋翼中心(葉轂)的距離為,旋翼中心、雷達以及地面所形成的夾角為,方位角為。

圖1 雷達和旋翼葉片的幾何關系圖Fig.1 Geometric relationship between radar and rotor blades
設雷達發射信號為線性調頻信號:

(1)

散射點模型下直升機旋翼葉片可等效為多個散射點,如圖2所示。

圖2 旋翼葉片散射點模型圖Fig.2 Scattering point model of rotor blades
旋翼葉片上散射點的回波信號進行脈沖壓縮后,可表示為

(2)
式中:為葉片數目;為旋翼上的第個葉片;為單個葉片上的散射點數目;為葉片上的第個散射點; 為對應散射點的散射系數;為觀測目標的時間;sinc(·)表示辛格函數;為線性調頻信號的帶寬; ()表示旋翼上對應散射點到雷達的距離;c為光速; 是旋翼葉片對應散射點的初相角;為發射信號的波長;′()表示脈壓后的噪聲。對應圖1,令方位角=0,有:
()=+ cos( +)cos
(3)

當旋翼葉片上各個散射點之間的距離 <2時,旋翼回波將出現閃爍現象,如圖3所示,單旋翼雙葉片散射點模型的閃爍現象仿真結果圖。圖3(a)表示=0、=0時葉片與雷達相對位置圖,圖3(b)為=03 m、=01 m、=120(每個葉片上的散射點數目)、=6 m、=12π rad/s、=0512 s、兩個葉片初始角度為0°和180°時通過短時傅里葉變換后的仿真結果。從圖3(b)可以看出,不加入噪聲時,時頻域中出現余弦包絡、零頻帶和閃爍帶,余弦包絡能量相對較弱。圖3(c)表示脈壓后SNR=10 dB時旋翼回波通過短時傅里葉變換后的仿真結果。從圖3(c)可以看出,當回波中存在噪聲時,信號中的余弦包絡被噪點覆蓋,難以辨識。

圖3 雙葉片散射點模型閃爍現象示意圖Fig.3 Schematic diagram of flashing phenomenon of double blades’ scattering point model
余弦包絡是各個葉片葉尖散射點回波的時頻結果,其微多普勒頻率可表示為

(4)


(5)



圖4 微動參數估計流程圖Fig.4 Flow chart of micro-motion parameters estimation
211 網絡結構
Zhang等提出一種用于圖像去噪的卷積神經網絡,稱為DnCNN神經網絡,該網絡包含17層,如圖5所示。輸入為帶噪灰度圖片=+,其中表示原始圖像,表示噪聲,第1層包含卷積層(Conv)和ReLu層,該卷積層包含64個卷積核且卷積核大小為3×3×1,ReLu代表一種非線性激活函數ReLu()=max(0,),用于避免因網絡層數增加造成的梯度消失等問題。第2層到第16層具有相同結構,包含卷積層(Conv)、批標準化(batch normalization,BN)層和ReLu層,卷積層包含64個卷積核且卷積核大小為3×3×64,在機器學習領域,有一個重要的獨立同分布(IID)假設,即訓練數據與測試數據滿足相同分布,BN層處于卷積層和ReLu層之間,可以使神經網絡每一層的輸入都保持相同的分布特性,較好地保證訓練出的模型在測試集上有一個很好的效果。第17層僅為一個卷積層,包含64個卷積核且卷積核大小為3×3×64。輸出為殘差圖,即中噪聲的分布圖片。整個網絡的損失函數定義為

圖5 DnCNN結構Fig.5 DnCNN structure

(6)

DnCNN最大的特點就是引入了BN層和殘差學習,BN層已經在前文介紹,這里不再贅述,所謂的殘差學習如圖6所示。其思想是不再通過神經網絡擬合潛在的映射(;)=(+;)=,而是通過神經網絡擬合出(;)=(+;)=。依據文獻[22],直接擬合函數(;)的難度更大,而且隨著網絡層數的增加,不僅會出現梯度消失和梯度彌散的問題,網絡模型也會退化,而擬合函數(;)的難度相對而言更小,網絡模型更加魯棒。

圖6 神經網絡映射方式Fig.6 Neural network mapping
2.1.2 網絡訓練方式
去噪和去閃爍網絡通過串行方式進行訓練,如圖7所示。本文不直接利用DnCNN訓練一個既可去噪又可去閃爍的網絡,因為根據文獻[22],輸入的訓練對如果表現特征差距較大,那么網絡將不易訓練。顯然,集合與集合的特征差距相比集合與集合的特征差距要小,所以分別訓練兩種模型。

圖7 網絡訓練方式示意圖Fig.7 Schematic diagram of network training
逆Radon變換是由平面中點的Radon變換推導出的,常用于微動目標參數的估計,其數學原理如下。
如圖8所示平面內任意一點(,)=(-)(-),記為,∠=,為該平面內的一條直線,表達式為

圖8 xoy平面內任意點示意圖Fig.8 Schematic diagram of any point in the xoy plane
=cos+sin
(7)
點沿該直線的Radon變換可表示為

(8)


(9)
式中:=cos;=sin??梢?經過逆Radon變換,余弦曲線聚焦為一點,該點坐標為(,)。







(10)

(11)



(12)

(13)
將以(,)為圓心,半徑為10個單位內的|(,)|置為0,計算|(,)|并確定其坐標(,),由式(10)和式(11)估計出第二個葉片的微動參數。
將以(,)為圓心,半徑為10個單位內的|(,)|置為0,計算|(,)+1|以及||(,)+1|-|(,)||。
(1) ||(,)|-|(,)+1||≤delta·|(,)|,那么求出|(,)+1|的坐標(+1,+1),并依據式(10)和式(11)得出第+1個葉片的微動參數,更新值為=+1,重復步驟3。

利用黃金分割法的頻率估計方法具有估計精度高、迭代次數少、計算量小的優點。下面對估計精度與迭代次數的關系進行分析。


(14)
通過式(1)即可得到黃金分割法所需的迭代次數為

(15)
若利用遍歷搜索,對于轉動頻率的搜索區間[,]和搜索精度,需要的搜索次數為

(16)
若設置頻率搜索區間為[2π,14π],估計精度為0001時,利用式(15)可以得到黃金分割法搜索需要的迭代次數僅為22次。當估計精度要求相同,通過遍歷搜索時,利用式(16)可以得到需要的搜索次數為37 699次??梢钥闯?基于黃金分割法的頻率搜索方法在保證估計精度的同時有效減少了搜索次數,運算速度大大提高。
訓練集和測試集數據在64位Window10系統,英特爾Core i7-8750H 2.2 GHz CPU,內存8 GB設備上,通過Matlab 2019生成;神經網絡的訓練在同樣設備配置下,采用Pytorch框架,GPU版本實現;旋翼目標微動參數的估計在同樣設備配置下,通過Matlab 2019得到。
本文所采用的雷達信號是模擬信號,帶寬=1 MHz、采樣頻率=2 MHz、脈沖重復頻率PRF=8 000 Hz、脈沖寬度=100 μs、觀測時間=0.512 s、載頻=1 GHz、雷達反射截面積RCS=1、光速c=3×10m/s,考慮實際情況,設置脈壓后SNR為5~15 dB之間任意數。在散射點模型下生成旋翼目標的雷達回波并作短時傅里葉變換獲得訓練集、、各400張,測試集40張,每張圖片的像素值為180×180,圖片類型為灰度圖,圖片塊的像素值大小為40×40,每次投入網絡訓練數為128個。為確保無重復參數作用,導致數據集失效,任意數的生成采用洗牌算法,具體仿真參數如表1和表2所示。

表1 微動參數

表2 短時傅里葉變換參數
采用正交初始化的方式設置初始權值,優化器采用Adam優化器。在訓練過程中,前15輪學習率為0.000 1,后15輪學習率為0.000 01。為了使數據得到增強,會隨機對圖片塊作翻轉、平移等操作。
基于深度學習的圖片預處理結果
圖9(a)、圖10(a)和圖11(a)為隨機選取的3張測試集圖片,其雷達直升機旋翼回波模型仿真參數如表3所示。圖9(b)、圖10(b)和圖11(b)分別為圖9(a)、圖10(a)和圖11(a)經過去噪網絡以后的結果。可以看出,去噪網絡有效去除了噪點的影響,恢復出了圖片中潛在的時頻特征,最外側余弦包絡也被很好地保留。圖9(c)、圖10(c)和圖11(c)分別為圖9(b)、圖10(b)和圖11(b)經過去閃爍網絡以后的結果??梢钥闯?閃爍帶和零頻帶被消除,而且余弦包絡能量被增強。

圖9 預處理結果(1)Fig.9 Preprocessing results (1)

圖10 預處理結果(2)Fig.10 Preprocessing results (2)


圖11 預處理結果(3)Fig.11 Preprocessing results (3)

表3 雷達直升機旋翼回波模型仿真參數
常規逆Radon變換的結果
圖12的3個分圖分別表示圖9(a)、圖10(a)和圖11(a)不做任何處理直接通過逆Radon變換得出的結果??梢钥闯?無論旋翼目標的微動參數如何變化,整個圖片僅會中心位置處出現一個聚焦點,因此不對閃爍條件下的時頻結果進行預處理,直接使用逆Radon變換是無法估計出旋翼的微動參數的。

圖12 預處理結果(4)Fig.12 Preprocessing results (4)
基于改進逆Radon變換的參數估計結果


表4 旋翼長度估計結果

圖13 預處理結果(5)Fig.13 Preprocessing results (5)

表5 旋翼目標轉動角頻率估計結果及誤差

表6 旋翼葉片初相估計結果

表7 旋翼葉片數估計結果

表8 不同算法運算速度比較
以旋翼葉片長度的均值、轉動角頻率的誤差、各個葉片初相角的估計值、葉片數量的估計值作為度量標準,從表4~表7的仿真結果可以看出,閃爍現象下,本文方法能夠很好地估計出旋翼目標微動參數,精度較高。對比表8和表9中改進逆Radon變換與傳統逆Radon變換參數估計時間和兩種參數搜索方法的迭代次數,也驗證了改進方法的時效性。

表9 不同算法迭代次數比較
旋翼目標特征參數的提取對機型的穩健識別意義重大。一般地,對閃爍現象下微動參數的提取方法研究不多,為此本文結合深度學習和變換域方法各自的優勢,首先進行時頻圖預處理,然后進行參數估計,主要有以下結論:一是訓練的去噪網絡和去閃爍網絡魯棒性較強,從測試集結果來看,針對不同的時頻圖類型,都能夠顯著消除噪點、零頻帶和閃爍帶,為后續參數估計打下了基礎;二是改進逆Radon變換方法能在較高的精度要求下估計出旋翼目標的微動參數,時效性強。
實際情況中,目標所處的噪聲環境十分復雜,DnCNN網絡的去噪能力和去閃爍能力會受到限制,不利于后續的參數估計。如何構建一種具有去除實際噪聲的網絡模型是下一步研究的重點。