一種高分辨雷達空間目標跟蹤方法研究

劉政瑋, 陳 映, 魯耀兵

(北京無線電測量研究所, 北京 100854)

0 引 言

高分辨雷達通過發射寬帶信號獲取目標的高分辨一維距離像,從而實現目標的分類與識別[1-2]。由于其高距離分辨率,高分辨雷達在彈道導彈防御和空間目標監視等領域受到了越來越多的關注[3-4]。但是由于寬帶檢測、跟蹤等關鍵問題尚未完全解決,現役的高分辨雷達只能采用窄帶-寬帶信號交替發射的工作方式,即主要通過窄帶信號實現目標跟蹤,而寬帶信號則用于目標特征提取[5]。這種工作方式不僅增加了雷達系統設計的復雜性,還會造成雷達資源的浪費。因此,實現高分辨雷達目標跟蹤對空間目標跟蹤有著極其重要的現實意義。

空間目標運動在笛卡爾坐標系中更易于描述,在該坐標系描述的目標運動特征也更符合目標的本質運動規律,因此目標的運動方程通常建立在笛卡爾坐標系中,而雷達通常將目標量測描述為球坐標系下的距離、方位、俯仰[6]。由于目標狀態方程和雷達方程通常建立在不同的坐標系下,目標的跟蹤方法基本可以劃分為3類[7]:混合坐標系的目標跟蹤方法、笛卡爾坐標系下的目標跟蹤方法以及球坐標系下的目標跟蹤方法?；旌献鴺讼迪碌哪繕烁櫡椒ㄊ悄壳氨粡V泛應用的濾波方法。大多數用于空間目標跟蹤的非線性濾波技術,如擴展Kalman濾波器(extended Kalman filter, EKF)[8-10]、無跡Kalman濾波器(unscented Kalman filter, UKF)[11-12]、容積Kalman濾波器(cubature Kalman filter, CKF)[13]等均屬于混合坐標系下的目標跟蹤方法。這類方法的缺點在于嚴重依賴坐標轉換及非線性函數的線性化方式。笛卡爾坐標系的目標跟蹤方法將球坐標下的測量轉換到笛卡爾坐標系,使得轉換后的量測及觀測方程表現為一種偽線性形式[14],在笛卡爾坐標系完成目標狀態預測和狀態更新。由于轉化后的量測是有偏差的,會使得跟蹤精度下降,因此必須使用各種補償方法[15]。此時狀態更新中非線性處理體現于偽量測對應的誤差協方差矩陣的計算。球坐標系下的目標跟蹤方法則是將目標的運動方程和觀測方程在球坐標系下進行建模。顯然,這類建模的優點在于測量模型是線性、非耦合的;量測噪聲符合高斯分布假設;量測中的多普勒信息可被直接利用[16]等。缺點是球坐標系下的目標運動方程推導困難,推導后的運動方程是高度非線性的,各狀態間也是高度耦合的。復雜運動的彈道目標、空間目標,球坐標系下的目標運動學方程的非線性程度遠遠高于笛卡爾坐標系下的非線性程度,以至于量測方程的非線性優勢被抵消,導致濾波精度下降[17]。

對于同一空間目標,高分辨雷達帶寬越大、波束寬度越寬,即雷達的距離分辨率越高、角度分辨率越低,量測的非線性就越高。隨著量測非線性程度的增大,會出現一類特殊的非線性問題[18-19]:量測空間內呈現高斯狀的噪聲分布在更新過程中被轉換到狀態空間內時,其不確定性區域會變為強烈的非高斯分布,這種現象在高分辨雷達跟蹤彈道目標、空間目標時經常出現。當高分辨雷達在跟蹤空間目標時出現上述特殊的非線性問題時,EKF、UKF等濾波方法會出現跟蹤精度下降,甚至濾波器發散的問題。因此,本文首先針對這一類特殊非線性問題進行了分析,并提出一種基于中間過渡狀態的Kalman濾波方法。改進的Kalman濾波方法將量測空間內的目標位置量(距離、方位、俯仰)及其速度量(徑向速度、方位角速率、俯仰角速率)作為中間過渡狀態,則Kalman濾波更新過程被轉換到量測空間中進行,實現了量測方程的線性化。仿真實驗結果驗證了改進濾波方法可有效提高高分辨雷達跟蹤空間目標的精度。

1 空間目標跟蹤方法

在空間目標跟蹤中,通常使用狀態量的狀態方程描述目標運動特征,用量測方程描述雷達探測目標的位置信息,兩者通常建立在不同的坐標系中。同時,由于空間目標的動力學特性復雜,因此雷達對空間目標的跟蹤涉及坐標系統、運動建模和濾波算法等方面的內容。

1.1 雷達坐標系

雷達探測到的位置信息是對空間目標進行跟蹤的基礎,現有的單脈沖雷達多提供在雷達站球坐標系下描述的徑向距離、方位和俯仰信息。在雷達站球坐標系下,空間目標的運動方程比較復雜,為簡化跟蹤過程,通常在直角坐標系下描述運動方程。根據實際應用情況,本文在雷達北天東直角坐標系中分析空間目標運動情況,在雷達球坐標系中探測目標徑向距離、角度信息,兩種坐標系的相互關系如圖1所示。

圖1 雷達站北天東坐標系與雷達站球坐標系Fig.1 Radar station north sky east coordinate system and radar station spherical coordinate system

雷達站北天東坐標系是一種站心坐標系,是以雷達中心為坐標原點Or,OrXr和OrZr軸是地球參考橢球的切線,分別指向北和東,而OrYr軸垂直于當地水平面向上。

雷達站球坐標系原點為雷達站,R為斜距,方位角A為正北軸OrXr順時針旋轉到目標與雷達站在水平面上的投影角度,俯仰角E定義為目標與雷達站連線與水平面的夾角。

1.2 空間目標運動模型及觀測模型

目前,空間目標跟蹤常用的3種地球模型分別是平地球模型、球地球模型和橢球地球模型[20],在實際應用中,可以根據空間目標的飛行距離相對于地球半徑的大小選擇合適的模型。對于長時間的空間運動目標跟蹤需要考慮由地球形狀的不規則引起的攝動,因此采用橢球地球模型可以獲得更準確的重力加速度模型ag:

(1)

在雷達北天東直角坐標系中分析空間目標運動受力情況更加直接,過程噪聲也更易于描述,因此大多數文獻在直角坐標系中建立目標的運動方程。在雷達北天東坐標系中,設空間目標的狀態變量為x=[x,vx,y,vy,z,vz],則空間目標運動狀態微分方程可描述為[21]

(2)

(3)

式中:ae、ac為雷達北天東坐標系中由地球自轉引起的牽連加速度矢量和柯氏加速度矢量,ad為氣動阻力加速度矢量。3項加速度矢量的具體表達式可參考文獻[20]。

通常在雷達測量坐標系下獲得目標量測Z=[R,A,E],因此目標的量測方程為

(4)

式中:vR、vA、vE為相互獨立的零均值高斯分布。

1.3 傳統的空間目標濾波方法

目前空間目標跟蹤方法中應用最為廣泛的是混合坐標系下的目標跟蹤方法,其因為直觀的運動方程及量測方程建模方法,受到眾多學者的關注[22-23]。其中,EKF通過泰勒級數展開的方法,獲得非線性系統的線性近似,進而采用卡爾曼濾波處理非線性系統的濾波問題。UKF、CKF則是通過確定性采樣點的變換求出后驗概率密度函數的均值和協方差,避免了對非線性函數的近似。

另外一種將非線性過程線性化的方式是將目標的量測值轉換到直角坐標系,在直角坐標系中進行目標狀態的預測更新,這種方法稱為轉換量測卡爾曼濾波器(converted measurement Kalman filter, CMKF)。傳統的CMKF存在偏差,可以通過將其乘以一個偏移因子實現去偏,但該方法存在兼容性問題[24],文獻[25]提出的改進方法對其計算結果進行了修正。這些算法在計算轉換誤差統計特性時均以量測值為條件,會造成增益矩陣與量測誤差相關,最終導致估計結果有偏。因此,有學者提出以目標預測信息為條件計算轉換誤差特性的無偏CMKF,從而克服濾波參數與量測誤差之間的相關性[26-27]。但是,無論是以量測值為條件還是以預測值為條件的轉換方法都會造成轉換后的量測噪聲是非高斯的,且在各坐標方向上互耦。

球坐標系下的目標跟蹤方法在量測空間內選取目標運動狀態,同樣可以實現量測方程的線性化。這類方法需要直接在球坐標系下描述目標運動方程,但是對于空間運動目標,其在運動過程中受到各種力的作用,在球坐標系中難以直接進行建模。在實際應用中,當目標運動的假設條件發生變化時,需要重新對其進行建模分析。另一方面,在目標的運動方程轉換到球坐標系時,通常不考慮過程噪聲的轉換,這時極易出現模型失配。

1.4 高分辨雷達跟蹤空間目標的非線性問題

高分辨雷達在跟蹤空間目標時,受到雷達帶寬、波束寬度、信噪比等因素影響,高分辨雷達的測距精度高,測角精度較低。當笛卡爾坐標系中的狀態使用來自球坐標系的非線性量測更新時,量測的不確定性區域會呈現出非常明顯的非高斯分布。這種分布呈現出一種曲線形狀(在二維空間內類似于香蕉形狀,在三維空間內類似于隱形眼鏡形狀),因此將這一分布稱為“隱形眼鏡(contact lens, CL)”分布,并將這一類非線性問題簡稱為CL問題[28]。本文以極坐標系到二維直角坐標系轉換為例,研究狀態空間內的量測分布特性及高分辨雷達量測非線性的產生機理,雷達站球坐標系到雷達北天東坐標系的轉換可由此類推。

令(x0,y0)為雷達直角坐標系下目標的真實位置坐標,(r0,θ0)為雷達極坐標系下目標的真實位置,則目標在雷達極坐標中的距離和方位量測值為

(5)

(6)

將其在目標位置真值(x0,y0)處進行泰勒展開得

(7)

由式(7)可以看出,在雷達直角坐標系下,目標的量測值不再符合高斯分布,且其真實分布復雜,難以獲得。同時,式(7)表明直角坐標系下目標量測分布與目標的真實位置及測距和測角噪聲有關。

圖2從幾何角度描述了這種量測的不確定性區域的形成原因。對于遠程空間目標,在雷達測距精度不變的情況下,若雷達探測角度σθ誤差增加,則跨角度誤差范圍rσθ相應增大,由于其遠遠大于測距誤差σr,最終導致量測的不確定性區域逐漸呈現為曲線狀。

圖2 量測非線性CL問題Fig.2 CL problem of measurement nonlinearity

Lerro等[29]提出一種非線性度量指標來描述量測在坐標轉換中偏離高斯分布的程度:

(8)

假設高分辨雷達對空間目標的作用范圍為2 000 km,圖3描述了不同β值對量測非線性的影響程度。圖中藍色點跡為104個蒙特卡羅仿真點跡,綠色實線為根據直角坐標系下量測值泰勒展開得到的3σ誤差橢圓理論值,紅色實線為根據蒙特卡羅仿真點跡得到的3σ誤差橢圓真實值。由圖3(a)可以看出,β值越大,高分辨雷達對空間目標量測的距離噪聲誤差越小。測角噪聲增大時,由量測高斯假設得到的3σ誤差橢圓并不能體現直角坐標系下量測分布的真實特性,會導致經典非線性濾波器的跟蹤性能下降。從圖3(b)可以看出,3σ誤差橢圓可以很好地描述量測的真實分布,各非線性濾波器仍適用,但此時雷達帶寬不符合高分辨雷達帶寬。

圖3 不同β值對量測非線性的影響Fig.3 Influence of different β value’s on measurement nonlinearity

當空間目標跟蹤中出現強非線性CL問題時,傳統的非線性濾波器的跟蹤精度會急劇下降,甚至出現濾波發散。在圖4中,黑色區域表示量測不確定性區域和預測狀態不確定性區域的交集區。在濾波的初始階段,目標狀態估計不夠準確,此時交集區域明顯是彎曲且非高斯的,如圖4(a)所示?；诟咚辜僭O的EKF濾波器會出現一致性問題,EKF濾波器極有可能出現發散。當目標狀態估計準確時,交集區可以近似地視作符合高斯分布,如圖4(b)所示,此時EKF可以較好地跟蹤目標,但由于其在線性化時只保留一階偏導項,跟蹤精度較差。在這一過程中,由于UKF、CKF通過采樣點逼近量測不確定區域,因此其可以較好地實現目標跟蹤。與EKF、UKF、CKF不同,CMKF是通過增加量測不確定性區域的“厚度”,即通過犧牲一定的測距精度來保證濾波器的一致性,使得交集區可以近似為高斯分布,如圖4(c)所示,因此CMKF可以應對CL問題,但其測距誤差會急劇增大[30]。

圖4 濾波器性能幾何圖解Fig.4 A geometric illustration diagram of filter performance

2 改進的Kalman濾波器

2.1 矢量分析法求解不同坐標系下的目標狀態轉換

(9)

球坐標系中,目標的距離、方位、俯仰為

(10)

令ex、ey、ez分別為直角坐標系中OrXr、OrYr、OrZr軸的單位矢量,在球坐標系中取eR為沿徑向方向的單位矢量,eA為垂直于徑向方向并指向方位角A增大一方的單位矢量,eE定義與eA同理。因此,目標徑向矢量可表示為

R=ReR=RcosEcosAex+RsinEey+RcosEsinAez

(11)

由此可得,球坐標系中各單位向量為

(12)

所以

(13)

所以,球坐標系中速度量為

(14)

將式(12)代入得

(15)

綜上所述,直角坐標系中的目標速度為

(16)

球坐標系中,目標徑向速度、方位角速率、俯仰角速率分別為

(17)

2.2 改進的Kalman濾波器

高分辨雷達在跟蹤空間目標時易出現非線性問題——CL問題,主要原因在于球坐標系中的量測轉換到直角坐標系時,其不確定性區域呈現出非常強烈的非高斯分布。因此,改進的Kalman濾波器將濾波器的更新過程轉換到量測空間內進行,選取量測量(距離、方位、俯仰)及其一階導數(徑向速度、方位角速率、俯仰角速率)作為目標在量測空間內的狀態,對目標進行更新,這樣在量測分布符合高斯假設的同時,量測方程也相應地變為線性方程。

為方便后續方法的比較,假設空間目標的非線性系統狀態估計的一般描述為

(18)

則改進的Kalman濾波器計算步驟如下:

步驟 1預測。假設預測濾波器輸入為k時刻濾波結果xk=[xk,vxk,yk,vyk,zk,vzk]和Pk,通過下式獲得狀態預測結果xk+1|k和Pk+1|k。

xk+1|k=f(xk)

(19)

(20)

(21)

式中:Qk為過程噪聲wk的協方差。

容積準則能夠將笛卡爾坐標系下的高斯分布的運動狀態轉換為某個多維幾何體容積的計算,具有較高的計算效率,能夠獲得較高的數值精度。因此,改進的Kalman濾波器根據三階球面徑向準則將狀態空間內的預測協方差Pk+1|k轉換為量測空間內的協方差PMk+1|k。

通過下式對xk+1|k采樣,得到2n個采樣點:

(22)

式中:ei表示第i個元素為1的單位向量。

采樣點經過式(10)和式(17)轉換后的值為Ci(i=1,2,…,2n)。

則量測空間內中間狀態的協方差PMk+1|k為

(23)

(24)

步驟 3更新。由狀態xMk+1|k及其協方差PMk+1|k得到量測空間內的狀態估計結果xMk+1PMk+1。

由于更新過程在量測空間進行,因此原本的非線性量測矩陣h(xk)變為線性矩陣Hz,即

(25)

量測空間內目標運動狀態更新過程為

(26)

(27)

xMk+1=xMk+1|k+KMk+1(zk+1-HzxMk+1|k)

(28)

PMk+1=(I6-KMk+1Hz)PMk+1|k

(29)

步驟 4狀態估計結果提取。根據式(9)、式(16)及三階球面徑向準則將量測空間內的狀態估計結果xMk+1和PMk+1轉換得到k時刻直角坐標系中的目標濾波估計結果xk+1和Pk+1,這一步可認為是中間過渡狀態轉換的逆過程。

改進的Kalman濾波器與CKF均采用了容積變換準則,但是二者有本質的不同。由于CKF在混合坐標系中進行目標的預測更新,失去了Kalman濾波器本身具有的一致性和無偏性。改進的Kalman濾波器則是在預測后將狀態空間內的目標預測值通過容積準則映射為量測空間內的中間過渡狀態,從而得到線性的Kalman更新過程,更新后再將中間過渡狀態逆變換到直角坐標系內,得到目標狀態的一致性估計。同時在轉換過程中,不同于CKF僅僅將空間目標的六維運動狀態轉換為量測空間內的三維位置預測信息,改進的Kalman濾波器選取了量測空間內目標的六維運動信息,進一步保證了轉換過程的一致性。另一方面,相比CKF需要對復雜空間目標的運動方程進行采樣,改進的Kalman濾波器能夠在保證濾波一致性的同時,減少通過容積變換求取目標預測值及協方差的步驟,目標的狀態估計是通過線性Kalman濾波器進行的,因此計算量會減少。

3 仿真實驗

3.1 仿真場景

仿真的彈道目標射程為3 700 km,在雷達北天東坐標系下,仿真生成的彈道軌跡如圖5所示,其中,“°”和“△”分別表示彈道目標發射點位置和落點位置。

圖5 雷達北天東坐標系彈道目標軌跡圖Fig.5 Radar north sky east coordinate system trajectory diagram of ballistic target

雷達的探照范圍為2 000 km,數據率為10 Hz,雷達跟蹤彈道目標時長為100 s。在雷達球坐標系下,導彈的徑向距離、方位角和俯仰角隨時間的變化關系分別如圖6、圖7、圖8所示,其中藍色段為彈道目標軌跡,紅色段為雷達探測弧段。

圖6 目標距離隨時間的變化關系Fig.6 Change of target distance with time

圖7 目標方位角隨時間的變化關系Fig.7 Change of target azimuth with time

圖8 目標俯仰角隨時間的變化關系Fig.8 Change of target pitch angle with time

3.2 高分辨雷達跟蹤彈道目標濾波仿真

基于第3.1節仿真給定的彈道目標運動軌跡,本文采用典型高分辨雷達,測距誤差為0.1 m。針對不同的測角誤差,本文通過50次蒙特卡羅求均方根誤差(root mean square error, RMSE)平均值的方法,進行UKF、CKF、EKF、兩種無偏CMKF分別為無偏CMKF(unbiased CMKF, UCMKF[31])、基于預測信息的CMKF(prediction position CMKF, PRE_CMKF),以及本文改進Kalman濾波器的性能對比。濾波器均采用兩點初始化方法,其他仿真參數設置如表1所示。

表1 跟蹤場景

在仿真場景1中,雷達跟蹤彈道目標時的測距誤差較大,測角精度高,量測的非線性程度低。在這種情形下,通過圖9可以看出,UKF、CKF、EKF及本文改進的Kalman濾波器的性能相近,而UCMKF的測距誤差較大,PRE_CMKF相較于UCMKF濾波性能有所提升。

圖9 仿真場景1濾波結果Fig.9 Simulation scenario 1: filtering results

隨著雷達性能提升,高分辨雷達的測距精度提高,圖10為仿真場景2時目標的狀態估計結果。從圖10中可以看出,高分辨雷達測角精度較高時,由量測轉換引起的非線性程度較低,本文提出的改進Kalman濾波器和UKF、CKF測距、測角精度相近。而EKF濾波誤差急劇增大,這是由于EKF在完成系統模型的線性化近似時忽略了高階項,所以導致EKF估計的狀態后驗分布產生較大誤差。而此時兩種CMKF雖然可以保證目標的穩定跟蹤,但它們的測距誤差較大。

圖10 仿真場景2濾波結果Fig.10 Simulation scenario 2: filtering results

由于當前高分辨雷達測角技術發展的限制,高分辨雷達跟蹤彈道目標時測量角度誤差較大,因此仿真場景3中假設方位角誤差為σA=0.5°。在高分辨雷達中,隨著測角誤差的增大,由量測值轉換引起的CL問題的非線性程度增大。從圖11中可以看出,EKF會出現濾波發散的情形。此時,兩類CMKF濾波器通過犧牲一定程度的測距精度,保證了濾波器的一致性,其仍然能夠跟蹤目標。本文提出的Kalman濾波器測距精度和UKF、CKF基本一致,但由于測角誤差增大,UKF、CKF的方位角精度急劇下降,在彈道目標防御過程中將會導致一系列嚴重的后果,而本文提出的改進Kalman濾波器的角度RMSE仍可收斂至較小的測角誤差。因此,當雷達球坐標系量測轉換到直角坐標系中出現CL分布時,本文提出的改進Kalman濾波器比UKF、CKF具有更優的濾波性能,彈道目標跟蹤精度更高。

圖11 仿真場景3濾波結果Fig.11 Simulation scenario 3: filtering results

為了比較各濾波器的計算復雜度,圖12列出了不同濾波器單次狀態估計的平均時間損耗。從圖中可以看出,基于中間過渡狀態的Kalman濾波器在保證濾波精度的同時,可以避免對復雜的運動方程進行采樣,計算復雜度比UKF、CKF更小。

圖12 平均時間消耗Fig.12 Average time consumption

4 結 論

本文研究了高分辨雷達跟蹤空間目標時,量測值由于坐標轉換在狀態空間內不確定性區域呈現出嚴重非高斯狀態的問題。此時,傳統的非線性濾波器的狀態估計精度會急劇下降。本文提出一種基于中間過渡狀態的Kalman濾波器,在直角坐標系下進行目標狀態預測,避免了球坐標系下對目標運動的復雜建模;預測后選取量測空間內目標位置和速度量及其協方差作為中間過渡狀態,在量測坐標系下進行目標狀態線性更新;更新后再將中間過渡狀態映射到直角坐標系內,從而得到目標狀態的一致性估計。仿真結果表明,相對于現有的非線性濾波器,改進算法可有效提高高分辨雷達對空間目標的狀態估計精度,這將對雷達系統快速準確形成目標軌跡有著重要作用。