基于LiDAR的青藏高原山坡小流域地形提取方法

4

(1.河海大学 水文水资源学院，南京 210098； 2.水利部海河水利委员会 水文局，天津 300170； 3.南京水利科学研究院水工水力学研究所，南京 210029； 4.河海大学 水文水资源与水利工程科学国家重点实验室，南京 210098； 5.中国科学院南京土壤研究所 土壤物理与盐渍土研究室，南京 210008)

1 研究背景

数字高程模型(Digital Elevation Model，DEM)使用有限的数据点近似模拟地表形态，是地貌、水文、土壤、生态等学科获取地形信息的数据基础[1]。在实际应用中，DEM的水平分辨率会对地形特征分析的结果产生影响，低分辨率DEM存在较多的地形信息遗失，使描述对象失真[2-4]。为保证成果精度，土壤制图[5]、水土流失评价[6-7]等需要精确地表输移路径的领域对DEM分辨率要求极高。

青藏高原是目前水文研究的热点地区。受地表地形起伏、植被覆盖、云层等的影响[8]，目前广泛使用的SRTM和ASTER GDEM这些DEM数据在该地区存在较大误差[9-11]。此外，SRTM和ASTER GDEM数据的空间分辨率仅为1″(约30 m)，这样的分辨率难以满足小流域水文与地形研究的需要[12]。因此，青藏高原区的水文研究需要高精度DEM数据的获取方式。然而，DEM分辨率的提升虽然可以提高所代表地形信息精度，但同样会增加相关水文模型的计算时间，所以水文研究往往需要选取能适当保障地形信息精度的最优分辨率[12]。

激光探测与测量系统(Light Detection And Ranging，LiDAR)作为一种先进的三维空间信息获取手段，是目前获取高精度、高分辨率DEM的主要方法[13-16]。本文选取拉萨河夺底沟小流域一处山坡为研究对象，通过LiDAR点云生成DEM，使用地形剖面法与SRTM数据及ASTER GDEM数据对比，评价LiDAR DEM质量。此外，通过比较不同分辨率LiDAR DEM提取的坡度、平面曲率、剖面曲率，研究高山区地形属性与分辨率的关系，选取最适合水文研究的DEM分辨率。

2 研究资料与方法

2.1 研究区概况

研究区选自拉萨市区以北的夺底乡境内，位于91°11′07″E—91°11′47″E，29°46′19″N—29°46′43″N之间，面积约0.46 km2，属拉萨河流域夺底沟子流域的一级子流域。该区域地形起伏变化明显，东北高，西南低，内有一处山脊将其划分为2道沟谷，海拔分布见图1。研究区土壤类型为砂土，地表覆盖为高山草甸和裸岩，近年无大规模人类活动或地质作用等改变地表形态的事件发生。由于无高株植被和人工建筑物等的遮挡且地形稳定，适合采集LiDAR点云进行地形分析，并与近20 a来获取的遥感地形数据进行对比研究。

图1 研究山坡位置及现场Fig.1 Map and scenic picture of the hillside slope

2.2 地形数据

2.2.1 LiDAR地形数据

LiDAR利用光速恒定原理获取目标物距离，即使用传感器发射激光束至目标物表面，激光束经反射后由传感器接收，系统精确记录发射与接收激光束的时刻，使用式(1)计算得到激光束到达的目标点与传感器间的斜距R。

(1)

式中：c为光速；t为接收激光束时刻与发射激光束时刻的时间差。使用测得的斜距R结合扫描仪坐标、传感器高度以及扫描角度信息可以计算得到激光束到达的每一个地面点的x，y，z三维空间坐标，大量离散分布的三维点组成地形点云。DEM数据可通过地形处理软件对点云数据进行插值处理得到。

在本研究中，地形测量选用Leica-HDS8800扫描测量系统，如图2所示。该系统最远扫描距离可达2 000 m，扫描水平视场角360°，垂直视场角80°，2 000 m处精度50 mm。研究区范围内最低点云密度为4.53 point/m2。使用I-Site Studio软件对LiDAR扫描获取的点云进行坐标位置纠正、边缘拼接等预处理，再由ArcGIS 10.1软件对点云生成不规则三角网(TIN)，使用线性插值得到1 m分辨率DEM栅格数据，通过对栅格数据进行子流域提取获得试验区完整地形。

图2 Leica-HDS8800扫描测量系统扫描场景Fig.2 Terrain scanning using Leica-HDS880 scanning measurement system

2.2.2 SRTM及ASTER GDEM数据

SRTM数据于2000年由航天飞机使用干涉雷达测绘获得，最新发布的4.1版本由国际热带农业中心(CIAT)使用新型插值算法得到。本研究使用的1″和3″分辨率SRTM数据下载自美国地质调查局网站(http://earthexplorer.usgs.gov/)。

ASTER GDEM数据由高分辨率卫星成像设备获取，覆盖全球地面面积的99%，高于SRTM的80%。V2版本的ASTER GDEM数据下载自地理空间数据云(http://www.gscloud.cn/)，分辨率为1″。

2.3 研究方法

2.3.1 不同数据源DEM对比分析

由于拥有试验区大量实地观测资料，使用与实际地形要素对比的方法评价DEM质量。研究区中部凸起的山脊及两侧的山谷是当地最明显的地形特征(见图1)。在研究中，选取横跨该山脊的上、下2道剖面(A-A′和B-B′，见图1)进行不同数据源与分辨率的DEM质量分析，并比较这些DEM对该地形的还原程度。在2道剖面中，偏低位置的A-A′剖面经过山脊末端，偏高位置的B-B′剖面经过山脊中部。此处选取的DEM包括1，30，90 m共3种分辨率的LiDAR DEM数据，1″和3″共2种分辨率的SRTM数据，以及1″分辨率的ASTER GDEM数据。为与另2种数据准确对比，对1″分辨率的ASTER DEM使用双线性插值得到90 m分辨率DEM数据加入分析。

通过对比3种数据源DEM数据展示的2处地形剖面，论证使用LiDAR技术获取DEM数据的可靠性，以及SRTM与ASTER DEM数据还原青藏高原山区地形的能力。

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

目前用来对不同分辨率DEM进行评价的指标主要为坡度、剖面曲率和平面曲率。坡度使用ArcGIS 10.1软件提取，剖面曲率和平面曲率根据Schmidt等[17]提出的方法计算。使用局部方差均值定量描述地形信息的尺度效应，方差均值越大表明地形的局部起伏变化越剧烈[18]。针对某个特定栅格，以其为中心构建分析窗口，求取分析窗口内各栅格对应地形属性的方差作为该中心栅格的局部方差。最后求得所有栅格对应局部方差的平均值即为局部方差均值。具体计算公式为：

(2)

(3)

本文对使用LiDAR点云生成的1 m分辨率DEM进行双线性插值，生成2，5，10，15，20，25，30 m 7种分辨率的DEM，加上1 m分辨率DEM共得到8种不同分辨率基于LiDAR数据的DEM。提取这8种DEM的坡度、平面曲率、剖面曲率3种地形特征，判断分辨率变化对地形特征的影响，并使用局部方差均值分析不同地形特征时的最优分辨率。

3 结果与分析

3.1 不同数据源DEM对比分析

对7种不同数据源或分辨率DEM的高程和坡度信息进行统计，统计结果见表1。根据统计数据，在30 m分辨率下，不同数据源提取的地形存在较大差异，使用SRTM和ASTER GDEM数据提取的地形高程较LiDAR地形抬升明显。此外，对比高程和坡度统计指标可知， SRTM数据的坡度最大值和最小值大于LiDAR数据，但二者平均坡度相似；ASTER GDEM地形较LiDAR地形平坦。由于研究区范围小，90 m LiDAR，SRTM，ASTER GDEM数据分别只有65，57，63个栅格，区域边缘粗糙，且栅格划分位置不同，导致这3种低分辨率DEM漏失部分边缘高程信息情况不一，高程体现的地形平坦化信息不足，但是对比坡度可以发现，90 m分辨率地形相比30 m分辨率存在较严重的平坦化现象。

表1 不同数据源DEM高程及坡度统计Table 1 Statistics of elevation and slope gradient of DEMs from diverse sources

图3 A-A′和B-B′地形剖面DEM对比Fig.3 Comparison of different DEMs at two topographic profiles

以研究区中部山脊划分出的2道山谷作为评价DEM地形的标准，对比分析如图3所示。1 m分辨率的LiDAR地形与30 m分辨率LiDAR地形大致相似。在30 m分辨率下，SRTM地形在两侧边缘处与LiDAR大致相似，但是SRTM地形的2个剖面都仅显示一个山谷，其中A-A′剖面的谷底异常平缓，B-B′剖面在山谷右侧存在一小块不合理平缓区域，这种误差极可能是将谷底数据缺失区插值成平面造成的。ASTER GDEM高程远高于另2种DEM，在A-A′剖面处山脊不明显，存在过高估计山谷谷底的现象，这与SRTM相似。ASTER GDEM数据显示B-B′剖面存在多个拐点，但是仅有一个明显下凹的山谷，这类误差可能是受当地常年存在的云雾影响[8]。3种数据源的90 m分辨率地形变化均与各自数据源的30 m分辨率数据相似，但都无法准确体现山脊的存在。通过对比，90 m分辨率的地形数据均无法准确反映地形特征，30 m分辨率的ASTER GDEM和SRTM精度同样不能满足地形研究的需要。精度达到30 m以内的LiDAR数据能够清楚地展现2道山谷的地形剖面，并且有效地还原研究区的实际地形。

图4 坡度与DEM分辨率关系Fig.4 Relation between the resolution of DEM and slope gradient

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

使用8种不同水平分辨率LiDAR DEM提取坡度的对比结果如图4所示。图4(a)中箱型图的5处节点自下而上表示DEM坡度的最小值、1/4值、中值，3/4值和最大值。DEM最大坡度由1 m分辨率时的85.51°降低至30 m分辨率时的41.46°，当分辨率>10 m时最低坡度自0°开始升高。频率分布在25%～75%区间时(图4(a)中箱型图的箱体部分)，坡度范围由1 m分辨率的24.14°～38.90°减小到30 m分辨率的27.20°～32.72°。随着分辨率降低，极陡坡和极缓坡面积均减少，最大坡度降低迅速，最小坡度上升较慢，各点的坡度向中值靠拢，坡度的变化范围减小。DEM分辨率越低，提取的山坡地形忽略的微地形信息越多，坡度越趋近一致。从图4(b)可知,2 m分辨率时DEM的平均坡度最高，且略大于1 m和5 m分辨率的平均坡度，当分辨率>5 m时平均坡度迅速减小，整体地形出现坦化。图4(c)展示了坡度的局部方差均值与DEM分辨率的关系，局部方差均值随着分辨率的降低呈指数型减小，因此在研究区，越高分辨率的DEM保留的地形坡度信息越多，当分辨率<5 m时局部方差均值较大。根据以上信息可以判断，分辨率<5 m的研究区DEM保留了较完善的坡度信息，更适合用于当地的坡度研究。

平面曲率与剖面曲率随DEM分辨率的变化较坡度更为明显，如图5所示，随着分辨率降低，曲率变化范围迅速减小。平均剖面曲率为负值说明该山坡整体为凹型山坡，但随着DEM分辨率的降低平均剖面曲率趋近于0，剖面形状下凹不再明显。分辨率越低，平面曲率均值越小，1 m和2 m分辨率时均值为正值，分辨率>5 m时均值为负值，因此实验山坡地形在高分辨率时总体呈发散坡型，低分辨率时总体呈收敛坡型，可见DEM分辨率对山坡数字地形有极大影响。使用局部方差均值显示1 m分辨率DEM包含的曲率信息远超其他分辨率。

图5 地形曲率与DEM分辨率关系Fig.5 Relation between terrain curvature and resolution of DEM

8种分辨率LiDAR DEM坡度、平面曲率、剖面曲率的结果显示，分辨率降低会造成青藏高原高山地区DEM数字地形起伏减小，保留的地形信息迅速减少。

4 结 论

本文对比了不同数据源DEM，并分析了不同分辨率下LiDAR DEM提取的地形特征。结论如下：

(1)90 m分辨率DEM数据在以本文研究区为代表的青藏高原部分高山区存在漏失地形要素的现象，不能满足地形研究的需要。在30 m分辨率DEM数据中，SRTM和ASTER GDEM等开源DEM数据过高估计山谷谷底高程，漏失地形要素。使用这些数据得到的数字地形存在失真现象，高原地形研究需要高精度的地形获取手段。

(2)青藏高原山坡陡峭，地形复杂，使用低分辨率DEM提取地表形态会有大量地形信息损失，并改变山坡的整体坡型。坡度相关的研究最好使用分辨率不超过5 m的DEM，曲率研究对DEM分辨率的要求更高。

(3)LiDAR技术通过激光扫描能准确获取地形特征，提供保留充足地形信息的高分辨率DEM，避免低分辨率DEM数据易出现的地形平坦化现象，在高精度地形数据获取方面具有远大的应用前景。