GPS坐標(biāo)時間序列噪聲協(xié)方差陣的構(gòu)造及其方差估計

金航航，馬海翔，馬俊

(1.東陽市經(jīng)濟(jì)開發(fā)區(qū)規(guī)劃國土建設(shè)局，浙江 東陽　322100；　2. 中鐵十一局集團(tuán)第二工程有限公司，湖北 十堰　442000；3.武漢大學(xué)，湖北 武漢　430079)

1　引　言

GPS坐標(biāo)時間序列因其可以為地殼形變監(jiān)測以及一些地球物理現(xiàn)象(如板塊運動、應(yīng)變積累、火山變形、冰后回彈、沉降以及海平面變化)的解釋提供寶貴的基礎(chǔ)數(shù)據(jù)，因此被廣泛應(yīng)用于大地測量學(xué)及地球動力學(xué)研究。由于受多種因素(如多路徑效應(yīng)、衛(wèi)星鐘差、地理環(huán)境以及大氣延遲)的影響，GPS坐標(biāo)時間序列中包含了信號和誤差(噪聲)兩部分。噪聲嚴(yán)重影響了坐標(biāo)時間序列的穩(wěn)定以及測站運動參數(shù)(包括速度、周期振幅等)估值的精度及其不確定度。針對GPS坐標(biāo)時間序列建立適當(dāng)?shù)脑肼暷Ｐ停蓪崿F(xiàn)測量信號和噪聲的有效分離，獲得準(zhǔn)確的測站運動趨勢，對于合理解釋地殼形變特征，維護(hù)和更新地球參考框架以及探求和揭示大地構(gòu)造變形運動過程具有重要的意義[1]。

由于GPS坐標(biāo)時間序列中的有色噪聲統(tǒng)計特性各不相同，因此，分析測站坐標(biāo)時間序列時，需根據(jù)時間序列的噪聲特性構(gòu)造相應(yīng)的協(xié)方差陣。Agnew[2]直接給出了冪律噪聲的功率譜密度公式，Williams[3]給出了有色協(xié)方差矩陣的廣義形式，Zhang[4]等人提出了WN+FN噪聲的混合模型，并給出了閃爍噪聲的協(xié)方差矩陣形式。在此基礎(chǔ)上一些學(xué)者提出了噪聲方差的估計方法，如Williams結(jié)合極大似然估計與單純形搜索算法編寫了CTAS軟件[5]，被廣泛地應(yīng)用于GPS單分量時間序列的分析。Amiri-Simkooei提出利用最小二乘方差分量估計法(LS-VCE)分析GPS坐標(biāo)時間序列噪聲，該方法可以較為精確地估計出噪聲的方差[6]。

對于常見的白噪聲和隨機(jī)漫步噪聲，其協(xié)方差矩陣已有明確的表達(dá)式，而對于其他譜指數(shù)對應(yīng)的有色噪聲，尚無文獻(xiàn)給出其協(xié)方差矩陣表達(dá)式詳細(xì)的推導(dǎo)過程，均直接給出了不同有色噪聲模型的協(xié)方差陣。這不利于理解不同噪聲在時間域上的相關(guān)特性及其功率譜密度在頻域上的特性，影響對GPS坐標(biāo)時間序列所反映的測站運動特征的深入研究。此外，在利用極大似然估計法以及最小二乘方差估計法估計GPS坐標(biāo)時間序列的噪聲方差時，需要多次的迭代與搜索過程，因此運算需要較長的時間，尤其是對于跨度較長的時間序列。

由于GPS坐標(biāo)時間序列中噪聲及其功率譜的特性與其在分?jǐn)?shù)差分中的表現(xiàn)形式有關(guān)，因此本文簡述了分形噪聲理論及其功率譜的特點，在此基礎(chǔ)上詳細(xì)推導(dǎo)了有色噪聲的協(xié)方差表達(dá)式，論述了利用最小范數(shù)二次無偏估計法估計噪聲方差，最后通過實測數(shù)據(jù)展示了該方法的計算過程。

2　分形噪聲及其功率譜密度

Hosking認(rèn)為分形噪聲是白噪聲的(0.5-H)階差分。其中H是用來度量時間序列的相關(guān)性和趨勢強(qiáng)度的赫斯特指數(shù)，其取值范圍為(0，1)。分形噪聲的特性與H的取值范圍有關(guān)。當(dāng)H=0.5時為正態(tài)的高斯白噪聲。H≠0.5時，對應(yīng)非正態(tài)過程，即分形噪聲。當(dāng)0.5

xt=-dat

(1)

上式中{xt}為隨機(jī)過程，d=(1-B)d為差分算子，B為后向移位算子(Bxt=xt-1)，{at}為白噪聲，其均值與方差分別0和差分算子展開式如下：

(2)

由上式可知，ARIMA(0，d，0)模型為基本的分?jǐn)?shù)差分模型，又被稱為分?jǐn)?shù)差分白噪聲[6]，將其用滑動平均模型表示如下：

(3)

根據(jù)差分算子式(2)可得：

(4)

ARIMA(p，d，q)模型的功率譜可用下式表示：

(5)

S(ω)=[2sin(0.5ω)]-2d

(6)

其中{xt}的特性與d的取值范圍有關(guān)：當(dāng)0

3　GPS坐標(biāo)時間序列噪聲估計

3.1　噪聲的譜指數(shù)

許多地球物理信號的統(tǒng)計模型可以用冪律過程來描述，其功率譜可以用如下式子來表示[2]：

(7)

其中P0和f0為常數(shù)，f為空間或時間上的頻率，K為譜指數(shù)，其取值范圍為[-3，1]。具有不同譜指數(shù)的噪聲其特性也不相同。-1

(8)

由此可知，確定噪聲類型及其特性的關(guān)鍵是獲取噪聲譜指數(shù)，具體方法讀者可自行查閱相關(guān)文獻(xiàn)。

3.2　利用最小范數(shù)二次無偏估計法估計噪聲方差

為了獲得測站準(zhǔn)確的運動趨勢，減小速度估值的不確定度，就必須采取有效的措施減小噪聲對數(shù)據(jù)分析的影響。確定譜指數(shù)后可利用極大似然估計法或者方差-協(xié)方差分量估計法估計噪聲方差然后確定其對速度不確定度的影響。受篇幅所限，本文主要介紹利用最小范數(shù)二次無偏估計法分析GPS坐標(biāo)時間序列噪聲。

GPS單站、單分量的運動模型如下所示：

(9)

其中，y(ti)為GPS單站單分量的坐標(biāo)時間序列，ti為時間，x(1)為初始位置，x(2)為速率，x(2K-1)和x(2K)是調(diào)和函數(shù)的系數(shù)，q-1為時間序列中的周期項數(shù)，vi為白噪聲α與有色噪聲β的線性組合，其表現(xiàn)形式如下：

vi=σwαi+σKβi

(10)

上式中σw和σK分別為白噪聲和有色噪聲的振幅。式(9)的矩陣形式以及隨機(jī)模型如下：

y=Ax+ε

(11)

(12)

上式中，A是設(shè)計矩陣，ε是誤差，D(y)為協(xié)方差矩陣，I和QK分別為白噪聲和有色噪聲的協(xié)因數(shù)矩陣。當(dāng)前研究表明GPS測站的非線性運動中包含周年運動和半周年運動，因此A矩陣的表現(xiàn)形式如下：

(13)

進(jìn)一步寫成矩陣形式，如下式所示：

(14)

其中n為時間序列的長度。因此其協(xié)因數(shù)矩陣可以表示為：

QK=TCWTT

(15)

由于白噪聲協(xié)方差矩陣Cw=I，因此Qk=TTT。-3

(16)

其中：ψ0=1，

在精確地估計出各類噪聲方差之前，首先要合理地確定觀測值的權(quán)陣。由于白噪聲和有色噪聲的方差(即驗前方差)未知，即使確定時間序列的最佳噪聲組合，也無法精確地定權(quán)。最小范數(shù)二次無偏估計(MINQUE)適用于估計同類觀測值中不同因素的方差分量[8]，因此可用該方法估計GPS坐標(biāo)時間序列中白噪聲以及有色噪聲的方差分量，進(jìn)而準(zhǔn)確地定權(quán)。限于篇幅本文直接給出最小范數(shù)二次無偏估計的計算公式如下：

(17)

skl=tr(CQlCQk)

(18)

wk=yTCQkCyl，k=1，2

(19)

4　算例說明

本文采用位于美國南加州的AZU1站以及ROCK站的垂直方向的時間序列。數(shù)據(jù)由JPL提供，時間跨度為2010年1月～2016年9月，利用Matlab編程計算噪聲方差。在估計噪聲之間，我們探測并剔除了時間序列中的粗差，根據(jù)JPL給出的階躍估值對剔除粗差的時間序列進(jìn)行了改正。

當(dāng)前國內(nèi)外研究表明，GPS坐標(biāo)時間序列中的季節(jié)性信號主要為周年以及半周年信號[9]，且最佳噪聲模型為白噪聲+閃爍噪聲組合[10]。因此，函數(shù)模型中的周期信號僅包括周年以及半周年信號，取有色噪聲的譜指數(shù)K=-1。根據(jù)式(15)、式(16)構(gòu)造閃爍噪聲的協(xié)因數(shù)矩陣，再結(jié)合設(shè)計矩陣A構(gòu)成M以及C矩陣；然后根據(jù)式(18)和式(19)分別構(gòu)成S矩陣以及W矩陣，最后帶入式(16)中獲得白噪聲以及閃爍噪聲的方差估值，具體結(jié)果如表1所示。表1中WN和FN分別表示白噪聲以及閃爍噪聲方差的估值。此外，表1中還列出了利用LS-VCE法獲得的噪聲方差估值。由表1可以看出，MINQUE法的估計結(jié)果與LS-VCE相差較小，且閃爍噪聲的強(qiáng)度明顯大于白噪聲。因此在分析GPS坐標(biāo)時間序列時，應(yīng)采取有效措施降低噪聲，尤其是有色噪聲對測站運動速度估值的影響。

噪聲方差估值　　　　　　　　　　表1

5　結(jié)　語

本文根據(jù)時間序列分形噪聲理論及其功率譜密度討論了冪律噪聲的功率譜密度的特性，推導(dǎo)了GPS坐標(biāo)時間序列中冪律噪聲的協(xié)方差矩陣構(gòu)造過程。在此基礎(chǔ)上論述了利用最小二乘方差分量估計分析GPS坐標(biāo)時間序列噪聲的方法。由以上分析可知：

(1)由冪律噪聲協(xié)方差陣構(gòu)造過程可知，GPS坐標(biāo)時間序列中的冪律噪聲可以看作是將某一卷積函數(shù)與一白噪聲卷積無限求和形式截斷到有限和后得到的，時間序列跨度越長，相關(guān)冪律噪聲的協(xié)方差矩陣越準(zhǔn)確，因此在研究GPS測站運動特征時應(yīng)選取跨度較長的時間序列。

(2)確定GPS坐標(biāo)時間序列噪聲的譜指數(shù)后，可利用MINQUE法估計噪聲的方差。尤其是分析較多GPS測站以及跨度較長的時間序列時，有助于提高分析效率，節(jié)省運算時間。此外還應(yīng)采取有效措施降低有色噪聲的強(qiáng)度，減小其對測站運動速度估計結(jié)果的影響。

(3)對于分布范圍較大的GPS測站，由于其所在地區(qū)的地球物理環(huán)境不同，對于不同的GPS測站，其坐標(biāo)時間序列中的有色噪聲也不盡相同。應(yīng)估計在不同噪聲模型組合下的GPS坐標(biāo)時間序列的極大似然估計值，確定最佳噪聲組合模型及其噪聲方差。

(4)以往坐標(biāo)時間序列噪聲模型的建立普遍基于單一分量，忽略了水平和高程分量之間的相關(guān)性，致使構(gòu)造信號不能得到充分的解釋。因此，在以后的研究中，應(yīng)量化水平及高程方向的相互滲透作用程度，構(gòu)建精確的三維噪聲模型，獲取受噪聲影響較小的GNSS基準(zhǔn)站坐標(biāo)時間序列。

綜上所述，時間序列相關(guān)噪聲分析從分?jǐn)?shù)差分理論、統(tǒng)計學(xué)以及噪聲的潛在來源出發(fā)，研究有色噪聲的類型、強(qiáng)度、及與其他噪聲之間的關(guān)系，對于GPS坐標(biāo)時間序列的深入研究與應(yīng)用具有指導(dǎo)意義。