高超聲速機動目標天基跟蹤魯棒性濾波方法

魏世君，翟 光，孫一勇，畢幸子，汪宏昇

(1. 北京理工大學宇航學院，北京 100081；2. 中國科學院微小衛星創新研究院，上海 201203；3. 中國運載火箭技術研究院研究發展部，北京 100076)

0 引 言

高超聲速滑翔飛行器(Hypersonic glide vehicle, HGV)可在臨近空間內實現5倍以上聲速的機動突防飛行。為了實現對HGV的有效攔截，需首先完成對HGV飛行軌跡的準確跟蹤。受地球曲率、大氣散射等因素影響，地基雷達難以對HGV進行持續穩定的探測跟蹤。由于HGV飛行過程中紅外特性顯著，天基紅外傳感器在未來將成為HGV的有效探測手段。然而，HGV在高速飛行過程中不斷實施軌跡機動，由于機動幅值及起始時間等信息未知，采用傳統卡爾曼濾波等方法進行軌跡估計跟蹤時，極易導致濾波誤差發散并丟失目標。因此仍需基于天基紅外探測信息，設計針對目標未知機動的魯棒性濾波跟蹤算法，實現對目標飛行狀態的準確估計。

當前，針對HGV軌跡跟蹤問題，大部分研究工作采用的方法為強跟蹤、擴維卡爾曼濾波以及交互多模型濾波。文獻[10]采用多模型交互濾波，對Singer,CA,CT等運動學模型的跟蹤結果進行融合；文獻[11]在Singer模型和擴維卡爾曼濾波的基礎上，基于多個具有不同時間常數的子濾波器構建交互多模型濾波器，以適應目標的不同機動頻率；文獻[12-14]根據氣動特性和制導律推導了機動模型的參數取值，適用于特定的氣動外形以及制導律，但對于不同型號或采用不同制導律的目標需要重新進行推導；文獻[15]研究了與HGV類似的機動再入飛行器(MaRV)跟蹤問題，為彈道系數構建維納過程模型，通過動力學反解得到彈道系數的偽量測，從而對其進行最優估計，但該方法對動力學模型進行了高度簡化，難以應用于復雜模型。總之，在HGV的跟蹤問題中，通常采用維納過程、馬爾可夫過程等統計模型對機動加速度進行建模，并采用擴維卡爾曼濾波進行狀態估計，通過強跟蹤或交互多模型方法增強濾波器的魯棒性。然而，單一跟蹤模型往往存在選參問題；若采用交互多模型手段增強魯棒性，則又會產生較大的計算開銷，不利于實時跟蹤。此外，已有研究工作多基于傳統雷達測量，涉及天基觀測的研究工作尚為少見。

針對現有工作的不足，本文研究了天基紅外觀測背景下存在未知時變機動的HGV目標跟蹤問題，提出了一種基于機動觀測器補償的魯棒擴展卡爾曼濾波(Maneuver observer compensated robust extended Kalman filter, MOCREKF)。構造機動觀測器，對HGV目標氣動加速度進行實時估計，并推導了該機動觀測器的誤差界；利用機動觀測器的輸出信息修正擴展卡爾曼濾波的動力學預測步驟，解決模型失配的問題；為克服機動觀測器中低通環節帶來的時延誤差，在擴展卡爾曼濾波的更新步驟引入次優漸消因子，增強濾波器的魯棒性。通過HGV典型跟蹤場景進行仿真對比，驗證了MOCREKF相對于現有算法在跟蹤精度、魯棒性和實時性上的優勢。

1 HGV滑翔段動力學模型

(1)

式中：為引力加速度向量；為氣動加速度向量。

HGV滑翔飛行時間相對較短(約500～3000 s)，在引力加速度中僅考慮二體引力和攝動項，則有

(2)

在速度系下對氣動加速度進行解算，HGV的地球相對速度向量為

(3)

(4)

(5)

(6)

式中：為HGV質量；,分別為阻力系數和升力系數；為當地大氣密度；為地球相對速度的模長；為HGV參考面積；為HGV傾側角。記氣動加速度在J2000慣性系下的坐標分量為A,A,A，則有

(7)

式中：為速度系到慣性系的變換矩陣，根據HGV的位置向量和地球相對速度計算得到。至此，已得到式(1)中和的表達式。

氣動加速度向量即為HGV目標根據制導和控制指令生成的機動加速度，通過調整，可產生不同類型的滑翔彈道。對于非合作HGV目標，未知，在目標跟蹤過程中必須予以相應的處理。

2 天基紅外立體觀測模型

圖1 天基紅外觀測模型Fig.1 Model of space-based infrared measurement

2.1 量測方程

在不影響結論的前提下進行簡化，假定各衛星對視線角度,(=1,…,，為協同觀測的衛星數目，后文假定≥2)的測量信息被轉換至J2000慣性系下，則衛星的視線角度量測方程為

(8)

(9)

(10)

式(8)右側第一項對時間求導，可得到衛星的視線角度變化率量測方程

(11)

2.2 目標位置速度反解

為便于后續機動觀測器的設計，需根據測角信息反解目標的位置及速度。忽略式(8)右側的量測噪聲，不考慮衛星自身星歷誤差，綜合顆衛星的測量信息，可整理得到關于目標位置的方程組

·=

(12)

式中：系數矩陣及列矩陣的表達式為

(13)

(14)

方程組(12)中存在冗余測量信息，采用最小二乘法可求解得到受噪聲污染的目標位置

(15)

忽略量測方程(11)右側的噪聲項，同樣不考慮衛星自身星歷誤差，綜合顆衛星的視線角變化率測量信息，可整理得到關于目標速度的方程組

·=

(16)

式中：系數矩陣及列矩陣的表達式為

(17)

(18)

(19)

受量測噪聲的影響，從角度和角度變化率量測信息中求解出的目標位置和速度精度較差，僅用于構建機動觀測器。為了結合動力學模型對目標位置和速度進行最優估計，應當將式(8)和(11)給出的量測信息輸入至濾波器。

3 機動觀測器

目標實施機動時，若機動信息已知，則常規濾波器即可達到較好效果；若機動信息未知，跟蹤模型與目標實際運動模型的失配將導致濾波器跟蹤發散。針對目標的未知機動，本節設計一種機動觀測器，從量測中提取目標的機動估計信息，從而抑制濾波發散、提高跟蹤精度。

3.1 模型離散化

+1=()++

(20)

=()+

(21)

記為自然數集，為實數空間，以上兩式中：∈為時間步數；,,分別為時刻的狀態向量、量測向量及未知輸入向量；∈為輸入驅動矩陣；∈為噪聲驅動矩陣，取值與相同；∈為零均值高斯過程噪聲，反映建模及離散化誤差；∈為零均值高斯量測噪聲；()由動力學方程(1)離散化得到；()聯立了兩顆衛星如式(8)、(11)所示的量測方程。在此情況下，任意時刻的均趨近于常值=[05T]，其中為單位陣，為采樣周期。后文在必要時以常值代替，并省略單位陣的維度下標。后續機動觀測器和濾波器的設計及誤差分析均基于離散化的狀態空間模型(20)、(21)展開。

“放心吧，現在什么事情也不會有。”他說。“這房子是租下的，簽著租賃合同。在租賃合同期滿之前，從法律上講這房子是屬于咱們的。”

3.2 機動觀測器迭代步驟

為對目標的未知機動進行估計，構造機動觀測器

(22)

(23)

(24)

(25)

式中：∈為待定矩陣。將式(25)代入式(24)，合并和+1的同類項，整理得到

(26)

式中：

(27)

(28)

(29)

(30)

在遞推過程中，該機動觀測器僅依賴于當前量測。獲取后，等式右側的各個量均為已知，即可給出目標的未知輸入估計。此外，上述推導過程中不存在針對未知輸入的假設，因此該機動觀測器在使用中無需獲取目標機動模式的先驗信息，理論上可對目標的各類機動模式進行自適應估計。

3.3 機動觀測器誤差分析

本節將對機動觀測器的誤差傳播特性進行分析，確定增益矩陣的取值范圍，并給出機動觀測器的誤差界。將機動觀測器的誤差定義為

(31)

將式(29)代入上式，得到

(32)

(33)

結合式(20)和式(30)，上式可展開為

(34)

(35)

(36)

式中：

(37)

式(36)描述了機動觀測器的誤差傳播特性，易知項有界。故機動觀測器誤差收斂的必要條件為

|1-|<1 (=1,2,3)

(38)

上式給出了機動觀測器增益矩陣的取值范圍。

(39)

(40)

在式(38)的前提下，當→∞時，有

(41)

對于本文研究的問題，有以下條件成立：

(1) HGV在滑翔飛行過程中，其氣動加速度連續，未知輸入的變化率存在上界，即

(42)

(3) 存在正實數,，使得

(43)

依據條件(1)，對于式(41)右側第一項，有

(44)

因此

(45)

依據條件(2)和(3)，對于式(41)右側第二項，有

(46)

對于式(41)右側第三項，有

(47)

由于假定-1服從高斯分布，取其3上界，則

(48)

(49)

據此，由矩陣范數性質

(50)

可知，對于式(41)所示的機動觀測誤差，有

(51)

上式給出了機動觀測器的誤差界。

3.4 噪聲抑制

(52)

式中：為時間常數，決定了該環節的通頻帶寬=1。當過小時，慣性環節對噪聲分量的抑制作用不明顯；當過大時，該環節會產生較強的時滯，有可能導致濾波發散。受結構強度、熱防護等條件約束，HGV的機動頻率無法產生較大變化。在實際應用中，可根據HGV的常規機動頻率選取。

4 基于機動觀測器的魯棒擴展卡爾曼濾波

將機動觀測器的輸出作為未知機動估計補償至擴展卡爾曼濾波的預測步驟中，即可對式(20)、(21)所示的非線性系統進行狀態估計。當目標進行快速時變機動時，低通環節(52)將會產生不可忽略的時延誤差。為克服時延誤差、進一步提高濾波算法的魯棒性，在擴展卡爾曼濾波的更新步驟中引入次優漸消因子。MOCREKF主要分為初始化、狀態預測、漸消因子計算、狀態更新和觀測器更新五個步驟。

(1)初始化

(2)狀態預測

計算狀態預測

(53)

(3)漸消因子計算

計算漸消因子

(54)

(55)

(56)

(57)

(58)

(4)狀態更新

利用漸消因子對先驗方差進行調整，有

(59)

(60)

(61)

(62)

式中：為MOCREKF的濾波增益矩陣。

(5)觀測器更新

(63)

(64)

式(53)～(64)構成了MOCREKF的迭代流程，算法結構如圖2所示。初始化完畢后，重復執行步驟(2)～(5)，即可實現對HGV目標的軌跡跟蹤。不同于擴維卡爾曼濾波、交互多模型濾波等方法，MOCREKF無需使用目標機動的先驗信息，理論上對各類機動模式均具備良好的自適應能力，漸消因子的引入也使濾波器具備了較強的魯棒性。此外，MOCREKF僅在單個擴展卡爾曼濾波器的基礎上添加了觀測器更新步驟和漸消因子計算步驟，與多模型濾波方法相比具有更低的計算負擔。

圖2 MOCREKF算法流程Fig.2 Algorithm flowchart of MOCREKF

5 仿真校驗

現構建HGV目標的天基觀測仿真場景，對本文提出的MOCREKF算法進行性能驗證。根據HTV-2相關參數設計跳躍滑翔彈道和擬平衡滑翔彈道，如圖3所示，代表兩類典型機動模式。目標在(121°E, 52°N)上空約60 km高度處進入滑翔段，初始速度設置為6.5 km/s，在30～60 km高度范圍內向西南方向飛行。構建搭載高精度凝視紅外傳感器的低軌星座，采用Walker構型，衛星總數為200顆，均勻分布在20個500 km高度軌道上。由于單個衛星觀測時長受限，仿真過程中采用星間接力的方式進行持續跟蹤。協同觀測的衛星數目為=2，始終選取對目標可見且基線最佳的觀測組合，忽略目標捕獲過程中的誤差以及衛星自身的定軌誤差。依據現階段傳感器精度水平，單星方位角和俯仰角測量標準差選取為5×10rad，角度變化率測量標準差選取為5×10rad/s。濾波跟蹤總時長為1400 s，離散周期為=1 s，濾波器過程噪聲標準差置為0.01 m/s。

圖3 HGV滑翔段彈道Fig.3 Trajectory of HGV in glide phase

作為對比，采用強跟蹤濾波(Strong tracking filter, STF)、擴維卡爾曼濾波(Augmented state Kalman filter, ASKF)、交互多模型(Interactive multiple model, IMM)濾波以及本文提出的MOCREKF分別對HGV目標進行跟蹤。文獻[9]對標準卡爾曼濾波跟蹤機動目標時出現的發散現象進行了詳細的仿真和分析，故此處不再采用標準卡爾曼濾波進行對比。各跟蹤算法中，ASKF的加速度分量采用一階馬爾可夫模型描述，在不使用先驗信息的情況下將機動頻率參數設定為0.3和1，分別記為ASKF1和ASKF2。在交互多模型濾波中，設計五個擴維子模型，機動頻率參數在區間[0,1]內等間隔選取。在MOCREKF中，根據式(38)將觀測器增益矩陣選取為=diag(1, 1, 1)，低通環節時間常數選取為=20。各個濾波器的位置、速度初值根據初始時刻的測量信息求解得到，機動估計初值置為零。

5.1 場景一：擬平衡滑翔彈道跟蹤

在擬平衡滑翔過程中，目標機動加速度變化較為平緩。分別采用STF, ASKF1, ASKF2, IMM, MOCREKF對圖3中的擬平衡滑翔彈道進行跟蹤，在i7-1185G7, 32G RAM, MATLAB環境下執行1000次蒙特卡洛仿真，得到位置、速度、加速度估計的均方根誤差(RMSE)曲線，分別如圖4～圖6所示。繪制單次仿真中的觀測器輸出曲線，如圖7所示。統計各算法平均單步RMSE及平均單步運行耗時，見表1。

圖4 位置RMSE(場景一)Fig.4 Position RMSE (Scenario 1)

圖5 速度RMSE(場景一)Fig.5 Velocity RMSE (Scenario 1)

圖6 機動加速度RMSE(場景一)Fig.6 Maneuver acceleration RMSE (Scenario 1)

圖7 觀測器機動加速度估計(場景一)Fig.7 Maneuver acceleration estimation of the observer (Scenario 1)

表1 各濾波器平均RMSE及耗時(場景一)Table 1 Average RMSE and computing time of filters (Scenario 1)

觀察圖4、圖5可知，在不使用目標機動先驗信息時，難以對擴維卡爾曼濾波中的機動模型參數(機動頻率)進行合理取值，致使單模型ASKF1, ASKF2跟蹤發散。在對目標機動規律進行充分分析的前提下，方能在馬爾可夫過程、維納過程等機動模型中給出合理的參數取值，但在多數情況下目標的機動信息難以事先獲取。STF根據漸消因子對方差進行修正，強迫濾波器收斂；IMM濾波對多個機動模型的輸出結果進行融合，利用模型集對目標可能的機動模式進行覆蓋；MOCREKF在使用機動觀測器修正跟蹤模型的同時，進一步利用漸消因子增強濾波器的魯棒性。在圖4和圖5中，STF, IMM, MOCREKF均取得了良好的效果，收斂到了相似的誤差下界。但STF的速度估計誤差略高于IMM和MOCREKF。

觀察圖6可知，ASKF1, ASKF2在機動加速度方面也無法給出準確的估計。STF僅通過方差膨脹迫使濾波器收斂，不具備機動加速度估計的能力。盡管IMM的加速度估計誤差顯著低于ASKF1和ASKF2，但受模型集中失配模型的影響，其加速度的估計誤差仍然較高，且存在輕微的發散趨勢。與IMM相比，MOCREKF能夠給出精度更高、收斂性更為良好的加速度估計。觀察圖7可知，本文所設計的機動觀測器能夠較好地逼近機動加速度真值。

由表1可知，本文提出的MOCREKF在精度上優于ASKF和STF，與IMM相近。然而MOCREKF在計算耗時上遠低于IMM，這是因為MOCREKF未進行擴維處理，僅在標準EKF的基礎上添加了觀測器和漸消因子遞推步驟，其運行耗時與單模型擴維卡爾曼濾波器相近。IMM中包含五個擴維子模型以及模型交互步驟，故其運行耗時約為單模型擴維卡爾曼濾波器的五倍。

5.2 場景二：跳躍滑翔彈道跟蹤

與擬平衡滑翔相比，HGV在跳躍滑翔過程中的機動加速度變化較為頻繁。各濾波器配置同場景一，采用STF, ASKF1, ASKF2, IMM及MOCREKF對圖3中的跳躍滑翔彈道進行跟蹤，執行1000次蒙特卡洛仿真，得到位置、速度、加速度的RMSE曲線，分別如圖8～圖10所示。繪制單次仿真中的觀測器輸出曲線，如圖11所示。統計各算法平均單步RMSE，見表2。

觀察圖8、圖9可知，單模型濾波器ASKF1和ASKF2均由于模型失配而產生發散，無法給出準確的狀態及機動加速度估計。STF, IMM及本文提出的MOCREKF在位置和速度估計上具備較好的收斂性。觀察圖10可知，MOCREKF具有比IMM更高的機動加速度估計精度，這與場景一的結論具有一致性。觀察圖11可知，MOCREKF中的機動觀測器在跳躍滑翔機動模式下同樣能夠較好地逼近目標加速度真值。由表2數據可知，在場景二中MOCREKF的精度仍然高于STF，且與IMM保持在相同水平。

圖8 位置RMSE(場景二)Fig.8 Position RMSE (Scenario 2)

圖9 速度RMSE(場景二)Fig.9 Velocity RMSE (Scenario 2)

圖10 機動加速度RMSE(場景二)Fig.10 Maneuver acceleration RMSE (Scenario 2)

圖11 觀測器機動加速度估計(場景二)Fig.11 Maneuver acceleration estimation of the observer (Scenario 2)

表2 各濾波器平均RMSE及耗時(場景二)Table 2 Average RMSE and computing time of filters (Scenario 2)

前文機動觀測器設計過程中未使用HGV目標機動的任何先驗信息，其機動估計能夠自適應地逼近目標機動真值，故MOCREKF對HGV目標的機動模式不敏感。在場景一和場景二所代表的兩類典型機動模式下，MOCREKF無需調整濾波器結構及參數，均能夠以較低的計算開銷取得收斂性良好、精度較高的HGV目標狀態及機動加速度估計。

6 結 論

為解決存在未知時變機動的HGV目標軌跡跟蹤問題，設計機動觀測器估計目標的未知機動，并提出了一種MOCREKF。擬平衡滑翔和跳躍滑翔兩類典型機動場景的仿真分析表明：

(1) 本文所提MOCREKF中的機動觀測器可根據測量信息有效提取出目標的機動估計，與強跟蹤濾波、擴維卡爾曼濾波相比，MOCREKF在HGV目標跟蹤問題中具有更優的穩定性和精度；

(2) MOCREKF設計過程中未使用目標機動規律的先驗信息，在場景一和場景二中均表現出良好的性能，對目標的機動模式具有較強的魯棒性；

(3) MOCREKF的實時性良好，其計算負擔大致與單模型擴維卡爾曼濾波相同，遠低于相同精度水平的交互多模型濾波。