一种形态学滤波时间序列冰面湖提取方法

杨晓彤，赵彬如，邢喆，张峰，王朝阳

(国家海洋信息中心，天津 300171)

0 引言

格陵兰冰盖表面消融影响冰盖物质平衡，是导致全球海平面上升的重要因素[1]。冰面湖(supraglacial lakes，SGLs)作为冰盖表面消融的主要表现形式，在建立冰上和冰下水文系统之间的水文联系方面发挥着十分重要的作用[2]。冰面湖主要通过两种方式导致格陵兰冰盖的物质损失[3-4]。一方面它们通过影响地表融化速率，导致更多的消融;另一方面，由于水力压裂作用(hydrofracturing)，冰面湖可以将融水快速排放至冰盖底部，进而加速冰体运动，从而促进溢出冰川崩解。因此，冰面湖成为了研究格陵兰冰盖对气候变化响应的一个重要因子[5]，尤其对消融期内冰面湖的形成、演变和消失过程开展监测，有助于理解格陵兰冰盖表面融水的存储机制[6]，具有十分重要的科学意义。

格陵兰冰面湖实地调查难度大、危险性高，卫星遥感技术为格陵兰冰面湖动态监测提供了宏观、综合、动态的对地观测手段。冰面湖范围提取是开展冰面湖分布特征、面积和体积变化监测等工作的基础。已有的冰面湖范围提取研究大多基于Landsat、ASTER、MODIS等中低分辨率卫星开展[7-9]，然而，随着卫星遥感影像空间分辨率的不断提高，有必要基于时间序列高分辨率卫星遥感影像研究冰面湖范围的精细化提取方法，兼顾工作效率和提取精度，完成覆盖整个消融期的冰面湖信息提取，为研究冰面湖的时空变化规律、监测冰面湖水量变化等提供重要的基础数据支持。

常用的冰面湖提取方法有人工[10-13]和自动/半自动两种[14-17]，其中自动/半自动方法由于能够节省大量的人工和时间成本，得到了更加广泛的应用，尤其适用于时间序列影像的冰面湖提取问题[18-20]。冰面湖信息提取本质上是冰雪环境背景下的水体信息提取问题，大多数的自动/半自动方法基于归一化水体指数(normalized difference water index，NDWI)并设定合适的阈值来实现[21-23]。然而，将这一基本思路应用于时间序列高分辨率影像的冰面湖信息自动化、高效提取则仍需解决三方面问题。一是水体提取的阈值确定问题。目前的水体-非水体阈值确定方法大多较为复杂，如多阈值分割法[24]、OTSU法等，均需要依靠编程技术实现，且针对不同时相的数据，往往需要设置不同的提取阈值，工作效率较低。因此，有必要研制一种适用于整个消融期时间序列影像，并且简单、易实现的阈值确定方法，应用单一阈值完成时间序列影像上水体信息的高效、快速提取。二是冰面湖与其他水体信息的区分问题。随着气温的升高，冰盖表面融水量达到高峰，会形成冰面河，此外还有冰面融水、冰水混合物等区别于冰面湖的干扰信息，通过阈值很难将这些地物与冰面湖进行自动区分[25-28]。如何从提取的水体信息中自动获取冰面湖仍有待进一步的研究。针对完成水体提取的二值化影像，可将专门用于二值图像处理的数学形态学算法引入冰面湖提取研究[29]。数学形态学可以对遥感影像进行图像增强、图像分割、边缘检测、特征分析、几何分析等[30]，已得到了广泛应用[31-34]。然而，形态学算法的应用效果主要受到滤波方法的选择和滤波窗口大小的影响[35-36]。因此，将其应用于格陵兰冰面湖提取研究，还需对滤波方法的选择和滤波窗口大小的设置进行探讨。三是提取结果的精度评价问题。对于基于高分辨率影像的信息提取，尤其对于时间序列的提取结果，很难获取更高分辨率的影像用于精度评价。格陵兰冰盖属于实测困难区，很难通过现场验证的方法对提取结果进行精度评价。同时，人员很难到达，利用无人机等新型测绘手段则需要较高的人力和物力成本[37]。针对该问题，王辉等[38]提出了基于目视解译结果的面积统计对比法。该方法通过目视解译完成对冰面湖范围的屏幕数字化，然后通过计算冰面湖自动提取结果与目视解译结果的重叠区面积占目视解译结果的百分比来评价自动提取结果的准确性。然而，利用重叠区面积所占百分比，存在无法有效评价冰面湖自动提取结果与目视解译结果的真实偏差情况的可能。因此，有必要研究针对高分辨率时间序列影像的冰面湖提取结果的精度评价方法。

本研究围绕目前高分辨率时间序列影像冰面湖信息提取方法存在的上述问题，开展基于形态学滤波的时间序列冰面湖信息提取方法研究。以格陵兰东北部一个研究区为例，基于覆盖2017年整个消融期的九幅时间序列WorldView多光谱影像，研究适用于整个消融期的单一阈值水体信息提取方法，探讨基于开运算的冰面湖范围自动获取方法，并提出一种基于欧式距离的精度评价方法，通过定量评价自动提取结果与目视解译结果的偏差情况，确定最优的开运算滤波窗口，并最终得到满足精度要求的时间序列冰面湖信息提取结果，验证了研究方法的适用性、鲁棒性和高效性。

1 研究区与数据

1.1 研究区概况

研究区位于格陵兰岛东北部一个快速冰流区，经纬度范围为23°4′56″W～23°35′8″W、78°59′20″N～79°4′8″N，总面积约163.2 km2。2017年6月至8月，该区域内共有7个冰面湖在此形成、发育并消失。选择该区域开展冰面湖信息提取研究，能够验证本文提出的方法对于不同发育阶段冰面湖范围提取的适用性，具有一定的典型性。图1为2017年7月15日的WorldView 3遥感影像图，显示了在盛夏时期该区域内的冰面湖形态特征及分布情况。

图1 研究区概览

1.2 研究数据及预处理

研究数据包含9幅四波段WorldView卫星遥感影像(表1)，时相从2017年6月14日至2017年8月26日，基本覆盖了格陵兰岛的一个消融期，同时相邻日期影像的时间间隔不超过10 d，适用于监测冰面湖在整个消融期内的消长变化情况。

表1 本研究所用WorldView时间序列影像

影像预处理步骤主要包括辐射校正、FLAASH 大气校正和基于RPC模型的正射校正等，将全部影像设置为WGS84 UTM 27N投影。本文仅选择2 m分辨率的多光谱波段进行冰面湖信息提取研究，因此对于自主定位精度达到3.5 m左右的WorldView影像，仅需对不同时相的影像进行相对校正，消除影像之间的位置偏差，以更好地监测冰面湖的空间变化特征。

2 研究方法

基于形态学滤波的时间序列冰面湖信息提取的研究思路是：首先，研究基于归一化水体指数统计分析的单一阈值确定方法，适用于消融期内九幅时间序列影像上的水体提取，完成水体-非水体二值化分割；然后，基于开运算，完成基于多种滤波窗口大小的冰面湖自动提取；最后，提出一种基于欧式距离的冰面湖提取精度快速、自动化、定量评价方法，确定最优滤波窗口并得到最优的时间序列冰面湖信息提取结果，并应用目视判读对提出的精度评价方法和提取结果进行进一步验证。研究流程图见图2。

图2 技术流程图

2.1 改进型归一化水体指数

传统的归一化水体指数，基于水体在近红外波段具有极低的反射率，而在蓝绿波段却具有极高的反射率的特点，通过增强这一差异实现水体信息的有效提取。然而，NDWI更加适用于陆地表面的开阔水域提取，对于冰雪背景下的格陵兰地区，该指数的敏感性较低，文献[6]针对格陵兰的冰雪环境特点对NDWI进行了改进，改进后的归一化水体指数(modified normalized difference water index，MNDWI)有效提高了以冰雪为主背景的水体信息提取精度。

MNDWI选择蓝光波段和红光波段进行计算，适用于本研究的四波段WorldView多光谱影像，因此本文选择该指数结合单一阈值完成研究区的时间序列水体信息提取。

2.2 水体阈值确定

计算影像的MNDWI结果后，可通过设置合适的阈值完成水体信息提取。但是已有的阈值设定方法往往较复杂，且对于时间序列影像要耗费更多的时间和人力资源。本研究尝试通过采集研究区内的典型地物样本，并对它们的MNDWI数值进行统计分析，探求单一阈值对于时间序列影像水体提取的适用性。

首先，通过目视判读确定研究区内的典型地物，主要包括冰、水、冰水雪的混合物等，为了方便阈值确定，将研究区内地物确定为水体和非水体两大类。然后，采用人机交互的方式在九幅影像上选取水体和非水体典型样本，并统计其MNDWI数值，并绘制统计直方图。其中水体像元的选择顾及多种水深，以提高阈值对水体提取的准确度。对于典型地物样本，一般情况下会有非水体的MNDWI数值小于水体的MNDWI数值，因此本研究通过统计九幅影像中全部非水体样本中的MNDWI最大值以及全部水体样本中的MNDWI最小值，并选择其中的较小者作为适用于全局的水体提取阈值，则阈值t可表示为式(1)。

t=min(max(MNDWInonwater)，min(MNDWIwater))

(1)

式中：max(MNDWInonwater)表示全部非水体样本中MNDWI数值中的最大值；min(MNDWIwater)表示全部水体样本中MNDWI数值中的最小值。在实际操作中，将统计直方图中的前100个像元和后100个像元剔除，以排除异常值。

本研究选择较小的阈值t是为了尽可能减少水体提取的漏提问题。尽管较小的阈值可能导致一些伪水体信息的误提取，但水体信息的漏提将直接影响后续的冰面湖信息提取精度，而水体提取的错分问题可通过后续的形态学算法进行修正。

2.3 基于形态学滤波的冰面湖提取

基于确定的阈值完成水体-非水体二值化后，采用数学形态学滤波，完成时间序列冰面湖边界范围的提取。常用的形态学算法有腐蚀、膨胀、开运算和闭运算等。

本文采用开运算进行冰面湖边界提取，开运算即先对二值化影像进行腐蚀计算，然后再进行膨胀计算。选择开运算进行冰面湖提取，主要利用开运算的三方面优势[39-40]：一是消除亮度较高的细小区域，即去除椒盐噪声，在本研究中可用于剔除较小的阈值t引起的伪水体信息的误提问题，或者去除尚未形成冰面湖的其他细小水体；二是在纤细点分离物体，在本研究中可用于分割与冰面湖相连的冰面河，提高冰面湖边界的提取精度；三是可以在不明显改变其面积的情况下对较大物体进行边界平滑，在本研究中可以平滑冰面湖提取边界。

本文的开运算应用ENVI软件自动实现，需要设置的参数主要包括滤波窗口大小(kernel size)、滤波次数(cycles)、滤波器类型(style)、滤波核权重(editable kernel)。由于滤波窗口，即结构元素的大小是影响形态学滤波运算结果的直接因素，本研究将其他参数均设置为默认值，仅重点探讨滤波窗口大小的变化对冰面湖提取结果的影响。在ENVI中，滤波窗口大小可以设置为大于等于3的奇数。通常，随着滤波窗口的增大，影像滤波结果会出现过平滑等问题，导致提取结果精度的降低。因此，首先，将滤波窗口依次设置为3、5、7、9、11等数值，得到多个冰面湖提取结果；然后，通过目视判读，直接将出现过平滑的结果进行剔除，确定最优滤波窗口的大致范围；最后，对无法直接比较提取精度的结果进行精度评价，从而确定最优的滤波窗口。

2.4 精度评价方法

对于高分辨率的卫星遥感影像，很难获取更高分辨率的影像对信息提取结果进行精度验证，同时格陵兰冰盖属于实测困难区，限制了现场调查的可能，因此基于屏幕数字化的提取结果常常被视为地面真实值带入精度评价中。应用重叠面积比屏幕数字化结果面积作为精度评价指标，无法有效评价提取结果与地面真实结果的实际偏差，因为存在面积统计精度值较高，但提取结果与目视解译结果偏差很大的可能。如图3所示，红色边界为冰面湖信息提取边界，蓝色为屏幕数字化边界，则采用重叠区面积比屏幕数字化结果面积会得出提取精度较高的结论，然而由图3可看出，冰面湖提取结果与屏幕数字化结果之间在局部存在较大差异，提取结果无法满足精度要求。

图3 面积统计法进行精度评价示意图

因此，本文基于屏幕数字化的冰面湖边界，提出一种基于欧式距离的冰面湖信息提取精度的定量化评价方法，以实现对提取边界与地面偏差的度量。具体步骤如下。

步骤1：采用屏幕数字化方法对九幅影像上研究区内的冰面湖进行人工勾绘，得到九期冰面湖精度评价参考数据。由于对于冰面湖范围的界定存在目视解释较为主观的问题，因此采用三名课题组人员独立勾绘并对数字化结果进行交叉检查修改的方法，提高屏幕数字化结果的准确性，并作为地面真实值参与精度评价。课题组人员在勾绘冰面湖时遵照如下判定准则：形状为近圆形或者椭圆形的且面积较大的水体作为冰面湖，对于狭长的辫状水体或面积较小的细碎水体则认为是冰面河等其他开放水体；针对判定为冰面湖的水体，以冰与水的交界线作为其边界，冰与水的交界线在1∶2 000比例尺下进行判读。

步骤2：获取基于本文研究算法提取的冰面湖范围线的节点，作为精度评价的检查点。

步骤3：计算每个检查点到基于屏幕数字化提取的冰面湖边界的垂直距离，即欧式距离。

步骤4：为了排除极端值对精度评价结果的影响，选择通过统计每幅影像上全部检查点的欧式距离的中位数L和标准差s，来定量化评价提取结果与真实结果的吻合程度。图4选择了一个出现错分的冰面湖提取结果作为精度评价方法的示意图。从图4(a)可以看出，提取的冰面湖边界在局部与真实情况存在一定偏差。图4(b)为局部放大图，其中的黄色线段表示每个精度检查点至真实值的欧式距离，即提取结果与真实边界的实际偏离情况。从检查点的分布密度看，本文提出的精度评价方法可以通过对提取的冰面湖边界进行精细化评价，定量评估提取结果与地面真值的一致性。

图4 精度评价方法示意图

3 结果及分析

3.1 基于光谱分析的阈值确定结果

本研究从9幅影像上分别选择约2 000个水体和非水体的纯像元，绘制其MNDWI数值的统计直方图(图5)。从图5可以看出，对于本研究任一时相影像，大部分非水体像元的MNDWI数值都小于水体，表明MNDWI对于冰雪背景下的水体和非水体具有较高的区分度，可以直接应用于本文的水体提取研究。

图5 典型地物的MNDWI分布直方图

通过记录剔除异常值后，记录各时相下MNDWInonwater的最大值和MNDWIwater的最小值(表2)，从中找出max(MNDWInonwater)和min(MNDWIwater)，带入式(1)，得到适用于全部9幅影像的水体-非水体分割阈值。由表2可以看出，9幅影像中max(MNDWInonwater)的值为0.476，min(MNDWIwater)的值为0.421，则确定水体-非水体分割阈值为0.421。

表2 各时相影像水体和非水体像元MNDWI统计值

应用该阈值完成水体分割，得到二值化后的水体分割结果(图6)，图中白色部分代表提取出的水体范围。如前所述，选择一个较低的水体分割阈值，是为了避免水体像元的漏提问题，因为漏提会直接影响最终的冰面湖提取精度。

图6 水体阈值分割结果图

3.2 冰面湖提取结果

1)冰面湖信息提取初步结果确定。在ENVI中，将滤波窗口大小依次设置为3、5、7、9、11等数值对九幅影像执行开运算，并将计算结果转换成矢量格式的线状要素，即可得到冰面湖边界信息。然后，通过目视判读对提取结果进行初步分析。经实验，当滤波窗口设置为11时，多期影像中的多个冰面湖提取结果出现了过平滑现象，精度已无法满足实验要求。因此，将滤波窗口大小分别设置为3、5、7、9的9期36个提取结果作为冰面湖信息的初步提取结果。接下来将通过进一步的精度评价，从中确定最优的冰面湖提取结果和最优的滤波窗口大小。

2)基于欧式距离的精度评价及最优提取结果确定。采用本文提出的基于欧式距离的精度评价方法对36个提取结果进行精度评价，找出最优解。针对36个提取结果，首先提取用于精度评价的检查点，即提取的冰面湖边界线的节点；然后计算检查点到屏幕数字化得到的冰面湖范围线的欧式距离，并计算欧式距离结果的统计中位数L和标准差s，其中L代表冰面湖的提取误差，s表示误差的离散程度；最后，绘制误差线图，分析确定最优滤波窗口和最优提取结果。

图7为对滤波窗口大小依次设置为3、5、7、9得到的提取结果进行基于欧式距离的精度评价后绘制的误差线图。对于消融初期和末期，即6月14日、6月24日和8月26日，冰面湖提取误差L随着滤波窗口值的增大而显著增大。当滤波窗口设置为3时，提取结果与真实值之间的总体误差最小，即为当期影像的最优提取结果。由误差线图还可看出，当滤波窗口设置为3时，6月14日和6月24日检查点的欧式距离计算结果的L均小于1.5 m，且s也较小(小于0.2)，表明提取的所有冰面湖边界信息均有较高的精度。而对于8月26日的提取结果，其与真实结果的总体偏差略大，统计L为9.12 m，且s亦较大(1.7 m)，提示该期提取结果存在整体或局部误差较大的情况。

图7 基于不同滤波窗口值的欧式距离统计误差线图

对于气温较高的盛夏季节，即7月6日至8月19日，冰面湖提取误差L随着滤波窗口值的增大呈现显著减小并趋于稳定的变化规律。滤波窗口设置由3变为5时，冰面湖提取误差显著降低；当滤波窗口设置为7时，提取误差的降低速率减缓；当滤波窗口设置为9时，针对不同影像，提取误差有的略有降低，有的略有升高，提取精度总体仍优于滤波窗口设置为5时的提取结果，但其s较滤波窗口设置为7的提取结果均有增大趋势，表明此时冰面湖提取边界出现局部精度降低的情况。因此认为，针对7月6日至8月19日的影像，滤波窗口设置为7时，其对应的冰面湖提取结果精度最高。由误差线图还可看出，在最优滤波窗口下，除8月1日的提取结果与地面真值之间的总体误差仍较大(7.02 m)外，其他期影像的提取结果与地面真值之间的总体偏差均小于4.5 m，精度检查点的欧式距离计算结果的标准差也较小，提取结果的总体精度较高。

最后，计算9期最优提取结果全部检查点的欧式距离中位数的平均值，得到9期提取结果与真实结果的偏离情况平均小于3.82 m。对于空间分辨率2 m的WorldView多光谱影像，即偏差在2个像元以内，表明本研究确定的9期冰面湖信息提取结果具有较高的提取精度。

图8为经过精度评价确定的基于最优滤波窗口的时间序列冰面湖提取结果。由图8可以看出，8月1日和8月26日的提取结果与实际偏差较大的原因为单个湖泊出现明显的错分或者漏分导致。图9为7月15日冰面湖提取结果放大图，可以更加清晰地展示基于本文方法提取的冰面湖与地面真实情况的吻合程度。

图8 时间序列冰面湖信息提取结果

图9 2017年7月15日研究区冰面湖提取结果放大图

3)最优提取结果精度评价。为了进一步验证基于本文精度评价方法确定的提取结果的可靠性，采用传统人工目视判读的方法对9幅影像上出现的冰面湖进行逐一评价，将冰面湖提取不完整的视为漏分，计为漏分误差；将提取边界与实际情况一致性较差的视为错分，计为错分误差。

应用本文的方法在9幅影像上共累计提取了48个冰面湖边界线，此处仅以面积大于5 000 m2的43个冰面湖进行评价。由于面积较小的冰面湖一般处于形成初期或消失末期，其边界较难界定，因此不计入统计。目视评价结果见表3，冰面湖提取总体精度为88.37%，漏分误差为2.33%，错分误差为9.30%。

表3 基于目视判读的最优提取结果精度评价

本研究唯一出现漏分的是8月26日的一个冰面湖(图10(a))，经分析，导致漏分是因为应用确定的阈值进行水体提取时，出现了漏分问题(图10(b))。图10(c)～图10(f)为出现错分的冰面湖，4个边界提取结果的局部与冰面湖的实际边界存在差异，其中3个错分的冰面湖集中在8月1日的影像上，也和水体阈值分割结果相关。

图10 提取结果出现误差的冰面湖

4 结束语

本文选择格陵兰东北部一个快速冰流区作为实验区，基于MNDWI和单一阈值，通过引入形态学开运算，实现了2017年整个消融期内9幅WorldView多光谱影像上的冰面湖信息提取。通过基于欧式距离的自动化精度评价方法，确定了不同时相对应的最优开运算滤波窗口为3或7，并以此获得最佳的时间序列冰面湖信息提取结果，提取的冰面湖范围与屏幕数字化结果的偏差平均小于2个像元。进一步通过目视判读，验证最终确定的冰面湖信息提取的总精度达88.37%，其中漏分误差为2.33%，错分误差为9.30%，表明基于本文方法提取的时间序列冰面湖信息，其精度满足后续研究需求。

本文提出的水体阈值确定方法简单、易实现，通过对各时相影像上的典型水体和非水体样本的MNDWI值进行统计分析即可获得。精度评价结果表明，本研究确定的单一阈值具有较好的鲁棒性，能够满足多时相影像的水体-非水体二值化分割，适用于时间序列影像的冰面湖信息快速提取。事实上，非水体的MNDWI最大值和水体的MNDWI最小值均可作为水体-非水体的分割阈值，本文选择二者中较小的水体MNDWI最小值(0.421)作为阈值，是为了避免水体信息的漏提问题，因为漏提会直接影响后续冰面湖信息提取的精度。本研究中唯一个发生漏提的冰面湖，即是水体信息提取出现漏提导致的，表明水体信息提取过程中出现的误差会直接影响后续冰面湖信息的提取精度。

基于形态学开运算的冰面湖信息提取方法，通过设置合适的滤波窗口，可以得到满足精度要求的冰面湖提取结果。应用开运算进行时间序列冰面湖信息提取的优势主要有：算法运行效率高，对于时间序列影像信息提取问题，可以大大提高其工作效率；可以剔除由于选择较小阈值导致的水体提取结果中的椒盐噪声，提高提取结果精度；能够有效区分冰面湖和其他冰面水体，在提高冰面湖提取精度的同时得到边界平滑的冰面湖提取结果。

通过计算提取结果范围线所有节点至屏幕数字化冰面湖范围线的欧式距离的精度评价方法，实现对冰面湖信息提取结果的精细化定量评价，并通过精度评价结果快速获得最优的冰面湖提取结果，不仅适用于时间序列提取结果的精度评价，而且可用于确定最优的开运算滤波窗口。与基于目视判读的精度评价方法相比，该方法免去了逐一进行人工目视判读的时间，效率更高，评价结果也更加客观。本研究中，在消融初期(6月14日、6月24日)和末期(8月26日)，冰面湖提取误差随着滤波窗口值的增大而显著增大，最优滤波窗口值为3；对于消融高峰期的影像，在滤波窗口值由3增大到9的过程中，冰面湖提取误差随着滤波窗口值的增大而显著减小并趋于稳定，表明滤波窗口对提取精度的影响趋于稳定。在综合分析误差离散程度的基础上，确定最优滤波窗口值为7。

尽管应用单一阈值提取水体信息影响了本研究中个别时相影像的水体提取精度，但是基于单一阈值完成水体信息的快速提取，可以在一定程度上提高工作效率，尤其适用于长时间序列的冰面湖信息提取。此外，通过后续的形态学算法可以修正单一阈值导致的水体提取误差。

综上所述，本文提出的冰面湖提取方法具有易操作、工作效率较高等特点，尤其适用于时间序列影像上的冰面湖提取问题。针对本研究采用的覆盖整个消融期的9幅WorldView时间序列多光谱影像，方法的适用性较强，能够对研究区内一个消融期内发育各个阶段的冰面湖进行较精确的提取。后续，将进一步开展该方法在其他星源影像和其他地区的适用性研究，研制鲁棒性更好、精度更高的冰面湖信息提取方法，为开展冰面湖水量反演研究提供更精准的数据支撑。