设X₁,X₂,…,X_n为地震震级随机变量列,假设它们相互独立且服从同一分布,记其最大震级为M_n=max(X₁,X₂,…,X_n),若存在{a_n>0,b_n∈R}和非退化分布函数G(x),使得:

Pr{(M_n-b_n)/(a_n)≤x}→G(x),n→∞(1)

则G称为极值分布。Fisher和Tippett(1928)证明了极值分布的3大类型定理,指出G必属于Gumbel分布、Frechet分布、Weibull分布之一。广义极值理论将这3种分布类型统一为1个分布类型,由位置参数μ、尺度参数σ、形状参数ξ来表示(Coles,2001; 史道济,2006),如式(2)所示:

G(x)=exp{-[1+ξ((x-μ)/σ)]^-1/ξ},

(1+(ξ(x-μ))/σ>0)(2)

选定震级阈值u,若X_i>u,称Y=X_i-u为震级超出量。则超出量的分布函数可表示为:

F_u(y)=Pr{X-u≤Y〖JB<1|〗X>u}=

(F(u+y)-F(u))/(1-F(u))(y>0)(3)

对于充分大的震级阈值u,震级超出量Y的分布函数近似服从广义帕累托分布:

H(y)=1-(1+(ξy)/(σ^～))^-1/ξ(4)

式中:y>0,(1+(ξy)/(σ^～))>0, 且σ^～=σ+ξ(u-μ)。

广义帕累托分布与广义极值分布具有相同的形状参数ξ。当ξ≥0时,广义帕累托分布均没有上边界; 只有当ξ<0时,广义帕累托分布具有有限上界。

震级上限是指潜在震源区内可能发生的最大地震的震级,可以认为在潜在震源区内未来发生超过该震级地震的概率几乎为零(胡聿贤,1999)。当选定震级阈值u,且X>u时,则有:

Pr{X>x〖JB<1|〗X>u}=[1+ξ((x-u)/(σ^～))]^-1/ξ(5)

令ζ_u=Pr(X>u),表示为超过阈值的样本个数与总样本的比值,式(5)可写为:

Pr{X>x}=ζ_u[1+ξ((x-u)/(σ^～))]^-1/ξ(6)

考虑到潜在震源区的地震事件的震级必为有限值,因此只有当形状参数ξ<0时,广义帕累托分布具有有限上界,才符合震级存在上限的条件,震级上限可表示为:

x=u-σ^～/ξ(7)

利用最大似然估计方法可得到3个参数(ζ_u,σ^～,ξ)的估计值。根据Delta定理可知,3个参数的最大似然估计值近似符合正态分布,因此根据方差协方差矩阵可以获得置信水平为1-α的震级上限的置信区间如式(8)所示:

x^{^}_N±Z_1-α/2(Var(x^{^}_N))^1/2(8)

式中:Var(x^{^}_N)≈ x_N^TV x_N; V是(ζ^{^}_u,σ^{^},ξ^{^})参数估计的协方差矩阵; x_N^T为((x_N)/(ζ_u),(x_N)/(σ^～),(x_N)/(ξ))在(ζ^{^}_u,σ^{^/～},ξ^{^})处的值; Z_1-α/2为标准正态分布1-α/2的分位数(Coles,2001)。

本文数据来源主要为国家地震科学数据共享中心(http://data.earthquake.cn)的中国历史地震目录数据库以及国家地震台网地震目录数据库,黄玮琼等(1994)对中国大陆大部分地区的历史地震资料作了较详细的分析与研究,显示龙门山区域M_S≥43/4地震资料在1931年之后基本完整。因此,笔者收集了1930—2016年龙门山地区(101°～106°E,29.5°～33.5°N)M≥4.5的历史地震目录。经统计,1930年至今,龙门山地区发生过5次M≥7.0强震,其中在龙门山断裂带上发生过2次,分别是2008年汶川8.0级地震和2013年芦山7.0级地震; 北部岷山地区发生过3次,分别是1933年叠溪7.5级地震和1976年松潘2次7.2级地震。

在利用广义帕累托模型进行地震危险性估计之前需要尽量删除地震目录中的余震。目前国内外地震危险性分析研究中,考虑到计算的简便和工程应用上的可操作性,普遍采用Gardner和Knopoff(1974)提出的时空窗法来删除余震(徐伟进,2012)。在该方法中,空间窗口根据式(9)确定(Knopoff,Gardner,1972):

lgR=aM+b(9)

式中:a,b为固定参数值,在我国a,b的取值分别为0.5、-1.78(刘杰等,1996)。

时间窗口如表1所示。删除余震后,共收集到龙门山地震构造区内地震事件118个,对收集到的地震资料统计发现,龙门山地震构造区内以浅源地震为主,其中M≥7地震5次,6≤M<7地震14次,5≤M<6地震34次。从图1可看出,地震震中分布位置与断裂带的展布形态密切相关。

表1 不同震级地震的余震时间窗(Gardner,Knopoff,1974)
Tab.1 Aftershock time windows of different magnitude earthquakes(Garder,Knopoff,1974)

利用广义帕累托分布估计龙门山地区震级上限,首先需确定震级阈值u。阈值选取通过2种方法:(1)依据震级平均超出量分布函数图,震级平均超出量与震级阈值应满足线性关系;(2)随着阈值选取不断增大,形状参数ξ和修饰的尺度参数σ^～估计值应基本保持稳定不变(Coles,2001)。由图2可以看出,当震级阈值在(5.1,5.5)、(5.7,6.0)以及(6.3,6.5)范围内时,平均超出量近似呈线性分布。为了对阈值进行最合理的选取,还需要综合考虑样本数据的实际情况,阈值选取过大,超出量样本数过少,误差增大。结合图3分析,阈值在(5.1,5.5)范围内相应的形状参数和修饰的尺度参数的估计值基本保持不变,根据广义帕累托分布应用条件,选取尽可能大的阈值u=5.5。

图2 震级平均超出量分布函数图
Fig.2 Diagram of mean excess function of seismic magnitude

图3 参数估计值随震级阈值选取的变化图
Fig.3 Variation diagram of parameter estimates according to the seismic magnitude threshold variation

图4a是概率图,横坐标表示实际数据的累积概率,纵坐标表示选用分布模型的累积概率; 图4b是分位数图,横坐标表示所选分布模型的分位数,纵坐标表示实际数据利用经验分布函数计算得到的分位数。从图4a,b可看出,大部分点均沿直线分布,因此不能否认所拟合的模型。由于ξ的估计值为负,因此最大震级分布存在有限上界,故重现水平图中的曲线为凸曲线,渐近地趋于有限值(图4c)。概率密度图中的密度曲线估计与直方图也显示一定的拟合(图4d)。因此4个诊断图均支持拟合的广义帕累托分布模型。

图4 广义帕累托分布拟合诊断图
Fig.4 Diagnostic plots of the Generalized Pareto Distribution fitted to seismic magnitude

利用最大似然法估计广义帕累托分布的参数,得到修饰的尺度参数和形状参数的估计值为(σ^{^/～},ξ^{^})=(1.31,-0.47),其中形状参数估计的标准差为0.18,其95%的置信区间为[-0.82,-0.12],形状参数的最大似然估计与置信区间都显示为负值,这表明龙门山地区在研究时间尺度范围内震级超出量分布具有上限值,这也与潜在震源区一定具有震级上限的认识相符。

将阈值、形状参数与修饰的尺度参数估计值代入式(7)得到震级上限为8.3。对震级上限的误差进行分析,首先需计算参数估计的协方差矩阵:

V={0.001 4 0 0

0 0.113 4 -0.057 3

0 -0.057 3 0.033 3}(10)

将式(10)协方差矩阵代入式(8),得到震级上限的标准差为0.50,由上限估计值和标准差可知震级上限置信度为95%的置信区间为[7.3,9.3]。