張風(fēng)霜，占 偉

0 引言

在研究地殼形變問題時(shí)，常利用流動(dòng)GNSS復(fù)測資料計(jì)算指定區(qū)域的位移場、速度場和應(yīng)變場，如利用GNSS流動(dòng)資料研究中國大陸應(yīng)變場及塊體運(yùn)動(dòng)，取得了豐碩的成果 (江在森等，2003;王敏等，2003;楊國華等，2002;張培震等，2003，劉志廣等，2013)。但是流動(dòng)GNSS資料復(fù)測周期長，且其觀測資料的質(zhì)量、點(diǎn)位的穩(wěn)定性等各方面都與連續(xù)觀測資料存在較大的差距。隨著《中國地殼運(yùn)動(dòng)觀測網(wǎng)絡(luò)》二期工程的實(shí)施，GNSS連續(xù)觀測站由33個(gè)增至253個(gè)，由地方或其他部門建設(shè)的GNSS連續(xù)站也已超過千個(gè)，這使得利用GNSS連續(xù)觀測資料進(jìn)行形變場的分析成為可能。若發(fā)揮GNSS連續(xù)觀測資料時(shí)間上的優(yōu)勢，計(jì)算并分析月、季度、半年和年各個(gè)時(shí)段尺度下的形變場的動(dòng)態(tài)演變過程，可使GNSS連續(xù)觀測資料更有效地服務(wù)于震情會(huì)商和形變分析，因此利用GNSS連續(xù)觀測資料及時(shí)計(jì)算獲取各個(gè)時(shí)段尺度下穩(wěn)定可靠的速度場是非常必要的。GNSS連續(xù)觀測時(shí)間序列中除了線性成分外，還包含多種周期運(yùn)動(dòng)(江在森等，2006;王敏等，2005a;朱文耀等，2003)。白噪聲、閃爍噪聲和隨機(jī)漫步噪聲的噪聲性質(zhì)是GNSS連續(xù)觀測時(shí)間序列的基本特征 (黃立人，2006;黃立人，符養(yǎng)，2007;蔣志浩等，2010，夏峰等，2014)。如何用GNSS時(shí)間序列計(jì)算不同時(shí)段下的高精度速度場，獲取穩(wěn)定速度場與GNSS時(shí)間序列時(shí)段長度的關(guān)系是怎樣的?周期性非構(gòu)造運(yùn)動(dòng)、隨機(jī)漫步噪聲和閃爍噪聲等有色噪聲對(duì)不同時(shí)段長度的GNSS時(shí)間序列速度場的影響有多大?本文以上述問題為研究背景，采用小波分析方法和加權(quán)最小二乘法評(píng)定顧及周期性非構(gòu)造運(yùn)動(dòng)和不同噪聲影響的中國GNSS連續(xù)觀測網(wǎng)的月時(shí)段、半年時(shí)段、年時(shí)段的速度估值和精度，計(jì)算不同時(shí)間尺度下穩(wěn)定的GNSS連續(xù)觀測速度場。

1 GNSS數(shù)據(jù)處理

本文使用的GNSS數(shù)據(jù)是中國地殼運(yùn)動(dòng)觀測網(wǎng)絡(luò)253個(gè)基準(zhǔn)站運(yùn)行以來記錄的連續(xù)觀測資料，其中一期網(wǎng)絡(luò)基準(zhǔn)站33個(gè)，觀測時(shí)間為1999年初至2013年4月的數(shù)據(jù)，二期網(wǎng)絡(luò)基準(zhǔn)站220個(gè)，觀測時(shí)間為2010年至2013年4月的數(shù)據(jù)，參與解算的還有90個(gè)IGS基準(zhǔn)站在相應(yīng)時(shí)間段內(nèi)的觀測數(shù)據(jù)。GNSS連續(xù)站觀測值的數(shù)據(jù)處理采用GAMIT/QOCA軟件完成 (Herring et al.，2010)。數(shù)據(jù)處理的基本流程 (王敏等，2005b;王敏，2007)首先利用GAMIT獲得網(wǎng)絡(luò)工程253個(gè)連續(xù)站及90個(gè)IGS測站的單日松弛解，同步觀測的測站較多時(shí)，采用分區(qū)處理。完成GAMIT計(jì)算之后，利用QOCA軟件將計(jì)算所得的各區(qū)單日松弛解進(jìn)行綜合平差，在此基礎(chǔ)上通過IGS核心站求解相對(duì)于全球參考框架ITRF2005的相似變換7參數(shù)，從而獲得ITRF2005下的單日解，即GLOBK的NEU坐標(biāo)值，不同于通常的站心坐標(biāo)，其具體含義可參閱黃立人等 (2006)。對(duì)計(jì)算得到的各測站坐標(biāo)時(shí)間序列，在剔除粗差的基礎(chǔ)上，我們進(jìn)行了2001年昆侖山口西8.1級(jí)地震、2004年蘇門答臘9級(jí)地震、2008年汶川8.0級(jí)地震和2013年廬山7.0級(jí)地震同震位移改正的處理。

2 周期性運(yùn)動(dòng)對(duì)不同時(shí)段長度GNSS時(shí)間序列速度場的影響

通過剔除GNSS站坐標(biāo)時(shí)間序列的粗差，僅可以削弱觀測值中的偶然成分，還需要對(duì)信息進(jìn)一步分離。GNSS連續(xù)觀測時(shí)間序列中包含著多種有規(guī)律的、周期性成分，研究表明，在GNSS連續(xù)觀測時(shí)間序列中比較顯著且共性較好的周期成分是年周期和半年周期 (江在森等，2006)?？刹捎弥C波擬合法(王敏等，2005b;韓鵬等，2009;謝凡等，2011)剔除GNSS連續(xù)觀測時(shí)間序列中的年周期和半年周期成分。諧波擬合方法是在傅里葉級(jí)數(shù)及傅里葉變換的基礎(chǔ)上，把滿足一定條件的周期函數(shù)分解為組成它的各次諧波成分。對(duì)GNSS連續(xù)觀測時(shí)間序列去線性趨勢后得到的殘差時(shí)間序列的擬合函數(shù)f(t)通過傅里葉級(jí)數(shù)表示成如下形式:

式中，f(t)可以看成是這些正弦函數(shù)和余弦函數(shù)的“加權(quán)和”，其中m為諧波的次數(shù)，系數(shù)Am、Bm是權(quán)重因子，表示各次諧波對(duì)總序列的貢獻(xiàn)，可以通過最小二乘法計(jì)算得到。對(duì)GNSS連續(xù)觀測時(shí)間序列而言，m取2，即用2次諧波擬合就可以較好的改正GNSS時(shí)間序列中的年周期和半年周期成分。

對(duì)GNSS連續(xù)觀測時(shí)間序列中的年周期和半年周期成分用諧波擬合的方法進(jìn)行改正后，筆者采用加權(quán)最小二乘法計(jì)算并評(píng)定未剔除和剔除周期性運(yùn)動(dòng)后的測站年速度及精度。最小二乘準(zhǔn)則為(武漢大學(xué)測繪學(xué)院測量平差學(xué)科組，2009)

剔除年、半年周期后計(jì)算得到半年速度的精度和年速度一樣也有較大幅的提高。表2給出了2003～2012年 (共10年)哈爾濱站和泰安站以每半年為一個(gè)時(shí)段計(jì)算所得的東向速度對(duì)比。由表中可見，2003～2012年哈爾濱站每年上半年未剔除周期性運(yùn)動(dòng)的速度明顯小于下半年速度，而剔除周期性運(yùn)動(dòng)后各年份半年速度間的差異明顯減小。哈爾濱站和泰安站剔除年周期和半年周期后得到的半年時(shí)段速度的離散度相對(duì)未剔除的結(jié)果顯著減小，說明剔除年周期和半年周期后的半年時(shí)段速度更穩(wěn)定，可靠性更高。

剔除GNSS連續(xù)觀測時(shí)間序列中的年、半年周期，并計(jì)算相對(duì)于區(qū)域無旋轉(zhuǎn)參考基準(zhǔn) (楊國華等，2005;胡新康等，2007)的速度，筆者計(jì)算得到了2010～2013年中國大陸地殼運(yùn)動(dòng)觀測網(wǎng)絡(luò)工程GNSS連續(xù)觀測站的背景速度場 (圖1中藍(lán)色箭頭)，其中一期網(wǎng)絡(luò)33個(gè)、二期網(wǎng)絡(luò)220個(gè)GNSS連續(xù)站計(jì)算所采用時(shí)間序列的長度分別為1999～2013年和2010～2013年。為了反映GNSS連續(xù)觀測站速度場的年度差異運(yùn)動(dòng)，筆者利用2011年4月至2012年4月的GNSS連續(xù)觀測資料，計(jì)算了2011～2012年中國大陸GNSS連續(xù)觀測站速度場 (圖1中黃色箭頭)。由圖1可以發(fā)現(xiàn)，計(jì)算得到的GNSS連續(xù)觀測站背景速度場與我國地質(zhì)構(gòu)造格局較為相關(guān)，基本反映了中國大陸的形變背景。2011年日本9.0級(jí)地震后，2011～2012年速度場相對(duì)背景速度場的差異變化在一定程度上說明了GNSS連續(xù)觀測時(shí)間序列年速度場能夠顯示中國大陸地殼水平差異運(yùn)動(dòng)的空間分布信息。

由上述分析，剔除GNSS時(shí)間序列中的年周期和半年周期成分，可以得到穩(wěn)定且高精度的年時(shí)段和半年時(shí)段速度場。但是，用上述方法計(jì)算了2011年4月初至2013年3月底 (共24個(gè)月)GNSS連續(xù)觀測站各月時(shí)段的速度，卻發(fā)現(xiàn)各個(gè)月時(shí)段速度離散度大，很不穩(wěn)定，不同年份相對(duì)應(yīng)月份的速度差別也很大，因此要得到穩(wěn)定的月時(shí)段速度，還需對(duì)GNSS連續(xù)觀測時(shí)間序列進(jìn)行更深入的信息分離。

表1 未剔除和剔除周期后計(jì)算得到的GNSS測站年速度精度對(duì)比Tab.1 Precision comparison of the calculated annual velocity of GNSS stations after rejecting the cycle motion or not

表2 哈爾濱和泰安站半年時(shí)段東向速度離散度比較Tab.2 Comparison of the dispersion of velocity in east recorded by Harbin and Taian Stations in the period of half a year

?

3 小波分析剔除不同噪聲影響的速度估值

為進(jìn)一步分離GNSS站坐標(biāo)時(shí)間序列所包含的有色噪聲，我們采用小波變換的方法。小波分析可以提取時(shí)間序列的時(shí)頻特征，具有分析信號(hào)的局部特征和多尺度分辨率的特點(diǎn)，借助小波方法可以將不同特性的噪聲進(jìn)行分離與估計(jì) (飛思科技產(chǎn)品研發(fā)中心，2005)。

設(shè)f(t)是觀測點(diǎn)變形序列，采用Mallat算法，將信號(hào)分解成不同頻率成分:

其中，信號(hào)在空間Vj上的投影為

在空間Wj上的投影為

其中，Aj－1f(x)和Ajf(x)分別是信號(hào)f(x)的頻率不超過2－j+1和2－j的成分，而 Djf(x)是頻率介于 2－j與2－j+1之間的成分。

上述小波分解式實(shí)際可寫成如下的矩陣形式:

式中:H是尺度函數(shù)對(duì)應(yīng)的低通濾波器，H=(hk－2n);G是小波函數(shù)對(duì)應(yīng)的帶通濾波器，G=(gk－2n);Cj是2j在分辨率下的離散逼近;Dj是2j在分辨率下的離散細(xì)節(jié);{hk}k∈Z和 {gk}k∈Z是一對(duì)離散正交鏡像濾波器 (低通和高通)。式 (6)即為Mallat塔式分解算法。

小波變換要求數(shù)據(jù)是等間隔采樣的，而GNSS時(shí)間序列的日觀測值中有一定的數(shù)據(jù)缺失，我們采用最小二乘配置方法對(duì)GNSS時(shí)間序列進(jìn)行插值。由于最小二乘配置方法不像譜分析或是小波變換要求數(shù)據(jù)等間隔，其數(shù)據(jù)輸入可以較為松散，而其輸出可以根據(jù)需求決定是否等間隔 (武艷強(qiáng)等，2007)。這就決定了該方法可以作為譜分析、小波變換等方法的預(yù)處理過程。即先由最小二乘配置輸出濾波、插值后的等間隔值，然后再由其它方法進(jìn)行后繼處理。

武漢大學(xué)測繪學(xué)院測量平差學(xué)科組 (2009)詳細(xì)論述了最小二乘配置，其誤差方程為

式中，BZ為信號(hào)的系數(shù)陣 (通常觀測點(diǎn)部分為單位陣，推估點(diǎn)部分為零陣)，Z^為隨機(jī)信號(hào)估值(包括觀測點(diǎn)和推估點(diǎn))，G為經(jīng)典平差問題的系數(shù)陣，Y^為經(jīng)典平差問題需要求解的待定參數(shù)，L為觀測值，LZ為信號(hào)的虛擬觀測值 (一般為信號(hào)的數(shù)學(xué)期望)。

根據(jù)間接平差理論和矩陣反演理論對(duì)式 (7)進(jìn)行求解得到

式中，DX為已測點(diǎn)信號(hào)的協(xié)方差陣，DΔ為觀測值的協(xié)方差陣。

應(yīng)用最小二乘配置方法的前提是獲取信號(hào)的協(xié)方差陣，推薦使用如式 (9)的高斯經(jīng)驗(yàn)協(xié)方差函數(shù):

式中，t為時(shí)間間隔，A和k為待定參數(shù)。采用赫爾默特方差估計(jì)方法獲得參數(shù)A，而對(duì)于參數(shù)k，可以根據(jù)關(guān)注的信息頻段人為設(shè)定。如果k值較大，得到的信息就接近于原始序列;如果k值較小，得到的信息就是相對(duì)低頻的數(shù)據(jù)。

對(duì)最小二乘配置插值后的數(shù)據(jù)采用具有對(duì)稱性和正則性的sym8小波進(jìn)行小波分解 (張風(fēng)霜，楊國華，2011)。以2011～2013年哈爾濱站東向站坐標(biāo)時(shí)間序列小波分解的結(jié)果為例 (圖2)，隨著階數(shù)的增加，信號(hào)的細(xì)節(jié)部分從最高頻向低頻分解，當(dāng)階數(shù)為0時(shí)，信號(hào)為采樣頻率;當(dāng)階數(shù)為1時(shí)，頻率二等分，對(duì)于GNSS連續(xù)觀測時(shí)間序列，1階信號(hào)的周期為1～2 d，同樣，2階信號(hào)的周期為2～4 d，以此類推，當(dāng)階數(shù)為8時(shí)，8階信號(hào)的周期為128～256 d。研究表明，GNSS時(shí)間序列小波分解的1～3階的高頻信號(hào)主要是白噪聲，4、5階的高頻信號(hào)主要是閃爍噪聲，6、7階的高頻信號(hào)主要是隨機(jī)漫步噪聲 (楊國華等，2007;楊博等，2010)。

對(duì)哈爾濱站和泰安站已剔除年周期和半年周期的東向站坐標(biāo)時(shí)間序列，分別剔除白噪聲、閃爍噪聲和隨機(jī)漫步噪聲后再計(jì)算2011年4月初至2013年3月底 (共24個(gè)月)各月時(shí)段速度，結(jié)果見表3、4和圖3。由圖3可見，隨著白噪聲、閃爍噪聲和隨機(jī)漫步噪聲的逐步剔除，24個(gè)月各月的速度由離散逐步趨于穩(wěn)定。由表3、4可知，只剔除白噪聲，計(jì)算得到的哈爾濱站和泰安站24個(gè)月的各月速度離散度較大，分別為0.795和0.464，說明各月速度并不穩(wěn)定，進(jìn)一步剔除閃爍噪聲后，哈爾濱站和泰安站的月速度離散度分別減小至0.307和0.247，再進(jìn)一步剔除隨機(jī)漫步噪聲后，哈爾濱站和泰安站離散度分別進(jìn)一步減小至0.165和0.180，說明剔除白噪聲、閃爍噪聲、隨機(jī)漫步噪聲后得到的月時(shí)段速度較穩(wěn)定。

表3 哈爾濱站東向剔除不同頻段噪聲后月時(shí)段速度對(duì)比Tab.3 Comparison of monthly velocity in east after rejecting different noises recorded by Harbin Station

表4 泰安站東向剔除不同頻段噪聲后月時(shí)段速度對(duì)比Tab.4 Comparison of monthly period velocity of rejecting different noises of the east recorded by Taian Station

4 認(rèn)識(shí)與討論

通過剔除GNSS連續(xù)觀測時(shí)間序列的年周期和半年周期成分，可得到穩(wěn)定、高精度的年時(shí)段和半年時(shí)段速度場，但由于GNSS連續(xù)觀測時(shí)間序列的噪聲性質(zhì)，短時(shí)段穩(wěn)定速度場的獲取需要對(duì)GNSS連續(xù)觀測時(shí)間序列進(jìn)行濾波剔除有色噪聲。本文采用小波分析方法分別剔除了GNSS連續(xù)觀測時(shí)間序列中的白噪聲、閃爍噪聲和隨機(jī)漫步噪聲，并比較了逐步剔除不同頻段噪聲后計(jì)算得到的各月時(shí)段速度的精度和離散度，結(jié)果表明，剔除白噪聲、閃爍噪聲、隨機(jī)漫步噪聲后得到的月時(shí)段速度更穩(wěn)定可靠。利用本文的方法，實(shí)時(shí)計(jì)算并分析月、半年和年各個(gè)時(shí)段尺度下的形變場的動(dòng)態(tài)演變過程，有利于發(fā)揮GNSS連續(xù)觀測資料時(shí)間上的優(yōu)勢，從而使GNSS連續(xù)觀測資料更有效地服務(wù)于震情會(huì)商和形變分析。

飛思科技產(chǎn)品研發(fā)中心.2005.小波分析理論與實(shí)現(xiàn)［M］．北京:電子工業(yè)出版社.

韓鵬，黃清華，修濟(jì)剛.2009.地磁日變與地震活動(dòng)關(guān)系的主成分分析［J］．地球物理學(xué)報(bào)，52(6):1556－1563.

胡新康，王倩，馬青，等.2007.區(qū)域無整體旋轉(zhuǎn)基準(zhǔn)的研究與應(yīng)用［J］．大地測量與地球動(dòng)力學(xué)，27(2):52－60.

黃立人，符養(yǎng).2007.GPS連續(xù)觀測站的噪聲分析［J］．地震學(xué)報(bào)，29(2):197－202.

黃立人，高硯龍，任立生，等.2006.關(guān)于NEU(ENU)坐標(biāo)系統(tǒng)［J］．大地測量與地球動(dòng)力學(xué)，26(1):97－99.

黃立人.2006.GPS基準(zhǔn)站時(shí)間序列的噪聲特性分析［J］．大地測量與地球動(dòng)力學(xué)，26(2):31－33.

江在森，鄧志輝，丁鑒海，等.2006.GPS、衛(wèi)星遙感及地球變化磁場地震短期預(yù)測方法研究［M］．北京:地震出版社.

江在森，馬宗晉，張希，等.2003.GPS初步結(jié)果揭示的中國大陸水平應(yīng)變場與構(gòu)造變形［J］．地球物理學(xué)報(bào)，46(3):352－358.

蔣志浩，張鵬，秘金鐘，等.2010.顧及有色噪聲影響的CGCS2000下我國CORS站速度估計(jì)［J］．測繪學(xué)報(bào)，39(4):355－363.

劉志廣，楊博，盧雙印，等.2013.青藏高原中南部近期地殼水平形變［J］.大地測量與地球動(dòng)力學(xué)，33(3):16－20.

王敏，沈振康，牛之俊，等.2003.現(xiàn)今中國大陸的地殼運(yùn)動(dòng)與活動(dòng)塊型［J］．中國科學(xué)(D 輯)，33(增刊):21－32.

王敏，沈正康，董大南.2005a.非構(gòu)造形變對(duì)GPS連續(xù)站位置時(shí)間序列的影響和修正［J］．地球物理學(xué)報(bào)，48(5):1045－1052.

王敏，張祖勝，許明元.2005b.2000國家GPS大地控制網(wǎng)的數(shù)據(jù)處理和精度評(píng)估［J］．地球物理學(xué)報(bào)，48(4):817－823.

王敏.2007.GPS數(shù)據(jù)處理方面的最新進(jìn)展及其對(duì)定位結(jié)果的影響［J］．國際地震動(dòng)態(tài)，(7):3－8.

武漢大學(xué)測繪學(xué)院測量平差學(xué)科組.2009.誤差理論與測量平差基礎(chǔ)［M］．武漢:武漢大學(xué)出版社.

武艷強(qiáng)，江在森，楊國華.2007.最小二乘配置方法在提取GPS時(shí)間序列信息中的應(yīng)用［J］．國際地震動(dòng)態(tài)，(7):99－103.

夏峰，張銳，馮勝濤，等.2014.GNSS不同類型實(shí)驗(yàn)觀測墩對(duì)不同類型噪聲幅值影響［J］.華北地震科學(xué)，32(2):39－44.

謝凡，滕云田，胡星星.2011.數(shù)學(xué)形態(tài)濾波在地磁數(shù)據(jù)干擾抑制中的應(yīng)用［J］．地球物理學(xué)進(jìn)展，26(1):147－156.

楊博，張風(fēng)霜，韓月萍.2010.GPS連續(xù)站水平分量時(shí)間序列共模誤差識(shí)別的歐拉—小波法［J］．大地測量與地球動(dòng)力學(xué)，30(3):100－104.

楊國華，江在森，武艷強(qiáng)，等.2005.中國大陸整體無凈旋轉(zhuǎn)基準(zhǔn)及其應(yīng)用［J］．大地測量與地球動(dòng)力學(xué)，25(4):6－10.

楊國華，李延興，韓月萍，等.2002.由GPS觀測結(jié)果推導(dǎo)中國大陸現(xiàn)今水平應(yīng)變場［J］．地震學(xué)報(bào)，24(4):337－347.

楊國華，張風(fēng)霜，武艷強(qiáng)，等.2007.GPS基準(zhǔn)站坐標(biāo)分量噪聲的時(shí)間序列與分類特征［J］．國際地震動(dòng)態(tài)，(7):80－85.

張風(fēng)霜，楊國華.2011.GPS連續(xù)觀測資料小波變換質(zhì)量評(píng)價(jià)探討［J］．測繪科學(xué)，36(6):10 －12.

張培震，鄧啟東，張國民，等.2003.中國大陸的強(qiáng)震活動(dòng)與活動(dòng)地塊［J］．中國科學(xué)(D 輯)，33(增刊):12－20.

朱文耀，符養(yǎng)，李彥.2003.GPS高程導(dǎo)出的全球高程振蕩運(yùn)動(dòng)及季節(jié)變化［J］．中國科學(xué)，33(增刊):72－83.

Herring T.A.，King R.W.，McClusky S.C..GAMIT Reference Manual:GPS Analysis at MIT.Release 10.4［EB/OL］．(2010 －01 －01)［2014－07－19］．http://www.doc88.com/p－31570247490.html.