LURR方法是通过潮汐应力在某一时刻在断层面上引起的库仑应力的增加和减少,判断加载和卸载,使用地震目录计算加卸载响应比值,预测地震迫近程度(尹祥础等,1994)。当系统处于稳态时,加载响应与卸载响应相当,因而LURR=1; 当系统偏离稳态时,加载响应大于卸载响应,LURR>1。

依据 LURR的计算规则,编制完成可在Linux下运行的LURR模块(图2)Python程序。以地震目录为基础,利用LURR方法可执行程序,产出LURR时序曲线和异常的空间分布。可实现未来数天、数周、数月和数年多时间尺度的地震预测。

图2 LURR计算模块
Fig.2 The calculation module of LURR

Chen 等(2020)在国内外多个强震前,通过地震仪资料和GPS资料,发现地壳振动信号在特定频段出现显著增强的现象(图3a)。图3b中,红色为异常时段的分析结果,即震前5～10 d的频谱,蓝色为背景时段的分析结果,即震前2个月的频谱。由图可见,在10^-5～10^-4 Hz的频带范围内,震前5～10 d振幅比背景振幅明显增强(于晨等,2021)。

图3 地壳振动方法定位过程
Fig.3 Positioning process of the crustal vibration method

本模块以识别0.2 Hz、0.05 Hz和0.0 001 Hz 3个频段(可根据实际应用情况进行调节)地壳振动信号为目的进行开发,当增强幅度超过背景值20%时,则认为记录到地壳振动信号。采用Tanimoto等(2006)提出的单台连续波形地震信号源定位法,通过对多个台站的单台法定位结果综合分析,判定地震发生位置。将水平方向的2个信号投射到正交轴上,并按照顺时针旋转叠加,确定最大振幅对应的方向为信号源方位。以信号源方向为中心,构造角度为30°、半径为300 km的扇形区域中,当至少3个台站出现重叠区域时,该区域为潜在的地震源(图3c)。

按照上述思路,通过设定背景观测时段和变化时段,产出异常分布图和异常来源分布图。目前该模型已经实现川滇地区地震危险性的实时跟踪。

态矢量概念来源于统计物理学,是一种对连续场进行粗粒化描述的方法。尹祥础等(2004)引入态矢量方法对1976 年唐山大地震、1975年海城大地震等进行了研究,发现大地震前态矢量参数有显著的变化,因而认为态矢量的特征变化可能是一种地震前兆。

将研究区域剖分成n个子区域(图4),通过子区域内小震活动在一定时期内的变化预测地震。具体方法为:将子区域i内对应的地震活动参量V_i看做n维矢量的一个分量,可以得到任意t时刻的n维态矢量:

图4 研究区域的n个子区域划分
Fig.4 Division of the study area in n sub-regions

V_t=[v₁(t),v₂(t),…,v_n(t)] (1)

其增量为:

ΔV_t-Δt,t=V_t-V_t-Δt (2)

分别计算态矢量的模:

M=|V_t| (3)

V_t-Δt和V之间的夹角为:

增量的模为:

ΔM=|ΔV_t-Δt,t|(5)

态矢量V_t与单位矢量V_e之间的转角关系如下:

式中:单位矢量V_e表示各个分量相同的矢量。

本模块以地震目录数据为基础,产出态矢量时序曲线和对应数据。可实现未来数天、数周、数月和数年多时间尺度的地震预测。

地震综合概率模块是为了整合各种预测方法的预测结论,以风险概率的形式给出危险区域的空间分布和危险程度,推进地震综合预报由定性预报向定量预报的转变。它兼顾地震背景发生概率和每种预测方法的历史预测效能,采用贝叶斯公式计算综合概率(Field,2007):

P=1-(1-P₁)(1-P₂)(1-P₃)… (7)

式中:P为综合概率; P₁、P₂、P₃…为通过不同方法计算得到的对同一地点的地震概率。

在实际预测时,将各种方法评价所得的R值评分作为权重,综合概率P可改写为:

式中:N表示所选用的预测方法的个数; R_i是第i种方法的R值; c定义为:

式中: 是平均概率; 是平均R值。

本模块以地震目录和各种算法的预测结果为输入,能够产出地震综合概率预测分布图和对应网格概率值,实现不同时间尺度的向前地震预测功能。

按照PI方法(Rundle et al,2003),地震活动可以作为受稳恒速度连续驱动的阈值系统,系统中的能量耗散在长期平均的意义上接近稳定。PI算法假定中小地震活动异常是中等以上地震的前兆现象,考虑地震活动的增强与平静。具体步骤为:①对研究区进行空间网格划分,以落入每一网格中的地震事件个数构建时间序列; ②给定滑动时间窗,通过时空变化计算每个网格未来发生显著地震的可能概率,再减去地震背景概率,将发震概率高的网格以深浅不一的暖色块的形式呈现出来,即热点分布图。

PI方法使用Matlab语言编写,可在Linux下运行。以地震目录为输入数据,产出图像信息空间分布图,实现未来数天、数周、数月和数年多时间尺度的地震预测。

ETAS模型描述了地震活动时空丛集的自激发点过程(蒋长胜,庄建仓,2010),可以应用于研究区域地震活动时空特征和余震活动特征分析。

按照ETAS模型,每次地震都会触发余震且余震活动在时间域内遵循Omori-Utsu公式,单位时间t内的余震发生次数为:

n_i(t)=Kexp[α(M_i-M_C)]/(t-t_i+c)^p (10)

式中:M_C为下限震级; M_i、t_i分别为第i次地震的震级和发震时刻; K、α、c、p是常量。t时刻的地震发生率为:

本模块输入的是地震目录数据,可以产出未来一定时间、空间、震级窗内的预测结果。

RI是基于统计学原理构造的地震预测模型,通过学习历史地震数据,预测未来将要发生的给定震级范围内的地震次数。将研究区域按经纬度划分为规则网格,以网格为单位统计给定震级条件下的历史地震次数,根据已发生的地震次数预测未来地震的发生次数,最后累加得到研究区域的预测值(Rundle et al,2000)。通过遍历搜索适应研究区地震活动特征的最优模型参数,进行地震危险性预测,绘制RI预测热点分布图。

本模块以地震目录数据为输入,产出RI热点空间分布图,可实现未来数天、数周、数月和数年多时间尺度的地震预测。

SRI基于统计学模型,定量识别地震活动显著增强与减弱两种典型异常。计算方法为:分别给定背景时段:T_b=t₂-t₀,研究时段:T_w=t₂-t₁。根据背景时段T_b内发生的地震个数N,给定T_w时间内的平均地震发生率:

根据泊松分布,T_w时间内发生k次地震的概率为:

根据T_w内实际发生的地震数n,定义地震发生率指数SRI为:

本模块以地震目录数据为输入,可以产出未来数天、数周、数月和数年SRI预测的异常空间分布图。

根据Gutenberg-Richter定律,b值大小与应力水平呈反比,可以用来指示区域地震危险性水平,且一般选取低b值区域作为发生地震的危险区域(易桂喜等,2008)。

b值用最大似然法估算为:

式中: 为计算选用目录的平均震级; M_C为地震目录的最小完备震级; dM为地震记录的最小标度,一般为0.1,通常地震总数N大于30。

在实际预测中,对研究区域进行网格划分,然后依次以各个网格节点为中心,筛选给定边长的区域内、给定时段内的地震目录,计算各节点处的b值。本模块以地震目录数据为输入,可以产出不同时段b值异常的预测结果。

小震调制比是地震学预测中的常用方法之一。当构造应力接近临界状态时,地震活动会受到固体潮调制(韩颜颜等,2017)。

本模块首先将研究区域进行网格化,依次以各个网格节点为中心,筛选一定半径、时段内的地震目录,按照农历初一和初二(朔),初七至初九(上弦),十五至十七(望),廿二至廿五(下弦)作为调制时间,计算相应节点的调制比,产出不同时段小震调制比异常的空间分布图。

R值评分是我国地震工作者普遍使用的地震预报效能评价方法(许绍燮,1989; Yu et al,2022),

计算公式为:

式中:S_P和S_T分别表示预测区域面积和整个研究区域面积; H和N分别表示一定时间窗内命中的地震数和地震总数。

当R>0时,表示预报有效,且R值越大效果越好,最大值为1; R=0表示预报无意义; R<0表示预报无效。结合统计分布检验,利用报准、漏报地震的个数,可以计算出大于97.5%置信度的临界值R₀,当R>R₀时,其预测的有效性高于自然发生概率,该预测结果具有较高的效能(图5a)。

图5 Molchan图表检验原理及其与R值评分的关系
Fig.5 The principle of Molchan chart test and its relationship with the R-score

R值评分类似于Molchan图表法的一种特例,表达为:

L=1-τ-ν (17)

式中:τ和n分别表示误报率和漏报率。

对于空间占有的预测方法,异常识别阈值降低,预测区域面积增大,同时命中率增加,其变化即Molchan曲线(Molchan et al,2008)。如图5b所示,图中任意点P₀(τ₀,m₀),由于其与直线τ+m-1=0,即直线R=0的距离,结合R值计算公式可推导出,Molchan图表中每个点与直线τ+m-1=0的距离D与R值评分成正比,二者可以相互转换。以直线τ+m-1=0为分界线,其左下方区域的R>0,右上方区域的R<0。

假设地震的发生是相互独立的,一定时段内发生的地震个数符合n重伯努利分布(泊松分布的离散化表达),该时段内单个地震的发生概率用预测区域的空间占有率表示,相应的显著性水平满足以下公式:

式中:N表示地震总数; H表示报准地震的个数; α表示置信度。例如取α=5%,根据τ值反算伯努利分布下的报准率h,进而得出Molchan图表的显著性曲线。该方法还定义了概率增益参数(Zechar,Jordan,2008):

N-test是以目标地震数目与实发地震数是否一致的统计检验方法(Schorlemmer et al,2007)。根据泊松分布、泊松累积分布函数,即发生不超过ω个事件的概率,可表示为:

使用评分变量δ₁和δ₂来定量表达发生“至少”和“不超过”N_obs个地震事件的概率,计算公式为:

当δ₁<α_eff时,表明预测的地震数目过少; 当δ₂<α_eff时,表明预测的地震数目过多,其中α_eff=0.5α。

ROC检验与R值评分有共同之处,在表现形式上给出不同危险概率阈值下各命中率与失败率的比值来描述预测的优劣程度(Holliday et al,2005)。变量a表示报准地震的格点个数,即预报有震的格点个数; 变量b表示误报地震的格点个数,即预报无震的格点个数; 变量c表示有地震无异常的格点数,即未报有震的格点个数; 变量d表示无地震无异常的格点数,即未报无震的格点数。则a+c表示总地震个数,b+d表示无地震的格点个数。定义报准率H=a/[KG0](a+c),虚报率F=b/[KG0](b+d)。定义预测ROC折线与随机预测直线所包络的面积E_f为有效预测系数,E_f值越大,表明ROC曲线越靠近图的左上部分,即在相同的报准率下付出的“虚报率”越低,预测效果越好,反之预测效果越差(图6)。

图6 ROC评价方法原理
Fig.6 The principle of ROC evaluation method