Copula 应用（2）：优化藤 Copula 的桥梁结构系统地震易损性分析

---

引言

这是一篇发在《土木工程学报》上的文章，是高维Copula在土木工程方向的应用。复现见后文的代码。
为在桥梁结构系统易损性中更好地考虑不同构件个体之间的相关性，文章提出一种新的优化藤 Copula 分析桥梁结构系统易损性的方法。 该方法与现有的考虑构件类别之间相关性的桥梁系统易损性分析不同， 是将桥梁系统包含的各个构件的易损性作为藤Copula 函数中的一个变量。 通过不同构件个体的地震需求样本间的秩相关系数， 基于最短哈密顿路径方法筛选优化的藤 Copula 函数。 从而避免大量的联合概率密度函数的求解以及信息判定准则的试算。 同时， 提出优化藤Copula 函数是基于构件易损性补集的联合概率密度函数(PDF)， 通过其分析桥梁系统易损性仅需构建一个联合PDF， 可避免求解构件不同组合情况下较多的联合PDF。 选用四跨有桥台混凝土桥梁结构模型， 用桥梁剩余能力模型模拟桥梁的损伤， 用提出的方法分析具有桥墩桥台和支座共 10个构件个体相关性的桥梁系统易损性， 并与目前通常采用的考虑构件类别之间相关性的桥梁系统易损性分析以及Monte Carlo(MC)方法的结果对比， 验证提出方法的有效性。

---

基本原理

R语言代码片段

library(readxl)
library(ggraph)
library(readxl)
# 加载所需的库
library(raster)
library(rgdal)
library(ncdf4)

setwd(dir = "G:/data/")
#ata <- read.csv("中度损伤.csv", header = TRUE)

# 读取Excel文件
data <- read_excel("构件Pf汇总1.xlsx", sheet = "完全")

#先大概看下他们之间的相关系数（这里用kendall）
RR<-round(cor(data, method="kendall"), 3)#基本都是强正相关关系
R=1-RR
write.table(R, file = "corr完全.txt", sep = "\t", row.names = FALSE, col.names = FALSE)

a1<-fitalldist(data$a1)
a2<-fitalldist(data$a2)
a3<-fitalldist(data$a3)
a4<-fitalldist(data$b1)
a5<-fitalldist(data$b2)
a6<-fitalldist(data$c1)
a7<-fitalldist(data$c2)
a8<-fitalldist(data$S5)
a9<-fitalldist(data$S6)
a10<-fitalldist(data$S10)
a11<-fitalldist(data$M5)
a12<-fitalldist(data$M6)

a13<-fitalldist(data$A)
a14<-fitalldist(data$B)
a15<-fitalldist(data$C)


uvz <- cbind(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12) #12类构件 


stru<-muticop$structure
pair<-muticop$pair_copulas
muticopcdf<-pvinecop(uvz, muticop)
copulas<-summary(muticop)
Pfs=1-muticopcdf

write.table(Pfs, file = "PFS完全.txt", sep = "\t", row.names = FALSE, col.names = FALSE)
write.table(stru[["order"]], file = "stru完全.txt", sep = "\t", row.names = FALSE, col.names = FALSE)

copulas_df <- as.data.frame(copulas)
capture.output(copulas_df, file = "copulas完全.txt")

#write.table(copulas, file = "copulas轻微.txt", sep = "\t", row.names = FALSE, col.names = FALSE)

uvz1 <- cbind(a13,a14,a15,a8,a9,a10,a11,a12)  #3类构件


Pfs1=1-muticopcdf1
write.table(Pfs1, file = "PFS1完全.txt", sep = "\t", row.names = FALSE, col.names = FALSE)

plot(muticop)
contour(muticop)
#write.table(stru$struct_array, file = "stru_array.txt", sep = "\t", row.names = FALSE, col.names = FALSE)
Calc_Emp_Prob <- function(D) {
  # Length of data
  n = length(D)
  # Pre-assign probability array
  #P = zeros(n,1)
  P=numeric(n)
  # Loop through the data
  for (i in 1:n) {
    #P[i,1] = sum( D <= D[i] )
    P[i] = sum( D <= D[i] )
  }
  
  # Gringorten plotting position
  Y = (P - 0.44) / (n + 0.12)
  
}

fitalldist <- function(data,X1=100) {
  
  require(fitdistrplus)
  require(extRemes)
  require(ismev)
  require(SCI) 
  require(goft)
  require(gPdtest)
  require(actuar)
  require(evd)
  
  # change the same data value slightly
  
  # dup_data <- duplicated(data)
  # for (i in 1:length(dup_data)) {
  #   if (dup_data[i] == 0) {
  #     data[data=dup_data[i]] = data[data=dup_data[i]] + runif(data[data=dup_data[i]],0,0.001)
  #   } else {
  #     data[data=dup_data[i]] = data[data=dup_data[i]] + runif(data[data=dup_data[i]],-0.001,0.001)
  #   }
  # }
  
  Calc_Emp_Prob <- function(D) {
    # Length of data
    n = length(D)
    # Pre-assign probability array
    #P = zeros(n,1)
    P=numeric(n)
    # Loop through the data
    for (i in 1:n) {
      #P[i,1] = sum( D <= D[i] )
      P[i] = sum( D <= D[i] )
    }
    
    # Gringorten plotting position
    Y = (P - 0.44) / (n + 0.12)
    
  }
  
  if (all(is.na(data)) || all(data == 0)) {
    ep <- rep(NA, length(data))
  } else {
    # calculate empirical CDF
    emp <- Calc_Emp_Prob(data)
    
    # total 13 distributions
    distnames <- c("gamma","exponential","weibull","normal","logistic",
                   "log-normal","log-logistic","cauchy",
                   "generalized extreme value","generalized pareto",
                   "inverse gaussian")
    
    results <- data.frame('aic' = numeric(), 'aicc' = numeric(), 
                          'bic' = numeric(), 'p' = numeric())
    
    ep_mat <- matrix( NA, length(data), 11 )
    pdf_mat <- matrix(NA, length(data), 11 )
    V1_mat <- matrix(NA, length(X1), 11 )

结果和出图

如果需要更多关于Copula函数应用的联合重现期、同现重现期和条件重现期等的代码（MATLAB和R语言），如下图所示。详见链接（代码可直接运行，替换数据即可运行出结果）https://mp.weixin.qq.com/s/Eh982f_u0rXYoqEuW9-V8Q
欢迎关注公众号，获取更多科研代码和前沿论文资讯等相关内容

参考文献

[1]李幸钰,雷鹰,林树枝,等.优化藤Copula的桥梁结构系统地震易损性分析[J].土木工程学报,2023,56(11):43-55+136.DOI:10.15951/j.tmgcxb.22060665.
[2]王国庆. Copula函数及ANN方法在矮塔斜拉桥地震易损性分析中的应用[D].合肥工业大学,2023.DOI:10.27101/d.cnki.ghfgu.2022.001376.