坡度的工作原理

需要 Spatial Analyst 许可。

需要 3D Analyst 许可。

坡度工具可确定栅格表面每个像元处的陡度。坡度值越小，地势越平坦；坡度值越大，地势越陡峭。

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

以度和百分比计的坡度值比较。

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

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

计算方法和边缘效应

有两种方法可用于坡度计算。可以通过方法参数在执行平面或测地线计算之间进行选择。

对于平面方法，坡度按一个像元到与其直接相邻的像元方向上值的最大变化率进行测量。利用 2D 笛卡尔坐标系对投影平面执行计算。坡度值通过最大平均值法（Burrough，1998）来计算。

使用测地线方法时，通过将地球形状视为椭球体，在 3D 笛卡尔坐标系中执行计算。通过测量地形面与参考基准面之间的角度来计算坡度值。

可以使用 3 像元 x 3 像元的邻域（移动窗口）来执行平面和测地线计算。对于每个邻域，如果处理（中心）像元为 NoData，则输出也为 NoData。计算还要求至少 7 个与处理像元相邻的像元具有有效值。如果有效像元少于 7 个，则不执行计算，并且该处理像元处的输出将为 NoData。

输出栅格最外侧行列的像元为 NoData。这是因为沿着输入数据集边界，这些像元没有足够的有效相邻像元。

平面方法

对于每个像元，该工具可计算从该像元到与其相邻的像元方向上值的最大变化率。实际上，高程随着像元与其相邻的八个像元之间距离的变化而产生的最大变化率可用来识别自该像元开始的最陡坡降。

平面坡度算法

坡度取决于表面从中心像元开始在水平 (dz/dx) 方向和垂直 (dz/dy) 方向上的变化率（增量）。用来计算坡度的基本算法如下所示：

 slope_radians = ATAN ( √ ([dz/dx]2 + [dz/dy]2) )

坡度通常以度为单位来测量，其算法如下所示：

 slope_degrees = ATAN ( √ ([dz/dx]2 + [dz/dy]2) ) * 57.29578

注：

这里显示的值 57.29578 是对 180/pi 的计算结果进行截断而得到的值。

坡度算法也可以表示为：

 slope_degrees = ATAN (rise_run) * 57.29578

• 其中：

   rise_run = √ ([dz/dx]2 + [dz/dy]2]

中心像元及其相邻的八个像元的值确定水平增量和垂直增量。这些相邻的像元使用字母 a 至 i 进行标识，其中 e 表示当前正在计算坡向的像元。

表面分析扫描窗口

像元 e 在 x 方向上的变化率将通过以下算法进行计算：

  [dz/dx] = ((c + 2f + i)*4/wght1 - (a + 2d + g)*4/wght2) / (8 * x_cellsize)

• 其中：

  wght1 和 wght2 是有效像元的水平加权计数。

  例如，如果：

  • c、f 和 i 均有一个有效值，则 wght1 = (1+2*1+1) = 4。
  • i 为 NoData，则 wght1 = (1+2*1+0) = 3。
  • f 为 NoData，则 wght1 = (1+2*0+1) = 2。

  除邻域位置为 a、d 和 g 以外，类似的逻辑也适用于 wght2。

像元 e 在 y 方向上的变化率使用以下算法进行计算：

  [dz/dy] = ((g + 2h + i)*4/wght3 - (a + 2b + c)*4/wght4) / (8 * y_cellsize)

• 其中：

  wght3 和 wght4 与 [dz/dx] 计算中的概念相同。

平面坡度计算示例

例如，将计算如下所示的移动窗口内中心像元的坡度值。

坡度示例输入

中心像元 e 在 x 方向上的变化率为：

  [dz/dx] = ((c + 2f + i)*4/wght1 - (a + 2d + g)*4/wght2) / (8 * x_cellsize)          = ((50 + 60 + 10)*4/(1+2+1) - (50 + 60 + 8)*4/(1+2+1)) / (8 * 5)          = (120 - 118) / 40          = 0.05

像元 e 在 y 方向上的变化率为：

  [dz/dy] = ((g + 2h + i)*4/wght3 - (a + 2b + c)*4/wght4) / (8 * y_cellsize)          = ((8 + 20 + 10)*4/(1+2+1) - (50 + 90 + 50)*4/(1+2+1)) / (8 * 5)          = (38 - 190 ) / 40          = -3.8

代入 x 方向和 y 方向上的变化率，使用下列算法来计算中心像元 e 的坡度：

  rise_run = √ ([dz/dx]2 + [dz/dy]2)           = √ ((0.05)2 + (-3.8)2)           = √ (0.0025 + 14.44)           = 3.80032

  slope_degrees = ATAN (rise_run) * 57.29578                = ATAN (3.80032) * 57.29578                = 1.31349 * 57.29578                = 75.25762

像元e 的整型坡度值为 75 度。

坡度示例输出

测地线方法

测地线方法通过将地球形状视为椭球体，在 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( φ ) = a²/ √(a²cosφ²+b²sinφ²)
  • φ = 纬度
  • λ = 经度
  • h = 椭球体高度
  • a = 椭球体的长轴
  • b = 椭球体的短轴

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

坡度计算

测地线坡度是指地形面与椭球体表面之间形成的角度。对于任何平行于椭球体表面的表面，其坡度为 0。要计算每个位置处的坡度，可以使用最小二乘法 (LSM) 在每个处理像元周围拟合 3 像元 x 3 像元的邻域平面。LSM 中的最佳拟合可使实际 z 值与拟合 z 值之间的平方差 (dz_i) 的总和最小。有关示例，请参见下图。

最小二乘拟合示例

此处，平面表示为 z = Ax + By + C。对于每个像元中心，dz_i 是实际 z 值与拟合 z 值之间的差值。

如果 ∑⁹_i=1dz_i² 最小，则平面将实现最佳拟合。

对平面进行拟合后，将在像元位置处计算表面法线。在相同位置，还将计算与椭球体表面的切平面垂直的椭球体法线。

测地线坡度计算

可以通过椭球体法线与地形面法线之间的角度（此处表示为 β）来计算坡度（以度为单位）。在上图中，角度 α 为测地线坡度，根据相等几何原理，该角度与角度 β 相等。

要以高程增量百分比计算坡度，可以使用以下公式：

Slope_PercentRise = ATAN(β) * 100%

是否应使用表面参数工具？

如果输入栅格参数值（Python 中的 in_raster）为高分辨率值（其像元大小小于几米，或者特别是存在噪点），请考虑使用表面参数工具及其用户定义的邻域距离选项，而非此工具的直接 3 x 3 邻域。 使用较大邻域可以最大程度地减少噪点表面的影响。当使用高分辨率表面时，使用较大邻域也可以更好地表示地形和表面特征。

GPU 的使用

对于测地线方法，如果您的系统上已安装某些 GPU 硬件，则此工具能够显著提高性能。有关如何支持、配置以及启用此功能的详细信息，请参阅通过 Spatial Analyst 处理 GPU 一节。

参考文献

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

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

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

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