CN105139375A - 一种结合全球dem与立体视觉的卫星影像云检测方法 - Google Patents

Abstract

本发明公开了一种结合全球DEM与立体视觉的卫星影像云检测方法，利用密集匹配的DSM数据和现有的全球DEM数据SRTM，结合高程与影像灰度信息提取云区域。主要步骤如下：通过多影像多匹配基元方法获得模型匹配的M-DSM，基于特征线的方法配准M-DSM与SRTM，通过高程差值获得初始种子点，然后把种子点计算到像方，结合影像的灰度信息提取出云区域。本发明在有大面积雪的影像，能有效的避免雪的干扰，准确地检测出云区，提高了遥感影像的利用率。

Description

一种结合全球DEM与立体视觉的卫星影像云检测方法

技术领域

本发明属于遥感技术领域，涉及一种卫星影像云检测方法，尤其涉及一种结合全球DEM与立体视觉的卫星影像云检测方法。

背景技术

卫星影像中云的存在，导致影像的可用信息缺失，给地形的自动提取带来了很多错误。同时，云覆盖给影像处理和遥感影像应用带来很多困难。早期的云检测主要靠人工做出判读，这种方法不仅费时，而且对云量的判断带有很大的主观性。因此光学遥感影像中高效的云区域自动识别和检测方法成为光学遥感影像处理和应用中的研究热点。

传统的云检测分为利用单一时相的影像数据进行判定和利用多时相的影像数据进行判定两大类。在利用单时相影像数据提取云区域的阈值分割算法和分类算法中，检测结果与选用的阈值有很大关系。阈值选择的困难使得这些算法很难识别和检测出影像中的薄云，且无法排除雪域和雾霾的影响，因此对于复杂的地物场景会出现很多错误分类；结合云阴影进行云检测的方法则会受到耕地、灌溉、裸露地的影响；利用多时相的影像数据则会因为时间变化造成地物变化，或者不同时相的影像色调引起差异所引起的错误分类。

发明内容

针对传统监督变化检测方法的不足，本发明提出一种结合全球DEM与立体视觉的卫星影像云检测方法。

本发明所采用的技术方案是：一种结合全球DEM与立体视觉的卫星影像云检测方法，其特征在于，包括以下步骤：

步骤1：将遥感卫星获得的卫星影像(如天绘立体测绘卫星影像)密集匹配获得的M-DSM(MatchingDigitalSurfaceModel)；

步骤2：将获得的M-DSM与现有的全球DEM数据SRTM(ShuttleRadarTopographyMission)通过地形配准后进行高程比较得到DEM的差值图；

步骤3：在DEM差值图中去除误差大于200米的点，然后剩下的空间点统计出单位权中误差(RMSE)，将大于三倍单位权中误差的点通过卫星影像三线阵正视影像的RPC参数投影到卫星影像上，得到疑似云区域种子点的像片坐标；

步骤4：根据影像灰度信息对疑似云区域的种子点进行区域增长和分类；

步骤5：去除增长区域中种子点灰度值小于100的低亮度区域、纹理缺乏区域和雪区域，得到最终的云区域。

作为优选，步骤1中所述的密集匹配方法采用的是多影像多匹配基元法(MPM)，采用由粗到细(coarsetofine)的多级影像匹配策略，基于多影像的匹配算法，并结合多种匹配基元的匹配结果，所述的多种匹配基元包括特征点、格网点和特征线。

作为优选，步骤2的具体实现包括以下子步骤：

步骤2.1：分别在全球DEM数据SRTM与M-DSM上提取地形特征线；

步骤2.2：通过投票算法实现两者的地形特征线的整体套合；

步骤2.3：进行特征点局部的最小二乘相关匹配提高匹配的精度；

步骤2.4：使用地形匹配获得的同名点解算M-DSM与全球DEM数据SRTM之间的空间相似变换关系式，将全球DEM数据SRTM相对于M-DSM进行校准；

步骤2.5：校准后进行高程比较得到DEM的差值图。

作为优选，步骤2.2中所述的投票算法，投票是根据地理特征线上除开地形特征线末梢点外的每个点，计算与其在同一地形特征线上的下一个相邻点之间的连线与水平方向的夹角，设定角度阈值α，α≤5°，以M-DSM的地形特征线中的点所对应的角度为依据，在全球DEM数据SRTM的地形特征线中搜索与其角度差小于阈值的匹配点作为备选匹配点。例如，对于M-DSM的地形特征线上的点M，其对应了一个角度值α_M，计算其与全球DEM数据SRTM的地形特征线中的点(如点N)所对应的角度值α_N的差的绝对值，若|α_M-α_N|小于阈值，则将点N选为点M的匹配备选点。根据M与N点的坐标差值来投票，对于所有与点M相匹配的点，都按上述方法在其对应的投票区投上一票，即完成点M的投票。对于单个点来说，正确的投票只有一票，而错误的投票可能有多票，正确的投票信息会被错误的投票信息淹没。然而将所有点的投票信息汇总之后，正确的投票会积累起来，而错误的投票往往是随机的，总体上呈现出均匀分散的状态，正确的投票会凸显出来，从而得到正确的结果。最后得到M-DSM与全球DEM数据SRTM的整体坐标偏移。

作为优选，步骤2.4中所述使用地形匹配获得的同名点解算M-DSM与参考DEM之间的空间相似变换关系式；空间相似变换关系为式：

$\left[\begin{array}{c}{X}_{tp}\\ Y_{tp}\\ Z_{tp}\end{array}\right]=\mathrm{\λ}\·\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]\·\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]+\left[\begin{array}{c}\mathrm{\Δ}X\\ \mathrm{\Δ}Y\\ \mathrm{\Δ}Z\end{array}\right];$

其中假设M-DSM上某点的坐标为(X，Y，Z)，该点在参考DEM上对应的同名点的坐标为(X_tp,Y_tp,Z_tp)，λ是缩放系数，R是由角元素Φ，Ω，Κ的函数组成的方向余弦，a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为R的9个元素，(ΔX,ΔY,ΔZ)是坐标平移量；空间相似变换的求解采用间接平差法，求解参数的同时探测并剔除匹配结果中的粗差，最后平差解算得到的中误差用来代表地形匹配的精度。

作为优选，步骤3中所述的疑似云区域种子点是比较遥感影像密集匹配生成的M-DSM和步骤2.4校准后的全球DEM数据SRTM，提取其中高程差大于三倍单位权中误差(RMSE)的格网点，通过三线阵正视影像的RPC参数将这些点投影到影像上，得到种子点的像片坐标。

作为优选，步骤4中所述的区域增长的具体过程为：以种子点像素作为生长的起点，然后将种子点像素周围邻域中与种子像素有相同或相似性质的像素合并到种子像素所在的区域中，将这些新像素当作新的种子像素继续进行上面的操作，直到再没有满足条件像素可以包括进来；其中所述的相同或相似性质的像素是根据事先确定的生长或相似性准则来判定，所述的区域增长的合并准则是每次合并时的灰度差值是和初始的种子点像素进行比较小于设置的阈值(一般设置为10或者15)。

作为优选，步骤5中所述的云区域和雪区域的区分主要依据区域中种子点的个数与区域总像素个数的比，小于预设的阈值为(一般设置为0.1)时则认为区域不是云区域。

本发明相对现有技术，其有益效果为：云提取正确率高，在无大面积的雪干扰的情况下也可以达到95％的正确率；在有大面积雪的影像，能有效的避免雪的干扰。最终能获得有效剔除云等因素干扰的DSM，提高了DSM的质量和利用率。

附图说明

图1：是本发明实施例的总体流程图。

图2：是本发明实施例的多影像多基元法流程图。

图3：是本发明实施例的无云影像M-DSM。

图4：是本发明实施例的地形配准算法流程图。

图5：是本发明实施例的地形特征线的提取示例图，其中(a)为山谷线影像，(b)为山脊线影像。

图6：是本发明实施例的云区域自动识别流程图。

图7：是本发明实施例的典型无云影像(a)和有云影像(b)中种子点区域集合的灰度直方图示例。

图8：是本发明实施例的少量云分布的自动提取效果图，其中(a)为原始影像，(b)为自动匹配的M-DSM影像，(c)为提取的云区域，(d)为使用SRTM填补后的DSM影像。

图9：是本发明大面积云分布的自动提取效果图，其中(a)为原始影像，(b)为自动匹配的M-DSM影像，(c)为提取的云区域，(d)为使用SRTM填补后的DSM影像。

图10：是本发明有雪分布时云区域的自动提取效果图，其中(a)为原始影像，(b)为自动匹配的M-DSM影像，(c)为提取的云区域，(d)为使用SRTM填补后的DSM影像。

图11：是本发明云和雪混叠时云区域的在自动提取效果图，其中(a)为原始影像，(b)为自动匹配的M-DSM影像，(c)为提取的云区域，(d)为使用SRTM填补后的DSM影像。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

实施例用天绘卫星三线阵全色影像进行试验。天绘卫星是立体测绘卫星，具有前视、下视、后视的全色相机，地面分辨率为5m，同时具有2.0m的全色高分相机，10m分辨率的多光谱相机，包括红、绿、蓝、红外四个波段。

请见图1，本发明提供的一种结合全球DEM与立体视觉的卫星影像云检测方法，包括以下步骤：

步骤1：将遥感卫星获得的卫星影像(如天绘立体测绘卫星影像)密集匹配获得的M-DSM(MatchingDigitalSurfaceModel)；

请见图2，是本发明实施例的多影像多基元法流程图，本实施例的密集匹配方法采用的是多影像多匹配基元法(MPM)，采用由粗到细(coarsetofine)的多级影像匹配策略，基于多影像的匹配算法，并结合多种匹配基元的匹配结果，所述的多种匹配基元包括特征点、格网点和特征线。请见图3，是本发明实施例获得的M-DSM。

步骤2：将获得的M-DSM与现有的全球DEM数据SRTM(ShuttleRadarTopographyMission)通过地形配准后进行高程比较得到DEM的差值图；

请见图4，其具体实现包括以下子步骤：

步骤2.1：分别在全球DEM数据SRTM与M-DSM上提取地形特征线；

步骤2.2：通过投票算法实现两者的地形特征线的整体套合；

步骤2.3：进行特征点局部的最小二乘相关匹配提高匹配的精度；

步骤2.4：使用地形匹配获得的同名点解算M-DSM与全球DEM数据SRTM之间的空间相似变换关系式，将全球DEM数据SRTM相对于M-DSM进行校准；

步骤2.5：校准后进行高程比较得到DEM的差值图。

请见图5，本实施例采用基于地形表面流水物理模拟分析的方法提取山脊线和山谷线作为地形特征线，本发明采用一种投票算法实现地形特征线的整体套合，实现粗匹配。这种投票算法是霍夫变换的一种改进，投票是根据地理特征线上每个点(地形特征线末梢点除外)，计算与其在同一地形特征线上的下一个相邻点之间的连线与水平方向的夹角，设定一个数值较小的角度阈值(如5°)，以M-DSM的地形特征线中的点为依据，在全球DEM数据SRTM的地形特征线中搜索与其角度差小于阈值的匹配点。例如，对于M-DSM的地形特征线上的点M，其对应了一个角度值α_M，计算其与参考DEM的地形特征线中的点(如点N)所对应的角度值α_N的差的绝对值，若|α_M-α_N|小于阈值，则将点N选为点M的匹配备选点。根据M与N点的坐标差值来投票，对于所有与点M相匹配的点，都按上述方法在其对应的投票区投上一票，即完成点M的投票。对于单个点来说，正确的投票只有一票，而错误的投票可能有多票，正确的投票信息会被错误的投票信息淹没。然而将所有点的投票信息汇总之后，正确的投票会积累起来，而错误的投票往往是随机的，总体上呈现出均匀分散的状态，正确的投票会凸显出来，从而得到正确的结果。最后得到M-DSM与参考DEM的整体坐标偏移。最后再使用最小二乘方法得到高精度的匹配结果。本发明将全球DEM数据SRTM相对于M-DSM进行校准的方法是首先使用地形匹配获得的同名点解算M-DSM与全球DEM数据SRTM之间的空间相似变换关系式。空间相似变换关系为：

$\left[\begin{array}{c}{X}_{tp}\\ Y_{tp}\\ Z_{tp}\end{array}\right]=\mathrm{\λ}\·\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]\·\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]+\left[\begin{array}{c}\mathrm{\Δ}X\\ \mathrm{\Δ}Y\\ \mathrm{\Δ}Z\end{array}\right]$

其中假设M-DSM上某点的坐标为(X，Y，Z)，该点在参考DEM上对应的同名点的坐标为(X_tp,Y_tp,Z_tp)，λ是缩放系数，R是由角元素Φ，Ω，Κ的函数组成的方向余弦，a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为R的9个元素，(ΔX,ΔY,ΔZ)是坐标平移量；空间相似变换的求解采用间接平差法，求解参数的同时还可以探测并剔除匹配结果中的粗差，最后平差解算得到的中误差可以用来代表地形匹配的精度。

步骤3：在DEM差值图中去除误差大于200米的点，然后剩下的空间点统计出单位权中误差(RMSE)，将大于三倍单位权中误差的点通过卫星影像三线阵正视影像的RPC参数投影到卫星影像上，得到疑似云区域种子点的像片坐标；

配准M-DSM与全球DEM数据SRTM之后，通过比较M-DSM与校准后的全球DEM数据SRTM之间的差值得到疑似云区域种子点。因为不同地形的高程范围变化很大，所以用统一的阈值会出现误判断的问题。一种可取的自适应算法是把M-DSM与校准后的全球DEM数据SRTM进行比较，首先根据输入的一个初始高程阈值剔除误差较大的点，根据剩下的空间点统计出单位权中误差(RMSE)：

${\mathrm{\&delta;}}_z=\sqrt{\frac{\underset{i=1}{\overset{N}{\&Sigma;}}{(Z_i-Z_{ref})}^2}{N}},if\left|Z_i-Z_{ref}\right|<Z_{threshold};$

其中，N为总点数，Z_i为第i个点的高程，Z_ref为第i个点在校准后的参考DEM上对应的高程，Z_threshold输入的初始高程阈值。新的自适应阈值则按照下式确定：

Z′_threshold＝3·δ_z

疑似云区域种子点是比较M-DSM产品和校准后的全球DEM数据SRTM，提取其中高程差大于自适应阈值的格网点，通过三线阵正视影像的RPC参数将这些点投影到影像上，得到种子点的像片坐标。

步骤4：去除增长区域中种子点灰度值小于100的低亮度区域、纹理缺乏区域的种子点，根据影像灰度信息对剩余疑似云区域的种子点进行区域增长和分类；请见图6，是本实施例的云区域自动识别流程图；包括以下步骤：

步骤4.1：在将候选云区域种子点投影到三线阵正射影像后，进行影像分析，去除去除增长区域中种子点灰度值小于100的低亮度区域、纹理缺乏区域的种子点。低亮度纹理缺乏区域中的种子点可以通过统计种子点集合的灰度直方图予以清除，如图7所示，是本发明实施例的典型无云影像(a)和有云影像(b)中种子点区域集合的灰度直方图示例，主要过程为判断影像是否含有云，如果存在云，再除去低亮度缺乏纹理区域的种子点。

步骤4.2：进行种子点区域增长。基本思想是将具有相似性质的像素集合起来构成区域。区域增长的具体过程为：以种子点像素作为生长的起点，然后将种子点像素周围邻域中与种子像素有相同或相似性质的像素(根据某种事先确定的生长或相似性准则来判定)合并到种子像素所在的区域中，将这些新像素当作新的种子像素继续进行上面的操作，直到再没有满足条件像素可以包括进来。本发明区域增长的合并准则是每次合并时的灰度差值是和初始的种子点像素进行比较小于设置的阈值。本发明所述的区域增长的合并准则是每次合并时的灰度差值是和初始的种子点像素进行比较小于设置的阈值(一般设置为10或者15)。

步骤5：对增长区域进行云雪分类，得到最终的云区域。

在提取种子点时，云区域中会出现大量种子点，而雪区域中则只会有少量种子点。这样通过这些种子点进行区域增长后，提取的云区域中种子点在整个区域中所占比例很大，而雪区域中种子点在整个区域中所占比例很小。因此可以统计提取的区域中种子点的比例值即可将云和雪区分开来，比值为区域中种子点的个数与区域总像素个数的比，小于阈值(如0.1)可以认为区域不是云区域。

在含有大面积的积雪情况下，因为在全色影像上，云和雪的灰度很相似，区域增长后的云区域会受雪的干扰。最后通过高程差值集合灰度差值区域增长，精化最后的云区域。

请见图8、图9和图10，可以看出本发明方法对少量云分布、大面积云分布、有雪分布时云区域以及云和雪混叠时云区域的云检测都有很准确的检测结果。可以发现，本发明方法能较好解决云检测是如何应对雪的干扰，利用立体视觉结合全球DEM创新地较好地解决了云检测问题。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (10)

1.一种结合全球DEM与立体视觉的卫星影像云检测方法，其特征在于，包括以下步骤：

步骤1：将遥感卫星获得的卫星影像密集匹配获得的M-DSM；

步骤2：将获得的M-DSM与现有的全球DEM数据SRTM通过地形配准后进行高程比较得到DEM的差值图；

步骤3：在DEM差值图中去除误差大于200米的点，然后将剩下的空间点统计单位权中误差(RMSE)，将大于三倍单位权中误差的点通过卫星影像三线阵正视影像的RPC参数投影到卫星影像上，得到疑似云区域种子点的像片坐标；

步骤4：去除增长区域中种子点灰度值小于100的低亮度区域、纹理缺乏区域的种子点，根据影像灰度信息对剩余疑似云区域的种子点进行区域增长和分类；

步骤5：去除雪区域，得到最终的云区域。

2.根据权利要求1所述的结合全球DEM与立体视觉的卫星影像云检测方法，其特征在于：步骤1中所述的密集匹配方法采用的是多影像多匹配基元法(MPM)，采用由粗到细(coarsetofine)的多级影像匹配策略，基于多影像的匹配算法，并结合多种匹配基元的匹配结果，所述的多种匹配基元包括特征点、格网点和特征线。

3.根据权利要求1所述的结合全球DEM与立体视觉的卫星影像云检测方法，其特征在于，步骤2的具体实现包括以下子步骤：

步骤2.1：分别在全球DEM数据SRTM与M-DSM上提取地形特征线；

步骤2.2：通过投票算法实现两者的地形特征线的整体套合；

步骤2.3：进行特征点局部的最小二乘相关匹配提高匹配的精度；

步骤2.4：使用地形匹配获得的同名点解算M-DSM与全球DEM数据SRTM之间的空间相似变换关系式，将全球DEM数据SRTM相对于M-DSM进行校准；

步骤2.5：校准后进行高程比较得到DEM的差值图。

4.根据权利要求3所述的结合全球DEM与立体视觉的卫星影像云检测方法，其特征在于：步骤2.2中所述的投票算法，投票是根据地理特征线上除开地形特征线末梢点外的每个点，计算与其在同一地形特征线上的下一个相邻点之间的连线与水平方向的夹角，设定角度阈值α，α≤5°，以M-DSM的地形特征线中的点所对应的角度为依据，在全球DEM数据SRTM的地形特征线中搜索与其角度差小于阈值的匹配点作为备选匹配点。

5.根据权利要求3或4所述的结合全球DEM与立体视觉的卫星影像云检测方法，其特征在于：步骤2.4中所述使用地形匹配获得的同名点解算M-DSM与参考DEM之间的空间相似变换关系式；空间相似变换关系为式：

$\left[\begin{array}{c}{X}_{tp}\\ Y_{tp}\\ Z_{tp}\end{array}\right]=\mathrm{\&lambda;}\&CenterDot;\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]\&CenterDot;\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]+\left[\begin{array}{c}\mathrm{\&Delta;}X\\ \mathrm{\&Delta;}Y\\ \mathrm{\&Delta;}Z\end{array}\right];$

其中假设M-DSM上某点的坐标为(X，Y，Z)，该点在参考DEM上对应的同名点的坐标为(X_tp,Y_tp,Z_tp)，λ是缩放系数，R是由角元素Φ，Ω，Κ的函数组成的方向余弦，a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为R的9个元素，(ΔX,ΔY,ΔZ)是坐标平移量；空间相似变换的求解采用间接平差法，求解参数的同时探测并剔除匹配结果中的粗差，最后平差解算得到的中误差用来代表地形匹配的精度。

6.根据权利要求1所述的结合全球DEM与立体视觉的卫星影像云检测方法，其特征在于：步骤3中所述的疑似云区域种子点是比较遥感影像密集匹配生成的M-DSM和步骤2.4校准后的全球DEM数据SRTM，提取其中高程差大于三倍单位权中误差(RMSE)的格网点，通过三线阵正视影像的RPC参数将这些点投影到影像上，得到种子点的像片坐标。

7.根据权利要求1所述的结合全球DEM与立体视觉的卫星影像云检测方法，其特征在于，步骤4中所述的区域增长的具体过程为：以种子点像素作为生长的起点，然后将种子点像素周围邻域中与种子像素有相同或相似性质的像素合并到种子像素所在的区域中，将这些新像素当作新的种子像素继续进行上面的操作，直到再没有满足条件像素可以包括进来；本发明所述的区域增长的合并准则是每次合并时的灰度差值是和初始的种子点像素进行比较小于设置的阈值。

8.根据权利要求7所述的结合全球DEM与立体视觉的卫星影像云检测方法，其特征在于：所述的设置的阈值为10或15。

9.根据权利要求1所述的结合全球DEM与立体视觉的卫星影像云检测方法，其特征在于：步骤5中所述的云区域和雪区域的区分主要依据区域中种子点的个数与区域总像素个数的比，小于预设的阈值时则认为区域不是云区域。

10.根据权利要求9所述的结合全球DEM与立体视觉的卫星影像云检测方法，其特征在于：所述的设置的阈值为0.1。