基于随机森林算法的非定常气动力建模研究

徐旺丁， 张 兵， 王华毕

（合肥工业大学 机械工程学院，合肥230009）

1 引 言

非定常气动力建模从理论研究到基于CFD（Computational Fluid Dynamics）技术的非定常流场数值模拟，再到基于CFD技术的非定常气动力降阶模型，是气动弹性力学研究的重要组成部分和关键技术之一。基于CFD技术的气动弹性分析对气体流动的刻画越细致，则气动弹性系统的维数就越高，从而导致气动弹性模拟计算时间耗费巨大，特别是涉及气动弹性优化问题时，有时能达到难以接受的地步。因此，发展非定常气动力降阶模型ROM（Reduced Order Model）已成为一个热门的研究领域。基于气动力降阶模型的气动弹性分析，可以快速求解该气动弹性系统对任意输入的响应，仅需利用CFD求解器求解给定状态下流场的输入输出特性进行构造，无需对程序进行大的改动。

目前，基于系统辨识的模型降阶方法为非定常气动力建模提供了新思路。系统辨识根据系统的输入输出数据来确定描述系统行为的数学模型，其基本原理和依据是现代系统辨识理论，对系统进行分析的主要问题是根据输入信号和系统的特性来确定输出信号。在气动弹性系统中，较为常见的基于系统辨识的降阶方法有 Volterra级数法［1，2］、神经网络模 型［3］、ARMA（Auto-Regressive Moving Average）模 型［4，5］和 支 持 向 量 机 SVM （Support Vector Machines）方法［6］等。非定常气动力 ROM和结构模型可以耦合得到整个气动弹性系统ROM，因此可以方便地用于气动弹性领域的颤振分析［7，8］。

Volterra／ROM在模拟较强非线性的气动弹性现象时会遇到模型维数灾难的问题；而ARMA／ROM是线性系统辨识模型，存在难以对强非线性特性进行建模的理论限制，且需要反复调制输入信号来确定激励特征频率范围；以神经网络等为代表的响应面模型则难以建立传递函数或状态空间模型，且易出现过学习或欠学习的困难，限制了其应用范围［9］；支持向量机在其训练阶段运算量非常大，特别是对于非线性气动弹性分析，需要采用的支持向量机数目非常多，造成计算量急剧增大。

随机森林RF（Random Forest）算法是一种比较新的机器学习算法，其计算开销小，具有较好的大数据集回归预测效果。根据大数据集理论，随机森林不存在过拟合的问题［10，11］，在很多现实任务中展现出强大的性能，是代表集成学习技术水平的方法［12］，该方法在航空航天领域已有较好的应用［13，14］。本文尝试将RF算法应用于建立非定常气动力模型ROM，并将所得到的模型命名为RF／ROM。RF／ROM也是一种基于系统辨识的气动力建模方法，采用CFD方法计算训练信号的非定常气动力，仅以一组训练样本对RF模型进行训练，采用简谐信号对该RF模型进行测试，最后将测试结果与CFD仿真结果对比，检验RF／ROM对非定常气动力的预测精度和效率，并将RF／ROM应用于颤振边界以及极限环振荡（Limit Cycle Oscillation，LCO）特性的预测。

2 随机森林的基本原理

随机森林是由多颗决策树组成的集成学习模型，在处理回归问题时，随机森林以回归树为基学习器。随机森林是在Bagging的基础上引入了随机属性选择的思想，回归效果优于Bagging。根据大数据集理论，随机森林不存在过拟合的问题［10］。

令fi（X）（i＝1，2，…，m）表示m 个子模型，δi表示fi（X）的权值，那么集成模型fE（X）可表示为

随机森林算法如图1所示，先采用Bootstrap方法从训练样本集TTrain＝｛（Xk，Yk）中按照随机的、有放回的和重新选择的方式，产生m个相互独立的样本子集Ti＝ ｛（Xk，Yk）（i＝1，2，…，m）。对于任意的Ti，NB（NB＜N）个样本服从统一分布且相互独立。然后，在每个Ti上都建立一个回归树子模型，并通过随机属性选择的方式分裂根结点或父结点。对于任意样本X，m个子模型将产生m 个预测值P1，P2，…，Pm，最后将m个预测值去平均值后作为随机森林的预测值PE。

从式（2）可以看出，所有子模型的权值相同，δi＝1／m（i＝1，2，…，m），但实际上子模型的样本子集不同，准确度也不同。如果给模型分配相同的权值，则会降低随机森林的预测精度。因此需要对子模型进行融合处理，本文采用GEM（Generalized Ensemble Method）算法［15］融合子模型。该算法是以fE（X）的均方误差最小为目标函数，作为约束条件，通过子模型的偏差构造相关矩阵，从而直接求解权值向量。GEM算法具体过程如下。

图1 随机森林算法Fig.1 Diagram of random forest algorithm

令g（X）代表真实函数，di（X）＝fi（X）－g（X）代表子模型与真实函数的偏差函数，则集成模型还可以表示为

设A为m×m的相关矩阵，且

式中E为单位矩阵。对有限个训练样本近似估计，则

则集成模型的误差可表示为

最优权值δ可通过求解最优化问题式（7）得到，

由拉格朗日定理可知，

式中λ为拉格朗日乘子。获得最优条件：

式中 1v＝［1，1，…，1］T。解得最优解为

3 基于RF／ROM的气动弹性建模

3.1 结构模型

为验证RF／ROM，选择 NACA0012和 NACA64A010翼型模型，模型结构示意图如图2所示。结构动力学系统包含绕弹性轴的俯仰运动α（顺时针为正）和y方向的沉浮运动h（向下为正）两个自由度。结构动力学方程如下。

式中 Sα＝mxαb为对弹性轴的静矩，fy为沉浮方向的外力合力，Mz为对z轴的合力矩。定义无量纲时间τ＝ωαt，二维翼型气动弹性方程可以写成无量纲形式为

式中 无量纲速度v＊＝V／（ωαbμ1／2），V 为来流速度，质量比μ＝m／（πρb2），h＝h／b为结构响应广义位移，rα为翼型围绕弹性轴的回转半径，ωh和ωα分别为弯曲和扭转模态的固有频率，CL和CM分别为升力系数和俯仰运动力矩系数。

3.2 基于CFD技术的ROM训练过程

RF／ROM的训练过程主要是基于CFD求解器的耦合计算和设计来训练输入。对于CFD耦合计算，非定常气动力的计算需要花费大量的计算成本，引入降阶模型用于非定常流场的模拟可以提高气动弹性计算的效率。在气动弹性系统中，结构状态值u和v作为输入值，广义气动力系数向量f作为输出，采样时间序列数据（u，v，f）可以由CFD耦合求计算。如两个自由度的气动弹性系统的输入为（h，α），输出为CL和CM。为了构建图1所示的基于RF的ROM，需要精心准备采样数据作为训练样本集。一个高质量的RF／ROM不仅需要精确的CFD训练计算，还要求训练输入对系统的动态特性进行充分激励。因此，需要精心设计训练输入信号以保证激励具有足够的频带宽度。图3为本文设计的以一种过滤的高斯白噪声信号作为输入训练信号，因为白噪声信号可以很好地代表非线性系统的自然动态特性。该信号是通过采用数字生成的高斯分布随机时间序列来创建，并对其进行滤波处理，可以看出，该训练信号既含有大振幅成分，也含有小振幅成分，频谱也比较宽，适合作为RF／ROM的训练输入信号。

由于非定常流场时间具有延迟效应，所以可以选择xi＝［（ui，vi），（ui－1，vi－1），…，（ui－r，vi－r），fi－1，…，fi－s）作为输入信号，r和s是由用户选择的时间延迟。因此，可以从非定常CFD模拟的结果构建RF／ROM 的新训练样本集（xi，fi）。

图2 二维翼型气动弹性模型Fig.2 Diagram of two-DOF aeroelastic model

3.3 RF／ROM的构建

基于RF的气动弹性ROM的构建流程如图4所示，共有5个主要步骤。（1）设计多组输入信号，并通过非定常CFD求解器计算相应的非线性气动力响应；（2）将输入信号和非定常气动力响应组合成训练样本集（xi，fi）；（3）执行 RF算法获得辨识模型RF／ROM，并用测试集验证模型的精度；（4）将RF／ROM应用于预测不同工况下的非定常气动力；（5）耦合结构动力学系统，根据图4的虚线路径，将RF／ROM用于气动弹性响应的预测。

4 算例验证

为了验证气动力RF／ROM的精确性和高效性，应用3.2节设计的过滤的高斯白噪声信号位移输入对CFD求解器进行激励，选择r＝5和s＝4为RF／ROM的准备训练样本集，根据3.3节流程建立RF／ROM。本文主要做了两方面工作，一是对NACA0012翼型的不同马赫数工况下的颤振边界进行了计算，并将RF／ROM的计算结果和CFD仿真结果以及试验结果进行对比；二是将RF／ROM用于NACA64A010翼型的跨音速极限环颤振特性预测，并与CFD仿真结果对比。

4.1 颤振边界预测

图3 训练输入信号Fig.3 Training input signal

图4 RF／ROM构建流程Fig.4 Workflow of RF-based ROM

颤振边界算例验证选用的NACA0012机翼结构模型、来流Ma、Re以及颤振边界实验数据来自文献［16］，选 择 0.45Ma，0.61Ma，0.71Ma，0.77Ma和0.8Ma共5个工况进行颤振边界分析计算，初始迎角为0°。模型结构参数为m＝87.07kg，c＝0.4064m，Kh＝38821N／m，Kα＝3928N／m，a＝0，xα＝0，Iα＝3.68kg·m2。

在0.71Ma工况下，训练输出的结果如图5所示，可以看出，RF／ROM预测升力系数和力矩系数的精度很高。为进一步验证RF／ROM的正确性，使用简谐信号输入下俯仰运动的气动力系数作为测试信号，图6为在最后一个周期中，ROM预测结果和CFD计算结果的对比。

使用已建立的RF／ROM可以方便地开展高效颤振边界预测，利用变密度的方式调整来流动压，可以快速预测气动弹性系统的颤振边界。如在Ma＝0.71工况下，计算动压q分别设为0.95qe，0.85qe和0.9qe时（qe为颤振动压的实验值），翼型的结构响应如图7所示。图7（a）中广义位移的幅值逐渐增大，说明此时结构振动是发散的，即已经发生颤振；图7（b）中广义位移的幅值逐渐减小，说明此时结构振动是收敛的，尚未达到颤振边界；图7（c）中广义位移表现为等幅振荡，正好位于颤振临界边界，此时的动压为计算颤振动压，相应的振动频率为计算颤振频率。

图5 NACA0012翼型CFD直接模拟和RF／ROM预测的输出对比Fig.5 Results of NACA0012airfoil comparison between CFD and RF／ROM

对5个工况下的翼型颤振边界进行计算，计算结果和CFD直接模拟计算以及文献［16］的试验值进行对比，如图8所示。可以看出，RF／ROM的预测结果和CFD直接模拟计算的结果吻合良好，表明RF／ROM预测颤振边界特性具有较高的精度。然而，模拟值和实验值匹配不佳的主要原因是部分结构参数，如重心位置等对颤振结果存在显著的影响，且实测数据存在着一定的误差，因此颤振实验值的离散度较大［17］。

在上述算例中，其中一种马赫数工况下，CFD直接模拟颤振特性通常需要进行4～5次的CFD／CSD耦合计算方可得到颤振临界点，每次CFD／CSD耦合计算需要计算5～8个周期，约1250～2000个时间步，共需5000～10000个时间步。而对于RF／ROM方法，在某种马赫数工况下，只需要进行一次训练过程，约500个时间步，除此之外的RF／ROM计算时间可以忽略不计。显然，在定马赫数工况下，基于RF／ROM的颤振计算效率提高了约10～20倍。气动网格数目越多，则CFD计算量越大，RF／ROM方法的效率越明显。

图6 最后一个周期气动力系数对比Fig.6 Comparison of aerodynamic coefficients for the last cycle

4.2 LCO预测

由于非线性因素，弹性结构在流场中的运动形式通常表现为极限环振荡。针对因跨音速气动力非线性而产生的极限环型颤振分析研究［18］，选择NACA64A010翼型作为结构模型，结构参数为xα＝0.25，r2α＝0.75，a＝－0.6，μ＝75，ωh／ωα＝0.5。先运用CFD方法计算该模型的LCO特性，然后建立RF／ROM，并开展LCO特性计算，将两种方法的计算结果进行对比，验证RF／ROM预测LCO的正确性和计算精度。

在 Ma＝0.8，v＊＝0.7工况下，给定初始沉浮运动和俯仰运动速度值均为0.1，CFD耦合计算LCO特性和RF／ROM预测LCO特性的结果对比如图9所示。

为了进一步验证RF／ROM预测LCO特性的性能，计算了马赫数为0.76，0.8，0.825，0.85和0.875共5种情况下的气动弹性系统的LCO特性，所有马赫数情况下均取v＊＝0.75。运用CFD方法和RF／ROM方法计算的LCO振幅如图10所示，并引入参考文献［19］的CFD和SVM／ROM计算结果作为参考和对比。可以看出，RF／ROM和CFD计算的结果很一致，平均误差在1 0%以内。CFD方法捕获到LCO振幅需要大约1小时，而RF／ROM方法预测相同的LCO特性不超过5分钟。故RF／ROM的计算效率较高，这对实时飞行模拟和飞行控制器的设计有重要意义。

图7 0.71Ma，不同动压下的NACA0012翼型结构响应Fig.7 0.71Ma，NACA0012airfoil structure response under different dynamic pressure

图8 NACA0012翼型的颤振边界计算结果Fig.8 Flutter boundary of the NACA0012airfoil

图9 NACA64A010翼型在 Ma＝0.8，v＊＝0.7工况下的LCO响应对比Fig.9 Comparison of LCO response at Ma＝0.8，v＊＝0.7

图10 NACA64A010翼型在不同马赫数工况下LCO振幅对比Fig.10 Comparison of LCO amplitude with Mach number

5 结 语

基于CFD技术，将随机森林算法引入非定常气动力降阶模型建模研究领域，详细说明了RF／ROM的构建过程，采用所构建的RF／ROM实现对二维翼型模型的颤振边界和LCO特性的预测，并与CFD直接耦合计算的结果进行对比。结果表明，RF／ROM可以用于非定常流场中的气动弹性特性预测，计算效率也提高了10～20倍。如何进一步提高ROM的精度和效率，特别是对具有变参数的模型，将会成为未来气动力建模研究领域的重点。