2. 中国科学院大学, 北京市玉泉路19号甲, 100049;
3. 国家卫星气象中心, 北京市中关村南大街46号, 100081
Montagner等[1]发现,在2011年Tohoku-Oki MW9.0地震主震发生后P波到达之前,KAM台超导重力仪记录到量级在0.1 μGal的瞬时重力变化,证明重力观测网可以缩短地震预警时间,并且更快更准确地估计震级大小。Imanishi等[2]对2003年Tokachi-Oki MW8.0地震期间3个近震台站的超导重力仪数据进行分析,发现3个台站全都记录到0.1 μGal量级的同震重力效应。雷湘鄂等[3]对秘鲁MS7.9地震后武汉基准台超导重力仪数据进行分析,成功观测到0 S0~0 S32的全部基频振荡。由此可知,在地震发生后,连续重力观测资料包含了大量的有效信息,能够清晰地记录到同震及震后变化。其中,同震信号是永久变形引起的重力效应和地震波(体波)的叠加,对于这样的混频非稳态信号,如何快速、准确地实现两者的有效分离对于快速反演震源参数及地震预警具有重要意义。针对这一问题,本文将正则化平滑因子方法引入连续重力观测数据处理中,以实现信号的简单快速分离。
1 数据处理与方法 1.1 地震概况及数据情况选取2011-03-11 Tohoku-Oki MW9.0地震、2014-02-12新疆于田MW6.9地震及2014-10-07云南景谷MW6.6地震,对地震前后的近震台站数据进行筛选。由于3次地震的信号能量都很大,有些仪器出现记录中断、错误等现象,因此选取数据质量较好的3个近震台站数据作为研究对象(表 1),台站及震中位置见图 1。
连续重力观测资料包含了固体潮汐信号、同震重力效应及地震波信号,为有效提取地震信号,本文采用国际地球潮汐研究中心推荐的Tsoft预处理软件[4]将原始观测资料中由仪器突跳和电脉冲等产生的尖峰、扰动等错误信号删除。选用合成潮汐法,首先利用理论公式计算出观测点的理论潮汐值,然后从观测数据中扣除理论潮汐值,进而得到同震信号,结果见图 2。从图中可以看到,合成潮汐法很好地去除了原始观测数据中的低频固体潮汐信号,得到的同震信号结果只包含了同震重力效应和地震波。
同震信号包含了同震重力效应和地震波,由于同震重力效应是由永久变形和密度变化引起的,表现为非稳态信号,使得同震信号(如图 2中最上边的点线所示)也是非稳态信号。当需要利用同震重力效应反演震源参数时,地震波会影响其结果的准确性;反之,当需要利用地震波反演断层破裂过程时,同震重力效应也会对结果产生影响。因此,快速准确地实现两者的有效分离对于反演震源参数和破裂过程及进行地震预警都有重要意义。但是,在进行信号分离时,如果直接采用最小二乘法拟合,参数估计时就会受到突跳(同震重力效应)的影响,导致拟合得到的周期信号失真。为解决这一问题,本文引入了正则化平滑因子方法[5-8]。
将非稳态的同震重力效应表示为:
$ \boldsymbol{z}_{\text {coseismic }}=\boldsymbol{H} \theta+\boldsymbol{\varepsilon} $ | (1) |
式中,H∈Rn×m为观测矩阵,θ∈Rm为回归参数,ε为观测误差,n为观测数据采样点的个数,m为回归参数的个数。选用正则化最小二乘准则方法求解回归参数的估计值
$ \hat{\theta}_{\lambda}=\arg \min\limits_{\theta}\left\{\|\boldsymbol{H} \theta-\boldsymbol{z}\|^{2}+\lambda^{2}\left\|D_{d}(\boldsymbol{H} \theta)\right\|^{2}\right\} $ | (2) |
式中,z为重力残差时间序列,λ为正则化平滑因子,Dd为第d次导数操作的离散逼近。式(2)在普通最小二乘方法的基础上增加了一个条件,其目的是使式(2)的解在满足‖Hθ-z‖取得极小值的同时,‖Dd(Hθ)‖的取值也尽可能小。在常见的正则化方法中,正则化参数的选择是难点,而本文要解决的物理问题具有一定的特殊性:地震波的传播频率固定,与震级、震中距等无关,因此可以很容易地选定截止频率,进而确定正则化参数,也解决了在应用正则化方法时平滑因子的选取受人为因素影响的问题。
式(2)的解可表示为:
$ \begin{aligned} &\hat{\theta}_{\lambda}=\left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{H}+\lambda^{2} \boldsymbol{H}^{\mathrm{T}} {D}_{{d}} \boldsymbol{H}\right)^{-1} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{z}\\ &\boldsymbol{\hat{z}}_{\text {cosessmic }}=\boldsymbol{H} \hat{\theta}_{\lambda} \end{aligned} $ | (3) |
式中,
图 3为KAM台的观测数据应用正则化平滑因子方法分离得到的结果。从图 3(a)中不仅可以清晰地看到MW9.0主震引起的地震波信号,同时也可以清晰地分辨出多次强余震引起的地震波。由于短时间内连续发生多次强余震,分离得到的同震重力效应(图 3(b))实际上是由多次地震引起的共同效应,当信号趋于稳定后进行定量估算,计算结果见表 2。
图 4为于田台的连续重力观测资料分离得到的结果,图 5为孟连台的连续重力观测资料分离得到的结果。与Tohoku-Oki地震相比,于田地震和景谷地震的震级较小,主震发生后短时间内都没有发生较大的余震,因此信号分离的效果非常明显。由于在地震波的冲击下,弹簧重力仪有可能出现掉格效应,为了判断仪器在震后是否出现了掉格,同时也为了验证结果的准确性,将计算结果与理论模拟值进行对比,其中,理论模拟所采用的震源参数和断层模型来源于GCMT计划[9-10]。从表 2中可以看到,KAM台的计算结果与理论模拟结果有较大差别,这是因为理论模拟结果为Tohoku-Oki地震主震引起的同震重力效应,而本文的计算结果是由Tohoku-Oki地震主震及余震共同引起的同震重力效应;于田台的计算结果与理论模拟结果吻合较好;对于孟连台,其理论模拟结果太小,易受震源参数和断层模型误差的影响,可靠性不高,因此这里没有给出其理论模拟的结果。
综上所述,利用正则化平滑因子方法能在快速有效地获取高频地震波的同时,完整地保留同震重力效应。因此,连续重力观测资料可以作为传统地震学观测仪器的有力补充,并为准确获取大震震源参数及震源破裂过程提供一种新的技术手段。
3 结语本文选取了Tohoku-Oki地震KAM台、新疆于田地震于田台及云南景谷地震孟连台的连续重力观测资料,通过引入正则化平滑因子方法,快速、有效地获取了同震重力效应和地震波信号,并得到以下结论:
1) 连续重力观测资料可以清晰地记录同震及震后变化,其中同震信号混叠了同震重力效应和地震波信号,两种信号相互叠加、相互影响,需要进行信号分离;
2) 采用正则化平滑因子方法可以快速有效地分离同震重力效应和地震波信号,根据实际物理问题本身的特点,可以快速确定正则化因子,分离结果受人为因素影响较小;
3) 对于大震引起的近震台站同震重力效应,本文的计算结果与理论模拟的结果具有较好的一致性。
致谢: 感谢重力与形变数据共享分中心提供数据;感谢德国波茨坦地球科学中心汪荣江老师提供QSSPSTATIC程序。
[1] |
Montagner J P, Juhel K, Barsuglia M, et al. Prompt Gravity Signal Induced by the 2011 Tohoku-Oki Earthquake[J]. Nature Communications, 2016(7): 13349
(0) |
[2] |
Imanishi Y, Sato T, Higashi T, et al. A Network of Superconducting Gravimeters Detects Submicrogal Coseismic Gravity Changes[J]. Science, 2004, 306(5695): 476-478 DOI:10.1126/science.1101875
(0) |
[3] |
雷湘鄂, 许厚泽, 孙和平. 利用超导重力观测资料检测地球自由振荡[J]. 科学通报, 2002, 47(18): 1432-1436 (Lei Xiang'e, Xu Houze, Sun Heping. The Free Oscillation of the Earth Detected by Superconducting Gravity Observations[J]. Chinese Science Bulletin, 2002, 47(18): 1432-1436)
(0) |
[4] |
Vauterin P. Tsoft: Graphical and Interactive Software for the Analysis of Earth Tide Data[C]. 13th Int Sympos on Earth Tides, Brussels, 1998
(0) |
[5] |
Tarvainen M P, Ranta-Aho P O, Karjalainen P A. An Advanced Detrending Method with Application to HRV Analysis[J]. IEEE Transactions on Biomedical Engineering, 2002, 49(2): 172-175 DOI:10.1109/10.979357
(0) |
[6] |
Gersch W. New Directions in Time Series Analysis[M]. New York: Springer, 1993
(0) |
[7] |
Karjalainen P A. Regularization and Bayesian Methods for Evoked Potential Estimation[D]. Kuopio: University of Kuopio, 1997
(0) |
[8] |
詹金刚, 王勇, 史红岭, 等. 应用平滑先验信息方法移除GRACE数据中相关误差[J]. 地球物理学报, 2015, 58(4): 1135-1144 (Zhan Jingang, Wang Yong, Shi Hongling, et al. Removing Correlative Errors in GRACE Data by the Smoothness Priors Method[J]. Chinese Journal of Geophysics, 2015, 58(4): 1135-1144)
(0) |
[9] |
Dziewonski A M, Chou T A, Woodhouse J H. Determination of Earthquake Source Parameters from Waveform Data for Studies of Global and Regional Seismicity[J]. Journal of Geophysical Research:Solid Earth, 1981, 86(B4): 2825-2852 DOI:10.1029/JB086iB04p02825
(0) |
[10] |
Ekstr m G, Nettles M, Dziewonski A M. The Global CMT Project 2004-2010:Centroid-Moment Tensors for 13017 Earthquakes[J]. Physics of the Earth and Planetary Interiors, 2012, 200-201: 1-9 DOI:10.1016/j.pepi.2012.04.002
(0) |
2. University of Chinese Academy of Sciences, A19 Yuquan Road, Beijing 100049, China;
3. National Satellite Meteorological Center, 46 South-Zhongguancun Street, Beijing 100081, China