表面参数的工作原理

表面参数工具确定栅格表面的参数,例如坡向、坡度和曲率。

坡向

坡向表面参数用于确定下坡坡度所朝的方向。 输出栅格中各像元的值可指示出各像元位置处表面所朝向的罗盘方向。 将按照顺时针方向进行测量,角度范围介于 0(正北)到 360(仍是正北)之间,即完整的圆。 不具有下坡方向的平坦区域将赋值为 -1。

坡向

下图显示的是输入高程数据集和输出坡向栅格。

坡向输出示例

坡向应用

通过坡向表面参数类型,可执行以下操作:

  • 在寻找最适合滑雪的山坡的过程中,查找某座山所有朝北的坡。
  • 在统计各地区生物多样性的研究中,计算某区域中各个位置的日照强度。
  • 作为判断最先遭受洪流袭击的居住区位置研究的一部分,在某山区中查找所有朝南的山坡,从而判断出雪最先融化的位置。

测地线坡向计算

某一位置的测地线坡向是相对于北向在与椭球体表面相切的平面(下图中的蓝色平面)上测量的下坡表面的角度方向 α。

测地线组成部分说明

要计算每个位置处的坡向,需要使用最小二乘法 (LSM) 将二次或四次平面拟合到邻域像元。 将基于该表面计算像元位置处的表面法线。 在相同位置,还将计算与椭球体表面的切平面垂直的椭球体法线。

由于将椭球体表面的切平面视为参考平面,因此表面法线将投影到平面上。 最后,通过按顺时针方向测量北向与平面法线的投影之间的角度 α 来计算测地线坡向(请参见上图)。

斜率

坡度表面参数指栅格表面的每个像元处的陡度。 坡度值越小,地势越平坦;坡度值越大,地势越陡峭。

输出坡度栅格可使用两种单位进行计算:度和百分比(高程增量百分比)。 将高程增量百分比视为高程增量除以水平增量后再乘以 100,这样可以更好地理解高程增量百分比。 请参考下面的三角形 B。 当角度为 45 度时,高程增量等于水平增量,所以高程增量百分比为 100%。 如三角形 C 所示,当坡度角接近直角(90 度)时,高程增量百分比开始接近无穷大。

坡度和百分比
比较以度为单位和以百分比为单位的坡度值。

坡度最常用在高程数据集处理中,如下图所示。 较陡的坡度在输出坡度栅格上以暗棕色阴影显示。

坡度输出示例

该工具可与其他类型的连续数据(如人口)配合使用,用来识别值的急剧变化。

测地线坡度计算

测地线坡度是指地形面与椭球体表面之间形成的角度。 对于任何平行于椭球体表面的表面,其坡度为 0。 要计算每个位置处的坡度,需要使用最小二乘法 (LSM) 将二次或四次平面拟合到邻域像元。 将基于该表面计算像元位置处的表面法线。 在相同位置,还将计算与椭球体表面的切平面垂直的椭球体法线。 可以通过椭球体法线与地形面法线之间的角度来计算坡度(以度为单位)。 此角度等于地形面与椭球体表面之间的夹角。

表面曲率概述

曲率是表面参数类型的集合,用于描述表面的形状,通常沿表面上与通过该表面的平面相交形成的线。 从概念上讲,几何曲率会找到最佳拟合圆(密切圆)以近似曲线在任何点的形状。 曲率是该圆的半径的倒数 (1/r)。 线越直,与较大圆的拟合度越高,生成的曲率越小,而线越弯曲,与较小圆的拟合度越高,生成的曲率越大(由 Crane 于 2018 年提出)。

曲率为相反相切圆

剖面(法向坡度线)曲率

剖面(法向坡度线)曲率表面参数用于沿坡度线测量几何法向曲率。 有时称为剖面曲率,可以将其可视化为通过表面的垂直(剖面)横截面的形状。 如下图所示,垂直平面沿橙色线切入表面,如果切掉,将显示为表面的横截面剖面。

剖面(法向坡度线)曲率平面

该平面由两个矢量定义,黄色箭头表示梯度方向或坡度线箭头,红色箭头表示表面法向。 将通过组合这些红色和黄色矢量来定义橙色平面及其与表面的橙色相交线。 剖面曲率沿橙色平面中的橙色线(法向坡度线)计算。

在此处将采用由 Minár 等人于 2020 年提出的“法向坡度线”术语,以最大程度地减少歧义以及与早期术语的混淆。

通常应用此曲率来表征流体在重力作用下在表面上的向下加速和减速。 速度越高,水流可以携带和移动材料数量越大,加速区域变成腐蚀区域,减速区域变成沉积区域。

在下图中,椎体山脊周围具有高凸出剖面(法向坡度线)曲率的区域显示为紫色。 椎体底部具有高凹入剖面(法向坡度线)曲率的区域显示为橙色。 具有小曲率值的区域是透明的。

火山锥的剖面曲率
将显示基于 5 米分辨率 DSM 和 35 米邻域距离(15 x 15 像元窗口)、自适应邻域和二次表面拟合计算的剖面(法向坡度线)曲率。

此曲率的结果不同于先前曲率工具的剖面曲率输出。 稍后将提供有关剖面曲率和剖面(坡度线)曲率之间差异的说明。

用于计算剖面(法向坡度线)曲率的公式如下:

剖面(法向坡度线)曲率方程

  • 其中:

    KP = 剖面(法向坡度线)曲率

    z = f(x,y)

切向(法向等值线)曲率

切向(法向等值线)曲率表面参数用于测量垂直于坡度线且与等值线相切的几何法向曲率。 将其称为切向曲率,因为它用于测量与等值线相切的曲率。 将其描述为法向等值线 (Minár et al., 2020) 因为形成计算曲率时随沿的紫色线的紫色切割平面由蓝色等值线矢量和红色表面法向矢量定义。

切向(法向等值线)曲率平面

通常应用切向(法向等值线)曲率来表征流经表面的流体的地形汇聚和分散。

在下图中,椎体山脊和朝向您的山脊周围具有高凸出剖面(法向等值线)曲率的区域显示为蓝色。 这些是分散流的区域。 椎体内部具有高凹入切向(法向等值线)曲率的区域以红色显示汇聚流。 具有小曲率值的区域是透明的。

火山锥的切向曲率
将显示基于 5 米分辨率 DSM 和 35 米邻域距离(15 x 15 像元窗口)、自适应邻域和二次表面拟合计算的切向(法向等值线)曲率。

用于计算切向(法向等值线)曲率的公式如下:

切向(标准等高线)曲率方程

  • 其中:

    KT = 切向(法向等值线)曲率

    z = f(x,y)

平面(投影等值线)曲率

平面(投影等值线)曲率表面参数测量沿等值线的曲率。 其有时也称为等值线曲率或水平曲率。 投影等值线曲率沿蓝色等值线测量,其中水平平面与表面相交。

平面(投影等值线)曲率

用于计算平面(投影等值线)曲率的公式如下:

平面(投影等值线)曲率方程

  • 其中:

    KPC = 平面(投影等值线)曲率

    z = f(x,y)

下图显示了沿紫色线测量的切向(法向等值线)曲率与沿蓝色等值线测量的平面(投影等值线)曲率之间的差异。

切向和平面(投影等值线)曲率平面

等值线测地挠率

等值线测地挠率表面参数测量沿等值线的坡度角变化率。

用于计算等值线测地挠率的公式如下:

等值线测地挠率方程

  • 其中:

    τ = 等值线测地挠率

    z = f(x,y)

平均曲率

平均曲率表面参数用于测量表面的总曲率。 计算最小曲率和最大曲率的平均值即可获得平均曲率。 在数学上,相当于剖面(法向坡度线)和切向(法向等值线)曲率的平均值。 下图显示了剖面(法向坡度线)切割平面(橙色)和切向(法向等值线)切割平面(紫色)。

剖面和切向曲率平面

剖面(法向坡度线)和切向(法向等值线)曲率分别用于测量特定方向上的凸度和凹度;而平均曲率用于描述表面的固有凸度或凹度,与方向或重力影响无关。 其符号(正或负)无法明确指示是凸度还是凹度(极值除外),因为表面可以在一个方向上凹入,而在另一个方向上凸出。 高正值表示最大剥蚀区域,高负值表示最大累积区域 (Minár et al., 2020)。

用于计算平均曲率的公式如下:

平均曲率方程

  • 其中:

    KM = 平均曲率

    z = f(x,y)

高斯曲率

高斯曲率表面参数用于测量表面的总曲率。 计算最小曲率和最大曲率的结果即可获得高斯曲率,并且可能为负值和正值。 正值表示该像元处的表面为凸面,负值表示其为凹面。 值为 0 说明表面是平的。

用于计算高斯曲率的公式如下:

高斯曲率方程

  • 其中:

    KG = 高斯曲率

    z = f(x,y)

Casorati 曲率

Casorati 曲率表面参数用于测量表面的总曲率。 该值可为零或任意其他正数。 高正值表示在多个方向上的尖锐弯曲区域。

用于计算 Casorati 曲率的公式如下:

Casorati 曲率方程

  • 其中:

    KC = Casorati 曲率

    z = f(x,y)

基本曲率和组合曲率类型

切向(法向等值线)曲率、剖面(法向坡度线)曲率和等值线测地挠率被视为基本曲率类型,因为其他曲率可以表示为它们的组合。 根据 Minár 等人 (2020) 提出的术语,其被描述为基本三重奏。

除了上面为平均曲率、高斯曲率和 Casorati 曲率提供的表达式之外,这些曲率也可以作为基本三重奏的组合进行计算。

用于计算平均曲率的公式如下:

平均曲率组合方程

  • 其中:

    KM = 平均曲率

    KT = 切向(法向等值线)曲率

    KP = 剖面(法向坡度线)曲率

    z = f(x,y)

用于计算高斯曲率的公式如下:

高斯曲率组合方程

  • 其中:

    KG = 高斯曲率

    KT = 切向(法向等值线)曲率

    KP = 剖面(法向坡度线)曲率

    τ = 等值线测地挠率

    z = f(x,y)

用于计算 Casorati 曲率的公式如下:

Casorati 曲率组合方程

  • 其中:

    KC = Casorati 曲率

    KM = 平均曲率

    KG = 高斯曲率

    z = f(x,y)

与传统曲率工具算法的比较

表面参数工具使用与曲率工具不同的曲率算法及其计算中的测地线数学,因此不应直接比较这两个工具的输出。 表面参数曲率类型(剖面(法向坡度线)和切线(法向等值线))是真实的几何曲率 (Minár et al. 2020)。 表面参数工具的平均曲率是该点最小和最大曲率的平均值。 曲率工具的剖面和平面类型是方向导数,实际上并不能测量某个位置表面的几何曲率(Zeverbergen 和 Thorne 于 1987 年提出)。 表面参数工具中剖面(法向坡度线)的符号(正或负)与曲率工具中的剖面曲率相反。 表面参数工具会在测地线空间中进行计算,而曲率工具将使用平面坐标和数学进行计算。 表面参数工具可以适合二次曲面或四次曲面,曲率工具仅支持四次运算。

邻域距离

邻域距离值是从当前处理像元中心到正交邻域中心的地图距离。 较小邻域距离可捕获地表中的更多局部变化,从而获得较小地表要素的特征。 高程数据的分辨率越高,邻域距离越大可能越合适,因为数据中的微尺度误差(噪声)不能反映出感兴趣地貌过程,或者因为感兴趣地貌在更大距离上更容易识别。

下方示例使用 5 米分辨率数字表面模型 (DSM),在剖面(法向坡度线)曲率结果中显示出明显的噪点和条带伪影。 第一个图像使用默认的 3 x 3 窗口或 5 米邻域距离,第二个图像为 9 x 9 像元窗口或 20 米邻域距离,第三个图像使用 15 x 15 像元窗口或 35 米邻域距离。 在本例中,随着邻域距离的增加,景观中最重要或最主要的要素变得更加清晰,而噪点和条带伪影则更不明显。 虽然较大的邻域距离总是会导致较少的噪点,但最合适的距离将取决于数据的像元大小和应用的重要地貌要素的大小。

邻域距离的剖面示例
将显示邻域距离为 3 的火山锥的剖面(法向坡度线)斜率。 紫色区域为高凸出曲率,橙色区域为高凹入曲率。

最小邻域距离等于输入栅格的像元大小。 最大邻域距离等于像元大小的七倍,从而形成 15 x 15 像元窗口。 任何大于七倍像元大小的指定距离将始终使用 15 x 15 像元窗口。

邻域距离和移动窗口的像素数之间的关系
将显示邻域距离(橙色线)和移动窗口的像素数之间的关系。 像元大小为 10 米时;10 米的邻域距离将使用 3 x 3 像元窗口(此为默认设置),20 米的邻域距离将使用 5 x 5 像元窗口,而 30 米的邻域距离将使用 7 x 7 像元窗口。

如果指定的邻域距离不会导致像元大小的间隔为奇数,它将向上舍入到像元大小的下一个间隔。 例如,在上图中,如果指定的邻域距离为 25 米,它将向上舍入到像元大小的下一个间隔,即 30 米(像元大小的三倍),从而得到 7 x 7 像元窗口。

如果高程数据的空间分辨率比分析感兴趣地形所需的空间分辨率更加精细,则“邻域窗口”选项的替代选项是将数据重采样或聚合为更适合应用的更大像元大小。

表面参数计算对像元大小和邻域距离非常敏感。 Wilson (2018) 和 Minár 等人 (2020) 提供了关于该主题的大量研究的最新有效汇总。

自适应邻域

如果选中,使用自适应邻域参数会更改用于计算表面参数的邻域距离(窗口大小或面积),以更好地捕获地表中的相关变化。 该工具将根据邻域内所有像元的值计算与平均高程 (DEV) 的局部偏差(Wilson 和 Gallant 于 2000 年提出)来自动确定适当的窗口大小。 它尝试在使用最大的窗口大小的同时,最大程度地减少表面变化 (James et al., 2014)。 使用的最大窗口大小在邻域距离参数中指定。

当使用固定邻域计算表面参数时,将使用邻域内的所有像元值。 在使用自适应邻域计算表面参数时,将仅使用邻域的九个像元(外部正交和对角像元以及中心处理像元)。

说明了在自适应邻域计算中涉及哪个像元
点表示使用自适应邻域时使用 7 x 7 窗口计算表面参数时使用的像元中心。

当从高分辨率 DEM 分析具有各种大小的 terrain 要素的景观(例如具有小溪谷或河道的起伏群山)时,自适应邻域非常有用。 在这种情况下,对于小溪谷,可以使用较小的邻域距离(例如 1 米),而对于山丘则可以使用 10 或 15 米的较大邻域距离。

在下图中,较小的邻域适用于小溪和崖边,较大的邻域适用于从山丘到平原的过渡,更大的邻域适用于近似平坦的同构高原。

自适应邻域大小

邻域距离的边缘效应

当没有足够的信息可用于计算时,将向输出的外边周围的像元分配 NoData。

使用自适应邻域选项时,输出栅格的范围将在其外边周围缩减一个像元。

当使用大于输入像元大小的固定邻域距离时,将根据使用的邻域距离缩减输出栅格范围。 缩减量可以计算为(以像素为单位的窗口宽度 - 1)/ 2

例如,如果邻域距离导致使用 7 x 7 像元窗口,则输出栅格将在其外边周围缩减三个像元。

二次和四次

可以将两种类型的局部表面拟合到邻域窗口,即二次和四次。 对于大多数数据和应用,建议使用默认值二次。

二次表面是点的最小二乘拟合,并不完全通过所有点。 因此,使用二次表面可以最大程度地降低噪声表面数据(例如高分辨率激光雷达表面)的影响。 这将为所有表面参数生成更具代表性的结果,这在计算曲率时尤其重要。

当指定的邻域大小大于像元大小以及使用自适应邻域选项时,应使用二次曲面。

四次表面可精确拟合邻域像元中的数据。 此选项适用于高精度输入表面(无随机噪声)。 如果邻域距离大于输入栅格像元大小,则四次表面类型将失去精度优势,因此应将邻域距离保留为默认值(等于像元大小)。

测地线坐标转换

表面参数工具通过将地球形状视为椭球体,在 3D 地心坐标系(也称为地心地固 (ECEF) 坐标系)中执行计算。 计算结果不受数据集的投影方式影响。 如果在空间参考中定义了 z 单位,此方法将使用输入栅格的 z 单位。 如果输入的空间参考未定义 z 单位,则需要使用 z 单位参数进行定义。

ECEF 坐标系是以地心为原点的 3D 右手笛卡尔坐标系,其中所有位置均由 X、Y 和 Z 坐标表示。 有关以地心坐标系表示的目标位置 T 的示例,请参见下图:

测地线系统中的笛卡尔坐标说明

将表面栅格从输入坐标系转换到 3D 地心坐标系中。

测地线计算使用基于其测地坐标(纬度 φ、经度 λ、高度 h)计算的 X, Y, Z 坐标。 如果输入表面栅格的坐标系为投影坐标系 (PCS),则需要先将栅格重新投影到地理坐标系 (GCS),其中每个位置均具有测地坐标。 然后将其转换为 ECEF 坐标系。 高度 h(z 值)是以椭球体表面为参考的椭球体高度。 请参见下图。

椭球体高度

要从测地坐标(纬度 φ、经度 λ、高度 h)转换为 ECEF 坐标,请使用以下公式:

X = (N(φ) + h) * cos(φ) * cos(λ)

Y = (N(φ) + h) * cos(φ) * sin(λ)

Z = (b2 / a2 * N(φ) + h) * sin(φ)

  • 其中:

    N(φ) = a2 / √( a2 * cos(φ)2 + b2 * sin(φ)2)

    φ = 纬度

    λ = 经度

    h = 椭球体高度

    a = 椭球体的长轴

    b = 椭球体的短轴

在以上公式中,椭球体高度 h 以米为单位。 如果以任何其他单位指定输入栅格的 z 单位,则该单位将内部转换为米。

推荐阅读

要更深入地了解表面分析方法及其应用,请参考以下参考资料。 此外,Hengl 和 Reuter (2008) 以及 Wilson (2018) 均提供了上述以及更多地形分析技术及其应用的综合目录。 Minár 等人 (2020) 提供了关于陆地表面曲率的前期工作的全面汇总和比较,并明确定义了许多类型的曲率。

参考资料

B. Hofmann-Wellenhof, H. Lichtenegger and J. Collins, 2001. GPS - theory and practice. Section 10.2.1. p. 282.

Burrough, P. A., and McDonell, R. A., 1998. Principles of Geographical Information Systems (Oxford University Press, New York), 190 pp.

Crane K., 2018. Discrete Differential Geometry: An Applied Introduction. Notices of the AMS, Communication. https://www.cs.cmu.edu/~kmcrane/Projects/DDG/paper.pdf

David Eberly 1999. Least Squares Fitting of Data (Geometric Tools, LLC), pp. 3.

E.J.Krakiwsky, and D.E.Wells, 1971. Coordinate Systems In Geodesy (GEODESY AND GEOMATICS ENGINEERING, UNB), LECTURE NOTES, No16, 1971, pp. 18-38

Hengl T. and Reuter H. 2008. Geomorphometry Concepts, Software, Applications. Elsevier.

James D.E., M.D. Tomer, S.A. Porter. 2014. Trans-scalar landform segmentation from high-resolution digital elevation models. 海报发表于:ESRI 年度用户大会;2014 年 7 月;加利福尼亚州圣地亚哥。

Lancaster, P. and Šalkauskas, K. Curve and Surface Fitting: An Introduction. London: Academic Press, 1986.

Marcin Ligas, and Piotr Banasik, 2011. Conversion between Cartesian and geodetic coordinates on a rotational ellipsoid by solving a system of nonlinear equations (GEODESY AND CARTOGRAPHY), Vol. 60, No 2, 2011, pp. 145-159

Minár, J., Evans, I. S., & Jenčo, M. (2020). A comprehensive system of definitions of land surface (topographic) curvatures, with implications for their application in geoscience modelling and prediction. Earth-Science Reviews, 103414. https://doi.org/10.1016/j.earscirev.2020.103414

Wilson J.P and Gallant, J.C. (Eds.) 2000. Terrain Analysis: Principles and Applications. John Wiley & Sons, Inc.

Wilson J.P 2018. Environmental Application of Digital Terrain Modeling. John-Blackwell, Inc.

Zevenbergen, L. W., and C. R. Thorne. 1987. Quantitative Analysis of Land Surface Topography. Earth Surface Processes and Landforms 12: 47–56.