ArcGIS環(huán)境下基于DEM的地形指標(biāo)和水文信息提取分析

周林虎，張秉來，曹榮泰，張 建，祁兆鑫，劉宇平，胡榮澤，張?chǎng)╆?/p>

(中國電建集團(tuán)青海省電力設(shè)計(jì)院有限公司，青海 西寧 810008）

0 引言

近年來，“3S技術(shù)”不斷發(fā)展成熟，被越來越多的工程設(shè)計(jì)和勘察單位所應(yīng)用，導(dǎo)致地質(zhì)勘察與測(cè)量鄰域也進(jìn)入了數(shù)字時(shí)代，尤其是地理信息系統(tǒng)(geographic information systems，GIS)和數(shù)字高程模型(digital elevation model，DEM)在地質(zhì)勘察方面得到了廣泛應(yīng)用[1]。地理信息，即與空間地理分布有關(guān)的信息，具有區(qū)域性、多維性和動(dòng)態(tài)性等特點(diǎn)，是指地表物體和環(huán)境固有的數(shù)量、質(zhì)量、分布特征，以及聯(lián)系和規(guī)律的數(shù)字、文字、圖形、圖像等的總稱[2]。DEM即指利用有限的高程數(shù)據(jù)實(shí)現(xiàn)地形三維可視化模擬，是一種實(shí)體地面模型，為數(shù)字地形模型(digital terrain model，DTM)的一個(gè)分支[3]。

基于DEM的地表特征提取和洪水淹沒模擬技術(shù)主要是根據(jù)DEM提供的高程信息，對(duì)研究區(qū)的高程、坡度變率、坡向變率、地形起伏度、地面粗糙度等基本地表特征信息，以及山脊線、山谷線和山頂點(diǎn)等地形特征參數(shù)進(jìn)行提取和分析，以專題圖的形式進(jìn)行矢量或柵格數(shù)據(jù)的輸出，該技術(shù)不僅可以對(duì)研究區(qū)地形地貌分布情況與變化情況進(jìn)行描述和分析，還可以對(duì)區(qū)域水土流失、滑坡、崩塌和泥石流等地質(zhì)災(zāi)害進(jìn)行評(píng)估，提出合理的避讓和防治措施。

另外，通過模擬洪災(zāi)發(fā)生時(shí)的洪水淹沒范圍和深度，可預(yù)先獲知和劃分區(qū)域洪水受災(zāi)程度，為洪災(zāi)評(píng)估、預(yù)警和防治提供參考[3]。對(duì)洪水淹沒情況進(jìn)行分析亦是防汛抗洪工作的重要內(nèi)容，是當(dāng)前GIS 技術(shù)在水利應(yīng)用鄰域的研究熱點(diǎn)。早期的洪水淹沒分析大多以水動(dòng)力學(xué)模型為主，但水動(dòng)力學(xué)模型具有較強(qiáng)的針對(duì)性，且構(gòu)建復(fù)雜、所需參數(shù)多，無法進(jìn)行廣泛推廣。隨著對(duì)洪水淹沒分析的深入研究，不少學(xué)者發(fā)現(xiàn)洪水淹沒程度主要取決于洪水特征以及地形因子，用DEM 能夠較好地反映地表起伏狀態(tài)[4]。特別是近些年，隨著暴雨和強(qiáng)降雨的不斷加劇，我國多地洪澇災(zāi)害數(shù)量增多、程度加重，因此對(duì)區(qū)域洪水進(jìn)行模擬，科學(xué)、準(zhǔn)確地預(yù)測(cè)洪水淹沒范圍和洪水深度分布對(duì)防洪減災(zāi)具有重要的現(xiàn)實(shí)意義[5-6]。

基于此，以青海省某一區(qū)域?yàn)檠芯繀^(qū)，利用ArcGIS軟件提取和分析了該區(qū)域基本地表特征信息，生成了相關(guān)專題圖件；并在此基礎(chǔ)上進(jìn)行了研究區(qū)洪水淹沒分析，為該區(qū)域的城鎮(zhèn)發(fā)展、農(nóng)業(yè)規(guī)劃、資源開發(fā)、災(zāi)害預(yù)測(cè)及治理等提供指導(dǎo)性建議。

1 研究區(qū)概況

本研究區(qū)海拔為2 071.96～2 344.55 m，高低起伏較大，最大高差達(dá)272.59 m。區(qū)內(nèi)地層以黃土和泥巖為主，滑坡、泥石流和崩塌等地質(zhì)災(zāi)害發(fā)育較多，尤其6～9月份降雨較多，易發(fā)生地質(zhì)災(zāi)害。該區(qū)域地層裸露較多，地表植被相對(duì)較少，為水土流失較嚴(yán)重的地區(qū)。長期而頻繁的水土流失、滑坡和泥石流等地質(zhì)災(zāi)害的發(fā)生，使得該區(qū)域土地退化嚴(yán)重，生態(tài)環(huán)境不斷惡化，嚴(yán)重影響和阻礙了該區(qū)域工程設(shè)施的建設(shè)和社會(huì)經(jīng)濟(jì)的可持續(xù)發(fā)展，因此對(duì)該地區(qū)地形地貌特征進(jìn)行分析，對(duì)地質(zhì)災(zāi)害進(jìn)行調(diào)查和評(píng)價(jià)，并提出合理的防治和治理措施具有重要意義。

2 傾斜攝影測(cè)量技術(shù)

1)優(yōu)點(diǎn)和應(yīng)用前景

在以往的工程地質(zhì)災(zāi)害勘察工作中，地質(zhì)災(zāi)害調(diào)查分析多以勘察員的工程經(jīng)驗(yàn)判斷為主，受現(xiàn)場(chǎng)視野的局部性、地形、現(xiàn)場(chǎng)復(fù)雜條件、勘察手段和人員安全因素等限制，勘察人員往往無法到達(dá)現(xiàn)場(chǎng)，野外實(shí)施過程中地質(zhì)災(zāi)害點(diǎn)要素描寫欠全面、規(guī)范，災(zāi)害成因分析不足，資料收集不夠完整，地質(zhì)災(zāi)害勘察結(jié)果只是定性分析得出的結(jié)果，而無法得出定量分析結(jié)果。利用無人機(jī)傾斜攝影測(cè)量技術(shù)，通過獲得的航空影像和DEM數(shù)據(jù)，可準(zhǔn)確地對(duì)研究區(qū)坡度、坡向、流向和流量等地表特征進(jìn)行分析和計(jì)算，從而調(diào)查地質(zhì)災(zāi)害點(diǎn)區(qū)域自然環(huán)境特征和孕育地質(zhì)災(zāi)害的地質(zhì)背景，查明地質(zhì)災(zāi)害的空間分布現(xiàn)狀和強(qiáng)度等發(fā)育特征，得出地基穩(wěn)定性評(píng)價(jià)，節(jié)約成本，提質(zhì)增效。

2)關(guān)鍵技術(shù)難點(diǎn)

近幾年，傾斜攝影測(cè)量技術(shù)得到了迅速發(fā)展，傾斜攝影測(cè)量技術(shù)能夠獲取建筑物、樹木等地理實(shí)體的紋理細(xì)節(jié)，不但豐富了影像數(shù)據(jù)源信息，同時(shí)高冗余度的航攝影像重疊為高精度的影像匹配提供了條件，使得基于人工智能的三維實(shí)體重建成為了可能。分層顯示技術(shù)、紋理映射技術(shù)成為傾斜攝影測(cè)量和建模的關(guān)鍵支撐點(diǎn)，極大地提升了三維建模的效率，同時(shí)也降低了建模生產(chǎn)成本。目前，基于傾斜攝影測(cè)量成果的應(yīng)用還比較少，因此，應(yīng)用創(chuàng)新點(diǎn)挖掘還需要逐步深入。同時(shí)不能局限于傾斜攝影測(cè)量技術(shù)，還需要研究其與雷達(dá)、紅外、多光譜、高光譜等多種傳感器的結(jié)合，集成在更小的無人機(jī)上拓寬攝影測(cè)量技術(shù)的應(yīng)用范圍。未來，基于點(diǎn)云數(shù)據(jù)計(jì)算的大規(guī)模三維數(shù)據(jù)生產(chǎn)使得工程測(cè)量、三維建模等工作進(jìn)行變革，開啟三維遙感的新時(shí)代。

3 地形指標(biāo)提取

地形指標(biāo)主要包括坡度變率、坡向變率、地形起伏度和地面粗糙度等4個(gè)基本地形指標(biāo)，即最基本的自然地理要素，同時(shí)也是對(duì)人類生產(chǎn)、生活影響最大的自然要素[7]。地形特征值廣泛應(yīng)用于諸多研究和應(yīng)用鄰域，地形指標(biāo)的提取對(duì)水土流失、土地利用、土地資源評(píng)價(jià)、城市規(guī)劃等方面的研究起著重要作用。

3.1 坡度變率

坡度變率，又名剖面曲率，是指坡度在微分空間上的變化率，地表面高程相對(duì)于水平面變化的二階導(dǎo)數(shù)，即坡度之坡度(slope of slope，SOS)[8]。其能很好地反應(yīng)出地面高程變化率的大小，從而反應(yīng)出剖面曲率信息。根據(jù)ArcGIS表面分析工具，提取出的坡度和坡度變率如圖1所示。

圖1 研究區(qū)坡度和坡度變率

對(duì)于坡度和坡度曲率的研究和分析不僅能直觀判斷和分析研究區(qū)地形高低起伏情況、微地形地貌，對(duì)于土地利用規(guī)劃和資源開發(fā)具有參考價(jià)值，還能作為衡量和評(píng)價(jià)區(qū)域物質(zhì)能量轉(zhuǎn)移、水土流失、土壤侵蝕以及洪澇災(zāi)害的重要基礎(chǔ)數(shù)據(jù)，對(duì)于地質(zhì)災(zāi)害的評(píng)價(jià)和防治具有重要意義[9]。在我國，坡度為25°是水土保持法規(guī)定的臨界坡度[10]，即坡度大于25°的地區(qū)，其水土保持能力相對(duì)較弱，相同條件下發(fā)生水土流失、土壤侵蝕、崩塌、滑坡和泥石流等地質(zhì)災(zāi)害的概率亦較大。通過分析研究區(qū)坡度數(shù)據(jù)分布情況，得到區(qū)內(nèi)坡度大于25°的區(qū)域面積占總面積的49.86%，由此可知該流域內(nèi)一半左右區(qū)域坡度較大，水土保持能力較弱，且結(jié)合現(xiàn)場(chǎng)植被分布情況可知，該部分區(qū)域植被覆蓋率較低，水土流失情況嚴(yán)重，區(qū)內(nèi)存在較多沖溝，加之豐水季節(jié)雨水較多，因此判定該區(qū)域發(fā)生地質(zhì)災(zāi)害的概率較大。為防止水土流失進(jìn)一步加劇，需采取專門的水土保持措施。通過坡度變率圖可知，圖中紅色高值即坡度圖中坡度變化較顯著、落差較大的區(qū)域，即地勢(shì)急劇變化的區(qū)域。通過計(jì)算得到，該區(qū)域坡度變率為80～90的區(qū)域占比82.75%，見表1所列，說明該區(qū)域地形起伏較大，坡度變化較顯著，多為山地、丘陵和溝谷地貌，不適宜發(fā)展大規(guī)模機(jī)械化種植業(yè)，可發(fā)展林業(yè)、山地畜牧業(yè)以及林下種植業(yè)。

表1 不同坡度變率面積占比

3.2 坡向變率

坡向變率，又名平面曲率，是指坡向在微分空間上的變化率，對(duì)提取的坡向數(shù)據(jù)再進(jìn)行一次坡度提取，即坡向之坡度(slope of aspect，SOA)，其可以很好地反映等高線的彎曲程度[8]。坡向和坡向變率提取方法如下：①通過坡向工具提取坡向數(shù)據(jù)DEM1；②提取坡向變率SOA1，即對(duì)坡向數(shù)據(jù)DEM1提取坡度，得到坡向變率SOA1；③提取反地形數(shù)據(jù)FDEM通過柵格計(jì)算器工具，輸入計(jì)算式為FDEM=2 344.55-DEM，2 344.55 m即本流域最大高程；④提取反地形坡向變率SOA2，即對(duì)FDEM提取反地形坡向FDEM1，然后對(duì)坡向數(shù)據(jù)FDEM1提取坡度，得到反地形坡向變率SOA2；⑤計(jì)算SOA。使用柵格計(jì)算器工具，輸入式為SOA=[(SOA1+SOA2)-ABS(SOA1-SOA2)]/2，即可得到無誤差的DEM坡向變率。

提取出的坡向和坡向變率如圖2～圖4所示，通過SOA可知，其與坡度變率分布情況基本相似，但是通過計(jì)算不同坡度變率區(qū)域面積占比情況得出，坡向變率為0～10和80～90的區(qū)域面積占比為8.29%和88.68%，均高于坡度變率為0～10和80～90的面積占比，而坡向變率為10～80的區(qū)域面積總體占比低于坡度變率為10～80的面積，見表2所列。

表2 不同坡向變率面積占比

圖2 研究區(qū)DEM1和SOA1

圖3 研究區(qū)FDEM1和SOA2

圖4 研究區(qū)SOA

3.3 地形起伏度

地形起伏度是描述某一特定區(qū)域內(nèi)高程極差(即最高點(diǎn)和最低點(diǎn)之間高差)的地形指標(biāo)[11-12]。鄰域范圍可根據(jù)研究區(qū)具體情況和實(shí)際需求進(jìn)行設(shè)置，鄰域范圍越大，求得的起伏度亦越大。本項(xiàng)研究將鄰域范圍設(shè)為12×12，而研究區(qū)像元柵格總數(shù)為200 970 931，即將研究區(qū)分為1 395 631個(gè)矩形柵格，求解這些矩形柵格像元內(nèi)的像元值極差，統(tǒng)計(jì)結(jié)果即該區(qū)域內(nèi)的相對(duì)高差(即起伏度)。通過ArcGIS焦點(diǎn)統(tǒng)計(jì)工具可直接生成起伏度成本數(shù)據(jù)，如圖5所示。由該圖可知，地形起伏度范圍為0～118.75 m，其中起伏度為0～13.19 m的區(qū)域面積占98.946 3%，而隨著地形起伏度不斷增大，所占面積比例也越來越小，見表3所列。

表3 不同地形起伏度面積占比

圖5 研究區(qū)地形起伏度

3.4 地面粗糙度

地面粗糙度是反映區(qū)域表面形態(tài)的一個(gè)宏觀指標(biāo)，指某一特定區(qū)域內(nèi)地球表面積與其投影面積之比[13]。根據(jù)地面粗糙度的定義，求每個(gè)柵格單元的表面積與其投影面積之比。如圖6所示，假如△ABC是一個(gè)柵格單元的縱剖面，a為此柵格單元的坡度，則AB面的面積(AB×AC)為此柵格的表面積，AC面的面積(AC×AC)為此柵格的投影面積(也既是此柵格分辨率平方)，此柵格單元的地面粗糙度M=“AB面的面積”/“AC面的面積”=(AC×AB)/(AC×AC)=AB/AC=1/cosa。根據(jù)ArcGIS表面分析工具計(jì)算出研究區(qū)地面粗糙度如圖7所示，可知研究區(qū)地面粗糙度為1～602.26，但大部分區(qū)域介于1～67.81，隨著地面粗糙度不斷增大，所占面積比例亦逐漸減少，這和地形起伏度表現(xiàn)出相同特征，見表4所列。

表4 不同地面粗糙度面積占比

圖6 △ABC

圖7 地面粗糙度

3.5 山脊線和山谷線

特征地形要素主要是指對(duì)地形在地表的空間分布特征具有控制作用的點(diǎn)、線或面狀要素[14]。特征地形要素構(gòu)成地表地形與起伏變化的基本框架，與地形指標(biāo)的提取主要采用小范圍鄰域分析不同的是，特征地形要素的提取更多地應(yīng)用較為復(fù)雜的技術(shù)方法，如山谷線、山脊線等的提取采用了全局分析法，成為柵格數(shù)據(jù)地學(xué)分析中很具特色的數(shù)據(jù)處理內(nèi)容[15]。

山脊線和山谷線的提取步驟為：①點(diǎn)擊初始DEM數(shù)據(jù)，使用焦點(diǎn)統(tǒng)計(jì)工具，設(shè)置統(tǒng)計(jì)類型為平均值MEAN，鄰域的類型為矩形(也可以為圓)，鄰域的大小為12×12(根據(jù)實(shí)際情況而定)，則可得到一個(gè)鄰域?yàn)?2×12的矩形的平均值數(shù)據(jù)層，記為MEAN，如圖8(a)所示；②使用柵格計(jì)算器計(jì)算，輸入公式ZFDX=DEM-MEAN，即可求出正負(fù)地形分布區(qū)域，如圖8(b)所示；③采用柵格計(jì)算器工具，輸入公式shanji=(ZFDX＞0)&(SOA＞80)，即可提取出山脊線，如圖9(a)所示；④同理，使用柵格計(jì)算器工具，輸入公式shangu=(ZFDX＜0)&(SOA＞80)，即可求出山谷線如圖9(b)所示。通過圖9可看出研究區(qū)地形地貌的山梁、山坡、溝谷和洼地等分布情況，從而直觀分析出區(qū)內(nèi)匯流和分流情況，為后續(xù)洪水淹沒分析結(jié)果的可靠性和真實(shí)性提供參考。

圖8 MEAN和ZFDX

圖9 研究區(qū)山脊線和山谷線

4 洪水淹沒分析

判別給定降水條件下某一區(qū)域的潛在淹沒區(qū)，既可以規(guī)避盲目建設(shè)所帶來的雨洪淹沒風(fēng)險(xiǎn)，也能為生態(tài)保護(hù)與修復(fù)提供可選的項(xiàng)目突破點(diǎn)，因而在區(qū)域規(guī)劃設(shè)計(jì)中具有重要作用。洪水淹沒分析一般分為無源淹沒和有源淹沒2種分析方法，洪水淹沒分析是進(jìn)行水文預(yù)測(cè)預(yù)報(bào)、洪水災(zāi)害評(píng)估的一項(xiàng)重要內(nèi)容，也是GIS地形三維仿真系統(tǒng)中的一個(gè)重要功能[16-17]，前者是指低于某一高程的所有柵格均被劃入淹沒區(qū)，相當(dāng)于區(qū)域大面積均勻降水，所有低洼處都有可能積水成災(zāi)；后者考慮了柵格之間的物理連通性，因而某些柵格受到地形阻擋而不會(huì)被淹沒，且除降水外還考慮地表水，如河流、水庫等的外溢現(xiàn)象[18]。與有源淹沒相比，無源淹沒方法更為簡(jiǎn)單，且在受地表水流動(dòng)影響較小的地區(qū)能實(shí)現(xiàn)良好的模擬效果，因而得到廣泛應(yīng)用。由于本項(xiàng)研究區(qū)周圍未發(fā)現(xiàn)能形成洪水的水庫和河流，洪水基本上由當(dāng)?shù)貜?qiáng)降雨引起，因此洪水淹沒模擬采用無源淹沒方法。

由于本項(xiàng)研究區(qū)高程差較大，高低起伏較顯著，為直觀顯示出不同淹沒深度條件下洪水淹沒區(qū)域和面積占比，本次分析設(shè)置了10 m、50 m、100 m、150 m和200 m等5個(gè)梯度。在實(shí)際分析和模擬中，則應(yīng)根據(jù)實(shí)際高程數(shù)據(jù)和地形情況設(shè)置合理的淹沒深度進(jìn)行分析。表5為不同淹沒深度條件下淹沒面積占比，由該表可知，隨著淹沒深度不斷增大，淹沒面積亦顯著增加。當(dāng)淹沒模擬深度為10 m時(shí)，淹沒面積僅為0.62%，該區(qū)域主要集中于研究區(qū)南部一小部分區(qū)域，隨著淹沒模擬深度增大100 m，淹沒區(qū)域逐漸擴(kuò)散至南部大部分區(qū)域和西北部，淹沒面積達(dá)到35.38%，由此亦可分析出本項(xiàng)研究區(qū)高程由低至高依次為南部、西北部和中東部，與實(shí)際地形相符，如圖10所示。

表5 不同淹沒深度條件下淹沒面積占比

圖10 不同淹沒深度條件下淹沒面積示意圖

值得討論的是，本項(xiàng)洪水淹沒分析對(duì)于相互聯(lián)通的小面積區(qū)域具有良好的模擬效果，但對(duì)于高程變化較大的大面積區(qū)域，由于有許多分水嶺和匯水溝的存在，一般其內(nèi)部并非完全貫通，因此不能形成完整的水域，與實(shí)際情況會(huì)略有偏差。

5 結(jié)語

基于ArcGIS的3D Analyst、表面分析和柵格計(jì)算分析模塊，提取和分析了研究區(qū)坡度、坡度變率、坡向、坡向變率、地形起伏度和地面粗糙度等地形變化指標(biāo)，以及山脊線和山谷線等特征地形要素，從而得到了區(qū)內(nèi)地形高低起伏情況和微地形地貌特征，在此基礎(chǔ)上通過ArcGIS的柵格計(jì)算工具，提取了不同深度洪水情況下的淹沒面積和洪水位，結(jié)合區(qū)域地形情況，驗(yàn)證了洪水淹沒分析的真實(shí)性和可靠性。該方法不需要大量的現(xiàn)場(chǎng)勘察和搜集資料工作即可得到研究區(qū)整個(gè)地形地貌、高程分布情況和洪水淹沒情況，節(jié)省了大量時(shí)間和人力。