单纯把“中心”相关计算改为“偏心”相关计算只是减少了一定相关窗口的操作,并没有减少计算量。而对于同一位置的Fw和Gw的所有点的相关计算可以根据卷积分定理和傅里叶变换(FFT)的相关定理得到:
(2.8)
其中,*为取共轭,F(U,V)是f(m,n)的傅里叶变换,G(U,V)是g(m,n)的傅里叶变换,只需对 作傅里叶逆变换,得到相关函数 (m,n),然后检测 (m,n)的峰值,即可得到对应粒子的位移。
基于FFT的互相关算法的计算量可以分为两个方面计算:
1) 一幅速度矢量场图像需要算出的位移s(i,j)的数目为: (2.9)
2) 因为快速傅里叶变换的运算量为O(Nlog2N)次。所以用式(2.8)获得计算一次位移所需的相关函数的计算量为:
(2.10)
这样,总的计算次数为C=Cs×CF,上节所述的图像所需的计算次数为O(108) ,又降低了两个数量级。所以采用基于FFT的互相关算法比普通算法快很多,所需的限制条件只是流体的最大速度在横向和纵向都不能超过wx/3和wy/3,这一点在实际操作中是可以满足的。因此,基于以上原因,作者选择此算法完成DPIV系统处理软件,以完成细水雾雾场的测量。
2.5 亚像素拟合
如果计算仅到此为止,则所获得的数据均为整数值,也即此时算法对小于一个像素的位移不响应。因此系统的位移的动态范围限制在1~wx(或wy)像素这一狭窄的范围内,而且得到的速度值皆为整数,量化效应明显,这种精度往往在实际应用中是不够的。为了获得“亚像素”精度,有必要对结果进一步拟合,得到最大相关峰值在两像素间的精确位置。考虑到示踪粒子成像时亮度的分布为高斯分布:
(2.11)
所以本文将采用高斯曲线拟合峰值。设峰值所在点的x轴上的坐标为xp,相关值为1,则其横向相邻两像素在x轴上的坐标分别为xp-1和xp+1,对应点的相关值为0、2,则在横向拟合,可得极值点的位置为:
(2.12)
对最大相关峰所在像素点进行横向、纵向、斜向(此时分母中的2改为 )四个方向进行拟合,最终确定最大峰值的精确位置。
2.6 数据后处理-误差矢量剔除
由于此算法中流体的最大速度在横向和纵向都不能超过wx/3和wy/3,所以在测量流场速度时,不仅要注意窗口大小的选取,还需在测量之后将超出限制范围的矢量剔除。
2.7 算法的计算机实现
本文采用矩阵预算语言MATLAB7.7.0作为软件开发和运行平台实现了上述算法,结合2.3和2.4两节的内容,FFT的互相关算法的步骤如图2.3(其中t为时间间隔):
图2.3 基于FFT的相关算法的步骤
下面介绍一下具体的运算过程(图2.4):首先输入需要进行互相关计算的连续两帧图像,然后选择两幅图的查询计算窗口大小和移动步长,接着对两幅图进行FFT快速傅里叶变换并得到结果,然后对结果在频率进行复数共轭相乘计算,对结果进行反傅里叶变化,求出相关函数值,并将零点平移到谐线中心,就是将FFT变换得到的相关谱结果的图像调整一下,因为FFT变换的结果是将低频分量放在了4个角上,而高频分量放在了正中间,在matlab中采用fftshift函数即可。我们对FFT的结果进行寻找,查找出模最大的点,该点与中心点的距离就是位移量。其中傅里叶变换及反变换可以通过matlab中的fft2函数和ifft2获得,共轭可通过conj函数完成。接下来用高斯曲线拟合精确确定相关峰的位置及位移,然后重复进行相关运算,直至全幅图像计算完毕,根据时间间隔获得流场速度,画出速度矢量图。 基于数字粒子图像处理的细水雾雾场速度测量(5):http://www.youerw.com/huaxue/lunwen_3049.html