CN111881816B - 一种长时序的河湖围埂养殖区域监测方法 - Google Patents

Abstract

本发明涉及一种长时序的河湖围埂养殖区域监测方法，先根据需求获取河湖遥感影像，在对遥感影像统一影像分辨率的基础上，计算遥感影像的云含量，云含量小于5％的影像为有效观察，同时构建其云模板Cloud_mask；采用迭代的阈值分割方法获取水体斑块，计算水体监测频率分布图；接着在河湖管理范围矢量裁剪的基础上，对河湖水体监测频率分布图构建水体监测斑块，根据河湖周长、面积特征，将水体监测斑块分为河湖主体和周边鱼塘；提取河湖主体水体频率分布图中频率变化大的区域为围埂区域，周边鱼塘为养殖区域；将提取的河湖围埂和养殖区域进行矢量化并编号，添加属性信息。该方法能够快速准确的识别河湖围埂养殖区域，准确率高。

Description

一种长时序的河湖围埂养殖区域监测方法

技术领域

本发明涉及一种长时序的河湖围埂养殖区域监测方法，属于水利、河湖水域岸线监管技术领域。

背景技术

随着水产品市场需求大幅增加和社会经济快速发展，河湖围埂用地、养殖用地面积迅速增长，大面积的河湖沼泽等已被围埂、养殖用地取代，加剧了沿湖生态系统的脆弱性。围埂养殖用地评估对做好养殖污染防治及围埂围网拆除工作，保护生态环境至关重要。对湖库围埂养殖区域进行监控分析，也有利于加强政府对河湖周边地区线管基础数据的掌握，为科学规划并实现对湖库资源的可持续发展提供重要的基础信息支撑。

传统的现场调查手段受限于人力物力，很难满足大范围的围埂养殖区监测及实时信息更新的需求。遥感影像因其覆盖范围面积大，能真实记录和反映地表信息，获取周期短等特点，已成为土地利用变化监测、资源管理等方面的重要数据源。将遥感技术应用于湖库围埂养殖的自动提取和变化识别，相比传统的实地调查统计方法具有实时性强、准确度高、空间可表达等优势。目前对围埂养殖区域的提取方法有目视解译法、面向对象的信息提取方法等，前者需要依靠人工解译，后者主要适用于近海养殖区域。河湖围埂养殖区域在影像上所表现的复杂性和不均匀性，造成该目标识别中的困难。利用遥感影像监测河湖周边围埂养殖用地多为采用两个时相高分辨率遥感影像进行围埂养殖区域监测，由于围埂相对于湖库水面高差小，受水位波动影响大，容易出现误判，往往需要现场调研进行判断。受卫星过境重返时间限制，难以获取合适时间的高分辨率遥感影像，且经济成本较高。此外，单一时相中等分辨率遥感影像识别结果难以区分水稻田和围埂养殖鱼塘。

对此，本发明提出了一种适用于中等分辨率遥感影像长时序监测的围埂养殖区域监测方法，可快速准确的识别河湖围埂养殖区域。

发明内容

本发明针对现有技术存在的不足，提供了一种长时序的河湖围埂养殖区域监测方法，能够快速准确的识别河湖围埂养殖区域。该方法主要包括4个步骤：首先根据需求获取河湖遥感影像，在对遥感影像进行正射校正、几何配准和重采样统一影像分辨率的基础上，计算遥感影像的云含量，云含量(即遥感影像中云覆盖影响区域面积占影像总面积的比例)小于5％的影像为有效观察，同时构建其云模板Cloud_mask(其中“1”表示云覆盖影响区域，“0”表示无云覆盖影响区域)；然后，采用迭代的阈值分割方法获取水体斑块，计算水体监测频率分布图；接着在河湖管理范围矢量裁剪的基础上，对河湖水体监测频率分布图构建水体监测斑块，根据河湖周长、面积特征，将水体监测斑块分为两部分：河湖主体和周边鱼塘；提取河湖主体水体频率分布图中频率变化大的区域为围埂区域，周边鱼塘为养殖区域；最后，将提取的河湖围埂和养殖区域进行矢量化并编号，添加属性信息。具体技术方案如下：

一种长时序的河湖围埂养殖区域监测方法，包括以下步骤：

步骤一、获取河湖遥感影像；

步骤二、对不同类型的遥感影像进行处理，然后重采样统一不同类型遥感影像的空间分辨率；

步骤三、对不同类型的遥感影像构建Cloud_mask；

步骤四、在构建Cloud_mask的基础上，筛选出云含量小于5％的遥感影像作为有效观察；

步骤五、取每景遥感影像的近红外波段和绿色波段，计算每一景影像的归一化水体指数NDWI值；其中，Band2表示绿光波段的DN值，Band4表示近红外波段的DN值，NDWI值的范围为-1～1；

NDWI＝(Band4-Band2)/(Band2+Band4)；

NDWI_new＝(NDWI+1)*127.5；

其中，NDWI_new表示修改后的归一化水体指数；

NDWI_new取值范围为0～255，并去除归一化水体指数NDWI_new中的云覆盖影像区域，即令Cloud_mask覆盖像素值为0；

步骤六、采用迭代的最大类间方差阈值分割方法获取水体斑块；

步骤七、在对所有有效观察遥感影像提取水体的基础上，计算水体监测频率分布图；

步骤八、在河湖管理范围矢量裁剪的基础上，对河湖水体监测频率分布图构建水体监测斑块；

步骤九、根据河湖周长、面积特征，将水体监测斑块分为河湖主体和周边鱼塘；将水体频率大于0.1的水体构建水体频率斑块，并对其进行编号，计算每一个斑块的面积，在去除湖面面积的基础上，得到湖泊周边离散的养殖区域，将其存储为带坐标“.GIF”文件；

步骤十、对河湖主体水体频率分布图计算3邻域标准差，去除最外层轮廓线，得到湖体内部的标准差较大的区域为围埂区域，将其存储为带坐标“.GIF”文件；

步骤十一、将提取的河湖围埂和养殖区域文件进行矢量化，并对边缘进行简化，得到面图层文件并编号，计算面积，中心点位置的属性信息。

上述技术方案的进一步优化，在步骤二中，对不同类型的遥感影像进行处理的方式包括正射校正、几何配准。

上述技术方案的进一步优化，在步骤三中，对于通过哨兵二号卫星或Landsat4-8卫星获取的遥感影像和数据文件确定云含量，并构建Cloud_mask，Cloud_mask中像素值为1表示云，像素值为0表示无云；

对于通过高分一号卫星、高分六号卫星、HJ-1A/1B卫星获取的遥感数据，利用波段运算，取阈值分割得到厚云范围，然后采用结构元素B对厚云范围进行膨胀，构建Cloud_mask，Cloud_mask中像素值为1表示云，像素值为0表示无云；其中，Cloud表示厚云范围，T表示阈值；

Cloud＝(Band1+Band2+Band3+Band4)>T；

Cloud_mask＝Cloud⊕B；

其中，Band1表示蓝光波段的DN值，Band2表示绿光波段的DN值，Band3表示红光波段的DN值，Band4表示近红外波段的DN值；

上述技术方案的进一步优化，采用迭代的最大类间方差阈值分割方法获取水体斑块的方法包括如下步骤：

1)、在取初始分割阈值为128的基础上，进行水体和陆地二值分割，得到初始水体斑块；

2、构建水体斑块缓冲区，使得缓冲区的面积为水体斑块的两倍；

3、采用最大类间方差法计算缓冲区NDWI_new影像的分割阈值，再次进行水体和陆地二值分割，直到前后两次分割，水体面积一致。

上述技术方案的进一步优化，在步骤七中，计算水体监测频率分布图的方法为：

其中，Fre为每一个像素点的水体有效观察频率，m为每一个像素观测为水体的总次数，M为每一个像素的总有效观察次数，EO为所有像素每一次有效观察组成的矩阵，W为观测为水体所有像素组成的矩阵；i为正整数，取值为1,2…m；j为正整数，取值为1,2…M。

本发明的有益效果：

所述长时序的河湖围埂养殖区域监测方法极大的减少人工解译工作量，能够快速准确的识别河湖围埂养殖区域，实验湖泊的准确率超过90％，实施效果好，应用价值高，具有极大的推广价值。

附图说明

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

图2为去除云后的NDWI影像图；

图3为初始水体和陆地二值分割结果图；

图4为水体斑块缓冲区；

图5为缓冲区NDWI影像图；

图6为河湖水体监测频率分布图；

图7为养殖区域监测结果图；

图8为河湖主体计算像素3邻域标准差结果图；

图9为围埂区域监测结果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例1

如图1所示，所述长时序的河湖围埂养殖区域监测方法的具体流程如下：

1、首先根据需求获取河湖遥感影像。包括但不限于Sentinel-2卫星地表反射率数据、Landsat系列卫星(Landsat4卫星、Landsat5卫星、Landsat7卫星、Landsat8卫星)地表反射率数据、高分一号卫星(GF-1)1级数据、高分六号卫星(GF-6)1级数据、HJ-1A/1B卫星2级数据等。

2、判断卫星影像数据类型，针对不同类型的卫星影像数据进行不同的处理：对16米分辨率的GF-1卫星、GF-6卫星影像数据进行辐射定标、正射校正、几何配准；对30米分辨率的HJ-1A/1B卫星影像数据进行几何配准；对Landsat7卫星影像数据进行条带修复。

因为不同类型的卫星影像空间分辨率不同，为综合利用多源卫星数据，最后重采样统一不同类型影像的空间分辨率。

3、判断输入卫星影像数据是否均是同一类型(如Sentinel-2卫星、Landsat系列卫星、GF-1卫星、GF-6卫星、HJ-1A/1B卫星)，如果卫星影像数据分辨率不一致，以最高分辨率影像数据为基准，将相对较低分辨率卫星遥感数据以最邻近插值法重采样成最高分辨率影像数据。

4、对卫星影像数据进行云检测，构建Cloud_mask。针对不同类型的卫星影像数据进行不同的处理：对于Landsat系列卫星遥感影像根据卫星提供的“pixel_qa.GIF”确定云含量，并构建Cloud_mask(像素值为1表示云，0表示无云)；对于Sentinel-2卫星获取的数据也根据卫星提供的文件确定云含量，并构建Cloud_mask(像素值为1表示云，0表示无云)。

5、由于高分一号、高分六号、HJ-1A/1B没有提供可确定云含量的文件，考虑到厚云在高分一号、高分六号、HJ-1A/1B遥感影像中各波段的像素值均较高(GF-1不小于650，GF-1不小于2000，HJ-1A/1B不小于255)以及对运算速度的需要和精度要求不高，根据下式阈值分割得到厚云范围Cloud，并去除零星散点：

Cloud_{GF1、HJ-1A/1B}＝(Band1+Band2+Band3+Band4)/4>T1，

Cloud_GF6＝Band7>T2∩(Band1+Band2+Band3+Band4)/4>T1,

其中，Band1表示蓝光波段的DN值，Band2表示绿光波段的DN值，Band3表示红光波段的DN值，Band4表示近红外波段的DN值，Band7表示紫波段的DN值。T1、T2为阈值，可人工选取合适的阈值(如T1、T2取值GF-1时不小于650，取值GF-1时不小于2000，取值HJ-1A/1B时不小于255)，也可采用最大类间方差法自动计算。

6、考虑到云阴影，采用结构元素B对厚云范围Cloud进行膨胀，构建Cloud_mask(像素值为1表示云，0表示无云)。

Cloud_mask＝Cloud⊕B；

其中，结构元素B是遥感数字图像处理，数学形态学中的结构元素矩阵。

7、在构建Cloud_mask的基础上，筛选出影像清晰，云含量Area_cloud小于5％的遥感影像作为有效观察。

8、取每景遥感影像的近红外波段和绿色波段，计算每一景影像的NDWI值(范围为-1～1)；

NDWI＝(Band4-Band2)/(Band2+Band4)；

其中，Band4为近红外波段反射率值(或DN值)，Band2为绿光波段反射率值(或DN值)。

采用下式将NDWI的取值(-1～1)修改为NDWI_new(范围：0～255)，同时令Cloud_mask覆盖像素值为0(如图2所示)。

NDWI_new＝(NDWI+1)*127.5；

9、采用迭代的阈值分割方法获取水体斑块：

1)、在取初始分割阈值为128的基础上，进行水体和陆地二值分割，水体为1，陆地为0，得到初始水体斑块(如图3所示)；

2)、构建水体斑块缓冲区，使得缓冲区的面积为水体斑块的两倍(如图4所示)；

为了快速计算水体和非水体二值分割的最佳阈值，降低干扰，通过上一次分割水体的结果外扩一倍的范围构建缓冲区，下一次只计算缓冲区内NDWInew影像的分割阈值。以此迭代，直到前后两次分割，水体面积一致；此时的阈值最接近最佳分割阈值，而采用这个阈值得到水体分割结果也是最好的。

3)、计算分割阈值，可以采用最大类间方差法自动计算缓冲区NDWI_new影像(如图5所示)的分割阈值，也可采用直方图反复试验目视解译；

4)、根据步骤3确定阈值，再次进行水体和陆地二值分割，直到前后两次分割，水体面积一致。

10、在对所有有效观察遥感影像提取水体的基础上，计算水体监测频率分布图；

其中，Fre为每一个像素点的水体有效观察频率，m为每一个像素观测为水体的总次数，M为每一个像素的总有效观察次数，EO为所有像素每一次有效观察组成的矩阵，W为观测为水体所有像素组成的矩阵；i为正整数，取值为1,2…m；j为正整数，取值为1,2…M。

11、在河湖管理范围矢量裁剪的基础上，对河湖水体监测频率分布图构建水体监测斑块(如图6所示)。

12、根据河湖周长、面积特征，将水体监测斑块分为两部分：河湖主体和周边鱼塘；将水体频率大于0.1的水体构建水体频率斑块，并对其进行编号，计算每一个斑块的面积，在去除湖面面积的基础上，得到湖泊周边离散的养殖区域(如图7所示)，将其存储为带坐标“.GIF”文件。

13、对河湖主体水体频率分布图计算3邻域标准差(如图8所示)，去除最外层轮廓线，得到湖体内部的标准差较大(即频率变化大)的区域为围埂区域(如图9所示)，将其存储为带坐标“.GIF”文件。

14、将提取的河湖围埂和养殖区域文件进行矢量化，并对边缘进行简化，得到面图层文件并编号，计算面积，中心点位置等属性信息。

在上述实施例中，在本发明中，“长时序”顾名思义是长时间序列，通俗的将是采用很多个时间获取的遥感影像进行监测，例如实验湖泊，本发明采用了2017年4月-2020年4月卫星过境获取的上百景所有卫星遥感影像。

传统的现场调查手段受限于人力物力，很难满足大范围的围埂养殖区监测及实时信息更新的需求。遥感影像因其覆盖范围面积大，能真实记录和反映地表信息，获取周期短等特点，已成为土地利用变化监测、资源管理等方面的重要数据源。将遥感技术应用于湖库围埂养殖的自动提取和变化识别，相比传统的实地调查统计方法具有实时性强、准确度高、空间可表达等优势。目前对围埂养殖区域的提取方法有目视解译法、面向对象的信息提取方法等，前者需要依靠人工解译，后者主要适用于近海养殖区域。河湖围埂养殖区域在影像上所表现的复杂性和不均匀性，造成该目标识别中的困难。利用遥感影像监测河湖周边围埂养殖用地多为采用两个时相高分辨率遥感影像进行围埂养殖区域监测，由于围埂相对于湖库水面高差小，受水位波动影响大，容易出现误判，往往需要现场调研进行判断。受卫星过境重返时间限制，难以获取合适时间的高分辨率遥感影像，且经济成本较高。此外，单一时相中等分辨率遥感影像识别结果难以区分水稻田和围埂养殖鱼塘。

对此，本发明提出了一种适用于中等分辨率遥感影像的长时序的河湖围埂养殖区域监测方法，极大的减少人工解译工作量，能够快速准确的识别河湖围埂养殖区域，经过大量的试验，实验湖泊的准确率超过90％，即提取的围埂、养殖区域的点位经过验证超过90％是对的。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种长时序的河湖围埂养殖区域监测方法，其特征在于包括以下步骤：

步骤一、获取河湖遥感影像；

步骤二、对不同类型的遥感影像进行处理，然后重采样统一不同类型遥感影像的空间分辨率；

步骤三、对不同类型的遥感影像构建Cloud_mask；

步骤四、在构建Cloud_mask的基础上，筛选出云含量小于5％的遥感影像作为有效观察；

步骤五、取每景遥感影像的近红外波段和绿色波段，计算每一景影像的归一化水体指数NDWI值；其中，Band2表示绿光波段的DN值，Band4表示近红外波段的DN值，NDWI值的范围为-1～1；

NDWI＝(Band4-Band2)/(Band2+Band4)；

NDWI_new＝(NDWI+1)*127.5；

其中，NDWI_new表示归一化水体指数；

NDWI_new取值范围为0～255，并去除归一化水体指数NDWI_new中的云覆盖影像区域，即令Cloud_mask覆盖像素值为0；

步骤六、采用迭代的最大类间方差阈值分割方法获取水体斑块；

步骤七、在对所有有效观察遥感影像提取水体的基础上，计算水体监测频率分布图；

步骤八、在河湖管理范围矢量裁剪的基础上，对河湖水体监测频率分布图构建水体监测斑块；

步骤九、根据河湖周长、面积特征，将水体监测斑块分为河湖主体和周边鱼塘；将水体频率大于0.1的水体构建水体频率斑块，并对其进行编号，计算每一个斑块的面积，在去除湖面面积的基础上，得到湖泊周边离散的养殖区域，将其存储为带坐标“.GIF”文件；

步骤十、对河湖主体水体频率分布图计算3邻域标准差，去除最外层轮廓线，得到湖体内部的标准差较大的区域为围埂区域，将其存储为带坐标“.GIF”文件；

步骤十一、将提取的河湖围埂和养殖区域文件进行矢量化，并对边缘进行简化，得到面图层文件并编号，计算面积，中心点位置的属性信息。

2.根据权利要求1所述的一种长时序的河湖围埂养殖区域监测方法，其特征在于：在步骤二中，对不同类型的遥感影像进行处理的方式包括正射校正、几何配准。

3.根据权利要求1所述的一种长时序的河湖围埂养殖区域监测方法，其特征在于：在步骤三中，对于通过哨兵二号卫星或Landsat4-8卫星获取的遥感影像和数据文件确定云含量，并构建Cloud_mask，Cloud_mask中像素值为1表示云，像素值为0表示无云；

对于通过高分一号卫星、高分六号卫星、HJ-1A/1B卫星获取的遥感数据，利用波段运算，取阈值分割得到厚云范围，然后采用结构元素B对厚云范围进行膨胀，构建Cloud_mask，Cloud_mask中像素值为1表示云，像素值为0表示无云；其中，Cloud表示厚云范围，T表示阈值；

Cloud＝(Band1+Band2+Band3+Band4)>T；

Cloud_mask＝Cloud⊕B；

其中，Band1表示蓝光波段的DN值，Band2表示绿光波段的DN值，Band3表示红光波段的DN值，Band4表示近红外波段的DN值；

4.根据权利要求1所述的一种长时序的河湖围埂养殖区域监测方法，其特征在于：采用迭代的最大类间方差阈值分割方法获取水体斑块的方法包括如下步骤：

1)、在取初始分割阈值为128的基础上，进行水体和陆地二值分割，得到初始水体斑块；

2)、构建水体斑块缓冲区，使得缓冲区的面积为水体斑块的两倍；

3)、采用最大类间方差法计算缓冲区NDWI_new影像的分割阈值，再次进行水体和陆地二值分割，直到前后两次分割，水体面积一致。

5.根据权利要求1所述的一种长时序的河湖围埂养殖区域监测方法，其特征在于：在步骤七中，计算水体监测频率分布图的方法为：

其中，Fre为每一个像素点的水体有效观察频率，m为每一个像素观测为水体的总次数，M为每一个像素的总有效观察次数，EO为所有像素每一次有效观察组成的矩阵，W为观测为水体所有像素组成的矩阵；i为正整数，取值为1,2…m；j为正整数，取值为1,2…M。

Images