CN112926408B - 一种基于数字高程模型的冰川崩解前沿自动提取方法 - Google Patents

Abstract

一种基于数字高程模型的冰川崩解前沿自动提取方法，首先根据山体阴影可视化技术创建数字高程模型的山体阴影地形图，再根据山体阴影地形图计算出视觉粗糙度；其次根据视觉粗糙度和数字高程模型中所有像素的对应高程值计算出高程地形特征图；再次根据基于博弈论的分类方法对高程地形特征图中所有像素进行海洋、冰盖两类分类，海洋、冰盖两类即I类与II类；最后根据分类结果提取边界线并矢量化。该方法通过创建高程地形特征图来有效地放大了冰盖和海洋两类之间的对比度，解决了在光学和合成孔径雷达遥感影像上类别混淆的难题；针对高程地形特征图的影像特点，设计了基于博弈论的两类分类方法，高效地区分冰盖和海洋两类地物。

Description

一种基于数字高程模型的冰川崩解前沿自动提取方法

技术领域

本发明属于遥感影像处理领域，特别涉及一种基于数字高程模型的冰川崩解前沿自动提取方法。

背景技术

冰川崩解前沿指的是冰盖和海洋的分界线，作为冰盖动力学的一个重要参数，它的变化揭示了冰川或冰架对不断变化的海洋环流、变暖的气候以及自身结构条件变化的动态响应，也影响了未来冰川的物质平衡变化。由于极地冰川特殊的地理位置和环境特点，人迹罕至，实地测量难以实现，因此对于大范围的极地冰川崩解前沿提取，主要依赖于遥感数据和遥感技术。目前提取冰川崩解前沿的遥感数据主要是基于光学和合成孔径雷达(SAR)影像，即利用冰盖和海洋不同地物类型的光谱特征或后向散射回波特点的不同，进行不同类别的区分。

基于光学或者SAR影像提取冰川崩解前沿的半自动以及自动方法主要包括两类，首先是基于边缘提取算子的冰川崩解前沿提取方法，常用的边缘提取算子有Sobel、Canny等。其次是基于影像分类的方法将数据分为冰盖和海洋两大类，其中冰盖一类包括注出冰川和冰架，而海洋一类则包括海水、海冰、冰山、混杂冰等地物类别，进而提取两类别之间的边界。基于边缘提取算子的方法通常无法自动提取完整的冰川崩解前沿，需要一定的后处理去除提取出的杂乱线特征；基于分类的方法面临的主要问题就是冰盖和海洋交接位置处不同类别地物特征在遥感影像上的对比度不足，比如在光学遥感影像上，海冰会存在与冰川或冰架非常相似的图谱特征；在SAR遥感影像上，冰川或冰架的后向散射回波强度一般来讲要高于海水，当冰川或冰架表面有一定程度的融水或者薄雪覆盖时，冰川或冰架表面的后向散射回波强度就会明显降低，另一方面来讲，当海水表面存在海冰，或者风浪时都会增大海水表面的粗糙度，进而增大其雷达回波强度，造成冰盖和海洋的困难区分。目前，最广泛采用的冰川崩解前沿提取方法仍然是手工提取。

发明内容

为了解决现有技术中存在的技术问题，本发明提供了一种基于数字高程模型的冰川崩解前沿自动提取方法，该方法通过创建高程地形特征图增强冰盖与海水之间的对比度，再利用一种创新的基于博弈论的分类方法对创建的高程地形特征图进行自动、高效分类，分为冰盖和海水两类。

实现本发明目的的技术方案为，一种基于数字高程模型的冰川崩解前沿自动提取方法，该方法至少包括如下步骤：

(1)、根据山体阴影可视化技术创建数字高程模型的山体阴影地形图，再根据山体阴影地形图计算出视觉粗糙度；

(2)、根据视觉粗糙度和数字高程模型中所有像素的对应高程值计算出高程地形特征图；

(3)、根据基于博弈论的分类方法对高程地形特征图中所有像素进行海洋、冰盖两类分类，海洋、冰盖两类即I类与II类；

(4)、根据分类结果提取边界线并矢量化。

对上述技术方案的进一步改进为：步骤(1)中所述山体阴影可视化技术采用双光源入射生成两副山体阴影地形图，且光源之间入射方向互为相反方位角。

且步骤(1)中所述视觉粗糙度R(n)通过寻找山体阴影地形图中n×n的窗口中像素值最大和最小的像素点B_max，B_min得到，具体为R(n)＝B_max-B_min。

且将所述相反方位角入射光源生成的两副山体阴影地形图的视觉粗糙度进行比较，取数值最大的视觉粗糙度进行计算高程地形特征图。

且其特征在于、步骤(2)中所述高程地形特征图通过如下步骤获得：

(2.1)、取数字高程模型中全部像素高程值的百分之三十分位数q₁与I类的临界高程值γ₁中的最小值α，取数字高程模型中全部像素高程值的百分之七十分位数q₂与II类的临界高程值γ₂中的最大值β，即α＝min{H_q1,γ₁},β＝max{H_q2,γ₂}；

(2.2)、根据视觉粗糙度R、数字高程模型中像素的高程值H、α和β计算出像素对应的高程地形特征值，具体通过如下公式获得：

式中，Median{H,α,β}为高程值H、α和β中的中位数；

(2.3)、重复步骤(2.2)直至所有的像素高程地形特征值被计算出；

(2.4)、根据高程地形特征值得到高程地形特征图，即将高程地形特征值中小于其他99％特征值的像素值赋值为0，高程地形特征值中大于其他99％特征值的像素值赋值为255，其他像素根据高程地形特征值的大小从小到大依次在(0,255)范围内赋值。

所述I类的临界高程值γ₁为30，所述II类的临界高程值γ₂为60。

且步骤(3)所述基于博弈论的分类方法为：将所有像素都被视为一个理性玩家，每个玩家的策略集为L＝{0，1，2}，其中0表示未分类的，1表示类I，2表示类II，对于每个策略，对应的策略成本为E_i(x_i)。考虑到N个玩家的所有策略，总策略成本E(x)是每个玩家的决策的成本值E_i(x_i)的总和，通过最小化总体成本E(x)，可以得到所有像素的最优策略集合；

步骤(3)所述高程地形特征图中所有像素通过如下步骤分类：

(3.1)、将像素的高程值H小于α的像素标记为I类；将像素的高程值H大于β的像素标记为II类；

(3.2)、将其他未分类的像素，即像素的高程值H大于α且小于β的像素按照使得其自身的策略成本E_i(x_i)最小进行标记，具体E_i(x_i)通过如下公式获得：

式中，i为任意像素，j为i的邻域内的像素，x_i为像素i选择的策略值，x_j为i的邻域内的像素j所选择的策略值，r为固定的惩罚值即为255，策略集合{0,1,2}代表未分类、I类、II类标记；f(i)为像素i的高程地形特征值，f(j)为像素j的高程地形特征值，N_i为像素i邻域内的其他像素，min||f(i)-f(j)||为像素i的高程地形特征值f(i)与像素i邻域像素j的高程地形特征值f(j)之间差值的最小值；

步骤(3.3)、重复步骤(3.2)直至所有的点被分类。

且所述边界线通过如下方式获得：将不同类别的相邻像素作为边界区域，再将边界区域转换为1像素宽的边界线，然后将1像素宽的边界线从栅格图像转换成矢量文件。

由上述技术方案可知，本发明提供的方法通过创建山体阴影地形图，再通过山体阴影地形图计算出视觉粗糙度，再根据视觉粗糙度和数字高程模型中所有像素的对应高程值计算出高程地形特征图，来增强冰盖和海水之间的对比度。

且基于博弈论的分类方法对高程地形特征图中所有像素进行海洋、冰盖两类分类，通过最小化总体成本，来增加区分冰盖和海水的效率。

且将第一类和第二类之间的边界提取并矢量化，可以直观的展示提取结果，方便工作人员的观看。

附图说明

图1是本发明的基于数字高程模型的冰川崩解前沿自动提取方法的12m分辨率TanDEM-X DEM数据生成山体阴影地形图。

图2是本发明的基于数字高程模型的冰川崩解前沿自动提取方法的步骤2中高程地形特征图。

图3是本发明的基于数字高程模型的冰川崩解前沿自动提取方法的步骤3中分类结果图。

图4是本发明的基于数字高程模型的冰川崩解前沿自动提取方法的步骤4中最终结果图。

具体实施方式

下面结合附图和实施例对本发明进行详细具体说明，本发明的内容不局限于以下实施例。

一种基于数字高程模型的冰川崩解前沿自动提取方法，该方法包括如下步骤：

(1)、参考图1，根据山体阴影可视化技术创建数字高程模型的山体阴影地形图，再根据山体阴影地形图计算出视觉粗糙度；

在本实施例中，山体阴影可视化技术假想一个光源在某个方向和某个太阳高度的模拟下，对数字高程模型中每个像素进行赋值，具体通过如下公式来生成一副0-255的灰度图。

H＝255.0×((cos(Z_rad)×cos(S_rad))+(sin(Z_rad)×sin(S_rad)×cos(A_{z_ra}d-A_{s_ ra}d)))

式中，H为计算出的山体阴影值，Z_rad为太阳天顶角的弧度数，S_rad为某一点的坡度弧度数，A_{z_rad}为太阳光线方向角的弧度数，A_{s_rad}为某一点的坡向弧度数，当H＜0时,则令H＝0。

视觉粗糙度R(n)通过寻找山体阴影地形图中n×n的窗口中像素值最大和最小的像素点B_max，B_min得到，具体为R(n)＝B_max-B_min。

因为单入射光源存在方向斜坡阴影会导致信息损失，所以在另一些实施例中会采用双光源入射生成两副山体阴影地形图，且光源之间入射方向互为相反方位角，即当方位角<180°时，反方位角＝此方位角+180°，当方位角>180°时，反方位角＝此方位角-180°。再根据两副山体阴影地形图通过上述步骤寻找到各自对应的视觉粗糙度，最后将所述两副山体阴影地形图的视觉粗糙度进行比较，取数值最大的视觉粗糙度进行后续的计算，从而选取的视觉粗糙度出自信息损失较少的山体阴影地形图，以此来提高后续计算的精度。

(2)、参考图2，根据视觉粗糙度和数字高程模型中所有像素的对应高程值计算出高程地形特征图，具体包括如下步骤；

(2.1)、取数字高程模型中全部像素高程值的百分之三十分位数q₁与I类的临界高程值γ₁中的最小值α，取数字高程模型中全部像素高程值的百分之七十分位数q₂与II类的临界高程值γ₂中的最大值β，即α＝min{H_q1,γ₁},β＝max{H_q2,γ₂}；

且I类的临界高程值γ₁为30，所述II类的临界高程值γ₂为60。

本实施例中冰盖被定义为I类，包括冰川和冰架；海洋被定义为II类，包括海水、海冰、冰山和混杂冰。

(2.2)、根据视觉粗糙度R、数字高程模型中像素的高程值H、α和β计算出像素对应的高程地形特征值，具体通过如下公式获得：

式中，Median{H,α,β}为高程值H、α和β中的中位数；

(2.3)、重复步骤(2.2)直至所有的像素高程地形特征值被计算出；

(2.4)、根据高程地形特征值得到高程地形特征图，即将高程地形特征值中小于其他99％特征值的像素值赋值为0，高程地形特征值中大于其他99％特征值的像素值赋值为255，其他像素根据高程地形特征值的大小从小到大依次在(0,255)范围内赋值。

数字高程模型中像素的噪声值会使山体阴影地图显得更粗糙，而海水即II类的噪声一般比冰川即I类的粗糙度大，同样海水的高程值一般比冰川小得多，而且在海水和冰川交界的地方高程值会有一个跳变，因为冰川一般会在交界处崩解形成断崖，即海水的粗糙度大于冰川，海水的高程值小于冰川，为了进一步扩大二者之间的差异，用粗糙度比上高程值，这样就将二者进一步的区分开了，然后根据高程特征值创建特征图，就能有效的增强I类与II类之间的对比度。

(3)、参考图3，根据基于博弈论的分类方法对高程地形特征图中所有像素进行海洋、冰盖两类分类，基于博弈论的分类方法为：将所有像素都被视为一个理性玩家，每个玩家的策略集为L＝{0，1，2}，其中0表示未分类的，1表示类I，2表示类II，对于每个策略，对应的策略成本为E_i(x_i)。考虑到N个玩家的所有策略，总策略成本E(x)是每个玩家的决策的成本值E_i(x_i)的总和，通过最小化总体成本E(x)，可以得到所有像素的最优策略集合；

(3.1)、将像素的高程值H小于α的像素标记为I类；将像素的高程值H大于β的像素标记为II类；

(3.2)、将其外他未分类的像素，即像素的高程值H大于α且小于β的像素按照使得其自身的策略成本E_i(x_i)最小进行标记，具体E_i(x_i)通过如下公式获得：

式中，i为任意像素，j为i的邻域内的像素，x_i为像素i选择的策略值，x_j为i的邻域内的像素所选择的策略值，r为固定的惩罚值即为255，策略集合{0,1,2}代表未分类、I类、II类标记；f(i)为像素i的高程地形特征值，f(j)为像素j的高程地形特征值，N_i为像素i邻域内的其他像素，min||f(i)-f(j)||为像素i的高程地形特征值f(i)与像素i邻域像素j的高程地形特征值f(j)之间差值的最小值；

步骤(3.3)、重复步骤(3.2)直至所有的点被分类。

参考图4，根据分类结果提取边界线并矢量化，具体为：将不同类别的相邻像素作为边界区域，再将边界区域转换为1像素宽的边界线，然后将1像素宽的边界线从栅格图像转换成矢量文件，该矢量文件为最终的冰川崩解前沿定位结果，将边界线提取并矢量化，可以直观的展示提取结果，方便工作人员的观看。

Claims (8)

1.一种基于数字高程模型的冰川崩解前沿自动提取方法，其特征在于，该方法至少包括如下步骤：

(1)、根据山体阴影可视化技术创建数字高程模型的山体阴影地形图，再根据山体阴影地形图计算出视觉粗糙度；

(2)、根据视觉粗糙度和数字高程模型中所有像素的对应高程值计算出高程地形特征图；取数字高程模型中全部像素高程值的百分之三十分位数q₁与I类的临界高程值γ₁中的最小值α，取数字高程模型中全部像素高程值的百分之七十分位数q₂与II类的临界高程值γ₂中的最大值β，即α＝min{H_q1,γ₁},β＝max{H_q2,γ₂}；

(3)、根据基于博弈论的分类方法对高程地形特征图中所有像素进行海洋、冰盖两类分类，海洋、冰盖两类即I类与II类；

步骤(3)所述高程地形特征图中所有像素通过如下步骤分类：

(3.1)、将像素的高程值H小于α的像素标记为I类；将像素的高程值H大于β的像素标记为II类；

(3.2)、将其外他未分类的像素，即像素的高程值H大于α且小于β的像素按照使得其自身的策略成本E_i(x_i)最小进行标记，具体E_i(x_i)通过如下公式获得：

式中，i为任意像素，j为i的邻域内的像素，x_i为像素i选择的策略值，x_j为i的邻域内的像素所选择的策略值，r为固定的惩罚值即为255，策略集合{0,1,2}代表未分类、I类、II类标记；f(i)为像素i的高程地形特征值，f(j)为像素j的高程地形特征值，N_i为像素i邻域内的其他像素，min||f(i)-f(j)||为像素i的高程地形特征值f(i)与像素i邻域像素j的高程地形特征值f(j)之间差值的最小值；

步骤(3.3)、重复步骤(3.2)直至所有的点被分类；

(4)、根据分类结果提取边界线并矢量化。

2.根据权利要求1所述的基于数字高程模型的冰川崩解前沿自动提取方法，其特征在于：步骤(1)中所述山体阴影可视化技术采用双光源入射生成两副山体阴影地形图，且光源之间入射方向互为相反方位角。

3.根据权利要求2所述的基于数字高程模型的冰川崩解前沿自动提取方法，其特征在于：步骤(1)中所述视觉粗糙度R(n)通过寻找山体阴影地形图中n×n的窗口中像素值最大和最小的像素点B_max，B_min得到，具体为R(n)＝B_max-B_min。

4.根据权利要求2或3所述的基于数字高程模型的冰川崩解前沿自动提取方法，其特征在于：将所述相反方位角入射光源生成的两副山体阴影地形图的视觉粗糙度进行比较，取数值最大的视觉粗糙度进行计算高程地形特征图。

5.根据权利要求1所述的基于数字高程模型的冰川崩解前沿自动提取方法，其特征在于、步骤(2)中所述高程地形特征图通过如下步骤获得：

(2.1)、取数字高程模型中全部像素高程值的百分之三十分位数q₁与I类的临界高程值γ₁中的最小值α，取数字高程模型中全部像素高程值的百分之七十分位数q₂与II类的临界高程值γ₂中的最大值β，即α＝min{H_q1,γ₁},β＝max{H_q2,γ₂}；

(2.2)、根据视觉粗糙度R、数字高程模型中像素的高程值H、α和β计算出像素对应的高程地形特征值，具体通过如下公式获得：

式中，Median{H,α,β}为高程值H、α和β中的中位数；

(2.3)、重复步骤(2.2)直至所有的像素高程地形特征值被计算出；

(2.4)、根据高程地形特征值得到高程地形特征图，即将高程地形特征值中小于其他99％特征值的像素值赋值为0，高程地形特征值中大于其他99％特征值的像素值赋值为255，其他像素根据高程地形特征值的大小从小到大依次在(0,255)范围内赋值。

6.根据权利要求5所述的基于数字高程模型的冰川崩解前沿自动提取方法，其特征在于：所述I类的临界高程值γ₁为30，所述II类的临界高程值γ₂为60。

7.根据权利要求6所述的基于数字高程模型的冰川崩解前沿自动提取方法，其特征在于，步骤(3)所述基于博弈论的分类方法为：将所有像素都被视为一个理性玩家，每个玩家的策略集为L＝{0，1，2}，其中0表示未分类的，1表示类I，2表示类II，对于每个策略，对应的策略成本为E_i(x_i)；考虑到N个玩家的所有策略，总策略成本E(x)是每个玩家的决策的成本值E_i(x_i)的总和，通过最小化总体成本E(x)，可以得到所有像素的最优策略集合。

8.根据权利要求1所述的基于数字高程模型的冰川崩解前沿自动提取方法，其特征在于：所述边界线通过如下方式获得：将不同类别的相邻像素作为边界区域，再将边界区域转换为1像素宽的边界线，然后将1像素宽的边界线从栅格图像转换成矢量文件。

Images