Source code for openquake.hazardlib.gsim.shahjouei_pezeshk_2016

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2013-2017 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:'ShahjoueiPezeshk2016'.
"""
from __future__ import division

import numpy as np

from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, PGV, SA


[docs]class ShahjoueiPezeshk2016(GMPE): """ Implements GMPE developed by Alireza Shahjouei and Shahram Pezeshk. Published as "Alternative Hybrid Empirical Ground‐Motion Model for Central and Eastern North America Using Hybrid Simulations and NGA‐West2 Models", 2016, Bulletin of the Seismological Society of America, vol. 106, no. 2, 734 - 754. """ #: Supported tectonic region type is 'stable continental region' #: equation has been derived from data from Eastern North America (ENA) # 'Introduction', page 735. DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.STABLE_CONTINENTAL #: Supported intensity measure types are spectral acceleration, #: and peak ground acceleration. See Table 7 on page 743 DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([ PGA, PGV, SA ]) #: An orientation-independent alternative to :attr:`AVERAGE_HORIZONTAL`. #: Defined at Boore et al. (2006, Bull. Seism. Soc. Am. 96, 1502-1511) #: and is used for all the NGA GMPEs. See page 742. #: :attr:'~openquake.hazardlib.const.IMC.RotD50'. DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.RotD50 #: Supported standard deviation types is total. #: See equation 4 and 5, page 744. #: We use aleatory uncertainty as the total #: standard deviation since page 745 states, #: The epistemic uncertainty for an individual GMM is infrequently employed... DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([ const.StdDev.TOTAL ]) #: No site parameters are needed. The GMPE was developed for hard-rock site # with Vs30 >= 3000 m/s (NEHRP site class A) only. Page 734. REQUIRES_SITES_PARAMETERS = set() #: Required rupture parameters are magnitude (eq. 4, page 742). REQUIRES_RUPTURE_PARAMETERS = set(('mag', )) #: Required distance measure is Rjb (eq. 3 page 742). REQUIRES_DISTANCES = set(('rjb', )) #: GMPE not tested against independent implementation so raise #: not verified warning non_verified = True
[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. """ # Extracting dictionary of coefficients specific to required # intensity measure type. C = self.COEFFS[imt] imean = (self._compute_magnitude(rup, C) + self._compute_attenuation(rup, dists, imt, C) + self._compute_distance(rup, dists, imt, C)) mean = np.log(10.0 ** (imean)) istddevs = self._get_stddevs(C, stddev_types, rup, imt, num_sites=len(dists.rjb)) stddevs = np.log(10.0 ** np.array(istddevs)) return mean, stddevs
def _get_stddevs(self, C, stddev_types, rup, imt, num_sites): """ Return standard deviations as defined in eq. 4 and 5, page 744, based on table 8, page 744. Eq. 5 yields std dev in natural log, so convert to log10 """ stddevs = [] for stddev_type in stddev_types: sigma_mean = self._compute_standard_dev(rup, imt, C) sigma_tot = np.sqrt((sigma_mean ** 2) + (C['SigmaReg'] ** 2)) sigma_tot = np.log10(np.exp(sigma_tot)) stddevs.append(sigma_tot + np.zeros(num_sites)) return stddevs def _compute_magnitude(self, rup, C): """ Compute the first term of the equation described on p. 742: "c1 + (c2 * M) + (c3 * M**2) " """ return C['c1'] + (C['c2'] * rup.mag) + (C['c3'] * (rup.mag ** 2)) def _compute_attenuation(self, rup, dists, imt, C): """ Compute the second term of the equation described on p. 742: " [(c4 + c5 * M) * min{ log10(R), log10(60.) }] + [(c6 + c7 * M) * max{ min{ log10(R/60.), log10(120./60.) }, 0.}] + [(c8 + c9 * M) * max{ log10(R/120.), 0}] " """ vec = np.ones(len(dists.rjb)) a1 = (np.log10(np.sqrt(dists.rjb ** 2.0 + C['c11'] ** 2.0)), np.log10(60. * vec)) a = np.column_stack([a1[0], a1[1]]) b3 = (np.log10(np.sqrt(dists.rjb ** 2.0 + C['c11'] ** 2.0) / (60. * vec)), np.log10((120. / 60.) * vec)) b2 = np.column_stack([b3[0], b3[1]]) b1 = ([np.min(b2, axis=1), 0. * vec]) b = np.column_stack([b1[0], b1[1]]) c1 = (np.log10(np.sqrt(dists.rjb ** 2.0 + C['c11'] ** 2.0) / (120.) * vec), 0. * vec) c = np.column_stack([c1[0], c1[1]]) return (((C['c4'] + C['c5'] * rup.mag) * np.min(a, axis=1)) + ((C['c6'] + C['c7'] * rup.mag) * np.max(b, axis=1)) + ((C['c8'] + C['c9'] * rup.mag) * np.max(c, axis=1))) def _compute_distance(self, rup, dists, imt, C): """ Compute the third term of the equation described on p. 742: " c10 * R " """ return (C['c10'] * np.sqrt(dists.rjb ** 2.0 + C['c11'] ** 2.0)) def _compute_standard_dev(self, rup, imt, C): """ Compute the the standard deviation in terms of magnitude described on page 744, eq. 4 """ sigma_mean = 0. if isinstance(imt, (PGA, SA)): psi = -6.898E-3 else: psi = -3.054E-5 if rup.mag <= 6.5: sigma_mean = (C['c12'] * rup.mag) + C['c13'] elif rup.mag > 6.5: sigma_mean = (psi * rup.mag) + C['c14'] return sigma_mean #: Equation coefficients, described in Table 2 on pp. 1865 COEFFS = CoeffsTable(sa_damping=5, table=""" IMT c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 SigmaReg SigmaPar pgv -2.3891 1.259 -0.07901 -2.9386 0.3034 -0.00929 -0.04605 -2.7548 0.3467 -0.0007623 -4.598 -0.0554 0.978 0.663 0.1 0.288 pga -0.3002 0.5066 -0.04526 -3.224 0.2998 -1.283 0.1045 -3.0856 0.2778 -0.0007711 3.81 -0.041 0.876 0.611 0.194 0.373 0.010 -0.3472 0.4838 -0.04093 -3.0832 0.2712 -0.9676 0.04983 -2.9695 0.2693 -0.0006695 -4.434 -0.056 0.982 0.664 0.132 0.281 0.020 0.832 0.1934 -0.0206 -3.1134 0.2786 -1.133 0.05994 -3.5023 0.2901 -0.0005857 -4.412 -0.0559 0.983 0.665 0.0928 0.281 0.030 1.185 0.1064 -0.01423 -3.1029 0.2792 -1.078 0.05239 -3.5722 0.2865 -0.000622 -4.353 -0.0577 1 0.676 0.0833 0.277 0.040 1.246 0.08986 -0.01268 -3.0785 0.2773 -0.9743 0.0416 -3.5083 0.2769 -0.0006818 -4.303 -0.0577 1.01 0.688 0.0798 0.279 0.050 1.1793 0.1037 -0.01321 -3.0488 0.2744 -0.8635 0.03077 -3.3986 0.2659 -0.0007439 -4.266 -0.0578 1.03 0.701 0.0776 0.272 0.075 0.8045 0.1866 -0.01788 -2.9697 0.266 -0.6122 0.007491 -3.0852 0.2391 -0.0008801 -4.214 -0.0561 1.03 0.721 0.0738 0.252 0.100 0.35 0.2871 -0.02381 -2.894 0.2576 -0.4123 -0.01012 -2.7947 0.2163 -0.0009848 4.201 -0.0565 1.05 0.732 0.0717 0.265 0.150 -0.5264 0.4782 -0.03519 -2.761 0.2426 -0.1319 -0.03338 -2.3312 0.1818 -0.001125 4.239 -0.0559 1.04 0.724 0.0716 0.276 0.200 -1.2884 0.6413 -0.04486 -2.6504 0.2301 0.04637 -0.0469 -1.9927 0.1576 -0.001209 4.325 -0.056 1.03 0.715 0.0743 0.258 0.250 -1.9422 0.7789 -0.05295 -2.5573 0.2196 0.1631 -0.05478 -1.7399 0.1398 -0.001258 4.438 -0.0537 1.02 0.712 0.0779 0.268 0.300 -2.5071 0.8961 -0.05976 -2.478 0.2107 0.2407 -0.05919 -1.547 0.1265 -0.001286 4.571 -0.0511 1.01 0.718 0.0815 0.284 0.400 -3.436 1.085 -0.07059 -2.3495 0.1961 0.3244 -0.06197 -1.2793 0.1085 -0.001304 -4.872 -0.047 0.987 0.725 0.0876 0.34 0.500 -4.1699 1.231 -0.07878 -2.251 0.1849 0.3544 -0.06046 -1.1111 0.09757 -0.001294 -5.211 -0.0442 0.981 0.736 0.0923 0.357 0.750 -5.4797 1.482 -0.09245 -2.0865 0.1659 0.3284 -0.04979 -0.9131 0.0857 -0.001219 -6.154 -0.0384 0.967 0.76 0.0991 0.374 1.000 -6.3464 1.641 -0.1006 -1.9931 0.1546 0.253 -0.03709 -0.8641 0.08405 -0.001123 -7.174 -0.0314 0.933 0.77 0.102 0.392 1.500 -7.4087 1.823 -0.1093 -1.9162 0.1438 0.09019 -0.01551 -0.92 0.09103 -0.0009407 -9.253 -0.0227 0.883 0.776 0.105 0.426 2.000 -8.0057 1.916 -0.113 -1.9173 0.1418 -0.03828 -0.001252 -1.0327 0.1016 -0.0007926 -11.22 -0.0184 0.857 0.778 0.106 0.44 3.000 -8.5793 1.985 -0.1146 -2.0184 0.1499 -0.1744 0.009393 -1.2453 0.1214 -0.0005919 14.38 -0.0189 0.859 0.777 0.107 0.58 4.000 -8.8246 1.99 -0.1131 -2.1475 0.1635 -0.1844 0.003919 -1.3849 0.1357 -0.0004855 16.19 -0.016 0.83 0.766 0.107 0.589 5.000 -8.9855 1.975 -0.1105 -2.2496 0.1764 -0.1043 -0.01187 -1.4511 0.1446 -0.0004439 16.71 -0.0153 0.826 0.766 0.107 0.631 7.500 -9.3927 1.925 -0.1032 -2.3572 0.1973 0.3465 -0.07832 -1.3728 0.149 -0.0005176 14.58 -0.0143 0.815 0.762 0.113 0.721 10.000 -9.735 1.879 -0.09666 -2.4139 0.2117 1.01 -0.1678 -1.0631 0.137 -0.000742 11.23 -0.017 0.822 0.752 0.14 0.739 """)