4

(1.河海大學(xué) 水文水資源學(xué)院，南京 210098； 2.水利部海河水利委員會(huì) 水文局，天津 300170； 3.南京水利科學(xué)研究院水工水力學(xué)研究所，南京 210029； 4.河海大學(xué) 水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210098； 5.中國(guó)科學(xué)院南京土壤研究所 土壤物理與鹽漬土研究室，南京 210008)

1 研究背景

數(shù)字高程模型(Digital Elevation Model，DEM)使用有限的數(shù)據(jù)點(diǎn)近似模擬地表形態(tài)，是地貌、水文、土壤、生態(tài)等學(xué)科獲取地形信息的數(shù)據(jù)基礎(chǔ)[1]。在實(shí)際應(yīng)用中，DEM的水平分辨率會(huì)對(duì)地形特征分析的結(jié)果產(chǎn)生影響，低分辨率DEM存在較多的地形信息遺失，使描述對(duì)象失真[2-4]。為保證成果精度，土壤制圖[5]、水土流失評(píng)價(jià)[6-7]等需要精確地表輸移路徑的領(lǐng)域?qū)EM分辨率要求極高。

青藏高原是目前水文研究的熱點(diǎn)地區(qū)。受地表地形起伏、植被覆蓋、云層等的影響[8]，目前廣泛使用的SRTM和ASTER GDEM這些DEM數(shù)據(jù)在該地區(qū)存在較大誤差[9-11]。此外，SRTM和ASTER GDEM數(shù)據(jù)的空間分辨率僅為1″(約30 m)，這樣的分辨率難以滿足小流域水文與地形研究的需要[12]。因此，青藏高原區(qū)的水文研究需要高精度DEM數(shù)據(jù)的獲取方式。然而，DEM分辨率的提升雖然可以提高所代表地形信息精度，但同樣會(huì)增加相關(guān)水文模型的計(jì)算時(shí)間，所以水文研究往往需要選取能適當(dāng)保障地形信息精度的最優(yōu)分辨率[12]。

激光探測(cè)與測(cè)量系統(tǒng)(Light Detection And Ranging，LiDAR)作為一種先進(jìn)的三維空間信息獲取手段，是目前獲取高精度、高分辨率DEM的主要方法[13-16]。本文選取拉薩河奪底溝小流域一處山坡為研究對(duì)象，通過(guò)LiDAR點(diǎn)云生成DEM，使用地形剖面法與SRTM數(shù)據(jù)及ASTER GDEM數(shù)據(jù)對(duì)比，評(píng)價(jià)LiDAR DEM質(zhì)量。此外，通過(guò)比較不同分辨率LiDAR DEM提取的坡度、平面曲率、剖面曲率，研究高山區(qū)地形屬性與分辨率的關(guān)系，選取最適合水文研究的DEM分辨率。

2 研究資料與方法

2.1 研究區(qū)概況

研究區(qū)選自拉薩市區(qū)以北的奪底鄉(xiāng)境內(nèi)，位于91°11′07″E—91°11′47″E，29°46′19″N—29°46′43″N之間，面積約0.46 km2，屬拉薩河流域奪底溝子流域的一級(jí)子流域。該區(qū)域地形起伏變化明顯，東北高，西南低，內(nèi)有一處山脊將其劃分為2道溝谷，海拔分布見(jiàn)圖1。研究區(qū)土壤類(lèi)型為砂土，地表覆蓋為高山草甸和裸巖，近年無(wú)大規(guī)模人類(lèi)活動(dòng)或地質(zhì)作用等改變地表形態(tài)的事件發(fā)生。由于無(wú)高株植被和人工建筑物等的遮擋且地形穩(wěn)定，適合采集LiDAR點(diǎn)云進(jìn)行地形分析，并與近20 a來(lái)獲取的遙感地形數(shù)據(jù)進(jìn)行對(duì)比研究。

圖1 研究山坡位置及現(xiàn)場(chǎng)Fig.1 Map and scenic picture of the hillside slope

2.2 地形數(shù)據(jù)

2.2.1 LiDAR地形數(shù)據(jù)

LiDAR利用光速恒定原理獲取目標(biāo)物距離，即使用傳感器發(fā)射激光束至目標(biāo)物表面，激光束經(jīng)反射后由傳感器接收，系統(tǒng)精確記錄發(fā)射與接收激光束的時(shí)刻，使用式(1)計(jì)算得到激光束到達(dá)的目標(biāo)點(diǎn)與傳感器間的斜距R。

(1)

式中：c為光速；t為接收激光束時(shí)刻與發(fā)射激光束時(shí)刻的時(shí)間差。使用測(cè)得的斜距R結(jié)合掃描儀坐標(biāo)、傳感器高度以及掃描角度信息可以計(jì)算得到激光束到達(dá)的每一個(gè)地面點(diǎn)的x，y，z三維空間坐標(biāo)，大量離散分布的三維點(diǎn)組成地形點(diǎn)云。DEM數(shù)據(jù)可通過(guò)地形處理軟件對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行插值處理得到。

在本研究中，地形測(cè)量選用Leica-HDS8800掃描測(cè)量系統(tǒng)，如圖2所示。該系統(tǒng)最遠(yuǎn)掃描距離可達(dá)2 000 m，掃描水平視場(chǎng)角360°，垂直視場(chǎng)角80°，2 000 m處精度50 mm。研究區(qū)范圍內(nèi)最低點(diǎn)云密度為4.53 point/m2。使用I-Site Studio軟件對(duì)LiDAR掃描獲取的點(diǎn)云進(jìn)行坐標(biāo)位置糾正、邊緣拼接等預(yù)處理，再由ArcGIS 10.1軟件對(duì)點(diǎn)云生成不規(guī)則三角網(wǎng)(TIN)，使用線性插值得到1 m分辨率DEM柵格數(shù)據(jù)，通過(guò)對(duì)柵格數(shù)據(jù)進(jìn)行子流域提取獲得試驗(yàn)區(qū)完整地形。

圖2 Leica-HDS8800掃描測(cè)量系統(tǒng)掃描場(chǎng)景Fig.2 Terrain scanning using Leica-HDS880 scanning measurement system

2.2.2 SRTM及ASTER GDEM數(shù)據(jù)

SRTM數(shù)據(jù)于2000年由航天飛機(jī)使用干涉雷達(dá)測(cè)繪獲得，最新發(fā)布的4.1版本由國(guó)際熱帶農(nóng)業(yè)中心(CIAT)使用新型插值算法得到。本研究使用的1″和3″分辨率SRTM數(shù)據(jù)下載自美國(guó)地質(zhì)調(diào)查局網(wǎng)站(http://earthexplorer.usgs.gov/)。

ASTER GDEM數(shù)據(jù)由高分辨率衛(wèi)星成像設(shè)備獲取，覆蓋全球地面面積的99%，高于SRTM的80%。V2版本的ASTER GDEM數(shù)據(jù)下載自地理空間數(shù)據(jù)云(http://www.gscloud.cn/)，分辨率為1″。

2.3 研究方法

2.3.1 不同數(shù)據(jù)源DEM對(duì)比分析

由于擁有試驗(yàn)區(qū)大量實(shí)地觀測(cè)資料，使用與實(shí)際地形要素對(duì)比的方法評(píng)價(jià)DEM質(zhì)量。研究區(qū)中部凸起的山脊及兩側(cè)的山谷是當(dāng)?shù)刈蠲黠@的地形特征(見(jiàn)圖1)。在研究中，選取橫跨該山脊的上、下2道剖面(A-A′和B-B′，見(jiàn)圖1)進(jìn)行不同數(shù)據(jù)源與分辨率的DEM質(zhì)量分析，并比較這些DEM對(duì)該地形的還原程度。在2道剖面中，偏低位置的A-A′剖面經(jīng)過(guò)山脊末端，偏高位置的B-B′剖面經(jīng)過(guò)山脊中部。此處選取的DEM包括1，30，90 m共3種分辨率的LiDAR DEM數(shù)據(jù)，1″和3″共2種分辨率的SRTM數(shù)據(jù)，以及1″分辨率的ASTER GDEM數(shù)據(jù)。為與另2種數(shù)據(jù)準(zhǔn)確對(duì)比，對(duì)1″分辨率的ASTER DEM使用雙線性插值得到90 m分辨率DEM數(shù)據(jù)加入分析。

通過(guò)對(duì)比3種數(shù)據(jù)源DEM數(shù)據(jù)展示的2處地形剖面，論證使用LiDAR技術(shù)獲取DEM數(shù)據(jù)的可靠性，以及SRTM與ASTER DEM數(shù)據(jù)還原青藏高原山區(qū)地形的能力。

2.3.2 不同分辨率LiDAR地形特征分析

目前用來(lái)對(duì)不同分辨率DEM進(jìn)行評(píng)價(jià)的指標(biāo)主要為坡度、剖面曲率和平面曲率。坡度使用ArcGIS 10.1軟件提取，剖面曲率和平面曲率根據(jù)Schmidt等[17]提出的方法計(jì)算。使用局部方差均值定量描述地形信息的尺度效應(yīng)，方差均值越大表明地形的局部起伏變化越劇烈[18]。針對(duì)某個(gè)特定柵格，以其為中心構(gòu)建分析窗口，求取分析窗口內(nèi)各柵格對(duì)應(yīng)地形屬性的方差作為該中心柵格的局部方差。最后求得所有柵格對(duì)應(yīng)局部方差的平均值即為局部方差均值。具體計(jì)算公式為：

(2)

(3)

本文對(duì)使用LiDAR點(diǎn)云生成的1 m分辨率DEM進(jìn)行雙線性插值，生成2，5，10，15，20，25，30 m 7種分辨率的DEM，加上1 m分辨率DEM共得到8種不同分辨率基于LiDAR數(shù)據(jù)的DEM。提取這8種DEM的坡度、平面曲率、剖面曲率3種地形特征，判斷分辨率變化對(duì)地形特征的影響，并使用局部方差均值分析不同地形特征時(shí)的最優(yōu)分辨率。

3 結(jié)果與分析

3.1 不同數(shù)據(jù)源DEM對(duì)比分析

對(duì)7種不同數(shù)據(jù)源或分辨率DEM的高程和坡度信息進(jìn)行統(tǒng)計(jì)，統(tǒng)計(jì)結(jié)果見(jiàn)表1。根據(jù)統(tǒng)計(jì)數(shù)據(jù)，在30 m分辨率下，不同數(shù)據(jù)源提取的地形存在較大差異，使用SRTM和ASTER GDEM數(shù)據(jù)提取的地形高程較LiDAR地形抬升明顯。此外，對(duì)比高程和坡度統(tǒng)計(jì)指標(biāo)可知， SRTM數(shù)據(jù)的坡度最大值和最小值大于LiDAR數(shù)據(jù)，但二者平均坡度相似；ASTER GDEM地形較LiDAR地形平坦。由于研究區(qū)范圍小，90 m LiDAR，SRTM，ASTER GDEM數(shù)據(jù)分別只有65，57，63個(gè)柵格，區(qū)域邊緣粗糙，且柵格劃分位置不同，導(dǎo)致這3種低分辨率DEM漏失部分邊緣高程信息情況不一，高程體現(xiàn)的地形平坦化信息不足，但是對(duì)比坡度可以發(fā)現(xiàn)，90 m分辨率地形相比30 m分辨率存在較嚴(yán)重的平坦化現(xiàn)象。

表1 不同數(shù)據(jù)源DEM高程及坡度統(tǒng)計(jì)Table 1 Statistics of elevation and slope gradient of DEMs from diverse sources

圖3 A-A′和B-B′地形剖面DEM對(duì)比Fig.3 Comparison of different DEMs at two topographic profiles

以研究區(qū)中部山脊劃分出的2道山谷作為評(píng)價(jià)DEM地形的標(biāo)準(zhǔn)，對(duì)比分析如圖3所示。1 m分辨率的LiDAR地形與30 m分辨率LiDAR地形大致相似。在30 m分辨率下，SRTM地形在兩側(cè)邊緣處與LiDAR大致相似，但是SRTM地形的2個(gè)剖面都僅顯示一個(gè)山谷，其中A-A′剖面的谷底異常平緩，B-B′剖面在山谷右側(cè)存在一小塊不合理平緩區(qū)域，這種誤差極可能是將谷底數(shù)據(jù)缺失區(qū)插值成平面造成的。ASTER GDEM高程遠(yuǎn)高于另2種DEM，在A-A′剖面處山脊不明顯，存在過(guò)高估計(jì)山谷谷底的現(xiàn)象，這與SRTM相似。ASTER GDEM數(shù)據(jù)顯示B-B′剖面存在多個(gè)拐點(diǎn)，但是僅有一個(gè)明顯下凹的山谷，這類(lèi)誤差可能是受當(dāng)?shù)爻Ｄ甏嬖诘脑旗F影響[8]。3種數(shù)據(jù)源的90 m分辨率地形變化均與各自數(shù)據(jù)源的30 m分辨率數(shù)據(jù)相似，但都無(wú)法準(zhǔn)確體現(xiàn)山脊的存在。通過(guò)對(duì)比，90 m分辨率的地形數(shù)據(jù)均無(wú)法準(zhǔn)確反映地形特征，30 m分辨率的ASTER GDEM和SRTM精度同樣不能滿足地形研究的需要。精度達(dá)到30 m以內(nèi)的LiDAR數(shù)據(jù)能夠清楚地展現(xiàn)2道山谷的地形剖面，并且有效地還原研究區(qū)的實(shí)際地形。

圖4 坡度與DEM分辨率關(guān)系Fig.4 Relation between the resolution of DEM and slope gradient

3.2 不同分辨率LiDAR地形特征分析

使用8種不同水平分辨率LiDAR DEM提取坡度的對(duì)比結(jié)果如圖4所示。圖4(a)中箱型圖的5處節(jié)點(diǎn)自下而上表示DEM坡度的最小值、1/4值、中值，3/4值和最大值。DEM最大坡度由1 m分辨率時(shí)的85.51°降低至30 m分辨率時(shí)的41.46°，當(dāng)分辨率>10 m時(shí)最低坡度自0°開(kāi)始升高。頻率分布在25%～75%區(qū)間時(shí)(圖4(a)中箱型圖的箱體部分)，坡度范圍由1 m分辨率的24.14°～38.90°減小到30 m分辨率的27.20°～32.72°。隨著分辨率降低，極陡坡和極緩坡面積均減少，最大坡度降低迅速，最小坡度上升較慢，各點(diǎn)的坡度向中值靠攏，坡度的變化范圍減小。DEM分辨率越低，提取的山坡地形忽略的微地形信息越多，坡度越趨近一致。從圖4(b)可知,2 m分辨率時(shí)DEM的平均坡度最高，且略大于1 m和5 m分辨率的平均坡度，當(dāng)分辨率>5 m時(shí)平均坡度迅速減小，整體地形出現(xiàn)坦化。圖4(c)展示了坡度的局部方差均值與DEM分辨率的關(guān)系，局部方差均值隨著分辨率的降低呈指數(shù)型減小，因此在研究區(qū)，越高分辨率的DEM保留的地形坡度信息越多，當(dāng)分辨率<5 m時(shí)局部方差均值較大。根據(jù)以上信息可以判斷，分辨率<5 m的研究區(qū)DEM保留了較完善的坡度信息，更適合用于當(dāng)?shù)氐钠露妊芯俊?/p>

平面曲率與剖面曲率隨DEM分辨率的變化較坡度更為明顯，如圖5所示，隨著分辨率降低，曲率變化范圍迅速減小。平均剖面曲率為負(fù)值說(shuō)明該山坡整體為凹型山坡，但隨著DEM分辨率的降低平均剖面曲率趨近于0，剖面形狀下凹不再明顯。分辨率越低，平面曲率均值越小，1 m和2 m分辨率時(shí)均值為正值，分辨率>5 m時(shí)均值為負(fù)值，因此實(shí)驗(yàn)山坡地形在高分辨率時(shí)總體呈發(fā)散坡型，低分辨率時(shí)總體呈收斂坡型，可見(jiàn)DEM分辨率對(duì)山坡數(shù)字地形有極大影響。使用局部方差均值顯示1 m分辨率DEM包含的曲率信息遠(yuǎn)超其他分辨率。

圖5 地形曲率與DEM分辨率關(guān)系Fig.5 Relation between terrain curvature and resolution of DEM

8種分辨率LiDAR DEM坡度、平面曲率、剖面曲率的結(jié)果顯示，分辨率降低會(huì)造成青藏高原高山地區(qū)DEM數(shù)字地形起伏減小，保留的地形信息迅速減少。

4 結(jié) 論

本文對(duì)比了不同數(shù)據(jù)源DEM，并分析了不同分辨率下LiDAR DEM提取的地形特征。結(jié)論如下：

(1)90 m分辨率DEM數(shù)據(jù)在以本文研究區(qū)為代表的青藏高原部分高山區(qū)存在漏失地形要素的現(xiàn)象，不能滿足地形研究的需要。在30 m分辨率DEM數(shù)據(jù)中，SRTM和ASTER GDEM等開(kāi)源DEM數(shù)據(jù)過(guò)高估計(jì)山谷谷底高程，漏失地形要素。使用這些數(shù)據(jù)得到的數(shù)字地形存在失真現(xiàn)象，高原地形研究需要高精度的地形獲取手段。

(2)青藏高原山坡陡峭，地形復(fù)雜，使用低分辨率DEM提取地表形態(tài)會(huì)有大量地形信息損失，并改變山坡的整體坡型。坡度相關(guān)的研究最好使用分辨率不超過(guò)5 m的DEM，曲率研究對(duì)DEM分辨率的要求更高。

(3)LiDAR技術(shù)通過(guò)激光掃描能準(zhǔn)確獲取地形特征，提供保留充足地形信息的高分辨率DEM，避免低分辨率DEM數(shù)據(jù)易出現(xiàn)的地形平坦化現(xiàn)象，在高精度地形數(shù)據(jù)獲取方面具有遠(yuǎn)大的應(yīng)用前景。