一种综合阈值分割和分水岭变换算法的冰面湖提取方法与流程

本发明属于地球地质研究领域，特别涉及一种综合阈值分割和分水岭变换算法的冰面湖提取方法。

背景技术：

冰面湖是位于冰盖(冰川)表面的湖泊，是因冰盖(冰川)表面差异消融形成的临时性积水洼地。冰面湖多出现在规模较大的冰盖(冰川)消融区表面，如格陵兰冰盖、南极大陆冰盖、珠穆朗玛峰北侧绒布冰川消融区。每年消融期，当冰盖(冰川)表面温度达到融点后形成融水，融水产生后在地形等因素作用下输送，形成规模庞大、结构复杂的冰盖表面水文系统，这一系统包括冰面湖(supraglaciallake)、冰面河(supraglacialstream)、锅穴(moulin)和冰裂隙(crevasse)等水文要素。冰面湖作为连接冰盖(冰川)地上与地下水文系统的重要桥梁，是反映冰盖(冰川)与气候变化之间关系的重要指示因素。以格陵兰冰盖为例，近20年来格陵兰冰盖的物质损失是造成海平面上升的主要原因，而造成格陵兰冰盖物质损失的主要原因之一就是冰盖表面消融。然而针对格陵兰、南极、珠穆朗玛峰等研究区，开展冰面水文系统实地调查的难度大、危险性高，卫星遥感技术作为一种宏观、综合、动态的对地观测手段，已成为冰面水文系统研究的新契机、新领域。依托卫星遥感技术研究冰面湖的提取方法，能够为监测冰面湖的演化过程、反演冰面湖水深和体积提供重要基础数据支撑，对于研究冰盖(冰川)表面融水的形成、存储、输送和释放机制具有重要意义。

目前，基于遥感技术提取冰面湖范围的方法主要有人工手动勾绘和自动阈值提取两种。人工手动勾绘是通过人机交互的方式，通过屏幕数字化开展。该方法能够在一定程度上确保冰面湖提取的准确度，但是该方法效率极低，并且提取结果主要取决于操作人员的专业水平，提取结果缺乏客观性。此外，应用该方法进行冰面湖的长时序、大面积观测则要消耗更多的人力和物力。因此，基于遥感影像进行冰面湖的自动化提取是研究的重点领域。该方法的主要思路是通过波段比值，归一化水体指数(normaliseddifferencewaterindex,ndwi)计算等方法，综合分析冰面湖的光谱特点，设置合适阈值,通过阈值分割区分冰面水体和非水体进而提取冰面湖。该方法面临的两个主要挑战：一是确定一个精确的阈值用于区分水体和非水体象元；二是从水体象元中区分冰面湖与其他冰面水文要素，如冰面湖，冰面融雪等。为了得到更加精确的阈值，otsu阈值确定法等各类阈值方法被引入冰面湖提取研究中；对于冰面湖与其他水体象元的区分问题，有些学者采用多阈值分割方法，应用不同的阈值对冰面湖与其他冰面水文要素进行区分。上述方法大多需要研究者具备较高的编程水平，并且针对不同的影像和研究区需要重复进行复杂的分析和试验，方法的鲁棒性较差。此外，即使设置合适的阈值，仍不可避免对非水体象元的误提取。

技术实现要素：

本发明的目的是提出一种综合阈值分割和分水岭变换算法的冰面湖提取方法，该方法阈值确定方法简单易实现，鲁棒性好，可从水体提取结果中剔除误分象元，并进一步区分冰面湖与冰面水文要素。

如上构思，本发明的技术方案是：一种综合阈值分割和分水岭变换算法的冰面湖提取方法，其特征在于：包括如下步骤：

①遥感影像预处理；

②改进型归一化水体指数计算，应用蓝光波段和红波段进行归一化指数计算，具体计算公式如下：

③基于改进型归一化水体指数的计算结果，在影像上选择水体和非水体的典型样本；

④将样本统计信息导入统计分析软件，对采集的水体，非水体典型样本的mndwi的数值进行直方图统计分析；

⑤选择步骤④记录的两个数值中较小的一个作为最终的阈值t，即：

t＝min(mndwi水体min，mndwi非水体max)；

⑥应用阈值t进行阈值分割，提取水体象元，对归一化水体指数计算结果进行二值化，其中水体象元赋值为1，非水体象元赋值为0，具体判断标准如下：

其中p为mndwi影像的象元值，t为阈值；

⑦应用基于影像梯度图的分水岭变换算法对二值化后的影像进行影像分割，得到矢量格式的影像分割结果；

⑧应用gis软件中的面转线工具，将分割后得到的面状矢量数据转为线状数据，然后从中提取出冰面湖的线状矢量范围，即可得到冰面湖的边界范围。

进一步，所述步骤①首先应用卫星遥感影像自带的rpc参数对影像进行正射校正，然后应用fastline-of-sightatmosphericanalysisofhypercubes(flaash)方法对影像进行大气校正，最后对影像进行投影设置，具体实施步骤如下：

a.影像正射校正：在通用遥感影像处理软件中打开待校正的遥感影像，选择影像自带的rpc参数，设置参考数字高程模型，根据设置影像的输出分辨率，采用重采样方法完成影像正射校正；

b.影像大气校正：在通用遥感影像处理软件中应用fastline-of-sightatmosphericanalysisofhypercubes(flaash)方法对影像进行大气校正；

c.影像投影设置：根据研究区所在位置，在gis或者遥感软件中对影像进行投影设置。

进一步，所述步骤②的具体实施步骤是：

a.在遥感软件中应用波段运算工具，编写mndwi计算公式，

(float(bblue)-float(bred))/(float(bblue)+float(bred))

b.将影像代入公式，即可计算得到mndwi计算结果。

进一步，所述步骤③的具体实施步骤是：

a.选择通用遥感软件中的roi选择工具，设置水体和非水利两类样本；

b.采用人机交互的方法，在遥感影像上选择水体和非水体的典型象元；c.

应用统计分析工具，计算得到水体和非水体样本的统计信息。

进一步，所述步骤④的具体实施步骤是：在统计典型样本的mndwi指数值时，将统计直方图中水体和非水体各自的前100个象元和后100个象元剔除，然后分别记录剩余所有水体典型样本中mndwi的最小值(mndwi水体min)和剩余所有非水体典型样本的mndwi最大值(mndwi非水体max)。

进一步，所述步骤⑥的具体实施步骤是：

a.在遥感软件中应用波段运算工具，编写二值化计算公式，如:

b1>t

该公式默认将所有值大于t的象元赋值为1，小于t的赋值为0；

b.将影像代入公式，即可计算得到二值化结果。

进一步，所述步骤⑦的具体实施步骤是：

a.计算影像梯度图，即通过边缘检测得到各个波段的梯度影像，然后将多个波段的梯度影像融合成一个单波段梯度图；

b.基于尺度参数对影像梯度值进行替换，首先需要得到梯度图的累计直方图，累计直方图横轴为梯度值，纵轴为梯度值从0开始梯度值范围内所有像元个数占总像元个数的百分比，即0％-100％。通过设置尺度参数，将小于该参数对应的梯度值gt1的像元值均替换为gt1；

c.分水岭变换，将改进后的梯度图经分水岭变换，即可得到连续、封闭的对象边界。

进一步，所述步骤⑦梯度影响融合的方法有欧式距离法，梯度值求和法和梯度值求最大值法。

进一步，所述重采样方法主要包括最近邻法、双线性插值法和三次卷积法。

进一步，所述水体和非水体的典型象元各选择2000-5000个。

本发明具有如下的优点和积极效果：

1、本发明提出的综合阈值分割和分水岭变换算法的冰面湖提取方法能够自动化、高精度的完成冰面湖的自动化提取；

2、本发明中设计的阈值确定方法简单易实现，能够有效区分水体与非水体象元，不需要编程即可快速获取阈值。并且该阈值分割方法适用于各类影像和研究区，鲁棒性较好，解决了现有技术中存在的阈值确定方法繁琐、迁移性差的问题。

3、针对水体提取后仍存在的非水体象元，以及冰面湖与其他冰面水文要素难以区分的问题，本发明通过引入分水岭变换算法，可充分利用冰面湖的几何形态特征，在保留精细化的冰面湖边界的基础上，将冰面湖与其他冰面水文要素进行有效区分，并对非水体象元进行了有效剔除，大大降低了冰面湖提取结果对阈值精确性的依赖程度，提高了基于遥感影像提取冰面湖的精度。

附图说明

图1是本发明的技术流程图；

图2是基于影像梯度图的分水岭变换算法的影像分割技术流程图。

具体实施方式

下面结合附图和具体实施方式对本技术方案的具体实施方式进行详细说明。

如图1所示，一种综合阈值分割和分水岭变换算法的冰面湖提取方法，包括如下步骤：

步骤1：影像预处理。首先应用卫星遥感影像自带的rpc参数对影像进行正射校正，然后应用fastline-of-sightatmosphericanalysisofhypercubes(flaash)方法对影像进行大气校正，消除大气对影像光谱特征的影响，最后对影像进行投影设置。具体实施步骤如下：

步骤1-1：影像正射校正。在通用遥感影像处理软件中打开待校正的遥感影像，选择影像自带的rpc参数，设置参考数字高程模型，根据设置影像的输出分辨率，选择重采样方法完成影像正射校正。所述重采样方法主要包括最近邻法、双线性插值法、三次卷积法等，本实施例采用双线性插值法；

步骤1-2：影像大气校正。在通用遥感影像处理软件中应用fastline-of-sightatmosphericanalysisofhypercubes(flaash)方法对影像进行大气校正。根据研究区的实际情况，输入影像中心点坐标、传感器类型、卫星过境时间、区域高程等信息；然后选择合适的大气模型，由于冰面湖一般位于两极地区和海拔较高的喜马拉雅山脉等，故大气模型可以选择极地冬季(saw)或中纬度冬季(mlw)，选择合适的气溶胶模型，气溶胶模型一般包括无气溶胶、城市气溶胶、乡村气溶胶、海洋气溶胶和对流层气溶胶模型等；最后根据影像波段情况选择大气反演波段，完成影像大气校正；

步骤1-3：影像投影设置。根据研究区所在位置，在gis或者遥感软件中对影像进行投影设置。

步骤2：改进型归一化水体指数(modifiednormalizeddifferencewaterindex，mndwi)计算。已有研究已证明，针对以冰雪为背景的卫星遥感影像特点对传统的归一化水体指数进行改进，应用蓝光波段和红波段进行归一化指数计算，能够更有效的区分水体与非水体。具体计算公式如下：

具体实施步骤如下：

步骤2-1：在遥感软件中应用波段运算工具，编写mndwi计算公式，如:

(float(bblue)-float(bred))/(float(bblue)+float(bred))，

步骤2-2：将影像代入公式，即可计算得到mndwi计算结果。

步骤3：基于改进型归一化水体指数的计算结果，在影像上选择水体和非水体的典型样本。由于不同水深的水体存在光谱差异，导致指数计算结果差异较大，因此在采集水体典型样本时应顾及不同的水深范围。对于非水体区域，各地物类型如冰、雪、裸岩等都应采集一定数量的典型样本。具体实施步骤如下：

步骤3-1：选择通用遥感软件中的roi选择工具，设置水体和非水利两类样本；

步骤3-2：采用人机交互的方法，在遥感影像上选择水体和非水体的典型象元。对于水体样本，由于不同水深的水体存在光谱差异，导致指数计算结果差异较大，因此在采集水体典型样本时应顾及不同的水深范围。对于非水体样本，影像内除水体以外的其他各地物类型如冰、雪、裸岩等都应采集。水体和非水体的典型样本各选择2000-5000个象元即可；

步骤3-3：应用roi统计工具，得到两类样本的统计信息，包括最大值，最小值，平均值，累积频率等。

步骤4：将样本统计信息导入统计分析软件，对采集的水体，非水体典型样本的mndwi的数值进行直方图统计分析。为了消除异常象元对统计结果的影响，在统计典型样本的指数值时，将统计直方图中水体和非水体各自的前100个象元和后100个象元剔除，然后分别记录剩余所有水体典型样本中mndwi的最小值(mndwi水体min)和剩余所有非水体典型样本的mndwi最大值(mndwi非水体max)。

步骤5：为了保证尽可能多的水体信息被提取出来，选择上述记录的两个数值中较小的一个作为最终的阈值t，即：

t＝min(mndwi水体min，mndwi非水体max)

步骤6：应用阈值t进行阈值分割，提取水体象元。具体方法为对归一化水体指数计算结果进行二值化，其中水体象元赋值为1，非水体象元赋值为0。具体判断标准如下：

其中p为mndwi影像的象元值，t为阈值。

具体实施步骤如下：

步骤6-1：在遥感软件中应用波段运算工具，编写二值化计算公式，如:

b1>t

该公式默认将所有值大于t的象元赋值为1，小于t的赋值为0；

步骤6-2：将影像代入公式，即可计算得到二值化结果。

步骤7：应用分水岭变换算法对二值化后的影像进行影像分割，得到矢量格式的影像分割结果。分水岭变换算法，亦可称为基于影像梯度图的分水岭变换算法。该算法的基本思想是：首先计算影像的梯度图，然后将梯度图像的每一个梯度值看成一个海拔高度，则整个影像可以看成一个拓扑地形。梯度图中的每一个局部极小值及其影响区域构成一个汇水范围，即“集水盆”。假设水位从“集水盆”最低点开始上涨，当水位上升到相邻“集水盆”的界线时，则停止上涨。控制水位上涨的相邻“集水盆”之间的界线，即“分水岭线”也就是分割对象的边界。具体实现步骤如图2所示。

步骤7-1：计算影像梯度图，即进行边缘检测。可选用roberts算子，sobel算子等。通过边缘检测得到各个波段的梯度影像，然后将多个波段的梯度影像合成一个单波段梯度图。梯度融合的方法主要有欧式距离法，梯度值求和法和梯度值求最大值法等；

步骤7-2：基于尺度参数对影像梯度值进行替换。此步骤首先需要得到梯度图的累计直方图，累计直方图横轴为梯度值，纵轴为梯度值从0开始梯度值范围内所有像元个数占总像元个数的百分比，即0％-100％。通过设置尺度参数，将小于该参数对应的梯度值gt1的像元值均替换为gt1；

步骤7-3：分水岭变换。改进后的梯度图经分水岭变换，即可得到连续、封闭的对象边界。

步骤8：应用gis软件中的面转线(featuretopolygon)工具，将分割后得到的面状矢量数据转为线状数据，然后从中提取出冰面湖的线状矢量范围，即可得到冰面湖的边界范围。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。