基于Matlab的PIV图像处理
本⽂为实验流体⼒学课程作业，需⽤Matlab处理两幅PIV图像来计算该流动的速度场，基于速度获得流线与涡量场信息。两幅图像如下。本⽂根据所查资料，⾃⼰编写了Matlab代码来处理，但Matlab⼯具箱中就有处理的PIVLab软件，处理效果也更好，因此本⽂仅提供⼀种思路⽅法。
⽬录
1 背景知识
1.1 PIV测速原理
PIV全称Particle Image Velocimetry（平⾯粒⼦成像测速），通过激光器和其他光学元件得到⽚光源，⽤其照亮掺有⽰踪粒⼦的待测流体的待测⽚层（尽量使该⽚层内流体的流速都在该⽚层内），利⽤数码照相机连续拍摄被照明的区域，得到⼀系列⼀定时间间隔下⽰踪粒⼦的位置，使⽤图像处理技术计算相同粒⼦的位移，得到当地流速。
流场中某⼀⽰踪粒⼦在⼆维平⾯上运动，其在x、y两个⽅向上的位移随时间的变化为x(t)、y(t)，是时间t的函数。那么，当两次曝光时间间隔∆t⾜够⼩时，该⽰踪粒⼦所在处的质点的⼆维流速可以表⽰如下公式：
式中:u与v是沿x⽅向与y⽅向的瞬时速度，∆x/∆t和∆y/∆t是沿x⽅向与y⽅向的平均速度,∆t是测量的时间间隔。
1.2 互相关理论
因为粒⼦⼩，所以很难从两幅图像中分辨同⼀个粒⼦，也就⽆法获得所需的相对位移。因此利⽤互相关分析理论。图像采集系统获得的每⼀对图像都是从相同的空间位置上得到的，且曝光的时间间隔可以作为已知参数。流场中的⽰踪粒⼦反射来⾃⽚光源的光线，每⼀粒⼦上反射的光强信号与其空间位置成单⼀映射，这就形成光强信号与空间位置的函数映射关系，使⽤互相关分析⽅法可以确定两幅图像之间的对应关系。
在对采集图像进⾏分析时，⾸先需要明确⼀个概念 “判读区” ：它是指图像中⼀定位置取⼀定尺⼨的⽅形图，通过对判读区进⾏信号处理，就可以获取速度。假设系统在以及这两个时刻分别获取图1和图2，在图1和图2中相同位置获取两个同样尺⼨⼤⼩的判读区f(m,n)以及g(m,n)，(m,n)表⽰ f与g分别在图1与图2中的相对位置，对f与g进⾏处理就可以获得此判读区对应位移s。
⽤傅⽴叶变换实现互相关分析，直接利⽤互相关函数的定义进⾏互相关函数的计算，计算量是⾮常巨⼤的，其⼀维的计算复杂性为O(N2)，计算量会随着序列长度的增⼤⽽急剧增长。为了提⾼运算效率，⼀般通过⼆维快速Fourier变换(FFT)在PIV技术中实现互相关函数的计算。
p(x,y)、q(x,y)与有如下公式：
在以上的公式中，表⽰互相关函数的最⼤值，i和j 表⽰互相关函数取最⼤值是相应的坐标，、、
和分别表⽰互相关函数数组中与最邻近的周围四个⽹格点上的互相关函数的数值。上述公式类似于⼀种在图像处理算法中通常采⽤的邻域运算，即互相关函数峰值位置不但和本⾝有关，⽽且还受到其相邻⽹格点的互相关函数数值的影响。
2 代码实现
2.1 预处理
读⼊图像，设置判断窗⼝，初始化速度。
clc ; clear all;
I1 = imread('xxx');
I2 = imread('xxx');
[row,col]=size(I1);
wsz=32;%设置判断窗⼝
stp=wsz/2;
cycl=floor((col-wsz)/stp);
cycw=floor((row-wsz)/stp);
velocity=zeros((cycw+1),(cycl+1),2);
2.2 速度计算
根据每个分区计算速度并保存
for m=1:1:(cycw+1)
for n=1:1:(cycl+1)
image1=I1((1+stp*(m-1)):(stp*(m-1)+wsz),(1+stp*(n-1)):(stp*(n-1)+wsz));
image2=I2((1+stp*(m-1)):(stp*(m-1)+wsz),(1+stp*(n-1)):(stp*(n-1)+wsz));
[x,y]=velo_calculate(image1,image2,stp); %计算粒⼦的u、v速度分量
velocity(m,n,1)=x;  %将粒⼦的u、v速度分量分别顺序存进velocity
velocity(m,n,2)=y;
end
end
其中的velo_calculate函数如下，使⽤快速傅⾥叶变换和⾼斯拟合修正。
注意这⾥并没有得到真正的速度，还需除以两幅图像之间的时间差，因此这⾥默认两幅图像的时间差为单位时间，但是实际中是远⼩于单位时间的，后续除去即可。
function[x,y]=velo_calculate(i1,i2,st)
r = fftshift(ifft2(fft2(i1).*conj(fft2(i2))));%计算互相关度
[x,y]=find(r==max(max(r)));
try  %⾼斯拟合修正
detax=(log(r(x-1,y))-log(r(x+1,y)))/(2*log(r(x-1,y))-4*log(r(x,y))+2*log(r(x+1,y)));
detay=(log(r(x,y-1))-log(r(x,y+1)))/(2*log(r(x,y-1))-4*log(r(x,y))+2*log(r(x,y+1)));
catch exception
detax=zeros(size(x));
detay=zeros(size(y));
end
x=x-st-1;
y=y-st-1;
d2=x.*x+y.*y;
[d2min,wh]=min(d2); %若互相关矩阵中有多个元素同为最⼩，则取位移量最⼩的那⼀个
x=x(wh)+detax(wh);
y=y(wh)+detay(wh);
temp=(-y);
y=x;
x=temp;
if (x^2+y^2)^0.5>8 %后处理1，Vector Validation，根据结果筛除过⼤的速度，速率上限选为8  x=0;
y=0;
end
2.3 速度修正
velocity=velo_correction(velocity,cycw,cycl,3);%速度修正，3为修正次数
其中，velo_correction函数如下：
function[velo]=velo_correction(v,cw,cl,t)
fontweight属性boldfor ce=1:t %进⾏t次拟合修正速度
for m=2:1:cw
for n=2:1:cl
refu=(v((m-1),(n-1),1)+2*v((m-1),n,1)+v((m-1),(n+1),1)+2*v(m,(n-1),1)+2*v(m,(n+1),1)+v((m+1),(n-1),1)+2*v((m+1),n,1)+v((m+1),(n+1),1))/12;
refv=(v((m-1),(n-1),2)+2*v((m-1),n,2)+v((m-1),(n+1),2)+2*v(m,(n-1),2)+2*v(m,(n+1),2)+v((m+1),(n-1),2)+2*v((m+1),n,2)+v((m+1),(n+1),2))/12;
%            refu=(velocity((m-1),(n-1),1)+velocity((m-1),n,1)+velocity((m-1),(n+1),1)+velocity(m,(n-1),1)+velocity(m,(n+1),1)+velocity((m+1),(n-1),1)+velocity((m+1), %            refv=(velocity((m-1),(n-1),2)+velocity((m-1),n,2)+velocity((m-1),(n+1),2)+velocity(m,(n-1),2)+velocity(m,(n+1),2)+velocity((m+1),(n-1),2)+velocity((m+1),            if abs(v(m,n,1)-refu)>(0.5*abs(refu))||abs(v(m,n,2)-refv)>(0.5*abs(refv))  %不合理标准选为任⼀⽅向速度与周围的平均值相差超过50%
v(m,n,1)=refu;
v(m,n,2)=refv;
end
end
end
m=1;
for n=2:cl
refu=(v(m,(n-1),1)+v(m,(n+1),1)+v((m+1),(n-1),1)+v((m+1),n,1)+v((m+1),(n+1),1))/5;
refv=(v(m,(n-1),2)+v(m,(n+1),2)+v((m+1),(n-1),2)+v((m+1),n,2)+v((m+1),(n+1),2))/5;
if abs(v(m,n,1)-refu)>(0.5*abs(refu))||abs(v(m,n,2)-refv)>(0.5*abs(refv))  %不合理标准选为任⼀⽅向速度与周围的平均值相差超过50%
v(m,n,1)=refu;
v(m,n,2)=refv;
end
end
m=cw+1;
for n=2:cl
refu=(v((m-1),(n-1),1)+v((m-1),n,1)+v((m-1),(n+1),1)+v(m,(n-1),1)+v(m,(n+1),1))/5;
refv=(v((m-1),(n-1),2)+v((m-1),n,2)+v((m-1),(n+1),2)+v(m,(n-1),2)+v(m,(n+1),2))/5;
if abs(v(m,n,1)-refu)>(0.5*abs(refu))||abs(v(m,n,2)-refv)>(0.5*abs(refv))  %不合理标准选为任⼀⽅向速度与周围的平均值相差超过50%
v(m,n,1)=refu;
v(m,n,2)=refv;
end
end
n=1;
for m=2:cw
refu=(v((m-1),n,1)+v((m-1),(n+1),1)+v(m,(n+1),1)+v((m+1),n,1)+v((m+1),(n+1),1))/5;
refv=(v((m-1),n,2)+v((m-1),(n+1),2)+v(m,(n+1),2)+v((m+1),n,2)+v((m+1),(n+1),2))/5;
if abs(v(m,n,1)-refu)>(0.5*abs(refu))||abs(v(m,n,2)-refv)>(0.5*abs(refv))  %不合理标准选为任⼀⽅向速度与周围的平均值相差超过50%
v(m,n,1)=refu;
v(m,n,2)=refv;
end
end
n=cl+1;
for m=2:cw
refu=(v((m-1),(n-1),1)+v((m-1),n,1)+v(m,(n-1),1)+v((m+1),(n-1),1)+v((m+1),n,1))/5;
refv=(v((m-1),(n-1),2)+v((m-1),n,2)+v(m,(n-1),2)+v((m+1),(n-1),2)+v((m+1),n,2))/5;
if abs(v(m,n,1)-refu)>(0.5*abs(refu))||abs(v(m,n,2)-refv)>(0.5*abs(refv))  %不合理标准选为任⼀⽅向速度与周围的平均值相差超过50%
v(m,n,1)=refu;
v(m,n,2)=refv;
end
end
end
velo=v;
end
2.4 计算涡量
根据涡量表达式，采⽤中⼼差分的⽅法计算涡量的⼤⼩，并把每个点的涡量⼤⼩与周围⼋个点的涡量的平均值进⾏⽐较，误差超过
50%的点位，直接取平均值。涡量处理有多种⽅法，此⽅法处理效果不是很好。
vorticity=vortex_calculate(velocity,cycw,cycl);
其中，vortex_calculate函数如下：
function[vor]=vortex_calculate(velocity_,cw,cl)
u=velocity_(:,:,1);
v=velocity_(:,:,2);
dv=zeros(cw+1,cl+1);
du=zeros(cw+1,cl+1);
vor=zeros(cw+1,cl+1);
for m=1:1:cw+1
for n=2:1:cl
dv(m,n)=(v(m,n+1)-v(m,n-1))/2;
end
dv(m,1)=(v(m,2)-v(m,1));
dv(m,cl+1)=v(m,cl+1)-v(m,cl);
end
for m=2:1:cw
for n=1:1:cl+1
du(m,n)=(u(m+1,n)-u(m-1,n))/2;
end
du(1,n)=(u(2,n)-u(1,n));
du(cw+1,n)=u(cw+1,n)-u(cw,n);
end
for m=1:1:cw+1
for n=1:1:cl+1
vor(m,n)=(dv(m,n)-du(m,n));
if abs(vor(m,n))<0.2||n>60
vor(m,n)=0;
end
end
end
vor=flipud(-vor);
for m=2:1:cw
for n=2:1:cl
refu=(vor((m-1),(n-1))+vor((m-1),n)+vor((m-1),(n+1))+vor(m,(n-1))+vor(m,(n+1))+vor((m+1),(n-1))+vor((m+1),n)+vor((m+1),(n+1)))/8;            if abs(vor(m,n)-refu)>(0.3*abs(refu))  %不合理标准选为任⼀⽅向速度与周围的平均值相差超过50%
vor(m,n)=refu;
end
end
end
for m=1:1:cw+1
for n=1:1:cl+1
if abs(vor(m,n))<0.1||n>60
vor(m,n)=0;
end
end
end
2.5 绘制图像
根据得到的速度和涡量绘制流场、流线和涡量图。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。