烟线仪不仅能用来看旋涡，还能做POD分解！
“眼见为实”。看得见的东西能给人留下深刻印象，所以烟线这个“老”技术依然能够在现代流体力学、空气动力学教学与科研中站住脚。其实，烟线不仅能做定性研究，也可以做定量研究。如果使用得当，该技术能揭示出流场中许多关键信息。本文中，我们将利用30度迎角羽毛球周边流场的32幅烟线照片为例，给感兴趣的读者介绍如何进行图像定量分析。
我们将利用正交分解（Proper Orthogonal Decomposition，POD）方法来提取湍流大尺度结构的特征。关于POD分解的数学基础，读者可以参阅Kutz[1]。主要过程包括1准备工作、2求平均、3给每幅图像去平均、4构建相关系数矩阵，求矩阵的特征值与特征向量、5将每幅图像向特征向量投影，获得投影系数（模态系数），最后6利用人工选择的若干投影系数及对应的特征向量重构图像，获得所谓的低阶图像。
1. 准备工作:
点击这里进入下载页面来下载全部32图像文件，储存在同一目录下(比如d:\smokewire_images)。打开Matlab或Octave软件，进入图像目录并创建脚本。键入
cd d:\smokewire_imagesedit compute_pod.m
在compute_pod.m脚本文件中键入下列初始化设置
image_width=160; %pixelsimage_height=100; %pixelsnumber_of_images=32; %number of images to analyzekey_modes=16; %number of modes of interestesData_mean=zeros(image_height,image_width); % initialize a matrix for the time-mean image
显示原始图像：
figure,for image_no = 1:number_of_imagesA=imresize(imread(['img_' num2str(image_no)],'jpeg'),[image_height,image_width]);Image_R=A(:,:,1); %only use the R information out of RGBif image_no < 7subplot(3,2,image_no),imshow(Image_R);title(['图像' num2str(image_no)]);endData_mean=Data_mean+double(Image_R);endprint -djpeg raw_images
2. 获取平均图像:
Data_mean=Data_mean/number_of_images;Image_mean=int8(Data_mean);figure, imshow(Image_mean);title('图像均值')print -djpeg mean_image
3. 图像去均值：
下面为每幅流显图像减去均值，这样可以突出旋涡引起变化的这一部分。这一部分也是我们下一步进行正交分解的对象。同时，我们建立一个容纳所有图像的大矩阵B，每一行存放一个图像，总行数为number_of_images，总列数为image_width*image_height。注意计算时需要使用双精度double格式，图像显示时数值则需要调整为整数int8格式。
B=zeros(number_of_images,image_height*image_width);figure,for image_no=1:number_of_images% load image file, image will be load into a row*column*3 arrays, each2D array represents RGB information,respectively.A=imresize(imread(['img_' num2str(image_no)],'jpeg'),[image_height,image_width]);Image_R=A(:,:,1); %use only the R information out of RGBsImage_mean_subtracted=int8(Image_R)-int8(Image_mean); %subtract the mean imageif image_no<7subplot(3,2,image_no),imshow(Image_mean_subtracted);title(['去均值后的图像' num2str(image_no)]);endData_mean_subtracted=double(Image_R)-double(Image_mean); %subtract the mean image data, forming a de-meaned data arraya=reshape(Data_mean_subtracted,1,image_height*image_width); %reshape image data to a single columnB(image_no,:)=a; % store each image data in a big array B, each image as a columnendprint -djpeg raw_images_de_meaned
在去除均值以后的图像中（下图），羽毛球变得非常模糊，但其轮廓在每幅图中依稀可见，且颜色深浅不同。这说明每次曝光强度都有所不同。不同的曝光强度将使每幅图的明暗出现微弱变化。这与闪光灯的特性有关，也与每次曝光过程中烟的浓度有关：不同浓度形成的光散射不同。
4. 构建相关系数矩阵并计算特征值、特征向量:
我们利用经典正交分解方法来构建相关矩阵C：
矩阵C为对称矩阵，行数与列数均为 image_width*image_height 。本例中C的大小为16000*16000。利用eigs命令计算C矩阵主要的特征值lambda与特征向量V。因总特征向量数量为16000，我们仅仅关注对应特征值最大的key_modes个特征向量以减少计算量。流体力学领域同行通常称这些特征向量代表着流场变化的“主要模态”。
C=B'*B; % construct a core correlation matrix[V,D]=eigs(C,key_modes,'lm'); % compute the eigenvectors and eigenvalues of the core matrixlambda=sqrt(diag(D)); % pick the eigenvalues out of the diagonals of the eigenvalue matrix
显示主要模态：
在脚本中键入下面命令显示前9阶模态
figure,for mode=1:9mode_shape=int8(reshape(V(:,mode)*lambda(mode),image_height,image_width)); %reshape the eigenvectors from a single column into image matrixsubplot(3,3,mode),imshow(mode_shape); %show eigenvectors as imagetitle(['模态' num2str(mode)]);endprint -djpeg mode_shapes
需要注意的是：（1）模态1,2中羽毛球的轮廓依然可见，而在其他模态中几乎看不见羽毛球轮廓。这说明模态1,2很可能代表了每次曝光的光强变化；（2）模态3-9剪切层中有大尺度明暗区域，尤其模态模态3、4明暗区域左右错开，错开距离约1/4波长。说明模态3、4很可能对应着向下游传播的大尺度结构。
读者可以检验模态之间的正交性：
V(:,1)'*V(:,1)V(:,1)'*V(:,2)
上面第一行的结果为1，第二行为0。说明模态之间的确是正交的。用代数的语言来说: 上面这些模态图像形成一组标准正交基。
4投影：不同模态对每幅图像的贡献
下一步，我们以图像5为例，考察不同模态对该图的贡献。我们将图像5投影到各个模态上（既然各个模态是正交的，我们就可以像下图那样：想象有这样一个多维的坐标系，每一个轴就是一个模态，图像在每个模态上的投影的结果就是它在这个模态上的坐标值，流体力学领域俗称模态系数，代表了这个模态对该图的贡献量）。
投影的方式为(图像向量a尺寸为1*16000)：
proj_a=a*V;计算前6个图像模态系数完整的程序为：
figure,for image_no=1:6A=imresize(imread(['img_' num2str(image_no)],'jpeg'),[image_height,image_width]);Image_R=A(:,:,1);Image_mean_subtracted=int8(Image_R)-int8(Image_mean);Data_mean_subtracted=double(Image_R)-double(Image_mean);a=reshape(Data_mean_subtracted,1,image_height*image_width);proj_a=a*V; %project image onto the eigenvectors, obtain a 1x32(number_of_images) array, representing the contributions from each mode to this imagesubplot(3,4,image_no*2-1),imshow(Image_mean_subtracted); % show image(de-meaned)title(['去均值后的图像' num2str(image_no)]);subplot(3,4,image_no*2),bar(proj_a(1:8)), yticklabels({''}),axis([09-15001500]); %show modal contribution to this imageylabel(['图像' num2str(image_no) '系数']); xlabel(['模态']);endprint -djpeg mode_coefficients
5. 图像重构与低阶模型：
就像一个向量可以用坐标值与坐标轴单位向量表示一样
图像也可以用模态系数与模态表示：
proj_a*V';利用全部、及部分模态来重构图像5的完整程序为：
figure,for image_no=5A=imresize(imread(['img_' num2str(image_no)],'jpeg'),[image_height,image_width]);Image_R=A(:,:,1);Image_mean_subtracted=int8(Image_R)-int8(Image_mean);Data_mean_subtracted=double(Image_R)-double(Image_mean);a=reshape(Data_mean_subtracted,1,image_height*image_width);subplot(3,2,1),imshow(int8(reshape(a,image_height,[]))); title(['去掉均值的图像5']); %show original image 5proj_a=a*V; %projecting image 5 to each modeA2=proj_a*V'; %reconstruct image 5 using all the modessubplot(3,2,2),imshow(int8(reshape(A2,image_height,[])));title(['全部模态重构的图像5']); % show reconstructed image 5using all modesA2=proj_a(1:2)*V(:,1:2)'; %reconstruct image 5 using modes 1,2subplot(3,2,3),imshow(int8(reshape(A2,image_height,[])));title(['模态1,2重构的图像5']);% show reconstructed image 5using modes 1,2A2=proj_a(3:8)*V(:,3:8)'; %reconstruct image 5 using modes 3-8subplot(3,2,4),imshow(int8(reshape(A2,image_height,[])));title(['模态3-8重构的图像5']);% show reconstructed image 5A2=proj_a(3:6)*V(:,3:6)'; %reconstruct image 5 using modes 3-6subplot(3,2,5),imshow(int8(reshape(A2,image_height,[])));title(['模态3-6重构的图像5']);% show reconstructed image 5A2=proj_a(3:4)*V(:,3:4)'; %reconstruct image 5 using modes 3,4subplot(3,2,6),imshow(int8(reshape(A2,image_height,[])));title(['模态3,4重构的图像5']);% show reconstructed image 5endprint -djpeg mode_reconstruction
图中可见：（1）前16个模态重构的图像已经与原图相差无几；2）1,2模态重构图像突出背景光的变化，与流动关系较小；（3）3,4模态重构出的图像已经可以大致描述大尺度结构的位置，5-8模态则添加了其他细节。模态3与4的模态系数可作为判断最大尺度运动出现、传播的标准，可被称为是一对可靠的低阶流动模型（low-order flow model）。
总结：
烟线流场图像不仅能用于定性研究，也可用于定量分析。我们希望读者能从本文得到启发，将正交分解方法、流场可视化方法应用于更广阔的空间去。欢迎大家利用本文中的matlab脚本程序。最后再次提醒读者点点击这里进入下载页面下载图像数据及程序。也欢迎您关注航华微信公众号，我们是骄傲的国产湍流仪器制造商。
文献：
[1] Kutz JN (2013) Data-Driven Modeling & Scientific Computation, Oxford University Press，第2.5节
