JLCORS網絡參考站坐標時間序列噪聲類型分析

劉俊清，肖輝峰，張晨俠，丁　廣

(1. 吉林大學地球探測科學與技術學院，吉林 長春 130026； 2. 吉林省地震局，吉林 長春 130117； 3. 北京麥格天寶科技股份有限公司，北京 100029)

?

JLCORS網絡參考站坐標時間序列噪聲類型分析

劉俊清1,2，肖輝峰3，張晨俠2，丁廣2

(1. 吉林大學地球探測科學與技術學院，吉林 長春 130026； 2. 吉林省地震局，吉林 長春 130117； 3. 北京麥格天寶科技股份有限公司，北京 100029)

利用吉林省連續運行參考站網絡(JLCORS)近3年的觀測資料，使用GAMIT/GLOBK軟件進行后處理，解算時用吉林省周邊IGS及陸態網絡站點作為控制點，獲得41個站單日松弛解的坐標時間序列。用最大似然法擬合經典的噪聲模型，根據噪聲模型是否能夠取得似然函數的最大值，獲得不同噪聲模型及其組合的參數估計。結果表明，JLCORS全部參考站坐標時間序列噪聲含有冪律噪聲及白噪聲，可用二者的組合噪聲模型來描述，不含帶通濾波噪聲及各階高斯-馬爾科夫噪聲。冪律噪聲模型部分的譜指數集中分布在-2附近，表明冪律噪聲的特例隨機游走噪聲占有主導地位。

CORS；IGS；陸態網；噪聲模型

地球動力學是建立在觀測試驗基礎上的一門學科，觀測數據不可避免地存在誤差、不穩定和不確定等非地殼構造因素的影響[1-3]，誤差性質及分離方法的研究一直是地學工作者的一項重要工作。由于地殼構造活動是一個非常緩慢的過程，地殼監測時間序列跨度很長，通常是幾年甚至幾十年，這就意味著觀測值中含有大量的與時間相關的誤差，而且誤差源也可能隨著時間變化[4-6]。如在多年的觀測過程中，儀器老化、更換及觀測墩標不穩定等，均與時間具有相關性，都會在觀測數據時間序列中產生誤差。隨著研究的深入，逐漸發現這類與時間相關的誤差的一些統計規律，可用白噪聲(white noise，WN)和冪率噪聲(power law noise，PLN)的組合模型進行很好的模擬[4]。白噪聲可以通過增加觀測次數來降低影響，冪率噪聲卻無法消除,只能通過噪聲相對準確的統計關系，在數據處理時建立數學模型來消除。

目前，很多國家建立了國家級、區域級的CORS網絡，主要用于測繪工作中。建站標準及其選址要求都低于地殼運動研究的標準，大部分站使用土層標墩建設規范，甚至建在了建筑物頂部。這樣的觀測環境會給GNSS觀測帶來很多誤差信號，對不同用途的觀測要求，需要相應地選擇計算模型消除誤差，這就需要對誤差性質進行研究，獲得其統計表達式，然后引入數據處理軟件。在信號研究中，誤差信號采用噪聲模型描述，本文應用最大似然法對吉林省連續運行參考站網絡(Jilin Continuously Operating Reference Station，JLCORS)參考站坐標時間序列噪聲性質進行研究。

一、數據資料及噪聲分析方法

1. JLCORS概況

CORS系統理論源于20世紀80年代，經過多年的發展，現在連續運行參考站系統(CORS)能夠常年連續不斷地運行，服務于高精度中短期天氣狀況的數值預報、變形監測、地震監測、地球動力學等[7]。它不僅滿足了各種測繪、基準的需求，還滿足了地震地殼形變監測需求。目前，很多國家建立了國家級、區域級的CORS網絡[8]，JLCORS由吉林省測繪地理信息局建設(如圖1所示)，圖中黑色正三角表示設計參考站位置。JLCORS是我國地區(省)級的連續運行的參考站網絡系統，它將在全省范圍內建立的永久性參考站通過網絡互聯，構成新一代的網絡化的大地測量系統，將GNSS技術綜合應用于吉林省的大地測量、工程測量、氣象監測、地震監測、地面沉降監測及城市地理信息系統等領域，2010年12月已經建成49個站并投入運行。

2. GNSS后處理軟件

GAMIT/GLOBK是國際上三大高精度GNSS后處理軟件之一，由美國麻省理工學院(MIT)和斯克里普斯海洋研究所(SIO)聯合開發，主要應用于分析研究地殼形變、高精度GPS測量數據處理等領域。GAMIT軟件處理雙差觀測量，采用最小二乘算法進行參數估計。GAMIT軟件主要由以下幾個模塊構成：ARC(軌道積分模塊)、MODEl(組成觀測方程)、SINCI N(單差自動修復周跳)、DBCI N(雙差自動修復周跳)、CVIEW(人工交互式修復周跳)、SOI VE(最小二乘解算模塊)、DFMRG(數據融合模塊)、FXDRV(生成批處理文件)、GI OBK(運用卡爾曼濾波進行網平差)等。

圖1　JLCORS基準站分布

3. 時間序列噪聲分析方法

時間序列是動態測試中被觀測量在一定的幅值范圍內隨時間出現的一系列隨機數據。頻譜分析和波形分析是最重要和最基本的方法。頻譜分析可求得時間序列的幅值譜、相位譜、功率譜和各種譜密度等，波形分析也即時間域分析，與頻譜分析可用傅里葉變換相互轉換。功率譜技術和最大似然法分別在頻率域和時間域對大地測量手段時間序列的噪聲分析表明，時間序列包含的噪聲類型非常復雜，如白噪聲、冪律噪聲、高斯-馬爾可夫噪聲、帶通濾波噪聲等，可以是僅含有一種噪聲，也可以是多種噪聲的組合。相應地用這些噪聲的數學模型擬合時間序列來求解噪聲參數。

最大似然估計(MLE)[9]可以同時獲得時間相關的噪聲模型的所有參數，經過試驗證實，MLE方法不但可以最大限度地發現一列隨機數據包含的噪聲，而且能準確計算出這些噪聲的功率譜，與譜分析所得的結果能夠精確符合。本次研究應用MLE方法進行時間序列噪聲分析。

(1) 最大似然估計

最大似然估計的公式為

(1)

在實際中，C是一種或多種噪聲統計模型的組合，可以用下式表示

(2)

式中，σ為噪聲模型的幅度；R為各噪聲模型的協方差矩陣。

式(2)取自然對數后可得

(3)

對σ求微分可得

(4)

通過上述變換，可以把多維問題轉化為一維問題，在解似然方程中，采用上述單純形法。

(2) 噪聲模型

根據地球物理觀測所包含噪聲統計規律，構造出如下模型，在噪聲分析中擬合時間序列噪聲。

① 冪律噪聲

因許多地球物理過程可以用一個冪律過程來描述，GNSS時間序列噪聲也具有這樣的冪律性質，即噪聲功率譜Px(f)與其對應的噪聲頻率f之間存在某種冪次關系

(5)

(6)

式中，Px為功率譜；f為噪聲頻率；P0、f0為常數；κ為譜指數；σpl為噪聲幅度。譜指數分布范圍通常為(-3，1),穩定過程譜指數為(-1，1)，不穩定過程譜指數為(-3，-1)。隨機過程分為穩定過程和不穩定過程，不穩定過程在低頻部分通常有較大的功率，而在高頻部分有負的譜指數，通常負譜指數分布在開區間(-3，-1)內，其中白噪聲的譜指數κ=0，閃爍噪聲的譜指數κ=-1、κ=1。類隨機游走過程的譜指數分布在區間(1，3)內,如典型的隨機游走噪聲譜指數κ=2。

② 高斯-馬爾可夫噪聲

一階高斯-馬爾可夫噪聲統計模型可用下式表示

(7)

一階高斯馬爾可夫信號可用下面微分方程表示

(8)

式中，Px為功率譜；f為噪聲頻率；P0、f0為常數；κ為譜指數。

③ 帶通噪聲

帶通信號是周期過程，在地球物理現象中常常是由季節變化引起，公式如下

(9)

式中，P為功率譜；f為噪聲頻率；f1為常數；κ為譜指數；σbp為噪聲幅度。

二、數據處理結果及噪聲性質分析

在獲得JLCORS坐標時間序列時，使用的數據策略是首先在GAMIT部分獲得所有站點的單日松弛解。在計算中選用IGS站(CHAN、BJFS、KHRJ、AIRA等)和陸態網站(HLAR、SUIY、JIXN、CHUN等)作為控制點，適當緊約束控制站坐標而對JLCORS站給予較松弛的約束。在GLRED模塊，先獲得所有點的時間序列，對數據結果進行檢查，利用GLOBK將SPOAC給出的全球單日松弛解和計算所得的區域單日松弛解進行綜合平差計算，在此基礎上通過IGS核心站求解相對于全球參考框架ITRF2008的相似變換參數，從而獲得ITRF2008下的單日解。最終獲得JLCORS參考站坐標時間序列，分別為E(東向)、N(北向)，U(垂直方向)，如圖2所示。

數據處理采用CATS程序[10]，其核心算法是參數估計中的最大似然估計。計算時，擬合多個經典噪聲統計模型到觀測時間序列，通過構造和求解似然函數，獲得最佳的噪聲模型，而且可以同時求解噪聲模型中的所有待估參數。計算是一個迭代的過程，為提高處理速度，CATS在計算過程中引入1個標量angle，即一個幫助計算的輔助量，沒有物理意義，通過angle的增加計算出似然函數的最大似然值，然后計算噪聲模型中各個參數的值。

圖2　JLCORS部分參考站(DGAN、DHUA、EDAO、FJTN)坐標時間序列

通過對冪律噪聲、閃爍噪聲、隨機游走噪聲及它們各自和白噪聲的組合進行似然函數值的計算，隨機游走噪聲和白噪聲的組合能取得最大的似然函數值。表1是觀測站N-S向坐標時間序列迭代后求得的噪聲組合模型中各個參數的估計值。表中σpl表示利用式(5)的冪律噪聲模型計算的噪聲幅度，σbp表示利用式(9)的帶通濾波模型計算的噪聲幅度，SD是各模型計算結果的標準差。表中帶通濾波噪聲很多站通過迭代后，似然函數不收斂，沒有計算結果，用“—”表示，其中有幾個站有迭代結果，其標準差也非常大，顯然是錯誤的結果。這樣的情況說明序列中不含帶通濾波噪聲，這類噪聲是由于季節更替產生的，存在于大多地球物理觀測資料中，在JLCORS參考站中卻沒有解算結果，所有站坐標時間序列均受白噪聲的影響。由于冪律噪聲模型覆蓋的噪聲類型范圍很廣，所有站結果可見主要含有冪律噪聲。

表1　JLCORS參考站噪聲參數分析結果

用冪律噪聲加白噪聲的組合噪聲模型擬合時間序列可獲得最大的似然函數值，同時獲得每個參考站所有分向冪律噪聲譜指數，分布區間為[-2.32，-1.71]，譜指數集中在κ=-2附近(如圖3所示)，此時為冪律噪聲的特例——隨機游走噪聲。由此可見GNSS時間序列噪聲占主導地位的應該是隨機游走噪聲。研究地殼運動的基巖站，隨機游走噪聲的幅度很小，較短的時間序列很難監測到，特別是在全球網解算數據時很大的共模誤差會將其淹沒[11]。本次研究的參考站網為區域級別，空間相關的誤差不會產生，有利于計算隨機游走噪聲的幅度。

圖3　JLCORS東西方向譜指數和噪聲幅度的關系

三、結束語

GNSS觀測資料結果解算后在結果中包含很多誤差信息，如GNSS整個觀測系統的誤差、解算過程的模型誤差、觀測墩的非構造誤差等。在信號處理過程中，這些信號可看作是隨機噪聲。本文通過最大似然法進行參數估計，根據JLCORS參考站坐標時間序列與經典噪聲模型擬合，估計各個噪聲模型的參數，來確定觀測結果的誤差信息可用哪種噪聲模型來描述。該計算過程使用噪聲分析軟件CATS來實現，經過分析表明，所有參考站均可用冪律噪聲和白噪聲的組合模型來描述，根據譜指數分布區間，冪律噪聲中的特殊情況之一的隨機游走噪聲占主導地位。也就是說JLCORS參考站坐標時間序列的誤差信息可用隨機游走噪聲和白噪聲的組合來描述。

[1]王敏, 沈正康, 董大南. 非構造形變對 GPS 連續站位置時間序列的影響和修正[J]. 地球物理學報, 2005, 48(5): 1045-1052.

[2]田云鋒, 沈正康. GPS 坐標時間序列中非構造噪聲的剔除方法研究進展[J]. 地震學報, 2009, 31(1): 68-81.

[3]黃立人, 符養. GPS 連續觀測站的噪聲分析[J]. 地震學報, 2007, 29(2): 197-202.

[4]LANGBEIN J. Noise in Two-color Electronic Distance Meter Measurements Revisited[J]. Journal of Geophysical Research: Solid Earth (1978—2012), 2004, 109(B4).DOI:10.1029/2003JB002819.

[5]LANGBEIN J, JOHNSON H. Correlated Errors in Geodetic Time Series: Implications for Time-dependent Deformation[J]. Journal of Geophysical Research: Solid Earth (1978—2012), 1997, 102(B1): 591-603.

[6]MAO A, HARRISON C G, DIXON T H. Noise in GPS Coordinate Time Series[J]. Journal of Geophysical Research: Solid Earth (1978—2012), 1999, 104(B2): 797-816.

[7]杜向鋒, 張興福, 張永毅. CORS 測量成果轉換的一步法及其精度分析[J]. 測繪通報, 2015(7): 23-26.

[8]鄭立平. 港口鎮單參考站 CORS 系統建設及應用[J]. 測繪通報, 2011(1): 43-45.

[9]莫惠棟. 最大似然法及其應用[J]. 遺傳, 1984, 6(5): 42-48.

[10]WILLIAMS S D. CATS: GPS Coordinate Time Series Analysis Software[J]. GPS solutions, 2008, 12(2): 147-153.

[11]BOS M, FERNANDES R, WILLIAMS S, et al. Fast Error Analysis of Continuous GPS Observations[J]. Journal of Geodesy, 2008, 82(3): 157-166.

Analysis of the Noise of Coordinate Time Series from JLCORS Network Reference Stations

LIU Junqing,XIAO Huifeng,ZHANG Chenxia,DING Guang

劉俊清，肖輝峰，張晨俠,等.JLCORS網絡參考站坐標時間序列噪聲類型分析[J].測繪通報，2016(10)：1-5.DOI:10.13474/j.cnki.11-2246.2016.0316.

2015-12-15

地震科技星火計劃(XH14016Y)

劉俊清(1977—)，男，博士生，高級工程師，主要從事地震火山監測預報研究。E-mail: lunwen_98@126.com

P228

B

0494-0911(2016)10-0001-05