趙文豪,閆 利,張?jiān)粕?陳斯煬

(1.武漢大學(xué) 測(cè)繪學(xué)院, 湖北 武漢 430079; 2.國(guó)家基礎(chǔ)地理信息中心, 北京 100830； 3.中南大學(xué) 地球科學(xué)與信息物理學(xué)院, 湖南 長(zhǎng)沙 410083 )

隨著激光掃描技術(shù)和多視影像密集匹配方法的發(fā)展，獲取高密度的點(diǎn)云越來越容易，為了減少冗余和提高數(shù)字高程模型(DEM)處理、網(wǎng)絡(luò)傳輸?shù)葢?yīng)用的效率，有必要對(duì)DEM進(jìn)行綜合以獲取多尺度DEM數(shù)據(jù)庫(kù)[1]。目前的DEM綜合方法按照DEM存儲(chǔ)方式可以大致分為2類：面向三角網(wǎng)(TIN)的綜合方法和面向格網(wǎng)DEM的綜合方法[2]。其中3維Douglas-Pucker (DP) 算法是由費(fèi)立凡等人根據(jù)經(jīng)典2維DP算法發(fā)展而來的DEM綜合算法，最早用于處理離散點(diǎn)，也就是TIN表達(dá)形式的DEM[3-4]。由于原始算法處理海量數(shù)據(jù)時(shí)效率較低，何津在利用3維DP算法綜合規(guī)則格網(wǎng)DEM時(shí)，利用行列掃描的方式來減少計(jì)算量[5]。該方法大幅提高了效率，但可能會(huì)出現(xiàn)漏掉重要地形點(diǎn)的情況。根據(jù)朱雪堅(jiān)等人的研究發(fā)現(xiàn)，基準(zhǔn)面的選擇和掃描方向的選擇對(duì)結(jié)果存在較大影響[6]。Ma等人認(rèn)為高程剖面線可以在一定程度上反映DEM的全局特征，因此提出一種4方向剖面線DP簡(jiǎn)化方法來綜合DEM[7]。該方法生成剖面線時(shí)依賴于VIP特征[8]，若VIP特征保留較少，4個(gè)方向的掃描難以確保掃描到全局的DEM特征。

針對(duì)以上問題，本文提出一種基于多方向剖面線簡(jiǎn)化的DEM綜合方法。首先基于多方向剖面線DP簡(jiǎn)化算法提取DEM關(guān)鍵特征點(diǎn)，然后構(gòu)建三角網(wǎng)再內(nèi)插規(guī)則格網(wǎng)DEM，以達(dá)到DEM綜合的目的。

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

為了驗(yàn)證本文算法，本文采用如圖1所示的兩組由航空激光點(diǎn)云數(shù)據(jù)內(nèi)插的格網(wǎng)DEM進(jìn)行實(shí)驗(yàn)，內(nèi)插過程采用surfer 8.0實(shí)現(xiàn)，選擇克里金插值算法生成間距為2 m格網(wǎng)DEM。如圖1(a)所示的數(shù)據(jù)1，格網(wǎng)大小為1 416像素×943像素，高程范圍從1 125～1 623 m；如圖1(b)所示的數(shù)據(jù)2，格網(wǎng)大小為913像素×555像素，高程范圍是從1 500～2 115 m，兩組數(shù)據(jù)地形起伏較大。

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

2 基于DP算法的DEM簡(jiǎn)化原理及其改進(jìn)

2.1 DP算法簡(jiǎn)介

DP算法，是Douglas和Pucker在1973年提出的一種將二維曲線化簡(jiǎn)的算法[3]。DP算法主要有2種形式，其一為“錨點(diǎn)前進(jìn)法”，另外一種為“分而治之法”，后一種方法效率更高[4]。因此本文采用后種方法，圖2是其原理示意圖。如圖2所示的一條曲線，將這個(gè)序列中的首點(diǎn)A和末點(diǎn)J連成一條直線，然后計(jì)算序列中間各個(gè)點(diǎn)到這條直線的距離di。如果最大距離大于預(yù)先設(shè)置的閾值T(該閾值與簡(jiǎn)化后DEM要保留的精度相關(guān)，閾值越小，保留點(diǎn)越多，本文實(shí)驗(yàn)取經(jīng)驗(yàn)值50)，設(shè)最大距離對(duì)應(yīng)的點(diǎn)為S，然后對(duì)端點(diǎn)AS和SJ之間的點(diǎn)序列進(jìn)行相同迭代操作，如果最大距離小于T則停止。最后將插入的端點(diǎn)與A、J連成折線替代原始曲線，從而起到簡(jiǎn)化的作用。

圖2 二維DP算法原理圖

2.2 基于多方向剖面線的特征提取

文獻(xiàn)[7]以VIP特征為起點(diǎn)，向東南西北4個(gè)方向發(fā)射射線與DEM有效邊界相交形成剖面線，然后進(jìn)行DP簡(jiǎn)化。如果為了更好地覆蓋整個(gè)DEM區(qū)域，最直接的方法是將DEM有效區(qū)域的每行的端點(diǎn)看成一條剖面線，然后采用2維DP算法簡(jiǎn)化，這類似于傳統(tǒng)3維DP算法按照行列掃描的方式。為了獲取更多方向的剖面線，本文不針對(duì)原始DEM直接運(yùn)算，而是圍繞原始DEM中心位置，按照一定角度間隔對(duì)原始DEM進(jìn)行旋轉(zhuǎn)。對(duì)于給定的角度θi，采用式(1)計(jì)算每個(gè)點(diǎn)旋轉(zhuǎn)以后的坐標(biāo)。

(1)

(2)

其中(X,Y)和(X′,Y′)分別是旋轉(zhuǎn)前后的DEM水平坐標(biāo)，(XC,YC)為旋轉(zhuǎn)前DEM中心位置水平坐標(biāo)。

提取特征后，采用式(2)對(duì)特征提取結(jié)果進(jìn)行變換，以便于將特征提取結(jié)果重新納入到原始DEM所在的坐標(biāo)系，也便于融合所有旋轉(zhuǎn)DEM提取的特征結(jié)果。圖3(a)DEM為圖1(a)所示DEM逆時(shí)針旋轉(zhuǎn)15°的結(jié)果。進(jìn)行DEM旋轉(zhuǎn)后，只需要沿著新DEM水平方向生成剖面線即可實(shí)現(xiàn)多方向的剖面線生成。剖面線生成過程中，依次掃描新DEM的每一行，其起始點(diǎn)和終點(diǎn)為DEM有效區(qū)域。圖3(b)為圖3(a) 所示剖面線的高程剖面圖。對(duì)于生成的每一條剖面線，采用如圖2所示的2維DP算法進(jìn)行簡(jiǎn)化以獲得特征點(diǎn)，并保留每個(gè)特征點(diǎn)對(duì)應(yīng)的距離，作為每個(gè)特征點(diǎn)的興趣值R，以用于后續(xù)處理。

圖3 DEM旋轉(zhuǎn)結(jié)果與對(duì)應(yīng)剖面線圖

當(dāng)旋轉(zhuǎn)間隔較小時(shí)，相鄰剖面線獲取的特征點(diǎn)可能重疊，因此本文在完成多方向DP簡(jiǎn)化獲得特征序列后，根據(jù)每個(gè)特征點(diǎn)記錄的興趣值，在局部進(jìn)行極值抑制。具體過程為循環(huán)選取每個(gè)特征點(diǎn)，然后尋找其水平距離半徑r(本文所有r取經(jīng)驗(yàn)值3 m)范圍的所有特征點(diǎn)集S，如果當(dāng)前特征點(diǎn)的興趣值R為S中的最小值或者最大值，則保留當(dāng)前點(diǎn)，否則舍棄。經(jīng)過抑制局部極值，局部范圍內(nèi)只保留一個(gè)特征明顯的點(diǎn)。由于提取的特征點(diǎn)為離散3維點(diǎn)，不能用原始DEM搜索其鄰域點(diǎn)。為了快速搜索每個(gè)點(diǎn)周圍鄰域內(nèi)的點(diǎn)，本文采用K-D樹進(jìn)行搜索[9]。

綜上所述，本文基于多方向剖面線簡(jiǎn)化的特征提取過程可以歸納為圖4所示的流程。

如圖5所示為按照30°角度間隔對(duì)圖1(a)所示DEM進(jìn)行多方向特征提取的結(jié)果。圖5(a) 為多方向剖面線DP簡(jiǎn)化后直接疊加的結(jié)果，一共有32 739個(gè)點(diǎn)，進(jìn)行局部抑制極值后，保留6 458個(gè)點(diǎn)，結(jié)果如圖5(b)所示。

圖4 多方向剖面線簡(jiǎn)化特征提取流程圖

圖5 特征點(diǎn)提取結(jié)果示意

獲得DEM的特征點(diǎn)后，將這些點(diǎn)構(gòu)建三角網(wǎng)，然后再采用線性內(nèi)插的方式生成10 m大小格網(wǎng)的DEM，即為本文方法DEM綜合的結(jié)果。

3 實(shí)驗(yàn)結(jié)果及分析

3.1 評(píng)價(jià)指標(biāo)

獲得綜合后的規(guī)則格網(wǎng)DEM，分別采用高程中誤差、平均坡度對(duì)結(jié)果進(jìn)行定量評(píng)價(jià)；采用等高線套合對(duì)綜合以后的DEM進(jìn)行可視化定性分析[10]，指標(biāo)定義如下：

1) 高程中誤差：通過比較簡(jiǎn)化以后DEM上每個(gè)格網(wǎng)點(diǎn)高程與原始DEM對(duì)應(yīng)位置的高程計(jì)算而來。

2) 坡度均值差：指簡(jiǎn)化以后的DEM坡度平均值，與原始DEM的坡度平均值進(jìn)行比較的差值。差值越大說明與原始的DEM吻合越低，則DEM簡(jiǎn)化過程中，特征保留效果越差。

3) 等高線套合：分別對(duì)比簡(jiǎn)化后DEM生成等高線與原始等高線套合在原始DEM上的吻合程度，從整體上定性分析綜合前和綜合后的等高線。吻合程度越高，則表明綜合后特征保留越好。

3.2 結(jié)果與分析

利用DP算法進(jìn)行DEM特征點(diǎn)提取，傳統(tǒng)做法多采用0°和90°方向聯(lián)合提取，本文通過一定間隔角度旋轉(zhuǎn)，可以進(jìn)行多方向剖面線DP簡(jiǎn)化。本文對(duì)不同間隔角度進(jìn)行了實(shí)驗(yàn)，最后選取15°、30°、45°間隔的結(jié)果進(jìn)行定量和定性評(píng)價(jià)。

1) 定量評(píng)價(jià)結(jié)果。 高程中誤差和坡度均值差如表1所示。從表1結(jié)果可以看出高程中誤差隨著旋轉(zhuǎn)角度間隔的增大而增大。結(jié)果表明多方向的DP算法的精度比傳統(tǒng)的算法要高。表1中坡度均值差也表明多方向簡(jiǎn)化結(jié)果比傳統(tǒng)兩個(gè)方向的結(jié)果精度高。

表1 定量評(píng)價(jià)結(jié)果

2) 定性評(píng)價(jià)結(jié)果。 圖6是傳統(tǒng)方法和本文方法重構(gòu)DEM等高線套合結(jié)果；圖7為等高線套合局部放大結(jié)果(藍(lán)色為傳統(tǒng)方法重構(gòu)DEM的等高線，紅色為原始的DEM等高線，綠色、紫色、淡紫色分別為45°、30°、15°多方向的DP算法重構(gòu)DEM等高線)。從所得的等高線套合情況可以看出藍(lán)色等高線偏離紅色等高線最多，其余顏色的等高線與紅色等高線走勢(shì)更加吻合。從圖7的局部放大結(jié)果中也可以明顯看出藍(lán)色與紅色的吻合沒有其余顏色與紅色吻合得好。

圖6 等高線套合效果

圖7 等高線套合局部放大結(jié)果

實(shí)驗(yàn)結(jié)果表明本文所提出的多方向DP算法內(nèi)插的DEM比傳統(tǒng)的DP算法內(nèi)插的DEM的精度更好。傳統(tǒng)的DP算法內(nèi)插時(shí)只考慮了2個(gè)方向，而本文方法通過格網(wǎng)DEM的旋轉(zhuǎn)可以獲取多方向剖面線的整合結(jié)果。故而用多方向算法得到的重構(gòu)DEM比起傳統(tǒng)算法得到的重構(gòu)DEM更接近原始DEM。

4 結(jié) 論

針對(duì)DEM綜合的問題，本文提出一種利用多方向剖面線DP算法簡(jiǎn)化提取特征點(diǎn)的方法，從實(shí)驗(yàn)結(jié)果可看出，不管是從高程誤差統(tǒng)計(jì)還是坡度均值上，都可以看出本文方法得到的DEM與原始DEM地形形態(tài)更加吻合，能更好地保持地形形態(tài)特征。

本文在DEM綜合時(shí)，只考慮了特征點(diǎn)，進(jìn)一步將加入線特征，使得綜合后的DEM能夠更好地保持原始DEM特征。

[1] ZHOU Q, CHEN Y. Generalization of DEM for terrain analysis using a compound method [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2011,66(1),38-45.

[2] 董有福, 湯國(guó)安.利用地形信息強(qiáng)度進(jìn)行DEM地形簡(jiǎn)化研究[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2013,38(3):353-357.

[3] DOUGLAS D H, PEUCKER T K., Algorithms for the reduction of the number of points required to represent a digitized line or its caricature [J]. The Canadian Cartographer, 1973, 10(2):112-122 .

[4] 費(fèi)立凡,何津,馬晨燕,等.3維Douglas-Peucker算法及其在DEM自動(dòng)綜合中的應(yīng)用研究[J].測(cè)繪學(xué)報(bào), 2006(3): 278-284.

[5] 何津,費(fèi)立凡.再論三維Douglas-Peucker算法及其在DEM綜合中的應(yīng)用[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2008 (2):160-163.

[6] 朱雪堅(jiān),葉遠(yuǎn)智,湯國(guó)安.運(yùn)用三維Douglas-Peucker算法提取DEM地形特征[J].測(cè)繪通報(bào),2014(3):118-121

[7] MA T, CHEN Y, HUA Y, et al. DEM generalization with profile simplification in four directions[J]. Earth Science Informatics, 2017, 10(1): 29-39.

[8] Chen Z T, GUEVARA J A. Systematic selection of very important points (VIP) from digital terrain model for constructing triangular irregular networks[C]. In Auto-carto, 1987, 8, pp. 50-56.

[9] ARYA S, MOUNT D M. Approximate Nearest Neighbor Searching [C], Proc. 4th Annu. ACM-SIAM Sympos. on Discrete Algorithms(SODA’93), 1993, 271-280.

[10] 吳艷蘭,胡海,胡鵬,等.數(shù)字高程模型誤差及其評(píng)價(jià)的問題綜述[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2011, 36(5): 568-574.