Bucher和Bourgund[7]使用不带交叉项的二次多项式拟合极限状态函数，并将响应面法应用于结构可靠度分析。Rajashekhar和Ellingwood[8]通过自适应迭代进一步对响应面法进行改进。Kaymaz和McMahon[9]赋予极限状态函数绝对值小的样本点更高的权重，采用加权回归的方法确定不含交叉项的二次多项式响应面函数中的待定系数。赵洁和吕震宙[10]提出一种新的权重矩阵构造方法，并采用线性多项式作为响应面函数。Nguyen等[11]提出双加权响应面法，该方法考虑了2个权重因子(样本点到真实响应面的距离和样本点到验算点的距离)，并且第1次迭代响应函数采用一次多项式，后续迭代中采用含交叉项的二次多项式。李贵杰等[12]在考虑联合概率密度函数值和极限状态函数值样本权数的加权线性响应面法基础上，利用鞍点逼近方法近似其累积分布函数和概率密度函数后再进行可靠性分析。钟宏林等[13]在加权响应面法的基础上，为了避免回归矩阵的病态，提出了新的权重表达式。李广博[14]采用Fourier正交基神经网络响应面代替传统二次多项式，并结合加权回归求解未知系数。张吉宁等[15]在加权非线性响应面法的基础上，通过向量梯度投影方法改进样本点的选择。张学刚[16]在迭代过程中采用混合响应面法，第1次迭代响应面函数采用一次多项式，后续迭代中采用不含交叉项的二次多项式，在迭代收敛后运用样本重用思想选取已有样本点并对其进行赋权来拟合含交叉项的二次多项式响应面函数。贾长安[17]将连续插值取样和加权响应面法相结合提出了一种随机变量相关情形下的可靠度计算方法。侯振兴[18]构造了一种指数型线性权数形式，同时将均匀试验设计法应用于加权线性响应面法。Lu等[19]提出加权回归的极值响应面法进行结构动态模糊可靠性分析。洪林雄等[20]提出一种考虑样本点与设计点距离和样本点极限状态函数值大小的混合加权方法，并采用变向搜索算法寻找设计点。效率和精度的平衡兼顾是响应面法研究中的难点。目前, 应用最为广泛的响应面表达式是不含交叉项的二次多项式，相较于含有交叉项的二次多项式，需要的样本点少，但是计算精度较低。此外，已有的响应面法大多在迭代中都只采用新样本点的信息，而忽视已有的样本点信息，造成大量资源的浪费。
针对上述不足，本文提出一种改进的加权响应面结构可靠度计算方法，运用加权回归方法并结合文献[21]中提出的样本重复使用思想，在加速迭代收敛、提高计算效率的同时，可以较为精确地计算结构可靠度。首先，同时考虑样本点与验算点距离、极限状态函数值、联合概率密度函数值3个权重因子对各样本点赋权；其次，在迭代过程中重复利用之前产生的样本点以更新不含交叉项的二次多项式响应面函数；最后，在迭代收敛后，再次利用已有的样本点信息，采用加权回归方法拟合含交叉项的二次多项式响应面函数并计算结构可靠度。本文方法考虑了不同的权重因子，通过加权回归、样本多次重用和最后采用含交叉项的二次多项式拟合响应面函数，兼顾了响应面法的计算效率和精度，算例的计算结果验证了所提方法的可行性和优越性。
1 加权响应面法 在传统的响应面法中，一般采用最小二乘法求解响应面函数中的待定系数。而加权响应面法通过对各样本点赋权，基于加权回归的统计思想求解响应面函数中的待定系数，如式(1)所示。通过赋予更接近失效面的样本点更高的权重，可以达到使越接近失效面的样本点在拟合响应面函数时起更重要作用的目的，从而提高拟合精度[9]。
 | (1) |
式中: A 为回归矩阵; W 为权重矩阵; g 为样本点的极限状态函数值; b 为响应面函数中的待定系数。
响应面函数不含交叉项时，回归矩阵的表达式为
 | (2) |
响应面函数含交叉项时，回归矩阵的表达式为
 | (3) |
非线性加权响应面法的流程如图 1所示[22]。
 |
| 图 1 非线性加权响应面法流程[22] Fig. 1 Flowchart of weighted nonlinear response surface method[22] |
| 图选项 |
2 改进加权响应面法及其算法流程 2.1 响应面表达式的选择 不含交叉项的二次多项式是目前应用最为广泛的响应面函数表达式，在一定程度上可以反映极限状态函数的非线性程度，应用方便，只需较少的结构分析次数，但是对于非线性程度较大的极限状态函数，其计算误差较大，同时计算精度较低。虽然含有交叉项的二次多项式的计算精度较高，但因响应面函数中待定系数增多，使得求解需要的样本点增多，计算效率较低。本文结合这2种响应面函数表达式，在迭代过程中，采用不含交叉项的二次多项式拟合响应面函数，而在迭代收敛后采用含有交叉项的二次多项式拟合响应面函数，根据最终得到的含有交叉项的二次多项式计算可靠度。
不含交叉项的响应面函数表达式为
 | (4) |
式中：x =[x1, x2, …, xn]为组成极限状态函数的n个随机变量；b0、bi、bn+i为待定系数，共2n+1个。
含交叉项的响应面函数表达式为
 | (5) |
式中：b0、bi、bij为待定系数，共
2.2 权重矩阵的构造 样本点的权重决定了该样本点在加权回归分析中的作用，相对越重要的点，应该被赋予更大的权重。为了提高理想样本点并弱化不理想样本点在拟合响应面函数中所起的作用，对各样本点合理赋权，提出一种同时考虑样本点与验算点距离、极限状态函数值、联合概率密度函数值3个权重因子的权重矩阵构造方法。
在结构可靠性分析中，希望响应面函数能更好地逼近失效面，故第1个权重因子用来表征该样本点处极限状态函数绝对值|g(x)|的大小，|g(x)|越小则表明该样本点越靠近失效面。该权重表达式赋予接近失效面的样本点较大的权重，远离失效面的样本点较小的权重，第j个样本点的第1个权重因子表达式为
 | (6) |
式中：gbest为已有样本点中极限状态函数绝对值的“最佳”值，即除去0以外的最小值；g(xj)为各样本点处的极限状态函数值。
 | (7) |
式中：wgj为各样本点处第1个权重因子。
为了使响应面函数在概率上更好地替代失效面，第2个权重因子用来表征该样本点处联合概率密度函数值f(x)的大小，目的是赋予联合概率密度函数值较大的样本点更大的权重，反之赋予样本点的权重越小，第j个样本点的第2个权重因子表达式如下：
 | (8) |
式中：fbest为已有样本点中联合概率密度函数值的“最佳”值，即最大值；f(xj)为各样本点处联合概率密度函数值。
 | (9) |
式中：wfj为各样本点处第2个权重因子。
第3个权重因子用来表征该样本点与验算点之间的距离。该权重表达式赋予更接近验算点的样本点更大的权重，使其在回归分析中起到更大的作用。参考文献[20]中的方法，权函数以样条函数的形式给出：
 | (10) |
式中：dj为各样本点与验算点的距离；x*(k)为第k次迭代中的验算点；xj为样本点。
 | (11) |
式中：dmax为已有样本点中样本点与验算点的最大距离。
 | (12) |
 | (13) |
式中：wdj为各样本点处第3个权重因子。
取3个权重的平均加权，得到最终各样本点处的混合权重：
 | (14) |
以各样本点的混合权重为对角线元素构造权重矩阵：
 | (15) |
2.3 算法流程 步骤1????第1次迭代(k=1)采用经典响应面法，初始迭代一般取均值点μx=[μx1, μx2, …, μxn]为抽样中心，计算得到新抽样中心x(1)和可靠度指标β(1)。
步骤2????在第k次迭代中，采用Bucher方法选取初始样本点，即x(k-1)=(x1(k-1), x2(k-1), …, xn(k-1))和(x1(k-1), x2(k-1), …, xi(k-1)±f(k)σi, …, xn(k-1))(i=1, 2, …, n)，并计算极限状态函数的真实值g(x), σi为随机变量的标准差。x(k-1)为上一次迭代中计算得到的抽样中心，f(k)为第k次迭代的插值系数。
步骤3????用第k次迭代新增的2n+1个样本点及第k-1次迭代的(k－1)×(2n+1)个样本点，共计m=k×(2n+1)个样本点，由式(2)构造回归矩阵 A 。
步骤4????根据各个样本点的极限状态函数值、联合概率密度函数值以及与验算点的距离计算3个权重因子，由式(14)确定各个样本点的混合权重并由式(15)构造权重矩阵 W 。
步骤5????进行加权回归分析，由式(1)拟合不含交叉项的响应面函数，求解响应面函数中的待定系数。
步骤6????根据不含交叉项的响应面函数，采用一次二阶矩法计算验算点x*(k)和可靠度指标β(k)。
步骤7????用均值点和验算点进行线性插值，采用插值法计算下一次迭代的抽样中心x(k)：
 | (16) |
步骤8????判断
步骤9????迭代收敛后，在已有的样本点中选取权重最大的nw个样本点(nw≥含交叉项响应面函数中未知系数的个数)，由式(3)和式(15)根据已有样本点信息重新构造回归矩阵 A 和权重矩阵 W 后，再由式(1)加权回归拟合含交叉项的响应面函数并计算可靠度。
改进加权响应面法流程如图 2所示。
 |
| 图 2 改进加权响应面法流程 Fig. 2 Flowchart of improved weighted response surface method |
| 图选项 |
3 算例 3.1 算例1 指数型极限状态函数表达式为
 | (17) |
式中：x1、x2均服从标准正态分布。
图 3为改进加权响应面法的求解过程。经过3次迭代达到理想的收敛，迭代收敛后用含交叉项的响应面函数再次拟合，对其求解得到最终的可靠性指标。
 |
| 图 3 算例1求解过程 Fig. 3 Solving process of Example 1 |
| 图选项 |
用蒙特卡罗法、经典响应面法、线性加权响应面法、单加权响应面法、双加权响应面法和改进加权响应面法求解得到的计算结果如表 1所示。由表可以看出，经典响应面法所需结构分析次数最多，误差较大；由于该算例非线性程度较高，线性加权响应面法虽然快速收敛，但计算误差最大，为7.178 36%；单加权响应面法一定程度上减少了结构分析次数，但计算误差比经典响应面法更大；双加权响应面法相比经典响应面法在效率和精度上都有所提升；改进加权响应面法只需15次结构分析，仅为经典响应面法结构分析次数的一半，且计算误差最小，仅为1.849 81%。可以看出，由于改进加权响应面法同时考虑了样本点与验算点距离、极限状态函数值、联合概率密度函数值3个权重因子并对样本点多次重用，在第2次迭代后既快速收敛，同时精度达到了较理想的状态。
表 1 算例1计算结果比较 Table 1 Comparison of calculation results of Example 1
| 方法 | 失效概率 | 误差/% | 结构分析次数 | 计算时间/s |
| 蒙特卡罗法 | 0.003 622 | 0 | 10 000 000 | |
| 经典响应面法 | 0.003 369 | 6.985 09 | 30 | 0.061 |
| 线性加权响应面法 | 0.003 362 | 7.178 36 | 20 | 0.023 |
| 单加权响应面法 | 0.003 366 | 7.067 92 | 20 | 0.026 |
| 双加权响应面法 | 0.003 393 | 6.322 37 | 20 | 0.034 |
| 改进加权响应面法 | 0.003 689 | 1.849 81 | 15 | 0.048 |
表选项
3.2 算例2 带有集中力的悬臂梁如图 4所示。弹性模量E、截面惯性矩I以及施加的载荷力P均为服从正态分布的独立随机变量：E~N(2×107, 0.5×107)kN/m2，I~N(1×10－4, 0.1×10－4)m4，P~N(8, 2.5)kN。长度为常量L=5 m。考虑该悬臂梁的最大变形作为极限状态，建立极限状态函数如下：
 | (18) |
 |
| 图 4 带有集中力的悬臂梁 Fig. 4 A cantilever beam with concentrated force |
| 图选项 |
采用改进加权响应面法计算该悬臂梁的失效概率，迭代过程如图 5所示。其他对比方法同算例1，对比表 2中各方法的计算结果，发现改进加权响应面法计算得到的失效概率与蒙特卡罗模拟107次得到的理论解基本一致，说明改进加权响应面法在计算精度方面更具优势，失效概率的计算误差仅为0.134 25%。在效率方面，改进加权响应面法结构分析次数最少，只需21次。算例结果说明，在实际工程应用中改进加权响应面法依然具有较大的优势。
 |
| 图 5 算例2求解过程 Fig. 5 Solving process of Example 2 |
| 图选项 |
表 2 算例2计算结果比较 Table 2 Comparison of calculation results of Example 2
| 方法 | 失效概率 | 误差/% | 结构分析次数 | 计算时间/s |
| 蒙特卡罗法 | 0.005 959 | 0 | 10 000 000 | |
| 经典响应面法 | 0.005 691 | 4.497 40 | 56 | 0.026 |
| 线性加权响应面法 | 0.005 691 | 4.497 40 | 28 | 0.021 |
| 单加权响应面法 | 0.005 689 | 4.530 96 | 49 | 0.031 |
| 双加权响应面法 | 0.005 861 | 1.644 57 | 21 | 0.030 |
| 改进加权响应面法 | 0.005 967 | 0.134 25 | 21 | 0.052 |
表选项
3.3 算例3 十杆桁架结构如图 6所示。该桁架在节点5与节点6处简支，节点2与节点4处承受相同大小的竖直向下载荷P=105 lb(1 lb=0.453 592 37 kg)，水平杆和垂直杆的杆长均为L=360 in(1 in=0.025 4 m)。结构材料为铝合金，弹性模量E=107 psi(1 psi=6.895 kPa)。水平杆、垂直杆以及倾斜杆取自3种不同横截面积的铝合金杆件，基本随机变量为3种杆件结构的横截面积：S1~N(13, 1.3)in2，S2~N(2, 0.2)in2，S3~N(9, 0.9)in2(1 in2=0.092 903 04 m2)。节点2处的允许位移为dallow=3.5 in，则该十杆桁架结构的极限状态函数表达式为
 | (19) |
 |
| 图 6 十杆桁架结构 Fig. 6 Ten-bar truss structure |
| 图选项 |
式中：S =[S1, S2, S3]，S1、S2、S3分别为水平杆、垂直杆以及倾斜杆的横截面积；dreal为节点2处的实际位移值。
图 7为单次ANSYS仿真计算结果，节点2处的实际位移为0.096 495 m，即3.799 016 in。图 8为改进加权响应面法的求解过程。利用改进加权响应面法对十杆桁架求解可靠度，并和蒙特卡罗法、经典响应面法、线性加权响应面法、单加权响应面法、双加权响应面法的结果进行对比，如表 3所示。其中，蒙特卡罗法的计算结果由文献[23]提供，作为精确解。改进加权响应面法的结构分析次数较经典响应面法显著减少，并且计算结果的精度更高，可以看出对于非线性较强的结构其优势较显著。线性加权响应面法、单加权响应面法、双加权响应面法在一定程度上也降低了结构分析次数，但误差均大于改进加权响应面法。该工程算例的计算结果进一步验证了改进加权响应面法在计算效率和精度上的优越性。
 |
| 图 7 十杆桁架ANSYS仿真结果 Fig. 7 ANSYS simulation results of ten-bar truss |
| 图选项 |
 |
| 图 8 算例3求解过程 Fig. 8 Solving process of Example 3 |
| 图选项 |
表 3 算例3计算结果比较 Table 3 Comparison of calculation results of Example 3
| 方法 | 失效概率 | 误差/% | 结构分析次数 |
| 蒙特卡罗法[23] | 0.844 20 | 0 | 100 000 |
| 经典响应面法 | 0.836 41 | 0.922 77 | 196 |
| 线性加权响应面法 | 0.879 58 | 4.199 10 | 21 |
| 单加权响应面法 | 0.838 17 | 0.714 29 | 28 |
| 双加权响应面法 | 0.836 32 | 0.933 43 | 49 |
| 改进加权响应面法 | 0.848 89 | 0.555 56 | 21 |
表选项
4 结论 针对响应面法的计算效率和精度难以平衡兼顾的问题，提出一种基于改进加权响应面的结构可靠度计算方法，经案例验证表明：
1) 提出了综合考虑3个权重因子的权重矩阵构造方法，有效选取了样本点，提高了理想样本点，并弱化了不理想样本点在拟合响应面函数中所起的作用。
2) 在迭代过程中，多次重用已知样本点更新响应面函数，既没有浪费已有样本点中的有用信息，又不会引入劣质样本点而降低响应面的拟合精度，充分利用了已知样本点信息，提高了计算效率。
3) 迭代收敛后，选取权重较大的样本点并采用加权回归方法拟合含交叉项的响应面函数，再次重复使用已知样本点，在不增加结构分析次数的同时提高了计算精度。
参考文献
| [1] | 徐铭阳, 刘林江. 可靠度理论及一次二阶矩法概述[J]. 山西建筑, 2016, 42(8): 67-68. XU M Y, LIU L J. Illustration on reliability theory and first-order second-moment method[J]. Shanxi Architecture, 2016, 42(8): 67-68. DOI:10.3969/j.issn.1009-6825.2016.08.035 (in Chinese) |
| [2] | BREITUNG K. Asymptotic approximations for multinormal integrals[J]. Journal of Engineering Mechanics, 1984, 110(3): 357-366. DOI:10.1061/(ASCE)0733-9399(1984)110:3(357) |
| [3] | 张明. 结构可靠度分析——方法与程序[M]. 北京: 科学出版社, 2009: 68-77. ZHANG M. Structural reliability analysis: Methods and procedures[M]. Beijing: Science Press, 2009: 68-77. (in Chinese) |
| [4] | 吕震宙, 宋淑芳, 李洪双, 等. 结构机构可靠性及可靠性灵敏度分析[M]. 北京: 科学出版社, 2009: 98-104. LV Z Z, SONG S F, LI H S, et al. Reliability and reliability sensitivity analysis of structural mechanism[M]. Beijing: Science Press, 2009: 98-104. (in Chinese) |
| [5] | ROBERT C P, CASELLA G. Monte Carlo statistical methods[M]. 2nd ed.Berlin: Springer, 2004: 79-122. |
| [6] | 赵国藩. 工程结构可靠性理论与应用[M]. 大连: 大连理工大学出版社, 1996: 145-158. ZHAO G F. Reliability theory and its applications for engineering structures[M]. Dalian: Dalian University of Technology Press, 1996: 145-158. (in Chinese) |
| [7] | BUCHER C G, BOURGUND U. A Fast and efficient response surface approach for structural reliability problems[J]. Structural Safety, 1990, 7(1): 57-66. DOI:10.1016/0167-4730(90)90012-E |
| [8] | RAJASHEKHAR M R, ELLINGWOOD B R. A new look at the response surface approach for reliability analysis[J]. Structural Safety, 1993, 12(3): 205-220. DOI:10.1016/0167-4730(93)90003-J |
| [9] | KAYMAZ I, MCMAHON C A. A response surface method based on weighted regression for structural reliability analysis[J]. Probabilistic Engineering Mechanics, 2005, 20(1): 11-17. DOI:10.1016/j.probengmech.2004.05.005 |
| [10] | 赵洁, 吕震宙. 隐式极限状态方程可靠性分析的加权响应面法[J]. 机械强度, 2006, 28(4): 512-516. ZHAO J, LV Z Z. Response surface method for reliability analysis of implicit limit state equation based on weighted regression[J]. Journal of Mechanical Strength, 2006, 28(4): 512-516. DOI:10.3321/j.issn:1001-9669.2006.04.009 (in Chinese) |
| [11] | NGUYEN X A, SELLIER A, DUPRAT F, et al. Adaptive response surface method based on a double weighted regression technique[J]. Probabilistic Engineering Mechanics, 2009, 24(2): 135-143. DOI:10.1016/j.probengmech.2008.04.001 |
| [12] | 李贵杰, 吕震宙, 赵新攀. 基于加权线性响应面的结构可靠性估计的鞍点逼近方法[J]. 机械强度, 2010, 32(6): 917-921. LI G J, LV Z Z, ZHAO X P. Saddlepoint approximation of structural reliability on weighted linear response surface method[J]. Journal of Mechanical Strength, 2010, 32(6): 917-921. (in Chinese) |
| [13] | 钟宏林, 吴剑国, 王恒军. 可靠性分析的双加权响应面法[J]. 浙江工业大学学报, 2010, 38(2): 218-221. ZHONG H L, WU J G, WANG H J. Reliability analysis of response surface method based on a double weighted regression technique[J]. Journal of Zhejiang University of Technology, 2010, 38(2): 218-221. DOI:10.3969/j.issn.1006-4303.2010.02.021 (in Chinese) |
| [14] | 李广博. Fourier正交基神经网络加权响应面法的结构可靠性分析[D]. 长春: 吉林大学, 2014: 37-41. LI G B. Structural reliability analysis based on Fourier orthogonal neural network weighted response surface method[D]. Changchun: Jilin University, 2014: 37-41(in Chinese). |
| [15] | 张吉宁, 郭书祥, 唐承, 等. 基于向量投影取样的改进加权响应面法[J]. 科学通报, 2017, 62(17): 1854-1860. ZHANG J N, GUO S X, TANG C, et al. An improved weighted response surface method based on vector projection sampling[J]. Chinese Science Bulletin, 2017, 62(17): 1854-1860. (in Chinese) |
| [16] | 张学刚. 一种改进响应面法结构可靠度计算方法[J]. 机械强度, 2018, 40(6): 1382-1388. ZHANG X G. Structural reliability analysis based on an improvement of the response surface method[J]. Journal of Mechanical Strength, 2018, 40(6): 1382-1388. (in Chinese) |
| [17] | 贾长安. 基于随机变量相关和失效模式相关的结构可靠性算法研究[D]. 西安: 西安电子科技大学, 2018: 34-41. JIA C A. The study on structural reliability algorithm based on random variable correlation and failure mode correlation[D]. Xi'an: Xidian University, 2018: 34-41(in Chinese). |
| [18] | 侯振兴. 基于响应面法的结构可靠度计算方法研究[D]. 秦皇岛: 燕山大学, 2019: 20-23. HOU Z X. Research on structural reliability calculation method based on response surface method[D]. Qinhuangdao: Yanshan University, 2019: 20-23(in Chinese). |
| [19] | LU C, FENG Y W, FEI C W. Weighted regression-based extremum response surface method for structural dynamic fuzzy reliability analysis[J]. Energies, 2019, 12(9): 1-16. |
| [20] | 洪林雄, 李华聪, 彭凯, 等. 基于高效搜索方法的可靠性分析改进响应面法[J]. 北京航空航天大学学报, 2020, 46(1): 95-102. HONG L X, LI H C, PENG K, et al. Improved response surface method of reliability analysis based on efficient search method[J]. Journal of Beijing University of Aeronautics and Astronautics, 2020, 46(1): 95-102. (in Chinese) |
| [21] | 冯欢欢, 蒋向华. 样本重复使用失效响应曲线分析结构可靠度方法[J]. 航空动力学报, 2013, 28(10): 2228-2234. FENG H H, JIANG X H. Sample points reused failure response curve method for structural reliability analysis[J]. Journal of Aerospace Power, 2013, 28(10): 2228-2234. (in Chinese) |
| [22] | 赵洁. 机械可靠性分析的响应面法研究[D]. 西安: 西北工业大学, 2006: 23-27. ZHAO J. Research on response surface method of mechanical reliability analysis[D]. Xi'an: Northwestern Polytechnical University, 2006: 23-27(in Chinese). |
| [23] | CHOI S K, CANFIELD R A, GRANDHI R. Reliability-based structural design[M]. Berlin: Springer, 2007: 138-141. |
