1 引　言

光斑质心估计是激光制导^[1-2]、光学探测^[3-4]、光电传感^[5-6]、计算机视觉等诸多领域的一项关键技术^[7-9]。许多的光斑质心估计算法被提出，如Fu等^[10]通过质心加权卡尔曼滤波器处理数字图像，有效地实现了视觉目标跟踪。质心加权的卡尔曼滤波算法是一种通过对点的灰度值进行加权来获得质心位置的方法。该算法仅对具有高信噪比、高对比度和均匀灰度分布的图像具有较高精度。因此，它很少用于高噪声场合。Bu等^[11]将传统的椭圆拟合方法发展为自适应圆−椭圆拟合方法，以提高树木直径的测量精度。采用自适应圆−椭圆拟合算法来估计质心位置，需要先检测边缘。因此，其最终的精度将受到边缘检测算法，如罗伯茨（Roberts）、索贝尔（Sobel）、普利维特（Prewitt）和拉普拉斯（Laplacian）的限制。边缘检测的误差会随着噪声的增加而增大，从而导致算法的精度下降。灰度重心法的基本原理是根据阈值分割出来的光斑灰度分布，按权重质心坐标作为跟踪点来进行激光光斑的定位，可以实现亚像素精度。在噪声存在的情况下，由于其对噪声敏感，该方法将导致平滑轮廓的测量结果显示为粗糙。

除上述的这几种方法外，还有基于数学形态学的质心估计方法^[12]、基于动态梯度的质心估计方法^[13]，以及基于边缘亚像素检测的质心估计方法^[14]等。

由此可见，各种算法在精确定位方面都存在一定问题。对于高精度挠度测量的场景，这些方法都不太适用。目前，基于函数拟合的方法，如点扩散函数（point diffusion function，PSF）模型拟合，可以实现精确定位。考虑到激光光斑的累积分布函数更接近高斯分布，本文提出了基于高斯拟合的光斑中心定位算法。自适应圆−椭圆拟合和质心加权卡尔曼滤波算法已在工程应用中广泛使用，尤其是在光学传感领域。因此，本文选择这两种算法与基于高斯拟合光斑中心提取方法进行比较。

1 研究背景介绍

要实现基于视觉的方法对混凝土结构挠度等参数的测量，首先要确定混凝土结构表面的“观察点”。对混凝土结构变化的测量都是以1个或几个观察点的位移来实现的。从数字图像处理的技术角度来看，混凝土表面的图像是一种非常特殊的图像形式，几乎没有显著的颜色特征，也没有典型的几何结构特征，只包含了一些不规则的简单纹理特征，如图1所示。在这种情况下，如何获得稳定可靠的“观察点”成为一个难题。

图 1. 混凝土表面图像

Fig. 1. Concrete surface image

下载图片 查看所有图片

在相关的研究工作中，有学者采用了一种折衷的办法：通过在混凝土表面设置“标靶”，从而回避了这一问题。采用的标靶形式大多为2种，如图2所示。

图 2. 观察点标靶形式

Fig. 2. Observation point target form

下载图片 查看所有图片

这样就把确定观察点简化为确定由特殊几何形状（直线、圆、矩形）构成的图形质心问题。但在检测前必须预设好标靶，而检测系统只能针对特定的标靶位置进行检测，这就在很大程度上牺牲了系统检测的方便性和灵活性。

采用激光指示的方式（如图3所示）可确定一个固定的“参考点”，从而实现对混凝土表面观测目标点相对于参考点的位移及形变的测量。

图 3. 激光指示光斑

Fig. 3. Laser indication spot

下载图片 查看所有图片

为了保证测量精度，本文将就激光指示“参考点”的精确、稳定定位展开研究。主要从图像分析的角度出发，通过图像处理技术提高定位精度，完成对参考点的精确定位。本文采用模拟的激光光斑，从噪声和对比度两方面比较了几种方法的定位精度。文中分析的仅仅是理想激光光斑，仅考虑噪声和对比度的影响。由于在实际场合中情况比较复杂，激光光斑受表面粗糙度和平整度，甚至环境的影响，实际光斑与理想的模拟光斑将有较大不同。另外，实际光斑的质量还与激光器有关。因此，采用理想的模拟光斑来定位，其结果将与实际情况存在差异。

2 基于高斯拟合的单激光光斑定位算法研究

2.1 基于高斯拟合光斑算法的推导

由于光学系统的衍射，光斑的边缘会产生模糊现象。在数学上，这个过程可表示为一个边缘函数ε（x）和一个点扩散函数h（x）的卷积^[15]。结果可表示为g（x）。g（x）是一个累积分布函数（cumulative distribution function，CDF），在应用中不够方便，需要运用其他简单且形状类似的函数来代替。由于激光光斑的CDF更趋近高斯分布^[16]，因此选择高斯函数作为拟合函数。

拟合函数的表达式为

$ I\left(x\text{，}y\right)=W\cdot \mathrm{exp}\left\{-\left[\frac{{\left(x-{x}_{0}\right)}^{2}}{{\sigma }_{1}^{2}}+\frac{{\left(y-{y}_{0}\right)}^{2}}{{\sigma }_{2}^{2}}\right]\right\} $ (1)

式中：I（x，y）为光斑在该截面（x，y）的光强； ${W}$为光强峰值； $ {\mathrm{\sigma }}_{1}，{\mathrm{\sigma }}_{2} $为2个方向上的标准差； $({x}_{0}，{y}_{0})$为光斑的中心位置。

光斑中心位置可由光强峰值位置来确定。本文就采用光强峰值所在的位置作为光斑中心位置。对上式两边取对数可得

$ z=a{x}^{2}+b{y}^{2}+cx+dy+f $ (2)

式中， $ z $为计算得到的像素灰度值的自然对数。

$ \left\{ \begin{gathered}a=\displaystyle-\frac{1}{{\sigma }_{1}^{2}} \\ b=\displaystyle-\frac{1}{{\sigma }_{2}^{2}} \\ c=\displaystyle\frac{2{x}_{0}}{{\sigma }_{1}^{2}} \\ d=\displaystyle\frac{2{y}_{0}}{{\sigma }_{2}^{2}} \\ f=\displaystyle {\rm{ln}}\;H-{x}_{0}^{2}/{\sigma}_{1}^{2}-{y}_{0}^{2}/{\sigma }_{2}^{2} \\ z={\rm{ln}}\;I(x,y) \\ \end{gathered}\right. $ (3)

取残差

$ {\varepsilon }_{i}=(a{x}_{i}^{\mathrm{'}2}+b{y}_{i}^{\mathrm{'}2}+c{x}_{i}^{\mathrm{'}}+d{y}_{i}^{\mathrm{'}}+f)-{z}_{i}^{\mathrm{'}} $ (4)

式中： ${(x}_{i}^{\mathrm{'}}，{y}_{i}^{\mathrm{'}})$∈ ${E}$，是构成光斑像素的坐标； ${z}_{i}^{\mathrm{'}}$为上述像素所在位置对应的像素强度的自然对数；a，b，c，d，f可通过最小平方法确定

$ \left[\begin{array}{ccccc}{x}_{1}^{\mathrm{'}2}& {y}_{1}^{\mathrm{'}2}& {x}_{1}^{\mathrm{'}}& {y}_{1}^{\mathrm{'}}& 1\\ {x}_{2}^{\mathrm{'}2}& {y}_{2}^{\mathrm{'}2}& {x}_{2}^{\mathrm{'}}& {y}_{2}^{\mathrm{'}}& 1\\ {x}_{3}^{\mathrm{'}2}& {y}_{3}^{\mathrm{'}2}& {x}_{3}^{\mathrm{'}}& {y}_{3}^{\mathrm{'}}& 1\\ \vdots & \vdots & \vdots&\vdots&\vdots \\ {x}_{n}^{\mathrm{'}2}& {y}_{n}^{\mathrm{'}2}& {x}_{n}^{\mathrm{'}}& {y}_{n}^{\mathrm{'}}& 1\end{array}\right]\left[\begin{array}{c}a\\ b\\ c\\ d\\ f\end{array}\right]=\left[\begin{array}{c}{z}_{1}^{\mathrm{'}}\\ {z}_{2}^{\mathrm{'}}\\ {z}_{3}^{\mathrm{'}}\\ \vdots \\ {z}_{n}^{\mathrm{'}}\end{array}\right] $ (5)

最小二乘拟合算法是获取参数的最常用方法之一。给定一组图像坐标 $\left( {{x_1}，{y_1}} \right)， \cdots ， \left( {{x_n}，{y_n}} \right)$，假定这些像素构成了光斑的轮廓，则可以使用最小二乘拟合算法来近似由这些像素拟合。该算法可以计算光斑的每个参数，且它的计算速度很快。

2.2 非线性参数的克拉美−罗下界证明

由于 $ z $是一个多元非线性函数，参数估计过程也必须是非线性的。在估计理论和统计中，克拉美−罗下界（Cramer-Rao lower bound，CRLB）^[17]定义了一个无偏估计的方差的下界。只有当CRLB收敛时， $ z $的估计才是可解的。接下来，将证明参数的CRLB是存在的。 $ {z}_{i,j} $表示计算得到的像素灰度值的自然对数， ${z}_{i,j}^{\mathrm{'}}$表示该点观察到的像素灰度值的自然对数。 ${p}=\{{a}，{b}， {c}， {d}， {f}\}$为未知参数。图像噪声被认为主要是加性高斯白噪声。

参数的似然函数为

$ L({z}^{\mathrm{'}},p)=\prod _{i\text{，}j}\frac{1}{\sqrt{2\text{π} {\sigma }^{2}}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{1}{2{\sigma }^{2}}{({z}_{i,j}^{\mathrm{'}}-{z}_{i,j})}^{2}\right] $ (6)

等式（6）两边同时取对数

$ \mathrm{l}\mathrm{n}({z}^{\mathrm{'}}\text{,}p)=-\frac{N}{2}\mathrm{l}\mathrm{n}\left(2\text{π} {\sigma }^{2}\right)-\frac{1}{2{\sigma }^{2}}\sum _{i,j}{({z}_{i,j}^{\mathrm{'}}-{z}_{i,j})}^{2} $ (7)

求方程（7）两边p的二阶导数

$ \begin{split} \displaystyle \frac{{\partial }^{2}\mathrm{ln}\;L \left({z}^{\mathrm{'}}\text{,}p\right)}{\partial {p}^{2}} = & \displaystyle \frac{1}{{ \displaystyle \sigma }^{2}}\sum _{i,j}\frac{{\partial }^{2}{z}_{i,j}}{\partial {p}^{2}}\left({z}_{i,j}^{\mathrm{'}} - {z}_{i,j}\right) - \\ &\frac{1}{{\sigma }^{2}} \sum _{i,j}{ \left(\frac{\partial {z}_{i,j}}{\partial p}\right)}^{2} \end{split} $ (8)

由于 $E\left({z}_{i,j}^{\mathrm{'}}\right)={z}_{i,j}$，式（8）的数学期望为

$ E\left[\frac{{\partial }^{2}\mathrm{l}\mathrm{n}\;L({z}^{\mathrm{'}}\text{,}p)}{\partial {p}^{2}}\right]=-\frac{1}{{\sigma }^{2}}\sum _{i,j}{\left(\frac{\partial {z}_{i,j}}{\partial p}\right)}^{2} $ (9)

对于满足公式（9）的任何无偏估计方差，其CRLB可定义为

$ \mathrm{v}\mathrm{a}\mathrm{r}\left(\hat{p}\right){\text{≥}} \displaystyle \frac{1}{-E\left[\displaystyle \frac{{\partial }^{2}\mathrm{l}\mathrm{n}\;L({z}^{\mathrm{'}}\text{,}p)}{\partial {p}^{2}}\right]} $ (10)

2.3 基于莱文贝格−马夸特的非线性参数优化

从CRLB的求解过程可以看出，所有未知参数的导数都是连续且存在的，并且拟合过程完全是非线性的。因此，基于莱文贝格−马夸特（Levenberg-Marquardt，LM）^[18]算法的最小二乘法被用来处理该拟合问题。计算值和采样值之间的残差为

$ r\left(p\right)=z-{z}_{i}^{\mathrm{'}} $ (11)

式中：z表示计算得到的像素灰度值的自然对数； ${z}_{i}^{\mathrm{'}}$为上述像素所在位置对应的像素强度的自然对数。本研究想要最小化 $ \parallel r\left(p\right)\parallel $，或者等价地找到

$ {p}^{\mathrm{*}}=\mathrm{a}\mathrm{r}\mathrm{g}\;\underset{p}{{\rm{min}}}\left\{R\right(p\left)\right\} $ (12)

式中， $R\left(p\right)=\displaystyle\frac{1}{2}{\sum }_{i}^{N}{\left[{r}_{i}\right(p\left)\right]}^{2}=\frac{1}{2}r\left(p\right)^{{\rm{T}}}r\left(p\right)$，式（12）的泰勒公式展开可表示为

$ r(p+h)=r\left(p\right)+{\boldsymbol{J}}\left(p\right)h+O(\parallel h{\parallel }^{2}) $ (13)

式中： ${\boldsymbol{J}}\in {{\boldsymbol{R}}}^{m\times n}$为雅可比矩阵； $ m $为采集数据的个数； $ n $为未知参数个数。这个雅可比矩阵展开式表示为：

$ {J}_{11}=\frac{\partial {r}_{1}}{\partial a}=-{x}^{2} $ (14)

$ {J}_{12}=\frac{\partial {r}_{1}}{\partial b}=-{y}^{2} $ (15)

$ {J}_{13}=\frac{\partial {r}_{1}}{\partial c}=-x $ (16)

$ {J}_{14}=\frac{\partial {r}_{1}}{\partial d}=-y $ (17)

$ {J}_{15}=\frac{\partial {r}_{1}}{\partial f}=-1 $ (18)

将式（13）代入 $ R\left(p\right) $

$ R(p+h)\simeq L\left(h\right)=R\left(p\right)+{h}^{{\rm{T}}}{{\boldsymbol{J}}}^{{\rm{T}}}r+0.5{h}^{{\rm{T}}}{{\boldsymbol{J}}}^{{\rm{T}}}{\boldsymbol{J}}h $ (19)

$ L\left(h\right) $的梯度和黑森量分别为

$ {L}^{\mathrm{'}}\left(h\right)={{\boldsymbol{J}}}^{{\rm{T}}}r+{{\boldsymbol{J}}}^{{\rm{T}}}{\boldsymbol{J}}h\text{，}{L}^{\mathrm{'}\mathrm{'}}\left(h\right)={{\boldsymbol{J}}}^{{\rm{T}}}{\boldsymbol{J}} $ (20)

${{L}}^{\mathrm{'}\mathrm{'}}\left({h}\right)$独立于h。如果J是满秩，它是对称的。换句话说， ${R}({p}+{h})$将有q个极小值，可以通过求解如下式子得到

$ h=-{\left({{\boldsymbol{J}}}^{{\rm{T}}}{\boldsymbol{J}}\right)}^{-1}{{\boldsymbol{J}}}^{{\rm{T}}}r $ (21)

通过确定迭代收敛方向的阻尼参数μ并将其代入等式（21），得到LM方法的迭代步长

$ h=-{({{\boldsymbol{J}}}^{{\rm{T}}}{\boldsymbol{J}}+\mu {\boldsymbol{I}})}^{-1}{{\boldsymbol{J}}}^{{\rm{T}}}r $ (22)

式中，I是与 ${{\boldsymbol{J}}}^{{\rm{T}}}{\boldsymbol{J}}$秩相同的单位矩阵。

停止标准是

$ (\parallel h\parallel{\text{≤}} {\varepsilon }_{1})\;{{\rm{and}}}\;[\parallel {r}_{{\rm{new}}}-r\parallel{\text{≤}} (\parallel r\parallel +{\varepsilon }_{2}\left)\right] $ (23)

式中， $ {\varepsilon }_{1} $和 $ {\varepsilon }_{2} $为小正数，由用户本身选择。

3 算法性能研究

3.1 实验设备及假定

实验过程中所用实验平台的参数见表1。

在本节中，测试了所研究的方法在模拟图像及真实图像上的应用结果，并与文献中一些主要方法的研究结果进行了比较。测试中，仅限于搜索那些圆心位于图像内的圆（如果估计的圆心位于图像外，则立即丢弃该圆）。如果像素小于30（T_min=30），则停止圆检测任务。讨论时，仅考虑半径大于5个像素且被检测到的圆周长大于40％的内部像素数的圆。

3.2 算法性能测试

为了更好地说明算法的性能，对3种算法进行了噪声测试，结果如图4所示。

图 4. 3种算法噪声鲁棒性测试结果

Fig. 4. Noise robustness test results of the three algorithms

下载图片 查看所有图片

在上述的测试中，将所研究的方法与自适应圆–椭圆拟合和质心加权卡尔曼滤波算法进行了比较。为了解决模拟问题，首先需要人工生成一个具有固定能量的光斑点，该光斑为半径是50 pixel的圆。

图4结果显示，在3种算法中光斑中心提取的均方根误差（root mean square error，RMSE）^[19]均随着信噪比的增加呈指数递减。这说明，各算法的提取精度均有明显提高。在信噪比为40 dB时，高斯拟合光斑定位算法中心提取的RMSE几乎衰减为0，而自适应圆−椭圆拟合算法的RMSE为50 dB，质心加权卡尔曼滤波算法的RMSE为60 dB。

3.3 算法鲁棒性测试

考虑到真实环境中存在噪声，在模拟对比度实验中，将图像的信噪比设置为50 dB，进行了噪声鲁棒性的测试，仿真结果如图5所示。

图 5. 3种算法的低对比度测试结果

Fig. 5. Low contrast test results of the three algorithms

下载图片 查看所有图片

光斑中心提取的RMSE随着对比度的降低而增大。也就是说，随着对比度的降低，3种方法的提取精度会变差。在对比度降至原图像的25%之前，本文算法的鲁棒性优于其他算法。但当对比度继续下降时，这一优势将消失。

为了评估原始图像受噪声直接影响时对整体工作流程的影响，进行了1组参数化仿真。图6展示了高斯白噪声作用下的测试性能。通过使用MATLAB函数中的Imnoise 添加零均值噪声，方差范围从0.01（1％附加噪声）到0.3（30％附加噪声）。噪声的增减会导致检测到的圆的数量减少，特别是对于沿圆周具有较低梯度的那些圆。值得注意的是，当附加噪声不超过27%时，测得的错误光斑数为0；当附加噪声接近30%时，将检测到（平均）1个错误光斑。

图 6. 模拟光斑图像中检测到的对象

Fig. 6. Simulation of an object detected in a spot image

下载图片 查看所有图片

4 结 论

本研究初步验证了将激光指示作为“参考点”的挠度测量的可行性。利用模拟的激光光斑，从噪声和对比度的影响来对比了几种方法的定位精度。实验结果表明，依据点扩散函数的模型拟合可以实现精确定位，说明激光光斑的累积分布函数可以较好地采用高斯拟合的方法进行中心定位。文中分析的仅仅是模拟光斑，也只是考虑了噪声和对比度的影响。在实际场合中，由于激光光斑受表面粗糙度、投射角度、距离远近、平整度等环境的影响较大，因此，采用理想的模拟光斑进行定位的计算结果将会与实际情况存在差异。后续需对该问题作进一步深入研究。