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

Abstract

本发明提出了一种基于小波去噪优化的飞行试验数据缺失值填补方法，用于解决现有技术中存在的低信噪比条件下飞行试验数据缺失值填补精度较低的技术问题。实现步骤为，对存在缺失值的飞行试验数据f(n)进行局部线性回归填补；对初填补后的完整飞行试验数据g(n)进行小波去噪；获取飞行试验数据缺失值填补后的飞行试验数据y(n)；本发明采用多重填补的思路，采用稳定性较好的局部线性回归填补方法得到缺失值的粗估计结果后，再利用小波去噪结合分段三次Hermite插值方法得到最终缺失值的精估计结果，解决了现有技术在低信噪比条件下飞行试验数据缺失值填补精度较低的问题。

Description

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

技术领域

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

背景技术

飞行试验数据是指航空装备在试飞过程在产生的实测数据，由试飞过程中航空装备上传感器测得的速度、位置、高度、加速度、压力等内传感器数据，以及雷达等外测量设备的测量数据组成，通常呈现为与时间高度相关的时间序列数据。基于飞行试验数据的后续分析研究是航空装备型号研制和应用技术研究的必需手段和关键环节。然而，在实际的数据采集以及存储过程中，由于多种因素的干扰，难免存在数据缺失的情况，造成飞行试验数据不连续，从而影响后续的数据分析处理。因此，研究飞行试验数据中缺失值的填补方法具有重要的应用意义。

目前飞行试验数据缺失值处理主要有以下三个思路：

第一种是简单删除处理，飞行试验数据中缺失值通常为空值，简单删除处理将缺失值删掉来达到对缺失值进行处理的目的，这种方法当飞行试验数据中缺失值为非完全随机缺失且缺失比例较小时是有效且便利的，但当不满足上述条件时会严重降低飞行试验数据的准确性和连续性，影响后续数据分析处理。

第二种是插值填补方法，其中比较常用的三次样条插值方法、分段三次Hermite插值方法等，这类方法利用飞行试验数据中缺失数据时序值前后的未缺失数据的时序值以及数据值对缺失数据值进行估计，这类方法对平滑数据中的缺失值有着良好的填补效果，但对于实际飞行试验中产生的含噪飞行试验数据中的缺失值填补精度相对较差。

第三种是基于构建填补模型的填补方法，越复杂的模型带来的是越庞大的计算量，其中比较常用的有基于回归模型的填补方法、基于ARMA模型的填补方法、基于LSSVM的填补方法、基于BP神经网络的填补方法等。由于飞行试验数据量庞大且模型机理复杂，构建复杂拟合模型很难保证拟合精度的稳定性，为了兼顾准确性和填补效率，因此针对飞行试验数据比较成熟的是基于局部线性回归模型的数据缺失值填补方法，通过飞行试验数据缺失值附近的未缺失数据构建局部线性回归模型得到缺失值的估计填补值。具体填补步骤如下：首先，构建每个缺失数据的时序值周围最近的L个未缺失数据的时序值为辅助变量向量x，以各时序值的飞行试验数据为目标变量向量y的线性回归填补模型；然后，运用最小二乘法计算线性回归模型的回归系数α和β；最后将飞行试验数据中缺失数据的时序值带入线性回归模型中，通过回归系数进行计算，即可得到飞行试验数据缺失值的填补估计值。这种方法有着较好的填补稳定性，然而在用于飞行试验数据值与时间之间存在着较小的相关关系时，可能导致填补值偏离数据真实值，除此之外该方法在处理低信噪比数据时，由于低信噪比数据局部呈现的严重不平稳现象，导致填补估计值会损失数据的局部趋势信息。

发明内容

本发明的目的在于针对上述现有技术的不足，提出了一种基于小波去噪优化的飞行试验数据缺失值填补方法，用于解决现有技术中在低信噪比条件下对飞行试验数据中缺失值的填补精度较低的技术问题。

为实现上述技术目的，本发明采取的技术方案包括如下步骤：

(1)对存在缺失值的飞行试验数据f(n)进行局部线性回归填补：

(1a)初始化长度为N的待处理飞行试验数据f(n)中包含Q个缺失数据的时序值n_m＝{n_m(1),n_m(2),...,n_m(q),...,n_m(Q)}和K个未缺失数据的时序值n_o＝{n_o(1),n_o(2),...,n_o(k),...,n_o(K)}，并构建每个缺失数据时序值n_m(q)周围最近的L个未缺失数据时序值x＝[n_o(k),n_o(k+1),...,n_o(k+L-1)]为辅助变量向量，以x中各时序值的飞行试验数据值y＝[f(n_o(k)),f(n_o(k+1)),...,f(n_o(k+L-1))]为目标变量向量的线性回归填补模型，其中，N＝Q+K；

(1b)采用最小二乘法通过辅助变量向量x以及目标变量向量y计算线性回归填补模型的回归系数α、β，并通过α、β对缺失数据时序值集合n_m中每个缺失数据时序值n_m(q)的飞行试验数据值进行粗估计，得到n_m的数据估计值集合f₁(n_m)＝{f₁(n_m(1))，f₁(n_m(2))，...，f₁(n_m(q))，...，f₁(n_m(Q))}，最后利用f₁(n_m)对待处理飞行试验数据f(n)中的Q个缺失数据进行填补，得到长度为N的初填补后的完整飞行试验数据g(n)，其中f₁(n_m(q))＝α*n_m(q)+β；

(2)对初填补后的完整飞行试验数据g(n)进行小波去噪：

(2a)对初填补后的完整飞行试验数据g(n)进行分解尺度为J的小波分解，得到第J尺度的小波近似系数序列a_J以及每一尺度的小波细节系数序列d_j，j＝1,2,...,J；

(2b)通过每一尺度的小波细节系数序列d_j计算每一尺度的小波去噪阈值Th_j，并利用Th_j对每一尺度的小波细节系数序列d_j进行阈值判别，将d_j中绝对值小于Th_j的点置零处理，绝对值大于Th_j的点利用去噪函数处理，得到去噪后的每一尺度的小波细节系数序列d_j′，利用d_j′和第J尺度的小波近似系数序列a_J进行小波重构，得到长度为N的去噪后完整飞行试验数据w(n)；

(3)获取飞行试验数据缺失值填补后的飞行试验数据y(n)：

(3a)依次获取缺失数据时序值集合n_m中每一个缺失数据时序值n_m(q)前一个未缺失数据的时序值n_o(k)和后一个未缺失数据的时序值n_o(k+1)；

(3b)根据每一个缺失数据时序值n_m(q)前一个未缺失数据的时序值n_o(k)和其去噪后数据值w(n_o(k))以及后一个未缺失数据的时序值n_o(k+1)和其去噪后数据值w(n_o(k+1))，利用分段三次Hermite插值方法依次对缺失数据时序值集合n_m中每一个缺失数据时序值n_m(q)的飞行试验数据值进行精估计，得到缺失数据时序值集合n_m的数据最终估计值集合w₁(n_m)＝{w₁(n_m(1))，w₁(n_m(2))，...，w₁(n_m(q))，...，w₁(n_m(Q))}；

(3c)利用缺失数据时序值集合n_m的数据最终估计值集合w₁(n_m)对待处理飞行试验数据f(n)中的缺失数据进行填补，得到长度为N的飞行试验数据缺失值填补后的飞行试验数据y(n)。

本发明与现有技术相比，具有如下优点：

1、本发明采用的多重填补思路首先利用局部线性回归填补方法填补结果稳定的特点对缺失值进行粗估计，得到初填补后的完整飞行试验数据；其次利用分段三次Hermite插值方法对于平滑数据填补精度较高的特点，对基于小波去噪平滑后的完整飞行试验数据中缺失值进行精估计，得到最终的缺失值填补结果。与现有技术相比，有效提高了飞行试验数据中缺失值的填补精度。

2、本发明采用小波去噪的方法对初填补后的完整飞行试验数据进行平滑处理，一方面利用小波变换的多分辨分析特点和表征数据局部信息的能力，得到更好的去噪平滑结果；另一方面，通过局部线性回归填补法保证了飞行试验数据的完整性，弥补了小波去噪在应用于含缺失数据时小波重构会出现失真的不足，与现有技术相比，提高了小波去噪精度的同时，更利于进一步对缺失值的精估计。

附图说明

图1为本发明的实现流程图；

图2为用本发明和线性回归填补法以及分段三次Hermite插值方法对长度为20的连续型缺失值的填补结果对比图；

图3为用本发明和线性回归填补法以及分段三次Hermite插值方法在不同信噪比条件下对长度为20的连续型缺失值填补并进行1000次蒙特卡洛试验的平均RMSE对比图。

具体实施方式

下面结合附图和具体实施例，对本发明作进一步详细描述，需要说明的是本发明符合专利法第二条第二项的规定，同时不属于专利法第二十五条所规定的不授予专利权的主体。

参照图1，本发明包括如下步骤：

步骤1，对存在缺失值的飞行试验数据f(n)进行线性回归填补：

步骤1a)初始化长度为N的待处理飞行试验数据f(n)中包含Q个缺失数据的时序值n_m＝{n_m(1),n_m(2),...,n_m(q),...,n_m(Q)}和K个未缺失数据的时序值n_o＝{n_o(1),n_o(2),...,n_o(k),...,n_o(K)}，并构建每个缺失数据时序值n_m(q)周围最近的L个未缺失数据时序值x＝[n_o(k),n_o(k+1),...,n_o(k+L-1)]为辅助变量向量，以x中各时序值的飞行试验数据值y＝[f(n_o(k)),f(n_o(k+1)),...,f(n_o(k+L-1))]为目标变量向量的线性回归填补模型，其中，N＝Q+K；

该步骤中待处理飞行试验数据可以包含航空装备上传感器测得的速度、位置、高度、加速度、压力等内传感器数据，以及雷达等外测量设备的测量数据，本实例中选取的是某传感器的时序测量数据；待处理飞行试验数据f(n)中包含的Q个缺失数据可以是随机型缺失数据，也可以是连续型缺失数据，本实例中选取的是长度为Q的连续缺失数据，缺失数据值为空值；

步骤1b)采用最小二乘法通过辅助变量向量x以及目标变量向量y计算线性回归填补模型的回归系数α、β，并通过α、β对缺失数据时序值集合n_m中每个缺失数据时序值n_m(q)的飞行试验数据值进行粗估计，得到n_m的数据估计值集合f₁(n_m)＝{f₁(n_m(1))，f₁(n_m(2))，...，f₁(n_m(q))，...，f₁(n_m(Q))}，最后利用f₁(n_m)对待处理飞行试验数据f(n)中的Q个缺失数据进行填补，得到长度为N的初填补后的完整飞行试验数据g(n)，其中f₁(n_m(q))＝α*n_m(q)+β；

构建的线性回归填补模型中辅助变量向量为x，目标变量向量为y，线性回归填补模型假设目标变量向量和辅助变量向量满足y＝αx+β的关系，即假设这L个未缺失数据时序值x＝[n_o(k),n_o(k+1),...,n_o(k+L-1)]和以及x中各时序值的飞行试验数据值y＝[f(n_o(k)),f(n_o(k+1)),...,f(n_o(k+L-1))]误差最小的满足该线性回归方程，以此来计算回归系数α、β。这个问题可以用最小二乘法解决，即目标变量向量值与辅助变量的回归模型估计值的误差平方和最小：

令

将该问题转换为求二元函数的极小值，即令

和

即可得到α、β的值，计算公式如下：

在确定线性回归填补模型中的回归系数α、β后，即可将缺失数据时序值集合n_m中每个缺失数据时序值n_m(q)带入线性回归填补模型中，计算得到缺失数据时序值集合n_m中每个缺失数据时序值n_m(q)的飞行试验数据值的粗估计值f₁(n_m(q))＝α*n_m(q)+β；

最后将缺失数据时序值集合n_m的数据估计值集合f₁(n_m)依次替换掉飞行试验数据中的缺失值，完成初填补，得到长度为N的初填补后的完整时序数据g(n)；

步骤2，对初填补后的完整飞行试验数据g(n)进行小波去噪：

步骤2a)对初填补后的完整飞行试验数据g(n)进行分解尺度为J的小波分解，得到第J尺度的小波近似系数序列a_J以及每一尺度的小波细节系数序列d_j，j＝1,2,...,J；

该对初填补后的完整时序数据g(n)进行分解尺度为J的小波分解采用逐层分解的方式，通过由小波变换基函数确定的高通分解滤波器G和低通分解滤波器H，每一尺度的近似系数序列和细节系数序列通过上一尺度的近似系数序列分别与低通分解滤波器和高通分解滤波器卷积得到，最终分解得到第J尺度的小波近似系数序列a_J以及每一尺度的小波细节系数序列d_j，具体计算公式如下：

其中H和G分别为小波变换中的低通分解滤波器和高通分解滤波器，由小波变换基函数确定，

表示卷积操作。

由于小波变换基函数在处理数据时各有特点，当利用小波变换进行数据去噪时，需要选用合适的小波函数作为小波变换基函数，目前常用于去噪的小波变换基函数有dbN系列小波函数、symN系列小波函数、haar小波等，本实例选用但不限于db8小波函数作为小波变换基函数；

小波去噪中小波分解尺度与最终的去噪结果息息相关，通常分解尺度过多，对于各层小波系数序列进行阈值处理会造成信息丢失严重，而分解尺度过少则会降低小波去噪效果，因此小波分解尺度通常为3-8层，本实例中小波分解尺度选用但不限于5层。

步骤2b)通过每一尺度的小波细节系数序列d_j计算每一尺度的小波去噪阈值Th_j，并利用Th_j对每一尺度的小波细节系数序列d_j进行阈值判别，将d_j中绝对值小于Th_j的点置零处理，绝对值大于Th_j的点利用去噪函数处理，得到去噪后的每一尺度的小波细节系数序列d_j′，利用d_j′和第J尺度的小波近似系数序列a_J进行小波重构，得到长度为N的去噪后完整飞行试验数据w(n)；

步骤2b-1)计算每一尺度的小波去噪阈值Th_j，计算公式为：

其中σ_j是第j尺度小波细节系数序列d_j的标准方差值，

Z_j表示第j尺度小波细节系数序列d_j的长度；

步骤2b-2)得到去噪后的每一尺度的小波细节系数序列d_j′，计算公式如下：

该步骤中利用Th_j逐点对每一尺度的小波细节系数序列d_j进行阈值判别，将d_j中绝对值小于Th_j的点置零处理，绝对值大于Th_j的点利用去噪函数处理，本实例选用但不限于软阈值小波去噪函数，得到去噪后的每一尺度的小波细节系数序列d_j′；

步骤2b-3)利用d_j′和第J尺度的小波近似系数序列a_J进行小波重构，计算公式如下：

其中H^*和G^*分别为高通滤波器H和低通滤波器G的对偶算子；

步骤3，获取飞行试验数据缺失值填补后的飞行试验数据y(n)：

(3a)获取缺失数据时序值集合n_m中每一个缺失数据时序值n_m(q)前一个未缺失数据的时序值n_o(k)和后一个未缺失数据的时序值n_o(k+1)；

步骤3b)根据每一个缺失数据时序值n_m(q)前一个未缺失数据的时序值n_o(k)和其去噪后数据值w(n_o(k))以及后一个未缺失数据的时序值n_o(k+1)和其去噪后数据值w(n_o(k+1))，利用分段三次Hermite插值方法依次对缺失数据时序值集合n_m中每一个缺失数据时序值n_m(q)的飞行试验数据值进行精估计，得到缺失数据时序值集合n_m的数据最终估计值集合w₁(n_m)＝{w₁(n_m(1))，w₁(n_m(2))，...，w₁(n_m(q))，...，w₁(n_m(Q))}；

该步骤中具体实现如下：

步骤3b-1)计算时序值n_m(q)的填补系数α_k(n_m(q))、α_k+1(n_m(q))、β_k(n_m(q))、β_k+1(n_m(q))：

步骤3b-2)计算每一个缺失数据时序值n_m(q)的飞行试验数据精估计值w₁(n_m(q))：

步骤3c)利用缺失数据时序值集合n_m的数据最终估计值集合w₁(n_m)对待处理飞行试验数据f(n)中的缺失数据进行填补，得到长度为N的飞行试验数据缺失值填补后的飞行试验数据y(n)。

以下结合仿真对本发明技术效果进行描述

1仿真条件及内容

仿真条件：Matlab R2018a，Window 10

硬件条件：Intel(R)Core(TM)i7-10750H CPU

数据：某飞行试验数据集中的一组时序测量数据，数据长度为7584；

仿真1：对本发明处理连续型缺失值的效果进行描述，采用本发明和线性回归填补法以及分段三次Hermite插值方法，分别对待处理飞行试验数据中长度为20的连续型缺失值进行填补，最终三种方法的填补结果如图2所示。

仿真2：对本发明在不同信噪比条件下缺失值填补性能进行描述，采用本发明和线性回归填补法以及分段三次Hermite插值方法，分别对不同信噪比条件下的待处理飞行试验数据中长度为20的连续型缺失值进行填补，并进行1000次蒙特卡洛试验，得到不同信噪比条件下三种方法的RMSE对比结果如图3所示。

2仿真结果分析

图2中，本发明对于飞行试验数据中缺失值的填补结果更接近于数据理想值，表明本发明有着比现有常用方法更好的缺失值填补效果。

图3中，本发明相较于现有常用缺失值处理方法在不同信噪比条件下，有着更低的填补误差，表明本发明有着比现有常用方法更好的缺失值填补精度。

以上描述仅是本发明的一个具体实例，并未构成对本发明的任何限制，显然对于本领域的专业人员来说，在了解了本发明的内容和原理后，都可能在不背离本发明原理、结构的情况下，进行形式和细节上的各种修改和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。

Claims (6)

1.一种基于小波去噪优化的飞行试验数据缺失值填补方法，其特征在于，包括以下步骤：

(1)对存在缺失值的飞行试验数据f(n)进行局部线性回归填补：

(1a)初始化长度为N的待处理飞行试验数据f(n)中包含Q个缺失数据的时序值n_m＝{n_m(1),n_m(2),...,n_m(q),...,n_m(Q)}和K个未缺失数据的时序值n_o＝{n_o(1),n_o(2),...,n_o(k),...,n_o(K)}，并构建每个缺失数据时序值n_m(q)周围最近的L个未缺失数据时序值x＝[n_o(k),n_o(k+1),...,n_o(k+L-1)]为辅助变量向量，以x中各时序值的飞行试验数据值y＝[f(n_o(k)),f(n_o(k+1)),...,f(n_o(k+L-1))]为目标变量向量的线性回归填补模型，其中，N＝Q+K；

(1b)采用最小二乘法通过辅助变量向量x以及目标变量向量y计算线性回归填补模型的回归系数α、β，并通过α、β对缺失数据时序值集合n_m中每个缺失数据时序值n_m(q)的飞行试验数据值进行粗估计，得到n_m的数据估计值集合f₁(n_m)＝{f₁(n_m(1))，f₁(n_m(2))，...，f₁(n_m(q))，...，f₁(n_m(Q))}，然后利用f₁(n_m)对待处理飞行试验数据f(n)中的Q个缺失数据进行填补，得到长度为N的初填补后的完整飞行试验数据g(n)，其中f₁(n_m(q))＝α*n_m(q)+β；

(2)对初填补后的完整飞行试验数据g(n)进行小波去噪：

(2a)对初填补后的完整飞行试验数据g(n)进行分解尺度为J的小波分解，得到第J尺度的小波近似系数序列a_J以及每一尺度的小波细节系数序列d_j，j＝1,2,...,J；

(2b)通过每一尺度的小波细节系数序列d_j计算每一尺度的小波去噪阈值Th_j，并利用Th_j对每一尺度的小波细节系数序列d_j进行阈值判别，将d_j中绝对值小于Th_j的点置零处理，绝对值大于Th_j的点利用去噪函数处理，得到去噪后的每一尺度的小波细节系数序列d_j′，利用d_j′和第J尺度的小波近似系数序列a_J进行小波重构，得到长度为N的去噪后完整飞行试验数据w(n)；

(3)获取飞行试验数据缺失值填补后的飞行试验数据y(n)：

(3a)依次获取缺失数据时序值集合n_m中每一个缺失数据时序值n_m(q)前一个未缺失数据的时序值n_o(k)和后一个未缺失数据的时序值n_o(k+1)；

(3b)根据每一个缺失数据时序值n_m(q)前一个未缺失数据的时序值n_o(k)和其去噪后数据值w(n_o(k))以及后一个未缺失数据的时序值n_o(k+1)和其去噪后数据值w(n_o(k+1))，利用分段三次Hermite插值方法依次对缺失数据时序值集合n_m中每一个缺失数据时序值n_m(q)的飞行试验数据值进行精估计，得到缺失数据时序值集合n_m的数据最终估计值集合w₁(n_m)＝{w₁(n_m(1))，w₁(n_m(2))，...，w₁(n_m(q))，...，w₁(n_m(Q))}；

(3c)利用缺失数据时序值集合n_m的数据最终估计值集合w₁(n_m)对待处理飞行试验数据f(n)中的缺失数据进行填补，得到长度为N的飞行试验数据缺失值填补后的飞行试验数据y(n)。

2.根据权利要求1所述的基于小波去噪优化的飞行试验数据缺失值填补方法，其特征在于，步骤(1b)中采用最小二乘法通过辅助变量向量x以及目标变量向量y计算线性回归填补模型的回归系数α、β，计算公式为：

其中，

表示辅助变量向量x的均值，

表示目标变量向量的y均值，x_l表示x中的第l个元素，y_l表示y中的第l个元素。

3.根据权利要求1所述的基于小波去噪优化的飞行试验数据缺失值填补方法，其特征在于，步骤(2a)中所述的得到第J尺度的小波近似系数序列a_J以及每一尺度的小波细节系数序列d_j，其计算如下：

对初填补后的完整时序数据g(n)进行分解尺度为J的小波分解采用逐层分解的方式，通过由小波变换基函数确定的高通分解滤波器G和低通分解滤波器H，每一尺度的近似系数序列和细节系数序列通过上一尺度的近似系数序列分别与低通分解滤波器和高通分解滤波器卷积得到，最终分解得到第J尺度的小波近似系数序列a_J以及每一尺度的小波细节系数序列d_j，具体计算公式如下：

其中H和G分别为小波变换中的低通分解滤波器和高通分解滤波器，由小波变换基函数确定，

表示卷积操作。

4.根据权利要求1所述的基于小波去噪优化的飞行试验数据缺失值填补方法，其特征在于，步骤(2b)中所述的每一尺度的小波去噪阈值Th_j，以及去噪后每一尺度的小波细节系数序列d_j′，计算公式分别为：

其中σ_j是第j尺度小波细节系数序列d_j的标准方差值，

Z_j表示第j尺度小波细节系数序列d_j的长度。

5.根据权利要求1所述的基于小波去噪优化的飞行试验数据缺失值填补方法，其特征在于，步骤(2b)中所述的长度为N的去噪后完整飞行试验数据w(n)，计算公式如下：

其中H^*和G^*分别为高通滤波器H和低通滤波器G的对偶算子。

6.根据权利要求1所述的基于小波去噪优化的飞行试验数据缺失值填补方法，其特征在于，步骤(3b)中所述的利用分段三次Hermite插值方法依次对缺失数据时序值集合n_m中每一个缺失数据时序值n_m(q)的飞行试验数据值进行精估计，其实现步骤为：

(3b1)计算时序值n_m(q)的填补系数α_k(n_m(q))、α_k+1(n_m(q))、β_k(n_m(q))、β_k+1(n_m(q))：

(3b2)计算每一个缺失数据时序值n_m(q)的飞行试验数据精估计值w₁(n_m(q))：

Images