张雪松,吴 楠,王 锋,李东泽

(1. 战略支援部队信息工程大学,河南 郑州 450001;2. 中国人民解放军32035部队,陕西 西安 710000;3. 军事科学院 国防科技创新研究院,北京 100071)

各类新型导弹、助推滑翔飞行器等兵器具有机动能力强、飞行空域大、打击精度高等优势,被广泛运用于实际作战,并得到国内外学者的广泛关注和研究[1-2]。各类飞行器通过搭载的传感器获取拦截器的测量信息进行灵活控制和机动调整,不仅丰富了探测手段、拓展了探测范围,还增强了突防能力,是机动突防和反拦截的发展方向[3-5]。研究人员采用单飞行器搭载红外传感器对拦截弹进行探测,并在弱可观条件下完成对拦截弹主动段的弹道估计,是飞行器信息处理技术研究的关键环节之一。

单飞行器搭载红外探测器探测拦截弹主动段时,同一时刻只能获得两个角度测量值,对拦截弹定位来说属于不完备观测[6]。主动段运动规律复杂,在弱可观测条件的限制下,运动模型复杂会导致模型不可辨识,运动模型简单会带来较大辨识误差[7-8]。解决此类问题有两种思路,一是融合其他先验信息以弥补观测信息的不完备,二是构建主动段的模板信息[9]。

目前,基于模板信息是解决问题的重要方法,主要包括标准弹道模板、推力加速度模板和动力特征模板。基于标准弹道模板或推力加速度模板等导弹的先验信息来提高弹道估计精度,主要受限于模板库的准确性和完备程度,而且运算时会出现不收敛的问题,在实际应用中较为困难[10-11]。动力特征模板基于同型号导弹主动段推力加速度曲线相同的事实,是同型号导弹的本质特征,体现了导弹的运动特性。在动力特征模板已知的条件下,导弹推力随时间的函数是确定的,未知的是所采取的飞行程序角模型[12]。因此,对拦截弹的弹道估计,主要解决的问题是对不同飞行程序角模型的辨识,以及拦截弹位置、速度的准确估计。

因此,在已有文献的研究基础上,本文针对单飞行器对拦截弹主动段弹道估计问题,通过引入拦截弹动力特征模板,将飞行程序角模型设计参数与发射方位角增广为目标状态分量,基于程序角模型和动力特征模板构建非线性MMAE模型集,研究采用MMAE-UKF算法对拦截弹的程序角模型和设计参数进行辨识,同时对主动段运动状态进行估计。最后,通过仿真实验验证了算法的有效性。

1 主动段弹道估计模型

1.1 动力特征模板

对于特定型号的拦截弹,部署区域已知,推进剂质量、发动机性能、燃烧效率和助推时长相对固定,因此推力加速度有着固定的变化规律[13]。也就是说当某型号拦截弹总体参数(包括质量、推力等特性参数)已知时,推力随时间的函数是确定的,未知的是不同弹道采用的程序角φ(t)和发射方位角A0。同时,在实际中为了减少控制设计的复杂度,程序角时间节点也是固定的。图1为某型单级助推拦截弹从不同发射点飞向预测拦截点的弹道和推力加速度曲线,可见对于同种类型的导弹,即使弹道不同,其推力加速度曲线基本相同,动力特征模板正是基于此提出的,这是同种类型导弹的本质特点,具有更强的实用性。

图1 同种类型拦截弹弹道及其推力加速度比较Fig.1 Trajectory comparison and thrust acceleration comparison of same type interceptor

1.2 主动段动力学模型

在拦截弹的助推段,总加速度a由以下几项组成,包括推力加速度aT、气动力加速度aA、引力加速度ag、离心力加速度acf和科氏力加速度aco,即

a=aT+aA+ag+acf+aco

(1)

(2)

(3)

式中,R[]为初等变换矩阵,φ、ψ、γ分别为俯仰角、偏航角和滚转角,通常偏航角和滚转角为0。

(4)

式中,A为发射方位角,φ、λ为发射点纬度、经度,可根据作战区域先验信息获得。

(5)

式中,ν、σ、θ分别为倾侧角、航迹偏航角、速度倾角,通常倾侧角、航迹偏航角为0。

(6)

对式(6)进行离散化处理,可以得到主动段运动模型:

Xk+1=ΦkXk+Ukak

(7)

1.3 程序角模型

对于一级助推的拦截弹来说,助推段主要包括垂直上升段(t0

1)直线型程序角

导弹在程序转弯段是直接给出直线斜率形式变化的程序角。程序角变化规律为

(8)

2)攻角控制型程序角

拦截弹在程序转弯段采用攻角控制程序角变化,可分为指数型攻角和三角函数型攻角。

(9)

式中,θ(t)为速度倾角;α(t)为攻角。

当α(t)为指数型攻角时

α(t)=-4αmz(1-z)

(10)

式中,z=e-0.28(t-t1),αm为攻角绝对值的最大值。

当α(t)为三角函数型攻角

α(t)=-αmsin2(f(t))

(11)

式中,

(12)

采用该种程序角的设计参数为αm。

3)抛物线型程序角

工程上设计出一种实用的抛物线形式的程序角,变化规律为

(13)

其中,

(14)

式中,C为选择的某一常数,是程序角的设计参数。

1.4 飞行器探测模型

利用红外探测器进行探测时,观测方程为角度观测量(A,E)与目标状态Xk间的关系。

Zk(A,E)=h(Xk)+wk

(15)

(16)

由拦截弹位置矢量可以计算红外探测器的两个角度测量值:定义方位角A为位置矢量在xoy平面内的投影与oy轴的夹角,仰角E为位置矢量与xoy平面的夹角。则位置矢量与红外探测器角度测量值的关系为

(17)

式(16)～(17)构成了观测方程(12)的解析表达式。

2 主动段弹道估计多模型滤波器设计

MMAE算法可以对参数不确定模型的状态进行估计,具有模型匹配程度高、参数辨识准确、易于工程实现的优点。在本文中,算法的核心步骤有以下3步:一是根据拦截弹的动力特征模型和程序角设计参数构建模型集;二是基于各个模型分别对拦截弹运动状态、程序角设计参数和发射方位角进行滤波初始化;三是基于各个模型分别设计一组单元滤波器,各个滤波器并行计算,获得各自状态估计值、程序角设计参数和发射方位角。通过各个单元滤波器的测量残差计算得到相应滤波器的权值(预测概率),能反映各个模型与实际系统的一致程度,从而实现拦截弹程序角模型确定及弹道参数的辨识,同时运动状态估计值是各个并行滤波器状态估计值的加权和。

2.1 MMAE模型集构建

在本文中,假设拦截弹的程序角模型为1.3节中的4种模型之一,并根据每种程序角的设计参数,与拦截弹的动力特征模型组合,构成MMAE算法的所需模型集β={βj|j=1,2,3,4},βj表示模型j。

2.2 滤波初始化方法

1)运动状态初始化

由于红外探测器仅能测得视线矢量的两个角度,无法测量距离,无法对目标形成三维立体定位,需要通过引入其他约束条件进行初始化定位。如图2所示,假设拦截弹在P处被飞行器S搭载的红外探测器首次探测,此时拦截弹距离地面很近,可认为位于地球椭球表面,视线射线与地球椭球表面的交点便是目标位置,可以利用红外探测器首次测量数据的两个角度测量值对目标进行定位。

图2 飞行器与拦截弹的空间几何关系Fig.2 Space geometric relationship between aircraft and interceptor

在飞行器视场坐标系中,设三角形OSP中∠OSP=S,则

(18)

其中,

(19)

(20)

(21)

式中,A为观测方位角,E为观测俯仰角。

根据余弦定理有

d2=f2+b2-2bfcos(S)

(22)

整理可得

f2+(-2bcos(S))f+(b2-d2)=0

(23)

(24)

转换到地固坐标系中:

(25)

(26)

由此,式(18)～式(26)便构成了拦截弹初始位置的定位方程,可以求得初始位置作为滤波初始状态估计值。根据拦截弹初始位置的定位方程以及测量数据的均值和协方差,采用无迹变换的方法可以得到滤波初始协方差矩阵。

2)弹道参数初始化

t=R/(VT+VM)

(27)

2.3 MMAE-UKF算法流程

MMAE-UKF算法主要分为3个步骤,即并行UKF滤波、模型概率概新和状态及弹道参数估计输出。

1)并行UKF滤波

由于各模型的滤波算法相同,均为UKF滤波算法,因此这里仅给出第j个模型βj的UKF滤波算法,设k时刻的状态为

(28)

具体算法流程如下:

a.状态预测

(29)

b.协方差预测

(30)

c.采用修正比例采样规则构造Sigma采样点χk+1/k(i)和相应权重Wi

(31)

d.基于观测方程的非线性传播

(32)

e.计算观测值的预估值和协方差矩阵

(33)

f.计算新息矢量和新息协方差矩阵

(34)

g.状态估计与协方差更新

(35)

2) 模型概率更新

采用Gaussian density函数计算在k+1时刻模型j为匹配模型的似然概率为

(36)

根据似然概率进行模型匹配概率的更新,各模型的预测概率为

(37)

式中,C为归一化系数。

(38)

由上述公式可知,各个模型的预测概率之和为1,且由于各模型之间的程序角设计参数的差异,导致单元滤波器对模型匹配的程度不同,因而滤波得到的预测概率存在差异,即模型集中能够与实际弹道采用的程序角模型相匹配的,预测概率快速收敛至1,而不能够匹配的,预测概率收敛于0,据此可以判断实际弹道所采用的程序角模型,对应单元滤波器输出的程序角设计参数、发射方位角估计值即为所求。

3)运动状态及弹道参数估计

根据多个模型的估计输出,则估计的总输出为

(39)

总的协方差矩阵为

(40)

由模型概率更新结果确定实际拦截弹所采用的程序角设计模型,则对应滤波器输出的程序角设计参数、发射方位角即为弹道参数辨识结果。

图3 算法的整体流程Fig.3 The overall flow of the algorithm

3 仿真实验与结果分析

为验证算法有效性,本文构建三维仿真场景,一方面验证对于不同程序角模型的辨识以及对增广参数估计的能力,另一方面验证算法对于运动状态的滤波性能,并与传统的增广状态的滤波算法进行了对比。

3.1 仿真场景

根据飞行器动力学方程和数值积分算法,优化一条飞行器滑翔弹道。拦截弹动力特征已知,飞行程序角节点为t1=2 s,t2=20 s,即t1时刻前垂直上升,t1～t2进行程序转弯。发射点地理坐标为19.6°E,0.4°N,高程为5 m,发射方位角为-117.65°,分别采用1.3节的程序角模型优化出一条飞向预测拦截点的弹道,4种程序角模型的参数取值如表1所示。

构建三维场景如图4所示,采用不同程序角模型的拦截弹射程-高度曲线如图5所示。飞行器搭载红外探测器对拦截弹进行探测,并假设拦截弹发射后即被飞行器发现,根据飞行器和拦截弹的相对关系,计算出助推段观测数据,并加入测量噪声,生成测量数据,数据率为T=0.1 s。

图5 射程-高度曲线Fig.5 Range-altitude curve

3.2 仿真环境及参数设置

实验平台如下:移动工作站,操作系统为Windows 10,CPU为Intel Core(TM) i7-9850(2.60 GHz),内存为32 GB,显卡为NVIDIA Quadro RTX 3000;编程环境为Matlab R2020b。

本文算法的滤波初始化采用2.2节方法,利用首次测量数据得到运动状态的初始估计值为

X0=[6 008 434.5,2 139 507.3,44 229.417,0,0,0]T

初始协方差设为

P0=diag([1e4,1e4,1e4,1e2,1e2,1e2])

程序角模型设计参数初值为直线型程序角0.05,指数型程序角、三角函数型程序角12.5,抛物线型程序角0.7。首次探测到拦截弹时,飞行器位置为[6 155 618.1 1 780 988.1 -0.000 460 42]T,速度为2 752 m/s,拦截弹初始化位置为[6 008 434.5,2 139 507.3,44 229.417]T,飞行平均速度取750 m/s,经计算发射方位角初值为115.1°。上述状态分量的过程噪声方差设计采用文献[16]的方法。传统的增广状态的滤波算法采用UKF,将时变参数程序角φ(t)与发射方位角A0作为增广状态分量进行滤波,采用同样的滤波初始化方法。观测角度标准差为1×10-4rad。

3.3 仿真结果及分析

1)程序角模型辨识

图6 直线型程序角辨识结果Fig.6 The identified result for flight program of linear form

图7 指数型程序角辨识结果Fig.7 The identified result for flight program of exponential form

图8 三角函数型程序角辨识结果Fig.8 The identified result for flight program of trigonometric form

图9 抛物线型程序角辨识结果Fig.9 The identified result for flight program of parabolic form

从上面的辨识结果可以看出,对于4种不同程序角模型的拦截弹助推段弹道,MMAE-UKF算法可以有效实现对拦截弹程序角模型的辨识。如图6a)所示,由于直线型程序角模型采用俯仰角变化率直接控制俯仰角的变化,与指数型程序角模型、三角函数型程序角模型及抛物线型程序角模型差异较大,能很快地被排除或者确认,因此在拦截弹577.4 s进入程序转弯后匹配概率便不断上升,在583 s前便很快被辨识匹配。对于指数型程序角模型、三角函数型程序角模型来说,两者属于攻角控制型程序角模型,对于程序俯仰角控制具有一定的相似性,如图7a)、图8a)所示,在辨识过程中两者在程序转弯时初始时会彼此互相影响,且同时受到抛物线型程序角模型的影响,但两者分别能在586.5 s和588 s前辨识出所匹配的程序角模型。对于抛物线型程序角来说,如图9a)所示,主要受两种攻角控制型程序角模型影响,且受三角函数型程序角模型影响较大,因此匹配概率在585 s才开始缓慢上升,约在590 s才完全匹配。

从4种程序角模型对应的拦截弹攻角来分析,可以更为显著地说明产生这种结果的原因,如图10所示。可以看出,在577.4 s进入程序转弯后,直线型程序角模型与其余模型的攻角差异明显,这也导致弹道有明显差异,因而能够较为快速地辨识匹配。同理,两种攻角控制型程序角模型虽然具有相同的最大负攻角参数,但两者到达最大负攻角的时间是不同的,因此虽然两者在辨识的起始阶段互有影响,但仍然能够明显辨识匹配。而抛物线型程序角模型的攻角在起始阶段介于指数型程序角模型与三角函数型程序角模型之间,且到达最大负攻角的时间与三角函数型程序角模型相近,因此在辨识起始时受三角函数型程序角模型影响较大,匹配概率上升缓慢。

图10 四种程序角模型对应的攻角Fig.10 Angle of attack corresponding to four flight program

对于设计参数估计来说,如图6b)、图7b)、图8b)和图9b)所示,指数型程序角、三角函数型程序角以及抛物线型程序角辨识结果较好,三者均能在585 s前实现参数估计的收敛,且参数估计值稳定收敛于真值附近。而直线型俯仰角虽然能够辨识出设计参数,但在585 s后参数才收敛,且参数估计值不够平稳,这是因为直线型设计参数值较小,在滤波过程中受噪声影响明显。

传统增广状态滤波算法直接以程序角这一时变参数作为增广状态分量进行滤波,这里以直线型程序角模型对应的弹道为例进行分析,其余几种程序角模型的弹道类似,辨识结果如图11所示。图中黑色实线为利用本文算法辨识出的设计参数重建后得到的程序角曲线,黑色虚线为传统增广状态滤波算法的辨识结果,红色虚线为程序角的真实值。可以看出,与真实程序角曲线相比,传统增广状态滤波算法的辨识结果精度较低,而本文算法采用对设计参数这一固定值进行辨识,能得到较为准确的设计参数估计值,重建后能得到比较准确的程序角变化曲线。

图11 程序角辨识结果对比Fig.11 The identified results comparison of flight program angle

2)发射方位角辨识

由于拦截弹的发射点和预测拦截点都是固定的,采用4种程序角模型的弹道,其发射方位角是相同的。这里仅对采用直线型程序角模型的弹道发射方位角进行辨识分析,其余几种程序角模型的弹道类似。如图12所示,本文算法和传统增广状态滤波算法对发射方位角的辨识估计结果总体上比较准确,但本文算法收敛速度更快,辨识更为准确。本文算法在583 s时滤波估计结果开始平稳且收敛于真值,但在初始阶段发射方位角估计值产生剧烈震荡,这是因为拦截弹在垂直上升段和程序转弯起始的一个阶段射向不明显,因此辨识估计结果较差。

图12 发射方位角辨识结果Fig.12 The identified results for launching azimuth

3)运动状态估计

对于运动状态估计,这里以直线型程序角模型对应的弹道进行分析,其余几种程序角模型的弹道类似。采用均方根误差RMSE表示跟踪过程中的误差变化。定义如下:

(41)

式中,xi(k)和xi(k|k)分别表示第i次仿真中第k时刻真值和估计值,N为蒙特卡洛仿真次数。

以x轴方向状态估计均方根误差为例,其他方向基本相同,进行500次蒙特卡罗仿真,结果如图13和图14所示。实线描述本文算法的弹道估计结果,虚线描述基于传统增广状态滤波算法的弹道估计结果。

图13 X方向位置均方根误差Fig.13 The RMSE of position in X direction

图14 X方向速度均方根误差Fig.14 The RMSE of velocity in X direction

从图13和图14可以看出,与传统增广状态滤波算法相比,本文算法估计精度更高。本文算法滤波稳定时的位置均方根误差约为5 m,而传统增广状态滤波算法的位置均方根误差约为10 m,本文算法的位置均方根误差比传统增广状态滤波算法减小50%;本文算法速度估计的均方根误差约为1.2 m/s,而传统增广状态滤波算法的速度均方根误差约为1.8 m/s,本文算法速度估计的均方根误差约比传统增广状态滤波算法减小30%。此外,本文算法在滤波时能更快到达稳定状态。

4 结束语

本文以单飞行器搭载红外探测器对拦截弹主动段探测跟踪为背景,基于动力特征模板和程序角模型构建非线性MMAE模型集,提出一种可以对拦截弹主动段程序角模型及设计参数进行辨识,同时对运动状态进行估计的MMAE-UKF算法。仿真结果表明,本文算法能够快速准确对拦截弹主动段程序角模型及设计参数进行辨识。同时,与传统的增广状态的滤波算法相比,运动状态估计误差显著减小,滤波稳定速度大幅度提高。后续,可以在此基础上进行弹道外推,对拦截弹的预测拦截点进行初步估计,评估拦截弹中末制导所需修正误差的量级大小,对飞行器规避拦截具有理论研究意义。