Analysis of anomalous environmental responses to the 5.5 earthquake in Pingyuan County, Dezhou in august 2023 based on GNSS
-
摘要: 为了探测2023年8月6日山东德州平原县5.5级地震引起的异常环境响应，本文基于震中东北方向约26 km的GNSS观测数据，提出滑动时变频率方法并引入四分位距(interquartile range，IQR)法对地震引发的异常环境响应进行分析. 结果显示，地震发生后N和E方向坐标速度时间序列在9~10 s处出现显著异常，垂向在16 s时出现轻微异常. 滑动时变频率方法能够在速度时间序列异常跳变发生前5~10 s探测到明显变化，显示了更高的灵敏度. 研究还发现，在地震发生前一天，电离层在震中上空出现了异常扰动，异常值达到4 TECU. 研究表明：GNSS技术能有效探测地震过程中的异常环境响应，为地震监测与预警提供了新的视角和工具.
-
关键词:
- 全球导航卫星系统(GNSS) /
- 滑动时变频率 /
- 地震异常 /
- 四分位距(IQR)法
Abstract: To investigate the abnormal environmental response caused by the magnitude 5.5 earthquake in Pingyuan County, Shandong Province, on August 6th, 2023, this study based on the Global Navigation Satellite System (GNSS) observation data located about 26 km northeast of the epicenter, proposes a moving time-varying frequency method and incorporates the interquartile range (IQR) method to analyze the abnormal environmental responses triggered by the earthquake. The results indicate significant anomalies in the N and E directional coordinate velocity time series at 9 to 10 seconds after the earthquake, with a minor anomaly in the U direction at 16 seconds. Moreover, the sliding time-varying frequency method detects noticeable changes 5 to 10 seconds before the abnormal jump in the velocity time series, demonstrating higher sensitivity. The study also discovered abnormal disturbances in the ionosphere above the epicenter the day before the earthquake, with the anomaly reaching 4 TECU. This research demonstrates the effectiveness of GNSS technology in detecting abnormal environmental responses during earthquakes, offering a new perspective and tools for earthquake monitoring and early warning. -
0. 引 言
地震作为一种突发性、具有破坏力的自然灾害，其预警与监测一直是地球科学研究的重点和难点[1-2]. 传统的地震监测技术，如地震波分析、地质调查等，虽然在地震研究领域发挥了重要作用[3]，但在预警方面存在一定局限性. 近年来，GNSS技术在地震监测领域展现出了独特的优势[4-6]，特别是在捕捉地震前的异常环境变化方面[7]. 这项技术通过监测大气、电离层和地表的细微变化，为地震预测提供了新的思路和方法. 当今地震研究领域，GNSS的应用已经成为一项关键技术[8-10]. 张一彬等[11]研究2019年5月26日秘鲁北部Mw8.0地震发现，震前3天内出现了电离层赤道异常消失、电子含量显著减少(高达10 TECU)以及特定时段内沿纬度分布的双峰曲线特征消失和震中附近低谷现象. 翟笃林等[12]基于2008至2019年中国区域7个MS≥6.0地震的地基GPS 电离层电子总含量(total electronic content, TEC)观测数据分析，研究发现地震前2~6天电离层TEC出现负异常扰动，震时及震后2~7天内出现正异常扰动，异常扰动与震中距离和震级相关. GNSS技术不仅用于监测地震引起的电离层扰动，还能检测到海啸引起的电离层变化[13-14]，其高采样率观测提供了高频速度信息[15]. 文献[16]也综合考察了太阳活动、对流层效应以及地震和海啸活动对电离层的综合影响. 然而，如何从大量的GNSS观测数据中准确提取地震信号，以及如何利用这些数据进行有效的地震分析和预警，仍然是一个极具挑战性的问题[4,17].
本文以2023年德州平原县地震为背景，提出滑动时变频率方法对三维坐标速度时间序列进行分析，引入四分位距(interquartile range，IQR)法对三维坐标速度时间序列和震中上空电子含量时间序列进行分析，旨在深入探讨基于GNSS观测站数据探测地震引起的异常环境响应.
1. 数据和方法
1.1 研究区域和数据
2023年8月6日2时33分59秒，在山东省德州市平原县(37.16°N，116.34°E)发生5.5级地震，震源深度10 km，多地震感明显.
GNSS观测站为架设在混凝土房顶的临时GNSS观测站，完整记录了地震期间的数据，数据采样率为1 s，截止高度角为10°. 如图1所示，GNSS观测站的位置为37.26°N、116.61°E，在震中的东北方向，距震中的地面距离约25.6 km，距离震源约为28 km.
电离层数据来自欧洲定轨中心(Center for Orbit Determination in Europe, CODE)提供的全球电离层格网(global ionospheric map, GIM)资料(数据网址: ftp://ftp.unibe.ch/aiub)，其中包含格网点上的TEC. 数据时段为2023年7月6日至2023年8月13日. GIM数据的时间分辨率为2 h，空间分辨率为5°×2.5°.
为排除太阳活动对电离层的异常扰动，本文对中科院空间中心提供的2023年7月6日至2023年8月13日的太阳射电通量F10.7、太阳黑子数(sunspot number, SSN)及地磁指数Ap(数据网址: http://www.sepc.ac.cn)进行分析.
1.2 研究方法
1) IQR法
当使用IQR法[18]来处理数据时，首先需要计算数据集的上下四分之一位数(Q1和Q3). 这两个值可以帮助我们了解数据的分布情况，特别是在处理异常值和离群点时非常有效.
以下是计算IQR法的具体步骤：
步骤1：将数据集按照从小到大的顺序排列.
步骤2：计算第一四分位数(Q1). Q1是中位数与最小值之间的中位数，表示了前25%的数据的位置. 如果数据集的大小为奇数，Q1即为中位数左侧一半数据的中位数；如果数据集的大小为偶数，Q1即为中位数左侧一半数据的中位数.
步骤3：计算第三四分位数(Q3). Q3是中位数与最大值之间的中位数，表示了前75%的数据的位置. 同样，如果数据集的大小为奇数，Q3即为中位数右侧一半数据的中位数；如果数据集的大小为偶数，Q3即为中位数右侧一半数据的中位数.
步骤4：计算IQR. IQR是Q3与Q1的差值，它表示了数据集的中间50%范围内的变异程度. 数学公式为
$$ \mathrm{I}\mathrm{Q}\mathrm{R}=Q_3-Q_1 $$ (1) 步骤5：检测异常值. 在IQR法中，通常将小于
$Q_1-1.5*\mathrm{I}\mathrm{Q}\mathrm{R}$ 或大于$Q_3+1.5*\mathrm{I}\mathrm{Q}\mathrm{R}$ 的数据点视为异常值或离群点.2)频谱分析方法
本文使用的频谱分析方法为快速傅里叶变换(fast Fourier transform，FFT)，FFT能高效地将时域中的信号转换为频域，揭示其构成频率的复合结构，在信号处理、图像处理、通信等领域中有广泛的应用. 以下是FFT的具体步骤和公式描述.
步骤1：获取离散信号. FFT适用于离散信号，因此首先需要获得一个离散的时域信号，表示为
$ x\left[n\right] $ ，其中$ n $ 表示离散时间点.步骤2：计算复数指数项. FFT利用复数的旋转特性，计算复数指数项
$ {W}_{N}^{kn} $ ，其中N表示信号长度，k表示频率索引，n表示时间索引. 复数指数项的数学表达式为$$ {W}_{N}^{kn}={\text{e}}^{-\mathrm{j}\frac{2{\text{π} }}{N}kn} $$ (2) 步骤3：执行蝶形运算. FFT算法通过蝶形运算将时域信号分解为不同频率分量. 蝶形运算是一种迭代计算，它将信号分为两部分，然后进行乘法和加法运算. 这个步骤的关键是利用了复数指数项的周期性和对称性.
步骤4：重排结果. 在进行蝶形运算后，需要将计算得到的频域分量按照频率顺序重新排列，以获得正确的频谱表示.
步骤5：重复步骤3~4. FFT算法会递归地执行蝶形运算和频谱重排，直到得到所有频率分量的幅度和相位信息.
步骤6：计算频谱. 最终，通过FFT算法得到的频率分量的幅度和相位信息构成了信号的频谱表示. 频率分量的幅度表示了不同频率在信号中的强度，而相位表示了信号在不同频率上的相对偏移.
FFT的离散形式可以用以下公式表示，其中
$ X\left[k\right] $ 表示频域信号在频率索引$ k $ 处的分量值$$ X\left[k\right]=\sum _{n=0}^{N-1}x\left[n\right]\cdot {W}_{N}^{kn} $$ (3) 2. 结果与分析
为了使用GNSS技术探测2023年德州平原县5.5级地震导致的异常，我们参考震中与GNSS观测站的距离，截取了地震发生时前20 s和后80 s的时间序列速度数据，使用IQR法和滑动时变频率方法进行分析，探测地震发生前后的异常环境响应.
2.1 IQR法探测速度时间序列异常值
对原始坐标时间序列进行一阶导数求解，得到测站3个方向的运动速度，采用IQR法对速度进行分析，探测异常跳变. 结果图2~4所示.
从图2的数据分析来看，在地震发生后的第9 s和第10 s，观察到了N方向速度时间序列的显著异常. 这一时段内的速度超过了5.6 mm/s的阈值，约为正常速度的两倍. 与此同时，震后速度时间序列的振动频率也出现了显著变化. 在地震发生后的第10 s之后，速度时间序列逐渐回归到了稳定状态.
根据图3的结果，E方向的速度时间序列在地震后的第10 s显示出了显著的异常. 异常值超出了2.2 mm/s的阈值下限. 同时，震后速度时间序列的震动频率经历了显著的变化.
由图4可知，U方向速度时间序列在震后16 s时探测到异常，但异常值较小. 在地震发生后，U方向速度时间序列可以观察到明显的震动频率变化.
2.2 速度时间序列频率变化研究
在第2.1节中发现震后的速度时间序列变化频率出现了明显的变化，因此，本文提出滑动时变频率方法对三个方向的速度时间序列进行分析，每10个历元为1组，每1个历元为步长向后滑动，分析速度时间序列的频率随时间的变化，结果如图5~7所示.
由图5可得，从频率变化角度分析，N方向速度时间序列从第5~6 s开始出现明显变化，较速度时间序列的异常跳变提前5 s左右被探测到. 到第14~15 s左右逐渐恢复正常. 上述现象与图2中的结果基本吻合.
由图6可得，在地震发生时，E方向速度时间序列便可观察到明显的频率变化，比速度时间序列的异常跳变值提前出现10 s左右. 在第18 s左右逐渐恢复正常. 在第68 s左右又出现短暂的频率异常现象. 上述现象与图3中的结果基本吻合.
由图7可知，在地震发生后第6 s左右，U方向速度时间序列的变化频率出现明显变化，比速度时间序列的跳变值提前出现10 s左右. 而速度时间序列的异常跳变并不明显.
综上所述，从速度时间序列的变化频率上可以更加提前地观察地震引起的速度异常变化，且异常现象非常明显.
2.3 电离层异常探测
从CODE提供的GIM资料中提取出距离震中最近的GNSS TEC值，利用双线性内插得到震中处的TEC值. 利用滑动四分位距法，以30天同一时刻的观测值为背景值，1.5倍IQR为限差，对汤加火山爆发中心位置前30天的电离层异常进行探测. 探测结果如图8所示.
由图8(a)可知，在地震前的第26天左右、第16天左右以及第11天左右出现TEC的高峰值，在地震前第28天左右、第9天左右以及第1天左右出现了TEC低谷值.
由图8(b)可知，在地震前的第23天、第19天、第16天、第10天均出现明显的TEC正异常，异常值最高达4 TECU. 地震前的第28天、第9天、第4天及第1天，出现了明显的负TEC异常，异常最大可达5 TECU. 地震发生后，基本没有TEC异常现象出现.
然后，利用IQR法分析了地震前30天的太阳射电通量F10.7、地磁Ap指数和太阳黑子数SSN的异常信息，以排除太阳活动和地磁活动引起的TEC异常. 分析结果如图9所示.
由图9可知，Ap指数在地震前第26~28天和第4~8天有明显的正异常，在震前第18~22天有明显的负异常；F10.7在震前第24~30天有明显的负异常，在震前第6~19天出现明显的正异常；SSN的变化趋势与F10.7的变化趋势相似，在震前第11~21天出现明显的正异常现象，通过与TEC时间序列进行对比，发现震前0~4天太阳活动和地磁活动都未出现明显异常，而TEC出现了明显的负异常，因此，我们认为震前0~4天的TEC异常可能与此次地震有关.
根据震中处TEC异常扰动探测的结果，震前0~4天的TEC异常扰动尤为明显，因此我们对震前第1天的亚洲及周边区域TEC异常分布进行研究，结果如图10所示.
由图10可得，震前1天的0~8 h亚洲及周边区域的电离层没有出现特别明显的电离层异常现象，从震前1天的10 h开始，出现了较为明显的电离层异常现象，震前10~16 h，震中同纬度上出现了明显的电离层负异常现象，震前12 h在30°N~60°N出现了大范围的负异常现象，从震前12 h开始，电离层异常范围逐渐缩小.
3. 结束语
本文提出滑动时变频率方法对三维坐标速度时间序列进行分析，引入IQR法对三维坐标时间序列、三维坐标速度时间序列和震中上空电子含量时间序列进行分析，得出如下结论：
1)在地震发生后的第9 s、10 s左右，N和E方向坐标速度时间序列可见明显异常，最大异常超过正常值的50%，U方向速度时间序列在震后16 s时探测到轻微异常.
2)利用滑动时变频率方法对速度时间序列进行分析， N、E、U三个方向分别可以提前5 s、6 s、10 s左右探测到地震引起的速度时间序列异常.
3)从时间维度分析，在地震的前1天，地震引起震中上空的电离层出现明显的异常扰动，异常值为4 TECU. 从空间维度分析，从震前1天的10 h开始，出现了较为明显的电离层负异常现象，震前12 h电离层负异常的范围达到最大，从震前12 h开始，电离层异常范围逐渐缩小.
-
[1] 杨剑, 张宇, 胡良晨, 等. 2017年九寨沟M7.0地震前临震电离层异常研究[J]. 地震, 2022, 42(4): 100-110. [2] 蔡华, 李子申, 王敏, 等. 汶川Mw7.9与日本Mw9.0地震同震电离层扰动研究[J]. 武汉大学学报(信息科学版), 2013, 38(6): 716-719,724. [3] 张怀, 聂兆生, 刘刚, 等. 高频BDS/GPS融合PPP确定2021年玛多M_W7.4地震三要素[J]. 地震研究, 2023, 46(2): 216-225. [4] BESHR A, ZARZOURA F. Using artificial neural networks for GNSS observations analysis and displacement prediction of suspension highway bridge[J]. Innovative infrastructure solutions, 2021, 6(1): 1-15. DOI: 10.1007/s41062-021-00458-4
[5] TU R, ZHANG R, ZHANG P, et al. Integration of single-frequency GNSS and strong-motion observations for real-time earthquake monitoring[J]. Remote sensing, 2018, 10(6): 886. DOI: 10.3390/rs10060886
[6] WANG J, JIANG W, LI Z, et al. A new multi-scale sliding window LSTM framework (MSSW-LSTM): a case study for GNSS time-series prediction[J]. Remote sensing, 2021, 13(16): 3328. DOI: 10.3390/rs13163328
[7] 马玉, 祝芙英. 基于GPS TEC的尼泊尔M_W7.8地震同震电离层扰动研究[J]. 大地测量与地球动力学, 2020, 40(9): 957-961,969. [8] GAO Z, LI Y, SHAN X, et al. Earthquake magnitude estimation from high-rate GNSS data: a case study of the 2021 Mw 7.3 Maduo earthquake[J]. Remote sensing, 2021, 13(21): 4478. DOI: 10.3390/rs13214478
[9] 姚宜斌, 翟长治, 孔建, 等. 2015年尼泊尔地震的震前电离层异常探测[J]. 测绘学报, 2016, 45(4): 385-395. [10] 姜卫平, 许才军, 李志伟, 等. 利用空间观测技术研究青海玛多7.4级地震孕育发生变形时空特征[J]. 地球物理学报, 2022, 65(2): 495-508. [11] 张一彬, 钟慧鑫, 张亭, 等. 2019-05-26秘鲁北部地震前低纬电离层变化分析[J]. 大地测量与地球动力学, 2023, 43(2): 135-140. [12] 翟笃林, 祝芙英, 林剑, 等. 基于地基GPS-TEC的中国区域地震电离层效应研究[J]. 中国地震, 2020, 36(4): 857-871. [13] 唐龙, 郭博峰, 李哲. 利用日本GPS网探测2011年Tohoku海啸引发的电离层扰动[J]. 地球物理学报, 2017, 60(2): 507-513. [14] LARSON K, LAY T, YAMAZAKI Y, et al. Dynamic sea level variation from GNSS: 2020 Shumagin earthquake tsunami resonance and hurricane Laura[J]. Geophysical research letters, 2021, 48(4): e2020GL091378. DOI: 10.1029/2020GL091378
[15] XU P L, DU F, SHU Y M, et al. Regularized reconstruction of peak ground velocity and acceleration from very high-rate GNSS precise point positioning with applications to the 2013 Lushan Mw6.6 earthquake[J]. Journal of geodesy, 2021, 95(1): 1-22. DOI: 10.1007/s00190-020-01449-6
[16] SUN Y. GNSS brings us back on the ground from ionosphere[J]. Geoscience letters, 2019, 6(1): 1-9. DOI: 10.1186/s40562-019-0144-0
[17] KING L S, UNWIN M, RAWLINSON J, et al. Towards a topographically-accurate reflection point prediction algorithm for operational spaceborne GNSS reflectometry-development and verification[J]. Remote sensing, 2021, 13(5): 1031. DOI: 10.3390/rs13051031
[18] EH N. Elementary statistical methods[J]. Nature, 1934, 133(3349): 9. DOI: 10.1038/133009c0
-
期刊类型引用(1)
1. 董龙飞，李艳艳，龙海立. 联合CEEMD和改进修正能量比的高频GNSS同震位移P波震相识别方法. 全球定位系统. 2025(01): 25-32+59 . 
本站查看
其他类型引用(0)
下载:
计量
- 文章访问数: 210
- HTML全文浏览量: 83
- PDF下载量: 25
- 被引次数: 1
