数据预处理之剔除异常值及平滑处理

figure(1); % 第一张图 subplot(3,2,1); plot(t,y); axis([-10 10 -1 6]); grid on title('无噪声信号'); subplot(3,2,2); plot(t,y1,'-*'); axis([-10 10 -1 6]); grid on title('含噪声信号'); subplot(3,2,3); plot(t,y2,'-*'); axis([-10 10 -1 6]); grid on title('3 点单纯移动平均'); subplot(3,2,4); plot(t,y3,'-*'); axis([-10 10 -1 6]); grid on title('9 点单纯移动平均'); subplot(3,2,5); plot(t,y4,'-*'); axis([-10 10 -1 6]); grid on title('3 点线性加权平滑'); subplot(3,2,6); plot(t,y5,'-*'); axis([-10 10 -1 6]); grid on title('9 点线性加权平滑'); figure(2); % 第二张图
min ( yi' k yi k )2
k
通常采用“五点二次平滑” （n=5， k=-2,-1,0,1,2）
2 2 ( yi k A0 A1k A2k ) 0 k 2 2 2 ( yi k A0 A1k A2k )k 0 k 2 2 2 2 ( yi k A0 A1k A2k )k 0 k 2
数据预处理——剔除异常值及平滑处理 测量数据在其采集与传输过程中， 由于环境干扰或人为因素有可 能造成个别数据不切合实际或丢失，这种数据称为异常值。为了恢复 数据的客观真实性以便将来得到更好的分析结果， 有必要先对原始数 据剔除异常值。 另外，无论是人工观测的数据还是由数据采集系统获取的数据， 都不可避免叠加上“噪声”干扰（反映在曲线图形上就是一些“毛刺 和尖峰” ） 。为了提高数据的质量，必须对数据进行平滑处理（去噪声 干扰） 。 （一）剔除异常值。 注：若是有空缺值，或导入 Matlab 数据显示为“NaN” （非数） ， 需要忽略整条空缺值数据，或者填上空缺值。 填空缺值的方法，通常有两种：A. 使用样本平均值填充；B. 使 用判定树或贝叶斯分类等方法推导最可能的值填充（略） 。 一、基本思想： 规定一个置信水平， 确定一个置信限度， 凡是超过该限度的误差， 就认为它是异常值，从而予以剔除。 二、常用方法：拉依达方法、肖维勒方法、一阶差分法。 注意：这些方法都是假设数据依正态分布为前提的。 1. 拉依达方法（非等置信概率） 如果某测量值与平均值之差大于标准偏差的三倍， 则予以剔除。
yi'
n 1 yi 1 2n 1 k n
优点：方法简单，计算方便。 缺点： 方法产生误差会造成信号失真； 前后各 n 个数据无法平滑。 适用性：适用于变化缓慢的数据。 注：n 越大平滑效果越好，但失真也越大。 例2. “9 点单纯移动平均”平滑滤波 代码：
% 建立“n 点单纯移动平均”的滤波函数 % 注意函数要单独保存为与函数名同名的.m 文件 function Y=smooth_data(y,n) m=length(y); j=1; for i=(n-1)/2+1:(m-(n-1)/2) p=i-(n-1)/2;
ˆn xn1 ( xn1 xn2 ) x
比较判别：
ˆn W xn x
注：该方法的特点是 （1）适合于实时数据采集与处理过程； （2）精度除了与允许误差限的大小有关外，还与前两点测量值的 精确度有关； （3）若被测物理量的变化规律不是单调递增或单调递减函数，这
一方法将在函数的拐点处产生较大的误差，严重时将无法使用。 （二）数据的平滑处理 对于一组测量数据（xi，yi） i=1,…,n，不要直接就想着求出的拟 合多项式的线性参数，而是要先平滑处理去掉“噪声” 。平滑处理在 科学研究中广泛使用，它可以减少测量中统计误差带来的影响，尤其 被用于无法利用多次重复测量来得到其平均值的情况和当ｙ i 随ｘi 有徒然变化的那些测量段。 1. “ （2n+1 点）单纯移动平均”平滑滤波 取出以 yi 为中心的前后各 n 个数据（yi-n, …,yi-1,yi,…yi+n）求 平均值代替 yi，即
subplot(3,2,5); plot(t,y8,'-*'); axis([-10 10 -1 6]); grid on title('3 点 rloess 平滑'); subplot(3,2,6); plot(t,y9,'-*'); axis([-10 10 -1 6]); grid on title('9 点 rloess 平滑');
t=-10:0.5:10; n=length(t); y=5./(1+t.^2); % 原始测试数据 % 给测试数据加上噪声干扰
y1=y+0.5*(0.5-rand(1,n)); % 调用多个滤波函数作滤波处理
y2=smooth(y1,3); y3=smooth(y1,9); y4=smooth(y1,3,'lowess'); y5=smooth(y1,9,'lowess'); y6=smooth(y1,3,'loess'); y7=smooth(y1,9,'loess'); y8=smooth(y1,3,'rloess'); y9=smooth(y1,9,'rloess');
xi x n S x
例1. 利用肖维勒方法对下列数据的异常值（2.5000）进行剔除： 1.5034 1.5062 1.5034 1.5024 1.4985 2.5000 1.5007
1.5067 1.4993
1.4969
上述数据保存于文件 erro.dat 代码：
x=load('error.dat'); n=length(x); subplot(2,1,1); plot(x,'o'); title('原始数据') axis([0,n+1,min(x)-1,max(x)+1]); w=1+0.4*log(n); yichang = abs(x-mean(x)) > w*std(x); % 若用拉依达方法，把 w 改成 3 即可，但本组数据将不能成功剔除异常值。 x(yichang)=[]; save errornew.dat x -ASCII subplot(2,1,2); plot(x,'rs'); title('异常值剔除后数据'); axis([0,n+1,min(x)-1,max(x)+1]);
xi x 3S x
2 1 n 1 n 2 其中，x xi 为样本均值，S x 为样本的标准偏 ( x x ) i n i 1 n 1 i 1
1
差。 注：适合大样本数据，建议测量次数≥50 次。 代码实例（略） 。 2. 肖维勒方法（等置信概率） 在 n 次测量结果中，如果某误差可能出现的次数小于半次时， 就予以剔除。 这实质上是规定了置信概率为 1-1/2n，根据这一置信概率，可计 算出肖维勒系数，也可从表中查出，当要求不很严格时，还可按下列 近似公式计算：
五点二次平滑权重系数表： 归一系数 y-2’ y-1’ y0’ y1’ y2’
' Fra Baidu bibliotek0
y-2 31 9 -3 -5 3
y-1 9 13 12 6 -5
y0 -3 12 17 12 -3
y1 -5 6 12 13 9
y2 3 -5 -3 9 31
35 35 35 35 35
1 3 y2 12 y1 17 y0 12 y1 3 y2 35
运行结果： Figure 1
Figure 2
4. 用“smoothts 函数” （盒子法、高斯窗法、指数法）平滑滤波 调用格式： output = smoothts(input) output = smoothts(input, ‘b’, wsize) % 盒子法 % 高斯窗方法
output = smoothts(input, ‘g’, wsize, stdev) output = smoothts(input, ‘e’, n) % 指数法
n 1 0.4ln(n)
Tab1. 肖维勒系数表 n 3 4 5 6 7 8 9 10 11 12
ωn 1.38 1.53 1.65 1.73 1.80 1.86 1.92 1.96 2.00 2.03 n 13 14 15 20 30 40 50 100 200 500
ωn 2.07 2.10 2.13 2.24 2.39 2.49 2.58 2.81 3.02 3.20 如果某测量值与平均值之差的绝对值大于标准偏差与肖维勒系 数之积，则该测量值被剔除。
subplot(3,2,1); plot(t,y); axis([-10 10 -1 6]); grid on title('无噪声信号'); subplot(3,2,2); plot(t,y1,'-*'); axis([-10 10 -1 6]); grid on title('含噪声信号'); subplot(3,2,3); plot(t,y6,'-*'); axis([-10 10 -1 6]); grid on title('3 点二次加权平滑'); subplot(3,2,4); plot(t,y7,'-*'); axis([-10 10 -1 6]); grid on title('9 点二次加权平滑');
q=i+(n-1)/2; Y(j)=sum(y(p:q))/n; j=j+1; end end % 主程序 clc clear t=-15:0.5:15; n=length(t); Y=5./(1+t.^2); % 原始测试数据 % 给测试数据加上噪声干扰 % 调用函数作 9 点滤波处理
y=Y+(0.5-rand(1,n)); y1=smooth_data(y,9);
3. 用“smooth 函数”平滑滤波 调用格式： Z = smooth(Y, span, method) 说明： Z： 平滑后的数据向量 Y： 被平滑的数据向量 span： 平滑点数，缺省为 5 点 method ：平滑方法，缺省为移动平滑，其它还有 ‘moving’ —— Moving average (default)单纯移动平均 ‘lowess’ —— Lowess (linear fit)线性加权平滑 ‘loess’ ——Loess (quadratic fit)二次加权平滑 'sgolay' —— Savitzky-Golay 'rlowess' ——Robust Lowess (linear fit) 'rloess' ——Robust Loess (quadratic fit) 例3. 用 matlab 自带的平滑函数作平滑滤波实例。 代码：
运行结果：x =
1.5034 1.5062 1.5034 1.5067 1.4993 1.4969 y= 1.5034 1.5062 1.4993 1.4969 1.5024 1.4985 2.5000 1.5007
1.5034
1.5024
1.4985
1.5007
1.5067
3. 一阶差分法（预估比较法） 用前两个测量值来预估新的测量值， 然后用预估值与实际测量值 比较，若大于事先给定的允许差限值，则剔除该测量值。 预估值
plot(1:n,Y,1:n,y,'-o',5:n-4,y1,'-*'); legend('无噪声','含噪声','9 点平滑后');
运行结果：
2.
“加权移动平均”平滑滤波
加权的基本思想： 作平均的区间内中心处数据的权值最大，愈远离中心处 的数据权值越小小。这样就减小了对真实信号本身的平滑作 用。 权重系数可以采用最小二乘原理，使平滑后的数据以最 小均方差逼近原始数据。即令