基于高分四号数据的多时相多通道云检测算法

王纪辉，李峰，鹿明， 马骏，郭毅

1. 河南大学 软件学院，开封 475100 2. 中国空间技术研究院 钱学森空间技术实验室，北京 100094 3. 澳大利亚西悉尼大学，悉尼 NSW 2150

1 引言

高分四号(GF-4)卫星是一颗地球同步轨道高分辨率的光学遥感卫星，配置了一个全色多光谱传感器(panchromatic multispectral sensor，PMS)和一个中红外传感器(infrared sensor，IRS)。PMS传感器的空间分辨率为50 m，重访时间为20 s，单幅成像幅宽超过500 km。IRS传感器的空间分辨为400 m，重访时间为20 s，单幅成像幅宽超过400 km[2]。GF-4可以对地面进行连续观测，提供高时间分辨率的连续图像，为观测地面移动目标变化提供有利条件，在减灾、气象、林业、地震等方面发挥着重要作用[3]。

在利用遥感图像获取地表信息的过程中，云会阻碍地表与传感器之间的辐射能量传输，导致获取到地表信息存在误差，从而降低数据的利用率。因此为了提高遥感图像数据的利用率和精度，云的检测和去除是非常重要的,遥感图像的云检测一直是一个非常活跃的研究领域[4]。

云检测方法大致上可以分为两类。一类为基于单时相遥感图像云检测方法，其中包括基于光谱分析的方法，利用云与典型地表的光谱特征差异来进行云检测。文献[5]利用单时相的高分一号数据来进行云检测，但是该方法仅考虑了裸土和植被这两种下垫面没有考虑冰雪的情况。文献[4]利用单时相的GF-4多光谱影像来进行云检测，但是该方法对于薄云的检测能力较差且较难区分冰雪与云像元。还有基于深度学习的云检测方法，卷积神经网络(convolutional neural networks，CNN)、生成式对抗神经网络(generative adversarial networks，GAN)等深度学习算法被应用于遥感图像云检测领域。文献[6]将影像的特征信息输入到支持向量机分类器进行分类，从而获得资源三号卫星影像的云检测结果，但是该方法步骤较多、网络复杂且无法实现像素级的云检测。另一类为基于多时相遥感图像云检测方法，这种方法通过计算多张不同时间获取的卫星影像之间的差异，将这些图像之间光谱特征发生显著变化的像元标记为云像元，来实现遥感图像的云检测。文献[7]提出了一种多时相多通道阈值结合的风云四号气象卫星云检测方法，但是GF-4影像的波段较少，该方法不适用于GF-4影像的云检测。文献[8]提出了一种不依赖于热红外波段增强型多时相云检测算法，解决了单时相云检测中难以区分高亮地表和云的问题，有效提高了云检测的精度，但是该方法需要一景无云的影像作为参考，存在一定的局限性，文献[2]提出一种多时相RTD云检测算法来实现GF-4影像的云检测，但是该算法在下垫面为冰雪时云检测精度较低。

云与像元在可见光波段的表观反射率相近，仅依靠可见光波段难以区分云与冰雪像元。而与其他地球同步轨道卫星相比，GF-4卫星影像的波段较少，只有全色、红、绿、蓝、近红外和中红外这6个波段。这导致了单时相云检测算法由于缺乏热红外波段难以检测出冰雪像元，而多时相云检测算法又依赖于云的移动来区分云与冰雪像元，这很大程度上影响了云检测的精度，进而影响了GF-4影像的下一步利用率和精度。

根据GF-4卫星可以提供高时间分辨率连续的可见光影像和中红外影像[9]这一特征，本文提出了一种多时相多通道云检测算法。该算法利用云与典型地表的光谱特征分析来识别出潜在云像元，使用多时相影像像元之间的光谱差异来识别出移动的云像元，利用中红外波段反演出地表的亮度温度识别出冰雪像元。通过大量试验和验证结果表明该算法对GF-4号卫星影像具有较高的云检测精度，特别是在区分云与冰雪像元的方面优于其他算法。

2 研究区与研究数据

2.1 研究区介绍

为了验证算法的适用性，选取了3个研究区域，研究区域一为中国海南省及其相邻海域，位于106°44′(E)～111°37′(E)，17°42′(N)～22°31′(N)，该区域包括海面、山地、河道丰富的地表类型；选择2017-07-27，当日薄云较多。研究区域二为中国辽宁省和吉林省部分地区，位于120°14′(E)～127°21′(E)，37°38′(N)～44°10′(N)，该区域的覆盖地表类型大多为积雪覆盖的裸土，还包括森林和结冰的河流、冰雪等多种地表类型；选择2017-01-17，当日大面积厚云覆盖于冰雪表面。研究区域三为中国安徽省和江苏省部分地区，位于115°68′(E)～121°51′(E)，29°47′(N)～35°26′(N)，该区域包含大面积积雪覆盖的城区；选择2017-01-08，当日大面积薄雪覆盖在城区。丰富的地表覆盖类型可以充分验证本文算法对于不同下垫面云的检测精度。

2.2 研究数据

本文采用的研究数据为GF-4卫星全色多光谱传感器获取的50m分辨率的多光谱影像和中红外传感器获取的400 m分辨率的中红外影像。GF-4卫星获取的多光谱影像和中红外影像的波段分布以及相关参数如表1所示[10]。其中，10组不同地区不同时间的多光谱和中红外影像作为参考影像包含裸土、水体、植被和冰雪等地表覆盖类型，用于确定光谱分析的阈值参数。每组有3张地点相同且时间为同一天的影像，用于多时相云检测。3组海南、辽宁和安徽研究区域的多光谱和中红外影像，作为测试影像用于验证本文云检测算法的精度。

表1 GF-4卫星传感器波段分布以及相关参数

3 云检测算法

云检测算法主要分为4个步骤：①对GF-4卫星多光谱影像和中红外影像进行预处理；②利用可见光近红外波段中云和其他典型地表的光谱差异识别出潜在云像元；③利用多时相可见光影像蓝光波段之间的差异识别出地表上移动的云像元；④利用中红外影像反演的地表亮度温度去除冰雪像元。GF-4卫星影像云检测算法流程如图1所示。

图1 GF-4卫星影像云检测算法流程Fig.1 GF-4 satellite image cloud detection algorithm flow

3.1 预处理

预处理模块主要分为4个步骤：辐射定标、多时相可见光波段配准、多光谱波段间配准以及多光谱和中红外波段配准。

(1)辐射定标

对GF-4号多光谱影像和中外影像进行辐射定标，将原始DN值转换为表观辐亮度[11]：

L=Gain×DN+offset

(1)

式中：L为表观辐亮度；DN为像元灰度值；Gain为定标斜率；offset为绝对定标系数偏移量。GF-4号卫星影像自带的XML文件中的辐射定标系数与实际值差别较大，可以从中国资源卫星应用中心(http:∥www.cresda.com/CN/)发布的GF-4绝对定标系数文件中获取得到。将上一步计算得到的表观辐亮度转换成表观反射率[12]。

(2)

式中：ρ为表观反射率；D为日地天文单位距离；ESUN为大气层顶的平均太阳光谱辐照度；θ为太阳天顶角，可以从GF-4卫星影像自带的XML文件中获得。

(2)多时相可见光波段配准

GF-4卫星采用的是面阵成像机制对地凝视观测[13]，使得L1级的相邻时间相同地点的多光谱影像通过简单的平移与旋转就能实现配准。采用基于相位相关和重采样的亚像素配准算法[14]来实现多时相可见光波段的亚像素级配准。首先基于相位相关实现像素级的粗定位，然后在粗定位点1.5倍的邻域范围内利用二维快速傅里叶变换实现100倍上采样，最后基于相位相关实现重采样范围内亚像素级的细定位。

(3)多光谱波段间配准

GF-4卫星的L1级产品虽然经过初步的几何校正，但是谱段间还是存在一定的偏差，本文采用多时相可见光波段配准的算法进行配准。

(4)多光谱和中红外波段配准

GF-4号卫星的多光谱影像的空间分辨率为50 m，而中红外影像的空间分辨率为400 m，因此，需要先将影像重采样至相同分辨率。中红外影像与多光谱影像是由不同的传感器拍摄的，在地理位置上还存在着偏差。使用GF-4影像自带的RPC文件对处理后的GF-4卫星多光谱影像与中红外影像进行波段合成[15]，使影像重投影到同一个投影坐标，并使用3次卷积内插法将中红外影像重采样至50 m分辨率，生成一个具有6个波段的新影像，新生成影像的6个波段分别为全色、蓝光、绿光、红光、近红外、中红外，然后以第3个波段为参考图像，第6个波段为配准图像使用ENVI[16]的配准方法进行配准。

3.2 单时相影像识别潜在云像元

根据云与典型地表的光谱分析可以看出，云在可见光和近红外波段的表观反射率要高于大多数典型地表。因此，利用云与典型地表在可见光和近红外波段的光谱差异来识别潜在云像元。本文通过蓝光波段检测、“白度”指数检测、HOT检测、NDWI检测和NDVI检测这5种检测方法来检测潜在云像元，只有同时满足以上5种检测方法的像元才会被认定为潜在云像元。

(1)蓝光波段检测

云在蓝光波段具有很高的表观反射率，一般都高于其他典型地表。因此，可以根据蓝光波段的表观反射率来识别基于蓝光波段检测的潜在云像元，检测方法为：

Rband2>0.15

(3)

式中：Rband2为蓝光波段的表观反射率。

(2)“白度”指数检测

部分地物在蓝光波段也有较高的反射率，但是云相较于其他典型地物在红、绿、蓝3个可见光波段的反射率变化趋于平缓，因此，文献[17]提出了“白度”指数测试的方法用于云检测。通过“白度”指数来识别基于“白度”指数的潜在云像元，检测方法如下：

Raver=(Rband2+Rband3+Rband4)÷3

(4)

(5)

式中：Rband3为绿光波段表观反射率；Rband4为红光波段表观反射率；Raver为红、绿、蓝光3个波段表观反射率的均值。

(3)HOT检测

雾霾最优化转换(haze optimized transformation，HOT)最初是用来检测卫星影像中的雾霾，在可见光波段中大多数晴空像元的蓝光与红光波段的表观反射率存在线性关系，当受到雾霾影响时线性关系会发生变化[18]。因此，文献[19]将HOT检测应用于区分云与雾霾，取得了很好的效果。通过HOT值来识别基于HOT检测的潜在云像元，检测方法为：

Rband2-0.5×Rband4>0.15

(6)

(4)NDWI检测

归一化水指数(normalized difference water index，NDWI)一般是用来提取影像中的水体信息，水体在绿光波段的表观反射率较高，而在近红外波段却随着波长的增加而下降[20]。因此，水体的NDWI值一般较高，通过NDWI值来识别基于NDWI检测的潜在云像元，检测方法为：

(7)

式中：Rband5为近红外波段的表观反射率。

(5)NDVI检测

归一化植被指数(normalized difference vegetation index，NDVI)一般是用来检测卫星影像中的植被覆盖情况[21]。植被在红光波段的表观反射率较低，而在近红外波段却迅速增加。因此，可以通过NDVI值来识别基于NDVI检测的潜在云像元，检测方法为：

(8)

3.3 多时相影像识别移动的云像元

GF-4卫星获取多光谱影像的重访时间为20 s，因此，便于获取同一天相邻时间内同一地点的多幅影像。GF-4卫星影像的空间分辨率为50 m，在GF-4多光谱影像的重访时间内可以合理地假设地表没有发生改变。在蓝光波段云的表观反射率一般高于其他典型地表，而云阴影的表观反射率明显低于大部分的典型地表，当云发生移动的时候云阴影也会随之发生移动。因此，当同一位置像元的表观反射率发生较大增加时，可以认定此像元为移动的云像元或者云阴影像元。通过3幅时间相邻预处理配准后的影像之间的蓝光波段表观反射率差值来检测移动的云像元并且去除云阴影像元，检测方法如下：

(9)

(10)

3.4 利用中红外影像去除冰雪像元

在可见光近红外波段，冰雪与云的表观反射率无论是在数值上还是在变化的趋势上都非常接近，仅利用可见光近红外波段很难去除冰雪像元。文献[22]利用绿光和短波红外的表观反射率计算归一化雪盖指数(normalized difference snow index，NDSI)来去除冰雪像元，而GF-4卫星仅有可见光近红外和中红外波段，缺少短波红外和热红外波段。本文通过GF-4中红外影像来反演地表的亮度温度，根据亮度温度设立阈值来区分云与冰雪像元。亮度温度(brightness temperature，BT),是与传感器观测地表物体辐射出射度相等的黑体温度，在数值上等于辐射温度[23]。其计算方法为：

(11)

式中：h=6.626 2×10-34J·s为普朗克常数；T为亮度温度；c=2.997 93×108m/s为光速；k=1.380 6×10-23J/K为玻尔兹曼常数；λ为波长。

3.5 云检测算法的阈值

固定阈值很难精准识别云像元，这是由于云像元的表观反射率会受到地表类型的影响，特别是薄云像元，不同的地表类型表观反射率会发生很大变化。本文将10组参考影像裁切成大小为512×512像素的下垫面分别为裸土、水体、植被和冰雪的影像，每种下垫面各10张影像，通过试验验证确定了如表2所示的云检测阈值。

表2 不同地表类型的云检测阈值

3.6 云检测算法的精度评价

为了对云检测的精度精确评估，本文采用了对比试验和定量分析的方法。

(1)传统单时相云检测算法

本文采用的传统单时相云检测算法是一种基于光谱分析的阈值云检测算法[4]，基于云与典型地表的光谱特征，使用光谱差异检测出潜在云像元，然后根据有云像元和无云像元的光谱变化率差异计算云概率。

(2)SVM算法

支持向量机(support vector machine，SVM)算法是一种高精度的机器学习分类算法[24]，适用于云检测领域。在实现SVM算法的过程中，使用本文算法相同的参考图像获取SVM算法最佳参数。通过测试线性核函数、多项式核函数、径向核函数和sigmod核函数的云分类效果，核函数选择径向核函数、核函数系数选择0.2、惩罚参数选择100时云分类效果最佳。

(3)RTD算法

实时差分(real time difference，RTD)算法是一种高精度的多时相云检测算法[2]，通过计算一对实时影像之间的差异来获得移动的云，以便去除高亮地表和冰雪像元。

(4)定量分析

为了对云检测的精度进行定量分析，本文采用云检测领域通常用到的准确识别率、误检率和漏检率[25]来对结果进行量化分析。

(12)

(13)

(14)

式中：correct为云像元的准确识别率，true为影像中被正确检测为云像元的个数，miss为影像中云像元被检测为晴空像元的个数，commission为云像元的误检率，false为影像中晴空像元被检测为云像元的个数，omission为云像元的漏检率。

4 结果分析

为了验证本文算法的适用性，选择了海南、辽宁和安徽3个不同的研究区域进行实验，采取对比试验的方式验证算法的正确性，将本文算法的检测结果与其他3个云检测算法进行对比并进行定量分析，其中包括传统单时相云检测算法、SVM算法， RTD云检测算法。

4.1 海南研究区域

海南研究区域有丰富的水体和薄云，便于测试本文算法对水体上的薄云的检测能力，检测结果如图2所示。

图2 海南研究区域云检测结果Fig.2 Cloud detection results in Hainan study area

图2(a)为海南研究区域截取大小为2048×2048像素的红、绿、蓝波段组成的RGB影像，图2(b)为预处理后的中红外影像，图2(c)(d)(e)分别为单时相算法云检测结果和SVM算法云检测结果和RTD算法云检测结果,图2(e)为多时相多通道算法云检测结果。为了更好地评价算法的适用性，在可见光影像中随机选取了2 000个像素点对算法进行定量分析。4种算法的定量分析结果如表3所示。通过分析结果对比发现，准确识别率都达到了90%以上，但其他3种算法漏检率比多时相多通道算法高，这是由于在可见光波段水体表面薄云的表观反射率与水体差别不大，导致其他3种算法没有将部分水体表面的薄云检测出来。而RTD算法的误检率相较于本文算法低，这是由于RTD算法通过多时相云检测将海岸线边缘的高亮地表检测出来，而本文算法将少部分高亮地表检测为云像元。

表3 海南研究区域云检测结果定量分析

4.2 辽宁研究区域

辽宁研究区域中含有大面积冰雪覆盖的裸土和森林以及结冰的河流，便于检测本文算法去除冰雪的能力，检测结果如图3所示。

图3 辽宁研究区域云检测结果Fig.3 Cloud detection results in Liaoning study area

从检测结果可以看出多时相多通道算法能很好地区分出云与冰雪像元，而单时相算法却将部分被积雪覆盖的地表以及结冰的河道误判成云，这是由于在可见光近红外波段冰雪与云的表观反射率差别不大，而SVM算法和RTD算法将部分裸土检测为云，部分薄云没有检测出来。对辽宁研究区域进行定量分析，分析结果如表4所示。通过对分析结果对比发现，多时相多通道算法的云检测正确率高于其他3种算法，而且误检率也低于其他3种算法。

表4 辽宁研究区域云检测结果定量分析

4.3 安徽研究区域

安徽研究区域中包括由大量被积雪覆盖的城区和部分覆盖在积雪上面的云，便于验证本文算法检测冰雪表面云的能力，检测结果如图4所示。

图4 安徽研究区域云检测结果Fig.4 Cloud detection results in Anhui study area

从检测结果上看多时相多通道算法能够很好地检测出积雪表面的云像元，而其他3种算法却将部分积雪覆盖的城区检测为云像元，这是由于两者在可见光近红外波段的表观反射率相近。SVM算法将积雪表面覆盖薄云像元检测为冰雪像元，导致部分薄云没有检测出来。对安徽研究区域进行定量分析，分析结果如表5所示，通过分析结果对比发现，多时相多通道算法的识别准确率高于单时相算法和SVM算法，而且误检率和漏检率均低于单时相算法和SVM算法。

表5 安徽研究区域云检测结果定量分析

5 结论

GF-4卫星是中国第一个高分辨率地球同步轨道卫星，GF-4卫星影像波段较少且缺少热红外波段，导致云检测过程中去除冰雪像元面临着较大困难。本研究根据GF-4卫星的特点，提出一种多时相多通道的云检测算法，采用光谱分析的方法[26]来检测影像中的潜在云像元，并采用多时相数据光谱差异的方法来检测影像中移动的云像元，结合中红外影像来去除影像中的冰雪像元。该算法在海南、辽宁和安徽3个研究区域进行了测试，并与传统单时相云检测算法、SVM算法和RTD算法的检测结果进行对比，结果表明，本文算法在不同下垫面的云检测准确识别率均优于其他3种算法，特别是对于薄云和冰雪像元有更好的检测精度。

虽然本文算法在GF-4卫星影像的云检测方面取得了较高的精度，但是也有一定的局限性。一方面是GF-4的多光谱影像和中红外影像是由不同传感器获取的，两者的分辨率和几何位置不同，配准的难度较大，而配准的精度影响冰雪像元的去除，从而导致云检测精度的下降。另一方面是在高亮地表的检测方面误检率较高。未来可以考虑更加精确的配准方法来实现GF-4可见光波段和中红外波段高精度的配准以及降低高亮地表的误检率，进一步完善GF-4卫星影像的云检测。