CN106020215A - 一种基于单步预测矩阵博弈的近距空战自主决策方法 - Google Patents

Abstract

一种基于单步预测矩阵博弈的近距空战自主决策方法，该方法的步骤如下：步骤1：搭建六自由度非线性无人作战飞机控制律结构；步骤2：初始化矩阵博弈棋局；步骤3：根据博弈棋局进行单步预测计算；步骤4：计算支付函数矩阵；步骤5：通过极大极小算法进行策略选择；步骤6：更新六自由度飞机运动学和动力学方程；步骤7：判断是否达到空战终止条件。本发明的优点是相比于三自由度质点模型更具有实际应用价值。同时，将现有的基于机动动作库的矩阵博弈方法更改为基于指令模型的机动库，只需要单步步长的预测，有效减小了决策时间，满足了空中验证的实时性需求，能更好地适应复杂动态的战场环境变化，提高了无人作战飞机在近距格斗中的作战能力。

Description

一种基于单步预测矩阵博弈的近距空战自主决策方法

【技术领域】

本发明是一种基于单步预测矩阵博弈的近距空战自主决策方法，属于无人机自主空战决策领域。

【背景技术】

随着无人作战飞机的性能和武器的提升，无人作战飞机充分利用信息技术革命时代的各种技术成果，进一步向更高性能和更高自主作战能力的全新信息化武器系统方向发展。

按照空战作战距离的分类，可以分为超视距空战和近距空战两个空战阶段，由于导弹的命中概率、敌机的电子干扰等问题，无人作战飞机在超视距攻击时可能会出现攻击失败的情况，此时双方就需要进入近距格斗作战阶段。作为现代空战的基本形式之一，近距格斗需要格斗双方根据自身状态和外部信息来进行机动决策，决策原则是相对于敌机获得有利的占位，使敌机进入我方飞机的攻击包络范围中。本发明旨在提高无人作战飞机的近距空战决策能力，使其可以实时地进行离散指令机动的选择，从而具备等同于有人战斗机上飞行员的决策能力。

目前，空战自主决策方法主要包括专家系统法、神经网络法、微分对策法和矩阵博弈等。其中矩阵博弈的核心思想就是对博弈树节点的估值过程与对博弈树搜索过程的结合，通过对博弈树搜索找出当前的最优策略，通过对博弈树进行极大极小的搜索可以达到这一目的。

在现有无人作战飞机空战自主决策研究中，矩阵博弈所采用的无人作战飞机模型均是较为简化的三自由度质点模型，而且其决策所需的待选机动库是根据空战战术动作而建立的，在进行博弈计算的过程中往往需要较大的计算量。然而，真实的无人作战飞机在空间中的运动还应当包括姿态角的变化，并且由于战场态势的快速改变，决策所产生的战术动作往往得不到有效执行就需要切换下一个动作。本发明即面向更为实际的六自由度非线性无人作战飞机模型，同时为满足战场态势的快速变化需求，提出了一种基于单步预测矩阵博弈的近距空战自主决策方法，以提高决策过程的实时性和有效性。

【发明内容】

1、发明目的：

本发明提出了一种基于单步预测矩阵博弈的近距空战自主决策方法，其目的是提供一种更具有实际应用价值的无人作战飞机空战自主决策方法，旨在保证决策正确性和科学性的同时，有效提高决策的实时性，以及对战场态势快速变化的适应能力，从而提高无人作战飞机的作战能力。

2、技术方案：

本发明针对采用六自由度非线性无人作战飞机模型，开发了一种基于单步预测矩阵博弈的近距空战自主决策方法，该方法的实现步骤如下：

步骤一：搭建六自由度非线性无人作战飞机控制律结构

(1)无人作战飞机配平

设置飞机的油门杆配平推力δ_Ttrim，升降舵配平偏角δ_etrim，副翼配平偏角δ_atrim，方向舵配平偏角δ_rtrim在给定的配平高度h_trim，配平迎角α_trim和飞行配平速度V_trim的条件下，使得无人作战飞机所受合力和合力矩为零，从而保持平飞状态，即确定飞机在配平状态下的参数。

(2)航迹倾斜角自动驾驶仪设计

在配平状态下，通过引入俯仰角速度q，飞机机体轴法向过载n_z，滚转角φ，航迹倾斜角μ的反馈信号，设计航迹倾斜角保持的自动驾驶仪，将控制指令以升降舵偏角指令δ_e的形式给到飞机模型，进而实现飞机的航迹倾斜角对航迹倾斜角指令信号μ_com的跟踪。

(3)自动油门系统设计

在配平状态下，为了在飞机爬升时对其动力进行补偿，通过引入飞行速度V的反馈信号，将飞行速度指令V_com。给到飞机的油门杆推力输入δ_T，从而实现飞行速度对飞行速度指令V_com的跟踪保持。

(4)滚转角自动驾驶仪设计

在配平状态下，通过引入滚转角速度p，偏航角速度r，滚转角φ，迎角α，和机体轴侧向过载n_y作为反馈信号，同时引入方向舵偏角指令δ_r，将滚转角的控制指令给到副翼偏角指令δ_a，设计滚转角保持的自动驾驶仪，从而实现对滚转角指令信号φ_com的跟踪。

步骤二：初始化矩阵博弈棋局

(1)初始化空战双方初始飞行状态

在本发明中，设定博弈的双方分别为红方和蓝方，对空战双方的初始飞行状态进行初始化设置，包括红方和蓝方，飞机的位置信息和姿态信息，即

X＝[x_g y_g h_g φ θ ψ V α β p q r]

其中，X为状态量的集合，x_g为飞机的北向坐标，y_g为飞机的东向坐标，h_g为飞机的高度坐标，φ为飞机的滚转角，θ为飞机的俯仰角，ψ为飞机的偏航角，V为飞行速度，α为飞机的迎角，β为飞机的侧滑角，p为滚转角速度，q为俯仰角速度，r为偏航角速度。

(2)初始化离散指令机动库

本发明中，为了满足自主决策实时性的需要，采取按照飞机机动指令离散信号的形式构建博弈所需的待选机动动作库的方法。认为空战过程中飞行速度V保持不变，同时飞行速度指令V_com保持不变，将纵向的航迹倾斜角指令信号μ_com和横侧向的滚转角指令信号φ_com进行组合，从而形成可以执行平飞、爬升、俯冲、协调转弯、爬升转弯和俯冲转弯等常见战术机动的机动指令信号。

本发明中，离散指令机动库为

Lib_RED ⁽ⁱ⁾＝[μ_RED ⁽ⁱ⁾ _com,φ_RED ⁽ⁱ⁾ _com] (1)

Lib_BLUE ^(j)＝[μ_BLUE ^(j) _com,φ_BLUE ^(j) _com] (2)

其中，Lib_RED为候选指令机动库集合，i表示第i个机动指令组合，μ_RED ⁽ⁱ⁾ _com表示第i组航迹倾斜角指令，φ_RED ⁽ⁱ⁾ _com表示第i组滚转角指令，Lib_BLUE为候选指令机动库集合，j表示第j个机动指令组合，μ_BLUE ^(j) _com表示第j组航迹倾斜角指令，φ_BLUE ^(j) _com表示第j组滚转角指令。

(3)初始化候选机动指令的博弈棋局

本发明中，博弈棋局表示在进行决策时，红蓝双方飞机的所有候选机动指令组合的集合，以矩阵的形式给出。认为红方一共有m种候选机动组合，蓝方一共有n种候选机动组合，则候选机动的矩阵为

步骤三：根据博弈棋局进行单步预测计算

(1)读入当前时刻的飞行状态

读取红蓝双方飞机当前时刻的飞行状态，以进行单步预测的计算，读取的飞行状态信息为

Y＝[x_g0 y_g0 h_g0 V₀ μ₀ χ₀]

其中，Y为状态量的集合，x_g0为飞机的初始北向坐标，y_g0为飞机的初始东向坐标，h₀为飞机的初始高度坐标，V₀为飞机的初始气流速度，μ₀为飞机的初始航迹倾斜角，χ₀为飞机的初始航迹方位角。

(2)读入矩阵棋局中的候选机动指令

读取步骤二中候选机动指令的博弈棋局，分别得到该矩阵中红蓝双方不同的机动指令组合，即

{Lib_BLUE,Lib_RED}＝{[μ_BLUEcom,φ_BLUEcom],[μ_REDcom,φ_REDcom]} (4)

其中，Lib_RED为候选指令机动库集合，Lib_BLUE为候选指令机动库集合，μ_REDcom表示红方的航迹倾斜角指令，φ_REDcom表示红方的滚转角指令，μ_BLUEcom表示蓝方航迹倾斜角指令，φ_BLUEcom表示蓝方滚转角指令。

(3)进行单步预测微分方程计算

考虑到飞机使用的是航迹倾斜角和滚转角的自动驾驶仪行驶结构以及速度保持的自动油门系统。本发明中，对无人机的飞行状态建立预测模型进行预测，建立飞机的预测航迹角变化的微分方程为

$\left\{\begin{array}{c}\stackrel{\·}{\mathrm{\μ}}=\frac{1}{T_{\mathrm{\μ}}}({\mathrm{\μ}}_{com}-\mathrm{\μ})\\ n_f=\frac{1}{{\mathrm{cos\φ}}_{com}}(\frac{V\stackrel{\·}{\mathrm{\μ}}}{g}+cos\mathrm{\μ})\\ \stackrel{\·}{\mathrm{\χ}}=\frac{g}{v\cos \mathrm{\μ}}{n}_f{\mathrm{sin\φ}}_{com}\end{array}\right.---\left(5\right)$

建立飞机的预测位置信息的微分方程为

$\left\{\begin{array}{c}{\stackrel{\·}{x}}_g=Vcos\mathrm{\μ}cos\mathrm{\χ}\\ {\stackrel{\·}{y}}_g=Vcos\mathrm{\μ}sin\mathrm{\χ}\\ {\stackrel{\·}{h}}_g=Vsin\mathrm{\μ}\end{array}\right.---\left(6\right)$

其中，μ为航迹倾斜角，为航迹倾斜角的微分，T_μ为航迹角响应的一阶等效模型的时间常数，μ_com为航迹倾斜角指令信号，n_f为飞机航迹坐标下的法向过载，φ_com为滚转角指令信号，g为重力加速度，为航迹方位角的微分，χ为航迹方位角，x_g为飞机的北向坐标，y_g为飞机的东向坐标，h_g为飞机的高度坐标，为飞机的北向坐标的微分量，为飞机的东向坐标的微分量，为飞机的高度坐标的微分量。

(4)通过数值积分进行单步预测结果的计算

将读入的初始飞行状态数据和单步预测的微分方程中的微分量进行数值积分，从而计算得出单步预测的飞行数据结果，数值积分的步长为t_s，单步决策步长为ΔT(时间长度)。

步骤四：计算支付函数矩阵

(1)建立支付函数模型

在本发明中，空战决策的进行需要建立一个可以有效反应红蓝双方态势的博弈支付函数，即对空战双方的态势进行定量的数值描述。由于到双机之间的角度关系和距离关系是当前态势的影响因素。

角度因素的评分函数为

$S_a=1-\frac{w_R+w_B}{180}---\left(7\right)$

其中，S_a为角度评分函数，w_R为红方无人作战飞机与两架无人作战飞机连线的夹角，w_B为蓝方无人作战飞机与两架无人作战飞机连线的夹角。

R＝(x_b-x_r)²+(y_b-y_r)²+(h_b-h_r)² (8)

$w_R=\frac{(x_b-x_r){\mathrm{cos\μ}}_r{\mathrm{cos\χ}}_r+(y_b-y_r){\mathrm{cos\μ}}_r{\mathrm{sin\χ}}_r+(h_b-h_r){\mathrm{sin\μ}}_r}{R}---\left(9\right)$

$w_B=\frac{(x_b-x_r){\mathrm{cos\μ}}_b{\mathrm{cos\χ}}_b+(y_b-y_r){\mathrm{cos\μ}}_b{\mathrm{sin\χ}}_b+(h_b-h_r){\mathrm{sin\μ}}_b}{R}---\left(10\right)$

其中，R为红蓝双方飞机之间的相对距离，x_r，x_b分为红方和蓝方的北向坐标，y_r，y_b分为红方和蓝方的东向坐标，h_r，h_b分为红方和蓝方的高度坐标，μ_r，μ_b分别为红方和蓝方的航迹倾斜角，χ_r，χ_b分别为红方和蓝方的航迹方位角。

距离因素的评分函数为

$S_r=e^{-(R-\stackrel{\‾}{r})/K}---\left(11\right)$

因此可以建立博弈的支付函数为

$S=S_a{S}_r=(1-\frac{w_R+w_B}{180})\left(e^{-(R-\stackrel{\‾}{r})/K}\right)---\left(12\right)$

其中，S为博弈支付函数，S_r为距离评分函数，为平均机炮攻击范围，K为灵敏度系数。

本发明中，所采用的支付函数计算结果值越大对红方越为有利，对蓝方越为不利。反之值越小对红方越为不利，对蓝方越为有利。

(2)计算支付函数并形成评分矩阵

对步骤三中按照博弈棋局计算得到的单步预测结果进行量化评分，将单步预测得到的结果信息代入支付函数的公式中进行计算，按照博弈支付函数的评分方式进行评分。将计算得到的评分结果按照原先的指令组合位置，以矩阵的形式储存为

其中，Score表示评分矩阵，S(m,n)表示对博弈矩阵中m行n列的指令单步预测结果的支付函数值。

步骤五：通过极大极小算法进行策略选择

极大极小算法是当代机器博弈搜索的基础，其始终是站在博弈一方的立场上对博弈系统的棋局进行估值，在进行估值时，将有利于这一方的棋局给予一个较高价值的数值，不利于这一方的给予一个较低价值的数值，博弈双方优劣势不明显的棋局给予一个中间价值的数值。

在本发明中，以红方为我方，根据步骤四中的评分矩阵计算结果，首先从每一列中选取支付函数值最小的一项，即对红方最为不利，对蓝方最为有利的指令组合情况，构成由极小值组成的行向量Score_min为

Score_min＝[S_min(1) S_min(2) … S_min(n)]_1×n (14)

其中，Score_min为由每列最小值组成的行向量，S_min(n)表示第n列的最小值。

然后再从每列的最小值中选取其中的最大值，即这些不利情况中对红方最为有利，对蓝方最为不利的指令组合情况，即在由极小值构成的行向量中找出其中的最大值Score_min,max为

Score_min,max＝max(S_min(1) S_min(2) … S_min(n)) (15)

记下最终选取的机动组合在矩阵中所在的位置，该位置所对应的机动动作即为最终决策产生的机动指令。

步骤六：更新六自由度飞机运动学和动力学方程

(1)计算飞机的输入量指令信号

将步骤五中决策产生的机动指令信号给到飞机的自动驾驶仪结构中，根据步骤一的控制律结构，计算出六自由度飞机模型所需的输入量

U＝[δ_T δ_e δ_a δ_r]

其中，U为控制量的集合，δ_T为油门杆推力输入，δ_e为升降舵偏角指令，δ_a为副翼偏角指令，δ_r为方向舵偏角指令。

(2)更新六自由度无人机的飞行状态

空气动力和力矩的计算公式为

$\left\{\begin{array}{c}L=C_L{\mathrm{QS}}_w=\frac{1}{2}{C}_L{\mathrm{\ρV}}^2{S}_w\\ D=C_D{\mathrm{QS}}_w=\frac{1}{2}{C}_D{\mathrm{\ρV}}^2{S}_w\\ Y=C_Y{\mathrm{QS}}_W=\frac{1}{2}{C}_Y{\mathrm{\ρV}}^2{S}_w\\ \stackrel{\‾}{L}=C_l{\mathrm{QS}}_{w}b=\frac{1}{2}{C}_l{\mathrm{\ρV}}^2{S}_{w}b\\ M=C_m{\mathrm{QS}}_w{c}_A=\frac{1}{2}{C}_m{\mathrm{\ρV}}^2{S}_w{c}_A\\ N=C_n{\mathrm{QS}}_{w}b=\frac{1}{2}{C}_n{\mathrm{\ρV}}^2{S}_{w}b\end{array}\right.---\left(17\right)$

其中，L为升力，D为阻力，Y为侧力，为绕机体轴x轴的转动力矩，M为绕机体轴y轴的转动力矩，N为绕机体轴z轴的转动力矩。C_L为无量纲的升力系数，C_D为无量纲的阻力系数，C_Y为无量纲的侧力系数，C_l为无量纲的滚转力矩系数，C_m为无量纲的俯仰力矩系数，C_n为无量纲的偏航力矩系数。Q为动压，S_w为机翼参考面积，ρ为空气密度，V为飞行速度，b为机翼展长，c_A机翼平均几何弦长。

六自由度飞机的导航方程组为

$\left\{\begin{array}{c}{\stackrel{\·}{x}}_g=u\cos \mathrm{\θ}\cos \mathrm{\ψ}+v(\sin \mathrm{\φ}\sin \mathrm{\θ}\cos \mathrm{\ψ}-\cos \mathrm{\φ}\sin \mathrm{\ψ})+w(\sin \mathrm{\φ}\sin \mathrm{\ψ}+\cos \mathrm{\φ}\sin \mathrm{\θ}\cos \mathrm{\ψ})\\ {\stackrel{\·}{y}}_g=u\cos \mathrm{\θ}\sin \mathrm{\ψ}+v(\sin \mathrm{\φ}\sin \mathrm{\θ}\sin \mathrm{\ψ}+\cos \mathrm{\φ}\sin \mathrm{\ψ})+w(-\sin \mathrm{\φ}\cos \mathrm{\ψ}+\cos \mathrm{\φ}\sin \mathrm{\θ}\sin \mathrm{\ψ})\\ {\stackrel{\·}{h}}_g=u\sin \mathrm{\θ}-v\sin \mathrm{\φ}\cos \mathrm{\θ}-w\cos \mathrm{\φ}\cos \mathrm{\θ}\end{array}\right.---\left(18\right)$

六自由度飞机的运动方程组为

$\left\{\begin{array}{c}\stackrel{\·}{\mathrm{\φ}}=p+(rcos\mathrm{\φ}+qsin\mathrm{\φ})tan\mathrm{\θ}\\ \stackrel{\·}{\mathrm{\θ}}=q\cos \mathrm{\φ}-rsin\mathrm{\φ}\\ \stackrel{\·}{\mathrm{\ψ}}=\frac{1}{\cos \mathrm{\θ}}(rcos\mathrm{\φ}+qsin\mathrm{\φ})\end{array}\right.---\left(19\right)$

六自由度飞机的力方程组为

$\left\{\begin{array}{c}\stackrel{\·}{V}=\frac{u\stackrel{\·}{u}+v\stackrel{\·}{v}+w\stackrel{\·}{w}}{V}\\ \stackrel{\·}{\mathrm{\α}}=\frac{u\stackrel{\·}{w}-w\stackrel{\·}{u}}{u^2+w^2}\\ \stackrel{\·}{\mathrm{\β}}=\frac{\stackrel{\·}{v}V-v\stackrel{\·}{V}}{V^2\cos \mathrm{\β}}\\ F_x=T+Lsin\mathrm{\α}-Ycos\mathrm{\α}sin\mathrm{\β}-Dcos\mathrm{\α}cos\mathrm{\β}\\ F_y=Ycos\mathrm{\β}-Dsin\mathrm{\β}\\ F_x=-L\cos \mathrm{\α}-Ysin\mathrm{\α}sin\mathrm{\β}-Dsin\mathrm{\α}cos\mathrm{\β}\end{array}\right.---\left(20\right)$

六自由度飞机的力矩方程组为

$\{\begin{array}{c}\stackrel{\·}{p}=\frac{1}{I_x{I}_z-{I_{xz}}^2}\[I_z\stackrel{\‾}{L}+I_{xz}N+(I_x-I_y+I_z)I_{xy}pq+(I_y{I}_z-{I_z}^2+{I^2}_{xz})qr\]\\ \stackrel{\·}{q}=\frac{1}{I_y}\[M-I_{xz}(p^2-r^2)\]\\ \stackrel{\·}{r}=\frac{1}{I_x{I}_y-{I_{xz}}^2}\end{array}---\left(21\right)$

其中，x_g为飞机的北向坐标，y_g为飞机的东向坐标，h_g为飞机的高度坐标，为飞机的北向坐标的微分量，为飞机的东向坐标的微分量，为飞机的高度坐标的微分量。u为机体轴系x轴方向的速度，v为机体轴系y轴方向的速度，w为机体轴z轴方向的速度；V为飞行速度，为飞行速度的微分量，α为飞机的迎角，为迎角的微分量，β为侧滑角，为侧滑角的微分量。F_x为机体x轴方向的合外力，F_y为机体y轴方向的合外力，F_z为机体z轴方向的合外力；φ为滚转角，θ为俯仰角，ψ为偏航角，为滚转角的微分量，为俯仰角的微分量，为偏航角的微分量；p为滚转角速度，q为俯仰角速度，r为偏航角速度，为滚转角速度的微分量，为俯仰角速度的微分量，为偏航角速度的微分量。I_x为飞机绕体轴x轴的转动惯量，I_y为飞机绕体轴y轴的转动惯量，I_z为飞机绕体轴z轴的转动惯量，I_xz为转动惯量积。

步骤七：判断是否达到空战终止条件

(1)判断是否达到开火条件

本发明假设机载航炮有效射程为r_m，攻击范围为以飞机的机体轴x轴为轴线，Γ为半顶角的类锥形空间区域为机载航炮的有效攻击范围如图5所示。当目标方位角(攻击目标与无人作战飞机的连线与x轴间的夹角)小于Γ，同时目标与攻击方的相对距离小于机炮射程时，满足开火条件，仿真结束。

(2)判断是否达到仿真时长

本发明可以设定空战仿真时间长度为T_max，如果当前仿真时间t大于仿真时间长度T_max是，则空战仿真过程结束；否则返回步骤三。整个空战仿真的流程图如附图1所示。

3、优点及效果：

本发明提出了一种基于单步预测矩阵博弈的近距空战自主决策方法。该方法主要针对六自由度非线性无人作战飞机模型进行空战决策，相比于三自由度质点模型更具有实际应用价值；同时，采用单步预测的矩阵博弈方法，将现有的基于机动动作库的矩阵博弈方法更改为基于指令模型的机动库，只需要单步步长的预测，从而有效减小了决策时间，满足了空中验证的实时性需求，可以更好地适应复杂动态的战场环境变化，进而有效提高了无人作战飞机在近距格斗中的作战能力。

【附图说明】

图1空战仿真流程图。

图2航迹倾斜角自动驾驶仪结构。

图3自动油门系统结构。

图4滚转角自动驾驶仪结构。

图5航炮攻击范围示意图。

图6航迹倾斜角响应曲线。

图7速度保持响应曲线。

图8滚转角响应曲线。

图9空战仿真飞行曲线。

图中标号及符号说明如下：

ΔT——决策周期步长

N——不满足条件(否)

Y——满足条件(是)

T_max——最大仿真时长

μ——航迹倾斜角

φ——滚转角

μ_com——航迹倾斜角指令信号

Δn_z——滚转角补偿的过载指令

q——俯仰角速度

n_z——法向过载

V_com——飞行速度指令

V——飞机速度

φ_com——滚转角指令信号

μ——航迹倾斜角

p——滚转角速度

r——偏航角速度

n_y——侧向过载

α——迎角

cos——余弦函数

sin——正弦函数

r_m——飞机的机炮射程

Γ——机炮攻击范围半顶角

deg——角度单位

t——时间

s——时间单位

x——北向坐标

y——东向坐标

h——高度坐标

m——距离单位

【具体实施方式】

下面通过一个具体的非线性六自由度无人机模型的近距空战格斗场景算例来验证本发明所提出的空战决策方法的有效性。由于F-16战斗机模型的各项数据比较公开，因此，本算例选用F-16的飞机模型作为空战飞机模型。实验计算机配置为Pentium处理器，2.50Ghz主频，1G内存，软件为MATLAB 2004版本。

步骤一：搭建六自由度非线性无人作战飞机控制律结构

(1)无人作战飞机配平

本发明中，飞机模型的采样周期为T＝0.01s，单位是秒，选用的配平高度h_trim＝6000m，配平速度为V_trim＝180m/s，配平迎角为α_trim＝3.4034deg，经过配平之后，飞机的俯仰角为θ＝3.4034deg，飞机的配平推力大小δ_Ttrim＝2085.1lbs，飞机的升降舵配平偏转角δ_etrim＝-2.2057deg，副翼配平偏角δ_atrim＝0deg，方向舵配平偏角δ_rtrim＝0deg，其余飞机的状态量的初始值都为零。其中，m为长度单位，deg为角度单位，lbs为推力单位。

(2)航迹倾斜角自动驾驶仪设计

设计如图2所示的航迹倾斜角自动驾驶仪，其中，俯仰角速度反馈的增益k_q＝0.5，航迹角增益为k_μ＝0.23，PI控制器为俯仰角速度陀螺传感器的传递函数为法向过载传感器传递函数为滚转角补偿器为

此时，给定单位航迹倾斜角指令信号，此时飞机航迹倾斜角响应曲线如图6所示。

(3)自动油门系统设计

设计如图3所示的自动油门系统，其中的PI控制器为k_V＝4.5+1/s，油门增益为

此时，给定单位航迹倾斜角指令信号的同时，给定V_com＝200m/s的飞行速度指令，此时飞机速度响应曲线如图7所示。

(4)滚转角自动驾驶仪设计

设计如图4所示的滚转角自动驾驶仪，其中，滚转角增益为k_φ＝2，交联增益为k_ARI＝-0.12，增益k₁＝0.5，k₂＝0.1，k₃＝0.45，超前滞后环节为高通滤波器为侧向过载滤波器迎角传感器传递函数滚转角速度陀螺传递函数为偏航角速度陀螺传递函数为

此时，给定单位滚转角的指令信号，此时飞机滚转角响应曲线如图8所示。

步骤二：初始化矩阵博弈棋局

(1)初始化空战双方初始飞行状态

在本算例中，红方飞机的初始位置信息为(x_g，y_g，h_g)＝(0,0,1500)，单位为米；初始姿态信息为(φ,θ,ψ,α,β)＝(0,3.4034,-90,3.4034,0)，单位为度；初始速度为V＝180m/s，并保持不变；初始欧拉角信息(p,q,r)＝(0,0,0)，单位为度每秒。蓝方飞机的初始位置信息为(x_g，y_g，h_g)＝(3000,3000,1500)，单位为米；初始姿态信息为(φ,θ,ψ,α,β)＝(0,3.4034,90,3.4034,0)，单位为度；初始速度为V＝165m/s；初始欧拉角信息(p,q,r)＝(0,0,0)，单位为度每秒。

(2)初始化离散指令机动库

为了体现不同飞机的机动特性，分别赋予红方和蓝方不同的候选指令机动库。红方的速度指令为V_com＝180m/s，蓝方的速度指令为V_com＝165m/s。

红方和蓝方的离散指令机动库分别为

${\mathrm{Lib}}_{RED}=\[{\mathrm{\μ}}_{com},{\mathrm{\φ}}_{com}\]=\left[\begin{array}{cc}0& 0\\ 2& 0\\ 5& 0\\ -2& 0\\ -5& 0\\ 0& -45\\ 0& 45\\ 2& 45\\ 2& -45\\ -2& 45\\ -2& -45\\ 5& 45\\ 5& -45\\ -5& 45\\ -5& -45\end{array}\right],{\mathrm{Lib}}_{BLUE}=\[{\mathrm{\μ}}_{com},{\mathrm{\φ}}_{com}\]=\left[\begin{array}{cc}0& 0\\ 2& 0\\ 5& 0\\ -2& 0\\ -5& 0\\ 0& -30\\ 0& 30\\ 2& 30\\ 2& -30\\ -2& 30\\ -2& -30\\ 5& 30\\ 5& -30\\ -5& 30\\ -5& -30\end{array}\right]$

其中，角度的单位为度。

(3)初始化候选机动指令的博弈棋局

将红方和蓝方的离散指令机动库中的指令进行组合，由于红方和蓝方分别有15中候选机动指令，因此形成一个m×n＝15×15的指令博弈棋局。

步骤三：根据博弈棋局进行单步预测计算

按照步骤三中的流程进行计算，其中T_μ＝2.5，单位为秒，重力加速度g＝9.8，单位为米每二次方秒。单步决策步长为ΔT＝1，单位为秒，数值积分的步长t_s＝0.1，单位为秒。

步骤四：计算支付函数矩阵

按照步骤四中的流程进行计算，其中机炮平均范围单位为米。灵敏度系数K＝1000。

步骤五：通过极大极小算法进行策略选择

按照步骤五的计算步骤进行计算，选择出红方和蓝方的待选机动指令。

步骤六：更新六自由度飞机运动学和动力学方程

F-16战斗机的固有参数设置如下：飞机机体轴的转动惯量分别为I_x＝12874.8kg·m²，I_y＝75673.6kg·m²，I_z＝85552.1kg·m²，惯量积I_xz＝1331.4kg·m²，I_xy＝0，I_zy＝0，其单位为千克乘以米的平方。

F-16战斗机的控制输入约束设置如下：δ_T的允许范围为1000～19000lbs,其变化速率的允许范围为-10000～10000lbs/s；δ_e的允许范围为-25～25deg，其变化速率的允许范围为-60～60deg/s；δ_a的允许范围为-21.5～21.5deg，其变化速率的允许范围为-80～80deg/s；δ_r的允许范围为-30～30deg，其变化速率的允许范围为-120～120deg/s。其中，lbs为推理的单位，deg/s为角速度单位，为度每秒。

步骤七：判断是否达到空战终止条件

机炮的有效射程r_m＝1000，为单位米，以飞机的机体轴x轴为轴线的类锥形空间区域半顶角Γ＝30，单位为度。

Claims (2)

1.一种基于单步预测矩阵博弈的近距空战自主决策方法，其特征在于：该方法的实现步骤如下：

步骤1：搭建六自由度非线性无人作战飞机控制律结构

步骤1.1无人作战飞机配平

设置飞机的油门杆配平推力δ_Ttrim，升降舵配平偏角δ_etrim，副翼配平偏角δ_atrim，方向舵配平偏角δ_rtrim在给定的配平高度h_trim，配平迎角α_trim和飞行配平速度V_trim的条件下，使得无人作战飞机所受合力和合力矩为零，从而保持平飞状态，即确定飞机在配平状态下的参数；

步骤1.2航迹倾斜角自动驾驶仪设计

在配平状态下，通过引入俯仰角速度q，飞机机体轴法向过载n_z，滚转角φ，航迹倾斜角μ的反馈信号，设计航迹倾斜角保持的自动驾驶仪，将控制指令以升降舵偏角指令δ_e的形式给到飞机模型，进而实现飞机的航迹倾斜角对航迹倾斜角指令信号μ_com的跟踪；

步骤1.3自动油门系统设计

在配平状态下，为了在飞机爬升时对其动力进行补偿，通过引入飞行速度V的反馈信号，将飞行速度指令V_com；给到飞机的油门杆推力输入δ_T，从而实现飞行速度对飞行速度指令V_com的跟踪保持；

步骤1.4滚转角自动驾驶仪设计

在配平状态下，通过引入滚转角速度p，偏航角速度r，滚转角φ，迎角α，和机体轴侧向过载n_y作为反馈信号，同时引入方向舵偏角指令δ_r，将滚转角的控制指令给到副翼偏角指令δ_a，设计滚转角保持的自动驾驶仪，从而实现对滚转角指令信号φ_com的跟踪；

步骤2：初始化矩阵博弈棋局

步骤2.1初始化空战双方初始飞行状态

设定博弈的双方分别为红方和蓝方，对空战双方的初始飞行状态进行初始化设置，包括红方和蓝方，飞机的位置信息和姿态信息，即

X＝[x_g y_g h_g φ θ ψ V α β p q r]

其中，X为状态量的集合，x_g为飞机的北向坐标，y_g为飞机的东向坐标，h_g为飞机的高度坐标，φ为飞机的滚转角，θ为飞机的俯仰角，ψ为飞机的偏航角，V为飞行速度，α为飞机的迎角，β为飞机的侧滑角，p为滚转角速度，q为俯仰角速度，r为偏航角速度；

步骤2.2初始化离散指令机动库

采取按照飞机机动指令离散信号的形式构建博弈所需的待选机动动作库的方法；认为空战过程中飞行速度V保持不变，同时飞行速度指令V_com保持不变，将纵向的航迹倾斜角指令信号μ_com和横侧向的滚转角指令信号φ_com进行组合，从而形成执行平飞、爬升、俯冲、协调转弯、爬升转弯和俯冲转弯的战术机动的机动指令信号；

离散指令机动库为

${{\mathrm{Lib}}_{RED}}^{\left(i\right)}=\[{\mathrm{\μ}}_{RED}{{}^{\left(i\right)}}_{com},{\mathrm{\φ}}_{RED}{{}^{\left(i\right)}}_{com}\]---\left(1\right)$

${{\mathrm{Lib}}_{BLUE}}^{\left(j\right)}=\[{\mathrm{\μ}}_{BLUE}{{}^{\left(j\right)}}_{com},{\mathrm{\φ}}_{BLUE}{{}^{\left(j\right)}}_{com}\]---\left(2\right)$

其中，Lib_RED为候选指令机动库集合，i表示第i个机动指令组合，表示第i组航迹倾斜角指令，表示第i组滚转角指令，Lib_BLUE为候选指令机动库集合，j表示第j个机动指令组合，表示第j组航迹倾斜角指令，表示第j组滚转角指令；

步骤2.3初始化候选机动指令的博弈棋局

博弈棋局表示在进行决策时，红蓝双方飞机的所有候选机动指令组合的集合，以矩阵的形式给出；认为红方一共有m种候选机动组合，蓝方一共有n种候选机动组合，则候选机动的矩阵为

步骤3：根据博弈棋局进行单步预测计算

步骤3.1读入当前时刻的飞行状态

读取红蓝双方飞机当前时刻的飞行状态，以进行单步预测的计算，读取的飞行状态信息为

Y＝[x_g0 y_g0 h_g0 V₀ μ₀ χ₀]

其中，Y为状态量的集合，x_g0为飞机的初始北向坐标，y_g0为飞机的初始东向坐标，h₀为飞机的初始高度坐标，V₀为飞机的初始气流速度，μ₀为飞机的初始航迹倾斜角，χ₀为飞机的初始航迹方位角；

步骤3.2读入矩阵棋局中的候选机动指令

读取步骤2中候选机动指令的博弈棋局，分别得到该矩阵中红蓝双方不同的机动指令组合，即

{Lib_BLUE,Lib_RED}＝{[μ_BLUEcom,φ_BLUEcom],[μ_REDcom,φ_REDcom]} (4)

其中，Lib_RED为候选指令机动库集合，Lib_BLUE为候选指令机动库集合，μ_REDcom表示红方的航迹倾斜角指令，φ_REDcom表示红方的滚转角指令，μ_BLUEcom表示蓝方航迹倾斜角指令，φ_BLUEcom表示蓝方滚转角指令；

步骤3.3进行单步预测微分方程计算

对无人机的飞行状态建立预测模型进行预测，建立飞机的预测航迹角变化的微分方程为

$\left\{\begin{array}{c}\stackrel{\·}{\mathrm{\μ}}=\frac{1}{T_{\mathrm{\μ}}}({\mathrm{\μ}}_{com}-\mathrm{\μ})\\ n_f=\frac{1}{{\mathrm{cos\φ}}_{com}}(\frac{V\stackrel{\·}{\mathrm{\μ}}}{g}+cos\mathrm{\μ})\\ \stackrel{\·}{\mathrm{\χ}}=\frac{g}{v\cos \mathrm{\μ}}{n}_f{\mathrm{sin\φ}}_{com}\end{array}\right.---\left(5\right)$

建立飞机的预测位置信息的微分方程为

$\left\{\begin{array}{c}{\stackrel{\·}{x}}_g=Vcos\mathrm{\μ}cos\mathrm{\χ}\\ {\stackrel{\·}{y}}_g=Vcos\mathrm{\μ}sin\mathrm{\χ}\\ {\stackrel{\·}{h}}_g=Vsin\mathrm{\μ}\end{array}\right.---\left(6\right)$

其中，μ为航迹倾斜角，为航迹倾斜角的微分，T_μ为航迹角响应的一阶等效模型的时间常数，μ_com为航迹倾斜角指令信号，n_f为飞机航迹坐标下的法向过载，φ_com为滚转角指令信号，g为重力加速度，为航迹方位角的微分，χ为航迹方位角，x_g为飞机的北向坐标，y_g为飞机的东向坐标，h_g为飞机的高度坐标，为飞机的北向坐标的微分量，为飞机的东向坐标的微分量，为飞机的高度坐标的微分量；

步骤3.4通过数值积分进行单步预测结果的计算

将读入的初始飞行状态数据和单步预测的微分方程中的微分量进行数值积分，计算得出单步预测的飞行数据结果，数值积分的步长为t_s，单步决策步长为ΔT；

步骤4：计算支付函数矩阵

步骤4.1建立支付函数模型

对空战双方的态势进行定量的数值描述；由于到双机之间的角度关系和距离关系是当前态势的影响因素；故

角度因素的评分函数为

$S_a=1-\frac{w_R+w_B}{180}---\left(7\right)$

其中，S_a为角度评分函数，w_R为红方无人作战飞机与两架无人作战飞机连线的夹角，w_B为蓝方无人作战飞机与两架无人作战飞机连线的夹角；

R＝(x_b-x_r)²+(y_b-y_r)²+(h_b-h_r)² (8)

$w_R=\frac{(x_b-x_r){\mathrm{cos\μ}}_r{\mathrm{cos\χ}}_r+(y_b-y_r){\mathrm{cos\μ}}_r{\mathrm{sin\χ}}_r+(h_b-h_r){\mathrm{sin\μ}}_r}{R}---\left(9\right)$

$w_B=\frac{(x_b-x_r){\mathrm{cos\μ}}_b{\mathrm{cos\χ}}_b+(y_b-y_r){\mathrm{cos\μ}}_b{\mathrm{sin\χ}}_b+(h_b-h_r){\mathrm{sin\μ}}_b}{R}---\left(10\right)$

其中，R为红蓝双方飞机之间的相对距离，x_r，x_b分为红方和蓝方的北向坐标，y_r，y_b分为红方和蓝方的东向坐标，h_r，h_b分为红方和蓝方的高度坐标，μ_r，μ_b分别为红方和蓝方的航迹倾斜角，χ_r，χ_b分别为红方和蓝方的航迹方位角；

距离因素的评分函数为

$S_r=e^{-(R-\stackrel{\‾}{r})/K}---\left(11\right)$

因此建立博弈的支付函数为

$S=S_a{S}_r=(1-\frac{w_R+w_B}{180})\left(e^{-(R-\stackrel{\‾}{r})/K}\right)---\left(12\right)$

其中，S为博弈支付函数，S_r为距离评分函数，为平均机炮攻击范围，K为灵敏度系数；

步骤4.2计算支付函数并形成评分矩阵

对步骤3中按照博弈棋局计算得到的单步预测结果进行量化评分，将单步预测得到的结果信息代入支付函数的公式中进行计算，按照博弈支付函数的评分方式进行评分；将计算得到的评分结果按照原先的指令组合位置，以矩阵的形式储存为

其中，Score表示评分矩阵，S(m,n)表示对博弈矩阵中m行n列的指令单步预测结果的支付函数值；

步骤5：通过极大极小算法进行策略选择

根据步骤4中的评分矩阵计算结果，首先从每一列中选取支付函数值最小的一项，即对红方最为不利，对蓝方最为有利的指令组合情况，构成由极小值组成的行向量Score_min为

Score_min＝[S_min(1) S_min(2) … S_min(n)]_1×n (14)

其中，Score_min为由每列最小值组成的行向量，S_min(n)表示第n列的最小值；

然后再从每列的最小值中选取其中的最大值，即这些不利情况中对红方最为有利，对蓝方最为不利的指令组合情况，即在由极小值构成的行向量中找出其中的最大值Score_min,max为

Score_min,max＝max(S_min(1) S_min(2) … S_min(n)) (15)

记下最终选取的机动组合在矩阵中所在的位置，该位置所对应的机动动作即为最终决策产生的机动指令；

步骤6：更新六自由度飞机运动学和动力学方程

步骤6.1计算飞机的输入量指令信号

将步骤5中决策产生的机动指令信号给到飞机的自动驾驶仪结构中，根据步骤1的控制律结构，计算出六自由度飞机模型所需的输入量

U＝[δ_T δ_e δ_a δ_r]

其中，U为控制量的集合，δ_T为油门杆推力输入，δ_e为升降舵偏角指令，δ_a为副翼偏角指令，δ_r为方向舵偏角指令；

步骤6.2更新六自由度无人机的飞行状态

空气动力和力矩的计算公式为

$\left\{\begin{array}{c}L=C_L{\mathrm{QS}}_w=\frac{1}{2}{C}_L{\mathrm{\ρV}}^2{S}_w\\ D=C_D{\mathrm{QS}}_w=\frac{1}{2}{C}_D{\mathrm{\ρV}}^2{S}_w\\ Y=C_Y{\mathrm{QS}}_w=\frac{1}{2}{C}_Y{\mathrm{\ρV}}^2{S}_w\\ \stackrel{\‾}{L}=C_l{\mathrm{QS}}_{w}b=\frac{1}{2}{C}_l{\mathrm{\ρV}}^2{S}_{w}b\\ M=C_m{\mathrm{QS}}_w{c}_A=\frac{1}{2}{C}_m{\mathrm{\ρV}}^2{S}_w{c}_A\\ N=C_n{\mathrm{QS}}_{w}b=\frac{1}{2}{C}_n{\mathrm{\ρV}}^2{S}_{w}b\end{array}\right.---\left(17\right)$

其中，L为升力，D为阻力，Y为侧力，为绕机体轴x轴的转动力矩，M为绕机体轴y轴的转动力矩，N为绕机体轴z轴的转动力矩；C_L为无量纲的升力系数，C_D为无量纲的阻力系数，C_Y为无量纲的侧力系数，C_l为无量纲的滚转力矩系数，C_m为无量纲的俯仰力矩系数，C_n为无量纲的偏航力矩系数；Q为动压，S_w为机翼参考面积，ρ为空气密度，V为飞行速度，b为机翼展长，c_A机翼平均几何弦长；

六自由度飞机的导航方程组为

$\left\{\begin{array}{c}{\stackrel{\·}{x}}_g=u\cos \mathrm{\θ}\cos \mathrm{\ψ}+v(\sin \mathrm{\φ}\sin \mathrm{\θ}\cos \mathrm{\ψ}-\cos \mathrm{\φ}\sin \mathrm{\ψ})+w(\sin \mathrm{\φ}\sin \mathrm{\ψ}+\cos \mathrm{\φ}\sin \mathrm{\θ}\cos \mathrm{\ψ})\\ {\stackrel{\·}{y}}_g=u\cos \mathrm{\θ}\sin \mathrm{\ψ}+v(\sin \mathrm{\φ}\sin \mathrm{\θ}\sin \mathrm{\ψ}+\cos \mathrm{\φ}\sin \mathrm{\ψ})+w(-\sin \mathrm{\φ}\cos \mathrm{\ψ}+\cos \mathrm{\φ}\sin \mathrm{\θ}\sin \mathrm{\ψ})\\ {\stackrel{\·}{h}}_g=u\sin \mathrm{\θ}-v\sin \mathrm{\φ}\cos \mathrm{\θ}-w\cos \mathrm{\φ}\cos \mathrm{\θ}\end{array}\right.---\left(18\right)$

六自由度飞机的运动方程组为

$\left\{\begin{array}{c}\stackrel{\·}{\mathrm{\φ}}=p+(rcos\mathrm{\φ}+qsin\mathrm{\φ})tan\mathrm{\θ}\\ \stackrel{\·}{\mathrm{\θ}}=q\cos \mathrm{\φ}-rsin\mathrm{\φ}\\ \stackrel{\·}{\mathrm{\ψ}}=\frac{1}{\cos \mathrm{\θ}}(rcos\mathrm{\φ}+q\sin \mathrm{\φ})\end{array}\right.---\left(19\right)$

六自由度飞机的力方程组为

$\left\{\begin{array}{c}\stackrel{\·}{V}=\frac{u\stackrel{\·}{u}+v\stackrel{\·}{v}+w\stackrel{\·}{w}}{V}\\ \stackrel{\·}{\mathrm{\α}}=\frac{u\stackrel{\·}{w}-w\stackrel{\·}{u}}{u^2+w^2}\\ \stackrel{\·}{\mathrm{\β}}=\frac{\stackrel{\·}{v}V-v\stackrel{\·}{V}}{V^2\cos \mathrm{\β}}\\ F_x=T+Lsin\mathrm{\α}-Ycos\mathrm{\α}sin\mathrm{\β}-Dcos\mathrm{\α}cos\mathrm{\β}\\ F_y=Ycos\mathrm{\β}-Dsin\mathrm{\β}\\ F_x=-L\cos \mathrm{\α}-Ysin\mathrm{\α}sin\mathrm{\β}-Dsin\mathrm{\α}cos\mathrm{\β}\end{array}\right.---\left(20\right)$

六自由度飞机的力矩方程组为

$\left\{\begin{array}{c}\stackrel{\·}{p}=\frac{1}{I_x{I}_z-{I_{xz}}^2}\[I_z\stackrel{\‾}{L}+I_{xz}N+(I_x-I_y+I_z)I_{xz}pq+(I_y{I}_z-{I_z}^2+{I^2}_{xz})qr\]\\ \stackrel{\·}{q}=\frac{1}{I_y}\[M-I_{xz}(p^2-r^2)\]\\ \stackrel{\·}{r}=\frac{1}{I_x{I}_y-{I_{xz}}^2}\end{array}\right.---\left(21\right)$

其中，x_g为飞机的北向坐标，y_g为飞机的东向坐标，h_g为飞机的高度坐标，为飞机的北向坐标的微分量，为飞机的东向坐标的微分量，为飞机的高度坐标的微分量；u为机体轴系x轴方向的速度，v为机体轴系y轴方向的速度，w为机体轴z轴方向的速度；V为飞行速度，为飞行速度的微分量，α为飞机的迎角，为迎角的微分量，β为侧滑角，为侧滑角的微分量；F_x为机体x轴方向的合外力，F_y为机体y轴方向的合外力，F_z为机体z轴方向的合外力；φ为滚转角，θ为俯仰角，ψ为偏航角，为滚转角的微分量，为俯仰角的微分量，为偏航角的微分量；p为滚转角速度，q为俯仰角速度，r为偏航角速度，为滚转角速度的微分量，为俯仰角速度的微分量，为偏航角速度的微分量；I_x为飞机绕体轴x轴的转动惯量，I_y为飞机绕体轴y轴的转动惯量，I_z为飞机绕体轴z轴的转动惯量，I_xz为转动惯量积；

步骤7：判断是否达到空战终止条件

7.1判断是否达到开火条件

假设机载航炮有效射程为r_m，攻击范围为以飞机的机体轴x轴为轴线，Γ为半顶角的类锥形空间区域为机载航炮的有效攻击范围如图5所示；当目标方位角小于Γ，同时目标与攻击方的相对距离小于机炮射程时，满足开火条件，仿真结束；

7.2判断是否达到仿真时长

设定空战仿真时间长度为T_max，如果当前仿真时间t大于仿真时间长度T_max是，则空战仿真过程结束；否则返回步骤3。

2.根据权利要求1所述的一种基于单步预测矩阵博弈的近距空战自主决策方法，其特征在于：在步骤4中，所采用的支付函数计算结果值越大对红方越为有利，对蓝方越为不利；反之值越小对红方越为不利，对蓝方越为有利。