小波变换轮廓术的测量系统的光路原理和傅里叶变换轮廓术相同^{[2, 4]}。成像装置获得的变形条纹图一般可表示为：

式中，$a(x, y)$为背景分量；$b(x, y)$为条纹可见度；${f_0}$为空间载波频率；$\varphi (x, y)$为物体高度引起的调制相位；$c(x, y)$为随机噪声。其中，${f_0}$必须满足下式才能恢复相位^[11]：

条纹图一维连续小波变换为^[12]：

式中，$a$为尺度因子；$b$为平移因子；$f(x)$为取条纹图一行的信号；$\psi (x)$为母小波；“$*$”表示共轭。

Morlet小波在空域和频域中具有良好的局部性，文献[13]通过对几种小波实验对比，指出复Morlet相对于其他小波在相位恢复中更具优势，故本文选用复Morlet作为母小波，其表示如下：

式中，${f_b}$为带宽参数；${f_c}$为小波中心频率参数。

通过使用复Morlet小波对条纹图信号进行连续小波变换，可得到小波变换系数$W(a, b)$，其实部记为${\rm{RE}}[W(a, b)]$，虚部为${\rm{IM}}[W(a, b)]$，则幅值和相位分别为^{[4, 13]}：

由式(3)可以看出，小波变换系数$W(a, b)$为信号和小波序列函数内积的结果，它反映了信号和小波序列函数的相似程度。因此，当信号的局部频率和相应尺度的小波函数振荡频率相近或者相同时，$W(a, b)$较大，即相应的幅值$A(a, b)$较大。小波脊即为在小波平移方向上，不同位置的幅值最大值的连线。脊处对应的尺度即为最佳伸缩尺度，对应的相位即为包含物体高度信息的调制相位。

条纹图受到噪声影响后，小波变换系数幅值会产生多个局部极大值点，若直接采用模极大值提取相位，易造成脊点的误选，导致小波脊定位不准确。根据条纹图瞬时频率的连续特性，文献[9]提出了基于代价函数的小波脊提取算法，该算法通过计算最小代价值来确定小波脊。

代价函数定义为：

式中，$b$为平移因子；$\phi (b)$为位置$b$处的尺度因子；$S[\phi (b), b]$表示在尺度为$\phi (b)$，位置为$b$处的小波变换系数幅值；${C_0}$和${C_1}$分别表示小波变换系数幅值和尺度参数曲线梯度的权值系数。

在实际应用中，处理的对象为数字图像，位置是离散化的，因此需要对式(7)进行离散化，虽然原理上可计算任意位置的代价值，但实际上这是不现实也是不必要的。式(7)离散化形式为：

式中，$L$表示小波变换系数幅值矩阵的宽度；$b$表示平移因子。

由式(8)可以看出，其选取的小波脊代价函数较小且较为光滑，避免了因噪声产生较大幅值造成的脊线不光滑。

如式(8)所示，在一般研究中，都默认设置${C_0}$和${C_1}$的值为1。小波变换系数幅值$S[\phi (b), b]$通常在0~1之间。小波变换尺度$\phi (b)$无论选择步长为1还是步长为指数增长形式^[10]，$\phi (b) - \phi (b - 1)$一般大于1。因此通常有：

由式(8)选择的小波脊过于偏向平滑，受到噪声影响后，易误选使脊线更加平滑的噪点作为脊点。因此，本文提出一种新的结合空域噪声信息和Feature Scaling方法的代价函数，其定义为：

式中，${\rm{Max}}s$和${\rm{Min}}s$分别表示候选脊点对应小波变换系数幅值矩阵除第一列；${\left| {S[\phi (b), b]} \right|^2}$的最大值和最小值；${\rm{Max}}b$和${\rm{Min}}b$分别表示候选脊点对应的尺度因子矩阵梯度${\left| {\phi (b) - \phi (b - 1)} \right|^2}$的最大值和最小值；$K$为$\lambda $的调整因子，一般为${10^n}$，n为整数，使得$K\lambda \geqslant 1$；$\lambda $为噪声程度系数，具体判断标准可根据实际情况定义使用，本文中取噪声数与总像素比，其表达式为：

式中，${\rm{Count}}$为检测的噪声总数；${\rm{Pixel}}\_{\rm{num}}$为变形条纹的总像素数。

由式(10)可以看出，首先使用Feature Scaling方法对小波变换系数幅值和尺度参数曲线梯度两项特征因子进行特征缩放，权衡两项特征因子的影响比重，然后利用空域噪声信息对尺度参数曲线梯度因子进行补偿。在较低噪声情况下，小波脊受噪声的影响较小，对尺度参数曲线梯度因子补偿较小，使两项因子权重等阶，既使用模极大值提取小波脊，又利用小波脊线的光滑性过滤噪声。随着噪声的增加，噪声会产生较多属于极值的小波变换系数幅值，更易干扰小波脊的正确提取，此时增加对尺度参数曲线梯度因子的补偿，加大该因子的权重，即增强利用小波脊线光滑特性滤除噪声，从而正确定位小波脊。

本文结合条纹图灰度的连续特性，选取噪声检测中常用的基于梯度的判别法来判断噪点^[15]，从而计算$\lambda $值。噪点分析示意图如图 1所示，相关公式为：

式中，$L(c)$为噪点梯度和；$c$表示检测点${x_0}$的邻接区域；${x_0}(c)$表示检测点${x_0}$的灰度值；${x_k}(c)$表示${x_0}$邻接点的灰度值。

如图 1a所示，若$L(c)$大于某个阈值，可判定${x_0}$为噪点，对于边缘像素点，如图 1b所示，$L(c)$远小于设定阈值。

图  1  噪点分析示意图

获取$\lambda $步骤如下：

1) 提取变形条纹图的一行。

2) 读取第$j$个像素点，其中$1 \leqslant j \leqslant n - 2$，$n$为条纹图一行的长度。设读取的像素点为A，其灰度值为${\rm{pixel}}[i][j]$。A后第2个像素点$B$的灰度值为${\rm{pixel}}[i][j + 2]$，计算A与B的灰度差值$d = $ $abs({\rm{pixel}}[i][j] - {\rm{pixel}}[i][j + 2])$。

3) 分别计算A和B的噪点梯度和$L1$、$L2$。

4) 判断是否$d < d\_{\rm{leftthreshold}}$，若否，执行步骤6)，若是，执行步骤5)。

5) 判断是否$L1 > {\rm{Threshold}}$且$L2 > {\rm{Threshold}}$，若是，则A、B都为噪点，噪点数${\rm{Count}}$加2。计算第$j + 1$个像素点的梯度和$L3$，若$L3 > {\rm{Threshold}}$，则为噪点，${\rm{Count}}$加1。$j$加2，返回执行步骤2)。

6) 判断是否$d > d\_{\rm{rightthreshold}}$，若是，则执行步骤7)，否则$j$加1，返回执行步骤2)。

7) 判断是否$L1 > L2$且$L1 > {\rm{Threshold}}$，若是，则$A$为噪点，${\rm{Count}}$加1，$j$加2，返回执行步骤2)。若否，执行步骤8)。

8) 判断是否$L1 \leqslant L2$且$L2 > {\rm{Threshold}}$，若是，则$B$为噪点，${\rm{Count}}$加1，$j$加2，返回执行步骤2)。

9) 重复上述步骤，直到检测完变形条纹图最后一行。计算$\lambda {\rm{ = }}\frac{{{\rm{Count}}}}{{{\rm{Pixel}}\_{\rm{num}}}}$。

本文中，阈值依据基准光栅来确定，$d\_{\rm{leftthreshold}}$和$d\_{\rm{rightthreshold}}$分别选取基准光栅中$d$的最小值和最大值。设$L{(c)_{\max }}$为基准光栅中噪点梯度和的最大值，由于受到物体高度的调制，变形光栅中的噪点梯度和必然大于$L{(c)_{\max }}$，因此本文设${\rm{Threshold}}$为$2L{(c)_{\max }}$。

小波脊提取步骤如下：

1) 分别计算$\lambda $值和$K$值。

2) 提取变形条纹图的一行对其进行连续小波变换，获得小波变换系数矩阵${\rm{Coefficient}}(m, n){\rm{ }}$，由式(5)可求得幅值矩阵${\rm{Modulus}}(m, n)$。其中，$m$为尺度$a$的维数，$n$为条纹图一行的长度。

3) 沿平移因子方向对${\rm{Modulus}}(m, n)$第$i$列做尺度-模值曲线，求取曲线局部极大值点和拐点，并将大于均值的点作为候选脊点。

4) 重复步骤3)直到$i = n$，获得候选脊点矩阵$\boldsymbol{\rm{CandidateRidgePoints}}$，进而获得对应的小波变换幅值矩阵$\boldsymbol{\rm{Can\_Magnitudes}}$和尺度因子矩阵$\boldsymbol{\rm{Can\_Scales}}$。

5) 根据$\boldsymbol{\rm{Can\_Magnitudes}}$和$\boldsymbol{\rm{Can\_Scales}}$分别计算${\rm{Min}}s$、${\rm{Max}}s$、${\rm{Min}}b$、${\rm{Max}}b$。

6) 设置和$\boldsymbol{\rm{CandidateRidgePoints}}$相同大小的代价值矩阵$\boldsymbol{\rm{Costs}}$，并将第一列代价值设置为0。

7) 从第二列开始，分别选取$\boldsymbol{\rm{CandidateRidgePoints}}$第$j$列每个候选脊点$i$，根据式(10)代价函数以及$\boldsymbol{\rm{Can\_Magnitudes}}$和$\boldsymbol{\rm{Can\_Scales}}$，依次求得列$j - 1$各个候选脊点$m$到选取脊点$i$的代价函数值$\operatorname{Cos} {\rm{t}}(m)$，并求最小代价函数值${\rm{Cost}}(i, j) = \min \{ {\rm{Cost}}(m, j - 1){\rm{ + Cost}}(m)\} $，将${\rm{Cost}}(i, j)$放入代价值矩阵$\boldsymbol{\rm{Costs}}$对应位置中。

8) 重复步骤7)到最后一列处理完毕，选取最后一列最小值的点作为脊点，并反推出所有脊点，得到小波脊。

9) 重复步骤2)~8)，得到整个变形条纹图的小波脊，并根据式(6)求解出小波脊相位。

为验证本文所提算法的有效性和鲁棒性，本文使用两种不同形状的物体来进行模拟实验验证。并对变形条纹图分别添加方差为0.02、0.04、0.06、0.08、0.1的斑点噪声来验证所提方法的抗噪性。本文使用常用的评价指标-均方根误差RMSE(root mean square error)来计算算法在不同噪声强度的重建误差，RMSE的表达式为：

式中，${h_f}(x, y)$和${h_r}(x, y)$分别表示模拟物体标准高度信息和重建物体高度信息；$m \times n$为条纹图像尺寸大小。

模拟的被测物体由Matlab函数库中的peaks函数生成。物体的高度设为$h = 3.7 * {\rm{peaks}}$，模拟条纹图大小为512×512 pixels，其中标准正弦光栅周期为16 pixel，模拟物体和在方差为0.1的噪声下变形条纹图分别如图 2、图 3所示。

图  2  模拟物体

图  3  含噪变形条纹图

首先，本文选取中值滤波技术对含噪的变形光栅图进行滤波处理，设定其滤波窗口为3*3。为了便于对比分析，分别采用了模极大值法、原代价函数法和本文算法对模拟变形条纹图进行小波脊相位提取。并对提取的截断相位都采用相同方法进行相位展开。选取的母小波为复Morlet小波，带宽参数为0.5，小波中心频率参数为1，尺度因子为1~64，步长为1。对于每种强度的噪声，分别使用3种方法模拟10次。在不同噪声情况下，各算法获得的10次RMSE均值如表 1所示，RMSE曲线如图 6所示，原代价函数法和本文算法在方差为0.1的噪声下重建形状分别如图 4、图 5所示。

图  6  各算法在不同噪声下模拟10次RMSE曲线图

图  4  原代价函数法

图  5  本文算法

可以看出，随着噪声的增加，3种方法RMSE值都在增加，本文算法RMSE较小，明显优于模极大值法和原代价函数法，具有更好的抗噪能力。本文算法在噪声低的情况下明显优于其他两种方法，随着噪声的增加，模极大值法波动增大，误差极不稳定，原代价函数法和文本算法则波动较小，较为稳定。噪声增加后，本文算法与原代价函数法差距逐渐减小，但本文算法仍优于其他两种方法。

为进一步验证本文算法的有效性和鲁棒性，使用高度变化较大的圆锥体作为验证对象，模拟的圆锥体的表达式为：

模拟条纹图大小为250×250 pixels，标准正弦光栅周期为16 pixel，模拟物体和在方差为0.1的噪声下变形条纹图分别如图 7、图 8所示。

图  7  模拟物体

图  8  含噪变形条纹图

其他步骤和参数设置与peaks模拟相同，3种方法在不同噪声情况下模拟10次，均方根误差RMSE均值如表 2所示，RMSE曲线如图 11所示，原代价函数法和本文算法在方差为0.1的噪声下重建形状分别如图 9、图 10所示。

图  11  各算法在不同噪声下模拟10次RMSE曲线图

图  9  原代价函数法

图  10  本文算法

可以看出，本文算法明显优于其他两种方法。在较低噪声情况下，模极大值法波动较大，原代价函数法和本文算法波动较小，较为稳定，随着噪声的增加，原代价函数出现波动情况，本文算法仍相对稳定，这是由于圆锥物体高度变化较大，原代价函数提取的小波脊趋于过于平滑，造成小波脊提取不准确。

本文提出了一种新的代价函数提取小波脊的算法。相比原代价函数，本文首先使用Feature Scaling方法对小波变换系数幅值和尺度参数曲线梯度的权值进行调整，然后结合空域噪声信息，对尺度参数曲线梯度进行补偿。避免了在低噪声情况下提取小波脊时，脊线过于趋向平滑，以及在高噪声情况下，脊线波动较大造成小波脊选取不准确的问题。实验表明，无论对于高度变化平滑还是陡峭的物体，本文算法具有良好的抗噪性、鲁棒性以及稳健性。