基于B样条曲线的异常数据剔除与修复

摘要

针对大量的复杂的靶场观测数据，通过构造初始拟合数据，利用B样条曲线的方法构造递推模型，使用基于样条平滑方法估计的判断门限对双向检验的结果数据是否异常值进行判定，并且对满足修复条件的数据进行拟合修复，当双向检验的结果不同时，通过构造内推模型来进一步检验。实例分析表明：本文提出的方法相对其他方法能更有效地剔除异常数据，通过数据分段处理能更有效地检验那些可能产生阶段性跳跃的数据，使得模型具有更好的稳定性、更广的适用性和更高的异常数据剔除率，具有较高的实际应用价值，是一套切实可行的异常数据剔除与修复方法

背景

开始时论文写得不太尽如人意，导师确实给了很多的建议和帮助。我的专业是计算数学，所以开始也上了些这方向的课程，例如数值计算方法、线性方程与非线性方程的迭代求解等，由于导师的项目中进行可靠性分析时需要对数据进行预处理，因此我就负责做数据预处理这块，以前都是用matlab编程较多，自从接触了Python之后就感觉很不错，当然入过的坑也比较多。

1 初始拟合数据的判定

设目标参数的离散观测量为$x_k(k=1,2,...,n)$，其对应的时间序列为$t_k(k=1,2,...,n)$，首先构造外推模型，选取合适数量的无异常数据作为起始拟合数据。当采样频率足够高时，相邻观测量之间的差距就会很小，那么可用三阶差分模型对连续四点是否可用作起始拟合点作出判断，差分就是把一个非平稳序列转变成平稳序列，如果序列蕴涵着显著的线性趋势，那么一阶差分就可以实现趋势平稳，如果序列蕴涵着趋势，那么通常低阶（二阶或者三阶）差分就可以提取曲线趋势的影响。差分模型的阶数并不是越高越好，过度的差分会造成有用信息的浪费。如果采样频率并不是足够高时，那么可用二阶差分模型对连续三点是否可用作起始拟合数据作出判断。
假设有连续四点的观测量为$x_k,x_{k+1},x_{k+2},x_{k+3}$，它的三阶差分为：

$$\Delta_k=x_k-3x_{k+1}+3x_{k+2}-x_{k+3} \tag{1}$$ 如果$|\Delta_k<\lambda|$ ，那么就判定连续四点数据中无异常数据，可作为起始拟合数据，其中$\lambda$为判断门限。如果不满足判断条件，那么则往后推进继续选取连续四点进行检验，直到满足条件，即可作为起始拟合数据。判断门限的选取基于以下原则：一般残差$\delta_k \sim N(\mu,\sigma^2)$，$\mu$为残差的均值，$\sigma$为残差的标准差，根据正态分布的原则$3\sigma$ ，判断门限一般取测量精度的3-5倍，这是数据质量较好的情况，如果测量数据质量较差，丢点较多时，可以适当放宽门限，至于取的多少倍是最适合的，可以通过实例对比来验证。由于一般精度未知，可以用样条平滑方法对进行估计，计算公式为: $$\hat\sigma=\sqrt{\frac{\sum(\delta_k-\frac{\sum\delta_k}{n})^2}{n-1}} \tag{2}$$ 式中$\delta_k=y_k-\widetilde y_k$,$n$为目标参数测量数量，$\widetilde y_k$为目标参数平滑值通过三阶差分判断连续四点是否有异常数据，然后以此作为起始拟合数据，再利用滑窗连续检验，对后一个数据进行检验。

2 基于B样条方法的递推模型

B样条曲线是在贝塞尔曲线基础上推广得来的，可以看看样条曲线的详细介绍。
假设给定$n+1$个控制点$P_0,P_1,...,P_n$和一个节点向量$U=\{u_0,u_1,...,u_m\}$,$u$是一个$m+1$个非减元素的集合，$u_0\leq u_1\leq ...\leq u_m$，$u_i$为节点，且$m=n+p+1$（*$p$为B样条曲线的次数）

$p$次B样条曲线的总方程为:

$$S(t)=\sum_{i=0}^n P_iN_{i,p}(t) \tag{3}$$
其中$P_i$是控制曲线的特征点，$N_{i,p}(u)$为$p$次B样条基函数。B样条基函数$N_{i,p}(u)$递归定义为： $$\begin{split} N_{i,0}(u)&= \begin{cases}1,\quad u_i\leq u\leq u_{i+1} \\ 0,\quad other \end{cases} \\ N_{i,p}(u)&=\frac{u-u_i}{u_{i+p}-u_i}N_{i,p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{p+1}}N_{i+1,p-1}(u)\end{split}\tag{4}$$

2.1 正向递推模型

首先是正向递推，拟合点数为4，拟合次数为$P_a$，一般B样条曲线选用3次就足够，本文节点的选择使用均匀分布的。正向递推公式:

$$\hat x_k=S(t_k)=\sum_{i=0}^3P_iN_{i,p_a}(t_k)\quad k=i+1,i+2,......\tag{5}$$
其中节点向量
$U={u_0,u_1,...,u_{3+p_a+1}}$,$t_{i-3}=u_0\leq u_1\leq u_2,...，\leq u_{3+p_a+1}=t_i$
当拟合点数为4时，那么就构造基于3次B样条曲线的递推模型，只要目标参数不是发生异动，那么目标参数运动方程在短时间一般不会超过三阶。首先计算3次递推模型，随后进行逐次降次，直至衰减为线性方程，这样就得到了3个结果分别为$\hat x_k^{(3)},\hat x_k^{(2)},\hat x_k^{(1)}$ 判定条件为：
如果$|\hat x_k^j-x_k|<\lambda(j=1,2,3)$,也就是说开始3次，然后逐次降次，只要有任一次递推结果满足判定条件，那么就可以认为$x_k$ 正常，然后把该点加入起始拟合数据，并同时把第一点拟合数据剔除外推模型。否则认为数据可能异常，此时$x_k$被标记，不能被用于构造外推模型。此时应该继续判定下一点，当再一次有数据被判定为正常时，记录此刻$|\hat x_k^j-x_k|$最小时的拟合次数$P_a(min)$。

2.2 逆向递推模型

然后是逆向递推，拟合点数为4，拟合次数为$P_b$，逆向递推公式为:

$$\hat x_k=S(t_k)=\sum_{i=0}^3 P_iN_{i,P_b}(t_k)\quad k=i^\prime-1,i^\prime-2,......\tag{6}$$
其中节点向量
$U={u_0,u_1,...,u_{3+p_b+1}}$,$t_{i^\prime-3}=u_0\leq u_1\leq u_2,...，\leq u_{3+p_a+1}=t_i^\prime$
逆向递推的方法和正向递推类似，这里不做重复介绍，当逆向递推再次数据被判定为正常时，记录此刻$|\hat x_k^j-x_k|$最小时的拟合次数$P_b(min)$。

3 结果分析

经过双向检验后的数据有如下几种情况:
（1） 如果双向检验都判定为异常数据：
如果 $P_a(min)\neq P_b(min)$
那么此点为异常数据，需剔除。
如果 $P_a(min)=P_b(min)=P$
那么此点也为异常数据，但是满足修复条件，给予修复。但是如果P=2时，则只允许中间存在一个异常数据并且数据无丢失，否则就无法确保拟合模型符合数据规律，不能进行数据修复。数据修复公式为：

$$\hat x_k=S(t_k)=\sum_{i=0}^7P_iN_{i,p}(t_k)\quad k=i+1,i+2,....i^\prime-1\tag{7}$$
其中节点向量
$U={u_0,u_1,...,u_{7+p+1}}$, $t_{i-3}=u_0\leq u_1\leq u_2,...，\leq u_{7+p+1}=t_{i+3}$
（2） 如果正向检验或逆向检验判定为异常数据（也就是一个判定为正常，一个判定为异常）时，说明此点可能发生阶跃等目标运动轨迹异常现象，如果是发生目标运动轨迹阶跃这种情况，那么就会存在连续多点检验异常，这里我们设定最多连续异常数据数量$n_{rmax}$（其值一般依据工程背景进行确定），如果连续异常数据达到了限定值，那么就要对数据进行分段处理，然后分别进行递推。因此对于连续异常少于$n_{rmax}$的数据并不能判定其为异常数据，需进一步检验。这里通过重新构造内推模型，内推模型比外推模型要更加精确。假设单向异常数据点为$x_{k^{\prime\prime}}$，选取该点逆向连续两个正常数据，正向连续两个正常数据构造3次B样条曲线内推模型。公式为： $$\hat x_k=S(t_{k^{\prime\prime}})=\sum_{i=0}^3P_iN_{i,p}(t_{k^{\prime\prime}})\tag{8}$$
其中节点向量
$U={u_0,u_1,...,u_{3+p+1}}$, $t_{k^{\prime\prime}-u}=u_0\leq u_1\leq u_2,...，\leq u_{m}=t_{k^{\prime\prime}+v}$
同样$\lambda$仍为判断门限，$|\hat x_{k^{\prime\prime}}^j-x_{k^{\prime\prime}}|$，只要有任一阶结果满足判断条件，则判定$x_{k^{\prime\prime}}$为正常数据。
（3） 如果双向检验都判定其为正常，则认为该数据正常。
在整个剔除与修复过程中，控制参数的设定也需特别注意，例如允许最多连续异常数据数量$n_{rmax}$，允许最多连续丢点数量$n_{lmax}$的设定，这些都要依据工程背景进行设定，还有判断门限的设定可以根据实际结果进行动态调整。

4 算法思想及实现框架

• 根据所要检测的数据质量使用公式（2）来设定判断门限，后续可以根据实际运算结果来进行修正。还有连续异常数据的判断次数的允许最大值和连续丢失数据的数量的允许最大值的设定，这些可以依据工程背景进行合理的设定。
• 通过正向连续4点检验和逆向连续4点检验进行初始数据判断，从而构造初始拟合递推模型。
• 然后根据公式（5）与公式（6）对数据进行是否异常判断，如果降次到线性模型后还是不满足判定条件则认为该点可能异常，标记此点，往后递推进行下一点检验。
• 如果再一次有数据判断为正常数据后，记录此时拟合外推结果中与目标结果最接近时拟合次数，同时此点需要加入递推模型，相应的剔除递推模型中的第一点数据，然后重新构造递推模型。
• 当连续丢点数量或连续异常数量超过了设定的最大值，那么此时的外推模型失效，需要进行数据分段处理，重新进行外推模型的构造。
• 双向检验均判定为异常的数据则判定为异常，双向检验均判定为正常的数据则判定为正常，一侧判定为正常一侧判定为异常的数据则需要使用公式（8）进一步的检验来判定。而其中对满足修复条件的异常数据可以使用公式（7）给予修复处理。
  程序流程图如下: