2. 河南测绘职业学院, 河南 郑州 451464
2. Henan College of Surveying and Mapping, Zhengzhou 451464, China
地表瞬态形变主要源于火山喷发、地下水变迁及地下断层的无震蠕滑[1-3]。无震蠕滑通常表现为几天至几个月的缓慢滑移，这种慢滑移事件常伴随断层破裂，同时释放出巨大能量[4-5]。该能量不断积聚于周围岩层相对闭锁的区域，为地震的发生提供应力积累。因此无震蠕滑可能成为地震发生的前兆，准确提取蠕滑信息有助于区域地震危险性评估及震后地壳形变机制的研究。
国外科研人员已经开展了瞬态蠕滑探测的相关研究[6-10]。文献[11]将GNSS时间序列描述为由地表瞬态形变、参考框架误差、共性误差和噪声等信息构成的模型，并结合卡尔曼滤波估计瞬态蠕滑信息。文献[12]利用卡尔曼滤波结合主成分分析(principal component analysis，PCA)进行时空滤波，提高GNSS观测数据的信噪比，建立基于一阶马尔可夫理论的无震瞬态蠕滑信息检测模型。文献[13]考虑到共模误差的多样性，利用多频段奇异谱分析探测瞬态形变信息。文献[14]结合GNSS数据和气象资料使用双曲正切函数计算地表瞬态形变。国内学者也做了大量关于无震蠕滑信号检测的研究[15-18]。目前，研究人员多是将整个区域作为整体，滤波后利用所有测站间信息的相关性来探测信号，如果蠕滑信息的相关性较差，就会影响信号的探测效果。
地壳形变信号常淹没于GNSS的背景噪声中，采用合理的方法进行时空滤波有利于瞬态形变信息的提取。PCA是目前较为常用的时空滤波方法[19-20]，但研究表明该方法存在过度滤波现象[21-22]。而独立成分分析(independent component analysis，ICA)克服了该缺点并且能够有效估计高阶项信息。
本文基于GNSS坐标时间序列波动的相对强度指数提出一种蠕滑形变信息的自动检测方法。首先采用ICA进行时空滤波，提高坐标序列的信噪比；而后计算单一测站时间序列波动的相对强度指数及峰度值，通过累积分布函数将其转换为蠕滑信号概率，进而探测断层蠕滑事件。该方法单独计算各测站异常信号的发生概率，克服了测站分布及信号相关性的影响。
1 方法 1.1 独立成分分析独立成分分析是一种盲信号分离(blind source separation，BSS)方法，能有效提取坐标时间序列中的高阶项信息，恢复原始信号[23-24]。假设观测数据由n个GNSS连续观测站的m天单日解组成，l(l≤min(n, m))为分离出的源信号数量。ICA数学模型如式(1)所示
(1)
式中，X、A、s分别表示观测矩阵、混叠矩阵和源信号矩阵。混叠矩阵A的列表示信号的空间响应，源信号矩阵s的行表示信号的时间响应。
滤波后坐标时间序列如式(2)所示
(2)
式中，Ak、Sk分别表示第k个主分量对应的空间响应和时间响应；dn×m表示滤波后的时间序列。
1.2 相对强度指数文献[25]首次将相对强度指数引入股票市场，监测股值的异常震荡。本文利用滤波后GNSS坐标时间序列的相对强度指数判别异常波动，进而探测地壳瞬态蠕滑信息。分别计算各测站单日解时间序列的上升U(ti)(di>di－1)和下降趋势D(ti)(di < di－1)。
趋势上升时
(3)
(4)
趋势下降时
(5)
(6)
式中，di表示滤波后所求测站第i天的坐标值。式(2)中矩阵dn×m表示所有测站滤波后的坐标序列矩阵，该矩阵的行和列分别表示测站及观测值时间序列，因此di即为dn×m中所求测站对应行的第i列数值。U(ti)、D(ti)分别表示该测站第i天的上升和下降偏移值，如为上升趋势，则上升偏移值为U(ti)=di－di－1，下降偏移值为D(ti)=0；反之，下降偏移值为D(ti)=di－1－di，上升偏移值U(ti)=0。每个测站均可求得一个上升偏移值时间序列和下降偏移值时间序列。式(3)—(6)显示如果地表没有明显位移，U(ti)、D(ti)均接近于0，当出现蠕滑现象时，U(ti)或D(ti)异常，明显大于0。
式(7)采用移动平均值计算相对强度(relative strength, RS)。n表示移动周期，周期过长将导致相对强度时间序列过于平滑，无法准确跟踪较短周期的瞬态蠕滑。如果移动周期较短，则对测站位移过于敏感，易受非构造等外界因素影响，造成异常值过多、难以分辨蠕滑信号。通常情况下，蠕滑过程持续数周至几个月，而长期的慢滑移有时可持续数年之久，本文主要探测持续时间较短的瞬态蠕滑信号，对于长期的慢滑移事件还需要其他专业方法进一步研究[26-27]。参照文献[28]的设置，本文以Akutan和四川省为例将移动周期n设置为3周(21 d)。利用式(8)进一步获取相对强度指数(relative strength index, RSI)后，参考文献[25, 28]的设置，将大于60或小于40的RSI值称为异常波动，可能为断层的瞬态蠕滑信号
(7)
(8)
峰度通常用来描述一组数据分布的极端性，峰度值越大，该组数据的极端性越强。累积分布函数(cumulative distribution function, CDF)能够完整描述一个实数变量的概率分布。本文利用累积分布函数将峰度值转化为概率分布，即得到各历元处异常信号出现的概率。首先将RSI时间序列以50为界分成两组，依次剔除每组远离中心的极值，获取峰度值时间序列，如式(9)所示；而后利用累积分布函数将峰度值时间序列转化为瞬态蠕滑事件发生的概率，如式(10)所示。计算过程中为保证峰度值序列的稳定性，要求RSI时间序列的历元数不少于200[28]
(9)
(10)
式中，r(下标r含义相同)、i分别表示RSI和第i个历元；K、μ、σ分别表示峰度值、数学期望和标准差。
2 仿真试验本文设计了两个仿真试验。第1个试验验证PCA、ICA是否存在过度滤波现象。利用Fakenet软件模拟生成与陆态网测站水平相当的共模误差、白噪声、闪烁噪声和随机游走噪声合成6个测站连续3 a坐标时间序列，如图 1所示。分别利用PCA和ICA提取共模误差，如图 2所示，计算两种方法滤波前后时间序列的标准差，见表 1。表中显示两种方法滤波后各测站标准差的平均值明显降低，表明两种方法均具有较好的滤波效果。但利用PCA滤波后的标准差低于滤波前未加入共模误差的噪声序列的标准差，表明利用PCA滤波时可能滤除了部分本应留存的信号，因此产生了过度滤波现象。这与文献[21]的研究成果一致。
|
| 图 1 模拟测站时间序列 Fig. 1 Time series of simulative stations |
|
| 图 2 ICA、PCA提取的CME Fig. 2 CME extracted by ICA and PCA |
| 滤波前 | 滤波后 | ||||
| 噪声+CME | 噪声 | PCA | ICA | ||
| 标准差 | 2.602 5 | 1.722 6 | 1.571 3 | 2.158 7 | |
第2个试验验证本文所述蠕滑信息探测方法的可行性。同样利用Fakenet软件模拟5个测站连续500 d的时间序列。该序列由与陆态网测站水平相当的共模误差、噪声(白噪声+闪烁噪声+随机游走噪声)和线性趋势项组成。同时在第1、3、4、5这4个测站中加入了从第201 d起连续25 d的蠕滑信号，合成地表位移时间序列。蠕滑信号与合成的地表位移时间序列如图 3、图 4所示。
|
| 图 3 蠕滑信号 Fig. 3 Creep signal |
|
| 图 4 地表位移时间序列 Fig. 4 Displacement time series |
图 5为消除线性趋势利用ICA滤波后的时间序列。滤波后各测站时间序列的振幅减小，信噪比明显提高。利用式(8)、式(10)分别计算该模拟区域GNSS网相对强度指数和蠕滑形变的发生概率，如图 6、图 7所示。图 6显示多数历元的相对强度指数在50附近波动，少数历元波动幅度略大，但只有第201 d突破了前文所给的阈值(60/40)。同样，图 7中多数历元信号异常的概率为0，少数历元检测到微小的异常信号，但概率均未超过0.2。第201 d的概率突然增大，接近0.4，且明显高于其他历元，表明这一天信号异常。这与本文加入蠕滑信号的起始时间吻合。
|
| 图 5 滤波后时间序列 Fig. 5 Time series after filtering |
|
| 图 6 相对强度指数 Fig. 6 Relative strength index |
|
| 图 7 蠕滑概率 Fig. 7 Creep probability |
图 8所示为5个仿真测站第201 d的CDF值，图中显示第1、3、4、5这4个测站均探测出异常信号，且第1、5测站异常值较大。第2测站没有探测出异常信号，这与仿真时没有加入蠕滑信号的先验条件一致。第3、4测站虽然加入了蠕滑信号，但探测出的异常概率较小，可能与该测站模拟的背景噪声有关。通过大量试验模拟不同信噪比的地表位移时间序列，结果表明当蠕滑信号至少与噪声水平相当时，可通过相对强度指数来探测无震蠕滑信息。
|
| 图 8 仿真测站CDF值 Fig. 8 CDF values of simulation stations |
3 实例分析 3.1 Akutan瞬态无震蠕滑信息探测
Akutan位于美国阿拉斯加—阿留申火山带，该区域受太平洋板块向欧亚板块俯冲挤压影响，构造运动强烈，火山和地震灾害频发。
测站分布与滤波后N方向坐标时间序列如图 9、图 10所示。利用分布于Akutan地区PBO网络7个GNSS基准站。数据计算异常信号发生概率的时间序列(图 11)，图 11显示2008年2月13日和2009年1月5日概率值突出，均超过0.4，表明这两个时间可能发生了地表异常形变。研究发现2008年初Akutan地区产生了超过100 d的无震蠕滑[29-30]。由于Akutan火山十分活跃，此次蠕滑事件可能与火山岩的强烈运动有关。2009年1月5日异常概率虽然较高，但是否为蠕滑事件尚待进一步研究。
|
| 图 9 Akutan地区GNSS测站分布 Fig. 9 Distribution of GNSS stations in Akutan |
|
| 图 10 滤波后N方向时间序列 Fig. 10 N time series after filtering |
|
| 图 11 异常信号概率 Fig. 11 Probability of abnormal signals |
3.2 四川省地表位移异常信号探测
四川省位于我国川滇地震带，地质条件复杂，地震活动频繁，是国内外学者开展地震研究的重要地区。本文选取四川省境内陆态网18个基准站2011—2017年的观测数据开展区域性异常信号的探测研究。测站分布及主要地震的震源机制如图 12所示。
|
| 图 12 测站分布及主要地震震源机制 Fig. 12 Stations distribution and main earthquake focal mechanism |
使用GAMIT/GLOBK软件解算并获取各测站ITRF2008坐标。利用QOCA软件扣除线性趋势，年、半年周期项，仪器和天线变更等引起的阶跃项等因素的影响，滤波后N方向时间序列如图 13所示。最后计算该地区7年来异常信号出现的概率，如图 14所示。
|
| 图 13 滤波后N方向时间序列 Fig. 13 N time series after filtering |
|
| 图 14 异常信号概率 Fig. 14 Probability of abnormal signals |
由于四川省地质活动频繁，因此图 14中各历元信号异常的概率较为“杂乱”。2011年3月8日和2017年12月14日信号突出。哈佛CMT目录显示这两个信号出现前后2—21 d云南盈江和重庆武隆分别发生Mw5.5和Mw4.9地震。由此推断这两个信号可能为震前应力积累造成的地表异常移动和地震发生后引起的震后余滑。异常信号概率较高表明这些地震震级大、破坏力强、地表位移明显。
值得注意的是，如果将该区域概率阈值降低至0.14，能够探测出另外2个较为明显的异常信号。两个信号均出现在震级较高(云南腾冲Mw5.1地震、四川九寨沟Mw6.5地震)的地震发生后5—30 d之内，推断这些地震发生后可能产生了震后余滑，引起了相应地区的地表位移。上述现象是否与地震有关，有待进一步研究。信号出现时间与发震时间见表 2。
| 异常信号 | 发震日期 | 经度/(°) | 纬度/(°) | 震级/Mw | 发震地点 |
| 2011-03-08 | 2011-03-10 | 97.98 | 24.64 | 5.5 | 云南 |
| 2011-09-05 | 2011-08-09 | 98.73 | 24.98 | 5.1 | 云南 |
| 2017-08-13 | 2017-08-08 | 103.89 | 33.21 | 6.5 | 四川 |
| 2017-12-14 | 2017-11-23 | 107.94 | 29.40 | 4.9 | 重庆 |
4 结论
本文基于GNSS坐标时间序列异常波动的相对强度指数提出一种蠕滑形变信息的自动检测方法。由于PCA存在过度滤波现象，且不能估计高阶项信息，因此采用ICA方法进行时空滤波。该方法在提高坐标序列信噪比的同时克服了测站分布及信号相关性的影响，为我国地壳形变监测提供一种可行的方法。仿真试验的结果表明，当蠕滑信号至少和噪声水平相当时，可有效探测地表瞬态位移。根据Akutan地区PBO 7个测站GNSS坐标数据的计算结果探测到2008年2月13日的蠕滑信号。由于该地区火山活动频繁，因此推断此次蠕滑事件可能与火山岩的强烈运动有关。对四川地区陆态网18个GNSS观测站资料处理后发现2011年3月8日和2017年12月14日出现两次高概率异常信号。两个信号出现在地震前后2—21天，由此推断可能与震前应力积累造成的地表异常移动和震后余滑密切相关。降低异常信号概率的阈值之后，探测出另外两个异常信号，均出现在较大地震发生后不久，这些现象是否与地震的产生有关尚待进一步研究。该探测方法能够有效探测区域断层蠕滑现象，为今后断层参数的估计和断层滑移时空分布反演提供重要的参考信息。但也存在一定的不足，由于该方法是通过计算各测站每个历元发生异常的CDF值来探测异常信号，而异常信号的CDF值并不具有连续性，因此尚无法判断其持续时间。今后将通过研究异常信号前后各历元异常概率的变化特征及其之间的相关性来进一步探索该问题的解决方法。
| [1] | FENG Lujia, NEWMAN A V. Constraints on continued episodic inflation at Long Valley Caldera, based on seismic and geodetic observations[J]. Journal of Geophysical Research, 2009, 114(B6): B06403. DOI:10.1029/2008JB006240 |
| [2] | KING N E, ARGUS D, LANGBEIN J, et al. Space geodetic observation of expansion of the San Gabriel Valley, California, aquifer system, during heavy rainfall in winter 2004-2005[J]. Journal of Geophysical Research, 2007, 112(B3): B03409. DOI:10.1029/2006JB004448 |
| [3] | MILLER M M, MELBOURNE T, JOHNSON D J, et al. Periodic slow earthquakes from the Cascadia subduction zone[J]. Science, 2002, 295(5564): 2423. DOI:10.1126/science.1071193 |
| [4] | NISHIMURA T. Short-term slow slip events along the Ryukyu trench, southwestern Japan, observed by continuous GNSS[J]. Progress in Earth and Planetary Science, 2014(1): 22. DOI:10.1186/s40645-014-0022-5 |
| [5] | MITSUI Y. Interval modulation of recurrent slow slip events by two types of earthquake loading[J]. Earth, Planets and Space, 2015, 67(1): 56. |
| [6] | OKUTANI T, IDE S. Statistic analysis of swarm activities around the Boso Peninsula, Japan:slow slip events beneath Tokyo Bay?[J]. Earth, Planets and Space, 2011, 63(5): 419–426. DOI:10.5047/eps.2011.02.010 |
| [7] | RIEL B, SIMONS M, AGRAM P, et al. Detecting transient signals in geodetic time series using sparse estimation techniques[J]. Journal of Geophysical Research:Solid Earth, 2014, 119(6): 5140–5160. DOI:10.1002/2014JB011077 |
| [8] | WALLACE L M, BEAVAN J, BANNISTER S, et al. Simultaneous long-term and short-term slow slip events at the Hikurangi subduction margin, New Zealand:implications for processes that control slow slip event occurrence, duration, and migration[J]. Journal of Geophysical Research:Solid Earth, 2012, 117(B11): B11402. DOI:10.1029/2012JB009489 |
| [9] | MAVROMMATIS A P, SEGALL P, JOHNSON K M. A decadal-scale deformation transient prior to the 2011 Mw9.0 Tohoku-Oki earthquake[J]. Geophysical Research Letters, 2014, 41(13): 4486–4494. DOI:10.1002/2014GL060139 |
| [10] | HOLTKAMP S, BRUDZINSKI M R. Determination of slow slip episodes and strain accumulation along the Cascadia margin[J]. Journal of Geophysical Research:Solid Earth, 2010, 115(B4): B00A17. DOI:10.1029/2008JB006058 |
| [11] | OHTANI R, MCGUIRE J J, SEGALL P. Network strain filter:a new tool for monitoring and detecting transient deformation signals in GPS arrays[J]. Journal of Geophysical Research:Solid Earth, 2010, 115(B12): B12418. DOI:10.1029/2010JB007442 |
| [12] | JI K H, HERRING T A. A method for detecting transient signals in GPS position time-series:smoothing and principal component analysis[J]. Geophysical Journal International, 2013, 193(1): 171–186. DOI:10.1093/gji/ggt003 |
| [13] | WALWER D, CALAIS E, GHIL M. Data-adaptive detection of transient deformation in geodetic networks[J]. Journal of Geophysical Research:Solid Earth, 2016, 121(3): 2129–2152. DOI:10.1002/2015JB012424 |
| [14] | CHAMOLI A, LOWRY A R, JEPPSON T N. Implications of transient deformation in the northern basin and range, western United States[J]. Journal of Geophysical Research:Solid Earth, 2014, 119(5): 4393–4413. DOI:10.1002/2013JB010605 |
| [15] |
徐克科, 伍吉仓, 吴伟伟.
基于GNSS时空数据的瞬态无震蠕滑信息检测[J]. 地球物理学报, 2015, 58(7): 2330–2338.
XU Keke, WU Jicang, WU Weiwei. Detection of transient aseismic slip signals from GNSS spatial-temporal data[J]. Chinese Journal of Geophysics, 2015, 58(7): 2330–2338. |
| [16] |
徐克科, 伍吉仓, 王成.
利用GNSS位移时空序列进行断层无震蠕滑特征分析[J]. 武汉大学学报(信息科学版), 2015, 40(9): 1247–1252.
XU Keke, WU Jicang, WANG Cheng. Analysis of fault aseismic slip feature based on GNSS displacement time-space series[J]. Geomatics and Information Science of Wuhan University, 2015, 40(9): 1247–1252. |
| [17] |
赵斌, 王东振, 谭凯, 等.
顾及黏弹性松弛效应的2015年尼泊尔Mw7.9地震震后早期余滑模型[J]. 大地测量与地球动力学, 2018, 38(3): 239–243.
ZHAO Bin, WANG Dongzhen, TAN Kai, et al. Afterslip model following the 2015 Nepal Mw7.9 earthquake considering viscoelastic relaxation displacements[J]. Journal of Geodesy and Geodynamics, 2018, 38(3): 239–243. |
| [18] |
谭凯, 乔学军, 杨少敏, 等.
汶川地震GPS形变约束的破裂分段特征及滑移[J]. 测绘学报, 2011, 40(6): 703–709.
TAN Kai, QIAO Xuejun, YANG Shaomin, et al. Rupture characteristic and slip constrained by GPS coseismic deformation induced by the Wenchuan earthquake[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(6): 703–709. |
| [19] | Dong D, FANG P, Bock Y, et al. Spatiotemporal filtering using principal component analysis and Karhunen-Loeve expansion approaches for regional GPS network analysis[J]. Journal of Geophysical Research:Solid Earth, 2006, 111(B3): B03405. DOI:10.1029/2005JB003806 |
| [20] | SHEN Yunzhong, LI Weiwei, XU Guochang, et al. Spatiotemporal filtering of regional GNSS network's position time series with missing data using principle component analysis[J]. Journal of Geodesy, 2014, 88(1): 1–12. DOI:10.1007/s00190-013-0663-y |
| [21] | MING Feng, YANG Yuanxi, ZENG Anmin, et al. Spatiotemporal filtering for regional GPS network in China using independent component analysis[J]. Journal of Geodesy, 2017, 91(4): 419–440. DOI:10.1007/s00190-016-0973-y |
| [22] |
姜卫平, 周晓慧.
澳大利亚GPS坐标时间序列跨度对噪声模型建立的影响分析[J]. 中国科学:地球科学, 2014, 44(11): 2461–2478.
JIANG Weiping, ZHOU Xiaohui. Effect of the span of Australian GPS coordinate time series in establishing an optimal noise model[J]. Science China:Earth Sciences, 2014, 44(11): 2461–2478. |
| [23] |
张倩倩, 归庆明, 宫轶松.
卫星多故障探测和识别的独立分量分析法[J]. 测绘学报, 2017, 46(6): 698–705.
ZHANG Qianqian, GUI Qingming, GONG Yisong. Multiple satellite faults detection and identification based on the independent component analysis[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(6): 698–705. DOI:10.11947/j.AGCS.2017.20160079 |
| [24] |
罗飞雪, 戴吾蛟, 唐成盼, 等.
参考经验模态分解-独立分量分析及其在GPS多路径误差处理中的应用[J]. 测绘学报, 2012, 41(3): 366–370.
LUO Feixue, DAI Wujiao, TANG Chengpan, et al. EMD-ICA with reference signal method and its application in GPS multipath[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(3): 366–370. |
| [25] | WILDER J W. New concepts in technical trading systems[M]. London: Trend Research, 1978. |
| [26] | OCHI T, KATO T. Depth extent of the long-term slow slip event in the Tokai district, central Japan:a new insight[J]. Journal of Geophysical Research:Solid Earth, 2013, 118(9): 4847–4860. DOI:10.1002/jgrb.50355 |
| [27] | OZAWA S. Long-term slow slip events along the Nankai trough subduction zone after the 2011 Tohoku earthquake in Japan[J]. Earth, Planets and Space, 2017, 69: 56. DOI:10.1186/s40623-017-0640-4 |
| [28] | CROWELL B W, BOCK Y, LIU Zhen. Single-station automated detection of transient deformation in GPS time series with the relative strength index:a case study of Cascadian slow slip[J]. Journal of Geophysical Research:Solid Earth, 2016, 121(12): 9077–9094. DOI:10.1002/2016JB013542 |
| [29] | JI K H, HERRING T A. Transient signal detection using GPS measurements:transient inflation at Akutan Volcano, Alaska, during early 2008[J]. Geophysical Research Letters, 2011, 38(6): L06307. DOI:10.1029/2011GL046904 |
| [30] | JI K H, YUN S H, RIM H. Episodic inflation events at Akutan Volcano, Alaska, during 2005-2017[J]. Geophysical Research Letters, 2017, 44(16): 8268–8275. DOI:10.1002/2017GL074626 |