基于转置卷积神经网络的翼型结冰冰形图像化预测方法*

何 磊，钱炜祺，易 贤，王 强，张显才

(1. 中国空气动力研究与发展中心 计算空气动力研究所, 四川 绵阳 621000；2. 中国空气动力研究与发展中心 低速空气动力研究所, 四川 绵阳 621000)

结冰问题一直以来是影响航空飞行器飞行安全的重要隐患之一[1]。飞机结冰后，机翼、尾翼和舵面上的积冰破坏了物面附近流场，严重影响飞机气动性能和飞行性能，增加飞行风险，危害飞行安全。我国幅员辽阔，地势和气象条件复杂多变，先后发生了多起因结冰导致的坠机事故[2]。由于飞机结冰与气象条件密切相关，若能根据气象条件实时预测可能发生的结冰外形情况，不仅可以为飞行安全包线确定提供依据，为飞行员提供结冰信息和安全预警信息，也能支撑飞机的防除冰系统设计和智能结冰系统研制。因此，国内外进行了大量飞机结冰研究。

长期以来，开展结冰研究主要依靠结冰风洞试验[3]、飞行试验[4]和CFD模拟[5-6]。欧美一些发达国家纷纷建造结冰风洞并开发结冰数值模拟软件[7-8]，获取了大量飞机结冰数据[9]。这些数据被很多结冰研究作为验证或建模样本。国内结冰实验设备和数值模拟研究方面起步较晚，但中国空气动力研究与发展中心、北京航空航天大学、西北工业大学、南京航空航天大学等单位都在积极开展相关研究[10]。

随着机器学习算法的发展，国内外许多专家学者将机器学习方法应用到飞机结冰研究领域，包括自组织特性映射，径向基函数、BP神经网络、概率神经网络等，用于解决结冰严重性探测[11]、结冰位置和体积探测[12-13]、结冰冰形预测[14]、结冰后飞机气动特性影响[15-16]等问题。在结冰冰形预测方面，人工神经网络类算法应用较为广泛和成功，取得了许多有价值的成果。

结冰冰形预测研究内容主要包括对冰形曲线进行描述和建立预测模型两个方面。Ogretim等[17]最早提出一种基于神经网络的翼型结冰预测方法：首先需要进行坐标变换得到新的冰形曲线，将坐标变换后的冰形曲线看作一个非周期的复杂信号，展开为傅里叶级数形式，将冰形与飞行状态和气象参数之间映射关系的建模转化为了傅里叶系数与飞行状态和气象参数之间的建模；然后利用神经网络进行建模和预测，模型输入为来流速度、温度、液态水含量、平均水滴等效直径和结冰时间5个参数，输出为冰形曲线傅里叶级数的正余弦系数；该方法主要针对NACA00XX 系列翼型，在坐标变换算法方面有一定局限性。潘环等[18]也用类似方法对冰形预测的建模与方法进行了研究，并增加了相对湿度和攻角2项参数。Chang等[14]和李珺[19]改进了Ogretim等的方法，分别提出用小波包变换方法和多值变量拟合函数替代傅里叶变换对冰形曲线的描述。2种方法的冰形预测效果不仅取决于神经网络模型的设计和训练效果，也取决于对冰形的描述能力。另外，李珺仅仅采用了NASA格林结冰研究中心的80组风洞结冰数据，样本偏少，限制了神经网络的训练效果。鲍雨晨等[20]基于BP神经网络研究了通用飞机机翼的冰形预测方法，该方法对翼型特征点上的结冰厚度进行预测，再通过计算公式转换为冰形特征点坐标，对坐标点连线可得冰形轮廓；该研究为冰形预测提供了不错的实践途径，但研究中仅标记了14个机翼特征点，对冰形轮廓的描述能力稍显不足；另外研究中仅包含43组训练样本，而输入的实验变量包含空速、姿态角、总温等10个参数，训练样本数量明显偏少，训练获得的模型可能无法应对实验变量值变化稍大的情况。

近年来，深度学习作为机器学习的重要分支发展迅速，并继承了神经网络对非线性关系的描述能力，在众多领域得到了广泛应用。本文的主要研究工作就是基于深度学习技术构建翼型结冰冰形的预测模型，实现对冰形的图像化预测能力。

1 翼型结冰冰形预测问题分析

对翼型结冰冰形进行智能预测需要建立预测模型，其核心就是描述影响飞机结冰因素和冰形之间的映射关系，如图1所示。

图1 翼型结冰预测问题Fig.1 Problem of airfoil ice accretion prediction

影响飞机结冰冰形的参数主要有2类：大气环境参数和飞行状态参数。对于飞行状态参数，主要考虑飞行攻角α、飞行速度v；而大气环境参数，主要考虑液态水含量(Liquid Water Content, LWC)、水滴平均直径(Median Volumetric Diameter, MVD) 、环境温度T、结冰时长t。

以往采用的描述映射关系的方法，如Kriging、BP神经网络、RBF神经网络等，输出冰形曲线都是数值类型的，需要通过翼面特征点法向的冰形值去描述冰形曲线，或通过类似傅里叶变换的方法对冰形曲线做进一步处理，这些方法都要求翼面同一位置处法线方向冰形厚度只能存在单值，从而无法解决复杂冰形在翼面同一位置处法线方向冰形厚度存在多值的问题，如图2所示[15]。

图2 复杂冰形多值问题Fig.2 Multivalued problem of complex ice shape

为克服这一问题，提出采用图像方式对冰形进行描述，图像方式比传统数值方式更直观，描述能力更强，因此建模的关键就在于设计直接将图像作为输出的模型架构。

2 翼型结冰冰形图像化预测模型设计

根据上述对结冰预测问题的分析，提出如图3所示的建模和预测框架。主要包括数据获取、数据预处理、模型设计、模型训练、模型预测等环节，框架的特点是将图像化的冰形作为模型训练标签和预测的输出。

图3 建模与预测框架Fig.3 Framework of modeling and prediction

2.1 数据规范

翼型结冰冰形预测模型的输入输出数据格式是模型结构和接口设计的关键。对于输入参数，即影响飞机结冰的飞行状态和大气环境因素，其值存在符号和数量级的差异，因此需要通过归一化处理将其数值范围调整到一致的范围之内。对于输出的冰形图像，考虑将其灰度化，既能满足对冰形几何特征描述的需求，又可避免多通道带来的额外计算量。冰形图像尺寸设置为512×256，每个像素点的取值在0～255之间，实际运算中将其归一化到0～1之间。

如图4所示，为避免由于图像中冰形的坐标轴取值范围不一致对模型带来的计算误差，冰形灰度图像宽度对应的x坐标轴取值范围固定为[-0.2, 1.0]，图像高度对应的y坐标轴取值范围固定为[-0.2, 0.2]，坐标轴取值范围可根据实际情况调整。

图4 结冰翼型图像坐标轴范围固定Fig.4 Fixed the coordinate range of iced airfoil image

2.2 模型结构

卷积神经网络(Convolutional Neural Network，CNN)是深度神经网络的典型结构，卷积操作为卷积神经网络提供了强大的图像特征提取能力。转置卷积(又名反卷积)是卷积操作的相反过程，可以对卷积操作提取的编码器中的特征进行解码。利用多个转置卷积操作可以实现生成图片的目的，深度卷积对抗生成网络的生成器就是使用多个转置卷积层对随机噪声值进行操作以实现生成完整的图片功能。

因此基于转置卷积神经网络设计了如图5所示冰形预测模型结构。输入的结冰条件参数通过2个全连接网络层映射到更大尺寸的数据结构，2个全连接网络层均使用ReLU激活函数。接着通过5个转置卷积网络层将低维特征向量向高维特征空间映射。除了预测模型的输出层外，其他的转置卷积网络层后都连接了批标准化(Batch Normalization， BN)网络层[21]和ReLU激活函数层，批标准化技术通过规范化手段，使得每一层神经网络的输入保持相同分布，以保证梯度传播到每一层，避免出现梯度消失现象。从结冰条件参数到冰形图像是一个典型的非线性映射问题，ReLU激活函数层的功能就是为神经网络提供非线性映射功能，从而提升冰形预测模型的预测能力，并使用Dropout网络层随机丢弃神经网络单元，提高泛化能力；输出层在转置卷积层后连接一个sigmoid激活函数层。

图5 翼型冰形预测模型网络结构Fig.5 Network structure of airfoil ice shape prediction model

模型采用二元交叉熵(binary cross entropy)损失函数，为了防止模型出现过拟合，进一步提升泛化性能，损失函数中加入L2正则化惩罚项，损失函数如式(1)所示。式中，右边第一项为交叉熵损失项，右边第二项为L2正则化项，它表示了模型的复杂度。

(1)

式中：yi表示类别i的真实标签；pi表示模型计算出类别i的概率值；N表示训练样本总数；k表示类别数，二元交叉熵中k=2；w表示神经网络的权重；λ表示L2正则化率；n表示整个网络的神经元总数。

3 实验数据生成

考虑以某运输机尾翼翼型(NACA0012)为实验对象，通过数值模拟方法生成实验所需冰形数据。

3.1 结冰参数范围确定

对于飞行参数取值范围，考虑固定飞行攻角条件下飞行速度的取值，所研究的运输机最小飞行速度约为268 km/h(74.44 m/s)、巡航飞行速度为550 km/h(152.77 m/s)。

对于气象参数取值范围，有研究表明：飞机最容易发生结冰的温度范围为0～-20 ℃，特别是在-2～-10 ℃ 范围内遭遇结冰的次数通常最多，而-2～-8 ℃是飞机发生强烈结冰的主要温度范围[22]。结冰中常见的过冷水滴平均直径在20～40 μm之间[19]，如《中国民用航空规章第 25 部运输类飞机适航标准》(CCAR 25部)附录C中通常考虑的过冷水滴尺寸就在15～40 μm范围。液态水含量是影响结冰最重要的因素之一，影响云中液态水含量的因素较多，随水汽凝结和降水而变化，云层中不同位置液态水含量也不相同。根据CCAR 25部附录C中给出的大气约束条件，通常液态水含量的考虑范围为0.2～0.8 g/m3。

综上分析并结合具体研究需求，确定的结冰参数取值范围如表1所示。

表1 结冰参数取值范围

3.2 数值模拟方法验证

本节对典型翼型结冰进行数值模拟，并将冰形计算结果与风洞试验结果进行对比，以验证数值模拟方法的可靠性。风洞试验数据取自美国NASA格林结冰研究中心和 FAA 威廉·J.休斯技术中心联合发起的现代翼型项目的翼型结冰资料[15]，该资料记录了3种典型飞机翼型在多种结冰条件下的冰形和气动特性风洞试验数据。

3.2.1 数值计算

翼型采用格林结冰研究中心的现代翼型项目中记录的商用喷气飞机(business jet)翼型[24-25]，如图6所示，该翼型是现代商用飞机的典型翼型。

图6 商用飞机翼型Fig.6 Airfoil of business jet

共安排2组实验，结冰条件参数包括液态水含量、水滴平均直径、环境温度、结冰时长、飞行速度、飞行攻角。2组实验结冰条件参数取值如表2所示，分别对应结冰资料中在结冰研究风洞的试验序号“run213”和“run214”[15]，除结冰时长不同外，其他条件一致。

表2 验证算例计算条件

翼型结冰数值计算采用中国空气动力研究与发展中心的计算方法和计算软件irc2d，主要包括流场计算、过冷水滴运动计算、结冰计算和物面外形更新4个步骤[16, 26]。如图7所示，计算采用C型网格，在网格生成阶段已经加入攻角。

图7 结冰计算网格(α=6.2°)Fig.7 Computational grid of ice accretion(α=6.2°)

3.2.2 结果验证

表3是风洞试验的冰形几何特征参数值，图8分别给出了2组实验数值计算冰形和风洞试验冰形的对比情况。对于冰体轮廓，两组实验的数值计算结果与风洞试验结果基本一致，冰角厚度、角度也符合较好；但对于冰体在翼型物面的结冰上极限和下极限位置，计算结果和试验结果对比稍有差异，两组实验中计算冰形上极限稍大，下极限稍小；实验1的驻点结冰厚度符合较好，实验2的驻点结冰厚度略有差异。总体而言，计算冰形和试验冰形在冰体轮廓、冰形体积和主要特征方面符合较好，说明所采用的翼型结冰数值模拟方法是可靠的。

表3 冰形几何特征参数值

(a) 实验1(a) Experiment 1

(b) 实验2(b) Experiment 2图8 计算结果与风洞试验结果对比Fig.8 Comparison of computational and wind tunnel test results

3.3 数据生成

使用上节所述irc2d翼型结冰数值仿真软件计算NACA0012翼型在α=2°的情况下的结冰外形。计算采用C型网格，如图9所示。

图9 结冰计算网格(α=2°)Fig.9 Computational grid of ice accretion(α=2°)

计算共获得11 200组训练样本数据，对应飞行速度为70 m/s, 80 m/s, 90 m/s, 100 m/s, 110 m/s, 125 m/s, 140 m/s, 150 m/s；温度为-20 ℃, -17 ℃, -14 ℃, -11 ℃, -8 ℃, -6 ℃, -4 ℃, -2 ℃；水滴直径为20 μm, 30 μm, 40 μm, 50 μm, 60 μm；液态水含量为0.2 g/m3, 0.35 g/m3, 0.5 g/m3, 0.65 g/m3, 0.8 g/m3；结冰时长为6 min, 9 min, 12 min, 15 min, 18 min, 20 min, 22.5 min的全组合。

共获得768组验证样本数据，对应飞行速度为85 m/s, 98 m/s, 120 m/s, 145 m/s；温度为-18 ℃, -10 ℃, -5 ℃, -3 ℃；水滴直径为23 μm, 32 μm, 52 μm；液态水含量为0.3 g/m3, 0.4 g/m3, 0.6 g/m3, 0.7g/m3；结冰时长为7 min, 10 min, 12.5 min, 18.5 min的全组合。

获得冰形样本数据后，按照建模框架所述数据规范将训练样本集转化为如图10所示的灰度图像。

图10 翼型结冰冰形灰度图像Fig.10 Gray scaled images of iced airfoil shape

4 仿真与结果

根据图5所示神经网络模型建立冰形预测模型，并使用上节生成的11 200组训练样本对网络模型进行训练，768组验证样本对训练好的网络模型进行测试。由于实验数据的飞行攻角固定值为2°，可以只考虑其他5个影响因素，因此网络输入层尺寸为5。

4.1 预测模型训练

论文基于开源深度学习架构tensorflow的keras高级接口实现所提的翼型结冰预测模型的构建、训练和预测。预测模型训练参数设置如下：训练优化方法选择Adam[27]，其中学习率(learning rate)的初始值为0.001，参数β1设置为0.9，β2设置为0.999，为稳定训练过程使用学习率衰减，设置值为1.0×10-8；分批大小为20，即每批输入20组结冰条件参数数据和冰形图像标签进行训练；迭代次数为200。

模型采用CPU训练模式，用于建模和训练的计算机配置为：Intel Core i7-7700、3.6 GHz、4核8线程CPU、16 GB内存。基于以上训练参数，对11 200个样本的训练耗时约为88.2 h。训练过程中，样本的损失函数值随迭代次数的变化曲线如图11所示。从图11中可以看出，模型在约前10次迭代中，损失函数值迅速趋于收敛，之后随着迭代次数的增加，有继续小幅下降的趋势，最终稳定在0.007 8附近。

图11 训练过程Fig.11 Training history

4.2 结果分析

模型训练完成后，将训练样本集和验证样本集的结冰条件参数输入预测模型，预测其对应的翼型结冰冰形图像，测试模型的预测能力，检验预测模型的训练和泛化效果。训练样本集预测耗时576 s(单个样本耗时约0.048 s)，验证样本集预测耗时38 s(单个样本预测耗时约0.049 s)。

模型预测输出结果为类似图10的灰度图像，为便于观察和对比分析，利用图像技术对结果进行了加入原始翼型物面、颜色替换、剪裁等后处理操作。

图12展示了3组典型训练冰形的模型预测结果，以及与CFD计算冰形进行对比的情况。从图12中可见，在冰体轮廓、结冰上下极限位置、冰角厚度和角度等主要几何特征方面，3组实验预测结果，以及与CFD计算结果都基本重合，符合较好，说明模型在训练样本集上表现良好。

(a) α=2°, v=70 m/s, T=-20 ℃, MVD=30 μm,LWC=0.5 g/m3, t=20 min

(b) α=2°, v=80 m/s, T=-20 ℃, MVD=20 μm,LWC=0.8 g/m3, t=18 min

(c) α=2°, v=140 m/s, T=-17 ℃, MVD=60 μm,LWC=0.8 g/m3, t=15 min图12 典型训练冰形预测结果Fig.12 Prediction of typical training ice shape

图13展示了4组典型测试冰形的模型预测结果，以及与CFD计算冰形进行对比的情况。对于冰体轮廓，4组实验的模型预测结果与CFD计算结果基本重合，符合较好。冰形冰角角度、结冰上下极限位置等几何特征参数也都符合较好。对于结冰厚度这一特征，部分实验模型预测结果与CFD计算结果略有差异，且这种差异随着冰形厚度的增加而变大。但总体而言，模型预测冰形和CFD计算冰形在冰体轮廓、冰形体积和主要特征方面符合较好。可见，预测模型的整体泛化性能较好，但对结冰较厚情况的泛化性能仍有进一步提升空间。

(a) α=2°, v=85 m/s, T=-10 ℃, MVD=35 μm,LWC=0.7 g/m3, t=12.5 min

(b) α=2°, v=98 m/s, T=-18 ℃, MVD=35 μm,LWC=0.4 g/m3, t=18.5 min

(c) α=2°, v=145 m/s, T=-18 ℃, MVD=52 μm,LWC=0.6 g/m3, t=10 min

(d) α=2°, v=145 m/s, T=-18 ℃, MVD=23 μm,LWC=0.7 g/m3, t=12.5 min 图13 典型测试冰形预测结果Fig.13 Prediction of typical testing ice shape

5 结论

本文针对飞机结冰问题，开展了翼型结冰冰形预测方法研究，提出了建模和预测框架，设计了基于深度神经网络的图像化冰形预测模型，用于建立冰形与飞行状态参数、气象参数之间的映射关系，主要考虑了飞行速度、攻角、大气中液态水含量、水滴平均直径、结冰时的温度、结冰时长等多物理参数对冰形的影响。

以某运输机水平尾翼(NACA0012翼型)为对象，利用CFD数值模拟生成的冰形作为训练和验证样本，对所建立的冰形预测模型进行了训练和测试。结果表明：

1)提出的翼型结冰冰形图像化预测方法是可行的，预测冰形与CFD数值计算的冰形在冰形轮廓、结冰上下极限、上下冰角位置、结冰厚度等主要几何特征参数方面都符合较好，但在结冰较厚情况下，模型泛化性能还可进一步提高。

2)虽然模型训练耗时较多，但预测模型训练完成后，便能够快速预测一定范围内的冰形，且计算速度快(单工况计算耗时约50 ms)，能适应机载要求；另外在模型训练过程中,若使用图形处理器计算能大幅减少训练时间。

总之，深度学习方法在飞机结冰研究领域具有很强的应用前景，下一步考虑将风洞试验数据加入样本集对预测模型进行训练，增强模型的工程实用性。同时，在现有研究基础上，继续基于深度学习方法研究飞机结冰后对气动特性的影响。