机器学习数据融合方法在火箭子级栅格舵气动特性建模应用中的比较研究

2022-07-14许晨舟杜涛韩忠华昝博文牟宇张津泽

实验流体力学 2022年3期

关键词：栅格力矩气动

许晨舟，杜涛，韩忠华,，昝博文，牟宇，张津泽

1.西北工业大学航空学院，西安 710072

2.翼型、叶栅空气动力学国家级重点实验室，西安 710072

3.北京宇航系统工程研究所，北京 100076

0 引言

飞行器气动设计的核心工作之一是提供准确、完备的气动特性数据，其中最为重要的是冻结外形气动特性数据，国外一般称为“设计气动数据库”。然而，对于流动条件变化复杂的高超声速飞行器，建立设计气动数据库通常需要开展规模庞大和复杂的风洞试验和CFD 计算，导致气动研制经费和进度难以控制。因此工程应用迫切需要发展新的方法降低气动数据库建立成本。

近年来，随着数字化飞行器设计技术的发展，研究人员在人工智能技术与空气动力学结合方面开展了诸多有益的尝试与探索。其中具有代表性的机器学习方法已经成为最主流的人工智能实现方式，被广泛应用于飞行器多源气动数据融合与建模。机器学习方法可以快速预测不同工况下飞行器的气动特性，对建立完备的气动数据库、降低飞行器设计成本以及提升飞行器综合性能具有重要意义。

变可信度融合模型是机器学习方法中的代表性方法。它是一种针对小样本数据的监督式机器学习方法，其核心思想是先利用大量的低可信度（成本低、精度差）数据去近似气动特性的变化趋势，再融合高可信度（成本高、精度好）数据进行修正，实现对未知状态下气动特性的高精度预测。目前常用的变可信度融合模型大致可以分为以下4 类：

1）基于标度函数的修正模型

基于标度函数的修正模型是在低可信度模型基础上，以相加、相乘或两者混合的方式引入标度函数，从而构建高可信度模型的近似模型。Haftka和Chang 等最早提出了乘法标度函数的修正方法，通过缩放比例系数可使低可信度模型在局部近似于高可信度分析的结果。Alexandrov 等将其与置信域法结合，提出了一阶加法和乘法标度函数，并应用于翼型和机翼的优化设计之中，显著降低了计算成本。舒乐时等和Park 等建立了更符合全局模型变化趋势的差异函数，提高了变可信度模型的预测精度。Wang提出了一种基于最小二乘法的加法标度函数构造方法，提高了模型的泛化能力。杜涛等提出了一种针对变可信度模型的自适应预处理技术，解决了加法标度函数修正模型在原始数据绝对值较小情况下出现数据失真的问题。韩忠华等改善了混合标度函数修正方法，并引入了梯度等信息，进一步提高了模型的预测精度。

2）空间映射模型

该方法最早由Bandler 在1994年提出，Bandler建立了高、低设计空间的映射关系，使低可信度函数的局部特性能够逼近高可信度函数。而后Robinson等提出了一种改进的空间映射方法，并应用于扑翼的气动建模与优化中。Jonsson 等也将此方法应用于跨声速机翼的气动优化设计中。空间映射模型最主要的特点是允许高、低可信度模型的设计空间维数不一致，并且适用于高可信度样本较少的高维建模问题。

3）Co-Kriging 类模型

Co-Kriging 是在Kriging 模型基础上发展起来的一类针对多可信度数据的高斯回归模型，被认为是一种浅层机器学习方法。其基于贝叶斯理论，通过建立自回归模型将多种精度的数据进行融合，利用交叉协方差来衡量不同精度数据间的相关性。Kennedy 和O’Hagan首次将其应用于工程科学领域，提出了经典的KOH 模型。之后，该模型的算法理论被不断改进，包括采用多项式缩放系数、引入更多不同可信度层数据、建模样本非嵌套化等。韩忠华和 Zimmermann 等提出了一种能够避免复杂交叉协方差计算的Co-Kriging 模型。此外，韩忠华等还提出了一种通用性更强、建模更加简洁的分层Kriging 模型（Hierarchical kriging，HK）。该模型的基本思想是以大量的低可信度数据建立全局趋势模型，通过少量的高可信度样本进行修正，最终建立起高精度的预测模型。HK 方法的后续改进研究可以参阅文献[6]和[8]。

4）多可信度神经网络模型

多可信度神经网络模型（Multi-Fidelity Neural Network，MFNN）采用最新发展的深度学习技术，先根据低可信度样本建立神经网络模型，再将低可信度模型的预测值和高可信度样本同时作为高可信度神经网络的输入，最后根据输出的预测值与高可信度样本标签值之间的差异来构建损失函数，迭代训练模型参数。Meng 等在MFNN 模型的基础上提出了多可信度贝叶斯神经网络，该方法可提供模型预测值的均方误差。何磊等将MFNN 网络应用于亚声速翼型的气动特性建模。张鑫帅等则将MFNN 网络应用于翼型与翼身组合体的气动优化设计。

虽然变可信度融合模型理论与算法在近十年已取得了长足的进步，但现有方法的泛化能力仍明显不足，大多数模型还停留在实验室研究阶段，缺乏工程研制中复杂应用场景的适应性研究。

2019年，在我国首次开展的在役火箭子级落区控制项目中，由于飞行操态较多，跨越马赫数范围广，气动研制工作面临经费和进度的双重压力。本文作者所在团队首次尝试引入变可信度融合模型技术，在解决了4 项关键算法和技术问题、提升了模型在复杂工程问题中的适应性，利用加法标度函数修正模型和HK 模型学习了项目部分风洞试验及CFD计算结果后，预测了全部工况的气动特性，按时建立了飞行气动数据库，并大幅降低了风洞试验和数值仿真的需求。2019年7月，飞行搭载试验取得成功，预测结果的正确性获得了验证。

本文是上述研究工作的继续。在应用中作者发现相关函数和模型的选择对预测结果的正确性和质量有较大影响，因此，本文将深入研究影响机理，进一步发掘变可信度融合模型在工程复杂场景中应用的鲁棒性和适应性。除文献[2]中的加法标度函数修正模型和HK 模型外，本文还将额外对比2 种变可信度融合模型—Co-Kriging 模型和MFNN 模型。这4 种变可信度融合模型几乎代表了当前主流的气动融合建模思路。通过对建模结果的对比分析，讨论不同相关函数和模型的工程适应性，并给出变可信度融合模型的使用建议。

1 变可信度模型理论与算法

1.1 基本假设

对于一个维的建模问题，假设存在通过不同来源或不同精度的分析模型得到的高、低可信度样本数据集：

以及对应的响应值：

其中，下标“1”表示高可信度，“2”表示低可信度；和分别为高、低可信度样本点个数， R为实数集。

变可信度模型构建中所需要解决的核心问题是如何有效地融合高、低精度分析模型的数据，使预测结果能更逼近真实的高精度结果。下面将介绍4 种常见的变可信度融合模型理论。

1.2 加法标度函数修正模型

基于加法标度函数修正的变可信度融合模型（Additive Scaling Function Based Multi-fidelity Surrogate Model，AS-MFS）是一种以低可信度模型为基础、通过加法标度函数的方式引入高可信度数据辅助构建高可信度分析模型的近似模型，具体的公式推导可参见文献[8]。在AS-MFS 中，标度函数（）定义为高可信度样本点处高、低可信度分析模型输出响应间的差值：

由于概念直观且模型简洁，AS-MFS 已成为标度函数修正模型中最常用的方法，被广泛应用于工程问题的建模和设计过程中。

1.3 Co-Kriging 模型

本文选用的Co-Kriging 模型来自文献[29]。假设模型预估值为已知样本响应值的线性加权：

为了降低模型的复杂程度，文献[29]假设高、低可信度数据具有相同的过程方差。同时为了描述不同可信度数据之间的相关性，额外引入了一个模型超参数[0，1]，于是通过拉格朗日乘子法可将公式（6）转为矩阵形式：

经过推导可得，Co-Kriging 模型的预估值表达式为：

1.4 分层Kriging 模型

韩忠华等于2012年提出了分层Kriging 模型，它是对传统Kriging 模型理论的进一步发展。首先，针对低可信度数据集建立Kriging 模型：

其次，假设在高可信度样本点处存在静态随机过程：

其中，Z（x）表示均值为0、差为的静态随机过程。这里低可信度Kriging 模型的预估值ˆ（）与常数缩放因子的乘积代替了原先的常数回归项，成为高可信度模型的全局趋势模型。与Kriging 模型类似，可推得HK 模型的均方误差为：

并满足如下无偏估计的条件：

最终，HK 模型的预估值可表示为：

1.5 多可信度神经网络模型

除了基于贝叶斯理论的变可信度模型，近年来新发展了一种基于深度学习的多可信度神经网络模型（MFNN），其模型结构如图1所示，具体推导过程可以参阅文献[31]。

图1 多可信度神经网络模型结构示意图Fig.1 Schematic of multi-fidelity neural network

假设高、低可信度数据之间的关系可以表示为：

其中，（）表示数据从低可信度层到高可信度层的映射关系。一般认为在高、低可信度数据层均存在线性和非线性相关性，于是映射函数（）可拆为两部分：

为了描述高、低可信度数据间线性和非线性相关性的程度，引入超参数缩放因子，于是MFNN 模型的预估值表达式为：

2 在栅格舵气动建模与预测上的应用和比较

2.1 模型训练样本集

预测对象为在役火箭的一子级，如图2所示，其级间段设计了4 片对称分布的栅格舵，对再入过程的俯仰、偏航和滚转通道进行气动操纵。本文将以俯仰通道气动特性为研究对象，输入变量为马赫数、攻角（）和俯仰舵偏角。

图2 火箭一子级落区控制的栅格舵方案示意图[2]Fig.2 Illustration of grid fins in the rocket first stage landing area control project[2]

项目以风洞试验（图3）为主，辅以少量CFD 计算（图4）。由于进度和经费的原因，风洞试验和CFD 计算的结果并不充分。亚声速和跨声速条件的试验成本低、进度快，均开展了试验；马赫数大于4 的状态未开展试验，最终仅开展了所需工况的54%，如表1所示。理论上，风洞试验状态比CFD 计算更加接近飞行器表面真实的流动状态，因此，本文将风洞试验数据视为高可信度数据，而CFD 计算结果是低可信度辅助数据。

图3 火箭一子级带栅格舵模型在风洞中试验情况[2]Fig.3 Experimental model of rocket first stage with grid fins in the wind tunnel test[2]

图4 火箭一子级带栅格舵外形数值仿真的计算网格[2]Fig.4 Computational grids of CFD simulation for rocket first stage with grid fins[2]

表1 俯仰特性风洞试验情况[2]Table 1 Fundamental state for wind tunnel experiment on pitching moment characteristics[2]

亚声速段时，火箭一子级带栅格舵外形的俯仰力矩在0°攻角附近为静稳定配平点，而超声速段变为不稳定配平点（图5），跨声速段容易出现显著的壅塞现象。栅格舵的俯仰力矩表现出远比轴向力和法向力复杂的变化特性。这种特性能够代表工程应用的复杂场景，因此本文将俯仰力矩特性的气动建模与预测作为主要研究对象。

图5 火箭一子级带栅格舵外形跨声速段和超声速段俯仰力矩特性Fig.5 Pitching moment characteristics of rocket first stage with grid fins in the transonic and supersonic regimes

从气动数据集中筛选出不同舵偏角、马赫数、攻角（输入数据）与俯仰力矩系数（输出数据），得到不同舵偏角下的310 组风洞试验数据和96 组CFD 计算数据，其空间分布情况如图6所示，4 个舵偏角从上往下依次为–20°、–10°、0°和10°。从图中可以看出：高可信度风洞试验数据（绿色球体）主要集中在跨声速段，缺少高超声速段数据（马赫数4 以上）；而CFD 计算数据（红色方块）覆盖了全马赫数工况，但缺少大攻角下（攻角15°以上）俯仰力矩数据。

图6 栅格舵俯仰力矩系数数据集空间分布情况Fig.6 Pitching moment data set of rocket first stage with grid fins

2.2 模型训练的关键问题

本文采用4 种变可信度融合模型对俯仰力矩系数进行建模与预测。下面先对模型训练过程中出现的一些关键问题进行分析和讨论。

除了MFNN 模型，其余3 种变可信度融合模型均基于高斯过程回归理论。高斯过程回归—本文特指Kriging 模型—是一种使用高斯过程先验对数据进行插值分析的模型。它通常包含不超过两层的非线性特征变换，因而被认为是一种浅层的机器学习模型。文献[2]中，我们发现相关函数R 和模型超参数训练对Kriging 模型的预测效果影响很大。相关函数又称为“核函数”，是一类仅与样本点空间距离相关的光滑函数，满足距离为0 时等于1、距离无穷大时等于0，且相关性随着距离的增大而减小。对于一个维的建模问题，相关函数一般具有如下的形式：

其中，、x为设计空间中任意两个不同位置，为待定的模型超参数。一般来说，带有的相关函数可以更好地反映不同变量在该维空间上的变化尺度，通过训练可以提高Kriging 模型对非线性函数的拟合能力。另外，还有一类不带的相关函数，如薄板样条函数，它们往往对非线性函数的拟合较差，而对线性函数的拟合较好。

目前常用的带的相关函数包括三次样条函数和高斯指数函数。本节将着重考察这两种相关函数对预测效果的影响。

三次样条函数的解析表达式为：

其中，ξ=θ|x-x|，=1，2，，。该函数二阶可导，在光滑性和鲁棒性方面能取得很好的折中。

高斯指数函数的解析表达式为：

其中，p为相关函数光滑程度的各向异性参数。p=2时，相关函数无穷阶次可导；p ＜2时，相关函数一阶可导。

图7、8 分别为超参数取不同值时三次样条函数和高斯指数函数的变化情况。从图中可以看到，随着的增大，相关函数曲线逐渐变“陡峭”，这意味着较大的会缩小作用范围。取值相同时，对比三次样条函数和高斯指数函数（p = 2），可以发现：三次样条函数随样本空间距离的变化更剧烈，并且当样本空间距离大到一定程度时，二者的相关性消失。文献[35]将这类样本局部性影响较大而全局性影响较小的相关函数称为“紧支撑性函数”。图9展示了= 1 时高斯指数函数取随参数p 的变化情况。可以看到，随着参数p 的增大，相关函数在样本距离较小的情况下过渡得更为平缓。为了保持这种光滑的过渡性，本文限制高斯指数函数参数p 的取值范围为[1.5，2]，且参数p 与超参数都可以通过最大似然估计的方法优化得到。

图7 不同超参数下的三次样条函数变化Fig.7 Influence of θ on the cubic correlation function

图8 不同超参数下的高斯指数函数变化（p = 2）Fig.8 Influence of θ on the Gauss correlation function with p = 2

图9 不同参数p 时的高斯指数函数变化（θ = 1）Fig.9 Influence of index p on the Gauss correlation function with θ = 1

上文通过理论分析解释了两种相关函数的差异性，下面将其应用于栅格舵的俯仰力矩参数建模问题，考察不同相关函数对预测效果的影响情况。本文以舵偏角=0°时的俯仰力矩系数建模为例。需要说明的是，当俯仰舵无偏转时，共有100 组风洞试验数据和24 组CFD 计算数据；模型的第一维输入是马赫数，第二维输入是攻角。

为在同一标准下进行对比，首先将两种相关函数的超参数优化范围均限制在[0.001，1.0]。图10 为HK 模型采用不同相关函数的预测结果对比。从图中可以看到，采用紧支撑性较强的三次样条函数时，预测曲面变得不光滑，在离样本点较远的区域出现波浪型“隆起”现象；而采用高斯指数函数时，整个预测曲面在设计空间内的变化都较为光滑。

图10 基于不同相关函数的HK 模型对俯仰力矩系数建模结果对比（超参数优化范围受限制）Fig.10 Comparison of predicted pitching moment coefficients via HK model based on different correlation functions （with optimization boundaries for hyperparameters）

然后，解除对超参数优化范围的限制，基于最大似然估计理论对超参数进行彻底优化，得到优化后的超参数及俯仰力矩系数预测结果，如图11 所示。从图中可以看到：采用三次样条作为相关函数时，优化后超参数的第一维数值较大，反映在预测曲面上，俯仰力矩系数预测值沿马赫数方向出现了过拟合现象；而采用高斯指数函数时，预测模型未出现明显的过拟合现象。造成这一现象的原因是：三次样条函数的紧支撑性使得其在较大时样本的影响域很小（图7），一旦预测点稍微偏离样本点，就无法获取样本响应值的相关信息，预测值将迅速靠近模型自身的常数回归项。在HK 模型的外插区域，由于缺少高可信度风洞数据，采用高斯指数函数的模型预测曲面出现了小范围凹坑与鼓包。根据马赫数无关性原理，该现象不符合实际的物理规律。这也从侧面反映出：并不是模型的超参数优化得越彻底，预测效果就会越好；而是需要结合建模问题本身，给出合理的优化空间。

图11 基于不同相关函数的HK 模型对俯仰力矩系数建模结果对比（超参数优化范围不作限制）Fig.11 Comparison of predicted pitching moment coefficients via HK model based on different correlation functions （with no optimization boundaries for hyperparameters）

本文将超参数的优化范围限制在[0.001，1.0]，再对比不同变可信度高斯回归模型采用这两种相关函数的俯仰力矩系数预测曲线，截取了攻角为10°的结果，如图12 所示。从图中可以看到：采用三次样条作为相关函数时，3 种变可信度模型的预测曲线均出现了不同程度的“振荡”现象；而采用高斯指数函数时，模型预测曲线变化较为光滑，更符合物理规律。图12 充分显示了相关函数对预测结果的重要性。在后续的分析中，选用高斯指数函数作为AS-MFS 模型、 Co-Kriging 模型和HK 模型的相关函数。

图12 采用两种相关函数的不同变可信度模型对俯仰力矩系数建模的结果对比Fig.12 Comparison of predicted pitching moment coefficients via different variable-fidelity model based on the two correlation functions

如何评估最终建立的变可信度融合模型的预测精度，是研究人员最关心、也是其在工程问题中取得成功应用的关键性问题之一。本文沿用文献[2]提出的检验方法：交叉验证与物理规律检验相结合。交叉验证，即屏蔽一小部分试验数据，将其设定为未知，然后与模型预测进行比较。气动特性数据并非毫无规律可言，而是需要满足客观存在的物理规律，例如Oswatitsch 高超马赫数无关性原则。这些规律可以帮助研究人员对气动模型进一步确认预测结果。

鉴于项目数据空间分布的特殊性，高可信度风洞试验结果仅分布在马赫数4 以下的状态（图6）。本文将马赫数区间分成两段，分别采用交叉验证和物理规律检验的方法。具体而言：

1）对马赫数小于4 的内插区域，由于存在高可信度样本数据，因此随机屏蔽一小部分风洞试验数据点作为测试集，计算相关系数（r）、相对均方根误差（RRMSE）和相对最大绝对误差（RMAE） 3 个指标来表征模型的预估精度。其中，r和RRMSE 用来表征模型的全局预测精度，且r越接近1、RRMSE越接近0，模型全局预测越精确；RMAE 用来表征模型的局部预测精度，且RMAE 越接近0，模型局部预测越精确。3 个指标的计算公式如下：

其中，N 为测试集中的样本个数。需要说明的是，N 的值不宜过大，否则会导致需要被屏蔽的高可信度风洞试验数据过多，从而减少参与气动建模的高可信度样本数据，造成模型预测精度的损失。综合考虑后，从310 组高可信度样本里抽取了30 组样本作为测试集。

2）对于马赫数大于4 的外插区域，由于缺乏风洞试验数据，无法采用交叉验证方法进行精度评估，因而采用物理规律进行检验。需要说明的是，火箭一子级带栅格舵外形的飞行搭载试验取得了成功，其在再入启控的过程中经历了外插段的高马赫数区间，也从侧面间接说明了采用变可信度融合模型进行外插的正确性和有效性。后续研究工作中将发展对外插结果的直接验证方法。

2.3 模型的预测结果及对比

本节采用2.1 小节介绍的4 种变可信度融合模型分别对–20°、–10°、0°和10°舵偏角下栅格舵的俯仰力矩系数进行建模。表2为各舵偏角下的样本个数分布情况。

表2 不同舵偏角下样本分布情况Table 2 Data sets for different elevator deflections of grid fins

图13 为采用高斯指数相关函数的AS-MFS 模型、Co-Kriging 模型、HK 模型和MFNN 模型在马赫数[0.4，4]区间内对俯仰力矩系数的内插结果对比。从图中可以看到，AS-MFS 模型、Co-Kriging 模型和HK 模型的预测结果较为接近，均比较光滑，各舵偏角对应的预测曲面近似平行且相互隔开，不存在交叉现象。MFNN 模型预测的俯仰力矩系数在跨声速、大攻角区域有较高的隆起现象；在跨声速区，–20°和–10°舵偏角下的俯仰力矩系数有明显的交叉现象。图14 为马赫数为3.5 时采用交叉验证方法检验不同模型预测结果的对比情况。从图中可以看到，Co-Kriging 模型和HK 模型的预测值与测试点吻合最好，且二者精度较为接近。AS-MFS 模型的预测精度稍弱于前两种模型。由于AS-MFS、Co-Kriging 模型和HK 模型为插值模型，因此预测曲线均精确通过了训练样本；而MFNN 的预测结果较差，且由于它是拟合模型，因此未精确通过训练样本。此外，表3还展示了采用交叉验证方法计算的模型精度指标对比情况。可以看到，在–20°、0°和10°舵偏角下，Co-Kriging 模型的3 项精度指标均优于其他3 种模型。在–10°舵偏角下，Co-Kriging 模型除RRMSE 指标稍弱于HK 模型外，其他2 项精度指标仍为最佳。这表明Co-Kriging 模型在栅格舵俯仰力矩系数的内插预测上具备一定的精度优势。HK 模型的精度指标仅次于Co-Kriging 模型，也表现出了较强的预测能力。在4 种模型中，MFNN 模型预测结果的3 项指标均最差，俯仰力矩系数的预测误差较大。

图13 不同变可信度模型对俯仰力矩系数建模内插结果对比Fig.13 Comparison of predicted pitching moment coefficients via different VFMs

表3 针对俯仰力矩系数建模不同模型精度指标对比Table 3 Comparison of different accuracy indicators for pitching moment coefficients via different VFMs

图14 不同变可信度模型在马赫数为3.5 时俯仰力矩系数预测结果的交叉验证对比Fig.14 Comparison of cross validation for predicted pitching moment coefficients via different VFMs at Mach 3.5

图15 为4 种变可信度模型的外插预测结果对比，其中单独截取了马赫数为5 和7 下的俯仰力矩系数随攻角变化结果进行展示。从图中可以看到，AS-MFS 模型和MFNN 模型的外插曲面最光滑，相比之下，Co-Kriging 模型和HK 模型外插曲面的光滑程度略微下降，存在少量凹坑。在马赫数为5 的状态下，AS-MFS 模型预测的0°和–10° 舵偏角下的俯仰力矩系数曲线在攻角14°附近有部分重叠，HK 模型的俯仰力矩系数曲线在2 个舵偏角下则存在交叉现象。Co-Kriging 模型和MFNN 模型不同舵偏角下的俯仰力矩系数曲线均无交叉现象。在马赫数为7 的状态下，HK 模型预测的2 个舵偏角下的俯仰力矩系数曲线仍然交叉，而其余3 种模型的预测结果较好，符合基本物理规律。需要注意的是，在0°攻角、0°舵偏角状态下，俯仰力矩系数理论上过零点，但AS-MFS 模型和Co-Kriging 模型预测值在马赫数为5 和7 状态下整组预测曲线均向下略微偏移，俯仰力矩系数未能过零点；而HK 模型和MFNN 模型预测值在该状态时经过了零点，符合客观物理规律。考虑到HK 模型预测结果存在交叉现象，因此可以认为MFNN 模型在俯仰力矩系数外插时的预测效果更好。