探讨流场数据填补方法与GappyPOD方法有机融合

　　摘要：本征正交分解方法是一种数据驱动的流场特征信息提取技术，可以按能量大小给出流场的结构模态，并且可以通过较少阶模态叠加获得高阶数据的近似描述。将该方法结合一组线性方程，即构成GappyPOD方法，可实现对缺失流场的重构。论文从样本数据已知和样本数据不完全已知两种情况展开，对GappyPOD方法在缺失流场数据填补方面的应用进行了研究。首先，对于样本数据已知的情况，研究了样本参数范围内任意参数值下的重构能力；其次，对于样本数据不完全已知的情况，研究了模态阶数、缺失率以及样本数对缺失流场重构精度的影响规律。结果表明，GappyPOD方法可以高效再现参数范围内的任意完整流场数据。但是，对于缺失率较高的样本集，需要适当增加样本数以提高重构精度。
　　关键词：GappyPODPOD数据驱动流体力学流场数据填补方法流场重构非定常流动
　　随着计算机技术和流场测量技术的高速发展，研究者获得流场数据越来越便利，也越来越精细〔1，2〕。与此同时，面对海量数据，高效地处理分析数据，获取有用信息，以支撑未来飞行器设计发展，逐渐成为了研究热点。目前，基于数据驱动的数据挖掘和机器学习等技术成为了处理和分析数据的主要途径。基于特征提取的本征正交分解方法作为一种高效的降阶方法，已在流场数据信息提取〔3，4〕和非定常流场动力学分析〔5，6，7〕等方面得到了广泛应用。
　　POD方法〔8〕的主要思路是寻找一组最佳的标准正交基，使得样本数据在该标准正交基上的投影依次迅速递减，截取投影较大或者说包含能量较高的前几阶模态，从而可以用较少的基展开获得较高阶数据的近似描述。该方法已成功应用于诸多领域，如图像处理，信号分析和数据压缩等〔9〕。在流体力学领域，POD方法主要有两方面的应用。一是用于特征模态的提取。Lumley〔10〕首先将POD方法引入了湍流领域，用于提取湍流的主要特性或拟序结构。魏斌斌等〔11〕将POD方法引入高超声速脉动压力数据处理中，发展了基于POD方法的高超声速边界层转捩判定方法。李卓越等〔12〕将POD方法应用于实验测量数据的分解，分析了受激后向台阶分离流的流场特性。二是用于构建气动力降阶模型。Deane等〔13〕结合Galerkin投影的方法，成功使用基于POD方法的降阶模型预测了槽道流和圆柱绕流的极限环行为。基于POD方法的降阶模型与结构模型进行耦合，可以实现气动弹性的快速预测〔14〕，如颤振特性分析〔15，16，17〕、极限环预测〔18，19〕以及气动弹性主动控制律设计〔20，21〕等。在保证计算精度的前提下，极大地提高了计算效率。另外，POD方法在气动外形优化〔22〕、气动外形反设计〔23〕和流动控制〔24〕等方面也得到了应用。
　　Everson等〔25〕针对不完整样本或缺失流场，在基本POD方法的基础上，发展了GappyPOD方法。给定一组POD模态，一个不完整流场快照可以通过求解一组线性方程组得以重构。并且，如果用于求解POD模态的样本存在数据缺失，那么可以通过迭代的方法实现缺失流场重构。该方法已成功应用于不完整人脸图像的重构。对于流体力学问题，几种情况可能导致数据的不完整。首先，实验中可能只得到了某些点的测量值，GappyPOD方法就可以通过这些已知点重构出完整流场的数据〔26〕。其次，流场的某些数据是未知的。例如，已知一组翼型形状和其相对应的流场解，给定一个新的翼型形状，这时，通过GappyPOD方法可以快速估算出新翼型对应的流场解〔27〕。另外，GappyPOD方法还可以应用于翼型反设计问题〔28，29〕。本文在上述应用的基础上针对GappyPOD方法在数据填补应用中可能会遇到的数据集缺失情况以及重构过程中的参数选择问题开展研究，以期为GappyPOD方法在缺失流场数据填补方面的应用提供较为全面的参考经验。
　　本文以非定常圆柱绕流流场作为样本数据，针对样本数据已知和不完全已知两种情况，探究了模态阶数、缺失率以及样本数对缺失流场重构精度的影响规律。结果表明，两种情况下，GappyPOD方法均能以较少阶模态高精度再现完整流场。另外，在样本缺失率较高的情况下，可适当增加样本数以提高重构精度。
　　1、GappyPOD方法
　　1。1适用于缺失流场重构的GappyPOD方法
　　GappyPOD方法是POD方法的一种变形，用于填补流场中的缺失数据。在获得某个参数范围内的一系列采样解后，采用GappyPOD方法可以填补该参数范围内任意参数值下的缺失数据。
　　首先建立一个标签向量，用以描述所需重构的流场中每个位置上的数据是已知还是缺失。对于流场解Uk，其相应的标签向量nk定义如下：如果Uki（Uki代表流场解向量Uk中的第i个元素）缺失或不正确，则nki0；如果Uki已知，则nki1。为方便起见，将解向量Uk中缺失或不正确的元素通过点乘（nk，Uk）inkiUki置为零。内积定义为（u，v）n〔（n，u），（n，v）〕，相应的模定义为（vn）2（v，v）n。
　　设｛i｝mi1为一组样本集｛Ui｝mi1的POD基，样本集中的所有快照均完整已知。设g为样本参数范围内具有一部分缺失数据的流场解向量，其对应的标签向量为n。假设g为需要重构或填补的解向量，且该解向量的流场特征可以由上述样本集的POD基表示，则由前p阶POD基表示的重构解向量可表示为
　　gi1pbii（1）
　　定义原始流场解向量与重构解向量的误差为
　　Egg2n（2）
　　通过保证该误差最小以求得模态系数。为保证E最小，可以将式（2）对各阶模态系数bi求导，从而得到如下线性方程组
　　Mbf（3）
　　式中，Mij（i，j）n，fi（g，i）n。求解式（3）得到b，并带入式（1），得到重构流场解向量g。最后，将原始解向量g中的缺失点用重构解向量g中的相应数据代替，即可得到完整的流场解，即，当ni0时，gigi。
　　1。2适用于不完整样本集的GappyPOD方法
　　上述GappyPOD方法可以延伸到样本集不完全已知的情况，在这种情况下，POD基可以通过迭代的方法得到，该方法由Everson和Sirovich提出。设一组不完整的样本集｛gk｝mk1，相对应的标签矩阵为｛nk｝mk1。首先将样本中的缺失元素用其余样本在该点的平均值代替，如下
　　hki（0）｛gki，gi，nki1nki0
　　式中
　　gi1Pik1mgki，Pik1mnki
　　并且，令hk（l）为原始流场解向量hk第l次内迭代的猜想解向量，由这些猜想解向量构成的样本集可以得到一组POD基，这组POD基将用于下一步迭代中不完整快照的更新。由l0开始，迭代步骤总结如下
　　（1）对样本集｛hk（l）｝mk1进行基本POD模态分析，得到当前迭代步的POD基向量｛k（l）｝mk1。
　　（2）如适用于缺失流场重构的GappyPOD方法所述，采用前pn阶POD基向量重构，以填补样本中的缺失元素。当前迭代步的临时重构解向量为
　　hk（l）i1pbki（l）i（l）
　　（3）该临时重构解向量用于填补原始解向量中的缺失元素，构成下一迭代步的临时猜想解向量，如下
　　hki（l1）｛hki（l），hki（l），nki1nki0
　　（4）令ll1，进入下一步迭代。
　　循环上述迭代过程直到达到最大迭代步骤或达到收敛要求，该迭代过程称为内迭代。令pn1，重复上述内迭代过程，进入下一循环的主迭代。Gunes等〔30〕提出，可以由直接投影得到的模态系数代替上述方法中的模态系数，这样省去了最小二乘的计算过程，提高计算效率。Venturi等〔31〕提出第n次主迭代中的POD模态阶数均可取为pn，且给出了外迭代的收敛机制为：当KpnKp0达到最大值时，外迭代停止，其中，
　　Kpni1N2n，j，i
　　内迭代的收敛机制为
　　1pni1pnn，j，in，j1，in，j，i
　　一般取0。1即可达到较高的收敛精度。
　　2、算例验证
　　以Re100的二维圆柱绕流为研究算例，选择中心格式有限体积法数值求解NS方程，得到样本数据。选取足够大的计算域，从而减小边界条件对流场计算结果的影响。本文的计算域如图1所示，入口边界距离圆心20D，出口边界距离圆心40D。将计算得到的升力系数幅值和阻力系数的均值与文献已有计算结果〔32，33，34〕进行了对比，见表1。结果表明，本文采用的数值方法有很高的精度，是可靠的。
　　图1计算网格
　　表1Re100圆柱绕流的升力系数幅值和阻力系数均值对比
　　2。1适用于缺失流场重构的GappyPOD方法
　　为了验证适用于缺失流场重构的GappyPOD方法的可行性，本文采用上述圆柱绕流构造缺失流场算例。CFD计算中一个周期保留120个样本点，文中每间隔2个点取1个点，即一共有40个采样点，该40个采样点用来构造POD基。选取该120个样本点中的第20个样本点作为缺失流场，该点不在40个采样点内，假设缺失流场中只有圆柱表面的200个点的数据已知。文中以横向速度场为例展示，其余变量具有相同的规律。
　　首先，对40个采样点样本集进行POD分析，得到用于重构缺失流场的POD基。如图2和3所示，分别为特征值曲线和前6阶POD模态。从图中可以看出，特征值是成对出现的，且各对间近似为指数递减的趋势。模态图关于x轴成对称或反对称，前两阶模态反映了流场的大尺度特征。
　　图2特征值曲线
　　图3前6阶POD模态
　　对缺失流场的横向速度场进行填补，图4为使用不同阶POD模态的缺失数据填补结果与CFD计算结果对比图。可以看出，只使用前2阶模态的填补结果再现了流场的大尺度特征，缺失流动的小尺度特征。使用前4阶POD模态，小尺度特征得以体现。当模态阶数增加到10阶时，GappyPOD的填补结果几乎与CFD计算结果完全吻合。为了更直观地显示填补结果与CFD计算结果之间的误差，文中定义结果误差为
　　error1Ni1NgGappyPODigCFDi
　　式中，N表示流场中缺失元素的个数。图5为填补的结果误差随所使用的POD模态数的变化曲线，可以看出，随着使用POD模态阶数的增加，填补结果逐渐趋于收敛。
　　图4GappyPOD填补结果与CFD计算结果对比
　　图4GappyPOD填补结果与CFD计算结果对比
　　图5误差随重构所用POD模态数的变化曲线
　　2。2适用于不完整样本集的GappyPOD方法
　　为了构造有数据缺失的流场，文中在CFD计算得到的完整流场中任意舍弃某些离散点上的数据。为了量化流场中数据的缺失程度，文中定义未知的离散点个数占总离散点个数的百分比为缺失率。采用随机抽样的方式从流场离散点中抽取相应比例的点数，舍弃其真实值，作为缺失点。如果取所用样本集的缺失率为某一值，则样本集中各快照的缺失率在该值附近有微小变化。文中采用各快照的缺失率分布的均值为，方差为的正态分布，其中的值等于样本集的缺失率。
　　仍选择上述Re100的二维圆柱绕流作为研究算例。下文首先以样本个数40，空间缺失率20的算例进行方法的结果展示。各快照缺失率正态分布的均值为0。2，方差为0。01。为了更直观地显示GappyPOD方法的填补结果与CFD计算结果的误差，文中定义误差为
　　error21mj1m1Nji1NjgGappyPODigCFDi
　　式中，Nj代表第j张快照中缺失元素点的个数。图6为该误差随重构所用POD模态阶数的变化曲线，可以看出，随着模态阶数的增加，误差单调递减，且最终收敛于接近0的值。使用10阶模态，即外迭代步达到10步，即可得到较精确的结果。
　　图6误差随重构所用POD模态阶数的变化曲线
　　图7为各外迭代步得到的特征值曲线与原始完整样本的特征值曲线对比图。可以看出，随着迭代步数的增加，迭代步的特征值曲线逐渐收敛到原始样本的特征值曲线，且低阶模态的特征值几乎与真实特征值完全吻合。但是随着迭代步数的增加，收敛速度逐渐减慢，且高阶模态的特征值与真实特征值相差较大。也就是说，当迭代步刚开始时，GappyPOD方法首先捕捉到流场的大尺度特征，随着迭代步数的增加，逐渐捕捉到流场的细小特征，对流场的重构也更精细。但是样本中离散点数据的缺失，导致某些流场特征的丢失，不能全完复现原始流场的特征。
　　图7各迭代步的特征值曲线与原始样本特征值曲线对比图
　　为了测试该方法在不同缺失率以及不同样本数情况下的重构能力，文中考察了缺失率固定为0。2，0。5和0。7，而样本数分别取40，80和120的情况，以及样本数固定为40，80和120，而缺失率分别取0。2，0。5和0。7情况下的重构误差随所用POD模态阶数的变化，如图8所示。可以看出，随着样本缺失率的增加，不同样本数下的重构误差随POD模态阶数变化的收敛速度随之降低，且能够达到的收敛精度也随之降低。当缺失率较低时，不同样本数下的重构精度相差不大；当缺失率较大时，增加样本数可以提高重构精度。
　　图8不同缺失率不同样本数的误差随重构所用POD模态阶数的变化曲线
　　图9给出的是缺失率为0。7时，分别选取样本数为40，80和120，缺失流场的填补结果与CFD的计算结果对比图，图中所示结果的计算误差均已收敛。可以看出，当缺失率为0。7时，只选取40张样本得到的重构流场较为模糊，而当样本数增加到80或者120时，重构流场的清晰度明显提高，且与CFD结果的吻合度较高。
　　图9缺失率为0。7的GappyPOD填补结果与CFD计算结果对比
　　3、结论
　　本节对适用于缺失流场重构的GappyPOD方法和适用于不完整样本集的GappyPOD方法进行了研究，以Re100的二维圆柱绕流作为研究算例，对上述两种方法的可行性进行了探究。得出如下结论：
　　（1）对于已知一组POD模态，只有某一流场存在部分数据缺失的情况，适用于缺失流场重构的GappyPOD方法可以高效再现完整流场数据。随着重构所用POD模态阶数的增加，重构误差快速收敛，并且收敛于较高的重构精度。
　　（2）对于用于求解POD模态的样本集各快照均存在部分数据缺失的情况，适用于不完整样本集的GappyPOD方法可以高效再现整个样本集数据。随着重构所用POD模态阶数的增加，重构误差快速收敛，并且收敛于较高的重构精度。
　　（3）对于用于求解POD模态的样本集各快照均存在部分数据缺失的情况，当数据缺失率较低时，使用较少的样本数即可达到较高的重构精度；当数据缺失率较高时，使用的样本数较低会导致重构精度不高，可增加样本数以提高重构精度。