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

摘要

一种临近空间高超声速滑翔弹头跳跃弹道预测方法，本发明涉及一种飞行器弹道预测方法。本发明是要解决现有方法对机动目标弹道预测精度低的问题，而提供一种临近空间高超声速滑翔弹头跳跃弹道预测方法。它按下述步骤实现一、建立高超声速滑翔弹头的弹道方程；二、设计实时跟踪高超声速滑翔弹头运动轨迹的卡尔曼滤波器；三、根据跟踪结束时刻的高超声速滑翔弹头的位置、速度和加速度，结合弹道方程估算高超声速滑翔弹头飞行时的攻角和滚转角，在随后的飞行时间内临近空间高超声速滑翔弹头进行等攻角和等滚转角飞行，应用弹道方程向下一时刻循环递推计算，得到一定时间之后的高超声速滑翔弹头的弹道预测值。属于目标跟踪技术领域。

主权项

一种临近空间高超声速滑翔弹头跳跃弹道预测方法，其特征在于它按以下步骤实现：一、高超声速滑翔弹头的弹道方程的建立；建立临近空间高超声速滑翔弹头的弹道方程，包括气动力计算方程、加速度计算方程、速度和位置计算方程、速率和弹道角计算方程，以及姿态角计算方程；综合计算方程，得到高超声速滑翔弹头质心运动速度和位置方程如式(1)～(3)所示：x·I=vxIv·xI=-μxIrI3+aaxI---(1)]]>y·I=vyIv·yI=-μyIrI3+aayI---(2)]]>z·I=v·zIv·zI=-μzIrI3+aazI---(3)]]>其中，xI、yI和zI代表高超声速滑翔弹头在地心惯性坐标系中的位置，代表高超声速滑翔弹头的位置到地球中心的距离,和则代表高超声速滑翔弹头在地心惯性坐标系中的速度，aaxI、aayI和aazI代表气动力产生的高超声速滑翔弹头机动加速度在地心惯性坐标系中的分量，μ为引力参数；二、高超声速滑翔弹头的跟踪滤波器设计；设计高超声速滑翔弹头的运动轨迹跟踪卡尔曼滤波器，实时地跟踪并估计高超声速滑翔弹头的位置、速度和加速度；其中，所述高超声速滑翔弹头的运动轨迹跟踪卡尔曼滤波器的状态一步预报方程为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步的状态预报值，Pi(k)代表滤波器i第k步的状态估计误差协方差矩阵，Pi(k+1/k)代表滤波器i从第k步到第k+1步的状态预报误差协方差矩阵，Qi为滤波器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的测量矩阵Hi＝[1 0 0],i＝x,y,z  (7)Ri为滤波器i测量误差协方差矩阵，I代表单位矩阵，Ki(k+1)代表滤波器i第k+1步的滤波增益矩阵，ηi(k+1)代表滤波器i第k+1步的测量信息，Pi(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时刻的位置xI、yI和zI、kT时刻的速度vxI、vyI和vzI代入式(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]]>cy为升力系数，cx为阻力系数，S为特征参考面积；q＝ρv2/2  (9)其中，q为动压头，v为高超声速滑翔弹头的速率，ρ为空气密度，它随着飞行高度h变化，可查标准的大气密度表得到；高超声速滑翔弹头飞行高度h用下式计算:h=xI2+yI2+zI2-Re---(10)]]>其中，xI、yI和zI代表高超声速滑翔弹头在地心惯性坐标系中的位置，Re＝6378km为地球平均半径；设在飞行中侧滑角β＝0，因此气动力侧向力Rz1＝0；(二)计算高超声速滑翔弹头的速率高超声速滑翔弹头的速率为v=vxI2+vyI2+vzI2---(11)]]>其中，vxI、vyI和vzI代表高超声速滑翔弹头在地心惯性坐标系中的速度；将高超声速滑翔弹头的速度矢量投影到目标惯性坐标系vxt0、vyt0与vzt0，即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)]]>其中，gxI、gyI与gzI代表重力加速度在地心惯性系中的投影；步骤二具体为：(一)在探测点惯性坐标系的三个轴上，均用Singer模型来分别描述高超声速滑翔弹头的加速度分量的变化，即a·x=-λxax+wxa·y=-λyay+wya·z=-λzaz+wz---(17)]]>其中，λx、λy和λz分别代表在探测点惯性坐标系下oFx轴、oFy轴和oFz轴方向上目标机动时间常数的倒数，wx、wy和wz分别代表探测点惯性坐标系下oFx轴、oFy轴和oFz轴方向上的零均值高斯白噪声；加速度ax、ay和az中既包含气动力产生的加速度aaxI、aayI和aazI，又包含重力产生的加速度gxI、gyI和gzI；探测点惯性坐标系的oFx轴、oFy轴和oFz轴与地心惯性系的三个轴oIxI、oIyI和oIzI分别平行，气动力产生的加速度在地心惯性系中的投影aaxI、aayI和aazI就是此加速度在探测点惯性系中的投影；(二)在探测点惯性坐标系下建立起如下三组运动方程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代表高超声速滑翔弹头在探测点惯性坐标系中的位置，和则代表高超声速滑翔弹头速度在探测点惯性坐标系中的投影；因为探测点惯性坐标系的oFx轴、oFy轴和oFz轴与地心惯性系的三个轴oIxI、oIyI和oIzI分别平行，所以vx=vxIvy=vyIvz=vzI---(21)]]>而由于探测点惯性坐标系只是相对于地心惯性系探测点惯性坐标系在oIyI轴方向上作了平移，所以有x=xIy=yI-Rez=zI---(22)]]>对于系统(18)，定义x轴子系统状态变量Xx＝[x vx ax]T，则得到如下状态方程X·x=FxXx+Gxwx---(23)]]>其中，Fx=01000100-λx,Gx=001---(24)]]>对式(23)以一定的周期Δt进行离散化后，得到从第k步到第k+1步的状态方程为Xx(k+1)＝ΦxXx(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轴子系统状态变量Xy＝[y vy ay]T，则得到如下状态方程X·y=FyXy+Gywy---(28)]]>其中，Fy=01000100-λy,Gy=001---(29)]]>对式(28)以一定的周期Δt进行离散化后，得到从第k步到第k+1步的状态方程为Xy(k+1)＝ΦyXy(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轴子系统状态变量Xz＝[z vz az]T，则得到如下状态方程X·z=FzXz+Gzwz---(33)]]>其中，Fz=01000100-λz,Gz=001---(34)]]>对式(28)以一定的周期Δt进行离散化后，得到从第k步到第k+1步的状态方程为Xz(k+1)＝ΦzXz(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步的状态预报值，Pi(k)代表滤波器i第k步的状态估计误差协方差矩阵，Pi(k+1/k)代表滤波器i从第k步到第k+1步的状态预报误差协方差矩阵，Qi为滤波器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),即Hi＝[1 0 0],i＝x,y,zRi为滤波器i测量误差协方差矩阵，I代表单位矩，Ki(k+1)代表滤波器i第k+1步的滤波增益矩阵，Pi(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，探测点惯性坐标系中的位置与地心惯性系中的位置关系为[xI yI zI]T＝[x y z]T+[0 Re 0]T  (38)其中，Re＝6378km为地球平均半径；根据跟踪结束时刻的高超声速滑翔弹头的加速度估计值[ax ay az]T，由于探测点惯性坐标系的三个与地心惯性系的三个轴平行，可得[axI ayI azI]T＝[ax ay az]T  (39)其中，axI、ayI和azI代表高超声速滑翔弹头的加速度在地心惯性系中的投影用下式计算出当前高超声速滑翔弹头受到的气动力：R=m(axIayIazI-gxIgyIgzI)---(40)]]>在15Mach左右的时候，升阻比大约为3:1，由于Rz1＝0，可得Ry1=310||R||---(41)]]>根据跟踪结束时刻的高超声速滑翔弹头的速度估计值[vx vy vz]T，由于探测点惯性坐标系的三个与地心惯性系的三个轴平行，可得[vxI vyI vzI]T＝[vx vy vz]T  (42)根据跟踪结束时高超声速滑翔弹头的速度[vxI vyI vzI]T，用式(12)和(13)计算其弹道角θ和ψv；步骤三中的步骤(三二)具体如下：根据跟踪结束时高超声速滑翔弹头的速度[vxI vyI vzI]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)]]>其中，Ry1由式(41)估算出来，Rz1＝0，用式(46)求得sinγ和cosγ，用下式确定当前的滚转角γ＝tan2‑1(sinγ,cosγ)  (47)。