1 引言

在遥感影像获取过程中,快速的云检测能够及时提供准确的云覆盖比例以评价影像质量,进而星上预筛选含云量较高的劣质图像,可以减小下传的数据量,提高遥感数据的传输效率,节约传输和储存等成本,提高目标检测效果。因此遥感图像云检测已成为光学遥感图像处理过程中亟待解决的问题之一^[1-2]。目前最常用的云检测算法是光谱阈值法或利用遥感图像的纹理信息进行云区和下垫面的判别。

纹理特征不仅反映了图像的灰度统计信息,还反映了地物本身的结构特征和地物的空间排列关系,因此可以作为区分不同地物的依据^[3]。典型的纹理特征分析方法有灰度共生矩阵法、分形维数法、空间变化检测法等。如Solvsteen^[4]提出的空间变化检测法用于海洋背景的检测,但是仅适合于简单下垫面上的云检测,同时计算量较大;陶淑苹等^[5]采用小波尺度共生矩阵(SCM)提取纹理特征,并利用基于角二阶矩(ASM)和熵的双判别方式辨别云层,取得了一定的成果,但受到传感器类型不同、影像分辨率不同等因素的限制,用于资源三号卫星影像时仍存在较高的误判率;单娜等^[6]提出高准确度的云检测算法,算法选用的纹理特征为分形维数和角二阶矩,在实现上应用一种基于树状判别结构的快速算法,可以显著提高算法的运行效率。只利用纹理特征的云检测算法大多需要提取图像的特征或进行复杂的决策判断,其计算量比较大,并不适用于星载高光谱图像的快速云检测。

对于高光谱遥感图像,由于光谱波段较多,常常采用光谱阈值法作为云区和下垫面的决策判别依据。光谱阈值法综合考虑了几种光谱谱段进行阈值检测,最后再判断某个给定的像素是否是云^[7-8]。如冯书谊等^[9]针对地球观测卫星-1(EO-1)Hyperion传感器的遥感图像数据采用多谱段阈值的云检测方法,取得了较好的效果。光谱阈值法的缺点是识别精度不高,很难将薄云与下垫面区分开来。

而将光谱阈值法与纹理特征相结合的高效快速算法鲜有报道。同时在现有的云检测算法中,大多只是区分云区和下垫面,对于厚云和薄云的决策判决研究较少^[10]。郁文霞等^[11]在进行云自动检测中提出根据云、雪表面的不规则度计算分形维数,区分云和雪;刘湘航等^[12]通过人机交互的方式对已检测出的云、雪区域进行手动剪切子区,然后进行子图分形维数的计算实现云的识别。传统的分形维数计算方法需要对子图进行多次遍历,每次遍历包含大量复杂运算,导致算法运行效率低。同时现有的算法均采用固定的分形维数阈值,其缺点是不能自适应图像的变化,不能有效地应用到不同卫星的遥感图像上。

本文在光谱反射特性的基础上结合纹理特性,提出基于光谱分析和动态分形维数的云检测算法。光谱阈值因卫星传感器的不同而不同,以EO-1卫星的Hyperion传感器拍摄的高光谱遥感图像为例,以光谱阈值法为主,辅以动态分形维数法。所提算法的重点在于解决光谱分析云检测精度不高和利用纹理特征速度过慢的问题。所提算法同时增加了对薄云区的检测,实验结果表明,所提算法既可达到星载高光谱图像快速云检测的要求,又极大地提高了遥感图像云检测的精度,还能够有效地检测厚云区和薄云区。

2 基于多谱段辐射量特性的遥感图像云检测

文献[ 9]针对EO-1卫星Hyperion传感器的遥感图像数据提出了基于反射率特性的高光谱遥感图像云检测方法,该方法根据各类地物在不同光谱下的反射率特征,从星载高光谱图像中选取5个特定波长的图像数据,然后依照图像地物反射率及阈值判定法区分出高云区和中低云区,最后对二者进行合并,得到目标云区。参照其方法步骤,根据光谱曲线,利用地物的辐射量特性区分厚云区、薄云区和下垫面。在文献[ 9]的基础上,研究不同下垫面的含云遥感图像,增加了对厚云区和薄云区的检测与分析。为下文叙述方便,简称为多光谱辐射特性(MRC)算法。

图1(a)为由EO-1卫星的Hyperion传感器于2017年1月17日拍摄于挪威的遥感图像,并截取400 pixel×256 pixel大小的子图。图中含有厚云、薄云、海洋、陆地和雪,用遥感图像处理平台(ENVI)提取各个地物的光谱曲线,如图1(b)所示。图中各曲线由上到下依次为厚云、薄云、雪、陆地和海洋。

由光谱曲线分析可以看出,在光谱仪的光谱覆盖范围(355~2577 nm)内,厚云的辐射量最大,因此在大多数图像中均表现为亮色调区域,在1043.590 nm(第88波段)达到局部峰值;而薄云和海洋的光谱曲线重合度较高,不易区分,但是二者在蓝绿可见光波段的辐射量有一定的偏差,可以作为分类判别依据,如在波段457.34 nm(第11波段)处,辐射量2500 W/m²可以作为薄云区和下垫面的阈值。

通过对含云遥感图像光谱曲线辐射特性的分析,并根据特定谱段下地物辐射量的差异,采用光谱辐射量阈值法进行云检测。算法步骤如下:

1) 选择1043.590 nm,辐射量大于3000 W/m²可判断为厚云区。其他部分待定。

图 1. 含厚云、薄云、雪、海洋、陆地的遥感图像及其光谱曲线。(a)遥感图像;(b)光谱曲线

Fig. 1. Remote sensing image with thick clouds, thin clouds, snow, sea, and land and their spectral curves. (a) Remote sensing image; (b) spectral curves

下载图片 查看所有图片

2) 待定部分选择457.34 nm,辐射量大于2500 W/m²暂定为薄云区。

3) 暂定为薄云区部分选择548.92 nm(第20波段),辐射量大于2200 W/m²可判断为薄云区。其他部分即为下垫面。

3 基于动态分形维数和辐射量特性的云检测算法

考虑到薄云和下垫面光谱反射值较为接近,仅仅依靠光谱曲线检测精度不够高,细节部分以及薄云区可能会检测不出。而图像的纹理特征不仅反映了图像的灰度统计信息,还反映了地物本身的结构特征和地物空间排列的关系,可以作为区分不同地物的依据。

选择云的分形维数作为云的纹理特性。分形维数值用于表征图像表面的不规则度,其值越大,表明图像表面越复杂,对应的纹理越粗糙;反之则表明图像表面越光滑。遥感图像中地面目标纹理细节丰富,其分形维数值偏大;云影像灰度值变化平缓,纹理细节少,其分形维数值偏小^[13-14]。因此可以利用分形维数值作为辅助区分云区和地物的指标^[15]。

传统利用分形维数作为分类指标的方法中,一般都采用先验固定的分形维数作为分类依据^[16]。先验固定的分形维数是基于前人的研究结果,将云的分形维数作为云固有不变的属性,对所有的遥感图像均采用先验固定值。其缺点是适应性差,分类效果不好,难以发挥纹理特征的优势^[17-18]。而本研究提出的动态分形维数法则利用当前遥感图像的信息,得到该幅遥感图像的含云子图,计算含云子图的分形维数,从而将得到的云的分形维数的动态范围作为阈值的判断依据,相较先验固定的分形维数法,所提方法具有自适应性。为下文叙述方便,简称为动态分形维数与辐射特性(DFD_RC)算法。

3.1 分形维数的计算及分析

遥感图像是地面三维空间在二维表面的投影,遥感图像的结构信息量可以通过图像上灰度的空间变化来反映。采用由Sarkar^[19]提出的计盒维数法计算分形维数。把二维遥感图像想象成三维空间中的曲面,(x,y)表示二维空间的像素坐标,z为图像灰度值。将M×M大小的图像(x,y)空间按照尺寸r分割成许多r×r的小区域。在每一个小区域上是一列r×r×h的盒子,其中盒子的高度h可通过图像的总灰度级G来确定,即

h=rM·G。(1)

设第(p,q)块区域中图像灰度级最小值与最大值分别落在第K和L个盒子中,则覆盖第(p,q)个区域所需的盒子数为

nr(p,q)=L-K+1。(2)

统计覆盖整个图像所需的盒子总数N_r为

Nr=∑p,qnr(p,q),(3)

分形维数值为

D=limlnNrln(1/r)。(4)

3.2 分形维数动态范围的确定

为了得到云的分形维数的动态范围,需要首先确定含云子图。不同的子图尺寸会对云检测的准确性带来一定的影响,若子图的尺寸过大,则包含的信息较多,影像中可能同时含有云类和地物类,或者含有不同类型的云。由于图像的特征是全图的综合反应,此时的子图不能作为一个云类或地物类的单元,且子图尺寸越大时,全图的取样数就越少,会导致准确率受到影响^[20-21]。若子图的尺寸过小,则所含信息量小,不足以体现云和地物的纹理区别,甚至出现部分云和地物的子图表现几乎一致,人眼也难以区分的情况^[22-23]。通过大量实验选择尺寸为50 pixel×32 pixel的子图进行分割,分割效果最佳,即分割成8×8个子图。

在分割原图的基础上,利用二均值聚类的方法将含云子图提取出来。二均值聚类即按照下述步骤将图像分成目标(云区)和背景(非云区)两部分。具体步骤如下:

1) 将图像等分成8×8个子图。

2) 求出图像最小灰度值R_min和最大灰度值R_max,计算初始阈值T₀=(R_min+R_max)/2。

3) 根据阈值将图像分为云区和非云区两部分,并求出两部分均值,即

RO=∑R(i,j)>TkR(i,j)·N(i,j)∑R(i,j)>TkN(i,j),(5)RG=∑R(i,j)≤TkR(i,j)·N(i,j)∑R(i,j)≤TkN(i,j),(6)

式中:R_O为目标灰度均值;R_G为背景灰度均值;R(i,j)为图像上(i,j)点的灰度值;N(i,j)为(i,j)点的权重系数;T_k为第k次迭代的阈值。

4) 重新选择阈值T_k+1=(R_O+R_G)/2。

5) 对整幅图像和每个子图分别循环步骤3)和4),直到T_k=T_k+1。

按照上面的步骤,求出整幅图像的阈值T_all和子图的阈值T_i(i=1,2,…,64),并作初步判断,

Ti>λ×Tall,cloud⁃containingsubgraphTi≤λ×Tall,non⁃cloud⁃containingsubgraph,(7)

式中:系数λ控制的是含云子图的个数和云的分形维数确信度之间的关系。λ越大,提取的含云子图数量越少,但是得到当前云图分形维数的动态范围的确信度越高。由大量实验可知,λ=1.5时可以较好地衡量二者之间的关系。

按照3.1节中分形维数的计算方法求出上述含云子图的分形维数,从而确定每幅图片云区分形维数的动态范围(d_min,d_max)。

3.3 云检测算法流程

基于动态分形维数和辐射量特性相结合的云检测综合算法的流程图如图2所示。由图1(b)可以看出:在光谱1043.590 nm处,厚云的辐射量最大,达到局部峰值,在大多数图像中均表现为亮色调区域;而薄云和海洋的光谱曲线重合度较高,不易区分。可以选择合适的阈值作为区分厚云的依据(如图2所示,可以选择阈值y₁=800 W/m²,y₂=3000 W/m²)。薄云和下垫面在蓝绿可见光波段的辐射量有一定的偏差,可以作为分类判别依据^[24]。如在波段457.34 nm处,辐射量阈值y₃=2100 W/m²可以作为薄云区和下垫面的阈值。

对于不同来源的遥感图像数据需要人工参与调节3个参数,但是仅需要微调即可。因为通过分析不同图像的光谱曲线可以看出:地物的光谱辐射特性是地物的固有属性,不因遥感数据采集来源的不同而改变^[25]。同时选择的波段均是能明显区分不同地物特征的波段,阈值的选择范围较广,对不同的遥感数据仅需要微调光谱阈值即可。

图 2. 基于动态分形维数和辐射量特性的云检测算法流程图

Fig. 2. Flow chart of cloud detection algorithm based on dynamic fractal dimension and radiation characteristics

下载图片 查看所有图片

4 实验和结果分析

4.1 云检测的样本选择

由于遥感图像尺寸过大,缩略显示无法看清云层纹理细节,为了观测检测结果,截取图中典型云域实际尺寸显示其检测过程及结果。选择EO-1卫星的Hyperion传感器拍摄的遥感图像,如图3所示。每幅遥感图像的大小均为400 pixel×256 pixel,其详细信息如表1所示。

图 3. 遥感云检测样本图像。(a)样本A;(b)样本B;(c)样本C;(d)样本D;(e)样本E;(f)样本F;(g)样本G;(h)样本H

Fig. 3. Sample images of remote sensing cloud detection. (a) Sample A; (b) sample B; (c) sample C; (d) sample D; (e) sample E; (f) sample F; (g) sample G; (h) sample H

下载图片 查看所有图片

4.2 算法比较及分析

为了定量评价所提云检测方法,参考文献[ 26]采用人工目视解译方式勾取的云区域作为真值。勾绘准则为:云边界模糊区域勾绘精度小于3 pixel,云边界清晰区域勾绘精度小于1 pixel^[26],如图4所示。根据参考文献[ 26]的方法得到所有样本真值云图的云含量占比,如表2所示。其中,云含量占比的统计方法为云的像素个数与整幅图像像素个数的比值^[26],云含量为厚云含量和薄云含量的总和。

图5和图6分别表示MRC算法和DFD_RC算法对图3样本图像灰度图的检测结果,其中白色表示厚云(灰度值为255),灰色表示薄云(灰度值为150),黑色表示下垫面(包括陆地和海洋,灰度值为0)。两种算法云含量的检测结果如表3所示,分别统计了两种算法所检测出的云含量占比、厚云含量占比以及薄云含量占比。分类指标如表4所示,分别统计了两种算法对于云含量检测的召回率、漏警率、精确率与虚警率。

图 4. 真值云图。(a)样本A;(b)样本B;(c)样本C;(d)样本D;(e)样本E;(f)样本F;(g)样本G;(h)样本H

Fig. 4. True value cloud maps. (a) Sample A; (b) sample B; (c) sample C; (d) sample D; (e) sample E; (f) sample F; (g) sample G; (h) sample H

下载图片 查看所有图片

图 5. MRC算法检测结果。(a)样本A;(b)样本B;(c)样本C;(d)样本D;(e)样本E;(f)样本F;(g)样本G;(h)样本H

Fig. 5. Detection results of MRC algorithm. (a) Sample A; (b) sample B; (c) sample C; (d) sample D; (e) sample E; (f) sample F; (g) sample G; (h) sample H

下载图片 查看所有图片

图 6. DFD_RC算法检测结果。(a)样本A;(b)样本B;(c)样本C;(d)样本D;(e)样本E;(f)样本F;(g)样本G;(h)样本H

Fig. 6. Detection results of DFD_RC algorithm. (a) Sample A; (b) sample B; (c) sample C; (d) sample D; (e) sample E; (f) sample F; (g) sample G; (h) sample H

下载图片 查看所有图片

从实验结果可以看出,当下垫面是海洋且有少量厚云时(即样本A、B),两种算法检测结果相差不大。但是当厚云含量增加时(样本C),MRC算法将不能识别出大面积的厚云,漏警率大大增加,而DFD_RC算法的检测结果则大大降低漏警率。原因是MRC算法只依赖于光谱阈值,而阈值的选择因图像的不同会有差异,因此MRC算法的稳健性低于DFD_RC算法。

当下垫面含有陆地或雪,且有薄云存在时(样本D、E、F、G、H),MRC算法将薄云区域全部判定为陆地,漏警率很大,而DFD_RC算法能较为准确地检测出薄云。原因在于MRC算法仅利用地物的光谱信息,而薄云和陆地的光谱特性较为接近,很难用单一的光谱阈值将薄云和陆地区分开;而DFD_RC算法在利用光谱阈值的同时又考虑云的纹理特征,极大地提高了遥感图像云检测的精度。

同时,所提算法可以快速进行云检测。传统的分形维数计算方法需要对子图进行多次遍历,每次遍历包含大量复杂运算,导致算法运行效率低。而所提算法则是根据图2的算法流程图,首先根据光谱特性确定厚云区和下垫面;只是在薄云区和下垫面部分利用分形维数作为辅助分析;极大地减少了遍历次数,从而减少了计算量,提高算法运行效率。

为进一步分析两种算法的差异,图7和图8展示了结果差异较大部分的局部放大影像。图像特征场景分别为:下垫面为海洋并且同时含有厚云和薄云(样本B);下垫面为陆地并且同时含有厚云和薄云(样本H);以薄云为主的特征场景(样本G);以厚云为主的特征场景(样本C)。从图7(a)、(b)与图8(a)、(b)的比较结果中可以看出,所提云检测算法得到的云边界清晰,能检测出厚云和薄云得到的而MRC算法在某些厚云薄云交错区域的漏警率很大,原因是MRC算法只依赖于光谱阈值,而阈值的选择因图像的不同会有些许差异。图7(c)与图8(c)是含有大量薄云的情况,通过比较可以看出,MRC算法不能识别出大量的薄云,而所提算法能较为准确地识别出薄云。但是所提算法仍有一些不足,如图8(c)方框所示,对于透明度较高的薄云,可以看见其下垫面,但往往会归于地物。图7(d)与图8(d)是含有大量厚云(层云、卷云)的情况,通过比较可以看出,MRC算法将在阴影区域的大量卷云识别成薄云,而且识别精度不高,存在大量漏检现象,而所提算法则极大地改善了这一情况,能较为准确地识别厚卷云。

图 7. MRC算法的部分放大图像。(a)样本B;(b)样本H;(c)样本G;(d)样本C

Fig. 7. Partially enlarged images of MRC algorithm. (a) Sample B; (b) sample H; (c) sample G; (d) sample C

下载图片 查看所有图片

图 8. DFD_RC算法的部分放大图像。(a)样本B;(b)样本H;(c)样本G;(d)样本C

Fig. 8. Partially enlarged images of DFD_RC algorithm. (a) Sample B; (b) sample H; (c) sample G; (d) sample C

下载图片 查看所有图片

5 结论

以EO-1卫星Hyperion传感器的高光谱遥感图像为例,介绍了光谱阈值判别法和纹理特征分析法的基本原理,分析比较了各自的优缺点。在光谱反射特性的基础上引入动态分形维数,提出了将光谱阈值和纹理特征相结合的云检测算法。算法以光谱阈值法为主,辅以动态分形维数法,既达到星载高光谱图像快速云检测的要求,又极大地提高了遥感图像云检测的精度,还能够有效地检测厚云区和薄云区。但是所提算法也存在一定的不足,对透明度较高薄云的正确识别有待提高。