飛行器軌跡的半?yún)?shù)估計(jì)方法

李振興,李 冬,劉 學(xué)

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

0 引言

飛行器試驗(yàn)鑒定一般采用光學(xué)、雷達(dá)等多臺(tái)套精密測(cè)量設(shè)備對(duì)飛行器跟蹤測(cè)量,然后對(duì)所有設(shè)備的測(cè)量數(shù)據(jù)作融合處理,獲得飛行器軌跡的高精度估計(jì),用于分離飛行器的制導(dǎo)工具誤差,進(jìn)而評(píng)定飛行器的命中精度[1-3]。隨著飛行器制導(dǎo)元器件精密度的不斷提高,為保證制導(dǎo)工具誤差的有效分離,對(duì)軌跡估計(jì)的精度要求越來(lái)越高。設(shè)備測(cè)量誤差導(dǎo)致軌跡的估計(jì)值偏離真實(shí)值,是影響軌跡估計(jì)精度的關(guān)鍵因素,研究能夠最大限度抑制測(cè)量誤差的軌跡估計(jì)方法對(duì)于提高飛行器軌跡的估計(jì)精度有著重要意義。

飛行器軌跡估計(jì)的精度主要受2種測(cè)量誤差影響,第1種誤差是由設(shè)備內(nèi)部熱噪聲、環(huán)境擾動(dòng)等引起的隨機(jī)誤差,在零值附近隨機(jī)波動(dòng),一般根據(jù)統(tǒng)計(jì)規(guī)律融合多個(gè)測(cè)量數(shù)據(jù)對(duì)其作補(bǔ)償抑制。第2種誤差是系統(tǒng)誤差,是由設(shè)備標(biāo)校不準(zhǔn)確、時(shí)間不同步、大氣折射修正殘差等引起的測(cè)量偏差。飛行器軌跡估計(jì)通常根據(jù)系統(tǒng)誤差的變化規(guī)律進(jìn)行參數(shù)化建模,例如建立常值模型、線性模型、周期模型等,采用參數(shù)估計(jì)方法同時(shí)估計(jì)軌跡和分離系統(tǒng)誤差,從而降低系統(tǒng)誤差對(duì)軌跡估計(jì)精度的影響[4-6]。然而,在實(shí)際測(cè)量中,受多種復(fù)雜誤差源的綜合影響,系統(tǒng)誤差往往不具有規(guī)律性,無(wú)法用參數(shù)模型對(duì)其進(jìn)行建模描述,而如果不考慮這種不可參數(shù)建模的系統(tǒng)誤差,勢(shì)必會(huì)降低軌跡估計(jì)的精度。

針對(duì)無(wú)法利用參數(shù)估計(jì)方法分離不可參數(shù)建模系統(tǒng)誤差的問(wèn)題,統(tǒng)計(jì)學(xué)界提出了半?yún)?shù)估計(jì)方法[7-9],用參數(shù)模型表示待估計(jì)的真實(shí)信號(hào),用非參數(shù)表示不可參數(shù)建模的系統(tǒng)誤差,建立既含有參數(shù)又含有非參數(shù)的半?yún)?shù)回歸模型,采用半?yún)?shù)估計(jì)方法同時(shí)獲得參數(shù)分量和非參數(shù)分量的估計(jì),可實(shí)現(xiàn)真實(shí)信號(hào)與系統(tǒng)誤差的有效分離,從而提高真實(shí)信號(hào)的估計(jì)精度。許多學(xué)者將半?yún)?shù)估計(jì)方法應(yīng)用于人造衛(wèi)星的軌道確定[10-12]。文獻(xiàn)[10]建立基于補(bǔ)償最小二乘的半?yún)?shù)聯(lián)合定軌模型,同時(shí)估計(jì)衛(wèi)星軌道參數(shù)和模型誤差。文獻(xiàn)[11]針對(duì)天基測(cè)量確定初軌問(wèn)題,提出了基于半?yún)?shù)模型的廣義正則化最小二乘估計(jì)方法,既有效抑制了不可參數(shù)建模系統(tǒng)誤差,又提高了初軌確定的穩(wěn)定性。文獻(xiàn)[12]理論上證明了半?yún)?shù)回歸補(bǔ)償最小二乘估計(jì)法優(yōu)于經(jīng)典最小二乘法。文獻(xiàn)[13]將半?yún)?shù)估計(jì)方法應(yīng)用于慣導(dǎo)離心機(jī)試驗(yàn)中的加速度計(jì)誤差模型辨識(shí),可同時(shí)辨識(shí)加速度計(jì)誤差模型參數(shù)和系統(tǒng)誤差。目前將半?yún)?shù)估計(jì)方法應(yīng)用于飛行器軌跡估計(jì)問(wèn)題的研究較少,文獻(xiàn)[14]利用多項(xiàng)式調(diào)制函數(shù)表示系統(tǒng)誤差,由軌跡估計(jì)后的殘差來(lái)估計(jì)多項(xiàng)式調(diào)制函數(shù)的系數(shù),該方法本質(zhì)上仍然是參數(shù)估計(jì)方法,多項(xiàng)式調(diào)制函數(shù)很難對(duì)實(shí)際應(yīng)用中的復(fù)雜系統(tǒng)誤差進(jìn)行精確表示,因此該方法對(duì)軌跡估計(jì)精度的提升效果有限。

文中針對(duì)飛行器軌跡估計(jì)的不可參數(shù)建模測(cè)量系統(tǒng)誤差的抑制問(wèn)題,提出一種半?yún)?shù)估計(jì)方法。該方法建立既包含參數(shù)待估量又包含非參數(shù)待估量的半?yún)?shù)回歸模型,參數(shù)待估量為飛行器軌跡的樣條表示系數(shù),非參數(shù)待估量為測(cè)量系統(tǒng)誤差,采用參數(shù)估計(jì)和非參數(shù)估計(jì)交替迭代的方式求解半?yún)?shù)回歸模型,參數(shù)估計(jì)利用非線性最小二乘法估計(jì)軌跡的樣條表示系數(shù),非參數(shù)估計(jì)利用小波分析方法從殘差數(shù)據(jù)中分離系統(tǒng)誤差。最后從理論上分析了半?yún)?shù)估計(jì)方法的精度,并通過(guò)仿真實(shí)驗(yàn)對(duì)文中方法的估計(jì)精度進(jìn)行了驗(yàn)證。

1 飛行器軌跡的參數(shù)估計(jì)

1.1 軌跡建模

(1)

根據(jù)式(1),飛行器軌跡的樣條表示模型為

x(t)=B(t)α

(2)

1.2 測(cè)量模型

(3)

(4)

式(4)中:β=(β1,β2,…,βp)T為p維系統(tǒng)誤差參數(shù)向量。

發(fā)射坐標(biāo)系向地心坐標(biāo)系的旋轉(zhuǎn)矩陣記為C,地心坐標(biāo)系向測(cè)站坐標(biāo)系的旋轉(zhuǎn)矩陣記為G,測(cè)站和發(fā)射原點(diǎn)在地心坐標(biāo)系下的矢量分別記為rs0=(xs0,ys0,zs0)T和rf0=(xf0,yf0,zf0)T,則飛行器的位置和速度從發(fā)射坐標(biāo)系到測(cè)站坐標(biāo)系的轉(zhuǎn)換關(guān)系表示為[1]

(5)

令y(t)為t時(shí)刻的所有測(cè)量參數(shù)組成的向量;根據(jù)式(3)和式(5),y(t)可表示為軌跡x(t)的函數(shù):

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

(6)

式(6)中:s(t)為測(cè)量的系統(tǒng)誤差向量,ε(t)為測(cè)量的隨機(jī)誤差向量,其協(xié)方差矩陣記為Σ(t)。將軌跡的參數(shù)模型(2)、系統(tǒng)誤差的參數(shù)模型(4)代入測(cè)量模型(6),可得到如下以樣條系數(shù)α和系統(tǒng)誤差參數(shù)β為未知參數(shù)的測(cè)量模型:

y(t)=z(t,α,β)+ε(t)

其中,z(t,α,β)=f(B(t)α)+s(t,β)。

1.3 參數(shù)估計(jì)

飛行器軌跡測(cè)量的采樣時(shí)刻記為t1,t2,…,tn,軌跡估計(jì)是根據(jù)所有采樣時(shí)刻測(cè)量數(shù)據(jù)y(t1),y(t2),…,y(tn)的融合值給出軌跡X=(x(t1),x(t2),…,x(tn))T的最優(yōu)估計(jì),從而獲得飛行器高精度的位置和速度。考慮到不同測(cè)站測(cè)量精度不一致,需利用測(cè)量隨機(jī)誤差的協(xié)方差矩陣對(duì)測(cè)量模型作歸一化處理[6],以實(shí)現(xiàn)不同測(cè)站多源數(shù)據(jù)的最佳融合:

Y=(Σ-1/2(t1)y(t1),…,Σ-1/2(tn)y(tn))T

Z(α,β)=(Σ-1/2(t1)z(t1,α,β),…,Σ-1/2(tn)z(tn,α,β))T

E=(Σ-1/2(t1)ε(t1),…,Σ-1/2(tn)ε(tn))T

其中,E服從均值為零協(xié)方差矩陣為單位矩陣的高斯分布,這樣便可建立以樣條系數(shù)和系統(tǒng)誤差參數(shù)為待估參數(shù)的非線性回歸模型:

Y=Z(α,β)+E

采用非線性最小二乘法估計(jì)樣條系數(shù)和系統(tǒng)誤差參數(shù):

可采用高斯牛頓迭代法[1]求解上述非線性最小二乘問(wèn)題。飛行器軌跡的估計(jì)為

其中,B=(B(t1),B(t2),…,B(tn))T。

2 飛行器軌跡的半?yún)?shù)估計(jì)

2.1 半?yún)?shù)估計(jì)

受測(cè)量設(shè)備狀態(tài)、環(huán)境因素的影響,測(cè)量系統(tǒng)誤差往往呈現(xiàn)復(fù)雜特性。例如圖1所示的某型光學(xué)設(shè)備俯仰角測(cè)量的系統(tǒng)誤差,對(duì)于這種復(fù)雜的系統(tǒng)誤差無(wú)法建立常值模型、線性模型、周期模型等參數(shù)模型,因此不能采用參數(shù)估計(jì)方法同時(shí)估計(jì)軌跡的樣條表示系數(shù)和分離系統(tǒng)誤差。以下將不可參數(shù)建模的系統(tǒng)誤差用非參數(shù)表示,建立既包含軌跡表示參數(shù)又包括非參數(shù)的半?yún)?shù)回歸模型,采用半?yún)?shù)估計(jì)方法求解半?yún)?shù)回歸模型,同時(shí)獲得軌跡和系統(tǒng)誤差的估計(jì)結(jié)果。

圖1 光學(xué)設(shè)備俯仰角測(cè)量的系統(tǒng)誤差

考慮不可參數(shù)建模系統(tǒng)誤差的影響,建立如下測(cè)量模型:

y(t)=f(B(t)α)+s(t)+ε(t)

利用采樣時(shí)刻t1,t2,…,tn的測(cè)量數(shù)據(jù)建立歸一化的半?yún)?shù)回歸模型:

Y=F(Bα)+S+E

(7)

其中

F(Bα)=(Σ-1/2(t1)f(B(t1)α),…,Σ-1/2(tn)f(B(tn)α))T

S=(Σ-1/2(t1)s(t1),…,Σ-1/2(tn)s(tn))T

模型(7)的待估量除了包括參數(shù)α,還包括非參數(shù)S,半?yún)?shù)估計(jì)方法需同時(shí)獲得參數(shù)α和非參數(shù)S的估計(jì)。

測(cè)量隨機(jī)誤差和系統(tǒng)誤差會(huì)在軌跡估計(jì)后的殘差數(shù)據(jù)中有所體現(xiàn)。其中,系統(tǒng)誤差通常表現(xiàn)為變化平緩的低頻信號(hào),而隨機(jī)誤差為變化劇烈的高頻信號(hào),可根據(jù)2種誤差的不同頻率特性采用小波分析方法[16]將兩者進(jìn)行分離,從而提取殘差中的低頻信號(hào),將其作為系統(tǒng)誤差的近似值,再將扣除系統(tǒng)誤差后的測(cè)量數(shù)據(jù)用于軌跡估計(jì),可進(jìn)一步改善軌跡估計(jì)的精度。

半?yún)?shù)估計(jì)采用軌跡估計(jì)和系統(tǒng)誤差分離交替迭代的方式實(shí)現(xiàn),軌跡估計(jì)采用非線性最小二乘法,軌跡估計(jì)完成后,采用小波分析方法對(duì)殘差數(shù)據(jù)作系統(tǒng)誤差分離,利用分離結(jié)果對(duì)測(cè)量數(shù)據(jù)作修正,由修正后的測(cè)量數(shù)據(jù)再進(jìn)行軌跡估計(jì),進(jìn)一步修正軌跡估計(jì)結(jié)果,這一過(guò)程反復(fù)迭代,直到殘差數(shù)據(jù)中不含系統(tǒng)誤差為止。半?yún)?shù)估計(jì)方法的原理如圖2所示。

圖2 半?yún)?shù)估計(jì)方法的原理圖

具體步驟如下:

步驟4檢驗(yàn)殘差ΔY是否為白噪聲,如果是,則轉(zhuǎn)到步驟6,否則,令j=j+1,轉(zhuǎn)到步驟5。

步驟5利用Daubechies7小波從殘差ΔY中分離出系統(tǒng)誤差ΔS(j)和測(cè)量噪聲,令Y=Y-ΔS(j),i=0,ω=1,轉(zhuǎn)到步驟2,重新利用高斯牛頓迭代法估計(jì)α。

2.2 估計(jì)精度分析

于是

即

BTAT(Y-F(Bαp))=0

(8)

(9)

由式(8)和式(9)可得:

于是,參數(shù)估計(jì)方法的軌跡估計(jì)精度為

B(BTATAB)-1BT+QQT

(10)

同理可得:

B(BTATAB)-1BT

(11)

即半?yún)?shù)估計(jì)方法得到的軌跡精度優(yōu)于忽略系統(tǒng)誤差的參數(shù)估計(jì)方法的精度。

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

利用仿真實(shí)驗(yàn)檢驗(yàn)半?yún)?shù)估計(jì)方法的性能。飛行器的真實(shí)軌跡由6自由動(dòng)力學(xué)仿真模型產(chǎn)生,時(shí)長(zhǎng)為94 s,采樣率為20 Hz。5臺(tái)雷達(dá)(測(cè)量數(shù)據(jù)為距離和距離變化率)和2臺(tái)光學(xué)設(shè)備(測(cè)量數(shù)據(jù)為方位角和俯仰角)對(duì)飛行器跟蹤測(cè)量。利用測(cè)量設(shè)備的站址和飛行器的軌跡真值生成測(cè)量真值,在測(cè)量真值的基礎(chǔ)上通過(guò)疊加測(cè)量誤差生成仿真測(cè)量數(shù)據(jù),其中,距離、方位角和俯仰角的測(cè)量誤差含有不可參數(shù)建模系統(tǒng)誤差的實(shí)測(cè)值,距離變化率測(cè)量誤差為均方根為0.05 m/s的隨機(jī)誤差。

分別采用參數(shù)估計(jì)方法和半?yún)?shù)估計(jì)方法對(duì)飛行器軌跡進(jìn)行估計(jì),其中,參數(shù)估計(jì)方法忽略距離、方位角和俯仰角的不可參數(shù)建模系統(tǒng)誤差的影響,而半?yún)?shù)估計(jì)方法對(duì)系統(tǒng)誤差進(jìn)行分離。

通過(guò)比較軌跡估計(jì)值與真實(shí)值的差異得到軌跡的估計(jì)誤差,由估計(jì)誤差的大小分析估計(jì)方法的精度。t時(shí)刻的位置和速度估計(jì)誤差分別定義為

圖3和圖4分別是位置和速度估計(jì)誤差隨時(shí)間變化的曲線圖。可見(jiàn),半?yún)?shù)估計(jì)方法相對(duì)參數(shù)估計(jì)方法關(guān)于位置和速度的估計(jì)精度都有較大幅度提高,參數(shù)估計(jì)方法對(duì)位置和速度估計(jì)的平均誤差分別為1.29 m和0.095 m/s,而半?yún)?shù)估計(jì)方法對(duì)位置和速度估計(jì)的平均誤差分別為0.62 m和0.055 m/s,與參數(shù)估計(jì)方法相比,半?yún)?shù)估計(jì)方法關(guān)于位置和速度的估計(jì)誤差分別減少了51.9%和42.1%。

圖3 位置估計(jì)誤差

圖4 速度估計(jì)誤差

圖5是其中1臺(tái)雷達(dá)的距離測(cè)量誤差及半?yún)?shù)估計(jì)方法獲得的系統(tǒng)誤差估計(jì)結(jié)果。圖6是半?yún)?shù)估計(jì)方法分離出的隨機(jī)誤差。可見(jiàn),系統(tǒng)誤差估計(jì)值的變化平緩,頻率較低,與測(cè)量誤差的變化趨勢(shì)相吻合,分離出的隨機(jī)誤差頻率較高,接近于白噪聲。圖7是其中1臺(tái)光學(xué)設(shè)備的俯仰角測(cè)量誤差及系統(tǒng)誤差估計(jì)結(jié)果。圖8是對(duì)應(yīng)的隨機(jī)誤差分離結(jié)果,可見(jiàn),半?yún)?shù)估計(jì)方法對(duì)光學(xué)設(shè)備系統(tǒng)誤差的分離效果也較好。圖7和圖8中的誤差曲線在45～75 s出現(xiàn)空白,這是因?yàn)樵摃r(shí)段飛行器進(jìn)入云層,光學(xué)設(shè)備丟失目標(biāo),沒(méi)有獲得測(cè)量數(shù)據(jù)。

圖5 雷達(dá)距離測(cè)量誤差及系統(tǒng)誤差估計(jì)結(jié)果

圖6 雷達(dá)距離測(cè)量隨機(jī)誤差的分離結(jié)果

圖7 光學(xué)設(shè)備俯仰角測(cè)量誤差及系統(tǒng)誤差估計(jì)結(jié)果

圖8 光學(xué)設(shè)備俯仰角測(cè)量隨機(jī)誤差的分離結(jié)果

4 結(jié)論

1) 提出了一種基于半?yún)?shù)回歸的飛行器軌跡估計(jì)新方法,可同時(shí)估計(jì)軌跡參數(shù)和分離不可參數(shù)建模的測(cè)量系統(tǒng)誤差,有效抑制了測(cè)量系統(tǒng)誤差對(duì)軌跡估計(jì)的影響,提高了軌跡估計(jì)精度。

2) 得到了半?yún)?shù)估計(jì)方法估計(jì)精度的理論公式,從理論上證明了半?yún)?shù)估計(jì)方法的精度優(yōu)于傳統(tǒng)的參數(shù)估計(jì)方法。

3) 通過(guò)仿真實(shí)驗(yàn)檢驗(yàn)了半?yún)?shù)估計(jì)方法的軌跡估計(jì)精度,結(jié)果表明,半?yún)?shù)估計(jì)方法對(duì)飛行器位置和速度的估計(jì)精度相對(duì)參數(shù)估計(jì)方法都有大幅度的提高,系統(tǒng)誤差分離結(jié)果準(zhǔn)確可靠,有效解決了不可參數(shù)建模系統(tǒng)誤差的處理難題,在飛行器試驗(yàn)鑒定領(lǐng)域有重要應(yīng)用價(jià)值。