基于小波去噪优化的飞行试验数据缺失值填补方法

1.本发明属于数据处理技术领域，涉及一种飞行试验数据的缺失值填补方法，具体涉及一种基于小波去噪优化的飞行试验数据缺失值填补方法。


背景技术：

2.飞行试验数据是指航空装备在试飞过程在产生的实测数据，由试飞过程中航空装备上传感器测得的速度、位置、高度、加速度、压力等内传感器数据，以及雷达等外测量设备的测量数据组成，通常呈现为与时间高度相关的时间序列数据。基于飞行试验数据的后续分析研究是航空装备型号研制和应用技术研究的必需手段和关键环节。然而，在实际的数据采集以及存储过程中，由于多种因素的干扰，难免存在数据缺失的情况，造成飞行试验数据不连续，从而影响后续的数据分析处理。因此，研究飞行试验数据中缺失值的填补方法具有重要的应用意义。
3.目前飞行试验数据缺失值处理主要有以下三个思路：
4.第一种是简单删除处理，飞行试验数据中缺失值通常为空值，简单删除处理将缺失值删掉来达到对缺失值进行处理的目的，这种方法当飞行试验数据中缺失值为非完全随机缺失且缺失比例较小时是有效且便利的，但当不满足上述条件时会严重降低飞行试验数据的准确性和连续性，影响后续数据分析处理。
5.第二种是插值填补方法，其中比较常用的三次样条插值方法、分段三次hermite插值方法等，这类方法利用飞行试验数据中缺失数据时序值前后的未缺失数据的时序值以及数据值对缺失数据值进行估计，这类方法对平滑数据中的缺失值有着良好的填补效果，但对于实际飞行试验中产生的含噪飞行试验数据中的缺失值填补精度相对较差。
6.第三种是基于构建填补模型的填补方法，越复杂的模型带来的是越庞大的计算量，其中比较常用的有基于回归模型的填补方法、基于arma模型的填补方法、基于lssvm的填补方法、基于bp神经网络的填补方法等。由于飞行试验数据量庞大且模型机理复杂，构建复杂拟合模型很难保证拟合精度的稳定性，为了兼顾准确性和填补效率，因此针对飞行试验数据比较成熟的是基于局部线性回归模型的数据缺失值填补方法，通过飞行试验数据缺失值附近的未缺失数据构建局部线性回归模型得到缺失值的估计填补值。具体填补步骤如下：首先，构建每个缺失数据的时序值周围最近的l个未缺失数据的时序值为辅助变量向量x，以各时序值的飞行试验数据为目标变量向量y的线性回归填补模型；然后，运用最小二乘法计算线性回归模型的回归系数α和β；最后将飞行试验数据中缺失数据的时序值带入线性回归模型中，通过回归系数进行计算，即可得到飞行试验数据缺失值的填补估计值。这种方法有着较好的填补稳定性，然而在用于飞行试验数据值与时间之间存在着较小的相关关系时，可能导致填补值偏离数据真实值，除此之外该方法在处理低信噪比数据时，由于低信噪比数据局部呈现的严重不平稳现象，导致填补估计值会损失数据的局部趋势信息。


技术实现要素：

7.本发明的目的在于针对上述现有技术的不足，提出了一种基于小波去噪优化的飞行试验数据缺失值填补方法，用于解决现有技术中在低信噪比条件下对飞行试验数据中缺失值的填补精度较低的技术问题。
8.为实现上述技术目的，本发明采取的技术方案包括如下步骤：
9.(1)对存在缺失值的飞行试验数据f(n)进行局部线性回归填补：
10.(1a)初始化长度为n的待处理飞行试验数据f(n)中包含q个缺失数据的时序值nm＝{nm(1),nm(2),...,nm(q),...,nm(q)}和k个未缺失数据的时序值no＝{no(1),no(2),...,no(k),...,no(k)}，并构建每个缺失数据时序值nm(q)周围最近的l个未缺失数据时序值x＝[no(k),no(k 1),...,no(k l-1)]为辅助变量向量，以x中各时序值的飞行试验数据值y＝[f(no(k)),f(no(k 1)),...,f(no(k l-1))]为目标变量向量的线性回归填补模型，其中，n＝q k；
[0011]
(1b)采用最小二乘法通过辅助变量向量x以及目标变量向量y计算线性回归填补模型的回归系数α、β，并通过α、β对缺失数据时序值集合nm中每个缺失数据时序值nm(q)的飞行试验数据值进行粗估计，得到nm的数据估计值集合f1(nm)＝{f1(nm(1))，f1(nm(2))，...，f1(nm(q))，...，f1(nm(q))}，最后利用f1(nm)对待处理飞行试验数据f(n)中的q个缺失数据进行填补，得到长度为n的初填补后的完整飞行试验数据g(n)，其中f1(nm(q))＝α*nm(q) β；
[0012]
(2)对初填补后的完整飞行试验数据g(n)进行小波去噪：
[0013]
(2a)对初填补后的完整飞行试验数据g(n)进行分解尺度为j的小波分解，得到第j尺度的小波近似系数序列aj以及每一尺度的小波细节系数序列dj，j＝1,2,...,j；
[0014]
(2b)通过每一尺度的小波细节系数序列dj计算每一尺度的小波去噪阈值thj，并利用thj对每一尺度的小波细节系数序列dj进行阈值判别，将dj中绝对值小于thj的点置零处理，绝对值大于thj的点利用去噪函数处理，得到去噪后的每一尺度的小波细节系数序列dj′
，利用dj′
和第j尺度的小波近似系数序列aj进行小波重构，得到长度为n的去噪后完整飞行试验数据w(n)；
[0015]
(3)获取飞行试验数据缺失值填补后的飞行试验数据y(n)：
[0016]
(3a)依次获取缺失数据时序值集合nm中每一个缺失数据时序值nm(q)前一个未缺失数据的时序值no(k)和后一个未缺失数据的时序值no(k 1)；
[0017]
(3b)根据每一个缺失数据时序值nm(q)前一个未缺失数据的时序值no(k)和其去噪后数据值w(no(k))以及后一个未缺失数据的时序值no(k 1)和其去噪后数据值w(no(k 1))，利用分段三次hermite插值方法依次对缺失数据时序值集合nm中每一个缺失数据时序值nm(q)的飞行试验数据值进行精估计，得到缺失数据时序值集合nm的数据最终估计值集合w1(nm)＝{w1(nm(1))，w1(nm(2))，...，w1(nm(q))，...，w1(nm(q))}；
[0018]
(3c)利用缺失数据时序值集合nm的数据最终估计值集合w1(nm)对待处理飞行试验数据f(n)中的缺失数据进行填补，得到长度为n的飞行试验数据缺失值填补后的飞行试验数据y(n)。
[0019]
本发明与现有技术相比，具有如下优点：
[0020]
1、本发明采用的多重填补思路首先利用局部线性回归填补方法填补结果稳定的特点对缺失值进行粗估计，得到初填补后的完整飞行试验数据；其次利用分段三次hermite
插值方法对于平滑数据填补精度较高的特点，对基于小波去噪平滑后的完整飞行试验数据中缺失值进行精估计，得到最终的缺失值填补结果。与现有技术相比，有效提高了飞行试验数据中缺失值的填补精度。
[0021]
2、本发明采用小波去噪的方法对初填补后的完整飞行试验数据进行平滑处理，一方面利用小波变换的多分辨分析特点和表征数据局部信息的能力，得到更好的去噪平滑结果；另一方面，通过局部线性回归填补法保证了飞行试验数据的完整性，弥补了小波去噪在应用于含缺失数据时小波重构会出现失真的不足，与现有技术相比，提高了小波去噪精度的同时，更利于进一步对缺失值的精估计。
附图说明
[0022]
图1为本发明的实现流程图；
[0023]
图2为用本发明和线性回归填补法以及分段三次hermite插值方法对长度为20的连续型缺失值的填补结果对比图；
[0024]
图3为用本发明和线性回归填补法以及分段三次hermite插值方法在不同信噪比条件下对长度为20的连续型缺失值填补并进行1000次蒙特卡洛试验的平均rmse对比图。
具体实施方式
[0025]
下面结合附图和具体实施例，对本发明作进一步详细描述，需要说明的是本发明符合专利法第二条第二项的规定，同时不属于专利法第二十五条所规定的不授予专利权的主体。
[0026]
参照图1，本发明包括如下步骤：
[0027]
步骤1，对存在缺失值的飞行试验数据f(n)进行线性回归填补：
[0028]
步骤1a)初始化长度为n的待处理飞行试验数据f(n)中包含q个缺失数据的时序值nm＝{nm(1),nm(2),...,nm(q),...,nm(q)}和k个未缺失数据的时序值no＝{no(1),no(2),...,no(k),...,no(k)}，并构建每个缺失数据时序值nm(q)周围最近的l个未缺失数据时序值x＝[no(k),no(k 1),...,no(k l-1)]为辅助变量向量，以x中各时序值的飞行试验数据值y＝[f(no(k)),f(no(k 1)),...,f(no(k l-1))]为目标变量向量的线性回归填补模型，其中，n＝q k；
[0029]
该步骤中待处理飞行试验数据可以包含航空装备上传感器测得的速度、位置、高度、加速度、压力等内传感器数据，以及雷达等外测量设备的测量数据，本实例中选取的是某传感器的时序测量数据；待处理飞行试验数据f(n)中包含的q个缺失数据可以是随机型缺失数据，也可以是连续型缺失数据，本实例中选取的是长度为q的连续缺失数据，缺失数据值为空值；
[0030]
步骤1b)采用最小二乘法通过辅助变量向量x以及目标变量向量y计算线性回归填补模型的回归系数α、β，并通过α、β对缺失数据时序值集合nm中每个缺失数据时序值nm(q)的飞行试验数据值进行粗估计，得到nm的数据估计值集合f1(nm)＝{f1(nm(1))，f1(nm(2))，...，f1(nm(q))，...，f1(nm(q))}，最后利用f1(nm)对待处理飞行试验数据f(n)中的q个缺失数据进行填补，得到长度为n的初填补后的完整飞行试验数据g(n)，其中f1(nm(q))＝α*nm(q) β；
[0031]
构建的线性回归填补模型中辅助变量向量为x，目标变量向量为y，线性回归填补模型假设目标变量向量和辅助变量向量满足y＝αx β的关系，即假设这l个未缺失数据时序值x＝[no(k),no(k 1),...,no(k l-1)]和以及x中各时序值的飞行试验数据值y＝[f(no(k)),f(no(k 1)),...,f(no(k l-1))]误差最小的满足该线性回归方程，以此来计算回归系数α、β。这个问题可以用最小二乘法解决，即目标变量向量值与辅助变量的回归模型估计值的误差平方和最小：
[0032][0033]
令将该问题转换为求二元函数的极小值，即令和即可得到α、β的值，计算公式如下：
[0034][0035][0036]
在确定线性回归填补模型中的回归系数α、β后，即可将缺失数据时序值集合nm中每个缺失数据时序值nm(q)带入线性回归填补模型中，计算得到缺失数据时序值集合nm中每个缺失数据时序值nm(q)的飞行试验数据值的粗估计值f1(nm(q))＝α*nm(q) β；
[0037]
最后将缺失数据时序值集合nm的数据估计值集合f1(nm)依次替换掉飞行试验数据中的缺失值，完成初填补，得到长度为n的初填补后的完整时序数据g(n)；
[0038]
步骤2，对初填补后的完整飞行试验数据g(n)进行小波去噪：
[0039]
步骤2a)对初填补后的完整飞行试验数据g(n)进行分解尺度为j的小波分解，得到第j尺度的小波近似系数序列aj以及每一尺度的小波细节系数序列dj，j＝1,2,...,j；
[0040]
该对初填补后的完整时序数据g(n)进行分解尺度为j的小波分解采用逐层分解的方式，通过由小波变换基函数确定的高通分解滤波器g和低通分解滤波器h，每一尺度的近似系数序列和细节系数序列通过上一尺度的近似系数序列分别与低通分解滤波器和高通分解滤波器卷积得到，最终分解得到第j尺度的小波近似系数序列aj以及每一尺度的小波细节系数序列dj，具体计算公式如下：
[0041][0042][0043]
其中h和g分别为小波变换中的低通分解滤波器和高通分解滤波器，由小波变换基函数确定，表示卷积操作。
[0044]
由于小波变换基函数在处理数据时各有特点，当利用小波变换进行数据去噪时，需要选用合适的小波函数作为小波变换基函数，目前常用于去噪的小波变换基函数有dbn
系列小波函数、symn系列小波函数、haar小波等，本实例选用但不限于db8小波函数作为小波变换基函数；
[0045]
小波去噪中小波分解尺度与最终的去噪结果息息相关，通常分解尺度过多，对于各层小波系数序列进行阈值处理会造成信息丢失严重，而分解尺度过少则会降低小波去噪效果，因此小波分解尺度通常为3-8层，本实例中小波分解尺度选用但不限于5层。
[0046]
步骤2b)通过每一尺度的小波细节系数序列dj计算每一尺度的小波去噪阈值thj，并利用thj对每一尺度的小波细节系数序列dj进行阈值判别，将dj中绝对值小于thj的点置零处理，绝对值大于thj的点利用去噪函数处理，得到去噪后的每一尺度的小波细节系数序列dj′
，利用dj′
和第j尺度的小波近似系数序列aj进行小波重构，得到长度为n的去噪后完整飞行试验数据w(n)；
[0047]
步骤2b-1)计算每一尺度的小波去噪阈值thj，计算公式为：
[0048][0049]
其中σj是第j尺度小波细节系数序列dj的标准方差值，zj表示第j尺度小波细节系数序列dj的长度；
[0050]
步骤2b-2)得到去噪后的每一尺度的小波细节系数序列dj′
，计算公式如下：
[0051][0052]
该步骤中利用thj逐点对每一尺度的小波细节系数序列dj进行阈值判别，将dj中绝对值小于thj的点置零处理，绝对值大于thj的点利用去噪函数处理，本实例选用但不限于软阈值小波去噪函数，得到去噪后的每一尺度的小波细节系数序列dj′
；
[0053]
步骤2b-3)利用dj′
和第j尺度的小波近似系数序列aj进行小波重构，计算公式如下：
[0054][0055]
其中h
*
和g
*
分别为高通滤波器h和低通滤波器g的对偶算子；
[0056]
步骤3，获取飞行试验数据缺失值填补后的飞行试验数据y(n)：
[0057]
(3a)获取缺失数据时序值集合nm中每一个缺失数据时序值nm(q)前一个未缺失数据的时序值no(k)和后一个未缺失数据的时序值no(k 1)；
[0058]
步骤3b)根据每一个缺失数据时序值nm(q)前一个未缺失数据的时序值no(k)和其去噪后数据值w(no(k))以及后一个未缺失数据的时序值no(k 1)和其去噪后数据值w(no(k 1))，利用分段三次hermite插值方法依次对缺失数据时序值集合nm中每一个缺失数据时序值nm(q)的飞行试验数据值进行精估计，得到缺失数据时序值集合nm的数据最终估计值集合w1(nm)＝{w1(nm(1))，w1(nm(2))，...，w1(nm(q))，...，w1(nm(q))}；
[0059]
该步骤中具体实现如下：
[0060]
步骤3b-1)计算时序值nm(q)的填补系数αk(nm(q))、α
k 1
(nm(q))、βk(nm(q))、β
k 1
(nm(q))：
[0061][0062]
步骤3b-2)计算每一个缺失数据时序值nm(q)的飞行试验数据精估计值w1(nm(q))：
[0063][0064]
步骤3c)利用缺失数据时序值集合nm的数据最终估计值集合w1(nm)对待处理飞行试验数据f(n)中的缺失数据进行填补，得到长度为n的飞行试验数据缺失值填补后的飞行试验数据y(n)。
[0065]
以下结合仿真对本发明技术效果进行描述
[0066]
1仿真条件及内容
[0067]
仿真条件：matlab r2018a，window 10
[0068]
硬件条件：intel(r)core(tm)i7-10750h cpu
[0069]
数据：某飞行试验数据集中的一组时序测量数据，数据长度为7584；
[0070]
仿真1：对本发明处理连续型缺失值的效果进行描述，采用本发明和线性回归填补法以及分段三次hermite插值方法，分别对待处理飞行试验数据中长度为20的连续型缺失值进行填补，最终三种方法的填补结果如图2所示。
[0071]
仿真2：对本发明在不同信噪比条件下缺失值填补性能进行描述，采用本发明和线性回归填补法以及分段三次hermite插值方法，分别对不同信噪比条件下的待处理飞行试验数据中长度为20的连续型缺失值进行填补，并进行1000次蒙特卡洛试验，得到不同信噪比条件下三种方法的rmse对比结果如图3所示。
[0072]
2仿真结果分析
[0073]
图2中，本发明对于飞行试验数据中缺失值的填补结果更接近于数据理想值，表明本发明有着比现有常用方法更好的缺失值填补效果。
[0074]
图3中，本发明相较于现有常用缺失值处理方法在不同信噪比条件下，有着更低的填补误差，表明本发明有着比现有常用方法更好的缺失值填补精度。
[0075]
以上描述仅是本发明的一个具体实例，并未构成对本发明的任何限制，显然对于本领域的专业人员来说，在了解了本发明的内容和原理后，都可能在不背离本发明原理、结构的情况下，进行形式和细节上的各种修改和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。