部分變量誤差模型的極大似然估計

于子龍, 胡友健, 汪奇生，肖本林

(1.中國地質(zhì)大學(xué)(武漢) 地理與信息工程學(xué)院,湖北 武漢 430078；2.湖北工業(yè)大學(xué) 土木建筑與環(huán)境學(xué)院，湖北 武漢 430068)

在測量數(shù)據(jù)處理中有些平差問題需要考慮系數(shù)矩陣的誤差，如線性回歸、坐標(biāo)轉(zhuǎn)換等，相應(yīng)的平差模型稱為變量誤差模型(errors in variables, EIV)，EIV模型的參數(shù)估計采用總體最小二乘法[1]，近年來受到測量領(lǐng)域的廣泛關(guān)注和研究，建立了相應(yīng)的平差理論并得到了較多的研究成果[1-3]，從EIV模型觀測值等精度情況下的總體最小二乘算法擴展到加權(quán)總體最小二乘算法[4]、非線性模型的總體最小二乘法[5,6]、觀測向量與系數(shù)矩陣元素相關(guān)的總體最小二乘算法[7]、觀測值含有粗差的穩(wěn)健總體最小二乘法[8]、平差模型附有不等式約束的總體最小二乘算法[9]，以及貝葉斯估計的總體最小二乘算法[10]等。由于總體最小二乘最初的解算采用的是奇異值分解法，不利于測量領(lǐng)域的應(yīng)用，特別是測量數(shù)據(jù)處理中通常要考慮觀測值不等精度的情況。因此，關(guān)于EIV模型的加權(quán)總體最小二乘算法得到較多研究，相應(yīng)提出一些迭代算法[11-22]。其中文獻(xiàn)[15]提出了一種算法迭代格式與最小二乘間接平差類似的算法，算法簡單，容易實現(xiàn)，而且迭代效率也較高。

由于EIV模型需要考慮系數(shù)矩陣誤差，但實際平差模型中系數(shù)矩陣構(gòu)造各不相同，可能同時存在隨機元素和常數(shù)項的情況。常規(guī)的EIV模型加權(quán)總體最小二乘算法對系數(shù)矩陣常數(shù)項構(gòu)造零權(quán)進(jìn)行處理增加了參數(shù)估計的個數(shù)。文獻(xiàn)[18]提出了部分變量誤差(partial EIV)模型，只提取系數(shù)矩陣的隨機元素項進(jìn)行處理，較大程度的減少了參數(shù)估計的個數(shù)，簡化了平差模型，但文獻(xiàn)[18]提出的迭代算法迭代格式復(fù)雜且計算效率較低[19-22]。文獻(xiàn)[19]提出了改進(jìn)的迭代算法，但算法仍然復(fù)雜需要進(jìn)一步簡化[20]，而且沒有驗證算法的效果，通過實驗發(fā)現(xiàn)其計算效率仍然較低。文獻(xiàn)[20]針對這一問題，通過對平差模型的部分元素進(jìn)行移項，通過兩次間接平差，提出一種新的算法(文中稱為算法1)，但新算法仍然存在收斂速度較慢的問題，通過進(jìn)一步分析與文獻(xiàn)[15]算法的關(guān)系，得到參數(shù)的另一種表達(dá)形式，提出一種迭代效率較高的算法(文中稱為算法2)。但文獻(xiàn)[20]提出的改進(jìn)算法并沒有從理論上進(jìn)行嚴(yán)密推導(dǎo)，而且沒有進(jìn)一步分析文獻(xiàn)[19]算法的關(guān)系。文獻(xiàn)[21]根據(jù)非線性平差提出一種新的算法，算法迭代格式也有待進(jìn)一步簡化。文獻(xiàn)[22]采用矩陣微分運算和矩陣反演推導(dǎo)一種新的算法，算法迭代格式其實與文獻(xiàn)[15]算法也類似，但公式推導(dǎo)過程較為復(fù)雜。王樂洋還對PEIV模型的其他擴展算法和應(yīng)用進(jìn)行大量研究[23-26]。

鑒于此，為了進(jìn)一步簡化PEIV模型加權(quán)總體最小二乘算法公式，提高PEIV模型的解算效率，本文參照文獻(xiàn)[10]的思路，采用極大似然估計方法進(jìn)一步推導(dǎo)PEIV模型的高效迭代算法，該算法的迭代格式與文獻(xiàn)[15]算法類似，易于理解和程序?qū)崿F(xiàn)。并將本文算法與文獻(xiàn)[18]算法、文獻(xiàn)[19]算法以及文獻(xiàn)[15]算法進(jìn)行對比和分析，驗證算法之間的等價性，最后采用兩個算例進(jìn)一步驗證本文算法的有效性和可行性。

1 PEIV平差模型

PEIV函數(shù)模型為：

(1)

(2)

則平差模型的隨機模型可以表示為：

(3)

不考慮系數(shù)矩陣元素與觀測值之間的相關(guān)性，隨機模型可以進(jìn)一步表示為：

(4)

(5)

文獻(xiàn)[18]提出的PEIV模型通過對系數(shù)矩陣的變換，只提取含有誤差的項來進(jìn)行處理，減少平差時待估參數(shù)的個數(shù)，但其提出的迭代算法效率并不高，收斂較慢且算法迭代格式復(fù)雜，F(xiàn)ang采用貝葉斯估計對EIV模型進(jìn)行參數(shù)求解[10]，提出了迭代算法，但經(jīng)過實驗驗證發(fā)現(xiàn)其迭代效率依然不高。因此，本文根據(jù)其思路，采用極大似然估計對PEIV模型參數(shù)求解進(jìn)行公式推導(dǎo)，并通過進(jìn)一步推導(dǎo)提出迭代效率較高的算法。

2 PEIV模型的極大似然估計

極大似然估計的原理，假設(shè)觀測值服從正態(tài)分布，構(gòu)建觀測值的似然函數(shù)，使得似然函數(shù)取最大值，通過對似然函數(shù)求極值的方法來估計待求參數(shù)，根據(jù)式(2)、式(3)可知觀測值的期望和方差可以表示為：

E(y)=(βT?Im)(h+BUa)=vec-1(h+BUa)β,

(6)

式中:Ua為系數(shù)矩陣中隨機元素a的期望，vec-1表示的是矩陣?yán)边\算的逆運算即重構(gòu)矩陣，根據(jù)式(6)，可構(gòu)造觀測值的似然函數(shù)：

(7)

式中:含有兩個待求參數(shù)β、Ua，極大似然估計的準(zhǔn)則為對似然函數(shù)取最大值，即將分別對兩個參數(shù)進(jìn)行求導(dǎo)，并令其等于零。為了推算方便，對式(7)求對數(shù)仍然是等價的，進(jìn)一步可以表示為：

lnf(l|β,Ua)=

(8)

將式(8)分別對兩個參數(shù)求導(dǎo)并進(jìn)行相應(yīng)簡化，可以表示為：

(9)

(10)

另式(9)等于0，可得到參數(shù)β的表達(dá)式：

β=((vec-1(h+BUa))TPyvec-1(h+BUa))-1×

(11)

參照文獻(xiàn)[10]Fang的推導(dǎo)，本文不再贅述，令式(10)等于0，可導(dǎo)出觀測值誤差和系數(shù)矩陣元素誤差之間的關(guān)系：

(12)

相應(yīng)的顧及到ey+(βT?Im)Bea=y-(βT?Im)(h+Ba)=y-Aβ，則可進(jìn)一步表示為：

ey=Qy(Qy+(βT?Im)BQaBT(βT?Im))-1×(y-Aβ),ea=-QaBT(β?Im)(Qy+(βT?Im)BQaBT×(βT?Im))-1(y-Aβ).

(13)

式(13)的表達(dá)式與文獻(xiàn)[10]Fang的推導(dǎo)是等價的，根據(jù)文獻(xiàn)[10]，則可以將迭代算法(本文算法1)表示為:

由于算法的推導(dǎo)和表達(dá)式與Fang基本相同，差別僅在于本文所述的是PEIV模型而Fang所述的是EIV模型，兩者可以進(jìn)行轉(zhuǎn)換，因此將上述方法仍稱為Fang法，其迭代格式與最小二乘相似，易于編程實現(xiàn)，但通過實驗分析，F(xiàn)ang法的迭代效率較低，收斂較慢。通過上述推導(dǎo)可以發(fā)現(xiàn)，式(13)中僅僅用到系數(shù)矩陣元素的誤差向量，而觀測值的誤差向量表達(dá)式并沒有用來進(jìn)一步推導(dǎo)，鑒于此，本文對式(13)進(jìn)行進(jìn)一步推導(dǎo)。結(jié)合上述推導(dǎo)，式(1)的第一式可以表示為：

y=(βT?Im)(h+BUa)+ey.

(14)

將式(13)中的ey帶入到式(14)中：

y=vec-1(h+BUa)β+Qy(Qy+(βT?Im)×

(βT?Im)BQaBT(β?Im))-1(y-Aβ).

(15)

再將式(15)代入到式(11)中：

BQaBT(β?Im))-1(y-Aβ).

(16)

(17)

(18)

式(18)的參數(shù)求解公式，相對于式(11)利用了更多的信息，經(jīng)過實驗分析迭代效率有較大提高，另外，式(8)的迭代格式與文獻(xiàn)[15] Jazaeri的EIV模型迭代算法是等價的，而本文基于PEIV模型更具一般性。

PEIV模型極大似然估計算法(本文算法2，為了方便在后文敘述中無特殊區(qū)分，所指的本文算法皆為本文算法2)的迭代步驟為：

PEIV模型的單位權(quán)方差的評定公式為[18-22]：

(19)

(20)

3 本文算法與文獻(xiàn)中算法的比較分析

3.1 與Xu算法的比較分析

由文獻(xiàn)[18]Xu提出的PEIV模型，相應(yīng)的構(gòu)建了PEIV模型的目標(biāo)函數(shù)：

(21)

y}=0.

(22)

進(jìn)一步可得到Xu算法的迭代公式：

β=(Nh+NB+NhB+NBh)-1(uh+uB).

(23)

Xu在對式(22)第二式的β求解時，沒有采用矩陣的拉直運算和重構(gòu)矩陣運算，實際上式(22)第二式目標(biāo)函數(shù)對β求導(dǎo)的表達(dá)式與本文式(9)是等價的，僅為表達(dá)形式不一樣。其中式(23)第二式中(Nh+NB+NhB+NBh)=(vec-1(h+BUa))TPyvec-1(h+BUa)，(uh+uB)=(vec-1(h+BUa))TPyy，因此，Xu算法與本文所述的Fang算法是完全等價的。

3.2 與Shi算法的比較分析

在文獻(xiàn)[19]中Shi提出了PEIV模型的改進(jìn)算法，根據(jù)式(21)的目標(biāo)函數(shù)對β求導(dǎo)并令其等于0，可得到參數(shù)β的表達(dá)式：

(24)

(25)

對比式(24)、式(25)與式(23)，Shi提出的表達(dá)式要比Xu的簡化，但進(jìn)一步比較可以發(fā)現(xiàn)，Shi提出的算法與本文所述的Fang算法是等價的，在迭代效率上并沒有提高。對比式(24)與本文的式(11)可以發(fā)現(xiàn)，兩個表達(dá)式是一致的。對比式(25)與本文的式(13)第二式，很明顯也是完全一致的。因此，雖然Shi提出的迭代算法公式要比Xu的算法簡單，但本質(zhì)上并無差別，只是公式的表達(dá)式不一樣而已。Shi提出的算法與Xu的算法都與本文所述的Fang算法等價，經(jīng)過實驗分析發(fā)現(xiàn)其迭代效率也基本一致。

3.3 與Jazaeri算法的比較分析

文獻(xiàn)[15]中，Jazaeri提出的迭代加權(quán)總體最小二乘算法是基于EIV模型，其構(gòu)造的目標(biāo)函數(shù)為：

(26)

(27)

vec(EA)=-QA(β?Im)Q-1(y-Aβ).

(28)

3.4 與Wang算法的比較分析

文獻(xiàn)[20]中，Wang通過采用兩次間接平差原理分別求解平差參數(shù)和系數(shù)矩陣元素，提出了一種算法(文中稱為算法1)，平差參數(shù)和系數(shù)矩陣元素的求解公式為：

(29)

Pa)-1(BT(β?Im)Py(y-(βT?Im)h)+aPa)=

a+QaBT(β?Im)(Qy+(βT?Im)×

BQaBT(β?Im))-1(y-Aβ)=a-ea.

(30)

4 算例分析

4.1 算例1：二維相似坐標(biāo)轉(zhuǎn)換

為了說明算法的適應(yīng)性，采用一個二維相似坐標(biāo)轉(zhuǎn)換的算例來進(jìn)行驗證，選用文獻(xiàn)[15]中的坐標(biāo)轉(zhuǎn)換數(shù)據(jù)如表1所示，表中有5個具有原始坐標(biāo)和目標(biāo)坐標(biāo)兩套坐標(biāo)系統(tǒng)的二維點，其相應(yīng)的權(quán)值為：

Pxy=diag([18.781 7,6.377 4,12.648 9,17.476 9,22.272 6,23.982 3,13.680 4,3.465 6,3.732 4,6.437 7]);PXY=diag([9.831 6,5.535 7,12.736 9,12.009 9,10.181,11.366 1,11.147,5.883 4,9.832 2,7.567 8]).

用表1中的坐標(biāo)數(shù)據(jù)來計算二維相似坐標(biāo)轉(zhuǎn)換參數(shù)，根據(jù)平差原理可知其也是一個典型的求解部分誤差變量誤差模型參數(shù)的問題。同樣采用本文所述的Fang算法、本文算法、Jazaeri算法、Xu算法和Shi算法5種方法進(jìn)行坐標(biāo)參數(shù)求解，二維相似坐標(biāo)轉(zhuǎn)換的參數(shù)有4個，分別設(shè)為ξ1、ξ2、ξ3、ξ4，算法的收斂條件為δ=10-12，各算法估計的坐標(biāo)轉(zhuǎn)換參數(shù)結(jié)果列于表2，為了進(jìn)一步分析各算法的計算效率，統(tǒng)計各個算法計算一次的迭代次數(shù)以及重復(fù)1 000次計算的總耗時和平均耗時列于表3，繪制幾種算法在每次迭代過程中前后兩次參數(shù)的差值范數(shù)于圖1。為了比較各算法解算結(jié)果的差異，統(tǒng)計本文算法與其他算法參數(shù)估計結(jié)果的差值列于表4。

圖1 各算法的迭代過程

表1 原始坐標(biāo)系和目標(biāo)坐標(biāo)系數(shù)據(jù) m

表2 坐標(biāo)轉(zhuǎn)換各算法估計結(jié)果

表3 坐標(biāo)轉(zhuǎn)換各算法計算效率

表4 本文算法與其他算法估計結(jié)果的差異

從表2中可以發(fā)現(xiàn)，本文算法與其他幾種算法得到的結(jié)果是一致的，這說明本文算法也能解算PEIV模型參數(shù)估計問題，從表3的計算一次的迭代數(shù)中可以發(fā)現(xiàn)，本文算法與Jazaeri算法的迭代次數(shù)一致且同樣最少，其他3種算法迭代次數(shù)相同且較多，從1000次重復(fù)計算的總耗時和平均耗時可以看出，本文算法與Jazaeri算法接近且較其他3種算法要少，圖1的參數(shù)迭代過程也清晰的表明本文算法與Jazaeri算法的收斂較其他3種算法要快，從表4的參數(shù)估計結(jié)果差異中可以看出，本文算法與其他算法解算結(jié)果差異最大都在10-12，少于本文迭代的收斂閾值，可認(rèn)為各算法的解算結(jié)果是一致的，對比本文算法與其他算法的差異可以發(fā)現(xiàn)，與Jazaeri算法差異為0，這也進(jìn)一步驗證了本文算法與Jazaeri算法完全等價的結(jié)論；與Fang算法差異和Shi算法差異完全相同，這也說明了前文推導(dǎo)中兩種算法是等價的結(jié)論；與Xu算法差異雖然與其他兩種算法差異的結(jié)果不同，但差異值很小，可能是Xu算法采用的大量矩陣離散計算沒有采用其他算法的矩陣?yán)焙椭貥?gòu)而造成的差異，但都在本文的迭代閾值之內(nèi)。綜上說明本文算法在計算結(jié)果上能得到正確的值且計算效率具有一定的優(yōu)勢，進(jìn)一步驗證了本文算法的有效性和可行性。

4.2 算例2：二維仿射坐標(biāo)轉(zhuǎn)換

為了進(jìn)一步驗證本文算法的有效性，采用一個模擬的二維仿射坐標(biāo)轉(zhuǎn)換算例，算例引用文獻(xiàn)[12]，13個具有原始坐標(biāo)和目標(biāo)坐標(biāo)真值的二維點列于表5，二維仿射坐標(biāo)轉(zhuǎn)換的6個轉(zhuǎn)換參數(shù)的真值分別為a1=0.9，b1=-0.8，c1=1，a2=0.6，b2=0.7，c2=5，根據(jù)表5的坐標(biāo)真值模擬添加誤差進(jìn)行實驗，本文在原始坐標(biāo)和目標(biāo)坐標(biāo)添加的誤差協(xié)因素陣分別設(shè)為：

表5 二維仿射坐標(biāo)轉(zhuǎn)換數(shù)據(jù)

Qxy=I2?(0.005×diag([1,3,6,1,1,8,4,3,6,5,4,5,2]))，

QXY=I2?(0.005×diag([1,2,3,1,5,4,2,7,2,1,8,3,6])).

根據(jù)誤差協(xié)因素陣共進(jìn)行5 000次模擬，分別采用本文所述的Fang算法、本文算法、Jazaeri算法、Xu算法和Shi算法5種方法進(jìn)行坐標(biāo)參數(shù)求解，算法的收斂條件為δ=10-12，統(tǒng)計每一次模擬各種算法解算的結(jié)果，將各種方法獲得的參數(shù)估計均值列于表6，統(tǒng)計各算法的平均迭代次數(shù)和迭代計算耗時列于表7，繪制5 000次實驗中各算法的迭代時間于圖2，統(tǒng)計本文算法與其他算法參數(shù)估計結(jié)果的平均差值列于表8。

表6 二維仿射坐標(biāo)轉(zhuǎn)換各算法估計結(jié)果的均值

表7 二維仿射坐標(biāo)轉(zhuǎn)換各算法計算效率

表8 本文算法與其他算法估計結(jié)果的平均差異

圖2 5 000次試驗各算法的迭代時間

從表6的各算法計算結(jié)果中可以看出，本文算法和其他幾種算法得到的結(jié)果一致，而且二維仿射轉(zhuǎn)換參數(shù)的均值與真值比較接近，這說明本文推導(dǎo)的算法能得到正確的結(jié)果。從表7的算法計算效率中可以發(fā)現(xiàn)，本文算法與Jazaeri算法所需的平均迭代次數(shù)為9.077，相對于其他3種算法迭代收斂速度要快，這說明本文算法要優(yōu)于其他3種算法，從迭代耗時上可以看出，本文算法與Jazaeri算法相對于其他3種算法耗時要少，而且本文算法比Jazaeri算法迭代收斂要快，這可能是因為本文算法采用的PEIV模型較Jazaeri算法采用的EIV模型在處理系數(shù)矩陣常數(shù)項時更方便，從圖2的5000次試驗各個算法的迭代時間中也可以說明本文算法相對于其他幾種算法迭代耗時較少，迭代效率較高。從表8的本文算法與其他算法解算結(jié)果的平均差異中可以進(jìn)一步說明，本文算法與Jazaeri算法解算的結(jié)果完全一致，在迭代公式上僅針對系數(shù)矩陣的處理上有差異，本文采用的是PEIV模型，而Jazaeri采用的是EIV模型，對于系數(shù)矩陣含有大量常數(shù)項的情況，本文算法處理起來更方便。對于其他3種算法，本文算法與它們的差異都在迭代收斂閾值之內(nèi)，另外還可以發(fā)現(xiàn)與Fang算法差異和與Shi算法差異是完全相同的，這也驗證了前文推導(dǎo)證明兩種算法是完全等價的，僅在算法迭代公式的表述上存在一些差異。

5 結(jié)束語

PEIV模型只提取系數(shù)矩陣含誤差的項進(jìn)行處理，減少了常規(guī)EIV模型系數(shù)矩陣的估計參數(shù)個數(shù)，一定程度上簡化了平差模型，但Xu提出的算法和Shi改進(jìn)的算法計算效率較低且算法較復(fù)雜，本文基于極大似然估計原理推導(dǎo)了一種新的PEIV模型算法，該算法的迭代格式與最小二乘類似，易于理解和程序?qū)崿F(xiàn)，通過與Xu算法、Shi算法以及Jazaeri算法的對比和分析，驗證了算法之間的等價性，最后通過兩個算例進(jìn)行了分析驗證，結(jié)果表明本文算法能得到與其他算法一致的結(jié)果，且算法迭代格式簡單，計算效率較高。