在基于性能的地震工程理论框架下，桥梁的地震可靠度分析对于桥梁的地震损伤评估和地震风险分析都具有重要意义^[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)分别为结构加速度、速度和位移响应向量；υ=[υ₁, υ₂, …, υ_d1]^T和χ=[χ₁, χ₂, …, χ_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Δω，Δω为频率区间；{X_l, Y_l}为2m个标准正交随机变量；S_U(ω_l, t)为随机地震动的时-频演化功率谱，可以表示为

$ S_{U}(\omega, t)=|A(\omega, t)|^{2} S(\omega) $

(3)

其中A(ω, t)为地震动的时-频调制函数。

根据谱表示理论，随机变量{X_l, Y_l}(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为地震动持时，Z_T为桥梁结构抗震能力，其通常可以由规范规定或者根据Pushover分析进行确定。

由于式(6)中结构可靠度的求解非常困难，为了对其进行求解，这里通过等价极值原理^[19]，将桥梁结构的地震可靠度转化为

$ R=\operatorname{Pr}\left\{Z_{\mathrm{ext}}<Z_{T}, \forall t \in[0, T]\right\} $

(7)

式中Z_ext为桥梁结构非线性地震响应的极值，可以定义为

$ 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)

式中p_Zext(z)为桥梁非线性地震响应极值Z_ext的概率密度，即桥梁结构的地震响应极值分布。同时，桥梁结构的失效概率则为

$ P_{\mathrm{f}}=1-R=\int_{Z_{T}}^{\infty} p_{Z_{\mathrm{ext}}}(z) \mathrm{d} z $

(11)

综上所述，在桥梁地震可靠度分析中，最为重要的就是对其地震响应极值分布的估计。一旦桥梁地震响应的极值分布确定，其地震可靠度和失效概率就可以通过简单的数值积分进行求解。此外，对于工程中的大跨度复杂桥梁结构，都是高可靠性结构，在设计地震作用下，失效的概率通常都非常小。为了对其地震可靠度进行估计，需要对桥梁结构地震响应极值分布的尾部进行精确计算。鉴于此，本文提出一种针对复杂非线性结构的地震响应极值分布估计的高效数值方法。

2 结构地震可靠度分析的最大熵方法 2.1 结构地震响应极值分布的最大熵模型

桥梁非线性地震响应的极值Z_ext在桥梁结构参数和地震动不确定性的影响下，通常是一个正的随机变量，为简单起见，这里将其表示为Z。根据最大熵原理，随机变量Z的无偏概率分布可以在分数矩约束下，通过熵的最大化进行求解。因此，根据最大熵原理桥梁非线性地震响应极值分布p_Zextz可以通过以下约束非线性优化问题求解得到^[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)$与真实值p_Z(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)

因此，p_Zext(z)的估计只需要对式(16)进行求解，得到分数指数向量α和拉格朗日乘子λ，并将其代入式(14)，即可得到桥梁非线性地震响应极值分布的无偏估计。

目前式(16)的求解普遍都是直接采用单纯形算法进行迭代求解，所得到的都是局部解，求解过程对α和λ的初值具有很强的依赖性，而且求解过程中经常存在数值不稳定的问题。为了克服这一问题，对式(16)中的优化问题进行快速求解，并保证求解过程的稳定性，本文提出一种新的基于似然函数的求解方法。

2.2 基于似然函数的最大熵求解方法

根据分数矩的定义，对于随机变量Z，其第α_i阶分数矩为

$ \mu^{\alpha_{i}}=\int_{Z} g_{i}(z) p_{Z}(z) \mathrm{d} z $

(17)

式中g_i(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)

其中G_i(z)=∫_Zg_i(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^α₁, μ_Z^α₂, …, μ_Z^α_M]。

由于式(20)建立的初衷是求解整数矩约束下的最大熵方法的拉格朗日乘子，即当α₁, α₂, …, α_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)即可得到桥梁结构地震响应极值分布p_Zz的无偏估计：

$ p_{Z}(z, \boldsymbol{\alpha})=\underset{\alpha}{\arg }\{\max \{L(M, \boldsymbol{\alpha})\}\} $

(22)

相比于式(16)，式(22)只需要对M个参数α₁, α₂, …, α_M进行确定，因此其计算更为高效。虽然式(22)中的优化问题可以采用单纯形迭代、遗传算法等迭代方法进行求解，但是这些方法往往都容易陷入局部极值。为了对分数指数向量进行求解，这里通过对α进行离散，再进行随机组合，从而根据式(22)选出α最优解，最终得到结构地震响应极值分布的最佳无偏估计。求解结构地震响应极值分布的主要步骤如下：1)给定分数阶矩阶数M和分数指数α的潜在区间[α_min, α_max]以及增量Δα；2)定义离散的分数指数{α₁, α₂, …, α_Q}，其中Q≥M，且α_Q=α_max；3)从分数指数集{α₁, α₂, …, α_Q}中无放回地随机抽取M个元素生成分数指数向量α，然后将得到的所有可能的α作为一个集合，构成分数指数向量集ψ_M, ψ_M中分数指数向量的个数N_{Q, 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)中分数矩最大熵模型转化为了N_{Q, M}个线性方程组的求解，因此其计算效率大大提高。另外，由于分数指数是在其整个潜在区间[α_min, α_max]离散，因此该方法得到的是全局最优解。

2.3 桥梁地震响应极值分数矩估计

由于式(2)是地震动$\ddot u$_g(t)的近似表达式，为了使得频率截断误差足够小，m通常需要取到500~1 000项，这不仅使得在地震动模拟时需要涉及到处理高维随机变量，而且在结构非线性地震响应极值分布的分数矩估计时需要在一个超高维的概率空间中进行。为了减小概率空间的维度，进而减小结构动力分析的计算量，这里采用Liu等^[21]近年来提出的基于正交函数的方法对随机变量{X_l, Y_l}进行模拟。标准随机变量样本的生成如下：

$ 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)，地震动模拟中的高维随机变量就被减少到只有两个基本随机变量，因此地震动的不确定性将完全由随机变量χ₁和χ₁所表征。此时，概率空间维度为d=d₁+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)

式中：p_X(X)为随机向量X的联合概率密度函数；X_i为积分点, 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所示。

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所示。

采用OpenSees作为分析平台，建立桥梁非线性有限元分析模型。主梁本文采用基于位移的梁柱单元并结合弹性截面进行模拟，不考虑主梁的非线性行为。桥墩作为地震中的易损构件，强震作用下，墩顶和墩底有可能进入屈服，产生塑性铰，本文采用基于力的非线性梁柱单元进行模拟，并通过纤维截面模拟桥墩的非线性行为。对于主梁和桥台伸缩缝位置处，强震作用下其可能发生严重碰撞，本文采用接触单元对梁-台之间的碰撞行为进行模拟，碰撞单元材料采用Hertz-damp模型^[25]，从而充分考虑碰撞过程中的刚度变化和能量耗散。在桥台位置处，分别放置了两个盆式橡胶支座对结构变形进行约束，支座是地震中的易损构件，其通常都会产生很强的非线性，这里采用零长度单元结合硬化材料对其非线性行为进行模拟。桥梁结构有限元模型如图 2所示。

3.2 桥梁结构不确定性

桥梁结构的不确定性通常主要包括结构材料、几何尺寸、边界条件以及阻尼比等参数。既有研究表明，混凝土抗压强度、钢筋屈服强度以及桥梁结构阻尼比等参数对桥梁结构的地震需求具有显著影响。本文在参考既有研究的基础上^[26-27]，选出7个影响桥梁结构地震既有最为显著的参数作为考虑桥梁结构不确定性的随机变量，其概率分布见表 1。

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控制地震动的形状，取值由桥位处所在场地确定。

此外，地震动演化功率谱模型中S₀(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)

式中：a_max表示地震动的峰值加速度；ω₀、ξ₀、a和b分别为反映场地条件的地震动参数，这里分别取a=3、b=0.35、c=6、d=2、γ=2.75。

该算例桥梁抗震设防烈度为VIII度，故a_max取为196 cm/s²，地震动演化功率谱如图 3所示。采用ICLHS方法生成400个随机样本，并根据式(2)即可生成模拟的地震记录样本。图 4给出了地震动样本平均功率谱(power spectral density, PSD)与目标值的比较，可以看到二者具有较好的一致性。另外，图 5给出了两条典型地震样本，模拟地震动均值和均方值与目标值的比较如图 6所示。从图 5中可以看到，模拟的地震动记录不仅具有显著的非平稳性，而且其峰值都在0.2g附近，与目标值较为吻合。从图 6中可见，模拟地震动均值、均方值与目标值也都十分吻合，由此说明了地震动降维模拟的有效性。

3.4 桥梁结构地震可靠度分析方法验证

根据桥梁结构随机参数概率分布和地震动随机参数，采用2.3节中的ICLHS方法生成400个积分点，通过非线性时程分析对结构地震响应进行求解即可得到400个桥梁结构地震响应极值。图 7以其中一个样本为例，给出了地震作用下1号墩的墩底弯矩-曲率滞回曲线。从图 7中可以看到，在地震作用下桥梁结构产生了显著的非线性变形。

为了对本文方法进行验证，说明其在复杂结构非线性地震可靠度分析中的适用性。这里以1号墩的墩顶位移为例，图 8给出了1号墩的墩顶位移响应极值分数矩与MCS(10⁴次)计算结果的对比。

从图 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)以及对数正态分布拟合得到的结果。

从图 9中可以看到，采用本文方法获得的桥墩非线性地震响应极值分布与MCS结果具有较好的一致性，即使在小失效概率水平下(P_f < 10^-3)，二者也十分吻合。这说明了采用该方法能够通过较少的几百次非线性动力分析，对大规模复杂非线性结构地震响应极值分布进行精确估计。同时，从图 9中还可以看到，不论极值分布的主体还是尾部分布，采用核密度估计和对数正态分布拟合的结果与蒙特卡洛模拟结果相比都存在较大误差，特别是对于尾部分布，二者间的误差非常大，由此说明了采用核密度估计和对数正态分布拟合都无法对小失效概率水平下结构的地震可靠度进行有效估计。

通常，对于非参数概率模型，其鲁棒性通常较差，如核密度估计，其结果都严重依赖于带宽的选择。为了说明2.2节中基于似然函数的最大熵模型求解方法的鲁棒性，图 10给出了采用不同分数指数增量Δα获得的桥梁非线性地震响应极值分布的对比。从图 10中可以看到，当分数指数增量Δα=0.1时，该方法获得的极值分布与蒙特卡洛模拟结果吻合得最好，二者基本完全重合。当Δα=0.2、0.25、0.3时，该方法获得的结果与蒙特卡洛模拟结果仍然一致，这说明了本文提出的基于似然函数的最大熵模型求解方法具有较好的鲁棒性。值得注意的是当Δα=0.2、0.25、0.3时，在尾部(10^-4数量级)，本文方法获得的超越概率曲线与MCS结果出现了非常细微的偏差。这主要是由于本文所采用的算例为实际的高墩连续刚构桥，模型较为复杂，动力分析计算量较大，因此在蒙特卡洛模拟时只进行了10⁴次抽样。在10^-4数量级，蒙特卡洛模拟的结果自身就存在一定的误差，因此由于Δα的不同带来的这一点偏差基本是可以忽略不计的。同时，这也说明只要分数指数增量Δα的取值小于0.3，2.2节中基于似然函数的方法就能够对最大熵模型进行精确求解。

4 结论

1) 建立了基于分数阶矩最大熵原理的桥梁非线性地震可靠度分析方法，通过与蒙特卡洛模拟结果的对比说明了本文方法能够通过少量的非线性动力分析实现小失效概率水平下结构非线性的地震响应极值分布的精确求解，该方法能够为大规模复杂结构非线性地震可靠度分析提供一种有效途径。

2) 提出了一种基于似然函数的最大熵模型求解方法，能够将最大熵模型求解过程中的无约束极值问题转化为线性方程组的求解，从而克服了采用单纯形方法求解最大熵分布迭代过程不稳定的问题。

3) 核密度估计和对数正态分布对桥梁非线性地震响应极值分布的主体具有一定近似能力，但是二者都无法准确模拟结构极值分布的尾部分布，而最大熵分布不仅能够模拟结构非线性地震响应极值分布的主体，而且对于尾部分布也具有非常好的近似能力。