DEM 分辨率对提取东北地区坡度的影响研究1

刘佳轩 高 瑞 师黎静

1）中国地震局工程力学研究所，地震工程与工程振动重点实验室，哈尔滨 154200 2）黑龙江科技大学，建筑工程学院，哈尔滨 154200

引言

相比钻孔数据受经费、环境及钻孔技术等的限制，基于地形坡度的大区域场地分类方法得到大量应用。目前国内外研究多基于Wald 等（2006）建立的坡度−VS30模型进行场地类别划分（史大成，2009；Thompson 等，2014；吴效勇，2019），其中坡度来源主要是基于30″分辨率数字高程数据（Digital Elevation Model，DEM）的提取。然而，提取的坡度受DEM 分辨率、算法及数据类型（网格点类型）等因素影响，且分辨率是最主要的影响参数。国内外大量研究表明（Walker 等，1999；Mukherjee 等，2013；Grohmann，2015；刘利峰，2018；师动等，2020），提取的坡度随着DEM 分辨率的降低而减小，即坡度均值与标准差明显降低，且不同地区分辨率的影响不同，这对建立坡度-剪切波速、坡度-覆盖层厚度等关系模型具有重要影响，因此应选择合适分辨率的DEM 进行数字地形分析（Thornley，1976）。李昕蕾（2020）通过对新疆地区不同分辨率DEM 坡度的提取，发现利用精度高或大区域小比例尺数据进行场地分类的效果较差，而30″分辨率效果较好；乔之轩（2019）研究发现，高分辨率DEM 提取的坡度具有更精细的局部变化，但有时会掩盖区域坡度的整体变化趋势和特征。受计算机能力与硬件显示等条件限制，高分辨率DEM往往不能得到良好应用，有时需采用与研究比例尺相适应的分辨率。部分研究者在提取不同分辨率下的坡度时考虑地貌条件，发现不同地貌下的坡度尺度具有明显差异性（Wang 等，2012；刘晓等，2017；土祥等，2018）。目前，对我国东北地区考虑地貌条件下的坡度研究相对较少，且不同研究区域可能搜集到不同分辨率DEM，需对提取的坡度进行统一，王英等（2019）、刘飞等（2019）以分形理论为基础，建立不同分辨率下DEM 坡度转换模型，分别以泾河区域和四川丘陵地区某流域为研究对象，研究成果不适用于我国东北地区。

基于上述问题，本文以美国地质勘探局（United States Geological Survey，USGS）发布的我国东北地区1″、3″、30″分辨率DEM 数据为基础，分别提取研究区域内公里网格点的百分比坡度，分析不同分辨率下坡度差异原因。区分平原、丘陵和山地后进行坡度-频率分布曲线统计，分析不同地貌单元下的坡度随分辨率变化特征。通过线性函数模型、多项式函数模型及幂函数模型，拟合不同分辨率坡度转换关系，并对比拟合优度R2和坡度分级平均相对差值Δmean，选取不同DEM 分辨率得到坡度最佳拟合公式。

1 研究区域及数据

1.1 研究区域

研究区域主要包括黑龙江省、吉林省、辽宁省及内蒙古自治区部分地区，地理坐标为116°E～135°E，41°N～53°N。东西向最大横距约1 441 km，南北向最大纵距约1 656 km，面积约1 450 000 km2。东北地区海拔高度为−275～2 692 m，整个地势大致分为3 环，外围由黑龙江、乌苏里江、鸭绿江等流域构成，整体地势较低；中部为由大兴安岭、小兴安岭及长白山构成的山地与丘陵，整体地势较高；内部为由松花江、嫩江等流域构成的东北平原，整体地势较低。

1.2 DEM 高程数据

研究数据来源于USGS 官网发布的DEM，原始数据采样间隔分别为地球等角坐标系的1″（约30 m）、3″（约90 m）和30″（约1 km），图1 所示为3 种分辨率下高程数据模型，不同分辨率下DEM 表现出的整体地势结构基本相同，高分辨率DEM 对研究区域的描述更详细，具有更精确的高程数据范围及更精细的地形刻画，随着分辨率的降低，区域内高程数据产生一定程度的过滤与平滑，使区域整体高程数据范围减小，但整体地势结构更突出。

图1 东北地区不同分辨率下高程数据模型图Fig. 1 Elevation maps of 1 arc second, 3 arc second, and 30 arc second resolution in the Northeast, China

2 坡度提取

2.1 提取方法

目前最常用的坡度提取方法为拟合曲面法（Olaya，2009），计算式为：

式中，slope 表示坡度，本文使用广泛应用于大区域场地分类的百分比坡度；dz/dx，dz/dy分别表示计算点处2 个正交水平方向高程梯度变化率，无量纲，考虑复杂地貌坡度影响时，dz/dx，dz/dy优选算法建议使用三阶反距离权算法，即Horn 算法（Horn，1981）。

图2 所示为计算高程梯度变化率dz/dx，dz/dy的3×3 像元窗口，Horn 算法中高程梯度变化率dz/dx，dz/dy计算如下：

式中，w为像元宽度，即分辨率。

以1″和30″分辨率DEM 坡度计算过程为例，说明不同分辨率下坡度提取差异。当分辨率为1″时，式（2）和式（3）是对图2 代表的约90 m×90 m 范围内高程数据点的计算，最终通过代表中间栅格数据的高程梯度变化率得到坡度。当分辨率为30″时，式（2）和式（3）是对图2 代表的约3 km×3 km 范围内高程数据点的计算，进而计算得到坡度。虽然1″、30″分辨率下计算中心点可以相同，但坡度体现计算面积平均特性，随着分辨率的增大，栅格计算面积变小，因此高程数据有所差异，导致不同分辨率下的坡度有时存在较大差距，且区域附近坡度变化越大，可能导致不同分辨率下坡度差异越大。

图2 计算高程梯度变化率的像元窗口Fig. 2 3×3 unit schematic diagram of slope calculation point

2.2 提取结果

图3 所示为通过式（1）及DEM 数据计算得到的不同分辨率下百分比坡度图，由图3 可知，高分辨率下的坡度图对地表地形的描述更精细，坡度特征表现层次更丰富；随着DEM 分辨率的降低，坡度值域逐渐减小，部分地区坡度趋于0，整体地形结构更突出，但缺少地形细节的表达。通过GIS 软件建立与不同分辨率下坡度图相同的渔网式公里网格点，并提取网格点对应的坡度。

图3 东北地区不同分辨率下坡度图Fig. 3 Slope maps of 1 arc second, 3 arc second, and 30 arc second resolution in the Northeast, China

图4 所示为东北地区不同DEM 分辨率下坡度统计曲线，由图4（a）可知，高分辨率下坡度累计频率曲线较低分辨率更平缓，即高分辨率下坡度整体相对偏高；对于坡度值域而言，高分辨率下坡度值域范围较低分辨率广，1″、3″、30″分辨率下坡度分别约为0.6、0.4、0.2 m/m。由图4（b）可知，随着分辨率的降低，坡度均值和标准差明显减小，且在高分辨率下减小速率较快，在低分辨率下减小速率较慢，标准差的降低说明了坡度值域整体向均值附近靠近。图4 基本表现出东北地区整体坡度随着DEM 分辨率的降低而减小的现象。

图4 不同DEM 分辨率下的坡度统计Fig. 4 Slope statistics of different DEM resolutions

3 不同地貌坡度统计分析

区分不同地貌条件，对研究区域网格采样点坡度进行统计分析，研究不同地貌下坡度随分辨率变化特征，为拟合公式按照不同地貌条件提供理论基础。

图5 所示为我国东北地区地貌分布，平原、丘陵和山地地貌占比分别约为33.3%、12.2%、43.1%，其他地貌占比约为11.4%。通常，平原地貌较平坦或起伏较小；丘陵地貌由多个低矮山丘组成，具有一定起伏；山地地貌由众多山脉区域组成，地面起伏较大。对于海拔高度而言，通常山地海拔较高，丘陵次之，平原较低。对于坡度而言，通常山地坡度较大、较陡，丘陵次之，平原坡度较小、较缓。

图5 地貌分布图Fig. 5 Landform distribution map

图6 所示为不同分辨率下不同地貌单元坡度-频率分布曲线。由图6 可知，平原在坡度较小时频率分布较大，山地在坡度较大时频率分布相对较大；随着坡度的增大，平原坡度-频率分布曲线下降速率较山地快，而丘陵坡度-频率分布曲线变化特征处于两者中间；随着分辨率的减小，平原坡度-频率分布曲线峰值位置迅速向低坡度段移动，峰值区域变窄，坡度减小至0.05 m/m 以下，而山地坡度-频率分布曲线峰值位置向低坡度移动速度相对较缓，峰值区域变窄程度相对较小，坡度减小至0.15 m/m 以下，丘陵坡度-频率分布曲线变化特征处于平原与山地之间，坡度减小至0.1 m/m 以下。各地貌单元坡度减小程度不同，因此需考虑地貌条件，建立不同DEM 分辨率的坡度转换关系。

图6 不同分辨率下不同地貌单元坡度-频率分布曲线Fig. 6 Slope-frequency distribution curves of different topography at different resolutions

4 转换关系拟合回归

4.1 拟合模型

本文采用常用函数模型对不同DEM 分辨率提取的坡度进行拟合，各函数模型如下：

线性函数模型：

式中，Si表示分辨率为i（其中i=1″、3″、30″）时计算得到的百分比坡度（m/m）；aj、bj（j=1、2、3）、c1、c2、d1分别表示不同函数模型的回归参数。为使拟合公式计算结果更合理，需考虑边界效应，即当1″、3″分辨率下坡度趋于0 时，30″下坡度也应趋于0。因此，当Si=0 时，S30=0，即a1、d1、c2=0。

4.2 拟合回归

基于上述函数模型，分别对不同地貌单元1″、3″分辨率下坡度与30″分辨率下坡度进行拟合，坡度转换关系参数如表1、2 所示。由表1、2 可知，平原回归关系参数b最小，而山地回归关系参数b最大，丘陵回归关系参数b处于二者之间，基本符合随着分辨率减小，平原坡度减小程度较大，而山地坡度减小程度相对较小的变化规律；多项式函数及幂函数模型拟合参数具有类似规律。

表1 不同地貌单元回归模型拟合参数及评价参数（1″与30″分辨率下）Table 1 Fitting and evaluation parameters of regression models for different landforms (1 arc second and 30 arc second)

5 回归方程评价方法

5.1 拟合优度

表2 不同地貌单元回归模型拟合参数及评价参数（3″与30″分辨率下）Table 2 Fitting and evaluation parameters of regression models for different landforms (3 arc second and 30 arc second)

拟合优度R2计算如下：式中，SSR为回归平方和，SST为总平方和，ESS为误差平方和，TSS为总离差平方和。

由表1、2 可知，幂函数模型转换的拟合优度R2最小，即幂函数模型拟合效果最差；线性函数模型拟合优度R2略小于多项式函数模型；整体来说，多项式函数模型拟合结果最优。

5.2 坡度分级平均相对差值

高分辨率DEM 可获得高分辨率坡度，进而获得高分辨率等效剪切波速或覆盖层厚度，但场地分类主要使用30″分辨率DEM 提取的坡度，分类方法基本为建立坡度与等效剪切波速的分级对应关系，因此还需考虑转换前后30″坡度分级结果。将坡度按0.01 m/m 宽度范围进行分级统计，以30″分辨率DEM 提取的坡度分级作为参考值，分别计算不同拟合模型对高分辨率坡度转换后的坡度分级与参考值的平均相对差值Δmean，平均相对差值越接近于0，说明该拟合公式向30″坡度转化越好，平均相对差值越接近于1，说明该拟合公式向30″坡度转化越差，平均相对差值计算如下：

式中，Pi,k表示在i（i=1″、3″）分辨率下转换后第k组分段坡度，n表示坡度分级级数。

不同地貌单元不同分辨率下坡度转换结果如图7 所示，由图7 可知，平原1″、3″与30″真值接近，基本将分布在0～0.2 m/m 的坡度转换为0～0.03 m/m；丘陵1″向30″转换时，线性函数模型转换与多项式函数模型转换结果较好，而3″向30″转换时多项式函数模型转换结果较好；山地线性函数模型转换效果较好，而多项式函数模型与幂函数模型转换效果较差。由表1、2 可知，线性函数模型分级平均相对差值Δmean最小，即仅考虑坡度分级时，线性函数拟合为最佳拟合。

图7 不同地貌单元不同分辨率下坡度转换结果Fig. 7 Conversion results of DEM slopes with different resolutions under different landforms

5.3 最佳拟合模型

由于分别考虑拟合优度R2和分级平均相对差值Δmean，各模型在坡度转换上均有优劣，因此，以比例参数K综合考虑R2、Δmeans，计算公式如下：

K值越大，说明拟合模型拟合效果越好。

根据表1、2 中的比例参数K，最终确定不同地貌单元不同分辨率下坡度转换公式，如表3 所示。平原和山地最佳坡度拟合模型均为线性函数模型；丘陵1″与30″分辨率最佳坡度拟合模型为线性函数模型，3″与30″分辨率最佳坡度拟合模型为多项式函数模型。

表3 不同地貌单元不同分辨率下坡度转换公式Table 3 The slope conversion formula of different resolution DEM for different landforms

6 结论

本文以我国东北地区1″、3″、30″分辨率DEM 数据为基础，分别提取研究区域公里网格点百分比坡度，分析不同分辨率下坡度提取差异原因。区分平原、丘陵和山地进行坡度-频率分布曲线统计，分析不同地貌单元坡度变化特征。通过线性函数模型、多项式函数模型及幂函数模型拟合不同分辨率坡度转换关系，并计算拟合公式拟合优度R2及坡度分级平均相对差值Δmean，得出以下结论：

（1）不同分辨率提取的坡度不同，主要原因是坡度计算方法体现计算面积平均特性，随着分辨率的增大，栅格计算面积变小，高程数据有所差异，导致不同分辨率下坡度存在较大差异。由于坡度计算平均特性，区域附近坡度变化越大，可能导致不同分辨率下获得的坡度差异越大。

（2）我国东北地区坡度随DEM 分辨率的降低而减小，高分辨率下坡度图对地表地形的描述更精细，坡度特征表现层次更丰富，而低分辨率下整体地形结构更突出，但缺少地形细节的表达。随着分辨率的降低，我国东北地区坡度均值、方差及变化范围减小，坡度-频率分布曲线向低坡度段移动。

（3）随着分辨率的降低，不同地貌单元坡度-频率分布曲线变化趋势相同，变化程度不同。各地貌单元坡度-频率分布曲线均向低坡度方向移动，且峰值区域变窄，其中平原地貌变化趋势性显著，丘陵次之，山地最小。因此，建立坡度转换关系时需区分不同地貌条件。

（4）给出不同地貌单元下东北地区不同分辨率坡度转换模型，为利用30″分辨率的坡度进行大区域场地分类提供转换关系。平原和山地最佳坡度拟合模型均为线性函数模型；丘陵1″与30″分辨率最佳坡度拟合模型为线性函数模型，3″与30″分辨率最佳坡度拟合模型为多项式函数模型。

（5）本文给出的最佳坡度转换公式主要应用于工程抗震领域中的场地分类。