安宁河–则木河–小江断裂带区域构造地震发生的时空概率分布估计-凯发国际一触即发

1. 引言

在川滇地区，存在着多个重要的活动断裂带，如鲜水河、龙门山、安宁河–则木河–小江、红河以及丽江–小金河断裂带。其中，安宁河–则木河–小江断裂带位于“川滇菱形地块”的东侧边界区域。安宁河断裂带是一条东西走向的逆冲断裂带，主要由北侧断块和南侧断块组成，两侧断块相对位移显著，主要沿着安宁河流域东西走向延伸，贯穿了云南省安宁市和昆明市一带，全长约160公里，总体走向为北南向。而则木河断裂带则连接安宁河断裂带与小江断裂带，全长约150公里，主要倾向为东北向，局部为西南向。小江断裂带位于四川省西南部和云南省北部，主要沿着大凉山地区东西走向延伸，东起四川省泸定县，西至云南省盐津县，全长约280公里，该断裂带主要以左旋走滑为主。

研究安宁河–则木河–小江断裂带区域的滑动速率、断层深度以及潜在强震风险，对于深入了解地震发生规律，提高地震安全性分析的科学水平至关重要[1]。

目前，地震研究和预测已成为一个多学科的综合性问题。许多学者结合统计学、地质学的方法和模型，对地震发生的时间、地点和危险性进行了预测，并提供了其他重要资料。这种综合性方法的运用，为我们更好地理解地震活动提供了有力的工具，并且为地震灾害防范和减轻工作提供了重要的科学依据。在强震危险性预测中，近年来许多学者(wensnosky, s. t. 1986; nishenko, s. p. et al., 1987；闻学泽，1990)利用地震原地复发的条件概率模型，在概率意义下预测不同时间范围内的强震风险，以确定未来不同时间范围内强震风险区段[2]。程万正根据断裂带出现的地震活动性、跨断层地壳形变和微观异常现象分析了强震预测及预测应用问题，对经验性预测进行了辩证的思考[3]。唐昌荣、黄祖智等人利用概率复发模型预测未来30年各活动段强震复发可能性，并分析了鲜水河断裂带恰叫段，乾宁段和色拉哈段，及安宁河断裂带西昌–冕宁段和野鸡洞段更易发生7级地震[4]。闻学泽通过小江带地震潜势概率预估研究，初步建立了该断裂带上地震发生时间区间统计模型，并介绍了以各剖面上1次地震规模作为背景来预估平均复发区间数学方法，基于此，进一步提出了一种适用于准时间可预测复发的断裂带上进行地震潜势实时概率分析的新方法，并利用该方法估算了小江断裂带在1991年至2005年间各段地震复发的条件概率和概率增益，得出该断裂带东川南部到嵩明段和澄江到华宁段将来复发强，大震概率较大[5]。

在地质统计学的应用中，指示克里格被广泛采纳为一种方法，但在观察数据中常常出现某些特别的值，而这些数据并不遵循对数正态的分布规律，这直接妨碍了变异函数的稳定性。倘若采用参数地质统计的技术手段，就必须排除这些独特的数据值并对数据执行非线性转化，目的是让数据的概率分布遵循正态关系，然而，这可能会对数据变量在空间内的变化产生不良的真实影响。在非参数地质统计学领域，指示克里格法被认为是一种处理偏数据的有力策略，能够在不需要移除实际存在的重要特异值的前提下处理各种现象，并成功地降低特异值对变异函数稳健型的效果。指示克立格方法是一个被广泛采用的非参数地质统计技术，它的名称源于它将区域变量的研究转化为对相关指令函数的探讨。关于指示不克里格方法的学术研究和实践应用，国内和国际学术人士已经进行了大量的研究和尝试，但大多数的研究焦点集中在单元指示克里格在某个尺度上的具体应用上。

本文选取中国安宁河–则木河–小江断裂带及其周边区域1983年至今的地震事件作为研究对象，从发震时刻间隔的角度，利用指示克里格方法对安宁河–则木河–小江断裂带区域构造地震发生的时空概率分布进行估计，对分析所得的时空概率分布估计进行辨证的思考。

2. 指示克里格基本原理

在地质属性空间分布预测方面，对连续型数据进行编码是一种最早应用的方法。二进制转换在克里格插值中起着关键作用，通过合理的二进制转换，能够根据阈值z k 原始数据可以转换为适合插值分析的格式，如下式：

i( x; z k )={ 1, ifz( x )≥ z k 0, ifz( x )≤ z k (1)

由(1)式可知，指示变换是一种对原始数据进行的非线性变换，这一转换可以改善数据的适用性，包括消除偏斜性、降低方差、处理异常值、增强空间相关性，以提高克里格插值结果的准确性和可靠性。在给定z的条件下，i( x h, z k ) 和i( x; z k ) 分别代表向量h分隔的两个区域化变量，于是，可以将指示变差函数及实验指示变差函数定义如下：

r( h; z k )= e{ [ i( x h; z k )−i( x; z k ) ] } 2 (2)

r( h; z k )= 1 2 n ( h ) [ i( x h; z k )−i( x; z k ) ] 2 (3)

在式中，n ( h ) 表示在距离h处的数据对数量。式(3)表明，可利用指示函数i( x α : z k ) 计算距离为h的数据点对之间的实验指示变差，而不是用样本值x α 本身。因此，指示变差函数r( h; z k ) 对z ( x ) 有极强的稳健性，意味着即使数据中存在极高或极低的值，也不需要常规的删除处理，这样就极大化的保留了数据信息的完整性。

首先收集了研究区域s内的n个区域化变量的观测值，{ z ( x α ) , x α ∈s,α=1,2,⋯,n } ，接着，根据给定的l个阈值{ z k ,k=1,2,⋯,l } 后，将这些观测值转换为指示函数域，表示为{ i( x α ; z k ),α=1,2,⋯,n } 。利用这些指示值，可以解出指示克里金方程组，从而获得模拟网格中各节点位置处的条件累积分布函数，该函数描述了每个位置处的不确定性概率模型[6] [7]，如下式：

f x ( z k |n )=prob{ z ( x ) ≤ z a | a ( x α ) ,α=1,2,⋯,n } (4)

e-型估值通过对条件累积分布函数进行积分来实现，具体公式如下：

∫ −∞ ∞ zd f x * ( z|n ) ≈ ∑ i=1 k z i [ f x * ( z|n ) ] − f x * ( z i−1 |n ) (5)

其中z min < z 0 < z 1 <⋯< z k−1 < z k < z max 表示对区间[ z min , z max ] 的阈值划分z i ∈[ z i−1 , z i ] ，i=1,2,⋯,k 。e-型估值之所以被采用，是因为它不仅满足最小二乘标准，而且能够计算衡量条件不确定性的条件方差σ 2 ( μ ) ，从而评估未知点的不确定性。此外，e型估值的误差平方的条件期望值达到了最小。具体地说，可以利用公式(4)建立的条件累积分布直方图来计算区域化变量在指定块段或点内大于或小于任意指定阈值的概率，而公式(5)则可以得到区域化变量在指定块段或点的平均值估计[8]。这些结果对于地质学或地理信息系统等领域的数据分析和预测非常有用，因为它们可以帮助我们更好地理解和评估未知点的值以及与之相关的不确定性。

指示克立格法以指定阈值对数据进行二元指示变换，有助于减轻偏态分布和特异值对变异函数稳健性的影响。借助指示变异函数，能够获取时间和空间上不同阈值下震级的指示变异函数模型及相关参数。指示半变异函数可遵循不同类型的模型，如球状模型、指数模型以及高斯模型。块金值co是由测量误差与空间变异性综合造成的，它体现了区域化变量的随机性与异质性；基台值 sill代表了由结构性和随机性变化组成的总变化，它代表了变量在空间上的异质性。co/sill值体现出空间自相关性的强弱，一般将其划分为小于25%、25%~75%和大于75%三个等级，分别与较强、中等和很弱的空间自相关性相对应。当在整个尺度上具有恒定变异时，它的数值会接近100%。变程是指各变量在空间上的异质性的尺度函数，它是度量各变量在变程范围内变化的最大值，并在变程范围内表现出空间上的自相关。区域化变量的空间变异性是由结构和随机性两个方面综合作用的结果，其中，非人为结构性是导致区域化变量空间连续的主要原因，而人为随机性是导致区域变量空间自相关性减弱的主要原因[9]。

3. 指示变异函数模型构建

3.1. 数据来源及预处理

本文采用1983年7月15日~2023年12月31日中国地震台网()地震目录数据。地震事件序列具有波动性和非线性特点，包括发震时间、经纬度和震级等因素。因此，通过数据预处理，可以将地震事件序列确定为特定空间和一定震级下的数据，即安宁河–则木河–小江断裂带及其周边区域的经纬度范围内(24.0˚n~30.0˚n，101.0˚e~104.0˚e) m_s 3.0及以上的时间序列。

如表1所示，1983~2023年间安宁河断裂带附近发生的震级在3级及以上的地震次数为310次，其中3~4.5级地震发生的概率最大，占总数的91.29%；4.5~6级地震次之，占总数的8.06%；6~7级地震发生的概率最小，只占总数的0.65%。

table 1. distribution characteristics by magnitude

表1. 按震级分布特征

根据地震发生的时间间隔与发震顺序画出时间间隔关系曲线，如图1，可以看出发震间隔曲线图中有两段时间间隔很小的近乎直线段，分别在发震的第11,342天和第14,297天左右，对应时间段约为2014年8月13日左右和2022年9月15日左右，也就是云南鲁甸6.8级地震和四川甘孜州泸定县6.8级地震的时间段，说明大多数弱震是伴随着一次强震发生而发生的。

figure 1. time series of earthquakes of magnitude 3.0 and above in the aninghe-zemuhe-xiaojiang fault zone

图1. 安宁河–则木河–小江断裂带区域3.0级及以上地震时序图

3.2. 地震活动性研究

为了使地震带数据分析结果精确度更高，本文利用古登堡李希特的地震重复律关系(g-r关系)来确定出了各地震带发震数据的最小完备震级[8]。g-r关系式如下：

lgn=a−bm

式中m是震级，n是震级 > m的地震次数，a、b为常数。针对地震带的1983年7月到2023年10月的发震数据，通过最小二乘法拟合出g-r关系式，得到图2所示的结果。

基于拟合结果示意图图5，可以确定安宁河–则木河–小江断裂带发震数据的最小完备震级是ms = 3.0。为了分析结果的一致性和可比性，断裂带去除小于震级m = 3.0的发震数据。最终整理得到安宁河–则木河–小江断裂带发震数据共310条。进一步处理后再得到地震带的发震强度(震级)、发震时刻(间隔)以及发震频率(次数)序列的数据，并按照发震震级设置值(m3.0, m (4.5, 6.0), m6.0)，最后各得到三段数据。

figure 2. the fitting result of g-r relationship of anninghe-zemuhe-xiaojiang fault zone

图2. 安宁河–则木河–小江断裂带g-r关系式拟合结果图

4. 时空概率分布指示克里格计算

4.1. 发震时间概率分布估计

时间维度上，按分组计算半变异函数后分别对其进行模型拟合，模型结果见图3。在时间维度中，各阈值变异函数模型的参数值见表2，从中可以看出各模型块金值co都较小，说明本研究尺度上由于采样误差、短距离的变异、随机和固有变异引起时间上的变异不大。

(a) 3.0级及以上拟合模型 (b) 4.5级及以上拟合模型

(c) 6.0级及以上拟合模型

figure 3. semi-variance function fitting model under different thresholds in time dimension

图3. 时间维度中不同阈值下的半变异函数拟合模型

table 2. theoretical model of variation function for each threshold group in time dimension

表2. 时间维度中各个阈值组的变异函数理论模型

由表2可知，在时间维度中，阈值为3.0级的分组符合高斯模型，而4.5和6.0级的分组皆符合稳定模型，且三个模型的co/sill值分别为1.45%、0%和0%均小于25%，表现出较强的空间自相关性，即空间变异主要由结构型因素作用。在变程上，阈值为3.0级、4.5级和6.0级分组分别为3923天、801天和582天，说明在阈值为3.0级的分组中，地震发生的3923天里与其他3.0级及以上地震的发生存在一定的相关性，一次地震可能是未来十一年内地震发生的因素之一；在阈值为4.5级的分组中，4.5级及以上的地震发生的801天范围内，与其他4.5级及以上地震的发生存在一定的相关性；在阈值为6.0级的分组中，6.0级及以上的地震发生的522天范围内，与发其他6.0级及以上地震发生存在一定的相关性。

将变异函数模型的参数导入到arcgis中，带入10%测试数据进行指示克里金插值，绘制相应阈值的概率曲线图，同时画出真实数据对比情况，红色虚线为预测值分界线，红色虚线前为模型拟合情况，虚线后为预测值，预测结果见图4。

从图4中的概率曲线中可以清楚看到，根据图中预测值分界线可知，阈值为3.0级和4.5级分组以2023年5月12日为界限2023年10月12日以后为预测区域；阈值为6.0分组以2022年9月5日天为界限，2022年9月5日以后为预测区域。从图中可以看出数据预测存在一个明显的预测概率曲线变化区间范围，且概率曲线在达到峰值后逐渐下降，后趋于稳定值本文将根据以上的模型拟合及预测结果，与真实数据对比，讨论预测区间内发生地震的情况。

图4(a)中，3.0级阈值组预测概率曲线变化区间范围为第2023年10月12日至2025年2月22 (共计652天)，预测概率大于0.8的范围为第2024年6月15日至2025年2月28日(共计259天)，预测区间占变化区间的39.7%。说明未来几年内安宁河–则木河断裂带发生3.0级以上地震的概率较大。此外，可以看出，阈值为3.0级模型对于未来发震概率估计结果为：2025年1月30日发震概率达到最大值0.92，其前后141天的概率皆超过0.85。

(a) 3.0级及以上

(b) 4.5级及以上

(c) 6.0级及以上

figure 4. probability estimation of time dimensions under different thresholds

图4. 不同阈值下时间维度的概率估计

图4(b)中，4.5级阈值组预测概率曲线变化区间范围为第2023年10月12日至2025年2月22 (共计652天)，预测概率大于0.7的范围为第2023年10月15至2024年5月11天(共计209天)，预测区间占变化区间的32.1%，说明未来几年内安宁河–则木河断裂带有一定概率发生4.5级以上地震，但概率偏低。此外可以看出，阈值为4.5级模型对于未来发震概率估计结果为：2024年1月15日发震概率最大值0.78，其24天的概率皆超过0.7。

图4(c)中，6.0级阈值组预测概率曲线未能给出有效预测区间，究其原因可能是由于数据量过少而导致模型信息量提取不够。

总的来说，在阈值为3.0级与4.5级分组的模型预测上，其准确度较高，说明指示克里金在对未来发震时间的区间预测上有着一定的效果。

4.2. 空间维度的概率分布估计

而对于空间维度，模型的拟合需根据不同方向的半变异函数综合考虑，本次研究在综合空间维度12个方向上的各向异性后得到了最终的指示变差函数，见图5。在空间维度中，各阈值变异函数模型的参数值见和表3。从中可以看出各模型块金值co都较小，说明本研究尺度上由于采样误差、短距离的变异、随机和固有变异引起空间上的变异不大。

图5给出了各阈值组中各向异模型与综合半变异函数模型。由表3可知，在空间维度中，阈值为3.0分组符合稳定模型，4.5级和6.0级的分组符合高斯和稳定模型。其中，阈值为3.0级、4.5级和6.0级分组的模型的co/sill值分别为0.00%、20.65%和57.35%，3.0级、4.5级的co/sill值小于25%，表现出较强的空间自相关性，而6.0级分组的模型的co/sill值介于25%~75%，表现出中等的空间自相关性。而在变程上，阈值为3.0级、4.5级和6.0级分组分别为108.78 km、131.43 km、212.23 km，说明在阈值为3.0级的分组中，地震发生的108.78 km范围内，与其他3.0级以上地震的发生存在一定相关性；在阈值为4.5级的分组中，4.5级及以上的地震发生的183.22 km范围内，与其他4.5级及以上地震的发生存在一定相关性；在阈值为6.0级的分组中，6.0级及以上的地震发生的212.23 km范围内，与其他6.0级及以上地震的发生存在一定相关性。

(a) 3.0级及以上各向异性(左)与综合模型(右)

(b) 4.5级及以上各向异性(左)与综合模型(右)

(c) 6.0级及以上各向异性(左)与综合模型(右)

figure 5. anisotropic models and synthetic semi-variogram models under different thresholds in spatial

图5. 空间维度中不同阈值下的各向异模型与综合的半变异函数模型

table 3. theoretical model of variation function for each threshold group in spatial dimension

表3. 空间维度中各个阈值组的变异函数理论模型

将空间维度的变异函数模型的参数导入到arcgis中，带入10%测试数据进行指示克里金插值，进行指示克里金插值，得到各阈值组的发震概率分布图，同时画出真实数据对比图，见图6~8。

figure 6. comparison of probability spatial distribution and data of grade 3.0 and above

图6. 3.0级及以上的概率空间分布与数据对比

figure 7. comparison of probability spatial distribution and data of level 4.5 and above

图7. 4.5级及以上的概率空间分布与数据对比

figure 8. comparison of probability spatial distribution and data of level 6.0 and above

图8. 6.0级及以上的概率空间分布与数据对比

3.0级及以上的概率空间分布与位置对比图见图6。从图中可以看出，在3.0级及以上的分组中，指示克里金很好地用发震概率的方式估计了发震位置，测试点基本出现在高概率区域。

4.5级及以上的概率空间分布与位置对比图见图7(b)。对比图6~8，可以看出阈值为4.5的空间概率图的高概率位置更加具体，从图中可以看出，在4.5级及以上的分组中，高概率区域主要在曲靖–昆明、华平–攀枝花和石棉–冕宁一带。

6.0级及以上的概率空间分布与位置对比图见图8从图中可以看出，高概率区域主要集中在曲靖、攀枝花、九龙一带。在6.0级及以上的分组中，由于数据量减少到13条，指示克里金法给出的高概率区域较4.5级分组更大，然而这种现象是样本数据量过小导致的，而不是真实的发生概率区域大。其高概率中心仍在4.5级分组区域中。在准确度上将测试数据基本出现在高概率区域。

以真实数据为基准，对阈值为3.0、4.5和6.0分别定义预测概率大于0.92、0.89和0.94为高概率区、反之为低概率区。表4是空间维度不同阈值条件下的高低概率区面积统计情况。由表可以看出，在空间维度上，阈值为3.0、4.5、6.0级的高概率区域栅格数为702、609和105，占总区域的比例为30.64%、26.76%、和7.8%。即随着阈值增大，各维度高概率区明显减少，这与断裂带地震发生的聚集型特点相符。

table 4. the proportion of probability interval of different threshold conditions in space dimension

表4. 空间维度下不同阈值条件的概率区间比例构成

对比可知，在空间维度上，发震高概率区则主要集中在(101.5˚e, 26˚n)至(103.5˚e, 30˚n)直线型断裂带上，当阈值为3.0级、4.5级和6.0级时，其区域范围内的发震数据为总发震数据的85.09%、76.07%、95.11%，主要分布在石棉县、巧家县、曲靖市、攀枝花市、九龙县，其余均演化为低概率区。

总的来看，不同维度上，大阈值的高概率区域均被包含在上一级小阈值的高概率区域，并且在空间分布上具有明显的相似性，空间维度上尤其是高震级发生的概率区主要集中在曲靖、石棉、攀枝花一带。对比分析不同阈值发震概率图，可以发现高阈值的高概率区包含在小阈值的高概率区，其分布区域随阈值的减小而扩散。

在基于指示克里金的时空概率估计上，3.0级阈值组模型对于未来发震的估计结果为：2025年1月30日前后50天里，在石棉县–冕宁县–宁南县–巧家县一带发生大于3.0级地震的概率为0.822；4.5级阈值组模型对于未来发震的估计结果为2024年1月15日前后17天里，在陆良县–曲靖市、九龙县–石棉县、华坪县–宁蒗彝族自治县一带发生大于4.5级地震的概率为0.70；阈值为6.0组模型对于未来发震未能做出有效估计。但给出了6.0级及以上地震高概率发震区域为九龙县、曲靖市、寻甸回族彝族自治县一带。

5. 结论

指示克里金法可以分析地震目录数据的震级变化，地震系列的震级可以看作是一个区域化的变量，有随机成分和结构成分，用变异图函数表征，并能进行推断和建模。从时间维度上看，虽然从技术上讲，根据拟合的变异图，可以对已知的最后一次地震的后续地震进行模拟，但外推的有效性会随着外推距离的增大而减小。而在空间维度，指示克里金的内插结果是比较准确的。最后，本章利用指示克里金法对未来地震的发生做出了估计，将指示克里金时间与空间分布估计结果做乘积，得到时空概率估计。根据结果，有以下几点结论：1) 3.0级阈值组模型对于未来发震的时空概率估计结果为：2025年1月30日，在石棉县–冕宁县–宁南县–巧家县一带发生大于3.0级地震的概率最大，为0.822，其前后50天概率大于0.804；2) 4.5级阈值组模型对于未来发震的估计结果为：2024年1月15日，在陆良县–曲靖市、九龙县–石棉县、华坪县–宁蒗彝族自治县一带发生大于4.5级地震的概率最大，为0.63，其前后17天概率大于0.616；3) 6.0级阈值组模型对于未来发震时间未能做出有效估计，但给出了6.0级及以上地震高概率发震区域为九龙县、曲靖市、寻甸回族彝族自治县一带。

本文开展了对突发地震的预测研究，主要研究方向是将地震预测任务同地质统计学相结合，旨在研究出一种更为精准，更为全面，且数据获取不过多耗费人力的预测方式。随着研究内容的深入，发现本文研究依然存在着一些不足，期待后续的进一步分析与拓展，研究展望如下：

(1) 尝试用条件地质统计模拟的方法来扩充指示克里金的方法单一性，并采用更丰富的数据来进行研究。未来能够更全面有效地对发震概率进行估计。

(2) 根据地震目录数据计算地震预测因子，同时找到更多与地震发生相关的变量，从而训练出拟合度更好的模型。

参考文献