CN106055885A - 基于过采样投影近似基追踪无人机飞行数据异常检测方法 - Google Patents

Abstract

本发明提供了一种基于过采样投影近似基追踪无人机飞行数据异常检测方法。本发明引入过采样和信号处理中的信号子空间投影近似思想，通过对过采样后数据投影近似子空间基向量方向的估计和追踪匹配，检测数据流中异常，同时利用子空间方向对飞行模式切换不敏感特点，抑制飞行模式对异常检测结果的影响。从而以此为基础提出了一种无人机飞行数据在线异常检测方法的框架。所提出的检测框架，在飞行数据在线异常检测中，消耗同等μs级计算时间下，较基于Online MD、BN、CCA和KOAD四种在线异常检测方法，误检率FPR降低53.82％以上，检测准确率AUC分数提高5.28％以上，达到0.9836，接近理论值1。表明本发明的方法可有效解决飞行数据在线异常检测问题。

Description

基于过采样投影近似基追踪无人机飞行数据异常检测方法

技术领域

本发明涉及一种基于过采样投影近似基追踪的无人机飞行数据异常检测方法，属于无人机飞行数据异常检测方法技术领域。

背景技术

无人机(Unmanned Aerial Vehicle,UAV)的部署和应用，需以保证对自身、其它飞行器、地面设施及人员有足够的安全性为前提。传感器是无人机感知自身飞行状态和与外界环境交互的关键部件，但容易发生故障，从而导致飞行失控等安全事故。同时，受限于无人机体积、重量和功耗约束，很难在传感器***中采用物理冗余技术，因而无人机在可靠性和安全性方面相比有人机存在较大差距。为了及时隔离故障并为后续任务的重新规划等应急响应提供决策依据，无人机传感器故障检测技术成为不可回避的挑战。

无人机传感器，如气压高度计、GPS传感器和姿态指示计等，发生性能退化或失效时，会在对应飞行数据中表现异常。异常是远离其他数据而疑为不同机制产生的数据。异常检测是一种发现不符合期望行为的数据或模式的技术，已广泛应用于网络流量检测和***状态监测等领域，近年来逐渐成为无人机故障检测的一种新颖方法。相比传统基于模型或知识的方法，异常检测不需要复杂的精确动力学模型建模和验证，并且适合广泛存在的多样和未知故障类型。

由于缺少可以依靠的其他感知手段，故通过发现飞行数据内部存在关联性数据的相关性破裂，已成为现有无人机飞行数据异常检测的有效思路。飞行数据间的关联性主要是控制与反馈数据关联和同类数据间关联。此外，飞行数据呈现多维数据流和随飞行模式动态变化的特性，为了满足机载处理的实时性要求，数据流的异常检测方法应具备计算轻量级的特点。因此，目前无人机异常检测领域的学者主要研究针对数据流的轻量级在线异常检测方法，以发现数据流的相关性破裂，并判断其中每个输入向量的异常程度。根据有无使用带标签的训练数据，现有方法可分为有监督和无监督方法两大类。

有监督方法通过有监督训练建立的数据模型都存在一定局限性，如限定为横向动力学或假设爬升和下降的概率相等，以及依赖可靠的先验传感器故障概率。但是真实的无人机飞行包含完整的横向和纵向动力学特性，并且不同飞行模式的切换概率和传感器真实故障概率很难获得。所以，目前有监督训练的数据模型，随着飞行模式切换，容易发生失效，从而导致异常检测准确率降低。

无监督方法相比于有监督方法，具有不依赖离线训练的数据模型的优势。由于传感器输出的正常数据与相关性破裂后的异常数据在统计意义下存在分布差异，因而可以用基于距离的无监督方法度量差异程度，判断出数据流的异常。然而，马氏距离数值上计算不稳定，此外基于距离的方法对于飞行模式切换下的多分布数据，容易误检，导致异常检测的准确率降低，引发频繁误警问题。

可见，现有方法尚未有效地解决飞行模式切换导致的分析方法失配和异常检测准确率低问题。异常检测方法本质上需要具备能刻画原始数据结构和评价该结构差异程度的能力。数据子空间是原始数据空间向低维映射得到的子空间。子空间方向是其典型特征之一，为原始数据向子空间投影的基向量。当原始数据空间出现异常数据时，数据子空间方向会随之发生角度变化。此外，飞行模式切换时保持着相关性的数据的变化模式基本一致，其数据子空间方向也基本保持相同角度，因而飞行模式切换对未发生相关性改变的数据的子空间方向影响较小。然而，数据流应用中，异常和正常数据子空间方向较小的差异，导致两者可分性降低。并且传统子空间基向量的求解方法复杂度高，并不适合机载在线应用。

发明内容

本发明的目的是为了解决上述现有技术存在的问题，即针对飞行模式切换导致的无人机飞行数据在线异常检测准确率低的问题。进而提供一种基于过采样投影近似基追踪的无人机飞行数据异常检测方法。

本发明的目的是通过以下技术方案实现的：

一种基于过采样投影近似基追踪无人机飞行数据异常检测方法，步骤如下：

(一)定义

定义1：飞行数据流模型A是以一定采样频率，连续不断，随时间增长且只能读取一次的p维时间序列；

$A=\{{\stackrel{\→}{a}}_1,{\stackrel{\→}{a}}_2,\mathrm{...},{\stackrel{\→}{a}}_t,{\stackrel{\→}{a}}_{t+1},\mathrm{...}\}---\left(1\right)$

其中，是飞行数据流模型A在t时刻的输入向量，p为所监控参数总数，也是向量维度，a_t,j∈R表示参数j在时刻t的值，j∈[1,p]分别对应所监控的不同飞行参数，如气压高度、气压升降速度、垂向速度、俯仰角、组合高度、GPS垂向速度、GPS高度；

定义2：观察矩阵H_p×m，无限的数据流无法全部存储和处理，可通过保持固定窗口长度的滑窗存储部分数据，并作归一化以统一数据尺度，从而建立的观察矩阵H_p×m：

$H_{p\×m}=\left\{{\stackrel{\→}{a}}_{t-m}\mathrm{......}{\stackrel{\→}{a}}_{t-1}\right\}---\left(2\right)$

其中，H_p×m维度为p×m，列数m为滑窗所记录的过去时间点长度，行数p为所记录的飞行参数数量；

定义3：观察矩阵子集H_r×m和飞行数据流子集H_r×m为从观察矩阵H_p×m中抽取的相关的观察矩阵子集，数据子集如高度类数据子集，垂向速度类数据子集，H_r×m中行数r为子集维度且r＜p，列数m为所记录的过去时间点长度，是观察矩阵子集对应的输入向量子集，r为子集维度且r＜p；

定义4：飞行数据流异常检测，根据观察矩阵子集H_r×m，判断t时刻输入向量子集是否为异常；

$Anomaly\left(t\right)=(H_{r\×m},{\stackrel{\→}{a}}_t(r\left)\right)---\left(3\right)$

(二)异常检测总体框架

无人机飞行数据在线异常检测框架，框架包括滑窗与归一化、观测矩阵子集抽取和过采样投影近似基追踪三个部分，输入为多维的飞行数据流输出为时刻t的异常报警；

(1)滑窗与归一化是为了建立定义2所示的观测矩阵，其中，滑窗从无限的飞行数据流中存储并更新预设宽度的飞行数据；归一化主要考虑到原始飞行数据中，不同量纲的参数相互间不存在可比性；

(2)飞行数据中相关的参数组成的数据子集是本方法的检测对象，需要从观测矩阵中抽取相关子集；

(3)过采样投影近似基追踪在线异常检测算法，当原始数据空间出现异常数据时，数据子空间方向会随之发生角度变化，OSPABP实质上通过度量子空间方向变化程度实现异常检测；

(三)滑窗与归一化

为消除不同参数的量纲并考虑数据的统计分布，将输入飞行数据输入向量利用式(4)进行Z-score变换，从而将滑窗中的数据归一化为均值为0、方差为1的数据，输入向量的Z-score变化为

$Z\left({\stackrel{\→}{a}}_t\right)=\{Z({\stackrel{\→}{a}}_{t,1},H_1^T),\mathrm{...},Z({\stackrel{\→}{a}}_{t,p},H_p^T)\}---\left(4\right)$

其中为向量均值，σ_x为向量方差，a_t,j为飞行数据向量中的参数，为观测矩阵H_p×m中参数j所对应的过去m点数据；

(四)观测矩阵子集抽取

从观察矩阵中抽取相关的子集H_r×m，作为异常检测算法的输入，观察矩阵子集中相关飞行参数的选取，可对历史的飞行数据进行相关性分析确定，选取相互间具有较大相关系数的参数构成子集，如相关系数大于0.997；

(五)过采样投影近似基追踪

将原始的多维飞行数据观测矩阵，处理为多个相关的子集矩阵H_r×m后，可对观测矩阵中的子集矩阵采用过采样投影近似基追踪方法跟踪和匹配t时刻输入数据带来的数据子空间方向变化，从而检测时刻t的飞行数据是否存在异常；

OSPABP方法包括两个步骤：第一步，使用过采样PCA放大子集输入向量异常效果，并实现原始数据向子空间转换；第二步，采用投影近似子空间估计方法，近似求解子空间方向，并追踪和匹配子空间方向变化，从而实现飞行数据流异常检测；

(六)过采样投影近似基追踪方法

采用投影近似子空间估计方法，以递推投影近似形式估计子空间，使计算复杂度下降两个数量级，简化子空间基向量求解，从而满足机载处理计算实时性要求；

(七)过采样PCA

假设过采样t时刻的目标实例次，则相应过采样PCA方程为：

${\mathrm{\Σ}}_{\tilde{A}}{\tilde{u}}_t=\mathrm{\λ}{\tilde{u}}_t---\left(5\right)$

其中，过采样后数据矩阵为的均值为μ，协方差矩阵为为特征向量；为原始数据矩阵，每行x_i代表i时刻的p维空间数据实例，n为实例数量，x_t为t时刻新的目标实例；协方差矩阵随着每一个新的目标实例x_t而变化；

标准PCA可看成判断投影数据有最大方差的子空间，可通过最小化数据重构误差J(U)求解；

$\underset{U\∈R^{p\×k},\left|\right|U\left|\right|=I}{\min }J\left(U\right)=\underset{i=1}{\overset{n}{\Σ}}\left|\right|{\stackrel{\‾}{x}}_i-{\mathrm{UU}}^T{\stackrel{\‾}{x}}_i|{|}^2---\left(6\right)$

其中U为包含k个显著特征向量的子空间矩阵，为去均值的第i个原始数据，为去均值的第i个原始数据向特征向量子空间投影，是从投影子空间重构估计得到的第i个用估计为前一个时刻估计的特征向量空间，从而，式(6)转化为：

$\underset{U\∈R^{p\×k},\left|\right|U\left|\right|=I}{\min }J\left(U_t\right)=\underset{i=1}{\overset{t}{\Σ}}\left|\right|{\stackrel{\‾}{x}}_i-U_i{y}_i|{|}^2---\left(7\right)$

同样，如果过采样目标实例次，在线过采样PCA的最小化数据重构误差为：

$\underset{U\∈R^{p\×k},\left|\right|U\left|\right|-I}{\min }J\left(\tilde{U}\right)\≈\mathrm{\β}\left(\underset{i=1}{\overset{n}{\Σ}}\right||{\stackrel{\‾}{x}}_i-\tilde{U}{y}_i|{|}^2)+\left|\right|{\stackrel{\‾}{x}}_t-\tilde{U}{y}_t|{|}^2---\left(8\right)$

其中，q为过采样数量与原始数据大小的比例，n为实例数量，所以β是过采样后，对历史数据信息的权重因子，因而，通过求解式(8)，就能得到原始数据过采样后的数据子空间矩阵

(八)投影近似基追踪

采用投影近似信号子空间估计方法，近似获取过采样下数据的子空间并跟踪和匹配t时刻下的投影近似基；

如果满足以下条件，式(8)中在线过采样PCA的最小化数据重构误差可以最小化，

$U\left(t\right)=C_{xy}\left(t\right)C_{yy}^{-1}\left(t\right)---\left(9\right)$

其中，

$C_{yy}\left(t\right)=\underset{i=1}{\overset{t}{\Σ}}{\mathrm{\β}}^{t-i}y\left(i\right)y^T\left(i\right)={\mathrm{\βC}}_{yy}(t-1)+y\left(t\right)y^T\left(t\right);$

当每次计算只考虑数据的一个投影近似基时，且设参考的投影近似子空间和对应的最大特征值、显著投影近似基向量分别为U_s、d_k(s)和u_k(s)，简化式(9)为

$u_k\left(t\right)=u_k\left(s\right)+\frac{1}{d_k\left(t\right)}\[x_k\left(t\right)-u_k\left(s\right)y_k\left(t\right)\]y_k{\left(t\right)}^T---\left(10\right)$

其中，d_k(t)＝βd_k(s)+y_k(t)²，k∈[1,l]，l为子空间维数，x_k(t)为去均值后的原始数据输入向量；

求得第一个投影近似基u₁(t)后，可由式(11)更新x_k+1(t)，用于求解下一个投影近似基；

x_k+1(t)＝x_k(t)-u_k(t)y_k(t) (11)

最终得到的t时刻所有u_k(t)张成的投影近似空间U(t)，其中最大特征值对应的投影近似基为u₁(t)；式(10)所得到U(t)中的投影近似基并不正交，最后利用Gram-Schmidt算法，正交化所有投影近似基；

利用估计的最显著投影近似基方向u₁(t)和参考的最显著投影近似基方向u₁(s)，分别向参考投影近似空间U_s和估计的投影近似空间U(t)投影，最后，用所得到的夹角平均，判断t时刻输入向量的异常程度；

$s_t=1-\left(\right|\left|U_1^T\right(t\left)u_1\right(s\left)\right||+|\left|U_S^T{u}_1\right(t\left)\right|\left|\right)/2---\left(12\right)$

s_t代表了t时刻输入向量对数据的投影近似子空间方向影响，值越大则输入向量的异常程度越大；

令去均值后的经归一化处理后飞行数据流子集为式(10)中x_k(t)，u_k(t)即为t时刻到来后所估计的第k个投影近似基，式(11)更新x_k+1(t)，用于求解下一个投影近似基，其中，参考数据的投影近似空间可在检测开始时，由预定义的初始值，利用滑窗获取的观察矩阵子集通过式(10)在线训练得到，最终式(12)所得到的异常分数s_t代表t时刻飞行数据流子集的异常程度，在线异常检测中异常阈值可以设置为过去m个异常分数的3δ统计。

本发明的技术效果：

针对无人机飞行数据异常检测的准确率问题，本发明提出了一种基于过采样投影近似基追踪(OverSampling Projection Approximation Basis Pursuit，OSPABP)的无监督在线异常检测方法，并据此提出无人机飞行数据在线异常检测框架。首先，消除飞行数据流的量纲，抽取相关的飞行数据观测矩阵子集，然后使用过采样PCA放大子集输入向量异常效果，以提高异常和正常数据间可分性，并实现原始数据向子空间转换。在此基础上，引入投影近似思想简化子空间基向量求解，并通过对数据投影近似子空间基向量的追踪和匹配检测飞行数据流异常，同时利用其对飞行模式切换的不敏感性抑制飞行模式造成的误检测，从而改善和提高了检测准确率。

附图说明

图1为无人机飞行数据在线异常检测框图。

图2为在线异常检测的异常分数示意图。

图3为公开数据在线异常检测ROC曲线对比图。

具体实施方式

下面将结合附图对本发明做进一步的详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式，但本发明的保护范围不限于下述实施例。

如图1所示，本实施例所涉及的一种基于过采样投影近似基追踪无人机飞行数据异常检测方法，步骤如下：

(一)定义

定义1：飞行数据流模型A是以一定采样频率，连续不断，随时间增长且只能读取一次的p维时间序列。

$A=\{{\stackrel{\→}{a}}_1,{\stackrel{\→}{a}}_2,\mathrm{...},{\stackrel{\→}{a}}_t,{\stackrel{\→}{a}}_{t+1},\mathrm{...}\}---\left(1\right)$

其中，是飞行数据流模型A在t时刻的输入向量，p为所监控参数总数，也是向量维度，a_t,j∈R表示参数j在时刻t的值，j∈[1,p]分别对应所监控的不同飞行参数，如气压高度、气压升降速度、垂向速度、俯仰角、组合高度、GPS垂向速度、GPS高度等。

定义2：观察矩阵H_p×m。无限的数据流无法全部存储和处理，可通过保持固定窗口长度的滑窗存储部分数据，并作归一化以统一数据尺度，从而建立的观察矩阵H_p×m：

$H_{p\×m}=\left\{{\stackrel{\→}{a}}_{t-m}\mathrm{......}{\stackrel{\→}{a}}_{t-1}\right\}---\left(2\right)$

其中，H_p×m维度为p×m，列数m为滑窗所记录的过去时间点长度，行数p为所记录的飞行参数数量。

定义3：观察矩阵子集H_r×m和飞行数据流子集H_r×m为从观察矩阵H_p×m中抽取的相关的观察矩阵子集。数据子集如高度类数据子集，垂向速度类数据子集等。H_r×m中行数r为子集维度且r＜p，列数m为所记录的过去时间点长度。是观察矩阵子集对应的输入向量子集，r为子集维度且r＜p。

定义4：飞行数据流异常检测。根据观察矩阵子集H_r×m，判断t时刻输入向量子集是否为异常。

$Anomaly\left(t\right)=(H_{r\×m},{\stackrel{\→}{a}}_t(r\left)\right)---\left(3\right)$

(二)异常检测总体框架

无人机飞行数据在线异常检测框架，如图1虚线方框所示。框架包括滑窗与归一化、观测矩阵子集抽取和过采样投影近似基追踪(OSPABP)三个部分，输入为多维的飞行数据流输出为时刻t的异常报警。

(1)滑窗与归一化是为了建立定义2所示的观测矩阵，其中，滑窗从无限的飞行数据流中存储并更新预设宽度的飞行数据；归一化主要考虑到原始飞行数据中，不同量纲的参数相互间不存在可比性，将影响后续异常检测结果，归一化到同一数量级后，可解决原始参数间量纲不统一问题。

(2)飞行数据中相关的参数组成的数据子集是本实施例提出方法的检测对象，因此，需要从观测矩阵中抽取相关子集(如高度类参数子集、垂向速度类参数子集等)。

(3)过采样投影近似基追踪(OSPABP)在线异常检测算法。当原始数据空间出现异常数据时，数据子空间方向会随之发生角度变化。因而，OSPABP实质上通过度量子空间方向变化程度实现异常检测，同时利用子空间方向对飞行模式切换的不敏感性，抑制飞行模式切换造成的误检测，改善检测准确率。然而，数据流应用中，异常和正常数据子空间方向较小的差异，导致两者可分性降低。并且传统求解代表子空间方向的基向量的方法复杂度高，并不适合机载在线应用。因此，采用OSPABP方法放大输入数据异常效果，然后近似估计、跟踪和匹配代表子空间方向的数据投影近似基，从而简化子空间基向量求解，并实现异常检测。

(三)滑窗与归一化

为消除不同参数的量纲并考虑数据的统计分布，将输入飞行数据输入向量利用式(4)进行Z-score变换，从而将滑窗中的数据归一化为均值为0、方差为1的数据。输入向量的Z-score变化为

$Z\left({\stackrel{\→}{a}}_t\right)=\{Z({\stackrel{\→}{a}}_{t,1},H_1^T),\mathrm{...},Z({\stackrel{\→}{a}}_{t,p},H_p^T)\}---\left(4\right)$

其中为向量均值，σ_x为向量方差。a_t,j为飞行数据向量中的参数，为观测矩阵H_p×m中参数j所对应的过去m点数据。

(四)观测矩阵子集抽取

从观察矩阵中抽取相关的子集H_r×m，作为异常检测算法的输入。观察矩阵子集中相关飞行参数的选取，可对历史的飞行数据进行相关性分析确定。选取相互间具有较大相关系数的参数构成子集，如相关系数大于0.997。

(五)过采样投影近似基追踪

将原始的多维飞行数据观测矩阵，处理为多个相关的子集矩阵H_r×m后。可对观测矩阵中的子集矩阵采用本实施例提出的过采样投影近似基追踪(OSPABP)方法跟踪和匹配t时刻输入数据带来的数据子空间方向变化，从而检测时刻t的飞行数据是否存在异常。

OSPABP方法包括两个步骤。第一步，为了提高数据流中异常和正常数据间的可分性以及获取数据子空间，使用过采样PCA放大子集输入向量异常效果，并实现原始数据向子空间转换。第二步，为满足飞行数据流处理计算轻量级需求，采用投影近似子空间估计方法，近似求解子空间方向，并追踪和匹配子空间方向变化，从而实现飞行数据流异常检测。下面以实例方式对实际无人机飞行数据进行分析。

(六)过采样投影近似基追踪(OSPABP)方法

本节详细介绍过采样投影近似基追踪(OSPABP)方法，分过采样PCA和投影近似基追踪两部分进行论述。

为获取数据的子空间，需采用子空间映射方法。主成分分析(PrincipalComponent Analysis,PCA)是将原始数据通过线性变换，获取数据子空间的方法。然而在不断更新的数据流中，单个异常数据对数据子空间变化程度影响轻微，导致PCA获取的异常数据子空间和正常数据子空间差异较小，不易区分。因此本实施例使用过采样PCA，通过复制目标数据，然后向子空间映射的同时，放大了异常数据对数据子空间影响，提高异常和正常数据间可分性，从而能更有效的评估异常数据带来数据子空间变化。

为了获得子空间方向，需求解子空间基向量，而传统求解方法，如特征值分解和奇异值分解方法，需对协方差矩阵进行分解，过高的计算复杂度(O(N³)，N为数据长度)影响算法的实时处理能力，因而并不适合机载在线应用。因此，本实施例采用投影近似子空间估计方法，以递推投影近似形式估计子空间，使计算复杂度下降两个数量级，简化子空间基向量求解，从而满足机载处理计算实时性要求。

(七)过采样PCA

假设过采样t时刻的目标实例次，则相应过采样PCA方程为：

${\mathrm{\Σ}}_{\tilde{A}}{\tilde{u}}_t=\mathrm{\λ}{\tilde{u}}_t---\left(5\right)$

其中，过采样后数据矩阵为的均值为μ，协方差矩阵为为特征向量。为原始数据矩阵，每行x_i代表i时刻的p维空间数据实例，n为实例数量，x_t为t时刻新的目标实例。协方差矩阵随着每一个新的目标实例x_t而变化，如果采用传统的求解特征值和特征向量方法，每次新的目标数据到来，都需要重复计算，导致较高的计算开销。

标准PCA可看成判断投影数据有最大方差的子空间，可通过最小化数据重构误差J(U)求解。

$\underset{U\∈R^{p\×k},\left|\right|U\left|\right|=I}{\min }J\left(U\right)=\underset{i=1}{\overset{n}{\Σ}}\left|\right|{\stackrel{\‾}{x}}_i-{\mathrm{UU}}^T{\stackrel{\‾}{x}}_i|{|}^2---\left(6\right)$

其中U为包含k个显著特征向量的子空间矩阵，为去均值的第i个原始数据，为去均值的第i个原始数据向特征向量子空间投影，是从投影子空间重构估计得到的第i个用估计为前一个时刻估计的特征向量空间。从而，式(6)转化为：

$\underset{U\∈R^{p\×k},\left|\right|U\left|\right|=I}{\min }J\left(U_t\right)=\underset{i=1}{\overset{t}{\Σ}}\left|\right|{\stackrel{\‾}{x}}_i-U_i{y}_i|{|}^2---\left(7\right)$

同样，如果过采样目标实例次，在线过采样PCA的最小化数据重构误差为：

$\underset{U\∈R^{p\×k},\left|\right|U\left|\right|-I}{\min }J\left(\tilde{U}\right)\≈\mathrm{\β}\left(\underset{i=1}{\overset{n}{\Σ}}\right||{\stackrel{\‾}{x}}_i-\tilde{U}{y}_i|{|}^2)+\left|\right|{\stackrel{\‾}{x}}_t-\tilde{U}{y}_t|{|}^2---\left(8\right)$

其中，q为过采样数量与原始数据大小的比例，n为实例数量。所以β是过采样后，对历史数据信息的权重因子。因而，通过求解式(8)，就能得到原始数据过采样后的数据子空间矩阵

(八)投影近似基追踪

为了解决传统子空间基向量的求解方法无法满足无人机飞行数据流在线异常检测对计算轻量级和实时性要求，本实施例引入投影近似信号子空间估计方法以简化式(8)的求解。

子空间估计是一种估计信号子空间特征值和特征向量的方法，在信噪比估计、图像压缩等现代信号处理应用中扮演重要角色。投影近似子空间估计方法，以递推投影近似形式估计子空间，使计算复杂度下降两个数量级，因而可应用于实时性要求较高的信号处理领域。本实施例采用投影近似信号子空间估计方法，近似获取过采样下数据的子空间并跟踪和匹配t时刻下的投影近似基。

如果满足以下条件，式(8)中在线过采样PCA的最小化数据重构误差可以最小化，

$U\left(t\right)=C_{xy}\left(t\right)C_{yy}^{-1}\left(t\right)---\left(9\right)$

其中，

$C_{yy}\left(t\right)=\underset{i=1}{\overset{t}{\Σ}}{\mathrm{\β}}^{t-i}y\left(i\right)y^T\left(i\right)={\mathrm{\βC}}_{yy}(t-1)+y\left(t\right)y^T\left(t\right).$

当每次计算只考虑数据的一个投影近似基时，且设参考的投影近似子空间和对应的最大特征值、显著投影近似基向量分别为U_s、d_k(s)和u_k(s)，简化式(9)为

$u_k\left(t\right)=u_k\left(s\right)+\frac{1}{d_k\left(t\right)}\[x_k\left(t\right)-u_k\left(s\right)y_k\left(t\right)\]y_k{\left(t\right)}^T---\left(10\right)$

其中，d_k(t)＝βd_k(s)+y_k(t)²，k∈[1,l]，l为子空间维数，x_k(t)为去均值后的原始数据输入向量。

求得第一个投影近似基u₁(t)后，可由式(11)更新x_k+1(t)，用于求解下一个投影近似基。

x_k+1(t)＝x_k(t)-u_k(t)y_k(t) (11)

最终得到的t时刻所有u_k(t)张成的投影近似空间U(t)。其中最大特征值对应的投影近似基为u₁(t)。式(10)所得到U(t)中的投影近似基并不正交，最后利用Gram-Schmidt算法，正交化所有投影近似基。

利用估计的最显著投影近似基方向u₁(t)和参考的最显著投影近似基方向u₁(s)，分别向参考投影近似空间U_s和估计的投影近似空间U(t)投影，如式(12)所示。最后，用所得到的夹角平均，判断t时刻输入向量的异常程度。

$s_t=1-\left(\right|\left|U_1^T\right(t\left)u_1\right(s\left)\right||+|\left|U_S^T{u}_1\right(t\left)\right|\left|\right)/2---\left(12\right)$

s_t代表了t时刻输入向量对数据的投影近似子空间方向影响，值越大则输入向量的异常程度越大。

令去均值后的经归一化处理后飞行数据流子集为式(10)中x_k(t)，u_k(t)即为t时刻到来后所估计的第k个投影近似基。式(11)更新x_k+1(t)，用于求解下一个投影近似基。其中，参考数据的投影近似空间可在检测开始时，由预定义的初始值，利用滑窗获取的观察矩阵子集通过式(10)在线训练得到。最终式(12)所得到的异常分数s_t代表t时刻飞行数据流子集的异常程度。在线异常检测中异常阈值可以设置为过去m个异常分数的3δ统计。

综上，本实施例所提出的过采样投影近似基追踪(OSPABP)方法旨在提高异常和正常数据子空间的可分性，并且采用一种计算轻量级的方式求解和匹配子空间方向，从而实现数据流的异常检测，并满足飞行数据处理的计算实时性要求。

实验验证

为验证本实施例提出方法的可行性和对真实数据的适用性，使用明尼苏达大学无人机实验室真实飞行数据进行实验验证。实验计算平台为Dell Inspiron7447，CPU为IntelCore [email protected]，软件采用Matlab2015。

将提出的过采样投影近似基追踪(OSPABP)方法用于飞行数据在线异常检测，并与基于在线马氏距离(Online Mahalanobis Distance,Online MD)、贝叶斯网络(Bayesiannetwork，BN)、典型相关性分析(Canonical Correlation Analysis，CCA)和KOAD(Kernel-based Online Anomaly Detection)的四种在线异常检测方法进行对比，验证提出方法的可行性和适用性。

图2给出针对明尼苏达大学公开数据集中高度类子集(GPS高度alt、气压高度h和导航高度navalt)在线异常检测的异常分数示意图。横坐标是采样点，纵坐标是高度值/异常分数。其中，图2a为经过滑窗和归一化处理后的高度类相关数据子集，可看出无人机起飞后在空中起伏飞行。数据共计13212个时刻采样点。其中在1145-1735时刻气压高度h和GPS高度alt分别出现明显漂移异常，在2056-2288时刻飞行参数GPS高度alt出现漂移异常，从而与气压高度h和导航高度navalt趋势不一致。图2b和图2c分别为对比方法中检测效果最好的Online MD方法和本实施例方法(OSPABP)针对数据流中每一时刻数据所在线计算的异常分数值。由图2可知，本实施例提出的方法准确检测出1145-1735时刻(图中竖点画线标记区域1)和2056-2288(图中竖点画线标记区域2)时刻存在的异常。虽然，在4494时刻、6248时刻、7651时刻、10950时刻、12430时刻左右，随着飞行模式的切换，出现少量误检测。但是本实施例的方法相比对比方法中检测效果最好的Online MD方法，异常分数受飞行模式切换的影响明显小，误检的正常点更少。

图3给出本实施例方法(OSPABP)、Online MD、BN、CCA和KOAD对公开数据集中高度类子集在线异常检测的ROC曲线图。从图3可以看出，相比其他四种方法，本实施例方法的ROC曲线(左侧星型线)更为陡峭，说明本实施例的方法对正常数据保持较小的误检率FPR下，对异常数据的检测率TPR迅速趋近理论值1，即本实施例的方法对飞行数据的异常检测效果最优。

表1给出本实施例的方法和对比方法的在线异常检测AUC分数统计和本实施例方法AUC分数改善效果。本实施例的方法AUC分数达到0.9836，相比基于Online MD、BN、CCA和KOAD的方法分别提高了5.28％、20.13％、17.90％和32.65％。表明本实施例的方法对真实无人机数据中的异常数据和正常数据的检测准确度最高，接近理论值1，检测效果最好。

表1公开数据在线异常检测AUC分数对比

表2给出本实施例的方法和对比方法的在线异常检测误检率FPR统计。本实施例的方法在检测率TPR分别为90％、95％、97％时，误检率FPR仅为2.49％、3.25％和10.82％。相比其他四种方法中表现最好的基于Online MD的方法，误检率FPR分别相对降低了86.63％、84.38％和53.82％。

表2公开数据在线异常检测误检率FPR对比

为进一步验证本实施例的方法是否满足无人机飞行数据流在线异常检测对计算和存储要求，本实施例对涉及方法的计算复杂和存储复杂度进行分析，并测量方法针对高度类子集数据每个时刻所处理需要的时间。表3为本实施例的方法(OSPABP)和对比方法对每一时刻多维输入向量的计算和存储复杂度，以及处理时间统计。其中m为观察矩阵宽度，r为输入向量子集维度，l为数据子空间维度，l<＝r<<m。

表3算法复杂度和输入向量处理时间统计

由表3可知，本实施例方法的存储复杂度为O(rl)，与m无关，只与r和l有关，而r和l的值远小于m的值，实验中，m为500，r为3，l为1。所以，本实施例方法的存储复杂度最小，满足无人机飞行数据实时处理存储要求。本实施例的方法计算复杂度为O(rl+rl²)，也与m无关，只大于基于BN的方法，复杂度相对较小。对于每采样时刻输入向量的处理时间上，本实施例的方法与对比方法同为μs级，满足无人机飞行数据实时处理时间要求。

综上所述，通过对公开无人机真实飞行数据的处理，验证了本实施例的方法对真实飞行数据异常检测的适用性。实验表明本实施例所提出的方法，在消耗同等级别计算时间下，相对基于Online MD、BN、CCA、KOAD等方法误检测率FPR得到明显降低，准确率AUC分数得到有效提高。

以上所述，仅为本发明较佳的具体实施方式，这些具体实施方式都是基于本发明整体构思下的不同实现方式，而且本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (1)

1.一种基于过采样投影近似基追踪的无人机飞行数据异常检测方法，其特征在于，

(一)定义

定义1：飞行数据流模型A是以一定采样频率，连续不断，随时间增长且只能读取一次的p维时间序列；

$A=\{{\stackrel{\&RightArrow;}{a}}_1,{\stackrel{\&RightArrow;}{a}}_2,...,{\stackrel{\&RightArrow;}{a}}_t,{\stackrel{\&RightArrow;}{a}}_{t+1},...\}---\left(1\right)$

其中， 是飞行数据流模型A在t时刻的输入向量，p为所监控参数总数，也是向量维度，a_t,j∈R表示参数j在时刻t的值，j∈[1,p]分别对应所监控的不同飞行参数，如气压高度、气压升降速度、垂向速度、俯仰角、组合高度、GPS垂向速度、GPS高度；

定义2：观察矩阵H_p×m，无限的数据流无法全部存储和处理，可通过保持固定窗口长度的滑窗存储部分数据，并作归一化以统一数据尺度，从而建立的观察矩阵H_p×m：

$H_{p\&times;m}=\{{\stackrel{\&RightArrow;}{a}}_{t-m}......{\stackrel{\&RightArrow;}{a}}_{t-1}\}---\left(2\right)$

其中，H_p×m维度为p×m，列数m为滑窗所记录的过去时间点长度，行数p为所记录的飞行参数数量；

定义3：观察矩阵子集H_r×m和飞行数据流子集H_r×m为从观察矩阵H_p×m中抽取的相关的观察矩阵子集，数据子集如高度类数据子集，垂向速度类数据子集，H_r×m中行数r为子集维度且r＜p，列数m为所记录的过去时间点长度，是观察矩阵子集对应的输入向量子集，r为子集维度且r＜p；

定义4：飞行数据流异常检测，根据观察矩阵子集H_r×m，判断t时刻输入向量子集是否为异常；

$Anomaly\left(t\right)=(H_{r\&times;m},{\stackrel{\&RightArrow;}{a}}_t(r\left)\right)---\left(3\right)$

(二)异常检测总体框架

无人机飞行数据在线异常检测框架，框架包括滑窗与归一化、观测矩阵子集抽取和过采样投影近似基追踪三个部分，输入为多维的飞行数据流输出为时刻t的异常报警；

(1)滑窗与归一化是为了建立定义2所示的观测矩阵，其中，滑窗从无限的飞行数据流中存储并更新预设宽度的飞行数据；归一化主要考虑到原始飞行数据中，不同量纲的参数相互间不存在可比性；

(2)飞行数据中相关的参数组成的数据子集是本方法的检测对象，需要从观测矩阵中抽取相关子集；

(3)过采样投影近似基追踪在线异常检测算法，当原始数据空间出现异常数据时，数据子空间方向会随之发生角度变化，OSPABP实质上通过度量子空间方向变化程度实现异常检测；

(三)滑窗与归一化

为消除不同参数的量纲并考虑数据的统计分布，将输入飞行数据输入向量利用式(4)进行Z-score变换，从而将滑窗中的数据归一化为均值为0、方差为1的数据，输入向量的Z-score变化为

$Z\left({\stackrel{\&RightArrow;}{a}}_t\right)=\{Z({\stackrel{\&RightArrow;}{a}}_{t,1},H_1^T),...,Z({\stackrel{\&RightArrow;}{a}}_{t,p},H_p^T)\}---\left(4\right)$

其中 为向量均值，σ_x为向量方差，a_t,j为飞行数据向量中的参数，为观测矩阵H_p×m中参数j所对应的过去m点数据；

(四)观测矩阵子集抽取

从观察矩阵中抽取相关的子集H_r×m，作为异常检测算法的输入，观察矩阵子集中相关飞行参数的选取，可对历史的飞行数据进行相关性分析确定，选取相互间具有较大相关系数的参数构成子集，如相关系数大于0.997；

(五)过采样投影近似基追踪

将原始的多维飞行数据观测矩阵，处理为多个相关的子集矩阵H_r×m后，可对观测矩阵中的子集矩阵采用过采样投影近似基追踪方法跟踪和匹配t时刻输入数据带来的数据子空间方向变化，从而检测时刻t的飞行数据是否存在异常；

OSPABP方法包括两个步骤：第一步，使用过采样PCA放大子集输入向量异常效果，并实现原始数据向子空间转换；第二步，采用投影近似子空间估计方法，近似求解子空间方向，并追踪和匹配子空间方向变化，从而实现飞行数据流异常检测；

(六)过采样投影近似基追踪方法

采用投影近似子空间估计方法，以递推投影近似形式估计子空间，使计算复杂度下降两个数量级，简化子空间基向量求解，从而满足机载处理计算实时性要求；

(七)过采样PCA

假设过采样t时刻的目标实例次，则相应过采样PCA方程为：

${\mathrm{\&Sigma;}}_{\tilde{A}}{\tilde{u}}_t=\mathrm{\&lambda;}{\tilde{u}}_t---\left(5\right)$

其中，过采样后数据矩阵为 的均值为μ，协方差矩阵为 为特征向量；为原始数据矩阵，每行x_i代表i时刻的p维空间数据实例，n为实例数量，x_t为t时刻新的目标实例；协方差矩阵随着每一个新的目标实例x_t而变化；

标准PCA可看成判断投影数据有最大方差的子空间，可通过最小化数据重构误差J(U)求解；

$\underset{U\&Element;R^{p\&times;k},\left|\right|U\left|\right|=I}{min}J\left(U\right)=\underset{i=1}{\overset{n}{\&Sigma;}}\left|\right|\stackrel{\&OverBar;}{x_i}-{\mathrm{UU}}^T\stackrel{\&OverBar;}{x_i}|{|}^2---\left(6\right)$

其中U为包含k个显著特征向量的子空间矩阵，为去均值的第i个原始数据，为去均值的第i个原始数据向特征向量子空间投影，是从投影子空间重构估计得到的第i个用估计 为前一个时刻估计的特征向量空间，从而，式(6)转化为：

$\underset{U\&Element;R^{p\&times;k},\left|\right|U\left|\right|=I}{min}J\left(U_t\right)=\underset{i=1}{\overset{t}{\&Sigma;}}\left|\right|\stackrel{\&OverBar;}{x_i}-U_i{y}_i|{|}^2---\left(7\right)$

同样，如果过采样目标实例次，在线过采样PCA的最小化数据重构误差为：

$\underset{U\&Element;R^{p\&times;k},\left|\right|U\left|\right|=I}{min}J\left(\tilde{U}\right)\&ap;\mathrm{\&beta;}\left(\underset{i=1}{\overset{n}{\&Sigma;}}\right||\stackrel{\&OverBar;}{x_i}-\tilde{U}{y}_i|{|}^2)+\left|\right|\stackrel{\&OverBar;}{x_t}-\tilde{U}{y}_t|{|}^2---\left(8\right)$

其中，q为过采样数量与原始数据大小的比例，n为实例数量，所以β是过采样后，对历史数据信息的权重因子，因而，通过求解式(8)，就能得到原始数据过采样后的数据子空间矩阵

(八)投影近似基追踪

采用投影近似信号子空间估计方法，近似获取过采样下数据的子空间并跟踪和匹配t时刻下的投影近似基；

如果满足以下条件，式(8)中在线过采样PCA的最小化数据重构误差可以最小化，

$U\left(t\right)=C_{xy}\left(t\right)C_{yy}^{-1}\left(t\right)---\left(9\right)$

其中，

$C_{yy}\left(t\right)=\underset{i=I}{\overset{t}{\&Sigma;}}{\mathrm{\&beta;}}^{t-i}y\left(i\right)y^T\left(i\right)={\mathrm{\&beta;C}}_{yy}(t-1)+y\left(t\right)y^T\left(t\right);$

当每次计算只考虑数据的一个投影近似基时，且设参考的投影近似子空间和对应的最大特征值、显著投影近似基向量分别为U_s、d_k(s)和u_k(s)，简化式(9)为

$u_k\left(t\right)=u_k\left(s\right)+\frac{1}{d_k\left(t\right)}\&lsqb;x_k\left(t\right)-u_k\left(s\right)y_k\left(t\right)\&rsqb;y_k{\left(t\right)}^T---\left(10\right)$

其中，d_k(t)＝βd_k(s)+y_k(t)²，k∈[1,l]，l为子空间维数，x_k(t)为去均值后的原始数据输入向量；

求得第一个投影近似基u₁(t)后，可由式(11)更新x_k+1(t)，用于求解下一个投影近似基；

x_k+1(t)＝x_k(t)-u_k(t)y_k(t) (11)

最终得到的t时刻所有u_k(t)张成的投影近似空间U(t)，其中最大特征值对应的投影近似基为u₁(t)，式(10)所得到U(t)中的投影近似基并不正交，最后利用Gram-Schmidt算法，正交化所有投影近似基；

利用估计的最显著投影近似基方向u₁(t)和参考的最显著投影近似基方向u₁(s)，分别向参考投影近似空间U_s和估计的投影近似空间U(t)投影，最后，用所得到的夹角平均，判断t时刻输入向量的异常程度；

$s_t=1-\left(\right|\left|U_1^T\right(t)u_1\left(s\right)\left|\right|+\left|\right|U_S^T{u}_1\left(t\right)\left|\right|)/2---\left(12\right)$

s_t代表了t时刻输入向量对数据的投影近似子空间方向影响，值越大则输入向量的异常程度越大；

令去均值后的经归一化处理后飞行数据流子集为式(10)中x_k(t)，u_k(t)即为t时刻到来后所估计的第k个投影近似基，式(11)更新x_k+1(t)，用于求解下一个投影近似基，其中，参考数据的投影近似空间可在检测开始时，由预定义的初始值，利用滑窗获取的观察矩阵子集通过式(10)在线训练得到，最终式(12)所得到的异常分数s_t代表t时刻飞行数据流子集的异常程度，在线异常检测中异常阈值可以设置为过去m个异常分数的3δ统计。