偏差改正的Partial EIV模型方差分量估計(jì)

王樂(lè)洋，溫貴森

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

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

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

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

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

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

1.1 Partial EIV模型解算

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

函數(shù)模型

(1)

隨機(jī)模型

(2)

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

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

(3)

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

(4)

(5)

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

e=f(x)-l

(6)

(7)

(8)

(9)

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

(10)

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

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

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

(11)

(12)

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

(13)

(14)

(15)

(16)

(17)

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

(18)

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

(19)

(20)

2 基于偏差改正的方差分量估計(jì)

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

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

(21)

(22)

對(duì)總體最小二乘解算可以得到

(23)

(24)

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

(25)

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

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

(26)

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

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

(27)

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

(28)

繼而有

(29)

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

(30a)

進(jìn)而得到

(30b)

(31a)

(31b)

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

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

3 算例與分析

3.1 算例1 直線擬合

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

(32)

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

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

(33)

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

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

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

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

表3 算例1參數(shù)估值的標(biāo)準(zhǔn)差

3.2 算例2 二維仿射變換

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

(34)

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

將式(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è)目標(biāo)坐標(biāo)系坐標(biāo)之間相互獨(dú)立且等精度，原始坐標(biāo)同一點(diǎn)中x、y坐標(biāo)獨(dú)立，不同點(diǎn)間坐標(biāo)相關(guān)，對(duì)應(yīng)協(xié)因數(shù)陣中非對(duì)角線元素為

(37)

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

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

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

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

表5 算例2原始坐標(biāo)與目標(biāo)坐標(biāo)方差分量估值

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

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

在算例2中，由于觀測(cè)量誤差較小，在構(gòu)造總體最小二乘平差模型時(shí)，系數(shù)矩陣誤差很小，TLS與LS結(jié)果相差較小；文獻(xiàn)[28]指出系數(shù)矩陣的信噪比達(dá)到一定量級(jí)時(shí)，參數(shù)估值的相對(duì)偏差很小，LS與TLS結(jié)果必然沒(méi)有差別。考慮一個(gè)四參數(shù)坐標(biāo)轉(zhuǎn)換模型

(38)

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

(39)

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

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

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

3.4 算例分析

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

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

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

4 結(jié)束語(yǔ)

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