引言

三维重建在航空观测、计算机视觉等领域中占有重要地位，如地表地貌的立体测绘、多目立体视觉中目标识别与三维表面恢复，目标运动参数获取等技术^[1,2]，大多依赖于物体表面的纹理特征进行识别，但是实际任务中会经常遇到低纹理甚至无纹理、表面光滑易产生高反光的物体，传统的三维重建算法针对此类物体会出现大面积的数据空洞，无法获得完整的表面，即便是基于结构光的重建算法在遇到大面积耀光时也无法得到精确的深度信息，故开展针对低纹理高反光物体表面的三维重构算法具有重要意义.

近年来的研究表明，通过观测物体表面反射光中的偏振特性，可以获得物体的形状信息，且偏振所具有的“弱光强化，强光弱化”特性，使得即便外界光源不理想，仍然能得到较为清晰的偏振信息，从而进行三维表面恢复.三维重建的方法大致可以分为两种，一种是利用三维激光扫描仪等直接获取测区及环境的三维点云，这种方法精度较高，但是作业时间长，设备昂贵，受限于物体尺寸和材质；另一种是基于视觉的重建方法，多目立体视觉通过提取多张影像的特征点并匹配同名像点^[3]，然后利用影像间的几何关系重建三维表面^[4]，明暗恢复法（Shape From Shading, SFS）^[5]利用单张影像的明暗变化获取物体形状信息，认为灰度值与表面不同位置反射光强弱有关，光度立体法（Photometric Stereo, PS）^[6]通过改变入射光源得到多幅明暗不同的影像，利用目标表面同一个像素点在不同影像中的亮度信息得到几何信息，但是这些方法都要求已知光源的方向与强度等，不能很好的适用于大多数情况.基于偏振信息的三维重构方法不依赖物体的纹理特征，同时又能很好地抑制耀光，对低纹理、高反光的物体也能获得完整的三维表面形态.Koshikawa^[7]最早将偏振视觉引入三维重建领域；S.Rahmann和N. Canterakis等^[8]针对绝缘体反射光的偏振信息首次提出了完整的偏振三维重建算法；Ghosh等^[9]利用圆偏光分离了漫反射与镜面反射，完成了物体材质和形状的同时估计；Miche^[10]设计一种特定的视觉系统，可实现偏振图像采集与处理，通过积分过程实现三维重建；Kadambi等^[11]融合了Kinect得到的粗深度图和偏振图像得到的法向量，得到了较高精度的三维重建结果.该文章出发点与本文较为相似，但本文使用较少的偏振图像，对方位角处理和积分算法更简单高效，并且将结果与扫描仪真值进行了对比分析.国内近年来也开展了偏振三维重建技术的研究，平茜茜等^[12]结合双目立体视觉，提出偏振双目视觉三维重建算法，并获得目标的世界坐标系下绝对数据；岳春敏等^[13]提出利用解析反射光偏振图像重建透明物体三维表面的方法，并推导偏振度与表面法向量得函数关系，构建三维重建算法；彭群聂等^[14]针对视觉任务中表面光滑低纹理的目标，提出多尺度Shapelets算子获取有效深度信息，恢复三维形状.

针对偏振重建时存在的方位角歧义，本文先利用边缘传播算法进行粗纠正，再通过比较粗糙深度图的方位角偏差，重新对偏振信息获得的方位角进行估值；针对目前尚未有完善的梯度积分算法，本文采用融合深度图的表面三维重建算法，重建出较高精度的物体表面.

1 偏振三维重建

1.1　偏振三维重建原理

1809年，法国科学家马吕斯在双折射晶体的实验中发现了光的偏振现象，进一步证实光是横波^[15].偏振光可利用Stokes矢量来描述，使用四个参量表示光的强度和偏振信息，得到Stokes矢量需要三个及以上偏振角度的光强信息（本文采用四个角度来测量），在成像系统前放置线偏振片，通过调节偏振片的角度θpolar，光强I(θpolar,φ)与角度的关系为正弦函数如式（1）所示.

其中，光矢量沿垂直方向上的振动占优势，光强为Imax，与其垂直方向处于劣势，光强为Imin，φ为偏振方位角，偏振度（DoP, Degree of Polarization） 为：

偏振相位角（AoP, Angle of Polarization）为：

其中，I0, I45, I90, I135表示0°、45°、90°和135°方向上线偏振光的强度.

三维重建的关键在于获得物体表面准确的法向量，法向量使用其天顶角θ和方位角α表达.对于镜面非导电物体，天顶角与偏振度的关系为：

其中n为物体表面的折射系数.对于金属导电物体，天顶角与偏振度的关系为：

其中n*=n2(1+κ2)，κ为材料的衰减系数.图1为两类物体天顶角与偏振度的仿真曲线.

图 1 偏振度与天顶角仿真曲线

Fig.1 Simulation of degree of polarization and zenith angle

方位角与偏振方位角的关系为：

用n⃗表示物体表面某处的法向量，可表示为：

其中f(x,y)为物体曲面方程，zx(x,y)与zy(x,y)分别为曲面方程对x,y轴上的偏导，即物体表面的梯度信息.对物体表面不同形状，不同曲率位置的法向量进行积分，便能得到物体的相对深度信息，通过深度信息的解算，最终获得物体表面的三维形状.

1.2　偏振三维重建存在的问题

图2以桌球为例，从左到右显示了利用偏振三维重建的过程，首先利用拍摄的四张偏振影像图2（a）求得目标物的偏振度图2（b）和偏振方位角图2（e），其次由偏振度和偏振方位角求得物体表面法向量的天顶角图2（c）和方位角图2（f），最终由获得的法向量图2（d）得到物体表面梯度并重建出三维表面图2（g）.图2中重建算法采用经典的梯度积分算法，利用Neumann边界条件进行泊松重建^[16]，但重建的结果并非是我们期望得到的球形.在利用偏振信息直接进行物体表面重建时，实际上存在一些问题：

图2 桌球偏振三维重建

Fig.2 Shape from polarization of snooker

（1）方位角歧义：偏振方位角解算时φ 与φ+π的正切值是一样的，反正切得到的值域为-π2,π2，而实际需要的方位角值域为-π,π，即偏振获得的方位角部分区域存在π歧义；

（2）图像噪声：实际拍摄受周围环境和相机系统噪声的影响，偏振图像无法避免噪声的存在.由梯度到三维表面是一个积分过程，对图像噪声的积分会引起累加性的误差，目前尚未有完善的积分算法来处理这类问题.

（3）曲面初始偏差：在未知初始边界的条件下，一般的积分算法只能设定或预估初始值，从而导致积分结果与真实值间存在系统性偏差.

2 整体流程与数据预处理

2.1　整体算法处理过程

针对偏振三维重建时存在的问题，本文整体处理的结构图如图3所示.

图3 整体处理结构图

Fig.3 Overall processing flowchart

（1）利用偏振相机获得的偏振图像与深度相机获得的粗糙深度图配准，得到目标物的粗糙法向量及深度图；

（2）利用粗糙深度图求得的目标法向量来纠正偏振图像获得的偏振方位角，得到准确的目标物表面法向量的方位角；

（3）利用已纠正的偏振法向量与粗糙深度图积分融合，获得更高精度的物体三维表面.

图4以杯子为例，显示本文整体算法的处理流程.图4（a）为由深度相机获得的粗糙深度图；图4（b）由偏振相机获得的0°、45°、90°、135°偏振图片；图4（c）为1.1节的处理结果，直接由偏振求解法向量，并用经典的积分算法重建，其结果存在方位角歧义，直接导致形状畸变；图4（d）为经过3.1节方位角去歧义处理的法向量积分重建结果；图4（e）为经过3.2节融合已纠正的偏振法向量信息和粗糙深度信息的三维重建结果.从图4中可以看出，粗糙深度图融入偏振信息后可以获得精确的物体三维表面.

图4 杯子三维重建过程

Fig.4 The 3D reconstruction process of a cup

2.2　粗糙深度图获取

本文利用Astra获得粗糙深度图作为先验信息，Astra是深圳奥比中光科技有限公司开发的一款3D相机，采用红外散斑结构光的方式来获取景深，其拍摄深度范围0.8～8 m、深度图分辨率为640×480@30FPS、精度1 m:±1-3 mm、视场角H:58.4∘, V:45.5∘.

图5（a）为利用Astra拍摄的一组深度图，其中正中央放置着待测的目标物（杯子），背面为挡板.对二维深度数据处理，采用基于视场角的方法可以获得点云数据，具体解算式为：

其中，x,y为投影坐标系中的x,y坐标；X,Y为世界坐标系中的X,Y坐标；Xres,Yres为深度图的分辨率（本文采用640×480）；FOV-h为水平方向视场角；FOV-v为竖直方向视场角；d为深度值.图5（b）为点云图计算结果，图5 （c）为所要提取的目标物点云图，图5 （d）为目标物三维显示.由于深度相机视场较大且拍摄距离需要大于0.8 m，因此无法获取近距离高密度的三维点云数据，只能得到粗糙的目标表面.

图5 Astra 3D相机数据获取

Fig.5 Astra 3D camera data acquisition

2.3　偏振相机标定

本文使用Lucid Phoenix系列偏振相机获取偏振图像，该相机采用覆盖防反射材料来抑制闪光和重影的气隙纳米线栅作为线偏振阵列层，偏振片按四个角度（90°、45°、135°和0°）放置于单个像元上，每四个像元一组作为一个计算单元.相机图像分辨率2448×2048@24FPS、精度为1 m:±1~3 mm.

为了将Lucid偏振相机获取的图片与Astra获取的深度图精确配准，有必要进行相机标定（Camera Calibration），从而获取相机内外参，本节以偏振相机为例来描述相机标定过程.Lucid偏振相机由于添加偏振片阵列光学结构，导致相机模型与原始的小孔成像模型有差异，但由于实际仅选取1/4的单张偏振影像进行配准，从而忽略组合单元的影响.本文采用张正友标定方法^[17]，该方法是一种有物理棋盘格辅助的相机参数检校，其优点是标定简单、精度较高.

对于相机成像，二维像点与三维坐标点之间的转化可表示为：

其中p=[u,v,1]T为二维相机平面像素坐标；P=[X,Y,Z,1]T为三维世界坐标点；s为世界坐标系到图像坐标系的尺度因子；K为相机内参矩阵；R为旋转矩阵；t为平移向量.以平面标定板上左上角的点为世界坐标系的原点，其XOY平面与标定板平面重合，则标定板上点的Z=0，公式（9）可以转化为

其中(u0,v0)为像主点坐标，γ为径向畸变参数，fx,fy将相机的焦距f变换为在x,y方向上像素度量表示，以上五个参数即为所要求解的相机内参.通过标定板角点检测可获得多组方程组，再加上旋转矩阵的约束条件，及对多组解的最大似然估计可获得较高的精度的标定参数，图6为相机内外参标定.

图 6 相机标定

Fig.6 Camera calibration

最终获得偏振相机的内外参矩阵分别为：

2.4　图像配准

同时得到偏振相机与深度传感器参数之后，需要将两类影像数据进行配准，主要原因是两类图像素之间不是一一对应关系.具体可通过图7观测，空间中任一点P在深度图与偏振图上的投影位置，会根据两类传感器之间的姿态与内参的不同而不同，且分别投影在点Pdep与点Ppol上.配准的最终目的是，将深度图上任一点Pdep，对应地赋予给自身对应点Ppol的偏振图像上.解算思路是将深度图中某像素投影到空间中，若该空间点被偏振相机观测到，可以将其投影到偏振图像坐标系中，最终将该点的深度值赋值到偏振图的同名像点上.

图7 深度图与偏振图配准示意图

Fig.7 Alignment of the depth map and polarization map

设偏振相机与深度传感器经过相机标定后的内参分别为Kpol,Κdep，并且在统一的参考坐标系（Coordinate Datum）下的其他外方位元素分别为Rpol,tpol与Rdep,tdep.设点P在偏振相机坐标系下的空间坐标为Ppol，在深度坐标系下的空间坐标为Pdep（均为非齐次坐标），则该空间点在两类影像上的投影坐标分别为:

由于点Ppol, Pdep在空间中为同一点P，所以之间的转换可以通过旋转矩阵R与平移向量t表达，即可以转换深度图像至偏振相机观测的视角下

由式（11）与式（12）可得到

结合相机内参，在图像像素坐标系下，式（13）可表示为：

对相同的目标物，深度相机获得的像素尺寸一般小于偏振相机获得的像素尺寸，需要进一步做插值处理，由此可将两幅不同相机的图像进行配准.

3 偏振信息与粗糙深度图融合

针对偏振三维重建存在的问题：方位角存在π歧义；噪声对法向量积分的误差累计影响；获得曲面的初始偏差.本文利用粗糙的物体深度图为基础，来融合两种图像信息，使获得的物体三维表面具有粗糙深度图的大体轮廓，并具备偏振获取的详细信息.

在此研究中有以下假设，（1）实验在无偏光的环境中进行；（2）物体材质为非导电体，且整个物体表面材质相近；（3）忽略物体与周边环境的多次交叉反射，天顶角的计算主要针对单次反射的结果.

3.1　方位角纠正

在1.3节中有论述过，部分物体表面偏振方位角与实际三维物体的表面法向量的方位角存在π歧义.由这个问题出发，首先采用边缘传播算法获得初步的纠正结果，边缘传播算法假定物体表面整体遵循圆周周期变化，由此预估图像四周的方位角分布如图8所示，从外围选取种子点向里迭代，当所在点的方位角处于该区域给定的方位角范围内时，则认定为准确值，否则利用邻近点插值代替.

图8 方位角的分布范围与数值大小

Fig.8 Azimuth distribution range and numerical value

该算法对于非对称的复杂物体容易出现误判，因此本文在此基础上进一步利用粗糙深度图的方位角来纠正，具体计算方式为：

其中φc0，φd分为边缘传播纠正、深度图计算得到的方位角，α为0,1间的系数，α取值大则会误判噪声点，取值小则会影响整体结果准确性，本文实验选用3/4.

3.2　物体表面重建

由法向量重建出物体三维表面是一个积分过程，除了需要边界条件限制，图像里的噪声、畸变等都会导致积分结果的错误，如何找到一种鲁棒的积分算法一直是三维重建里的难题.本文利用粗糙的深度信息作为整理框架，并融入偏振细节信息.为此首先定义两个最小化的能量项^[18]：点误差Ep和法向量误差En.

点误差定义为所有的优化三维点和深度图测量三维点平方距离的总和：

其中Pi为优化后第i个点的位置，Pid为深度图中对应的i点.当Ep达到最小值时，优化后的三维点与深度图的三维点一致.实际上三维点在x, y的方向上是固定的，因此式（16）可以转换为z轴深度图上表示：

其中fx, fy为相机在像素点上的焦距.

法向量表示垂直于空间曲面点的任一切线方向，因此可利用切线与法线垂直的关系来表述法线误差：

其中Nipc为从偏振图像获得，已方位角去奇异的法向量，图像x, y方向的切线向量可由偏微分方程表示：

最终的优化方程可以表述为：

其中k为0,1 间的参数，决定深度信息与偏振信息三维重建的比重，当k取值为0时，优化的结果等同于纠正后的偏振重建结果，当k取值为1时，优化的结果等同于深度图的结果.而且该方程组为超定的线性系统，可用最小二乘法求得最优解.

4 实验结果与分析

针对低纹理乃至无纹理的目标物，本文选用合成树脂光滑表面（桌球）、陶瓷磨砂表面（杯子）、陶瓷光滑表面（两个不同形状的花瓶）作为验证实验，并利用扫描仪得到的结果作为真实值来评估实验结果的准确性.

图9为法向量处理的结果.偏振获取的信息处理最终得到物体表面的法向量，而法向量的重建结果受目前积分算法不完善的限制，因此通常也以法向量的结果图来表示物体处理的优劣.图中第一列为实物图；第二列为扫描仪获得的结果作为真值；第三列为本文算法对偏振获得的方位角去奇异处理后的结果；第四列为文献[19]算法得到的结果；第五列为文献[20]算法得到的结果；第六列为文献[21]算法得到的结果.结果与真值的误差以MAE（Mean Absolute Error）表示，即真值法向量与求解法向量之间夹角绝对值的平均值，单位为度.从结果可以看出，在物体垂直观测视角的切面上，磨砂、光滑表面均会出现明显的噪声，这是由于在这个区域偏振度较小，噪声的影响更明显.作为对比的三类结果中，文献[19]获得的结果与真值较为接近，但其实验依赖于对光源的确定，而且不同的光源参数值结果相差很大.从MAE值比较可以看出本文算法获得的法向量与真值较为接近.

图9 法向量结果

Fig.9 Normal vector results

图10为物体表面重建的结果.图中第一列为扫描仪获得的结果作为真值；第二列为Astra获得的粗糙深度图；第三列多视角处理的结果，利用商业软件Agisoft PhotoScan处理同一物体多张影像，在最高质量的模式下配准图像并生成密集点云，其中桌球选取10张配准0张无三维结果，杯子选取10张配准4张，花瓶1选取10张配准5张，花瓶2选取10张配准7张；第四列为融合深度图的重建结果，即本文重建的最终结果；第五列为文献[16]经典的梯度积分重建算法得到的结果；第六列为文献[22] Shapelet算法得到的结果.结果误差以MAE表示重建点深度与真实点深度误差绝对值的平均值，单位为mm.由于积分算法在未设定边界初始值时，重建结果会存在一定初始偏差，可能会影响MAE的评判结果，因此本文引入相关系数r（Correlation Coefficient）来辅助评判结果与真实值的相似程度.MAE越高说明结果整体与真值相差较大，r越高说明重建形状与真值越接近.文中选取的低纹理物体，多视角结果并不理想，多视角三维重建依赖于目标物的纹理特征，对于低纹理乃至无纹理物体提取的特征点少，无法生成准确密集的目标点云.文献[16]算法的重建结果存在较大畸变，尤其在桌球上明显有错；文献[22] Shapelet算法能避免方位角歧义问题，但对噪声敏感，中间部分容易出现毛刺，当调整天顶角参数时，会引起部分过饱和；深度相机只能获取粗糙深度图，尤其在花瓶2里，丝带纹理基本没有，只存在大体轮廓，本文算法重建结果依赖于偏振信息获取物体丝带立体纹理.从三维形状和误差值比较可以看出本文算法获得的物体三维表面与真值较为接近.

图10 表面重建结果

Fig.10 Surface reconstruction results

5 结论

以粗糙深度图作为先验信息，可以有效地处理偏振三维重建时出现的一些问题，包括方位角π歧义，积分受图像噪声影响，重建初始偏差，最终可以获得精度较高的物体三维表面.使得利用偏振手段，能很好地处理表面光滑、纹理单一的物体表面三维重建.

在天顶角求解过程中，垂直于观测方向的区域，偏振度较低，且容易受噪声影响，所求解的天顶角存在较大的误差，这部分的误差会直接影响物体的形状，后续可以考虑引入融合算法等提升重建的精度；对于结构比较复杂的物体，粗糙深度图与偏振图像的精确匹配直接影响重建的精度，需要有更高精度的初始深度图，以及优化的匹配算法，从而能处理结构更为复杂的物体.

References