李 冬，劉 學(xué)

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

在航天測(cè)控領(lǐng)域，通常利用光學(xué)、雷達(dá)等多個(gè)測(cè)量設(shè)備對(duì)飛行器跟蹤測(cè)量，再由采集到的測(cè)量數(shù)據(jù)給出飛行器高精度的軌跡參數(shù)估計(jì)。在此跟蹤測(cè)量過(guò)程中經(jīng)常出現(xiàn)不完全測(cè)量的情況，即測(cè)量設(shè)備由于工作狀態(tài)不理想或受到環(huán)境干擾，在某些時(shí)段丟失大量測(cè)量數(shù)據(jù)，導(dǎo)致無(wú)法解算完整的軌跡參數(shù)，直接影響了后續(xù)的精度評(píng)估工作。如何利用不完全測(cè)量數(shù)據(jù)獲得精度盡量高的軌跡參數(shù)是測(cè)控?cái)?shù)據(jù)處理亟須解決的難題。

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

本文針對(duì)不完全測(cè)量條件下的飛行器軌跡參數(shù)估計(jì)問(wèn)題，提出了一種基于稀疏優(yōu)化的軌跡參數(shù)估計(jì)新方法。該方法將軌跡參數(shù)用B樣條函數(shù)表示，結(jié)合不完全測(cè)量數(shù)據(jù)建立軌跡參數(shù)估計(jì)的稀疏表示尋優(yōu)模型，通過(guò)樣條節(jié)點(diǎn)的稀疏優(yōu)化最大限度壓縮樣條節(jié)點(diǎn)數(shù)，從而降低參數(shù)空間的維數(shù)，緩解模型的病態(tài)性。

1 飛行器軌跡參數(shù)的稀疏表示尋優(yōu)模型

(1)

(2)

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

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

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

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服從多維標(biāo)準(zhǔn)正態(tài)分布。對(duì)X的準(zhǔn)確估計(jì)采用如下非線性最小二乘估計(jì)：

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

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

(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

這樣軌跡參數(shù)X的估計(jì)轉(zhuǎn)化為樣條節(jié)點(diǎn)T和樣條系數(shù)β的估計(jì)：

(5)

待估參數(shù)T和β的數(shù)量顯著小于X的數(shù)量。

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

(6)

模型(6)在估計(jì)樣條系數(shù)β的同時(shí)，需確定樣條節(jié)點(diǎn)序列T，模型的非線性程度非常高，計(jì)算量太大，難以獲得全局最優(yōu)解?？梢韵扔纱致缘能壽E參數(shù)預(yù)先確定樣條節(jié)點(diǎn)T*，然后再估計(jì)樣條系數(shù)。事實(shí)上，樣條節(jié)點(diǎn)數(shù)決定了樣條系數(shù)的數(shù)量，從而決定了待估參數(shù)的數(shù)量，因此，T*的選取應(yīng)在保證較小的擬合誤差的前提下，盡量壓縮其數(shù)量，以改善估計(jì)模型的病態(tài)性。第2節(jié)將給出一種基于稀疏優(yōu)化的樣條節(jié)點(diǎn)選取方法，能夠大幅度減少樣條節(jié)點(diǎn)數(shù)。確定T*后，只需估計(jì)樣條系數(shù)，即求解優(yōu)化問(wèn)題：

(7)

2 基于稀疏優(yōu)化的樣條節(jié)點(diǎn)選取

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

2.1 樣條節(jié)點(diǎn)的稀疏優(yōu)化模型

給定區(qū)間[a,b]上的擴(kuò)充分劃η:τ-n<τ-n+1

(8)

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

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

(b) 一階導(dǎo)數(shù)S(1)(t)(b) The first order derivative S(1)(t)

(c) 二階導(dǎo)數(shù)S(2)(t)(c) The second order derivative S(2)(t)

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

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

則D中的非零元對(duì)應(yīng)于S(n)(t)的不連續(xù)點(diǎn)，即有效節(jié)點(diǎn)。為壓縮有效節(jié)點(diǎn)數(shù)，應(yīng)在保證S(t)對(duì)r和v的擬合誤差不超過(guò)允許值的條件下，使D中的非零元素盡量少，即建立如下稀疏優(yōu)化模型：

(9)

其中，δ1和δ2為設(shè)定的位置和速度擬合誤差門(mén)限，δ1>0，δ2>0。

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

步驟1：取較為密集的等距樣條節(jié)點(diǎn)序列η。

步驟3：對(duì)η′按距離聚類(lèi)，把每一類(lèi)各節(jié)點(diǎn)間的中點(diǎn)加入節(jié)點(diǎn)序列η中，轉(zhuǎn)到步驟2。

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

2.2 稀疏優(yōu)化模型的求解方法

為求解模型(9)，需要獲得D關(guān)于θ的函數(shù)表達(dá)式。事實(shí)上，D是關(guān)于θ的線性函數(shù)，可將模型(9)轉(zhuǎn)化為線性凸優(yōu)化問(wèn)題，便能夠利用多項(xiàng)復(fù)雜度的凸優(yōu)化方法完成模型求解。為確定D關(guān)于θ的函數(shù)表達(dá)式，給出如下定理：

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

(10)

其中，

(11)

下面采用數(shù)學(xué)歸納法完成該定理的證明。

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

(12)

當(dāng)k=1時(shí)，由式(8)和式(12)得：

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

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

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

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

定理得證。

□

由定理得知：

而

于是對(duì)于i=1,…,N-1，有：

(13)

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

但上述最小化l0范數(shù)為NP難解問(wèn)題，可用l1范數(shù)代替l0范數(shù)，轉(zhuǎn)化為線性凸優(yōu)化問(wèn)題：

(14)

可利用具有多項(xiàng)式復(fù)雜度的凸優(yōu)化方法[14]求解。

3 仿真實(shí)驗(yàn)與結(jié)果分析

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

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

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

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

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

圖4 x方向的優(yōu)化樣條節(jié)點(diǎn)和位置擬合誤差Fig.4 Optimized knots of spline and fitting errors of position for x-coordinate

圖5 x方向的優(yōu)化樣條節(jié)點(diǎn)和速度擬合誤差Fig.5 Optimized knots of spline and fitting errors of velocity for x-coordinate

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

圖6和圖7分別是場(chǎng)景1位置和速度估計(jì)誤差的對(duì)比圖，可以看出，本文方法在不完全測(cè)量段落對(duì)軌跡參數(shù)的估計(jì)精度相對(duì)方法1和方法2都有大幅度的提高，方法1在不完全測(cè)量段落的位置和速度平均誤差分別為1.33 m和0.240 m/s，方法2的平均誤差為3.91 m和0.573 m/s，而本文方法的平均誤差只有0.38 m和0.067 m/s。場(chǎng)景2位置和速度估計(jì)誤差的對(duì)比見(jiàn)圖8和圖9，可見(jiàn)，本文方法相對(duì)方法1和方法2仍然能夠大幅度提高不完全測(cè)量段落的軌跡參數(shù)估計(jì)精度，本文方法在不完全測(cè)量段落的位置和速度平均誤差為0.74 m和0.117 m/s，而方法1的平均誤差達(dá)到了2.73 m和0.471 m/s，方法2的平均誤差達(dá)到了5.74 m和0.823 m/s。本文方法、方法1和方法2在不完全測(cè)量段落待估參數(shù)的數(shù)量分別為17、29和20。本文方法相對(duì)方法1能夠提高估計(jì)精度，是因?yàn)橥ㄟ^(guò)樣條節(jié)點(diǎn)的稀疏優(yōu)化大幅度減少了待估參數(shù)的數(shù)量，緩解了模型的病態(tài)性；方法2待估參數(shù)的數(shù)量雖然也較少，但用多項(xiàng)式表示軌跡參數(shù)引入了較大的截?cái)嗾`差，估計(jì)精度低。場(chǎng)景2的估計(jì)誤差相對(duì)場(chǎng)景1較大，這是由于場(chǎng)景2在不完全測(cè)量段落的測(cè)量數(shù)據(jù)數(shù)量比場(chǎng)景1更少。

圖6 場(chǎng)景1位置估計(jì)誤差Fig.6 Errors of position estimates in scenario 1

圖7 場(chǎng)景1速度估計(jì)誤差Fig.7 Errors of velocity estimates in scenario 1

圖8 場(chǎng)景2位置估計(jì)誤差Fig.8 Errors of position estimates in scenario 2

圖9 場(chǎng)景2速度估計(jì)誤差Fig.9 Errors of velocity estimates in scenario 2

4 結(jié)論

研究不完全測(cè)量條件下的飛行器軌跡參數(shù)估計(jì)問(wèn)題，提出了一種基于稀疏優(yōu)化的軌跡參數(shù)估計(jì)新方法。利用B樣條函數(shù)給出軌跡參數(shù)的稀疏表示，建立了估計(jì)軌跡參數(shù)的尋優(yōu)模型；依據(jù)樣條函數(shù)高階導(dǎo)數(shù)在節(jié)點(diǎn)處的不連續(xù)性，將表示樣條的節(jié)點(diǎn)劃分為有效節(jié)點(diǎn)和無(wú)效節(jié)點(diǎn)兩類(lèi)，模型待估參數(shù)的數(shù)量取決于有效節(jié)點(diǎn)的數(shù)量；采用稀疏優(yōu)化方法最大限度壓縮樣條的有效節(jié)點(diǎn)數(shù)，從而減少待估參數(shù)的數(shù)量，改善模型的病態(tài)性。仿真結(jié)果表明，在不完全測(cè)量段落，稀疏優(yōu)化方法相對(duì)已有的基于軌跡動(dòng)力特征的估計(jì)方法能夠大幅度提高軌跡參數(shù)的估計(jì)精度。本文的方法沒(méi)有考慮測(cè)量設(shè)備的系統(tǒng)誤差，下一步將研究不完全測(cè)量條件下的軌跡參數(shù)和系統(tǒng)誤差的聯(lián)合估計(jì)方法，實(shí)現(xiàn)系統(tǒng)誤差自校準(zhǔn)。