Python輔助下的坡度提取不確定性分析

吳 餒，李朝奎，方 軍，褚 楠，吳柏燕

(1. 湖南科技大學地理空間信息技術國家地方聯合工程實驗室，湖南 湘潭 411201； 2. 湖南科技大學地理空間信息湖南省工程實驗室，湖南 湘潭 411201)

滑坡和泥石流是主要的地質災害形式，其基本形成條件離不開地形地貌[1-2]，而坡度是地貌學中描述地貌形態的兩個重要因子之一[3]，對地表過程有強烈影響，具有重要的地理學意義[4]。

在區域尺度的研究中，坡度通常基于DEM提取[5]，DEM只是真實地表形態的近似表達[6]，而不同分辨率DEM數據提取的坡度精度存在差異，地形起伏等因素導致差異增大。因此，了解坡度提取誤差的成因、大小和空間分布規律有利于改進坡度提取方法，提高坡度提取的精度。

國內許多學者在基于DEM提取坡度的不確定性分析中作了大量的研究。趙牡丹等選取陜北黃土丘陵溝壑區、黃土梁峁丘陵區、黃土低丘區及黃土破碎塬區4種地貌類型區作為試驗樣區，采用高精度1∶1萬比例尺5 m分辨率的DEM為基準數據，應用比較分析、回歸分析和相關分析等方法研究黃土高原不同分辨率DEM提取地面平均坡度精度的量化估算方法[7]；湯國安等選擇陜北黃土高原6個典型地貌類型區為試驗樣區，采用野外實測及高精度的1∶1萬比例尺DEM為基準數據，研究柵格分辨率及地形粗糙度對DEM所提取地面平均坡度精度的影響[8]；劉敏等以陜西省671幅1∶5萬比例尺，25 m分辨率的DEM數據為信息源，分析DEM分辨率對提取坡度信息的不確定性影響[9]。

上述研究多從理論層面研究坡度提取的不確定性，但很少考慮整個計算過程的方法誤差，以及可能存在的人為誤差。本文針對現有坡度提取方法，運用以Python為開發語言的地理建模技術，從實踐層面研究DEM分辨率對坡度提取的影響，整合分析過程，提高分析效率，減少人為誤差。研究結論對完善現有坡度提取理論具有顯著的理論意義。

1 坡度計算及不確定性分析方法

1.1 基于DEM的坡度計算

地表面上任一點的坡度是指過該點的切平面與水平地面的夾角[10]，表示地表面在該點的傾斜程度。因此，坡度是一個既有大小又有方向的量，即矢量。在數值上等于過該點的地表微分單元的法矢量n與地表曲面函數Z軸的夾角,如圖1所示。用公式表示為

圖1 地表單元坡度示意圖

(1)

式中，n為地面一點的微分單元的法矢量；z為Z軸分量；|n|與|z|分別為各自方向的模。地面上每一點都有坡度，它是一個點上的微分概念，是地表曲面函數z=f(x,y)在東西、南北方向上的高程變化率的函數。因此，進行坡度提取常采用簡化差分公式計算

(2)

式中，fx為X方向的高程變化率；fy為Y方向的高程變化率。由式(2)可知，fx與fy的求解是獲取坡度值的關鍵。由于DEM數據是離散的，曲面函數f(x,y)未知，故坡度的提取通常是在3×3的柵格分析窗口中(如圖2所示)，采用幾何平面來擬合或差分計算的方法進行。

圖2 三階差分

基于DEM格網計算fx、fy可以有以下3種不同的計算方法

(3)

(4)

(5)

式中，d表示DEM的分辨率，即柵格格網的距離大小；fx與fy的含義與式(2)同。從式(3)—式(5)中可以看出，fx與fy的值主要受格網大小d的影響，因此坡度的提取精度也將明顯受到格網分辨率的限制。為了提高坡度計算精度，式(5)常被用于求解fx與fy。

1.2 坡度計算的不確定性分析方法

Python作為一種通用的開源編程語言，易于學習和使用，在ArcGIS地理處理中用作腳本語言。而ArcGIS中的Arcpy站點包給Python處理地理數據提供了基礎[11-12]。本文基于Arcpy站點包編寫Python腳本程序對DEM數據進行批量重采樣、批量坡度提取及批量掩膜提取。利用numpy包對數據進行整理，scipy包對數據進行科學計算并分析，matplotlib包將數據分析結果制成圖表。

在Model Builder中構建模型工具。Model Builder也叫模型構建器，是ArcGIS提供的一種通過構造工作流與處理腳本的圖形化建模工具。利用模型構建器將需要用到的工具連接，前一個工具的輸出作為后一個工具的輸入，形成工作流，從而不用頻繁調用工具。本文使用Model Builder構建柵格河網模型計算溝壑密度。為了便于操作，將腳本工具和模型工具都集成到ArcGIS的工具箱中統一管理與調用。

2 試驗及結果分析

2.1 試驗數據及試驗流程

采用平均坡度來表征區域坡度的一般水平[13]。以天津薊縣為試驗區，地貌類型如圖3所示。根據地貌類型將其劃分為侵蝕中低山、侵蝕丘陵、侵蝕堆積山間盆地、平原洼地、洪積沖積傾斜平原、沖積平原、沖積洪積傾斜平原7個地貌類別區域。

以30 m空間分辨率的DEM數據為基礎數據，對上述7類地貌類型區坡度提取的不確定性進行分析，具體試驗流程如下：

(1) 在批量重采樣Python腳本代碼中，調用arcpy.Resample_management()函數，使用循環將30 m分辨率DEM重采樣成40、50、60、70、80、90、100、110、120 m共10組不同分辨率的DEM數據。采樣結果如圖4所示。

注：1 miles=1.609 344 km圖3 天津薊縣地貌類型

圖4 不同分辨率DEM

(2) 在批量坡度提取Python腳本中，調用arcpy.Slope()函數，使用循環功能對上述重采樣的DEM數據分別提取坡度值，如圖5所示。調用掩膜提取函數arcpy.ExtractByMask()提取出每種地貌類型區域的平均坡度，結果見表1。

圖5 不同分辨率DEM提取的坡度

(°)

(3) 運用分辨率與平均坡度分析模型進行分析計算，得到以DEM分辨率為橫坐標、7種地貌類型區域的平均坡度值為縱坐標的散點圖，并對散點圖進行擬合，通過回歸算法求解得到回歸方程的系數及常數項。

(4) 建立不同地貌類型區域的4種地形因子，即地形曲率、地表起伏度、地面粗糙度及溝壑密度，以此作為坡度提取的地形因子并構建對應的模型工具。分別擬合上述回歸方程中求得的常數項，建立平均坡度與DEM分辨率、區域地貌特征的關系。以溝壑密度模型提取為例，其提取模型及過程[14-15]如圖6所示，其他因子提取與此相同。

圖6 溝壑密度提取模型

2.2 試驗結果分析

以DEM的分辨率為自變量，以不同地貌區域內的平均坡度值為因變量，定量分析平均坡度的變異特征。采用分辨率與平均坡度模型計算功能，其計算結果如圖7所示。

圖7 不同地貌單元的分辨率與平均坡度回歸關系

由圖7可知，平均坡度值與分辨率呈現很強的線性關系，且是正相關關系。因此以線性方程來描述這種關系，方程式如下

y=ax+b

(6)

式中，x為DEM分辨率；y為平均坡度；a為回歸方程系數；b為常數項。則7種地貌類型區域內的平均坡度與分辨率的回歸關系計算結果見表2。

表2 不同地貌單元內平均坡度與分辨率的回歸關系

從表2可見，除了侵蝕中低山地貌類型，其他地貌類型的曲率基本一致，故取其他6種地貌類型曲率的平均值作為研究區的平均曲率，即a=-0.012 4，分別對4種地形因子模型進行計算，其計算結果如圖8所示。

由圖8可見，不同地貌單元內的地形因子溝壑密度值與相應的回歸方程的常數項b存在顯著的二次函數關系，其二次回歸方程如下

b=0.004 7x2-0.095 5x+1.247 1

(7)

式中，x為不同地貌類型所對應的溝壑密度值；b為不同地貌類型區域內平均坡度與分辨率線性回歸方程的常數項。由此得到試驗區域內平均坡度的變異特征如下

F(g，x)=-0.012 4g+(0.004 7x2-0.095 5x+1.247 1)

(8)

式中，g為DEM空間分辨率；F(g，x)為對應的平均坡度值；x為溝壑密度值。

圖8 不同地形因子與常數項b的關系

假定以分辨率g=30 m的DEM提取的平均坡度值為參考M，則M可以表示為

M(x)=-0.372+0.004 7x2-0.095 5x+1.247 1

(9)

式中，x的含義與式(7)、式(8)同。則以其他分辨率提取的平均坡度的誤差E估算模型如下

E=M(x)-F(g,x)= -0.372+0.012 4g

(10)

式中，g為DEM空間分辨率；E為試驗區域內平均坡度提取誤差。可見坡度的提取平均誤差只與DEM的空間分辨率有關。

3 結 論

針對現有基于DEM的坡度提取方法，通過地理建模技術，從實踐層面研究DEM分辨率對坡度提取的影響，整合分析過程，提高分析效率，減少人為誤差。以天津薊縣為研究區，采用ArcGIS平臺和Python腳本語言功能分析了基于DEM的坡度提取方法的不確定性。研究結果表明：

(1) 研究區在不同地貌類型下提取的平均坡度與DEM空間分辨率之間呈現很強的線性關系，平均坡度隨著DEM分辨率的減小而減小，但衰減速率基本不變。

(2) 不同地貌單元平均坡度與DEM分辨率之間回歸方程的常數項與溝壑密度呈顯著的二次函數特征。

(3) 坡度提取的精度與DEM的分辨率呈正相關，即坡度的提取精度隨DEM分辨率的升高而升高。

(4) 基于Python的地理建模技術能有效地整合坡度提取精度分析內容，極大地提高坡度提取的不確定性分析效率。

研究結論進一步驗證了當前基于DEM的坡度提取理論方法的實用性，是對現有理論方法的有效補充，具有一定的理論價值。

參考文獻：

[1] 劉成，徐剛，楊娟，等.激發坡面泥石流發生的坡度因素探討——以重慶北碚地區為例[J].中國地質災害與防治學報，2006，17(4)：165-168.

[2] 郭芳芳，楊農，孟暉，等.地形起伏度和坡度分析在區域滑坡災害評價中的應用[J].中國地質，2008，35(1)：131-143.

[3] 楊勤科，賈大韋，李銳，等.基于DEM的坡度研究——現狀與展望[J].水土保持通報 ,2007，27(1)：146-150.

[4] 李天文，劉學軍，湯國安.地形復雜度對坡度坡向的影響[J].山地學報，2004，22(3)：272-277.

[5] 馬思煜，王春梅，土祥，等.坡度隨DEM分辨率變化的地統計學解釋——以黃土丘陵溝壑區為例[J].地下水，2015，37(5)：278-280.

[6] 陳楠.基于信息損失量的坡度精度分析[J].地球信息科學，2014，16(6)：852-858.

[7] 趙牡丹，陳正江，晉銳.DEM提取地面平均坡度誤差的量化模擬[J].水土保持通報，2002，22(3)：55-57.

[8] 湯國安，趙牡丹，李天文，等.DEM提取黃土高原地面坡度的不確定性[J].地理學報，2001，58(6)：824-830.

[9] 劉敏，湯國安，王春，等.DEM提取坡度信息的不確定性分析[J].地球信息科學，2007，19(2)：65-69.

[10] 湯國安，李發源，劉學軍.數字高程模型教程[M].3版.北京:科學出版社，2016.

[11] 劉雄，李宇，陳剛.Arcpy支持下的線要素幾何自動接邊與屬性檢查[J].測繪通報,2016(8)：114-117.

[12] 余詠勝,彭艷麗,尹言軍,等.基于Arcpy的影像地圖自動處理技術研究[J].測繪通報,2015(3)：82-85.

[13] 胡卓瑋，李洋，王志恒.基于DEM的四川省低山丘陵區坡度提取不確定性分析[J].山地學報，2012，30(5)：636-640.

[14] 師長興，周園園，范小黎，等.利用DEM進行黃河中游河網提取及河網密度空間差異分析[J].測繪通報,2012(10)：24-27.

[15] 劉學軍，盧華興,卞璐，等.基于DEM的河網提取算法的比較[J].水利學報，2006，37(9):1134-1141.

歡迎訂閱《測繪通報》

《測繪通報》創刊于1955年，是由國家測繪地理信息局主管、中國地圖出版社(測繪出版社)主辦的反映我國測繪地理信息科技發展現狀和指導全國測繪地理信息生產業務的綜合性、技術性刊物。現為《中文核心期刊要目總覽》收錄期刊、中國科技核心期刊(中國科技論文統計源期刊)、中國科學引文數據庫(CSCD)來源期刊、中國學術期刊網絡出版總庫(CAJD)全文收錄期刊、中國學術期刊綜合評價數據庫(CAJCED)核心統計源期刊、《中國學術期刊文摘》收錄期刊。曾被評為第二、三屆國家期刊獎百種重點期刊，多次榮獲全國優秀測繪地理信息期刊獎。

《測繪通報》的主要內容包括大地測量、全球導航衛星系統(GNSS)、攝影測量、遙感(RS)、地圖制圖、地理信息系統(GIS)、工程測量、礦山測量、地籍測繪、海洋測繪、測繪儀器、信息傳輸、圖形圖像處理等方面的新成果和新技術；行業管理、科研、教學、生產的先進經驗；計算機、通信等相關理論技術在測繪地理信息領域里的應用及測繪地理信息科技在國家經濟建設各個方面的應用；國內外測繪地理信息學術動態及有關測繪地理信息科技信息等。

《測繪通報》的主要欄目有學術研究、技術交流、測繪地理信息論壇、經驗介紹、測量員之窗、國外測繪地理信息、行業管理、行業研究、行業調查、知識窗、新書介紹、企業之窗、測繪地理信息教學等。

《測繪通報》為月刊，2018年定價：35.00元，郵發代號：2-223。

編輯部地址：北京市西城區三里河路50號，郵編：100045，訂閱電話：010-68531192(金老師)，68531317(傳真)。

網址：http:∥tb.sinomaps.com