0.   引　言

随着数值仿真、高性能计算等技术的发展，计算流体力学（Computational Fluid Dynamics, CFD）方法在航空、航天等领域得到了广泛应用^[1]。然而，由于知识不完备、试验样本匮乏和建模过程的模型简化，导致CFD模拟分析中不可避免地存在不确定性^[2]。例如，CFD模拟采用湍流模型计算机翼表面流场参数，但湍流模型参数往往来源于有限次的试验标定、经验假设^[3]，致使湍流模型参数具有不确定性，从而导致CFD模拟结果的不确定性。

传统CFD模拟不确定性量化方法将不确定性因素假设为概率分布或区间。Du等^[4]考虑马赫数和翼型指定位置升力系数服从均匀分布，利用效用理论（utility theory）构建RAE2822翼型气动优化问题的目标函数，获得了具有最小平均阻力系数和最小平均阻力系数标准差的翼型气动性能CFD模拟结果。Zhao等^[5]考虑k-ω SST湍流模型参数为区间数，分析了湍流模型参数对高速飞行器外表面压力系数和热通量的灵敏度。但是，上述方法并不能综合量化CFD模拟中的不确定性因素。

证据理论方法是对经典概率论和区间理论的扩展，采用更通用的框架描述不确定性信息。由于证据理论能够有效地表征随机和认知不确定性，且在满足一定条件时可与概率分布、区间数、模糊数和粗糙集等进行等效转化，近年来被广泛应用于结构不确定性量化和可靠性分析中^[6]。鉴于证据理论独特的优势，本文将利用证据理论量化CFD模拟中的不确定性。证据理论框架下的CFD模拟不确定性量化包含两大部分：输入参数不确定性量化、输出响应不确定性评估。前者将CFD模型输入参数表征为证据结构，后者又被称为不确定性传播。然而利用传统CFD模型进行证据变量传播的时间成本难以接受，因此需要开发证据理论框架下的代理模型主动学习策略实现证据理论框架下的CFD模型的不确定性传播。

本文提出了一种证据理论框架下主动学习代理模型驱动的CFD模拟不确定性量化方法，采用最优最大最小距离（optimization-based max-min distance, OMD）策略生成空间分布良好的候选样本点，设计了动态熵权-TOPSIS主动学习函数，综合了探索性、开发性和鲁棒性3个混合加点准则，兼顾了代理模型训练过程中的全局探索、局部开发和鲁棒性。此外，传统的主动学习代理模型在没有获得准确的全局解时进行不确定性传播，此时得到的结果误差较大且易产生错误的加点方向。为此，本文提出了基于证据信息的全局收敛准则，第一阶段只加点而不进行不确定性传播。加点范围为输入变量的最大取值范围，通过优化计算输出响应区间的中心和半径，当二者分别收敛时进入第二阶段。第二阶段结合了Hartley测度^[7]和Jousselme距离^[8]证据不确定性度量指标，以构成复合收敛准则，作为判断是否结束加点的依据。

最后，针对NASA SC(2)0410超临界机翼跨声速飞行流场CFD模拟，本文利用证据变量量化来流参数和湍流模型参数的不确定性，通过主动学习策略构建CFD代理模型，计算输出响应升阻比，并展示了证据理论相较其他不确定性表征数学方法的优越性。

1.   理论基础

1.1   证据理论

证据理论（evidence theory）最初被用于描述事件发生的可能性^[9]。记${\boldsymbol{\varOmega}} = \{ {{\boldsymbol{x}}_1},{{\boldsymbol{x}}_2}, \cdots ,{{\boldsymbol{x}}_N}\}$为辨识框架，${{\boldsymbol{x}}_i}{\text{ }}(i = 1,2, \cdots ,N)$表示事件。证据理论将事件发生的可能组合表示为幂集2^Ω：

式中，$\varnothing$表示空集，幂集中每一个元素或事件发生可能称为焦元（focal element, FE），记做A。各焦元发生的可能性用基本概率分配（basic probability assignment，BPA）$m(A)$表示，且满足以下条件：

与概率论不同的是，证据理论用信度函数$Bel(A)$和似真函数$Pl(A)$度量事件A发生可能性的下界和上界：

从而，$\left[ {Bel(A),Pl(A)} \right]$组成了关于事件$A$的不确定性大小的区间，且$Pl(A) - Bel(A)$表示事件A的认知不确定性大小。

1.2   代理模型

代理模型是用于替代复杂仿真模型的近似数学模型，通过少量已知的输入输出样本构建而成。在证据理论框架下，CFD模拟的不确定性传播涉及大量的极值优化计算，这些优化计算需与原始仿真模型进行大量迭代运算，造成巨大的计算负担。因此，构建CFD模拟的代理模型可以显著减少对复杂仿真的调用，从而有效地进行不确定性传播。常见的代理模型包括高斯过程代理模型^[10]、径向基函数代理模型^[11]和支持向量回归代理模型^[12]等。为比较不同模型的传播效果，本文将对这3种模型的不确定性量化结果进行定量对比分析。下面先对3种代理模型的数学理论进行介绍。

高斯过程（Gaussian process, GP）代理模型是一种基于贝叶斯统计的方法。它假设所有的函数点都是从一个相同的多元正态分布中抽取的，这使得不仅可以预测输入变量与输出变量之间的关系，还可以量化预测的不确定性，其假设函数的先验分布是高斯分布，即每一点的函数值均服从多元高斯分布。对于任意的输入样本集${\boldsymbol{X}} = {{\boldsymbol{x}}_1},{{\boldsymbol{x}}_2}, \cdots ,{{\boldsymbol{x}}_n}$，对应的输出${\boldsymbol{Y}} = {{\boldsymbol{y}}_1}, {{\boldsymbol{y}}_2}, \cdots ,{{\boldsymbol{y}}_n}$服从以下联合高斯分布：

式中：$\mu ({\boldsymbol{X}})$是均值函数，通常可以假设为零（$\mu ({\boldsymbol{X}}) = 0$）或基于数据拟合得到的线性模型；$\Sigma ({\boldsymbol{X}})$是协方差矩阵，由核函数$k({{\boldsymbol{x}}_i},{{\boldsymbol{x}}_j})$确定。核函数描述了输入空间中点与点之间的相似度，常用的核函数包括线性核、多项式核和指数型核等。由于平方指数核能够捕捉连续空间中平滑函数特性，并且拟合得到的代理模型足够光滑且无限可导，所以本文选取平方指数核函数作为计算协方差矩阵的核函数，即：

式中：${{\boldsymbol{x}}_i}$和${{\boldsymbol{x}}_j}$表示输入空间中的两个点，$\left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_j}} \right\|$表示样本点${{\boldsymbol{x}}_i}$和${{\boldsymbol{x}}_j}$之间的欧氏距离；$l$是长度尺度参数，可以控制函数的平滑程度，其值越大，生成的函数越平滑；${\sigma ^2}$表示振幅参数，其决定函数输出的变化幅度。对于新的输入点${{\boldsymbol{x}}^*}$，高斯过程预测输出${{\boldsymbol{y}}^*}$的分布为：

预测均值${\mu _ * }$和方差$\sigma _*^2$表示如下：

式中，${\boldsymbol{k}}({{\boldsymbol{x}}^*},{{\boldsymbol{x}}^*})$是新输入点的先验协方差，表明在该新输入点上的先验不确定性；${{\boldsymbol{k}}_*}$是新试验点${{\boldsymbol{x}}^*}$与训练数据点集X之间的协方差向量，${\boldsymbol{K}}$是训练数据点集的协方差矩阵，${\boldsymbol{I }}$表示单位矩阵，$\sigma _n^2$表示观测噪声的方差。

径向基函数（radial basis function, RBF）代理模型的任何一个预测点的函数值可以表示为所有中心点函数值的加权和，一般可以表示为：

式中，${\boldsymbol{x}}$是预测点，${{\boldsymbol{x}}_i}$是第$i$个中心点，${\lambda _i}$是待定系数，$\phi$是径向基函数，通常采用如下形式的高斯函数表示：

式中：$r = \left\| {{\boldsymbol{x}} - {{\boldsymbol{x}}_i}} \right\|$，表示预测点到中心点的欧氏距离；$\beta$表示形状参数，控制函数的平滑度。

支持向量回归（support vector regression, SVR）代理模型的基本思想可以理解为，求解一个函数$f(x) = {{\boldsymbol{w}}^{\rm{T}}}\phi ({\boldsymbol{X}}) + b$，使得对于所有训练数据$({{\boldsymbol{x}}_i},{y_i})$，误差$\left| {{y_i} - f({{\boldsymbol{x}}_i})} \right|$都不超过给定的误差阈值。其中${\boldsymbol{w}}$和$b$分别表示两个模型参数，$\phi ({\boldsymbol{X}})$是将输入数据映射到高维空间的函数。SVR的目标是最小化模型的复杂度，即${\left\| {\boldsymbol{w}} \right\|^2}$，同时保证大部分数据点的预测误差小于$\varepsilon$。因此通常通过求解以下优化问题实现：

式中：${\xi _i}$和$\xi _i^*$是松弛变量，用来处理不可能严格在边际$\epsilon$内的数据点；$C$是一个正则化参数，用来平衡模型的复杂度和训练误差。

2.   证据理论框架下CFD不确定性量化

CFD模拟结果的不确定性通常由模型形式不确定性（湍流模型、化学反应模型、状态方程等）、模型参数不确定性（模型参数、来流条件、外形参数等）、数值解偏差和模型预测以及后处理误差等引起。从不确定性来源角度来看，可以将上述不确定性分为随机不确定性和认知不确定性^[13]。

本文选用证据理论量化CFD模拟中的不确定性。证据变量是一种在连续空间取值的证据理论表述形式^[14]，证据理论采用了专家打分的思想，在试验数据的邻域划分焦元区间，并通过专家评判给每个焦元进行BPA赋值。例如，通过4次实验室和大气近地面试验测量，卡门常数$\kappa$的取值为0.39、0.41、0.43和0.45^[15]，由于实验数据较少，工程人员可以根据经验推测将其作为证据的区间焦元中心，其卡门常数的取值落在4个可能的区间上，从而得到关于卡门常数的辨识框架，如图1所示。

图  1  卡门常数的辨识框架

Figure  1.  Framework of discerning Karman constant

下载: 全尺寸图片 幻灯片

由于CFD模拟中变量的不确定性通常表现为在其名义值附近小范围的波动，即区间形式，所以本文中用于数值模拟不确定性传播的输入变量的辨识框架仅考虑区间形式。

一般来说，用于量化CFD模拟参数的不确定性证据变量由区间及其对应的BPA表示，包括贝叶斯证据变量、嵌套证据变量和重叠证据变量3种类型。值得注意的是，嵌套和重叠的证据变量可以转换为Bayesian证据变量，详见文献[16]。

针对图1所示的辨识框架，根据专家经验确定了不同焦元的BPA后得到的证据结构即为Bayesian证据变量。而输入证据变量的不确定性可通过CFD传播至输出响应，并利用响应的累积信度函数（cumulative belief function, CBF）和累积似真函数（cumulative plausibility function, CPF）量化响应结果的不确定性。假设函数$y = g({{\boldsymbol{x}}_1},{{\boldsymbol{x}}_2}, \cdots ,{{\boldsymbol{x}}_n})$表示从n维输入证据变量${{\boldsymbol{x}}_i}{\text{ }}(i = 1,2, \cdots ,n)$到输出响应的映射函数。其中证据变量${{\boldsymbol{x}}_i}$的焦元表示为$[{{x}_{i,{j_i}}}{\text{] }}({j_i} = 1,2, \cdots ,{n_{{{{\boldsymbol{x}}}_i}}})$。证据变量的不确定性通过在输入联合焦元构成的证据空间内进行极值运算得到输出响应焦元，实现传播。具体而言，输出的响应焦元$[{\underline {y}_k},{\bar y_k}]{\text{ }}(k = 1,2, \cdots , {n_{{{\boldsymbol{x}}_1}}} \times {n_{{{\boldsymbol{x}}_2}}} \times \cdots \times {n_{{{\boldsymbol{x}}_n}}})$通过求解${n_{{{\boldsymbol{x}}_1}}} \times {n_{{{\boldsymbol{x}}_2}}} \times \cdots \times {n_{{{\boldsymbol{x}}_n}}}$次优化模型得到，即：

式中，${\underline y _k}$和${\bar y_k}$表示第$k$个响应焦元的上界和下界。响应焦元$[{\underline y _k},{\bar y_k}]$的基本概率分配通过式（13）计算：

式中${m_{i \uparrow {{\boldsymbol{\varOmega}} _{1,2, \cdots, {n_{{{\boldsymbol{x}}_1}}} \times {n_{{{\boldsymbol{x}}_2}}} \times \cdots \times {n_{{{\boldsymbol{x}}_n}}}}}}}$表示将第i个变量的识别框架${\boldsymbol{\varOmega}} _i$映射到联合识别框架${\boldsymbol{\varOmega}}_{1,2,\cdots,n}$空扩展y运算。得到输出响应的证据结构后，分别通过累积信度函数和累积似真函数量化输出的不确定性^[17]：

式中：$\delta$表示狄利克雷函数，当${\bar y_k} \leqslant y$或者${\underline y _k} \leqslant y$时，函数值取1，否则为0；$y$表示响应范围${\varOmega _y}$内的离散变量。

3.   面向全局近似的主动学习代理模型

主动学习加点策略直接影响代理模型建模效率和建模精度，因此加点准则至关重要，本文采用多属性混合加点准则进行新增样本点选择。所提出的证据理论框架下不确定性传播代理模型的建模策略基本思想如下：首先，在CFD模拟输出证据空间内进行少量的均匀采样，选取初始训练样本点构建初始高斯过程代理模型。然后，对每个训练样本点的未探索程度、非线性程度和鲁棒性进行评估，并利用这些属性的熵信息依据TOPSIS决策^[18]模型确定权重矩阵，以识别敏感样本区域。接着，采用委员会查询策略（query by committee, QBC）添加新的样本点。

由于初始代理模型精度较低，直接利用其进行不确定性传播结果不准确，因此本文提出基于证据信息的全局收敛准则以提高传播效率。首先，在第一阶段的迭代过程中，计算全局优化的最大值和最小值组成第一阶段的迭代特征向量。当迭代特征向量收敛后，表明代理模型的全局精度已足够，此时才进行不确定性传播。随着迭代继续进行，进一步通过计算每次迭代得到的证据结构的Hartley测度和相邻两次迭代证据结构的Jousselme距离组成新的迭代特征向量。最后，通过比较迭代特征向量的相对修正结果作为终止准则，从而实现利用少量的试验点实现CFD模拟不确定性的高效传播。

3.1   混合加点准则

对于一个n维响应函数$y = g({X_1},{X_2}, \cdots ,{X_n})$，根据Loeppky等^[19]提出的初始样本点数经验公式$m = 4n$，利用拉丁超立方采样（Latin hypercube sampling, LHS）生成$m$个初始样本并得到对应的输出，从而获得初始样本集合${\boldsymbol{C}}=\left\{\left({\boldsymbol{x}}_{j}, {\boldsymbol{y}}_{j}\right) \big| {\boldsymbol{y}}_{j}=g\left({\boldsymbol{x}}_{j}\right),\right. \left.j=1,2, \cdots, m\right\}$和初始的高斯过程代理模型。

3.1.1   基于Voronoi单元相对大小的探索性指标

样本点的探索性是量化已知样本点附近样本密度的指标，是对模型训练点的未知程度的量化。较低的探索性指标表明训练点附近的采样数量较低，因此需要重点考虑增加采样点。Voronoi图^[20]是一种被广泛应用于计算机图形学、路径规划等领域的数学结构，它基于一组给定的点在平面（空间）内划分区域，使得每个点的区域内的任何位置都比其他点的区域更靠近该点。因此，Voronoi单元的相对大小成为量化每个样本点未探索区域大小的理想指标。然而，Voronoi划分通常会产生不规则的多边体，这使得计算Voronoi单元的相对大小较为困难，传统方法多采用LHS生成大量候选样本点，通过比例近似来估计Voronoi单元的相对大小。针对高维问题，这种方法不仅需要产生大量的候选样本点，而且生成的候选点邻近性较高，导致Voronoi单元的评估效率和精度较低，并且会降低后续新增样本点的质量。因此本文提出最优最大最小距离策略来迭代生成空间分布良好的候选样本点。该策略通过最大化候选样本点与已知训练样本的距离来生成候选样本。对于初始训练集$\mathcal{D}(0) = \{ {{\boldsymbol{x}}_1},{{\boldsymbol{x}}_2}, \cdots ,{{\boldsymbol{x}}_{{n_{\mathcal{D}(0)}}}}\}$，对其进行归一化处理得到$\mathcal{D}'(0) = \{ {{\boldsymbol{x}}'_1},{{\boldsymbol{x}}'_2}, \cdots ,{{\boldsymbol{x}}'_{{n_{\mathcal{D}(0)}}}}\}$，设置每次迭代生成的候选样本点数目为${N_{{\mathrm{Cand}}}}$，本文取${N_{{\mathrm{Cand}}}} = 80$，其值随着问题维数的增加可适当增大。同时，令${\mathcal{L}'_C}(0) = \varnothing$，${\mathcal{D}_{{\mathrm{known}}}} = \mathcal{D}'(t - 1)\bigcup {{{\mathcal{L}'}_C}(t - 1)}$，通过在归一化空间${\mathbb{I}_X}$内迭代求解对优化问题${D_{\min }}({\boldsymbol{x}}') = \mathop {\min }\limits_{{\boldsymbol{x}}' \in {\mathbb{I}_X},{{{\boldsymbol{x}}'}_{i}} \in {\mathcal{D}_{{\mathrm{known}}}}} \left\| {{\boldsymbol{x}}' - {{{\boldsymbol{x}}'}_{{\mathrm{known}}}}} \right\|$${N_{{\mathrm{Cand}}}}$次后，将得到的${N_{{\mathrm{Cand}}}}$个样本点重映射到原来的样本空间后，得到候选样本点集${\mathcal{L}_C}(t)$。图2给出了提出的候选样本点策略与传统的LHS方法的采样结果。

图  2  提出的候选样本点生成策略效果图

Figure  2.  Results of the proposed candidate sample generation strategy

下载: 全尺寸图片 幻灯片

从图2中可以看出，所提方法能在少量的样本点下近似出Voronoi单元的相对面积，并且生成的样本点在每个Voronoi单元内都具有较低的邻近性，为后续新增样本点的选择做下了铺垫。在得到每个已知训练点的Voronoi单元和候选样本点集后，通过将候选样本集${\mathcal{L}_C}$拆分为$m$个子集，每个子集${{\boldsymbol{Q}}_i}{\text{(}}{{\boldsymbol{x}}_i}{\text{), }}i = 1,2, \cdots ,m$的元素由候选样本集中最邻近第$i$个训练样本${{\boldsymbol{x}}_i}$的点组成。因此，每个训练点的探索指标${\mathcal{P}_{\rm{ue}}}({{\boldsymbol{x}}_i})$表示为：

式中，${N_{\boldsymbol{Q}}}_i$和${N_{{\mathcal{L}_C}}}$分别表示集合${{\boldsymbol{Q}}_{\boldsymbol{i}}}$和${\mathcal{L}_C}$的基数。

3.1.2   基于局部线性近似误差的开发性指标

样本点的开发性是对样本点附近的特征进行描述，如非线性程度、预测方差和累积变化误差等，因此该属性是对模型训练点的已知信息的量化。由于证据理论下不确定性传播涉及大量极值优化问题，而最优解往往出现在非线性程度区域，因此本文利用每个训练点的局部线性近似误差${\mathcal{P}_{\rm{Ep}}}({{\boldsymbol{x}}_i})$作为样本点的开发性指标。具体计算步骤如下：首先，定义每个训练点${{\boldsymbol{x}}_i}$的邻域输入集${\mathcal{D}_{{\mathrm{neg}}}}({{\boldsymbol{x}}_i}) = \{ {{\boldsymbol{x}}_{i1}},{{\boldsymbol{x}}_{i2}}, \cdots ,{{\boldsymbol{x}}_{iv}}\} \subset \mathcal{D}$，该邻域集由训练样本集中离训练点${{\boldsymbol{x}}_i}$最近的$v = 2n$个训练样本组成。训练点${{\boldsymbol{x}}_i}$处的梯度估计$\hat {\boldsymbol{g}}({{\boldsymbol{x}}_i})$利用最小二乘法求解一个超定方程组获得：

式中，${{\boldsymbol{y}}_{ij}}$表示集合$\mathcal{C}=\left\{\left({\boldsymbol{x}}_{i}, {\boldsymbol{y}}_{i}\right)\right\}_{i\in[1, 2, \cdots ,m]}$中与${{\boldsymbol{x}}_{ij}}$对应的元素。求解得到的训练点${{\boldsymbol{x}}_i}$处的近似梯度表示为$\hat {\boldsymbol{g}}({{\boldsymbol{x}}_i}) = [\hat {\boldsymbol{y}}_i^1,\hat {\boldsymbol{y}}_i^2, \cdots ,\hat {\boldsymbol{y}}_i^n]$，从而训练点${{\boldsymbol{x}}_i}$的开发指标${\mathcal{P}_{\rm{Ep}}}({{\boldsymbol{x}}_i})$可通过局部线性近似值与真实映射函数值之间的误差表示，即：

式中${{\boldsymbol{x}}_{ij}}$表示邻域集${\mathcal{D}_{{\mathrm{neg}}}}({{\boldsymbol{x}}_i})$的第$j$个元素。

3.1.3   基于留一交叉验证误差的鲁棒性指标

样本点的鲁棒性是指模型缺失某个样本点后受影响的水平。当缺失某个样本点后，计算得到的留一交叉验证误差较低，说明该样本点不会显著影响模型精度，因此说明该样本点周围的插值是鲁棒的。每个已知训练点的鲁棒性指标留一交叉验证（leave-one-out cross validation, LOOCV）误差表示为${e_{\rm{LOOCV}}}({{\boldsymbol{x}}_i})$，具体计算公式如下：

式中，$\hat G({{\boldsymbol{x}}_i})$表示利用已知训练点构建的代理模型，${\left\{ {{{\hat G}_{ - i}}} \right\}_{i \in [1,2, \cdots ,m]}}$表示将训练点${{\boldsymbol{x}}_i}$从训练集中除去后构建的代理模型。

3.2   动态熵权-TOPSIS下的敏感样本点确定

在计算得到已知训练点的探索性、开发性和鲁棒性指标后，可以利用这3个指标确定敏感样本点。为避免迭代过程过度依赖于某一指标带来的信息，本文提出动态熵权-TOPSIS主动学习策略用于确定敏感训练样本点。该策略在获得每个训练点的探索性、开发性和鲁棒性指标后，首先构建如下归一化决策矩阵${{\boldsymbol{D}}_{\rm{norm}}}$：

作为TOPSIS多属性决策的关键步骤，如何构造标准化加权向量是关键，对于每一个指标，低熵值表明该指标具有较高的确定性，参考可信性高。根据这一性质，提出基于熵信息动态权重${\boldsymbol{\omega}} = [{w_1},{w_2},{w_3}]$的计算方法，在第$t$轮迭代中，每个指标的熵信息可由式（21）计算：

式中：$k$表示迭代过程的训练样本点总数，${a_{ij}}$表示决策矩阵第$i$行第$j$列对应的元素。基于熵信息，第$t$轮迭代过程的不同指标的权重定义如式（22）：

在为每个属性赋予相应权重之后，将所有属性值都达到最优的情况定义为理想解，而所有属性值都处于最劣的情况定义为负理想解。接着，计算每个样本点属性向量与理想解和负理想解的欧氏距离，从而得到每个样本点的相对接近度。这一接近度是样本点属性向量到理想解的距离与到理想解和负理想解距离之和的比值。最后，选择样本点的相对接近度较大的点作为敏感训练样本点。

3.3   基于证据信息的两段式判断准则

获得敏感训练样本点后，新增样本点则在敏感样本点对应的Voronoi单元区域内产生，通过QBC策略从候选样本点集中产生新的试验样本点${{\boldsymbol{x}}_{\rm{new}}}$。利用探索性指标计算过程中对候选样本集的拆分结果，采用最远距离策略选择新试验点${{\boldsymbol{x}}_{\rm{new}}}$如下：

从上式中可以看出，当敏感样本点对应的Voronoi单元区域内存在候选样本点时，选择该集合中距离敏感样本点最远的候选样本作为新增样本点。随着迭代的进行，部分非线性程度较高的区域的样本点将变得较为密集，这可能导致敏感区域内无可用的候选样本点。在这种情况下，本文会选择整个候选样本集中距离敏感样本点最近的点作为新增样本点。由于采用OMD策略生成的候选样本点具有良好的空间分布特性，新增的样本点将仍与已有样本点集保持较低的邻近性。

由于初始代理模型精度较低，导致直接利用初始代理模型进行不确定性传播的结果不准确，因此，本文将代理模型训练到一定精度时再进行传播，以提高传播效率；接着再利用证据信息作为整个传播过程的终止准则。

本文两段迭代判断准则都采用连续迭代特征向量相对修正误差作为判断依据，具体表示如下：

式中，${{\mathcal{{K}}}_t} = [I_t^c,I_t^r]$表示迭代特征向量，${\boldsymbol{\tau}} = {\tau _s} \times [1,1]$表示相对修正容差向量。${\tau _s}$表示容差阈值，通常取很小的正数值，在本文中取0.0001。

在构建初始近似代理模型阶段，通过在证据空间内对迭代更新得到的代理模型进行全局优化计算，得到目标响应的极值响应区间${I_t} = [{\hat G_{\min }}({\boldsymbol{x}}),{\hat G_{\max }}({\boldsymbol{x}})],$${\boldsymbol{x}} \in \mathbb{X}$；然后将其转化为中心-半径的形式表示，其中，区间中心表示为$I_t^c = \dfrac{{{{\hat G}_{\min }}({\boldsymbol{x}}) + {{\hat G}_{\max }}({\boldsymbol{x}})}}{2}$，半径表示为$I_t^r = \dfrac{{{{\hat G}_{\max }}({\boldsymbol{x}}) - {{\hat G}_{\min }}({\boldsymbol{x}})}}{2}$。此时，判断准则中的迭代特征向量由区间中心和半径组成。

在得到具有一定精度的代理模型后，则可以利用所构建的代理模型进行证据变量的不确定性传播，当所有的输入焦元利用模型进行预测得到响应焦元后，则可以构建得到目标响应的证据结构。本文采用Jousselme距离和Hartley测度定义基于证据信息的全局收敛准则，Jousselme距离和Hartley测度定义如下：

令${{\boldsymbol{m}}_1}$和${{\boldsymbol{m}}_2}$分别表示两个证据结构，两个证据结构的Jousselme距离指标定义为：

式中，${\boldsymbol{J}}$表示维度为${2^{\boldsymbol{\varOmega}} } \times {2^{\boldsymbol{\varOmega}} }$的Jaccard冲突矩阵，具体定义如下：

令${\boldsymbol{m}}$表示输出响应辨识框架$Y$对应的基本概率分配集合，该证据结构的Hartley测度定义为：

式中，$\mu (A)$是对命题A的信息量大小的量化函数，当响应焦元在实数域$R$表现为区间$\left[ {\underline y ,\bar y} \right]$形式时，此时，$\mu (A) = \bar y - \underline y$。

在证据理论框架下，相邻两次迭代得到的证据结构之间Jousselme距离是对迭代结果相似度的量化，单个证据结构的Hartley测度是对该证据结构包含的不确定性信息的量化。随着迭代过程的进行，传播得到的证据结构趋于稳定，不会出现新增信息，此时Jousselme距离和Hartley测度组成的迭代特征向量也稳定收敛。因此，将相邻两次迭代结果的Jousselme距离表示为${E_{{\text{Jd}}}}$，每次迭代结果的Hartley测度表示为${e_{{\text{non}}}}$，同时采用式（24）所示的判断准则，此时迭代特征向量表示为${\mathcal{K}_t} = [{E_{{\text{Jd}}}},{e_{{\text{non}}}}]$。当判断准则满足时，得到证据变量的传播结果。

4.   算例分析

4.1   数值算例

为了说明所提出算法的性能优势，首先选用修正Branin函数开展数值算例分析。修正Branin函数^[21]具有两个局部最小值和一个全局最小值，其具体表达式如下：

式中，${\bar X_1}$、${\bar X_2}$分别表示对两个独立的证据变量${X_1}$和${X_2}$做如下尺度变换：

式中，修正后的变量${\bar X_1}$和${\bar X_2}$的辨识框架范围均为$[0,1]$，该函数的三维示意图如图3所示。

图  3  放缩后的Branin函数的示意图

Figure  3.  Schematic diagram of rescaled Branin function

下载: 全尺寸图片 幻灯片

假设每个证据变量的证据结构内包含8个区间焦元，则该算例共包含64个联合焦元，其具体的BPA结构如表1所示。

表  1  证据变量${X_1}$和${X_2}$的BPA结构

Table  1.  BPA structure for evidence variables

${X_1}$			${X_2}$
区间焦元	BPA		区间焦元	BPA
[0,0.125]	0.020		[0,0.125]	0.050
[0.125,0.250]	0.090		[0.125,0.250]	0.050
[0.250,0.375]	0.163		[0.250,0.375]	0.150
[0.375,0.500]	0.247		[0.375,0.500]	0.200
[0.500,0.625]	0.239		[0.500,0.625]	0.200
[0.625,0.750]	0.141		[0.625,0.750]	0.150
[0.750,0.875]	0.083		[0.750,0.875]	0.100
[0.875,1.000]	0.017		[0.875,1.000]	0.100

下载: 导出CSV 

| 显示表格

图4给出了本文所提主动学习策略的示例。首先在证据空间内利用拉丁超立方采样选择8个初始样本点，利用提出的OMD候选样本点生成策略迭代产生候选样本点集，然后利用各Voronoi单元内候选样本点比例评估各训练样本的全局探索性，最后利用基于近似梯度的局部线性近似误差评估局部开发性和留一交叉验证误差评估的鲁棒性。图中每个训练样本对应的玉玦图中的每一环表示不同指标的相对大小，在利用动态熵权-TOPSIS决策后，得到敏感样本点及其对应Voronoi单元，选择单元内距离敏感样本点最远的样本点作为新增样本点。

图  4  提出的主动学习策略的示例

Figure  4.  An example for the proposed strategy

下载: 全尺寸图片 幻灯片

为了说明当满足主动学习判断准则时，此时训练得到的代理模型在整个证据空间内具有较好的精度，图5（a）展示了满足主动学习阶段判断准则时，样本点在证据空间内的分布情况，图5（b）给出了此时构建得到的代理模型的相对误差热力图。可以看出，此时在整个空间内的相对误差都已经达到一个较小的值，并且新增的样本点都能很好地遍布在空间区域内，特别是利用所提出的学习策略得到的样本点在非线性程度较高的区域也具有良好的分布性。因此，当满足第一段主动学习的判断准则时，此时利用代理模型进行传播得到的结果已经具有一定的可参考性。并且这种两阶段式迭代判断准则能够一定程度地减少迭代过程证据变量传播所需的极值优化时间，这种策略在输入变量的联合焦元数量规模比较大的情况下更具优势。

图  5  提出的主动学习策略的细节和效果

Figure  5.  Details and effect of the proposed active learning strategy

下载: 全尺寸图片 幻灯片

当迭代48次后，第一段主动学习过程的判断准则得到满足，再迭代13次后，基于证据信息的收敛准则得到满足。为了对比不同代理模型在本文所提出的主动策略下的建模效果，利用本文所提学习策略学习得到的样本集建模不同的代理模型，并利用不同的模型传播得到输出响应的证据结构。图6中对比了采用相同训练不同代理模型后传播得到的输出响应不确定性表征结果，在图中利用输出响应的累积信度函数和累积似真函数进行定量表征。

图  6  不同代理模型传播得到的CBF和CPF

Figure  6.  Results of CBF and CPF achieved by different surrogate models

下载: 全尺寸图片 幻灯片

表2中给出了不同结果与参考方法之间的定量对比结果，可以看出，高斯过程代理模型的精度比其他模型的精度更具优势，这得益于所提学习策略学习得到的采样点在空间内良好的分布特性和学习函数能够平衡不同样本指标带来的信息。

表  2  不同代理模型定量对比结果

Table  2.  Quantitative comparison results of different surrogate models

方法	与参考解之间的Jousselme距离	Hartley测度
参考方法	—	0.7604
基于GP	0.9997	0.7605
基于RBF	0.9551	0.7548
基于SVR	0.8982	0.7369

下载: 导出CSV 

| 显示表格

4.2   超临界机翼CFD模拟

以剖面翼型为NASA SC(2)0410的超临界机翼风洞缩比模型跨声速流场CFD数值模拟为例，说明本文所提算法的有效性以及证据理论在量化CFD模拟不确定性上的重要作用。由于该机翼模型具有独特的气动特性，使其在提高飞行效率方面具有显著优势，因此被广泛应用于关于超临界机翼的研究中。Ghanooni等^[22]在低雷诺数条件下对不同倾斜角度的小翼对该型号超临界机翼气动性能的影响进行了系统研究，结果表明，可以通过减小小翼翼尖涡流的大小和强度实现增升减阻。机翼通常由多个翼型剖面构成，本文以大展弦比、中等后掠双段机翼为例进行不确定性分析，图7展示了剖面翼型为NASA SC(2)0410的后掠机翼三视图和机翼结构变量参数^[30]。为了方便展示，图中右侧两幅视图的左侧的视图放大了2倍。半展长为800 mm，翼根剖面翼型弦长为200 mm，平均气动弦长为130.64 mm，下反角为3.5°。图8（a） 展示了NASA SC(2)0410超临界机翼的远场网格划分结果。该网格先以结构网格形式进行划分，为了兼容CFD模拟求解器的非结构化网格需求，其后又利用网格划分软件将其转换为非结构网格进行求解。图8（b）是对机翼进行网格划分的结果。远场尺寸为12000 mm×7000 mm×5000 mm，单元数共计2563241个，其中机翼表面采用O-block拓扑网格划分，共计46034个单元。

图  7  剖面翼型为NASA SC(2)0410的超临界机翼三视图

Figure  7.  The orthographic projection of the supercritical wing profile of NASA SC (2) 0410

下载: 全尺寸图片 幻灯片

图  8  NASA SC(2)0410超临界机翼表面及远场网格划分

Figure  8.  The grid setup for NASA SC(2)0410 supercritical wing and its farfield

下载: 全尺寸图片 幻灯片

在对网格划分质量的评估中，网格划分质量用${y^ + }$表示，一般要求机翼表面${y^ + } < 1$以满足线性黏性底层要求^[23]。机翼表面${y^ + }$分布情况如图9所示，最大值为0.036，远小于1，说明网格划分质量较高。

图  9  NASA SC(2)0410超临界机翼表面$y^+$值分布

Figure  9.  The $y^+$ distribution of NASA SC(2)0410 supercritical wing

下载: 全尺寸图片 幻灯片

在本算例中，利用ANSYS/Fluent 2021R1的压力-速度耦合求解器和Spalart-Allmaras （S-A）湍流模型，模拟NASA SC(2)0410超临界机翼跨声速飞行状态时的表面流场，并求解关键气动性能升阻比${{{C_L}} \mathord{\left/ {\vphantom {{{C_L}} {{C_D}}}} \right. } {{C_D}}}$。其可以表示为：

式中，${g_{{\text{CFD}}}}$表示CFD仿真模型，${{\boldsymbol{\theta}} _1}$和${{\boldsymbol{\theta}} _2}$分别表示来流参数和湍流模型参数。本算例中，来流参数${{\boldsymbol{\theta}} _1}$包括来流马赫数$Ma$、迎角α、侧滑角β和飞行海拔高度H，${{\boldsymbol{\theta}} _2}$为S-A湍流模型参数，其取值源于圆管流、槽道流、平板流等简单流动标定或依据专家经验给出，其试验标定值和经验取值范围^[24]如表3所示。

表  3  S-A湍流模型参数

Table  3.  Parameters of S-A turbulence model

湍流模型参数	试验标定值[24]	经验取值范围
${c_{b1}}$	0.1355	[0.12893, 0.137]
${c_{b2}}$	0.622	[0.60983, 0.6875]
${\sigma _{\tilde \upsilon }}$	0.667	[0.6, 1]
${c_{v1}}$	7.1	[6.9, 7.3]
${c_{w2}}$	0.3	[0.055, 0.3525]
${c_{w3}}$	2	[0.175, 2.5]
$\kappa$	0.41	[0.38, 0.46]

下载: 导出CSV 

| 显示表格

此外，由于跨声速飞行过程可能出现的大气紊流以及飞行过程中飞机受扰动出现的沉降，在满足跨声速马赫数和不失速条件下，表4中给出了来流参数的取值范围^[25]。

表  4  来流参数

Table  4.  Incoming flow parameters

来流参数	经验取值范围
Ma	[0.82, 0.849]
$\alpha$/(°)	[3, 3.2]
$\beta$/(°)	[0, 1]
H/m	[9990, 10000]

下载: 导出CSV 

| 显示表格

以来流参数马赫数的不确定性为例，图10展示了不同来流马赫数下翼根剖面绝对压力的变化。随着来流马赫数增加，上方翼面的低压区（$\leqslant 7 \times 10^4 \;{\mathrm{Pa}}$）和下方翼面的高压区（$\geqslant 8 \times 10^4 \;{\mathrm{Pa}}$）扩大，而升阻比为上下表面压力差和阻力系数的函数。因此，来流参数马赫数对升阻比的预测将产生不可忽略的影响，所以在对机翼气动性能进行评估的过程中，考虑各种参数的不确定性对结果进行可信度评估至关重要。

图  10  翼根剖面绝对压力随来流马赫数变化

Figure  10.  The variation of absolute pressure contour on the wing root with freestream Mach number

下载: 全尺寸图片 幻灯片

为了说明证据理论在少量样本下对输入不确定性表征的优势，基于表3和表4的模型参数经验取值范围，在表5中给出了来流和湍流模型参数的几组可能取值，用于模拟在实际风洞试验和试飞试验中只能获取有限个来流、湍流模型参数样本的情况，分析不同不确定性表征方法对升阻比预测结果的影响。

表  5  来流和湍流模型参数取值

Table  5.  The incoming flow and turbulence model’s parameters provided by wind-tunnel experiments

模型参数	参数取值
Ma	0.824, 0.8315, 0.8385, 0.8455
$\alpha$/(°)	3.025, 3.075, 3.125, 3.175
$\beta$/(°)	0.125, 0.375, 0.625, 0.875
H/m	9991, 9993.5, 9996.5, 9999
${c_{b1}}$	0.13, 0.132, 0.134, 0.136
${c_{b2}}$	0.6182, 0.6375, 0.6575, 0.6775
${\sigma _{\tilde \upsilon }}$	0.65, 0.75, 0.85, 0.95
${c_{v1}}$	6.95, 7.05, 7.15, 7.25
${c_{w2}}$	0.1052, 0.2055, 0.304
${c_{w3}}$	1.875, 2.125, 2.375
$\kappa$	0.39, 0.41, 0.43, 0.45

下载: 导出CSV 

| 显示表格

在表5所示的这种少量样本条件下，随机不确定性框架通常利用分布函数方法拟合经验概率密度函数来表征输入来流和湍流模型参数的不确定性^[26]，然而这种估计方法依赖于样本规模、统计方法和参数估计方法。核密度估计是一种用于估计随机变量的平滑概率密度函数的非参数方法，并可以通过Bootstrap方法计算估计结果的置信区间，因此被广泛应用于试验样本的随机不确定性表征，但是其表征结果的置信度与样本数量具有强相关性^[27]。以迎角为例，利用核密度估计方法估计得到的迎角不确定性表征结果如图11所示。

图  11  迎角的核密度估计及置信区间

Figure  11.  Angle of attack kernel density estimation and confidence intervals

下载: 全尺寸图片 幻灯片

从图11中可以看出，拟合得到的经验概率密度的置信区间较宽，且取值范围会存在超出经验取值范围的情况，从而违背经验信息。此时进行不确定性传播会导致评估结果可信度低，失去参考价值。

本文通过证据变量表征并量化输入来流与湍流模型参数的不确定性，然而确定证据变量焦元的BPA取值作为利用证据理论进行不确定性量化的核心步骤，一直以来都是一个悬而未决的问题。梁宵等^[28]提出对于无物理意义、通过试验直接标定的参数，在考虑不同专家工程经验和知识不完备的情况下，可以利用证据理论的多证据融合规则得到认知不确定性量化的结果，从而获取不同参数赋值区间的基本概率分配值。Tonon^[29]提出了基于参数离散化的方法，将概率盒变量转换为证据结构，从而获得区间焦元BPA。本文采用Tonon提出的参数离散化BPA确定方法，再结合文献[30]对来流参数和湍流参数的不确定性表征结果，通过将输入不确定性信息的度量转化处理后得到来流参数和湍流模型参数的证据变量表征结果，如表6所示。表6共给出11个湍流模型参数和来流参数的证据结构，通过焦元的笛卡尔积运算后得到联合焦元总数，共计${4^9} \times {3^2} = 2\;359\;296$。

表  6  证据理论框架下来流和湍流模型参数不确定性量化结果

Table  6.  Uncertainty quantification of incoming flow and turbulence model’s parameters under the framework of evidence theory

Ma			$\alpha$/(°)			$\beta$/(°)			H/m			${c_{b1}}$			${c_{b2}}$
FE	BPA		FE	BPA		FE	BPA		FE	BPA		FE	BPA		FE	BPA
[0.820, 0.828]	0.34		[3.00,3.05]	0.25		[0.00,0.25]	0.16		[9990,9992]	0.2		[0.12893,0.131]	0.13		[0.60983,0.6275]	0.4
[0.828, 0.835]	0.22		[3.05,3.10]	0.25		[0.25,0.50]	0.34		[9992,9995]	0.3		[0.131,0.133]	0.37		[0.6275,0.6475]	0.2
[0.835, 0.842]	0.22		[3.10,3.15]	0.25		[0.50,0.75]	0.34		[9995,9998]	0.3		[0.133,0.135]	0.37		[0.6475,0.6675]	0.2
[0.842, 0.849]	0.22		[3.15,3.20]	0.25		[0.75,1.00]	0.16		[9998,10000]	0.2		[0.135,0.137]	0.13		[0.6675,0.6875]	0.2
${\sigma _{\tilde \upsilon }}$			${c_{v1}}$			${c_{w2}}$			${c_{w3}}$			$\kappa$
FE	BPA		FE	BPA		FE	BPA		FE	BPA		FE	BPA
[0.6,0.7]	0.7		[6.9,7.0]	0.1		[0.055,0.1555]	0.2		[1.75,2.00]	0.4		[0.38,0.40]	0.14
[0.7,0.8]	0.1		[7.0,7.1]	0.4		[0.1555,0.2555]	0.3		[2.00,2.25]	0.4		[0.40,0.42]	0.36
[0.8,0.9]	0.1		[7.1,7.2]	0.4		[0.2555,0.3525]	0.5		[2.25,2.50]	0.2		[0.42,0.44]	0.36
[0.9,1.0]	0.1		[7.2,7.3]	0.1		—	—		—	—		[0.44,0.46]	0.14

下载: 导出CSV 

| 显示表格

进一步地，利用本文所提出的方法构建证据理论框架下高斯过程主动学习代理模型。首先，选择44个初始训练样本构建初始高斯过程代理模型。调用64次CFD模型使输出响应升阻比的最大范围区间中点和半径收敛。利用式（12）进行证据变量的不确定性传播，得到输出响应升阻比的焦元和BPA值，通过式（25）和式（27）计算Hartley测度和Jousselme距离，共加点183个后达到收敛指标。模型收敛后，采用留一交叉验证相对误差对建立的代理模型的精度进行评估，此时得到的交叉验证相对误差为$1.391\;7 \times {10^{ - 5}}$，表明训练得到的模型具有较高的精度。图12展示了计算得到输出响应升阻比的证据结构。从图中可以看出，利用证据理论得到的输出响应的证据结构中，每一个响应的可能取值范围都给出了似真度和信度量化不确定性。相对于传统的概率框架下得到的累积概率分布，本文方法能够提供更具参考价值的结果。

图  12  证据理论框架下输出响应升阻比

Figure  12.  The output response lift-to-drag ratio under the framework of evidence theory

下载: 全尺寸图片 幻灯片

5.   结　论

本文针对CFD模拟过程存在的多类型不确定因素，利用证据理论统一量化其随机与认知不确定性，并提出方法解决证据理论框架下不确定性传播效率低的问题，同时应用于翼型升阻比预测结果不确定性量化，很好地实现了CFD模拟可信度评估。主要结论如下：

1） 建立了一种证据理论框架下的CFD不确定性量化与传播方法。利用证据变量量化输入不确定性，利用累积信度函数和累积似真度函数量化输出变量的不确定性，能够更好地评估CFD模拟结果的可信度。

2） 提出了一种适用于证据变量传播的主动学习策略。该策略使用OMD策略生成候选样本，并使用基于最大距离原则的QBC策略选择新的试验样本。同时，该策略能够综合全局探索性、局部开发性和稳健性指标提供的信息，用以选择敏感样本点，最后基于Hartley测度和Jousselme距离作为整个迭代过程终止的判断准则。

3） 通过数值算例和超临界机翼算例说明了该方法可以在试验点数较少的非线性空间中选择新样本构建代理模型，提出的方法对于构建可信的CFD模拟不确定性传播方法起到了理论支撑作用。

然而在本研究中，不同证据变量的传播是在相互独立的假设下进行的。在相互依赖的条件下证据变量的高效传播是一个巨大的挑战，也是未来研究的一个方向。此外，在证据理论框架内探索不确定性的反向传播也值得进一步研究。