2016, Vol. 59
地震空间分布特征是地震学研究的一个重要方面.在区域地震预测和地震危险性分析中，建立经验地震空间活动性模型是一个主要的工作环节，这需要计算地震在空间格点的活动频度，采用的方法主要以核密度估计法为主(如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-functionRipley′s K function(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) |
参照附录A的推导，可以导出Thomas model的理论K-function公式为
|
(7) |
对于广义Thomas model的理论K-function公式为
|
(8) |
同样对于(4)式逆幂律模型，我们无法代入参数直接求得其理论K-function，故采用数值计算方法求取.
特别地对于空间泊松点过程，其K-function公式可写为
|
(9) |
给定一个空间点过程，将采用统计学方法回归上述模型中的参数.在本研究中依然采用最大似然法求取模型中的参数.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)准则来判定模型的拟合优度:
|
其中，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 |
表 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地震是如此的).
|
表 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为华北地震区两套地震目录计算得到的参数和非参数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次理论丛集过程计算求得的.
|
图 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 China In 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个地震区的Poisson模型及各丛集模型的参数(表 2).在这5个地震区只分析了Thomas模型和广义Thomas模型两个丛集模型.逆幂律模型由于其庞大的计算量并且从华北地震区的拟合情况看该模型表现最差，极有可能并不适合其他地震区，因此在该5个地震区未分析逆幂律模型.
|
|
表 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 |
通过比较每个地震区各模型的AIC值发现，在5个地震区中，Poisson模型的拟合优度是最差的，这说明了各地震区地震的空间分布都是丛集的，从而也可以认为我国大陆地区地震的空间分布是丛集的.
同样，在这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值最小是一致的.
|
图 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 regions The 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. |
文中根据中国大陆的地震目录资料，采用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. IEEE Transactions on Automatic Control , 19 (6) : 716-723. DOI:10.1109/TAC.1974.1100705 | |
| Archibald E E A. 1948. Plant populations:I. A new application of Neyman's contagious distribution. Annals of Botany , 12 (47) : 221-235. | |
| Besag J E. 1977. Comment on "modelling spatial patterns" by B. D. Ripley. Journal of the Royal Statistical Society, 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. Journal of Geophysical Research , 106 (B5) : 8699-8711. DOI:10.1029/2000JB900269 | |
| Daley D J, Vere-Jones D. An Introduction to the Theory of Point Processes. Berlin: Springer-Verlag, 1988 . | |
| Frankel A. 1995. Mapping seismic hazard in the central and eastern United States. Seismological Research Letters , 66 (4) : 8-21. DOI:10.1785/gssrl.66.4.8 | |
| 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. DOI:10.1785/0120120052 | |
| Ogata Y, Katsura K. 1991. Maximum likelihood estimates of the fractal dimension for random spatial patterns. Biometrika , 78 (3) : 463-474. DOI:10.1093/biomet/78.3.463 | |
| Ripley B D. 1976. The second-order analysis of stationary point processes. Journal of Applied Probability , 13 (2) : 255-266. DOI:10.1017/S0021900200094328 | |
| Ripley B D. 1977. Modelling spatial patterns (with discussion). Journal of the Royal Statistical Society, Series B , 39 : 172-212. | |
| Stoyan D, Stoyan H. 1996. Estimating pair correlation functions of planar cluster processes. Biometrical Journal , 38 (3) : 259-271. DOI:10.1002/(ISSN)1521-4036 | |
| Tanaka U, Ogata Y, Stoyan D. 2008. Parameter estimation and model selection for Neyman-Scott point processes. Biometrical Journal , 50 (1) : 43-57. DOI:10.1002/(ISSN)1521-4036 | |
| Thomas M. 1949. A generalization of Poisson's binomial limit for use in ecology. Biometrika , 36 (1-2) : 18-25. DOI:10.1093/biomet/36.1-2.18 | |
| Veen A, Schoenberg F P. 2006. Assessing spatial point process models using weighted K-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. Acta Seismologica Sinica (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. DOI:10.1198/016214502760046925 | |
| 徐伟进. 2012. 地震危险性分析中地震时空统计分布模型研究[博士学位论文]. 北京:中国地震局地球物理研究所. | |
| 徐伟进, 高孟潭. 2012. 空间光滑地震活动性模型中光滑函数的比较研究. 地震学报 , 34 (2) : 244–256. | |