CN102628676A - 一种光学三维测量中的自适应窗口傅里叶相位提取法 - Google Patents

Abstract

一种光学三维测量中的自适应窗口傅里叶相位提取法，是对灰度正弦光栅条纹图进行处理：(1)将一行条纹信号进行经验模态分解得到一系列本征模函数。(2)计算每个本征模函数的瞬时频率及边际谱，通过分析边际谱来选择能够准确描述条纹信号变化的瞬时频率。(3)确定条纹信号的背景分量。(4)用所确定的瞬时频率自适应定位局部平稳区域并据此计算高斯窗口的尺度因子。(5)将原信号去除背景分量后，对其进行自适应窗口傅里叶变换得到包裹相位分布。逐行进行上述步骤可得条纹图的全局包裹相位，展开后得到绝对相位并转换为高度信息，从而实现对物体的三维测量。本发明自适应性强，对面形复杂或含有陡峭边缘物体的测量精度有较大提高。

Description

一种光学三维测量中的自适应窗口傅里叶相位提取法

技术领域

本发明属于三维信息重构的领域。基于灰度正弦光栅投影，采用希尔伯特-黄变换法分析变形条纹图像并求得条纹的瞬时频率及背景分量，用所求的瞬时频率自适应计算窗口傅里叶变换的窗口尺度因子，通过对去除背景的条纹信号直接做傅里叶变换或自适应窗口傅里叶变换从而得到精确的包裹相位，这种方法有利于提高对形貌复杂物体或携带陡峭边缘物体的测量精度。

背景技术

光学三维测量技术能够准确获得物体的三维面形数据，可用于三维模型重建、物体表面轮廓测量、工业环境中的尺寸和形位参数的检测等，因此它在虚拟现实、影视特技、医学整形和美容、工业产品的外观设计、艺术雕塑和文物保护等领域都有广阔的应用前景。

光栅投影法是一种重要的三维测量技术，通过向物体表面投射正弦光栅，将物体的高度信息以相位的形式隐含在光栅中，利用CCD获得物体表面的光栅条纹图像，并使用特定算法对条纹图像进行处理，提取其中的相位，从而建立物体的三维信息。

常用的求解条纹图像相位的方法有基于时间域的方法和基于变换域的方法。基于变换域的方法由于只需采集一幅变形条纹图即可完成相位的测量，有利于实现物体的动态测量，因此而被广泛研究和应用，其中包括傅里叶变换法、小波变换法、S变换法、窗口傅里叶变换法等。傅里叶变换是一种全局的信号分析工具，不能提取局部信号的特征，只适用于平稳信号。近年来，各种时频分析技术被广泛研究以期能够准确获取条纹图像的细节相位信息。连续小波变换由于具有多分辨率分析的特性被引入光学三维测量领域，通过检测小波变换最大脊处的相位，得到变形条纹图的包裹相位。但是这种方法有个前提条件，即被检测的相位必须是近似线性变化且变化缓慢，否则该方法的理论推导是不成立的。同时，小波母函数及相关参数的选取没有成熟的理论依据，也为小波变换轮廓术的广泛应用带来了限制。S变换与小波变换求解条纹图包裹相位的原理高度相似，故它也不可避免的受到与小波变换同样的条件限制。窗口傅里叶变换也具有较好的时频分析性能，但窗口尺寸即窗口尺度因子的自适应选取一直是研究的热点。现有的选取窗口尺度因子的方法大都是先通过连续小波变换或S变换检测最大脊处的尺度因子，将该尺度因子直接作为窗口傅里叶的尺度因子，或将该尺度因子取倒数作为条纹信号的瞬时频率，再通过相关算法计算窗口的尺度因子。这些方法需要事先设定基函数及经验值，因此自适应性不好，同时也同样受到被检测相位前提条件的限制，使求解的尺度因子或瞬时频率过平滑化，而不能很好的描述条纹信号的局部特征，从而导致这些方法只能测量表面变化相对平滑的物体，而对表面变化复杂或含有陡峭边缘物体的测量会受到较大限制。

发明内容

技术问题：针对光学三维测量中窗口傅里叶变换求取条纹图包裹相位的精确性和自适应性等问题，本发明在自适应进行窗口傅里叶变换及提高包裹相位的求解精度等方面给出了解决办法。本方法借助希尔伯特-黄变换，精确且自适应地求解能准确描述条纹图变化情况的瞬时频率，从而自适应计算窗口傅里叶变换的窗口尺度，同时，在不需要额外计算的情况下可以有效去除条纹图的背景分量从而大大减少窗口傅里叶变换处理过程中零频频谱的干扰。本方法能大大提高对复杂面形或含有陡峭边缘物体的测量精度。

技术方案：一种光学三维测量中的自适应窗口傅里叶相位提取法，步骤如下：

步骤1：将灰度正弦条纹图投影到被测物体表面，用CCD对被测物体表面进行拍摄，得到一幅宽度为c、高度为l的变形条纹图像g(y，x)：

g(y，x)＝A(y，x)+B(y，x)cos[2πf₀x+φ(y，x)]，

其中，A(y，x)是背景分量，B(y，x)是物体表面反射率，f₀是正弦条纹基频，φ(y，x)是待求的包裹相位分布，(y，x)表示变形条纹图像中各像素点的二维坐标，取值范围分别为：1≤y≤l，1≤x≤c，B(y，x)cos[2πf₀x+φ(y，x)]为基频分量，此处将每个像素点看作一个信号；

步骤2：令y＝1，n＝1，用经验模式分解法即EMD对g(y，x)进行分解，方法如下：

步骤2.1：将g(y，x)记为一行信号g(x)，其中x仍然满足1≤x≤c，该行信号的相位为φ(x)，寻找g(x)的极大值点和极小值点，分别采用公知的三次样条插值法对这些极大值点和极小值点进行插值，然后将这些值连接得到极大值包络线g_max(x)和极小值包络线g_min(x)；

步骤2.2：用原始信号g(x)减去极大值包络线g_max(x)与极小值包络线g_min(x)的平均值，得到h(x)：

$h\left(x\right)=g\left(x\right)-\frac{g_{\max }\left(x\right)+g_{\min }\left(x\right)}{2};$

步骤2.3：分别计算h(x)的平均幅度mean_amp(x)和包络幅度env_amp(x)：

$\mathrm{mean}\_\mathrm{amp}\left(x\right)=\frac{|g_{\max }\left(x\right)+g_{\min }\left(x\right)|}{2},$

$\mathrm{env}\_\mathrm{amp}\left(x\right)=\frac{|g_{\max }\left(x\right)-g_{\min }\left(x\right)|}{2};$

步骤2.4：如果h(x)同时满足以下三个条件，则得到一个本征模函数IMF，记为s_n(x)，且s_n(x)＝h(x)，同时将s_n(x)从原信号g(x)中分离，得到新的信号g′(x)＝g(x)-s_n(x)，否则，直接将h(x)从原信号g(x)中分离，得到新的信号g′(x)＝g(x)-h(x)，所述的三个条件是：

条件1：mean_amp(x)/env_amp(x)＜0.5，

条件2：满足不等式mean_amp(x)/env_amp(x)＜0.5的像素个数占同一整行像素总数的比例小于0.05，

条件3：极大值和极小值的个数之和与h(x)过零点的个数相等或至多相差1个，

步骤2.5：若步骤2.4得到的g′(x)也同时满足以上三个条件，令g(x)＝g′(x)，n＝n+1，返回步骤2.1；否则分解停止，并将g′(x)记为res(x)，得到最后的分解结果为：

$g\left(x\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{s}_n\left(x\right)+\mathrm{res}\left(x\right),$

其中n为一个IMF即s(x)的序数，1≤n≤N，N为IMF的总个数；

步骤3：通过希尔伯特-黄变换，确定能够准确描述一行条纹信号变化规律的瞬时频率，具体过程如下：

步骤3.1：对第n个IMF即s_n(x)做希尔伯特变换，得到：

${s_n}^{\′}\left(x\right)=s_n\left(x\right)*\frac{1}{\mathrm{\πx}}=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{\∞}\frac{s_n\left(\mathrm{\τ}\right)}{x-\mathrm{\τ}}\mathrm{d\τ},$

其中“*”为卷积运算符，τ为积分变量，s_n′(x)为希尔伯特变换的结果，分别构造s_n(x)的解析信号z_n(x)：

其中i为虚部单位，

为该解析信号z_n(x)的模值，

为解析信号z_n(x)的相位；

步骤3.2：对第n个IMF求其瞬时频率f_n(x)：

对第n个IMF求其边际谱M_n(f)：

$M_n\left(f\right)={\∫}_{-\∞}^{\∞}\left|{\mathrm{\λ}}_n\left(x\right)e^{i2\mathrm{\π}{\∫}_{-\∞}^{\∞}{f}_n\left(x\right)\mathrm{dx}}\mathrm{dx}\right|;$

步骤3.3：确定含基频分量最多的IMF的序数K：

$K=E\left\{\right|f_n^{M_{\max }}-f_0\left|\right\};$

其中

为第n个IMF的边际谱最大值M_max所对应的频率值，E{}为取

的最小值所对应的IMF的序数，

如果IMF所对应的

最小值个数为1个以上，则从具有

最小值的IMF中选择出具有最大的边际谱最大值M_max的IMF，并以具有最大的边际谱最大值M_max的IMF所对应的序数为所求的K值，

选取f_K(x)作为能够准确描述该行条纹信号变化规律的瞬时频率；

步骤4：确定一行条纹信号的背景分量，具体过程如下：根据步骤3.2已得的每个IMF的瞬时频率，求出每个IMF的瞬时频率均值，并从中找到最小的瞬时频率均值，然后确定最小的瞬时频率均值所对应的IMF序数，记为K_b，则以步骤2.5所得的第K_b+1个IMF到第N个IMF以及残余分量之和即

为该行条纹信号的背景分量组合；

步骤5：根据一个局部平稳区域内的瞬时频率最大值小于2倍最小值这一属性，用步骤3.3确定的瞬时频率f_K(x)自适应定位一行条纹信号的局部平稳区域，具体过程如下：

步骤5.1：瞬时频率向量f_K(x)为[f_K(1)，f_K(2)，……f_K(c)]，其中元素全部为不小于0的瞬时频率，对f_K(x)求转置向量，得到

为 $\left[\begin{array}{c}{f}_K\left(1\right)\\ f_K\left(2\right)\\ \·\·\·\·\·\·\\ f_K\left(c\right)\end{array}\right];$

步骤5.2：令

中每个元素分别减去整个2×f_K(x)向量，得到一个c×c的方阵F：

$F=\left[\begin{array}{cccc}{f}_K\left(1\right)-2\×f_K\left(1\right)& f_K\left(1\right)-2\×f_K\left(2\right)& \·\·\·\·\·\·& f_K\left(1\right)-2\×f_K\left(c\right)\\ f_K\left(2\right)-2\×f_K\left(1\right)& f_K\left(2\right)-2\×f_K\left(2\right)& \·\·\·\·\·\·& f_K\left(2\right)-2\×f_K\left(c\right)\\ \·\·\·\·\·\·& \·\·\·\·\·\·& \·\·\·\·\·\·& \·\·\·\·\·\·\\ f_K\left(c\right)-2\×f_K\left(1\right)& f_K\left(c\right)-2\×f_K\left(2\right)& \·\·\·\·\·\·& f_K\left(c\right)-2\×f_K\left(c\right)\end{array}\right],$ 令方阵内所有负值为0，并找到方阵中全部元素为0的零对角线，得到：

$F=\left[\begin{array}{cccc}0& f_K\left(1\right)-2\×f_K\left(2\right)& \·\·\·\·\·\·& f_K\left(1\right)-2\×f_K\left(c\right)\\ f_K\left(2\right)-2\×f_K\left(1\right)& 0& \·\·\·\·\·\·& f_K\left(2\right)-2\×f_K\left(c\right)\\ \·\·\·\·\·\·& \·\·\·\·\·\·& \·\·\·\·\·\·& \·\·\·\·\·\·\\ f_K\left(c\right)-2\×f_K\left(1\right)& f_K\left(c\right)-2\×f_K\left(1\right)& \·\·\·\·\·\·& 0\end{array}\right];$

步骤5.3：对方阵F建立坐标系，坐标原点为F左上角第一个元素，其坐标值记为(1，1)，横坐标的坐标方向为方阵的行方向，横坐标取值范围为1～c，纵坐标的坐标方向为方阵的列方向，纵坐标的取值范围仍为1～c，

在方阵F中寻找出以方阵F的零对角线的一部分或全部为对角线的最大全零子方阵，从方阵F的坐标原点开始，沿零对角线方向，依次将这些最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标记录下来，依次用(p₁，p₁)，(p₂，p₂)，……，(p_r，p_r)，(c，c)表示，其中(p₁，p₁)为第一个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标，依次排列，(p_r，p_r)为倒数第二个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标，(c，c)为最后一个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标，

最后，该行条纹信号的平稳区域划分情况为：(1，……，p₁)，(p₁+1，……，p₂)，……，(p_r+1，……，c)；

步骤6：将原条纹信号g(x)减去步骤4中确定的背景分量，得到步骤1中去除A(y，x)后的基频分量g₀(x)，即为g₀(x)＝B(x，y)cos[2πf₀x+φ(x，y)]；

步骤7：判断步骤5确定的局部平稳区域个数是否为1，若“是”，则直接对g₀(x)做快速傅里叶变换得到傅里叶频谱，记作F_f，若“否”，则对g₀(x)进行自适应高斯窗口傅里叶变换，具体过程为：

步骤7.1：对g₀(x)做自适应高斯窗口傅里叶变换：

${\mathrm{WF}}_a^b\left[g_0\left(x\right)\right]={\∫}_{-\∞}^{\∞}\left[g_0\left(x\right)W_a^b\left(x\right)\exp (-\mathrm{j\ωx})\right]\mathrm{dx},$

其中b为平移因子，b依次取值为1、2、3、…、c，a为每一个局部平稳区域内每个像素点对应的高斯窗口的尺度因子，a＝L/4，其中L为相应的局部平稳区域的长度值，单位为像素， $W_a^b\left(x\right)=\frac{1}{a}W\left(\frac{x-b}{a}\right)=\frac{1}{\sqrt{2\mathrm{\π}}a}\exp [-\frac{1}{2}{\left(\frac{x-b}{a}\right)}^2]$ 为高斯窗函数，

对一行条纹信号所有的信号依次进行自适应高斯窗口傅里叶变换后，得到由

构成的一个二维复数方阵，大小为c×c，二维复数方阵每一行的元素为每一个窗口内信号的频谱，二维复数方阵一共有c行，代表b从1取值到c即一共有c个窗口的信号的频谱；

步骤7.2：将步骤7.1所得的二维复数方阵按列叠加，得到总频谱，记作F_f；

步骤8：根据F_f确定基频频谱的范围，并提取出来，记为F_f0；对F_f0求傅里叶逆变换，并根据对F_f0求傅里叶逆变换所得的结果，求得相位角即包裹在[-π，π]之间的相位φ(x)，令y＝y+1，返回步骤2，若y＞l，则进入步骤9；

步骤9：对包裹相位进行展开得到绝对相位，根据经典光栅投影的相位到高度的转换公式，最终求得测量物体的三维信息。

有益效果：

相比相移法，本方法不会受到设备关于正弦性、相移的精度和速度等条件限制，且仅需要一幅变形条纹图即可完成包裹相位的提取，有利于实现动态测量。与其他基于变换域的现有技术相比，本发明具有如下优点：

首先，本发明采用希尔伯特-黄变换法求解条纹信号的瞬时频率，相比于现有常用的小波变换最大脊法或者S变换最大脊法，避免了被测相位必须线性逼近且变化缓慢的条件限制，得到的瞬时频率更加真实的描述了变形条纹信号的变化，通过本方法所得的瞬时频率确定的窗口尺度因子能真实的根据信号的局部特征进行自适应变化，从而达到窗口傅里叶精确提取局部细节相位的目的；

其次，本发明提出的自适应定位局部平稳区域的方法，不需要事先确定任何经验值，也不需要迭代计算，效率高、速度快，大大提高了整个自适应窗口傅里叶相位提取的处理效率；

最后，本发明在求解条纹图的瞬时频率时，可以很方便的“顺便”去除条纹图的背景分量，大大减少了提取基频分量时背景分量所带来的零频频谱干扰，当确定的窗口尺度因子较小时，这一干扰将会非常严重，因此本发明这一处理会有很大的优势。

综上，本发明仅用一幅变形条纹图，便能够准确地获取被测物体的包裹相位，且克服了传统方法对具有复杂面形或陡峭边缘的物体测量不高的缺点。

附图说明

图1是本发明整个过程的流程图。

图2是CCD采集到的以塑料泡沫为被测物体的变形条纹图像。

图3是步骤2中对一行条纹信号进行经验模式分解的具体过程流程图。

图4是图2中直线所代表的一行条纹信号及其EMD分解结果。

图5是图4中相应IMF的瞬时频率。

图6是图4中相应IMF的边际谱。

图7是步骤5的示意图。

图8是对图2中直线所在行条纹信号所定位的局部平稳区域。

图9是去掉背景分量的变形条纹图。

图10是所求的尺度因子分布图。

图11是所求的包裹相位图。

图12是最终恢复的相位图。

具体实施方式

下面结合附图对本发明的具体实施方式做进一步说明。在Windows操作系统下选用VC++6.0作为编程工具，对CCD采集到的变形条纹图像进行处理。该实例采用塑料泡沫作为被测物体，最终得到比较精确的含有三维信息的全场包裹相位分布和。

图1是本发明的整个过程的流程图。

光学三维测量中条纹图进行包裹相位提取时，针对现有的时频分析方法对表面复杂或含有陡峭边缘的物体的测量精度较低等问题，本发明提出基于希尔伯特-黄变换的多尺度窗口傅里叶相位提取法。首先采用经验模态分解法逐行分解采集到的变形条纹图像并进行希尔伯特-黄变换处理，逐行分析确定能够描述变形条纹变化情况的瞬时频率，同时确定条纹图的背景分量。将原条纹图去除背景分量。根据确定的瞬时频率逐行定位每一行信号的局部平稳区域，若任意一行信号的局部平稳区域个数为1，则对该行信号直接进行傅里叶变换然后提取基频分量的频谱。若任意一行信号的局部平稳区域个数不为1，根据所确定的局部平稳区域的长度，计算每一个局部平稳区域内每个信号相应的高斯窗函数的尺度因子。由不同的尺度因子引导不同尺寸的高斯窗，对相应的信号逐点进行由高斯窗函数引导的窗口傅里叶变换，每个像素点的变换结果都是一个一维复数数组即傅里叶频谱，将基频频谱提取出来。对一行信号逐个信号进行相同的处理，所有的处理结果依次按行排列，得到一个二维复数矩阵。根据高斯窗函数的时频统计特性，将所有的提取出的基频频谱按列相加，得到最终的一行信号的基频频谱。对相加后的一行信后的基频频谱进行傅里叶反变换，得到一行信号的基频分量的相位。对条纹图逐行进行处理后，最终可得整幅条纹图的全场包裹相位分布。由于所采用的方法不会使所求的瞬时频率被局部平滑化，故根据其所确定的窗口尺寸能够较好的反映条纹信号的局部信息，从而能够精确地测量条纹图细节处的相位。同时，由于背景分量被去除，故在基频提取的过程中极大地减少了零频频谱的干扰，有利于基频频谱的精确提取。

本发明的具体实施步骤如下：

步骤1：将灰度正弦条纹图投影到被测物体表面，用CCD对被测物体表面进行拍摄，得到一幅宽度为c、高度为l的变形条纹图像g(y，x)：

g(y，x)＝A(y，x)+B(y，x)cos[2πf₀x+φ(y，x)]，

其中，A(y，x)是背景分量，B(y，x)是物体表面反射率，f₀是正弦条纹基频，φ(y，x)是待求的包裹相位分布，(y，x)表示变形条纹图像中各像素点的二维坐标，取值范围分别为：1≤y≤l，1≤x≤c，B(y，x)cos[2πf₀x+φ(y，x)]为基频分量，此处将每个像素点看作一个信号。

图2为所得的变形条纹图像，大小为1020×1020，图中直线为y＝170时的一行示例条纹信号g(x)。

步骤2：令y＝1，n＝1，用经验模式分解法即EMD对g(y，x)进行分解，方法如下：

步骤2.1：将g(y，x)记为一行信号g(x)，其中x仍然满足1≤x≤c，该行信号的相位为φ(x)，寻找g(x)的极大值点和极小值点，分别采用公知的三次样条插值法对这些极大值点和极小值点进行插值，然后将这些值连接得到极大值包络线g_max(x)和极小值包络线g_min(x)。此处也可选择其他插值法，具体的选取办法可参考文献“N.E.Huang，Z.Shen，and S.R.Long et al，“The empirical mode decompositionand the Hilbert spectrum for nonliner and non-stationary time series analysis，”J.Proc.Lond.A.Great Britain，454，903-995(1998)”，及“G.Rilling，P.Flandrin and P.Goncalves，“Onempirical mode decomposition and its algorithm，”in IEEEEURASIP Workshop onNonlinear Signal andImage Processing，NSIP-03，GRADO(I)(IEEE，2003)”；

步骤2.2：用原始信号g(x)减去极大值包络线g_max(x)与极小值包络线g_min(x)的平均值，得到h(x)：

$h\left(x\right)=g\left(x\right)-\frac{g_{\max }\left(x\right)+g_{\min }\left(x\right)}{2};$

步骤2.3：分别计算h(x)的平均幅度mean_amp(x)和包络幅度env_amp(x)：

$\mathrm{mean}\_\mathrm{amp}\left(x\right)=\frac{|g_{\max }\left(x\right)+g_{\min }\left(x\right)|}{2},$

$\mathrm{env}\_\mathrm{amp}\left(x\right)=\frac{|g_{\max }\left(x\right)-g_{\min }\left(x\right)|}{2};$

步骤2.4：如果h(x)同时满足以下三个条件，则得到一个本征模函数IMF，记为s_n(x)，且s_n(x)＝h(x)，同时将s_n(x)从原信号g(x)中分离，得到新的信号g′(x)＝g(x)-s_n(x)，否则，直接将h(x)从原信号g(x)中分离，得到新的信号g′(x)＝g(x)-h(x)，所述的三个条件是：

条件1：mean_amp(x)/env_amp(x)＜0.5；

条件2：满足不等式mean_amp(x)/env_amp(x)＜0.5的像素个数占同一整行像素总数的比例小于0.05；

条件3：极大值和极小值的个数之和与h(x)过零点的个数相等或至多相差1个；

步骤2.5：若步骤2.4得到的g′(x)也同时满足以上三个条件，令g(x)＝g′(x)，n＝n+1，返回步骤2.1；否则分解停止，并将g′(x)记为res(x)，得到最后的分解结果为：

$g\left(x\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{s}_n\left(x\right)+\mathrm{res}\left(x\right),$

其中n为一个IMF即s(x)的序数，1≤n≤N，N为IMF的总个数。

图3为一行信号分解过程的流程图。这里以变形条纹图的一行变形信号为例，即对如图2中直线所代表信号g(x)进行处理。在处理过程中，三个条件的阈值0.5、0.05和0.05为建议阈值，实际应用过程中可根据需要适当微调整，调整的前提是使分解彻底，即剩余分量不能再分离出IMF，同时又不能过分解。图4为原信号及其被分解后的结果，从中可看见，IMFs随着尺度的变化由高到底依次排列，第一个IMF为高频噪声分量，res(x)的尺度最大，趋近于0，它能够描述一行信号的整体趋势。

步骤3：通过希尔伯特-黄变换，确定能够准确描述一行条纹信号变化规律的瞬时频率，具体过程如下：

步骤3.1：对第n个IMF即s_n(x)做希尔伯特变换，得到：

${s_n}^{\′}\left(x\right)=s_n\left(x\right)*\frac{1}{\mathrm{\πx}}=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{\∞}\frac{s_n\left(\mathrm{\τ}\right)}{x-\mathrm{\τ}}\mathrm{d\τ},$

其中“*”为卷积运算符，τ为积分变量，s_n′(x)为希尔伯特变换的结果，分别构造s_n(x)的解析信号z_n(x)：

其中i为虚部单位，

为该解析信号z_n(x)的模值，

为解析信号z_n(x)的相位；

步骤3.2：对第n个IMF求其瞬时频率f_n(x)：

对第n个IMF求其边际谱M_n(f)：

$M_n\left(f\right)={\∫}_{-\∞}^{\∞}\left|{\mathrm{\λ}}_n\left(x\right)e^{i2\mathrm{\π}{\∫}_{-\∞}^{\∞}{f}_n\left(x\right)\mathrm{dx}}\mathrm{dx}\right|,$ (1≤n≤N)；

步骤3.3：确定含基频分量最多的IMF的序数K：

$K=E\left\{\right|f_n^{M_{\max }}-f_0\left|\right\},$ (1≤n≤N)，

其中为第n个IMF的边际谱最大值M_max所对应的频率值，E{}为取

的最小值所对应的IMF的序数，

如果IMF所对应的

最小值个数为1个以上，则从具有

最小值的IMF中选择出具有最大的边际谱最大值M_max的IMF，并以具有最大的边际谱最大值M_max的IMF所对应的序数为所求的K值。

选取f_K(x)作为能够准确描述该行条纹信号变化规律的瞬时频率。

如图5为图4中每个IMF所得的瞬时频率，图6为相应的边际谱。图6中，由于条纹的基频是0.05，可见IMF2的边际谱最大值所对应的频率值与基频0.05的差值最小，故IMF2为含有基频分量最多的分量，IMF-f₂则被选定为能够准确描述该行条纹信号变化规律的瞬时频率。

步骤4：确定一行条纹信号的背景分量，具体过程如下：根据步骤3.2已得的每个IMF的瞬时频率，求出每个IMF的瞬时频率均值，并从中找到最小的瞬时频率均值，然后确定最小的瞬时频率均值所对应的IMF序数，记为K_b，则以步骤2.5所得的第K_b+1个IMF到第N个IMF以及残余分量之和即

为该行条纹信号的背景分量组合。

步骤5：根据一个局部平稳区域内的瞬时频率最大值小于2倍最小值这一属性，用步骤3.3确定的瞬时频率f_K(x)自适应定位一行条纹信号的局部平稳区域，具体过程如下：

步骤5.1：瞬时频率向量f_K(x)为[f_K(1)，f_K(2)，……f_K(c)]，其中元素全部为不小于0的瞬时频率，对f_K(x)求转置向量，得到

为 $\left[\begin{array}{c}{f}_K\left(1\right)\\ f_K\left(2\right)\\ \·\·\·\·\·\·\\ f_K\left(c\right)\end{array}\right];$

步骤5.2：令

中每个元素分别减去整个2×f_K(x)向量，得到一个c×c的方阵F：

$F=\left[\begin{array}{cccc}{f}_K\left(1\right)-2\×f_K\left(1\right)& f_K\left(1\right)-2\×f_K\left(2\right)& \·\·\·\·\·\·& f_K\left(1\right)-2\×f_K\left(c\right)\\ f_K\left(2\right)-2\×f_K\left(1\right)& f_K\left(2\right)-2\×f_K\left(2\right)& \·\·\·\·\·\·& f_K\left(2\right)-2\×f_K\left(c\right)\\ \·\·\·\·\·\·& \·\·\·\·\·\·& \·\·\·\·\·\·& \·\·\·\·\·\·\\ f_K\left(c\right)-2\×f_K\left(1\right)& f_K\left(c\right)-2\×f_K\left(2\right)& \·\·\·\·\·\·& f_K\left(c\right)-2\×f_K\left(c\right)\end{array}\right],$

令方阵内所有负值为0，并找到方阵中全部元素为0的零对角线，得到：

$F=\left[\begin{array}{cccc}0& f_K\left(1\right)-2\×f_K\left(2\right)& \·\·\·\·\·\·& f_K\left(1\right)-2\×f_K\left(c\right)\\ f_K\left(2\right)-2\×f_K\left(1\right)& 0& \·\·\·\·\·\·& f_K\left(2\right)-2\×f_K\left(c\right)\\ \·\·\·\·\·\·& \·\·\·\·\·\·& \·\·\·\·\·\·& \·\·\·\·\·\·\\ f_K\left(c\right)-2\×f_K\left(1\right)& f_K\left(c\right)-2\×f_K\left(1\right)& \·\·\·\·\·\·& 0\end{array}\right];$

步骤5.3：对方阵F建立坐标系，坐标原点为F左上角第一个元素，其坐标值记为(1，1)，横坐标的坐标方向为方阵的行方向，横坐标取值范围为1～c，纵坐标的坐标方向为方阵的列方向，纵坐标的取值范围仍为1～c，

在方阵F中寻找出以方阵F的零对角线的一部分或全部为对角线的最大全零子方阵，从方阵F的坐标原点开始，沿零对角线方向，依次将这些最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标记录下来，依次用(p₁，p₁)，(p₂，p₂)，……，(p_r，p_r)，(c，c)表示，其中(p₁，p₁)为第一个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标，依次排列，(p_r，p_r)为倒数第二个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标，(c，c)为最后一个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标，

最后，该行条纹信号的平稳区域划分情况为：(1，……，p₁)，(p₁+1，……，p₂)，……，(p_r+1，……，c)。

图7为自适应定位一行条纹信号局部平稳区域的直观描述，图8为图3中对直线代表信号局部平稳区域的划分。

步骤6：将原条纹信号g(x)减去步骤4中确定的背景分量，得到步骤1中去除A(y，x)后的基频分量g₀(x)，即为g₀(x)＝B(x，y)cos[2πf₀x+φ(x，y)]。

图9为去除背景分量后的条纹图。

步骤7：判断步骤5确定的局部平稳区域个数是否为1，若“是”，则直接对g₀(x)做快速傅里叶变换得到傅里叶频谱，记作F_f，若“否”，则对g₀(x)进行自适应高斯窗口傅里叶变换，具体过程为：

步骤7.1：对g₀(x)做自适应高斯窗口傅里叶变换：

${\mathrm{WF}}_a^b\left[g_0\left(x\right)\right]={\∫}_{-\∞}^{\∞}\left[g_0\left(x\right)W_a^b\left(x\right)\exp (-\mathrm{j\ωx})\right]\mathrm{dx},$

其中b为平移因子，b依次取值为1、2、3、…、c，a为每一个局部平稳区域内每个像素点对应的高斯窗口的尺度因子，a＝L/4，其中L为相应的局部平稳区域的长度值，单位为像素， $W_a^b\left(x\right)=\frac{1}{a}W\left(\frac{x-b}{a}\right)=\frac{1}{\sqrt{2\mathrm{\π}}a}\exp [-\frac{1}{2}{\left(\frac{x-b}{a}\right)}^2]$ 为高斯窗函数，

对一行条纹信号所有的信号依次进行自适应高斯窗口傅里叶变换后，得到由

构成的一个二维复数方阵，大小为c×c，二维复数方阵每一行的元素为每一个窗口内信号的频谱，二维复数方阵一共有c行，代表b从1取值到c即一共有c个窗口的信号的频谱；

步骤7.2：将步骤7.1所得的二维复数方阵按列叠加，得到总频谱，记作F_f。

步骤8：根据F_f确定基频频谱的范围，并提取出来，记为F_f0；对F_f0求傅里叶逆变换，并根据对F_f0求傅里叶逆变换所得的结果，求得相位角即包裹在[-π，π]之间的相位φ(x)，令y＝y+1，返回步骤2，若y＞l，则进入步骤9。

步骤9：对包裹相位进行展开得到绝对相位，此处相位展开的方法采用质量图相位展开法。最后，根据经典光栅投影的相位到高度的转换公式求得测量物体的三维信息。高度转换公式如下：

其中，l、d是测量系统的几何参数，l是投影仪到测量平面的距离，d是CCD摄像头到投影仪的距离，

表示相位变化量，

为展开相位结果，

为初始相位结果，由测量参考面决定，ω₀为投影光栅的角频率，由系统标定可得。

如图10为尺度因子的全场分布图，可见用本方法能够较好地描述物体的变化规律。图11为对一整幅变形条纹图进行本发明方法处理后所得的全场包裹相位分布图。图12为所测物体的相位变化分布图，即先将全场包裹相位用质量图法进行相位展开，然后用同样方法得到没有被物体调制的参考平面条纹图的绝对相位图，将包含物体三维信息的变形条纹图的展开相位与参考平面条纹图的展开相位相减即可得到相位变化分布图。用本方法所测的包裹相位，在洪水法进行相位展开时，由于物体边缘相位测量的精度较高，基本没有误差传递，而其他方法会有大片的误差传递现象。故本方法对复杂面形或含有陡峭边缘物体的测量时精度会有具有较大的提高。

Claims (1)

1.一种光学三维测量中的自适应窗口傅里叶相位提取法，具体步骤如下：

步骤1：将灰度正弦条纹图投影到被测物体表面，用CCD对被测物体表面进行拍摄，得到一幅宽度为c、高度为l的变形条纹图像g(y，x)：

g(y，x)＝A(y，x)+B(y，x)cos[2πf₀x+φ(y，x)]，

其中，A(y，x)是背景分量，B(y，x)是物体表面反射率，f₀是正弦条纹基频，φ(y，x)是待求的包裹相位分布，(y，x)表示变形条纹图像中各像素点的二维坐标，取值范围分别为：1≤y≤l，1≤x≤c，B(y，x)cos[2πf₀x+φ(y，x)]为基频分量，此处将每个像素点看作一个信号；

步骤2：令y＝1，n＝1，用经验模式分解法即EMD对g(y，x)进行分解，方法如下：

步骤2.1：将g(y，x)记为一行信号g(x)，其中x仍然满足1≤x≤c，该行信号的相位为φ(x)，寻找g(x)的极大值点和极小值点，分别采用公知的三次样条插值法对这些极大值点和极小值点进行插值，然后将这些值连接得到极大值包络线g_max(x)和极小值包络线g_min(x)；

步骤2.2：用原始信号g(x)减去极大值包络线g_max(x)与极小值包络线g_min(x)的平均值，得到h(x)：

$h\left(x\right)=g\left(x\right)-\frac{g_{\max }\left(x\right)+g_{\min }\left(x\right)}{2};$

步骤2.3：分别计算h(x)的平均幅度mean_amp(x)和包络幅度env_amp(x)：

$\mathrm{mean}\_\mathrm{amp}\left(x\right)=\frac{|g_{\max }\left(x\right)+g_{\min }\left(x\right)|}{2},$

$\mathrm{env}\_\mathrm{amp}\left(x\right)=\frac{|g_{\max }\left(x\right)-g_{\min }\left(x\right)|}{2};$

步骤2.4：如果h(x)同时满足以下三个条件，则得到一个本征模函数IMF，记为s_n(x)，且s_n(x)＝h(x)，同时将s_n(x)从原信号g(x)中分离，得到新的信号g′(x)＝g(x)-s_n(x)，否则，直接将h(x)从原信号g(x)中分离，得到新的信号g′(x)＝g(x)-h(x)，所述的三个条件是：

条件1：mean_amp(x)/env_amp(x)＜0.5，

条件2：满足不等式mean_amp(x)/env_amp(x)＜0.5的像素个数占同一整行像素总数的比例小于0.05，

条件3：极大值和极小值的个数之和与h(x)过零点的个数相等或至多相差1个，

步骤2.5：若步骤2.4得到的g′(x)也同时满足以上三个条件，令g(x)＝g′(x)，n＝n+1，返回步骤2.1；否则分解停止，并将g′(x)记为res(x)，得到最后的分解结果为：

$g\left(x\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{s}_n\left(x\right)+\mathrm{res}\left(x\right),$

其中n为一个IMF即s(x)的序数，1≤n≤N，N为IMF的总个数；

步骤3：通过希尔伯特-黄变换，确定能够准确描述一行条纹信号变化规律的瞬时频率，具体过程如下：

步骤3.1：对第n个IMF即s_n(x)做希尔伯特变换，得到：

${s_n}^{\′}\left(x\right)=s_n\left(x\right)*\frac{1}{\mathrm{\πx}}=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{\∞}\frac{s_n\left(\mathrm{\τ}\right)}{x-\mathrm{\τ}}\mathrm{d\τ},$

其中“*”为卷积运算符，τ为积分变量，s_n′(x)为希尔伯特变换的结果，分别构造s_n(x)的解析信号z_n(x)：

其中i为虚部单位，

为该解析信号z_n(x)的模值，

为解析信号z_n(x)的相位；

步骤3.2：对第n个IMF求其瞬时频率f_n(x)：

对第n个IMF求其边际谱M_n(f)：

$M_n\left(f\right)={\∫}_{-\∞}^{\∞}\left|{\mathrm{\λ}}_n\left(x\right)e^{i2\mathrm{\π}{\∫}_{-\∞}^{\∞}{f}_n\left(x\right)\mathrm{dx}}\mathrm{dx}\right|;$

步骤3.3：确定含基频分量最多的IMF的序数K：

$K=E\left\{\right|f_n^{M_{\max }}-f_0\left|\right\};$

其中

为第n个IMF的边际谱最大值M_max所对应的频率值，E{}为取

的最小值所对应的IMF的序数，

如果IMF所对应的

最小值个数为1个以上，则从具有

最小值的IMF中选择出具有最大的边际谱最大值M_max的IMF，并以具有最大的边际谱最大值M_max的IMF所对应的序数为所求的K值，

选取f_K(x)作为能够准确描述该行条纹信号变化规律的瞬时频率；

步骤4：确定一行条纹信号的背景分量，具体过程如下：根据步骤3.2已得的每个IMF的瞬时频率，求出每个IMF的瞬时频率均值，并从中找到最小的瞬时频率均值，然后确定最小的瞬时频率均值所对应的IMF序数，记为K_b，则以步骤2.5所得的第K_b+1个IMF到第N个IMF以及残余分量之和即

为该行条纹信号的背景分量组合；

步骤5：根据一个局部平稳区域内的瞬时频率最大值小于2倍最小值这一属性，用步骤3.3确定的瞬时频率f_K(x)自适应定位一行条纹信号的局部平稳区域，具体过程如下：

步骤5.1：瞬时频率向量f_K(x)为[f_K(1)，f_K(2)，……f_K(c)]，其中元素全部为不小于0的瞬时频率，对f_K(x)求转置向量，得到

为 $\left[\begin{array}{c}{f}_K\left(1\right)\\ f_K\left(2\right)\\ \·\·\·\·\·\·\\ f_K\left(c\right)\end{array}\right];$

步骤5.2：令中每个元素分别减去整个2×f_K(x)向量，得到一个c×c的方阵F：

$F=\left[\begin{array}{cccc}{f}_K\left(1\right)-2\×f_K\left(1\right)& f_K\left(1\right)-2\×f_K\left(2\right)& \·\·\·\·\·\·& f_K\left(1\right)-2\×f_K\left(c\right)\\ f_K\left(2\right)-2\×f_K\left(1\right)& f_K\left(2\right)-2\×f_K\left(2\right)& \·\·\·\·\·\·& f_K\left(2\right)-2\×f_K\left(c\right)\\ \·\·\·\·\·\·& \·\·\·\·\·\·& \·\·\·\·\·\·& \·\·\·\·\·\·\\ f_K\left(c\right)-2\×f_K\left(1\right)& f_K\left(c\right)-2\×f_K\left(2\right)& \·\·\·\·\·\·& f_K\left(c\right)-2\×f_K\left(c\right)\end{array}\right],$

令方阵内所有负值为0，并找到方阵中全部元素为0的零对角线，得到：

$F=\left[\begin{array}{cccc}0& f_K\left(1\right)-2\×f_K\left(2\right)& \·\·\·\·\·\·& f_K\left(1\right)-2\×f_K\left(c\right)\\ f_K\left(2\right)-2\×f_K\left(1\right)& 0& \·\·\·\·\·\·& f_K\left(2\right)-2\×f_K\left(c\right)\\ \·\·\·\·\·\·& \·\·\·\·\·\·& \·\·\·\·\·\·& \·\·\·\·\·\·\\ f_K\left(c\right)-2\×f_K\left(1\right)& f_K\left(c\right)-2\×f_K\left(1\right)& \·\·\·\·\·\·& 0\end{array}\right];$

步骤5.3：对方阵F建立坐标系，坐标原点为F左上角第一个元素，其坐标值记为(1，1)，横坐标的坐标方向为方阵的行方向，横坐标取值范围为1～c，纵坐标的坐标方向为方阵的列方向，纵坐标的取值范围仍为1～c，

在方阵F中寻找出以方阵F的零对角线的一部分或全部为对角线的最大全零子方阵，从方阵F的坐标原点开始，沿零对角线方向，依次将这些最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标记录下来，依次用(p₁，p₁)，(p₂，p₂)，……，(p_r，p_r)，(c，c)表示，其中(p₁，p₁)为第一个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标，依次排列，(p_r，p_r)为倒数第二个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标，(c，c)为最后一个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标，

最后，该行条纹信号的平稳区域划分情况为：(1，……，p₁)，(p₁+1，……，p₂)，……，(p_r+1，……，c)；

步骤6：将原条纹信号g(x)减去步骤4中确定的背景分量，得到步骤1中去除A(y，x)后的基频分量g₀(x)，即为g₀(x)＝B(x，y)cos[2πf₀x+φ(x，y)]；

步骤7：判断步骤5确定的局部平稳区域个数是否为1，若“是”，则直接对g₀(x)做快速傅里叶变换得到傅里叶频谱，记作F_f，若“否”，则对g₀(x)进行自适应高斯窗口傅里叶变换，具体过程为：

步骤7.1：对g₀(x)做自适应高斯窗口傅里叶变换：

${\mathrm{WF}}_a^b\left[g_0\left(x\right)\right]={\∫}_{-\∞}^{\∞}\left[g_0\left(x\right)W_a^b\left(x\right)\exp (-\mathrm{j\ωx})\right]\mathrm{dx},$

其中b为平移因子，b依次取值为1、2、3、…、c，a为每一个局部平稳区域内每个像素点对应的高斯窗口的尺度因子，a＝L/4，其中L为相应的局部平稳区域的长度值，单位为像素， $W_a^b\left(x\right)=\frac{1}{a}W\left(\frac{x-b}{a}\right)=\frac{1}{\sqrt{2\mathrm{\π}}a}\exp [-\frac{1}{2}{\left(\frac{x-b}{a}\right)}^2]$ 为高斯窗函数，

对一行条纹信号所有的信号依次进行自适应高斯窗口傅里叶变换后，得到由

构成的一个二维复数方阵，大小为c×c，二维复数方阵每一行的元素为每一个窗口内信号的频谱，二维复数方阵一共有c行，代表b从1取值到c即一共有c个窗口的信号的频谱；

步骤7.2：将步骤7.1所得的二维复数方阵按列叠加，得到总频谱，记作F_f；

步骤8：根据F_f确定基频频谱的范围，并提取出来，记为F_f0；对F_f0求傅里叶逆变换，并根据对F_f0求傅里叶逆变换所得的结果，求得相位角即包裹在[-π，π]之间的相位φ(x)，令y＝y+1，返回步骤2，若y＞l，则进入步骤9；

步骤9：对包裹相位进行展开得到绝对相位，根据经典光栅投影的相位到高度的转换公式，最终求得测量物体的三维信息。

Images