中国大陆地震空间分布模型检验

徐伟进, 李雪婧

中国地震局地球物理研究所， 北京 100081



中国大陆地震空间分布模型检验

徐伟进, 李雪婧

中国地震局地球物理研究所， 北京 100081

本文以中国大陆地震目录为基础资料，以泊松模型为零假设模型，并将Neyman-Scott空间丛集过程的各子模型设为检验模型，采用K-function 点过程分析法和最大似然估计法计算各模型参数，并以AIC准侧判定模型的拟合优度，来检验中国大陆地震空间分布模型.检验结果表明：泊松模型的拟合优度最差，说明地震在空间的分布不是完全随机的；广义Thomas模型的拟合优度最好，说明地震的空间分布是丛集的，可用由两个高斯核组成的广义Thomas模型较好地描述.研究结果还表明，同一研究区内，采用不同时段具有不同最小完整起始震级的地震目录计算得到的地震空间分布的丛集尺度几乎不变，这意味着地震空间丛集尺度不受小地震的控制，且可能与研究区的断层规模有关.

泊松模型； Neyman-Scott丛集模型； 地震空间分布； 中国大陆

1 引言

地震空间分布特征是地震学研究的一个重要方面.在区域地震预测和地震危险性分析中，建立经验地震空间活动性模型是一个主要的工作环节，这需要计算地震在空间格点的活动频度，采用的方法主要以核密度估计法为主(如Frankel，1995； Console and Murru, 2001; Zhuang et al.,2002)，而核密度估计中核函数的选择是尤为重要的，核函数是描述地震空间点分布的数学模型.Frankel(1995)在美国中东部地震区划中采用的地震空间模型为高斯模型，Zhuang 等(2002)基于ETAS模型的随机除丛方法中采用的也是高斯模型；Cao等(1996)在对美国加州的概率地震危险性分析中采用的是幂律模型；Woo(1996)则采用了地震分形分布模型.徐伟进和高孟潭(2012a)系统比较了上述三个模型，认为地震空间分布模型的选择对计算地震空间格点发生频度的影响是显著的.因此在计算地震空间活动性模型时地震空间分布理论模型的选择是非常重要的.

徐伟进(2012)的研究表明地震空间点分布格局是丛集的，即地震在空间上是一个丛集点过程.我们希望找到合适的数学模型去描述这一空间点过程.文中我们将采用目前应用最为广泛的Neyman-Scott(Archibald,1948)空间丛集过程及相似的点过程模型来描述地震空间分布，通过系统比较Neyman-Scott空间丛集过程的各子模型来选择合适的地震空间分布理论模型.

对于一组空间点，需要根据这组数据来估算各个模型的参数，本文将采用基于二阶统计量K-function(Ripley 1976，1977)的Maximum Palm Likelihood Estimation (MPLE)(Ogata and Katsura, 1991)方法估计模型参数，并根据AIC准则(Akaike, 1974)判断模型的优劣.事实上K-function(Ripley 1976，1977)本身也被广泛应用于点过程的丛集性分析与探测.Veen 和Schoenberg (2005)曾经采用K-function对南加州地区的地震空间分布做了研究.因此我们采用K-function和Maximum Palm Likelihood Estimation相结合的手段进行模型参数的估计、模型比较及模型的选择.

2 模型

2.1 泊松模型

空间点格局分析的零假设一般为随机分布模型，即泊松模型.泊松模型中地震在空间中的位置是完全随机的，并且每个地震事件的位置完全独立.地震在空间格点单位面积的发生频度是相同的，在一个统计单元内发生n个事件的概率密度函数可写为

(1)

其中a为统计单元面积，λ为事件单位面积的密度.λ=NR/AR,R为研究区域，NR为研究区内事件个数，AR为研究区面积.

2.2 Neyman-Scott 空间丛集模型

丛集分布是自然界中最为广泛的一种空间分布格局，Neyman-Scott模型是应用比较广泛的描述空间点分布格局的概率模型(Archibald,1948).根据Tanaka等(2008)的研究，Neyman-Scott空间丛集过程可以这样来描述：假设空间上有一个平稳泊松随机点过程，其强度为μ，我们把这些点称之为母点，每个母点随机产生M个独立、同分布的子点，M的分布由概率pm=P{M=m},m=0,1,…给出，其对应的均值记做ν.子点各向同性地分布在母点的周围，而距离是独立且同分布的，其概率密度函数为fτ(r)，其中τ代表概率密度函数中的一系列参数.Neyman-Scott点过程即是由这些子点组成的，其强度为λ=μν.

在下文研究中，参数θ=(μ,ν,τ)被看成是模型中的参数，将采用统计回归方法对各参数进行估计.

Neyman-Scott空间丛集过程中几个典型的模型如下：

Thomas process or Gaussian type model (Thomas,1949):

(2)

该模型最早由Thomas(1949)提出，Daley和Vere-Jones(1988)给出了距离分布的概率密度函数，即(2)式.这一模型中，子点以二维高斯分布分撒在母点周围.

广义Thomas model：

参考Stoyan 和Stoyan (1996)及Tanaka等(2008)的研究，还考虑了广义托马斯丛集模型，其距离分布概率密度函数为

(3)

该模型由两个不同丛集大小和距离尺度的二维高斯模型组合而成，其中σ1,σ2分别为两个高斯分布的参数，a称其为混合率.两个高斯模型母点的强度分别为μ1,μ2，母点周围子点的平均个数分别为ν1=ν2=ν，并且μ=μ1+μ2，μ1=aμ，μ2=(1-a)μ.

Inverse-power type model:

逆幂律模型的距离分布概率密度函数为

(4)

其中p>1,c>0，分别表示距离衰减的程度和尺度.

3 K-function空间点格局分析

3.1 经验K-function

Ripley′sKfunction(Ripley，1976)是最为熟知的用来分析空间点过程的工具.其广泛应用于各种空间点过程的模式分析.空间K-function的使用前提是点过程是一个平稳且各向同性的空间点过程，然而在实际应用中研究者们很难判断一个点过程是否是各向同性的.本文中我们将要估计的地震丛集分布尺度远小于研究区的分布尺度，因此可以将地震事件的空间分布近似看成是平稳且各向同性的.

本文通过非参数法求取了研究区的K-function，并与模型的理论K-function做了对比研究.K-function的估算公式为

(5)

其中，为dij第i个点与第j个点的距离；I为示性函数，dij小于t时，I等于1，否则为0；w为边界校正项，当以li为圆心，以t为半径的圆面完全处在统计区域内时，w为1，否则w为在统计区域内的面积与圆的面积的比值.

为了更加直观的表现出空间点在不同距离尺度上的分布特征，Besag(1977)根据K-function提出了L-function，其估计式为

(6)

3.2 理论K-fucntion

参照附录A的推导，可以导出Thomas model的理论K-function公式为

(7)

对于广义Thomas model的理论K-function公式为

(8)

同样对于(4)式逆幂律模型，我们无法代入参数直接求得其理论K-function，故采用数值计算方法求取.

特别地对于空间泊松点过程，其K-function公式可写为

(9)

4 模型参数估计及最优模型选择

4.1 模型参数估计

给定一个空间点过程，将采用统计学方法回归上述模型中的参数.在本研究中依然采用最大似然法求取模型中的参数.Ogata和Katsura (1991)导出了被称之为Palm likelihood function的公式：

(10)

其中，W为研究区域，往往把研究区域定为规则的矩形，N(W)为研究区域内点的个数，R是一个在计算Palm似然函数前已经确定的值，R的取值一般定为矩形短边的1/2.关于详细的参数估计及Palm likelihood function计算详见附录A.

对(4)式逆幂律模型，很难直接写出其似然函数式中积分的解析式，因此只能通过数值积分计算.

4.2 模型选择

基于上述介绍的方法，根据研究区中的地震目录，可以估计出Poisson模型、Thomas模型、广义Thomas模型及逆幂律模型中的各参数值，并采用Akaike Information Criterion(AIC)(Akaike, 1974)准则来判定模型的拟合优度:

AIC=2k-2ln(L)

其中，k为模型中参数的个数，L为最大似然函数值.一般认为AIC值越小，模型越合理.下文中AIC值均为各模型真实AIC值减去泊松模型AIC值后的结果，为的是与泊松模型进行比较，若AIC值小于0，则说明该模型优于泊松模型，AIC越小表示模型的拟合优度越好.

5 研究区及地震目录

本研究中选取了东北、华北、华南、川滇、青藏高原、新疆地震区和青藏高原地震区西端6个地区作为研究区域(见图1).选取6个研究区中1975年来M≥3.0地震和1970年来M≥4.0地震作为研究对象，其中东北地震区和华南地震区由于地震样本较少，只选取了1975年来M≥3.0地震.地震目录根据Mignan等(2013)的研究结果选取，地震目录中的余震采用Gardner和Knopoff(1974)的方法去除.

图1 研究区及地震空间分布Fig.1 Study regions and spatial distribution of earthquakes, black rectangles are study regions and black dots are earthquakes

6 结果分析

表1为根据华北地震区的地震目录采用上述参数估计方法得到的参数值及对数最大似然函数值和AIC值.每个模型的AIC值均减去Poisson模型的AIC值.通过观察表1中最后一列的AIC值发现Poisson模型的AIC值相对很大，这说明地震空间分布是丛集的.比较三个丛集模型，发现广义Thomas模型的AIC最小，这说明了本文验证的三个丛集模型中广义Thomas模型是最优的.此外，通过比较分别采用M≥3.0地震和M≥4.0地震计算得到的广义Thomas模型的σ1,σ2，发现采用两套地震目录计算得到的地震空丛集尺度是基本相同的，优势丛集尺度(若a>0.5则优势丛集尺度为σ1，否则为σ2)均为27 km左右，这在一定程度上暗示着地震空间丛集尺度是不受震级影响的(至少对于华北地区M≥3.0和M≥4.0地震是如此的).

图2为华北地震区两套地震目录计算得到的参数和非参数K-function的比较.图2(a，c)中红色实线是采用非参数法根据地震目录直接计算得到的K-function随距离变化的曲线，黑色实线为理论Poisson模型的K-function值，黑色虚线为Poisson模型的K-function值95%置信区间，可以发现R在200 km以内采用观测数据计算得到的K-function值均在黑色实线和虚线之上，这说明时间地震事件在空间的分布确实是呈现丛集的.比较三种丛集模型的曲线发现逆幂律丛集模型(图2(a，c)绿色曲线)与实际计算曲线吻合最差，这与表中三种丛集模型的AIC值中逆幂律模型AIC值最大是一致的.Thomas模型的曲线(图2(a，c)中紫色曲线)好于逆幂律模型.与实际数据吻合最好的是广义Thomas模型，从图中可以看出在200 km以内，广义Thomas模型(图中蓝色曲线)与实际数据吻合的非常好，有的部分几乎重合.图2(b，d)为实际观测数据测算的曲线及广义Thomas模型曲线及其95%置信区间曲线，我们发现在根据实际观测数据计算得到的K-function值完全落在广义Thomas模型95%的置信区间之内.95%置信区间是根据采用估计的参数模拟1000次理论丛集过程计算求得的.

此外，还采用上文介绍的方法估计了川滇、东北、华南、新疆地震区西端及青藏高原地震区等5个地震区的Poisson模型及各丛集模型的参数(表2).在这5个地震区只分析了Thomas模型和广义Thomas模型两个丛集模型.逆幂律模型由于其庞大的计算量并且从华北地震区的拟合情况看该模型

表1 根据华北地震区M≥3.0地震和M≥4.0地震计算得到的各模型参数及对数似然函数值和AIC值Table 1 MPLE parameters and AIC values of each model for M≥3.0 and M≥4.0 catalogs in Northern China

图2 华北地震区实测及各丛集模型的L-function曲线(a)、(c)中红色实线为根据实际观测数据计算所得，蓝色曲线为广义Thomas模型，紫色曲线为Thomas模型，绿色曲线为逆幂律模型，黑色实线及虚线分别为Poisson模型曲线及其95%置信区间； (b)、(d)中红色实线为根据实际观测数据计算所得，蓝色实线及虚线为广义Thomas模型曲线及其95%置信区间.Fig.2 The graphs of L(r)s across r for Northern ChinaIn graph (a) and (c), the red lines are estimated L(r) sacross r for M≥3.0 and M≥4.0 catalogs. The blue lines, purple lines and green lines are the best-fit lines between L(r)s and r of General Thomas model, Thomas model and Inverse-power model, respectively.The black solid lines are for Poisson model, and the black dash lines are 95% confidence limits. In graph (b) and (d), the red lines are estimated L(r)s across r for M≥3.0 and M≥4.0 catalogs. The blue solid lines are for General Thomas model, and the blue dash lines are 95% confidence limits.

表现最差，极有可能并不适合其他地震区，因此在该5个地震区未分析逆幂律模型.

通过比较每个地震区各模型的AIC值发现，在5个地震区中，Poisson模型的拟合优度是最差的，这说明了各地震区地震的空间分布都是丛集的，从而也可以认为我国大陆地区地震的空间分布是丛集的.

表2 川滇、东北、华南、新疆地震区西端及青藏高原地震区等5个地震区各模型参数及对数似然函数值和AIC值Table 2 MPLE parameters and AIC values of each model for M≥3.0 and M≥4.0 catalogs in Northeastern China, Southern China, Sichuan-Yunnan, Tibet and Western of Xinjiang seismic regions

图3 东北、华南、川滇、青藏高原地震区及新疆地震区西端等5个地震区实际数据测算及不同模型的L-function曲线其中红色实线是根据时间观测数据计算得到的L-function曲线，黑色实线为Poisson模型，绿色实线为Thomas模型，蓝色实线为广义Thomas模型.Fig.3 The graphs of L(r)s across r for Northeastern China, Southern China, Sichuan-Yunnan, Tibet and Western of Xinjiang seismic regionsThe red lines are estimated L(r)s for M≥3.0 and M≥4.0 catalogs, the blue lines, green lines and black lines are L(r)s of General Thomas model, Thomas model and Poisson model, respectively.

同样，在这5个地震区中广义Thomas模型的AIC值均为最小，这说明该模型的拟合优度最好.通过比较M≥3.0和M≥4.0两套目录计算得到的广义Thomas模型的σ1,σ2，发现在每个地震区内采用这两套地震目录计算得到的优势地震丛集尺度差别不显著，这也再次印证了在同一研究区内地震空间丛集的丛集尺度可能不受震级影响.

图3为川滇、东北、华南、新疆地震区西端及青藏高原地震区等5个地震区的K-function曲线，其中红色实线是根据时间观测数据计算得到的K-function曲线，紫色实线为Poisson模型，蓝色虚线为Thomas模型，蓝色实线为广义Thomas模型.通过观察发现，每个地震区现采用M≥3.0和M≥4.0两套目录计算得到的K-function均与理论广义Thomas模型的K-function吻合最好，这与表2中广义Thomas模型AIC值最小是一致的.

7 结论与讨论

文中根据中国大陆的地震目录资料，采用K-function点过程分析法和最大似然估计法，分析了中国大陆6个地区的地震空间分布特征，得出如下有益结论：

在我国大陆地区地震空间分布是丛集的.并且丛集模式可以用由两个高斯核组成的广义Thomas模型进行描述.这一研究结果为地震预测、地震危险性分析及其他地震学研究中经验空间地震活动性模型的计算提供了一定的理论依据.

同一研究区内，采用不同时段具有不同最小完整起始震级的地震目录进行计算得到的地震丛集尺度是大致相同的，这说明了中国大陆地震空间丛集分布的尺度可能是不受小地震影响的，该尺度可能与研究区内断层规模有关.从物理上讲一个丛集区就是一个断层带或者应力集中区，大量地震的震中分布在断层面上，因此只要地震事件足够多，无论采用多大震级的地震进行统计得出只是同一处断层带的尺度，故计算结果是相当的.在不同研究区内地震空间丛集尺度不同.在西部地区，计算得到的地震空间丛集尺度相对较大，这可能是与西部地区地质环境复杂，地震发生强度大频度高有关.在东部地区由于地震强度相对较小，频度也相对低的多，因此计算的地震空间丛集尺度就相对较小.总之，各研究区地震空间丛集尺度的不同反应了地震地质环境的不同，其更深层次的地震地质学意义需要更进一步的研究.

从表1和表2中各地区丛集尺度σ值的大小可以看出该值并不是研究区内断层的实际规模尺度.地震的空间丛集尺度显示的是断层带上破裂相对集中地区的大小，即地震事件初始破裂点的聚集地区.因此地震空间丛集尺度为概率地震危险性分析中潜在震源区几何尺度的确定以及地震预测中地震危险区的划定提供了一定的定量依据.

从统计学上讲，本文的结论是预料之中的，根据中心极限定理，对于任何一个空间点过程，即使不符合高斯过程，都可以构造一个由多个高斯核组成的模型去逼近这个点过程，当然，模型会显得很复杂，不利于研究结果的应用.在科学研究中我们期望用简单的数学模型去描述复杂的物理过程.本文中采用含有两个高斯核的混合模型就能很好地去拟合实际的地震空间过程，这个结果是比较理想的.

附录A：Palm likelihood function

给定一个空间点过程，将采用统计学方法回归上述模型中的参数.在本研究中依然采用最大似然法求取模型中的参数.Ogata和Katsura (1991)导出了被称之为Palm likelihood function的公式：

(A1)

其中，W为研究区域，往往把研究区域定为规则的矩形，N(W)为研究区域内点的个数，R是一个在计算Palm似然函数前已经确定的值，R的取值一般定为矩形短边的1/2.

在计算Palm似然函数时用到了Palm intensitiesλo.据Daley和Vere-Jones (1988)研究：

(A2)

其中，λ即为点过程的强度，K(r)为K-function.对于Neyman-Scott过程，其二阶统计量K-function可用下式求取(Stoyan and Stoyan 1994)：

(A3)

其中：

(A4)

根据(A2)(A3)(A4)式，(A1)可重写为

(A5)

对于Thomas model和广义Thomas model，(A4)式的积分我们可以写出其解析式.对于Thomas model：

(A6)

对于广义Thomas model：

(A7)

将(A6)式(A7)式带入(A5)式，则Thomas model和广义Thomas model的Palm似然函数可写为：

Thomas model：

(A8)

广义Thomas model：

(A9)

变化各个参数值，使似然函数取最大值时，即求得参数的最优估计值.

Akaike H. 1974. A new look at the statistical model identification.IEEETransactionsonAutomaticControl, 19(6): 716-723.

Archibald E E A. 1948. Plant populations: I. A new application of Neyman′s contagious distribution.AnnalsofBotany, 12(47): 221-235.

Besag J E. 1977. Comment on “modelling spatial patterns” by B.D. Ripley.JournaloftheRoyalStatisticalSociety, Series B, 39: 193-195.

Cao T Q, Petersen M D, Reichle M S. 1996. Seismic hazard estimate from background seismicity in Southern California.Bull.Seismol.Soc.Am., 86(5): 1372-1381.

Console R, Murru M. 2001. A simple and testable model for earthquake clustering.JournalofGeophysicalResearch, 106(B5): 8699-8711, doi: 10.1029/ 2000JB900269.

Daley D J, Vere-Jones D. 1988. An Introduction to the Theory of Point Processes. Berlin: Springer-Verlag.

Frankel A. 1995. Mapping seismic hazard in the central and eastern United States.SeismologicalResearchLetters, 66(4): 8-21.

Gardner J K, Knopoff L. 1974. Is the sequence of earthquakes in Southern California, with aftershocks removed, Poissonian?.Bull.Seismol.Soc.Am., 64(5): 1363-1367.

Mignan A, Jiang C, Zechar J D, et al. 2013. Completeness of the mainland China earthquake catalog and implications for the setup of the china earthquake forecast testing center.Bull.Seismol.Soc.Am., 103(2A): 845-859.

Ogata Y, Katsura K. 1991. Maximum likelihood estimates of the fractal dimension for random spatial patterns.Biometrika, 78(3): 463-474.

Ripley B D. 1976. The second-order analysis of stationary point processes.JournalofAppliedProbability, 13(2): 255-266.

Ripley B D. 1977. Modelling spatial patterns (with discussion).JournaloftheRoyalStatisticalSociety,SeriesB, 39: 172-212.

Stoyan D, Stoyan H. 1996. Estimating pair correlation functions of planar cluster processes.BiometricalJournal, 38(3): 259-271.

Tanaka U, Ogata Y, Stoyan D. 2008. Parameter estimation and model selection for Neyman-Scott point processes.BiometricalJournal, 50(1): 43-57.

Thomas M. 1949. A generalization of Poisson′s binomial limit for use in ecology.Biometrika, 36(1-2): 18-25.

Veen A, Schoenberg F P. 2006. Assessing spatial point process models using weightedK-functions: analysis of California earthquakes.∥ Baddeley A, Gregori P, Mateu J, et al., eds. Case Studies in Spatial Point Process Modeling. NY: Springer, 293-306.

Woo G. 1996. Kernel estimation methods for seismic hazard area source modeling.Bull.Seismol.Soc.Am., 86(2): 353-362.

Xu W J. 2012. Study on Space-time Statistical Distribution Models of Earthquakes in Seismic Hazard Analysis [Ph. D. thesis] (in Chinese). Beijing: Institute of Geophysics, China Earthquake Administration.

Xu W J, Gao M T. 2012. Comparison of the smoothing functions in the analysis of spatially smoothed seismicity.ActaSeismologicaSinica(in Chinese), 34(2): 244-256.

Zhuang J C, Ogata Y, Vere-Jones D. 2002. Stochastic declustering of space-time earthquake occurrences.J.Am.Stat.Assoc., 97(458): 369-380.

附中文参考文献

徐伟进. 2012. 地震危险性分析中地震时空统计分布模型研究 [博士学位论文]. 北京: 中国地震局地球物理研究所.

徐伟进, 高孟潭. 2012. 空间光滑地震活动性模型中光滑函数的比较研究. 地震学报, 34(2): 244-256.

(本文编辑 汪海英)

Statistical test of spatial distribution model of the Earthquakes in Chinese Mainland

XU Wei-Jin, LI Xue-Jing

InstituteofGeophysics，ChinaEarthquakeAdministrator，Beijing100081，China

In this paper, we test the spatial distribution model of the earthquakes in Chinese Mainland. We take Poisson model as null hypothesis and test Neyman-Scott and similar point processes. TheK-function point processes analysis tool and maximum likelihood estimation method are used to calculate parameters of each model. We use the Akaike Information Criterion (AIC) to compare the goodness-of-fit of these models. The results show that the goodness-of fit of Poisson model is the worst compared with that of any other model, while the goodness-of-fit of General Thomas model is the best. It can be confirmed that the spatial distribution of earthquakes is indeed significantly clustered (General Thomas model) rather than completely random (Poisson model). The best fit is got by the General Thomas model which contains two Gaussian kernels. The scale of cluster of spatial distribution of earthquakes does not change for each seismic sub-catalog in the same research region. This means that the scale of cluster is not controlled by small earthquakes, and that may be related with the scale of faults in research regions.

Poisson model; Neyman-Scott cluster model; Seismic spatial distribution; Chinese Mainland

10.6038/cjg20160911.

中国地震局地震行业专项(201408014)和中国地震局地球物理研究所基本业务专项(DQJB12B26)联合资助.

徐伟进，男，1982年生，副研究员，主要从事地震活动性及地震危险性分析方法研究.E-mail:wjxuwin@163.com

10.6038/cjg20160911

P315

2015-11-05，2016-07-26收修定稿

徐伟进, 李雪婧. 2016. 中国大陆地震空间分布模型检验. 地球物理学报，59(9):3260-3268,

Xu W J, Li X J. 2016. Statistical test of spatial distribution model of the Earthquakes in Chinese Mainland.ChineseJ.Geophys. (in Chinese),59(9):3260-3268,doi:10.6038/cjg20160911.