新疆北天山地区破年变异常识别及其效能检验①

颜 龙， 史勇军， 闫 玮， 雷 晴， 朱治国

(新疆维吾尔自治区地震局，新疆 乌鲁木齐 830011)

地震孕育过程要经历弹性变形和非弹性变形两个重要阶段，其中非弹性变形可能引起岩石物理性质的显著变化以及伴随地壳形变、地下水等非线性、非均匀的变化[1-3]。这些低频非线性过程大都受气压负载、地下水负载、地表气温以及地壳基岩热胀冷缩等物理效应的共同作用，尤其以季节冷热变化带来的周期成分尤为显著，并表现为年周期震荡运动。已有研究表明，地震发生前不仅可能有短临异常出现，当应力积累到一定程度，也可能会出现大量破年变异常，这些破年变异常在地震预测过程中可能具有中长期尺度的指示意义[4]。例如，米亚齐金等利用弹性空间孕震模型分析前兆异常在孕震空间尺度与震级之间的关系[5]，认为地震发生前不仅会有短临异常出现，当应力积累到一定程度，也会出现趋势性转折和破年变异常。刘耀炜等[6]提出用月均值计算月变差率方法来定量描述破年变异常,以此方法确定本年相对上年在同一月的变化程度。刘冬英等[7]针对大华北地区地下流体资料的转折时间点和形态特征进行分析，并指出异常集中分布的区域可能是当前或今后中强地震发生的主体区域。以往研究中，多是通过对数据曲线形态的经验辨别，将不同学科、不同手段的观测资料进行联合对比[8-10]。但不同地震前的异常测项不同，各测项的异常形态、幅度和持续时间也不尽相同，异常的判断存在明显的主观性。本文中拟对新疆北天山地区具有年动态的定点形变资料和地下流体资料进行系统梳理，开展破年变异常定量识别研究，提取出异常信息，并在此基础上分析新疆北天山地区破年变异常的时空演化特征，为客观识别出地震异常信息提供参考依据。

1 研究区域概况

新疆地球物理观测台站主要沿南北天山地震带布设，在昆仑山、阿尔泰地震带上也分布有个别台站，以北天山地震带中段和南天山地震带中、西段分布最为密集，北天山地震带历史上中强震多发，是新疆地震监测的重中之重，该区域主要以形变和地下流体等多项监测为主，广泛分布于乌鲁木齐—奎屯—温泉线(图1)。

图1 新疆北天山地区地球物理观测点分布图Fig.1 Distribution map of precursor observation point in North Tianshan area，Xinjiang

为了进一步定量分析北天山地区地球物理观测资料中具有年变异常的测项及其预报效能，本文中基于地球物理观测资料预报效能评估结果，选取预报效能B级以上的资料，包括5个台站的10个定点形变和10个台站的10个地下流体测项。其中A类资料9项(水平摆3项、水位2项、水氡3项、流量1项)，B类资料11项(水管仪3项、洞体应变3项、水平摆1项、水位2项、水氡1项，流量1项)(表1)。

表1 新疆北天山地区地球物理观测测项预报效能评估表

2 研究方法

根据不同地球物理观测资料变化特征选取适用的年变拟合方法，即制作标准年变模型，然后利用R值方法进行预报效能检验，提取异常阈值。首先采用破年变异常自动识别软件ADDv1.exe，该软件主要利用4种数学模型(矩平法、傅里叶滑动法、最小二乘法、分段最小二乘法)对具有年动态的观测数据制作 “标准年变”模型，即年变拟合曲线，以此获取观测数据清晰可识别的年变形态，然后根据拟合曲线和观测数据的残差判断异常，即残差序列中高于阈值部分，再将破年变异常的测项进行R值检验，提取异常指标。

R值评分方法由许绍燮提出[11]，该方法在时间序列中的研究较多，20世纪90年代以来，该方法多被用于对历年年度预测效能的检验，之后又被多次完善[12-14]，在检验前兆异常和地震预报效能得到了普遍应用和认可。

(1)

式中，当R>0表示预报有效，R值越大效果越佳，R=1表示预测完全正确，所有地震全部落入危险区，且不存在任何虚报。通过对一段时间不同目标震级地震进行若干次R值动态检验，可以得到各测项异常的最佳阈值。

3 破年变异常检验示例

根据表1中新疆北天山地区地球物理观测测项，本文中将分2个震级档分别进行检验。测点震中距300 km范围内的MS5.0地震和震中距500 km范围内的MS≥6.0地震，各地球物理观测数据均截止到2019年12月31日。

3.1 阜康洞体应变EW分量(矩平法)

阜康洞体应变自观测以来积累近10年的观测资料，其NS分量无显著年变动态，而EW分量年变形态较显著。图2(a)对阜康洞体应变EW分量进行去趋势处理，黑色实线为EW分量预处理后的观测曲线，红色虚线为通过矩平法计算出的年变拟合曲线。而图2(d)中，色标差异则代表预报效能大小，黑色圆点代表对应最佳阈值和预测时间的最优R值，通过对测点周边MS≥6.0地震进行检验，可以看出阜康洞体应变EW分量异常出现后102 d为优势对应时间，最佳阈值为1.2，对应的预报效能R=0.418 3。图2(c)为观测值和拟合值的残差绝对值曲线，红色虚线代表残差序列所带入阈值，从全部筛选的MS≥6.0地震检验结果来看，检验效果不是很理想，仅报准2016年12月8日呼图壁MS6.2地震和2017年8月9日精河MS6.6地震[15-16]，其原因可能是阜康洞体应变EW分量长时间受天池库容变化引起的载荷变化影响，剔除明确的干扰时段后，在有效的预报时间内仅发生了4次MS≥6.0地震，该资料积累时间较短也是预报效能较低的可能原因之一。

3.2 温泉水平摆倾斜EW分量(傅里叶滑动法)

温泉水平摆倾斜EW分量数据存在11个完整的年变，年变振幅大小虽有小幅变化，但周期形态较为稳定[17]。通过傅里叶滑动方法拟合标准年变选取周期T=365 d，采用3T窗长滑动，拟合结果如图3(a)，可以看出该拟合结果和原始曲线较贴合。利用R-TT曲线(效能-预测时间/阈值)分析发现温泉水平摆倾斜EW分量对应周边300 km内的MS5.0地震具有较好的映震效能，但无优势集中区域，地震分布如图3(b)。图3(d)的检验结果中，当阈值取1.5时，预报优势时效为82 d，对应的最佳预报效能R=0.802 5，通过了预报效能检验。图3(c)中红色虚线上的绿色方块代表超过残差阈值的异常对应时段，可以看出温泉水平摆倾斜EW分量破年变异常总共出现了9次异常时段，在此期间总共发生了9次MS5.0地震，报准了其中的5次地震，说明温泉水平摆倾斜EW分量在出现破年变异常后发生地震的可能性较大。

图2 阜康洞体应变EW分量年变识别结果(a) 年变曲线拟合图 (b) 震例空间分布图 (c) 异常自动识别结果 (d) R-TT曲线Fig.2 Annual change anomaly recognition results of EW component of cavern strain in Fukang

图3 温泉水平摆倾斜EW分量年变识别结果(a) 年变曲线拟合图 (b) 震例空间分布图 (c) 异常自动识别结果 (d) R-TT曲线Fig.3 Annual change anomaly recognition results of EW component of horizontal pendulum tilt in Wenquan

3.3 榆树沟洞体应变NS分量(最小二乘法)

榆树沟洞体应变NS分量数据累计超过12年，且年变幅周期和振幅均无明显规律。由于该点NS分量数据年变无显著的分段点，因此利用最小二乘法和分段最小二乘法年变拟合后效果一致。利用年变拟合残差对该测项500 km内的MS6.0地震进行检验如图4(d)，R-TT曲线显示，当阈值取0.6时，预报优势时效为62 d内，对应的最佳预报效能R=0.4612。从图4(c)中看出，资料的异常时间段较多，有效预报时间段内发生28次MS6.0地震，但仅报准4次，并且这4次报准的地震并非发生在最大异常之后。从空间分布上来看，4次地震都在该测点西侧且无优势区域，该测项的破年变异常的信度较低。

图4 榆树沟洞体应变NS分量年变识别结果(a) 年变曲线拟合图 (b) 震例空间分布图 (c) 异常自动识别结果 (d) R-TT曲线Fig.4 Annual change anomaly recognition results of NS component of cavern strain in Yushugou

3.4 新10井水位(分段最小二乘法)

新10井水位自观测以来无明显年变规律。2000年以来的变化量均不相同，多处数据还存在相位改变的的情况[18]，但在2011年前后存在年变幅差别，2011年后整体年变的振幅明显增大。根据新10井水位的变化特征，可采用分段曲线拟合，在每段区间上进行局部最小二乘拟合，然后通过F检验判断分段点的显著性，在显著性水平α>0.2的水平时，即为分段点。图5(a)可以看出，年变拟合曲线较好的适应原始数据2011年前后的年变幅变化，分段处理可有效减少拟合曲线与原始数据的残差序列。利用残差序列对区域300 km范围内的MS5.0地震进行预报效能检验如图5(d)所示，预报效能最高为R=0.592 7，阈值为2.2。有效预报时间段内总共发生10次MS5.0地震，报准7次，291 d内的预测效果最好。综合考虑整体预报效能和周边地震分布特征，新10井水位破年变异常发生后对周边MS5.0地震的预报效能较好。

图5 新10井水位年变识别结果(a) 年变曲线拟合图 (b) 震例空间分布图 (c) 异常自动识别结果 (d) R-TT曲线Fig.5 Annual change anomaly recognition results of Xin 10 well

4 检验结果

参照上述示例，分别利用4种拟合方法对表1中的测项进行年变拟合，并利用残差序列定量化提取异常指标：对应时间、预报效能、最佳阈值(表2)，主要包括北天山区域10个定点形变的20分量测项和10个地下流体测项，其中阜康洞体应变NS分量由于年变不显著，暂不做分析。

按照对应的空间(天山以北、周边)划分，地震优势对应空间为天山以北范围的测项有23项(占研究测项总数的76.7%)，无优势对应区域的7项(占研究测项总数的23.3%)。按照对应的时间划分，优势对应时间段为短期的测项有22项(占研究测项总数的73.3%)，优势对应时间段为中期的6项(占研究测项总数的20.0%)，优势对应时间段为长期的2项(占研究测项总数的6.7%)。按照对应的震级划分，破年变异常优势对应为MS≥6.0地震的测项有22项(占研究测项总数的73.3%)，优势对应为MS5.0地震的8项(占研究测项总数的26.7%)。按照测项所适应的拟合方法划分，适应矩平法拟合的测项有3项(占研究测项总数的10.0%)，适应傅里叶滑动法拟合的17项(占研究测项总数的56.7%)，适应最小二乘法拟合的4项(占研究测项总数的13.3%)，适应分段最小二乘法拟合的6项(占研究测项总数的20.0%)。从整体研究结果看，多数地球物理观测测项对北天山地区中强震预测效果较好，对应时间以中短期居多。定点形变测项较适应选择最小二乘法和分段最小二乘法进行年变拟合，地下流体测项较适用傅里叶滑动法进行年变拟合。

表2 新疆北天山地区地球物理观测测项破年变异常预报效能评估与预报指标

5 结论与讨论

破年变异常作为地震预测中常见的异常形态，如果只利用单学科单手段进行预报，往往可能出现虚报漏报，开展多学科破年变异常特征研究可以更好地为各前兆异常信度进行相互佐证，大大消除干扰，降低前兆异常判定的复杂度。通过对北天山区域的定点形变和地下流体中具有年动态的观测资料进行破年变异常分析，选用4种不同的拟合方法进行年变周期拟合，得到正常年变时间拟合序列，从分级震例统计角度，得到最佳预测效能的参数组合，即异常最佳阈值和预测时间，并尝试探讨该区域地球物理观测资料破年变异常的时空特征。

(1) 在识别破年变异常过程中，针对地球物理资料各自变化特征需要选取不同年变拟合方法。矩平法适合积累时间较短年变形态变化相对较少的资料，此方法对大多已积累10年以上年变形态复杂多样的观测资料来说，拟合度相对过低，得到残差序列取值过大，且不利于异常的提取。最小二乘拟合法对年变周期比较规律的资料拟合度较高，但对于振幅多变的观测资料拟合度就较差，若振幅变化存在明显前后分段现象，则可使用分段最小二乘法寻找出显著分段点，这样得到的拟合结果更接近于原始观测数据。傅里叶滑动法在观测数据较长年变形态比较复杂且存在非规律的振幅改变时拟合效果最佳，一般选取1年为周期，3年窗长进行滑动，这样可以较好避免某年年变幅异常带来的拟合不稳定。

(2) 从异常对应时间来看，北天山地区地球物理定点观测资料的破年变异常大多为中短期异常，这与以往分析人员利用数据曲线形态的经验辨别的破年变异常震例结果较为一致，即破年变异常对应的地震以半年内的居多。说明破年变异常一定程度上反映了区域构造应力状态与异常所处区域整体动力学背景的调整和演化特征，此类异常对区域内中强地震的对应性较好。

(3) 从对应空间上来看，主要优势对应区域为北天山地区。其原因是所选地球物理观测测项主要沿北天山地震构造带分布，该区域中强地震频发。因此，对北天山区域地震检验效果较好，虽然部分南天山地震也在所选空间范围内，但从检验结果上来看，对南天山区域的预报效能较差。

北天山地区沿构造具有明显的分带性，强震构造带之间力学性质不尽相同，不同测项对不同构造带上的地震预报效能也各不相同。定量化破年变异常识别不仅能直观的对各测项整体预报效能进行评价，还能克服人为因素影响，更为客观和科学的显示定量化结果，具有一定优势，但单从震级和震中距两个方面作为约束条件，不同的测项和不同方向分量对周边地震的反应都存在其特异性，需要结合具体的地质构造背景等进一步深入研究。