[发明专利]一种临近空间高超声速滑翔弹头跳跃弹道预测方法

权利要求书

1.一种临近空间高超声速滑翔弹头跳跃弹道预测方法，其特征在于它按以下步骤实现：

一、高超声速滑翔弹头的弹道方程的建立；

建立临近空间高超声速滑翔弹头的弹道方程，包括气动力计算方程、加速度计算方程、速度和位置计算方程、速率和弹道角计算方程，以及姿态角计算方程；

综合计算方程，得到高超声速滑翔弹头质心运动速度和位置方程如式(1)～(3)所示：

x·I=vxIv·xI=-μxIrI3+aaxI---(1)]]>

y·I=vyIv·yI=-μyIrI3+aayI---(2)]]>

z·I=v·zIv·zI=-μzIrI3+aazI---(3)]]>

其中，x_I、y_I和z_I代表高超声速滑翔弹头在地心惯性坐标系中的位置，代表高超声速滑翔弹头的位置到地球中心的距离,和则代表高超声速滑翔弹头在地心惯性坐标系中的速度，a_axI、a_ayI和a_azI代表气动力产生的高超声速滑翔弹头机动加速度在地心惯性坐标系中的分量，μ为引力参数；

二、高超声速滑翔弹头的跟踪滤波器设计；

设计高超声速滑翔弹头的运动轨迹跟踪卡尔曼滤波器，实时地跟踪并估计高超声速滑翔弹头的位置、速度和加速度；

其中，所述高超声速滑翔弹头的运动轨迹跟踪卡尔曼滤波器的状态一步预报方程为

X‾i(k+1)=ΦiX^i(k)Pi(k+1/k)=ΦiPi(k)ΦiT+Qi,i=x,y,z---(4)]]>

上式中，代表滤波器i第k步的状态估计值，代表滤波器i从第k步到第k+1步的状态预报值，P_i(k)代表滤波器i第k步的状态估计误差协方差矩阵，P_i(k+1/k)代表滤波器i从第k步到第k+1步的状态预报误差协方差矩阵，Q_i为滤波器i模型预报误差协方差矩阵,

Φi=exp(FiΔt)=1Δt1/λi2(e-λit+λiΔt-1)011/λi1-e-λiΔt)00e-λiΔt---(5)]]>

所述运动轨迹跟踪卡尔曼滤波器的测量修正方程为

Ki(k+1)=Pi(k+1/k)HiT[HiPi(k+1/k)HiT+Ri]-1X^i(k+1)=X‾i(k+1)+Ki(k+1)[ηi(k+1)-HiX‾i(k+1)]Pi(k+1)=[I-Ki(k+1)Hi]Pi(k+1/k),i=x,y,z---(6)]]>

其中，上标^T代表矩阵转置运算，上标^-1代表矩阵求逆运算，滤波器i的测量矩阵

H_i＝[1 0 0],i＝x,y,z(7)

R_i为滤波器i测量误差协方差矩阵，I代表单位矩阵，K_i(k+1)代表滤波器i第k+1步的滤波增益矩阵，η_i(k+1)代表滤波器i第k+1步的测量信息，P_i(k+1)代表滤波器i第k+1步的状态估计误差协方差矩阵；

三、高超声速滑翔弹头的弹道预测；

(三一)根据步骤二中运动轨迹跟踪卡尔曼滤波器跟踪结束kT时刻的高超声速滑翔弹头的加速度，其中，k为大于0的整数，T为计算周期，一般可取T＝0.01s，应用步骤一中的临近空间高超声速滑翔弹头的加速度计算方程求出跟踪结束kT时刻临近空间高超声速滑翔弹头的气动力产生的加速度，并根据跟踪结束kT时刻高超声速滑翔弹头的速度，应用步骤一中的临近空间高超声速滑翔弹头的弹道角计算方程计算其弹道角；

(三二)根据计算出来的kT时刻临近空间高超声速滑翔弹头的气动力产生的加速度和弹道角,结合步骤一中的临近空间高超声速滑翔弹头的气动力计算方程和姿态角计算方程估算跟踪结束kT时刻高超声速滑翔弹头飞行时的攻角和姿态角；

(三三)设在随后的kT到(k+1)T，(k+2)T,...,(k+N)T时间内，N为大于2的正整数，临近空间高超声速滑翔弹头进行等攻角和等滚转角飞行，则纵向平面内呈现出上下跳跃式的飞行弹道；同时，高超声速滑翔弹头采用BTT控制方式，侧滑角β等于零，通过控制滚转角在其体系侧向输出一个机动加速度进行转弯，即滚转角指令γ_c保持常值；

(三四)根据等攻角、零侧滑角和等滚转角飞行条件，计算kT时刻气动力产生的加速度、kT时刻的位置x_I、y_I和z_I、kT时刻的速度v_xI、v_yI和v_zI代入式(1)～(3)，经过一步数值积分运算后得到高超声速滑翔弹头在(k+1)T时刻在地心惯性坐标系中的速度和位置；

计算高超声速滑翔弹头在(k+1)T时刻受到的气动力，计算高超声速滑翔弹头在(k+1)T时刻受到的气动力加速度，再把(k+1)T时刻的气动力加速度以及(k+1)T时刻的位置和速度代入式(1)～(3),计算出(k+2)T时刻的位置和速度，以此类推，计算出直到(k+N)T时刻的位置和速度；

步骤一具体为：

(一)计算高超声速滑翔弹头气动力

在滑翔阶段，高超声速滑翔弹头受到的气动力在体系内的投影按照下式计算：

Rx1=-cxqSRy1=cysign(α)qSRz1=0---(8)]]>

其中，α为攻角，

sign(α)=1,α>0-1,α<0]]>

c_y为升力系数，c_x为阻力系数，S为特征参考面积；

q＝ρv²/2(9)

其中，q为动压头，v为高超声速滑翔弹头的速率，ρ为空气密度，它随着飞行高度h变化，可查标准的大气密度表得到；高超声速滑翔弹头飞行高度h用下式计算:

h=xI2+yI2+zI2-Re---(10)]]>

其中，x_I、y_I和z_I代表高超声速滑翔弹头在地心惯性坐标系中的位置，R_e＝6378km为地球平均半径；

设在飞行中侧滑角β＝0，因此气动力侧向力R_z1＝0；

(二)计算高超声速滑翔弹头的速率

高超声速滑翔弹头的速率为

v=vxI2+vyI2+vzI2---(11)]]>

其中，v_xI、v_yI和v_zI代表高超声速滑翔弹头在地心惯性坐标系中的速度；

将高超声速滑翔弹头的速度矢量投影到目标惯性坐标系v_xt0、v_yt0与v_zt0，即

vxt0vyt0vzt0=C0→t0vxIvyIvzI---(12)]]>

其中，

C0→t0=-10001000-1]]>

(三)计算高超声速滑翔弹头的弹道仰角θ和弹道偏角ψ_v：

θ=sin-1(vyt0/v)ψv=-tan-1(vzt0/vxt0)---(13)]]>

高超声速滑翔弹头的姿态角ψ和γ按照下式计算：

其中，为俯仰角，ψ为偏航角，γ为滚转角，γ_c为滚转角指令；

(四)计算高超声速滑翔弹头机动加速度在地心惯性坐标系中的投影

aaxIaayIaazI=C0→t0TCt0→1TRx1/mRy1/mRz1/m---(15)]]>

其中，

m为高超声速滑翔弹头质量；

(五)计算高超声速滑翔弹头的速度和位置

高超声速滑翔弹头质心运动速度和位置方程如式(1)～(3)所示,其中包含了高超声速滑翔弹头受到的重力加速度，即

gxI=-μxIrI3gyI=-μyIrI3gzI=-μzIrI3---(16)]]>

其中，g_xI、g_yI与g_zI代表重力加速度在地心惯性系中的投影；

步骤二具体为：

(一)在探测点惯性坐标系的三个轴上，均用Singer模型来分别描述高超声速滑翔弹头的加速度分量的变化，即

a·x=-λxax+wxa·y=-λyay+wya·z=-λzaz+wz---(17)]]>

其中，λ_x、λ_y和λ_z分别代表在探测点惯性坐标系下o_Fx轴、o_Fy轴和o_Fz轴方向上目标机动时间常数的倒数，w_x、w_y和w_z分别代表探测点惯性坐标系下o_Fx轴、o_Fy轴和o_Fz轴方向上的零均值高斯白噪声；加速度a_x、a_y和a_z中既包含气动力产生的加速度a_axI、a_ayI和a_azI，又包含重力产生的加速度g_xI、g_yI和g_zI；探测点惯性坐标系的o_Fx轴、o_Fy轴和o_Fz轴与地心惯性系的三个轴o_Ix_I、o_Iy_I和o_Iz_I分别平行，气动力产生的加速度在地心惯性系中的投影a_axI、a_ayI和a_azI就是此加速度在探测点惯性系中的投影；

(二)在探测点惯性坐标系下建立起如下三组运动方程

x·=vxv·x=axa·x=-λxax+wx---(18)]]>

y·=vyv·y=aya·y=-λyay+wy---(19)]]>

z·=vzv·z=aza·z=-λzaz+wz---(20)]]>

其中，x、y和z代表高超声速滑翔弹头在探测点惯性坐标系中的位置，和则代表高超声速滑翔弹头速度在探测点惯性坐标系中的投影；因为探测点惯性坐标系的o_Fx轴、o_Fy轴和o_Fz轴与地心惯性系的三个轴o_Ix_I、o_Iy_I和o_Iz_I分别平行，所以

vx=vxIvy=vyIvz=vzI---(21)]]>

而由于探测点惯性坐标系只是相对于地心惯性系探测点惯性坐标系在o_Iy_I轴方向上作了平移，所以有

x=xIy=yI-Rez=zI---(22)]]>

对于系统(18)，定义x轴子系统状态变量X_x＝[x v_x a_x]^T，则得到如下状态方程

X·x=FxXx+Gxwx---(23)]]>

其中，

Fx=01000100-λx,Gx=001---(24)]]>

对式(23)以一定的周期Δt进行离散化后，得到从第k步到第k+1步的状态方程为

X_x(k+1)＝Φ_xX_x(k)+ω_x(k)(25)

其中，

Φx=exp(FxΔt)=1Δt1/λx2(e-λxΔt+λxΔt-1)011/λx(1-e-λxΔt)00e-λxΔt---(26)]]>

ω_x(k)代表x轴子系统状态噪声向量，为零均值高斯白噪声向量；

高超声速滑翔弹头的位置由地面目标跟踪装置测量得到，则测量方程为

η_x＝x+υ_x(27)

其中,υ_x代表x轴子系统测量噪声向量，为零均值高斯白噪声向量；

同理，对于系统(19)，定义y轴子系统状态变量X_y＝[y v_y a_y]^T，则得到如下状态方程

X·y=FyXy+Gywy---(28)]]>

其中，

Fy=01000100-λy,Gy=001---(29)]]>

对式(28)以一定的周期Δt进行离散化后，得到从第k步到第k+1步的状态方程为

X_y(k+1)＝Φ_yX_y(k)+ω_y(k)(30)

其中，

Φy=exp(FyΔt)=1Δt1/λy2(e-λyΔt+λyΔt-1)011/λy(1-e-λyΔt)00e-λyΔt---(31)]]>

ω_y(k)代表y轴子系统状态噪声向量，为零均值高斯白噪声向量；

测量方程为

η_y＝y+υ_y(32)

其中,υ_y代表y轴子系统测量噪声向量，为零均值高斯白噪声向量；

对于系统(20)，定义z轴子系统状态变量X_z＝[z v_z a_z]^T，则得到如下状态方程

X·z=FzXz+Gzwz---(33)]]>

其中，

Fz=01000100-λz,Gz=001---(34)]]>

对式(28)以一定的周期Δt进行离散化后，得到从第k步到第k+1步的状态方程为

X_z(k+1)＝Φ_zX_z(k)+ω_z(k)(35)

其中，

Φz=exp(FzΔt)=1Δt1/λz2(e-λzΔt+λzΔt-1)011/λz(1-e-λzΔt)00e-λzΔt---(36)]]>

ω_z(k)代表z轴子系统状态噪声向量，为零均值高斯白噪声向量；

测量方程为

η_z＝z+υ_z(37)

其中,υ_z代表z轴子系统测量噪声向量，为零均值高斯白噪声向量；

把探测点惯性坐标系三个轴上的子系统，即式(23)和(27)构成的x轴子系统，式(28)和(32)构成的y轴子系统，以及式(33)和(37)构成的z轴子系统，分别设计Kalman滤波器，其预报方程为公式(4),即

X‾i(k+1)=ΦiX^i(k)Pi(k+1/k)=ΦiPi(k)ΦiT+Qi,i=x,y,z]]>

上式中，代表滤波器i第k步的状态估计值，代表滤波器i从第k步到第k+1步的状态预报值，P_i(k)代表滤波器i第k步的状态估计误差协方差矩阵，P_i(k+1/k)代表滤波器i从第k步到第k+1步的状态预报误差协方差矩阵，Q_i为滤波器i模型预报误差协方差矩阵,

滤波器的测量修正方程为公式(6),即

Ki(k+1)=Pi(k+1/k)HiT[HiPi(k+1/k)HiT+Ri]-1X^i(k+1)=X‾i(k+1)+Ki(k+1)[ηi(k+1)-HiX‾i(k+1)]Pi(k+1)=[I-Ki(k+1)Hi]Pi(k+1/k),i=x,y,z]]>

其中，上标^T代表矩阵转置运算，上标^-1代表矩阵求逆运算，滤波器i的测量矩阵为式(7),即

H_i＝[1 0 0],i＝x,y,z

R_i为滤波器i测量误差协方差矩阵，I代表单位矩，K_i(k+1)代表滤波器i第k+1步的滤波增益矩阵，P_i(k+1)代表滤波器i第k+1步的状态估计误差协方差矩阵,η_i(k+1)代表滤波器i第k+1步的测量信息，即η_x(k+1)＝x(k+1)+υ_x(k+1)，η_y(k+1)＝y(k+1)+υ_y(k+1),η_z(k+1)＝z(k+1)+υ_z(k+1)；

当没有地面目标跟踪装置测量信息时，只运行预报方程(4)，有地面目标跟踪装置测量信息时，则同时运行预报方程(4)和测量修正方程(6)；

步骤三中的步骤(三一)具体如下：

首先，根据三个子滤波器的估计结果得到跟踪结束时刻kT时空间高超声速滑翔弹头在探测点惯性系中的位置[x y z]^T，探测点惯性坐标系中的位置与地心惯性系中的位置关系为

[x_I y_I z_I]^T＝[x y z]^T+[0 R_e 0]^T(38)

其中，R_e＝6378km为地球平均半径；

根据跟踪结束时刻的高超声速滑翔弹头的加速度估计值[a_x a_y a_z]^T，由于探测点惯性坐标系的三个与地心惯性系的三个轴平行，可得

[a_xI a_yI a_zI]^T＝[a_x a_y a_z]^T(39)

其中，a_xI、a_yI和a_zI代表高超声速滑翔弹头的加速度在地心惯性系中的投影

用下式计算出当前高超声速滑翔弹头受到的气动力：

R=m(axIayIazI-gxIgyIgzI)---(40)]]>

在15Mach左右的时候，升阻比大约为3:1，由于R_z1＝0，可得

Ry1=310||R||---(41)]]>

根据跟踪结束时刻的高超声速滑翔弹头的速度估计值[v_x v_y v_z]^T，由于探测点惯性坐标系的三个与地心惯性系的三个轴平行，可得

[v_xI v_yI v_zI]^T＝[v_x v_y v_z]^T(42)

根据跟踪结束时高超声速滑翔弹头的速度[v_xI v_yI v_zI]^T，用式(12)和(13)计算其弹道角θ和ψ_v；

步骤三中的步骤(三二)具体如下：

根据跟踪结束时高超声速滑翔弹头的速度[v_xI v_yI v_zI]^T，用式(11)计算出其速率v，再利用式(9)求得当前的动压头q；

用下式计算出当前的升力系数为

cy=Ry1qS---(43)]]>

攻角为10°时，升力系数为0.37，所以用下式估算当前的攻角：

α=cy0.371057.3---(44)]]>

高超声速滑翔弹头采用倾斜转弯控制方式，即侧滑角β＝0，然后由式(14)中的前两式可以计算出俯仰角和偏航角ψ；

估算出当前的滚转角γ：由式(40)已经计算出气动力在惯性系的投影，把它投影到目标体坐标系，即

Rx1Ry1Rz1=Ct0→1C0→t0R]]>

上式进一步可写作

令

则根据式(45)可得

b3b2-b2b3sinγcosγ=Ry1Rz1---(46)]]>

其中，R_y1由式(41)估算出来，R_z1＝0，用式(46)求得sinγ和cosγ，用下式确定当前的滚转角

γ＝tan2^-1(sinγ,cosγ)(47)。