Nanosecond pulse laser damage characteristics of sub-wavelength gratings polarizers
-
摘要: 偏振片在诸多光学系统中有着重要的应用。亚波长介质光栅可用作正入射偏振片，在高能激光系统中有着广泛的应用前景。为了探究波长为1 064 nm的纳秒脉冲激光对于亚波长全介质光栅的诱导损伤特性，使用了粒子群优化算法结合严格耦合波分析设计了光栅的几何参数，计算表明亚波长光栅偏振片在入射光波长1 064 nm附近带宽0.5 nm内，平均消光比为1 500。使用了紫外曝光配合离子束刻蚀的工艺制备了HfO2光栅，并对其纳秒脉冲激光损伤阈值进行了测试。测试结果表明S光损伤阈值约为P光损伤阈值的5倍，且都大于5 J/cm2。结果表明亚波长全介质光栅偏振片可广泛用于正入射激光系统中。Abstract: Polarizers are important components for many optical systems. Sub-wavelength dielectric gratings can be used as polarizers under normal incidence, and these have broad application prospects in high-energy laser systems. In order to investigate the induced damage characteristics of nanosecond pulse laser with wavelength of 1 064 nm for subwavelength all dielectric gratings, the particle swarm optimization algorithm was used to design the geometric parameters of the grating with rigorous coupled wave analysis. The calculation results show that the subwavelength grating polarizer has 0.5 nm bandwidth and a theoretical extinction ratio of 1 500 at wavelength of 1 064 nm. This grating was fabricated by ultraviolet(UV) exposure and ion beam etching, and its nanosecond pulsed laser damage threshold was tested. The test results show that the laser induced damage threshold of S polarized light is about 5 times of that of P polarized light, and all of them are greater than 5 J/cm2. The results indicate that the sub-wavelength all dielectric grating polarizers can be widely used in laser system under normal incidence.
-
引言
单目光栅投影三维测量系统由相机、投影仪和计算机组成。投影仪向待测物体投影一组编码图案，图案受到物体表面形貌的调制产生形变，再由相机捕获；通过解码调制后的编码图案，获得待测物体表面上各点在相机和投影仪的图像坐标，最后根据光学三角法计算待测物表面各点在世界坐标系下的三维坐标[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%。使用上述标定方法，有效地提高了单目光栅三维测量系统的标定精度。
-
![]()
图 1 不同偏振光正入射时光栅0级反射效率
Figure 1. 0 order efficiency of different polarized lights at normal incidence
![]()
图 3 理论计算光栅0级次透射光谱
Figure 3. Theoretical calculation of 0 order transmission spectra of grating
![]()
图 6 实测光谱与理论计算光谱对比
Figure 6. Comparison between measured spectra and theoretical calculation spectra
-
[1] 李淑军, 姜会林, 朱京平, 等.偏振成像探测技术发展现状及关键技术[J].中国光学, 2013, 6(6):803-809. http://d.old.wanfangdata.com.cn/Periodical/zggxyyygxwz201306001 LI Shujun, JIANG Huilin, ZHU Jingping, et al.Development status and key technologies of polarization imaging detection[J].Chinese Optics, 2013, 6(6):803-809. http://d.old.wanfangdata.com.cn/Periodical/zggxyyygxwz201306001
[2] KIMEL A V, KIRILYUK A, RASING T.Femtosecond opto-magnetism:ultrafast laser manipulation of magnetic materials[J].Laser & Photonics Reviews, 2007, 1(3):275-287. http://cn.bing.com/academic/profile?id=b570b0e8e7514b78445b0f529dcd9651&encoded=0&v=paper_preview&mkt=zh-cn
[3] 马科斯·玻恩, 埃米尔·沃尔夫.光学原理[M].7版, 杨葭荪, 译.北京: 电子工业出版, 2006: 699-701. MAX B, EMIL W.Principles of optics[M].7th ed.Beijing: Electronie Industry Press, 2006: 699-701.
[4] EKINCI Y, SOLAK H H, DAVID C, et al.Bilayer Al wire-grids as broadband and high-performance polarizers[J].Optics Express, 2006, 14(6):2323-2334. doi: 10.1364/OE.14.002323
[5] TAMADA H, DOUMUKI T, YAMAGUCHI T, et al.Al wire-grid polarizer using the s-polarization resonance effect at the 0.8-microm-wavelength band[J].Optics Letters, 1997, 22(6):419-421. doi: 10.1364/OL.22.000419
[6] 唐晋发, 顾培夫, 刘旭, 等.现代光学薄膜技术[M], 杭州, 浙江大学出版社, 2007:167-171. TANG Jinfa, GU Peifu, LIU Xu, et al.Modern optical thin film technology[M].Hangzhou:Zhejiang University Press, 2007:167-171.
[7] LEE K J, CURZAN J, SHOKOOH-SAREMI M, et al.Resonant wideband polarizer with single silicon layer[J].Applied Physics Letters, 2011, 98(21):211112. doi: 10.1063/1.3594244
[8] LEE K J, GIESE J, AJAYI L, et al.Resonant grating polarizers made with silicon nitride, titanium dioxide, and silicon:design, fabrication, and characterization[J].Optics Express, 2014, 22(8):9271-9281. doi: 10.1364/OE.22.009271
[9] WANG S S, Magnusson R.Theory and applications of guided-mode resonance filters[J].Applied Optics, 1993, 32(14):2606-2613. doi: 10.1364/AO.32.002606
[10] MOHARAM M G, GAYLORD T K.Rigorous coupled-wave analysis of planar-grating diffraction[J].Journal of the Optical Society of America, 1981, 71(7):811-818. doi: 10.1364/JOSA.71.000811
[11] MOHARAM M G, POMMET D A, GRANN E B, et al.Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings:enhanced transmittance matrix approach[J].Journal of the Optical Society of America A, 1995, 12(5):1077-1086. doi: 10.1364/JOSAA.12.001077
[12] 王旭迪, 徐向东, 刘颖, 等.HfO2薄膜的离子束刻蚀特性研究[J].光学精密工程, 2004, 12(5):454-458. doi: 10.3321/j.issn:1004-924X.2004.05.002 WANG Xudi, XU Xiangdong, LIU Ying, et al.Ion beam etching of HfO2 film[J].Optics & Precision Engineering, 2004, 12(5):454-458. doi: 10.3321/j.issn:1004-924X.2004.05.002
[13] 徐向东, 刘颖, 邱克强, 等.HfO2顶层多层介质膜脉宽压缩光栅的离子束刻蚀[J].物理学报, 2013, 62(23):165-172. http://d.old.wanfangdata.com.cn/Periodical/wlxb201323025 XU Xiangdong, LIU Ying, QIU Keqiang, et al.Ion beam etching of multilayer dielectric pulse compressor gratings with top layer of HfO2[J].Acta Physica Sinica, 2013, 62(23):165-172. http://d.old.wanfangdata.com.cn/Periodical/wlxb201323025
-
期刊类型引用(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)
下载:
计量
- 文章访问数: 866
- HTML全文浏览量: 274
- PDF下载量: 52
- 被引次数: 16
陕公网安备 61011302001501号 