基于核函數(shù)最大值點(diǎn)跟蹤的彈道實(shí)時(shí)融合方法

王奕波,劉雪峰,李 玫

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

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

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

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

1.1 基于核函數(shù)最大值點(diǎn)的估計(jì)方法

(1)

證畢。

(2)

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

1.2.1 跟蹤方法描述

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

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

(3)

(4)

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

1.2.2 跟蹤有效性分析

(5)

(6)

(7)

(8)

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

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

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

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

2 實(shí)驗(yàn)分析

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

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

①目標(biāo)真實(shí)飛行軌跡。選用3條航天試驗(yàn)任務(wù)的歷史彈道(每條彈道總計(jì)15 920幀)反演目標(biāo)的真實(shí)飛行軌跡。

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

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

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

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

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

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

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

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

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

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

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

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

2.3 對(duì)比實(shí)驗(yàn)

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)束語

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