王奕波,劉雪峰,李 玫

(中國人民解放軍63801部隊,四川 西昌 615042)

由于系統(tǒng)誤差和隨機誤差的存在[1],航天測量任務(wù)中任何設(shè)備輸出的目標測量數(shù)據(jù)都不可能完全準確,而獲取目標的真實飛行彈道,對于實時引導(dǎo)數(shù)據(jù)生成[2]、設(shè)備測量性能分析[3]和目標飛行控制[4]等工作尤為重要,因此一般采用相關(guān)的數(shù)據(jù)處理方法,對不同設(shè)備的測量數(shù)據(jù)進行融合,從而實現(xiàn)目標真實彈道的估計。目前的彈道實時融合方法[5]主要有3種：基于選優(yōu)的方法[6]、樣條約束(spline restraint,SR)法[7]、無跡卡爾曼濾波(unscented Kalman filter,UKF)法[8]。選優(yōu)法的基本思路是采用如神經(jīng)網(wǎng)絡(luò)[9]等方法對各設(shè)備測量數(shù)據(jù)的可靠性進行量化,選取可靠性最高的數(shù)據(jù),該方法存在的問題是其估計結(jié)果總是某個設(shè)備的測量數(shù)據(jù),因此估計精度對設(shè)備本身的測量精度、測站分布依賴性較強,況且可靠性量化不可能完全正確;SR方法利用一段時間內(nèi)所有設(shè)備的多幀測量數(shù)據(jù),采用樣條函數(shù)對彈道位置進行估計,但由于該方法本質(zhì)上是對目標真實位置在多幀測量數(shù)據(jù)下的最小二乘(least squares,LS)估計[10],因此對異常測量數(shù)據(jù)的抗干擾性差;UKF方法對于高斯噪聲過程[11]具有較強的濾波性能,但依賴理論彈道方程[12],其估計精度易受目標真實飛行軌跡偏離理論彈道的影響,并且其估計精度對初值選取較為敏感。

本文基于測量誤差主要為高斯噪聲這一基本依據(jù)[13],證明了以服從高斯分布的有限樣本構(gòu)造的高斯核密度估計(kernel density estimation,KDE)函數(shù)[14]的最大值點為產(chǎn)生該組樣本高斯分布期望的無偏估計量這一結(jié)論,提出了一種基于高斯KDE函數(shù)最大值點跟蹤的目標彈道實時融合方法,從而將彈道融合問題轉(zhuǎn)化為求該函數(shù)最大值點的問題,并通過理論分析與對比實驗證明了本文方法的優(yōu)勢。

1 基于高斯核密度估計函數(shù)最大值點跟蹤的方法

1.1 基于核函數(shù)最大值點的估計方法

(1)

證畢。

(2)

1.2 核函數(shù)最大值點跟蹤

1.2.1 跟蹤方法描述

求取ρ(X)最大值點的直接思路是先求出ρ(X)的全部駐點,函數(shù)值最大的駐點則為最大值點。但ρ(X)的駐點約束方程ρ(X)=0為超越方程,無法直接求得精確解,因此只能采用一些數(shù)值分析方法求得近似解。傳統(tǒng)方法是采用迭代的思路,即將每組測量數(shù)據(jù)作為初值進行迭代運算,取核函數(shù)值最大的迭代結(jié)果作為最終估計結(jié)果。該方法對樣本點有一定的聚類[15]能力,能有效求出核函數(shù)值最大的位置,但存在耗時較長的問題：一方面是由于該方法采用相鄰兩步迭代結(jié)果的歐式距離閾值作為迭代終止判據(jù),對于不同的測量數(shù)據(jù)分布,其迭代步數(shù)是不定的,較小的閾值能夠保證收斂精度,但以犧牲時間成本為代價;另一方面由于該方法需要對每個都進行迭代運算,其時間消耗與設(shè)備數(shù)量N相關(guān)性較大。本文將在2.2節(jié)的實驗中進一步說明該問題。

針對上述問題,本文提出了一種核函數(shù)最大值點跟蹤的方法：首先根據(jù)每個設(shè)備當(dāng)前幀與前一幀的測量數(shù)據(jù)求出若干平移矢量;再以每個平移矢量為基準對上一幀的估計結(jié)果進行平移采樣,從而得到若干采樣結(jié)果;取當(dāng)前幀核函數(shù)值最大的若干采樣結(jié)果作為核函數(shù)最大值點的近鄰,采用一階泰勒估計的方法,計算出若干核函數(shù)最大值點的估計位置;最后取核函數(shù)值最大的位置,作為當(dāng)前幀的估計結(jié)果。具體步驟如下。

(3)

(4)

式中：r≠l,l=1,2,…,S。

1.2.2 跟蹤有效性分析

(5)

(6)

(7)

(8)

表1 PQ的計算結(jié)果Table 1 Calculation results of PQ

1.2.3 時間復(fù)雜度分析

下面對本文跟蹤方法的耗時情況進行分析,設(shè)求取一次由N組測量數(shù)據(jù)構(gòu)造的核函數(shù)值的平均耗時為t(N),本文方法耗時主要在步驟③、④、⑤中求核函數(shù)值的部分,步驟③的耗時為Nt(N);步驟④需計算S個點的一階和二階方向?qū)?shù),共9個方向,其耗時為9St(N);步驟⑤的耗時為St(N)。因此,本文方法總的理論耗時為(N+10S)t(N)。而經(jīng)典迭代法的耗時為TNt(N),T為平均迭代步數(shù),最后選取最大核函數(shù)值的迭代結(jié)果時需求取全部迭代結(jié)果的核函數(shù)值,耗時為Nt(N),因此總的耗時為N(T+1)t(N)。由于T≥1,S為常數(shù),當(dāng)N足夠大時,本文方法的耗時相對較低,且隨N的增大,本文方法相對經(jīng)典迭代法節(jié)省的時間成本越多。

需要說明的是,在飛行初始段,參與測量的設(shè)備數(shù)量較少,可以運用經(jīng)典迭代法,不會導(dǎo)致較高的時間消耗;但一段時間后,當(dāng)設(shè)備數(shù)量達到一定程度(由上文分析,在工程應(yīng)用中可設(shè)定當(dāng)設(shè)備數(shù)量達到5個時),采用本文跟蹤方法,可以同時保證估計精度和實時性,將在第2.2節(jié)的實驗中進一步證明。

2 實驗分析

2.1 條件與參數(shù)說明

本文實驗數(shù)據(jù)的生成方式如下：

①目標真實飛行軌跡。選用3條航天試驗任務(wù)的歷史彈道(每條彈道總計15 920幀)反演目標的真實飛行軌跡。

②測站分布。假設(shè)發(fā)射點的坐標為(0 m,0 m,0 m),根據(jù)我國領(lǐng)土的實際范圍,在以發(fā)射點坐標為基準,從(-2×106m,-2×106m,-2×103m)到(2×106m,2×106m,2×103m)的范圍內(nèi)按均勻分布隨機產(chǎn)生N個坐標模擬設(shè)備分布情況。

④異常測量數(shù)據(jù)。由于異常測量數(shù)據(jù)的統(tǒng)計特性難以預(yù)知,以目標真實位置為中心、從(0 km,-180°,-90°)到(10 000 km,180°,90°)的范圍內(nèi),按均勻分布隨機生成Z=N-M個離群點。

⑤硬件環(huán)境。計算機平臺內(nèi)存1 GB、CPU主頻3.0 GHz(雙核心)。

⑥其他實驗條件。下述所有實驗統(tǒng)計的各項指標(誤差均值、幀平均耗時等)均為對3條彈道分別運行100次實驗(總計4 776 000幀)的平均結(jié)果;為保證不同方法誤差比較的同一性,幀誤差、誤差均值、誤差均方差基于發(fā)射點進行計算;在實驗中任意時刻的測量協(xié)方差計算方法采用經(jīng)典的變量差分法[16];由于存在很多種選優(yōu)法,下述實驗中假定采用的選優(yōu)法是完全理想的,即每一幀都選取測量誤差最小的數(shù)據(jù),在此基礎(chǔ)上與本文方法進行比較,以說明本文方法的高精度優(yōu)勢。

2.2 跟蹤參數(shù)S的選取及有效性驗證

本文跟蹤方法涉及的參數(shù)S影響耗時,如果S設(shè)定得過高則耗時較長,過低則導(dǎo)致近鄰丟失,從而降低跟蹤精度。為對跟蹤精度進行評價,本文計算不同N和S條件下跟蹤結(jié)果與目標真實位置的平均偏差(設(shè)置離群點數(shù)量Z=0),結(jié)果如表2所示。表中,R為距離,A為方位角,E為俯仰角。

表2 不同N和S條件下的跟蹤結(jié)果與真實位置平均偏差Table 2 MD between the tracking results and the pinpoint under different N and S values

表2的計算結(jié)果表明,當(dāng)S≥6時估計精度基本不隨S的增大而提高,此時估計誤差是較低的(在第2.3節(jié)的對比實驗中進一步證明)。根據(jù)1.2.2節(jié)的分析,只要S足夠大就能保證至少存在一個平移采樣結(jié)果為KDE函數(shù)最大值點的近鄰,由于近鄰本身是從核函數(shù)值最大的S個平移采樣結(jié)果中選取的,再增加S的值只會將核函數(shù)值較小的采樣結(jié)果包含進來,而按照KDE函數(shù)的分布特性,核函數(shù)值越小與核函數(shù)最大值點的位置就越遠,因此再增加S的值不會顯著提高估計精度。而根據(jù)第1.2.3節(jié)的分析,S對耗時的影響較大,因此在工程應(yīng)用中不宜取較大的值。當(dāng)N=20和S=6的條件下幀平均耗時約為3.367 m/s(約297 Hz/s,比目前航天測量任務(wù)的數(shù)據(jù)融合頻率高一個數(shù)量級)。較小的S即可保證必存在一個采樣結(jié)果為最大值點的近鄰,跟蹤精度就可以保證,而S的取值不必隨N的增加而增大,因此同時保證了實時性。

為進一步證明本文方法的有效性,在N取值不同的條件下與經(jīng)典迭代法進行比較,結(jié)果如表3所示。實驗中設(shè)定迭代收斂條件為相鄰兩步迭代結(jié)果的歐式距離小于1 m,該條件下迭代法的距離收斂精度近似為1 m,角度收斂精度近似為0.001°。

表3表明,在不同的設(shè)備數(shù)量條件下,本文方法與經(jīng)典迭代法的估計精度近似相等,說明本文跟蹤方法能夠有效估計核函數(shù)最大值點。另外,本文對這2種方法在不同設(shè)備數(shù)量條件下的幀平均耗時情況進行了統(tǒng)計,如圖1。當(dāng)N較小時,2種方法的幀平均耗時均較低,但隨N的增大,由于經(jīng)典迭代法需將每組測量數(shù)據(jù)作為迭代初值進行運算,其時間消耗與N的相關(guān)性較大;隨N的增大,本文方法與經(jīng)典迭代法的幀平均耗時之差也越大,說明隨設(shè)備數(shù)量增多,本文方法相對經(jīng)典迭代法節(jié)省的時間成本越多。

表3 與經(jīng)典迭代法的估計精度比較Table 3 Comparison of MD with the classic iterative method

圖1 與經(jīng)典迭代法的時間消耗情況比較Fig.1 Comparison of the time-consuming test with the classic iterative method

2.3 對比實驗

①抗干擾性測試。該項測試用于評價離群點對本文方法估計精度的影響。在給定N=20的條件下,設(shè)定不同的離群點數(shù)量Z,統(tǒng)計本文方法、選優(yōu)法、UKF方法(參數(shù)選取σ點數(shù)量為3,α=0.5,β=2,λ=0)、SR方法的誤差均值,結(jié)果如圖2所示。

圖2 抗干擾性實驗的對比結(jié)果Fig.2 Comparison of the noise-tolerate capability test with other methods

由圖2可知,當(dāng)不存在離群點時,4種方法均能保持較高的估計精度,而一旦出現(xiàn)離群點,SR方法的估計精度迅速下降,原因是該方法將當(dāng)前幀及歷史多幀的測量數(shù)據(jù)采用LS方法進行估計,而LS方法的數(shù)學(xué)本質(zhì)是求解一個與所有測量數(shù)據(jù)距離之和最小的點作為估計結(jié)果,因此SR方法易受離群點位置及數(shù)量的影響,且離群點分布越遠,對估計結(jié)果的影響越大,這與本文方法恰好相反。選優(yōu)法只有在完全理想的情況下,即選優(yōu)完全正確的情況下,其估計精度才和本文方法相當(dāng)。UKF方法具有一定的抗干擾性能,這是由于該方法利用了歷史先驗信息,即理論彈道方程,因此具有一定的抗噪能力,但由于該方法依賴于理論彈道方程,如果目標飛行偏離理論彈道,UKF方法的估計精度會迅速降低,如圖3所示,以存在俯仰偏差為例。而本文方法不依賴理論彈道方程,不受目標真實飛行是否偏離理論彈道的影響。

圖3 UKF方法對理論彈道方程的敏感度測試Fig.3 Initial value sensitivity test for UKF

②初值敏感度測試。本文方法和UKF方法同屬于跟蹤方法,該項測試用于評價第一幀數(shù)據(jù)對這2種方法跟蹤精度的影響。在實際航天發(fā)射任務(wù)中的初始段,只有少量設(shè)備(如遙測雷達)參與測量,因此假設(shè)第一幀數(shù)據(jù)只由一臺設(shè)備產(chǎn)生,在第一幀數(shù)據(jù)與目標位置真值存在不同偏差的條件下,比較本文方法和UKF方法的跟蹤性能(設(shè)置設(shè)備數(shù)量N=20,離群點數(shù)量Z=0),如圖4所示。

圖4表明,當(dāng)初值存在一定偏度時(與目標真實位置距離偏差為10 000 m,方位角、俯仰角偏差均為0.2°),本文方法和UKF方法均具有一定的收斂性能。而在相同的初值條件下,本文方法的收斂精度更高,如表4所示。表中,ξ為不同的初值偏度系數(shù),所模擬的距離、方位、俯仰初值偏度步進分別為4 000 m,0.08°,0.08°,即所模擬的距離、方位、俯仰初值分別為對應(yīng)真值加上4 000ξ,0.08ξ和0.08ξ。本文方法和UKF方法的共同點是均采用了多點采樣的方式,UKF方法對采樣點進行加權(quán)求和,利用當(dāng)前幀實測結(jié)果對加權(quán)求和結(jié)果進行修正,因此UKF方法的收斂精度與初值偏度具有較大的相關(guān)性。而本文跟蹤方法是在多個采樣點中選取一個目標真實位置的近鄰進行估計,根據(jù)第1.2.2節(jié)的分析,在若干采樣點中至少存在1個近鄰的概率與前一幀的估計結(jié)果并無較大關(guān)系,因此本文方法具有較高的收斂精度。

表4 不同初值偏度條件下本文方法與UKF方法的誤差均值Table 4 Comparison of MD with UKF under different initial values

圖4 本文方法和UKF方法的跟蹤情況示例Fig.4 Comparison of MD with UKF under specific initial value

③伸縮性能測試。該項測試用于評價本文方法在不同設(shè)備數(shù)量和不同測站分布條件下的估計精度,設(shè)置離群點數(shù)量Z=0,統(tǒng)計不同設(shè)備數(shù)量(如圖5所示)和不同測站分布(如圖6所示,設(shè)置設(shè)備數(shù)量N=20)條件下本文方法、理想選優(yōu)法、UKF方法、SR方法的估計精度。

圖5表明,一方面本文方法在不同設(shè)備數(shù)量條件下的估計精度總是高于UKF方法和SR方法,而選優(yōu)法在完全理想的情況下,估計精度才和本文方法相當(dāng);另一方面,本文方法的數(shù)據(jù)曲線較為平穩(wěn),說明設(shè)備數(shù)量對本文方法的影響較小。圖6表明,在不同測站分布條件下,本文方法的估計精度高于其他3種方法,選優(yōu)法依賴于設(shè)備測量誤差,而測量誤差大小與目標距離相關(guān),因此對測站分布較為敏感。綜上所述,本文方法具有較強的伸縮性能。

圖5 本文方法、理想選優(yōu)法、UKF方法、SR方法在不同設(shè)備數(shù)量條件下的估計精度Fig.5 Comparison of MD with other methods under different N values

圖6 本文方法、理想選優(yōu)法、UKF方法、SR方法在不同測站分布條件下的估計精度Fig.6 Comparison of MD with other methods under different distribution

④穩(wěn)定性測試。該項測試用于評價本文方法估計數(shù)據(jù)的穩(wěn)定性,在給定設(shè)備數(shù)量N=20和離群點數(shù)量Z=3的條件下,統(tǒng)計本文方法、理想選優(yōu)法、UKF方法、SR方法的誤差均方差如表5所示。表5的計算結(jié)果表明,在相同的實驗條件下,本文方法的估計穩(wěn)定性優(yōu)于UKF方法和SR方法,選優(yōu)法必須完全理想,其估計穩(wěn)定性才和本文方法相當(dāng)。估計穩(wěn)定性是彈道估計平滑程度的重要指標,該實驗表明本文方法估計的目標飛行軌跡相比UKF方法和SR方法更加平滑,抖動較小。

表5 本文方法和其他三種方法的誤差均方差對比Table 5 Comparison of MSE with other methods

3 結(jié)束語

目前各種火箭彈道實時融合方法難以在估計精度、估計穩(wěn)定性、抗干擾性能、伸縮性能、泛化性能、時間消耗這幾個方面均能夠滿足工程應(yīng)用。針對這些局限性,本文提出了一種基于高斯KDE函數(shù)最大值點跟蹤的彈道實時融合方法。通過理論分析及實驗證明,本文跟蹤方法是有效的,與UKF方法、SR方法、選優(yōu)法相比,具有較高的估計精度、較強的抗干擾性、伸縮性能和穩(wěn)定性,其時間消耗完全滿足工程應(yīng)用。本文方法模型簡潔,涉及參數(shù)少、易于調(diào)優(yōu),且不依賴于理論彈道方程、歷史數(shù)據(jù)等任何先驗知識,具有較強的泛化性能。下一步將對本文方法在彈道事后數(shù)據(jù)處理問題中的應(yīng)用進行研究。