Source code for openquake.hazardlib.gsim.campbell_1997

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2013-2019 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:`Campbell1997`
"""
import numpy as np

from openquake.hazardlib.gsim.base import GMPE
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA


[docs]class Campbell1997(GMPE): """ Implements GMPE (PGA) by Campbell, Kenneth W. "Empirical near-source attenuation relationships for horizontal and vertical components of peak ground acceleration, peak ground velocity, and pseudo-absolute acceleration response spectra." Seismological research letters 68.1 (1997): 154-179. """ #: Supported TRT active...we specify active_shallow_crust DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST #: Supported intensity measure types are PGA, PGV, PSA, but we only define #: PGA because this is the only IMT used by an implemented model (09/18) DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([ PGA ]) #: Supported intensity measure component is the horizontal component DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.AVERAGE_HORIZONTAL #: Supported standard deviation type is only total, see equation 4, pg 164 DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([ const.StdDev.TOTAL, ]) #: Requires vs30 REQUIRES_SITES_PARAMETERS = set(('vs30',)) #: Required rupture parameters are magnitude and top of rupture depth REQUIRES_RUPTURE_PARAMETERS = set(('mag', 'rake')) #: Required distance measure is closest distance to rupture. In the #: publication, Rseis is used. We assume Rrup=Rseis, justified by #: our calculations matching the verification tables REQUIRES_DISTANCES = set(('rrup', )) #: Verification of the mean value was done by digitizing Figs. 9 and 10 #: using Engauge Digitizer. The tests check varied magnitude, distance, #: vs30, and faulting type. Maximum error was ~1.3%. OpenQuake trellis #: plots match these figures. Also tested against a matlab implementation #: (web.stanford.edu/~bakerjw/GMPEs/C_1997_horiz.m), which also has no #: verification tables.
[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. """ R = (dists.rrup) M = rup.mag # get constants Ssr = self.get_Ssr_term(sites.vs30) Shr = self.get_Shr_term(sites.vs30) rake = rup.rake F = self.get_fault_term(rake) # compute mean mean = -3.512 + (0.904 * M) - (1.328 * np.log(np.sqrt(R**2 + (0.149 * np.exp(0.647 * M))**2))) \ + (1.125 - 0.112 * np.log(R) - 0.0957 * M) * F \ + (0.440 - 0.171 * np.log(R)) * Ssr \ + (0.405 - 0.222 * np.log(R)) * Shr stddevs = self.get_stddevs(mean, stddev_types) return mean, stddevs
[docs] def get_fault_term(self, rake): """ Returns coefficient for faulting style (pg 156) """ rake = rake + 360 if rake < 0 else rake if (rake >= 45) & (rake <= 135): f = 1. elif (rake >= 225) & (rake <= 315): f = 0.5 else: f = 0. return f
[docs] def get_Ssr_term(self, vs30): """ Returns site term for soft rock (pg 157) """ return (vs30 >= 760) & (vs30 < 1500)
[docs] def get_Shr_term(self, vs30): """ Returns site term for hard rock (pg 157) """ return vs30 >= 1500
[docs] def get_stddevs(self, mean, stddev_types): """ Returns the standard deviations from mean (pg 164; more robust than estimate using magnitude) """ mean = np.exp(mean) sigma = 0.39 + np.zeros(mean.shape) sigma[mean < 0.068] = 0.55 idx = np.logical_and(mean >= 0.068, mean <= 0.21) sigma[idx] = 0.173- 0.140 * np.log(mean[idx]) stddevs = [] for stddev in stddev_types: if stddev == const.StdDev.TOTAL: stddevs.append(sigma) return stddevs