Source code for openquake.hazardlib.gsim.sera_amplification_models
# -*- 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/>."""Implements SERA site amplification modelsclass:`PitilakisEtAl2018`, `PitilakisEtAl2020`, `Eurocode8Amplification`,`Eurocode8AmplificationDefault`,`SandikkayaDinsever2018`"""importnumpyasnpimportcopyfromscipy.constantsimportgfromopenquake.baselib.generalimportCallableDictfromopenquake.hazardlib.gsim.baseimportGMPE,CoeffsTable,registryfromopenquake.hazardlib.imtimportPGA,SA,from_stringfromopenquake.hazardlibimportconst,contexts# Pitilakis GMPE Wrapperuppernames='''DEFINED_FOR_INTENSITY_MEASURE_TYPESDEFINED_FOR_STANDARD_DEVIATION_TYPESREQUIRES_SITES_PARAMETERSREQUIRES_RUPTURE_PARAMETERSREQUIRES_DISTANCES'''.split()CONSTANTS={"F0":2.5,"kappa":5.0,"TA":0.03}IMLS=[0.,0.25,0.5,0.75,1.,1.25]MEAN,SIGMA,INTER,INTRA=0,1,2,3get_amplification_factor=CallableDict()
[docs]@get_amplification_factor.add("base")defget_amplification_factor_1(kind,F1,FS,s_s,s_1,sctx,ec8=None):""" 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)forec8binnp.unique(sctx.ec8_p18):ec8=ec8b.decode('ascii')ifec8=="A":# Amplification factors are 1continueidx=sctx.ec8_p18==ec8bifnp.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.25ub=s_ss>1.25f_ss[lb]=FS[ec8][0]f_ls[lb]=F1[ec8][0]f_ss[ub]=FS[ec8][-1]f_ls[ub]=F1[ec8][-1]forjinrange(1,len(IMLS)-1):jdx=np.logical_and(s_ss>=IMLS[j],s_ss<IMLS[j+1])ifnotnp.any(jdx):continuedfs=FS[ec8][j+1]-FS[ec8][j]dfl=F1[ec8][j+1]-F1[ec8][j]diml=IMLS[j+1]-IMLS[j]f_ss[jdx]=FS[ec8][j]+(s_ss[jdx]-IMLS[j])*(dfs/diml)f_ls[jdx]=F1[ec8][j]+(s_ss[jdx]-IMLS[j])*(dfl/diml)f_s[idx]=f_ssf_l[idx]=f_lsreturnf_s,f_l
[docs]@get_amplification_factor.add("euro8")defget_amplification_factor_2(kind,F1,FS,s_s_rp,s_1_rp,sctx,ec8):""" 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.fors_cinnp.unique(ec8):ifs_c==b"A":continueidx=ec8==s_cifnotnp.any(idx):continueifs_cin(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])elifs_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.))elifs_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:passreturnf_s,f_l
[docs]@get_amplification_factor.add("euro8default")defget_amplification_factor_3(kind,F1,FS,s_s_rp,s_1_rp,sctx,ec8=None):""" Returns the default amplification factor dependent upon the site class """f_s=np.ones(sctx.ec8.shape)f_l=np.ones(sctx.ec8.shape)forkeyinEC8_FS_default:idx=sctx.ec8==keyf_s[idx]=EC8_FS_default[key]f_l[idx]=EC8_FL_default[key]returnf_s,f_l
[docs]defget_amplified_mean(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"instr(imt)orimt.period<=CONSTANTS["TA"]:# PGA or v. short period accelerationreturnnp.log(s_s/CONSTANTS["F0"])mean=np.copy(s_s)t_c=s_1/s_st_b=t_c/CONSTANTS["kappa"]t_b[t_b<0.05]=0.05t_b[t_b>0.1]=0.1t_d=2.0+np.zeros_like(s_1_rp)idx=s_1_rp>0.1t_d[idx]=1.0+(10.*s_1_rp[idx])idx=np.logical_and(CONSTANTS["TA"]<imt.period,t_b>=imt.period)mean[idx]=(s_s[idx]/(t_b[idx]-CONSTANTS["TA"]))*\
((imt.period-CONSTANTS["TA"])+(t_b[idx]-imt.period)/CONSTANTS["F0"])idx=np.logical_and(t_c<imt.period,t_d>=imt.period)mean[idx]=s_1[idx]/imt.periodidx=t_d<imt.periodmean[idx]=t_d[idx]*s_1[idx]/(imt.period**2.)returnnp.log(mean)
[docs]defget_ec8_class(vsh,h800):""" Method to return the vector of Eurocode 8 site classes based on Vs30 and h800 """ec8=np.array([b"A"foriinrange(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"returnec8
[docs]classPitilakisEtAl2018(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 """kind="base"#: Supported tectonic region type is undefined (applies to any)DEFINED_FOR_TECTONIC_REGION_TYPE=""#: Supported intensity measure types are not setDEFINED_FOR_INTENSITY_MEASURE_TYPES={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 typeDEFINED_FOR_STANDARD_DEVIATION_TYPES={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 laterREQUIRES_RUPTURE_PARAMETERS={'mag'}#: Required distance metrics will be set by the GMPEsREQUIRES_DISTANCES=set()#: Defined reference velocity is 800 m/sDEFINED_FOR_REFERENCE_VELOCITY=800.0def__init__(self,gmpe_name,reference_velocity=None,**extra_args):ifisinstance(gmpe_name,str):self.gsim=registry[gmpe_name](**extra_args)else:# An instantiated class is passed as an argumentself.gsim=copy.deepcopy(gmpe_name)ifreference_velocity:self.rock_vs30=reference_velocityelse:self.rock_vs30=self.DEFINED_FOR_REFERENCE_VELOCITYfornameinuppernames:setattr(self,name,frozenset(getattr(self,name)|getattr(self.gsim,name)))
[docs]defcompute(self,ctx:np.recarray,imts,mean,sig,tau,phi):""" 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. """# Get PGA and Sa (1.0) from GMPEctx_r=copy.copy(ctx)ctx_r.vs30=np.full_like(ctx_r.vs30,self.rock_vs30)rock=contexts.get_mean_stds(self.gsim,ctx_r,[PGA(),SA(1.0)])pga_r=rock[0,0]s_1_rp=rock[0,1]s_s_rp=CONSTANTS["F0"]*np.exp(pga_r)s_1_rp=np.exp(s_1_rp)# Get the short and long period amplification factorsifself.kind=='euro8':ec8=get_ec8_class(ctx.vs30,ctx.h800)else:ec8=Nonef_s,f_l=get_amplification_factor(self.kind,self.F1,self.FS,s_s_rp,s_1_rp,ctx,ec8)s_1=f_l*s_1_rps_s=f_s*s_s_rp# NB: this is wasteful since means are computed and then discardedout=contexts.get_mean_stds(self.gsim,ctx_r,imts)form,imtinenumerate(imts):# Get the mean ground motion using the design code spectrummean[m]=get_amplified_mean(s_s,s_1,s_1_rp,imt)sig[m]=out[1,m]tau[m]=out[2,m]phi[m]=out[3,m]
# 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]classPitilakisEtAl2020(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. """# 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]classEurocode8Amplification(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. """kind="euro8"#: Supported tectonic region type is undefined (applies to any)DEFINED_FOR_TECTONIC_REGION_TYPE=""#: Supported intensity measure types are not setDEFINED_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 typeDEFINED_FOR_STANDARD_DEVIATION_TYPES=set([const.StdDev.TOTAL])#: Required site parameters will be set be selected GMPESREQUIRES_SITES_PARAMETERS={'vs30','h800'}#: Required rupture parameter is magnitude, others will be set laterREQUIRES_RUPTURE_PARAMETERS={'mag'}#: Required distance metrics will be set by the GMPEsREQUIRES_DISTANCES=set()#: Defined reference velocity is 800 m/sDEFINED_FOR_REFERENCE_VELOCITY=800.0def__init__(self,gmpe_name,reference_velocity=800.0,**kwargs):super().__init__(gmpe_name=gmpe_name,**kwargs)self.rock_vs30=reference_velocityifreference_velocityelse\
self.DEFINED_FOR_REFERENCE_VELOCITYfornameinuppernames:setattr(self,name,frozenset(getattr(self,name)|getattr(self.gsim,name)))
# Default short period amplification factors defined by Eurocode 8 Table 3.4EC8_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.4EC8_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]classEurocode8AmplificationDefault(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 """kind="euro8default"#: Required site parameters are the EC8 site class, everything else will#: be set be selected GMPESREQUIRES_SITES_PARAMETERS={'ec8'}
# Sandikkaya & DinseverREGION_SET=["USNZ","JP","TW","CH","WA","TRGR","WMT","NWE"]def_get_basin_term(C,ctx,region=None):""" Get basin amplification term """returnC["b2"]*np.log(ctx.z1pt0)
[docs]defget_site_amplification(C,psarock,ctx,ck):""" Returns the site amplification model define in equation (9) """vs30_s=np.copy(ctx.vs30)vs30_s[vs30_s>1000.]=1000.fn_lin=(C["b1"]+ck)*np.log(vs30_s/760.)fn_z=_get_basin_term(C,ctx)fn_nl=C["b3"]*np.log((psarock+0.1*g)/(0.1*g))*\
np.exp(-np.exp(2.0*np.log(ctx.vs30)-11.))returnfn_lin+fn_z+fn_nl
[docs]defget_stddevs(phi_0,C,tau,phi,psa_rock,vs30,imt):""" Returns the standard deviation adjusted for the site-response model """ysig=np.copy(psa_rock)ysig[ysig>0.35]=0.35ysig[ysig<0.005]=0.005vsig=np.copy(vs30)vsig[vsig>600.0]=600.0vsig[vsig<150.]=150.sigma_s=C["sigma_s"]*C["c0"]*(C["c1"]*np.log(ysig)+C["c2"]*np.log(vsig))ifphi_0:phi0=phi_0[imt]['value']+np.zeros(vs30.shape)else:# In the case that no input phi0 is defined take 'approximate'# phi0 as 85 % of phiphi0=0.85*phiphi=np.sqrt(phi0**2.+sigma_s**2.)return[np.sqrt(tau**2+phi**2),tau,phi]
[docs]classSandikkayaDinsever2018(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 undefinedDEFINED_FOR_TECTONIC_REGION_TYPE=""#: Supported intensity measure types are not setDEFINED_FOR_INTENSITY_MEASURE_TYPES={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 typeDEFINED_FOR_STANDARD_DEVIATION_TYPES={const.StdDev.TOTAL}#: Required site parameters will be set be selected GMPESREQUIRES_SITES_PARAMETERS={'vs30','z1pt0'}#: Required rupture parameter is magnitude, others will be set laterREQUIRES_RUPTURE_PARAMETERS={'mag'}#: Required distance metrics will be set by the GMPEsREQUIRES_DISTANCES=set()#: Defined reference velocity is 800 m/sDEFINED_FOR_REFERENCE_VELOCITY=760.0def__init__(self,gmpe_name,reference_velocity=760.,region=None,phi_0=None,**kwargs):super().__init__(**kwargs)ifisinstance(gmpe_name,str):self.gsim=registry[gmpe_name](**kwargs)else:# An instantiated class is passed as an argumentself.gsim=copy.deepcopy(gmpe_name)# Define the reference velocity - set to 760. by defaultself.rock_vs30=reference_velocityifreference_velocityelse\
self.DEFINED_FOR_REFERENCE_VELOCITYfornameinuppernames:setattr(self,name,frozenset(getattr(self,name)|getattr(self.gsim,name)))stddev_check=(const.StdDev.INTER_EVENTinself.DEFINED_FOR_STANDARD_DEVIATION_TYPES)and\
(const.StdDev.INTRA_EVENTinself.DEFINED_FOR_STANDARD_DEVIATION_TYPES)ifnotstddev_check:raiseValueError("Input GMPE %s not defined for inter- and intra-""event standard deviation"%str(self.gsim))ifisinstance(phi_0,dict):# Input phi_0 modeliphi_0={from_string(key):{'value':phi_0[key]}forkeyinphi_0}self.phi_0=CoeffsTable.fromdict(iphi_0)else:# No input phi_0 modelself.phi_0=None# Regionalisation of the linear site term is possible# check if region is in the set of supported terms and# raise error otherwiseifregionisnotNone:ifregioninREGION_SET:self.region="ck{:s}".format(region)else:raiseValueError("Region must be one of: %s"%" ".join(REGION_SET))else:self.region=region
[docs]defcompute(self,ctx:np.recarray,imts,mean,sig,tau,phi):""" Returns the mean and standard deviations """ctx_r=copy.copy(ctx)ctx_r.vs30=np.full_like(ctx_r.vs30,self.rock_vs30)rock=contexts.get_mean_stds(self.gsim,ctx_r,imts)form,imtinenumerate(imts):psarock=np.exp(rock[MEAN][m])C=self.COEFFS_SITE[imt]ifself.region:ck=self.COEFFS_REG[imt][self.region]else:ck=0.0mean[m]=rock[MEAN][m]+get_site_amplification(C,psarock,ctx,ck)t=rock[INTER][m]p=rock[INTRA][m]sig[m],tau[m],phi[m]=get_stddevs(self.phi_0,C,t,p,psarock,ctx.vs30,imt)
