随着计算机性能的不断提高，高精度仿真分析被广泛应用于工程实践，能够提高生产效率，降低生产成本[1-2]；同时，高精度仿真模型也变得更加复杂，模型离散程度的增加，求解器保真度的增加，求解精度的提高等都使得高精度仿真仍然非常耗时．当优化问题涉及高精度仿真时，优化算法调用大量仿真模型的代价仍然很高．代理模型(surrogate model)是一种近似模型，能够建立设计变量和目标函数间的关系，以代替真实模型．代理模型的精度对研究的结果有重要影响，当对某一特定问题建立代理模型时，模型精度主要与代理模型技术和样本分布情况有关．克里金 (Kriging) 模型能够预测未知样本的不确定性，是应用最广泛的代理模型．样本在设计空间的分布由试验设计(design of experiments，DoE)确定，合理的样本分布有助于获得更多的信息，为建模提供基础．根据样本生成方式的不同，DoE可以分为一次性采样(one-shot sampling)和序列采样(sequential sampling)[3]．一次性采样根据采样准则确定样本在空间的分布，一般须要同时满足空间填充和投影特性[4]，代表性的方法有拉丁超立方设计(Latin hypercube design，LHD)、优化拉丁超立方设计(optimized LHD，OLHD)和平移传播拉丁超立方设计(translational propagation LHD，TPLHD)[5]等．一次性采样存在两个问题，即样本数量须要预先设置和没有充分利用已有样本信息．a. 在实际工程问题中，设计人员无法预知合适的样本数量．样本数量过多，须要进行更多的仿真计算，造成资源浪费；样本数量过少，建立的代理模型精度不足，无法对研究的问题提供准确预测．b. 已有样本及其响应值建立的代理模型能够提供预测信息与不确定性信息，在不确定性大的区域增加样本，能够提高模型精度．序列采样方法利用已有样本的空间分布信息和代理模型提供的信息迭代的增加样本，直至满足精度或者达到计算资源上限．已有研究表明，序列采样方法建立的代理模型精度高于相同样本数量的一次性采样[6-8]．文献[9]提出一种累积误差(accumulative error，ACE)序列采样方法，在非线性区域选择更多样本．文献[10]提出一种空间填充与交叉验证权衡的序列采样方法(space-filling cross-validation tradeoff，SFCVT)，通过最大化约束的交叉验证误差预测值确定新样本．文献[11]提出一种序列采样方法CV-Voronoi，利用Voronoi多边形分割设计空间，利用交叉验证(cross validation，CV)确定敏感多边形进行采样，测试结果表明该方法的性能优于ACE和SFCVT．文献[12]提出一种K折交叉验证误差的CV-Voronoi方法，提高了采样效率．文献[13]利用序列采样方法建立稳态视觉诱发电位的代理模型并进行优化，提出的方法比拉丁超立方设计和正交设计方法更高效．文献[14]回顾了建立全局代理模型的序列采样方法．现有大部分序列采样方法在一次迭代过程中仅增加一个样本，当存在并行计算资源时，每次迭代过程选择多个样本进行并行计算，可以极大缩短建模时间[15]．本研究提出一种基于空间分割的并行序列采样方法(parallel sequential sampling based on design space partition，PSS-DSP)，利用Voronoi图分割设计空间，得到一系列子空间；利用权重的误差指标和稀疏性指标判断子空间的样本需求度；根据采样顺序同时考虑与该子空间中已有样本以及与本次迭代中已加入样本的空间距离，按最大最小准则选择新增样本．基于空间分割的并行序列采样方法能够根据计算资源在每次迭代过程中选择任意数量的样本，考虑新增样本间的空间距离，避免样本聚集，提高序列采样效率．1 并行序列采样方法一种好的序列采样方法应该协同考虑全局探索与局部搜索特性，全局探索指在样本分布稀疏区域增加样本，局部搜索指在非线性程度较大区域增加样本，两种策略都有助于提高代理模型精度．基于空间分割的并行序列采样方法采用Voronoi图对设计空间进行分割，每个子空间含有一个样本，尺寸较大的子区域表示样本分布较稀疏．每个样本的交叉验证误差(CV error，CVE)作为子区域的非线性程度，较大的交叉验证误差意味着该样本存在与否对模型精度有很大影响，进而能够推测该区域的非线性程度较高．稀疏性指标和误差指标通过权重和的形式确定目标子空间，在每个目标子空间中选择新样本．1.1　基于Voronoi图的空间分割假设一个n维设计空间含有m个样本，X=[x1,x2,…,xm]⊆Rn，Voronoi图根据m个样本将整个空间分割为m个子空间，C=(C1,C2,…,Cm)．每个子空间Ci含有一个样本xi，该子空间中任意点与xi的距离小于其与其他m-1个样本的距离．子区域Ci的定义方式为dom(xi, xj)={x∈Rn|‖x-xi‖≤‖x-xj‖};                                            Ci=∩xj∈X\xidom(xi, xj), (1)式中：•表示两点间的欧氏距离；dom(xi,xj)表示以样本xi和xj连线的垂直平分线为边界的封闭超平面，超平面内部的点与xi的距离更近．分割后的子空间边界是不规则的，用理论方法不容易描述．文献[11]采用蒙特卡罗法描述多边形边界，具体步骤为：a. 在设计空间生成大量随机样本；b. 计算每个样本与已有样本的欧氏距离；c. 与样本xi距离最近的点则位于子空间Ci中．图1给出使用该策略得到的Voronoi图．为方便对比理论方法与统计方法得到的子空间大小，已有样本是规则分布的，由图中四个红色方块表示．不同颜色圆点代表生成的大量随机样本，红色实线代表不同子空间的边界线，能够看出不同颜色区域的差别与理论边界线基本重合．10.13245/j.hust.220607.F001图1利用统计方法得到的子空间利用子空间中随机样本的数量表示子空间尺寸，样本数量与其大小成正比例．表1给出了图1四个子空间中随机样本数量与其面积的关系．二维空间的面积为1，随机样本数量为2 000．表中给出每个子空间的理论面积(S0)、实际样本数量(m)和实际面积(S)．四个子空间的理论面积都是0.25，能够看出统计方法得到的实际面积与理论面积基本相同．10.13245/j.hust.220607.T001表1基于蒙特卡罗法得到的子空间尺寸i样本坐标mS0S1(0.25，0.75)5210.250.260 52(0.75，0.75)5000.250.250 03(0.25，0.25)4800.250.240 04(0.75，0.25)4990.250.249 51.2　目标子空间确定通过协同考虑每个子空间的稀疏性指标和样本误差指标判断目标子空间的样本需求度，确定须要增加样本的目标子空间．稀疏性通过子空间尺寸体现，子空间尺寸越大表明其稀疏程度越高，在该空间增加样本能够提高样本密度，减小不确定性．样本误差通过交叉验证误差体现，对于典型的留一交叉验证(leave-one-out CV，LOOCV)，一个样本的交叉验证误差可以表示为eCV,i=y(xi)-ŷ-i(xi)     (   i=1, 2,…,m), (2)式中：y(xi)为样本xi的真实响应；ŷ-i(xi)为xi不作为训练样本时建立的代理模型对xi的预测值；m为当前样本数量．一个样本的交叉验证误差越大，表明该样本的缺失会显著影响代理模型在该区域的精度．可以进一步推测该区域的非线性程度较高，或者样本分布较稀疏，须要在该样本附近区域增加样本．考虑稀疏性指标和样本误差指标的样本需求度函数可表示为Di=Ae˜CV,i+(1-A)V˜i       (  A∈[0,    1]), (3)式中：Di，e˜CV,i，V˜i分别为第i个子空间的样本需求度、归一化交叉验证误差、归一化稀疏度；A为两种因素间的权重，当A 0.5时，样本误差指标占有更大权重，反之，稀疏性指标占有更大权重．通常来说，减小交叉验证误差能够直接提高代理模型的精度，因此占有更高的权重．一个子空间的样本需求度越大，说明其内部已有样本的交叉验证误差较大，或者其尺寸较大，或者二者的权重和较大，须要优先在该空间中增加样本．e˜CV,i和V˜i的表达式为e˜CV,i= (ei-emin)/(emax-emin)    ; (4)V˜i=Vi-VminVmax-Vmin=Ni-NminNmax-Nmin      , (5)式中：emax和emin为m个交叉验证误差的最大值和最小值；Vmax和Vmin为m个子空间尺寸的最大值和最小值，由每个子空间中的随机样本数量表示；Ni为第i个子空间中的随机样本数量．1.3　新增样本确定在目标子空间确定后，须要在对应子空间中增加样本．为了获得更多信息，一个子空间仅增加一个样本[11]．当前大部分序列采样方法在一次迭代中仅增加一个样本；当存在并行计算资源时，可以同时计算多个样本的响应值，提高采样效率．本研究提出一种并行样本选择策略，根据计算资源选择新增样本数量，同时考虑新增样本间的空间关系，避免样本聚集．假设现有计算资源允许同时计算d个新样本的响应值，即一次迭代可以选择d个新样本．根据式(3)得到样本需求度最大的d个子空间，每个子空间选择一个新样本，新样本从该空间的随机样本中选择．新样本的选择方式为xinew=max(‖xj-xi‖)          (i=1,   j=1, 2,…,Ni);    max(min(‖xj-(xi,x1new,x2new,…,xi-1new)‖))(  i1), (6)式中：xinew为第i个子空间的新增样本；xi为该子空间内已有样本；xj为第j个随机样本．对于第一个新增样本(i=1)，仅考虑该空间中与已有样本(xi)距离最大的随机样本．对于第i个(i1)新增样本，要计算该子空间中每个随机样本与该空间已有样本(xi)和前i-1个新增样本之间的最小距离，选择最小距离最大的随机样本作为该子空间的新增样本．通过最大化后加入样本与本次迭代中先增加样本间的最小距离，能够避免各子空间独立增加新样本时的样本聚集情况．1.4　算法主要流程图2给出采样方法的主要流程，具体步骤如下．a. 利用优化拉丁超立方设计产生少量初始样本，计算样本的真实响应值，将样本及其响应值放入数据库；b. 由式(2)计算数据库中每个样本的留一交叉验证误差并根据式(4)归一化；c. 利用式(1)分割设计空间，利用蒙特卡罗法计算每个子空间尺寸并由式(5)归一化；d. 根据式(3)计算子空间的样本需求度；e. 根据计算资源确定d个目标子空间；选择新样本，并行计算样本响应值，加入数据库；f. 当代理模型满足精度或达到计算代价时，采样过程结束，建立全局代理模型，否则转到步骤b，继续进行采样．10.13245/j.hust.220607.F002图2基于空间分割的并行序列采样流程2 结果和讨论在此对提出的并行序列采样方法进行测试，验证其有效性．首先介绍参数设置，然后对权重系数A进行分析，最后与其他采样方法在标准测试函数中进行对比，验证基于空间分割的并行序列采样方法的有效性．2.1　参数设置使用Kriging模型建立代理模型，利用SURROGATES Toolbox工具箱[16]实现建模．Kriging模型采用0阶多项式、高斯核函数；超参数θ的初值为[1，1，…，1]n，取值范围为[1×10-3，20]n．当n3时，A取值为1；当n≥3时，A取值为0.5．测试集由7个函数组成，变量数量从2到6，初始样本数量为5n．对于并行采样，主要探讨每次迭代选择3个和5个样本的情况，因此增加的样本数量为15的倍数．采样方法在每个函数独立执行10次，得到统计结果．测试函数的相关信息见表2，函数的具体表达式参见文献[17]．10.13245/j.hust.220607.T002表2测试函数相关信息序号测试函数n初始样本数量(5n)增加样本数量设计空间F1Easom21060[0，7]nF2Peaks21060[-5，5]nF3Shubert21060[1,3]nF4Hartman3315120[0，1]nF5Shekel420180[3,5]nF6Ackley525240[-2，2]nF7Hartman6630300[0，1]n在相同函数评估次数下，对比建立的代理模型精度．一种序列采样方法的效果越好，建立的代理模型精度越高．测试标准为均方根误差(root mean square error，ERMS)，ERMS是一种全局精度标准，ERMS=1nt∑k=1nt(y(xk)-ŷ(xk))2,式中：y(xk)和ŷ(xk)分别为测试样本xk的真实响应值与预测响应值；nt为测试样本的数量，本研究所有函数的测试样本数量为5 000．2.2　参数分析基于空间分割的并行序列采样方法通过权衡考虑样本稀疏性指标和样本误差指标判断子空间的样本需求度，二者通过系数A体现不同的权重．本节通过对比A的不同取值对基于空间分割的并行序列采样方法的影响，得到A的较优取值．设置A的取值分别为0.0，0.2，0.5，0.8和1.0．当A=1.0时，采样方法只考虑样本误差指标的影响；当A=0.0时，只考虑样本稀疏性指标的影响；当A=0.5时，两种因素起到的影响相同．图3给出权重系数A的不同取值在测试函数中得到的均方根误差，每种方法在每个测试函数独立10.13245/j.hust.220607.F003图3权重系数A的变化对测试问题均方根误差的影响测试10次得到平均值．基于空间分割的代理模型并行序列采样方法分别进行1点、3点和5点采样，表明在一次迭代过程中选择1个、3个和5个新增样本．对于函数F1，F2和F7，随着A的增加，均方根误差保持下降，当A=0.8时结果基本稳定，表明增加样本误差指标权重能够提高代理模型精度．对于函数F3，1点和3点采样受A值影响较大，5点采样的结果受A值影响不明显．对于函数F4~F6，均方根误差随着A的增加，有先减小后增大的趋势，且大部分在A=0.5时取得最好结果．这说明样本误差指标与样本稀疏性指标共同影响目标子空间的确定，进而影响新增样本的选择．对于前三个二维问题，A倾向于取较大值(A≥ 0.8)，能够迅速减小非线性区域的不确定性．对于函数F4～F6，给予稀疏性指标和不确定性指标相同的权重(A=0.5)，能够得到比较好的结果．对于函数F7，虽然A=1.0时结果最好；相比于A取值0.0和0.2，A取值0.5和0.8时的结果也比较接近，仍然可以将其归到取值大于0.5的范围内．因此，当问题变量数小于3(n3)时，A取值为1.0；当n≥3时，A取值为0.5．本研究提出一种并行采样策略，在目标子空间确定后，考虑新增样本加入顺序，避免样本聚集．通过对比不采用该策略的并行加点策略，验证其有效性．并行加点的数量为3个和5个，对比策略在每个目标子空间中选择与该空间已有样本距离最大的随机样本，各个目标子空间互不影响．图4给出两种并行加点策略(3点和5点)在测试函数随迭代次数(I)得到的均方根误差变化情况．整体上看，两种策略没有明显区别．在迭代前期有一定差别，随着样本数量的增加，得到的均方根误差差别不大．从3点采样的迭代曲线看，随机选择策略在前四个函数上表现较好，提出的策略在剩下问题上表现较好．对于5点采样的迭代过程，能够看出提出的策略在迭代初期的收敛速度更快．该策略在选择新增样本时，同时考虑与已有样本和本次迭代中已加入样本的空间距离，避免新增样本在各子空间边缘聚集，更有效地体现了空间填充特性．这对于昂贵问题的建模更有优势，当计算资源很有限时，该策略能够利用更少的样本得到精度更高的代理模型．10.13245/j.hust.220607.F004图4不同新增样本选择策略得到的均方根误差收敛曲线2.3　对比测试对比基于空间分割的代理模型并行序列采样方法与其他采样方法对建立的代理模型精度的影响见表3．PSS-DSP-1，PSS-DSP-3和PSS-DSP-5表示每次迭代过程中增加1个、3个和5个样本，后两种变体为并行采样．CV-Voronoi为一种基于空间分割的序列采样方法，当PSS-DSP-1的权重系数A=1.0时，基于空间分割的代理模型并行序列采样方法与该方法相同．均方误差(mean squared error，MSE)方法在Kriging模型均方根误差最大处增加样本．OLHD为一次性采样方法，按照一定的数量间隔生成样本，建立代理模型并计算均方根误差，每个样本数量进行10次独立测试．对于每个测试函数，序列采样方法的10组初始样本相同，测试样本相同并且只有一组．10.13245/j.hust.220607.T003表3不同采样方法在测试函数中得到的结果对比(x¯±s)测试函数mPSS-DSP-1PSS-DSP-3PSS-DSP-5CV-VoronoiMSEOLHDEasom (F1)704.587 7×10-3±7.496 4×10-3±5.130 4×10-3±5.018 2×10-3±3.281 5×10-2±1.503 8×10-2±1.522 2×10-33.439 4×10-31.267 4×10-31.203 3×10-39.019 8×10-36.044 8×10-3Peaks (F2)700.177 3±0.270 8±0.160 0±0.158 9±0.638 8±0.601 4±5.885 7×10-22.967 3×10-13.812 1×10-25.472 5×10-21.265 9×10-11.096 1×10-1Shubert (F3)702.052 7±1.893 9±1.896 2±1.901 8±1.806 7±0.868 4±4.503 1×10-14.617 7×10-13.129 6×10-15.369 0×10-11.648 3×10-19.456 6×10-2Hartman3 (F4)1352.089 4×10-2±1.903 5×10-2±1.888 3×10-2±2.220 0×10-2±2.573 1×10-2±2.666 3×10-2±4.162 1×10-35.959 2×10-35.589 3×10-32.480 0×10-34.509 8×10-35.966 9×10-3Shekel (F5)2000.109 3±0.111 0±0.110 8±0.208 1±0.235 4±0.155 8±9.799 4×10-31.129 7×10-21.078 0×10-23.001 0×10-22.865 7×10-22.333 5×10-2Ackley (F6)2650.434 9±0.434 7±0.433 9±0.466 8±0.453 7±0.452 4±1.532 2×10-25.215 6×10-37.731 8×10-31.844 8×10-21.305 1×10-29.785 9×10-3Hartman6 (F7)3300.111 6±0.114 1±0.125 7±0.104 6±0.229 3±0.138 6±9.206 5×10-31.046 2×10-21.432 0×10-25.381 6×10-31.813 5×10-21.117 2×10-2图5给出不同采样方法建立的Kriging模型在测试函数上得到的均方根误差变化情况，能够看到随着样本数量增加，均方根误差逐渐减小，表明增加样本能够提高代理模型精度．均方误差方法在所有测试函数中得到的模型精度都比较差．在二维问题中观察样本分布，发现很多样本都位于设计空间边缘，不利于提高代理模型精度．除了函数F3，大部分序列采样方法得到的结果好于一次性采样方法(OLHD)，表明序列采样方法能够有针对性地确定样本的分布．对比PSS-DSP-1和CV-Voronoi，两种方法在二维问题上表现相同，在函数F4上的结果也基本符合．PSS-DSP-1在F5和F6两个函数中表现更好，由于考虑了样本分布，空间填充特性更好，收敛速度和结果都比较好．CV-Voronoi在函数F7上表现更好．10.13245/j.hust.220607.F005图5不同采样方法在测试函数上均方根误差收敛情况对比相比于单点序列采样方法，基于空间分割的代理模型并行序列采样方法的并行采样策略能够大大提高收敛速度．由于一次迭代过程中分别增加多个样本，均方根误差能够快速减小．表3给出六种采样方法在相同样本数量下得到的均方根误差统计结果，可以看到CV-Voronoi和OLHD分别在两个和一个函数上表现最好，基于空间分割的代理模型并行序列采样方法三种变体在剩下四个函数中表现最好．PSS-DSP-5在两个函数中表现最好，在提升模型精度和建模时间两方面都体现出很大的优势．通过对比，能够看到基于空间分割的代理模型并行序列采样方法是一种很有竞争力的并行序列采样方法．3 结语本研究提出一种基于空间分割的代理模型并行序列采样方法，在分割设计空间的基础上，协同考虑样本误差指标与稀疏性指标，确定须要增加样本的目标子空间；在并行采样过程中，考虑新增样本间的空间距离，避免聚集情况的发生．通过在测试函数中对比不同权重系数的结果，给出了A的较优取值．对比提出的并行采样与随机采样策略结果，发现提出的策略能够在采样前期得到较好的结果．最后，通过与一次性采样和若干序列采样方法对比，验证了提出的并行采样方法在采样效率和建模结果方面都有一定的优势．