双目视觉因模拟人眼完成对客观三维场景的还原，已经在机器人智能视觉、三维测量、产品零件检测、医学外科手术导航等领域得到广泛应用^[1-3]。由于其具有结构简单、操作简便、测量范围大、灵活性强等优点，目前已成为三维测量领域中的研究热点。所涉及的关键技术包括摄像机标定、图像预处理、特征点识别、立体匹配、三维重建等^[4]。但是，随着研究的不断深入，双目视觉的局限性也逐渐突显出来。首先，测量精度不高，当测量环境复杂或者被测工件表面形态复杂时，特征点的识别和匹配就会面临很大的挑战；其次，计算速度慢，很难保障实时性，双目视觉三维测量技术以像素为单位进行处理，计算量庞大，对于复杂三维形貌的重建更是耗时长，难以实现实时测量；最后，系统对硬件要求高，由于双目视觉以像素为单位进行处理，所以测量精度在很大程度上依赖于硬件的水平。

为解决这些问题，本文提出一种三维测量方案，用光敏位置探测器(position sensitive detector，PSD)来替代传统双目视觉中所采用的电荷耦合元件(charge coupled device，CCD)相机，通过2个不同视角的PSD对激光点的实时捕捉来实现三维测量，从而省去了传统双目视觉技术中的特征点识别和立体匹配环节，极大地简化了计算模型，为实时测量提供了可能。由于与传统双目视觉在测量形式上产生了根本变化，已有的标定方案不再适用，所以需要探索针对双PSD测量系统的标定方案。

1 PSD的工作原理

PSD是一种对其光敏面上的入射光斑重心位置敏感的光电器件，其测量原理基于横向光电效应，具有分辨率高、响应速度快、传输速度快等优点^[5-7]。

本研究采用二维PSD，其结构如图 1所示。当激光器照向PSD光敏面时，会产生光电流，由于接触层是一层均匀的电阻层，所以被两端电极收集到的光电流将与入射点到电极的距离成反比。由此，入射光斑P点与4个电极上的电流的关系如式(1)所示，其参考坐标系以PSD光敏面的中心点为坐标原点。

$ \left\{ \begin{array}{l} u = \frac{L}{2} \cdot \frac{{\left( {{I_1} + {I_4}} \right) - \left( {{I_2} + {I_3}} \right)}}{{{I_1} + {I_2} + {I_3} + {I_4}}},\\ v = \frac{L}{2} \cdot \frac{{\left( {{I_3} + {I_4}} \right) - \left( {{I_1} + {I_2}} \right)}}{{{I_1} + {I_2} + {I_3} + {I_4}}}. \end{array} \right. $

(1)

其中：L为光敏面边长, I₁、I₂、I₃、I₄分别为4个电极上的光电流。

2 系统测量原理

本研究提出的基于双PSD的视觉系统，就是通过左右2个PSD对同一个激光点的测量数据来还原激光点的真实坐标，即还原工件的三维信息^[8]。在测量过程中，2个PSD同一时刻测得的光敏面坐标值都与该激光点本身的世界坐标一一对应，这种映射关系是由系统本身唯一确定的。因此，系统的标定工作就是建立起左右PSD的光敏面坐标和激光点的世界坐标之间的转换模型，即获得系统的结构参数。理想的双PSD视觉系统成像模型如图 2所示。

以左PSD为例，o-uv为左PSD光敏面坐标系，O_c-X_cY_cZ_c为左镜头坐标系，平面X_cO_cY_c平行于平面uov，且左光轴Z_c交平面uov于(u₀, v₀)。O_w-X_wY_wZ_w为世界坐标系。激光点P在左PSD光敏面上的坐标为P₁。

由此，P点左PSD光敏面坐标与世界坐标之间的关系如式(2)所示，该表达式其实就是基于小孔成像原理的线性映射关系^[9]。

$ {Z_c}\left[ {\begin{array}{*{20}{c}} u\\ v\\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{f_x}}&0&{{u_0}}&0\\ 0&{{f_y}}&{{v_0}}&0\\ 0&0&1&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{R}}&\mathit{\boldsymbol{T}}\\ 0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{X_w}}\\ {{Y_w}}\\ {{Z_w}}\\ 1 \end{array}} \right]. $

(2)

其中：(X_w, Y_w, Z_w)为激光点P点的世界坐标，(u, v)为P点在PSD光敏面上的坐标，f_x、f_y为镜头的焦距，(u₀, v₀)为镜头光轴在PSD光敏面上的初始坐标值(这4个参数即为系统的内参)，R、T分别为世界坐标系到镜头坐标系的旋转矩阵和平移矩阵(即系统的外参)。

3 标定方案 3.1 方案1：改进的Faugeras标定+ LM优化

本方案先利用改进的Faugeras标定方法求出系统的内外参数初值，再对该初值进行LM算法优化。

1) 求取内外参数初值。

令

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} f&0&{{u_0}}\\ 0&f&{{v_0}}\\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{R}}&\mathit{\boldsymbol{T}} \end{array}} \right] = }\\ {\left[ {\begin{array}{*{20}{c}} {{m_{11}}}&{{m_{12}}}&{{m_{13}}}&{{m_{14}}}\\ {{m_{21}}}&{{m_{22}}}&{{m_{23}}}&{{m_{24}}}\\ {{m_{31}}}&{{m_{32}}}&{{m_{33}}}&{{m_{34}}} \end{array}} \right].} \end{array} $

(3)

将式(3)展开变换可得：

$ \left\{ \begin{array}{l} {m_{11}}{x_w} + {m_{12}}{y_w} + {m_{13}}{z_w} + {m_{14}} - {m_{31}}{x_w}u - \\ \;\;\;{m_{32}}{y_w}u - {m_{33}}{z_w}u = {m_{34}}u,\\ {m_{21}}{x_w} + {m_{22}}{y_w} + {m_{23}}{z_w} + {m_{24}} - {m_{31}}{x_w}v - \\ \;\;\;\;{m_{32}}{y_w}v - {m_{33}}{z_w}v = {m_{34}}v. \end{array} \right. $

(4)

若已知n个空间点的世界坐标，代入式(4)，则能得到2n个方程组

$ {\mathit{\boldsymbol{A}}_1}{\mathit{\boldsymbol{X}}_1} + {\mathit{\boldsymbol{A}}_2}{\mathit{\boldsymbol{X}}_2} = 0. $

(5)

其中：

$ {\mathit{\boldsymbol{A}}_1} = \left[ {\begin{array}{*{20}{c}} {{x_{w1}}}&{{y_{w1}}}&{{z_{w1}}}&1&0&0&0&0&{ - {u_1}}\\ 0&0&0&0&{{x_{w1}}}&{{y_{w1}}}&{{z_{w1}}}&1&{ - {v_1}}\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots \\ {{x_{wn}}}&{{y_{wn}}}&{{z_{wn}}}&1&0&0&0&0&{ - {u_n}}\\ 0&0&0&0&{{x_{wn}}}&{{y_{wn}}}&{{z_{wn}}}&1&{ - {v_n}} \end{array}} \right], $

$ {\mathit{\boldsymbol{A}}_2} = \left[ {\begin{array}{*{20}{c}} { - {u_1}{x_{w1}}}&{ - {u_1}{y_{w1}}}&{ - {u_1}{z_{w1}}}\\ { - {v_1}{x_{w1}}}&{ - {v_1}{y_{w1}}}&{ - {v_1}{z_{w1}}}\\ \vdots&\vdots&\vdots \\ { - {u_n}{x_{wn}}}&{ - {u_n}{y_{wn}}}&{ - {u_n}{z_{wn}}}\\ { - {v_n}{x_{wn}}}&{ - {v_n}{y_{wn}}}&{ - {v_n}{z_{wn}}} \end{array}} \right], $

$ {\mathit{\boldsymbol{X}}_1} = {\left[ {\begin{array}{*{20}{c}} {{m_{11}}}&{{m_{12}}}&{{m_{13}}}&{{m_{14}}}&{{m_{21}}}&{{m_{22}}}&{{m_{23}}}&{{m_{24}}}&{{m_{34}}} \end{array}} \right]^{\rm{T}}}, $

$ {\mathit{\boldsymbol{X}}_2} = {\left[ {\begin{array}{*{20}{c}} {{m_{31}}}&{{m_{32}}}&{{m_{33}}} \end{array}} \right]^{\rm{T}}}. $

由于R为正交矩阵，则‖X₂‖=1，由此构造性能指标函数：

$ J = {\left\| {{\mathit{\boldsymbol{A}}_1}{\mathit{\boldsymbol{X}}_1} + {\mathit{\boldsymbol{A}}_2}{\mathit{\boldsymbol{X}}_2}} \right\|^2} + \lambda \left( {{{\left\| {{\mathit{\boldsymbol{X}}_2}} \right\|}^2} - 1} \right). $

(6)

其中：λ为所有实数。为使J最小，将J分别对X₁、X₂求偏导数，并令其等于0，则有：

$ \left\{ \begin{array}{l} \frac{{\partial J}}{{\partial {\mathit{\boldsymbol{X}}_1}}} = 2\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_1}{\mathit{\boldsymbol{X}}_1} + 2\mathit{\boldsymbol{A}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}_2}{\mathit{\boldsymbol{X}}_2} = 0,\\ \frac{{\partial J}}{{\partial {\mathit{\boldsymbol{X}}_2}}} = 2\mathit{\boldsymbol{A}}_2^{\rm{T}}{\mathit{\boldsymbol{A}}_2}{\mathit{\boldsymbol{X}}_2} + 2\mathit{\boldsymbol{A}}_2^{\rm{T}}{\mathit{\boldsymbol{A}}_1}{\mathit{\boldsymbol{X}}_1} + 2\lambda {\mathit{\boldsymbol{X}}_2} = 0. \end{array} \right. $

(7)

由此可求出X₁、X₂，即可得到模型的内外参数。此方法更好地保证了旋转矩阵R的正交性，比线性标定法精度更高^[10]。

2) LM优化。

假设需要优化的模型为：

$ {y_i} = {f_i}\left( {y,\theta } \right) + {e_i},i = 1,2, \cdots ,n. $

(8)

其中：f_i(y, θ)为式(2)所反映的函数关系，θ=[θ₁, θ₂, …, θ_m]^T为待优化的系统内外参数，{e_i}为白噪声序列，且n>m。

则最小二乘性能指标函数为：

$ \begin{array}{*{20}{c}} {J = \sum\limits_{i = 1}^n {{{\left[ {{y_i} - {f_i}\left( {y,\theta } \right)} \right]}^2}} = }\\ {{{\left[ {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{F}}\left( \mathit{\boldsymbol{\theta }} \right)} \right]}^{\rm{T}}}\left[ {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{F}}\left( \mathit{\boldsymbol{\theta }} \right)} \right].} \end{array} $

(9)

取极小，则令$\frac{{\partial J}}{{\partial \mathit{\boldsymbol{\theta }}}} = 0 $，得到：

$ \frac{{\partial J}}{{\partial \mathit{\boldsymbol{\theta }}}} = 2\frac{{\partial {{\left[ {\mathit{\boldsymbol{Y}} - F\left( \mathit{\boldsymbol{\theta }} \right)} \right]}^{\rm{T}}}}}{{\partial \mathit{\boldsymbol{\theta }}}}\left[ {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{F}}\left( \mathit{\boldsymbol{\theta }} \right)} \right] = 0. $

(10)

将F(θ)在θ_k附近用Taylor级数展开，可得到：

$ \mathit{\boldsymbol{F}}\left( \mathit{\boldsymbol{\theta }} \right) \approx \mathit{\boldsymbol{F}}\left( {{\mathit{\boldsymbol{\theta }}_k}} \right) + \frac{{\partial \mathit{\boldsymbol{F}}\left( \mathit{\boldsymbol{\theta }} \right)}}{{\partial \mathit{\boldsymbol{\theta }}}}\left| {_{\mathit{\boldsymbol{\theta }} = {\mathit{\boldsymbol{\theta }}_k}}} \right.\Delta \mathit{\boldsymbol{\theta }} = \mathit{\boldsymbol{F}}\left( {{\mathit{\boldsymbol{\theta }}_k}} \right) + g\left( {{\mathit{\boldsymbol{\theta }}_k}} \right)\Delta \mathit{\boldsymbol{\theta }}. $

(11)

其中：

$ g\left( {{\mathit{\boldsymbol{\theta }}_k}} \right) = \frac{{\partial \mathit{\boldsymbol{F}}\left( \mathit{\boldsymbol{\theta }} \right)}}{{\partial \mathit{\boldsymbol{\theta }}}}\left| {_{\mathit{\boldsymbol{\theta }} = {\mathit{\boldsymbol{\theta }}_k}}} \right.. $

代入到式(10)，可得

$ \Delta \mathit{\boldsymbol{\theta }} = - \frac{1}{2}{\left[ {g{{\left( {{\mathit{\boldsymbol{\theta }}_k}} \right)}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{\theta }}_k}} \right)} \right]^{ - 1}}\frac{{\partial J}}{{\partial \mathit{\boldsymbol{\theta }}}}\left| {_{\mathit{\boldsymbol{\theta }} = {\mathit{\boldsymbol{\theta }}_k}}} \right.. $

(12)

上述方法为Gauss-Newton法，虽然它的收敛速度很快，但无法满足Jacobi矩阵奇异的情况。本研究采用的LM算法在此基础之上，引入阻尼因子μ_k，解决了Jacobi矩阵奇异时的收敛问题^[11-12]。

则有：

$ \Delta \mathit{\boldsymbol{\theta }} = - \frac{1}{2}{\left[ {g{{\left( {{\mathit{\boldsymbol{\theta }}_k}} \right)}^{\rm{T}}}g\left( {{\theta _k}} \right) + {\mu _k}\mathit{\boldsymbol{I}}} \right]^{ - 1}}\frac{{\partial J}}{{\partial \mathit{\boldsymbol{\theta }}}}\left| {_{\mathit{\boldsymbol{\theta }} = {\mathit{\boldsymbol{\theta }}_k}}} \right.. $

(13)

其中I为单位矩阵。LM算法在一般条件下都能保证收敛到极小值。

3.2 方案2：BP神经网络

一个典型的三层BP神经网络的基本结构如图 3所示，它由输入层、输出层和隐层组成。理论上，一个简单的三层BP神经网络在隐层的神经元个数可以任意设定的前提下，能够以任意精度逼近任意连续函数^[13]。BP神经网络的学习算法是基于最小二乘法的基本思想，通过梯度下降法使得网络的均方误差达到最小值。学习过程主要分为正向传播和反向权值调整两个阶段^[14]。

假设神经网络输入层、隐层、输出层的神经元个数分别为L、M、N，训练总的样本数为Q，输入层神经元i与隐层神经元j之间的连接权值为ω_ij，隐层神经元k与隐层神经元t之间的连接权值为ω_kt。

1) 正向传播。

假设单个输入样本为p，p=[x₁, x₂, …, x_L]，则隐层第j个神经元的输入为

$ {\rm{ne}}{{\rm{t}}_j} = \sum\limits_{i = 1}^L {{\omega _{ij}}{x_i}} - {\theta _j}. $

(14)

输出为

$ {y_j} = f\left( {{\rm{ne}}{{\rm{t}}_j}} \right) = f\left( {\sum\limits_{i = 1}^L {{\omega _{ij}}{x_i}} - {\theta _j}} \right). $

(15)

其中：θ_j为阈值, j=1, 2, …, M.

输出层第t个神经元的输入为：

$ {\rm{ne}}{{\rm{t}}_t} = \sum\limits_{k = 1}^M {{\omega _{kt}}{y_k}} - {\theta _t}. $

(16)

输出为

$ {O_t} = f\left( {{\rm{ne}}{{\rm{t}}_t}} \right) = f\left( {\sum\limits_{k = 1}^M {{\omega _{kt}}{y_k}} - {\theta _t}} \right). $

(17)

其中t=1, 2, …, N.

2) 权值反向修正。

假设期望输出为q，则单个输入样本p的误差为：

$ {J_p} = \frac{1}{2}\sum\limits_{t = 1}^N {{{\left( {{q_t} - {O_t}} \right)}^2}} . $

(18)

整个样本的误差为：

$ {J_1} = \frac{1}{2}\sum\limits_{l = 1}^Q {\sum\limits_{t = 1}^N {{{\left( {q_t^l - O_t^l} \right)}^2}} } . $

(19)

BP算法采用梯度下降法来调整神经网络的权值使得误差达到最小值，权值的修正方向沿着J₁梯度的方向：

$ \begin{array}{*{20}{c}} {\Delta {\omega _{kt}} = - \eta \frac{{\partial {J_P}}}{{\partial {\omega _{kt}}}} = - \eta \frac{{\partial {J_P}}}{{\partial {O_t}}}\frac{{\partial {O_t}}}{{\partial {\omega _{ij}}}} = }\\ {\eta \sum\limits_{t = 1}^N {\left( {{q_t} - {O_t}} \right){{f'}_2}\left( {{\rm{ne}}{{\rm{t}}_t}} \right){y_k}} ,} \end{array} $

(20)

$ \begin{array}{*{20}{c}} {\Delta {\omega _{ij}} = - \eta \frac{{\partial {J_P}}}{{\partial {\omega _{ij}}}} = - \eta \frac{{\partial {J_P}}}{{\partial {O_t}}}\frac{{\partial {O_t}}}{{\partial {y_j}}}\frac{{\partial {y_j}}}{{\partial {\omega _{ij}}}} = }\\ {\eta \sum\limits_{t = 1}^N {\left( {{q_t} - {O_t}} \right){{f'}_2}\left( {{\rm{ne}}{{\rm{t}}_t}} \right){\omega _{kt}}{{f'}_1}\left( {{\rm{ne}}{{\rm{t}}_j}} \right){x_i}} .} \end{array} $

(21)

其中：η为学习速率，f₁为隐层的传递函数，f₂为输出层的传递函数。

结合式(19)—(21)，可以推导得出各层权值的修正公式：

$ \left\{ \begin{array}{l} {\omega _{ij}}\left( {z + 1} \right) = {\omega _{ij}}\left( z \right) + \Delta {\omega _{ij}},\\ {\omega _{kt}}\left( {z + 1} \right) = {\omega _{kt}}\left( z \right) + \Delta {\omega _{kt}}. \end{array} \right. $

(22)

在整个神经网络的训练过程中，权值不断地修正使得误差函数沿负梯度方向下降，当整个样本误差函数的值满足初始设定的精度要求或者达到最大循环次数时修正停止，即训练结束^[15]。

4 实验验证

双PSD视觉系统的视野范围选定在60 mm×60 mm×40 mm范围内。PSD选用光敏面大小为9 mm ×9 mm的二维PSD探测器。镜头选用固定焦距工业镜头，根据所测工件的大小及镜头视野范围选定镜头焦距为12 mm。PSD和镜头通过C镜头座连接。数据采集选用USB数据采集卡，8路差分输入。激光器采用二极管点状可控激光器，波长为660 nm。最终搭建出的系统如图 4所示。

标定过程利用三维导轨的移动来获取激光点的实际三维坐标，采用镜头移动，激光器相对固定的标定方式，可以有效地避免因平台表面不平整所带来的误差。激光点在空间40 mm×40 mm×20 mm的范围内，每隔10 mm进行采集，最终得到75个空间点在左右2个PSD上的检测数据。

4.1 改进的Faugeras标定+LM优化

利用采集到的75组数据，采用本研究提出的方案1对系统进行标定，最终标定结果内部参数如表 1所示，外部参数如表 2所示。

重新采集Z=0平面上169个激光点的数据，利用标定出的参数还原激光点的世界坐标，得到169组误差值，分别对坐标x、y、z方向上的误差R_x、R_y、R_z进行曲面拟合，得到的曲面如图 5所示。

从图 5可以看出，已标定区域40 mm×40 mm的误差相对较小且方差小，而未标定的区域则误差较大且方差大。

4.2 BP神经网络

BP神经网络的输入层神经元个数为4个，分别是左右2个PSD的靶面坐标，输出层的神经元个数为3个，即世界坐标。训练样本采用与方案1同样的75组数据。经实验验证，选择隐层神经元8个，训练函数选取Bayesian Regularization函数，训练效果最佳。同理，得到Z=0平面上的误差拟合曲面如图 6所示。

从图 6可以看出，采用神经网络标定法得到的模型，x、y方向上的误差较为理想，而z方向上的误差较大。

5 结论

本文提出一种基于双PSD的三维测量方法，利用2个PSD从不同的角度捕获、跟踪激光点的方式来还原工件的三维信息，省去了传统双目视觉中的特征点识别和匹配部分，极大简化了计算模型。本文提出的2种标定方案，改进的Faugeras标定+LM优化和BP神经网络方法均取得了较好的精度，其中：2种方法都得到了较好的误差特性；BP神经网络标定法建模简单、快速，而且在y、z方向上误差特性均优于改进的Faugeras标定+LM优化方法；改进的Faugeras标定+LM优化方法在误差纠正之后，得到的精度更高，且更稳定，更能符合系统高精度的要求。