Introducción

El análisis probabilista de peligro sísmico es un procedimiento que permite evaluar la tasa media de excedencia de medidas de intensidad del movimiento sísmico del suelo. Es una herramienta ampliamente utilizada en la ingeniería sísmica y se ha aplicado en la elaboración de mapas de peligro sísmico asociados a un cierto periodo de retorno y en el cálculo de espectros de peligro uniforme. ^{Esteva (1967)} y ^{Cornell (1968)} concibieron originalmente el marco metodológico del análisis probabilista de peligro sísmico incorporando las incertidumbres asociadas con la magnitud y localización de las fuentes sísmicas, así como con los parámetros que caracterizan el movimiento del suelo en un sitio de interés. Empleando el teorema de la probabilidad total, desarrollaron una formulación univariada para calcular la probabilidad o la tasa media de excedencia de medidas de intensidad escalares tales como la aceleración máxima del suelo o la aceleración espectral. Dichas formulaciones se expresan en términos de la probabilidad condicional de excedencia de una medida de intensidad escalar dada una cierta magnitud y una localización de la fuente sísmica, y de la distribución de probabilidad marginal de éstas.

La tasa de excedencia de las medidas de intensidad del movimiento del suelo se emplea en el análisis probabilista de demanda sísmica, el cual tiene por objeto estimar la tasa de excedencia de la respuesta estructural (ver e.g. ^{Shome, 1999}; ^{Moehle y Deierlein, 2004}; ^{Conte y Zhang, 2007}; ^{Barbosa, 2011}). En el análisis probabilista de demanda sísmica se considera que una medida de intensidad eficiente es aquella que permite caracterizar con menor dispersión la respuesta estructural y, por lo tanto, estimar con mejor precisión la probabilidad de excedencia de la respuesta condicional a dicha medida de intensidad (^{Luco y Cornell, 2007}). Se ha encontrado que el uso de una medida de intensidad escalar, como por ejemplo la aceleración espectral para el periodo fundamental, resulta ser un estimador ineficiente en ciertos casos (^{Bazzurro y Cornell, 2002}; ^{Baker y Cornell, 2005}; ^{Baker, 2007}; ^{Faggella et al., 2013}). Esto ocurre cuando la respuesta estructural está correlacionada significativamente con un conjunto de medidas de intensidad del movimiento sísmico, como puede ser el caso de edificios altos, puentes de gran longitud y plataformas marinas cuya respuesta no necesariamente está gobernada por el modo fundamental, o el de construcciones cercanas a una fuente sísmica y puentes cuya respuesta puede depender del movimiento sísmico del terreno en dos o tres direcciones ortogonales, además de los efectos de amplificación de sitio y de la duración de la excitación, o de estructuras que exhiban una respuesta altamente no-lineal. La tasa de excedencia de la respuesta estructural debe ser calculada en tales casos utilizando tasas de excedencia de medidas de intensidad vectoriales. El análisis probabilístico multivariado de peligro sísmico provee el marco metodológico para la estimación de la tasa de excedencia de medidas de intensidad vectoriales considerando su distribución de probabilidad conjunta. A diferencia de un enfoque univariado, en el que las tasas de excedencia de un conjunto de medidas de intensidad, como por ejemplo las aceleraciones espectrales para distintos periodos estructurales, son calculadas por separado sin tomar en cuenta la dependencia estadística entre ellas, el cálculo de la tasa de excedencia de medidas de intensidad vectoriales toma en cuenta, mediante una formulación multivariada, su distribución de probabilidad conjunta.

^{Bazzurro (1998)} propuso una metodología para el análisis multivariado de peligro sísmico. Esta metodología es una extensión del enfoque univariado y, aun cuando se puede generalizar a un cierto número de medidas de intensidad, fue ilustrada para el caso bivariado de dos aceleraciones espectrales. La densidad de probabilidad bivariada se expresa como el producto de la densidad marginal de una aceleración espectral y la densidad condicional de la segunda aceleración dada la primera. La metodología supone que la distribución conjunta de los logaritmos de las aceleraciones espectrales es de tipo normal o gaussiana, por lo que la densidad de probabilidad condicional y las marginales son desde luego de tipo log-normal. Esta suposición implica adicionalmente que la dependencia estadística entre los logaritmos de las aceleraciones espectrales puede ser caracterizada apropiadamente mediante el coeficiente de correlación lineal. Bajo esta perspectiva, algunos trabajos se han dedicado al análisis y modelado de la correlación lineal entre aceleraciones espectrales usando datos de registros sísmicos (^{Baker y Jayaram, 2008}; ^{Jayaram et al., 2011}). La metodología de ^{Bazzurro (1998)} fue empleada por ^{Bazzurro y Cornell (2002)} para evaluar la tasa de excedencia de la distorsión máxima de entrepiso de un modelo de edificio de acero de 20 niveles (^{Shome y Cornell, 1999}; ^{Gupta y Krawinkler, 1999}). Los resultados obtenidos mostraron una reducción de la tasa de excedencia en comparación con la calculada mediante el enfoque univariado. Aplicando la formulación de ^{Bazzurro (1998)}, ^{Rajeev et al. (2008)} evaluaron un procedimiento para seleccionar la combinación óptima de aceleraciones espectrales en el periodo del modo fundamental y en un modo superior de manera que se minimice la incertidumbre en la estimación de la tasa de excedencia aplicando una técnica de bootstrap a resultados de análisis de distorsión de columnas en un edificio de concreto reforzado. Con el fin de simplificar el análisis multivariado y de facilitar su aplicación en la práctica, se han propuesto formulaciones para evaluar de manera indirecta el peligro sísmico de medidas de intensidad vectoriales a través de la combinación de resultados de análisis univariados de medidas de intensidad escalares (^{Bazzurro et al. 2009}; ^{Koranghi et al., 2016}).

En este trabajo se presenta una formulación multivariada para el análisis probabilístico de peligro sísmico de medidas de intensidad vectoriales usando cópulas. La formulación se basa en la propuesta por ^{Bazurro (1998)}, sin embargo, a diferencia de aquella, la tasa de excedencia se expresa en forma genérica en términos de la cópula de una medida de intensidad vectorial, sin suponer un tipo de distribución de probabilidad multivariada en particular. Las cópulas son funciones que contienen la información sobre la estructura o patrón de dependencia estadística de las componentes de un vector aleatorio y definen su distribución de probabilidad conjunta (^{Nelsen, 2006}). La ventaja es que, en lugar de suponer a priori una estructura de dependencia particular y estimar los parámetros que la caracterizan, por ejemplo, el coeficiente de correlación lineal en el caso gaussiano, este enfoque permite estimar la cópula que representa apropiadamente la estructura de dependencia del vector aleatorio y, en conjunto con las distribuciones marginales, definir el modelo de distribución multivariada. Las primeras aplicaciones de la teoría de cópulas se produjeron en disciplinas como las matemáticas financieras y la actuaría (^{Embrechts et al., 2002}). Posteriormente, se han aplicado también en la ingeniería civil para la caracterización multivariada de fenómenos naturales (^{Salvadori et al., 2007}), el cálculo de contornos de peligro ambiental para el diseño de estructuras marinas (^{Montes-Iturizaga y Heredia-Zavoni, 2015}, ²⁰¹⁶, ²⁰¹⁷), el análisis de demanda sísmica de múltiples parámetros de respuesta en edificios de concreto reforzado (^{Goda y Tesfamariam, 2015}) y estudios de confiabilidad estructural (^{Tang et al., 2013a}, ^b).

El presente artículo inicia con una exposición de los conceptos básicos para el modelado multivariado empleando cópulas. Estos conceptos se utilizan luego para desarrollar la formulación multivariada del análisis probabilístico de peligro sísmico de medidas de intensidad vectoriales. La formulación se aplica a un ejemplo de análisis de la tasa de excedencia de la distorsión máxima de entrepiso usando el modelo de regresión de ^{Bazzurro y Cornell (2002)} para un edificio de acero de 20 niveles. En dicho análisis, se consideró una medida de intensidad vectorial compuesta por las aceleraciones espectrales para los dos primeros periodos modales de la estructura. Se presenta el procedimiento detallado de estimación de la cópula del vector de aceleraciones espectrales mediante estimadores basados en el τ de Kendall (^{Genest y Favre, 2007}), el cual es una medida de asociación que permite caracterizar apropiadamente el patrón de dependencia entre variables aleatorias. Para la estimación de la cópula, se utilizaron datos de registros sísmicos de la zona de transición de la Ciudad de México. Este tipo de excitación genera una contribución importante de la respuesta asociada con el segundo modo de la estructura. Se describe enseguida el modelo de ley de atenuación utilizado para definir las distribuciones de probabilidad marginales de las aceleraciones espectrales. Se analizan los resultados del cálculo de la tasa de excedencia empleando la formulación multivariada y se examinan las diferencias con los resultados de un análisis univariado que sólo considera la aceleración espectral para el modo fundamental de la estructura. El trabajo concluye con unos comentarios finales y se apuntan algunos temas de interés para posibles trabajos futuros.

Modelado multivariado empleando cópulas

Las cópulas son funciones que definen rigurosamente la distribución de probabilidad conjunta o multivariada de un grupo de variables aleatorias 𝑋=(𝑋₁,…,𝑋_𝑛). Estas funciones permiten ensamblar la distribución conjunta 𝐹_𝑋(𝑥₁,…,𝑥_𝑛) a partir de las distribuciones marginales 𝐹_𝑖(𝑥_𝑖), 𝑖=1,⋯,𝑛,, de cada una de las variables de la siguiente manera (^{Nelsen, 2006}),

donde, 𝐶(𝑢₁,…𝑢_𝑛) es la cópula y 𝑈_𝑖= 𝐹_𝑖(𝑋_𝑖) son variables uniformemente distribuidas en [0,1] La ecuación 1 se conoce como el teorema de ^{Sklar (1959)}). Dada la distribución multivariada, 𝐹_𝑋(𝑥₁,…,𝑥_𝑛),, se puede deducir entonces la siguiente expresión para la función de densidad de probabilidad conjunta,

en la que fi(xi) es la función de densidad de probabilidad marginal de 𝑋_𝑖 y 𝑐(𝑢₁,…𝑢_𝑛) es la densidad de la cópula,

Tal como se observa en la ecuación 2, la densidad de probabilidad conjunta es igual a la densidad de la cópula multiplicada por el producto de las funciones de densidad marginales de las variables. Por lo tanto, toda la información sobre el patrón de dependencia estadística o estructura de asociación entre las variables está contenida íntegramente en la cópula puesto que las densidades marginales, fi(xi) , que caracterizan individualmente a cada una de las variables, no contienen dicha información.

Existe una gran variedad de tipos o familias de cópulas. Entre las familias comúnmente utilizadas en diversas disciplinas se encuentran las cópulas arquimedianas, tales como las cópulas de Frank, Gumbel y Clayton, reconocidas por la facilidad con que pueden ser construidas y por la gran variedad de patrones de dependencia que permiten reproducir (^{Salvadori et al., 2007}), así como la cópula Gaussiana con la que se pueden modelar estructuras de dependencia definidas mediante matrices de correlación lineal. Estas familias de cópulas se presentan en la tabla 1 para el caso bivariado en que 𝑛=2..

Para la estimación estadística de los parámetros de las cópulas se cuenta con métodos de máxima verosimilitud (^{Genest y Favre, 2007}). En el caso de cópulas de un sólo parámetro, como las cópulas bivariadas de Frank, Gumbel, Clayton y Gaussiana mostradas en la tabla 1, se pueden usar estimadores definidos en términos de medidas de asociación como el 𝜏 de Kendall o el p de Spearman (^{Lebrun y Dutfoy, 2009}). Estas medidas son estimadores cuantitativos del grado de asociación estadística entre dos variables y cumplen con las siguientes propiedades: (1) varían en el intervalo [-1,1]; (2) son iguales a cero si las variables son estadísticamente independientes; y (3) son invariantes bajo transformaciones monótonas de las variables. El coeficiente de correlación lineal, que ha sido utilizado ampliamente en diversas disciplinas, principalmente debido a la facilidad con que puede ser estimado y por su utilidad en la construcción de modelos multivariados de tipo gaussiano, no es rigurosamente una medida de asociación. Como se sabe, dicho coeficiente es una medida de la dependencia lineal entre dos variables y permite modelar distribuciones conjuntas de tipo gaussiano en las que la estructura de dependencia se puede caracterizar apropiadamente mediante la correlación lineal. Sin embargo, no cumple con la tercera propiedad de las medidas de asociación. El cálculo del coeficiente de correlación lineal involucra el uso de las distribuciones marginales de probabilidad de las variables, las cuales no contienen información sobre la asociación estadística entre éstas. Por otro lado, una de las propiedades de las cópulas es que son invariantes ante transformaciones monotónicas de las variables. Dado que como se observa en la ecuación 2, la estructura de asociación o de dependencia entre las variables depende exclusivamente de la cópula, una medida apropiada de dicha asociación o dependencia debe exhibir también la misma propiedad invariante de las cópulas.

En este trabajo se utiliza el 𝜏 de Kendall como medida de asociación para la estimación de las cópulas. El 𝜏 de Kendall se define como la diferencia entre las probabilidades de concordancia y de discordancia de dos vectores estadísticamente independientes e idénticamente distribuidos, (𝑋₁, 𝑌₁) y (𝑋₂, 𝑌₂) ,

Si (X _i , Y _i ), i=1,…,n es una muestra aleatoria de un par de variables (X,Y), el 𝜏 de Kendall se estima como,

donde, P _n y Q _n y son el número de pares de datos concordantes y discordantes en la muestra (^{Genest y Favre, 2007}). Sean (R _i ,S _i ) ,i=1,…,n, los rangos de las observaciones en la muestra aleatoria (X _i , Y _i ). Puesto que (X _i - X _j ) (Y _i , Yj) > 0 sí y solo sí (R _i - R _j ) (S _i , Sj) > 0 , el 𝜏 de Kendall es entonces una función de los rangos y del mismo modo que éstos, se mantiene invariante bajo transformaciones monótonas de las variables. Usando la ecuación 5, se puede demostrar que cuando n→∞ , τ^ es un estimador asintótico centrado del 𝜏 de Kendall dado por ^{Genest y Favre, (2007)},

donde 𝐶(𝑢,𝑣) es la cópula del par de variables (𝑋, 𝑌), 𝑢=𝐹_𝑋(𝑥), 𝑣=𝐹_𝑌(𝑦).. Mediante la ecuación 6, se pueden obtener relaciones funcionales entre el ?? de Kendall y el parámetro θ de la cópula, tal como se muestran en la tabla 1. El estimado del parámetro de una cópula se puede entonces obtener como θ^=g(τ^) , donde θ=g(τ) es la relación funcional que se muestra en la tabla 1. ^{Genest y Favre (2007)} han desarrollado formulaciones para calcular los estadísticos (media y varianza) del estimador θ^=gτ^.

Análisis multivariado de peligro sísmico

Consideremos un 𝑋=(𝑋₁,…,𝑋_𝑛) de medidas de intensidad del movimiento sísmico del suelo requeridas para caracterizar una demanda de capacidad o de respuesta de una estructura. De acuerdo con la ecuación 2, la función de densidad de probabilidad fXx1,…xn|m,r se puede escribir en términos de la densidad de la cópula correspondiente y de las densidades marginales condicionales 𝑀=𝑚 y 𝑅=𝑟

donde, ui=Fixi|m,r 𝑖=1,⋯,𝑛, y Fixi|m,r y fixi|m,r son las distribuciones y las densidades de probabilidad condicionales de las medidas de intensidad, respectivamente. Nótese que no solamente los argumentos de la densidad condicional de la cópula, cu1,…un|m,r , dependen de M y R, sino también la propia función.

La densidad de la tasa media de eventos (𝑥₁,…𝑥_𝑛) considerando una fuente sísmica se define como (^{Bazurro, 1998}; ^{Bazzurro y Cornell, 2002}),

donde, λo es la tasa de ocurrencia de eventos de magnitud mayor que un límite inferior. La ecuación 8 es una extensión de la expresión propuesta por Esteva (1970) para el análisis univariado. Sustituyendo la ecuación 7 en la ecuación 8, se obtiene,

Cuando se consideran N fuentes sísmicas, la densidad de la tasa media de eventos (𝑥₁,…𝑥_𝑛) se calcula mediante la suma de las contribuciones de cada una de ellas,

donde, λi , 𝑖=1,⋯,𝑁,, son las tasas de ocurrencia de sismos de magnitud mayor que una cierta magnitud mínima en cada fuente.

Tal como se muestra en la ecuación 9, la caracterización multivariada de las medidas de intensidad se introduce en el análisis a través de la densidad de la cópula cu1,…un|m,r condicional a M y R. Supongamos ahora que es de interés evaluar la tasa media de excedencia de una cierta demanda de capacidad o de respuesta sísmica γ de una estructura. Si Pγ>z|x1,…,xn denota la probabilidad de excedencia, condicional (𝑥₁,…,𝑥_𝑛), dicha tasa es

Consideremos como ejemplo el caso de la estimación de la tasa media de excedencia de la distorsión máxima de entrepiso, γmax , de un edificio alto de periodo largo en el que la contribución del segundo modo de vibrar sea relevante. 𝑆_𝑎1 y 𝑆𝑎2 las seudo-aceleraciones espectrales en los dos primeros periodos de vibración de la estructura, T1 y T2 . De acuerdo con la ecuación 7, la densidad de probabilidad conjunta de Sa1 y Sa2 , condicional a M=m y R=r es igual a

En la ecuación 12, cSau, v|m,r denota la densidad de la copula condicional de las seudo-aceleraciones Sa1 y Sa2 ,

Y FSa1a1|m,r , FSa2a2|m,r , fSa1a1|m,r , fSa2a2|m,r , son las funciones de distribución y de densidad de probabilidad marginales de dichas variables condicionales a M=m y R=r , respectivamente. La densidad de la tasa media de eventos de seudo-aceleraciones espectrales a1,a2 está dada por:

En caso de que se consideren N fuentes sísmicas, νSaa1,a2 está dada por

donde, λi es la tasa de ocurrencia de eventos correspondiente a cada fuente. La tasa de excedencia de la distorsión máxima de entrepiso es entonces,

Ejemplo de aplicación

La formulación multivariada se aplicó con fines de ilustración a la estimación de la tasa media de excedencia de la distorsión de entrepiso máxima en edificios. Se consideró un modelo de regresión de la distorsión de entrepiso máxima de un edificio de acero de 20 niveles en función de las aceleraciones espectrales asociadas con las frecuencias de los dos primeros modos de vibrar (^{Bazurro y Cornell, 2002}). Para los fines de este ejemplo, el edificio se sometió a un tipo de excitación que genere una contribución importante de la respuesta asociada con el segundo modo y pueda entonces apreciarse la diferencia entre los resultados del análisis multivariado y del univariado.

El modelo de ^{Bazurro y Cornell (2002)} fue desarrollado mediante análisis dinámicos no-lineales paso a paso de un marco de acero del edificio de 20 niveles que se muestra en la figura 1, el cual fue empleado por ^{Shome y Cornell (1999)} en un estudio de demanda símica. Se trata de un marco interior de un edificio con planta regular que tiene cinco crujías cuyos claros son de 6.10 m; la altura del primer entrepiso es de 5.50 m y las alturas de los entrepisos superiores son de 4 m. Las frecuencias y periodos de los dos primeros modos de vibración son f1 = 0.25 Hz ( T1 = 4 s) y f2 = 0.75 Hz ( T2 = 1.33 s); para mayor información sobre la estructura se puede consultar el trabajo de ^{Gupta y Krawinkler (1999)}.

Se usaron registros de sismos de subducción de estaciones ubicadas en la zona II de transición de la Ciudad de México. Los nombres de las estaciones y la relación de sismos considerados se indican en la tabla 2. Los registros corresponden a sismos con magnitudes M _w entre 5 y 8.1, y distancias epicentrales comprendidas entre 283 km y 421 km. Los espectros de respuesta de los registros seleccionados se muestran en la figura 2 para una fracción de amortiguamiento crítico igual a 5%. En dicha figura se indican los dos primeros periodos del modelo de edificio de acero. Dada la forma que guardan los espectros en la zona de transición, se observan mayores aceleraciones espectrales para la segunda forma modal del edificio en relación con el modo fundamental; el periodo de la segunda forma modal se encuentra en la región de mayor demanda de aceleración de la zona de transición de la Ciudad de México. Para el modo fundamental de vibración del edificio se pueden alcanzar aceleraciones espectrales de hasta 25 gal, mientras que para la segunda forma modal éstas son superiores a los 100 gal.

Para la caracterización multivariada del peligro sísmico se utilizaron los datos de las aceleraciones espectrales Sa1 y Sa2 , asociadas con los periodos modales T1 y T2 , obtenidas de los registros. La gráfica de dispersión de los rangos normalizados de las aceleraciones espectrales Sa1 y Sa2 se muestra en la figura 3. Se observa que la relación de dependencia entre ambas variables no es estrictamente lineal. Con base en esta muestra se obtuvo un estimado del τ de Kendall igual a 0.622 usando el paquete Copula de R (^{R Core Team, 2020}; ^{Hofert et al., 2020}). Este valor representa un grado de asociación significativo entre las variables. El valor positivo indica que cuando una de las variables se incrementa (o disminuye) es más probable que la otra también se incremente (o disminuya).

Los parámetros de las cópulas de Frank, Gumbel, Clayton y Gaussiana se estimaron en función del valor estimado del τ de Kendall mediante las relaciones funcionales mostradas en la tabla 1. La tabla 3 muestra los valores de los parámetros estimados y sus correspondientes desviaciones estándar. Los coeficientes de variación de los estimados son del orden de 10 a 20%. Estos valores dependen, como en cualquier problema de inferencia estadística, del tamaño de la muestra de datos y, desde luego, se reducirían en caso de contar con un mayor número de registros. Para los fines de este estudio se considera que dichos valores son aceptables. Se realizaron pruebas de bondad de ajuste para cópulas usando el estadístico de Cramér-von Mises Sn  basado en la cópula empírica y el estadístico SnB basado en la transformación de Rosenblatt. Una descripción detallada de estos estadísticos se presenta en el apéndice A. Los valores p de las pruebas de bondad se obtienen usando un método de bootstrap propuesto por ^{Genest et al. (2009)}. Para el cálculo de los estadísticos y de los valores p de estas pruebas se usó el paquete Copula de R (^{R Core Team, 2020}; ^{Hofert et al., 2020}). Los resultados de las pruebas de bondad de ajuste se muestran en la tabla 3. El modelo de Gumbel proporciona en ambas pruebas los menores valores de los estadísticos y en consecuencia los mayores valores p. Si bien el modelo Gaussiano no se rechaza, presenta valores p menores que los de las copulas de Gumbel y de Frank. El modelo de Clayton se rechaza con un nivel de significación igual a 5%.

Se aplicaron adicionalmente los criterios de Akaike y bayesiano para selección de cópulas, los cuales están basado en teoría de la información (^{Akaike, 1973}; ^{Schwarz, 1978}). Estos criterios permiten comparar puntajes entre diferentes modelos y se selecciona aquel que resulta en un menor puntaje. La formulación para la aplicación de estos criterios se presenta a en el Apéndice A. Los criterios AIC y BIC (por sus siglas en inglés) miden la diferencia entre la parsimonia de un modelo, representada por el número de parámetros del modelo, y su bondad de ajuste, representada por la función de verosimilitud. Los criterios favorecen aquellos modelos con mayor bondad de ajuste y penalizan aquellos que tienen un mayor número de parámetros. Se presupone que un mejor modelo es aquel que con menos parámetros puede representar mejor los datos. Dada la formulación lineal con la que se definen estos criterios (ver Apéndice A), el puntaje que resulta de estos criterios puede ser un valor positivo o negativo dependiendo del valor de la función de verosimilitud y del número de parámetros. Independientemente de esto, el criterio de selección consiste en escoger, entre un grupo de modelos candidatos, aquel que tiene el menor puntaje. Dado que en este trabajo los modelos de cópulas tienen sólo un parámetro, la comparación entre los modelos depende exclusivamente del valor de la función de verosimilitud. Los resultados que se muestran en la tabla 4 indican que con el modelo de Gumbel se obtiene el menor puntaje. Los modelos de Frank y Gaussiano obtienen valores similares sin que se muestre una clara diferencia a favor de uno de ellos. Claramente el modelo de Clayton no resulta favorecido, observación que coincide con los resultados de las pruebas de bondad de ajuste. Se seleccionó por tanto el modelo de Gumbel como el que mejor representa el patrón de dependencia entre los datos de aceleraciones espectrales. La densidad de la cópula de Gumbel se muestra en la figura 4. En la figura 5 se muestran iso-contornos de la cópula de Gumbel comparados con la cópula empírica de los datos de aceleraciones espectrales. La cópula empírica sólo depende de los datos de una muestra; se define en términos de un estadístico de orden y es una representación no paramétrica del patrón de dependencia entre las variables. La comparación con la cópula empírica que se muestra en la figura 5 sugiere que el modelo de Gumbel seleccionado representa los datos adecuadamente.

Las distribuciones marginales de las aceleraciones espectrales se modelaron con base en la ley de atenuación espectral estimada por ^{Reyes et al. (2002)} y el uso de cocientes espectrales de respuesta. Empleando registros de sismos de subducción de la estación CU en suelo firme, ^{Reyes et al. (2002)} estimaron mediante un método de regresión bayesiana los coeficientes de la ley de atenuación espectral,

basada en la propuesta por ^{Joyner y Boore (1988)}; SarT es la aceleración espectral en suelo firme en cm/s², T es el periodo estructural en segundos, m es la magnitud de momento, r es la distancia más corta del sitio al área de ruptura en km y αi(T) son los coeficientes del modelo. La razón de emplear esta ley de atenuación, obedeció a que está calibrada para obtener aceleraciones espectrales en roca (estación CU) para sismos de la fuente de subducción, lo cual sirve de base para obtener la aceleración espectral tanto en suelo blando como en suelo de transición a través de cuantificar el efecto de sitio mediante un cociente espectral de respuesta. En este trabajo se consideró que R= 300 km y se adoptaron los estimados de los coeficientes del modelo para el componente N-S, los cuales se muestran en la tabla 5 para los dos periodos modales. Las desviaciones estándar del modelo de ley de atenuación son σln⁡SarT1=4= 0.603 y σln⁡SarT2=1.33= 0.495. En el modelo de ^{Reyes et al. (2002)}, la desviación estándar no es función de m y r , y depende solamente del periodo estructural.

La distribución de probabilidad de las aceleraciones espectrales en la zona de transición se obtuvo a partir de la distribución en suelo firme mediante el uso de cocientes espectrales de respuesta. Un cociente espectral es la relación entre un espectro de respuesta en un suelo blando y el espectro de respuesta en suelo firme. Los cocientes espectrales se pueden estimar mediante espectros de amplitudes de Fourier de los registros en suelo bando y suelo firme. Suponiendo que el espectro de amplitudes de Fourier es una estimado de la función de densidad espectral, el cociente espectral se puede interpretar como una función de transferencia empírica. La técnica está basada en un comportamiento predominante lineal del suelo y en un modelado unidimensional. Una vez que se cuenta con el espectro de amplitudes de Fourier en suelo blando se puede calcular el espectro de respuesta aplicando teoría de vibraciones aleatorias. Alternativamente, los cocientes espectrales pueden estimarse también calculando el cociente entre los espectros de respuesta en un suelo blando y suelo firme. La aproximación que se obtiene en este caso es similar a la del empleo de espectros de Fourier (^{Rosenblueth y Arciniega, 1992}). Los cocientes espectrales H(T) se estimaron con la relación de los espectros de respuesta de las estaciones de la zona de transición (tabla 2) y el espectro correspondiente de la estación CU. Las figuras 6 y 7 muestran los cocientes espectrales en función del periodo para las direcciones N-S y E-W; se muestra también el valor medio de los cocientes espectrales. Se puede notar que las mayores amplificaciones promedio en ambas direcciones horizontales se presentan en periodos alrededor de 1 s, cercanos al periodo del segundo modo de la estructura. Los valores máximos de los cocientes espectrales son del orden de 6 a 10 en las direcciones N-S y E-W, respectivamente. Para efectos de este estudio se consideró el mayor cociente espectral en ambas direcciones para el periodo del segundo modo, el cual corresponde a un valor H4= 5.48 , cuyo valor asociado para el periodo del primer modo es H1.33= 1.47.

El valor medio del logaritmo de la aceleración espectral en el sitio, dados M y R , se calcula entonces como

En el caso de la desviación estándar , σ[ln⁡ SaT]=σ[ln⁡SarT] . Las distribuciones y densidades de probabilidad marginales de las seudo-aceleraciones se pueden calcular entonces como,

en donde, Φ(∙) y ϕ(∙) son las funciones de distribución y de densidad de probabilidad normal estándar, respectivamente, considerando para la aceleración espectral una distribución log-normal.

Empleando las ecuaciones 19 y 20, y el modelo seleccionado de cópula de Gumbel, se pueden calcular las funciones de densidad de probabilidad conjuntas de las aceleraciones espectrales de acuerdo con la ecuación 12. La figura 8 muestra las densidades de probabilidad conjuntas de Sa1 y Sa2 para magnitudes Mw =7 y 8. Se observa una asociación estadística significativa entre las aceleraciones espectrales que se deriva del valor muestral del τ de Kendall. La influencia de la magnitud se evidencia en intervalos de valores notablemente mayores de aceleración espectral para contornos con la misma densidad de probabilidad. Cabe destacar que aun cuando no se ha adoptado un modelo de cópula cuya forma funcional dependa de la magnitud, el argumento ( u, v) de la cópula sí depende de la magnitud a través de las funciones de distribución de las aceleraciones espectrales, tal como se muestra en la ecuación 13.

En el análisis de peligro sísmico, se ha usado ampliamente la ley de ^{Gutenberg y Richter (1944)}, mediante la cual la tasa anual de ocurrencia sismos de magnitud mayores que M en una fuente se puede expresar como,

donde, α y β son los parámetros del modelo. ^{Cornell y Vanmarcke (1969)} propusieron un modelo de tasa de excedencia que toma en cuenta un límite superior Mu para las magnitudes,

En la ecuación 22, λo es la tasa de excedencia de la magnitud Mo . La función de densidad de magnitud queda definida entonces por,

En este estudio se usaron los valores λo= 4.79 (1/año) y β= 2 de acuerdo con el Manual de Obras Civiles, Diseño por Sismo (^{CFE, 2015}).

Para la distorsión máxima de entrepiso γmax , se empleó el modelo de regresión desarrollado por ^{Shome y Cornell (1999)} considerando la contribución de las aceleraciones espectrales (Sa ₁ , Sa ₂ ) asociadas con las frecuencias del modo fundamental y del segundo modo de vibrar del edificio, por lo que la ecuación está dada por,

expresada en función de dichas aceleraciones espectrales de los dos primeros modos, f1=1T1= 0.25 Hz y f2=1T2= 0.75 Hz. La desviación estándar del modelo de regresión es, σln⁡γmax= 0.23. La probabilidad de excedencia de la respuesta γmax en la ecuación 16, dado que Sa1f1=a1 y Sa2f2=a2 , se calcula suponiendo una distribución log-normal, como sigue,

Con fines comparativos, se realizó adicionalmente un análisis de peligro sísmico univariado considerando como medida de intensidad escalar la aceleración espectral asociada con el periodo fundamental de la estructura. Análogamente a la ecuación 16, en este caso la tasa media de excedencia de la distorsión máxima de entrepiso está dada por,

La tasa media de eventos de densidad de aceleración espectral a1 , se calcula de acuerdo con la ecuación 14, como,

El modelo de regresión desarrollado por ^{Shome y Cornell (1999)} para la distorsión máxima de entrepiso γmax del marco de 20 niveles en función de la aceleración espectral Sa1 , está dado por

con una desviación estándar σln⁡γmax= 0.37; la dispersión de este modelo es mayor que la del modelo de regresión expresado en función de las dos aceleraciones espectrales. Se tiene entonces, que

Los resultados del cálculo de las tasas anuales de excedencia para diferentes distorsiones de entrepiso máximas se muestran en la tabla 6, mientras que en la figura 9 se muestran las curvas de tasas de excedencia correspondientes. El cálculo se realizó con un programa desarrollado con el lenguaje R (^{R Core Team, 2020}). En este ejemplo, el empleo de la formulación multivariada produjo estimados significativamente menores. La tasa de excedencia calculada con el análisis univariado es del orden de tres a cuatro veces mayor que la tasa que se obtiene con la formulación multivariada para distorsiones entre 0.01 y 0.02. Para distorsiones mayores, la diferencia entre las tasas se incrementa aún más. Estos resultados reflejan la influencia que la aceleración espectral asociada al segundo modo tiene en la caracterización de la respuesta no lineal en este caso. El origen fundamental de estas diferencias se puede identificar en las ecuaciones 16 y 26 para el cálculo de las tasas de excedencia. La tasa de excedencia en la formulación univariada, se puede escribir como,

Por otro lado, la tasa de excedencia de la respuesta en la formulación multivariada se puede escribir a partir de la ecuación 14 como sigue,

donde

Dado que el producto cSau, v|m,rfSa2a2|m,r en la ecuación 32 representa la densidad de probabilidad de Sa2 dado Sa1, M=m , y R=r , el término P-Sa2γmax>z |a1;m,r se puede interpretar entonces como el valor medio de la probabilidad de excedencia con respecto de Sa2 , condicional a Sa1 , para cada par de valores m , r . Este término permite que en el análisis se incorpore: (1) el efecto del patrón de dependencia estadística entre Sa1 y Sa2 en la caracterización de la respuesta no-lineal; y (2) la variación de dicho patrón de dependencia con las diferentes combinaciones de valores de M y R . En cambio, el término Pγmax>z |a1 en la formulación univariada sólo caracteriza la dependencia de la respuesta con respecto a Sa1 . En el caso de este ejemplo, se muestra que dichas diferencias son significativas. Desde luego, en caso de que la respuesta estructural dependiera ligeramente de la segunda aceleración espectral, los resultados de ambas formulaciones serían similares. Esto ocurre en general cuando el modo fundamental gobierna la respuesta no-lineal de la estructura.

Los resultados del ejemplo de aplicación muestran que considerar una medida de intensidad escalar puede resultar en estimaciones de la tasa de excedencia de la respuesta considerablemente distintas de las que se obtienen cuando la respuesta se caracteriza apropiadamente, según sea el caso, en función de un conjunto de medidas de intensidad. Una mejor estimación de la tasa de excedencia de la respuesta estructural impacta a su vez en una mejor estimación del riesgo sísmico. En consecuencia, el empleo de la formulación multivariada para el análisis de demanda sísmica puede mejorar las decisiones de diseño, evaluación, reparación o reforzamiento, con base en una mejor estimación de pérdidas esperadas adoptando el enfoque del análisis basado en desempeño del Pacific Earthquake Engineering Research Center (PEERC).

Conclusiones

Se ha presentado una formulación multivariada para el análisis probabilístico de peligro sísmico de medidas de intensidad vectoriales empleando cópulas. La formulación es una generalización de la propuesta por ^{Bazurro (1998)} pues expresa la tasa de excedencia en forma genérica en términos de la cópula de una medida de intensidad vectorial, sin suponer que ésta tiene un tipo de distribución de probabilidad conjunta en particular. En lugar de suponer a priori una cierta estructura de dependencia y estimar los parámetros que la caracterizan, por ejemplo, el coeficiente de correlación lineal en el caso gaussiano, este enfoque permite estimar la cópula que represente apropiadamente el patrón de dependencia entre las variables y, en conjunto con las distribuciones marginales, definir el modelo de distribución de probabilidad multivariada. La formulación propuesta se aplicó a un ejemplo de análisis de la tasa de excedencia de la distorsión máxima de entrepiso usando el modelo de regresión de ^{Bazzurro y Cornell (2002)} para un edificio de acero de 20 niveles. Se consideró una medida de intensidad vectorial compuesta por las aceleraciones espectrales para los dos primeros periodos modales de la estructura. El modelo de cópula del vector de aceleraciones espectrales se estimó usando datos de registros sísmicos de la zona de transición de la Ciudad de México. Este tipo de excitación genera una contribución importante de la respuesta asociada con el segundo modo de la estructura. Los parámetros del modelo de distribución de probabilidad marginal de las aceleraciones espectrales se obtuvieron con base en la ley de atenuación espectral de ^{Reyes et al. (2002)} para suelo firme y el empleo de cocientes espectrales. Empleando la formulación multivariada se estimaron valores de la tasa de excedencia de la distorsión de entrepiso máxima que son significativamente menores que los obtenidos mediante un análisis univariado usando como medida de intensidad escalar la aceleración espectral del modo fundamental. La tasa de excedencia calculada con el análisis univariado fue del orden de tres a cuatro veces mayor que la tasa obtenida con la formulación multivariada para distorsiones máximas entre 0.01 y 0.02; para distorsiones mayores, la diferencia entre las tasas se incrementa aún más. Estos resultados ilustran que, tal como ha sido documentado en otros estudios sobre análisis probabilista de demanda sísmica, el uso de una medida de intensidad escalar para el cálculo de la tasa de excedencia de la respuesta estructural no resulta ser apropiado en el caso de estructuras cuya respuesta está correlacionada con varias medidas de intensidad; por ejemplo, estructuras cuyo comportamiento depende de la contribución de modos superiores o de múltiples componentes y parámetros del movimiento sísmico del suelo. El análisis multivariado de peligro sísmico permite considerar la distribución de probabilidad conjunta del vector de medidas de intensidad, a diferencia de una formulación univariada en donde cada medida de intensidad escalar es considerada por separado sin tomar en cuenta su dependencia estadística con las demás.

Dos aspectos del ejemplo de análisis multivariado son relevantes a destacar en comparación con un análisis univariado. En primer lugar, comparando las desviaciones estándar de los modelos de regresión para la distorsión máxima de entrepiso, se observa que la incertidumbre en la predicción de la respuesta es menor al considerar las aceleraciones espectrales asociadas a los dos primeros modos que cuando ésta se expresa sólo en términos de la aceleración espectral para el modo fundamental (^{Shome y Cornell, 1999}; ^{Bazzurro y Cornell, 2002}). De acuerdo con la definición de eficiencia de una medida de intensidad planteada por ^{Luco y Cornell (2007)}, el empleo del modelo bivariado es más eficiente en el sentido de que una menor cantidad de registros y análisis estructurales serían necesarios para caracterizar con una adecuada precisión la respuesta en comparación con el caso univariado. En segundo lugar, el análisis bivariado incorpora explícitamente el efecto del patrón de dependencia entre las aceleraciones espectrales en el análisis de demanda sísmica, empleando el modelo de regresión de Bazzurro y Cornell (1999) que se basa en análisis no lineales, así como la variación de dicho patrón de dependencia con la magnitud M y la distancia R del sitio a la fuente. Incorporar dichos efectos a través de un enfoque univariado puede resultar en un procedimiento complejo e impráctico. Una opción sería realizar el análisis probabilista de demanda sísmica empleando modelos de medidas de intensidad escalares expresados como combinaciones de un conjunto de medidas de intensidad que caractericen la respuesta en diferentes periodos modales (^{Luco y Cornell, 2007}). La construcción de dichos modelos es un intrincado procedimiento que requiere caracterizar apropiadamente no sólo la contribución de modos superiores a la respuesta no lineal de interés, sino también la relación entre los estimados de la respuesta considerando varios modos de vibrar y considerando solo el modo fundamental. La formulación que aquí se ha presentado establece un procedimiento explícito y riguroso para un análisis multivariado, evitando el empleo de procedimientos que intentan caracterizar formas espectrales e incorporar la contribución de los modos superiores en función de la respuesta del modo fundamental.

El cálculo de la tasa de excedencia de las medidas de intensidad en el análisis de peligro sísmico y de la tasa de excedencia de un parámetro de respuesta en el análisis probabilista de demanda sísmica, son dos de los pasos fundamentales para el análisis de riesgo sísmico de acuerdo con el marco metodológico del análisis sísmico basado en desempeño propuesto por el Pacific Earthquake Engineering Research Center (PEERC). Los resultados del ejemplo de aplicación muestran que, en ciertos casos, el empleo de una formulación univariada considerando una medida de intensidad escalar puede resultar en estimaciones de la tasa de excedencia de la respuesta considerablemente distintas de las que se obtienen cuando la respuesta se caracteriza de manera más eficiente mediante un conjunto de medidas de intensidad. En tales casos, las decisiones de diseño, evaluación, reparación o reforzamiento, dirigidas a mitigar el riesgo sísmico y alcanzar niveles apropiados de confiabilidad, debieran estar basada en los resultados de un enfoque de análisis multivariado. Algunos tópicos que no han formado parte del alcance de este trabajo pueden ser abordados en estudios futuros. Trabajos futuros pueden enfocarse a la caracterización multivariada de peligro sísmico para el análisis probabilista de demanda sísmica de edificios altos de concreto reforzado y de obras de infraestructura como puentes, taludes y plataformas marinas, al estudio de la influencia del componente vertical del movimiento sísmico en la tasa de excedencia de la respuesta, y al análisis de construcciones sometidas a efectos direccionales de movimientos sísmicos de fuentes cercanas. La estimación de los parámetros de las cópulas depende, como en cualquier problema de inferencia estadística, del tamaño de la muestra de datos. Dado que se pueden calcular intervalos de confianza de los parámetros estimados, en estudios futuros se puede evaluar el efecto de la incertidumbre de la estimación estadística de las cópulas en la evaluación de la tasa de excedencia de la respuesta.