一种面向高分五号高光谱影像的云检测算法

田静国, 范磊*, 于峻川, 容俊, 黄非, 崔艳芳

(1.航天宏图信息技术股份有限公司, 北京 100095; 2.中国自然资源航空物探遥感中心, 北京 100083)

高光谱光学影像成像过程易受到云的干扰,使得高光谱遥感数据中存在大量云覆盖区域,影响了遥感影像的质量,从而降低了影像的数据利用率。去云已成为影像预处理的一个重要步骤,而其中云检测是高光谱影像预处理的首要前提[1-6]。

目前,遥感影像云检测的方法大体可以分为纹理分析法、机器学习方法和阈值法。纹理分析法利用云的形状纹理进行云检测,如利用灰度共生矩阵、分形维数[7-8]等,但该方法仅仅使用纹理特征,很少涉及云的物理遥感特征,其在高光谱数据应用中的稳定性较低、精度有限。机器学习方法主要利用聚类分析或者监督分类的方法进行云检测,该方法通过已有的样本数据综合云的光谱、空间特征,构建特定的分类器,如Liu等[9]利用随机森林对葵花8号卫星影像云检测,提高了相对官方算法提高了5%的精度;Shao等[10]建立神经网络分类器,能够实现后厚云和薄云的识别;Cilli等[11]采用支持向量机(support vector machine,SVM)结合云样本,生成云检测模型实现云检测;Li等[12]、Segal-Rozenhaimer等[13]基于深度学习的方法云检测,并取得了较高的精度,但是该方法自动化程度较低,需要人为不断更新样本,进行大量样本训练用以提高检测精度,且该算法计算复杂,不利于快速提取和工程化应用。阈值法是利用云在可见光近红外波段的高反射和热红外波段的低温等物理特征,划定阈值进行云检测,如Huang等[14]以林覆盖区域为背景,通过简单阈值条件提取出无云像素点作为真实地面参考,并构建“光谱-温度”空间,根据云层与森林覆盖地表在该空间的差异,提出了一种适用于森林覆盖区域的云检测方法[14];Luo等[15]利用蓝(blue)、红(red)、近红外(near infrared,NIR)以及短波红外波段的简单阈值,实现云检测;Zhu等[16]提出掩码函数(function of mask,Fmask)算法,针对云层不同的物理、光谱和温度特征构建云概率影像,实现云检测,但是这些方法主要基于表观反射率或地表反射率,需要对原始数字量化值(digital number,DN)数据先进行辐射定标、大气校正等复杂的预处理步骤,不便工程化运用,且很少结合云的视觉色彩特征,云检测精度相对较低。

目前,以上研究方法多是面向多光谱影像,而高光谱相对于多光谱拥有高光谱分辨率、图谱合一等突出特点,上述方法在一定程度上不适用于高光谱影像。作为中国首颗全谱段高光谱卫星的GF5[17],该卫星的高光谱数据被应用于农业遥感、大环境监测等多个领域,但是针对GF5的云检测,特别是利用视觉特征结合光谱特征的云检测较少,故现综合利用云的视觉色彩特征和光谱特征,提出一种面向GF5原始高光谱影像处理简单、准确度高、便于工程化应用的高光谱云检测算法。

1 材料与方法

1.1 GF5高光谱数据

高分五号卫星(GF5卫星)于2018年5月9日发射成功,其上搭载的高光谱相机光谱分辨率为5 nm和10 nm,光谱范围为400～2 500 nm,是世界上首颗实现对大气和陆地综合观测的全谱段高光谱卫星,可满足环境综合监测等方面的迫切需求,是中国实现高光谱分辨率对地观测能力的重要标志[17]。GF5轨道标称技术指标见表1,可见短波红外高光谱相机(advanced hyperspectral imager,AHSI)技术指标见表2。

表2 可见短波红外高光谱相机指标

1.2 云检测算法

首先获得高光谱数字量化值(DN)影像,根据云的物理和视觉特性优选红绿蓝近红外等波段;根据优选波段的DN影像构建色调-饱和度-明度(hue-saturation-value,HSV)色彩空间,计算归一化植被指数(NDVI)、归一化水分指数(NDWI)、归一化雪指数(NDSI)、归一化建筑指数(NDBI)、霾优化变换(HOT)等参量;基于HSV色彩空间,分析云的视觉色彩特征,划定疑似云的色彩空间范围,并提取疑似云,获得疑似云影像和疑似非云影像;利用计算的NDVI、NDWI、NDSI等归一化指数以及云的物理特征,排除植被、水体、土壤、雪等地物的干扰,并构建概率影像;以视觉特征构建的疑似非云掩膜影像为基础,结合概率影像,两景影像取交集获得疑似非云概率影像和疑似云概率影像;通过自适应阈值划定算法,判定云分割阈值,超过该阈值的像元为云像元,从而得到云影像;对云影像进行云区域整合和优化,通过形态学和过滤的方法,去除碎斑块并填补细小空洞,最终得到云检测成果。具体流程图如图1所示。

1.2.1 云检测波段优选

GF5高光谱影像共有330波段,波段数据众多,需要对云的光谱特征进行分析,优选出对云敏感的波段析[18-24]。基于GF5高光谱的原始数字量值(DN)影像进行云光谱特征分析。由表2可知,高光谱一级数据的量化位数为12 bit,所以DN的最大值为4 096,将典型地物DN除以4 096,对DN进行归一化处理,云与植被、水体、雪、裸土地物DN归一化曲线分布如图2所示。

由图2可知,其中在400～900、1 200～1 300、1 500～1 700、2 100～2 200 nm云与其他地物对比最为明显。在400～900 nm为可见光-近红外波段,不同地物的反射特征不同,其中水体和植被DN值较低,裸土和雪其次,云的DN值最高,在此光谱范围内,晴空下各个地物具有较低的DN,DN归一化值在0～0.4,但云的DN值在400～700 nm都高于0.9,甚至接近与1,明显高于其他地物的DN。在1 200～1 300 nm范围内,水体的DN最低,接近为0,其次为植被、裸土、雪,DN在0.2左右,云的DN最高,在0.35以上。1 500～1 700 nm光谱范围内,水体DN最低,接近0,其次为植被,在0.15左右,再次为裸土和雪,阈值范围为0.2～0.3,雪的DN最高,在0.3以上。2 100～2 200 nm范围内,水体最低,植被在0.1左右,裸土在0.2左右,雪在0.3左右,雪最高,在0.4左右。

根据云和地物的DN光谱曲线分析,优选出480、560、660、860、1 250、1 600、2 200 nm 6个对云敏感的波段,用于云检测。

1.2.2 HSV色彩空间

HSV是一种将RGB色彩模型中的点在圆柱坐标系中进行展示的方法,这种表示法比基于笛卡尔坐标系的几何结构RGB更加直观,更能量化视觉特征,有利于云的识别[25-26]。

在色调-饱和度-明度(HSV)中,色调(H)是色彩的基本属性,用于量化颜色色彩的种类,取值范围为0°～360°,饱和度(S)是指色彩的纯度,越高色彩越纯,低则逐渐变灰,取0～100%的数值,明度(V)取值0～100%。

HSV可由RGB(即高光谱影像的670、560、460 nm的归一化DN)波段值转化构建,具体公式如下。

(1)

(2)

V=max

(3)

R=R′/212

(4)

G=G′/212

(5)

B=B′/212

(6)

式中:R′、G′、B′为高光谱影像中波长为670、560、460 nm处的DN;max为R、G、B中的最大值;min为R、G、B中的最小值。

云的高光谱影像的RGB真彩色合成和HSV假彩色合成如图3所示。

图3 GF5高光谱影像HSV合成影像以及云影像

由图3可知,HSV色彩空间中云呈现紫色、蓝色,薄云呈现为黄色、红色。其他地物如植被、裸土、城镇等呈现亮绿色、深绿色和浅蓝色,云阴影呈现黑色。经过大量的对比实验,归纳GF5高光谱云的HSV阈值范围为0°

1.2.3 云检测遥感指数

选取6种遥感指数,分别为NDVI、NDWI、NDSI、NDBI、HOT指数和SR。具体公式如表3所示。其中,NDVI是归一化植被指数,可用于检测植被生长状态、植被覆盖度和消除部分辐射误差等,云的NDVI多为负值,而植被的一半为正值,且随覆盖度增大而增大;NDWI为归一化水指数,主要用于提取水体,一般的水体的NDWI为正值,NDWI值越大代表是水体的可能性越高;NDSI为归一化雪指数,该指数基于积雪光谱特征提出的,主要应用与积雪检测和提取;NDBI为归一化建筑指数,用于建筑物提取;HOT指数用于表征雾霾的严重程度,值越大,雾霾相对越厚。

表3 云检测遥感指数

1.2.4 云概率影像

根据NDVI、NDWI等归一化参量以及云的光谱特征,排除植被、水体、建筑、土壤、雪等地物在云检测过程中的干扰,排除植被、水体、土壤的云概率为P1,排除雪、冰、建筑物的云概率P2,并构建面向云检测的概率影像Pcloud。示例影像见图4。

图4 GF5高光谱云概率影像

图5 基于数学形态学的云区优化整合

当NDVI<0.8且NDWI<0.8时,P1为

P1=ρ2 200ρ1 250

(7)

式(7)中:ρ2 200、ρ1 250分别为DN的归一化值。NDVI和NDWI见表3。

当HOT>0时,P2为

P2=HOT|(1-NDSI)|SR|(1-NDBI)|

(8)

式(8)中:HOT、NDSI、NDBI、SR计算公式见表3,然后,根据每个像元的P1、P2值计算每个像元的云概率像元值,进而得到整幅影像的云概率像元值,即云概率影像。具体计算公式为

Pcloud=P1P2

(9)

1.2.5 云检测动态阈值分割

云概率阈值是基于影像本身的动态云概率阈值进行云检测,即将云概率影像中的大于某一个阈值的像元,判定为云。根据影像的HSV色彩空间获取的疑似非云层区域M2,再叠加上Pcloud云概率影像分析,通过大量统计,得到云检测阈值为非层云区域的云概率82.5%。

1.2.6 云区域优化整合

通过云检测阈值分割之后,会存在大量碎斑噪声和云层间隙。需要在阈值分割的基础上,尽可能消除云检测噪声、填补云间隙,并优化云的轮廓,使云分布连通和分布合理。故采用过滤处理和数学形态算法。

过滤处理是解决分类图像中出现的孤岛问题。过滤处理使用斑点分组方法来消除这些被隔离的分类像元。类别筛选方法通过分析周围的4个或8个像元,判定一个像元是否与周围的像元同组。如果一类中被分析的像元数少于输入的阈值,这些像元就会被从该类中删除。

数学形态算法是指以形态为基础的一种数学分析工具,可用于视觉图像分析,常用于二值图像的分析和处理。主要使用了腐蚀、膨胀2种数学形态算法[33]。

对GF5高光谱影像云检测结果的优化整合步骤如下。

步骤1在阈值分割的结果图6(a)的基础上,对按照3×3的尺度的腐蚀(对应GF5的影像为9个像素),消除结果如图6(b)所示。

图6 高光谱影像云检测结果

步骤2在上一步的基础,进行5×5尺度(对的GF5的影像为25个像素)的膨胀,结果如图6(c)所示。

步骤3在上一步基础上,进行过滤处理,过滤128个像素以下的碎斑,填补128个像素以下的空隙,结果如图6(d)所示。

1.2.7 云区域优化整合

采用查准率(p)和查全率(r)两个指标对云检查结果进行评价[34]。公式如下。

(10)

(11)

式中:Pt为实际是云,并且能够被正确识别为云的像元数目;Pf为不是云,但被误识别为云的像元数目;Pn为实际是云,但是没有识别出来是云的像元数目。

2 结果与分析

2.1 云检测结果展示

选取典型地区、不同时期的多景GF5影像,按照本文算法进行云检测,结果如图6和图7所示。

图7 判定无云的情况

由图6可知,算法对森林、戈壁、农田、冰山等场景下有很好的适用性,能够很好地将云和冰/雪区别,但是在检测薄云和碎云方面,有一定的漏检,在一些高亮地物有个别错检的情况。如图7所示为被判定为无云的案例。对于不含云的场景如水域、植被、冰雪等地区中,在无明显云特征的情况下,算法可成功地将此类影像判定为无云。

2.2 定量评价

选取5组数据对算法进行定量评价,这5组数据包含了森林、高原、平原、沙漠、城镇等多个场景,覆盖西北地区、青藏地区、东北地区、华北地区、南方地区多个区域,每组数据约30景GF5影像,可以有效代表中国区的典型地貌,评价结果具有普遍代表性。

在评价过程中,使用查准率和查全率两个指标,云层分布的真实是由人工目视解译后勾画出来的。为了保证算法运行速度的客观性,算法选取以每景数据的中心像元为中心的1 000×1 000的影像块,进行算法运行和定量评价。

算法的运行环境为一台CPU为AMD 4800H (主频2.9 GHz),运行内存在16 GB, 512 GB固态的个人电脑,操作系统为Windows10 家庭版。算法程序由基于GDAL和numpy的python3.8开发而成。算法的输入为GF5-AHSI的VN 和SW 原始DN 值tiff影像以及人工勾画的云矢量文件,经过波段筛选、无效值去除、构建HSV色彩空间、构建云概率影像、动态阈值分割、云区域形态学优化整合、精度评价等步骤后,生成云掩膜文件,格式为geotiff,云检测矢量文件,格式为shapefile。

定量评价统计表如表4所示。

表4 高光谱影像云检测算法精度评价

2.3 结果分析

通过对不同场景、不同地区、不同时间的云检测结果进行分析,如图6和图7以及表4所示,本文算法对于含云大多数场景,云检测结果精度都比较高,云提取的轮廓准确,云检测查准率在87.22%～98.17%,查全率在82.85%～98.34%;算法能够有效排除雪/冰的影响,算法精度受云形态的影响较小;算法执行效率高,对单景的高光谱运行速度约为25 s。但是该算法在云斑面积较小时会出现漏检现象;在薄云区域,由于下垫面光谱信息占主导地位,云的光谱信息较少,会出现漏检;对于少量高亮的区域,土壤、建筑有误检的情况。

3 结论

综合以上分析,得到以下结论。

(1)提出的色彩空间特征和光谱特征相结合的云检测算法,能够对不同场景下GF5高光谱影像实现云检测,检测准确率在87%以上,可有效排除雪、冰、植被等地物的干扰。

(2)算法执行速度快,且不需要对影像数据做复杂的预处理操作,便于工程化应用。