一种检测山地国产高分辨率卫星影像云与云阴影的方法

1.本发明涉及遥感卫星影像的计算机处理技术领域，具体而言，涉及一种检测山地国产高分辨率卫星影像云与云阴影的方法。


背景技术：

2.准确检测光学卫星影像中云及其阴影像元的位置是开展遥感影像应用的重要前提。然而，云的易变性、复杂性及其与雪的光谱相似性使得光学卫星影像中云及其阴影检测成为一大难点。此外，山地的地形起伏导致云阴影发生与山地阴影的混合和形变，进一步影响了山地遥感影像云及其阴影检测的难度。
3.近年来，我国国产资源环境卫星迅猛发展，海量高分辨率国产卫星数据的获取迫切需要高分辨率、自动化的云监测方法。然而，针对仅有可见光到近红外波段、缺少短波红外波段的大多数国产卫星ccd影像，如cbers-02b, cbers-04, hj-1a/b ccd1/2, gf-1/2/6等而言，现有的诸多成熟的云检测方法均难以适用，山地面积占我国国土面积的近70%，我国国产卫星影像的云及其阴影检测，必须需要考虑山地的特殊性。
4.然而，此类云监测算法极少考虑山地问题。为解决上述缺少短波红外和热红外波段国产高分辨率山地卫星影像云及其阴影检测问题，本发明提出了一种检测山地国产高分辨率卫星影像云与云阴影的方法。


技术实现要素：

5.本发明的目的在于解决上述问题，提供一种能够实现高效率、高精度，缺少短波红外波段和热红外波段的山地国产高分辨率卫星影像云与云阴影自动检测方法。
6.本发明的实施例通过以下技术方案实现：一种检测山地国产高分辨率卫星影像云与云阴影的方法，包括如下步骤：s1.基于像元亮度值对所采集的遥感图像进行云雪提取；s2.利用时间窗口内的多时相影像进行时相合成，消除云的干扰，得到无云参考影像；s3.基于所述无云参考影像以及像元的蓝波段反射率进行云雪区分；s4.计算表示目标影像纹理特征的区域协方差矩阵，基于区域协方差矩阵距离，得到云影像三维空间坐标；s5.基于所述云影像三维空间坐标计算云影像遮挡光线，得到云阴影位置。
7.根据一种优选实施方式，步骤s1具体包括：s11.对蓝波段、绿波段以及红波段的亮度均值进行归一化，选用亮度值较高的像元，进行厚云和积雪检测，表达式如下：
上式中，mean vis表示亮度均值归一化，band1表示蓝波段的反射率，band2表示绿波段的反射率，band3表示红波段的反射率，wt表示权重，band i{i=1,2,3}分别表示不同波段的反射率；s12.通过寻找投影方向突出薄云和雾霾的像元亮度值，并抑制无云区域的像元亮度值，进行薄云和雾霾检测，表达式如下：上式中，hottest表示雾霾优化转换。
8.根据一种优选实施方式，步骤s2具体包括：s21.利用时间窗口内的多时相影像进行时相合成，得到参考影像；s22.采用中值滤波算法对所述参考影像进行残云检测，表达式如下：上式中，表示第j个像元在合成日期x的蓝波段反射率值，median表示中值滤波算子，t1表示常数，表示第个像元在合成日期x的蓝波段反射率值；s23.采用时间序列滤波算法重构残云区域的像元，得到无云参考影像，表达式如下：上式中，表示第j个像元在合成日期x的反射率值，表示像元的原始数据，ci表示滤波系数，n表示滑动窗口内的像元个数并等于平滑窗口的大小2m+1。
9.根据一种优选实施方式，步骤s3具体包括：若一个像元的蓝波段反射率与无云参考影像相比满足如下表达式，则判定该像元为云，表达式如下：上式中，表示给定像元在日期d的蓝波段反射率，表示对应的无云参考影像相应的蓝波段反射率，表示无云参考影像到待检测影像之间的天数，表示待检测影像所处月份的实际天数，表示常数。
10.根据一种优选实施方式，步骤s4具体包括：s41.计算表示目标影像纹理特征的区域协方差矩阵，表达式如下：
上式中，r表示给定的区域，n表示区域内的像元个数，zk表示r内的d维特征点，表示特征点的均值，t表示转置；s42.计算区域协方差矩阵距离，表达式如下：上式中，表示参考影像协方差矩阵和检测的目标影像协方差矩阵的距离，表示协方差矩阵和的广义特征值，i=1k n；s43.计算sobel响应算子和laplace响应算子，得到八维灰度特征向量，表达式如下：上式中，表示卫星的四个波段的反射率，表示植被指数，表示影像蓝波段的sobel算子滤波器特征响应，表示影像蓝波段的laplace算子滤波器特征响应，表示灰度共生矩阵的局部均质性指数，表示像元的位置。
11.根据一种优选实施方式，步骤s5具体包括：基于云像元三维空间坐标计算云像元遮挡光线，表达式如下：上式中，表示云阴影的三维空间坐标，表示水平地表条件下的云阴影投影坐标，表示给定的云位置三维空间坐标。
12.根据一种优选实施方式，所述方法还包括，获取区域所有日期影像的参考云掩膜结果，对检测结果进行验证。
13.根据一种优选实施方式，所述参考云掩膜生成的步骤如下：采用面向对象的影像分类算法将影像分割为对象，进行云、雪和干净地表的区分；对区分结果中错分、漏分的像元进行改正，获取最终的参考云掩膜。
14.根据一种优选实施方式，所述采用面向对象的影像分类算法将影像分割为对象具体包括：将影像分割为规则瓦片，并采用随机函数获取每一期影像中的一个瓦片作为验证区域。
15.本发明实施例的技术方案至少具有如下优点和有益效果：（1）本发明所提供的一种检测山地国产高分辨率卫星影像云与云阴影的方法，利用影像的时相信息和影像自身空间纹理信息，来补充光谱特征空间，可提高云与地表光谱相似目标的可区分性；（2）本方法利用影像的时相信息和云的移动性特征，构建时相光谱差值规则，通过分析多时相影像的
光谱变化特征来描述光谱变化检测阈值，在时相分析基础上，增加信息维度，利用影像的纹理信息来表述影像中云、干净地表的空间特征，引入区域协方差矩阵技术，将影像光谱、空间纹理等多维信息进行综合，实现了多维信息的辅助云检测；（3）本方法考虑了复杂山地对云阴影的形变影响，解决了山地区域云阴影检测问题；（4）本方法对海量国产卫星历史数据以及在轨的卫星数据的云检测、卫星遥感数据应用等具有十分重要的意义。
附图说明
16.图1为本发明方法流程示意图；图2中第一行和第二行为月合成影像和待检测影像，第三行为光谱测试结果；图3中图3（a）为最大ndvi合成方法获得的无云参考影像，图3（b）为最小蓝波段合成方法获得的无云参考影像，图3（c）为本发明综合规则合成方法获得的无云参考影像；图4中图4（a）为三个典型区合成影像的光谱特征，图4（b）为中值滤波法的云识别结果，图4（c）为s-g滤波算法重构后的影像特征；图5 为图2待检测影像时相上下文信息监测结果；图6 为基于区域协方差矩阵的云检测流程；图7为云和积雪区域协方差矩阵距离直方图；图8 第一行为图2待检测影像基于区域协方差矩阵的云检测结果，第二行为对应的区域协方差矩阵距离；图9 为山地条件下太阳、传感器、云及其阴影的几何关系图；图10为不同山地地形云阴影投影的提取结果；图11 为贡嘎山地区全部影像云检测结果与参考掩膜结果的精度对比。
具体实施方式
17.为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
18.实施例1经申请人研究发现，云和雪在可见光波段的光谱特征相似，通常在该波段范围内难以有效区分。然而，雪在短波红外波段中反射率明显低于云，因此可以根据该波段的光谱特征来区别云和雪。
19.近年来，根据云的光谱特性和不同传感器的光谱特征，国内外已相继发展了多种云检测方法。在诸多云及其阴影的检测方法中，短波红外波段（1.5~2.2μm）和热红外波段（》7μm）是进行云雪区分以及云顶高度判别的重要谱段。通过设定相应规则，可以将云和非云区下垫面进行区分。
20.早期基于光谱特性的云检测方法以单一经验阈值法为主，然而由于阈值的确定带有一定的主观性，同时对先验知识要求也很高，而且在时间和区域上存在局限性，经验阈值法往往不具有较好的普适性。
21.近年来，国内利用云、冰雪等在更多波长范围的光谱特性，研究了多波段信息综合
的规则阈值法和光谱分类法。规则阈值法根据云在特定波段的光谱特性，通过设定各种规则和阈值特征来区分云和非云下垫面。规则阈值法中，常用的波段包括可见光波段、近红外波段、短波红外波段等。
22.卫星影像光谱分类法从模式识别的角度出发，通过云像元的特征提取及分类器的设计进行云检测。在分类器的选择上，常见的分类器包括支持向量机、神经网络、回归树等。
23.此外，云的纹理特征在云识别中也起着重要作用。云的纹理属于随机纹理特征。其纹理基元虽然变化较为复杂，但又区别于其它下垫面物体。
24.最后，近年来也发展了基于时相信息的云检测算法，但其仍然主要依赖短波红外波段的特性。
25.因此，本发明实施例提供一种检测山地国产高分辨率卫星影像云与云阴影的方法，参考图1所示，步骤具体如下：步骤一、基于云的光谱特性，通过光谱规则进行云和积雪的范围提取：首先，面向像元亮度值较高的厚云和积雪区域，基于亮度测试法进行厚云和积雪检测。需要说明的是，亮度测试法主要基于云和积雪在可见光波段的“亮-白”特性，对蓝波段、绿波段以及红波段的亮度均值进行归一化，选用亮度值较高的像元，进行厚云和积雪检测，表达式如下：上式中，mean vis表示亮度均值归一化，band1表示蓝波段的反射率，band2表示绿波段的反射率，band3表示红波段的反射率，wt表示权重，band i{i=1,2,3}分别表示不同波段的反射率。
26.对薄云像元而言，其既包含云的信息，也包含地表信息；因薄云像元亮度低于厚云像元，而裸岩、积雪等亮地表在遥感影像上也表现出与云相类似的光谱特征，因此无法通过亮度进行有效检测。因此，针对薄云像元，本发明实施例采用雾霾优化转换hottest进行薄云检测，即通过寻找投影方向以突出薄云和雾霾的像元亮度值，并抑制无云区域的像元亮度值，表达式如下：上式中，band1表示蓝波段的反射率，band3表示红波段的反射率。需要说明的是，雾霾优化转换hottest方法的基本假设是晴空条件下，大多地表像元反射率在可见光的不同波段相关性很高，而云以及薄云的出现将导致这种相关性发生较大变化。
27.步骤二、进行时相上下文信息测试分析：首先，获取无云参考影像。在本发明实施例中，采用时相合成法进行区域无云参考影像的获取，具体如下：基于云的位置多变特征，利用时间窗口内多期局部有云的多时相影像进行时相合成，消除云的干扰；在合成方法的选择方面，本发明实施例考虑不同的合成方法对植被特性和水体特性的适应性，提出一种考虑水体的综合合成规则，表达式如下：
上式中，上式中，ndvi
t
表示考虑水体的植被指数信息结果，表示考虑水体的蓝波段反射率值，water test表示水体特性，band4表示近红外波段的反射率，可以理解的是，该综合规则能够在水体和云的时间序列观测中提取出水体信息，而同时也能够保留最佳的植被观测。
28.进一步地，由于合成的参考影像中不可避免地会有残云存在，因此需要对合成的参考影像做进一步的检测。在本发明实施例的一种实施方式中，采用中值滤波算法对所述参考影像进行残云检测，表达式如下：上式中，表示第j个像元在合成日期x的蓝波段反射率值，median表示中值滤波算子，t1表示常数，以确保当整个时间序列中没有残云时，干净地表不被错误的剔除，表示第个像元在合成日期x的蓝波段反射率值。
29.需要说明的是，该算法基于云的运动性，假设云不会在同一个位置长期存在且合成的参考影像中的大部分云已被剔除，因此，残云的反射率会高于整个时间序列反射率的中值。
30.进一步地，在采用中值滤波算法获得合成影像中的残云检测结果后，采用s-g时间序列滤波算法重构残云区域的像元，得到无云参考影像，表达式如下：上式中，表示第j个像元在合成日期x的反射率值，表示像元的原始数据，ci表示滤波系数，n表示滑动窗口内的像元个数并等于平滑窗口的大小2m+1。
31.在该重构方法中，主要包括两个步骤，残云区域的光谱线性内插和基于s-g滤波算法的反射率重构。需要说明的是，s-g滤波算法是一个简单的多项式滤波器，它的设计思想是通过寻找合适的滤波系数以保护高阶距，即在对基础函数进行近似时，采用高阶多项式而非常数窗口，实现时间滑动窗口内的最小二乘拟合。
32.进一步地，利用云的位置易变性和积雪位置的固定性，采用时相上下文信息测试分析方法进行云和积雪的区分。基于所述无云参考影像以及像元的蓝波段反射率进行云雪区分，具体如下：时相检测模型：若一个像元的蓝波段反射率与无云参考影像相比满足如下表达式，则判定该像元为云，表达式如下：
上式中，表示给定像元在日期d的蓝波段反射率，表示对应的无云参考影像相应的蓝波段反射率，表示无云参考影像到待检测影像之间的天数，可以从合成影像中像元的合成日期与待检测日期之间的差值求取，表示待检测影像所处月份的实际天数，表示常数，以确保给定像元与无云参考影像蓝波段间的光谱差异阈值，可通过分析已知地物类型光谱反射率差值确定。
33.步骤三、空间上下文信息测试分析：针对步骤一和步骤二分析中未被检测出来的残留云像元，本发明实施例采用空间上下文信息测试方法进一步检测。计算表示目标影像纹理特征的区域协方差矩阵，基于区域协方差矩阵距离，得到云影像三维空间坐标，具体如下：采用区域协方差矩阵，构建区域纹理描述算子进一步地处理被错误检测和未能正确检测的云像元。需要说明的是，本方法利用影像的时相信息和影像自身空间纹理信息，来补充光谱特征空间，可提高云与地表光谱相似目标的可区分性。在本发明实施例中，为了计算影像的区域协方差矩阵，设i为多光谱的卫星遥感影像，f为从i中提取的w
×h×
d维的特征图像，表达式如下：上式中，函数表示卫星遥感影像反射率及其各种映射（例如梯度、倒数、滤波器响应等），表示影像图像坐标系的行列号坐标。在本发明实施例中，对于给定矩形区域，设为r内的d维特征点，可通过的协方差矩阵来表示区域r，表达式如下：上式中，r表示给定的区域，n表示区域内的像元个数，zk表示r内的d维特征点，表示特征点的均值，t表示转置；由于区域协方差矩阵的距离不属于欧式空间，因此在本发明实施例中，采用如下表达式表示两个区域协方差矩阵距离，表达式如下：上式中，表示参考影像协方差矩阵和检测的目标影像协方差矩阵的距离，表示协方差矩阵和的广义特征值，i=1k n。
34.本发明实施例针对仅具有可见光波段和近红外波段的大多数国产卫星影像，采用其四个波段的反射率数据作为光谱灰度特征，同时计算对应的植被指数作为光谱输入特征之一，为了引入影像的纹理特征，分别计算sobel响应算子和laplace响应算子，得到八维灰度特征向量，表达式如下：上式中，表示卫星的四个波段的反射率，表示植被指数，
表示影像蓝波段的sobel算子滤波器特征响应，表示影像蓝波段的laplace算子滤波器特征响应，表示灰度共生矩阵的局部均质性指数，表示像元的位置。
35.综上，本方法利用影像的时相信息和云的移动性特征，构建时相光谱差值规则，通过分析多时相影像的光谱变化特征来描述光谱变化检测阈值，在时相分析基础上，增加信息维度，利用影像的纹理信息来表述影像中云、干净地表的空间特征，引入区域协方差矩阵技术，将影像光谱、空间纹理等多维信息进行综合，实现了多维信息的辅助云检测。
36.步骤四、山地云阴影检测：需要说明的是，在云阴影的检测识别中，地形效应经常被忽略，相比平原地区，山区云投影到地面的距离以及投影形状受地形起伏的影响。特殊情况下，当云很低或观测方向垂直于云及其阴影时，云阴影有时还会被云所遮挡。本方法考虑了复杂山地对云阴影的形变影响，解决了山地区域云阴影检测问题。
37.本发明实施例在给定云位置三维空间坐标的情况下，计算云像元遮挡光线的表达式如下：上式中，表示云阴影的三维空间坐标，表示水平地表条件下的云阴影投影坐标，表示给定的云位置三维空间坐标。根据高程数据中各点的坐标、高程值带入上式即可计算获得与上式匹配的云阴影的位置。
38.步骤五、精度评价：本发明实施例采用监督分类和手工编辑相结合的方式获取区域所有日期影像的参考云掩膜结果，并在此基础上对检测结果进行验证。
39.参考云掩膜生成的步骤如下：采用面向对象的影像分类算法将影像分割为对象，然后进行云、雪和干净地表的区分；最后，采用手工编辑的方式对区分结果中错分、漏分的像元进行改正，获取最终的参考云掩膜。
40.为了降低手工编辑的工作量，在本发明实施例中，将影像分割为规则瓦片，并采用随机函数获取每一期影像中的一个瓦片作为验证区域。
41.以下结合附图和以中国国产环境减灾卫星（型号：hj-1a/b）影像为具体实施例对本发明作进一步具体描述：参考图1，图1示出了本发明总体技术路线流程。步骤一利用hj-1a/b影像可见光波段的亮度特征及统计特征，基于亮度测试法和雾霾优化转换将厚云、薄云和积雪等一通提取，参考图2第三行。
42.进一步地，开展无云参考影像合成，参考图3所示的合成影像效果。通过对比可以发现，本发明所提供的合成算法能够综合最大植被指数法和最小蓝波段法的优势，达到更好的融合效果。但融合影像中仍然后少量残云存在，因此继续利用中值滤波算法和s-g时间序列滤波算法进行时间序列的影像滤波和残云区域地表反射率的重构。从图4给出的三个典型区反射率重构先后的效果可以看出，中值滤波法能够成功识别区域中的云像元，重构
后的影像没有了残云的干扰，且影像的过渡性较好。
43.进一步地，进行时相上下文分析，基于云位置的多变特征，将待检测影像的蓝波段反射率与构建的无云参考影像蓝波段反射率进行相减；参考图5所示的时空上下文信息监测结果，由于云的反射率较高，因此无云地表差值后反射率接近于0，而云覆盖区域的差值反射率仍表现出较高的特征。
44.进一步地，区域协方差矩阵算法确定，区域协方差矩阵技术流程参见图6和图7。首先，分别采用sobel算子和laplace算子构建hj-1a/b影像的纹理特征，综合四个波段反射率、植被指数以及纹理特征，构建各待检测云影像和无云参考影像的特征波段。其次，设置一个3
×
3的像元矩阵区域位于待检测影像和无云参考影像。进一步地，分别计算两幅影像的区域协方差矩阵响应。然后计算待检测影像和无云参考影像的区域协方差矩阵距离，参考图8所示。
45.需要说明的是，区域协方差矩阵距离可根据分析已知影像中的云和积雪的值进行确定。当薄云出现时，区域协方差矩阵的增加非常显著。rcm分析后，根据云和积雪的区域协方差矩阵直方图高斯分布规律确定划分阈值，进行云雪区分。进一步地，在本发明实施例中，所有的云检测像元在8个联通方向向外扩展1个像元以消除云边缘像元的影响。
46.进一步地，进行云阴影检测。参见图9示出的山区阴影。若根据高程数据计算获得的交点位置唯一，则该交点即为云阴影的匹配点。
47.由于实际光线在投射过程中会被目标像元前方山体高程较高的坡面所遮挡，而计算获得的交点个数为多个则说明光线穿过所有的点，此时选择交点中高程值高同时云至阴影间光线段距离最短的交点作为阴影匹配点。
48.针对hj-1a/b没有热红外波段作为云顶高程判断的辅助，云高度无法直接计算获得，本发明实施例基于云高度变化通常在地表0.5 km至12 km的经验区间，首先计算获得投影的潜在范围；由于大气散射在短波波段要强于长波波段，且阴影区地表反射率主要是散射辐射的贡献，因此波长较长的近红外波段与波长较短的蓝波段的比值可以作为云阴影的光谱指数，结合潜在投影范围进行云阴影检测。
49.最后进行精度验证，参见图10所示，图10示出了参考云掩膜和本发明实施例所提供方法的云检测结果的云量散点图。参考图11，总体上，本方法在云检测精度上较为准确，r2达到0.9，出线性回归方程的截距较小，为3.29%，斜率为0.94。
50.本方法对海量国产卫星历史数据以及在轨的卫星数据的云检测、卫星遥感数据应用等具有十分重要的意义。
51.以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。