CN113191989B - 一种针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法 - Google Patents

Abstract

本发明涉及一种针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法，包括步骤：对遥感影像数据进行预处理；低频信息提取；对低频信息进行形态学处理。本发明的有益效果是：引入高通滤波处理，使淤泥质海岸海域侧和陆域测频率信息的差异显著增大，有利于水边线位置的精确度；针对遥感影像高通滤波后的像素值进行直方图统计，通过构建自适应低频信息阈值区间算法，获得了能够随影像自适应变化调整的更加稳健的阈值范围，解决了传统阈值法提取水边线中阈值范围过度依赖影像本身，阈值可迁移性差的问题；与传统边缘检测法提取水边线相比，本发明利用低频信息的边界替代图像的边缘，有效地解决了传统边缘检测法提取水边线结果破碎、不连续的缺陷。

Description

一种针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法

技术领域

本发明属于遥感影像处理领域，尤其涉及一种针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法。

背景技术

海洋是维持陆地生态系统平衡和稳定的生态屏障，陆地是海洋开发和保护的重要依托。海洋和陆地的交界线被称之为海岸线，海岸线对于海洋资源利用、船舶航行以及沿海城市规划等都具有重要的意义，是海岸带的重要标志。传统的海岸线测定主要是依靠人工实地测量，但受地理环境条件限制，许多地方无法到达，并且实地测量工作强度大、效率低，难以实现长时间动态监测。近年来，遥感技术凭借其数据获取容易、数据信息量大、高时空分辨率等优点逐渐成为海岸线测定的重要手段。

由于遥感影像瞬时成像的特点，从遥感影像中直接提取的大多是瞬时水边线，即影像成像时海水与陆地的交界线，而并非实际意义上的海岸线。虽然水边线不是海岸线，但在大多数海岸线提取方法中，海岸线是在提取出水边线后对水边线进行一系列修正后得到的，因此水边线提取是海岸线提取的基础。

目前常用的基于遥感的水边线自动解译方法有边缘检测法、阈值法等。边缘检测法得到的水边线通常存在破碎、不连续的问题，阈值法的关键是阈值的确定，但阈值往往存在较大的不确定性，如不同时相、不同区域以及不同传感器影像的阈值往往差别较大，阈值的可扩展性差，不利于水边线的智能化提取。

淤泥质海岸的近岸水体浑浊，悬浮泥沙含量高；陆侧潮间带普遍存在水斑，岸滩含水量较高，且存在混合像元，导致淤泥质海岸水陆两侧的遥感特征差异小，增加了水边线提取的难度。遥感影像作为一种二维信号，具有信号的高、低频特征。影像的高频信号是影像中强度变化剧烈的区域，通常对应于影像的边缘、细节等部分，如水陆分界区域的水边线。影像的低频信号则是影像中强度变化平缓的区域，对应于影像中的大片均匀色块区域，如空地、大片植被以及水体内部等。高通滤波处理能够突出高频信息，抑制低频信息，从而达到增大高低频信息间差异的目的，因此，水边线的遥感提取通过引入高通滤波方法能够增大淤泥质海岸水陆两侧的遥感特征差异，有助于实现水边线的高精度、智能化提取。

发明内容

本发明的目的是克服现有技术中的不足，提供一种针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法。

这种针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法，包括以下步骤：

步骤1、对遥感影像数据进行预处理操作，将原始影像像素值校正为实际地表反射率值，得到地表反射率数据；

步骤2、低频信息提取；

步骤2.1、对地表反射率数据进行高通滤波：

上式中，m表示遥感影像中像元水平方向的位置；n表示遥感影像中像元垂直方向的位置；H(m,n)表示遥感影像中(m,n)处像元经过高通滤波之后的像素值；DN(m,n)表示遥感影像中(m,n)处像元的原始像素值；F表示高通滤波器；i表示滤波器中的元素水平方向的位置；j表示滤波器中元素垂直方向的位置；d为大于等于1的整数；

步骤2.2、对高通滤波后的图像像素值进行直方图统计；

步骤2.3、计算得到低频信息阈值区间，构建自适应的阈值区间算法，根据阈值区间自适应提取低频信息，低频信息值为1，高频信息值为0；

步骤3、对低频信息进行形态学处理，处理完成后得到水陆分离的二值图像；

步骤3.1、利用形态学开运算填充低频信息中的孤立小点；所述形态学开运算为先进行腐蚀运算，再进行膨胀运算；

步骤3.2、利用形态学闭运算去除低频信息中的细小孔洞；所述形态学闭运算为先进行膨胀运算，再进行采用腐蚀运算填充细小孔洞；

步骤3.3、利用膨胀运算弥补低频信息边界与影像实际水边线之间的像素误差；

步骤4、对水陆分离的二值图像进行后处理：将水陆分离的二值图像转为矢量格式，按照矢量格式的多边形面积及属性进行筛选：

(area≥0.05)and(GRIDCODE＝1)

上式中，area表示多边形面积，单位为平方千米，GRIDCODE代表多边形属性，取值为1表示水体；

采用人机交互解译进行人工检核，确定最终水边线。

作为优选，步骤2.2具体包括以下步骤：

步骤2.2.1、将直方图划分为N个像素区间，求得直方图的间隔S：

上式中，H_max为高通滤波后像素值的最大值，H_min为高通滤波后像素值的最小值，N为像素区间个数；

步骤2.2.2、统计各像素区间的占比，按降序排列；

步骤2.2.3、筛选出百分比最高的前B个像素区间，确定直方图中百分比最高的前1％个像素区间；B的值由如下公式确定：

上式中，

为向下取整符号，A为直方图中百分比最高的前1％个像素区间。

作为优选，步骤2.3具体包括以下步骤：

步骤2.3.1、假设前B个像素区间中最小像素值为H_{B_min}，则低频信息阈值区间表示为如下公式：

[H_{B_min},H_{B_min}+B×S]

上式中，H_{B_min}为前B个像素区间中最小像素值；B为百分比最高的像素区间个数，S为直方图的间隔；

步骤2.3.2、构建自适应的阈值区间算法为：

(H ge H_{B_min})and(H le(H_{B_min}+B×S))

上式中，H表示高通滤波后的像素值，ge表示大于等于，le表示小于等于，H_max为高通滤波后像素值的最大值，H_min为高通滤波后像素值的最小值；B为百分比最高的像素区间个数，S为直方图的间隔。

作为优选，步骤1中遥感影像数据为光学遥感影像数据，光学遥感影像数据含有近红外波段。

作为优选，步骤1中对遥感影像数据进行预处理的操作包括几何纠正、辐射定标和大气校正。

作为优选，步骤2.1中用于进行高通滤波的滤波器的尺寸大小为(2d+1)×(2d+1)，其中d为大于等于1的整数；滤波器中所有元素值之和为0，滤波器中心元素值的绝对值大于周围元素值的绝对值。

作为优选，步骤3.1和步骤3.2中低频信息中的孤立小点为小于等于10个像素的低频信息图斑，细小孔洞为低频信息中由高频信息造成的小于等于10个像素的孔洞。

作为优选，步骤3.1和步骤3.2中的腐蚀运算和膨胀运算分别为：

腐蚀运算：

上式中，C表示低频信息提取结果，

为腐蚀运算符号，D表示形态学运算卷积核,m表示影像中像元水平方向的位置，n表示影像中像元垂直方向的位置，(D)_mn表示当卷积核D的核中心在影像中(m,n)处时卷积核D在影像中所处的位置；

膨胀运算：

上式中，C表示低频信息提取结果，

为膨胀运算符号，D表示形态学运算卷积核,m表示影像中像元水平方向的位置，n表示影像中像元垂直方向的位置，(D)_mn表示当卷积核D的核中心在影像中(m,n)处时卷积核D在影像中所处的位置；

其中，形态学运算卷积核D为

本发明的有益效果是：本发明引入高通滤波处理，使淤泥质海岸海域侧和陆域测频率信息的差异显著增大，有利于水边线位置的精确度；本发明针对遥感影像高通滤波后的像素值进行直方图统计，通过构建自适应低频信息阈值区间算法，获得了能够随影像自适应变化调整的更加稳健的阈值范围，解决了传统阈值法提取水边线中阈值范围过度依赖影像本身，阈值可迁移性差的问题；与传统边缘检测法提取水边线相比，本发明利用低频信息的边界替代图像的边缘，有效地解决了传统边缘检测法提取水边线结果破碎、不连续的缺陷。

附图说明

图1为本发明的整体流程图；

图2为本发明实施例中2020年2月22日九段沙区域Landsat8遥感影像；

图3为本发明实施例中短波红外影像数据预处理结果图；

图4为本发明实施例中影像数据高通滤波处理结果图；

图5为本发明实施例中高通滤波结果灰度直方图；

图6为本发明实施例中低频信息处理结果图；

图7为本发明实施例中水陆分离的二值图像图；

图8为本发明实施例中九段沙区域水边线结果图；

图9为阈值区间示意图。

具体实施方式

下面结合实施例对本发明做进一步描述。下述实施例的说明只是用于帮助理解本发明。应当指出，对于本技术领域的普通人员来说，在不脱离本发明原理的前提下，还可以对本发明进行若干修饰，这些改进和修饰也落入本发明权利要求的保护范围内。

淤泥质海岸由于受到海域侧水体高悬浮泥沙浓度和近岸陆地侧岸滩高含水量的影响，海域侧和陆地侧的光谱信息存在光谱混合现象，难以准确区分水体和陆地，增加了水边线遥感提取的不确定性，水边线提取难度大。本发明针对以上问题，提出一种针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法。本发明提供的针对复杂淤泥质海岸水边线遥感提取的一整套流程化方法能够获得比传统边缘检测法和阈值法更为精确且连续的水边线结果，实现了提高水边线提取精度和提取效率的目的。

实施例1：

如图1所示，一种针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法，包括以下步骤：

步骤1、对光学遥感影像数据进行几何纠正、辐射定标和大气校正等预处理操作，将原始影像像素值校正为实际地表反射率值，得到地表反射率数据；充分结合水体对于近红外波段的敏感性，光学遥感影像数据含有近红外波段，可根据实际应用需求选择包含近红外波段的光学遥感影像数据。

步骤2、低频信息提取；

步骤2.1、对地表反射率数据进行高通滤波，突出图像中的高频信息分量，抑制低频信息分量，突出高频信息、低频信息之间的差异：

上式中，m表示遥感影像中像元水平方向的位置；n表示遥感影像中像元垂直方向的位置；H(m,n)表示遥感影像中(m,n)处像元经过高通滤波之后的像素值；DN(m,n)表示遥感影像中(m,n)处像元的原始像素值；F表示高通滤波器；i表示滤波器中的元素水平方向的位置；j表示滤波器中元素垂直方向的位置；d为大于等于1的整数，可用于确定滤波器的大小；进行高通滤波的滤波器的尺寸大小为(2d+1)×(2d+1)，其中d为大于等于1的整数；滤波器中所有元素值之和为0，滤波器中心元素值的绝对值大于周围元素值的绝对值；

步骤2.2、高通滤波后像素值的直方图呈现出明显的单峰特征，单峰部分为光学遥感影像低频信息，可以利用双阈值进行提取，低频信息提取的双阈值需具有自适应性，即阈值随影像变化而自适应变化调整，避免阈值因对于影像过度依赖而造成阈值可迁移性差的缺陷，对高通滤波后的图像像素值进行直方图统计；

步骤2.2.1、将直方图划分为N个像素区间，求得直方图的间隔S：

上式中，H_max为高通滤波后像素值的最大值，H_min为高通滤波后像素值的最小值，N为像素区间个数；

步骤2.2.2、统计各像素区间的占比，按降序排列；

步骤2.2.3、筛选出百分比最高的前B个像素区间，确定直方图中百分比最高的前1％个像素区间；B的值由如下公式确定：

上式中，

为向下取整符号，A为直方图中百分比最高的前1％个像素区间；

步骤2.3、计算得到低频信息阈值区间，构建自适应的阈值区间算法，根据阈值区间自适应提取低频信息，低频信息值为1，高频信息值为0；

步骤2.3.1、假设前B个像素区间中最小像素值为H_{B_min}，则低频信息阈值区间表示为如下公式：

[H_{B_min},H_{B_min}+B×S]

上式中，H_{B_min}为前B个像素区间中最小像素值；B为百分比最高的像素区间个数，S为直方图的间隔；

步骤2.3.2、构建自适应的阈值区间算法为：

(H ge H_{B_min})and(H le(H_{B_min}+B×S))

上式中，H表示高通滤波后的像素值，ge表示大于等于，le表示小于等于，H_max为高通滤波后像素值的最大值，H_min为高通滤波后像素值的最小值；B为百分比最高的像素区间个数，S为直方图的间隔；

步骤3、对低频信息进行形态学处理，处理完成后得到水陆分离的二值图像；

步骤3.1、利用形态学开运算填充低频信息中的孤立小点(小于等于10个像素的低频信息图斑)；所述形态学开运算为先进行腐蚀运算，再进行膨胀运算；

步骤3.2、利用形态学闭运算去除低频信息中的细小孔洞(由高频信息造成的小于等于10个像素的孔洞)；所述形态学闭运算为先进行膨胀运算，再进行采用腐蚀运算填充细小孔洞；

步骤3.3、利用膨胀运算弥补低频信息边界与影像实际水边线之间的像素误差；

其中：

腐蚀运算：

上式中，C表示低频信息提取结果，

为腐蚀运算符号，D表示形态学运算卷积核,m表示影像中像元水平方向的位置，n表示影像中像元垂直方向的位置，(D)_mn表示当卷积核D的核中心在影像中(m,n)处时卷积核D在影像中所处的位置；

膨胀运算：

上式中，C表示低频信息提取结果，

为膨胀运算符号，D表示形态学运算卷积核,m表示影像中像元水平方向的位置，n表示影像中像元垂直方向的位置，(D)_mn表示当卷积核D的核中心在影像中(m,n)处时卷积核D在影像中所处的位置；

其中，形态学运算卷积核D为

步骤4、对水陆分离的二值图像进行后处理，后处理操作包括但不限于按面积筛除、人工目视解译等方法，可根据影像实际处理结果按需选择，最终得到水边线结果：将水陆分离的二值图像转为矢量格式，按照矢量格式的多边形面积及属性进行筛选：

(area≥0.05)and(GRIDCODE＝1)

上式中，area表示多边形面积，单位为平方千米，GRIDCODE代表多边形属性，取值为1表示水体；

采用人机交互解译进行人工检核，确定最终水边线。

实施例2：

将实施例1的针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法应用于九段沙区域；九段沙是位于长江口外南北槽之间的拦门沙河段，是典型的淤泥质海岸区域；使用九段沙区域2020年2月22日Landsat8 OLI遥感影像进行瞬时水边线提取。

(1)所采用的遥感影像数据为光学遥感影像数据，所使用传感器为Landsat8卫星所搭载的OLI传感器，获得如图2所示2020年2月22日九段沙区域Landsat8遥感影像，影像云量为0.01％；

(2)对光学遥感影像数据进行图像预处理，包括几何纠正、辐射定标及大气校正，得到地表反射率数据，为便于存储和运算，将地表反射率数据扩大10000倍；获得如图3所示对短波红外影像数据的预处理结果；

(3)本实施例中取d＝1，确定高通滤波器大小为3×3，所使用的高通滤波器表示如下:

对经过预处理的光学遥感影像数据采用高通滤波法，突出图像中的高频分量，抑制低频分量，得到如图4所示影像数据的高通滤波处理结果；

(4)取N＝256，将直方图划分为256个像素区间，高通滤波后的最大像素值H_max和最小像素值H_min分别为21886和-11461，代入公式

可求得直方图的间隔S＝131，按照上述参数对光学遥感影像数据的高通滤波结果进行直方图统计，得到如图5所示高通滤波结果的灰度直方图；

(5)确定直方图中百分比最高的前B个像素区间，由公式

计算得出A的结果为2，当A为偶数时，B＝A+1，因此需确定直方图中百分比最高的前3个像素区间，这3个像素区间中的最小像素值为-195，因此根据自适应的阈值区间公式[H_{B_min},H_{B_min}+B×S]，对应图9的阈值区间示意图，可计算出低频信息的阈值区间为[-195,198]，按照阈值区间算法进行低频信息提取，提取公式为：

(H ge-195)and(H le 198)

提取得到如图6所示低频信息；图9中H_min为高通滤波后的最小像素值，H_{B_min}为前B个像素区间中最小像素值；N为像素区间个数；S为直方图的间隔；

(6)对低频信息进行形态学处理，先后经过开运算、闭运算和膨胀运算，得到如图7所示的水陆分离的二值图像；

(7)将水陆分离的二值图像转为矢量面进行编辑，并进行人工检核，得到如图8所示最终水边线结果。

Claims (8)

1.一种针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法，其特征在于，包括以下步骤：

步骤1、对遥感影像数据进行预处理操作，将原始影像像素值校正为实际地表反射率值，得到地表反射率数据；

步骤2、低频信息提取；

步骤2.1、对地表反射率数据进行高通滤波：

上式中，m表示遥感影像中像元水平方向的位置；n表示遥感影像中像元垂直方向的位置；H(m,n)表示遥感影像中(m,n)处像元经过高通滤波之后的像素值；DN(m,n)表示遥感影像中(m,n)处像元的原始像素值；F表示高通滤波器；i表示滤波器中的元素水平方向的位置；j表示滤波器中元素垂直方向的位置；d为大于等于1的整数；

步骤2.2、对高通滤波后的图像像素值进行直方图统计；

步骤2.3、计算得到低频信息阈值区间，构建自适应的阈值区间算法，根据阈值区间自适应提取低频信息，低频信息值为1，高频信息值为0；

步骤3、对低频信息进行形态学处理，处理完成后得到水陆分离的二值图像；

步骤3.1、利用形态学开运算填充低频信息中的孤立小点；所述形态学开运算为先进行腐蚀运算，再进行膨胀运算；

步骤3.2、利用形态学闭运算去除低频信息中的细小孔洞；所述形态学闭运算为先进行膨胀运算，再进行采用腐蚀运算填充细小孔洞；

步骤3.3、利用膨胀运算弥补低频信息边界与影像实际水边线之间的像素误差；

步骤4、对水陆分离的二值图像进行后处理：将水陆分离的二值图像转为矢量格式，按照矢量格式的多边形面积及属性进行筛选：

(area≥0.05)and(GRIDCODE＝1)

上式中，area表示多边形面积，单位为平方千米，GRIDCODE代表多边形属性，取值为1表示水体；

采用人机交互解译进行人工检核，确定最终水边线。

2.根据权利要求1所述针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法，其特征在于，步骤2.2具体包括以下步骤：

步骤2.2.1、将直方图划分为N个像素区间，求得直方图的间隔S：

上式中，H_max为高通滤波后像素值的最大值，H_min为高通滤波后像素值的最小值，N为像素区间个数；

步骤2.2.2、统计各像素区间的占比，按降序排列；

步骤2.2.3、筛选出百分比最高的前B个像素区间，确定直方图中百分比最高的前1％个像素区间；B的值由如下公式确定：

上式中，

为向下取整符号，A为直方图中百分比最高的前1％个像素区间。

3.根据权利要求1所述针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法，其特征在于，步骤2.3具体包括以下步骤：

步骤2.3.1、假设前B个像素区间中最小像素值为H_{B_min}，则低频信息阈值区间表示为如下公式：

[H_{B_min},H_{B_min}+B×S]

上式中，H_{B_min}为前B个像素区间中最小像素值；B为百分比最高的像素区间个数，S为直方图的间隔；

步骤2.3.2、构建自适应的阈值区间算法为：

(H ge H_{B_min})and(H le(H_{B_min}+B×S))

上式中，H表示高通滤波后的像素值，ge表示大于等于，le表示小于等于，H_max为高通滤波后像素值的最大值，H_min为高通滤波后像素值的最小值；B为百分比最高的像素区间个数，S为直方图的间隔。

4.根据权利要求1所述针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法，其特征在于：步骤1中遥感影像数据为光学遥感影像数据，光学遥感影像数据含有近红外波段。

5.根据权利要求1所述针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法，其特征在于：步骤1中对遥感影像数据进行预处理的操作包括几何纠正、辐射定标和大气校正。

6.根据权利要求1所述针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法，其特征在于：步骤2.1中用于进行高通滤波的滤波器的尺寸大小为(2d+1)×(2d+1)，其中d为大于等于1的整数；滤波器中所有元素值之和为0，滤波器中心元素值的绝对值大于周围元素值的绝对值。

7.根据权利要求1所述针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法，其特征在于：步骤3.1和步骤3.2中低频信息中的孤立小点为小于等于10个像素的低频信息图斑，细小孔洞为低频信息中由高频信息造成的小于等于10个像素的孔洞。

8.根据权利要求1所述针对近岸复杂环境下淤泥质海岸的水边线遥感提取方法，其特征在于，步骤3.1和步骤3.2中的腐蚀运算和膨胀运算分别为：

腐蚀运算：

上式中，C表示低频信息提取结果，

为腐蚀运算符号，D表示形态学运算卷积核,m表示影像中像元水平方向的位置，n表示影像中像元垂直方向的位置，(D)_mn表示当卷积核D的核中心在影像中(m,n)处时卷积核D在影像中所处的位置；

膨胀运算：

上式中，C表示低频信息提取结果，

为膨胀运算符号，D表示形态学运算卷积核,m表示影像中像元水平方向的位置，n表示影像中像元垂直方向的位置，(D)_mn表示当卷积核D的核中心在影像中(m,n)处时卷积核D在影像中所处的位置；

其中，形态学运算卷积核D为

Images