CN113177473B - 遥感影像水体自动化提取方法和装置 - Google Patents
遥感影像水体自动化提取方法和装置 Download PDFInfo
- Publication number
- CN113177473B CN113177473B CN202110471313.5A CN202110471313A CN113177473B CN 113177473 B CN113177473 B CN 113177473B CN 202110471313 A CN202110471313 A CN 202110471313A CN 113177473 B CN113177473 B CN 113177473B
- Authority
- CN
- China
- Prior art keywords
- water body
- remote sensing
- pixel
- sensing image
- index
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/13—Satellite images
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/168—Segmentation; Edge detection involving transform domain methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/50—Depth or shape recovery
- G06T7/507—Depth or shape recovery from shading
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/10—Image acquisition
- G06V10/12—Details of acquisition arrangements; Constructional details thereof
- G06V10/14—Optical characteristics of the device performing the acquisition or on the illumination arrangements
- G06V10/143—Sensing or illuminating at different wavelengths
Abstract
本发明公开了一种遥感影像水体自动化提取方法和装置，属于水体监测领域。该方法首先获取卫星遥感影像及成像参数并进行预处理，其次基于影像的光谱、梯度及统计特征等，构建水体NDWI、MSWI、ESWI、Yellow、Canny等特征信息；利用NDWI和p(NIR)特征识别纯净水体，然后耦合特征信息，综合去除云层、建筑和山体阴影的错分；构造黄度特征以及特征水体样本补充高含沙量水体及其他异常水体；最后通过边界处理及最小图斑调整，获得普适水体的空间分布。本发明解决了因建筑、山体等阴影以及云层导致的水体错分问题，并解决了对高泥沙含量水体、富营养化等异常水体的漏分问题，克服传统单一阈值、复杂地物背景下水体信息提取的局限性。
Description
技术领域
本发明涉及水体监测领域，特别是指一种遥感影像水体自动化提取方法和装置。
背景技术
水体信息在水利规划，水资源调查、保护与监测，洪涝灾害应急救助，水储量动态变化等方面是重要的基础地理信息数据，掌握水体信息的覆盖范围在水资源的动态监测，在保证水源安全、合理规划及保护等方面发挥着巨大的作用。
遥感在大尺度空间范围的水体信息提取方面具有强大的优势，如短周期观测、覆盖范围广、空间分辨率及尺度信息丰富、光谱信息量大等，其获取的水体信息在全国湖泊水域及水质量监测、城市黑臭水体巡查、入海河口浑浊水体监测、地表水域与地下水影响机制研究、洪水淹没范围及面积等方面是基础性数据。
水体在遥感卫星影像上的背景信息简单，特征明显且易分辨。常规获取水体分布信息的方法主要包括两类，一类是采用人工目视解译方法，解译水体特征，获取水域空间分布，手动勾绘水域边界，得到水体空间矢量。另外一类是采用光谱特征及其差异并结合人工解译的半自动化提取方法，如单波段阈值法、谱间关系法及图像分类法等。单波段阈值法为采用单一波段阈值分割后提取；谱间关系法为采用多个波段，借助差值、比值等模型，每个波段采用单一阈值分割提取；图像分类为借助分类器将影像分为水体和非水体对象，将水体与非水体信息分离，取水体分类对象边界作为水域的分布范围。
采用人工目视解译方法提取水体的方法，优点是在水体背景信息简单 (或其他地物对水体边界判断不造成干扰)、范围较小、水体分布较为集中的条件下，提取速度快，精度相对有保证，无需复杂操作，简便易行。该方法的缺点和不足也显而易见，在水域分布较为分散、水域特征较为复杂的条件下，勾绘工作量相对较大、目视识别相对较为困难，容易有遗漏、边界误判的错误。同时不同的专业人员容易造成边界界定不统一，存在标准不一导致的水域边界错判误差。
采用光谱特征及其差异并结合人工解译的半自动化提取方法主要包括单波段阈值法、谱间关系法、图像分类法及其组合法等，采用不同光谱对比分析增强水体特征信息，将水体与非水体进行区分，该方法容易受影像的成像参数(包括卫星传感器的光谱范围及波段数)及影像信噪比的影响；其次，容易导致灰度值与水体接近的地物误提，在水体中混淆非水体信息，不合理的阈值容易造成水域边界不拟合或过拟合的极端现象，同时阈值因影像的不同而不同，普适性较低；此外，水体与陆地交界区域，水动力、水环境、水底质等条件受区域性影响较大，特别是湿地滩涂、浅滩等区域，该方法对高泥沙含量水体、富营养水体的识别度不高；同时未考虑建筑物阴影、山体阴影以及云层的影响，导致部分区域水体提取精度不高。
发明内容
为解决上述技术问题，本发明提供一种遥感影像水体自动化提取方法和装置，本发明解决了因建筑、山体等阴影以及云层导致的水体错分问题，并解决了对高泥沙含量水体、富营养化等异常水体的漏分问题，克服传统单一阈值、复杂地物背景下水体信息提取的局限性。
本发明提供技术方案如下：
一种遥感影像水体自动化提取方法，所述方法包括：
S1：获取遥感影像并对所述遥感影像进行预处理；其中，所述遥感影像包括红波段、绿波段、蓝波段和近红外波段的光谱信息；
S2：利用遥感影像的归一化水体指数以及近红外波段的光谱信息提取纯净水体；
S3：利用遥感影像的蓝波段的光谱信息提取云层信息，并去除云层信息对纯净水体的干扰；
S4：将遥感影像变换到HIS空间，利用遥感影像的HIS信息、归一化水体指数、改进的阴影水体指数以及增强阴影水体指数从纯净水体中剔除山体阴影；
S5：根据纯净水体中各个图斑的边界强度从纯净水体中剔除建筑阴影；
S6：利用遥感影像的黄度指数提取高泥沙含量水体；
S7：选取异常水体特征样本，根据随机森林分类法分离提取异常水体；
S8：将纯净水体、高泥沙含量水体和异常水体合并，并进行修正，得到水体空间分布和面积。
进一步的，所述S2包括：
S21：计算遥感影像的各个像元的归一化水体指数；
NDWI＝(p(Green)-p(NIR))/(p(Green)+p(NIR))；
其中，NDWI为归一化水体指数，p(Green)和p(NIR)分别为遥感影像中各个像元在绿波段和近红外波段的像元值；
S22：将遥感影像的各个像元的归一化水体指数与设置的纯净水体阈值进行比较，将归一化水体指数大于纯净水体阈值的像元划分为纯净水体；
S23：将遥感影像中各个像元在近红外波段的像元值按照从低到高排序，将从低到高第90％位置的像元值作为补充水体阈值；
S24：将遥感影像中各个像元在近红外波段的像元值与所述补充水体阈值进行比较，将近红外波段的像元值小于补充水体阈值的像元划分为纯净水体；
S25：将所述纯净水体阈值缩小一定数值，将遥感影像的各个像元的归一化水体指数与缩小的纯净水体阈值进行比较，将归一化水体指数大于缩小的纯净水体阈值的像元划分为纯净水体。
进一步的，所述S3包括：
S31：设置高亮地物的值域范围[p1,p2]；
S32：获取遥感影像中高亮地物的各个像元在蓝波段的像元值的平均值p；
S33：将所述平均值P与p1、p2进行比较，若所述平均值P小于p1，则将P1作为云层分割阈值，若所述平均值P大于p2，则将P2作为云层分割阈值，若所述平均值P∈[p1,p2]，则将所述平均值P作为云层分割阈值；
S34：将遥感影像中各个像元在蓝波段的像元值与所述云层分割阈值进行比较，将蓝波段的像元值大于云层分割阈值的像元划分为云层；
S35：将包含在所述纯净水体内部且大小小于200个像元的图斑划分为纯净水体；
S36：将包含在所述纯净水体内部的云层划分为纯净水体。
进一步的，所述S4包括：
S41：将所述遥感影像变换到HIS空间，将H＜0.085，并且I＜0.115，并且S＜0.441，并且NDWI＜0.35的像元划分为山体阴影；
S42：计算遥感影像的各个像元的改进的阴影水体指数；
MSWI＝(p(Blue)-p(NIR))/(p(NIR))；
其中，MSWI为改进的阴影水体指数，p(Blue)为遥感影像中各个像元在蓝波段的像元值；
S43：将遥感影像的各个像元的改进的阴影水体指数与设置的第一阴影阈值进行比较，将改进的阴影水体指数小于第一阴影阈值的像元划分为山体阴影；
S44：计算遥感影像的各个像元的增强阴影水体指数；
ESWI＝(p(Blue)-p(Green))/(p(NIR)+p(NIR))；
其中，ESWI为增强阴影水体指数；
S45：将遥感影像的各个像元的增强阴影水体指数与设置的第二阴影阈值进行比较，将增强阴影水体指数小于第二阴影阈值的像元划分为山体阴影；
S46：将所述山体阴影从所述纯净水体中剔除。
进一步的，所述S5包括：
S51：计算纯净水体中各个图斑的边界指数特征border index；
其中，bv为图斑的边界周长，Iv为图斑的长度，Wv为图斑的宽度；
S52：将边界指数特征大于4并且面积小于1500个像元的图斑划分为建筑阴影；
S53：使用Canny边缘检测算子对各个图斑进行处理，并计算各个图斑在Canny特征层的均值mean canny；
S54：将mean canny大于设置的边缘检测阈值并且面积小于1000个像元的图斑划分为建筑阴影；
S55：计算纯净水体中各个图斑的全波段标准差mean_std；
mean_std＝(sd1+sd2+...sdN)/N
其中，sd1,sd2,…,sdN为图斑各个波段的标准差，N为波段数目；
S56：将全波段标准差大于30并且面积小于1500个像元的图斑划分为建筑阴影；
S57：将所述建筑阴影从所述纯净水体中剔除。
进一步的，所述S6包括：
S61：计算遥感影像中各个像元的黄度指数Yellow，并设置高泥沙含量水体的值域范围[p6,p7]；
Yellow＝p(Green)+p(Red)-2p(Blue)
其中，p(Red)为遥感影像中各个像元在红波段的像元值；
S62：将遥感影像中各个像元的黄度指数与所述p6、p7进行比较，将黄度指数∈[p6,p7]的像元划分为疑似图斑；
S63：计算各个疑似图斑内各个像元的归一化水体指数的均值mean NDWI、各个疑似图斑的边界指数特征以及各个疑似图斑的亮度值 brightness；
brightness＝(meanB1+meanB2+…+meanBN)/N
其中，meanB1,meanB2,…,meanBN为高泥沙含量水体图斑各个波段的亮度均值，N为波段数目；
S64：将mean NDWI大于设定的阈值P8且边界指数特征小于设定的阈值P9的疑似图斑作为高泥沙含量水体图斑，将brightness小于设定的阈值P10的疑似图斑作为高泥沙含量水体图斑；
S65：将与所述纯净水体的边界指数特征大于0的高泥沙含量水体图斑作为高泥沙含量水体。
进一步的，所述S8包括：
S81：将纯净水体、高泥沙含量水体和异常水体合并，得到水体分布区域；
S82：对所述水体分布区域进行形态学腐蚀和膨胀操作；
S83：将水体分布区域中面积小于设定的面积阈值的图斑删除；
S84：对水体分布区域进行统计分析，得到水体空间分布及面积。
一种遥感影像水体自动化提取装置，所述装置包括：
影像获取模块，用于获取遥感影像并对所述遥感影像进行预处理；其中，所述遥感影像包括红波段、绿波段、蓝波段和近红外波段的光谱信息；
纯净水体提取模块，用于利用遥感影像的归一化水体指数以及近红外波段的光谱信息提取纯净水体；
云层干扰去除模块，用于利用遥感影像的蓝波段的光谱信息提取云层信息，并去除云层信息对纯净水体的干扰；
山体阴影剔除模块，用于将遥感影像变换到HIS空间，利用遥感影像的HIS信息、归一化水体指数、改进的阴影水体指数以及增强阴影水体指数从纯净水体中剔除山体阴影；
建筑阴影剔除模块，用于根据纯净水体中各个图斑的边界强度从纯净水体中剔除建筑阴影；
高泥沙含量水体提取模块，用于利用遥感影像的黄度指数提取高泥沙含量水体；
异常水体提取模块，用于选取异常水体特征样本，根据随机森林分类法分离提取异常水体；
修正模块，用于将纯净水体、高泥沙含量水体和异常水体合并，并进行修正，得到水体空间分布和面积。
进一步的，所述纯净水体提取模块包括：
NDWI计算单元，用于计算遥感影像的各个像元的归一化水体指数；
NDWI＝(p(Green)-p(NIR))/(p(Green)+p(NIR))；
其中，NDWI为归一化水体指数，p(Green)和p(NIR)分别为遥感影像中各个像元在绿波段和近红外波段的像元值；
纯净水体划分单元，用于将遥感影像的各个像元的归一化水体指数与设置的纯净水体阈值进行比较，将归一化水体指数大于纯净水体阈值的像元划分为纯净水体；
补充水体阈值确定单元，用于将遥感影像中各个像元在近红外波段的像元值按照从低到高排序，将从低到高第90％位置的像元值作为补充水体阈值；
第一补充单元，用于将遥感影像中各个像元在近红外波段的像元值与所述补充水体阈值进行比较，将近红外波段的像元值小于补充水体阈值的像元划分为纯净水体；
第二补充单元，用于将所述纯净水体阈值缩小一定数值，将遥感影像的各个像元的归一化水体指数与缩小的纯净水体阈值进行比较，将归一化水体指数大于缩小的纯净水体阈值的像元划分为纯净水体。
进一步的，所述云层干扰去除模块包括：
高亮地物设置单元，用于设置高亮地物的值域范围[p1,p2]；
均值计算单元，用于获取遥感影像中高亮地物的各个像元在蓝波段的像元值的平均值p；
云层分割阈值确定单元，用于将所述平均值P与p1、p2进行比较，若所述平均值P小于p1，则将P1作为云层分割阈值，若所述平均值P大于p2，则将P2作为云层分割阈值，若所述平均值P∈[p1,p2]，则将所述平均值P作为云层分割阈值；
云层划分单元，用于将遥感影像中各个像元在蓝波段的像元值与所述云层分割阈值进行比较，将蓝波段的像元值大于云层分割阈值的像元划分为云层；
第一干扰去除单元，用于将包含在所述纯净水体内部且大小小于200 个像元的图斑划分为纯净水体；
第二干扰去除单元，用于将包含在所述纯净水体内部的云层划分为纯净水体。
进一步的，所述山体阴影剔除模块包括：
第一山体阴影划分单元，用于将所述遥感影像变换到HIS空间，将H＜0.085，并且I＜0.115，并且S＜0.441，并且NDWI＜0.35的像元划分为山体阴影；
MSWI计算单元，用于计算遥感影像的各个像元的改进的阴影水体指数；
MSWI＝(p(Blue)-p(NIR))/(p(NIR))；
其中，MSWI为改进的阴影水体指数，p(Blue)为遥感影像中各个像元在蓝波段的像元值；
第二山体阴影划分单元，用于将遥感影像的各个像元的改进的阴影水体指数与设置的第一阴影阈值进行比较，将改进的阴影水体指数小于第一阴影阈值的像元划分为山体阴影；
ESWI计算单元，用于计算遥感影像的各个像元的增强阴影水体指数；
ESWI＝(p(Blue)-p(Green))/(p(NIR)+p(NIR))；
其中，ESWI为增强阴影水体指数；
第三山体阴影划分单元，用于将遥感影像的各个像元的增强阴影水体指数与设置的第二阴影阈值进行比较，将增强阴影水体指数小于第二阴影阈值的像元划分为山体阴影；
山体阴影剔除单元，用于将所述山体阴影从所述纯净水体中剔除。
进一步的，所述建筑阴影剔除模块包括：
边界指数特征计算单元，用于计算纯净水体中各个图斑的边界指数特征borderindex；
其中，bv为图斑的边界周长，Iv为图斑的长度，Wv为图斑的宽度；
第一建筑阴影划分单元，用于将边界指数特征大于4并且面积小于 1500个像元的图斑划分为建筑阴影；
anny边缘检测单元，用于使用Canny边缘检测算子对各个图斑进行处理，并计算各个图斑在Canny特征层的均值mean canny；
第二建筑阴影划分单元，用于将mean canny大于设置的边缘检测阈值并且面积小于1000个像元的图斑划分为建筑阴影；
全波段标准差计算单元，用于计算纯净水体中各个图斑的全波段标准差mean_std；
mean_std＝(sd1+sd2+...sdN)/N
其中，sd1,sd2,…,sdN为图斑各个波段的标准差，N为波段数目；
第三建筑阴影划分单元，用于将全波段标准差大于30并且面积小于 1500个像元的图斑划分为建筑阴影；
建筑阴影剔除单元，用于将所述建筑阴影从所述纯净水体中剔除。
进一步的，所述高泥沙含量水体提取模块包括：
黄度指数计算单元，用于计算遥感影像中各个像元的黄度指数Yellow，并设置高泥沙含量水体的值域范围[p6,p7]；
Yellow＝p(Green)+p(Red)-2p(Blue)
其中，p(Red)为遥感影像中各个像元在红波段的像元值；
疑似图斑划分单元，用于将遥感影像中各个像元的黄度指数与所述p6、 p7进行比较，将黄度指数∈[p6,p7]的像元划分为疑似图斑；
高泥沙含量水体特征计算单元，用于计算各个疑似图斑内各个像元的归一化水体指数的均值mean NDWI、各个疑似图斑的边界指数特征以及各个疑似图斑的亮度值brightness；
brightness＝(meanB1+meanB2+…+meanBN)/N
其中，meanB1,meanB2,…,meanBN为高泥沙含量水体图斑各个波段的亮度均值，N为波段数目；
高泥沙含量水体图斑确定单元，用于将mean NDWI大于设定的阈值 P8且边界指数特征小于设定的阈值P9的疑似图斑作为高泥沙含量水体图斑，将brightness小于设定的阈值P10的疑似图斑作为高泥沙含量水体图斑；
高泥沙含量水体确定单元，用于将与所述纯净水体的边界指数特征大于0的高泥沙含量水体图斑作为高泥沙含量水体。
进一步的，所述修正模块包括：
合并单元，用于将纯净水体、高泥沙含量水体和异常水体合并，得到水体分布区域；
形态学修正单元，用于对所述水体分布区域进行形态学腐蚀和膨胀操作；
面积修正单元，用于将水体分布区域中面积小于设定的面积阈值的图斑删除；
统计单元，用于对水体分布区域进行统计分析，得到水体空间分布及面积。
本发明具有以下有益效果：
本发明针对常规人工识别和勾绘方法的耗费人力多、精度得不到保证、效率低的问题，可提高水体识别的精度和效率；对常规单一算法的水体区域漏分、非水体区域错分、特殊水体未分以及水体边界拟合程度不够等缺陷，融合了多种水体信息层及判别算法，同时去除多种干扰信息，减少了因建筑、山体等阴影以及云层导致的水体错分问题。
本发明可自动化提取常规纯净水体，减少常规人工判别勾绘水体边界的边界误差不确定、水体遗漏以及耗时、费力等缺点，同时兼顾高含沙量水体、富营养化等特征水体的提取，提高水体空间分布的准确性及全面性。在提高水体信息提取速度的同时，保证了数据精度，在湖库敞开水面面积、河流流域范围、洪水灾后监测以及海水与陆地水边线提取等方面具有广泛的应用前景。
附图说明
图1为本发明的遥感影像水体自动化提取方法的流程图；
图2为原始遥感影像与去除非成像区后的影像对比图；
图3为提取的NDWI特征示意图。
图4为云层信息分布示意图；
图5为去除云层干扰后的纯净水体示意图；
图6为提取的高泥沙含量水体示意图；
图 7为水体空间分布示意图；
图8为本发明的感影像水体自动化提取装置的示意图。
具体实施方式
为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。
实施例1：
本发明实施例提供一种遥感影像水体自动化提取方法，如图1所示，该方法包括：
S1：获取遥感影像并对遥感影像进行预处理；其中，遥感影像包括红波段、绿波段、蓝波段和近红外波段的光谱信息。
本步骤用于获取包含待提取水体的卫星遥感影像并获取卫星遥感影像的基本参量信息，包括空间尺度(分辨率)、波段数量、光谱范围、投影坐标、非成像区值等，遥感影像至少包含“Red(红)、Green(率)、Blue (蓝)、Near Infrared(近红外)”4个光谱波段。该遥感影像既可以为原始影像，也可以为采用经过辐射定标、大气校正、几何校正以及其他影像增强处理后的高级别影像。
本步骤还可以同时获取水体待检测区域矢量，若影像与水体待检测区域矢量投影坐标一致，则不做坐标转换处理；若不一致，则进行坐标转换。
利用待检测区域矢量边界对遥感影像进行裁剪，获取待检测区域的遥感影像；若整景影像为待检测区域，则依据影像的基本参量信息去除非成像区域，或将非成像区域归为非水体区域。
S2：利用遥感影像的归一化水体指数以及近红外波段的光谱信息提取纯净水体。
具体实现方法包括：
S21：计算遥感影像的各个像元的归一化水体指数；
NDWI＝(p(Green)-p(NIR))/(p(Green)+p(NIR))；
其中，NDWI为归一化水体指数，p(Green)和p(NIR)分别为遥感影像中各个像元在绿波段和近红外波段的像元值。
S22：将遥感影像的各个像元的归一化水体指数NDWI与设置的纯净水体阈值ndwi_corewater进行比较，将归一化水体指数NDWI大于纯净水体阈值ndwi_corewater的像元划分为纯净水体，其他像元划分为未分类对象。
S23：将遥感影像中各个像元在近红外波段的像元值p(NIR)按照从低到高排序，将从低到高第90％位置的像元值作为补充水体阈值NIR_ve。
S24：将遥感影像中各个像元在近红外波段的像元值p(NIR)与补充水体阈值NIR_ve进行比较，将近红外波段的像元值p(NIR)小于补充水体阈值NIR_ve的像元划分为纯净水体，其他像元划分为未分类对象。
S23和S24用于对纯净水体进行遗漏补充。
S25：将纯净水体阈值缩小一定数值(例如ndwi_corewater-0.1)，将遥感影像的各个像元的归一化水体指数NDWI与缩小的纯净水体阈值ndwi_corewater-0.1进行比较，将归一化水体指数NDWI大于缩小的纯净水体阈值ndwi_corewater-0.1的像元划分为纯净水体，其他像元划分为未分类对象。
本步骤利用NDWI阈值内缩，去除纯净水体中的混分水体。
本发明的S2中，使用NDWI提取纯净水体(即可信度较高的水体)，并利用p(NIR)及NDWI阈值内缩法补充纯净水体信息，提取的纯净水体较为准确。
S3：利用遥感影像的蓝波段的光谱信息提取云层信息，并去除云层信息对纯净水体的干扰。
云层在p(Blue)波段具有高反射特性，其值一般较大，因此可以根据蓝波段的光谱信息提取云层信息，具体实现方法包括：
S31：设置高亮地物的值域范围[p1,p2]。
S32：获取遥感影像中高亮地物的各个像元在蓝波段的像元值p(Blue) 的平均值p。
S33：将平均值P与p1、p2进行比较，若平均值P小于p1，则将P1 作为云层分割阈值，若平均值P大于p2，则将P2作为云层分割阈值，若平均值P∈[p1,p2]，则将平均值P作为云层分割阈值。
S34：将遥感影像中各个像元在蓝波段的像元值p(Blue)与云层分割阈值进行比较，将蓝波段的像元值大于云层分割阈值的像元划分为云层，其他像元划分为未分类对象。
S35：将包含在纯净水体内部且大小小于200个像元的图斑划分为纯净水体。
S36：将包含在纯净水体内部的云层划分为纯净水体。
判断图斑和云层是否位于纯净水体内部时，可以根据图斑/云层与纯净水体边界特征(Rel.border.to)进行界定，该边界特征反映了图斑/云层与纯净水体的边界重叠度。
其中，water为纯净水体的类别名称，L(polygon,waterN)为polygon对象(即图斑/云层)与纯净水体中的waterN对象的公共边长； Npolygon(distance,water)为与polygon对象相邻的所有water对象，bpolygon为 polygon对象的周长。
当Rel.border.to water大于0.7时，认为图斑/云层位于纯净水体内部。
本发明利用p(Blue)及动态阈云层分割阈值取云层空间分布信息，并利用图斑面积与边界特征，去除云层干扰及对造成的水体混淆。
S4：将遥感影像变换到HIS空间，利用遥感影像的HIS信息、归一化水体指数、改进的阴影水体指数以及增强阴影水体指数从纯净水体中剔除山体阴影。
山体阴影与水体在光谱特征上具有相似性，均为光谱吸收性、低反射性，本发明采用多特征耦合去除山体阴影的影响，具体实现方式为：
S41：将遥感影像变换到HIS空间，将H＜0.085，并且I＜0.115，并且 S＜0.441，并且NDWI＜0.35的像元划分为山体阴影。
S42：计算遥感影像的各个像元的改进的阴影水体指数MSWI。
MSWI＝(p(Blue)-p(NIR))/(p(NIR))。
其中，MSWI为改进的阴影水体指数，p(Blue)为遥感影像中各个像元在蓝波段的像元值。
S43：将遥感影像的各个像元的改进的阴影水体指数MSWI与设置的第一阴影阈值P3(P3一般取值2.5)进行比较，将改进的阴影水体指数 MSWI小于第一阴影阈值P3的像元划分为山体阴影，大于等于P3的像元划分为水体。
S44：计算遥感影像的各个像元的增强阴影水体指数ESWI。
ESWI＝(p(Blue)-p(Green))/(p(NIR)+p(NIR))。
其中，ESWI为增强阴影水体指数。
S45：将遥感影像的各个像元的增强阴影水体指数ESWI与设置的第二阴影阈值P4进行比较，将增强阴影水体指数ESWI小于第二阴影阈值 P4的像元划分为山体阴影，大于等于P4的像元划分为水体。
S46：将山体阴影从纯净水体中剔除。
本发明中，构建HIS空间变换、结合MSWI、ESWI谱间关系特征等协同去除山体阴影。
S5：根据纯净水体中各个图斑的边界强度从纯净水体中剔除建筑阴影。
水体与建筑物的边缘强度存在差异，水体边缘强度较小，建筑及建筑阴影的边缘强度较大，本发明采用多特征耦合去除建筑阴影的影响，具体实现方式为：
S51：计算纯净水体中各个图斑的边界指数特征border index。
其中，bv为图斑的边界周长，Iv为图斑的长度，Wv为图斑的宽度。
S52：将边界指数特征border index大于4并且面积小于1500个像元的图斑划分为建筑阴影。
S53：使用Canny边缘检测算子对各个图斑进行处理，并计算各个图斑在Canny特征层的均值mean canny。
S54：将mean canny大于设置的边缘检测阈值P5(一般取值为0.2) 并且面积小于1000个像元的图斑划分为建筑阴影。
S55：计算纯净水体中各个图斑的全波段标准差mean_std。
mean_std＝(sd1+sd2+...sdN)/N
其中，sd1,sd2,…,sdN为图斑各个波段的标准差，N为波段数目。
S56：将全波段标准差大于30并且面积小于1500个像元的图斑划分为建筑阴影。
S57：将建筑阴影从纯净水体中剔除。
本发明中，利用水体图斑面积、边界指数、Canny边缘强度以及全波段标准差等去除建筑阴影，去除阴影误提及水体混淆。
S6：利用遥感影像的黄度指数提取高泥沙含量水体。
泥沙含量高的水体，黄度指数yellow高，但低于黄土地的黄度指数，因此采用yellow特征提取高泥沙含量水体，具体实现方式为：
S61：计算遥感影像中各个像元的黄度指数Yellow，并设置高泥沙含量水体的值域范围[p6,p7]，一般P6为80，P7为100。
Yellow＝p(Green)+p(Red)-2p(Blue)
其中，p(Red)为遥感影像中各个像元在红波段的像元值。
S62：将遥感影像中各个像元的黄度指数与所述p6、p7进行比较，将黄度指数∈[p6,p7]的像元划分为疑似图斑。
S63：计算各个疑似图斑内各个像元的归一化水体指数的均值mean NDWI、各个疑似图斑的边界指数特征以及各个疑似图斑的亮度值 brightness。
brightness＝(meanB1+meanB2+…+meanBN)/N
其中，meanB1,meanB2,…,meanBN为高泥沙含量水体图斑各个波段的亮度均值，N为波段数目。
S64：将mean NDWI大于设定的阈值P8且边界指数特征小于设定的阈值P9的疑似图斑作为高泥沙含量水体图斑，将brightness小于设定的阈值P10的疑似图斑作为高泥沙含量水体图斑。
S65：将与所述纯净水体的边界指数特征大于0的高泥沙含量水体图斑作为高泥沙含量水体。
本发明将黄度指数Yellow∈[p6，p7]，且mean NDWI大于P8，且边界指数小于P9范围内的对象作为高泥沙含量水体。其次将亮度值 brightness小于P10的高泥沙含量水体，采用边界指数特征(计算方法参见S51)，取高泥沙含量水体与纯净水体边界指数大于0的归为高泥沙含量水体类，边界指数小于0的归为未分类。
本发明利用黄度指、亮度特征、NDWI均值特征及边界指数特征等协同提取高泥沙含量水体，以增加该类型的特征水体。
S7：选取异常水体特征样本，根据随机森林分类法分离提取异常水体。
本步骤用于提取富营养化水体、黑臭水体、水色异常水体等异常水体，提取时，选取该部分水体特征样本，采用随机森林分类法，设置特征样本，分类提取该部分水体图斑。
S8：将纯净水体、高泥沙含量水体和异常水体合并，并进行修正，得到水体空间分布和面积。
具体方法包括：
S81：将纯净水体、高泥沙含量水体和异常水体合并，得到水体分布区域。
S82：对水体分布区域进行形态学腐蚀和膨胀操作。
获取的水体分布区域为单排像元、单列像元等形状，采用数学形态学法的腐蚀和膨胀算法，去除该部分误差，具体的定义为：
腐蚀定义为：
膨胀定义为：
S83：将水体分布区域中面积小于设定的面积阈值的图斑删除。
本步骤用于去除水域最小图斑，部分水体所包含像元数较少，设置最小包含像元数据min_pixel(即面积阈值)，将面积小于min_pixel的图斑去除，删除过碎图斑。
S84：对水体分布区域进行统计分析，得到水体空间分布及面积。
本发明首先获取卫星遥感影像及成像参数并进行预处理，其次基于影像的光谱、梯度及统计特征等，构建水体NDWI、MSWI、ESWI、Yellow、 Canny等特征信息；利用NDWI和p(NIR)特征识别纯净水体，然后耦合特征信息，综合去除云层、建筑和山体阴影的错分；构造黄度特征以及特征水体样本补充高含沙量水体及其他异常水体；最后通过边界处理及最小图斑调整，获得普适水体的空间分布。
本发明针对常规人工识别和勾绘方法的耗费人力多、精度得不到保证、效率低的问题，可提高水体识别的精度和效率；对常规单一算法的水体区域漏分、非水体区域错分、特殊水体未分以及水体边界拟合程度不够等缺陷，融合了多种水体信息层及判别算法，同时去除多种干扰信息，减少了因建筑、山体等阴影以及云层导致的水体错分问题。
本发明可自动化提取常规纯净水体，减少常规人工判别勾绘水体边界的边界误差不确定、水体遗漏以及耗时、费力等缺点，同时兼顾高含沙量水体、富营养化等特征水体的提取，提高水体空间分布的准确性及全面性。在提高水体信息提取速度的同时，保证了数据精度，在湖库敞开水面面积、河流流域范围、洪水灾后监测以及海水与陆地水边线提取等方面具有广泛的应用前景。
下面通过一个具体的示例对本发明进行详细说明：
1、获取卫星遥感影像及成像参数，至少包括空间分辨率(spatial resolution)、波段数量(bands)、光谱范围、投影坐标(spatialreference)、非成像区值(DataIgnoreValue)等。
其中，“spatial resolution＝2m”、“bands＝4”、“band1＝p(Blue),band2 ＝p(Green),band3＝p(Red),band4＝p(NIR)”,“spatialreference＝WGS84 UTM51N”，“DataIgnoreValue＝0”。
2、获取水体提取的目标区域，可选择影像的成像区或者由ROI感兴趣区裁剪得到，其中ROI感兴趣区域与非成像区的空间关系自动进行判断，以两者的交集区域作为水体信息的提取区域。
如图2所示，图2中(a)为原始遥感影像，(b)为原始遥感影像去除非成像区后的影像。
3、计算水体信息提取的特征波段数据，至少包括NDWI、MSWI、 ESWI、Yellow、Canny等特征。
图3为提取的NDWI特征。
4、依据云层的高亮特征提取可能造成影响的云层分布特征，如图4 所示。
5、利用NDWI和p(NIR)特征识别纯净水体，并利用步骤4中云层空间分布信息，去除部分云层对水体范围的影响，如图5所示。
6、利用HIS变换、MSWI、ESWI、边界指数、canny、面积及标准差等特征去除步骤5中误提的建筑、山体阴影。
7、利用Yellow指数、边界指数、亮度指数、标准差等特征，提取高泥沙含量水体，如图6所示。
8、选择其他类型水体样本，采用随机森林分类法提取异常水体。
9、将纯净水体、高泥沙含量水体和异常水体合并，并进行水体图斑边界调整及最小图斑去除，获得最终的水体空间分布数据。
图7为最终的水体空间分布，图7中，(a)为水体空间分布整体图， (b)、(c)为水体空间分布局部图
实施例2：
本发明实施例提供一种遥感影像水体自动化提取装置，如图8所示，该装置包括：
影像获取模块1，用于获取遥感影像并对遥感影像进行预处理；其中，遥感影像包括红波段、绿波段、蓝波段和近红外波段的光谱信息。
纯净水体提取模块2，用于利用遥感影像的归一化水体指数以及近红外波段的光谱信息提取纯净水体。
云层干扰去除模块3，用于利用遥感影像的蓝波段的光谱信息提取云层信息，并去除云层信息对纯净水体的干扰。
山体阴影剔除模块4，用于将遥感影像变换到HIS空间，利用遥感影像的HIS信息、归一化水体指数、改进的阴影水体指数以及增强阴影水体指数从纯净水体中剔除山体阴影。
建筑阴影剔除模块5，用于根据纯净水体中各个图斑的边界强度从纯净水体中剔除建筑阴影。
高泥沙含量水体提取模块6，用于利用遥感影像的黄度指数提取高泥沙含量水体。
异常水体提取模块7，用于选取异常水体特征样本，根据随机森林分类法分离提取异常水体。
修正模块8，用于将纯净水体、高泥沙含量水体和异常水体合并，并进行修正，得到水体空间分布和面积。
其中，前述的纯净水体提取模块包括：
NDWI计算单元，用于计算遥感影像的各个像元的归一化水体指数。
NDWI＝(p(Green)-p(NIR))/(p(Green)+p(NIR))。
其中，NDWI为归一化水体指数，p(Green)和p(NIR)分别为遥感影像中各个像元在绿波段和近红外波段的像元值。
纯净水体划分单元，用于将遥感影像的各个像元的归一化水体指数与设置的纯净水体阈值进行比较，将归一化水体指数大于纯净水体阈值的像元划分为纯净水体。
补充水体阈值确定单元，用于将遥感影像中各个像元在近红外波段的像元值按照从低到高排序，将从低到高第90％位置的像元值作为补充水体阈值。
第一补充单元，用于将遥感影像中各个像元在近红外波段的像元值与补充水体阈值进行比较，将近红外波段的像元值小于补充水体阈值的像元划分为纯净水体。
第二补充单元，用于将纯净水体阈值缩小一定数值，将遥感影像的各个像元的归一化水体指数与缩小的纯净水体阈值进行比较，将归一化水体指数大于缩小的纯净水体阈值的像元划分为纯净水体。
前述的云层干扰去除模块包括：
高亮地物设置单元，用于设置高亮地物的值域范围[p1,p2]。
均值计算单元，用于获取遥感影像中高亮地物的各个像元在蓝波段的像元值的平均值p。
云层分割阈值确定单元，用于将平均值P与p1、p2进行比较，若平均值P小于p1，则将P1作为云层分割阈值，若平均值P大于p2，则将 P2作为云层分割阈值，若平均值P∈[p1,p2]，则将平均值P作为云层分割阈值。
云层划分单元，用于将遥感影像中各个像元在蓝波段的像元值与云层分割阈值进行比较，将蓝波段的像元值大于云层分割阈值的像元划分为云层。
第一干扰去除单元，用于将包含在纯净水体内部且大小小于200个像元的图斑划分为纯净水体。
第二干扰去除单元，用于将包含在纯净水体内部的云层划分为纯净水体。
山体阴影剔除模块包括：
第一山体阴影划分单元，用于将遥感影像变换到HIS空间，将H＜0.085，并且I＜0.115，并且S＜0.441，并且NDWI＜0.35的像元划分为山体阴影。
MSWI计算单元，用于计算遥感影像的各个像元的改进的阴影水体指数。
MSWI＝(p(Blue)-p(NIR))/(p(NIR))。
其中，MSWI为改进的阴影水体指数，p(Blue)为遥感影像中各个像元在蓝波段的像元值。
第二山体阴影划分单元，用于将遥感影像的各个像元的改进的阴影水体指数与设置的第一阴影阈值进行比较，将改进的阴影水体指数小于第一阴影阈值的像元划分为山体阴影。
ESWI计算单元，用于计算遥感影像的各个像元的增强阴影水体指数。
ESWI＝(p(Blue)-p(Green))/(p(NIR)+p(NIR))。
其中，ESWI为增强阴影水体指数。
第三山体阴影划分单元，用于将遥感影像的各个像元的增强阴影水体指数与设置的第二阴影阈值进行比较，将增强阴影水体指数小于第二阴影阈值的像元划分为山体阴影。
山体阴影剔除单元，用于将山体阴影从纯净水体中剔除。
建筑阴影剔除模块包括：
边界指数特征计算单元，用于计算纯净水体中各个图斑的边界指数特征borderindex。
其中，bv为图斑的边界周长，Iv为图斑的长度，Wv为图斑的宽度。
第一建筑阴影划分单元，用于将边界指数特征大于4并且面积小于 1500个像元的图斑划分为建筑阴影。
anny边缘检测单元，用于使用Canny边缘检测算子对各个图斑进行处理，并计算各个图斑在Canny特征层的均值mean canny。
第二建筑阴影划分单元，用于将mean canny大于设置的边缘检测阈值并且面积小于1000个像元的图斑划分为建筑阴影。
全波段标准差计算单元，用于计算纯净水体中各个图斑的全波段标准差mean_std。
mean_std＝(sd1+sd2+...sdN)/N
其中，sd1,sd2,…,sdN为图斑各个波段的标准差，N为波段数目。
第三建筑阴影划分单元，用于将全波段标准差大于30并且面积小于 1500个像元的图斑划分为建筑阴影。
建筑阴影剔除单元，用于将建筑阴影从纯净水体中剔除。
高泥沙含量水体提取模块包括：
黄度指数计算单元，用于计算遥感影像中各个像元的黄度指数Yellow，并设置高泥沙含量水体的值域范围[p6,p7]。
Yellow＝p(Green)+p(Red)-2p(Blue)
其中，p(Red)为遥感影像中各个像元在红波段的像元值。
疑似图斑划分单元，用于将遥感影像中各个像元的黄度指数与所述p6、 p7进行比较，将黄度指数∈[p6,p7]的像元划分为疑似图斑。
高泥沙含量水体特征计算单元，用于计算各个疑似图斑内各个像元的归一化水体指数的均值mean NDWI、各个疑似图斑的边界指数特征以及各个疑似图斑的亮度值brightness。
brightness＝(meanB1+meanB2+…+meanBN)/N
其中，meanB1,meanB2,…,meanBN为高泥沙含量水体图斑各个波段的亮度均值，N为波段数目。
高泥沙含量水体图斑确定单元，用于将mean NDWI大于设定的阈值 P8且边界指数特征小于设定的阈值P9的疑似图斑作为高泥沙含量水体图斑，将brightness小于设定的阈值P10的疑似图斑作为高泥沙含量水体图斑。
高泥沙含量水体确定单元，用于将与所述纯净水体的边界指数特征大于0的高泥沙含量水体图斑作为高泥沙含量水体。
修正模块包括：
合并单元，用于将纯净水体、高泥沙含量水体和异常水体合并，得到水体分布区域。
形态学修正单元，用于对水体分布区域进行形态学腐蚀和膨胀操作。
面积修正单元，用于将水体分布区域中面积小于设定的面积阈值的图斑删除。
统计单元，用于对水体分布区域进行统计分析，得到水体空间分布及面积。
本发明实施例所提供的装置，其实现原理及产生的技术效果和前述方法实施例相同，为简要描述，装置实施例部分未提及之处，可参考前述方法实施例1中相应内容。所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，前述描述的装置和单元的具体工作过程，均可以参考上述方法实施例中的对应过程，在此不再赘述。
需要说明的是，本说明书上述所述的装置或者系统根据相关方法实施例的描述还可以包括其他的实施方式，具体的实现方式可以参照方法实施例的描述，在此不作一一赘述。本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于硬件+程序类、存储介质+程序实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
上述对本说明书特定实施例进行了描述。其它实施例在所附权利要求书的范围内。在一些情况下，在权利要求书中记载的动作或步骤可以按照不同于实施例中的顺序来执行并且仍然可以实现期望的结果。另外，在附图中描绘的过程不一定按照要求示出的特定顺序或者连续顺序才能实现期望的结果。在某些实施方式中，多任务处理和并行处理也是可以的或者可能是有利的。
上述实施例阐明的系统、装置、模块或单元，具体可以由计算机芯片或实体实现，或者由具有某种功能的产品来实现。一种典型的实现设备为计算机。具体的，计算机例如可以为个人计算机、膝上型计算机、车载人机交互设备、蜂窝电话、相机电话、智能电话、个人数字助理、媒体播放器、导航设备、电子邮件设备、游戏控制台、平板计算机、可穿戴设备或者这些设备中的任何设备的组合。
为了描述的方便，描述以上装置时以功能分为各种模块分别描述。当然，在实施本说明书一个或多个时可以把各模块的功能在同一个或多个软件和/或硬件中实现，也可以将实现同一功能的模块由多个子模块或子单元的组合实现等。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。
本领域技术人员也知道，除了以纯计算机可读程序代码方式实现控制器以外，完全可以通过将方法步骤进行逻辑编程来使得控制器以逻辑门、开关、专用集成电路、可编程逻辑控制器和嵌入微控制器等的形式来实现相同功能。因此这种控制器可以被认为是一种硬件部件，而对其内部包括的用于实现各种功能的装置也可以视为硬件部件内的结构。或者甚至，可以将用于实现各种功能的装置视为既可以是实现方法的软件模块又可以是硬件部件内的结构。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
在一个典型的配置中，计算设备包括一个或多个处理器(CPU)、输入/ 输出接口、网络接口和内存。
还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、商品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、商品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法或者设备中还存在另外的相同要素。
本领域技术人员应明白，本说明书一个或多个实施例可提供为方法、系统或计算机程序产品。因此，本说明书一个或多个实施例可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本说明书一个或多个实施例可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本说明书一个或多个实施例可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本说明书一个或多个实施例，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本底和远程计算机存储介质中。
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本说明书的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述并不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围。都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。
Claims (8)
1.一种遥感影像水体自动化提取方法，其特征在于，所述方法包括：
S1：获取遥感影像并对所述遥感影像进行预处理；其中，所述遥感影像包括红波段、绿波段、蓝波段和近红外波段的光谱信息；
S2：利用遥感影像的归一化水体指数以及近红外波段的光谱信息提取纯净水体；
S3：利用遥感影像的蓝波段的光谱信息提取云层信息，并去除云层信息对纯净水体的干扰；
S4：将遥感影像变换到HIS空间，利用遥感影像的HIS信息、归一化水体指数、改进的阴影水体指数以及增强阴影水体指数从纯净水体中剔除山体阴影；
S5：根据纯净水体中各个图斑的边界强度从纯净水体中剔除建筑阴影；
S6：利用遥感影像的黄度指数提取高泥沙含量水体；
S7：选取异常水体特征样本，根据随机森林分类法分离提取异常水体；
S8：将纯净水体、高泥沙含量水体和异常水体合并，并进行修正，得到水体空间分布和面积；
其中，所述S5包括：
S51：计算纯净水体中各个图斑的边界指数特征border index；
其中，bv为图斑的边界周长，Iv为图斑的长度，Wv为图斑的宽度；
S52：将边界指数特征大于4并且面积小于1500个像元的图斑划分为建筑阴影；
S53：使用Canny边缘检测算子对各个图斑进行处理，并计算各个图斑在Canny特征层的均值mean canny；
S54：将mean canny大于设置的边缘检测阈值并且面积小于1000个像元的图斑划分为建筑阴影；
S55：计算纯净水体中各个图斑的全波段标准差mean_std；
mean_std＝(sd1+sd2+...sdN)/N
其中，sd1,sd2,…,sdN为图斑各个波段的标准差，N为波段数目；
S56：将全波段标准差大于30并且面积小于1500个像元的图斑划分为建筑阴影；
S57：将所述建筑阴影从所述纯净水体中剔除；
所述S6包括：
S61：计算遥感影像中各个像元的黄度指数Yellow，并设置高泥沙含量水体的值域范围[p6,p7]；
Yellow＝p(Green)+p(Red)-2p(Blue)
其中，p(Red)为遥感影像中各个像元在红波段的像元值；
S62：将遥感影像中各个像元的黄度指数与所述p6、p7进行比较，将黄度指数∈[p6,p7]的像元划分为疑似图斑；
S63：计算各个疑似图斑内各个像元的归一化水体指数的均值mean NDWI、各个疑似图斑的边界指数特征以及各个疑似图斑的亮度值brightness；
brightness＝(meanB1+meanB2+…+meanBN)/N
其中，meanB1,meanB2,…,meanBN为高泥沙含量水体图斑各个波段的亮度均值，N为波段数目；
S64：将mean NDWI大于设定的阈值P8且边界指数特征小于设定的阈值P9的疑似图斑作为高泥沙含量水体图斑，将brightness小于设定的阈值P10的疑似图斑作为高泥沙含量水体图斑；
S65：将与所述纯净水体的边界指数特征大于0的高泥沙含量水体图斑作为高泥沙含量水体。
2.根据权利要求1所述的遥感影像水体自动化提取方法，其特征在于，所述S2包括：
S21：计算遥感影像的各个像元的归一化水体指数；
NDWI＝(p(Green)-p(NIR))/(p(Green)+p(NIR))；
其中，NDWI为归一化水体指数，p(Green)和p(NIR)分别为遥感影像中各个像元在绿波段和近红外波段的像元值；
S22：将遥感影像的各个像元的归一化水体指数与设置的纯净水体阈值进行比较，将归一化水体指数大于纯净水体阈值的像元划分为纯净水体；
S23：将遥感影像中各个像元在近红外波段的像元值按照从低到高排序，将从低到高第90％位置的像元值作为补充水体阈值；
S24：将遥感影像中各个像元在近红外波段的像元值与所述补充水体阈值进行比较，将近红外波段的像元值小于补充水体阈值的像元划分为纯净水体；
S25：将所述纯净水体阈值缩小一定数值，将遥感影像的各个像元的归一化水体指数与缩小的纯净水体阈值进行比较，将归一化水体指数大于缩小的纯净水体阈值的像元划分为纯净水体。
3.根据权利要求2所述的遥感影像水体自动化提取方法，其特征在于，所述S3包括：
S31：设置高亮地物的值域范围[p1,p2]；
S32：获取遥感影像中高亮地物的各个像元在蓝波段的像元值的平均值p；
S33：将所述平均值P与p1、p2进行比较，若所述平均值P小于p1，则将P1作为云层分割阈值，若所述平均值P大于p2，则将P2作为云层分割阈值，若所述平均值P∈[p1,p2]，则将所述平均值P作为云层分割阈值；
S34：将遥感影像中各个像元在蓝波段的像元值与所述云层分割阈值进行比较，将蓝波段的像元值大于云层分割阈值的像元划分为云层；
S35：将包含在所述纯净水体内部且大小小于200个像元的图斑划分为纯净水体；
S36：将包含在所述纯净水体内部的云层划分为纯净水体。
4.根据权利要求3所述的遥感影像水体自动化提取方法，其特征在于，所述S4包括：
S41：将所述遥感影像变换到HIS空间，将H＜0.085，并且I＜0.115，并且S＜0.441，并且NDWI＜0.35的像元划分为山体阴影；
S42：计算遥感影像的各个像元的改进的阴影水体指数；
MSWI＝(p(Blue)-p(NIR))/(p(NIR))；
其中，MSWI为改进的阴影水体指数，p(Blue)为遥感影像中各个像元在蓝波段的像元值；
S43：将遥感影像的各个像元的改进的阴影水体指数与设置的第一阴影阈值进行比较，将改进的阴影水体指数小于第一阴影阈值的像元划分为山体阴影；
S44：计算遥感影像的各个像元的增强阴影水体指数；
ESWI＝(p(Blue)-p(Green))/(p(NIR)+p(NIR))；
其中，ESWI为增强阴影水体指数；
S45：将遥感影像的各个像元的增强阴影水体指数与设置的第二阴影阈值进行比较，将增强阴影水体指数小于第二阴影阈值的像元划分为山体阴影；
S46：将所述山体阴影从所述纯净水体中剔除。
5.根据权利要求4所述的遥感影像水体自动化提取方法，其特征在于，所述S8包括：
S81：将纯净水体、高泥沙含量水体和异常水体合并，得到水体分布区域；
S82：对所述水体分布区域进行形态学腐蚀和膨胀操作；
S83：将水体分布区域中面积小于设定的面积阈值的图斑删除；
S84：对水体分布区域进行统计分析，得到水体空间分布及面积。
6.一种遥感影像水体自动化提取装置，其特征在于，所述装置包括：
影像获取模块，用于获取遥感影像并对所述遥感影像进行预处理；其中，所述遥感影像包括红波段、绿波段、蓝波段和近红外波段的光谱信息；
纯净水体提取模块，用于利用遥感影像的归一化水体指数以及近红外波段的光谱信息提取纯净水体；
云层干扰去除模块，用于利用遥感影像的蓝波段的光谱信息提取云层信息，并去除云层信息对纯净水体的干扰；
山体阴影剔除模块，用于将遥感影像变换到HIS空间，利用遥感影像的HIS信息、归一化水体指数、改进的阴影水体指数以及增强阴影水体指数从纯净水体中剔除山体阴影；
建筑阴影剔除模块，用于根据纯净水体中各个图斑的边界强度从纯净水体中剔除建筑阴影；
高泥沙含量水体提取模块，用于利用遥感影像的黄度指数提取高泥沙含量水体；
异常水体提取模块，用于选取异常水体特征样本，根据随机森林分类法分离提取异常水体；
修正模块，用于将纯净水体、高泥沙含量水体和异常水体合并，并进行修正，得到水体空间分布和面积；
其中，所述建筑阴影剔除模块包括：
边界指数特征计算单元，用于计算纯净水体中各个图斑的边界指数特征borderindex；
其中，bv为图斑的边界周长，Iv为图斑的长度，Wv为图斑的宽度；
第一建筑阴影划分单元，用于将边界指数特征大于4并且面积小于1500个像元的图斑划分为建筑阴影；
anny边缘检测单元，用于使用Canny边缘检测算子对各个图斑进行处理，并计算各个图斑在Canny特征层的均值mean canny；
第二建筑阴影划分单元，用于将mean canny大于设置的边缘检测阈值并且面积小于1000个像元的图斑划分为建筑阴影；
全波段标准差计算单元，用于计算纯净水体中各个图斑的全波段标准差mean_std；
mean_std＝(sd1+sd2+...sdN)/N
其中，sd1,sd2,…,sdN为图斑各个波段的标准差，N为波段数目；
第三建筑阴影划分单元，用于将全波段标准差大于30并且面积小于1500个像元的图斑划分为建筑阴影；
建筑阴影剔除单元，用于将所述建筑阴影从所述纯净水体中剔除；
所述高泥沙含量水体提取模块包括：
黄度指数计算单元，用于计算遥感影像中各个像元的黄度指数Yellow，并设置高泥沙含量水体的值域范围[p6,p7]；
Yellow＝p(Green)+p(Red)-2p(Blue)
其中，p(Red)为遥感影像中各个像元在红波段的像元值；
疑似图斑划分单元，用于将遥感影像中各个像元的黄度指数与所述p6、p7进行比较，将黄度指数∈[p6,p7]的像元划分为疑似图斑；
高泥沙含量水体特征计算单元，用于计算各个疑似图斑内各个像元的归一化水体指数的均值mean NDWI、各个疑似图斑的边界指数特征以及各个疑似图斑的亮度值brightness；
brightness＝(meanB1+meanB2+…+meanBN)/N
其中，meanB1,meanB2,…,meanBN为高泥沙含量水体图斑各个波段的亮度均值，N为波段数目；
高泥沙含量水体图斑确定单元，用于将mean NDWI大于设定的阈值P8且边界指数特征小于设定的阈值P9的疑似图斑作为高泥沙含量水体图斑，将brightness小于设定的阈值P10的疑似图斑作为高泥沙含量水体图斑；
高泥沙含量水体确定单元，用于将与所述纯净水体的边界指数特征大于0的高泥沙含量水体图斑作为高泥沙含量水体。
7.根据权利要求6所述的遥感影像水体自动化提取装置，其特征在于，所述纯净水体提取模块包括：
NDWI计算单元，用于计算遥感影像的各个像元的归一化水体指数；
NDWI＝(p(Green)-p(NIR))/(p(Green)+p(NIR))；
其中，NDWI为归一化水体指数，p(Green)和p(NIR)分别为遥感影像中各个像元在绿波段和近红外波段的像元值；
纯净水体划分单元，用于将遥感影像的各个像元的归一化水体指数与设置的纯净水体阈值进行比较，将归一化水体指数大于纯净水体阈值的像元划分为纯净水体；
补充水体阈值确定单元，用于将遥感影像中各个像元在近红外波段的像元值按照从低到高排序，将从低到高第90％位置的像元值作为补充水体阈值；
第一补充单元，用于将遥感影像中各个像元在近红外波段的像元值与所述补充水体阈值进行比较，将近红外波段的像元值小于补充水体阈值的像元划分为纯净水体；
第二补充单元，用于将所述纯净水体阈值缩小一定数值，将遥感影像的各个像元的归一化水体指数与缩小的纯净水体阈值进行比较，将归一化水体指数大于缩小的纯净水体阈值的像元划分为纯净水体。
8.根据权利要求7所述的遥感影像水体自动化提取装置，其特征在于，所述云层干扰去除模块包括：
高亮地物设置单元，用于设置高亮地物的值域范围[p1,p2]；
均值计算单元，用于获取遥感影像中高亮地物的各个像元在蓝波段的像元值的平均值p；
云层分割阈值确定单元，用于将所述平均值P与p1、p2进行比较，若所述平均值P小于p1，则将P1作为云层分割阈值，若所述平均值P大于p2，则将P2作为云层分割阈值，若所述平均值P∈[p1,p2]，则将所述平均值P作为云层分割阈值；
云层划分单元，用于将遥感影像中各个像元在蓝波段的像元值与所述云层分割阈值进行比较，将蓝波段的像元值大于云层分割阈值的像元划分为云层；
第一干扰去除单元，用于将包含在所述纯净水体内部且大小小于200个像元的图斑划分为纯净水体；
第二干扰去除单元，用于将包含在所述纯净水体内部的云层划分为纯净水体。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110471313.5A CN113177473B (zh) | 2021-04-29 | 2021-04-29 | 遥感影像水体自动化提取方法和装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110471313.5A CN113177473B (zh) | 2021-04-29 | 2021-04-29 | 遥感影像水体自动化提取方法和装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113177473A CN113177473A (zh) | 2021-07-27 |
| CN113177473B true CN113177473B (zh) | 2021-11-16 |
Family
ID=76925879
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110471313.5A Active CN113177473B (zh) | 2021-04-29 | 2021-04-29 | 遥感影像水体自动化提取方法和装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113177473B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113822255B (zh) * | 2021-11-24 | 2022-03-01 | 腾讯科技（深圳）有限公司 | 一种水体识别方法及相关装置 |
Citations (19)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102930496A (zh) * | 2012-10-29 | 2013-02-13 | 南京信息工程大学 | 基于tm影像的水体信息提取方法 |
| US9131642B2 (en) * | 2011-05-13 | 2015-09-15 | Hydrobio, Inc. | Method and system to control irrigation across large geographic areas using remote sensing, weather and field level data |
| CN105046087A (zh) * | 2015-08-04 | 2015-11-11 | 中国资源卫星应用中心 | 一种遥感卫星多光谱影像的水体信息自动提取方法 |
| CN105139015A (zh) * | 2015-07-24 | 2015-12-09 | 河海大学 | 一种遥感图像水体提取方法 |
| CN105590316A (zh) * | 2015-12-11 | 2016-05-18 | 中国测绘科学研究院 | 面向对象的高分辨率遥感影像阴影提取方法 |
| CN107085851A (zh) * | 2017-05-26 | 2017-08-22 | 环境保护部卫星环境应用中心 | 一种提取水域岸线的方法及装置 |
| CN107103295A (zh) * | 2017-04-20 | 2017-08-29 | 苏州中科天启遥感科技有限公司 | 光学遥感影像云检测方法 |
| CN107862667A (zh) * | 2017-11-23 | 2018-03-30 | 武汉大学 | 一种基于高分辨率遥感影像的城市阴影检测与去除方法 |
| CN108169142A (zh) * | 2017-12-20 | 2018-06-15 | 环境保护部卫星环境应用中心 | 基于遥感影像的水色异常快速定位方法和装置 |
| CN108647738A (zh) * | 2018-05-17 | 2018-10-12 | 中国科学院遥感与数字地球研究所 | 基于多指数的全球尺度遥感影像水体智能提取方法 |
| CN109300133A (zh) * | 2018-11-19 | 2019-02-01 | 珠江水利委员会珠江水利科学研究院 | 一种城市河网区水体提取方法 |
| CN109376600A (zh) * | 2018-09-20 | 2019-02-22 | 中国农业大学 | 多光谱遥感影像综合特征云检测方法及装置 |
| CN109544558A (zh) * | 2018-09-27 | 2019-03-29 | 浙江工业大学 | 一种城市复杂环境下阴影与水体分离方法 |
| CN109815894A (zh) * | 2019-01-23 | 2019-05-28 | 中国石油大学(华东) | 一种针对哨兵2a影像的建筑物阴影提取处理方法 |
| CN110210438A (zh) * | 2019-06-10 | 2019-09-06 | 南京林业大学 | 北方土石山区水土流失监测土地利用/覆盖分类方法 |
| CN111398176A (zh) * | 2020-03-13 | 2020-07-10 | 生态环境部卫星环境应用中心 | 基于像元尺度特征的水体水色异常遥感识别方法和装置 |
| CN111415357A (zh) * | 2020-03-19 | 2020-07-14 | 长光卫星技术有限公司 | 一种基于彩色影像的便携式阴影提取方法 |
| CN111931709A (zh) * | 2020-09-17 | 2020-11-13 | 航天宏图信息技术股份有限公司 | 遥感影像的水体提取方法、装置、电子设备及存储介质 |
| CN112613545A (zh) * | 2020-12-17 | 2021-04-06 | 中国农业科学院农业资源与农业区划研究所 | 一种基于多光谱遥感数据提取水体的方法 |
-
2021
- 2021-04-29 CN CN202110471313.5A patent/CN113177473B/zh active Active
Patent Citations (19)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9131642B2 (en) * | 2011-05-13 | 2015-09-15 | Hydrobio, Inc. | Method and system to control irrigation across large geographic areas using remote sensing, weather and field level data |
| CN102930496A (zh) * | 2012-10-29 | 2013-02-13 | 南京信息工程大学 | 基于tm影像的水体信息提取方法 |
| CN105139015A (zh) * | 2015-07-24 | 2015-12-09 | 河海大学 | 一种遥感图像水体提取方法 |
| CN105046087A (zh) * | 2015-08-04 | 2015-11-11 | 中国资源卫星应用中心 | 一种遥感卫星多光谱影像的水体信息自动提取方法 |
| CN105590316A (zh) * | 2015-12-11 | 2016-05-18 | 中国测绘科学研究院 | 面向对象的高分辨率遥感影像阴影提取方法 |
| CN107103295A (zh) * | 2017-04-20 | 2017-08-29 | 苏州中科天启遥感科技有限公司 | 光学遥感影像云检测方法 |
| CN107085851A (zh) * | 2017-05-26 | 2017-08-22 | 环境保护部卫星环境应用中心 | 一种提取水域岸线的方法及装置 |
| CN107862667A (zh) * | 2017-11-23 | 2018-03-30 | 武汉大学 | 一种基于高分辨率遥感影像的城市阴影检测与去除方法 |
| CN108169142A (zh) * | 2017-12-20 | 2018-06-15 | 环境保护部卫星环境应用中心 | 基于遥感影像的水色异常快速定位方法和装置 |
| CN108647738A (zh) * | 2018-05-17 | 2018-10-12 | 中国科学院遥感与数字地球研究所 | 基于多指数的全球尺度遥感影像水体智能提取方法 |
| CN109376600A (zh) * | 2018-09-20 | 2019-02-22 | 中国农业大学 | 多光谱遥感影像综合特征云检测方法及装置 |
| CN109544558A (zh) * | 2018-09-27 | 2019-03-29 | 浙江工业大学 | 一种城市复杂环境下阴影与水体分离方法 |
| CN109300133A (zh) * | 2018-11-19 | 2019-02-01 | 珠江水利委员会珠江水利科学研究院 | 一种城市河网区水体提取方法 |
| CN109815894A (zh) * | 2019-01-23 | 2019-05-28 | 中国石油大学(华东) | 一种针对哨兵2a影像的建筑物阴影提取处理方法 |
| CN110210438A (zh) * | 2019-06-10 | 2019-09-06 | 南京林业大学 | 北方土石山区水土流失监测土地利用/覆盖分类方法 |
| CN111398176A (zh) * | 2020-03-13 | 2020-07-10 | 生态环境部卫星环境应用中心 | 基于像元尺度特征的水体水色异常遥感识别方法和装置 |
| CN111415357A (zh) * | 2020-03-19 | 2020-07-14 | 长光卫星技术有限公司 | 一种基于彩色影像的便携式阴影提取方法 |
| CN111931709A (zh) * | 2020-09-17 | 2020-11-13 | 航天宏图信息技术股份有限公司 | 遥感影像的水体提取方法、装置、电子设备及存储介质 |
| CN112613545A (zh) * | 2020-12-17 | 2021-04-06 | 中国农业科学院农业资源与农业区划研究所 | 一种基于多光谱遥感数据提取水体的方法 |
Non-Patent Citations (5)
| Title |
|---|
| Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery;Hanqiu hu;《International journal of remote sensing》;20070222;第27卷(第14期);1-9 * |
| 国产高分系列遥感影像水体提取方法研究;刘双童;《中国优秀博硕士学位论文全文数据库(硕士) 基础科学辑》;20200415(第04期);全文 * |
| 基于GF-2高分辨率遥感影像的水体提取方法研究;刘双童 等;《全球定位系统》;20181231;第43卷(第6期);1-7 * |
| 基于对象及隶属规则的海岸水边线提取方法;毕京鹏 等;《海岸工程》;20191231;第38卷(第4期);1-14 * |
| 基于高分遥感影像的黄土高原地区水体高精度提取;孙娜 等;《国土资源遥感》;20171231;第29卷(第4期);1-6 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113177473A (zh) | 2021-07-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108253943B (zh) | 一种基于卫星遥感影像的赤潮浒苔一体化监测方法 | |
| US20210199579A1 (en) | Method and system for urban impervious surface extraction based on remote sensing | |
| WO2017071160A1 (zh) | 一种大幅面遥感图像海陆分割的方法及系统 | |
| CN106650812B (zh) | 一种卫星遥感影像的城市水体提取方法 | |
| US20050175253A1 (en) | Method for producing cloud free and cloud-shadow free images | |
| CN111797712B (zh) | 基于多尺度特征融合网络的遥感影像云与云阴影检测方法 | |
| CN103778627B (zh) | 一种基于sar图像的海域溢油检测方法 | |
| CN111027446B (zh) | 一种高分辨率影像的海岸线自动提取方法 | |
| WO2023039959A1 (zh) | 一种基于金字塔机制的遥感影像海洋与非海区域分割方法 | |
| CN108764132B (zh) | 一种湖泊湿地遥感图像误差检测方法 | |
| CN107358161B (zh) | 一种基于遥感影像分类的海岸线提取方法及系统 | |
| Cai et al. | Study on shadow detection method on high resolution remote sensing image based on HIS space transformation and NDVI index | |
| CN113033315A (zh) | 一种稀土开采高分影像识别与定位方法 | |
| Wang et al. | Extracting oil slick features from VIIRS nighttime imagery using a Gaussian filter and morphological constraints | |
| CN111597930A (zh) | 一种基于遥感云平台的海岸线提取方法 | |
| CN111310771B (zh) | 遥感影像的道路图像提取方法、装置、设备及存储介质 | |
| CN113177473B (zh) | 遥感影像水体自动化提取方法和装置 | |
| Kadhim et al. | Shadow detection from very high resoluton satellite image using grabcut segmentation and ratio-band algorithms | |
| Li et al. | Hybrid cloud detection algorithm based on intelligent scene recognition | |
| CN109785318B (zh) | 基于面线基元关联约束的遥感图像变化检测方法 | |
| CN113284066B (zh) | 遥感影像自动云检测方法和装置 | |
| CN111046783A (zh) | 一种改进分水岭算法的斜坡地质灾害边界提取方法 | |
| El Rai et al. | Integrating deep learning with active contour models in remote sensing image segmentation | |
| Zhang et al. | Dynamic Threshold Selection for the Classification of Large Water Bodies within Landsat-8 OLI Water Index Images | |
| CN114862883A (zh) | 一种目标边缘提取方法、图像分割方法及系统 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |