遺傳算法和最小二乘匹配相結合的DEM匹配方法

楊容浩,岑敏儀,張同剛,楊佳

(1.西南交通大學土木工程學院,四川成都610031;2.成都理工大學地球科學學院測繪工程系,四川成都610059;3.成都理工大學核技術與自動化學院,四川成都610059)

滑坡、泥石流等固體地球自然災害已成為全球性問題,這些地表自然災害通常都伴隨著巨大的地表變化,有效確定地表變化范圍和幅度是進行災害評估、防治決策的重要依據。

數字高程模型(DEM,digital elevation model)是實際地球表面的理想表達工具,可以通過航空攝影測量、衛星遙感、激光掃描、地面測量和地形圖數字化等技術適時、準確獲取,將同一地區不同時相、采用不同方式獲取的DEM進行匹配,是進行地表變形監測的有效手段[1]。

無控制DEM匹配技術是解決DEM匹配中建立控制點成本高,精度和效率低,周期長等問題的有效途徑,已開始用于數字地面信息融合,地表變形探測與分析,三維空間數據的絕對定向等方面[2-10]。目前無控制DEM匹配技術可以分為兩大類,一類是基于特征提取的匹配方法,另一類是基于整體點的匹配方法[2]。

基于特征提取的匹配方法受地面特征(如等高線、地性線等)的提取精度影響較大,匹配精度不穩定,且地面特征的提取需要花費較長的時間,應用受到較大限制。傳統基于整體點的匹配方法是建立在最小二乘匹配(LSM,least squares matching)基礎上的,這類方法具有較高的計算冗余度,匹配精度較高,還可以填補基于特征提取的匹配方法特征信息不明顯對象難于匹配的空白[9-11]。

目前,對基于整體點的匹配方法的研究主要集中在基準DEM與待匹配DEM之間的對應關系建立、提高算法運行效率和差異探測能力[2-5]等方面。如Besl等[12]提出的最近點迭代(lterative closest points,ICP)算法,Rosenholm等[8]提出的最小高差(least Z-difference,LZD)算法和張同剛等[13]提出的最小法向距離(least normal distance,LND)算法等都是對基準DEM與待匹配DEM之間對應關系建立進行的研究;張同剛[6]、Besl[11]等在他們的文獻中分別提出了采用截尾最小二乘估計和K—D樹的策略來提高算法的運行效率;Pilgrim等[8-9]提出的M—LZD和Li[14]等提出的LMS—LZD以及張同剛[6]等提出的對LZD的改進方法等,都是對算法的差異探測能力的研究。

但是,由于無控制DEM匹配是一個復雜的多維、多峰最優化問題,而最小二乘法是一種局部最優搜索策略。傳統基于整體點的匹配方法收斂效果要受匹配初始狀態的影響,當待匹配DEM相對于基準DEM變形較大時,很容易收斂到錯誤的局部最優解甚至不收斂[2],因此,有必要尋求一種方法來拉大基于整體點的匹配方法的收斂范圍,以此來提高該類方法的適用范圍和自動化程度。

遺傳算法(GA,genetic algorithms)是一種常用的全局最優化方法,該方法具有自組織、自適應、自學習性和本質并行性等特點,對求解多維、多峰、全局優化問題相對于最小二乘法等傳統優化算法有明顯的優越性[15]。

本文擬研究基于GA和 LSM相結合的無控制DEM匹配方法,利用GA的全局最優搜索能力,為LSM提供合適的初始解,并充分發揮LSM精度高、速度快的優點。

1 遺傳算法

GA是由美國Michigan大學的John Holland與其同事、學生們在20世紀60年代末期到70年代初期研究形成的一個較完整的理論和方法,主要通過模擬生物進化的機制來構造人工系統的模型[15],其基本算法流程如圖1所示。目前,針對GA中的各個步驟,均有不同的改進,應用中可以根據實際需要進行選擇。

圖1 遺傳算法基本流程式

2 新匹配方法

2.1 匹配模型

傳統整體點無控制DEM匹配方法一般基于以下模型[2]:

式中:wi——取值0或1,用來處理DEM 表面沒有覆蓋相同區域的問題=[pix,piy,piz]T——待匹配DEM 上一點的坐標向量=[qix,qiy,qiz]——通過某種規則在基準DEM 上建立的的對應點坐標向量 ;s——縮放系數;R——旋轉矩陣=[tx,ty,]T——平移向量,i=1,2,…,N,N為基準模型包含的點數。目標函數的實際意義為基準DEM與待匹配DEM對應點之間的距離平方和最小。

由于GA是一種全局最優化方法,采用以上傳統匹配模型會出現以下2個問題:

(1)算法容易搜索到錯誤的全局極小點,即s=0,(為基準DEM 上任一點)或wi全部為0。

(2)由于目標函數連續,所以對匹配效果敏感,收斂速度快,但容易收斂到局部極值點。

為了解決這2個問題,設計匹配模型:

式中:N——基準模型點數;d(i)——基準模型上第i點與待匹配模型對應點之間接近程度的量,由以下函數給出

其中,j=1,2,…,Ng-2,dz(i)=|z基(i)-z待(i)|表示基準模型上第i點與待匹配模型對應點的Z坐標之差的絕對值,其對應關系建立方法與LZD方法相同,對于不能在兩模型上找到相同X,Y坐標的點,其dz(i)=∞;Ng為對基準模型與待匹配模型對應點之間接近程度進行等級劃分的等級數;gr和dg分別為對基準模型與待匹配模型對應點之間接近程度進行等級劃分的劃分標準和對應等級的接近程度表示量。

這里,dbz為對應點Z坐標之差參考量,通過對基準模型以初始預匹配精度值pφ,pω,pκ,ps,px,py,pz為參數進行變換,然后與原始基準模型進行對應點Z坐標差求解,dbz為絕對值最大的Z坐標之差(不包括dz(i)=∞)。其實際意義為:如果待匹配模型與基準模型之間的對應點Z坐標絕對差均小于dbz,則可認為已經收斂到精度要求范圍。

2.2 匹配方法

由于GA是一種具有定向制導的隨機搜索算法[15],其優點在于能夠以較快的速度找到求解問題的近似全局最優解,但在最優解附近會收斂很慢,甚至無法達到較高的精度。而LSM在有合適初始解的情況下,能夠快速、高精度的收斂到最優解。因此,本文設計一種將GA和LSM相結合的無控制DEM匹配方法,算法流程如圖2所示。

圖2 本文方法流程圖

從圖2可看出,本文方法是以LSM作為GA迭代的終止條件,而GA的作用是為LSM提供合適的初始匹配狀態。其中,LSM計算不是在GA每代搜索結束后,而是在GA連續搜索Nu代后才進行,目的在于避免GA沒有搜索到全局最優解附近而頻繁進行LSM計算浪費時間。

3 實驗分析

3.1 實驗準備

為了驗證本文方法的有效性,設計了以下2組實驗。第1組實驗采用模擬數據,原因在于模擬數據匹配前后基準DEM和待匹配DEM數據之間的點對應關系為已知,能夠有效、準確評價匹配效果。實驗目的在于對本文提出的模型和方法進行性能指標測試,并與傳統模型和方法進行對比研究。第2組實驗采用某泥石流溝區不同時期的實際DEM數據,通過對匹配前后提取的山脊(谷)線匹配情況的分析,進一步驗證本文方法的有效性。

2組實驗數據均為格網數據,地面分辨率10 m,Z軸方向坐標轉換為“格”單位,每格10 m。第1組實驗基準DEM為某山區DEM(圖3a),大小為101格 ×108格,待匹配DEM為在基準DEM上截取大小為91格 ×101格,添加均值為0,方差為1的高斯噪聲并進行平移、旋轉、縮放以后生成。第2組實驗基準DEM為某泥石流溝區1987年DEM(圖3b),大小為230格×383格,待匹配DEM為該地區1957年DEM,大小為149格×283格。

圖3 實驗數據

采用MathWorks公司的MATLAB 7.01軟件作為試驗程序開發和執行平臺,操作系統為MS Windows XP,基本硬件條件為 Intel Core Duo T2300雙核處理器,2.00GB內存,1.66GHz主頻。

第1組實驗性能評價指標包括速度、精度和拉入范圍。其中,速度指標用完成一次匹配所花費時間(time)衡量;精度指標用匹配完成后基準DEM與待匹配DEM之間對應點距離的均方根誤差(RMSE)衡量;拉入范圍指標為相應方法能夠正確完成匹配的情況下,待匹配DEM相對于基準DEM最大的變形參數。

結合實際應用需要,將3個平移參數的期望精度(px,py,pz)設置為0.01格,縮放系數的期望精度(ps)為0.01,旋轉角度的期望精度(pφ,pω,pκ)為 1′。遺傳算法種群大小N為100,初始種群范圍tx,ty,tz同為[0,30],s為[0.5,2],φ,ω,κ同為[-30°,30°],接近程度等級劃分級數Ng為10,GA迭代終止條件為達到最大遺傳代數(60代)或者連續10代最優目標值變化小于1%。采用Rosenholm[8]提出的方法對待匹配DEM表面點在基準DEM表面尋找對應點。由于GA是一種隨機搜索算法,并不能保證百分百收斂到全局最優解,因此,每組數據重復測試10次。

3.2 實驗設計

第1組實驗分為以下幾個部分進行。

(1)基于傳統模型(式(1))的GA匹配方法性能測試。模擬變形參數為 φ=ω=κ=5k(°),s=0.2k+1,tx=ty=tz=k+1(格)。當k=1,2,…,10時,100次實驗結果表明,傳統模型無法避開收斂到s=0的錯誤全局最優極值處(或其附近)。

(2)LSM方法性能測試。模擬變形參數為φ=ω=κ=10k(°),s=0.01k+1,tx=ty=tz=0.1k(格),其中,k=1時,匹配前RMSE為 4.96格,k=2時,匹配前RMSE為9.68格。迭代終止條件為連續2次迭代求得的轉換參數之差小于限差(平移量為0.01格,縮放系數為0.001,角度為1″)或者達到最大的迭代次數50。

結果顯示,ICP算法在k=1時,RMSE為9.23格,已經發散;LZD算法在k＜23時收斂,最快3.98 s,最慢40.02 s,平均17.15 s,且收斂時間是隨k的增大而增加;RMSE不超過1.01格。

(3)基于本文模型式(2)的GA匹配方法性能測試。模擬變形參數同1,實驗結果見表1。

表1 基于本文模型的GA匹配方法效果

(4)本文方法匹配性能測試。模擬變形參數同(1),每組參數重復運行10次,參數Nu=3,LZD迭代終止條件與(2)同,但最大迭代次數為20,且只有在其終止時迭代次數小于20的情況下,才認為收斂。結果顯示,在k＜15時,每次都能收斂,RMSE與(2)測試LZD算法一致。收斂時間情況見圖4。第2組實驗為采用本文方法對如圖3b所示實際泥石流溝區間隔30 a的兩期DEM進行匹配實驗,并對匹配前后DEM提取山脊線和山谷線,匹配效果見圖5所示。

圖4 GA+LZD方法收斂時間

3.3 實驗分析

3.3.1 第1組實驗 實驗部分(1)結果易見,基于傳統模型的GA匹配方法容易陷入錯誤的全局最優極值處,與前述理論分析結論一致。實驗部分(2)結果表明,ICP算法不適合用于只有部分重疊區域的3D表面匹配,與文獻[3]結論一致。

圖5 實際DEM數據特征線匹配效果

比較實驗部分(2)和(3)結果,可以看出,基于本文模型的GA匹配方法相較于LZD算法拉入范圍提高明顯,但精度整體不如LZD算法,最佳收斂效果與LZD算法接近,平均收斂效果和最差收斂效果都明顯低于 LZD算法;而收斂最快速度和LZD算法接近,整體上要慢得多。

比較實驗部分(2),(3)和(4)結果,可以看出,本文提出的基于GA和LZD相結合的匹配方法,能夠保持LZD算法高精度和GA拉入范圍大的特點,收斂效率相對于GA匹配方法提高明顯,其平均收斂時間和LZD算法相當,少數情況下最慢收斂時間明顯多于LZD算法,但大多數情況下只是略低于LZD算法,表明該方法有較高的穩定性和收斂效率。

從圖5也能看出,在有噪聲影響的情況下,待匹配DEM與基準DEM之間不具有一致的等高線和特征點(山頂點和谷底點),本文方法收斂時,也能收斂到正確的位置。

3.3.2 第2組實驗 從圖5a,5b均可以看出,57 a和87 a的兩期DEM 在匹配前山谷(脊)線位置差異明顯,平移、旋轉和縮放產生的位置差異占主導地位。而從圖5中匹配完成后的57 a山谷(脊)線和87 a山谷(脊)線匹配情況可以看出,山谷(脊)線在絕大多數地方都能較好的吻合,特別是大多數的山脊(谷)線分叉點(如圖中圓圈標記處(未完全標記))都能夠處于重合位置,表明本文方法能夠較好的完成匹配。對于少量不重合的地方(如圖中方框標記處(未完全標記)),大部分是由于泥石流沖刷和堆積產生,也有部分是由于山脊(谷)線提取的不確定性因素產生[16]。

4 結論

基于GA的無控制DEM匹配,經實驗證實,采用等級劃分的DEM匹配模型,可以解決常規匹配模型易收斂到縮放系數s為0的錯誤全局最優極值處(或其附近)的問題;該方法在拉入范圍方面,相對于常規的LZD算法有明顯的改進,但匹配精度和收斂效率要較LZD算法低。

將GA和LZD算法相結合的無控制DEM匹配方法,融合了GA拉入范圍大和LZD算法收斂精度高的優點,并且GA能夠以較快的速度收斂到全局最優解附近,為LZD算法提供合適的初始匹配狀態,因此也有較高的收斂效率。

由于GA的收斂效果要受初始種群范圍設置的影響,因此,本文方法還不是一種絕對意義上的全局最優化方法,今后的研究要研究具有普適意義的初始種群設置方法。另外,GA的編/解碼策略、遺傳算子的選擇和參與匹配計算的采樣點的選擇等對匹配效率影響也比較大,因此,今后的研究工作中也要對這些方面進行關注。

從以上仿真實驗和實際DEM數據實驗結果容易看出,本文方法能夠準確、快速地將待匹配DEM和基準DEM的統一到同一坐標系下,受噪聲(較小變形)影響較小。因此,本文方法可用于泥石流災害監測、土壤侵蝕、滑坡變形預報等領域,這也是本文后續工作的重點。

[1] 馮義從,岑敏儀,張同剛.用于無控制DEM匹配與差異探測監測泥石流災害的方法研究[J].水土保持通報,2007,27(1):74-77.

[2] 張同剛,岑敏儀,馮義從.用于無控制DEM匹配的LZD和ICP算法的比較[J].中國圖象圖形學報,2006,11(5):714-719.

[3] 張同剛,岑敏儀,吳興華.基于差分模型的無控制DEM差異探測方法[J].西南交通大學學報,2006,41(1):91-96.

[4] Salvi J,Matabosch C,Fofi D,et al.A review of recent range image registration methods with accuracy evaluation[J].Image and Vision Computing,2007,25(5):578-596.

[5] Akca D,Gruen A.Recent Advances in Least Squares 3D Surfaces M atching[C]//Optical 3-D Measurement Techniques VII,Vienna,Austria,2005(II):197-206.

[6] 張同剛,岑敏儀,秦軍,等.多時相DEM 匹配探測泥石流地表變形的新方法[J].水土保持通報,2006,26(3):96-100.

[7] Zhang T,Cen M.Robust DEM co-registration method for terrain changes assessment using least trimmed squares estimator[J].Advances in Space Research,2008,41(11):1827-1835.

[8] Rosenholm D,Torelegard K.Three-dimensional absolute orientation of stereo models using digital elevation models[J].Photogrammetric Engineering and Remote Sensing,1988,54(10):1385-1389.

[9] Pilgrim L J.Robust estimation applied to surface matching[J].ISPRSJournal of Photogrammetry and Remote Sensing,1996,51:243-257.

[10] Pilgrim L J.Surface matching and difference detection without the aid of control points[J].Survey Review,1996,33(259):291-304.

[11] 熊興華,錢曾波,王任享.遺傳算法與最小二乘法相結合的遙感圖像子像素匹配[J].測繪學報,2001,30(1):54-59.

[12] Besl P J,Mckay N D.A Method for Registration of 3-D Shapes[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1992,14(2):239-256.

[13] 張同剛,岑敏儀,吳興華.無控制DEM匹配的最小法向距離算法[J].自然科學進展,2006,16(7):868-873.

[14] Li Z,Xu Z,Cen M,et al.Rubust surface matching for automated detection of local deformations using least-median-of-squares estimator[J].Photogrammetric Engineering and Remote Sensing,2001,67(11):1283-1292.

[15] 王小平,曹立明.遺傳算法:理論、應用與軟件實現[M].西安:西安交通大學出版社,2002:20-49.

[16] 湯國安,楊昕.ArcGIS地理信息系統空間分析實驗教程[M].北京:科學出版社,2006:446-452.