# -*- 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
""")