基于QOCA軟件的GNSS連續站中非構造噪聲剔除方法研究

紀海源,何遠梅,郝 明

(1.陜西工業職業技術學院,陜西 咸陽 712000; 2.中國有色金屬工業西安勘察設計研究院有限公司,西安 710000; 3.中國地震局第二監測中心,西安 710043)

0 引言

連續GNSS變形監測是地殼形變監測的主要方法。研究發現,在連續站位移時間序列中存在著非地殼內部動力因素造成的非構造噪聲,這種非構造噪聲由很多因素引起,包括數據處理軟件中模型沒有消除的、地表質量負荷、隨著季節性下雨與下雪等負荷及觀測墩的季節性熱脹冷縮等。姚宜斌等[1]指出,GLOBK 軟件假定測站運動呈線性變化,與現實中觀測到的測站運動不符,可能會引入粗差。QOCA軟件在數據處理中能夠兼顧測站運動的非線性變化,運用QOCA軟件剔除GNSS連續站位移序列中的非構造噪聲,可為地殼的形變研究提供重要的基礎數據。QOCA(Quasi-Observation Combination)是空間大地測量資料處理軟件,可對各種松弛約束的大地測站坐標及速度解(如準觀測值)進行聯合解算,獲取地殼形變信息,不僅能夠對全球定位系統及衛星激光測距等空間大地擬觀測值進行解算,還可以與陸地大地測量擬觀測值(電子激光測距、三角測量、水準測量等)進行聯合解算,能夠聯合合成孔徑雷達數據(SAR)、重力、地震活動及地表螺變等數據進行解算。通過計算獲得測站坐標、測站速度、網參數及同震與震后形變參數,探測出包含在擬觀測值中的異常值,采取穩健性分析直接處理噪聲時間序列。

1 非構造噪聲

GAMIT/GLOBK處理得到的GNSS連續站位移時間序列是相對精確的,GAMIT中所加載的各種改正模型在很大程度上降低了噪聲的影響,但是不能完全消除,故這種位移序列結果中除了包含真實的地殼構造形變外,還包含了一些干擾噪聲,圖1為連續站GSJN的位移時間序列,在NEU三個方向明顯能看見周期性成分,這不是由地殼內部動力因素造成的。在GNSS監測過程中產生的干擾統稱為非構造噪聲,這種非構造噪聲由很多因素引起,包括GAMIT中模型沒有消除的地表質量負荷,隨著季節性下雨、下雪的負荷,觀測墩的季節性熱脹冷縮等[1]。地殼的運動不可能呈周期性變化(一般為線性),故去除其中的非構造噪聲信息對于有效運用GNSS數據研究構造形變場來說至關重要。

圖1 位移時間序列Fig.1 Displacement time series

2 analyze_tseri與pca模塊的數據分析原理

analyze_tseri與pca是QOCA軟件的主要組成模塊,對非構造噪聲的分析有重要的作用。

analyze_tseri為時間序列分析(time series analysis),是一種動態數據處理及統計方法,廣泛應用于多領域,可將來自地表與地下的策動源產生的確定性地形變模式(pattern)以參數化方式估計出來。把時間序列中未知的信號看作一隨機過程,利用自回歸分析、最大熵譜估計、ARMA模型擬合等統計方法把隱藏的信號分離出來,通常將地形變時間序列用數學模型表示為:

x(t)=x0+v·(t-t0)+∑[Sisin(ωit)

pca(principal component analysis)稱為主分量分析、經驗正交函數分析(EOF)、特征向量分析,能夠把隨時間變化的臺站網時間序列分解成時間域的主分量及空間域的特征向量,這種分解按照每個主分量貢獻的功率(能量)來排列,一般前幾個主分量的變化特征可以將區域性時間變化特征最大限度地表現出來,這幾個主分量對應的特征向量反映了這些時間變化強弱的空間分布[2]。如果把臺站網中每個臺站的時間序列排列起來,每個臺站時間序列為一列,可表示為如下m×n(通常情況下,m>n)的數字矩陣:

可以由奇異值分析分解成如下形式:X=UΠV,式中,U為m×m正交歸一矩陣,為m×n準對角矩陣,V為n×n正交歸一矩陣。地形變分析中最常見的情況是m>n,且X的秩為n。這時X的方差矩陣可表為:

3 QOCA軟件剔除非構造噪聲

為了剔除GNSS位移時間序列中的非構造噪聲,必須獲取這種噪聲的時間序列,分析其中包含的信息。根據黃立人等[3]的研究資料,把噪聲按照功率譜密度P(f)與其對應的噪聲頻率f之間的關系P(f)∝(f)k分類為:若k=0則對應的是白噪聲,與頻率沒有多大關系,各種頻率下噪聲大小一樣。若k=1,則對應的是閃爍噪聲,功率譜的密度與頻率成反比,有一定的波動。若k=2,則對應的噪聲稱為隨機漫步噪聲,有一定的隨機性,波動性大。噪聲性質的確定關鍵是獲得噪聲序列并求出它的功率譜密度(簡稱功率譜)與譜指數。對于連續的位移時間序列分量來說,可采用模型消減這種噪聲影響[3]。為了得到觀測噪聲的時間序列,對位移分量每日解觀測序列建立參數模型:

式中,ti(i=1…N)為以年為單位的時間,在這里以每日的單日解構成時間序列。a為初始位置,b為速率,s1,s2,s3、s4分別為年、半年周期項系數,g為偏移,h為震后速率變化,k為震后速率衰減(指數模型),H為階梯(heavisidestep)函數,cme為共模誤差,υ為誤差。以上數學模型可通過QOCA軟件的模塊來實現。

3.1 analyze_tseri分離線性成分

使用analyze_tseri模塊對GLOBK處理結果的36個連續站進行位移時間序列分析,擬合其中的線性項、周年、半周年項等參數,得到擬合殘差。

analyze_tseri輸入文件準備。jz_unf.drv:驅動文件,這是analyze_tseri模塊最重要的一個文件,包含所有輸入輸出文件的路徑、名字及各種選擇項的參數。analyze_tseri會根據驅動文件中的參數來運行。文件由用戶建立,采用自由格式,主要文件包括coord.soln(臺站先驗坐標和速度文件)、stxyz.list(輸入數據文件的列表文件)、jzsite.list(臺站文件)、jz_tseri.para(估計參數文件)。其中,jz_tseri.para文件給出了每個臺站要估計的參數及每個參數的情況,這是程序最關鍵的文件,參數擬合與估計的好壞全看它給的是否合理。

運行analyze_tseri指令,給出驅動文件名字(analyze_tseri analyze_tseri.drv)。

結果文件。jz_unf.out:輸出結果文件。頭文件中必須給出輸出文件名及路徑。這個文件包含了所有估計參數的結果,其中以大寫字母開頭的一般是估計參數解。小寫字母行的意義如下:lsq_velo(速度項解)、lsq_coor(臺站坐標解)、coor_xyz與velo_xyz(坐標與速度估計在xyz坐標系中的解)、jz_unf.resi(殘差文件)。jz_unf.resi給出了第一個殘差文件的名字及路徑,考慮為后面的區域濾波做準備,必須輸出殘差文件。圖2是根據參差文件jz_unf.resi畫的殘差時間序列。

3.2 pca分離共模誤差

由analyze_tseri處理得到的殘差圖2可以看出,analyze_tseri所得的殘差時間序列中仍可見周期性成分及不規則波動。根據Wdowinski S等[4]的研究表明,殘差中包含的這類非構造噪聲稱為共模誤差(common modeerror,即cme),這類噪聲存在區域相關性,來源尚不清楚,有可能來自衛星軌道誤差、水體及大氣質量負荷、參考框架定義的不確定性等[5],此外還存在單站誤差,源于天線墩不穩定、多路徑效應等[6]。Dong等[5]采用pca方法對臺站網中是否存在cme制定了標準,如果大于50%的站出現明顯的空間響應,且這種模式特征值與特征值總和大于1%,按照這個標準一般構成cme的主體部分是排在前面的幾個主成分。SCIGN通過pca分析的殘差坐標序列文件表明,cme具有空間特征,以同一方式在一定的區域內影響各站點。因為這種影響不是地殼構造活動造成的,故有必要從坐標時間序列中剔除[7-9]。采用pca模塊對analyze_tseri的殘差文件得到的殘差坐標序列進行主成分分析。

pca分析文件準備。pca_jz.drv:pca的驅動文件,裝載所有輸入輸出文件的路徑及名字及各種選擇項的參數,這是一個必須由用戶建立的文件,有3種類型,即regional filtering、tectonic signal、seasonal signal,pca的內插方式是不一樣的。對于regional filtering工作類型,程序要求內插盡量不影響由沒有空缺觀測的臺站給出的前1～2個主分量,對于后面兩種工作類型,程序要求內插盡量不影響由沒有空缺觀測的臺站給出的主要構造運動模式即季節性變化模式。處理選擇regional filtering,需準備臺站先驗坐標及速度文件coord.soln,臺站文件site_list,輸入數據文件的列表文件pca_jz.list。

運行pca指令,給出驅動文件名字pca pca_jz.drv即可。

輸出結果如下:

pca_jz.out:輸出結果文件,給出輸出文件名及路徑,主要儲存各個分量分解的本征值及其功率。

pca_jz.cpt:主分量文件,給出主分量文件的名字及路徑。輸出12個主分量(單位a,m)按東西、南北、垂直向依次排列如下:

* time cpt1 cpt2 cpt3 cpt4 cpt5 cpt6 cpt7 cpt8 cpt9 cpt10 cpt11 cpt12

2011.00100000 -0.004025 -0.001041 -0.000234 0.001173 0.001535 0.000044 -0.005732 -0.002654 -0.000331 -0.000613

-0.002197 0.000410

pca_jz.seign:輸出特征向量文件,給出空間響應本征矢文件名字及路徑。輸出的前10個特征向量(無單位)按東西、南北、垂直向依次排列:

* long lati site cpt1 cpt2 cpt3 cpt4 cpt5 cpt6 cpt7 cpt8 cpt9 cpt10

104.6046 35.5543 GSDX_GNSS 0.758543 0.201866 0.131942 -0.020168 0.021891 -0.002359 -0.057869 -0.030266 0.434962 0.214931

pca_jz.res:濾波后殘差序列輸出文件,必須給出輸出文件名及路徑。通過pca的主分量分析得到了兩個重要文件,即pca_jz.cpt、pca_jz.seign,包含主要的共模誤差,可作為analyze_tseri的輸入文件。

3.3 analyze_tseri剔除非構造噪聲

經過pca主成分分析得到了主分量文件pca_jz.cpt及本征矢文件pca_jz.seign,將其加入到analyze_tseri的驅動文件中,可以剔除pca分離出來的共模誤差。對剔除共模誤差后的位移時間序列中的周年項、半周年項等非構造噪聲參數信息進行擬合,從位移序列中剔除。QOCA會根據剔除干凈的位移時間序列來擬合一個最終速度,這個速度即GNSS連續站在扣除非構造噪聲后相對真實的線性運動速度。analyze_tseri處理應注意在驅動文件jz_flt.drv中的關鍵字cme_correction file行后面添加濾波文件pca_jz.seign pca_jz.cpt。

3.4 QOCA處理結果分析

經過QOCA軟件對GNSS連續站的處理,已經將大量的非構造噪聲剔除,得到需要的各種估計參數及速度文件。表1是從最終的結果文件jz_flt.out中提取的部分連續站的速度。

表1 剔除非構造噪聲后的部分連續站速度統計Tab.1 Partial continuous station velocity statistics after removing unstructured nois

圖3是用連續站GSJN剔除線性成分與共模誤差后的殘差文件畫的,從圖中明顯能看到比濾波前的波動小,但是出現個別點誤差的增大,這可能是因為個別單日解的誤差較大。通過pca方法得到的結果與Dong等[7]的結果基本一致。

圖3 剔除線性成分和共模誤差后的GSJN站殘差Fig.3 Residual of GSJN station after removing linear component and common mode error

將剔除非構造噪聲前后的速度進行對比,表2列出了部分連續站的數據,在計算速度差均方根時采用全部36個測站,經過對比發現,在南北方向上最大速度差為0.6 mm,最小為0 mm,速度差均方根為0.26 mm。在東西方向上最大差0.5 mm,最小差0 mm,速度差均方根0.19 mm。在垂直方向上最大-1.8 mm,最小為-0.3 mm,速度差均方根0.64 mm。水平方向上的速度差均方根為0.33 mm,可見,垂直方向上速度差均方根是水平方向上的1.97倍。這符合GNSS高程測量精度小于水平方向上的規律。QOCA處理結果在水平方向還不是太明顯,而在垂直方向上對提高GNSS高程測量有很大的作用。

表2 剔除非構造噪聲前后部分連續站速度對比Tab.2 Velocity comparison of partial continuous stations before and after removing unstructured noise

4 結論

使用QOCA軟件的analyze_tseri與pca模塊對36個GNSS連續站位移時間序列進行分析處理,剔除連續站中的共模誤差,獲取了位移時間序列中的周期性參數信息,將其中的周期性等非構造噪聲剔除,令GNSS連續站水平運動速度的精度提高了約0.3 mm/yr,垂直方向運動速度的精度提高了約0.6 mm/yr,得到相對真實的地殼形變信息,為地殼的形變研究提供了重要的基礎數據。