基于改进统计矩点估计法和最大熵原理的结构整体可靠度分析

王 涛，李正良,2，范文亮,2

(1. 重庆大学土木工程学院，重庆 400045；2. 山地城镇建设与新技术教育部重点实验室(重庆大学)，重庆 400045)

由于材料、荷载不可避免的随机性，结构可靠度分析成为结构设计的基础和重要保障。经过长期的发展，结构构件可靠度方法已形成了数值近似法(如一次二阶矩法[1])、随机模拟法(如Monte Carlo 法[2])、代理模型法(如响应面法[3]) 和数值积分法(如矩方法[4]) 四大类方法，其基本理论已经相对成熟且能够应用于工程实际[1]。然而，结构整体可靠度分析仍在发展阶段。大体上，结构整体可靠度方法可分为三类，即Monte Carlo 法[5]、失效模式识别法[5]以及等价描述方法[6-8]。Monte Carlo 法由于其计算昂贵难以应用于大型复杂结构系统，常作为一种校准法；而经典的失效模式识别法则存在组合爆炸和相关失效问题[6]；对于复杂工程结构进行整体可靠度分析，Monte Carlo 法计算成本难以接受，失效模式识别法应用困难，相比而言，由于等价描述方法将结构整体可靠度用一个等价的功能函数描述，有效地避免了组合爆炸与相关失效问题，形成了一种较为实用可行的思路。

功能函数经过等价描述后，其函数形式变得较为复杂，由于基于验算点的方法处理函数形式复杂的功能函数会存在多验算点问题，采用此类方法进行整体可靠度分析可能致使其计算精度不理想[9]。将等价描述后的功能函数与不依赖验算点的方法结合进行结构整体可靠度分析是一种行之可效的思路。矩方法是获得结构可靠度的一种有效途径，既不依赖于验算点，也无需复杂的积分计算和求导运算，因而得到广泛应用[10-15]。其主要思路是先获取功能函数的统计矩，再而基于统计矩信息求解功能函数的概率密度函数，最后通过对功能函数的概率密度函数积分获取失效概率或可靠指标。

功能函数的统计矩估计通常通过数值积分方法实现，其主要包括降维近似统计矩点估计法[10-13]、稀疏网格配点法[14]以及容积积分法[15]。其中，研究者对降维近似统计矩点估计法进行了大量的研究，通过引入不同的降维近似模型，分别发展了单变量、双变量和三变量降维近似统计矩点估计法[10-13]。对于结构整体可靠度问题，等价功能函数蕴含了多个失效模式的相关性，变量间的交互影响一般较为复杂，采用单变量与双变量降维近似统计矩点估计法往往不能达到满意的精度[13]。研究发现，三变量降维近似统计矩点估计法能够较为准确地考虑结构整体可靠度问题中各随机变量的交互影响，适用于结构整体可靠度的评估[13]。但三变量降维近似统计矩点估计法所需的结构分析次数较多，如何在保证精度的前提下改善其计算效率有待进一步提高。

一般而言，在获取功能函数的统计矩后，采用Pearson 分布[16]、立方正态[4]、鞍点近似[13]、移位广义对数分布[17]以及最大熵原理[18]等方法可获得功能函数的概率密度函数。其中，最大熵原理具有统一形式的最大熵概率密度函数，能够较为精确地拟合单峰与多峰概率密度函数，被认为是拟合概率密度函数最具有无偏估计的一类方法[18]，因而学者们将该方法广泛应用于工程结构可靠度分析[18-22]。

为此，本文在结构整体可靠度问题分类的基础上给出其对应等价功能函数的统一描述，进而提出了将有效维度两步分析法与共轭无迹变换法结合的改进统计矩点估计法，并结合最大熵原理发展了一类兼顾精度与效率的结构整体可靠度分析方法。

1 结构整体可靠度问题描述

根据多失效模式产生的来源不同，可将结构整体可靠度问题分为两类：结构经典整体可靠度问题和结构一般整体可靠度问题[23]。结构经典整体可靠度问题主要研究理想弹塑性结构倒塌或形成可变机构的概率，即结构经典整体可靠度问题的失效模式由结构倒塌这一物理机制进行识别[7]；而结构一般整体可靠度问题的失效模式可以由失效准则直接确定。例如考察某多层结构(层数n≥2)的可靠度问题，任意某层的位移大于某限值时即认为结构失效，则该结构体系具有n个失效判别准则，相应地对应着n个失效模式[23]。因此，结构一般整体可靠度的多失效模式源于可靠度问题本身。

1.1 结构经典整体可靠度

当理想弹塑性结构承受完全相关的随机荷载作用时，随机向量Θ={ΘL,ΘS}T，其中ΘS和ΘL分别表示结构的随机参数和荷载的随机参数。定义ΘL,1为参考荷载，则ΘL可表示为：

其中，r(ΘL)为荷载向量。根据非线性发展过程[24]，弹塑性结构整体可靠度可等效为结构的极限承载力Fmax大于所施加荷载ΘL,1的概率。因此，结构经典整体可靠度对应的功能函数可等价描述为：

对于非完全相关荷载下的结构经典整体可靠度分析，详见文献[25]。

1.2 结构一般整体可靠度

当失效模式已知或已识别出时，此类结构整体可靠度问题称为结构一般整体可靠度。对于结构一般整体可靠度问题，其涉及的各失效模式及其逻辑关系均可根据预定的失效准则直接确定。

根据等价极值事件[8]，串联系统、并联系统及混联系统这三类逻辑关系的结构一般整体可靠度问题的等价功能函数可表示为：

式中：Gi(·)为串联系统或并联系统中第i个单元的功能函数；Gij(·)为第i个串联子系统中第j个单元的功能函数。

综上，结构经典整体可靠度与结构一般整体可靠度对应的等价功能函数可统一表示为：

2 基于有效维度两步分析法和共轭无迹变换法的改进统计矩点估计法

2.1 变量的独立标准正态变换

对于随机向量Θ={Θ1, Θ2,···,ΘN}，通过引入Nataf 变换[26]，可将其转换到独立标准正态空间，得到新的随机向量U={U1,U2,···,UN}。于是，两类结构整体可靠度问题对应的功能函数可改写为：

式中，N-1(·)表示Nataf 变换的逆变换。

2.2 传统的三变量降维近似统计矩点估计法

根据文献[11]，Z的原点矩MZ,k可近似计算如下：

虽然式(7)将高维积分转换为多个低维积分使得结构分析次数得到缩减，但三维和二维积分所需的结构分析次数仍较多。据此，本文提出有效维度两步分析法，并结合共轭无迹变换法进一步提高式(7)的计算效率。

2.3 有效维度两步分析法

其中，ɛ为衡量影响程度的阈值，本文中取ɛ=5%[15]。

采用有效维度分析法进行交互项评估仅需2N+2 次结构分析，具有较高的效率。然而，该方法存在交互项误判问题[21]。鉴于此，本文首先通过有效维度分析法进行交互项预判，在此基础上引入一次交叉项判别准则[12]进行二次判定。定义Us和Ut的交互项判定的示性函数Qh,2(Us,Ut)为[12]：

本文在有效维度分析法的基础上引入二次判定，故将此方法称为有效维度两步分析法。

2.4 改进统计矩点估计法

相比于传统的高斯求积公式，研究者近年来在非线性滤波领域发展的共轭无迹变换法具备更高的计算效率[27-29]。本文引入9 阶代数精度的共轭无迹变换法进行分量函数积分计算，其基本原理与过程如下。

其次，通过Isserlis 定理[27]，可建立如下矩约束方程：

式中：Uj为随机变量；为求积节点的坐标；ν为随机变量数目；κj为非负的整数且1≤κ1+κ2+···+κν≤9；n为除零点外求积节点的数量。求解式(19)，可获取缩放变量ri和权系数wi。特别的，当ν=2 或3 时，缩放变量ri和权系数wi如表1所示。

表1 二维及三维系统共轭无迹变换方法对应的ri 和wiTable 1 ri and wi in conjugate unscented transformation method for two-dimensional and three-dimensional systems

据 此，令ui={us,i,ut,i}或ui={us,i,ut,i,uv,i}，ui与αi分别为二维或三维共轭无迹变换方法的节点与权系数，则式(7)中的二维及三维积分可由下式确定：

其中，N1=4N(N-1)(N-2)/3。

对于式(7)中的一维积分，直接采用高斯求积公式进行计算：

式中：uq,i和ωi分别为Gauss-Hermite 积分的求积节点和权系数；本文取m=5。

对于式(7)中的二维积分，采用下式计算：

对于式(7)中三维积分，若Qh(Us,Ut,Uv)=1，则：

若Qh(Us,Ut,Uv)=0，则存在如下三种情况：

1) 当任意一对变量的交叉项不存在，如Qh,2(Us,Ut)=0 时，经推导，式(7)中三维积分可简化为：

2) 当任意两对变量的交叉项不存在，如Qh,2(Us,Ut)=0 且Qh,2(Us,Uv)=0 时，经推导，式(7)中三维积分可简化为：

3) 当三对变量的交叉项均不存在，经推导，式(7)中三维积分可简化为：

3 基于最大熵原理的结构整体可靠度分析

3.1 最大熵原理

结构响应统计矩获取后，通过最大熵原理进行结构整体可靠度评估[18-19]。以功能函数Z的熵取最大值为目标函数，以功能函数Z的前4 阶原点矩MZ,k(k=0,1,···,4)为约束条件，建立优化模型为：

对熵函数H引入Lagrange 乘子，构造Lagrange函数，得到：

通过求解非线性方程可获取Lagrange 系数λ0, λ1, ···, λ4，从而估计响应函数Z的密度函数：

则结构的失效概率与可靠指标为：

3.2 算法流程

基于本文方法的结构整体可靠度分析基本步骤如下：

1) 通过式(5)将结构整体可靠度问题的功能函数统一表示。

2) 通过式(6)将任意分布转换到标准独立正态空间。

3) 根据式(7)计算功能函数的前四阶原点矩：

a) 根据式(21)和式(22)分别计算单变量函数和双变量函数的积分值；

b) 根据Qh(Us,Ut,Uv)的结果通过式(23)～式(26)计算三变量函数的积分值；

c) 将单变量函数、双变量函数及三变量函数的积分值代入式(7)中，计算功能函数的前四阶矩。

4) 基于求得的功能函数前四阶矩，按照式(27) ～式(30)计算结构的失效概率和可靠指标。

4 算例分析

本节通过结构一般整体可靠度与结构经典整体可靠度两个算例验证本文方法的计算效率、精度以及适用性，分别将本文方法分别与Monte Carlo 模拟方法(MCS)、宽界限法(WBM)[30]、窄界限法(NBM)[31]、5 点Gauss-Hermite 积分的双变量降维近似统计矩点估计法(BDRM)[11]及5 点Gauss-Hermite 积分的三变量降维近似统计矩点估计法(TDRM)[11]计算的统计矩及可靠度进行了对比。本文将MCS 计算结果视为标准解，其他各方法的相对误差为：

式中：Value为对应方法的计算结果；MCV为Monte Carlo 方法计算的结果。

4.1 算例1. 结构一般整体可靠度问题

考察一个如图1 所示的六层钢筋混凝土框架结构，框架的跨度l=7.5 m，层高h=3 m，梁与柱的截面尺寸分别为300 mm×400 mm 和500 mm×500 mm，水平荷载F可表示为：

图1 六层钢筋混凝土框架结构Fig. 1 Six story reinforced concrete frame structure

考虑梁柱的弹性模量Eb和Ec，竖向荷载P1和P2以及水平荷载F0为随机变量，各随机变量的统计特性如表2 所示。考虑结构体系的最大层间位移角超过某限值为失效状态，该可靠度问题可以描述为：

表2 算例1 中随机变量的统计特征Table 2 Statistical characteristics of random variables for Example 1

式中：Xj(j=1,2,···,6)代表第j层与第j-1 层的相对位移； φB为层间位移角限值，本文中 φB取值为1/50。

该问题可归类为结构一般整体可靠度问题，根据式(6)，其单一等价功能函数可以表示为：

通过有效维度两步分析法易知，Qh(U1,U3,U4) =Qh(U2,U3,U4) =Qh(U3,U4,U5) =Qh,2(U3,U4) =0，进而，可计算Z的前四阶统计矩和可靠指标，其计算结果如表3 和表4 所示，表中ɛ表示该方法与Monte Carlo 方法的相对误差。

表3 算例1 中极限状态函数前四阶矩结果Table 3 Results of first four order moments of limit state function for Example 1

表4 算例1 中各类方法的可靠指标与有限元分析次数Table 4 Reliability index and number of finite element analysis of each method for Example 1

由表3 可以发现，TDRM 和本文方法估计的前四阶矩均具有较高的精度，其最大相对误差低于2%；虽然BDRM 亦能精确地估计前两阶矩，然而其估计的三阶矩具有6.32%的相对误差。就可靠指标的计算精度与效率而言，由表4 易知：采用经典的WBM 仅能给出较大的可靠指标界限值，其上、下界限可靠指标均具有较大的相对误差；虽然NBM 计算的可靠指标的界限值较窄，但其仍难以定量评估结构整体可靠度；通过BDRM估计的可靠指标相对误差达到12.5%，亦难以精确评估结构整体可靠度；幸运的是，本文方法和TDRM 能够精确地估计可靠指标，其最大误差低于1%，但本文方法所需的有限元分析次数为650 次，相较于TDRM 所需的1528 次有限元分析，本文方法计算效率提高2.35 倍左右。综合计算精度与效率，推荐采用本文方法进行结构一般整体可靠度的评估。

4.2 算例2. 结构经典整体可靠度问题

考察如图2 所示的理想弹塑性单层单跨钢框架结构的整体可靠度问题，其中，材料的弹性模量E为2.0×105MPa，各构件截面积均为0.09 m2。考虑节点荷载F0以及各构件截面抗弯承载力M1、M2和M3为相互独立的随机变量，其统计特性如表5 所示。当钢框架结构变成可变机构时认为结构失效。

图2 单层单跨的理想弹塑性钢框架结构 /mFig. 2 One-story and one-bay perfectly elastoplastic steel frame structure

表5 算例2 中随机变量的统计特征Table 5 Statistical characteristics of random variables for Example 2

根据问题描述可知，该问题为结构经典整体可靠度问题。根据式(6)，该问题的单一等价功能函数可表示为：

通过有效维度两步分析法可得，Qh(U1,U2,U3) =Qh(U1,U3,U4) =Qh,2(U1,U3) = 0，进而，采用各类方法计算出Z的前四阶统计矩和可靠指标如表6 和表7 所示。

表6 算例2 中极限状态函数前四阶矩结果Table 6 Results of first four order moments of limit state function for Example 2

表7 算例2 中各类方法的可靠指标与有限元分析次数Table 7 Reliability index and number of finite element analysis of each method for Example 2

通过表6 和表7 可以看出，采用BDRM 估计的三阶矩和可靠指标的相对误差均较大，其三阶矩的相对误差高达38.97%；采用经典的WBM 和NBM 仅能给出可靠指标的界限值难以精确量化，且需进行结构失效模式识别步骤；而采用TDRM及本文方法估计的统计矩和可靠指标均能够得到较为精确的结果，其相对误差均小于2%，但本文方法的计算效率较TDRM 提升2.5 倍多，效果显著。

5 结论

结构整体可靠度的定量评估是可靠度领域的基本问题之一。本文结合有效维度两步分析法和共轭无迹变换法，发展了高效的改进统计矩点估计法；进而，结合最大熵原理和结构整体可靠度等价描述法，提出了适用于结构经典和一般整体可靠度的分析方法。结果表明：

(1) 本文方法对两类结构整体可靠度问题均表现出较高的精度水平。

(2) 相比于传统的三变量降维近似统计矩点估计法，本文方法能在保证精度的前提下具有更高的效率；

(3) 双变量降维近似统计矩点估计法对于两类整体可靠度问题均具有较高的相对误差，难以对结构整体可靠度问题进行精确评估。