0 引　言

翼型设计是飞机机翼、涡扇发动机叶片、风力机叶片、赛车稳定尾翼等设计的基础。针对航空航天、能源动力、交通运输等不同领域的应用，目前工程上已发展了NACA系列、RAE系列、S系列、DU系列、NPU-WA系列等翼型族^[1-3]。如何从这些翼型族中有效地总结出物理规律，形成知识，以更好地指导工程应用及翼型族的不断优化完善，在理论和工程上都有重要的研究价值。近年来，数据挖掘和数据库知识发现（Data Mining & Knowledge Discovery in Database，简称DM&KDD）技术得到了快速发展，该技术融合交叉了统计学、模式识别、人工智能、符号推理、机器学习、数据库以及高性能计算等技术，包括特征化与区分、聚类、分类与回归、关联分析及异常点检测等挖掘模式，在经济、商业、金融、天文等行业得到了成功应用，尤其是在故障诊断、生产优化、丰富知识库、决策支持等方面展现出强大优势，在国际上掀起了研究热潮 ^[4-5]。该技术在空气动力学研究中也已得到了应用，如飞行器气动布局设计知识挖掘^[6]、流场结构特征提取^[7-8]、基于飞行试验数据的气动建模和气动参数辨识结果可信度确认^[9]以及翼型气动系数预测等^[10]方面。文献[11-13]将该技术用于翼型气动设计参数空间的挖掘分析，采用显著变量识别、总变差分析（ANOVA）、自组织神经网络等数据挖掘技术对设计空间进行知识挖掘，但在一定的局限性，一是仅针对单个翼型或是一类翼型进行分析，所得到的研究结论适用范围相对较窄；二是主要进行设计变量的参数影响分析，对翼型参数本身内蕴的一些规律挖掘较少。本文从工程翼型库出发，基于CST参数化方法实现了数据异常翼型的检测，揭示了不同翼型CST参数内蕴的一些规律；采用级差分析、SOM自组织映射、Apriori算法等知识挖掘方法分析了CST参数对典型工况翼型气动特性的影响规律，并建立了预测模型。

1 工程翼型库翼型数据的清洗与参数化

选取伊利诺伊大学香槟分校UIUC的工程翼型库^[14-15]为研究对象。该翼型库中初始包含1550组常用翼型数据。但其中存在部分翼型重复、部分翼型数据异常的情况，因此，在对翼型库中翼型进行分析之前需要对这些翼型数据进行数据清洗。

首先通过几何数据直接对比的方法去掉了重复和数据不完整的翼型；然后对翼型进行重建，一方面通过比较重建前后翼型差异的方法可以识别局部数据点异常和数据过于稀疏的翼型，另一方面，基于不同翼型重建参数的对比也可以进一步识别出翼型库中的重复翼型。

这里采用的翼型重建方法，也称参数化方法，为美国波音公司Kulfan等提出的基于型函数/类函数变换的参数化方法（Class Function/Shape Function Transformation，CST），该方法参数具有明确的几何意义、控制参数少、适应性强、精度较好^[16]。CST方法中，曲线表示为：

${y_u} = C_{{N_2}}^{{N_1}}(x)\;{S_u}(x) + x\;{y_{_{TEu}}} $

(1)

${y_l} = C_{{N_2}}^{{N_1}}(x)\;{S_l}(x) + x\,{y_{_{TEl}}} $

(2)

${S_u}(x) = \sum\limits_{i = 0}^{{N_u}} {{A_{ui}}\frac{{{n_u}!}}{{i!({n_u} - i)!}}{x^i}{{(1 - x)}^{{n_u} - i}}} $

(3)

${S_l}(x) = \sum\limits_{j = 0}^{{N_l}} {{A_{lj}}\frac{{{n_l}!}}{{j!\,\,({n_l} - j)!}}{x^j}{{(1 - x)}^{{n_l} - j}}} $

(4)

其中：下标u、l分别代表翼型上、下表面；x为横坐标；y_u、y_l为上下表面纵坐标；y_TEu、y_TEl为后缘点纵坐标； $C_{{N_2}}^{{N_1}}(x) = {x^{{N_1}}}{(1 - x)^{{N_2}}}$ 为类别函数；S(x)为形状函数，采用Bernstein多项式形式；n_u、n_l为上下翼面Bernstein多项式的次数。对一般翼型而言，N₁ = 0.5，N₂ = 1.0，翼型即可用Bernstein多项式的系数向量：A_ui，( i = 1,···,n_u)、A_lj，(j = 1,···,n_l)来描述，这些系数称为翼型设计参数，通过拟合逼近翼型几何坐标来得到。

图1给出了翼型库中三种翼型“dga1182”、“saratov”和“fx84w150”重建前后的结果对比，其中“Original”表示原始翼型，“CST”表示采用CST方法的重建结果。可以看到，对前两种翼型，由于原始翼型本身数据存在异常，使得重建结果和原始翼型有较大差异，而对于“fx84w150”翼型，原始数据完整光滑，重建前后符合较好。因此，采用这一方法来进行翼型数据异常的识别是可行的。

采用上述方法对原始翼型数据进行数据清洗后，最终可得到1499组可用的翼型。同时，对这1499组翼型，由于已经用CST方法进行了参数化重建，所以每个翼型都可用上下翼面的A_ui、A_li（i = 1,···,n）来描述，这里取n = 7，共14个参数。下面先对这些参数蕴含的一些规律进行挖掘。

2 翼型CST参数规律挖掘 2.1 CST参数取值分布规律分析

针对上述1499组CST参数，首先对各参数的取值范围进行统计分析，得到各参数取值的概率密度分布函数如图2。可以看到，各参数的取值分布近似为正态分布形式，其原因可用中心极限定理来进行定性解释^[17-18]，即将翼型几何形状视为随机变量，不同领域专家针对各自应用场景来总结提炼出不同的翼型几何形状，这些工作同时也具有较强的独立性，使得翼型几何参数在一定程度上具备了随机变量的多样性和独立性，同时，针对一族翼型的不断完善也可视为大量微小因素作用的叠加，从而使得翼型库中的翼型几何参数近似满足中心极限定理的成立条件，A_li和A_ui近似满足正态分布。

针对各CST参数都取概率密度最大的值，生成一个翼型，如图3中的“Mean airfoil”示，在翼型库中与其欧式距离最接近的是“sd7080”翼型，二者的比较示于图3。这一翼型并不是气动性能最好的翼型，但应该是比较“中庸”的翼型，可适应大多数的工程应用需求，但这一推测尚需后续工程实践检验。

2.2 CST参数之间关联规律挖掘

进一步采用开源数据挖掘平台WEKA（Waikato Environment for Knowledge Analysis）^[19]对各CST参数之间的关联关系进行挖掘，该软件由新西兰Waikato大学开发，集成了大量数据预处理、关联规则挖掘、分类、聚类等算法，目前在科学研究和工业领域都得到了较成功的应用。该软件给出的A_ui（i = 1,···,7）之间、A_li（i = 1,···,7）之间的关联关系分析结果如图4所示，从图中可以看到，A_l2（简记为“L2”，下同）与A_l4、A_l3与A_l5、A_l4与A_l6、A_u2与A_u4、A_u3与A_u5、A_u4与A_u6之间存在一定的弱相关性，计算出相关系数约为0.8。软件同时也计算了上下翼面参数A_ui（i = 1,···,7）与A_li（i = 1,···,7）之间的相关系数，值较低，在图中略去。

2.3 CST参数聚类分析

接下来针对翼型的CST参数对其进行聚类分析。对于对称翼型，其CST参数满足A_li = − A_ui, i = 1,···,7，因此可首先通过分析CST参数识别出翼型库中的对称翼型，共108个翼型；后面只需对剩余的1391组翼型进行聚类。

为了克服翼型厚度对聚类结果的影响，将这1391组翼型的上翼面最大厚度取为相同值，对y坐标进行规范化，然后针对变换后的翼型重新计算CST参数用于聚类。采用WEKA软件中的Kmeans聚类算法，选取聚类数K = (2,3,···,12)，得到聚类误差E（所有点到其聚类中心的距离和）随K变化曲线如图5。计算曲线的导数，选取导数变化趋于平缓的点（即K = 5），作为最终聚类数。图6给出了此时聚类中心对应的翼型形状，可以看到，不同聚类中，上翼面曲线整体形状差异不大，仅峰值点位置有所差异，但下翼面曲线则分为不同类别：

1）曲线整体为凹曲线，如图6中的“Group 1”；

2）曲线在前缘附近为凸曲线，后段近似为直线，如图6中的“Group 2”；

3）曲线存在拐点，前段为凸曲线，后段为凹曲线，但尾部没有出现明显向内反凹的曲线，如图6中的“Group 3”；

4）曲线存在拐点，前段为凸曲线，后段为凹曲线，但尾部出现了向内反凹的曲线，如图6中的“Group 4”、“Group 5”，二者整体形状一致，但凸曲线、内凹曲线的峰值点位置与幅值差异较显著。

这一聚类结果与工程上将翼型分为“平凸翼型”、“双凸翼型”、“凹凸翼型”等类别^[20]是基本一致的。

3 翼型CST参数与气动特性关联规律挖掘

翼型CST参数与气动特性之间的关联规律挖掘是翼型气动知识挖掘的重要内容。本文采用两类方法来分析UIUC工程翼型库中翼型CST参数与气动特性的关联规律。一类是因果逻辑的思路，采用基于试验设计的极差分析方法；另一种是数据约简归纳的思路，采用自组织映射方法（SOM）和基于Aprior算法的关联关系挖掘方法。翼型的气动特性以马赫数0.2、雷诺数2×10⁶、0°～25°攻角下的气动特性为例，重点关注最大升力系数C_L,max、最大升力系数对应攻角α_CL,max、升力线斜率C_Lα、零升阻力系数C_D0、最大升阻比K_max、最大升阻比对应攻角α_K,max等6个气动特性指标参数。采用本单位自主研发的流体计算软件MBNS2D^[21]，计算求解雷诺平均N-S方程，湍流模型为k-ω SST两方程湍流模型。计算时间推进采用LU-SGS方法，对流项采用Roe通量差分格式，扩散项和源项离散采用中心差分格式。该计算软件已通过了多个算例的考核，广泛应用于工程实践。针对本文中的翼型，计算网格为373×121的C型结构网格（如图7），翼型升阻力系数典型状态计算结果示于图7右下角。

3.1 基于均匀试验设计的极差分析方法

在参数影响分析过程中采用14参数3水平51样本的均匀试验设计方法^[22]，根据图2中每个参数的取值范围，同时为了避免生成的翼型出现上下翼面交叉的情况，选取水平值为：A_l1的水平值为（−0.3, −0.18, −0.06），其余A_li（i = 2,···,7）的水平值为（−0.18, −0.06, 0.06）；所有A_ui（i = 1,···,7）的水平值为（0.1, 0.22, 0.34）。参照文献[23]的U₅₁(3¹⁴)均匀设计表，生成51个样本的计算网格进行流场和气动特征参数计算，然后进行极差分析^[24]。图8中分别给出了不同CST参数对应的三个典型气动参数（最大升力系数C_L,max、升力线斜率C_Lα、零升阻力系数C_D0）计算结果极大值和极小值，图中横轴对应不同参数，纵轴对应气动参数的最大值（“Max”）和最小值（“Min”），二者之差即为该CST参数的极差，其值越大，表示气动参数受该参数的影响越大。表1给出了具体的极差计算结果，从结果中可以看到：当前工况下，对最大升力系数C_L,max、最大升力系数对应攻角α_CL,max、升力线斜率C_Lα、零升阻力系数C_D0、最大升阻比K_max等5个气动特性指标参数影响最大的CST参数分别是A_u1、A_u1、A_l7、A_u1、A_l7，对最大升阻比对应攻角α_K,max影响最大的CST参数是A_l1和A_l7。

3.2 自组织映射方法（SOM）

3.1节对试验设计的51组样本进行了马赫数0.2、雷诺数2×10⁶、0°～25°工况的气动特性计算。本小节对UIUC翼型库清洗后的1499组翼型进行了计算，并对结果进行自组织映射分析。自组织映射分析是人工神经网络中的一种无监督学习算法，在保留原数据特征的前提下将高维数据有序映射到SOM成分图中，各成分图相同坐标位置对应同一SOM神经元，通过比较各图中相同区域颜色模块（即云图，图9）的相似程度，可以直观地分析变量间相互作用关系，该方法不仅可以探索自变量对函数的影响规律，还能够探索相同自变量对应的多个函数之间的相互影响关系^[6,11]。

对于SOM成分图中相同区域颜色模块相似程度的分析，通常依靠的是主观判断，具有较大的不确定性。本节对其进行了定量处理，将各成分二维云图值通过矩阵拉直的方法进行向量化，然后再计算所得向量之间的相关系数。n维向量x和y的相关系数r(x,y)可以定义为：

$ r({\boldsymbol{x}},{\boldsymbol{y}}) = \frac{{{s_{xy}}}}{{{s_x}\,{s_y}}} $

(5)

其中， ${s_{xy}} = \dfrac{1}{{n - 1}}\displaystyle\sum\limits_{i = 1}^n {\left[ {({x_i} - \bar x)({y_i} - \bar y)} \right]}$ ， ${s_x} = \sqrt {\dfrac{1}{{n - 1}}\displaystyle\sum\limits_{i = 1}^n {{{({x_i} - \bar x)}^2}} }$ ， ${s_y} = \sqrt {\dfrac{1}{{n - 1}}\displaystyle\sum\limits_{i = 1}^n {{{({y_i} - \bar y)}^2}} }$ ， $\bar x = \dfrac{1}{n}\displaystyle\sum\limits_{i = 1}^n {{x_i}}$ ， $\bar y = \dfrac{1}{n}\displaystyle\sum\limits_{i = 1}^n {{y_i}}$ ， ${x_i}$ 和 ${y_i}$ 分别为向量x和y中的各个元素。

从而可以通过相关系数的大小来对云图的相似程度进行定量化的判断。基于这一方法，从图中发现：

1）自变量之间：L2与L4、L4与L6、U3与U5之间的相关系数超过0.95，相应云图之间有较强的正相关，这与图4中的结果是一致的。

2）气动特性参数之间：C_L,max与K_max正相关；α_CL,max与C_Lα、C_D0与C_Lα有一定负相关。

3）自变量与气动特性参数之间：C_L,max与L7、K_max与L7正相关；C_Lα与U5负相关；C_D0和U2有一定的正相关；C_Lα与U3有一定的负相关。其中部分结果与表1结果一致，部分在表1中体现不显著，这与表1中的参数取值范围有限，而不同参数区间影响规律不同有关。

3.3 基于Aprior算法的关联关系挖掘

本小节利用WEKA中的Apriori算法对上述1499组翼型数据集进行关联关系的挖掘。Apriori算法利用频繁项集所有非空子集是频繁项集这一性质，通过连接和剪枝可以高效地找出数据集中的频繁项集，从而发现相应的关联规则，关联规则是“如果X，那么Y”的形式。要得到有用的规则，就需要计算支持度s和置信度γ。

$ s\left( {X \to Y} \right) = P\left( {X \cup Y} \right) $

(6)

$ \gamma \left( {X \to Y} \right) = P\left( {Y|X} \right) $

(7)

支持度表示规则出现的概率，而置信度表示规则正确的概率。关联规则的支持度以及置信度越大，该规则实用性越高。

对1499组翼型数据集分别对CST参数与气动特性参数之间的关联关系进行挖掘，主要结果如表2所示。从表中可以看到，可以给出有较大影响的参数，同时还能够给出多个参数与气动特性之间的关联关系。以C_D0为例，从关联规则可以看出参数U1、U2与C_D0关联较为紧密，与表1中的影响规律基本符合。同时，通过挖掘出的关联规则能得出多个参数与气动特性参数取值范围的关联关系，这是其他方法不具备的优势。

综上，三种关联关系分析方法之间各有优缺点。极差分析方法能给出各CST参数影响的显著程度，但对数据要求较为严格，需在分析前对样本进行试验设计；SOM方法和Aprior算法对数据样本没有严格要求，但只能识别出关联关系，不能给出显著程度。SOM方法的特色在于能较为直观地给出自变量之间、因变量之间、以及自变量与因变量之间的影响关系，Aprior算法则能进一步给出多个自变量与因变量之间的关联关系。

4 翼型气动特性建模预测

基于上面的影响因素分析和1499组翼型典型工况的气动特性计算结果，本节建立该典型工况下的翼型气动特性预测模型，即输入为翼型的14个CST参数，输出翼型典型工况下的最大升力系数、升力线斜率、零升阻力及最大升阻比。建模时将1499组翼型分为1300组翼型组成的训练集和199组翼型组成的测试集，在训练集上进行模型训练，然后对测试集中的翼型进行预测，通过计算预测值与标签值（CFD计算结果）之间的平均相对误差来衡量所构建模型的预测精度。平均相对误差定义为：

$ {\varepsilon _{{\rm{avg}}}} = \frac{{\displaystyle\sum\limits_{i = 1}^N {\left| {y_p^i - y_l^i} \right|} }}{{\displaystyle\sum\limits_{i = 1}^N {\left| {y_l^i} \right|} }} \times 100 $

(8)

其中，N为样本数， $ y_p^i $ 为预测值， $ y_l^i $ 为标签值。

采用两种模型架构，一种是支持向量机（SVM）模型，建模算法采用的是WEKA软件中的SMO回归算法；另一种是基于Tensorflow的深度神经网络模型，包含14个节点的输入层、10层16个节点的隐藏层以及4个节点的输出层，隐藏层的节点采用ReLU激活函数，选择预测值与标签值的均方误差作为网络损失函数，采用结合Adam方法的小批量梯度下降法作为优化训练方法，分批大小设为8，Adam初始学习率设为0.001，迭代轮数设为1000。

图10给出了训练集上深度神经网络模型的收敛过程。以最大升力系数和零升阻力为例，利用训练完的网络在训练集和测试集上的预测结果如图11以及表3所示，其中也给出了与SVM训练和预测结果的对比。从图中可以看到，相比于SVM在测试集上的预测结果，深度网络的预测值与相应的测试标签值的分布更接近于45°斜对角线，表明深度网络比SVM具有更好的预测能力。从表3也可以看出，无论是在训练集还是在测试集上，深度神经网络预测值和标签值之间的平均相对误差远小于SVM的平均相对误差。

从上述建模结果可以看到，对于工程翼型库，采用深度神经网络模型能对其气动特性进行较为有效的建模，建模精度和泛化能力优于SVM模型。同时，所建立的模型蕴含了工程翼型库的信息，基于此模型进行翼型的反设计可以得到更具有工程可实现性的设计结果。

5 小结与展望

本文以伊利诺伊大学香槟分校UIUC的工程翼型库为研究对象，在对翼型进行CST参数化的基础上对CST参数内蕴规律、CST参数与典型工况气动特性的关联规律进行了挖掘，同时也对关联挖掘方法的各自特点进行了探讨，最后建立了气动特性与CST参数之间的深度神经网络模型。主要结论如下：

1）通过翼型CST参数化后重建结果与原始翼型的比较，可以识别存在异常的翼型数据。

2）UIUC工程翼型库中的翼型CST参数取值分布近似为正态分布形式；翼型库CST参数A_l2与A_l4、A_l3与A_l5、A_l4与A_l6、A_u2与A_u4、A_u3与A_u5、A_u4与A_u6之间存在一定的弱相关性；CST参数聚类分析的结果与工程上“平凸翼型”、“双凸翼型”、“凹凸翼型”等分类结果基本一致。

3）可采用极差分析、SOM自组织映射、Apriori算法等方法挖掘CST参数与典型工况气动特性之间的关联关系。极差分析能给出各CST参数影响的显著程度，但在分析前需对样本进行试验设计；SOM和Apriori算法对数据样本没有严格要求，但只能识别出关联关系，不能给出显著程度。此外，在SOM分析中，可以通过分析矩阵拉直向量的相关系数大小来对成分图的相似程度进行定量判断。

4）采用深度神经网络建立的CST参数与典型工况气动特性之间的预测模型在拟合能力和泛化能力上优于SVM模型，同时，神经网络模型中蕴含了目前工程翼型库中的信息，基于此模型进行翼型的反设计可以得到更具有工程可实现性的设计结果。

下步工作主要考虑在以下方面深入开展：

1）目前UIUC翼型库中的翼型均为比较成熟的翼型，近年来工程上发展的一些新型翼型没有包含在内。下步将对翼型库不断扩充，完善相应的挖掘分析方法，不断丰富和深化对翼型性能的认识。

2）本文只针对工程翼型进行了低速典型工况下的气动特性分析，重在建立方法。后续需对翼型先按大致用途分类，如分为厚翼型和薄翼型，或者分为高升力翼型和超临界翼型等，然后再重点针对低速或亚跨超气动特性进行规律挖掘分析。

3）用户使用数据挖掘软件不可能一劳永逸地从数据中获取知识。知识挖掘过程需要专家经验、数据与知识挖掘工具的深度结合，有时甚至需要针对特定的挖掘需求对数据和挖掘工具进行扩充。因此，如何在知识挖掘中融入专家经验、背景知识、约束、规则，以及如何有效理解数据挖掘的结果，都是后续需要深入研究的问题。

4）后续可在本文所发展的数据挖掘结果的基础上进一步提炼翼型参数与气动特性之间的因果关系，以揭示翼型参数对气动特性的影响规律，更好地指导和改进翼型设计。