POINT ESTIMATION FOR STATISTICAL MOMENTS BASED ON THE BIVARIATE DIMENSION-REDUCTION MODEL AND KRIGING APPROXIMATION
-
摘要: 对复杂随机系统进行统计矩分析时，双变量降维近似模型一定程度上可以缓解“维数灾难”。但当系统维数较高时，双变量分量函数较多，计算量仍然较大。为此，该文将降维近似和Kriging代理模型有机结合起来，提出了一类高效、合理的改进点估计法。充分考虑函数逼近和数值积分中积分点的特点，提出了“米”字形的选点策略，并基于此发展了双变量分量函数的Kriging近似模型；将此近似模型用于原函数和矩函数的双变量降维近似模型中双变量分量函数的近似，分别建立了基于原函数近似和矩函数近似的统计矩改进点估计法；通过多个算例对该文提出方法进行了效率和精度的分析。算例分析结果表明：基于“米”字形选点策略的双变量分量函数的Kriging近似具有较高的精度；相比于已有的基于双变量降维近似模型的统计矩点估计法，建议方法仅需较少的结构分析即可达到与已有方法相当的精度，能更好地体现精度和效率的平衡。Abstract: When estimating statistical moments for complex random systems, the bivariate dimension-reduction method can alleviate the curse of dimension to some extent. However, there are many bivariate component functions for high-dimension random systems. It makes the estimation unfeasible. An efficient point estimation for moments is proposed, in which the dimension-reduction model is modified by the Kriging approximation model. Considering the characteristics of function approximation and the abscissas of numerical integration, a "star" shape point-selection strategy is proposed. Based on this point-selection strategy, a Kriging approximation model of the bivariate component function is developed. By replacing each bivariate component function in the bivariate dimension-reduction model for the original function or its moment function with their corresponding Kriging approximations, two modified methods for the moment estimation are presented. The efficiency and accuracy of the proposed methods are verified by several examples. The results show that the Kriging approximation of the bivariate component function based on the "star" shape point-selection strategy has higher accuracy. Correspondingly, the statistical moment estimation of the proposed methods has comparable accuracy with the existing methods, while fewer function evaluations are required.
-
点估计法[1-3]是随机系统常用分析方法，其计算简便有效且可直接与矩方法结合进行结构可靠度分析[4-5]。点估计法中多维数值积分不可避免，往往遭遇维数灾难问题。高效的选点策略与降维近似是目前常用的缓解维数灾难的方法。基于稀疏点集的点估计法[6]是高效选点策略中具有相对较广适用范围与较高精度的方法，但对于交叉项影响较显著的情形精度会有所减弱[6-7]。文献[8]的比较研究表明将多维函数近似为一系列较低维数分量函数组合的降维近似模型在有效缓解维数灾难的同时可保持较好的精度。因此，近年来基于降维近似模型的统计矩估计引起了广泛的关注。对于一个多变量函数，通常可有乘积降维近似[9-10]和加和降维近似[11]两种策略。Rosenblueth[2]首先引入了单变量乘积降维近似模型计算了多变量函数的统计矩。多变量乘积降维近似中，由于具有相同变量的分量函数通过累乘组合在一起，其多维积分并不能转化为低维积分，因此，除单变量乘积降维近似模型外，鲜有乘积降维近似模型用于统计矩估计；然而，随着随机系统维数的增加，单变量乘积降维近似模型并不能保证合理的精度。相比较而言，加和降维近似模型中，分量函数通过加和的方式相组合，其多维积分必定可转化为一系列的低维积分，因此，加和降维近似模型在统计矩估计中应用更为广泛。Zhao等[12]首先将多变量函数的单变量降维近似模型用于统计矩估计；随后Rahman等[13]亦提出了基于单变量降维近似模型的统计矩估计方法，并对单变量降维近似模型进行了详细的理论分析；李洪双等[14]采用了在文献[13]的基础上引入了Nataf变换，拓展了其适用范围。由于仅涉及到多个单变量函数的概率积分，基于单变量降维近似模型的点估计法具有极高的计算效率，然而对于交叉项作用较复杂随机系统，其精度亦不够理想。为改善精度，Xu等[15]进一步提出了针对矩函数进行降维多变量加和降维近似的统计矩估计方法；Huang等[16]和Yu等[17]将矩函数的双变量降维近似分别与Rosenblatt变换和Nataf变换相结合，发展了较为实用的统计矩估计方法。相比于矩函数，原函数更为简单，Fan等[18]提出了针对原函数采用多变量加和降维近似的直接积分法，在不增加计算量的基础上，进一步提高了统计矩估计的精度。不难发现，随着降维近似模型截断维数的增加，模型的精度会逐步改善，但分量函数会逐步增多，进而导致计算量激增。因此，将截断维数取为2是兼顾精度和效率的实用选择[16-18]。然而，即使对于双变量加和降维近似模型，一旦随机系统的维数较高，双变量分量函数将急剧增加，其计算效率并不总是很理想。事实上，双变量分量函数的积分中，各积分点的权重系数是不同的，若允许权重较小的计算点的响应存在一定的误差，则可针对分量函数引入较为成熟的代理模型以改善计算效率。
为此，本文基于函数近似和数值积分中求积节点的特性，提出了Kriging近似模型训练点集确定的“米”字型策略，并基于此建立了双变量分量函数的Kriging近似，将其与现有的基于双变量降维近似模型的统计矩估计方法相结合，建立了更为高效的统计矩点估计法。
1 基于双变量降维近似模型的统计矩点估计法
通常，基于双变量降维近似模型的统计矩计算包括三个步骤：1)将分布类型不同的随机变量X向独立标准正态空间转换，得到新的随机向量U=(U1,U2,···,Un)T和响应函数h(U)；2)利用双变量降维近似模型将多维响应函数分解为二维及以下分量函数的线性组合；3)基于双变量降维近似模型进行响应函数的统计矩计算。
1.1 变量的独立化正态变换
随机向量X中的变量概率信息(如分布类型和相关信息等)常常不同，在随机分析和可靠度分析中常将随机向量X向独立标准正态空间转换，得到新的随机向量U=(U1,U2,···,Un)T。于是，响应函数Z=g(X)可改写为：
Z=g(X)=g(XN,XR)=g(N−1(UN),R−1(UR))≜ (1) 式中：XR为X中已知联合概率分布的变量组成的向量；XN为剩余变量组成的向量，既包括独立变量，亦包括已知边缘密度函数和相关信息的变量；U={UN,UR}。N−1(·)、R−1(·)分别表示Nataf变换[19]与Rosenblatt变换[20]的逆变换。
1.2 基于双变量降维近似模型的统计矩估计
由概率论可知，随机响应函数Z的b阶统计矩矩MZ,b的表达式如下：
{M_{Z,b}} = E[ {{{\left( {h( {{U}} ) - a} \right)}^b}} ] (2) 式中，当b=1时，a=0，MZ,1表示响应函数的均值；当b>1时，a=MZ,1，MZ,b为响应函数的b阶中心矩。
根据降维近似的对象不同，基于双变量降维近似模型的统计矩估计可以分为两类：1) 基于矩函数[h(U)−a]b的双变量降维近似模型的统计矩估计；2) 基于原函数h(U)的双变量降维近似模型的统计矩估计。
1.3 基于矩函数双变量降维近似模型的统计矩估计
若引入[h(U)−a]b的双变量降维近似模型，并进行期望运算，可给出MZ,b的近似计算表达式如下[16]：
\begin{split} {M_{Z,b}} \approx &{\sum\limits_{i < j} {E[ {( {h( {{U_i},{U_j},{{\underline {{u}}}_{ij,{\rm c}}}} ) - a} )} } ^b} - \\ & ( {n - 2} )\sum\limits_{k = 1}^n {E[ {{{( {h( {{U_k},{{\underline {{u}}}_{k,{\rm c}}}} ) - a} )}^b}} ]} + \\ & \frac{{( {n - 1} )( {n - 2} )}}{2}{( {h( {{{{u}}_{\rm c}}} ) - a} )^b}\approx \\ & {\sum\limits_{i < j} \!{\sum\limits_{{l_i} = 1}^r \!{\sum\limits_{{l_j} = 1}^r {{w_{i,{l_i}}}} } {w_{j,{l_j}}}( {h( {{u_{i,{l_i}}},{u_{j,{l_j}}},{{\underline {{u}}}_{ij,{\rm c}}}} )\! -\! a} )^b} } \!- \\ & ( {n - 2} )\sum\limits_{k = 1}^n {\sum\limits_{{l_k}{\rm{ = }}1}^r {{w_{k,{l_k}}}} } {( {h( {{u_{k,{l_k}}},{{\underline {{u}}}_{k,{\rm c}}}} ) - a} )^b} + \\ & \frac{{( {n - 1} )( {n - 2} )}}{2}{( {h( {{{{u}}_{\rm c}}} ) - a} )^b} \end{split} (3) 式中：uc=T−1(μX)为参考点坐标，μX为X的均值；
{{\underline {{u}}} _{ij,c}} 为除变量Ui、Uj以外其余变量对应的参考点坐标；{{\underline {{u}}} _{k,c}} 为除变量Uk以外其余变量对应的参考点坐标；{u_{i,{l_i}}}{\rm{ = }}\sqrt {\rm{2}} {y_{{\rm{GH}},{l_i}}} 、{w_{i,{l_i}}}{\rm{ = }}{{{\omega _{{\rm{GH}},{l_i}}}} / {\sqrt \pi }} 分别是以标准正态密度函数为权函数的高斯求积分公式的求积节点和权系数，其中{y_{{\rm{GH}},{l_i}}} 和{\omega _{{\rm{GH}},{l_i}}} 为一维Gauss-Hermite求积公式的节点及对应的权系数。显然，双变量分量函数
h( {{U_i},{U_j},{{\underline {{u}}}_{ij,c}}} ) 可以有效地考虑Ui与Uj的交互影响。为简便，该方法记为D-2。1.4 基于原函数双变量降维近似模型的统计矩估计
若先引入h(U)的双变量降维近似模型，然后再对[h(U)−a]b进行期望运算，可给出MZ,b的另一种近似计算表达式[18]：
\begin{split} {M_{Z,b}} \approx & E\left[ {\left( {\sum\limits_{i < j} {h( {{U_i},{U_j},{{\underline {{u}}}_{ij,{\rm c}}}} )} \! - \!\left( {n \!- \!2} \right)\sum\limits_{k = 1}^n {h( {{U_k},{{\underline {{u}}}_{k,{\rm c}}}})} \!+} \right.} \right. \\& {\left. {\left. { \frac{{\left( {n - 1} \right)\left( {n - 2} \right)}}{2}h\left( {{{{u}}_{\rm c}}} \right)} \right) - a} \right]^b} \approx \\& \sum\limits_{{l_1} = 1}^r {\sum\limits_{{l_2} = 1}^r \cdots } \sum\limits_{{l_n} = 1}^r {\prod\limits_{q = 1}^n {{w_{q,{l_q}}}} \cdot \Bigg[ {\sum\limits_{i < j} {h( {{u_{i,{l_i}}},{u_{j,{l_j}}},{{\underline {{u}}}_{ij,{\rm c}}}} )} } \Bigg.} - \\& \left( {n - 2} \right)\sum\limits_{k = 1}^n {h( {{u_{k,{l_k}}},{{\underline {{u}}}_{k,{\rm c}}}} }) \! +\! \frac{{\left( {n - 1} \right)( {n \!- \!2} )}}{2}h( {{{{u}}_{\rm c}}} ) {\Bigg. { \!-\! a} \Bigg]^b} \end{split} (4) 为简便，该方法记为N-2。
2 改进的双变量降维近似统计矩估计
在双变量降维近似统计矩估计中，无论是D-2方法还是N-2方法，其计算量主要体现于双变量分量函数在积分节点处的函数值计算，对于隐式的g(X)，其函数值需通过结构分析获得。实用中，通常取r=7，则任一双变量分量函数均涉及7×7=49次的结构分析。若能在保持计算精度相当的前提下减少原函数的调用次数，则可以有效改善计算效率，较好地兼顾精度和效率的平衡。
事实上，在双变量分量函数的49个求积节点中，权系数的差异较为明显，权系数较小节点的函数值即使存在着一定的误差，总体计算结果的误差仍是可以接受的。为此，本文将任一双变量分量函数的49个求积节点分为计算节点和近似节点，其中计算节点的函数值通过函数调用或结构分析获得，近似节点的函数值则由基于计算节点确定的双变量分量函数的近似模型计算得到，此处近似模型采用Kriging模型。由于近似节点的函数值不再涉及原函数调用或结构分析，将有效改善计算效率。
2.1 Kriging近似模型
对于一个n维随机系统U={U1, U2, ···,Un}，Kriging模型
\tilde h({{U}}) 可表示为[21]：\tilde h({{U}}) = f({{U}}) + m({{U}})\;\qquad\qquad\qquad (5) 式中：
f({{U}}) 为回归函数；m({{U}}) 为零均值、σ方差的高斯随机过程。给定训练集合u0={u(1), u(2),···, u(k)}及对应响应集合z0={z1), z(2),···, z(k) }，该高斯过程中任意两点z(i)和z(j)的协方差可表示为：{\rm{cov}} ( {{{z}}({u^{(i)}}),{{z}}({u^{(j)}})} ) = {\sigma ^2}{{{R}}_{{\theta}} }( {{{{u}}^{(i)}},{{{u}}^{(j)}}} ) (6) 式中，
{{R}}_{{\theta}} ^{}\left( \cdot \right) 是带参数θ的相关函数方程，文中采用如下形式：{{{R}}_{{\theta}} }( {{{{u}}^{(i)}},{{{u}}^{(j)}}} ){\rm{ = }}\exp \left( { - \sum\limits_{l = 1}^n {{\theta _l}\left| {{{{u}}_l}^{(i)}\left. { - {{{u}}_l}^{(j)}} \right|} \right.} } \right) (7) 上述参数θ和σ可由最大似然法估计，其表达式为：
{{\theta}} {\rm{ = }}( {{{{F}}^{\rm T}}{{{R}}_{{\theta}} }{{F}}} ){{{F}}^{\rm T}}{{{R}}_{{\theta}} ^{ - 1}}{{{z}}_0}\quad\qquad\;\; (8) \sigma = \frac{1}{k}{\left( {{{{z}}_0} - F{{\theta}} } \right)^{\rm T}}{R_{{\theta}} ^{ - 1}}\left( {{{{z}}_0} - F{{\theta}} } \right) (9) 其中，F={f(u(1)), f(u(2)),···, f(u(k)) }。
由此，可得
h({{U}}) 的Kriging模型\tilde h({{U}}) 为：\tilde h({{U}}) = f({{U}}) + {\left( {t({{U}},{{{u}}_0})} \right)^{\rm T}}{\left( {{R}} \right)^{ - 1}}({{{z}}_0} - {{H}}({{{u}}_0})) (10) 其中：
{{H}}({{{u}}_0}) = {[h({{{u}}^{(1)}}),h({{{u}}^{(2)}}),\cdots,h({{{u}}^{(k)}})]^{\rm T}}\qquad\qquad (11) {{R}} \!\!=\! \!\left[ \!\!\!{\begin{array}{*{20}{c}} {R_{{\theta}} ^{}({{{u}}^{(1)}},{{{u}}^{(1)}})}&{\cdots}&{R_{{\theta}} ^{}({{{u}}^{(1)}},{{{u}}^{(k)}})} \\ \vdots &{R_{{\theta}} ^{}({{{u}}^{(i)}},{{{u}}^{(i)}})}& \vdots \\ {R_{{\theta}} ^{}({{{u}}^{(k)}},{{{u}}^{(1)}})}&{\cdots}&{R_{{\theta}} ^{}({{{u}}^{(k)}},{{{u}}^{(k)}})} \end{array}} \!\!\!\right] (12) t( {{{U}},{{{u}}_0}} ) \!= \!{[ \!{{{{R}}_{{\theta}} }( {{{U}},{{{u}}^{( 1 )}}} ), {{{R}}_{{\theta}} }( {{{U}},{{{u}}^{( 2 )}}} ), \cdots ,{{{R}}_{{\theta}} }( {{{U}},{{{u}}^{( k )}}} )} \!]^{\rm T}} (13) 2.2 改进的双变量降维近似统计矩估计
2.2.1 双变量分量函数的Kriging模型
1)计算节点的选取策略
图1(a)给出了以标准正态密度函数为权函数的二维高斯求积分公式的49个积分节点，记为集合DOE。理论上，定义域内任意点皆可以作为Kriging模型的训练点集，即计算节点。然而，统计矩估计中仅需要双变量分量函数在49个积分节点的函数值，因此，计算节点可在此49个节点中选取，且应遵循如下准则：
① 具有较大权系数的计算节点函数值误差对整体积分的影响往往较大，因此计算节点优先选择权系数较大的积分节点，换言之，计算节点应尽量靠近中心区域；
② 函数近似中，内插精度往往高于外插精度，因此计算节点应包含一定的外围节点；
③ 当涉及对称型分布的随机变量，其参考点坐标分量为0，二维积分节点中的轴点和一维积分节点重合，为充分利用这些节点，计算节点宜优先选择轴点。
2)双变量分量函数的Kriging近似
首先，针对集合DOEY中所有的计算节点
{{{u}}_0} ，通过调用函数或结构分析确定双变量函数h(Ui,Uj,uij,c)在计算节点处的值{z_0} ；然后，由式(10)可得到h(Ui,Uj,uij,c)的Kriging近似模型为：\begin{split} & \tilde h( {{U_i},{U_j},{{\underline {{u}}}_{ij,c}}} ){\rm{ = }}f({U_i},{U_j}) + \\&\qquad {( {t( {{U_i},{U_j},{{{u}}_0}} )} )^{\rm T}}{( {{R}} )^{ - 1}}( {{{{z}}_0} - {{H}}( {{{{u}}_0}} )} ) \end{split} (14) 通常，回归函数
f\left( \cdot \right) 取为0阶多项式函数，即f( {{U_i},{U_j}} ) = 1 。2.2.2 改进的双变量降维近似统计矩估计
相应于常规的基于双变量降维近似模型的统计矩估计，结合双变量分量函数Kriging近似的改进统计矩估计亦分为两类：基于矩函数双变量降维和Kriging近似的统计矩估计、基于原函数双变量降维和Kriging近似的统计矩估计。
基于上述考虑，本文给出的计算节点选取方案如图1(b)所示。为简便，可将计算节点的集合记为DOEY，其形状呈“米”字形，故称为“米”字形策略。
1)基于矩函数双变量降维和Kriging近似的统计矩估计
将式(14)得到的双变量分量函数的Kriging近似代入式(3)，可得到改进的基于矩函数双变量降维近似模型的统计矩估计：
\begin{split} & {M_{Z,b}} \approx {\sum\limits_{i < j} {\sum\limits_{{l_i} = 1}^r {\sum\limits_{{l_j} = 1}^r {{w_{i,{l_i}}}} } {w_{j,{l_j}}}( {\tilde h( {{u_{i,{l_i}}},{u_{j,{l_j}}},{{\underline {{u}}}_{ij,{\rm c}}}} ) - a} )^b} } - \\& \;\;\;\;\;\;\;( {n - 2} )\sum\limits_{k = 1}^n {\sum\limits_{{l_k}{\rm{ = }}1}^r {{w_{k,{l_k}}}} } {( {h( {{u_{k,{l_k}}},{{\underline {{u}}}_{k,{\rm c}}}} ) - a} )^b} + \\& \;\;\;\;\;\;\;\frac{{( {n - 1} )( {n - 2} )}}{2}{( {h( {{{{u}}_{\rm c}}} ) - a} )^b} \\[-16pt] \end{split} (15) 为简便，该方法记为米D-2。
2)基于原函数双变量降维和Kriging近似的统计矩估计
将式(14)得到的双变量分量函数的Kriging近似函数代入式(4)中，可得到改进的基于原函数双变量降维近似模型的统计矩估计，其表达式如下：
\begin{split} & {M_{Z,b}} \approx \sum\limits_{{l_1} = 1}^r {\sum\limits_{{l_2} = 1}^r \cdots } \sum\limits_{{l_n} = 1}^r {\prod\limits_{q = 1}^n {{w_{q,{l_q}}}} \cdot \left[ {\sum\limits_{i < j} {\tilde h( {{u_{i,{l_i}}},{u_{j,{l_j}}},{{\underline {{u}}}_{ij,{\rm c}}}} )} - } \right.} \\& \;\;\quad ( {n - 2} )\sum\limits_{k = 1}^n {h( {{u_{k,{l_k}}},{{\underline {{u}}}_{k,{\rm c}}}} )} + \frac{{( {n - 1} )( {n - 2} )}}{2}h( {{{{u}}_{\rm c}}} ){\left. { - a} \right]^b} \end{split} (16) 为简便，该方法记为米N-2。
不难发现，与D-2、N-2方法相比，米D-2和米N-2方法亦充分考虑了交叉项的影响，只不过是真实双变量分量函数被合适的Kriging函数近似代替。
2.3 算法步骤
综合上述过程，改进的双变量降维近似统计矩估计实现过程如下：
1)根据已知概率分布信息，由式(1)将随机向量X转化为标准正态向量U。
2)确定各双变量分量函数在计算节点处的函数值，并基于此由式(14)确定各双变量分量函数的Kriging近似。
3)确定各单变量分量函数在求积节点处的函数值，由各双变量分量函数的Kriging近似确定其在二维积分节点的函数值，将上述函数值以及参考点处的函数值代入式(15)或式(16)即可得到响应Z的b阶统计矩。
3 算例分析
为了验证改进的双变量降维近似统计矩估计方法的精度和效率，此处通过三个算例(1个数值算例和2个工程算例)开展了详细的分析。为评估不同方法的精度，定义相对误差如下：
\varepsilon {\rm{ = }}\frac{{\left| {{\text{各方法计算值}} - {\text{标准解}}} \right|}}{{{\text{标准解}}}} \times 100 (17) 其中，标准解为多维Gauss-Hermite数值积分迭代收敛的结果或3×107个样本点的蒙特卡罗法的计算结果。
算例1. 数值算例
假设响应Z的表达式为：
g\left( {{X}} \right) = 200 - {X_1^3} - 3.5{X_2^3} + {X_3^3} + {X_1}{X_2} (18) 式中：X1为对数正态分布变量，均值为22、标准差为2；X2为正态分布变量，均值为10、标准差为0.9；X3为极值I型变量，均值为2、标准差为0.6。
将X变换为独立标准正态向量可得：
\begin{split} Z =& h\left({{U}} \right) = 200 - {( {{F_{X_1}^{ - 1}}( {\varPhi ( {{U_1}} )} )} )^3} - \\& 3.5{( {0.9{U_2} + 10} )^3} + {( {{F_{{X_3}}^{ - 1}}( {\varPhi ( {{U_3}} )} )} )^6} + \\& {F_{{X_1}}^{ - 1}}( {\varPhi ( {{U_1}} )} )( {0.9{U_2} + 10} ) \end{split} (19) 式中：
\varPhi \left( \cdot \right) 为标准正态变量的累积分布函数；{F_{{X_i}}^{ - 1}}\left( \cdot \right) 为Xi的累积分布函数{F_{{X_i}}}\left( \cdot \right) 的反函数。对于上述三变量函数，参考点为uc=T−1(μX)={0.0454,0, 0.1773}，即u12,c=0.1773、u13,c=0、u23,c=0.0454，其双变量降维近似模型包含3个双变量分量函数。对于双变量分量函数h(U1,U2,u12,c)，根据图1(b)所示的米字形计算节点及其函数值，由式(8)和式(9)可确定参数θ=−0.6527和σ=6.5159×103，然后由式(14)确定的Kriging近似函数
\tilde h( {{U_1},{U_2},{{\underline {{u}}}_{12,{\rm c}}}} ) 计算h(U1,U2,u12,c)在近似节点的函数值。图2括号中的数值为Kriging近似函数在近似节点处的值与相应的实际函数值的相对误差。不难发现，对于大多数的近似节点，Kriging近似模型具有非常高的精度；即使个别近似节点的误差较大，但由于其权重较小，在统计矩的贡献很小。因此，基于米字形选点方案确定的Kriging模型是合理、有效的。类似地，可确定另外两个双变量分量函数的Kriging近似函数\tilde h( {{U_2},{U_3},{{\underline {{u}}}_{23,{\rm c}}}} ) 和\tilde h( {{U_1},{U_3},{{\underline {{u}}}_{13,{\rm c}}}} ) 。![]()
图 2 h(U1,U2,u12,c)的选点方案和Kriging近似模型的精度Figure 2. Scheme for selecting points and accuracy of Kriging approximation model of h(U1,U2,u12,c)分别采用D-2、N-2和米D-2、米N-2计算g(X)的前四阶统计矩，标准解为多维Gauss-Hermite数值积分迭代收敛的结果，结果如表1所示。从表1可以发现：1) 4种方法均具有很高的精度，D-2和N-2两者精度相当，米D-2和米N-2两者精度大体相当但米N-2略优，且后两者较前两者精度略有下降，体现了Kriging近似误差的影响，但结果的高精度表明了引入Kriging近似的可行性和有效性，此外，从低阶矩到高阶矩，各方法误差呈增大趋势；2) 函数调用次数上，D-2和N-2两者相同，米D-2和米N-2两者相同，且由于引入了Kriging近似，后两者的函数调用次数较前两者显著减少，由154次降低至82。具体而言，对于双变量分量函数h(U1,U2,u12,c)，需要计算
h( {{u_{1,{l_1}}},{u_{2,{l_2}}}, {{\underline {{u}}}_{12,{\rm c}}}} ) (l1=1,2,···,7; l2=1,2,···,7)，共49个函数值，但是当u2,l2=0=u13,c时，h(u1,l1,0,u12,c)= h(u1,l1, u13,c,u12,c)，即图2的所有节点中，U1轴上的节点与U1的单变量分量函数的节点相同，不需重复计算。因此，在D-2和N-2中，h(U1,U2,u12,c)需调用函数49−7=42次；在米D-2和米N-2中，h(U1,U2,u12,c)需调用函数25−7=18次。不难发现，若双变量分量函数中某一变量的参考点坐标为0(通常对应着具有对称分布的随机变量)，另一变量的7个轴点可不用计算；若双变量分量函数中两个变量的参考点坐标均为0，则两个变量的13个轴点可不用计算；此外，对于参考点坐标为0的变量的单变量分量函数，7个节点的中间节点坐标与参考点相同，亦不需重复计算。于是，D-2和N-2的函数调用次数为1+7+7+6+42+42+49=154；米D-2和米N-2的函数调用次数为1+7+7+6+18+18+2=82。表 1 统计矩计算结果对比Table 1. Comparison of statistical moments方法 μZ mZ,2 函数调用
次数N数值/
(×104)相对误差
ε/(%)数值/
(×106)相对误差
ε/(%)标准解 −1.4069 − 9.9240 − 73 D-2 −1.4069 0 9.9240 0.00 154 N-2 −1.4069 0 9.9240 0.00 154 米D-2 −1.4069 0 9.9229 0.01 82 米N-2 −1.4069 0 9.9229 0.01 82 方法 mZ,3 mZ,4 函数调用
次数N数值/
(×1010)相对误差
ε/(%)数值/
(×1014)相对误差
ε/(%)标准解 −2.3721 − 4.0460 − 73 D-2 −2.3721 0.00 4.0465 0.01 154 N-2 −2.3721 0.00 4.0460 0.00 154 米D-2 −2.3756 0.15 4.0542 0.20 82 米N-2 −2.3734 0.10 4.0529 0.17 82 算例2. 工程算例
考察图3所示承受轴向压缩载荷F的环形截面柱，其功能函数为[22]：
\begin{split} Z = & g\left( {E,F,L,D,T} \right) {\rm{= }} \\& \frac{{{\pi ^2}E}}{{{L^2}}}\left\{ {\frac{\pi }{{64}}[ {{{\left( {D + T} \right)}^4}} ] - {D^4}} \right\} - F \end{split} (20) 式中：E、F、L、D、T分别是材料弹性模量、轴荷载、柱高、内径和厚度；E、F为对数正态变量，均值为2.1×1011 Pa和6000 N，标准差为0.4×1011 Pa和400 N；D、T、L是独立的正态变量，均值分别为0.03、0.006和2.5，标准差分别为0.001、0.0004和0.09。
分别采用D-2、N-2和米D-2、米N-2 四种方法计算式(20)所示功能函数的前四阶矩，标准解为多维Gauss-Hermite数值积分迭代收敛的结果，结果如表2所示。从表2可以发现：1)高阶矩计算结果中，D-2较N-2精度偏低，米D-2较米N-2精度偏低，米D-2和米N-2分别较D-2、N-2略高，其原因在于建议两方法中降维近似误差和Kriging近似误差的相互抵消；2)函数调用次数上，D-2和N-2相同，米D-2和米N-2相同，且由于引入了Kriging近似后两者的计算效率提升1倍多，效果显著。
表 2 统计矩计算结果对比Table 2. Comparison of statistical moments方法 μZ mZ,2 函数调用
次数N数值/
(×103)相对误差
ε/(%)数值/
(×107)相对误差
ε/(%)标准解 8.2704 − 1.1957 − 75 D-2 8.2704 0 1.1954 0.03 449 N-2 8.2704 0 1.1954 0.03 449 米D-2 8.2703 0 1.1955 0.02 209 米N-2 8.2703 0 1.1955 0.02 209 方法 mZ,3 mZ,4 函数调用
次数N数值/
(×1010)相对误差
ε/(%)数值/
(×1014)相对误差
ε/(%)标准解 2.9301 − 5.6011 − 75 D-2 2.8603 2.38 5.1838 7.45 449 N-2 2.9005 1.01 5.5069 1.68 449 米D-2 2.8990 1.06 5.2897 5.56 209 米N-2 2.9328 0.09 5.6234 0.40 209 算例3. 工程算例
考察如下的磁盘极限承载力[23]：
\begin{split} Z = & \left( {MUF,UTS,\delta ,N,R,{R_0}} \right) = \\& \sqrt {\frac{{\left( {MUF} \right)\left( {UTS} \right)}}{{3 \cdot 385.82 \cdot \delta \left[ {N\dfrac{{2\pi }}{{60}}} \right]( {{R^3} - {R_0^3}} )\left( {R - {R_0}} \right)}}} \end{split} (21) 式中：MUF、UTS、δ、N、R、R0 分别是材料利用率、极限拉伸强度、密度、转子转速、外半径和内半径；上述各变量均为独立正态分布变量，其均值分别0.8358、22 000、0.29、21 000、24和8，标准差分别为0.2、5000、0.005、1000、0.5和0.3。
分别采用D-2、N-2和米D-2、米N-2 四种方法计算式(21)所示功能函数的前四阶矩，标准解为3×107个样本点的蒙特卡罗法的计算结果，结果如表3所示。从表3可以发现：1) D-2和N-2两者精度相当，米D-2较米N-2略低，且前两者在二、三阶矩精度较后两者略好，但在四阶矩精度略差，因此，建议方法中的降维近似误差和Kriging近似误差的叠加效应和抵消效应对于各阶矩并不总是保持一致；2)建议方法的函数调用次数从现有方法的577次降为217，计算效率显著提高。
表 3 统计矩矩计算结果对比Table 3. Comparison of statistical moments方法 μZ mZ,2 函数调用
次数N数值/
(×10−5)相对误差
ε/(%)数值/
(×10−11)相对误差
ε/(%)标准解 2.2985 − 1.0817 − 3×107 D-2 2.2985 0 1.0819 0.018 577 N-2 2.2985 0 1.0819 0.018 577 米D-2 2.2985 0 1.0813 0.037 217 米N-2 2.2985 0 1.0813 0.037 217 方法 mZ,3 mZ,4 函数调用
次数N数值/
(×10−18)相对误差
ε/(%)数值/
(×10−22)相对误差
ε/(%)标准解 −3.2444 − 3.8767 − 3×107 D-2 −3.3335 2.75 3.8335 1.11 577 N-2 −3.3335 2.75 3.8335 1.11 577 米D-2 −3.3728 3.96 3.8547 0.57 217 米N-2 −3.3471 3.17 3.8981 0.55 217 进一步对比上述3个算例，随机系统维数从3上升至6时，现有方法的函数调用次数从154上升至577，而本文建议方法则从82上升至217。事实上，由算例1的效率分析可以发现基于双变量降维近似模型的统计矩估计方法的函数调用次数与参考点坐标有着较大的关系。考察n变量函数，对于参考点坐标均为0和参考点坐标均不为0两种极端情况，N-2(或D-2)与米N-2(或米D-2)的函数调用次数表达式为：
{N_1} = 49\frac{{n\left( {n - 1} \right)}}{2} + 7n + 1\qquad\; (22) {N_2} = 25\frac{{n\left( {n - 1} \right)}}{2} + 7n + 1\qquad\;\; (23) {N_3} = \left( {49 - 13} \right)\frac{{n\left( {n - 1} \right)}}{2} + 6n + 1 (24) {N_4}{\rm{ = }}\left( {25 - 13} \right)\frac{{n\left( {n - 1} \right)}}{2} + 6n + 1\;\; (25) 式中：N1、N2分别为参考点坐标均不为0时N-2与米N-2的函数调用次数；N3、N4分别为参考点坐标均为0时N-2与米N-2的函数调用次数。图4给出了2种情况下各方法函数调用次数随维数的变化规律，虽然米N-2的函数调用次数与N-2具有类似的增长规律，但总的分析次数显著减少。若定义：
{\gamma _1} = \frac{{{N_2}}}{{{N_1}}} (26) {\gamma _2} = \frac{{{N_4}}}{{{N_3}}} (27) 显然，γ1和γ2分别反映了两种极端情况下米N-2与N-2的函数调用次数的比值。图5给出了γ1和γ2随维数的变化规律，不难发现γ1约为52%，γ2约为35%。
易知，对于参考点坐标不全为0的一般情况，米N-2与N-2的函数调用次数的比值应在35%~52%。换言之，相比于已有方法，文中建议方法可提升1倍的计算效率，能更好地兼顾精度和效率。
4 结论
为了进一步提高基于双变量降维近似模型的统计矩估计方法的效率，本文发展了基于“米”字形选点策略的双变量函数的Kriging近似模型，并结合现有的矩函数或原函数双变量降维近似模型，提出了更好地兼顾精度和效率的改进统计矩估计方法。文中分析结果表明：
(1)本文建议方法(米N-2和米D-2)的计算效率相比于现有方法(N-2和D-2)显著提高，函数调用次数约为现有方法的50%左右；
(2)建议方法的误差主要包括降维近似误差和Kriging近似误差，两种误差有时叠加、有时抵消；总体而言，建议方法与现有方法的精度大体相当；
(3)建议的米N-2和米D-2中，米N-2具有相对更好的精度，可优先采用。
-
![]()
图 2 h(U1,U2,u12,c)的选点方案和Kriging近似模型的精度
Figure 2. Scheme for selecting points and accuracy of Kriging approximation model of h(U1,U2,u12,c)
表 1 统计矩计算结果对比
Table 1 Comparison of statistical moments
方法 μZ mZ,2 函数调用
次数N数值/
(×104)相对误差
ε/(%)数值/
(×106)相对误差
ε/(%)标准解 −1.4069 − 9.9240 − 73 D-2 −1.4069 0 9.9240 0.00 154 N-2 −1.4069 0 9.9240 0.00 154 米D-2 −1.4069 0 9.9229 0.01 82 米N-2 −1.4069 0 9.9229 0.01 82 方法 mZ,3 mZ,4 函数调用
次数N数值/
(×1010)相对误差
ε/(%)数值/
(×1014)相对误差
ε/(%)标准解 −2.3721 − 4.0460 − 73 D-2 −2.3721 0.00 4.0465 0.01 154 N-2 −2.3721 0.00 4.0460 0.00 154 米D-2 −2.3756 0.15 4.0542 0.20 82 米N-2 −2.3734 0.10 4.0529 0.17 82 
下载: 导出CSV
表 2 统计矩计算结果对比
Table 2 Comparison of statistical moments
方法 μZ mZ,2 函数调用
次数N数值/
(×103)相对误差
ε/(%)数值/
(×107)相对误差
ε/(%)标准解 8.2704 − 1.1957 − 75 D-2 8.2704 0 1.1954 0.03 449 N-2 8.2704 0 1.1954 0.03 449 米D-2 8.2703 0 1.1955 0.02 209 米N-2 8.2703 0 1.1955 0.02 209 方法 mZ,3 mZ,4 函数调用
次数N数值/
(×1010)相对误差
ε/(%)数值/
(×1014)相对误差
ε/(%)标准解 2.9301 − 5.6011 − 75 D-2 2.8603 2.38 5.1838 7.45 449 N-2 2.9005 1.01 5.5069 1.68 449 米D-2 2.8990 1.06 5.2897 5.56 209 米N-2 2.9328 0.09 5.6234 0.40 209
下载: 导出CSV
表 3 统计矩矩计算结果对比
Table 3 Comparison of statistical moments
方法 μZ mZ,2 函数调用
次数N数值/
(×10−5)相对误差
ε/(%)数值/
(×10−11)相对误差
ε/(%)标准解 2.2985 − 1.0817 − 3×107 D-2 2.2985 0 1.0819 0.018 577 N-2 2.2985 0 1.0819 0.018 577 米D-2 2.2985 0 1.0813 0.037 217 米N-2 2.2985 0 1.0813 0.037 217 方法 mZ,3 mZ,4 函数调用
次数N数值/
(×10−18)相对误差
ε/(%)数值/
(×10−22)相对误差
ε/(%)标准解 −3.2444 − 3.8767 − 3×107 D-2 −3.3335 2.75 3.8335 1.11 577 N-2 −3.3335 2.75 3.8335 1.11 577 米D-2 −3.3728 3.96 3.8547 0.57 217 米N-2 −3.3471 3.17 3.8981 0.55 217
下载: 导出CSV
-
[1] Evans D H. An application of numerical integration techniques to statistical tolerancing [J]. Technimetrics, 1967, 9(3): 441 − 456.
[2] Rosenblueth E. Point estimates for probability moments [J]. Proceedings of the National Academy of Sciences of the United States of America, 1975, 72(10): 3812 − 3814. doi: 10.1073/pnas.72.10.3812
[3] Rosenblueth E. Two-point estimates in probabilities [J]. Applied Mathematical Modeling, 1981, 5(2): 329 − 335.
[4] Zhao Y G, Ono T. Moment methods for structural reliability [J]. Structural Safety, 2001, 23(1): 47 − 75. doi: 10.1016/S0167-4730(00)00027-8
[5] Zhao Y G, Ang A H-S. System reliability assessment by method of moments [J]. Journal of Structural Engineering, 2003, 129(10): 1341 − 1349. doi: 10.1061/(ASCE)0733-9445(2003)129:10(1341)
[6] He J, Gao S B, Gong J H. Gong, A sparse grid stochastic collocation method for structural reliability analysis [J]. Structural Safety, 2014, 51: 29 − 34. doi: 10.1016/j.strusafe.2014.06.003
[7] Zhang W B, Kang Z. Robust shape and topology optimization considering geometric uncertainties with stochastic level set perturbation [J]. International Journal for Numerical Methods in Engineering, 2016, 110(1): 31 − 56.
[8] 范文亮, 李正良, 王承启. 多变量函数统计矩点估计法的性能比较[J]. 工程力学, 2012, 29(11): 1 − 11. doi: 10.6052/j.issn.1000-4750.2010.09.0698 Fan Wenliang, Li Zhengliang, Wang Chengqi. Comparison of point estimate methods for probability moments of multivariate function [J]. Engineering Mechanics, 2012, 29(11): 1 − 11. (in Chinese) doi: 10.6052/j.issn.1000-4750.2010.09.0698
[9] Tunga M A, Demiralp M. A Factorized high dimensional model representation on the partitioned random discrete data [J]. Applied Numerical Analysis & Computational Mathematics, 2004, 1(1): 231 − 241.
[10] Tunga M A , Demiralp M. A factorized high dimensional model representation on the nodes of a finite hyperprismatic regular grid [J]. Applied Mathematics & Computation, 2005, 164(3): 865 − 883.
[11] Zhang X, Pandey M D. Structural reliability analysis based on the concepts of entropy, fractional moment and dimensional reduction method [J]. Structural Safety, 2013, 43: 28 − 40. doi: 10.1016/j.strusafe.2013.03.001
[12] Zhao Y G, Ono T. New point estimates for probability moments [J]. Journal of Engineering Mechanics, 2000, 126(4): 433 − 436. doi: 10.1061/(ASCE)0733-9399(2000)126:4(433)
[13] Rahman S, Xu H. A univariate dimension-reduction method for multi-dimensional integration in stochastic mechanics [J]. Probabilistic Engineering Mechanics, 2004, 19(4): 393 − 408. doi: 10.1016/j.probengmech.2004.04.003
[14] 李洪双, 吕震宙, 袁修开. 基于Nataf变换的点估计法[J]. 科学通报, 2008(6): 13 − 18. Li Hongshuang, Lü Zhenzhou, Yuan Xiukai. Nataf transformation based point estimate method [J]. Chinese Science Bulletin, 2008(6): 13 − 18. (in Chinese)
[15] Xu H, Rahman S. A generalized dimension-reduction method for multidimensional integration in stochastic mechanics [J]. International Journal for Numerical Methods in Engineering, 2004, 61(12): 1992 − 2019. doi: 10.1002/nme.1135
[16] Huang X Z, Zhang Y M. Reliability–sensitivity analysis using dimension reduction methods and saddlepoint approximations [J]. International Journal for Numerical Methods in Engineering, 2012, 93(8): 857 − 886.
[17] Yu X H, Lu D G. An advanced point estimate method for uncertainty and sensitivity analysis using nataf transformation and dimension-reduction integration [J]. Numerical Methods for Reliability and Safety Assessment, 2015: 215 − 239.
[18] Fan W L, Liu R Y, Ang A H-S, et al. A new point estimation method for statistical moments based on dimension-reduction method and direct numerical integration [J]. Applied Mathematical Modelling, 2018, 62: 664 − 679. doi: 10.1016/j.apm.2018.06.022
[19] Der Kiureghian A, Liu P L. Structural reliability under incomplete probability information [J]. Journal of Engineering Mechanics, 1986, 112(1): 85 − 104. doi: 10.1061/(ASCE)0733-9399(1986)112:1(85)
[20] Rosenblatt M. Remarks on a multivariate transformation [J]. Annals of Mathematical Statistics, 1952, 23(3): 470 − 472. doi: 10.1214/aoms/1177729394
[21] Wen Z X, Pei H Q, Liu H. et al. A sequential Kriging reliability analysis method with characteristics of adaptive sampling regions and parallelizability [J]. Reliability Engineering & System Safety, 2016, 153: 170 − 179.
[22] Chowdhury R, Rao B N, Prasad A M. High-dimensional model representation for structural reliability analysis [J]. Communications in Numerical Methods in Engineering, 2009, 25(4): 301 − 337. doi: 10.1002/cnm.1118
[23] Wang L P, Beeson D, Wiggs G. Efficient and accurate point estimate method for moments and probability distribution estimation [C]// 10th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference. New York: American Institute of Aeronautics and Astronautics, 2004: 693 − 703.
-
期刊类型引用(3)
1. 刘丞，范文亮，余书君，李正良. 基于主动学习Kriging模型的改进一次可靠度方法. 工程力学. 2024(02): 35-42 . 
本站查看
2. 洪兆溪，谭建荣，何利力，胡炳涛，张志峰，宋秀菊，冯毅雄. 高维混合不确定条件下基于降维积分的产品关键结构可靠性优化设计. 计算机集成制造系统. 2024(12): 4152-4167 . 百度学术
3. 韩忠皓，张德权，杨美德，赵海文. 基于Kriging模型的双变量降维方法. 河北工业大学学报. 2023(06): 46-52+59 . 百度学术
其他类型引用(4)
计量
- 文章访问数: 401
- HTML全文浏览量: 114
- PDF下载量: 42
- 被引次数: 7
