Design and application of distributed Raman thermometry system based on multi-modal Gaussian fitting
-
摘要:
针对温感区域长度小于空间分辨率导致温度测量不准确的问题，提出一种多点温度校正方案，并介绍了一种可应用于中小尺度剖面温度分布测量的装置和方法。采用基于多峰高斯拟合的方案，实现了多个测量不准确温度点的同时校正，提高了系统的温度准确度，设计了分布式拉曼测温系统样机，通过将光纤部分缠绕在PPR管上实现了中小尺度剖面温度分布测量。实验装置使2 m和1.5 m温感区域的测量准确度提高到1 ℃；测得冬季黄河万家寨水库冰盖厚度约为42.33 cm，与人工测量结果相差约为1.67 cm；测得36 m深河水剖面的温度约在0.32 ℃～0.90 ℃范围内波动。实验结果表明，采用该分布式拉曼测温系统能够在野外大范围监测河道流水剖面的温度分布。
Abstract:Aiming at the problem that the length of the temperature sensing area is smaller than the spatial resolution, which leads to inaccurate temperature measurement, a multi-point temperature correction scheme was proposed, and a device and method that could be applied to the measurement of temperature distribution in medium and small scale sections were introduced. Using the scheme based on multi-modal Gaussian fitting, the simultaneous correction of multiple inaccurate temperature points was realized, and the temperature accuracy of the system was improved. A prototype of the distributed Raman thermometry system was designed, and the measurement of temperature distribution in medium and small scale sections was realized by wrapping the fiber around the PPR tube. The experimental device improved the measurement accuracy of 2 m and 1.5 m temperature sensing areas to 1 °C, the thickness of the ice sheet of the Wanjiazhai Reservoir on the Yellow River in winter was measured to be 42.33 cm, which was 1.67 cm different from the artificial measurement results, and the measured temperature of 36 m deep river water profile fluctuated in the range of 0.32 ℃ ~ 0.90 ℃. The experimental results show that the distributed Raman thermometry system can monitor the temperature distribution of the river flow profile on a large scale in the field.
-
引言
单目光栅投影三维测量系统由相机、投影仪和计算机组成。投影仪向待测物体投影一组编码图案，图案受到物体表面形貌的调制产生形变，再由相机捕获；通过解码调制后的编码图案，获得待测物体表面上各点在相机和投影仪的图像坐标，最后根据光学三角法计算待测物表面各点在世界坐标系下的三维坐标[1-3]。光栅投影三维测量技术具有非接触、测量速度快、自动化等优点，近年来在逆向工程、质量检测、文物保护等领域得到广泛应用。
其中相机和投影仪的标定是系统工作的基本步骤与关键环节，其精度直接决定了系统的测量精度。
圆形特征对噪声的抑制性较强，识别率较高[4]，且在投影仪标定时，标志点中心处的相位调制度较好，故圆形阵列标定板常作为三维测量系统的标定物。圆形特征的提取在图像处理领域研究较早，对图像中圆形及椭圆的识别具有较高的准确度与速度[5-7]，标定场景中标定板图像特征简单，标志点与非标志点在轮廓的闭合性、轮廓长度、圆度、对比度等特征上区分明显，综合使用以上方式可以准确快速提取标志点图像[8，10，12]。
圆形标志点像素坐标的获取精度主要由边缘检测精度决定，朱统晶等基于正交傅里叶—马林距的算法来获取亚像素边缘[9]，张虎等人通过边缘梯度方向的灰度插值拟合高斯曲线进行亚像素边缘定位[10]。除传统的插值、拟合、基于距等方法外，Agustin等人提出了一种基于边缘局部区域效应的亚像素边缘提取方法[11]。此外，标志点图像坐标需要与空间坐标匹配，夏仁波等人通过标志点之间的几何位置关系，对标志点进行了排序和匹配[12]，以实现相机的自动标定。
投影仪使用逆向相机模型，使用相移轮廓术进行投影仪像素编码，通过相位建立标志点的投影仪图像坐标与相机图像坐标的对应关系[13]，进而与标志点的空间坐标相匹配来进行标定。光栅条纹会受到物体表面调制及环境光干扰等影响，引起相位误差。通过补偿、校正等方式可以减小相位误差，毛翠丽等人提出多频法和反向相位误差补偿结合的高效方法[14]，郑东亮等人提出双四步相移的方法[15]，有效减小非线性引起的相位误差。标志点圆心坐标往往是亚像素级别的，圆心处的相位值需要通过插值获得，Jakob Wilm等人使用高斯径向基函数方法进行插值[16]，进一步提高圆心处相位获取精度与稳定性。
本文考虑边缘变化剧烈程度和噪声影响，根据边缘特征灵活提出了基于拟合及局域区域效应的亚像素边缘检测算法，获得高精度标志点圆心图像坐标；同时提供一种基于透视变换图像矫正的简易快速标志点匹配方法。投影仪标定中，除使用双四步相移减小误差之外，针对插值算法，提出一种基于径向基的线性插值方式。
1 系统模型
单目光栅投影三维测量系统模型如图1所示。其中，世界坐标系
${{{O}}_{{w}}} - {{{X}}_{{w}}}{{{Y}}_{{w}}}{{{Z}}_{{w}}}$ 是物体在空间中所建立的绝对坐标系，以标定时的标定物为基础而建立，重建的点云通过该坐标系来描述；相机坐标系${O_{{c}}} - {{{X}}_{{c}}}{{{Y}}_{{c}}}{{{Z}}_{{c}}}$ 是以相机镜头光心为原点建立的坐标系；相机图像坐标系${{{o}}_{{c}}} - {{{u}}_{{c}}}{{{v}}_{{c}}}$ 是以图像的左上角为原点建立的二维坐标系；投影仪坐标系与投影仪图像坐标系同理。在相机的小孔成像模型中，相机的成像过程可以简化为透视投影模型：
$${{s}}\left[ {\begin{array}{*{20}{c}} {{{{u}}_{{c}}}} \\ {{{{v}}_{{c}}}} \\ 1 \end{array}} \right] = {{A}}\left[ {\begin{array}{*{20}{c}} {{{{R}}_{3 \times 3}}}&{{{{T}}_{3 \times 1}}} \\ {{0_{1 \times 3}}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{{X}}_{{w}}}} \\ {{{{Y}}_{{w}}}} \\ {{{{Z}}_{{w}}}} \\ 1 \end{array}} \right] = {{{M}}_{{c}}}\left[ {\begin{array}{*{20}{c}} {{{{X}}_{{w}}}} \\ {{{{Y}}_{{w}}}} \\ {{{{Z}}_{{w}}}} \\ 1 \end{array}} \right]$$ (1) 式中：
${{s}}$ 为尺度因子；${{A}}$ 为相机的内参矩阵；${{{R}}_{3 \times 3}}$ 为世界坐标系到相机坐标系的正交旋转矩阵；${{{T}}_{3 \times 1}}$ 为平移向量；${{{M}}_c}$ 是3×4大小的相机投影变换矩阵。使用逆向相机模型对投影仪进行数学建模，世界坐标系下物点${P_{{W}}}\left( {{X_{{W}}},{Y_{{W}}},{Z_{{W}}}} \right)$ 到投影仪图像${P_{{p}}}\left( {{u_{{p}}},{v_{{p}}}} \right)$ 的转换关系同相机一致，描述为$${{s}}\left[ {\begin{array}{*{20}{c}} {{{{u}}_{{p}}}} \\ {{{{v}}_{{p}}}} \\ 1 \end{array}} \right] = {{{M}}_{{p}}}\left[ {\begin{array}{*{20}{c}} {{{{X}}_{{w}}}} \\ {{{{Y}}_{{w}}}} \\ {{{{Z}}_{{w}}}} \\ 1 \end{array}} \right]$$ (2) ${{{M}}_{{c}}}$ 、${{{M}}_{{p}}}$ 进一步表示为$${{{M}}_{{c}}} = \left( {\begin{array}{*{20}{c}} {{{m}}_{11}^{{c}}}&{{{m}}_{12}^{{c}}}&{{{m}}_{13}^{{c}}}&{{{m}}_{14}^{{c}}} \\ {{{m}}_{21}^{{c}}}&{{{m}}_{22}^{{c}}}&{{{m}}_{23}^{{c}}}&{{{m}}_{24}^{{c}}} \\ {{{m}}_{31}^{{c}}}&{{{m}}_{32}^{{c}}}&{{{m}}_{33}^{{c}}}&{{{m}}_{34}^{{c}}} \end{array}} \right)$$ (3) $${{{M}}_{{p}}} = \left( {\begin{array}{*{20}{c}} {{{m}}_{11}^{{p}}}&{{{m}}_{12}^{{p}}}&{{{m}}_{13}^{{p}}}&{{{m}}_{14}^{{p}}} \\ {{{m}}_{21}^{{p}}}&{{{m}}_{22}^{{p}}}&{{{m}}_{23}^{{p}}}&{{{m}}_{24}^{{p}}} \\ {{{m}}_{31}^{{p}}}&{{{m}}_{32}^{{p}}}&{{{m}}_{33}^{{p}}}&{{{m}}_{34}^{{p}}} \end{array}} \right)$$ (4) 通过系统标定，可得投影变换矩阵
${{{M}}_{{c}}}$ 和${{{M}}_{{p}}}$ ，由（1）式～（4）式可求得物体表面的世界坐标[3]为$$\;\!\! \begin{split} \left(\!\! {\begin{array}{*{20}{c}} {{{{X}}_{{w}}}} \\ {{{{Y}}_{{w}}}} \\ {{{{Z}}_{{w}}}} \end{array}}\!\! \right) = & {\left(\!\! {\begin{array}{*{20}{c}} {{{m}}_{11}^{{c}} - {{m}}_{31}^{{c}}{{{u}}_{{c}}}}&{{{m}}_{12}^{{c}} - {{m}}_{32}^{{c}}{{{u}}_{{c}}}}&{{{m}}_{13}^{{c}} - {{m}}_{33}^{{c}}{{{u}}_{{c}}}} \\ {{{m}}_{21}^{{c}} - {{m}}_{31}^{{c}}{{{v}}_{{c}}}}&{{{m}}_{22}^{{c}} - {{m}}_{32}^{{c}}{{{v}}_{{c}}}}&{{{m}}_{23}^{{c}} - {{m}}_{33}^{{c}}{{{v}}_{{c}}}} \\ {{{m}}_{21}^{{p}} - {{m}}_{31}^{{p}}{{{v}}_{{p}}}}&{{{m}}_{22}^{{p}} - {{m}}_{32}^{{p}}{{{v}}_{{p}}}}&{{{m}}_{23}^{{p}} - {{m}}_{33}^{{p}}{{{v}}_{{p}}}} \end{array}} \!\!\right)^{ - 1}} \times\\ &\left( {\begin{array}{*{20}{c}} {{{m}}_{14}^{{c}} - {{m}}_{34}^{{c}}{{{u}}_{{c}}}} \\ {{{m}}_{24}^{{c}} - {{m}}_{34}^{{c}}{{{v}}_{{c}}}} \\ {{{m}}_{24}^{{p}} - {{m}}_{34}^{{p}}{{{v}}_{{p}}}} \end{array}} \right)\end{split}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!$$ (5) 由上述公式可以看出，相机和投影仪的标定精度直接决定了三维扫描系统的工作精度。
2 相机标定
使用圆形阵列标定板进行相机标定时，需要拍摄不同方位的标定板，并建立各个标志点的世界坐标与相机图像坐标的对应关系。标志点圆心的图像坐标直接参与标定计算，提取的精度决定相机标定的精度。
2.1 标志点边缘提取
标志点图像特征简单，由黑色背景和白色圆形填充组成。由方形孔径采样定理[17]，标志点由一系列离散的灰度阵列描述，设对应标志点内、外的像素灰度值分别为
${{{f}}_{{A}}}$ 、${{{f}}_{{B}}}$ ，当像素点内有边缘穿过时，像素落在标志点内外两个部分，该像素的灰度值可以表示为$${{{f}}_{{E}}} = {{{f}}_{{A}}} \times {{{S}}_{{A}}} + {{{f}}_{{B}}} \times {{{S}}_{{B}}}$$ (6) 式中：
${{{S}}_{{A}}}$ 、${{{S}}_{{B}}}$ 分别为该像素点在标志点内、外所占的面积，设像素宽、高均为1，则有${{{S}}_{{A}}} + {{{S}}_{{B}}} = 1$ 。圆形标志点通过投影变换成像时，往往退化成椭圆，而且由于相机的离焦等影响，实际空间中剧变的灰度经过光学成像成为渐变的形式[17-18]。这种变化可以看作是系统点扩展函数与图像卷积的结果。图像边缘灰度值变化表征为高斯分布，标志点边缘退化模型可以用典型的二维高斯点扩展函数描述为
$${\rm{PSF}}(x,y) = \frac{1}{{2{\rm{{\text{π}} }}{\sigma _{{p}}}^2}}\exp \left( { - \left( {{x^2} + {y^2}} \right)/2{\sigma _{{p}}}^2} \right)$$ (7) 式中：
${{{\sigma }}_{{p}}}$ 为点扩展函数方差，方差越大，代表像素点灰度值受周围像素影响越大。边缘的精确提取还受到噪声的影响，通过分析标志点内外的灰度分布均匀区域的图像，相机拍摄到图像包含随机噪声[19]，信噪比（signal noise ratio，SNR）为18 dB，SNR定义为
$${\rm{SNR}} = 20{\log _{10}}\frac{{{k}}}{{{{{\sigma }}_{{r}}}}}$$ (8) 式中:
${{k}}$ 为图像的对比度；${{{\sigma }}_{{r}}}$ 为随机噪声的标准差。根据上述分析，采用图像模拟方式生成标志点图像[4]。为直观体现边缘退化情况，本文提出了四分区的模拟标志点图像进行边缘检测验证，如图2（a）所示，四个分区1~4分别为理想边缘图像、添加随机噪声的边缘图像、使用点扩展函数与边缘图像卷积后添加随机噪声、使用点扩展函数与边缘图像卷积。
根据文献[10]所提供的方法，将边缘的二维高斯曲面拟合转化为一维的高斯曲线拟合，其亚像素边缘提取分为4个步骤：
1） 提取包含椭圆的感兴趣区域（region of interesting，ROI），使用Canny算子定位像素级的单边缘；
2） 边缘点梯度方向的确定以及基于梯度方向的区域划分；
3） 梯度方向的灰度插值以及灰度差求取；
4） 高斯曲线拟合。
步骤3中求出的灰度差分曲线可以使用高斯曲线描述，曲线极值点即梯度方向上的边缘，描述为
$${{{N}}_{{{k,\sigma ,\mu }}}}{{(x) = k}}\frac{{{1}}}{{\sqrt {{{2{\text{π}} }}{{{\sigma }}^{{2}}}} }}{{{\exp}}{\left({{ - }}\frac{{{{{{(x - \mu )}}}^{{2}}}}}{{{{2}}{{{\sigma }}^{{2}}}}}\right)}}$$ (9) 式中：
${{N}}({{x}})$ 为梯度方向上插值点x的灰度差分值；${{k,\sigma ,\mu }}$ 分别为高斯曲线系数、方差和均值，其中$k$ 对应曲线极值点大小，${\rm{\mu}} $ 对应极值点位置。对（9）式两边取对数，等式变成关于x的多项式：
$${{{{\ln}}(N(x)) = {{\ln}}}}\Bigg(k\frac{{{1}}}{{\sqrt {{{2{\text{π}} }}{{{\sigma }}^{{2}}}} }}{{{\exp}}{\left({{ - }}\frac{{{{{{(x - \mu )}}}^{{2}}}}}{{{{2}}{{{\sigma }}^{{2}}}}}\right)}}\Bigg){{ = }}{{{p}}_{{1}}}{{{x}}^{{2}}}{{ + }}{{{p}}_{{2}}}{{x + }}{{{p}}_{{3}}}$$ (10) 当存在2个以上有效灰度插值时，拟合多项式可以得到
${{\mu}} = {{ - {{{p}}_2}} / {2{{{p}}_1}}}$ 。使用文献[11]中的方法对模拟生成的标志点图像进行亚像素边缘提取，结果如图2（b）所示。可以看到，对于图中区域2中的边缘，该方法提取效果并不理想。如图3所示， 以区域2中的点
${{{P}}_{({{x}},{{y}})}}$ 为例，该点处的边缘变化剧烈且灰度差分分布曲线受噪声影响，使用最小二乘法拟合时求得高斯曲线定点偏差较大。![]()
图 3${{P}}$ 点局部像素分布及高斯曲线拟合示意图Figure 3. Point P partial pixel distribution and schematic of Gaussian curve fitting实际使用中，文献[10]中的方法在
${{{\sigma }}_{{p}}}$ 比较大的情况下表现较好，但在${{{\sigma }}_{{p}}}$ 较小的情况下表现较差。计算边缘梯度方向上插值点的灰度差分值，本文取阈值为10，当差分值中大于阈值的个数小于4的时候，判断边缘变化剧烈。变化剧烈的边缘处，使用基于边缘灰度值分布局部区域效应的方法[11]，替换高斯曲线拟合方法，进行亚像素边缘提取，如图4所示。![]()
图 4 二阶曲线模拟边缘示意图及局部区域效应边缘提取Figure 4. Schematic of second order curve simulated edge and partial area effect edge extraction图4（a）中的边缘用虚线表示，二次曲线可以表示在点
${{{P}}_{{(x,y)}}}$ 附近的边缘走向，以${{{P}}_{{(x,y)}}}$ 为原点建立坐标轴，列出曲线方程为$${{y = {{a}} + {{b}}x + {{c}}}}{{{x}}^{{2}}}$$ (11) 求解系数a、b、c后，可以计算曲线边缘点
${{{P}}_{{(x,y)}}}$ 处的亚像素位置。在已知边缘像素点及周围像素灰度值的情况下，以该点为中心选择周围像素建立方程与灰度值的对应关系，可以求解方程系数。首先进行梯度方向区域划分，划分标准是边缘点在灰度梯度方向上的直线斜率
${{{k}}_{\rm{l}}}$ 的大小。当$\left| {{{{k}}_{\rm{l}}}} \right| < 1$ 时，为k1区域，否则为k2区域。已知点${{{P}}_{{(x,y)}}}$ 属于k2区域，选择该点周围$5 \times 3$ 区域的像素，当属于k1区域时，取$3 \times 5$ 区域的像素进行计算。如图4（a）所示，区域分为3列，计算每一列的灰度值之和，以${{{S}}_{{L}}}{\text{、}}{{{S}}_{{M}}}{\text{、}}{{{S}}_{{R}}}$ 表示为$$\left\{ \begin{array}{l} {{{S}}_{{L}}}{{ = }}\displaystyle\sum\limits_{{{n = y - 2}}}^{{{y + 2}}} {{{f(x - 1,n) = 5}}{{{f}}_{{B}}}} {{ + (}}{{{f}}_{{A}}}{{ - }}{{{f}}_{{B}}}{{)L}}\\ {{{S}}_{{M}}}{{ = }}\displaystyle\sum\limits_{{{n = y - 2}}}^{{{y + 2}}} {{{f(x,n) = 5}}{{{f}}_{{B}}}} {{ + (}}{{{f}}_{{A}}}{{ - }}{{{f}}_{{B}}}{{)M}}\\ {{{S}}_{{R}}}{{ = }}\displaystyle\sum\limits_{{{n = y - 2}}}^{{{y + 2}}} {{{f(x + 1,n) = 5}}{{{f}}_{{B}}}} {{ + (}}{{{f}}_{{A}}}{{ - }}{{{f}}_{{B}}}{{)R}} \end{array} \right. \quad \quad \quad\;\; $$ (12) 式中：
$L{\text{、}}M{\text{、}}R$ 表示区域中每一列在标志点内部，即边缘内侧的面积，带入曲线方程使用积分求得$$\;\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\left\{ \begin{array}{l} {{L = }}\displaystyle\int_{{{ - }}{{{3}} / {{2}}}}^{{{ - }}{{{1}} / {{2}}}} \Bigg({{{a + bx + c}}{{{x}}^{{2}}}} {{ + }}\frac{{{5}}}{{{2}}}\Bigg){{{\rm{d}}x = a - b + }}\frac{{{{13}}}}{{{{12}}}}{{c + }}\frac{{{5}}}{{{2}}}\\ {{M = }}\displaystyle\int_{{{ - }}{{{1}} / {{2}}}}^{{{{1}} / {{2}}}} \Bigg({{{a + bx + c}}{{{x}}^{{2}}}} {{ + }}\frac{{{5}}}{{{2}}}\Bigg){{{\rm{d}}x = a + }}\frac{{{1}}}{{{{12}}}}{{c + }}\frac{{{5}}}{{{2}}}\\ {{R = }}\displaystyle\int_{{{{1}} / {{2}}}}^{{{{3}} / {{2}}}} \Bigg({{{a + bx + c}}{{{x}}^{{2}}}} {{ + }}\frac{{{5}}}{{{2}}}\Bigg){{{\rm{d}}x = a + b + }}\frac{{{{13}}}}{{{{12}}}}{{c + }}\frac{{{5}}}{{{2}}} \end{array} \right.\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!$$ (13) 根据3个等式，可以求得曲线方程系数为
$$\left\{ \begin{array}{l} {{c = }}\dfrac{{{{{S}}_{{L}}}{{ + }}{{{S}}_{{R}}}{{ - 2}}{{{S}}_{{M}}}}}{{{{2(}}{{{f}}_{{A}}}{{ - }}{{{f}}_{{B}}}{{)}}}}\\ {{b = }}\dfrac{{{{{S}}_{{R}}}{{ - }}{{{S}}_{{L}}}}}{{{{2(}}{{{f}}_{{A}}}{{ - }}{{{f}}_{{B}}}{{)}}}}\\ {{a = }}\dfrac{{{{2}}{{{S}}_{{M}}}{{ - 5(}}{{{f}}_{{A}}}{{ + }}{{{f}}_{{B}}}{{)}}}}{{{{2(}}{{{f}}_{{A}}}{{ - }}{{{f}}_{{B}}}{{)}}}}{{ - }}\dfrac{{{1}}}{{{{12}}}}{{c}} \end{array} \right.\quad \quad \quad \quad \quad \quad \quad\quad\; $$ (14) 实际使用中，
${{{f}}_{{A}}}$ 、${{{f}}_{{B}}}$ 使用靠近所求边缘内外、侧附近像素点均值求得。使用该方法求得亚像素边缘如图4（b）所示，边缘获取精度大幅度提高。对比试验分别用方法1) 基于Canny算子边缘像素提取，方法2) 基于文献[10]中高斯曲线拟合的亚像素边缘提取，方法3) 本文所使用方法，对实际拍摄所得的标志点图像进行处理及椭圆拟合（标志点边缘及椭圆焦点由方法3得出）。由于实际拍摄的标志点无法得到圆心的准确位置，使用残差的方法进行对比，残差由以下方式给出：首先用边缘点进行椭圆拟合，得到椭圆方程；计算所有参与拟合的点到椭圆方程曲线的距离，距离均值则为残差。如表1所示，方法2的平均残差为0.1477，方法3的平均残差为0.0871，本文方法对比方法2，残差精度提高41%。
表 1 拟合残差比较Table 1. Comparison of fitting residual error标志点 方法 残差 标志点 方法 残差 
1 0.5295 
1 0.5589 2 0.1621 2 0.1428 3 0.0789 3 0.0857 
1 0.5353 
1 0.5881 2 0.1302 2 0.1731 3 0.0694 3 0.1012 
1 0.5522 
1 0.5512 2 0.1301 2 0.1481 3 0.0914 3 0.0958 2.2 标志点匹配
对标志点进行排序，建立标志点的图像坐标与其世界坐标一一的对应关系。当相机正视标定板时，标志点横平竖直地排列，容易进行排序匹配，透视变换可以将图像投影到一个新的视平面，使用透视变换矩阵将倾斜的标定板进行矫正得到
$$\left[ {\begin{array}{*{20}{c}} {{{x'}}}\\ {{{y'}}}\\ {{1}} \end{array}} \right]{{ = A}}\left[ {\begin{array}{*{20}{c}} {{x}}\\ {{y}}\\ {{1}} \end{array}} \right]{{,A = }}\left[ {\begin{array}{*{20}{c}} {{a_{{{11}}}}}&{{a_{{{12}}}}}&{{a_{{{13}}}}}\\ {{a_{{{21}}}}}&{{a_{{{22}}}}}&{{a_{{{23}}}}}\\ {{a_{{{31}}}}}&{{a_{{{32}}}}}&{{a_{{{33}}}}} \end{array}} \right]$$ (15) 式中：
$\left( {{{x}},{{y}},1} \right)$ 为图像坐标；$({{x}}',{{y}}',1)$ 为矫正后的图像坐标；A为透视矩阵，通常令${{{a}}_{33}}$ 值为1。透视矩阵为8自由度，因此透视矩阵只需要4个对应的点就可以解出。透视变换矫正步骤如下：1） 提取到图像中标志点后，通过面积筛选出定位标志点，如图5（a）所示；
![]()
图 5 透视变换标志点匹配过程及相机重投影误差Figure 5. Process of mark points matching by perspective transformation and re-projection error of camera2） 通过定位标志点之间的距离，首先为每个定位标志点进行区分，预算出每个定位标志点透视变换矫正后的图像坐标；
3） 求出透视矩阵；
4） 使用（16）式求出标定板图像透视变换矫正后的图像坐标为
$$ \left\{ \begin{array}{l} {{x}}' = \dfrac{{{a_{11}} \times {{x}} + {a_{21}} \times {{y}} + {a_{31}}}}{{{a_{13}} \times {{x}} + {a_{23}} \times {{y}} + {a_{33}}}} \\ {{y}}' = \dfrac{{{a_{12}} \times {{x}} + {a_{22}} \times {{y}} + {a_{32}}}}{{{a_{13}} \times {{x}} + {a_{23}} \times{{ y}} + {a_{33}}}} \end{array} \right. $$ (16) 5） 此时矫正后的标志点横平竖直排列，如图5（b）所示，可以进行排序和匹配。
完成上述操作后，可以使用张氏标定法对相机进行快速标定，相机标定的重投影误差均值为0.0524像素，误差分布如图5（d）所示。
3 投影仪标定
投影仪使用相移轮廓术进行像素编码[20]，相移条纹图片的生成公式为
$$ \begin{split} &{{{I}}_{{n}}}({{x}},{{y}}) = {{{I}}_{{A}}}({{x}},{{y}}) + {{{I}}_{{B}}}({{x}},{{y}})\cos \left( {{{\varphi }}({{x}},{{y}}) + \frac{{2{{{\text{π}}}}({{n}} - 1)}}{{{N}}}} \right),\\ &{{n}} = 1,2,\cdots,{{N}} \end{split} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!$$ (17) 式中：
${{{I}}_{{n}}}$ 为图片在点$({{x}},{{y}})$ 处的灰度值；${{{I}}_{{A}}}$ 为条纹图案的背景光强；${{{I}}_{{B}}}$ 为条纹调制强度；${{\varphi}} ({{x}},{{y}})$ 为待求相位，${{\varphi}} ({{x}},{{y}}) = {{2{{\pi}} {{x}}} / {{p}}}$ ，${{p}}$ 为条纹周期；N为相移步数；n表示第几次相移。实际使用中，生成的图像从投影出来到相机捕获的过程还会引入非线性的误差，这种非线性往往不能用Gamma变换等来描述。系统的灰度响应如图6（a）所示，这种非线性导致正弦条纹变形，从而引起周期性的相位误差，如图6（b）。本文使用文献[15]中双四步相移方式进行误差抑制。
在投影仪标定时，分别投影上述两组水平、垂直条纹图像序列并使用相机捕获，另拍摄一张不包含条纹的图片。使用多频条纹，基于时域方式进行相位展开。条纹图像序列分别解码得到水平、垂直的绝对相位，由标志点圆心的绝对相位值可以计算出标志点的投影仪图像坐标[13]，如图7所示。
拍摄的不包含条纹的标定板图片可以提取亚像素精度的标志点圆心图像坐标。将标定板的绝对相位看做关于图像坐标的二维函数，通过插值方式获取圆心的绝对相位值，常用插值方式为双线性插值，公式为
$${{\theta }}({{x}},{{y}}) \approx \left[ {\begin{array}{*{20}{c}} {1 - {{{\lambda}} _{{x}}}}&{{{{\lambda}} _{{x}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\theta}} ({{{x}}_0},{{{y}}_0})}&{{{\theta }}({{{x}}_0},{{{y}}_1})} \\ {{{\theta}} ({{{x}}_1},{{{y}}_0})}&{{{\theta}} ({{{x}}_1},{{{y}}_1})} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {1 - {{{\lambda}} _{{y}}}} \\ {{{{\lambda}} _{{y}}}} \end{array}} \right]$$ (18) 式中：
$\theta (x,y)$ 表示该点处插值后的相位值，${{\theta}} ({{{x}}_0},{{{y}}_0})$ 表示图7中对应已知点的相位值，${{{\lambda}} _{{x}}} = {{x}} - \left[ {{x}} \right]$ ，${{{\lambda}} _{{y}}} = {{y}} - \left[ {{y}} \right]$ ，$\left[ {{x}} \right]{\text{、}}\left[ {{y }}\right]$ 表示取整数部分。使用双线性差值的方式，进行投影仪标定，重投影误差均值为0.1582，误差分布如图8（a）所示，分布特征明显。文献[16]中提到使用基于径向基函数的插值方法，将插值问题转化成曲面重构问题，使用该方法求出标志点圆心的投影仪图像坐标，标定后的重投影误差均值为0.1602，误差分布如图8（b）所示。考虑到存在噪声等因素影响，本文提出一个基于径向基的插值方式为$$ {{r}}({{{x}}_{{i}}},{{{y}}_{{j}}}) = \left\| {({{x}},{{y}}) - ({{{x}}_{{i}}},{{{y}}_{{j}}})} \right\|,\;{{i}},{{j}} \in ( - 3, - 2,\cdots ,3)\;\;\;\;$$ (19) 式中：
${{r}}({{{x}}_{{i}}},{{{y}}_{{j}}})$ 是点$({{{x}}_{{i}}},{{{y}}_{{j}}})$ 到待插值点的欧氏距离，遍历待插值点周围已知点，本文取所有欧氏距离小于2.6的已知点，令${{{\sigma}} ^2}{\rm{ = }}\sqrt 2$ ，经式（20）变换后参与计算得$$R({{{x}}_{{i}}},{{{y}}_{{j}}}) = {{{e}}^{ - \frac{{{{{r}}^2}({{{x}}_{{i}}},{{{y}}_{{j}}})}}{{{{{\sigma }}^2}}}}}$$ (20) $${{\theta}} ({{x}},{{y}}) = \sum {\Bigg(\frac{{{{R}}({{{x}}_{{i}}},{{{y}}_{{j}}})}}{{{{{S}}_{{R}}}}} \times {{\theta }}} ({{{x}}_{{i}}},{{{y}}_{{j}}})\Bigg)$$ (21) $${{{S}}_{{R}}} = \sum {{{R}}({{{x}}_{{i}}},{{{y}}_{{j}}})} $$ (22) 该方法是线性插值的方式，将已知点到待插值点径向距离的高斯变换作为权值，使用
${{{S}}_{{R}}}$ 保证权值和为1。使用本文方法对投影仪进行标定，重投影误差均值为0.1203，误差分布如图8（c）所示，对比使用双线性插值和基于径向基函数的插值方法，标定精度分别提高23.9%、24.9%。4 实验及结果
搭建如图9所示单目光栅投影三维测量系统对本文算法进行验证。使用PointGrey公司FL3-U3-13Y3M-C COMS相机，分辨率1280×1024，相机搭配Computar 8 mm定焦镜头；投影仪为TI公司DLP LightCrafter 4500数字微镜（DMD）投影仪，本征分辨率为912×1140，偏轴率为100%；使用9行11列圆形阵列标定板，标志点直径18 mm，中心距离45 mm。使用本文标定方法进行系统标定，结果如表2所示。
表 2 系统标定结果Table 2. System calibration results相机参数 投影仪参致 内参矩阵 1730.1713 0 628.4165 1122.7932 0 421.2725 0 1730.1724 518.3202 0 1122.7905 584.0842 0 0 1 0 0 1 径向畸变参数 k1=−0.0907
k2=0.2018k1=0.0542
k2=−0.1328旋转参数 2.0565 1.9864 −0.177 2.2336 2.1582 −0.2147 平移参致/mm −102.5556 −39.5949 378.6716 −78.2600 −138.9417 382.8964 单目光栅投影三维测量系统完成标定后，对标准块进行测量，测量距离为1050 mm。如图10所示，导出点云到Geomagic Studio软件中进行分析，本文方法测得标准块正面有效点云30 290个。将左中右三个平面进行拟合，标准差分别为0.221 mm、0.232 mm、0.210 mm，系统相对精度约为1:5000。
使用单目光栅投影三维测量系统进行人脸扫描，重建效果较好，如图11所示。
![]()
图 11 三维测量系统人脸重建结果Figure 11. Reconstruction results of face with 3D measurement system5 结论
本文研究了单目光栅投影三维测量系统标定的两个关键步骤。提出了基于拟合和局部区域效应相结合的亚像素边缘检测方法，对标志点边缘进行高精度提取，标志点边缘平均残差为0.0871。对比基于拟合方法提取的边缘残差，精度提高了41%，并提供了一种快速简单的、基于透视矫正的标志点自动化匹配方法，相机标定重投影误差均值为0.0524像素。投影仪标定中，分析相位误差产生原因并使用双四步相移减小误差。针对插值算法，提出一种基于径向基的线性插值方式，投影仪标定重投影误差均值为0.1203像素，对比双线性插值方式，精度提高23.9%。使用上述标定方法，有效地提高了单目光栅三维测量系统的标定精度。
-
![]()
图 2 标准温度与光纤测量温度之间的线性关系
Figure 2. Linear relationship between standard temperature and fiber-measured temperature
![]()
图 3 高斯拟合处理前后的温度分布
Figure 3. Temperature distribution before and after Gaussian fitting
![]()
图 8 冰盖垂直方向上的温度分布
Figure 8. Temperature distribution in vertical direction of ice sheet
-
[1] 付辉, 杨开林, 王涛, 等. 河冰水力学研究进展[J]. 南水北调与水利科技,2010,8(1):14-18. FU Hui, YANG Kailin, WANG Tao, et al. Progress in the study of river ice hydraulics[J]. South-to-North Water Transfers and Water Science & Technology,2010,8(1):14-18.
[2] 滕晖, 邓云, 黄奉斌, 等. 水库静水结冰过程及冰盖热力变化的模拟试验研究[J]. 水科学进展,2011,22(5):720-726. TENG Hui, DENG Yun, HUANG Fengbin, et al. Experimental study on the simulation of freezing processes in calm waters and thermal changes on reservoir ice cover[J]. Advances in Water Science,2011,22(5):720-726.
[3] 闫慧荣, 刘文涛, 冯民权, 等. 水库冰冻期冰内温度的模拟[J]. 武汉理工大学学报(信息与管理工程版),2010,32(5):738-741. YAN Huirong, LIU Wentao, FENG Minquan, et al. Distribution of ice temperature in freeze reservoir[J]. Journal of Wuhan University of Technology (Information & Management Engineering),2010,32(5):738-741.
[4] CUI L Q, QIN J M, DENG X. Freshwater ice thickness apparatus based on differences in electrical resistance and temperature[J]. Cold Regions Science and Technology,2015,119:37-46. doi: 10.1016/j.coldregions.2015.07.009
[5] SOTO M A, NANNIPIERI T, SIGNORINI A, et al. Raman-based distributed temperature sensor with 1 m spatial resolution over 26 km SMF using low-repetition-rate cyclic pulse coding[J]. Optics Letters,2011,36(13):2557-2559. doi: 10.1364/OL.36.002557
[6] SELKER J S, THÉVENAZ L, HUWALD H, et al. Distributed fiber-optic temperature sensing for hydrologic systems[J]. Water Resources Research,2006,42(12):235-241.
[7] ZHANG L, FENG X, ZHANG W, et al. Improving spatial resolution in fiber Raman distributed temperature sensor by using deconvolution algorithm[J]. Chinese Optics Letters,2009,7(7):560-563. doi: 10.3788/COL20090707.0560
[8] BAZZO J P, PIPA D R, MARTELLI C, et al. Improving spatial resolution of Raman DTS using total variation deconvolution[J]. IEEE Sensors Journal,2016,16(11):4425-4430. doi: 10.1109/JSEN.2016.2539279
[9] 石希, 夏军强, 孙健. 基于热红外遥感影像的河流水温反演方法比较——以长江上游流域为例[J]. 湖泊科学,2022,34(1):307-319. doi: 10.18307/2022.0125 SHI Xi, XIA Junqiang, SUN Jian. Comparison of methods to derive river water temperature using thermal infrared imagery: a case study of the Upper Yangtze River catchment[J]. Journal of Lake Sciences,2022,34(1):307-319. doi: 10.18307/2022.0125
[10] TYLER S W, HOLLAND D M, ZAGORODNOV V, et al. Using distributed temperature sensors to monitor an Antarctic ice shelf and sub-ice-shelf cavity[J]. Journal of Glaciology,2013,59(215):583-591. doi: 10.3189/2013JoG12J207
[11] KOBS S, HOLLAND D M, ZAGORODNOV V, et al. Novel monitoring of Antarctic ice shelf basal melting using a fiber-optic distributed temperature sensing mooring[J]. Geophysical Research Letters,2014,41(19):6779-6786. doi: 10.1002/2014GL061155
[12] 王玎睿, 邓霄, 张均, 等. 面向冰盖剖面的高空间分辨率分布式光纤测温系统设计[J]. 应用光学,2021,42(5):941-948. doi: 10.5768/JAO202142.0508003 WANG Dingrui, DENG Xiao, ZHANG Jun, et al. High spatial resolution distributed optical fiber temperature measurement system for ice cover profile[J]. Journal of Applied Optics,2021,42(5):941-948. doi: 10.5768/JAO202142.0508003
[13] 王伟杰. 基于拉曼散射的分布式光纤测温系统设计及优化[D]. 济南: 山东大学, 2013. WANG Weijie. Design and optimization of distributed fiber temperature sensor based on Raman scattering[D]. Jinan: Shandong University, 2013.
[14] SAXENA M K, RAJU S D V S J, ARYA R, et al. Empirical mode decomposition-based detection of bend-induced error and its correction in a Raman optical fiber distributed temperature sensor[J]. IEEE Sensors Journal,2016,16(5):1243-1252. doi: 10.1109/JSEN.2015.2499242
[15] 方一竹, 余果, 李海涛, 等. 基于多峰高斯模型的天然气产量趋势预测[J]. 天然气勘探与开发,2018,41(3):65-69. FANG Yizhu, YU Guo, LI Haitao, et al. Prediction of natural-gas production trend based on multi-peak Gauss model[J]. Natural Gas Exploration and Development,2018,41(3):65-69.
[16] 宁枫, 朱永, 崔海军, 等. 一种提高分布式光纤测温系统空间分辨率的线性修正算法[J]. 光子学报,2012,41(4):408-413. doi: 10.3788/gzxb20124104.0408 NING Feng, ZHU Yong, CUI Haijun, et al. A linear correcting algorithm for improving space resolution of distributed optical fiber Raman temperature measurement system[J]. Acta Photonica Sinica,2012,41(4):408-413. doi: 10.3788/gzxb20124104.0408
[17] BRIGGS M A, LAUTZ L K, MCKENZIE J M, et al. Using high-resolution distributed temperature sensing to quantify spatial and temporal variability in vertical hyporheic flux[J]. Water Resources Research,2012,48(2):670-677.
-
期刊类型引用(7)
1. 刘鸿智，童景琳，梁科. 基于光栅投影图像的硬质合金刀具磨损监测. 计算机仿真. 2024(04): 279-283 . 
百度学术
2. 甘勇，王昭太，朱明凯，文龙. 基于局部RANSAC的结构光标定方法. 应用光学. 2024(06): 1252-1260 . 本站查看
3. 严飞，祁健，刘银萍，吴迪，于强，刘佳. 基于格雷码的分区间相位展开方法. 应用光学. 2023(01): 79-85 . 本站查看
4. 郑小钰，董祉序，杨赫然，孙兴伟，刘寅. 基于结构光测量技术的DMD自适应掩膜生成. 电子测量与仪器学报. 2023(07): 148-155 . 百度学术
5. 邓春华，陈鑫源，汪普庆. 机器人视觉导航传感器光栅投射误差校正系统. 激光杂志. 2022(05): 182-186 . 百度学术
6. 赵婷婷，白福忠，徐永祥，高晓娟. 一种基于参数标定的投影条纹周期简易校正方法. 应用光学. 2021(01): 119-124 . 本站查看
7. 余宗璞，杨懿昕，王俨铮，刘雪岩，周平，周光泉. 三维人脸测量与分割系统研究. 应用光学. 2021(04): 664-670 . 本站查看
其他类型引用(9)
下载:
计量
- 文章访问数: 347
- HTML全文浏览量: 168
- PDF下载量: 37
- 被引次数: 16
陕公网安备 61011302001501号 