# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2014-2020 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:`YuEtAl2013Ms`, :class:`YuEtAl2013MsTibet`,
:class:`YuEtAl2013MsEastern`, :class:`YuEtAl2013MsStable`
:class:`YuEtAl2013Mw`, :class:`YuEtAl2013MwTibet`,
:class:`YuEtAl2013MwEastern`, :class:`YuEtAl2013MwStable`
"""
import numpy as np
from scipy.constants import g
from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, PGV, SA
[docs]def gc(coeff, mag):
"""
Returns the set of coefficients to be used for the calculation of GM
as a function of earthquake magnitude
:param coeff:
A dictionary of parameters for the selected IMT
:param mag:
Magnitude value
:returns:
The set of coefficients
"""
if mag > 6.5:
a1ca = coeff['ua']
a1cb = coeff['ub']
a1cc = coeff['uc']
a1cd = coeff['ud']
a1ce = coeff['ue']
a2ca = coeff['ia']
a2cb = coeff['ib']
a2cc = coeff['ic']
a2cd = coeff['id']
a2ce = coeff['ie']
else:
a1ca = coeff['a']
a1cb = coeff['b']
a1cc = coeff['c']
a1cd = coeff['d']
a1ce = coeff['e']
a2ca = coeff['ma']
a2cb = coeff['mb']
a2cc = coeff['mc']
a2cd = coeff['md']
a2ce = coeff['me']
return a1ca, a1cb, a1cc, a1cd, a1ce, a2ca, a2cb, a2cc, a2cd, a2ce
[docs]def rbf(ra, coeff, mag):
"""
Calculate the median ground motion for a given magnitude and distance
:param ra:
Distance value [km]
:param coeff:
The set of coefficients
:param mag:
Magnitude value
:returns:
"""
a1ca, a1cb, a1cc, a1cd, a1ce, a2ca, a2cb, a2cc, a2cd, a2ce = gc(coeff, mag)
term1 = a1ca + a1cb * mag + a1cc * np.log(ra + a1cd*np.exp(a1ce*mag))
term2 = a2ca + a2cb * mag
term3 = a2cd*np.exp(a2ce*mag)
return np.exp((term1 - term2) / a2cc) - term3
[docs]def fnc(ra, *args):
"""
Function used in the minimisation problem.
:param ra:
Semi-axis of the ellipses used in the Yu et al.
:returns:
The absolute difference between the epicentral distance and the
adjusted distance
"""
#
# epicentral distance
repi = args[0]
#
# azimuth
theta = args[1]
#
# magnitude
mag = args[2]
#
# coefficients
coeff = args[3]
#
# compute the difference between epicentral distances
rb = rbf(ra, coeff, mag)
t1 = ra**2 * (np.sin(np.radians(theta)))**2
t2 = rb**2 * (np.cos(np.radians(theta)))**2
xx = ra * rb / (t1+t2)**0.5
return xx-repi
[docs]def get_ras(repi, theta, mag, coeff):
"""
Computes equivalent distance
:param repi:
Epicentral distance
:param theta:
Azimuth value
:param mag:
Magnitude
:param coeff:
GMPE coefficients
"""
rx = 100.
ras = 200.
#
# calculate the difference between epicentral distances
dff = fnc(ras, repi, theta, mag, coeff)
while abs(dff) > 1e-3:
# update the value of distance computed
if dff > 0.:
ras = ras - rx
else:
ras = ras + rx
dff = fnc(ras, repi, theta, mag, coeff)
rx = rx / 2.
if rx < 1e-3:
break
return ras
[docs]class YuEtAl2013Ms(GMPE):
"""
Implements the Yu et al. (2013) GMPE used for the calculation of the 2015
version of the national seismic hazard maps for China. Note that magnitude
supported is Ms.
"""
#: Supported tectonic region type is active shallow crust
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST
#: Supported intensity measure types are peak ground velocity and
#: peak ground acceleration
DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([
PGA,
PGV,
SA
])
#: Supported intensity measure component is geometric mean (supposed)
DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.AVERAGE_HORIZONTAL
#: Supported standard deviation types is total
DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([
const.StdDev.TOTAL
])
#: No site parameters required
REQUIRES_SITES_PARAMETERS = set()
#: Required rupture parameter is magnitude
REQUIRES_RUPTURE_PARAMETERS = {'mag'}
#: Required distance measures are epicentral distance and azimuth
REQUIRES_DISTANCES = {'repi', 'azimuth'}
[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.
"""
# Check that the requested standard deviation type is available
assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
for stddev_type in stddev_types)
#
# Set parameters
mag = rup.mag
epi = dists.repi
theta = dists.azimuth
#
# Set coefficients
coeff = self.COEFFS[imt]
a1ca, a1cb, a1cc, a1cd, a1ce, a2ca, a2cb, a2cc, a2cd, a2ce = \
gc(coeff, mag)
#
# Get correction coefficients. Here for each site we find the
# the geometry of the ellipses
ras = []
for epi, theta in zip(dists.repi, dists.azimuth):
res = get_ras(epi, theta, mag, coeff)
ras.append(res)
ras = np.array(ras)
rbs = rbf(ras, coeff, mag)
#
# Compute values of ground motion for the two cases. The value of
# 225 is hardcoded under the assumption that the hypocentral depth
# corresponds to 15 km (i.e. 15**2)
mean1 = (a1ca + a1cb * mag +
a1cc * np.log((ras**2+225)**0.5 +
a1cd * np.exp(a1ce * mag)))
mean2 = (a2ca + a2cb * mag +
a2cc * np.log((rbs**2+225)**0.5 +
a2cd * np.exp(a2ce * mag)))
#
# Get distances
x = (mean1 * np.sin(np.radians(dists.azimuth)))**2
y = (mean2 * np.cos(np.radians(dists.azimuth)))**2
mean = mean1 * mean2 / np.sqrt(x+y)
if imt.name == "PGA":
mean = np.exp(mean)/g/100
elif imt.name == "PGV":
mean = np.exp(mean)
else:
raise ValueError('Unsupported IMT')
#
# Get the standard deviation
stddevs = self._compute_std(coeff, stddev_types, len(dists.repi))
#
# Return results
return np.log(mean), stddevs
def _compute_std(self, C, stddev_types, num_sites):
return [np.ones(num_sites)*C['sigma']]
#: Coefficient table
COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT a b c d e ua ub uc ud ue ma mb mc md me ia ib ic id ie sigma
PGA 4.1193 1.656 -2.389 1.772 0.424 7.8269 1.0856 -2.389 1.772 0.424 2.2609 1.6399 -2.118 0.825 0.465 6.003 1.0649 -2.118 0.825 0.465 0.5428
PGV -1.2581 1.932 -2.181 1.772 0.424 3.013 1.2742 -2.181 1.772 0.424 -3.1073 1.9389 -1.945 0.825 0.465 1.3087 1.2627 -1.945 0.825 0.465 0.6233
""")
[docs]class YuEtAl2013MsTibet(YuEtAl2013Ms):
#: Supported tectonic region type is Tibetan plateau
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST
#: Coefficient table
COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT a b c d e ua ub uc ud ue ma mb mc md me ia ib ic id ie sigma
PGA 5.4901 1.4835 -2.416 2.647 0.366 8.7561 0.9453 -2.416 2.647 0.366 2.3069 1.4007 -1.854 0.612 0.457 5.6511 0.8924 -1.854 0.612 0.457 0.5428
PGV -0.1472 1.7618 -2.205 2.647 0.366 3.9422 1.1293 -2.205 2.647 0.366 -2.9923 1.7043 -1.696 0.612 0.457 1.0189 1.0902 -1.696 0.612 0.457 0.6233
""")
[docs]class YuEtAl2013MsEastern(YuEtAl2013Ms):
#: Supported tectonic region type is eastern part of China
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.STABLE_CONTINENTAL
#: Coefficient table
COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT a b c d e ua ub uc ud ue ma mb mc md me ia ib ic id ie sigma
PGA 4.5517 1.5433 -2.315 2.088 0.399 8.1259 0.9936 -2.315 2.088 0.399 2.7048 1.518 -2.004 0.944 0.447 6.3319 0.9614 -2.004 0.944 0.447 0.5428
PGV -0.8349 1.8193 -2.103 2.088 0.399 3.3051 1.1799 -2.103 2.088 0.399 -2.6381 1.8124 -1.825 0.944 0.447 1.6376 1.1546 -1.825 0.944 0.447 0.6233
""")
[docs]class YuEtAl2013MsStable(YuEtAl2013Ms):
#: Supported tectonic region type is stable part of China
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.STABLE_CONTINENTAL
#: Coefficient table
COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT a b c d e ua ub uc ud ue ma mb mc md me ia ib ic id ie sigma
PGA 5.5591 1.1454 -2.079 2.802 0.295 8.5238 0.6854 -2.079 2.802 0.295 3.9445 1.0833 -1.723 1.295 0.331 6.187 0.7383 -1.723 1.295 0.331 0.5428
PGV 0.2139 1.4283 -1.889 2.802 0.295 3.772 0.8786 -1.889 2.802 0.295 -1.3547 1.3823 -1.559 1.295 0.331 1.5433 0.9361 -1.559 1.295 0.331 0.6233
""")
[docs]class YuEtAl2013Mw(YuEtAl2013Ms):
"""
This is a modified version of the original Yu et al. (2013) that supports
the use of Mw rather than Ms. The Mw to Ms conversion equation used is the
one proposed by Cheng et al. (2017). Note that this version does not
propagate the uncertainty related to the magnitude conversion process.
"""
[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.
"""
# Check that the requested standard deviation type is available
assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
for stddev_type in stddev_types)
#
# Set parameters
magn = rup.mag
epi = dists.repi
theta = dists.azimuth
#
# Convert Mw into Ms
if magn < 6.58:
mag = (magn - 0.59) / 0.86
else:
mag = (magn + 2.42) / 1.28
#
# Set coefficients
coeff = self.COEFFS[imt]
a1ca, a1cb, a1cc, a1cd, a1ce, a2ca, a2cb, a2cc, a2cd, a2ce = \
gc(coeff, mag)
#
# Get correction coefficients. Here for each site we find the
# the geometry of the ellipses
ras = []
for epi, theta in zip(dists.repi, dists.azimuth):
res = get_ras(epi, theta, mag, coeff)
ras.append(res)
ras = np.array(ras)
rbs = rbf(ras, coeff, mag)
#
# Compute values of ground motion for the two cases. The value of
# 225 is hardcoded under the assumption that the hypocentral depth
# corresponds to 15 km (i.e. 15**2)
mean1 = (a1ca + a1cb * mag +
a1cc * np.log((ras**2+225)**0.5 +
a1cd * np.exp(a1ce * mag)))
mean2 = (a2ca + a2cb * mag +
a2cc * np.log((rbs**2+225)**0.5 +
a2cd * np.exp(a2ce * mag)))
#
# Get distances
x = (mean1 * np.sin(np.radians(dists.azimuth)))**2
y = (mean2 * np.cos(np.radians(dists.azimuth)))**2
mean = mean1 * mean2 / np.sqrt(x+y)
if imt.name == "PGA":
mean = np.exp(mean)/g/100
elif imt.name == "PGV":
mean = np.exp(mean)
else:
raise ValueError('Unsupported IMT')
#
# Get the standard deviation
stddevs = self._compute_std(coeff, stddev_types, len(dists.repi))
#
# Return results
return np.log(mean), stddevs
[docs]class YuEtAl2013MwTibet(YuEtAl2013Mw):
#: Supported tectonic region type is Tibetan plateau
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST
#: Coefficient table
COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT a b c d e ua ub uc ud ue ma mb mc md me ia ib ic id ie sigma
PGA 5.4901 1.4835 -2.416 2.647 0.366 8.7561 0.9453 -2.416 2.647 0.366 2.3069 1.4007 -1.854 0.612 0.457 5.6511 0.8924 -1.854 0.612 0.457 0.5428
PGV -0.1472 1.7618 -2.205 2.647 0.366 3.9422 1.1293 -2.205 2.647 0.366 -2.9923 1.7043 -1.696 0.612 0.457 1.0189 1.0902 -1.696 0.612 0.457 0.6233
""")
[docs]class YuEtAl2013MwEastern(YuEtAl2013Mw):
#: Supported tectonic region type is eastern part of China
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.STABLE_CONTINENTAL
#: Coefficient table
COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT a b c d e ua ub uc ud ue ma mb mc md me ia ib ic id ie sigma
PGA 4.5517 1.5433 -2.315 2.088 0.399 8.1259 0.9936 -2.315 2.088 0.399 2.7048 1.518 -2.004 0.944 0.447 6.3319 0.9614 -2.004 0.944 0.447 0.5428
PGV -0.8349 1.8193 -2.103 2.088 0.399 3.3051 1.1799 -2.103 2.088 0.399 -2.6381 1.8124 -1.825 0.944 0.447 1.6376 1.1546 -1.825 0.944 0.447 0.6233
""")
[docs]class YuEtAl2013MwStable(YuEtAl2013Mw):
#: Supported tectonic region type is stable part of China
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.STABLE_CONTINENTAL
#: Coefficient table
COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT a b c d e ua ub uc ud ue ma mb mc md me ia ib ic id ie sigma
PGA 5.5591 1.1454 -2.079 2.802 0.295 8.5238 0.6854 -2.079 2.802 0.295 3.9445 1.0833 -1.723 1.295 0.331 6.187 0.7383 -1.723 1.295 0.331 0.5428
PGV 0.2139 1.4283 -1.889 2.802 0.295 3.772 0.8786 -1.889 2.802 0.295 -1.3547 1.3823 -1.559 1.295 0.331 1.5433 0.9361 -1.559 1.295 0.331 0.6233
""")