張同剛，岑敏儀，任自珍，金國清

(1. 西南交通大學地球科學與環(huán)境工程學院，四川 成都 610031； 2. 西南交通大學高速鐵路運營安全空間信息技術國家地方聯(lián)合工程實驗室，四川 成都 610031； 3. 中鐵第五勘測設計院集團有限公司，北京 100000)

?

M-LZD算法中表面變形的可探測性分析

張同剛1，2，岑敏儀1，2，任自珍1，2，金國清3

(1. 西南交通大學地球科學與環(huán)境工程學院，四川 成都 610031； 2. 西南交通大學高速鐵路運營安全空間信息技術國家地方聯(lián)合工程實驗室，四川 成都 610031； 3. 中鐵第五勘測設計院集團有限公司，北京 100000)

摘要：現有的基于LZD方法的DEM變形探測方法都是將表面變形作為粗差，采用不同的穩(wěn)健估計方法來探測。粗差具有發(fā)現和定位能力是粗差探測的前提，而現有方法均沒有進行粗差發(fā)現和定位能力分析。本文首先采用判斷矩陣對其進行了理論分析，然后對結果在不同類型地面的DEM進行了試驗。根據試驗結果，任意2個判斷向量線性無關，判斷矩陣中沒有零元素，相關數等于多余觀測個數。試驗結果表明，DEM表面上無論多少粗差均具有發(fā)現和定位能力；基于LZD方法的DEM表面變形探測方法在理論上是完備的，可以正確地發(fā)現和定位粗差。

關鍵詞：粗差發(fā)現；粗差定位；判斷矩陣；相關數；LZD；DEM

無控制點DEM匹配差異探測是當前地學領域內的熱點和難點問題。隨著地面和機載激光掃描技術的發(fā)展，可以定期或根據需要對興趣地點進行重復掃描，獲取高分辨率的地面三維模型[1]。提取的多時相的DEM包含了豐富的地面變化信息，這是將其應用于滑坡等固體地球自然災害監(jiān)測、防治的關鍵問題[2-3]。不借助地面控制點完成不同時期DEM的配準來探測地面變形的技術為該問題提供了自動化的解決途徑。

當前領域內最有代表性的DEM配準算法是最小高差算法(least Z-difference，LZD)，該方法采用最小二乘技術，能夠在不借助控制點的情況下完成2期DEM的配準。在此方法的基礎上，為了在有變形的情況下精確完成配準并確定地面變形的范圍和幅度，Pilgrim利用穩(wěn)健估計中的M估計來替代最小二乘，提出M-LZD算法，能夠在不超過16%的情況下完成匹配[4]。更復雜的LMS估計和LTS估計也被引入到該領域，如LMS-LZD[5]和LTS-LZD[6]算法；為了進一步提升算法性能，表面局部不變特征[7]和變形量的空間關系[8]也被引入到該領域。隨著算法性能的提升，算法也越來越復雜。

利用穩(wěn)健估計替代最小二乘算法，在匹配過程

中將表面變形作為粗差來對待，是這類方法共同的基本思想。從這個角度來看，探測表面變形的過程實際上就是匹配過程中觀測量粗差的探測過程。根據測量誤差處理理論，觀測量是否具備粗差發(fā)現和定位能力是任何粗差探測方法的前提[9-12]，對于不具備粗差定位能力的情況，任何方法準確定位其位置的概率不會超過50%[9]。這是關系到這一類基于穩(wěn)健估計的LZD算法探測到的表面變形可靠性的基礎性問題。文獻[13]采用判斷矩陣對DEM中任意一個粗差能否發(fā)現和定位問題進行了一些探索，但DEM變形是一個區(qū)域，與之對應的粗差肯定不止一個。

考慮到LMS估計和LTS估計算法模型非常復雜，本文采用判斷矩陣中的相關數[14]，以這類算法中具有代表性的M-LZD算法為依托，對多個粗差發(fā)現和定位能力進行理論分析和模擬試驗研究。

一、M-LZD算法模型

LZD算法采用剛體轉換模型，采用七參數來描述表面間的轉換關系[15]，即表面轉換關系由3個旋轉參數(θx,θy,θz)、3個平移參數(Δx，Δy，Δz)和縮放系數組成，其算法模型可以表示如下

V=AX+dzw

(1)

(2)

式中，A為系數陣，并且列滿秩；V是改正數向量；X為待求轉換參數；dz為Z坐標差；Fxi,Fyi表示第i點的X和Y方向的坡度，xi,yi表示i點的平面坐標；w為對應觀測量的權，由M估計確定。

為方便描述起見，采用單變量來代替A中的元素，這樣A可以簡寫成

(3)

二、判斷矩陣和相關數

根據LZD算法的系數陣A，通過簡單的矩陣行變換可以得到判斷矩陣J[10]

(4)

(5)

判斷矩陣中6個判斷向量(J1-J6)可以求得6個粗差發(fā)現和定位相關數(以下簡稱“相關數”)(η1～η6)。根據文獻[9]

1. 判斷向量的線性相關性分析

判斷矩陣可以表達成判斷向量的形式

(6)

再令

(7)

將式(5)中各項代入J1—J7，并合并同類項系數，得

(8)

(9)

(10)

(11)

(12)

(13)

從式(8)—式(13)可以看出，判斷向量J1—J6的形式非常類似，因此討論6個判斷向量兩兩之間的線性相關性只需要討論任意2個判斷向量之間的線性相關性即可。以J1、J2為例進行討論。

假設

sJ1+tJ2=0(s、t為常數)

(14)

將式(8)、式(9)代入式(14)，并簡化系數得

(15)

根據式(6)和式(7)將上式展開，得

(16)

兩兩作差，并代入各單變量的具體形式，得

(17)

式(17)成立的條件可以根據系數全部為0或不全部為0分成2種情況來分析。顯然當系數全部為0時，判斷向量J1、J2線性無關；其他情況下，J1、J2線性相關。下面著重討論系數不全為0的情況，由于式(17)是和的形式，其成立的條件非常復雜，很難進行系統(tǒng)完善的分析。但顯然式(18)是式(17)成立的條件中的一個。

(18)

根據式中符號的物理含義，Fx和Fy是DEM表面上X、Y方向上的坡度。式(18)中前2條比較簡單，物理意義很明確，分別是DEM表面為X、Y方向的坡度分別相等，從DEM表面特征來看，也就是該表面在X、Y方向上是規(guī)則曲面。如果DEM確實為規(guī)則曲面，那么LZD數學模型中系數陣A(式(2))列秩虧，式(1)沒有唯一解，這種情況與LZD方法的要求矛盾。式(18)中的其他3個條件的物理含義不太明確，不便于分析，但是前2個條件不成立，式(18)肯定不成立。

對于式(17)中在系數不全為0時，其他成立條件很難通過嚴格的數學公式來進行分析。結合實際DEM特點對此進行討論。式(17)中包括DEM的平面坐標和表面坡度。對于實際的規(guī)則格網DEM而言，一方面DEM表面坡度與地形相關，沒有固定的變化規(guī)律；另一方面，平面坐標X、Y具有很強的規(guī)律。再考慮到式(17)中各個條件相互獨立，這樣式(17)在系數不全為零時通常不成立。

對于實際DEM而言，其包含的數據量一般非常大，這樣式(17)中包括的條件就非常多，由于這些條件相互獨立，因此要使得它們同時成立的情況非常罕見。因此，對于實際DEM匹配而言，其判斷向量J1、J2線性無關。類似的，可以證明J2和J3，J3和J4，J4和J5，J5和J6，J1和J6也線性無關。

2. 粗差發(fā)現和定位相關數分析

根據式(8)—式(13)，判斷矩陣的6個判斷向量(J1-J6)形式一致，對于其對應的6個相關數(η1-η6)分析可以任取其中的一個來進行分析。

(18)

三、試驗

為了對以上理論進行驗證，并兼顧試驗的代表性，選擇各種類型地形的DEM來進行模擬試驗(如圖1所示)。這些DEM是實際地形數據，為規(guī)則格網數據。格網間距為10m，大小為120×100像素，在試驗中作為基準DEM。待匹配DEM通過對基準DEM旋轉和平移后，并添加隨機誤差和變形得到。根據M-LZD算法的探測能力，添加變形為15%。試驗中添加的隨機誤差和粗差均服從正態(tài)分布N(0,σ)。

圖1　試驗模型

試驗中在DEM上均勻選取6個點作為必要觀測量，分別進行計算相應的判斷矩陣和相關數。限于篇幅，表1僅給出了其中一個判斷矩陣及其相關數。所有的試驗結果中，任意2個判斷向量均線性無關，這表明DEM表面存在任意一個粗差都具備發(fā)現和定位能力；在判斷矩陣中沒有出現零元素的情況，相關數均等于多余觀測量的個數，這表明最多可以發(fā)現粗差的個數為多余觀測數的一半，可以定位的粗差比可以發(fā)現粗差個數少一個。

表1　判斷矩陣和相關數

根據上面的試驗結果，在實際應用中只要DEM表面變形的比例不超過50%，均可以發(fā)現和定位。本文6組數據的測試結果也表明，添加的15%的變形均被準確定位，因此M-LZD算法中DEM表面變形具備可探測性，算法確定的變形是可靠的。

四、結論

采用判斷矩陣對基于LZD算法的DEM表面粗差發(fā)現和定位能力進行理論分析，并進行了試驗驗證，結果表明DEM表面具有粗差發(fā)現和定位能力，因此基于LZD方法將表面變形作為粗差來處理的方法在理論上是完備的，可以正確地發(fā)現和定位粗差。

采用各種類型的DEM進行了一系列的模擬試驗，不同DEM上的試驗結果一致，并與理論分析相符。試驗結果表明任意2個判斷向量之間均是線性無關的；判斷矩陣中沒有出現零元素，相關數與多余觀測數相等。

參考文獻：

[1]武永斌, 盧小平, 陳曦東, 等. 機載LiDAR鐵路測繪關鍵技術及應用[J]. 測繪通報, 2015(9): 64-67.

[2]張同剛, 岑敏儀, 馮義從, 等. 采用截尾最小二乘估計的DEM匹配方法[J]. 測繪學報, 2009, 38(2): 144-151.

[3]汪燕麟, 殷義程, 施昆. 地震災區(qū)中地面三維激光掃描測繪技術的應用方法分析[J]. 測繪通報, 2015(6): 65-68.

[4]PILGRIMLJ.RobustEstimationAppliedtoSurfaceMatching[J].ISPRSJournalofPhotogrammetryandRemoteSensing, 1996(51)：243-257.

[5]LIZ,XUZ,CENM,etal.RobustSurfaceMatchingforAutomatedDetectionofLocalDeformationsUsingLeast-Median-of-SquaresEstimator[J].PhotogrammetricEngineering&RemoteSensing, 2001, 67(11): 1283-1292.

[6]ZHANGT,CENM.RobustDEMCo-registrationMethodforTerrainChangesAssessmentUsingLeastTrimmedSquaresEstimator[J].AdvancesinSpaceResearch, 2008, 41(11): 1827-1835.

[7]SHARPGC,LEESW,WEHEDK.ICPRegistrationUsingInvariantFeatures[J].IEEETransactionsonPatternAnalysisandMachineIntelligence, 2002, 24(1): 90-102.

[8]LIUY,WEIB.DevelopingStructuralConstraintsforAccurateRegistrationofOverlappingRangeImages[J].RoboticsandAutonomousSystems, 2004, 47(1): 11-30.

[9]李德仁, 袁修孝. 誤差處理與可靠性理論[M]. 武漢: 武漢大學出版社, 2002: 235-292.

[10]岑敏儀, 卓健成, 李志林, 等. 判斷觀測值粗差能否發(fā)現和定位的一種驗前方法[J]. 測繪學報, 2003, 32(2): 134-138.

[11]岑敏儀, 顧利亞, 李志林, 等. 基于判斷矩陣的觀測量粗差發(fā)現和定位相關性分析[J]. 測繪學報, 2005, 34(1): 24-29.

[12]戴激光, 王新亮, 郭一洋. 一種穩(wěn)健的傾斜航空影像粗差探測方法[J]. 測繪通報, 2014(5): 28-31.

[13]ZHANG T, CEN M, REN Z, et al. Ability to Detect and Locate Gross Errors on DEM Matching Algorithm[J]. International Journal of Digital Earth, 2010, 3(1): 72-82.

[14]顧利亞, 岑敏儀, 李志林. 粗差發(fā)現和定位能力與相關系數的關系[J]. 武漢大學學報(信息科學版)，2005, 30(7): 621-624.

[15]ROSENHOLM D, TORELEG?RD K. Three-dimensional Absolute Orientation of Stereo Models Using Digital Elevation Models[J]. Photogrammetric Engineering & Remote Sensing, 1988, 54(10): 1385-1389.

Analysis of Surface Deformation Detectability Based on M-LZD Algorithm

ZHANG Tonggang，CEN Minyi，REN Zizhen，JIN Guoqing

收稿日期：2015-06-23； 修回日期： 2016-01-04

基金項目：中國鐵路總公司重大科研項目(2014G004-B)；中國鐵路總公司科研試驗任務(Z2013-038-5)；長江學者和創(chuàng)新團隊發(fā)展計劃(IRT13092)

作者簡介：張同剛(1977—)，男，博士，副教授，研究方向為數字攝影測量與模式識別。E-mail：swjtuztg@gmail.com

中圖分類號：P237

文獻標識碼：B

文章編號：0494-0911(2016)06-0010-04

引文格式: 張同剛，岑敏儀，任自珍，等． M-LZD 算法中表面變形的可探測性分析[J]． 測繪通報，2016 ( 6) : 10-13． DOI: 10． 13474 /j． cnki． 11-2246． 2016． 0179．