飞行器气动外形直接决定其飞行性能与飞行品质。气动外形与性能之间存在复杂关系，并受多种约束条件限制，使得气动外形设计成为整个飞行器设计流程中最为困难的环节。气动外形数字化设计技术充分利用计算流体力学、数值优化方法，在提高设计效率的同时，有助于获得更优异的气动性能。气动外形数字化设计技术包括优化设计方法和反设计方法。优化设计方法鲁棒性好，对设计经验依赖程度低，但依赖于数值优化算法的特性，如全局优化算法需要大量的流场分析，计算效率低，而基于梯度的优化设计结果强烈依赖于初始构型在设计空间中的位置。

反设计方法通过改变气动外形，使得其表面压力分布与给定的目标压力分布尽量吻合，该方法通常计算效率较高，在给定优秀目标压力分布后能够很快得到满意的气动外形。Lighthill^[1]最早提出了基于保角变换的翼型反设计方法。之后，Takanashi^[2]发展了余量修正法及其改进方法^[3-4]，该方法是工程常用的反设计方法之一。此类经典的反设计方法要求给出合理的目标压力分布，否则无法收敛。Bui-Thanh等^[5]利用本征正交分解(proper orthogonal decomposition, POD)降阶方法发展了Gappy POD翼型反设计方法。白俊强等^[6]对Gappy POD进行了改进，引入了迭代校正法提高反设计精度。Gappy POD反设计方法对目标压力分布合理性没有要求，计算效率高，但是不能处理各种几何和气动约束。

近年来，机器学习技术不断发展并被应用到翼型反设计领域。单志辉^[7]基于高斯过程回归模型进行了翼型反设计研究，压力分布和外形数据用于模型训练，以预测目标压力分布对应的翼型。Sun等^[8]基于人工神经网络开展了翼型反设计研究。Sekar等^[9]利用深度卷积神经网络将翼型和压力分布的图像作为学习对象，实现了压力分布图像特征提取，从而计算获得翼型外形。目前，基于神经网络，尤其是深度神经网格的反设计方法还存在模型训练复杂、训练时间长等问题。

生成拓扑影射(generative topographic mapping, GTM)^[10]是机器学习领域的一种非线性隐变量模型，能实现高维数据的非线性降维，广泛应用于数据分析及可视化分析^[11-13]。GTM能够尽可能地保持数据在高维空间上原有的拓扑相对关系，并且建立了高维与低维空间的双向映射。Viswanath等^[14-15]利用气动外形与气动力组成的高维数据作为GTM模型训练数据，进行了翼型的气动力优化设计，是GTM模型在气动设计领域的首次应用。

本文针对反设计方法目标压力分布难以给定、对设计经验依赖度高的问题，利用GTM模型的特点，建立了翼型设计变量与压力分布的映射模型，并结合全局优化算法实现了翼型压力分布高效反设计。

1 基于生成拓扑映射的气动反设计方法 1.1 生成拓扑映射模型

GTM能够实现高维数据到低维空间的映射，并保持原有的拓扑相对关系。设有D维数据集T={t₁, t₂, …, t_N}，但数据本质是L维的(L < D)。映射函数y(x, W)将L维空间的K个隐变量点x={x₁, x₂, …, x_K}映射到D维数据集T。图 1表示一个二维隐空间上的9个隐变量点，通过映射函数y(x, W)被映射到图右侧所示的三维空间上。

文献[10]给出了GTM模型详细的描述，采用中心位于y(x, W)且方差为β的高斯函数作为D维数据集T的分布

(1)

对x积分可得：

(2)

式中，p(x)为先验分布

(3)

表示均匀分布在隐变量空间内的K个点上。

由此得到

(4)

通过最大对数似然函数可求得参数W和β

(5)

分布p(t|W, β)为中心固定的带约束混合高斯模型。矩阵W是隐空间点向高维空间映射函数的权重参数矩阵，其维度为D×M, M为基函数个数。GTM算法可归纳为：已知隐变量空间中x的分布和数据空间的数据集合T，使用最大似然法求解参数W, β。常采用期望最大化(expectation maximization, EM)算法训练GTM模型参数^[16-17]。

1.2 基于GTM的气动反设计方法

基于GTM的气动反设计方法流程如图 2所示。

首先在根据气动外形参数化方法确定的设计变量确定设计空间。然后在设计空间内对设计变量进行抽样，得到样本集S，并计算对应外形的流场，这一步骤为实验设计。若设计变量个数为m，样本点数目为n，则S为n×m的矩阵。利用流场求解器得到样本点i对应的压力分布，若绕气动外形的压力分布点数为k，则得到维度为n×k的压力分布数据集C_p。将样本集与压力分布数据集合并，作为GTM模型的训练高维数据X_train。X_train是维度为n×(m+k)的矩阵，每一行由一个样本点及其对应的压力分布数据组成。

GTM模型的参数包括训练次数、隐空间维度L、隐空间网格点数K、基函数个数M。GTM模型训练过程一般在10步左右收敛^[18]。GTM模型常用于高维数据可视化，因此隐空间维度经常取为2。另外隐空间网格点数K随隐空间维度L呈指数增长，即K=N_mesh^L，N_mesh为一个维度上的网格划分数目。过多的隐空间网格点数会导致训练时间长或内存溢出等问题，一般L不大于4。基函数的个数一般取为N_mesh/2。

GTM模型训练完成后，算法建立了高维空间与隐空间的映射关系。隐空间的维度即为L维，给定隐空间的任意位置x，其中x是一个长度为L的向量，每个分量在±1范围内。利用GTM反映射关系可以得到对应高维空间的数据，即设计变量与压力分布组成的数据。因此，对于隐空间中的每一个点，可以通过GTM映射得到一个气动外形及其对应的压力分布。

定义表面压力分布与目标压力分布之间的差异为E_p

(6)

式中: N为表面压力分布点数; c_pⁱ为表面第i个点的压力值; c_pTⁱ为目标压力分布的第i个值。

本文采用遗传算法在隐空间进行寻优，目标函数为E_p。GTM 模型充当类似代理模型的作用。遗传算法的设计变量为隐空间变量，通过调用GTM模型，得到对应的压力分布和气动外形参数化变量。通过公式(6)可以得到隐空间变量对应的适应度值，供遗传算法执行交叉变异等操作。通过遗传算法得到最接近目标压力分布的气动外形参数化变量，进一步由参数化方法得到对应的气动外形。

隐空间维度较低，因此在种群数目及迭代步数较少的情况下即可高效找到全局最优解。最后通过CFD校验设计结果。若不满足设计目标，可以进行设计样本的更新及GTM模型的更新，或者调整目标压力分布，进行新一轮的设计。

2 模型参数影响分析

GTM模型参数对模型精度有重要影响，进而影响气动外形反设计结果，有必要对模型参数的影响规律进行研究。以NACA2412翼型为基本翼型，采用类函数/型函数参数化方法(class function/shape function transformation, CST)^[19]，翼型上下表面型函数阶数均取6，共14个设计变量，设计变量上下限取为基本翼型参数的±30%。计算状态为：α=2.2, Ma=0.2, Re=4.5×10⁶。采用拉丁超立方抽样方法(Latin hypercube sampling, LHS)取50个样本点，并利用Xfoil求解器得到翼型压力分布。GTM模型隐空间维度、隐空间网格点数、基函数个数取值如表 1所示。

在隐空间均匀取2 500个样本点，利用GTM映射回高维数据空间，得到翼型的设计变量与模型预测的压力分布。利用Xfoil求解器计算得到翼型的真实压力分布，并根据(5)式定义计算预测压力分布的误差，结果在表 1中列出。从表中可以看出，模型预测的误差最小值、均值、标准差等总体变化趋势为随隐空间维度的增加而减小，随隐空间网格点数与基函数个数的增加而减小。采用4维隐空间，20⁴个隐空间网格点、10个基函数训练GTM模型，选取预测误差最大的一组压力分布，如图 4所示。测试表明GTM模型能够建立高精度的映射关系，满足压力分布反设计的要求。综合考虑GTM精度、计算机内存、训练时间等因素，本文选取的隐空间维度为4，网格点数为20⁴，基函数个数为10。

3 翼型反设计算例 3.1 低雷诺数翼型压力反设计

本节根据1.2节的设计框架开展翼型压力反设计研究。基准翼型为NACA2412，设计状态为α=2.2, Ma=0.2, Re=4.5×10⁶，压力分布通过Xfoil求解得到。为验证方法的稳健性，给出不同的目标压力分布形态。遗传算法的种群数目取为20，迭代步数取为30步。在算法得到最优解之后，利用求解器进行验证。

图 4给出了基本翼型形状及其压力分布，以及输入的目标压力分布、求解器验证得到的设计压力分布，其中目标压力分布共给出了4种不同的形式。本文方法不需要给出绕翼型一周的完整压力分布，也不要求压力分布具有实际物理意义。给出的目标压力分布约在弦长40%~50%位置处。

图 4a~4b)给出的压力分布较为合理，设计压力分布与目标压力分布几乎重合，达到了设计预期。图 4c)~4d)给出的压力分布形态不合理，尤其是图 4d)“W”形的压力分布形式显然不具有物理意义，设计方法仍然能够找到与目标压力分布最为接近的结果，证明了设计方法的稳健性。总之，从翼型压力分布反设计的结果来看，设计得到的压力分布都很好地趋近于目标压力分布。

3.2 跨声速翼型压力反设计

本节以NACA0010翼型为初始翼型，进行压力分布反设计研究。设计状态为α=0, Ma=0.76, Re=7.7×10⁶。翼型压力分布利用PMB3D求解器得到^[20]，湍流模型为Spalart-Allmaras(SA)模型。计算网格如图 5所示，网格单元数为21 285，物面第一层网格数满足y⁺=1，远场距离翼型30倍弦长。

翼型采用CST参数化方法，设计变量数目为14。采用拉丁超方方法选取50个样本点，并计算对应的压力分布。给定的目标压力分布参考了超临界翼型的压力分布，翼型上表面压力分布平缓，翼型下表面前半部分压力分布保持一定前加载，其值约束在-0.3左右，靠近后缘部分保持翼型后加载。第1轮反设计结果如图 6所示，设计压力分布已经靠近目标压力分布，但仍有一定差距。将第1轮的设计变量作为初始值，进行第2轮设计。第2轮的设计结果如图 7所示，设计压力分布在翼型上表面、下表面前部基本与目标值重合，整体更加接近目标值。为更好逼近目标压力分布，尝试进行第3轮反设计，如图 8所示，设计压力分布整体上更好地趋近于目标压力分布，达到了理想设计结果。表 2给出了3轮设计翼型的气动特性，从表中看出，随着设计迭代，翼型升力系数不断提升，阻力系数不断减小，翼型升阻力特性提升明显。

由于目标翼型与基本翼型外形差异较大，本次设计进行了3轮设计，不断更新设计空间，取得了理想的设计效果。反设计问题中，无法确定事先给定的设计空间是否包含目标气动外形。本次设计也很好反映了实际设计过程的这一特点。即使设计空间未能包含目标压力分布对应的气动外形，本方法依然表现出很好的稳健性，通过不断调整设计空间，能够达到理想外形与压力分布。

3.3 三维层流短舱反设计

为了进一步验证本文方法适应性，选取客机发动机短舱模型为研究对象，以提高短舱外表面层流区域为目标进行压力分布反设计。选取的短舱模型长度5.54 m，最大直径3.719 m，模型左右对称，上下非对称。计算网格采用结构网格，网格量为508万。设计状态为马赫数Ma=0.85，雷诺数Re=3.1×10⁶，自由来流湍流度T_u=0.1%，攻角α=2°，侧滑角β=0°，湍流采用SST模型求解，转捩模型采用γ-Re_θ模型。采用FFD参数化方法单独对短舱外形进行参数化。选取位于最左侧的截面为设计剖面，如图 9所示，将图中所示红色圆点的沿短舱径向坐标变化量作为设计变量，共9个设计变量。

采用拉丁超方方法选取50个样本点，获取短舱流场信息，并提取短舱设计剖面外形和对应的压力分布作为训练样本。为提高短舱外表面层流区域范围，提高10%~30%弦长范围内的压力梯度，给定的目标压力分布如图 10所示。

图中的设计压力分布是经CFD验证的，结果表明反设计的压力分布很好地与目标压力分布吻合，达到了设计意图。设计后的剖面前缘半径增大，后段相对厚度有所减小。图 11给出了基本外形和反设计得到外形外表面对应的摩阻分布。从图中可知，层流区域延长了约弦长长度的4.4%。

图 12给出了基准短舱与设计短舱的摩阻云图，从图中可以明显看出经过设计的剖面处层流区域有明显延长。

4 结论

本文结合GTM机器学习方法与全局优化算法，发展了一种简单、鲁棒、高效的气动外形目标压力反设计方法。

1) GTM模型能够建立气动外形及其压力分布数据在低维隐空间的准确映射关系，借助GTM模型，能够准确预测压力分布对应的气动外形；

2) GTM模型将高维压力分布数据降低到低维隐空间，有利于全局优化算法高效寻优。另外设计目标定义灵活，不要求压力分布物理上存在，不要求给出绕气动外形一周的全部压力分布，极大提升了反设计方法的鲁棒性，使得设计问题定义更加灵活，使用更加方便。

3) 本文方法需要样本点较少，计算效率高，并且不局限于二维设计，在三维外形的设计中同样具有很好的适应性。在下一步研究中，本方法将用于复杂三维气动外形设计，进一步验证方法的实用性。