基于Matlab的PIV图像处理

最新推荐文章于 2023-12-27 14:22:02 发布

qq_45885422

最新推荐文章于 2023-12-27 14:22:02 发布

阅读量1.6w

点赞数 41

文章标签： matlab 图像处理

本文链接：https://blog.csdn.net/qq_45885422/article/details/121996247

版权

本文为实验流体力学课程作业，需用Matlab处理两幅PIV图像来计算该流动的速度场，基于速度获得流线与涡量场信息。两幅图像如下。本文根据所查资料，自己编写了Matlab代码来处理，但Matlab工具箱中就有处理的PIVLab软件，处理效果也更好，因此本文仅提供一种思路方法。

1 背景知识

1.1 PIV测速原理

PIV全称Particle Image Velocimetry（平面粒子成像测速），通过激光器和其他光学元件得到片光源，用其照亮掺有示踪粒子的待测流体的待测片层（尽量使该片层内流体的流速都在该片层内），利用数码照相机连续拍摄被照明的区域，得到一系列一定时间间隔下示踪粒子的位置，使用图像处理技术计算相同粒子的位移，得到当地流速。

流场中某一示踪粒子在二维平面上运动,其在x、y两个方向上的位移随时间的变化为x(t)、y(t)，是时间t的函数。那么,当两次曝光时间间隔∆t足够小时，该示踪粒子所在处的质点的二维流速可以表示如下公式：

$u=\frac{dx(t)}{dt}\approx \frac{x(t+\Delta t)-x(t)}{\Delta(t) }=\frac{\Delta x}{\Delta t}$

$v=\frac{dy(t)}{dt}\approx \frac{y(t+\Delta t)-y(t)}{\Delta(t) }=\frac{\Delta y}{\Delta t}$

式中:u与v是沿x方向与y方向的瞬时速度，∆x/∆t和∆y/∆t是沿x方向与y方向的平均速度, ∆t是测量的时间间隔。

1.2 互相关理论

因为粒子小，所以很难从两幅图像中分辨同一个粒子，也就无法获得所需的相对位移。因此利用互相关分析理论。图像采集系统获得的每一对图像都是从相同的空间位置上得到的，且曝光的时间间隔可以作为已知参数。流场中的示踪粒子反射来自片光源的光线，每一粒子上反射的光强信号与其空间位置成单一映射,这就形成光强信号与空间位置的函数映射关系，使用互相关分析方法可以确定两幅图像之间的对应关系。

在对采集图像进行分析时，首先需要明确一个概念 “判读区” :它是指图像中一定位置取一定尺寸的方形图，通过对判读区进行信号处理，就可以获取速度。假设系统在 $t_{0}$ 以及 $t_{0}+\Delta t$ 这两个时刻分别获取图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)与 $r_{pq(\tau _{x},\tau _{y})}$ 有如下公式：

$\Delta x=i+\frac{ln(r_{pq(i-1,j))})-ln(r_{pq(i+1,j))})}{2ln(r_{pq(i-1,j))})-4ln(r_{pq(i,j))})+2ln(r_{pq(i+1,j))})}$

$\Delta y=j+\frac{ln(r_{pq(i,j-1))})-ln(r_{pq(i,j+1))})}{2ln(r_{pq(i,j-1))})-4ln(r_{pq(i,j))})+2ln(r_{pq(i,j+1))})}$

在以上的公式中，表示互相关函数的最大值，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)
for 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),n,1)+velocity((m+1),(n+1),1))/8;
%             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),n,2)+velocity((m+1),(n+1),2))/8;
            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 绘制图像

根据得到的速度和涡量绘制流场、流线和涡量图。

%绘制流场
figure(1);
ax1=axes('linewidth',3, 'box', 'on', 'FontSize',12);
set(gcf,'color','white');
set(ax1,'YTick', []);
set(ax1,'XTick', []);

for m=1:1:(cycw+1)
    for n=1:1:(cycl+1)
        quiver((stp*n)*0.09,stp*(cycw-m+1)*0.09,velocity(m,n,1),velocity(m,n,2),'Color','k','MaxHeadSize',1,'AutoScaleFactor',1);
        hold on;
    end
end
xlabel('X(mm)','FontWeight','bold'); 
ylabel('Y(mm)','FontWeight','bold');
axis([0 143 0 106])

% 绘制流线
figure(2);
ax1=axes('linewidth',1,'box','on','FontSize',12);
set(gcf,'color','white');
set(ax1,'YTick', [0:20:106]);
set(ax1,'XTick', [0:20:143]);
xlabel('X(mm)','FontWeight','bold'); 
ylabel('Y(mm)','FontWeight','bold');
axis([0 143 0 106])
streamlines=flipud(imresize(velocity,1.456));
streamslice(streamlines(:,:,1),streamlines(:,:,2),5);

%绘制涡量
figure(3)
ax1=axes('linewidth',1,'box','on','FontSize',12);
set(gcf,'color','white');
xlabel('X(mm)','FontWeight','bold'); 
ylabel('Y(mm)','FontWeight','bold');
F=imresize(vorticity,1.456);
contourf(F,10)
colorbar;
xlabel('X(mm)','FontWeight','bold'); 
ylabel('Y(mm)','FontWeight','bold');

绘制得到的图像如下：

3 完整代码

以下只列出主程序代码，各调用函数已经在上部分中给出。

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);

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

velocity=velo_correction(velocity,cycw,cycl,3);%速度修正，3为修正次数

%绘制流场
figure(1);
ax1=axes('linewidth',3, 'box', 'on', 'FontSize',12);
set(gcf,'color','white');
set(ax1,'YTick', []);
set(ax1,'XTick', []);

for m=1:1:(cycw+1)
    for n=1:1:(cycl+1)
        quiver((stp*n)*0.09,stp*(cycw-m+1)*0.09,velocity(m,n,1),velocity(m,n,2),'Color','k','MaxHeadSize',1,'AutoScaleFactor',1);
        hold on;
    end
end
xlabel('X(mm)','FontWeight','bold'); 
ylabel('Y(mm)','FontWeight','bold');
axis([0 143 0 106])

% 绘制流线
figure(2);
ax1=axes('linewidth',1,'box','on','FontSize',12);
set(gcf,'color','white');
set(ax1,'YTick', [0:20:106]);
set(ax1,'XTick', [0:20:143]);
xlabel('X(mm)','FontWeight','bold'); 
ylabel('Y(mm)','FontWeight','bold');
axis([0 143 0 106])
streamlines=flipud(imresize(velocity,1.456));
streamslice(streamlines(:,:,1),streamlines(:,:,2),5);

%计算涡量
vorticity=vortex_calculate(velocity,cycw,cycl);

figure(3)
ax1=axes('linewidth',1,'box','on','FontSize',12);
set(gcf,'color','white');
xlabel('X(mm)','FontWeight','bold'); 
ylabel('Y(mm)','FontWeight','bold');
F=imresize(vorticity,1.456);
contourf(F,10)
colorbar;
xlabel('X(mm)','FontWeight','bold'); 
ylabel('Y(mm)','FontWeight','bold');