1 引言
面对大多数复杂的工程结构设计问题，为获取设计参数与结构响应之间的关系，工程技术人员往往采用代理模型。实践表明，代理模型具有精度高、计算量小，低误差等特点[1]，可广泛应用于航空航天等领域。现有较为常用代理模型主要有：响应面法[2]、神经网络法[3]、支持向量机[4]和Kriging模型[5-7]等。由于Kriging模型具有良好处理非线性预测能力，Kriging模型已成为最具代表的代理模型之一。
为完善Kriging模型，研究者已经进行了大量的研究工作。增强梯度Kriging模型[8]，优点在于通过利用梯度信息来提高Kriging模型的计算精度。Kennedy等人将协同Kriging模型[9]推广到工程领域，该模型通过用容易抽样的量替换较困难抽样的量来进行辅助预测。Joseph等人提出了盲Kriging模型[10-12]，通过贝叶斯法识别趋势函数。Kwon等人提出了趋势Kriging模型[13]，但趋势函数求解过程不利于工程实际问题的求解。
本文提出了一种改进的Kriging模型。与以往的Kriging模型不同的是：
1) 在生成初始样本后，本文通过响应面法确定趋势函数。考虑到实际工程问题，本文使用一至三阶不含交叉项的多项式函数作为备选趋势函数，通过回归分析计算出趋势函数的相关系数，选取趋势拟合程度较高的备选趋势函数的阶数作为Kriging模型整体趋势函数的阶数。
2) 在确定趋势项的阶数后，将整体趋势函数中所忽略的交叉项重新加入Kriging模型中，通过遗传算法重新筛选基函数，排除对模型计算精度无影响或者影响极小的基函数项，保留影响较大基函数项，降低计算量的同时提高模型的精度。
2 本文方法
经典Kriging模型包含整体趋势函数和随机函数2个部分。已有研究表明，在样本数量足够多的情况下，整体趋势函数对Kriging模型的预测精度影响较小[9-13]。然而足够多的样本必然带来计算成本的增加，不利于工程实际。因此，在样本数量有限的情况下，正确选取整体趋势函数有助于提高Kriging模型的预测精度。本文首先利用响应面方法确定Kriging模型整体趋势函数的阶次，然后通过遗传算法对整体趋势函数中的各项进行优化选取。
2.1 初始样本点
利用拉丁超立方抽样方式生成初始样本，并求得对应样本点的响应值。由于在最佳基函数选择中需要对全基函数进行计算对比，因此初始样本点的数量一定要大于等于全基函数的项数，初始样本点的数量为：
k1≥(n+p)!/n!×p!
(1)
式(1)中：k1为初始样本点数量；n是函数的维数；p是函数的阶数。
2.2 响应面法
响应面函数设计思路应当形式简单，并使待定系数尽量少，以减小结构分析的工作量，因此在建立整体趋势函数时忽略交叉项。
考虑到实际工程问题需求，本文建立一至三阶不含交叉项的多项式作为备选趋势函数。一阶、二阶、三阶的响应函数为：
(2)
式(2)中，a、bi、ci和di为待定系数。
2.3 趋势判别
通过回归分析，计算出所有备选趋势函数的相关系数，选取趋势拟合程度最高的趋势备选函数的阶数作为Kriging模型中趋势函数的阶数。相关系数为：
(3)
式(3)中： yi是样本点处的真实响应；
是真实响应值的期望；
是备选趋势函数。
选取相关系数最趋近于1的备选趋势函数的阶数作为Kriging模型中趋势函数的阶数。
2.4 最佳基函数选择
对于一个二维三阶函数，它的全基函数为，共计10个基函数项
假设真实函数为
则
均为不必要基函数项，而基函数中不必要的项不但增加计算成本，甚至可能降低模型的质量。因此本文采用了遗传算法，通过排除Kriging代理模型趋势函数中不必要的项，降低计算量的同时提高模型的精度。
Martin等[14]开发的二进制编码非支配排序遗传算法，可以搜索基函数的最优子集。二进制编码的生成一个(n+p)!/n!×p!位编码，其中n是函数的维数，p是函数的阶数。全基函数中的每个多项式项都被分配给不同的编码。是否使用指定的多项式项由二进制数决定，其中1表示选择多项式项作为候选项，0表示没有选择多项式项。可定义目标函数为均方根误差最小，来比较每个基函数的适应值，通过迭代选择，直到满足收敛条件，即可从全基函数项中找到基函数的最佳子集。
2.5 主动学习函数及收敛条件
当基函数选取后，利用主动学习函数逐步增加具有重要意义的样本，进而更新已构建的Kriging模型，直至满足一定收敛条件。常见的学习函数有EEF函数、U函数和H函数。实践表明[15]，U函数和其他学习函数相比，在保持较高精度的同时，所需的样本数量较少，因此本文采用U函数，U函数的表达式为：
(4)
式(4)中：
是Kriging代理模型的预测值；
是Kriging代理模型的预测值标准差。
学习函数的停止条件(即Kriging模型的收敛条件)为：
min.(U(X))>2
(5)
2.6 基本流程
本文方法的基本流程如下：
1) 使用拉丁超立方抽样方式生成初始样本；
2) 通过响应面法求得一至三阶响应面函数，作为备选趋势函数；
3) 对所有备选趋势函数进行回归分析，分别计算对应R2值并判别，选取趋势拟合程度较高的备选趋势函数的阶数作为Kriging模型整体趋势函数的阶数；
4) 使用遗传算法搜索趋势函数的基函数最优子集，目标函数为均方根误差最小；
5) 通过主动学习函数训练已构建的Kriging模型，直至满足收敛条件。
3 算例
本文给出3个数值算例，包括三角函数、二维函数和三维函数，以验证本文提出的改进 Kriging模型(improved kriging model，IKM)的计算效率与计算精度。采用误差均方根衡量计算精度，采用总样本数量衡量计算效率，将IKM与经典Kriging模型(classical kriging model，CKM)、泛Kriging模型(universal kriging model，UKM)和趋势Kriging模型(trended kriging model，TKM)进行对比。计算表格中的UKM1代表一阶UKM、UKM2代表二阶UKM、TKM1代表一阶TKM、TKM2代表2阶TKM，k是样本的总数量。
k=k1+k2
(6)
式(6)中：k1是初始样本点的数量;k2是由主动学习函数筛选来的后续样本点数量。
3.1 测试函数1
测试函数1来源于文献[13]，有：
(7)
图1为测试函数1的真实响应值的等高线图，图2为本文方法给出的响应值等高线图。对比图1、图2可知，IKM的响应值与真实函数的响应值总体趋势吻合较好。
图1 测试函数1真实响应值等高线图
Fig.1 Real contour of test function 1
图2 基于IKM的测试函数1响应值等高线图
Fig.2 Contour of test function 1 using IKM
趋势函数计算结果见表1，IKM与其他模型的总样本数及对比见表2。通过表1可知，当代理模型整体趋势函数阶数等于2时，趋势拟合程度最高。由表2可知，本文算法需要12个初始样本和3个后续样本点，本文方法的计算效率和其他算法的计算效率相同，计算精度有显著提升。
表1 趋势函数的计算结果(测试函数1)
Table 1 Results of trend function (test function 1)
参数数值取值域x1∈[1,2], x2∈[1.5, 3]一阶趋势函数k1=10, R2=0.998 38二阶趋势函数k1=12, R2=0.998 41三阶趋势函数k1=15, R2=0.997 61整体趋势函数的阶数2趋势函数最佳基函数项[x21,x22,x1,x2,1]
表2 各种计算结果(测试函数1)
Table 2 Comparison of results (test function 1)
模型样本kRMSECKM[13]150.869 846UKM1[13]150.857 204UKM2[13]150.774 842TKM1[13]150.857 206TKM2[13]150.774 840IKM(本文)12+30.751 709
3.2 测试函数2
测试函数2来源于文献[13]，有：
0.1(x1-3)2+0.1(x2-2)2
(8)
图3为测试函数2的真实响应值的等高线图，图4为本文方法给出的响应值等高线图。对比图3、图4可知，IKM的响应值与真实函数的响应值总体趋势吻合较好。
图3 测试函数2真实响应值等高线图
Fig.3 Real contour of test function 2
图4 基于IKM的测试函数2响应值等高线图
Fig.4 Contour of test function 2 using IKM
趋势函数计算结果见表3，IKM与其他模型的总样本数见表4。由表4可知，本文算法需要20个初始样本和5个后续样本点，计算效率比其他算法略低，但计算精度有显著提升。
表3 趋势函数的计算结果(测试函数2)
Table 3 Results of trend function (test function 2)
参数数值取值域x1∈[-5,10], x2∈[0,15]一阶趋势函数k1=20, R2=0.398 98二阶趋势函数k1=20, R2=0.890 51三阶趋势函数k1=20, R2=0.989 10确定整体趋势函数的阶数3趋势函数最佳基函数项[x31,x32,x21x2,x1x22,x21,x22,x1,x2,1]
表4 各种计算结果(测试函数2)
Table 4 Comparison of results (test function 2)
模型样本kRMSECKM[13]20108.415 40UKM1[13]20113.073 43UKM2[13]2093.582 74TKM1[13]20113.105 26TKM2[13]2093.577 40 IKM(本文)20+590.777 90
3.3 测试函数3
测试函数3来源于文献[16]，有：
(9)
趋势函数的计算结果见表5，IKM与其他模型的总样本数见表6。通过表5可知，当代理模型整体趋势函数阶数等于3时，趋势拟合程度最高。由表6可知，本文算法需要40个初始样本和5个后续样本点，虽然计算效率比其他算法略低，但是计算精度有显著提升。
表5 趋势函数的计算结果(测试函数3)
Table 5 Results of trend function (test function 3)
参数数值取值域x1∈[-5,10], x2∈[0,15], x3∈[-5,10]一阶趋势函数k1=35, R21=0.439 819二阶趋势函数k1=35, R22=0.869 01三阶趋势函数k1=40, R23=0.997 90确定整体趋势函数的阶数3趋势函数最佳基函数项[x31,x32,x21x2,x22x3,x21,x22,x23,x1,x2,x3,1]
表6 各种计算结果(测试函数3)
Table 6 Comparison of results (test function 3)
模型样本kRMSECKM[16]3511.269 9UKM1[16]3511.592 8UKM2[16]359.517 5 IKM(本文)40+58.734 4
4 结论
在Kriging模型构建过程中，首先利用不含交叉项的多项式响应面法确定趋势函数的阶次，然后将趋势函数中忽略的交叉项重新加入Kriging模型的整体趋势函数中，通过遗传算法筛选趋势函数的基函数项，排除整体趋势函数中不必要的项，降低Kriging模型的计算量的同时提高模型的计算精度。
算例结果表明：由于本文考虑了整体趋势函数，在计算量增加不多的情况下，计算精度明显高于其他方法。然而值得注意的是，为了便于最佳基函数的选取，本文选取的初始样本点数量等于全基函数的项数，这将导致样本点数量随着变量个数的增加出现急剧增加的现象。对于这一问题，需对基函数选取方法进行完善，将在后续工作中进行改进。
[1] 穆雪峰,姚卫星,余雄庆,等.多学科设计优化中常用代理模型的研究[J].计算力学学报,2005,22(05):608-612.
Mu X F，Yao W X，Yu X Q，etal.A survey of surrogate models used in MDO[J].Chinese Journal of Computational Mechanics，2005,22(05):608-612.
[2] 赵维涛，邱志平.基于合理子域的改进响应面方法[J].力学学报，2014，46(03):409-416.
Zhao W T，Qiu Z Z.An improved response surface method based on the reasonable subdomain[J].Chinese Journal of Theoretical and Applied Mechanics，2014，46(03):409-416.
[3] 孙大鑫,牛余宝,常浩.基于集成神经网络的航空维修差错风险评估模型[J].兵器装备工程学报，2009,30(08)：86-88.
Sun D X，Niu Y B，Chang J.Aviation maintenance errors risk assessment model based on integrated neural network[J].Journal of Ordnance Equipment Engineering，2009,30(08)：86- 88.
[4] 谭志杨,林友红,张小海,等.基于粗糙集的支持向量机回归建模及应用[J].兵器装备工程学报,2010,31(08)：53-55.
Tan Z Y，Lin Y H，Zhang X H.[J].Support vector machine regression modeling and application based on Rough Set[J].Journal of Ordnance Equipment Engineering，2010,31(08)：53-55.
[5] 黄晓旭,陈建桥.基于主动学习Kriging模型的可靠度分析[J].固体力学学报，2016,37(02):172-180.
Huang X H，Chen J Q.Reliability analysis based on active learning Kriging model[J].Chinese journal of solid mechanics，2016,37(02):172-180.
[6] 林智文,葛建立,张鸿浩,等.迫击炮炮箍位置优化设计[J].兵器装备工程学报,2017,38(01):79-82,87.
Lin Z W，Ge J L，Zhang H H，etal.Optimization design for gun hoop position of mortar[J].Journal of Ordnance Equipment Engineering，2017,38(01):79-82,87.
[7] 王浩，张海瑞,王尧,等.基于主动学习Kriging的飞行器射程评估[J].兵器装备工程学报,2019,40(05):56-60.
Wang H，Zhang H R，Wang Y，etal.Vehicle range assessment based on active learning Kriging[J].Journal of Ordnance Equipment Engineering，2019，40(05)：56-60.
[8] Laurenceau J，Sagaut P.Building efficient response surfaces of aerodynamic functions with Kriging and Cokriging[J].American Institute of Aeronautics and Astronautics Journal，2008,46(02):498-507.
[9] Kennedy M，O’nagan Y.Predicting the output from a complex computer code when fast approximations are available[J].Geoderma,2008,146(01):391-396.
[10] Joseph V R，Hung Y，Suduanto A.Blind Kriging：A new method for developing metamodels[J].Journal of Mechanical Design，2008,130(03)：350-353.
[11] 韩忠华.Kriging模型及代理优化算法研究进展[J].航空学报,2016,37(11):3197-3225.
Han Z H.Kriging surrogate model and its application to design optimization：A review of recent progress[J].Acta Aeronautica et Astronautica Sinica，2016,37(11):3197-3225.
[12] Zhao L,Choi K K，Lee I.Metamodeling method using dynamic Kriging for design optimization[J].American Institute of Aeronautics and Astronautics Journal，2014，52(10)：2313-2327.
[13] Kwon H I,Choi S I.A trended kriging model with R2 indicator and application to design optimization[J].Aerospace Science and Technology.2015,43:111-125.
[14] Martin J D，Simpson T W.Use of kriging models to approximate deterministic computer models[J].American Institute of Aeronautics and Astronautics Journal,2005，4,(43)：853-863.
[15] 赵维涛,陈欢,祁武超.基于主动学习Kriging模型的结构多失效模型可靠度计算[J].计算力学学报,2020，37(01)：8-13.
Zhao W T，Chen H，Qi W C.Structural reliability calculation for multiple failure modes based on an active learning Kriging model[J].Chinese Journal of Computational Mechanics，2020，37(01)：8-13.
[16] Lee H，Lee D J，Kwon H.Development of an optimized trend kriging model using regression analysis and selection process for optimal subset of basis functions[J].Aerospace Science and Technology,2018，77:273-285.