利用光流场进行物体变形相位测量的方法

摘要

本发明公开了利用光流场进行物体变形相位测量的方法，包括以下步骤：步骤一：根据光流场理论的基本假设，获得两帧连续图像之间的运动矢量的两个分量u和v；步骤二：由窗口傅里叶变换法提取第一帧图像条纹的横向频率f

说明书

技术领域

本发明涉及一种利用光流场进行物体变形相位测量的方法。

背景技术

光学干涉测量方法因其独特的优势在物体的三维形貌和变形测量中占有十分重要的地位，如电子散斑干涉、全息干涉、云纹干涉法等。在这些方法中需要解调出条纹图中每一点的相位信息才可以实现变形测量。解调相位的方法主要分为时间相位测量法和空间相位测量法两大类。时间相位测量法，如时间相移法(TemporalPhase-shiftingMethod,TPM)等等，具有位移测量精度高，结构简单等优点，但是该类方法需要的条纹图较多，耗时较长，一般情况下只适用于静态测量；空间相位测量法是通过在空间上获取更多信息来提取相位的，如空间相移法(SpatialPhase-shiftingMethod,SPM)和Fourier变换法(FourierTransformMethod,FTM)等等，可实现对动态过程相位的测量。但是，前者的光路较为复杂，且精度较低；后者的处理过程比较繁琐，很难实现条纹处理的自动化。同时，以上方法对于位移场的测量均没有考虑到时间参量，因而对动态场的定量测量就无法获得满意的结果。针对此问题，有研究工作者提出了时间序列相位法(TimeSequencePhaseMethod,TSPM)，将时间参量引入散斑干涉当中，成功地提取出了相位信息，能适应大变形、时变测量。然而该方法在位移方向性上无法给出确定值，故只能用于单调的时变场中。除此之外，以上方法都需要进行相位解包络操作，这无疑增加了相位提取的运算量和误差。同时，由于在条纹稠密区域无法判断每个条纹的细节，故以上方法在条纹密集处解调的相位误差都比较大。

20世纪50年代，Gribson等研究人员基于运动结构重建原理(StructureFromMotion,SFM)第一次提出了光流(opticalflow)的概念：在随时间变化的二维图像序列中，运动的三维场景由于亮度模式的变化而产生的流动即为光流。从此以后，国内外学者陆续提出了不同的光流方法，在精度，鲁棒性和实时性等多方面均取得了重大突破。1997年，国内研究人员提出可以将光流场分析法引入光学干涉计量当中，为光干涉动态测量提供了新的思路，但是其并未做深入的探讨。2011年，J.Vargas等人通过运用光流场理论和SPT算符成功提取出了单幅图像相位，但是该方法仍需要相位解包等一系列复杂操作，且不适用于实现动态形变测量。

发明内容

本发明的目的就是为了解决上述问题，提供一种利用光流场进行物体变形相位测量的方法，只需在空间域中用两幅条纹图就可以完成变形相位的提取，无需转换到频域和相位解包络操作，适用于动态测量，且在条纹越密集的区域提取相位的效果越好，完全克服了现有相位解调方法中条纹过密时解调误差较大的缺点。

为了实现上述目的，本发明采用如下技术方案：

利用光流场进行物体变形相位测量的方法，包括以下步骤：

步骤一：根据光流场理论的基本假设，获得两帧连续图像之间的运动矢量的两个分量u和v；

步骤二：由窗口傅里叶变换法提取第一帧图像条纹的横向频率f_x0和纵向频率f_y0；

步骤三：利用所述步骤一得到的运动矢量的两个分量u和v及所述步骤二得到条纹的横向频率f_x0和纵向频率f_y0计算得到相位的改变量。

所述步骤一的具体方法为：运用公式：

uⁿ⁺¹＝uⁿ-I_x[(I_xuⁿ+I_yvⁿ+I_t)]/(α²+I_x²+I_y²)和

vⁿ⁺¹＝vⁿ-I_y[(I_xuⁿ+I_yvⁿ+I_t)]/(α²+I_x²+I_y²)迭代得到，

其中，I_x、I_y和I_t分别为图像I对于x，y，t的偏导数，I为t时刻像素点(x，y)处灰度值，uⁿ⁺¹和vⁿ⁺¹为在图像的每点第n+1次迭代后所得到的运动矢量的两个分量u和v的值，即两幅图像之间的运动矢量场；uⁿ和vⁿ为在图像的每点第n次迭代后所得到的运动矢量的两个分量u和v的值；α为平滑参数。

所述步骤二的结果由公式求得，(f_x0,f_y0)为使|SI(x，y，ξ，η)|取最大值的(ξ,η)，其中，$SI(\mu ,v;\xi ,\eta )=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}I(x,y)g(x-\mu ,y-v)\exp (-j\xi x-j\eta y)dxdy,$(ξ,η)为空间一点(x,y)对应的频域坐标，SI(μ,ν；ξ,η)是I(x,y)经过窗口傅里叶变换后对应频域的结果，I(x,y)为一幅输入的二维图像；g(x,y)为窗口函数，取高斯函数σ_x和σ_y为高斯分布的标准差，其取值决定着高斯窗口的大小。

所述步骤三的具体方法为：将运动矢量的两个分量u和v，条纹的横向频率f_x0和纵向频率f_y0带入公式为点(x₀，y₀)在Δt时间内相位的改变量，求得两帧连续图像的相位该变量。

本发明的有益效果：

本发明只需在空间域中用两幅条纹图就可以完成变形相位的提取，无需转换到频域和相位解包络操作，适用于动态测量，且在条纹越密集的区域提取相位的效果越好，完全克服了现有相位解调方法中条纹过密时解调误差较大的缺点，过程简单、方便，同时该方法也为光干涉动态测量提供了新的途径。

附图说明

图1为原始干涉条纹模拟图；

图2(a)为附加相位为后的条纹模拟图,图2(b)为附加相位后的条纹模拟图；

图3为附加相位为时所估算的在y＝127处x方向光流场横截面图；

图4为对附加相位提取的归一化模拟结果；

图5(a)为对附加相位提取的归一化理论结果，图5(b)为对附加相位提取的归一化模拟结果；

图6(a)为实验获取的初始条纹图，图6(b)为实验获取的附加π/7相位后的条纹图；

图7为实验结果图。

具体实施方式

下面结合附图与实施例对本发明作进一步说明。

利用光流场进行物体变形相位测量的方法，包括以下步骤：

步骤一：根据光流场理论的基本假设，获得两帧连续图像之间的运动矢量的两个分量u和v；

步骤二：由窗口傅里叶变换法提取第一帧图像条纹的横向频率f_x0和纵向频率f_y0；

步骤三：利用所述步骤一得到的运动矢量的两个分量u和v及所述步骤二得到条纹的横向频率f_x0和纵向频率f_y0计算得到相位的改变量。

1运动矢量的计算

设在t时刻像素点(x,y)处灰度值为I(x,y,t)，在(t+Δt)时刻该点运动到新的位置(x+Δx,y+Δy)，灰度值记为I(x+Δx,y+Δy)。根据图像一致性假设，即图像沿着运动轨迹亮度保持不变，满足$\frac{dI(x,y,t)}{dt}=0,$则

I(x，y，t)＝I(x+Δx，y+Δy，t+Δt)(1)

设u和v分别为该点光流矢量沿x和y方向的两个分量，且将I(x+Δx,y+Δy,t+Δt)用泰勒公式展开，得到：

$I(x+\Delta x,y+\Delta y,t+\Delta t)=I(x,y,t)+\frac{\partial I}{\partial x}\Delta x+\frac{\partial I}{\partial y}\Delta y+\frac{\partial I}{\partial t}\Delta t---\left(2\right)$

式(2)忽略二阶以上的高次项，有

$\frac{\partial I}{\partial x}\Delta x+\frac{\partial I}{\partial y}\Delta y+\frac{\partial I}{\partial t}\Delta t=0---\left(3\right)$

由于Δt→0，式(3)变为

$\frac{\partial I}{\partial x}\frac{\Delta x}{\Delta t}+\frac{\partial I}{\partial y}\frac{\Delta y}{\Delta t}+\frac{\partial I}{\partial t}=0---\left(4\right)$

I_xu+I_yv+I_t＝0(5)

其中I_x、I_y和I_t分别为图像I对于x，y，t的偏导数，它们的值可以用图像序列中相邻图像目标像素点的一阶差分估计得到。公式(5)即为光流基本等式，由于只有一个方程，只能求得光流沿梯度方向上的值I_x、I_y和I_t，而u,v的解是非唯一的，故需要加入约束条件才能求解速度矢量。根据Horn-Schunck算法计算光流运动场。

Horn-Schunck算法的基本思想是在求解光流时要求光流本身尽可能的平滑。所谓平滑，就是在给定的邻域内尽可能的小，即

$\min \{{\left[\frac{\partial u}{\partial x}\right]}^2+{\left[\frac{\partial u}{\partial y}\right]}^2+{\left[\frac{\partial v}{\partial x}\right]}^2+{\left[\frac{\partial v}{\partial y}\right]}^2\}---\left(6\right)$

综合式(5)和式(6)，Horn-Schunck算法将光流(u,v)的计算归结为式(7)的变分问题。

$\min \{\int \int {\left(I_{x}u+I_{y}v+I_t\right)}^2+{\alpha }^2[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial u}{\partial y}\right)}^2+{\left(\frac{\partial v}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2\left]\right\}---\left(7\right)$

可以得到相应的欧拉-拉格朗日方程，并利用高斯-赛德尔方法进行求解，得到图像上每点第(n+1)次迭代估计的(uⁿ⁺¹,vⁿ⁺¹)为：

uⁿ⁺¹＝uⁿ-I_x[(I_xuⁿ+I_yvⁿ+I_t)]/(α²+I_x²+I_y²)(8)

vⁿ⁺¹＝vⁿ-I_y[(I_xuⁿ+I_yvⁿ+I_t)]/(α²+I_x²+I_y²)(9)

根据公式(8)和(9)迭代所得到的u和v即为两帧图像之间的运动矢量场的两个分量。

2条纹频率的计算

为提取条纹频率，引入以下窗口傅里叶运算：

$SI(\mu ,v;\xi ,\eta )=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}I(x,y)g(x-\mu ,y-v)\exp (-j\xi x-j\eta y)dxdy---\left(10\right)$

式(10)中，(ξ,η)为空间一点(x,y)对应的频域坐标，SI(μ,ν；ξ,η)是I(x,y)经过窗口傅里叶变换后对应频域的结果。g(x,y)为窗口函数，取高斯函数σ_x和σ_y为高斯分布的标准差，其取值决定着高斯窗口的大小。

$(f_{x0},f_{y0})=\arg \underset{\xi ,\eta }{\max }\left(\right|SI\left(x,y,\xi ,\eta \right)\left|\right)---\left(11\right)$

其中(f_x0,f_y0)为使|SI(x,y,ξ,η)|取最大值的(ξ,η)，同时也是干涉条纹图中每一点的横向与纵向条纹频率。

故条纹频率f_x0和f_y0的提取方法可以表述为：首先要将条纹图的条纹频率(ξ,η)在其取值范围内均匀离散化，然后把这些值分别带入|SI(x,y,ξ,η)|，其中能使|SI(x,y,ξ,η)|取最大值的点(f_x0,f_y0)即为公式(11)的解，也就是条纹频率。

3物体变形相位的计算

干涉条纹图一般可以表示为：

其中a(x,y,t)是背景光强，b(x,y,t)是条纹图的幅值，为条纹图的相位信息。将函数I(x,y,t)在点(x0,y0)做泰勒级数展开，取一级近似可得：

一般地认为a和b为常数。根据光流场理论的假设和(1)-(5)式，光流场基本等式(5)可表示为

即

其中为点(x₀,y₀)在Δt时间内相位的改变量。由于其中f_x0和f_y0代表了点(x₀,y₀)处的横向和纵向条纹频率。故式(15)又可写为

化简式(16)可得

将根据式(8)和式(9)采用迭代的方法得到的光流场运动矢量u和v以及公式(11)得到的条纹横向和纵向频率f_x0和f_y0带入公式(17)即可求得干涉条纹图的相位的变化量

本专利基于光流场理论，提出了一种由两幅条纹图像解调出物体变形相位的新方法。首先根据光流场理论获得两帧连续图像之间的光流矢量场，然后运用双幅图像中初始图像的条纹频率与该矢量场进行计算来获得物体的全场变形相位分布。本发明给出了光流场相位计算方法，对其进行了模拟和实验验证。模拟结果和实验结果表明，该算法能够将直观的运动矢量场和变形相位的提取联系起来，能够解调出物体变形相位信息，不需要相位解包络运算，过程简单方便，并且在条纹越密集处提取的变形相位信息越准确。该方法为计算物体全场变形相位分布和动态测量提供了新的途径。

4模拟和实验结果及分析

4.1模拟结果及分析

根据公式(10)，采用Matlab软件生成一幅干涉条纹图作为原始图像，大小为256*256像素，如图1所示，其中a＝0，b＝1，相位分布为

根据采样定理，一幅图像所能表达的最大条纹频率不超过0.5，故将频率(ξ,η)均匀离散化的范围设定为[-0.5,0.5]，抽样间隔取Δf＝0.01。积分窗口的σ_x和σ_y均取10个像素大小。在提取完模拟图的条纹频率(f_x0,f_y0)后，再对条纹图分别加入的固定相位和的变形相位。之所以选择变形相位的大小为是因为当两幅图像中某点的变形相位值大于π时，该点会出现黑白条纹的突变现象，从而导致所估算的运动矢量u、v以及最后所求的相位的值为零。故在本方法所能测得的变形相位的取值范围为[0,π]。当变形相位时，其变化的范围约为(0,π/3)，符合本方法的条件限制要求。附加相位后所得到的条纹图分别为图2(a)和图2(b)所示。

运用光流场理论分别对图1与图2(a)和图2(b)进行两幅图之间运动场的测量，光流场迭代次数为1000次，其中图1与图2(a),即附加了相位后所估算的x方向光流运动场在y＝127处截面图如图3所示。令两幅图像之间的时间变化量Δt＝1，最后根据公式(17)得到附加的相位，即变形相位信息，其理论和模拟的归一化结果对比如图4、图5(a)、图5(b)所示。

由图1、图2(a)和图2(b)可以看到，所有干涉条纹图的中央均为一个较大的圆斑，并且越靠近中心部分条纹越稀疏，其条纹频率也就越低。由于在条纹越密集处运动矢量越小，故在两幅图像之间所估计的光流运动场的大小应是由中心向两边单调递减的。但是正如图3所示，本算法在条纹频率较低处估算的光流场会产生较大的误差，尤其是在圆斑中心处所估算的运动场为0，所以在点a1(x＝80)和a2(x＝176)之间所估算的运动场会有较大的起伏，而在x<80和x>176处所估算的运动场较为准确。由于上述误差的存在，故在两模拟图的中央会出现一个“坑”一样的凹陷，如模拟结果中图4和图5(b)所示。经数据分析，a1和a2处所对应的条纹频率绝对值为0.3，即当条纹频率的绝对值大于0.3时所估算的光流运动场较为理想。在图4中条纹频率绝对值大于0.3处所提取的相位趋于一个平面，即附加的相位为一个常量，这与理论值是相符的。同时，根据图5(a)和图5(b)在条纹频率绝对值大于0.3处的理论与模拟结果对比可以看出，该方法提取的相位是准确的。

4.2实验数据处理结果及分析

采用迈克尔逊干涉获取干涉条纹，由CCD采集初始图像，图像大小为512×512像素，如图6(a)所示。用相移器对干涉条纹图附加π/7的固定相位，获得的条纹图如图6(b)所示。取Δt＝1，得到的实验结果如图7所示。

从图7可以看到，由于在实验当中图像之间无法保持各个像素点之间运动高度一致且有噪声因素的干扰，在条纹稀疏处估计的运动场误差较大，故在中央凹陷的周边会出现较大的起伏，而在条纹频率绝对值大于0.3的区域所测得的数据趋于平稳，即附加的相位趋于一个常量，故该结果表明此算法可行。通过模拟和实验结果可以证明，该方法在条纹越密集处获得的相位越准确，完全克服了现有的相位解调方法中条纹越密误差越大的缺陷。此外，若处理的图像为数幅连续变化的图像，只需令每两帧图像之间的时间变化量Δt取两帧图像之间的时间间隔即可完成干涉条纹图相位差的动态定量提取。

5结论

本专利提出了一种基于光流场理论的解调双幅图像之间相位差的算法，该方法根据光流场理论中的Horn-Schunck算法获取两幅图像之间的运动场，通过窗口傅里叶变换算法提取出初始图像条纹频率并与运动场进行计算，进而求解出两幅图像之间相位的变化量。这种方法只需在空间域当中用两幅条纹图就可以完成变形相位的提取，免去了转换到频域和相位解包络操作，过程简单、方便，对全息条纹变形场、电子散斑条纹场均适应。模拟和实验结果表明该算法对变化范围在[0,π]之间的变形相位的计算是准确的，且在条纹越密集的区域提取的相位信息越准确。此外，只要得知在视频序列中两帧图像间的时间间隔即可完成光干涉动态测量变形相位的提取。该方法为光干涉动态测量提供了新的途径。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。