基于混合随机相依SFE-Copula模型的地铁施工诱发临近建筑的可靠性分析

董 超， 刘文黎， 王彦玉

(华中科技大学 土木工程与力学学院， 湖北 武汉 430074)

在地铁建设过程中，如何保证地铁本体及周围构筑物的安全稳定性尤为重要。在城市地铁建设过程中，新建地铁难免在既有建筑物的附近区域施工，导致这些临近建筑物出现变形、倾斜、扭曲等现象，甚至造成不可恢复的损坏[1,2]。因此，为确保地铁施工期间临近既有建筑物的安全性和可靠性，有必要对地铁开挖引起的建筑物扰动情况进行分析，并对临近建筑物的安全可靠性进行评估。

在地铁盾构施工对临近建筑物变形研究方面，最著名的为Peck公式，在一定的时期内，经验法是初步分析建筑变形的有效方法。Maleki认为数值分析方法在对地铁隧道施工对建筑物安全问题的分析上面更具优势[3]。三维有限元方法也被应用到地铁施工对临近建筑物桩基的影响中[4]。骆建军等对大跨度地铁盾构施工临近建筑物的安全保护问题提出了相应的监测和保护措施[5]。虽然国内外相关研究已经取得了许多有意义的研究成果，但地铁盾构施工建筑物的变形问题是一个涉及地质、环境等多参数相互作用的复杂问题，然而这些参数的不确定性在一定程度上限制了Peck公式、解析法等经典方法的使用[6]。

地下工程研究的各参数间存在一定的相依性，为确保可靠度分析的精度，需要构建指标之间的联合分布函数，从而需要大量数据以保证联合分布函数的精度，导致其在实际工程中实行困难。近些年来，在数学和金融领域应用较多的Copula函数广泛应用于多维参数联合概率分布建模中。Copula函数能够提供一种灵活的方式来表征多维变量数据之间的非线性相依性，现已广泛应用于金融和精算应用[7]、水文[8]和岩土工程[9, 10]等领域。

为精确描述地铁施工对临近建筑物的影响，本文采用有限元法构建地铁-建筑物相互作用模型，引用随机有限元法输入外界不确定性因素用以研究其对模型响应的影响，并选用四种不同的Copula函数构建临近建筑物的最大沉降Smax和最大倾斜率Imax二维联合分布模型，然后利用AIC(Akaike Information Criterion)和BIC(Bayesian Information Criterion)准则识别出最优的Copula模型，并使用最优Copula函数模型生成足够的数据用于蒙特卡罗模拟，从而得出地铁施工作用下临近建筑物的失效概率和可靠度。

1 SFE-Copula可靠性分析模型构建

本文分析了地铁施工对临近建筑物的安全可靠性影响，主要思路为：(1)通过ABAQUS软件，构建地铁和建筑的相互作用模型，精确描述地铁施工对临近建筑物的影响规律；(2)通过自行编制GUI命令，构建ABAQUS随机有限元模型，将周围环境的不确定性因素作为输入参数导入有限元模型中，从而模拟生成180个FE(Footloose Entrepreneurs)模拟数据；(3)采用4种不同的Copula函数构建Smax和Imax的二维Copula模型，并运用AIC和BIC准则确定最优Copula函数，以实现在构建模型中考虑参数之间的相依性；(4)通过确定的失效函数和蒙特卡洛法得出建筑物在周围环境不确定性因素影响下的失效概率，从而得到建筑物的可靠度指标，并通过建筑物的重要程度和可靠度指标β评估建筑物的安全风险等级。

1.1 使用ABAQUS建立随机有限元模型

传统有限元法主要关注确定性问题的求解，但是，岩土工程中存在多重不确定性因素[11]，本文选用随机有限元法来描述这些不确定性问题。

本文采用随机有限元法对隧道施工诱发周围建筑物安全风险的可靠度进行分析。基于ABAQUS软件的随机有限元模型构建过程如下所述。构建模型的核心思想为参数的数据更新，即编制一个Python程序[12]用于更新输入参数。其步骤如下：

(1)用Python程序编写ABAQUS图形用户界面(Graphical User Interface，GUI)；

(2)定义输入参数的分布形式，生成输入参数，在ABAQUS中执行确定性有限元分析程序，并将结果存储在输出的.odb文件中，重复此过程，产生与试验次数相同数量的.odb文件；

(3)编制另一个Python脚本从.odb文件中提取输入随机变量的值

1.2 Copula理论

Copula方法基于Sklar定理，是一种用于构建多变量联合分布模型的常用方法。Sklar定理为：设F(x1,x2,…,xn)为边缘分布函数的联合分布函数F(x1),F(x2),…,F(xn)。存在对于所有真实的n维Copula，x1,x2,…,xn[13]。

F(x1,x2,…,xn)=

C(F1(x1),F2(x2),…,Fn(xn))

(1)

式中：C为F(x1),F(x2),…,F(xn)的Copula函数。Sklar定理将多维参数联合分布的构建过程，分为两个步骤：(1)从样本数据估计边缘分布的统计参数；(2)最优Copula函数识别(本文选用AIC或BIC准则)。根据Sklar定理，两个随机变量x1和x2的两参数联合分布如下：

F(x1,x2) =C(F1(x1),F2(x2);θ)

=C(μ1,μ2;θ)

(2)

式中：C(μ1,μ2;θ)为Copula函数中的待定参数，θ为x1和x2之间的相关性。

1.2.1Copula函数选择

Copula函数可以用来表征两个随机变量之间的关联性。通常，采用不同Copula函数表征的相依性结构有很大不同，本文选择的Copula函数如表1所示[14]。选择这四种Copula函数的主要原因是，这四个Copula函数是常用的Copula，并且与几个典型的Copula族相关联。其中，Gaussian是椭圆族Copula函数，Clayton，Gumbel和Frank Copulas是常用的阿基米德Copula函数。

表1 本文采用的Copula函数类型

注：Φρ为正态分布的联合概率密度函数；Φ-1(·)为标准正态分布函数的反函数

1.2.2Copula函数相关参数的估计

Copula参数θ可以通过Pearson线性相关系数或秩相关系数(如Kendall相关系数)确定[13]。由于Pearson系数在线性条件下能保持不变，但在非线性条件下不稳定， Kendall秩相关τk用来度量两个随机变量的依赖性程度，可以克服这一问题[9]，故而本文采用τk确定Copula参数。对于具有N个值的两变量(x1,x2)的非参数估计τk可由式(3)计算：

(3)

这里，sign(·)定义为：

i,j=1,2,…,N

(4)

根据Copula理论，τk可以用Copula函数表示为：

(5)

因此，与给定的Copula函数相关联的Copula参数θ可以通过求解上述积分方程确定。

1.2.3 最优Copula函数识别

在确定Copula参数θ之后，可以唯一地得到Copula联合分布函数和Copula联合概率密度分布函数(见表1)。然后，可以使用Akaike信息准则(AIC)[13]和贝叶斯信息准则(BIC)[9]的计算值区分识别最优Copula函数，公式如下：

(6)

(7)

(8)

式中：k为Copula参数的数量，本文中k值为1；N为相关数据的样本大小。当AIC和BIC取值最小时所对应的Copula函数为最优Copula函数。

1.2.4 二维Copula函数的模拟方法

Copula函数不仅在构造变量的联合概率分布函数方面具备便利性及灵活性，并且还能模拟服从构造概率分布函数的相关非正态变量。利用Copula函数进行参数模拟的步骤如下：(1)模拟服从给定Copula函数的相关标准均分布变量U=(U1,U2)T；(2)根据等概率变换原则将U=(U1,U2)T映射为相关非正态变量X=(X1,X2)T。

1.3 建筑物风险评价和保护管理

许多既有建筑物存在一定程度的性能老化使得建筑承载能力降低。在隧道掘进过程中可能对现有建筑物造成一些潜在的结构损伤。为以最小的经济代价保证隧道掘进区既有建筑物的安全，必须找到工程安全和经济成本的平衡点，建立建筑物安全等级与保护措施之间的对应关系。

1.3.1 建筑物安全等级

由于地面位移会对建筑物造成较大的影响，地表沉降是隧道挖掘期的主要关注指标。参考文献[15, 16]中关于隧道掘进对建筑物的影响研究，本文将靠近建筑物附近的最大地表沉降Smax以及最大建筑物倾斜Imax作为隧道开挖对建筑物安全影响的评价指标，并对其进行相应分级，其中A等级对建筑物影响最小，D等级与之相反，如表2所示。

表2 建筑物安全等级标准

1.3.2 建筑物风险状态评估

确定建筑物的风险状态水平，有助于保护受地铁开挖影响下的既有建筑物。本文将系统的失效概率作为判断建筑风险状态的基础。为研究开挖地铁隧道对临近建筑物安全风险的影响，本文采用蒙特卡洛模拟法确定系统的失效概率，计算公式如下：

Pf=P(g1(Smax)<0∪g2(Imax)<0)

(9)

(10)

(11)

β=Φ-1(1-Pf)

(12)

根据标准[17]通过建筑物的重要程度来确定不同建筑物的可靠性指数β，如表3所示。

表3 建筑物重要度等级标准

构建系统化和结构化的风险评估程序，是实现系统风险准确评估的关键[18]。本文基于Eskesen等[19]提出的风险水平判断矩阵构建建筑物安全风险可靠性等级评价标准，其中在隧道掘进作用下建筑物的潜在风险分为四个不同等级，即“Ⅰ(安全)，Ⅱ(低风险)，Ⅲ(中风险)，Ⅳ(高风险)”，如表4所示。

表4 建筑物安全风险可靠性等级评价标准

2 指标体系的建立

本文依托武汉轨道交通二号线越江段隧道工程汉口段，根据盾构掘进方向自DK12+400(汉口江滩)—DK11+739(江汉路站接收井)，全长约660 m。隧道沿线6栋主要建筑物及其平面位置如图1所示。

图1 临近建筑物平面位置

2.1 建筑调查

建筑物在隧道掘进过程中所受影响属于建筑物属性与盾构施工的综合作用效果。本文用以下四个属性来体现6栋建筑的差异性：(1)建筑距地铁开挖轴线的水平距离；(2)地铁埋深；(3)建筑物层高；(4)建筑物结构形式。建筑物特征属性值如表5所示。

表5 各栋建筑物的特征属性

2.2 输入指标体系

在地铁盾构掘进过程中，隧道和建筑物之间的相互作用还受到许多环境因素的影响。因此，本文在模型中引入大量的输入变量以精确描述隧道和建筑物的相互作用机理。在参考相关专家经验和相关文献[15, 20, 21]的基础上，共选取了11个独立随机变量作为模型的输入参数，其分布形式均采用正态分布，各变量的均值、参数变异率见表6，而其他的输入参数则被假定为定值。

表6 输入变量分布形式

3 隧道建设有限元模型建立

本节建立了六个相关建筑物的确定性和随机有限元模型。在FE模型中考虑了周围环境因素，通过11个输入变量的变化，在计算结果中反映了隧道引起的建筑物位移的概率和可靠性。

3.1 确定性有限元模型建立

本文利用ABAQUS模拟隧道开挖对周围建筑物的影响，土体采用Drucker-Prager本构模型，混凝土和钢筋采用弹性本构。模型系数设置见表7。

表7 模型相应部件的材料属性

通过确定性有限元分析和随机有限元分析两个过程来构建如图2所示的3D有限元模型。得到共包含47940个单元，尺寸为50 m×40 m×60 m的模型。其中，盾构机长度、内径及厚度分别为9，6.2，0.4 m；衬砌管片的内径及管片厚度分别为6，0.3 m；注浆区域的内径和注浆加固土厚度分别为6.6，0.2 m；边界条件为固定xz底面单元，限制yz侧面的x方向位移，限制xy侧面的z方向位移。

图2 隧道有限元模型/m

以有限元模型中隧道掘进诱发的1#～6#建筑物地表位移来表征隧道开挖对建筑物的扰动情况，其具体计算为用有限元模拟结束时的总地表位移减去建筑物自重导致的位移差值。图3为1#建筑物下隧道掘进模型的模拟结果。

图3 有限元模型中的位移

3.2 随机有限元模型建立

本文将11个参数作为独立的随机变量输入到模型中，通过上述方法，对1#～6#建筑物逐一将有限元模型重复运行提取出180组输出参数(即Smax，Imax)。为保证有限元模型的鲁棒性，本文在输入参数导入过程中剔除其异常值。通过随机有限元模型，生成6个不同建筑物的Smax，Imax数据的盒须图，如图4所示。盒须图反映了随机FE模型的180组模拟结果，也是下一步基于Copula模拟所采用的数据基础。

图4 1#～6#建筑物的Smax与Imax盒须图

4 基于Copula理论的联合分布函数模拟

由于需要大样本的模拟数据来确保蒙特卡洛法的准确性(当故障概率的准确度为0.01%时，需要样本大小为104的数据)，显然这个数据量难以通过随机有限元法得到。因此，为实现对工程安全性的客观评价，本文基于Copula理论构建Smax，Imax的二维联合分布函数，并采用蒙特卡洛法生成模拟数据用于得到系统的失效概率和可靠度。另外，为分析不同Copula函数对模型的影响，本文利用AIC和BIC准则识别出最优Copula函数，并用其结果评价建筑物的安全性。

4.1 基于Copula理论联合分布函数构建

由于利用Copula理论构建的联合分布函数，需要确定模拟样本数据的边缘分布类型，本文在利用Copula理论构建Smax，Imax的联合分布函数之前，先对其边缘分布类型进行测试检验。通过正态检验发现，Smax和Imax的边缘分布符合Gaussian分布。

4.1.1 最优Copula函数的确定

最优Copula函数的识别可分为两个子步骤[9]。

Step 1：获取Copula参数。基于Smax，Imax的模拟样本，可以确定Copula参数。

对于Gaussian Copula，Kendall秩相关系数由式(3)计算，可以从MATLAB中使用τk=corr([Smax,Imax],‘type’,‘Kendall’)获得。然后，通过式(13)[9]计算Gaussian Copula参数ρ：

(13)

例如，基于1#建筑物中Smax，Imax的180个模拟数据，得到的τk为0.2912，表明Smax，Imax存在显著的相关性，因此不能忽略Smax，Imax之间的参数相关性，证明了利用Copula函数构建参数的二维联合分布函数的合理性。由式(13)，Gaussian Copula参数ρ计算为0.441631。对于其他三个Copula函数，表1中所示的相应Copula函数C(μ1,μ2)可以代入式(5)，并通过使用二分法[9]确定。获得三个相应参数θ，如表8所示。

表8 1#～6#建筑物的Copula参数计算值

Step 2：通过AIC和BIC准则识别最优Copula函数。在得到Copula参数之后，可以确定Copula函数及其密度函数的形式，如表1。此后，利用式(6)，(7)可得到AIC和BIC的计算值并用来识别出最优Copula函数。

基于Smax，Imax的180个模拟数据，可以获得所选四个Copula函数的AIC和BIC值，如表8所示，具有最小AIC和BIC值的Copula函数形式是最优Copula函数。从表8可以看出对于1#，2#，3#，6#建筑物，Gumbel Copula是最优Copula函数，对于4#，5#建筑物，Gaussian Copula是最优Copula函数。综上，对于不同的建筑物其Smax和Imax的依赖关系并不一致，因此在构建Copula函数模型时，需要在不同的Copula函中选择出最优Copula函数来计算。

4.1.2 二维联合分布函数模拟

如上所述，通过Copula方法可以生成数据样本，本文基于随机有限元模拟的Smax，Imax数据，分别用Clayton，Gaussian，Gumbel，Frank Copula函数构建Copula函数模型。

下面以1#建筑为例，给出不同Copula函数构造的Smax，Imax联合分布函数的数据模拟结果，如图5所示。图5为模拟次数N=1000时模拟的Smax，Imax的Copula函数生成数据(蓝色空心圆点)以及随机有限元数据(红色实心圆点)的散点图。可以看出：对于1#建筑，Copula函数模拟数据较好地包含了原始Smax，Imax数据点的分布。从AIC和BIC准则计算结果来看，Gumbel Copula是构建Copula函数模型的最优Copula函数。

图5 1#建筑物有限元模拟与Copula模拟对比

4.2 不同Copula函数计算结果比较

由于不同的Copula函数描述变量间依赖性的规律并不相同，因此本文从联合概率密度函数和条件累积分布函数这两个角度，对采用的四个Copula函数的计算结果进行比较分析。

4.2.1 联合概率密度函数

以1#建筑为例，图6给出四种Copula函数构造的Smax，Imax联合概率密度函数的曲线图。由于Copula函数的不同结构性质，在相同的Kendall秩相关系数(τk=0.1912)和相同的边缘分布函数的情况下，不同的Copula函数构造的联合概率密度函数不同，主要表现在两个方面：(1)不同的Copula函数构造的联合概率密度函数，其等概率密度线形状不同，Gaussian和Frank Copula函数是对称Copula函数，Clayton和Gumbel Copula是不对称的Copula函数；(2)最大联合概率密度函数值也不同，Clayton Copula函数具有最大的联合概率密度函数值，其次依次是Gumbel，Frank，Gaussian Copula函数。

图6 1#建筑物4种Copula函数联合概率密度函数曲线

4.2.2 条件累积分布函数

根据表2，当Smax，Imax的值较大时，建筑物的潜在风险较大。因此，对具有较高Smax，Imax值的条件累积分布函数需要特别关注。同时，为了定量比较不同Copula函数之间的差异，以1#建筑为例，在图7中比较了四个不同Copula函数相关的条件累积分布函数。

当Smax的取值范围为[μSmax+σSmax,μSmax+3σSmax]时，其中μSmax，σSmax分别为Smax的期望和方差，由Clayton，Gaussian，Gumbel，Frank Copula函数构造的条件累积分布函数得到的Imax值的期望分别为0.5172，0.4534，0.5170，0.5179，说明当Smax取较大值时，Imax的值也较大。当Imax的取值范围为[μImax+σImax,μImax+3σImax]时，其中μImax，σImax分别为Imax的期望和方差，由四个Copula函数构造的条件累积分布函数得到的Imax值的期望分别为0.0257，0.0234，0.0275，0.0268，说明当Imax取较大值时，Smax的值也较大。

图7 1#建筑物4种Copula函数获取条件累积分布函数

由图7可知，由Gaussian Copula函数构建的参数依赖关系与其他函数不同，利用AIC和BIC准则判断发现Gumbel函数为最优Copula函数，说明利用Copula函数构建参数相关性时，不同的Copula函数的结果可能相差很大；在传统的Nataf模型中，存在对变量的依赖结构服从Gaussian分布的潜在假设，而模拟结果证明这种假设是与客观事实存在差异的。一般来说，不同的Copula函数构建的相关依赖结构不同，表明寻找最佳拟合Copula函数来精确表征变量的依赖关系结构是非常重要的。

5 建筑物安全可靠性分析

建筑物安全状况分为如表4所示的四个等级，通过对每个建筑使用不同的Copula函数得出样本大小为106的蒙特卡洛模拟结果，计算出建筑物处于“B(正常)”“C(差)”“D(濒危)”状态的失效概率，如表9，表9还列出了当Smax，Imax的相互独立的结果用于对比分析。由表9可知随着建筑状态由从B变为D，失效概率逐渐减小，且Smax，Imax相互独立时的失效概率大于考虑参数相依性的失效概率。

表9 不同的Copula模型及相互独立下建筑的失效概率

根据AIC和BIC准则，对1#，2#，3#，6#建筑物，采用Gumbel Copula得到的失效概率作为参考结果，对4#，5#建筑物，采用Gaussian Copula的失效概率作为参考结果。

根据式(12)，建筑可靠性指标β可以通过失效概率来计算。根据建筑物的重要程度(表3)和建筑物安全风险可靠性等级评价标准(表4)，表10总结了6座建筑的风险等级。由表10可知，1#建筑的风险等级最高，其安全性需要在地铁施工中重点关注；其次是4#建筑，再次是3#，6#建筑，最次是2#，5#建筑。

表10 6座建筑物的风险评估结果

6 结 论

武汉地铁二号线线路附近存在6栋既有建筑物，需要研究地铁施工对其安全性的影响。本文提出了SFE-Copula模型，首先采用随机有限元模型构建了隧道-建筑相互作用的数值模型，然后通过Copula函数构建了Smax，Imax的联合分布函数，最后采用不同Copula函数模型生成的数据分析评估了这6栋建筑的安全等级。

在面对数据量不足的情况下，本文引入Copula函数来构建地铁盾构隧道施工对临近建筑物的安全可靠性的联合分布函数的构造问题，解决了使用传统方法对地下工程参数相依性建模时由于数据需求量不能满足的问题而影响其实际应用。由于Copula函数的引入，使得联合分布函数的精度提升，以用于评估邻近建筑在隧道挖掘期间的潜在风险，对于临近建筑物的安全保护问题具较高的参考意义和价值，实现以最小的经济代价提高地铁施工的安全性。

然而，本文中所使用的参数依赖于有限元模型的输出，由于其只是对真实物理系统进行的模拟，无法完全还原工程的实际情况，因此，结果存在与实际建筑物的安全可靠性不完全一致的可能性。在以后的研究中，若能获得实际的参数，将会对实际工程具有更高的指导意义。但本文依然为在数据量不足的情况下进行相依性建模的问题提供了一个新的解决思路。