Source code for openquake.hazardlib.gsim.yu_2013

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2014-2025 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
    """
    a1ca = np.where(mag > 6.5, coeff['ua'], coeff['a'])
    a1cb = np.where(mag > 6.5, coeff['ub'], coeff['b'])
    a1cc = np.where(mag > 6.5, coeff['uc'], coeff['c'])
    a1cd = np.where(mag > 6.5, coeff['ud'], coeff['d'])
    a1ce = np.where(mag > 6.5, coeff['ue'], coeff['e'])
    a2ca = np.where(mag > 6.5, coeff['ia'], coeff['ma'])
    a2cb = np.where(mag > 6.5, coeff['ib'], coeff['mb'])
    a2cc = np.where(mag > 6.5, coeff['ic'], coeff['mc'])
    a2cd = np.where(mag > 6.5, coeff['id'], coeff['md'])
    a2ce = np.where(mag > 6.5, coeff['ie'], 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 = {PGA, PGV, SA}

    #: Supported intensity measure component is geometric mean (supposed)
    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.GEOMETRIC_MEAN

    #: Supported standard deviation types is total
    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {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 compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        """
        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.compute>`
        for spec of input and result values.
        """
        for m, imt in enumerate(imts):
            coeff = self.COEFFS[imt]
            a1ca, a1cb, a1cc, a1cd, a1ce, a2ca, a2cb, a2cc, a2cd, a2ce = gc(
                coeff, ctx.mag)

            # Get correction coefficients. Here for each site we find the
            # the geometry of the ellipses
            ras = []
            for epi, mag, theta in zip(ctx.repi, ctx.mag, ctx.azimuth):
                res = get_ras(epi, theta, mag, coeff)
                ras.append(res)
            ras = np.array(ras)
            rbs = rbf(ras, coeff, ctx.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(ctx.azimuth)))**2
            y = (mean2 * np.cos(np.radians(ctx.azimuth)))**2
            mean_ = mean1 * mean2 / np.sqrt(x+y)
            if imt.string == "PGA":
                mean_ = np.exp(mean_)/g/100
            elif imt.string == "PGV":
                mean_ = np.exp(mean_)
            else:
                raise ValueError('Unsupported IMT')

            mean[m] = np.log(mean_)
            sig[m] = coeff['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 compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        """
        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.compute>`
        for spec of input and result values.
        """
        # Convert Mw into Ms
        mag = np.where(ctx.mag < 6.58,
                       (ctx.mag - 0.59) / 0.86,
                       (ctx.mag + 2.42) / 1.28)
        for m, imt in enumerate(imts):
            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, mg, theta in zip(ctx.repi, mag, ctx.azimuth):
                res = get_ras(epi, theta, mg, 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(ctx.azimuth)))**2
            y = (mean2 * np.cos(np.radians(ctx.azimuth)))**2
            mean_ = mean1 * mean2 / np.sqrt(x+y)
            if imt.string == "PGA":
                mean_ = np.exp(mean_)/g/100
            elif imt.string == "PGV":
                mean_ = np.exp(mean_)
            else:
                raise ValueError('Unsupported IMT')
            mean[m] = np.log(mean_)
            sig[m] = coeff['sigma']


[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
     """)