王樂洋，陳 濤

1. 東華理工大學(xué)測(cè)繪工程學(xué)院，江西 南昌 330013; 2. 自然資源部環(huán)鄱陽湖區(qū)域礦山環(huán)境監(jiān)測(cè)與治理重點(diǎn)實(shí)驗(yàn)室，江西 南昌 330013; 3. 武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢 430079

在大地測(cè)量數(shù)據(jù)處理領(lǐng)域中，測(cè)量平差理論是在基于加性誤差的線性或非線性高斯-馬爾可夫(Gauss-Markov,GM)[1-2]模型基礎(chǔ)上發(fā)展起來的。隨著現(xiàn)代測(cè)繪技術(shù)的高速發(fā)展，以電磁波為代表的觀測(cè)值的隨機(jī)誤差表現(xiàn)為乘性誤差或加乘性混合誤差[3-4]，例如LiDAR觀測(cè)值的隨機(jī)誤差表現(xiàn)為加乘性混合誤差，電磁波測(cè)距(electronic distance measurement,EDM)、全球?qū)Ш叫l(wèi)星系統(tǒng)(GNSS)和甚長(zhǎng)基線干涉測(cè)量(very long baseline interferometry,VLBI)基線精度計(jì)算中，乘常數(shù)和加常數(shù)本質(zhì)上也是加乘性混合誤差[3-7]。傳統(tǒng)的加性誤差模型平差理論不再適用于新型加乘性混合誤差模型，若直接采用基于加性誤差的平差理論進(jìn)行解算，模型不能正確反映觀測(cè)值隨機(jī)誤差的特征，此時(shí)求得的結(jié)果將不可靠，因此需要針對(duì)加乘性混合誤差模型進(jìn)行更深入的研究來豐富現(xiàn)代大地測(cè)量數(shù)據(jù)處理理論。

目前，加乘性混合誤差模型的解法包括擬極大似然法和最小二乘法[5]。雖然擬極大似然法是廣義線性模型或乘性誤差模型參數(shù)估計(jì)的主要方法，但它存在以下問題[5]：①?zèng)]有明確定義的目標(biāo)函數(shù)；②擬似然解是由一個(gè)非線性方程組決定的，可能存在多組解，無法確定最似然解。乘性誤差模型的最小二乘法最早是由文獻(xiàn)[8]提出，文獻(xiàn)[5]在乘性誤差模型的基礎(chǔ)上最早給出了加乘性混合誤差模型的函數(shù)模型和隨機(jī)模型，并基于最小二乘原理推導(dǎo)了最小二乘解法和加權(quán)最小二乘迭代解法，此外還通過公式推導(dǎo)將文獻(xiàn)[8]中乘性誤差模型的偏差改正加權(quán)最小二乘迭代解法推廣到加乘性混合誤差模型，但該方法是通過參數(shù)估值的迭代公式根據(jù)誤差傳播定律求得參數(shù)估值的精度信息，由于在進(jìn)行泰勒函數(shù)展開時(shí)忽略了求偏導(dǎo)數(shù)的二階及其高階項(xiàng)，故只能達(dá)到一階精度。文獻(xiàn)[9]給出了乘性誤差和加性誤差彼此獨(dú)立且各自具有相同方差情況下的偏差改正加權(quán)最小二乘迭代解法，但該方法也只能求得參數(shù)估值的一階精度。因此需要針對(duì)加乘性混合誤差模型參數(shù)估值的精度評(píng)定問題作進(jìn)一步的研究。

無跡變換(unscented transformation,UT)方法是一種確定性采樣策略的近似非線性函數(shù)概率密度分布法，最早由文獻(xiàn)[10—11]在非線性Kalman濾波領(lǐng)域中首次提出，其通過選擇特定的采樣策略來近似非線性函數(shù)的概率分布，無須求導(dǎo)、無須了解非線性函數(shù)的構(gòu)成，通過加權(quán)的方法便可計(jì)算非線性函數(shù)的均值和協(xié)方差陣。經(jīng)過20多年的發(fā)展，不同的學(xué)者提出了多種采樣策略，擴(kuò)展了其應(yīng)用范圍。文獻(xiàn)[12—15]提出了一種適應(yīng)性強(qiáng)的比例對(duì)稱采樣策略，相應(yīng)的UT即為SUT。文獻(xiàn)[13]證明了通過SUT法計(jì)算的非線性函數(shù)的均值和協(xié)方差陣均可以達(dá)到二階精度。文獻(xiàn)[16]證明了將SUT法應(yīng)用在測(cè)量非線性函數(shù)誤差傳播中是可行和有效的，并通過與二階近似函數(shù)法進(jìn)行對(duì)比，證明了SUT法求得的參數(shù)估值及其精度信息均可以達(dá)到二階精度。文獻(xiàn)[17]使用SUT法研究了變量誤差(errors-in-variables,EIV)模型的參數(shù)估計(jì)和精度評(píng)定問題。文獻(xiàn)[18]研究了部分變量誤差(partial errors-in-variables,Partial EIV)模型方差分量估計(jì)精度評(píng)定問題。文獻(xiàn)[19]研究了Okada模型中震源參數(shù)估值對(duì)滑動(dòng)分布反演中格林函數(shù)矩陣的影響分析并進(jìn)行精度評(píng)定。文獻(xiàn)[20]研究了Okada模型的震源參數(shù)的精度評(píng)定問題。EIV模型、Partial EIV模型及Okada模型中的參數(shù)估值和觀測(cè)值之間為一個(gè)復(fù)雜的非線性迭代過程，加乘性混合誤差模型中參數(shù)估值和觀測(cè)值之間也是一個(gè)非常復(fù)雜的非線性迭代過程，因此，SUT法可以較好地用于加乘性混合誤差模型的精度評(píng)定。

綜上所述，針對(duì)已有加乘性混合誤差模型精度評(píng)定方法只能達(dá)到一階精度，且加乘性混合誤差模型中參數(shù)估值與觀測(cè)值為一個(gè)復(fù)雜的非線性關(guān)系，無法通過求導(dǎo)運(yùn)算和公式推導(dǎo)進(jìn)行精度評(píng)定的問題，本文使用一種無須求導(dǎo)、無須了解非線性函數(shù)構(gòu)造的SUT法對(duì)其進(jìn)行精度評(píng)定，旨在求得參數(shù)估值的二階精度信息；最后，通過算例驗(yàn)證本文將SUT法應(yīng)用于加乘性混合誤差模型精度評(píng)定問題中的可行性及優(yōu)勢(shì)。

1 加乘性混合誤差模型及其解法

1.1 加乘性混合誤差模型

加乘性混合誤差模型可以表示為[5]

y=f(β)⊙(1+εm)+εa

(1)

式中，y∈Rn×1表示觀測(cè)值向量；f(·)表示未知參數(shù)的函數(shù)；β∈Rt×1表示未知參數(shù)向量；⊙表示Hadamard積符號(hào)；1∈Rn×1表示元素全為1的列向量；εm∈Rn×1表示服從正態(tài)分布的乘性誤差列向量；εa∈Rn×1表示服從正態(tài)分布的加性誤差列向量。

本文研究的函數(shù)f(·)是β的線性函數(shù)形式[4-7,9]，則式(1)可表示為

y=(Aβ)⊙(1+εm)+εa

(2)

式中，A∈Rn×t表示系數(shù)矩陣，且A=[a1,a2,…,at](ai∈Rn×1)。

E(y)=Aβ

(3)

D(y)=E[(y-E(y))(y-E(y))T]=

(4)

由式(4)可知，在加乘性混合誤差模型中，觀測(cè)值的協(xié)方差陣D(y)是參數(shù)估值β的非線性函數(shù)，這明顯區(qū)別于傳統(tǒng)的加性誤差模型，即觀測(cè)值的協(xié)方差陣D(y)與參數(shù)估值β無關(guān)，因此傳統(tǒng)的加性誤差模型平差理論不再適用于新型的加乘性混合誤差模型。

1.2 加乘性混合誤差模型的最小二乘解

對(duì)式(2)應(yīng)用最小二乘準(zhǔn)則，可得目標(biāo)函數(shù)為

min:F1(β)=(y-Aβ)T(y-Aβ)

(5)

根據(jù)求解目標(biāo)函數(shù)自由極值的原理，將函數(shù)F1(β)對(duì)β求偏導(dǎo)，并令其為0，可得加乘性混合誤差模型的最小二乘解為

(6)

(7)

1.3 加乘性混合誤差模型的加權(quán)最小二乘迭代解

對(duì)式(2)應(yīng)用加權(quán)最小二乘準(zhǔn)則，可得目標(biāo)函數(shù)為

(8)

根據(jù)求解目標(biāo)函數(shù)自由極值的原理，將函數(shù)F2(β)對(duì)β求偏導(dǎo)，并令其為0，可得

(9)

由上述分析可知，由于加乘性混合誤差模型中觀測(cè)值的權(quán)是參數(shù)估值的非線性函數(shù)，因此通過式(9)無法得到參數(shù)估值的解析解，只能使用數(shù)值逼近的方法進(jìn)行迭代求得數(shù)值解，進(jìn)而可以得到加乘性混合誤差模型的加權(quán)最小二乘迭代解為

ATP(y)iy+

(10)

由文獻(xiàn)[5]可得加權(quán)最小二乘迭代解的單位權(quán)中誤差和一階協(xié)方差陣分別為

(11)

(12)

1.4 加乘性混合誤差模型的偏差改正加權(quán)最小二乘迭代解

文獻(xiàn)[5]通過公式推導(dǎo)證明了加權(quán)最小二乘解的偏差等于觀測(cè)值的權(quán)對(duì)參數(shù)估值求偏導(dǎo)數(shù)項(xiàng)，因此，在加權(quán)最小二乘解的基礎(chǔ)上刪去觀測(cè)值的權(quán)對(duì)參數(shù)估值求偏導(dǎo)數(shù)項(xiàng)(即式(10)右邊第二項(xiàng))得到了參數(shù)估值的二階無偏估計(jì)，對(duì)應(yīng)的偏差改正加權(quán)最小二乘解為

(13)

式(13)中觀測(cè)值的權(quán)是參數(shù)估值的非線性函數(shù)，而參數(shù)估值又是觀測(cè)值的權(quán)的非線性函數(shù)，因此式(13)也為一個(gè)迭代解，其迭代格式為

(14)

由文獻(xiàn)[5]可得偏差改正加權(quán)最小二乘迭代解的單位權(quán)中誤差和一階協(xié)方差陣分別為

(15)

(16)

由式(16)可以看出，文獻(xiàn)[5]是通過參數(shù)估值的迭代公式和協(xié)方差傳播率來獲取參數(shù)估值的精度信息，由于忽略了求偏導(dǎo)數(shù)的二階及其高階項(xiàng)，故只能達(dá)到一階精度。

綜上所述，本文將文獻(xiàn)[5]中加乘性混合誤差模型的偏差改正加權(quán)最小二乘迭代解法總結(jié)為算法1，具體步驟如下。

(1) 輸入觀測(cè)值y，系數(shù)矩陣A，計(jì)算最小二乘估值作為迭代初值

(17)

(2) 計(jì)算第i次迭代觀測(cè)值的權(quán)P(y)i

(18)

(19)

(20)

(3) 由式(14)計(jì)算加乘性混合誤差模型的偏差改正加權(quán)最小二乘迭代解。

(5) 由式(15)和式(16)分別計(jì)算偏差改正加權(quán)最小二乘迭代解的單位權(quán)中誤差和一階協(xié)方差陣。

2 加乘性混合誤差模型精度評(píng)定的SUT法

由算法1可以看出，偏差改正加權(quán)最小二乘迭代解法中每次迭代的參數(shù)估值與觀測(cè)值的函數(shù)關(guān)系相同，因此在經(jīng)過第i次迭代后，最終得到的參數(shù)估值的迭代過程可以表示為

(21)

式(21)表明最終得到的參數(shù)估值與觀測(cè)值之間為一個(gè)復(fù)雜的非線性函數(shù)關(guān)系，此時(shí)若通過近似函數(shù)法來進(jìn)行精度評(píng)定，必然需要復(fù)雜的求導(dǎo)運(yùn)算。為了避開復(fù)雜的求導(dǎo)運(yùn)算和公式推導(dǎo)，同時(shí)得到參數(shù)估值的二階精度信息，本文在偏差改正加權(quán)最小二乘迭代解法中引入無須求導(dǎo)、無須了解非線性函數(shù)構(gòu)成的SUT采樣法進(jìn)行精度評(píng)定，旨在獲得參數(shù)估值的二階精度信息。由于式(21)是一個(gè)復(fù)雜的非線性函數(shù)，無法將其具體表達(dá)式表示出來，但其本質(zhì)上也是一個(gè)非線性函數(shù)，因此，不失一般性，現(xiàn)從理論上證明通過SUT法求得一般的函數(shù)Y=F(y)的因變量的均值和協(xié)方差陣均可以達(dá)到二階精度，具體如下。

在式(21)中，因變量Y對(duì)自變量y的Jacobian矩陣和Hessian矩陣可表示為[23]

(22)

(23)

(24)

進(jìn)一步對(duì)函數(shù)Y=F(y)進(jìn)行二階泰勒展開可得

(25)

式中，Δy=y-E(y)；G=[ΔyTb1ΔyΔyTb2Δy…ΔyTbtΔy]T。

根據(jù)期望和協(xié)方差陣的定義[21]，對(duì)式(25)進(jìn)行求期望和協(xié)方差陣運(yùn)算可得

(26)

(27)

式中，?表示Kronecker積[24]。

由式(26)和式(27)可以看出，若通過傳統(tǒng)泰勒級(jí)數(shù)展開的近似函數(shù)方法進(jìn)行運(yùn)算，需要求得Jacobi矩陣a和Hessian矩陣c，但當(dāng)非線性函數(shù)復(fù)雜且非線性強(qiáng)時(shí)，很難進(jìn)行求取。因此，為了確保因變量的均值和協(xié)方差陣在誤差傳播中能達(dá)到二階精度，并且避免求導(dǎo)，利用SUT法求取因變量的均值和協(xié)方差陣，對(duì)于函數(shù)Y=F(y)，SUT法精度評(píng)定的具體步驟為[12-13]。

(1) 基于自變量y的先驗(yàn)均值E(y)和協(xié)方差陣D(y)依據(jù)式(28)生成采樣點(diǎn)集{χi}(i=0,1,…,2n)

(28)

(2) 將采樣點(diǎn)集{χi}(i=0,1,…,2n)代入函數(shù)Y=F(y)中，求得對(duì)應(yīng)的解集

(29)

(3) 加權(quán)計(jì)算因變量Y的均值E(Y)SUT和協(xié)方差陣D(Y)SUT

(30)

(31)

將通過SUT法求得的因變量Y的均值E(Y)SUT和協(xié)方差陣D(Y)SUT，即式(30)和式(31)具體展開可得[16]

(32)

(33)

(34)

(35)

由于在SUT法中，α=10-3，γ=2，κ=0，因此式(35)進(jìn)一步可表示為

(36)

對(duì)比式(26)與式(34)可知，兩者表達(dá)式的形式完全相同，表明利用SUT法求得函數(shù)Y=F(y)的均值可以精確地達(dá)到二階精度；同時(shí)，對(duì)比式(27)與式(36)可知，兩者表達(dá)式的形式也完全相同，表明利用SUT法求得函數(shù)Y=F(y)的協(xié)方差陣也可以精確地達(dá)到二階精度。因此，表明本文利用SUT法求解加乘性混合誤差模型，所求得的參數(shù)估值及其協(xié)方差陣均能達(dá)到二階精度。

綜上所述，本文將加乘性混合誤差模型精度評(píng)定的SUT法總結(jié)為算法2，具體步驟如下。

p=0,1,…,2n

(37)

(4) 由式(38)和式(39)計(jì)算參數(shù)估值的單位權(quán)中誤差和標(biāo)準(zhǔn)差

(38)

(39)

由上述步驟可以看出，SUT法能夠解決復(fù)雜或難以實(shí)現(xiàn)的求導(dǎo)計(jì)算，具有操作簡(jiǎn)單、易于實(shí)現(xiàn)等特點(diǎn)，且利用SUT法求解加乘性混合誤差模型能夠有效避免復(fù)雜的求導(dǎo)運(yùn)算，所求得的參數(shù)估值及其協(xié)方差陣均能達(dá)到二階精度。

3 算例與分析

為驗(yàn)證本文提出方法在加乘性混合誤差模型中精度評(píng)定方面的性能及探討在大地測(cè)量領(lǐng)域中的應(yīng)用，本文共設(shè)計(jì)了3個(gè)算例：通過算例1來展示本文方法在GPS高程擬合中的應(yīng)用價(jià)值，算例2和算例3來展示本文方法在數(shù)字地面模型(DTM)中的應(yīng)用價(jià)值，并驗(yàn)證本文方法在精度評(píng)定方面的優(yōu)勢(shì)。分別使用最小二乘解法、加權(quán)最小二乘迭代解法、文獻(xiàn)[5]中的偏差改正加權(quán)最小二乘迭代解法(即算法1)及加乘性混合誤差模型精度評(píng)定的SUT法(即算法2)進(jìn)行對(duì)比解算以驗(yàn)證本文方法的可行性及優(yōu)勢(shì)，4種方案對(duì)應(yīng)的方法分別列于表1中。

表1 4種方案及其對(duì)應(yīng)的方法

3.1 算例1

為初步驗(yàn)證本文方法在加乘性混合誤差模型中的可行性和正確性，算例1采用文獻(xiàn)[30]中的GPS高程擬合算例，該算例中高程真值與離某一坐標(biāo)原點(diǎn)的距離除以100滿足以下函數(shù)模型

y=10+x+2x2+2x3

(40)

假設(shè)高程真值受加乘性混合誤差的影響，具體可表示為[7]

Y=y⊙(1+εm)+εa

(41)

式中，Y表示31維高程觀測(cè)值向量；y表示31維高程真值向量；1表示元素全為1的31維列向量；εm表示獨(dú)立且服從正態(tài)分布的31維乘性誤差列向量；εa表示獨(dú)立且服從正態(tài)分布的31維加性誤差列向量。

根據(jù)文獻(xiàn)[5—7]，本算例將乘性誤差的標(biāo)準(zhǔn)差設(shè)置為0.05，加性誤差的標(biāo)準(zhǔn)差設(shè)置為0.15 m，單位權(quán)中誤差設(shè)置為σ0=0.3 m，隨機(jī)誤差由Matlab中的mvnrnd函數(shù)生成。表2列出了地面點(diǎn)的高程模擬值，為了說明高程值受加乘性混合誤差的影響程度，將不含隨機(jī)誤差時(shí)的高程真值及受加乘性混合誤差干擾的高程值繪于圖1中。

表2 地面點(diǎn)的高程值

圖1 不含隨機(jī)誤差時(shí)的高程真值及受加乘性混合誤差干擾的高程值Fig.1 The true elevations without random error and the elevations disturbed by mixed additive and multiplicative random error

由圖1可知，當(dāng)高程受到加乘性混合誤差的干擾時(shí)，觀測(cè)值會(huì)受到一定的影響，進(jìn)而與原始的高程真值產(chǎn)生了偏離。

圖2 WLS法和bcWLS法的流程Fig.2 Flowchart of WLS method and bcWLS method

圖3 SUT法的流程Fig.3 Flowchart of SUT method

表3 算例1中LS法、WLS法、bcWLS法及SUT法求得的參數(shù)估值、估值與真值之間的二范數(shù)和單位權(quán)中誤差

由表3的結(jié)果可知，LS法由于在平差時(shí)未考慮觀測(cè)值的權(quán)，因此求得的二范數(shù)結(jié)果最大，且單位權(quán)中誤差也嚴(yán)重偏離真值；WLS法由于在LS法的基礎(chǔ)上考慮了觀測(cè)值的權(quán)，因此求得的二范數(shù)及單位權(quán)中誤差結(jié)果均優(yōu)于LS法；bcWLS法由于在WLS法的基礎(chǔ)上減去了其偏差，因此求得的二范數(shù)的結(jié)果優(yōu)于WLS法，表明文獻(xiàn)[5]中方法的有效性。對(duì)比bcWLS法和SUT法可以看出，兩者求得的二范數(shù)結(jié)果一樣，而文獻(xiàn)[5]已通過嚴(yán)格的數(shù)學(xué)公式推導(dǎo)證明了bcWLS法求得的參數(shù)估值能達(dá)到二階精度，說明本文通過SUT法求得的參數(shù)估值也能達(dá)到二階精度，表明本文方法的可行性。

表4 算例1中LS法、WLS法、bcWLS法及SUT法求得的參數(shù)估值的標(biāo)準(zhǔn)差

由表4求得的標(biāo)準(zhǔn)差結(jié)果可知，LS法求得的結(jié)果最大，WLS法的結(jié)果小于LS法，而bcWLS法的結(jié)果略小于WLS法，精度最高，這與文獻(xiàn)[5]得到的結(jié)論一致，進(jìn)一步表明在平差時(shí)考慮觀測(cè)值權(quán)的必要性。然而，由于bcWLS法是基于參數(shù)估值的最后一步迭代表達(dá)式，通過協(xié)方差傳播率求得參數(shù)估值的精度信息，因此只能達(dá)到一階精度。對(duì)比bcWLS法和SUT法求得的結(jié)果可以看出，SUT法求得的標(biāo)準(zhǔn)差最小，精度最高，而上文已證明利用SUT法求得加乘性混合誤差模型參數(shù)估值的協(xié)方差陣可以達(dá)到二階精度，說明通過SUT法求得的參數(shù)估值的二階精度信息的情形下中誤差是減小的，這與文獻(xiàn)[28—29]得到的結(jié)論一致。此外，若參數(shù)估值的二階精度信息情形下中誤差是增大的，此時(shí)參數(shù)估值的二階標(biāo)準(zhǔn)差將大于bcWLS法求得的標(biāo)準(zhǔn)差，而WLS法求得的標(biāo)準(zhǔn)差也大于bcWLS法求得的標(biāo)準(zhǔn)差，說明通過WLS法求得的標(biāo)準(zhǔn)差為二階精度，但這顯然是不合理的，因此進(jìn)一步說明通過SUT法求得參數(shù)估值的二階精度信息的情形下中誤差是減小的，從而驗(yàn)證了本文方法的可行性和優(yōu)勢(shì)。

3.2 算例2

為進(jìn)一步驗(yàn)證本文方法在加乘性混合誤差模型中的可行性和正確性，算例2是一個(gè)DTM模型算例，該模型通過有限的地形高程數(shù)據(jù)實(shí)現(xiàn)對(duì)地面地形的數(shù)字化模擬，被廣泛用于工程、軍事、遙感與制圖、地理分析等方面，例如在工程項(xiàng)目中的挖填方計(jì)算、線路勘測(cè)設(shè)計(jì)、航天遙感數(shù)字圖像定量解釋中的應(yīng)用等[31-32]。目前，機(jī)載LiDAR被廣泛用于DTM數(shù)據(jù)的獲取[33-34]，但是以LiDAR技術(shù)為數(shù)據(jù)源構(gòu)建DTM模型都是基于加性誤差來進(jìn)行處理的[35-37]，然而LiDAR數(shù)據(jù)的噪聲已被證明具有乘性噪聲性質(zhì)[38-40]，且文獻(xiàn)[5—7,9]已將LiDAR數(shù)據(jù)的誤差當(dāng)作加乘性混合誤差進(jìn)行處理，因此本算例將仿真設(shè)計(jì)受加乘性混合誤差干擾的LiDAR DTM模型試驗(yàn)，以更好地驗(yàn)證本文方法的優(yōu)勢(shì)。本算例模擬的含4個(gè)未知參數(shù)的加乘性混合誤差的DTM模型觀測(cè)方程如下

H(x,y)=h(x,y)⊙(1+εm)+εa

(42)

(43)

式中，(x,y)表示地面點(diǎn)坐標(biāo)，橫坐標(biāo)x和縱坐標(biāo)y都在0～100 m之間，以步長(zhǎng)為2 m取值；H(x,y)表示DTM模型高程觀測(cè)值；h(x,y)表示DTM模型高程真值；1表示元素全為1的列向量；εm表示獨(dú)立且服從正態(tài)分布的乘性誤差列向量；εa表示獨(dú)立且服從正態(tài)分布的加性誤差列向量；βi表示未知待求參數(shù)；函數(shù)fi(·)(i=1,2,3,4)的具體形式為

(44)

根據(jù)文獻(xiàn)[38—40]中對(duì)LiDAR數(shù)據(jù)的理論和實(shí)際測(cè)試及已有關(guān)于加乘性混合誤差模型的文獻(xiàn)[5—7]，本算例將加性誤差的標(biāo)準(zhǔn)差固定為σa=0.15 m，將乘性誤差的標(biāo)準(zhǔn)差分別設(shè)置為σm=0.005、σm=0.05和σm=0.1這3組值進(jìn)行對(duì)比解算。為了說明DTM模型受加乘性混合誤差的影響程度，將不含隨機(jī)誤差時(shí)的DTM模型繪于圖4中，將受到加乘性混合誤差干擾的DTM模型分別繪于圖5、圖6和圖7中。令單位權(quán)中誤差σ0=0.3 m[5-7]，分別使用表1中的4種算法計(jì)算，將參數(shù)真值、參數(shù)估值及其與真值的二范數(shù)和單位權(quán)中誤差列于表5，參數(shù)估值的標(biāo)準(zhǔn)差列于表6。

通過對(duì)比圖4、圖5、圖6和圖7可以發(fā)現(xiàn)，隨機(jī)誤差對(duì)高程產(chǎn)生了較大的影響，并且隨著乘性誤差變大，DTM模型的高程也嚴(yán)重偏離真值，因此針對(duì)加乘性混合誤差的DTM模型需要進(jìn)行更深入的研究來豐富現(xiàn)代大地測(cè)量數(shù)據(jù)處理理論。

圖4 不含誤差時(shí)的數(shù)字地形模型Fig.4 The digital terrain model without random errors

圖5 乘性誤差σm=0.005時(shí)的數(shù)字地形模型Fig.5 The digital terrain model in the case σm=0.005

圖6 乘性誤差σm=0.05時(shí)的數(shù)字地面模型Fig.6 The digital terrain model in the case σm=0.05

圖7 乘性誤差σm=0.1時(shí)的數(shù)字地面模型Fig.7 The digital terra in model in the case σm=0.1

表5 算例2中LS法、WLS法、bcWLS法及SUT法求得的參數(shù)估值、估值與真值之間的二范數(shù)和單位權(quán)中誤差

表6 算例2中LS法、WLS法、bcWLS法及SUT法求得的參數(shù)估值的標(biāo)準(zhǔn)差

通過比較表5中4種算法求得的參數(shù)估值與真值之間的二范數(shù)和單位權(quán)中誤差結(jié)果可以看出，相比于LS法，WLS法、bcWLS法及SUT法求得的參數(shù)估值更接近于參數(shù)真值，且單位權(quán)中誤差也更接近真值，表明考慮觀測(cè)值的權(quán)能夠得到較優(yōu)的結(jié)果。同時(shí)，相比于WLS法，bcWLS法能夠在一定程度上減小參數(shù)估值的偏差得到較優(yōu)的參數(shù)估值，表明文獻(xiàn)[5]中方法的有效性。通過對(duì)比bcWLS法和SUT法的結(jié)果可以看出，SUT法能夠得到與文獻(xiàn)[5]中方法相同的參數(shù)估值結(jié)果，這與算例1得到的結(jié)論一致，進(jìn)一步表明本文通過SUT法求得的參數(shù)估值能達(dá)到二階精度。此外，對(duì)比不同乘性誤差求得的結(jié)果可以發(fā)現(xiàn)，當(dāng)乘性誤差σm=0.005時(shí)，4種算法所得參數(shù)估值與真值較為接近；當(dāng)乘性誤差σm=0.05和σm=0.1時(shí)，4種算法所得參數(shù)估值逐漸偏離參數(shù)真值，表明乘性誤差的大小會(huì)影響所求得的參數(shù)估值結(jié)果。

通過比較表6中4種算法求得的參數(shù)估值的標(biāo)準(zhǔn)差可以看出，LS法求得的結(jié)果最大，WLS法求得的結(jié)果遠(yuǎn)小于LS法，bcWLS法求得的結(jié)果略小于WLS法的結(jié)果，而SUT法求得的結(jié)果最小，精度最高，這與算例1得到的結(jié)論一致，進(jìn)一步表明文獻(xiàn)[5]中方法只能求得參數(shù)估值的一階精度，而本文方法能求得參數(shù)估值的二階精度。同時(shí)，對(duì)比不同乘性誤差求得的結(jié)果可以發(fā)現(xiàn)，隨著乘性誤差變大，4種算法所得的參數(shù)估值的精度也逐漸變差，但SUT法總能求得最優(yōu)的參數(shù)估值的精度信息，進(jìn)一步表明將SUT法與加乘性混合誤差模型結(jié)合是可行且有效的。

3.3 算例3

為進(jìn)一步驗(yàn)證本文方法在實(shí)測(cè)數(shù)據(jù)中的適用性及優(yōu)勢(shì)，算例3通過在網(wǎng)站http:∥srtm.csi.cgiar.org/srtmdata/提供的LiDAR型DTM模型數(shù)據(jù)下載了塔里木盆地864個(gè)點(diǎn)的三維坐標(biāo)數(shù)據(jù)進(jìn)行解算[7,29]。同算例2一樣，受加乘性混合誤差干擾的DTM模型觀測(cè)方程見式(42)，不同的是DTM模型函數(shù)，算例3采用文獻(xiàn)[7,9,29,41]中的DTM模型函數(shù)擬合為

h(x,y)=β1+β2x+β3y+β4xy+β5x2+β6y2

(45)

式中， (x,y)表示地面點(diǎn)坐標(biāo)；h(x,y)表示DTM模型高程值；βi表示未知待求參數(shù)。

目前，由于已有關(guān)于加乘性混合誤差模型的研究成果較少，表現(xiàn)為加乘性混合誤差特征的觀測(cè)值未得到足夠重視，且實(shí)際觀測(cè)數(shù)據(jù)中加性誤差和乘性誤差的分析是很復(fù)雜的，需要進(jìn)行專門的論文研究，因此本算例暫不討論分離實(shí)際觀測(cè)數(shù)據(jù)中的加性誤差和乘性誤差，而是將獲取的實(shí)測(cè)數(shù)據(jù)視為真值通過模擬誤差進(jìn)行解算。為了模擬出實(shí)際地形，本算例將乘性誤差的標(biāo)準(zhǔn)差設(shè)置為0.002，加性誤差的標(biāo)準(zhǔn)差設(shè)置為σa=0.15 m，單位權(quán)中誤差同算例1和2一樣。將不含隨機(jī)誤差時(shí)的DTM模型繪于圖8中，受到加乘性混合誤差干擾的DTM模型繪于圖9中。由算例1和2的結(jié)果可知，LS法和WLS法的結(jié)果不足以反映參數(shù)估值及其真實(shí)精度，因此算例3不再使用LS法和WLS法求解參數(shù)估值及其精度信息，僅通過bcWLS法與SUT法進(jìn)行對(duì)比來驗(yàn)證本文方法的可行性及優(yōu)勢(shì)。針對(duì)求解過程中出現(xiàn)的病態(tài)奇異情況，采用嶺估計(jì)法[42-43]來解決法矩陣求逆不穩(wěn)定的問題，所得的參數(shù)估值的精度信息被列于表7，擬合得到的DTM模型如圖10和圖11所示。在本算例中，將獲取的DTM模型高程真值與模型反演值之間的均方根誤差(root mean square error,RMSE)的大小來比較兩種方法的優(yōu)勢(shì)，具體可表示為

(46)

式中，n表示觀測(cè)值個(gè)數(shù)；htrue,i表示獲取的DTM模型高程真值；cmod,i表示模型反演的高程值。

圖8 不含誤差時(shí)的數(shù)字地形模型Fig.8 The digital terrain model without random errors

圖9 受加乘性混合誤差干擾后的數(shù)字地形模型Fig.9 The digital terrain model disturbed by mixed additive and multiplicative random error

圖10 不含誤差時(shí)的數(shù)字地形模型(點(diǎn))及bcWLS法求得的數(shù)字地面模型(圈)Fig.10 The digital terrain model without random errors (point) and the digital terrain model obtained by the bcWLS method (circle)

由圖8可以看出，本文獲取的DTM模型為一個(gè)陡峭的曲面，高程最高點(diǎn)和最低點(diǎn)的偏差為18.892 8 m；同時(shí)，通過對(duì)比圖8和圖9中的DTM模型可以看出，盡管本算例加入的乘性誤差的標(biāo)準(zhǔn)差僅為0.002，但其對(duì)DTM模型的高程產(chǎn)生了明顯的偏差。此外，由圖10和圖11可以看出，兩種方法均能較好地?cái)M合原始的DTM模型，但通過仔細(xì)對(duì)比也可以發(fā)現(xiàn)，SUT法較bcWLS法能夠更好地?cái)M合原始的DTM模型，表明本文方法在實(shí)測(cè)數(shù)據(jù)處理中的優(yōu)勢(shì)。

圖11 不含誤差時(shí)的數(shù)字地形模型(點(diǎn))及SUT法求得的數(shù)字地面模型(圈)Fig.11 The digital terrain model without random errors (point) and the digital terrainmodel obtained by the SUT method (circle)

對(duì)比表7中的均方根誤差可以看出，SUT法求得的結(jié)果較bcWLS法在數(shù)值上更小，表明本文算法能夠更好地?cái)M合原始的DTM模型。同時(shí)，對(duì)比參數(shù)估值的標(biāo)準(zhǔn)差可以看出，兩種方法均能求得較優(yōu)的參數(shù)估值的精度信息，但由于SUT法為二階精度評(píng)定方法，因此求得的參數(shù)估值的精度明顯優(yōu)于bcWLS法，進(jìn)一步表明本文方法在實(shí)測(cè)數(shù)據(jù)處理中的優(yōu)勢(shì)。

表7 算例3中bcWLS法和SUT法求得的參數(shù)估值的標(biāo)準(zhǔn)差及均方根誤差

4 結(jié) 論

本文首先闡述了加乘性混合誤差模型的統(tǒng)計(jì)性質(zhì)及其誤差結(jié)構(gòu)，然后分析了已有加乘性混合誤差模型的參數(shù)估計(jì)方法能達(dá)到二階無偏，但精度評(píng)定方法只能達(dá)到一階精度，且無法通過近似函數(shù)法來獲取參數(shù)估值的二階精度信息，最后通過結(jié)合一種無須求導(dǎo)、無須了解非線性函數(shù)構(gòu)成的SUT法對(duì)加乘性混合誤差模型進(jìn)行精度評(píng)定，并設(shè)計(jì)了具體的算法流程。

通過對(duì)本文算例中的各方法對(duì)比分析得出，當(dāng)乘性誤差較小時(shí)，本文算法求得的參數(shù)估值的二階精度信息略優(yōu)于已有方法求得的參數(shù)估值的一階精度信息；而當(dāng)乘性誤差較大時(shí)，此時(shí)通過SUT法求得的參數(shù)估值的二階精度信息明顯優(yōu)于已有方法；此外，將SUT法應(yīng)用于實(shí)測(cè)DTM模型中，可以得到比已有方法更小的均方根誤差結(jié)果，且參數(shù)估值的精度更高，表明將SUT方法用于解決加乘性混合誤差模型的精度評(píng)定問題是可行且有效的。

本文僅是在理論方面對(duì)加乘性混合誤差模型精度評(píng)定問題進(jìn)行探討，并從理論上針對(duì)加乘性混合誤差模型已有解法的不足給出了一種新的算法，為今后研究LiDAR數(shù)據(jù)、SAR圖像等[44-45]實(shí)際觀測(cè)數(shù)據(jù)的乘性誤差特性，以及結(jié)合地空觀測(cè)的實(shí)踐與應(yīng)用提供堅(jiān)實(shí)的理論基礎(chǔ)，而如何將本文方法與方差分量估計(jì)法[46]結(jié)合進(jìn)而切實(shí)運(yùn)用于實(shí)際中是接下來要研究的內(nèi)容。

致謝：感謝許光煜和鄒傳義提出的寶貴建議和幫助。