1 引言

卫星影像在成像过程中会受到云层遮挡,这会对遥感分析、动态监测等研究过程产生重要影响^[1-3],同时这也是造成地表信息分类、变化检测以及其他定量信息提取精度降低的主要原因之一^[4-7]。因此,在利用遥感数据提取地表、海洋及大气信息之前进行云识别具有重要意义。

目前,已经发展了多种云检测方法,其中的固定阈值法是应用比较广泛、成熟的云检测方法,其原理是利用云与典型地物在可见光、近红外波段(NIR)的反射率差异及在热红外波段的亮度温度差异来识别云像元。该方法具有容易实现且效率和精度相对较高的特点^[8-9]。目前已经发展了多种典型的阈值法云检测方法,如:国际人造卫星云气候学计划(ISCCP)法^[10-12]、CLAVR (the NOAA Cloud Advanced Very High Resolution Radiometer) 法^[13]和APPOLLO (the AVHRR Processing scheme Over Clouds, Land and Ocean) 法^[14]。传统的阈值法主要利用云在可见光、近红外波段的反射率高于大部分典型地物反射率的特点,在不同波段设定固定的阈值来实现云与大部分地表的分离。但由于地表结构的复杂性以及云的多样性,很难找到合适的阈值将所有类型的地表和云分离开,特别是对容易受地表影响的薄云、碎云、云边缘以及沙漠、裸土、岩石等高反射率地表上空的云进行识别比较困难。当卫星传感器波段较少时,难以提供足够的信息来确定云与典型地表的辐射差异,导致云像元的识别精度较低。针对该问题,Sun等^[15]提出了一种先验地表反射率数据库支持的动态阈值云检测算法(UDTCDA),该算法以现有的高质量地表反射率产品作为支撑,辅助识别有云像元,其优点是:1)由于背景反射率已知,阈值的设定有了针对性,可以有效提高不同地表背景上空云(尤其是薄、碎云)的识别精度;2)仅利用少数可见光及近红外波段便可实现遥感影像的云识别,在波段较少的卫星传感器的云识别中具有独特优势。由于大部分高分辨率卫星传感器的波段较少,且波段集中在可见光到近红外等波长范围内,云检测可利用的波段范围有限。因此,UDTCDA在高分辨率影像云检测中具有显著优势。UDTCDA使用当前中分辨率成像光谱仪(MODIS)提供的地表反射率数据来支撑云识别,空间分辨率为500 m,这对于高空间分辨遥感数据来说较为粗糙,原有的算法是将待检测像元的空间分辨率重采样到与先验支撑数据一致以后再进行云检测。这种方法不仅降低了待检测图像的空间分辨率,使检测结果的空间范围与实际的云覆盖范围有较大差异,而且会造成像元信息丢失,产生较大的识别误差。本文根据高空间分辨率的特点,采用逐像元空间地理坐标配准的方法实现了高空间分辨率影像与支撑数据的配准,提高了该算法在高分辨率影像云检测中的应用。将该算法用于空间分辨率较高(2~50 m)的高分1号(GF-1)、高分2号(GF-2)、高分4号(GF-4)和资源3号(ZY3)多光谱影像进行云检测实验。

2 数据源与预处理

2.1 高分辨率卫星数据

本课题组选择空间分辨率较高的ZY3和高分系列卫星(GF-1、GF-2和GF-4)数据开展云识别实验。ZY3于2012年1月9日成功发射,是中国首个民用高分辨率光学传输型立体测图卫星,搭载有一台用于地物解译与应用研究的空间分辨率为6 m的多光谱相机,包含蓝、绿、红和近红外4个波段,重访周期为5 d,相机刈幅宽为51 km。GF-1、GF-2和GF-4是高分辨率对地观测系统的组成部分,它们的发射实现了亚米级高空间分辨率与高时间分辨率的有机结合^[16]。GF-1搭载了两台全色-多光谱相机(PMS1、PMS2),它们的空间分辨率分别为2 m和8 m,刈幅宽为60 km;4台多光谱中分辨率宽幅相机(WFV1、WFV2、WFV3、WFV4)的空间分辨率均为16 m,刈幅宽均为800 km。GF-2搭载的两台全色-多光谱相机(PMS1、PMS2)的空间分辨率分别为1 m和4 m,刈幅宽均为45 km。GF-4搭载了一台凝视相机,包含可见光-近红外波段(VNIR)和中红外波段(MWIR),空间分辨率分别为50 m和400 m,刈幅宽为400 km。

2.2 地表反射率数据集

本课题组使用美国国家航空航天局(NASA)提供的MOD09地表反射率产品(MOD09A1)构建先验地表反射率数据集支持云检测。MOD09A1是由MOD09逐日地表反射率产品合成得到的,可以在各个波段上较精确地估计地物光谱的反射率,整体精度为±(0.005+0.05ρ),其中ρ为地表反射率^[17-19]。本课题组获取2012年覆盖中国全境的所有MOD09A1数据后,分别对其进行几何校正、数据拼接和裁剪等基本数据预处理。为了降低云以及阴影的影响,采用次最小值合成方法合成得到包含蓝、绿、红和近红外4个波段的月合成地表反射率数据集^[20]。该地表反射率数据集的空间分辨率为500 m,统一投影为阿尔伯斯投影,坐标系为WGS84。图1为2012年10月覆盖中国全境的地表反射率假彩色合成图像(RGB:213),可以看出,该数据的空间连续性较好,云覆盖较少,整体质量较高。

图 1. 2012年10月地表反射率数据集假彩色合成图像

Fig. 1. False color composite image of surface reflectance in October, 2012

下载图片 查看所有图片

3 改进的UDTCDA

作为一种基于光谱反射率差异的识别算法,云识别算法的关键在于阈值的确定。UDTCDA利用可靠的先验地表反射率数据集为云检测提供了真实下垫面的地表反射率信息,然而卫星影像观测到的是表观反射率,由于大气-地表系统间复杂的相互作用,导致这两种反射率之间存在明显差异。研究表明,在特定的观测条件和大气环境背景下,卫星测得的表观反射率与真实的地表反射率具有明显的函数关系^[21-22]。因此,该算法综合考虑了卫星观测几何条件和气溶胶等大气背景条件等因素的影响,利用太阳光谱内卫星信号的第二代模拟(6S)模型模拟传感器地表反射率与表观反射率之间的量化关系,通过最小二乘法拟合得到地表反射率与表观反射率极大值之间的函数关系,构建出不同光谱通道的动态阈值云检测模型。该方法能够有效降低混合像元的影响,提高对不同地区上空薄、碎云的识别能力,降低高亮地表上云识别的不确定性,具有较高的云识别精度。

3.1 动态阈值模型的构建

本课题组首先利用6S模型模拟了不同观测和大气条件下,各高空间分辨率传感器在可见光到近红外光谱通道上的表观反射率和地表反射率之间的量化关系。如果待检测像元的表观反射率大于地表反射率在不同观测和大气条件下模拟的表观反射率的极大值,则认为该像元为潜在云像元。因此,表观反射率的极大值可以表示为地表反射率与观测几何的函数:

ρ*=aρ+bcosθυcosθs+c,(1)

式中:ρ^*为模拟得到的云的表观反射率阈值;ρ为地表反射率,由先验地表反射率计算得到;θ_υ和θ_s分别为太阳天顶角和卫星天顶角;a、b、c为模型中的三个系数,可通过最小二乘法拟合得到。只有在4个通道都被判定为潜在云像元时,该像元才最终被确定为云像元。表1给出了不同陆地观测卫星不同通道的动态阈值计算模型参数。

3.2 光谱响应误差的校正

光谱响应函数描述的是传感器某一通道对特定波谱区间的响应能力。随着影像空间分辨率的提高,地物的光谱特征更加丰富^[23]。图2列出了MODIS与GF-1、GF-2、GF-4和ZY3在可见光和近红外波段的光谱响应函数曲线,可以看出,MODIS的光谱响应函数与高分辨率卫星(GF-1、GF-2、GF-4、ZY3)的光谱响应函数在各个波段上的差异较为明显,这主要是由传感器本身设计和光谱通道设置的差异造成的。

为了降低光谱响应差异对云识别的影响,提高识别的整体精度,本课题组利用美国地质调查局(USGS)地物波谱库中植被和土壤的光谱数据对不同传感器间的光谱响应误差进行校正^[24-25]。本课题组共获取了29条实测地物的光谱数据,其中包括16条植被和13条土壤的光谱数据,根据各传感器的光谱参数构建MODIS与高分辨率卫星地表反射率之间的转化统计模型。回归分析结果表明,MODIS与高分辨率卫星地表反射率之间存在较强的线性关系:

ρ'=mρM+n,(2)

图 2. 国产卫星传感器与MODIS传感器光谱响应函数的对比。(a)红光波段中光谱响应函数对比;(b)绿光波段中光谱响应函数对比;(c)蓝光波段中光谱响应函数对比;(d)近红外波段中光谱响应函数对比

Fig. 2. Comparison of spectral response function between MODIS and different domestic satellite sensors. (a) The comparison of spectral response functions in red channel; (b) the comparison of spectral response functions in green channel; (c) the comparison of spectral response functions in blue channel; (d) the comparison of spectral response functions in near infrared (NIR) channel

下载图片 查看所有图片

式中:ρ'为校正后的地表反射率;ρ_M先验地表反射率;m和n分别为系数。表2给出了国产陆地观测系列卫星(GF-1 PMS、GF-1 WFV、GF-2、GF-4、ZY-3)与MODIS不同光谱通道的地表反射率转换参数。

3.3 像元匹配方法改进

UDTCDA识别结果的空间分辨率依赖于地表反射率数据集。由于待检测影像与地表反射率数据集具有不同的空间分辨率,因此,实现准确的影像匹配是提高云识别精度的重要一步。先验地表反射率数据集的空间分辨率为500 m,而本课题组用到的高分辨率数据则主要为4~50 m,它们在像元尺度上具有较大差距,若通过重采样将两幅影像的空间分辨率统一到500 m,就会导致高空间分辨率像元中的部分云信息损失,增加了识别难度。为了解决这一问题,本课题组根据高空间分辨率像元中心坐标与先验地表反射率像元进行地理配准,选取中心坐标处的先验地表反射率像元值为先验地表反射率(如图3所示)。该方法不仅可以保持识别结果的空间分辨率,还可以保留高空间分辨率影像中丰富的信息。以GF-2 PMS数据为例,当地表反射率数据集的空间分辨率为500 m时,一个地表反射率像元范围对应125×125个GF-2像元。当该范围内存在云时,这些GF-2像元中就会存在反射率较高的有云像元,基于先验地表反射率利用动态阈值模型就能够实现有云像元与晴空像元的区分。

图 3. GF-2 PMS与地表反射率数据集像元配准示意图

Fig. 3. Sketch map of pixel registration of GF-2 PMS and surface reflectance dataset

下载图片 查看所有图片

4 实验结果与分析

本课题组选择了15景GF-1 PMS影像、82景GF-1WFV影像、3景GF-2影像、37景GF-4影像和38景ZY3影像进行云识别试验。选取的实验数据尽可能包含不同类型的云(如厚云、薄云和碎云等),同时覆盖多种下垫面类型,如植被、土壤、沙漠和人工地表等。使用改进后的UDTCDA和UDTCDA分别进行云识别实验,并采用遥感目视解译方法对云识别结果进行对比和验证。

4.1 实验结果

选取改进后的UDTCDA对厚云、薄云和碎云这三种云类型的识别结果进行分析(如图4~6所示)。图4(a)、(b)分别为原始假彩色合成影像(RGB:近红外-红-绿)和云识别结果;图4(c)、和(d)分别为局部地区假彩色合成影像和云识别结果,其中白色表示云,黑色表示为背景。厚云在遥感影像中表现为大面积的亮斑,反射率整体较高,识别较为容易。从图4中可以看出:改进后的UDTCDA对厚云的识别精度较高,误判和漏判的情况较少(如图4(b)所示);厚云边缘的薄云由于反射率与地表反射率的差异小,往往会产生漏判。由局部识别结果可以看出,改进后的UDTCDA不仅能够以较高的精度实现厚云的检测,对厚云边缘的薄云也具有较好的检测效果,提高了对于厚云的整体识别准确度(如图4(d)所示)。

图 4. 改进后的UDTCDA对ZY3卫星数据中厚云的识别结果。(a) ZY3卫星原始假彩色合成影像;(b) ZY3云检测结果;(c) ZY3局部假彩色合成影像;(d) ZY3局部云检测结果

Fig. 4. Thick cloud detection results from the improved UDTCDA for ZY3 images. (a) The original false color composite image of ZY3; (b) the corre-sponding cloud detection result of ZY3; (c) the local standard false color composite image of ZY3; (d) the corresponding cloud detection result of local ZY3 image.

下载图片 查看所有图片

图5为改进后的UDTCDA对GF-1 WFV中薄云的识别结果。由于薄云的透过性高,其反射率容易受下垫面地物类型的影响,导致云与下垫面之间的反射率差异较小,难以通过固定阈值来区分有云像元与晴空像元。图5的识别结果显示,先验地表反射率支持的云识别方法可以有效降低下垫面的影响,能将大部分薄云识别出来(如图5(b)所示)。从图5(c)、(d)所示的局部区域中可以看出,虽然改进后的UDTCDA对薄云能够达到较好的识别效果,但随着云层厚度的降低,漏判的情况增加。

图 5. 改进后的UDTCDA对GF-1 WFV卫星数据中薄云的识别结果。(a) GF-1 WFV卫星原始假彩色合成影像;(b) GF-1 WFV云检测结果;(c) GF-1 WFV局部假彩色合成影像;(d) GF-1 WFV局部云检测结果

Fig. 5. Thin cloud detection results from the improved UDTCDA for GF-1 WFV images. (a) The original false color composite image of ZY3; (b) the corresponding cloud detection result of GF-1 WFV; (c) the local standard false color composite image of GF-1 WFV; (d)the corresponding cloud detection result of local GF-1 WFV image

下载图片 查看所有图片

图6为改进后的UDTCDA对GF-1 PMS影像中碎云的识别结果。碎云在卫星影像中表现为面积小、分布广泛的团块状云层,其空间尺度与空间分布给云识别造成了很大困难。如图6的识别结果所示,在对配准方式进行改进后,云识别结果与原始假彩色影像中的云层分布整体吻合度较高,对植被与裸地上空的云能达到较好的识别效果,漏判与误判的情况较少(如图6(b)所示)。从图6(c)、(d)中可以看出,改进后的UDTCDA不仅提高了碎云的识别效果,而且对面积小、云层薄的碎云也有较高的识别精度。

图 6. 改进后的UDTCDA对GF-1 PMS卫星数据中碎云的识别结果。(a) GF-1 PMS卫星原始假彩色合成影像;(b) GF-1 PMS云检测结果;(c) GF-1 PMS局部假彩色合成影像;(d) GF-1 PMS局部云检测结果

Fig. 6. Broken cloud detection results from the improved UDTCDA for GF-1 PMS images. (a) The original false color composite image of GF-1 PMS; (b) the corresponding cloud detection result of GF-1 PMS; (c) the local standard false color composite image of GF-1 PMS; (d) the corresponding cloud detection result of local GF-1 PMS image

下载图片 查看所有图片

4.2 精度验证

为了定量评价改进后的UDTCDA在不同高分辨率卫星上的适用性,利用遥感目视解译的方法从原始假彩色影像中提取云信息作为真实信息,通过构建误差矩阵的方式(表3),计算并验证该算法云识别结果的精度。同时选择总体精度A_O、生产者精度A_P、用户精度A_U、错分误差E_C和漏分误差E_O共5个评价指标来描述该算法的准确性和不确定性^[26]。表3中:N为像元总数; ρ11为改进后UDTCDA准确识别的云像元数量; ρ1+为改进后UDTCDA识别的所有云像元数量;ρ₁₂为改进后UDTCDA误提的云像元数量;ρ₊₁为目视解译结果中真实的云像元数量;ρ₂₂为改进后UDTCDA准确识别的晴空像元数量。

获取不同传感器影像中1000×1000大小的子图像对云识别结果进行展示,并使用上述方法对检测结果的精度进行验证。选取的影像包含下垫面为植被、城镇、裸土和沙漠等不同类型地表的厚云、薄云和碎云等多种云类型。图7展示了不同传感器不同下垫面对不同类型云层的识别结果,每一组图像从左到右依次是原始标准假彩色合成影像、改进后的UDTCDA和UDTCDA的云检测结果。从图7可以看出,UDTCDA能够识别不同传感器影像中的大部分云信息,但由于该算法利用重采样对影像与先验地表反射率数据进行配准,降低了原始影像的空间分辨率,导致原始影像的光谱信息丢失。从GF-4云识别结果可以看出(图7(m)~(p)),UDTCDA对中等空间分辨率影像中云的整体识别效果较好,但与改进后的UDTCDA相比仍存在少量漏判(图7 (p))与误判(图7(n))的情况。当将UDTCDA应用于GF-2影像云识别时,识别结果的空间连续性大幅降低,云层边缘锯齿化严重,云层细节信息大量丢失。尤其是对于云量较小的碎云,漏判情况比较严重(图7(k)、(l))。对比不同空间分辨率影像的云识别结果可以发现,影像的空间分辨率越高,UDTCDA识别结果的漏判情况越严重。然而,改进后的UDTCDA能够较好地解决这一问题,云识别结果的空间连续性和细节得到了显著提高,有效地降低了云识别的不确定性。

图 7. 典型区域的云识别结果

Fig. 7. Cloud detection results over different typical areas

下载图片 查看所有图片

图7(a)~(d)为ZY3影像在不同地表类型上的云识别结果,可以看出:改进后的UDTCDA对植被上空厚云、薄云和碎云的识别结果与原始影像中云层的分布具有较好的一致性(如图7(a)和(b)所示),而且使用先验地表反射率数据集可以有效降低高亮地表的影响,使高亮地表上空的厚云、碎云(如图7(c)所示)以及薄云(如图7(d)所示)的识别效果整体较好,但薄云仍存在漏判的情况。整体来看,改进后的UDTCDA对于ZY3多光谱传感器的适用性整体较高,平均识别精度达到了96.32%(表4),漏分和错分误差整体分别低于8%和9%。图7(e)~(h)为GF-1 WFV影像中植被、裸土和沙漠上空的厚云、薄云和碎云的识别结果。从识别结果可以看出:改进后的UDTCDA对于厚云和碎云的整体识别效果较好;而对于薄云,特别是高亮地表上空的薄云,漏判情况较多,平均生产者精度和用户精度分别为88.18%和91.38%。图7(i)、(j)为GF-1 PMS影像的云识别结果。从植被与城镇上空的厚云、薄云和碎云的识别结果可以看出,改进后的UDTCDA对厚云与碎云的识别精度较高(如图7(i)所示),并且能够以较高的精度识别植被上空的薄碎云,但仍存在少量漏判和错判的情况,如图7(j)所示,漏分误差为14.31%,错分误差为12.94%。图7 (k)、(l)为GF-2影像中植被与裸土上空厚云和薄云的识别结果。可以看出,改进后的UDTCDA对稀薄云层识别时,会出现明显的漏判和误判,如图7(k)所示,漏分误差和错分误差分别为22.40%和20.98%,总体精度为83.11%。图7(m)~(p)为GF-4影像中植被、城镇、裸土和沙漠上空不同类型云的识别结果。改进后的UDTCDA对植被上空云的识别效果最好,总体精度为95%(如图7(m)所示)。但对于部分薄云仍存在漏判和误判的情况。整体来看,改进后的UDTCDA对GF-4云识别的适用性较高,对不同类型下垫面上的云检测均能够达到较高的识别精度,平均总体精度为93.92%。

5 结论

为了降低混合像元对云识别的影响和重采样对原始影像造成的损失,本课题组在UDTCD的基础上采用空间配准的方式实现了非同源影像的配准,既实现了卫星影像与地表反射率数据集的准确配准,又保留了原始影像的空间尺度。然后,将改进后的UDTCDA与原始UDTCDA的云检测结果进行对比。

结果表明,改进后的UDTCDA在国产系列卫星云识别中能够取得较好的识别效果。与UDTCDA的云检测结果相比,改进后的UDTCDA算法大幅度提高了碎云、薄云和云边缘的识别准确度。精度验证结果表明,改进后的算法能够较好地适用于不同国产陆地观测卫星的云识别,对ZY3、GF-1 WFV、GF-1 PMS、GF-2和GF-4卫星的平均云检测精度分别为96.32%、94.04%、96.31%、85.98%和93.92%,整体云识别精度大于85%,漏判和错判的情况相对较少,两种误差整体分别低于9.57%和10.40%。改进后的云识别算法对不同下垫面上空的不同云类型均能够取得较好的识别效果,对植被和湿地等低反射率地表上空的云层具有整体较高的识别精度(A_O>88%);同时,对于高亮地表(裸土、城镇和沙漠等)上空的云识别平均精度约为83%,漏判情况整体较少。验证和对比结果表明,改进后的UDTCDA对高分辨率卫星数据具有较强的适应性,可以用于该系列数据的云检测业务化运行。但仍存在以下不足:1)受限于地表反射率数据集的精度,由于自然或人为原因造成的地表反射率变化明显的地区或具有崎岖复杂地形的地区,算法的适用性降低;2)数据间的配准误差,在河流、海岸线与陆地交界等地区,容易产生误判现象。