目標微動參數估計的曲線跟蹤算法

李 飛 糾 博 邵長宇 劉宏偉

（西安電子科技大學 雷達信號處理國家重點實驗室，陜西 西安710071）

引 言

目標或目標組成部分的振動或轉動會對雷達回波的頻譜產生調制，在目標多普勒頻率上產生邊帶（微多普勒頻率），美國海軍實驗室的Victor C.Chen教授將這種現象稱為微多普勒效應，Victor C.Chen教授還研究了雷達振動和轉動點目標對單頻信號的微多譜勒效應［1］.Sparr等人分析了APY-6SAR數據，指出振動對雷達回波的相位調制是正弦調制［2］.由于目標微動對雷達回波的頻譜產生調制，所以可以利用調制的多普勒譜估計目標微動參數，為目標分類識別提供有效的目標特征信息［3-5］.因此，目標微多普勒特性研究具有廣泛的工程應用價值.

微多普勒估計等價于多分量非平穩信號的瞬時頻率 估 計 （Instantaneous Frequency Estimation，IFE）問題，現有IFE的方法［6-12］可以分為兩類：基于濾波器的估計方法和基于時頻分析的估計方法.

基于濾波器的IFE算法原理是根據接收信號調整濾波器結構，利用鎖相環（Phase Locked Loop，PLL）跟蹤信號相位的變化.這種方法的缺點是瞬時頻率估計精度受噪聲影響比較大且不能準確跟蹤瞬時頻率快速變化的信號.時頻分析將信號的時間和頻率信息表現在2D平面（時頻域）中，在時頻域的IFE算法是一種有效解決微多普勒估計問題的方法.Henry K.Kwok提出一種基于自適應短時傅里葉變換 （Adaptive Short-Time Fourier Transform，ASTFT）和隱馬爾科夫模型 （Hidden Markov Model，HMM）的IFE算法［6］，該方法解決了單分量信號IFE問題；但是多分量信號的時頻曲線在時頻域中是混疊在一起的，所以多分量信號的IFE問題就是如何將各個分量時頻曲線分離.近期有部分文獻利用通信領域中維特比算法（Viterbi Algorithm，VA）在時頻域中將各個分量信號的時頻曲線分離［7-9］.文獻［7］首先用S-method時頻分析方法得到多分量信號的時頻分布，然后通過VA在時頻域將多條時頻曲線分離得到各個分量信號的IFE；但是VA并未考慮微多普勒的正弦調制，這就使得VA在處理時頻曲線交點時容易產生錯誤.文獻［10］應用霍夫變換（Hough Transformation，HT）估計微動參數，HT在參數域進行多維搜索提取曲線；而HT算法計算復雜度會隨參數估計精度的提高而增加.

針對窄帶目標微多普勒估計即多分量正弦調制信號的IFE問題，提出一種基于曲線跟蹤（Curve Tracking，CT）的IFE算法.該方法將CT算法中的系統模型設定為正弦模型，更有利于跟蹤時頻曲線隱狀態的變化，而且CT算法采用最近鄰數據關聯算法（Near Neighbor Data Association NNDA）與擴展Kalman濾波器對時頻曲線進行分離和平滑，能夠保證算法的實時性.

1 信號模型

本節以旋轉為例分析目標微多普勒，系統幾何示意圖如圖1所示.假設在窄帶條件下，目標回波可簡單表示為

式中：Ai為目標點反射系數；Ri（tm）表示散射點P至雷達的距離為

式中：R0為旋轉中心O到雷達的距離；ν和a分別為雷達與目標之間的相對速度和加速度；ωi表示其旋轉角速度；ri，θi代表散射點P的旋轉半徑和初始相位.

圖1 系統幾何示意圖

由上述分析可得目標多普勒與微多普勒分別為

由式（3）可以看出：目標與雷達相對平動會使多普勒產生平移與線性變化；而微多普勒頻率fdmi（tm）是以目標旋轉頻率為頻率的余弦函數，即散射點P回波信號相位受正弦調制（Sinusoidal Modulation，SM），且最大微多普勒頻率變化量為.當目標進行振動或者擺動等微動形式時，目標回波信號也是正弦調制信號［1］.

運用短時傅里葉變換（Short-Time Fourier Transform，STFT）可以得到信號的時頻分布.STFT的思想就是在傅里葉變換中引入一個窗函數，它的數學定義為

式中：＊代表復數共軛；g（t）是一個窗函數.

目標與雷達相對平動會引起時頻曲線的平移與旋轉.本文假設目標與雷達相對平動參數已知或者已做運動補償，平動補償后回波信號為

對式（6）進行STFT可得時頻分布

式中Tw為窗長.

由回波信號的時頻分布可以提取微多普勒對應的時頻曲線峰值點，得到單分量信號的IFE；當回波信號包含多個分量信號（多個微動散射點）時，不能保證每次提取的峰值點對應同一分量信號的時頻曲線，因此在時頻域進行IFE所面臨的問題是如何將各個分量信號時頻曲線正確分離，第二節將詳細介紹基于CT的IFE算法.

2 基于CT的多分量信號IFE算法

當目標包含多個微動部件時，微多普勒的估計相當于多分量信號的IFE.假設目標有M個旋轉散射點（旋轉中心相同），且相對平動已經補償，則由式（6）可得目標的M個旋轉散射點回波信號為

式中：Ai，ri，ωi，θi分別表示第i個 散射點回波幅度、旋轉半徑、旋轉角速度和初始相位.

提出的基于曲線跟蹤IFE算法是在回波信號的時頻域完成對多分量信號s（tm）的瞬時頻率的估計.由式（7）可得信號時頻譜模值TF（tm，fd）

式中：fdmi（tm）為第i個散射點的微多普勒.通過式（10）利用CLEAN的思想提取峰值得到時頻域的峰值曲線（n），i＝1，2，…，M.

式中：TF0（n，k）＝TF（n，k），TF（n，k）為TF（tm，fd）的離散形式，n，k分別為慢時間采樣標號和頻率采樣標號，n＝1，2，…，N，2Δ＋1為信號時頻譜的寬度，N為時頻譜的慢時間采樣點數.但是由此得到的信號的時頻信息是混疊在一起的，本文采用NNDA算法將各分量信號s（tm）的時頻曲線在時頻域內分離；最后通過傳統的平滑濾波對分離的時頻曲線進行平滑濾波再處理，得到更為精確的IFE.

CT算法包括數據關聯算法和濾波算法.數據關聯算法以濾波算法中的預測數據點與下一時刻各個數據點的統計距離為準則，將各個分量信號的時頻曲線關聯、分離.濾波算法是利用時頻域峰值數據估計各個分量信號時頻曲線的隱狀態，通過隱狀態的變化跟蹤各個時頻曲線，最后對其進行平滑.由式（10）可知各個分量的時頻曲線f＾i（n），i＝1，2，…，M是正弦曲線，所以可將濾波算法的運動模型設為正弦模型，以便更精確地描述時頻曲線隱狀態變化.

2.1 CT濾波、平滑算法

濾波和預測是跟蹤系統中最基本的因素.濾波和預測的目的是估計當前和未來時刻目標的運動狀態.常用的濾波與預測方法有線性自回歸濾波，最小二乘濾波，Kalman濾波，擴展Kalman濾波等.本文所用的濾波算法為擴展Kalman.擴展Kalman是利用線性化技巧將非線性濾波問題轉化為一個近似的線性濾波問題，其中最常用的線性化方法是泰勒級數展開，在此采用一階擴展Kalman濾波算法［13-14］.

本文以旋轉目標為例，由式（9）可知回波信號的時頻曲線按正弦規律變化，采用正弦模型的擴展Kalman能更準確地描述時頻曲線隱狀態變化趨勢.因篇幅有限，不再詳細介紹擴展Kalman濾波平滑算法的數學結構，擴展Kalman濾波器詳細設計請參考文獻［14］，在此僅給出擴展Kalman濾波器的狀態方程、觀測方程和雅克比矩陣.狀態方程為

式中：Δt為采樣間隔；下標k表示時刻k；過程噪聲qk-1為零均值協方差矩陣為Qk-1的高斯噪聲，狀態向量xk為

式中：θk，ωk，ak分別表示k時刻正弦函數的相位、角速度和幅度值.

觀測方程為

式中：yk表示觀測；觀測噪聲rk為零均值方差為σk的高斯噪聲.

雅克比矩陣 Hx（m，k）為

式中：m為狀態均值向量.

2.2 數據關聯算法

數據關聯算法以濾波算法中的預測數據點與下一時刻各個數據點的統計距離為準則，將各個分量信號的時頻曲線關聯、分離.文中采用計算簡單的NNDA算法對各分量信號時頻曲線進行關聯、分離.從第3節實驗結果可以發現該數據關聯算法能夠有效地將時頻域中的多條時頻曲線分離.NNDA算法思想是：以預測狀態得到的觀測為橢圓中心（該橢圓大小由門限γ確定），判定接收的觀測值是否落入該橢圓內.若落入橢圓內的測量值只有一個，則該測量值被用于更新狀態；但若有一個以上回波落在此橢圓內，此時選取統計距離最小的測量值用于更新狀態［13］.NNDA算法準則的數學表達式為

式中νk和Sk分別為新息和新息協方差.

綜上所述，本文所提出的基于CT的多分量信號IFE算法流程圖如圖2所示.

圖2 基于CT的多分量信號IFE算法流程圖

3 仿真實驗

該部分分別用仿真數據與電磁仿真數據驗證本文提出的IFE算法的有效性.

3.1 三分量信號IFE仿真實驗

假設場景中包含三個旋轉散射點，各個散射點微動參數如表1所示.

表1 微動參數

圖3 回波信號時頻分布

對回波信號做STFT，得到目標散射點的時頻分布，如圖3所示.根據式（10）提取出的3條時頻曲線峰值如圖4（a）所示.由圖4（a）可以看出：3條時頻曲線峰值混疊在一起，而且由于時頻曲線一個周期內存在多個交點的影響，使得在曲線交點處提取的峰值點偏差較大，連續性也較差.本文提出的基于CT的IFE算法首先利用數據關聯算法將圖4（a）中混疊在一起的時頻曲線關聯、分離，然后經過平滑濾波處理使分離后的數據點在交點處的連續性得到改善.分離后的時頻曲線峰值點如圖4（b）所示，由此可以看出NNDA算法能夠有效地將各個分量信號的時頻曲線關聯、分離；平滑后的瞬時頻率估值曲線如圖3中曲線所示，而且三分量平滑后的時頻曲線與信號的時頻分布相吻合.

根據平滑后的微多普勒可估計各個分量信號對應的微動參數：旋轉半徑、旋轉頻率和初相，如表2所示.由表2可知：該算法能夠精確估計旋轉頻率，三分量的旋轉半徑和初相估計平均誤差是真實值的2%和5.8%，本文所提出的基于曲線跟蹤的IFE算法能夠用于精確估計多分量信號的瞬時頻率和微動參數.

表2 微動參數估值

圖5給出加性高斯白噪聲背景下，信噪比RSN＝10dB時瞬時頻率分離結果和平滑濾波結果，表3給出了微動參數的估值.信噪比的定義為

由于噪聲影響，提取的峰值點中含有噪點，如圖5（a）所示，本文提出的CT算法仍能有效分離和準確估計時頻曲線，如圖5（b）所示.圖6給出了VA在RSN＝10dB時時頻曲線分離結果，對應的微動參數估計結果如表4所示.由圖6可以看出：由于曲線交點與噪點的影響，VA算法提取的分量1時頻曲線有較大的形變；在時頻曲線前兩個周期內分量2與分量3的幅度較小，而且分量3有兩處頻率發生突變.比較表3與表4，可知在RSN＝10dB時，CT算法與VA算法旋轉頻率估值準確，CT算法對三分量的旋轉半徑和初相估計較VA算法稍好.

圖6 RSN＝10dB VA算法IFE估計結果

表3 RSN＝10dB CT算法微動參數估計結果

表4 RSN＝10dB VA算法微動參數估計結果

當HT算法參數搜索間隔較小時，估計的微動參數如表5所示，增大HT算法搜索間隔，微動參數估值如表6所示.由于條件Ⅰ中參數搜索間隔較條件Ⅱ小，所以條件Ⅰ下微動參數估計精度高，但是耗時較長約10s；條件Ⅱ耗時雖少約1.6s，但是半徑估計精度變差.而對同樣的數據處理，本文在保證參數估計精度前提下提高了算法的實時性，耗時約1.2s.

表5 RSN＝10dB HT微動參數估計Ⅰ

表6 RSN＝10dB HT微動參數估計Ⅱ

為了充分驗證算法的抗噪聲性能，圖7給出400次蒙特卡洛實驗三分量的旋轉半徑估計均值與標準差隨信噪比變化曲線.由圖7可知三分量的旋轉半徑估計均值與標準差隨著信噪比的增加得到改善.當RSN＜5dB時，由于提取峰值中噪點的原因，導致該算法無法正確提取時頻曲線.當5dB＜RSN＜10dB時，分量3的估計均值較真值小，這是由于CT算法受噪聲影響，正弦曲線峰值偏小引起的.當RSN≥10dB時，算法受噪聲影響逐漸減弱，各個分量的半徑估計均值趨于真值，且估計的標準差進一步降低.蒙特卡洛實驗結果說明本文算法在RSN≥10dB時能夠估計各個分量的微動參數.

圖7 旋轉半徑估計均值與標準差

3.2 電磁仿真信號IFE仿真實驗

該小節旨在利用電磁仿真數據驗證本文所提出的基于CT的IFE算法的有效性.錐體目標尺寸如圖8實線所示，錐體高0.96m，錐底半徑0.25m.載頻fc＝9GHz，雷達脈沖重頻fPR＝500Hz，雷達視線俯仰角45°，方位角90°；錐體擺動角為10°，擺動頻率1Hz.

圖8 錐體目標尺寸圖

在雷達視線內存在兩個微動散射點，利用STFT對回波數據進行時頻分析，可得目標微多普勒時頻譜，回波時頻譜以及CT算法提取的時頻曲線如圖9所示.圖9中提取的時頻曲線與信號時頻曲線重合說明本文算法能夠精確地提取目標微多普勒.由提取的微多普勒估計目標散射點的微動周期：T1＝1s，T2＝1s.

圖9 無噪聲時頻分布與IFE比較

加性高斯白噪聲背景下，RSN＝5dB時IFE結果與時頻分布比較如圖10所示.對最終的平滑濾波后的瞬時頻率估值做傅里葉變換，可以得到目標散射點的微動周期：T1＝1s，T2＝1s.實驗結果說明本文提出的算法能夠有效將各個分量信號的時頻曲線分離，通過濾波平滑精確估計瞬時頻率，并且具有一定的抗噪聲性能.

圖10 RSN＝5dB時頻分布與IFE比較

4 結 論

目標微動參數估計的研究受到越來越多的關注.本文通過分析微動對目標回波的影響，說明振動、旋轉以及擺動等微動對回波的調制為正弦調制；應用STFT得到目標的時頻分布，以此為基礎通過CT算法將各個分量的時頻曲線分離，并通過濾波平滑得到精確的瞬時頻率估值，進而估計目標的微動參數.本文針對簡單微動（旋轉、振動以及擺動等）將CT算法中的系統模型設定為正弦模型；當散射點微動形式復雜時，如何建立系統模型并能準確估計微多普勒和微動參數是本文的后續工作.

［1］CHEN V C，LI F Y，HO S S.Micro-Doppler effect in radar phenomenon，model and simulation study ［J］.IEEE Trans Aerospace and Electronic System，2006，42（1）：2-21.

［2］SPARR T，KRANE B.Time-frequency analysis of vibrating targets in airborne SAR systems［J］.IET Radar，Sonar ＆ Navigation，2003，150（3）：173-176.

［3］楊立明，曹祥玉.直升機旋翼對回波的調制效應分析［J］.電波科學學報，2002，17（1）：93-96.YANG Liming，CAO Xiangyu.Analysis of the backscattered wave from a helicopter rotor［J］.Chinese Journal of Radio Science，2002，17（1）：93-96.（in Chinese）.

［4］高紅衛，謝良貴，文樹梁，等.基于微多普勒特征的真假目標雷達識別研究［J］.電波科學學報，2008，23（4）：775-780.GAO Hongwei，XIE Lianggui，WEN Shuliang，et al.Research on radar target identification ofwarhead and decoys based on micro-Doppler signature［J］.Chinese Journal of Radio Science，2008，23（4）：775-780.（in Chinese）

［5］李彥兵，杜 蘭，劉宏偉，等.基于信號特征譜的地面運動目標分類［J］.電波科學學報，2011，26（4）：641-648.LI Yanbing，DU Lan，LIU Hongwei，et al.Moving ground targets classification based on eigenvalue sequence of echo signal［J］.Journal of Electronisc＆Information Technology，2011，26（4）：641-648.（in Chinese）

［6］KWOK H K.JONES D L.Improved instantaneous frequency estimation using an adaptive short-time Fourier transform［J］.IEEE Trans Signal Processing，2000，48（10）：2964-2972.

［7］THAYAPARAN T，STANOVICB L，DJUR＇OVICB I，et al.Micro-Doppler-based target detection and feature extraction in indoor and outdoor environments［J］.Franklin Inst，2008，345（6）：700-722.

［8］STANK＇OVIC L，DJUR＇OVICI，OHSUMI A，et al.Instantaneous frequency estimation by using Wigner distribution and Viterbi algorithm ［C］// ICASSP.Hong Kong，April 6-10，2003：121-124.

［9］DJUR＇OVICI.Viterbi algorithm for chirp-rate and instantaneous frequency estimation［J］.Franklin Inst，2011，345（6）：700-722.

［10］ZHANG Q，YEO T S，TAN H S，et al.Imaging of a moving target with rotating parts based on the Hough transform［J］.IEEE Trans Geosci Remote Sens，2008，46（1）：291-299.

［11］SUCIC V，SAULIG N，BOASHASH B，et al.Estimating the number of components of a multicomponent nonstaionary signal using the short-term timefrequency Rényi entropy［J］.EURASIP Journal on Advances in Signal Processing，2011：125.doi：10.1186/1687-6180-2011-125.

［12］BOASHASH B.Estimating and interpreting the instantaneous frequency of a signal-part II：Algorithms and Applications［J］.Proc IEEE，1992，80（4）：521-538.

［13］HUE C，CADRE J，PEREZ P.Sequential Monte Carlo methods for multiple target tracking and data fusion［J］.IEEE Trans Signal Processing，2002，50（2）：309-325.

［14］BAR-SHALOM Y，LI X R，KIRUBARAJAN T.Estimation with Applications to Tracking and Navigation［M］.New York：John Wiley and Sons，2001.