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

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

(1. 湖南科技大學(xué)地理空間信息技術(shù)國(guó)家地方聯(lián)合工程實(shí)驗(yàn)室，湖南 湘潭 411201； 2. 湖南科技大學(xué)地理空間信息湖南省工程實(shí)驗(yàn)室，湖南 湘潭 411201)

滑坡和泥石流是主要的地質(zhì)災(zāi)害形式，其基本形成條件離不開(kāi)地形地貌[1-2]，而坡度是地貌學(xué)中描述地貌形態(tài)的兩個(gè)重要因子之一[3]，對(duì)地表過(guò)程有強(qiáng)烈影響，具有重要的地理學(xué)意義[4]。

在區(qū)域尺度的研究中，坡度通?；贒EM提取[5]，DEM只是真實(shí)地表形態(tài)的近似表達(dá)[6]，而不同分辨率DEM數(shù)據(jù)提取的坡度精度存在差異，地形起伏等因素導(dǎo)致差異增大。因此，了解坡度提取誤差的成因、大小和空間分布規(guī)律有利于改進(jìn)坡度提取方法，提高坡度提取的精度。

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

上述研究多從理論層面研究坡度提取的不確定性，但很少考慮整個(gè)計(jì)算過(guò)程的方法誤差，以及可能存在的人為誤差。本文針對(duì)現(xiàn)有坡度提取方法，運(yùn)用以Python為開(kāi)發(fā)語(yǔ)言的地理建模技術(shù)，從實(shí)踐層面研究DEM分辨率對(duì)坡度提取的影響，整合分析過(guò)程，提高分析效率，減少人為誤差。研究結(jié)論對(duì)完善現(xiàn)有坡度提取理論具有顯著的理論意義。

1 坡度計(jì)算及不確定性分析方法

1.1 基于DEM的坡度計(jì)算

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

圖1 地表單元坡度示意圖

(1)

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

(2)

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

圖2 三階差分

基于DEM格網(wǎng)計(jì)算fx、fy可以有以下3種不同的計(jì)算方法

(3)

(4)

(5)

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

1.2 坡度計(jì)算的不確定性分析方法

Python作為一種通用的開(kāi)源編程語(yǔ)言，易于學(xué)習(xí)和使用，在ArcGIS地理處理中用作腳本語(yǔ)言。而ArcGIS中的Arcpy站點(diǎn)包給Python處理地理數(shù)據(jù)提供了基礎(chǔ)[11-12]。本文基于Arcpy站點(diǎn)包編寫(xiě)Python腳本程序?qū)EM數(shù)據(jù)進(jìn)行批量重采樣、批量坡度提取及批量掩膜提取。利用numpy包對(duì)數(shù)據(jù)進(jìn)行整理，scipy包對(duì)數(shù)據(jù)進(jìn)行科學(xué)計(jì)算并分析，matplotlib包將數(shù)據(jù)分析結(jié)果制成圖表。

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

2 試驗(yàn)及結(jié)果分析

2.1 試驗(yàn)數(shù)據(jù)及試驗(yàn)流程

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

以30 m空間分辨率的DEM數(shù)據(jù)為基礎(chǔ)數(shù)據(jù)，對(duì)上述7類地貌類型區(qū)坡度提取的不確定性進(jìn)行分析，具體試驗(yàn)流程如下：

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

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

圖4 不同分辨率DEM

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

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

(°)

(3) 運(yùn)用分辨率與平均坡度分析模型進(jìn)行分析計(jì)算，得到以DEM分辨率為橫坐標(biāo)、7種地貌類型區(qū)域的平均坡度值為縱坐標(biāo)的散點(diǎn)圖，并對(duì)散點(diǎn)圖進(jìn)行擬合，通過(guò)回歸算法求解得到回歸方程的系數(shù)及常數(shù)項(xiàng)。

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

圖6 溝壑密度提取模型

2.2 試驗(yàn)結(jié)果分析

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

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

由圖7可知，平均坡度值與分辨率呈現(xiàn)很強(qiáng)的線性關(guān)系，且是正相關(guān)關(guān)系。因此以線性方程來(lái)描述這種關(guān)系，方程式如下

y=ax+b

(6)

式中，x為DEM分辨率；y為平均坡度；a為回歸方程系數(shù)；b為常數(shù)項(xiàng)。則7種地貌類型區(qū)域內(nèi)的平均坡度與分辨率的回歸關(guān)系計(jì)算結(jié)果見(jiàn)表2。

表2 不同地貌單元內(nèi)平均坡度與分辨率的回歸關(guān)系

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

由圖8可見(jiàn)，不同地貌單元內(nèi)的地形因子溝壑密度值與相應(yīng)的回歸方程的常數(shù)項(xiàng)b存在顯著的二次函數(shù)關(guān)系，其二次回歸方程如下

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

(7)

式中，x為不同地貌類型所對(duì)應(yīng)的溝壑密度值；b為不同地貌類型區(qū)域內(nèi)平均坡度與分辨率線性回歸方程的常數(shù)項(xiàng)。由此得到試驗(yàn)區(qū)域內(nèi)平均坡度的變異特征如下

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

(8)

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

圖8 不同地形因子與常數(shù)項(xiàng)b的關(guān)系

假定以分辨率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為試驗(yàn)區(qū)域內(nèi)平均坡度提取誤差?？梢?jiàn)坡度的提取平均誤差只與DEM的空間分辨率有關(guān)。

3 結(jié) 論

針對(duì)現(xiàn)有基于DEM的坡度提取方法，通過(guò)地理建模技術(shù)，從實(shí)踐層面研究DEM分辨率對(duì)坡度提取的影響，整合分析過(guò)程，提高分析效率，減少人為誤差。以天津薊縣為研究區(qū)，采用ArcGIS平臺(tái)和Python腳本語(yǔ)言功能分析了基于DEM的坡度提取方法的不確定性。研究結(jié)果表明：

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

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

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

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

研究結(jié)論進(jìn)一步驗(yàn)證了當(dāng)前基于DEM的坡度提取理論方法的實(shí)用性，是對(duì)現(xiàn)有理論方法的有效補(bǔ)充，具有一定的理論價(jià)值。

參考文獻(xiàn)：

[1] 劉成，徐剛，楊娟，等.激發(fā)坡面泥石流發(fā)生的坡度因素探討——以重慶北碚地區(qū)為例[J].中國(guó)地質(zhì)災(zāi)害與防治學(xué)報(bào)，2006，17(4)：165-168.

[2] 郭芳芳，楊農(nóng)，孟暉，等.地形起伏度和坡度分析在區(qū)域滑坡災(zāi)害評(píng)價(jià)中的應(yīng)用[J].中國(guó)地質(zhì)，2008，35(1)：131-143.

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

[4] 李天文，劉學(xué)軍，湯國(guó)安.地形復(fù)雜度對(duì)坡度坡向的影響[J].山地學(xué)報(bào)，2004，22(3)：272-277.

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

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

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

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

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

[10] 湯國(guó)安，李發(fā)源，劉學(xué)軍.數(shù)字高程模型教程[M].3版.北京:科學(xué)出版社，2016.

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

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

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

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

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

歡迎訂閱《測(cè)繪通報(bào)》

《測(cè)繪通報(bào)》創(chuàng)刊于1955年，是由國(guó)家測(cè)繪地理信息局主管、中國(guó)地圖出版社(測(cè)繪出版社)主辦的反映我國(guó)測(cè)繪地理信息科技發(fā)展現(xiàn)狀和指導(dǎo)全國(guó)測(cè)繪地理信息生產(chǎn)業(yè)務(wù)的綜合性、技術(shù)性刊物?，F(xiàn)為《中文核心期刊要目總覽》收錄期刊、中國(guó)科技核心期刊(中國(guó)科技論文統(tǒng)計(jì)源期刊)、中國(guó)科學(xué)引文數(shù)據(jù)庫(kù)(CSCD)來(lái)源期刊、中國(guó)學(xué)術(shù)期刊網(wǎng)絡(luò)出版總庫(kù)(CAJD)全文收錄期刊、中國(guó)學(xué)術(shù)期刊綜合評(píng)價(jià)數(shù)據(jù)庫(kù)(CAJCED)核心統(tǒng)計(jì)源期刊、《中國(guó)學(xué)術(shù)期刊文摘》收錄期刊。曾被評(píng)為第二、三屆國(guó)家期刊獎(jiǎng)百種重點(diǎn)期刊，多次榮獲全國(guó)優(yōu)秀測(cè)繪地理信息期刊獎(jiǎng)。

《測(cè)繪通報(bào)》的主要內(nèi)容包括大地測(cè)量、全球?qū)Ш叫l(wèi)星系統(tǒng)(GNSS)、攝影測(cè)量、遙感(RS)、地圖制圖、地理信息系統(tǒng)(GIS)、工程測(cè)量、礦山測(cè)量、地籍測(cè)繪、海洋測(cè)繪、測(cè)繪儀器、信息傳輸、圖形圖像處理等方面的新成果和新技術(shù)；行業(yè)管理、科研、教學(xué)、生產(chǎn)的先進(jìn)經(jīng)驗(yàn)；計(jì)算機(jī)、通信等相關(guān)理論技術(shù)在測(cè)繪地理信息領(lǐng)域里的應(yīng)用及測(cè)繪地理信息科技在國(guó)家經(jīng)濟(jì)建設(shè)各個(gè)方面的應(yīng)用；國(guó)內(nèi)外測(cè)繪地理信息學(xué)術(shù)動(dòng)態(tài)及有關(guān)測(cè)繪地理信息科技信息等。

《測(cè)繪通報(bào)》的主要欄目有學(xué)術(shù)研究、技術(shù)交流、測(cè)繪地理信息論壇、經(jīng)驗(yàn)介紹、測(cè)量員之窗、國(guó)外測(cè)繪地理信息、行業(yè)管理、行業(yè)研究、行業(yè)調(diào)查、知識(shí)窗、新書(shū)介紹、企業(yè)之窗、測(cè)繪地理信息教學(xué)等。

《測(cè)繪通報(bào)》為月刊，2018年定價(jià)：35.00元，郵發(fā)代號(hào)：2-223。

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

網(wǎng)址：http:∥tb.sinomaps.com