数据滤波

数据滤波（精选八篇）

数据滤波篇1

在便携式的电子类产品中,触摸屏由于其轻便、占用空间少、方便灵活等优点,已经逐渐取代键盘成为嵌入式计算机系统的输入设备。基于触摸屏的输入系统实际上是由触摸屏、触摸屏控制器、微控制器及其相应软件构成的,本文从手机系统中触摸屏硬件和软件接口入手,分析触摸屏的驱动,采样,滤波及校准原理,及其软件实现方法。

2 触摸屏系统组成

3 采样坐标数据和滤波

当触摸屏被点击后,数据会从DOUT得到,但由于触摸的时间长短我们并不能控制,以及触摸的点可能有变动,所以必须在短时间内采样多点,来判断该点是否是有效点及对SPI口采样的点进行滤波。滤波可以通过延时检测坐标点采样点的方式。

采样点算法如下:

触摸屏接口开始对坐标X和Y的模拟量进行采样,根据试验选取适合的的采集次数。这里使用9次采集,分别记入到ptx[TouchSample]和pty[TouchSample]数组中,TouchSample为采集次数。

为了减少运算量,将ptx[]和pty[]分别分三组取平均值,存储在px[3]和py[3]中。

以处理X坐标为例:

计算以上三组数据的差值:

然后对上述差值取绝对值,所得结果简称绝对差值:

判断上述计算的色对差值是否都超过差值门限,如果这3个绝对差值都超过门限值,判定这次采样点为野点,抛弃采样点,程序返回等待下次采样。其中的差值门限需要根据试验测试得到,本文取值为2。找出其中绝对差值最小的2组数据,再将它们作平均,同时赋值给X坐标,及为最终的采样数据。

采样数据得到后,具体滤波算法如下所述:

假设三次连续采样时刻为T1、T2、T3(T3>T2>T1),采样间隔为10ms。由于采样间隔远小于人的反应时间,所以在前两种操作模式下,如果采样点有效,将T1和T3时刻的采样值作平均。其平均值和T2时刻的采样值比较一般不会大于某个门限,否则判定此次采样点为野点。而对于采样点有很大的跳变。跳变过程中的数据是不稳定的,虽然记入了数据,但被判定成无效的采样点,所以需要在程序中定义一个静态数组x[2]记录相邻的两次采样数据。只有当前后数据持续稳定一段时间,才认为这时的采样点有效。程序中使用的间隔门限经多次实验,设置为2。

在本系统中软件开启一个采样任务,当触摸笔按下后,产生中断,软件处理中断同时新建一个采样10ms定时器来进行多次采样,当发现采样的点大于两个点时开始进行滤波处理,当触摸笔抬起时关闭定时器。

4 坐标校准

系统中采样3点校准。(XR0,YR0),(XR1,YR1),(XR2,YR2)为采样出的物理坐标,(X0,Y0),(X1,Y1),(X2,Y2)为LCD上实际的坐标并且与其物理坐标一一对应,通过下面两个方程解出相应的参数序列(A,B,C,D,E,F,K),该序列及为触摸屏的校准序列,K为公分母。

5 实验结果

在手机软件中,我们以LCD实际坐标点(50,20),(60,80),(128,160)为例,采样的坐标为(2695,3515),(2509,2657),(1201,1513),根据上述校准算法得出的该LCD的校准序列(A,B,C,D,E,F K))为:(-230608,0,836135192,0,-299952,1140234200,4288897),经运算确认,可以得到LCD的比较正确的坐标点。

6 结束语

触摸屏的驱动和滤波是现在很多手机开发中的一个重要模块,本文中的采样和滤波算法在手机开发中取得了良好的效果,在手写识别中,也有较高的识别率,同时整个系统有较强的可移植性,可以方便用户移植到其他平台。

摘要：针对当前手机开发中的触摸屏模块,该文介绍了触摸屏坐标采样,滤波,以及校准的算法。实验表明这些算法能够获取准确的坐标数据,和手写识别效果。

关键词：触摸屏,数据采样,数据滤波

参考文献

[1]周立功.ARM微控制器基础与实战[M].北京航天航空大学出版社,2003

[2]刘锬.触摸屏校准程序设计[N].计算机世界报,2006-04-10.

[3]杨国栋,沈培宏.触摸屏技术发展简介[J].光电技术.2002,43(1):49-55.

[4]王卫京.认识触摸屏[J].微型计算机,2000(12):88-91.

数据滤波篇2

卡尔曼滤波在变形监测数据处理中的应用

通过对某大坝动态变形监测数据的卡尔曼滤波处理和结果分析,并考虑大坝水位和温度影响建立与水位和温度因子有关的`模型,同多元线性回归分析方法和多项式拟合方法对比,表明卡尔曼滤波模型能够实时,快速地处理大量动态变形数据,并能有效地改善动态变形监测数据的精度.

作者：王琦孙华李伟华王晓强 Wang Qi Sun Hua Li Weihua Wang Xiaoqiang 作者单位：中国地质大学,测绘工程系,武汉,430074刊名：工程地球物理学报英文刊名：CHINESE JOURNAL OF ENGINEERING GEOPHYSICS年，卷(期)：6(5)分类号：P258关键词：卡尔曼滤波大坝变形监测数据处理多元线性回归

数据滤波篇3

(1) 基于数学形态学方法的滤波方法[2,3]。此类方法是利用特定的形态学算子对点云数据进行操作, 分类出地面点和地物点, 从而滤波。此类方法具有成熟的理论依据, 且应用于城市地区时滤波效果较好;但滤波结果过分依赖形态学窗口尺寸, 且对非平坦地区的滤波效果不理想。

(2) 基于坡度的滤波方法[4,5]。此类方法首先选取种子点, 再通过比较种子点与邻近点之间的高程变化程度, 即坡度, 来滤除非地面点。基于坡度的滤波算法原理简单, 易于实施, 适应性较好, 滤波结果的精度较高。但是这种方法需要逐点计算坡度值, 由于在离散的点空间搜索邻近节点比较困难, 所以其计算量较大, 速度慢;并且对高植被覆盖区域和断裂地形区域滤波困难。

(3) 迭代加密滤波方法[6,7]。此类算法是选取种子点生成趋势面, 根据每个点到趋势面的残差, 得到点再次拟合的权重, 迭代得到最终的趋势面, 即真实DEM。迭代加密方法适用于密集林区, 但算法计算量较大, 且在有陡坡、低矮植被、断裂线的地区滤波效果较差。

(4) 基于表面的滤波方法[8,9]。此类方法的基本思想是把地面看做是连续且平缓变化的表面, 用带限制条件的参数曲面 (如平面、二次曲面、Snake样条曲面等) 滤波。基于表面的滤波方法能够抑制粗差, 滤波效果较好, 但过于强调适用地形的平缓性, 计算量较大, 且需根据地形频繁的调整阈值。

(5) 聚类分割方法[10]等。聚类分割算法是一种基于特征方式 (如几何、光学、统计等) 对激光脚点进行识别的方法, 根据原始数据组织方式 (如规则格网、TIN、K-D树等) 的不同, 利用图像分割的方法对点云数据进行分割, 再利用分割结果进行聚类分析, 来完成数据滤波。基于聚类分割的滤波方法在设计思想上与其他几种方法存在明显区别, 它不依赖于几何假设描述地形或地物形态, 对特殊的地形数据 (如较大的房屋群、植被) 有较好的滤波作用。

基于最小二乘法的滤波方法是一种迭代加密滤波方法, 由Pfeifer最先提出的[7]。其基本原理是利用最小二乘法拟合点集, 根据残差确定再次拟合曲面时点的权重, 迭代得最终曲面即为真实DEM, 参与拟合的点分为地面点, 其余点为非地面点, 从而完成滤波。此类方法将所有原始点赋予相同的“权”值, 因此不具有先验知识, 地形变化剧烈时滤波效果较差[3]。

针对最小二乘法的不足, 文章提出了一种改进算法。对离散点进行腐蚀, 将腐蚀后的高程与原始高程相比, 选取符合阈值条件的点作为种子点, 并根据腐蚀前后的高程差分配不同的权重, 利用最小二乘原理拟合曲面, 此时的曲面为粗略的地面, 再计算所有点的拟合高程, 与实际高程比较, 将小于高差阈值的点加入种子点, 分配权重并重新拟合曲面, 迭代至种子点集不发生变化。实验结果表明算法能够滤除地物点, 保留地面点, 滤波效果较好。

1 最小二乘法滤波原理

首先将原始数据划分为块, 假设区块内所有的点具有相同权重, 将区块内所有点拟合为一个曲面:

依据最小二乘原理, 使得点的拟合高程值Zi与实际高程值Z之差的平方, 即拟合残差vi之和最小, 式 (1) 中a0至a5为系数, (xi, yi) 为i点的坐标。

初次拟合时, 地物点和地面点是等权的, 一般情况下地面点较多, 地物点较少, 所以地面点与拟合曲面间的距离小于地物点与曲面间的距离。利用残差值确定下一次曲面拟合时点的权值, 残差越大, 权值越小。定权函数为:

式 (2) 中, vi为第i点拟合残差值, pi是第i点权值, a, b为定权参数, g为门限, 均由经验值确定。如果点的残差值极小, 则点是地面点, 再次曲面拟合时权重设为1, 残差值极大则为地物点, 再次曲面拟合时权值设为0。点的残差值除了极大和极小之外, 存在一些中间状态, 可能是低小的地物, 也可能是凸起的地面, 引入参数w, 当残差值处于 (g, g+w) 时, 赋予一定的权值, 权值在0~1之间, 权值越小, 对下次曲面拟合的作用越小。通过多次迭代, 最终的拟合曲面就是DEM模型。

最小二乘算法是在所有的点具有相同的“权”值的条件下, “平差迭代”从而滤波。算法适用于平坦城区、密集林区, 但因为所有原始点“等权”, 不具有先验知识, 故算法用于陡坡地区时效果较差。本文对此处不足进行了改进。

2 改进的最小二乘滤波算法原理

首先利用腐蚀算子选取种子点。腐蚀算子是图像处理中常用的算子, 其定义如式 (3) 。

式 (3) 中A为图像灰度值的集合, B为结构算子, 利用B对A的腐蚀运算, 就是得到所有B完全包含在A中时B的原点位置的集合[11]。直观来说, 对于A中某像素点Ai, 腐蚀是取n张图中对应图像灰度值的最小值作为Ai的值, 其中n为结构算子B的元素个数, n张图是对于B的每个元素B (i, j) , 将原图反向平移B (i, j) , 对应的每个像素点的值减去 (B (i, j) -1) , 平移而得。

本算法的腐蚀运算处理的是离散点的高程值, 没有规则结构, 且不用结构算子, 对离散点p (x, y, z) 的高程z先进行腐蚀, 操作如式 (4) 。

式 (4) 中, (xp, yp, zp) 是脚点p邻域内的点的坐标值, W是窗口边长, 设置W大于最大建筑物尺寸, hp_erosion是腐蚀过的点的高程值。通过对点的腐蚀运算, 点的高程值取为以本点为中心边长为W的窗口内点高程的最小值zi。将此时高程值hp_erosion与原始高程值z做差, 绝对值小于高程阈值的点作为种子点。因为腐蚀处理的窗口边长大于最大建筑物的边长, 所以建筑物和植被点能够全被滤除, 且由于窗口太大, 与邻近点高程相差稍大的地面点也被滤除, 此时的种子点均为地面点。

利用种子点拟合曲面, 如式 (5) 。

式 (5) 中a0至a5为系数, (x, y) 为点的坐标。拟合残差V由式 (6) 可得

式 (6) 中,

每个种子点对拟合曲面的影响是不同的, 权重的分配与腐蚀前后的高差变化成反比, 如式 (11) 。

权阵如式 (12) 。

根据最小二乘原理VTPV=min, 可解得曲面方程系数如式 (13) 。

此时得到的拟合曲面是粗糙的, 再计算网格内所有点的拟合高程值, 设置高差阈值h, 拟合点高程与实际高程之差小于h的点判断为地面点, 加入种子点集;大于h的点可能为地面点, 也可能为地物点, 再次拟合曲面后进行判别。根据高差设置点再次拟合的权重, 利用此时的种子点集重新拟合曲面, 再判断点, 迭代直到种子点集不再变化, 此时得到的种子点集即为地面点。算法流程图如图1。

算法步骤如下;

(1) 对所有离散点进行腐蚀运算, 得到腐蚀后的点的高程, 记为h_erosion;

(2) 计算所有点腐蚀后的高程与原始高程之差, 记为diff_erosion;

(3) 比较diff_erosion与特定阈值thre_erosion的大小, 若小于高差, 则将点选为种子点, 若大于高差, 则点为待判断点, 拟合曲面后进行判别;

(4) 判断是否存在新的点加入种子点集, 若存在则对种子点分配权重, 否则输出种子点集;

(5) 计算所有离散点到拟合面距离的残差V;

(6) 比较残差V与特定阈值V_thre的大小, 若小于阈值, 则将点加入种子点, 大于阈值则为待判断点, 再次拟合曲面后判别;回到步骤 (4) 。

3 实验过程与结果分析

以MATLAB为平台, 利用ISPRS提供的点云数据进行实验。ISPRS提供了多个数据样本点集[1,13], 每个点集的参考数据都事先人工做了分类, 标明了地面点或非地面点, Sample54共有点8 608个, 点密度0.173 2/m2, 点间距2.402 8/m, 范围为185×267 (m×m) 高程范围为228~294 m, 地面点比例为0.468 2, 地势坡度较大, 如图2所示。左图中颜色越深, 点的高程值越大, 其中右图为左图红色方框内的切面, 由图可见, 区域的高程变化剧烈。故利用Sample54实验, 可以很好的说明本文算法。

3.1 实验过程

实验处理的数据为.txt文件, 读入数据, 并显示原始点云插值图如图3 (a) 。对散点进行腐蚀, 选取处理后高差符合条件的点为种子点, 种子点内插为粗略DEM如图3 (b) 。利用最小二乘法拟合种子点为曲面, 判断散点与曲面的距离是否小于阈值, 将符合条件的点加入种子点集, 迭代得散点图如3 (c) 所示, 内插为DEM如图3 (d) 。

3.2 实验结果与分析

2003年, ISPRS挑选了8种滤波方法进行定量测试与定性分析[13]。这8种算法分别由Elmqvist, Sohn, Roggero, Brovelli, Wack, Axelsson, Sithole, Kraus提出 (表格1和图4中以提出者名字代指算法) 。Elmqvist提出了一种基于活动轮廓法的滤波方法, 其基本思想是利用高程的活动轮廓逼近局部最小值, 从而得到地面点。Sohn提出了渐进不规则三角网加密算法, 利用最小描述长度判断其余点是否为地面点, 重复进行向下加密及向上加密。Roggero提出了一种坡度滤波方法, 通过比较相邻点的坡度来滤波。Brovelli提出了线性曲线内插法滤波方法, 算法基于样条差值和区域生长。Wack改良了分层最小区域滤波方法, 利用拉普拉斯 (LOG) 操作元在由粗到细的网格内滤除高程急剧变化的点。Axelsson利用局部最低点作为种子点构建初始三角网, 渐进加密滤波。Sithole提出一种改进的坡度滤波算法, 根据坡度值图设置不同的斜率阈值。Kraus提出的分等级稳健内插法的滤波效果, 根据点到DEM的距离分配权值, 迭代加密。本文与这8种算法进行对比分析。

文章采用George Sithole等人[4]提出的评价指标。假定:A表示参考数据中地面点的总数, B表示参考数据中非地面点的总数, C表示被误分类为非地面点的地面点的总数, D表示被误分类为地面点的非地面点的总数。一类误差 (Ι类误差) 为地面点误判为非地面点的比例:C/A, 二类误差 (П类误差) 为非地面点误判为地面点的比例:D/B, 总误差为所有误判点的比例: (C+D) / (A+B) 。

本算法与原最小二乘算法和这8种算法误差结果相比得表1 (原最小二乘算法是本文编写) , 由表1可知, 本文算法一类误差为1.56%, 只有62个地面点被误判为非地面点, 说明大部分地面点都得到了正确判别, 一类误差最小的是Roggero提出的改进坡度滤波方法, 仅为1.01%, 最大的为Brovelli提出的线性曲线内插法滤波算法, 有近一半的地面点被误判为非地面点。相对于其他几类算法, 本文算法的一类误差较小;本文算法二类误差为4.63%, 有214个非地面点被误判为地面点, 大部分非地面点得到滤除, 二类误差最小的是Wack改良的分层最小区域滤波算法, 仅为0.48%, 几乎所有的非地面点都被滤除, 最大的是本文编写的最小二乘算法, 为9.30%, 其次是Roggero提出的改进坡度滤波算法, 为8.37%, ;本文算法总误差为3.21%, 比Axelsson算法的3.23%仅少0.02%, 但远小于其他几类算法, 尤其与原最小二乘算法的7.59%相比, 误差较小, 滤波效果较好。

将本算法得到的结果可视化, 并与原最小二乘算法和这8种算法相比较得到图4, 其中, 蓝色点表示一类误差点, 蓝色点越多, 被错分的地面点越多, 红色点表示二类误差点, 红色点越多, 被错分的非地面点越多, 浅灰色表示地面点, 深灰色表示非地面点。

图4 (a) 是由Elmqvist提出的基于活动轮廓法的滤波结果, 右上角蓝色点较多, 说明地面点被错分的数量较大;图4 (b) 为Sohn提出的渐进不规则三角网加密算法的滤波结果, 蓝色点主要集中在右侧, 错分的点不多;图4 (c) 由Roggero提出的改进坡度滤波方法所得, 蓝色点较少, 即错分的地面点较少, 但红色点较多, 有大量错分的非地面点;图4 (d) 为Brovelli提出的线性曲线内插法滤波所得, 蓝色点较多, 地面点被误判为非地面点的数目较大, 有近一半的地面点被错误分类;图4 (e) 是Wack改良的分层最小区域滤波方法所得, 几乎没有红色点, 说明二类误差较小, 非地面点得到很好的滤除;图4 (f) 是Axelsson提出的渐进不规则三角网加密算法的滤波结果, 蓝色点较少, 地面点得到很好的保留;图4 (g) 是由Sithole提出的改进的坡度滤波算法的滤波结果, 蓝色点较多, 地面点的分类误差较大;图4 (h) 是Kraus提出的分等级稳健内插法的滤波效果, 非地面点得到很好的滤除;图4 (i) 是本文利用MATLAB编写的最小二乘算法的滤波效果, 由图可知, 地面点没有得到很好的保留, 非地面点没有很好的滤除, 有大量点被误判, 且在右上角地势变化剧烈的区域错分点尤多, 可见原始最小二乘算法对地势变化剧烈的区域滤波效果不理想;图4 (j) 为本文滤波方法所得, 地面点得到很好的保留, 非地面点得以滤除, 图4 (i) 中的错分现象得到明显改进, 滤波效果较好。

4 结束语

利用形态学腐蚀算子进行种子点的选取, 能够得到真实地面点, 获取粗略的DEM, 根据点腐蚀前后的高程差分配权重, 使得首次拟合曲面具有先验知识, 从而适用于地形变化剧烈的区域。实验证明算法的滤波性能较好。算法的阈值选取与地形有关, 需要进一步的探讨。

摘要：最小二乘算法在“等权”的条件下进行, 不具有先验知识, 不适用地形复杂的区域。针对此问题, 提出了一种改进的基于最小二乘的激光雷达 (light detection and ranging, Li DAR) 数据滤波方法。首先对离散点高程进行腐蚀操作, 将得到的高程与原始高程相比, 选取小于阈值条件的点作为种子点。其次根据腐蚀前后的高程差值赋予点不同的权重, 利用最小二乘原理拟合这些种子点, 计算所有点的拟合值与实际高程的差;特定阈值比较并将小于阈值的点加入种子点, 迭代至地面点集不再变化。利用国际摄影测重与遥感学会 (international society for photogrammetry and remote sensing, ISPRS) 提供的数据进行实验, 滤波总误差小于经典算法。实验结果表明算法能够较好的滤除地物, 适用于地形复杂的区域。

滤波算法在测井数据预处理中的应用篇4

测井数据的准确性是保证测井解释结果可靠的前提。在测井数据处理的逐点计算中, 计算机是极严格地按照深度连续驻点取出各个采样点的测井数据来定量计算的。因此, 对测井曲线深度和幅度的准确性更有十分严格的要求。然而, 由于野外测井作业和测井环境的许多随机因素的影响, 即使采用数控测井及严格的技术措施, 同一口井各测井曲线之间的深度一致性也往往难以实现, 各测井数据曲线幅度也不可避免地要受到许多非地层的测量因素的影响。因此, 测井资料预处理是测井解释与数据处理的一项重要工作, 它是保证测井解释与数据处理结果精度的极重要前提。

在碳酸盐岩的声波测井中, 由于声波探头与井壁的随机碰撞干扰, 或在缝洞孔隙和裂缝发育的地层中声波经过多次反射折射, 使测出的声波曲线上出现许多毛刺干扰。显然, 用这些具有毛刺干扰的测井曲线作数字处理, 会给计算的地质参数带来很大的误差;统计起伏或毛刺干扰严重的曲线, 根本不能直接用作数字处理。因此, 必须设法把这些与地层性质无关的毛刺干扰滤掉, 只保留曲线上反映地层特性的有用成分。本文主要讨论的是直线滑动平均法、五点三次平滑法和局部线性估计光滑算法。

1 算法原理

1.1 直线滑动平均法

该方法主要根据对某点临近的采样点的波幅来对该点进行波幅修正, 从而达到对波形进行去噪的目的。一般取5个临近的点, 平均法的基本计算公式如下:

其中:x为采样数据;y为平滑处理后的数据;m为数据点的个数;N为所取的平均点个数;h为加权平均因子

对于简单平均法, 即:

直线滑动平均法就是利用最小二乘法对离散数据进行线形平滑的方法。5点直线滑动平均法的计算公式为:

公式中:i=3, 4, …m-2

1.2 五点三次平滑法

五点三次平滑法是利用最小二乘法原理对离散数据进行三次最小二乘多项式平滑的方法, 五点三次平滑法计算公式为:

公式中:i=3, 4, …m-2

1.3 基于局部线性估计的光滑处理算法

这种光滑处理算法是一种基于临近点的局部线性加权估计算法, 光滑处理后的数据是建立在对原有的每一个数据点的局部线性回归上的。Lowess (局部加权散点平滑法, 采用线性最小二乘法和一阶多项式拟合得到的数据进行替换) 由cleveland在1985年提出, 是加权回归中最精巧的形式, 方法类似于滑动平均技术, 是在指定的窗口之内, 每一点的数值都用窗口内临近的数据进行加权回归得到的, 回归方程可用线性的或者二次的。如果在指定的窗口宽度之内, 拟进行平滑的数据点两侧的进行平滑的数据点是相等的, 则为对称Lowess, 如果两侧数据点不等, 则为非对称Lowess。

利用Lowess方法进行数据平滑处理的重点参数在于窗口宽度的选择, 窗口宽度过大, 将使得光滑描点涵盖的历史数据过多, 降低最新数据信息对平滑值的影响。反之, 过窄的窗口宽度使“平滑”后的数据并不平滑。本实验中, 选择的是窗口大小为5。简而言之就是在最小二乘法的基础上增加一次滑动窗口, 每个窗口根据实际情况确定大小, 就是在计算每个光滑点的同时, 用它附近所有窗口点的大小值。因此每个点都有一个权值, 离当前点越远它的权值越小, 权值根据欧式距离的计算公式为:

其中δmax为离当前点最远的点到当前点的距离

2 编程实现

MATLAB可以对通过编制的M文件直接进行运行和调试, 由于M文件是解释性的程序语言, 所以M文件无论从形式、结构、语法规则等方面都比一般计算机高级语简单、易写、易读得多, 而且程序容易调试, 人机交互性强。MATLAB的一系列工具箱所包含的全部函数都是M文件编写的。工具箱可以灵活编程实现所需的要求;并且MATLAB语言提供了一套强大的图形程序, 为计算过程和结果的可视化提供了极佳的手段。

因此, 我们就是采用MATLAB下的lowess函数,

function c=lowess (x, y, span, method, robust, iter, ws)

根据上节讨论, 一般来说, Lowess方法包括以下步骤:

Step1, 计算指定窗口内各个数据点的初始权重, 权重函数一般表达为数值之间欧氏距离比值的立方函数;

Step2, 利用初始权重进行回归估计, 利用估计式的残差定义稳健的权函数, 计算新的权重;

Step3, 利用新的权重重复步骤2, 不停地修正权函数, 第N步收敛后可根据多项式和权重得到任意点的光滑值。

利用MATLAB可以很方便实现这个过程, 而其余两种方法, 我们只需将算法编制后制作成相应的M文件, 方便调用, 最后我们根据PLOT () 函数对结果绘图。

3 算法验证

由上述图片可知, 五点平滑法应用比较简单, 但是平滑后的数据没有保持原来曲线特征;而五点三次法保留了曲线的特征, 但是又不够光滑;lowess方法在保持曲线特征的同时又光滑, 能够抵御异常点引起的偏差, 但是缺点在于计算复杂, 而且在多于2 000个数据点以上这些方法都不合适, 就要使用小波等信号处理工具。在实际工程应用中, 要根据现实状况选择比较合适的方法。

参考文献

[1]孙苗钟.基于MATLAB的振动信号平滑处理方法[J].电子测量技术, 2007 (6) .

[2]王济, 胡晓.MATLAB在振动信号处理中的应用[M].北京:中国水利水电出版社, 2006.

[3]雍世和.测井数据处理与综合解释[M].东营:石油大学出版社, 2002.

[4]刘伍丰, 何旭春.数据平滑处理算法的编程[J].软件时空, 2007 (3) .

[5]代广进, 侯正信.小波域信号去噪算法[J].电子测量技术, 2005 (6) .

数据滤波篇5

实验的目的是把传统的基于线性滤波器的FBG解调系统受环境温度变化带来的影响降低,从而提高系统的稳定性与可靠性。对系统改进后实现基于GPRS网络的数据传输系统应用[1,7,8]。使基于线性滤波器的传感解调系统达到远程、智能、高效的目的。实验中采用GPRS网络对解调系统所测物理量进行远程实时在线监测,从而可以快捷、方便地得出FBG传感器所处环境的温度值。以实现解调系统的远程实时智能测量。

1 基于线性滤波器的FBG传感解调系统

宽带光源发出的光经过隔离器后经过耦合器1到达FBG,FBG产生的反射光波长会随外界温度的变化而改变。反射光再经耦合器1后到达耦合器2,耦合器2将FBG的反射光等功率分成两路光,S1作为参考光进行直接探测,S2经过线性滤波器作为探测光,对采集到的光信号进行光电转换变为电信号,再经过AD采集模块采样[9]。最后经过有线传送到计算机。实现数据的显示与存储。解调原理图如图1所示。

利用光谱分析仪测量线性滤波器的透过率曲线,光谱图如图2所示,在1 540.0～1 541.5 nm的波长范围内,光功率与波长近似成线性关系。本实验是利用此波段的光信号来进行解调实验,实验中采用中心反射波长为1 540.0 nm的光纤Bragg光栅。

两路光信号S1和S2通过光电探测器变为相应的电压信号U1和U2,电压信号值U1和U2的大小反映了S1和S2两路光信号光功率的大小,利用数据采集器对电压信号U1和U2进行采样,将两路光电压值进行比较得B=U1/U2,如下式

式中,S(λ)为光源功率,T1(λ)为耦合器1的透过率,T2(λ)为耦合器2的透过率,R(λ)为FBG的反射率,D(λ)为线性滤波器的透过率。由于S(λ)、T1(λ)、R(λ)、T2(λ)为常数,则式(1)可简化为

由式(2)可知,比值的变化只与线性滤波器的透过率D(λ)有关,当FBG所处环境温度不变时,FBG反射波长一定,此时线性滤波器的透过率理论上应为一常数,这时比值应该为一个常数。

而实际上,当FBG所处环境温度不变时,线性滤波器所处温度变化时,比值发生了一定的变化,这时观察两路光S1和S2在线性滤波器所处不同温度下电压值的变化情况及比值变化情况,如表1所示。

从表1中可知,比值随着线性滤波器所处环境温度的变化而发生了一定的变化。因此线性滤波器所处的环境温度将对解调系统的稳定性和准确性有很大影响。

2 加温控装置后解调系统的稳定性测量

由表1可知,线性滤波器的透过率受其所处环境温度影响较大,因此实验设计了一套恒温系统,将线性滤波器置于恒温系统中,这样线性滤波器的透过率就不会随着外界环境温度的变化而变化。实验将恒温系统的温度值控制在40℃±0.5℃来减小线性滤波器受环境温度变化带来的影响。从而使得该解调系统的测量值更加精确,通过改进使系统的性能得到极大提高,恒温系统原理框图如图3所示。

恒温系统采用DS18B20温度传感器来测量线性滤波器所处环境的温度值,当线性滤波器所处环境温度低于预设的温度值时,由ARM7控制PMOS管导通实现对加热布供电,对恒温系统进行加热,使线性滤波器所处环境温度基本恒定在预设的温度值。当线性滤波器所处环境温度高于预设的温度值时,由于没有设计降温系统,此时温控系统将失去作用。因此,为了避免出现环境温度大于恒温系统预设温度值,应该将预设温度值设置为大于环境温度,试验中选择了40℃。对线性滤波器加上恒温系统后当FBG所处环境温度不变时,每隔半小时对两路光信号S1和S2的电压值进行一次数据采集,采集数据如表2所示。

从表2中可以看出,在加入温控装置后,两路光信号的比值基本恒定在一个常数的状态下,这对于解调系统的稳定性有了极大地提高。此时整套解调系统只有FBG的反射光波长受环境温度变化而变化,文中实验正是利用此变化来测量环境温度。

3 基于GPRS的现场数据采集系统的设计

基于GPRS的远程监控系统是集现代无线通信技术、信号采集技术以及计算机网络技术为一体的现代监控系统,已经成为现代监控领域发展的新趋势。该系统通过GPRS网络与Internet相结合,实现了计算机与远程ARM监控终端的通信。整个系统由监控中心上位机软件、无线通信链路和监控终端(FBG传感解调设备)三部分组成[1]。监控中心上位机软件完成与远程FBG传感解调设备连接及对设备的控制,对远程FBG传感解调设备传回来的数据进行分析、显示、存储等操作。远程FBG传感解调设备主要用来采集现场温度、图像或视频信息,向监控中心上位机软件传送相关信息。GPRS网络起桥梁作用,完成与Internet网络的信息交互[3]。数据采集系统原理框图如图4所示。

基于GPRS网络的现场数据采集系统由基于线性滤波器的FBG传感解调设备采集现场的两路光信号的能量值进行比值,由SIMCOM公司的一款三频段GSM/GPRS模块SIM300模块通过GPRS网络将比值信息、图片或图像信息传送到监控中心上位机软件,SIM300内部集成了TCP/IP协议栈、具有双音频通道、并具有低功耗模式等优点,因此特别适合于开发数据传输设备。上位机软件根据SIM300通过GPRS网络传回的比值信息与实验阶段存储的标准比值与温度关系数据库信息比较,最后得出现场的温度值。当温度值超过规定的值时,监控中心上位机软件会下发报警信息,并通过GPRS网络传送给数据采集系统的中央处理器(ARM7-LPC2364)[2,6],中央处理器控制报警器进行报警,并通知摄像头采集现场的图片或视频信息,通过GPRS网络传送到监控中心的上位机进行显示,监控中心的工作人员可以根据传回的图片或视频信息进行相应的控制。并可通过

模块的音频功能与设备现场工作人员进行语音通话。最终完成对设备的远程、实时、智能化控制。

远程监控中心上位机部分介绍:FBG传感解调系统与监控中心主站的上位机软件以TCP/IP方式进行连接通信,上位机软件实时显示监控环境内的温度值,显示监控终端设备传回的图片信息或视频信息[5]。上位机提供存储信息功能,此外还提供报警功能。当探测点温度超过设定值时监测系统的CPU输出报警信号,也可通过上位机向CPU通过GPRS网络传送新的报警信号、传送摄像或拍照命令,达到永久在线、智能监测的目的。监控中心上位机软件原理图如图5所示。

监控中心上位机软件采用美国NI公司推出的LabWindows/CVI6.0软件编写,该软件具有使用简单、灵活、易于上手等优点,特别适合用于工业控制、通信测试、设计验证等环境。该上位机软件通过TCP/IP方式与Internet进行连接,对现场数据采集系统进行控制与状态显示,同时具有存储操作日志、绘制温度曲线、显示及存储视频或图片信息等功能。

4 总结

由于线性滤波器本身的透过率受所处环境温度的影响而产生一定的变化,进而引起两路光信号S1与S2光功率比值的不稳定。实验中,通过引入温控装置使得线性滤波器处在恒温条件下来消除这种不稳定对解调系统的准确性带来的影响,从而极大地提高了系统的稳定性。在此基础上设计了解调系统的远程、实时、智能化监测系统,采用GPRS无线网络与Internet相结合的方式作为连接远程监测终端和监控中心上位机的传输网络,取代了传统的有线传输,实现了远程实时智能化监测和快速、精确解调的目的。系统可用于对特殊环境下的温度、压力等应变量的测量,具有一定的实用价值。

参考文献

[1]廖胜,胡文齐.基于ARM和GPRS远程监控系统的研究[D].北京:北京邮电大学,2008.

[2]周立功,王祖麟,陈明计,等.ARM嵌入式系统基础教程[M].北京:北京航空航天大学出版社,2008.

[3]韩基荣.基于GPRS远程监控系统的设计[D].哈尔滨:哈尔滨工程大学,2009.

[4]卢挺.基于GPRS的路灯远程监控系统的设计[D].大连:大连理工大学,2009.

[5]张毅刚,乔立岩.虚拟仪器软件开发环境:LabWindows/CVI6.0编程指南[M].北京:机械工业出版社,2002.

[6]谭浩强.C程序设计[M].3版.北京:清华大学出版社,2005.

[7]艾谦.光纤光栅在电力系统测温中的应用研究[D].武汉:武汉理工大学,2006.

[8]任建新,熊亮,张鹏.基于GPRS的油井远程监控系统设计[J].测控技术,2010,29(8):98-101.

数据滤波篇6

地球介质中的许多岩石都呈现各向异性的特征,并且地层的各向异性系数和快横波方位能够用来指示地层的裂隙和地应力信息,有着重要的应用。测量的方法原理是根据各向异性地层中偶极声源激发的快慢弯曲波具有波形相似的特征,对所有接收阵列的快慢波形进行匹配并构造目标函数,利用反演方法计算地层的各向异性系数和快横波方位[1]。但偶极阵列声波数据往往会被噪声干扰而使我们计算的各向异性产生错误的结果,因此对偶极声波数据的滤波处理是非常重要的,本文主要介绍一种新的滤波方法,与传统时域滤波相比它是对频率波数域进行滤波,极大的提高了偶极声波数据的质量。

1 频率波数域相干滤波

1.1 方法原理

对于阵列声波数据中的一些低信噪比的信号,处理通常是比较困难的。Schmitt首次针对井中声场采用频率波数域的分析方法[2],我们利用阵列数据的相干性实现了一种新的阵列波形处理技术。我们首先将原始数据转换到频率波数域,然后计算一个相干函数,并将计算得到的相干函数与频率波数域中数据进行卷积,最后对需要的相干信号进行逆变换就可以获得相干滤波后的阵列数据,可由如下(1)表示[3]。

$X_{c f i l} (k, w) = X (k, w) C o h (k, w) (1)$

式(1)中Xcfil(k,w)是频率波数域相干滤波的波形数据,X(k,w)是对一个二维的声波波形时间序列做二维傅里叶变换得到的频率波数域的波形数据,如式(2)所示[3]。

$X (k, w) = \iint \begin{matrix} X (x, t) \end{matrix} e^{- i (k z - w t)} (2)$

而Coh(k,w)是我们构造的相干函数,如偶极子中的弯曲波,我们可以构造一个相干函数如式(3)[3]。

$C o h (k, w) = \frac{∥ \sum_{n = 1}^{Ν} X_{n}^{*} (w) e^{i k (n - 1) d} ∥}{\sqrt{Ν \sum_{n = 1}^{Ν} X_{n}^{*} (w) X_{n} (w)}} (3)$

式(3)中Xn(w)是对从阵列数据中各道进行傅里叶变换获得的谱阵列数据进行的近似,*号代表对数据进行复共轭变换。式(3)的本质是对频率波数域平面阵列数据进行相似性叠加。对频率波数域的相干滤波图如图1所示,其中纵轴为波数,横轴为频率,图中亮点是我们需要的有用信号。

1.2 处理实例

在地层进行声波测井时,沿着地层传播的声波信号常常会被别的不同路径的声波信号干扰。例如,在随钻测井中,声波信号可能会沿着仪器本身传播并对地层信号造成严重的干扰。在套管井中,当套管与水泥胶结不好时,沿着套管传播的套管波可能变得很大。在这些条件下,都可以利用频率波数域相干滤波的方法来获得地层信号的有用信息。我们采用实际阵列数据接收到的偶极波形数据进行频率波数域滤波如下图所示。图2是原始波列,可以看出阵列波形之间的差异比较大

我们对这组数据进行频率波数域相干滤波后对其进行二维傅里叶逆变换得到滤波后的数据如图2所示。从中可以看出进行频率波数域相干滤波后的数据已经得到了很好的改善。

2 各向异性的计算

传统的滤波方法就是对偶极横波阵列数据做带通滤波,效果不是很理想,本文基于频率波数域滤波的方法对偶极横波阵列数据做相干滤波,得到比以往带通滤波更高质量的声波数据,并将滤波后的声波数据用于各向异性的计算。用阵列偶极横波数据计算地层的各向异性通常是一个二维反演求全局极小值的过程,如今我们经常用的计算地层各向异性的软件是由美国Baker Hughes公司研发的eXpress软件,它采用的是模拟退火的方法求二维平面内的全局极小,模拟退火的方法虽然能够很好的求出地层各向异性大小和快波方位,但它花费的时间过长,而别的处理软件虽然时间花费较少,但处理结果的精度很低。本文基于极快速模拟退火的理论,用频率波数域滤波后的偶极阵列声波数据进行各向异性的计算,不但达到了很高的精度而且使计算时间大大缩短,具有重要的应用意义,极快速模拟退火主要是取一维的Cauchy分布,用不同的初始温度来反映不同参数对退火时间敏感性的不同,并采用改进的算法和给出冷却进度表来对模拟退火算法进一步提速[4],虽然处理速度较快但各向异性的精度有待提高,因此我们采用频率波数域滤波后的偶极阵列声波数据进行各向异性的计算,结果如图4所示。

图4显示的是一个实际的油田数据计算的各向异性,下部各向异性较大的地方裂隙比较发育,开发后产气,其中右边两道中黑色的线是采用eXpress处理软件计算出的各向异性的大小和快波方位曲线,而红色的曲线是我们将相干滤波后的阵列声波数据用极快速模拟退火的方法计算出的各向异性大小和快波方位曲线,可以看出两者有很高的相似性,而我们处理各向异性所花费的时间非常短,比eXpress软件快了好几倍。具有重要的应用意义。

3 结论

通过将频率波数域滤波后的偶极声波数据用极快速模拟退火的方法计算结果可以得出如下结论:

(1)对偶极阵列声波数据进行频率波数域相干滤波可以大大提高波形的相似性和质量,对我们后面的进一步处理有重要的意义。

(2)对进行频率波数域滤波后的数据采用极快速模拟退火的方法来计算地层各向异性的大小和快波方位不但具有很高的精度而且大大减少了计算时间。

参考文献

[1]苏远大,乔文孝,陈雪莲,等.正交偶极声波测井资料在评价地层各向异性中的应用.石油物探,2005;44(4):409—412

[2] Schmitt D P,Bouchon M.Full wave acoustic logging:Synthetic mi-cro-seismogram and frequency-wavenumber analysis.Geophysics,1985;50:1756—1778

[3] Tang X M,Land S T X.Method for coherence-filtering of acoustic ar-ray signal.U.S0287830,2006-12—21

数据滤波篇7

关键词：鲁棒自适应卡尔曼滤波,平稳时间序列,陀螺漂移,趋势项,建模

惯导系统作为精确制导武器上核心系统,其精度直接影响到武器的打击精度。要提高惯导系统的精度,目前最有效的方法是对误差进行标定并加以补偿。当飞行器竖立在发射台上时会受到各种干扰因素的影响,比如阵风干扰,补偿这种干扰现阶段比较有效的作法是,首先建立起阵风干扰下平台漂移误差模型,然后选用合适的滤波算法对数据进行处理,以期达到较好的效果。但是在实际计算中,普通的卡尔曼滤波在数据处理方面是基于小误差情况,也就是平台在工作过程中没有受到大的影响,这种条件并不适用于本文所研究的情况,所以必须要选用一种收敛效果好的滤波算法,鲁棒自适应卡尔曼滤波就满足这种要求。它能够在线估计干扰误差,并且在建模时就考虑了模型本身所带来的误差并进行补偿,使数据处理达到较理想的效果。

近年来,时序分析方法被越来越多地应用到惯导系统漂移误差建模中, 本文首先利用ARMA(2,1)模型建立起阵风干扰下平台漂移模型,然后利用鲁棒自适应卡尔曼滤波对数据进行处理取得了比较好的效果。

1 数据分析

1.1 数据采集

香农定理:采样频率必须大于信号最高频率的两倍。这一定理是信号采集的依据,本文所采样的间隔为0.1 s,采样时长为60 s,这样就得到了600个数据,如图1所示。

1.2 数据平稳性检验

陀螺在生产安装过程中,由于各种因素的影响会产生安装误差,这就导致了陀螺漂移数据表现出一定的趋向性,所以需要对其进行平稳性检验。

常用的检验方法有参数和非参数分析法两类,本文采用逆序检验法属于后者,检验结果表明,陀螺仪的漂移数据存在趋势项,不满足平稳性要求。

2 陀螺漂移数据建模

2.1 趋势项的提取

对于趋势项:

$\begin{array}{l} D_{1} (t) = β_{0} + β_{1} t + β_{2} t^{2} + β_{3} t^{3} + β_{4} t^{- 1} + \\ β_{5} t^{- 2} + β_{6} t^{\frac{1}{2}} + β_{7} t^{\frac{- 1}{2}} (1) \end{array}$

在实际测试中的趋势不一定如此复杂,通常取两三项即可,采用逐步回归法进行因子舍取,

令t=1,2,…,N,可得矩阵方程

x=βA+ε;

$x = [\begin{matrix} x_{1} & x_{2} & \dots & x_{Ν} \end{matrix}]^{Τ}$ ;

$β = [\begin{matrix} β_{0} & β_{1} & \dots & β_{7} \end{matrix}]^{Τ}$ ;

$\begin{array}{l} A = \\ [\begin{matrix} 1 & 1 + \frac{1}{Ν} & (1 + \frac{1}{Ν})^{2} & \dots & (1 + \frac{1}{Ν})^{- 1 / 2} \\ 1 & 1 + \frac{2}{Ν} & (1 + \frac{2}{Ν})^{2} & \dots & (1 + \frac{2}{Ν})^{- 1 / 2} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 1 & 1 + \frac{Ν}{Ν} & (1 + \frac{Ν}{Ν})^{2} & \dots & (1 + \frac{Ν}{Ν})^{- 1 / 2} \end{matrix}] 。 \end{array}$

参数β的最小二乘估计式为:

$β = (A^{Τ} A)^{- 1} A^{Τ} x$

结果表明,漂移信号含有一阶线性趋势项

2.2 正态性检验

正态性检验就是检验时序的3阶矩(偏态系数ζ)和4阶矩(峰态系数v)是否满足正态随机变量的概率密度函数,即ζ=0,v=3,经检验陀螺漂移数据满足正态性要求。

2.3 平稳时间序列建模

本文采用ARMA(2,1)模型对漂移数据建模其算法如下:

$[\begin{array}{l} φ_{1} \\ φ_{2} \\ ⋮ \\ φ_{n} \end{array}] = [\begin{array}{l} θ_{1} \\ θ_{2} \\ ⋮ \\ θ_{n} \end{array}] + [\begin{matrix} 1 & 0 & \dots & 0 \\ - θ_{1} & 1 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ - θ_{n - 1} & - θ_{n - 2} & \dots & 1 \end{matrix}] [\begin{array}{l} Ι_{1} \\ Ι_{2} \\ ⋮ \\ Ι_{n} \end{array}] (4)$

$[\begin{array}{l} Ι_{n + 1} \\ Ι_{n + 2} \\ ⋮ \\ Ι_{n + m} \end{array}] = [\begin{matrix} Ι_{n} & Ι_{n - 1} & Ι_{n - 2} & \dots & Ι_{n - m + 1} \\ Ι_{n + 1} & Ι_{n} & Ι_{n - 1} & \dots & Ι_{n - m + 2} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ Ι_{n + m - 1} & Ι_{n + m - 2} & Ι_{n + m - 3} & \dots & Ι_{n} \end{matrix}] [\begin{array}{l} θ_{1} \\ θ_{2} \\ ⋮ \\ θ_{m} \end{array}] (5)$

经过matlab计算,得到如下的模型:

$x_{(t)} = 0.133 9 \bar{x}_{(t - 1)} - 0.681 7 \bar{x}_{(t - 2)} + a_{(t)} + 0.802 7 a_{(t - 1)} (6)$

3 鲁棒自适应卡尔曼滤波

3.1 滤波方程建立

基于上述模型写出如下状态空间模型:

${\begin{cases} X_{(t + 1)} = Φ X_{(t)} + Γ W_{(t)} \\ Ζ_{(t)} = Η X_{(t)} + V_{(t)} \end{cases} (7)$

式(7)中,W(k),V(k)的统计特性为:E(W(t))=0,E(V(t))=0;E(V(t)VT(t))=Qtδtj,E(W(t)WT(t))=Rtδtj,E(V(t)WT(t))=0。δtj是荻拉克利函数。

$δ_{t j} = {\begin{matrix} 0 ‚ & t \neq j ； \\ 1 ‚ & t = j 。 \end{matrix}$

系统的状态为: $X_{(t)} = [\overset{⌢}{x}_{(t)} \overset{⌢}{x}_{(t - 1)}]^{Τ}$ ,过程噪声为:W(t)=[a(t)a(t-1)]T。

$\begin{array}{l} 式中的系数为 ∶ Φ = [\begin{matrix} 0.133 9 & - 0.681 7 \\ 1 & 0 \end{matrix}] ‚ Γ = \\ [\begin{matrix} 1 & 0.802 7 \\ 0 & 1 \end{matrix}] 。 \end{array}$

系统的输出为Z(t)=x(t),则输出方程系数H=[1 0]。

在这个模型中,把需要辨识的漂移角设为一个状态变量,然后通过参数辨识的方法把它估计出来。

3.2 鲁棒自适应卡尔曼滤波

考虑上述的线性离散随机系统,假设其是带模型误差系统,即真实系统是如式(8)。

${\begin{cases} X (t + 1) = (Φ + Δ Φ) X (t) + Γ w (t) \\ Ζ (t) = (Η + Δ Η) X (t) + v (t) \end{cases} (8)$

式(8)中状态转移矩阵误差ΔΦ和观测矩阵误差ΔH是未知的,但Φ和H是已知的。因此,既使噪声统计是已知的,对带模型误差系统进行常规卡尔曼滤波也会使滤波器性能变坏。为了补偿模型误差,将真实系统写为:

${\begin{cases} X (t + 1) = Φ X (t) + ε (t) \\ Ζ (t) = Η X (t) + η (t) \end{cases} (9)$

式(9)中定义虚拟模型噪声ε(t)和虚拟观测噪声η(t)为

${\begin{cases} ε (t) = Δ Φ X (t) + Γ (t) w (t) \\ η (t) = Δ Η X (t) + v (t) \end{cases} (10)$

虚拟噪声ε(t)补偿了状态模型误差ΔΦX(t),虚拟噪声η(t)补偿了观测方程误差ΔHX(t)。通常模型误差ΔΦ,ΔH相对于Φ,H而言是较小的,因此可近似假设虚拟噪声ε(t),η(t)是带未知时变噪声统计的相互独立的白噪声。

$\begin{array}{l} E ε (t) = q (t), cov [ε (t), ε (j)] = Q (t) ； \\ E η (t) = r (t), cov [η (t), η (j)] = R (t) 。 \end{array}$

对真实系统的滤波问题转化为带未知时变噪声系统的自适应卡尔曼滤波问题,得到如下鲁棒自适应卡尔曼滤波公式:

给定初值: $\overset{⌢}{X} (0 | 0) ‚ Ρ (0 | 0) ‚ r (0) ‚ Q (0) ‚ \overset{⌢}{q} (0) ‚ b 。$

dt=(1-b)/(1-bt+1) (11)

$\begin{array}{l} \overset{⌢}{X} (t + 1 | t + 1) = \overset{⌢}{X} (t + 1 | t) + Κ (t + 1) \times \\ ε (t + 1) (12) \end{array}$

$\begin{array}{l} \overset{⌢}{X} (t + 1 | t) = Φ \overset{⌢}{X} (t | t) + Γ \overset{⌢}{q} (t) (13) \\ ε (t + 1) = Ζ (t + 1) - Η \overset{⌢}{X} (t + 1 | t) - \overset{⌢}{r} (t) (14) \\ Κ (t + 1) = Ρ (t + 1 | t) Η^{Τ} \\ [Η Ρ (t + 1 | t) Η^{Τ} + \overset{⌢}{R} (t)]^{- 1} (15) \\ Κ (t + 1) = Ρ (t + 1 | t) Η^{Τ} \\ [Η Ρ (t + 1 | t) Η^{Τ} + \overset{⌢}{R} (t)]^{- 1} (16) \\ Ρ (t + 1 | t) = Φ Ρ (t | t) Φ^{Τ} + Γ \overset{⌢}{Q} (t) Γ^{Τ} (17) \\ Ρ (t + 1 | t + 1) = [Ι - Κ (t + 1) Η] Ρ (t + 1 | t) (18) \\ \overset{⌢}{q} (t + 1) = (1 - d_{t}) \overset{⌢}{q} (t) + \\ d_{t} [\overset{⌢}{X} (t + 1 | t + 1) - Φ \overset{⌢}{X} (t | t)] (19) \\ \overset{⌢}{Q} (t + 1) = (1 - d_{t}) \overset{⌢}{Q} (t) + \\ d_{t} [Κ (t + 1) ε (t + 1) ε^{Τ} (t + 1) Κ^{Τ} (t + 1) + \\ Ρ (t + 1 | t + 1) - Φ Ρ (t | t) Φ^{Τ}] (20) \\ \overset{⌢}{r} (t + 1) = (1 - d_{t}) \overset{⌢}{r} (t) + \\ d_{t} [Ζ (t + 1) - Η \overset{⌢}{X} (t + 1 | t)] (21) \\ \overset{⌢}{R} (t + 1) = (1 - d_{t}) \overset{⌢}{R} (t) + d_{t} [ε (t + 1) \\ ε^{Τ} (t + 1) - Η Ρ (t + 1 | t) Η^{Τ}] (22) \end{array}$

其中,T为转置,I为单位矩阵,b是遗忘因子(通常取值为0.95≤b≤0.995)。以上就构成了自适应卡尔曼滤波器,交替运用式(11)～式(22)就可计算出状态和噪声统计特性的估计。

3.3 数据处理结果

将处理后的平稳时间序列作为系统输出输入到上述滤波器中,交替运用公式可得出如图3所示的数据处理结果。

从图3中可以看出,图3所得到数据明显比图2平滑了许多,为了更好F 说明问题,表1给出了数据处理前后的方差值,方差大时说明数据较分散,方差小时说明数据较集中。

从表1中可以看出,滤波后数据方差比滤波小了一个数量级,说明滤波效果较好,所得的数据可信,鲁棒自适应卡尔曼滤波在陀螺大干扰情况下能够取得很好的效果。

4 结论

本文首先利用ARMA(2,1)模型建立起阵风干扰下平台漂移模型,并且将其用状态空间的形式描述,通过模型适用性检验得知这个模型能够满足数据处理要求。在对数据处理时,考虑到普通卡尔曼滤波对大干扰情况下数据处理效果不好而且需要预知系统干扰的方差等信息,本文选用的鲁棒自适应卡尔螺滤波可以在线估计系统噪声的方差和协法差,而且在对大干扰数据处理的情况下取得了较好的效果,实际实验数据表明,经过滤波后的数据的方差可以比原始数据方差小一个数量级,满足了标定系数计算精度要求。

参考文献

[1]杨叔子,吴雅,等.时间序列分析的工程应用.武汉:华中理工大学出版社,1992

[2]武丽花,凌林本.三浮陀螺仪漂移模型的建立及MATLAB仿真.中国惯性技术学报,2004;12(6):75—78

[3]李纪光,刘锡敬.MATLAB在液浮陀螺随机漂移建模中的应用.中国舰船研究,2008;4(3):50—53

[4]李璐,孙纯祥.三浮陀螺随机漂移数据建模.压电与声光,2011;6(3):378—381

数据滤波篇8

关键词：小波分析,Kalman滤波,GPS,数据去噪

1 基于小波分析的GPS数据去噪方法

1.1 小波的种类

1)Haar小波。Haar函数是小波分析中最早用到的一个具有紧支撑的正交小波函数,也是最简单的一个小波函数,它是支撑在t∈[0,1]范围内的单个矩形波。Haar函数定义如下

2)Daubechie小波(dbN)。在dbN中,N代表次序。这族小波包括Haar小波,写作dbl,最简单的小波可以想象是最早期的。它们是紧支撑小波并且在给定的支持宽度拥有最高数目的消失矩。他们具有正交,双正交和紧致性。例如dbl或者Haar小波,db4小波,它们的消失矩为N。小波函数φ(t)和尺度函数Υ(t)中的支撑取为2N-1,φ(t)的消失矩为N。除N=1外,dbN不具有对称性(即非线性相位)且没有明确的表达式,但转换函数h的平方模是很明确的。

3)Symlet小波(symN)。Symlet小波函数是Daubechies提出的近似对称的小波函数,它是对db函数的一种改进。Symlets小波系通常表示为SymN。 Daubechies希望在保持小波足够简单的基础上修改他的小波对称性。方法是以z变换的形式|m0(ω)|2,设它为W,因此,可以用不同的方法分解 $W ‚ W (z) = U (z) \bar{U (1 / z)}, W$ 不等于1的根组成一对。如果选择所有的根都严格小于1的U,就得到dbN小波,U滤波器是最小相位滤波器;如果选择其他的,就得到更为对称的滤波器,这就是Symlets滤波器。Symlets还有与其他dbN相似的属性。他们在给定的支持宽度下都拥有很高的消失矩。N可以上2,3……。他们具有正交性,双正交性和紧致性。消失矩的长度变为N。

4)Coiflet小波(coifN)。Coiflet小波函数是由Daubechies构造的小波函数,它包括coifN(N=1,2,3,4,5)一系列。coifN具有比dbN更好的对称性。从支撑长度的角度看,coifN具有和db3N及sym3N相同的支撑长度;从消失矩的数目来看,coifN具有和db2N及sym2N相同的消失矩数目(见表1)。

1.2 GPS信号结构和噪声特性

1)GPS信号结构。

每个GPS卫星均发射可区分的信号,包括载波信号、P码(Y码)、C/A码和数据码(D码)等多种信号分量,所有这些信号分量都是在同一个基本频率f0=10.23 MHz的基础上产生。其中的P码和C/A码统称为测距码,载波信号包括L1载波和L2载波,频率为1 575.42 MHz和1 227.60 MHz。L1信号是由两个正交载波分量组成,同相分量调制有10.23 MHz钟速的精密伪随机噪声信号(P码),90°相移分量上调制有1.023 MHz的粗/捕获伪随机噪声信号(C/A码),L2信号上仅调制有P码,L1,L2信号均还连续调制SOBPS的导航电文,所有信号使用同一频率相干导出,L1和L2的载波均是正弦波,其卫星发射的信号可表示为 L1:SL1=APP(t)D(t)cos(2πf1t+

φ1)+AcG(t)sin(2πf1t+φ1), (2)

L2:SL2=BPP(t)D(t)cos(2πf2t+φ2). (3)式中:AP、BP和Ac分别为P码和C/A码的振幅; f1、f2为L1和L2的载波频率; φ1、φ2为L1和L2载波信号的初相; P(t)、G(t)、D(t) 分别为P码、C/A码和导航电文。

2)GPS信号的噪声特性。

GPS技术与伪距码是分不开的,GPS接收机利用伪距码接收伪距信号,辨别GPS卫星。精密GPS相位测量中利用伪距码解算相位的模糊度。伪距测量中有很多误差源,如电离层折射、对流层折射、多路径效应、钟差影响和接收机噪声等,关于对流层折射和电离层折射已有许多技术或模型来消除和改正,钟差影响也可以用一定模型来较好地消除,但多路径效应不同,接收机不同或地方不同所受多路径影响都不一样。多路径效应在GPS数据处理中不容忽视,它与接收机噪声是GPS精密定位的重要误差源。由于GPS接收机观测噪声和多路径延迟误差是GPS相位观测值中的一种噪声信号,从小波检测信号的角度分析,可以认为它没有幅值上的突变,是平稳信号。在GPS信号中,信号的有用频率可理解为信号发射频率、数据采样频率和数据解算频率等,至于哪种信号才是我们需要的有用信号,取决于我们准备从数据中提取何种信息,那么该信息的频率范围即为我们所感兴趣的,即有用的频率,该频率范围以外的信息都可视为无用信息或称为噪声。从纯频率的角度看,各种观测噪声和多路径效应的信号都有一定的频率范围,在GPS观测数据序列中的有用信号和噪声的时频特性通常是不一样的。有用信号在时域和频域上是局部化的,表现为低频特性或平稳的信号,而观测噪声和多路径效应在时频空间中的分布是全局性的,在整个观测时域内处处存在,频域上表现为高频信号。

1.3 小波去噪的方法

1)小波分解与重构法。

小波分解与重构法去噪的基本思想类似于FFT方法中多通道的滤波去噪,所不同的是小波方法是将信号变换到小波域,而FFT是将信号变换到频域。对于确定性噪声在特定的情况下,由先验知识可事先确定噪声的频率范围,且信号和噪声的频带不相互重叠时,用该方法就特别有效。应用小波分解与重构法的具体步骤是:根据需要,将含噪信号在某一尺度下分解到不同的频带内,然后再将噪声所处的频带置零(或直接提取有用信号所在的频带),进行小波重构,这样就可达到去噪的目的。

2)非线性小波变换阈值法。

当信号和噪声的频带相互重叠时,用小波分解和重构的方法将不能得到令人满意的去噪效果。在这种情况下,非线性小波变换阈值法去噪就可以有效地解决这个问题。非线性小波变换阈值法也称为“小波收缩”(wavelet shrinkage),该方法主要针对信号中混有白噪声的情况,是当前各种小波去噪方法中应用最为广泛的一种。

3)平移不变量小波去噪法。

用阈值法虽然能取得很好的去噪效果,但在有些情况下,不连续点附近的信号会在一个特定的目标水平上下跳变。可以通过平移不变量小波去噪法加以有效地抑制。平移不变量小波去噪方法是在阈值法基础上的改进,若将Donoho的非线性小波变换阈值法称为“第一代去噪方法”,则平移不变量小波去噪法就可称为“第二代去噪方法”。由于一个信号中可能包含若干不连续点,它们之间会相互干扰,即对于一个不连续点的最佳平移可能是另一个不连续点的最差平移,因此,我们通常采用n次循环平移,并将每次平移去噪后的结果进行平均,这得到所谓“平移-去噪-平均”的平移不变量小波去噪法。

4)小波变换模极大值去噪法。

信号的奇异点就是指信号中的突变点,它往往包含了信号的重要特征信息。我们可以观察不同尺度间小波变换模极大值变化的不同规律,去除那些幅度随尺度的增加而减小的点(对应噪声的极值点),保留那些幅度随尺度增加而增大的点,然后再由保留的模极大值点进行重建,即可以达到去噪的目的。

1.4 小波去噪的步骤

1.4.1 阈值法去噪过程的3个步骤

1) 信号的小波分解:

选择一个小波函数,并确定分解的层次,然后进行小波分解的计算。

2)小波分解高频系数的阈值量化:

对高频系数选择合适的阈值进行阈值量化处理。

3) 一维小波重构:

根据小波分解的底层低频系数和各层高频系数进行一维小波重构。

1.4.2 小波分析进行阈值处理方法

1)默认阈值去噪处理。

该方法利用ddencmp生成信号的默认阈值,然后利用函数wdencmp进行去噪处理。

2)给定阈值去噪处理。

在实际的去噪处理过程中,阈值往往通过经验获得,目前这种阈值比默认阈值的可信度高,在进行阈值量化处理时可用函数wthresh。

3)强制去噪处理。

该方法是将小波分解结构中的高频系数全部置为0,即滤掉所有高频部分,然后对信号进行小波重构。这种方法比较简单,且去噪后的信号比较平滑,但是容易丢失信号中的有用成分。

2 基于小波分析的Kalman滤波在GPS数据去噪中的方法研究

2.1 Kalman滤波的基本原理

卡尔曼(Kalman)滤波是由匈牙利数学家卡尔曼博士(Kalman R E)为满足计算机用于人造地球卫星定轨和导航等计算的要求,在1960年提出的一种新的线性滤波模型,通常称它为卡尔曼滤波方法。这种滤波方法很适合处理动态系统的观测数据,因而广泛的用于动态测量,飞行器导航,导弹制导和工业自动化等。

卡尔曼滤波是一种具有无偏性的递推线性最小方差估计,即估计误差的均值或数学期望为零。在计算方法上,卡尔曼滤波采用递推形式,即在t-1时刻估值的基础上,利用t时刻的观测值,递推得到t时刻的状态估值。由于一次仅处理一个时刻观测值,无需存储先前的观测数据,因而计算量大大减少。卡尔曼滤波能有效的消除干扰噪声,获得逼近真实情况的有用信息。同时Kalman滤波是一套计算机实现的实时递推算法,它所处理的对象是随机信号,利用系统噪声和观测噪声的统计特性,以系统的观测量作为滤波器的输入,以所要估计值(系统的状态和参数)作为滤波器的输出,滤波器的输入与输出之间是由时间更新和观测更新算法联系在一起的,根据系统方程和观测方程估计出所有需要处理的信号,所以Kalman滤波是一种最优估计方法。

2.2 基于小波分析的自适卡尔曼滤波方法

经典最优卡尔曼滤波的缺点和局限性是要求精确已知系统的数学模型和噪声统计。然而大多数情况下系统模型的噪声统计是部分已知、近似已知或者完全未知的。应用不精确或错误的模型和噪声统计设计卡尔曼滤波系统可能导致滤波系统性能变坏,甚至使滤波发散。为了解决这个矛后,产生了自适应卡尔曼滤波。自适应卡尔曼滤波的原理是:在进行状态滤波时,利用观测数据提供的信息,在线估计未知噪声的统计,构成状态和噪声统计估计的两段互耦自适应卡尔曼滤波算法。Sage和Husa的自适应卡尔曼滤波算法可在互耦估计状态进行噪声统计(均值和方差阵),算法简单且具有良好的性能,因而被人们广泛应用。它由互耦的常规卡尔曼滤波算法和噪声统计估计器组合而成。Sage和Husa的噪声统计估计器是次优无偏极大后验估值器,适合估计未知常的噪声统计。考虑时变线性离散随机系统的状态方程x(t+1)=Φx(t)+w(t), (4)

式中:x(t+1)为状态变量,y(t)为观测变量,w(t)为系统噪声向量,v(t)为观测噪声向量。

采用的卡尔曼滤波器为

式中: $\hat{x}$ _filter(t)为滤波估计值, $\hat{x}$ _predicted(t)为预测估计值。

考虑到y(t)是一系列带噪声的信号,所以在把y(t)作为观测数据进入卡尔曼滤波系统之前,先将它进行小波去噪预处理,这样可以得到较为平稳的观测数据,再对其进行自适应卡尔曼滤波处理。小波变换实际上是将信号通过低通和高通两组滤波器,平滑分量主要反映信号的低频部分,细节分量主要反映信号的高频部分。首先对信号进行小波分解,一般地,噪声信号多包含在具有较高频率的细节中,从而,可以利用阈值选取方法对所分解的高频信号进行处理,然后对信号进行小波重构即可达到对信号降噪的目的。图1所示为基于小波分析的自适应Kalman滤波方法流程图。

3 结束语

GPS观测值的误差源复杂多样,通过选择合理的应对措施,可以有效的降低各种误差对数据质量的影响,选择合适的小波函数、小波阈值和分解层数对其进行去噪处理,可以明显提高精度;反之,则会降低精度。利用小波去噪法、自适应Kalman滤波法以及基于小波分析的自适应Kalman滤波法对GPS原始观测数据的处理,自适应Kalman滤波法具有不稳定性,对精度总体上有所提高,小波去噪法和基于小波分析的自适应Kalman滤波法对精度提高明显,且后者的稳定性更好。

参考文献

[1]黄声享,刘经南,柳响林.小波分析在高层建筑动态监测中的应用[J].测绘学报,2003,32(2):153-157.

[2]何秀凤,杨光.扩展区间Kalman滤波器及其在GPS/INS组合导航中的应用[J].测绘学报,2004,33(l):47-52.

[3]郭杭.迭代扩展卡尔曼滤波用于实时GPS数据处理[J].武汉测绘科技大学学报,1998,24(2):112-114.

[4]张双成,杨元喜,张勤.一种基于抗差自校正Kalman滤波的GPS导航算法[J].武汉大学学报(信息科学版),2005,30(10):881-884.

[5]黄声享,刘经南.GPS变形监测系统中消除噪声的一种有效方法[J].测绘学报,2002,31(2):104-107.

[6]朱长青,扬启和,朱文忠.基于小波变换特征的地貌影像纹理分析和分类[J].测绘学报,1996,21(4):252-256.

[7]耿则勋,钱曾波.遥感影像无损压缩编码的实验分析与改进[J].Acta Geodaetica et Cartographica Sinica,1996,25(4):262-265.

[8]柳林涛.小波基本理论及其在大地测量等领域中的应用[D].武汉:中国科学院测量与地球物理研究所,1999.

本文来自古文书网(www.gwbook.cn)，转载请保留网址和出处

窄带滤波技术02-02

混合滤波02-02