偏差改正的Partial EIV模型方差分量估計

王樂洋，溫貴森

1. 東華理工大學(xué)測繪工程學(xué)院，江西 南昌 330013； 2. 流域生態(tài)與地理環(huán)境監(jiān)測國家測繪地理信息局重點實驗室，江西 南昌 330013； 3. 江西省數(shù)字國土重點實驗室，江西 南昌 330013

在測量數(shù)據(jù)處理中，Gauss-Markov模型是常用的平差模型，模型解算常用的方法為最小二乘方法。然而在某些實際情況中，Gauss-Markov模型系數(shù)矩陣A的元素可能由某些含有觀測誤差的觀測值構(gòu)成，此時經(jīng)典的最小二乘平差方法不夠嚴密。總體最小二乘[1-4]是同時顧及了觀測向量與系數(shù)矩陣誤差的平差方法，是變量誤差(errors-in-variables，EIV)模型的嚴密估計方法。文獻[5—7]研究了總體最小二乘的迭代算法；文獻[8—14]研究了附有等式約束和不等式約束條件的總體最小二乘方法；文獻[15]研究了提高計算效率的總體最小二乘方法。部分變量誤差[16](partial error-in-variables，Partial EIV)模型是考慮了EIV模型系數(shù)矩陣中存在隨機元素和非隨機元素的情況，相比于EIV模型，Partial EIV模型更具有一般性和適用性，且在算法及應(yīng)用中都得到了廣泛的研究[17-25]。

由于EIV模型或Partial EIV模型受非線性的影響，總體最小二乘解算得到的參數(shù)估值及單位權(quán)方差估值是有偏的[26-27]。文獻[28]分析了系數(shù)矩陣誤差對EIV模型平差參數(shù)估值的影響。文獻[29]基于泰勒公式展開的二次項，得到了非線性函數(shù)的偏差公式。文獻[30]基于EIV模型，通過線性化處理得到了參數(shù)估值及參數(shù)估值一階近似協(xié)方差陣。文獻[27]對迭代的最后一步表達式進行線性化，得到了參數(shù)估值的一階近似協(xié)方差陣。文獻[16]根據(jù)非線性理論得到了參數(shù)估值偏差及一階近似協(xié)方差陣。文獻[31]通過迭代計算得到參數(shù)估值的協(xié)方差陣。文獻[32]推導(dǎo)了多元總體最小二乘算法并得到了參數(shù)估值的一階近似協(xié)方差陣。以上計算方法及通過協(xié)方差傳播得到的參數(shù)估值方差只計算到一階近似。文獻[33]使用SUT(scaled unscented transformation)方法計算得到總體最小二乘參數(shù)估值的二階近似協(xié)方差陣。文獻[34]推導(dǎo)了總體最小二乘精度評定的二階近似函數(shù)法，將參數(shù)估值精度提高至二階，進一步完善了總體最小二乘精度評定理論。

方差分量估計[35-39]是針對隨機模型不準確從而對權(quán)進行修正的驗后估計方法，方差分量估計方法通過預(yù)平差得到的信息根據(jù)一定的原則對觀測量的驗前方差或協(xié)方差進行估計、重新定權(quán)。文獻[40]分析了EIV模型的最小二乘方差分量估計。文獻[41]引入權(quán)修正因子，推導(dǎo)了Partial EIV模型的Helmert方差分量估計。文獻[42]得到了Partial EIV模型的非負最小二乘方差分量估計；然而以上總體最小二乘方差分量估計方法中并未考慮參數(shù)估值有偏性的影響。文獻[43]推導(dǎo)了EIV模型的最小范數(shù)二次無偏估計，分析了由參數(shù)估值偏差引起的方差分量估值的偏差，推導(dǎo)了方差分量估值偏差的表達式。

基于上述分析，本文考慮總體最小二乘估計參數(shù)估值的偏差，進行偏差改正得到偏差改正后的參數(shù)估值，更新由參數(shù)估值引起變化的中間變量，將偏差改正與方差分量歸為整體同時迭代計算，進行偏差改正后的方差分量估計及二階精度評定。通過算例試驗，將偏差改正后得到的結(jié)果與未進行偏差改正的結(jié)果進行對比分析，驗證本文方法的可行性。

1 Partial EIV模型參數(shù)偏差及精度評定

1.1 Partial EIV模型解算

針對EIV模型系數(shù)矩陣中存在固定元素及重復(fù)出現(xiàn)的隨機元素，文獻[16]將系數(shù)矩陣中的隨機元素分離，提出了Partial EIV模型，其數(shù)學(xué)模型為

函數(shù)模型

(1)

隨機模型

(2)

目前針對Partial EIV模型的解法較多，本文參考文獻[19]的求解思路，構(gòu)造拉格朗日條件極值函數(shù)

Φ=eTQ-1e+2λT(y-(βT?In)(h+Ba)+Ce)

(3)

對式(3)進行解算最終可得到參數(shù)估值的表達式為

(4)

(5)

由于總體最小二乘方法顧及了系數(shù)矩陣誤差，式(1)可以看作為非線性模型，Partial EIV模型解算得到的參數(shù)是有偏或近似無偏的[16]。將Partial EIV模型視為非線性模型得到

e=f(x)-l

(6)

(7)

(8)

(9)

通過迭代計算最終得到參數(shù)估值為

(10)

式(4)與式(10)是等價的，通過迭代計算都能得到相同的參數(shù)估值。

1.2 基于二階導(dǎo)數(shù)的偏差計算及精度評定

文獻[16]針對Partial EIV模型的非線性推導(dǎo)了參數(shù)估值的偏差表達式及參數(shù)估值的一階近似協(xié)方差陣。參數(shù)估值偏差表達式和協(xié)方差表達式分別為[16]

(11)

(12)

同理，式(4)通過協(xié)方差傳播律可以得到參數(shù)估值的一階近似協(xié)方差

(13)

(14)

(15)

(16)

(17)

式(15)可以得到參數(shù)估值的偏差

(18)

進而可以得到偏差改正后的參數(shù)估值

(19)

(20)

2 基于偏差改正的方差分量估計

隨機模型的驗后估計又稱方差分量估計，是針對平差前給定的隨機模型不準確而提出的方法。方差分量估計的基本思想是先對各類觀測值定初權(quán)，進行預(yù)平差，利用預(yù)平差后得到的信息，主要是各類觀測值的改正數(shù)，并根據(jù)得到的觀測值的改正數(shù)對各類觀測值的驗前方差和協(xié)方差進行估計，重新定權(quán)。

考慮式(2)隨機模型形式為

(21)

(22)

對總體最小二乘解算可以得到

(23)

(24)

式中，p在本文中取值為2，Qk的形式為

(25)

式中，k=1時，表示觀測量y對應(yīng)的初始協(xié)因數(shù)陣；k=2時，表示系數(shù)矩陣對應(yīng)的初始協(xié)因數(shù)陣。

本文以最小二乘方差分量估計方法為例，文獻[40]分析了總體最小二乘的最小二乘方差分量估計，將方差分量作為參數(shù)進行最小二乘平差解算得到

(26)

式中，vh表示取出對稱矩陣的上三角元素按照一定的順序排列成列向量。

文獻[36]通過公式轉(zhuǎn)換得到法矩陣N和列向量L的形式為

(27)

針對Partial EIV模型式(1)，模型的總觀測量個數(shù)為n+t，必要觀測量個數(shù)為t+m，在平差計算時是將系數(shù)矩陣含有觀測誤差的數(shù)據(jù)作為參數(shù)進行解算；根據(jù)非線性解算理論，參數(shù)β與系數(shù)矩陣觀測值a組成的向量x在進行式(10)的迭代計算時是存在偏差的，而參數(shù)估值影響觀測值改正數(shù)，從而影響方差分量估計。將偏差改正后的參數(shù)估值代入方差分量估計中，即存在

(28)

繼而有

(29)

式中，下標bc表示偏差改正后對應(yīng)的值。

(30a)

進而得到

(30b)

(31a)

(31b)

將參數(shù)估值偏差改正與方差分量估計作為一個整體進行迭代計算，其迭代流程見圖1。

圖1 偏差改正方差分量估計迭代流程Fig.1 Flow chart of bias correction variance components estimation iteration

3 算例與分析

3.1 算例1 直線擬合

采用文獻[18]的數(shù)據(jù)，直線擬合的模型為

(32)

已知坐標觀測值(xi,yi)和相應(yīng)的權(quán)值(pxi,pyi)，見表1。

考慮觀測數(shù)據(jù)的個數(shù)，構(gòu)造Partial EIV模型的向量h與矩陣B的形式為

(33)

現(xiàn)分別使用最小二乘(LS)、總體最小二乘(TLS)、總體最小二乘的最小二乘方差分量估計方法[40](TLS-VCE)及本文方法(偏差改正的方差分量估計)計算，得到的參數(shù)估值、參數(shù)估值偏差及方差分量估值結(jié)果見表2。

表1 坐標觀測值及相應(yīng)的權(quán)值[18]

表2 算例1不同方法得到的參數(shù)估值、參數(shù)估值偏差及方差分量估值

從表2可以看出，經(jīng)過偏差改正后的方差分量估計得到的方差分量估值及參數(shù)估值都有所改變，這主要是由參數(shù)估值的偏差引起的，參數(shù)估值偏差越大，對結(jié)果的影響也越大。計算得到參數(shù)估值的標準差結(jié)果見表3。

表3 算例1參數(shù)估值的標準差

3.2 算例2 二維仿射變換

根據(jù)文獻[38—39]算例思想，在長寬為100 m的正方形區(qū)域內(nèi)，沿縱橫坐標方向每隔10 m選取一個點，共121個點，在正方形區(qū)域的所有點中均勻選取31個點作為公共點進行二維仿射變換。二維仿射變換模型可以表示為

(34)

式中，Δx、Δy為平移參數(shù)；ωx、ωy為旋轉(zhuǎn)參數(shù)；kx、ky為尺度參數(shù)；(xi，yi)、(Xi，Xi)分別為第i個公共點的原始坐標系及目標坐標系的坐標。

將式(34)改寫成矩陣形式得到

(35)

式中，a1=kxcosωx；b1=kysinωy；a2=-kxsinωx；b2=kycosωy，構(gòu)造Partial EIV模型的向量h與矩陣B的形式為

(36)

給定一組轉(zhuǎn)換參數(shù)真值，Δx=Δy=0，ωx=100，ωy=110，kx=1.01，ky=1.02，假設(shè)目標坐標系坐標之間相互獨立且等精度，原始坐標同一點中x、y坐標獨立，不同點間坐標相關(guān)，對應(yīng)協(xié)因數(shù)陣中非對角線元素為

(37)

式中，qij為協(xié)因數(shù)陣中對應(yīng)的元素值；dij為i、j點間的距離；d0為常數(shù)且取值為10 m。

現(xiàn)分別給原始坐標與目標坐標模擬均值為0，服從正態(tài)分布的隨機誤差。其中，原始坐標單位權(quán)中誤差為1 cm，目標坐標單位權(quán)中誤差為3 cm。分別使用最小二乘(LS)、總體最小二乘(TLS)、總體最小二乘的最小二乘方差分量估計[40](TLS-VCE)及本文方法計算，得到的參數(shù)估值、參數(shù)估值偏差及參數(shù)估值與真值的差值范數(shù)見表4。

表4 算例2不同方法得到的參數(shù)估值、偏差及差值范數(shù)

從表4可以看出，由于該算例中參數(shù)估值偏差較小，不同方法得到的參數(shù)估值相差較小，參數(shù)估值與真值的差值范數(shù)也較小，但方差分量估計方法可以得到與平差前給定的方差分量相近的估值；計算得到原始坐標與目標坐標的方差分量估值見表5。

表5 算例2原始坐標與目標坐標方差分量估值

從表5可以看出，方差分量估計方法得到原始坐標的觀測精度為2.994 689 cm，目標坐標的觀測精度為1.142 152 cm，與平差前給定的3 cm和1 cm接近；而由于該算例中參數(shù)估值偏差較小，經(jīng)過偏差改正后的方差分量估值與未經(jīng)過偏差改正的相差較小，但都接近于驗前給定的方差分量。

3.3 算例3 四參數(shù)坐標轉(zhuǎn)換

在算例2中，由于觀測量誤差較小，在構(gòu)造總體最小二乘平差模型時，系數(shù)矩陣誤差很小，TLS與LS結(jié)果相差較??；文獻[28]指出系數(shù)矩陣的信噪比達到一定量級時，參數(shù)估值的相對偏差很小，LS與TLS結(jié)果必然沒有差別?？紤]一個四參數(shù)坐標轉(zhuǎn)換模型

(38)

根據(jù)所提供坐標個數(shù)，得到向量h、B形式如下

(39)

表6 算例3坐標觀測值及協(xié)因數(shù)

表7 算例3不同方法得到參數(shù)估值、偏差及方差分量估值

從表7與表8可以看出，方差分量估計方法可以得到與平差前相近的方差分量，對參數(shù)估值也有修正作用，而加入偏差改正后的方差分量估計方法得到的參數(shù)估值與真值的差值范數(shù)最小，數(shù)值為0.269 620，更接近于參數(shù)真值。

3.4 算例分析

(1) 算例1中，方差分量估計方法可以得到不同類數(shù)據(jù)的方差分量估值，而考慮總體最小二乘參數(shù)估值的有偏性時，對參數(shù)估值進行偏差改正，將改正后的參數(shù)估值與方差分量估計結(jié)合成整體進行計算，可以得到修正后的參數(shù)估值與方差分量估值。方差分量估值影響參數(shù)估值的精度評定，從表3的參數(shù)估值標準差可以看出，偏差改正后的參數(shù)估值標準差小于未偏差改正的參數(shù)估值標準差，而二階方差在數(shù)值上要大于一階，說明在非線性模型的精度評定時應(yīng)盡可能考慮線性化過程中舍去的高次項信息。

(2) 算例2中，由于觀測數(shù)據(jù)的量級與觀測誤差的量級差別較大，系數(shù)矩陣誤差對平差結(jié)果的影響較小，相應(yīng)地將平差模型作為非線性模型時，得到的參數(shù)估值偏差很小，量級為10-7，表4所提供的參數(shù)估值結(jié)果都與LS結(jié)果相差較?。槐?中方差分量估計方法都可以得到接近于平差前給定的方差分量，說明了方差分量估計方法的有效性；而由于參數(shù)估值偏差較小的原因，經(jīng)過偏差改正后得到的方差分量估值與未經(jīng)過偏差改正得到的方差分量估值接近一致，這也說明偏差大小與數(shù)據(jù)初始精度及模型的非線性強度有關(guān)[44]。

(3) 算例3中，考慮觀測數(shù)據(jù)量級較小且觀測量的方差較大的四參數(shù)坐標轉(zhuǎn)換模型，表7中給出計算得到參數(shù)估值的偏差最大為0.017 362，參數(shù)的偏差百分比為1.34%，模型具有非線性形態(tài)且有必要考慮泰勒展開式的二階項[44]，表7提供的參數(shù)估值與真值的差值范數(shù)中，經(jīng)過偏差改正后的方差分量估計得到的參數(shù)估值最接近于真值，相比于不考慮偏差改正的方差分量估計方法，更說明了偏差改正的必要性；從表8中的參數(shù)估值標準差可以看出，加入偏差改正后得到一階近似標準差要比不加偏差改正的小，而考慮泰勒級數(shù)展開的二階項時，二階近似要大于一階近似標準差，在線性化過程中忽略高次項時可能會舍去某些信息，而考慮了部分隨機量的隨機性的方法在理論上更加嚴密。

4 結(jié)束語

非線性模型在線性近似過程中會帶來模型誤差，導(dǎo)致得到的最小二乘估計的參數(shù)估值是有偏的。本文考慮了Partial EIV模型參數(shù)估值的偏差并得到偏差改正后的參數(shù)估值，在方差分量估計中，參數(shù)估值對觀測值改正數(shù)等有直接影響，進而影響方差分量估計的準確性。因此，本文將參數(shù)估計、偏差改正和方差分量估計作為整體進行迭代計算，將偏差改正后的參數(shù)估值代入方差分量估計中，從而獲得更加可靠的方差分量估值與參數(shù)估值。根據(jù)不同的平差模型及數(shù)據(jù)結(jié)構(gòu)，當系數(shù)矩陣誤差與系數(shù)矩陣量級很大時，模型的非線性強度較低，參數(shù)估值的偏差必然很小，此時進行偏差改正得到的結(jié)果變化不大；當偏差較大時，進行偏差改正可以得到更好的結(jié)果；得到更加合理的方差分量估值對隨機模型進行修正，從而獲得更加合理的參數(shù)估值及精度信息，是對總體最小二乘理論的進一步完善。本文將偏差改正結(jié)合到方差分量估計中，然而并未對偏差量級對結(jié)果產(chǎn)生影響的程度進行分析，這也是今后將要開展的工作。