GF-4序列图像的云自动检测[PDF全文]
胡昌苗, 白洋, 唐娉
摘要: 以高分四号(GF-4)卫星L1级标准分幅数据产品提供高精度的云检测产品为目的,研究针对地球同步轨道卫星数据的云检测算法,改进自动阈值以适应同日不同时刻成像数据的辐射亮度与地表反射特性的变化差异。利用GF-4卫星凝视成像方式获取的同区域序列图像以及云在不同图像上的运动特性,结合自动阈值与Savitzky-Golay(SG)滤波修正检测结果中的误检。算法的两个关键预处理,一是通过自动的几何配准解决未经几何校正的分幅数据之间像素位置对应的问题,二是通过基于典型相关变换自动提取序列图像之间的伪不变特征点集,进而利用相对辐射归一减小了不同时刻成像数据之间的辐射差异。通过内蒙古自治区东部及长江中下游区域70余组数据对算法进行验证,整体上获得了稳定的结果与精度,并且基于序列图像的云检测算法在云边界、高亮地表及薄云区域的检测精度整体优于单幅自动阈值的检测结果。结果表明算法精度上满足GF-4云检测数据产品需求,且算法自动化程度高,便于工程化的数据生产。
关键词: GF-4     云检测     自动匹配     Savitzky-Golay滤波    
DOI: 10.11834/jrs.20186401    
收稿日期: 2016-11-11
中图分类号: TP701    文献标识码: A    
作者简介: 胡昌苗(1981— ),男,助理研究员,研究方向为遥感辐射传输及大气校正理论和应用。E-mail:hucm@radi.ac.cn
基金项目: 国家自然科学基金(编号:41601384);高分辨率对地观测系统重大专项(编号:11-Y20A05-9001-15/16,03-Y20A04-9001-15/16)
Automatic cloud detection for GF-4 series images
HU Changmiao, BAI Yang, TANG Ping
1. Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing 100094, China
2. University of Chinese Academy of Sciences, Beijing 100094, China
Abstract: The research on cloud detection is an important branch of remote sensing image research. In recent years, with the increasing number and type of remote sensing satellites, cloud detection based on reference map/sequence image has become a subject receiving close review in cloud detection. GF-4 is a geo-synchronous orbit satellite launched by China in December 2015. This satellite is equipped with a visible-light/near-infrared camera with a resolution of 50 m and has typical high-resolution and multi-spectral satellite data characteristics. GF-4 satellites have many common characteristics as meteorological satellites. These common features include the geostationary orbit, area array starring imaging, and the mid-infrared band. GF-4 has the capability to acquire the sequence data of the same area in a short time. This paper attached importance to the algorithm of automatic cloud detection for early GF-4 satellite data acquisition. The research is based on the same area of multiple GF-4 images. First, according to the characteristics of the image data and cloud in different images on the movement characteristics, this work performed automatic geometric registration and relative radiation normalization on multiple images. Then, the image was set to automatic threshold by cloud detection and was processed by the Savitzky–Golay filtering. Finally, this work implemented an automatic cloud detection algorithm for GF-4 sequence images. This paper selected 36 data in the eastern region of Inner Mongolia and 39 data in the middle and lower reaches of the Yangtze River for cloud testing to detect the practical feasibility of the new algorithm. The following preliminary conclusions were obtained. (1) The results of the cloud detection algorithm based on sequence image are superior to those of the single-image cloud detection in terms of overall accuracy. The main difference was observed in the image of the cloud boundary, highlighted surface, and thin cloud area. Through the experiment, this work showed that the algorithm has a high degree of automation and can satisfy the needs of engineering data. (2) Based on the single-image cloud detection, using the automatic threshold method can provide an overall stability for the test results. However, owing to the diversity of the cloud in the image, improving the accuracy is obviously difficult for the proposed algorithm. (3) The mid-infrared band data of GF-4 cannot be simply used for cloud detection due to differences in coverage area, spatial resolution, and acquisition time. The shortcoming of the current algorithm is that the acquisition time of the sequence data is extremely long. Eliminating the radiation difference between the data obtained in the morning and those obtained in the noon with the simple linear relation is difficult. Simultaneously, the follow-up research will focus on the systematic cloud detection accuracy evaluation method.
Key Words: GF-4    cloud detection    automatically matching    Savitzky-Golay filtering    
1、引 言

云在遥感图像中具有较高的辐射亮度,并且遮蔽地物,影响遥感数据的应用。对云像元进行标记,生成云掩膜被很多卫星数据产品所采用。云的研究一直是遥感领域的重点,尤其是在气象卫星数据研究领域(Stowe 等,1999Walder和Maclaren,2000Heidinger 等,2002),气象卫星通常在设计时便考虑到对云的监测,比如采用相对静止的地球同步轨道、相机波段设置在对云敏感的非大气窗口区域、且空间分辨率不高等。气象卫星云检测的常用算法有阈值法、统计法及基于辐射传输的方法等(Goodman和Henderson-Sellers,1988)。针对风云系列气象卫星,马芳等人(2007)综合利用多个波段进行云检测。文雄飞等人(2009)对FY-2C提出了一种基于云指数的云检测方法。刘健(2010)采用滑动分析区和嵌套分析区的方法改进FY-2C云检测的阈值提取方法。王根等人(2015)提出了基于最小剩余法的FY-3B云检测方法。卢晶等人(2015)提出了基于FY-3C卫星双氧通道的云检测方法。相对气象卫星云检测,高分辨率地表监测卫星由于波段设置、分辨率差异等因素,云检测方法更多样,Sedano等人(2011)利用同时相、同区域的多幅高分辨率遥感图像检测云,并通过基于目标的区域生长来修正云腌膜;Martinuzzi等人(2007)利用结合地物先验特征的自动云覆盖估计法对Landsat图像进行云检测。从模式识别角度,金炜等人(2010)提出一种基于密度聚类支持向量机的云检测方法。闫宇松和龙腾(2010)通过提取多种特征,采用最小距离分类器,实现实时的云检测;赵敏等人(2012)提出了基于支持向量机与无监督聚类相结合的云分类算法。基于单幅图像的云检测方法是目前应用主要的方法,算法稳定且易于业务化运行,缺点是容易出现误检与漏检的情况,尤其是高亮地表、积雪等。

近年来,随着遥感卫星数量与种类的快速增长,遥感图像获取变得越来越容易,基于参考底图/序列图像的云检测得到发展,该类方法利用云在不同图像上的运动特性,有效减少了高亮地表的误检。目前利用不同高分辨率遥感卫星数据构建的多套中国区域多季相近无云底图都可用于对应卫星数据的云检测。Savitzky-Golay(SG)滤波是一种常用的算法(Savitzky和Golay,1964),对序列图像进行滤波的同时还可去除云与云下阴影(辜智慧,2003Chen 等,2004)。尽管基于序列图像的处理可以获得更高的精度,但在高分辨率遥感卫星云检测的实际业务运行初期中却很少采用,原因是积累足够的数据构成参考底图或者序列图像通常是滞后于云检测业务处理,当传感器积累的数据足够多时才开始采用。

高分四号卫星(GF-4)是中国2015年12月发射的一颗地球同步轨道卫星,搭载空间分辨率为50 m的可见光近红外相机,地物细节特征丰富,具有典型的高分辨率多光谱卫星数据特性,同时GF-4卫星又具有气象卫星的很多特征,包括采用地球同步轨道、面阵凝视成像、具有中红外波段及短时间内获取同一区域的序列图像的能力等。GF-4数据的上述特征,使得其云检测研究具有挑战性。结合GF-4数据特性,研究卫星数据获取初期适用的自动云检测算法是本文研究的出发点。

基于GF-4数据特征,对于同区域的多幅图像,利用云在不同图像上的运动特性,在自动的几何配准与辐射相对归一的基础上,结合自动阈值与Savitzky-Golay滤波实现了一种GF-4序列图像的自动云检测算法。算法自动化程度高,可以满足工程化处理数据的需求。通过内蒙古东部区域36幅数据及长江中下游区域39幅数据对云检测算法进行了验证,结果表明利用GF-4卫星数据的序列图像来提高单幅图像云检测的精度是有效的。

2、数据和方法     (2.1) GF-4数据分析

GF-4由中国航天科技集团公司空间技术研究院(航天五院)研制,属于国家高分辨率对地观测系统重大专项天基系统中的系列卫星,设计寿命8年,于2015年12月29日发射入轨。GF-4运行在距地36000 km的地球同步轨道,定点位置105.6°E,采用面阵凝视方式成像,具备全色、多光谱和中波红外成像能力,有高时间分辨率和较高空间分辨率的优势。通过指向控制,GF-4可实现对中国及周边地区进行自由观测,也可采用凝视模式对幅宽400 km的固定区域进行持续观测。GF-4卫星载荷相机参数见表1

由于采用地球同步轨道,GF-4卫星具备短时间内获取同一地理区域多幅数据的能力。但是受相机硬件约束,GF-4无法像气象卫星一样24 h连续成像,相机每天工作时间平均为2 h,凝视模式下可对同一区域连续成像20景左右。GF-4在实际的灾害监测应用中,对灾害区域的多天定时连续观测数据通常具有较强的序列特性。对于可见光近红外数据,形成序列的图像可以有两种方式:(1) GF-4相机在凝视工作模式下获取的同一天、相近时间内获取的连续多幅图像组成的序列;(2) 非凝视模式下,不同天获取的、同一地区的多幅图像。在GF-4卫星的实际运营中,获取第2种非凝视模式构成的序列图像数据很容易。因为灾害监控作为GF-4的重要任务,当某一区域发生灾害时,GF-4可以在经纬度误差0.1°的范围内对固定区域进行多天、多次观测。GF-4数据发布的L1级数据包含的头文件中记录辐射定标参数及成像几何信息。L1级数据未经系统几何校正,但包含系统几何校正的有理函数参数文件。

表 1 GF-4载荷相机技术参数 Table 1 GF-4 payload technical indicators
    (2.2) GF-4云检测算法流程

云检测算法针对可见光近红外数据,算法流程主要包含两个步骤,分别为单幅图像云检测与序列图像云检测,具体流程图如图1

图 1 GF-4云检测算法流程 Figure 1 Flow chart of GF-4 cloud detection algorithm
    (2.3) 单幅图像云检测

阈值法是目前确定云盖量的常用方法。简单的固定阈值法被用于很多太阳同步轨道的卫星遥感数据,根据传感器不同增益状态设置经验阈值。利用波段差异增强云特征的云指数也常被用于确定阈值。复杂些的自动阈值法则通过分析直方图的基础上添加各种判断等方法获得相对精确的云阈值。

GF-4单幅图像云检测同样基于阈值法, 是根据GF-4数据特征定制的自动阈值法。利用可见光近红外各波段云和背景场之间的数值差别,反差越大越容易识别云。第2波段(0.52—0.60 μm)云也有较高反射率,但受大气散射影响明显,薄雾较多,云边界模糊,不利于云检测。第3波段(0.52—0.60 μm)绿色植被最低,水体其次,土壤和城镇反射率最高,云的反射率明显高于下垫面,并随着厚度、高度而变化。第4波段(0.63—0.69 μm)晴空数据具有较低的反射率,厚云则有高的反射率,同时该波段能够增强云和陆地之间的对比,是云检测的重要波段。第5波段(0.76—0.90 μm)云也有较高反射率,可是一些裸露地表、贫瘠土地以及植被也会有和云相似的反射率,难以单独进行云检测。用第3与第4波段进行基于自动阈值的云检测,策略是分别统计波段3与波段4的直方图,有云图像的直方图假定有两个波峰,根据高反射率所在的统计波峰确定两个波段云的自动阈值所在的大致范围,以该波峰为起点,向低反射率方向,以步长5个灰度值为单位,统计步长内像素点个数的增长率。当增长率大于一定阈值a时(比如a=0.2%),认为到达云和地表的分界点。由于云的波段性差异,波段3和波段4各自通过上述方法获取的阈值并不相同,进一步的修正是在两个阈值基础上小范围内寻找中间值,以两个波段统计的云覆盖百分比最接近的点作为最终单幅云检测的阈值。这样,确定阈值的关键参数只有a,调整a的数值大小可改变整体上检测出的云量的多少,以修正不同时期数据的检测结果。

由于GF-4是地球同步轨道卫星,数据获取在白天的任意时段,辐射值变化差异大。虽然GF-4的L1级数据包含定标参数,通过辐射预处理转换可得表观反射率或地表反射率,但并不适用固定阈值的云检测。同样原因,构造的云指数数值也是不稳定的。

自动阈值法简单稳定,满足大多数数据的云检测精度质量要求,但云在遥感图像中的表现是极端复杂的,对于检测结果不理想的数据,想要进一步提高精度,将自动阈值方法与其他方法相结合是获得更加精确的云检测结果的有效途径。例如结合多阈值启发式和人工神经网络对GMS卫星红外云图中尺度云系进行自动分割(师春香和瞿建华,2001),采用从纹理特征的进行云检测的方法(郁文霞 等,2006单娜 等,2009)等。本文选取了两种修正阈值检测结果的方法。

第1种修正是简单的利用形态学的方法进行云区整饰。用于处理个别数据云与非云的边界处产生大量的碎片,如图2(a)所示,通过去除云边界外小于一定像素数的孤立云区,并填充云边界内小于一定像素数的空洞整饰云边界,效果如图2(b)。云区整饰优化云掩膜的质量的同时大大减少了云区数量,减少了后续处理的计算量。需要注意的是,由于云的复杂性,实际的云分布是过度离散的,云区整饰反而会降低云检测的精度。

图 2 GF-4云区整饰示例 Figure 2 An example of GF-4 cloud area finishing

第2种修正是基于纹理的云检测方法。利用分形维数来区分高亮地表与云。分形维数反映物体纹理的复杂程度(Takayasu,1990),维数越大,分形越复杂。该法假设分辨率50 m GF-4数据地表纹理分形维数数值偏大;云层纹理细节的分形维数数值偏小。

具体方法是把2维图像视作3维空间中的一个表面 $\left( {x,y,f\left( {x,y} \right)} \right)$ ,其中 $f\left( {x,y} \right)$ 为图像 $\left( {x,y} \right)$ 位置处的灰度值,将M×N大小的图像按照尺度r分割成许多r×r的网格。在每个网格上,是一列r×r×h的盒子,单个盒子的高度h=r×G/M,其中G为总的灰度级数。设在第 $\left( {i,j} \right)$ 网格中图像灰度的最小值和最大值分别落在第k和第l个盒子中,则覆盖第 $\left( {i,j} \right)$ 网格中的图像灰度值所需的盒子数为 ${n_r}(i,j) = l - k + 1$ ,则整个图像所需的盒子数为 ${N_r} = \sum\limits_{i,j} {{n_r}(i,j)} $ ,分形维数值为 $D = \lim \displaystyle\frac{{\lg ({N_r})}}{{\lg (l/r)}}$

选取一组r以lg(l/r)为横坐标,lg(Nr) 为纵坐标绘制点聚图。利用最小二乘法拟合样本点,所得斜率即为此图像样本分形维数D

分形维数计算选取64×64大小的子图作为基本的处理单元,选择4、8、16等3个尺度因子(单娜 等,2009)。分形的计算量远大于普通的基于灰度共生矩阵的纹理计算,但可以得到对地物类别更加敏感的纹理图,有助于区分云、雪与高亮地表。

    (2.4) 可见光近红外序列图像云检测

GF-4成像的一个重要特征是可以获得相同地区同一天或者连续多天的序列图像。序列图像云检测利用云在不同图像上的运动特性来与地表区分。序列图像云检测的算法流程如图3

图 3 GF-4可见光近红外序列图像云检测算法流程 Figure 3 Flow chart of GF-4 visible light & near infrared series images cloud detection

可见光近红外序列图像的自动匹配是算法流程的关键技术。由于GF-4卫星具有静止卫星的面阵成像方式、良好的姿态稳定性等硬件条件保证了成像几何特性的一致性,使得L1级的序列图像之间通过简单的平移与旋转即可实现高精度的逐像素配准。以序列图像其中一幅图像作为基准,获取其余图像到该基准图像的平移与旋转关系实现对应,无需对图像本身进行变换与重采样。

序列图像自动匹配3个步骤:(1) SIFT (Scale-invariant Feature Transform)特征提取(Lowe, 2004, 2005);(2)特征匹配;(3)图像配准关系计算。

在特征匹配过程中,利用L1级数据里的有理函数模型,可确定待匹配的两幅图像之间大致的地理位置关系,这样便可以将每一个SIFT特征点的匹配限定在一个局部区域里,提高匹配速度(Shan等,2014)。在两幅图像上分别取模板窗口 $\left( {M \times M} \right)$ 和搜索窗口 $\left( {N \times N,N > M,M = 65,N = 165} \right)$ ,让模板窗口在搜索窗口内移动,计算每一个位置与搜索窗口的相似度,如果相似度最大值大于某一阈值T(T=0.65),则匹配成功。相似度判别准则采用归一化相关系数(NCC),具体计算采用NCC的加速算法Fast NCC(Lewis,1995)。特征匹配完成后,采用RANSAC (RANdom SAmple Consensus)算法剔除误匹配点。考虑到GF-4卫星L1级数据图像单幅像素尺寸达10240×10240,若整幅图像进行SIFT特征提取将会消耗大量的运算时间与资源,所以采用分块处理的策略,将整幅图像分成1024×1024的小块进行逐块自动匹配,单块匹配成功则停止块循环,节省匹配时间。

在配准关系计算过程中,利用两幅图像的特征匹配点集,采用最小二乘法拟合一阶多项式y=ax+b的旋转参数a与平移参数b

序列图像自动匹配全部完成后,序列图像中的像素便通过旋转与平移参数实现对应。之后,对于凝视模式成像的序列图像可以直接进行云检测。而对于非凝视模式的序列图像,由于图像成像时间的差异、相机成像状态不同(主要是成像积分时间的差异)等因素,使得各个序列图像的像素值存在整体上的差异。利用线性相对辐射归一技术消除这种差异,即利用一阶多项式y=ax+b描述两幅图像之间的辐射差异,将辐射差异分解为增益a与偏移b。以非凝视模式下序列图像的其中一幅作为基准图像,其余图像用线性相对辐射归一到该基准图像。拟合参数ab所需的伪不变特征地物点通过迭代加权的多源变化检测IR-MAD (Multivariate Alteration Detection)算法自动提取,拟合变换则采用正交变换(Hu和Tang,2011)。

根据序列图像数量的多少,采用不同的云检测方法,当序列数量小于10幅时,采用分级自动阈值进行云检测,当序列数量大于10幅时,采用SG滤波结合自动阈值的方法进行云检测。

分级自动阈值云检测的原理是利用云的运动特性,序列图像中某些像素的像素值只在部分图像上较高,而在其余图像上不高,这些便是受云影响的像素,将高像素值所在图像的位置处标记为云。所有图像中像素值均明显偏高的像素则判定为高亮地表。判定每个像素点是否为云的依据有两个,1是每幅图像中根据单幅云检测的自动阈值方法计算得到的云阈值,以及统计得到的可能的阴影区域阈值,2是该像素与对应位置所有图像的像素,排除每个图像中大于该图云阈值与小于该图像阴影阈值的像素后,像素中值的差值,记为a。单幅图像云检测结果中被检测为云的像素若a值小,且与该图云阈值接近,则认为是误检测的高亮地表,单幅图像云检测结果中被检测为地表的像素若a值大,且排除阴影像元,则将该像素识别为云,并根据差值大小在云阈值区间的位置,标记为薄云与厚云。

SG滤波通过滑动窗口多项式拟合来达到对序列数据进行平滑的目的(Savitzky和Golay,1964)。序列数为N,对其中长度为n=2m+1的子序列进行 $k\left( {k \leqslant n} \right)$ 阶多项式拟合可表示为

$f(t) = {a_0} + {a_1}{t^1} + \cdots + {a_k}{t^k} = \sum\limits_{i = 0}^{i = k} {{a_i}{t^i},} t \in [ - m,m]$

SG滤波过程是对序列中的某一点t0及其左右m邻域共n=2m+1个点 $\left( {{t_i},{y_i}} \right)$ $i \in \left[ { - m,m} \right]$ ,进行k $\left( {k \leqslant n} \right)$ 的多项式拟合,用拟合后的滑动窗口中心的数据 $\left( {{t_0},{y_{0'}}} \right)$ 置换原始时间序列中的数据 $\left( {{t_0},{y_0}} \right)$ ,然后向右移动窗口,使窗口中心移至序列中下一数据,重复上述过程,直到滑动窗口到达序列末尾。平滑窗口系数通过最小二乘求得。通过SG滤波便可以判断序列中每一个图像中的像素值在序列变化趋势中的异常,结合分级自动阈值后可以更加有效地区分运动的云区与不变的高亮地表,并标记出薄云与厚云。

由于分级自动阈值与SG滤波都是逐像元的处理,可能导致局部区域的云掩膜出现很多由多个像素组成的漏洞或者碎片,通过云区整饰提高最终云掩膜结果的成图质量。

    (2.5) 中波红外在云检测中的局限性

云在中波红外波段中的低温特性可用于云检测,但GF-4的中波红外数据在实际应用中有局限性:

首先,中波红外与可见光近红外的分辨率差异过大,一个400 m分辨率的中波红外像素对应一个8×8像素块、共64个50 m分辨率的可见光近红外像素。图像尺寸也不同,可见光近红外为10240×10240像素,中波红外为1204×1024像素。数据质量上,中波红外还存在明显的噪声,如图4中矩形框放大图。

其次,中波红外与可见光近红外之间的几何配准困难。两者覆盖地理范围不完全重合,图5为可见光近红外与中波红外根据头文件中的系统几何模型,按照地理坐标配准叠加显示的结果,A为中波红外,B为可见光近红外,可见中波红外数据的覆盖范围比可见光近红外要小,并且两者的地图中心位置也存在差异。有理函数模型无法实现地理位置的相对配准,图5中C与D分别为箭头白框处中波红外于可见光近红外局部放大数据,点a与点c为相同地物点,点b为可见光近红外上点a在中波红外上的地图坐标点,c与b之间的差异便是地理位置误差,需要进一步的修正,而两者的数据差异使得实现自动配准困难,难以适应工程化快速处理的需求。

最后,中波红外与可见光近红外成像时间存在差异。由于中波红外与可见光近红外被设计为成像时间固定间隔45 s。云的运动特性使得中波红外与可见光近红外数据中云的位置、形态发生变化,如图4(a)(b)所示,仔细比较可发现云的细节差异。

基于以上原因未将中波红外数据用于云检测。

图 4 可见光近红外与中波红外云区示例 Figure 4 GF-4 cloud in visual near infrared and medium wave infrared bands

图 5 GF-4可见光近红外与中波红外数据示例(A白色虚线框为中波红外,B为可见光近红外,C&D为白框中波红外与可见光近红外局部,点a与点c为相同地物点,点b为可见光近红外上点a在中波红外上的地图坐标点) Figure 5 GF-4 bands in visual near-infrared and medium wave infrared (medium wave infrared band in white dashed box A, visual near-infrared wave band B, C&D is a partial area of white box in A&B, a and c is the same point of earth’s surface, a and b is the same point of in visual near-infrared and medium wave infrared bands)
3、数据结果处理与分析     (3.1) GF-4实验数据

实验数据采用36幅内蒙古自治区东部区域数据及39幅长江中下游区域数据。受“厄尔尼诺现象”影响,2016年夏季中国内蒙古自治区区域和长江中下游区域夏季出现了典型灾害性气候。民政部国家减灾中心针对7月、8月份内蒙古东部发生的干旱现象,以及长江中下游进入主汛期出现的洪涝灾害,制定了受灾区域GF-4数据的获取计划,体现了GF-4在灾害监测方面的价值。长江中下游区域数据,获取时间为2016-06-17—2016-08-24,经纬度范围为111.3°E—115.9°E,28.3°N—31.0°N。内蒙古区域数据,获取时间为2016-05-10—2016-08-15,经纬度范围为114.6°E—121.7°E,44.0°N—45.5°N。实验数据云量统计见表2表3,其中数据编号根据获取时间排序,其中字体加黑并下划线的编号表示可以构成序列的图像。构成序列图像的条件是至少存在两幅图像,图像之间4个角点的经纬度坐标差异≤0.1°。内蒙古自治区东部区域有18幅图像构成多组序列,占总数的50%。长江中下游数据可以构成序列的有28幅,占总数的72%。可见GF-4获取数据具有明显的序列特性。

表 2 内蒙古东部区域GF-4 36景影像数据云量统计 Table 2 GF-4 data cloud cover statistics in east of neimenggu province (36 scenes)

表 3 长江中下游区域GF-4 39景影像数据云量统计 Table 3 GF-4 data cloud cover statistics in the middle and lower regions of the changjiang river (39 scenes)
    (3.2) 单幅图像云检测实验

单幅图像云检测,首先对实验数据进行自动阈值法云检测,然后进行云区整饰与基于纹理的修正。

自动阈值法云检测结果主要通过目视进行评价。由于云在图像中的表现极端复杂,目前还难以找到一种定量化的GF-4云检测精度比较算法。比较了多种基于统计的与构造云指数的检测方法后,最终形成的自动阈值云检测方法针对实验数据获得的结果稳定性最好,选取的阈值整体上最优,可满足业务化单幅云检测数据生产要求。

云区整饰与基于纹理的修正目的是提高单幅图像自动阈值云检测的精度,但实验结果中对精度的提升非常有限。实验数据中出现局部检测结果过离散的情况不多,其中极少区域是由于薄云边界过渡范围大导致,云区整饰是有效的,剩余的多是极小的碎云,云区整饰反而将面积只有几个像素的碎云去除了。纹理修正希望能去除部分高亮地表的误检区域,在部分实验数据中取得了很好的效果,但同时也将很多纹理信息丰富的云识别为地表。

对75幅数据的结果比较,整体上自动阈值结果与修正后结果差异很小。图6是单幅云检测修正前后云覆盖率比较图,纵轴为自动阈值结果,横轴为修正后结果,两种结果云盖率差异小。图7是单幅云检测示例,图7(a)上方局部区域有小范围极小的碎云,在云区整饰的过程中被错误删除了,同时图7(a)下方局部区域,一些误检测的细小高亮地表则会在云区整饰的过程中被修正。其中纹理信息强的高亮地表,纹理修正取得一定的效果,但亮度极高且纹理弱的地表仍然误识别为云。

通过单幅图像云检测的实验,在夏季干旱与汛期两组图像中云的复杂特性得到体现,认为基于自动阈值的云检测方法取得的结果适应性与稳定性最好,可用于业务化运行。而云区整饰与基于纹理的修正对复杂云的检测表现不稳定,对精度的提升不明显,可用于对某些类别的误检测结果进行修正。

图 6 单幅云检测修正前后云覆盖率比较图 Figure 6 Comparison of cloud coverage for single image’s cloud detection

图 7 GF-4单幅云检测实例 Figure 7 An example of GF-4 single image’s cloud detection
    (3.3) 序列图像云检测实验

序列图像云检测,先进行自动匹配与线性相对辐射归一等预处理,然后进行序列图像云检测验证。

自动配准是序列图像云检测预处理的关键。实验数据中最长序列数为凝视模式下获取的9幅图像,表4为该序列其中的8幅到另外1幅的配准参数,将拟合的一阶多项式系数转换成旋转与平移系数,可发现所有图像的旋转量都接近0值,说明序列图像可通过像素平移实现逐像素配准,这是GF-4凝视模式成像的特点,非凝视模式下、不同天获取的序列图像之间几何偏差一般大于凝视模式。需要说明的是,该基于仿射变换的自动配准只适用于L1级数据,而利用数据自带的有理函数模型进行系统几何校正后,序列数据的配准关系会变复杂。

表 4 GF-4序列配准参数示例 Table 4 An example of registration parameters for a GF-4 series images

线性相对辐射归一对于非凝视模式获取的序列图像是必要的预处理步骤。由于GF-4采用地球同步轨道,同一区域不同天获取的数据获取时刻差异大,辐射亮度差异大,传感器增益不同,输出值不同。实验表明,线性相对辐射归一可有效降低辐射差异,保障后续的云检测的稳定。局限性是由于不同的地表类型在不同时刻的反射特性差异,线性相对辐射归一减小了大部分的地表类型的辐射差异,对于个别辐射变化大的地物效果有限。图8为内蒙古东部局部区域,两幅获取时间分别为早上8点与中午12点数据云检测结果,图8中的部分水体及周围由于在不同时刻的辐射差异大于整体辐射差异,而被误检测为云。

图 8 数据获取时刻差异的影响示例 Figure 8 An example of effects at differet acquisition time

序列图像云检测实验多采用分级自动阈值,随着GF-4数据量的增加,结合SG滤波能得到更好的云检测结果。序列图像云检测结果保存为8位图像,0值表示无云区域,数值1为单幅云检测结果,数值2与3为基于序列检测后的云区,云掩膜数据示例如图9。序列图像云检测结果一般多于单幅云检测,多出的部分以过渡型云边界部分为主(图10)。序列图像云检测在云区边缘检测精度高,云区边缘亮度相对较低,利用基于单幅图像阈值检测难以区分,尤其是在薄云区域,亮度上还包含部分地表的贡献,亮度较低。薄云的检测主要依赖地表参照,通过序列图像可检测出。序列图像云检测可以减小高亮地表造成的误检,如图11,高亮的水体与薄云区域通过序列分析得以区分。

精度验证方面,通过目视评价,序列图像云检测结果整体上优于单幅云检测的结果,主要体现在云边界、高亮地表及薄云区域。从所有序列云检测数据中挑选出18幅高质量的结果数据,并人工逐幅修正了其中的误检区域后,将其作为检验数据,验证单幅云检测的精度。方法是利用像素统计,统计对应数据的单幅云检测结果中检测正确的云区像素数占检验数据中所有云区像素数量的百分比,18幅图像的最低值为78%,最大值93%,平均值为88%。

图 9 序列图像云掩膜标记示例(黑色0值为无云,深灰色1值为单幅云检测结果,浅灰色2值为厚云,白色3值为薄云) Figure 9 An example of cloud mask’s value (black 0 value is cloudless, dark gray 1 value is the results of single image’s cloud detection, light gray 2 value for the thick cloud, and white 3 value for the thin cloud)

图 10 序列图像云检测示例 Figure 10 An example of series images’ cloud detection

图 11 序列图像云检测示例 Figure 11 An example of series images’ cloud detection
4、结 论

本文针对GF-4地球同步轨道卫星数据,研究生产L1级数据云掩膜产品的问题,提出了一种序列图像云检测算法,通过内蒙古东部及长江中下游区域的多组数据实验,初步得到以下规律性的认识:

(1) 基于序列图像的云检测算法结果精度整体上优于基于单幅图像云检测,主要体现在云边界、高亮地表及薄云区域。

(2) 基于单幅图像云检测采用自动阈值法可以获得整体稳定的检测结果,由于云在图像中表现出的复杂性,复杂的修正算法对精度的提升有限。

(3) GF-4中波红外数据由于覆盖区域、空间分辨率、获取时间等差异难以用于云检测。

目前算法存在的不足主要是获取时刻差异过大的序列数据上面,早上与正午获取的数据之间的辐射差异难以利用简单的线性关系消除,造成云检测精度变差。另外配准与相对辐射校正过程占据大部分算法耗时,算法有待优化。

GF-4云检测算法仅用于数据获取初期,在其8年的设计寿命期间,随着数据的不断获取,积累的不同季相、不同时相、不同时刻的历史数据将会降低云检测的难度的同时提高检测精度。系统的云检测精度评价方法将是后续研究的重点。

参考文献
[1] Chen J, Jönsson P, Tamura M, Gu Z H, Matsushita B, Eklundh L. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter[J]. Remote Sensing of Environment, 2004, 91 (3/4) : 332 –344. DOI: 10.1016/j.rse.2004.03.014
[2] Goodman A H, Henderson-Sellers A. Cloud detection and analysis: A review of recent progress[J]. Atmospheric Research, 1988, 21 (3/4) : 203 –228. DOI: 10.1016/0169-8095(88)90027-0
[3] 辜智慧. 2003. 中国农作物复种指数的遥感估算方法研究——基于SPOT/VGT多时相NDVI遥感数据. 北京: 北京师范大学 Gu Z H. 2003. A study of calculating multiple cropping index of crop in China using SPOT/VGT multi-temporal NDVI Data. Beijing: Beijing Normal University
[4] Heidinger A K, Anne V R, Dean C. Using MODIS to estimate cloud contamination of the AVHRR data record[J]. Journal of Atmospheric and Oceanic Technology, 2002, 19 (5) : 586 –601. DOI: 10.1175/1520-0426(2002)019<0586:UMTECC>2.0.CO;2
[5] Hu C M, Tang P. Automatic algorithm for relative radiometric normalization of data obtained from Landsat TM and HJ-1A/B charge-coupled device sensors[J]. Journal of Applied Remote Sensing, 2011, 6 (1) : 063509 . DOI: 10.1117/1.JRS.6.063509
[6] 金炜, 俞建定, 符冉迪, 岑雄鹰, 尹曹谦. 利用密度聚类支持向量机的气象云图云检测[J]. 光电子•激光, 2010, 21 (7) : 1079 –1082. Jin W, Yu J D, Fu R D, Cen X Y, Yin C Q. Meteorological imagery cloud detection using density clustering support vector machine[J]. Journal of Optoelectronics•Laser, 2010, 21 (7) : 1079 –1082.
[7] Lewis J P. 1995. Fast normalized cross-correlation. http://scribblethink.org/Work/nvisionInterface/nip.pdf[2016-11-11]
[8] 刘健. FY-2云检测中动态阈值提取技术改进方法研究[J]. 红外与毫米波学报, 2010, 29 (4) : 288 –292. Liu J. Improvement of dynamic threshold value extraction technic in FY-2 cloud detection[J]. Journal of Infrared and Millimeter Waves, 2010, 29 (4) : 288 –292.
[9] Lowe D G. Distinctive image features from scale-invariant keypoints[J]. International Journal of Computer Vision, 2004, 60 (2) : 91 –110. DOI: 10.1023/B:VISI.0000029664.99615.94
[10] Lowe D. 2005. Demo code for detecting and matching SIFT features, Version 4, July 6, 2005, http://www.nexoncn.com/read/3470e251dac51e3cf7289acb.html
[11] 卢晶, 薛胜军, 韩阳, 张夏琨, 孙晓娟, 王志伟, 田伟, 杨润芝. 基于风云3C卫星双氧通道的云检测算法[J]. 科学技术与工程, 2015, 15 (2) : 179 –182. Lu J, Xue S J, Han Y, Zhang X K, Sun X J, Wang Z W, Tian W, Yang R Z. A cloud detection algorithm based on FY-3C satellite double O2 channels [J]. Science Technology and Engineering, 2015, 15 (2) : 179 –182.
[12] 马芳, 张强, 郭妮, 张杰. 多通道卫星云图云检测方法的研究[J]. 大气科学, 2007, 31 (1) : 119 –128. Ma F, Zhang Q, Guo N, Zhang J. The study of cloud detection with multi-channel data of satellite[J]. Chinese Journal of Atmospheric Sciences, 2007, 31 (1) : 119 –128. DOI: 10.3878/j.issn.1006-9895.2007.01.12
[13] Martinuzzi S, Gould W A and Ramos Gonzalez O M. 2007. Creating cloud-free landsat ETM+ data sets in tropical landscapes: cloud and cloud-shadow removal. General Technical Report IITF-GTR-32, United States Department of Agriculture: 1–18 [DOI: 10.2737/IITF-GTR-32]
[14] Savitzky A, Golay M J E. Smoothing and differentiation of data by simplified least squares procedures[J]. Analytical Chemistry, 1964, 36 (8) : 1627 –1639. DOI: 10.1021/ac60214a047
[15] Sedano F, Kempeneers P, Strobl P, Kucera J, Vogt P, Seebach L, San-Miguel-Ayanz J. A cloud mask methodology for high resolution remote sensing data combining information from high and medium resolution optical sensor[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66 (5) : 588 –596. DOI: 10.1016/j.isprsjprs.2011.03.005
[16] 单娜, 郑天垚, 王贞松. 快速高准确度云检测算法及其应用[J]. 遥感学报, 2009, 13 (6) : 1138 –1155. Shan N, Zheng T Y, Wang Z S. High-speed and high-accuracy algorithm for cloud detection and its application[J]. Journal of Remote Sensing, 2009, 13 (6) : 1138 –1155. DOI: 10.3321/j.issn:1007-4619.2009.06.012
[17] Shan X J, Tang P, Hu C M. An automatic geometric precision correction system based on hierarchical registration for HJ-1 A/B CCD images[J]. International Journal of Remote Sensing, 2014, 35 (20) : 7154 –7178. DOI: 10.1080/01431161.2014.967884
[18] 师春香, 瞿建华. 用神经网络方法对NOAA-AVHRR资料进行云客观分类[J]. 气象学报, 2002, 60 (2) : 250 –255. Shi C X, Qu J H. Cloud classification for NOAA-AVHRR data by using a neural network[J]. Acta Meteorologica Sinica, 2002, 60 (2) : 250 –255. DOI: 10.11676/qxxb2002.031
[19] Stowe L L, Davis P A, McClain E P. Scientific basis and initial evaluation of the CLAVR-1 global clear/cloud classification algorithm for the advanced very high resolution radiometer[J]. Journal of Atmospheric and Oceanic Technology, 1999, 16 (6) : 656 –681. DOI: 10.1175/1520-0426(1999)016<0656:SBAIEO>2.0.CO;2
[20] Takayasu H. 1990. Fractals in the physical sciences. Manchester: Manchester University Press
[21] Walder P, Maclaren I. Neural network based methods for cloud classification on AVHRR images[J]. International Journal of Remote Sensing, 2000, 21 (8) : 1693 –1708. DOI: 10.1080/014311600209977
[22] 王根, 华连生, 刘惠兰, 张苗苗. 基于最小剩余法的FY-3B/IRAS资料云检测研究[J]. 红外, 2015, 36 (9) : 15 –20, 29. Wang G, Hua L S, Liu H L, Zhang M M. Study of FY-3B/IRAS data cloud detection based on minimum residual method[J]. Infrared, 2015, 36 (9) : 15 –20, 29.
[23] 文雄飞, 董新奕, 刘良明. " 云指数法”云检测研究[J]. 武汉大学学报(信息科学版), 2009, 34 (7) : 838 –841. Wen X F, Dong X Y, Liu L M. Cloud index method for cloud detection[J]. Geomatics and Information Science of Wuhan University, 2009, 34 (7) : 838 –841.
[24] 闫宇松, 龙腾. 遥感图像的实时云判技术[J]. 北京理工大学学报, 2010, 30 (7) : 817 –821. Yan Y S, Long T. Real-time cloud detection in optical remote sensing image[J]. Transactions of Beijing Institute of Technology, 2010, 30 (7) : 817 –821.
[25] 郁文霞, 曹晓光, 徐琳, BencherkeiM. 遥感图像云自动检测[J]. 仪器仪表学报, 2006, 27 (6S) : 2184 –2186. Yu W X, Cao X G, Xu L, Bencherkei M. Automatic cloud detection for remote sensing image[J]. Chinese Journal of Scientific Instrument, 2006, 27 (6S) : 2184 –2186.
[26] 赵敏, 张荣, 尹东, 王奎. 一种新的可见光遥感图像云判别算法[J]. 遥感技术与应用, 2012, 27 (1) : 106 –110. Zhao M, Zhang R, Yin D, Wang K. Cloud classification algorithm for optical remote sensing image[J]. Remote Sensing Technology and Application, 2012, 27 (1) : 106 –110. DOI: 10.11873/j.issn.1004-0323.2012.1.106