等價條件平差模型的方差-協(xié)方差分量最小二乘估計方法

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

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

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

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

1 等價條件閉合差的方差分量最小二乘估計

1.1 等價條件平差模型

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

(1)

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

(2)

(3)

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

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

(4)

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

(5)

(6)

(7)

(8)

(9)

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

(10)

分析可知，式(8)—式(10)含有大量的高維矩陣迭代運算(r2階)，因而需進一步簡化矩陣運算?？紤]到矩陣拉直算子vec(·)、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)有文獻已證明殘差型VCE方法之間的等價性，并形成了一套較為完整的理論框架。通過概括平差因子矩陣R也可以證明本文LSV-ECM方法與殘差型VCE方法的等價性。對于概括平差模型，LSV-ECM法與Helmert型通用VCE方法[7](簡稱通用Helmert法)等價；對于間接平差模型，LSV-ECM法與LS-VCE法[9,11]、基于等效殘差的LS-VCE法[15]等價；對于具有參數(shù)的條件平差模型，LSV-ECM法與文獻[10]相應(yīng)算法等價。為節(jié)省篇幅，這里僅給出LSV-ECM法與通用Helmert法[7]等價關(guān)系的證明過程。

(14)

(15)

(16)

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

表2 不同VCE方法的計算復(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法進行方差分量估計和參數(shù)平差計算。其中測角和測邊單位權(quán)中誤差的迭代初值均取先驗中誤差，迭代終止條件為測角和測邊的單位權(quán)中誤差相等。

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

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

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

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

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

2.2 GNSS站坐標時序建模

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

GNSS站坐標時序的觀測方程和隨機模型分別為[11]

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

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

(17)

(18)

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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

表5 站點運動速度和不確定度