基于連續波雷達的人體目標微多普勒頻率估計方法

丁一鵬，柳潤金，許雪梅

(中南大學物理與電子學院，湖南長沙，410083)

連續波多普勒雷達在人體目標探測方面[1-2]具有獨特優勢，它可以抑制靜止雜波，穿透墻壁等障礙物，且檢測條件不受天氣、光線影響，因此被廣泛應用于反恐、邊境安防、災難救援等重要領域。微多普勒頻率作為雷達目標所具有的獨特特征，可用于目標運動狀態的識別[3-4]。因此，微多普勒頻率的精確估計對提高雷達探測和目標識別的準確性至關重要。為了更精確地識別感興趣目標的散射部位(例如手臂和腿)的運動，需要先提取這些部位的基本回波分量。然而，人體柔性關節和復雜的運動模式使得其雷達回波動態范圍大，時頻譜帶相互交疊，人體各部位的瞬時微多普勒頻率很難被精確估計。此外，人體各部位微多普勒回波為非平穩隨機信號，如何準確擬合回波局部特性并精確估計目標微多普勒頻率，是檢測算法設計中需要考慮的問題。

學者們針對微多普勒頻率估計技術進行了大量研究。最常用的分析方法是線性時頻分析技術，如短時傅里葉變換[5-6]和S變換[7]，這些技術在分析處理多分量的信號時可免受交叉項的干擾且易于實現，但時頻分辨率的折中處理使其在雷達目標探測領域的應用受到限制。傳統的非參數信號分離算法如經驗模式分解(EMD)[8-9]、主成分分析(PCA)[10]和Hilbert-Huang 變換[11]等都是基于回波頻率特性而不是目標散射特性，這會給分離后的頻率估計結果帶來很大的偏差。參數信號分離算法[12]的參數搜索過程則通常涉及大量計算。稀疏分解算法[13]近年來在語音信號、地震信號和雷達信號的分離和去噪研究中也取得了較大進展，然而，算法所需原子庫的完備性和原子匹配會對算法的收斂性造成影響。借鑒圖論中的邊緣檢測算法思想，利用維特比算法[14]可以對瞬時頻率軌跡進行跟蹤，但該算法只能準確估計單分量信號頻率。改進的維特比算法[15]增加了局部斜率差懲罰函數，然而，在多散射點條件下該算法仍表現出穩健性不足、估計精度較低的特點。結合航跡關聯思想，通過設置隸屬度懲罰函數[16]可以提高改進型維特比算法對多分量信號頻率估計的精度，但在人體目標特定散射部位的微多普勒頻率估計中會出現路徑分叉問題。

本文作者基于三次動態指數平滑預測方法獲取時頻分布中最佳路徑的非線性變化趨勢，將其表示為新型懲罰函數的形式并對傳統維特比算法進行改進，從而增強最佳路徑的識別能力；為了更好地擬合人體散射部位的局部特性，借鑒Hough變換的頻率估計思想[17]，用動態平滑系數代替傳統指數平滑法的靜態平滑系數；在人體微多普勒頻率估計的應用場景中，該方法以歷史路徑為參考，并通過動態確定候選路徑的搜索范圍以提高最佳路徑的尋找效率。

1 維特比瞬時頻率估計算法

維特比瞬時頻率估計算法的基本思想來源于圖論和數字圖像處理中的邊緣跟蹤問題[14]。對信號進行時頻分析后可以得到二維時頻分布，利用維特比算法從時頻分布中搜索最佳路徑從而提取信號瞬時頻率主要基于以下2個假設：1)每個時刻的瞬時頻率對應的時頻點的幅值要盡可能大；2)兩相鄰時刻的瞬時頻率變化不會太大，瞬時頻率曲線比較平滑。

對信號進行時頻分析后，可得到N行M列的時頻分布，將時頻分布數據中時刻m的時頻幅值T排列成非遞增序列：

式中：j= 1,2,…,N；m= 1,2,…,M；kj為m時刻的第j個時頻點。

采用如下算式定義懲罰函數f(·)：

考慮分別位于相鄰時刻m和m+ 1的2個時頻點(m,kl)和(m,kj)， 其 中，kl= 1,2,…,N且kj=1,2,…,N，則懲罰函數g(·)可定義為

式中：c1為懲罰因子，可以反映g(·)懲罰函數值對總懲罰值的貢獻度；Δ1為2個連續點之間瞬時頻率變化的最大期望值。路徑的總懲罰值計算方法如下：

式中：p為總懲罰值。最終的最優路徑即為使懲罰值最小的時頻曲線Kmin：

為了對多分量信號瞬時頻率進行估計且增強算法穩健性，改進的維特比算法[16]設置了如下隸屬度懲罰函數：

式中：f′k,i為通過多點多步預測的瞬時頻率；fk,i為需要進行隸屬度判別的瞬時頻率；u為常數。則式(5)所示路徑總懲罰值可表示為

傳統維特比算法在單分量信號和頻譜混疊不嚴重的多分量信號瞬時頻率估計方面表現良好，然而，傳統維特比算法是基于時頻分布中脊線的幅值定義f(·)懲罰函數，其所定義的g(·)懲罰函數僅通過限制2個連續采樣點絕對頻率變化值來提取路徑最短的平滑脊線，而不考慮脊線的變化趨勢。由于采樣間隔太小，懲罰函數的約束是較弱的局部約束[18]，在面臨時頻分量交疊處的路徑搜索時，會產生路徑分叉問題。傳統維特比算法路徑分叉示意圖如圖1 所示，其中Im(t)和In(t)分別為時頻域交疊的2 個分量信號，gm(k4,k5)為由Im(t)分量上的點4 和5 的頻率變化而引起的g(·)懲罰值；gn(k4,k8)為由In(t)分量上的點4和8的頻率變化而引起的g(·)懲罰值。對于待估分量Im(t)，點1至點7為最佳路徑。 然而， 當gm(k4,k5)≈gn(k4,k8)且點8的幅值比點5的幅值更大時，傳統維特比算法會將干擾分量In(t)上的點8 誤選為待估分量Im(t)的路徑，后續的路徑搜索可能完全偏離待估分量，從而將點1至點4和8至點10誤搜索為最佳路徑，產生路徑分叉問題。在隸屬度懲罰函數[14]的設置過程中，所采用的預測方法是選擇局部若干時刻的瞬時頻率對當前時刻瞬時頻率進行線性預測而不是基于全部歷史數據，故該方法不適用于非平穩信號序列的預測。此外，在整個時頻平面范圍內搜索候選路徑會使算法計算量增大。

為解決上述問題，提高微多普勒頻率估計精度和最佳路徑的搜索效率，本文采用一種基于三次動態指數平滑預測的具有高搜索效率的新型維特比算法用于估計人體目標微多普勒頻率。

2 基于三次動態指數平滑預測的新型維特比算法用于人體微多普勒頻率估計

2.1 人體步態微多普勒回波模型

為了保持實驗模型所固有的動態特性，采用基于簡單運動學方法設計的經典Boulic人體模型[19]進行人體微多普勒回波分析。當連續波多普勒雷達發射天線發射電磁載波時，運動的目標對載波信號產生頻率調制，接收天線接收到的人體基帶回波模型r(t)可用式(9)表示，Boulic 模型將人體目標分為若干散射部分，回波信號r(t)即人體各散射部位(如頭部、軀干和腿部)的回波之和。

式中：Nh為人體散射部位的數量；c和fc分別為光速和雷達載波頻率；θi為目標運動速度與雷達視線方向的夾角；Ri(0)為第i個散射中心與雷達的初始視距；vi(t)為第i個散射部位運動速度。第i個散射部位的微多普勒頻率fdi(t)可以表示為

微多普勒頻率是由人體運動模式和柔性關節的復雜相互作用對雷達回波產生的非線性頻率，在時頻譜上表現為整體多普勒頻率周圍的旁瓣和邊帶效應。

2.2 微多普勒頻率估計算法

人體正常行走時具有大慣量、低速度的特點，且其微多普勒頻率曲線為類正弦曲線?；诖?，本文考慮將可靠性和實時性高的三次指數平滑預測技術應用于人體目標微多普勒頻率估計。下面介紹基于動態指數平滑預測模型的新型維特比算法的原理和構建步驟。

三次指數平滑預測是一種序列分析和預測方法，能將非線性時間數據序列的數量差異抽象化，對序列進行修勻從而消除不規則和隨機擾動，顯示出預測對象變動的基本趨勢。利用傳統維特比算法在短時傅里葉變換后的時頻分布中搜索初始路徑，并計算得到三次指數平滑預測的歷史路徑平滑值：

獲取歷史路徑的平滑值后，可得三次指數平滑法的預測系數λ，γ和η的計算公式[20]如下：

為增強傳統維特比算法懲罰函數的最佳路徑選擇能力，減少頻率模糊區域的微多普勒瞬時頻率分叉問題，將新型維特比算法的懲罰函數μ(·)定義如下：

式中：ε為三次指數平滑基于全部歷史數據平滑值的短期預測值，且ε=λ+γ+η；Δ2為設定的頻率閾值，在人體步態微多普勒頻率估計的應用場景中，Δ2一般設置為3～5 Hz；c2為懲罰因子，用于表征μ(·)懲罰函數對總懲罰值的貢獻度，當預測值與當前估計值距離較小時，懲罰值較小，該預測值所在路徑屬于最佳路徑的可能性較大；若路徑在頻率模糊區域向其他路徑分叉，遠遠偏離預測路徑，則將產生較大的μ(·)懲罰值，使之被選為最佳路徑的可能性減小。因此，這樣定義的罰函數既考慮了時頻脊線的變化趨勢，也加強了確定最佳路徑需要的整體約束，可對當前搜索路徑進行修正。

將懲罰函數μ(mi,kj)代入式(5)，有

以p(mi+1,kj)為總懲罰值對時頻分布進行搜索，則以αml為模型參數的微多普勒頻率估計模型F(αml,t)表達式如下：

下面根據Hough變換頻率估計的思想[17]確定動態平滑系數αml。調整平滑系數αml，以F(αml,t)為頻率估計模型對雷達回波信號r(t)進行解調，并將傅里葉變換應用于解調后的信號，得到如下表達式：

在理想情況下，若待估計分量的微多普勒頻率fp(t)≈F(αml,t)，則解調后該分量將成為恒定信號，其能量將以脈沖的形式在頻域中累積。對于估計模型以外的其他分量，其能量則不會在頻域中收斂，因此，當S(f)的信號能量收斂到基帶的最大值時，可以得到感興趣的人體散射部位的微多普勒頻率fp'(t)≈F(α′ml,t)，α′ml即為最佳模型參數?；谥笖灯交A測模型的新型維特比算法構建流程如圖2所示。

將窄帶濾波器應用于具有最佳模型參數的S(α′ml,f)，將當前所估計分量從原始信號中分離出來后，可對其他部位的微多普勒頻率進行估計：

式中：s(t)為被解調分量；FB(f)為濾波器的頻率響應；IFFT為逆傅里葉算子；sp(t)為由F(α′ml,t)分離的分量。

在路徑搜索過程中，并不需要搜索整個時頻平面上的所有點，只有在能量收斂區域附近的搜索才有價值，其他搜索點會耗費更多處理時間。由于正常人體行走時的運動為慣性運動，其雷達回波具有局部靜態特性，且兩相鄰時刻的瞬時頻率變化不會非常劇烈，瞬時頻率曲線應該比較平滑。因此，通過式(6)確定某一時刻的歷史路徑Kpi后，以[Kpi-δ,Kpi+δ]mi作為參考間隔，[Kpi-δ,Kpi+δ]mi+1為搜索區間來進行候選路徑的動態搜索，參考區間和搜索區間長度均為2δ+ 1，若時頻分布中的極差為R，則δ可估計為R/2。在一般情況下，人體正常行走時微多普勒頻率估計的δ不超過5，則在確定從mi到mi+1的路徑時，此搜索方法只需要不超過(2δ+ 1)2即121次搜索。

新型維特比算法用于人體特定散射部位微多普勒頻率估計的步驟如下：在三次動態指數平滑預測的基礎上獲取最佳路徑的非線性變化趨勢，并依據g(·)懲罰函數形式定義新懲罰函數μ(·)，對當前搜索的最佳路徑進行修正從而提高傳統懲罰函數的最佳路徑識別能力，抑制傳統維特比算法在頻率模糊區域的路徑分叉，同時動態調整候選路徑的搜索范圍，提高最優路徑的尋找效率。

3 仿真實驗與分析

為了驗證所提出的基于三次動態指數平滑預測的維特比算法用于人體目標微多普勒頻率估計的有效性和優越性，采用經典的Boulic人體模型作為分析標準，通過MATLAB 軟件進行仿真實驗，對比不同算法的相應估計結果，并進行分析。所應用的Boulic 模型中的人體不同部位的尺寸見表1，這些數據是通過對1 個身高約1.78 m 的男性目標進行實驗測量獲得的。雷達探測場景示意圖如圖3(a)所示，連續波雷達所使用的載波頻率為2.4 GHz，發射單元和接收單元間隔距離約為6.5 cm。目標以0.5 m/s 的恒定徑向速度向雷達走去，每0.005 s采集1次回波信號，接收回波的頻譜如圖3(b)所示。

表1 人體目標各部位尺寸Table 1 Dimensions of different parts for human targets

在基于微多普勒頻率對目標進行識別和分類的實際應用中，目標四肢的微多普勒信息是區分人類和其他動物的典型特征。由于四肢運動的動態范圍大且其微多普勒頻率形狀不規則，最難精確估計，下面分別采用基于短時傅里葉變換的峰值檢測法、傳統維特比算法和本文所提出算法依次對人體左腿和左臂的微多普勒頻率進行估計，結果分別如圖4和圖5所示。在傳統算法的頻率估計過程中，其他部位的雷達回波會對估計結果產生較大干擾，產生路徑分叉問題，故在圖4 和圖5中用細實線標示出其他干擾散射部位的實際微多普勒頻率以用于對比。

由圖4(a)和圖5(a)可見，無論是否處于頻譜交疊區，峰值檢測法所估計的左腿微多普勒瞬時頻率曲線都顯示出來自其他分量的干擾。這是因為峰值檢測法按能量脊線提取頻率，而干擾分量的能量大于感興趣部位分量的能量，使得所估計的頻率曲線往其他分量的頻率曲線跳變。

由圖4(b)和圖5(b)可見：與峰值檢測法相比，傳統維特比算法在一定程度上減小了微多普勒頻率估計誤差。然而，該方法在頻率模糊區域仍存在路徑分叉問題。這是因為在這些區域能量分布密集，傳統維特比算法的f(·)懲罰函數值相近，而g(·)通過僅限制2個連續采樣點絕對頻率變化值來提取路徑最短的平滑脊線，不考慮脊線的變化趨勢，導致維特比算法搜索的路徑分叉至干擾分量頻率曲線上懲罰值最小的點。

由圖4(c)和圖5(c)可見：新型維特比算法所得結果與實際目標微多普勒頻率幾乎相同。不同算法所得手臂微多普勒頻率的均方根誤差和處理時間對比見表2。由表2 可見，與峰值檢測法和傳統維特比算法相比，本文所提方法的均方根誤差顯著降低。

表2 不同算法所得微多普勒頻率的均方根誤差及處理時間對比Table 2 Comparison of micro-Doppler frequency root mean square error and processing time of different algorithms

4 結論

1)基于三次指數平滑預測獲取傳統維特比算法歷史路徑的變化趨勢，并定義新型維特比算法的懲罰函數μ(·)，可增強傳統維特比算法懲罰函數的最佳路徑選擇能力，抑制路徑分叉，提高連續波雷達頻率模糊區域的微多普勒瞬時頻率估計精度。

2)采用Hough變換頻率估計思想，以F(αml,t)為頻率估計模型對雷達回波信號進行解調，依據雷達回波的局部特性動態調整平滑系數αml；同時，動態調整維特比算法候選路徑的搜索范圍，提高了最優路徑的尋找效率。

3)與基于短時傅里葉變換的峰值檢測算法和傳統維特比算法相比，基于三次指數平滑預測的新型維特比算法得到的人體正常行走時左臂微多普勒頻率估計值的均方根誤差顯著降低。