Source code for openquake.hazardlib.gsim.campbell_bozorgnia_2008

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2013-2016 GEM Foundation
#
# OpenQuake is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# OpenQuake is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>.

"""
Module exports :class:`CampbellBozorgnia2008`, and 
:class:'CampbellBozorgnia2008Arbitrary'
"""
from __future__ import division

import numpy as np
from math import log, exp
from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, PGV, PGD, CAV, SA


[docs]class CampbellBozorgnia2008(GMPE): """ Implements GMPE developed by Kenneth W. Campbell and Yousef Bozorgnia, published as "NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5 % Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10s" (2008, Earthquake Spectra, Volume 24, Number 1, pages 139 - 171). This class implements the model for the Geometric Mean of the elastic spectra. Included in the coefficient set are the coefficients for the Campbell & Bozorgnia (2010) GMPE for predicting Cumulative Absolute Velocity (CAV), published as "A Ground Motion Prediction Equation for the Horizontal Component of Cumulative Absolute Velocity (CSV) Based on the PEER-NGA Strong Motion Database" (2010, Earthquake Spectra, Volume 26, Number 3, 635 - 650). """ #: Supported tectonic region type is active shallow crust DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST #: Supported intensity measure types are spectral acceleration, peak #: ground velocity, peak ground displacement and peak ground acceleration #: Additional model for cumulative absolute velocity defined in #: Campbell & Bozorgnia (2010) DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([ PGA, PGV, PGD, CAV, SA ]) #: Supported intensity measure component is orientation-independent #: average horizontal :attr:`~openquake.hazardlib.const.IMC.GMRotI50` DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.GMRotI50 #: Supported standard deviation types are inter-event, intra-event #: and total, see section "Aleatory Uncertainty Model", page 147. DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([ const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT ]) #: Required site parameters are Vs30, Vs30 type (measured or inferred), #: and depth (km) to the 2.5 km/s shear wave velocity layer (z2pt5) REQUIRES_SITES_PARAMETERS = set(('vs30', 'z2pt5')) #: Required rupture parameters are magnitude, rake, dip, ztor REQUIRES_RUPTURE_PARAMETERS = set(('mag', 'rake', 'dip', 'ztor')) #: Required distance measures are Rrup and Rjb. REQUIRES_DISTANCES = set(('rrup', 'rjb'))
[docs] def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types): """ See :meth:`superclass method <.base.GroundShakingIntensityModel.get_mean_and_stddevs>` for spec of input and result values. """ # extract dictionaries of coefficients specific to required # intensity measure type and for PGA C = self.COEFFS[imt] C_PGA = self.COEFFS[PGA()] # compute median pga on rock (vs30=1100), needed for site response # term calculation # For spectral accelerations at periods between 0.0 and 0.25 s, Sa (T) # cannot be less than PGA on soil, therefore if the IMT is in this # period range it is necessary to calculate PGA on soil if isinstance(imt, SA) and (imt.period > 0.0) and (imt.period < 0.25): get_pga_site = True else: get_pga_site = False pga1100, pga_site = self._compute_imt1100(C_PGA, sites, rup, dists, get_pga_site) # Get the median ground motion mean = (self._compute_magnitude_term(C, rup.mag) + self._compute_distance_term(C, rup, dists) + self._compute_style_of_faulting_term(C, rup) + self._compute_hanging_wall_term(C, rup, dists) + self._compute_shallow_site_response(C, sites, pga1100) + self._compute_basin_response_term(C, sites.z2pt5)) # If it is necessary to ensure that Sa(T) >= PGA (see previous comment) if get_pga_site: idx = mean < np.log(pga_site) mean[idx] = np.log(pga_site[idx]) stddevs = self._get_stddevs(C, sites, pga1100, C_PGA['s_lny'], stddev_types) return mean, stddevs
def _compute_imt1100(self, C, sites, rup, dists, get_pga_site=False): """ Computes the PGA on reference (Vs30 = 1100 m/s) rock. """ # Calculates simple site response term assuming all sites 1100 m/s fsite = (C['c10'] + (C['k2'] * C['n'])) * log(1100. / C['k1']) # Calculates the PGA on rock pga1100 = np.exp(self._compute_magnitude_term(C, rup.mag) + self._compute_distance_term(C, rup, dists) + self._compute_style_of_faulting_term(C, rup) + self._compute_hanging_wall_term(C, rup, dists) + self._compute_basin_response_term(C, sites.z2pt5) + fsite) # If PGA at the site is needed then remove factor for rock and # re-calculate on correct site condition if get_pga_site: pga_site = np.exp(np.log(pga1100) - fsite) fsite = self._compute_shallow_site_response(C, sites, pga1100) pga_site = np.exp(np.log(pga_site) + fsite) else: pga_site = None return pga1100, pga_site def _compute_magnitude_term(self, C, mag): """ Returns the magnitude scaling factor (equation (2), page 144) """ fmag = C['c0'] + C['c1'] * mag if mag <= 5.5: return fmag elif mag > 6.5: return fmag + (C['c2'] * (mag - 5.5)) + (C['c3'] * (mag - 6.5)) else: return fmag + (C['c2'] * (mag - 5.5)) def _compute_distance_term(self, C, rup, dists): """ Returns the distance scaling factor (equation (3), page 145) """ return (C['c4'] + C['c5'] * rup.mag) * \ np.log(np.sqrt(dists.rrup ** 2. + C['c6'] ** 2.)) def _compute_style_of_faulting_term(self, C, rup): """ Returns the style of faulting factor, depending on the mechanism (rake) and top of rupture depth (equations (4) and (5), pages 145 - 146) """ frv, fnm = self._get_fault_type_dummy_variables(rup.rake) if frv > 0.: # Top of rupture depth term only applies to reverse faults if rup.ztor < 1.: ffltz = rup.ztor else: ffltz = 1. else: ffltz = 0. return (C['c7'] * frv * ffltz) + (C['c8'] * fnm) def _get_fault_type_dummy_variables(self, rake): """ Returns the coefficients FRV and FNM, describing if the rupture is reverse (FRV = 1.0, FNM = 0.0), normal (FRV = 0.0, FNM = 1.0) or strike-slip/oblique-slip (FRV = 0.0, FNM = 0.0). Reverse faults are classified as those with a rake in the range 30 to 150 degrees. Normal faults are classified as having a rake in the range -150 to -30 degrees :returns: FRV, FNM """ if (rake > 30.0) and (rake < 150.): return 1., 0. elif (rake > -150.0) and (rake < -30.0): return 0., 1. else: return 0., 0. def _compute_hanging_wall_term(self, C, rup, dists): """ Returns the hanging wall scaling term, the product of the scaling coefficient and four separate scaling terms for distance, magnitude, rupture depth and dip (equations 6 - 10, page 146). Individual scaling terms defined in separate functions """ return (C['c9'] * self._get_hanging_wall_distance_term(dists, rup.ztor) * self._get_hanging_wall_magnitude_term(rup.mag) * self._get_hanging_wall_depth_term(rup.ztor) * self._get_hanging_wall_dip_term(rup.dip)) def _get_hanging_wall_distance_term(self, dists, ztor): """ Returns the hanging wall distance scaling term (equation 7, page 146) """ fhngr = np.ones_like(dists.rjb, dtype=float) idx = dists.rjb > 0. if ztor < 1.: temp_rjb = np.sqrt(dists.rjb[idx] ** 2. + 1.) r_max = np.max(np.column_stack([dists.rrup[idx], temp_rjb]), axis=1) fhngr[idx] = (r_max - dists.rjb[idx]) / r_max else: fhngr[idx] = (dists.rrup[idx] - dists.rjb[idx]) / dists.rrup[idx] return fhngr def _get_hanging_wall_magnitude_term(self, mag): """ Returns the hanging wall magnitude scaling term (equation 8, page 146) """ if mag <= 6.0: return 0. elif mag >= 6.5: return 1. else: return 2. * (mag - 6.0) def _get_hanging_wall_depth_term(self, ztor): """ Returns the hanging wall depth scaling term (equation 9, page 146) """ if ztor >= 20.0: return 0. else: return (20. - ztor) / 20.0 def _get_hanging_wall_dip_term(self, dip): """ Returns the hanging wall dip scaling term (equation 10, page 146) """ if dip > 70.0: return (90.0 - dip) / 20.0 else: return 1.0 def _compute_shallow_site_response(self, C, sites, pga1100): """ Returns the shallow site response term (equation 11, page 146) """ stiff_factor = C['c10'] + (C['k2'] * C['n']) # Initially default all sites to intermediate rock value fsite = stiff_factor * np.log(sites.vs30 / C['k1']) # Check for soft soil sites idx = sites.vs30 < C['k1'] if np.any(idx): pga_scale = np.log(pga1100[idx] + (C['c'] * ((sites.vs30[idx] / C['k1']) ** C['n']))) - np.log(pga1100[idx] + C['c']) fsite[idx] = C['c10'] * np.log(sites.vs30[idx] / C['k1']) + \ (C['k2'] * pga_scale) # Any very hard rock sites are rendered to the constant amplification # factor idx = sites.vs30 >= 1100. if np.any(idx): fsite[idx] = stiff_factor * log(1100. / C['k1']) return fsite def _compute_basin_response_term(self, C, z2pt5): """ Returns the basin response term (equation 12, page 146) """ fsed = np.zeros_like(z2pt5, dtype=float) idx = z2pt5 < 1.0 if np.any(idx): fsed[idx] = C['c11'] * (z2pt5[idx] - 1.0) idx = z2pt5 > 3.0 if np.any(idx): fsed[idx] = (C['c12'] * C['k3'] * exp(-0.75)) *\ (1.0 - np.exp(-0.25 * (z2pt5[idx] - 3.0))) return fsed def _get_stddevs(self, C, sites, pga1100, sigma_pga, stddev_types): """ Returns the standard deviations as described in the "ALEATORY UNCERTAINTY MODEL" section of the paper. Equations 13 to 19, pages 147 to 151 """ std_intra = self._compute_intra_event_std(C, sites.vs30, pga1100, sigma_pga) std_inter = C['t_lny'] * np.ones_like(sites.vs30) stddevs = [] for stddev_type in stddev_types: assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES if stddev_type == const.StdDev.TOTAL: stddevs.append(self._get_total_sigma(C, std_intra, std_inter)) elif stddev_type == const.StdDev.INTRA_EVENT: stddevs.append(std_intra) elif stddev_type == const.StdDev.INTER_EVENT: stddevs.append(std_inter) return stddevs def _compute_intra_event_std(self, C, vs30, pga1100, sigma_pga): """ Returns the intra-event standard deviation at the site, as defined in equation 15, page 147 """ # Get intra-event standard deviation at the base of the site profile sig_lnyb = np.sqrt(C['s_lny'] ** 2. - C['s_lnAF'] ** 2.) sig_lnab = np.sqrt(sigma_pga ** 2. - C['s_lnAF'] ** 2.) # Get linearised relationship between f_site and ln PGA alpha = self._compute_intra_event_alpha(C, vs30, pga1100) return np.sqrt( (sig_lnyb ** 2.) + (C['s_lnAF'] ** 2.) + ((alpha ** 2.) * (sig_lnab ** 2.)) + (2.0 * alpha * C['rho'] * sig_lnyb * sig_lnab)) def _compute_intra_event_alpha(self, C, vs30, pga1100): """ Returns the linearised functional relationship between fsite and pga1100, determined from the partial derivative defined on equation 17 on page 148 """ alpha = np.zeros_like(vs30, dtype=float) idx = vs30 < C['k1'] if np.any(idx): temp1 = (pga1100[idx] + C['c'] * (vs30[idx] / C['k1']) ** C['n']) ** -1. temp1 = temp1 - ((pga1100[idx] + C['c']) ** -1.) alpha[idx] = C['k2'] * pga1100[idx] * temp1 return alpha def _get_total_sigma(self, C, std_intra, std_inter): """ Returns the total sigma term as defined by equation 16, page 147 This method is defined here as the Campbell & Bozorgnia (2008) model can also be applied to the "arbitrary" horizontal component definition, in which case the total sigma is modified. """ return np.sqrt(std_intra ** 2. + std_inter ** 2.) COEFFS = CoeffsTable(sa_damping=5, table="""\ imt c0 c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 k1 k2 k3 c n s_lny t_lny s_lnAF c_lny rho cav -4.354 0.942 -0.178 -0.346 -1.309 0.087 7.24 0.111 -0.108 0.362 2.549 0.090 1.277 400 -2.690 1.000 1.88 1.18 0.371 0.196 0.300 0.089 0.735 pgd -5.270 1.600 -0.070 0.000 -2.000 0.17 4.00 0.000 0.000 0.000 -0.820 0.300 1.000 400 0.000 2.744 1.88 1.18 0.667 0.485 0.300 0.290 0.174 pgv 0.954 0.696 -0.309 -0.019 -2.016 0.17 4.00 0.245 0.000 0.358 1.694 0.092 1.000 400 -1.955 1.929 1.88 1.18 0.484 0.203 0.300 0.190 0.691 pga -1.715 0.500 -0.530 -0.262 -2.118 0.17 5.60 0.280 -0.120 0.490 1.058 0.040 0.610 865 -1.186 1.839 1.88 1.18 0.478 0.219 0.300 0.166 1.000 0.010 -1.715 0.500 -0.530 -0.262 -2.118 0.17 5.60 0.280 -0.120 0.490 1.058 0.040 0.610 865 -1.186 1.839 1.88 1.18 0.478 0.219 0.300 0.166 1.000 0.020 -1.680 0.500 -0.530 -0.262 -2.123 0.17 5.60 0.280 -0.120 0.490 1.102 0.040 0.610 865 -1.219 1.840 1.88 1.18 0.480 0.219 0.300 0.166 0.999 0.030 -1.552 0.500 -0.530 -0.262 -2.145 0.17 5.60 0.280 -0.120 0.490 1.174 0.040 0.610 908 -1.273 1.841 1.88 1.18 0.489 0.235 0.300 0.165 0.989 0.050 -1.209 0.500 -0.530 -0.267 -2.199 0.17 5.74 0.280 -0.120 0.490 1.272 0.040 0.610 1054 -1.346 1.843 1.88 1.18 0.510 0.258 0.300 0.162 0.963 0.075 -0.657 0.500 -0.530 -0.302 -2.277 0.17 7.09 0.280 -0.120 0.490 1.438 0.040 0.610 1086 -1.471 1.845 1.88 1.18 0.520 0.292 0.300 0.158 0.922 0.100 -0.314 0.500 -0.530 -0.324 -2.318 0.17 8.05 0.280 -0.099 0.490 1.604 0.040 0.610 1032 -1.624 1.847 1.88 1.18 0.531 0.286 0.300 0.170 0.898 0.150 -0.133 0.500 -0.530 -0.339 -2.309 0.17 8.79 0.280 -0.048 0.490 1.928 0.040 0.610 878 -1.931 1.852 1.88 1.18 0.532 0.280 0.300 0.180 0.890 0.200 -0.486 0.500 -0.446 -0.398 -2.220 0.17 7.60 0.280 -0.012 0.490 2.194 0.040 0.610 748 -2.188 1.856 1.88 1.18 0.534 0.249 0.300 0.186 0.871 0.250 -0.890 0.500 -0.362 -0.458 -2.146 0.17 6.58 0.280 0.000 0.490 2.351 0.040 0.700 654 -2.381 1.861 1.88 1.18 0.534 0.240 0.300 0.191 0.852 0.300 -1.171 0.500 -0.294 -0.511 -2.095 0.17 6.04 0.280 0.000 0.490 2.460 0.040 0.750 587 -2.518 1.865 1.88 1.18 0.544 0.215 0.300 0.198 0.831 0.400 -1.466 0.500 -0.186 -0.592 -2.066 0.17 5.30 0.280 0.000 0.490 2.587 0.040 0.850 503 -2.657 1.874 1.88 1.18 0.541 0.217 0.300 0.206 0.785 0.500 -2.569 0.656 -0.304 -0.536 -2.041 0.17 4.73 0.280 0.000 0.490 2.544 0.040 0.883 457 -2.669 1.883 1.88 1.18 0.550 0.214 0.300 0.208 0.735 0.750 -4.844 0.972 -0.578 -0.406 -2.000 0.17 4.00 0.280 0.000 0.490 2.133 0.077 1.000 410 -2.401 1.906 1.88 1.18 0.568 0.227 0.300 0.221 0.628 1.000 -6.406 1.196 -0.772 -0.314 -2.000 0.17 4.00 0.255 0.000 0.490 1.571 0.150 1.000 400 -1.955 1.929 1.88 1.18 0.568 0.255 0.300 0.225 0.534 1.500 -8.692 1.513 -1.046 -0.185 -2.000 0.17 4.00 0.161 0.000 0.490 0.406 0.253 1.000 400 -1.025 1.974 1.88 1.18 0.564 0.296 0.300 0.222 0.411 2.000 -9.701 1.600 -0.978 -0.236 -2.000 0.17 4.00 0.094 0.000 0.371 -0.456 0.300 1.000 400 -0.299 2.019 1.88 1.18 0.571 0.296 0.300 0.226 0.331 3.000 -10.556 1.600 -0.638 -0.491 -2.000 0.17 4.00 0.000 0.000 0.154 -0.820 0.300 1.000 400 0.000 2.110 1.88 1.18 0.558 0.326 0.300 0.229 0.289 4.000 -11.212 1.600 -0.316 -0.770 -2.000 0.17 4.00 0.000 0.000 0.000 -0.820 0.300 1.000 400 0.000 2.200 1.88 1.18 0.576 0.297 0.300 0.237 0.261 5.000 -11.684 1.600 -0.070 -0.986 -2.000 0.17 4.00 0.000 0.000 0.000 -0.820 0.300 1.000 400 0.000 2.291 1.88 1.18 0.601 0.359 0.300 0.237 0.200 7.500 -12.505 1.600 -0.070 -0.656 -2.000 0.17 4.00 0.000 0.000 0.000 -0.820 0.300 1.000 400 0.000 2.517 1.88 1.18 0.628 0.428 0.300 0.271 0.174 10.00 -13.087 1.600 -0.070 -0.422 -2.000 0.17 4.00 0.000 0.000 0.000 -0.820 0.300 1.000 400 0.000 2.744 1.88 1.18 0.667 0.485 0.300 0.290 0.174 """)
[docs]class CampbellBozorgnia2008Arbitrary(CampbellBozorgnia2008): """ Implements the Campbell & Bozorgnia (2008) GMPE as modified to represent the arbitrary horizontal component of ground motion, instead of the Rotationally Independent Geometric Mean (GMRotI) originally defined in the paper. """ #: Supported intensity measure component is arbitrary horizontal #: :attr:`~openquake.hazardlib.const.IMC.HORIZONTAL`, DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.HORIZONTAL def _get_total_sigma(self, C, std_intra, std_inter): """ Returns the total sigma term for the arbitrary horizontal component of ground motion defined by equation 18, page 150 """ return np.sqrt(std_intra ** 2. + std_inter ** 2. + C['c_lny'] ** 2.)