馬超,李斐,2,張勝凱*,雷錦韜,鄂棟臣,郝衛(wèi)峰,張卿川

1 武漢大學(xué)中國(guó)南極測(cè)繪研究中心,武漢　430079 2 武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,武漢　430079

?

顧及共性誤差的南極半島地區(qū)連續(xù)GPS站坐標(biāo)時(shí)間序列分析

馬超1,李斐1,2,張勝凱1*,雷錦韜1,鄂棟臣1,郝衛(wèi)峰1,張卿川1

1 武漢大學(xué)中國(guó)南極測(cè)繪研究中心,武漢430079 2 武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,武漢430079

共性誤差是區(qū)域連續(xù)GPS網(wǎng)絡(luò)中存在的一種與時(shí)空相關(guān)的主要誤差源.對(duì)區(qū)域網(wǎng)絡(luò)使用空間濾波能夠有效減小這種誤差并改善坐標(biāo)時(shí)間序列的精度.本文使用GAMIT/GLOBK 10.5軟件解算南極半島區(qū)域11個(gè)GPS測(cè)站2010—2014的數(shù)據(jù)，各測(cè)站獨(dú)立估計(jì)的單天解組成坐標(biāo)時(shí)間序列.利用堆棧法(stacking)、主分量分析法(principal component analysis,簡(jiǎn)稱PCA)和Karhunen-Loeve展開法(KLE)對(duì)其進(jìn)行空間濾波，在利用PCA和KLE方法進(jìn)行空間濾波時(shí)，使用前兩個(gè)主分量作為共有模式分量.結(jié)果表明：三種方法均能有效提取南極半島區(qū)域共性誤差，PCA方法的濾波效果明顯優(yōu)于stacking方法，略微優(yōu)于KLE方法；空間濾波能夠有效減小殘差時(shí)間序列的振幅、功率和RMS；對(duì)坐標(biāo)時(shí)間序列進(jìn)行空間濾波能夠有效降低其線性項(xiàng)和周期項(xiàng)誤差，從而提高線性項(xiàng)和周期項(xiàng)估計(jì)的精度.共性誤差的譜分析結(jié)果顯示在U方向上存在9.4天、13.7天等短周期的信號(hào)，推測(cè)共性誤差組成成分中可能含有與海潮相關(guān)的誤差源.

南極半島;GPS;時(shí)間序列;共性誤差

1　引言

全球連續(xù)GPS網(wǎng)絡(luò)在不同的空間和時(shí)間尺度上為形變模式的監(jiān)測(cè)提供了統(tǒng)一框架.連續(xù)GPS測(cè)量能夠在全球板塊運(yùn)動(dòng)、板塊邊界區(qū)域地殼形變以及局部同震變形等一系列空間尺度上揭示形變模式.形變模式的量化有助于從冰后回彈、季節(jié)性形變、瞬變運(yùn)動(dòng)等不同時(shí)間尺度來(lái)研究地球物理過(guò)程.在區(qū)域連續(xù)GPS網(wǎng)絡(luò)中，存在著一種與時(shí)空相關(guān)的誤差，稱為共性誤差(也稱共模誤差，Common Mode Error).共性誤差的有效提取與剔除對(duì)于提高區(qū)域連續(xù)GPS坐標(biāo)時(shí)間序列精度、改進(jìn)形變模式分析具有重要的作用.

Dong等(2006)使用PCA和KLE聯(lián)合濾波方法對(duì)南加利福尼亞綜合GPS網(wǎng)絡(luò)(SCIGN)5年單天坐標(biāo)時(shí)間序列進(jìn)行濾波，比stacking方法更有效地提取了共性誤差.伍吉倉(cāng)等(2008)和胡守超等(2009)也對(duì)該區(qū)域的共性誤差進(jìn)行了提取，前者使用了PCA和KLE法，后者還對(duì)stacking、PCA和KLE三種方法進(jìn)行了比較.蔣志浩等(2010)采用PCA方法提取了CGCS2000框架下國(guó)家連續(xù)進(jìn)行參考站系統(tǒng)(CORS)網(wǎng)絡(luò)坐標(biāo)時(shí)間序列中的共性誤差，其結(jié)果表明使用空間濾波能夠有效降低N、E、U三個(gè)方向的平均坐標(biāo)重復(fù)性和U方向周年項(xiàng)振幅.Serpelloni等(2013)利用800多個(gè)連續(xù)GPS站點(diǎn)的坐標(biāo)時(shí)間序列來(lái)研究歐洲地中海地區(qū)的垂向變形速率，利用PCA方法估計(jì)并剔除了坐標(biāo)時(shí)間序列中的共性誤差，使位移數(shù)據(jù)的信噪比有了明顯提高.Shen等(2014)針對(duì)不連續(xù)的GPS時(shí)間序列提出了改進(jìn)的PCA方法，從中國(guó)地殼運(yùn)動(dòng)觀測(cè)網(wǎng)絡(luò)(CMONOC)不完整的時(shí)間序列中提取了共性誤差.雖然目前有關(guān)共性誤差的研究已日趨成熟，但仍然存在一些問(wèn)題值得探討，例如當(dāng)測(cè)站較少時(shí)，PCA或KLE方法得到的高階主分量模式是否可以作為共有模式來(lái)提取共性誤差.

一直以來(lái)，南極區(qū)域連續(xù)GPS網(wǎng)絡(luò)由于受各種不利因素影響，測(cè)站稀少、數(shù)據(jù)質(zhì)量差且空間分布十分不均勻.2007/08年國(guó)際極地年間，由多國(guó)合作開始在南極實(shí)施POLENET計(jì)劃(The Polar Earth Observing Network)，該計(jì)劃的主要任務(wù)是在極地布設(shè)GPS與地震儀等應(yīng)用于地球科學(xué)研究的常年連續(xù)觀測(cè)儀器(鄂棟臣和張勝凱，2006).近年來(lái)，隨著POLENET計(jì)劃的實(shí)施和南極IGS測(cè)站的不斷擴(kuò)充，南極半島地區(qū)的連續(xù)GPS站在數(shù)量和空間分布上有了很大大改善，相對(duì)南極其他區(qū)域更為密集，研究該地區(qū)GPS網(wǎng)絡(luò)的共性誤差成為可能.

2　共性誤差提取方法

目前提取共性誤差一般使用空間濾波方法，主要包括stacking、PCA和KLE方法.提取共性誤差需要?dú)埐顣r(shí)間序列，各測(cè)站的殘差時(shí)間序列由坐標(biāo)時(shí)間序列剔除線性項(xiàng)和周期項(xiàng)得到，公式如下：

vi=y(ti)-(a+bti+csin(2πti)+dcos(2πti)

(1)

其中，vi為第i個(gè)歷元的觀測(cè)噪聲，也就是殘差；y(ti)為站坐標(biāo)時(shí)間序列；ti代表第i個(gè)歷元時(shí)間，以小數(shù)年的形式表示；a代表初始?xì)v元測(cè)站坐標(biāo)，b代表線性速度，c、d代表全年周期運(yùn)動(dòng)系數(shù)，e、f代表半年周期運(yùn)動(dòng)系數(shù)，g表示歷元Tg處由各種原因引起的階躍式偏移量(如天線位置變化，同震位移等)；H表示Heaviside階梯函數(shù).

2.1堆棧法(stacking)

假設(shè)一個(gè)區(qū)域網(wǎng)絡(luò)有n個(gè)測(cè)站，站坐標(biāo)時(shí)間序列長(zhǎng)度為m天，建立一個(gè)m×n的矩陣X，其元素x(ti,sj)代表第i天第k個(gè)測(cè)站的坐標(biāo)殘差值，那么該網(wǎng)絡(luò)的共性誤差可以表示為

(2)

其中，δi,k是第i天第k個(gè)測(cè)站的中誤差(Nikolaidis,2002).由(2)式可知，該方法提取的共性誤差即為各測(cè)站殘差的加權(quán)平均值，故也稱其為加權(quán)平均法.

堆棧法假設(shè)共性誤差在空間均勻分布，能夠?qū)?shù)百公里范圍GPS網(wǎng)絡(luò)有很好的近似.但隨著空間范圍擴(kuò)大，共性誤差會(huì)逐漸減弱至消失.

2.2PCA方法(PCA)

PCA方法是一種廣義的空間濾波方法，在同震變形的研究中該方法已成功用于大地測(cè)量數(shù)據(jù)的分解(Scherneck et al.,2000)和區(qū)域?yàn)V波(Johansson et al.,2002).某一區(qū)域殘差時(shí)間序列X為m×n矩陣，可定義其協(xié)方差陣B為

(3)

B是一個(gè)n×n的對(duì)稱矩陣，對(duì)B進(jìn)行正交分解可得

(4)

其中，VT是一個(gè)行向量正交的n×n特征向量矩陣,Λ是具有r個(gè)非零特征值的對(duì)角矩陣(n≥r).對(duì)于大多數(shù)大地測(cè)量數(shù)據(jù)，矩陣B是滿秩的，即n=r.由線性代數(shù)可知，任何秩為n的矩陣都由n維正交基擴(kuò)充得到.因此，可以使用正交基V來(lái)擴(kuò)充矩陣X：

(5)

即

2.1兩組患者臨床療效對(duì)比 觀察組患者治療顯效29例,有效8例,無(wú)效3例,總有效37例,總有效率為92.5%;對(duì)照組患者治療顯效25例,有效5例,無(wú)效10例,總有效30例,總有效率為75.0%。觀察組總有效率明顯高于對(duì)照組,差異具有統(tǒng)計(jì)學(xué)意義(P<0.05)。

(6)

其中，A是m×n矩陣，vk,j是第j個(gè)特征向量第k個(gè)元素.由于正交矩陣的性質(zhì)，(5)式可以變?yōu)?/p>

(7)

即

(8)

(6)式中ak(t)稱為矩陣X的第j個(gè)主分量，代表時(shí)間變化；vk(s)是其對(duì)應(yīng)特征向量，代表對(duì)主分量ak(t)的空間響應(yīng).上述分解方法稱為經(jīng)驗(yàn)正交函數(shù)(EOF)分析，也稱為主分量分析(PCA)(Preisendorfer,1988).將特征向量按照特征值大小進(jìn)行降序排列，排在前面的幾個(gè)主分量對(duì)殘差時(shí)間序列X的貢獻(xiàn)最大，通常與共源時(shí)間函數(shù)有關(guān)；而排在后面的主分量通常與當(dāng)?shù)鼗騻€(gè)別測(cè)站的影響有關(guān).那么，PCA方法的共性誤差可以定義為

(9)

其中，p為定義共性誤差的主分量個(gè)數(shù).由(9)式可知，PCA方法定義的共性誤差在不同的測(cè)站上是不同的，它舍棄了堆棧法空間均勻分布的假設(shè)，而讓數(shù)據(jù)本身去揭示共性誤差的空間分布.

2.3Karhunen-Loeve展開法(KLE)

該方法與PCA方法類似，只是將協(xié)方差矩陣B標(biāo)準(zhǔn)化得到相關(guān)系數(shù)矩陣C，其元素：

(10)

其中，δ是殘差時(shí)間序列矩陣X的標(biāo)準(zhǔn)差向量，即

(11)

相關(guān)系數(shù)矩陣C也可以正交分解為

(12)

ΛC、W分別是C的特征值組成的對(duì)角矩陣和特征向量組成的正交矩陣.

進(jìn)而得到主分量矩陣A：

A=XW,

(13)

同理，Karhunen-Loeve展開法的共性誤差可以定義為

(14)

3　共性誤差的提取

3.1GPS數(shù)據(jù)處理

為了提取和研究南極半島區(qū)域的共性誤差，本文選取了該區(qū)域11個(gè)GPS連續(xù)跟蹤站，包括5個(gè)IGS測(cè)站和6個(gè)POLENET測(cè)站，時(shí)間跨度為2010—2014年(5年)，測(cè)站分布如圖1所示.采用GAMIT/GLOBK10.5對(duì)GPS數(shù)據(jù)進(jìn)行解算和平差.GAMIT基線解算策略如下：解算類型為松弛解(RELAX)，同時(shí)解算衛(wèi)星軌道、測(cè)站坐標(biāo)；觀測(cè)值使用無(wú)電離層組合、自動(dòng)修復(fù)周跳模式(LC_AUTCLN)；測(cè)站N、E、U方向先驗(yàn)坐標(biāo)約束：IGS站約束為0.05 m、0.05 m、0.10 m，其他站約束為100 m、100 m、100 m；截止高度角為15°；歷元間隔為30 s；對(duì)流層延遲改正使用Saastamoinen模型；潮汐改正使用FES2004模型；太陽(yáng)光壓攝動(dòng)改正使用BERNE模型；參考框架為ITRF2008；其他均采用GAMIT10.5默認(rèn)設(shè)置.

圖1　南極半島區(qū)域GPS測(cè)站分布圖，圓點(diǎn)表示IGS站，三角形表示POLENET站Fig.1　The distribution of GPS stations in Antarctic Peninsula,the dots represent the IGS stations,and the triangles represent the POLENET stations

3.2殘差時(shí)間序列

得到單天坐標(biāo)時(shí)間序列之后，通過(guò)公式(1)給出的模型來(lái)剔除趨勢(shì)項(xiàng)和周期項(xiàng)等，從而得到單天殘差時(shí)間序列.具體方法是：

1)將GLOBK平差結(jié)果(org文件)利用sh_plotcrd 的-mb命令生成各測(cè)站的mb文件，再利用sh_mb2cats命令將mb文件轉(zhuǎn)換成cats格式文件.

2)將cats格式文件輸入到CATS軟件，噪聲模型選擇白噪聲+閃爍噪聲+隨機(jī)漫步噪聲(李偉偉等,2014)，估計(jì)各測(cè)站的6個(gè)系數(shù)(公式(1)中a、b、c、d、e、f).

3)將各測(cè)站的系數(shù)帶入公式(1)，由于在平差中使用地震重命名文件改正了階躍項(xiàng)，所以階躍項(xiàng)不再考慮，從而得到各測(cè)站的單天殘差時(shí)間序列.

3.3異常值剔除和時(shí)間序列插值

時(shí)間序列中的異常值會(huì)對(duì)時(shí)間序列分析造成偏差，所以有必要對(duì)異常值進(jìn)行剔除.首先對(duì)坐標(biāo)時(shí)間序列使用三倍四分位數(shù)差法則進(jìn)行孤立值剔除(Nikolaidis,2002)，然后剔除線性項(xiàng)和周期項(xiàng)，得到殘差時(shí)間序列.由于殘差時(shí)間序列是在0附近波動(dòng)，因此采用閾值法，超過(guò)某一閾值的的殘差被當(dāng)作粗差予以剔除.剔除標(biāo)準(zhǔn)是：N、E、U三個(gè)分量的中誤差閾值分別設(shè)為50、50、100 mm，殘差閾值設(shè)為100、100、300 mm.殘差時(shí)間序列中還存在間隙，無(wú)法直接提取共性誤差，需要進(jìn)行插值.本文參考Dong 等(2006)的插值方法，結(jié)合實(shí)際情況(測(cè)站較少)對(duì)其插值方法進(jìn)行了相應(yīng)修改，以減少迭代次數(shù).具體方法是：首先使用三點(diǎn)拉格朗日法(間隙小于等于3天)和空間平均法(間隙大于三天)對(duì)站坐標(biāo)時(shí)間序列進(jìn)行初步插值；然后使用PCA方法對(duì)初步插值的時(shí)間序列進(jìn)行主分量分析，取前三個(gè)主分量與缺失測(cè)站空間響應(yīng)乘積的和作為該測(cè)站新的插值；迭代這一過(guò)程直至前、后結(jié)果之差的平均值小于0.01 mm.

3.4共性誤差提取

堆棧法一般不考慮測(cè)站的特異性，所有測(cè)站的殘差按中誤差加權(quán)平均即可得到共性誤差，接下來(lái)主要介紹PCA和KLE方法的提取過(guò)程.為了使不同主分量之間便于比較，我們把每個(gè)特征向量除以其絕對(duì)值最大元素，得到標(biāo)準(zhǔn)化空間特征向量，也稱標(biāo)準(zhǔn)化空間響應(yīng)，正值代表對(duì)主分量的積極響應(yīng)，負(fù)值代表對(duì)主分量的消極響應(yīng).將主分量乘以對(duì)應(yīng)標(biāo)準(zhǔn)化因子，可得到標(biāo)準(zhǔn)化的主分量.

使用PCA方法對(duì)11個(gè)測(cè)站的殘差時(shí)間序列進(jìn)行主分量分析，分別得到了N、E、U方向第一、二、三的標(biāo)準(zhǔn)化主分量時(shí)間序列及其對(duì)應(yīng)的標(biāo)準(zhǔn)化空間向量，如圖2中a、b、c所示.通過(guò)對(duì)比三個(gè)方向的標(biāo)準(zhǔn)化空間向量，發(fā)現(xiàn)N方向第三主分量(圖2c)OHI2、OHI3、ROTH三個(gè)測(cè)站具有異常大的空間特征向量，其絕對(duì)值比其他測(cè)站平均大了約80%，且OHI2、OHI3與ROTH空間特征向量的方向相反，在E方向第三主分量模式以及U方向第二、三主分量模式也發(fā)現(xiàn)了類似的現(xiàn)象.考慮到OHI2、OHI3兩個(gè)測(cè)站位于南極半島東北端，而ROTH位于南極半島西南部的阿德萊德島，在空間上與其他測(cè)站的距離都相距較遠(yuǎn)，因此這三個(gè)測(cè)站有可能存在與其他測(cè)站不一致的強(qiáng)烈局部效應(yīng)，這里先將它們從GPS網(wǎng)絡(luò)中剔除，以便進(jìn)行下一步工作.利用KLE方法進(jìn)行主分量分析，得到的N、E、U三個(gè)方向的第三主分量模式(圖2d)與PCA第三主分量模式十分相似，從某種程度上也驗(yàn)證了PCA方法的適用性.

圖2　11個(gè)測(cè)站N、E、U方向的標(biāo)準(zhǔn)化主分量序列(上)及其對(duì)應(yīng)的標(biāo)準(zhǔn)化空間響應(yīng)(下)，其中a、b、c列為PCA分析得到的第一、二、三主分量模式，d列為KLE方法分析的第三主分量模式.Fig.2　The scaled principal components and normalized spatial eigenvectors of 11 stations in north,east and up component,which a,b,c are the first,second,third mode of PCA solution respectively,and d is the third mode of KLE solution

利用PCA和KLE方法對(duì)剩余8個(gè)測(cè)站重新進(jìn)行主分量分析，得到N、E、U方向的第一、二標(biāo)準(zhǔn)化主分量序列及其標(biāo)準(zhǔn)化空間響應(yīng)，如圖3所示.在剔除OHI2、OHI3、ROTH之后，兩種方法的結(jié)果中N、E、U三個(gè)方向的第一標(biāo)準(zhǔn)化主分量序列振幅都有所減小，但標(biāo)準(zhǔn)化空間響應(yīng)更加一致.在PCA方法得到的第二主分量模式中，南極半島東部和西部測(cè)站呈現(xiàn)出方向相反的空間響應(yīng)，KLE方法也發(fā)現(xiàn)類似的現(xiàn)象，說(shuō)明南極半島可能存在著對(duì)東部和西部影響相反的某種誤差源，再加上該模式中位于內(nèi)陸的FONP測(cè)站的標(biāo)準(zhǔn)化空間響應(yīng)遠(yuǎn)小于其他沿海測(cè)站，本文認(rèn)為這種現(xiàn)象可能是該區(qū)域海洋潮汐未完全改正或改正過(guò)度引起的.

圖3　8個(gè)測(cè)站N方向的標(biāo)準(zhǔn)化主分量(上)及其對(duì)應(yīng)的標(biāo)準(zhǔn)化空間向量(下)，其中a、c列為PCA方法得到的第一、二主分量模式，b、d列為KLE方法得到的第一、二主分量模式Fig.3　The scaled principal components and its normalized spatial eigenvectors of 8 stations in north component,which a,c are the first,second mode of PCA solution respectively,and b,d are the first,second mode of KLE solution respectively

對(duì)共性誤差的定義，尚無(wú)明確的標(biāo)準(zhǔn)和共識(shí).Dong等(2006)利用共有模式來(lái)定義共性誤差，即某一主分量模式中大多數(shù)測(cè)站(50%)具有明顯(>25%)的標(biāo)準(zhǔn)化空間響應(yīng)，且該模式的特征值超過(guò)了所有特征值總和的1%，可認(rèn)為是共有模式.根據(jù)這個(gè)定義，前幾階主分量模式都有可能是共有模式，而高階模式通常只與部分測(cè)站有關(guān)并且可能反映局部效應(yīng)，因此不會(huì)成為共有模式.為了確定PCA和KLE方法提取共性誤差要使用的共有模式主分量個(gè)數(shù)p，本文參考Dong等(2006)共有模式標(biāo)準(zhǔn)，對(duì)各主分量模式從標(biāo)準(zhǔn)化空間響應(yīng)和特征值兩方面分別進(jìn)行顯著性檢驗(yàn).在標(biāo)準(zhǔn)化空間響應(yīng)方面，由圖3也可以看出，利用PCA和KLE方法得到的第一、二主分量均能夠滿足50%以上的測(cè)站標(biāo)準(zhǔn)化空間響應(yīng)大于25%的條件，而更高階主分量并不滿足該條件，這里不再列出.在特征值方面，表1給出了第一、二主分量特征值占總特征值的百分比，可以看出PCA和KLE分析得到的第一、二主分量模式均滿足特征值百分比大于1%的條件.為了更加客觀、準(zhǔn)確地提取共性誤差，本文還分析了各主分量對(duì)殘差的貢獻(xiàn)，主分量模式特征值的百分比反映了該主分量模式對(duì)整個(gè)殘差時(shí)間序列的貢獻(xiàn)，圖4給出了兩種方法得到的各主分量特征值的累計(jì)百分比，兩種方法特征值累計(jì)百分比都在80%左右，且前兩個(gè)主分量的特征值遠(yuǎn)大于其他主分量的特征值，由此可見(jiàn)前兩階主分量對(duì)整個(gè)殘差時(shí)間序列的貢獻(xiàn)十分可觀.總體來(lái)說(shuō)，在PCA和KLE方法中第一、二主分量模式都可以作為共有模式.將PCA和KLE方法獲取的第一、二主分量序列及其對(duì)應(yīng)空間響應(yīng)分別代入公式(9)和(14)，其中p=2，就得到了相應(yīng)的共性誤差.

表1　第一、二主分量特征值占總特征值的百分比

由于堆棧法提取的共性誤差是空間均勻分布，而PCA方法和KLE方法得到的共性誤差在空間上并非均勻分布，為了便于與堆棧法進(jìn)行比較，對(duì)后兩者得到的共性誤差在空間上求平均，從而得到了與堆棧法一致的共模時(shí)間函數(shù)，即一維共性誤差.圖5對(duì)三種方法提取的一維共性誤差進(jìn)行了比較，可以看出三種方法提取的一維共性誤差隨時(shí)間變化趨勢(shì)十分相近，特別是PCA和KLE兩種方法的結(jié)果幾乎一致，水平(N、E)方向的一維共性誤差最大量級(jí)約為7 mm，振幅絕對(duì)值的平均量級(jí)約為1.5 mm；而垂直方向一維共性誤差的最大量級(jí)約為24 mm，振幅絕對(duì)值的平均量級(jí)約為5.2 mm.三種空間濾波方法振幅絕對(duì)值的平均量級(jí)差別很小(<0.1 mm)，PCA結(jié)果比KLE 和stacking的結(jié)果要稍微大一些.

圖4　各主分量特征值占總特征值的累計(jì)百分比，左為PCA分析的結(jié)果，右為KLE分析的結(jié)果Fig.4　Cumulative percentage of PC eigenvalues：(left)PCA result；(right)KLE result

圖5　分別使用stacking、PCA和KLE方法提取的共模時(shí)間函數(shù)(一維共性誤差).其中，灰色實(shí)線代表stacking提取結(jié)果，黑色圓點(diǎn)代表PCA提取結(jié)果，黑色”+”號(hào)代表KLE提取結(jié)果Fig.5　Comparison of the common mode time function extracted by stacking、PCA and KLE approach.The gray solid line represents stacking extraction results,black dots represent PCA extraction results,black “+” represents the KLE extraction result

圖6　stacking、PCA和KLE三種方法剔除共性誤差后的L1范數(shù)離散度序列比較Fig.6　Comparison of the daily L1 scatters of the regionally filtered time series by the stacking,PCA and KLE approach

4　共性誤差的剔除

4.1從殘差時(shí)間序列中剔除共性誤差

從殘差時(shí)間序列X中減去共性誤差ε，對(duì)于堆棧法，由于其共性誤差在空間均勻分布，公式如下：

(15)

對(duì)于PCA方法和KLE方法，由于其共性誤差并非均勻分布，公式如下：

(16)

為了統(tǒng)計(jì)和對(duì)比各方法共性誤差的剔除效果并減弱個(gè)別測(cè)站異常殘差值的影響，對(duì)剔除前后殘差的絕對(duì)值求和再除以測(cè)站數(shù)，得到單天L1范數(shù)離散度序列；對(duì)剔除前后的殘差平方和除以測(cè)站數(shù)再開平方，得到單天L2范數(shù)離散度序列，殘差的振幅定義為L(zhǎng)1范數(shù)的離散度，殘差的功率定義為L(zhǎng)2范數(shù)離散度的平方.圖6給出了不同方法剔除共性誤差前后殘差L1范數(shù)離散度(即殘差振幅)的比較，由圖6可知PCA和KLE方法濾波后的殘差振幅明顯小于Stacking方法的結(jié)果，說(shuō)明PCA和KLE方法對(duì)南極半島GPS數(shù)據(jù)的適應(yīng)性更好.N、E、U三個(gè)方向?yàn)V波后殘差L1范數(shù)離散度序列中冬季(南極)的振幅明顯大于夏季的，這可能與受溫度影響的誤差源有關(guān)，例如低層大氣環(huán)流、小尺度溫度擾動(dòng)、天線熱噪聲以及熱彈性張力(Prawirodirdjo et al.,2006)等.

表2比較了不同方法剔除共性誤差前后殘差的平均振幅、平均功率(對(duì)歷元求平均)和平均RMS(對(duì)測(cè)站求平均)，表3統(tǒng)計(jì)了殘差平均振幅、平均功率和平均RMS的下降幅度.由表2、3可知，三種空間濾波方法在濾波后殘差的平均振幅、平均功率以及平均RMS都有了明顯的下降，PCA方法的效果最好，其次是KLE方法.

4.2從站坐標(biāo)時(shí)間序列中剔除共性誤差

從站坐標(biāo)時(shí)間序列中剔除共性誤差的公式與從殘差時(shí)間序列中剔除共性誤差類似，只需將公式(15)、(16)中殘差時(shí)間序列替換為坐標(biāo)時(shí)間序列，共性誤差采用PCA方法提取的結(jié)果.圖7給出了剔除共性誤差前、后各測(cè)站坐標(biāo)時(shí)間序列，其中各測(cè)站縱軸相對(duì)坐標(biāo)是指相對(duì)于各測(cè)站的2010年1月1日初始坐標(biāo)的坐標(biāo).由圖7可以看出，濾波后的坐標(biāo)時(shí)間序列(黑點(diǎn))相比濾波前的序列更加收斂，其中U方向最為明顯，說(shuō)明使用PCA方法對(duì)坐標(biāo)時(shí)間序列進(jìn)行空間濾波能夠有效降低各站點(diǎn)坐標(biāo)時(shí)間序列的不確定性.

表2　不同方法濾波后殘差平均振幅、平均功率和平均RMS的比較

表3　不同方法濾波后殘差的平均振幅、平均功率和平均RMS的降幅比較

利用CATS軟件重新估計(jì)空間濾波前、后的坐標(biāo)時(shí)間序列的線性項(xiàng)和周期項(xiàng)參數(shù)，噪聲模型選擇白噪聲+閃爍噪聲+隨機(jī)漫步噪聲，各測(cè)站的參數(shù)估計(jì)結(jié)果在表4中給出.表4中各測(cè)站時(shí)間序列參數(shù)的誤差均有明顯的下降，表明坐標(biāo)時(shí)間序列剔除共性誤差后能夠大大減小估計(jì)線性項(xiàng)、周期項(xiàng)的誤差，從而提高坐標(biāo)時(shí)間序列的精度，使原來(lái)淹沒(méi)在誤差里的信號(hào)更加明顯.

圖7　PCA濾波前(灰點(diǎn))、后(黑色點(diǎn))的各測(cè)站相對(duì)坐標(biāo)時(shí)間序列Fig.7　The unfiltered and filtered position time series by PCA

站名濾波狀態(tài)速度(mm·a-1)周年項(xiàng)振幅(mm)周年項(xiàng)相位(rad)半周年振幅(mm)半周年項(xiàng)相位(rad)CAPF濾波前5.91±2.066.56±13.760.92±0.322.11±3.17-0.75±0.71濾波后6.15±0.875.56±4.940.80±0.162.40±1.49-0.16±0.27DUPT濾波前14.35±1.352.50±3.16-0.78±0.512.21±2.090.18±0.18濾波后14.56±0.492.65±1.18-0.58±0.161.56±0.500.32±0.20FONP濾波前17.31±1.742.94±5.260.09±0.583.57±4.600.10±0.36濾波后17.27±0.532.84±1.510.19±0.183.22±1.20-0.01±0.12HUGO濾波前6.56±1.825.74±9.660.62±0.301.70±2.06-1.30±0.73濾波后6.06±0.864.73±3.460.61±0.151.60±0.831.47±0.33PALM濾波前8.76±1.342.55±3.31-0.80±0.512.09±1.930.98±0.44濾波后8.85±0.402.53±0.93-0.74±0.141.89±0.490.83±0.14PALV濾波前9.48±1.352.47±3.19-0.99±0.532.29±2.131.03±0.40濾波后8.96±0.402.52±0.90-0.83±0.142.01±0.520.81±0.13ROBI濾波前9.65±1.961.45±2.740.76±1.292.70±3.630.61±0.50濾波后9.15±0.591.71±0.960.70±0.332.68±1.070.43±0.15VNAD濾波前7.43±1.352.39±3.07-0.14±0.542.32±2.150.10±0.39濾波后7.72±0.562.31±1.18-0.12±0.212.19±0.790.20±0.16

5　共性誤差的譜分析

共性誤差作為一種區(qū)域時(shí)空相關(guān)的誤差，并不是特指某一種誤差，它是某一區(qū)域大部分站點(diǎn)都存在的一種或多種誤差的集合，通過(guò)空間濾波，我們僅能得到這些誤差的總和，卻無(wú)法確定它由哪些誤差構(gòu)成以及哪些原因造成的.對(duì)時(shí)間序列進(jìn)行譜分析，可以對(duì)不同周期的信號(hào)進(jìn)行分離，有助于分析和解釋共性誤差的構(gòu)成.由于前面三種濾波方法中PCA濾波后的殘差離散度功率降幅最大，這里選用PCA方法得到的共性誤差進(jìn)行分析.

為了便于在時(shí)間域進(jìn)行分析，對(duì)該方法提取的二維共性誤差在空間域取平均(即在測(cè)站間取平均)，得到一維共性誤差，然后再進(jìn)行傅里葉變換，得到N、E、U三個(gè)分量的振幅譜.為了便于統(tǒng)計(jì)周期，以周期取代頻率作為X軸坐標(biāo)，單位為天，如圖8所示.我們發(fā)現(xiàn)N、E、U三個(gè)方向的最大振幅的周期分別出現(xiàn)在89.1、170.8和227.7天，對(duì)應(yīng)的最大振幅分別為0.56、0.46和1.71 mm·a-1，而且89.1天的周期同時(shí)存在于N、E、U三個(gè)方向，170.8天的周期同時(shí)存在于E、U兩個(gè)方向，227.7天的周期存在于N、U兩個(gè)方向，另外僅U方向存在周期50天振幅1.5 mm·a-1的信號(hào)，這些周期由于振幅足夠大，可以認(rèn)為是在共性誤差中的主要成分.其中部分周期可能是由于參考框架、大氣質(zhì)量負(fù)荷、未建模的衛(wèi)星軌道、EOP和交點(diǎn)周期(Wdowinski et al.,1997)等不準(zhǔn)確而引入的系統(tǒng)誤差，這些周期的具體成因目前尚無(wú)法解釋，還需要進(jìn)一步的研究和探討.對(duì)于U方向，還存在周期為9.4天(振幅0.69 mm·a-1)和13.66天(振幅0.90 mm·a-1)的短周期信號(hào)，考慮到海洋潮汐中存在9.1、9.6和13.7天(例如Mf分潮、半月太陰潮周期為13.7天)等周期，推測(cè)共性誤差組成成分中可能含與海潮相關(guān)的誤差源，該誤差源可能是由南極半島區(qū)域海洋潮汐未建模或模型不準(zhǔn)確引起的.

6　結(jié)論

本文使用了stacking、PCA和KLE三種空間濾波方法對(duì)南極半島區(qū)域GPS網(wǎng)絡(luò)坐標(biāo)時(shí)間序列進(jìn)行區(qū)域?yàn)V波，其中PCA和KLE方法使用了前兩階主分量.通過(guò)比較發(fā)現(xiàn)三種方法在南極半島區(qū)域均能有效提取共性誤差，但PCA方法的濾波效果明顯優(yōu)于stacking，且略微優(yōu)于KLE方法.使用PCA方法進(jìn)行空間濾波后，殘差時(shí)間序列的振幅、功率和RMS均有大幅下降，坐標(biāo)時(shí)間序列的線性項(xiàng)和周期項(xiàng)誤差也明顯減小，說(shuō)明進(jìn)行空間濾波能夠有助于提高坐標(biāo)時(shí)間序列參數(shù)估計(jì)的精度.

圖8　PCA方法提取共性誤差的幅值譜分析Fig.8　The amplitude spectrum analysis of the common mode error extracted by PCA method

通過(guò)對(duì)共性誤差進(jìn)行分析，本文總結(jié)了共性誤差的一些性質(zhì)：共性誤差是一種與區(qū)域空間尺度相關(guān)的誤差，隨著區(qū)域范圍的擴(kuò)大，原本的部分共性誤差因不適用整個(gè)區(qū)域會(huì)變?yōu)榫植啃?yīng)引起的特性誤差；共性誤差并非恒定不變，而是隨時(shí)間不斷發(fā)生變化；共性誤差并非均勻分布，因?yàn)橥灰蛩匾鸬恼`差對(duì)不同站點(diǎn)的影響大小并不相同.共性誤差的譜分析結(jié)果表明其U方向存在明顯的9.4、13.7天短周期信號(hào)，可能與南極半島區(qū)域海潮未建?；蚰Ｐ筒粶?zhǔn)確有關(guān).

References

Dong D,Fang P,Bock Y,et al.2006.Spatiotemporal filtering using principal component analysis and Karhunen-Loeve expansion approaches for regional GPS network analysis.Journal of Geophysical Research:Solid Earth (1978—2012),111(B3):B03405.E D C,Zhang S K.2006.Infrastructure and progress of international Antarctic geodetic reference framework.Journal of Geodesy and Geodynamics (in Chinese),26(2):104-108.

Hu S C,Wu J C,Sun Y F.2009.Comparison among three spatiotemporal filtering methods for regional GPS networks analysis.Journal of Geodesy and Geodynamics (in Chinese),29(3):95-99.

Jiang Z H,Zhang P,Mi J Z,et al.2010.Velocity estimation on the colored noise properties of CORS network in China based on the CGCS2000 frame.Acta Geodaetica et Cartographica Sinica (in Chinese),39(4):355-363.Johansson J M,Davis J L,Scherneck H G,et al.2002.Continuous GPS measurements of postglacial adjustment in Fennoscandia 1.Geodetic results.Journal of Geophysical Research:Solid Earth (1978—2012),107(B8):ETG 3-1-ETG 3-27.

Li W W,Ju X L,Shen Y Z,et al.2014.Vertical deformation analysis of Antarctic GNSS stations combined using GNSS and GRACE data.Chinese Journal of Polar Research (in Chinese),26(2):238-243.Nikolaidis R.2002.Observation of geodetic and seismic deformation with the Global Positioning System[Ph.D.thesis].Santiago:University of California.

Prawirodirdjo L,Ben-Zion Y,Bock Y.2006.Observation and modeling of thermoelastic strain in Southern California Integrated GPS Network daily position time series.Journal of Geophysical Research:Solid Earth (1978—2012),111(B2):B02408.Preisendorfer R W.1988.Principal Component Analysis in Meteorology and Oceanography.Amsterdam:Elsevier.Scherneck H G.1991.A parametrized solid earth tide model and ocean tide loading effects for global geodetic baseline measurements.Geophysical Journal International,106(3):677-694.

Serpelloni E,Faccenna C,Spada G,et al.2013.Vertical GPS ground motion rates in the Euro-Mediterranean region:New evidence of velocity gradients at different spatial scales along the Nubia-Eurasia plate boundary.Journal of Geophysical Research:Solid Earth,118(11):6003-6024.

Shen Y Z,Li W W,Xu G C,et al.2014.Spatiotemporal filtering of regional GNSS network′s position time series with missing data using principle component analysis.Journal of Geodesy,88(1):1-12.

Wdowinski S,Bock Y,Zhang J,et al.1997.Southern California Permanent GPS Geodetic Array:Spatial filtering of daily positions for estimating coseismic and postseismic displacements induced by the 1992 Landers earthquake.Journal of Geophysical Research:Solid Earth (1978—2012),102(B8):18057-18070.Wu J C,Sun Y F,Liu C G.2008.Extraction of common mode errors for continuous GPS networks and deformation analysis.Journal of Geodesy and Geodynamics (in Chinese),28(4):97-101.

附中文參考文獻(xiàn)

鄂棟臣,張勝凱.2006.國(guó)際南極大地參考框架的構(gòu)建與進(jìn)展.大地測(cè)量與地球動(dòng)力學(xué),26(2):104-108.

胡守超,伍吉倉(cāng),孫亞峰.2009.區(qū)域GPS網(wǎng)三種時(shí)空濾波方法的比較.大地測(cè)量與地球動(dòng)力學(xué),29(3):95-99.

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

李偉偉,鞠曉蕾,沈云中等.2014.聯(lián)合GNSS和GRACE數(shù)據(jù)分析南極跟蹤站的徑向形變.極地研究,26(2):238-243.

伍吉倉(cāng),孫亞峰,劉朝功.2008.連續(xù)GPS站坐標(biāo)序列共性誤差的提取與形變分析.大地測(cè)量與地球動(dòng)力學(xué),28(4):97-101.

(本文編輯胡素芳)

The coordinate time series analysis of continuous GPS stations in the Antarctic Peninsula with consideration of common mode error

MA Chao1,LI Fei1,2,ZHANG Sheng-Kai1*,LEI Jin-Tao1,E Dong-Chen1,HAO Wei-Feng1,ZHANG Qing-Chuan1

1 Chinese Antarctic Center of Surveying and Mapping,Wuhan University,Wuhan 430079,China 2 State Key Laboratory of Information Engineering in Surveying,Mapping and Remote Sensing, Wuhan University,Wuhan 430079,China

The spatially and temporally correlated common mode errors (CME)are the dominant error sources in daily GPS solutions in the regional CGPS network.Spatial filtering is an effective way to improve the precision of coordinate time series for regional CGPS networks by reducing these errors.The data of 11 GPS stations from 2010 to 2014 in Antarctic Peninsula are processed with GAMIT/GLOBK10.5 software.The set of individually estimated daily positions then make up the position time series for every stations.Three filtering approaches including stacking,principal component analysis (PCA)and Karhunen-Loeve expansion (KLE)were applied to daily coordinate time series of Antarctic Peninsula GPS network from 2010 to 2014.The first two principal components in PCA and KLE were used as common modes.The filtering results show that all of the three methods can effectively extract the common mode errors in the Antarctic Peninsula,but the results of PCA are much better than those of stacking,and slightly better than those of KLE.We demonstrate that spatial filtering can effectively reduce the amplitude and power of the residual time series and effectively reduce the errors of linear term and periodic term in coordinate time series,so as to improve the accuracy of their estimations.The spectrum results of common mode errors from PCA analysis show that 9.4-day and 13.7-day period signals in the vertical component might be related to the ocean tide model used in the analysis.

Antarctic Peninsula;GPS;Time series;Common mode errors

馬超,李斐,張勝凱等.2016.顧及共性誤差的南極半島地區(qū)連續(xù)GPS站坐標(biāo)時(shí)間序列分析.地球物理學(xué)報(bào)，59(8):2783-2795,

10.6038/cjg20160806.

Ma C,Li F,Zhang S K,et al.2016.The coordinate time series analysis of continuous GPS stations in the Antarctic Peninsula with consideration of common mode error.Chinese J.Geophys.(in Chinese),59(8):2783-2795,doi:10.6038/cjg20160806.

國(guó)家自然科學(xué)基金重點(diǎn)項(xiàng)目(41531069)，國(guó)家重大科學(xué)研究973計(jì)劃(2012CB957701，2013CBA01804)，國(guó)家自然科學(xué)基金項(xiàng)目(41174019,41176173)，南北極環(huán)境綜合考察及資源潛力評(píng)估專項(xiàng)(CHINARE2015)聯(lián)合資助.

馬超,男，1990年生，山東滕州人，現(xiàn)為武漢大學(xué)碩士研究生，主要從事極地大地測(cè)量學(xué)的研究.E-mail：mac1207@whu.edu.cn

張勝凱，男，1977年生，山東東明人，博士，副教授，研究方向?yàn)槟蠘OGNSS應(yīng)用研究.E-mail：zskai@whu.edu.cn

10.6038/cjg20160806

P228

2015-12-02，2016-01-14收修定稿