基于貝葉斯濾波算法的警戒相控陣雷達目標跟蹤時間資源優化分配算法

李紀三， 侯嬌， 劉溶， 任淵

(南京船舶雷達研究所， 江蘇 南京 210003)

0 引言

相控陣雷達通過改變信號單元發射的電磁波相位來實現不同方向的波束指向，在計算機控制下，在方位和仰角上可達到微秒量級的波束捷變，具有掃描方式靈活、工作波形豐富、控制參數可變、工作模式可選等優點[1-2]，通過對相控陣雷達的優化控制，包括任務合理規劃、資源有效調度以及工作參數優化選擇等，可實現多功能多任務并能釋放雷達作戰潛能[3-5]。

雷達的主要任務是在各種電磁環境下探測、跟蹤和識別空中及海面目標，其中跟蹤任務的時間資源分配是資源調度的重要方面，跟蹤任務消耗的時間資源與跟蹤數據率呈正比，跟蹤數據率是指單位時間內對目標的采樣次數，是采樣時間間隔的倒數。跟蹤的數據率影響跟蹤精度和穩定性，數據率越高，跟蹤的穩定性和跟蹤精度越高，但是消耗的時間資源也越多。跟蹤任務的時間資源分配是指在保證一定精度和穩定性的前提下，使得跟蹤花費的時間資源最少。根據當前跟蹤的狀況，自適應地確定目標的跟蹤數據率是優化分配時間資源的有效途徑。文獻[6]介紹了基于跟蹤殘差自適應數據率跟蹤算法，該算法通過殘差與量測噪聲方差間的比值對下一時刻采樣周期進行調整，當目標機動時該算法可自適應地減小采樣周期，但是當目標無機動性時采樣周期卻會不斷增大。van Keuk等[7]研究了采樣周期與信噪比、跟蹤穩定度以及資源消耗的關系，提出了基于跟蹤誤差、相對跟蹤精度和波束寬度的采樣間隔計算公式。上面兩個文獻奠定了后續研究基礎。Benoudnine等[8]基于交互多模型濾波算法實現了自適應采樣間隔計算，下個采樣間隔是各模型計算采樣間隔的加權平均。文獻[9]引入了目標威脅等級，作為調整數據率的參考值。文獻[10]將van Keuk等[7]的自適應準則擴展到擴展卡爾曼濾波算法、不敏卡爾曼濾波算法、交互多模型濾波算法和粒子濾波算法。

以上研究把跟蹤精度作為調整數據率的依據和準則，隨著相控陣雷達體制的發展，和差測角波束因其具有實現簡單和精度高的特點已得到廣泛應用，從濾波算法角度，量測誤差相對于傳統機掃雷達大大減小，跟蹤的主要矛盾變成跟蹤穩定性，即跟蹤波束能否罩住目標。

基于殘差的自適應跟蹤算法具有運算量小、工程實現簡單等優點，該算法為了實時響應目標的機動，通常利用當前采樣點的殘差計算下個周期的采樣間隔。在導彈攻擊或者飛機逃離過程中會采取多種樣式的機動方式，在蛇形機動反向機動階段殘差會下降到很小的值，此時會誤判為機動結束而增加跟蹤數據率，進而造成目標逃出雷達的跟蹤波門，造成失跟。

相對于傳統基于殘差的變數據率跟蹤算法把機動引起的殘差看作一個確定變量，本文基于貝葉斯濾波算法，把殘差看成隨機變量，并用相應的概率密度加以描述。在目標勻速直線運動階段，殘差服從正態分布，通過對殘差進行拉依達準則判決來抑制量測噪聲起伏產生的數據率抖動問題。在機動階段，殘差服從修正的瑞利分布，對數據率采用指數平滑法進行平滑，通過在數據率增加和減小時段采用不同的平滑因子，在數據率下降時采用較大的平滑因子，即在計算下個采樣間隔時本周期的殘差中所占比重較大，降低噪聲引起的機動誤判和快速響應目標的機動，在數據率上升時采用較小的平滑因子，即歷史殘差所占比重大，從而減少誤判機動結束所造成的目標失跟。

1 貝葉斯濾波原理

貝葉斯原理是利用已知信息來構造系統狀態變量的后驗概率密度，即用系統狀態轉移方程來計算狀態的先驗概率密度，再利用觀測值對先驗概率密度函數進行修正，從而得到狀態變量的后驗概率密度[11]。

設隨機、離散系統S的狀態空間模型為

(1)

式中：Xk為k時刻(k≥1)的系統狀態向量；F(·)為系統狀態轉移模型；Uk為系統隨機噪聲；Zk為k時刻的系統觀測向量；H(·)為系統觀測模型；Vk為隨機觀測噪聲。

求解狀態后驗分布P(Xk|Zk)是貝葉斯濾波的關鍵問題，可由如下貝葉斯濾波方程遞推求解[12-14]：

1)預測方程(Chapman-Kolmogorov方程)

(2)

2)更新方程

(3)

若系統S是線性的，系統噪聲Uk和觀測噪聲Vk為零均值高斯白噪聲(即服從高斯分布)，其協方差矩陣記為Var(U)=Q，Var(V)=R，則基于貝葉斯濾波算法的卡爾曼濾波算法[15]可表述為

步驟1一步預測：

(4)

Pk|k-1=APk-1|k-1AT+Q，

(5)

步驟2觀測更新：

Kk=Pk|k-1CT(CPk|k-1CT+R)-1，

(6)

(7)

Pk|k=Pk|k-1-KkCPk|k-1,

(8)

經考察得知，Pk|k-1的變化趨勢是非負遞減的，當k趨于無窮時，系統達到穩態存在極限：

(9)

式中：P為協方差矩陣；K為增益矩陣；M為預測方差矩陣。

當目標運動方程采用常速度模型時，得到穩態下的常增益濾波器為

(10)

(11)

2 常增益濾波器的殘差計算

1)一個加速度為u的勻加速模型

X(k)=AX(k-1)+Gu，

(12)

2)測量方程

Y(k)=CX(k)+ω(k),

(13)

式中：Y(k)為k時刻的量測值；ω(k)為量測噪聲，并滿足高斯分布。

3)濾波方程

(14)

4)狀態預測方程

(15)

5)量測預測方程

(16)

將濾波的誤差定義為真值減去估計值，有

(17)

于是，當系統穩定之后，k時刻的誤差應該與k-1時刻的誤差相同，即

(18)

(19)

對于勻加速直線模型，不考慮量測噪聲ω(k)下，由純加速度u引起的位置誤差為

(20)

同理可推導出由u引起的位置殘差[16]為

e(k)=T2u/β.

(21)

3 自適應數據率算法

3.1 基于殘差的采樣周期計算

如果誤差在加速度增加時要自適應地保持恒定，則由(20)式就得出采樣周期應按與加速度的平方根呈反比的方式減小，文獻[6]給出了周期遞推的計算公式：

(22)

式中：D(k)為跟蹤的數據率，是采樣周期T(k)的倒數；e0(·)為歸一化殘差。(22)式可展開為

(23)

跟蹤殘差包含系統的過程噪聲(也稱為模型噪聲)和量測噪聲，過程噪聲是由選擇的目標運動動態方程與實際目標運動過程不符產生的，例如目標當前在做圓周運動，而用直線運動方程來表述，則會產生誤差。

在目標做勻速直線運動、沒有機動情況時過程噪聲為0，殘差主要是由量測噪聲貢獻的，由于量測噪聲服從正態分布，e0(1)，…，e0(k)相互獨立且服從正態分布N(μ,σ2)，密度函數[17]為

(24)

式中：x為跟蹤殘差；σ為標準差；μ為均值。

圖1所示為在目標做勻速運動時按照(22)式直接計算出的數據率。從圖1中可觀察到，由于量測噪聲的正態分布特性，噪聲的起伏引起虛警(目標沒有機動，誤判為目標機動，從而增加數據率)。

圖1 勻速直線運動的數據率Fig.1 Data rate of uniform linear motion

針對圖1由觀察噪聲引起的數據率起伏問題，本文引入統計決策中2σ準則[18]:

P(μ-2σ

(25)

式中：P(·)為概率密度函數。由(25)式可知，如果當前歸一化殘差小于2σ，則認為殘差是由量測噪聲起伏引起的，數據率保持不變；如果當前歸一化殘差大于2σ，則認為殘差是由目標的機動引起的，需要根據殘差調整數據率。

3.2 目標機動模型

目標的機動可由以下模型進行描述：勻速模型、勻加速模型、Jerk模型、Singer模型和當前統計模型等。在這些模型中，當前統計模型較為典型，該模型假定加速度服從修正的瑞利分布，其均值是當前加速度，其均值與方差之間的關系可以用來建立機動加速度的均值和方差自適應算法。修正的瑞利密度函數[19]可描述如下。

1)當目標的當前加速度為正時，概率密度函數為

(26)

式中：umax為已知目標加速度的上限，umax>0；u的均值和方差為

(27)

(28)

2)當目標的當前加速度為負時，概率密度函數為

(29)

式中：u-max為已知目標加速度的負下限，u-max<0；u的均值和方差為

(30)

(31)

3)當目標在機動時，加速度u為修正的瑞利分布，分布的起伏特性決定了加速度引起的殘差T2u/β(見(21)式)，也服從瑞利分布。如果通過(22)式直接計算數據率，則會引起數據率的劇烈振蕩，不利于穩定跟蹤，特別是在機動階段的反向機動時刻，即由正負加速度相互切換時會出現殘差為0的情況，此時按照公式要增加采樣間隔，會引起目標的失跟，為此利用指數平滑法對數據率進行平滑。

一次指數平滑法[20]的公式為

Ds,k=ρDk+(1-ρ)Ds,k-1,

(32)

式中：ρ為平滑系數；Ds,k為時刻k的平滑值；Dk為時刻k的實際值；Ds,k-1為時刻k-1的平滑值。(22)式代入(32)式，有

(33)

(34)

從(34)式中可以看出，當前的數據率與第1次的數據率和殘差序列以及平滑系數ρ有關。k時刻計算數據率時，ρ越大則當前的殘差在計算中占比越大，ρ越小則以前的殘差占比越大。因此可以通過實時調整ρ的取值來調節當前殘差占的比重：在機動開始時刻(即殘差突然變大)采用較大的值，使得系統響應迅速；在殘差變小時采用較小的值，防止機動未結束、誤判為結束而采用較小數據率，造成目標的失跟。

3.3 變數據率跟蹤算法流程

本文算法是在直角坐標系中進行殘差計算，雷達的量測在球坐標系中進行，轉換到直角坐標系時存在非線性問題。若利用過程噪聲對殘差進行歸一化，則過程噪聲計算誤差較大，而文獻[6]介紹的算法利用0.1 n mile的常量對殘差進行歸一化會引起較大的誤差。本文算法在直角坐標系下濾波，在量測坐標系下計算殘差，利用雷達本身的量測噪聲對殘差進行歸一化，具體實施步驟[21]如下：

步驟1對k時刻狀態進行預測:

(35)

式中：xp(k)、yp(k)、zp(k)分別為x軸、y軸、z軸的位置預測值；xs(k)、ys(k)、zs(k)分別為x軸、y軸、z軸的位置平滑值；vs，x(k)、vs，y(k)、vs，z(k)分別為x軸、y軸、z軸的速度平滑值；T(k)為時刻t(k)與t(k+1)間的采樣周期。

步驟2直角坐標系下的預測值轉換為球坐標系：

(36)

式中：rp(k)、βp(k)、θp(k)分別為k時刻預測值在極坐標下的距離、仰角和方位。

步驟3雷達在(βp(k),θp(k))照射后得到量測值(rm(k),βm(k),θm(k))，將量測值轉為直角坐標系[21]：

(37)

式中：xm(k)、ym(k)、zm(k)分別為x軸、y軸、z軸在k時刻的量測值。

步驟4在直角坐標系下對k時刻目標的位置進行濾波：

(38)

步驟5根據k時刻的量測值和預測值計算k時刻的殘差，殘差e(k)定義為t(k)時刻測量值與預測值之差，即

(39)

式中:er(k)、eβ(k)、eθ(k)分別為距離殘差、方位殘差和仰角殘差。

步驟6殘差歸一化：

e0(k)=max(eβ(k)/σβ,eθ(k)/σθ)，

(40)

式中:σβ、σθ分別為仰角量測標準差和方位量測標準差。

步驟7根據k時刻的殘差計算k+1時刻的采樣周期，計算下個周期的數據率

(41)

步驟8對計算的數據率進行平滑，

D(k)=ρD(k)+(1-ρ)D(k-1)，

(42)

式中：

(43)

(43)式表明：當e0(k)>1時，要殘差變大需要降數據率，此時平滑因子取值較大，能夠較快地響應目標的機動，防止目標失跟；當e0(k)=1時，保持數據率不變；當e0(k)<1時，要增加數據率，此時選取的平滑因子較小。為了防止誤判，即機動沒有結束而增大數據率造成失跟，此時采用這種策略會在某種程度上犧牲一些時間資源，保證了目標的穩定跟蹤。

4 仿真實驗

為驗證本文所提算法的有效性，與常規天線周期作為數據率的定數據率Cohen跟蹤算法作比較，利用300次蒙特卡洛算法計算得到的均方根誤差(RMSE)作為濾波性能的衡量標準，對比分析兩種算法在目標作勻速直線運動和機動情況下的跟蹤性能。

4.1 實驗1

仿真的雷達參數為：雷達位于坐標原點(0 km，0 km)處，量測精度為距離60 m，方位0.25°，仰角0.25°.

仿真場景如圖2所示，飛機從方位45°、距離120 km、高度15 km，以速度830 km/h沿直線飛向雷達，勻速直線飛行40 s后開始作蛇形機動，轉彎向心加速度10g、轉彎半徑6.9 km，先向右轉彎90°，然后向左轉彎180°，再向右轉彎90°，最后沿著直線以速度830 km/h勻速直線飛行40 s結束。

圖2 目標飛行軌跡Fig.2 Target flight trajectory

跟蹤的殘差為預測值和量測值的差，殘差反映了模型預測的準確性和對目標的跟蹤效果，用來作為調整數據率的依據。對于相控陣雷達跟蹤和搜索(TAS)，采用了單波束進行跟蹤，為了提高測量精度，波束寬度設計得通常不寬。當跟蹤的殘差大于波束寬度的一半后就會造成失跟，雖然可以用更寬的波束進行重新捕獲，但是捕獲花費的時間資源是普通跟蹤的數十倍，失去了資源優化調度的意義。圖3和圖4為兩種算法進行50次仿真的分布。從圖3和圖4中可以看出：在0～40 s和100～140 s階段是勻速直線運動階段，兩種算法的殘差分布相當；在50～100 s階段蛇形機動階段，本文算法的跟蹤殘差要小。

圖3 Cohen算法的跟蹤殘差Fig.3 Tracking residual of Cohen algorithm

圖4 本文算法的跟蹤殘差Fig.4 Tracking residual of the proposed algorithm

為了更明顯地對比兩種算法的性能，計算了殘差的方差，結果如圖5所示。從圖5中可以看出本文算法優于傳統算法，本文算法的跟蹤穩定性高。

圖5 跟蹤殘差對比Fig.5 Comparison between tracking residuals of traditional algorithm and the proposed algorithm

由于目標主要在方位上進行機動，比較了方位上的跟蹤精度，跟蹤的誤差分布如圖6和圖7所示，可見誤差分布和殘差的分布相似，如果采用常系數阿爾法- 貝塔濾波器，則誤差和殘差有固定關系，誤差是殘差的一半。從圖6和圖7中也可以看出這種趨勢。

圖6 Cohen算法的跟蹤誤差Fig.6 Tracking error of Cohen algorithm

圖7 本文算法的跟蹤誤差Fig.7 Tracking error of the proposed algorithm

方位跟蹤精度如圖8所示。由圖8可見，本文算法的精度稍優，傳統算法隨著蛇形機動的方向改變，精度發生了起伏，主要原因是反向機動時殘差變小，跟蹤的數據率增大，但此時機動加速度大小沒有變化，采樣周期的變大造成了跟蹤精度的下降。

圖8 方位RMSE對比Fig.8 Comparison between bearing RMSEs of traditional algorithm and the proposed algorithm

圖9和圖10為兩種算法采樣周期的分布圖。由圖9和圖10可以看出：在反向機動的過程中由于沒有采用平滑因子，50 s、70 s和100 s附近是殘差最大的時刻，采樣周期最小，但是采樣周期變大；采用本文算法在目標的蛇形機動過程中，采樣周期仍采用較低的值；本文算法全程平均采樣周期為1.28 s，傳統算法平均采樣周期為1.79 s，精度和穩定性增加的代價是時間資源增多。

圖9 Cohen算法的跟蹤采樣周期Fig.9 Sampling period of Cohen algorithm

圖10 本文算法的跟蹤采樣周期Fig.10 Sampling period of the proposed algorithm

4.2 實驗2

目標的初始位置為：x軸25 km，y軸50 km，指向為90°，初始速度為0 m/s，加速度為0g. 雷達坐標在(0 km，0 km)位置。目標的運動情況為：經過30 s，以10 m/s2縱向加速度達到300 m/s的速度；從40 s開始持續40 s，以-20 m/s2的橫向加速度進行轉彎；從200 s開始持續15 s，以20 m/s2的橫向加速度進行轉彎，調整姿態，飛向雷達；從300 s開始持續30 s，以30 m/s2的橫向加速度進行轉180°，轉彎半徑約為2 km；從400 s開始持續10 s，以-90 m/s2的橫向加速度進行180°轉彎，轉彎半徑約為1 km；從450 s開始持續30 s，以10 m/s2的縱向加速度加速到600 m/s；從500 s開始持續22.5 s，以-90 m/s2的橫向加速度進行180°轉彎，轉彎半徑約為2 km；從620 s開始持續60 s，以-10 m/s2的縱向加速度，減速停止，目標運動軌跡如圖11所示。

圖11 目標飛行軌跡Fig.11 Targer flight trajectory

實驗2相對于實驗1機動更豐富，更能反映雷達的跟蹤情況。圖12所示為跟蹤的殘差對比。由圖12可見，在弱機動階段兩種算法的殘差水平基本相當，在突然機動階段由于對目標的運動描述模型與目標的實際模型相差很大，因此按照理論模型對目標的預測位置與目標的實際位置相差很大，殘差很大，如圖12中的223 s處和600 s處。殘差超過波束寬度的一半就容易造成目標失跟，圖12中傳統算法的殘差方位最大值達到1.8°，本文算法最大值為1.2°.

圖12 殘差RMSE對比Fig.12 Comparison between residual RMSEs of traditional algorithm and the proposed algorithm

圖13為兩種算法方位的跟蹤誤差對比，可見兩種算法跟蹤誤差曲線與殘差一致，只是幅度不同，在機動時誤差大，在弱機動時誤差小。圖14為兩種算法平均數據率，可見本文算法采用平滑因子，在采樣周期減少階段采用較大的平滑因子，因此采用周期不會突然增大，保證了目標的穩定跟蹤。

圖13 方位RMSE對比Fig.13 Comparison between bearing RMSEs of traditional algorithm and the proposed algorithm

圖14 跟蹤數據率對比Fig.14 Comparison between tracking data rates of traditional algorithm and the proposed algorithm

5 結論

本文針對警戒相控陣雷達目標跟蹤的時間資源優化分配問題，基于貝葉斯濾波算法提出一種基于殘差的跟蹤時間資源自適應分配算法，在直角坐標下進行濾波在球坐標系下計算殘差，并利用量測精度對殘差進行歸一化，引進平滑因子并在采樣間隔上升和下降時期采用不同取值。得到主要結論如下：

1)本文算法比傳統算法計算歸一化殘差簡單，避免了坐標轉換的非線性問題。

2)平滑因子能夠有效抑制目標在蛇形機動反向機動階段因殘差減小導致的采樣間隔擴大引起的失跟問題。