CN107608335A

CN107608335A - 一种无人机飞行控制系统故障检测与故障分离的数据驱动方法

Info

Publication number: CN107608335A
Application number: CN201710825789.8A
Authority: CN
Inventors: 钟麦英; 李明虎; 赵岩; 周东华; 陈杰
Original assignee: Shandong University of Science and Technology
Current assignee: Shandong University of Science and Technology
Priority date: 2017-09-14
Filing date: 2017-09-14
Publication date: 2018-01-19
Anticipated expiration: 2037-09-14
Also published as: CN107608335B

Abstract

本发明涉及一种无人机飞行控制系统故障检测与故障分离的数据驱动方法。本发明以无人机正常和故障情况下的实际飞行数据为基础，利用正常数据建立基于多元统计分析的数据模型，获取T²和SPE两种监测统计量的阈值；然后利用故障数据计算故障发生后各采样点T²和SPE两种统计量，并与阈值比较从而获得故障检测结果，进而采用贡献图的方法，计算故障发生后各变量对两种统计结果的贡献值，进而得到累积贡献图。最终利用贡献图和累积贡献图，确定故障变量，获得故障分离结果。该故障检测与分离方法能够有效解决实际飞行系统中不同类型故障混叠问题，为无人机系统设计及故障定位与故障传播方法分析提供一种参考依据。

Description

一种无人机飞行控制系统故障检测与故障分离的数据驱动方法

技术领域

本发明涉及一种无人机飞行控制系统故障检测与故障分离的数据驱动方法，用于从数据驱动角度对无人机飞行控制系统故障进行定性及定量分析，属于无人机系统故障诊断技术领域。

背景技术

目前针对无人机飞行控制系统的故障诊断研究成果，多数是采用基于解析模型的方法进行研究。基于解析模型的方法，最主要特点是离不开精确的数学模型，而模型的不确定性使得故障诊断的难度增大。并且，目前多数研究使用的无人机模型为线性的，真实无人机的飞行过程是非线性的，而且受到大气扰动、外界干扰等一系列因素的影响，这使得故障诊断的难度进一步增大。而数据驱动的故障诊断方法不需要获取系统的精确模型，对于模型不确定的小型无人机进行故障诊断具有更大的优势。

目前故障诊断的数据驱动方法主要包括基于多元统计分析的方法、基于机器学习的方法以及信号处理的方法。其中，基于机器学习的故障诊断方法具有强自适应性，鲁棒性和抗干扰性等优点，但是存在模型结构的物理意义不明，算法复杂、收敛慢，在线实现困难，对复杂系统学习速度慢等不可避免的缺点。基于信号处理的故障诊断方法具有灵敏度高，克服噪声能力强等优点，但是存在计算量大，难以在线辨识，在线进行故障诊断等缺点。以上两种方法均难以用在模型不确定的小型无人机的故障检测与故障分离上，而动态主元分析作为基于多元统计分析的故障诊断方法中的一种，不需要对系统的结构和原理有深入的了解，完全基于系统运行过程中传感器的测量数据，而且适用于动态系统，算法简单，易于实现。因此，本文选择主元分析方法作为无人机飞行控制系统的故障检测与分离方法进行分析研究。

基于主元分析的故障诊断方法在工业过程中的应用比较成熟，在无人机上的应用才刚刚起步。白志强所撰论文“飞行控制系统故障检测研究与仿真软件开发[D].西北工业大学,2006.”将主元分析应用到无人机执行机构的故障诊断中，不仅能够检测执行机构卡死故障，而且能够检测舵面损伤故障即部分失效故障。邱宗江等所撰论文“无人机PCA故障检测与诊断技术研究[J].计算机工程与应用,2013,49(4):262-266.”将主元分析应用到无人机传感器的故障诊断中，能够实现传感器故障的故障检测和分离。可以发现目前的研究仅限于对执行机构或者传感器进行故障诊断，没有对执行机构以及传感器进行全面的故障诊断并分析其故障检测与分离，且大多针对的是无人机线性模型。尚未有文献或专利提出针对非线性无人机飞行控制系统进行全面的故障检测与故障分离方法的分析，故该课题有待深入研究。

发明内容

本发明要解决的技术问题是：通过从故障发生后获得的贡献图以及累积贡献图，对所发生的故障进行准确分离。利用多元统计分析中动态主元分析方法建立正常情况的数据模型，通过T²和SPE两种监测统计量实现在线故障检测；计算出故障数据中各变量对T²和SPE统计量的贡献值，进而得到贡献图以及累积贡献图，最终完成对系统的故障分离。

本发明解决上述技术问题采用的技术方案为：一种无人机飞行控制系统故障检测与故障分离的数据驱动方法，包括下列步骤：

步骤1：采集无人机正常平稳飞行以及各种故障情况下的输入输出数据，对需要采集的无人机输入输出变量介绍如下：

无人机执行机构包含升降舵、方向舵、副翼和油门杆这四个执行机构，其舵偏角分别用δ_e、δ_r、δ_a和δ_p表示。因此无人机输入数据包含以上4个变量。

无人机传感器主要包括攻角传感器、侧滑角传感器、空速表、高度表、垂直陀螺以及角速率陀螺，输出数据依次为攻角α、侧滑角β、空速V、高度H，俯仰角θ、偏航角ψ、滚转角φ、滚转角速度p、俯仰角速度q、偏航角速度r。因此无人机输出数据包含以上10个变量。

步骤2：将步骤1采集到的数据整理为矩阵形式，假设有N次测量，将每次采样点获得的无人机输入输出数据整理为向量x，假设x中包含m个变量，那么形成的数据矩阵为X：X＝[x₁，x₂，…,x_N]∈R^m×N，并利用移动均值法对数据进行平滑预处理；

步骤3：由步骤2得到平滑预处理的原始量测数据后，为了获得更多具有相关性的数据信息，需要掌握当前时刻数据与过去s个时刻数据的关系，s为采样数据间最大持续影响的采样次数，即动态阶数。在这里，我们通过平行分析法来确定动态阶数s。

步骤4：根据步骤3确定的动态阶数s，把前面s个时刻的量测数据的每个变量当成新的增广阵的变量，建立数据增广矩阵X_a。对该数据增广矩阵X_a进行动态主元分析，建立动态主元分析模型，为后续的故障检测与分离提供依据。动态主元分析方法的基本原理如下：

将原始的系统测量数据整理为数据矩阵，比如行数据代表所有的测量变量，列数据代表所有的测量时刻；然后把前面S个时刻的量测数据的每个变量当成新的增广阵的变量，构建新的增广数据矩阵；最后通过对数据矩阵进行线性变换(类似于线性投影)，将原始数据空间分解为两个正交的子空间，分别是主元子空间和残差子空间。主元子空间中有少部分不相关的重要变量可以代表系统的大部分信息，残差子空间的变量代表着系统的噪声误差等信息，这样就获得了一个基于正常数据的动态主元模型了。这个动态主元模型反映了系统正常运行时数据的结构关系，当系统处于正常状态时，当前时刻获取的测量数据经过动态主元分析分解后得到的得分向量将会投影在主元子空间里，如果当前时刻获得测量数据与建立的主元模型不匹配时，那么就可以推断系统出现了异常或者故障。

步骤5：根据步骤4建立的动态主元模型，通过计算新数据在主元子空间变化的程度以及在残差子空间变化的程度来检测系统是处于正常状态还是出现故障。为了提高故障检测与分离的准确率，采用将Hotelling的T²统计量和平方预测误差SPE统计量两种故障检测方法配合使用的方法，利用飞机正常飞行数据，计算两种监测统计量的阈值。其中：

T²统计量表示每个采样在变化趋势和幅值上偏离主元模型的程度，用来检验系统新的量测数据是否依然在主元确定的控制域内，是模型内部变化的一种测度。

SPE统计量表示采样向量在残差子空间上投影变化的大小，用来检验原始数据中没有被主元模型所解释的变化，是模型外部变化的一种测度。

步骤6：根据步骤4建立的动态主元模型，计算故障发生后各采样点的T²和SPE两种统计量，分别与步骤5得到的阈值比较。当T²统计量或者SPE统计量至少有一个会超出阈值时，即认为发生故障，从而实现对故障的检测；

步骤7：在步骤6完成故障检测的基础上，为了对故障进行分离，采用贡献图的方法进行分析。具体方法是：在故障发生后的每一个时刻，计算出每个变量对T²和SPE统计量的贡献值，并通过除以当前时刻的相应统计量进行归一化，从而获得在各个采样点各个变量对统计量的贡献率，最终将其以变量为纵坐标，时间为横坐标作图，即可得到贡献图；

步骤8：根据步骤7得到的贡献率，将故障发生后一段时间内各变量对统计量的贡献值进行累加，获得各变量对统计量的累积贡献率，以变量为横坐标作柱状图，即为累积贡献图；

步骤9：根据贡献图和累积贡献图理论：在一般情况下，故障变量对统计量的贡献率大于其它变量的贡献率，而故障变量的累积贡献率更是远大于其它变量的累计贡献率。则可在步骤7和步骤8的基础上，通过贡献图和累积贡献图确定故障变量，实现故障的分离；

进一步地，步骤2中所述平滑预处理的移动均值法为：

x_i＝(x_i+x_i-1+x_i-2+…+x_i-9)/10

为了较大限度地平滑由随机性所带来的偏差，这里取当前时刻的测量数据x_i以及过去9个时刻的测量数据x_i-1、x_i-2、…、x_i-9共10个测量数据进行平滑处理。

进一步地，步骤3所述的进行平行分析的具体方法为：

首先，取动态阶数s＝0，对正常飞行的数据矩阵X进行归一化，使其均值为零并具有单位方差，并计算归一化后数据矩阵的协方差阵：

对S进行奇异值分解，以获得特征值组成的对角矩阵Λ＝diag(λ₁,…λ_m)以及相应的特征向量组成的矩阵P，其中λ₁≥…≥λ_m≥0。将P中前k个特征向量组成的矩阵定义为载荷矩阵P_pc∈R^m×k，剩余的特征向量组成的矩阵定义为载荷矩阵P_res∈R^m×(m-k)。对于一个采样向量x∈R^m可以向载荷向量进行投影，获得x在每个特征向量方向上的得分t＝P_pc ^Tx。

然后，取s阶时间滞后情况下总的线性关系的个数r(s)＝0，判断第j＝m×(s+1)个主元是否代表着一个线性关系。是则令j＝j-1，r(s)＝r(s)+1，重复该判断；否，则进行下一步。

最后，计算第s阶包含的线性关系的个数r_new(0)＝0，并判断r_new(s)<0是否成立，若成立，则当前s即为所求的动态阶数，平行分析结束；反之，取s＝s+1，并把前面s个时刻的量测数据的每个变量当成新的增广阵的变量，建立数据增广矩阵X_a，代替原始数据矩阵X重新进行整个平行分析过程。

进一步地，步骤4所述的建立数据的动态主元模型具体步骤为：

首先，构建数据增广矩阵X_a。

其中N₁≥k+N。

其次，对矩阵X_a进行归一化处理，得到均值为零，方差为1的标准化后的数据矩阵X_S。

再次，对X_S的协方差矩阵进行奇异值分解，以获得特征值组成的对角矩阵Λ以及相应的特征向量组成的矩阵P：

最后，利用累积方差贡献率的方法来确定主元个数并获得负荷矩阵P。当累积方差贡献率CPV达到90％时，计算出来的主元可以综合系统的有效信息，此时的k即为主元个数。

Λ_pc＝diag(λ₁,…,λ_k),Λ_res＝diag(λ_k+1,…,λ_m)；P＝[P_pc P_res],P_pc∈R^m×k,P_res∈R^m ^×(m-k)。

进一步地，步骤5所述的计算T²和SPE两种监测统计量阈值的具体求解方法为：

置信度为α时T²统计量阈值求解公式如下：

其中，F_α代表置信度为α时的F分布。

置信度为α时SPE统计量阈值求解公式如下：

其中，C_α为置信度为α时的标准正态分布的分位数，

进一步地，步骤6所述的故障数据各采样点T²和SPE两种统计量具体求解方法为：

T²＝t_newΛ^-1t_new ^T＝x_new ^Tp_pcΛ_pc ^-1p_pc ^Tx_new

SPE＝e_newe_new ^T＝x_new(I-P_pcP_pc ^T)x_new ^T

其中，x_new为当前采样点的测量数据向量。

此时，便可以通过以下逻辑规则进行故障检测：

故障

正常。

进一步地，步骤7所述的获取贡献图的具体步骤为：

首先，在采样时刻j，计算变量i对T²和SPE的贡献值：

T² _ij＝x_j ^TM^(1/2)ξ_iξ_i ^TM^(1/2)x_j,M＝P_pcΛ_pc ^-1P_pc ^T

其中，x_j是采样时刻j的测量数据，ξ_i是单位矩阵的第i列，e_ij是采样时刻j变量i的残差，x_ij是采样时刻j变量i的测量值，是采样时刻j变量i的估计值。

然后，计算变量i对T²和SPE的累积贡献值：

由于测量数据矩阵的增广矩阵X包含时间滞后s，因此m×(s+1)个贡献值属于原始的m个变量，将属于同一个变量的s+1个贡献值累加起来作为这个原始变量的贡献值：

SPE^ori _i,j＝SPE_i,j+SPE_(i+m),j+…+SPE_(i+s*m),j,i＝1,…m

最后，采用归一化的方法，获得每个变量对统计量的贡献率。

无人机飞行控制系统是动态的，T²和SPE统计量的大小在每次采样都会发生变化。在每一个时刻各个变量的贡献值可以通过除以当前时刻的相应统计量进行归一化，从而获得每个变量对统计量的贡献率。

在j时刻变量i对T²统计量的贡献率为：

在j时刻变量i对SPE统计量的贡献率为：

CR_SPE_i,j＝SPE_ori _i,j/SPE_j

进一步地，步骤8所述的获取累积贡献图的具体步骤为：

为了更好的分离故障变量，我们需要得到故障发生后一段时间内各变量对统计量的累积贡献率，其可以通过将故障发生时间段内各个变量的贡献率分别进行累加得到：

其中，t_f是故障发生的时刻，n_f是故障发生后总共的采样数,是变量i对T²统计量的累积贡献率，SPE^CCR _i是变量i对SPE统计量的累积贡献率；

然后以变量为横坐标作累积贡献率的柱状图，即为累积贡献图。

进一步地，步骤9所述的获取故障变量的具体步骤为：

首先，根据平均值原理计算累积贡献图中的故障控制线，具体值为：

然后，在T²统计量累积贡献图以及SPE统计量累积贡献图中，将超过控制线的变量作为可能的故障变量。下面将对故障分离结果进行分类讨论，给出如下故障分离判别准则；

(1)当可能的故障变量为单个变量时，即可判定这个单一变量为根本的故障变量。

(2)当可能的故障变量为多个变量时，结合贡献图中故障发生初始时刻多个可能故障变量对统计量贡献率的变化情况来推断根本的故障变量。

(3)当T²统计量累积贡献图以及SPE统计量累积贡献图中故障变量一致时，即可判定这个变量为根本的故障变量。

(4)当T²统计量累积贡献图以及SPE统计量累积贡献图中故障变量不一致时，由于SPE统计量很难反应系统的动态特性，基于SPE统计量贡献图的故障分离结果大部分都不准确，此时以T²统计量贡献图的故障分离结果为准，以SPE统计量累积贡献图的故障分离结果为参考。

本发明与现有技术相比的优点在于：

(1)本发明首次将动态核主元应用于无人机飞行控制系统的故障诊断中，不仅无需已知无人机的解析模型、不受无人机机型限制，而且能获得很好的故障分离结果，为后续的故障机理分析奠定基础。

(2)本发明通过故障分离的T²和SPE统计量贡献图来对飞行控制系统故障检测、故障分离进行分析，可以获得各种故障情况下受故障变量影响较大的其它变量，其分析结果为无人机飞行控制系统进行故障容错控制提供一种参考依据。

附图说明

图1为无人机飞行控制系统框图；

图2为本发明提出的无人机飞行控制系统故障检测与故障分离的数据驱动方法流程图；

图3为本发明所采用的无人机纵向飞行控制系统闭环仿真平台Simulink框图；

图4为不同动态阶数下的平行分析结果图；

图5为特征值及累积方差贡献率结果图，其中：图5(a)为特征值结果图；图5(b)为累积方差贡献率结果图；

图6为无人机飞行控制系统T²和SPE统计量故障检测图，其中：图6a为角速率陀螺偏差故障T²统计量故障检测；图6b为角速率陀螺偏差故障SPE统计量故障检测；图6c为油门杆部分失效故障T²统计量故障检测；图6d为油门杆部分失效故障SPE统计量故障检测；

图7为无人机飞行控制系统T²和SPE统计量贡献图；

图8为无人机飞行控制系统T²和SPE统计量累积贡献图。

具体实施方式

本发明无人机飞行控制系统框图如图1所示，主要由控制器、执行器、无人机机体和传感器组成，控制指令通过控制器作用于执行器，通过控制输入引起无人机机体位置、速度和姿态等状态变量的变化，经过传感器进行量测并将量测变量反馈给控制器，控制器再根据控制指令和量测变量产生控制输入。无人机飞行控制系统故障检测与故障分离的数据驱动方法流程图如图2所示，具体方法如下：

步骤2：将步骤1采集到的数据整理为矩阵形式，假设有N次测量，将每次采样点获得的无人机输入输出数据整理为向量x，假设x中包含m个变量，那么形成的数据矩阵为X：

X＝[x₁，x₂，…,x_N]∈R^m×N

然后，利用移动均值法对数据进行平滑预处理：

x_i＝(x_i+x_i-1+x_i-2+…+x_i-9)/10

步骤3：由步骤2得到平滑预处理的原始量测数据后，为了获得更多具有相关性的数据信息，需要掌握当前时刻数据与过去s个时刻数据的关系，s为采样数据间最大持续影响的采样次数，即动态阶数。在这里，我们通过以下的平行分析方法来确定动态阶数s：

a.令s取值为0；

b.取正常飞行的数据矩阵X并对其进行归一化，使其均值为零并具有单位方差，然后计算归一化数据矩阵X的协方差阵：

c.对S进行奇异值分解，以获得特征值组成的对角矩阵Λ＝diag(λ₁,…λ_m)以及相应的特征向量组成的矩阵P，其中λ₁≥…≥λ_m≥0。将P中前k个特征向量组成的矩阵定义为载荷矩阵P_pc∈R^m×k，剩余的特征向量组成的矩阵定义为载荷矩阵P_res∈R^m×(m-k)。对于一个采样向量x∈R^m可以向载荷向量进行投影，获得x在每个特征向量方向上的得分t＝P_pc ^Tx

d.令j＝m×(s+1)并取s阶时间滞后情况下总的线性关系的个数r(s)＝0；

e.确定第j个主元是否代表着一个线性关系，是则继续进行下一步，否则跳转到步骤g；

f.令j＝j-1，取r(s)＝r(s)+1，重复步骤e；

g.计算第s阶包含的线性关系的个数

h.如果r_new(s)<0，跳转到步骤j，否则继续进行下一步；

i.令s＝s+1，把原始量测数据X前面s个时刻的量测数据的每个变量当成新的增广阵的变量，建立数据增广矩阵X_a，代替原始数据矩阵X跳转到步骤b；

j.停止计算，当前s即为所求的动态阶数。

步骤4：根据步骤3确定的动态阶数s，把前面s个时刻的量测数据的每个变量当成新的增广阵的变量，建立数据增广矩阵X_a。对该数据增广矩阵X_a进行如下动态主元分析，以建立动态主元分析模型：

a.构建数据增广矩阵X_a。

其中N₁≥k+N。

b.对矩阵X_a进行归一化处理，得到标准化后的数据矩阵X_S，其均值为零，方差为1。

c.对X_S的协方差矩阵进行奇异值分解，以获得特征值组成的对角矩阵Λ以及相应的特征向量组成的矩阵P：

d.利用累积方差贡献率的方法来确定主元个数并获得负荷矩阵P。当累积方差贡献率CPV达到90％时，计算出来的主元可以综合系统的有效信息，此时的k即为主元个数：

Λ_pc＝diag(λ₁,…,λ_k),Λ_res＝diag(λ_k+1,…,λ_m)；P＝[P_pc P_res],P_pc∈R^m×k,P_res∈R^m ^×(m-k)

步骤5：根据步骤4建立的动态主元模型，通过计算新数据在主元子空间变化的程度以及在残差子空间变化的程度来检测系统是处于正常状态还是出现故障。为了提高故障检测与分离的准确率，本发明配合使用Hotelling的T²统计量和平方预测误差SPE统计量两种故障检测方法，利用飞机正常飞行数据，计算T²和SPE两种监测统计量的阈值。其中：

T²统计量测量了新数据在主元子空间的变化情况，表示每个采样在变化趋势和幅值上偏离主元模型的程度，用来检验系统新的量测数据是否依然在主元确定的控制域内，因此T²统计量表征了一种模型内部变化的一种测度。它的定义为：

T²＝x^TP_pcΛ^-1 _pcP_pc ^Tx

正常情况下的数据服从一个多元正态分布，因此T²统计量服从一个F分布：

那么就可以得到置信度为α的T²统计量控制上限即阈值为：

平方预测误差(SPE)又称为Q统计量，测量了采样向量在残差子空间上投影变化的大小，表示原始数据中没有被主元模型所解释的变化，是模型外部的一种测度。它的定义为：

令那么c可以近似看作服从标准正态分布：

对于置信度α，其分位数为C_α，那么相应的SPE统计量控制上限即阈值为：

步骤6：根据步骤4建立的动态主元模型，可知当故障发生后，通过故障飞行数据可以计算故障发生后各采样点的T²和SPE两种统计量：

T²＝t_newΛ^-1t_new ^T＝x_new ^Tp_pcΛ_pc ^-1p_pc ^Tx_new

SPE＝e_newe_new ^T＝x_new(I-P_pcP_pc ^T)x_new ^T

其中，x_new为当前采样点的测量数据向量。

此时，便可以通过与步骤5得到的阈值比较，以下面的逻辑规则进行故障检测：

T²统计量和SPE统计可以联合起来进行故障检测，当故障发生后，T²统计量或者SPE统计量至少有一个会超出阈值，这就意味着有故障发生，从而实现对故障的及时检测。

步骤7：在步骤6完成故障检测的基础上，为了对故障进行分离，在故障发生后的每一个时刻，计算出每个变量对T²和SPE统计量的贡献值。在采样时刻j，变量i对T²和SPE的贡献值的计算公式如下：

T² _ij＝x_j ^TM^(1/2)ξ_iξ_i ^TM^(1/2)x_j，M＝P_pcΛ_pc ^-1P_pc ^T

由于测量数据矩阵的增广矩阵X包含时间滞后s，因此m×(s+1)个贡献值属于原始的m个变量，画出原始m个变量对故障的贡献图更具有实际意义。所以将属于同一个变量的s+1个贡献值累加起来作为这个原始变量的贡献值。

SPE^ori _i,j＝SPE_i,j+SPE_(i+m),j+…+SPE_(i+s*m),j,i＝1,…m

因为无人机飞行控制系统是动态的，T²和SPE统计量的大小在每次采样都会发生变化。在每一个时刻各个变量的贡献值可以通过除以当前时刻的相应统计量进行归一化，从而获得每个变量对统计量的贡献率。

在j时刻变量i对T²统计量的贡献率为：

在j时刻变量i对SPE统计量的贡献率为：

CR_SPE_i,j＝SPE^ori _i,j/SPE_j

最后，将贡献率以变量为纵坐标，时间为横坐标作图，即可得到贡献图；

步骤8：将故障发生后一段时间内各变量对统计量的贡献值进行累加，获得各变量对统计量的累积贡献率，获得累积贡献图。

为了更好的分离故障变量，故障发生后一段时间内各变量对统计量的累积贡献率也可以通过将各个变量的贡献率分别进行累加得到：

其中，t_f是故障发生的时刻，n_f是故障发生后总共的采样数，是变量i对T²统计量的累积贡献率，SPE^CCR _i是变量i对SPE统计量的累积贡献率；

步骤9：在一般情况下，故障变量对统计量的贡献率大于其它变量的贡献率，而故障变量的累积贡献率更是远大于其它变量的累计贡献率。那么在步骤7和步骤8的基础上，即可通过贡献图和累积贡献图确定故障变量，实现故障的分离

实施例

本发明以一种无人机纵向飞行控制系统为例对数据驱动的故障检测与分离进行分析，首先对该无人机纵向飞行控制系统进行说明：

输入变量为升降舵偏转量δ_e和油门杆偏转量δ_p，分别用来实现对速度通道和高度通道的控制；传感器配置为攻角传感器、空速管、俯仰陀螺、垂直陀螺和气压式高度表，分别量测无人机运动过程中攻角α、空速V、俯仰角速率q、俯仰角θ及高度H，并假设系统传感器无硬件冗余，即一个物理量仅由一个传感器量测。

考虑到无人机在飞行过程中不可避免受到大气紊流、阵风等扰动因素的影响，建立如下变化风场中的无人机纵向运动方程为：

其中：m为机体质量，g为重力加速度；α为攻角；P为发动机推力，Q为空气阻力，Y为升力，M_z为俯仰力矩；I_z为飞机沿机体系俯仰轴的转动惯量；分别为水平风速与垂直风速沿前进方向的风速梯度，w_h表示垂直风速的大小。

考虑风场中最常见的两种形式为大气紊流和突变，因此本文对这两种风场形式进行仿真。在大气扰动模块的设计当中，分别采用紊流Dryden模型和突变全波长(1-consine)模型对进行仿真。

Dryden模型利用白噪声通过成形滤波器仿真大气紊流。通过对紊流频谱进行分解，可以得到生成前向紊流成形滤波器和垂向紊流速度成形滤波器的参数如下：

其中，L_x,L_y是前向和垂向的紊流尺度，σ_x,σ_y为前向和垂向的紊流强度，V为空速。

全波长(1-cos)型对突变建模，表示为：

其中，V_w为突变强度，d_m为突变尺度，x为突变在空间域的长度。

取系统状态变量x(t)＝[V α q θ H]^T，控制输入量u(t)＝[δ_e δ_p]^T，量测变量y(t)＝[V_m α_m q_m θ_m H_m]^T，风干扰f(t)表示传感器或作动器故障，u(t)得到无人机飞行控制系统模型如下：

其中，w(t)、v(t)分别为系统噪声与量测噪声，F(x(t))、B(x(t))、G(x(t))、C(t)的计算方法如下：

其中，S_w为机翼面积，ρ为空气密度，为平均气动弦长，K表示发动机推力与油门杆偏转量为比例系数；C_x0、C_y0、均为气动导数，由无人机在参考飞行状态下风动实验的结果得到；m_y0、表示俯仰力矩导数，由无人机尺寸等参数确定。

根据如上所述纵向模型搭建无人机纵向飞行控制系统闭环仿真平台如图3所示，用于产生正常和故障情况下的飞行数据。以角速率陀螺偏差故障和油门杆部分失效故障为例分别进行仿真实验，对两种故障情况下的故障检测与分离情况进行分析。

下面对该无人机系统的故障问题，采用数据驱动的方法进行故障检测与分离：

步骤1：首先设定仿真飞行时间300s，采样周期T＝0.01s，前200s飞行数据作为训练数据，用于平行分析、建立主元模型以及计算阈值，后100s飞行数据作为测试数据，用于进行故障检测与分离。在250s时注入故障；系统初始状态和目标状态一致，这样使无人机处于稳态飞行状态；系统噪声、量测噪声以及风扰动均设为高斯白噪声。

实验采集的无人机输入输出数据列在表1中：

表1无人机的输入和输出

变量编号	变量类型	变量描述
			1	传感器输出	攻角α
2	传感器输出	空速V
			3	传感器输出	高度H
4	传感器输出	俯仰角θ
			5	传感器输出	俯仰角速率q
6	执行器输入	油门杆偏转量δ_p
			7	执行器输入	升降舵偏转量δ_e

对于传感器故障，以角速率陀螺偏差故障为例：

在0至250s正常飞行，250s至300s发生角速率陀螺偏差故障，故障形式及故障注入时间如下式所示：

对于执行机构故障，以油门杆部分失效故障为例：

在0至250s正常飞行，250s至300s发生油门杆部分失效故障，故障形式及故障注入时间如下式所示：

步骤2：将数据整理为矩阵形式，利用移动均值法对数据进行平滑预处理。其中每个采样点的测量数据整理为向量x＝[α,V,H,θ,q,δ_p,δ_e]^T。

角速率陀螺偏差故障下的平滑数据矩阵为：

油门杆部分失效故障下的平滑数据矩阵为：

由于采用同一仿真模型，前200s的训练数据是相同的，后100s的测试数据因故障设置不同而不同。故而两种故障情况下的平行分析、主元模型以及阈值均是一致的，下述步骤3、4、5中不再重复。

步骤3：通过平行分析确定系统的动态阶数s。

平行分析的结果如图4所示，小于1的特征值的个数被看做线性关系的个数，那么图中低于红线阈值的特征值的个数就是所有线性关系的个数。当s＝0时，r(0)＝4，r_new(0)＝4，说明该系统具有4个零阶线性关系；当s＝1时，r(1)＝10，r_new(1)＝2，说明该系统具有2个一阶线性关系，；当s＝2时，r_new(2)＝1，能够得到一个二阶线性关系；当s＝3时，r_new(3)＝0，不能得到新的动态线性关系，因此取s＝2比较合适。

步骤4：根据动态阶数s构建数据增广矩阵X_a，建立正常飞行数据的动态主元模型。

由s＝2构建增广矩阵如下：

对X_a归一化后得到X_S，并对X_S的协方差矩阵进行奇异值分解得到：

特征值及累积方差贡献率如图5所示，可见当累积方差贡献率达到90％时，主元个数k＝8。由此得到主元模型为：

步骤5：计算T²和SPE两种监测统计量的阈值。

将上述数据带入T²和SPE两种监测统计量的阈值计算公式中，得到J_th,SPE＝5.73。

步骤6：通过故障数据可以计算故障发生后各采样点的T²和SPE两种统计量，并结合步骤5得到的阈值，对系统进行故障检测。

角速率陀螺偏差故障以及油门杆部分失效故障下的T²和SPE两种统计量的故障检测结果分别如图6a、6b、6c、6d所示，可见在两种故障发生后，T²统计量和SPE统计量均迅速超过阈值，实现了对故障的快速检测，检测正确率列表如下：

表2故障检测正确率

故障	T²统计量	SPE统计量
			角速率陀螺偏差故障	100％	98.22％
油门杆部分失效故障	96.38％	97.88％

由表可知，无人机故障检测的数据驱动方法是准确可行的。

步骤7：获取故障发生后各变量对T²和SPE两种统计量的贡献值，并通过除以当前时刻的相应统计量进行归一化，从而获得在各个采样点各个变量对统计量的贡献率，最终得到角速率陀螺偏差故障以及油门杆部分失效故障下的T²和SPE两种统计量贡献图如图7所示。

步骤8：将故障发生后一段时间内各变量对统计量的贡献值进行累加，获得各变量对统计量的累积贡献率，得到角速率陀螺偏差故障以及油门杆部分失效故障下的T²和SPE两种统计量累积贡献图如图8所示。

步骤9：通过步骤7、8得到的贡献图和累积贡献图进行故障分离，以确定故障变量。分离结果如下表所示：

表3故障分离结果

从表中可以看出，角速率陀螺偏差故障发生后，由T²统计量与SPE统计量判断得到的故障分离结果均只有一个，并且不同，根据步骤9中判别准则(4)：当T²统计量累积贡献图以及SPE统计量累积贡献图中故障变量不一致时，由于SPE统计量很难反应系统的动态特性，基于SPE统计量贡献图的故障分离结果大部分都不准确，此时以T²统计量贡献图的故障分离结果为准。可以认定故障分离结果为故障5，与故障源一致。具体分析如下：

角速率陀螺作为传感器，在故障发生后对无人机飞行控制系统的影响很小，故障难以在系统内部传播，因此由T²统计量与SPE统计量判断得到的故障分离结果均只有一个。而T²统计量用来检验系统新的量测数据是否依然在主元确定的控制域内，是模型内部变化的一种测度，能够较好地反应系统的动态特性；SPE统计量用来检验原始数据中没有被主元模型所解释的变化，是模型外部的一种测度，难以反应系统的动态特性。而无人机的动态特性相当显著，因此T²统计量的分离结果可信度较高，而SPE统计量的分离结果可信度较低，可以判定故障源为故障5。

油门杆部分失效故障发生后，由T²统计量与SPE统计量判断得到的故障分离结果均有分别为一个、两个，但只有一个相同，根据步骤9中判别准则(3)：当T²统计量累积贡献图以及SPE统计量累积贡献图中故障变量一致时，即可判定这个变量为根本的故障变量。可以认定故障分离结果为故障6，与故障源一致。具体分析如下：

油门杆作为执行器，在故障发生后对无人机飞行控制系统的影响较大，故障得以在系统内部传播，因此由图8可见T²统计量得到的故障分离结果虽然可以确定为故障6，但其他变量的故障累积贡献率相对角速率陀螺故障时依然较大，而SPE统计量判断得到的故障分离结果更是有故障3、故障6两个。但是依照准则(3)的一致原则，仍可以判定故障源为故障6。

由此可知，本发明所提出的无人机故障分离的数据驱动方法是有效的。

Claims

1.一种无人机飞行控制系统故障检测与故障分离的数据驱动方法，其特征在于，包括下列步骤：

无人机执行机构包含升降舵、方向舵、副翼和油门杆这四个执行机构，其舵偏角分别用δ_e、δ_r、δ_a和δ_p表示，因此无人机输入数据包含以上4个变量；

无人机传感器主要包括攻角传感器、侧滑角传感器、空速表、高度表、垂直陀螺以及角速率陀螺，输出数据依次为攻角α、侧滑角β、空速V、高度H，俯仰角θ、偏航角ψ、滚转角φ、滚转角速度p、俯仰角速度q、偏航角速度r，因此无人机输出数据包含以上10个变量；

步骤3：由步骤2得到平滑预处理的原始量测数据后，为了获得更多具有相关性的数据信息，需要掌握当前时刻数据与过去s个时刻数据的关系，s为采样数据间最大持续影响的采样次数，即动态阶数，在这里，通过平行分析法来确定动态阶数s；

步骤4：根据步骤3确定的动态阶数s，把前面s个时刻的量测数据的每个变量当成新的增广阵的变量，建立数据增广矩阵X_a，对该数据增广矩阵X_a进行动态主元分析，建立动态主元分析模型，为后续的故障检测与分离提供依据；

步骤5：根据步骤4建立的动态主元模型，通过计算新数据在主元子空间变化的程度以及在残差子空间变化的程度来检测系统是处于正常状态还是出现故障，为了提高故障检测与分离的准确率，采用将Hotelling的T²统计量和平方预测误差SPE统计量两种故障检测方法配合使用的方法，利用飞机正常飞行数据，计算两种监测统计量的阈值，其中：

T²统计量表示每个采样在变化趋势和幅值上偏离主元模型的程度，用来检验系统新的量测数据是否依然在主元确定的控制域内，是模型内部变化的一种测度；

SPE统计量表示采样向量在残差子空间上投影变化的大小，用来检验原始数据中没有被主元模型所解释的变化，是模型外部变化的一种测度；

步骤6：根据步骤4建立的动态主元模型，计算故障发生后各采样点的T²和SPE两种统计量，并分别与步骤5得到的阈值比较，当且仅当T²统计量或者SPE统计量均在相应的阈值范围内时，认为系统正常，反之故障，从而实现了对故障的检测；

步骤7：在步骤6完成故障检测的基础上，为了对故障进行分离，采用贡献图的方法进行分析，具体方法是：在故障发生后的每一个时刻，计算出每个变量对T²和SPE统计量的贡献值，并通过除以当前时刻的相应统计量进行归一化，从而获得在各个采样点各个变量对统计量的贡献率，最终将其以变量为纵坐标，时间为横坐标作图，即可得到贡献图；

步骤9：根据贡献图和累积贡献图理论：在一般情况下，故障变量对统计量的贡献率大于其它变量的贡献率，而故障变量的累积贡献率更是远大于其它变量的累计贡献率，则可在步骤7和步骤8的基础上，通过贡献图和累积贡献图确定故障变量，实现故障的分离。

2.根据权利要求1所述的无人机飞行控制系统故障检测与故障分离的数据驱动方法，其特征在于：步骤2中所述平滑预处理的移动均值法为：

x_i＝(x_i+x_i-1+x_i-2+…+x_i-9)/10

3.根据权利要求1所述的无人机飞行控制系统故障检测与故障分离的数据驱动方法，其特征在于：步骤3所述的进行平行分析的具体方法为：

对S进行奇异值分解，以获得特征值组成的对角矩阵Λ＝diag(λ₁,…λ_m)以及相应的特征向量组成的矩阵P，其中λ₁≥…≥λ_m≥0，将P中前k个特征向量组成的矩阵定义为载荷矩阵P_pc∈R^m×k，剩余的特征向量组成的矩阵定义为残差矩阵P_res∈R^m×(m-k)，对于一个采样向量x∈R^m可以向载荷向量进行投影，获得x在每个特征向量方向上的得分t＝P_pc ^Tx；

然后，取s阶时间滞后情况下总的线性关系的个数r(s)＝0，判断第j＝m×(s+1)个主元是否代表着一个线性关系，是则令j＝j-1，r(s)＝r(s)+1，重复该判断；否，则进行下一步；

4.根据权利要求1所述的无人机飞行控制系统故障检测与故障分离的数据驱动方法，其特征在于：步骤4所述的建立数据的动态主元模型具体步骤为：

首先，构建数据增广矩阵X_a，

其中N₁≥k+N，

其次，对矩阵X_a进行归一化处理，得到均值为零，方差为1的标准化后的数据矩阵X_S；

Λ＝diag(λ₁,…λ_m)，λ₁≥…≥λ_m≥0

最后，利用累积方差贡献率的方法来确定主元个数并获得负荷矩阵P，当累积方差贡献率CPV达到90％时，计算出来的主元可以综合系统的有效信息，此时的k即为主元个数，

<mrow> <mi>C</mi> <mi>P</mi> <mi>V</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>100</mn> <mo>*</mo> <mo>&lsqb;</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>k</mi> </munderover> <msub> <mi>&lambda;</mi> <mi>i</mi> </msub> <mo>/</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msub> <mi>&lambda;</mi> <mi>i</mi> </msub> <mo>&rsqb;</mo> <mi>%</mi> </mrow>

5.根据权利要求1所述的无人机飞行控制系统故障检测与故障分离的数据驱动方法，其特征在于：步骤5所述的计算T²和SPE两种监测统计量阈值的具体求解方法为：

置信度为α时T²统计量阈值求解公式如下：

<mrow> <msub> <mi>J</mi> <mrow> <mi>t</mi> <mi>h</mi> <mo>,</mo> <msup> <mi>T</mi> <mn>2</mn> </msup> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <mi>k</mi> <mrow> <mo>(</mo> <mi>N</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <mi>N</mi> <mo>-</mo> <mi>k</mi> </mrow> </mfrac> <msub> <mi>F</mi> <mi>&alpha;</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>,</mo> <mi>N</mi> <mo>-</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow>

其中，F_α代表置信度为α时的F分布；

置信度为α时SPE统计量阈值求解公式如下：

<mrow> <msub> <mi>J</mi> <mrow> <mi>t</mi> <mi>h</mi> <mo>,</mo> <mi>S</mi> <mi>P</mi> <mi>E</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>Q</mi> <mi>&alpha;</mi> </msub> <mo>=</mo> <msub> <mi>&theta;</mi> <mn>1</mn> </msub> <msup> <mrow> <mo>&lsqb;</mo> <mfrac> <mrow> <msub> <mi>C</mi> <mi>&alpha;</mi> </msub> <msqrt> <mrow> <mn>2</mn> <msub> <mi>&theta;</mi> <mn>2</mn> </msub> <msup> <msub> <mi>h</mi> <mn>0</mn> </msub> <mn>2</mn> </msup> </mrow> </msqrt> </mrow> <msub> <mi>&theta;</mi> <mn>1</mn> </msub> </mfrac> <mo>+</mo> <mn>1</mn> <mo>+</mo> <mfrac> <mrow> <msub> <mi>&theta;</mi> <mn>2</mn> </msub> <msub> <mi>h</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <msub> <mi>h</mi> <mn>0</mn> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <msup> <msub> <mi>&theta;</mi> <mn>1</mn> </msub> <mn>2</mn> </msup> </mrow> </mfrac> <mo>&rsqb;</mo> </mrow> <mfrac> <mn>1</mn> <msub> <mi>h</mi> <mn>0</mn> </msub> </mfrac> </msup> </mrow>

其中，C_α为置信度为α时的标准正态分布的分位数，

6.根据权利要求1所述的无人机飞行控制系统故障检测与故障分离的数据驱动方法，其特征在于：步骤6所述的故障数据各采样点T²和SPE两种统计量具体求解方法为：

T²＝t_newΛ^-1t_new ^T＝x_new ^Tp_pcΛ_pc ^-1p_pc ^Tx_new

SPE＝e_newe_new ^T＝x_new(I-P_pcP_pc ^T)x_new ^T

其中，x_new为当前采样点的测量数据向量；

此时，便可以通过以下逻辑规则进行故障检测：

SPE＞J_th,SPE或故障

SPE<J_th,SPE与正常。

7.根据权利要求1所述的无人机飞行控制系统故障检测与故障分离的数据驱动方法，其特征在于：步骤7所述的获取贡献图的具体步骤为：

首先，在采样时刻j，计算变量i对T²和SPE的贡献值：

T² _ij＝x_j ^TM^(1/2)ξ_iξ_i ^TM^(1/2)x_j,M＝P_pcΛ_pc ^-1P_pc ^T

其中，x_j是采样时刻j的测量数据，ξ_i是单位矩阵的第i列，e_ij是采样时刻j变量i的残差，x_ij是采样时刻j变量i的测量值，是采样时刻j变量i的估计值；

然后，计算变量i对T²和SPE的累积贡献值：

SPE^ori _i,j＝SPE_i,j+SPE_(i+m),j+…+SPE_(i+s*m),j,i＝1,…m

最后，采用归一化的方法，获得每个变量对统计量的贡献率；

无人机飞行控制系统是动态的，T²和SPE统计量的大小在每次采样都会发生变化，在每一个时刻各个变量的贡献值可以通过除以当前时刻的相应统计量进行归一化，从而获得每个变量对统计量的贡献率；

在j时刻变量i对T²统计量的贡献率为：

在j时刻变量i对SPE统计量的贡献率为：

CR_SPE_i,j＝SPE^ori _i,j/SPE_j。

8.根据权利要求7所述的无人机飞行控制系统故障检测与故障分离的数据驱动方法，其特征在于：步骤8所述的获取累积贡献图的具体步骤为：

9.根据权利要求1所述的无人机飞行控制系统故障检测与故障分离的数据驱动方法，其特征在于：步骤9所述的获取故障变量的具体步骤为：

<mrow> <msub> <mi>J</mi> <msup> <mi>T</mi> <mn>2</mn> </msup> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>m</mi> </mfrac> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msub> <msup> <mi>T</mi> <msup> <mn>2</mn> <mrow> <mi>C</mi> <mi>C</mi> <mi>R</mi> </mrow> </msup> </msup> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>J</mi> <mrow> <mi>S</mi> <mi>P</mi> <mi>E</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>m</mi> </mfrac> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msub> <msup> <mi>SPE</mi> <mrow> <mi>C</mi> <mi>C</mi> <mi>R</mi> </mrow> </msup> <mi>i</mi> </msub> </mrow>

然后，在T²统计量累积贡献图以及SPE统计量累积贡献图中，将超过控制线的变量作为可能的故障变量，下面将对故障分离结果进行分类讨论，给出如下故障分离判别准则；

(1)当可能的故障变量为单个变量时，即可判定这个单一变量为根本的故障变量；

(2)当可能的故障变量为多个变量时，结合贡献图中故障发生初始时刻多个可能故障变量对统计量贡献率的变化情况来推断根本的故障变量；

(3)当T²统计量累积贡献图以及SPE统计量累积贡献图中故障变量一致时，即可判定这个变量为根本的故障变量；