李勝天，徐貴興

（江西省地質(zhì)局地理信息工程大隊，江西南昌 330001）

0 引言

數(shù)字高程模型（digital elevation model，DEM）的數(shù)據(jù)結(jié)構(gòu)簡單、容易管理、便于存貯，許多國家的DEM 數(shù)據(jù)都是以規(guī)則格網(wǎng)的數(shù)據(jù)矩陣形式提供的，因此，基于規(guī)則格網(wǎng)DEM 提取等高線、山脊線、山谷線一直是研究人員的研究方向[1]。山脊線、山谷線是地形骨架結(jié)構(gòu)的重要組成部分[2]，常應(yīng)用于地貌表達以及三維地形分析。建立格網(wǎng)DEM 的空間內(nèi)插方法可以分為分塊內(nèi)插法、整體內(nèi)插法、逐點內(nèi)插法3 類。其建立一般包括數(shù)據(jù)取樣、數(shù)據(jù)內(nèi)插和數(shù)據(jù)精度分析等步驟[3]。整體內(nèi)插法的原理簡單，適用于總體預(yù)測，但不適用于復(fù)雜地形，對于復(fù)雜地形，常使用分塊內(nèi)插法[4-5]。逐點內(nèi)插法是應(yīng)用最為廣泛的方法，包括距離倒數(shù)乘方法、克里金法，在樣本點分布均勻且數(shù)量較多時均可取得較好的結(jié)果。當(dāng)實際地形的坡度較大時，采用距離倒數(shù)乘方法較為適宜[6]。提取山脊線、山谷線常用的算法有3 類。一是幾何分析法，它利用地表幾何形態(tài)的實際規(guī)律，通過一定的算法（如曲率判別法、斷面極值法和骨架化法）提取地形特征點，其中斷面極值法與曲率判別法均是利用極值提取地形特征點。斷面極值法一般采用兩個方向的斷面，并不能提取出所有地形特征點[7-9]。骨架化法假設(shè)山脊線、山谷線兩側(cè)地形對稱，提取中心軸線，與實際不符，有很大程度的近似性，因此所得結(jié)果不盡如人意[10-11]。二是圖像處理法，它將DEM 數(shù)據(jù)作為灰度圖像處理，通過數(shù)學(xué)形態(tài)學(xué)的“序貫減薄”運算得到骨架線[12]，該方法可利用移動窗口法提取特征點，但是將特征點連成線非常困難。三是水文分析法，它利用地表水流的客觀規(guī)律對山脊線、山谷線進行提取，常用D8 流向法計算流向，在帶有匯聚型河流的多山地形能生成較好的結(jié)果。但在平坦地區(qū)不易確定流向，不適用于含泛濫平原和濕地的多變地形。對格網(wǎng)DEM 構(gòu)建常用的5 種算法，即距離倒數(shù)乘方法、自然鄰域法、趨勢面分析法、克里金法、樣條函數(shù)法進行研究，并分析它們在本溪關(guān)門山地區(qū)的適用性與準(zhǔn)確性，選出最適合的插值方法——距離倒數(shù)乘方法和自然鄰域法用于格網(wǎng)DEM 的構(gòu)建?；诟窬W(wǎng)DEM，利用水文分析法和曲率判別法兩種方法提取山脊線、山谷線，將提取結(jié)果進行對比分析，并研究和實現(xiàn)了兩種山脊線、山谷線提取算法的融合。

1 格網(wǎng)DEM 構(gòu)建及山脊線、山谷線提取方法與技術(shù)路線

本溪關(guān)門山地區(qū)，采用不同的空間插值方法提取山脊線、山谷線，通過GeoScene Pro 建模和二次開發(fā)的形式，分別完成格網(wǎng)DEM 的構(gòu)建、等高線的繪制、等高線曲率的計算、山脊線和山谷線的自動化提取。同時分析不同插值算法的適宜性，以及不同山脊線、山谷線提取算法的優(yōu)缺點。其技術(shù)路線步驟如下：

第一步：格網(wǎng)DEM 的構(gòu)建，利用5 種空間插值算法進行柵格插值，并通過交叉驗證選出最合適的方法用于本溪關(guān)門山地區(qū)格網(wǎng)DEM的構(gòu)建。

第二步：等高線的自動繪制，基于格網(wǎng)DEM 數(shù)據(jù)，生成等高線。

第三步：等高線曲率的計算，基于等高線數(shù)據(jù)，計算等高線曲率。

第四步：曲率判別法的實現(xiàn)，基于等高線曲率，完成山谷線及山脊線的提取。

第五步：水文分析法的實現(xiàn)，基于水文分析法，完成山谷線及山脊線的提取。

第六步：水文分析法與曲率判別法提取結(jié)果的對比分析以及兩種算法的融合。

2 空間插值算法實現(xiàn)與分析

空間插值算法構(gòu)建格網(wǎng)DEM有5種方法，包括距離倒數(shù)乘方法、克里金法、薄板樣條函數(shù)法、趨勢面分析、自然鄰域法。在此基礎(chǔ)上，對插值結(jié)果進行對比分析。在進行插值之前，首先將樣本點分為兩部分，其一為訓(xùn)練點，其二為測試點，樣本點共有37 660 個，訓(xùn)練點為1 000 個，測試點為36 660 個，其中訓(xùn)練點用于檢測插值質(zhì)量，測試點用于格網(wǎng)DEM 的構(gòu)建。在插值成功之后，將5種插值結(jié)果進行3次對比，即將樣本點隨機分配3次，進行獨立試驗，互不干擾。每次對比的測試點與訓(xùn)練點的分配比例相同，但分布各不相同，均是隨機的。根據(jù)測試點的5 種插值算法的插值計算結(jié)果，將訓(xùn)練點所在位置的值分別提取至訓(xùn)練點屬性表，共分別提取5 次。之后便可計算每種插值算法在每個訓(xùn)練點處的偏差。基于此屬性進行統(tǒng)計分析，分別計算均值、最大值、最小值、標(biāo)準(zhǔn)差、方差。統(tǒng)計表見表1。

從表1、圖1 可以看出，本實驗區(qū)（本溪關(guān)門山地區(qū)）應(yīng)用的5 種插值算法中，距離倒數(shù)乘方法、自然鄰域法的均方根誤差為4.533 0 和3.339 3，均值為-0.162 2 和-0.073 8，插值結(jié)果最好，最能反映真實地形?？死锝鸱ê蜆訔l函數(shù)法稍差，均方根誤差為5.603 2和5.748 9，均值分別為-0.293 4、-0.238 1，插值結(jié)果不如距離倒數(shù)乘方法和樣條函數(shù)法，但差距不大，也能反映真實地形。趨勢面分析法插值結(jié)果極差，標(biāo)準(zhǔn)差為123.728 0，偏差較大，只能反映整體趨勢，不能反映真實地形，不適用于DEM的建立。

圖1 構(gòu)建網(wǎng)格的5種空間插值算法插值結(jié)果Figure 1. Interpolation results of five spatial interpolation algorithms for grid construction

表1 三次插值結(jié)果對比Table 1. comparison of the results of cubic interpolation

3 山脊線、山谷線提取模型的建立

3.1 曲率計算模型

等高線曲率與GeoScene Pro 中的曲率有根本上的不同，等高線曲率是指地表面上通過某一點的等高面與地表的交線在地形表面水平方向上的凹凸性。

地形表面S的曲面方程為z=f(x,y)，點P為曲面S上任意一點，平面α為過點P的水平面，α與S的交線為c1過點P的等高線；a為等高線在c1點P的切向量，b為點P的坡向向量。b⊥a平面β、γ，分別過a、b垂直于水平面XOY并與曲面S相交于曲線c2、c3。曲面S：z=f(x,y)的各階偏導(dǎo)數(shù)用以下符號表示：

GeoScene Pro 中提供了平面曲率的計算工具，該曲率是垂直于坡度方向的曲率。從圖2中的數(shù)學(xué)關(guān)系來看，平面曲率就是曲線c2在點P的二階導(dǎo)數(shù)，它的計算公式為：

在GeoScene Pro 中，用最大平均值法計算坡度值。GeoScene Pro中的坡度可利用式（1）中的變量進行表達，其計算公式如下：

地形曲面S：z=f(x,y)的任意一條等高線方程可表示為f(x,y)=c（c為表示任意高程值的常數(shù)），由曲率計算公式可得等高線曲率計算公式為：

通過綜合公式（2）~（4），可以發(fā)現(xiàn)平面曲率、等值線曲率、坡度三者之間具有一定的關(guān)系，通過平面曲率與坡度便可計算等值線曲率柵格。

推導(dǎo)出的公式如下：

曲率判別法的實現(xiàn)思路是：首先獲取等值線曲率柵格，其次基于曲率柵格獲取地形的特征點，通過一定的分類閾值提取山脊點與山谷點。

3.2 水文分析模型

水文分析法的基本原理是：對于山脊線而言，它也是分水線，分水線所在的柵格區(qū)域在一定程度上不存在水流，即流量累積柵格的對應(yīng)值為零。通過地表徑流模擬計算可以得到流量累積柵格，并提取零累積值區(qū)域，即分水線所在區(qū)域，也就得到了山脊線[13]。若以某一相對高程面為對稱基準(zhǔn)面，求原始地形的對稱地形，則山谷線與山脊線位置也互換。因此，可以通過建立的反地形提取模型，求出對稱地形（反地形），那么，正地形山谷線就可以在反地形中利用提取山脊線的方法進行提取[14]。

（1）已填洼DEM。水文分析法中已填洼DEM 是指不存在小洼地，大多數(shù)小洼地是DEM 生成過程中帶來的數(shù)據(jù)錯誤所致，如果未進行填充，則生成的水系網(wǎng)絡(luò)可能會呈現(xiàn)不連續(xù)性。常用的除去洼地的方法是把洼地像元值加高至周圍最低像元值[15]。填洼操作不僅能去除小洼地，同時也可以去除錯誤的峰。峰是一種偽像元，其高程高于所預(yù)期的高程，是反地形中的小洼地。

（2）反地形DEM。傳統(tǒng)的水文分析法在提取山谷線過程中，反地形的計算是基于原始DEM 數(shù)據(jù)的，利用已填洼DEM 數(shù)據(jù)提取反地形，從而消除反地形中的峰和小洼地，進而消除山谷線提取結(jié)果的局部錯誤值。反地形的計算公式為：

式中：Abs為絕對值運算，Dem為已填洼DEM 數(shù)據(jù)；Dmax與Dmin分別為原始DEM 的最大高程值和最小高程值。模型中山脊線、山谷線提取閾值的確定和反地形的計算極為重要。對匯流累積量零值數(shù)據(jù)進行重分類時，屬性值越接近1 越可能是山脊線或山谷線的位置，分界閾值的確定尤為重要。根據(jù)原始DEM 建立的山體陰影（地貌暈渲圖），提取等值線，通過綜合3D 顯示的衛(wèi)星影像三方面的因素，輔助判讀，經(jīng)過反復(fù)驗證與試驗，得到山谷線提取的最佳閾值為0.510 2，山脊線提取的閾值為0.561 4。

3.3 融合算法模型

融合算法實現(xiàn)思路是對水文分析法和曲率判別法的一種集成，利用兩種算法進行加權(quán)評價，分值越高，代表兩種算法均認(rèn)可提取結(jié)果的準(zhǔn)確，分值越低，則與之相反。具體的實現(xiàn)過程是：首先利用創(chuàng)建的等值線曲率計算模型或腳本，獲取等值線曲率；同時，截取水文分析模型中的上半部分，獲取蓄積柵格數(shù)據(jù)（流量）和正反地形；然后基于等值線曲率，利用條件函數(shù)與閾值獲取初始山脊線曲率、初始山谷線曲率?；谒姆治霁@取初始山脊線流量、山谷線流量；再者，基于上一步的初始山脊線、山谷線，使用分位數(shù)法將以上初始數(shù)據(jù)各自分為8個等級，等級越高，代表正確的可能性越大。最后，將拉伸過的初始山脊線曲率與初始山脊線流量數(shù)據(jù)進行相加，將結(jié)果拉伸到0～1，越接近于1，結(jié)果是山脊線的可能性越大。山谷線與之相同。

4 山脊線和山谷線提取結(jié)果的對比分析

4.1 水文分析法與曲率判別法的對比分析

采用兩種山脊線、山谷線的提取算法，一是水文分析法，二是曲率判別法。運用兩種算法進行了建模實現(xiàn)與開發(fā)實現(xiàn)，并對兩種算法的提取山脊線、山谷線的結(jié)果分別進行對比分析，兩種算法各有優(yōu)缺點。首先建立了對比分析模型，將水文分析法與曲率判別法提取山脊線、山谷線的結(jié)果進行分析，然后在本模型中將該對比分析功能進行了實現(xiàn)。對比分析模型建立的原理主要是地圖代數(shù)運算與柵格局域運算，得到的分析結(jié)果柵格具有4 個值，即為4 個類型，其中“0”值代表兩種算法均為提取的區(qū)域，而并不是山脊線、山谷線區(qū)域；“1”值代表曲率分析法提取的單獨區(qū)域，此區(qū)域水文分析法未提取出山脊線、山谷線；“2”值代表水文分析法提取的單獨區(qū)域，此區(qū)域曲率分析法未提取出山脊線、山谷線。

利用水文分析法與曲率判別法分別提取山脊線的結(jié)果（柵格）對比如圖2 所示，提取山脊線的像元數(shù)如表2所示。利用水文分析法與曲率判別法分別提取山谷線的結(jié)果（柵格）對比如圖3 所示，提取山谷線的像元數(shù)如表3所示。

表2 兩種算法提取山脊線的像元數(shù)對比Table 2. Comparison of pixel elements of ridge lines extracted by two algorithms

表3 兩種算法提取山谷線的像元數(shù)對比Table 3. Comparison of pixel elements of valley lines extracted by two algorithms

圖2 兩種算法提取山脊線的結(jié)果對比Figure 2. Comparison of ridge lines extracted by two algorithms

圖3 兩種算法提取山谷線的結(jié)果對比Figure 3. Comparison of valley lines extracted by two algorithms

4.1.1 兩種算法提取山脊線結(jié)果的對比分析

綜合圖2 和表2 的對比結(jié)果，可以看出兩種算法提取山脊線均取得了很好的效果，在整體上都是正確的。相比而言，曲率判別法提取的范圍較為全面，但連續(xù)性較差，水文分析法提取山脊線的連續(xù)性較好，但局部區(qū)域未能提取。山脊線的提取結(jié)果受限于多方面的因素，首先是DEM 的質(zhì)量，其次是填洼、流向、流量、坡度、曲率等中間數(shù)據(jù)生成算法的影響。以流向計算為例，在GeoScene Pro中采用的是D8流向法，即八方向法，該方法并非最優(yōu)算法，但相對而言卻是效率高、精度高的算法。最后是人為干涉因素，包括各類閾值的設(shè)定，需要人為輔助判斷。

4.1.2 兩種算法提取山谷線結(jié)果的對比分析

綜合圖3 和表3 的對比結(jié)果，可以看出水文分析法與曲率判別法在山谷線的提取方面各有優(yōu)缺點，首先，二者提取的結(jié)果都是正確的，但都有不足之處。水文分析法提取山谷線的連續(xù)性效果較好，在主山谷線處提取結(jié)果突出，正確率較高，且少有斷裂現(xiàn)象發(fā)生，但是在分支處的提取結(jié)果較差，部分局部山谷線未得到提取。曲率判別法提取山谷線的整體性較好，局部的山谷線也得到了提取，但連續(xù)性較差，有斷裂情況發(fā)生，局部存在錯誤點，需要進行眾數(shù)濾波算法處理錯誤值。山谷線的提取結(jié)果同樣受限于多方面的因素，與山脊線的影響因素相同，首先是DEM的質(zhì)量，其次是流向、流量提取算法，最后是各類閾值的設(shè)定。

4.2 融合算法與水文分析法、曲率判別法的對比分析

通過表2 和表3 中的數(shù)據(jù)來看，水文分析法和曲率判別法單獨提取區(qū)域占總柵格區(qū)域的2.50%～3.56%，將該部分?jǐn)?shù)據(jù)通過圖3 和圖4 進行分析對比，可以看出：水文分析法提取的單獨部分主要位于主體特征線，曲率判別法單獨提取的部分大部分是水文分析法未提取的局部山脊線，因此均是正確的。總的來說，二者各有長短，優(yōu)勢互補。水文分析法的主體提取結(jié)果較好，連續(xù)性強，斷裂線少，但局部提取結(jié)果差；而曲率判別法局部提取結(jié)果較好，但連續(xù)性較差。因此將二者算法進行融合，得到的結(jié)果相對而言是更為全面而準(zhǔn)確的。研究了兩種算法的融合，利用水文分析提取的山脊線、山谷線結(jié)果輔助曲率判別法的閾值設(shè)定，并將兩種算法的提取結(jié)果進行疊加，利用眾數(shù)濾波剔除局部錯誤值。由兩種方法融合后提取結(jié)果可以明顯看出，融合算法的提取結(jié)果較水文分析法和曲率判別法有所提高，無論是整體還是局部的山脊線與山谷線都得到了提取。3 種結(jié)果的對比如圖4所示。

圖4 曲率判別法、水文分析法、融合算法分別提取的山谷線、山脊線對比圖Figure 4. Comparison of valley lines and ridgelines extracted by curvature discrimination method, hydrologic analysis method and fusion algorithm

5 結(jié)論與展望

5.1 結(jié)論

基于GeoScene Pro 建模和腳本，實現(xiàn)了DEM 構(gòu)建、等高線繪制、等高線曲率計算、山脊線與山谷線的提取等功能。主要研究結(jié)果：

（1）完成了格網(wǎng)DEM 的構(gòu)建，包括5 種空間插值算法：距離倒數(shù)乘方法、薄板樣條函數(shù)法、趨勢面分析法、自然鄰域法、克里金法。5 種算法中，自然鄰域法和距離倒數(shù)乘方法對本溪關(guān)門山地區(qū)的適宜性最佳。

（2）實現(xiàn)了山脊線、山谷線的自動提取，完成了水文分析法與曲率判別法兩種算法對山脊線與山谷線的自動提取。提取結(jié)果顯示水文分析法連續(xù)性較好，但存在局部特征線未提取，曲率判別法整體性較好，提取結(jié)果較為全面，但存在特征線斷裂情況。

（3）研究了水文分析法與曲率判別法的融合，在兩種算法的基礎(chǔ)上實現(xiàn)了融合算法對山脊線與山谷線的自動提取。并對水文分析法、曲率判別法以及二者融合算法提取結(jié)果進行了對比分析，建立了對比分析模型，將分析結(jié)果量化顯示。研究發(fā)現(xiàn)融合算法兼具二者之長，提取的山脊線與山谷線不僅連續(xù)性好，而且更加完整、全面。

（4）提出了基于坡度與平面曲率計算柵格等高線曲率的方法，并實現(xiàn)了柵格形式等高線曲率的計算，同時也實現(xiàn)了矢量形式等高線曲率的計算。通過對比發(fā)現(xiàn)，柵格形式的等高線曲率計算速度較快，且精度與準(zhǔn)確性更高。

5.2 展望

對柵格空間插值方法、等高線的曲率計算方法、山脊線與山谷線的提取方法做了初步探索，存在以下問題有待進一步研究：在山脊線、山谷線的提取方面，仍需驗證更多的提取算法，并進行對比分析；水文分析法與等高線曲率判別法的融合策略有待進一步的加深。在今后的學(xué)習(xí)和工作中，需要不斷探索，吸取更多經(jīng)驗，綜合多領(lǐng)域的技術(shù)，繼續(xù)加深研究，提高算法的效率與準(zhǔn)確性，進一步優(yōu)化“格網(wǎng)DEM 的構(gòu)建及山脊線、山谷線提取模型”。