强约束复杂任务条件下的三维剖面优化设计方法、系统及
介质
技术领域
[0001] 本
发明涉及
飞行器动
力学、弹道规划与制导领域，具体涉及一种
强约束复杂任务条件下的三维剖面优化设计方法、系统及介质，可广泛应用于高超声速飞行器、载人飞船等飞行器的强约束复杂任务弹道生成，为其弹道规划和制导方法提供
支撑。
背景技术
[0002] 目前弹道优化设计方法主要包括两大类，即直接法和间接法。直接法通常是指通过数值积分或其他搜索
算法，直接找到满足问题的最优解。在利用直接法进行弹道优化求解的计算中，以高斯伪谱法、Radau伪谱法等应用较为著名，比如Yu Zhengshi针对火星再入轨迹优化问题，将伪谱法应用到了考虑干扰和初始状态偏差轨迹优化问题中；而Ma Lin则将高斯伪谱法应用到月球上升段常推力下的燃油轨迹优化设计；针对禁飞区和
航路点等复杂约束条件下的弹道优化问题，Zhao Jiang基于分段高斯伪谱法优化的策略进行求解等。这类方法虽然都已在弹道优化问题中取得了较好的解，但由于伪谱法本身求解时需要对整个求解空间的所有状态量和控制量同时进行优化计算，求解
精度高但是计算效率低，有时甚至无法保证问题的解稳定收敛。间接法一般是通过建立必要的微分方程、目标函数以及约束条件，并基于一定的假设和简化，解析推导得出满足问题的最优表达式。对于
滑翔飞行器弹道优化设计而言，由于问题本身较为复杂，在进行大量的假设和简化后获得的解析解通常适应性较差。因此，本发明主要基于直接法开展三维剖面的弹道优化设计。
[0003] 高超声速滑翔飞行器弹道优化设计问题本质上属于非线性复杂多约束条件下的最优化问题。高度作为纵向运动状态的典型代表，通过建立其与时间或
能量的函数关系剖面后可基本确定速度、速度倾
角等纵向运动量以及动压和热流等影响飞行器关键结构性能考核的约束指标。同时，航路点和禁飞区的典型路径约束已成为当前滑翔飞行器飞行试验和机动性能考核的基本任务需求。因此，研究给定高度剖面约束下航路点和禁飞区约束的弹道优化设计方法具有重要的理论意义和应用价值，尤其是脱离事先给定参考
攻角剖面约束下实现控制能力按需分配的弹道优化设计方法对于充分发挥飞行器机动潜能更是具有极大地促进作用。目前，关于非线性最优设计问题尚未有统一的解决算法，一般做法是通过将连续问题离散化后转换为非线性规划问题进行求解。序列二次规划方法是当前用于求解复杂约束非线性规划问题的公认求解最有效的算法之一，具有类似
牛顿
迭代算法的超一次收敛速率和
稳定性，这也是本发明将其选为三维剖面弹道优化设计的重要原因。
发明内容
[0004] 本发明要解决的技术问题：给定高度剖面约束，并考虑航路点和禁飞区等复杂飞行任务需求，兼顾飞行过程中的动压、过载以及驻点热流等过程约束、控制量大小和变化速率约束以及初终端条件而开展的滑翔段可行弹道设计问题。为了充分发挥滑翔飞行器的机动潜能，要求脱离参考攻角剖面，即产生的三维弹道应能实现控制能力的按需分配。
[0005] 为了解决上述技术问题，本发明采用的技术方案为：
[0006] 一种强约束复杂任务条件下的三维剖面优化设计方法，实施步骤包括：
[0007] 1)综合考虑多种约束条件，构建满足给定高度剖面约束的侧向
升阻比走廊边界；
[0008] 2)初始化侧向剖面参数，并利用侧向升阻比走廊边界进行归一化得到初始化侧向剖面；
[0009] 3)根据给定高度剖面和初始化侧向剖面解算参考弹道；
[0010] 4)建立考虑复杂多约束的非线性规划模型，基于参考弹道进行求解得到强约束复杂任务条件下的三维剖面。
[0011] 可选地，步骤4)的详细步骤包括：
[0012] 4.1)输入强高度剖面约束、复杂飞行任务，提取需要参数信息，包括：参考高度剖面、初始终端状态条件以及过程约束；
[0013] 4.2)设计侧向剖面的初值生成器，对侧向剖面参数进行初始化产生初始侧向剖面；
[0014] 4.3)根据给定的高度剖面约束和初始侧向剖面，求解纵向剖面获得参考三维剖面；
[0015] 4.4)判定解算得到的参考三维剖面是否满足给定的约束要求，满足则跳转执行下一步；否则，检查约束条件和联系方程是否合理，调整侧向剖面或相应的约束条件重新跳转执行步骤4.2)直至产生的三维剖面满足约束要求；
[0016] 4.5)建立融入航路点、禁飞区以及目标点的混合性能指标，联立三维剖面和相关约束条件形成非线性弹道规划的数学求解模型；
[0017] 4.6)利用非线性序列二次规划算法求解非线性弹道规划的数学求解模型，并输出求解得到的参考弹道；
[0018] 4.7)判断求解得到的参考弹道是否满足任务要求，若满足任务要求，则输出参考弹道；否则通过非线性规划求解器自动调整剖面参数，返回4.2)继续进行迭代计算，直到求解得到的参考弹道满足任务要求。
[0019] 可选地，步骤4.2)的详细步骤包括：
[0020] 4.2.1)根据飞行任务，判定初终端符号；根据确定的初终端符号，建立相应的侧向和纵向走廊边界；
[0021] 4.2.2)假定从滑翔起点到终点间的飞行路径为理想大圆弧，则根据式(4-1)所示球面几何三角变换可求得需要的总飞行航程Stogo；
[0022] stogo＝cos-1(sinφfsinφ+cosφfcosφcos(λf-λ)) (4-1)
[0023] 式(4-1)中，φf表示目标点地心纬度，φ表示飞行器当前地心纬度，λ表示飞行器当前经度，λf表示目标点经度；
[0024] 4.2.3)记阻力
加速度走廊的最大最小边界分别为Dconmax和Dconmin，根据式(4-2)确定给定初始阻力加速度剖面权重系数ωD时对应的阻力加速度剖面D(e)；
[0025] D(e)＝ωD(Dconmax(e)-Dconmin(e))+Dconmin(e) (4-2)
[0026] 式(4-2)中，D(e)为给定初始阻力加速度剖面权重系数ωD时对应的阻力加速度剖面，为，Dconmax(e)和Dconmin(e)分别为加速度走廊的最大最小边界，ωD为给定初始阻力加速度剖面权重系数；
[0027] 4.2.4)为了快速求出式(4-2)对应的剖面航程，将式(4-2)进行在积分区间内离散，并假设各个剖面点之间采用线性连接，求得此时相邻两点间对应的阻力加速度剖面系数如式(4-3)所示；
[0028]
[0029] 式(4-3)中，aD(i)表示第i点到第(i+1)点间直线连接的阻力加速度剖面斜率，bD(i)表示第i点到第(i+1)点间直线连接的阻力加速度剖面常数项，D(i+1)表示第(i+1)点的阻力加速度，D(i)第i点的表示阻力加速度，e(i+1)表示第(i+1)点的归一化能量值，e(i)表示第i点的归一化能量值；
[0030] 4.2.5)记各阻力加速度点间对应的剖面航程为SD(i)，则可求得剖面对应的总航程ST如式(4-4)所示；
[0031]
[0032] 式(4-4)中，Nint表示离散点数目，SD(i)表示从第i点到第(i+1)点对应的剖面航程，aD(i)表示第i点到第(i+1)点间直线连接的阻力加速度剖面斜率，bD(i)表示第i点到第(i+1)点间直线连接的阻力加速度剖面常数项，e(i+1)表示第(i+1)点的归一化能量值，e(i)表示第i点的归一化能量值；
[0033] 4.2.6)迭代给定初始阻力加速度剖面权重系数ωD，当剖面对应的总航程ST、总飞行航程Stogo成立时即找到初始阻力加速度剖面；
[0034] 4.2.7)利用侧向走廊大小边界，通过插值获得对应的侧向剖面大小；根据初始判定的初终端符号对侧向剖面的符号进行调整，初始正、终端负时将求得的侧向指令大小剖面最小值点为边界，且按照初始正、终端负确定侧向剖面的符号；然后再从找到的初始剖面中反求出需要的初始迭代点，并进行归一化，从而得到需要的初值迭代点，产生初始侧向剖面。
[0035] 可选地，步骤4.5)的详细步骤包括：将侧向升阻比剖面离散化为一系列参考点，各参考点间进行线性连接产生标准侧向指令剖面；当确定纵、侧向剖面后，将其代入侧向运动模型中积分求解余下的侧向运动状态量；设定目标函数为与终端目标点的距离最小，则可建立如式(4-5)所示的非线性弹道规划的数学求解模型；
[0036] min J＝min(|λ(ef)-λf|+|φ(ef)-φf|)
[0037]
[0038] 式(4-5)中，J表示性能指标函数，λ(ef)和φ(ef)分别表示终端归一化能量ef对应的经度和纬度，λf和φf分别表示目标点经度和纬度；h(e)表示高度剖面，fh(e)表示给定的强高度剖面函数，L表示升力加速度，υ表示
倾侧角，g为
重力加速度，V是飞行器速度，r为地心距， 表示弹道状态量的导数，x(e)表示弹道状态量，u(e)表示控制变量，u(e)＝[D(e),Lz(e)]T，e表示归一化能量，f(x(e),u(e),e)表示用归一化能量表示的运动
状态方程，e0表示初始归一化能量，ef表示终端归一化能量， 表示能量e时刻的驻点热流
密度，q(e)表示能量e时刻的动压，L(e)表示能量e时刻的升力加速度，D(e)表示能量e时刻的阻力加速度，g0为海平面重力加速度， 以及nmax分别为最大热流密度、最大动压及最大过载。
[0039] 可选地，步骤1)中构建满足参考高度剖面约束的侧向升阻比走廊边界时，参考高度剖面约束的函数表达式如式(1-1)所示，侧向升阻比走廊边界的变化区间的函数表达式如式(1-2)所示；
[0040] h(e)＝(hf-h0)(e-e0)2+h0 (1-1)
[0041] 式(1-1)中，h(e)表示参考高度剖面，h0表示滑翔段初始高度，hf表示滑翔段终端高度，e表示归一化能量，e0表示归一化初始能量；
[0042] Lzαmin≤|Lz|≤Lzαmax (1-2)
[0043] 式(1-2)中，Lz表示侧向升阻比走廊边界，Lzαmin表示最小攻角对应的侧向升阻比走廊边界，Lzαmax表示最大攻角对应的侧向升阻比走廊边界；Lzαmin和Lzαmax的函数表达式如式(1-3)所示；
[0044]
[0045] 式(1-3)中，Lzαmin表示最小攻角对应的侧向升阻比走廊边界，Lzαmax表示最大攻角对应的侧向升阻比走廊边界，LD表示总升阻比，Lu(e)表示纵向升力加速度关于能量的函数，D表示阻力加速度。
[0046] 可选地，步骤2)的详细步骤包括：
[0047] 2.1)记待优化侧向剖面参数的个数为N，假设侧向剖面是通过侧向剖面参数点之间线性连接得到，得侧向剖面各连接点的函数表达式如式(2-1)所示；
[0048]
[0049] 式(2-1)中，aLz(i)表示第i个连接点的连线斜率，bLz(i)表示第i个连接点的连线常数项，Lz(i+1)表示第i+1个连接点的侧向升阻比，Lz(i)表示第i个连接点的侧向升阻比，e(i)表示第i个连接点的归一化能量，e(i+1)表示第i+1个连接点的归一化能量；
[0050] 2.2)根据式(2-2)利用侧向升阻比走廊边界对侧向剖面参数进行归一化；
[0051]
[0052] 式(2-2)中，xz(N)表示归一化侧向升阻比，Lz表示侧向升阻比，Lzdn(N)表示侧向升阻比走廊下边界，Lzup(N)表示侧向升阻比走廊上边界；
[0053] 2.3)通过搜索归一化侧向升阻比xz(N)，当其确定的如式(2-1)所示的剖面满足给定的约束方程同时使性能指标达到最小时，即可将其作为得到的初始化侧向剖面。
[0054] 此外，本发明还提供一种强约束复杂任务条件下的三维剖面优化设计系统，包括计算机设备，该计算机设备被编程或配置以执行前述强约束复杂任务条件下的三维剖面优化设计方法的步骤，或该计算机设备的存储介质上存储有被编程或配置以执行前述强约束复杂任务条件下的三维剖面优化设计方法的
计算机程序。
[0055] 此外，本发明还提供一种计算机可读存储介质，该计算机可读存储介质上存储有被编程或配置以执行前述强约束复杂任务条件下的三维剖面优化设计方法的计算机程序。
[0056] 和
现有技术相比，本发明具有下述优点：本发明针对高超声速飞行器滑翔段强约束复杂任务下的三维剖面设计问题，采用侧向优先设计策略，综合兼顾飞行器纵、侧向任务需求，通过建立高度剖面约束与三维剖面的映射关系，并构建融入航路点和禁飞区约束的综合性能指标，从而将复杂三维剖面问题转换为单剖面搜索问题进行求解，从而简化了强约束复杂任务条件下的三维剖面优化设计过程，使得计算更加简单和高效。此外，本发明还可进一步在建立考虑复杂多约束的非线性规划模型时采用初值生成器，以确保利用序列二次规划算法快速获得收敛解，为滑翔飞行器弹道规划与制导方法提供了有效技术支撑。
附图说明
[0058] 图2为本发明实施例的算法流程。
[0059] 图3为本发明实施例的侧向剖面的
跟踪结果。
[0060] 图4为本发明实施例的纵向剖面跟踪结果。
[0061] 图5为本发明实施例的参考高度剖面跟踪结果。
[0062] 图6为本发明实施例的航向角跟踪结果。
[0063] 图7为本发明实施例的攻角曲线。
[0064] 图8为本发明实施例的倾侧角曲线。
[0065] 图9为本发明实施例的航路点约束的优化轨迹。
[0066] 图10为本发明实施例的航路点约束的侧向剖面。
[0067] 图11为本发明实施例的禁飞区约束的优化轨迹。
[0068] 图12为本发明实施例的禁飞区约束的侧向剖面。
[0069] 图13为本发明实施例的两个航路点约束的优化轨迹。
[0070] 图14为本发明实施例的两个航路点约束的三维剖面。
具体实施方式
[0071] 如图1所示，本实施例强约束复杂任务条件下的三维剖面优化设计方法的实施步骤包括：
[0072] 1)综合考虑多种约束条件，构建满足给定高度剖面约束的侧向升阻比走廊边界；
[0073] 2)初始化侧向剖面参数，并利用侧向升阻比走廊边界进行归一化得到初始化侧向剖面；
[0074] 3)根据给定高度剖面和初始化侧向剖面解算参考弹道；
[0075] 4)建立考虑复杂多约束的非线性规划模型，并基于参考弹道进行求解得到强约束复杂任务条件下的三维剖面。
[0076] 本实施例步骤1)中构建满足参考高度剖面约束的侧向升阻比走廊边界时，参考高度剖面约束的函数表达式如式(1-1)所示，侧向升阻比走廊边界的变化区间的函数表达式如式(1-2)所示；
[0077] h(e)＝(hf-h0)(e-e0)2+h0 (1-1)
[0078] 式(1-1)中，h(e)表示参考高度剖面，h0表示滑翔段初始高度，hf表示滑翔段终端高度，e表示归一化能量，e0表示归一化初始能量；式(1-1)即采用
顶点式抛物线作为参考高度剖面，既满足初终端条件约束，同时使初始点为最大高度点，确保弹道高度平滑降低，相应的参考速度可以通过能量定义求得；
[0079] Lzαmin≤|Lz|≤Lzαmax (1-2)
[0080] 式(1-2)中，Lz表示侧向升阻比走廊边界，Lzαmin表示最小攻角对应的侧向升阻比走廊边界，Lzαmax表示最大攻角对应的侧向升阻比走廊边界；Lzαmin和Lzαmax的函数表达式如式(1-3)所示；
[0081]
[0082] 式(1-3)中，Lzαmin表示最小攻角对应的侧向升阻比走廊边界，Lzαmax表示最大攻角对应的侧向升阻比走廊边界，LD表示总升阻比，Lu(e)表示纵向升力加速度关于能量的函数，D表示阻力加速度。
[0083] 基于式(1-1)，修正后的准平衡滑翔条件的函数表达式如式(1-4)所示；
[0084]
[0085] 式(1-4)中，L表示升力加速度，υ表示倾侧角，g表示引力加速度，V表示飞行器速度，r表示地心距，Lhu(e)表示纵向升力加速度关于能量的函数。
[0086] 要想确保式(1-4)恒成立，飞行器升力加速度最小值必须大于等于式(1-4)确定的升力加速度，故满足条件的最小攻角αegmin如式(1-5)：
[0087]
[0088] 式(1-5)中，CL(αegmin,Ma)表示升力系数，αegmin表示最小攻角，Ma表示
马赫数，m表示飞行器
质量，L表示升力加速度，υ表示倾侧角，q表示动压，Sref表示飞行器参考面积。因此，飞行器的最小允许攻角αfeamin如式(1-6)：
[0089] αfeamin＝max(αegmin,αmin) (1-6)
[0090] 式(1-6)中，αegmin表示飞行器的最小允许攻角，αmin为给定的最小攻角。
[0091] 根据侧向指令的定义，可求得其可达边界的绝对值如式(1-7)：
[0092]
[0093] 式(1-7)中，Lz表示侧向升阻比，LD表示升阻比，υ表示倾侧角，L表示升力加速度。当参考高度曲线确定后根据给定的过程约束和攻角变化范围，可确定纵向飞行走廊。
[0094] 本实施例中，步骤1)综合考虑多种约束条件包括综合考虑热流密度、动压、过载、攻角、倾侧角等约束条件。滑翔飞行器在飞行过程中受到热流密度、过载、动压等复杂过程约束限制表示如式(1-8)：
[0095]
[0096] 式(1-8)中， qmax以及nmax分别为最大热流密度、最大动压及最大过载，kh为常系数，ρ为大气密度，g0为海平面重力加速度，V是飞行器速度，D为
气动阻力加速度，L为气动升力加速度，D、L计算公式分别如式(1-9)：
[0097]
[0098] 式(1-9)中，CD和CL分别表示
气动阻力和升力系数， 表示动压头，Sref为特征面积，M为飞行器质量。fD(α,V,h)表示阻力加速度函数，阻力加速度函数是关于攻角α、飞行器速度V、飞行器当前最小可行高度的函数；fL(α,V,h)表示升力加速度函数，升力加速度函数是关于攻角α、飞行器速度V、飞行器当前最大可行高度的函数。
[0099] 气动阻力加速度D用能量E相关的气动阻力加速度D(E)可表达如式(1-10)：
[0100]
[0101] 式(1-10)中， 以及nmax分别为最大热流密度、最大动压及最大过载，Dq(E)以及Dn(E)分别为最大热流密度、最大动压及最大过载工况下的气动阻力加速度D，阻力加速度系数CD(α,Ma)和升阻比LD(α,Ma)是关于攻角α和马赫数Ma的函数，Sr为特征面积，M为飞行器质量，m为常数，值为3.15，kh为常系数，V是飞行器速度，g0为海平面重力加速度。则有气动阻力加速度D(E)如式(1-11)：
[0102]
[0103] 式(1-10)中，Dmax表示 Dq(E)以及Dn(E)中的最小值， Dq(E)以及Dn(E)分别为最大热流密度、最大动压及最大过载工况下的气动阻力加速度D。
[0104] 攻角α、倾侧角υ的约束可表达如式(1-12)：
[0105]
[0106] 式(1-12)中，αmin,αmax分别为给定的最小攻角和最大攻角，υmin,υmax分别为给定的最小倾侧角和最大倾侧角。
[0107] 在给定的能量条件下，不同的攻角对应不同的最大阻力加速度Dmax和升阻比LD，根据平衡滑翔假设平衡滑翔条件式(1-13)：
[0108]
[0109] 式(1-13)中，g为重力加速度，V是飞行器速度，r为地心距，L为气动升力，θ为速度倾角，υ为飞行器倾侧角。由于滑翔段速度倾角一般保持在0附近，故通常取cosθ＝1。
[0110] 记Dconmax和Dconmin分别表示给定参考高度的纵向走廊最大、最小边界值，Dαmax和Dαmin分别表示根据定义求得的最大和最小可达阻力加速度边界，则给定参考高度下飞行器的最大和最小允许阻力加速度如式(1-14)：
[0111]
[0112] 式(1-14)中，Dconmax和Dconmin分别表示给定参考高度的纵向走廊最大、最小边界值，Dαmax和Dαmin分别表示根据定义求得的最大和最小可达阻力加速度边界，Dhmax和Dhmin分别表示给定参考高度下飞行器的最大和最小允许阻力加速度。
[0113] 从而对侧向指令Lz，有式(1-15)：
[0114]
[0115] 式(1-15)中，Lz表示侧向升阻比，LD表示升阻比，Lhu(e)表示纵向升力加速度关于能量的函数，D表示阻力加速度。
[0116] 因为高度、速度已知，所以式(1-15)的变量只有攻角。令攻角的取值范围为α∈[αfeamin,αmax]，同时阻力加速度满足式(1-14)。记Lzαmin表示最小攻角对应的侧向升阻比，Lzαmax表示最大攻角对应的侧向升阻比，则有式(1-3)，则侧向指令走廊大小的变化区间有式(1-2)，因此结合式(1-2)和式(1-3)可确定此时对应的侧向指令走廊边界的绝对值变化范围。
[0117] 本实施例中，步骤2)的详细步骤包括：
[0118] 2.1)记待优化侧向剖面参数的个数为N，假设侧向剖面是通过侧向剖面参数点之间线性连接得到，得侧向剖面各连接点的函数表达式如式(2-1)所示；
[0119]
[0120] 式(2-1)中，aLz(i)表示第i个连接点的连线斜率，bLz(i)表示第i个连接点的连线常数项，Lz(i+1)表示第i+1个连接点的侧向升阻比，Lz(i)表示第i个连接点的侧向升阻比，e(i)表示第i个连接点的归一化能量，e(i+1)表示第i+1个连接点的归一化能量；
[0121] 2.2)根据式(2-2)利用侧向升阻比走廊边界对侧向剖面参数进行归一化；
[0122]
[0123] 式(2-2)中，xz(N)表示归一化侧向升阻比，Lz表示侧向升阻比，Lzdn(N)表示侧向升阻比走廊下边界，Lzup(N)表示侧向升阻比走廊上边界；
[0124] 2.3)通过搜索归一化侧向升阻比xz(N)，当其确定的如式(2-1)所示的剖面满足给定的约束方程同时使性能指标达到最小时，即可将其作为得到的初始化侧向剖面。本实施例中，给定的约束方程具体可参见下文的式(4-5)。
[0125] 本实施例采用Lockheed-Martin公司1999年设计的高超声速滑翔飞行器概念模型CAV-H为对象开展仿真实验。本实施例中驻点热流密度、动压和过载约束分别为qmax＝200kpa，nmax＝3g，攻角和倾侧角的幅值约束分别为α∈[10°,20°]，σ∈[-80°,80°],所得侧向可行走廊和初始侧向剖面如图3所示。
[0126] 本实施例，步骤3)的详细步骤包括：
[0127] 3.1)记Lz(e)和h(e)分别表示初始化侧向剖面和参考高度剖面，确定函数如式(3-1)所示；
[0128]
[0129] 式(3-1)中，D(e)表示阻力加速度剖面，Lz(e)表示侧向剖面，L表示升力速度，υ表示倾侧角。将气动加速度中与气动系数无关的项记为Aconst且函数如式(3-2)所示；
[0130]
[0131] 式(3-2)中，m表示飞行器质量，q表示动压，Sref表示参考面积；
[0132] 3.2)将式(3-2)代入式(3-1)化解
整理得式(3-3)：因为已经将飞行器的气动参数拟合为关于攻角和马赫数的二元函数，因此当高度、速度已知时可令升力系数CL、阻力系数CD的函数表达式如式(3-4)所示；
[0133]
[0134] 式(3-3)中，CL表示升力系数，CD表示阻力系数，Lz表示侧向升阻比，υ表示倾侧角，L表示升力加速度；
[0135]
[0136] 式(3-4)中，CL表示升力系数，CD表示阻力系数，Cli,Cdi,i＝1,2,3为与攻角无关的飞行器气动系数；
[0137] 3.3)将式(3-4)代入式(3-3)化解整理得式(3-5)和式(3-6)：
[0138]
[0139]
[0140] 式(3-5)和式(3-6)中，CL表示升力系数，CD表示阻力系数，Cli,Cdi,i＝1,2,3为与攻角无关的飞行器气动系数，Lz表示侧向升阻比，fα1～fα5为中间变量；求解(3-5)的一元四次方程将产生四个不同的根,根据实际物理意义和弹道特性，飞行器在给定的攻角变化区间内最多存在两个可行的实根,因此其余两个根一般为复数或极大和极小的数；同时为了确保阻力加
速度曲线的平滑过渡，当已知前一时刻的值后，根据(3-5)求解得到的两个可行攻角就只能有一个是可行解；当攻角αr确定后，结合设计的高度曲线可以唯一得到此时对应的升阻比 和阻力加速度D*，从而可得参考倾侧角υr如式(3-7)，可反求出此时对应的速度倾角θ*如式(3-8)，
[0141]
[0142] 式(3-7)中，Lz(e)表示e时刻的侧向加速度剖面值， 为结合设计的高度曲线可以唯一得到此时对应的升阻比；
[0143] θ＝-arcsin(Dh′) (3-8)
[0144] 式(3-8)中，θ表示结合设计的高度曲线可以唯一得到此时对应的速度倾角θ*，D表示阻力加速度，h′表示高度关于能量的导数；
[0145] 3.4)至此纵向运动状态r,V以及θ和控制量α和υ已唯一确定，其中，r表示地心距，V表示飞行器速度，θ表示当地速度倾角，α表示攻角，υ表示倾侧角，然后通过数值积分降阶运动模型求解余下的航向角σ、经度λ以及纬度φ等3个侧向运动状态，从而得到设计剖面对应的所有参考弹道状态。
[0146] 如图2所示，本实施例步骤4)的详细步骤包括：
[0147] 4.1)输入强高度剖面约束、复杂飞行任务，提取需要参数信息，包括：参考高度剖面、初始终端状态条件以及过程约束；
[0148] 4.2)设计侧向剖面的初值生成器，对侧向剖面参数进行初始化产生初始侧向剖面；
[0149] 4.3)根据给定的高度剖面约束和初始侧向剖面，求解纵向剖面获得参考三维剖面；
[0150] 4.4)判定解算得到的参考三维剖面是否满足给定的约束要求，满足则跳转执行下一步；否则，检查约束条件和联系方程是否合理，调整侧向剖面或相应的约束条件重新跳转执行步骤4.2)直至产生的三维剖面满足约束要求；
[0151] 4.5)建立融入航路点、禁飞区以及目标点的混合性能指标，联立三维剖面和相关约束条件形成非线性弹道规划的数学求解模型；
[0152] 4.6)利用非线性序列二次规划算法求解非线性弹道规划的数学求解模型，并输出求解得到的参考弹道；
[0153] 4.7)判断求解得到的参考弹道是否满足任务要求，若满足任务要求，则输出参考弹道；否则通过非线性规划求解器自动调整剖面参数，返回4.2)继续进行迭代计算，直到求解得到的参考弹道满足任务要求。
[0154] 本实施例中，步骤4.2)的详细步骤包括：
[0155] 4.2.1)根据飞行任务，判定初终端符号；根据确定的初终端符号，建立相应的侧向和纵向走廊边界；
[0156] 4.2.2)假定从滑翔起点到终点间的飞行路径为理想大圆弧，则根据式(4-1)所示球面几何三角变换可求得需要的总飞行航程Stogo；
[0157] stogo＝cos-1(sinφfsinφ+cosφfcosφcos(λf-λ)) (4-1)
[0158] 式(4-1)中，φf表示目标点地心纬度，φ表示飞行器当前地心纬度，λ表示飞行器当前经度，λf表示目标点经度；
[0159] 4.2.3)记阻力加速度走廊的最大最小边界分别为Dconmax和Dconmin，根据式(4-2)确定给定初始阻力加速度剖面权重系数ωD时对应的阻力加速度剖面D(e)；
[0160] D(e)＝ωD(Dconmax(e)-Dconmin(e))+Dconmin(e) (4-2)
[0161] 式(4-2)中，D(e)为给定初始阻力加速度剖面权重系数ωD时对应的阻力加速度剖面，为，Dconmax(e)和Dconmin(e)分别为加速度走廊的最大最小边界，ωD为给定初始阻力加速度剖面权重系数；
[0162] 4.2.4)为了快速求出式(4-2)对应的剖面航程，将式(4-2)进行在积分区间内离散，并假设各个剖面点之间采用线性连接，求得此时相邻两点间对应的阻力加速度剖面系数如式(4-3)所示；
[0163]
[0164] 式(4-3)中，aD(i)表示第i点到第(i+1)点间直线连接的阻力加速度剖面斜率，bD(i)表示第i点到第(i+1)点间直线连接的阻力加速度剖面常数项，D(i+1)表示第(i+1)点的阻力加速度，D(i)第i点的表示阻力加速度，e(i+1)表示第(i+1)点的归一化能量值，e(i)表示第i点的归一化能量值；
[0165] 4.2.5)记各阻力加速度点间对应的剖面航程为SD(i)，则可求得剖面对应的总航程ST如式(4-4)所示；
[0166]
[0167] 式(4-4)中，Nint表示离散点数目，SD(i)表示从第i点到第(i+1)点对应的剖面航程，aD(i)表示第i点到第(i+1)点间直线连接的阻力加速度剖面斜率，bD(i)表示第i点到第(i+1)点间直线连接的阻力加速度剖面常数项，e(i+1)表示第(i+1)点的归一化能量值，e(i)表示第i点的归一化能量值；
[0168] 4.2.6)迭代给定初始阻力加速度剖面权重系数ωD，当剖面对应的总航程ST、总飞行航程Stogo成立时即找到初始阻力加速度剖面；
[0169] 4.2.7)利用侧向走廊大小边界，通过插值获得对应的侧向剖面大小；根据初始判定的初终端符号对侧向剖面的符号进行调整，初始正、终端负时将求得的侧向指令大小剖面最小值点为边界，且按照初始正、终端负确定侧向剖面的符号；然后再从找到的初始剖面中反求出需要的初始迭代点，并进行归一化，从而得到需要的初值迭代点，产生初始侧向剖面。
[0170] 本实施例中，步骤4.5)的详细步骤包括：将侧向升阻比剖面离散化为一系列参考点，各参考点间进行线性连接产生标准侧向指令剖面；当确定纵、侧向剖面后，将其代入侧向运动模型中积分求解余下的侧向运动状态量；设定目标函数为与终端目标点的距离最小，则可建立如式(4-5)所示的非线性弹道规划的数学求解模型；
[0171] min J＝min(|λ(ef)-λf|+|φ(ef)-φf|)
[0172]
[0173] 式(4-5)中，J表示性能指标函数，λ(ef)和φ(ef)分别表示终端归一化能量ef对应的经度和纬度，λf和φf分别表示目标点经度和纬度；h(e)表示高度剖面，fh(e)表示给定的强高度剖面函数，L表示升力加速度，υ表示倾侧角，g为重力加速度，V是飞行器速度，r为地心距， 表示弹道状态量的导数，x(e)表示弹道状态量，u(e)表示控制变量，u(e)＝[D(e),Lz(e)]T，e表示归一化能量，f(x(e),u(e),e)表示用归一化能量表示的运动状态方程，e0表示初始归一化能量，ef表示终端归一化能量， 表示能量e时刻的驻点热流密度，q(e)表示能量e时刻的动压，L(e)表示能量e时刻的升力加速度，D(e)表示能量e时刻的阻力加速度，g0为海平面重力加速度， qmax以及nmax分别为最大热流密度、最大动压及最大过载。
[0174] 本实施例中将侧向升阻比剖面离散化为一系列参考点，各参考点间进行线性连接产生标准侧向指令剖面；当确定纵、侧向剖面后，将其代入侧向运动模型中积分求解余下的侧向运动状态量；值得注意的是，此时的高度已经已知，相应的速度和速度倾角也都可以唯一确定。因此，积分侧向运动模型产生的误差，理论上只有忽略地球自转产生的。由此，设定目标函数为与终端目标点的距离最小，则可建立如式(4-5)所示的最优控制问题模型。
[0175] 当考虑航路点和禁飞区任务时，记航路点的坐标为(λw,φw,ew)，禁飞区中心坐标为(λN,φN,RN)，则航路点约束可表示为式(4-6)：
[0176] JW＝|λ(ew)-λw|+|φ(ew)-φw|<εw (4-6)
[0177] 式(4-6)中，JW表示考虑航路点约束的性能指标，(λw,φw,ew)为记航路点的坐标，λ(ew)表示归一化能量ew时刻的经度，φ(ew)表示归一化能量ew时刻的纬度，φw表示给定的航路点纬度，λw表示给定的航路点经度，εw表示过航路点的精度要求。
[0178] 为了快速确定优化得到的轨迹是否通过航路点，可以将式(4-6)
修改为式(4-7)：
[0179] JW＝|φ(λw)-φw|<εw (4-7)
[0180] 式(4-7)中各参量与式(4-6)中相同。即将同一经度下两者的纬度差作为飞行轨迹与航路点的距离差。
[0181] 对于禁飞区，将禁飞区约束转换为航路点约束进行弹道优化设计。因此，考虑航路点和禁飞区约束的总性能指标可表示为式(4-8)：
[0182] J＝|λ(ef)-λf|+|φ(ef)-φf|+ωWJW+ωNJN (4-8)
[0183] 式(4-8)中，J表示考虑航路点和禁飞区约束后的复合性能指标，λ(ef)表示终端归一化能量ef对应的经度，ef表示终端归一化能量，λf表示设定的终端经度，φ(ef)表示终端归一化能量ef对应的纬度，φf表示设定的终端纬度，JW表示航路点约束对应的性能指标，JN表示禁飞区约束对应的性能指标，ωW和ωN分别表示航路点约束和禁飞区约束的惩罚系数，两者值越大，说明对其约束要求越严格。
[0184] 为检验本实施例方法的有效性，本实施例中还进一步设计了三维跟踪制导律。
[0185] 三维跟踪制导律根据采用纵侧向分别设计
控制器进行跟踪控制。因为纵向升力加速度分量Lcosυ和侧向加速度分量Lsinυ完整地包含了制导所需的控制量攻角和倾侧角，同时这两个方向的加速度分量容易通过飞行器自身携带的加表敏感得到。因此分别选择其作为纵侧向中间制导指令。在纵向，有高度动力学方程如式(5-1)：
[0186]
[0187] 式(5-1)中，表示高度对时间的导数，V表示飞行器速度，θ表示速度倾角。
[0188] 对方程两边继续求关于时间的二阶导数，同时将忽略地球自转后的动力学方程代入，整理化简得式(5-2)：
[0189]
[0190] 式(5-2)中，表示高度对时间的二阶导数， 表示飞行器速度的导数，V表示飞行器速度，表示速度倾角的导数，θ表示速度倾角，D表示阻力加速度，g表示引力加速度，r表示地心距，υ表示倾侧角，L表示升力加速度。
[0191] 由此，可设计一个纵向的二阶PD控制律对参考高度进行跟踪如式(5-3)：
[0192]
[0193] 式(5-3)中，h表示高度，表示高度对时间的导数，表示高度对时间的二阶导数，href表示参考高度， 表示参考高度的导数， 表示参考高度的二阶导数，ξh,ωh分别为阻尼系数和振荡
频率。将式(5-2)代入式(5-3)，整理得式(5-4)：
[0194]
[0195] 式(5-4)中，各变量含义参见式(5-2)和式(5-3)。
[0196] 对于侧向，通过设计一阶PD控制器跟踪三维剖面解算得到的航向角剖面产生侧向制导指令如式(5-5)：
[0197]
[0198] 式(5-5)中，σ表示航向角， 表示航向角的导数，σref表示参考航向角， 表示参考航向角的导数，kσ表示调节系数，ω表示振荡频率。同样，将忽略地球自转后的动力学方程代入式(5-5)，整理化简后可得需要的侧向加速度为式(5-6)：
[0199]
[0200] 式(5-6)中，L表示升力加速度，υ表示倾侧角，V表示飞行器速度，φ表示飞行器当前地心纬度，θ表示速度倾角，r为地心距，σref表示参考航向角， 表示参考航向角的导数，kσ表示调节系数，ω表示振荡频率。
[0201] 结合式(5-4)和式(5-6)，可求得需要的总升力加速度L*为式(5-7)：
[0202]
[0203] 式(5-7)中，L表示升力加速度，υ表示倾侧角。
[0204] 根据升力加速度的定义式，可反求出此时对应的升力系数 为式(5-8)：
[0205]
[0206] 式(5-8)中，L*表示需要的总升力加速度，m表示飞行器质量，q表示动压，Sref表示参考面积。从而利用拟合的升力函数关系式或气动系数表可以反求出需要的攻角α*。相应的制导需要的倾侧角υ*为式(5-9)：
[0207]
[0208] 式(5-9)中，L表示升力加速度，υ表示倾侧角。
[0209] 值得注意的是，在利用式(5-4)和式(5-6)求解的中间制导指令时需要将设计的参考剖面转换为关于时间的函数关系。最终，通过前述设计的三维跟踪制导律，证明了本实施例方法基于三维剖面产生的弹道落点可行性和有效性。
[0210] 为进一步检验本实施例方法基于三维剖面产生的弹道落点可行性和有效性，本实施任意选取侧向走廊内一条参考剖面基准，通过跟踪其对应的所有纵向阻力加速度剖面弹道进行仿真实验。本实施例设置滑翔起点经纬度(0°,0°)，初始射向为正东方向，初始高度为50.73km，初始速度为6380m/s，要求交班高度为30km，交班速度为2500m/s。图3、图4、图5、图6分别为阻力加速度剖面、侧向剖面、参考高度剖面以及航向角曲线的跟踪与规划对比结果，其目的在于对提出的三维剖面弹道进行可行性验证。图7、图8分别为攻角和倾侧角曲线。图9、图10分别为航路点约束的优化轨迹和侧向剖面。图11、图12分别为禁飞区约束的优化轨迹和侧向剖面。图13、图14分别为两个航路点约束的优化轨迹和三维剖面。通过对比可以明显看出，采用本实施例的强约束复杂任务条件下的三维剖面优化设计方法获得的滑翔机动目标
覆盖范围相比传统方法得到了有效释放。
[0211] 综上所述，本实施例的强约束复杂任务条件下的三维剖面优化设计方法无需事先优化设计参考攻角和倾侧角，采用侧向优先设计策略，综合兼顾飞行器纵、侧向任务需求，通过建立高度剖面约束与三维剖面的映射关系，并构建融入航路点和禁飞区约束的综合性能指标，从而将复杂三维剖面问题转换为单剖面搜索问题进行求解。此外，设计了初值生成器，确保利用序列二次规划算法快速获得收敛解。与现有方法相比，该方法解决了攻角和倾侧角同时优化调整下的强约束复杂任务下弹道设计问题，充分发挥了飞行器的机
动能力，适应了广域大机动目标。
[0212] 此外，本实施例还提供一种强约束复杂任务条件下的三维剖面优化设计系统，包括计算机设备，该计算机设备被编程或配置以执行前述强约束复杂任务条件下的三维剖面优化设计方法的步骤，或该计算机设备的存储介质上存储有被编程或配置以执行前述强约束复杂任务条件下的三维剖面优化设计方法的计算机程序。
[0213] 此外，本实施例还提供一种计算机可读存储介质，该计算机可读存储介质上存储有被编程或配置以执行前述强约束复杂任务条件下的三维剖面优化设计方法的计算机程序。
[0214] 以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
