利用GPS和GRACE分析四川地表垂向位移变化

引用本文

丁一航, 黄丁发, 师悦龄, 等. 2018. 利用GPS和GRACE分析四川地表垂向位移变化. 地球物理学报, 61(12): 4777-4788, doi: 10.6038/cjg2018L0654.

Ding Y H, Huang D F, Shi Y L, et al. 2018. Determination of vertical surface displacements in Sichuan using GPS and GRACE measurements. Chinese J. Geophys. (in Chinese), 61(12): 4777-4788, doi: 10.6038/cjg2018L0654.

利用GPS和GRACE分析四川地表垂向位移变化

丁一航, 黄丁发, 师悦龄, 姜中山, 陈婷

西南交通大学地球科学与环境工程学院, 成都 611756

收稿日期 2017-10-17, 2018-06-07 收修定稿

基金项目: 国家自然科学基金（41374032）资助

第一作者简介: 丁一航, 男, 硕士研究生, 研究方向为GPS导航与定位.E-mail:2011301610267@whu.edu.cn

通讯作者: 黄丁发, 男, 博士, 教授, 主要从事GNSS地基增强系统、地壳形变与大地测量反演研究.E-mail:dfhuang@swjtu.edu.cn

摘要：陆地水储量的季节性变化是导致地表周期性负荷形变位移的主要因素，有效地剔除地表位移中的陆地水储量影响，是获取地壳构造垂向运动的必要过程.四川地处青藏高原东边缘，地形分区明显，境内以长江水系为主，水资源丰富，研究四川地区地表负荷形变位移，有助于分析陆地水储量的时空分布特性及地壳构造形变信息.本文利用研究区域内59个CORS站的GPS观测数据，计算了CORS站点的垂向位移，并将其与GRACE所得相应结果进行对比分析.结果显示，GPS和GRACE所得垂向位移时间序列的振幅大小整体相符，但存在明显的相位差.GPS站点振幅最大值为12.7 mm，对应HANY站，最小值为1.5 mm，对应SCMX站.GRACE所得的地表垂向位移振幅大小均为3~4 mm，且最大位移集中出现在7-9月份；而GPS站点出现最大位移的月份和地形相关，东部盆地、西北部高原和南部山地分别出现在7-8月份、10-11月份和10月份.GPS站点时间序列中的周年项与陆地水的季节性变化强相关，为了讨论陆地水储量对GPS站点位移的影响，本文利用改进的总体经验模态分解方法（MEEMD：Modified Ensemble Empirical Mode Decomposition），从GPS垂向位移时间序列中提取出周年项及约2年的年际变化项.发现利用MEEMD获取的周年项改正原始GPS时间序列时可使其WRMS（Weight Root Mean Square）减少量减小约26%，结果优于最小二乘拟合方法提取的GPS周年项改正效果，验证了MEEMD方法在GPS坐标时间序列处理中的可行性及有效性.

关键词: 地表垂向位移陆地水储量连续运行参考站时间序列分析改进的总体平均经验模态分解

Determination of vertical surface displacements in Sichuan using GPS and GRACE measurements

DING YiHang, HUANG DingFa, SHI YueLing, JIANG ZhongShan, CHEN Ting

Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, Chengdu 611756, China

Abstract: Seasonal variations of terrestrial water storage (TWS) mainly result in periodical fluctuations of crustal deformation, eliminating the effects from TWS on surface displacement is accordingly conducive to highlight regional tectonic movements. Sichuan Province lies close to the eastern margin of the Qinghai-Tibetan Plateau and has a rising terrain feature from its east to the west. Moreover, abundant river systems, such as Yangtze River system, provide Sichuan with adequate water resource. We combine GPS-derived and GRACE-derived vertical displacement time series to explore the spatial-temporal distribution characteristics of TWS and the crustal tectonic movements in Sichuan. The result shows that the amplitudes of vertical displacements derived by those two ways are accorded very well, but there exist phase delays. The GPS-determined amplitudes show a maximum of 12.7 mm at station HANY and minimum of 1.5 mm at station SCMX, respectively. Due to a low spatial resolution, GRACE-derived amplitudes present comparable magnitude of 3~4 mm and the largest displacements are in between July and September. However, the time when the largest displacements of GPS happen correspond with topography, the eastern basin and hills, northwestern plateau and southern mountains possess its own largest displacement in July to August, October to November and October, separately. Seasonal TWS changes are relevant to GPS-derived vertical displacement variations. So, in order to illustrate the effect of TWS changes on GPS sites' movements, we use the Modified Ensemble Empirical Mode Decomposition (MEEMD) to extract annual term and a 2-year periodic signal from GPS vertical displacement time series. Consequently, we find that when using MEEMD-derived annual terms to correct GPS time series, the WRMS (Weight Root Mean Square) reduction can decrease roughly 26%, which is superior to least square method, suggesting that the MEEMD method is feasible in dealing with GPS time series.

Keywords: Earth surface vertical displacement Terrestrial water storage CORS Time series analysis Modified EEMD

0 引言

GPS坐标时间序列包含陆地水迁徙引起的周期性地壳垂向形变信号(Wahr et al., 1998; Tesmer et al., 2011).如何利用GRACE时变重力场信息, 有效地去除GPS坐标时间序列中由陆地水迁徙所引起的非构造形变影响, 最大限度地凸显真实的构造运动信息, 成为近年来GPS地壳形变监测研究领域的研究热点(刘任莉等, 2013).

Davis等(2004)利用GRACE时变重力场研究了亚马逊流域陆地水储量变化, 并将其与GPS坐标时间序列进行比较, 发现GRACE所得陆地水储量变化与GPS垂向时间序列具有很好的一致性.Fu等(2013)利用GPS和GRACE观测资料分别比较了亚马逊流域和亚洲南部陆地水储量变化对各自区域内地表水平位移的影响, 发现GPS和GRACE计算所得的地表水平位移之间也具有很强的相关性.Fu等(2012)利用GPS和GRACE数据研究了尼泊尔和喜马拉雅地区的长期垂向形变, 并发现当利用GRACE计算得到的负荷形变位移改正GPS坐标时间序列时, 可使GPS站点残差时间序列的WRMS(Weighted Root Mean Square)减少量减小约45%.Wahr等(2013)利用GPS水平位移, 结合负荷弹性形变理论, 研究了负荷源的地域分布情况及量级大小.Pan等(2016)利用主成分分析方法提取GPS观测数据中的共模误差, 并将其与GRACE计算的地表垂向位移进行比较, 发现两者具有很好的一致性, 说明GPS垂向坐标时间序列中, 共模误差主要是由于陆地水负荷的季节性变化所致; 王林松等(2014), 比较了由GPS和GRACE计算得到的喜马拉雅地区负荷变化; 其他学者, 结合GPS和GRACE数据, 分别对比研究了山西、滇西地区的负荷形变情况, 均得出GPS和GRACE结果符合较好的结果(刘任莉等, 2013;盛传贞, 2014).

在分析GPS和GRACE垂向位移时间序列时, 通常采用最小二乘拟合方法, 这样的数据处理策略, 往往会忽略时间序列中蕴含的其他细节信号.Bryant等(2014)通过实验证明, 总体/补充的总体平均经验模态分解方法(EEMD/CEEMD: Ensemble/Complement Ensemble Empirical Mode Decomposition)可有效提取GRACE所得的陆地水储量变化时间序列中的周年变化.Pan等(2015)利用EEMD方法, 提取出了GPS站点的周年变化及其年际变化, 发现GPS站点垂向位移时间序列中的年际变化可能与大气负荷有关, 并得出在剔除周年变化和年际变化后, 原始序列WRMS减少量明显减小.

本文利用四川省59个CORS站的GPS观测资料, 计算站点的垂向位移时间序列, 并与GRACE得到的垂向位移时间序列进行对比, 以分析陆地水储量的时空分布特性及地表运动情况; 利用改进的总体经验模态分解方法(MEEMD: Modified Ensemble Empirical Mode Decomposition)(郑近德等, 2013), 分别从GPS观测值和GRACE观测值中提取出了周年信号与年际信号, 并分析了其年际信号的成因, 为合理解释陆地水循环周期和分析地表质量迁移提供重要依据.最后利用MEEMD提取的GPS站点垂向位移周年信号, 分析了汶川地震前后, 沿龙门山断裂带部分站点的运动速率变化情况.

1 数据处理 1.1 GPS观测数据及预处理

本文采用四川省地震局和中国大陆环境监测网络(CMONOC: Crustal Movement Observation Network of China)提供的四川省境内共59个CORS站的GPS观测数据, 站点分布如图 1所示.GPS数据采样间隔为30 s, 数据解算使用麻省理工学院提供的GAMIT软件, 解算策略如下:使用消电离层无关组合LC消除1阶电离层效应; 高阶电离层误差使用CODE(Center for Orbit Determination in Europe)提供的单日电离层产品进行改正; 绝对天线相位中心校正文件使用从IGS(International GNSS Service)网站获取的结果; 海潮负荷效应的改正使用了FES2004海洋潮汐模型; 先验对流层延迟投影函数使用VMF1(Vienna Mapping Function); 大气负荷和非潮汐海洋影响认为可以忽略.在获得单日松弛约束解后, 利用GLOBK滤波得到单日解的站点时间序列, 最后选择中国境内和周边13个IGS站点作为联系站点, 使用赫尔默特6参数转换模型将站点时间序列结果转换到IGS08框架下.

图 1 四川省CORS站分布图 SCB指四川盆地, SGB指松潘—甘孜地块, SYR指川滇地区. Fig. 1 Distribution of Sichuan CORS sites SCB, SGB and SYR are Sichuan Basin, Songpan-Garze Block and Sichuan-Yunnan Region, respectively.

1.2 GRACE数据与处理

本文GRACE时变重力场采用美国德克萨斯大学空间研究中心(CSR: Center of Space Research)提供的RL05模型, 时间跨度为2005年1月—2015年9月, 共112个月数据.缺省数据采用三次样条插值补齐后, 共获得129个月数据.GRACE数据预处理时, 用卫星激光测距(SLR)获得的C₂₀项替换GRACE数据的C₂₀项(Cheng and Tapley, 2004), 并向球谐系数中加入Swenson等(2008)计算的地心改正项, 再使用去相关滤波(P3M8), 结合350 km高斯滤波, 削弱球谐系数的高阶项噪声(Swenson and Wahr, 2002).按式(1)计算GRACE探测到的地表垂向位移(王林松等, 2014):

(1)

其中, dr(θ, φ)为地表径向位移; h_l和k_l为l阶负荷勒夫数; P_{(l, m)}(cosθ)为完全规格化缔合勒让德函数; C_lm, S_lm为去均值的月重力场球谐系数, W_l为高斯平滑核函数.

2 GPS与GRACE垂向位移时间序列对比分析 2.1 最小二乘方法提取季节性位移

为了获取准确的季节性位移, 在最小二乘拟合中, 需要考虑去除坐标时间序列的线性趋势, 以及地震和仪器更换导致的跳变.以成都(CHDU)站为例, GPS垂向坐标时间序列拟合结果见图 2所示, 图 2a去掉了线性趋势项, 图 2b为去趋势和跳变结果.拟合公式如下:

(2)

图 2 成都(CHDU)站垂向坐标位移时间序列 (a)时间序列去趋势结果; (b)时间序列去除跳变项结果. Fig. 2 Vertical displacement time series of CHDU (a) The detrended time series; (b) The detrended time series with jumps removed.

其中y₀为观测序列, t_i为观测时间, H(t_i-T_gj)为由于地震和仪器更换导致的跳变项, e(t_i)为残差.振幅为, 周年相位为tan^-1(d₁/c₁).

由于GPS站点时间序列中的周年项与陆地水负荷的季节性变化具有很强的相关性, 因此, 对比了GPS与GRACE的垂向位移时间序列周年项拟合振幅, 结果如图 3所示, 图 3a中箭头指向垂向位移最大时所在月份, 箭头向右指向1月, 月份按顺时针增加.结合图 3a与图 3c可以看出, GPS与GRACE垂向位移振幅整体符合较好; 图 3a还显示出, GPS站点垂向位移最大值出现的时间与地形有关.

图 3 GPS与GRACE垂向位移对比结果及GPS站点水平向与垂向运动趋势 (a)和(c)是GPS与GRACE的周年垂向位移时间序列振幅对比, (a)中箭头指向站点位移最大时所在月份(箭头向右指向1月, 月份按顺时针增大), 箭头长短代表振幅大小.(b)是GPS站点垂向位移速率, 绿色和白色框标记内容表示速率较大站所在区域.(d)是GPS水平速度场. Fig. 3 Comparison between GPS-derived and GRACE-derived vertical and horizontal displacements (a) and (c) Denotes the comparison of the annual amplitudes between GPS and GRACE vertical displacement time series. The length of arrows shows the value of annual amplitudes (mm) and the arrows point to the month when the largest displacements appear (clockwise from the right); (b) Presents the trends of the vertical displacements of GPS sites, Green and white annotation indicates the location of sites with large velocities; (d) Denotes the GPS measured horizontal velocity field.

四川省横跨我国第一、第二两级阶梯, 地势由东往西急剧抬升.川西高山高原属青藏高原东南缘, 整体隆起; 东部为海拔约200~750 m的山脉环绕盆地丘陵区域; 川西南, 亦称攀西地区, 主要由山地和丘陵构成.复杂的地形地貌, 使得四川的气候呈现多样化, 表现在大面积区域内地带性气候类型被地形气候类型所取代.四川盆地为亚热带地带性气候类型, 而山地气候则表现出不同区域的垂直分异结果.气候的多变性, 主要体现在温度和降雨的区域分化上.川西、四川盆地和攀西地区, 全年高温月份都集中在夏季, 7—9月份, 海拔越高, 全年均温越低; 降雨量则呈由东向西减少趋势, 盆地夏季降雨丰富.川西北高山高原高寒气候, 年降水量较少(陈青青等, 2017).降雨量是陆地水储量变化的主要因素来源之一, 地表垂向位移主要受陆地水储量影响, 故推测, 导致最大地表位移出现在不同月份的原因, 四川盆地主要是由于夏季降雨等导致的.而川西和攀西山地区由于温度较低, 全年降雨量低于四川盆地, 所以其地表垂向位移出现的月份主要受温度影响, 包括基岩的热胀冷缩、冰雪消融、冰川覆盖等影响.

据图 3a所示, 东部四川盆地(SCB)GPS站点垂向位移最大值出现在7—8月份, 南部川滇地区(SYR)出现在约10月份, 西北部则为10—11月份.而GRACE由于其空间分辨率较低, 整体出现在7—9月份.图 3b是GPS站点垂向位移变化的趋势图, 可以看出, 四川盆地西缘, BEIC、MYAN、YAAN、QLAI等沿龙门山断层、北川—映秀断裂带分布的站点, 运动趋势较明显, 平均速率超过2.5±0.25 mm·a^-1, 其中SCMX站最大达到7.98±0.6 mm·a^-1, BEIC为-4.85±0.32 mm·a^-1.南部SCML、XICH等站运动趋势也较为突出.具体数值结果如表 1所示.

表 1 四川省CORS站点信息 Table 1 Relevant information of Sichuan CORS sites

站点	经度(°E)	纬度(°N)	时间段	振幅/(mm)	速率/(mm·a^-1)	GPS振幅/(mm)	GPS速率/(mm·a^-1)
站点	经度(°E)	纬度(°N)	时间段	GPS/GRACE	GPS/GRACE	East/North	East/North
ANYU	105.3304	30.1038	2010.88—2015.81	5.3/3.79	0.66/-0.26	1.22/0.75	35.71/-8.36
BEIC	104.439	31.7895	2010.80—2015.63	1.63/3.25	-4.85/-0.21	0.84/1.99	39.4/-8.88
CAXI	105.9252	31.7396	2010.90—2015.81	5.47/3.15	-0.51/-0.23	0.96/0.89	34.19/-9.48
CHDU	104.0645	30.6394	2006.00—2015.74	5.45/3.73	-0.23/-0.22	3.06/0.78	35.25/-8.72
GANL	102.77	28.9591	2010.84—2015.80	5.82/4.63	0.47/-0.17	1.36/0.74	37.12/-10.13
GYAO	102.2666	27.7962	2007.83—2009.99	7.51/5.29	6.64/-0.14	0.83/0.69	36.24/-17.42
HANY	102.6324	29.3492	2010.88—2015.77	12.73/4.47	3.51/-0.17	0.66/3.06	36.97/-9.4
HONY	102.5482	32.7924	2010.70—2015.80	3.24/2.98	3.51/-0.17	1.32/1.32	44.37/-8.22
JIGE	105.5595	32.2054	2010.74—2015.81	5.15/3.01	0.40/-0.22	0.85/0.77	33.44/-8.97
JYAN	104.5448	30.3885	2006.00—2015.74	4.18/3.78	0.36/-0.23	0.71/0.35	34.37/-9.1
KAIY	102.112	27.8821	2007.83—2010.9	9.63/5.28	-14.32/-0.13	0.69/0.50	45.8/-17.11
LESH	103.7553	29.5648	2006.25—2015.74	4.28/4.21	-0.21/-0.21	1.13/0.62	35.34/-8.8
LUZH	105.4136	28.8724	2006.25—2015.74	3.37/4.23	-0.04/-0.27	0.76/0.92	35.62/-9.43
MAEK	102.2287	31.8991	2010.70—2015.84	5.15/3.37	1.88/-0.16	1.14/0.75	42.45/-10.68
MEIG	103.134	28.335	2008.27—2015.72	4.95/4.86	-0.75/-0.18	1.13/0.91	37.13/-9.69
MNIU	101.7445	30.5999	2006.67—2012.83	6.30/3.99	-1.15/-0.13	1.05/0.27	38.56/-9.41
MOXI	102.1151	29.6497	2010.88—2015.77	6.21/4.40	3.27/-0.14	0.96/1.59	39.51/-10.87
MYAN	104.7264	31.4399	2006.25—2015.74	3.74/3.36	-0.55/-0.22	1.14/0.64	35.03/-8.8
NEIJ	105.1186	29.62	2006.25—2015.74	4.23/4.00	0.35/-0.26	0.84/0.29	34.7/-8.67
PIWU	104.5471	32.4142	2010.88—2013.88	4.06/3.01	-2.68/-0.20	1.53/1.33	39.37/-7.31
PIXI	103.757	30.9102	2006.00—2015.74	3.29/3.66	-0.66/-0.21	0.71/0.84	35.27/-9.94
QICH	105.2275	32.5939	2010.88—2015.77	5.55/2.89	-0.61/-0.20	0.60/0.54	37.02/-5.15
QIME	101.185	31.0229	2006.74—2012.83	5.10/3.83	-1.33/-0.12	0.94/0.67	40.31/-10.66
QLAI	103.306	30.3544	2006.00—2015.74	4.06/3.94	-2.88/-0.19	1.47/0.62	35.04/-9.25
RENS	104.1029	30.2004	2006.00—2015.74	4.56/3.91	-0.32/-0.22	1.05/0.71	35.09/-8.9
ROXI	104.4338	29.4579	2007.37—2015.74	3.95/4.16	-0.91/-0.24	1.16/0.40	34.53/-8.23
SCBZ	106.7446	31.8408	2012.09—2015.87	5.11/3.05	0.81/-0.24	1.20/0.43	34.21/-9.37
SCDF	101.1227	30.9779	2012.00—2015.87	4.26/3.86	0.72/-0.12	0.76/1.13	43.14/-12.23
SCGU	105.8997	32.4248	2014.96—2015.87	7.30/2.91	1.97/-0.21	2.26/0.65	35.93/-11.46
SCGY	105.8523	32.4389	2012.00—2014.96	4.55/2.91	-1.15/-0.21	1.28/0.10	34.0/-8.04
SCGZ	100.0188	31.6103	2012.00—2015.87	3.88/3.62	-0.51/-0.09	0.80/0.59	47.64/-9.59
SCJL	101.5004	29.0081	2012.00—2015.87	5.47/4.82	0.19/-0.11	1.35/1.12	39.63/-17.83
SCJU	104.5157	28.1792	2012.00—2015.87	4.45/4.65	-1.06/-0.24	1.18/1.32	34.23/-8.03
SCLH	100.672	31.3878	2012.00—2015.87	5.08/3.69	0.2/-0.11	0.71/0.57	45.76/-10.9576
SCLT	100.2181	29.9916	2012.36—2015.87	5.47/4.45	1.00/-0.06	1.11/1.23	44.11/-15.0
SCMB	103.5336	28.8404	2012.00—2015.87	6.24/4.56	-0.68/-0.2	1.30/0.53	35.28/-9.35
SCML	101.2753	27.9294	2012.00—2015.87	9.41/5.44	-6.07/-0.1	1.85/0.82	49.44/-19.3
SCMN	102.1738	28.3328	2012.00—2015.87	5.73/5.04	-1.14/-0.14	1.36/0.34	38.84/-16.20
SCMX	103.8496	31.6712	2012.00—2015.87	1.55/3.34	7.98/-0.2	1.63/0.68	44.46/-6.91
SCNC	105.8822	30.9793	2012.00—2015.87	5.58/3.42	0.77/-0.25	1.17/0.83	34.76/-9.73
SCNN	102.7164	27.057	2012.00—2015.87	7.5/5.54	-0.59/-0.16	0.96/0.68	37.30/-16.83
SCPZ	101.7439	26.5032	2012.00—2015.87	6.26/6.08	-0.95/-0.13	0.63/0.89	35.71/-17.64
SCSM	102.3537	29.2294	2012.00—2015.87	7.0/4.57	1.18/-0.15	1.67/0.63	38.25/-12.05
SCSN	105.5621	30.5076	2012.00—2015.87	5.36/3.62	0.38/-0.26	1.59/0.30	33.97/-8.94
SCSP	103.5824	32.6484	2012.00—2015.87	2.59/2.89	0.23/-0.19	1.35/1.20	41.23/-9.85
SCTQ	102.765	30.0736	2012.00—2015.87	3.76/4.12	-0.32/-0.17	1.35/1.73	35.63/-9.4
SCXC	99.8032	28.9374	2012.00—2015.87	6.15/5.12	0.9/-0.03	1.02/0.55	41.35/-16.88
SCXD	102.437	28.3004	2012.00—2015.87	5.97/5.01	-1.32/-0.15	0.80/0.96	38.58/-14.11
SCXJ	102.3721	31.0004	2012.00—2015.87	4.37/3.75	-0.04/-0.16	1.25/0.70	40.67/-5.2
SCYX	102.5115	28.6509	2012.00—2015.87	4.79/4.82	-0.56/-0.15	1.16/1.28	37.84/-13.70
SCYY	101.5126	27.4319	2012.00—2015.87	5.95/5.66	-0.87/-0.11	0.97/1.74	38.55/-17.51
TAGO	101.5258	30.3261	2006.66—2012.82	5.42/4.14	-1.32/-0.12	0.74/0.52	42.54/-16.66
WARI	101.1254	30.8812	2006.66—2012.82	7.63/3.91	-0.12/-0.11	1.02/0.51	44.05/-14.82
WCXI	103.4626	30.985	2010.88—2015.76	2.51/3.66	1.48/-0.2	2.42/0.64	42.55/-13.23
WENC	103.6001	31.4741	2010.88—2015.75	4.18/3.44	-2.72/-0.2	1.66/0.50	51.53/-24.92
XICH	102.3479	27.8094	2007.83—2015.72	5.66/5.27	-1.4/-0.14	1.82/0.26	36.06/-16.98
YAAN	103.0106	29.981	2006.22—2015.74	4.58/4.13	-0.18/-0.6	1.72/0.72	36.55/-9.24
YANT	105.3859	31.2208	2010.92—2015.80	5.69/3.38	0.5/-0.24	1.05/0.78	34.54/-7.85
YBIN	104.5962	28.7983	2006.24—2015.74	4.12/4.39	-0.9/-0.24	1.31/0.80	34.16/-9.19
ZHJI	104.5458	31.0063	2006.00—2015.74	3.74/3.54	-0.53/-0.23	0.83/0.51	35.05/-8.71
ZTON	105.1691	31.6505	2010.71—2015.81	5.09/3.24	1.97/-0.23	0.33/1.52	28.51/-10.77

表 1 四川省CORS站点信息 Table 1 Relevant information of Sichuan CORS sites

图 3d反映的是GPS站点在IGS08框架下的水平速度场, 可以看出, 四川整体呈现右旋运动趋势, 受华南块体影响, 四川盆地内站点的位移速度, 略小于川西高原及攀西地区.值得一提的是, WENC(汶川)呈现出往北方向运动趋势, 导致其水平位移与其他站点方向不一致, 本文推测是震后断层活动局部变形引起的, 其具体原因还有待进一步研究.

2.2 MEEMD方法提取季节性位移

GPS站点垂向位移时间序列主要由趋势项和非线性项组成, 非线性项主要包括周年项和年际项, 这些成分的准确计算对确定GPS站点的垂向运动速率具有很大影响.最小二乘拟合算法主观性较强, 通常仅能得到时间序列中具有明确周期的周年和半周年项, 不能全面的反映坐标时间序列原有的非平稳特性.为了克服上述缺点, 本文引入了非平稳时间序列处理MEEMD方法, 并将其应用到GPS垂向时间序列的周期信号提取中.MEEMD不仅保证了原始时间序列的完整性, 还能提取出GPS垂向时间序列中的非线性分量.

MEEMD方法是以总体完备经验模态分解方法为基础, 利用排列熵作为信号的随机性判别标准的信号分解方法.MEEMD方法避免了EEMD/CEEMD方法中不必要的集成平均, 使得分解结果中更小几率出现虚假信号.MEEMD实施步骤如下(郑近德等, 2013):

(1) 利用CEEMD分解方法, 得到IMF1.

(2) 利用排列熵检验IMF1的随机性, 熵值阈值取H_p=0.6, 嵌入维度取m=6, 时间延迟取λ=1.

(3) 如果IMF1的熵值H>H_p, 则说明该IMF1具有随机分布特性, 继续利用CEEMD算法提取下一个IMF, 重复步骤1—3.直到对于某个IMF, H < H_p, 即IMF不具有随机分布特性, 则自该IMF以下的IMF, 直接采取EMD算法进行提取.

以成都站(CHDU)为例, 对其进行MEEMD分解, 结果如图 4所示.其中IMF6主要包含半周年项, IMF7主要包含周年项, Trend为趋势项, Res为各IMF叠加值与原始数据之间的差值.利用快速傅里叶变换计算IMF6—IMF8的频率, 发现成都站GPS垂向位移时间序列中主要包含了约2.01年的年际变化对应图中的IMF8, 利用MEEMD分析GLDAS-Noah水文模型和GRACE垂向位移, 发现其分别包含约1.8和1.9年的陆地水循环周期, 说明GPS站点中约2.01年的年际变化主要与陆地水循环有关.后文为了比较MEEMD和最小二乘拟合方法, 将IMF6—IMF9叠加作为周年项.

图 4 成都(CHDU)站GPS垂向位移时间序列的MEEMD分解结果 Fig. 4 Decomposition results of GPS vertical displacement time series in CHDU by using MEEMD

选取时间跨度较长的8个站点, 如JYAN、LESH等, 利用MEEMD分别对其GPS和GRACE数据进行分解, 结果如图 5.可以看出, 利用MEEMD方法分解得到的站点GPS/GRACE时间序列周年项更能反映GPS/GRACE数据的时变特点, 更加符合站点的真实运动情况.

图 5 四川省部分站点GPS和GRACE垂向位移时间序列MEEMD分解所得周年项结果 Fig. 5 Annual terms of GPS and GRACE vertical displacement time series extracted by MEEMD in representative sites

2.3 MEEMD与最小二乘拟合周年项结果对比

为了比较MEEMD方法与最小二乘拟合方法提取信号的能力, 本文利用WRMS(Weighted Root Mean Square)减少量WRMS_reduction作为指标, 分别计算GPS垂向位移时间序列与MEEMD和最小二乘拟合提取的周年项之间的相关性, WRMS_reduction计算公式如下(Fu et al., 2012):

(6)

图 6是MEEMD和最小二乘拟合计算的WRMS_reduction结果对比.图 6a是利用MEEMD计算的WRMS_reduction结果, 其均值为26%, 最大值为48%, 最小值为7%;图 6b则是最小二乘拟合计算的WRMS_reduction结果, 其均值为18%, 最大为45%, 最小为2%.两种方法计算的WRMS_reduction最大和最小值均分别对应HANY站和PIWU站.比较图 6a和图 6b, 利用MEEMD得到的WRMS_reduction明显大于最小二乘拟合方法的结果, 进一步计算可知, MEEMD方法计算的WRMS_reduction比最小二乘拟合计算结果平均要大40%, 说明利用MEEMD提取的周年项相比于最小二乘拟合结果更加反映原始数据GPS站点垂向坐标时间序列的周年特性.

图 6 MEEMD方法与最小二乘拟合方法提取的GPS垂向位移周年项的WRMS_reduction比较 (a) MEEMD提取的周年项计算的WRMS_reduction结果; (b)最小二乘拟合提取的周年项计算结果. Fig. 6 Comparison of WRMS_reduction of the GPS′s vertical displacement annual terms derived by MEEMD method and conventional fitting method (a) The WRMS after removing the annual term derived by MEEMD from GPS vertical displacement time series; (b) The result with the annual term removed, which is derived by least square method.

GPS站点垂向位移时间序列周年项和年际项的提取, 会影响其运动趋势和地震导致的跳变项的计算精度.以CHDU站为例, 探究汶川地震对研究区域内GPS站点运动的影响.首先利用最小二乘方法计算GPS站点的趋势项及跳变项, 并将其从GPS垂向位移时间序列中剔除; 然后从剩余的时间序列中, 直接减掉利用MEEMD提取的周年项结果; 最后再以发生地震的时间为节点, 对余下时间序列残差进行多线段拟合, 再次提取上述原因导致的趋势项和跳变项, 将第二次得到的趋势项和跳变项与之前的提取结果相叠加, 得到最终较为精确的趋势项和跳变项数值.CHDU站多线段拟合结果如图 7所示.

图 7 成都(CHDU)GPS垂向时间序列残差多线段拟合结果 Fig. 7 Multi-linear fitting results from the residual of vertical displacement time series of CHDU

选取数据记录时间跨越汶川地震前后的GPS站点, 采取前述MEEMD结合最小二乘拟合方法提取站点运动趋势与地震导致的位移, 得到更为精确的GPS站点运动信息, 提取结果见表 2.从表 2中看出, 地震前后, GPS站点垂向运动速率均有明显变化, TAGO、WARI、MNIU等海拔较高站点受局部形变影响, 位移速率有所增大, 其他站点的速率基本呈现减小趋势, 具体原因有待进一步研究.

表 2 受汶川地震影响的站点运动信息 Table 2 Information of sites affected by the Wenchuan earthquake

从表 2还可得到, 站点受汶川地震导致的地表垂向和水平向位移大小, 与GPS观测站点距汶川的距离成反比.PIXI、CHDU、QLAI、YAAN、MYAN等站点由于汶川地震导致的垂向位移均偏大, 其中: PIXI最大达到了-86.8 mm, 其E向位移更是达到-566 mm.QLAI、YAAN站的垂向位移要大于CHDU、ZHJI等站, 而水平位移则相反, 且与其他大多数站不同的是, 这两个站点的垂向位移要显著大于水平位移.本文推测是由于这两个GPS站点位于龙门山断裂带上, 而龙门山断裂带构造运动较为活跃所致.

图 8显示的是GPS站点垂向位移运动速率(趋势)内插图, 从图中可以得出, 四川盆地内的站点运动速率更多地呈现减小趋势.图 8a是2008年汶川地震前站点的运动速率图, 图 8b图则是地震后的运动速率图, 其中红色箭头所示站点为内插所用站点.可以看出, 汶川地震前后, 部分GPS站点的运动速率发生变化, 其中最为明显的是ROXI站, 其运动速率从震前的-3.5 mm·a^-1到震后的-1.6 mm·a^-1.

图 8 部分站点汶川地震前后速率变化图 (a)震前运动速率; (b)震后运动速率; 红色箭头为站点运动速率变化的站点. Fig. 8 Partial sites′ rates before and after Wenchuan earthquake (a) and (b) denote the rates before and after the earthquake, separately.Red arrows are the sites whose rates changed.

上述分析结果表明, 从原始时间序列中剔除利用MEEMD方法提取的周年项, 然后以地震发生时间为节点进行多线段拟合, 提取出地震前后站点运动速率的变化, 相比于整体的运动趋势, 更能反映局部时间内, 站点的运动情况, 有利于分析地震前后地表的运动情况.

3 结论

本文对比分析了四川省境内59个CORS站的GPS观测数据与GRACE得到的垂向位移时间序列; 利用MEEMD方法从GPS垂向坐标时间序列中提取出了年际变化和周年变化, 并以WRMS作为判别标准比较了MEEMD方法和最小二乘方法对原始GPS时间序列拟合的效果.

受局部地形与气候影响, 四川地区, 部分GPS站点垂向位移周年振幅较GRACE稍大, 但整体相符, 说明陆地水储量变化是导致四川地区地表位移呈现出季节性变化特点的主要原因.GPS振幅最大值为12.7 mm, 对应HANY站, 最小值为1.5 mm, 对应SCMX站, 而GRACE由于空间分辨率较低原因, 研究区域内, 垂向位移振幅均约为3~4 mm.研究结果还发现, GPS站点发生最大位移的时间与地形强相关, 东部盆地、西北部高原和南部山地垂向位移分别在7—8月份、10—11月份和10月份达到最大.

为了讨论陆地水储量对GPS站点垂向位移的影响, 利用MEEMD方法对GPS垂向位移时间序列进行分解, 提取了更为准确的周年项与年际变化, 并发现周期约为2.01年的年际变化与陆地水储量循环周期相关.进一步, 当利用MEEMD获取的周年项改正GPS垂向位移时间序列时, 可使其WRMS减少量减小约26%, 而最小二乘拟合方法提取的GPS周年项改正结果仅减小约18%, 验证了MEEMD方法在GPS时序处理中的有效性.最后发现, 利用MEEMD结合多线段拟合的方法, 可从GPS垂向位移时间序列中获得更加清晰的地表位移量及运动趋势, 为研究地震前后地壳运动情况提供参考.

致谢感谢四川省地震局和中国陆态网提供的四川省CORS站GPS观测数据, 感谢德克萨斯州立大学提供的GRACE L2数据, 本文的图片由GMT软件生成, 在此一并表示感谢.

References

Chen Q Q, Tang Z Y, Yang L, et al. 2017. Analysis of air temperature and precipitation in Sichuan. Journal of Chengdu University of Information Technology (in Chinese), 32(2): 200-207.

Cheng M K, Tapley B D. 2004. Variations in the Earth's oblateness during the past 28 years. Journal of Geophysical Research:Solid Earth, 109(B9): B09402. DOI:10.1029/2004JB003028

Davis J L, Elósegui P, Mitrovica J X, et al. 2004. Climate-driven deformation of the solid Earth from GRACE and GPS. Geophysical Research Letters, 31(24): L24605. DOI:10.1029/2004GL021435

Fu Y N, Freymueller J T. 2012. Seasonal and long-term vertical deformation in the Nepal Himalaya constrained by GPS and GRACE measurement. Journal of Geophysical Research:Solid Earth, 117(B3): B03407. DOI:10.1029/2011JB008925

Fu Y N, Argus D F, Freymueller J T, et al. 2013. Horizontal motion in elastic response to seasonal loading of rain water in the Amazon Basin and monsoon water in Southeast Asia observed by GPS and inferred from GRACE. Geophysical Research Letters, 40(23): 6048-6053. DOI:10.1002/2013GL058093

Liu R L, Li J C, Jiang W P, et al. 2013. Comparing vertical surface displacements using GRACE and GPS over Shanxi province. Geomatics and Information Science of Wuhan University (in Chinese), 38(4): 426-430.

Loomis B D, Luthcke S B. 2014. Optimized signal denoising and adaptive estimation of seasonal timing and mass balance from simulated GRACE-Like regional mass variations. Advances in Adaptive Data Analysis, 6(1): 145003. DOI:10.1142/S1793536914500034

Pan Y J, Shen W B, Ding H, et al. 2015. The quasi-biennial vertical oscillations at global gps stations:Identification by ensemble empirical mode decomposition. Sensors, 15(10): 26096-26114. DOI:10.3390/s151026096

Pan Y J, Shen W B, Hwang C, et al. 2016. Seasonal mass changes and crustal vertical deformations constrained by GPS and GRACE in northeastern Tibet. Sensors, 16(8): 1211. DOI:10.3390/s16081211

Sheng C Z, Gan W J, Liang S M, et al. 2014. Identification and elimination of non-tectonic crustal deformation caused by land water from GPS time series in the western Yunnan province based on GRACE observation. Chinese Journal of Geophysics (in Chinese), 57(1): 42-52. DOI:10.6038/cjg20140105

Swenson S, Wahr J. 2002. Methods for inferring regional surface-mass anomalies from Gravity Recovery and Climate Experiment (GRACE) measurements of time-variable gravity. Journal of Geophysical Research:Solid Earth, 107(B9): 2193. DOI:10.1029/2001JB000576

Swenson S, Chambers D, Wahr J. 2008. Estimating geocenter variations from a combination of GRACE and ocean model output. Journal of Geophysical Research:Solid Earth, 113(B8): B08410. DOI:10.1029/2007JB005338

Tesmer V, Steigenberger P, Van Dam T, et al. 2011. Vertical deformations from homogeneously processed GRACE and global GPS long-term series. Journal of Geodesy, 85(5): 291-310. DOI:10.1007/s00190-010-0437-8

Wahr J, Molenaar M, Bryan F. 1998. Time variability of the Earth's gravity field:Hydrological and oceanic effects and their possible detection using GRACE. Journal of Geophysical Research:Solid Earth, 103(B12): 30205-30229. DOI:10.1029/98JB02844

Wahr J, Khan S A, Van Dam T, et al. 2013. The use of GPS horizontals for loading studies, with applications to northern California and southeast Greenland. Journal of Geophysical Research:Solid Earth, 118(4): 1795-1806. DOI:10.1002/jgrb.50104

Wang L S, Chen C, Zou R, et al. 2014. Using GPS and GRACE todetect seasonal horizontal deformation caused by loading of terrestrial water:A case study in the Himalayas. Chinese Journal of Geophysics (in Chinese), 57(6): 1792-1804. DOI:10.6038/cjg20140611

Zheng J D, Cheng J S, Yang Y. 2013. Modified EEMD algorithm and its applications. Journal of Vibration and Shock (in Chinese), 32(21): 21-26, 46.

陈青青, 汤志亚, 杨玲, 等. 2017. 四川气温和降水量特征分析. 成都信息工程大学学报, 32(2): 200-207.

刘任莉, 李建成, 姜卫平, 等. 2013. 联合GRACE与GPS比较山西省垂向地表形变. 武汉大学学报:信息科学版, 38(4): 426-430.

盛传贞, 甘卫军, 梁诗明, 等. 2014. 滇西地区GPS时间序列中陆地水载荷形变干扰的GRACE分辨与剔除. 地球物理学报, 57(1): 42-52. DOI:10.6038/cjg20140105

王林松, 陈超, 邹蓉, 等. 2014. 利用GPS与GRACE监测陆地水负荷导致的季节性水平形变:以喜马拉雅山地区为例. 地球物理学报, 57(6): 1792-1804. DOI:10.6038/cjg20140611

郑近德, 程军圣, 杨宇. 2013. 改进的EEMD算法及其应用研究. 振动与冲击, 32(21): 21-26, 46. DOI:10.3969/j.issn.1000-3835.2013.21.004


地球物理学报 2018, Vol. 61 Issue (12): 4777-4788	PDF