时序形变资料的多核函数滤波方法研究及应用

梁洪宝, 刘雪龙， 郭炳辉

(中国地震局第一监测中心，天津 300180)



时序形变资料的多核函数滤波方法研究及应用

梁洪宝, 刘雪龙， 郭炳辉

(中国地震局第一监测中心，天津 300180)

针对多年时序形变观测资料有效信息提取复杂的问题，对基于多核函数的滤波方法进行研究，得到以下有益结论：(1)当核函数指数为0.5，光滑因子为0.003时，10天及以上核点间隔的滤波模型单位权中误差最小；(2)核点间隔控制滤波信息频谱的高低, 间隔越大频谱信息越低，反之则频谱信息越高；(3)因数据缺失部分造成核点减少，当连续减少2个以上时滤波失败，当连续减少2个时数据缺失部分滤波出现失真，当减少1个时滤波效果不受影响；(4)通过对GPS时序资料、定点形变时序资料和非构造形变时序资料的滤波应用，获取不同频谱的信息，验证了本文方法的稳定性和可靠性。

多核函数； 时序形变资料； 核点间隔； 滤波

0 引言

为监测中国大陆断层的构造活动，自1976年唐山7.8级地震以来，在重点断层带建立了定点形变地震观测台站，如唐山地震台至今积累了三十多年的单天形变观测资料[1]；上世纪90年代，随着GPS技术兴起并成功应用于地壳形变监测中，我国启动了国家重大科学工程——中国地壳运动观测网络[2](简称“网络工程”)，建立了25个GPS基准站(后来陆续增加到30个)，并于1998年开始运行；2007年启动的国家重大科技基础设施——中国大陆构造环境监测网络[3](简称“陆态网络”)，对网络工程GPS基准站进行加密，共建成260个测站，并于2010年下半年开始运行，至今积累了数年的单天观测资料。众所周知，上述形变观测资料均含有大气、非潮汐海洋、积雪和土壤湿度等非构造形变信息，对于地壳形变分析而言，需要对其含有的非构造形变信息进行剔除[4]。根据多年非构造地球物理模型可以获取目标区域的非构造形变序列，以此对观测资料进行修正。

时序形变原始观测资料不能直接用于形变分析，需对其含有的有效信息进行提取，即根据研究需要对包含多种信息的原始资料进行滤波，以得到较为客观、清晰的信息。时序资料的滤波方法有多种[5]，但能完全满足需求的方法并不多，如样条函数方法和最小二乘配置法从理论上讲均可以满足要求，但在实际应用中却存在计算不稳定、耗时较长等问题[6]。利用高斯型函数作为多核函数的核函数对GPS时间序列进行滤波，其计算过程较为简易[7]。 因此，本文在此基础上探讨基于指数型核函数的多核函数模型对时间跨度较长的时序形变资料的滤波效果，确定其较优的光滑因子，并将其应用于上述三种时序形变资料的滤波处理中。

1 多核函数模型及滤波原理

多核函数法是由美国Hardy教授于1977年提出并建议用于地壳垂直形变分析的一种纯数学的方法，其基本思想是任何一个圆滑的数学表面总可以用一系列有规则的数学表面的总和以任意精度逼近[8-9]。其数学表达式为：

(1)

式中:Q为核函数;a为待定系数;m为核函数个数。核函数可任意选用，一般表示为：

(2)

式中:δ2称为光滑因子，用来对核函数进行调整;k为非零实数，一般取0.5、-0.5或1.5。

根据多核函数法的基本思想，假设时序形变的变化是由多条曲线组成的一条圆滑曲线，因此可以将其应用于以时间变化的一维自变量的领域，函数表达式可以描述为：

(3)

相应的核函数可以描述为：

(4)

1.1 光滑因子的确定

以长春IGS站2005-2011年垂向去线性趋势运动序列为例(以下类同)，核点间隔分别为10、20、30、60、90和120天，光滑因子在0～0.5取值，步长取值为0.000 5，滑动计算1 000次，并以验后单位权中误差为评价指标，探讨不同核点间隔下光滑因子与单位权中误差间的关系，结果如图1所示。从图1可以看出，在理想的单位权中误差下，不同的核点间隔光滑因子的取值均有一个相对稳定的区间，这个区间随核点间隔的增大而增大。核点间隔为10天时，光滑因子稳定区间为[0.001,0.005],为了应用简便快捷，本文取中间值0.003，以满足核点间隔为10天及以上所有情况。

1.2 核点间隔与滤波信息的关系

确定光滑因子为0.003后，分别取核点间隔为20、30、60、90、120和150天，探讨不同核点间隔与滤波信息的关系，结果如图2所示。从图2可以看出，核点间隔大小控制滤波信息的频谱，间隔越大频谱信息越低，间隔越小频谱信息越高，因此在实际应用中要根据研究对象合理确定核点间隔，以获取所需信息。

1.3 数据连续性缺失对滤波效果的影响

因各种因素影响，在实际监测工作中观测资料难免出现连续性缺失。为稳健起见，将长春站2006年10—12月的连续观测数据删除，以1.1节确定的光滑因子为0.003，核点间隔分别为20、30、39、46、60和70天，探讨数据连续性缺失对滤波效果的影响，结果如图3所示。与图2对比，从图3可以看出，核点间隔为20天和30天时数据缺失部分滤波出现错误；核点间隔为39天和46天时数据缺失部分滤波出现失真；剩余两种情况滤波正常。经过分析，图3(a)和(b)中数据缺失的同时也造成核点分别减少了4个和3个，图3(c)和(d)中核点丢失了2个，图3(e)和(f)中核点丢失了1个，由此可以认为，核点连续丢失数不能超过2个。当观测数据连续性缺失时若连续丢失2个以上核点则滤波失败；若连续丢失2个核点则数据缺失部分滤波出现失真；若丢失1个核点则滤波效果不受影响。

图1 不同核点间隔下光滑因子与单位权中误差间的关系Fig.1 Relationship between the smooth factor and mean square error of unit weight with different kernel point intervals

图2 核点间隔与滤波信息的关系Fig.2 Relationship between the kernel point interval and filtering information

图3 数据缺失时不同核点间隔的滤波效果Fig.3 The filtering effect of different kernel point intervals in the case of when missing data missing

2 时序形变资料的滤波应用

根据上述滤波方法及参数设置，分别对GPS时序资料、定点形变时序资料和非构造形变时序资料进行滤波应用。

2.1 GPS时序资料

随着网络工程和陆态网络项目的实施，连续GPS站产出了多年时序资料，如图4(黑色圆点表示计算值，红色曲线表示本文方法的滤波值，下图表示方法相同)所示，测站分布如图5所示。KMIN(昆明)站积累了16年左右的坐标序列，从图4(a)中滤波曲线(核点间隔为90天)可以看出，该站水平方向主要以线性运动为主，垂向以周年运动为主。为探测该站的精细运动，对其线性运动进行剔除，结果如图4(b)所示。通过滤波曲线可以看出，该站自2004年开始南向运动加剧，积累的大量能量经2004年12月26日苏门答腊9.3级地震影响得到一定程度的释放；自2007年开始南北向运动出现“闭锁”，由此积累的能量经2008年5月12日汶川8.0级地震的影响再次进行释放，构造运动经过后续四年左右时间的调整期逐渐恢复到原来的运动特征。其东西向也因受汶川地震影响出现西向运动加剧的特征，一直持续至今，垂向运动受地震影响较小；除构造运动外，滤波曲线还显示出年周期的非构造运动，垂向最为显著。

图4 KMIN站坐标序列Fig.4 Coordinate series of KMIN station

图5 地震活动与GPS测站分布Fig.5 Seismic activity and the GPS station distribution

基线序列是用来分析块体间构造运动的拉张或挤压特征，应变序列是用来分析区域块体的受力应变特征。基于GPS连续站序列可以计算测站基线和块体应变序列[12]，计算结果和滤波结果(核点间隔为180天)分别如图6和图7所示。图6中LHAZ(拉萨)—DLHA(德令哈)—DXIN(鼎新)三站位置近似处于一条直线上，并北东向横跨青藏高原东部，基线缩短速率拉萨—鼎新约为22 mm/a，拉萨—德令哈约为17 mm/a，德令哈—鼎新约为5 mm/a。这说明青藏高原东部缩短速率由南向北递减，这也与青藏高原南部受印度洋北向推挤、北部受塔里木块体和鄂尔多斯块体阻挡的基本受力格局相吻合，期间受2001年11月14日昆仑山8.1级地震的影响，拉萨—德令哈和德令哈—鼎新分别有10 mm的缩短和伸长。图7中由YANC(盐池)-DXIN(鼎新)-XNIN(西宁)三站得到该区域的应变系列，该区域地处青藏高原东北缘，是印度与欧亚两大板块碰撞作用由近南北方向向北东、东方向转换的重要场所，也是中国大陆东西及南北构造结合部位和重要的构造转换区域。图7中其东西向和南北向均处于持续挤压状态特征也说明了这一点。2001年11月14日昆仑山8.1级地震和2010年4月14日7.1级地震对该区域均有一定影响，分别表现出东西方向和南北方向的拉张变形特征。

图6 GPS基线序列Fig.6 GPS baseline series

图7 GPS应变序列Fig.7 GPS strain series

2.2 定点形变时序资料滤波应用

唐山地震台始建于1978年，位于唐山市路南区复兴路原第十中学院内，是唐山地区唯一的大地形变台站。1976年唐山7.8级地震后，地壳升降和错动的巨大形变大体上都以发震断层为交界面，由于深部形变与强烈地震波的共同作用，在震区形成了大量的总体走向为N30°～40° E呈雁状排列的地裂缝，其中发震断层附近尤为明显。唐山地震台正处于地裂缝最为发育的典型地段，台站布设有水准和基线测线各4条，长度24～48 m不等，以不同角度与断层交汇。测站分布如图8所示，基线测量顺序是J1、J2、J3、J4、J1，水准测量顺序是S1、S2、S3、S4、S1。水准S1-S2和基线J3-J4观测结果如图9所示，从图中可以看出，通过不同核点间隔(分别为30天、182天和365天)可以获取不同频谱的信息，除趋势性变化外，还可以清晰展现短期动态变化。在高差曲线变化中曲线上升表示断裂正断活动，下降则表示逆断活动；在基线曲线变化中，上升则表示断裂拉张运动，下降则表示挤压运动。

图8 唐山台观测点分布Fig.8 Distribution of observation points at Tangshan station

图9 断层高差与基线序列Fig.9 The fault height difference and baseline series

2.3 非构造形变时序资料滤波应用

由于垂直形变观测资料是多种信息的综合体现，且较水平运动弱。在进行地壳垂直形变分析时，必须对其含有的非构造形变信息进行剔除，这种非构造形变主要是由非潮汐海洋、大气、积雪和土壤湿度等非构造地球物理因素产生的。本文计算了KMIN(昆明)GPS站的非构造垂向形变，非构造地球物理模型的选择和荷载形变计算方法见文献[13]，荷载形变与滤波结果(核点间隔为60天)如图10所示。从图中滤波曲线可以清晰看出，土壤湿度荷载影响最大，大气荷载影响次之，非潮汐海洋荷载影响再次之，积雪影响最小。滤波结果还显示出非构造影响的周年特征，与图4中KMIN站垂向周年运动特征相吻合。基于滤波结果的连续性非构造变化，对GPS垂向运动进行修正效果优于基于原始结果的修正，具体结果拟另文介绍。

3 结语

(1) 本文对多核函数时序滤波模型进行论述，并对多核函数滤波原理及其滤波特性进行探讨。认为当k=0.5、δ2=0.003时，单位权中误差较小；核点间隔控制滤波信息的频谱高低，间隔越大频谱信息越低，间隔越小频谱信息越高，实际应用中应根据不同研究对象合理选择核点间隔；另外，当观测数据出现缺失时，若连续丢失2个以上核点则滤波失败，若连续丢失2个核点则数据缺失部分滤波出现失真，若丢失1个核点则滤波效果不受影响，因此当数据连续性缺失时，应合理确定核点间隔，确保滤波信息可靠。

(2) 利用本文建立的多核函数滤波模型对GPS时序资料(包括站坐标序列、基线序列和应变序列)、定点形变资料(包括高差序列和基线序列)和非构造形变序列(包括非潮汐海洋、大气、土壤湿度和积雪荷载序列)进行滤波应用，获取其不同频谱的有效信息，验证其滤波效果的可靠性和稳定性，为对该资料的日常跟踪与分析提供了便利。

图10 KMIN站垂向非构造荷载形变序列Fig.10 The vertical non-tectonic deformation series of KMIN station

References)

[1] 李文静,杨国华,武艳强.地震前后唐山地震台地形变数据频谱特征分析[J].地震,2009,29(2)：141-146.

LI Wen-jing,YANG Guo-hua,WU Yan-qiang.Autoregression and Spectrum Analysis of Tangshan Deformation Data before and after Earthquakes[J].Earthquake,2009,29(2)：141-146.(in Chinese)

[2] 牛之俊,马宗晋,陈鑫连,等.中国地壳运动观测网络[J].大地测量与地球动力学,2002,22(3)：88-93.

NIU Zhi-jun,MA Zong-jin,CHEN Xin-lian,et al.Crustal Movement Observation Network of China[J].Journal of Geodesy and Geodynamics,2002,22(3)：88-93.(in Chinese)

[3] 李强,游新兆,杨少敏,等.中国大陆构造变形高精度大密度GPS监测——现今速度场[J].中国科学D辑：地球科学,2012,42(5)：629-632.

LI Qiang,YOU Xin-zhao,YANG Shao-min,et al.A Precise Velocity Field of Tectonic Deformation in China as Infered from Intensive GPS Observation[J].Sci China：Earth Sci,2012,42(5):629-632.(in Chinese)

[4] 张飞鹏，董大南，程宗颐.利用GPS监测中国地壳的垂向季节性变化[J].科学通报,2002,47(18)：1370-1379.

ZHANG Fei-peng,DONG Da-nan,CHENG Zong-yi.Seasonal Vertical Crustal Motions in China Detected by GPS[J].Chinese Science Bulletin,2002,47(21)：1772-1779.(in Chinese)

[5] 张恒璟，程鹏飞.基于经验模式分解的滤波去噪方法研究进展[J].测绘通报，2014(2)：1-4.

ZHANG Heng-jing,CHENG Peng-fei.Filtering and De-noising Based on Empirical Mode Decomposition Methods and Issues[J].Bulletin of Surveying and Mapping,2014(2)：1-4.(in Chinese)

[6] 武艳强,江在森,杨国华.最小二乘配置方法在提取GPS时间序列信息中的应用[J].国际地震动态，2007,343(7)：99-103.

WU Yan-qiang,JIANG Zai-sen,YANG Guo-hua.The Application of Least Square Collocation in Obtaining Information from GPS Time Series[J].Recent Developments in World Seismology,2007,343(7)：99-103.(in Chinese)

[7] 杨博,张风霜,韩月萍.多核函数法在GPS时序资料处理中的应用[J].大地测量与地球动力学,2010,30(4):137-141.

YANG Bo,ZHANG Feng-shuang,HAN Yue-ping.Method of Multi-kernel Function and Its Application in GPS Time Series Data Processing[J].Journal of Geodesy and Geodynamics,2010,30(4):137-141.(in Chinese)

[8] 黄立人,刘天奎.几种地壳垂直运动分析方法的数值实验[J].地壳形变与地震,1992,12(1):29-35.

HUANG Li-ren,LIU Tian-kui.A Numerical Experiment on Some Methods Used in Crustal Vertical Movement Analysis[J].Crustal Deformation and Earthquake,1992,12(1):29-35.(in Chinese)

[9] 陶本藻,王新洲,于正林,等.用于垂直形变模型的多面函数拟合法的试验研究[J].地壳形变与地震，1992,12(1):1-13.

TAO Ben-zao,WANG Xin-zhou,YU Zheng-lin,et al.An Investigation on Multiquadric Fitting Applied to Modeling of Vertical Crustal Movements[J].Crustal Deformation and Earthquake,1992,12(1):1-13.(in Chinese)

[10] 马洪滨,董仲宇.多面函数水准高程拟合中光滑因子求定方法[J].东北大学学报:自然科学版，2008,29(8)：1176-1178.

MA Hong-bin,DONG Zhong-yu.Calculation of Smooth Factor in GPS Level/ Elevation Fitting by Hardy Function[J].Journal of Northeastern University:Natural Science,2008,29(8)：1176-1178.(in Chinese)

[11] 黄继磊,罗志清.GPS高程拟合中一种多面函数平滑因子求法的研究[J].科学技术与工程,2013,13(6)：1157-1560.HUANG Ji-lei,LUO Zhi-qing.The Study on Method of Acquire Smooth Factor of Multi-quadric Function in GPS Height Fitting[J].Science Technology and Engineering,2013,13(6)：1157-1560.(in Chinese)

[12] 朱爽，周伟.甘肃岷县漳县6.6级地震前后区域地壳形变分析[J].地震工程学报，2015，37(3)：731-738.

ZHU Shuang,ZHOU Wei.Regional Crustal Deformation Analysis before and after the Minxin-Zhangxian Earthquake[J].China Earthquake Engineering Journal,2015,37(3):731-738.(in Chinese)

[13] 梁洪宝，刘志广，黄立人，等.非构造形变对中国大陆GNSS基准站垂向周期运动的影响[J].大地测量与地球动力学，2015，35(4)：46-50.

LIANG Hong-bao,LIU Zhi-guang,HUANG Li-ren,et al.Effects of Non-tectonic Deformation on the Vertical Periodic Motion of GNSS Reference Stations in China[J].Journal of Geodesy and Geodynamics，2015，35(4)：46-50.(in Chinese)

Study and Application of Multi-kernel Function Filtering Method in Time-series Deformation Data Processing

LIANG Hong-bao, LIU Xue-long, GUO Bing-hui

(FirstCrusta1MonitoringandApplicationCenter,CEA,Tianjin300180,China)

Due to observations of environmental impact carried out by monitoring stations on the China Mainland, we need to deal with the daily data for many years for daily tracking purposes. Existing filtering methods include the wavelet method, least-squares co-location method, Gaussian weighted average method, and Vondark method, etc., but for daily tracking these methods prove to be inconvenient. Aimed at solving the problem of extracting efficient information from sequential deformation observation data over some years, a filtering method, based on multi-kernel function, is studied in this paper. Taking continuous GPS station vertical sequential observation data as an example, we discuss the parameters for the multi-kernel function and their physical meaning. Conclusions are as follows: (1) When the kernel function index is 0.5 and the smoothing factor 0.003, the mean square error of unit weight of the filtering model with a kernel point interval of more than 10 days, is the least. (2) The kernel point interval controls the level of the filtering information frequency spectrum, the larger the interval, the lower the spectrum information; the smaller the interval, the higher the spectrum information. (3) Sometimes kernel points are lost because of missing data. When more than two continuous points are lost, the filtering fails; when two continuous points are lost, part of the filtering waves are distorted because of the missing data; when just one point is lost, the filtering effect is not affected. (4) From the filtering application in the GPS time-series data, the fixed-point deformation time-series data, and the non-tectonic deformation time-series data, information on different spectra are obtained and the stability and reliability of the method verified. This provides a more convenient way to daily process time-series observation data from a large number of stations.

multi-kernel function; time-series deformation data; kernel point interval; filtering

2015-04-23

中国地震局地震科技星火计划项目(XH15060Y)；科技部基础性工作专项(2015FY210400)；中国地震局青年震情跟踪课题(2015010211)

梁洪宝(1984-)，男，汉族，研究生，工程师，从事GNSS高精度解算与形变分析工作。E-mail：lhb131421@126.com。

P228.4

A

1000-0844(2016)05-0808-07

10.3969/j.issn.1000-0844.2016.05.808