時序形變資料的多核函數(shù)濾波方法研究及應用

梁洪寶, 劉雪龍， 郭炳輝

(中國地震局第一監(jiān)測中心，天津 300180)

?

時序形變資料的多核函數(shù)濾波方法研究及應用

梁洪寶, 劉雪龍， 郭炳輝

(中國地震局第一監(jiān)測中心，天津 300180)

針對多年時序形變觀測資料有效信息提取復雜的問題，對基于多核函數(shù)的濾波方法進行研究，得到以下有益結論：(1)當核函數(shù)指數(shù)為0.5，光滑因子為0.003時，10天及以上核點間隔的濾波模型單位權中誤差最小；(2)核點間隔控制濾波信息頻譜的高低, 間隔越大頻譜信息越低，反之則頻譜信息越高；(3)因數(shù)據缺失部分造成核點減少，當連續(xù)減少2個以上時濾波失敗，當連續(xù)減少2個時數(shù)據缺失部分濾波出現(xiàn)失真，當減少1個時濾波效果不受影響；(4)通過對GPS時序資料、定點形變時序資料和非構造形變時序資料的濾波應用，獲取不同頻譜的信息，驗證了本文方法的穩(wěn)定性和可靠性。

多核函數(shù)； 時序形變資料； 核點間隔； 濾波

0 引言

為監(jiān)測中國大陸斷層的構造活動，自1976年唐山7.8級地震以來，在重點斷層帶建立了定點形變地震觀測臺站，如唐山地震臺至今積累了三十多年的單天形變觀測資料[1]；上世紀90年代，隨著GPS技術興起并成功應用于地殼形變監(jiān)測中，我國啟動了國家重大科學工程——中國地殼運動觀測網絡[2](簡稱“網絡工程”)，建立了25個GPS基準站(后來陸續(xù)增加到30個)，并于1998年開始運行；2007年啟動的國家重大科技基礎設施——中國大陸構造環(huán)境監(jiān)測網絡[3](簡稱“陸態(tài)網絡”)，對網絡工程GPS基準站進行加密，共建成260個測站，并于2010年下半年開始運行，至今積累了數(shù)年的單天觀測資料。眾所周知，上述形變觀測資料均含有大氣、非潮汐海洋、積雪和土壤濕度等非構造形變信息，對于地殼形變分析而言，需要對其含有的非構造形變信息進行剔除[4]。根據多年非構造地球物理模型可以獲取目標區(qū)域的非構造形變序列，以此對觀測資料進行修正。

時序形變原始觀測資料不能直接用于形變分析，需對其含有的有效信息進行提取，即根據研究需要對包含多種信息的原始資料進行濾波，以得到較為客觀、清晰的信息。時序資料的濾波方法有多種[5]，但能完全滿足需求的方法并不多，如樣條函數(shù)方法和最小二乘配置法從理論上講均可以滿足要求，但在實際應用中卻存在計算不穩(wěn)定、耗時較長等問題[6]。利用高斯型函數(shù)作為多核函數(shù)的核函數(shù)對GPS時間序列進行濾波，其計算過程較為簡易[7]。 因此，本文在此基礎上探討基于指數(shù)型核函數(shù)的多核函數(shù)模型對時間跨度較長的時序形變資料的濾波效果，確定其較優(yōu)的光滑因子，并將其應用于上述三種時序形變資料的濾波處理中。

1 多核函數(shù)模型及濾波原理

多核函數(shù)法是由美國Hardy教授于1977年提出并建議用于地殼垂直形變分析的一種純數(shù)學的方法，其基本思想是任何一個圓滑的數(shù)學表面總可以用一系列有規(guī)則的數(shù)學表面的總和以任意精度逼近[8-9]。其數(shù)學表達式為：

(1)

式中:Q為核函數(shù);a為待定系數(shù);m為核函數(shù)個數(shù)。核函數(shù)可任意選用，一般表示為：

(2)

式中:δ2稱為光滑因子，用來對核函數(shù)進行調整;k為非零實數(shù)，一般取0.5、-0.5或1.5。

根據多核函數(shù)法的基本思想，假設時序形變的變化是由多條曲線組成的一條圓滑曲線，因此可以將其應用于以時間變化的一維自變量的領域，函數(shù)表達式可以描述為：

(3)

相應的核函數(shù)可以描述為：

(4)

1.1 光滑因子的確定

以長春IGS站2005-2011年垂向去線性趨勢運動序列為例(以下類同)，核點間隔分別為10、20、30、60、90和120天，光滑因子在0～0.5取值，步長取值為0.000 5，滑動計算1 000次，并以驗后單位權中誤差為評價指標，探討不同核點間隔下光滑因子與單位權中誤差間的關系，結果如圖1所示。從圖1可以看出，在理想的單位權中誤差下，不同的核點間隔光滑因子的取值均有一個相對穩(wěn)定的區(qū)間，這個區(qū)間隨核點間隔的增大而增大。核點間隔為10天時，光滑因子穩(wěn)定區(qū)間為[0.001,0.005],為了應用簡便快捷，本文取中間值0.003，以滿足核點間隔為10天及以上所有情況。

1.2 核點間隔與濾波信息的關系

確定光滑因子為0.003后，分別取核點間隔為20、30、60、90、120和150天，探討不同核點間隔與濾波信息的關系，結果如圖2所示。從圖2可以看出，核點間隔大小控制濾波信息的頻譜，間隔越大頻譜信息越低，間隔越小頻譜信息越高，因此在實際應用中要根據研究對象合理確定核點間隔，以獲取所需信息。

1.3 數(shù)據連續(xù)性缺失對濾波效果的影響

因各種因素影響，在實際監(jiān)測工作中觀測資料難免出現(xiàn)連續(xù)性缺失。為穩(wěn)健起見，將長春站2006年10—12月的連續(xù)觀測數(shù)據刪除，以1.1節(jié)確定的光滑因子為0.003，核點間隔分別為20、30、39、46、60和70天，探討數(shù)據連續(xù)性缺失對濾波效果的影響，結果如圖3所示。與圖2對比，從圖3可以看出，核點間隔為20天和30天時數(shù)據缺失部分濾波出現(xiàn)錯誤；核點間隔為39天和46天時數(shù)據缺失部分濾波出現(xiàn)失真；剩余兩種情況濾波正常。經過分析，圖3(a)和(b)中數(shù)據缺失的同時也造成核點分別減少了4個和3個，圖3(c)和(d)中核點丟失了2個，圖3(e)和(f)中核點丟失了1個，由此可以認為，核點連續(xù)丟失數(shù)不能超過2個。當觀測數(shù)據連續(xù)性缺失時若連續(xù)丟失2個以上核點則濾波失敗；若連續(xù)丟失2個核點則數(shù)據缺失部分濾波出現(xiàn)失真；若丟失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 數(shù)據缺失時不同核點間隔的濾波效果Fig.3 The filtering effect of different kernel point intervals in the case of when missing data missing

2 時序形變資料的濾波應用

根據上述濾波方法及參數(shù)設置，分別對GPS時序資料、定點形變時序資料和非構造形變時序資料進行濾波應用。

2.1 GPS時序資料

隨著網絡工程和陸態(tài)網絡項目的實施，連續(xù)GPS站產出了多年時序資料，如圖4(黑色圓點表示計算值，紅色曲線表示本文方法的濾波值，下圖表示方法相同)所示，測站分布如圖5所示。KMIN(昆明)站積累了16年左右的坐標序列，從圖4(a)中濾波曲線(核點間隔為90天)可以看出，該站水平方向主要以線性運動為主，垂向以周年運動為主。為探測該站的精細運動，對其線性運動進行剔除，結果如圖4(b)所示。通過濾波曲線可以看出，該站自2004年開始南向運動加劇，積累的大量能量經2004年12月26日蘇門答臘9.3級地震影響得到一定程度的釋放；自2007年開始南北向運動出現(xiàn)“閉鎖”，由此積累的能量經2008年5月12日汶川8.0級地震的影響再次進行釋放，構造運動經過后續(xù)四年左右時間的調整期逐漸恢復到原來的運動特征。其東西向也因受汶川地震影響出現(xiàn)西向運動加劇的特征，一直持續(xù)至今，垂向運動受地震影響較小；除構造運動外，濾波曲線還顯示出年周期的非構造運動，垂向最為顯著。

圖4 KMIN站坐標序列Fig.4 Coordinate series of KMIN station

圖5 地震活動與GPS測站分布Fig.5 Seismic activity and the GPS station distribution

基線序列是用來分析塊體間構造運動的拉張或擠壓特征，應變序列是用來分析區(qū)域塊體的受力應變特征。基于GPS連續(xù)站序列可以計算測站基線和塊體應變序列[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(西寧)三站得到該區(qū)域的應變系列，該區(qū)域地處青藏高原東北緣，是印度與歐亞兩大板塊碰撞作用由近南北方向向北東、東方向轉換的重要場所，也是中國大陸東西及南北構造結合部位和重要的構造轉換區(qū)域。圖7中其東西向和南北向均處于持續(xù)擠壓狀態(tài)特征也說明了這一點。2001年11月14日昆侖山8.1級地震和2010年4月14日7.1級地震對該區(qū)域均有一定影響，分別表現(xiàn)出東西方向和南北方向的拉張變形特征。

圖6 GPS基線序列Fig.6 GPS baseline series

圖7 GPS應變序列Fig.7 GPS strain series

2.2 定點形變時序資料濾波應用

唐山地震臺始建于1978年，位于唐山市路南區(qū)復興路原第十中學院內，是唐山地區(qū)唯一的大地形變臺站。1976年唐山7.8級地震后，地殼升降和錯動的巨大形變大體上都以發(fā)震斷層為交界面，由于深部形變與強烈地震波的共同作用，在震區(qū)形成了大量的總體走向為N30°～40° E呈雁狀排列的地裂縫，其中發(fā)震斷層附近尤為明顯。唐山地震臺正處于地裂縫最為發(fā)育的典型地段，臺站布設有水準和基線測線各4條，長度24～48 m不等，以不同角度與斷層交匯。測站分布如圖8所示，基線測量順序是J1、J2、J3、J4、J1，水準測量順序是S1、S2、S3、S4、S1。水準S1-S2和基線J3-J4觀測結果如圖9所示，從圖中可以看出，通過不同核點間隔(分別為30天、182天和365天)可以獲取不同頻譜的信息，除趨勢性變化外，還可以清晰展現(xiàn)短期動態(tài)變化。在高差曲線變化中曲線上升表示斷裂正斷活動，下降則表示逆斷活動；在基線曲線變化中，上升則表示斷裂拉張運動，下降則表示擠壓運動。

圖8 唐山臺觀測點分布Fig.8 Distribution of observation points at Tangshan station

圖9 斷層高差與基線序列Fig.9 The fault height difference and baseline series

2.3 非構造形變時序資料濾波應用

由于垂直形變觀測資料是多種信息的綜合體現(xiàn)，且較水平運動弱。在進行地殼垂直形變分析時，必須對其含有的非構造形變信息進行剔除，這種非構造形變主要是由非潮汐海洋、大氣、積雪和土壤濕度等非構造地球物理因素產生的。本文計算了KMIN(昆明)GPS站的非構造垂向形變，非構造地球物理模型的選擇和荷載形變計算方法見文獻[13]，荷載形變與濾波結果(核點間隔為60天)如圖10所示。從圖中濾波曲線可以清晰看出，土壤濕度荷載影響最大，大氣荷載影響次之，非潮汐海洋荷載影響再次之，積雪影響最小。濾波結果還顯示出非構造影響的周年特征，與圖4中KMIN站垂向周年運動特征相吻合。基于濾波結果的連續(xù)性非構造變化，對GPS垂向運動進行修正效果優(yōu)于基于原始結果的修正，具體結果擬另文介紹。

3 結語

(1) 本文對多核函數(shù)時序濾波模型進行論述，并對多核函數(shù)濾波原理及其濾波特性進行探討。認為當k=0.5、δ2=0.003時，單位權中誤差較小；核點間隔控制濾波信息的頻譜高低，間隔越大頻譜信息越低，間隔越小頻譜信息越高，實際應用中應根據不同研究對象合理選擇核點間隔；另外，當觀測數(shù)據出現(xiàn)缺失時，若連續(xù)丟失2個以上核點則濾波失敗，若連續(xù)丟失2個核點則數(shù)據缺失部分濾波出現(xiàn)失真，若丟失1個核點則濾波效果不受影響，因此當數(shù)據連續(xù)性缺失時，應合理確定核點間隔，確保濾波信息可靠。

(2) 利用本文建立的多核函數(shù)濾波模型對GPS時序資料(包括站坐標序列、基線序列和應變序列)、定點形變資料(包括高差序列和基線序列)和非構造形變序列(包括非潮汐海洋、大氣、土壤濕度和積雪荷載序列)進行濾波應用，獲取其不同頻譜的有效信息，驗證其濾波效果的可靠性和穩(wěn)定性，為對該資料的日常跟蹤與分析提供了便利。

圖10 KMIN站垂向非構造荷載形變序列Fig.10 The vertical non-tectonic deformation series of KMIN station

References)

[1] 李文靜,楊國華,武艷強.地震前后唐山地震臺地形變數(shù)據頻譜特征分析[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監(jiān)測——現(xiàn)今速度場[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監(jiān)測中國地殼的垂向季節(jié)性變化[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].國際地震動態(tài)，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] 楊博,張風霜,韓月萍.多核函數(shù)法在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] 黃立人,劉天奎.幾種地殼垂直運動分析方法的數(shù)值實驗[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] 陶本藻,王新洲,于正林,等.用于垂直形變模型的多面函數(shù)擬合法的試驗研究[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] 馬洪濱,董仲宇.多面函數(shù)水準高程擬合中光滑因子求定方法[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高程擬合中一種多面函數(shù)平滑因子求法的研究[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級地震前后區(qū)域地殼形變分析[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