李廷豪，仲濟(jì)強(qiáng)

(南京海河測繪科技有限公司，江蘇 南京 210008)

0 引言

數(shù)字高程模型(DEM)的出現(xiàn)使得地形的表達(dá)更加的多元化，人們可以更加快速、簡便、精確地了解地形信息。數(shù)字高程模型是地面高程信息的數(shù)字表示形式，對于地學(xué)領(lǐng)域的研究具有廣泛的應(yīng)用價(jià)值[1]。與傳統(tǒng)的地形圖相比，DEM可以用于工程規(guī)劃建設(shè)、景觀分析、日照分析等。山脊線和山谷線代表了地形最為重要的結(jié)構(gòu)信息，在水文分析、地形重建等方面具有重要的意義[2]。隨著地理信息系統(tǒng)的快速發(fā)展，提取地形線的方法也越來越多，提取的準(zhǔn)確性和便捷性等方面也得到了一定的提升。主要方法有：基于圖像處理技術(shù)的方法[3]、基于地形表面幾何分析的方法[4-5]、基于地表流水物理模擬的方法[6]以及基于地形表面幾何分析與地表流水物理模擬相結(jié)合的方法[7]。不同的提取方法具有不同的優(yōu)勢以及缺陷，本文基于地表流水物理模擬方法結(jié)合數(shù)理統(tǒng)計(jì)方法來提取山脊線與山谷線，以期得到更加精準(zhǔn)的提取效果。

1 基于地表流水物理模擬方法

水從高處向低處流，在山谷處流量大，形成匯水，在山脊處水量小，形成分水等，因此山谷線與山脊線也被稱為合水線和分水線。通過分析地表水流的這種自然規(guī)律，結(jié)合合水線和分水線的匯水量的累計(jì)特征，即合水線(山谷線)的匯水量單調(diào)遞增，而分水線(山脊線)上的匯水量是不累積的，實(shí)現(xiàn)對地形線進(jìn)行提取，稱為基于地表流水物理模擬的地形特征提取方法。

其基本原理是首先從每個(gè)點(diǎn)8個(gè)方向找出坡度下降最大的方向，然后根據(jù)高程從大到小依次計(jì)算每個(gè)點(diǎn)的匯水量，其中區(qū)域匯水量最大的點(diǎn)為匯水線的終結(jié)點(diǎn)，從終結(jié)點(diǎn)出發(fā)，找出與該匯水線的終結(jié)點(diǎn)相鄰近的點(diǎn)，并且該點(diǎn)的匯水量遠(yuǎn)遠(yuǎn)小于終結(jié)點(diǎn)的匯水量，同時(shí)大于其他點(diǎn)的匯水量，確定這個(gè)點(diǎn)為新的匯水線上的終結(jié)點(diǎn)，直到判定點(diǎn)的所有鄰近點(diǎn)的匯水量差異較小，即可得到矢量化的匯水線數(shù)據(jù)。通過計(jì)算流向匯水線的點(diǎn)構(gòu)成的區(qū)域的邊界線得到匯水區(qū)域的邊界線。大量研究表明，這種方法提取的地形線能保持較好的連貫性，但提取結(jié)果的質(zhì)量與鄰域窗口的大小有較大的關(guān)系?；诖?，本文以該方法為基礎(chǔ)，提出采用統(tǒng)計(jì)學(xué)中的均值變點(diǎn)分析法選取最佳鄰域窗口，以提高提取的準(zhǔn)確性。

2 基于均值變點(diǎn)分析法的鄰域窗口確定

正負(fù)地形是兩種相反的地形，分別高于或者低于鄰區(qū)?；诘乇砹魉锢砟M的地形特征提取方法需要對原始DEM進(jìn)行正負(fù)地形的提取。在提取正負(fù)地形時(shí)首先需要進(jìn)行鄰域分析，計(jì)算鄰域均值DEM，鄰域窗口的大小直接影響均值DEM的計(jì)算，從而也就影響地形線的提取結(jié)果[8]。鄰域分析的基本原理是以待計(jì)算柵格為中心，向其周圍擴(kuò)展一定范圍，基于這些擴(kuò)展柵格數(shù)據(jù)進(jìn)行函數(shù)運(yùn)算，從而得到此柵格的值。這里采用鄰域分析的焦點(diǎn)統(tǒng)計(jì)法，通過訪問柵格中的每個(gè)像元，根據(jù)識別出的鄰域范圍來計(jì)算高程平均值，得到均值DEM。

由于鄰域窗口的大小直接關(guān)系到地形線的提取效果，同時(shí)顧及山脊線和山谷線的特點(diǎn)，本文提出采用均值變點(diǎn)分析法來選取合適的窗口。均值變點(diǎn)分析法是一種處理非線性數(shù)據(jù)的數(shù)理統(tǒng)計(jì)方法，特別對只有一個(gè)變點(diǎn)的情況更為有效。對于一非線性的系統(tǒng)輸出序列，若在某未知時(shí)刻或空間位置突然發(fā)生變化，這一點(diǎn)就稱為變點(diǎn)[9]。主要計(jì)算過程為：

設(shè)有樣本序列M，令i=2，3，…，N，第i個(gè)樣本將序列分成兩段：X1，X2，…，Xi-1與Xi，Xi+1，...，XN，分別統(tǒng)計(jì)兩段樣本的算術(shù)平均值：

(1)

(2)

計(jì)算兩段樣本的離差平方和之和為：

(3)

計(jì)算樣本的算術(shù)平均值與樣本的離差平方和：

(4)

(5)

計(jì)算期望值：

(6)

由于變點(diǎn)可以使S-Sn的差距變大，因此可以根據(jù)計(jì)算S-Sn值變化來判斷變點(diǎn)存在的位置。

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

3.1 正負(fù)地形提取

結(jié)合地形起伏程度等相關(guān)因素，本文選取安徽省黃山市黃山區(qū)進(jìn)行研究，黃山區(qū)位于東經(jīng)117°50′與118°21′，北緯30°00′與30°32′之間。黃山區(qū)整體地勢起伏較大，其中最高峰海拔達(dá)到1 864 m，最低海拔為85 m。本文將地理空間數(shù)據(jù)云上下載的SRTM 90 m分辨率的黃山區(qū)數(shù)字高程數(shù)據(jù)經(jīng)過裁剪作為實(shí)驗(yàn)DEM數(shù)據(jù)。首先對原始DEM數(shù)據(jù)進(jìn)行正負(fù)地形的提取，以獲得兩種相反的地形來對兩種相反的地形線分別進(jìn)行提取分析。在進(jìn)行正負(fù)地形提取時(shí)首先要進(jìn)行鄰域分析來獲取均值DEM數(shù)據(jù)，其中鄰域窗口的選擇至關(guān)重要，鄰域窗口選取過大或者過小都影響著山脊線和山谷線的提取，而均值變點(diǎn)分析法的使用則較好地解決了窗口大小的選取問題。不同鄰域窗口下均值變點(diǎn)分析法統(tǒng)計(jì)結(jié)果見表1。

表1 不同鄰域窗口下均值變點(diǎn)分析法統(tǒng)計(jì)結(jié)果Tab.1 Statistical results of mean change point analysis under different neighborhood windows

如表1所示，為確定最佳鄰域窗口，將鄰域窗口設(shè)為2×2一直到30×30依次進(jìn)行運(yùn)算，得到不同窗口下的均值DEM。再通過前文公式分別計(jì)算樣本離差平方和與兩段樣本的離差平方和。通過計(jì)算可以得到S=20.345 6，進(jìn)一步計(jì)算在不同鄰域窗口大小下兩者的差值。

如圖1所示，將格網(wǎng)大小與均值DEM進(jìn)行擬合發(fā)現(xiàn)隨著格網(wǎng)面積的變化，均值DEM也單調(diào)變化，剛開始均值DEM上升速度較慢，而到了一定點(diǎn)，增長速率急劇上升，這個(gè)點(diǎn)稱為變點(diǎn)。通過分析S-Sn的變化趨勢來確定變點(diǎn)的具體位置，發(fā)現(xiàn)當(dāng)鄰域窗口為22×22時(shí)，S-Sn的值達(dá)到最大，因此在進(jìn)行鄰域分析的時(shí)候，將22×22窗口作為最佳鄰域窗口來得到最合適的均值DEM，以達(dá)到更好的提取效果，提高結(jié)果的準(zhǔn)確性與精密性。

圖1 格網(wǎng)大小與均值DEM及S-Sn的關(guān)系Fig.1 Relationship between grid size， mean DEM and S-Sn

通過以上方法計(jì)算均值DEM，將原始DEM與均值DEM作差可以得到研究區(qū)域的概略DEM數(shù)據(jù)，對概略DEM數(shù)據(jù)進(jìn)行兩次重分類可以得到研究區(qū)域的正負(fù)地形結(jié)果。為了驗(yàn)證采用均值變點(diǎn)分析法所得到的正負(fù)地形具有更好的效果，本文選取最佳鄰域窗口大小的同時(shí)選取其他兩組不同的窗口大小進(jìn)行實(shí)驗(yàn)，獲得不同的正負(fù)地形結(jié)果進(jìn)行對比分析，圖2(a)、(b)、(c)分別是鄰域窗口為15×15、22×22、28×28時(shí)得到的正地形結(jié)果圖，圖2(d)、(e)、(f)為對應(yīng)的負(fù)地形結(jié)果圖。從圖中可以發(fā)現(xiàn)，不同鄰域窗口下的正負(fù)地形有所差異，當(dāng)鄰域窗口較小時(shí)，破碎點(diǎn)比較多，當(dāng)鄰域窗口比較大時(shí)，提取的正負(fù)地形又比較過度化，而當(dāng)鄰域窗口適中時(shí)，提取的正負(fù)地形比較完整。因此本文選取22×22的鄰域窗口具有一定的可行性。

圖2 不同鄰域窗口下的正負(fù)地形Fig.2 Positive and negative terrain under different neighborhood windows

3.2 山脊線、山谷線提取

與山脊線不同，在進(jìn)行山谷線的提取時(shí)，首先要提取原始DEM的反地形DEM數(shù)據(jù)，反地形與正負(fù)地形概念不同，它不是客觀存在的，只是基于一個(gè)基準(zhǔn)面將原始地形來進(jìn)行翻轉(zhuǎn)而得到的地形。由如圖3所示，基于ArcMap對原始DEM以及反地形DEM進(jìn)行洼地填充、水流方向提取、匯流累積量，得到原始DEM匯流累計(jì)結(jié)果與反地形DEM匯流累計(jì)結(jié)果，將兩個(gè)結(jié)果分別與正負(fù)地形進(jìn)行相交可以得到研究區(qū)域的山脊線與山谷線柵格結(jié)果，對柵格結(jié)果進(jìn)行自動(dòng)矢量化得到的矢量結(jié)果圖。

圖3 山脊線、山谷線矢量化結(jié)果Fig.3 Vectorization results of ridge line and valley line

4 結(jié)束語

鄰域窗口的選擇對山脊線與山谷線的提取具有重要的影響，因此本文在基于地表流水物理模擬的基礎(chǔ)之上，通過數(shù)理統(tǒng)計(jì)方法選取最佳鄰域窗口來進(jìn)行山脊線與山谷線的提取。由實(shí)驗(yàn)分析可以得到，不同的鄰域窗口所得到的正負(fù)地形有所差異，通過對均值變點(diǎn)分析法的原理介紹以及應(yīng)用，得到最佳鄰域窗口，該方法提高了地形線提取的準(zhǔn)確性，同時(shí)這種方法提取的自動(dòng)化有待進(jìn)一步研究。