张申华,杨延西,秦峤孟
(1. 西安理工大学 自动化学院,陕西 西安 710048;2. 安康学院 电子与信息工程学院,陕西 安康 725000;3. 陕西省复杂系统控制与智能信息处理重点实验室,陕西 西安 710048)
由于具有非接触、高精度和快速的优点,基于正弦光栅条纹投影的相位轮廓术目前被广泛研究。正弦光栅被数字光处理(Digital Light Procession, DLP)投影仪投影至被测物体的轮廓上,经物体外形调制后发生形变,并被工业相机采集[1-8]。计算出捕获条纹的相位,借助标定出的参数,最终可以计算出物体的三维外形轮廓尺寸。然而在测量过程中,受到环境及工业相机自身的影响,采集到的光栅图像存在噪声。噪声的存在降低了图像的质量,影响后续相位的高精度提取,并最终导致三维测量准确度降低。因此,有必要对噪声进行处理。在实际应用中,由于采集到的光栅图像噪声未知,从而给噪声处理带来了困难。
虽然在数字图像处理领域很多研究人员对去噪问题进行了大量研究[5-7],对于条纹图像噪声的处理,主要采用中值滤波的方法[8]和偏微分方程(Partial Differential Equation, PDE)方法[9-11]。中值滤波方法实现简单,但是对噪声的抑制效果有限。基于偏微分方程的滤波方法是利用偏微分方程进行数值运算,逐个像素修改灰度值,以达到滤除噪声的目的。为了滤除噪声,Wang提出了一种方向性的相干增强扩散滤波方法[9],该方法沿条纹水平方向及垂直方向滤波。文献[10]提出了一种基于二阶偏微分方程模型的方向性滤波方法,根据电子散斑干涉条纹的方向滤波。Zhou对文献[10]中的模型进行了改进,提出了一种自适应的方向性偏微分方程模型[11]。总的来说,基于PDE扩散的滤波方法存在以下不足:(1)该类方法基于扩散的思想,逐点处理每个像素点,而相位轮廓术一般需要处理多帧光栅图像,噪声处理过程耗时;(2)数值计算中的迭代次数、时间步长等参数通常需要反复实验才能确定。
基于噪声水平估计的方法,是盲噪声处理的重要研究方向[12-18]。该方法一般利用高斯分布模型对未知噪声进行建模,并获取噪声的方差值。Donoho利用离散小波变换对信号进行三层分解,并取对角线系数的中值作为噪声方差的估计值[12]。基于主成分分析(Principal Component Analysis,PCA)的方法[13-14]利用单幅图像生成样本块库,并通过样本块库估算出噪声方差。这种方法在高光谱图像去噪[15]、视频噪声估计[16]、中子图像去噪[17]及特征提取[19]等方面得到了应用。但是基于PCA的方法估计出来的噪声水平常常偏低[20],导致利用估计值进行去噪的效果不够理想。
在现有的研究中,对于图像中的未知噪声,一般是进行高斯建模[20-22],基于此,本文借鉴现有文献的处理方法,对光栅图像中的噪声进行高斯建模,并引入PCA方法对噪声水平进行估计,并依据估计结果完成噪声处理。首先,针对PCA估计出的方差水平偏低的问题,利用残差模型,提出一种真值图像和噪声图像分离方法,从含噪光栅图像中分离出噪声图像,然后对噪声图像进行PCA估计。其次,提出了一种基于相图的噪声滤除方法。该方法利用估计出的噪声方差值,在展开的相位图上进行高斯滤波,避免了对多帧光栅条纹分别进行滤波处理,降低了处理的数据量,提高了噪声处理的时效性。
在光栅投影测量系统中,DLP投影仪投影到被测物体上的光栅条纹被工业相机采集,用符号f表示。根据残差模型,f中包含噪声n和未被污染真值条纹图像r,用公式表示为
式(1)中的噪声n为 均值为零、方差为的高斯白噪声。。根据Liu的分析[14],在某一向量U的方向上有
为了滤除条纹中的噪声,需要首先估算出
其中Cf表示f的协方差矩阵,Cr表示r的协方差矩阵,λmin表示求取协方差矩阵的最小特征值。
在式(2)中,对于细节信息不丰富的图像,在信息冗余特性的假设下,r的数据主要分布在有限主成分分量上。在此前提下,
由此,可以估算出噪声的方差值
研究组成员利用PCA方法对仿真光栅条纹进行了噪声方差估计,结果如图1所示。由图1(b)中的比值曲线可以看出,PCA方法估计出的方差普遍偏小(图1中 σ表示噪声方差的真值,表示估计值)。
图1 PCA方法对仿真噪声的估计Fig. 1 The estimation of simulation noise by PCA method
若选择数据的协方差矩阵大小为9×9,其各主成分占比如图2所示,由图中占比可以看出,后7个主成分占比小且极为接近,若仅选择占比最小的主成分对应的特征值作为噪声方差,实际上人为降低了噪声的成分,这是PCA方法估计结果偏小的原因。
图2 主成分分量占比Fig. 2 Proportion of principal component
针对PCA估计方法的不足,本文对噪声估计方法进行改进。对式(1)中的f进行特征值分解,可以得到
式中,U表示特征值矩阵 Λ对应的特征向量。符号T表示转置操作。Λ为对角矩阵,表示为
其中 λ1,···,λM表示特征值且以降序排列,其大小代表了各主成分的能量贡献度,M表示f的秩。
令各主成分的能量贡献率表示为
若令前P项特征值(P<M)对应的向量为RP,则真值光栅图像可表示为
根据式(1),噪声图像可以用残差模型表示为
根据文献[12]的证明,n的协方差矩阵最大特征值收敛于,也就是最终的噪声方差估计公式为:
为了实现式(9)中真值光栅图像与噪声图像的分离,需要合理确定P值。假设工业相机的信噪比为ddB,也就是信号的功率Wo与噪声功率Wn的比例关系为
根据式(11)所表示的比例关系,且所有主成分的能量贡献值为,又知真值图像的能量贡献与噪声图像的能量贡献比值与Wo/Wn相等。由此,噪声部分的能量贡献值可以表示为
若某一成分的特征值λi满足
则可以将该主成分分量划归到真实光栅图像r,否则,该主成分分量划归到噪声图像n。
根据式(13)将满足条件的最大i值,设定为P值,即可实现式(9)中噪声图像与真值图像的分离。
为了降低数据处理量,减少处理时间,本文提出了基于相位图的高斯去噪方法。在相移法测量中,通常需要投影和采集多幅光栅图像。以4步相移法为例,投影的4幅光栅条纹为
其中A表示平均分量,B表示调制度,ϕ表示相位。
若工业相机采集到多帧光栅图像后,再对每一幅图像进行去噪,无疑会大大增加要处理的数据量和时间。为此,本文将光栅图像上的噪声处理转换到相图上进行。投影的光栅条纹的光栅频率分别为
本文采用三频外差法[23]展开包裹相位,需要, 其包裹相位分别用 ψ1,ψ2,ψ3表示。通常,取s=64, 首先由ψ1和 ψ2做外差,得到频率为的 包裹相位 ψ12,再由 ψ2和 ψ3做外差,得到频率为的包裹相位 ψ23,然后由 ψ12和 ψ23做外差,得到频率为1的相位ψ123,最后,利用式(15)和式(16),通过中间展开相位 φ■s+1即 可得到ψ1的最终展开相位ϕ。
在式(15)和式(16)中,相位ψ1,ψ12,ψ123取值在(−π,π)之 间,而和 ϕ的取值范围一般远大于该取值,因此,可以对式(15)和式(16)做如下近似:
在仅考虑零均值的高斯噪声影响下,时域上光栅条纹的噪声方差 σn和 展开相位图 ϕ的方差σ1及 相移步长N、调制度B之间关系为[24]
根据方差的性质,依据式(17)和式(18)所表示的比例关系,相位和 ψ123的方差分别为
图3为本文搭建的光栅投影测量平台结构图,其主要由投影仪、工业相机及便携式计算机组成,主要参数如表1所示。实验首先投影12帧
正弦光栅条纹图像,并由同步信号触发控制工业相机完成光栅图像的采集。利用采集的光栅条纹图像,对本文所提方法和对比方法进行了测试。
图3 测量系统结构图Fig. 3 Framework of measuring system
表1 仪器设备型号和参数Tab. 1 The instrument types and parameters
在测量平台上采集到的12帧光栅图像中选取一帧,如图4(彩图见期刊电子版)所示。利用本文提出的方法和对比方法进行噪声方差估计后,分别进行噪声处理。各种对比方法的噪声估计值如表2所示。
图4 采集的光栅图像Fig. 4 Captured grating image
表2 各种方法的噪声方差估计值Tab. 2 The estimated noise variances of different methods
为了验证本文方法的估计效果,利用表2的估计值对图4的光栅图像进行噪声滤除,效果如图5(彩图见期刊电子版)所示。对12帧光栅图像采用文献[10]的方法分别进行方向性PDE滤波后,提取出的相位图如图5(a)所示。实验中时间步长设置为 ∆t=0.1s, 迭代次数N=75(实验参数非最优设置)。对光栅图像采用中值滤波后,利用相移法得到的相位图如图5(b)所示,实验中滤波模板大小设置为3×3。首先采用文献[12]提出的方法估计出噪声方差,然后,利用估计值对光栅图像进行高斯滤波处理,处理后的相位图如图5(c)所示,实验中先采用db2小波对光栅图像进行三层分解,然后求取第三层分量的对角系数,计算出其中间值与0.674 5的比值,最后将计算出的比值作为噪声方差的估计值。采用文献[13]中提出的PCA方法对单帧光栅图像估计出噪声方差后,再利用估计值对采集到的每帧光栅图像进行高斯滤波,结果如图5(d)所示。利用本文方法先估计出单张条纹图像的方差值,再进行相图滤波后提取的相位结果如图5(e)所示。图5(f)是投影相移步长为20的正弦条纹,然后利用相移法和三频外差法获得的相位图。
对于图4中红色虚线标定的一行区域,以二十步相移法获取的相位为基准,采用其他方法计算出的相位与基准相位取差值,相位差值曲线如图6(彩图见期刊电子版)所示。由相位差曲线可知,和其他对比方法相比,本文方法处理后的相位更加逼近真实相位。本文将二十步相移法获取的相位作为测量体的真值相位并将其作为基准相位。这样设置的依据如下:根据式(19),增大相移步数N能够有效降低相位差的方差值,并且在现有的文献研究中,一般采用二十步相移法[25-26]或四十八步相移法[27]作为测量体的真值,因此,本文采用同样的处理方法获取测量真值。
为了定量评估本文所提方法的性能,定义最大相位误差(Maximum Phase Error, MPE)和均方根误差(Root Mean Square Error, RMSE)作为评价指标,其中 MPE表征相位值偏离真值的最大幅度, RMSE表示相位和真值相位的逼近程度,定义分别为:
式中 ϕ′表示采用各种方法处理后的相位,表示被测体的真值相位,L表示相位点个数。
图5 几种不同方法的滤波效果Fig. 5 Filtering effect obtained by different methods
图6 各种方法的相位差Fig. 6 Phase errors obtained by different methods
任意选取测量体的第150行相位,计算出各种方法的M PE、R MSE及执行时间,结果如表3和表4所示。由表3和表4的数据可以看出,在最短的运行时间内,本文方法的 RMSE值较中值滤波、PDE、小波和传统PCA相比,分别下降了38.6%、88.5%、42.6%和43.1%。表明本文方法的相位更加接近测量体的真值相位,对噪声抑制效果最好。需要说明的是,受采集角度和光源等因素影响,光栅图像存在阴影区域,导致区域内的相位信息丢失。对阴影区域采用PDE扩散滤波后,错误的像素信息会经过扩散传导至正常的光栅区域,导致测量结果发生错误。如表3所示,PDE方法的相位最大偏离幅度达9.234 3 rad,因此,当测量体轮廓不连续时,PDE方法噪声处理效果不好,同时,耗时也最久。对于最大相位误差值,中值滤波方法的 MPE 值为0.385 6 rad,小于本文方法的0.700 3 rad。这是由于测量过程中光栅图像存在孤立噪声点,而中值滤波在处理这类噪声时的性能要优于本文方法,但是中值滤波的相位平均运算会导致在物体轮廓不连续区域产生明显的错误相位(如图5(b)中面具左眼区域)。
表3 几种算法的量化指标对比Tab. 3 Comparison of quantitative indicators for different methods
表4 几种算法的执行时间对比Tab. 4 Comparison of computing time for different methods (s)
由相移法提取相位,一般通过求反正切函数计算出包裹相位,且至少需要3帧光栅图像才能计算出包裹相位。本文方法在估计噪声方差时,为了提高时效性,任意选取其中的一帧光栅图像估计方差值。这是由于在测量平台上,光栅图像在数秒内甚至更短的时间内完成投影和采集,噪声方差这一统计值在较短的时间内能够维持稳定。
光栅条纹投影测量过程中,由于采集到的光栅条纹图像存在噪声,且噪声未知,从而导致提取的相位精度降低,甚至出现错误。针对这一情况,本文提出了一种快速盲去噪方法。首先,利用残差模型和相机参数,提出了一种光栅图像的噪声和真值图像分离方法,并引入PCA对噪声图像进行噪声水平估计,克服了传统PCA方法估计结果偏低的问题。为了减少噪声处理的数据量,提出了基于相位图的高斯去噪方法,利用式(19)~式(21)所描述的时域噪声方差与相位图方差间的关系,将在时域上的光栅图像去噪转换至在相位图上进行。转换后,需要处理的光栅图像帧数减少到3帧,处理数据量大为减少,处理时间变短。测试实验结果表明,和对比方法相比,本文方法可在最短的执行时间内,使光栅图像相位的RMSE值最高下降88.5%,所提方法对噪声的抑制最为明显,处理后的相位图更加接近真值相位。
需要指出的是,本文方法将随机选取的一帧光栅条纹图像的噪声方差估计值作为所有采集图片的噪声方差值。这样处理的前提条件是:在光栅条纹图像投影测量过程中,采集到的每帧光栅图像的噪声方差值应基本保持稳定,这也是本文方法有效的条件。