CN114935277A - 一种滑翔增程制导炮弹理想弹道的在线规划方法 - Google Patents

Abstract

本发明为一种滑翔增程制导炮弹理想弹道的在线规划方法。包括如下步骤：建立地面、平动、弹体、弹道和速度坐标系；建立3dof制导炮弹运动方程；进行风洞试验得到气动参数；根据打击要求建立弹道规划问题的性能指标，建立约束条件；在发射后的初始上升段或助推段开始进行理想弹道的规划，获取并处理此时刻的卫星数据，将数据作为起点代入制导炮弹运动方程；利用自适应hp方法重构网格，对网格数目和插值多项式阶次进行优化，利用伪谱法离散最优控制问题，转化为非线性规划问题，求解得到制导炮弹质心运动的理想弹道曲线、速度曲线、弹道倾角曲线、过载曲线和舵偏角曲线。本发明对制导炮弹在线规划理想弹道前受到的干扰不敏感，且打击精度高。

Description

一种滑翔增程制导炮弹理想弹道的在线规划方法

技术领域

本发明属于制导弹道领域，具体涉及一种滑翔增程制导炮弹理想弹道的在线规划方法。

背景技术

制导炮弹的弹道规划问题在工程技术领域有着广泛的应用，理想弹道指的是制导炮弹从运动起点到终点的满足打击精度和力度要求的轨迹。理想弹道的在线规划是减轻火控***负担和弹道跟踪难度的关键。发射时制导炮弹射角、偏角、初速等参数的偏差，低空上升段的干扰，都会使开始在线规划弹道的时刻制导炮弹的相关运动参数千差万别。滑翔增程制导炮弹理想弹道的在线规划问题是近年来研究的热点。

弹道规划问题本质上是带有状态约束和控制约束的满足某项指标最优的最优控制问题，随着***和约束条件越来越复杂，通常通过数值法求解，有间接法和直接法两种。

文献“高超声速跳跃-滑翔弹道方案设计及优化”(赵吉松等，固体火箭技术，第32卷第2期，第123～126页，2009年)采用间接法将弹道优化问题转化为两点边值问题，并运用共轭梯度法进行解算。文献“基于间接法的上升段轨迹优化方法研究”(吴嘉梁，导航定位与授时，第3卷第2 期，第14～19页，2016年)对于固体火箭上升段轨迹优化问题，采用间接法在满足攻角约束的条件下获得真实大气环境下的最优轨迹。间接法主要是应用Pontryagin极值原理求解最优控制问题的一阶必要条件，得到哈密尔顿多点边值问题，通过数值计算使得***状态、控制输入和参数满足两点边值问题，从而得到最优解。但Pontryagin极值原理推导最优一阶必要条件时引入了协态变量，对初值敏感，局限性较大，且处理复杂最优控制问题时的两点边值难以求解，在实际的工程应用中受到限制。

文献“高超声速无动力远程滑翔飞行器多约束条件下的轨迹快速生成”(雍恩米等，宇航学报，第1期，第46～52页，2008年)使用p法更新网格，改变插值多项式维数实现收敛，并基于高斯伪谱法提出一种可以生成初值的串行求解策略来应对复杂多约束条件。直接法包括打靶法和配点法，通过参数化方法将无线维最优化问题转化为有限维非线性规划问题，通过非线性规划问题求解器直接得到原最优控制问题要求的控制变量，从而解算出最优理想弹道轨迹。打靶法只离散控制变量，不能求解复杂的最优控制问题，因此很少被使用；配点法离散控制变量和状态变量，其中伪谱法因有较好的收敛性而广泛应用，结合p法的伪谱法因节点固定，对于不连续或非光滑的问题难以求解，且过多提高维数时会带来维数灾难。

现在的滑翔增程制导炮弹，是通过火控***的控制完成对目标的准确射击，尤其是舰载滑翔增程制导炮弹，在发射前通过舰炮火控***完成理想弹道的规划再进行射击，实现对岸远距离精确打击的作战任务，通常对落点距目标的圆周误差(CEP)有准确要求，譬如实施例要求对地静止指标 CEP小于20m，射速不小于5发/分钟。根据作战需求舰上备弹数通常达到几百至上千的数量级别，在作战任务紧急时，需要舰上所有舰炮以允许的最大射速进行射击，并完成精确打击，这种情况下，对每一枚发射的制导炮弹都进行理想弹道的规划给火控***带来负担。

发明内容

本发明的目的在于提供一种滑翔增程制导炮弹理想弹道的在线规划方法。

实现本发明目的的技术解决方案为：一种滑翔增程制导炮弹理想弹道的在线规划方法，包括如下步骤：

步骤(1)：建立地面Axyz、平动Oxyz、弹体Ox₁y₁z₁、弹道Ox₂y₂z₂和速度Ox₃y₃z₃坐标系；

步骤(2)：建立3dof制导炮弹运动方程；

步骤(3)：进行风洞试验得到气动参数；

步骤(4)：根据打击要求建立弹道规划问题的性能指标，根据炮台的发射设施能力、弹身的强度限制和打击目标的力度要求建立约束条件；

步骤(5)：在发射后的初始上升段或助推段开始进行理想弹道的规划，获取并处理此时刻的卫星数据，将数据作为起点代入制导炮弹运动方程；

步骤(6)：利用自适应hp方法重构网格，对网格数目和插值多项式阶次进行优化，利用伪谱法离散最优控制问题，转化为非线性规划问题，利用非线性规划求解器进行求解，得到制导炮弹质心运动的理想弹道曲线、速度曲线、弹道倾角曲线、过载曲线和舵偏角曲线。

进一步的，步骤(1)的“建立地面Axyz、平动Oxyz、弹体Ox₁y₁z₁、弹道Ox₂y₂z₂和速度Ox₃y₃z₃坐标系”具体方法为：

地面坐标系固连地面，原点A在发射时弹体的质心上，Ax轴是连接发射装置和目标的水平线，指向目标，Ay轴沿垂线向上，Az轴与Ax和Ay垂直；

O为弹体瞬时质心，Ox、Oy、Oz分别和Ax、Ay、Az平行且同向；Ox₁轴与纵轴重合指向弹头为正，Oy₁轴位于纵向平面，与Ox₁垂直，向上为正，Oz₁轴垂直于Ox₁y₁平面；

Ox₂轴与制导炮弹质心的瞬时速度V重合，Oy₂轴位于含V的铅垂面内与Ox₂轴垂直，Oz₂轴垂直于Ox₂轴与Oy₂轴；Ox₃轴与V重合，Oy₃轴位于纵向对称面内与Ox₃轴垂直，Oz₃轴垂直于 Ox₃y₃平面。

进一步的，步骤(2)建立3dof制导炮弹运动方程，具体方法为：

忽略制导炮弹的滚转，弹道坐标系下制导炮弹质心的动力学方程组为：

地面坐标系下制导炮弹质心的运动学方程为：

质心质量变化方程和控制关系方程为：

式中，m是制导炮弹的质量，

为质量变化率，仅在发动机工作期间即助推期不为零，助推期的开始时间和结束时间取决于制导炮弹本身，m_c为制导炮弹的载药量和发动机燃烧时间有关，是一个具体数值，V是制导炮弹质心的瞬时速度，

是速度变化率，g是重力加速度，P是发动机平均推力，仅存在于助推期，X是制导炮弹所受空气阻力，Y是升力，Z是侧向力，攻角α是速度在弹体坐标系Ox₁y₁z₁上的投影与轴的Ox₁夹角，α_c(t)是控制律作用下攻角α随时间变化的取值，侧滑角β是速度与弹体坐标系Ox₁y₁z₁的夹角，β_c(t)是控制律作用下侧滑角β随时间变化的取值， f_α、f_β是基于瞬时平衡假设展开得到的与气动参数、制导炮弹重心、压心和控制舵的压心有关的攻角、侧滑角和两个舵偏角的关系函数，制导炮弹飞至弹道顶点时才展开控制舵进行控制，速度倾斜角γ_v是速度坐标系的Oy₃轴与弹道坐标系的Ox2y2之间的夹角，弹道倾角θ是V与水平面之间的夹角，

是弹道倾角的变化率，弹道偏角

是V在水平面上的投影和Ox轴之间的夹角，

是弹道倾角的变化率，(x,y,z)为制导炮弹质心在地面坐标系的瞬时位置，

是制导炮弹质心在地面坐标系x轴、y轴、z轴的位置变化率。

进一步的，步骤(3)进行风洞试验得到气动参数，具体方法为：

将滑翔增程制导炮弹模型置于风洞中，根据气体流动与制导炮弹模型的作用，得到多组不同马赫数对应的零升阻力系数、升力系数导数、俯仰舵和偏航舵零升阻力系数、俯仰舵和偏航舵升力系数导数、弹翼组合压心位置、控制舵压心位置的数据，并选取诱导阻力系数和控制舵诱导阻力系数。

进一步的，步骤(4)具体为：

使射程最大时，选取制导炮弹飞行终点的x轴坐标作为终值型性能指标，表达式如下：

J＝-x(t_f)

其中，J代表性能指标，t_f表示飞完理想弹道所需的时间，x(t_f)表示弹道终点在地面坐标系Ax 轴上的坐标；

设置制导炮弹起点参数值为起点约束，设置飞行过程中状态变量和控制变量的最大值和最小值为过程约束，制导炮弹落地时的落角和落速应不小于用户设定的终点约束。

进一步的，步骤(5)具体为：

北斗卫星基于WGS-84坐标系下表示制导炮弹的位置P点(x_ε,y_ε,z_ε)转换为大地坐标系中的位置数据

式中，a表示WGS-84椭球长半轴，数值为6378137.0m；e为WGS-84第一偏心率，平方为 0.0066943799；N为椭球卯酉圈曲率半；WGS-84坐标系的原点w为地球质心，wz轴指向椭球北极，wx轴指向赤道平面和格林尼治子午面的交点，wy轴在赤道平面内和xwz构成右手***坐标系；大地坐标系的椭球短半轴和地球自转轴wz重合，经度λ是P点所在的椭球子午面和格林尼治子午面的夹角，纬度

是过P点的椭球法线和赤道平面的夹角，高度h是沿P点椭球法线到椭球面的距离；

制导炮弹在大地坐标系的坐标转换为在北天东坐标系中的坐标(x_n,y_n,z_n)：

式中，(x_ε0,y_ε0,z_ε0)是北天东坐标系原点在WGS-84坐标系里的位置，北天东坐标系的原点n 为发射装置点，xnz为原点所在地的大地水平面，nx轴指向地球北，nz轴指向地球东，ny轴垂直于xnz面指向空中；

北天东坐标系到射击坐标系(x,y,z)为：

式中，α为北天东坐标系下的射向，射击坐标系即地面坐标系。

进一步的，步骤(6)具体为：

根据3、步骤(2)的动力学、运动学方程组和质量变化方程组成得到最优控制的运动状态方程约束

根据步骤(4)得到性能指标J和约束条件，最终确定最优控制问题的描述模型如下：

式中，J代表性能指标，Φ是终值型性能指标，t_f表示飞完理想弹道所需的时间，

分别表示所有状态变量、所有控制变量构成的向量，4、控制变量为两个舵偏角关于时间的导数，舵偏角的初始值为0°，积分求得舵偏角数值，

表示

对时间的导数，f是一个映射，通过映射用

和

表示

q为

和

的函数，

表示弹道终点的状态向量，

表示弹道起点的状态向量，φ是边界事件约束函数，包括对起点状态向量x(t₀)和终点状态向量x(t_f)的约束；c₁和c₂为等式和不等式约束条件。

进一步的，步骤(6)中的“利用自适应hp方法重构网格，对网格数目和插值多项式阶次进行优化，利用伪谱法离散最优控制问题，转化为非线性规划问题”具体包括如下步骤：

步骤(61)：进行网格初始划分；设定最大容许误差ε、最大容许曲率r_max、初始网格数目和拉格朗日插值多项式的阶次；

将时间区间[t₀,t_f]划分成K个网格子区间[t_k-1,t_k],k＝1,…,K，引入新的时间变量τ∈[-1,1]，通过下式将原变量t转换为新的变量τ，即可将单区间非线性最优控制问题转化为多区间非线性最优控制问题：

步骤(62)：计算当前各网格子区间内的LGR插值点的位置、LGR权重和Radau伪谱微分矩阵，结合拉格朗日插值多项式近似表达出最优控制模型，将其转化为非线性规划问题；

设第k个网格子区间[t_k-1,t_k]内的状态变量和控制变量表示为X^k(τ)和U^k(τ)，N_k+1个离散点

和末点

构成了区间τ∈[-1,1]内的N_k+2个拉格朗日插值点，其中

设

为N_k+2阶拉格朗日插值基函数，

为N_k+1阶拉格朗日插值基函数，表达式为：

状态变量X^(k)在N_k+2个离散点上用

近似为：

式中，

是状态变量X^k在离散点

处的取值；

控制变量U^(k)在N_k+1个离散点上用

近似为：

式中，

是状态变量U^k在离散点

处的取值；

对状态变量X^(k)求导得到近似状态方程的导数为：

令

表示第k个子区间的 Radau伪谱微分近似矩阵，其具体表达式为：

式中，g^(k)(τ)是P_n(τ)的函数，P_n(τ)是勒让德正交多项式，微分表达式为：

则N_k+1个插值点上的状态方程约束转化为代数方程约束，表达为：

式中F_s ^(k)是关于状态变量和控制变量的状态方程函数向量；

用LGR积分近似表达博尔查型性能指标里的积分项，近似为：

式中，

是状态向量和控制向量的函数，

为第k个子区间的离散点

对应的LGR权重，具体形式为：

最终，将最优控制问题离散为非线性规划问题，表达为：

其中，φ是边界事件约束函数，C₁和C₂为等式和不等式约束条件，

是第k个网格的末点的状态变量，

是第k+1个网格的起点的状态变量；

步骤(63)计算第k个网格子区间的最大相对误差

若

则进行下一步，否则跳至步骤(65)；

设第k个网格的N_k+1个采样点为：

状态约束方程和过程约束方程在采样点上的残差表示为：

其中

和

是残差，

是状态变量和控制变量在采样点

上的取值；

取当前区间的最大残差为最大相对误差：

步骤(64)：计算当前网格的相对曲率r^(k)，若r^(k)≤r_max，则重新计算插值多项式的阶次；若 r^(k)＞r_max，则将当前网格细分，计算新插值点的位置；然后再进行步骤(62)到步骤(64)，直至

相对曲率r^(k)的表达式为：

式中，

为第k个网格内所有LGR点的最大曲率和平均曲率；

各点曲率

表达式为：

式中，

和

分别表示第k个网格子区间内的最大相对误差

对应的LGR点处的状态函数的一阶导函数和二阶导函数；

若增加多项式的阶数，则增加的数量ΔN可以表示为：

其中，ceil(·)表示向上取整，B是值大于0且可调的正整数；

细化网格，则该区间重新划分的子区间数量n_k由下式求得：

其中，λ是值大于0的整数调节因子；

步骤(65)：令k＝k+1，对下一个网格进行步骤(64)和步骤(65)，直至所有网格子区间的

均满足要求，则输出结果，退出程序。

一种滑翔增程制导炮弹，包括弹载计算机，弹载计算机采用权利要求1-8任一项所述的方法进行滑翔增程制导炮弹理想弹道的在线规划。

本发明与现有技术相比，其显著优点在于：

(1)相比于其他满足打击目标要求的理想弹道规划方法，本发明在每一个制导炮弹均通过弹载计算机进行在线的弹道规划，本发明提供的在线理想弹道规划方法减轻了火控***的负担。

(2)针对制导炮弹在初始上升段受到阵风等干扰，到达实际弹道顶点前偏离理想弹道太远的问题，相比于由火控***规划的理想弹道对应的控制律不适用干扰，导致无法击中目标的情况，利用本发明提供的方法在助推段进行在线规划，在线规划时干扰已经改变了制导炮弹的位置和运动状态信息，以此为基础规划的理想弹道与实际弹道误差较小，减小了制导炮弹跟踪理想弹道的难度，能精确击中目标。

(3)相对于间接法带来的初值敏感和两点边值难以求解的问题，本发明提供的技术方法依赖于伪谱法，在连续时间段内用牛顿迭代法取得离散点，并选取合适的多项式阶次再进行方程近似，避开了间接法的两个难点，使伪谱法在求解无穷光滑问题时收敛性好且易于计算。

(4)相对于基于p法的伪谱法无法求解不连续或不光滑问题的问题，本发明提供的技术方法基于自适应hp法更新网格，令网格内采样点上的误差和曲率都满足要求，既保留了h法能捕捉解的非光滑性的特点，也保持p法求解时的指数级收敛速度，自适应hp方法提高了非线性规划问题的计算效率和解的精度。

附图说明

图1为本发明建立的弹体坐标系和速度坐标系间的角度关系图。

图2为本发明建立的平动坐标系和弹道坐标系间的角度关系图。

图3为本发明建立的弹道坐标系和速度坐标系间的角度关系图。

图4为在线规划的制导炮弹理想弹道的纵向飞行曲线图。

图5为在线规划的制导炮弹理想弹道的横向飞行曲线图。

图6为在线规划的制导炮弹理想弹道的速度曲线图。

图7为在线规划的制导炮弹理想弹道的弹道倾角曲线图。

图8为在线规划的制导炮弹理想弹道的过载曲线图。

图9为在线规划的制导炮弹理想弹道的俯仰舵舵偏角角度变化图。

图10为在线规划的制导炮弹理想弹道的偏航舵舵偏角角度变化图。

图11为阵风干扰时基于离线规划和在线规划的制导炮弹理想弹道图。

具体实施方式

下面结合附图对本发明作进一步详细描述。

为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。

本发明基于最优控制理论，考虑了初始发射偏差和在线规划弹道前制导炮弹受到的干扰，通过自适应伪谱法，能够有效在线规划满足落点、落速和落角要求的理想弹道，为研究滑翔增程制导炮弹理想弹道的在线规划提供了一种高效准确的计算方法。一种滑翔增程制导炮弹理想弹道的在线规划方法，具体包括如下步骤：

步骤(1)的“建立地面Axyz、平动Oxyz、弹体Ox₁y₁z₁、弹道Ox₂y₂z₂和速度Ox₃y₃z₃坐标系”具体方法为：

地面坐标系固连地面，原点A在发射时弹体的质心上，Ax轴是连接发射装置和目标的水平线，指向目标，Ay轴沿垂线向上，Az轴与Ax和Ay垂直；

O为弹体瞬时质心，Ox、Oy、Oz分别和Ax、Ay、Az平行且同向；Ox₁轴与纵轴重合指向弹头为正，Oy₁轴位于纵向平面，与Ox₁垂直，向上为正，Oz₁轴垂直于Ox₁y₁平面；

Ox₂轴与制导炮弹质心的瞬时速度V重合，Oy₂轴位于含V的铅垂面内与Ox₂轴垂直，Oz₂轴垂直于Ox₂轴与Oy₂轴；Ox₃轴与V重合，Oy₃轴位于纵向对称面内与Ox₃轴垂直，Oz₃轴垂直于 Ox₃y₃平面。

平动坐标系、弹体坐标系、弹道坐标系和速度坐标系间的角度关系如图1-3所示。

步骤2，建立3dof制导炮弹运动方程，具体方法为：

忽略制导炮弹的滚转，弹道坐标系下制导炮弹质心的动力学方程组为：

地面坐标系下制导炮弹质心的运动学方程为：

质心质量变化方程和控制关系方程为：

式中，m是制导炮弹的质量，

为质量变化率，仅在发动机工作期间即助推期不为零，助推期的开始时间和结束时间取决于制导炮弹本身，m_c为制导炮弹的载药量和发动机燃烧时间有关，是一个具体数值，V是制导炮弹质心的瞬时速度，

是速度变化率，g是重力加速度，P是发动机平均推力，仅存在于助推期，X是制导炮弹所受空气阻力，Y是升力，Z是侧向力，攻角α是速度在弹体坐标系Ox₁y₁z₁上的投影与轴的Ox₁夹角，α_c(t)是控制律作用下攻角α随时间变化的取值，侧滑角β是速度与弹体坐标系Ox₁y₁z₁的夹角，β_c(t)是控制律作用下侧滑角β随时间变化的取值， f_α、f_β是基于瞬时平衡假设展开得到的与气动参数，制导炮弹重心、压心和控制舵的压心有关的攻角、侧滑角和两个舵偏角的关系函数，制导炮弹飞至弹道顶点时才展开控制舵进行控制，速度倾斜角γ_v是速度坐标系的Oy₃轴与弹道坐标系的Ox₂y₂之间的夹角，弹道倾角θ是V与水平面之间的夹角，

是弹道倾角的变化率，弹道偏角

是V在水平面上的投影和Ox轴之间的夹角，

是弹道倾角的变化率，(x,y,z)为制导炮弹质心在地面坐标系的瞬时位置，

是制导炮弹质心在地面坐标系x轴、y轴、z轴的位置变化率。

步骤(3)进行风洞试验得到气动参数，具体方法为：

将滑翔增程制导炮弹模型置于风洞中，根据气体流动与制导炮弹模型的作用，得到多组不同马赫数对应的零升阻力系数、升力系数导数、俯仰舵和偏航舵零升阻力系数、俯仰舵和偏航舵升力系数导数、弹翼组合压心位置、控制舵压心位置的数据，并选取诱导阻力系数和控制舵诱导阻力系数。

步骤(4)具体为：

工程应用的很多情况需要指定射程，因此这里以最短飞行时间为要求，选取制导炮弹飞行时间 t作为终值型性能指标，表达式如下：

J＝t_f

其中，J代表性能指标，t_f表示飞完理想弹道所需的时间。

设置制导炮弹的初始发射速度为790m/s，弹道倾角为55°，位置为(0km,0km,0km),设置飞行过程中状态变量和控制变量的最大值和最小值为过程约束，本实施例要求飞行过程中攻角、侧滑角和舵偏角的角度不能超过6°，舵偏角每秒钟变化的角度不能超过2°，过载绝对值不能超过3，设置制导炮弹终点参数值为终点约束，制导炮弹落地时的落角应不小于65°或为指定落角，落速应不小于200m/s。

步骤(5)具体为：

北斗卫星基于WGS-84坐标系下表示制导炮弹的位置P点(x_ε,y_ε,z_ε)转换为大地坐标系中的位置数据

式中，a表示WGS-84椭球长半轴，数值为6378137.0m；e为WGS-84第一偏心率，平方为 0.0066943799；N为椭球卯酉圈曲率半径。WGS-84坐标系的原点w为地球质心，wz轴指向椭球北极，wx轴指向赤道平面和格林尼治子午面的交点，wy轴在赤道平面内和xwz构成右手***坐标系；大地坐标系的椭球短半轴和地球自转轴wz重合，经度λ是P点所在的椭球子午面和格林尼治子午面的夹角，纬度

是过P点的椭球法线和赤道平面的夹角，高度h是沿P点椭球法线到椭球面的距离。

制导炮弹在大地坐标系的坐标转换为在北天东坐标系中的坐标(x_n,y_n,z_n)：

式中，(x_ε0,y_ε0,z_ε0)是北天东坐标系原点在WGS-84坐标系里的位置，北天东坐标系的原点n为发射装置点，xnz为原点所在地的大地水平面，nx轴指向地球北，nz轴指向地球东，ny轴垂直于xnz 面指向空中。

北天东坐标系到射击坐标系(x,y,z)为：

式中，α为北天东坐标系下的射向，射击坐标系即为步骤1中的地面坐标系。

步骤(6)具体为：

根据步骤(2)动力学、运动学方程组和质量变化方程组成得到最优控制的运动状态方程约束

根据步骤(4)得到性能指标J和约束条件，最终确定最优控制问题的描述模型如下：

J代表性能指标，Φ是终值型性能指标，t_f表示飞完理想弹道所需的时间，

分别表示所有状态变量、所有控制变量构成的向量，4、控制变量为两个舵偏角关于时间的导数，舵偏角的初始值为0°，积分求得舵偏角数值，

表示

对时间的导数，f是一个映射，通过映射用

和

表示

q为

和

的函数，

表示弹道终点的状态向量，

表示弹道起点的状态向量，φ是边界事件约束函数，包括对起点状态向量x(t₀)和终点状态向量x(t_f)的约束；c₁和c₂为等式和不等式约束条件。

步骤(6)中的“利用自适应hp方法重构网格，对网格数目和插值多项式阶次进行优化，利用伪谱法离散最优控制问题，转化为非线性规划问题”具体包括如下步骤：

步骤(61)：进行网格初始划分；设定最大容许误差ε为0.0001，最大容许曲率r_max为1，初始化网格数目为10、拉格朗日插值多项式的阶次为4。

将时间区间[t₀,t_f]划分成K个网格子区间[t_k-1,t_k],k＝1,…,K，引入新的时间变量τ∈[-1,1]，通过下式将原变量t转换为新的变量τ，即可将单区间非线性最优控制问题转化为多区间非线性最优控制问题

步骤(62)：计算当前各网格子区间内的LGR插值点的位置、LGR权重和Radau伪谱微分矩阵，结合拉格朗日插值多项式近似表达出最优控制模型，将其转化为非线性规划问题；

设第k个网格子区间[t_k-1,t_k]内的状态变量和控制变量表示为X^k(τ)和U^k(τ)，N_k+1个离散点

和末点

构成了区间τ∈[-1,1]内的N_k+2个拉格朗日插值点，其中

设

为N_k+2阶拉格朗日插值基函数，

为N_k+1阶拉格朗日插值基函数，表达式为：

状态变量X^(k)在N_k+2个离散点上用

近似为：

式中，

是状态变量X^k在离散点

处的取值；

控制变量U^(k)在N_k+1个离散点上用

近似为：

式中，

是状态变量U^k在离散点

处的取值；

对状态变量X^(k)求导得到近似状态方程的导数为：

令

表示第k个子区间的 Radau伪谱微分近似矩阵，其具体表达式为：

式中，g^(k)(τ)是P_n(τ)的函数，P_n(τ)是勒让德正交多项式，微分表达式为：

则N_k+1个插值点上的状态方程约束转化为代数方程约束，表达为：

式中F_s ^(k)是关于状态变量和控制变量的状态方程函数向量；

用LGR积分近似表达博尔查型性能指标里的积分项，近似为：

式中，

是状态向量和控制向量的函数，

为第k个子区间的离散点

对应的LGR权重，具体形式为：

最终，将最优控制问题离散为非线性规划问题，表达为：

其中，φ是边界事件约束函数，C₁和C₂为等式和不等式约束条件，

是第k个网格的末点的状态向量，

是第k+1个网格的起点的状态向量；

步骤(63)计算第k个网格子区间的最大相对误差

若

则进行下一步，否则跳至步骤(65)；

设第k个网格的N_k+1个采样点为：

状态约束方程和过程约束方程在采样点上的残差表示为：

其中

和

是残差，

是状态变量和控制变量在采样点

上的取值；

取当前区间的最大残差为最大相对误差：

步骤(64)：计算当前网格的相对曲率r^(k)，若r^(k)≤r_max，则重新计算插值多项式的阶次；若 r^(k)＞r_max，则将当前网格细分，计算新插值点的位置；然后再进行步骤(62)到步骤(64)，直至

相对曲率r^(k)的表达式为：

式中，

为第k个网格内所有LGR点的最大曲率和平均曲率；

各点曲率

表达式为：

式中，

和

分别表示第k个网格子区间内的最大相对误差

对应的LGR点处的状态函数的一阶导函数和二阶导函数；

若增加多项式的阶数，则增加的数量ΔN可以表示为：

其中，ceil(·)表示向上取整，B是值大于0且可调的正整数；

细化网格，则该区间重新划分的子区间数量n_k由下式求得：

其中，λ是值大于0的整数调节因子；

步骤(65)：令k＝k+1，对下一个网格进行步骤(64)和步骤(65)，直至所有网格子区间的

均满足要求，则输出结果，退出程序。

令某型号制导炮弹以-65°打击(65km，0m，0m)的目标，设制导炮弹以初速790m/s，射角 55°发射，发射时为0时刻，上升到弹道顶点的时间为57.8s，弹载计算机在t₀＝20s时刻接收到制导炮弹的参数并开始规划弹道。制导炮弹前20s飞行无干扰时，此刻的参数信息应为：位置(8585m,10170m,0m)、速度622.94m/s、弹道倾角44.114°、弹道偏角为0°。假定炮弹从发射到t₀时刻经历了未知干扰，此刻状态出现了偏差，任意选择四组弹道参数进行仿真，参数如表1所示，距离的单位为(m)，速度的单位为(m/s)，角度的单位为(°)。仿真结果如图4-图10所示。图 4-图10分别表示在线规划的方案弹道的轨迹和运动参数曲线，可以看到以表1四组参数为起点在线规划的弹道均在满足各变量约束的前提下完成了打击目的。

表1随机选取的弹道规划点参数

令弹丸以v＝790m/s，θ＝55°发射，以落角为-65°打击x＝40km处的目标。在18s～28s时加入 0～50m/s方向与弹速相反的阵风，用自适应hp拉道伪谱法规划弹道。第一组试验使用离线规划的控制律，第二组在t＝30s基于本发明的技术方法在线规划弹道。图11为两组的轨迹，其中第一组的落点为35.41km，与目标距离误差太大，打击失败，第二组完成了落角和定点的要求。

一种滑翔增程制导炮弹，包括弹载计算机，弹载计算机采用权利要求1-8任一项所述的方法进行滑翔增程制导炮弹理想弹道的在线规划。

Claims (9)

1.一种滑翔增程制导炮弹理想弹道的在线规划方法，其特征在于，包括如下步骤：

步骤(1)：建立地面Axyz、平动Oxyz、弹体Ox₁y₁z₁、弹道Ox₂y₂z₂和速度Ox₃y₃z₃坐标系；

步骤(2)：建立3dof制导炮弹运动方程；

步骤(3)：进行风洞试验得到气动参数；

步骤(4)：根据打击要求建立弹道规划问题的性能指标，根据炮台的发射设施能力、弹身的强度限制和打击目标的力度要求建立约束条件；

步骤(5)：在发射后的初始上升段或助推段开始进行理想弹道的规划，获取并处理此时刻的卫星数据，将数据作为起点代入制导炮弹运动方程；

步骤(6)：利用自适应hp方法重构网格，对网格数目和插值多项式阶次进行优化，利用伪谱法离散最优控制问题，转化为非线性规划问题，利用非线性规划求解器进行求解，得到制导炮弹质心运动的理想弹道曲线、速度曲线、弹道倾角曲线、过载曲线和舵偏角曲线。

2.根据权利要求1所述的方法，其特征在于，步骤(1)的“建立地面Axyz、平动Oxyz、弹体Ox₁y₁z₁、弹道Ox₂y₂z₂和速度Ox₃y₃z₃坐标系”具体方法为：

地面坐标系固连地面，原点A在发射时弹体的质心上，Ax轴是连接发射装置和目标的水平线，指向目标，Ay轴沿垂线向上，Az轴与Ax和Ay垂直；

O为弹体瞬时质心，Ox、Oy、Oz分别和Ax、Ay、Az平行且同向；Ox₁轴与纵轴重合指向弹头为正，Oy₁轴位于纵向平面，与Ox₁垂直，向上为正，Oz₁轴垂直于Ox₁y₁平面；

Ox₂轴与制导炮弹质心的瞬时速度V重合，Oy₂轴位于含V的铅垂面内与Ox₂轴垂直，Oz₂轴垂直于Ox₂轴与Oy₂轴；Ox₃轴与V重合，Oy₃轴位于纵向对称面内与Ox₃轴垂直，Oz₃轴垂直于Ox₃y₃平面。

3.根据权利要求2所述的方法，其特征在于，步骤(2)建立3dof制导炮弹运动方程，具体方法为：

忽略制导炮弹的滚转，弹道坐标系下制导炮弹质心的动力学方程组为：

地面坐标系下制导炮弹质心的运动学方程为：

质心质量变化方程和控制关系方程为：

式中，m是制导炮弹的质量，

为质量变化率，仅在发动机工作期间即助推期不为零，助推期的开始时间和结束时间取决于制导炮弹本身，m_c为制导炮弹的载药量和发动机燃烧时间有关，是一个具体数值，V是制导炮弹质心的瞬时速度，

是速度变化率，g是重力加速度，P是发动机平均推力，仅存在于助推期，X是制导炮弹所受空气阻力，Y是升力，Z是侧向力，攻角α是速度在弹体坐标系Ox₁y₁z₁上的投影与轴的Ox₁夹角，α_c(t)是控制律作用下攻角α随时间变化的取值，侧滑角β是速度与弹体坐标系Ox₁y₁z₁的夹角，β_c(t)是控制律作用下侧滑角β随时间变化的取值，f_α、f_β是基于瞬时平衡假设展开得到的与气动参数、制导炮弹重心、压心和控制舵的压心有关的攻角、侧滑角和两个舵偏角的关系函数，制导炮弹飞至弹道顶点时才展开控制舵进行控制，速度倾斜角γ_v是速度坐标系的Oy₃轴与弹道坐标系的Ox₂y₂之间的夹角，弹道倾角θ是V与水平面之间的夹角，

是弹道倾角的变化率，弹道偏角

是V在水平面上的投影和Ox轴之间的夹角，

是弹道倾角的变化率，(x,y,z)为制导炮弹质心在地面坐标系的瞬时位置，

是制导炮弹质心在地面坐标系x轴、y轴、z轴的位置变化率。

4.根据权利要求3所述的方法，其特征在于，步骤(3)进行风洞试验得到气动参数，具体方法为：

将滑翔增程制导炮弹模型置于风洞中，根据气体流动与制导炮弹模型的作用，得到多组不同马赫数对应的零升阻力系数、升力系数导数、俯仰舵和偏航舵零升阻力系数、俯仰舵和偏航舵升力系数导数、弹翼组合压心位置、控制舵压心位置的数据，并选取诱导阻力系数和控制舵诱导阻力系数。

5.根据权利要求4所述的方法，其特征在于，步骤(4)具体为：

使射程最大时，选取制导炮弹飞行终点的x轴坐标作为终值型性能指标，表达式如下：

J＝-x(t_f)

其中，J代表性能指标，t_f表示飞完理想弹道所需的时间，x(t_f)表示弹道终点在地面坐标系Ax轴上的坐标；

设置制导炮弹起点参数值为起点约束，设置飞行过程中状态变量和控制变量的最大值和最小值为过程约束，制导炮弹落地时的落角和落速应不小于用户设定的终点约束。

6.根据权利要求5所述的方法，其特征在于，步骤(5)具体为：

北斗卫星基于WGS-84坐标系下表示制导炮弹的位置P点(x_ε,y_ε,z_ε)转换为大地坐标系中的位置数据

式中，a表示WGS-84椭球长半轴，数值为6378137.0m；e为WGS-84第一偏心率，平方为0.0066943799；N为椭球卯酉圈曲率半；WGS-84坐标系的原点w为地球质心，wz轴指向椭球北极，wx轴指向赤道平面和格林尼治子午面的交点，wy轴在赤道平面内和xwz构成右手***坐标系；大地坐标系的椭球短半轴和地球自转轴wz重合，经度λ是P点所在的椭球子午面和格林尼治子午面的夹角，纬度

是过P点的椭球法线和赤道平面的夹角，高度h是沿P点椭球法线到椭球面的距离；

制导炮弹在大地坐标系的坐标转换为在北天东坐标系中的坐标(x_n,y_n,z_n)：

式中，(x_ε0,y_ε0,z_ε0)是北天东坐标系原点在WGS-84坐标系里的位置，北天东坐标系的原点n为发射装置点，xnz为原点所在地的大地水平面，nx轴指向地球北，nz轴指向地球东，ny轴垂直于xnz面指向空中；

北天东坐标系到射击坐标系(x,y,z)为：

式中，α为北天东坐标系下的射向，射击坐标系即地面坐标系。

7.根据权利要求6所述的方法，其特征在于，步骤(6)具体为：

根据3、步骤(2)的动力学、运动学方程组和质量变化方程组成得到最优控制的运动状态方程约束

根据步骤(4)得到性能指标J和约束条件，最终确定最优控制问题的描述模型如下：

式中，J代表性能指标，Φ是终值型性能指标，t_f表示飞完理想弹道所需的时间，

分别表示所有状态变量、所有控制变量构成的向量，控制变量为两个舵偏角关于时间的导数，舵偏角的初始值为0°，积分求得舵偏角数值，

表示

对时间的导数，f是一个映射，通过映射用

和

表示

q为

和

的函数，

表示弹道终点的状态向量，

表示弹道起点的状态向量，φ是边界事件约束函数，包括对起点状态向量x(t₀)和终点状态向量x(t_f)的约束；c₁和c₂为等式和不等式约束条件。

8.根据权利要求7所述的方法，其特征在于，步骤(6)中的“利用自适应hp方法重构网格，对网格数目和插值多项式阶次进行优化，利用伪谱法离散最优控制问题，转化为非线性规划问题”具体包括如下步骤：

步骤(61)：进行网格初始划分；设定最大容许误差ε、最大容许曲率r_max、初始网格数目和拉格朗日插值多项式的阶次；

将时间区间[t₀,t_f]划分成K个网格子区间[t_k-1,t_k],k＝1,…,K，引入新的时间变量τ∈[-1,1]，通过下式将原变量t转换为新的变量τ，即可将单区间非线性最优控制问题转化为多区间非线性最优控制问题：

步骤(62)：计算当前各网格子区间内的LGR插值点的位置、LGR权重和Radau伪谱微分矩阵，结合拉格朗日插值多项式近似表达出最优控制模型，将其转化为非线性规划问题；

设第k个网格子区间[t_k-1,t_k]内的状态变量和控制变量表示为X^k(τ)和U^k(τ)，N_k+1个离散点

和末点

构成了区间τ∈[-1,1]内的N_k+2个拉格朗日插值点，其中

设

为N_k+2阶拉格朗日插值基函数，

为N_k+1阶拉格朗日插值基函数，表达式为：

状态变量X^(k)在N_k+2个离散点上用

近似为：

式中，

是状态变量X^k在离散点

处的取值；

控制变量U^(k)在N_k+1个离散点上用

近似为：

式中，

是状态变量U^k在离散点

处的取值；

对状态变量X^(k)求导得到近似状态方程的导数为：

令

表示第k个子区间的Radau伪谱微分近似矩阵，其具体表达式为：

式中，g^(k)(τ)是P_n(τ)的函数，P_n(τ)是勒让德正交多项式，微分表达式为：

则N_k+1个插值点上的状态方程约束转化为代数方程约束，表达为：

式中F_s ^(k)是关于状态变量和控制变量的状态方程函数向量；

用LGR积分近似表达博尔查型性能指标里的积分项，近似为：

式中，

是状态向量和控制向量的函数，

为第k个子区间的离散点

对应的LGR权重，具体形式为：

最终，将最优控制问题离散为非线性规划问题，表达为：

其中，φ是边界事件约束函数，C₁和C₂为等式和不等式约束条件，

是第k个网格的末点的状态变量，

是第k+1个网格的起点的状态变量；

步骤(63)计算第k个网格子区间的最大相对误差

若

则进行下一步，否则跳至步骤(65)；

设第k个网格的N_k+1个采样点为：

状态约束方程和过程约束方程在采样点上的残差表示为：

其中

和

是残差，

是状态变量和控制变量在采样点

上的取值；

取当前区间的最大残差为最大相对误差：

步骤(64)：计算当前网格的相对曲率r^(k)，若r^(k)≤r_max，则重新计算插值多项式的阶次；若r^(k)＞r_max，则将当前网格细分，计算新插值点的位置；然后再进行步骤(62)到步骤(64)，直至

相对曲率r^(k)的表达式为：

式中，

为第k个网格内所有LGR点的最大曲率和平均曲率；

各点曲率

表达式为：

式中，

和

分别表示第k个网格子区间内的最大相对误差

对应的LGR点处的状态函数的一阶导函数和二阶导函数；

若增加多项式的阶数，则增加的数量ΔN可以表示为：

其中，ceil(·)表示向上取整，B是值大于0且可调的正整数；

细化网格，则该区间重新划分的子区间数量n_k由下式求得：

其中，λ是值大于0的整数调节因子；

步骤(65)：令k＝k+1，对下一个网格进行步骤(64)和步骤(65)，直至所有网格子区间的

均满足要求，则输出结果，退出程序。

9.一种滑翔增程制导炮弹，其特征在于，包括弹载计算机，弹载计算机采用权利要求1-8任一项所述的方法进行滑翔增程制导炮弹理想弹道的在线规划。

Images