0 引言

跨断层形变数据分析时通常将两期观测值之差超过3倍中误差、较长时段出现增值(或减值)、断层活动突破近几年的周期变化等作为断层活动异常信息的确定原则，分别适合于提取突变、趋势转折等短临或中短期异常(龚复华等，1998)。周硕愚(1990)提出建立正常动态基准，以观测值偏离动态基值超过一定限度作为异常定量识别，适用于尖点、台阶等异常识别(陈兵等，2000)。郭良迁等(2004)利用华北地区的断层形变资料求解形变异常参数和应变累计率，提取中期和中长期前兆异常信息和应变信息，研究其区域断层形变对强震的响应能力。这些异常判定指标和方法在分析预报中起到了一定作用，但却具有很强的主观性和经验性，缺乏对方法系统性的评价检验。

跨断层形变观测获得的时间序列具有非线性、非平稳、多层次特征，增加了异常分析和检测的难度。前人在分析非平稳人类心率数据时提出了启发式分割算法(BG算法)，用于处理非线性、非平稳序列的均值异常点，并已在医学、气象和形变等领域得到一定的应用(Bernaola-Galván et al，2001；封国林等，2005；姚宜斌等，2019)。本文采用BG算法对唐山台跨断层形变观测时间序列进行均值突变检测，以此为基础讨论唐山台跨断层形变的阶段性变化特征，并进行预报效能评估。

1 资料概况

唐山地震台被1976年唐山7.8级地震后留下的SN向断裂穿过，地面有明显的水平和垂直错动。台站跨唐山大型断裂带东支断裂中的Ⅴ号断层(唐山断裂)，该断裂与其上部地表主破裂带产状一致，上下贯通，近乎直立(谢觉民等，1997)，断层倾角约72°(郝书检等，1998)。台站在大震出露的断裂两侧布设有5个水准、基线公用观测墩(图 1)。其中，3(Ⅰ)号、4(Ⅲ)号水准(基线)点位于断层上盘，1(Ⅳ)号、2(Ⅱ)号水准(基线)点位于断层下盘。基线和水准测线各4条，斜交测线J Ⅰ-Ⅳ(S3-1)与断层夹角约6°，J Ⅲ-Ⅱ(S4-2)与断层夹角约126°。

本文选取4条跨断层短基线测段(J Ⅰ-Ⅱ、J Ⅰ-Ⅳ、J Ⅲ-Ⅱ、J Ⅲ-Ⅳ)和4条跨断层短水准测段(S3-2、S3-1、S4-2、S4-1)1984年1月1日~2019年2月28日的形变观测资料进行分析。按照距台站100km以内M_S2.0~4.0和M_S4.0~5.0、200km以内M_S5.0~6.0、300km以内M_S6.0~7.0的原则选择地震，并参考相关文献(周海涛等，2009)进行补充。

2 数据处理方法 2.1 资料预处理

对原始观测资料预处理，包括剔除粗差和明显异常干扰项，对观测数据进行合并、线性插值等以获得连续等间隔观测时间序列。考虑到唐山台跨断层形变观测期次为每天一期，将同一天发生的多个小震计为一个发震日期。

2.2 断层运动分量计算

为探究断层运动状态与唐山台及其周围地区典型地震之间的关系，计算断层运动三分量累积参量(周海涛等，2009；张希等，2012；李文静，2014)，即走滑分量a、张压分量b、垂向分量h。

2.3 BG算法

启发式分割算法(BG算法)是将一个长度为N的跨断层形变观测时间序列x(t)，从左到右分别计算每个时刻i左右两边的平均值μ₁(i)、μ₂(i)及标准偏差s₁(i)、s₂(i)，定义其合并偏差为

$ {s_D}\left(i \right) = \left({\frac{1}{{{N_1}}} + \frac{1}{{{N_2}}}} \right) \times \sqrt {\frac{{\left({{N_1} - 1} \right) \times {s_1}{{\left(i \right)}^2} + \left({{N_2} - 1} \right) \times {s_2}{{\left(i \right)}^2}}}{{{N_1} + {N_2} - 2}}} $

(1)

其中，N₁、N₂分别为i左边和右边部分的时间序列长度。

使用t检验的统计量T(i)量化表示i时刻左右两部分均值的差异

$ T\left(i \right) = \left| {\frac{{{\mu _1}\left(i \right) - {\mu _2}\left(i \right)}}{{{s_D}\left(i \right)}}} \right| $

(2)

对x(t)中每个点重复上述计算过程，得到检验统计值序列T(t)，T越大则表示该时刻左右两部分均值相差越大。计算T(t)最大值T_max的统计显著性

$ P\left({{T_{\max }}} \right) = {\rm{Prod}}\left({T \le {T_{\max }}} \right) $

(3)

P(T_max)表示随机过程中取得T值小于T_max的概率。一般情况下

$ P\left({{T_{\max }}} \right) \approx {\left({1 - {I_{\frac{v}{{\left({v + T_{\max }^2} \right)}}}}\left({{\delta _v}, \delta } \right)} \right)^\eta } $

(4)

由蒙特卡洛模拟得到η=4.19lnN-11.54，δ=0.4，v=N-2，I_x(a, b)为不完全β函数。

设定一个临界值P₀，如果P(T_max)≥P₀，则于该点将x(t)分割为两段均值有一定差异的子序列，否则不分割。对新得到的2个子序列分别重复上述步骤，如果子序列有P(T_max)≥P₀，且子序列与其左右相邻的子序列间均值的差异程度满足上述条件，则对子序列进行分割，否则不分割。如此重复至所有子序列不可分割为止(子序列长度小于等于最小分割尺度l₀)。

通过上述步骤，将原时间序列分割为若干不同均值的子序列，各子序列分别代表不同层次的信息，分割点则为均值突变点。通过调整值P₀和l₀的大小可以改变检测的尺度和精度(封国林等，2005；龚志强等，2006)。

2.4 速率阶段性变化识别

本文将临界值P₀设为95%。为获得不同尺度的检测结果，经比较选择10年(l₀≈3650)、3年(l₀≈1095)和1个月(l₀≈30)分别作为最小分割尺度。采用BG算法对上述8个原始形变观测时间序列和3个断层运动分量时间序列(表 1)分别进行均值突变检测，将唐山台跨断层形变观测时间序列划分为不等长度的子序列，计算突变前后子序列均值变化量，探究唐山台跨断层形变的阶段性变化特征。

2.5 预报效能评估

将检测结果与台站及其周围地区典型震例的对应情况进行统计，发现当唐山台跨断层形变观测时间序列最小分割尺度设为10年、3年和1个月时，检测到的突变数量分别与M_S5.0及以上、M_S4.0~5.0和M_S2.0~4.0地震发震总次数相当，其对应的最优发震间隔约为5年、2年和3周，故以此作为预报有效时间，假设突变均为前兆异常，统计突变距离发震日的时间间隔，计算虚报率a、漏报率b和R值(罗兰格，2015)，并据此对唐山台跨断层形变观测时间序列进行预报效能评估(表 1)。计算虚报率a、漏报率b和R值的计算公式分别为

$ a = \frac{{n_0^1}}{{{N^1}}} $

(5)

$ b = \frac{{n_1^0}}{{{N^1} + n_1^0}} $

(6)

$ R = 1 - a - b - p $

(7)

其中，n₀¹为预报有震而未发生地震次数；n₁⁰为预报无震而发生地震次数；N¹为预报有震总次数。p为n次预报均发生地震的自然概率，p值通常很小，可以忽略不计。R＞0表示预报有效，且R值越大效果越好，最大值为1；R=0表示预报无意义；R＜0则预报无效。

3 计算结果分析

由断层运动三分量累积参量(图 2(b))可以看出，唐山台所跨断层总体呈右旋压性逆断运动；1996年断层左旋逆转，压性逆断活动减弱；2003年恢复右旋，张性转折；2013年至今断层呈微弱左旋压性运动，逆断活动近年来有所加速。

3.1 10年尺度检测结果

最小分割尺度为10年时，采用BG算法检测到的均值突变次数(4~6次)与唐山台及其周围地区M_S5.0及以上地震发震次数(4次)相当。从突变前后子序列均值变化情况看，基线和水准变化幅度相当，自0.14mm至3.61mm不等(图 3(a))。基线变化基本可以分为4个阶段，1984~1998年以伸长为主，突变前后均值变化量为0.55~2.20mm；1999~2008年以缩短为主，均值变化量为0.41~1.18mm，明显低于上一阶段；2008~2012年恢复以伸长为主，均值变化量略小于第一阶段，为0.28~1.95mm；2012年至今以更低速率缩短，均值变化量仅为0.18~0.33mm。水准曲线自观测以来一直呈下降趋势，突变前后均值变化量为0.51~3.61mm；但2000~2004年各测线存在陆续小幅上升(0.14~0.58mm)。结合典型震例来看，台站及其周围地区M_S5.0及以上地震发生前5年内，唐山台基线和水准一般各出现2项或2项以上突变。

断层运动三分量突变前后均值变化情况也显示唐山台所跨断层以右旋压性逆断活动为主(图 2(b)、3(b))，其走滑、张压和垂向运动分量的均值变化区间依次增大，分别为0.34~0.8mm、0.23~1.42mm和0.58~2.82mm。但是，2000~2004年断层左旋张性正断突变，与水准观测同期的变化相对应。其中，水平面走滑和张压分量的变化量较大，分别为0.97mm和0.99mm；垂向分量变化量较小，仅为0.18mm，趋势近乎水平。说明在此期间唐山断裂周边应力场发生了重大转折(周海涛等，2009)，可能与2006年河北文安M_S5.2地震有关。

定量统计结果(表 1)显示，各测项突变后发生地震的概率为50%~100%。突变距发震日的时间间隔自51天至2274天不等，最短为水准测线(S3-1)1991年4月9日突变后51天，即1991年5月30日，河北唐山发生M_S5.1地震(图 3(a))。

以5年作为预报有效时间，计算得到各测项的虚报率和漏报率分别在0~50.0%和0~33.3%之间，R值评分介于0.17~1.00之间。

在基线原始形变观测时间序列中，J Ⅲ-Ⅳ虚报率和漏报率最低，R值评分最高；J Ⅰ-Ⅳ虚报率和漏报率最高，R值评分最低。在水准原始形变观测时间序列中，S4-1虚报率和漏报率最低，R值评分最高；S3-2和S4-2的R值评分最低。结合测线布设情况(图 1)可以看出，J Ⅲ-Ⅳ(S4-1)对M_S5.0及以上地震的预报效能最好，J Ⅰ-Ⅳ(S3-1)的预报效能最差。

在断层运动三分量中，垂向分量的R值评分最高；走滑分量的虚报率和漏报率最高，R值评分最低。走滑分量对M_S5.0及以上地震的预报效能明显低于垂向分量和张压分量。

3.2 3年尺度检测结果

最小分割尺度设为3年检测到的突变次数(17~21次)与M_S4.0~5.0地震的发震次数(20次)相当，对应较好(表 1)。基线和水准原始形变观测时间序列突变检测结果的阶段性特征与10年尺度上的结果一致。其中，基线突变前后子序列的均值变化量为0.05~1.49mm(图 4(a))，大致可分为4个阶段，1984~1998年以伸长为主，均值变化量为0.1~1.36mm；1999~2008年伸缩交替，均值变化量为0.13~1.49mm；2008~2011年恢复以伸长为主，均值变化量为0.23~0.75mm；2012年至今伸缩交替，均值变化量为0.05~0.79mm，幅度明显减小。水准均值变化量略大于基线，为0.07~2.33mm，曲线整体呈下降趋势；但2000~2004年各测线曾陆续小幅上升(0.41~0.67mm)。结合震例来看，台站及其周围地区M_S4.0~5.0地震发生前2年内，一般先后出现2项以上的基线或水准突变。

由断层走滑分量(图 2(b)、3(b)、4(b))可以看出，1984~1996年断层以右旋运动为主，突变前后均值变化量为0.16~0.78mm。其中，1990年前后断层短期内左旋转向，可能与1990年唐山M_S4.5地震有关，1991年恢复右旋后发生唐山M_S4.8和M_S5.1地震；1997~2002年以左旋为主，均值变化量为0.37~0.81mm；2003~2013年继续以右旋运动为主，均值变化量为0.29~0.63mm。但是，2004年和2005年断层曾出现2次左旋转向突变，随后2006年发生河北文安M_S5.2地震。2013年以来，均值变化量为0.05~0.18mm，左旋走滑运动微弱。

张压分量显示1984~2002年断层以压性活动为主，突变前后均值变化量为0.1~0.67mm。1985年、1994年和2001年出现张性突变，其中1994年突变前后均值变化量为0.52mm，可能与1995年唐山M_S5.1地震有关，另外2次突变均值变化量较小，分别为0.09mm和0.17mm。2002年开始张压交替，断层整体呈弱压性运动。其中，压性突变均值变化量为0.08~0.72mm；张性突变变化量为0.13~1.26mm，且幅值依次递减。另外，统计发现2003~2016年张压交替突变时间间隔约为2年，自17个月至35个月不等，与张跃刚等(2013)对唐山余震区2004年以来M_L≥4.0地震的准周期性研究结论基本一致。

垂向分量显示1984~2002年断层上盘相对上升，以逆断为主，突变前后均值变化量为0.27~1.94mm。1997年出现1次小幅正断突变，均值变化量为0.28mm，可能与1998年唐山M_S4.7地震有关。2002~2007年垂向分量曲线起伏变化，总体趋势较为平稳。其中，逆断突变均值变化量为0.35~0.47mm，正断突变均值变化量为0.30~0.69mm。2010年恢复逆断，均值变化量为0.18~1.12mm，且近年来活动速率有所增加。整体而言，垂向分量突变检测结果与同时段张压分量显示的断层运动趋势基本一致。

定量统计(表 1)发现，各测项突变后发生地震的概率为50%~67%。突变距发震日的时间间隔自8天至1048天不等，最短(S3-2)为1985年5月14日突变后8天，即1985年5月22日，在距唐山台36.7km处发生了M_S4.3地震(图 4(a))。

以2年作为有效时间，计算得虚报率和漏报率分别为41.2%~60.0%和30.0%~37.5%，R值评分介于0.03~0.26之间。在基线原始形变观测时间序列中，J Ⅰ-Ⅱ的R值评分最高；J Ⅲ-Ⅱ虚报率和漏报率最高，R值评分最低。在水准原始形变观测时间序列中，S3-2虚报率和漏报率最低，R值评分最高；S3-1的R值评分最低。结合测线布设情况(图 1)可以看出，J Ⅰ-Ⅱ(S3-2)对M_S4.0~5.0地震的预报效能最好，J Ⅲ-Ⅱ(S4-2)的预报效能最差，且水准的预报效能整体上优于基线。

在断层运动三分量中，张压分量的R值评分最高；走滑分量的虚报率和漏报率最高，R值评分最低。因此，可以认为走滑分量的预报效能最差。

3.3 1个月尺度检测结果

最小分割尺度设为1个月检测到的突变次数(667~701次)与M_S2.0~4.0地震的发震次数(730次)相当。据统计(表 1)，各测项突变后发生地震的概率为55.6%~59.8%。突变距发震日的时间间隔自0(当天)至34天不等。突变当天发震次数占全部对应次数的6.7%~11.4%，1周内发震次数占全部对应次数的51.4%~65%，2周内发震次数占全部对应次数80%以上(图 5)。也就是说，在1个月尺度上进行突变检测，突变后有50%以上的概率会在唐山台及其周围100km范围内发生M_S2.0~4.0小震，但如果突变后2周内没有地震发生，基本可以认为这是一次虚报。

以3周作为预报有效时间，计算得虚报率和漏报率分别为41.6%~46.2%和32.4%~34.2%，R值评分介于0.20~0.26之间。在基线原始形变观测时间序列中，J Ⅰ-Ⅱ的R值评分最高；J Ⅲ-Ⅱ的R值评分最低。在水准原始形变观测时间序列中，S4-2的R值评分最高；S4-1的R值评分最低。考虑突变距离发震日的时间间隔，结合测线布设情况(图 1)，认为J Ⅰ-Ⅱ(S3-2)对M_S2.0~4.0地震的预报效能最好，水准的预报效能整体优于基线。在断层运动三分量中，垂向分量的R值评分最高，预报效能最好。

4 结论与讨论

本文采用BG算法对唐山台跨断层形变观测时间序列进行突变检测，减少了跨断层形变异常分析中的主观性和经验性因素的影响。根据检测结果，计算突变前后子序列均值变化情况，探讨唐山台跨断层形变的阶段性变化特征，并结合台站及其周围地震典型地震进行预报效能评估。结果显示：

(1) 唐山台跨断层形变阶段性变化特征显著。断层运动三分量表明，唐山台所跨断层1984年以来整体呈右旋压性变化趋势；1996年、2002年和2013年前后断层走滑运动转向；2002年和2010年前后断层张压运动分段特征明显，体现了唐山台所跨断裂及其周边地区应力场在上述时间段的重大变化。

(2) 原始形变观测时间序列和断层三分量形变时间序列在最小分割尺度设为10年、3年和1个月时，检测到的突变次数分别与M_S5.0及以上、M_S4.0~5.0和M_S2.0~4.0地震发震次数相当，各测项突变后对应发生地震的概率均不低于50%，突变距发震日的最短时间间隔自当天到2个月不等。

(3) 原始形变观测时间序列和断层三分量时间序列的预报效能相当。从单条测线角度进行比较，J Ⅲ-Ⅳ(S4-1)与M_S5.0及以上地震的对应情况较好，J Ⅰ-Ⅱ(S3-2)与M_S4.0~5.0和M_S2.0~4.0地震的对应情况较好；从观测手段角度进行比较，水准的预报效能整体上优于基线；从断层运动特征角度进行比较，垂向分量和张压分量的预报效能整体上优于走滑分量。

本文假设检测到的突变均为前兆异常，是对唐山跨断层形变观测资料进行的回溯性的预报效能评估，结果可供唐山地区后续地震预报研究参考。需要指出的是：①BG算法基于整个时间序列计算突变，当时间序列积累延长，突变检测结果可能发生改变。且由于唐山台及其附近M_S5.0及以上地震和M_S4.0~5.0地震数量较少，本文相关的计算结果仅作为参考，不具有统计学意义。②BG算法可以较好地检测均值变异点，但跨断层形变的突变表现形式复杂多样，因此，可以考虑与窗函数、滤波、距平等方法结合，提高异常检测和识别能力。③在数据积累充足的情况下，该方法理论上可以应用于任一时间序列的均值突变检测，但存在一定滞后性，要真正做到数据实时处理和前兆异常提取，还有很多问题需要解决。④BG算法基于一定的数学理论基础，应用于实际的跨断层形变时间序列检测，尚缺乏物理机理层面的理论支撑，且地震的孕育发生过程复杂，跨断层形变表现形式多样，如何结合具体情况进行算法的优化改进，仍需进一步研究。