本文為實驗流體力學課程作業，需用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');

轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/385450.html

標籤：其他

上一篇：實戰內容（13）- Invalid audio stream. Exactly one MP3 audio stream is required.

下一篇：在 Python 中使用 OpenCV 制作簡單影像影片

基于Matlab的PIV影像處理

目錄