遥感影像阴影自动扩充提取算法
2019-03-29薛理杨树文马吉晶贾鑫闫如柳
薛理， 杨树文， 马吉晶， 贾鑫， 闫如柳
(1.兰州交通大学测绘与地理信息学院，兰州 730070；2.甘肃省地理国情监测工程实验室，兰州 730070)
0 引言
目前常用的高空间分辨率遥感影像阴影提取方法主要可分为2类： 一类是基于模型的[1-2]，使用该类方法需要传感器、太阳入射角和数字表面模型(digital surface model,DSM)等参数，同时需要一定的先验知识，因此只针对特定的影像，局限性较大； 另一类是基于阴影特征的[3-5]，该类方法通过分析阴影区与非阴影区的共性和差异性来提取阴影区。目前，第二类方法使用较为普遍，学者们已在该方向上提出了多种识别算法，如张云飞等[6]提出利用二维熵自动确定图像分割阈值； 虢建宏等[7]提出了一种多波段检测阴影的方法； 于东方等[8]提出利用势函数拟合灰度直方图进而获得阈值的阴影自动检测算法； 韩青松[9]提出了约束灰度范围的Otsu算法和基于散度差的一维Otsu算法及其快速递推算法，并将其拓展到二维空间，实现遥感图像阈值分割； 高贤君等[10]改进了Otsu阈值算法，可以自动获取各特征的合适阈值，实现阴影的自动检测； 李轶鲲等[11]提出基于多峰直方图的遥感影像阴影检测方法； 赵显富等[12]提出一种基于彩色模型的遥感影像阴影检测方法； 位明露等[13]结合HSV变换和区域生长原理，提出了一种改进的阴影检测方法； 杨猛等[14]提出一种结合特征分量构建和多尺度分割面向对象的阴影检测方法； 池毓锋[15]等利用Landsat8 OLI多光谱影像，结合地物对阳光反射的亮度曲线特征，优选波段组合，通过c！函数实现山区阴影提取。上述方法都取得了不错的效果，然而大部分方法通过分析整幅影像或部分影像，从而获得阈值，并没有考虑到每个阴影的特征，从而造成分割阈值稍大就会出现误提取，分割阈值稍小又会出现提取效果不理想等问题。为了避免误提，通常只能将大部分阴影区域提取出来，而靠近阴影边界的区域提取较为困难。由此作者通过分析各种影像近红外波段阴影边界的像素值特征，发现阴影局部的像素值变化率波动较小而边界的像素值变化率存在一个波峰，于是利用该特征对初提取的阴影进行扩充，提出一种遥感影像阴影自动扩充提取算法。
1 阴影边界特征分析
阴影是由于光源被遮挡而形成的，分为本影和半影。光线被完全遮挡的区域为本影，部分遮挡的区域为半影[3]，如图1所示。相关研究表明，当光源无穷远时半影区很小。通过分析发现遥感影像中多数阴影区和非阴影区边界较为明显，由于太阳距地球十分遥远，因此本文不考虑半影区。
图1 阴影分类Fig.1 Shadow classification
通过对各类遥感影像阴影边界进行统计分析，发现近红外波段中阴影边界的像素值变化较为明显。为了进一步扩大明显性，本文计算了阴影区到非阴影区像素值的变化率，即
(1)
式中：vi为第i个像素的变化率；li+1为第i+1个像素的值；li为第i个像素的值； ││代表取绝对值。
根据上述公式对QuickBird、资源三号(ZY3)、WorldView和高分一号(GF-1)影像阴影边界区域像素值变化率进行计算，累计求和并取平均值。由于阴影内部像素值变化不太明显，而阴影边界周围阴影区域的像素值变化较为明显，阴影边界处变化最剧烈。非阴影区域在阴影边界周围的像素值变化较为剧烈，这与阴影有关。但如果离阴影边界更远的非阴影区像素值变化剧烈，这可能是物体边界等现象引起的，与阴影无关。因此以阴影边界周围2个非阴影像素值为例，其像素值变化率见表1。
表1 阴影像素值变化率Tab.1 Change rate table of shadow pixel value
根据表1得出： ①从阴影区到非阴影区的过程中，像素值的变化先急剧增加后急剧减小，在阴影和非阴影边界处有一个波峰，在0.448左右； ②阴影区内部的像素值变化率很小，变化率在0.140左右； ③非阴影区的像素值变化率也在0.140左右，但都小于相应边界处波峰值。
通过以上分析可知： 根据近红外波段的像素值变化率在阴影边界处存在波峰和阴影区及非阴影区内变化率较小的特点，建立边界判断准则，可以对初步阴影提取区域进行扩充，提高提取效果。
2 算法描述与算法分析
根据上述原理，本文提出一种遥感影像阴影自动扩充提取算法，根据需要分为2种情况： 一种为半自动提取； 另一种为自动提取。
如果提取精度要求高，数据种类少，需要分析该类影像阴影内部区域近红外波段的像素值变化率，得出变化率的范围，进而得出阈值α。具体方法为： 随机选取20幅不同区域、不同时间的该类影像，通过人工计算每幅影像部分阴影区域内靠近阴影边界的近红外波段5个像素值变化率。统计所有计算得出的像素值变化率，并对其排序。为了避免特殊值的影响，去掉0.5%的较大值，将剩下像素值变化率的最大值设为阈值α。同时也可根据经验，人为设定阈值，不需要对像素值变化率进行统计。通过大量实验得出： 阈值比阴影区域内离阴影边界相差一个像素的像素值变化率(以下简称SPR)稍大较为合适，该方法适用于大部分影像。例如根据表1，ZY3影像的SPR为0.045，则阈值可以设为0.06。WorldView的SPR为0.04，则阈值可以设为0.05。GF-1的SPR为0.099，则阈值可以设为0.1。该设定阈值的方法适合于大部分影像，但对某些影像不适合，如QuickBird的SPR为0.13，如果根据上述获得阈值的方法，阈值会设为0.13～0.29之间的值，但通过实验发现该范围内的变化率并不适合，会出现大量非阴影区被误提取为阴影区的现象。对于该情况，可以采用第一种方法，通过统计和计算的方式获得阈值，也可以减小阈值，通过实验得出较为合适的阈值。本文通过实验得出QuickBird的阈值设为0.07较为合适。
通过对多幅GF-1，GF-2和ZY3等数据分析发现，同一类影像其阴影内部像素值变化率范围相差不大，因此得出阈值α可应用于同一类影像。
如果精度要求不高，数据种类繁多，数据量大，可以不用考虑阴影内部情况，在阴影初步提取的基础上，只需根据阴影边界的像素值变化率存在波峰的特征来扩充阴影初步提取结果。
2.1 阴影扩充
阴影扩充是根据阴影边界判断准则，沿着某一方向由阴影内部向外部进行扩充，直到遇到阴影边界为止。根据阴影提取的种类不同，阴影边界判断分为2种情况： ①如果为半自动提取，当像素值变化率小于前一个像素值变化率，且前一个像素值变化率大于阈值α时，则前一个像素值为阴影边界； ②如果为自动提取，当像素值变化率小于前一个像素值变化率时，则前一个像素值为阴影边界。阴影扩充算法流程如图2所示。
图2 阴影扩充流程Fig.2 Flow chart for shadow expansion
上述阴影扩充算法流程主要包括2个步骤： 首先，获得阴影扩充方向dir，同时获得该方向的阴影边界； 然后将其位置存储到栈stack中。该步骤分半自动提取和自动提取2种情况：
1)半自动提取。从stack中，取出栈顶数据，作为种子，向阴影dir方向扩充，计算像素值变化率vi和vi-1(vi-1为dir方向上vi的前一个像素的像素值变化率)。如果vi-1
2)自动提取。从stack中，取出栈顶数据，作为种子，向阴影dir方向扩充，计算像素值变化率vi和vi-1(vi-1为dir方向上vi的前一个像素的像素值变化率)。如果vi-1
2.2 阴影提取
2.2.1 半自动提取
首先，对影像做统计分析，获得阴影区域内近红外波段变化率阈值α； 然后，遍历阴影边界，根据阈值α和边界判断准则由阴影内部向外部扩充； 最后，遍历上一步阴影结果的左右边界，再向阴影左右方向进行扩充得到最终结果。半自动提取流程如图3所示。
图3 半自动提取流程Fig.3 Flow chart for semiautomatic extraction
具体步骤如下：
1)分析此类遥感影像阴影区域内近红外波段的像素值变化情况，确定该类阴影区域像素值变化率阈值α。
2)获得该类影像阴影初步提取结果，遍历阴影区域，获得阴影左、右、上、下边界位置，分别存放到栈stackL，stackR，stackT和stackB中，作为种子。由于获得阴影左、右、上、下边界位置的方法类似，因此以左、上边界为例。
左边界： 从图像左上角开始遍历阴影初步检测结果，如果Pi(第i个像素)为非阴影区域，Pi+1(第i+1个像素)为阴影区域，，则Pi+1为阴影的左边界，将i+1保存到stackL中，继续遍历像素，直到阴影初步检测结果遍历完为止。
上边界： 从图像左上角第一个像素开始遍历阴影初步检测结果，如果Pi(第i个像素)为非阴影区域，Pi+r(Pi+r代表整幅影像第i+r个像素，即Pi对应的下一行像素，r代表每一行的像素数量)为阴影区域，则Pi+r为阴影的上边界，将i+r保存到stackT中，继续遍历像素。由于遍历像素时，不能超出阴影初步检测结果的总像素数量m，即i+r≤m，i≤m-r，m-r为第n-1行(其中n为阴影初步检测结果的总行数)的最后一个像素。因此遍历像素，直到第n-1行最后一个像素为止。
3)按照2.1节中给出的阴影扩充算法，依次从左、右、上、下4个方向对阴影进行扩充。
4)遍历第3)步的阴影提取结果，获得阴影左、右边界，分别存放到栈stackL和stackR中，按照第3)步的方法依次向左、右2个方向扩充，算法结束。
对于某一类影像，其阴影内近红外波段像素值变化率类似，因此对于同一类影像，第1)步只需一次。通过实验得出，阴影扩充4次后，阴影到非阴影的过渡区域虽然大部分被提取出来，但仍然会有少部分未被提取。为了便于理解，以图的形式进行说明。如图4所示，图4(a)中I区为阴影初始提取区域，S区为实际阴影区域，I区被包含于S区。图4(b)为阴影初始提取区域左右扩充结果，SL区为左扩充区域，SR区为右扩充区域。图4(c)为上下扩充结果，ST区为上扩充区域，SB区为下扩充区域。通过图4(c)可以发现上下左右扩充之后，阴影区域4个顶角区域仍然没有被检测。图4(d)为在图4(c)的基础上获得阴影左右边界，并再次左右扩充后的结果，可以发现阴影4个顶角区域也被检测出。由于阴影区域大都为规则的多边形，因此扩充6次就可达到满意的效果，即左、右、上、下、左、右，再次扩充效果就不明显。
(a) 初始阴影(b) 左右扩充(c) 上下扩充(d) 左右扩充
图4扩充过程
Fig.4Expansionprocess
2.2.2 自动提取
自动提取和半自动提取过程类似，不同之处在于： ①自动提取无需分析遥感影像阴影区域内近红外波段的像素值变化情况，同时也无需获得阴影内部阈值α； ②自动提取步骤中过渡区域判断时无需考虑内部阈值α。因此本文不再展示自动提取流程，自动提取具体步骤如下：
1)获得该类影像阴影初步提取结果，遍历阴影初步提取结果，获得阴影左、右、上、下边界，分别存放到栈stackL，stackR，stackT和stackB中，作为种子。
2)按照2.1节中给出的算法，依次从左、右、上、下4个方向对阴影进行扩充。
3)遍历第2)步的阴影扩充结果，获得阴影左、右边界，分别存放到栈stackL和stackR中，按照第2)步的方法依次向左、右2个方向扩充，算法结束。
2.3 算法分析与验证
分别选取GF-1,QuickBird和ZY3影像数据进行实验，对近红外波段使用文献[11]的方法提取阴影，其结果作为阴影初步提取结果。本文不仅与文献[11]算法进行对比，而且也与人工经验阈值提取结果进行了对比。对比结果如图5、图6和图7所示，为便于观察阴影用黑色表示，背景用白色表示。图5(a)和图6(a)分别为GF-1和QuickBird城市高大建筑物影像，图7(a)为ZY3大面积水体、建筑物和部分植被的影像； 图5(b),6(b)和7(b)分别为文献[11]算法所得结果； 图5(c),6(c)和7(c)分别为本文自动提取结果； 图5(d),6(d)和7(d)分别为本文半自动提取结果，实验中GF-1，QuickBird和ZY3的阈值α分别为0.1，0.07和0.06； 图5(e),6(e)和7(e)分别为人工经验提取结果。
(a) 原图像(b) 阴影初步提取结果(c) 本文自动提取结果
(d) 本文半自动提取结果(e) 人工经验提取结果
图5高分一号数据阴影提取结果对比
Fig.5ComparisonofresultsofGF-1shadowextraction
(a) 原图像(b) 阴影初步提取结果(c) 本文自动提取结果
(d) 本文半自动提取结果(e) 人工经验提取结果
图6QuickBird数据阴影提取结果对比
Fig.6ComparisonofresultsofQuickBirdshadowextraction
(a) 原图像(b) 阴影初步提取结果(c) 本文自动提取结果
(d) 本文半自动提取结果(e) 人工经验提取结果
图7资源三号数据阴影提取结果对比
Fig.7ComparisonofresultsofZY3shadowextraction
由图5—7可知，图5(b)，6(b)和7(b)中阴影基本被提取出来，但阴影提取信息不太完整，过渡区很少被提取出来。图5(c)，6(c)和7(c)是本文自动提取算法的结果，阴影过渡区被提取了出来，与实际阴影区域更为吻合，并且阴影区域较为连贯。图5(d)，6(d)和7(d)为本文半自动提取算法的结果，其过渡区的提取比自动提取更加完整，阴影区域也更加连贯。图5(e)，6(e)和7(e)为人工经验阈值提取结果，使用逻辑与运算去除水体。可以看出阴影大部分被提取出来，但存在大量碎部和中空，且阴影未被完整提出。
对比本文算法结果和人工经验阈值提取结果得出： 虽然人工阈值提取能够大致提取出阴影，但存在大量碎部，且单个阴影区域提取不完整，而本文算法阴影提取较为完整。产生这种结果的原因为： 人工提取是根据整幅影像的直方图来确定阈值，从而无法考虑到每个阴影区域，造成阴影碎部较多，无法完整提取。而本文阴影过渡区提取算法能够遍历每个阴影区域，进而考虑阴影边界的像素值变化率特征对每个区域进行扩充，因此提取更完整。
对比本文自动提取结果和半自动提取结果得出： 半自动提取的阴影比自动提取更加完整，精度更高。产生这种结果的原因，可以以举例的方式进行说明，如图8所示。
图8 阴影边界像素Fig.8 Pixels of shadow border
图8中，A—F为某影像阴影边界周围的一段连续的像素，对应的像素值变化率分别为0.02，0.1，0.07，0.19，0.17和0.11。其中A—D为阴影区域像素，E和F为非阴影区域像素。假设阴影初步检测结果只检测到像素A，而像素B，C，D并没有被检测为阴影，即为漏检像素。如果按照自动提取从像素A开始对阴影进行扩充，但扩充到像素C时，由于0.07<0.1，则停止扩充，所以只能检测到像素B，像素C和D依然没有被扩充为阴影。如果按照半自动提取从像素A开始对阴影进行扩充，假设阈值α设为0.12。当扩充到像素C时，虽然0.07<0.1，但0.1<0.12，因此，将像素C标记为阴影，继续扩充。当扩充到像素E时，0.17<0.19，0.19>0.12，则停止扩充，此时像素B，C和D被重新检测为阴影。因此半自动提取精度要高于自动提取。
上述半自动提取方法中，当阈值α<0.1时，半自动提取的结果和自动提取类似。当阈值α≥0.19时会将非阴影区误检测为阴影区。因此阈值的设置需要按照本文方法进行确定，不能按照自己的意愿随意选取。
为了更加客观地分析结果，采用文献[12]提出的阴影检测准确率C，总体精度A和漏检率L指标对检测结果进行定量分析,即
(2)
(3)
(4)
式中:MFP为非阴影被误检测为阴影的像素数；MFN为阴影被检测为非阴影的像素数；MTP为正确检测的像素数。
根据上述阴影评价指标，本文将检测后的结果进行统计，结果如表2所示。
表2 阴影检测结果统计Tab.2 Statistical table of shadow test results
从表2中容易得出： ①初步提取结果、人工经验提取结果、自动提取结果和半自动提取结果的准确率依次增大，漏检率依次降低； ②初步提取结果总体精度最低，其准确率和总体精度大致相等，这表明初步提取结果非阴影被误检测为阴影的像素数很少，这与大部分阈值提取算法为了减少误提取而缩小阈值的做法相一致； ③自动提取结果准确率、总体精度和漏检率相对初步提取结果有明显改善； ④GF-1阴影初步提取结果的准确率和总体精度较低，同时本文算法提取结果的准确率和总体精度也较低。
通过上述分析得出，本文算法能够明显改善阴影初步提取结果，且提取结果的准确率优于人工经验提取结果，总体精度和人工经验提取结果相差不大，且漏检率低于人工经验提取结果。本文半自动提取结果的准确率高于自动提取结果，其总体精度和自动提取相差不大，其漏检率要低于自动提取结果。因此本文提取算法能够代替人工经验提取结果，更有效地实现了影像数据的批量处理和实际应用。
3 结论
本文提出了一种遥感影像阴影自动扩充提取算法，利用近红外波段阴影内部及边界特性，判断阴影边界，进而对阴影初步提取结果进行扩充，从而能够针对单个阴影区域，在一定程度上避免了只考虑全局或局部和避免误提而导致阴影提取不完整的现象。同时该算法也能够实现自动提取，便于对影像阴影进行批量提取，以减少人力，且思路清晰，易于理解。实验表明本文算法能明显改善阴影提取效果。
该算法的不足之处在于： ①当阴天拍摄时，影像较暗，导致阴影边界特征不明显，此时需要增加检测条件； ②若选择的像素值变化率阈值不合理时，会出现误提的现象； ③由于该算法是在阴影初步提取结果的基础上进行阴影提取，因此初步提取结果中若存在完全没有提取的阴影区域，本文算法依然无法提取出来。
