Source code for openquake.hazardlib.gsim.chiou_youngs_2014
# -*- coding: utf-8 -*-# vim: tabstop=4 shiftwidth=4 softtabstop=4## Copyright (C) 2012-2023 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:`ChiouYoungs2014` :class:`ChiouYoungs2014Japan` :class:`ChiouYoungs2014Italy` :class:`ChiouYoungs2014Wenchuan` :class:`ChiouYoungs2014PEER` :class:`ChiouYoungs2014NearFaultEffect`"""importnumpyasnpfromopenquake.baselib.generalimportCallableDictfromopenquake.hazardlib.gsim.baseimportGMPE,CoeffsTablefromopenquake.hazardlib.gsim.abrahamson_2014importget_epistemic_sigmafromopenquake.hazardlibimportconstfromopenquake.hazardlib.imtimportPGA,PGV,SACONSTANTS={"c2":1.06,"c4":-2.1,"c4a":-0.5,"crb":50.0,"c8a":0.2695,"c11":0.0,"phi6":300.0,"phi6jp":800.0}def_get_centered_cdpp(clsname,ctx):""" Returns the centred dpp term (zero by default) """ifclsname.endswith("NearFaultEffect"):returnctx.rcdppreturnnp.zeros(ctx.rrup.shape)def_get_centered_z1pt0(clsname,ctx):""" Get z1pt0 centered on the Vs30- dependent average z1pt0(m) California and non-Japan regions """ifclsname.endswith("Japan"):mean_z1pt0=(-5.23/2.)*np.log(((ctx.vs30**2.)+412.39**2.)/(1360**2.+412.39**2.))returnctx.z1pt0-np.exp(mean_z1pt0)#: California and non-Japan regionsmean_z1pt0=(-7.15/4.)*np.log(((ctx.vs30)**4.+570.94**4.)/(1360**4.+570.94**4.))returnctx.z1pt0-np.exp(mean_z1pt0)def_get_centered_ztor(ctx):""" Get ztor centered on the M- dependent avarage ztor(km) by different fault types. """# Strike-slip and normal faultingmean_ztor=np.clip(2.673-1.136*np.clip(ctx.mag-4.970,0.,None),0.,None)**2# Reverse and reverse-oblique faultingrev=(30.<=ctx.rake)&(ctx.rake<=150.)mean_ztor[rev]=np.clip(2.704-1.226*np.clip(ctx.mag[rev]-5.849,0.0,None),0.,None)**2returnctx.ztor-mean_ztordef_get_ln_y_ref(ctx,C):""" Get an intensity on a reference soil. Implements eq. 13a. """# Reverse faulting flagFrv=1.if30<=ctx.rake<=150else0.# Normal faulting flagFnm=1.if-120<=ctx.rake<=-60else0.# A part in eq. 11mag_test1=np.cosh(2.*np.clip(ctx.mag-4.5,0.,None))# Centered DPPcentered_dpp=0# Centered Ztorcentered_ztor=0dist_taper=np.fmax(1-(np.fmax(ctx.rrup-40,np.zeros_like(ctx))/30.),np.zeros_like(ctx))dist_taper=dist_taper.astype(np.float64)ln_y_ref=(# first part of eq. 11C['c1']+(C['c1a']+C['c1c']/mag_test1)*Frv+(C['c1b']+C['c1d']/mag_test1)*Fnm+(C['c7']+C['c7b']/mag_test1)*centered_ztor+(C['c11']+C['c11b']/mag_test1)*np.cos(np.radians(ctx.dip))**2# second part+C['c2']*(ctx.mag-6)+((C['c2']-C['c3'])/C['cn'])*np.log(1+np.exp(C['cn']*(C['cm']-ctx.mag)))# third part+C['c4']*np.log(ctx.rrup+C['c5']*np.cosh(C['c6']*np.clip(ctx.mag-C['chm'],0,None)))+(C['c4a']-C['c4'])*np.log(np.sqrt(ctx.rrup**2+C['crb']**2))# forth part+(C['cg1']+C['cg2']/(np.cosh(np.clip(ctx.mag-C['cg3'],0,None))))*ctx.rrup# fifth part+C['c8']*dist_taper*np.clip((ctx.mag-5.5,0)/0.8,0.,1.)*np.exp(-1*C['c8a']*(ctx.mag-C['c8b'])**2)*centered_dpp# sixth part# + C['c9'] * Fhw * np.cos(math.radians(ctx.dip)) *# (C['c9a'] + (1 - C['c9a']) * np.tanh(ctx.rx / C['c9b']))# * (1 - np.sqrt(ctx.rjb ** 2 + ctx.ztor ** 2)# / (ctx.rrup + 1.0)))returnln_y_refdef_get_mean(ctx,C,ln_y_ref,exp1,exp2):""" Add site effects to an intensity. Implements eq. 13b. """eta=epsilon=0.ln_y=(# first line of eq. 12ln_y_ref+eta# second line+C['phi1']*np.log(ctx.vs30/1130).clip(-np.inf,0)# third line+C['phi2']*(exp1-exp2)*np.log((np.exp(ln_y_ref)*np.exp(eta)+C['phi4'])/C['phi4'])# fourth line - removed# fifth line+epsilon)returnln_y
[docs]defget_basin_depth_term(clsname,C,centered_z1pt0):""" Returns the basin depth scaling """ifclsname.endswith("Japan"):returnC["phi5jp"]*(1.0-np.exp(-centered_z1pt0/CONSTANTS["phi6jp"]))returnC["phi5"]*(1.0-np.exp(-centered_z1pt0/CONSTANTS["phi6"]))
[docs]defget_directivity(clsname,C,ctx):""" Returns the directivity term. The directivity prediction parameter is centered on the average directivity prediction parameter. Here we set the centered_dpp equal to zero, since the near fault directivity effect prediction is off by default in our calculation. """cdpp=_get_centered_cdpp(clsname,ctx)ifnotnp.any(cdpp>0.0):# No directivity termreturn0.0f_dir=np.exp(-C["c8a"]*((ctx.mag-C["c8b"])**2.))*cdppf_dir*=np.clip((ctx.mag-5.5)/0.8,0.,1.)rrup_max=ctx.rrup-40.rrup_max[rrup_max<0.0]=0.0rrup_max=1.0-(rrup_max/30.)rrup_max[rrup_max<0.0]=0.0returnC["c8"]*rrup_max*f_dir
get_far_field_distance_scaling=CallableDict()
[docs]@get_far_field_distance_scaling.add("CAL")defget_far_field_distance_scaling_1(region,C,mag,rrup):""" Returns the far-field distance scaling term - both magnitude and distance - for California and other regions """# Get the attenuation distance scalingf_r=(CONSTANTS["c4a"]-CONSTANTS["c4"])*np.log(np.sqrt(rrup**2.+CONSTANTS["crb"]**2.))# Get the magnitude dependent termf_rm=C["cg1"]+C["cg2"]/np.cosh(np.clip(mag-C["cg3"],0.0,None))returnf_r+f_rm*rrup
[docs]@get_far_field_distance_scaling.add("JPN")defget_far_field_distance_scaling_2(region,C,mag,rrup):""" Returns the far-field distance scaling term - both magnitude and distance - for Japan """# Get the attenuation distance scalingf_r=(CONSTANTS["c4a"]-CONSTANTS["c4"])*np.log(np.sqrt(rrup**2.+CONSTANTS["crb"]**2.))# Get the magnitude dependent termf_rm=(C["cg1"]+C["cg2"]/np.cosh(np.clip(mag-C["cg3"],0.0,None)))*rrup# Apply adjustment factor for Japanf_rm[(mag>6.0)&(mag<6.9)]*=C["gjpit"]returnf_r+f_rm
[docs]@get_far_field_distance_scaling.add("ITA")defget_far_field_distance_scaling_3(region,C,mag,rrup):""" Returns the far-field distance scaling term - both magnitude and distance - for Italy """# Get the attenuation distance scalingf_r=(CONSTANTS["c4a"]-CONSTANTS["c4"])*np.log(np.sqrt(rrup**2.+CONSTANTS["crb"]**2.))# Get the magnitude dependent termf_rm=(C["cg1"]+C["cg2"]/np.cosh(np.clip(mag-C["cg3"],0.0,None)))*rrup# Apply adjustment factor for Italyf_rm[(mag>6.0)&(mag<6.9)]*=C["gjpit"]returnf_r+f_rm
[docs]@get_far_field_distance_scaling.add("WEN")defget_far_field_distance_scaling_4(region,C,mag,rrup):""" Returns the far-field distance scaling term - both magnitude and distance - for Wenchuan """# Get the attenuation distance scalingf_r=(CONSTANTS["c4a"]-CONSTANTS["c4"])*np.log(np.sqrt(rrup**2.+CONSTANTS["crb"]**2.))# Get the magnitude dependent termf_rm=(C["cg1"]+C["cg2"]/np.cosh(np.clip(mag-C["cg3"],0.0,None)))*rrup# Apply adjustment factor for Wenchuanreturnf_r+(f_rm*C["gwn"])
[docs]defget_geometric_spreading(C,mag,rrup):""" Returns the near-field geometric spreading term """# Get the near-field magnitude scalingreturnCONSTANTS["c4"]*np.log(rrup+C["c5"]*np.cosh(C["c6"]*np.clip(mag-C["chm"],0.0,None)))
[docs]defget_hanging_wall_term(C,ctx):""" Returns the hanging wall term """fhw=np.zeros(ctx.rrup.shape)idx=ctx.rx>=0.0ifnp.any(idx):fdist=1.0-(np.sqrt(ctx.rjb[idx]**2.+ctx.ztor[idx]**2.)/(ctx.rrup[idx]+1.0))fdist*=C["c9a"]+(1.0-C["c9a"])*np.tanh(ctx.rx[idx]/C["c9b"])fhw[idx]+=C["c9"]*np.cos(np.radians(ctx.dip[idx]))*fdistreturnfhw
[docs]defget_linear_site_term(clsname,C,ctx):""" Returns the linear site scaling term """ifclsname.endswith("Japan"):returnC["phi1jp"]*np.log(ctx.vs30/1130).clip(-np.inf,0.0)returnC["phi1"]*np.log(ctx.vs30/1130).clip(-np.inf,0.0)
[docs]defget_ln_y_ref(clsname,C,ctx):""" Returns the ground motion on the reference rock, described fully by Equation 11 """region=get_region(clsname)delta_ztor=_get_centered_ztor(ctx)return(get_stress_scaling(C)+get_magnitude_scaling(C,ctx.mag)+get_source_scaling_terms(C,ctx,delta_ztor)+get_hanging_wall_term(C,ctx)+get_geometric_spreading(C,ctx.mag,ctx.rrup)+get_far_field_distance_scaling(region,C,ctx.mag,ctx.rrup)+get_directivity(clsname,C,ctx))
[docs]defget_magnitude_scaling(C,mag):""" Returns the magnitude scaling """f_m=np.log(1.0+np.exp(C["cn"]*(C["cm"]-mag)))f_m=CONSTANTS["c2"]*(mag-6.0)+\
((CONSTANTS["c2"]-C["c3"])/C["cn"])*f_mreturnf_m
[docs]defget_nonlinear_site_term(C,ctx,y_ref):""" Returns the nonlinear site term and the Vs-scaling factor (to be used in the standard deviation model """vs=ctx.vs30.clip(-np.inf,1130.0)f_nl_scaling=C["phi2"]*(np.exp(C["phi3"]*(vs-360.))-np.exp(C["phi3"]*(1130.-360.)))f_nl=np.log((y_ref+C["phi4"])/C["phi4"])*f_nl_scalingreturnf_nl,f_nl_scaling
[docs]defget_phi(C,mag,ctx,nl0):""" Returns the within-event variability described in equation 13, line 3 """phi=C["sig3"]*np.ones(ctx.vs30.shape)phi[ctx.vs30measured]=0.7phi=np.sqrt(phi+((1.0+nl0)**2.))mdep=C["sig1"]+(C["sig2"]-C["sig1"])*np.clip(mag-5.,0.,1.5)/1.5returnmdep*phi
[docs]defget_source_scaling_terms(C,ctx,delta_ztor):""" Returns additional source scaling parameters related to style of faulting, dip and top of rupture depth """f_src=np.zeros_like(ctx.mag)coshm=np.cosh(2.0*np.clip(ctx.mag-4.5,0.,None))# Style of faulting termpos=(30<=ctx.rake)&(ctx.rake<=150)neg=(-120<=ctx.rake)&(ctx.rake<=-60)# reverse faulting flagf_src[pos]+=C["c1a"]+(C["c1c"]/coshm[pos])# normal faulting flagf_src[neg]+=C["c1b"]+(C["c1d"]/coshm[neg])# Top of rupture termf_src+=(C["c7"]+(C["c7b"]/coshm))*delta_ztor# Dip termf_src+=((CONSTANTS["c11"]+(C["c11b"]/coshm))*np.cos(np.radians(ctx.dip))**2.0)returnf_src
[docs]defget_stddevs(clsname,C,ctx,mag,y_ref,f_nl_scaling):""" Returns the standard deviation model described in equation 13 """ifclsname=='ChiouYoungs2014PEER':# the standard deviation, which is fixed at 0.65 for every sitereturn[0.65*np.ones_like(ctx.vs30),0,0]# Determines the nonlinear term described in equation 13, line 4nl0=f_nl_scaling*(y_ref/(y_ref+C["phi4"]))# Get between and within-event variabilitytau=get_tau(C,mag)phi_nl0=get_phi(C,mag,ctx,nl0)# Get total standard deviation propagating the uncertainty in the# nonlinear amplification termsigma=np.sqrt(((1.0+nl0)**2.)*(tau**2.)+phi_nl0**2.)return[sigma,np.abs((1+nl0)*tau),phi_nl0]
[docs]defget_stress_scaling(C):""" Returns the stress drop scaling factor """returnC["c1"]
[docs]defget_tau(C,mag):""" Returns the between-event variability described in equation 13, line 2 """# eq. 13 to calculate inter-event standard errormag_test=np.clip(mag-5.0,0.,1.5)returnC['tau1']+(C['tau2']-C['tau1'])/1.5*mag_test
[docs]defget_mean_stddevs(name,C,ctx):""" Return mean and standard deviation values """# Get ground motion on reference rockln_y_ref=get_ln_y_ref(name,C,ctx)y_ref=np.exp(ln_y_ref)# Get the site amplification# Get basin depthdz1pt0=_get_centered_z1pt0(name,ctx)# for Z1.0 = 0.0 no deep soil correction is applieddz1pt0[ctx.z1pt0<=0.0]=0.0f_z1pt0=get_basin_depth_term(name,C,dz1pt0)# Get linear amplification termf_lin=get_linear_site_term(name,C,ctx)# Get nonlinear amplification termf_nl,f_nl_scaling=get_nonlinear_site_term(C,ctx,y_ref)# Add on the site amplificationmean=ln_y_ref+(f_lin+f_nl+f_z1pt0)# Get standard deviationssig,tau,phi=get_stddevs(name,C,ctx,ctx.mag,y_ref,f_nl_scaling)returnmean,sig,tau,phi
[docs]classChiouYoungs2014(GMPE):""" Implements GMPE developed by Brian S.-J. Chiou and Robert R. Youngs Chiou, B. S.-J. and Youngs, R. R. (2014), "Updated of the Chiou and Youngs NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra, Earthquake Spectra, 30(3), 1117 - 1153, DOI: 10.1193/072813EQS219M """adapted=False# overridden in acme_2019#: Supported tectonic region type is active shallow crustDEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.ACTIVE_SHALLOW_CRUST#: Supported intensity measure types are spectral acceleration,#: peak ground velocity and peak ground accelerationDEFINED_FOR_INTENSITY_MEASURE_TYPES={PGA,PGV,SA}#: Supported intensity measure component is orientation-independent#: measure :attr:`~openquake.hazardlib.const.IMC.RotD50`,DEFINED_FOR_INTENSITY_MEASURE_COMPONENT=const.IMC.RotD50#: Supported standard deviation types are inter-event, intra-event#: and total, see chapter "Variance model".DEFINED_FOR_STANDARD_DEVIATION_TYPES={const.StdDev.TOTAL,const.StdDev.INTER_EVENT,const.StdDev.INTRA_EVENT}#: Required site parameters are Vs30, Vs30 measured flag#: and Z1.0.REQUIRES_SITES_PARAMETERS={'vs30','vs30measured','z1pt0'}#: Required rupture parameters are magnitude, rake,#: dip and ztor.REQUIRES_RUPTURE_PARAMETERS={'dip','rake','mag','ztor'}#: Required distance measures are RRup, Rjb and Rx.REQUIRES_DISTANCES={'rrup','rjb','rx'}#: Reference shear wave velocityDEFINED_FOR_REFERENCE_VELOCITY=1130def__init__(self,sigma_mu_epsilon=0.0,**kwargs):super().__init__(sigma_mu_epsilon=sigma_mu_epsilon,**kwargs)self.sigma_mu_epsilon=sigma_mu_epsilon
[docs]defcompute(self,ctx:np.recarray,imts,mean,sig,tau,phi):""" See :meth:`superclass method <.base.GroundShakingIntensityModel.compute>` for spec of input and result values. """name=self.__class__.__name__# reference to page 1144, PSA might need PGA valuepga_mean,pga_sig,pga_tau,pga_phi=get_mean_stddevs(name,self.COEFFS[PGA()],ctx)form,imtinenumerate(imts):ifrepr(imt)=="PGA":mean[m]=pga_meanmean[m]+=(self.sigma_mu_epsilon*get_epistemic_sigma(ctx))sig[m],tau[m],phi[m]=pga_sig,pga_tau,pga_phielse:imt_mean,imt_sig,imt_tau,imt_phi= \
get_mean_stddevs(name,self.COEFFS[imt],ctx)# reference to page 1144# Predicted PSA value at T ≤ 0.3s should be set equal to the value of PGA# when it falls below the predicted PGAmean[m]=np.where(imt_mean<pga_mean,pga_mean,imt_mean) \
ifrepr(imt).startswith("SA")andimt.period<=0.3 \
elseimt_meanmean[m]+=(self.sigma_mu_epsilon*get_epistemic_sigma(ctx))sig[m],tau[m],phi[m]=imt_sig,imt_tau,imt_phi
[docs]classChiouYoungs2014Japan(ChiouYoungs2014):""" Regionalisation of the Chiou & Youngs (2014) GMPE for use with the Japan far-field distance attuation scaling and site model """
[docs]classChiouYoungs2014Italy(ChiouYoungs2014):""" Adaption of the Chiou & Youngs (2014) GMPE for the the Italy far-field attenuation scaling, but assuming the California site amplification model """
[docs]classChiouYoungs2014Wenchuan(ChiouYoungs2014):""" Adaption of the Chiou & Youngs (2014) GMPE for the Wenchuan far-field attenuation scaling, but assuming the California site amplification model. It should be note that according to Chiou & Youngs (2014) this adjustment is calibrated only for the M7.9 Wenchuan earthquake, so application to other scenarios is at the user's own risk """
[docs]classChiouYoungs2014PEER(ChiouYoungs2014):""" This implements the Chiou & Youngs (2014) GMPE for use with the PEER tests. In this version the total standard deviation is fixed at 0.65 """#: Only the total standars deviation is definedDEFINED_FOR_STANDARD_DEVIATION_TYPES={const.StdDev.TOTAL}#: The PEER tests requires only PGADEFINED_FOR_INTENSITY_MEASURE_TYPES={PGA}
[docs]classChiouYoungs2014NearFaultEffect(ChiouYoungs2014):""" This implements the Chiou & Youngs (2014) GMPE include the near fault effect prediction. In this version, we add the distance measure, rcdpp for directivity prediction. """#: Required distance measures are RRup, Rjb, Rx, and RcdppREQUIRES_DISTANCES={'rrup','rjb','rx','rcdpp'}
[docs]classChiouYoungs2014ACME2019(ChiouYoungs2014):""" Implements a modified version of the CY2014 GMM. Main changes: - Hanging wall term excluded - Centered Ztor = 0 - Centered Dpp = 0 """adapted=True
[docs]defcompute(self,ctx:np.recarray,imts,mean,sig,tau,phi):""" See :meth:`superclass method <.base.GroundShakingIntensityModel.compute>` for spec of input and result values. """form,imtinenumerate(imts):C=self.COEFFS[imt]# intensity on a reference soil is used for both mean# and stddev calculations.ln_y_ref=_get_ln_y_ref(self.__class__.__name__,ctx,C)# exp1 and exp2 are parts of eq. 12 and eq. 13,# calculate it once for both.exp1=np.exp(C['phi3']*(ctx.vs30.clip(-np.inf,1130)-360))exp2=np.exp(C['phi3']*(1130-360))mean[m]=_get_mean(ctx,C,ln_y_ref,exp1,exp2)