陈志强1,郑史雄1,2
(1.西南交通大学 土木工程学院,成都 610031; 2.陆地交通地质灾害防治技术国家工程实验室(西南交通大学),成都 610031)
摘要:
为了对复杂非线性桥梁结构进行随机地震作用下的动力可靠度分析,基于最大熵原理,构建了一种高效的桥梁非线性地震可靠度分析方法。首先阐述了桥梁结构地震可靠度分析与极值分布估计的关系,将桥梁结构的地震可靠度分析转化为对应的极值分布估计；然后建立了复杂桥梁结构非线性地震响应极值分布估计的最大熵方法,针对现有的最大熵分布求解过程受初值影响较大、不容易收敛等问题,提出一种基于似然函数的最大熵分布求解方法,从而实现了桥梁结构地震响应极值分布和地震可靠度的高效求解；最后以一座典型的高墩大跨连续刚构桥为例,通过蒙特卡洛模拟验证所提方法的精度和效率,并将结果与核密度估计和对数正态分布拟合的计算结果进行对比分析。结果表明:基于似然函数的最大熵分布求解方法不仅能够获得全局最优解,而且求解过程不受初值的影响、具有较好的数值稳定性,所提方法能够对随机地震作用下桥梁结构地震响应的极值分布和动力可靠度进行精确估计；核密度估计和对数正态分布无法对小失效概率水平下结构地震响应的极值分布进行估计,建议采用最大熵方法对桥梁结构进行地震可靠度分析。
关键词: 桥梁 地震可靠度 非线性 最大熵原理 极值分布 分数矩
DOI：10.11918/202005018
分类号:U442.5+5
文献标识码:A
基金项目:国家自然科学基金(U1434205)
Improved maximum entropy method for seismic reliability analysis of bridges under uncertainty
CHEN Zhiqiang1,ZHENG Shixiong1,2
(1.School of Civil Engineering, Southwest Jiaotong University, Chengdu 610031, China; 2. National Engineering Laboratory for Technology of Geological Disaster Prevention in Land Transportation(Southwest Jiaotong University), Chengdu 610031, China)
Abstract:
For the dynamic reliability analysis of complex nonlinear bridge structures under stochastic ground motions, an efficient seismic reliability analysis method was developed based on the principle of maximum entropy. First, the relationship between the seismic reliability and extreme value distribution (EVD) of bridges was clarified. The seismic reliability of bridge structures was transformed into the corresponding EVD estimation. Then, the maximum entropy method (MEM) for estimating the EVD of nonlinear seismic responses of complex bridge structures was established. In view of the fact that the iterative solution of the existing MEM is greatly affected by the initial value and is not easy to converge, a method based on likelihood function was proposed to solve the EVD and seismic reliability of bridge structures more efficiently. Finally, a typical high-pier and long-span continuous rigid frame bridge was taken as an example. The accuracy and efficiency of the proposed method were verified via Monte Carlo simulation, and the results were compared with the results of kernel density estimation (KDE) and lognormal distribution fitting. Results show that the proposed MEM based on likelihood function could obtain the global optimal solution of EVD, and the solution process was not affected by the initial value and had good numerical stability. The method could accurately estimate the EVD and dynamic reliability of complex bridge structures under stochastic earthquakes. Both KDE and lognormal distribution failed to estimate the EVD of the structural seismic responses at small failure probability levels. It is thereby recommended that MEM is adopted for the seismic reliability analysis of bridge structures.
Key words: bridge seismic reliability nonlinearity maximum entropy principle extreme value distribution fractional moments
陈志强, 郑史雄. 随机桥梁地震可靠度分析的改进最大熵方法[J]. 哈尔滨工业大学学报, 2021, 53(9): 79-87. DOI: 10.11918/202005018.
CHEN Zhiqiang, ZHENG Shixiong. Improved maximum entropy method for seismic reliability analysis of bridges under uncertainty[J]. Journal of Harbin Institute of Technology, 2021, 53(9): 79-87. DOI: 10.11918/202005018.
基金项目 国家自然科学基金(U1434205) 作者简介 陈志强(1993—)，男，博士研究生; 郑史雄(1965—)，男，教授，博士生导师 通信作者 郑史雄，zhengsx@home.swjtu.edu.cn 文章历史 收稿日期: 2020-05-05
Abstract Full text Figures/Tables PDF
随机桥梁地震可靠度分析的改进最大熵方法
陈志强1, 郑史雄1,2
1. 西南交通大学 土木工程学院，成都 610031;
2. 陆地交通地质灾害防治技术国家工程实验室(西南交通大学)，成都 610031
收稿日期: 2020-05-05
基金项目: 国家自然科学基金(U1434205)
作者简介: 陈志强(1993—)，男，博士研究生; 郑史雄(1965—)，男，教授，博士生导师
通信作者: 郑史雄，zhengsx@home.swjtu.edu.cn
摘要: 为了对复杂非线性桥梁结构进行随机地震作用下的动力可靠度分析，基于最大熵原理，构建了一种高效的桥梁非线性地震可靠度分析方法。首先阐述了桥梁结构地震可靠度分析与极值分布估计的关系，将桥梁结构的地震可靠度分析转化为对应的极值分布估计；然后建立了复杂桥梁结构非线性地震响应极值分布估计的最大熵方法，针对现有的最大熵分布求解过程受初值影响较大、不容易收敛等问题，提出一种基于似然函数的最大熵分布求解方法，从而实现了桥梁结构地震响应极值分布和地震可靠度的高效求解；最后以一座典型的高墩大跨连续刚构桥为例，通过蒙特卡洛模拟验证所提方法的精度和效率，并将结果与核密度估计和对数正态分布拟合的计算结果进行对比分析。结果表明: 基于似然函数的最大熵分布求解方法不仅能够获得全局最优解，而且求解过程不受初值的影响、具有较好的数值稳定性，所提方法能够对随机地震作用下桥梁结构地震响应的极值分布和动力可靠度进行精确估计；核密度估计和对数正态分布无法对小失效概率水平下结构地震响应的极值分布进行估计，建议采用最大熵方法对桥梁结构进行地震可靠度分析。
关键词: 桥梁 地震可靠度 非线性 最大熵原理 极值分布 分数矩
Improved maximum entropy method for seismic reliability analysis of bridges under uncertainty
CHEN Zhiqiang1, ZHENG Shixiong1,2
1. School of Civil Engineering, Southwest Jiaotong University, Chengdu 610031, China;
2. National Engineering Laboratory for Technology of Geological Disaster Prevention in Land Transportation(Southwest Jiaotong University), Chengdu 610031, China
Abstract: For the dynamic reliability analysis of complex nonlinear bridge structures under stochastic ground motions, an efficient seismic reliability analysis method was developed based on the principle of maximum entropy. First, the relationship between the seismic reliability and extreme value distribution (EVD) of bridges was clarified. The seismic reliability of bridge structures was transformed into the corresponding EVD estimation. Then, the maximum entropy method (MEM) for estimating the EVD of nonlinear seismic responses of complex bridge structures was established. In view of the fact that the iterative solution of the existing MEM is greatly affected by the initial value and is not easy to converge, a method based on likelihood function was proposed to solve the EVD and seismic reliability of bridge structures more efficiently. Finally, a typical high-pier and long-span continuous rigid frame bridge was taken as an example. The accuracy and efficiency of the proposed method were verified via Monte Carlo simulation, and the results were compared with the results of kernel density estimation (KDE) and lognormal distribution fitting. Results show that the proposed MEM based on likelihood function could obtain the global optimal solution of EVD, and the solution process was not affected by the initial value and had good numerical stability. The method could accurately estimate the EVD and dynamic reliability of complex bridge structures under stochastic earthquakes. Both KDE and lognormal distribution failed to estimate the EVD of the structural seismic responses at small failure probability levels. It is thereby recommended that MEM is adopted for the seismic reliability analysis of bridge structures.
Keywords: bridge seismic reliability nonlinearity maximum entropy principle extreme value distribution fractional moments
在基于性能的地震工程理论框架下，桥梁的地震可靠度分析对于桥梁的地震损伤评估和地震风险分析都具有重要意义[1]。一旦桥梁的地震可靠度确定，其地震风险可以通过失效概率对不同地震灾害水平进行积分计算[2]。因此，从结构工程的观点来看，桥梁地震风险分析中最为重要的部分就是桥梁的地震可靠度分析。
在实际工程中，地震动和桥梁结构本身通常都存在显著的不确定性。原则上来说，在桥梁结构地震可靠度分析中最为合理的应该同时考虑这两类不确定性，然而由于动力可靠度理论的不完善，现有的桥梁地震可靠度分析中通常采用随机振动理论和随机有限元理论两类方法分别对地震动的不确定性和桥梁结构参数的不确定性单独进行考虑[3]。当考虑地震动的不确定性时，桥梁的地震可靠度分析主要采用随机振动理论进行[3-5]，然而现有的随机振动理论主要是针对线弹性结构。虽然也有学者建立了针对非线性结构的随机振动分析方法，如路径积分法[6]、矩近似法[7]、等效线性化方法[8]以及概率密度演化理论[9-10]等，但是这些方法通常都只适用于简单结构，大规模、强非线性结构的随机地震响应分析仍然是一个巨大挑战。当考虑桥梁结构参数的不确定性时，桥梁的地震可靠度分析通常采用随机有限元方法[11-14]。然而，现有的方法不仅无法考虑地震动的不确定性，而且在复杂结构二阶统计量的计算上都还存在较大难度[3]。另外，对于非线性桥梁结构而言，其地震响应具有显著的非高斯性，仅采用一阶和二阶统计矩对其进行地震可靠度分析精度还远远不够。当同时考虑桥梁结构和地震动的不确定性时，桥梁地震可靠度分析方法的唯一可行方法只有基于抽样的数值模拟方法，如蒙特卡洛模拟(Monte Carlo simulation, MCS)[15]、子集模拟[16]以及重要性抽样[17]等。这些方法虽然理论上可以获得结构地震可靠度的精确结果，但是其所需要的样本数量通常都与失效概率成反比。在小失效概率水平下，采用这些抽样数值方法对桥梁结构进行地震可靠度分析时所需要的动力分析计算量非常大，同样很难适用于大跨度斜拉桥、悬索桥以及高墩连续刚构桥等复杂桥梁结构的地震可靠度分析。
最近几年，随着分数阶矩的出现，基于最大熵原理的结构动力可靠度分析获得了研究者极大的关注[14, 18-20]。文献[18]分别采用单变量降维模型、双变量降维模型以及稀疏网格积分对结构动力系统的分数矩进行估计，然后将其与最大熵模型相结合，建立了结构非线性地震可靠度分析的最大熵方法。文献[14]采用乘法形式的降维模型对结构地震响应的矩进行估计，对一座典型高墩桥梁进行了地震可靠度分析。文献[19]以最大熵原理为基础，提出了结构动力可靠度分析中的非等权重拟蒙特卡洛方法。文献[20]采用拉丁超立方抽样对结构地震响应分数矩进行估计，并将其与最大熵原理相结合，建立了近断层脉冲地震作用下高墩大跨连续刚构桥的地震可靠度分析方法。
由于分数阶矩中含有高阶中心矩信息，通过少数几阶分数矩就可以对结构地震响应极值分布进行准确模拟，因此基于分数矩的最大熵原理在桥梁地震可靠度分析中具有明显的优势。然而，由于现有的分数矩最大熵原理需要对多个变量进行迭代求解，在采用最大熵原理求解结构地震响应极值分布时求解过程受初值影响非常大，而且经常面临着优化过程不稳定的问题[18]。为了克服这一问题，本文提出了一种基于似然函数的最大熵分布模型，并将该模型应用于随机地震作用下，考虑非线性的随机桥梁结构地震可靠度分析，建立了基于最大熵原理的桥梁非线性地震可靠度分析方法。最后，以一座典型的高墩大跨连续刚构桥为例，通过蒙特卡洛模拟对所提出的方法进行验证。
1 桥梁地震可靠度问题阐述对于一个考虑结构参数和地震动双重不确定性的N自由度非线性动力系统，地震作用下，其动力方程可表示为
$\begin{aligned}&\boldsymbol{M}(\boldsymbol{\upsilon}) \ddot{\boldsymbol{U}}(t)+\boldsymbol{C}(\boldsymbol{\upsilon}) \dot{\boldsymbol{U}}(t)+\boldsymbol{G}[\boldsymbol{U}(t), \dot{\boldsymbol{U}}(t), \boldsymbol{\upsilon}]= \\&\quad-\boldsymbol{M}(\upsilon) \boldsymbol{I}\ddot{\boldsymbol{u}}_{g}(\boldsymbol{\chi}, t)\end{aligned}$ (1)
式中：M和C分别为质量矩阵和阻尼矩阵；$\ddot{\boldsymbol{U}}$(t)、$\mathit{\boldsymbol{\dot U}}$(t)和U(t)分别为结构加速度、速度和位移响应向量；υ=[υ1, υ2, …, υd1]T和χ=[χ1, χ2, …, χd2]T分别为桥梁结构和地震动随机参数向量；G为非线性恢复力向量，其通常为$\mathit{\boldsymbol{\dot U}}$(t)、U(t)和υ的函数；$\mathit{\boldsymbol{\ddot u}}$g(χ, t)为随机地面运动，通常可以采用谱表示法[21]或者Karhunen-Loeve (K-L)展开[22]进行模拟。
鉴于谱表示法具有简单、高效、易于实现等优点，因此本文采用谱表示法对地震动进行模拟。根据非平稳随机过程谱表示理论，一维地震动的表达式[21]通常为
$\ddot{u}_{g}(t)=\sum\limits_{l=1}^{m} \sqrt{2 S_{U}\left(\omega_{l}, t\right) \Delta \omega}\left[X_{l} \cos \left(\omega_{l} t\right)+Y_{l} \sin \left(\omega_{l} t\right)\right]$ (2)
式中：ωl=lΔω，Δω为频率区间；{Xl, Yl}为2m个标准正交随机变量；SU(ωl, t)为随机地震动的时-频演化功率谱，可以表示为
$S_{U}(\omega, t)=|A(\omega, t)|^{2} S(\omega)$ (3)
其中A(ω, t)为地震动的时-频调制函数。
根据谱表示理论，随机变量{Xl, Yl}(l=1, 2, …, m)应满足[21]:
$\begin{aligned}&E\left[X_{l}\right]=E\left[Y_{l}\right]=0, E\left[X_{l} Y_{r}\right]=0 \\&E\left[X_{l} X_{r}\right]=E\left[Y_{l} Y_{r}\right]=\delta_{r l}\end{aligned}$ (4)
其中: E[·]表示期望算子，δrl为Kronecker-Delta函数。
将模拟的地震动代入式(1)中的动力学方程，并采用数值方法对式(1)进行求解，由此可以得到桥梁结构地震可靠度分析所关心的物理量(如墩顶位移、支座相对位移等)，其通常可以表示为
$Z(t)=\boldsymbol{H}(\boldsymbol{\upsilon}, \boldsymbol{\chi}, \boldsymbol{d}, t)$ (5)
式中: d为确定性参数向量，H(·)为确定性算子。
对于通常所关心的首次超越地震可靠度，其可以定义为
$R=\operatorname{Pr}\left\{Z(t)<Z_{T}, \forall t \in[0, T]\right\}$ (6)
式中：Pr表示概率；T为地震动持时，ZT为桥梁结构抗震能力，其通常可以由规范规定或者根据Pushover分析进行确定。
由于式(6)中结构可靠度的求解非常困难，为了对其进行求解，这里通过等价极值原理[19]，将桥梁结构的地震可靠度转化为
$R=\operatorname{Pr}\left\{Z_{\mathrm{ext}}<Z_{T}, \forall t \in[0, T]\right\}$ (7)
式中Zext为桥梁结构非线性地震响应的极值，可以定义为
$Z_{\text {ext }}=\underset{t \in(0, T]}{\rm{Ext}}\{Z(t)\}=\varphi(\upsilon, \chi, T)$ (8)
其中$\mathop {{\mathop{\rm Ext}\nolimits} }\limits_{t \in (0, T]} $表示随机过程极值算子。若关心的是桥梁地震响应的绝对最大值，此时式(8)可以写为
$Z_{\mathrm{ext}}=\max \limits_{t \in(0, T]}\{|Z(t)|\}$ (9)
桥梁结构的地震可靠度则定义为
$R=\int_{0}^{Z_{T}} p_{Z_{\text {ext }}}(z) \mathrm{d} z$ (10)
式中pZext(z)为桥梁非线性地震响应极值Zext的概率密度，即桥梁结构的地震响应极值分布。同时，桥梁结构的失效概率则为
$P_{\mathrm{f}}=1-R=\int_{Z_{T}}^{\infty} p_{Z_{\mathrm{ext}}}(z) \mathrm{d} z$ (11)
综上所述，在桥梁地震可靠度分析中，最为重要的就是对其地震响应极值分布的估计。一旦桥梁地震响应的极值分布确定，其地震可靠度和失效概率就可以通过简单的数值积分进行求解。此外，对于工程中的大跨度复杂桥梁结构，都是高可靠性结构，在设计地震作用下，失效的概率通常都非常小。为了对其地震可靠度进行估计，需要对桥梁结构地震响应极值分布的尾部进行精确计算。鉴于此，本文提出一种针对复杂非线性结构的地震响应极值分布估计的高效数值方法。
2 结构地震可靠度分析的最大熵方法 2.1 结构地震响应极值分布的最大熵模型桥梁非线性地震响应的极值Zext在桥梁结构参数和地震动不确定性的影响下，通常是一个正的随机变量，为简单起见，这里将其表示为Z。根据最大熵原理，随机变量Z的无偏概率分布可以在分数矩约束下，通过熵的最大化进行求解。因此，根据最大熵原理桥梁非线性地震响应极值分布pZextz可以通过以下约束非线性优化问题求解得到[18-19]：
$\begin{aligned}&\max J(z)=-\int_{Z} p_{Z}(z) \ln p_{Z}(z) \mathrm{d} z \\&\text { s.t. }\left\{\begin{array}{l}\int_{Z} p_{Z}(z) \mathrm{d} z=1 \\\mu_{Z}{ }^{\alpha_{i}}=\int_{Z} z^{\alpha_{i}} p_{Z}(z) \mathrm{d} z, i=1,2, \cdots, M\end{array}\right.\end{aligned}$ (12)
式中μZαi为随机变量Z的第αi阶分数矩，J(z)为随机变量Z的信息熵，定义为
$J(z)=-\int_{Z} p_{Z} z) \ln p_{Z} z) \mathrm{d} z$ (13)
通过采用拉格朗日乘数法对式(12)进行求解，可以得到桥梁地震响应极值分布的一般形式[18]为
$\hat{p}_{Z}(z)=\exp \left(-\sum\limits_{i=0}^{M} \lambda_{i} z^{\alpha_{i}}\right)$ (14)
式中λi(i=0, 1, …, M)为M+1个拉格朗日乘子。
根据概率密度函数积分等于1的正则化条件可以得到
$\alpha_{0}=0, \lambda_{0}=\ln \left[\int_{Z} \exp \left(-\sum\limits_{i=1}^{M} \lambda_{i} z^{\alpha_{i}}\right) \mathrm{d} z\right]$ (15)
通过引入Kullback-Leibler(K-L)散度表征桥梁地震响应极值分布估计值$\hat{p}_{Z}(z)$与真实值pZ(z)之间的偏差，式(12)中的约束优化问题可以转变为以下无约束优化问题[19]：
$\min I(\boldsymbol{\alpha}, \boldsymbol{\lambda})=\ln \left[\int_{Z} \exp \left(-\sum\limits_{i=1}^{M} \lambda_{i} z^{\alpha_{i}}\right) \mathrm{d} z\right]+\sum\limits_{i=1}^{M} \lambda_{i} \mu_{Z}{ }^{\alpha_{i}}$ (16)
因此，pZext(z)的估计只需要对式(16)进行求解，得到分数指数向量α和拉格朗日乘子λ，并将其代入式(14)，即可得到桥梁非线性地震响应极值分布的无偏估计。
目前式(16)的求解普遍都是直接采用单纯形算法进行迭代求解，所得到的都是局部解，求解过程对α和λ的初值具有很强的依赖性，而且求解过程中经常存在数值不稳定的问题。为了克服这一问题，对式(16)中的优化问题进行快速求解，并保证求解过程的稳定性，本文提出一种新的基于似然函数的求解方法。
2.2 基于似然函数的最大熵求解方法根据分数矩的定义，对于随机变量Z，其第αi阶分数矩为
$\mu^{\alpha_{i}}=\int_{Z} g_{i}(z) p_{Z}(z) \mathrm{d} z$ (17)
式中gi(z)=zαi。将式(14)代入式(17)，并采用分部积分法可以得到
$\mu^{\alpha_{i}}=\sum\limits_{i=1}^{M} \lambda_{i} E\left[G_{i}(z) {g^{\prime}}_{j}(z)\right]=\sum\limits_{i=1}^{M} \lambda_{i} \omega_{i j}$ (18)
其中Gi(z)=∫Zgi(z)dz，而
$\omega_{i j}=E\left[\frac{\alpha_{j}}{\alpha_{i}+1} z^{\alpha_{i}+\alpha_{j}}\right]$ (19)
根据式(18)、(19)，拉格朗日乘子向量可以通过以下线性方程组进行确定[18]：
$\lambda=\boldsymbol{\omega}^{-1} \boldsymbol{\mu}_{Z}$ (20)
式中: ω=[ωij]，μZ=[μZα1, μZα2, …, μZαM]。
由于式(20)建立的初衷是求解整数矩约束下的最大熵方法的拉格朗日乘子，即当α1, α2, …, αM已知时，拉格朗日乘子λ可以通过式(20)中的线性方程组进行确定。但是在分数矩最大熵原理中，分数指数向量α是未知的。因此，根据式(16)和式(20)，随机变量Z概率分布的求解即转化为分数指数向量的确定。为了求解分数指数向量α，如果采用最大似然估计的思想，则可以通过定义似然函数L(M, α)[23]：
$L(M, \boldsymbol{\alpha})=\frac{1}{n_{s}} \sum\limits_{j=1}^{n_{s}} \ln \hat{p}_{Z}\left(z_{j}\right)-\frac{M}{n_{s}}$ (21)
根据最大似然原理，则可以使L(M, α)达到最大的分数指数向量$\mathit{\boldsymbol{\tilde \alpha }}$作为α的无偏估计。得到$\mathit{\boldsymbol{\tilde \alpha }}$之后，将其代入式(14)即可得到桥梁结构地震响应极值分布pZz的无偏估计：
$p_{Z}(z, \boldsymbol{\alpha})=\underset{\alpha}{\arg }\{\max \{L(M, \boldsymbol{\alpha})\}\}$ (22)
相比于式(16)，式(22)只需要对M个参数α1, α2, …, αM进行确定，因此其计算更为高效。虽然式(22)中的优化问题可以采用单纯形迭代、遗传算法等迭代方法进行求解，但是这些方法往往都容易陷入局部极值。为了对分数指数向量进行求解，这里通过对α进行离散，再进行随机组合，从而根据式(22)选出α最优解，最终得到结构地震响应极值分布的最佳无偏估计。求解结构地震响应极值分布的主要步骤如下：1)给定分数阶矩阶数M和分数指数α的潜在区间[αmin, αmax]以及增量Δα；2)定义离散的分数指数{α1, α2, …, αQ}，其中Q≥M，且αQ=αmax；3)从分数指数集{α1, α2, …, αQ}中无放回地随机抽取M个元素生成分数指数向量α，然后将得到的所有可能的α作为一个集合，构成分数指数向量集ψM, ψM中分数指数向量的个数NQ, M=$\left({\begin{array}{*{20}{l}} Q\\ M \end{array}} \right) = \frac{{Q!}}{{M!(Q - M)!}}$；4)从集合ψM中选择第r个元素α(r)，并通过式(18)~(20)求解拉格朗日乘子λ；5)将所求得的α(r)和λ代入式(14)从而得到结构非线性地震响应极值分布$\hat{p}_{Z}(z)$的初步值，并将结果代入式(21)，确定出其对应的似然函数；6)重复步骤4、5，确定出分数指数向量集中的每一个元素对应的似然函数值，选择其中使得似然函数L(M, α)最大的α和λ作为分数指数向量和拉格朗日乘子的最优解；7)将所得的分数指数向量α和拉格朗日乘子λ的最优解代入式(14)，即可得到结构非线性地震响应极值分布的无偏估计。
可以看到，通过步骤1~7，式(16)中分数矩最大熵模型转化为了NQ, M个线性方程组的求解，因此其计算效率大大提高。另外，由于分数指数是在其整个潜在区间[αmin, αmax]离散，因此该方法得到的是全局最优解。
2.3 桥梁地震响应极值分数矩估计由于式(2)是地震动$\ddot u$g(t)的近似表达式，为了使得频率截断误差足够小，m通常需要取到500~1 000项，这不仅使得在地震动模拟时需要涉及到处理高维随机变量，而且在结构非线性地震响应极值分布的分数矩估计时需要在一个超高维的概率空间中进行。为了减小概率空间的维度，进而减小结构动力分析的计算量，这里采用Liu等[21]近年来提出的基于正交函数的方法对随机变量{Xl, Yl}进行模拟。标准随机变量样本的生成如下：
$X_{l}=\varGamma\left(\operatorname{cas}\left(l \chi_{1}\right)\right), Y_{l}=\varGamma\left(\operatorname{cas}\left(l \chi_{2}\right)\right)$ (23)
式中：cas=sin(·)+cos(·)为Hartley正交基，Γ(·)为随机映射算子。
根据式(23)，地震动模拟中的高维随机变量就被减少到只有两个基本随机变量，因此地震动的不确定性将完全由随机变量χ1和χ1所表征。此时，概率空间维度为d=d1+2在结构地震响应极值分布分数矩的估计时，只需要在低维概率空间中进行即可。
对于结构非线性地震响应极值分布的第α阶分数矩，其通常可以估计如下：
$\mu^{\alpha}=\int_{X} z(\boldsymbol{X})^{\alpha} p_{X}(\boldsymbol{X}) \mathrm{d} \boldsymbol{X}=\sum\limits_{i=1}^{N} \widetilde{\omega}_{i} z\left(\overline{\boldsymbol{X}}_{i}\right)^{\alpha}$ (24)
式中：pX(X)为随机向量X的联合概率密度函数；Xi为积分点, i=1, 2, …, N; $\tilde \omega $i为积分权重，对于等权重抽样，$\tilde \omega $i=1/N。
为了确定估计中的积分点，文献[24]通过中心化的L2偏差对样本点集进行筛选，提出了一种改进的伪相关性折减拉丁超立方抽样方法(improved spurious correlation reduced LHS, ICLHS)。该方法在非线性结构动力响应的分数矩估计中取得了较好的效果，本文采用该方法对结构非线性地震响应极值分布的分数矩进行估计。
2.4 桥梁地震可靠度分析采用2.3节中的ICLHS方法生成桥梁结构和地震随机参数样本点，并通过式(2)对地震动进行模拟，从而形成结构-地震动样本对。然后采用非线性时程分析对结构地震响应进行求解，从而得到结构地震响应极值的分数矩。最后采用2.2节中所提出的最大熵方法对结构地震响应极值分布进行求解，即可得到桥梁结构非线性地震响应极值分布。再根据桥梁结构不同损伤极限状态的定义，通过对极值分布进行数值积分即可得到桥梁结构的地震可靠度。桥梁结构地震可靠度分析流程如图 1所示。
Fig. 1
图 1 桥梁地震可靠度分析流程图 Fig. 1 Flow chart of bridge seismic reliability analysis 3 算例分析 3.1 桥梁概况及非线性有限元模型以一座已经建成的高墩连续刚构桥为例，说明本文方法在复杂非线桥梁结构地震可靠度分析中的应用。该桥梁跨径布置为90 m+170 m+90 m，桥梁总长350 m，左右两个桥墩墩高分别为110.5 m和126.49 m，为典型的山区不规则高墩桥梁。桥梁上部主梁采用单箱单室预应力混凝土箱型梁，材料为C50预应力混凝土，主跨的跨中和桥墩的墩顶处截面梁高分别为3.7 m和10 m。为了适应主梁的变形，桥墩采用空心薄壁墩，墩身材料采用C40混凝土，桥梁结构详细布置如图 2所示。
Fig. 2
图 2 桥梁结构有限元模型示意(m) Fig. 2 Schematic diagram of finite element model of the bridge(m) 采用OpenSees作为分析平台，建立桥梁非线性有限元分析模型。主梁本文采用基于位移的梁柱单元并结合弹性截面进行模拟，不考虑主梁的非线性行为。桥墩作为地震中的易损构件，强震作用下，墩顶和墩底有可能进入屈服，产生塑性铰，本文采用基于力的非线性梁柱单元进行模拟，并通过纤维截面模拟桥墩的非线性行为。对于主梁和桥台伸缩缝位置处，强震作用下其可能发生严重碰撞，本文采用接触单元对梁-台之间的碰撞行为进行模拟，碰撞单元材料采用Hertz-damp模型[25]，从而充分考虑碰撞过程中的刚度变化和能量耗散。在桥台位置处，分别放置了两个盆式橡胶支座对结构变形进行约束，支座是地震中的易损构件，其通常都会产生很强的非线性，这里采用零长度单元结合硬化材料对其非线性行为进行模拟。桥梁结构有限元模型如图 2所示。
3.2 桥梁结构不确定性桥梁结构的不确定性通常主要包括结构材料、几何尺寸、边界条件以及阻尼比等参数。既有研究表明，混凝土抗压强度、钢筋屈服强度以及桥梁结构阻尼比等参数对桥梁结构的地震需求具有显著影响。本文在参考既有研究的基础上[26-27]，选出7个影响桥梁结构地震既有最为显著的参数作为考虑桥梁结构不确定性的随机变量，其概率分布见表 1。
表 1
表 1 桥梁结构的随机参数概率分布 Tab. 1 Probability distribution of random variables for bridge structures 模型参数 分布类型 参数均值 变异系数 桥墩混凝土弹性模量 正态 32.5 GPa 0.08
桥墩混凝土抗压强度 正态 32 MPa 0.135
钢筋弹性模量 正态 200 GPa 0.033
钢筋屈服强度 对数正态 335 MPa 0.05
主梁混凝土弹性模量 正态 35 500 MPa 0.08
支座摩擦系数 正态 0.06 0.5
阻尼比 正态 0.05 0.2
表 1 桥梁结构的随机参数概率分布 Tab. 1 Probability distribution of random variables for bridge structures
3.3 随机地震动功率谱模型为了考虑地震动的时-频非平稳性，这里采用文献[28]中演化功率谱(evolutionary power spectral density, EPSD)模型生成地震记录，从而对结构进行地震可靠度分析，其演化功率谱密度函数表示为
$\begin{gathered}S(\omega, t)=A(t)^{2} \frac{\omega_{e}(t)^{4}+4 \omega_{e}(t)^{2} \xi_{e}(t)^{2} \omega^{2}}{\left[\omega^{2}-\omega_{e}(t)^{2}\right]^{2}+4 \omega_{e}(t)^{2} \xi_{e}(t)^{2} \omega^{2}} \times \\\frac{\omega^{4}}{\left[\omega^{2}-\omega_{f}(t)^{2}\right]^{2}+4 \omega_{f}(t)^{2} \xi_{f}(t)^{2} \omega^{2}} \times S_{0}(t)\end{gathered}$ (25)
式中：A(t)为地震动非平稳调制函数，本文将其取为
$A(t)=\left[\frac{t}{c} \exp \left(1-\frac{t}{c}\right)\right]^{d}$ (26)
其中c控制地震动峰值到达时间，d控制地震动的形状，取值由桥位处所在场地确定。
此外，地震动演化功率谱模型中S0(t)为地震动谱强度因子，将其取为
$S_{0}(t)=\frac{\bar{a}_{\max }^{2}}{\bar{\gamma}^{2}\left[{\rm{ \mathit{ π} }} \omega_{e}(t)\left(2 \xi_{e}(t)+\frac{1}{2 \xi_{e}(t)}\right)\right]}$ (27)
$\omega_{e}(t)=\omega_{0}-a \frac{t}{T}, \xi_{e}(t)=\xi_{0}+b \frac{t}{T}$ (28)
式中：amax表示地震动的峰值加速度；ω0、ξ0、a和b分别为反映场地条件的地震动参数，这里分别取a=3、b=0.35、c=6、d=2、γ=2.75。
该算例桥梁抗震设防烈度为VIII度，故amax取为196 cm/s2，地震动演化功率谱如图 3所示。采用ICLHS方法生成400个随机样本，并根据式(2)即可生成模拟的地震记录样本。图 4给出了地震动样本平均功率谱(power spectral density, PSD)与目标值的比较，可以看到二者具有较好的一致性。另外，图 5给出了两条典型地震样本，模拟地震动均值和均方值与目标值的比较如图 6所示。从图 5中可以看到，模拟的地震动记录不仅具有显著的非平稳性，而且其峰值都在0.2g附近，与目标值较为吻合。从图 6中可见，模拟地震动均值、均方值与目标值也都十分吻合，由此说明了地震动降维模拟的有效性。
Fig. 3
图 3 地震演化功率谱 Fig. 3 EPSD of ground motions Fig. 4
图 4 模拟地震动功率谱与目标值比较 Fig. 4 Comparison between simulated PSD and target values Fig. 5
图 5 典型地震记录样本 Fig. 5 Typical samples of ground motions Fig. 6
图 6 模拟地震记录地震动均值和均方值与目标值比较 Fig. 6 Comparison of simulated mean and mean square values of ground motions with target values 3.4 桥梁结构地震可靠度分析方法验证根据桥梁结构随机参数概率分布和地震动随机参数，采用2.3节中的ICLHS方法生成400个积分点，通过非线性时程分析对结构地震响应进行求解即可得到400个桥梁结构地震响应极值。图 7以其中一个样本为例，给出了地震作用下1号墩的墩底弯矩-曲率滞回曲线。从图 7中可以看到，在地震作用下桥梁结构产生了显著的非线性变形。
Fig. 7
图 7 桥梁结构非线性地震响应典型样本 Fig. 7 Typical samples of structural nonlinear seismic responses 为了对本文方法进行验证，说明其在复杂结构非线性地震可靠度分析中的适用性。这里以1号墩的墩顶位移为例，图 8给出了1号墩的墩顶位移响应极值分数矩与MCS(104次)计算结果的对比。
Fig. 8
图 8 桥梁非线性地震响应分数矩 Fig. 8 Fractional moments of bridge nonlinear seismic responses 从图 8中可以看到，采用ICLHS方法估计的分数矩与MCS的结果十分吻合，特别是在-2到1阶内，二者相对误差不到2%，由此说明了ICLHS方法在复杂非线性地震响应的分数阶矩时具有较高的精度和效率。另外，从图 8中还可以看到当分数阶数较小时，桥梁结构地震响应极值分布分数矩的相对误差较小，而当分数阶数较高时，分数矩相对误差较大，由此说明了采用低阶分数矩能够更为准确地模拟结构地震响应极值分布。
若以桥梁结构非线性地震响应极值分数矩作为约束，采用2.2节中的方法对最大熵模型进行求解即可得到结构非线性地震响应极值分布。图 9给出了采用本文方法获得的1号墩墩顶位移非线性地震响应极值分布的概率密度函数(probability density function, PDF)、累积密度分布(cumulative density function, CDF)以及超越概率(probability of exceedance, POE)与MCS结果的对比，同时，图 9中还给出了采用核密度估计(kernel density estimation, KDE)以及对数正态分布拟合得到的结果。
Fig. 9
图 9 1号墩的墩顶位移地震响应极值分布 Fig. 9 EVD of seismic responses of top displacement of pier 1 从图 9中可以看到，采用本文方法获得的桥墩非线性地震响应极值分布与MCS结果具有较好的一致性，即使在小失效概率水平下(Pf < 10-3)，二者也十分吻合。这说明了采用该方法能够通过较少的几百次非线性动力分析，对大规模复杂非线性结构地震响应极值分布进行精确估计。同时，从图 9中还可以看到，不论极值分布的主体还是尾部分布，采用核密度估计和对数正态分布拟合的结果与蒙特卡洛模拟结果相比都存在较大误差，特别是对于尾部分布，二者间的误差非常大，由此说明了采用核密度估计和对数正态分布拟合都无法对小失效概率水平下结构的地震可靠度进行有效估计。
通常，对于非参数概率模型，其鲁棒性通常较差，如核密度估计，其结果都严重依赖于带宽的选择。为了说明2.2节中基于似然函数的最大熵模型求解方法的鲁棒性，图 10给出了采用不同分数指数增量Δα获得的桥梁非线性地震响应极值分布的对比。从图 10中可以看到，当分数指数增量Δα=0.1时，该方法获得的极值分布与蒙特卡洛模拟结果吻合得最好，二者基本完全重合。当Δα=0.2、0.25、0.3时，该方法获得的结果与蒙特卡洛模拟结果仍然一致，这说明了本文提出的基于似然函数的最大熵模型求解方法具有较好的鲁棒性。值得注意的是当Δα=0.2、0.25、0.3时，在尾部(10-4数量级)，本文方法获得的超越概率曲线与MCS结果出现了非常细微的偏差。这主要是由于本文所采用的算例为实际的高墩连续刚构桥，模型较为复杂，动力分析计算量较大，因此在蒙特卡洛模拟时只进行了104次抽样。在10-4数量级，蒙特卡洛模拟的结果自身就存在一定的误差，因此由于Δα的不同带来的这一点偏差基本是可以忽略不计的。同时，这也说明只要分数指数增量Δα的取值小于0.3，2.2节中基于似然函数的方法就能够对最大熵模型进行精确求解。
Fig. 10
图 10 分数指数增量Δα对极值分布的影响 Fig. 10 Effect of fractional exponential increment Δα on EVD 4 结论1) 建立了基于分数阶矩最大熵原理的桥梁非线性地震可靠度分析方法，通过与蒙特卡洛模拟结果的对比说明了本文方法能够通过少量的非线性动力分析实现小失效概率水平下结构非线性的地震响应极值分布的精确求解，该方法能够为大规模复杂结构非线性地震可靠度分析提供一种有效途径。
2) 提出了一种基于似然函数的最大熵模型求解方法，能够将最大熵模型求解过程中的无约束极值问题转化为线性方程组的求解，从而克服了采用单纯形方法求解最大熵分布迭代过程不稳定的问题。
3) 核密度估计和对数正态分布对桥梁非线性地震响应极值分布的主体具有一定近似能力，但是二者都无法准确模拟结构极值分布的尾部分布，而最大熵分布不仅能够模拟结构非线性地震响应极值分布的主体，而且对于尾部分布也具有非常好的近似能力。
参考文献
[1] GHOBARAH A. Performance-based design in earthquake engineering: state of development[J]. Engineering Structures, 2001, 23(8): 878. DOI:10.1016/S0141-0296(01)00036-0
[2] GHOSH S, ROY A, CHAKRABORTY S. Support vector regression based metamodeling for seismic reliability analysis of structures[J]. Applied Mathematical Modelling, 2018, 64: 584. DOI:10.1016/j.apm.2018.07.054
[3] XU J, FENG D C. Seismic response analysis of nonlinear structures with uncertain parameters under stochastic ground motions[J]. Soil Dynamics and Earthquake Engineering, 2018, 111: 149. DOI:10.1016/j.soildyn.2018.04.023
[4] 欧进萍, 王光远. 结构随机振动[M]. 北京: 高等教育出版社, 1998.
OU Jinping, WANG Guangyuan. Random vibration of structures[M]. Beijing: Higher Education Press, 1998.
[5] 林家浩, 张亚辉. 随机振动的虛拟激励法[M]. 北京: 科学出版社, 2004.
LIN Jiahao, ZHANG Yahui. Pseudo excitation method of random vibration[M]. Beijing: Science Press, 2004.
[6] DI-MATTEO A, DI-PAOLA M, PIRROTTA A. Path integral solution for nonlinear systems under parametric Poissonian white noise input[J]. Probabilistic Engineering Mechanics, 2016, 44: 89. DOI:10.1016/j.probengmech.2015.09.020
[7] CRANDALL S H. Non-Gaussian closure for random vibration of non-linear oscillators[J]. International Journal of Non-Linear Mechanics, 1980, 15(4/5): 303.
[8] ATALIK T S, UTKU S. Stochastic linearization of multi-degree-of-freedom non-linear systems[J]. Earthquake Engineering & Structural Dynamics, 1976, 4(4): 411.
[9] LI J, CHEN J B. Probability density evolution method for dynamic response analysis of structures with uncertain parameters[J]. Computational Mechanics, 2004, 34(5): 400. DOI:10.1007/s00466-004-0583-8
[10] 陈建兵, 李杰. 结构随机地震反应与可靠度的概率密度演化分析研究进展[J]. 工程力学, 2014, 31(4): 1.
CHEN Jianbing, LI Jie. probability density evolution method for stochastic seismic response and reliability of structures[J]. Engineering Mechanics, 2014, 31(4): 1.
[11] KLEIBER M, HIEN T D. The stochastic finite element method: basic perturbation technique and computer implementation[M]. New York: John Wiley & Sons, 1992.
[12] NAESS A, MOE V. Efficient path integration methods for nonlinear dynamic systems[J]. Probabilistic Engineering Mechanics, 2000, 15(2): 221. DOI:10.1016/S0266-8920(99)00031-4
[13] GHANEM R, SPANOS P D. A stochastic Galerkin expansion for nonlinear random vibration analysis[J]. Probabilistic Engineering Mechanics, 1993, 8(3/4): 255.
[14] ZHANG J, BI K, ZHENG S, et al. Seismic system reliability analysis of bridges using the multiplicative dimensional reduction method[J]. Structure and Infrastructure Engineering, 2018, 14(11): 1455. DOI:10.1080/15732479.2018.1450428
[15] SHINOZUKA M. Monte Carlo solution of structural dynamics[J]. Computers & Structures, 1972, 2(5/6): 855.
[16] AU S K, BECK J L. Estimation of small failure probabilities in high dimensions by subset simulation[J]. Probabilistic Engineering Mechanics, 2001, 16(4): 263. DOI:10.1016/S0266-8920(01)00019-4
[17] MELCHERS R E. Importance sampling in structural systems[J]. Structural Safety, 1989, 6(1): 3. DOI:10.1016/0167-4730(89)90003-9
[18] XU J. A new method for reliability assessment of structural dynamic systems with random parameters[J]. Structural Safety, 2016, 60: 130. DOI:10.1016/j.strusafe.2016.02.005
[19] XU J, ZHANG W, SUN R. Efficient reliability assessment of structural dynamic systems with unequal weighted quasi-Monte Carlo simulation[J]. Computers & Structures, 2016, 175: 37.
[20] CHEN Z Q, ZHENG S X, ZHOU Q, et al. Extreme value distribution and dynamic reliability estimation of high-pier bridges subjected to near-fault impulsive ground motions[J]. Advances in Structural Engineering, 2020, 23(7): 1367. DOI:10.1177/1369433219894245
[21] LIU Z, LIU W, PENG Y. Random function based spectral representation of stationary and non-stationary stochastic processes[J]. Probabilistic Engineering Mechanics, 2016, 45: 115. DOI:10.1016/j.probengmech.2016.04.004
[22] PHOON K K, HUANG H W, QUEK S T. Simulation of strongly non-Gaussian processes using Karhunen-Loeve expansion[J]. Probabilistic Engineering Mechanics, 2005, 20(2): 188. DOI:10.1016/j.probengmech.2005.05.007
[23] ALIBRANDI U, MOSALAM K M. Kernel density maximum entropy method with generalized moments for evaluating probability distributions, including tails, from a small sample of data[J]. International Journal for Numerical Methods in Engineering, 2018, 113(13): 1904. DOI:10.1002/nme.5725
[24] 陈志强, 郑史雄, 张宁, 等. 非平稳地震作用下强非线性结构动力可靠度分析[J]. 振动与冲击, 2020, 39(19): 121.
CHEN Zhiqiang, ZHENG Shixiong, ZHANG Ning, et al. Dynamic reliability assessment of nonlinear structures under non-stationary ground motions[J]. Journal of Vibration and Shock, 2020, 39(19): 121.
[25] MUTHUKUMAR S. A contact element approach with hysteresis damping for the analysis and design of pounding in bridges[D]. Atlanta: Georgia Institute of Technology, 2003
[26] 成虎. 氯离子侵蚀钢筋混凝土桥墩地震易损性研究[D]. 大连: 大连理工大学, 2019
CHENG Hu. Seismic fragility analyses of reinforced concrete bridge piers under chloride-induced corrosion[D]. Dalian: Dalian University of Technology, 2019
[27] 宋帅, 钱永久, 钱聪. 桥梁地震需求中随机参数的重要性分析方法研究[J]. 工程力学, 2018, 35(3): 106.
SONG Shuai, QIAN Yongjiu, QIAN Cong. Research on methods for importance analysis of random parameter in bridge seismic demand[J]. Engineering Mechanics, 2018, 35(3): 106.
[28] CACCIOLA P, DEODATIS G. A method for generating fully non-stationary and spectrum-compatible ground motion vector processes[J]. Soil Dynamics and Earthquake Engineering, 2011, 31(3): 351. DOI:10.1016/j.soildyn.2010.09.003
