递推最小二乘算法

递推最小二乘算法（精选7篇）

递推最小二乘算法 篇1

摘要：递推最小二乘RLS (Recursive of Least Square) 算法是自适应滤波算法中的精确分析算法。它具有收敛速率快, 精确度高等特点, 但是发现目前RLS算法多用于对一维信号的去噪处理。使用递推最小二乘 (RLS) 算法对二维图像进行去噪, 从处理一维信号变成处理二维图像信号, 需要对RLS算法进行改进。先迭代得到滤波器参数, 形成3×3滤波掩模, 再改进算法对图像进行滤波;同时与常数比率维纳滤波和自相关函数的维纳滤波算法的去噪效果进行对比。结论证明在对图像进行较严重的模糊和加噪处理后, 其他两种算法对图像的还原能力差, 而递推最小二乘自适应滤波 (RLS) 算法具有优良的图像去噪性能。

关键词：自适应滤波,递推最小二乘算法 (RLS) ,图像去噪,MATLAB仿真

0 引言

自适应滤波器是能够根据环境的改变, 使用自适应算法来改变滤波器的参数和结构这样一类滤波器。这方面一些有价值的早期研究工作是由联邦德国的卡尔斯卢埃技术学院和斯坦福大学的研究人员完成的, 他们从1959 年起就开始研究自适应模式识别系统。在20 世纪60 年代中叶, 鲁丁 ( Rudin) 发表了有关自适应滤波器进展现状的出色综述及其在自动均衡应用方面的早期参考资料。最小均方 ( LMS) 自适应算法的进展源于斯坦福大学的模式识别研究, 并于1967 年首次由威德洛 ( Widrow) 发表在自适应阵列的论文中。接着在1971 年他又提出自适应滤波器方面的报告。自适应滤波器的特点是[1]: 滤波器的参数可以自动地按照某种准则调整到最佳滤波; 实现时不需要任何关于信号和噪声的先验统计知识, 尤其当输入统计特性变化时, 自适应滤波器都能调整自身的参数来满足最佳滤波的需要。常常将这种输入统计特性未知, 调整自身的参数到最佳的过程称为“学习过程”。将输入信号统计特性变化时, 调整自身的参数到最佳的过程称为“跟踪过程”, 因此自适应滤波器具有学习和跟踪的性能。

自适应滤波器发展很快, 已广泛地用于系统模型识别, 通信信道的自适应均衡, 雷达与声纳的波束形成, 减少或消除心电图中的周期干扰, 噪声中信号的检测、跟踪、增强和线性预测等[2 - 7]。由于最小均方自适应算法采用梯度搜索法, 使得收敛速度比其它算法快, 所以现在这一算法已广泛地用于计算自适应滤波器的权系数[8 - 10]。

目前采用的较广泛的自适应算法主要有递归最小二乘 ( RLS) 算法、Widrow-Hoff最小均方 ( LMS) 算法、格型滤波器算法和无限冲激响应 ( ⅡR) 算法等[1]。递推最小二乘 ( RLS) 滤波器是一种由LMS算法演变而来的顺序衰退滤波器, 其基本原理是使预测误差最小化。由于RLS算法的初始收敛速率远大于LMS算法, 受到遗忘因子 λ 的影响远比 μ 对LMS算法的影响小得多。

虽然自适应滤波器有了长足发展和广泛应用, 但是在图像处理尤其图像去噪方面的应用十分稀少, 而经过研究发现, 在没有图像统计特性先验知识的情况下, 自适应滤波器能够自动跟踪并抑制图像中的相关噪声, 对既存在模糊[11]又存在噪声的图像去噪具有较大优势, 尤其图像已被噪声污染较重, 维纳滤波器等常规滤波器已无法还原图像的情况下, 自适应滤波器常常有较好的图像还原效果。因此, 本文使用递推最小二乘 ( RLS) 自适应算法对图像进行去噪, 同时与常数比率维纳滤波和自相关函数的维纳滤波算法的去噪效果进行对比, 结论证明RLS算法在对图像的去噪效果方面优于其他两种算法。

1 最小二乘自适应滤波

最小二乘LS ( Least Square) 准则[1]为:

式中, ξ ( n) 是误差信号的平方和; ej是j时刻的误差信号, 即:

其中, dj是j时刻的期望信号, Xj是j时刻的输入信号构成的向量, W表示滤波器的权系数构成的向量。通过选择W, 使 ξ ( n) 取得最小值的滤波称为最小二乘LS滤波, 而满足E[ej2]取得最小值的滤波称为最小均方误差LMS ( Least Mean Square) 滤波。和LMS滤波相比, LS滤波对非平稳信号的适应性要强许多, 这是由于LS滤波总是采用新的准则, 在每一个时刻对所有已输入信号而言, 重新评估使其误差的平方和最小, 因此具有更精确的含义, 属于精确分析法。而LMS滤波是以集合平均为基础的, 属于统计分析的方法。

1. 1 最小二乘滤波

已知n个数据{ x ( 1) , x ( 2) , …, x ( n) } , 采用M个权的FIR滤波器对数据进行滤波, 假设期望信号为d ( i) , 如图1 所示。

图 1 M 个权的 FIR 滤波器

滤波器的输出是对期望信号d ( i) 的估计, 即:

其中, wk ( i) , k = 1, 2, …, M为FIR滤波器在i时刻的M个系数值 ( 说明滤波器的系数可以变化) , 它是一个M维的向量, 记为wM ( i) =[w1 ( i) , w2 ( i) , …, wM ( i) ]T。同理, 输入信号也是一个M维的向量, xM ( i) =[x1 ( i) , x2 ( i) , …, xM ( i) ]T。n时刻, 估计误差为:

误差信号的平方加权和为:

式中, 由于数据长度有限, 对观测区间以外的数据所做的约定不同, 当i的取值范围不同时, 得到不同的 ξ ( n) 。这里采用前加窗法, 约定:

得到:

式中, λ 称为遗忘因子, 主要用于增加新数据的权重, 以增强对非平稳信号的适应性, 一般取值在0. 95 至0. 995 之间。ξ ( n) 为误差信号的能量。

要使 ξ ( n) 取得最小值, 满足:

式 ( 7) 成立的w ( n) 就是wM ( n) 的最小二乘估计, 记为^wLS ( n) 。

1. 2 递推最小二乘法 ( RLS)

这里仅对经典的一次完成的最小二乘算法作了介绍, 在具体使用时, 这种算法不仅占用内存大, 而且不能实时实现, 因此, 现在介绍经常采用的递推最小二乘算法[12]。

递推最小二乘法RLS的基本思想就是新的估计值是在老的估计值的基础上修正而成的, 为了分析简单, 便设 θ 为一向量, 且 θ 仅与当前观测值有关, 则新的估计值。

最小二乘递推算法的一种形式为:

也可以设:

因此可以得到递推最小二乘的另一种形式为:

式中, z ( k) 为k时刻期望信号, a ( k) 为k时刻输入信号, P ( k) 、q ( k) 为中间量。

2 适用于图像去噪的RLS算法设计

2. 1 基本思想

从上述理论得到的RLS算法是适用于一维信号去噪的算法, 现在要对RLS算法进行改进, 以适应对二维图像信号进行去噪。

由式 ( 7) 得到的w ( n) 记为, 即滤波参数, 是一个列向量, 把赋给一个3 × 3 的滤波掩模如图2 所示, 需要把滤波掩模先变成一个列向量, 然后才能赋值, 用这个3 × 3 的掩模充当滤波参数, 对二维图像进行滤波。

图 2 3 × 3 的滤波掩模

另一方面, 更重要的是, 处理一维信号的RLS算法在计算机程序化设计中只需要一层循环, 而对于二维图像, 需要逐行扫描或逐列扫描的方法进行滤波, 这就要把RLS程序化算法的一层循环改进成两层的循环嵌套, 以实现用RLS算法对图像信号进行滤波去噪。

2. 2 具体实现方法与核心编程

( 1) 对一幅图像 ( 大小为512 × 512 像素分辨率为72dpi的8 位无符号整形灰度图像) 用fspecial函数[11] ( len = 7, theta =45°) 进行模糊, 又对图像加均值为0 方差为0. 01的高斯噪声。

( 2) 把处理后的图像g1 赋给一个变量, 程序为:

RLS滤波后的图像如图6 所示。

3 仿真实验与结果对比

本仿真实验的硬件条件为ASUS A43S系列计算机一台, 软件为Math Works公司的MATLAB R2011b软件, 实验所用图像数据来源于文献[11], 图像特征为: 大小512 × 512 像素, 分辨率72dpi的8 位无符号整形灰度图像, 并且改变图像特征 ( 如分辨率) 对本文三种滤波算法的图像滤波效果基本无影响。

采用递推最小二乘 ( RLS) 算法、常数比率维纳滤波和自相关函数的维纳滤波三种算法对模糊和加噪后的图像进行滤波, 并对这三种滤波算法的效果进行对比。

仿真实验流程如图3 所示。

具体仿真实验步骤为:

(1) 图4为原图像。

( 2) 用fspecial函数[11] ( len = 7, theta = 45°) 对原图像进行模糊, 又对图像加均值为0 方差为0. 01 的高斯噪声, 模糊且加噪后的图像如图5 所示。

( 3) 采用3 × 3 滤波掩模的递推最小二乘 ( RLS) 算法对模糊和加噪后的图像 ( 图5) 进行滤波, 滤波后的图像如图6 所示。

( 4) 对图6 图像计算均方误差 ( MSE) 和峰值信噪比 ( PSNR) 。

( 5) 用常数比率维纳滤波算法对模糊且加噪的图像 ( 图5) 进行滤波, 并计算此种方法的MSE和PSNR, 滤波后的图像如图7 所示。

( 6) 用自相关函数维纳滤波算法对模糊和加噪的图像 ( 图5) 进行滤波并计算MSE和PSNR, 滤波后的图像如图8 所示。

三种算法处理图像后的MSE和PSNR值见表1 所示。

对比图6图7和图8可看出图6即递推最小二乘 ( RLS) 算法对图像的去噪效果最好, 图8 即采用自相关函数维纳滤波算法滤波后的图像可大致看出人物轮廓, 图7 采用常数比率维纳滤波算法滤波后的图像已无法辨认图像内容, 而表1的数据更加精确地证明了这一直观判断。

由表1 看出递推最小二乘 ( RLS) 滤波算法在三种滤波算法中的MSE最小且PSNR最大, 所以递推最小二乘 ( RLS) 滤波算法对进行较严重模糊和加入高斯噪声的图像具有优良的去噪性能。

4 结语

针对目前RLS算法多用于对一维信号去噪处理的现状, 对RLS算法进行改进, 从程序化算法的一层循环改进为两层循环嵌套, 实现了对二维图像的去噪功能, 由以上MATLAB仿真实验可得出在对图像较严重的模糊和加噪处理后, 递推最小二乘 ( RLS) 滤波算法相对于常数比率维纳滤波和自相关函数的维纳滤波算法具有优良的二维图像去噪性能, 但是就递推最小二乘 ( RLS) 滤波算法本身滤波效果来看, 还有很大改进空间; 且由于RLS滤波算法应用于图像处理采用的是像素点信息作为输入参数, 其造成矩阵运算计算量过大、处理时间长等缺点, 下一步将在这些方面进行改进。

参考文献

[1]丁玉美, 阔永红, 高新波.数字信号处理—时域离散随机信号处理[M].西安:西安电子科技大学出版社, 2002.

[2]Yong Yin, Nenghai Yu Ph D, Weijie Dong.A research of UWB rake receiver based on novel RLS adaptive algorithm[J].Journal of Electronics (China) , 2006, 3:341-345.

[3]吴雪芬, 陈功, 朱锡芳.基于小波和RLS的锂电池激光在线测厚[J].通信技术, 2012 (7) :90-93.

[4]Xiangyuan Bu, Yaoyao Yang.Adaptive equalization using a new RLS algorithm combined with the coordinate descent method[J].Advanced Technology in Teaching, 2013, 163:631-638.

[5]王金明, 张雄伟.一种基于自适应模糊滤波的语音增强方法[J].解放军理工大学学报:自然科学版, 2003 (1) :17-20.

[6]Zakharov Y V, White G P, Liu J.Low-Complexity RLS algorithms using dichotomous coordinate descent iterations[J].IEEE Transactions on Signal Processing, 2008, 56 (7) :3150-3161.

[7]陶青, 谢勤岚.基于RLS算法的脑电信号的去噪[J].现代电子技术, 2006, 11 (226) :143-145.

[8]Bharath Sridhar, I Akram Sheriff, Dr K A Narayanan Kutty, et al.Comparison of cascaded LMS-RLS, LMS and RLS adaptive filters in Non-Stationary environments[J].Novel Algorithms and Techniques in Telecommunications and Networking, 2010:495-499.

[9]Rongshan Yu, Chi Chung Ko.Lossless compression of digital audio using cascaded RLS-LMS prediction[J].IEEE Transactions on Speech and Audio Processing, 2003, 6.

[10]Farhangboroujeny B.Variable-step size LMS algorithm:new developments and experiments[C].IEE Proceedings Visual Image Signal Processing, 1994, 5:311-317.

[11]Rafael C Gonzalez.数字图像处理 (MATLAB版) [M].阮秋奇, 等译.北京:电子工业出版社, 2005.

[12]Muneyasu M, Uemoto E.A novel 2-D adaptive filter based on the 1-D RLS algorithm[C]//Proceedings of the 1997 IEEE Intertional Symposium on Circuits and Systems.Hong Kong, 1997:2317-2320.

递推最小二乘算法 篇2

最小二乘算法是系统辨识中用得最广泛的估计方法之一。标准的递推最小二乘算法是通过极小化关于输入输出数据的一个二次准则函数,即极小化估计残差平方和而得到的算法[1]。本文基于卡尔曼滤波原理[2],从一个新的角度来推导递推最小二乘算法:定义一个参数估汁误差协方差矩阵,通过极小化该协方差矩阵而推导出递推最小二乘算法。

1最小二乘辨识算法

考虑下列单输入单输出(SISO)系统的参数估计问题,

A(z)y(t)=B(z)u(t)+v(t)(1)

(1)式中{u(t)}和{y(t)}分别是模型的输入和输出序列,{v(t)}为零均值、方差为σ2的不相关随机白噪声序列A(z)和B(z)均为单位后移算z-1的多项式[z-1y(t)=y(t-1)],且

A(z)=1+a1z-1+a2z-2+…+anaz-n,

B(z)=b1z-1+b2z-2+…+bnbz-nb。

定义参数向量θ和信息向量φ(t)分别为

θ:=[a1,a2,…,ana,b1,b2,…,bnb]T∈Rn,

φ(t):=[-y(t-1),-y(t-2),…,-y(t-na),u(t-1),u(t-2),…,u(t-nb)]T∈Rn。

那么式(1)写可成下列最小二乘辨识模型,

y(t)=φT(t)θ+v(t)(2)

本文的目标是,利用系统的输入输出数据{u(t),y(t)}或{y(t),φ(t)},基于卡尔曼滤波原理,极小化参数估计误差协方差阵,而推导辨识参数向量θ的递推最小二乘算法。

很多文献给出了计算θ估计白$\widehat{\theta }$(t)的递推最小二乘算法[1,3]:

$\widehat{\theta }(t)=\widehat{\theta }(t-1)+L(t)[y(t)-\phi {}^{\text{Τ}}(t)\widehat{\theta }(t-1)]$ (3)

L(t)=P(t-1)φ(t)[1+φT(t)P(t-1)φ(t)]-1(4)

P(t)=[I-L(t)φT(t)]P(t-1)(5)

其中$\widehat{\theta }$(t)为t时刻的估计,L(t)称为增益向量,P(t)称为协方差矩阵。

2算法推导

2.1单输入单输出系统的递推最小二乘算法

仍考虑上述单输入单输出系统式(1)或式(2)。参照式(3),设t寸刻的参数估计$\widehat{\theta }$(t)足由前一时刻参数估计$\widehat{\theta }(t-1)$加上一个补偿项得到,且具有下列形,

$\widehat{\theta }(t)=\widehat{\theta }(t-1)+L(t)[y(t)-\phi {}^{\text{Τ}}(t)\widehat{\theta }(t-1)]$ (6)

那么问题就是如何找到增益向量L(t)∈Rn。方法如下,定义参数估计误差为

$\widehat{\theta }$(t):$=\widehat{\theta }(t)-\theta$ (7)

将式(6)代入式(7),并利用式(2)展开可得

$\widehat{\theta }(t)=[{\rm I}-L(t)\phi {}^{\text{Τ}}(t)]\widehat{\theta }(t-1)+L(t)v(t)$ (8)

(8)式中I是一个适当维数的单位阵.目标是确定一个最优增益向量L(t)使参数估计误差$\widehat{\theta }$(t)最小。由于v(t)是白噪声,假设L(t)和φ(t)与v(t)独立,对式(8)两边取期望得$\text{E}[\widehat{\theta }(t)]=0(\text{E}$为期望算子)。定义参数估计误差协方差阵.

P(t):$=\text{E}[\widehat{\theta }(t)\widehat{\theta }{}^{\text{Τ}}(t)]$ (9)

将式(8)代入式(9)可得

P(t)=[I-L(t)φT(t)]P(t-1)×

[I-L(t)φT(t)]T+L(t)σ2LT(t)(10)

因为P(t)为非负定矩阵,将P(t)配成下列形式:

把P(t)看作待定增益向量L(t)的函数。通过极小化估计误差协方差矩P(t),可求得最优增益向量L(t).上式中矩阵P(t)包含3项,第一、二项与L(t)无关,第二项中,因为P(t-1)是非负定矩阵,所以σ2+φT(t)P(t-1)φ(t)≥0。如果选择增益L(t)使得第三项为零,即取

L(t)=P(t-1)φ(t)[σ2+φT(t)P(t-1)φ(t)]-1(12)

那么协方差阵P(t)最小,也就是参数估计误差最小。

$\text{E}[\parallel \widehat{\theta }(t)\parallel {}^2]=tr[{\rm P}(t)]=\min$。

将式(12)代入式(11)得到

P(t)=P(t-1)-P(t-1)φ(t)×

[σ2+φT(t)P(t-1)φ(t)]-1φT(t)P(t-1)=

[I-L(t)φT(t)]P(t-1)(13)

式(6),式(12)和式(13)构成了基于卡尔曼原理的单输入单输出系统的最小方差递推最小二乘算法它与标准递推最小二乘算法式(3)-式(5)的差别在于式(12)-式(13)中取σ2=1.

下面讨论一种特殊情况。对于确定性系统,噪声方差σ2=0,则算法式(6),式(12)-式(13)退化为

$\widehat{\theta }(t)=\widehat{\theta }(t-1)+L(t)[y(t)-\phi {}^{\text{Τ}}(t)\widehat{\theta }(t-1)]$,(14)

L(t)=P(t-1)φ(t)[φT(t)P(t-1)φ(t)]-1(15)

P(t)=[I-L(t)φT(t)]P(t-1)(16)

为了防止式(15)中分母为零,可加上一个小常数ε(如取ε=10-6),则算法式(14)-式(16)可修改为

$\widehat{\theta }(t)=\widehat{\theta }(t-1)+L(t)[y(t)-\phi {}^{\text{Τ}}(t)\widehat{\theta }(t-1)]$ (17)

$L(t)=\frac{{\rm P}(t-1)\phi (t)}{\epsilon +\phi {}^{\text{Τ}}(t){\rm P}(t-1)\phi (t)},0＜\epsilon \le 1$ (18)

P(t)=[I-L(t)φT(t)]P(t-1)(19)

在后面的例子中,将研究ε取不同值时,参数估计精度的变化。

这是文献[2]状态估计算法推导过程,推导递推最小二乘参数估计的一种新的简便方法。从上述推导看:这种参数状态估计误差具有最小方差性质。也可理论上证明这个递推最小二来算法的参数估计误差收敛于零。

2.2多输入多输出系统的递推最小二乘算法

上述方法也叮推广到多输入多输出(MIMO)系统.考虑下列多输入多输出系统,

A(z)y(t)=B(z)u(t)+v(t)(20)

A(z):=I+A1z-1+A2z-2+…+Anaz-na,

B(z):=B1z-1+B2z-2+…+Bnbz-nb,

u(t)∈Rr为系统输入向量,y(t)∈Rm为系统输出向量,u(t)∈Rm为零均值、协方差阵为Rv=E[v(t)vT(t)]∈Rm×m的随机噪声向量。

定义参数矩阵ϑ和信号向量φ(t)如下

ϑT:=[A1,A2,…,Ana,B1,…,Bnb]∈Rm×m,φ(t):=[-yT(t-1),…,-yT(t-ma),tT(t-1),…,uT(t-nb)]T∈Rn。

则式(20)可等价写为

y(t)=ϑTϕ(t)+v(t)(21)

参照式(6),由上式可知ϑ的递推估计应取为下列形式:

$\widehat{\vartheta }(t)=\widehat{\vartheta }(t-1)+L(t)[y{}^{\text{Τ}}(t)-\text{ϕ}{}^{\text{Τ}}(t)\widehat{\vartheta }(t-1)]$ (22)

则参数估计误差可写为

$\begin{array}{l}\widehat{\vartheta }(t):=\widehat{\vartheta }(t)-\text{ϑ}=[{\rm I}-L(t)\text{ϕ}{}^{\text{Τ}}(t)]\times \\ \widehat{\vartheta }(t-1)+L(t)v{}^{\text{Τ}}(t).\end{array}$

定义参数估计误差协方差短阵P(t):$=E[\widehat{\vartheta }(t)\times \widehat{\vartheta }{}^{\text{Τ}}(t)]$。令ϵ:=tr[Rv]=E[‖v(t)‖2].采用与上节相同的方法,可以得到多输入多输出系统的最小方差递推最小二乘算法:

$\widehat{\vartheta }(t)=\widehat{\vartheta }(t-1)+L(t)[y{}^{\text{Τ}}(t)-\text{ϕ}{}^{\text{Τ}}(t)\widehat{\vartheta }(t-1)]$ (23)

$L(t)=\frac{{\rm P}(t-1)\text{ϕ}(t)}{ϵ+\text{ϕ}{}^{\text{Τ}}(t){\rm P}(t-1)\text{ϕ}(t)}$ (24)

P(t)=[I-L(t)ϕT(t)]P(t-1)](25)

如果式(24)中取ϵ=1,就得到多输入多输出系统的标准递推最小二乘算法.

3仿真例子

例1 考虑下列随机系统,

A(z)y(t)=B(z)u(t)+v(t),

A(z)=1+a1z-1+a2z-2=1-1.60z-1+0.80z-2,

B(z)=b1z-1+b2z-2=0.412z-1+0.309z-2,

θ=[a1,a2,b1,b2]T=[-1.60,0.80,0.412,0.309]T。

仿真时,输入u(t)采用零均值单位方差不相关可测随机序列,v(t)采用零均值方差为σ2的白噪声序列.当噪声方差为σ2=0.102和σ2=1.002时,对应的噪信比分别为δns=14.26%和δns=142.63%.应用提出的算法式(6),式(12)-式(13)估计这个系统参数,不同噪声方差下的参数估计及其误差如表1所示,参数估计误差δ:$=\parallel \widehat{\theta }(t)-\theta \parallel /\parallel \theta \parallel$随t变化曲线如图1所示.

当噪声方差为σ2=0.102,噪信比为δns=14.26%时,利用算法式(17)-式(19)来估计这个系统参数,不同ϵ下的信真结果如表2所示.

从表1和图1可知,在不同噪声水平下,噪信比越小(噪声方差小),参数估计精度高,参数估计误差随数据长度增加总的趋势是减小的。从表2可知,在同一噪声水平下,算法参数估计对ϵ取值大小不敏感。

例2 考虑下列确定性系统,即例v(t)=0情形。输入v(t)采用零均值单位方差不相关随机序列,应用算法式(17)-式(19)估计其参数,不同ϵ值时的仿真结果如表3所示。

4结论

本文基于卡尔曼滤波原理,通过极小化参数估计误差协方差矩阵的方法,分别推导了单输入单输出系统和多输入多输出系统的递推最小二乘算法,并用仿真例子研究了参数估计的性质。

摘要：基于卡尔曼滤波原理,通过极小化参数估计误差协方差矩阵,导出了递推最小二乘辨识算法。仿真例子说明了算法的有效性。

关键词：最小二乘,参数估计,参数估计误差,协方差

参考文献

[1]谢新民,丁锋.自适应控制系统.北京:清华大学出版社,2002

[2]str om KJ.Introduction to stochastic control theory.NewYork:Ac-ademic Press,1970

递推最小二乘算法 篇3

三维激光扫描测量技术自出现以来就以快速性, 不接触性, 穿透性, 实时、动态、主动性, 高密度、高精度, 数字化、自动化等特性广泛应用到测绘领域, 通过主动测量方式获取观测目标对象表面丰富的信息, 不仅包含目标对象的三维坐标, 还包括物体表面的反射率、纹理等信息, 这些信息具有数据量大、种类多、速度快等特性, 在点云数据中以模拟的形式体现出来。由于数据量非常庞大, 在计算机中存储、运算需要大量的存储空间, 且调用数据的速度比较慢, 且通过模拟点形式体现目标形态, 因此大多数情况下目标对象的矢量化信息不能直接读取和加以利用[1]。

目前针对点云滤波的方法概括起来主要是以下几种:最大局部坡度滤波方法、增大窗口高程阈值滤波方法、移动曲面拟合滤波方法、区域增长滤波方法、数学形态学滤波方法等[2]。其中, 坡度理论的滤波方法设定的阀值需要大量的人工试验, 移动曲面拟合法在地势平坦的区域取得很好的滤波效果, 在地形起伏较大的区域效果较差。

本文采用移动窗口的最小二乘曲面拟合算法, 在地形数据滤波中表现良好, 具有适用于各种特征地形, 运行稳健、分类精度高等特点, 在地形起伏变化大的地区也可以得到良好的滤波效果, 不足之处是需要多次迭代。

1 移动窗口的最小二乘曲面拟合算法

移动窗口的最小二乘法曲面拟合算法是一种逼近理论, 是将采样数据通过格网化处理, 选取一定格网数据, 利用最小二乘曲面拟合法对对采样数据进行拟合计算, 得到采样数据拟合值, 满足拟合曲面在取样处的拟合值与实际值之差的平方和最小。

Axelsson等人提出了使用TIN作为参考DTM的移动窗口法, 这里使用六个已知点解算出的二次曲面代替TIN作为拟合的局部参考DTM[7]。

1.1 基本原理

移动窗口的最小二乘曲面拟合法算法的基本原理是在一定区域内, 使用给定的拟合曲面作为参考数字地面模型, 对区域内的点云数据进行判断。通过比较实际数据点高程值和参考的数字地面模型值的大小, 来实现数据的滤波, 保留阀值内的实际地形点, 滤除大于阀值的地物点。该算法在满足以下两个假设前提的条件下才能够成立:第一是在一定的区域范围内高程值最低点必须为地形点;第二是地面的点云数据在一定的区域范围内的分布情况是符合二次曲面分布的, 而其他的地物点则不在这个曲面上, 并且高于该曲面。在以上两个假设成立的前提下, 使用一个固定的窗口获取其范围内的六个高程最低点值, 并作为地形点解算出对应的曲面拟合方程作为参考的数字地面模型。将此作为依据, 对其他数据进行处理, 通过移动窗口的多次迭代计算, 逐步缩减窗口的大小和高程阀值的大小, 滤波剔除掉原始数据中的地物点, 最终获得全部是地形点的数据点集。

1.2 滤波应用

移动窗口最小二乘曲面拟合需要将点云数据进行格网化处理, 即按照一定间隔进行格网划分, 如图1所示, 这里使用长宽都为20m的格网对所有原始数据进行格网化处理。设定初始格网的大小为2×3的格网 (40m×60m) , 从首选窗口中确定六个格网中每个格网的最低点高程值, 并利用这六个点进行二次曲面拟合方程计算, 作为局部参考的数字地面模型, 将该窗口 (2×3个格网) 中的每个点云数据带入该曲面方程, 判断由拟合方程计算出来的值与实际值之间的差是否大于阀值, 若大于阀值, 则认为该点是地物点, 进行剔除, 若小于阀值, 认为该点是地面点, 予以保留。

判断完毕后, 移动到下一个六个格网处重复上述的判断, 按照顺序依次移动窗口, 直到窗口到达最后六个格网单元结束。初次滤波后, 得到初始的地形点集Points1。对得到的点集Points1重复进行移动窗口最小二乘曲面拟合, 此时缩小网格单元尺寸, 减小高程阀值, 得到新的地形点集Points2。往返多次, 最终可以得到较好的DTM。

对于格网单元内没有点的处理:把一个格网平面覆盖在点云数据上, 大多数的格网内是有数据的, 而且点云的个数很多, 但对于点云区域的边界或是由于建筑物的遮挡, 或是水面的吸收, 有些格网中没有点云数据, 这会导致曲面拟合方程无法构建或是出错。此时需要设置一个虚拟的最低点, 参与曲面方程的构建。

若格网中有两个格网内没有点, 则可以用上述的方法用其周围的一个点的高程作为一个虚拟点的高程, 但是有的时候, 格网单元内没有点的个数很可能大于两个, 此时, 要用其周围点的高程做虚拟高程, 拟合的面可能会成为平面, 与实际地形不符, 所以, 在这里当格网内无点个数大于两个的时候直接进行下一个格网的拟合。

1.3 算法的实现

该算法的具体实现过程如下:

(1) 首先通过对原始的点云数据的预处理进行极大极小点云的剔除;

(2) 在平面坐标系内, 找出两坐标的最大、最小值, 然后根据X、Y坐标的最大、最小值将原始点云数据进行格网划分, 根据设计大小的格网单元窗口对预处理后的点云数据进行最小二乘曲面拟合的迭代计算, 获取初始化的地面点集合。

(3) 通过缩减移动窗口的大小和高程阀值, 对初始化的地面点集进行第二次的移动窗口最小二乘曲面拟合迭代计算, 得到新的地形点集。

(4) 根据实际的滤波精度需要, 对新得到的数据点集进行多次重复迭代计算, 直到得到满意的结果为止。

1.4 试验分析

本文以徕卡地面三维激光扫描测量系统Scan Station2获取的地形点云数据为研究对象, 实验点云总数183155个, 覆盖面积17143.67478m2, 见图3。

由图3可见, 本试验区域范围内植被比较低矮, 地面具有一定的坡度, 故进行滤波时对滤波格网和高程阀值的设置不宜过大, 第一次设置滤波格网间距为8m, 高程阀值为1.5, 对点云进行滤波处理, 从图4中的点云模型可以看到, 滤波效果比较明显, 大部分非地面点已经被过滤掉, 很好地保持了地形特征。图5为经过移动窗口最小二乘曲面拟合滤波调整后, 设置滤波格网间距为4m, 高程阀值为0.2的模型显示, 该图中过滤了大部分第一次滤波中没有过滤掉的低矮植被。图6为经过二次移动窗口最小二乘曲面拟合滤波调整后, 设置滤波格网间距为1m, 高程阀值为0.05的模型显示, 从图6中可以明显地看出, 和原始点云相比, 植被覆盖及其他大量噪声点已经被剔除, 滤波效果较为理想。但是, 在用移动窗口最小二乘曲面拟合滤波处理前需要对极低点进行剔除, 否则, 在有极低点存在的那个窗口内, 会将大部分的地面点滤除, 这样就不能够得到较好的滤波效果。

经过三次移动窗口设置不同窗口大小和阀值, 采用最小二乘曲面拟合算法对点云进行滤波处理, 对处理后的点云数据进行精简与DEM模型重建, 从上述重建的模型上来看, 对于非地面的分类选择与剔除已基本达到要求, 且很好地保留了地形特征, 由于该算法对窗口大小和阀值设置还有一定的经验存在, 故主观性比较大, 需要对数据进行多次处理, 才能达到理想的效果。从图4与图6模型对比分析可以看出, 经过多次迭代滤波处理后, 地面点中的非地面点都被滤除了, 滤波效果比较理想。

2 点云滤波的质量评价[1]

2.1 定性评价

由于现实世界中, 相邻点之间的空间结构上有一个明确的联系, 然而, 在分类时, 由于对这种联系的未知, 导致在进行分类时不能充分利用好点之间的联系, 产生分类错误。如果只考虑将点云数据分为两类, 一是表示地形特征的数据;一是表示地物信息的数据。我们也可以分为两种分类错误, Ⅰ型错误 (将地面点分为地物点的分类错误) 和Ⅱ型错误 (将地物点分为地面点的分类错误) 两种类型[1]。

对点云数据用上述移动窗口最小二乘曲面拟合滤波算法, 通过对处理后的地面点集生成的模型的分析, 获得对滤波质量的评价, 从模型上的分析可以看出, 对于该地形区域, 利用移动窗口最小二乘曲面拟合的方法滤波效果较好, 基本上过滤掉地物点, 能很好保留地形特征。

2.2 定量评价

对于点云数据滤波质量的定量评价, 本文基于以下假设成立时, 采用了如下的评价指标:假定a表示的是分类处理时数据集中的地面点的总数;b表示的是分类数据集中非地面点总数;c表示分类过程中将地面点误分为地物点的总数;d表示分类时将地物点误分为地面点的总数。根据以上假设, 定义了以下三个评价指标:。

在实际情况下, 地形都是比较复杂的, 因此为了得到高精度的数字地形模型, 通常在进行滤波处理时, 首先要根据实际需要选择滤波参数, 达到所需的滤波要求。如对滤波结果的比较, 选择较小的一类误差, 但也不能忽略二类误差。总误差是衡量两类错误总数的指标, 所以通过总误差的比对, 选择更适合的数据。

本文通过植被比较密集的区域数据点分类及滤波处理, 得到一系列结果 (如表1所示) , 并通过所得到的结果对上述滤波方法进行定量评价:从一类误差的比较中可以看出, 移动窗口最小二乘拟合方法得一类误差较大, 说明在滤波处理时, 被误分为地物点的地面点较多;从二类误差的比较中可以看出, 移动窗口最小二乘拟合方法的二类误差较小, 说明在滤波的过程中被误分为地面点的地物点较少。

3 结论

对于点云数据的滤波处理, 一直是点云数据处理研究的热点和难点, 本文利用MATLAB语言编程实现了基于移动窗口的最小二乘曲面拟合滤波算法, 对试验数据进行了滤波处理, 同时结合徕卡测量系统自身配备的点云数据处理软件cyclone, 对点云数据进行预处理、点云拼接、点云精简、点云分割、平滑处理等操作, 取得了一定的效果。基于移动窗口的最小二乘曲面拟合滤波算法在点云数据滤波中表现较好, 具有适用于各种特征地形, 运行稳健、分类精度高等特点, 在地形复杂且连续变化的地方能够获得很好的滤波效果, 该算法的不足之处就是对数据进行滤波处理时需要对数据进行多次的迭代计算。

本文所做实验只是对该特定的区域进行的, 通过对滤波后的实验数据进行定量分析, 综合评价了上述滤波方法, 但是从总体上来说, 由于这种方法本身还存在一定的局限性, 因而还不能完全将地面点与地物点分离开来, 因此还需要在今后的研究中进一步的改进。

摘要：点云滤波是地面三维激光扫描数据处理的重点和热点问题, 也是点云数据处理的一项关键工作。本文在阐述了目前几种主要的点云滤波方法的基础上, 提出了基于移动窗口最小二乘曲面拟合算法, 详细介绍了该算法的基本原理、实现过程及其应用。该算法具有拟合曲面精度高、光滑性好等特性, 经试验数据分析, 该算法在点云数据率处理中取得了良好的滤波效果, 能够有效对点云进行分类, 去除点云噪声, 很好保持了地形特征, 获得高精度建模数据。在实际工程中得到了广泛应用。

递推最小二乘算法 篇4

伴随着多媒体和微电子技术的迅速发展,视频图像无线传输已成为移动业务新的增长亮点H264和MPEG-4等流行的国际视频压缩标准运用了块离散余弦变换(BDCT)、运动补偿、变字长编码等技术,能够在高压缩率下获得高质量图像。然而在易出错的无线信道中传输视频图像时,单个比特误码可能导致单个或连续宏块丢失的现象,严重影响了视频图像的传输质量。解码端的差错掩盖(errorconcealment)[1]技术是一种有效对抗信道干扰的方法。

差错掩盖是指解码器试图用主观可以接受的、近似原始质量的视频数据来遮盖差错受损数据,而不需要从编码器得到额外的信息。它利用视频图像在空间和时间上的高度相关性来掩盖差错带来的不良影响。当视频图像受损区域为帧内编码时,可用其周围的空域信息或频域信息进行插值掩盖当受损区域为帧间编码时,则主要采用时域掩盖。视频图像中的I帧采用帧内编码,其恢复图像的质量将影响以其为参考的其他帧,因此空间差错掩盖具有重要的意义。另一方面,在图像中物体运动剧烈区域,时间域差错掩盖效果不好,这时可以考虑用空间差错掩盖算法来恢复丢失像素点的信息。

由于已有成熟的丢失块检测方法,本文在丢失块确定的情况下,提出了针对I帧的基于最小二乘误差算法的空域掩盖算法。Xin.Li等人已经提出了从八个方向对丢失像素点预测,补偿的算法,但该算法是基于点序列预测的,对八个方向的最小二乘误差运算需要十分复杂的矩阵运算,因此该算法很难应用于实际图像掩盖当中本文首先用改进的算法,判断丢失像素点中是否含有边缘,对不包含在边缘中的像素点,本文认为是平稳点,采用双线形插值的算法进行恢复。对于包含在边缘中像素点,采用多方向的点序列预测算法,根据边缘方向,修正各方向的权重系数,最后根据各方向的预测值和权重,得到最终的恢复像素值。

1 改进的Sobel算法检测边缘

经典的Sobel算法的边缘检测,是在图像空间利用两个方向模板与图像进行邻域卷积来完成的。这两个模板一个是检测水平边缘,一个是检测垂直边缘。算法的基本原理是:由于图像边缘附近的亮度变化较大,所以可以把那些在邻域内,灰度变化超过某个适当阈值的像素点当作边缘点。

图像的边缘除了有垂直和水平两个方向以外,还有其他的边缘方向,为此,添加了其他六个方向的模板矩阵,如图1所示。

如图2所示,对每一个丢失像素点,利用Sobel算子对丢失块周围已经正确接收的像素点作边缘检测。将模板的中心像素与对应方向的像素点重合,把模板内的系数与其对应的图像像素值相乘作卷积,把得到的数值与预设的阈值Th作比较,如果大于该阈值Th,则认为存在该方向的边缘。对不含有边界的丢失像素点,本文采用边界插值的方式作掩盖。对在边界处的非平稳像素点,将对其采用多方向的最小二乘误差算法作预测,并根据它的边界情况,重新估算各方向的权重值,以求得到更准确的估计。

2 空域掩盖算法

为了提高掩盖算法的效率,对平稳像素点和非平稳像素点采用不同的掩盖方式。

2.1 边缘插值算法

当丢失像素点是平稳的,采用边缘插值的掩盖算法[4],把丢失像素点f(x,y)周围的块边界正常像素点的加权平均值来恢复。具体运算公式(1)如下:

如图3所示,其中d(i,j)※(x,y)为像素点f(x,y)到上,下,左,右的边界距离。

2.2 多方向的最小二乘(LS)掩盖算法

根据丢失像素点周围正确像素的分布情况,本文建立最小二乘预测模型,如图4所示。把丢失像素点周围正确的像素点作为输入,通过组成抽头输入u(i),u(i-1),…,u(i-M+1)与相应的抽头权重ω0,ω1,……,ωM-1之间的内积,并用d(i)作为期望响应,将估计误差或残差e(i)定义为期望响应d(i)和滤波器输出y(i)之间的差值,即e(i)=d(i)-y(i)

在本文中,把要修复像素点周围的正确数据看作是滤波器的输入,包括已经恢复出来的像素点,也认为是正确的数据滤波器的输出作为已修复的数据,根据最小二乘滤波器正则方程的矩阵形式[5],得到如下模型:

(1)在实验中,取M=7,把丢失像素点周围M×M范围内正确的像素点作为输入,包括已经恢复出来的数据。如图5所示,输入像素点的u(i),u(i-1),…,u(i-M+1)的M×M时间平均自相关矩阵为:

(2)抽头输入u(i),u(i-1),…,u(i-M+1)与期望响应d(i)之间的M×1时间平均互相关向量为:

(3)最小二乘滤波器的M×1抽头权向量为:

根据线性最小二乘滤波器正则方程的矩形形式,可得抽头权向量表达式为:

得到抽头系数以后,根据式(2)可以算出像素点的预测值。

最小二乘滤波器模型只针对平稳像素点可以作出准确的预测,由于边界情况的存在,如果只做单一方向的预测,将会产生较大的误差。因此,要对非平稳的像素点作八个方向的扫描和预测,通过加权平均的补偿方式得到最终的恢复值。算法进行过程中,(2)式是一个串行的扫描过程。在同一个丢失块中,后恢复的点由前面已恢复的点预测得出。实验证明,通过八方向的预测值的加权,可以得到更加准确的预测值。八个扫描顺序如图6所示。

根据像素点所处丢失像素点中的相对位置,赋予各方向的权重值。根据Xin.Li的算法,如图7所示,按照丢失像素点与周围正确点的相对位置,可以得出D0,D1,D2,D3三个数值,分别是丢失像素点跟上,下,左,右最近的正确像素点的距离。再根据式(8)和表1,可以得到各方向权重的初识值ωnk

按照通过Sobel算法得出来的像素点的边界情况,对各方向的权重作出修改,根据实验结果,具体如表2:

其中“++”和“+”代表增加该方向权重的程度,“--”和“-”代表减少该方向权重的程度。

最后,把经过修正的ωnk(k=1,2……8)代入式(9),得出最终恢复值Yn。

3 实验结果

为了验证本文算法的有效性,选择了Lena图和标准图像序列Football的I帧(第6帧和第30帧)进行试验。实验中,丢失像素点采用和两种,图8是离散丢失率23%的Lena图。图9是离散丢失率是22%的Football图像第6帧,图10是丢失率为10%的连续丢失图,为Football图像的第30帧。图10(b),图10(c),图10(d)是分别采用了线性插值算法,Xin.Li[2]算法以及本文所提出来的算法得到的结果。

由这些图中可以看出,无论是离散的,还是连续的块丢失,对比于传统的线性插值算法,Xin.Li的算法和本文的算法对丢失像素点的主边缘和纹理都有很好的恢复,平均有10%左右的提高。但Xin.Li的算法没有本文算法的效率高,也没有本文算法的效果好,如表3所示。Xin.Li[2]的算法由于要涉及大量的矩阵运算,运行时间要远远大于本文算法的时间,很难应用于实际的图像处理当中。本文算法考虑了边界因素,对处于平稳区域内的直接使用线性插值的算法,而对于非平稳区域,则采用最小二乘误差预测的办法,并根据边界因素,调整了各方向预测的权重系数,从而在降低了时间复杂度的同时,提高了图像的质量。

4 结束语

本文针对Xin.Li算法的不足,提出了一种有效的空域差错掩盖算法。算法根据丢失像素点中是否含有边界,把丢失像素点分为平稳像素点和非平稳像素点。对于不含有边缘的平稳像素点,直接采用线性插值的算法;对于含有边界的非平稳像素点,采用最小二乘算法,利用周围的像素点,建立最小二乘滤波器模型,然后通过多方向的扫描方式,根据边界情况,调整各方向的权重系数,最后用各方向预测值价权求和的方法,得出预测值。实验结果表明,本文算法能够获得优秀的图像质量,在降低了复杂度的同时适用于视频图像的实时传输

参考文献

[1]Wang Y,Zhu Q F,Shaw L.Maximally smooth image recovery in trans-form coding[J].IEEE Transactions on Communication,1993,41(10):1544-1551.

[2]Li Xin,Orchard M.Novel sequential error-concealment techniquesusing orientation adaptive interpolation[J].IEEE Transactions onCircuits and Systems for Video Technology,2002,12(10):857-864.

[3]Wilson Kwok,Huifang Sun.Multi-directional interpolation for spatialerror concealment[J].IEEE Trans.Consumer Electron 1993,39(3):455-460.

[4]Wang Y,Zhu Q F.Error control and concealment for video communica-tion:a review[J]//Proceedings of IEEE Multimedia Signal Processing,1998,5:974-997.

递推最小二乘算法 篇5

统计学习理论是由V.Vapnik等人建立的一种专门研究小样本情况下机器学习规律的理论,为解决有限样本学习问题提供了一个统一的框架,支持向量机就是在这一理论基础上发展而来的一种新的通用学习方法[1]。支持向量机通过结构风险最小化原理来提高泛化能力,较好地解决了小样本、非线性、高维数等实际问题,已在模式识别、信号处理、函数拟合等领域得到了成功应用,成为机器学习领域的研究热点之一。然而,用支持向量机求解分类问题,其本质上是一个二次规划问题,当样本数据量较大时,常规的数值优化算法及软件很难实现二次规划问题的求解。因此,运行时间和计算内存是支持向量机求解的主要瓶颈[2]。而作为支持向量机的改进算法,最小二乘支持向量机在运行速度上有了很大提高,同时减少了求解所需的计算资源,而其准确率并未明显下降,因此在模式识别和非线性函数拟合上得到了很好的应用。和其它算法一样,支持向量机的性能也依赖于模型的参数,研究人员对支持向量机的参数选择已作了很多研究,但至今还未提出明确的理论依据。如何实现模型的参数优化成为提高支持向量机学习性能和泛化能力的主要问题之一。本文结合J.A.K.Suyken等人开发的基于MATLAB语言的最小二乘支持向量机工具箱(LS-SVMlab),利用遗传算法实现了对支持向量机模型的参数优化,并基于一组木材表面颜色特征得到了很好的实验效果。

1 最小二乘支持向量机

支持向量机的基本思想是通过非线性变换将输入空间变换到一个高维空间,然后在这个新空间中求取最优线性分类超平面。支持向量机具有全局最优性和很强的泛化能力,能够很好地解决小样本学习问题,因此在模式识别、函数拟合等领域得到了广泛的应用。但同时,由于统计学习理论尚有很多问题需要进一步研究,在支持向量机结构、参数和核函数选择上还未有明确的理论指导。

最小二乘支持向量机LS-SVM(Least Squares Support Vector Machine)是在标准支持向量机上的一种扩展,由J.A.K.SuyKens和J.Vandewalle提出[3]。它采用最小二乘线性系统误差平方和作为损失函数,将求解过程变成了解一组等式方程,加快了求解速度,求解所需的计算资源较少,在模式识别和非线性函数拟合的应用中取得了很好的效果。

设训练样本D={ (xk,yk)| k=1,2,…,N},其中xk∈Rn为输入数据,yk∈R是输出类别。在权ω空间(原始空间)中的最小二乘支持向量机分类问题可以描述为求解以下问题:

$\underset{\omega ,b,e}{{\displaystyle \min }}\Phi (\omega ,b,e)=\frac{1}{2}\omega {}^{{\rm T}}\omega +\frac{1}{2}\gamma {\displaystyle \sum _{k=1}^{{\rm N}}e}{}_k^2(1)$

约束条件:yk[ωTφ(xk)+b]=1-ekk=1,...,N (2)

定义拉格朗日函数:

L(ω,b,e;α)=Φ(ω,b,e)-${\displaystyle \sum _{k=1}}$αk{yk[ωTφ(xk)+b]-1+ek} (3)

其中,拉格朗日乘子αk∈R。对上式进行优化:

$\{\begin{array}{l}\frac{\partial L}{\partial \omega }=0\to \omega ={\displaystyle \sum _{k=1}^{{\rm N}}\alpha }{}_{k}y{}_k\phi (x{}_k)\\ \frac{\partial L}{\partial b}=0\to {\displaystyle \sum _{k=1}^{{\rm N}}\alpha }{}_{k}y{}_k=0\\ \frac{\partial L}{\partial e{}_k}=0\to \alpha {}_k=\gamma e{}_k\\ \frac{\partial L}{\partial \alpha {}_k}=0\to y{}_k[\omega {}^{{\rm T}}\phi (x{}_k)+b]-1+e{}_k=0\end{array}(4)$

上式可化为求解下面的矩阵方程:

$\left[\begin{array}{cccc}{\rm I}& 0& 0& -{\rm Z}{}^{{\rm T}}\\ 0& 0& 0& -Y{}^{{\rm T}}\\ 0& 0& \gamma {\rm I}& -{\rm I}\\ {\rm Z}& Y& {\rm I}& 0\end{array}\right]\left[\begin{array}{l}\omega \\ b\\ e\\ \alpha \end{array}\right]=\left[\begin{array}{l}0\\ 0\\ 0\\ \rho {}_1\end{array}\right](5)$

即$\left[\begin{array}{cc}0& -Y{}^{{\rm T}}\\ Y& {\rm Z}{\rm Z}{}^{{\rm T}}+\gamma {}^{-1}{\rm I}\end{array}\right]\left[\begin{array}{l}b\\ \alpha \end{array}\right]=\left[\begin{array}{l}0\\ \rho {}_1\end{array}\right]$。

其中,Z=[φ(x1)Ty1,φ(x2)Ty2,…,φ(xN)TyN]T,Y=[y1,y2,…,yN]T,ρ1=[1,1,…,1]T,e=[ e1,e2,…,eN]T,α=[α1,α2,…,αN]T,同时将Mercer条件代入到Ω=ZZT,可得:

Ωkl=ykylφ(xk)Tφ(xl)=ykylΨ(xk,xl) (6)

因此,式(1)的分类问题通过式(5)和式(6)的线性问题得到最终解,而不是解二次规划问题。核函数可采用高斯径向基核函数、多项式核函数、多层感知器核函数和线性核函数等多种符合Mercer条件的核函数[3,4]。常用的核函数的表示形式如表1所示[5]。

最小二乘支持向量机将不等式约束转化为等式约束,将误差平方和损失函数作为训练集的经验损失,这样就将二次规划问题转化成线性方程组的求解,简化了计算复杂性,提高了支持向量机求解问题的速度和收敛精度。但同时也必须指出,由于每个样本数据对分类期都有贡献,LS-SVM失去了标准SVM稀疏性的优点。

2 最小二乘支持向量机的参数优化

2.1 支持向量机的参数选择问题

采用支持向量机求解模式识别或函数拟合问题需要选择一个核函数,由于不同的核函数构成的分类器的性能不同,因此选择哪种核函数对于支持向量机的设计是非常重要的。研究表明,惩罚因子和核函数的参数是影响支持向量机性能的主要原因,其中核函数的参数主要影响样本数据在高维特征空间中分布的复杂程度,而惩罚因子的作用是在确定的特征空间中调节支持向量机的置信范围和经验风险的比例[6]。所以,要想获得泛化能力良好的支持向量机,首先要选择合适的核函数参数将数据映射到合适的特征空间,然后针对该确定的特征空间寻找合适的惩罚因子以使支持向量机的置信范围和经验风险具有最佳比例。

作为支持向量机的一种扩展,最小二乘支持向量机同样存在参数选择问题。目前应用较多的是采用网格搜索和交叉验证相结合的支持向量机参数优化算法,该算法由初始给定的参数出发,在给定的参数范围内采用网格搜索算法和交叉验证寻找最优参数[7]。首先网格搜索算法选择要优化的参数对(γ,ρ),然后用交叉验证法对目标函数(如均方差最小)进行寻优,直至找到最佳的参数对,使交叉验证的精度最高。其中,网格搜索算法要将惩罚因子γ和核函数参数形成M×N个(γ,ρ)组合(其中M和N分别为γ和ρ的个数),分别训练不同的支持向量机,估计其学习精度,从而在M×N个(γ,ρ)的组合中得到学习精度最高的一个组合作为最优值。其计算量较大,尤其是训练样本集较大时搜索过程非常费时,因此这种方法在使用上受到一定程度的限制。在最优化方法众多的算法中,遗传算法以其强大的全局搜索能力、并行性、高效性的优点得到了广泛的应用。本文正是利用遗传算法的优点优化最小二乘支持向量机的参数,得到了较好的结果。

2.2 用遗传算法优化最小二乘支持向量机参数

遗传算法是一种借鉴生物界自然选择和自然遗传机制的随机搜索算法,能够在搜索过程中自动获取和积累有关搜索空间的知识,并自适应地控制搜索过程以求得最优解。它模拟自然选择和遗传中发生的复制、交叉和变异等现象,从任一初始种群出发,通过随机选择、交叉和变异操作,产生一群更适应环境的个体,使群体进化到搜索空间中越来越好的区域,这样一代一代不断繁衍下去,最后收敛到一群最适应环境的个体,求得问题的最优解[8]。

用遗传算法优化最小二乘支持向量机参数主要过程如下:

步骤1:设置初始值,如遗传算法的初始种群规模、最大遗传代数T、交叉概率、变异概率等。

步骤2:对要优化的参数根据其设定的范围进行二进制编码,随机产生初始种群。染色体为各参数二进制顺序排列组成,长度即为各参数二进制长度之和。设置遗传代数计数器t=0。

步骤3:计算种群中各个个体的适应度。这里将最小二乘支持向量机的分类正确率作为目标函数值,即个体的适应度,个体对应的参数的分类正确率越高,则该个体的适应度越大。

步骤4:根据个体适应度,按照一定规则(这里采用轮盘赌法)从当前种群中选出个体进入下一代。

步骤5:选择群体中的两个个体x1、x2作为父体以某个概率(交叉概率)进行交叉操作,产生两个新个体。这里采用单点交叉,交叉概率设为0.8。

步骤6:随机选择种群中的个体以一定的概率(变异概率)进行变异操作,通过随机改变个体中某些基因而产生新个体。变异概率设为0.05。

步骤7:终止条件判断。若t≤T,则转到步骤2;若t>T或平均适应度值变化持续小于某一常数超过一定代数,则所得到的具有最大适应的个体作为最优解输出,算法终止。

步骤8:对得到的最优解译码,得到优化的参数。

3 实验及结果

这里选择了东北常见的五种树种的木材样本,建立了包含1000幅图像的图像库。其中包括白桦、红松、落叶松、水曲柳和柞木五种树种的径切和弦切样本,共十类,每类各100个样本。提取了这些木材图像的19个L*a*b*颜色空间的颜色特征(包括L*、a*、b*分量、色相角Ag*、色饱和度C*、色差ΔE*、色相差ΔH*七个量的均值与标准差,以及前5个分量与其均值的差的均值),以此作为支持向量机的输入。选择每类样本中的30个作为参数优化及最终分类的训练样本,另外的30个作为参数优化时的测试样本,剩余的40个作为最终分类的测试样本。实验基于MATLAB R2007a实验平台,所用计算机配置为:P4 3.00 CPU、256M 内存、80G 硬盘、中文Windows XP 操作系统。

本文采用LS-SVMlab工具箱中提供的最小二乘支持向量机模型,将用遗传算法优化最小二乘支持向量机参数的方法与此工具箱中网格搜索优化参数的方法进行比较,每种方法分别运行10次,初始参数值均在设定范围内随机产生,最后将优化得到的参数构造支持向量机模型进行分类比较,比较结果如表2所示。

从表2的结果可以看出,用遗传算法进行优化得到的总体结果优于网格搜索的优化方法得到的结果,并且在运行时间上大大缩短了。由于篇幅有限,仅列出其它三种核函数参数优化的平均分类正确率和优化时间作为比较(10次),如表3所示。表3的结果中,采用线性核函数和多项式核函数用遗传算法进行参数优化在结果上优于网格搜索的优化方法得到的结果,而采用遗传算法也实现了对多层感知器核函数参数的优化,10次的平均分类正确率达到了80%以上。从各种核函数的结果对比也可以看出,采用不同的核函数对支持向量机分类器的分类性能也有很大差异,其中采用高斯径向基核函数和多层感知器核函数的结果都比较好,而采用线性核函数的分类结果并不理想(仅有70%左右),这也充分说明了核函数选择的重要性。

4 结 论

本文提出了一种采用遗传算法优化最小二乘支持向量机参数的方法。对于支持向量机参数的选择,现在也只能采用各种优化算法进行优化,而各种算法的寻优性能不同,采用遗传算法优化支持向量机参数的方法充分利用了遗传算法全局寻优和很强的搜索能力的优点,试验表明其在优化结果和时间上均优于采用网格搜索算法,取得了较为满意的效果。

参考文献

[1]张学工.关于统计学习理论与支持向量机[J].自动化学报,2000,26(1):32-42.

[2]杜树新,吴铁军.模式识别中的支持向量机方法[J].浙江大学学报:工学版,2003,37(5):521-527.

[3]Suyken J A K,Vandewalle J.Least squares support vector machineclassifiers[J].Neural Processing Letters,1999,9(3):293-300.

[4]朱家元,陈开陶,张恒喜.最小二乘支持向量机算法研究[J].计算机科学,2003,30(7):157-159.

[5]邓乃扬,田英杰.数据挖掘中的新方法——支持向量机[M].北京:科学出版社,2004.

[6]Vapnik V.Statistical Learning Theory[M].New York:Wiley-Inter-science,1998.

[7]王兴玲,李占斌.基于网格搜索的支持向量机核函数参数的确定[J].中国海洋大学学报,2005,35(5):859-862.

递推最小二乘算法 篇6

近年来基于导频的信道估计算法得到了广泛研究,R.Ne等[1]研究了慢时变信道环境下最优导频的放置方案,指出慢时变或非时变信道环境下,导频应该等间隔分布。良好的信道估计算法需要同时具有很好的估计精度和良好的跟踪能力。从实现准则来分,信道估计方法可分为最小均方误差(MMSE)[2,3]、最大似然估计(MLE)[4]、最小平方误差(LS)[5]等。MMSE算法[2]性能最好,可以得到精确的信道估计,但是MMSE算法需要知道信道的统计信息,并且具有很高的计算复杂度。线性最小均方误差(LMMSE)[3]信道估计器,是根据最佳低阶秩理论推导出的一种低复杂度的信道估计器,该估计器在信噪比(SNR)较低时性能受限,并且仍然需要已知信道的统计特性和信噪比。最大似然(ML)[4]信道估计算法,无须知道信道的频率相关特性,但是需要已知信道最大时延扩展和信噪比,且对有效信道长度具有敏感性。因此实际系统中常采用LS[5]算法,LS算法计算简单,不需要任何先验信道信息,但是LS算法的性能较差。为了提高LS算法的性能,[6]提出了基于离散傅立叶变换(DFT)的变换域信道估计算法,基于DFT的信道估计算法计算简单,不需要任何信道先验信息,性能优于LS算法。

本文首先分析了OFDM系统中子载波间干扰产生的原因和特性[7],接着分析了基于最小二乘LS的信道估计技术的性能,然后提出了同时考虑信道噪声和载波间干扰的OFDM信道估计算法。根据Tikhonov正则化原理[8],导出了总体最小二乘正则化解法的计算公式。该算法根据子载波间干扰特性,避免了一般总体最小二乘[9]方法带来的设计矩阵结构上的缺陷,并通过迭代计算进一步去除载波间干扰和噪声对系统性能的影响。理论分析和仿真结果表明该算法不仅可以有效的提高系统的误码率性能和信道估计的精度,而且对有限长度的CP仍然能保证通信质量,从而提高了频率利用率。

1 OFDM系统模型

1.1 信道模型

信道为瑞利多径衰落信道,且假设在一个OFDM符号间隔T内不变。在第n个OFDM符号间隔内,离散信道脉冲响应函数可表示为:

其中L为多径延迟扩展用采样时间归一化后的时间单位数,hl为第l条路径信号的衰落因子且为零均值的复高斯随机变量;τl第l条路径的多径时延,fDl为第l条路径的多谱勒频移,它会引起接收信号的载波间干扰(ICI);δ(t)为Kronecker delta函数。对于最大有效延迟数为L的信道,采用抽头系数为复系数的FIR滤波器作为信道的模型,其冲激响应的系数h=[h0,h1,…,hL]T。

1.2 OFDM信号模型

基于导频插入的典型的OFDM系统基带模型如图1所示,基于导频信道估计的OFDM系统首先在映射器中根据系统的调制方式对二进制信息进行符号映射分组并调制,串并变换后以固定周期在所有子载波上插入导频,形成频域序列X[k]。长度为N的序列X[k]被调制到工作子载波上,经过逆离散傅里叶变换(IFFT)模块被转化为时域序列x[n]。

为了消除码间干扰ISI,每帧OFDM符号后要插入保护间隔,保护间隔的长度Ng,Ng应大于等于信道冲击相应的最大延时,保护间隔可以由循环前缀组成,插入循环前缀后,此时OFDM系统的输出信号表示为:

发送信号通过带有加性白高斯噪声(AWGN)的频率选择性信道进行传输。在接收端,串并转换除去CP后,FFT之前的时域接收信号向量yt=[yi,0,yi,1,…,yi,N-1]T可以表示为:

式中G,A和B都是N×N矩阵,它们的组成如下面三个式子所示:

G是循环矩阵;A的最后g列为零。(4)式右边第二项和第三项分别代表ICI和ISI。如果g≥L,A和B为零矩阵,这时不存在ICI和ISI,均衡可以简化为简单的频域单抽头均衡。在CP不足时,残余ISI导致的当前符号内的ICI和前一符号带来的ISI会影响系统性能。这时,需要采取措施来消除这些干扰。

经过FFT变换,接收信号的频域表达式为:

式中Y(k)表示接收端第k个子载波的接收信号值,H(k)表示第k个子载波的信道频率特性值,w(k)表示频域AWGN,I(k)表示第k个子载波上干扰项。

2 最小二乘信道估计性能分析

LS信道估计是最简单的一种信道估计方法。设发送的导频信号X=[X0,X1,…,XP-1],从接收信号的导频位置获得导频符号Y=[Y0,Y1,…,YP-1],依据最小二乘准则:得到导频位置的信道特性:

在OFDM系统中,LS估计的MSE性能取决于导频所占用的频点及功率,当不考虑虚载波时,采用[10]中等间隔等功率的最优导频时,LS估计的MSE为:

式中σ2表示加性高斯白噪声信号功率,Ep为分配给导频的总功率。从式(4),(8)和(10)可以看到,LS信道估计结果仅考虑噪声而没有抑制干扰。如果能够采取一定的方法减小LS信道估计结果中噪声和干扰,那么LS信道估计的MSE可以得到降低。

3 正则总体最小二乘信道估计算法

针对上述ICI和ISI问题,接收端信号经过FFT后可以表示为:

式中W表示导频信号产生的数据矩阵,包括了I-CI和ISI干扰项后,不在为常数对角阵;E表示信道的不确定性对数据矩阵扰动;e表示噪声。

总体最小二乘模型的正则化准则[7]为:

其中α是正则化参数。

构造拉格朗日目标函数:

λ为NP×1拉格朗日乘算子矢量。对式(14)中Φ(E,e,λ)求偏导,并令偏导数等于0得:

求解上面(15)和(16)式得:

将(18)和(19)式代入约束方程得:

由(17)和(19)式得:

将(20)代入(21)化简得:

其中:

下面给出利用(22)、(23)式迭代计算的正则总体最小二乘信道估计算法,具体计算步骤如下:

步骤1用最小二乘算法估计信道的冲激响应,再用信道冲激响应的估计值分别构造出G,A和B;

步骤2初始化数据矩阵W和H(0),按照式(23)计算中间变量β;

步骤3按式(22)、(23)计算第1次迭代结果H(1)和β(1);

步骤4重复步骤3求第i次迭代结果:H(i)和β(i);

步骤5判断||H(i)-H(i-1)||<δ(δ为阈值),若成立则停止迭代;反之,重复步骤4。

上面迭代过程可以得到很精确的信道估计,其性能可由其均方误差来衡量:

4 性能仿真

为了评估文中提出的综合考虑了噪声和干扰的正则总体最小二乘RTLS信道估计算法的有效性,并与LS算法作比较。进行了系统的仿真,仿真参数如下:OFDM系统载频为2.0 GHz,子载波数N=128,16QAM调制,子载波间隔为:7.8125k Hz。各径时延依次为:0,2,4,8和12μs,延迟功率谱服从指数分布。先考察循环前缀充足CP=16时的仿真,即不存在存在I-CI和ISI干扰时改进算法有效性。下面给出仿真图如图2所示:

图2给出了两种算法的MSE随信噪比SNR变化的曲线。因为LS算法的性能受限于噪声,从图2可以看出,在低信噪比时,LS算法抑制噪声的能力比较差所以其MSE性能也差;而RTLS的性能有显著的提高。随着信噪比的增加,RTLS算法与LS算法性能越来越接近,因为随着噪声噪声功率的降低,噪声的影响越来越小。

系统BER随信噪比SNR变化的曲线如图3所示。从图3的BER比较图中也可得出与图2相同的结果,RTLS算法更能抑制噪声。在信噪比低于10d B时,RTLS算法有超过5d B的性能提高,信噪比越低性能提高越明显。

为了验证RTLS算法对干扰的稳健性,在两种循环前缀不足的情况下,系统BER随信噪比SNR变化的比较曲线如图4所示。从图4的比较结果可以看出,LS算法受ICI干扰影响较大,CP长度越小影响越大。因为RTLS总体考虑了噪声和干扰,ICI干扰对其影响不明显。随着SNR的增加,干扰是降低估计精度的主要原因,采用干扰抑制RTLS算法性能比LS算法有约6 d B的性能提高。理论上总体最小二乘算法比LS算法性能有10~15d B的性能提高[7],所以RTLS可大大提高估计的精度。

5 结束语

信道估计是实现OFDM系统性能的关键。由于总体最小二乘算法可以同时兼顾噪声和干扰,本文通过对频率选择性衰落信道下的OFDM系统ISI和ICI产生机理分析,避免了一般总体最小二乘算法设计矩阵结构上的缺陷。根据Tikhonov正则化原理,提出了总体最小二乘正则化信道估计算法。仿真结果都验证了该算法的有效性,尤其在存在干扰时改进明显,比普通最小二乘算法估计性能提高5d B以上。因此对噪声和干扰综合考虑可以提升信道估计算法鲁棒性。

参考文献

[1]NEGI R,CIOFFI J.Pilot tone selection for ehannel estimation in amobile OFDM system[J].IEEE Trans oll Consumer Electronics.1998,44(3):1122-1128.

[2]Coleri Sinem,Ergen Mustafa,Puf Anuj,Bahai Ahmad.Channelestimation techniques based on pilot arrangement in OFDM systems[J].IEEE Trans.on Broadcasting,2002,48(3):223-229.

[3]Y.Li,L.J.Cimini,and N.R.Sollenberger,Robust channel estimation forOFDM systems with rapid dispersive fading channels[J].IEEE Trans.Commun.,1998,46(7):902-915.

[4]Deneire L,Vandenameele P,Van der Perre L,et al.AlowcomplexityMLchannel estimator for OFDM[J].IEEE Transactions on Communications,2003,51(2):135-140.

[5]J.-C.Lin,Least-squarechannelestimationformobileOFDMcommunicationon time-varying frequency-selective fading channels[J].IEEE Trans.Veh.Technol.,2008,57(6):3538-3550.

[6]ZhaoY,HuangA.A novel channel estimation method for OFDMmobilecommunication systems based on pilot signals and transform-domainprocessing[C].ProceedingsofIEEE47thVehicularTechnologyConference.IEEE,1997:2089-2093.

[7]ZhongWei,MaoZhigang.EfficientTimeDomainResidualISICancellationfor OFDMBased WLANSystems[J].IEEE,2006,52(2):321-326.

[8]AmirBeck,AharonBen-Tal.OntheSolutionoftheTikhonovRegularizationoftheTotalLeastSquaresProblem[J].SIAMJ.OPTIM,2006,17(1):98-118.

[9]S.Chandrasekaran,MGu,A.h.Sayed,and K.E.Schubert.The degeneratebounded errors-in-variables model[J].SIAM J.Matrix.Anal.Appl.,2001,23(1):138-166.

递推最小二乘算法 篇7

目前,常规的电流、电压互感器的二次信号一般以额定电压AC 110 V、电流5 A的形式直接进入保护与测量装置。随着智能电网概念的提出,电流、电压互感器的直接数字化和网络化已成为必然趋势[1,2,3]。直接将二次信号转换为与一次信号成比例关系的数字量[2]是一种自然的做法。但这种方式的数字量频率高,信息量大,对接收端要求较高。如果对采样信息进行预处理,以电气相量的形式传输,则可以大大提高传输效率。

温漂是电子测量电路中都存在的普遍问题。在数字互感器中,温漂不仅会造成模拟运算放大器零点漂移,还会造成数字电路中关键部件——晶振频率的改变。前者直接叠加到A/D测量结果中,后者还影响采样周期,导致计算误差。而且,当温漂严重时还影响设备通信。

为此,本文针对交流信号直接采样样本,采用最小二乘法,采用优化滤波器计算工频相量,并通过波特率自适应法获得温漂程度,自动补偿温度造成的采样周期的偏移,并消除因温度造成的影响。

1 基于最小二乘法测量相量

设电气输入信号为:

$u(t)={\rm P}{}_0\text{e}{}^{-\lambda t}+{\displaystyle \sum _{k=1}^{{\rm M}}{\rm P}}{}_k\sin (k\omega t+\theta {}_k)(1)$

式中:P0为直流分量;λ为直流分量衰减时间常数;Pk为第k次谐波分量峰值,k=1,2,…,M;θk为第k次谐波相位角;ω=2πf,为基波角频率,f为电网频率。

对于测量与保护,最关心的是工频相量P1和θ1。

函数u(t)含有噪声时,与采样值u(k)之间必然存在误差e(k),根据函数理论,如果函数e2(t)在区间[0,T]内存在极大值或极小值,则函数在极值处的偏导数必然为0。经推导,此处即:

$u(t)=u(k)(2)$

式中:t∈[0,T],T为工频函数周期;k=0,1,…,N-1,为采样次数。

式(2)就是使误差最小的判据,简称最小二乘法判据。如果将式(1)中的P0e-λt按泰勒级数展开并取前2项,Pksin(kωt+θk)按三角函数展开,则式(1)可表示为如下形式:

$\begin{array}{l}u(t)={\rm P}{}_0-{\rm P}{}_0\lambda t+{\displaystyle \sum _{k=1}^{{\rm M}}\sin }k\omega t{\rm P}{}_k\cos \theta {}_k+\\ {\displaystyle \sum _{k=1}^{{\rm M}}\cos }k\omega t{\rm P}{}_k\sin \theta {}_k(3)\end{array}$

式(3)中,如果采样周期为τ,且τ=T/N,则sin kωt,cos kωt是可提前计算的常数,只有P0,λ,Pkcos θk,Pksin θk等是未知数。经过连续N次采样后,根据式(2),将得到N个方程,可表示为:

$\begin{array}{l}[\begin{array}{cccc}1& 1& \cdots & 1\\ \tau & 2\tau & \cdots & {\rm N}\tau \\ \sin \omega \tau & \sin \omega 2\tau & \cdots & \sin \omega {\rm N}\tau \\ \cos \omega \tau & \cos \omega 2\tau & \cdots & \cos \omega {\rm N}\tau \\ ⋮& ⋮& & ⋮\\ \sin {\rm M}\omega \tau & \sin {\rm M}\omega 2\tau & \cdots & \sin {\rm M}\omega {\rm N}\tau \\ \cos {\rm M}\omega \tau & \cos {\rm M}\omega 2\tau & \cdots & \cos {\rm M}\omega {\rm N}\tau \end{array}〗{}^{\text{Τ}}\cdot \\ [\begin{array}{c}{\rm P}{}_0\\ -\lambda {\rm P}{}_0\\ {\rm P}{}_1\cos \theta {}_1\\ {\rm P}{}_1\sin \theta {}_1\\ ⋮\\ {\rm P}{}_{{\rm M}}\cos \theta {}_{{\rm M}}\\ {\rm P}{}_{{\rm M}}\sin \theta {}_{{\rm M}}\end{array}〗=[\begin{array}{c}u(1)\\ u(2)\\ ⋮\\ u({\rm N})\end{array}〗(4)\end{array}$

如果将式(4)中系数矩阵用A表示,未知变量矩阵用X表示,采样矩阵用U表示,则在采样周期τ已定的情况下,A为N×2(M+1)阶的常数矩阵,X为2(M+1)×1阶的列矩阵,U为N×1阶的列矩阵。如果令2(M+1)=N,则A就是一个方阵,且最高可识别谐波倍数M=N/2-1,正好符合香农采样定理。根据矩阵理论,系数矩阵A存在逆矩阵A-1的条件为:|A|≠0,且AA-1=I。I为单位矩阵。为简便起见,此处不证明|A|≠0,只用AA-1=I来验证。

因此,滤波方程式(4)表示为:

$\mathit{X}=\mathit{A}{}^{-1}\mathit{U}(5)$

式中:A-1即为最小二乘滤波器。则

$\begin{array}{l}{\rm P}{}_1=\sqrt{({\rm P}{}_1\cos \theta {}_1){}^2+({\rm P}{}_1\sin \theta {}_1){}^2}(6)\\ \theta {}_1=\arctan \frac{{\rm P}{}_1\sin \theta {}_1}{{\rm P}{}_1\cos \theta {}_1}(7)\end{array}$

式中:P1cos θ1和P1sin θ1分别为列矩阵X的第3行和第4行元素,前者为基波相量的实部,后者为虚部。

从式(6)、式(7)可以看到,基波分量与直流无关,即A/D转换器的零漂及参考基准对交流量无影响。

2 优化的最小二乘滤波器

由于我们只关心工频成分,理论上每周期只需采样3次即可,但A-1的元素是带小数位的实数,且采样周期有偏差。如果每周期采样4次,则采样方程为:

$[\begin{array}{cccc}1& \tau & \sin \omega \tau & \cos \omega \tau \\ 1& 2\tau & \sin \omega 2\tau & \cos \omega 2\tau \\ 1& 3\tau & \sin \omega 3\tau & \cos \omega 3\tau \\ 1& 4\tau & \sin \omega 4\tau & \cos \omega 4\tau \end{array}$

$[\begin{array}{c}{\rm P}{}_0\\ -\lambda {\rm P}{}_0\\ {\rm P}{}_1\cos \theta {}_1\\ {\rm P}{}_1\sin \theta {}_1\end{array}$

$=[\begin{array}{c}u(1)\\ u(2)\\ u(3)\\ u(4)\end{array}〗(8)$

式中:τ=0.005 s;ω=2πf,f=50.0 Hz。

经过计算,得到的最小二乘滤波器A-1为:

$\mathit{A}{}^{-1}=[\begin{array}{cccc}\hfill 1.5& \hfill -1.0& \hfill 1.5& \hfill -1.0\\ \hfill -100.0& \hfill 100.0& \hfill -100.0& \hfill 100.0\\ \hfill 0.0& \hfill 0.5& \hfill -1.0& \hfill 0.5\\ \hfill 0.5& \hfill -1.0& \hfill 0.5& \hfill 0.0\end{array}〗(9)$

为了验证A-1是A的逆矩阵,可将τ=0.005 s,ω=2πf,f=50.0 Hz代入A中,可得:AA-1=I。因此,此处A-1为采样系数矩阵A的逆矩阵。

采用如式(9)所示的最小二乘滤波器,则采样时刻为第5 ms、第10 ms、第15 ms和第20 ms。在采用微处理器(MCU)的数字互感器中,采样时刻的定时一般来自于时钟中断,时钟是数字电路的基准,因此,时钟精度是一个必须要考虑的问题。众所周知,温度低则时钟慢,温度高则时钟快。因此,必须研究温度对本算法的影响。

3 温度对最小二乘法的影响

对于采用交流信号直接测量的算法,最小二乘法要求采样周期是固定的。而要实现准时采样,单片机必须通过定时器实现,而定时器的时钟基准来自单片机的晶振。一旦晶振因为温度原因而偏移,采样周期就发生了改变。为了研究采样周期对测量结果的影响,本文进行了仿真计算。

首先,设电网交流信号为u(t)=sin ωt=sin 2πft,f=50.0 Hz。仿真程序根据采样周期由-10%到+10%的变化量,分别根据式(6)、式(7)计算基波幅值和相位角,得到如表1和表2所示结果。

表1中的采样周期缩短,即采样频率加快,会导致采样提前结束,属欠采样。环境温度上升,晶振频率加快时即出现这一现象。表2中的采样周期加长,即采样速率变慢,将导致采样信号样本超过20 ms时间窗,属过采样。环境温度下降,晶振频率变慢时即出现这一现象。从表1和表2看,采样周期变化带来的计算误差有一定的规律,基本呈线性关系。因此,如果能得到采样周期的偏移量,可按表1、表2对计算结果进行修正。当然,如果采样周期能根据温度变化自动补偿,则计算误差自然不会产生。后者是最佳解决方案。下面提出解决这一问题的方法。

4 波特率自适应法校正测量误差

对MCU来说,如果没有直接测量温度或没有外部参考,是不可能得知温度变化的。但通信信道为此提供了可能。因为MCU的测量结果要通过通信的方式传送到控制中心,而正常工作情况下,控制中心一般处于室内,环境相对稳定。要保持设备的正常通信,两端的波特率必须完全一致,如果一端设备的波特率偏移严重,通信将不能正常进行。由温度引起的波特率偏移是一个缓慢的过程,伴随着波特率的偏移,通信误码率会逐渐增大,MCU通过监测误码率便可以得知温度发生了变化,然后修正波特率发生器的值,使误码率回归到正常范围。波特率发生器实际上是一个时钟基准。在MCU中,波特率发生器与晶振可以是1个,也可以是2个。以下只考虑同一个晶振情况下的解决办法。

4.1 波特率自适应法

对于具有通用异步收发器(UART)接口的MCU来说,错误信息有“帧错误”和“奇偶错误”。UART一帧信息一般可包括1个起始位、8个数据位、1个奇偶校验位、1个停止位,共11位。发送和接收端均根据自己的时间标尺(波特率周期)度量每位的时间长度,当本端时钟变慢时,时间标尺拉长。当一帧信息的第11位即停止位到来时,本端仲裁电路认为尚未到达,将会把已到达的停止位当做普通的第10位,虽然随后能采到停止位,但前10位的累积延时此时可能已导致了数据错位,现象表现为“奇偶错误”。

当本端时钟变快时,时间标尺变短,当一帧信息的第11位即停止位尚未到达时,本端仲裁电路认为时间到,将会把已到达的第10位误认为停止位。如果恰好第10位为“1”,则本端电路不报“帧错误”,但可能因累积产生了“奇偶错误”;但如果第10位为“0”,则立即产生“帧错误”,同时可能还伴随着“奇偶错误”。

以上是两端通信时钟失步时一帧信息传递过程中可能发生的现象。由于温度变化导致的这一现象是一个缓慢的过程,MCU可以监测到这一现象。由于波特率发生器的数值是可调的,MCU可以微调这个数值,与远端继续保持同步。按这个办法即实现了波特率的自动调整。

如果按这个微调量的百分比将交流采样的定时时钟同步调整,则采样时钟便得到了校正。既然采样时钟已经根据温度变化经过了调整,则自然不会产生因采样周期变化而造成的误差。因此,便通过波特率自适应实现了采样周期的自我校正,并因此保证了测量精度。

4.2 波特率自适应法校正测量误差举例

以PIC16F688单片机为例。设晶振频率fFosc为16 MHz,通信波特率BAUD0为38.4 kbit/s,采样周期为5 ms。以下步骤显示了通过波特率自适应法实现采样周期的自动调整进而避免测量误差的过程。

1)设定时器1(TIMER1)作为波特率发生器和采样时钟的共同基准(C语言格式):

T1CON=0x30;//设置滴答长度为2 μs

2)将波特率初始值(BAUD0)设定为38.4 kbit/s,将采样周期设定为5 ms(C语言格式):

void setbaud(void){

BRG0=25;//波特率寄存器初始值

SYNC=0;//异步模式

BRG16=0;//波特率控制寄存器8位模式

BRGH=1;//波特率加倍

SPBRG=BRG0;//BRG0=1000/38.4-1

}

void interrupt sample(void){//定时采样

TMR1IF=0;//清除TIMER1溢出标志

TMR1L=0x3c;//每5 ms产生一次溢出中断

TMR1H=0xf6;//0x10000-0xf63c=2 500

//2 500×2 μs=5 ms

AD_sample();//每5 ms启动一次A/D采样

Millisecond=millisecond+5;

if(millisecond>=1 000){//毫秒

Millisecond=0;

Second++;//秒

if (second>=60){

Second=0;

Minute++;}//分

}

if (FERR||RX9D) {ERROR++;}//帧或奇偶错误

}

3)MCU不断监测通信错误率,通过自动调整波特率获得波特率微调量kbaud_adj。

当大于门槛值,则kbaud_adj增加1或减去1,当ERROR回归正常范围后,波特率寄存器的值已变为kBRG1,其中:

KBRG1=KBRG0±kbaud_adj (10)

根据PIC16F688数据手册,调整后的波特率BAUD1(式中记为KBAUD1)为:

${\rm K}{}_{\text{B}\text{A}\text{U}\text{D}1}=\frac{f{}_{\text{F}\text{o}\text{s}\text{c}}}{16({\rm K}{}_{\text{B}\text{R}\text{G}1}+1)}(11)$

如果kbaud_adj≠0,则KBAUD1一定不等于KBAUD0。而且从式(11)还看到,波特率与晶振频率fFosc成正比关系。而此时控制中心的波特率并没有变化,新的波特率KBAUD1之所以不等于KBAUD0,是因为fFosc因温度发生了变化。

4)调整采样周期:

根据PIC16F688数据手册,将定时器1(TIMER1)的方式控制字T1CON设为“0x30”后,则TIMER1的时钟来源先将fFosc进行4分频,再8分频,分频系数为32,则一个滴答(tick)周期ttick为2 μs(ttick=32/fFosc=32/16=2)。

如果TIMER1产生5 ms的时钟中断,则TIMER1初始值为十六进制的“f63c”,产生一次溢出的间隔次数为2 500,乘以滴答周期ttick,正好为5 ms。即采样周期τ为:

τ=(65 536-tTIMER1)ttick (12)

由于晶振频率fFosc与波特率成正比,通过式(11)可以得到fFosc的偏移量为:

$f{}_{\text{o}\text{f}\text{f}}=\frac{{\rm K}{}_{\text{B}\text{A}\text{U}\text{D}1}-{\rm K}{}_{\text{B}\text{A}\text{U}\text{D}0}}{{\rm K}{}_{\text{B}\text{A}\text{U}\text{D}0}}\times 100％(13)$

如果将fFosc变化前后用0和1区分,则

fFosc1=fFosc0(1+foff) (14)

因此,随着温度变化,晶振频率fFosc也变化,如不修正必然产生误差。

要保证采样周期不变,必须按式(12)调整tTIMER1。如果tTIMER1的调整正好抵消fFosc的变化,则采样周期是精确的,也就无计算误差。

5 最小二乘滤波器对零漂的抑制

当A/D转换器因温度而发生零点漂移时,相当于在每个A/D转换结果后叠加了一个额外数值,即相当于每个采样样本u(k)叠加了一个直流成分Δ。如果采用本文的最小二乘滤波器,则此时的滤波器方程(5)成为如下形式:

$[{\rm P}{}_0-\lambda {\rm P}{}_0{\rm P}{}_1\cos \theta {}_1{\rm P}{}_1\sin \theta {}_1]{}^{\text{Τ}}=$

$[\begin{array}{c}u(1)+\Delta \\ u(2)+\Delta \\ u(3)+\Delta \\ u(4)+\Delta \end{array}$

设输入函数u(t)为式(13),当Δ按基波峰值的2%～10%变化时,根据式(11)和式(12)得到的仿真计算结果见表3。

表3说明最小二乘滤波器对A/D零漂具有完全的抑制作用。

6 实验结果

本互感器的设计精度为5P。实验按照电流互感器国标《GB 1208—2006》,针对互感器的电子测量部分进行了测试。升温速率为每5 min上升1 ℃,输入的二次电流为额定电流的40%,即2 A。恒温箱为GHX-800,电流源为继电保护测试仪TW30AE。通过后台主机阴极射线管(CRT)显示测量结果。主机与测量单元的连接采用RS-232光纤转换器,通信速率为38.4 kbit/s,通信周期为2 s。结果见表4。

从实验结果可以看到,温度达到50 ℃后,通信依然保持正常,且测量精度没有明显变化。电子测量单元的精度达到了1级。

7 结语

当数字电流、电压互感器采用最小二乘算法直接测量时,温度引起的晶振偏移对算法有影响,但通过本文所述的波特率自适应法可以修正这一因素,而A/D的零点漂移则对测量结果无任何影响。波特率自适应法是一种不需要直接测温而对温漂进行修正的间接方法,对其他的具有通信接口的数字化测量也有一定的参考价值。

参考文献

[1]袁宇波,卜强生,包玉树,等.电子式互感器数据分析系统的设计与应用.电力系统自动化,2009,33(20):78-82.YUAN Yubo,BU Qiangsheng,BAO Yushu,et al.Design andapplication of a data analyzing systemof electronic transformers.Automation of Electric Power Systems,2009,33(20):78-82.

[2]刘忠战,莫卫东.自励源电子式电流互感器研发中的关键技术.电力系统自动化,2009,33(6):67-69.LI U Zhongzhan,MO Weidong.The key technology of ECTwith self-excitation power supply.Automation of Electric PowerSystems,2009,33(6):67-69.