Source code for openquake.hazardlib.gsim.abrahamson_2014
# -*- coding: utf-8 -*-# vim: tabstop=4 shiftwidth=4 softtabstop=4## Copyright (C) 2014-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:`AbrahamsonEtAl2014` :class:`AbrahamsonEtAl2014RegCHN` :class:`AbrahamsonEtAl2014RegJPN` :class:`AbrahamsonEtAl2014RegTWN`"""importcopyimportnumpyasnpfromscipyimportinterpolatefromopenquake.hazardlib.gsim.baseimportGMPE,CoeffsTable,add_aliasfromopenquake.hazardlibimportconstfromopenquake.hazardlib.imtimportPGA,PGV,SAMETRES_PER_KM=1000.0#: equation constants (that are IMT independent)CONSTS={'n':1.5,# m2 specified at page 1032 (top)'m2':5.00,# h1, h2, h3 specified at page 1040 (top)'h1':+0.25,'h2':+1.50,'h3':-0.75}def_get_phi_al_regional(C,mag,vs30measured,rrup):""" Returns intra-event (Phi) standard deviation (equation 24, page 1046) """s1=np.ones_like(mag)*C['s1e']s2=np.ones_like(mag)*C['s2e']s1[vs30measured]=C['s1m']s2[vs30measured]=C['s2m']phi_al=s1+(s2-s1)/2.*(mag-4.)phi_al[mag<4]=s1[mag<4]phi_al[mag>=6]=s2[mag>=6]returnphi_aldef_get_phi_al_regional_JPN(C,mag,vs30measured,rrup):""" Returns intra-event (Tau) standard deviation (equation 26, page 1046) """phi_al=np.ones((len(vs30measured)))idx=rrup<30phi_al[idx]*=C['s5']idx=(rrup<=80.)&(rrup>=30.)phi_al[idx]*=C['s5']+(C['s6']-C['s5'])/50.*(rrup[idx]-30.)idx=rrup>80phi_al[idx]*=C['s6']returnphi_aldef_get_basic_term(C,ctx):""" Compute and return basic form, see page 1030. """# Fictitious depth calculationc4m=C['c4']-(C['c4']-1.)*(5.-ctx.mag)c4m[ctx.mag>5.]=C['c4']c4m[ctx.mag<4]=1.R=np.sqrt(ctx.rrup**2.+c4m**2.)# basic formbase_term=C['a1']*np.ones_like(ctx.rrup)+C['a17']*ctx.rrup# equation 2 at page 1030after=ctx.mag>=C['m1']within=(ctx.mag>=CONSTS['m2'])&(ctx.mag<C['m1'])before=ctx.mag<CONSTS['m2']base_term[after]+=(C['a5']*(ctx.mag[after]-C['m1'])+C['a8']*(8.5-ctx.mag[after])**2.+(C['a2']+C['a3']*(ctx.mag[after]-C['m1']))*np.log(R[after]))base_term[within]+=(C['a4']*(ctx.mag[within]-C['m1'])+C['a8']*(8.5-ctx.mag[within])**2.+(C['a2']+C['a3']*(ctx.mag[within]-C['m1']))*np.log(R[within]))base_term[before]+=(C['a4']*(CONSTS['m2']-C['m1'])+C['a8']*(8.5-CONSTS['m2'])**2.+C['a6']*(ctx.mag[before]-CONSTS['m2'])+C['a7']*(ctx.mag[before]-CONSTS['m2'])**2.+(C['a2']+C['a3']*(CONSTS['m2']-C['m1']))*np.log(R[before]))returnbase_termdef_get_derivative(C,sa1180,vs30):""" Returns equation 30 page 1047 """derAmp=np.zeros_like(vs30)n=CONSTS['n']c=C['c']b=C['b']idx=vs30<C['vlin']derAmp[idx]=(b*sa1180[idx]*(-1./(sa1180[idx]+c)+1./(sa1180[idx]+c*(vs30[idx]/C['vlin'])**n)))returnderAmpdef_get_faulting_style_term(C,ctx):""" Compute and return faulting style term, that is the sum of the second and third terms in equation 1, page 74. """# this implements equations 5 and 6 at page 1032. f7 is the# coefficient for reverse mechanisms while f8 is the correction# factor for normal rupturesf7=C['a11']*np.clip(ctx.mag-4.,0.,1.)f8=C['a12']*np.clip(ctx.mag-4.,0.,1.)# ranges of rake values for each faulting mechanism are specified in# table 2, page 1031return(f7*((ctx.rake>30)&(ctx.rake<150))+f8*((ctx.rake>-150)&(ctx.rake<-30)))def_get_hanging_wall_term(C,ctx):""" Compute and return hanging wall model term, see page 1038. """Fhw=np.zeros_like(ctx.rx)Fhw[ctx.rx>0]=1.# Taper 1T1=_hw_taper1(ctx)# Taper 2T2=_hw_taper2(ctx)# Taper 3T3=_hw_taper3(ctx)# Taper 4T4=_hw_taper4(ctx)# Taper 5T5=_hw_taper5(ctx)Fhw[ctx.dip==90.]=0.# Finally, compute the hanging wall termreturnFhw*C['a13']*T1*T2*T3*T4*T5def_get_inter_event_std(C,mag,sa1180,vs30):""" Returns inter event (tau) standard deviation (equation 25, page 1046) """tau_al=C['s3']+(C['s4']-C['s3'])/2.*(mag-5.)tau_al[mag<5.]=C['s3']tau_al[mag>=7.]=C['s4']tau_b=tau_altau=tau_b*(1+_get_derivative(C,sa1180,vs30))returntaudef_get_intra_event_std(region,C,mag,sa1180,vs30,vs30measured,rrup):""" Returns Phi as described at pages 1046 and 1047 """phi_al=_get_phi_al_regional_JPN(C,mag,vs30measured,rrup)ifregion=='JPN'else_get_phi_al_regional(C,mag,vs30measured,rrup)derAmp=_get_derivative(C,sa1180,vs30)phi_amp=0.4idx=phi_al<phi_ampifnp.any(idx):# In the case of small magnitudes and long periods it is possible# for phi_al to take a value less than phi_amp, which would return# a complex value. According to the GMPE authors in this case# phi_amp should be reduced such that it is fractionally smaller# than phi_alphi_amp=0.4*np.ones_like(phi_al)phi_amp[idx]=0.99*phi_al[idx]phi_b=np.sqrt(phi_al**2-phi_amp**2)phi=np.sqrt(phi_b**2*(1+derAmp)**2+phi_amp**2)returnphidef_get_regional_term(region,C,imt,vs30,rrup):""" In accordance with Abrahamson et al. (2014) we assume California as the default region. """ifregion=='TWN':vs30star=_get_vs30star(vs30,imt)returnC['a31']*np.log(vs30star/C['vlin'])+C['a25']*rrupelifregion=='CHN':returnC['a28']*rrupelifregion=='JPN':# See page 1043f3=interpolate.interp1d([150,250,350,450,600,850,1150,2000],[C['a36'],C['a37'],C['a38'],C['a39'],C['a40'],C['a41'],C['a42'],C['a42']],fill_value='extrapolate',kind='linear')returnf3(vs30)+C['a29']*rrupelse:# Californiareturn0.def_get_site_response_term(C,imt,vs30,sa1180):""" Compute and return site response model term see page 1033 """# vs30 starvs30_star=_get_vs30star(vs30,imt)# compute the site termsite_resp_term=np.zeros_like(vs30)gt_vlin=vs30>=C['vlin']lw_vlin=vs30<C['vlin']# compute site response term for ctx with vs30 greater than vlinvs30_rat=vs30_star/C['vlin']site_resp_term[gt_vlin]=((C['a10']+C['b']*CONSTS['n'])*np.log(vs30_rat[gt_vlin]))# compute site response term for ctx with vs30 lower than vlinsite_resp_term[lw_vlin]=(C['a10']*np.log(vs30_rat[lw_vlin])-C['b']*np.log(sa1180[lw_vlin]+C['c'])+C['b']*np.log(sa1180[lw_vlin]+C['c']*vs30_rat[lw_vlin]**CONSTS['n']))returnsite_resp_termdef_get_soil_depth_term(region,C,z1pt0,vs30):""" Compute and return soil depth term. See page 1042. """# Get reference z1pt0z1ref=_get_z1pt0ref(region,vs30)# Get z1pt0z10=copy.deepcopy(z1pt0)z10/=METRES_PER_KM# This is used for the calculation of the motion on reference rockidx=z1pt0<0z10[idx]=z1ref[idx]factor=np.log((z10+0.01)/(z1ref+0.01))# Here we use a linear interpolation as suggested in the 'Application# guidelines' at page 1044# Above 700 m/s the trend is flat, but we extend the Vs30 range to# 6,000 m/s (basically the upper limit for mantle shear wave velocity# on earth) to allow extrapolation without throwing an error.f2=interpolate.interp1d([0.0,150,250,400,700,1000,6000],[C['a43'],C['a43'],C['a44'],C['a45'],C['a46'],C['a46'],C['a46']],kind='linear')returnf2(vs30)*factordef_get_stddevs(region,C,imt,ctx,sa1180):""" Return standard deviations as described in paragraph 'Equations for standard deviation', page 1046. """std_intra=_get_intra_event_std(region,C,ctx.mag,sa1180,ctx.vs30,ctx.vs30measured,ctx.rrup)std_inter=_get_inter_event_std(C,ctx.mag,sa1180,ctx.vs30)return[np.sqrt(std_intra**2+std_inter**2),std_inter,std_intra]def_get_top_of_rupture_depth_term(C,imt,ctx):""" Compute and return top of rupture depth term. See paragraph 'Depth-to-Top of Rupture Model', page 1042. """returnC['a15']*np.clip(ctx.ztor/20.0,None,1.)def_get_vs30star(vs30,imt):""" This computes equations 8 and 9 at page 1034 """# compute the v1 value (see eq. 9, page 1034)ifimt.string[:2]=="SA":t=imt.periodift<=0.50:v1=1500.0elift<3.0:v1=np.exp(-0.35*np.log(t/0.5)+np.log(1500.))else:v1=800.0elifimt.string=="PGA":v1=1500.0else:# This covers the PGV casev1=1500.0# set the vs30 star value (see eq. 8, page 1034)vs30_star=np.ones_like(vs30)*vs30vs30_star[vs30>=v1]=v1returnvs30_stardef_get_z1pt0ref(region,vs30):""" This computes the reference depth to the 1.0 km/s interface using equation 18 at page 1042 of Abrahamson et al. (2014) """ifregion=='JPN':return1./1000.*np.exp(-5.23/2.*np.log((vs30**2+412.**2.)/(1360.**2+412**2.)))return(1./1000.)*np.exp((-7.67/4.)*np.log((vs30**4+610.**4)/(1360.**4+610.**4)))def_hw_taper1(ctx):# Compute taper t1T1=(90.-ctx.dip)/45.0T1[ctx.dip<=30.]=60./45.returnT1def_hw_taper2(ctx):# Compute taper t2 (eq 12 at page 1039) - a2hw set to 0.2 as# indicated at page 1041a2hw=0.2T2=(1.+a2hw*(ctx.mag-6.5)-(1.-a2hw)*(ctx.mag-6.5)**2)T2[ctx.mag>6.5]=1.+a2hw*(ctx.mag[ctx.mag>6.5]-6.5)T2[ctx.mag<=5.5]=0.returnT2def_hw_taper3(ctx):# Compute taper t3 (eq. 13 at page 1039) - r1 and r2 specified at# page 1040T3=np.zeros_like(ctx.rx)r1=ctx.width*np.cos(np.radians(ctx.dip))r2=3.*r1idx=ctx.rx<r1T3[idx]=(np.ones_like(ctx.rx)[idx]*CONSTS['h1']+CONSTS['h2']*(ctx.rx[idx]/r1[idx])+CONSTS['h3']*(ctx.rx[idx]/r1[idx])**2)idx=(ctx.rx>=r1)&(ctx.rx<=r2)T3[idx]=1.-(ctx.rx[idx]-r1[idx])/(r2[idx]-r1[idx])returnT3def_hw_taper4(ctx):# Compute taper t4 (eq. 14 at page 1040)T4=np.zeros_like(ctx.rx)T4[ctx.ztor<=10.]=1.-ctx.ztor[ctx.ztor<=10.]**2./100.returnT4def_hw_taper5(ctx):# Compute T5 (eq 15a at page 1040) - ry1 computed according to# suggestions provided at page 1040T5=np.zeros_like(ctx.rx)ry1=ctx.rx*np.tan(np.radians(20.))idx=(ctx.ry0-ry1)<=0.0T5[idx]=1.idx=((ctx.ry0-ry1)>0.0)&((ctx.ry0-ry1)<5.0)T5[idx]=1.-(ctx.ry0[idx]-ry1[idx])/5.0returnT5def_get_sa_at_1180(region,C,imt,ctx):""" Compute and return mean imt value for rock conditions (vs30 = 1100 m/s) """# reference vs30 = 1180 m/svs30_1180=np.ones_like(ctx.vs30)*1180.# reference shaking intensity = 0ref_iml=np.zeros_like(ctx.vs30)# fake Z1.0 - Since negative it will be replaced by the default Z1.0# for the corresponding regionfake_z1pt0=np.ones_like(ctx.vs30)*-1return(_get_basic_term(C,ctx)+_get_faulting_style_term(C,ctx)+_get_site_response_term(C,imt,vs30_1180,ref_iml)+_get_hanging_wall_term(C,ctx)+_get_top_of_rupture_depth_term(C,imt,ctx)+_get_soil_depth_term(region,C,fake_z1pt0,vs30_1180)+_get_regional_term(region,C,imt,vs30_1180,ctx.rrup))
[docs]defget_epistemic_sigma(ctx):""" This function gives the epistemic sigma computed following USGS-2014 approach. Also, note that the events are counted in each magnitude and distance bins. However, the epistemic sigma is based on NZ SMDB v1.0 """n=2dist_func_5_6=np.where(ctx.rrup<=10,0.4*np.sqrt(n/11),np.where((ctx.rrup>10)&(ctx.rrup<30),0.4*np.sqrt(n/38),0.4*np.sqrt(n/94)))dist_func_6_7=np.where(ctx.rrup<=10,0.4*np.sqrt(n/2),np.where((ctx.rrup>10)&(ctx.rrup<30),0.4*np.sqrt(n/7),0.4*np.sqrt(n/13)))dist_func_7_above=np.where(ctx.rrup<=10,0.4*np.sqrt(n/2),np.where((ctx.rrup>10)&(ctx.rrup<30),0.4*np.sqrt(n/2),0.4*np.sqrt(n/4)))sigma_epi=np.where((ctx.mag>=5)&(ctx.mag<6),dist_func_5_6,np.where((ctx.mag>=6)&(ctx.mag<7),dist_func_6_7,dist_func_7_above))returnsigma_epi
[docs]classAbrahamsonEtAl2014(GMPE):""" Implements GMPE by Abrahamson, Silva and Kamai developed within the the PEER West 2 Project. This GMPE is described in a paper published in 2014 on Earthquake Spectra, Volume 30, Number 3 and titled 'Summary of the ASK14 Ground Motion Relation for Active Crustal Regions'. """#: Supported tectonic region type is active shallow crust, see title!DEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.ACTIVE_SHALLOW_CRUST#: Supported intensity measure types are spectral acceleration, peak#: ground velocity and peak ground acceleration, see tables 4#: pages 1036DEFINED_FOR_INTENSITY_MEASURE_TYPES={PGA,PGV,SA}#: Supported intensity measure component is orientation-independent#: average horizontal :attr:`~openquake.hazardlib.const.IMC.RotD50`,#: see page 1025.DEFINED_FOR_INTENSITY_MEASURE_COMPONENT=const.IMC.RotD50#: Supported standard deviation types are inter-event, intra-event#: and total, see paragraph "Equations for standard deviations", page#: 1046.DEFINED_FOR_STANDARD_DEVIATION_TYPES={const.StdDev.TOTAL,const.StdDev.INTER_EVENT,const.StdDev.INTRA_EVENT}#: Required site parameters are Vs30 and Z1.0, see table 2, page 1031#: Unit of measure for Z1.0 is [m]REQUIRES_SITES_PARAMETERS={'vs30','z1pt0','vs30measured'}#: Required rupture parameters are magnitude, rake, dip, ztor, and width#: (see table 2, page 1031)REQUIRES_RUPTURE_PARAMETERS={'mag','rake','dip','ztor','width'}#: Required distance measures are Rrup, Rjb, Ry0 and Rx (see Table 2,#: page 1031).REQUIRES_DISTANCES={'rrup','rjb','rx','ry0'}#: Reference rock conditions as defined at pageDEFINED_FOR_REFERENCE_VELOCITY=1180def__init__(self,sigma_mu_epsilon=0.0,region=None):self.region=regionassertself.regionin(None,'CHN','JPN','TWN'),regionself.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. """form,imtinenumerate(imts):C=self.COEFFS[imt]# compute median sa on rock (vs30=1180m/s). Used for site response# term calculationsa1180=np.exp(_get_sa_at_1180(self.region,C,imt,ctx))# For debugging purposes# f1 = _get_basic_term(C, ctx)# f4 = _get_hanging_wall_term(C, ctx)# f5 = _get_site_response_term(C, imt, ctx.vs30, sa1180)# f6 = _get_top_of_rupture_depth_term(C, imt, ctx)# f7 = _get_faulting_style_term(C, ctx)# f10 =_get_soil_depth_term(self.region, C, ctx.z1pt0, ctx.vs30)# fre = _get_regional_term(self.region, C, imt, ctx.vs30, ctx.rrup)# get the mean valuemean[m]=(_get_basic_term(C,ctx)+_get_hanging_wall_term(C,ctx)+_get_site_response_term(C,imt,ctx.vs30,sa1180)+_get_top_of_rupture_depth_term(C,imt,ctx)+_get_faulting_style_term(C,ctx)+_get_soil_depth_term(self.region,C,ctx.z1pt0,ctx.vs30))mean[m]+=_get_regional_term(self.region,C,imt,ctx.vs30,ctx.rrup)mean[m]+=(self.sigma_mu_epsilon*get_epistemic_sigma(ctx))# get standard deviationssig[m],tau[m],phi[m]=_get_stddevs(self.region,C,imt,ctx,sa1180)
