一种基于高斯混合模型的卫星影像自动云检测方法
【专利摘要】本发明公开了一种基于高斯混合模型的卫星影像自动云检测方法，包括：S1统计卫星影像的原始灰度直方图，预处理灰度直方图以消除干扰；S2对灰度直方图构建高斯混合模型；S3 确定分割云区域的灰度阈值；S4 采用灰度阈值分割云区域。本发明方法弥补了传统纹理分析法、同态滤波法、多光谱综合法的不足，可用于卫星影像的筛选与质量控制、云区自动提取、含云影像快速修补等工作。本发明基于单波段影像直方图，不受限于卫星光谱范围，同时适用于全色、多光谱影像；检测精度高，无需辅助信息和人工干预，计算速度快，可自动化处理。
【专利说明】
一种基于高斯混合模型的卫星影像自动云检测方法
技术领域
[0001] 本发明属于遥感图像处理技术领域，特别涉及一种基于高斯混合模型的卫星影像 自动云检测方法。
【背景技术】
[0002] 以Worldview系列、Pleiades系列、中国"高分"系列卫星为代表，高分辨率对地观 测技术迅速发展，卫星重访周期不断变小、卫星影像分辨率和幅宽不断提高、卫星影像的数 据量急剧增加，广泛应用于资源调查、防灾减灾等领域。然而在卫星影像数据处理过程中， 云的存在是一个不可忽视的问题，它不仅对地面场景形成遮挡，还在一定程度上改变了影 像的光谱特性，给影像匹配、匀色和镶嵌工作造成诸多不利影响。如何精准、快速地检测卫 星影像中的云层，从而更好地发挥卫星影像的利用价值，是目前摄影测量与遥感领域的热 点课题。
[0003] 当前较为常用的云检测方法主要有纹理分析法、同态滤波法、多光谱综合法等。但 是，考虑到高分辨率对地观测卫星影像波段少、光谱范围窄(0.45~0.90微米）、数据量大的 特点，现有云检测方法并不能完全适用于高分辨率对地观测卫星影像。纹理分析法基于统 计方法提取云区和非云区的空间特性，可有效识别大片层云，但难以识别纹理较强的卷积 云。同态滤波法不擅长处理厚云影像，且滤波的过程中会丢失一些有用信息，降低影像产品 质量。多光谱综合法利用云和晴空反射率的差异，但它要求传感器配有多个热红外波段，探 测波长需涵盖水或二氧化碳的吸收带，而对地观测卫星影像一般仅含有蓝、绿、红、近红外4 个波段，光谱范围一般为〇. 43~0.90mi，多光谱综合法的优势难以体现。
【发明内容】
[0004] 本发明的目的是针对高分辨率对地观测卫星影像，提供一种高效精确的基于高斯 混合模型的卫星影像自动云检测方法。
[0005] 为达到上述目的，本发明采用如下技术方案：
[0006] -种基于高斯混合模型的卫星影像自动云检测方法，包括步骤：
[0007] S1统计卫星影像的原始灰度直方图，预处理灰度直方图以消除干扰；
[0008] S2对灰度直方图构建高斯混合模型，本步骤进一步包括：
[0009] 2.1采用高斯混合模型表示灰度直方图函数；
[001 0] 2.2采用局部最大值法获得灰度直方图的峰值点，峰值点数即高斯分量数M;
[0011] 2.2使用最大类间方差法获取相邻峰值点间的谷值点；同时，将灰度直方图的两端 点也标记为谷值点，谷值点横坐标记为Vj，j = 1，2，. . .M+1;
[0012] 2.3初始化高斯混合模型的模型参数，获得各高斯分量的权重、均值、标准差的初 始值 rf、<、C其中，足。>义，^ =窆》)(.丫- m个高斯分量峰值点横坐标;h(x)表示灰度直方图中灰阶x对应的像素数；
[0013] 2.4采用最大似然法对模型参数进行迭代拟合，获得高斯混合模型；
[0014] S3确定分割云区域的灰度阈值，本步骤进一步包括：
[0015] 3.1根据预处理后灰度直方图的灰阶范围[Xmin，Xmax]确定云层和地物间的初始分 类阈值
[0016] 3.2设置第m个高斯分量的高频区间|>m-l . 3〇m，lim+l . 3〇m]，iim、〇m分别为第m个高斯 分量的均值和标准差，从子步骤2.4获得的高斯混合模型获得；
[0017] 3.3找出权重最大的高斯分量，记为口(叉|似，〇人）；若口(叉|似，(^)的均值似>￥加，以口 (x | iu，〇a)的左临界点(iu-K〇A)作为灰度阈值;否则，执行子步骤3 ? 4;
[0018] 3.4判断口(1|似，〇、)与当前高斯分量口(1|11(1，〇(1)的高频区间是否有交集印(1|11 (1， ~)初始化为P(x | iu，〇a)右侧相邻高斯分量;若有交集，令p(x |此，％)右侧相邻高斯分量为当 前高斯分量P (x | ，〇a)，执行本子步骤;否则，以p (x |此，〇a)右临界点ya+K〇a为灰度阈值；
[0019] 上述，若卫星影像为全色影像，系数K取2.5~3;若卫星影像为多光谱影像，系数K 取2~2.5;
[0020] S4采用灰度阈值分割云区域。
[0021] 步骤S1中所述的预处理灰度直方图以消除干扰，进一步包括：
[0022] 舍弃位于灰度直方图左右两端各占总像素数a%的像素，得灰度直方图主体区间， a%在0.01 %~1 %范围取值；
[0023] 对灰度直方图主体区间进行曲线平滑。
[0024] 子步骤2.1中，若获取的峰值点数大于5,仅保留灰度值最大的5个峰值点，则高斯 分量数为5。
[0025] 子步骤2.4具体为：
[0026]构建最大似然模型
，采用最大似然函数对各高斯 分量的权重、均值、标准差进行迭代，当满足收敛条件
吋，停止迭代；
[0027]其中，代表最大似然函数;<+1)、乂广>、〇，+1)分别为第t+1次迭代下第m个高斯 分量的权重、均值、标准差；/〇、分别为第t次迭代下第m个高斯分量的权重、均 值、标准差;卩;|表示第t次迭代下的第m个高斯分量。
[0028]本发明方法弥补了传统纹理分析法、同态滤波法、多光谱综合法的不足，可用于卫 星影像的筛选与质量控制、云区自动提取、含云影像快速修补等工作。
[0029] 本发明具有如下优点和有益效果：
[0030] (1)检测精度高，误判较少；同时适用于含云、无云影像。
[0031] (2)检测基于单波段影像直方图，不受限于卫星光谱范围，同时适用于全色、多光 谱影像。
[0032] (3)无需辅助信息和人工干预，计算速度快，可自动化处理。
[0033] (4)可扩展性强，当条件允许时，可将本发明方法作为初始判别准则，使用复合决 策的方式，取得更精准的云检测结果。
【附图说明】
[0034]图1为卫星影像及其对应的灰度直方图，其中，图（a)为理想的含云卫星影像及其 对应的灰度直方图，图（b)为含云量较少且云层较薄的含云卫星影像及其对应的灰度直方 图，图（c)为无云卫星影像及其对应的灰度直方图；
[0035]图2为本发明灰度阈值确定的流程图；
[0036]图3中，图（a)为本发明实施例的资源三号卫星全色含云影像，图（b)为云检测结 果，图（c)为影像灰度直方图及高斯混合模型；
[0037]图4中，图（a)为本发明实施例的资源三号卫星全色无云影像，图（b)为云检测结 果，图（c)为影像灰度直方图及高斯混合模型。
【具体实施方式】
[0038] 为了更好地理解本发明技术方案，下面将结合附图和实施例对本发明做进一步详 细说明。
[0039] 本【具体实施方式】在Microsoft Visual C++6.0环境下，由C语言编程实现，整个过 程可实现自动化处理。
[0040] 步骤1，直方图统计及预处理。
[0041] 读取卫星影像，统计卫星影像的原始灰度直方图，灰度直方图可反映图像灰度分 布的统计特征。由于地物本身辐射特征的复杂性，以及传感器CCD在某些情况下的异常(极 端）响应，直方图会受到"噪声"的干扰。因此在对灰度直方图进行分析前，需要对其进行必 要的预处理。
[0042]本步骤的【具体实施方式】如下：
[0043]步骤1.1，读取卫星影像，并统计卫星影像的原始灰度直方图。
[0044] 步骤1.2,舍弃位于原始灰度直方图左右两端各占总数a%的像素，得到灰度直方 图主体区间，将主体区间的灰阶范围即有效灰阶范围，记为[Xmin，Xm ax]，Xmin和)(_分别为像 素灰度值的极低值和极高值，a %在0.01 %~1 %范围取值。
[0045] 步骤1.3,采用大小为5X5的窗口对灰度直方图主体区间进行曲线平滑。
[0046] 下面步骤均基于预处理后的灰度直方图进行。
[0047]步骤2，建立高斯混合模型。
[0048]对于一幅理想的单波段含云卫星影像，其灰度直方图应该具有双峰，见图1(a)。但 对于含云量较少且云层较薄的卫星影像，灰度直方图在高亮区域并不会形成明显的波峰， 而是呈平缓延伸状，见图1(b)。另一方面，无云卫星影像没有表示云层的波峰，见图1(c)。综 合考虑，灰度直方图分布可看作由多个简单分布混合形成的多峰形态，此时灰度直方图函 数h(x)可以用高斯混合模型GMM(x)近似表示：
(1) (2)
[0051 ]式（1)~（2)中：
[0052] X表示灰阶；
[0053] m表示高斯分量的序号，M为高斯分量数；
[0054] p ( X I ，〇m)表示第m个高斯分量，( X I ，〇m)的权重；
[0055] ^和~分别表示第m个高斯分量的均值和标准差。
[0056] 设影像总像素数为N，则Tm满足： M
[0057] N 二 Yjm C3 ) m-\
[0058]本具体实施中，使用局部窗口法获取高斯混合模型的初始模型参数，所述的模型 参数包括高斯分量数M以及各高斯分量的权重均值、标准差^，并通过最大似然法(EM) 对模型参数进行迭代计算，最终得到稳定的高斯混合模型。
[0059]本步骤的【具体实施方式】如下：
[0060]步骤2.1，设置局部窗口，使用局部最大值法获取灰度直方图的峰值点，峰值点数 即高斯分量数M。若所获取峰值点数大于5，仅保留灰度值最大的5个峰值点，则高斯分量数M 为5。记峰值点横坐标为Pi，i = l，2，...M。
[0061 ]步骤2.2，使用最大类间方差法，获取相邻峰值点间的谷值点；同时，灰度直方图的 最左和最右两端点也被标记为谷值点。记谷值点横坐标为Vj，j = 1，2，. . .M+1。谷值点Vk、Vk+1 分别为位于峰值点Pk两侧，k=l，2,. . .M。
[0062] 步骤2.3,计算第m个高斯分量的权重、均值、标准差的初始值^、of :
[0063] ⑷
[0064] ^'〇, = S [/?'-v)(-v-/;.r)"] ( 5 ) K;;(m
[0065] zf) =艺./?(." (6)
[0066]式(4)~(6)中，Pm表示第m个高斯分量峰值点横坐标;h(x)表示灰度直方图，即灰 阶x对应的像素数。
[0067]步骤2.4，构建最大似然模型，迭代拟合模型参数，如下：
[0072] 式(7)~（10)中：
[0073] 代表最大似然函数，t表示迭代次数。
[0074]当满足下述收敛条件，停止迭代：
[0078]各模型参数趋于收敛后，高斯混合模型趋于稳定。
[0079]步骤3，高斯混合模型特征分析。
[0080] 在可见光及近红外波谱范围内，云层反射率较高，而多数地物反射率要低很多，表 现在卫星图像上，即云像元和地物像元的辐射特性具有明显差异。定量地分析各分量之间 的关系，主要包括步骤：
[0081] 步骤3.1，灰度直方图有效灰阶范围[Xmin，Xmax]内，表征地物的高斯分量主要集中 在灰度相对较低的区域，即灰度直方图的左侧;表征云层的高斯分量则主要集中在灰度直 方图的右侧。因此，将有效灰阶范围[X min，Xmax]的中间点￥^作为云层和地物间的初始分类 (14) 阈值：
[0083] 步骤3.2，表征地物和表征云层的高斯分量之间波峰相距较远，二者相对独立；而 当表征不同地物的高斯分量之间波峰距离较近，二者高度相干。所以，将[ym-l.30m，y m+l.3 0m]看作第m个高斯分量的高频区间。由高斯概率密度分布函数可知，任意像素分布在此高 频区间内的概率约为80 %。
[0084]步骤3.3，理想的灰度阈值应位于高斯混合模型表征地物的高斯分量右侧、表征云 层高斯分量左侧。基于此，将ym±K〇m#记为第m个高斯分量左右两侧的临界点。其中K值由影 像类型确定。
[0085]步骤4,灰度阈值确定，见图2。
[0086]步骤4.1，找出权重最大的高斯分量，记为p(x |iU，〇A)。若p(x |iU，〇A)的均值似> Yini，则认为该权重最大的高斯分量表征云，卫星影像含云量较大，以高斯分量P(X |iU，〇A)的 左临界点(似-Kw)作为灰度阈值;否则，认为该权重最大的高斯分量表征地物，卫星影像少 云或无云，此时执行步骤4.2。
[0087]步骤4 ? 2，判断表征地物的高斯分量p (x | iu，〇a)与当前高斯分量p (x | ，〇d)的高频 区间是否有交集，当前高斯分量P(X|Wc〇初始化为高斯分量P(x|iU，〇0的右侧相邻高斯 分量;若有交集，令P(X|Wc〇右侧相邻高斯分量为当前高斯分量口（1|如，〇。），重复步骤 4 ? 2;若无交集，以p (x | ，〇a)的右临界点ya+K〇a作为灰度阈值。
[0088]步骤5,基于灰度阈值分割的云区提取。
[0089]如果卫星影像是全色影像，系数K取2.5~3,得到灰度阈值Y;采用灰度阈值Y进行 分割得到卫星影像的高亮区域，即被判定为云区域。
[0090]如果卫星影像是多光谱影像，系数K取2~2.5,分别计算红、绿、蓝三个波段的灰度 阈值￥[?、￥(；、￥8，将灰度值同时高于灰度阈值￥[?、￥( ；、￥8的像素判断为初始云区像素。
【主权项】
1. 一种基于高斯混合模型的卫星影像自动云检测方法，其特征是，包括： Sl统计卫星影像的原始灰度直方图，预处理灰度直方图以消除干扰； S2对灰度直方图构建高斯混合模型，本步骤进一步包括： 2.1采用高斯混合模型表示灰度直方图函数； 2.2采用局部最大值法获得灰度直方图的峰值点，峰值点数即高斯分量数M; 2.2使用最大类间方差法获取相邻峰值点间的谷值点；同时，将灰度直方图的两端点也 标记为谷值点，谷值点横坐标记为Vj，j = 1，2，. . .M+1; 2.3初始化高斯混合模型的模型参数，获得各高斯分量的权重、均值、标准差的初始值)，Pm表示第m个 高斯分量峰值点横坐标;h(x)表示灰度直方图中灰阶X对应的像素数； 2.4采用最大似然法对模型参数进行迭代拟合，获得高斯混合模型； S3确定分割云区域的灰度阈值，本步骤进一步包括： ?V1相《??々卜捆辰mr直方图的灰阶范围[xmin，X max ]确定云层和地物间的初始分类阈 3.2设置第m个高斯分量的高频区间[1.3〇m，ym+1.3〇m]，y m、om分别为第m个高斯分量 的均值和标准差，从子步骤2.4获得的高斯混合模型获得； 3 · 3找出权重最大的高斯分量，记为p(x I μλ，σλ);若ρ(χ I μλ，σλ)的均值μλ>Υ?η?，以ρ(χ I μλ，〇λ)的左临界点(μλ-Κ〇λ)作为灰度阈值;否则，执行子步骤3.4; 3.4判断口(1|以、，似)与当前高斯分量口(1卜(1，〇(1)的高频区间是否有交集4(1|以 (1，〇(1)初 始化为Ρ(χ|μλ，〇0右侧相邻高斯分量;若有交集，令Ρ(χ|μα，~)右侧相邻高斯分量为当前高 斯分量P (XI ，ο。），执行本子步骤;否则，以P (XI ，ο。)右临界点μα+Κσα为灰度阈值； 上述，若卫星影像为全色影像，系数K取2.5~3;若卫星影像为多光谱影像，系数K取2~ 2.5； S4米用灰度阈值分割云区域。2. 如权利要求1所述的基于高斯混合模型的卫星影像自动云检测方法，其特征是： 所述的预处理灰度直方图以消除干扰，进一步包括： 舍弃位于灰度直方图左右两端各占总像素数a%的像素，得灰度直方图主体区间，a% 在0.01%~1%范围取值； 对灰度直方图主体区间进行曲线平滑。3. 如权利要求1所述的基于高斯混合模型的卫星影像自动云检测方法，其特征是： 子步骤2.1中，若获取的峰值点数大于5,仅保留灰度值最大的5个峰值点，则高斯分量 数为5。4. 如权利要求1所述的基于高斯混合模型的卫星影像自动云检测方法，其特征是： 子步骤2.4具体为： 构建最大似然模？1最大似然函数对各高斯分量 的权重、均值、标准差进行迭代，当满足收敛条件；止迭代； 其中，尺Γ代表最大似然函数；、σ/+η分别为第t+Ι次迭代下第m个高斯分 量的权重、均值、标准差；:rf、/4卩、分别为第t次迭代下第m个高斯分量的权重、均值、 标准差；"^|/0,4卩；)表示第〇欠迭代下的第 111个高斯分量。
【文档编号】G06T7/00GK105894520SQ201610259961
【公开日】2016年8月24日
【申请日】2016年4月25日
【发明人】康飞, 康一飞, 孙明伟, 王树根
【申请人】武汉大学