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

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

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

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

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

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

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

但是,由于無(wú)控制DEM匹配是一個(gè)復(fù)雜的多維、多峰最優(yōu)化問(wèn)題,而最小二乘法是一種局部最優(yōu)搜索策略。傳統(tǒng)基于整體點(diǎn)的匹配方法收斂效果要受匹配初始狀態(tài)的影響,當(dāng)待匹配DEM相對(duì)于基準(zhǔn)DEM變形較大時(shí),很容易收斂到錯(cuò)誤的局部最優(yōu)解甚至不收斂[2],因此,有必要尋求一種方法來(lái)拉大基于整體點(diǎn)的匹配方法的收斂范圍,以此來(lái)提高該類方法的適用范圍和自動(dòng)化程度。

遺傳算法(GA,genetic algorithms)是一種常用的全局最優(yōu)化方法,該方法具有自組織、自適應(yīng)、自學(xué)習(xí)性和本質(zhì)并行性等特點(diǎn),對(duì)求解多維、多峰、全局優(yōu)化問(wèn)題相對(duì)于最小二乘法等傳統(tǒng)優(yōu)化算法有明顯的優(yōu)越性[15]。

本文擬研究基于GA和 LSM相結(jié)合的無(wú)控制DEM匹配方法,利用GA的全局最優(yōu)搜索能力,為L(zhǎng)SM提供合適的初始解,并充分發(fā)揮LSM精度高、速度快的優(yōu)點(diǎn)。

1 遺傳算法

GA是由美國(guó)Michigan大學(xué)的John Holland與其同事、學(xué)生們?cè)?0世紀(jì)60年代末期到70年代初期研究形成的一個(gè)較完整的理論和方法,主要通過(guò)模擬生物進(jìn)化的機(jī)制來(lái)構(gòu)造人工系統(tǒng)的模型[15],其基本算法流程如圖1所示。目前,針對(duì)GA中的各個(gè)步驟,均有不同的改進(jìn),應(yīng)用中可以根據(jù)實(shí)際需要進(jìn)行選擇。

圖1 遺傳算法基本流程式

2 新匹配方法

2.1 匹配模型

傳統(tǒng)整體點(diǎn)無(wú)控制DEM匹配方法一般基于以下模型[2]:

式中:wi——取值0或1,用來(lái)處理DEM 表面沒(méi)有覆蓋相同區(qū)域的問(wèn)題=[pix,piy,piz]T——待匹配DEM 上一點(diǎn)的坐標(biāo)向量=[qix,qiy,qiz]——通過(guò)某種規(guī)則在基準(zhǔn)DEM 上建立的的對(duì)應(yīng)點(diǎn)坐標(biāo)向量 ;s——縮放系數(shù);R——旋轉(zhuǎn)矩陣=[tx,ty,]T——平移向量,i=1,2,…,N,N為基準(zhǔn)模型包含的點(diǎn)數(shù)。目標(biāo)函數(shù)的實(shí)際意義為基準(zhǔn)DEM與待匹配DEM對(duì)應(yīng)點(diǎn)之間的距離平方和最小。

由于GA是一種全局最優(yōu)化方法,采用以上傳統(tǒng)匹配模型會(huì)出現(xiàn)以下2個(gè)問(wèn)題:

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

(2)由于目標(biāo)函數(shù)連續(xù),所以對(duì)匹配效果敏感,收斂速度快,但容易收斂到局部極值點(diǎn)。

為了解決這2個(gè)問(wèn)題,設(shè)計(jì)匹配模型:

式中:N——基準(zhǔn)模型點(diǎn)數(shù);d(i)——基準(zhǔn)模型上第i點(diǎn)與待匹配模型對(duì)應(yīng)點(diǎn)之間接近程度的量,由以下函數(shù)給出

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

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

2.2 匹配方法

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

圖2 本文方法流程圖

從圖2可看出,本文方法是以LSM作為GA迭代的終止條件,而GA的作用是為L(zhǎng)SM提供合適的初始匹配狀態(tài)。其中,LSM計(jì)算不是在GA每代搜索結(jié)束后,而是在GA連續(xù)搜索Nu代后才進(jìn)行,目的在于避免GA沒(méi)有搜索到全局最優(yōu)解附近而頻繁進(jìn)行LSM計(jì)算浪費(fèi)時(shí)間。

3 實(shí)驗(yàn)分析

3.1 實(shí)驗(yàn)準(zhǔn)備

為了驗(yàn)證本文方法的有效性,設(shè)計(jì)了以下2組實(shí)驗(yàn)。第1組實(shí)驗(yàn)采用模擬數(shù)據(jù),原因在于模擬數(shù)據(jù)匹配前后基準(zhǔn)DEM和待匹配DEM數(shù)據(jù)之間的點(diǎn)對(duì)應(yīng)關(guān)系為已知,能夠有效、準(zhǔn)確評(píng)價(jià)匹配效果。實(shí)驗(yàn)?zāi)康脑谟趯?duì)本文提出的模型和方法進(jìn)行性能指標(biāo)測(cè)試,并與傳統(tǒng)模型和方法進(jìn)行對(duì)比研究。第2組實(shí)驗(yàn)采用某泥石流溝區(qū)不同時(shí)期的實(shí)際DEM數(shù)據(jù),通過(guò)對(duì)匹配前后提取的山脊(谷)線匹配情況的分析,進(jìn)一步驗(yàn)證本文方法的有效性。

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

圖3 實(shí)驗(yàn)數(shù)據(jù)

采用MathWorks公司的MATLAB 7.01軟件作為試驗(yàn)程序開(kāi)發(fā)和執(zhí)行平臺(tái),操作系統(tǒng)為MS Windows XP,基本硬件條件為 Intel Core Duo T2300雙核處理器,2.00GB內(nèi)存,1.66GHz主頻。

第1組實(shí)驗(yàn)性能評(píng)價(jià)指標(biāo)包括速度、精度和拉入范圍。其中,速度指標(biāo)用完成一次匹配所花費(fèi)時(shí)間(time)衡量;精度指標(biāo)用匹配完成后基準(zhǔn)DEM與待匹配DEM之間對(duì)應(yīng)點(diǎn)距離的均方根誤差(RMSE)衡量;拉入范圍指標(biāo)為相應(yīng)方法能夠正確完成匹配的情況下,待匹配DEM相對(duì)于基準(zhǔn)DEM最大的變形參數(shù)。

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

3.2 實(shí)驗(yàn)設(shè)計(jì)

第1組實(shí)驗(yàn)分為以下幾個(gè)部分進(jìn)行。

(1)基于傳統(tǒng)模型(式(1))的GA匹配方法性能測(cè)試。模擬變形參數(shù)為 φ=ω=κ=5k(°),s=0.2k+1,tx=ty=tz=k+1(格)。當(dāng)k=1,2,…,10時(shí),100次實(shí)驗(yàn)結(jié)果表明,傳統(tǒng)模型無(wú)法避開(kāi)收斂到s=0的錯(cuò)誤全局最優(yōu)極值處(或其附近)。

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

結(jié)果顯示,ICP算法在k=1時(shí),RMSE為9.23格,已經(jīng)發(fā)散;LZD算法在k＜23時(shí)收斂,最快3.98 s,最慢40.02 s,平均17.15 s,且收斂時(shí)間是隨k的增大而增加;RMSE不超過(guò)1.01格。

(3)基于本文模型式(2)的GA匹配方法性能測(cè)試。模擬變形參數(shù)同1,實(shí)驗(yàn)結(jié)果見(jiàn)表1。

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

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

圖4 GA+LZD方法收斂時(shí)間

3.3 實(shí)驗(yàn)分析

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

圖5 實(shí)際DEM數(shù)據(jù)特征線匹配效果

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

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

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

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

4 結(jié)論

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

將GA和LZD算法相結(jié)合的無(wú)控制DEM匹配方法,融合了GA拉入范圍大和LZD算法收斂精度高的優(yōu)點(diǎn),并且GA能夠以較快的速度收斂到全局最優(yōu)解附近,為L(zhǎng)ZD算法提供合適的初始匹配狀態(tài),因此也有較高的收斂效率。

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

從以上仿真實(shí)驗(yàn)和實(shí)際DEM數(shù)據(jù)實(shí)驗(yàn)結(jié)果容易看出,本文方法能夠準(zhǔn)確、快速地將待匹配DEM和基準(zhǔn)DEM的統(tǒng)一到同一坐標(biāo)系下,受噪聲(較小變形)影響較小。因此,本文方法可用于泥石流災(zāi)害監(jiān)測(cè)、土壤侵蝕、滑坡變形預(yù)報(bào)等領(lǐng)域,這也是本文后續(xù)工作的重點(diǎn)。

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

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

[3] 張同剛,岑敏儀,吳興華.基于差分模型的無(wú)控制DEM差異探測(cè)方法[J].西南交通大學(xué)學(xué)報(bào),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] 張同剛,岑敏儀,秦軍,等.多時(shí)相DEM 匹配探測(cè)泥石流地表變形的新方法[J].水土保持通報(bào),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] 熊興華,錢曾波,王任享.遺傳算法與最小二乘法相結(jié)合的遙感圖像子像素匹配[J].測(cè)繪學(xué)報(bào),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] 張同剛,岑敏儀,吳興華.無(wú)控制DEM匹配的最小法向距離算法[J].自然科學(xué)進(jìn)展,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] 王小平,曹立明.遺傳算法:理論、應(yīng)用與軟件實(shí)現(xiàn)[M].西安:西安交通大學(xué)出版社,2002:20-49.

[16] 湯國(guó)安,楊昕.ArcGIS地理信息系統(tǒng)空間分析實(shí)驗(yàn)教程[M].北京:科學(xué)出版社,2006:446-452.