融合GPS与强震仪数据实时监测瞬时地壳形变

陈向阳,于金池,葛 建,侯勇涛

(1. 南通职业大学建筑工程学院，江苏 南通 226007； 2. 南通大学交通与土木学院，江苏 南通 226019； 3. 上海华测导航技术股份有限公司，上海 201702)

近些年来，GPS被广泛地运用于长周期地壳形变、火山监测和台风预警[1-7]。这些监测主要依靠于毫米级的GPS坐标时间序列，与此同时，GPS站的采样间隔通常设置为30 s[8]。地震及其引起的火山爆发、滑坡和泥石流等对人类的生产和生活影响巨大，因此，快速准确地监测地震具有尤为重要的意义。随着GPS定位精度与方法的成熟，高采样率的GPS数据(如1、10、50 Hz)已经被用于监测由地震引起的瞬时地壳形变。高频GPS被成功运用于监测2002年Denali地震、2008年汶川地震和2011年日本宫城9.0级地震等[9-11]。之后，高频GPS结合其他手段也被成功应用于2003年Miyagi地震、2009年L’Aquila地震和2015年尼泊尔地震[12-13]。鉴于高频GPS与强震仪互有优缺点，因此可以融合GPS与强震仪的数据来实时监测地震引起的动态形变。本文以2016年意大利中部的Mw 6.0地震为例，探讨融合GPS数据和强震仪数据实时监测中等强度地震引起瞬时地壳形变的可行性与优势。

1 数据与方法

1.1 GPS数据

2016年8月24号，意大利中部发生了Mw 6.0级地震，震中位于Perugia，震源深度为8.1 km。目前，意大利境内已有完善的高密度GNSS观测网络，主要包括RING、ISPRA、DPC、Regione Lazio、Regione Abruzzo、Leica ITALPOS和Topcon NETGEO，使探测中等强度引起的同震动态形变成为可能。此次地震引起的同震形变成功地被近场的10 Hz高频GPS和200 Hz强震仪所记录。为了融合GPS与强震仪数据实时监测地震引起的地壳形变，笔者挑选了3组位置近似相同的GPS/强震仪组合进行分析(即RM33/MTER、GSA/INFN和AMT/AMAT)，GPS与强震仪之间的相对距离少于1.1 km。此次所用的高频GPS数据由INGV提供，强震仪数据由美国强震运动监测中心(CESMD)提供。如图1所示。

1.2 数据融合方法

鉴于强震仪记录的加速度二次积分求解的位移波形不理想，而高频GPS由于采样率的问题不能够更为细致地描述位移波形，许多学者尝试将GPS与强震仪数据结合来监测地震的同震动态形变。文献[14]于2011年提出了利用Kalman滤波将GPS与强震仪数据融合来获取准确和细致的位移波形。

将以上Kalman滤波应用于测站地表位移波形融合时，假设其在N、E、U 3个方向上均为一维的运动形式，可用一阶线性微分方程来描述其运动状态，即

(1)

对上述一阶线性微分方程离散化，其状态方程系数阵会发生变化，若离散数据采样间隔为τ，采用离散GPS位移观测量zk和强震仪加速度ak进行Kalman滤波融合，可建立如下的离散状态方程与观测方程

(2)

(3)

1.3 交叉小波变换

首先介绍连续小波变换(CWT)，一组时间序列Xn的连续小波变换如下

(4)

WXY=WXWY*

(5)

式中，Y*为WY的共轭矩阵。平均相位角可以更好地描述两个时间序列的相位关系。一组相位角(ai,i=1,2，…，n)的平均相位角可以表示为

(6)

对am进一步处理可以更直观地描述两组时间序列之间的相关性。交叉小波相似值可以用下式来计算

ρ=cosαm

(7)

式中，ρ=1表示两组时间序列完全相关；ρ=0表示两组时间序列没有明显的相关性；ρ=-1表示两组时间序列呈负相关。

2 结果与分析

2.1 GPS与强震仪位移波形对比

由于此次地震持续时间较短，仅选取震前10 s到震后60 s这个时间段的坐标时间序列进行分析。与GPS站不同的是，强震仪记录的往往是地震引起的加速度变化值，需要一个二次积分的过程才能够得到震时的位移波形。对强震仪记录的加速度进行二次积分后，利用SeismoSignal软件进行基线校正和低通滤波。图2为GPS位移波形和强震仪二次积分得到的位移波形的对比。由图2可知，GPS/强震仪组合RM33/MTER在水平分量上符合较好，两种手段获取的位移波形差别在2 cm以内；由于垂直方向的精度问题，GPS/强震仪组合在垂直方向上位移波形符合度较差。

2.2 GPS与强震仪数据融合

Kalman滤波可以利用GPS和强震仪的优点，进而补充GPS和强震仪各自的缺陷。融合后的位移波形可以更为细致地描述地震引起的动态地表形变，同时也可以探测地震引起的同震永久阶跃。以GPS/强震仪组合RM33/MTER为例，利用Kalman滤波进行数据的融合。从图3可以清楚看到，Kalman滤波获取的位移波形与GPS位移波形符合很好，与此同时，地震引起的同震阶跃在N方向上可以明显看到。GPS位移数据可以作为融合数据的整体趋势，而细节部分则用高采样率的强震仪数据进行描述。Kalman滤波获取的位移波形相对于GPS位移波形更为平滑，说明融合GPS与强震仪数据能够获取更为精确的位移波形。尽管GPS垂向方向精度较差，但是在强震仪数据的辅助下Kalman滤波仍然能够获取精确的位移波形。

2.3 交叉小波分析

为了进一步验证Kalman滤波位移波形与GPS、强震仪位移波形的相关性，笔者采用交叉小波变换进行任意2个位移波形之间的相关性分析。在交叉小波谱中，黑色粗线包围的区域表示通过了95%的显著性检验。两组时间序列的相对相位关系如交叉小波谱中的箭头所示(向右表示相位关系相同，向左表示相位关系相反，向上表示相位落后90°变化，向下表示相位超前90°变化)。

图4为GPS/强震仪组合RM33/MTER的强震仪与Kalman滤波位移波形、强震仪与GPS位移波形、GPS与Kalman滤波位移波形的交叉小波谱。对于GPS与强震仪位移波形的交叉小波谱(如图4(b)所示)，共振区域内的箭头几乎全部朝右，表明GPS与强震仪位移波形在地震发生期间呈明显的正相位关系。GPS与强震仪在共振区域的交叉小波相似值为0.99，表明二者在地震时呈明显的正相关性。对于强震仪、GPS与Kalman滤波位移波形的交叉小波谱(如图4(a)和(b)所示)，共振区域内的部分箭头向下有一定的偏移，说明Kalman滤波位移波形比GPS、强震仪位移波形的相位有一定的超前。图4(a)和(b)的共振区域的交叉小波相似值均在0.8以上，说明Kalman滤波波形与GPS、强震仪位移波形均呈较强的正相关性。对于3个组合，共振区域的小波相似值，GPS与强震仪位移波形的交叉小波相似值均在0.9以上，而Kalman滤波与GPS、强震仪位移波形的交叉小波相似值均在0.8以上。交叉小波分析表明，Kalman滤波获取的位移波形与GPS、强震仪位移波形均呈较强的正相关性，表明Kalman滤波可能获取准确的位移波形。

3 结 语

本文以2016年意大利Mw 6.0级地震为例，基于3组位置近似相同的GPS/强震仪组合，探讨了利用Kalman滤波融合GPS与强震仪数据监测中等强度地震引起的瞬时地壳形变的可行性与优势。结果表明Kalman滤波能够有机地结合GPS与强震仪数据的优点，融合后的位移波形精度更高，并能够对地震引起的动态地壳形变的细节部分进行较好地描述。然后利用交叉小波变换，对融合后的位移波形与GPS、强震仪的位移波形在时频域的相关性进行了分析。交叉小波分析表明，Kalman滤波位移波形与GPS、强震仪位移波形呈较强的正相关性，说明Kalman滤波确实获取了比较稳定可靠的位移波形。由于GPS与强震仪位置不完全一致，Kalman滤波获取的位移波形与GPS、强震仪位移波形会有略微的相位差异。