基于代理模型全局优化的自适应参数化方法

2020-11-05张伟高正红周琳夏露

航空学报 2020年10期

张伟，高正红，周琳，夏露

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

飞翼布局[1-4]因其具有良好的气动效率和隐身特性，被新一代军用飞机采用。相较于常规布局，由于尾翼等部件的缺失以及翼身融合的特征，翼型的作用对于飞翼布局更加重要。在气动特性方面，相较于传统翼型设计主要考虑升阻特性要求，飞翼气动布局则要求在翼型设计中增加量化的俯仰力矩以满足飞机纵向力矩配平的需要。在隐身特性方面，翼型的隐身设计是飞翼布局雷达散射截面(Radar Cross Section,RCS)减缩的重要手段之一[5]。因此，飞翼布局翼型设计需要在满足气动、隐身及控制等多学科性能要求下进行。

对于翼型设计而言，设计变量(Design Variable，DV)的选择直接影响设计结果[6]。在开展满足气动、隐身与控制等多学科影响的翼型设计时，需要在更宽范围内对设计变量的影响进行研究。本文以基于代理模型的全局优化设计(Surrogate Base Optimization，SBO)[5-12]方法为基础，开展了设计空间对翼型气动与隐身设计结果的影响分析，并提出了基于自适应参数化翼型气动隐身一体化优化设计方法。

影响气动外形优化结果的设计空间包括尺度和维度两个方面。为了解决空间尺度不足导致无法找到理想设计结果的问题，传统大多采用多轮优化方法[13]，即以上一轮优化的结果作为新的初始外形，扩大设计空间尺度再次进行优化，但这种“人在回路”的方法无疑会造成大量的计算资源消耗以及设计者工作量的增加。为此，王超等[6]在研究了设计变量每一维的尺度对设计空间及代理模型影响的基础上，提出了自适应设计空间尺度扩展方法，即在设计过程中选择性扩张设计空间尺度，有效提升了优化设计质量和效率。

然而在气动设计中，设计变量及其数量的选择对设计结果影响很大[14-15]。以往设计者多是根据经验或直接布置大量设计变量，但对于新型翼型而言，设计者一般缺乏必要的经验，而大幅增加设计变量，将对基于代理模型的全局优化造成“维度灾难”现象。

为了缓解设计空间维度扩张对优化算法的压力，设计者们开展了自适应参数化方法研究。Desideri[16]、Masters[17]等采用逐步加密设计变量的方法增加设计空间维度，但这样增加设计变量缺乏针对性，造成计算资源的浪费；Anderson[18]、Zingg等[19]以设计变量梯度信息作为敏感性信息，选择增加梯度较大的设计变量，但对于进化类算法而言，设计变量梯度难以直接获取，并且梯度信息只能反映设计变量在局部设计空间的敏感性。此外，设计空间维度扩张后，已有的低维空间产生的样本难以在高维空间内直接利用，重新在高维空间取样则会带来额外的计算负担，这更加限制了自适应参数化方法在全局优化算法中的应用。

针对现有的自适应参数化方法难以在代理模型全局优化中应用的问题，本文首先研究不同维度设计空间对翼型设计结果影响，对比分析了不同设计变量对翼型气动与隐身特性影响的敏感性，在此基础上提出了适用于代理模型全局优化方法的自适应参数加点方法：为了避免设计变量分布过度集中，引入密度阈值控制设计变量分布；针对维度扩展后低维样本无法利用的问题，应用节点插入算法将已有的低维样本转化为高维样本，避免了重新取样的工作量，实现了样本的高效配置。最后通过飞翼布局翼型气动隐身设计，并与不同设计空间维度设计结果进行对比，验证了该方法的可靠性。

1 基本方法

1.1 基本效应法

基本效应法能够凭借较少的计算量识别计算模型Y=f(X)中对输出量Y影响较大的输入变量X，该方法尤其适用于计算量较大的物理模型，是应用最为广泛的全局敏感性分析方法之一[20]。基本效应法通过OAT(One at A Time)[21]的方式计算任一维输入量Xi对输出量Y的影响，即任一维输入量在取值范围内每次变化Δi，其对输出量的敏感性信息计算方式为

EEi=(f(X1,X2,…，Xi+Δi，…，Xn)-

f(X1,X2,…，Xi，…，Xn))/Δi

(1)

文献[20]提出径向取样(Radial Sample)法,通过在整体设计空间内进行r次随机取样，可以得到每一维设计变量的全局敏感度信息SAi：

(2)

由式(2)可知，SAi的计算结果存在正负问题，可能导致计算结果不能充分反映设计变量的敏感性信息。因此，Campolongo等将式(2)转化为式(3),以消除正负号对结果的影响。

(3)

对于气动优化问题而言，通过优化过程获得样本中的设计变量对目标函数的敏感度分析，可以获得不同设计变量对目标特性影响的敏感度，并根据敏感度选择扩展设计空间维度。

1.2 Bspline 曲线参数化方法

Bspline[22]方法凭借灵敏的局部扰动特性，广泛应用于工业设计中。通过一组基函数及其系数实现对外形的拟合和扰动，基函数的系数称为控制点，其位置由横、纵坐标共同决定。在翼型参数化中通常采用三次样条曲线。

Bspline通常由一组u∈[0,1]的标量来组成。

(4)

(5)

(6)

式中：递增节点序列r={r0,r2，…,rl},ri≤ri+1，u={r0,…,r0,rk+1,…,rn-1,rl,…rl}，重节点r0、rl的数量均为k+1，其余节点一般均匀分布。

相较于其他经典参数化方法，如CST[23]、Hicks-Henne[24]等,Bspline曲线具备良好的局部扰动特性[25]。

由式(6)可以看出，Bspline曲线每段曲线的特性由相应节点矢量决定，因此Bspline曲线控制点的增加是从节点开始的。在节点区间[ri,ri+1]插入新的节点r，从而得到新的节点标量u1={r0,…,ri,r,ri+1,…，rl}，并由新的节点矢量确定新的基函数和控制点矢量[26]：

Pnew,j=

(7)

式中：αj=(t-tj)/(tj+k-1-tj)，这样一来，原来的样条曲线就可以用新的基函数与控制点矢量表示，并与初始曲线完全一致。控制点的增加仅影响其左右相邻控制点的位置。因此，Bspline曲线节点插入算法的特性能够实现低维样本在高维设计空间的重构，避免了重新取样，提高了样本利用率[26]。

1.3 二维矩量法

本文采用基于电场积分方程的二维矩量法[27]进行翼型RCS特性计算。

(8)

矩量法计入了各部分感应电流相互之间的影响，是麦克斯韦方程的精确求解方法。由于矩量法离散的面元数与照射波长成正比，对于飞行器这种大尺寸物体的高频计算，三维矩量法需要的计算量和存储量往往难以承受。而二维情况下计算量较小，虽然电磁散射的二维情况在实际中很少遇到，但是其对分析相应的三维问题具有很好的参考意义。

二维情况下，散射体为无限长柱体，雷达散射截面σ变为雷达散射宽度σ′：

(9)

式中：r为柱体到散射观察点的距离；ρ为表面电荷密度；Es和Hs均为目标散射电磁波在雷达处的电磁场强度。

三维情况下长度为有限值l的柱形结构的雷达散射截面σ与无限长相同结构的二维雷达散射宽度σ′之间有如下简单关系[27]：

(10)

式中：λ是照射波长。

在隐身数值计算中，EMSS公司发布的三维全波电磁仿真商用软件FEKO得到了较为广泛的认可。为了验证本文应用的二维矩量法的可靠性，下面以NACA65，3-018翼型为例，分别采用商用计算软件FEKO(三维矩量法)和二维矩量法程序，在入射频率分别为3 GHz、5 GHz、7 GHz，入射角度为±90°的条件下进行计算，其中，二维矩量法的计算结果经式(10)转化为三维结果，并与FEKO软件计算结果进行对比。FEKO计算模型为三维模型，如图1所示，模型展长为1 m。计算环境为Inter(R)Core(TM) i7-8700 CPU @3.2 GHz RAM 32 G。

图1 FEKO三维计算模型Fig.1 3D computation model of FEKO

图2为计算结果对比，其中横坐标为入射方位角θ。可以看出，二维矩量法程序与FEKO软件的计算结果吻合良好，计算精度满足要求。因此，在针对翼型的气动隐身协同设计中，本文应用二维矩量法程序进行翼型RCS特性的数值模拟。

图2 FEKO与二维矩量法计算结果对比Fig.2 Results comparison of FEKO and 2D moment method

2 设计空间维度对翼型气动隐身特性影响

2.1 不同设计空间维度设计结果对比

本文以相对厚度为18%的对称翼型NACA65，3-018为例[28]，选取典型跨声速设计状态马赫数Ma=0.70,升力系数CL=0.25。翼型气动特性要求在保持量化的抬头力矩时提高升阻特性；隐身特性要求在入射频率9 GHz下，前向±30°均值降低，入射角θ如图3所示。

图3 入射角示意图Fig.3 Sketch map of incident angle

设计目标为减小阻力系数CD和提高前向RCS特性，气动约束为力矩系数Cm不小于0.03，几何约束为翼型最大厚度thickmax不减小，设计模型为

objCD,RCS

(11)

在取样过程中，气动特性计算采用基于RANS方程的Cfl3d程序，隐身特性计算采用二维矩量法[27]。针对减阻和降低前向RCS设计目标，本文将两者进行正则化，然后采用线性加权的方式转换为一个目标进行优化搜索，权重比例为1∶1。

为了对比分析设计变量对目标特性的影响，本文首先采用设计变量均匀分布的方式，设计变量数(DV)即空间维度分别为12、18、24、30。利用拉丁超立方(LHS)方法在各设计空间内分别取样100个，加点500次。表1为不同设计空间维度下的设计结果，可以看出在满足力矩约束下，各组设计外形阻力特性和隐身特性都有很大的提升，图4为前向RCS对比。从图5对比中可以看出，低维设计空间(DV_12)收敛效率高，但低维空间内难以搜索到理想解，因此目标值很快陷入停滞；随着设计空间维度增加(DV_12→DV_24)，较之低维空间此时得到了更好的设计结果；当设计空间维度均匀增加30维(DV_30)时，由于设计变量分布缺乏针对性，且空间维度扩张导致代理模型精度降低，因此在相同的计算量下，得到的结果反而不如低维设计空间的结果。

图4 设计翼型前向RCS特性对比Fig.4 Forward RCS comparison of optimized airfoils

图5 目标值收敛历程Fig.5 Convergence history of fitness

表1 翼型计算结果Table 1 Computational results of airfoils

设计翼型外形和压力分布对比分别如图6、图7所示。各组设计结果外形接近：前缘半径都有所变小，具备明显的前缘正加载特征和中后部反加载，同时最大厚度位置后移。从结果对比中可知各组设计外形基本相似，说明优化搜索的方向是正确的，但设计变量的选择直接影响设计结果，简单地增加设计变量并不能确保改善设计结果。

图6 设计翼型外形对比Fig.6 Shape comparison of optimized airfoils

图7 设计翼型压力分布对比Fig.7 Pressure distribution comparison of optimized airfoils

2.2 设计变量关于气动隐身特性敏感性分析

设计空间由设计变量及其数量决定，为了分析设计空间内关于气动、隐身特性敏感的区域，本节以DV_12算例为例，对设计变量进行敏感性分析。

图8、图9分别为设计变量关于翼型阻力系数和前向RCS的敏感性评估结果，图中前6个为上表面设计变量，后6个为下表面变量。对于阻力特性而言，翼型上表面外形主要影响着激波的强弱和发展，因此这部分设计变量敏感性更强；对于隐身特性而言，上、下表面第1个设计变量的影响远大于其他变量，而这两个设计变量主要影响翼型前缘20%左右的区域，因此这部分区域对前向RCS特性影响更大。

图8 设计变量关于阻力系数敏感性Fig.8 Sensitivity to CD of design variables

图9 设计变量关于前向RCS敏感性Fig.9 Sensitivity to forward RCS of design variables

结合减阻和降低前向RCS的设计目标，为了得到更好的设计结果，需要在设计空间敏感区域扩展维度。而均匀分布的设计空间内，设计变量的分布缺乏针对性，为了得到理想的结果往往需要更多的设计变量，即更高的维度，而维度的扩展意味着设计空间呈指数增加，导致代理模型效率急剧降低，需要更多的样本才可能得到理想的设计结果。因此，在设计空间敏感的区域进行维度扩展，使得每个设计变量充分发挥作用，避免作用不大的设计变量的增加导致空间的过度增长，提高代理模型的效率的同时，得到更加理想的设计结果。

3 自适应参数化方法

对基于代理模型的全局气动优化而言，设计空间的构造对优化设计结果和效率有很大的影响，而设计者一般难以直接给出合适的设计空间。本文提出的自适应参数化方法，设计初始阶段只需凭借较少的设计变量定义低维设计空间，在优化过程中逐步在设计空间敏感区域增加维度，自适应构造设计空间，提升优化设计质量。

本文构建的自适应参数化方法基本框架如图10所示，具体步骤为

图10 自适应参数化方法流程Fig.10 Flowchart of adaptive parameterization method

步骤1初始设计空间确定，利用较少的设计变量确定低维设计空间。

步骤2初始样本库构建，利用LHS方法在初始设计空间内生成样本，并加入样本库中。

步骤3训练代理模型并优化，利用样本库内的样本训练代理模型，在当前设计空间内进行优化搜索，将搜索到的最优值进行评估，更新样本库和代理模型。

步骤4设计空间维度扩展判断，维度扩展需满足如下条件：① 目标值长期陷入停滞；② 代理模型精度满足要求。

本文定义目标函数值的更新效率Δt和代理模型精度f*，作为维度扩展的判断：

Δt=fgi-t+1-fgi

(12)

(13)

Tr=f*/t

(14)

本文设置目标函数值更新的阈值objtr，代理模型精度阈值c。

当Δt≤objtr，即目标函数优化效率不满足要求，以及Tr≤c，即代理模型精度满足要求的情况下，进行设计空间维度的扩展。

步骤5生成候补设计变量，进行维度扩展时，首先根据已有节点[ri,ri+1]生成新的节点r：

(15)

然后利用节点插入算法，在设计变量[Pi,Pi+1]间生成候补设计变量Pnew，在此过程中根据设计变量密度分布阈值dρ对候补设计变量进行第一次遴选：

Pnew=

(16)

式中：|Pnew-Pi|为设计变量横坐标距离。

将不满足密度分布阈值要求的候补变量去掉(即为空集[])，留下满足要求的变量进行下一步。

步骤6设计空间维度扩展，通过对候补设计变量敏感性分析进行第2次遴选，敏感性分析过程如第2节所示，将满足敏感度阈值的候补变量加入已有设计变量中，实现设计空间维度扩展。

步骤7确定新增设计变量取值范围，新增候补变量的取值范围VRnew由其相邻变量的取值范围VRi，以及设计变量横坐标距离|Pnew-Pi|共同决定，计算公式为

(17)

步骤8样本维度变换，结合新增候补设计变量的位置，同样利用节点插入算法将样本库内的低维样本转换为相应的高维样本，避免了在新的高维设计空间内重新取样的工作量。

步骤9继续进行优化，利用变换后的高维样本重新训练代理模型，在高维设计空间继续进行优化。

步骤10直至收敛，程序收敛条件为预设迭代数，满足上述条件即终止程序，否则转到步骤3直至满足停止条件。

4 NACA 65,3-018翼型气动隐身设计

采用本文的自适应参数化方法对NACA 65,3-018翼型重新进行优化设计，算例初始设置与DV_12保持一致，设计目标依旧为减阻和提高前向RCS特性，气动约束为力矩系数不小于0.03，几何约束为翼型最大厚度不减小，设计模型为

objCD，RCS

(18)

将设计结果与固定设计空间维度方法进行对比。

4.1 自适应参数化方法设计结果

表2为设计翼型与初始翼型的特性评估对比，相较于初始翼型，设计翼型阻力系数降低26 counts，力矩特性满足约束，隐身特性提升92%，气动与隐身特性均有明显提升。应用自适应参数化方法的优化收敛历程如图11所示，设计空间维度经历12-18-24-30的变化过程，在此过程中，目标函数最优值保持了较高的收敛效率，在120代左右达到收敛。

图11 自适应参数化方法目标值收敛历程Fig.11 Fitness convergence history of adaptive parameterization method

表2 自适应参数化方法设计翼型计算结果Table 2 Computational results of optimized airfoil of adaptive parameterization method

设计翼型与初始翼型的外形和压力分布对比分别如图12、图13所示。应用自适应参数化方法设计的翼型外形与第2节设计结果相似，前缘半径呈现了较为明显的“鹰嘴”特征，上表面趋于平坦，最大厚度位置后移，具备明显的前加载和中后部反加载；从压力分布对比可知，由于设计翼型的前缘外形呈现从钝头体到“鹰嘴”的变化，其头部尤其是最前缘的外形曲率变化较为剧烈，使设计翼型压力分布出现了较小的吸力峰，为了满足抬头力矩的约束，设计翼型压力分布前缘正升力区域扩大，中后部出现负升力区域，初始翼型的强激波削弱为两道弱激波。设计翼型的RCS特性如图14所示。

图12 设计翼型外形Fig.12 Shape of optimized airfoil

图13 设计翼型压力分布Fig.13 Pressure distribution of optimized airfoil

图14 设计翼型前向RCS特性Fig.14 Forward RCS of optimized airfoil

图15和图16展示了优化过程中节点和设计变量位置的变化。相较于初始，设计变量数量明显增加，并且设计变量的分布更加具有针对性。翼型前缘设计变量明显增加，有利于提升设计翼型的前向RCS特性；为了减弱激波，翼型上表面中部区域的设计变量分布更加密集；此外，翼型后缘下表面的变量也有增加，从最终外形可知，翼型后缘的设计变量使得其具有明显的后缘反加载特征，配合前缘的正加载，使设计翼型在满足力矩特性的约束下，获得更好的阻力和隐身特性。

图15 自适应参数化过程增加的节点Fig.15 Added knots in adaptive parameterization method

图16 初始与最终阶段设计变量横坐标位置对比Fig.16 Design variable position comparison between initial and final stages

4.2 设计结果对比

为了验证自适应参数化方法的可靠性，本文将其与固定设计空间维度的设计结果进行比较。图17和图18分别为各组设计翼型的外形和收敛历程对比，可以看出，各组设计的翼型外形特征相似。各组设计翼型的特性评估对比如表3所示，在满足力矩系数的约束下，应用自适应参数化方法调用的CFD计算次数更少，设计的翼型其阻力特性和隐身特性更优于固定设计空间维度的设计结果。