CN101799921A - 一种光学遥感图像云检测方法 - Google Patents

Abstract

本发明提供一种光学遥感图像云检测方法，包括下列步骤：1)将所述遥感图像划分为基本处理单元；2)根据灰度阈值将基本处理单元划分为可能为云和不可能为云两类；3)对于可能为云的基本处理单元，根据分形维数划分；4)对于分形维数划分结果中的待定单元根据角二阶矩判定为云层或地面目标物体。本发明将三个判别特征参数以合理的次序结合起来，不同的特征参数在检测过程中发挥不同的作用；从而不仅能够获得较高的分辨精确度，还能够在平均意义上显著提高方法运行效率。

Description

一种光学遥感图像云检测方法

技术领域

本发明涉及光学遥感图像处理技术领域，尤其涉及一种对遥感平台传感器获取的光学图像进行云片检测和云层相关参数提取的方法。

背景技术

随着科技的高速发展，卫星遥感技术进入了一个能快速及时地获取对地观测海量信息的新阶段。这些丰富的数据信息被广泛的应用在农业、气象、环境保护、资源调查等各个领域，有着重要的应用价值。在实际情况下，由于天气情况复杂多变，遥感卫星获取卫星图像，尤其是获取大范围的遥感影像时，多数情况下都会有云的存在，真正的无云晴空地面目标物体影像很难获取。遥感图像上，云的存在会遮蔽或者覆盖云层区域下面的真实的地面目标物体的信息。为了能够有效的从遥感图像数据中提取地面目标物体的信息，提高遥感卫星图像数据的可用性和利用率，必须对遥感图像中存在的云进行检测，采取必要的措施，以消除或降低云的干扰和影响。由于云的种类繁多，纹理呈现随机分布特性，其形态、特性等随时间和空间发生变化，同时还受季节和气候等多种条件的综合影响，即不同季节、不同空间高度、不同天气条件下的云的反射率和亮度都有所不同，这对用同一方法应用于不同条件下获取的遥感图像数据的云检测结果精度带来了考验。

另一方面，遥感平台搭载的传感器的类型不同，其成像原理也是不同的，获取的遥感图像数据也各有不同的性质和特点。针对不同类型的遥感图像数据，产生了很多种云检测方法。

对于单幅光学遥感图像，其所能够提供的信息量有限，不像多光谱遥感数据中有丰富的光谱信息可以使用。在单光谱遥感图像中对图像中的信息进行分类和提取要困难得多。这就对检测方法提出了很高的要求。而且，由于传感器获取的遥感数据基本为海量数据，遥感信息处理平台对处理效率的要求很高。否则，即使获得十分精确的判别结果，也没有实际的应用价值。因此，判别精确度和方法的执行效率都是评价方法可用性的重要标准。

现有的应用较广泛的方法按照是否需要辅助信息主要分为两类：只需图像本身数据信息的方法和利用历史数据信息的方法。利用历史数据信息的方法，通常是将之前在同一个地理位置上空获取的无云的遥感图像保存在遥感平台上，通过位置等信息标定后，比较对应像素点信息，进行判定。这种方法的优点是，不仅能够检测到是否有云层存在，并且可以直接用无云遥感图像上的信息点替换云层遮盖的区域，得到无云覆盖的图像数据。但是对于遥感平台有限的资源来说，存储海量的历史数据对平台造成极大的负担。对于没有历史数据的区域无法应用此方法。无需辅助信息的方法，又主要分为光谱分析和纹理分析两类。光谱分析利用云与地面物体在可见光和红外波段的反射率和温度的差异进行分辨。纹理分析利用物体的局部或整体有着某种程度的自相似性，通过提取目标物体的空间统计特性，将各类目标物体归类。

光谱空间信息主要是指图像的灰度值信息。云在可见光/近红外波段对于光线的反射率比大多数地面物体要强的多，在卫星遥感图像上表现为相对于地面目标有较高的灰度值。利用云的这个特性，可以使用阈值判断方法对卫星遥感图像进行粗略的云与非云的判别。这是较早的云层检测使用的方法，也是实现简单、效率较高的检测方法。图像光谱阈值的确定大体有两种方法：一是根据经验，通过对数据库中大量的历史样图进行统计，获取云层灰度值的经验阈值；二是利用图像本身的信息，确定阈值的选取。经验阈值的方法可以在统计意义上的接近最优，但是却没有考虑图像本身的信息特点，所以不能保证对每幅遥感图像都能得到很好的分辨结果。利用图像自身信息的方法，虽然充分利用了图像自身数据特点，却可能会遇到无法正确选取阈值的问题。光谱分析方法的优点是方法实现简单，能够获得很好的执行效率，在某些特定的场景下有非常好的判别结果，例如，黑暗背景下且无噪声干扰的遥感图像。但是其缺点是当图像中云层与地面背景的光谱值相差不明显时，判定结果精确度将大大降低。尤其是雪和沙漠等对光线反射率大的物体，很难进行区分；而且对阈值比较敏感，阈值的选择有一定的主观性，精确度容易受噪声的影响；同时，光谱分析以像素为处理基本单位，没有考虑图像的结构信息。

虽然云的纹理基元多变难测，属于随机纹理，但它仍有别于下垫面物体的纹理特征。基于云层纹理特征的分析利用云的局部和整体有着某种程度的相似性，从图像的统计特征出发，提取云与地球表面物体的纹理特征，对云和地球表面物体进行有效的区分。在云检测方法中最常用的纹理特征有分形维数和角二阶矩值。纹理分析方法虽然能得到比阈值法更高的判别精确度，但是这种方法要求足够的样本数据进行训练以得到不同类别物体的特征值区间，并且样本质量和数量直接影响训练区间的准确性，继而影响方法的判别精确度。而且纹理分析的方法忽视了云的独特光谱特征。因此，仍然存在一定数量的误判和漏判情况。

较高的判定精确度往往意味着较高的检测复杂度，在实际的应用中通常需要在精度和效率之间进行取舍，需要根据遥感图像数据的特点和应用场景等综合决定选择使用哪种检测方法。为了充分利用各种信息，提高云检测结果精确度，迫切需要一种将多种判别标准综合起来，充分利用遥感图像的光谱信息和纹理信息，从多个角度提取云层特征参数，对图像进行更加准确有效的云检测的方法。

发明内容

本发明要解决的技术问题是提供一种光学遥感图像云检测方法，以充分利用遥感图像的数据信息，从不同的判别空间提取目标的特征参数，从而对图像目标进行更准确有效的检测。

为实现上述目的，根据本发明的一个方面，提供了一种光学遥感图像云检测方法，包括下列步骤：

1)将遥感图像划分为基本处理单元；

2)根据灰度阈值将基本处理单元划分为可能为云和不可能为云两类；

3)对于可能为云的基本处理单元，根据分形维数划分；

4)对于分形维数划分结果中的待定单元根据角二阶矩判定为云层或地面目标物体。

其中，步骤3)包括：31)对于可能为云的基本处理单元中极有可能为云的基本处理单元，根据分形维数划分为云层单元和云层待定单元，该云层单元的判定结果为云层；

步骤4)包括：41)根据角二阶矩判定该云层待定单元为云层或地面目标物体。

其中，步骤3)包括：32)对于该可能为云的基本处理单元中既可能为云也可能为下垫面的基本处理单元，根据分形维数划分为下垫面单元和下垫面待定单元，该下垫面单元的判定结果为地面目标物体；

步骤4)包括：42)根据角二阶矩判定该下垫面待定单元为云层或地面目标物体。

其中，步骤2)进一步包括：

20)统计基本处理单元中大于灰度阈值的像素个数；

21)计算大于灰度阈值的像素个数所占百分比；

22)根据该百分比将基本处理单元划分为极有可能为云、不可能为云和既可能为云也可能为下垫面三类。

其中，灰度阈值根据经验确定或根据对历史图像样本进行统计获得。

其中，步骤3)可包括：

301)根据盒维数计算方法计算分形维数。

其中，步骤3)可包括：

302)根据基于二维灰度二值盒维数计算方法计算分形维数。

其中，基本处理单元的大小为64×64。

本发明提出的树形结构将三个判别特征参数以合理的次序结合起来。不同的特征参数在检测过程中发挥不同的作用，恰当的次序和优先级别可以保证各判别特征能够充分发挥优势同时互相弥补缺点，从而不仅能够有效的区分云区和非云下垫面，获得较高的分辨精确度，还能够在平均意义上显著提高方法运行效率。具体而言，本发明能够充分的利用遥感图像中光谱空间和位置空间的信息，实现云层的自动检测，且判别精确度高。对子图图像的区分可能不必计算所有的特征参数。在判别树中，叶子节点即为判别的最终结果，因此，只要沿着某一条路径走到叶子节点，即可得到判别结果。分形维数值和角二阶矩值的运算复杂度比较高，避免这两种运算将会使方法的效率明显提高。该树形结构实现复杂度低，基本满足实时性，尤其适用于遥感图像数据量非常大的应用场景。

附图说明

图1是本发明的云层检测方法的框图。

图2是根据本发明一个具体实施例的计算分形维数的三维空间示意图。

图3是根据本发明一个具体实施例的计算分形维数的点聚图。

图4是根据本发明一个具体实施例的计算灰度共生矩阵的示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对根据本发明一个实施例的光学遥感图像云检测方法进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明提出的云层检测方法需要进行光谱分析和纹理分析。传统光谱分析以像素点为最基本的处理单元，而纹理分析以图像中能体现纹理特征的最小图元为最基本的处理单元。应用本发明提供的方法对光学遥感图像进行处理时，以诸如64×64大小的子图为基本处理单元和结果判定单元，以下统一以基本处理单元描述。

根据本发明的一个具体实施例，将遥感图像划分成64×64大小的图像基本处理单元，对图像边界处不足64×64大小的单元不做检测判定处理。由于遥感卫星图像一般数据量巨大，以64×64大小的图像基本处理单元为单位能够完成对图像的充分细分，边界处小于64×64大小的单元可以忽略，对判定结果不造成影响。对于每个基本处理单元，按照树形判别结构所定义的判别准则完成对基本处理单元的判定工作。下面结合图1详细描述对每个基本处理单元进行处理的具体步骤。

对于每个待判定的图像子图，即基本处理单元，根据物体的反射特性：云对光线的反射比下垫面物体对光线的反射要强的多，在遥感卫星图像上表现为云层的灰度值明显高于下垫面物体的灰度值，将基本处理单元中的每个像素的灰度值与设定的灰度阈值进行比较，统计基本处理单元中大于灰度阈值的像素点的个数。其中，灰度阈值可以选择经验阈值，例如，可见光、近红外波段的遥感图像其云检测使用的经验阈值通常为200；也可以通过对遥感图像典型云层样本数据进行统计，选择适合的灰度阈值，确定方法通常有统计法，灰度直方图确定法等。本领域普通技术人员可以根据对算法精确度、复杂度和效率的要求选择合适的方法来确定该灰度阈值。

根据公式(1)计算基本处理单元中大于灰度阈值的像素点的百分比。

根据计算得到的基本处理单元中大于灰度阈值的像素点的百分比数值，对基本处理单元进行界定，将其划分为可能为云和不可能为云两类，其中可能为云的基本处理单元又具体分为极有可能为云以及既可能为云也可能为下垫面两类。边界值的确定可通过经验选取，例如，选择85％和15％作为边界值，如果cloud_percent＞85％则判定为极有可能为云，如果cloud_percent＜15％则判定为不可能为云，如果15％≤cloud_percent≤85％则判定为既可能为云也可能不为云。当然，也可以通过标定历史图像样本，然后统计cloud_percent的范围来获得边界值，例如，标定并选择典型的三类图像各20个样本，统计其cloud_percent的取值范围，根据统计值确定边界值。

本领域普通技术人员可以理解，除采用上述灰度阈值百分比进行判定的方法，还可以采用其他光谱判别方法，例如使用灰度平均值，灰度方差或者灰度熵等方法进行判别。

根据上述实施例，若根据大于灰度阈值的像素点的百分比将基本处理单元判定为不可能为云，则直接得到判定结果，该基本处理单元的最终判定结果为地面目标物体。

若根据大于灰度阈值的像素点的百分比将基本处理单元判定为极有可能为云，需要进行分形维数的判定，进一步从纹理分析的角度对其进行确认。如果该子图在纹理特征上也满足云层的特点，那么可以判定其为云层单元；如果不能满足云层纹理特征，但是从光谱阈值的角度它是云的可能性非常大，则为待定单元，称为云层待定单元，对该单元再进行角二阶矩值的判断，通过角二阶矩值减少漏判(是云但被判定为下垫面)的可能性。具体的，极有可能为云的基本处理单元的检测步骤为：

计算极有可能为云的基本处理单元的分形维数值D。

根据本发明一个具体实施例，采用盒维数计算方法计算基本处理单元的分形维数值。如图2所示，把二维灰度图像视作三维空间中的曲面(x，y，f(x，y))，其中f(x，y)为图像(x，y)位置处的灰度值。则计算分形维数D的具体执行步骤如下。

步骤1初始化尺度因子r，在本实施例中将尺度因子初始化为r＝4，本领域普通技术人员可以理解r只要在4≤r＜64范围内即可。然后执行步骤3。

步骤2如果r≥64(基本处理单元的大小)，则执行步骤4，否则r＝r×2，执行步骤3。

步骤3将基本处理单元分成r×r大小的子块，即64×64大小的图像被分割成许多r×r大小的网格。在每个r×r大小的网格上，假设有一列r×r×h的盒子，其中h表示一个盒子的高度。按照公式

可计算h。其中，G表示图像的总灰度级数，r是分割的尺度因子，

表示向下取整操作。由盒子的高度h根据公式256/h可以确定每个网格的上方的盒子的个数，进而可以确定每个灰度值落在第几个盒子里。

对图像中的每个网格做如下操作：

对于网格(m，n)找到网格中图像灰度最小值和最大值分别落在哪个盒子里，并记录在nr_min[m，n]和nr_max[m，n]中。根据公式(2)计算覆盖该网格中的图像灰度值所需的盒子数nr。

n_r＝nr_max[m，n]-nr_min[m，n]+1            (2)

根据公式(3)计算覆盖整个图像所需的盒子数。

$N_r=\underset{m,n}{\mathrm{\Σ}}{n}_r(m,n)---\left(3\right)$

执行步骤2。

步骤4以log(1/r)为横坐标，log(N_r)为纵坐标绘制点聚图，具体如附图3所示。利用最小二乘法拟合样本点，所得的斜率即为此图像样本的分形维数D。本领域的技术人员可以理解，除采用最小二乘法进行直线拟合以外，还可以通过使残差的最大绝对值为最小或使残差的绝对值之和为最小准则对直线进行拟合。

本领域的技术人员可以理解分形维数的计算方法，除了上述盒维数计算方法，还可以采用基于二维灰度二值盒维数计算方法等。

判断基本处理单元的分形维数值D在云层所在的区间范围内还是在地面目标物体的区间范围内。其中，区间范围通过对于典型的遥感图像云层样本和地面目标物体样本进行训练来获得，样本数量越多，质量越高，训练得到的边界值越精确。训练可以在边界值稳定时停止。

如果分形维数值在云层所在区间范围内，则其为云层单元，判定结果为云层。

如果分形维数值在地面目标物体所在区间范围内，则其为云层待定单元，调用角二阶矩值计算模块，根据如下具体步骤计算云层待定单元的角二阶矩值ASM。

首先将图像灰度级别由0-255压缩到0-31。

然后，如图4所示，统计图像中一个灰度为i而另一个灰度为j的两个相距为(Δx，Δy)的像素对出现的次数。得到灰度共生矩阵中位于(i，j)位置的元素p(i，j)的值。对于64×64大小的基本处理单元，通常Δx＝Δy，优选的取值范围是[2，6]。

P(i，j)＝#{(x1，y1)，(x2，y2)∈M×M|f(x1，y1)＝i，f(x2，y2)＝j}

其中#(A)表示集合A中的元素个数，显然P为32×32的矩阵。(x1，y1)，(x2，y2)为二维图像中像素点的位置，f(x1，y1)，f(x2，y2)分别为二维数字图像位置(x1，y1)，(x2，y2)处的像素值，M＝64是基本处理单元的大小。

最后，对灰度共生矩阵进行归一化处理，并根据公式(4)计算角二阶矩值ASM。

$\mathrm{ASM}=\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}{\left[P(i,j)\right]}^2---\left(4\right)$

判断云层待定单元的角二阶矩值在云层所在的区间范围内还是在地面目标物体的区间范围内。其中，区间范围通过对于典型的遥感图像云层样本和地面目标物体样本数据进行训练获得，样本数量越多，质量越高，训练得到的边界值越精确。训练可以在边界值稳定时停止。

如果角二阶矩值在云层所在区间范围内，则得到判定结果，该云层待定单元的判定结果为云层。如果角二阶矩值在地面目标物体所在区间范围内，则得到判定结果，该云层待定单元的判定结果为地面目标物体。

与极有可能为云的基本处理单元类似，对于判定为可能为云也可能为下垫面的基本处理单元，则如上所述计算基本处理单元的分形维数值D。判断基本处理单元的分形维数值在云层所在的区间范围内还是在地面目标物体的区间范围内。如果分形维数值在地面目标物体所在区间范围内，则该基本处理单元为下垫面单元，判定结果为地面目标物体。如果分形维数值在云层所在区间范围内，则该基本处理单元为待定单元，称为下垫面待定单元。如上所述计算下垫面待定单元的角二阶矩值ASM。判断下垫面待定单元的角二阶矩值在云层所在的区间范围内还是在地面目标物体的区间范围内。如果角二阶矩值在云层所在区间范围内，则得到判定结果，该下垫面待定单元的判定结果为云层。如果角二阶矩值在地面目标物体所在区间范围内，则得到判定结果，该下垫面待定单元的判定结果为地面目标物体。

根据上述判定结果分别统计被判定为云层和地面目标物体的基本处理单元数量，计算图像中云层基本单元所占总子图数的百分比，即得到图像的云层覆盖率。

本领域普通技术人员可以理解，利用本发明不仅能够获得遥感图像云层覆盖百分比，并且只需在方法的相应步骤进行少量修改即可获取云层覆盖在图像中的位置信息和边界信息等。

本发明提出的上述方法对于单光谱高分辨卫星遥感图像的云检测具有广泛的良好的通用性。实际测试结果表明，该云检测方法完全达到自动云检测的虚警概率小于5％，漏警概率小于10％的要求，并且能够获得很好的处理速度。

应该注意到并理解，在不脱离后附的权利要求所要求的本发明的精神和范围的情况下，能够对上述详细描述的本发明做出各种修改和改进。因此，要求保护的技术方案的范围不受所给出的任何特定示范教导的限制。

Claims (9)

1.一种光学遥感图像云检测方法，包括下列步骤：

1)将所述遥感图像划分为基本处理单元；

2)根据灰度阈值将所述基本处理单元划分为可能为云和不可能为云两类；

3)对于所述可能为云的基本处理单元，根据分形维数划分；

4)对于所述分形维数划分结果中的待定单元根据角二阶矩判定为云层或地面目标物体。

2.根据权利要求1所述的方法，其特征在于，

所述步骤3)包括：31)对于所述可能为云的基本处理单元中极有可能为云的基本处理单元，根据所述分形维数划分为云层单元和云层待定单元，所述云层单元的判定结果为云层；

所述步骤4)包括：41)根据所述角二阶矩判定所述云层待定单元为云层或地面目标物体。

3.根据权利要求1或2所述的方法，其特征在于，

所述步骤3)包括：32)对于所述可能为云的基本处理单元中既可能为云也可能为下垫面的基本处理单元，根据所述分形维数划分为下垫面单元和下垫面待定单元，所述下垫面单元的判定结果为地面目标物体；

所述步骤4)包括：42)根据所述角二阶矩判定所述下垫面待定单元为云层或地面目标物体。

4.根据权利要求1所述的方法，其特征在于，所述步骤2)包括：

20)统计所述基本处理单元中大于灰度阈值的像素个数；

21)计算所述大于灰度阈值的像素个数所占百分比；

22)根据所述百分比将所述基本处理单元划分为极有可能为云、不可能为云和既可能为云也可能为下垫面三类。

5.根据权利要求1或4所述的方法，其特征在于，所述灰度阈值根据经验值确定。

6.根据权利要求1或4所述的方法，其特征在于，所述灰度阈值是根据对历史图像样本进行统计获得。

7.根据权利要求1所述的方法，其特征在于，所述步骤3)包括：

301)根据盒维数计算方法计算所述分形维数。

8.根据权利要求1所述的方法，其特征在于，所述步骤3)包括：

302)根据基于二维灰度二值盒维数计算方法计算所述分形维数。

9.根据权利要求1所述的方法，其特征在于，所述基本处理单元的大小为64×64。

Images