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

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

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

?

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

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

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

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

關(guān)鍵詞：粗差發(fā)現(xiàn)；粗差定位；判斷矩陣；相關(guān)數(shù)；LZD；DEM

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

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

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

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

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

一、M-LZD算法模型

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

V=AX+dzw

(1)

(2)

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

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

(3)

二、判斷矩陣和相關(guān)數(shù)

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

(4)

(5)

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

1. 判斷向量的線性相關(guān)性分析

判斷矩陣可以表達(dá)成判斷向量的形式

(6)

再令

(7)

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

(8)

(9)

(10)

(11)

(12)

(13)

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

假設(shè)

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

(14)

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

(15)

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

(16)

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

(17)

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

(18)

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

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

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

2. 粗差發(fā)現(xiàn)和定位相關(guān)數(shù)分析

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

(18)

三、試驗

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

圖1　試驗?zāi)Ｐ?/p>

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

表1　判斷矩陣和相關(guān)數(shù)

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

四、結(jié)論

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

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

參考文獻(xiàn)：

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

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

[3]汪燕麟, 殷義程, 施昆. 地震災(zāi)區(qū)中地面三維激光掃描測繪技術(shù)的應(yīng)用方法分析[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]. 武漢: 武漢大學(xué)出版社, 2002: 235-292.

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

[11]岑敏儀, 顧利亞, 李志林, 等. 基于判斷矩陣的觀測量粗差發(fā)現(xiàn)和定位相關(guān)性分析[J]. 測繪學(xué)報, 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ā)現(xiàn)和定位能力與相關(guān)系數(shù)的關(guān)系[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版)，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)；中國鐵路總公司科研試驗任務(wù)(Z2013-038-5)；長江學(xué)者和創(chuàng)新團隊發(fā)展計劃(IRT13092)

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

中圖分類號：P237

文獻(xiàn)標(biāo)識碼：B

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

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