L1范數探測粗差失效的觀測量識別方法

閆廣峰，岑敏儀

1. 西南交通大學地球科學與環境工程學院，四川 成都 611756; 2. 高速鐵路運營安全空間信息技術國家地方聯合工程實驗室，四川 成都 610031; 3. 高速鐵路線路工程教育部重點實驗室，四川 成都 610031

觀測值中含有粗差時，采用最小二乘法求得的參數估值會嚴重偏離其真值，消除或削弱粗差的影響是獲取精確、可靠參數估值的必要前提。以最小二乘估計為基礎的粗差探測方法，無論是均值漂移模型[1-3]中統計檢驗量的建立，還是方差膨脹模型[4-7]中等價權的構造，都依賴最小二乘殘差或單位權中誤差。然而對最小二乘估計，粗差只能部分地反映在相應的觀測值殘差上，甚至含有粗差的觀測值不一定對應著大的改正數[8-9]，因此容易導致粗差的誤判。實際上，對線性最小二乘估計，最小二乘殘差是各個觀測值的線性函數，在數值上與各觀測值組成的條件方程閉合差相等但符號相反[10]。當觀測值含有粗差時，在相應的條件方程閉合差中，L1范數估計較LS估計更能集中反映粗差，因此探測粗差的能力更強。對L1范數估計在測量數據處理中的應用，國內外學者作了大量的研究工作。文獻[11]根據Bahadur型線性表達式導出了L1范數估計的方差-協方差矩陣，在此基礎上，文獻[12]導出了L1范數估計的巴爾達統計量并討論了其可靠性。文獻[13]采用線性規劃的單純形法進行攝影測量中相對定向的粗差檢測。文獻[14]采用單純形法討論了L1范數估計在測量控制網觀測值粗差探測中的應用。文獻[15]研究了秩虧高斯馬爾可夫模型下的L1范數估計單純形算法。文獻[16]研究了線性平差模型的L1范數估計的遞歸算法。文獻[17]通過對相關觀測值進行去相關化處理，將單純形法用于GNSS基線網的粗差探測。

L1范數估計在測量領域受到廣泛的關注和研究，正是因為它具有良好的抗差性質，不僅求得的未知參數估值少受甚至不受粗差觀測值的影響，而且粗差能夠集中反映在相應的閉合差中，從而有助于粗差的發現與定位。然而，在采用L1范數估計解決粗差探測問題的過程中發現，存在一類觀測值，雖然具有粗差發現和定位能力，但其含有的粗差無論是多大量級都不能準確定位。若測量系統存在L1抗差性失效點(robustness failpoint inL1-norm estimation,RFP-L1)，由于其在平差系統中的位置及是否含有粗差等情況都是未知的，因此無法確保基于L1范數的粗差探測結果的準確性和可靠性，因為只有當RFP-L1不含粗差時結果才會是可靠的，而這無疑會使得粗差定位結果帶有不確定性。對此，RFP-L1的識別及其粗差的探測將是解決該問題的關鍵[18]，因為只有確定測量系統中不存在RFP-L1，或存在時能夠準確識別其位置并對其是否含有粗差作出準確判斷，這種不確定性才能得以消除。考慮到最小二乘的選權迭代法(iteratively reweighed least squares,IRLS)和線性規劃的單純形法等兩種L1范數估計的求解方法中，IRLS方法要依賴最小二乘殘差確定初始迭代權，當初始值偏離真值較大時最終的估計結果將會極不可靠[8,16]，從而會引入新的不確定性因素。單純形法求解過程不需要進行最小二乘運算。因此，本文從觀測方程出發，結合單純形法L1范數平差算法，探索研究測量系統中RFP-L1觀測值的一般特征和識別方法。

1 條件方程與影響系數

平差問題的函數模型為

(1)

式中，L為觀測值向量；B為設計矩陣；X為未知參數向量；d為常數向量。

由式(1)，經矩陣初等變換運算或采用單純形法，可導出條件方程[19-20]

(2)

式中，[L1L2]T為觀測值向量；E為L1對應的系數矩陣且為單位矩陣；-J為L2對應的系數矩陣；W為閉合差向量。

實際上，式(2)中的J為粗差判斷方程的判斷矩陣。令A=E-J，A即為條件方程的系數矩陣。A矩陣中第j列向量Aj元素絕對值和ζj為

(3)

觀測值Lj中若含有大小為Δg的粗差，其對閉合差向量帶來的影響ΔWg為

ΔWg=AjΔg

(4)

從式(4)可以看出，向量Aj實質上就是Lj中粗差在閉合差向量中的投影向量。

設r個條件方程閉合差的絕對值和為Z，則由式(3)、式(4)可知，Lj中的粗差Δg反映在Z上的大小決定于ζj，當式(2)中Lj在L1時，ζj=1，粗差等量反映在Z中；當Lj在L2中時，ζj不一定為1，粗差會被放大或縮小后反映在Z中。ζj反映的是粗差Δg對函數Z影響程度的大小，為敘述方便，本文稱ζj為觀測值Lj的粗差對閉合差絕對值和(Z)的影響系數，簡稱影響系數。

若以殘差絕對值和最小為準則，采用單純形法對測量系統的觀測方程進行求解[14]，根據其解的結構可知，在單純形解中，有t個未知參數和r個觀測值殘差為基變量，剩余的t個觀測值殘差為取值為零的非基變量。這就說明，t個未知參數是由t個非基變量決定的[10]，可以由非基變量對應的觀測方程組求得，而余下的r個多余觀測值殘差可以通過t個未知參數計算得到，目標函數等于r個作為基變量的觀測值殘差的絕對值和[21]。由此可知，采用單純形法可以得到形如式(2)的條件方程，而此時各條件方程閉合差的絕對值和Z實質上就是采用單純形法取得最優解時的目標函數，對于給定的測量系統，要使得Z取得最小值，在偶然誤差服從正態分布時，不考慮粗差相互抵消的情況，Lj中的粗差對Z的貢獻應達到最小，即觀測值Lj的影響系數ζj取各種可能取值中的最小值ζjMIN。

2 最小影響系數與RFP-L1

為方便進一步分析最小影響系數與RFP-L1間的關系，結合式(2)引入RFP-L1，即具有粗差發現與定位能力的觀測值L′，無論其含有多大量級的粗差，單純形法求得的式(2)中，其始終位于L2，粗差無法被準確定位[19]。

下面就非RFP-L1和RFP-L1最小影響系數ζMIN的可能取值情況進行分析。

對觀測值Lj為非RFP-L1的情況，當Lj中的粗差足夠大時，單純形解中，其以多余觀測值的形式出現，從而粗差可以被準確發現和定位。由單純形解構造的條件方程中，Lj在L1中，其影響系數ζj取得最小值，且ζjMIN=1。

對觀測值Lj為RFP-L1的情況，顯然，單純形法得到的條件方程中，Lj必然位于L2。由式(2)構造適用于粗差探測問題的線性規劃的標準型[10]。

(1) 目標函數

(5)

(2) 約束條件

(6)

采取以下的原則選取初始基變量。

容易理解，如此選取的初始基變量V′實際上就是單純形法的最優解形式。

在單純形法求解過程中，檢驗數r為[13]

r=CGG-1N-CN

(7)

式中，CG=[1,1,…,1]和CN=[1,1,…,1]分別為基變量和非基變量對應的目標函數式(5)中的1×r和1×(2n-r)系數向量；G為r×r基變量系數矩陣，由式(6)中E和-E的列向量組成；N為r×(2n-r)非基變量系數矩陣，由式(6)中E、-E、J和-J的列向量組成。

(8)

(9)

(10)

式中，Jj為單純形法求得的條件方程式(2)系數矩陣中粗差觀測值Lj對應的列向量。

由于粗差觀測值Lj的影響系數ζj在由單純形解構造的條件方程中取得最小值ζjMIN，結合式(10)可知，RFP-L1觀測值Lj的最小影響系數

(11)

由以上分析可得，非RFP-L1和RFP-L1的最小影響系數分別具有以下特征：

(1) 非RFP-L1觀測值Lj的影響系數在其為多余觀測值時可取得最小值，最小影響系數ζjMIN=1。

(2) RFP-L1觀測值Lj的影響系數在其為必要觀測值時可取得最小值，最小影響系數ζjMIN<1。

在常見的測量問題中，高程控制網的函數模型(觀測方程)是線性的，其設計矩陣中僅含有±1和0，設計矩陣只與控制點個數、構網形式有關，而與控制點高程和空間位置無關[24]。然而，對測邊網、測角(或測方向)網、邊角網等平面控制網或三維控制網，其函數模型是非線性的，通常采用給定未知參數概略值線性化得到函數模型，其設計矩陣除了與控制點個數、網形結構有關外，還與控制點的空間位置有關。此外，對線性回歸、坐標轉換等問題的函數模型雖然是線性的，但其設計矩陣也同樣與樣本點的個數、空間分布及位置有關，以上這兩類測量系統的線性化(或線性)函數模型中未知參數前的系數不會再僅含有±1和0。由線性代數理論可知，對設計矩陣中僅含有±1、0的函數模型，經過矩陣初等變換得到的判斷矩陣J[20]中也僅會含有±1和0值，由此可以推斷，對具有粗差發現和定位能力的觀測值(判斷矩陣J中其對應列向量非零元素個數不小于2的觀測值)，其位于條件方程式(2)的L2時，該觀測值對應的列向量元素絕對值和始終不小于2，只有位于條件方程式(2)的L1時可取得影響系數的最小值。

綜合以上分析可以得出：

(1) 對觀測值具有粗差發現和定位能力且設計矩陣僅含±1和0的測量系統，所有觀測值的最小影響系數均等于1，觀測值中不存在RFP-L1。

(2) 對設計矩陣除含有±1和0(或不含±1和0)外還含有其他元素的測量系統，觀測值中可能會存在RFP-L1，觀測值Lj(j=1,2,…,n)的最小影響系數ζjMIN<1時，Lj為RFP-L1。

3 RFP-L1的識別

根據函數模型設計矩陣的數值特點不難判斷出測量系統是否會存在RFP-L1。對一個可能存在RFP-L1的測量系統，由于非RFP-L1和RFP-L1的最小影響系數呈現出明顯的差異性，據此也不難判斷觀測值是否為RFP-L1，而其中最為關鍵的是計算觀測值的最小影響系數ζMIN。粗差觀測值Lj的影響系數在利用單純形法得到的條件方程下可以取得最小值，考慮到間接平差函數模型的計算機自動建立過程比較簡便[24]，給出一種從觀測方程出發計算每個觀測值ζMIN的算法。

結合式(1)，設計的算法為：

(1) 對觀測值具有粗差發現和定位能力的測量系統，式(1)中設計矩陣B，如果其僅含有±1和0值，則所有觀測值的最小影響系數等于1，算法結束，否則執行步驟(2)。

(2) 根據式(1)構造線性規劃標準型，取目標函數為殘差絕對值和最小，采用單純形法求解條件方程，得到多余觀測值向量L1和必要觀測值向量L2，令i=1。

(3) 往L2中第i個觀測值L2i對應的約束條件常數項l2i中添加一個大粗差Δg(如1000倍觀測中誤差)。

(4) 采用單純形法得到條件方程，并根據式(3)計算觀測值L2i的最小影響系數ζiMIN。

(5) 如果i=t，算法結束，否則剔除在步驟(3)加入l2i的粗差，令s=i+1，并令i=s，執行步驟(3)。

通過以上步驟，可以求得必要觀測值向量L2中各個觀測值的最小影響系數，由此即可判斷L2中的觀測值L2i(i=1,2，…，t)是否為RFP-L1，若ζiMIN=1，觀測值為非RFP-L1；若ζjMIN<1，觀測值為RFP-L1。由于RFP-L1不會出現在單純形解的多余觀測中，因此多余觀測向量L1中的觀測值不會是RFP-L1、最小影響系數均等于1。

觀測值最小影響系數的計算步驟和RFP-L1判別方法，首先進行單純形求解，將觀測值分為L1和L2兩部分，然后逐個地把L2中的每一個觀測值作為研究對象，從局部來逐一考察并判別每一個觀測值是否為RFP-L1。為闡述方便，稱為局部分析識別法(local analysis identification method,LAIM)。

需要說明的是，前文的公式推導視觀測值為獨立且等權，雖然如此，相關結論及RFP-L1識別算法同樣適用于觀測值獨立但不等權的情況，只是在求解時目標函數采用加權殘差絕對值和最小。對于相關觀測的情況，可以首先采用Cholesky分解法[25]將相關觀測等價地轉換為獨立觀測，然后再利用LAIM進行求解。盡管本文的討論是基于單個粗差假設，不考慮粗差相互抵消的情況，相關結論同樣適用于多維粗差情形。此外，對于有t個必要觀測的測量系統，LAIM共需進行(t+1)次單純形法運算，計算效率主要決定于單純形算法的執行效率。

4 算例分析

為驗證設計矩陣中僅含“±1”和“0”的測量系統觀測值最小影響系數的數值特點，設計了水準網仿真試驗；為了進一步檢驗最小影響系數與RFP-L1觀測值的判別關系，并對LAIM方法的有效性進行評價，采用線性回歸測量數據進行仿真試驗。

4.1 水準網

如圖1所示為一水準網，B點高程已知，其余9個點高程待求，全網共計高差觀測值18個。

圖1 水準網Fig.1 Leveling network

觀測量對應的觀測方程為

Lm=Hj-Hi

(12)

式中，Lm為第m個高差觀測值(m=1,2,…,18)；下標i、j為水準點點名(A～J);H為水準點高程參數值或已知高程。

由式(12)可以看出，水準網的每個觀測方程中，未知參數的系數為常數1、-1或0，與各水準路線兩端水準點高程無關。以L1、L4、L5、L6、L9、L12、L14、L16、L17為必要觀測值，得到水準網的判斷矩陣J，列于表1。

表1 水準網判斷矩陣J

由表1可以看出，水準網的判斷矩陣J中僅包含1、-1和0，并且每個必要觀測值對應的列向量中非零元素個數至少為2，結合圖1的網形特點，每個水準點的自由度均不小于3，容易得出，水準網中的每個觀測值都有粗差發現和定位能力，各種必要觀測、多余觀測組合形式的判斷矩陣J中必要觀測值對應列向量非零元素個數始終≥2。由判斷矩陣J容易得到水準網的條件方程，可以發現，各個條件方程實質上就是獨立閉合環或附合路線，不難理解，對該水準網的各種L1、L2組合形式的條件方程，均反映的是高差觀測值之間簡單的圖形條件關系，在每個條件方程中，觀測值的系數為1、-1或0。因此，這類測量系統中，若觀測值Lm(m=1,2,…,n)具有粗差定位能力，且位于條件方程的L2中，則對應的系數列向量非零元素絕對值和必然≥2，即影響系數ζm≥2；如果位于L1中，則對應的系數列向量非零元素絕對值和恒等于1，即影響系數ζm=1。綜合這兩種情況可知，對水準網中具有粗差定位能力的觀測值Lm，最小影響系數ζmMIN=1。

4.2 線性回歸

表2為線性回歸y=2x+6的模擬數據，其中在y中添加了服從正態分布N(0,0.52)的隨機誤差。

表2 線性回歸模擬數據

對以上的線性回歸，有斜率a、截距b兩個待求參數。要確定a和b，至少在9組觀測值中選擇任意兩組來構成觀測方程進行解算。

采用LAIM方法求取各個觀測值的最小影響系數，列于表3。其中，ζmin為觀測值最小影響系數。

表3 觀測值的最小影響系數

從表3可以看出，L1—L8等8個觀測值的最小影響系數均等于1，而L9的最小影響系數小于1，由此不難判斷，L9為RFP-L1，而其他觀測值不是。為了驗證這一結論，設計隨機粗差試驗：選擇Lm(m=1,2，…，9)中的一個添加粗差，粗差大小介于5～20倍測量中誤差之間(2.5～10)，粗差的正負號隨機，遍歷所有觀測值，每個觀測值添加粗差的試驗進行100次，每次試驗采用以殘差絕對值和最小為目標函數的單純形法進行求解。分別統計各個觀測值的100次隨機粗差試驗中，單純形解中其作為多余觀測值、必要觀測值出現的次數，列于表4。

表4 隨機粗差試驗結果

從表4可以看出，L1—L8中的每個觀測值在添加粗差時，單純形解中，粗差觀測值都出現在多余觀測值中，從而粗差可以完全反映在相應的觀測值殘差上，粗差均可以定位；L9中無論添加小粗差還是大粗差，單純形解中都以必要觀測值出現，這樣粗差就被分配在7個多余觀測值殘差中，從而粗差無法準確定位。由此可以得出，L1—L9中L9為RFP-L1。為了檢驗傳統粗差處理方法是否能夠解決RFP-L1的粗差探測問題，進行以上試驗時，同時采用數據探測法和抗差最小二乘法分別進行粗差檢核和抗差估計，結果表明，對L1—L8，兩種方法可以識別到絕大多數的粗差，但對L9，對介于5～20倍中誤差的粗差，兩種方法均無法發現粗差。進一步地，往L9中加入更大量級的粗差進行試驗，結果表明，當粗差足夠大時，數據探測法可以探測到L9中的粗差，抗差最小二乘法也能夠識別L9中的粗差。雖然兩種傳統的粗差處理方法可以處理L9中較大量級的粗差，但均無法有效解決大小為5～20倍中誤差的粗差探測問題，而這樣量級的粗差卻是粗差探測中更為關注的。

為了進一步分析線性回歸算例中L1—L8為非RFP-L1，而L9為RFP-L1的現象，保持L1—L8各觀測值中的x和y不變，L9中的x依次取值23～49，對應y值取2x+6，并在y中加入服從N(0,0.52)的隨機誤差。采用LAIM方法計算各觀測值的最小影響系數，統計發現對L9的各種取值情況，L1—L8的最小影響系數均為1，L9中x取值與其對應的最小影響系數如圖2所示。同時，對L9的各種取值情況，均往L9的y中添加大小為10的粗差，并采用L1范數估計和LS估計方法求得線性回歸的斜率a和截距b，結果分別如圖3、圖4所示。

圖3 L9位置變化時的斜率Fig.3 Slope corresponding to different positions of L9

圖4 L9位置變化時的截距Fig.4 Intercept corresponding to different positions of L9

由圖2可以看出，L9中x的取值較小時，L9距離其他樣本點較近，樣本點的空間分布相對均勻，所有觀測值的最小影響系數均為1，測量系統中不存在RFP-L1。隨著L9中x取值的增大，L9遠離其他樣本點，其最小影響系數小于1，并逐漸減小，此時單純形解中，L9中的粗差將被縮放后分配在7個多余觀測值殘差中，粗差無法正確定位。而且當L9距離其他樣本點足夠遠時，對多余觀測值殘差進行假設檢驗甚至可能無法發現粗差。從圖3、圖4可以發現，L9中含有粗差，無論L9中x如何取值，LS估計求得的斜率a和截距b均嚴重偏離真實值；當L9不為RFP-L1時，L1范數估計求得的斜率a和截距b不受L9中粗差的影響，但隨著L9中x取值增大而變為RFP-L1，求得的斜率a和截距b嚴重偏離真實值。由此可見，對線性回歸問題，測量系統存在RFP-L1且含有粗差時，無論是LS估計還是L1范數估計，求得的未知參數均會受到粗差的較大影響；樣本點的空間分布是測量系統是否存在RFP-L1的重要影響因素，在采用L1范數估計進行粗差探測分析前，若沒進行RFP-L1的識別分析，粗差探測的結果是不可靠的。

綜合以上兩個算例的試驗分析結果可以得出，觀測值粗差是否能完全反映在條件方程閉合差函數Z中，與L1、L2的組合形式有關，L1范數估計的目標函數對取得最優解時L1、L2的組合形式起到了約束作用，而最小影響系數ζMIN刻畫了觀測值粗差投影到L1范數估計目標函數時變化程度的大小；若觀測值的最小影響系數小于1，則其為RFP-L1；LAIM可以有效識別出測量系統中的RFP-L1；對函數模型的設計矩陣僅含有1、-1和0的測量系統，不存在RFP-L1觀測值。

5 結 論

針對采用L1范數估計進行粗差探測分析時，存在有一類觀測值粗差始終無法準確定位的問題。本文由條件方程，推導出觀測值的影響系數計算式，得到最小影響系數大小與RFP-L1觀測值的關系，根據矩陣初等變換理論，獲得從設計矩陣判斷存在RFP-L1觀測值的一般性規律。通過算例驗證基于最小影響系數識別RFP-L1方法的有效性，并得到如下結論：

(1) 采用L1范數估計探測粗差時，若存在RFP-L1觀測值，則粗差會被分配到其他觀測值的殘差中，粗差將不能準確定位。因此，利用L1范數估計進行粗差檢驗時應先進行觀測值的RFP-L1識別，否則粗差探測的結果將會是錯誤的。

(2) 最小影響系數直觀地描述了觀測值粗差對L1范數估計目標函數的影響關系，在殘差絕對值和最小的準則下，非RFP-L1的最小影響系數為1，RFP-L1的最小影響系數小于1。

(3) 設計矩陣中僅含有±1和0的測量系統，具有粗差發現和定位能力的觀測值，最小影響系數均等于1，不存在RFP-L1問題。

篇幅所限，有關RFP-L1問題解決方法的研究將另文介紹。