CN110287898B - 一种光学卫星遥感影像云检测方法 - Google Patents

一种光学卫星遥感影像云检测方法 Download PDF

Info

Publication number
CN110287898B
CN110287898B CN201910568174.0A CN201910568174A CN110287898B CN 110287898 B CN110287898 B CN 110287898B CN 201910568174 A CN201910568174 A CN 201910568174A CN 110287898 B CN110287898 B CN 110287898B
Authority
CN
China
Prior art keywords
cloud
area
remote sensing
value
detection
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910568174.0A
Other languages
English (en)
Other versions
CN110287898A (zh
Inventor
闫小明
胡旭东
尹烁
姚燕
胡晓东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Suzhou Zhongketianqi Remote Sensing Technology Co ltd
Original Assignee
Suzhou Zhongketianqi Remote Sensing Technology Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Suzhou Zhongketianqi Remote Sensing Technology Co ltd filed Critical Suzhou Zhongketianqi Remote Sensing Technology Co ltd
Priority to CN201910568174.0A priority Critical patent/CN110287898B/zh
Publication of CN110287898A publication Critical patent/CN110287898A/zh
Application granted granted Critical
Publication of CN110287898B publication Critical patent/CN110287898B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/20—Image preprocessing
    • G06V10/26—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion
    • G06V10/267—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion by performing operations on regions, e.g. growing, shrinking or watersheds
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/56—Extraction of image or video features relating to colour
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/74—Image or video pattern matching; Proximity measures in feature spaces
    • G06V10/75—Organisation of the matching processes, e.g. simultaneous or sequential comparisons of image or video features; Coarse-fine approaches, e.g. multi-scale approaches; using context analysis; Selection of dictionaries
    • G06V10/758—Involving statistics of pixels or of feature values, e.g. histogram matching
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/13—Satellite images
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Multimedia (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Probability & Statistics with Applications (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Astronomy & Astrophysics (AREA)
  • Remote Sensing (AREA)
  • Evolutionary Biology (AREA)
  • Health & Medical Sciences (AREA)
  • Computing Systems (AREA)
  • Databases & Information Systems (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Image Analysis (AREA)

Abstract

本发明揭示了一种光学卫星遥感影像云检测方法，包括统计遥感影像灰度直方图并对其进行预处理；根据灰度直方图计算高斯混合模型分量和迭代最大类间方差阈值，进一步计算混合阈值；使用混合阈值对遥感影像进行分割，获得云初检图；获取云初检图中云区域的轮廓边界特征值和云阴影匹配特征值，进一步去除误检测的云区域，获得误检修正图；对云初检图进行漏检测薄云补偿，获得薄云补偿图；将误检修正图与薄云补偿图进行综合处理获得最终云检测图。本发明无需人工或其他辅助数据参与即可快速、准确地检测出云区域，并在检测过程中去除误检的云区域和对漏检的薄云区域进行补偿，提高云检测精度，适用于全色遥感影像和多光谱遥感影像。

Description

一种光学卫星遥感影像云检测方法
技术领域
本发明涉及光学卫星遥感影像处理技术领域，尤其是涉及一种光学卫星遥感影像云检测方法。
背景技术
光学卫星遥感影像具有广泛的应用领域，如应用于导航定位、环境保护和资源利用等，为陆地、大气和海洋等综合观测提供了有效的数据保障。光学卫星传感器在成像过程中容易受到多种因素影响，尤其是云层的影响。全球约50％以上的地球表面都被云层覆盖，遥感卫星影像拍摄地位于云层之上，导致遥感影像中覆盖大量的云区域，造成地物有效像素丢失，大大降低了遥感影像的质量。为使遥感影像满足智能化处理需求，需在预处理过程中检测遥感影像中的云区域并对其进行替换或者去除。
目前云检测方法多种多样，如物理阈值法、基于空间特性法、模式识别法，以及以上方法综合优化后得到的方法。其中，物理阈值法通过分析云区域的光谱特性，使用不同波段阈值组合进行云区分割，虽然该方法相对较简单，但难以对高亮地物和云进行区分，且容易造成薄云区域的漏检；基于空间特性法通过利用图像光谱亮度值空间变化的特性，使用纹理特征值和光谱变化量进行云区分割，然而对于纹理与云区纹理相似的地物容易误检，并且该方法实施时相对复杂；模式识别法通过提取影像特征，同时使用大批量人工标注数据进行训练，数据量大且费时、费力，无法满足云检测系统时间限制要求。另外，还有利用不同时相同一目标区域的多期数据进行变化监测，基于变化监测的结果进行云提取，此种方法虽然精度高，但对数据要求较高，影像需具备精确的地理位置信息或者匹配点信息。
发明内容
本发明的目的在于克服现有技术的缺陷，提供一种无需人工或其他辅助数据参与即可快速、准确地检测出云区域，并在检测过程中去除误检的云区域和对漏检的薄云区域进行补偿的光学卫星遥感影像云检测方法。
为实现上述目的，本发明提出如下技术方案：一种光学卫星遥感影像云检测方法，包括如下步骤：
S100，统计光学卫星遥感影像的灰度直方图，并对所述灰度直方图进行预处理；
S200，根据预处理后的灰度直方图计算高斯混合模型分量和迭代最大类间方差阈值，进一步根据所述高斯混合模型分量和迭代最大类间方差阈值计算混合阈值L；
S300，使用所述混合阈值L对原始光学卫星遥感影像进行分割，获得云初检图；
S400，获取所述云初检图中云区域的轮廓边界特征值和云阴影匹配特征值，并根据所述轮廓边界特征值和云阴影匹配特征值去除误检测的云区域，获得误检修正图；
S500，对所述云初检图进行漏检测薄云补偿，获得薄云补偿图；
S600，将所述误检修正图与薄云补偿图进行综合处理，获得最终云检测图。
优选地，在步骤S100中，对所述灰度直方图进行预处理包括如下步骤：
对所述灰度直方图进行截断处理，获得直方图有效区间为[X_min，X_max]，进一步通过一维平滑模板对截断后的灰度直方图进行一维卷积平滑处理。
优选地，所述一维平滑模板为[1/10，1/5，2/5，1/5，1/10]。
优选地，在步骤S200中，高斯混合模型分量和迭代最大类间方差阈值通过如下步骤获得：
步骤S201，设置K-means迭代算法的初始均值，并根据所述初始均值计算灰度直方图区间[X_min，X_max]各级数和各级数对应的频数与初始均值和初始均值对应的频数之间的欧式距离，根据欧式距离的大小对灰度直方图区间进行归类，并根据归类结果重新计算每一类的均值，使用重新计算后的每一类均值作为下一次迭代的初始值并继续迭代，直至两次迭代的均值之差小于预设限差值时停止迭代，此时迭代结果作为高斯混合模型参数的初始均值；
步骤S202，构建对数似然函数模型，并使用期望最大算法对高斯分量进行拟合，获得高斯混合模型分量X1_GMM，X2_GMM和X3_GMM；
步骤S203，使用迭代的最大类间方差算法计算迭代最大类间方差阈值X_OTSU。
优选地，在步骤S203中，迭代最大类间方差阈值X_OTSU通过如下步骤获得：
步骤S203a，将灰度直方图区间[X_min，X_max]分为[X_min，X]和[X，X_max]，分别计算两个区间的最大类间方差阈值T1和T2，其中，X为X_min和X_max的均值；
步骤S203b，判断T1与T2之差是否小于预设限差值，若是，则T1与T2的均值为迭代最大类间方差阈值X_OTSU，否则，将灰度直方图区间[X_min，X_max]更新为[T1，T2]，重复执行步骤S203b。
优选地，所述混合阈值L通过如下步骤获得：
计算迭代最大类间方差阈值X_OTSU与高斯混合模型分量X2_GMM和X3_GMM的欧式距离d1和d2；
若d1大于d2，则[X_OTSU，X2_GMM]的最大类间方差阈值为混合阈值L；否则，直方图区间[X_OTSU，X3_GMM]的最大类间方差阈值为混合阈值L。
优选地，在步骤S400中，轮廓边界特征值包括云轮廓多边形外角为凹角的数量和云轮廓多边形外角为锐角的数量，所述云轮廓多边形外角为凹角的数量和云轮廓多边形外角为锐角的数量通过如下步骤获得：
S401，使用多边形近似算法对云初检图中云区域进行多边形近似，获得云轮廓多边形；
S402，计算云轮廓多边形的外角和，并进一步根据所述外角和计算多边形的外角为凹角的数量count_i，其中，外角的角度通过如下公式计算：
li＝arctan(yi-yi-1)/(xi-xi-1)-arctan(yi+1-yi)/(xi+1-xi)
其中，(xi，yi)为外角对应的两条边的多边形顶点坐标；
凹角的数量count_i通过如下公式计算：
Figure BDA0002110096110000041
S403，统计云轮廓多边形的外角为锐角的数量Angel_i，其中，通过如下方式进行判断外角是否为锐角；以一个顶点为起点，其相邻左右两个端点为终点，分别形成两个向量，进一步计算两个向量的叉乘，若计算结果小于0，则为锐角。
优选地，误检测的云区域包括高亮山地雪地物，所述高亮山地雪地物通过如下步骤去除：判断count_i是否大于J1且Angel_i是否大于J2，若是，则云轮廓多边形为误检的高亮山地雪地物，进一步去除所述高亮山地雪地物，其中，J1为预设外角为凹角的数量，J2为预设外角为锐角的数量。
优选地，误检测的云区域包括高亮裸地，所述高亮裸地通过如下步骤去除：
计算云阴影相对于云的方向；
使用云阴影模板在云阴影相对于云的方向上迭代步长，统计云阴影模板在移动过程中覆盖区域的像素值均值，并在像素值均值最小时记录此时的步长，进一步统计每个云轮廓对应的步长并进行排序，取排序后步长的中值作为云阴影区域距离，统计在所述云阴影区域距离下云阴影模板覆盖原始遥感影像对应区域的像素值，取统计结果下分位五分之一的像素值作为云阴影阈值；
根据云阴影阈值确定云阴影匹配特征值，当云阴影区域距离下云阴影模板覆盖的小于云阴影阈值的像素个数是云阴影模板像素个数的一半以上时，云阴影匹配特征值为1，否则云阴影匹配特征值为0，对云阴影匹配特征值为0的云轮廓判定为误检的高亮裸地，进一步去除所述高亮裸地。
优选地，漏检测薄云补偿包括如下步骤：
步骤S501，降低混合阈值L到Ln，并使用混合阈值L和混合阈值Ln对光学卫星遥感影像进行阈值分割，将位于区间[Ln，L]的像素值预测为薄云区域，生成薄云预测图，其中，Ln＝0.9*L；
步骤S502，计算薄云预测图中每个薄云预测区域内的暗像素，生成暗像素图；
步骤S503，通过公式D＝N/S计算暗像素图中每个薄云预测区域内的暗像素密度，并在暗像素密度小于预设暗像素密度阈值时薄云预测区域为真实薄云区域，进一步对真实薄云区域进行补偿处理，其中，N为薄云预测区域内暗像素数，S为薄云预测区域面积。
优选地，所述暗像素通过如下步骤获得：将预设暗像素模板在薄云预测区域内移动，进一步将预设暗像素模板内最小值小于混合阈值Ln的像素作为暗像素。
本发明的有益效果是：
(1)本发明无需人工或其他辅助数据参与即可快速、准确地检测出遥感影像中的云区域，并在检测云区域的过程中去除误检的云区域，并对漏检的薄云区域进行补偿，最终获得精确的云检测图，适用于全色遥感影像和多光谱遥感影像；
(2)本发明通过云初检图中云区域的轮廓边界特征值和云阴影匹配特征值去除误检测的云区域，同时通过暗像素密度对误排除的薄云区域进行补偿，有效解决了光谱阈值法中对于高亮地物的误检和对于薄云漏检的问题，并且去除物检测的云区域和对漏检测的薄云区域进行补偿并不依赖于紧密像素亮度关系，可先对初检结果进行重采样，获得尺度较小的图像后再进行分割，能够有效提高云检测效率。
附图说明
图1是本发明的流程图示意图；
图2是本发明的方法流程图示意图；
图3是本发明的原始遥感影像示意图；
图4是本发明的云初检图示意图；
图5是本发明的误检修正图一示意图；
图6是本发明的误检修正图二示意图；
图7是本发明的暗像素密度图示意图；
图8是本发明的综合结果图示意图；
图9是本发明的最终云检测图示意图。
具体实施方式
下面将结合本发明的附图，对本发明实施例的技术方案进行清楚、完整的描述。
本发明所揭示的一种光学卫星遥感影像云检测方法，无需人工或其他辅助数据参与即可快速、准确地检测出遥感影像中的云区域，并在检测云区域的过程中去除误检的云区域，并对漏检的薄云区域进行补偿，最终获得精确的云检测图，适用于全色遥感影像和多光谱遥感影像。
结合图1和图2所示，为本发明所揭示的一种光学卫星遥感影像云检测方法，包括如下步骤：
步骤S100，统计光学卫星遥感影像的灰度直方图，并对所述灰度直方图进行预处理。
具体地，光学卫星遥感影像包括全色遥感影像和多光谱遥感影像，其中，多光谱遥感影像具有蓝(B)、绿(G)、红(R)和近红外(NIR)四个波段，蓝波段光谱特性良好，易于区分影像上的云和地物，适合进行云检测，全色遥感影像具有一个波段。本实施例中，多光谱遥感影像的第一波段，也即蓝(B)波段进行灰度直方图统计，全色遥感影像则统计单波段灰度直方图。
灰度直方图预处理包括截断和降噪处理。具体地，取灰度直方图左右两端各占a％的像素所在灰阶作为灰度直方图主体区间上下边界，其中，a为0.01～2。实施时，从下灰阶和上灰阶处对灰度直方图进行截断处理，获得灰度直方图有效区间[X_min，X_max]。进一步通过一维平滑模板对截断后的灰度直方图进行一维卷积平滑处理，实现对灰度直方图的降噪。
本实施例中，一维平滑模板为[1/10，1/5，2/5，1/5，1/10]。
步骤S200，根据预处理后的灰度直方图计算高斯混合模型分量和迭代最大类间方差阈值，进一步根据所述高斯混合模型分量和迭代最大类间方差阈值计算混合阈值L。
步骤S300，使用所述混合阈值L对原始光学卫星遥感影像进行阈值分割，获得云初检图。
具体地，高斯混合模型分量和迭代最大类间方差阈值通过如下步骤获得：
首先，设置K-means迭代算法的初始均值，并根据所述初始均值计算灰度直方图区间[X_min，X_max]各级数和各级数对应的频数与初始均值和初始均值对应的频数之间的欧式距离，根据欧式距离的大小对灰度直方图区间进行归类，并根据归类结果重新计算每一类的均值，使用重新计算后的每一类均值作为下一次迭代的初始值，继续迭代，直至两次迭代的均值之差小于预设限差值时停止迭代，此时K-means迭代算法的迭代结果作为高斯混合模型参数的初始均值，也即作为μk的初始值，μk为高斯混合模型第k个均值分量。
实施时，设置K-means迭代算法的初始均值为u1，u2，u3，初始值分别为1/4M，1/2M和3/4M，其中，M为灰度直方图预处理后的区间长度。进一步根据如下公式计算灰度直方图区间各级数和各级数对应的频数与初始均值和初始均值对应的频数之间的欧式距离：
d＝(p-ui)2+(fp-fi)2
其中，p为灰度直方图级数，p＝X_min，X_min+1，X_min+2，...，X_max，fp为灰度直方图第p级数对应的频数；ui为初始均值，fi为第i个初始均值对应的频数，i＝1，2，3。
对于三个类C1，C2，C3，直方图级数p取d最小的类集合Ci加入，计算分类完成的各个集合的均值作为下一次迭代的初始值，继续迭代，直到两次迭代的均值之差小于限差即可，限差取级数最小值1，迭代结果Ci中的级数均值作为高斯混合模型分量的初始均值μ1 (0)，μ2 (0)，μ3 (0)
其次，构建对数似然函数模型，并使用期望最大算法对高斯分量进行拟合，获得高斯混合模型分量X1_GMM，X2_GMM和X3_GMM。
具体地，对数似然函数模型如下所示：
Figure BDA0002110096110000081
其中，μk，αk，
Figure BDA0002110096110000082
分别表示等k个分模型对应的均值，权重和方差，
Figure BDA0002110096110000083
为当前模型参数下第j个观测数据来自第k个分模型的概率，K＝3。
使用期望最大算法对参数进行拟合，下一次迭代对数似然函数：
θi+1＝argmax(Q(θ，θi))
下一次迭代对数似然函数的参数解为：
Figure BDA0002110096110000091
其中，
Figure BDA0002110096110000092
Figure BDA0002110096110000093
N为灰度直方图总级数，xj为灰度直方图级数，yj为灰度直方图级数对应的频数，迭代终止条件为：
Figure BDA0002110096110000094
实施时，初始方差分量和权重分量根据K-means算法初始化的高斯混合模型均值分量得出：
Figure BDA0002110096110000095
Figure BDA0002110096110000096
使用上述公式迭代即可获得高斯混合模型三个均值分量X1_GMM，X2_GMM和X3_GMM。
最后，使用迭代的最大类间方差算法计算迭代最大类间方差阈值X_OTSU。
具体地，首先，将灰度直方图区间[X_min，X_max]二分为[X_min，X]和[X，X_max]，X为X_min和X_max的均值，分别计算两个区间的最大类间方差阈值T1和T2，计算公式如下；
σ＝pA pB(uA-u)2
其中，u表示迭代级数左侧直方图的均值，uA表示迭代级数右侧直方图的均值，pA表示迭代级数左侧直方图的权重，为左侧所有频数的和与所有频数和的比值，pB表示迭代级数右侧直方图的权重，为右侧侧所有频数的和与所有频数和的比值，当σ最大时对应直方图级数为最大类间方差阈值。
最后，若两个最大类间方差阈值T1与T2之差小于预设限差值时，则两个最大类间方差阈值的均值为迭代最大类间方差阈值X_OTSU，否则，将灰度直方图区间[X_min，X_max]更新为[T1，T2]，重复上一步，直至满足条件时迭代停止，实施时，预设限差值优选5，当然，可以根据实际需求进行设定。
进一步地，混合阈值L根据如下步骤获得：
首先，计算迭代最大类间方差阈值X_OTSU与高斯混合模型分量X2_GMM和X3_GMM的欧式距离d1和d2；
最后，判断d1是否大于d2。
当d1大于d2时，根据如下公式：
σ＝pA pB(uA-u)2
计算直方图区间[X_OTSU，X2_GMM]的最大类间方差阈值，该最大类间方差阈值为混合阈值L；
当d1小于或等于d2时，根据如下公式：
σ＝pA pB(uA-u)2
计算直方图区间[X_OTSU，X3_GMM]的最大类间方差阈值，该最大类间方差阈值为混合阈值L。
获得混合阈值L后进一步使用混合阈值L对原始光学卫星遥感影像(图3所示)进行阈值分割，获得云初检图，如图4所示。实施时，如果X_OTSU小于或等于500，则选择经验阈值作为最终的混合阈值L。
步骤S400，获取所述云初检图中云区域的轮廓边界特征和云阴影匹配特征，并根据所述轮廓边界特征和云阴影匹配特征去除误检测的云区域，获得误检修正图。
具体地，高亮山地雪地物的亮度与云区域的亮度基本相同，容易将高亮山地雪地物误检测为云区域，但高亮山地雪地物的边界与云区域的边界存在明显不同。高亮山地雪地物的边界呈现锥形状的边界点较多，并且边界凹凸度较大，而云区域的边界较为平滑，呈现锥形状的边界点较少。因此，可根据云区域的轮廓边界特征对误检测的高亮山地雪区域进行去除。
本实施例中，轮廓边界特征值包括云轮廓多边形外角为凹角的数量和云轮廓多边形外角为锐角的数量。
误检测的高亮山地雪地物通过如下步骤去除：
首先，使用多边形近似算法对云初检图中云区域进行多边形近似，获得云轮廓多边形；
其次，计算云轮廓多边形的外角和，并进一步根据所述外角和计算多边形的外角为凹角的数量count_i，其中，外角的角度通过如下公式计算：
li＝arctan(yi-yi-1)/(xi-xi-1)-arctan(yi+1-yi)/(xi+1-xi)
其中，(xi，yi)为外角对应的两条边的多边形顶点坐标，(xi-1，yi-1)为一条边的终点坐标，(xi+1，yi+1)为另一条边的终点坐标；
凹角的个数通过如下公式计算：
Figure BDA0002110096110000111
再其次，统计云轮廓多边形的外角为锐角的数量Angel_i，其中，外角锐角可通过如下方式进行判断；以一个顶点为起点，其相邻左右两个端点为终点，分别形成两个向量，进一步计算两个向量的叉乘，若计算结果小于0，则为锐角。
最后，判断count_i是否大于J1且Angel_i是否大于J2时，若是，则云轮廓多边形为误检的高亮山地雪地物，进一步对其进行标记，剔除误检轮廓后进行多边形轮廓重绘，生成误检修正图一，如图5所示，其中，J1为预设外角凹角数，J2为预设外角锐角数。
进一步地，高亮裸地区域的亮度与云区域的亮度也基本相同，容易将高亮裸地区域误检测为云区域，但高亮裸地区域不具有阴影，而云区域具有阴影，因此，可根据云阴影匹配特征来对误检测的高亮裸地区域进行去除。
误检测的高亮裸地区域通过如下步骤去除：
首先，计算云阴影相对于云的方向w。具体地，根据太阳方位角计算云阴影相对于云的方向w，
Figure BDA0002110096110000121
其中，
Figure BDA0002110096110000122
为原始影像当时拍摄条件下的太阳方位角。
其次，使用云阴影模板在云阴影相对云的方向w上迭代步长，统计云阴影模板在移动过程中覆盖区域的像素值均值，并在像素值均值最小时记录此时的步长，进一步统计每个云轮廓对应的步长并进行排序，取排序后的步长的中值作为云阴影区域距离，统计在所述云阴影区域距离下云阴影模板覆盖原始遥感影像对应区域的像素值，并取统计结果下分位五分之一的像素值作为云阴影阈值；
最后，根据云阴影阈值确定云阴影匹配特征值。当云阴影区域距离下云阴影模板覆盖的小于云阴影阈值的像素个数是云阴影模板像素个数的一半以上时，云阴影匹配特征值为1，否则云阴影匹配特征值为0，对云阴影匹配特征值为0的云轮廓认定为误检的高亮裸地，进行误检去除，生成误检修正图二，如图6所示。
进一步地，根据误检修正图一和误检修正图二可获得误检修正图。
S500，对所述云初检图进行漏检测薄云补偿，获得薄云补偿图。
具体地，薄云区域的亮度与高亮裸地基本相同，在进行初检时容易将薄云区域遗漏，导致生成的云初检图不准确，因此，需对云初检图中漏掉的薄云区域进行补偿。
由于薄云区域不能形成云阴影，并且薄云区域相对于自然地物较为光滑，亮度值偏高且几乎不存在暗像素，而自然地物则由于地物相互遮挡的关系，暗像素较多。因此，可以通过薄云区域暗像素先验知识对薄云进行漏检补偿。通过如下步骤对薄云漏检区域进行补偿：
首先，降低混合阈值L到Ln，其中，Ln＝0.9*L，进一步使用混合阈值L和混合阈值Ln对光学卫星遥感影像进行阈值分割，将位于区间[Ln，L]的像素值预测为薄云区域，生成薄云预测图；
其次，计算薄云预测图中薄云预测区域内的暗像素，生成暗像素图，如图7所示。暗像素通过如下方式计算：将预设暗像素模板在薄云预测区域内移动，进一步将预设暗像素模板内最小值小于混合阈值Ln的像素作为暗像素。实施时，可对薄云预测图中的薄云预测区域进行腐蚀膨胀，去除小面积的薄云预测区域，然后计算薄云预测区域内的暗像素。本实施例中，预设暗像素模板优选3*3的模板，当然，也可以选择其他尺寸的模板，如4*4等，可根据实际需求进行设定。
最后，计算每个薄云预测区域内的暗像素密度D。暗像素密度D通过如下公式计算：D＝N/S，其中，N为薄云预测区域内暗像素数，S为薄云预测区域面积。当暗像素密度小于预设暗像素密度阈值时薄云预测区域为真实薄云区域，进一步对真实薄云区域进行补偿处理，否则，薄云预测区域为真实地物区域。本实施例中，预设暗像素密度阈值为5％。
步骤S600，将所述误检修正图与薄云补偿图进行综合处理获得最终云检测图。
具体地，将误检修正图与薄云补偿图进行合并获得云检二值图，如图8所示，进一步使用形态学方法进行腐蚀、膨胀和腐蚀操作对云检二值图进行形态学处理，去除内部空洞和小面积碎班，获得最终云检测图，如图9所示。
本发明通过云初检图中云区域的轮廓边界特征和云阴影匹配特征值去除误检测的云区域，同时通过暗像素密度对初检测遗漏的薄云区域进行补偿，得到精确的云检二值图，有效解决了光谱阈值法中对于高亮地物的误检和对于薄云漏检的问题，并且去除误检测的云区域和对漏检测的薄云区域进行补偿并不依赖于紧密像素亮度关系，可先对初检结果进行重采样，获得尺度较小的图像，最后再进行分割，能够有效提高云检测效率。
本发明的技术内容及技术特征已揭示如上，然而熟悉本领域的技术人员仍可能基于本发明的教示及揭示而作种种不背离本发明精神的替换及修饰，因此，本发明保护范围应不限于实施例所揭示的内容，而应包括各种不背离本发明的替换及修饰，并为本专利申请权利要求所涵盖。

Claims (10)

1.一种光学卫星遥感影像云检测方法，其特征在于，包括如下步骤：
S100，统计光学卫星遥感影像的灰度直方图，并对所述灰度直方图进行预处理；
S200，根据预处理后的灰度直方图计算高斯混合模型分量和迭代最大类间方差阈值，进一步根据所述高斯混合模型分量和迭代最大类间方差阈值计算混合阈值L；
S300，使用所述混合阈值L对原始光学卫星遥感影像进行分割，获得云初检图；
S400，获取所述云初检图中云区域的轮廓边界特征值和云阴影匹配特征值，并根据所述轮廓边界特征值和云阴影匹配特征值去除误检测的云区域，获得误检修正图；
S500，对所述云初检图进行漏检测薄云补偿，获得薄云补偿图，其中，所述漏检测薄云补偿包括如下步骤：
步骤S501，降低混合阈值L到Ln，并使用混合阈值L和混合阈值Ln对光学卫星遥感影像进行阈值分割，将位于区间[Ln，L]的像素值预测为薄云区域，生成薄云预测图，其中，Ln＝0.9*L；
步骤S502，计算薄云预测图中每个薄云预测区域内的暗像素，生成暗像素图；
步骤S503，通过公式D＝N/S计算暗像素图中每个薄云预测区域内的暗像素密度，并在暗像素密度小于预设暗像素密度阈值时薄云预测区域为真实薄云区域，进一步对真实薄云区域进行补偿处理，其中，N为薄云预测区域内暗像素数，S为薄云预测区域面积；
S600，将所述误检修正图与薄云补偿图进行综合处理，获得云检测图。
2.根据权利要求1所述的光学卫星遥感影像云检测方法，其特征在于，在步骤S100中，对所述灰度直方图进行预处理包括如下步骤：
对所述灰度直方图进行截断处理，获得直方图有效区间为[X_min，X_max]，进一步通过一维平滑模板对截断后的灰度直方图进行一维卷积平滑处理。
3.根据权利要求2所述的光学卫星遥感影像云检测方法，其特征在于，所述一维平滑模板为[1/10，1/5，2/5，1/5，1/10]。
4.根据权利要求1所述的光学卫星遥感影像云检测方法，其特征在于，在步骤S200中，高斯混合模型分量和迭代最大类间方差阈值通过如下步骤获得：
步骤S201，设置K-means迭代算法的初始均值，并根据所述初始均值计算灰度直方图区间[X_min，X_max]各级数和各级数对应的频数与初始均值和初始均值对应的频数之间的欧式距离，根据欧式距离的大小对灰度直方图区间进行归类，并根据归类结果重新计算每一类的均值，使用重新计算后的每一类均值作为下一次迭代的初始值并继续迭代，直至两次迭代的均值之差小于预设限差值时停止迭代，此时迭代结果作为高斯混合模型参数的初始均值；
步骤S202，构建对数似然函数模型，并使用期望最大算法对高斯分量进行拟合，获得高斯混合模型分量X1_GMM，X2_GMM和X3_GMM；
步骤S203，使用迭代的最大类间方差算法计算迭代最大类间方差阈值X_OTSU。
5.根据权利要求4所述的光学卫星遥感影像云检测方法，其特征在于，在步骤S203中，迭代最大类间方差阈值X_OTSU通过如下步骤获得：
步骤S203a，将灰度直方图区间[X_min，X_max]分为[X_min，X]和[X，X_max]，分别计算两个区间的最大类间方差阈值T1和T2，其中，X为X_min和X_max的均值；
步骤S203b，判断T1与T2之差是否小于预设限差值，若是，则T1与T2的均值为迭代最大类间方差阈值X_OTSU，否则，将灰度直方图区间[X_min，X_max]更新为[T1，T2]，重复执行步骤S203a。
6.根据权利要求4所述的光学卫星遥感影像云检测方法，其特征在于，所述混合阈值L通过如下步骤获得：
计算迭代最大类间方差阈值X_OTSU与高斯混合模型分量X2_GMM和X3_GMM的欧式距离d1和d2；
若d1大于d2，则[X_OTSU，X2_GMM]的最大类间方差阈值为混合阈值L；否则，直方图区间[X_OTSU，X3_GMM]的最大类间方差阈值为混合阈值L。
7.根据权利要求1所述的光学卫星遥感影像云检测方法，其特征在于，在步骤S400中，轮廓边界特征值包括云轮廓多边形外角为凹角的数量和云轮廓多边形外角为锐角的数量，所述云轮廓多边形外角为凹角的数量和云轮廓多边形外角为锐角的数量通过如下步骤获得：
S401，使用多边形近似算法对云初检图中云区域进行多边形近似，获得云轮廓多边形；
S402，计算云轮廓多边形的外角和，并进一步根据所述外角和计算多边形的外角为凹角的数量count_i，其中，外角的角度通过如下公式计算：
li＝arctan(yi-yi-1)(xi-xi-1)-arctan(yi+1-yi)(xi+1-xi)
其中，(xi，yi)为外角对应的两条边的多边形顶点坐标；
凹角的数量count_i通过如下公式计算：
Figure FDA0004085238570000031
S403，统计云轮廓多边形的外角为锐角的数量Angel_i，其中，通过如下方式进行判断外角是否为锐角；以一个顶点为起点，其相邻左右两个端点为终点，分别形成两个向量，进一步计算两个向量的叉乘，若计算结果小于0，则为锐角。
8.根据权利要求7所述的光学卫星遥感影像云检测方法，其特征在于，误检测的云区域包括高亮山地雪地物，所述高亮山地雪地物通过如下步骤去除：判断count_i是否大于J1且Angel_i是否大于J2，若是，则云轮廓多边形为误检的高亮山地雪地物，进一步去除所述高亮山地雪地物，其中，J1为预设外角为凹角的数量，J2为预设外角为锐角的数量。
9.根据权利要求1所述的光学卫星遥感影像云检测方法，其特征在于，误检测的云区域包括高亮裸地，所述高亮裸地通过如下步骤去除：
计算云阴影相对于云的方向；
使用云阴影模板在云阴影相对于云的方向上迭代步长，统计云阴影模板在移动过程中覆盖区域的像素值均值，并在像素值均值最小时记录此时的步长，进一步统计每个云轮廓对应的步长并进行排序，取排序后步长的中值作为云阴影区域距离，统计在所述云阴影区域距离下云阴影模板覆盖原始遥感影像对应区域的像素值，取统计结果下分位五分之一的像素值作为云阴影阈值；
根据云阴影阈值确定云阴影匹配特征值，当云阴影区域距离下云阴影模板覆盖的小于云阴影阈值的像素个数是云阴影模板像素个数的一半以上时，云阴影匹配特征值为1，否则云阴影匹配特征值为0，对云阴影匹配特征值为0的云轮廓判定为误检的高亮裸地，进一步去除所述高亮裸地。
10.根据权利要求1所述的光学卫星遥感影像云检测方法，其特征在于，所述暗像素通过如下步骤获得：将预设暗像素模板在薄云预测区域内移动，进一步将预设暗像素模板内最小值小于混合阈值Ln的像素作为暗像素。
CN201910568174.0A 2019-06-27 2019-06-27 一种光学卫星遥感影像云检测方法 Active CN110287898B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910568174.0A CN110287898B (zh) 2019-06-27 2019-06-27 一种光学卫星遥感影像云检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910568174.0A CN110287898B (zh) 2019-06-27 2019-06-27 一种光学卫星遥感影像云检测方法

Publications (2)

Publication Number Publication Date
CN110287898A CN110287898A (zh) 2019-09-27
CN110287898B true CN110287898B (zh) 2023-04-18

Family

ID=68019319

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910568174.0A Active CN110287898B (zh) 2019-06-27 2019-06-27 一种光学卫星遥感影像云检测方法

Country Status (1)

Country Link
CN (1) CN110287898B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111024072B (zh) * 2019-12-27 2021-06-11 浙江大学 一种基于深度学习的卫星地图辅助导航定位方法
CN112102180B (zh) * 2020-08-21 2022-10-11 电子科技大学 一种基于Landsat光学遥感图像的云识别方法
CN112347956B (zh) * 2020-11-12 2022-05-06 上海交通大学 基于多无人机和机器视觉的云观测系统及方法
CN112634349A (zh) * 2020-12-22 2021-04-09 福建省星云大数据应用服务有限公司 一种基于遥感影像的茶园面积估计方法及系统
CN113436092B (zh) * 2021-06-16 2022-04-22 中国电子科技集团公司第五十四研究所 一种全色遥感图像云区域识别方法
CN113298836B (zh) * 2021-06-29 2023-03-14 天津市测绘院有限公司 一种顾及要素轮廓强度的遥感图像薄云去除方法及系统
CN113935917B (zh) * 2021-10-14 2024-06-18 中国石油大学(华东) 一种基于云图运算和多尺度生成对抗网络的光学遥感影像薄云去除方法
CN114332085B (zh) * 2022-03-11 2022-06-24 西安中科西光航天科技有限公司 一种光学卫星遥感影像检测方法
CN115512236B (zh) * 2022-10-13 2023-04-28 昆明理工大学 一种基于K-means ++的Himawari-8多光谱云检测方法及系统
CN116245757B (zh) * 2023-02-08 2023-09-19 北京艾尔思时代科技有限公司 多模态数据的多场景通用性遥感影像云修复方法和系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101894382A (zh) * 2010-07-23 2010-11-24 同济大学 一种集成LiDAR点云的卫星立体影像阴影计算方法
JP2014107589A (ja) * 2012-11-22 2014-06-09 Canon Inc 画像処理装置、画像処理方法、およびプログラム
CN105678777A (zh) * 2016-01-12 2016-06-15 武汉大学 一种多特征联合的光学卫星影像云与云阴影检测方法
CN105894520A (zh) * 2016-04-25 2016-08-24 武汉大学 一种基于高斯混合模型的卫星影像自动云检测方法
CN107564017A (zh) * 2017-08-29 2018-01-09 南京信息工程大学 一种城市高分遥感影像阴影检测及分割方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101894382A (zh) * 2010-07-23 2010-11-24 同济大学 一种集成LiDAR点云的卫星立体影像阴影计算方法
JP2014107589A (ja) * 2012-11-22 2014-06-09 Canon Inc 画像処理装置、画像処理方法、およびプログラム
CN105678777A (zh) * 2016-01-12 2016-06-15 武汉大学 一种多特征联合的光学卫星影像云与云阴影检测方法
CN105894520A (zh) * 2016-04-25 2016-08-24 武汉大学 一种基于高斯混合模型的卫星影像自动云检测方法
CN107564017A (zh) * 2017-08-29 2018-01-09 南京信息工程大学 一种城市高分遥感影像阴影检测及分割方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于高斯混合模型的遥感影像云检测技术;杨帆等;《南京林业大学学报(自然科学版)》;20180321(第04期);第138-144页 *

Also Published As

Publication number Publication date
CN110287898A (zh) 2019-09-27

Similar Documents

Publication Publication Date Title
CN110287898B (zh) 一种光学卫星遥感影像云检测方法
CN113592882B (zh) 一种基于无人机多源遥感的树冠提取方法
CN106780524B (zh) 一种三维点云道路边界自动提取方法
CN109191432B (zh) 基于域变换滤波多尺度分解的遥感图像云检测方法
CN107564017B (zh) 一种城市高分遥感影像阴影检测及分割方法
CN111027446B (zh) 一种高分辨率影像的海岸线自动提取方法
CN110443201B (zh) 基于多源影像联合形状分析与多属性融合的目标识别方法
CN111881816B (zh) 一种长时序的河湖围埂养殖区域监测方法
CN109829423B (zh) 一种结冰湖泊红外成像检测方法
CN111553922B (zh) 一种卫星遥感影像自动云检测方法
CN107992856B (zh) 城市场景下的高分遥感建筑物阴影检测方法
CN112099046B (zh) 基于多值体素模型的机载lidar三维平面检测方法
CN111079596A (zh) 高分辨率遥感影像的海上典型人造目标识别系统及方法
CN110428425B (zh) 一种基于海岸线矢量数据的sar图像海陆分离方法
CN107862271B (zh) 一种舰船目标的检测方法
CN110110687B (zh) 基于颜色信息和三维轮廓信息的树上水果自动识别方法
CN101114337A (zh) 一种地面建筑物识别定位方法
CN116704333B (zh) 一种基于激光点云数据的单株林木检测方法
CN111310771B (zh) 遥感影像的道路图像提取方法、装置、设备及存储介质
CN112669363A (zh) 城市绿地三维绿量测算方法
CN107194405B (zh) 一种交互式半自动高分辨率遥感影像建筑物提取的方法
CN106204596B (zh) 一种基于高斯拟合函数与模糊混合估计的全色波段遥感影像云检测方法
Omidalizarandi et al. Segmentation and classification of point clouds from dense aerial image matching
Aytekin et al. Automatic and unsupervised building extraction in complex urban environments from multi spectral satellite imagery
CN109785318B (zh) 基于面线基元关联约束的遥感图像变化检测方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant