1、Haze Tool User Guide 使用方式. 2 云层厚度检测. 2 HOT13 2 操作步骤. 3 HOT123 4 操作步骤. 4 案例. 5 云层厚度完善. 8 Haze perfection TM 9 操作步骤. 9 Haze perfection QB. 11 Maskandinterpolation . 12 案例. 13 Thresholdandinterpolation 13 fill_sink . 14 flatten_peak. 15 操作步骤. 15 云层去除. 17 Dark Substract. 17 操作步骤. 20 Hist match . 22 操作步骤
2、. 22 Cloud Point . 22 操作步骤. 23 案例. 28 鸣谢. 29 引文. 29 使用方式 该模块是在ENVI4.4中二次开发而成。使用时只需将文件置于*ITTIDL64productsenvi44save_add文件夹内,然后运行ENVI,Basic Tool菜单中会出现一个Haze tool按钮。绝对不要修改文件名,否则无效。 Haze tool主要由三部分组成:云层厚度检测(haze detection),云层厚度完善(haze perfection),云层去除(haze removal)。这三部分前后连贯,每一步的结果都会影响到最终的除云效果,而且在参数的选择上主
3、观性较大,这也是这个模块的缺点即不够自动化。 有问题联系作者请加QQ:27126797 或者e-mail联系: 除云案例请见作者博客: 云层厚度检测 一副影像中往往云层厚薄不一,因此不同厚度云覆盖下的地表需要恢复的强度不一;当然,完全遮挡地面的厚云不在考虑范围内。在该模块中云层厚度检测有两个指数:HOT13和HOT123。两者都是相对厚度检测,而不是绝对的光学厚度检测,HOT123是作者在HOT13的基础上的改进。 HOT13 HOT13由加拿大遥感研究中心的Zhang ying提出的(Zhang et al. 2002) (原称HOT,13是作者加上去的以示与HOT123的区别)。根据地物在
4、蓝色(TM1)和红色(TM3)波段的高度相关性,在特征空间里绝大部分像素分布在晴空线上。云的存在会使得云下地物的光谱偏离这条晴空线,云越厚,偏离越大。HOT13等于偏移距离。 晴空线 : b1sin-b3cos-a=0 (1) HOT= b1sin-b3cos-a (2) 是晴空线的倾角, a晴空线的截距, b1,b3 分别是蓝色(TM1)和红色(TM3)波段。通过人工选择无云区域的1、3波段回归,得到晴空线。 操作步骤 1、首先打开需要处理的影像,建立无云区域的ROI作为对照,并保存(之后其他操作仍然要用到). 2、Basic Tool-haze tool-haze detection-HO
5、T13 3、弹出对话框要求选择需要处理的影像,并且在Spectral subset中选择蓝色和红色两个波段。 4、选择一个之前就建立并且已经打开的无云区域ROI,点击OK即在内存中生成一个HOT13图像。为了减少内存占用,生成的HOT13是放大十倍的int格式。 HOT123 在很多情况下, 蓝色(TM1)和红色(TM3)高度相关,相关系数大于0.9。但是,当地物更加复杂的时候往往就不成立。比如除了植被,同时存在大量的土壤和水体,水体非常浑浊等等。在Quickbird等高分辨率影像中,除柏油和水泥以外的其他彩色地物也会降低HOT13的准确性。因此作者对其进行了改进,提出了HOT123,即利用可
6、见光波段(TM123)提取云层厚度。 HOT123=k1*b1+k2*b2+k3*b3-b K1k2k3b这4个参数的值使得|Mean_cloud-Mean_clear|/SD_clear 最小。Mean_cloud和Mean_clear分别是有云区域(选择的区域的云越厚越好)和无云区域的HOT123平均值,SD_clear是无云区域的HOT123标准差。 HOT123要满足有云区域和无云区域的厚度差值尽可能大,而又要使无云区域(背景)的方差尽可能小。 操作步骤 主要步骤和HOT13基本一致,不过波段选择不是蓝色和红色两个波段,而是蓝色、绿色、红色三个波段,即TM图像上的第一第二第三波段。最后
7、一步的不同之处是,HOT123需要选择两个ROI:无云区域ROI和厚云区域ROI(其中的厚云是指将地物完全遮挡住的云层,但是如果图像中没有此类厚云,也可以选择相对最厚的薄云来代替)。无云区域的选择要尽量涵盖各种土地利用类型,厚云区域无此要求。 为了减少内存占用,生成的HOT123是放大十倍的int格式。 案例 左边是431组合的假彩色图,中间是HOT123,右边是HOT13。根据作者处理大量影像后的经验,推荐使用HOT123。当没有厚云存在,云非常薄的情况下,有的时候可以选择HOT13。 云层厚度完善 虽然HOT13HOT123尽量地突出云层信息,抑制背景信息,但是仍然有不少无云区域的云层厚度
8、值偏离零很大,主要发生在水、浪花、土壤、雪、建筑区等地物上。以HOT13举例说明,在如下的13波段特征空间里(横坐标是第一波段,纵坐标是第三波段),无云区域用灰色表示,厚云区域用黑色表示;无云区域形状如同飞机。一翼代表偏小的云层厚度检测,一翼代表偏大的云层厚度检测。云层厚度的完善实际上就是把无云区域的这辆飞机的两翼去掉。 因此,当地物很复杂的时候就需要对HOT值进行一些额外处理,以修正这些偏差。当然,云层厚度完善中的几种方法都是非必须的,可以根据影像的特征酌情选用,也可不用。 Haze perfection QB适用于几乎所有光学遥感影像,如Quickbird等4波段影像(当然也适用于TM等更
9、多波段影像),这类影像会有几个可见光波段被云污染。Haze perfection TM 只适用于TM等含有中红外波段的影像,中红外波段往往很少受薄云影响。不管选择哪一个,处理后的结果都需要重新将无云区域ROI内的HOT平均值归零;因此,必须确保在进行云层厚度完善操作之前内存中有无云区域ROI的存在,否则最后将不会进行归零操作。 我们将以以下案例举例说明。左边是有云TM影像,中间是无云的TM影像(城市区域),右边是有云的QUICKBIRD影像。可见有很多特别暗的背景的云层厚度是负值,也有很多特别亮而不是云的背景。 Haze perfection TM 这种方法的灵感来自Liang shunlin
10、的聚类匹配方法(Liang et al. 2001; Liang et al. 2002)。作者稍作修改用于HOT值的完善上。该方法以利用TM457波段(忽略这些波段的云污染)进行非监督分类的结果作为输入参数,通过计算每类地物的平均HOT值并将其减去以达到归零的目的。 操作步骤 1、检查457波段是否受云的影响,如果影响可以忽略,则利用TM457波段进行非监督分类。类的数量看地物复杂程度随意决定,推荐20-50类。 2、选择需要处理的HOT影像(必须有无云区域的ROI)。 3、选择对应的非监督分类图 4、选择无云区域ROI 5、决定是否要对某些方差特别大的类型进行插值操作。如果不需要,则选择C
11、ancel;如果需要,则在相对应的类型前面打勾,然后点击OK。结果会在内存中生成。 以下是案例的结果,这步操作的效果取决于地物类型的光谱特征。当存在某几种地物类型,其457波段相似,而123可见光波段相差很大的时候,这一步操作就会对这几种地物类型失效。 Haze perfection QB 点击Haze Perfection QB,会弹出对话框要求选择要处理的HOT图像。然后,弹出一个多项选择的对话框,选择合适处理方法。可以选择一种,也可以同时选择多种,作者不推荐同时选择多种方法处理。主要有四种方法可供选择:maskandinterpolation、thresholdandinterpolat
12、ion、fill_sink、flatten_peak。其中fill_sink和flatten_peak非常吃内存,是haze tool最大能处理的影像的限制步骤。 Maskandinterpolation 选择一个手动勾勒的需要掩膜并插值的异常区域ROI进行插值。该操作是人工方法,操作简单,当异常较少的时候比其他haze perfection方法要快。 运行完毕后,将跳出对话框要求选择无云区域ROI进行归零操作以纠正改变后的HOT带来的偏差。 案例 如下,左边有个多边形区域被手工勾勒并且插值,右边是原始图像。 Thresholdandinterpolation 通过设定一个阈值范围minimu
13、m, maximum,将范围之外的HOT值掩模并进行插值。插值方式采用周围像素取平均的方式。这种方法较简单,因为HOT值偏小的地物受到云的影响依然可以有较大的HOT值。因此,这种方法现在作者很少采用,仅仅作为一种异常值去除的方法(去除厚云,和一些极端负值)。 运算需要输入4个参数:minimum(最小值),maximum(最大值),kernel size for interpolation (插值窗口大小),kernel size for smooth (平滑窗口大小)。前三个参数用于插值,最后一个参数用于结果图像的平滑滤波。默认情况下,minimum=整幅影像最小值,maximum=整幅影像
14、最大值,kernel size for interpolation=7, 不进行平滑滤波。 运行完毕后,将跳出对话框要求选择无云区域ROI进行归零操作以纠正改变后的HOT带来的偏差。 fill_sink 该方法来自于水文领域中对地形图的处理(Planchon and Darboux 2002):填洼算法。我们把HOT图像看作地形图,HOT值偏小的地物往往会在地形图上形成一个一个的洼地。洼地虽然比周围的值要小,一个洼地可以在海拔很低的地方,也可以在一座高山上。这也是作者很少利用之前的maskandinterpolation的原因。 运行完毕后,将跳出对话框要求选择无云区域ROI进行归零操作以纠正
15、改变后的HOT带来的偏差。 案例 如下图,左边是经过fill_sink处理的HOT,右边是未经过处理的HOT。很明显,fill_sink对暗斑去除很有效。 flatten_peak 除了HOT值偏小的地物,还有HOT值偏大的地物。一个像素的HOT值很大,有可能是该像素被云覆盖,也有可能是本身HOT值偏大。作者利用云层厚度具有渐变的特征而地物边界骤变的特征将这两者区别开来。 操作步骤 1、预处理。在弹出对话框中选择morph reconstruction。并在接下去弹出的参数对话框中输入一个合适整数n, 进行n次形态学侵蚀操作。n 的大小要保证绝大多数偏大的地物能够被侵蚀完全。过大,则徒耗计算量
16、。 预处理步骤在内存中生成两个图像max_change_of_successive_erosion和 total_change_after_reconstruction。这两个图像作为输入图像用于下一步真正的flatten_peak操作。 2、在预处理完后再重新flatten_peak,不过这回在弹出的对话框中选择choose previous results。并在接下来的对话框中依次选择预处理中生成的max_change_of_successive_erosion和total_change_after_reconstruction。 3、对max_change_of_successive_e
17、rosion和total_change_after_reconstruction分别设定最大阈值maxchange_limit和totalchange_limit。将同时大于这两阈值的区域进行掩模并插值。Kernel size for mask interpolation和kernel size for smooth和之前提到的意义相同。 4、运行完毕后,将跳出对话框要求选择无云区域ROI进行归零操作以纠正改变后的HOT带来的偏差。 当然,haze perfection操作的目的尽可能地突出云的信息,抑制背景信息。这里提供的几种方法的最终结果并不一定完美,有待将来开发更多更好的算法。 云层去除
18、 完成了云层厚度检测和云层厚度完善之后即可进行云层去除。这里提供了三种方法:Dark Substract (暗物质法), Hist match(直方图匹配法), Cloud Point (云点法)。Cloud Point是在暗物质法的基础上的改进。以上的几种方法都是按照云层厚度对图像进行分层除云处理。 Dark Substract Dark substract即HOT的提出这Zhang ying采用的方法(Zhang et al. 2002),主要是对HOT图像分层分割,每个像素各个波段的DN值减去该像素所在层各个波段的最小值。 然后考虑到最小值不够稳定,因此作者采用Dal Moro推荐的固定
19、百分点方法代替最小值作为下限(Dal Moro and Halounova 2007)。 作者不推荐使用这种除云方法,因为云层不仅能够增强较暗地物的亮度,也能降低较亮地物的亮度,从而缩小云下地物的对比度。云的模糊效果不仅仅是加性的,还有乘性的效果。而且云的厚度应该是渐变的,人为将其按照厚度分割成一层层会导致结果的分块效应。 操作步骤 1、确保无云区域ROI(必须)以及需要处理的区域ROI(可选)已经在内存中。然后选择需要除云的影像,并且在Spectral subset中选择需要除云的波段。作者推荐处理可见光波段和近红外波段。 2、选择经过云层厚度检测和云层厚度完善(可选)而得到的HOT图像。
20、3、选择无云区域ROI。 4、选择需要进行除云处理的区域,如果需要对整幅影像进行处理则选择Cancel。 5、对图像云层进行切割分别处理。minimum,maximum范围之外将不被处理,默认情况下,minimum=最小HOT,maximum=最大HOT。Increment表示切割云层的厚度,默认情况下increment=1,太大则有明显的分块效应,太小则每层的像素数量太少而结果不够稳定。Lowerandupperlimit是作为下限的固定百分点,默认情况下等于2%,0%则表示最小值作为下限。 6、处理结果自动输出到内存,并附带一张各层下限的对应波段的最小值曲线。曲线波动越小效果越好。如下图,
21、说明上一步定义的minimum不合适,minumum至少要大于-4。 Hist match Hist match是分层直方图匹配方法,即将HOT图像分层分割,然后以无云区域ROI作为参照对每一层进行直方图匹配操作。这种方法在地物分布较均匀的情况下效果较好。当然,结果也会由于分层的原因导致分块效应。 操作步骤 与Dark substract的操作步骤雷同。 Cloud Point Cloud point是dark substract的改进,相同之处是仍然要对HOT图像进行分层分割,不同之处是每层同时提取上限和下限两个值。通过上限和下限分别回归得到两条回归直线的交点即云点(代表最厚的云),如下图(
22、蓝色波段为例)。任何一个像素点都将通过云点进行点投影到HOT=0的竖直直线上实现除云处理。 a是虚拟的云点,是各层的上限回归直线和下限回归直线的焦点。b 点是除云前的像素,c是除云后的像素,c点的HOT值等于0. 操作步骤 前两步与Dark substract的操作步骤雷同。 3、选择有云区域ROI用于计算云点,ROI的选择直接决定了最终的除云质量。该ROI的建立的目的在于得到稳定的云点,因此要求上限点和下限点稳定而且回归直线的回归系数尽可能的大。所以手工建立的有云区域要求:1、不同厚度的云覆盖,该区域内的HOT值范围尽可能大,使有更多的下限和上线点用于回归。2、区域内的地物尽量分布均匀,如此
23、可以减少地物造成的影响。如果选择Cancel,表明利用整幅影像计算云点,不仅计算量大而且效果不如人工选择更合适的区域。 4、对图像云层进行切割分别处理。minimum,maximum范围之外将不被处理,默认情况下,minimum=最小HOT,maximum=最大HOT。Increment表示切割云层的厚度,默认情况下increment=1,。Lowerandupperlimit是作为下限的固定百分点,默认情况下等于2%,那么上限的固定百分点是98%。0%则表示最小值作为下限,最大值作为上限。 确定后,会输出各层的结果图: 各层的像素个数 各层的下限值 各层的上限值 各层的上限值减去下限值 6、
24、根据上一步的结果选择用于计算云点的HOT范围。1、要求范围内各层的像素不能太少(根据各层的像素个数图),2、要求各层的上限值减去下限值随着HOT的增加逐渐稳定地减少(根据各层的上限值减去下限值图)。如上面的例子,选择范围最小值不能是-3以下,最大值不能超过15。 以下案例说明了第3步,有云区域ROI选择的重要性。下图的有云区域的ROI选择比较合适。 下图的有云区域ROI选择不合适。造成的错误是:蓝色、绿色、红色波段随着HOT的增加对比度增加,这个显然与现实不符。 案例 Cloud point的结果 Dark substract的结果 鸣谢 南京师范大学蒋圣同学在模块测试中的无私帮助。 引文 D
25、al Moro, G., & Halounova, L. (2007). Haze removal for high-resolution satellite data: A case study. International Journal of Remote Sensing, 28, 2187-2205 Liang, S.L., Fang, H.L., & Chen, M.Z. (2001). Atmospheric correction of landsat ETM+ land surface imagery - Part I: Methods. Ieee Transactions on
26、 Geoscience and Remote Sensing, 39, 2490-2498 Liang, S.L., Fang, H.L., Morisette, J.T., Chen, M.Z., Shuey, C.J., Walthall, C.L., & Daughtry, C.S.T. (2002). Atmospheric correction of landsat ETM plus land surface imagery - Part II: Validation and applications. Ieee Transactions on Geoscience and Remo
27、te Sensing, 40, 2736-2746 Planchon, O., & Darboux, F. (2002). A fast, simple and versatile algorithm to fill the depressions of digital elevation models. Catena, 46, 159-176 Zhang, Y., Guindon, B., & Cihlar, J. (2002). An image transform to characterize and compensate for spatial variations in thin cloud contamination of Landsat images. Remote Sensing of Environment, 82, 173-187