王樂洋,溫貴森
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é)果進行對比分析,驗證本文方法的可行性。
針對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ù)估值。
文獻[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)
隨機模型的驗后估計又稱方差分量估計,是針對平差前給定的隨機模型不準確而提出的方法。方差分量估計的基本思想是先對各類觀測值定初權(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
采用文獻[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ù)估值的標準差
根據(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)過偏差改正的相差較小,但都接近于驗前給定的方差分量。
在算例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ù)真值。
(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ù)展開的二階項時,二階近似要大于一階近似標準差,在線性化過程中忽略高次項時可能會舍去某些信息,而考慮了部分隨機量的隨機性的方法在理論上更加嚴密。
非線性模型在線性近似過程中會帶來模型誤差,導(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)生影響的程度進行分析,這也是今后將要開展的工作。