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

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

(1.中國地質大學(武漢) 地理與信息工程學院,湖北 武漢 430078；2.湖北工業大學 土木建筑與環境學院，湖北 武漢 430068)

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

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

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

1 PEIV平差模型

PEIV函數模型為：

(1)

(2)

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

(3)

不考慮系數矩陣元素與觀測值之間的相關性，隨機模型可以進一步表示為：

(4)

(5)

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

2 PEIV模型的極大似然估計

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

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

(6)

式中:Ua為系數矩陣中隨機元素a的期望，vec-1表示的是矩陣拉直運算的逆運算即重構矩陣，根據式(6)，可構造觀測值的似然函數：

(7)

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

lnf(l|β,Ua)=

(8)

將式(8)分別對兩個參數求導并進行相應簡化，可以表示為：

(9)

(10)

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

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

(11)

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

(12)

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

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)的表達式與文獻[10]Fang的推導是等價的，根據文獻[10]，則可以將迭代算法(本文算法1)表示為:

由于算法的推導和表達式與Fang基本相同，差別僅在于本文所述的是PEIV模型而Fang所述的是EIV模型，兩者可以進行轉換，因此將上述方法仍稱為Fang法，其迭代格式與最小二乘相似，易于編程實現，但通過實驗分析，Fang法的迭代效率較低，收斂較慢。通過上述推導可以發現，式(13)中僅僅用到系數矩陣元素的誤差向量，而觀測值的誤差向量表達式并沒有用來進一步推導，鑒于此，本文對式(13)進行進一步推導。結合上述推導，式(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)的參數求解公式，相對于式(11)利用了更多的信息，經過實驗分析迭代效率有較大提高，另外，式(8)的迭代格式與文獻[15] Jazaeri的EIV模型迭代算法是等價的，而本文基于PEIV模型更具一般性。

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

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

(19)

(20)

3 本文算法與文獻中算法的比較分析

3.1 與Xu算法的比較分析

由文獻[18]Xu提出的PEIV模型，相應的構建了PEIV模型的目標函數：

(21)

y}=0.

(22)

進一步可得到Xu算法的迭代公式：

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

(23)

Xu在對式(22)第二式的β求解時，沒有采用矩陣的拉直運算和重構矩陣運算，實際上式(22)第二式目標函數對β求導的表達式與本文式(9)是等價的，僅為表達形式不一樣。其中式(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算法的比較分析

在文獻[19]中Shi提出了PEIV模型的改進算法，根據式(21)的目標函數對β求導并令其等于0，可得到參數β的表達式：

(24)

(25)

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

3.3 與Jazaeri算法的比較分析

文獻[15]中，Jazaeri提出的迭代加權總體最小二乘算法是基于EIV模型，其構造的目標函數為：

(26)

(27)

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

(28)

3.4 與Wang算法的比較分析

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

(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：二維相似坐標轉換

為了說明算法的適應性，采用一個二維相似坐標轉換的算例來進行驗證，選用文獻[15]中的坐標轉換數據如表1所示，表中有5個具有原始坐標和目標坐標兩套坐標系統的二維點，其相應的權值為：

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

圖1 各算法的迭代過程

表1 原始坐標系和目標坐標系數據 m

表2 坐標轉換各算法估計結果

表3 坐標轉換各算法計算效率

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

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

4.2 算例2：二維仿射坐標轉換

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

表5 二維仿射坐標轉換數據

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])).

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

表6 二維仿射坐標轉換各算法估計結果的均值

表7 二維仿射坐標轉換各算法計算效率

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

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

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

5 結束語

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