近年来，计算机视觉技术发展极为迅速，被大量的应用到工业测试领域。在大型挤压机活动部件五自由度实时监测系统^[1]中，稳定准确的求取激光束中心是其中一个关键步骤^[2]，寻求具有强鲁棒性和高滤波精度的滤波方法是实现这一目标的重要途径。

卡尔曼滤波算法是线性高斯系统的最优滤波器^[3-6]，它具有完备的理论基础和简便的算法，在计算机视觉领域有着广泛的应用^[7]。标准的卡尔曼滤波鲁棒性差，滤波精度取决于能否建立系统精确的数学模型和获得噪声的准确统计特性^[8]，在工业在线检测中，这些通常难以做到，需采用自适应滤波来代替常规卡尔曼滤波。目前常用的自适应滤波方法有强跟踪卡尔曼滤波、基于极大似然准则的自适应卡尔曼滤波、Sage-Husa自适应卡尔曼滤波等^[9-11]。强跟踪滤波器采用多重渐消因子削弱以往数据对模型误差的影响，提高了算法准确跟踪状态突变的能力；基于极大似然准则的自适应滤波器从系统量测量出现概率最大的角度利用滤波过程中获取的残差对系统过程噪声和观测噪声进行实时的估计和调整，从而提高滤波算法的精确度和鲁棒性；Sage-Husa自适应卡尔曼滤波也是一种在线估计过程噪声和观测噪声滤波算法，相对于基于极大似然准则的自适应滤波而言，改进在于引入了遗忘因子，进而增加了较新的观测值的权重，相应减小了较陈旧数据对参数估计的影响。这些算法中，Sage-Husa自适应卡尔曼滤波实用性较强，在工程实际中最为常用，但因为计算量较大，用于工业在线检测时，无法保证系统的实时性。文献[12]提到了一种简化的Sage-Husa算法，该算法假设过程噪声和观测噪声均是无偏的，且只对观测噪声进行在线估计，该算法在保证滤波精度的前提下减少了计算量，但仍需每次滤波都对噪声进行重新估计，也没有考虑野值的判定和处理，仍不能满足工业在线检测的要求。

针对上述问题，文中在国内外相关研究的基础上^[13-15]，对简化的Sage-Husa算法进行了进一步改进，该方法先对检测环境是否发生变化和新的观测值是否为野值做出判断，并分别采取不同的对策，且不需要每次都对观测噪声进行估计，算法的精度、鲁棒性和实用性比简化的Sage-Husa波算法更高，能满足工业在线检测的要求。

1 Sage-Husa算法及其简化

由卡尔曼滤波原理可知，线性离散系统的状态转移方程和观测方程分别为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{X}}\left( k \right) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {k,k - 1} \right)\mathit{\boldsymbol{X}}\left( {k - 1} \right) + }\\ {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left( {k,k - 1} \right)\mathit{\boldsymbol{W}}\left( {k - 1} \right),} \end{array} $

(1)

$ \mathit{\boldsymbol{Z}}\left( k \right) = \mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{X}}\left( k \right) + \mathit{\boldsymbol{V}}\left( k \right), $

(2)

式中：X(k)为系统状态向量；Z(k)为观测向量；Φ(k，k-1)为系统状态转移矩阵；Γ(k，k-1)为系统噪声矩阵；W(k)为过程噪声向量；V(k)为观测噪声向量。线性高斯系统的过程噪声W(k)和观测噪声V(k)可处理成具有时变均值和时变协方差阵的不相关高斯随机过程，即

$ E\left( {\mathit{\boldsymbol{W}}\left( k \right)} \right) = \mathit{\boldsymbol{q}}\left( k \right), $

(3)

$ E\left\{ {\left[ {\mathit{\boldsymbol{W}}\left( k \right) - \mathit{\boldsymbol{q}}\left( k \right)} \right]{{\left[ {\mathit{\boldsymbol{W}}\left( j \right) - \mathit{\boldsymbol{q}}\left( j \right)} \right]}^{\rm{T}}}} \right\} = \mathit{\boldsymbol{Q}}\left( k \right){\delta _{kj}}, $

(4)

$ \mathit{E}\left( {\mathit{\boldsymbol{V}}\left( k \right)} \right) = \mathit{\boldsymbol{r}}\left( k \right), $

(5)

$ E\left\{ {\left[ {\mathit{\boldsymbol{V}}\left( k \right) - \mathit{\boldsymbol{r}}\left( k \right)} \right]{{\left[ {\mathit{\boldsymbol{V}}\left( j \right) - \mathit{\boldsymbol{r}}\left( j \right)} \right]}^{\rm{T}}}} \right\} = \mathit{\boldsymbol{R}}\left( k \right){\delta _{kj}}, $

(6)

式中δ_kj为克罗内克函数。Sage-Husa自适应滤波算法是在噪声均值q(k)、r(k)协方差阵Q(k)、R(k)都未知的情况下，基于观测值来求取噪声统计特性的卡尔曼滤波器。标准的Sage-Husa算法由于计算较繁琐，实际应用过程中，观测噪声协方差R(k)滤波的影响更重要一些，通常可只认为R(k)为未知，且假设过程噪声和观测噪声均是无偏的，则Sage-Husa算法可简化为

$ \mathit{\boldsymbol{X}}\left( k \right) = \mathit{\boldsymbol{X}}\left( {k\left| {k - 1} \right.} \right) + \mathit{\boldsymbol{K}}\left( k \right)v\left( k \right), $

(7)

$ \mathit{\boldsymbol{X}}\left( {k\left| {k - 1} \right.} \right) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {k,k - 1} \right)\mathit{\boldsymbol{X}}\left( {k - 1} \right), $

(8)

$ \mathit{\boldsymbol{v}}\left( k \right) = \mathit{\boldsymbol{Z}}\left( k \right) - \mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{X}}\left( {k\left| {k - 1} \right.} \right), $

(9)

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}}\left( k \right) = \mathit{\boldsymbol{p}}\left( {k\left| {k - 1} \right.} \right){\mathit{\boldsymbol{H}}^{\rm{T}}}\left( k \right)\left[ {\mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{p}}\left( {k\left| {k - } \right.} \right.} \right.}\\ {{{\left. {\left. 1 \right){\mathit{\boldsymbol{H}}^{\rm{T}}}\left( k \right) + \mathit{\boldsymbol{R}}\left( k \right)} \right]}^{ - 1}},} \end{array} $

(10)

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{p}}\left( {k\left| {k - 1} \right.} \right) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {k\left| {k - 1} \right.} \right)p\left( {k - 1} \right){\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}\left( {k,k - 1} \right) + }\\ {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left( {k,k - 1} \right)\mathit{\boldsymbol{Q}}\left( {k - 1} \right){\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^{\rm{T}}}\left( {k,k - 1} \right)} \end{array} $

(11)

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{p}}\left( k \right) = \left[ {1 - \mathit{\boldsymbol{K}}\left( k \right)\mathit{\boldsymbol{H}}\left( k \right)} \right]\mathit{\boldsymbol{p}}\left( {k\left| {k - 1} \right.} \right)\left[ {1 - } \right.}\\ {{{\left. {\mathit{\boldsymbol{K}}\left( k \right)\mathit{\boldsymbol{H}}\left( k \right)} \right]}^{\rm{T}}} + \mathit{\boldsymbol{K}}\left( k \right)\mathit{\boldsymbol{R}}\left( {k - 1} \right){\mathit{\boldsymbol{K}}^{\rm{T}}}\left( k \right),} \end{array} $

(12)

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{R}}\left( k \right) = \left( {1 - {d_k}} \right)\mathit{\boldsymbol{R}}\left( {k - 1} \right) + {d_k}\left\{ {\left[ {1 - } \right.} \right.}\\ {\left. {\mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{K}}\left( {k - 1} \right)} \right]\mathit{\boldsymbol{v}}\left( k \right){\mathit{\boldsymbol{v}}^{\rm{T}}}\left( k \right){{\left[ {1 - \mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{K}}\left( {k - 1} \right)} \right]}^{\rm{T}}} + }\\ {\left. {\mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{p}}\left( {k - 1} \right){\mathit{\boldsymbol{H}}^{\rm{T}}}\left( k \right)} \right\}。} \end{array} $

(13)

式(13)即观测噪声估计方程，其中d_k=(1-b)/(1-b^k+1)，b为遗忘因子，取值范围通常为[0.95，0.99]。

2 简化Sage-Husa算法的改进

简化的Sage-Husa自适应滤波算法能在系统数学模型和噪声的统计特性不够精确的情况下达到较好的滤波效果，但是将其应用于工业环境下激光束中心定位时仍然存在明显的不足：简化的算法仍然需要每次对观测噪声协方差进行计算，降低了算法的实用性；对于工业环境下的在线检测，原始检测数据通常存在一定比例的野值，在使用自适应滤波算法时，如果不进行甄别，并将野值排除于滤波记忆之外，则反而可能造成滤波结果不稳定以及精度下降。

针对这些问题，引入一种改进的Sage-Husa算法，该方法先通过滤波发散判据对检测环境是否发生变化和新的观测值是否为野值做出判断，并分别采取不同的对策，从而提高算法的实用性、鲁棒性和精度。

2.1 滤波异常判断

滤波状态可通过实际余项与理论余项是否相符来判断，滤波发散判据为^[16]

$ {\mathit{\boldsymbol{v}}^{\rm{T}}}\left( k \right)\mathit{\boldsymbol{v}}\left( k \right) > \gamma tr\left( {\left[ {E\left( {\mathit{\boldsymbol{v}}\left( k \right){\mathit{\boldsymbol{v}}^{\rm{T}}}\left( k \right)} \right.} \right]} \right), $

(14)

式中：γ为储备系数，γ>1；tr为矩阵的迹；v(k)即kalman滤波器的新息向量。

由kalman滤波原理有

$ \begin{array}{*{20}{c}} {E\left( {\mathit{\boldsymbol{v}}\left( k \right){\mathit{\boldsymbol{v}}^{\rm{T}}}\left( k \right)} \right) = \mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{p}}\left( {k\left| {k - } \right.} \right.}\\ {\left. 1 \right){\mathit{\boldsymbol{H}}^{\rm{T}}}\left( k \right) + \mathit{\boldsymbol{R}}\left( k \right),} \end{array} $

(15)

把式(15)带入式(14)，滤波发散判据变为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{v}}\left( k \right){\mathit{\boldsymbol{v}}^{\rm{T}}}\left( k \right) > \gamma tr\left[ {\mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{p}}\left( {k\left| {k - } \right.} \right.} \right.}\\ {\left. {\left. 1 \right){\mathit{\boldsymbol{H}}^{\rm{T}}}\left( k \right) + \mathit{\boldsymbol{R}}\left( k \right)} \right],} \end{array} $

(16)

其中，v(k)在为获得一个新的观测向量后，按照式(9)计算。这样，滤波器工作是否正常可叙为以下假设检验问题

$ {\mathit{\boldsymbol{H}}_0}:\mathit{\boldsymbol{v}}\left( k \right){\mathit{\boldsymbol{v}}^{\rm{T}}}\left( k \right) < \gamma tr\left[ {\mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{p}}\left( {k\left| {k - 1} \right.} \right){\mathit{\boldsymbol{H}}^{\rm{T}}}\left( k \right) + \mathit{\boldsymbol{R}}\left( k \right)} \right]; $

$ {\mathit{\boldsymbol{H}}_1}:\mathit{\boldsymbol{v}}\left( k \right){\mathit{\boldsymbol{v}}^{\rm{T}}}\left( k \right) \ge \gamma tr\left[ {\mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{p}}\left( {k\left| {k - 1} \right.} \right){\mathit{\boldsymbol{H}}^{\rm{T}}}\left( k \right) + \mathit{\boldsymbol{R}}\left( k \right)} \right]。$

当H₀成立时，则说明滤波器工作正常，可让滤波器继续工作，而当H₁成立时则说明滤波器工作异常，需要进一步判断。

2.2 测量野值判断

前面已对滤波器工作是否正常做出了判断，而若滤波器工作异常，则还存在两种可能：新的观测值为野值；观测环境已经发生了变化，需对观测噪声R(k)进行重新估算。

由随机过程相关性理论可知，平稳自相关随机序列的每一步的取值虽然是随机的，但它的取值范围受限于自身序列的自相关性。基于这一原理，研究的方法在判定滤波器工作异常的情况下，通过一次延时判定来确定滤波异常的原因。

当2.1小节中H₁成立时，先不做如何处理，直接跳过k时刻的状态，依据k-1时刻的状态对k+1时刻进行估计，有

$ \mathit{\boldsymbol{X}}\left( {k + 1\left| {k - 1} \right.} \right) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {k + 1,k - 1} \right)\mathit{\boldsymbol{X}}\left( {k - 1} \right), $

(17)

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{v}}\left( {k + 1} \right) = \mathit{\boldsymbol{Z}}\left( {k + 1} \right) - \mathit{\boldsymbol{H}}\left( {k + 1} \right)}\\ {\mathit{\boldsymbol{X}}\left( {K + 1\left| {k - 1} \right.} \right)。} \end{array} $

(18)

此时，理论上有

$ \begin{array}{*{20}{c}} {E\left( {\mathit{\boldsymbol{v}}\left( {k + 1} \right){\mathit{\boldsymbol{v}}^{\rm{T}}}\left( {k + 1} \right)} \right) = }\\ {\mathit{\boldsymbol{H}}\left( {k + 1} \right)\mathit{\boldsymbol{p}}\left( {k + 1\left| {k - 1} \right.} \right){\mathit{\boldsymbol{H}}^{\rm{T}}}\left( {k + 1} \right) + \mathit{\boldsymbol{R}}\left( {k + 1} \right)。} \end{array} $

(19)

跳过k时刻状态后，k-1时刻与k+1时刻直接联系时，滤波器是否工作正常可叙为以下假设检验问题：

$ \begin{array}{l} {{\mathit{\boldsymbol{H'}}}_0}:\gamma tr\left[ {\mathit{\boldsymbol{v}}\left( {k + 1} \right){\mathit{\boldsymbol{v}}^{\rm{T}}}\left( {k + 1} \right) < \mathit{\boldsymbol{H}}\left( {k + 1} \right)\mathit{\boldsymbol{p}}\left( {\left. {k + 1} \right|} \right.} \right.\\ \left. {\left. {k - 1} \right){\mathit{\boldsymbol{H}}^{\rm{T}}}\left( {k + 1} \right) + \mathit{\boldsymbol{R}}\left( {k + 1} \right)} \right], \end{array} $

$ \begin{array}{l} {{\mathit{\boldsymbol{H'}}}_1}:\gamma tr\left[ {\mathit{\boldsymbol{v}}\left( {k + 1} \right){\mathit{\boldsymbol{v}}^{\rm{T}}}\left( {k + 1} \right) \ge \mathit{\boldsymbol{H}}\left( {k + 1} \right)\mathit{\boldsymbol{p}}\left( {\left. {k + 1} \right|} \right.} \right.\\ \left. {\left. {k - 1} \right){\mathit{\boldsymbol{H}}^{\rm{T}}}\left( {k + 1} \right) + \mathit{\boldsymbol{R}}\left( {k + 1} \right)} \right]。\end{array} $

当H′₀成立时，则说明跳过k时刻状态后，滤波器工作正常，k时刻状态判断为测量野值，直接用估计值代替k时刻状态，不把该观测数值引入滤波器，从而抑制野值对滤波精度的影响；而当H′₁成立时则说明对于新的观测值Z(k)和Z(k+1)，滤波器均工作异常，观测环境已发生变化，需要在k时刻即对R(k)进行重新计算。

3 实验及分析

下面通过模拟一种工业环境下计算机视觉激光束中心检测的典型工况来验证研究算法的有效性。

3.1 实验系统

实验系统由发射装置^[17]、运动平台和接收装置组成，如图 1所示。激光发射器为准直扩束激光器，波长650 nm、功率20 mw；接收装置中图像传感器分辨率为320×240像素，工业镜头为8 mm固定焦距镜头。图像处理及滤波算法是在PC机(主频2.1 G双核处理器、内存2G)上通过Visual C++ 6.0编程实现的。

3.2 实验过程

实验通过对光束相对位移检测结果滤波来验证算法的有效性。检测过程为：把装置第一次采集到的光斑中心位置设为基准参考点，通过运动平台把接收装置在水平和竖直方向各移动1 mm采集一次图像，通过中心定位算法对光斑相对位移进行测定，并对检测结果进行滤波。实验工况为：前20帧检测图像为检测装置状态良好的情况下采集的，后20帧检测图像为激光接收装置存在严重漏光的情况下采集的，并引入密度为0.01的椒盐噪声，另外，第7帧、第19帧以及第35帧图像引入密度为0.03的椒盐噪声模拟产生野值的环境。实验中采用的检测图像如图 2所示。

实验过程中，通过物/像比例尺标定^[18]得到单像素尺寸为0.343 1 mm×0.343 1 mm。正常的检测图像中心提取算法为曲线拟合法，干扰环境下检测图像和野值图像中心提取算法为Hough圆变换。经多次测定，在平面坐标系的x和y方向，正常检测图像的观测噪声幅值为0.12 mm，协方差约为0.006；干扰环境下检测图像观测噪声幅值为0.45 mm，协方差约为0.014；野值图像的测量偏差大于0.7 mm。过程噪声幅值为0.1 mm，协方差约为0.01，储备系数γ取值2.2。

滤波过程中，系统状态变量选为光斑的中心位置坐标(x、y)和速度$\dot x$、${\dot y}$则系统的状态转移方程为

$ \left[ {\begin{array}{*{20}{c}} {x\left( k \right)}\\ {\dot x\left( k \right)}\\ {y\left( k \right)}\\ {\dot y\left( k \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&T&0&0\\ 0&1&0&0\\ 0&0&1&T\\ 0&0&0&1 \end{array}} \right] \times \left[ {\begin{array}{*{20}{c}} {x\left( {k - 1} \right)}\\ {\dot x\left( {k - 1} \right)}\\ {y\left( {k - 1} \right)}\\ {\dot y\left( {k - 1} \right)} \end{array}} \right] + \mathit{\boldsymbol{Q}}\left( k \right), $

(20)

式中，T为每两帧图像间时间差，这里，可统一把T归为1 s。直接观测量为光斑的中心位置坐标(x、y)，则系统的观测方程为

$ \left[ {\begin{array}{*{20}{c}} {x\left( k \right)}\\ {y\left( k \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&0&0&0\\ 0&0&1&0 \end{array}} \right] \times \left[ {\begin{array}{*{20}{c}} {x\left( k \right)}\\ {\dot x\left( k \right)}\\ {y\left( k \right)}\\ {\dot y\left( k \right)} \end{array}} \right] + \mathit{\boldsymbol{R}}\left( k \right)。$

(21)

3.3 实验结果及分析

研究的滤波方法工作过程中，按照式(14)计算得到的实际余项与理论余项比的变化曲线如图 3所示，其中，第7帧、第19帧以及第35帧图像，由于野值的出现，实际余项与理论余项比有明显的增加，而按照(2.2)小节所述方法，直接跳过这些状态，联系之前与之后状态，则该异常情况消失，可据此判断这些观测值为野值。第22帧图像实际余项与理论余项比也有明显增加，这是因为第21帧图像开始的观测环境的变化这时在测量数据中得到了体现，此时，跳过这一状态，实际余项与理论余项仍然不符，需采用Sage-Husa算法对系统观测噪声进行重新估计。该组数据表明，文中方法采用的判据能准确的判断观测环境是否变化以及新观测值是否为野值。

激光中心定位测量误差、简化的Sage-Husa算法滤波后的定位误差以及文中方法滤波后的定位误差如图 4所示。

由图可知，2种滤波算法均能减小中心定位误差，且都能适应观测环境的变换，2种方法均是有效的。正常情况下文中方法和Sage-Husa滤波算法滤波精度相当，而一旦出现测量野值，算法能对其进行辨别，并排除在滤波记忆之外，因此滤波精度更高；另外，由于文中方法不需要每次对观测噪声进行估计，算法的实用性比简化的Sage-Husa波算法更好。

4 结论

1) 针对工业在线检测的特点及要求，对简化的Sage-Husa自适应滤波v算法进行了进一步改进，该方法先通过滤波发散判据判断滤波器工作状态是否异常，并进一步找出滤波异常的原因，分别采取不同的对策，从而提高滤波算法在工业检测中的性能。

2) 提出了优化算法在模拟工业环境下激光束中心定位中的应用，实验表明，该方法能对观测环境的变化和观测野值做出准确判断，与传统方法相比，精度较高、实用性好、鲁棒性强，适用于工业在线检测的随机误差滤波。