分析研究基于最小二乘原理

分 析 研 究 基于最小二乘原理的趋势项处理 及其 MA T L A B的实现 王 广斌 ， 刘义伦 ， 金 晓宏， 何玉辉 ( 中南大学 机 电工程 学院 ， 湖南 长沙 4 1 0 0 8 3 ) 摘要趋势项处理是信号检测技术中重要的环节， 本文利用最小二乘原理建立了趋势项 的一般模 型， 并利用 MAT L AB软件强大的矩阵计算功能， 对一个含有非线性趋势项的周期信号进行了趋势项消 除， 并得到了预期效果。 关键词 趋势项 ；最小 二乘 原理 ； MA T L AB 中图分类号T P 3 9 文献标识码A 文章编号1 0 0 3 8 8 8 4 ( 2 0 0 5 ) 0 5 0 0 0 40 5 0 前言 趋势项是振动周期比信号采样长度大 的频率成 分， 一般是由于测试 系统的某些原 因在 时间序列中 产生的一个线性 的或者慢变 的趋势误差 。它是 针对随机信号而言的， 但是在确定性信号 中， 趋势项 也常常存在。产生趋势项的原因很多， 主要有： 采样 时没有对原始信号进行适当的处理， 使采集 的信号 中含有周期 比采样长度长得多的低频成分 ； 由于外 界原因， 包括传感器或仪器的零 点飘移、 基础运动等 引起的信号波形偏 移 ； 由于操作 不当， 信号 经过 积分放大后产生的趋势项， 如零点未调准所产生的 常数， 经积分后成为一条直线： 低噪声经积分放大后 成为缓慢变化的趋势项 ； 在截取记录时， 样本的长度 选择不当。 趋势误差的存在， 会使相关函数、 功率谱函数处 理中出现变形， 甚至可能使低频时的谱 估计完全失 去真实性， 特别对振动或者冲击信号做长时间测量 时， 经常使用测量范围较大 、 测试性能较好的加速度 传感器 J ， 而这种 传感 器测 出 的信号是 加速 度 信 号， 在分析时往往要通过一次积分转化为速度( 用速 度有效值即烈度指标) 、 或两次积分转化 为位移( 用 收稿 日期 2 0 0 4 0 9 0 2 作者简介王广斌( 1 9 7 4一) ， 男， 河南林州人， 讲师， 现为中 南大学博士研究生， 主要研究方向为信号处理与大型回转机 械设备的故障诊断。 4 位移峰值指标) ， 来对设备的性能进行分析和寿命预 测， 在对含趋势项的信号进行积分时， 由于初始条件 不清楚， 积分后的波形会产生很大的畸变， 信号误差 急剧增大， 波形甚至完全扭曲， 使积分后的信号失去 了作为设备工况评判标准的价值 。因此， 在测试 信号分析中常常要消除趋势项， 而在实际工程中， 技 术人员往往因为趋势项拟合计算复杂， 常使用线性 趋势项来代替， 此种处理对 于精度要求不高的场合 是可以的， 但如果要对精密分析信号进行高精度分 析时， 必须用非线性趋势项才能取得很好的趋势项 消除效果。本文利用最小二乘原理建立了趋势项多 项式的一般模 型， 并用 MAT L A B程 序 比较简单 的 消除了信号中的非线性趋势项。 1 最小二乘法建模原理 消除趋势项的方法有多种， 视信号特征、 被测试 对象的物理模型等因素而定。对于随机信号和稳态 信号， 一般采用最小二乘方法， 它既可以消除呈线性 状态的基线偏移， 又可以消除具有高阶多项式 的趋 势项， 也是工程实际中常用的方法。其建模步骤为 ： 首先， 假设一趋势项多项式， 用最小二乘原理列出求 解方程 ； 其次， 用矩 阵法求 出趋势项系数矩 阵， 并得 出趋势项拟合曲线 ； 最后 ， 用原始信号渐去趋势项即 可得出有用信号 。 设 “ 是以 h为采样间隔的数据 采样序列 ( =1 ， 2 ， ， N) ， 现用 K 阶多项式 来拟合趋势项， 令 5 l ： 有 色设 备2 0 0 5 【 5 ) 维普资讯 分 析 研 究 U = b k ( n h ) ( = 1 ， 2 ， 3 ， N ) ( 1 ) 式中 6 一多项式的系数。 u 点的集合是 “ 中多项式元 素的估计 ， 根据 最小二乘原理， 定义中间函数 E( h) 为估计值与真 实值之间的误差 ： E ( h ) 7- ( “ 一 U ) = l “ 一b k ( ) k I ( 2 ) 误差 E( h) 按最 小二乘法求极 小值， 将 式 2对 b j 取偏导数， 并令其为零 ， 则有 ： 薏 = 砉 2 一 K ( n h _ ( n h = 。 ( 3 ) 整理后， 可得出 k+1 个方程如下： b ( n h ) J = fA n ( n h ) j J=0 ， 1 ， 2 ， ， K ( 4 ) 在式( 4 ) 中， 只要求出拟合趋势项系数 b ， 就可 以得 出趋势项的估计多项式， 但对于很大 的 K 值， 按照一般的代数方法计算很复杂， 也很容易出错， 如 果采用矩 阵的方法 ， 利用 MA TL A B运算工 具就 可 以方便快 捷的求 出趋势项 系数 b ， 从 而求得趋 势 项 。 “ “ ( 7 ) 以此类推， 当 K =S ( 0 4 S N 2 ) 时， 得趋势项 系数矩阵： 60 b1 62 ： b “ “ “ n sA h h “ h h ( 8 ) 式( 8 ) 就是趋势项系数矩 阵的模型， 将该模型代 入式( 1 ) ， 就可以得到趋势项多项式的一般模型。 2 趋势项多项式一般模型的递推求解 3 趋势项处理的 MA T L A B仿真 令 n= l ，则 ： 当 K=0时， 得趋势项系数矩阵： b 0 = N “ ( 5 ) 可以看到， 一阶趋势项 b 0 实际上就是均值。 当 K=1时， 得趋势项系数矩阵： ： 当 K=2时， 得趋势项系数矩 阵： f_6 。 N - 1 l 6 l =l 。 I L 6 z J l z 。 。 J 有色设置2 0 0 5 ( 5 ) MA TL A B语言是在 C语言和 F O R TR A N基础 上发展起来的更高级的程序设计语言， 是 1 9 8 4年 由 美国的 Ma t h Wo r k s 推 出的、 为科 学计算和 数据 分 析而专门设计的高级交互式软件包， 它具有强大的 数值计算功能， 可以进行矩阵、 数值插值、 曲线拟合 、 数值微分和数值 积分等运算， 现在 已经成为计算机 辅助设计与分析、 算法 研究和应用开发 的基本工具 和首要平 台。式 ( 8 ) 的数学运 算很 繁琐， 但是 使用 MAT L A B语言对趋势项进行多项式拟合， 则可以大 大减少计算量， 将复杂计算简单化。 现举一例说明之。首先用 MAT A L B构造一个 含线性及非线性 趋势项 的周期信号， 假设采样频率 F 和数据量 N 均为 2 0 0 0 ， 原始含噪声周期信号 为 z、含线性 趋势项的信号为 “ h含二 阶非 线性趋势 项的信号为 “ ， 则原程序如下 ： 5 2 3 4 ； ； 维普资讯 分 析 研 究 F N 2 0 00； = 0： N 一 1； =n F ； h=1 F ； z=8 s i n ( 2*5 0 P * ) +4 c o s ( 2*8 0 P * ) +r a n d n ( s i z e ( ) ) ； 1 = 1 0*+2 0； 1= z + 1 ； 2 =5 0* ： 2+1 0* +2 0；u 2 =z+ 2 ； 对于信号 即含 线性 趋势项 的信号来说， 在 MA TL A B里既可以用所建立的趋势项模型来实现， 也可以用命令 d e t r e n d ， 以如下方式实现： z=d e t r e n d ( 1 ) 但是， 当趋势项是非线性时， 用这个命令就不行 了， 而实际工程中的趋势项大部分是非线性 的， 只能 用最小二乘法建立的趋势项模型来 实现， 对上述信 号 2 对二阶趋势项消除实现程序如下； 含线性趋势项的信号 “ 0 J 4 =s u m( j 4 ) ； s s =s u m( 2 ) ； s 2 =s u m( 2 J ) ； s 3 =s u m( 2 ( J 2 ) ) ； A = N， h J 1 ， h*h J 2 ； J 1 ， h J 2 ， h*h* J 3 ；J 2 ， h J 3 ， h*h J 4 J ； E= s s ； s s 2 ； s 3 ； B=i n v ( A)*E； b d O=B( 1 ) ； b d1 =B( 2 ) ； b d2 =B( 3 ) ； =n F ； t r e nd 2= b dO+ b d1* + b d2* A 2： z 2= 2一t r e nd 2； 图 1 是 以上程序的趋势项消除计算结果， 其中 图 1 a 为对含线性趋势 项的信号 的趋势项消除， 从上到下依此为信号 的曲线 、 线性趋势项曲线、 消除趋势项后的信号 曲线； 图 1 b为对含非线性趋势 项的信号 的趋势 项消除， 从上到下依此为信号 的曲线 、 非线性趋势项 曲线、 消除趋势项后的信 号曲线。 l o 0 5 0 0 含非线性趋势项的信号 “ ： 5 0 0 l 0 0 0 l 5 o 0 20 0 0 0 笪呈 ： 塑塑 ： 一 8 0 0 5 o 0 l 0 0 0 l 5 o 0 消除趋势项后的信号 6 0 4 0 2 0 2 0 l 0 0 - 1 0 -2 0 5 0 0 l 0 0 0 l 5 o 0 20 o 0 信号 的非线性趋势项 5 0 0 l 0 0 0 l 5 o 0 20 o 0 消除趋势项后的信号 0 5 0 0 l 0 o 0 l 5 o 0 2 O 0 0 b二阶非线性趋势项消除 图 1 线性趋势项消除与二阶非线性趋势项消除 注 ： 横坐标为采样次序点， 纵坐标为振动值， m 。 可以看出， 本文对趋势项的拟合效果很好， 可以 将非线性趋势项很好的消除掉， 说 明所建立的非线 性趋势项模型 能够很好的反映信号的变化趋势。在 这个例子中， 信号趋势项是线性和二阶非线性 的， 对 于高阶的情况 ， 只要相应增加趋势项模 型的阶数即 可， 而且低阶趋势项模型是高阶趋势项模型的特例， 建立高阶趋势项模型 ， 则可 以使实际工 程中去趋势 项的效果更为理想。当然， 在实际应用中， 该模型的 阶数也不能取得太大， 否则， 误差会随之增大， 一般 取到三到四阶即可。 6 4 排烟风机振动信号中的非线性趋势 项处理 现以对回转窑排烟风机的振动信号进行趋势项 消除， 来检验本文的算法效果。 中铝某分公司氧化铝厂的回转窑排烟风机是氧 化铝生产中的大型关键设备， 其运行状态的好坏， 直 接影响氧化铝的生产和工人的生命财产安全 。在对 其建立适时监测 系统时， 如果不能有效消除趋势项 ， 就不能正确得 出正确 的监测量数值， 结果可能造成 有色 设备2 0 0 5 ( 5 ) 3 m U l l 2 m U l I II_ 加 0 如 加 加 m 0 m 加 维普资讯 分 析 研 究 误判， 最终导致生产中断或者人身财产损失。 在实际运行 中的风机振动信号中， 其趋势项一 般来说是非线性的， 而且其非线性的形式也很难用 一 定的多项式来表示， 并且随着设备的维护、 仪 器的 劣化、 环境的变化等 ， 其趋势项多项式的结构也在不 断的变化 ， 故上述 方法似不实用， 但是， 从实验可 知， 我们可以用固定的低阶非线性趋势项代替高 阶 1 5 0 1 o o 5 0 0 - 5 0 1 o o 非线性趋势项， 虽然不能消除之， 至少可以在较大程 度上减弱它的影响。由于风机振动信号的不确定性 和长期有效性， 我们在一组风 机振动信号 中加入 了 一 个三阶趋势项， 从 而构成信号 ， 然后对该信号 依照上面 的思想设计程序， 用线性消除法 、 二阶非线 性消除法消除趋势项。 图 2 表示为两种方法消除趋势项后的信号曲线。 含三阶非线性趋势项的风机振动信号 I ! ： “ llI J lr I nJ I 1 _ l i ill“14 I_ I I 丌 IM 1 l 。T 。I1 0 黼I 啊 姗删岍 I 蜘删舯啊虻 明口 胍 r l 邛 1 _ I 1 1 ” 坤 f r r 。 。 1 5 0 0 - 5 0 1 IJIIl IIII血 lI ik 11IIiIl- IIJ 山 IIl l ltl i il l J 1IfIII1 l 】IIIIf1Il 孵 1 l1 1 ll lll p p l 2 1 6 0 0 3 ， 而用二阶非线性趋势项消除法消除趋势项 后， 信号的均值和方差则为 0 5 2 8 8和 2 1 2 9 9 3 ，最 表 1 不同方法消除高阶趋势项后信号均值和方差比较 角色设 置2 0 0 5 ( 5 ) 5 结束语 对于信号 的预处理， 趋势项消除常常不被重视， 而它又在工程实际中真实存 在， 并对测试结果和分 析过程产生很 大影响， 利用 MA TL A B强大 的矩阵 运算功能和本文建立的趋势项多项式 的一般模型， 可使处理后的信号更加准确反映设备运行状况。 另外， 在某些情况下， 趋势项本身就是要知道的 结果， 反映了系统内部的因果关 系， 对设备故障及其 7 维普资讯 分 析 研究 发展有一定反映【 8 ， 通过对信号 中的趋势项提取并 建立模型， 就可以进行设备的运行状态预测和控制， 本文所建立的趋势项多项式一般模型， 也可 以用来 进行趋势预测 并有着广泛的应用场合。 参考文献 1 高频贤 趋势项对 时域参数识别的影响及 消除 J 振 动、 测试与 诊断 。 1 9 9 4 ， ( 6 ) ： 81 0 2 曾禹村， 张宝俊， 吴鹏翼 信号与系统 M 北京： 北京理工大学 出版社 1 9 9 2 北京： 冶金工业出版社， 1 9 9 2 4 寇惠， 原培新 故障诊断中的振动信号处理 M 北京： 冶金工业 出版社， 1 9 8 9 5 李少远， 曹保定， 孟昭终， 王仙琼 基于时间序列分析方法的预测 模型研究 J 河北工学院学报， 1 9 9 5 ， ( 3 ) ： 7 1 1 6 陈怀深， 吴大正， 高西全 MA T L A B及在电子信息课程中的应用 ( 第 2 版