CN116246272A - 针对国产卫星多光谱图像质量标记的云雪区分方法 - Google Patents
针对国产卫星多光谱图像质量标记的云雪区分方法 Download PDFInfo
- Publication number
- CN116246272A CN116246272A CN202211504485.9A CN202211504485A CN116246272A CN 116246272 A CN116246272 A CN 116246272A CN 202211504485 A CN202211504485 A CN 202211504485A CN 116246272 A CN116246272 A CN 116246272A
- Authority
- CN
- China
- Prior art keywords
- snow
- mountain
- cloud
- mask
- pixels
- 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.)
- Pending
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/70—Labelling scene content, e.g. deriving syntactic or semantic representations
-
- 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/20—Image preprocessing
- G06V10/24—Aligning, centring, orientation detection or correction of the image
-
- 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/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/764—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using classification, e.g. of video objects
-
- 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
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Abstract
本申请涉及遥感图像处理技术领域,提供一种针对国产卫星多光谱图像质量标记的云雪区分方法,目的是针对国产多光谱卫星数据生产质量标记产品应用,云与雪在可见光与近红外四波段多光谱图像中存在光谱不可分的问题,本申请提供了一种自动化的云与雪山检测算法流程。首先通过光谱阈值结合导向滤波边缘修正法进行云检测,并利用低分辨率逐月全球雪山覆盖分类数据判断是否存在雪山,然后利用低分辨率卫星的标准雪覆盖产品、DEM高程数据生成模拟雪山掩膜,最后根据模拟雪山掩膜提取云检测结果中误检的雪山。本申请通过利用多源参照数据从云检测结果中提取出误检测的雪山,涉及的关键技术采用成熟的算法实现,完全自动化且具备较高的鲁棒性与普适性。
Description
技术领域
本申请涉及遥感图像处理技术领域,具体涉及一种针对国产卫星多光谱图像质量标记的云雪区分方法。
背景技术
随着传感器硬件的进步,遥感卫星图像标准数据产品的质量近年来有了很大的提升,数据产品的定量化精度得到改善,同时出现了逐像素质量标记与特性的数据产品。目前科研使用较多的数据源包括MODIS、Landsat以及Sentinel传感器数据,都包含质量标记数据产品,用于逐个像素标记包括云、云下阴影、陆地、水体、积雪等类型,方便数据用户使用时根据使用场景过滤干扰像素。例如研究地表变化与分类需要过滤掉云与云下阴影。
生产遥感图像质量标记数据产品的核心算法是云与云下阴影检测,同时包含水体、积雪等类别的检测。目前国外相对成熟的卫星数据质量标记算法主要依赖定量化光谱阈值。在多光谱卫星数据中以美国的陆地资源系列卫星Landsat为代表,目前主要包括Landsat-4/5TM、Landsat-7ETM+以及Landsat-8OLI传感器数据,Landsat系列卫星空间分辨率30米,目前产品数据集基础是Collection 1/2,标准数据产品中包含QA质量标记波段数据,核心算法采用CFMask[1],该算法综合利用多光谱图像可见光、近红外、短波红外以及热红外波段联合阈值区分云、云下阴影、高亮地表、雪等,相比原来的ACCA算法在云与云下阴影标记的精度方面有显著提升。
中国目前是遥感卫星发射大国,中高分辨率的多光谱传感器是主要的成像载荷,并且近年来很多卫星数据标准产品中开始将逐像素质量标记数据纳入到标准流程中,研究针对最新的国产卫星数据质量标记算法开始受到国内学者的重视。目前国产卫星多光谱数据与国外卫星数据相比,主要的差距在于波段数量少,以及定量化精度不高。多光谱图像多为可见光到近红外的四波段数据。由于缺少中波红外波段,质量标记算法便无法利用高云在红外波段的低温特性,导致云与雪、沙漠等高亮地表在国产卫星多光谱数据中的光谱近似而难以有效区分,成为质量标记算法的主要误差来源。
国内研究者针对国产多光谱卫星质量标记算法的精度提升开展了专门研究,目前主要的思路有三个:一是结合图像处理算法改进精度,如采用超像素分割的修正算法提高云区边缘精度的手段,超像素是指具有相似纹理、颜色、亮度等特征的相邻像素构成的图像块,超像素生成算法主要包括分水岭算法、区域增长算法、基于图的图像分割算法和SLIC算法等,武汉大学王密等[2]使用SLIC算法进行国产卫星在轨云检测,提升厚云边缘精度。武汉大学沈焕峰等[3]利用导向滤波对云进行边缘修正,提高了国产GF-2卫星WFV多光谱图像边缘区域精度。二是结合多源辅助数据改进精度,比如基于云的运动特性,利用相近时相的多景数据获取无云地表参照数据进行云与云下阴影的检测。例如利用MODIS的高质量数据产品优化阈值提升检测精度等。三是利用深度学习算法通过大量数据样本的学习训练,使得模型具备与人工判读近似的、在RGB彩色图像中直接区分云、雪等目标。比如基于U-Net网络的Cloud-AttU[4]云检测方法基于特有特征进行识别,在冰雪干扰情况下也能达到很好的效果。这些研究都为国产多光谱卫星质量标记数据产品提供了很好的借鉴。
参考文献:
[1]Zhu,Zhe,and Curtis E.Woodcock.2012."Object-based cloud and cloudshadow detection in Landsat imagery."Remote Sensing of Environment 118:83-94.doi:10.1016/j.rse.2011.10.028.
[2]Wang Mi,Zhang Zhiqi,Dong Zhipeng,et al.Stream-computing Based HighAccuracy On-board Real-time Cloud Detection for High Resolution OpticalSatellite Imagery[J].Acta Geodaetica et Cartographica Sinica,2018,47(6):76-769.DOI:10.11947/jAGCS.2018.20170618.
[3]Zhiwei Li,Huanfeng Shen,Huifang Li,Guisong Xia,Paolo Gamba,Liangpei Zhang.Multi feature combined cloud and cloud shadow detection inGaoFen-1wide field of view imagery[J].Remote Sensing of Environment 191(2017)342-358.
[4]Guo,Yanan,Xiaoqun Cao,Bainian Liu,and Mei Gao.2020."CloudDetection for Satellite Imagery Using Attention-Based U-Net ConvolutionalNeural Network"Symmetry 12(6):1056.doi:10.3390/sym12061056.
发明内容
本申请实施例提供一种针对国产卫星多光谱图像质量标记的云雪区分方法,目的是针对国产卫星多光谱图像的数据预处理应用,特别是多光谱图像质量标记算法,提供一种自动的云与雪山检测算法流程。
本申请实施例提供一种针对国产卫星多光谱图像质量标记的云雪区分方法,包括:
云检测,对于输入的多光谱图像,通过紧缩光谱阈值的方式检测置信度高的厚云,利用大尺度导向滤波修正云区边缘,生成Byte型单波段云掩膜图像;
判断云掩膜图像中是否包含误检测的雪山,利用构建的低分辨率的、逐月份的全球雪山覆盖分类图,判断所述云掩膜图像是否包含雪山,在确定所述云掩膜图像中包含雪山的情况下,执行后续步骤;
参照数据准备,根据输入的多光谱图像的地理范围与成像时间、网络在线下载覆盖的地理范围、成像日期接近的低分辨率卫星标准雪覆盖产品,提取雪掩膜图像并进行裁切和重投影,获得与所述云掩膜图像相同投影、分辨率与像素尺寸的Byte型单波段雪掩膜图像;根据输入的多光谱图像的地理范围裁切与重投影DEM数字高程,获得与输入图像相同投影、分辨率与像素尺寸的单波段DEM图像;
计算模拟雪山掩膜,模拟雪山掩膜是与输入图像相同投影、分辨率与像素尺寸的Byte型单波段图像,对DEM图像进行二值分类划分为山区与非山区,统计雪掩膜图像中所有雪像素对应DEM图像山区像素值的平均高程值作为雪线值,DEM图像中大于雪线值的像素在模拟雪山掩膜中标记为雪山,DEM图像中小于0的所有像素在模拟雪山掩膜中标记为水体;
云掩膜中区分误检雪山,首先将云掩膜中标记为云,而对应模拟雪山掩膜中为雪山的像素修正为雪山,然后将云掩膜中标记为云,而对应模拟雪山掩膜中为山区,并且像素所在山区联通区域在最大半径200像素范围内与雪山连接的像素修正为雪山;
云掩膜中雪山区域修正,对于每一个雪山区域,判断其周围是否全为云像素,如果是则认为该雪山区域在原输入图像中被厚云覆盖,将该雪山联通区域重新标记为云,得到最终包含云与雪山的质量标记掩膜。
在一个实施例中,所述计算模拟雪山掩膜的方法,具体步骤是首先对DEM图像中大于0的所有像素使用OTSU大津法进行二值分类,在模拟雪山掩膜中分别标记山区与非山区,统计雪掩膜图像中雪像素数量占全部非填充像素的百分比,超过2%则进行雪山雪线检测,统计雪掩膜图像中所有雪像素对应DEM图像山区像素值的平均高程值作为雪线值,DEM图像中大于雪线值的像素在模拟雪山掩膜中标记为雪山,DEM图像中小于0的所有像素在模拟雪山掩膜中标记为水体;模拟雪山掩膜包含的类别有雪山、山区、非山区、水体、填充值,分别标记为不同值,目的是对应云掩膜中的云像素包含的雪山误检,假定误检的雪山只出现在雪山掩膜中的雪山与山区区域,不会出现在水体区域。
在一个实施例中,所述超过2%则进行雪山雪线检测,其中数值2%可根据参照数据使用的分辨率卫星标准雪覆盖产品数据源的差异进行调整。
在一个实施例中,所述将云掩膜中标记为云,而对应模拟雪山掩膜中为雪山的像素修正为雪山,具体步骤是:
对云掩膜中的所有云像素;
如果对应模拟雪山掩膜中的像素为雪山,则云掩膜中对应云像素修正为雪山,同时对模拟雪山掩膜中标记为雪山或山区,而在云掩膜中没有被标记为云的像素重新标记为非山区,此时云掩膜中剩余误检测的雪山区域仅存在于模拟雪山掩膜中标记为山区的部分像素。
在一个实施例中,所述然后将云掩膜中标记为云,而对应模拟雪山掩膜中为山区,并且像素所在山区联通区域在最大半径200像素范围内与雪山连接的像素修正为雪山,具体步骤是对模拟雪山掩膜中的每个联通区域按照最大半径200像素进行形态学膨胀,并且膨胀的范围限制在山区范围,将膨胀后的联通区域在云掩膜中的对应像素重新标记为雪山。
在一个实施例中,所述最大半径200像素可根据数据源的差异进行调整。
在一个实施例中,所述云掩膜中雪山区域修正,具体步骤是:
统计云掩膜中所有雪山像素联通区域,并对每个联通区域进行编号,然后遍历每一个雪山的联通区域,对每个雪山的联通区域的所有像素;
如果像素的相邻像素不全为雪山,则该像素为该雪山联通区域的边界点,对所有边界点像素统计其周围20个像素半径范围内是否全部包含云像素;
如果是,则重新标记该雪山联通区域内的所有像素为云区域。
在一个实施例中,所述20个像素的数值可根据数据源的差异进行调整。
本申请实施例提供的针对国产卫星多光谱图像质量标记的云雪区分方法,目的是针对国产多光谱卫星数据生产质量标记产品应用,云与雪在可见光与近红外四波段多光谱图像中存在光谱不可分的问题,本申请提供了一种自动化的云与雪山检测算法流程。首先通过光谱阈值结合导向滤波边缘修正法进行云检测,并利用低分辨率逐月全球雪山覆盖分类数据判断是否存在雪山,然后利用低分辨率卫星的标准雪覆盖产品、DEM高程数据生成模拟雪山掩膜,最后根据模拟雪山掩膜提取云检测结果中误检的雪山。本申请通过利用多源参照数据从云检测结果中提取出误检测的雪山,涉及的关键技术采用成熟的算法实现,完全自动化且具备较高的鲁棒性与普适性。
本申请与现有技术相比有如下特点:本申请提供了一种针对国产卫星多光谱图像质量标记算法中区分云与雪山的解决方案。算法完全自动化,整个过程无需人机交互,用户仅需对最终的检测结果进行简单的检查。涉及的关键步骤采用成熟的算法实现,具有较高的稳定性与适用性。对于国产卫星海量数据自动生产高精度质量标记数据产品,提供了关键的技术支撑。本技术是国产卫星数据标准产品生产研究中必要的数据预处理过程,通过引入多源参照数据来解决云与雪山在可见光与近红外波段的光谱不可分问题。相较于此前的国产卫星数据质量标记算法出现的雪山误检问题提供了一种可行的技术解决途径。
附图说明
为了更清楚地说明本申请或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍,显而易见地,下面描述中的附图是本申请的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1是本申请实施例提供的针对国产卫星多光谱图像质量标记的云雪区分方法的国产卫星遥感图像云与雪山检测流程图;
图2是本申请实施例提供的针对国产卫星多光谱图像质量标记的云雪区分方法的云检测中大尺度导向滤波作用示意图;
图3是本申请实施例提供的针对国产卫星多光谱图像质量标记的云雪区分方法的国产卫星遥感图像云与雪山检测涉及数据实例图;
图4是本申请实施例提供的针对国产卫星多光谱图像质量标记的云雪区分方法利用模拟雪山掩膜识别云掩膜中的雪山误检流程示意图。
具体实施方式
为使本申请的目的、技术方案和优点更加清楚,下面将结合本申请实施例中的附图,对本申请中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本申请一部分实施例,而不是全部的实施例。基于本申请中的实施例,本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例,都属于本申请保护的范围。
本技术的思想基于目前国产卫星多光谱遥感图像数据质量标记算法的两个局限特性:1)红外波段缺失,仅包含可见光与近红外波段,无法利用厚云在红外波段的低温特性进行厚云的检测,导致存在雪山与高亮地表被误检测为云;2)定量化程度不高,通过辐射定标将原始DN(Digital Number,遥感影像像元亮度值)转换到表观反射率的精度有限,难以仅靠光谱阈值有效划分云的精确边界。针对以上两点特性,结合现有相关研究成果,本技术设计一套完全自动化的国产卫星多光谱图像云与雪山检测算法流程,依靠多源参照数据区分雪山误检,以满足对海量国产卫星数据自动、快速生产高精度逐像素质量标记产品数据的需求。
下面结合实施例对本发明提供的针对国产卫星多光谱图像质量标记的云雪区分方法进行详细描述。
图1是本申请实施例提供的针对国产卫星多光谱图像质量标记的云雪区分方法的国产卫星遥感图像云与雪山检测流程图。
本申请实施例提供一种针对国产卫星多光谱图像质量标记的云雪区分方法,包括:
云检测,对于输入的多光谱图像,通过紧缩光谱阈值的方式检测置信度高的厚云,利用大尺度导向滤波修正云区边缘,生成Byte型单波段云掩膜图像;
判断云掩膜图像中是否包含误检测的雪山,利用构建的低分辨率的、逐月份的全球雪山覆盖分类图,判断云掩膜图像是否包含雪山,在确定云掩膜图像中包含雪山的情况下,执行后续步骤;
参照数据准备,根据输入的多光谱图像的地理范围与成像时间、网络在线下载覆盖的地理范围、成像日期接近的低分辨率卫星标准雪覆盖产品,提取雪掩膜图像并进行裁切和重投影,获得与云掩膜图像相同投影、分辨率与像素尺寸的Byte型单波段雪掩膜图像;根据输入的多光谱图像的地理范围裁切与重投影DEM(Digital Elevation Model,数字高程模型)数字高程,获得与输入图像相同投影、分辨率与像素尺寸的单波段DEM图像;
计算模拟雪山掩膜,模拟雪山掩膜是与输入图像相同投影、分辨率与像素尺寸的Byte(字节)型单波段图像,对DEM图像进行二值分类划分为山区与非山区,统计雪掩膜图像中所有雪像素对应DEM图像山区像素值的平均高程值作为雪线值,DEM图像中大于雪线值的像素在模拟雪山掩膜中标记为雪山,DEM图像中小于0的所有像素在模拟雪山掩膜中标记为水体;
云掩膜中区分误检雪山,首先将云掩膜中标记为云,而对应模拟雪山掩膜中为雪山的像素修正为雪山,然后将云掩膜中标记为云,而对应模拟雪山掩膜中为山区,并且像素所在山区联通区域在最大半径200像素范围内与雪山连接的像素修正为雪山;
云掩膜中雪山区域修正,对于每一个雪山区域,判断其周围是否全为云像素,如果是则认为该雪山区域在原输入图像中被厚云覆盖,将该雪山联通区域重新标记为云,得到最终包含云与雪山的质量标记掩膜。
在一个实施例中,计算模拟雪山掩膜的方法,具体步骤是首先对DEM图像中大于0的所有像素使用OTSU大津法进行二值分类,在模拟雪山掩膜中分别标记山区与非山区,统计雪掩膜图像中雪像素数量占全部非填充像素的百分比,超过2%则进行雪山雪线检测,统计雪掩膜图像中所有雪像素对应DEM图像山区像素值的平均高程值作为雪线值,DEM图像中大于雪线值的像素在模拟雪山掩膜中标记为雪山,DEM图像中小于0的所有像素在模拟雪山掩膜中标记为水体;模拟雪山掩膜包含的类别有雪山、山区、非山区、水体、填充值,分别标记为不同值,目的是对应云掩膜中的云像素包含的雪山误检,假定误检的雪山只出现在雪山掩膜中的雪山与山区区域,不会出现在水体区域。
在一个实施例中,超过2%则进行雪山雪线检测,其中数值2%可根据参照数据使用的分辨率卫星标准雪覆盖产品数据源的差异进行调整。
在一个实施例中,将云掩膜中标记为云,而对应模拟雪山掩膜中为雪山的像素修正为雪山,具体步骤是:
对云掩膜中的所有云像素;
如果对应模拟雪山掩膜中的像素为雪山,则云掩膜中对应云像素修正为雪山,同时对模拟雪山掩膜中标记为雪山或山区,而在云掩膜中没有被标记为云的像素重新标记为非山区,此时云掩膜中剩余误检测的雪山区域仅存在于模拟雪山掩膜中标记为山区的部分像素。
在一个实施例中,然后将云掩膜中标记为云,而对应模拟雪山掩膜中为山区,并且像素所在山区联通区域在最大半径200像素范围内与雪山连接的像素修正为雪山,具体步骤是对模拟雪山掩膜中的每个联通区域按照最大半径200像素进行形态学膨胀,并且膨胀的范围限制在山区范围,将膨胀后的联通区域在云掩膜中的对应像素重新标记为雪山。
在一个实施例中,最大半径200像素可根据数据源的差异进行调整。
在一个实施例中,云掩膜中雪山区域修正,具体步骤是:
统计云掩膜中所有雪山像素联通区域,并对每个联通区域进行编号,然后遍历每一个雪山的联通区域,对每个雪山的联通区域的所有像素;
如果像素的相邻像素不全为雪山,则该像素为该雪山联通区域的边界点,对所有边界点像素统计其周围20个像素半径范围内是否全部包含云像素;
如果是,则重新标记该雪山联通区域内的所有像素为云区域。
在一个实施例中,20个像素的数值可根据数据源的差异进行调整。
本申请中包括以下内容:
输入数据:国产多光谱图像,涉及国产军用、民用卫星所搭载的中高分辨率多光谱传感器获取的图像数据,例如GF系列、HJ系列、ZY系列、JB系列等卫星。获取的多光谱图像空间分辨率不低于30米,不包含中波红外或热红外波段,仅包含可见光到近红外的四个波段。输入数据已经通过定标参数由原始DN转换到表观反射率。
输入数据:全球逐月雪覆盖分类图,预先制作并作为参照数据集保存于服务器本地文件供程序调用时查询,制作该数据集可利用现有多种数据源及雪山、雪线数据产品,要求分辨率可大于500米,但不能低于2000米,云盖量为0或尽量低。
输入数据:低分辨率卫星标准雪覆盖产品,空间分辨率低于500米,昼间全天时成像,且同一地表重访周期不超过一周的卫星数据,比如目前使用最多的是MODIS系列卫星MOD10A1、MOD10A2标准雪覆盖产品。
MODIS积雪产品MOD10A1是基于Terra卫星上搭载的多波段传感器影像生成的,卫星的成像时间为地方时间10:30。积雪覆盖范围自动制图是使用卫星在MODIS波段4(0.545~0.565μm)和波段6(1.628~1.652μm)的反射系数所计算的归一化积雪指数(NSDI)方法,结合阈值测试以及MODIS云掩膜数据生成,二值日尺度积雪产品(“Snow Cover DailyTile”)的精度约94%。MOD10A1数据还提供了像元积雪面积比例产品,该产品是根据在阿拉斯加、加拿大和俄罗斯等地区建立的MODIS NDSI与由Landsat影像提取的真实像元的雪覆盖范围之间的线性关系计算生成,平均绝对误差<0.1,在尼泊尔-喜马拉雅地区的均方根误差为0.15,精度较高。MOD10A2则是进一步整合处理得到的8天合成数据产品,利用云的运动特性过滤掉了大部分厚云覆盖,具备更高的精度,但存在部分产品数据无法生产及检索不到的情况。
输入数据:全球DEM高程,空间分辨率不低于90米,比如全球SRTM数据。
处理步骤:云检测,通过紧缩光谱阈值的方式检测置信度高的厚云,紧缩阈值是通过缩小阈值范围,减少检测出的云与阴影数量,达到提高阈值法检测的精度,厚云检测使用HOT指数与VBR指数,检测结果以厚云为主。
HOT指数计算公式为:
HOT=B1-0.5×B3
VBR指数计算公式为:
其中B1、B2、B3分别为蓝波段、绿波段、红波段数据。云粗检测的阈值为HOT>0.2&&VBR>0.7。
云区大尺度导向滤波是利用滤波半径为500像素的大尺度导向滤波修正云区边缘,主要是将云区边缘漏检的薄云添加到云掩膜,适应由云粗检测可能导致的大面积云区漏检的修正。
通过最小化下面的窗口ωk的代价函数:
得到局部线性系数ak,bk的值。
导向滤波算法如下:
step1:meanL=fmean(L(x,y))
meanV=fmean(V(x,y))
corrL=fmean(L(x,y).*L(x,y))
corrLV=fmean(L(x,y).*V(x,y))
step2:varL=corrI-meanI.*meanI
covLV=corrIV-meanI.*meanV
step3:a=covLV./(varL+ε)
b=meanV-a.*meanL
step4:meana=fmean(a)
meanb=fmean(b)
其中fmean为中值滤波器,滤波半径为r,ε为归一化参数。为了适应大面积厚云区与地表区之间的过渡区域,即云区边界或薄云区,使用500像素的大滤波半径。随着滤波半径的增大会导致一般中值滤波器运算效率的急剧增大,采用和半径无关的Boxfilter快速滤波算法。使复杂度为O(MN)的求和、均值、方差等运算降低到近似于O(1)的复杂度。导向滤波的输出图像为修正的云区掩膜。
图2是云检测中大尺度导向滤波作用示意图,表明通过导向滤波图像(图2.c),得到的云掩膜(图2.d)相比之前(图2.b)更符合原始图像(图2.a)中云的真实边界。
处理步骤:判断是否存在雪山,利用全球逐月雪覆盖分类图判断云掩膜中是否存在雪山误检。
根据输入图像地理区域与成像时间,查询对应全球雪山覆盖分类图中该月份该区域包含雪山的雪山像素数占比是否超过2%,超过则进行后续的云雪区分。
如果判断云掩膜中不包含雪山误检测,则不进行后续的云雪区分处理,直接进行后续的质量标记处理流程。
处理中间结果:云掩膜,云掩膜中包含误检测的雪山,云掩膜是图像像素尺寸与输入图像一致的Byte型单波段图像,不同像素值标记不同类别,比如云像素被标记为5值,其清晰地表像素标记为1值,其填充像素标记为0值,掩膜图像像素尺寸与输入图像一致。云掩膜图像保存到硬盘。
图3为国产卫星遥感图像云与雪山检测涉及数据实例图。图3中的B为云掩膜的例子,其中黑圈内的雪山被识别为云而标记为5值,在图像中显示为白色。
处理步骤:自动查询下载,通过编程的方式实现自动查询当前处理数据需要的低分辨率雪覆盖产品数据,并下载到本地硬盘。目前国外MODIS系列卫星的MOD10A1、MOD10A2标准雪覆盖产品可以方便地通过HTTPS请求下载,比如通过美国冰雪数据中心(NSIDC)(http://nsidc.org/data/MOD10A1)下载MOD10A1逐日产品数据格式为hdf,数据投影为正弦地图投影,空间分辨率为500米,基本可保证对于输入的多光谱数据成像日期下能够查询到当天、或者前后1天对应的MOD10A1逐日产品数据。对于输入的多光谱数据地理覆盖区域有可能查询到2至4景国产多光谱图像逐日产品数据,利用命令行调用MODIS ReprojectionTool(MRT)工具的批处理接口实现多景数据的自动拼接,同时选择积雪覆盖率数据(FSC)波段进行hdf到tiff格式的转换。
由于MOD10A1逐日产品数据记录的是当天具体时刻的雪覆盖标记,同时也包含了该时刻云覆盖标记,存在雪山区域被厚云遮挡的情况,而8天合成数据产品MOD10A2很大程度消除了运动的厚云,且8天时间范围内雪山雪线的变化通常不大,所以相比逐日数据更适合本算法使用,但存在部分产品数据无法生产及检索不到的情况。所以数据获取的策略是优先下载MOD10A2,如果失败则下载MOD10A1数据。
处理步骤:裁切重投影,裁切重投影的目的是生成与输入的国产多光谱图像像素尺寸一致的参照数据,包括雪掩膜图像与DEM图像,具体步骤是:首先根据输入的国产多光谱图像信息确定地理范围的外包矩形,为了确保重投影后地理范围依然外包输入图像,适当扩大外包矩形;然后根据外包矩形范围裁切单波段的积雪覆盖图像与全球DEM高程,由于裁切的对象分辨率低于输入图像(MOD10A1积雪覆盖率数据分辨率500米,SRTM高程数据分辨率90米),将裁切的对象按照虚拟文件的方式存储在内存;最后对内存的裁切文件重投影到与输入图像相同的地图投影、分辨率,并裁切到与输入图像一致的地理范围与像素尺寸,保存输出结果的雪覆盖图像与DEM图像到外部文件。
处理中间结果:雪覆盖图像,输入来源于低分辨率雪覆盖产品,是和输入图像像素尺寸一致的Byte型单波段图像数据,雪覆盖图像保存到硬盘。
如果数据来源于MODIS系列卫星MOD10A1/MOD10A2数据,则雪覆盖图像为FSC波段,像素值0-100表示雪覆盖的比例,同时包含水体(239、237)、云(250)与填充值(255)标识。
图3中的C为雪覆盖图像的例子,数据来源于MODIS系列卫星MOD10A1数据FSC波段。
处理中间结果:DEM图像,输入来源于全球DEM高程产品,是和输入图像像素尺寸一致的Int16型单波段图像数据,DEM图像保存到硬盘。
如果数据来源于SRTM数据,则DEM图像像素值表示高度数值,没有高程的像素值为-32768,通常为海洋或大面积湖泊。
图3中的D为DEM高程图像的例子,数据来源于SRTM数据。
处理步骤:雪山模拟,目的是模拟雪山用于识别对应云掩膜中的云像素包含的雪山误检。误检的雪山最可能出现在雪山掩膜中的雪山与山区区域,不会出现在水体区域。雪山模拟的目标是生产包含雪山、山区、非山区、水体类别的雪山掩膜作为参照数据。
具体步骤首先是根据输入的国产多光谱图像生产相同像素尺寸的单波段Byte型图像作为模拟雪山掩膜图像,并根据输入的国产多光谱图像在掩膜中标记图像有效像素值为1,填充区域为0。
对于DEM图像中小于0的所有像素,如果在模拟雪山掩膜对应像素值为1,则重新标记为水体2。
对于DEM图像中像素值大于0,并且在模拟雪山掩膜对应像素值为1的所有像素,使用OTSU大津法进行二值分类,目的是在模拟雪山掩膜中分别标记山区与非山区。大津法是一种确定图像二值化分割阈值的算法,假设图像像素能够根据全局阈值,被分成背景(非山区C1)和目标(山区C2)两部分。然后采用最大类间方差的原理,计算该最佳阈值来区分这两类像素,使得分割后的两部分类间方差最大。
假定C1与C2两类像素均值分别为m1与m2,所有像素的均值为mg,每个像素被分为C1与C2的概率分别为p1与p2,则有:
类间方差为:
σ2=p1(m1-mg)2+p2(m2-mg)2
可得:
σ2=p1p2(m1-m2)2
OTSU阈值kOTSU使得上式最大化,小于kOTSU的像素划分为非山区C1,大于kOTSU的像素划分为山区C2,划分为山区的像素,雪山掩膜中像素值由1重新标记为山区值44。
统计雪掩膜图像中所有雪像素对应DEM图像山区像素值的平均高程值作为雪线阈值ksnow,DEM图像中大于雪线值的像素在模拟雪山掩膜中标记为雪山。具体步骤是首先对DEM图像中大于0的所有像素统计雪掩膜图像中雪像素数量占全部非填充像素的百分比,超过2%则进行雪山雪线检测,统计雪掩膜图像中所有雪像素对应DEM图像山区像素值的平均高程值作为雪线值,DEM图像中大于雪线值的像素在模拟雪山掩膜中标记为雪山。
其中数值2%可根据参照数据使用的分辨率卫星标准雪覆盖产品数据源的差异进行调整,这里的2%数值是针对逐日的MOD10A1雪覆盖产品制定,若采用MODIS系列卫星8天合成的MOD10A2数据,由于云覆盖率明显降低,则可调整为5%。由于固定时刻成像数据容易受厚云遮挡,虽然当前处理数据区域属于雪山存在范围,但有可能出现当天对应MOD10A1数据产品雪山区域被厚云完全覆盖的现象,此时程序将重新自动下载前/后一天的MOD10A1数据产品,直到检索到的雪山覆盖比例超过2%。考虑到雪山雪线的变化缓慢,前/后一周的MOD10A1数据产品都是适用的。
DEM图像中像素值大于雪线阈值ksnow,则在雪山掩膜中对应的像素值标记为雪山值4。通常雪线阈值ksnow大于山区阈值kOTSU,即在雪山掩膜中原本标记为山区的部分像素被标记为雪山。但如果出现雪线阈值ksnow小于山区阈值kOTSU的情况,则认为OTSU法计算的阈值kOTSU划分不准确,将雪山掩膜中的山区像素值44重新标记为1,仅保留雪线阈值ksnow标记的雪山值4。
最终得到的模拟雪山掩膜中,包含的标记值有雪山4、山区44、非山区1、水体2、填充值0。
需要说明的是,模拟雪山掩膜中标记的雪山距离输入的国产多光谱图像中的实际雪山尚有很大差距,主要原因有两个:一是标准雪覆盖产品分辨率低,MODIS数据分辨率仅为500米,且MOD10A1逐日数据容易受云影响,导致估计的雪线阈值精度有限;二是真实的雪线通常高程并不一致,山的阴面与阳面、以及气候、季风等因素导致实际的雪线在30米空间分辨率的图像中复杂多变。模拟雪山掩膜的作用是作为云掩膜中雪山误检的判断提供参照数据,如果云掩膜中标记为云的像素在模拟雪山掩膜中被标记为雪山,则认为该像素真实的情况大概率是雪山,而如果在模拟雪山掩膜中被标记为雪山的像素在云掩膜中没有被标记为云,则该像素的真实情况大概率不是雪山。
处理中间结果:模拟雪山掩膜,模拟雪山掩膜是与输入图像相同投影、分辨率与像素尺寸的Byte型单波段图像,模拟雪山掩膜包含的类别有雪山、山区、非山区、水体、填充值,分别标记为不同值,模拟雪山掩膜保存到硬盘。
模拟雪山掩膜用于识别对应云掩膜中的云像素包含的雪山误检,假定误检的雪山只出现在雪山掩膜中的雪山与山区(高概率出现在雪山,低概率出现在山区),不会出现在水体区域。
图3中的E为模拟雪山掩膜的例子,其中雪山区域被高亮显示。
处理步骤:云掩膜中区分误检雪山,目的是将云掩膜中标记为云像素,但根据模拟雪山掩膜高概率是雪山的像素修改标记为雪山。具体包含两次过滤步骤:
第一遍过滤是针对高概率误检的雪山像素,即在模拟雪山掩膜中标记为雪山的像素。对云掩膜中的所有云像素值5,如果对应模拟雪山掩膜中的像素为雪山值4,则云掩膜中对应云像素修正为雪山值4。同时,如果模拟雪山掩膜中标记为雪山的像素在云掩膜中没有被标记为云,则将模拟雪山掩膜中的该像素重标记为非山区值1,此时云掩膜中标记的雪山区域与模拟雪山掩膜中标记的雪山区域一致。如果模拟雪山掩膜中标记为山区的像素在云掩膜中的标记不为云,则将模拟雪山掩膜中的该像素重标记为非山区值1,此时云掩膜中剩余误检测的雪山区域仅存在于部分模拟雪山掩膜中标记为山区的像素。
参照图4,图4为本申请实施例提供的针对国产卫星多光谱图像质量标记的云雪区分方法利用模拟雪山掩膜识别云掩膜中的雪山误检流程示意图。图4中的云掩膜中区分误检雪山-1步骤显示了处理前后云掩膜与模拟掩膜标记值的变化。
第二遍过滤是对低概率误检的雪山像素,即在模拟雪山掩膜中标记为山区的像素。算法首先统计模拟雪山掩膜中所有标记为雪山值4的像素联通区域,并标号每个联通区域。为了避免联通区域数量过多、形态过离散,在统计并标号雪山联通区域之前,先对雪山像素进行整合修正,将小于20个联通像素的雪山重新标记为山区,将雪山联通区域内部小于20个像素的山区“漏洞”重新标记为雪山。然后遍历每个雪山联通区域,每个联通区域按照最大半径200像素进行形态学膨胀,并且膨胀的范围限制在山区值44所在范围。最后将膨胀后的联通区域在云掩膜中的对应像素重新标记为雪山值4。图4中的云掩膜中区分误检雪山-2步骤显示了处理前后云掩膜与模拟掩膜标记值的变化。
需要说明的是,对于输入的国产多光谱图像,实际存在雪山被厚云完全或部分遮盖的情况
处理步骤:云掩膜中雪山区域修正,目的是识别被厚云完全遮盖的雪山区域并重新修正为厚云。对于每一个雪山区域,判断其周围是否全为云像素,如果是则认为该雪山区域在原输入图像中被厚云覆盖,将该雪山联通区域重新标记为云。具体步骤是统计云掩膜中所有雪山像素联通区域,并对每个联通区域进行编号,然后遍历每一个雪山的联通区域,对每个雪山的联通区域的所有像素,如果像素的4邻域像素不全为雪山,则该像素为该雪山联通区域的边界点,对所有边界点像素统计其周围20个像素半径范围内是否全部包含云像素,如果是则重新标记该联通区域雪山所有像素为云区域。图4中的云掩膜中雪山区域修正步骤显示了处理前后云掩膜与模拟掩膜标记值的变化,最终得到的云与雪山掩膜正确标记了雪山,并修正了被厚云覆盖的雪山。
需要说明的是,对于输入的国产多光谱图像中存在雪山被厚云部分遮盖的情况,本算法无法区分厚云与雪山之间的分割边界,算法目前按照全不遮盖处理,即雪山上方覆盖的部分厚云由于无法区分而标记为雪山。
输出数据:云与雪山掩膜,是和输入图像像素尺寸一致的Byte型单波段图像数据,包含云、雪山、清晰地表、水体、填充值类别,不同类别标记为不同数值,这里填充值标记为0,清晰地表标记为1,水体标记为2,雪山标记为4,云标记为5。
图3中的F为云与雪山掩膜的例子,其中黑圈内的雪山基本被检测出来并标记为4,在图像中显示为蓝色。
需要说明的是,本算法的核心是解决雪山误检测的问题,本算法属于国产多光谱卫星质量标记算法流程中的一部分,本算法输出的云与雪山掩膜也是作为国产多光谱卫星质量标记算法流程中的中间数据存储在硬盘,后续算法还将进行云下阴影检测、云区修正、水体检测与修正等处理,所以云与雪山掩膜中的类别标记不代表质量标记最终结果。
本发明的C++算法实例已经在PC平台上实现,前期通过实验数据验证了算法的有效性与鲁棒性。
最后应说明的是:以上实施例仅用以说明本申请的技术方案,而非对其限制;尽管参照前述实施例对本申请进行了详细的说明,本领域的普通技术人员应当理解:其依然可以对前述各实施例所记载的技术方案进行修改,或者对其中部分技术特征进行等同替换;而这些修改或者替换,并不使相应技术方案的本质脱离本申请各实施例技术方案的精神和范围。
Claims (8)
1.一种针对国产卫星多光谱图像质量标记的云雪区分方法,其特征在于,包括:
云检测,对于输入的多光谱图像,通过紧缩光谱阈值的方式检测置信度高的厚云,利用大尺度导向滤波修正云区边缘,生成Byte型单波段云掩膜图像;
判断云掩膜图像中是否包含误检测的雪山,利用构建的低分辨率的、逐月份的全球雪山覆盖分类图,判断所述云掩膜图像是否包含雪山,在确定所述云掩膜图像中包含雪山的情况下,执行后续步骤;
参照数据准备,根据输入的多光谱图像的地理范围与成像时间、网络在线下载覆盖的地理范围、成像日期接近的低分辨率卫星标准雪覆盖产品,提取雪掩膜图像并进行裁切和重投影,获得与所述云掩膜图像相同投影、分辨率与像素尺寸的Byte型单波段雪掩膜图像;根据输入的多光谱图像的地理范围裁切与重投影DEM数字高程,获得与输入图像相同投影、分辨率与像素尺寸的单波段DEM图像;
计算模拟雪山掩膜,模拟雪山掩膜是与输入图像相同投影、分辨率与像素尺寸的Byte型单波段图像,对DEM图像进行二值分类划分为山区与非山区,统计雪掩膜图像中所有雪像素对应DEM图像山区像素值的平均高程值作为雪线值,DEM图像中大于雪线值的像素在模拟雪山掩膜中标记为雪山,DEM图像中小于0的所有像素在模拟雪山掩膜中标记为水体;
云掩膜中区分误检雪山,首先将云掩膜中标记为云,而对应模拟雪山掩膜中为雪山的像素修正为雪山,然后将云掩膜中标记为云,而对应模拟雪山掩膜中为山区,并且像素所在山区联通区域在最大半径200像素范围内与雪山连接的像素修正为雪山;
云掩膜中雪山区域修正,对于每一个雪山区域,判断其周围是否全为云像素,如果是则认为该雪山区域在原输入图像中被厚云覆盖,将该雪山联通区域重新标记为云,得到最终包含云与雪山的质量标记掩膜。
2.根据权利要求1所述的针对国产卫星多光谱图像质量标记的云雪区分方法,其特征在于,
所述计算模拟雪山掩膜的方法,具体步骤是首先对DEM图像中大于0的所有像素使用OTSU大津法进行二值分类,在模拟雪山掩膜中分别标记山区与非山区,统计雪掩膜图像中雪像素数量占全部非填充像素的百分比,超过2%则进行雪山雪线检测,统计雪掩膜图像中所有雪像素对应DEM图像山区像素值的平均高程值作为雪线值,DEM图像中大于雪线值的像素在模拟雪山掩膜中标记为雪山,DEM图像中小于0的所有像素在模拟雪山掩膜中标记为水体;模拟雪山掩膜包含的类别有雪山、山区、非山区、水体、填充值,分别标记为不同值,目的是对应云掩膜中的云像素包含的雪山误检,假定误检的雪山只出现在雪山掩膜中的雪山与山区区域,不会出现在水体区域。
3.根据权利要求2所述的针对国产卫星多光谱图像质量标记的云雪区分方法,其特征在于,所述超过2%则进行雪山雪线检测,其中数值2%可根据参照数据使用的分辨率卫星标准雪覆盖产品数据源的差异进行调整。
4.根据权利要求1所述的针对国产卫星多光谱图像质量标记的云雪区分方法,其特征在于,
所述将云掩膜中标记为云,而对应模拟雪山掩膜中为雪山的像素修正为雪山,具体步骤是:
对云掩膜中的所有云像素;
如果对应模拟雪山掩膜中的像素为雪山,则云掩膜中对应云像素修正为雪山,同时对模拟雪山掩膜中标记为雪山或山区,而在云掩膜中没有被标记为云的像素重新标记为非山区,此时云掩膜中剩余误检测的雪山区域仅存在于模拟雪山掩膜中标记为山区的部分像素。
5.根据权利要求1所述的针对国产卫星多光谱图像质量标记的云雪区分方法,其特征在于,
所述然后将云掩膜中标记为云,而对应模拟雪山掩膜中为山区,并且像素所在山区联通区域在最大半径200像素范围内与雪山连接的像素修正为雪山,具体步骤是对模拟雪山掩膜中的每个联通区域按照最大半径200像素进行形态学膨胀,并且膨胀的范围限制在山区范围,将膨胀后的联通区域在云掩膜中的对应像素重新标记为雪山。
6.根据权利要求5所述的针对国产卫星多光谱图像质量标记的云雪区分方法,其特征在于,所述最大半径200像素可根据数据源的差异进行调整。
7.根据权利要求1所述的针对国产卫星多光谱图像质量标记的云雪区分方法,其特征在于,所述云掩膜中雪山区域修正,具体步骤是:
统计云掩膜中所有雪山像素联通区域,并对每个联通区域进行编号,然后遍历每一个雪山的联通区域,对每个雪山的联通区域的所有像素;
如果像素的相邻像素不全为雪山,则该像素为该雪山联通区域的边界点,对所有边界点像素统计其周围20个像素半径范围内是否全部包含云像素;
如果是,则重新标记该雪山联通区域内的所有像素为云区域。
8.根据权利要求1所述的针对国产卫星多光谱图像质量标记的云雪区分方法,其特征在于,所述20个像素的数值可根据数据源的差异进行调整。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211504485.9A CN116246272A (zh) | 2022-11-28 | 2022-11-28 | 针对国产卫星多光谱图像质量标记的云雪区分方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211504485.9A CN116246272A (zh) | 2022-11-28 | 2022-11-28 | 针对国产卫星多光谱图像质量标记的云雪区分方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN116246272A true CN116246272A (zh) | 2023-06-09 |
Family
ID=86635417
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202211504485.9A Pending CN116246272A (zh) | 2022-11-28 | 2022-11-28 | 针对国产卫星多光谱图像质量标记的云雪区分方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN116246272A (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117173578A (zh) * | 2023-11-01 | 2023-12-05 | 长江水利委员会长江科学院 | 雪线高程检测方法、装置、计算机设备及存储介质 |
-
2022
- 2022-11-28 CN CN202211504485.9A patent/CN116246272A/zh active Pending
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117173578A (zh) * | 2023-11-01 | 2023-12-05 | 长江水利委员会长江科学院 | 雪线高程检测方法、装置、计算机设备及存储介质 |
CN117173578B (zh) * | 2023-11-01 | 2024-02-06 | 长江水利委员会长江科学院 | 雪线高程检测方法、装置、计算机设备及存储介质 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Shahtahmassebi et al. | Review of shadow detection and de-shadowing methods in remote sensing | |
Hirschmugl et al. | Single tree detection in very high resolution remote sensing data | |
CN109919875B (zh) | 一种高时频遥感图像特征辅助的居民地提取与分类方法 | |
CN110765934B (zh) | 一种多源数据融合的地质灾害识别方法 | |
CN108596103A (zh) | 基于最佳光谱指数选择的高分辨率卫星遥感影像建筑物提取方法 | |
CN111242224A (zh) | 一种基于无人机提取分类样本点的多源遥感数据分类方法 | |
CN113537018B (zh) | 基于多时相卫星遥感和无人机技术的水土保持监测方法 | |
Cai et al. | Study on shadow detection method on high resolution remote sensing image based on HIS space transformation and NDVI index | |
Li et al. | A comparative analysis of index-based methods for impervious surface mapping using multiseasonal Sentinel-2 satellite data | |
CN103839267A (zh) | 一种基于形态学建筑物指数的建筑物提取方法 | |
CN110988909A (zh) | 基于tls进行高寒脆弱区沙地植被的植被盖度测定方法 | |
CN107688777A (zh) | 一种协同多源遥感影像的城市绿地提取方法 | |
Chen et al. | A mathematical morphology-based multi-level filter of LiDAR data for generating DTMs | |
CN114022783A (zh) | 基于卫星图像的水土保持生态功能遥感监测方法和装置 | |
CN109472804A (zh) | 基于遥感影像的陆表水体提取方法和装置 | |
CN116246272A (zh) | 针对国产卫星多光谱图像质量标记的云雪区分方法 | |
Rashidi et al. | Extraction of ground points from LiDAR data based on slope and progressive window thresholding (SPWT) | |
CN116229459A (zh) | 国产卫星多光谱图像逐像素质量标记方法 | |
CN115760885B (zh) | 基于消费级无人机影像的高郁闭度湿地森林参数提取方法 | |
CN105681677B (zh) | 一种高分辨率光学遥感卫星相机最佳焦面确定方法 | |
CN115294183A (zh) | 一种基于多源遥感数据的碟型子湖水体时序提取方法 | |
CN112906528B (zh) | 一种融合多源卫星遥感数据的城市建筑物材质分类方法 | |
Grigillo et al. | Classification based building detection from GeoEye-1 images | |
KR102373278B1 (ko) | 지상기반 전천 영상자료를 이용한 야간 전운량 산출방법 | |
de Carvalho Jr et al. | Normalization of multi-temporal images using a new change detection method based on the spectral classifier |
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 |