基于压缩感知的高频响流场重构方法及其应用

陈子玉，周楷文，刘应征，温 新

(上海交通大学 机械与动力工程学院，上海 200240)

0 引 言

流场的时空变化包含了重要的流动信息，但实际测量却难以同时得到较高时间分辨率和空间分辨率的流场。粒子图像测速（PIV）[1-4]诞生于20世纪80年代，是一种可获得高空间分辨率的光学流场测试技术，但其时间分辨率通常限于15 Hz以下。尽管利用高速相机和激光器可以获得相对较高时间分辨率的PIV流场信息，但价格昂贵，且有时受实验条件限制，而实际应用中高频流场非常常见，因此有必要发展高频PIV流场重构方法。

为了提高PIV测量流场的时间分辨率，国内外诸多学者进行了相关研究。Druault等[5]对所测得的PIV数据进行本征正交分解，然后直接对模态系数进行样条插值，从而得到高时间分辨率的流场。这种方法鲁棒性好且易于实施，但这种方法本身依赖于高时间分辨率采样，频率提高有限。之后，许多研究者借用具有高时间分辨率的局部点测量技术，如热线、麦克风等，来提高PIV的时间分辨率。He等[6]在使用PIV 测量技术的同时，附加了34个局部点测量的高频采样传感器，结合线性随机估计和数据同化，重构出时间分辨率的湍流流场，但这种方法需要解控制方程，计算量较大，且在某些工况时，无法安装大量的高频采样传感器。Jin[7]等使用双向循环神经网络，并借助于高频局部点测量和PIV测量，恢复出时间分辨率的圆柱绕流流场，但这种方法仍需要消耗较大计算量，且加入的局部测量点同样会扰乱流场。黄路[8]针对非定常流场，结合互相关算法和卡尔曼滤波，利用仿真验证了时间分辨率实时自适应调整算法的有效性，但这种方法依赖于纳秒级分幅成像的硬件装置。因此，目前针对该问题的研究存在许多缺陷，如仅使用PIV测量的时间分辨率重构依赖于自身的高时间分辨率测量或依赖于昂贵的硬件，而结合侵入式点测量的重构在某些工况下甚至无法安装点测量仪器或会扰乱流场。

考虑上述问题，本文提出了一种基于压缩感知（CS）从亚采样的PIV数据中重构高时间分辨率流场的方法。该方法无需借助侵入式的高频点测量，使用压缩感知对亚采样的POD系数进行高分辨率重构，可实现12倍的频率放大，且不依赖于高时间分辨率的PIV测量，利用常用的PIV测量装置即可实现。该方法将随机采样得到的PIV数据进行本征正交分解，将所得到的时间系数利用压缩感知进行重构，继而得到高时间分辨率的系数，结合分解得到的空间模态，重构出时间分辨率的流场。本方法首先在较为简单的周期性实验数据中进行验证，并且为探索该方法的适用边界，将该方法进一步应用至非周期的复杂流场（不同直径的双圆柱尾涡流场）中。研究结果表明，该方法可以对周期性流场实现较好的重构，而对于非周期流场则表现较差。

1 流场重构原理

1.1 压缩感知概述

压缩感知理论[9-10]由Candès、陶哲轩及Donoho等提出。该理论指出，只要信号在某个变换域是稀疏的，那么就可以用一个与变换基不相关的观测矩阵将变换所得的高维信号投影到一个低维空间上，然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原信号[11]。由于压缩感知方法基于信号的随机采样，观测矩阵 ϕ定义了在高时间分辨率全流场F上的随机采样规则，基于这样的随机采样规则，亚采样的PIV数据G得以产生，可以定义为：

如式（3）所示，对随机采样的PIV数据G进行本征正交分解，可以获得空间模态 π和相应于数据G的亚采样时间系数Asa。在本文中，由于时间系数在离散余弦变换的变换空间中是稀疏的[12]，因此可将离散余弦变换（DCT）作为所选取的变换域。通过压缩感知，可以从亚采样时间系数Asa中重构出系数A，重构过程求解最优化问题：

其中，ψ为离散余弦变换，S为时间分辨率系数A在离散余弦变换中的系数。但是，求解上述问题是一个NP难问题[13]。Chen、Donoho和Saunders[14]指出，求解一个更加简单的L1优化问题会产生同等的解，使得原问题简化成一个线性规划问题。本文中，选择典型的基追踪算法，使用Matlab中的linprog函数进行求解。结合所得到的空间模态，时间分辨率的流场重建完成。更多细节可见Baraniuk等[13]、Zhao等[15]、Stankovic等[16]对压缩感知理论的论述。

1.2 本征正交分解概述

本征正交分解[17-18]是一种对流场降维，获取流场中主要流动成分的方法，已经得到了广泛应用。记g为瞬态流场数据集G在某时刻的数据，g可以表示为：

其中：表 示平均流场， πi表示空间模态，w表示模态总数，aisa表示相应于某一模态的时间系数。由于在该方法中，数据G为PIV采样数据的集合，因此aisa所记录的时刻点相应于数据G的采样时间。

时间系数asa可通过奇异值分解计算得到：

其中，C表示对流动数据内积 (G,G′)进行空间互相关计算得到的矩阵。将数据G映射到式（4）计算所得的时间系数asa上并进行归一化即可得到空间模态π 。特征值 λ 是逐渐递减的序列，表示模态 π所占有的能量。模态所占有的能量比例通常被记为

1.3 流场重构流程

如图1所示，首先对低时间分辨率随机采样的PIV数据G进行本征正交分解，得到空间模态及相应的时间系数，由于该时间系数与PIV数据G遵循相同的采样规则，因此式（1）中的采样矩阵即为时间系数的观测矩阵，选取离散余弦变换作为稀疏基，从而重构得到高时间分辨率的模态系数，结合前面所得到的空间模态，即可完成流场重构。

图1 流场重构流程Fig. 1 Flowchat of the flow field reconstruction method

2 算例验证与分析

2.1 振荡器流场的验证与分析

振荡器无需任何移动部件即可产生周期性的射流。在流动控制、换热等应用领域中，这种振荡器射流具有良好的鲁棒性[19]，但其振荡频率远高于PIV采样频率。针对这个问题，该时间分辨率流场重构方法首先应用于振荡器实验数据，来验证该方法的可靠性。PIV拍摄区域如图2所示。为验证重构准确性，使用TR-PIV获取真实流场信息，有关该实验的更多详细信息可见Wen等[20]的论文。在本文中，PIV采样数据为间隔采样时间超过0.444 s的随机采样，振荡器频率为2.687 Hz，而重构所得的高时间分辨率全流场数据采样频率为27 Hz。

图2 振荡器PIV拍摄区域示意图Fig. 2 Schematic of the PIV measurement region of the oscillator flow field

为了选择适当数量的模态，对PIV数据的POD模态累积能量谱进行分析，由图3所示，可以发现，在一共7208个模态中，随着模态数的增加，每一模态所占据的能量比例逐渐降低，前7个模态占据了振荡器射流约72%的能量，而前100个模态占据了约88%的能量，因此，在下面的分析中，首先选取前7个模态用于重构，继而选择前100个模态用于对比。

图3 振荡器流场POD模态能量谱Fig. 3 Energy spectrum of the POD modes of the oscillator flow field

图4对比了该模态系数在时间域上的真实值及重构值。低频采样数据G所得的模态系数用黑色菱形标出，可以看出该重构得到的模态系数与真实值有着较好的拟合。尽管相应于其他模态的系数这里未给出，但重构效果基本类似，不再赘述。

基于图4重构得到的模态系数，结合POD分解所得空间模态，完成高时间分辨率流场重构。随机选取某周期中未被采样的三张快照观察其重构效果（图5）。可见，7阶模态重构流场准确获取了真实流场的主模态相关信息，但由于所选模态较少，部分细节信息有所丢失。针对该问题，进一步增加重构阶数至100阶。当使用100阶模态进行重构时，重构流场的细节信息更为丰富。

图4 振荡器流场时间域上重构模态系数与真实值对比Fig. 4 Comparison of reconstructed POD coefficients with ground truth data in the temporal domain for the oscillator flow field

图5 振荡器流场非采样点重构与真实值对比Fig. 5 Comparison of the reconstructed flow field with the ground truth data at a non-sampling point for the flow oscillator

为量化重构准确率，误差由Ruscher等[21]的分析定义为：

其中，uriecon表示在一张PIV快照中空间i位置的重构值，utiruth表示空间i位置的真实参考值，和分别表示真实PIV数据F在时间域和空间域的最大值和最小值。最终每张快照的误差计算值对上式所有所得ei进行平均，继而排除采样点位置处的快照数据，对所得每张快照的误差进行平均。所得误差曲线如图6所示。当选取模态阶数由3阶升至20阶时，重构误差从2.8%降至2.7%，而随着重构阶数进一步增多至100阶时，误差稍微增加。经过分析，认为该误差的增加是由高阶模态包含更多噪声所导致的。

图6 振荡器流场重构误差随选取模态数量的变化Fig. 6 Error incurred in data reconstruction vs the number of POD modes employed of oscillator flow field

2.2 双圆柱流场的验证与分析

为了探索提出的方法在非周期流场中的重构效果，对不同大小的双圆柱尾涡流场原始PIV数据进行了处理分析（图7）。实验在循环水槽[22]中进行，实验的具体细节参见文献[23]。圆柱采样频率为250 Hz，因此本文用于重构的数据间隔采样大于0.048 s。

图7 双圆柱PIV拍摄区域示意图Fig. 7 Schematic of the PIV measurement region of double-cylinder wake flow

图8为 POD模态累积能量谱，可以发现，第一个模态占据了超过80%的能量，前7个模态占据了约92%的能量，因此，在本案例中，仍然选用前7个模态进行重构。

图8 双圆柱流场POD模态能量谱Fig. 8 Energy spectrum of the POD modes of the double-cylinder wake flow

相应于第一模态，图9对比了该模态系数的真实值及重构值，如图9（a）所示，由低频采样数据G所得的模态系数用蓝色五角星标出。在非采样点，重构得到的模态系数偏离于真实参考值，存在较多的高频噪声，而结合图9（b）所示，尽管在非采样点重构值偏离于真实值，但也是在真实值附近波动。这是由于L1范数无法区分稀疏系数尺度的位置，所以存在低尺度能量搬移到了高尺度的现象，从而容易出现一些人工效应。即，如图9（a）所示的一维信号高频振荡现象。整体上来说，重构信号在欧氏距离上仍然逼近原信号。相应于其他模态的系数未在图中示出，但重构效果基本类似，所以不再赘述。

图9 双圆柱流场时间域上重构模态系数与真实值对比Fig. 9 Comparison of the reconstructed POD coefficients with the ground truth data in the temporal domain for the double-cylinder wake flow

随机选取1张在非采样点的重构快照，见图10。可以发现虽然压缩感知方法可以重构流场主要特征，但是难以捕捉高频流场细节。相比于周期性的振荡器射流流场，误差显著增大。这说明，该方法对非周期性流场重构适应性较低。

图10 双圆柱流场非采样点7阶模态重构流场与真实流场对比Fig. 10 Comparison of the reconstrcucted flow field using leading 7 POD modes with the ground truth data at a nonsampling point for the double-cylinder wake flow

3 结 论

研究提出了一种基于压缩感知和POD的PIV流场高时间分辨率重构方法，对低频随机采样数据进行本征正交分解，然后对时间系数进行压缩感知重构，结合分解所得的空间模态，重构出时间分辨率的流场。利用该方法处理了周期性的亚采样振荡器射流实验流场和非周期性的双圆柱实验流场。结果表明，对于周期性流场，该重构方法的重构误差在3%以下，能够准确恢复时间分辨率的流场。而对于非周期性复杂流场，则出现较大的高频噪声。

该方法为一种纯数据驱动的方法，可以实现任意选择的高阶或低阶重构，而且无需任何接触式的传感器，不会扰乱流场，且计算量较低。该方法不依赖于高时间分辨率的PIV测量，利用常用的PIV测量装置即可实现，且可实现12倍的频率放大。因此，该方法对于周期性流场的时间分辨率重构是一种行之有效的方案。但目前该方法对非周期流场的重构适应性较低，需要做进一步研究。