一种地表类型数据支持的遥感图像云阴影检测方法与流程

本发明属于检测技术领域，具体涉及一种地表类型数据支持的遥感图像云阴影检测方法。

背景技术

云阴影的检测相比较云检测而言具有更大的挑战。目前，主要的云阴影检测方法分为阈值法、分类法以及基于观测几何法。其中，阈值法云阴影检测和云检测一样，是应用最为广泛的云阴影检测方法。

经验阈值法的基本思想是基于光谱差异分析，使用一个或多个波段或者衍生度量进行云阴影检测的。shusp等人早期针对航空影像的阴影问题提出一种阴影区域识别和亮度补偿方法；simpsonjj等人提出一种针对avhrr的多光谱算法，该方法用几何和光学约束组合，从逐像素几何和图像光谱分析方面检测云阴影；martinuzzis等人研发出一种简单的半自动化方法检测并修复etm+数据中存在的云和云阴影；luoy等人开发出一种改进的modis云和阴影掩膜阈值法，将所需波段重采样成一致的250m分辨率，达到了较高的精度。国内很多学者也对云阴影的识别有较为深入的研究。陈奋等人针对高分辨率影像提出一种人机交互的半自动云阴影去除的算法，该算法以人机交互的半自动化方式将遥感影像分为云、阴影、清晰地物以及过渡地区四部分，之后通过建立阴影地区和清晰区域的直方图映射对阴影地区进行补偿，最后再对边界进行处理；米雪婷等人提出一种基于多时相遥感数据的云阴影检测方法，主要思想是使用同一区域相近时期的晴空影像作为支持，当待检测的图像反射率低于一定的阈值，则被判别为云阴影；云雅等提出一种针对gf-1的云和云阴影检测算法，对于云阴影，该方法利用阴影光谱特征和云对阴影地理位置关系结合的方法对阴影进行检测；孙林等提出一种针对landsat-8fmask算法改进的的云阴影检测方法，先利用云阴影光谱特征阈值对云阴影区域有一个界定，之后使用改进的云高计算方法基于云和阴影几何关系对云阴影进行识别。

技术实现要素：

针对现有技术中存在的上述技术问题，本发明提出了一种地表类型数据支持的遥感图像云阴影检测方法，设计合理，克服了现有技术的不足，具有良好的效果。

为了实现上述目的，本发明采用如下技术方案：

一种地表类型数据支持的遥感图像云阴影检测方法，包括如下步骤：

步骤1：地表晴空像元库与阴影像元库的构建；

在全球选取耕地、森林、草地、灌木、湿地、人造地表以及裸地样本，构建晴空像元库和云阴影像元库；同时还要选取一定的虚阴影样本，保证阈值模拟和概率计算时虚阴影识的正确率；

步骤2：最优波段的选取；

不同地物识别阴影的最佳波段是不同的，为了选取各地表的最优波段，根据建立的先验像元库，模拟各个阈值范围所对应的阴影像元正确率和误判率，若某波段满足较高阴影正确率的同时还有较低的误判率，则选取该波段作为检测波段之一；

步骤3：云阴影概率图的生成；

云阴影概率图生成是利用各地表类型先验地表库计算各个波长处的云阴影概率，将得到的结果以s函数拟合出云阴影概率计算公式；地物需要多个波段进行检测，分别以每个波段拟合的函数的标准误差作为权重加权合成最终云阴影概率结果，同时得到云阴影二值结果。

优选地，在步骤1中，构建地表晴空像元库与阴影像元库时，必须满足以下条件：

首先，像元库中的像元须保证其中的像元为阴影或晴空且二者数量接近，当阴影上方有薄云覆盖时，不宜选取为像元库像元，保证像元库的正确性；其次，像元库必须保证有足够多有代表性的样本；最后，样本的选取需要考虑阴影的类型。

优选地，在步骤2中，为了找到每种地物的最佳波段选择，根据建立的先验像元库，以反射率从0到1以0.01为间隔变化，随阈值的变化分别计算对应像元库里的阴影像元正确率和晴空像元误判率；阴影像元正确率为阴影像元正确个数与阴影像元总个数的比值；晴空像元误判率为误判为阴影的晴空像元个数与晴空像元总个数的比值；

当某一地表类型存在符合阴影像元正确率和晴空像元误判率均分别高于0.95、低于0.1条件的阈值，则将此波段列入该地表类型的最佳波段之一。

优选地，在步骤3中，具体包括如下步骤：

步骤3.1：阴影像元正确率和晴空像元误判率同时为0时阈值的最大值为min1，同时为1时阈值的最小值为max1，表达式如下，其中scr表示阴影正确率，fr表示晴空错判率，t表示阈值：

miν1＝tmax{scr＝0,fr＝0}(1)；

mαx1＝tmin{scr＝1,fr＝1}(2)；

当阈值t∈[min1，max1]时，在该范围内以0.001为间隔依次遍历n个区间，统计计算每个0.001范围区间内云阴影概率，取每个区间右端值为该区间的阈值，因此得到每一阈值处对应云阴影概率，第i个区间云阴影概率计算如下：

其中，nshadow表示落在[miν1+(i-1)×0.001,miν1+i×0.001]内阴影像元的个数，ntotal表

示该区间内像元总个数，由此，计算出区间内阴影占所有像元的比例，进而得到了某一阈值处对应云阴影的概率值；

步骤3.2：根据公式(5)，计算得到基于地表类型支持的云阴影概率；

在此计算中，记概率为1时对应最大的一个阈值为min，概率为0时对应的一个最小阈值为max，当影像像元值小于min时，云阴影的概率为1，像元值大于max时，云阴影概率为0；当阈值在min和max之间时，根据统计，各个阈值对应的云阴影概率会呈现一定趋势，据统计该曲线拟合效果最佳的函数为sigmoid函数，该函数由以下公式定义：

经拟合统计，基于各地表类型所选最优波段，分别进行以上云阴影概率的模拟统计，得到单波段s型曲线拟合系数、拟合优度以及标准误差；

总云阴影概率计算公式，如公式(6)所示；

其中，a、b和x0分别为s型函数的拟合系数，n为地表类型i所需波段数，wi,j为地表类型i的j波段的权重；ρi,j为遥感影像表观反射率，通过公式(5)计算出在地表类型i的j波段云阴影概率；当表观反射率低于min时，云阴影概率为1，认为该像元确定为阴影像元；当表观反射率高于max时，云阴影概率为0，则认为该像元为晴空像元；

步骤3.3：根据公式(6)，计算得到各最佳波段的加权合成，最终得到云阴影检测概率结果。

本发明所带来的有益技术效果：

针对一般阈值法云阴影检测阈值设置无地表参考进而造成检测精度低的问题，本发明基于像元波谱库特征信息与地表类型相结合，对landsat8影像进行云阴影检测，本发明利用地表类型支持解决了不同地表类型下云阴影光谱特征不同带来的问题，并以一定反射率值步长增长的方式解决了云阴影阈值设置困难的问题(云阴影和地物波谱特征差异界限模糊，不易确定阈值)，因此该方法具有较好的云阴影检测效果，达到较高的云阴影检测精度，带来十分有益的云阴影检测效果。

附图说明

图1为不同下垫面(地表)、实阴影和虚阴影下云阴影的表观反射率示意图；图(a)、(b)和(c)分别为植被、人造地表和裸地分别在晴空、虚阴影覆盖和实阴影覆盖时的表观反射率曲线示意图；

图2为阴影像元正确率和晴空像元误判率的模拟趋势示意图；

图3为各地表类型不同波段云阴影概率随反射率值变化趋势示意图；

图4为lscsd算法与目视解译云阴影占比统计结果回归分析图；

图5为多景landsat8影像数据分布示意图；

图6为植被类地表检测结果示意图；

其中，下垫面类型(a)-耕地；(b)-森林；(c)-草地；(d)-灌木；(e)-湿地；

图7为非植被类地表检测结果示意图；

其中，下垫面类型(a)-人造地表；(b)-裸地。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

云阴影是高空的云阻挡了来自太阳的辐射，使地面较少或无法获得应有的辐射。云阴影检测比云检测更具有挑战性。不同的是，云阴影的表观反射率与地表类型有很大的关系。图1统计了不同下垫面、实阴影和虚阴影下云阴影的表观反射率。(a)、(b)和(c)分别为植被、人造地表和裸地分别在晴空、虚阴影覆盖和实阴影覆盖时的表观反射率曲线，由图可知，不同地表类型上的阴影波谱曲线与地表类型相关性很大，其趋势分布一致且整体偏低。云检测是需要运用反射率低的波段进行检测，保证与云的高反射有足够的差异，同理，云阴影需要选择波段差异大的高反射率波段进行检测，能够达到较好检测的效果。通过模拟的方法，利用一定判定条件选取最优波段进行检测。

云与典型地表有较大的反射率差异，容易以一定的阈值去判定，云阴影的反射随着下垫面的变化而变化，难以设定统一的阈值判定。阈值的设定通常有两种方法：基于图像的阈值确定法和经验阈值确定法。由于不同图像地物属性不同、影像中反射率分布不同，此二类方法在空间和时间上有一定的局限性。针对云阴影光谱的不确定性，利用各地表类型划分不同阈值进行判别是十分有效的方法。以耕地、森林、草地、灌木、湿地、人造地表以及裸地为例，对基于地表类型数据支持的云阴影检测展开讨论。

本发明的云阴影算法是基于landsat8数据的云阴影和晴空地表的反射波谱差异，云阴影反射率和地表的反射率在有些波段差异相对较大，因此，未用传统阈值确定的方法确定阈值，而是基于不同地表类型的云阴影和晴空像元库，将landsat8所有波段考虑在内，以一定步长增长的阈值模拟每个波段不同阈值的晴空正确率和阴影正确率，进而以单波段反射率存在的差异确定云阴影检测的波段以及阈值，并生成云阴影概率图。用此方法获得的阈值更加准确合理，并且有很高的适用性。

1、landsat8阴影像元库建立

landsat8像元库由云阴影像元和晴空像元组成。构建像元库时，是在全球选取40景影像，对其中的耕地、森林、草地、灌木、湿地、人造地表以及裸地进行人工目视采样。首先，像元库中的像元须保证其中的像元为阴影或晴空且二者数量接近，当阴影上方有薄云覆盖时，不宜选取为像元库像元，保证像元库的正确性；其次，像元库必须保证有足够多有代表性的样本，选取40景影像进行采样，保证了样本数量，同时在全球范围内选取，又保证了样本具有在不同区域的代表性，减小空间局限性，使样本更丰富全面；最后，样本的选取需要考虑阴影的类型，像元库的建立是为了更加全面的统计不同地表类型中阴影的特征和计算阈值，所以虚阴影也需要考虑在内，虚阴影既含有地表信息，又存在阴影特征，为了使模拟和阈值更加准确，在样本采样时，同样需要有虚阴影的像元，使阴影像元库包含的阴影特征更加完整，进而为之后的阴影阈值模拟计算提供精度保障。

2、最优波段选取

由于不同地表类型的云阴影反射率趋势与地表本身反射率相关，并且云阴影自身反射率较低且不恒定，无法通过波谱曲线找到合适的阈值去界定云阴影的范围。综合考虑了不同地物和云阴影波谱在可见光、近红外和短波红外的微小差异，以单波段检测方法进行云阴影检测的波段选取。

为了找到每种地物的最佳波段选择，根据建立的先验像元库，以反射率从0到1以0.01为间隔变化，随阈值的变化分别计算对应像元库里的阴影像元正确率(阴影像元正确个数与阴影像元总个数的比值)和晴空像元误判率(误判为阴影的晴空像元个数与晴空像元总个数的比值)。图2为阴影像元正确率和晴空像元误判率的模拟趋势示意图，当阈值很低时，阴影像元正确率和晴空像元误判率都为0，之后阈值开始识别出一些阴影，阴影像元正确率有很明显的提高，并且对于晴空像元的判别也较为准确，有着较低的晴空像元误判率，直到阴影像元正确率和晴空像元误判率的差异达到最大，随着阈值增大，阴影像元正确率趋近于1，而开始有更多的高反射率晴空被误判为阴影，晴空像元误判率提高并趋近于1。图2中，当某一地表类型存在符合阴影像元正确率和晴空像元误判率均分别高于0.95、低于0.1条件的阈值，则将此波段列入该地表类型的最佳波段之一。

3、基于地表类型支持云阴影概率图生成方法

通过波段选择后，选出每个地表类型识别阴影的最佳波段，记图2中阴影像元正确率和晴空像元误判率同时为0时阈值的最大值为min1，同时为1时阈值的最小值为max1，表达式如下，其中scr(shadowcorrectrate)表示阴影正确率，fr(faultrate)表示晴空错判率，t表示阈值：

miν1＝tmax{scr＝0,fr＝0}(1)；

mαx1＝tmin{scr＝1,fr＝1}(2)；

当阈值t∈[min1，max1]时，在该范围内以0.001为间隔依次遍历n个区间，统计计算每个0.001范围区间内云阴影概率，取每个区间右端值为代表该区间的阈值，因此得到每一阈值处对应云阴影概率，第i个区间云阴影概率计算如下：

其中，nshadow表示落在[miν1+(i-1)×0.001,miν1+i×0.001]内阴影像元的个数，ntotal表示该区间内像元总个数，由此，计算出区间内阴影占所有像元的比例，进而得到了某一阈值处对应云阴影的概率值。在此计算中，记概率为1时对应最大的一个阈值为min，概率为0时对应的一个最小阈值为max，当影像像元值小于min时，云阴影的概率为1，像元值大于max时，云阴影概率为0。当阈值在min和max之间时，根据统计，各个阈值对应的云阴影概率会呈现一定趋势，据统计该曲线拟合效果最佳的函数为sigmoid函数(s型曲线)，该函数由以下公式定义：

经拟合统计，基于各地表类型所选最优波段，分别进行以上云阴影概率的模拟统计，得到单波段s型曲线拟合系数、拟合优度(r²)以及标准误差(standarderrorofestimate，see)如图3和表1所示，各类地表的每一波段的拟合函数优度较高，与样本数据有很强的相关性，且标准误差很小。

表1各地表类型最优波段拟合方程系数、拟合优度以及标准误差

根据公式(5)，计算得到基于地表类型支持的云阴影概率；

总云阴影概率计算公式，如公式(6)所示；

其中，a、b和x0分别为s型函数的拟合系数(表1)，n为地表类型i所需波段数，wi,j为地表类型i的j波段的权重；ρi,j为遥感影像表观反射率，通过公式(5)计算出在地表类型i的j波段云阴影概率；当表观反射率低于min时，云阴影概率为1，认为该像元确定为阴影像元；当表观反射率高于max时，云阴影概率为0，则认为该像元为晴空像元。

不同的波段和不同的s曲线得到云阴影概率的结果是不同的，如单纯将各波段取并集，则会形成阴影误判误差的累积，导致概率结果不准确；如直接取交集，则会损失概率信息，使阴影漏判率增大。

根据公式(6)，计算得到各最佳波段的加权合成，最终得到云阴影检测概率结果。

因此，综合上述问题，根据拟合s函数的标准误差，对所涉及的波段赋予不同的权值，误差大的权值低，误差小的权值高，这样设置使概率结果更加准确可信。表2为各地表每一波段对应的权值，公式(5)和(6)为最终云概率计算公式，根据云阴影概率图很容易得到云阴影结果，以75％为界限，大于该概率的像元则被判定为云阴影像元，否则为晴空像元。

表2各波段云阴影概率min、max及权重分布

需要指出的是，因水体和海洋反射率低，其阴影目视判别困难，样本选取无法保证准确度，因此在这两种地表类型上空采用恒定阈值进行阴影判别。由波谱曲线可知，水体的反射率呈现随波长增大而减小的趋势，在深蓝波段反射率最大(大于0.1)，对阴影降低反射率影响的敏感性大，故当该波段反射率低于0.1时，被判定为含有阴影的像元。

此外，植被区域会出现植被类型变为水体类型、地表库和图像上的植被与水体存在一定边界误差以及将水体误判为阴影等的误差可能性。因此在得到阴影检测结果后，在植被类地表(耕地、森林、草地、灌木及湿地)的检测结果中，以ndvi小于0为判别条件去除水体误判为阴影的情况。

云阴影检测结果的定量验证仍采用之前云检测结果验证的方法，即对各类型在对应影像中随机选取6个500*500像素大小的样本区域进行目视解译云阴影区域工作，当云和云阴影同时存在时，将该类像元划分为云像元。之后将目视解译完成的结果作为真值，与lscsd(landcovercloudshadowdetection，地表类型云阴影检测)算法云检测结果进行比对。利用云阴影占比(csp，cloudshadowproportion)来定量评价lscsd算法精度。

首先对样本区域目视解译结果云阴影占比与lscsd算法结果云阴影占比进行回归分析，并计算出lscsd云阴影检测结果的rmse(rootmeansquareerror，均方根差)；之后利用crcs(correctrateofcloudshadow，云阴影正确率)和crs(correctrateofclear-sky，晴空正确率)两个指数对云检测结果进行云阴影像元和晴空像元详细的统计分析和精度评价；最后，统计所有样本点的总体正确率，即lscsd算法对云和晴空像元均正确识别的正确率(tcr)。

图4为lscsd算法与目视解译云阴影占比统计结果回归分析图，表3为各地表类型云阴影占比的rmse计算结果。可以看出各地表类型在与真值比较时，lscsd算法云阴影占比的总体趋势与目视解译结果保持一致，其总体rmse为3.37％，达到较高相关性。分析其误差，精度总体低估，但高估现象较云阴影占比明显，说明阴影的高估比云要多，这是由于阴影与地表呈现晴空时的波谱特性更具有相似之处。由图4和表3看出，湿地偏离参考线较多，低估现象较其他地表类型明显，其rmse为5.83％；其次误差较大的地表类型为耕地和森林，rmse分别为4.62％和4.54％；云阴影占比统计中误差最小的地表类型为人造地表，rmse为1.42％。

表3各地表类型云阴影占比rmse结果

该统计表明湿地云阴影占比统计误差最大，由于湿地存在各种程度的暗地表，在选取晴空像元样本时，包含了一些较暗地表的晴空，同时在选取云阴影像元时也选取了一些稍亮地表的虚阴影，因此相对于其他地表，湿地晴空像元库和阴影像元库的差别较小，因此lscsd算法在湿地的云阴影占比误差相对比其他地物大。人造地表由于本身反射率较高，与云阴影的光谱信息区分较大，因此云阴影占比统计最接近于真实值。总的来说，所有地表的云阴影占比误差普遍较低，精度较好。

对于landsatoli和tirs(thermalinfraredscanner，热红外扫描仪)数据，需要经公式(7～9)将影像本身的dn(digitalnumber，像元亮度值)值转化为具有明确物理意义表观辐亮度，进而计算对应表观反射率或亮温，之后再对待检测图像进行进一步处理。转化公式为：

lλ＝gain·qcal+bias(7)

其中，gain为增益值，bias为偏移值，qcal表示一级产品的dn值，lλ表示传感器入瞳出的光谱辐亮度，单位为w·m^-2·sr^-1·μm^-1；ρ表示为大气顶层表观反射率，d为日地平均距离，esunλ为大气顶层的太阳光谱辐照度，θs为太阳天顶角；t为星上亮温，单位为k，k1和k2分别为定标常量，对于landsat8tirs10来说，其值分别为774.885w·m^-2·sr^-1·μm^-1和1321.08k。

在精度验证中，在选取各类典型地表区域以及复杂地表区域在全球分布的多景landsat8影像数据分布如图5所示，由于北美洲和亚洲面积较大且地表类型相对多样，因此在该区域选取图像偏多，每一类型选取影像尽量保证均匀分布。进行上述预处理过程以及云检测后，对实验结果进行精度验证分析。共选取60个样本区域，以样本区域的假彩色合成图为基准勾画云图斑，对lscsd云检测结果进行定量验证。

水体本身反射率极低，无论是目视判别对比还是勾画阴影区域进行定量分析，均无法得到较为准确的判断，影响精度评价，因此，云阴影的精度验证暂且不考虑下垫面为水体类时的情况。云阴影的波谱特征受下垫面自身波谱特征的影响，因此在分析时直接将地表分为植被类地表与非植被类地表。图6和图7分别为植被类与非植被类云阴影检测结果目视解译对比图，其中，图中左侧为假彩色合成图像，中间为云阴影概率图(由黑色到白色，云阴影概率由0到1)，右侧为最终云阴影检测二值结果。为了更加直观、完整地表现出云阴影检测的结果，选取区域时，综合考虑实阴影、虚阴影，同时也就不同形态的阴影进行目视解译的对比。

植被类地表

图6为植被类地表检测结果示意图，地表类型包括耕地、森林、草地、灌木地及湿地。耕地情况较为复杂，地表类型由植被向裸地过度程度不一，同时，草地和灌木地也存在稀疏程度和季节变化导致的不同地域地表特性的差异。lscsd算法在图(a)裸地和植被覆盖不均一，以及图(d)中呈现较多裸地信息的情况下，均可以完整地识别出云阴影，有较强的空间连续性和适应性。这是由于在选取样本时，不仅考虑了植被茂盛时期，同时也就稀疏植被区或呈现出裸地特性的区域也考虑在内，得到的最佳波段中均含有近红外和短波红外，尤其是在近红外波段，植被和裸地均为高反射，和阴影有很大的区分性，因此能最大程度地降低地表类型不一引起的判别误差。

图(b)为裸地上阴影的识别结果。裸地一般反射率偏高，因此也是较容易判别阴影的地物类型之一。可以看出对于薄云产生的密集碎状阴影，包含实阴影和虚阴影，绝大部分的是可以被识别的，但在少量云和阴影混合的区域，受到高反射率云的影响。总体看来，该算法在裸地上空有较高的检测精度。

综合植被类地表和非植被类地表的目视解译情况，无论在不同植被覆盖度的植被类地表上，还是在高亮的非植被类地表上，lscsd算法在各类地表上均能表现出较好的识别能力，尤其是对各种形态虚阴影的识别。但也在某些地表类型区域自身呈现特殊的较低反射率时，lscsd算法可能会将其误判为阴影。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。