CN111176325A

CN111176325A - 一种吸气式高超声速无人机协调区域分析方法

Info

Publication number: CN111176325A
Application number: CN202010009031.9A
Authority: CN
Inventors: 王玉惠; 侯思远; 陈谋; 吴庆宪; 沈艺; 徐超; 李云鑫
Original assignee: Nanjing University of Aeronautics and Astronautics
Current assignee: Nanjing University of Aeronautics and Astronautics
Priority date: 2020-01-06
Filing date: 2020-01-06
Publication date: 2020-05-19
Anticipated expiration: 2040-01-06
Also published as: CN111176325B

Abstract

本发明公开了一种吸气式高超声速无人机协调区域分析方法，步骤：建立吸气式高超声速无人机的气动力矩模型；求取气动力矩之间的气动耦合矩阵；利用对角优势矩阵理论确定可忽略与不可忽略气动耦合的界限；建立吸气式高超声速无人机五自由度姿态运动模型，利用后向可达集的思想，结合水平集的求解方法对吸气式高超声速无人机的飞行安全区域进行求解；将气动耦合界限作为约束，对飞行安全区域中的状态点进行筛选，符合条件的状态点构成飞行协调区域。通过本发明得到的协调区域，该区域内无人机姿态三通道之间的气动耦合较小，变量之间基本协调，对保证无人机的安全协调飞行有着重要的意义。

Description

一种吸气式高超声速无人机协调区域分析方法

技术领域

本发明属于高超声速无人机领域，特别涉及了一种高超声速无人机协调区域分析方法。

背景技术

吸气式高超声速无人机因其极高的军事应用价值，成为各国技术发展的重点。同时，由于无人机在高超声速飞行时，复杂的气动特性使其通道之间存在耦合，并在飞行过程中不断累积，影响飞行安全。因此，开展高超声速无人机耦合的研究，对其飞行协调的区域进行限定对保证飞行安全具有重要的意义。

近年来，关于高超声速无人机耦合方面的研究越来越多。Zhen Wubin，针对高超声速飞行器的纵向动力学模型做出了相应的耦合分析，并基于求解的耦合关系设计了滑模控制器。在飞行器众多耦合中由气动力与气动力矩引起的气动交叉耦合最为明显，针对气动耦合的研究也已经取得很多成果。Lin Yang等人对大迎角导弹控制系统进行深入研究，从气动参数和稳定裕度出发，论述了气动耦合大迎角导弹控制系统的稳定性会产生不利影响，指出解耦的必要性。Wang bin利用单变量控制系统理论的方法，对于气动耦合的解耦条件进行限定，实现飞行器气动耦合模型的解耦。Zhang Keke对“十”字形布局飞行器进行了详细的气动耦合进行分析，并采用BP神经网络对其进行耦合补偿，达到良好的控制效果。

可见，针对高超声速无人机耦合方面的研究成果显著，这都为后续的研究提供了重要依据。然而，这些研究成果都存在着一定的缺陷，对飞行器气动交联耦合的分析工作更是少之又少，且缺少定量的分析。同时也没有一个确定的协调区域的概念，不利于高超声速飞行的飞行安全。

发明内容

为了解决上述背景技术提到的技术问题，本发明提出了一种吸气式高超声速无人机协调区域分析方法。

为了实现上述技术目的，本发明的技术方案为：

一种吸气式高超声速无人机协调区域分析方法，包括以下步骤：

(1)建立吸气式高超声速无人机的气动力矩模型；

(2)对气动力矩之间的气动耦合进行定量的分析，求取气动力矩之间的气动耦合矩阵；

(3)基于步骤(2)中的气动耦合矩阵，利用对角优势矩阵理论确定可忽略与不可忽略气动耦合的界限，为之后的协调区域的划分奠定基础；

(4)建立吸气式高超声速无人机五自由度姿态运动模型，利用后向可达集的思想，结合水平集的求解方法对吸气式高超声速无人机的飞行安全区域进行求解；

(5)将步骤(3)中的气动耦合界限作为约束，对步骤(4)中的飞行安全区域中的状态点进行筛选，符合条件的状态点即构成飞行协调区域。

进一步地，在步骤(1)中，吸气式高超声速无人机的气动力矩模型如下：

上式中，l为滚转力矩，m为俯仰力矩，n为偏航力矩；

为飞行动压；S为无人机机翼的参考面积；C_l为滚转力矩系数；C_m为俯仰力矩系数；C_n为偏航力矩系数。

进一步地，所述力矩系数C_l、C_m、C_n的形式如下：

上式中，β为侧滑角；V代表飞行速度；b为机翼的翼展长度；c为平均气动弦长；p,q,r分别表示无人机的滚转角速率、俯仰角速率和偏航角速率；δ_e,δ_a,δ_r分别表示飞行器左、右升降副翼舵和方向舵；C_l,β为基本滚转力矩系数；C_m,α为基本俯仰力矩系数；C_n,β为基本偏航力矩系数；

分别表示左、右升降副翼舵和方向舵引起的滚转力矩增量系数；C_l,q,C_l,r分别表示由俯仰角速率与偏航角速率导致的滚转力矩增量系数；

分别表示由左、右升降副翼舵和方向舵导致的俯仰力矩增量系数；C_m,q是由俯仰角速率引起的俯仰力矩增量系数；

分别表示由左、右升降副翼舵和方向舵引起的偏航力矩增量系数；C_n,q,C_n,r分别表示由俯仰角速率和偏航角速率引起的偏航增量系数。

进一步地，在步骤(2)中，所述气动耦合矩阵如下：

上式中，

为姿态角耦合度矩阵，

为姿态角速率耦合度矩阵，

为等效舵面耦合度矩阵；F描述的是吸气式高超声速飞行器的气动力矩系数模型；Ω＝[α,β,μ]^T表示无人机飞行的姿态角，α,β,μ分别表示飞行迎角、侧滑角和滚转角；ω＝[p,q,r]^T表示无人机姿态角速率，p,q,r分别表示滚转角速率、俯仰角速率与偏航角速率；

为三通道的等效舵面，

δ_γ,δ_ψ分别为俯仰舵、偏航舵和滚转舵；上标T表示转置。

进一步地，在步骤(3)中，对步骤(2)得到的气动耦合矩阵进行主对角元素优势处理，具体步骤如下：

(3-1)对姿态角耦合度矩阵

中的每个元素都除以对应行的主对角元素：

上式中，

为处理后的矩阵

中的第i行第j列的元素，a_ij为处理前的矩阵

中的第i行第j列的元素，a_ii为对应行的主对角元素；处理后得到矩阵

同理，对姿态角速率耦合度矩阵

和等效舵面耦合度矩阵

按照上述方法进行处理，得到矩阵

和

(3-2)定义耦合强弱的界限k_i：

上式中，k_i表示矩阵

第i行的非对角元素之和，将k_i≤1的耦合视为可忽略的气动耦合，k_i＞1的耦合视为不可忽略的气动耦合；

同理，矩阵

和矩阵

按照上述方法定义其每一行的k_i；

(3-3)取一飞行状态，保持飞行高度不变，仅改变迎角与飞行速度，分别研究矩阵

中k_i的变化。

进一步地，在步骤(4)中，所述吸气式高超声速无人机五自由度姿态运动模型：

上式中，α为迎角；β为侧滑角；p为滚转角速率；q为俯仰角速率；r为偏航角速率；V是无人机飞行速度；M为无人机质量；L为升力；Y为侧力；l为滚转力矩；m为俯仰力矩；n为偏航力矩；I_xx为机体轴x方向的转动惯量；I_yy机体轴y方向的转动惯量；I_yy为机体轴z方向的转动惯量；上方一点表示微分。

进一步地，在步骤(4)中，吸气式高超声速无人机的飞行安全区域的求解过程如下：

(4-1)取一平衡状态，以俯仰舵为控制变量，利用连续算法求取吸气式高超声速无人机关于迎角α和侧滑角β的平衡面；

(4-2)取一平衡状态，保持偏航舵与滚转舵不变，仅改变俯仰舵的偏转角，研究(4-1)中得到的平衡面上迎角与俯仰舵的关系；

(4-3)取一平衡状态，保持偏航舵与滚转舵不变，仅改变俯仰舵的偏转角，研究(4-1)中得到的平衡面上侧滑角与俯仰舵的关系；

(4-4)根据(4-1)求得的平衡面作为目标集，利用水平集思想对后向可达集求取，将此后向可达集作为吸气式高超声速无人机的飞行安全区域。

进一步地，在步骤(5)中，对于步骤(4)求出的飞行安全区域，利用步骤(3)中得到的可忽略与不可忽略气动耦合的界限进行筛选，将不可忽略气动耦合的状态点剔除，剩余的状态点满足气动耦合要求，这些状态点即构成了飞行器的飞行协调区域。

采用上述技术方案带来的有益效果：

(1)本发明所提出的利用对角优势理论的气动耦合矩阵处理方法，比以往的耦合分析方法更具动态性和整体性，可以反应出通道之间的整体耦合，而不局限于某些变量之间的耦合，更具有实际应用价值；

(2)本发明将后向可达集与气动耦合的约束相结合，求取高超声速无人机的协调区域既满足飞行的稳定性又保证了协调性，从多方面保证了飞行器的飞行安全；

(3)本发明所研究的无人机结构，为国际上最为经典的乘波体构型，因此研究结果更加符合实际应用，具有良好的前景。

附图说明

图1是本发明的方法流程图；

图2是吸气式高超声速无人机姿态角耦合度矩阵的k值变化图；其中，(a)、(b)、(c)依次代表俯仰通道、偏航通道和滚转通道；

图3是吸气式高超声速无人机姿态角角速率耦合度矩阵的k变换图；其中，(a)、(b)、(c)依次代表俯仰通道、偏航通道和滚转通道；

图4是吸气式高超声速无人机等效舵面耦合度矩阵的k变换图；其中，(a)、(b)、(c)依次代表俯仰通道、偏航通道和滚转通道；

图5是吸气式高超声速无人机俯仰舵面与迎角α平衡状态图；

图6是吸气式高超声速无人机俯仰舵面与侧滑角β平衡状态图；

图7是吸气式高超声速无人机迎角α与侧滑角β可达集状态图；

图8是吸气式高超声速无人机对安全飞行区域进行气动耦合约束筛选流程图；

图9是吸气式高超声速无人机经过气动耦合约束筛选的迎角α与侧滑角β协调集状态图。

具体实施方式

以下将结合附图，对本发明的技术方案进行详细说明。

本发明设计了一种吸气式高超声速无人机协调区域分析方法，如图1所示，步骤如下：

步骤1：建立吸气式高超声速无人机的气动力矩模型；

步骤2：对气动力矩之间的气动耦合进行定量的分析，求取气动力矩之间的气动耦合矩阵；

步骤3：基于步骤2中的气动耦合矩阵，利用对角优势矩阵理论确定可忽略与不可忽略气动耦合的界限，为之后的协调区域的划分奠定基础；

步骤4：建立吸气式高超声速无人机五自由度姿态运动模型，利用后向可达集的思想，结合水平集的求解方法对吸气式高超声速无人机的飞行安全区域进行求解；

步骤5：将步骤3中的气动耦合界限作为约束，对步骤4中的飞行安全区域中的状态点进行筛选，符合条件的状态点即构成飞行协调区域。

在本实施例中，上述步骤1采用如下优选方案实现：

本发明研究的吸气式高超声速无人机姿态三通道对应的气动力矩为滚转力矩l，俯仰力矩m与偏航力矩n，l,m,n的具体形式可以表示为：

其中，

为飞行动压；S为无人机机翼的参考面积；C_l为滚转力矩系数；C_m为俯仰力矩系数；C_n为偏航力矩系数。C_l,C_m,C_n的具体形式可以写为：

其中，α为飞行迎角；β为侧滑角；V代表飞行速度；b为机翼的翼展长度；c为平均气动弦长；p,q,r分别表示无人机的滚转角速率、俯仰角速率和偏航角速率；δ_e,δ_a,δ_r分别表示飞行器左、右升降副翼舵和方向舵；C_l,β为基本滚转力矩系数；C_m,α为基本俯仰力矩系数；C_n,β为基本偏航力矩系数；

分别表示由左、右升降副翼舵和方向舵引起的偏航力矩增量系数；C_n,q,C_n,r分别表示由俯仰角速率以及偏航角速率引起的偏航增量系数。涉及到的迎角α单位为度(°)，侧滑角α单位为弧度(rad)。

当吸气式高超声速无人机进行高超声速飞行时，由于其工作环境恶劣，这使得飞行器的空气动力学特性变的十分复杂，状态变量之间会出现强烈的非线性关系，从而导致俯仰、偏航与滚转通道之间的气动力矩出现相互耦合的情况，这就是高超声速无人机中的气动耦合。气动耦合的出现对于飞行器的飞行安全十分不利，特别是存在侧滑角和滚转角的情况下，飞行过程中极易出现不对称的随随机涡流现象，从而导致偏航通道与滚转通道上力矩的出现，如果对次效应不加以控制，该力矩会不断累积，可能会使得飞行器失去对偏航通道与滚转通道的控制，严重影响飞行器的飞行安全。因此对吸气式高超声速飞行器姿态系统中存在的气动耦合进行处理是十分必要的。

在本实施例中，上述步骤2采用如下优选方案实现：

对各通道之间的气动耦合进行定量的分析，借鉴动态方程的耦合分析方法，求得各气动力矩之间的耦合度矩阵，具体包括姿态角耦合度矩阵、姿态角速率耦合度矩阵以及等效舵面耦合度矩阵，从而实现对各通道之间存在的气动耦合的定量表示，为之后可忽略/不可忽略耦合界限的提出奠定基础。

吸气式高超声速无人机飞行器的气动力矩系数可以表示为C＝[C_l,C_m,C_n]^T，姿态角可以表示为Ω＝[α,β,μ]^T，姿态角速率可以表示为ω＝[p,q,r]^T，舵偏角可以表示为δ＝[δ_e,δ_a,δ_r]^T。由步骤1可知，气动力矩系数可以进一步写为关于姿态角、姿态角速率以及舵偏角的函数，具体可以表示为如下形式：

C＝F(Ω,ω,δ)

高超声速飞行器模型中的原始左、右升降副翼舵和方向舵并不能完全对应为俯仰、偏航与滚转通道的控制舵面，因此为了得到三通道各自对应的等效舵面，需要对操纵舵面的进行等效偏转策略的设计。具体策略是左、右升降舵全动提供高超声速无人机的俯仰力矩，依旧由方向舵提供偏航力矩，而用左、右升降舵差动来提供滚转力矩，因此可以得到三通道对应的等效舵面

即俯仰舵

偏航舵δ_γ和滚转舵δ_ψ，具体形式可写为：

δ_γ＝δ_r

δ_ψ＝(δ_e-δ_a)

因此，可以进一步将姿态三通道的气动力矩系数表达式写为:

C＝F(Ω,ω,δ_eq)

此时，俯仰通道对应的系统变量为迎角α，俯仰角速率q和俯仰舵

偏航通道对应的系统变量为侧滑角β，偏航角速率r以及偏航舵δ_γ；滚转通道对应的系统变量为滚转角μ，滚转角速率p和滚转舵δ_ψ。则若各通道的气动力矩只受到对应的变量影响，则三通道气动力矩之间相互独立，不存在耦合。但是由步骤1中表示的吸气式高超声速无人机气动力矩系数的模型明显可以看出，三通道气动力矩系数中都存在着非对应通道的变量，因此俯仰、偏航和滚转通道之间相互存在着不同程度的耦合。为了得到定量的耦合程度，借鉴动态方程的耦合分析方法，对三通道之间的耦合进行计算。具体过程为：

将气动力矩系数C分别对Ω,ω,δ_eq求取偏导数，可以得到以下三个Jacobian矩阵，即

具体表示为如下形式：

由上述三个的Jacobian矩阵，结合步骤1中的气动力矩模型，可将三通道气动力矩之间姿态角耦合度矩阵定义为

姿态角速度耦合度矩阵定义为

等效舵面耦合度矩阵定义为

具体可写为如下形式：

由Jacobian的含义知道，矩阵中的每一个元素都代表自变量与对应因变量之间的影响程度。因此

中每个元素表示对应的系统变量对该行通道气动力矩的影响大小，这与耦合的定义是一致的，因此这样定义耦合度矩阵是合理的。

若三通道气动力矩之间相互独立，互不影响，

都应为对角矩阵，只受到对应通道元素的影响，但由步骤1中描述的气动力矩系数模型易知，气动耦合力矩系数之间交叉含有非对应通道的元素，即矩阵

中存在着不为0的非对角元素，因此对三通道气动力矩之间存在着气动耦合。

本发明通过对吸气式高超声速无人机的气动力矩进行定性的分析，可以更加准确的表示出通道之间的相互耦合关系，对于无人机的协调飞行具有重要意义，同时该种耦合分析方法，同样适应于其他的强耦合系统。

在本实施例中，上述步骤3的优选实施方案如下：

步骤3-1)，对步骤2中得到的姿态耦合度矩阵

姿态角速率耦合度矩阵

以及等效舵面耦合度矩阵

中的每个元素的大小表示对应的自变量对该行通道气动力矩系数的影响程度。但此时的元素只表示每个变量的绝对影响的大小，并不能真正的反应出通道之间整体的耦合情况，因此需要对矩阵做进一步的主对角处理，其具体处理过程可以用矩阵

为例：

对姿态角耦合度矩阵

中的每个元素都除以对应行的主对角元素，具体处理方式如下：

矩阵

经过以上处理之后可以得到新的耦合度矩阵

此时

中的元素表示对应变量对该行通道气动力矩系数的相对影响大小，同样的对姿态角速率耦合度矩阵

和等效舵面耦合度矩阵

做主对角处理之后可以得到矩阵

和

的具体形式可写为：

步骤3-2)，由步骤3-1)中的分析易知，此时的矩阵

与

中各行的主对角元素都为1。以矩阵

为例，定义耦合强弱的界限k_i，具体可以写为：

实际上此时的k_i表示姿态角耦合度矩阵

第i行的非对角元素之和，按照主对角占优的理论，k_i的数值越大，说明该通道受到其他通道的影响就越大，耦合就越严重。当每行的k_i≤1时，可认为是该通道对应的元素对本通道影响占优势，反之该通道受其他通道元素的影响最大，视为不协调。则此时可将k_i≤1的耦合视为可忽略耦合，k_i＞1的耦合视为不可忽略的耦合。

根据步骤1中吸气式高超声速无人机气动力矩模型，选定此时的一飞行状态飞行高度H＝30000m,侧滑角β＝0°/s，仅改变迎角α与飞行速度V，研究姿态角度耦合度矩阵

中k_i的变化，由图2中的(a)可以看出对于姿态角耦合度矩阵

中俯仰通道对应的k值随着迎角α的增大会呈现出先减小后增大的趋势，不难看出此时k值都较大，存在着k＞1的情况；由图2中的(b)可以看出姿态角耦合度矩阵

中偏航通道对应的k值先减小后增大，且增大的速度比减小的迅速，在迎角α＝0°时取得最小值。由图2中的(c)可以看出，当飞行速度一定时，姿态角耦合度矩阵中滚转通道对应的k值基本在某个定值附近来回波动，但同样存在k＞1的情况。

步骤3-3)，根据步骤1中吸气式高超声速无人机气动力矩模型，选定此时的一飞行状态飞行高度H＝30000m,侧滑角β＝0°/s，仅改变迎角α与飞行速度V，研究姿态角速率耦合度矩阵

中k_i的变化。由图3中的(a)与(b)可以观察出，姿态角速率耦合度矩阵中俯仰通道和偏航通道对应k值存在这大于1的情况，而图3中的(c)描述的姿态角速率耦合度矩阵中滚转通道对应k小于1。

步骤3-4)，根据步骤1中吸气式高超声速无人机气动力矩模型，选定此时的一飞行状态飞行高度H＝30000m,侧滑角β＝0°/s，仅改变迎角α与飞行速度V，研究等效舵面耦合度矩阵

中k_i的变化。由图4中的(a)可以看出等效舵面耦合度矩阵

中俯仰通道对应的k值会随着迎角与飞行速度的增大而呈现上升的趋势，但是始终小于1；由图4中的(b)可以看出，矩阵

中偏航通道对应的k值会随着迎角与飞行速度的增大而减小，同样的k的数值始终小于1；由图4中的(c)可以看出，随着迎角的增大，等效舵面耦合度矩阵中对应滚转通道的k值会呈现先减小后增大的趋势，但是始终保持在1之下。

综合上述分析可知，姿态角耦合度矩阵中的k_i数值始终都大于1，为不可忽略的气动耦合。而姿态角速率耦合度矩阵和等效舵面耦合度矩阵中的k_i数值始终都小于1，为可忽略耦合。因此本发明将姿态角耦合度中的k_i作为步骤5中吸气式高超声速无人机气动耦合的约束，从而使得无人机在保持性能的同时可以保持较大的工作范围，使无人机具有良好的稳定性。

在本实施例中，上述步骤4的优选实施方案如下：

步骤4-1)，建立高超声速无人机姿态五自由度姿态模型，所述模型如下：

其中，V是无人机飞行速度；M为无人机质量；α为迎角；β为侧滑角；p为滚转角速率；q为俯仰角速率；r为偏航角速率；L为升力；Y为侧力；l为滚转力矩；m为俯仰力矩；n为偏航力矩；I_xx为机体轴x方向的转动惯量；I_yy机体轴y方向的转动惯量；I_yy为机体轴z方向的转动惯量。

对于上述的高超声速飞行器姿态模型，求解α,β可达集之前需要对其状态进行约束作为目标集，目标集可以由高超声速飞行器的平衡状态点表示。当高超声速飞行器处于较大攻角、与大马赫数飞行状态时，机身对于方向舵会产生一定的遮蔽作用，从而航向的稳定性受到一定的影响，会产生较大的激波阻力，进而导致方向舵的偏转效率的下降，当由左、右升降副翼舵联动等效的俯仰舵面偏转时，会引起阻力的变化进而导致侧滑角的出现，随后会导致滚转力矩的改变，飞行过程中不断的累积，最后甚至使得滚转力矩的数值大于左、右升降副翼舵面产生的力矩，因此出现舵面反操作的现象，进而会引起攻角的变化，这也将影响到两个机翼所受到的阻力，产生阻力差。由于机翼间阻力差的存在进而又会产生侧滑角，在相对较短的时间内飞行器发生攻角与侧滑角交替式强耦合效应，此耦合会使得侧滑角的不断加大，经过长时间的耦合积累严重影响着飞行器的姿态运动。因此在高超声速飞行器姿态运动过程中，保证α,β的协调是十分必要的。

步骤4-2)，对于步骤4-1)中表示的五自由度姿态非线性系统，令x＝[α,β,p,q,r]^T，

则姿态系统可表示为：

取一平衡点α＝4.372°，β＝0.235°，p＝0°/s，q＝0°/s，r＝0°/s为初始点x₀，此时俯仰舵

记为u₀。则此时

然后取u₀附近的一点为u₁，通过牛顿迭代法求得下一个平衡点x₁。

由

Df(x_0(k))(x_0(k+1)-x_0(k))+f(x_0(k))＝0

可得

x_0(k+1)＝x_0(k)-[Df(x_0(k))]^-1f(x_0(k))

其中

其中k＝0,1,2...为迭代的次数，当||x_0(k+1)-x_0(k)||＜10^-5时，新的平衡点为

x₁＝x_0(k+1)

接着可以求取平衡点x₁处的雅可比矩阵的特征值，然后根据特征值判断分叉点类型和稳定性情况，雅克比矩阵J为：

具体形式可以为：

步骤4-3)，由步骤4-2)中计算可得到从一个平衡点x₀,u₀出发，其中x₀具体为α＝4.372°，β＝0.235°，p＝0°/s，q＝0°/s，r＝0°/s，u₀具体为

经连续算法计算可得

的变化图，具体如图5所示，其中H为霍夫分叉点。

图5中实线代表稳定的平衡面(平衡点处的雅可比矩阵特征值全部分布在左半平面)，虚线代表不稳定的平衡面(平衡点处雅可比矩阵特征值不完全分布在左半平面)。由图5可知迎角α的平衡面大体分为三部分，其中当迎角位于5°左右回到了可控的区域范围，这是要主要研究的部分。当俯仰舵面为-11.5°时，会出现一个霍夫分叉点，平衡点的稳定性会发生改变，此时平衡点的雅可比矩阵的特征值会出现一对共轭复数根，其特征值λ为

经过该点之后，系统由不稳定的平衡面进入稳定的平衡面。同样的，当俯仰舵为12.5°时会出现另外一个霍夫分叉点，经过该点之后，由稳定的平衡面过渡到不稳定的平衡面。

步骤4-4)，由步骤4-2)中计算可得到从一个平衡点x₀,u₀出发，其中x₀具体为α＝4.372°，β＝0.235°，p＝0°/s，q＝0°/s，r＝0°/s，u₀具体为

经连续算法计算可得

的变化图，具体如图5所示，其中H为霍夫分叉点。

图6中实线代表稳定的平衡面(平衡点处的雅可比矩阵特征值全部分布在左半平面)，虚线代表不稳定的平衡面(平衡点处雅可比矩阵特征值不完全分布在左半平面)。由图6可以看出，对于侧滑角β会存在着许多交错的平衡面，但大多都是不稳定的，当偏航舵面位于[-11.5°,12.5°]之间会出现稳定的平衡面，该平衡面两端的霍夫分叉点附近可能会出现周期振荡，如出现机翼的摇晃等。

经连续算法计算出来的图5与图6描述的吸气式高超声速无人机的稳定平衡面上飞行状态是稳定的，因此可以将稳定的平衡面的状态点作为步骤4-5)中的目标集。

步骤4-5)根据步骤4-2)求得的迎角与侧滑角平衡面作为目标集，利用水平集思想对后向可达集求取。

状态变量构成的区域边界曲线形状在系统动态的映射函数f驱动下不断发生变化，在任意时刻满足水平集的方程:

其中，φ(x,t)为Lipschitz连续的水平集函数，为隐函数，有助于描述可达集的演化过程；

为φ(x,t)的梯度。设目标集J₀为有界开区域，其边界为

则目标集合的零水平集函数具有如下性质：

此时，目标集J₀可用其零水平集函数表示为：

J₀＝{x∈Rⁿ|φ(x,0)≤0}

要想求得目标集J₀在系统动态的映射函数f作用下的后向可达集，可以通过求解Hamilton-Jacobi得到。

Hamilton-Jacobi方程可以表示为：

其中，φ(x,0)＝φ(x),x∈Rⁿ，时间t是逆向的，Hamilton函数中H(x,p)的具体形式为：

上式中p∈Rⁿ，为哈密尔顿共态向量，具体表示为：

将使得H(x,p)取得最大的控制变量记作u^*，具体可写为：

将上式表示的最优控制变量u^*(x,p)代入式Hamilton-Jacobi方程求得目标集J₀的后向可达集P_τ(J₀)：

P_τ(J₀)＝{x∈Rⁿ|φ(x,-s)≤0}

其中s∈[0,τ]，s表示后向可达集到达目标集的时间。

Hamilton-Jacobi方程的数值计算只要包括三部分：空间导数D_xφ(x,t)，Hamilton函数H(x,p)与时间导数D_tφ(x,t)，其中

对于空间导数D_xφ(x,t)使用网络状态空间进行划分，采用定向逼近的方法求取。设划分的总格点数为n₀，第i(0＜i＜n₀)个格点对应的状态为x_i，则其左导数与右导数可以利用相邻的网络格点x_i-1,x_i+1得到具体表示为：

其中p^-、p⁺分别代表对D_xφ(x,t)导数的左、右逼近。

对于Hamilton函数使用Lax-Friedrichs格式进行逼近得到其近似值，具体为：

其中，h与Hamilton函数关于p的偏导数相关，对于状态空间的第i个维度的h_i可以写为如下形式：

其中，

与

分别代表D_xφ(x,t)在第i个维度上做左导数p^-和右导数p⁺的最小值与最大值，具体为

对于时间导数D_tφ(x,t)本发明采用收敛模式对后向可达集进行计算，从时间t＝0开始进行求解，直到H(x,p)≈0为止，此时的终止时间τ为收敛时间。

因为，此时考虑的主要是迎角与侧滑角之间的相互影响，即俯仰通道与偏航通道之间的气动耦合，则忽略滚转通道，约束滚转角速率p＝0°/s。由步骤4-1)描述的吸气式高超声速无人机五自由度姿态模型，此时的状态变量为α,β,q,r，控制变量为

则此时的Hamilton函数H(x,p)可以表示为：

其中，p₁,p₂,p₃,p₄分别为水平集函数φ(x,t)对于状态参数α,β,q,r的偏导数。

由上式可知，要想保证H(x,p)取得最大值只需要取得最优的控制变量

需要对H(x,p)求取

的一阶导数，可得

因此当

说明H(x,p)是关于

的增函数，此时取

当

说明H(x,p)是关于

的增函数，此时取

将最优的控制量

代入Hamilton函数，从从t＝0时对H(x,p)进行求解，直至H(x,p)≈0，由式(4.5)可知此时

意味着可达集停止增长，此时的解即为要求的后向可达集，其中α,β截面图如图7所示。

由图7可以观察出α,β的可达集是一个近似的三角区域，迎角α的最大值基本保持在8°左右，而此时的侧滑角β维持在0°上下。但是随着侧滑角β增大，迎角α达到的最大值呈现出整体下降的趋势，这种现象是因为α,β之间存在强烈的非线性耦合，这种特性会使得常规的气动舵面出现控制效能下降的现象，若α,β都出现过大的角度，此时需要很强的控制力矩去保持飞行器的飞行状态，一旦耦合累积超过一定的限度很可能出现舵面反转的现象，此时飞行器的飞行状态很难再控制回安全区域。

在本实施例中，上述步骤5的优选实施方案如下：

对于图7所表示的安全可达集区域上选取一状态点，H＝30000m，V＝3000m/s，α＝4°，β＝4°，p＝q＝r＝0°/s。由步骤2描述的耦合分析方法，对该状态的姿态角耦合度矩阵中关于α,β的部分为：

带入状态值上式可写为

由式上式可知姿态角耦合度矩阵中俯仰通道对应的k₁＝1.327＞1，偏航通道的对应的k₂＝1.024＞1。由以上分析可知，这说明俯仰通道受到侧滑角β影响强烈，偏航通道受到迎角α的影响剧烈，显然这种飞行状态存在着较强的气动耦合，属于飞行不协调的状态。因此对于图7描述的飞行区域要对不协调的飞行状态进行剔除，即将那些姿态角耦合度矩阵中存在k＞1情况的状态点剔除，剔除的具体流程如图8所示。

由图9所示，经过剔除之后，以上的区域均为飞行协调的状态，与图7相比可见，该区域子的最大迎角α为6°左右，而对于侧滑角β也只有一部分的状态超过了4°。大多的飞行状态维持在一个五边形的区域，明显小于图7区域，这也说明了飞行协调的要求比仅作出的安全性能要求要严格的多。

综上所述，本发明基于气动耦合的高超声速无人机协调区域如图9所示，在该区域内无人机姿态三通道之间的气动耦合较小，变量之间基本协调，对保证无人机的安全协调飞行有着重要的意义。

实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。