GPS坐標時間序列噪聲協方差陣的構造及其方差估計

金航航，馬海翔，馬俊

(1.東陽市經濟開發區規劃國土建設局，浙江 東陽　322100；　2. 中鐵十一局集團第二工程有限公司，湖北 十堰　442000；3.武漢大學，湖北 武漢　430079)

1　引　言

GPS坐標時間序列因其可以為地殼形變監測以及一些地球物理現象(如板塊運動、應變積累、火山變形、冰后回彈、沉降以及海平面變化)的解釋提供寶貴的基礎數據，因此被廣泛應用于大地測量學及地球動力學研究。由于受多種因素(如多路徑效應、衛星鐘差、地理環境以及大氣延遲)的影響，GPS坐標時間序列中包含了信號和誤差(噪聲)兩部分。噪聲嚴重影響了坐標時間序列的穩定以及測站運動參數(包括速度、周期振幅等)估值的精度及其不確定度。針對GPS坐標時間序列建立適當的噪聲模型，可實現測量信號和噪聲的有效分離，獲得準確的測站運動趨勢，對于合理解釋地殼形變特征，維護和更新地球參考框架以及探求和揭示大地構造變形運動過程具有重要的意義[1]。

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

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

由于GPS坐標時間序列中噪聲及其功率譜的特性與其在分數差分中的表現形式有關，因此本文簡述了分形噪聲理論及其功率譜的特點，在此基礎上詳細推導了有色噪聲的協方差表達式，論述了利用最小范數二次無偏估計法估計噪聲方差，最后通過實測數據展示了該方法的計算過程。

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

Hosking認為分形噪聲是白噪聲的(0.5-H)階差分。其中H是用來度量時間序列的相關性和趨勢強度的赫斯特指數，其取值范圍為(0，1)。分形噪聲的特性與H的取值范圍有關。當H=0.5時為正態的高斯白噪聲。H≠0.5時，對應非正態過程，即分形噪聲。當0.5

xt=-dat

(1)

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

(2)

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

(3)

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

(4)

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

(5)

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

(6)

其中{xt}的特性與d的取值范圍有關：當0

3　GPS坐標時間序列噪聲估計

3.1　噪聲的譜指數

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

(7)

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

(8)

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

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

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

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

(9)

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

vi=σwαi+σKβi

(10)

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

y=Ax+ε

(11)

(12)

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

(13)

進一步寫成矩陣形式，如下式所示：

(14)

其中n為時間序列的長度。因此其協因數矩陣可以表示為：

QK=TCWTT

(15)

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

(16)

其中：ψ0=1，

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

(17)

skl=tr(CQlCQk)

(18)

wk=yTCQkCyl，k=1，2

(19)

4　算例說明

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

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

噪聲方差估值　　　　　　　　　　表1

5　結　語

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

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

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

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

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

綜上所述，時間序列相關噪聲分析從分數差分理論、統計學以及噪聲的潛在來源出發，研究有色噪聲的類型、強度、及與其他噪聲之間的關系，對于GPS坐標時間序列的深入研究與應用具有指導意義。