摘要 2017 年 8 月 8 日中国四川省九寨沟县 Ms 7.0 级地震前, 在震源 200 km 范围内的九寨沟、松潘、平武、青川、茂县和汶川布设 AETA 系统的 6 个监测台。采用窗口为 30 日的滑动主成分分析(principal component analysis, PCA)方法, 对 6 个台站的电磁扰动进行特征分析, 探讨九寨沟 Ms7.0 级地震前后 AETA 系统观测的电磁场变化特征, 并由此提出一种新的临震特征, 即 AETA 电磁数据经滑动 PCA 方法处理后得到的“条带异常”。结果显示: 1)除九寨沟 Ms 7.0 级地震发震当日, 九寨沟监测台于震前 5 日至震后 15 日, 每天在同一时段(5:00am—6:00am)出现异常点, 形成“异常条带”, 异常程度在震前大体上呈递增变化, 震后大体上呈递减变化; 2)震源附近其他监测台也存在异常点, 但发布零散且异常程度较小; 3)九寨沟 Ms7.0 级地震“异常条带”消失后的两个月内, 未出现新的“异常条带”, 且无强震发生。
关键词 九寨沟地震; AETA; 电磁场; 滑动PCA方法
2017 年 8 月 8 日 21 时 19 分 46 秒, 四川省九寨沟县发生 Ms7.0 级地震, 截至 8 月 13 日 20 时, 导致25 人死亡, 525 人受伤[1]。
地震孕育及发生过程中, 震源及周围地区往往会提前数天或数小时出现不同程度的电磁辐射现象[2-4]。观测资料表明, 电磁异常是对短临地震反应最敏感的地震前兆之一[5-6]。因此, 该现象引起国内外学者的广泛关注, 许多学者针对某一特定频段的信号, 采用多种方法(如主成分分析法、滑动四分位法[13]和极化法[14-15]等)分析震前的标识异常[7-12], 期望在短临地震的预报方面有所突破。滑动四分位法的探测效果易受太阳地磁活动的影响[7], 极化法需要使用空间的 3 个分量, 所以这两种方法都具有一定的局限性。主成分分析法不受此局限, 作为一种常用的统计方法, 能够分离不同的信号, 有效地识别相对较弱的信号, 因此被引入地震电磁扰动分析中[13]。Chang等[16]将 PCA、滑动时窗法以及张小红等[17]提出的限差确定策略相结合, 提出滑动 PCA 方法, 并将其运用到 4 次强震(青海地震、本州地震、于田地震和尼泊尔地震)前后的全球电离层格网数据(GIM)分析中, 发现震前 10 日内皆出现短时的震前异常, 由此证明了该方法可用于地震前的异常监测。
我国目前的数字地震台网大多密度低, 台距较大[18], 布设在全国各地的电磁观测仪器也各不相同[19]。基于此背景, 本文使用的数据来源于我们设计的多分量地震监测系统 AETA[20]。AETA 可以同时监测电磁和地声两种信号, 已经在川滇等地区实现较大区域的密集布设, 且具备很好的一致性[21], 可以更好地监测地震前后的电磁信号。
本文将滑动 PCA 方法运用到 AETA 系统的低频电磁信号中, 对距九寨沟 Ms7.0 级地震震源 200km范围内的 6 个监测台数据进行分析, 探测九寨沟Ms7.0 级地震前后地表电磁场的变化规律。结果发现距, 震中 40km 的九寨沟监测台出现明显的异常条带, 我们认为该异常是一种较为明确的临震特征。
多分量地震监测系统(AETA)由数据处理终端、地声传感探头、电磁传感探头、监测数据云平台和数据分析系统组成, 可以感知来自地下的电磁扰动和地声信号。实时采集的数据通过互联网(有线或无线)传输到云平台, 进行存储、特征提取和异常分析。截至目前, 在全国范围内安装 AETA 系统约200 余台, 遍布河北、四川、云南、西藏、广东和台湾地区, 其中在四川布设设备数量达 93 台, 基本上覆盖四川全境重点区域。
九寨沟 Ms 7.0 级地震发生时, 距震源 200km 内共有 6 个 AETA 系统监测台, 分别位于九寨沟县、松潘县、平武县、青川县、茂县和汶川县(表 1)。图 1 展示本文使用的部分低频电磁信号的特征数据(EMD)。
表1 各监测台安装时间及地理位置
Table 1 Installation time and location of monitoring stations
监测台缩写安装日期地理位置震中距/km 九寨沟JZG2017-06-1033.25°N, 104.24°E40 松潘SP2017-06-1232.65°N, 103.60°E64 平武PW2017-06-0832.41°N, 104.55°E111 青川QC2017-06-0732.59°N, 105.23°E147 茂县MX2017-06-1331.69°N, 103.85°E167 汶川WC2017-06-1331.48°N, 103.59°E192
说明: 九寨沟Ms 7.0级地震震中位于33.20°N, 103.82°E。
本文运用 Chang 等[16]提出的滑动 PCA 方法, 并结合 AETA 系统低频电磁信号的特征数据, 进行地震预测分析。
设时间序列矩阵
(1)表示时长 N 日、每日数据时间间隔为 1 小时的原始数据。通过计算, 得到原始矩阵 X 的协方差矩阵 C, 求得协方差矩阵的特征值 λ1≥λ2≥…≥λn 及对应的特征向量矩阵E, 则
(2)
其中
, 
。选取前k个特征值, 使累计贡献率
≥85%,
。 (3)设置
,
, 计算重构矩阵
, (4)
选取重构矩阵 XP 的最后一列 Xref 为当日背景参考值。
为防止 AETA 系统在运行过程中可能受偶然环境干扰的影响, 并提升模型的鲁棒性, 我们设立合理的误差范围。
图1 本文所用部分特征数据
Fig. 1 Part of the characteristic data used in this paper
图2 2018年6月27日至2018年7月3日不同时间窗口下的评估指标(UTC +8)
Fig. 2 
in different windows from June 27, 2018 to July 3, 2018 (UTC +8)
首先计算当日实际值
与参考值
差值的绝对值, 按从小到大排列, 选取第 23 位(95.8%)[17], 得到公差
, 并依此进行异常判定。
计算异常限值:
(5)计算异常程度:
(6)
若
, 则认为存在正异常; 若
, 则认为存在负异常。
首先利用 PCA 方法构建主成分模型, 得到目标日期的参考背景值, 探测该日的异常程度, 然后将样本向后滑动 1 日, 继续计算, 直至计算完毕, 最后利用得到的异常值绘制热力图, 观察异常情况。
在实际情况中, 可能由于监测台断电、断网等原因, 出现数据短时空缺现象。如果样本数据存在不超过 7 日的短时空缺[17], 则使用相邻日期对应时间的数据进行填补, 保证数据结构及长度一致即可, 不影响预测结果。
为了确定滑动 PCA 方法在 AETA 电磁数据背景下的滑动窗口长度, 本文选取九寨沟监测台为期7 天(2018 年 6 月 27 日至 2018 年 7 月 3 日)的低频电磁信号特征数据, 滑动窗口长度分别为 20, 25, 27, 30, 32, 35 和 40日, 评估指标为
(7)距该时段前后两个月, 100km 范围内无 2.5 级以上地震发生(中国地震台网, http://news.ceic.ac.cn/)。评估结果如图 2 所示, 可以看出当滑动步长发生变化时, 评估指标整体变化趋势趋于一致。此外, 由表 2 的统计数据可知, 当滑动窗口为 30 日时, 平均评估指标最大, 因此选定 30 日作为本模型的滑动窗口长度。
图 3 表示各监测台 2017 年 8 月 1 日至 2017 年 8月 31 日电磁信号特征数据异常值的热力图。可以看出, 地震发生前 5 日, 九寨沟监测台 5—6 时出现逐渐增强的连续异常条带。地震发生前 2 日, 异常值大于 2, 明显高于历史背景中其他异常点; 地震发生时, 异常点位于 21—22 时, 表明此时的异常程度大于该日其他时间; 地震发生后, 5—6 时再次出现大体上逐渐减弱的连续异常条带, 持续 15 日; 震源附近其他监测台的异常程度较小, 且分布比较散乱, 没有构成条带状异常。
表2 2018年6月27日至2018年7月3日不同时间窗口下的平均评估指标(UTC+8)
Table 2 Average Prelin different windows from June 27, 2018 to July 3, 2018 (UTC +8) %
月/日20日25日27日30日32日35日40日 6/2781.8180.9179.3778.9679.2776.5175.27 6/2884.6885.8685.4286.7686.1286.2585.92 6/2986.4088.1688.1590.2489.5289.3991.60 6/3084.1381.3381.2479.5379.0979.2876.55 7/189.3289.4589.7990.3389.1989.6488.19 7/291.5993.1593.0093.2691.9492.3390.08 7/382.6386.6088.8688.6689.7089.8890.64 平均值85.7986.4986.5586.8286.4086.1885.47
图 4 为九寨沟监测台 2017 年 7 月 14 日至 11 月 12日的电磁信号异常值热力图。可以看出, 九寨沟监测台 8 月的条带状异常结束后, 未形成新的高幅值连续条带。
为了探究九寨沟 Ms 7.0 级地震前后“异常条带”这一临震特征的偶然性, 我们将滑动 PCA 方法运用到 AETA 系统其他台站的电磁数据中, 也发现类似的异常(如冕宁监测台)。表 3 为冕宁监测台 2017 年5 月 4 日至 2018 年 1 月 19 日, 距离冕宁监测台 150km 范围内的地震数据(来自中国地震台网 http:// news.ceic.ac.cn), 图 5 表示冕宁监测台的原始数据经滑动 PCA 方法处理后得到的异常值热力图, 可以看到在距离冕宁监测台 90~110km 的范围内, 条带状异常发生 3~4 天后, 皆会发生地震; 没有条带状异常时, 没有地震发生。
AETA 系统自 2016 年开始运行, 当电磁、地声信号出现异常时, 工作人员会向当地运维人员调查核实, 无一例是人为干扰所致。类似的震前异常现象在其他台站同样存在, 并取得切实的映震效果。所以数据是可靠的。我们也考察了空间活动的影响, 2017 年 8 月 1 日至 2017 年 8 月 14 日 九寨沟地区未发生磁暴现象(空间环境预报中心, http://www. sepc.ac.cn), 因此九寨沟台站监测到的异常也不太可能是空间活动造成的。在排除上述影响因素后, 我们认为该异常现象可能是附近地下电导率发生改变造所致, 并且很可能与数日后发生的地震活动有一定的关系。
图3 九寨沟Ms 7.0级地震震源附近监测台电磁信号的异常值热力图(UTC+8)
Fig. 3 Thermograms of abnormal value of electromagnetic signals of monitoring stations near Jiuzhaigou Ms 7.0 earthquake source (UTC+8)
图4 九寨沟监测台2017年7月14日至11月12日电磁信号的异常值热力图(UTC+8)
Fig. 4 Thermogram of abnormal value of electromagnetic signals of Jiuzhaigou monitoring station from July 14 2017 to November 12, 2017 (UTC+8)
表3 2017 年 5 月 4 日至 2018 年 1 月 19 日冕宁监测台 150 km 范围内地震情况(UTC +8)
Table 3 Earthquakes within 150 km of Mianning monitoring point from May 4, 2017 to January 19, 2018 (UTC+8)
发震时刻(UTC+8)震级(Ms)震中位置深度/km震中距/km 2017-08-30 21:46:173.027.90°N, 101.37°E8106.54 2017-09-12 18:40:103.227.92°N, 101.42°E13104.35 2017-09-12 19:25:593.027.91°N, 101.39°E15104.35 2017-09-12 19:26:404.427.93°N, 101.42°E13100.68 2017-09-21 10:03:123.027.93°N, 101.41°E15101.40 2017-10-18 02:54:203.328.33°N, 102.83°E1868.97 2017-10-25 22:29:373.027.90°N, 101.43°E13102.30
说明: 冕宁监测台位于28.55°N, 102.17°E。
图5 2017年5月4日至2018年1月19日冕宁监测台电磁信号的异常值热力图(UTC+8)
Fig. 5 Thermogram of abnormal value of electrimagnetic signals of Mianning monitoring station from May 4, 2017 to January 19, 2018 (UTC+8)
本文使用滑动 PCA 模型, 对九寨沟 Ms 7.0 级地震震源 200km 范围内的 AETA 系统 6 个监测台的低频电磁信号进行分析, 监测到地震前后电磁信号的明显异常。
1)在九寨沟 Ms7.0 级地震震前 5 日及震后 15日, 九寨沟监测台的地磁信号在同一时段出现异常点, 即条带状异常, 异常程度在震前大体上呈现递增变化, 震后大体上呈现递减变化。
2)震源附近其他监测台虽也存在异常点, 但分布零散, 异常程度较小。
3)在九寨沟 Ms 7.0 级地震异常条带消失后的两个月内, 九寨沟监测台未出现明显的条带异常, 且无强震发生。
我们认为, 如果监测台出现类似的条带异常, 则意味着具有高的地震风险, 震中可能位于此监测台附近。该异常表征具备很强的直观性。随着设备布设与震例积累, 在后续工作中, 我们将对此类异常与地震的时间、地点和震级的定量关系进行深入的研究。
致谢 AETA 实验得到中国地震局以及四川省、云南省、河北省地震局的大力支持, 尤其在九寨沟 Ms7.0 级地震后, 得到阿坝藏族羌族自治州和九寨沟县防震减灾局的支持和配合。AETA实验系统的研发得到深圳市深创谷技术服务有限公司的设计与制造服务支持。一并表示感谢。
参考文献
[1]Lei H, Wang X, Hou H, et al. The earthquake in Jiuzhaigou County of Northern Sichuan, China on August 8, 2017. Natural Hazards Journal of the Inter-national Society for the Prevention & Mitigation of Natural Hazards, 2018, 90(2): 1021-1030
[2]Molchanov O A, Mazhaeva O A, Protopopov M L, et al. Electromagnetic VLF radiation of seismic origin observed on the interkosmos-24 satellite. Geomagne-tizm I Aeronomiya, 1992, 32(6): 128-137
[3]Hirano T, Hattori K. ULF geomagnetic changes pos-sibly associated with the 2008 Iwate-Miyagi Nairiku earthquake. Journal of Asian Earth Sciences, 2011, 41(4/5): 442-449
[4]Karakelian D, Klemperer S L, Fraser-Smith A C, et al. Ultra-low frequency electromagnetic measurements associated with the 1998 Mw 5.1 San Juan Bautista, California earthquake and implications for mech-anisms of electromagnetic earthquake precursors. Tec-tonophysics, 2002, 359(1/2): 65-79
[5]Bleier T, Freund F. Earthquake. IEEE Spectrum, 2005, 42(12): 22-27
[6]Huang Q H. Retrospective investigation of geophysi-cal data possibly associated with the Ms 8.0 Wen-chuan earthquake in Sichuan, China. Journal of Asian Earth Sciences, 2011, 41(4/5): 421-427
[7]Zhang X, Zhang T, Pei J, et al. Time series modelling of syphilis incidence in China from 2005 to 2012. PLOS ONE, 2016, 11(2): e0149401
[8]Sivavaraprasad G, Ratnam D V. Performance evalua-tion of ionospheric time delay forecasting models using GPS observations at a low-latitude station. Advances in Space Research, 2017, 60(2): 475-490
[9]Chavez O, Pérez-Enríquez R, Cruz-Abeyro J A, et al. Detection of electromagnetic anomalies of three earthquakes in Mexico with an improved statistical method. Nat Hazards Earth Syst Sci, 2011, 11: 2021-2027
[10]Verma M, Steffen M, Denker C. Evaluating local correlation tracking using CO5BOLD simulations of solar granulation. Astronomy & Astrophysics, 2013, 555: A136
[11]Hattori K, Han P, Yoshino C, et al. Investigation of ULF seismo-magnetic phenomena in Kanto, Japan during 2000-2010: case studies and statistical stu-dies. Surveys in Geophysics, 2013, 34(3): 293-316
[12]Hayakawa M, Hobara Y, Ohta K, et al. The ultra- low-frequency magnetic disturbances associated with earthquakes. Earthquake Science, 2011, 24: 523-534
[13]邹斌, 郭金运, 常晓涛, 等. 基于主成分分析与滑动四分位法的震前 TEC 异常探测对比分析. 全球定位系统, 2016, 41(4): 63-69
[14]Li Q, Zhu P Y, Mamatemin A, et al. Detection of ULF electromagnetic emissions as a precursor to two earthquakes in China. Earthquake Science, 2011, 24 (6): 601-607
[15]Li Q, Yang X, Cai S P. Case study of applying polarization geomagnetic array date. Technology for Earthquake Disaster Prevention, 2015, 10(2): 412-417
[16]Chang X, Zou B, Guo J, et al. One sliding PCA method to detect ionospheric anomalies before strong earthquakes: cases study of Qinghai, Honshu, Hotan and Nepal earthquakes. Advances in Space Research, 2017, 59(8): 2058-2070
[17]张小红, 任晓东, 吴风波, 等. 震前电离层 TEC 异常探测新方法. 地球物理学报, 2013, 56(2): 441-449
[18]陈运泰. 地震预测——进展、困难与前景. 地震地磁观测与研究, 2007, 28(2): 1-24
[19]蔡润, 武震, 谭大诚, 等. 地震前的电磁异常综述. 华南地震, 2018, 38(1): 1-16
[20]王新安, 雍珊珊, 徐伯星, 等. 多分量地震监测系统 AETA 的研究与实现. 北京大学学报(自然科学版), 2018, 54(3): 487-494
[21]雍珊珊, 王新安, 庞瑞涛, 等. 多分量地震监测系统AETA 的感应式磁传感器磁棒研制. 北京大学学报(自然科学版), 2018, 54(3): 495-501
Research on Jiuzhaigou Ms 7.0 Earthquake Based on AETA Electromagnetic Disturbance
Abstract On August 8, 2017, an Ms7.0 earthquake hit Jiuzhaigou County, Sichuan Province, China. Before the earthquake, the AETA system was equipped with 6 monitoring stations in Jiuzhaigou County, Songpan County, Pingwu County, Qingchuan County, Maoxian County and Wenchuan County of Sichuan Province within a 200 km radius of the epicenter. With sliding principal component analysis (PCA) method, whose length of time window is 30 days, we analyzed the characteristics of electromagnetic disturbance of these six stations, and discussed the co- and post-seismic changes of the electromagnetic field observed by the AETA system, and proposed a new impending seismic feature, namely the “anomaly band”. The results show that except the day of Jiuzhaigou Ms7.0 earthquake, anomalous points appeared on monitoring station of Jiuzhaigou in the same period (5:00 am-6: 00 am) from the 5th day before the earthquake to the 15th day after the earthquake, and the anomalous value generally presented the trend of increase before the earthquake and decrease after the earthquake. Although anomalous points appeared on other monitoring stations, these points were scattered and the anomaly value was low. Within Two months after the abnormal band disappeared, there were no obvious band anomaly and no strong earthquakes.
Key words Jiuzhaigou earthquake; AETA; electromagnetic field; sliding principal component analysis method
doi: 10.13209/j.0479-8023.2019.108
深圳市科技创新委员会项目(JCYJ20170412151159461, KJYY20170721151955849)资助
收稿日期: 2018-11-29;
修回日期: 2019-01-17