引言

相比钻孔数据受经费、环境及钻孔技术等的限制,基于地形坡度的大区域场地分类方法得到大量应用。目前国内外研究多基于Wald等(2006)建立的坡度−VS30模型进行场地类别划分(史大成,2009Thompson等,2014吴效勇,2019),其中坡度来源主要是基于30″分辨率数字高程数据(Digital Elevation Model,DEM)的提取。然而,提取的坡度受DEM分辨率、算法及数据类型(网格点类型)等因素影响,且分辨率是最主要的影响参数。国内外大量研究表明(Walker等,1999Mukherjee等,2013Grohmann,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),计算式为:

$ {\rm{slope}}{=\tan}\left(\sqrt{({\rm{d}}z/{\rm{d}}x{)}^{2}{+(}{\rm{d}}z/{\rm{d}}x{)}^{2}}\right) $ (1)

式中,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计算如下:


图 2 计算高程梯度变化率的像元窗口 Fig. 2 3×3 unit schematic diagram of slope calculation point
$ {\rm{d}}z/{\rm{d}}x = \frac{{({z_3} + 2{z_6} + {z_9}) - ({z_{\text{1}}} + 2{z_{\text{4}}} + {z_{\text{7}}})}}{{8w}} $ (2)
$ {\rm{d}}z/{\rm{d}}y = \frac{{({z_1} + 2{z_2} + {z_3}) - ({z_7} + 2{z_8} + {z_9})}}{{8w}} $ (3)

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

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

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分辨率提取的坡度进行拟合,各函数模型如下:

线性函数模型:

$ {S_{30}} = {a_1} + {b_1}{S_{{i}}} $ (4)

多项式函数模型:

$ {S_{30}} = {a_2}{S_{{i}}} + {b_2}{S_{{i}}}^2 + {c_{\text{1}}}{S_{{i}}}^3 + {d_1} $ (5)

幂函数模型:

$ {S_{{\text{30}}}} = {a_{\text{3}}}{S_{{i}}}^{{b_3}} + {c_{\text{2}}} $ (6)

式中,Si表示分辨率为i(其中i=1″、3″、30″)时计算得到的百分比坡度(m/m);ajbjj=1、2、3)、c1c2d1分别表示不同函数模型的回归参数。为使拟合公式计算结果更合理,需考虑边界效应,即当1″、3″分辨率下坡度趋于0时,30″下坡度也应趋于0。因此,当Si=0时,S30=0,即a1d1c2=0。

4.2 拟合回归

基于上述函数模型,分别对不同地貌单元1″、3″分辨率下坡度与30″分辨率下坡度进行拟合,坡度转换关系参数如表12所示。由表12可知,平原回归关系参数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)
表 2 不同地貌单元回归模型拟合参数及评价参数(3″与30″分辨率下) Table 2 Fitting and evaluation parameters of regression models for different landforms (3 arc second and 30 arc second)
5 回归方程评价方法
5.1 拟合优度

拟合优度R2计算如下:

$ {R^2} = \frac{{SSR}}{{SST}} = 1 - \frac{{ESS}}{{TSS}} $ (7)

式中,SSR为回归平方和,SST为总平方和,ESS为误差平方和,TSS为总离差平方和。

表12可知,幂函数模型转换的拟合优度R2最小,即幂函数模型拟合效果最差;线性函数模型拟合优度R2略小于多项式函数模型;整体来说,多项式函数模型拟合结果最优。

5.2 坡度分级平均相对差值

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

$ {\varDelta _{{\rm{means}}}} = {{\left(\sum\limits_{k = 1}^n {\frac{{\left| {{P_{{{30, }}k}} - {P_{i,{{ }}k}}} \right|}}{{{P_{{{30, }}k}}}}} \right)} / n} $ (8)

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

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


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

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

$ K = \frac{{{R^{\text{2}}}}}{{{\varDelta _{{\rm{means}}}}}} $ (9)

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

根据表12中的比例参数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)本文给出的最佳坡度转换公式主要应用于工程抗震领域中的场地分类。

参考文献
李昕蕾, 2020. 基于多因素的区域场地分类方法研究. 哈尔滨: 中国地震局工程力学研究所.
Li X. L., 2020. Regional site classification method based on multiple factors. Harbin: Institute of Engineering Mechanics, China Earthquake Administration. (in Chinese)
刘飞, 范建容, 崔兆岩等, 2019. 基于DEM分形特征的坡度尺度变换模型[J]. 山地学报, 37(1): 129-136.
Liu F., Fan J. R., Cui Z. Y., et al., 2019. A model of re-scaling slope based on DEM fractal feature[J]. Mountain Research, 37(1): 129-136.
刘利峰, 2018. DEM分辨率对白马小流域坡度与坡向提取的影响[J]. 山西水土保持科技, (2): 10-12.
Liu L. F., 2018. DEM resolving capability influence on extract slope and aspect in Baima minor watershed[J]. Soil and Water Conservation Science and Technology in Shanxi, (2): 10-12. DOI:10.3969/j.issn.1008-0120.2018.02.005
刘晓, 赵荣, 梁勇等, 2017. 顾及地貌与DEM分辨率的坡度算法适应性研究[J]. 测绘科学, 42(3): 29-34.
Liu X., Zhao R., Liang Y., et al., 2017. Research of algorithms for deriving slope based on different geomorphic types and DEM resolution[J]. Science of Surveying and Mapping, 42(3): 29-34.
乔之轩, 2019. 基于DEM的多分辨率水文特征提取. 北京: 北京林业大学.
Qiao Z. X., 2019. Multi-resolution hydrological feature extraction based on DEM. Beijing: Beijing Forestry University. (in Chinese)
师动, 朱奇峰, 杨勤科等, 2020. DEM分辨率对坡度算法选择影响的分析[J]. 山地学报, 38(6): 935-944.
Shi D., Zhu Q. F., Yang Q. K., et al., 2020. DEM resolution influence on slope algorithm selection[J]. Mountain Research, 38(6): 935-944.
史大成, 2009. 基于GIS的场地分类新方法研究. 哈尔滨: 中国地震局工程力学研究所.
Shi D. C., 2009. Study on new methods of site classification based on GIS. Harbin: Institute of Engineering Mechanics, China Earthquake Administration. (in Chinese)
土祥, 王春梅, 庞国伟等, 2018. 黄土丘陵沟壑区坡度尺度效应空间分异分析[J]. 山地学报, 36(6): 964-972.
Tu X., Wang C. M., Pang G. W., et al., 2018. On the spatial differentiation of scale effect of slopes in Loess Hill and Gully Region, China[J]. Mountain Research, 36(6): 964-972.
王英, 龚家国, 贾仰文等, 2019. 基于不同分辨率DEM提取坡度值的转换关系研究[J]. 水利水电技术, 50(8): 45-51.
Wang Y., Gong J. G., Jia Y. W., et al., 2019. Study on conversion relationship of slope information extracted from different resolution DEM[J]. Water Resources and Hydropower Engineering, 50(8): 45-51.
吴效勇, 2019. 基于多因素的地震影响场研究. 北京: 中国地震局地震预测研究所.
Wu X. Y., 2019. Study on seismic influence field based on multi-factors. Beijing: Institute of Earthquake Forecasting, China Earthquake Administration. (in Chinese)
Grohmann C. H., 2015. Effects of spatial resolution on slope and aspect derivation for regional-scale analysis[J]. Computers & Geosciences, 77: 111-117.
Horn B. K. P., 1981. Hill shading and the reflectance map[J]. Proceedings of the IEEE, 69(1): 14-47. DOI:10.1109/PROC.1981.11918
Mukherjee S., Mukherjee S., Garg R. D., et al., 2013. Evaluation of topographic index in relation to terrain roughness and DEM grid spacing[J]. Journal of Earth System Science, 122(3): 869-886. DOI:10.1007/s12040-013-0292-0
Olaya V., 2009. Basic land-surface parameters[J]. Developments in Soil Science, 33: 141-169.
Thompson E. M., Wald D. J., Worden C. B., 2014. A Vs30 map for California with geologic and topographic constraints[J]. Bulletin of the Seismological Society of America, 104(5): 2313-2321. DOI:10.1785/0120130312
Thornley J. H. M., 1976. Mathematical models in plant physiology[J]. London:Academic Press: 86-110.
Wald D. J., Earle P. E., Lin K. W., et al., 2006. Challenges in rapid ground motion estimation for the prompt assessment of global urban earthquakes[J]. Bulletin of the Earthquake Research Institute, 81: 273-281.
Walker J. P., Willgoose G. R., 1999. On the effect of digital elevation model accuracy on hydrology and geomorphology[J]. Water Resources Research, 35(7): 2259-2268. DOI:10.1029/1999WR900034
Wang C. M., Yang Q. K., Guo W. L., et al., 2012. Influence of resolution on slope in areas with different topographic characteristics[J]. Computers & Geosciences, 41: 156-168.