王奕波,刘雪峰,李 玫

(中国人民解放军63801部队,四川 西昌 615042)

由于系统误差和随机误差的存在[1],航天测量任务中任何设备输出的目标测量数据都不可能完全准确,而获取目标的真实飞行弹道,对于实时引导数据生成[2]、设备测量性能分析[3]和目标飞行控制[4]等工作尤为重要,因此一般采用相关的数据处理方法,对不同设备的测量数据进行融合,从而实现目标真实弹道的估计。目前的弹道实时融合方法[5]主要有3种：基于选优的方法[6]、样条约束(spline restraint,SR)法[7]、无迹卡尔曼滤波(unscented Kalman filter,UKF)法[8]。选优法的基本思路是采用如神经网络[9]等方法对各设备测量数据的可靠性进行量化,选取可靠性最高的数据,该方法存在的问题是其估计结果总是某个设备的测量数据,因此估计精度对设备本身的测量精度、测站分布依赖性较强,况且可靠性量化不可能完全正确;SR方法利用一段时间内所有设备的多帧测量数据,采用样条函数对弹道位置进行估计,但由于该方法本质上是对目标真实位置在多帧测量数据下的最小二乘(least squares,LS)估计[10],因此对异常测量数据的抗干扰性差;UKF方法对于高斯噪声过程[11]具有较强的滤波性能,但依赖理论弹道方程[12],其估计精度易受目标真实飞行轨迹偏离理论弹道的影响,并且其估计精度对初值选取较为敏感。

本文基于测量误差主要为高斯噪声这一基本依据[13],证明了以服从高斯分布的有限样本构造的高斯核密度估计(kernel density estimation,KDE)函数[14]的最大值点为产生该组样本高斯分布期望的无偏估计量这一结论,提出了一种基于高斯KDE函数最大值点跟踪的目标弹道实时融合方法,从而将弹道融合问题转化为求该函数最大值点的问题,并通过理论分析与对比实验证明了本文方法的优势。

1 基于高斯核密度估计函数最大值点跟踪的方法

1.1 基于核函数最大值点的估计方法

(1)

证毕。

(2)

1.2 核函数最大值点跟踪

1.2.1 跟踪方法描述

求取ρ(X)最大值点的直接思路是先求出ρ(X)的全部驻点,函数值最大的驻点则为最大值点。但ρ(X)的驻点约束方程ρ(X)=0为超越方程,无法直接求得精确解,因此只能采用一些数值分析方法求得近似解。传统方法是采用迭代的思路,即将每组测量数据作为初值进行迭代运算,取核函数值最大的迭代结果作为最终估计结果。该方法对样本点有一定的聚类[15]能力,能有效求出核函数值最大的位置,但存在耗时较长的问题：一方面是由于该方法采用相邻两步迭代结果的欧式距离阈值作为迭代终止判据,对于不同的测量数据分布,其迭代步数是不定的,较小的阈值能够保证收敛精度,但以牺牲时间成本为代价;另一方面由于该方法需要对每个都进行迭代运算,其时间消耗与设备数量N相关性较大。本文将在2.2节的实验中进一步说明该问题。

针对上述问题,本文提出了一种核函数最大值点跟踪的方法：首先根据每个设备当前帧与前一帧的测量数据求出若干平移矢量;再以每个平移矢量为基准对上一帧的估计结果进行平移采样,从而得到若干采样结果;取当前帧核函数值最大的若干采样结果作为核函数最大值点的近邻,采用一阶泰勒估计的方法,计算出若干核函数最大值点的估计位置;最后取核函数值最大的位置,作为当前帧的估计结果。具体步骤如下。

(3)

(4)

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

1.2.2 跟踪有效性分析

(5)

(6)

(7)

(8)

表1 PQ的计算结果Table 1 Calculation results of PQ

1.2.3 时间复杂度分析

下面对本文跟踪方法的耗时情况进行分析,设求取一次由N组测量数据构造的核函数值的平均耗时为t(N),本文方法耗时主要在步骤③、④、⑤中求核函数值的部分,步骤③的耗时为Nt(N);步骤④需计算S个点的一阶和二阶方向导数,共9个方向,其耗时为9St(N);步骤⑤的耗时为St(N)。因此,本文方法总的理论耗时为(N+10S)t(N)。而经典迭代法的耗时为TNt(N),T为平均迭代步数,最后选取最大核函数值的迭代结果时需求取全部迭代结果的核函数值,耗时为Nt(N),因此总的耗时为N(T+1)t(N)。由于T≥1,S为常数,当N足够大时,本文方法的耗时相对较低,且随N的增大,本文方法相对经典迭代法节省的时间成本越多。

需要说明的是,在飞行初始段,参与测量的设备数量较少,可以运用经典迭代法,不会导致较高的时间消耗;但一段时间后,当设备数量达到一定程度(由上文分析,在工程应用中可设定当设备数量达到5个时),采用本文跟踪方法,可以同时保证估计精度和实时性,将在第2.2节的实验中进一步证明。

2 实验分析

2.1 条件与参数说明

本文实验数据的生成方式如下：

①目标真实飞行轨迹。选用3条航天试验任务的历史弹道(每条弹道总计15 920帧)反演目标的真实飞行轨迹。

②测站分布。假设发射点的坐标为(0 m,0 m,0 m),根据我国领土的实际范围,在以发射点坐标为基准,从(-2×106m,-2×106m,-2×103m)到(2×106m,2×106m,2×103m)的范围内按均匀分布随机产生N个坐标模拟设备分布情况。

④异常测量数据。由于异常测量数据的统计特性难以预知,以目标真实位置为中心、从(0 km,-180°,-90°)到(10 000 km,180°,90°)的范围内,按均匀分布随机生成Z=N-M个离群点。

⑤硬件环境。计算机平台内存1 GB、CPU主频3.0 GHz(双核心)。

⑥其他实验条件。下述所有实验统计的各项指标(误差均值、帧平均耗时等)均为对3条弹道分别运行100次实验(总计4 776 000帧)的平均结果;为保证不同方法误差比较的同一性,帧误差、误差均值、误差均方差基于发射点进行计算;在实验中任意时刻的测量协方差计算方法采用经典的变量差分法[16];由于存在很多种选优法,下述实验中假定采用的选优法是完全理想的,即每一帧都选取测量误差最小的数据,在此基础上与本文方法进行比较,以说明本文方法的高精度优势。

2.2 跟踪参数S的选取及有效性验证

本文跟踪方法涉及的参数S影响耗时,如果S设定得过高则耗时较长,过低则导致近邻丢失,从而降低跟踪精度。为对跟踪精度进行评价,本文计算不同N和S条件下跟踪结果与目标真实位置的平均偏差(设置离群点数量Z=0),结果如表2所示。表中,R为距离,A为方位角,E为俯仰角。

表2 不同N和S条件下的跟踪结果与真实位置平均偏差Table 2 MD between the tracking results and the pinpoint under different N and S values

表2的计算结果表明,当S≥6时估计精度基本不随S的增大而提高,此时估计误差是较低的(在第2.3节的对比实验中进一步证明)。根据1.2.2节的分析,只要S足够大就能保证至少存在一个平移采样结果为KDE函数最大值点的近邻,由于近邻本身是从核函数值最大的S个平移采样结果中选取的,再增加S的值只会将核函数值较小的采样结果包含进来,而按照KDE函数的分布特性,核函数值越小与核函数最大值点的位置就越远,因此再增加S的值不会显著提高估计精度。而根据第1.2.3节的分析,S对耗时的影响较大,因此在工程应用中不宜取较大的值。当N=20和S=6的条件下帧平均耗时约为3.367 m/s(约297 Hz/s,比目前航天测量任务的数据融合频率高一个数量级)。较小的S即可保证必存在一个采样结果为最大值点的近邻,跟踪精度就可以保证,而S的取值不必随N的增加而增大,因此同时保证了实时性。

为进一步证明本文方法的有效性,在N取值不同的条件下与经典迭代法进行比较,结果如表3所示。实验中设定迭代收敛条件为相邻两步迭代结果的欧式距离小于1 m,该条件下迭代法的距离收敛精度近似为1 m,角度收敛精度近似为0.001°。

表3表明,在不同的设备数量条件下,本文方法与经典迭代法的估计精度近似相等,说明本文跟踪方法能够有效估计核函数最大值点。另外,本文对这2种方法在不同设备数量条件下的帧平均耗时情况进行了统计,如图1。当N较小时,2种方法的帧平均耗时均较低,但随N的增大,由于经典迭代法需将每组测量数据作为迭代初值进行运算,其时间消耗与N的相关性较大;随N的增大,本文方法与经典迭代法的帧平均耗时之差也越大,说明随设备数量增多,本文方法相对经典迭代法节省的时间成本越多。

表3 与经典迭代法的估计精度比较Table 3 Comparison of MD with the classic iterative method

图1 与经典迭代法的时间消耗情况比较Fig.1 Comparison of the time-consuming test with the classic iterative method

2.3 对比实验

①抗干扰性测试。该项测试用于评价离群点对本文方法估计精度的影响。在给定N=20的条件下,设定不同的离群点数量Z,统计本文方法、选优法、UKF方法(参数选取σ点数量为3,α=0.5,β=2,λ=0)、SR方法的误差均值,结果如图2所示。

图2 抗干扰性实验的对比结果Fig.2 Comparison of the noise-tolerate capability test with other methods

由图2可知,当不存在离群点时,4种方法均能保持较高的估计精度,而一旦出现离群点,SR方法的估计精度迅速下降,原因是该方法将当前帧及历史多帧的测量数据采用LS方法进行估计,而LS方法的数学本质是求解一个与所有测量数据距离之和最小的点作为估计结果,因此SR方法易受离群点位置及数量的影响,且离群点分布越远,对估计结果的影响越大,这与本文方法恰好相反。选优法只有在完全理想的情况下,即选优完全正确的情况下,其估计精度才和本文方法相当。UKF方法具有一定的抗干扰性能,这是由于该方法利用了历史先验信息,即理论弹道方程,因此具有一定的抗噪能力,但由于该方法依赖于理论弹道方程,如果目标飞行偏离理论弹道,UKF方法的估计精度会迅速降低,如图3所示,以存在俯仰偏差为例。而本文方法不依赖理论弹道方程,不受目标真实飞行是否偏离理论弹道的影响。

图3 UKF方法对理论弹道方程的敏感度测试Fig.3 Initial value sensitivity test for UKF

②初值敏感度测试。本文方法和UKF方法同属于跟踪方法,该项测试用于评价第一帧数据对这2种方法跟踪精度的影响。在实际航天发射任务中的初始段,只有少量设备(如遥测雷达)参与测量,因此假设第一帧数据只由一台设备产生,在第一帧数据与目标位置真值存在不同偏差的条件下,比较本文方法和UKF方法的跟踪性能(设置设备数量N=20,离群点数量Z=0),如图4所示。

图4表明,当初值存在一定偏度时(与目标真实位置距离偏差为10 000 m,方位角、俯仰角偏差均为0.2°),本文方法和UKF方法均具有一定的收敛性能。而在相同的初值条件下,本文方法的收敛精度更高,如表4所示。表中,ξ为不同的初值偏度系数,所模拟的距离、方位、俯仰初值偏度步进分别为4 000 m,0.08°,0.08°,即所模拟的距离、方位、俯仰初值分别为对应真值加上4 000ξ,0.08ξ和0.08ξ。本文方法和UKF方法的共同点是均采用了多点采样的方式,UKF方法对采样点进行加权求和,利用当前帧实测结果对加权求和结果进行修正,因此UKF方法的收敛精度与初值偏度具有较大的相关性。而本文跟踪方法是在多个采样点中选取一个目标真实位置的近邻进行估计,根据第1.2.2节的分析,在若干采样点中至少存在1个近邻的概率与前一帧的估计结果并无较大关系,因此本文方法具有较高的收敛精度。

表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

③伸缩性能测试。该项测试用于评价本文方法在不同设备数量和不同测站分布条件下的估计精度,设置离群点数量Z=0,统计不同设备数量(如图5所示)和不同测站分布(如图6所示,设置设备数量N=20)条件下本文方法、理想选优法、UKF方法、SR方法的估计精度。

图5表明,一方面本文方法在不同设备数量条件下的估计精度总是高于UKF方法和SR方法,而选优法在完全理想的情况下,估计精度才和本文方法相当;另一方面,本文方法的数据曲线较为平稳,说明设备数量对本文方法的影响较小。图6表明,在不同测站分布条件下,本文方法的估计精度高于其他3种方法,选优法依赖于设备测量误差,而测量误差大小与目标距离相关,因此对测站分布较为敏感。综上所述,本文方法具有较强的伸缩性能。

图5 本文方法、理想选优法、UKF方法、SR方法在不同设备数量条件下的估计精度Fig.5 Comparison of MD with other methods under different N values

图6 本文方法、理想选优法、UKF方法、SR方法在不同测站分布条件下的估计精度Fig.6 Comparison of MD with other methods under different distribution

④稳定性测试。该项测试用于评价本文方法估计数据的稳定性,在给定设备数量N=20和离群点数量Z=3的条件下,统计本文方法、理想选优法、UKF方法、SR方法的误差均方差如表5所示。表5的计算结果表明,在相同的实验条件下,本文方法的估计稳定性优于UKF方法和SR方法,选优法必须完全理想,其估计稳定性才和本文方法相当。估计稳定性是弹道估计平滑程度的重要指标,该实验表明本文方法估计的目标飞行轨迹相比UKF方法和SR方法更加平滑,抖动较小。

表5 本文方法和其他三种方法的误差均方差对比Table 5 Comparison of MSE with other methods

3 结束语

目前各种火箭弹道实时融合方法难以在估计精度、估计稳定性、抗干扰性能、伸缩性能、泛化性能、时间消耗这几个方面均能够满足工程应用。针对这些局限性,本文提出了一种基于高斯KDE函数最大值点跟踪的弹道实时融合方法。通过理论分析及实验证明,本文跟踪方法是有效的,与UKF方法、SR方法、选优法相比,具有较高的估计精度、较强的抗干扰性、伸缩性能和稳定性,其时间消耗完全满足工程应用。本文方法模型简洁,涉及参数少、易于调优,且不依赖于理论弹道方程、历史数据等任何先验知识,具有较强的泛化性能。下一步将对本文方法在弹道事后数据处理问题中的应用进行研究。