基于不完全測量數據的飛行器軌跡參數估計*

李 冬，劉 學

(中國人民解放軍91550部隊， 遼寧 大連 116023)

在航天測控領域，通常利用光學、雷達等多個測量設備對飛行器跟蹤測量，再由采集到的測量數據給出飛行器高精度的軌跡參數估計。在此跟蹤測量過程中經常出現不完全測量的情況，即測量設備由于工作狀態不理想或受到環境干擾，在某些時段丟失大量測量數據，導致無法解算完整的軌跡參數，直接影響了后續的精度評估工作。如何利用不完全測量數據獲得精度盡量高的軌跡參數是測控數據處理亟須解決的難題。

飛行器軌跡參數估計通常采用基于樣條函數表示的數據融合方法[1-3]，在不完全測量情況下，該方法的估計模型呈現病態性，難以有效確定軌跡參數。文獻[4]利用誤差診斷方法，在完全測量段落獲取輔助測量設備低精度測元的系統誤差，扣除系統誤差后融合其他設備少量的高精度測元用于估計不完全測量段落的軌跡參數，獲得較高精度的估計結果，但該方法的前提是具有輔助測量設備并且該設備的測元完好，在實際應用中受到限制。文獻[5]利用飛行器軌跡的動力特性優化軌跡的表示樣條，在丟失部分測量數據的情況下，仍然能夠獲得較高精度的軌跡參數，但是該方法只適用于不完全測量段落較短的情況，當丟失數據的時間較長時，估計模型依然呈現病態性，估計結果并不理想。不完全測量問題本質上是病態逆問題，信號處理領域的稀疏優化方法被大量應用于病態逆問題的研究[6-9]，該方法充分利用先驗知識進行信號優化建模，得到信號的稀疏表示，通過確定少量參數即可實現真實信號的重構，從而縮小參數空間維數，將病態問題轉變為良態問題。值得注意的是，信號稀疏優化問題求解困難的根源在于信號的稀疏結構或支撐未知，很多方法的出發點都是要解決這一問題。而對于飛行器來說，其軌跡參數的稀疏結構在一定程度上是可建模的，可用樣條模型表示[2-3,10]。國內外一些學者將稀疏優化方法應用于樣條擬合問題的研究中。文獻[11]通過求解l1范數稀疏優化問題使擬合樣條的待求參數最小化，提高了擬合性能。文獻[12]建立了多分辨率樣條基函數集合，采用稀疏優化方法從集合中挑選出少量基函數，用其表示函數曲線能夠有效減少待求參數的數量。文獻[13]由插值樣條函數的l0范數稀疏優化獲得函數曲線的稀疏表示。

本文針對不完全測量條件下的飛行器軌跡參數估計問題，提出了一種基于稀疏優化的軌跡參數估計新方法。該方法將軌跡參數用B樣條函數表示，結合不完全測量數據建立軌跡參數估計的稀疏表示尋優模型，通過樣條節點的稀疏優化最大限度壓縮樣條節點數，從而降低參數空間的維數，緩解模型的病態性。

1 飛行器軌跡參數的稀疏表示尋優模型

(1)

(2)

y(t)=f(x(t))+ε(t)

其中，f(x(t))為測量真值，ε(t)為服從零均值多維正態分布的隨機誤差，其協方差矩陣記為Γ(t)。

測量數據的全部采樣時刻記為t1,t2,…,tM，對應的軌跡參數記為X=(x(t1),x(t2),…,x(tM))T。軌跡參數估計就是融合所有時刻的測量數據y(t1),y(t2),…,y(tM)給出X的高精度估計，從而確定飛行器完整的運動軌跡。對測量數據作歸一化處理，令

Y=(Γ-1/2(t1)y(t1),Γ-1/2(t2)y(t2),…,

Γ-1/2(tM)y(tM))T

E=(Γ-1/2(t1)ε(t1),Γ-1/2(t2)ε(t2),…,

Γ-1/2(tM)ε(tM))T

F(X)=(Γ-1/2(t1)f(x(t1)),Γ-1/2(t2)f(x(t2)),…,

Γ-1/2(tM)f(x(tM)))T

則有

Y=F(X)+E

(3)

其中，E服從多維標準正態分布。對X的準確估計采用如下非線性最小二乘估計：

在某些時刻點如果沒有足夠多的測量(3個位置測元和3個速度測元)解算6個軌跡參數，就出現了不完全測量的情況，由于測量的數量少于待估參數的數量，上述最小二乘估計模型是病態的，其解可能不唯一，即使能確定唯一解，但由于誤差傳播過大而使得求解結果在數值上不穩定，難以給出有效的軌跡參數估計。解決不完全測量問題一個重要途徑是通過物理機理分析、從歷史數據挖掘先驗信息等，尋求待估參數的稀疏表示，以降低模型解空間的維數，改善問題的病態性。通過對飛行器運動特性和大量實際軌跡數據的分析發現，采用B樣條函數表示軌跡參數能夠在保證精度的前提下大幅度減少待估參數。

將軌跡參數表示為如下n+1階B樣條函數：

(4)

則式(4)表示為：

x(t)=B(t,T)β

令B(T)=(B(t1,T),B(t2,T),…,B(tM,T))T，于是，式(3)變換為：

Y=F(B(T)β)+E

這樣軌跡參數X的估計轉化為樣條節點T和樣條系數β的估計：

(5)

待估參數T和β的數量顯著小于X的數量。

為進一步減少待估參數的數量，可將稀疏性約束(非零待估參數的數量最少)加入優化模型(5)中，構造如下稀疏表示尋優模型：

(6)

模型(6)在估計樣條系數β的同時，需確定樣條節點序列T，模型的非線性程度非常高，計算量太大，難以獲得全局最優解。可以先由粗略的軌跡參數預先確定樣條節點T*，然后再估計樣條系數。事實上，樣條節點數決定了樣條系數的數量，從而決定了待估參數的數量，因此，T*的選取應在保證較小的擬合誤差的前提下，盡量壓縮其數量，以改善估計模型的病態性。第2節將給出一種基于稀疏優化的樣條節點選取方法，能夠大幅度減少樣條節點數。確定T*后，只需估計樣條系數，即求解優化問題：

(7)

2 基于稀疏優化的樣條節點選取

樣條節點的選取應依據軌跡參數的動力學特性，在變化劇烈、可微性較差的軌跡段落(例如級間段)節點分布較為密集，在變化平緩的段落節點分布較為稀疏。遙測軌跡參數源于飛行器平臺系統的加速度表測量，其動力特征與飛行器的實際運動軌跡具有很好的一致性[5]，可作為確定樣條節點的粗略軌跡參數。下面只給出x方向的樣條節點確定方法，y和z方向的方法相同。

2.1 樣條節點的稀疏優化模型

給定區間[a,b]上的擴充分劃η:τ-n<τ-n+1

(8)

其中，bi為常值樣條系數，τi(i=1,2,…,N-1)稱為內節點序列。S(t)是[a,b]上的n-1階連續可微函數，n階導函數S(n)(t)為分段常值函數，其不連續點全部位于內節點序列中。關于S(n)(t)不連續的內節點稱為有效節點，連續的內節點稱為無效節點。無效節點是冗余的，即刪除無效節點后，不會改變樣條函數值。為減少樣條節點的數量，應去除無效節點，并且最大限度壓縮有效節點的數量。圖1為一個包含9個內節點的四階B樣條函數及其前三階導數的曲線圖，內節點中包括6個無效節點和3個有效節點，可以看出，三階導數為階梯函數，其不連續點位于3個有效節點處。

(a) 四階B樣條函數S(t)(a) Four order B-spline S(t)

(b) 一階導數S(1)(t)(b) The first order derivative S(1)(t)

(c) 二階導數S(2)(t)(c) The second order derivative S(2)(t)

(d) 三階導數S(3)(t)(d) The third order derivative S(3)(t)圖1 四階B樣條函數及其前三階導數Fig.1 A four order B-spline and its first three derivatives

θ=(b-n,b-n+1,…,bN-1)T

則D中的非零元對應于S(n)(t)的不連續點，即有效節點。為壓縮有效節點數，應在保證S(t)對r和v的擬合誤差不超過允許值的條件下，使D中的非零元素盡量少，即建立如下稀疏優化模型：

(9)

其中，δ1和δ2為設定的位置和速度擬合誤差門限，δ1>0，δ2>0。

為保證模型(9)存在可行解并且D具有稀疏性，應選取密集的初始樣條節點η，求解模型后，將無效節點從η中刪除，剩余的有效節點即為所需的優化樣條節點。本文采用節點加密方法確定初始樣條節點，步驟為：

步驟1：取較為密集的等距樣條節點序列η。

步驟3：對η′按距離聚類，把每一類各節點間的中點加入節點序列η中，轉到步驟2。

下面給出模型(9)的凸優化求解方法。

2.2 稀疏優化模型的求解方法

為求解模型(9)，需要獲得D關于θ的函數表達式。事實上，D是關于θ的線性函數，可將模型(9)轉化為線性凸優化問題，便能夠利用多項復雜度的凸優化方法完成模型求解。為確定D關于θ的函數表達式，給出如下定理：

定理n+1階B樣條函數S(t)的k階導數S(k)(t)(1≤k≤n)在[a,b]上滿足：

(10)

其中，

(11)

下面采用數學歸納法完成該定理的證明。

證明：Bi,n+1(t)滿足如下微分關系[2]：

(12)

當k=1時，由式(8)和式(12)得：

當t∈[a,b]時，B-n,n(t)≡0，BN,n(t)≡0，于是

定理成立。 假設當k=j時定理成立，當k=j+1時，由式(10)和式(12)得：

當t∈[a,b]時，B-n+j,n-j(t)≡0，BN,n-j(t)≡0，式(11)成立，則有：

Bi,n+1-(j+1)(t)

定理得證。

□

由定理得知：

而

于是對于i=1,…,N-1，有：

(13)

由式(11)和式(13)知，D與θ存在線性關系D=Cθ，其中C為(N-1)×(N+n)階的矩陣。這樣，模型(9)轉化為約束優化問題：

但上述最小化l0范數為NP難解問題，可用l1范數代替l0范數，轉化為線性凸優化問題：

(14)

可利用具有多項式復雜度的凸優化方法[14]求解。

3 仿真實驗與結果分析

飛行器的真實軌跡參數由動力學方程[15]積分產生，時長為94 s，軌跡中含有一個級間段，位于33.4 s附近。由9臺連續波雷達(測元為距離和徑向速度)和4臺光學設備(測元為方位角和俯仰角)對飛行器跟蹤測量，仿真測量數據由測量真值加隨機誤差構成，其中，連續波雷達的測距精度為8.3 m，測速精度為0.05 m/s，光學設備方位角和俯仰角的測量精度都為5″。設置兩個不完全測量場景：場景1的不完全測量段落為45 s—75 s的軌跡平穩段，該時間段只有2臺連續波雷達的2個距離和2個徑向速度測元，其余時間所有設備的測量數據完好；場景2的不完全測量段落也是45 s—75 s，該時間段只有1臺光學設備的1個俯仰角和1個方位角測元。

采用5階B樣條函數表示軌跡參數，利用2.1節中的節點加密方法確定x、y和z三個方向遙測軌跡的初始樣條節點，擬合誤差門限取為σ1=0.1 m，σ2=0.05 m/s，不完全測量段落的初始節點距取為0.2 s，其余段落的初始節點距取為 1 s。經過3～4次節點加密即可使擬合誤差小于門限值，x、y和z三個方向加密后的內節點數分別為237、239和213。圖2和圖3給出了x方向初始樣條節點的分布以及位置、速度的擬合誤差，可見，節點加密部分主要位于級間段，節點加密后的擬合誤差滿足門限設置要求。

圖2 x方向的初始樣條節點和位置擬合誤差Fig.2 Initial knots of spline and fitting errors of position for x-coordinate

由不完全測量段落(45 s—75 s)的初始樣條節點建立稀疏優化模型，模型中的擬合誤差門限取為δ1=0.1 m，δ2=0.05 m/s，利用CVX凸優化程序包[16]對模型進行求解，x、y和z三個方向剔除節點的|Cθ|閾值都取為1E-5，三個方向上優化后的樣條內節點數分別為95、96和62，明顯少于初始內節點數。圖4和圖5為x方向稀疏優化后的樣條節點分布以及位置、速度的擬合誤差，可以看出，不完全測量段落的樣條節點分布稀疏，擬合誤差未超出門限。

圖3 x方向的初始樣條節點和速度擬合誤差Fig.3 Initial knots of spline and fitting errors of velocity for x-coordinate

圖4 x方向的優化樣條節點和位置擬合誤差Fig.4 Optimized knots of spline and fitting errors of position for x-coordinate

圖5 x方向的優化樣條節點和速度擬合誤差Fig.5 Optimized knots of spline and fitting errors of velocity for x-coordinate

利用上述稀疏優化后的樣條節點結合測量數據估計軌跡參數。將本文方法的估計結果與文獻[5]提出的基于軌跡動力特征的估計方法(方法1)、基于多項式模型的估計方法(方法2)進行比較，其中，方法2采用三個多項式表示x、y和z三個方向的軌跡參數，由測量數據優化求解多項式系數，用以確定不完全測量段落的軌跡參數，方法2在完全測量段落的估計方法與方法1相同。

圖6和圖7分別是場景1位置和速度估計誤差的對比圖，可以看出，本文方法在不完全測量段落對軌跡參數的估計精度相對方法1和方法2都有大幅度的提高，方法1在不完全測量段落的位置和速度平均誤差分別為1.33 m和0.240 m/s，方法2的平均誤差為3.91 m和0.573 m/s，而本文方法的平均誤差只有0.38 m和0.067 m/s。場景2位置和速度估計誤差的對比見圖8和圖9，可見，本文方法相對方法1和方法2仍然能夠大幅度提高不完全測量段落的軌跡參數估計精度，本文方法在不完全測量段落的位置和速度平均誤差為0.74 m和0.117 m/s，而方法1的平均誤差達到了2.73 m和0.471 m/s，方法2的平均誤差達到了5.74 m和0.823 m/s。本文方法、方法1和方法2在不完全測量段落待估參數的數量分別為17、29和20。本文方法相對方法1能夠提高估計精度，是因為通過樣條節點的稀疏優化大幅度減少了待估參數的數量，緩解了模型的病態性；方法2待估參數的數量雖然也較少，但用多項式表示軌跡參數引入了較大的截斷誤差，估計精度低。場景2的估計誤差相對場景1較大，這是由于場景2在不完全測量段落的測量數據數量比場景1更少。

圖6 場景1位置估計誤差Fig.6 Errors of position estimates in scenario 1

圖7 場景1速度估計誤差Fig.7 Errors of velocity estimates in scenario 1

圖8 場景2位置估計誤差Fig.8 Errors of position estimates in scenario 2

圖9 場景2速度估計誤差Fig.9 Errors of velocity estimates in scenario 2

4 結論

研究不完全測量條件下的飛行器軌跡參數估計問題，提出了一種基于稀疏優化的軌跡參數估計新方法。利用B樣條函數給出軌跡參數的稀疏表示，建立了估計軌跡參數的尋優模型；依據樣條函數高階導數在節點處的不連續性，將表示樣條的節點劃分為有效節點和無效節點兩類，模型待估參數的數量取決于有效節點的數量；采用稀疏優化方法最大限度壓縮樣條的有效節點數，從而減少待估參數的數量，改善模型的病態性。仿真結果表明，在不完全測量段落，稀疏優化方法相對已有的基于軌跡動力特征的估計方法能夠大幅度提高軌跡參數的估計精度。本文的方法沒有考慮測量設備的系統誤差，下一步將研究不完全測量條件下的軌跡參數和系統誤差的聯合估計方法，實現系統誤差自校準。