CN102750701A - 针对Landsat TM和ETM图像的厚云及其阴影检测方法 - Google Patents

Abstract

本发明公开了一种针对Landsat TM和ETM图像的厚云及其阴影检测方法，其实现步骤为：（1）将输入图像分成16个子图集，并对其进行维纳滤波去噪和归一化；（2）对各子图集中的厚云和阴影进行粗检测，并从粗检测结果中选取基准对；（3）根据各基准对求最终基准对质心连线倾角和间距；（4）根据最终基准对质心连线倾角和间距对各子图集中厚云和阴影配对，并对未配对阴影和厚云进行补充检测；（5）对云影配对结果和补充检测结果进行相加，得到各子图集的最终检测结果子图；（6）顺序拼接各子图集的最终检测结果子图，获得最终检测结果图。本发明不需要辅助信息和人工参与，检测精确度高，可用于遥感图像变化检测、分类以及图像分割的预处理。

Description

针对Landsat TM和ETM图像的厚云及其阴影检测方法

技术领域

本发明属于光学遥感图像处理技术领域，涉及对专题制图仪TM图像和增强专题制图仪ETM图像的厚云及其阴影检测，可用于遥感图像变化检测、分类以及图像分割等工作的预处理。 

背景技术

随着空间技术的不断发展，利用中高分辨率的遥感图像进行陆地资源调查、变化检测、植被以及水域的监测等已经越来越普及。由于天气情况复杂多样，卫星在获取遥感图像时极少能遇到大范围区域内完全晴好无云的情况。云的存在使得人们无法从图像中获取云覆盖区之下的真实地表信息，这会对一些后续的分割、分类、变化检测等工作造成困难。为了消除由于云存在而造成的不利影响，提高遥感图像的可用性，人们提出了一些云检测和云去除的方法。 

云的形状、高度、厚度、成分及分布等情况各有不同，导致其对光学遥感图像产生的影响也具有较大差异。根据云是否完全遮挡地面信息，可将其分为薄云和厚云。薄云是指对太阳入射光和地表反射光只具有部分遮挡效果的云，在其覆盖区仍能反映出一定的地表信息，并且不产生明显的阴影。相对的，厚云是指对太阳入射光和地表反射光具有完全遮挡效果的云，在其覆盖区地表信息被完全屏蔽，并形成明显阴影，使阴影区的地面细节难以分辨。因此，对于薄云只需尽可能恢复或增强云覆盖区下的地面细节信息，即可达到消除其影响的目的，一般采用同态滤波法。而由于厚云遮蔽了在其覆盖下全部的地表信息，仅利用原始图像无法恢复，因而只能先检测出其位置、大小，再利用同一地区其他时期拍摄的无云图像进行像元替换，以达到消除厚云的目的。由于本发明只涉及厚云及其阴影的位置、大小的检测，因此对后续的像元替换不进行阐述。 

厚云的检测主要有灰度阈值法、空间纹理分析以及光谱分析三类方法。灰度阈值法是根据厚云灰度值比一般地物大这一特点用人工或自适应阈值将其检出，但冰雪也具有较大的灰度值，因此容易与云混淆而产生误检，同时检测效果很大程度受到图像本身和阈值选取的影响。空间纹理分析根据厚云与地物间的空间纹理差异进行区分， 最常采用的纹理特征有分形维数和角二阶矩值。虽然空间纹理分析检测效果较灰度阈值法好，但该方法通常需要训练样本，不但耗时且样本的数量和质量直接影响结果准确性。对于多光谱和高光谱数据，由于其包含大量光谱信息，因此广泛采用光谱分析法。普遍认为，厚云在可见光和近红外波段相对于其他地物具有较高反射率，而在热红外波段具有较低的亮温值。具体做法通常是在多个波段图像上通过分别设定阈值进行判断再集体决策，进而将厚云检测出来。中分辨率成像光谱仪MODIS图像具有36个波段，光谱分辨率较高，针对该图像的光谱分析法已经有大量比较成熟且简单高效的衍生算法。 

空间分辨率中高的专题制图仪TM图像，其光谱分辨率不高，只有7个波段，而增强专题制图仪ETM图像也只有8个波段。现有光谱分析法检测这两种图像中厚云的效果还不够理想，易受到与厚云具有类似光谱反射特性的冰雪干扰，因此常结合其他方法对检测结果进行改善。另外，由于专题制图仪或增强专题制图仪图像的第六波段比其他波段图像的空间分辨率低，光谱分析法使用起来非常不便。 

图像中有厚云存在时往往还伴有阴影，它会导致地表信息几乎无法分辨，因此在检测厚云的同时常常还将其阴影一同检测出来。2010年李炳燮、马张宝、齐清文等在文献“Landsat TM遥感影像中厚云和阴影去除”（遥感学报,2010,Vol.14,No.3,Pages:534-545）中通过对比分析同一地区不同时期图像中有云地域与无云地域间的光谱特征，提出了云地域和云阴影地域增强模型来检测厚云及阴影。该方法需要经过人工灰度匹配的一幅年份相近、相同季节的图像作为辅助，对辅助图像要求较高，并且需要人工参与。1996年Zhenlei Cai和Anthony T.S.在文献“Cloud detection and removal in satellite images for tropical regions”（1996 3rd International Conference on Signal Processing Proceedings,1996,Vol.2,Pages:1174-1177）中先对图像做非监督聚类，人工选择最佳聚类结果并从中找出云和阴影类，再通过云和阴影的位置关系进行云影配对以去除误检。该方法对聚类结果要求较高且需要人工参与，而且云和阴影的配对准则不完备导致配对稳定性不高，因而检测准确性受到影响。2010年Jianzhong Feng、Linyan Bai、Huajun Tang等在文献“A new context-based procedure for the detection and removal of cloud shadow from moderate-and-high resolution satellite data over land”（Proceedings 2010 IEEE International Geoscience and Remote Sensing Symposium2010,Pages:1718-1721）中针对FY-3-A图像和HJ-1-A图像特点，通过计算其中某个 特定可见光波段与近红外波段图像的比值，采用阈值分割将厚云和阴影粗略检测出来，然后利用太阳、卫星和厚云的空间关系在灰度直方图上反映的特征进行进一步确认。该方法在考虑光谱信息的同时还结合了成像时的空间关系以及像素邻域信息，但算法存在缩放效应，并且只适用于FY-3-A或HJ-1-A图像，不能直接用于专题制图仪TM图像或增强专题制图仪ETM图像的厚云及阴影检测。 

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种针对Landsat TM和ETM图像的厚云及其阴影检测方法，以消除对辅助信息的需求，提高自动化程度和检测准确度。 

本发明是这样实现的： 

实现本发明目的的技术思路是：根据专题制图仪图像或增强专题制图仪图像在除第六波段外的第一至第七波段厚云和冰雪的灰度值大，而阴影的灰度值小的特点，分别检测出厚云、冰雪和阴影；利用厚云及冰雪的饱和度极低这一性质检测出厚云和冰雪；根据厚云和阴影在各波段间的灰度值变化量小的特点检出云和阴影，从而排除冰雪。利用厚云及阴影的伴随关系对厚云和阴影进行配对，以除去误检，并通过对未配对成功的厚云和阴影进行补充检测，消除漏检。具体步骤包括如下： 

（1）输入一幅包含除第六波段外的第一至第七波段专题制图仪TM图像或增强专题制图仪ETM图像，以输入图像长的1/4和宽的1/4分别作为要获得子图像的长和宽，将输入图像分为互不重叠的16组子图像，用每一组子图像包含的6个波段图像Y_ij ^b构成一个子图集I_ij，其中i、j分别为行和列方向分成子图的序号，i＝{1,2,3,4}，j＝{1,2,3,4}，b为波段序号，b＝{1,2,3,4,5,6}； 

（2）将各子图集中的各波段子图像Y_ij ^b分别进行维纳滤波去噪，并归一化处理，得到归一化后子图像 

和归一化后子图集 

（3）将归一化后子图集 

所有波段图像 

的相同空间位置处的像素灰度值相加并求平均值，得到一幅灰度平均图像E_ij，用云阈值T₁对E_ij做二值化分割，即把E_ij中灰度值大于T₁的像素灰度值置为1，其余像素灰度值置为0，得到 

的云初检二值图CL1_ij；再对E_ij用阴影阈值T₂进行分割，把E_ij中灰度值小于T₂的像素灰度值置为0.5，其余像素灰度值置为0，得到 

的阴影初检二值图SH_ij，其中T₁=0.8，T₂=0.1； 

（4）对归一化后子图集 

所有波段图像 

的相同空间位置处的像素灰度值求方 差，得到一幅方差图像V_ij，用云影阈值T₃对V_ij进行二值化分割，把V_ij中灰度值小于T₃的像素灰度值置为1，其余像素灰度值置为0，得到 

的云影二值图CS'_ij，其中T₃=0.002； 

（5）将归一化后子图集 

中的 和 

分别作为蓝、绿和红三个颜色分量，得到一幅红绿蓝RGB彩色合成图像，然后将该图像从红-绿-蓝RGB色彩空间转换到色相-饱和度-明度HSV色彩空间，进而得到一幅饱和度图像S_ij；再用饱和度阈值T₄对S_ij进行二值化分割，把S_ij中灰度值小于T₄的像素灰度值置为1，其余像素灰度值置为0，得到 

的饱和度云检测图CL2_ij，其中T₄=0.02； 

（6）将归一化后子图集 

的云初检二值图CL1_ij、云影二值图CS'_ij以及饱和度云检测图CL2_ij这三者相同空间位置处的像素灰度值相乘，得到云种子图，并将该图像中的云种子像素点在饱和度图像S_ij中进行区域生长，得到 

的云粗检结果图Cd_ij，其中区域生长使用邻域判断阈值T，T=0.03； 

（7）将归一化后子图集 

的阴影初检二值图SH_ij和云影二值图CS'_ij相同空间位置处的像素灰度值相乘，得到阴影种子图，并将该图像中的阴影种子像素点在灰度平均图E_ij中进行区域生长，得到 

的阴影粗检结果图Sd_ij，其中区域生长使用邻域判断阈值T，T=0.03； 

（8）在8邻域条件下，用半径为两个像素的圆盘形结构元素分别对归一化后子图集 

的云粗检结果图Cd_ij和阴影粗检结果图Sd_ij做一次数学形态学闭运算，后将该两幅图中封闭区块面积小于8个像素的区块从图像中剔除，并对Cd_ij和Sd_ij的对应像素灰度值相加，得到 

的云影粗检测图CS_ij； 

（9）分别从每个归一化后子图集 

的云影粗检测图CS_ij中选取该 的基准对，共得到Ж个基准对： 

（9a）将归一化后子图集 的云影粗检测图CS_ij中包含的所有M个云块与所有N个阴影块一一组成云影对，共M×N个；如果M=0或N=0，则该归一化后子图集 

无基准对，否则统计每个云影对中云块的面积S_C和周长L_C，阴影块的面积S_S和周长L_S、云影对质心连线长度d以及倾角θ； 

（9b）判断每一个云影对是否满足1）式： 

$\left\{\begin{array}{c}100\≤S_C\≤900\\ 100\≤S_S\≤900\\ |S_C-S_S|\≤\frac{\mathrm{\α}}{2}(S_C+S_S)\\ |L_C-L_S|\≤\frac{\mathrm{\β}}{2}(L_C+L_S)\\ d\≤\mathrm{\γ}\sqrt{0.5(S_C+S_S)}\end{array}\right.;---1)$

若满足1）式，则将该云影对作为一个候选基准对，并在所有云影对判断完毕后进行步骤（9e），若没有一个云影对满足1）式，则进行步骤（9c），其中1）式中α为面积约束系数，β为周长约束系数，γ为间距约束系数，α＝0.3，β＝0.25，γ＝3； 

（9c）将1）式中面积约束系数α、周长约束系数β和间距约束系数γ同时各增大百分之一，再重新判断是否有云影对满足1）式，如有云影对满足1）式则将其作为一个候选基准对，并在所有云影对判断完毕后进行步骤（9e）；否则进行步骤（9d）； 

（9d）如果α<1、β＜1且γ<5，则转到步骤（9c）；否则，则该归一化后子图集 

无基准对； 

（9e）在所有得到的候选基准对中选择云影对质心连线长度d最小、阴影面积S_S最大的候选基准对作为归一化后子图集 的基准对； 

（9f）分别对16个归一化后子图集 

重复步骤（9a）至（9e），得到Ж个基准对，其中Ж≤16； 

（10）统计上述Ж个基准对的质心连线主角度和质心连线主长度，得到最终基准对倾角A和最终基准对间距D； 

（11）根据最终基准对倾角A和间距D对归一化后子图集 

中的所有云块和阴影块进行配对，得到 

的云影配对结果图P_ij； 

（12）将归一化后子图集 

的云影配对结果图P_ij与云影粗检测图CS_ij对应位置处的像素灰度值相减并求绝对值，得到 

的未配对云影块图像； 

（13）根据最终基准对倾角A和间距D对归一化后子图集 

的未配对云影块图像中的云块补充检测对应阴影块，得到 

的阴影补充检测结果图AD1_ij； 

（14）根据最终基准对倾角A和间距D对归一化后子图集 

的未配对云影块图像中的阴影块补充检测对应云块，得到 

的云补充检测结果图AD2_ij； 

（15）将归一化后子图集 

的云影配对结果图P_ij、云补充检测结果图AD1_ij和阴影补充检测结果图AD2_ij的对应像素相加，得到 

的最终检测结果子图CaS_ij； 

（16）将16个归一化后子图集 

对应的16个最终检测结果子图CaS_ij按顺序拼接，得到一幅完整的最终检测结果图。 

上述针对Landsat TM和ETM图像的厚云及其阴影检测方法，其中所述步骤（10）中统计上述Ж个基准对的质心连线主角度和质心连线主长度，按如下步骤进行： 

（10a）设角度阈值Θ=20，并设基准对的序号为Д，其中1≤Д≤Ж，此处令Д=1； 

（10b）选择第Д个基准对； 

（10c）将选中的基准对记为Δ，将Δ的质心连线倾角分别与其余Ж-1个基准对的质心连线倾角作差并取绝对值，得到Ж-1个基准对的Ж-1个差异角；如果这Ж-1个差异角中有超过(Ж-1)/2个小于角度阈值Θ，则将所有小于角度阈值Θ的差异角对应基准对的质心连线倾角与Δ的质心连线倾角相加并求平均值，得到质心连线主角度，并将该质心连线主角度作为最终基准对倾角A，将所有小于角度阈值Θ的差异角对应基准对的质心连线长度与Δ的质心连线长度相加并求平均值，得到质心连线主长度，并将该质心连线主长度作为最终基准对间距D；否则执行步骤（10d）； 

（10d）若Д<Ж，则将Д的值加1，并转到步骤（10b）；若Д=Ж，并仍未得到最终基准对倾角A和最终基准对间距D，则将角度阈值Θ的值增加5且令Д=1，并转到步骤（10b）。 

上述针对Landsat TM和ETM图像的厚云及其阴影检测方法，其中步骤（11）所述的根据最终基准对倾角A和间距D对归一化后子图集 

中的所有云块和阴影块进行配对，是根据云影对是否满足如下公式进行： 

$\left\{\begin{array}{c}A-\mathrm{\&delta;}\&le;\mathrm{\&theta;}\&le;A+\mathrm{\&delta;}\\ D-\mathrm{\&xi;}\&le;d\&le;D+\mathrm{\&xi;}\\ S_S\&le;a{S}_C\\ L_S\&le;b{L}_C\end{array}\right.;---2)$

若在归一化后子图集 

的所有M×N个云影对中有任何一个满足2）式，则将该云影对在一幅与 的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像的相应位置处保存，即将云对应像素的灰度值置为1，将阴影对应像素的灰度值置为0.5，并待所有云影对判断完毕后，将保存云影对的所有这些图像对应位置的像素灰度值相加，得到 

的云影配对结果图P_ij；若没有任何一个云影对满足2）式或 的云影粗检测图CS_ij中包含的云块个数M=0或阴影块个数N=0，则将一幅与 

的云影粗 检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像作为 的云影配对结果图P_ij，其中a为面积比系数，b为周长比系数，δ为角度松弛变量，ξ为间距松弛变量，a＝10，b＝5，δ＝20，ζ＝20。 

上述针对Landsat TM和ETM图像的厚云及其阴影检测方法，其中步骤（13）所述的根据最终基准对倾角A和间距D对归一化后子图集 

的未配对云影块图像中的云块补充检测对应阴影块，通过如下步骤进行： 

（13a）设M_X为归一化后子图集 

的未配对云影块图像中云块的个数，其中M_X≥0，若M_X等于0，则将一幅与 

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像作为 的阴影补充检测结果图AD1_ij；否则设m为 

的未配对云影块图像中云块的序号，其中1≤m≤M_X，此处令m＝1； 

（13b）根据第m个云块的质心像素坐标(X_cm，Y_cm)，利用3）式求出对应阴影块的质心像素坐标(X_sm，Y_sm)； 

$\left\{\begin{array}{c}{X}_{\mathrm{sm}}=X_{\mathrm{cm}}-D*\cos \left(A\right)\\ Y_{\mathrm{sm}}=Y_{\mathrm{cm}}-D*\sin \left(A\right)\end{array}\right.;---3)$

（13c）在 的灰度平均图E_ij上选取以像素坐标(X_sm，Y_sm)为中心、K×K大小的窗，然后将该窗口内灰度值最小的像素点作为种子点，利用区域生长法在E_ij中进行区域生长，得到补充检测的阴影块；如果补充检测的阴影块其面积和周长分别超过对应云块面积和周长的 

倍和ψ倍，则将该补充检测的阴影块剔除，否则保留，其中区域生长使用邻域判断阈值T，T=0.03， 

ψ＝5，K＝9； 

（13d）若m＜M_X，则将m的值加1，并转到步骤（13b）；若m＝M_X，则将得到的所有补充检测的阴影块与对应的云块在一幅与 

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像的相应位置处保存，即将云对应像素的灰度值置为1，将阴影对应像素的灰度值置为0.5，得到归一化后子图集 

的阴影补充检测结果图AD1_ij。 

上述针对Landsat TM和ETM图像的厚云及其阴影检测方法，其中步骤（14）所述的根据最终基准对倾角A和间距D对归一化后子图集 

的未配对云影块图像中的阴影块补充检测对应云块，通过如下步骤进行： 

（14a）设N_X为归一化后子图集 的未配对云影块图像中阴影块的个数，其中N_X≥0，若N_X等于0，则将一幅与 

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全 为0的空白图像作为 

的云补充检测结果图AD2_ij；否则设n为 

的未配对云影块图像中阴影块的序号，其中1≤n≤N_X，此处令n＝1； 

（14b）根据第n个阴影块的质心像素坐标(X_sn，Y_sn)，利用4）式求出对应云块质心像素坐标(X_cn， Y_cn)； 

$\left\{\begin{array}{c}{X}_{\mathrm{cn}}=X_{\mathrm{sn}}-D*\cos (A+180)\\ Y_{\mathrm{cn}}=Y_{\mathrm{sn}}-D*\sin (A+180)\end{array}\right.;---4)$

（14c）在 

的灰度平均图E_ij上选取以像素坐标(X_cn，Y_cn)为中心、K×K大小的窗，然后将该窗口内灰度值最大的像素点作为种子点，利用区域生长法在E_ij中进行区域生长，得到补充检测的云块；如果补充检测的云块其面积和周长分别超过对应阴影块面积和周长的 

倍和ψ倍，则同时剔除该补充检测的云块和它对应的阴影块，否则都予以保留，其中区域生长使用邻域判断阈值T，T=0.03， 

ψ＝5，K＝9； 

（14d）若n＜N_X，则将n的值加1，并转到步骤（14b）；若n＝N_X，则将得到的所有补充检测的云块与对应的阴影块在一幅与 

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像的相应位置处保存，即将云对应像素的灰度值置为1，将阴影对应像素的灰度值置为0.5，得到归一化后子图集 

的云补充检测结果图AD2_ij。 

本发明与现有的技术相比具有以下优点： 

a、本发明与现有的针对专题制图仪TM图像或增强专题制图仪ETM图像数据的光谱分析法相比，由于只利用了图像本身进行厚云及其阴影的检测，因此无需辅助信息；而且由于整个检测过程不需要人工参与，自动化程度较高； 

b、本发明与现有的针对专题制图仪TM图像或增强专题制图仪ETM图像数据的光谱分析法相比，由于利用了厚云及其阴影的空间上下文关系信息进行云影配对和补充检测，因而降低了误检和漏检的发生，检测准确度得到提高。 

附图说明

图1是本发明的整体流程图； 

图2是本发明的实验图像； 

图3是本发明的厚云及其阴影检测结果图。 

具体实施方式

参照图1，本发明的实现步骤如下： 

步骤1，输入图像，将其分为16组子图像，并构成16个子图集。 

一景完整的专题制图仪TM图像包括第一至第七波段共7个波段图像，而增强专题制图仪ETM图像包括第一至第八波段共8个波段图像。由于TM图像中的第六波段图像与其他波段图像分别率不同，而ETM图像中的第六和第八波段图像也与其他波段图像分辨率不同，因此本发明将一景专题制图仪TM图像或增强专题制图仪ETM图像中的第一至第五波段以及第七波段图像作为输入图像X，X={x₁,x₂,x₃,x₄,x₅,x₆}，其中x₁至x₅分别对应第一至第五波段图像，x₆对应第七波段图像。输入图像X的各波段图像x₁～x₆其图像尺寸大小相等，并且各波段图像相同空间位置处的像素对应相同的地理位置。 

以输入图像X中任意一个波段图像长的1/4和宽的1/4分别作为要获得子图像的长和宽，将X中的每一波段图像都分为互不重叠的16个子图像，其中同一空间位置对应的6个波段子图像视为一组，共得到16组子图像；将每一组子图像包含的6个波段子图像Y_ij ^b构成一个子图集I_ij，I_ij＝{Y_ij ¹,Y_ij ²，Y_ij ³，Y_ij ⁴，Y_ij ⁵，Y_ij ⁶}，其中i和j分别为行方向和列方向分成子图的序号，i＝{1,2,3,4}，j＝{1,2,3,4}，b为波段序号，b＝{1,2,3,4,5,6}。 

步骤2，将各子图集中的各波段子图像进行维纳滤波去噪，并归一化处理，得到归一化后子图像 

和归一化后子图集 

2.1）将子图集I_ij中每一幅子图像Y_ij ^b的灰度值区间从0~255转换到0~1之间，然后采用窗口大小为3×3像素的维纳滤波去噪处理，得到去噪后子图像 

2.2）对去噪后子图像 

按如下公式进行归一化处理，得到归一化后子图像 

并构成归一化后子图集 

即 ${\hat{I}}_{\mathrm{ij}}=\{{{\hat{Y}}_{\mathrm{ij}}}^1,{{\hat{Y}}_{\mathrm{ij}}}^2,{{\hat{Y}}_{\mathrm{ij}}}^3,{{\hat{Y}}_{\mathrm{ij}}}^4,{{\hat{Y}}_{\mathrm{ij}}}^5,{{\hat{Y}}_{\mathrm{ij}}}^6\},$

${{\hat{Y}}_{\mathrm{ij}}}^b=\frac{{{\stackrel{\&CenterDot;\&CenterDot;}{Y}}_{\mathrm{ij}}}^b-\min \left({{\stackrel{\&CenterDot;\&CenterDot;}{Y}}_{\mathrm{ij}}}^b\right)}{\max \left({{\stackrel{\&CenterDot;\&CenterDot;}{Y}}_{\mathrm{ij}}}^b\right)-\min \left({{\stackrel{\&CenterDot;\&CenterDot;}{Y}}_{\mathrm{ij}}}^b\right)},$

其中， 

表示取 

中最大的灰度值， 表示取 

中最小的灰度值。 

步骤3，将归一化后子图集 中所有波段图像的对应空间位置处的像素灰度值相加并求平均，得到 

的灰度平均图像： 

步骤4，分别用云阈值和阴影阈值对 

的灰度平均图像E_ij进行二值化分割，得到 

的云初检二值图CL1_ij和阴影初检二值图SH_ij。 

4.1）根据厚云在各波段图像上的灰度值处于灰度区间的极大值一端，经相加平均 后基本处于0.75到1之间的特性，设云阈值T₁的取值范围为0.8到0.85之间，本实例中取T₁=0.8； 

4.2）根据阴影在各波段图像上的灰度值处于灰度区间的极小值一端，经相加平均后基本处于0到0.2之间的特性，设阴影阈值T₂的取值范围为0.1到0.15之间，本实例中取T₂=0.1； 

4.3）用云阈值T₁对 的灰度平均图E_ij做二值化分割，即把E_ij中灰度值大于T₁的像素灰度值置为1，其余像素灰度值置为0，得到 的云初检二值图CL1_ij； 

4.4）用阴影阈值T₂对 

的灰度平均图E_ij做二值化分割，即把E_ij中灰度值小于T₂的像素灰度值置为0.5，其余像素灰度值置为0，得到 

的阴影初检二值图SH_ij。 

步骤5，由归一化后子图集 

中所有波段图像的对应空间位置处的像素灰度值求方差，得到 

的方差图像： 

步骤6，用云影阈值对 

的方差图像V_ij进行二值化分割，得到 

的云影二值图CS'_ij。 

6.1）根据实验统计得出的厚云及其阴影在各波段图像上的灰度值变化量极小的特性，设云影阈值T₃的取值范围为0.002到0.004之间，本实例中取T₃=0.002； 

6.2）用云影阈值T₃对 

的方差图像V_ij进行二值化分割，即把V_ij中灰度值小于T₃的像素灰度值置为1，其余像素灰度值置为0，得到 

的云影二值图CS'_ij。 

步骤7，由于第一波段图像的波谱范围为0.45~0.52um，对应蓝光波段，第二波段图像的波谱范围为0.52~0.60um，对应绿光波段，第三波段图像的波谱范围为0.63~0.69um，对应红光波段，因此将归一化后子图集 

中的 和 

分别作为红-绿-蓝RGB彩色图像中的蓝、绿和红三个颜色分量，得到一幅 

的红-绿-蓝RGB彩色合成图像。 

步骤8，将 

的红-绿-蓝RGB彩色合成图像的色彩空间由红-绿-蓝RGB转换到色相-饱和度-明度HSV，得到 

的饱和度图像： 

其中max(R,G,B)表示取组成该像素点的三个颜色分量中最大的灰度值，反之，min(R,G,B)表示取组成该像素点的三个颜色分量中最小的灰度值。 

步骤9，用 

的饱和度阈值对饱和度图像S_ij进行二值化分割，得到 

的饱和度云检测图CL2_ij。 

9.1）根据厚云的饱和度比一般地面目标都要低的特性，设饱和度阈值T₄的取值范围为0.02到0.03之间，本实例中取T₄=0.02； 

9.2）用饱和度阈值T₄对 

的饱和度图像S_ij进行二值化分割，即将S_ij中灰度值小于T₄的像素灰度值置为1，其余像素灰度值置为0，得到 

的饱和度云检测图CL2_ij。 

步骤10，将归一化后子图集 的云初检二值图CL1_ij、云影二值图CS'_ij以及饱和度云检测图CL2_ij的对应空间位置处的像素灰度值相乘，得到 

的云种子图。 

步骤11，利用区域生长法得到 

的云粗检结果图Cd_ij。 

11.1）设区域生长的邻域判断阈值T的取值范围为0.02到0.04之间，本实例中取T=0.03； 

11.2）把 的云种子图中像素灰度值为1的像素点作为云种子点； 

11.3）利用区域生长方法将每一个云种子点在 

的饱和度图像S_ij中进行区域生长，得到 的云粗检结果图Cd_ij： 

11.3a）任选一个云种子点； 

11.3b）判断以选中的云种子点在 

的饱和度图像S_ij中对应像素为中心的8邻域像素点分别与中心像素点的灰度值之差的绝对值是否满足小于邻域判断阈值T；如果满足，且该邻域像素点在云种子图中对应的像素点不是云种子点，则将该像素点作为一个新的云种子点加入云种子图；否则转到步骤11.3a），为了防止加入错误的云种子点，要求新加入的云种子点在 

的饱和度图像S_ij中对应像素的灰度值必须小于等于λ₁T₄，其中T₄为饱和度阈值，λ₁为T₄的约束系数，λ₁=3； 

11.3c）对云种子图中的每一个云种子点和新加入的云种子点重复步骤11.3a）和11.3b），直到不再有新的云种子点加入云种子点图为止，则经过这一过程加入了新的云种子点的云种子点图，即为所要求得的 

的云粗检结果图Cd_ij。 

步骤12，将归一化后子图集 

的阴影初检二值图SH_ij和云影二值图CS'_ij的对应空间位置处的像素灰度值相乘，得到 

的阴影种子图。 

步骤13，利用区域生长法得到 的阴影粗检结果图Sd_ij。 

13.1）设区域生长的邻域判断阈值T的取值范围为0.02到0.04之间，本实例中取T=0.03； 

13.2）把 

的阴影种子图中像素灰度值为0.5的像素点作为阴影种子点； 

13.3）利用区域生长方法将每一个阴影种子点在 的灰度平均图像E_ij中进行区域 生长，得到 

的阴影粗检结果图Sd_ij； 

13.3a）任选一个阴影种子点； 

13.3b）判断以选中的阴影种子点在 的灰度平均图像E_ij中对应像素为中心的8邻域像素点分别与中心像素点的灰度值之差的绝对值是否满足小于邻域判断阈值T；如果满足，且该邻域像素点在阴影种子图中对应的像素点不是阴影种子点，则将该像素点作为一个新的阴影种子点加入阴影种子图；否则转到步骤13.3a），为了防止加入错误的阴影种子点，要求新加入的阴影种子点在E_ij中对应像素的灰度值必须小于等于λ₂T₂，其中T₂为阴影阈值，λ₂为T₂的约束系数，λ₂=1； 

13.3c）对阴影种子图中的每一个阴影种子点和新加入的阴影种子点重复步骤13.3a）到13.3b），直到不再有新的阴影种子点加入为止，则经过这一过程加入了新的阴影种子点的阴影种子点图，即为所要求得的 

的阴影粗检结果图Sd_ij。 

步骤14，在8邻域条件下，用半径为两个像素的圆盘形结构元素分别对归一化后子图集 

的云粗检结果图Cd_ij和阴影粗检结果图Sd_ij做一次数学形态学闭运算，后将该两幅图中封闭区块面积小于8个像素的区块从图像中剔除，并对Cd_ij和Sd_ij的对应像素灰度值相加，得到 

的云影粗检测图CS_ij。 

步骤15，分别从每个归一化后子图集 

对应的云影粗检测图CS_ij中选取该 

的基准对，共得到Ж个基准对。 

15.1）归一化后子图集 的云影粗检测图CS_ij中包含M个灰度值为1的封闭区块和N个灰度值为0.5的封闭区块，它们分别代表云块和阴影块。将所有M个云块与所有N个阴影块一一组成云影对，共M×N个；如果M=0或N=0，则该归一化后子图集 

无基准对，否则统计每个云影对中云块的面积S_C和周长L_C，阴影块的面积S_S和周长L_S、云影对质心连线长度d以及倾角θ，其中S_C、S_S、L_C、 L_S和d都以像素为单位； 

15.2）根据厚云和其产生的阴影具有相似形状的特性，提出判断每一个云影对是否满足1）式： 

$\left\{\begin{array}{c}100\&le;S_C\&le;900\\ 100\&le;S_S\&le;900\\ |S_C-S_S|\&le;\frac{\mathrm{\&alpha;}}{2}(S_C+S_S)\\ |L_C-L_S|\&le;\frac{\mathrm{\&beta;}}{2}(L_C+L_S)\\ d\&le;\mathrm{\&gamma;}\sqrt{0.5(S_C+S_S)}\end{array}\right.;---1)$

若满足1）式，则将该云影对作为一个候选基准对，并在所有云影对判断完毕后进行步骤15.5），若没有一个云影对满足1）式，则进行步骤15.3），其中1）式中α为面积约束系数，β为周长约束系数，γ为间距约束系数，α＝0.3，β＝0.25，γ＝3； 

15.3）将1）式中面积约束系数α、周长约束系数β和间距约束系数γ同时各增大百分之一，再重新判断是否有云影对满足1）式，如有云影对满足1）式则将其作为一个候选基准对，并在所有云影对判断完毕后进行步骤15.5）；否则进行步骤15.4）； 

15.4）如果α<1、β＜1且γ<5，则转到步骤15.3）；否则，则该归一化后子图集 

无基准对； 

15.5）在所有得到的候选基准对中选择云影对质心连线长度d最小、阴影面积S_S最大的候选基准对作为归一化后子图集 

的基准对； 

15.6）分别对16个归一化后子图集 

重复步骤15.1）至15.5），得到Ж个基准对，其中Ж≤16； 

步骤16，统计上述Ж个基准对的质心连线主角度和质心连线主长度，得到最终基准对倾角A和最终基准对间距D。 

16.1）设角度阈值Θ=20，并设基准对的序号为Д，其中1≤Д≤Ж，此处令Д=1； 

16.2）选择第Д个基准对； 

16.3）将选中的基准对记为Δ，将Δ的质心连线倾角分别与其余Ж-1个基准对的质心连线倾角作差并取绝对值，得到Ж-1个基准对的Ж-1个差异角；如果这Ж-1个差异角中有超过(Ж-1)/2个小于角度阈值Θ，则将所有小于角度阈值Θ的差异角对应基准对的质心连线倾角与Δ的质心连线倾角相加并求平均值，得到质心连线主角度，并将该质心连线主角度作为最终基准对倾角A，将所有小于角度阈值Θ的差异角对应基准对的质心连线长度与Δ的质心连线长度相加并求平均值，得到质心连线主长度，并将该质心连线主长度作为最终基准对间距D；否则执行步骤16.4）； 

16.4）若Д<Ж，则将Д的值加1，并转到步骤16.2）；若Д=Ж，并仍未得到最终 基准对倾角A和最终基准对间距D，则将角度阈值Θ的值增加5且令Д=1，并转到步骤16.2）。 

步骤17，根据最终基准对倾角A和间距D对归一化后子图集 

中的所有云块和阴影块进行配对，得到 

的云影配对结果图P_ij。 

对于一幅图像中的每一个像素而言，拍摄时的太阳高度角和拍摄角度是固定不变的，因此图像中所有厚云及其产生的阴影具有一致的相对空间位置。根据这一特点，提出判断归一化后子图集 中的所有云影对是否满足如下公式： 

$\left\{\begin{array}{c}A-\mathrm{\&delta;}\&le;\mathrm{\&theta;}\&le;A+\mathrm{\&delta;}\\ D-\mathrm{\&xi;}\&le;d\&le;D+\mathrm{\&xi;}\\ S_S\&le;a{S}_C\\ L_S\&le;b{L}_C\end{array}\right.;---2)$

若在所有M×N个云影对中有任何一个满足2）式，则将该云影对在一幅与 

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像的相应位置处保存，即将云对应像素的灰度值置为1，将阴影对应像素的灰度值置为0.5，并待所有云影对判断完毕后，将保存云影对的所有这些图像对应位置的像素灰度值相加，得到 

的云影配对结果图P_ij；若没有任何一个云影对满足2）式或 

的云影粗检测图CS_ij中包含的云块个数M=0或阴影块个数N=0，则将一幅与 

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像作为 的云影配对结果图P_ij，其中a为面积比系数，b为周长比系数，δ为角度松弛变量，ξ为间距松弛变量，a＝10，b＝5，δ＝20，ζ＝20。 

步骤18，将归一化后子图集 

的云影配对结果图P_ij与云影粗检测图CS_ij对应位置处的像素灰度值相减并求绝对值，得到 

的未配对云影块图像，其中像素灰度值为1的封闭区域代表云块，像素灰度值为0.5的封闭区域代表阴影块。 

步骤19，根据最终基准对倾角A和间距D对归一化后子图集 

的未配对云影块图像中的云块补充检测对应阴影块，得到 

的阴影补充检测结果图AD1_ij。 

19.1）设M_X为归一化后子图集 

的未配对云影块图像中云块的个数，其中M_X≥0，若M_X等于0，则将一幅与 的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像作为 

的阴影补充检测结果图AD1_ij；否则设m为 

的未配对云影块图像中云块的序号，其中1≤m≤M_X，此处令m＝1； 

19.2）根据第m个云块的质心像素坐标(X_cm，Y_cm)，利用3）式求出对应阴影块的质心像素坐标(X_sm，Y_sm)； 

$\left\{\begin{array}{c}{X}_{\mathrm{sm}}=X_{\mathrm{cm}}-D*\cos \left(A\right)\\ Y_{\mathrm{sm}}=Y_{\mathrm{cm}}-D*\sin \left(A\right)\end{array}\right.;---3)$

19.3）在 

的灰度平均图E_ij上选取以像素坐标(X_sm，Y_sm)为中心、K×K大小的窗，然后将该窗口内灰度值最小的像素点作为种子点，利用区域生长法在E_ij中进行区域生长，得到补充检测的阴影块；如果补充检测的阴影块其面积和周长分别超过对应云块面积和周长的 

倍和ψ倍，则将该补充检测的阴影块剔除，否则保留，其中区域生长使用邻域判断阈值T，T=0.03， 

ψ＝5，K＝9； 

19.4）若m＜M_X，则将m的值加1，并转到步骤19.2）；若m＝M_X，则将得到的所有补充检测的阴影块与对应的云块在一幅与 

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像的相应位置处保存，即将云对应像素的灰度值置为1，将阴影对应像素的灰度值置为0.5，得到归一化后子图集 

的阴影补充检测结果图AD1_ij。 

步骤20，根据最终基准对倾角A和间距D对归一化后子图集 

的未配对云影块图像中的阴影块补充检测对应云块，得到 

的云补充检测结果图AD2_ij。 

20.1）设N_X为归一化后子图集 

的未配对云影块图像中阴影块的个数，其中N_X≥0，若N_X等于0，则将一幅与 

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像作为 

的云补充检测结果图AD2_ij；否则设n为 

的未配对云影块图像中阴影块的序号，其中1≤n≤N_X，此处令n＝1； 

20.2）根据第n个阴影块的质心像素坐标(X_sn，Y_sn)，利用4）式求出对应云块质心像素坐标(X_cn，Y_cn)； 

$\left\{\begin{array}{c}{X}_{\mathrm{cn}}=X_{\mathrm{sn}}-D*\cos (A+180)\\ Y_{\mathrm{cn}}=Y_{\mathrm{sn}}-D*\sin (A+180)\end{array}\right.;---4)$

20.3）在 

的灰度平均图E_ij上选取以像素坐标(X_cn，Y_cn)为中心、K×K大小的窗，然后将该窗口内灰度值最大的像素点作为种子点，利用区域生长法在E_ij中进行区域生长，得到补充检测的云块；如果补充检测的云块其面积和周长分别超过对应阴影块面积和周长的 

倍和ψ倍，则同时剔除该补充检测的云块和它对应的阴影块，否则都予以保留，其中区域生长使用邻域判断阈值T，T=0.03， 

ψ＝5，K＝9； 

20.4）若n＜N_X，则将n的值加1，并转到步骤20.2）；若n＝N_X，则将得到的所有补充检测的云块与对应的阴影块在一幅与 的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像的相应位置处保存，即将云对应像素的灰度值置为1，将阴 影对应像素的灰度值置为0.5，得到归一化后子图集 

的云补充检测结果图AD2_ij。 

步骤21，将归一化后子图集 

的云影配对结果图P_ij、云补充检测结果图AD1_ij和阴影补充检测结果图AD2_ij的对应空间位置处的像素灰度值相加，得到 

的最终检测结果子图CaS_ij。 

步骤22，将16个归一化后子图集 

对应的16个最终检测结果子图CaS_ij按顺序拼接，得到一幅完整的最终检测结果图CaS。 

本发明的效果可通过如下仿真实验结果进一步说明： 

1.实验数据 

实验所使用的输入图像如图2所示，其中图2（a）至图2（f）分别为除第六、第八波段外的第一至第七波段的增强专题制图仪ETM图像，该图像由一景L1G原始数据通过旋转得到，幅面均为5800×5600像素，灰度级为256级。 

2.厚云及其阴影检测内容及结果 

用本发明对图2所示的实验数据进行厚云及其阴影检测，得到如图3（a）~（e）所示的实验结果，其中图3（a）为云影粗检测图CS；图3（b）为云影配对结果图P；图3（c）为阴影补充检测结果图AD1；图3（d）为云补充检测结果图AD2；图3（e）为最终检测结果图CaS，图中白色部分代表云，灰色部分代表阴影，黑色部分代表地面。 

结论：通过将图2所示的实验图像与图3（e）所示的最终检测结果图进行对比观察发现，图2中包含的厚云及其阴影基本都被检测出来，并且几乎没有错误地将地面物体当作厚云及其阴影检测出来，可见本发明的检测准确性较高。 

Claims (5)

1.一种针对Landsat TM或ETM图像的厚云及其阴影检测方法，包括步骤如下：

（1）输入一幅包含除第六波段外的第一至第七波段专题制图仪TM图像或增强专题制图仪ETM图像，以输入图像长的1/4和宽的1/4分别作为要获得子图像的长和宽，将输入图像分为互不重叠的16组子图像，用每一组子图像包含的6个波段图像Y_ij ^b构成一个子图集I_ij，其中i、j分别为行和列方向分成子图的序号，i＝{1,2,3,4}，j＝{1,2,3,4}，b为波段序号，b＝{1,2,3,4,5,6}；

（2）将各子图集中的各波段子图像Y_ij ^b分别进行维纳滤波去噪，并归一化处理，得到归一化后子图像和归一化后子图集

（3）将归一化后子图集

所有波段图像的相同空间位置处的像素灰度值相加并求平均值，得到一幅灰度平均图像E_ij，用云阈值T₁对E_ij做二值化分割，即把E_ij中灰度值大于T₁的像素灰度值置为1，其余像素灰度值置为0，得到的云初检二值图CL1_ij；再对E_ij用阴影阈值T₂进行分割，把E_ij中灰度值小于T₂的像素灰度值置为0.5，其余像素灰度值置为0，得到

的阴影初检二值图SH_ij，其中T₁=0.8，T₂=0.1；

（4）对归一化后子图集

所有波段图像

的相同空间位置处的像素灰度值求方差，得到一幅方差图像V_ij，用云影阈值T₃对V_ij进行二值化分割，把V_ij中灰度值小于T₃的像素灰度值置为1，其余像素灰度值置为0，得到

的云影二值图CS'_ij，其中T₃=0.002；

（5）将归一化后子图集

中的

和

分别作为蓝、绿和红三个颜色分量，得到一幅红绿蓝RGB彩色合成图像，然后将该图像从红-绿-蓝RGB色彩空间转换到色相-饱和度-明度HSV色彩空间，进而得到一幅饱和度图像S_ij；再用饱和度阈值T₄对S_ij进行二值化分割，把S_ij中灰度值小于T₄的像素灰度值置为1，其余像素灰度值置为0，得到

的饱和度云检测图CL2_ij，其中T₄=0.02；

（6）将归一化后子图集的云初检二值图CL1_ij、云影二值图CS'_ij以及饱和度云检测图CL2_ij这三者相同空间位置处的像素灰度值相乘，得到云种子图，并将该图像中的云种子像素点在饱和度图像S_ij中进行区域生长，得到的云粗检结果图Cd_ij，其中区域生长使用邻域判断阈值T，T=0.03；

（7）将归一化后子图集

的阴影初检二值图SH_ij和云影二值图CS'_ij相同空间位置处的像素灰度值相乘，得到阴影种子图，并将该图像中的阴影种子像素点在灰度平均图E_ij中进行区域生长，得到

的阴影粗检结果图Sd_ij，其中区域生长使用邻域判断阈值T，T=0.03；

（8）在8邻域条件下，用半径为两个像素的圆盘形结构元素分别对归一化后子图集的云粗检结果图Cd_ij和阴影粗检结果图Sd_ij做一次数学形态学闭运算，后将该两幅图中封闭区块面积小于8个像素的区块从图像中剔除，并对Cd_ij和Sd_ij的对应像素灰度值相加，得到

的云影粗检测图CS_ij；

（9）分别从每个归一化后子图集的云影粗检测图CS_ij中选取该

的基准对，共得到Ж个基准对：

（9a）将归一化后子图集的云影粗检测图CS_ij中包含的所有M个云块与所有N个阴影块一一组成云影对，共M×N个；如果M=0或N=0，则该归一化后子图集无基准对，否则统计每个云影对中云块的面积S_C和周长L_C，阴影块的面积S_S和周长L_S、云影对质心连线长度d以及倾角θ；

（9b）判断每一个云影对是否满足1）式：

$\left\{\begin{array}{c}100\&le;S_C\&le;900\\ 100\&le;S_S\&le;900\\ |S_C-S_S|\&le;\frac{\mathrm{\&alpha;}}{2}(S_C+S_S)\\ |L_C-L_S|\&le;\frac{\mathrm{\&beta;}}{2}(L_C+L_S)\\ d\&le;\mathrm{\&gamma;}\sqrt{0.5(S_C+S_S)}\end{array}\right.;---1)$

若满足1）式，则将该云影对作为一个候选基准对，并在所有云影对判断完毕后进行步骤（9e），若没有一个云影对满足1）式，则进行步骤（9c），其中1）式中α为面积约束系数，β为周长约束系数，γ为间距约束系数，α＝0.3，β＝0.25，γ＝3；

（9c）将1）式中面积约束系数α、周长约束系数β和间距约束系数γ同时各增大百分之一，再重新判断是否有云影对满足1）式，如有云影对满足1）式则将其作为一个候选基准对，并在所有云影对判断完毕后进行步骤（9e）；否则进行步骤（9d）；

（9d）如果α<1、β＜1且γ<5，则转到步骤（9c）；否则，则该归一化后子图集

无基准对；

（9e）在所有得到的候选基准对中选择云影对质心连线长度d最小、阴影面积S_S最大的候选基准对作为归一化后子图集

的基准对；

（9f）分别对16个归一化后子图集

重复步骤（9a）至（9e），得到Ж个基准对，其中Ж≤16；

（10）统计上述Ж个基准对的质心连线主角度和质心连线主长度，得到最终基准对倾角A和最终基准对间距D；

（11）根据最终基准对倾角A和间距D对归一化后子图集

中的所有云块和阴影块进行配对，得到

的云影配对结果图P_ij；

（12）将归一化后子图集

的云影配对结果图P_ij与云影粗检测图CS_ij对应位置处的像素灰度值相减并求绝对值，得到

的未配对云影块图像；

（13）根据最终基准对倾角A和间距D对归一化后子图集

的未配对云影块图像中的云块补充检测对应阴影块，得到

的阴影补充检测结果图AD1_ij；

（14）根据最终基准对倾角A和间距D对归一化后子图集的未配对云影块图像中的阴影块补充检测对应云块，得到

的云补充检测结果图AD2_ij；

（15）将归一化后子图集

的云影配对结果图P_ij、云补充检测结果图AD1_ij和阴影补充检测结果图AD2_ij的对应像素相加，得到

的最终检测结果子图CaS_ij；

（16）将16个归一化后子图集

对应的16个最终检测结果子图CaS_ij按顺序拼接，得到一幅完整的最终检测结果图。

2.根据权利要求1所述的针对Landsat TM和ETM图像的厚云及其阴影检测方法，其中所述步骤（10）中统计上述Ж个基准对的质心连线主角度和质心连线主长度，按如下步骤进行：

（10a）设角度阈值Θ=20，并设基准对的序号为Д，其中1≤Д≤Ж，此处令Д=1；

（10b）选择第Д个基准对；

（10c）将选中的基准对记为Δ，将Δ的质心连线倾角分别与其余Ж-1个基准对的质心连线倾角作差并取绝对值，得到Ж-1个基准对的Ж-1个差异角；如果这Ж-1个差异角中有超过(Ж-1)/2个小于角度阈值Θ，则将所有小于角度阈值Θ的差异角对应基准对的质心连线倾角与Δ的质心连线倾角相加并求平均值，得到质心连线主角度，并将该质心连线主角度作为最终基准对倾角A，将所有小于角度阈值Θ的差异角对应基准对的质心连线长度与Δ的质心连线长度相加并求平均值，得到质心连线主长度，并将该质心连线主长度作为最终基准对间距D；否则执行步骤（10d）；

（10d）若Д<Ж，则将Д的值加1，并转到步骤（10b）；若Д=Ж，并仍未得到最终基准对倾角A和最终基准对间距D，则将角度阈值Θ的值增加5且令Д=1，并转到步骤（10b）。

3.根据权利要求1所述的针对Landsat TM和ETM图像的厚云及其阴影检测方法，其中步骤（11）所述的根据最终基准对倾角A和间距D对归一化后子图集

中的所有云块和阴影块进行配对，是根据云影对是否满足如下公式进行：

$\left\{\begin{array}{c}A-\mathrm{\&delta;}\&le;\mathrm{\&theta;}\&le;A+\mathrm{\&delta;}\\ D-\mathrm{\&xi;}\&le;d\&le;D+\mathrm{\&xi;}\\ S_S\&le;a{S}_C\\ L_S\&le;b{L}_C\end{array}\right.;---2)$

若在归一化后子图集的所有M×N个云影对中有任何一个满足2）式，则将该云影对在一幅与的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像的相应位置处保存，即将云对应像素的灰度值置为1，将阴影对应像素的灰度值置为0.5，并待所有云影对判断完毕后，将保存云影对的所有这些图像对应位置的像素灰度值相加，得到

的云影配对结果图P_ij；若没有任何一个云影对满足2）式或

的云影粗检测图CS_ij中包含的云块个数M=0或阴影块个数N=0，则将一幅与

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像作为的云影配对结果图P_ij，其中a为面积比系数，b为周长比系数，δ为角度松弛变量，ξ为间距松弛变量，a＝10，b＝5，δ＝20，ζ＝20。

4.根据权利要求1所述的针对Landsat TM和ETM图像的厚云及其阴影检测方法，其中步骤（13）所述的根据最终基准对倾角A和间距D对归一化后子图集的未配对云影块图像中的云块补充检测对应阴影块，通过如下步骤进行：

（13a）设M_X为归一化后子图集

的未配对云影块图像中云块的个数，其中M_X≥0，若M_X等于0，则将一幅与

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像作为

的阴影补充检测结果图AD1_ij；否则设m为

的未配对云影块图像中云块的序号，其中1≤m≤M_X，此处令m＝1；

（13b）根据第m个云块的质心像素坐标(X_cm，Y_cm)，利用3）式求出对应阴影块的质心像素坐标(X_sm，Y_sm)；

$\left\{\begin{array}{c}{X}_{\mathrm{sm}}=X_{\mathrm{cm}}-D*\cos \left(A\right)\\ Y_{\mathrm{sm}}=Y_{\mathrm{cm}}-D*\sin \left(A\right)\end{array}\right.;---3)$

（13c）在

的灰度平均图E_ij上选取以像素坐标(X_sm，Y_sm)为中心、K×K大小的窗，然后将该窗口内灰度值最小的像素点作为种子点，利用区域生长法在E_ij中进行区域生长，得到补充检测的阴影块；如果补充检测的阴影块其面积和周长分别超过对应云块面积和周长的

倍和ψ倍，则将该补充检测的阴影块剔除，否则保留，其中区域生长使用邻域判断阈值T，T=0.03，

ψ＝5，K＝9；

（13d）若m＜M_X，则将m的值加1，并转到步骤（13b）；若m＝M_X，则将得到的所有补充检测的阴影块与对应的云块在一幅与

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像的相应位置处保存，即将云对应像素的灰度值置为1，将阴影对应像素的灰度值置为0.5，得到归一化后子图集

的阴影补充检测结果图AD1_ij。

5.根据权利要求1所述的针对Landsat TM和ETM图像的厚云及其阴影检测方法，其中步骤（14）所述的根据最终基准对倾角A和间距D对归一化后子图集的未配对云影块图像中的阴影块补充检测对应云块，通过如下步骤进行：

（14a）设N_X为归一化后子图集

的未配对云影块图像中阴影块的个数，其中N_X≥0，若N_X等于0，则将一幅与

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像作为的云补充检测结果图AD2_ij；否则设n为

的未配对云影块图像中阴影块的序号，其中1≤n≤N_X，此处令n＝1；

（14b）根据第n个阴影块的质心像素坐标(X_sn，Y_sn)，利用4）式求出对应云块质心像素坐标(X_cn，Y_cn)；

$\left\{\begin{array}{c}{X}_{\mathrm{cn}}=X_{\mathrm{sn}}-D*\cos (A+180)\\ Y_{\mathrm{cn}}=Y_{\mathrm{sn}}-D*\sin (A+180)\end{array}\right.;---4)$

（14c）在

的灰度平均图E_ij上选取以像素坐标(X_cn，Y_cn)为中心、K×K大小的窗，然后将该窗口内灰度值最大的像素点作为种子点，利用区域生长法在E_ij中进行区域生长，得到补充检测的云块；如果补充检测的云块其面积和周长分别超过对应阴影块面积和周长的倍和ψ倍，则同时剔除该补充检测的云块和它对应的阴影块，否则都予以保留，其中区域生长使用邻域判断阈值T，T=0.03，

ψ＝5，K＝9；

（14d）若n＜N_X，则将n的值加1，并转到步骤（14b）；若n＝N_X，则将得到的所有补充检测的云块与对应的阴影块在一幅与

的云影粗检测图CS_ij尺寸相同、所有像素灰度值全为0的空白图像的相应位置处保存，即将云对应像素的灰度值置为1，将阴影对应像素的灰度值置为0.5，得到归一化后子图集

的云补充检测结果图AD2_ij。

Images