烟线照片也可以玩AI？
本文首先发表于航华微信公众号，请点击访问
烟线这种老技术也能玩人工智能吗？是的，科研人员可以利用一些新颖的算法来分析烟线流场可视化图像，从中提取流动的动态信息来研究湍流结构。本文以32幅羽毛球烟线照片为例，利用本征正交分解（POD）来为图像进行智能分类，排序，构建羽毛球尾流的“低阶流动”，该模型可刻画流场的主要动态特征。拍几张照片就能复现动态特征？听起来有些不可思议！具体的过程如下：
1 准备工作
本文为公众号文章“烟线仪不仅能用来看旋涡，还能做POD分解！”的后续，请有兴趣的读者利用matlab软件（或Octave软件）依次执行该文中的五个步骤。
2 考察关键模态的形态
公众号文章“烟线仪不仅能用来看旋涡，还能做POD分解！”中通过分析各个模态的形态发现模态3与4代表着向下游传播的大尺度旋涡结构。这样一对模态通常形态相似但空间上有一定的错位，错位距离大约是1/4个波长（见下图）。本文以这对最主要模态为基础，构建流场的动态模型。
3 构建表述动态特征的复平面
将每幅图像向模态3,4投影，获得一对模态系数a3与a4。于是，每一幅烟线照片在一个复平面（实部a3，虚部a4）内都有了一个自己的位置：坐标就是（a3，a4）注：图中实虚部标反。
for image_no=1:32A=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 imagea_3(image_no)=proj_a(3)/lambda(3);a_4(image_no)=proj_a(4)/lambda(4);end
将所有图对应的复数a3+i a4的模画出来；选择一个阈值（本例中阈值为所有图像中最大模值的50%）。如果某图a3+i a4之模高于阈值，代表该图有明显大尺度结构出现（模态3,4对该图的贡献较大），相反，则表示模态3,4对该图的贡献微小。
Modulus=sqrt(a_3.^2+a_4.^2);Max_modulus=max(Modulus);figure, axes('position',[.2 .2 .5 .6])bar(Modulus/Max_modulus); hold on;xlabel(['Image #'],'fontsize',14);ylabel(['$\sqrt{a_3^{*2}+a_4^{*2}}/max(\sqrt{a_3^{*2}+a_4^{*2}})$'],'fontsize',14,'interpreter', 'latex');plot([032.5],[11]*.5,'r')axis([032.501])print -djpeg modulus_mode34
上图中红线代表阈值。图中可见图像5,8,10等13幅图有较大的模。
figure, axes('position',[.2 .2 .5 .6])for i=1:32mode_modulus=sqrt(max(a_3(i).^2+a_4(i).^2));if mode_modulus>Max_modulus*.5plot(a_3(i),a_4(i),'ob','markerfacecolor','b'); hold on;text(a_3(i),a_4(i),num2str(i),'color','r','fontsize',12)endendr=Max_modulus*.75;for i=0:2:360plot(r*cos(i*3.14/180),r*sin(i*3.14/180),'.k');endaxis equal;axis([-11 -11]/2);plot([-11]/2,[00],'-k');plot([00],[-11]/2,'-k');ylabel(['a^*_4=a_4/\lambda_4'],'fontsize',14);xlabel(['a^*_3=a_3/\lambda_3'],'fontsize',14);print -djpeg complex_field
把这些幅值（模）较大的图标记在a3，a4复平面上(见下图)，其中蓝点代表着每幅流线照片在复平面中的位置，红色数字是照片编号。黑色的圆圈仅仅是为了辅助观察。流场在做着周期性变化，每幅图像都是周期性变化的一个瞬间，图像散落在一个圆圈上，每个瞬间都对应一个辐角atan(a4/a3)。本文仅仅使用了32个数据样本（图像），如果有大量的样本，可以预见它们将沿圆圈分布，填补空白空间。注：图中缺俩根号。
我们可以考察一下复平面内距离较近的两点，比如19与23。对比图19与图23（见下图），它们的确有很多相似之处。读者可以自己对比在复平面内距离较近的其他样本。这是一个有实用意义的聚类方法，非常适用于周期性强的问题（比如尾流等）。该方法的更多信息可见 Oberleithner et al. （2011）。
4 “复现”旋涡的动态发展
最后是本文的重头戏，我们将根据辐角atan(a4/a3)构建一个流场模型，来“复现”旋涡的动态发展过程。
% create movieimg_sequence=[2319152629102117];figure,vidObj = VideoWriter('shedding.avi');vidObj.FrameRate = 5;vidObj.Quality = 75;open(vidObj);axis tight manualset(gca,'nextplot','replacechildren');for j=1:length(img_sequence)A=imread(['img_' num2str(img_sequence(j))],'jpeg');A=fliplr(A);imshow(A); %hold on;% Write each frame to the file.currFrame = getframe(gcf);writeVideo(vidObj,currFrame);endclose(vidObj);
00:04
该视频中对应的图像辐角变化见下图。视频中可见剪切层大幅度振荡生成大尺度旋涡，这些旋涡随即”脱落“并向下游发展。随即拍摄的几幅照片的确可以被”智能“地合成为动态图像。更重要的是，一对数值a3和a4可以相对准确的描述复杂流动的核心动态特征！不可思议吧
00:04
参考文献：
Oberleithner, K.,Sieber, M., Nayeri, C.N., Paschereit, C.O., Petz, C., Hege, H.C., Noack, B.R., etal. (2011), “Three-dimensional coherent structures in a swirling jet undergoing vortex breakdown: Stability analysis and empirical mode construction”, Journal of Fluid Mechanics, Vol. 679, pp. 383–414.
