数据滤波

数据滤波（精选9篇）

数据滤波 篇1

摘要：本文介绍了点云滤波的基本原理, 对异常点检测问题的特点、分类及应用领域进行了阐述, 同时对异常点检测的各种算法进行了分类研究与深入分析, 最后指出异常点检测今后的研究方向。

关键词：点云,滤波,离群点

1 网格滤波问题

目前网格的光顺算法已经得到广泛研究。网格曲面光顺算法中最经典的算法是基于拉普拉斯算子的方法[1]。通过求取网格曲面的拉普拉斯算子, 并且对网格曲面迭代使用拉普拉斯算法, 可以得到平滑的网格曲面。这种算法的本质是求取网格曲面上某点及其临近点的中心点, 将该中心点作为原顶点的新位置。Jones等根据各顶点的邻域点来预测新顶点位置, 该方法的优点是不需要进行迭代计算。但是上述两种方法的缺点是经过平滑处理后, 得到的网格模型会比原来的网格模型体积变小, 并且新的模型会出现过平滑问题, 也就是原有的尖锐的特征会消失。为了克服这两个问题, Wu等提出一种基于梯度场的平滑方法, 该方法区别于前述的基于法向或顶点的平滑方法, 而是通过求解泊松方程来得到平滑的网格曲面。等提出一种保持原有特征的网格曲面滤波算法, 这种滤波方法的目的在于提高滤波后模型的可信度;Fan等提出一种鲁棒的保特征网格曲面滤波算法, 这种算法基于以下原则:一个带有噪声的网格曲面对应的本原的曲面应该是分片光滑的, 而尖锐特征往往在于多个光滑曲面交界处[2]。

2 点云滤波问题

以上网格曲面光顺算法都需要建立一个局部的邻域结构, 而点云模型中的各个点本身缺乏连接信息, 因此已有的网格光顺算法不能简单的推广到点云模型上来如果仅仅简单地通过最近邻等方式在点云数据中引入点与点之间的连接关系, 那么取得的光顺效果很差。所以, 相对于网格模型来说, 对点云模型进行滤波光顺比较困难, 而且现有针对点云模型的滤波算法也较少。

逆向工程中广泛采用的非接触式测量仪为基于激光光源的测量仪。其优点在于能够一次性采集大批量的点云数据, 方便实现对软质和超薄物体表面形状的测量真正实现“零接触力测量”。在激光测量仪数据采集过程中, 噪声产生的主要因素是被测物体的位置、物体表面的粗糙度和波纹等反射特性、物体颜色和对比度、环境光照条件和测量系统的误差等。

3 点云滤波方法

对于不同类型的点云数据有着相对应合适的滤波方法。对于扫描线式点云数据和阵列式点云数据来说, 现有的三维点云数据的滤波光顺算法可以从以下两个方面进行分类:根据噪声在各个方向上的扩散方式不同可分为各向同性和各向异性算法;根据算法的复杂度分析可分为基于Laplace算子的方法、基于最优化的方法和简单的非迭代方法。各向同性算法优点是算法简单, 但对噪声和模型的尖锐特征不能加以区分, 在去除噪声的同时, 尖锐特征不能得到保持。为了区分尖锐特征和噪声各向异性算法修改了扩散方程, 在去除噪声的同时, 可保持模型的尖锐特征。这些算法需要计算大量的模型结构信息以保持细节, 因此计算量非常大。其他可以采用的滤波算法有:最小二乘滤波、卡尔曼滤波和平滑滤波等。针对于散乱点云数据, 许多相关学者也进行了深入的研究, 提出了拉普拉斯 (Laplacian) 算子、平均曲率流、移动最小二次曲面等方法。

4 点云滤波中的离群点检测

离群点检测问题越来越受到重视, 出现了很多有效的算法, 本节将分别介绍基于统计的方法, 基于距离的方法和基于密度的方法。

基于统计方法的离群点检测基本思想为:对于整个点集来说, 假设其分布符合某种统计模型, 那么离群点就是点集中那些不符合该种统计模型的一些点。

这种基于统计的方法的前提假设是数据点集基本符合某种统计模型, 离群点严重地偏离这个统计模型。因此基于统计的方法具有很多缺点。首先, 这种算法鲁棒性不强, 均值$mu$和协方差矩阵$Sigma$的计算受离群点影响很大, 所以得到的模型不一定真实地反映数据分布情况。之后再用这个误差很大的分布去判断离群点, 必将导致判断的不准确。该方法的第二个缺点是对数据分布的先验知识过于依赖, 如果预先不能正确设定数据分布模型, 则检测结果会有极大的误差。

基于距离的离群点检测方法基本思想为, 通过一个点与其周围相邻点的距离来判断该点是不是离群点。这种判断方法的最基本假设是正常数据点周围存在许多距离较近的数据点;离群点远离它们的最近的邻居。

这种基于距离的方法当数据点集的密度十分不均匀时, 就会得到错误的检测结果。

考虑到基于距离的方法在数据点集的密度不均匀时容易导致离群点检测错误breunig等提出一种基于密度的离群点检测方法。该方法与基于距离的方法根本区别在于, 一个数据点的是否为离群点, 不再是一个布尔型的属性, 而是一个介于0和1之间的有理数, 这个数值越大, 该点就更可能是离群点。这个数值被称为“局部离群系数” (local outlier factor, LOF) 。

5 结论与展望

点云数据滤波问题在逆向工程、三维重建等问题中具有重要的意义。本文研究了点云数据滤波问题及离群点检测问题分析了离群点检测的三种方法。

点云数据的滤波和离群点检测问题还有很多难点没有克服, 未来的研究方向主要包括: (1) 根据目前的点云数据量越来越大, 考虑到现有计算机系统的效率, 寻找快速有效的滤波方法具有重要的意义。 (2) 随着kinect等设备的应用, 动态点云数据的滤波问题显得越来越重要。对动态点云数据开展滤波和离群点检测研究具有重要意义。 (3) 目前的滤波和离群点检测方法主要表现为无监督学习方法[3], 开展半监督和全监督的滤波和离群点检测算法研究有着广泛的应用前景。

参考文献

[1]Vollmer, J.and Mencl, R.and Mueller, H.Improved Laplacian smoothing ofnoisy surface meshes[J].Wiley OnlineLibrary, 1999, 10.

[2]罗大兵, 高明, 王培俊.逆向工程中数字化测量与点云数据处理[J].机械设计与制造, 2005, 4.

[3]董明晓, 郑康平.一种点云数据噪声点的随机滤波处理方法[J].中国图象图形学报:A辑, 2004, 11.

数据滤波 篇2

外测数据处理中小波滤波的鲁棒性算法研究

为消除外测数据处理中异常值和噪声信号对处理结果的影响,结合数据处理的实际,给出一种基于小波变换的`鲁棒性滤波算法.首先用移动中值滤波算法剔除原始数据中的异常值,然后采用小波系数去噪算法并结合经验雏纳周值滤波算法,抑制数据中的噪声.仿真计算及实际工程应用表明,该算法在保留特征段及有用信息的同时,有效地剔除了异常值,押制了噪声,具有很好的鲁棒性.

作 者：刘海波 LIU Hai-bo 作者单位：91550部队,辽宁大连,116023刊 名：飞行器测控学报 ISTIC英文刊名：JOURNAL OF SPACECRAFT TT&C TECHNOLOGY年，卷(期)：27(6)分类号：V557.5关键词：小波变换 经验维纳滤波 移动中值滤波 异常值

数据滤波 篇3

标准Kalman滤波是用于线性系统的最小均方意义上的最优状态估计。Kalman滤波对非平稳信号具有较强的估计能力,但标准最优Kalman滤波的缺点和局限性是要求知道系统的数学模型和噪声统计特性。所以,在单一尺度上,如果由于多种不确定因素干扰,导致使用不精确的模型和噪声统计,设计出的Kalman滤波器会导致滤波器性能变坏[1]。

多分辨率的数据进行多级别、多层次的处理能够获得更有价值的信息,而这种信息是单一尺度上所无法获得的。通过多尺度分解,在不同的尺度空间上分别刻画目标的特征属性再加以融合,可以有效地降低不确定因素的干扰,从而提高滤波器的滤波性能[2,3,4]。小波变换的多尺度特点非常适合多尺度信号的处理,特别是依靠小波变换的重构算法对低分辨率补充一定信息可以获得高分辨率近似补充。文献[5]中提出基于小波变换,寻找基于全局信息的一步预测值和估计误差协方差矩阵,通过综合不同尺度的信息来降低Kalman滤波对数学模型和噪声统计特性的依赖,这种方法比起单一尺度上的Kalman滤波效果是明显的,但计算量较大,由于算法侧重于横向数据更新,全局最优预测值和误差矩阵对滤波效果的影响变得较大,容易受到干扰。本文就此提出一种基于多尺度Kalman的融合滤波方法。充分利用初始估计序列来有效地降低向量和矩阵维数,减少运算量,把标准Kalman滤波只在单一尺度和时间轴上,对状态估计值和误差协方差进行数据更新,改进为基于小波变换在尺度轴和时间轴上,通过不同尺度的观测数据,纵向和横向双向数据更新,算法侧重于纵向数据更新,充分利用多尺度数据进行多层次融合处理,从而获得更好的滤波效果。

本文主要思想可以概括为:首先,在最细尺度N上通过标准Kalman滤波得到初始估计序列,然后通过小波变换将状态估计值和状态估计误差协方差的分块序列分解到最粗尺度,并根据相应尺度上的观测数据进行数据更新,然后通过小波重构,实现不同尺度上的数据更新,当纵向回归到原始尺度即最细尺度时,再跳至下一个数据块进行,直至滤波结束。在实际中,也可以利用此方法解决多传感器多分辨率数据的融合,更好地利用不同分辨率数据的互补信息,达到更佳的融合效果。

2 多尺度动态系统的描述及Kalman最优滤波

通常情况,动态系统状态方程在尺度N上的描述为

其中:k为离散时间变量,x(N,k)是状态向量,Φ(N,k)是系统状态转移矩阵,w(N,k)是系统噪声,且服从均值为零,方差为Q(N,k)的正态分布,初始值x(N0,)为一随机变量,其均值和方差是E[x(N,0)]=x0和E{[x(N,O)-x0][x(N,O)-x0]T}=P0。

动态系统量测方程在尺度N上的描述为

其中:C(N)系统观测阵,z(N,k)是对x(N,k)的观测值,观测噪声v(N,k)服从均值为零,方差为R(N,k)的正态分布,初始值x(N0,)与观测噪声v(N,k)以及系统噪声w(N,k)之间是互不相关的。

在尺度N上标准Kalman最优滤波的基本方程,如式(3)~式(7)所示:

3 基于小波变换的多尺度动态系统描述

我们由小波理论可知[6],通过一个脉冲响应为h(l)的低通滤波器可以从尺度i上获得粗尺度i-1上的平滑信号xL(i-,1k),脉冲响应为g(l)的高通滤波器可获得细节信号xH(i-,1k):

通过小波变换将状态方程从尺度i分解到粗尺度i-1,我们可以发现不同尺度之间的参数矩阵之间的内在联系,为了计算方便,一般我们取状态方程中的G(N,k)为单位阵,则i-1尺度上系统状态方程为

其中上标i表示结果由尺度i上小波分解得到的,比较式(10)和(11)可得式(12)~式(13):

由式(13)可推导出:

同样,用小波变换将量测方程从尺度i分解到尺度i-1,则i-1尺度上系统量测方程为

比较式(15)和式(16)可得式(17)~式(18):

从式(12)∼式(13)和式(17)∼式(18)可以看出,通过尺度i上的系统状态和量测方程的参数矩阵,可以很容易地获得尺度i-1上的相应的参数矩阵,研究标准Kalman最优滤波的基本方程可知,在尺度i-1上进行标准Kalman滤波估计还需要和Pk+1/k(i-1),我们将在多尺度Kalman滤波算法中给出计算方法。

4 基于小波变换的多尺度Kalman滤波

首先,假设一长度为T的信号,我们在单一尺度N上用标准Kalman滤波根据观测值进行估计。由于多种不确定因素干扰,滤波效果并不会很理想,在尺度N上,经Kalman滤波后,我们将得到估计序列和相应的估计误差协方差阵序列Pk/k(N),序列长度均为T,我们把得到的这一系列估计序列定义为初始估计序列,将这些序列分割成长度为M的数据块,其中T和M均设定为2的整数次幂,方便后续的小波变换。这样分割后的序列可表示为Xm(N)和Pm(N),则:

对于数据块Xm(N)中的每一个数据元素均来自于初始估计序列的估计值,由于我们要在尺度轴上对此估计值进一步修正,所以我们可以将此值看作一步预测状态估计值,为了保持与标准Kalman最优滤波方程中的符号的一致性,我们用符号来代替数据块Xm(N)的表示,此时是由初始状态估计值组成的数据块。

由式(8)~式(9)可知:经过从尺度N到尺度i的小波塔式分解,信号可分解为尺度i上的平滑信号和相应的各尺度l(i≤l≤N-)1的细节信号,这样我们把初始状态估计序列通过小波变换向尺度i分解,可得到尺度i上的平滑信号和相应各尺度l上的细节信号,如式(21)~式(22)所示:

对于初始状态估计序列,分解后不同尺度上的序列集合为

接下来,我们要获得相应尺度上的预测误差协方差阵。由估计误差协方差阵的定义可知:

其中:为尺度N上的一步预测误差,由于我们是对长度为M的数据块操作,所以此时估计误差协方差为一矩阵块。根据初始估计误差协方差阵序列,我们可以获得数据块Xm(N)的估计误差协方差阵为,如果用表示第i行第j列的元素,则:

同样,为了保持与标准Kalman最优滤波方程中的符号的一致性,我们用符号Pk/k-1(N)来代替数据块表示。这样对Pk/k-1(N)做小波变换实际上是一个二维离散小波变换[6],通过小波变换我们可以得到尺度i上的Pk/k-1(i)。

若在单一尺度上测量,我们可以通过下采样Zk(N)获得尺度i上的测量序列Zk(i),若是多传感器多分辨率系统,则尺度i上的测量值Zk(i)可以通过不同传感器的观测获得。在尺度i上,我们利用Kalman滤波方程对数据块进行数据更新。由式(17)、(18)、(23)和(26),我们可以得到尺度i上的Ck(i),Rk(i),和Pk/k-1(i),并通过式(27)~式(31)在尺度i上Kalman滤波:

数据更新结束后,通过小波重构我们可以得到尺度i+1上的以及Pk/k-1(i+1),在通过尺度i+1上观测值Zk(i+1)用同样的方法对序列进行数据更新,以此类推,直至尺度N,我们便得到了尺度N上的。同时在小波重构时,还可根据实际的研究情况,结合一定的小波域值去噪方法来去除部分噪声影响[7,8]。

5 仿真算例

实验一

系统在尺度4上的状态模型和尺度j(j=4,3,2,1)上的量测模型,如式(32)~式(33)所示:

其中:、w(,4k)~N[,0Q(,4k)]、v(j,k)~N[,0R(j,k)];式中Q(,4k)=diag1(1,),测量矩阵C(j,k)是单位矩阵,测量误差协方差矩阵R(j,k)为

尺度4为最细尺度,初始值x0=[10,]0T和P0=diag(4,4),选用d B4小波分解和重构。图1(a)信号在尺度4上的量测值(T=512,含两个状态分量)。图1(b)为最细尺度(i=4)上Kalman滤波结果。图1(c)为融合4层信息后的多尺度Kalman滤波结果(M=32),可以很直观地看出,多尺度Kalman滤波的有效性。

从表1中可以看出,融合的尺度越多,测量值的改进越明显,但当融合尺度数目过多时,测量值的改进程度明显减少且运算量随之增大,这是由于测量值提供的有用信息变少了,在粗尺度上的数据更新也就相应变少。从表2中可以看出,不同数据块大小对滤波效果也有明显影响,当数据块选取过大时状态估计误差协方差矩阵较大,数据之间的影响变大,导致滤波性能降低。所以,我们可以确定实际需求,确定适当融合尺度数目和数据块大小有效地提高滤波性能和效率。

实验二

继续采用式(32)∼式(33)的系统状态模型和量测模型,尺度4为最细尺度。假定此系统为多传感器多分辨率系统,我们用三个分辨率不同的传感器来观测系统信号,并进行数值仿真,三个传感器对信号观测的结果如图2(a)∼图2(c)所示,其信噪比(SNR)分别为13.4d B、9.2d B、6.7d B,图中为信号的第一个状态分量,初始值x0=[10,0]T和P0=diag(4,4)。这样式(27)的尺度i上的测量值Zk(i)可以通过不同传感器的观测获得,我们选用d B4小波分解和重构,传感器1、2、3的量测分别对应尺度2、3、4上的观测,将三个尺度上的数据融合后的结果图如图2(d)所示,其信噪比为19.3d B,可以看出基于多尺度Kalman的多传感器数据融合滤波算法具有良好的滤波效果,可用于多分辨率多传感器数据融合。

6 结论

本文通过深入分析基于小波变换的动态系统模型,提出一种基于多尺度Kalman的数据融合滤波的方法,利用小波的多尺度特点,把初始估计序列多尺度分解,并进行分层Kalman滤波估计,通过小波重构进行估计融合。该算法将小波多尺度分解和Kalman滤波结合起来,同时还有效地利用了多分辨率的数据信息。通过实验可以看出,该算法对实际中含较强噪声的动态系统状态估计效果较好,同时也能用于多分辨率多传感器数据融合。

摘要：本文通过分析基于小波变换的动态系统模型,提出一种基于小波多尺度的Kalman数据滤波方法,本文利用小波的多尺度特点,把初始估计序列多尺度分解,并在不同尺度层上进行Kalman滤波估计,再利用小波重构来融合各层的估计信息,把标准Kalman滤波只在单一尺度和时间轴上对状态估计值和误差协方差进行数据更新,改进为基于小波变换的尺度轴和时间轴上的双向数据更新,该算法将小波多尺度分解去噪和Kalman滤波相结合,对实际中含较强噪声的动态系统的状态估计效果较好。算法也可用于多分辨率多传感器数据融合。

关键词：多尺度,Kalman滤波,小波变换,数据融合

参考文献

[1]邓自立,许燕.基于Kalman滤波的白噪声估值理论[J].自动化学报,2003,29(1):21-23.DENG Zi-li,XU Yan.White Noise Estimation Theory Based on Kalman Filtering[J].Acta Automatica Sinica,2003,29(1):21-23.

[2]HONG L,CHEN G R,CHUI C K.A Filter-bank-based Kalman Filtering Technique for Wavelet Estimation and Decomposition of Random Signal[J].IEEE Trans.On Circuits and systems II:Analog and Digital Signal Processing,1998,45(2):237-241.

[3]WEN X B,TIAN Z,LIN W.Robust Kalman Filter and Smoothing Recursive Estimator for Multiscale Autoregressive Process[C]//ICSP’04.7th International Conf.Beijing:IEEE,2004,1(1):364-367.

[4]雷明,韩崇昭,元向辉,等.非线性系统的小波分频的扩展Kalman滤波[J].光电工程,2006,33(6):68-72.LEI Ming,HAN Chong-zhao,YUAN Xiang-hui,et al.Wavelet Decomposing based Subbandwise Extended Kalman Filtering for the Nonlinear System[J].Opto-Electronic Engineering,2006,33(6):68-72.

[5]文成林.多传感器单模型动态系统多尺度数据融合[J].电子学报,2001,29(3):341-345.WEN Cheng-lin.Multiscale Data Fusion for Multi-Sensor Single Model Dynamic Systems[J].Acta Electronica Sinica,2001,29(3):341-345.

[6]MALLAT S.A Wavelet Tour of Signal Processing(Second Edition)[M].San Diego:Academic Press,1999.

[7]彭玉华.一种改进的小波变换阈值去噪方法[J].通信学报,2004,25(8):119-123.PENG Yu-hua.An Improved Thresholding Method in Wavelet Transform Domain for Denosing[J].Journal on Communications,2004,25(8):119-123.

数据滤波 篇4

课程 数字图像处理 名称 数字图像变换

实验名称: 邻域平均法(box模板)和中值滤波处理

一、实验目的

图像变换是数字图像处理中的一种综合变换，如直方图变换、几何变换等。通过本实验，使得学生掌握两种变换的程序实现方法。

二、实验任务

请设计程序，分别用邻域平均法，其模板为：

111

111*1 9111 和中值滤波法对testnoise图像进行去噪处理（中值滤波的模板的大小也设为3×3）。

三、实验环境

本实验在Windows平台上进行，对内存及cpu主频无特别要求，使用VC或者MINGW（gcc）编译器均可。

四、设计思路

介绍代码的框架结构、所用的数据结构、各个类的介绍（类的功能、类中方法的功能、类的成员变量的作用）、各方法间的关系

试验要求中以给出大致的编程思路和源代码以及代码注释，只有黑框部分需要自己填写。在此不进行赘述。

五、具体实现

实现设计思路中定义的所有的数据类型，对每个操作给出实际算法。对主程序和其他模块也都需要写出实际算法。

注意：源代码中要加上注释。代码：(红色为重点代码)<邻域平均法>(3*3)/*------利用第一次实验课提供的 dhc.h 和 dhc.c文件以获取位图的高 宽 以及从文件头到实际的位图数据的偏移字节数，从而实现对位图实际数据的操作。------*/ #include #include #include #include “hdr.h” /*------定义结构指针------*/ struct bmphdr *hdr;//定义用于直方图变量

unsigned char *bitmap,*count,*new_color;/*------main()函数编写------*/ int main(){ //定义整数 i，j 用于函数循环时的，nr_pixels为图像中像素的个数

int i, j ,nr_pixels,nr_w,nr_h;//定义两个文件指针分别用于提取原图像的数据和生成直方图均衡化后的图像

FILE *fp, *fpnew;//定义主函数的参数包括：输入的位图文件名和输出的位图文件名,此处内容可以不要，在DOS下执行命令的时候再临时输入也可，为了方便演示，我这里直接把函数的参数确定了。

// argc=3;// argv[1]=“test.bmp”;// argv[2]=“testzf.bmp”;

//参数输入出错显示

/* if(argc!= 3){ printf(“please input the name of input and out bitmap filesn”);exit(1);}*/ // 获取位图文件相关信息

// hdr = get_header(argv[1]);hdr = get_header(“testnoise.bmp”);if(!hdr)exit(1);//以二进制可读方式打开输入位图文件

fp = fopen(“testnoise.bmp”, “rb”);if(!fp){

printf(“File open error!n”);

exit(1);} // 文件指针指向数据区域

fseek(fp, hdr->offset, SEEK_SET);

//计算位图像素的个数 nr_pixels = hdr->width * hdr->height;nr_w = hdr->width;nr_h = hdr->height;bitmap = malloc(nr_pixels);new_color = malloc(nr_pixels);count = malloc((nr_w+2)*(+nr_h+2));//读取位图数据到bitmap中

fread(bitmap, nr_pixels, 1, fp);fclose(fp);

//因为图像边缘无法使用邻域平均，所以根据邻近颜色填补图像的周围一圈，存入count[]数组中

//中心图像存入count[] for(i=nr_w+3;i<(nr_w+2)*(+nr_h+2)-nr_w-3;i++){

j=i/(nr_w+2);

if(i%(nr_w+2)!=0&&(i+1)%(nr_w+2)!=0)count[i]=bitmap[i-nr_w-1-2*j];} //填补第一排

for(i=1;i

count[i]=bitmap[i-1];} //填补最后一排

for(i=1;i

count[(nr_w+2)*(nr_h+1)+i]=bitmap[nr_w*(nr_h-1)+i-1];} //填补左边一排

for(i=0;i

count[i*(nr_w+2)]=count[i*(nr_w+2)+1];} //填补右边一排

for(i=0;i

count[(i+1)*(nr_w+2)-1]=count[(i+1)*(nr_w+2)-2];}

//邻域平均3*3 for(j=nr_w+3,i=0;j<(nr_w+2)*(+nr_h+2)-nr_w-3;j++){

if(j%(nr_w+2)!=0&&(j+1)%(nr_w+2)!=0)

new_color[i]=(count[j]+count[j-1]+count[j+1]+count[j-nr_w-2]+count[j-1-nr_w-2]+count[j+1-nr_w-2]+count[j+nr_w+2]+count[j-1+nr_w+2]+ count[j+1+nr_w+2])/9,i++;}

//结果存入bitmap[]中

for(i = 0;i < nr_pixels;i++)

bitmap[i]=new_color[i];

// 打开一个以输出文件名命名的文件,设为可写的二进制形式

fpnew = fopen(“test_lynoise.bmp”, “wb+”);

//由于位图文件的头部信息并没有因直方图均衡化而改变，因此输出图像的头部信息从原位图文件中拷贝即可：

fwrite(hdr->signature, 2, 1, fpnew);fwrite(&hdr->size, 4, 1, fpnew);fwrite(hdr->reserved, 4, 1, fpnew);fwrite(&hdr->offset, 4, 1, fpnew);fwrite(&hdr->hdr_size, 4, 1, fpnew);fwrite(&hdr->width, 4, 1, fpnew);fwrite(&hdr->height, 4, 1, fpnew);fwrite(&hdr->nr_planes, 2, 1, fpnew);fwrite(&hdr->bits_per_pixel, 2, 1, fpnew);fwrite(&hdr->compress_type, 4, 1, fpnew);fwrite(&hdr->data_size, 4, 1, fpnew);fwrite(&hdr->resol_hori, 4, 1, fpnew);fwrite(&hdr->resol_vert, 4, 1, fpnew);fwrite(&hdr->nr_colors, 4, 1, fpnew);fwrite(&hdr->important_color, 4, 1, fpnew);if(hdr->offset > 54)

fwrite(hdr->info,(hdr->offset54), 1, fpnew);

//直方图均衡化的数据(bitmap)赋值 fwrite(bitmap, nr_pixels, 1, fpnew);//关闭

fclose(fpnew);//释放内存（优化程序必需）free(hdr);free(bitmap);free(new_color);free(count);return 0;

得出实验结果图像后，比较这两种方法去噪的效果好坏，并分析具体原因。

通过比较，邻域平均法在降低噪声的同时，会使图像产生模糊；中值滤波在降低噪声的同时，能够较好的保持图像边缘，不过我感觉有点水粉画的效果，让图像有点失真。

完成上述工作后，使用程序进行验证分析：使用邻域平均法时，3×3和5×5模板大小对图像进行处理的效果有何差别？并分析原因。

程序及源代码都放在文件夹中，不放进报告了。上面是对比图，显然5×5模板比3×3模版模糊，因为使用邻域平均法时，模版尺寸越大，则图像模糊程度越大。

六、心得体会

这次实验是对两种去噪方法的比较。而在书本中，我们已经看过两种去噪方法的效果了，通过上机验证，果然如此。

基于卡尔曼滤波的陀螺仪数据处理 篇5

Kalman滤波是一种高性能的递归滤波器, 最大的特点是可以在一个不完全、甚至包含噪声的测量时间序列里估计出动态系统的状态。Wiener滤波是以当前和过去全部的观测值为依据, 对信号的当前状态进行估计, 它的解是以均方误差最小为目的的系统的传递函数或单位脉冲响应。而Kalman滤波器不需要过去全部的测量值。它是根据前一个估计值和最近一个测量值来估计信号的当前值, 它是用状态方程和递推方法进行估计的, 因而Kalman滤波对信号的平稳性和时不变性不做要求。[1]

1 Kalman滤波器

美国科学家Wiener和苏联科学家Kолмогоров等人是最早开始研究最佳线性滤波理论的, 后来, 人们把它称之为Wiener滤波理论。从理论上说, Wiener滤波有一最大缺点, 那就是它一定要用无限长的历史数据, 但是, 这在实时处理中往往是不能接受的。为了解决这一问题, 在60年代Kalman在滤波理论中引入了状态空间模型, 同时设计出了一套递推估计算法, 称为Kalman滤波理论。Kalman滤波是以最小均方误差为估计的最佳准则, 再寻求一套递推估计的算法, 从而实现状态的预测。它的基本思想为:采用信号与噪声的状态空间模型, 通过对前一时刻的估计值和现时刻的测量值完成对状态变量估计的更新, 求出现时刻的估计值。由于Kalman滤波对存储量和运算量的要求较小, 因此它适合于实时处理和单片机运算。

Kalman滤波是以线性代数和隐马尔可夫模型为基础的。可以用一个马尔可夫链表示它的最基本动态系统, 建立在一个含有高斯噪声的线性算子上。可以用一个实数向量来表示系统的状态。这个线性算子随着离散时间测量点的增加, 就会作用在当前状态上, 产生一个新的状态, 同时会把一些噪声带入, 并把系统的一些已知的控制器的控制信息加入。而另一个受噪声干扰的线性算子产生出这些隐含状态的可见输出。为了应用Kalman滤波器从一系列有噪声的观察数据中估计出被观察过程的内部状态, 必须要建立Kalman滤波框架下的模型。如下:

其中:Fk表示作用在xk-1上的状态变换模型;Bk表示作用在控制器向量uk上的输入-控制模型。wk表示过程噪声, 假定噪声的符合均值为零, 协方差矩阵为Qk的多元正态分布。

对于任意时刻k, 真实状态xk的一个测量zk满足下式:

其中Hk表示观测模型, 它把真实状态空间映射成观测空间, vk是观测噪声, 均值为零, 协方差矩阵为Rk, 且服从正态分布。

初始状态以及每一时刻的噪声{x0, w1, ..., wk, v1...vk}都认为是互相独立的。但是, 很多真实世界的动态系统都并不确切的符合这个模型;但是由于Kalman滤波器被设计在有噪声的情况下工作, 所以这个滤波器在这种情况下也能得到近似的结果, 所以这个滤波器的到了广泛的应用。

2 陀螺仪ADXRS150

ADXRS150是一款角速度范围为150°/S的MEMS角速度传感器。具有Z轴响应、宽频、抗高振动、噪音为0.05°/S/sqrt Hz, 2000g冲击耐受力, 温度传感器输出, 对精确应用绝对速率输出, 5V单电压操作等特点。由于其小而轻 (<0.15平方厘米, <0.5克) 常应用于辆底盘滚转传感, 惯性测量单元IMU, 平台稳定, 无人机控制, 弹道测量等场合。由于MEMS陀螺仪本身的特点, 其输出数据在长时间内温度漂移、零点漂移都相当严重, 必须使用滤波算法去除噪声干扰。由于加速度是一个瞬间变化的物理量, 对于这种长时间的温度漂移和零点漂移, 其他的滤波算法无法保证相应的及时性。所以只能选用具有估计特性的Kalman滤波器。[2]

3 Kalman滤波程序实现及处理效果

现设线性时变系统的离散状态防城和观测方程为:

其中:X (k) 和Y (k) 分别是k时刻的状态矢量和观测矢量;F (k, k-1) 为状态转移矩阵;U (k) 为k时刻动态噪声;T (k, k-1) 为系统控制矩阵;H (k) 为k时刻观测矩阵;N (k) 为k时刻观测噪声。则卡尔曼滤波的算法流程为:第一步, 预估计X (k) ;第二步, 计算预估计协方差矩阵;第三步, 计算Kalman增益矩阵;第四步, 更新估计;第五步, 计算更新后估计协防差矩阵;不断重复以上步骤即可。[3]以上运算中第一次计算中用到的第一次估计得数据对滤波的稳定会有一定影响, 选择好的话会有较好的滤波效果。图1为滤波处理后的数据和滤波前数据的比较。

4 结语

从滤波的图像可以看出, Kalman滤波能较好地滤除陀螺漂移中的随机噪声。该实验程序的算法部分的运算开销主要在矩阵的转置、相乘上, 序列的长度也会影响速度, 整体运算复杂度为О (n, 3) , 相对来说不高, 适用于目前流行的各种高速8位单片机使用。

参考文献

[1]郑天明.船舶智能监控导航技术研究.华南理工大学硕士论文.

[2]熊海林, 邓方林.陀螺静态漂移系数的两种估计方法[J].上海航天, 2002, 1:39-42.

数据滤波 篇6

在电网数据采集系统中,无论是测量用还是保护用,我们所关心的数据大多数都是工频分量,因此可以在采样之前将最高信号频率分量限制在一定频带之内,即限制输入信号的最高频率。这样一方面降低了对硬件的处理速度的要求;另一方面,根据奈奎斯特采样定律可知,当以频率fs对输入信号进行采样时,能够正确还原的信号的最高频率为fs/2,对于频率超过fs/2的输入信号,如果不在A/D转换之前加以滤除,都将以“合法”的身份呈现在0到fs/2频谱范围内,从而造成频率混叠,消除这种混叠现象的唯一办法是在信号到达A/D转换器之前先对其进行低通滤波,滤除其中的高频噪声和一些峰值噪声[1,2,3,4]。

1 无源低通滤波器

传统方式的低通滤波器设计常采用无源器件,利用电阻和电容来实现,图1(a)为简单低通滤波器的实现电路。它的电压传递函数为:

undefined。 (1)

截止频率fc为:

undefined。 (2)

其频率响应如图1(b)所示,虽然这种电路能执行基本的滤波功能,但在实际应用中存在如下几个缺点:①输出阻抗很高,负载效应大,当输出端有负载时,截止频率fc可能会改变;②最大增益将降低到小于1。

2 有源低通滤波器

与无源滤波器相比,有源滤波器有2个突出的优点[5]:①最大增益或传递函数的最大值可以大于1;②负载效应小,这意味着滤波器的输出响应基本上与滤波器驱动的负载无关。

2.1 单极点低通滤波器

图2(a)为单极点有源滤波器的实现电路。该滤波器的输入阻抗很高,等于运放的输入阻抗,对于CMOS运放,其输入阻抗在1013数量级;与无源滤波器相比,该滤波器电路的输出阻抗很低,典型值为几十欧姆至几百欧姆。它的电压传递函数为:

undefined。 (3)

截止频率fc为:

undefined。 (4)

其频率响应如图2(b)所示,在低于fc的频段内,直流增益A=1+R2/R1,也就是说该滤波器在滤除噪声的同时还能将其通带内的信号按指定比例放大,这一点很有实际意义。在实际应用中,为了减少输入偏置电流的影响,往往在同相输入端串联一个补偿电阻Rc[6,7]:

Rc=(R1‖R2)-rs 。 (5)

其中:rs为信号源阻抗。

2.2 双极点Butterworth低通滤波器

图3为双极点Butterworth低通滤波器,直接采用频域分析方法得到其传递函数为[8,9]:

undefined。 (6)

截止频率fc为:

undefined。 (7)

Butterworth低通滤波器是幅值最大的平坦滤波器,传递函数的设计应使其幅值在滤波器通带内尽可能地平坦。通过取传递函数对频率的导数,使导函数在通带中有尽可能多的等于0的点,这些等于0的点就是低通滤波器的零点频率。直流增益A=1+Rf/R1,但是为了使传递函数的幅值在滤波器通带内尽可能地平坦,A只能取undefined,所以在实际应用中存在很大的局限性。

3 设计实例

在微机综合保护系统设计中,为了提取工频分量,并滤除高次谐波分量的影响,需要在A/D转换之前添置模拟低通滤波器。本设计实例分别选取单极点有源滤波器和双极点Butterworth滤波器来实现所需的设计要求。

对于周波24点采样系统,即采样频率fs=1.2 kHz的数据采集系统来说,只要滤波器的截止频率fc满足fc≤600 Hz,即可滤除不需要的信号。为了达到更好的滤波效果,本实例选取fc=100 Hz,直流增益A的选取可以根据实际需要设定,在此为方便比较两种不同滤波器的性能,选取直流增益A=1.6。

3.1 选用单极点有源低通滤波器

由fc=100 Hz,A=1.6可得:

undefined

选C2=1 μF,则有:R2=1.59 kΩ≈1.6 kΩ,R1≈2.7 kΩ。在忽略信号源阻抗的前提下,补偿电阻Rc=R1‖R2≈1 kΩ。

该滤波电路的幅频特性和相频特性如图4所示。

3.2 选用双极点Butterworth低通滤波器

由fc=100 Hz,A=1.6可得:

undefined

选C=1 μF,Rf=6 kΩ,则有:R≈1.6 kΩ,R1=10 kΩ。该滤波电路的幅频特性和相频特性如图5所示,单位阶跃响应曲线图如图6所示。

假若fc=100 Hz保持不变,把直流增益A改为2.3,则易得在该条件下,单极点滤波器与Butterworth滤波器的单位阶跃响应曲线图,见图7。

分析比较单极点滤波器与双极点Butterworth滤波器的Bode图及单位阶跃响应曲线图可知:①单极点滤波器的系统稳定性更好,但幅频特性较差;②双极点Butterworth滤波器幅频特性较好,通带更为平坦,并且在过渡带的衰减比单极点滤波器快,只是该滤波器系统稳定性较差;③当增益A改变时,双极点低通滤波器的系统超调量增大,阶跃响应的振荡性变强,稳定性明显变差,而单极点低通滤波器的阶跃响应特性基本不受增益A的影响,仍然保持良好的系统稳定性。

综合上述分析可以得出:单极点低通滤波器在满足一定滤波效果的同时,又可通过调节电阻R1改变直流增益A的值,以达到我们所需的输出信号范围。凭借这一特点,单极点低通滤波器在电网数据采集中得到了广泛的应用。

4 结论

在数据采集系统中,模拟低通滤波器是关键的组成部分,用以滤除或衰减高频噪声的影响,提高系统的精度。本论文先从理论角度剖析了无源滤波器与两种常见有源滤波器的优缺点,然后给出两种常见有源滤波器的设计实例,通过详细地对比分析得出:单极点低通滤波器的适用范围更广,更具有现实应用价值。

参考文献

[1]Bonnie Baker.嵌入式系统中的模拟设计[M].李喻奎,译.北京:北京航空航天大学出版社,2006.

[2]彭永胜,王太勇.测控系统中的低通滤波技术[J].精密制造与自动化,2003(S1):108-110.

[3]唐振.模拟低通滤波器的设计[J].安徽水利水电职业技术学院学报,2008(1):56-57.

[4]黄知涛,郑龙席.高性能数据采集系统中信号的低通滤波原理及实践[J].测控技术,1999(5):55-56.

[5]Donald A Neamen.电子电路分析与设计[M].赵桂钦,卜艳萍,译.北京:电子工业出版社,2003.

[6]Theodore F Bogart Jr,Jeffrey S Beasley,GuillermoRico.电子器件与电路[M].蔡勉,王建明,孙兴芳,译.第6版.北京:清华大学出版社,2006.

[7]Li Zhongshen.New study on butterworth low pass filter[G]//Proceedings of the First International Symposiumon Test Automation&Instrumentation.Beijing:BeijingWorld Publishing Corporation,2006:1532-1535.

[8]李钟慎,洪健.基于改进型Butterworth传递函数的高阶低通滤波器的有源设计[J].电子测量与仪器学报,2008(1):86-89.

数据滤波 篇7

金属材料在外力作用下所表现出的各种特征, 如弹性、塑性、韧性、强度等统称为力学性能指标, 是机械设计的主要依据。获得力学性能指标的唯一可靠的方法是进行测试。拉伸试验法就是检验金属材料力学性能的一种极为重要的方法。由拉伸试验所确定的金属力学性能四大指标:抗拉强度δb、屈服强度δs、伸长率δ和断面收缩率ψ最具代表性, 是设计制造的主要依据[3]。

在金属拉伸实验中, 通过压力传感器将材料所受的拉力经A/D采样得到离散的载荷数据, 将这些离散数据拟合成载荷曲线, 并从曲线中按照一定的原则找到金属的性能指标。采样数据由于各种原因会迭加各种干扰, 为了满足测量数据的准确性和快速性, 必须采用适当的方法对采样数据进行快速数字滤波处理。目前常用的方法有傅里叶变换数字滤波、中值平均滤波、移动平均滤波等方法, 这些方法在处理金属拉伸实验数据时并不是最好的, 针对这一问题, 提出了采用小波变换对采样数据进行滤波处理的方法, 取得良好效果。

1 拉伸实验数据小波变换分析方法

1.1 拉伸实验测量数据

图1为温度与拉伸速度均恒定的条件下, 圆型低碳钢沿轴向拉伸时得到的位移-载荷典型曲线。曲线由4部分组成, (1) 弹性阶段:oa段, 位移与载荷呈线性关系。载荷在a点突变, 称为上屈服限, 记为σsu。 (2) 屈服阶段:ac段, 材料出现微量塑性变形, 即载荷增加不大, 但位移出现较大变化, 这一现象称为屈服现象。不计初始瞬时效应b点的最低点, bc段载荷最低点b′称下屈服限, 记为σsl。 (3) 塑变阶段:cd段, 随着塑性变形的增加金属材料产生硬化, 材料的抵抗力持续增加, 在d承受的载荷达到最大, 称为抗拉强度限, 记为σb。 (4) 局部变形阶段:de段, 材料出现“缩颈”现象, 随着位移的增加, 到达e点时材料被拉断, e点称为断裂强度, 记为σe。金属的主要性能指标的确定方法见表1[2]。

金属材料在拉伸实验过程中, 受机械振动、夹具滑溜等原因会使采样的数据迭加强烈的干扰信号, 通过对实际数据进行分析, 干扰信号相当于最大幅值为0.6KN的白噪声, 如图2 (a) 所示。由表1中计算金属力学性能指标的方法可知, 求取金属性能指标所需的几个数据均处于实验数据曲线的突变点位置, 为了精确得到这些数据, 需要选择合适的滤波方式。采用的滤波器不但要能滤除信号中的噪声, 最关键的问题是在对信号滤波后能保留采样数据曲线中的这些突变点特征。

1.2 滤波方法选择

根据干扰信号的统计特征和频谱分布的规律, 可以有很多种滤除信号噪声的方法。传统的方法是根据噪声能量一般集中于高频, 而信号频谱分布于一个有限区间的特点, 采用傅里叶变换将含噪的时域信号变换到频域, 然后采用低通滤波器进行滤波。这种方法在有用信号与噪声的频带相互分离时, 滤波效果十分明显;但当有用信号和噪声的频带相互重叠时, 比如噪声为白噪声, 滤波效果则不尽理想。因为低通滤波器在滤除噪声的过程中, 也同时将信号的突变部分变的模糊, 丢失了有用信号的特征;而高通滤波器虽然可以使突变部分得到清晰表现, 但高频量噪声却不能被很好滤除。因此, 传统的傅里叶变换滤波方法存在着保护信号局部性特征和滤除噪声之间的矛盾。

小波变换具有自适应的时-域局部化功能, 在信号的突变部分, 某些小波分量表现幅度突出, 它与噪声的高频部分的均匀表现正好形成明显的对比, 因此小波分解特别是正交小波分解能有效的区分信号中的突变部分和噪声, 达到消噪的效果, 为解决上述的矛盾提供了有力的工具。小波去噪之所以取得成功是因为小波变换具有以下重要特点: (1) 时频局部化特性, 小波变换可在时间轴上准确的定位信号的突变点位置; (2) 多分辨性质, 由于采用多尺度分解方法, 可以很好表现信号的非平稳特点, 如边缘、尖峰、突变点等, 有利于特征提取和保护; (3) 去相关性, 小波变换具有很强的去数据相关性, 使信号的能量集中于少数几个小波系数上, 而噪声能量分布于大部分小波系数上; (4) 小波基选择多样性, 有各种各样的小波基可供选择, 针对不同应用场合选用最合适的小波函数, 可以使去噪性能达到最优。

对图1 (b) 所示含噪信号采用小波变换滤波, 就是利用了小波变换优良的时域局部化能力, 既能滤除干扰, 又能保留突变信号。

1.3 小波阈值去噪方法

信号与噪声在小波域有不同的性态表现, 它们的小波系数随尺度变化的趋势不同。随着尺度的增加, 噪声系数的幅值很快衰减为0, 而有用信号系数的幅值基本保持不变。采用阈值滤波方法的实质在于减小甚至完全剔除由噪声产生的系数, 同时最大限度保留有用信号的系数, 最后由经过处理后的小波系数重构原信号, 得到真实信号的最优估计。“最优”的精确定义依赖于具体的应用要求。

小波变换阈值滤波的步骤一般为:

(1) 首先对一维含噪信号进行小波分解, 选择一个合适的小波基并确定具体的分解层数, 然后进行分解计算;

(2) 对小波分解高频系数进行阈值量化, 对各个分解尺度下的高频系数选择一个软阈值或硬阈值进行量化处理;

(3) 一维小波重构, 根据小波分解的最低层低频系数和经过量化处理后各层细节系数进行一维信号重构。

小波阈值滤波不仅能几乎完全滤除噪声, 而且可以很好地保留反映原始信号的特征信息, 如边缘及突变点信息, 因而具有良好的滤除噪声效果。事实上, 人们已经证明在均方误差意义上阈值法能得到信号的近似最优估计, 且采用软阈值所得到的估计信号至少与原信号同样光滑。

1.4 小波基选择

小波阈值去噪的滤波效果强烈依赖于所选的小波基。经大量分析知道, 如果用和信号形状相近的小波进行滤波时, 会得到较好的滤波效果。具有正交和紧支撑性质的Symlets小波比较适合实时小波滤波要求, 因此选择Symlets小波对实验数据进行滤波处理。在Sym N系列中经常用的是Sym4和Sym8, 而Sym8的正则性优于Sym4, 所以采用Sym8小波基。

2 采用小波变换方法对拉伸实验数据滤波

2.1 对采样数据的处理流程

如图1 (b) 是在拉伸实验过程中经A/D采集的离散数据, 系统采用CS5532作为外接A/D转换器, 该器件应用电荷平衡技术和极低噪声的可编程增益斩波稳定测量放大器, 内部集成了放大和滤波功能, 可得到高达24位分辨率的输出结果, 转换速率最高可达3200Hz。图1 (b) 中数据的采样频率f0=3KHz。无论采用傅里叶变换还是小波变换进行数字滤波都需对采样数据进行实时性处理。用滑动数据窗可解决实时性问题, 其软件实现流程如图2所示。

2.2 用傅里叶变换与小波变换对采样数据进行滤波的结果比较

如图1 (b) 所示, 拉伸实验开始时, 微处理器将采样数据依次放入预先定义的32位宽、长度为128的内存空间 (构成移动数据窗) , 数据采集达到128时采用傅里叶变换或小波变换对窗内数据进行数字滤波, 并输出最后一个滤波值。然后再采样新数据并更新内存空间, 并再次采用傅里叶变换或小波变换对数据窗中数据进行滤波。从而可得到如图1 (c) 所示的傅里叶变换滤波后波形及图1 (d所示的小波变换滤波后波形。

从图1 (c) 处理结果看, 经傅里叶变换滤波后, 含噪信号的滤波效果在信号光滑局域是相当不错的, 但在屈服阶段体现着信号几个关键突变信息也变得非常平坦 (采用中值平均滤波、移动平均滤波等方法时, 会使突变信息变得更加平坦, 这种结论比较容易理解, 不再重述) , 这使得这一阶段中需要的载荷数据采集产生较大误差。从图1 (d) 处理结果看, 经小波变换滤波后的去噪效果在光滑局域甚至弱于傅里叶变换滤波, 但在屈服阶段的突变信息几乎完全得以保留, 这是我们希望看到的结果。

对滑动数据窗中128个采样点进行离散小波分解时其最大分解尺度为7层。从分解得到的细节信号可以看出, 信号在第5层基本达到稳定状态。而随着分解层次的增加, 计算量也相应增加, 因此从尽可能提高运算速度方面考虑, 仅对第5层的尺度系数进行分析。根据小波阈值滤波方法, 将第5层以下的细节系数均置为0, 然后对由第5层低频信号及处理后细节系数对信号进行重构, 就能得到图1 (d) 所示波形。

数据的采样频率为3KHz, 根据Nyquist定理知道采样数据中包含的有效信号最高频率为1500Hz。因此对该信号进行小波分解时, 各层细节信号所占频带第1层为750~1500Hz;第2层为375~750Hz;第3层为187~375Hz;第4层为94~187Hz;第5层为47~94Hz。采用小波变换阈值滤波时, 将第5层以下的细节系数均置0, 这一过程的实质是强制滤除了47Hz以上的所有频率;而采用傅里叶变换滤波时, 先用离散傅里叶变换将时域信号转换至频域, 再将频谱中所有频率范围在47Hz~Nyquist频率与Nyquist频率~2953Hz内的频率全部滤除, 就能得到图1 (c) 所示波形。也就是说, 经傅里叶变换或小波变换滤波后得到的去噪信号中所含频段是相同的。所以, 图1 (c) 与图1 (d) 所示的结果是完全可以进行比较的, 比较结果能真实反映两种滤波方法的优劣。

3 结语

通过傅里叶变换与小波变换对同一含噪信号滤波进行分析, 小波变换的局部时-域分析能力既能滤除噪声, 又能保持信号的原貌。这对于保证拉伸实验屈服阶段信号的特征不丢失, 提高处于实验数据突变位置数据采样的精确度无疑是非常有效的, 这一结论在实践中也得到充分验证。

摘要：小波变换具有良好的时频局部化能力。当采用小波变换对金属拉伸实验数据进行数字滤波时, 不但能滤除干扰信号, 而且能很好的保留信号的突变部分, 体现信号原貌, 这一功能是傅里叶变换所无法比拟的。通过比较傅里叶变换与小波变换对金属拉伸实验数据的滤波结果, 表明采用小波变换处理含噪的实验数据, 具有良好效果。

数据滤波 篇8

关键词：S模式数据链,拉格朗日级数,残差-新息

在空中交通管制系统中,S模式数据链是空- 空数据链和空-地数据链的重要组成部分, 也是反映空中交通态势、保证飞行安全的一项重要手段。 随着国家在空管事业上的大力投入, 以及传统的航管雷达体制[1]存在覆盖范围小及运营成本高等诸多问题,发展具有自主知识产权的S模式数据链技术[2]( 如图1) 成为目前和未来几年中国内各相关研究机构的重要任务。 而在该领域中, 自适应滤波目标跟踪定位技术是关键, 该技术的合理运用,有助于滤除目标飞行器地理位置观测信息中的观测噪声, 提高目标飞行器的航迹稳定度, 对于确保空中交通情景态势的可靠性和准确性起到重要作用。

为适应实际应用需求,可靠性和实时性是行业内对于该项技术的两个核心要求,前者要求所采用的核心算法具有良好的鲁棒性(robust)和收敛性(convergence) ,后者则对算法的优化及复杂度的降低提出了很高要求。 目前,国内自适应滤波目标跟踪定位技术在S模式领域的研究和实际应用还处于起步阶段,而国外相关领域的研究和应用目前主要集中在基于卡尔曼滤波的核心算法上, 虽然能够达到较好的收敛特性和较小的残余方差,但运算复杂度较高,其算法实现过程中需要进行大量的矩阵运算( 即使采用行业内普遍的近似算法, 算法复杂度也是很可观的),这无疑给实时性、成本和功耗带来很大挑战。 即使是基于 α-β-γ 核心算法的研究,也主要集中在如何从残差(relics) 中提取 “ 新息”(innovation) 以及如何建模增益矩阵(Gain matrix) 上, 对该算法本身的稳定性、 初始状态和相关参数的选取与优化未作深入考虑,而这些方面对于飞行器安全是非常重要的。

本文技术探讨的目标有以下三点:

( 1 ) 设计一种适合空中交通管制系统的自适应( 航迹跟踪) 滤波器, 解决当前国内相关技术领域普遍存在的实时性与设计复杂度之间的矛盾;

( 2 ) 对飞行器的三维空间运动轨迹建模 “ 拉格朗日”( Lagrange ) 三阶级数展开, 设计一种兼顾实时性和预测滤波效果的算法模型, 并对该滤波器模型的稳定性、 初始状态、收敛性、滤波参数等进行深入分析;

( 3 ) 提出一种基于 “ 查表” 策略的 “ 残差- 新息” 估计方法, 在保证收敛性的同时兼顾了收敛速度, “ 表” 的维护更新就是SNR实时估计的过程, 可根据样值二阶数学期望的无偏估计统计得出。

1 算法模型设计

为了提高机载GNSS设备的可靠性和稳定性,DO -260A协议[3]明确规定, 具备ADS - B功能的S模式应答机必须使用自适应滤波跟踪算法对本机导航定位信息进行处理,实现对航迹的滤波(消除观测噪声影响)、预测( 丢点时) 和跟踪。 滤波跟踪算法可采用Kalman滤波[4,5],也可以采用 α-β-γ 滤波。

首先建立空间状态转移方程和观测方程如下:

状态转移方程模型的建立源自一个基本理论 ———“ 任何具有有理功率谱密度的随机信号均可以看作由一个白噪声激励某个物理网络形成”[4], 即Z域功率谱密度满足关系式 Φxx( z ) = σw2× B ( z ) × B ( z- 1) , 其中 σw2为激励白噪声的功率谱密度。 式(1) 中: 下标k表示时刻点,yk为状态变量,zk为观测变量, 矩阵Ak为状态转移矩阵,wk - 1为激励变量,Ck为观测矩阵,vk为观测噪声(噪声源)。

定义状态矢量 , 分别表示k时刻的距离、速度和加速度。 假定飞机运动模型(用状态转移方程描述) 为白噪声激励源下的匀加速运动模型, 飞机运动观测环境( 用观测方程描述) 为加性高斯白噪声环境, 利用“拉格朗日级数展开”思想,目标运动模型可视为距离sk的三阶幂级数展开,状态转移方程中wk - 1变量表征飞机机动, 可定义为目标加速度的高阶导数, 时间间隔为T ( 采样或内插周期) 的相邻两次 “ 采样点” 之间视飞行器为匀加速运动( 只要采样率足够, 则可以描述实际运动轨迹)。 定义该飞行器三维运动模型的开环预测矢量和闭环估计矢量为:

k时刻的增益矩阵Hk仅由参数 α、β 和 γ 决定, 且三者之间存在如下关系:

式(3) 中 γ =p ×q,q为滤波增益,p反映表征飞行器机动的 “新息”与观测噪声的 “信噪比”信息,可通过对最近N次残差的方差求数学期望获得。

由于并没有足够的先验统计知识能够区分开新息中激励和观测噪声( 因为残差的均方和平均值指标既反映了测量误差的大小,也反映了目标的机动情况), 应当根据实际观测环境中的噪声数值分析结果( 二阶矩分析或特征值分析) 对式(4) 中的q进行适当取值, 以获得最小均方误差。 从时域角度上看,q值控制跟踪滤波的收敛性;从频域角度看,q值影响滤波器的 “带宽”,也就是说,在机动( 激励信号) 较大时, 如果q取值较小, 滤波器等效带宽较窄,航迹关联会出现失真;反之如果噪声较大,而q取值较大, 滤波器等效带宽较宽, 则航迹关联 “ 抖动” 会变大。

2 关键参数分析

2 . 1 初始值选取

工程上一般用数据通道建立后接收到的机载导航信息作为滤波器初值,渡过滤波暂态后可输出稳态最优值(MSE准则下),如图2 所示。

2 . 2 q值选取与收敛性

自适应滤波采用了一步迭代结构, 是负反馈过程,如何保证滤波器的收敛是保证正常工作的最重要问题之一。 对式(2)中预测方程进行Z域变换:

第k时刻增益矩阵Hk是对角矩阵, 特征值 α、β 和γ 仅由参数p和q确定。 为保证自适应滤波器收敛, 要求式(5)表征的系统传递函数的 “全部复极点包含在单位圆内”[6]。 通常需要实时计算矩阵求逆和特征值,这将给工程应用带来较大困难。 以下给出工程应用时的一种较为简便的保证收敛性的处理方法。

计算信号与观测噪声的比值最为关键,根据Bernoulli大数定理对观测噪声进行数值分析,从大量的观测样本中可以通过二阶中心矩反映观测噪声强度, 再通过 “ 查找表”( 预置的q-SNR收敛性曲线, 可仿真获得) 的方法确定合适的自适应参数,从而保证算法的收敛性。

3 建模仿真

3 . 1 跟踪性能仿真

真实运动航迹为二次曲线, 前半程匀加速运动、 后半程匀速运动。 图3 仿真的是理想观测环境(无观测噪声)中自适应算法的跟踪性能。

从图3(a) 中可以看出, 航迹估计值和真实值重合,图3(b) 和图3(c) 的速度估值曲线和加速度估值曲线则反映了飞机的运动过程,图3(d)的均方误差曲线显示在初始时刻(由静止到启动)和由匀加速运动转为匀速运动的时刻, 会出现 “ 冲击”, 这与图3(b) 的速度估值曲线和图3(c)的加速度估值曲线的分析结果是吻合的。

减小q值, 滤波器等效带宽变小, 不能准确跟踪真实航迹, 在大机动情况下存在跟踪误差, 仿真结果如图4 所示。

3 . 2 滤波性能仿真

真实运动航迹为正弦曲线, 图5 仿真的是有噪环境下自适应算法的滤波跟踪结果,图5(a)对比了观测航迹和滤波跟踪航迹,图5(b)对比了真实航迹和滤波跟踪航迹。 仿真结果说明:通过选择合适的滤波器带宽,可以同时达到滤除观测噪声和跟踪目标机动的目的。

上述仿真结果验证了以下两个结论:

( 1 ) 自适应航迹跟踪滤波器的数学模型是三阶拉格朗日级数展开,因而在理想观测环境中可以由任意估计初值无误差地跟踪机动目标的真实二次运动轨迹(匀加速运动或者匀速运动)。

( 2 ) 从时域角度看, q值影响跟踪滤波的收敛性, 从频域角度看,q值影响滤波器带宽, 在不同噪声强度的观测环境下,q值的选取是否合适直接影响到航迹滤波和跟踪的效果。

在空中交通管制系统S模式数据链领域,国内外当前普遍采用的自适应航迹滤波跟踪技术主要是 “ 基于Kalman滤波模型的MSE准则算法” , 这也是FAA DO -181C 、 DO - 260A等适航标准推荐的信号处理方法。 与之相比, 本文论述关键技术在以下几个方面有其自身特点:

( 1 ) 实现复杂度方面: 相对Kalman滤波模型, 本文采用的算法建模复杂度较低, 这是因为Kalman滤波模型在运算过程中需要牵涉数次矩阵求逆,即使采用逼近简化算法,其运算量也相当可观,与之相比,采用滤波参数替代了 “增益矩阵”,在实现复杂度、运算资源以及系统功耗等方面具备优势。

( 2 ) 收敛性方面: 如果采用矩阵运算求解Kalman滤波模型,则不牵涉收敛性问题(也不存在收敛时间),如果采用传统LMS方法[4]迭代 “逼近”,则同样存在收敛性问题;本文提出的基于 “查表” 策略的 “ 残差- 新息” 估计方法, 在保证收敛性的同时兼顾了收敛速度, 与传统LMS方法解算Kalman滤波模型相比, 在保证收敛性和收敛速度方面具有优势。

数据滤波 篇9

机载激光扫描(Li DAR)系统应用广泛,可快速获取地表的三维点云数据[1],在获取的点云数据中,包含有地面点与地物点,地物点主要包括建筑物、植被、电力线、车辆等信息。为获取地面点得到DEM,需要对点云数据进行处理,滤掉地物信息,该项工作称为滤波。

机载Li DAR数据滤波作为目前Li DAR数据处理的主要问题,国内外已经有许多学者对其进行了研究,如基于多次回波的滤波算法[1~4]、数学形态学的滤波方法[5~8]、迭代线性最小二乘内插法[9]、基于多方向的地面滤波算法[10]、基于TIN的滤波算法[11]等。在多次回波滤波算法中,主要考虑了多次回波特性对地物进行滤除,文献[1,2]主要基于双次回波进行滤波,主要通过双次回波特性首先分离出植被点,再进行建筑物等地物的滤除;文献[3]对森林地区的多次回波进行分析并滤波;文献[4]则通过多次回波特性预先剔除部分建筑物与植被点信息,减少数据量。基于数学形态学的滤波方法主要通过数学形态学的膨胀与腐蚀运算进行滤波,该类算法关键在于结构元素的选取,文献[5]通过改进的形态开运算对未填补空白区域的格网化DSM进行滤波,通过迭代计算实现点云分类;文献[6]考虑实际地形变化,根据形态学知识判断点间的高差与距离来实现地面点与地物点的分类,但是对于高差变化较大或者陡峭的地形具有一定局限性;文献[7]通过改进的形态学梯度计算每个点的梯度,再基于梯度值选取特定的点进行迭代开运算,根据梯度直方图减少迭代的次数,通过判断每次开运算后点的高程与原高程的差值是否小于一定的阈值,逐步滤除非地面点;文献[8]在数学形态学基础上,对格网化DSM进行腐蚀运算获得标记图像,对其反复进行测地膨胀实现开重建,利用白顶帽变换得到n DSM,实现地物与地面点的分类;文献[9]则使用线性最小二乘法内插激光脚点,根据拟合残差不服从正态分布的原则,使用迭代线性最小二乘内插法滤波与内插,并逐渐逼近和修正地形,该算法在地形变化较大的区域滤波时间较长,效率较低;文献[10]提出了基于多方向的地面滤波算法,对Li DAR数据内插出格网点,在扫描线上自动选取最低点作为地面点,通过该点的多个扫描方向的坡度,每一个点与最近邻点的高程差,每一个点与邻域最低高程点的高程差一系列的相关关系来判别数据点的类型,在城市地区与森林地区起到一个比较好的滤波效果;文献[11]采用改进后的三角网迭代滤波算法进行滤波,以待处理格网为中心,将待处理格网最低点与领域格网最低点构成三角网,再对格网内的数据与所在三角形的三个顶点进行比较,实现地面点判别,算法的难点在于不同地形阈值的确定。在这些算法中,部分算法仅考虑了双次回波而忽略了三次及以上的回波,或者适用于森林地区具有一定的局限性,或者需要先进行格网化再进行滤波。

分析Li DAR点云数据的多次回波特性,考虑现有算法思想,提出了考虑多次回波的数学形态学机载Li DAR滤波方法,滤波对象为原始离散Li DAR点云数据,首先提取多次回波中的首次回波与中间次回波作为种子点,对其进行区域增长,滤除掉大部分建筑物与植被,再采用较小结构元素进行形态学滤波,得到地面点,通过内插即可生成DEM。通过实例验证了滤波方法的有效性,具有较高的精确度和较好的实用性。

1 研究方法

1.1 多次回波分析

机载Li DAR系统发射激光脉冲,遇到物体时即返回信号,返回的回波信息包括单次回波与多次回波。同一束激光脉冲只发生一次反射则为单次回波,发生多次反射则为多次回波,其中对于多次回波而言,返回的第一个回波信号称为首次回波,返回的最后一个回波信号称为尾次回波,其余的回波信号皆为中间次回波。在实际中,两次回波之间的距离达到一定间距时,第二次回波才有可能被探测到[3]。目前,现有机载Li DAR系统都具备多次回波信息记录的能力,目前ASPRS激光委员会已制定并不断完善点云数据格式,即LAS格式,LAS数据的主要格式见表1,其中,回波次数为当前脉冲产生的回波总数,回波号表示当前数据为第几次回波。若回波次数为3,回波号为2,表示当前数据为第二次回波,当前脉冲共产生了3次回波。

LAS数据具有多次回波信息记录,分析多回波信息可以判断被测目标的类型:

(1)对于植被密集区域:首次回波大多为树冠部分反射形成,最后一次回波一般是较树冠低一些的枝叶或地面的反射信号[3,12],同时树冠部分或者地面也易形成单次回波信息,但是对于植被区域而言,多次回波数据中的首次与中间次回波一定为植被点。

(2)对于建筑区域,由于建筑物表面较为规则,高差较小,通常产生单次回波,在建筑物边缘区域高差较大,易产生多次回波[4],多次回波中的首次回波通常为建筑物边缘的屋顶点或者墙壁点,中间次与末次回波则可能为地面点或者墙壁点。

通过以上分析,多次回波的首次与中间次回波通常为植被点或者建筑物边界点,对其进行区域增长可以滤除大部分的植被点或者大面积的建筑物。

1.2 区域增长

通过对点云数据的多次回波分析,将多次回波的首次与中间次回波数据作为种子点,执行区域增长算法,即可获取大部分的地物点(主要为较大的建筑或者密集植被区域部分点)。由于滤波方法直接基于原始离散点云数据进行处理,点云密度高、数据量大,在进行区域增长之前,需要对点云进行组织,方式通常为重采样[13],如通过使用一个二维的格网覆盖在原始点云上,判断其落入哪一个格网单元内,每一个格网单元只保留一个最低点,当格网单元内没有点落入时,搜寻邻域最临近点作为该格网单元的点。这种重采样方式的缺点在于:当一个格网存在多个点时,由于只保留一个点,造成了精度损失。

为保持原始点云数据原貌,基于重采样思想[13],设计了一种网格分块链表的方法建立空间索引:

(1)计算原始Li DAR数据的平均间距d,构造一个合适的二维矩阵,矩阵的单元格大小与d相当,对每一单元格设置一空的链表,表示当前格网无点云数据;

(2)遍历Li DAR数据点,判断点云落入的格网单元,并对该格网单元加入链表索引,反复执行该步骤,直至所有点云判断完毕。

建立空间索引后,将空间区域按照平面坐标划分为方形小块,并基于行、列值对小块进行索引,在后期区域增长与形态学滤波时,对某一参考点的搜索可通过设置搜寻半径,限制在该参考点附近的子块当中进行,极大地提高了数据处理效率。建立空间索引后,即可对点云数据执行区域增长算法[14]:

(1)遍历点云数据,获取多次回波中的首次与中间次回波数据,将其作为地物点种子点;

(2)将种子点的高程同邻域点高程进行比较,若邻域点与种子点的高差小于给定的阈值,则标记该邻域点为地物点,并将其作为种子点继续进行增长;

(3)反复执行第(2)步的操作,直至所有点已判别并分类完毕。

1.3 形态学滤波

形态学图象处理是基于集合论对图像进行分析的工具,它着重研究图像的几何结构,基本思想是用具有一定形态的结构元素去量度和提取图像中的对应形状以达到对图像分析和识别的目的[9]。数学形态学的基本运算主要为膨胀、腐蚀、开运算、闭运算,若以A、C表示灰值图,A(m,n)表示图A在(m,n)的灰度,B为结构矩阵,B(i,j)为其在(i,j)的值。基本运算定义如下[10]:

定义1

A被B膨胀,记为AB:

定义2

A被B腐蚀,记为:

式中,表示对所有的B(i,j)取[·]中值的最大值。表示对所有的B(i,j)取[·]中值的最小值,结构元素可以为一维直线或者其他形状。对图像A先腐蚀后膨胀则为开运算,先膨胀后腐蚀则为闭运算。

基于数学形态学的滤波方法关键在于结构元素的大小设置,同时还与Li DAR点云数据内地物分布有关。合适的结构元素窗口执行一次开运算可以滤除掉所有的地物点,但是固定窗口大小的结构元素很难滤除不同类型的地物,结构元素窗口过小,只能滤除掉较小地物点(比如车、单株树木、小建筑物、行人等),大面积的房屋则较难滤除,当结构元素窗口过大,则易造成过度滤除,丢失地形特征信息。由于在多次回波中,较大的建筑物或者植被区域,已经通过区域增长进行滤除,剩余的地物多为较小的建筑物或一些低矮地物,因此,采用较小的结构元素窗口进行形态学滤波,获取地面点,避免了较大结构元素进行形态学滤波的过度滤除问题。

1.4 DEM内插

经过区域增长与数学形态学滤波,建筑物与植被被滤除,得到地面点数据,此时地面点存在大量数据空洞,为得到DEM,还需要进行内插处理,逐点内插法[15]以原始离散点云数据为基础,能较好地保留局部地形,该方法以待内插的数据点为中心,使用一个曲面函数去拟合邻域数据点,逐点内插法拟合函数常采用二次曲面进行拟合:

求解出二次曲面参数后,代入待插点的坐标值,即可求出该点的高程。

2 数据试验与分析

选取城市区域约480m×300m大小的原始LAS点云数据对本文提出的滤波方法进行验证,实验数据采用LAS文件原始格式,LAS版本为1.2,区域总点数为255074个,平均每平方米点数为1.77,点云平均间距为0.75m。该数据最大回波次数为4,产生的多次回波数依次为233948、20603、520、3,主要为首次、二次与三次回波,而四次回波数据只有3个。根据原始Li DAR点云数据生成灰度图像,见图1,从图中可以看出,选取的实验区域被大量建筑物与植被覆盖,且建筑物之间也存在大量树木,建筑物与植被等地物较为复杂,选取的数据较为典型。

在多次回波分析中得出多次回波的首次与中间次回波为植被点或者建筑物边界点,首先提取多次回波数据中的首次与中间次回波数据,提取原则为该数据点的回波次数(Number of Returns)与回波号(Return Number)不相等,且回波次数不为1。根据提取的结果数据生成灰度图像,见图2。对比图1可以看出,多次回波数据中的首次与中间次回波数据主要集中在密集植被区域与建筑物边缘处,与分析一致。

以提取的多次回波数据中的首次与中间次回波数据为种子点,执行区域增长算法,经过测试,高程阈值设置为0.3m,遍历点云数据,根据高程一致性原则进行区域增长,得到地物点,见图3。与图1对比,从图3中可以看出,绝大部分建筑物与植被点已被正确提取出来。

对比图3与图1可以看出,还有部分地物点没有滤除,主要为部分植被点及低矮地物点,对区域增长滤波的结果再采用数学形态学的方法进行滤波,结构元素窗口设置为5m,滤波后,采用逐点内插法内插出DEM,生成的DEM见图4。

从图4中可以看出,建筑物与地面点已经被滤除,得到的DEM平滑,基本保持了地表形状,得到的DEM质量较高,说明提出的考虑多次回波影响的形态学滤波方法较为理想,具有一定的实用性。

3 结论