随着社会经济的发展，城市化进程不断加快，在人口和财富高度聚集的大城市，地震灾害往往造成严重后果，如大量人员伤亡、巨额经济损失和城市功能中断等，前不久发生的土耳其地震就是惨痛案例。党的十九大报告中明确提出要“健全公共安全体系，提升防灾减灾救灾能力”。提高建筑抗震韧性已被纳入住房和城乡建设部“十三五”规划，其中明确提出了建筑抗震韧性的设计目标，引导建筑行业健康发展。建筑抗震韧性是指建筑在设定水准地震作用后，维持与恢复原有建筑功能的能力，近年来逐渐成为地震工程和结构工程学科的研究热点。
本期重点介绍建筑抗震韧性评价流程中的样本矩阵扩充方法。
初学韧性评价的小伙伴可能会有个疑问，直接计算每一次时程分析的修复费用、修复时间及人员伤亡率再拟合对数正态分布曲线不就行了吗，为啥要进行样本扩充呢，怪麻烦的。当然啦，如果时程分析的数目足够多，这样做是完全没问题的。然而，非线性动力时程分析往往需要耗费大量的时间，一般只做少量地震动的时程分析。如果样本数目太少，拟合修复费用、修复时间及人员伤亡率对数正态分布曲线误差会很大，甚至拟合失败。这就需要对原始样本进行扩充，得到大量的模拟样本。图1为某工程20组地震动时程分析修复费用的概率分析，图2为将原始样本扩充为1000个样本后的修复费用概率分析，可见原始样本修复费用对数正态分布曲线拟合效果不太好，误差很大，而1000个模拟样本的修复费用对数正态分布曲线拟合效果很好。
图1 原始样本概率分析
图2 模拟样本概率分析
小伙伴们可能又会突发奇想，能不能直接根据原始时程分析修复费用的概率分布特征直接生成大量模拟样本呢，那样岂不方便？答案是NO，如果这样做就没有意义了，模拟样本的概率分布特征和原始样本完全一样。样本扩充必须在更早的阶段完成，由构件损伤状态计算修复费用、修改时间及人员伤亡率的方法还是比较复杂的，且多以表格形式给出，为了覆盖更多可能的路径，需要采用蒙塔卡罗方法。
《建筑抗震韧性评价标准》（GB /T 38591-2020）给出的流程是：提取时程分析工程需求参数（Engineering Demand Parameter，EDP，如层间位移角、楼面加速度等），扩充工程需求参数样本矩阵，再通过易损性数据库确定构件损伤状态，进而计算修复费用、修复时间及人员伤亡率。以层间位移角、楼面加速度作为工程需求参数进行样本扩充，简便易行，对于规则结构取得了很好的效果，极大地促进了抗震韧性评价方法在实际工程中的应用。然而，对于存在楼层变形不均匀、扭转不规则等情况的结构，同一楼层不同位置的构件变形是不同的，通过一个统一的层间位移角难以准确反映同楼层所有构件的变形。且对于没有层概念的空间结构，无法定义层间位移角。为了提升构件损伤状态判断的准确性，SAUSG-RES软件直接基于非线性时程分析结果判断构件损伤状态，再对构件损伤状态矩阵进行扩充，进而计算修复费用、修复时间及人员伤亡率。
无论是扩充工程需求参数矩阵，还是扩充构件损伤状态矩阵，样本扩充的算法是类似的，以下以工程需求参数矩阵扩充为例来详细介绍样本矩阵扩充算法。
每次时程分析可得到一组工程需求参数，将各次时程分析的工程需求参数组装为原始工程需求参数矩阵，每行表示一次时程分析的结果，每列表示一个特定的工程需求参数的取值。扩充后的工程需求参数与分析得到的工程需求参数应具有相同的联合分布，即扩充后的工程需求参数矩阵与时程分析得到的工程需求参数矩阵应具有相同的均值与方差。
图3 矩阵扩充过程
（1）
式中：U=[U1, U2…, Un]T是满足独立正态分布的 n 维向量；A为 n×n 的常数矩阵；B 为n×1的常数矩阵。
由均值矩阵与协方差矩阵的性质可知：
（2）
（3）
式中：MZ为Z的均值矩阵；ΣZZ为Z的协方差矩阵；MU 为 U 的均值矩阵；ΣUU为U的协方差矩阵。
由于U满足独立正态分布，故MU=0，ΣUU=I，代入式（2）和式（3）得：
由于扩充后的工程需求参数矩阵与原始工程需求参数矩阵具有相同的均值与方差，MZ =MY，ΣZZ=ΣYY，则式（1）转化为：
式中：当∑YY满秩时，L=(chol(∑YY)) T，其中chol()为cholesky分解。
独立正态分布变量的取值可由计算机的伪随机数生成。利用式（6）计算得到一组工程需求参数的对数取值，取指数后得到工程需求参数的取值，多次利用该式，将计算结果组装为矩阵，即可得到扩充后的工程需求参数分布的取值。
（7）
式中：p为∑YY的秩；Dpp通过分割协方差矩阵∑YY特征值的平方根矩阵得到；Lnp通过分割协方差矩阵∑YY特征向量矩阵得到。
对于多数工程，采用SAUSG-RES软件进行韧性评价时，构件数目远大于时程分析的地震动数目，且多数构件处于无损坏状态，构件损伤状态向量线性相关，协方差矩阵非满秩，Cholesky分解算法将变得不稳定，因此，采用FEMA P-58改进后的方法进行样本矩阵扩充，代码如下：
clear all;close all;clc;
clear;
num_realization = 1000;
EDPs =load('EDP.txt');
% 需求矩阵取自然对数
lnEDPs=log(EDPs);
[num_rec num_var]=size(lnEDPs)
%B=load('beta.txt'); %模型对数标准差βm 和地震动对数标准差βgm
B=zeros(2,num_var);
% 计算自然对数矩阵lnEDPs均值
lnEDPs_mean=mean(lnEDPs);
lnEDPs_mean=lnEDPs_mean';
% 计算自然对数矩阵 lnEDPs协方差矩阵
lnEDPs_cov=cov(lnEDPs);
% 计算协方差矩阵的秩
lnEDPs_cov_rank=rank(lnEDPs_cov);
% 放大方差
% sigma=sqrt(diag(lnEDPs_cov));
% R=lnEDPs_cov ./ (sigma*sigma');
% B=B';
% sigmap2=sigma.*sigma;
% sigmap2=sigmap2+(B(:,1).* B(:,1)); % 考虑βm
% sigmap2=sigmap2+(B(:,2).* B(:,2)); % 考虑βgm
% sigma=sqrt(sigmap2);
% lnEDPs_cov_inflated=R.*(sigma*sigma');
lnEDPs_cov_inflated=lnEDPs_cov;
% 计算协方差矩阵特征值
[L_total D2_total]=eig(lnEDPs_cov_inflated);
D2_total=eig(lnEDPs_cov_inflated);
% 分割 L_total 得 L_use，L_use是L_total的一部分
if lnEDPs_cov_rank >= num_var
L_use =L_total;
else
L_use =(L_total(:,num_var- lnEDPs_cov_rank+1:num_var));
end
% 分割 D2_total 得 D2_use，D2_use是D2_total的一部分
if lnEDPs_cov_rank >= num_var
D2_use =D2_total;
else
D2_use =D2_total(num_var- lnEDPs_cov_rank+1:num_var);
end
% 计算D2_use的平方根
D_use =diag((D2_use).^0.5);
% 生成标准正态分布随机数
if lnEDPs_cov_rank >= num_var
U = randn(num_realization,num_var) ;
else
U = randn(num_realization, lnEDPs_cov_rank) ;
end
U = U' ;
%计算Z矩阵
Lambda = L_use * D_use ;
Z = Lambda * U + lnEDPs_mean * ones(1,num_realization);
lnEDPs_sim_mean=mean(Z');
lnEDPs_sim_cov=cov(Z');
A=lnEDPs_sim_mean./lnEDPs_mean';
B=lnEDPs_sim_cov./lnEDPs_cov;
W=exp(Z);
W=W';
save EDP_Extend.txt -ascii W
%end
参考文献：
[1] GB /T 38591-2020, 建筑抗震韧性评价标准[S]. 北京: 中国标准出版社, 2020.
[2] Yang T Y, Moehle J, Stojadinovic B, et al. Seismic Performance Evaluation of Facilities: Methodology and Implementation [J]. Journal of Structural Engineering, 2009, 135(10): 1146-1154.
[3] FEMA. Seismic performance assessment of buildings: Volume 1-Methodology, Second Edition: FEMA P-58 [S]. Washington DC: Federal Emergency Management Agency, 2018.
[4] 潘鹏教授解读GB/T 38591—2020《建筑抗震韧性评价标准》. 清华大学互联网产业研究院. http://www.iii.tsinghua.edu.cn/info/1059/2520.htm
[5] 申洲洋, 潘鹏, 黄世敏. 基于抗震设防的建筑抗震韧性评价方法
联系客服
