劉志平，朱丹彤，余 航，張克非,2

1. 中國礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院, 江蘇 徐州 221116； 2. 皇家墨爾本理工大學(xué)空間科學(xué)研究中心，澳大利亞 維多利亞州 墨爾本 3001

隨著現(xiàn)代測(cè)量技術(shù)的發(fā)展，測(cè)量數(shù)據(jù)處理的對(duì)象已經(jīng)由單一同類觀測(cè)數(shù)據(jù)轉(zhuǎn)變?yōu)槎嘣础⒍囝惢蛲惗嘁蛩氐漠惙讲罱Y(jié)構(gòu)觀測(cè)數(shù)據(jù)，在這種情況下隨機(jī)模型的準(zhǔn)確性對(duì)參數(shù)解估計(jì)的統(tǒng)計(jì)特性具有重要影響[1-2]。方差-協(xié)方差分量估計(jì)能夠極大地改善隨機(jī)模型不精確的問題，從而成為國內(nèi)外學(xué)者關(guān)注的焦點(diǎn)，并形成了一整套的理論框架[3-8]。近年來主要理論研究包括：①以平差結(jié)果之一的殘差向量為基本輸入量的殘差型VCE方法，如Helmert法[1-2,7]、最小二乘方差分量估計(jì)法(簡稱LS-VCE)[9-11]、MINQUE法[12]、BIQUE法[13]、MLE-VCE法[14]、等效殘差法[15-16]及EIV模型的VCE方法[4,8]等；②以等價(jià)條件閉合差為基本輸入量的解析型VCE方法，如VCE-ECM法[17-18]。近年來的主要應(yīng)用研究有：文獻(xiàn)[19—20]研究了常規(guī)導(dǎo)線控制網(wǎng)中角度觀測(cè)值與距離觀測(cè)值的定權(quán)問題；文獻(xiàn)[21—22]探討了GNSS觀測(cè)值隨機(jī)模型的精化問題；文獻(xiàn)[23—24]研究了坐標(biāo)轉(zhuǎn)換中不同公共點(diǎn)坐標(biāo)的精度問題；文獻(xiàn)[11,25—26]開展了GNSS站坐標(biāo)時(shí)間序列的噪聲估計(jì)及評(píng)價(jià)研究。

大量研究表明，LS-VCE方法[9-10,15]可導(dǎo)出Helmert、MINQUE、MLE-VCE等方法的方差-協(xié)方差分量估計(jì)公式，具有最優(yōu)無偏特性。與此同時(shí)，基于等價(jià)條件平差模型的VCE-ECM法[17-18]以等價(jià)條件閉合差作為基本輸入量，摒棄了平差值求解或殘差估計(jì)過程。鑒于此，本文從方差-協(xié)方差分量最小二乘估計(jì)準(zhǔn)則和等價(jià)條件平差模型出發(fā)，利用等價(jià)條件閉合差的二次型構(gòu)造方差分量估計(jì)公式，進(jìn)而將方差分量估計(jì)方程變換為線性Gauss-Markov形式并建立了基本估計(jì)方程；然后，顧及矩陣?yán)?、半拉直算子和Kronecker積的性質(zhì)進(jìn)行化簡，導(dǎo)出了基于等價(jià)條件閉合差的方差-協(xié)方差分量最小二乘估計(jì)公式(least-square variance-covariance estimation based on the equivalent condition misclosure)，簡稱為LSV-ECM法。在此基礎(chǔ)上，證明了本文LSV-ECM法與殘差型VCE方法的等價(jià)性，定量分析了LSV-ECM法、Helmert法和LS-VCE法的計(jì)算復(fù)雜度，以邊角網(wǎng)平差和中國區(qū)域GNSS站坐標(biāo)時(shí)序建模實(shí)例驗(yàn)證本文方法的正確性和高效性。

1 等價(jià)條件閉合差的方差分量最小二乘估計(jì)

1.1 等價(jià)條件平差模型

設(shè)概括平差模型為

(1)

式中，A、[BTCT]均為行滿秩系數(shù)矩陣；W為具有參數(shù)的條件方程閉合差；Z為限制條件方程閉合差；DL為觀測(cè)值L的方差陣；V、x分別為待求的殘差與參數(shù)向量；下標(biāo)c、s分別為條件方程個(gè)數(shù)和限制條件方程個(gè)數(shù)；n、u分別為觀測(cè)數(shù)和參數(shù)個(gè)數(shù)；平差模型自由度r=(c+s)-u。

(2)

(3)

1.2 方差-協(xié)方差分量最小二乘估計(jì)方法

假設(shè)概括平差模型式(1)中隨機(jī)模型DL可表示為

(4)

式中，∑(·)表示累加運(yùn)算；k為方差-協(xié)方差分量的個(gè)數(shù)；Qi為第i個(gè)協(xié)因數(shù)分量矩陣。

(5)

(6)

(7)

(8)

(9)

式中，D為復(fù)制矩陣(duplication matrix)[9]，也即半拉直算子與拉直算子的映射矩陣

(10)

分析可知，式(8)—式(10)含有大量的高維矩陣迭代運(yùn)算(r2階)，因而需進(jìn)一步簡化矩陣運(yùn)算。考慮到矩陣?yán)彼阕觱ec(·)、Kronecker積?和矩陣跡算子tr(·)具有如下性質(zhì)[9]

vec(T)T·(U?S)·vec(VT)=tr(UVST)

(11)

(12)

(13)

參量模型條件平差具有參數(shù)的條件平差間接平差附有限制的間接平差B=0,C=0C=0A=-I,C=0A=-IQiAQiATHAQiATHTHQiHTHcQiHTcWWHWHWHcW+HsZDWADLATHADLATHTHDLHTHcDLHTc

現(xiàn)有文獻(xiàn)已證明殘差型VCE方法之間的等價(jià)性，并形成了一套較為完整的理論框架。通過概括平差因子矩陣R也可以證明本文LSV-ECM方法與殘差型VCE方法的等價(jià)性。對(duì)于概括平差模型，LSV-ECM法與Helmert型通用VCE方法[7](簡稱通用Helmert法)等價(jià)；對(duì)于間接平差模型，LSV-ECM法與LS-VCE法[9,11]、基于等效殘差的LS-VCE法[15]等價(jià)；對(duì)于具有參數(shù)的條件平差模型，LSV-ECM法與文獻(xiàn)[10]相應(yīng)算法等價(jià)。為節(jié)省篇幅，這里僅給出LSV-ECM法與通用Helmert法[7]等價(jià)關(guān)系的證明過程。

(14)

(15)

(16)

分析對(duì)比式(15)—式(16)和式(12)—式(13)可知，LSV-ECM方法可以導(dǎo)出通用Helmert法[7]，等價(jià)性得證。其次，殘差型VCE方法(如Helmert法、LS-VCE法等)均以殘差向量為輸入量，而本文LSV-ECM法以等價(jià)條件閉合差為輸入量，后者實(shí)現(xiàn)了平差值求解(殘差為平差結(jié)果之一)與隨機(jī)模型估計(jì)的分離。再次，等價(jià)條件閉合差維數(shù)低于殘差向量維數(shù)，且前者為不需求解的已知量，因此計(jì)算效率方面明顯優(yōu)于前者。

表2 不同VCE方法的計(jì)算復(fù)雜度

Tab.2 A comparison of computational complexities of different VCE methods

2 應(yīng)用結(jié)果及分析

2.1 邊角網(wǎng)平差

方案1：分別采用通用Helmert法[7]、LS-VCE法[9,11]、本文LSV-ECM法進(jìn)行方差分量估計(jì)和參數(shù)平差計(jì)算。其中測(cè)角和測(cè)邊單位權(quán)中誤差的迭代初值均取先驗(yàn)中誤差，迭代終止條件為測(cè)角和測(cè)邊的單位權(quán)中誤差相等。

方案2：在方案1的基礎(chǔ)上，剔除含有粗差的12號(hào)角度觀測(cè)值和16號(hào)邊長觀測(cè)值[18]，基于余下的16個(gè)觀測(cè)值進(jìn)行方差分量估計(jì)和參數(shù)平差計(jì)算。

按上述方案分別進(jìn)行計(jì)算，統(tǒng)計(jì)不同方法迭代收斂時(shí)的方差分量估值、參數(shù)估值及其中誤差。另外，為避免單次計(jì)算的隨機(jī)誤差，將上述方案重復(fù)計(jì)算100次，并統(tǒng)計(jì)LS-VCE法和本文LSV-ECM法相對(duì)于通用Helmert法的耗時(shí)比值Tratio，所得結(jié)果如表3所示。

表3 不同VCE方法的估計(jì)結(jié)果

從方案1方差分量和待估參數(shù)的估計(jì)結(jié)果看，3種結(jié)果完全相同，表明3種方法均能獲得一致的估計(jì)結(jié)果。同時(shí)比較3種方法的計(jì)算效率，通用Helmert法和LS-VCE的計(jì)算效率處于同一水平，本文所提出的LSV-ECM法的計(jì)算效率最高，其計(jì)算時(shí)間約為通用Helmert法的65%，較通用Helmert法和LS-VCE法均有明顯提升，從而驗(yàn)證本文方法的有效性。

方案2中同樣可以驗(yàn)證本文所提方法的正確性和有效性，不再贅述。同時(shí)，方案2在剔除兩個(gè)粗差觀測(cè)值以后，兩類觀測(cè)值的單位權(quán)方差、參數(shù)的估計(jì)精度均有較大提升，表明3種方差分量估計(jì)方法都無法抑制粗差影響，必須引入質(zhì)量控制方法進(jìn)行粗差識(shí)別與剔除。

2.2 GNSS站坐標(biāo)時(shí)序建模

為進(jìn)一步檢驗(yàn)本文LSV-ECM法的正確性和高效性，選用中國大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)(crustal movement observation network of China，CMONOC)中16個(gè)長期線性趨勢(shì)穩(wěn)定的GNSS基準(zhǔn)站2005.0014(DOY)～2015.0014(DOY)共10年的原始坐標(biāo)時(shí)間序列進(jìn)行噪聲特性分析(方差-協(xié)方差分量估計(jì))和三維運(yùn)動(dòng)速度估計(jì)(平差參數(shù)計(jì)算)，測(cè)站位置分布如圖1所示。

GNSS站坐標(biāo)時(shí)序的觀測(cè)方程和隨機(jī)模型分別為[11]

y(ti)=a+bti+ccos(2πti)+dsin(2πti)+

ecos(4πti)+fsin(4πti)+vi

(17)

(18)

式中，ti是以年為單位的時(shí)間序列歷元點(diǎn)；a為常數(shù)項(xiàng)；b為線性速度項(xiàng)；c、d組合表示全年性周期運(yùn)動(dòng)；e、f組合表示半年性周期運(yùn)動(dòng)；vi為殘差；σWN、σFN為所求的噪聲分量大?。籕WN、QFN分別為白噪聲和閃爍噪聲的協(xié)因數(shù)陣，具體形式見文獻(xiàn)[11]。

圖1 CMONOC所選測(cè)站分布Fig.1 Geographical distribution of selected stations in the CMONOC network

表4 噪聲分量(WN,FN)估計(jì)結(jié)果

注：高斯白噪聲(WN)的單位為mm；閃爍噪聲(FN)的單位為mm/a0.25

對(duì)比表4中通用Helmert法、LS-VCE法和LSV-ECM法的計(jì)算結(jié)果可知，3種方法所得的噪聲分量結(jié)果完全相同，驗(yàn)證了本文方法的正確性。同時(shí)，對(duì)比3種方法的計(jì)算效率可知，通用Helmert法與LS-VCE法基本一致，而本文LSV-ECM法有較大提高，三維方向計(jì)算時(shí)間比分別為通用Helmert法的74.0%、73.8%、72.1%，表明本文方法的高效性。

由表4噪聲分量的估計(jì)結(jié)果可知，白噪聲分量均遠(yuǎn)小于閃爍噪聲分量，表明有色噪聲為中國區(qū)域GNSS站坐標(biāo)時(shí)間序列的主要噪聲，在參數(shù)估計(jì)時(shí)若直接采用白噪聲模型會(huì)導(dǎo)致估計(jì)結(jié)果有偏，并會(huì)產(chǎn)生較高的虛假估計(jì)精度。另外，對(duì)比不同方向的噪聲分量估計(jì)結(jié)果，北方向和東方向的噪聲分量相差不大，且量級(jí)較小，94%和88%的測(cè)站白噪聲分量在1 mm以內(nèi)，50%的測(cè)站閃爍噪聲分量在3 mm/a0.25以內(nèi)，而豎直方向的噪聲分量遠(yuǎn)高于水平方向，這與現(xiàn)有結(jié)論相一致，約有56%的測(cè)站白噪聲分量在3 mm以內(nèi)，31%的測(cè)站閃爍噪聲分量在9 mm/a0.25以內(nèi)。

為進(jìn)一步分析噪聲大小與經(jīng)緯度之間的關(guān)系，繪制了白噪聲分量和閃爍噪聲分量隨經(jīng)緯度的變化情況。由于噪聲大小與緯度相關(guān)性較強(qiáng)、與經(jīng)度相關(guān)性較弱，為節(jié)省篇幅，此處僅討論噪聲大小與緯度變化的關(guān)系，見圖2。從圖2可以看出，水平方向的白噪聲和閃爍噪聲分量大小較為平穩(wěn)，豎直方向的噪聲大小波動(dòng)較大。從整體趨勢(shì)上看，白噪聲和閃爍噪聲大小隨緯度變化的趨勢(shì)較為明顯。具體表現(xiàn)為：噪聲大小隨緯度增加而逐漸減小，減少過程中有輕微波動(dòng)，且在中緯度地區(qū)有翹尾現(xiàn)象。該現(xiàn)象原因分析可能是GNSS的GDOP值隨緯度增大而降低[28]，從而使噪聲分量表現(xiàn)出隨緯度增加而減小的現(xiàn)象，而且噪聲分量在中緯度地區(qū)呈現(xiàn)一定的翹尾現(xiàn)象。

此外，為檢驗(yàn)方差-協(xié)方差分量估計(jì)結(jié)果的正確性和有效性，以CMONOC的公布數(shù)據(jù)(http:∥www.cgps.ac.cn/cgs/index.action)作為對(duì)比參考值，表5統(tǒng)計(jì)了中國區(qū)域GNSS站的三維運(yùn)動(dòng)速度和不確定度。由該表可知，中國區(qū)域不同站點(diǎn)水平方向均呈現(xiàn)東南方向運(yùn)動(dòng)趨勢(shì)，這與亞歐板塊的整體運(yùn)動(dòng)趨勢(shì)相一致，但豎直方向運(yùn)動(dòng)變化差異性較大，其中約40%的測(cè)站呈現(xiàn)下沉趨勢(shì)。對(duì)比本文結(jié)果與CMONOC的公布結(jié)果，除JIXN、HLAR站N方向，KMIN、GUAN站E方向，CHUN、BJFS站U方向以外，其余計(jì)算結(jié)果均在2倍中誤差范圍內(nèi)。此外，參照文獻(xiàn)[29]的計(jì)算結(jié)果，對(duì)比發(fā)現(xiàn)除CHUN站U方向外，其余測(cè)站的計(jì)算結(jié)果與本文結(jié)果具有較好的相符性，進(jìn)一步驗(yàn)證了本文方法的正確性。需要指出的是，CHUN站U方向的差異性可能是數(shù)據(jù)預(yù)處理和共模誤差提取策略不同所致。

圖2 噪聲分量大小與緯度的關(guān)系Fig.2 The relationship between noise components and latitude

3 結(jié) 論

本文基于等價(jià)條件平差模型和最小二乘準(zhǔn)則，利用等價(jià)條件閉合差的二次型構(gòu)建方差-協(xié)方差分量估計(jì)方程，并通過矩陣半拉直算子將其變換為線性Gauss-Markov形式，進(jìn)而顧及矩陣?yán)彼阕?、半拉直算子和Kronecker積運(yùn)算性質(zhì)，導(dǎo)出了基于等價(jià)條件平差模型的方差-協(xié)方差分量最小二乘估計(jì)公式，簡稱LSV-ECM法。該方法實(shí)現(xiàn)了平差值求解(殘差為平差結(jié)果之一)與隨機(jī)模型估計(jì)的分離，有效兼顧了等價(jià)條件閉合差(已知)和最小二乘的特性。邊角網(wǎng)平差實(shí)例的結(jié)果表明，本文的LSV-ECM法與通用Helmert法、LS-VCE法的結(jié)果完全相同，但計(jì)算效率更高，驗(yàn)證了該方法與殘差型VCE方法的等價(jià)性和計(jì)算高效性。

分析指出了方差-協(xié)方差分量估計(jì)方法相較于常規(guī)GNSS-MLE法進(jìn)行GNSS站坐標(biāo)時(shí)間序列噪聲分析的優(yōu)勢(shì)，并利用LSV-ECM法、通用Helmert法、LS-VCE法計(jì)算分析了中國區(qū)域16個(gè)GNSS站坐標(biāo)時(shí)間序列在WN+FN模型下的噪聲估計(jì)和站點(diǎn)速度。估計(jì)結(jié)果表明，有色噪聲是中國GNSS站坐標(biāo)時(shí)序的主要噪聲，且白噪聲和有色噪聲分量隨緯度增大而減小，速度估計(jì)結(jié)果與陸態(tài)網(wǎng)絡(luò)的公布結(jié)果基本一致，水平方向的整體運(yùn)動(dòng)趨勢(shì)呈東南方向，而豎直方向的運(yùn)動(dòng)趨勢(shì)差異性較大。因此，計(jì)算結(jié)果也進(jìn)一步驗(yàn)證了本文LSV-ECM法的正確性和高效性。

表5 站點(diǎn)運(yùn)動(dòng)速度和不確定度