紀(jì)海源,何遠(yuǎn)梅,郝 明

(1.陜西工業(yè)職業(yè)技術(shù)學(xué)院,陜西 咸陽 712000; 2.中國有色金屬工業(yè)西安勘察設(shè)計(jì)研究院有限公司,西安 710000; 3.中國地震局第二監(jiān)測(cè)中心,西安 710043)

0 引言

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

1 非構(gòu)造噪聲

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

圖1 位移時(shí)間序列Fig.1 Displacement time series

2 analyze_tseri與pca模塊的數(shù)據(jù)分析原理

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

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

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

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

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

3 QOCA軟件剔除非構(gòu)造噪聲

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

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

3.1 analyze_tseri分離線性成分

使用analyze_tseri模塊對(duì)GLOBK處理結(jié)果的36個(gè)連續(xù)站進(jìn)行位移時(shí)間序列分析,擬合其中的線性項(xiàng)、周年、半周年項(xiàng)等參數(shù),得到擬合殘差。

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

運(yùn)行analyze_tseri指令,給出驅(qū)動(dòng)文件名字(analyze_tseri analyze_tseri.drv)。

結(jié)果文件。jz_unf.out:輸出結(jié)果文件。頭文件中必須給出輸出文件名及路徑。這個(gè)文件包含了所有估計(jì)參數(shù)的結(jié)果,其中以大寫字母開頭的一般是估計(jì)參數(shù)解。小寫字母行的意義如下:lsq_velo(速度項(xiàng)解)、lsq_coor(臺(tái)站坐標(biāo)解)、coor_xyz與velo_xyz(坐標(biāo)與速度估計(jì)在xyz坐標(biāo)系中的解)、jz_unf.resi(殘差文件)。jz_unf.resi給出了第一個(gè)殘差文件的名字及路徑,考慮為后面的區(qū)域?yàn)V波做準(zhǔn)備,必須輸出殘差文件。圖2是根據(jù)參差文件jz_unf.resi畫的殘差時(shí)間序列。

3.2 pca分離共模誤差

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

pca分析文件準(zhǔn)備。pca_jz.drv:pca的驅(qū)動(dòng)文件,裝載所有輸入輸出文件的路徑及名字及各種選擇項(xiàng)的參數(shù),這是一個(gè)必須由用戶建立的文件,有3種類型,即regional filtering、tectonic signal、seasonal signal,pca的內(nèi)插方式是不一樣的。對(duì)于regional filtering工作類型,程序要求內(nèi)插盡量不影響由沒有空缺觀測(cè)的臺(tái)站給出的前1～2個(gè)主分量,對(duì)于后面兩種工作類型,程序要求內(nèi)插盡量不影響由沒有空缺觀測(cè)的臺(tái)站給出的主要構(gòu)造運(yùn)動(dòng)模式即季節(jié)性變化模式。處理選擇regional filtering,需準(zhǔn)備臺(tái)站先驗(yàn)坐標(biāo)及速度文件coord.soln,臺(tái)站文件site_list,輸入數(shù)據(jù)文件的列表文件pca_jz.list。

運(yùn)行pca指令,給出驅(qū)動(dòng)文件名字pca pca_jz.drv即可。

輸出結(jié)果如下:

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

pca_jz.cpt:主分量文件,給出主分量文件的名字及路徑。輸出12個(gè)主分量(單位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:輸出特征向量文件,給出空間響應(yīng)本征矢文件名字及路徑。輸出的前10個(gè)特征向量(無單位)按東西、南北、垂直向依次排列:

* 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的主分量分析得到了兩個(gè)重要文件,即pca_jz.cpt、pca_jz.seign,包含主要的共模誤差,可作為analyze_tseri的輸入文件。

3.3 analyze_tseri剔除非構(gòu)造噪聲

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

3.4 QOCA處理結(jié)果分析

經(jīng)過QOCA軟件對(duì)GNSS連續(xù)站的處理,已經(jīng)將大量的非構(gòu)造噪聲剔除,得到需要的各種估計(jì)參數(shù)及速度文件。表1是從最終的結(jié)果文件jz_flt.out中提取的部分連續(xù)站的速度。

表1 剔除非構(gòu)造噪聲后的部分連續(xù)站速度統(tǒng)計(jì)Tab.1 Partial continuous station velocity statistics after removing unstructured nois

圖3是用連續(xù)站GSJN剔除線性成分與共模誤差后的殘差文件畫的,從圖中明顯能看到比濾波前的波動(dòng)小,但是出現(xiàn)個(gè)別點(diǎn)誤差的增大,這可能是因?yàn)閭€(gè)別單日解的誤差較大。通過pca方法得到的結(jié)果與Dong等[7]的結(jié)果基本一致。

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

將剔除非構(gòu)造噪聲前后的速度進(jìn)行對(duì)比,表2列出了部分連續(xù)站的數(shù)據(jù),在計(jì)算速度差均方根時(shí)采用全部36個(gè)測(cè)站,經(jīng)過對(duì)比發(fā)現(xiàn),在南北方向上最大速度差為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高程測(cè)量精度小于水平方向上的規(guī)律。QOCA處理結(jié)果在水平方向還不是太明顯,而在垂直方向上對(duì)提高GNSS高程測(cè)量有很大的作用。

表2 剔除非構(gòu)造噪聲前后部分連續(xù)站速度對(duì)比Tab.2 Velocity comparison of partial continuous stations before and after removing unstructured noise

4 結(jié)論

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