采用ArcGIS平臺的地勢起伏度自動提取技術(shù)研究

王 康，何俊仕，于德浩，龍 凡，李 霞

(1.沈陽農(nóng)業(yè)大學(xué) 水利學(xué)院，遼寧 沈陽 110866;2.沈陽軍區(qū)司令部工程科研設(shè)計院，遼寧 沈陽 110162)

隨著3S技術(shù)的快速發(fā)展，利用DEM(digital elevation model)數(shù)字高程模型和遙感影像等數(shù)據(jù)源進(jìn)行地貌特征信息的提取已成為地貌研究的主要手段。其核心思想是根據(jù)相對起伏度或海拔高度來快速地確定基本地貌形態(tài)類型［1－4］。從理論上分析，只要DEM數(shù)據(jù)足夠精確，就可以很準(zhǔn)確地確定基本地貌形態(tài)。因此對于地勢起伏度的提取，關(guān)鍵在于一定統(tǒng)計窗口的確定。即確定面積有多大，才能恰到好處地反映山體的完整性，并具有一定范圍內(nèi)較強(qiáng)的代表性，亦即普適性［5］。劉振東［6］提出按照地貌發(fā)育的基本理論，存在一個使最大高差達(dá)到相對穩(wěn)定的最佳統(tǒng)計窗口，這個最佳統(tǒng)計窗口即為地勢起伏度的最佳統(tǒng)計單元。一般采取“拐點(diǎn)”理論［5］來確定窗口大小，這里所定義的“拐點(diǎn)”與數(shù)學(xué)上的拐點(diǎn)不同，這里“拐點(diǎn)”的定義為，高差由大變小的點(diǎn)為“拐點(diǎn)”。做統(tǒng)計時會發(fā)現(xiàn)這樣一條規(guī)律:那就是隨著統(tǒng)計窗口的擴(kuò)大，窗口中最高點(diǎn)與最低點(diǎn)的高差不斷增加，開始增加的很快，當(dāng)高差增加到一定值時，隨著統(tǒng)計窗口的再增加高程差增加不大，而這樣的點(diǎn)，則為要尋找的“拐點(diǎn)”。這個“拐點(diǎn)”所對應(yīng)的面積，即為最終的統(tǒng)計單元，而“拐點(diǎn)”所對應(yīng)的高差值即為地勢起伏高度。

利用ArcGIS軟件，加載Spatial Analyst模塊，借助Neighborhood Statistic功能，對工作區(qū)以3×3起始窗口，到90×90終止窗口，移動步距多為1，計算窗口內(nèi)最大高程值MAXIMUM和最小高程值MINIMUM［7］。然后利用 Raster Calculator工具計算專題層的MAXIMUM與MINIMUM的差值即為地勢起伏度層，其每個柵格的值即為以該柵格為中心的確定領(lǐng)域的地面起伏度值。最后，利用窗口大小與最大高程差值生成所需圖表文件，并依此尋求拐點(diǎn)位置。

上述方法的確可行，但是此方法需依次重復(fù)從3×3窗口到90×90窗口的領(lǐng)域計算等操作，該過程工作繁瑣，效率低，中間生成大量緩存文件，而且大量的繁瑣操作容易出現(xiàn)人為操作失誤，導(dǎo)致最終數(shù)據(jù)精度不可靠。

于德浩等［8］學(xué)者曾利用遙感特征指數(shù)對地表水體實(shí)現(xiàn)了自動化提取?；诖顺晒Π咐?，本文針對常規(guī)方法提取地勢起伏度操作繁瑣、效率低，以及提取過程生成大量緩存文件等問題，提出一種自動提取的方法。該方法對傳統(tǒng)“拐點(diǎn)”理論手工確定窗口大小的過程加以改進(jìn)，利用計算機(jī)編程實(shí)現(xiàn)不同統(tǒng)計單元起伏度值的自動計算，并依此尋求最佳統(tǒng)計單元，實(shí)現(xiàn)了地勢起伏度的自動提取。

1 地勢起伏度及其作用

地勢起伏度(Relief Amplitude)，也稱為地形起伏度或地勢能量、相對高度。在我國國土大地貌研究和1:400萬地貌圖工作中，將地勢起伏度定義為:一個山體的山脊最高點(diǎn)與順?biāo)鞣较蜃罱恿鞯牡谝粋€匯流處為基準(zhǔn)面的最大高差［9］。在實(shí)際操作中，涂漢明等［10］又將某點(diǎn)的地勢起伏度RAi定義為某一確定面積(A)中的最高點(diǎn)(Zimax)與最低點(diǎn)(Zimin)之高差。它是描述一個區(qū)域地形宏觀性指標(biāo)之一，在宏觀的區(qū)域內(nèi)反映地面的起伏特征。用公式表示如下:

式中:RAi指某區(qū)域的地勢起伏度;Zimax指該區(qū)域內(nèi)的最大高程值;Zimin則指該區(qū)域內(nèi)的最小高程值;i為自然數(shù)，指某一統(tǒng)計區(qū)域。

坡度和起伏度是地形描述中最常用的參數(shù)，它們能快速而直觀地反映地勢起伏特征。其中，坡度是劃分平原和非平原的重要依據(jù)之一;而地勢起伏度又可進(jìn)一步劃分為平原、臺地、丘陵、小起伏山地、中起伏山地、大起伏山地和極大起伏山地等7級［11］。地形起伏是構(gòu)造作用與地表剝蝕過程相互作用的結(jié)果，是地貌研究的必要手段［1］。地勢起伏度的提取為獲取地表信息，進(jìn)行地形定量化分析提供有效手段。

2 地勢起伏度自動提取的實(shí)現(xiàn)

2.1 自動計算3×3窗口到90×90窗口的最大值與最小值的差值

應(yīng)用ArcGIS的Python模塊進(jìn)行批處理，輸入代碼:

在工作空間目錄下新建一個文件夾命名為Output，用來保存輸出結(jié)果。

運(yùn)行結(jié)束之后會在output文件夾下面自動生成r_3至r_90的柵格數(shù)據(jù)集。操作界面如圖1所示。

圖1 Python操作界面

2.2 讀取所生成數(shù)據(jù)集的起伏度值(最大值屬性)

在事先準(zhǔn)備好的Excel中A1輸入“r_3”;B1輸入“MAXIMUM”。通過下拉生成88行如圖2所示。

將表格內(nèi)容復(fù)制，然后應(yīng)用ArcGIS10在Arc-Toolbox中選擇“獲取柵格屬性”工具，選擇“批處理”。添加行，添加至88行。將88行全部選中然后粘貼可簡化對設(shè)置每行的數(shù)據(jù)輸入和所獲取屬性等的復(fù)雜操作，之后在“環(huán)境…”中把當(dāng)前工作空間和臨時工作空間都設(shè)置為:E:/exe/output。

圖2 Excel所生成內(nèi)容

確定執(zhí)行后返回的MAXIMUM屬性將顯示在地理處理移動窗口中，見圖3。

圖3 消息窗口

右擊“消息”，將其復(fù)制，粘貼在Word中。

將第一行“消息”、第三行“開始時間……”和最后一行“成功在……”等內(nèi)容刪掉。將“執(zhí)行:GetRasterProperties E:/exe/output/”、 “ MAXIMUM”和“最大值 =”均替換為空，之后轉(zhuǎn)換成2列的表格。粘貼到Excel中，便可根據(jù)要求生成圖表來尋求“拐點(diǎn)”。

提取結(jié)果見圖4、圖5。從圖中均可看到在一定的范圍內(nèi)隨著移動窗口的擴(kuò)大，起伏度也在不斷擴(kuò)大，但當(dāng)窗口達(dá)到一定程度后起伏度增長已十分緩慢，之后隨著移動窗口的增大起伏度已趨于穩(wěn)定。這里的“拐點(diǎn)”即為所要尋求的最佳點(diǎn)——地勢起伏度。

通過對測試區(qū)提取結(jié)果進(jìn)行分類，劃出地勢起伏度分級圖，見圖6a，較符合實(shí)際情況;同時，對局部區(qū)域放大并結(jié)合等高線分析，見圖6b，地勢起伏度較大地區(qū)集中在等高線密集的山區(qū)地帶，而地勢起伏度較小的地區(qū)主要集中在等高線稀疏的河谷等平緩地帶，與實(shí)際情況基本吻合。

圖4 試驗(yàn)區(qū)I最佳統(tǒng)計單元提取結(jié)果

圖5 試驗(yàn)區(qū)II最佳統(tǒng)計單元提取結(jié)果

圖6 地勢起伏度提取結(jié)果

3 結(jié)論

本文提出了一種地勢起伏度自動提取的方法，對傳統(tǒng)的利用“拐點(diǎn)”理論確定最佳統(tǒng)計單元的過程加以改進(jìn)，通過編程實(shí)現(xiàn)了不同統(tǒng)計單元起伏度值的自動計算，并依此確定最佳統(tǒng)計單元，實(shí)現(xiàn)了地勢起伏度的自動提取。與常規(guī)方法比較，此方法具有如下優(yōu)點(diǎn):

(1)操作簡便，效率高:常規(guī)方法需要重復(fù)從3×3窗口到90×90窗口的鄰域計算等操作，操作復(fù)雜，工作量大，本方法通過Python模塊編程能夠?qū)崿F(xiàn)此過程的自動化運(yùn)行，減少了繁瑣而復(fù)雜的操作，效率得到了明顯提高。

(2)減少大量不必要緩存文件:常規(guī)方法計算中，首先求算分析窗口內(nèi)DEM數(shù)據(jù)的最大與最小高程值，分別記作專題層Max與Min，之后利用Raster Calculater工具計算專題層Max與Min的差值，得到地勢起伏度專題層。這意味著將要生成三倍于原始文件大小的緩存文件，如果重復(fù)計算3×3窗口到90×90窗口的鄰域計算則需要生成264倍于原始文件大小的緩存文件。本研究一步到位，能夠直接計算出起伏度值專題層，可使緩存文件縮小到常規(guī)方法的1/3。

該研究有利于地勢起伏度研究工作的快速開展，可為該研究提供有效的技術(shù)支撐，同時也必將促進(jìn)地勢起伏度研究工作的快速發(fā)展。此外，由于地勢起伏度值計算過程均是計算機(jī)自動計算，因此難免出現(xiàn)許多面積相對較小、很難在圖上得以顯示的破碎圖斑，這也是以后工作中有待提高的方面。

［1］A Kühni，O.A Pfiffener.The relief of the Swiss Alps and adjacent areas and structure:topographic analysis from a250-m DEM［J］.Geomorphology，2001，41:285－307.

［2］G Onorati，R Ventura，M Poscolieri，et al.The Digital Elevation Model of Italy for geomorphology and structural geology［J］.CATENA，1992，19(2):147 －178.

［3］閭國年，錢亞東，陳鐘明.基于柵格數(shù)字高程模型提取特征地貌技術(shù)研究［J］.地理學(xué)報，1998，53(6):562－569.

［4］Balázs，Székely，Dávid Karátson.DEM-based morphometry as a tool for reconstructing primary volcanic landforms:examples from the Borzsony Mountains，Hungary［J］.Geomorphology，2004，63:25 －37.

［5］涂漢明，劉振東.中國地勢起伏度最佳統(tǒng)計單元的求證［J］.湖北大學(xué)學(xué)報(自然科學(xué)版)，1990，12(3):266－271.

［6］劉振東.利用DTM編制小比例尺地勢起伏度圖的初步研究［J］.測繪學(xué)報，1990，19(1):57－62.

［7］中國科學(xué)院地理科學(xué)與資源研究所資源與環(huán)境信息系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室.中國1∶100萬數(shù)字地貌遙感解譯與集成流程［S］.(試行本).2005.

［8］于德浩，龍凡，鄧正棟，等.基于遙感特征指數(shù)的地表水體自動提取技術(shù)研究［J］.地球物理學(xué)進(jìn)展，2009(2):707－713.

［9］劉振東，涂漢明.中國地勢起伏度統(tǒng)計單元的初步研究［J］.熱帶地理，1989，9(1):31 －38.

［10］涂漢明，劉振明.中國地勢起伏度研究［J］.測繪學(xué)報，1991，20(4):311 －319.

［11］中科院地理科學(xué)與資源研究所(資源與環(huán)境信息系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室).中華人民共和國1:100萬數(shù)字地貌圖制圖規(guī)范［Z］.(征求意見稿).2005.