《粒子图像测速技术.doc》由会员分享，可在线阅读，更多相关《粒子图像测速技术.doc（6页珍藏版）》请在金锄头文库上搜索。
1、热灾害实验技术粒子图像测速技术 实验目的：测量流场的速度分布实验装置：摄像头，片光源，计算机实验原理与方法：互相关算法：根据DPIV互相关算法理论，过一采样窗口，同一特定位置，可顺序获取两幅数字图像。粒子的位移可从一幅图像到另一幅相对应的图像经互相关计算获得。两幅顺序获取图像中相同位置的两个查问窗口函数f(m，n)和g(m，n)，从信号系统的观点出发，g(m，n)可以看作是f(m，n)经线性转换后叠加以噪声而成，如图所示。1、 预先在流场中撒入一定浓度和大小的示踪粒子，用脉冲激光片光源照明二维流场，在垂直片光的方向上摄像，获得两次或者多次曝光的粒子图像。2、 对采取的图像使用Photoshop
2、进行处理使之变成长宽为256x256的灰度图像，使用matlab和互相关算法编程，对图像进行处理，得到二维流场中速度矢量分布图。实验结果：隔两帧得到的灰度图： 得到速度矢量图如图所示： 得到速度矢量图如图所示：源程序：function XNum,YNum,AutoFlag=g(a1,a2,OutFile)fid1=fopen(a1.jpg,r);fid2=fopen(a2.jpg,r);fid3=fopen(OutFile.jpg,w);XL=256;YL=256;frewind(fid1);Image1=fread(fid1,XL,YL);frewind(fid2);Image2=fread
3、(fid2,XL,YL);WinSizeX=64;WinSizeY=64;DeltaX=16;DeltaY=16;XNum = (XL-WinSizeX)/DeltaX;YNum = (YL-WinSizeY)/DeltaY;a=3.4968e+006 +4.7425e-010i;frewind(fid3);for j1=1:YNum+1; for i1=1:XNum+1; CorrWin1=Image1(i1-1)*DeltaX+1:(i1-1)*DeltaX+WinSizeX,(j1-1)*DeltaY+1:(j1-1)*DeltaY+WinSizeY); CorrWin2=Image2(
4、i1-1)*DeltaX+1:(i1-1)*DeltaX+WinSizeX,(j1-1)*DeltaY+1:(j1-1)*DeltaY+WinSizeY); F1=fft2(CorrWin1); F2=fft2(CorrWin2); S=F2.*conj(F1); Pks=ifft2(S); Pks=fftshift(Pks); Pks=Pks/a; if (i1=2) SURF(abs(Pks); pause end temp1,maxx=max(Pks); temp2,maxy=max(temp1); if temp2=0.0 | abs(maxx(maxy)-WinSizeX/2-1)
5、WinSizeX/3.0 | abs(maxy-WinSizeY/2-1)WinSizeY/3.0 VX(i1,j1)=0; VY(i1,j1)=0; Corr(i1,j1)=0; else TX=maxx(maxy); TY=maxy; Corr(i1,j1)=temp2;%-用高斯拟合峰值的精确位置- y0=Pks(TX-1,TY); %横向 y1=Pks(TX,TY); y2=Pks(TX+1,TY); CX1=(y0-y2)/(2.0*(y0+y2-2.0*y1); y0=Pks(TX,TY-1); %纵向 y1=Pks(TX,TY); y2=Pks(TX,TY+1); CY1=(y0
6、-y2)/(2.0*(y0+y2-2.0*y1); y0=Pks(TX-1,TY-1); %第一对角线 y1=Pks(TX,TY); y2=Pks(TX+1,TY+1); CX2=(y0-y2)/(2.828*(y0+y2-2.0*y1); CY2=(y0-y2)/(2.828*(y0+y2-2.0*y1); y0=Pks(TX+1,TY-1); %第二对角线 y1=Pks(TX,TY); y2=Pks(TX-1,TY+1); CX3=-(y0-y2)/(2.828*(y0+y2-2.0*y1); CY3=(y0-y2)/(2.828*(y0+y2-2.0*y1);%- VX(i1,j1)=T
7、X+CX1+CX2+CX3-WinSizeX/2-1; VY(i1,j1)=TY+CY1+CY2+CY3-WinSizeY/2-1; end X(i1,j1)=WinSizeX/2+(i1-1)*DeltaX; Y(i1,j1)=WinSizeY/2+(j1-1)*DeltaY; %if(abs(Corr(i1,j1)abs(a) % a=Corr(i1,j1); %end; fprintf(fid3,%8.2f,%8.2f,%12.6f,%12.6f,%8.2fn,X(i1,j1),Y(i1,j1),VX(i1,j1),VY(i1,j1),Corr(i1,j1); end %j1end%aquiver(X,Y,VX,VY);fclose(fid1);fclose(fid2);fclose(fid3);cx=5;cy=5;sigma=0;for j1=1:YNum+1, for i1=1:XNum+1, sx=cx*abs(VX(i1,j1)-cx)/(cx*cx+cy*cy); sy=cy*abs(VY(i1,j1)-cy)/(cx*cx+cy*cy); sigma=sigma+sqrt(sx*sx+sy*sy); endendsigma=sigma*100/(XNum+1)*(YNum+1); %百分比s,errmsg=sprintf(sigma=%fn,sigma);s
