許榮剛

數(shù)字高程模型作為地球表面地形的數(shù)字描述和模擬已成為空間數(shù)據(jù)基礎(chǔ)設(shè)施和“數(shù)字地球”的重要組成部分,目前廣泛用于自動(dòng)繪制等高線、制作坡度、坡向圖、立體透視圖、制作正射影像圖、立體匹配片、立體景觀圖、立體地形模型及地形圖修測等領(lǐng)域[1]。

1 基于LZD算法的DEM匹配基本原理

LZD算法基本思想是：先以兩個(gè)表面上的平面坐標(biāo)相同的點(diǎn)為對應(yīng)點(diǎn)(如果不存在對應(yīng)點(diǎn)就內(nèi)插一個(gè)臨時(shí)點(diǎn)),然后利用對應(yīng)點(diǎn)之間的z坐標(biāo)差(在DEM表面上就是高差)的平方和最小為原則來建立目標(biāo)方程,最后根據(jù)最小二乘原理來求解轉(zhuǎn)換參數(shù)向量,這組參數(shù)能夠拉近兩個(gè)表面。反復(fù)迭代上述過程,就可以正確完成匹配。LZD算法結(jié)合規(guī)則格網(wǎng)DEM自身數(shù)據(jù)結(jié)構(gòu)特點(diǎn)進(jìn)行匹配,適合以規(guī)則格網(wǎng)DEM形式表述的地形表面和其他對象三維表面的匹配[2]。

2 基于LZD算法的DEM匹配方法程序?qū)崿F(xiàn)

1)算法流程分析。

a.建立匹配模型間對應(yīng)關(guān)系：對于待匹配模型中的每一格網(wǎng)點(diǎn)(x,y,z),如果在基準(zhǔn)模型中能找到格網(wǎng)點(diǎn)(X,Y,Z)滿足(X,Y)與(x,y)相等,則這兩點(diǎn)相互對應(yīng);否則,根據(jù)(x,y),在基準(zhǔn)模型中內(nèi)插對應(yīng)點(diǎn)(x,y,Z)。

b.建立觀測方程：

c.求解轉(zhuǎn)換參數(shù)：利用最小二乘法原理,取初值為 ξ=η=α=0,m=1。由n(待匹配模型格網(wǎng)點(diǎn)數(shù))個(gè)觀測方程,建立最小二乘誤差方程,求解七個(gè)未知參數(shù),逐次迭代求解最優(yōu)值。

d.用所求參數(shù),對匹配模型進(jìn)行轉(zhuǎn)換：

2)算法模塊實(shí)現(xiàn)。

a.DEM數(shù)據(jù)處理。由于所給的DEM數(shù)據(jù)是規(guī)則格網(wǎng)數(shù)據(jù),只有高程序列、格網(wǎng)距離和起始坐標(biāo),首先需把一維高程序列轉(zhuǎn)換成二維矩陣。高程矩陣元素下標(biāo)與其X,Y的平面坐標(biāo)相對應(yīng)。

b.方程式的列出和求解。參數(shù)初始值選擇為無平移(X0=Y0=Z0=0)、無旋轉(zhuǎn)(ξ=η=α=0)和比例系數(shù)m=1。也即是說所有計(jì)算的轉(zhuǎn)換信息,將被作為誤差處理,對于精度較差的原始數(shù)據(jù),在此沒有對此方法的靈敏度進(jìn)行詳細(xì)的調(diào)查研究。

DEM中高程已知的這些點(diǎn),能夠列出其觀測方程并采用最小二乘原理通過迭代進(jìn)行計(jì)算,這些離散點(diǎn)的一階導(dǎo)數(shù)也都能夠通過立體模型來計(jì)算[4]。具體流程如圖1所示。

3 LZD算法在DEM匹配中的應(yīng)用研究

1)應(yīng)用實(shí)例的設(shè)計(jì)。

針對LZD方法的DEM匹配問題,進(jìn)行了模擬實(shí)驗(yàn)：

a.實(shí)驗(yàn)數(shù)據(jù)：采用的原始 DEM 大小為160×200,格網(wǎng)間距為xcell：4.522 6,ycell：4.402 5,其X軸,Y軸,Z值單位均為 m。

b.迭代限差：在匹配過程中,迭代終止條件由各參數(shù)相鄰迭代值之差來判定,即七個(gè)轉(zhuǎn)換參數(shù)的增量,七個(gè)轉(zhuǎn)換參數(shù)的限差值由各個(gè)觀測值的精度決定。在本篇論文中,各個(gè)參數(shù)的迭代限差如表1所示。在實(shí)際匹配應(yīng)用中,各個(gè)參數(shù)的限差選取則依靠其對應(yīng)的觀測精度、經(jīng)驗(yàn)值確定,在此不再做進(jìn)一步討論。

表1 參數(shù)限差表

2)實(shí)例設(shè)計(jì)：

a.待匹配模型建立：為了驗(yàn)證LZD算法的實(shí)際匹配能力,首先對三維模型即標(biāo)準(zhǔn)模型進(jìn)行了一定大小的平移、旋轉(zhuǎn)和縮放,來產(chǎn)生具有不同旋轉(zhuǎn)、平移和縮放系數(shù)的變換后模型即目標(biāo)模型。在本次實(shí)驗(yàn)中選擇的參數(shù)如表2所示。

表2 參數(shù)表

把上式兩側(cè)均左乘旋轉(zhuǎn)矩陣R的逆矩陣R-1,根據(jù)剛性旋轉(zhuǎn)矩陣R的性質(zhì)RT=R-1,所以R-1×R=RT×R=I,P′=RT×(P-T)/m＊。

這樣,通過上式可以產(chǎn)生具有不同平移、旋轉(zhuǎn)和縮放參數(shù)的目標(biāo)模型P′,然后按照本文所述方法進(jìn)行DEM匹配,可以直接求出平移參數(shù)、旋轉(zhuǎn)參數(shù)和縮放系數(shù),并可直接驗(yàn)證匹配是否正確。

b.模型插值：利用 MATLAB函數(shù) griddata,基于待匹配DEM,對基準(zhǔn)DEM進(jìn)行內(nèi)插,使標(biāo)準(zhǔn)模型的所有點(diǎn)在待匹配模型中都有點(diǎn)與之相對應(yīng)。之后便可利用最小二乘法進(jìn)行模型間的匹配參數(shù)求解[3]。

3)應(yīng)用效果分析。

LZD方法進(jìn)行DEM匹配的模型差值和計(jì)算效率分析：

迭代次數(shù)：20次。

匹配后兩個(gè)模型間的位置差：0.000 1,-0.028 0,0.000 6,1.000,-0.000,0.000,-0.000。

由上可知,迭代方法能夠收斂,模型能夠成功匹配。由于LZD算法建立對應(yīng)關(guān)系的準(zhǔn)則比較簡單,雖然迭代速度較慢,但因其算法簡單,計(jì)算量小,所以LZD算法是一種高效便捷的方法。

4)誤差分析。

我們在設(shè)計(jì)匹配實(shí)驗(yàn)時(shí),系數(shù)矩陣中包含的誤差是由于以下兩方面原因：首先,在大于兩格網(wǎng)距離的長度時(shí),一階導(dǎo)數(shù)(斜率)取的是離散的近似值;其次,一階導(dǎo)數(shù)計(jì)算時(shí),高程中通常含有誤差,觀測值也會受影響。在DEM插值中,在每次進(jìn)行迭代時(shí),不是在地表做插值而是在一個(gè)數(shù)學(xué)表面做插值。事實(shí)上,從估計(jì)量模型得到的這一額外偏差是很重要的。

在最小二乘法的應(yīng)用中,已假設(shè)觀測值之間相互獨(dú)立,盡管我們沒有在最小二乘估計(jì)中引入先驗(yàn)協(xié)方差陣,忽視在此實(shí)驗(yàn)中的協(xié)方差,可能是正確統(tǒng)計(jì)模型中出現(xiàn)的最大偏差[4]。

4 結(jié)語

由本次試驗(yàn)可得如下結(jié)論：1)LZD算法能夠成功進(jìn)行模型匹配,并具有較大的拉入范圍和較好的精度、速度。2)LZD算法不要求三維模型是規(guī)則格網(wǎng)的DEM,但規(guī)則格網(wǎng)不論在地面DEM中,還是在立體模型的數(shù)據(jù)簡化求導(dǎo)中和在DEM內(nèi)插方面都是一個(gè)優(yōu)勢。3)基于最小二乘原理的LZD算法,對于兩個(gè)待匹配表面,無需任何先驗(yàn)信息,待匹配模型無需任何特征提取。4)本文沒有對該方法的表面差異探測能力進(jìn)行探討,對于較差的初始參數(shù),其收斂性問題有待進(jìn)一步研究。

[1]張祖勛,張劍清.數(shù)字?jǐn)z影測量學(xué)[M].武漢：武漢大學(xué)出版社,1997.

[2]馮義從,岑敏儀.三維自由表面匹配及其應(yīng)用[J].測繪工程,2005,14(3)：36-40.

[3]郝文化.MATLAB圖形圖像處理應(yīng)用教程[M].北京：中國水利水電出版社,2004.

[4]柴登峰,舒 寧,張劍清.動(dòng)態(tài)DEM 匹配方法研究[J].國土資源遙感,2002,51(1)：38-42.

[5]余祖鋒,許才軍.MATLAB在測量中的應(yīng)用[J].東北測繪,2002(4)：43-46.