# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (c) 2013-2018 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:`HongGoda2007RotD100`.
"""
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 HongGoda2007(GMPE):
"""
Implements GMPE developed for RotD100 ground motion as defined by
Hong, H. P. and Goda, K. (2007), "Orientation-Dependent Ground Motion
Measure for Seismic Hazard Assessment", Bull. Seism. Soc. Am. 97(5),
1525 - 1538
This is really an experimental GMPE in which the amplification term
is taken directly from Atkinson & Boore (2006) rather than constrained
by the records themselves. There may exist a possible units issue as
the amplification function for AB2006 is in cm/s/s whereas the
GMPE here is given in g
"""
#: The supported tectonic region type is active shallow crust
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST
#: The supported intensity measure types are PGA, PGV, and SA
DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([
PGA,
PGV,
SA
])
#: The supported intensity measure component is RotD100
DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.RotD100
#: The supported standard deviations are total, inter and intra event, see
#: table 4.a, pages 22-23
DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([
const.StdDev.TOTAL,
const.StdDev.INTER_EVENT,
const.StdDev.INTRA_EVENT
])
#: The required site parameter is vs30, see equation 1, page 20.
REQUIRES_SITES_PARAMETERS = set(('vs30',))
#: The required rupture parameters are magnitude
REQUIRES_RUPTURE_PARAMETERS = set(('mag',))
#: The required distance parameter is 'Joyner-Boore' distance
REQUIRES_DISTANCES = set(('rjb',))
#: GMPE not tested against independent implementation
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.
Implements equation 14 of Hong & Goda (2007)
"""
C = self.COEFFS[imt]
C_PGA = self.COEFFS[PGA()]
C_AMP = self.AMP_COEFFS[imt]
# Gets the PGA on rock - need to convert from g to cm/s/s
pga_rock = self._compute_pga_rock(C_PGA, rup.mag, dists.rjb) * 980.665
# Get the mean ground motion value
mean = (self._compute_nonlinear_magnitude_term(C, rup.mag) +
self._compute_magnitude_distance_term(C, dists.rjb, rup.mag) +
self._get_site_amplification(C_AMP, sites.vs30, pga_rock))
# Get standard deviations
stddevs = self._get_stddevs(C, stddev_types, dists.rjb.shape)
return mean, stddevs
def _compute_pga_rock(self, C_PGA, mag, rjb):
"""
Returns the PGA (g) on rock, as defined in equation 15
"""
return np.exp(self._compute_linear_magnitude_term(C_PGA, mag) +
self._compute_simple_distance_term(C_PGA, rjb))
def _compute_linear_magnitude_term(self, C, mag):
"""
Computes the linear part of the magnitude term
"""
return C["b1"] + C["b2"] * (mag - 7.0)
def _compute_nonlinear_magnitude_term(self, C, mag):
"""
Computes the non-linear magnitude term
"""
return self._compute_linear_magnitude_term(C, mag) +\
C["b3"] * ((mag - 7.0) ** 2.)
def _compute_simple_distance_term(self, C, rjb):
"""
The distance term for the PGA case ignores magnitude (equation 15)
"""
return C["b4"] * np.log(np.sqrt(rjb ** 2. + C["h"] ** 2.))
def _compute_magnitude_distance_term(self, C, rjb, mag):
"""
Returns the magntude dependent distance term
"""
rval = np.sqrt(rjb ** 2. + C["h"] ** 2.)
return (C["b4"] + C["b5"] * (mag - 4.5)) * np.log(rval)
def _get_site_amplification(self, C_AMP, vs30, pga_rock):
"""
Gets the site amplification term based on equations 7 and 8 of
Atkinson & Boore (2006)
"""
# Get nonlinear term
bnl = self._get_bnl(C_AMP, vs30)
#
f_nl_coeff = np.log(60.0 / 100.0) * np.ones_like(vs30)
idx = pga_rock > 60.0
f_nl_coeff[idx] = np.log(pga_rock[idx] / 100.0)
return np.log(np.exp(
C_AMP["blin"] * np.log(vs30 / self.CONSTS["Vref"]) +
bnl * f_nl_coeff))
def _get_bnl(self, C_AMP, vs30):
"""
Gets the nonlinear term, given by equation 8 of Atkinson & Boore 2006
"""
# Default case 8d
bnl = np.zeros_like(vs30)
if np.all(vs30 >= self.CONSTS["Vref"]):
return bnl
# Case 8a
bnl[vs30 < self.CONSTS["v1"]] = C_AMP["b1sa"]
# Cade 8b
idx = np.logical_and(vs30 > self.CONSTS["v1"],
vs30 <= self.CONSTS["v2"])
if np.any(idx):
bnl[idx] = (C_AMP["b1sa"] - C_AMP["b2sa"]) *\
(np.log(vs30[idx] / self.CONSTS["v2"]) /
np.log(self.CONSTS["v1"] / self.CONSTS["v2"])) + C_AMP["b2sa"]
# Case 8c
idx = np.logical_and(vs30 > self.CONSTS["v2"],
vs30 < self.CONSTS["Vref"])
if np.any(idx):
bnl[idx] = C_AMP["b2sa"] *\
np.log(vs30[idx] / self.CONSTS["Vref"]) /\
np.log(self.CONSTS["v2"] / self.CONSTS["Vref"])
return bnl
def _get_stddevs(self, C, stddev_types, stddev_shape):
"""
Returns the standard deviations given in Table 2
"""
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(C["sigtot"] + np.zeros(stddev_shape))
elif stddev_type == const.StdDev.INTRA_EVENT:
stddevs.append(C['sig2'] + np.zeros(stddev_shape))
elif stddev_type == const.StdDev.INTER_EVENT:
stddevs.append(C['sig1'] + np.zeros(stddev_shape))
return stddevs
COEFFS = CoeffsTable(sa_damping=5, table="""\
imt b1 b2 b3 b4 b5 h sig1 sig2 sigtot
pga 1.365 0.349 0.000 -1.123 0.062 5.9 0.184 0.449 0.485
pgv 5.540 0.931 0.000 -0.866 -0.009 3.8 0.248 0.494 0.553
0.10 2.305 -0.084 -0.054 -1.461 0.167 8.2 0.218 0.467 0.515
0.15 2.605 -0.045 -0.044 -1.514 0.179 8.6 0.218 0.473 0.521
0.20 2.514 0.234 -0.053 -1.204 0.067 8.4 0.166 0.499 0.526
0.25 2.228 0.369 0.000 -1.118 0.057 6.9 0.170 0.495 0.523
0.30 1.762 0.515 0.000 -0.878 0.003 4.7 0.182 0.510 0.541
0.40 1.608 0.577 0.000 -0.898 0.012 4.9 0.234 0.528 0.577
0.50 1.713 0.837 0.000 -0.843 -0.041 6.0 0.216 0.542 0.584
0.60 1.451 0.924 -0.030 -0.755 -0.066 5.0 0.274 0.557 0.621
0.70 1.138 0.740 -0.093 -0.838 -0.014 4.2 0.320 0.581 0.664
0.80 0.781 0.549 -0.182 -0.834 0.005 3.4 0.316 0.591 0.670
0.90 0.763 0.484 -0.197 -0.960 0.052 4.0 0.324 0.594 0.677
1.00 0.763 0.359 -0.270 -1.024 0.067 4.9 0.326 0.593 0.677
1.10 0.827 0.596 -0.333 -0.819 -0.032 4.6 0.346 0.584 0.679
1.20 0.853 0.845 -0.328 -0.689 -0.093 4.3 0.358 0.578 0.680
1.30 0.682 0.921 -0.322 -0.634 -0.108 4.0 0.362 0.576 0.681
1.40 0.540 0.954 -0.303 -0.635 -0.103 3.8 0.362 0.574 0.678
1.50 0.433 1.005 -0.294 -0.617 -0.109 3.7 0.352 0.568 0.668
1.60 0.289 0.988 -0.302 -0.617 -0.103 3.6 0.364 0.566 0.673
1.70 0.102 0.976 -0.301 -0.611 -0.093 3.4 0.378 0.564 0.679
1.80 -0.098 0.965 -0.310 -0.588 -0.088 3.0 0.380 0.560 0.676
1.90 -0.216 0.936 -0.325 -0.601 -0.079 2.9 0.380 0.561 0.677
2.00 -0.379 0.693 -0.308 -0.759 -0.001 2.7 0.364 0.562 0.670
2.20 -0.549 0.643 -0.336 -0.776 0.011 2.5 0.378 0.570 0.684
2.40 -0.663 0.772 -0.325 -0.706 -0.016 2.3 0.402 0.570 0.697
2.60 -0.747 0.909 -0.302 -0.655 -0.039 2.3 0.404 0.575 0.703
2.80 -0.883 1.024 -0.259 -0.630 -0.045 2.2 0.414 0.584 0.716
3.00 -0.955 1.027 -0.265 -0.677 -0.029 2.3 0.420 0.594 0.728
""")
AMP_COEFFS = CoeffsTable(sa_damping=5, table="""\
imt blin b1sa b2sa
pgv -0.60000 -0.49500 -0.06000
pga -0.36100 -0.64100 -0.14400
0.02500 -0.33000 -0.62400 -0.11500
0.03125 -0.32200 -0.61800 -0.10800
0.04000 -0.31400 -0.60900 -0.10500
0.05000 -0.28600 -0.64300 -0.10500
0.06289 -0.24900 -0.64200 -0.10500
0.07937 -0.23200 -0.63700 -0.11700
0.10000 -0.25000 -0.59500 -0.13200
0.12500 -0.26000 -0.56000 -0.14000
0.15873 -0.28000 -0.52800 -0.18500
0.20000 -0.30600 -0.52100 -0.18500
0.25000 -0.39000 -0.51800 -0.16000
0.31250 -0.44500 -0.51300 -0.13000
0.40000 -0.50000 -0.50800 -0.09500
0.50000 -0.60000 -0.49500 -0.06000
0.62500 -0.67000 -0.48000 -0.03100
0.76923 -0.69000 -0.46500 -0.00200
1.00000 -0.70000 -0.44000 0.00000
1.58730 -0.72600 -0.39500 0.00000
2.00000 -0.73000 -0.37500 0.00000
3.12500 -0.74000 -0.33000 0.00000
4.00000 -0.74500 -0.31000 0.00000
5.00000 -0.75200 -0.30000 0.00000
""")
CONSTS = {"Vref": 760.0,
"v1": 180.0,
"v2": 300.0}