漸進性優(yōu)化的尺度自適應Cauchy穩(wěn)健估計模型及其應用

李加元，張永軍，艾明耀，胡慶武

武漢大學遙感信息工程學院，湖北 武漢 430079

穩(wěn)健估計是幾何處理與平差系統中的一個核心基礎問題[1-5]，在攝影測量與遙感、計算機視覺、機器人視覺等領域發(fā)揮著重要作用。它是一種同時進行幾何模型擬合與異常值檢測的技術，能夠從被污染的觀測值中正確估計出觀測值所滿足的幾何模型。由于實際觀測值中不可避免地存在粗差(異常值)，該技術已被廣泛應用于誤匹配剔除[6]、光束法平差[7-8]、點云配準[9-10]、同時定位與制圖(simultaneous localization and mapping,SLAM)[11]等攝影測量任務中?；贛估計的迭代加權最小二乘(iterative reweighted least squares,IRLS)[12-19]和隨機采樣一致性方法(random sample consensus,RANSAC)[20-27]為最常采用的兩類穩(wěn)健估計方法[28]。

IRLS又稱選權迭代法[12-13]，相比RANSAC類方法，其優(yōu)點是運算效率高并且能夠確保局部最優(yōu)解[28]，IRLS通常采用M估計作為穩(wěn)健核函數。在迭代過程中，IRLS的核心思想是基于觀測值殘差賦予粗差小權值，可靠觀測大權值，從而緩解粗差對平差系統的影響。IRLS的關鍵在于權函數的選取，常用的權函數包括Huber、Cauchy、Welsch、Tukey等M估計核[13]。文獻[15]基于驗后方差估計思想推導權函數。文獻[16]提出一種基于IGG(institute of geodesy and geophysics)權函數的迭代加權總體最小二乘法。文獻[18—19]分別從效率與穩(wěn)健性方面對傳統M估計進行改進。IRLS及現有變種都存在一個嚴重缺陷，即IRLS在理論上無法處理高于50%的粗差。一旦觀測值中粗差比率高于50%，此類方法完全失效。隨著自動化處理的普及，觀測值包含50%粗差的情形隨處可見，如影像特征匹配、點云特征匹配、SLAM閉環(huán)檢測等。因此，在這些任務中，對高粗差比率較為穩(wěn)健的RANSAC方法更受青睞。

RANSAC是一種假設檢驗技術[20]，它主要包括兩個步驟：①基于隨機最小子集的幾何模型估計；②基于一致性集的幾何模型驗證。RANSAC交替執(zhí)行上述兩個步驟直至滿足終止迭代條件，對應于最大一致性集的幾何模型被認為是正確解。RANSAC存在大量改進算法[21-27]。MLESAC(maximum-likelihood estimation by random sampling consensus)[21]基于概率統計思想，采用最大化似然性來取代RANSAC方法的最大化一致集大小。LO-RANSAC(locally optimized RANSAC)[22]和FLO-RANSAC(fixing locally optimized RANSAC)[23]通過增加局部優(yōu)化過程來減少所需的采樣次數。GC-RANSAC(graph-cut RANSAC)[24]通過引入圖割技術來進一步提高局部優(yōu)化的可靠性。RECON(residual consensus)[25]利用殘差一致性來消除RANSAC對內點閾值的依賴性。USAC(universal RANSAC)[26]將多種改進算法的特性納入一個通用框架內，包括工程及計算等方面的技巧。MAGSAC++(marginalizing sample consensus)[27]提出一種獨立于內點閾值的打分函數。但是，RANSAC類方法存在3大缺點：①采用最小集估計幾何模型僅能獲得近似解，最優(yōu)解需全部觀測值參與解算；②隨著粗差比率的增加，此類方法的時間復雜度呈指數增長，運行效率低下；③RANSAC類方法難以直接應用于測量平差系統中。

近年來，許多研究工作試圖突破上述M估計的瓶頸問題。例如，M估計被pbM估計器[29]和廣義pbM估計器[30]重新表述為基于投影的優(yōu)化問題，并基于變量誤差(EIV)模型進行求解。加權Q范數[31]估計器進一步提高了Q范數模型(0

本文旨在構建一種同時繼承IRLS的優(yōu)點(高效且最優(yōu))與RANSAC的優(yōu)點(能處理高于50%的粗差)的模型，提出了基于漸進優(yōu)化的尺度自適應Cauchy穩(wěn)健估計模型(注：本文所提尺度自適應模型具有通用性，并不僅限定于Cauchy估計，也適用于Huber、Tukey等M估計)。所提模型引入尺度因子來不斷過濾大殘差觀測值并調節(jié)Cauchy核函數的穩(wěn)健性，采用由粗到精的方式進行漸進優(yōu)化。與前期工作[33]的不同之處主要表現為兩個方面：①本文模型與統一代價函數的數學模型不同。本文模型是M估計模型的推廣，為非分段函數，不存在奇異點問題；②在模型優(yōu)化求解過程中，每次迭代利用控制參數剔除一部分大殘差觀測值，因此，在本文的優(yōu)化過程中，觀測值是動態(tài)變化的。而在統一代價函數模型中，盡管觀測值在某次迭代中獲取了很小的權值，但在后續(xù)迭代中有可能重新獲得較大權值，導致迭代震蕩現象。此外，本文還給出了其在經典攝影測量任務中的應用，包括匹配剔除、后方交會及點云配準。試驗結果表明，該模型明顯優(yōu)于基于M估計的IRLS模型和RANSAC類方法。

1 經典M估計與IRLS

M估計是最大似然估計的泛化形式，其數學表達如下(以線性回歸為例)

(1)

式中，ρ(vi)為穩(wěn)健核函數，vi=(f(xi,θ)-yi)為觀測值(xi,yi)的殘差；f(·)是觀測值所滿足的幾何模型(此處為線性模型)；θ為待求模型參數；N為觀測值數目。

表1 常用的M估計核函數ρ(v)和權函數w(v)(v為觀測值殘差，c為常量)Tab.1 Several commonly used M-estimators and their weight functions (v is the residuals of observations,c is a constant)

式(1)通常轉換為下述IRLS問題進行求解

(2)

(1) 變量估計：給定權值w(vi)t-1作為已知量(t為迭代計數器)，采用加權最小二乘法求解問題(3)

(3)

(3) 交替迭代執(zhí)行步驟(1)和(2)，直至滿足算法終止條件。

2 尺度自適應Cauchy穩(wěn)健估計模型

動機：經典Cauchy估計只能處理低粗差比率問題。如果粗差比率較高(>50%)，則初始估計(首次迭代)模型將大幅偏離真值，從而造成正確觀測值的殘差較大，權重很小。圖1(a)為經典Cauchy估計的權函數，可以看出，當殘差為20時，其權重僅為0.01。在如此苛刻的條件下，許多正確觀測值并沒有實際參與優(yōu)化過程。因此，理想的穩(wěn)健估計代價函數應該首先讓所有觀測值都參與到優(yōu)化過程，然后，隨著迭代的進行逐漸提高代價函數的穩(wěn)健性，最終變?yōu)榻浀涞腗估計函數。具體如圖1(b)所示，觀測值取得小權值(假設為0.01)所對應的殘差閾值應隨著迭代的進行不斷縮小。開始時所有觀測值都參與優(yōu)化，盡管初始模型可能不準確，但殘差較大的正確觀測值仍會對優(yōu)化過程做出較大貢獻。每次迭代后由于閾值減小，一部分具有最大殘差的觀測值將會被賦予小權值，這些具有最大殘差的觀測值通常都是粗差觀測。因此，隨著代價函數的變化，許多粗差觀測被剔除，真實的粗差比率會不斷降低。當代價函數變?yōu)榻浀銶估計時，真實粗差比率將會小于50%。

圖1 傳統Cauchy M估計與理想穩(wěn)健估計的權函數曲線對比Fig.1 Comparison of traditional M-estimation function and the ideal function for robust estimation

模型：為實現上述過程，本文提出了尺度自適應Cauchy穩(wěn)健估計模型。在經典Cauchy估計中引入尺度參數α來控制代價函數的穩(wěn)健性，其數學公式為

(4)

與其對應的權函數為

(5)

圖1(b)顯示了不同尺度下的權函數曲線。可以看到，當α取值越大時，越多的觀測值能夠參與優(yōu)化過程；反之，當α取值越小時，只有殘差很小的觀測值參與優(yōu)化，模型估計精度越高。

求解：所提模型為經典M估計的變種，可采用“由粗到精”的IRLS優(yōu)化策略。首先，將α初始化為大值(本文將初始尺度α0設為初始估計模型的最大殘差)，進行模型估計、粗差剔除和權值更新；然后，隨著迭代過程α不斷減小，并執(zhí)行上述步驟直至模型收斂。

(1)變量估計：給定尺度αt-1和觀測值權值w(vi,αt-1)t-1作為已知量，采用加權最小二乘法求解問題，見式(6)

(6)

(3)尺度更新：減小尺度αt=αt-1/τ，其中τ為步距常數(本文設為1.3)。

(4) 交替迭代執(zhí)行步驟(1)—(3)直至滿足終止條件。

可以看出，所提模型的優(yōu)化過程與傳統IRLS相似，僅多出粗差剔除與尺度更新過程。因而，其也具有傳統IRLS方法的優(yōu)點，如優(yōu)化效率高、具有最優(yōu)解、適用于平差系統等。

3 模型應用

穩(wěn)健估計在攝影測量與遙感任務中具有重要作用，本節(jié)將所提尺度自適應Cauchy模型應用于誤匹配剔除、相機位姿估計及點云配準中。

3.1 誤匹配剔除

yi=f(xi,θ)=Axi+t

(7)

式中，A為2×2仿射矩陣；t為2×1平移向量；θ=(A,t)為仿射變換模型參數。誤匹配剔除的尺度自適應Cauchy權函數為

(8)

3.2 后方交會

(9)

式中，P=K[R,T]為3×4相機投影矩陣；Pk(k=1,2,3)為P的第k行；θ=(R,T)為位姿參數。

后方交會的尺度自適應Cauchy權函數為

(10)

3.3 點云配準

Yi=f(Xi,θ)=RXi+T

(11)

式中，θ=(R,T)為配準參數；R為三維旋轉矩陣；T為三維平移向量。

誤匹配剔除的尺度自適應Cauchy權函數為

(12)

如果在模型f中加入尺度因子，該問題即變?yōu)閿z影測量中經典的絕對定向問題。

4 試驗結果與分析

采用模擬試驗和真實數據試驗來評估本文模型的穩(wěn)健性、精度及效率?；谡`匹配剔除、相機位姿估計及點云配準3個任務，將所提模型與Cauchy估計、Welsch估計、RANSAC、FLO-RANSAC及MAGSAC++進行對比分析。試驗選用均方根誤差(RMSE)、成功率及運行時間作為定量評價指標。每種方法的詳細配置信息(參數與代碼來源)見表2。所有方法的運行時間均在單核1.8 GB赫茲CPU i7-8550U上統計得到。

表2 對比方法參數設置與代碼實現Tab.2 The parameter settings and implement details of each compared method

4.1 誤匹配剔除

試驗結果如圖2所示，Cauchy和Welsch估計在粗差比率低于50%的情況下精度與速度均為最優(yōu)。但是，一旦粗差比率超過50%，其性能將急劇下降，成功率曲線現斷崖式下跌至0。這也與M估計在理論上只能處理低于50%的粗差觀測的結論相吻合。與M估計相比，RANSAC類方法(RANSAC、FLO-RANSAC和MAGSAC++)和本文模型則對高粗差比率較為穩(wěn)健，能處理高達90%的粗差觀測。由于RANSAC類方法無法獲得最優(yōu)解，RANSAC和FLO-RANSAC的模型估計精度相對較差(RMSE最大)。雖然MAGSAC++也是RANSAC的一種變體，但是其在算法中采用了M估計來提升模型精度。本文方法精度與M估計相當，運行效率僅次于M估計。在粗差比率為90%時，本文模型比RANSAC方法快將近3個數量級，比采用C++實現的MAGSAC++快將近2個數量級。

圖2 誤匹配剔除模擬試驗結果Fig.2 Comparison of mismatch removal results on the simulated data

真實試驗：選取6種類型(包括可見光-可見光、紅外-可見光、SAR-可見光、深度圖-可見光、地圖-可見光、夜光-白天)的多模態(tài)影像對進行試驗對比。本文采用RIFT算法提取初始匹配點集。由于多模態(tài)影像存在嚴重的非線性輻射畸變，初始匹配點的粗差比率較高(>80%)，傳統M估計無法成功解算，故只與其他3種方法進行對比，試驗結果見表3。本文模型同時取得了最佳模型精度與運行速度。模型精度相比排名第二的方法提升了15.6%；運行效率比采用C++實現的MAGSAC++快25倍，比RANSAC和FLO-RANSAC快350多倍。本文模型在每種影像對上的定性定量結果如圖3所示。

表3 誤匹配剔除真實試驗結果Tab.3 The results of mismatch removal on real data

注：n為初始匹配數量，γ為粗差比率，黃線表示正確匹配對。圖3 本文方法在多模態(tài)影像上的匹配結果Fig.3 Our feature matching results on the multi-modal image dataset

4.2 后方交會

圖4 后方交會模擬試驗結果Fig.4 Comparison of image orientation results on the simulated data

真實試驗：選用2張由SWDC傾斜相機所攝影像作為測試數據，影像大小為5406×7160像素，攝影比例尺約為1∶8000，相對航高接近700 m。其中，第1幅影像由中心垂直攝影相機拍攝所得，焦距為12 102.1像素；第2幅影像由傾斜攝影相機獲取，焦距為14 671.5像素。采用GPS-RTK技術在第1幅影像覆蓋測區(qū)內布設了36個控制點，第2幅影像測區(qū)布設了60個控制點。然后，由人工刺點方式獲取這些控制點所對應的像點坐標。然而，由于人工刺點錯誤，兩幅影像中正確像點觀測個數分別僅為12和15，即兩幅影像的粗差比率分別為66.7%和75%。實現結果見表4，本文方法再次在RMSE和運行時間兩個指標上取得了最優(yōu)，解算精度約為1個像素。

表4 后方交會真實試驗結果Tab.4 The results of image orientation on real data

4.3 點云配準

模擬試驗：模擬過程與誤匹配剔除類似，不同之處在于點云配準的觀測值為三維坐標點，并且?guī)缀文Ｐ蜑閮H包含三維旋轉和平移的剛體變換。三維點由正態(tài)分布N(0,1002)獲得，旋轉角在區(qū)間[-π/2,π/2]內隨機生成，平移量在區(qū)間[-100,100]內生成，高斯噪聲標準差為0.1 m。試驗結果如圖5所示，M估計(Cauchy和Welsch)在粗差比率超過40%后，性能開始急劇下降。反觀本文方法，即使在90%的粗差比率下，其解算成功率依然接近100%。

圖5 點云配準模擬試驗結果Fig.5 Comparison of point cloud registration results on the simulated data

真實試驗：從ETH公開激光點云數據集中選取3個點云匹配對進行試驗，每個匹配對間的真值剛體變換已知。這些激光點云的數據量非常大(千萬級點數)，因此，首先利用格網濾波方法進行降采樣，采樣間隔為0.1 m。然后，采用ISS(intrinsic shape signatures)[36]算法和FPFH(fast point feature histograms)[37]算法分別進行特征檢測與特征描述，并通過最近鄰匹配得到初始3D特征匹配點。由于激光點云缺乏紋理信息，3D特征描述符的顯著性比較差，所以初始匹配集中粗差比率非常高(本文試驗中>96%)。在進行算法對比試驗之前，首先基于支持線投票策略進行初步篩選，過濾掉部分粗差觀測。試驗結果見表5，可知，本文模型的配準精度達到了0.2 m(兩倍的點云分辨率)。除開預處理部分(降采樣、初始匹配、支持線投票)的耗時外，本文模型僅需0.02 s即能實現配準參數解算。詳細的定性定量結果見圖6。

表5 點云配準真實試驗結果Tab.5 The results of point cloud registration on real data

注：n為初始匹配數量，γ為粗差比率，綠線表示正確匹配對。圖6 本文方法在三維激光點云上的配準結果Fig.6 Our registration results on 3D LiDAR point clouds

4.4 參數τ敏感性分析

基于點云配準模擬試驗來分析步長參數τ的最佳取值，試驗過程與3.3節(jié)類似，區(qū)別在于步長參數τ的取值不同，試驗結果如圖7所示。

由圖7可知，當粗差比率小于80%時，τ在1.1～2.0間取值時均能取得100%的成功率；當粗差比率大于80%時，τ的取值越大，成功率越低。而算法所需的迭代次數與τ成反比，即τ越小，迭代次數越多，計算效率越低。綜合成功率和效率因素，τ的建議取值范圍為1.3～1.6，本文取1.3。

圖7 參數τ敏感性分析試驗結果Fig.7 Results of sensitivity analysis of parameter τ

4.5 模型缺陷

本文模型存在一個隱含假設前提，即粗差觀測近似滿足均勻或者隨機分布(非對抗性粗差，non-adversarial outliers)。如果該假設前提不成立，則本文模型可能失效。比如，當粗差觀測大于50%，并且它們之間也滿足另一個幾何模型時，本文模型的解算結果將為粗差觀測所滿足的幾何模型。因此，其不適用于多幾何模型同時估計的情形。再比如，當正確觀測與粗差觀測在空間分布上各自聚為一團，并且相距較遠時，那么，由于初始估計偏差過大可能導致模型無法收斂。解決該類問題只需提供待求模型的良好初值。

5 結 論

本文提出了一種漸進優(yōu)化的尺度自適應Cauchy穩(wěn)健估計模型。該模型引入了尺度控制參數來調節(jié)經典Cauchy核函數的穩(wěn)健性。此外，該參數還發(fā)揮內點閾值的作用。在優(yōu)化過程中，采用粗到精的迭代加權最小二乘法(IRLS)不斷提升核函數的穩(wěn)健性并過濾掉大殘差觀測值，降低真實粗差比率。采用大量模擬與真實數據對本文模型的正確性、解算精度、運行效率、通用性進行驗證，得到以下幾點結論：

與采用M估計作為核函數的IRLS方法相比，本文模型能處理高粗差比率的情形。與RANSAC類方法相比，本文模型能獲取最優(yōu)解，模型求解精度更高，并且運行效率提升了2～3個數量級。本文模型具有較好的通用性，在誤匹配剔除、后方交會和點云配準任務中均取得了很好的效果。綜上所述，筆者認為該模型具有很好的實用性，在大多數攝影測量、計算機視覺等任務中能夠取代已被廣泛應用的RANSAC方法。

因為本文模型與經典選權迭代法在形式上保持一致，筆者相信本文模型亦能用于平差系統中，進一步提升現有平差系統的抗差性能，這將是筆者下一步的研究重點。此外，IRLS方法比牛頓法的收斂速度慢，在高維度問題上尤為突出，后續(xù)筆者將進一步研究基于牛頓法的穩(wěn)健估計模型[38]，提升平差系統運行效率。