Source code for openquake.hazardlib.gsim.sera_amplification_models

# -*- 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/>.

"""
Implements SERA site amplification models class: `PitilakisEtAl2018`,
                                                 `PitilakisEtAl2020`,
                                                 `Eurocode8Amplification`,
                                                 `Eurocode8AmplificationDefault`,
                                                 `SandikkayaDinsever2018`
"""
import numpy as np
import copy
from scipy.constants import g
# from scipy.interpolate import interp1d
from openquake.hazardlib.gsim.base import (GMPE, CoeffsTable, registry)
from openquake.hazardlib.imt import PGA, SA, from_string
from openquake.hazardlib import const


# Pitilakis GMPE Wrapper
uppernames = '''
DEFINED_FOR_INTENSITY_MEASURE_TYPES
DEFINED_FOR_STANDARD_DEVIATION_TYPES
REQUIRES_SITES_PARAMETERS
REQUIRES_RUPTURE_PARAMETERS
REQUIRES_DISTANCES
'''.split()


[docs]class PitilakisEtAl2018(GMPE): """ Implements a site amplification model based on a design code approach, using the site characterisation and amplification coefficients proposed by Pitilakis et al. (2018) Pitilakis, K., Riga, E., Anastasiadis, A., Fotopoulou, S. and Karafagka, S. (2018) "Towards the revision of EC8: Proposal for an alternative site classification scheme and associated intensity dependent spectral amplification factors", Soil Dynamics & Earthquake Engineering, Care should be taken to note the following: 1. In the absence of a specific guidance from Eurocode 8 as to how the short period coefficient SS is determine from the UHS the choice is made to anchor the short period spectrum to PGA, with SS taken as being equal to 2.5 * PGA. This is implied by the Eurocode 8 decision to fix F0 to 2.5 and that the ground motion is fixed to SS / F0 for T -> 0 2. No uncertainty in amplification factor is presented in a code based approach and therefore the standard deviation of the original GMPE is retained. :param gmpe: Input ground motion prediction equation :param float rock_vs30: Reference shearwave velocity used for the rock calculation """ experimental = True #: Supported tectonic region type is undefined (applies to any) DEFINED_FOR_TECTONIC_REGION_TYPE = "" #: Supported intensity measure types are not set DEFINED_FOR_INTENSITY_MEASURE_TYPES = set((PGA, SA)) #: Supported intensity measure component is horizontal #: :attr:`~openquake.hazardlib.const.IMC.HORIZONTAL`, DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.HORIZONTAL #: Supported standard deviation type DEFINED_FOR_STANDARD_DEVIATION_TYPES = set((const.StdDev.TOTAL,)) #: Required site parameters are Vs30 and the Pitilakis et al (2018) site #: class (others will be added for the GMPE in question) REQUIRES_SITES_PARAMETERS = {'vs30', 'ec8_p18'} #: Required rupture parameter is magnitude, others will be set later REQUIRES_RUPTURE_PARAMETERS = {'mag'} #: Required distance metrics will be set by the GMPEs REQUIRES_DISTANCES = set() #: Defined reference velocity is 800 m/s DEFINED_FOR_REFERENCE_VELOCITY = 800.0 def __init__(self, gmpe_name, reference_velocity=None, **kwargs): super().__init__() if isinstance(gmpe_name, str): self.gmpe = registry[gmpe_name](**kwargs) else: # An instantiated class is passed as an argument self.gmpe = copy.deepcopy(gmpe_name) if reference_velocity: self.rock_vs30 = reference_velocity else: self.rock_vs30 = self.DEFINED_FOR_REFERENCE_VELOCITY for name in uppernames: setattr(self, name, frozenset(getattr(self, name) | getattr(self.gmpe, name)))
[docs] def get_mean_and_stddevs(self, sctx, rctx, dctx, imt, stddev_types): """ Returns the mean and standard deviations calling the input GMPE for the mean acceleration for PGA and Sa (1.0) on the reference rock, defining the amplification factors and code spectrum to return the mean ground motion at the desired period, before the calling the input GMPE once more in order to return the standard deviations for the required IMT. """ sctx_r = copy.copy(sctx) sctx_r.vs30 = self.rock_vs30 * np.ones_like(sctx_r.vs30) # Get PGA and Sa (1.0) from GMPE pga_r = self.gmpe.get_mean_and_stddevs(sctx_r, rctx, dctx, PGA(), stddev_types)[0] s_1_rp = self.gmpe.get_mean_and_stddevs(sctx_r, rctx, dctx, SA(1.0), stddev_types)[0] s_s_rp = self.CONSTANTS["F0"] * np.exp(pga_r) s_1_rp = np.exp(s_1_rp) # Get the short and long period amplification factors f_s, f_l = self.get_amplification_factor(s_s_rp, s_1_rp, sctx) s_1 = f_l * s_1_rp s_s = f_s * s_s_rp # Get the mean ground motion at the IMT using the design code spectrum mean = self.get_amplified_mean(s_s, s_1, s_1_rp, imt) # Call the original GMPE to return the standard deviation for the # IMT in question stddevs = self.gmpe.get_mean_and_stddevs(sctx, rctx, dctx, imt, stddev_types)[1] return mean, stddevs
[docs] def get_amplification_factor(self, s_s, s_1, sctx): """ Returns the short and long-period amplification factors given the input Pitilakis et al. (2018) site class and the short and long-period input accelerations """ f_s = np.ones(sctx.ec8_p18.shape, dtype=float) f_l = np.ones(sctx.ec8_p18.shape, dtype=float) for ec8b in np.unique(sctx.ec8_p18): ec8 = ec8b.decode('ascii') if ec8 == "A": # Amplification factors are 1 continue idx = sctx.ec8_p18 == ec8b if np.any(idx): s_ss = s_s[idx] f_ss = np.ones(np.sum(idx)) f_ls = np.ones(np.sum(idx)) lb = s_ss < 0.25 ub = s_ss > 1.25 f_ss[lb] = self.FS[ec8][0] f_ls[lb] = self.F1[ec8][0] f_ss[ub] = self.FS[ec8][-1] f_ls[ub] = self.F1[ec8][-1] for j in range(1, len(self.IMLS) - 1): jdx = np.logical_and(s_ss >= self.IMLS[j], s_ss < self.IMLS[j + 1]) if not np.any(jdx): continue dfs = self.FS[ec8][j + 1] - self.FS[ec8][j] dfl = self.F1[ec8][j + 1] - self.F1[ec8][j] diml = self.IMLS[j + 1] - self.IMLS[j] f_ss[jdx] = self.FS[ec8][j] + (s_ss[jdx] - self.IMLS[j]) *\ (dfs / diml) f_ls[jdx] = self.F1[ec8][j] + (s_ss[jdx] - self.IMLS[j]) *\ (dfl / diml) f_s[idx] = f_ss f_l[idx] = f_ls return f_s, f_l
[docs] def get_amplified_mean(self, s_s, s_1, s_1_rp, imt): """ Given the amplified short- and long-period input accelerations, returns the mean ground motion for the IMT according to the design spectrum construction in equations 1 - 5 of Pitilakis et al., (2018) """ if "PGA" in str(imt) or imt.period <= self.CONSTANTS["TA"]: # PGA or v. short period acceleration return np.log(s_s / self.CONSTANTS["F0"]) mean = np.copy(s_s) t_c = s_1 / s_s t_b = t_c / self.CONSTANTS["kappa"] t_b[t_b < 0.05] = 0.05 t_b[t_b > 0.1] = 0.1 t_d = 2.0 + np.zeros_like(s_1_rp) idx = s_1_rp > 0.1 if np.any(idx): t_d[idx] = 1.0 + (10. * s_1_rp[idx]) idx = np.logical_and(self.CONSTANTS["TA"] < imt.period, t_b >= imt.period) if np.any(idx): mean[idx] = (s_s[idx] / (t_b[idx] - self.CONSTANTS["TA"])) *\ ((imt.period - self.CONSTANTS["TA"]) + (t_b[idx] - imt.period) / self.CONSTANTS["F0"]) idx = np.logical_and(t_c < imt.period, t_d >= imt.period) if np.any(idx): mean[idx] = s_1[idx] / imt.period idx = t_d < imt.period if np.any(idx): mean[idx] = t_d[idx] * s_1[idx] / (imt.period ** 2.) return np.log(mean)
CONSTANTS = { "F0": 2.5, "kappa": 5.0, "TA": 0.03} IMLS = [0., 0.25, 0.5, 0.75, 1., 1.25] # Short period amplification factors defined by Pitilakis et al., (2018) FS = { "A": [1.0, 1.0, 1.0, 1.0, 1.0, 1.0], "B1": [1.3, 1.3, 1.2, 1.2, 1.2, 1.2], "B2": [1.4, 1.3, 1.3, 1.2, 1.1, 1.1], "C1": [1.7, 1.6, 1.4, 1.3, 1.3, 1.2], "C2": [1.6, 1.5, 1.3, 1.2, 1.1, 1.0], "C3": [1.8, 1.6, 1.4, 1.2, 1.1, 1.0], "D": [2.2, 1.9, 1.6, 1.4, 1.2, 1.0], "E": [1.7, 1.6, 1.6, 1.5, 1.5, 1.5]} # Long period amplification factors defined by Pitilakis et al., (2018) F1 = { "A": [1.0, 1.0, 1.0, 1.0, 1.0, 1.0], "B1": [1.4, 1.4, 1.4, 1.4, 1.3, 1.3], "B2": [1.6, 1.5, 1.5, 1.5, 1.4, 1.3], "C1": [1.7, 1.6, 1.5, 1.5, 1.4, 1.3], "C2": [2.1, 2.0, 1.9, 1.8, 1.8, 1.7], "C3": [3.2, 3.0, 2.7, 2.5, 2.4, 2.3], "D": [4.1, 3.8, 3.3, 3.0, 2.8, 2.7], "E": [1.3, 1.3, 1.2, 1.2, 1.2, 1.2]}
[docs]class PitilakisEtAl2020(PitilakisEtAl2018): """ Adaptation of the Pitilakis et al. (2018) amplification model adopting the revised FS and F1 factors proposed by Pitilakis et al., (2020) Pitilakis, K., Riga, E. and Anastasiadis, A. (2020), Towards the Revision of EC8: Proposal for an Alternative Site Classification Scheme and Associated Intensity-Dependent Amplification Factors. In the Proceedings of the 17th World Conference on Earthquake Engineering, 17WCEE, Sendai, Japan, September 13th to 18th 2020. Paper No. C002895. """ experimental = True # Short period amplification factors defined by Pitilakis et al., (2020) FS = { "A": [1.0, 1.0, 1.0, 1.0, 1.0, 1.0], "B1": [1.3, 1.3, 1.3, 1.2, 1.2, 1.2], "B2": [1.3, 1.3, 1.2, 1.2, 1.2, 1.1], "C1": [1.7, 1.7, 1.6, 1.5, 1.5, 1.4], "C2": [1.6, 1.5, 1.3, 1.2, 1.1, 1.0], "C3": [1.7, 1.6, 1.4, 1.2, 1.2, 1.1], "D": [1.8, 1.7, 1.5, 1.4, 1.3, 1.2], "E": [1.7, 1.6, 1.6, 1.5, 1.5, 1.4]} # Long period amplification factors defined by Pitilakis et al., (2020) F1 = { "A": [1.0, 1.0, 1.0, 1.0, 1.0, 1.0], "B1": [1.1, 1.1, 1.1, 1.1, 1.1, 1.1], "B2": [1.4, 1.4, 1.3, 1.3, 1.3, 1.3], "C1": [1.5, 1.5, 1.4, 1.4, 1.4, 1.4], "C2": [2.3, 2.2, 2.0, 1.9, 1.9, 1.8], "C3": [2.4, 2.3, 2.1, 2.0, 2.0, 1.9], "D": [4.0, 3.5, 3.0, 2.7, 2.4, 2.3], "E": [1.2, 1.1, 1.1, 1.1, 1.1, 1.1]}
[docs]class Eurocode8Amplification(PitilakisEtAl2018): """ Implements a general class to return a ground motion based on the Eurocode 8 design spectrum: CEN (2018): "Eurocode 8: Earthquake Resistant Design of Structures" Revised 2nd Draft SC8 PT1 - Rev 20 The potential notes highlighted in :class:`PitilakisEtAl2018` apply in this case too. """ experimental = True #: Supported tectonic region type is undefined (applies to any) DEFINED_FOR_TECTONIC_REGION_TYPE = "" #: Supported intensity measure types are not set DEFINED_FOR_INTENSITY_MEASURE_TYPES = set((PGA, SA)) #: Supported intensity measure component is horizontal #: :attr:`~openquake.hazardlib.const.IMC.HORIZONTAL`, DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.HORIZONTAL #: Supported standard deviation type DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([const.StdDev.TOTAL]) #: Required site parameters will be set be selected GMPES REQUIRES_SITES_PARAMETERS = {'vs30', 'h800'} #: Required rupture parameter is magnitude, others will be set later REQUIRES_RUPTURE_PARAMETERS = {'mag'} #: Required distance metrics will be set by the GMPEs REQUIRES_DISTANCES = set() #: Defined reference velocity is 800 m/s DEFINED_FOR_REFERENCE_VELOCITY = 800.0 def __init__(self, gmpe_name, reference_velocity=800.0, **kwargs): super().__init__(gmpe_name=gmpe_name) self.rock_vs30 = reference_velocity if reference_velocity else\ self.DEFINED_FOR_REFERENCE_VELOCITY for name in uppernames: setattr(self, name, frozenset(getattr(self, name) | getattr(self.gmpe, name)))
[docs] def get_mean_and_stddevs(self, sctx, rctx, dctx, imt, stddev_types): """ As with the :class:`PitilakisEtal2018`, the mean ground motion is determined by construction of the Eurocode 8 design spectrum from the short- and long-period acceleration coefficients amplified to the desired site class, with the standard deviations taken from the original GMPE at the desired IMT """ sctx_r = copy.copy(sctx) sctx_r.vs30 = self.rock_vs30 * np.ones_like(sctx_r.vs30) # Get PGA and Sa (1.0) from GMPE pga_r = self.gmpe.get_mean_and_stddevs(sctx_r, rctx, dctx, PGA(), stddev_types)[0] s_1_rp = self.gmpe.get_mean_and_stddevs(sctx_r, rctx, dctx, SA(1.0), stddev_types)[0] s_s_rp = self.CONSTANTS["F0"] * np.exp(pga_r) s_1_rp = np.exp(s_1_rp) ec8 = self.get_ec8_class(sctx.vs30, sctx.h800) f_s, f_l = self.get_amplification_factor(s_s_rp, s_1_rp, ec8, sctx) s_1 = f_l * s_1_rp s_s = f_s * s_s_rp mean = self.get_amplified_mean(s_s, s_1, s_1_rp, imt) stddevs = self.gmpe.get_mean_and_stddevs(sctx, rctx, dctx, imt, stddev_types)[1] return mean, stddevs
[docs] def get_amplification_factor(self, s_s_rp, s_1_rp, ec8, sctx): """ Returns the amplification factors based on the proposed EC8 formulation in Table 3.4 """ r_alpha = 1.0 - 2.0E3 * ((s_s_rp * g) / (sctx.vs30 ** 2)) r_beta = 1.0 - 2.0E3 * ((s_1_rp * g) / (sctx.vs30 ** 2)) f_s = np.ones(sctx.vs30.shape) f_l = np.ones(sctx.vs30.shape) vsh_norm = sctx.vs30 / 800. for s_c in np.unique(ec8): if s_c == b"A": continue idx = ec8 == s_c if not np.any(idx): continue if s_c in (b"B", b"C", b"D"): f_s[idx] = vsh_norm[idx] ** (-0.25 * r_alpha[idx]) f_l[idx] = vsh_norm[idx] ** (-0.7 * r_beta[idx]) elif s_c == b"E": f_s[idx] = vsh_norm[idx] ** (-0.25 * r_alpha[idx] * (sctx.h800[idx] / 30.) * (4.0 - (sctx.h800[idx] / 10.))) f_l[idx] = vsh_norm[idx] ** (-0.7 * r_beta[idx] * (sctx.h800[idx] / 30.)) elif s_c == b"F": f_s[idx] = 0.9 * (vsh_norm[idx] ** (-0.25 * r_alpha[idx])) f_l[idx] = 1.25 * (vsh_norm[idx] ** (-0.7 * r_beta[idx])) else: pass return f_s, f_l
[docs] @staticmethod def get_ec8_class(vsh, h800): """ Method to return the vector of Eurocode 8 site classes based on Vs30 and h800 """ ec8 = np.array([b"A" for i in range(len(vsh))], dtype="|S1") idx = np.logical_and(vsh >= 400., h800 > 5.) ec8[idx] = b"B" idx1 = np.logical_and(vsh < 400., vsh >= 250.) idx = np.logical_and(idx1, np.logical_and(h800 > 5., h800 <= 30.)) ec8[idx] = b"E" idx = np.logical_and(idx1, np.logical_and(h800 > 30., h800 <= 100.)) ec8[idx] = b"C" idx = np.logical_and(idx1, h800 > 100.) ec8[idx] = b"F" idx1 = vsh < 250. idx = np.logical_and(idx1, h800 <= 30.) ec8[idx] = b"E" idx = np.logical_and(idx1, np.logical_and(h800 > 30., h800 <= 100.)) ec8[idx] = b"D" idx = np.logical_and(idx1, h800 > 100.) ec8[idx] = b"F" return ec8
# Default short period amplification factors defined by Eurocode 8 Table 3.4 EC8_FS_default = { b"A": 1., b"B": 1.2, b"C": 1.35, b"D": 1.5, b"E": 1.7, b"F": 1.35 } # Default long period amplification factors defined by Eurocode 8 Table 3.4 EC8_FL_default = { b"A": 1., b"B": 1.6, b"C": 2.25, b"D": 3.20, b"E": 3.0, b"F": 4.0 }
[docs]class Eurocode8AmplificationDefault(Eurocode8Amplification): """ In the case that Vs30 and h800 are not known but a Eurocode 8 site class is otherwise determined then a set of default amplification factors are applied. This model implements the Eurocode 8 design spectrum """ experimental = True #: Required site parameters are the EC8 site class, everything else will #: be set be selected GMPES REQUIRES_SITES_PARAMETERS = set(('ec8',))
[docs] def get_mean_and_stddevs(self, sctx, rctx, dctx, imt, stddev_types): """ Returns the mean and standard deviations following the approach in :class:`Eurocode8Amplification` """ sctx_r = copy.copy(sctx) sctx_r.vs30 = self.rock_vs30 * np.ones_like(sctx_r.vs30) # Get PGA and Sa (1.0) from GMPE pga_r = self.gmpe.get_mean_and_stddevs(sctx_r, rctx, dctx, PGA(), stddev_types)[0] s_1_rp = self.gmpe.get_mean_and_stddevs(sctx_r, rctx, dctx, SA(1.0), stddev_types)[0] s_s_rp = self.CONSTANTS["F0"] * np.exp(pga_r) s_1_rp = np.exp(s_1_rp) f_s, f_l = self.get_amplification_factor(s_s_rp, s_1_rp, sctx) s_1 = f_l * s_1_rp s_s = f_s * s_s_rp mean = self.get_amplified_mean(s_s, s_1, s_1_rp, imt) stddevs = self.gmpe.get_mean_and_stddevs(sctx, rctx, dctx, imt, stddev_types)[1] return mean, stddevs
[docs] def get_amplification_factor(self, s_s_rp, s_1_rp, sctx): """ Returns the default amplification factor dependent upon the site class """ f_s = np.ones(sctx.ec8.shape) f_l = np.ones(sctx.ec8.shape) for key in EC8_FS_default: idx = sctx.ec8 == key if np.any(idx): f_s[idx] = EC8_FS_default[key] f_l[idx] = EC8_FL_default[key] return f_s, f_l
# Sandikkaya & Dinsever REGION_SET = ["USNZ", "JP", "TW", "CH", "WA", "TRGR", "WMT", "NWE"]
[docs]class SandikkayaDinsever2018(GMPE): """ Implements the nonlinear site amplification model of Sandikkaya & Dinsever (2018), see Sandikkaya, M. A. and Dinsever, L. D. (2018) "A Site Amplification Model for Crustal Earthquakes", Geosciences, 264(8), doi:10.3390/geosciences8070264 Note that the nonlinear amplification model has its own standard deviation, which should be applied with the phi0 model of the original GMPE. This is not defined for all GMPEs in the literature, nor is the retrieval of it consistently applied in OpenQuake. Therefore we allow the user to define manually the input phi0 model, and if this is not possible a "default" phi0 is taken by reducing the original GMPE's phi by 15 %. The amplification model is compatible only with GMPEs with separate inter- and intra-event standard deviation, otherwise an error is raised. :param gmpe: Input GMPE for calculation on reference rock and standrd deviation at the period of interest on surface rock :param phi_0: Single-station within-event standard deviation (as a period-dependent dictionary or None) :param str region: Defines the region for the region-adjusted version of the model """ experimental = True #: Supported tectonic region type is undefined DEFINED_FOR_TECTONIC_REGION_TYPE = "Active Shallow Crust" #: Supported intensity measure types are not set DEFINED_FOR_INTENSITY_MEASURE_TYPES = set((PGA, SA)) #: Supported intensity measure component is horizontal #: :attr:`~openquake.hazardlib.const.IMC.HORIZONTAL`, DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.HORIZONTAL #: Supported standard deviation type DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([const.StdDev.TOTAL]) #: Required site parameters will be set be selected GMPES REQUIRES_SITES_PARAMETERS = {'vs30', 'z1pt0'} #: Required rupture parameter is magnitude, others will be set later REQUIRES_RUPTURE_PARAMETERS = {'mag'} #: Required distance metrics will be set by the GMPEs REQUIRES_DISTANCES = set() #: Defined reference velocity is 800 m/s DEFINED_FOR_REFERENCE_VELOCITY = 760.0 def __init__(self, gmpe_name, reference_velocity=760., region=None, phi_0=None, **kwargs): super().__init__() if isinstance(gmpe_name, str): self.gmpe = registry[gmpe_name](**kwargs) else: # An instantiated class is passed as an argument self.gmpe = copy.deepcopy(gmpe_name) # Define the reference velocity - set to 760. by default self.rock_vs30 = reference_velocity if reference_velocity else\ self.DEFINED_FOR_REFERENCE_VELOCITY for name in uppernames: setattr(self, name, frozenset(getattr(self, name) | getattr(self.gmpe, name))) stddev_check = (const.StdDev.INTER_EVENT in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES) and\ (const.StdDev.INTRA_EVENT in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES) if not stddev_check: raise ValueError("Input GMPE %s not defined for inter- and intra-" "event standard deviation" % str(self.gmpe)) if isinstance(phi_0, dict): # Input phi_0 model iphi_0 = {} for key in phi_0: iphi_0[from_string(key)] = phi_0[key] self.phi_0 = CoeffsTable(sa_damping=5, table=iphi_0) else: # No input phi_0 model self.phi_0 = None # Regionalisation of the linear site term is possible # check if region is in the set of supported terms and # raise error otherwise if region is not None: if region in REGION_SET: self.region = "ck{:s}".format(region) else: raise ValueError("Region must be one of: %s" % " ".join(REGION_SET)) else: self.region = region
[docs] def get_mean_and_stddevs(self, sctx, rctx, dctx, imt, stddev_types): """ Returns the mean and standard deviations """ sctx_r = copy.copy(sctx) sctx_r.vs30 = self.rock_vs30 * np.ones_like(sctx_r.vs30) mean, stddevs = self.gmpe.get_mean_and_stddevs( sctx_r, rctx, dctx, imt, [const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT]) psarock = np.exp(mean) C = self.COEFFS_SITE[imt] if self.region: ck = self.COEFFS_REG[imt][self.region] else: ck = 0.0 ampl = self.get_site_amplification(C, psarock, sctx, ck) mean += ampl stddevs = self.get_stddevs(C, stddevs, psarock, sctx.vs30, imt, stddev_types) return mean, stddevs
[docs] def get_site_amplification(self, C, psarock, sites, ck): """ Returns the site amplification model define in equation (9) """ vs30_s = np.copy(sites.vs30) vs30_s[vs30_s > 1000.] = 1000. fn_lin = (C["b1"] + ck) * np.log(vs30_s / 760.) fn_z = C["b2"] * np.log(sites.z1pt0) fn_nl = C["b3"] * np.log((psarock + 0.1 * g) / (0.1 * g)) *\ np.exp(-np.exp(2.0 * np.log(sites.vs30) - 11.)) return fn_lin + fn_z + fn_nl
[docs] def get_stddevs(self, C, istddevs, psa_rock, vs30, imt, stddev_types): """ Returns the standard deviation adjusted for the site-response model """ tau, phi = istddevs ysig = np.copy(psa_rock) ysig[ysig > 0.35] = 0.35 ysig[ysig < 0.005] = 0.005 vsig = np.copy(vs30) vsig[vsig > 600.0] = 600.0 vsig[vsig < 150.] = 150. sigma_s = C["sigma_s"] * C["c0"] * (C["c1"] * np.log(ysig) + C["c2"] * np.log(vsig)) if self.phi_0: phi0 = self.phi_0[imt] + np.zeros(vs30.shape) else: # In the case that no input phi0 is defined take 'approximate' # phi0 as 85 % of phi phi0 = 0.85 * phi phi = np.sqrt(phi0 ** 2. + sigma_s ** 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(np.sqrt(tau ** 2. + phi ** 2.) + np.zeros(vs30.shape)) elif stddev_type == const.StdDev.INTRA_EVENT: stddevs.append(phi) elif stddev_type == const.StdDev.INTER_EVENT: stddevs.append(tau + np.zeros(vs30.shape)) return stddevs
COEFFS_SITE = CoeffsTable(sa_damping=5, table="""\ imt b1 b3 b2 sigma_s c0 c2 c1 pga -0.53307 -0.46412 0.02105 0.47096 1.24013 0.09542 -0.05865 0.010 -0.53307 -0.46412 0.02105 0.47096 1.24013 0.09542 -0.05865 0.025 -0.50842 -0.39040 0.02023 0.47508 1.24682 0.09906 -0.05951 0.040 -0.45025 -0.31255 0.01858 0.48906 1.33552 0.12324 -0.06481 0.050 -0.38023 -0.23187 0.02029 0.50412 1.67790 0.18762 -0.08741 0.070 -0.35050 -0.18413 0.02376 0.50892 1.57403 0.12994 -0.07910 0.100 -0.42752 -0.37652 0.03221 0.49777 1.52282 0.12604 -0.07408 0.150 -0.55919 -0.53679 0.03248 0.47977 1.31863 0.11085 -0.05612 0.200 -0.66730 -0.65710 0.02956 0.46896 1.21025 0.10065 -0.04777 0.250 -0.73135 -0.69189 0.02516 0.45698 1.13978 0.07837 -0.03958 0.300 -0.78840 -0.68208 0.03152 0.45065 1.05645 0.04621 -0.03245 0.350 -0.83320 -0.69252 0.03233 0.44141 1.01481 0.05533 -0.02765 0.400 -0.86810 -0.74537 0.03521 0.43589 1.00182 0.05914 -0.02363 0.450 -0.88575 -0.73547 0.03923 0.42954 0.94803 0.06557 -0.01790 0.500 -0.89944 -0.69269 0.04159 0.42699 0.94724 0.06067 -0.01710 0.600 -0.91493 -0.63480 0.04580 0.41593 0.95504 0.07576 -0.01606 0.700 -0.93236 -0.63204 0.04993 0.40303 1.01362 0.08323 -0.01527 0.750 -0.93217 -0.63780 0.04989 0.40219 1.03634 0.08203 -0.01622 0.800 -0.92975 -0.65092 0.05114 0.39766 1.05807 0.08385 -0.01434 0.900 -0.92777 -0.57775 0.05266 0.38861 1.11036 0.09388 -0.01658 1.000 -0.93815 -0.60041 0.05421 0.38150 1.16634 0.09095 -0.01502 1.200 -0.93377 -0.56801 0.05576 0.36982 1.29484 0.08078 -0.01434 1.400 -0.93847 -0.48684 0.05782 0.35868 1.32222 0.08353 -0.00681 1.600 -0.92242 -0.40484 0.05645 0.35713 1.30431 0.07158 -0.00268 1.800 -0.91608 -0.29053 0.05615 0.34643 1.35426 0.07341 0.00000 2.000 -0.90369 -0.18149 0.05307 0.34133 1.38763 0.06790 0.00000 2.500 -0.89442 -0.04175 0.05954 0.33960 1.41986 0.08582 0.00000 3.000 -0.87386 0.00000 0.05596 0.35349 1.37795 0.10208 0.00000 3.500 -0.85510 0.00000 0.05469 0.35286 1.34678 0.07501 0.00000 4.000 -0.84680 0.00000 0.05469 0.36845 1.25830 0.05876 0.00000 """) COEFFS_REG = CoeffsTable(sa_damping=5, table="""\ imt ckUSNZ ckJP ckTW ckCH ckWA ckGRTR ckWMT ckNWE pga -0.0302 0.0117 -0.0233 0.01580 0.10010 -0.0118 0.01720 0.0314 0.010 -0.0302 0.0117 -0.0233 0.01580 0.10010 -0.0118 0.01720 0.0314 0.025 -0.0303 0.0135 -0.0272 0.01500 0.10130 -0.0100 0.01740 0.0264 0.040 -0.0336 0.0298 -0.0394 0.01110 0.10590 -0.0148 0.01010 0.0178 0.050 -0.0400 0.0575 -0.0541 0.00990 0.10710 -0.0240 -0.00930 0.0038 0.070 -0.0346 0.0508 -0.0560 -0.00120 0.11190 -0.0190 -0.01140 -0.0206 0.100 -0.0287 0.0199 -0.0450 0.02200 0.12510 -0.0095 0.00840 -0.0222 0.150 -0.0187 -0.0228 -0.0114 0.01430 0.11050 0.0044 0.02580 -0.0307 0.200 -0.0196 -0.0439 0.0089 0.00560 0.11340 0.0133 0.03500 -0.0254 0.250 -0.0227 -0.0543 0.0222 0.00590 0.10160 0.0162 0.04800 0.0274 0.300 -0.0216 -0.0583 0.0300 -0.00003 0.08600 0.0153 0.05800 0.0407 0.350 -0.0187 -0.0583 0.0301 0.00250 0.08900 0.0135 0.05340 0.0650 0.400 -0.0239 -0.0544 0.0313 0.00800 0.09462 0.0070 0.05177 0.0728 0.450 -0.0254 -0.0502 0.0327 0.01420 0.09990 0.0041 0.05190 0.0798 0.500 -0.0322 -0.0461 0.0360 0.01560 0.10730 -0.0022 0.05530 0.0879 0.600 -0.0388 -0.0389 0.0356 0.01630 0.12090 -0.0125 0.05650 0.0978 0.700 -0.0411 -0.0333 0.0336 0.02200 0.12460 -0.0197 0.04830 0.1104 0.750 -0.0416 -0.0305 0.0339 0.02520 0.12240 -0.0269 0.04850 0.1166 0.800 -0.0436 -0.0289 0.0346 0.02970 0.12440 -0.0321 0.05120 0.1193 0.900 -0.0412 -0.0262 0.0289 0.03250 0.12390 -0.0408 0.05740 0.1303 1.000 -0.0397 -0.0195 0.0146 0.03750 0.12730 -0.0434 0.06730 0.1369 1.200 -0.0395 -0.0071 -0.0025 0.04630 0.13760 -0.0467 0.06680 0.0914 1.400 -0.0365 -0.0036 -0.0115 0.05740 0.13970 -0.0446 0.06400 0.0893 1.600 -0.0361 0.0073 -0.0188 0.06200 0.13190 -0.0473 0.06000 0.0914 1.800 -0.0307 0.0108 -0.0252 0.06090 0.13320 -0.0452 0.05230 0.1062 2.000 -0.0280 0.0129 -0.0328 0.05910 0.14080 -0.0445 0.04100 0.1092 2.500 -0.0336 0.0277 -0.0413 0.05880 0.14710 -0.0316 0.01970 0.0509 3.000 -0.0325 0.0369 -0.0579 0.05660 0.16790 -0.0268 0.01380 0.1050 3.500 -0.0272 0.0461 -0.0630 0.05250 0.14220 -0.0294 0.02160 0.1560 4.000 -0.0203 0.0503 -0.0641 0.05720 0.19450 -0.0242 0.01380 0.2198 """)