一種改進的無控制DEM匹配算法

劉志衛,朱建軍,左廷英,周 璀

(1.中南大學 地球科學與信息物理學院，湖南 長沙 410083;2.中南林業科技大學 理學院，湖南 長沙 410018)

隨著DEM獲取方式的發展，DEM的空間分辨率和時間分辨率越來越高，由于DEM獲取手段和獲取時間的不同，使得同一區域的不同時相DEM覆蓋區域不可能完全重合。傳統的方法要求找到兩組數據之間的多對同名點，通過坐標轉換尋求兩組數據之間的轉換參數，才能使兩組數據匹配到同一坐標系下。但是在高山區和分辨率較低的地圖上，很難獲取準確的地面控制點或同名特征點，因此，基于無控制點的DEM匹配具有廣泛的引用價值。此外，DEM無控制匹配不僅可以確定兩組數據之間的絕對定向參數[3]，將其匹配到同一坐標系下，也能根據不同時相的數據探測到地表一定程度的變形量[4-5]。

DEM無控制匹配最早是由Ebner和Mueller[6-7]提出的，其主要目的是用于立體模型的絕對定向。針對無控制DEM匹配問題，Rosenholm和Torlegard[8]提出的最小高程差(LZD,Least Square Z-axes Differences)算法，是為了尋找一種代替傳統的使用控制點進行絕對定向的方法，且獲得了比傳統的使用控制點方法更高的匹配精度；Zhang等人[5]基于LZD算法結合差分模型，通過對不同點進行自適應加權，實現了DEM表面變形量的自動探測；Karras[9]在LZD算法的基礎上，引入數據探測技術得到了穩健的LZD算法，能夠探測到一定程度的變形。但是，利用LZD算法進行匹配時，很難選取合適的初始轉換參數(包括尺度系數、平移參數和旋轉參數)，這會對算法的收斂速度和計算效率產生很大的影響。

最近鄰點迭代(ICP)算法是由Besl等[10]提出來的，該算法通過兩個點集任意對應點間的距離平方和最小為原則來求解轉換參數(3個旋轉參數和3個平移參數)，使得兩個表面的姿態更加接近。但是，傳統ICP算法的缺點是，只適用于存在明確對應關系的點集之間的定位。此外，該算法也需要消耗大量的計算時間[11-12]，為此許多學者都提出了改進方法。方邵江等人[13]提出了使用加權最小二乘進行配準；袁建英等人[14]提出改進的ICP算法，實現重合區域的快速自動定位，實現不同視下點云的快速精確配準。

依上述可以得出，ICP算法采用搜索距離最小的點作為對應點對，其計算距離比較復雜，這使得算法迭代收斂速度很快，幾次迭代計算就可以很好地接近真值，但是整體計算效率很低；而LZD算法由于建立對應關系的準則較簡單，計算量小，整體計算效率較高，但是對于姿態差異較大的模型迭代次數較多，收斂速度較慢。因此，本文在ICP算法的基礎上，提出一種融合ICP與LZD的改進無控制DEM匹配算法，并在傳統最小二乘估計的基礎上引入差分模型求抗差解，最后通過模擬實驗和實測數據對改進算法進行驗證，實驗結果表明，改進后的算法在配準準確的基礎上，不僅克服了傳統ICP算法的局限性，也解決了LZD匹配算法由于初值選擇問題所導致的收斂效率低的問題，提高了無控制DEM匹配的收斂速度和配準精度。

1 ICP和LZD匹配算法的基本原理

ICP(迭代最近點法)算法主要用于三維點云數據的配準問題，其主要思想是通過一定的方法獲取兩點云數據集間點與點之間的對應關系，并使所有對應點之間的距離最近，重復這個過程直到待匹配模型上所有點均找到對應點為止；然后再利用對應點來求解剛體轉換參數，本質問題就是求解對應點間的坐標轉換參數，使得兩點云集可以統一到同一坐標系的算法，實現點云表達的實物信息的融合，這樣反復迭代就可完成匹配。常用的方法：單位四元數法和SVD正交分解法。

LZD算法的基本思想是：先以兩表面上平面坐標相同的點位對應點(如果不存在就內插一個臨時點)，然后利用對應點之間的Z坐標差(在DEM表面上就是高差)的平方和最小為原則來建立目標方程，最后根據最小二乘原理來求解轉換參數向量，這組參數能夠拉近兩個表面。反復迭代上述過程，就可以正確完成匹配。

從上面的基本原理可以看出，兩個算法匹配的框架基本相似，它們共同的算法流程如下：

1)建立兩個表面上點之間的對應關系

(1)

其中，S={pi}為基準模型，M={qi}為待匹配模型；

2)根據對應關系，建立目標方程：

min∑wi‖pi-qi‖2.

(2)

3)根據不同的參數估計準則求解轉換參數(ICP算法采用單位四元素方法求解，LZD采用最小二乘原理求解)

4)根據求得的轉換參數更新待匹配模型；

5)判斷匹配是否完成，若不滿足條件，重復步驟(1)—(5)，直至滿足迭代條件結束。

從上面算法流程可以看出，LZD與ICP算法的核心差別就在于它們處理表面對象的策略不同，這就導致了它們建立點對應關系的算法不同。ICP算法利用的是三維表面點的空間距離最近建立點對應關系，這使得該算法迭代收斂速度很快，但是由于其建立點對應關系的計算量大，從而導致整體計算效率不高，此外，該算法也有一定的局限性：

1)要求目標數據集和參考數據集要具有明顯的特征，否則最終的匹配結果容易陷入局部最優；

2)目標數據集和參考數據集的對應近似點數要相等；而LZD算法是通過利用內插臨時對應點的方法來避免復雜的搜索過程，其建立的關系較為粗略，計算量較小，但是當待匹配DEM對之間的姿態相差較大時，模型之間轉換參數的初始值很難獲取，所以完成匹配需要的迭代次數較多，嚴重影響了算法的迭代收斂速度和計算效率。

2 融合ICP和LZD的改進無控制DEM匹配算法

2.1 算法的實現

ICP算法的實質是基于最小二乘的最優匹配方法，重復進行“確定對應點集-計算最優剛體變換”過程，直到迭代誤差足夠小或迭代次數大于預設的最大迭代次數為止。ICP算法中經常用到的“點到點”的四元素轉換參數法，其主要的求解步驟如下[15]：

1)計算基準DEM(標準模型)和變形DEM(待匹配模型)中任意對應點之間的距離(Di)(如式3)。

(3)

其中，(Xi1,Yi1,Zi1)為基準DEM中任意點的三維坐標，(Xi2,Yi2,Zi2)為待匹配模型DEM中對應點的三維坐標。

2)利用KD-Tree方法，在基準DEM中搜索與待匹配模型DEM的最近鄰點集，然后構造協方差矩陣。

3)根據2)中的協方差矩陣構造4×4的矩陣Q，然后根據矩陣Q的最大特征值對應的特征向量計算旋轉矩陣參數R，進而求得平移參數T(對初始變換參數進行更新)。

4)利用所求的旋轉矩陣參數和平移參數對待匹配DEM進行更新，重復步驟2～4，直到迭代誤差足夠小(θ<ε)或迭代次數大于預設的最大迭代次數為止。

將ICP匹配算法得到的6個轉換參數(尺度系數除外)作為LZD匹配的初始轉換信息進行匹配，具體流程如圖1所示。

基準模型與待匹配模型之間對應的數學轉換關系(見式4)：

(4)

其中，(XR,YR,ZR)是參考DEM的坐標，(XT,YY,ZT)是對參考數據修改后待匹配DEM的坐標；ΔX,ΔY,ΔZ,S,R分別為兩個DEM之間的平移參數、比例縮放系數和旋轉參數矩陣。

圖1 匹配算法具體流程

利用ICP算法求解參考DEM和經處理后的待匹配DEM之間的變換參數，并通過式(5)對待匹配DEM進行坐標變換。

(5)

其中，(XR1,YR1,ZR1)是經過ICP算法匹配后，與基準DEM接近的第一組近似坐標；(X1,Y1,Z1)是修改后的變形模型DEM的坐標；ΔX1，ΔY1，ΔZ1和R1分別表示基準DEM與處理后的待匹配DEM經ICP算法匹配所得到的平移和旋轉矩陣參數。

將經過ICP算法所得到的6個轉換參數(尺度系數除外)作為LZD匹配的初始參數，再采用LZD算法進行精確匹配運算，并通過式(6)進行兩模型間的轉換。

(6)

其中，(XR2,YR2,ZR2)是經過LZD算法匹配后，與參考DEM對應的第二組近似坐標；ΔX2,ΔY2,ΔZ2,S和R2分別表示參考DEM與經過ICP算法匹配后的近似模型之間的平移參數、尺度系數和旋轉矩陣。

將式(5)所得到的6個轉換參數(尺度因子除外)作為LZD匹配的初始值，然后，通過式(6)迭代求解參考DEM和待匹配DEM之間的7個精確轉換參數。因此，將式(5)、式(6)整理到一起可以得到最終的無控制匹配表達式(如式(7)所示)。

(7)

2.2 抗差最小二乘的引入

LZD匹配算法是通過參考模型與待匹配模型列立條件方程，根據最小二乘準則來求取模型之間轉換參數的最優估計值，估計準則如下：

(8)

而LZD算法通過插值建立的同名特征點可能存在“偽同名點”(由于錯誤測量或遮擋等原因造成)，從最小二乘的估計準則可以看出，在平差過程中，異常值(即“偽同名點”)對殘差平方和的影響較大，從而導致了最小二乘估計失去了對粗差的抵抗能力。文獻[16]指出了最小二乘估計的BP值為1/n，即數據中僅有一個非常極端的異常值會對最后的平差結果產生非常惡劣的影響，故而不穩健。因此，需要引入抗差最小二乘的方法來抑制異常值對參數估計的影響，從而獲得具有抗差性的參數估值。

Zhang等人[5]在LZD算法的基礎上，引入高程差分模型，通過對不同點進行自適應加權，實現了DEM表面變形量的自動探測，可以探測到大于50%的表面形變?；舅悸啡缦拢涸谶M行無控制匹配之前，將參與轉換參數計算的非邊緣高程點及其8鄰域范圍內的高程點按從左到右、自上到下的順序進行排列。

2)將兩模型間任一非邊緣點的高差dZi與其8個相鄰的高差均值dZm做差，得到(ΔdZ=dZi-dZm)，顯然，每一個高程值Zi與ΔdZ是一一對應的。根據新的統計量ΔdZ對每一個觀測值賦予不同的權值，以消除或削弱由于地表變形引起的匹配誤差。規則如下：

(9)

其中，Med為中位數。

根據式(9)，每個觀測值被賦予一個權0或1。只有權為1的觀測值參與匹配，其余觀測值在匹配過程中被剔除。通過上述方法進行確權，還存在很多孤立觀測量，即其本身的權值為1，而其相鄰的8個觀測量的權為0，這是由于觀測量中含有的表面變形被隨機誤差掩蓋，用傳統的最小二乘很難發現，但是可以很容易通過觀測值之間的相互關系發現。通過對式(9)中的權值進行適當調整，將孤立觀測點剔除，以提高最終的配準精度。

3 實驗和結果分析

3.1 模擬實驗

3.1.1 模擬實驗數據的選擇

為了對本文算法進行驗證，試驗數據采用張家界地區30 m分辨率的SRTM DEM數據為基礎來設計模擬試驗數據(如圖2所示)。其優點在于理論上完全正確匹配的情況下，匹配完成后與仿真模擬變換前同名點處處重合，這樣便于對本算法的實際試驗結果進行準確的評價和分析。

1)基準DEM。以張家界地區的SRTM DEM數據為基礎，格網間距為30 m，大小為128像素×128像素的子塊作為基準DEM(如圖2(a)圖所示)。

2)模擬變形。

①對基準DEM在高程上添加0～2 m的隨機誤差；

②然后按表1參數表進行旋轉、平移和縮放，產生與基準DEM對應的具有表面形變的待匹配模型。

表1 模擬轉換參數表

3.1.2 模擬實驗結果和分析

圖2中分別為從張家界SRTM DEM數據中截取的初始DEM數據(稱為基準DEM)和按表1中的轉換參數所得到的具有變形的待匹配模型匹配前的分布圖。分別采用基于四元素的經典ICP、ICP和LZD融合的算法，但是沒有考慮表面形變進行自適應加權和采用差分模型進行自適應加權情況下的本文算法對基準模型和待匹配模型進行無控制匹配。

圖3(a)為利用基于四元素經典ICP的配準效果圖、圖3(b)為ICP和LZD相結合算法的配準效果圖、圖3(c)是在ICP和LZD相結合的基礎上，引入了差分估計模型的配準效果圖。為了評定算法的配準結果，采用算法完成配準的迭代次數、完成迭代時間和平均誤差對比，在理論上完全正確匹配的情況下，匹配完成后與模擬變換前的點與點處處重合，因此，本文選擇將迭代結束后各對應點間配準的殘差平均值即“平均誤差”作為評定標準。

圖2 初始DEM數據和匹配前的分布圖(紅色:基準DEM;藍色:變形DEM)

圖3 不同匹配方法的點云數據(紅色:基準DEM點云數據;藍色:配準后的DEM點云數據)

從圖3中的匹配效果圖可以看出，傳統的經典ICP算法(圖3(b))，對于有尺度變換的DEM數據無法完成正確的匹配，匹配后的平均誤差也較大；但是在ICP的基礎上引入LZD算法(圖3(b))，可以實現有尺度變換的數據之間的匹配，且與本文算法(圖3(c))的匹配效果圖很接近，但是，表2中給出的迭代次數、完成迭代的時間和平均誤差的對比，可以看出，本文算法不僅比傳統的ICP算法得到了更精確的配準結果，并且與單一的LZD匹配算法相比，其配準的收斂速度和精度得到了大幅度的提高。此外，表3給出了不同算法匹配后所得到的模型之間的轉換參數與模擬參數之間的誤差比較，可以看出：經典的ICP算法對存在尺度變換的數據進行匹配時，效果最差；而LZD算法與ICP和LZD相結合的算法最后的匹配精度相當。這是由于這兩種算法匹配精度都是由LZD算法所決定，ICP算法只是起到加快收斂速度的作用。本文是是在最小二乘的基礎上引入了差分模型，所求的是轉換參數的抗差解。雖然本文算法相較于ICP+LZD算法在x，y軸平移參數上精度稍差，但是對高程和三軸旋轉參數的估計較精確，而旋轉參數誤差(乘性誤差)與平移參數誤差(加性誤差)相比， 前者對最后匹配精度的影響更大。因此，即使圖3(a)和圖3(b)在主觀視覺上效果很接近，通過轉換參數誤差和平均誤差的對比，本文算法更優。

本文算法匹配過程中，引入差分模型對高程點進行自適應加權處理，求得轉換參數的抗差解，很大程度上消除或削弱了由于地表形變(隨機誤差)所引起的配準誤差，使配準后的平均誤差由1.6 m提高到了0.7 m。實驗結果表明，本文算法不僅可以顧及兩組數據之間的尺度變換參數，且能夠削弱多時期數據之間的地表形變對配準結果的影響，具有較快的收斂速度和較高的配準精度。

3.2 無控制DEM匹配的應用

3.2.1 實驗數據

為了評價全球DEM數據的精度，通常需要將其與我國現有的地形數據轉換到同一坐標系下，保證數據之間具有相同的數學基礎，而投影轉換所需要的平移、旋轉和縮放系數屬于國家保密數據，僅僅通過現有的商業軟件并不能使兩組數據完全重疊，因此需要通過七參數轉換對預處理后的數據進一步校正。試驗數據為張家界地區的參考DEM數據和2003年獲取美國國家航空航天局(NASA)發布的SRTM1 DEM數據。圖4中分別為參考DEM數據和經過預處理后的SRTM1 DEM數據，格網間距均為30 m，數據大小均為128像素×128像素，且存在一定的重疊區域，由于獲取手段和地形條件問題，很難通過在圖上直接選取同名控制點進行數據之間的投影轉換，因此只能采用無控制DEM匹配方法將兩組數據轉換到同一坐標系下。

表2 DEM配準迭代次數和精度對比

表3 不同匹配方法的參數誤差對比

圖4 張家界試驗區DEM數據

從圖5(b)可以看出，在沒有進行精確配準之前，兩組數據的等高線之間存在著明顯的不重合，兩組數據之間存在著一定的旋轉、平移等轉換關系。為了驗證本文算法的實用性和可靠性，分別采用LZD算法和引入差分模型的本文算法對兩組DEM數據進行無控制匹配。

3.2.2 實驗結果

依據不同匹配算法得到的兩組DEM數據間的7個轉換參數，將SRTM1 DEM自動轉換到與參考DEM數據同一坐標系中，匹配后的兩個DEM的點云效果圖和等高線圖如圖6所示。從匹配后的點云效果圖來看，采用無控制DEM匹配方法基本上可以將兩組數據配準到同一坐標系下，兩種算法的匹配效果很接近；但是通過匹配后的等高線圖可以看出，采用本文算法所得到的兩組數據等高線總體吻合程度優于僅采用LZD算法的匹配結果。此外， SRTM1 DEM是通過雷達方式獲取的，考慮到合成孔徑雷達成像的特點，在地面起伏度較大的地區該DEM存在明顯的異常[17]，如圖6區域A和區域B所示。僅采用傳統LZD算法進行無控制匹配，不僅收斂速度很慢，由于異常數據的影響，無法獲取較精確的匹配結果；引入差分模型后的本文算法，降低了粗差(異常)點對配準結果的影響，從而使經過匹配后的兩組等高線數據重疊度較高，匹配結果更加準確、可靠。從匹配時間上看，LZD匹配算法需要迭代28次，時間消耗136 s，而本文算法只需要迭代8次，時間消耗39 s，收斂速度提高約70%，匹配效率更高。

圖5 參考DEM和SRTM1 DEM匹配前點云和等高線圖

圖6 不同算法匹配后點云效果圖及等高線比較圖

4 結 論

本文改進的無控制DEM匹配算法，首先通過經典ICP算法進行初始配準，并將所獲取的初始轉換參數(尺度參數除外)作為初始值，進一步采用LZD算法進行精確配準，提高了LZD匹配算法的收斂速度；此外，由于傳統最小二乘估計不具有抗差性，本文引入高程差分估計模型，在迭代過程中對不同精度的高程點進行自適應加權處理，消除或削弱了由于地表形變所引起的配準誤差，從而保證了算法的精度。依照實驗結果和對比分析可知，改進后的算法在配準效率、配準精度上較傳統算法都有大幅度的提高。最后，將本文算法應用于SRTM1 DEM與參考數據的無控制匹配實驗中，可以獲得較好的匹配結果，其匹配結果可以為DEM數據質量的評價和多源DEM數據的融合提供較好的基礎數據。

LZD算法是通過最小化三維表面所有對應點間的高差平方和來求解表面轉換參數，隨著DEM數據的全球化和多源數據的發展，為了將本文算法的配準結果應用于大范圍多源、多時相DEM數據的融合中，后續工作將在本文算法的基礎上，結合特征點檢測的思想對大范圍DEM數據之間的配準工作進行研究，并對其做進一步的完善。