钻孔应变观测技术的发展历史已经有40多年，自从1966年邢台发生两次大地震以后，国家充分意识到了大震的严重危害。由著名的地质学家李四光先生牵头成立的地震大队，倡导群防群测预防大震，之后各种预报方法不断推出，期间观测地形变应力的预报方法也进入了大家视野，并发挥着重要作用^[1-3]。经过众多学者专家的努力探索研究，陆续成功研制出了RZB-1型电容式钻孔应变仪、RZB-2型钻孔分量式应变仪、TJ-1型体积式应变仪、YRY-4型分量钻孔应变仪。目前，我国已建成143个形变观测台站，安装了150余套体积式和分量式钻孔应力应变仪，构成了应用于地震观测的地应变台网^[4]。

地震是地壳应力积累到一定程度的爆发，钻孔应变观测方法在间接反应地壳应力变化方面有一定的优势。由于观测仪器安装在深入地下100～200 m的深井中，这很好的消除了近地表各种非构造因素的影响，得到地壳应力变动的清晰信息^[5-6]。持续观测更是得到了台站附近区域形变应力场的实时变化情况。而如何将钻孔应变的的资料和优势应用于地震预报，也是地震人的长期课题，任重而道远。本文利用高台地震台的观测资料，研究了张掖M_S5.0地震前的异常变化情况。

高台地震台，位于青藏块体东北缘，祁连山北缘段裂带中段，线壮褶皱的地质构造非常发育，垂直差异运动强烈，有一条NWW走向的逆冲大断层。地质构造上位于合黎山—龙首山褶皱带南缘与走廊断陷的分界处，河西系高台—榆木山隆起与祁吕西褶皱外缘与古北西向构造斜撞复合部位。台基为深层大面积裸露的完整花岗岩，台站南的地表表层为第四纪亚砂土、粉砂、沙砾石，厚度约5 m，下部全为花岗岩。

高台钻孔应变井于2007年3月建成，井深50 m，岩石均为海西期花岗闪长岩。仪器使用YRY-4分量式钻孔应变仪，安装在裸岩上，基岩无裂隙，无破碎，无涌水和漏水现象，钻孔用标准钢套管保护，然后用水泥浆压力封井，使得井套管与岩壁紧密凝固，底部垫入30 mm小石子，探头周边均匀灌入水泥石英沙进行固定偶合。四个方向元件方位角1路-65°，2路-20°，3路25°，4路70°（如图1），2007年6月20日安装调试工作完成，仪器开始观测运行，21日产出数据。自架设观测至今，运行稳定，完整率均保持在99.5%以上，能够记录到清晰的应变固体潮和同震变化。

图1 高台台分量式钻孔应变元器件分布图Fig.1 Distribution of component drilling strain components of Gaotai Station

北京时间2019年9月16日20时48分，在（38.60 N，100.35 E）发生了张掖M_S5.0级地震，震源深度11 km，震中距离高台地震台及钻孔应变115 km，中国地震台网给出的震源机制中心解（如图2）。震后经过资料分析和研究，发现高台钻孔应变观测存在明显的前兆异常，现进行回溯性分析，对资料进行积累和总结。

图2 张掖Ms5.0级地震震中分布及震源中心解图Fig.2 Epicenter distribution and focal center solution of Zhangye MS5.0 Earthquake

2.1 资料可靠性分析

高台台YRY-4钻孔应变仪，2019年全年数据连续可靠无断记，使用四元件钻孔应变观测的自洽方程^[7]，

S₁ + S₃ =S₂ + S₄ =2A(ε₁ + ε₂)（1）

对2019年的资料进行自洽方程检验，S₁+S₃与S₂+S₄基本一致重合较好（如图3 ），很好地验证了高台地震台的YRY-4仪器稳定数据可靠。

对高台钻孔应变观测数据进行调和分析，主要提取固体潮的M2波进行分析，将2015—2019年的资料数据进行计算，潮汐因子2015年开始至2019年相对稳定（表1），从另一方面证明了高台钻孔应变仪仪器稳定，数据可靠。

图3 高台台YRY-4自洽曲线图Fig.3 Self consistent curve of YRY-4 at Gaotai Station

表1 2015—2019年高台钻孔应变M2波潮汐因子表Table1 The tidal factor of M2 wave of borehole strain at Gaotai Station from 2015 to 2019

2.2 异常分析

2.2.1 年变趋势异常分析

将2013—2019年的整点值进行下载，排除标定台阶和人为干扰等记录到的影响事件，得到的数据曲线为钻孔应变的年变趋势图，对比2013—2018年的年变趋势发现，2013—2018年数据都在4月下旬上升，9月上旬下降。2019年2月20日发现NE分量数据加速下降，一直持续至7月22日达到最低点-64 745.46×10^-8，超出同期2013—2018年年变均值54.38%，NW分量数据持续下降至2019年7月8日，达到-65 512.91×10^-8，转平后至2019年8月3日，达到最低值-65 663.8×10^-8，最大超出同期2013—2018年年变均值17.68%，后转折回升至2019年9月16日发震时，NE分量和NW分量形成破年变低值异常现象（如图4）。

图4 2013—2019年高台钻孔应变年变对比图Fig.4 Comparison chart of borehole strain of Gaotai Station from 2012 to 2019

2.2.2 主应变与主方位角异常变化分析

钻孔应变的数据处理通常是把四个直接观测值换算成三个间接观测值如（式3），

将换算后的间接观测值按主应变—主方向（ε₁，ε₂，φ）进行分解，（式3）可转化为（式4）

再利用三角公式分别求得主应变和主方向如（式5），

将2019年高台钻孔应变数据进行换算后，将（ε₁，ε₂，φ）进行求解，每月抽取一组数据进行作图发现，2019年1月1日之后，半长轴ε₂持续缓慢增大，半短轴ε₁加速增大张力增强，4月2日ε₁达到最大值26 814.1×10^-8，后张力逐步减弱，9月16日达到最小值6441.6×10^-8，张掖M_S5.0级地震后，应力释放ε₁逐步恢复，12月02日逐步恢复至1月的主应变大小。方位角φ年变化相对平稳，但2019年3月23日快速下降，3月28日达到最小值-64.6327°后发生转折，逐步恢复至平稳状态。ε₂方向正好与榆木山断裂NWW一致，自仪器安装以来稳步增大，ε₁方向正好与断裂垂直，主方向的大小变化也很好的反映了张掖M_S5.0级地震的应力积累和释放过程，很好地表现了地震孕育的前兆异常^[8] （如图8）。

图5 高台钻孔应变替换量应力变化图与替换量数值图Fig.5 Stress variation diagram and numerical replace diagram in borehole strain at Gaotai Station

2.2.3 差分变换异常分析

在钻孔应变观测中需要排除气压干扰造成的影响，对高台钻孔应变换算后的S_a数值进行差分处理，得到分钟值变化速率曲线，能清楚地发现钻孔应变观测中的高频信号变化情况（如图6），发现2019年4月—2019年7月，高频毛刺明显增多。

2.2.4 小波变换异常分析

小波方法最早起源于1980年，就是将信号分解成不同频率的分量，按相应的频率进行数据分析的方法^[9-10]。在研究时间和频率的局域变化，一般情况下使用Daubechies^[11]和Mallat^[12]的离散小波变换。在脉冲形态信号的提取中，我们对Sa数据进行DB小波分解，从第一层2~4 min，到第九层512~1024 min，发现d4~d7在4月以后毛刺逐渐增多，9月以后逐渐减少，这些大小毛刺随时间不断变化，反映了在4月至9月之间，地震前岩石的断裂和裂隙逐渐增大（如图7）。

图6 高台钻孔应变Sa、S13、S24分钟值差分曲线图Fig.6 Difference curve of Sa，S13，S24 minute value of borehole strain at Gaotai Station

图7 高台钻孔应变Sa DB小波分解4-7阶分解图Fig.7 The 4-7 order decomposition diagram of Gaotai Station borehole strain based on Sa DB wavelet decomposition

2.2.5 超限率异常分析

经过高通滤波和小波分解得到的观测值都是在0线的附近上下波动，当毛刺很小的时候就很难观察出来。将滤波后的信号时间序列记为X（i i=1， 2，...，N）N为数据的总个数，求时间序列的均值X—，得到标准差S_D。把S_D作为一个信号检测阈值，超出这个范围记做超限，单位时间内超限点的个数叫做数量超限率，记为Ron。将高台钻孔应变2019年的替换数据S_a、S₁₃、S₂₄进行超限分析后，发现4月至8月数目超限率Ron明显增多，8月开始明显减少，直至发震（如图8）。

图8 高台钻孔应变Sa超限分析图Fig.8 Sa transfinite analysis chart of drilling strain in Gaotai

2019年9月16日张掖M_S5.0级地震发生在全国重点危险区内，通过对高台钻孔应变的资料进行分析，在正常的年变状态下，NE分量和NW分量出现了明显的异常现象。

（1）NE分量和NW分量出现了破年变低值异常， NE分量较往年最低点达到-64 745.46×10^-8，超出同日2013年-2018年年变均值54.38%，NW分量数据最低点达到最低值-65 663.8×10^-8，最大超出同日2013—2018年年变均值17.68%，年变转折点滞后出现破年变。

（2）求解（ε_1，ε_2，φ）后发现，2019年1月1日之后，半长轴ε₂持续缓慢增大，半短轴ε₁加速增大张力增强，4月2日ε₁达到最大值26 814.1×10^-8，后张力逐步减弱， 9月16日达到最小值6441.6× 10^-8，张掖M_S 5.0级地震后，应力释放ε₁逐步恢复，12月02日逐步恢复至1月的应力大小。方位角φ年变化相对平稳，但2019年3月23日快速下降， 3月28日达到最小值-64.6327°后发生转折，逐步恢复至平稳状态，ε₂方向正好与榆木山断裂NWW一致，自仪器安装以来稳步增大，ε₁方向正好与断裂垂直，主方向的大小变化也很好的反映了张掖M_S 5.0级地震的应力积累和释放过程。

（3）替换量S_a、S₁₃、S₂₄的分钟值差分曲线在3～8月出现了毛刺增多现象，大小毛刺表现了岩石的破裂情况，DB离散小波分解后发现在4～7阶，在4～9月毛刺振幅变大，数量增多，地震发生后毛刺消失。求得超限点S_D后，进行数量超限发现， 4～7月数目超限率Ron明显增多，8月开始明显减少至发震。

另外，此次地震前山丹地震台电阻率在2016—2018年出现了EW分量的破年变异常^[13]。根据此次地震前这些异常信息的相互印证，综合分析，在张掖M_S5.0级地震前高台钻孔应变的异常信息可靠，可能是此次地震的前兆异常。

致谢：感谢应急管理部国家自然灾害研究院邱泽华研究员提供的钻孔应变分析软件对本文的技术支撑，唐磊副研究员的无私指导。