金航航，馬海翔，馬俊

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

1　引　言

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

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

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

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

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

Hosking認(rèn)為分形噪聲是白噪聲的(0.5-H)階差分。其中H是用來度量時(shí)間序列的相關(guān)性和趨勢(shì)強(qiáng)度的赫斯特指數(shù)，其取值范圍為(0，1)。分形噪聲的特性與H的取值范圍有關(guān)。當(dāng)H=0.5時(shí)為正態(tài)的高斯白噪聲。H≠0.5時(shí)，對(duì)應(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]，將其用滑動(dòng)平均模型表示如下：

(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)時(shí)間序列噪聲估計(jì)

3.1　噪聲的譜指數(shù)

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

(7)

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

(8)

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

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

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

GPS單站、單分量的運(yùn)動(dòng)模型如下所示：

(9)

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

vi=σwαi+σKβi

(10)

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

y=Ax+ε

(11)

(12)

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

(13)

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

(14)

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

QK=TCWTT

(15)

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

(16)

其中：ψ0=1，

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

(17)

skl=tr(CQlCQk)

(18)

wk=yTCQkCyl，k=1，2

(19)

4　算例說明

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

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

噪聲方差估值　　　　　　　　　　表1

5　結(jié)　語

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

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

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

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

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

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