元亞菲，張以文，孔慶珠，李朋朋，李志遠(yuǎn)，張雨陽

（1.中國礦業(yè)大學(xué) 環(huán)境與測繪學(xué)院，江蘇 徐州 221116；2.水發(fā)規(guī)劃設(shè)計有限公司，山東 濟(jì)南 250100）

0 引 言

近幾十年來，高精度、高時間分辨率的GNSS區(qū)域網(wǎng)建設(shè)取得了長足發(fā)展，為研究板塊運(yùn)動機(jī)制、建立和維護(hù)地球參考框架、揭示大地構(gòu)造形變過程等提供了重要的數(shù)據(jù)支撐。大量研究結(jié)果表明，由于GNSS站易受外部環(huán)境、自身觀測技術(shù)和地球構(gòu)造運(yùn)動等影響，GNSS站坐標(biāo)時間序列存在明顯的有色噪聲［1-2］，且不同方向上的有色噪聲存在一定相關(guān)性［3］。因此，研究GNSS站坐標(biāo)時間序列的最優(yōu)噪聲模型，確定不同方向噪聲的耦合關(guān)系，可以獲取更為準(zhǔn)確的運(yùn)動特征，進(jìn)而更加合理、全面地進(jìn)行地球物理現(xiàn)象研究。針對GNSS站坐標(biāo)時間序列的最優(yōu)噪聲模型和運(yùn)動特征，國內(nèi)外眾多學(xué)者進(jìn)行了深入分析，得到了一系列有益成果。如：A.Mao等［4］對全球23個GPS站三年的坐標(biāo)序列進(jìn)行分析，認(rèn)為白噪聲（white noise，WN）和閃爍噪聲（flicker noise，F(xiàn)N）的模型組合可以作為全球GPS站噪聲模型的最佳選擇；S.D.P.Williams等［5］利用極大似然估計法（MLE）對全球954個GPS站坐標(biāo)序列進(jìn)行分析，證明了WN+FN能更好地表示坐標(biāo)序列的噪聲模型；明鋒等［6］利用STL法提取了中國區(qū)域IGS站的長時間變化趨勢，并利用MLE法估計了各站點(diǎn)的最優(yōu)噪聲模型；馬俊等［7］采用最小范數(shù)二次無偏估計法（MINQUE）分析了中國區(qū)域20個IGS站的噪聲分布情況；He X等［8］采用MLE法估計了全球110個測站的最優(yōu)噪聲模型，并分析低頻噪聲對趨勢估值的長期影響。

已有GNSS站坐標(biāo)時間序列研究中，噪聲模型估計方法主要為極大似然估計法和方差協(xié)方差分量估計法（VCE）。極大似然估計法是噪聲估計最為常見的方法，該方法根據(jù)預(yù)先設(shè)置的噪聲模型和相應(yīng)殘差建立似然函數(shù)，并根據(jù)不同噪聲模型的極大似然值確定最優(yōu)的噪聲類型和噪聲強(qiáng)度［9-12］。但該方法計算復(fù)雜［13］，且無法直接給出參數(shù)的不確定度，只能通過相應(yīng)的簡化公式計算求得［6］，同時無法對噪聲強(qiáng)度的估值進(jìn)行質(zhì)量評價。方差協(xié)方差分量估計法以最小二乘估計為基礎(chǔ)，可以直接給出參數(shù)的不確定度，并且可以根據(jù)方差估計式的函數(shù)關(guān)系進(jìn)行噪聲強(qiáng)度估值的精度評價［14-16］。當(dāng)前噪聲模型的研究大多采用上述兩種方法，然而鮮有這兩種方法適用性的對比研究。同時，隨著GNSS觀測值不斷累積，時間跨度的增大必然進(jìn)一步提高最終的估計精度，同時也會降低計算效率，因此，最優(yōu)噪聲模型的估計還需同時顧及估計精度和估計效率。

基于此，本文針對極大似然估計法和方差分量估計法在GNSS站坐標(biāo)時間序列噪聲模型估計中的估計效率和估計精度，利用仿真數(shù)據(jù)對比極大似然估計法和兩種方差分量估計法（最小二乘方差分量估計法（LS-VCE）和最小范數(shù)二次無偏估計法（MINQUE））對不同噪聲強(qiáng)度、不同時間跨度和時間序列的適應(yīng)性。在此基礎(chǔ)上，利用本文確定的最優(yōu)方法對中國區(qū)域15個GNSS站的坐標(biāo)時間序列進(jìn)行噪聲方差估計，并利用線性回歸法分析不同方向噪聲間的相關(guān)性和交互關(guān)系。

1 GNSS站坐標(biāo)時間序列噪聲方法估計

GNSS站坐標(biāo)時間序列分為趨勢項、季節(jié)項和噪聲項，其觀測模型可表示為

式中：y t()i為ti歷元的坐標(biāo)值；a，b分別為站點(diǎn)坐標(biāo)的初始位置和線性速度；q為周期項的個數(shù)；fk為第k個周期項的頻率；ck，sk為相應(yīng)的周期項系數(shù)，用于表示各個周期項的振幅；εi為包含多路徑效應(yīng)、電離層高階項等多種未模型化誤差在內(nèi)的誤差項。一般情況下，GNSS站坐標(biāo)時間序列取年周期項和半年周期項，fk取值為1或2。

將式（1）寫成矩陣形式，可得GNSS站坐標(biāo)時間序列的函數(shù)模型和隨機(jī)模型，即

對于中國區(qū)域大多數(shù)GNSS站而言，最佳的坐標(biāo)時間序列隨機(jī)模型為白噪聲和閃爍噪聲［17］，即QWN=QFN。根據(jù)分?jǐn)?shù)階差分方程的性質(zhì)，閃爍噪聲可以表示為特定卷積函數(shù)與高斯白噪聲的卷積［18-19］，則閃爍噪聲和協(xié)因數(shù)矩陣為

式中，H為白噪聲到閃爍噪聲的卷積矩陣，其具體表達(dá)式為

根據(jù)式（2）～（6）構(gòu)建GNSS站坐標(biāo)時間序列的函數(shù)模型和隨機(jī)模型，即利用VCE法估計白噪聲和閃爍噪聲的方差和。LS-VCE法噪聲方差估計方程［20］為

其中，Qj為白噪聲或閃爍噪聲的協(xié)因數(shù)陣，Q0為噪聲強(qiáng)度已知的協(xié)因數(shù)陣，為殘差，且～

為了更加直觀地評定對噪聲方差估值的精度，可利用誤差傳播率求得LS-VCE法噪聲方差估值的協(xié)方差矩陣

①您是否忘記過服藥；②過去2周內(nèi)，有多少天忘記服藥；③治療期間，出現(xiàn)黃疸等癥狀或出現(xiàn)其他癥狀，是否告知醫(yī)師，自行減少用藥劑量，或停止用藥；④您外出時，是否忘記隨身攜帶藥物；⑤昨天是否用藥；⑥當(dāng)自覺高低血糖等癥狀是否好轉(zhuǎn)或消失，是否停止過用藥；⑦您是否覺得堅持治療是否存在困難。7～8分為優(yōu)良，6分以下為依從性差。

Q0=0時，即噪聲模型中不含有噪聲強(qiáng)度已知的部分時，式（7）等號右邊Wi可化簡為

矩陣二次型的協(xié)方差矩陣公式

可得

因此，DW可以表示為

2 仿真數(shù)據(jù)對比

極大似然估計法是GNSS站坐標(biāo)時間序列噪聲模型估計中最常用的方法，該方法利用殘差概率密度分布構(gòu)建似然函數(shù)，通過似然函數(shù)的極大值確定最佳噪聲分量強(qiáng)度，但極大似然估計法僅能得到噪聲分量的估計值，無法給出噪聲分量的估計精度，因此該方法的正確性需要進(jìn)一步檢驗［6］。此外，在估計噪聲模型時，時間序列長度是影響時間序列分析結(jié)果的一個重要影響因素，時間序列越長，參數(shù)估計的精度越高，但帶來的計算時間也急劇增加，因此最佳時間跨度需要考慮計算效率和計算準(zhǔn)確度兩方面的要求?；诖?，利用仿真數(shù)據(jù)對比不同時間長度、不同噪聲強(qiáng)度下MLE法、LS-VCE法和MINQUE法的估計準(zhǔn)確度和計算效率。

為使模擬數(shù)據(jù)更接近真實(shí)的時間序列，采用噴氣推進(jìn)實(shí)驗所（Jet Propulsion Laboratory，JPL）提供的BJFS站北方向的運(yùn)動參數(shù)，站點(diǎn)運(yùn)動參數(shù)如表1所示。此外，在噪聲模擬時，由于各分析中心沒有給出最終的噪聲模型，因此本文基于全球IGS站坐標(biāo)時間序列的分析結(jié)果［10］選擇白噪聲和閃爍噪聲作為噪聲模型，并設(shè)計如下方案。

表1 BJFS站點(diǎn)運(yùn)動參數(shù)Tab.1 Motion parameters of BJFS stations

（1）方案A。白噪聲方差取0.25 mm2，閃爍噪聲方差取1 mm2/a0.5。

（2）方案B。白噪聲方差取2.25 mm2，閃爍噪聲方差取9.75 mm2/a0.5。

（3）方案C。白噪聲方差取6.25 mm2，閃爍噪聲方差取25 mm2/a0.5。

采用上述方案生成7，9，11，13，15，17和19年的模擬時間序列，并分別采用MLE法、LS-VCE法和MINQUE法進(jìn)行噪聲方差估計，為使結(jié)果更具統(tǒng)計意義，每個方案模擬20次，并統(tǒng)計不同時間長度、不同噪聲強(qiáng)度下噪聲方差估計值和估計中誤差（僅限LS-VCE和MINQUE法）的平均值，統(tǒng)計結(jié)果如表2所示。

由表2可知，從整體上看，無論采用哪種方法，求得的白噪聲方差估計值相較于閃爍噪聲方差估計值更接近理論真值，白噪聲的平均偏差約為3%，而閃爍噪聲方差的平均偏差約為白噪聲平均偏差的3倍，該值也接近于預(yù)先設(shè)置的方差比值，因此閃爍噪聲估值較差的原因可能與閃爍噪聲強(qiáng)度較大而白噪聲強(qiáng)度較小有關(guān)。此外，3種方法得到的噪聲方差估值基本相等，噪聲強(qiáng)度較小時，MLE法的估計結(jié)果與真值相差較大，平均偏差約為2.6%，兩種方差分量估計方法所得結(jié)果幾乎相等，驗證了LS-VCE法和MINQUE法在方差分量估計方面的等價性［20］，總體而言，相較于MLE法偏差更小，其平均偏差約為1.9%。

表2 仿真時間序列噪聲方差統(tǒng)計結(jié)果Tab.2 Estimation resalts of the noise models using the simulated time series

從噪聲方差估值偏差的變化趨勢看，當(dāng)時間長度較短時，噪聲方差的估值偏差與噪聲方差呈現(xiàn)出一定的負(fù)相關(guān)性，即隨著噪聲強(qiáng)度提高，估計偏差逐漸減小，估計精度逐漸增加；當(dāng)時間序列超過一定長度時，噪聲方差的估值偏差與噪聲方差呈現(xiàn)出正相關(guān)性，即噪聲方差越大，估計偏差越大，估計精度越低。從估計中誤差角度看，噪聲方差的估計精度隨噪聲方差增加而降低。時間長度較短時，估計偏差表現(xiàn)出與估計中誤差完全相反的變化趨勢，這種不一致性說明在時間長度不足時估計結(jié)果與估計精度容易出現(xiàn)異常，因此，對GNSS站坐標(biāo)時間序列的噪聲模型分析必須保證足夠的時間長度。

噪聲強(qiáng)度相同時，噪聲方差的估計精度與時間長度呈現(xiàn)出明顯的正相關(guān)性，但時間長度超過一定值后，噪聲方差的估計精度并沒有隨時間長度的增加顯著提高，如采用19 a時間序列的噪聲強(qiáng)度估計值和估計精度相較于采用15 a時間序列的并沒有顯著提高。從理論上分析，時間長度越大，呈正態(tài)分布的白噪聲對時間序列的影響越小，越有利于有色噪聲特性的表達(dá)，噪聲強(qiáng)度的估計值也會更加精準(zhǔn)。但是時間長度對估計精度的提高有限，達(dá)到閾值后繼續(xù)提高時間長度只會增加計算時間成本，因此盲目提高時間長度并不會進(jìn)一步提高估計精度。

此外，比較3種方法的計算時間，MINQUE法與LS-VCE法單次迭代的計算復(fù)雜度近乎相等［16］，且在相同收斂條件下兩者的迭代次數(shù)沒有明顯區(qū)別，故MINQUE法與LS-VCE法的計算時間近乎相等。MLE法受限于噪聲矩陣的條件數(shù)，即使采用優(yōu)化算法，其計算效率仍然低于VCE法的，特別是在時間長度較長時，MLE法的計算時間為LS-VCE法的數(shù)十倍。因此，方差分量估計法更適用于噪聲模型估計，為便于比較，選擇LSVCE法作為后續(xù)噪聲模型的估計方法。同時，綜合考慮噪聲方差的估計精度和計算效率，采用15 a的GNSS站坐標(biāo)時間序列作為最佳時間長度。

3 GNSS站坐標(biāo)時間序列噪聲模型 估計

3.1 數(shù)據(jù)預(yù)處理

本文采用的GNSS站坐標(biāo)時間序列來源于中國地殼監(jiān)測網(wǎng)（CMONOC）提供的ITRF2014框架下的GNSS站單天坐標(biāo)數(shù)據(jù)。該數(shù)據(jù)在解算過程中采用雙差觀測值構(gòu)建觀測方程，利用GAMIT10.4進(jìn)行單日松弛解，解算出各站點(diǎn)的三維坐標(biāo)、速度和相關(guān)參數(shù)，然后將斯克里普斯軌道和永久陣列中心（Scripps Orbit and Permanet Array Center，SOPAC）提供的IGS全球解與松弛解進(jìn)行聯(lián)合處理，得到ITRF2014框架下的站坐標(biāo)。CMONOC提供原始時間序列（Raw）和去除粗差階躍項等的時間序列（Detrend）數(shù)據(jù)，為了盡可能保留數(shù)據(jù)的原始信息，選用原始坐標(biāo)時間序列作為原始數(shù)據(jù)。

為了更加準(zhǔn)確地估計站點(diǎn)坐標(biāo)時間序列的噪聲模型和三維運(yùn)動特征，綜合考慮計算效率和計算準(zhǔn)確度，結(jié)合仿真數(shù)據(jù)結(jié)果，選用具有15 a以上觀測值的GNSS站。同時為避免時間長度不一致導(dǎo)致地球物理因素對分析結(jié)果產(chǎn)生影響，選用2000—2014年共15 a的觀測值數(shù)據(jù)作為本文的數(shù)據(jù)源。

針對原始坐標(biāo)時間序列中采樣率不一致的問題，為減少數(shù)據(jù)缺失和數(shù)據(jù)插值策略影響結(jié)果，以連續(xù)缺失180 d和總?cè)笔?00 d為閾值進(jìn)行篩選，共選出15個具有相同接收機(jī)天線的高質(zhì)量GNSS觀測站。利用最小二乘法和自回歸模型對數(shù)據(jù)缺失點(diǎn)進(jìn)行插值［21-22］，同時采用IQR準(zhǔn)則對原始坐標(biāo)時間序列進(jìn)行粗差探測與剔除。為保證處理后坐標(biāo)時間序列的可靠性，在IQR粗差探測過程中多次迭代，直到無法探測出新的粗差為止。此外，為了消除地球質(zhì)量荷載等因素引起的共模誤差，采用相關(guān)加權(quán)疊加濾波法對上述坐標(biāo)時間序列進(jìn)行處理。以BJFS站為例，經(jīng)數(shù)據(jù)預(yù)處理后的N（北）方向，E（東）方向和U（豎直）方向的坐標(biāo)時間序列如圖1所示。

圖1 BJFS站預(yù)處理后的坐標(biāo)時間序列圖Fig.1 Coordinate time series diagram at BJFS after data processing

在上述基礎(chǔ)上，利用Lomb-Scarge法對所選站點(diǎn)的坐標(biāo)時間序列進(jìn)行頻譜分析，所得的疊加功率譜分析結(jié)果如圖2所示。結(jié)果表明，年周期和半年周期為中國區(qū)域所選站點(diǎn)的主要周期項，此外在3cpy和4cpy（cycles per year）附近出現(xiàn)一定的頻率譜峰，但不能體現(xiàn)在所有方向上，因此，選擇年周期項和半年周期項作為最終的周期項。

圖2 疊加功率譜Fig.2 Superimposed power spectrums

3.2 噪聲模型估計結(jié)果

相較于VCE法而言，MLE法雖在計算效率和計算準(zhǔn)確度方面均不如VCE法，但可以通過比較不同噪聲類型的似然值確定最佳噪聲類型和相應(yīng)的噪聲強(qiáng)度，而VCE法不能通過直接統(tǒng)計量確定噪聲類型。為了確定最佳的噪聲類型，本文采用w-test［3］進(jìn)行噪聲類型檢驗，w-test檢驗量近似服從正態(tài)分布，值越大閃爍噪聲的可接受度越高，檢驗結(jié)果如圖3所示。結(jié)果表明，除LUZH站U方向，XIAM站E方向（檢驗值分別為1.86和0.95）以外，其他站點(diǎn)坐標(biāo)時間序列的隨機(jī)性均符合白噪聲和閃爍噪聲的噪聲模型，故選擇基于白噪聲和閃爍噪聲的噪聲模型進(jìn)行噪聲方差的估計，估計結(jié)果如表3所示。

表3 中國區(qū)域GNSS站坐標(biāo)時間序列噪聲估計結(jié)果Tab.3 Estimation results of the noises in the GNSS coordinate time series in China

圖3 坐標(biāo)時間序列閃爍噪聲的w-test檢驗值Fig.3 Results of w-test for flicker noises in the coordinate time series

由表3可以看出，15個測站3個方向上的白噪聲強(qiáng)度均遠(yuǎn)遠(yuǎn)小于閃爍噪聲強(qiáng)度，表明中國區(qū)域GNSS站坐標(biāo)時間序列表現(xiàn)出明顯的有色噪聲特性。此外，豎直方向上的噪聲強(qiáng)度明顯高于水平方向的，造成這種現(xiàn)象的原因是：在坐標(biāo)解算過程中，水平方向分布的衛(wèi)星可以通過篩選，使其滿足圖形對稱性，而豎直方向分布的衛(wèi)星，受限于衛(wèi)星可見性，只能選擇地面以上的衛(wèi)星，無法滿足衛(wèi)星分布的對稱性，致使豎直方向的圖形強(qiáng)度較差。因此，豎直方向的噪聲水平高于水平方向的。

現(xiàn)有研究表明，水平方向和豎直方向的噪聲變化趨勢相同，為了更好地研究各個方向噪聲強(qiáng)度的相互關(guān)系，本文采用線性回歸分析了水平方向和豎直方向噪聲強(qiáng)度的線性關(guān)系，該結(jié)果可以為研究衛(wèi)星不對稱分布對定位結(jié)果的影響提供一定數(shù)據(jù)支持，同時為防止隨機(jī)抽樣的誤差影響回歸結(jié)果，對回歸模型進(jìn)行了顯著性水平為0.05的模型正確性檢驗，線性回歸結(jié)果見圖4，回歸結(jié)果見表4。

圖4為水平方向和豎直方向的白噪聲和閃爍噪聲的噪聲強(qiáng)度分布和擬合結(jié)果。由圖4可以看出，各噪聲間的擬合效果較好，除個別點(diǎn)外，大部分噪聲均位于擬合直線置信水平為95%的置信區(qū)間內(nèi)，結(jié)合表4中的擬合分位數(shù)，說明整體擬合效果較好。由表4可知，水平方向與豎直方向的噪聲水平存在一定相關(guān)性，相關(guān)水平為0.57～0.82，平均相關(guān)系數(shù)為0.68，其中E方向和U方向閃爍噪聲強(qiáng)度的相關(guān)性最大，相關(guān)系數(shù)為0.82，表明兩者間存在明顯相關(guān)性，而N方向和U方向的閃爍噪聲相關(guān)性最低，其相關(guān)系數(shù)僅為0.57。水平方向與豎直方向噪聲存在相關(guān)性的原因可能有兩點(diǎn)：（1）函數(shù)模型的不完善導(dǎo)致噪聲中含有未模型化的站點(diǎn)坐標(biāo)信息［23］（圖2）；（2）站點(diǎn)坐標(biāo)中含有大氣壓強(qiáng)、土壤濕度等地球物理因素引起的質(zhì)量荷載部分［24］，且由于中國區(qū)域GNSS站點(diǎn)跨度較大，常規(guī)方法無法完全剔除共模分量的影響，殘余的共模分量吸收到3個方向的噪聲中，進(jìn)一步提高了噪聲間的相關(guān)性。從擬合優(yōu)度看，水平方向上的噪聲可以引起豎直方向上的一部分噪聲，其中N方向的白噪聲和閃爍噪聲可以引起U方向53%和32%的相關(guān)噪聲，E方向可以引起U方向38%和67%的白噪聲和閃爍噪聲。

圖4 各方向噪聲分布和線性回歸結(jié)果Fig.4 Distribution and linear regression results of the noises in different directions

表4 各方向線性回歸結(jié)果Tab.4 Linear regression results of noises in different directions

4 結(jié) 論

（1）推導(dǎo)了LS-VCE法噪聲分量估計中的誤差，并通過仿真數(shù)據(jù)對比了極大似然估計法和兩種方差-協(xié)方差分量估計法對GNSS站坐標(biāo)時間序列噪聲模型估計的適用性。比較計算效率和計算準(zhǔn)確度后，確定了最佳噪聲模型估計法為LSVCE法，最佳坐標(biāo)時間序列長度為15 a。

（2）利用LS-VCE法對中國區(qū)域15個GNSS站坐標(biāo)時間序列進(jìn)行噪聲模型估計，分析結(jié)果表明，中國區(qū)域的閃爍噪聲強(qiáng)度遠(yuǎn)大于白噪聲強(qiáng)度，且豎直方向上的噪聲強(qiáng)度較高，水平方向上的噪聲強(qiáng)度較低?；谠肼暪烙嫿Y(jié)果，利用線性回歸法分析了水平方向噪聲與豎直方向噪聲的相關(guān)性，顯示水平方向噪聲能夠引起豎直方向32%以上的噪聲。