CN105894520A - 一种基于高斯混合模型的卫星影像自动云检测方法 - Google Patents

Abstract

本发明公开了一种基于高斯混合模型的卫星影像自动云检测方法，包括：S1统计卫星影像的原始灰度直方图，预处理灰度直方图以消除干扰；S2对灰度直方图构建高斯混合模型；S3 确定分割云区域的灰度阈值；S4 采用灰度阈值分割云区域。本发明方法弥补了传统纹理分析法、同态滤波法、多光谱综合法的不足，可用于卫星影像的筛选与质量控制、云区自动提取、含云影像快速修补等工作。本发明基于单波段影像直方图，不受限于卫星光谱范围，同时适用于全色、多光谱影像；检测精度高，无需辅助信息和人工干预，计算速度快，可自动化处理。

Description

一种基于高斯混合模型的卫星影像自动云检测方法

技术领域

本发明属于遥感图像处理技术领域，特别涉及一种基于高斯混合模型的卫星影像自动云检测方法。

背景技术

以WorldView系列、Pleiades系列、中国“高分”系列卫星为代表，高分辨率对地观测技术迅速发展，卫星重访周期不断变小、卫星影像分辨率和幅宽不断提高、卫星影像的数据量急剧增加，广泛应用于资源调查、防灾减灾等领域。然而在卫星影像数据处理过程中，云的存在是一个不可忽视的问题，它不仅对地面场景形成遮挡，还在一定程度上改变了影像的光谱特性，给影像匹配、匀色和镶嵌工作造成诸多不利影响。如何精准、快速地检测卫星影像中的云层，从而更好地发挥卫星影像的利用价值，是目前摄影测量与遥感领域的热点课题。

当前较为常用的云检测方法主要有纹理分析法、同态滤波法、多光谱综合法等。但是，考虑到高分辨率对地观测卫星影像波段少、光谱范围窄(0.45～0.90微米)、数据量大的特点，现有云检测方法并不能完全适用于高分辨率对地观测卫星影像。纹理分析法基于统计方法提取云区和非云区的空间特性，可有效识别大片层云，但难以识别纹理较强的卷积云。同态滤波法不擅长处理厚云影像，且滤波的过程中会丢失一些有用信息，降低影像产品质量。多光谱综合法利用云和晴空反射率的差异，但它要求传感器配有多个热红外波段，探测波长需涵盖水或二氧化碳的吸收带，而对地观测卫星影像一般仅含有蓝、绿、红、近红外4个波段，光谱范围一般为0.43～0.90μm，多光谱综合法的优势难以体现。

发明内容

本发明的目的是针对高分辨率对地观测卫星影像，提供一种高效精确的基于高斯混合模型的卫星影像自动云检测方法。

为达到上述目的，本发明采用如下技术方案：

一种基于高斯混合模型的卫星影像自动云检测方法，包括步骤：

S1统计卫星影像的原始灰度直方图，预处理灰度直方图以消除干扰；

S2对灰度直方图构建高斯混合模型，本步骤进一步包括：

2.1采用高斯混合模型表示灰度直方图函数；

2.2采用局部最大值法获得灰度直方图的峰值点，峰值点数即高斯分量数M；

2.2使用最大类间方差法获取相邻峰值点间的谷值点；同时，将灰度直方图的两端点也标记为谷值点，谷值点横坐标记为V_j，j＝1,2,...M+1；

2.3初始化高斯混合模型的模型参数，获得各高斯分量的权重、均值、标准差的初始值其中，P_m表示第m个高斯分量峰值点横坐标；h(x)表示灰度直方图中灰阶x对应的像素数；

2.4采用最大似然法对模型参数进行迭代拟合，获得高斯混合模型；

S3确定分割云区域的灰度阈值，本步骤进一步包括：

3.1根据预处理后灰度直方图的灰阶范围[X_min,X_max]确定云层和地物间的初始分类阈值

3.2设置第m个高斯分量的高频区间[μ_m-1.3σ_m,μ_m+1.3σ_m]，μ_m、σ_m分别为第m个高斯分量的均值和标准差，从子步骤2.4获得的高斯混合模型获得；

3.3找出权重最大的高斯分量，记为p(x|μ_λ,σ_λ)；若p(x|μ_λ,σ_λ)的均值μ_λ＞Y_ini，以p(x|μ_λ,σ_λ)的左临界点(μ_λ-Kσ_λ)作为灰度阈值；否则，执行子步骤3.4；

3.4判断p(x|μ_λ,σ_λ)与当前高斯分量p(x|μ_α,σ_α)的高频区间是否有交集，p(x|μ_α,σ_α)初始化为p(x|μ_λ,σ_λ)右侧相邻高斯分量；若有交集，令p(x|μ_α,σ_α)右侧相邻高斯分量为当前高斯分量p(x|μ_α,σ_α)，执行本子步骤；否则，以p(x|μ_α,σ_α)右临界点μ_α+Kσ_α为灰度阈值；

上述，若卫星影像为全色影像，系数K取2.5～3；若卫星影像为多光谱影像，系数K取2～2.5；

S4采用灰度阈值分割云区域。

步骤S1中所述的预处理灰度直方图以消除干扰，进一步包括：

舍弃位于灰度直方图左右两端各占总像素数a％的像素，得灰度直方图主体区间，a％在0.01％～1％范围取值；

对灰度直方图主体区间进行曲线平滑。

子步骤2.1中，若获取的峰值点数大于5，仅保留灰度值最大的5个峰值点，则高斯分量数为5。

子步骤2.4具体为：

构建最大似然模型采用最大似然函数对各高斯分量的权重、均值、标准差进行迭代，当满足收敛条件时，停止迭代；

其中，代表最大似然函数；σ_m ^(t+1)分别为第t+1次迭代下第m个高斯分量的权重、均值、标准差；分别为第t次迭代下第m个高斯分量的权重、均值、标准差；表示第t次迭代下的第m个高斯分量。

本发明方法弥补了传统纹理分析法、同态滤波法、多光谱综合法的不足，可用于卫星影像的筛选与质量控制、云区自动提取、含云影像快速修补等工作。

本发明具有如下优点和有益效果：

(1)检测精度高，误判较少；同时适用于含云、无云影像。

(2)检测基于单波段影像直方图，不受限于卫星光谱范围，同时适用于全色、多光谱影像。

(3)无需辅助信息和人工干预，计算速度快，可自动化处理。

(4)可扩展性强，当条件允许时，可将本发明方法作为初始判别准则，使用复合决策的方式，取得更精准的云检测结果。

附图说明

图1为卫星影像及其对应的灰度直方图，其中，图(a)为理想的含云卫星影像及其对应的灰度直方图，图(b)为含云量较少且云层较薄的含云卫星影像及其对应的灰度直方图，图(c)为无云卫星影像及其对应的灰度直方图；

图2为本发明灰度阈值确定的流程图；

图3中，图(a)为本发明实施例的资源三号卫星全色含云影像，图(b)为云检测结果，图(c)为影像灰度直方图及高斯混合模型；

图4中，图(a)为本发明实施例的资源三号卫星全色无云影像，图(b)为云检测结果，图(c)为影像灰度直方图及高斯混合模型。

具体实施方式

为了更好地理解本发明技术方案，下面将结合附图和实施例对本发明做进一步详细说明。

本具体实施方式在Microsoft Visual C++6.0环境下，由C语言编程实现，整个过程可实现自动化处理。

步骤1，直方图统计及预处理。

读取卫星影像，统计卫星影像的原始灰度直方图，灰度直方图可反映图像灰度分布的统计特征。由于地物本身辐射特征的复杂性，以及传感器CCD在某些情况下的异常(极端)响应，直方图会受到“噪声”的干扰。因此在对灰度直方图进行分析前，需要对其进行必要的预处理。

本步骤的具体实施方式如下：

步骤1.1，读取卫星影像，并统计卫星影像的原始灰度直方图。

步骤1.2，舍弃位于原始灰度直方图左右两端各占总数a％的像素，得到灰度直方图主体区间，将主体区间的灰阶范围即有效灰阶范围，记为[X_min,X_max]，X_min和X_max分别为像素灰度值的极低值和极高值，a％在0.01％～1％范围取值。

步骤1.3，采用大小为5×5的窗口对灰度直方图主体区间进行曲线平滑。

下面步骤均基于预处理后的灰度直方图进行。

步骤2，建立高斯混合模型。

对于一幅理想的单波段含云卫星影像，其灰度直方图应该具有双峰，见图1(a)。但对于含云量较少且云层较薄的卫星影像，灰度直方图在高亮区域并不会形成明显的波峰，而是呈平缓延伸状，见图1(b)。另一方面，无云卫星影像没有表示云层的波峰，见图1(c)。综合考虑，灰度直方图分布可看作由多个简单分布混合形成的多峰形态，此时灰度直方图函数h(x)可以用高斯混合模型GMM(x)近似表示：

$GMM\left(x\right)=\underset{m=1}{\overset{M}{\Σ}}\[{\mathrm{\τ}}_{m}p\left(x\right|{\mathrm{\μ}}_m,{\mathrm{\σ}}_m)\]---\left(1\right)$

$p\left(x\right|{\mathrm{\μ}}_m,{\mathrm{\σ}}_m)=\frac{1}{{\mathrm{\σ}}_m\sqrt{2\mathrm{\π}}}{e}^{-\frac{{(x-{\mathrm{\μ}}_m)}^2}{2{\mathrm{\σ}}_m^2}}---\left(2\right)$

式(1)～(2)中：

X表示灰阶；

m表示高斯分量的序号，M为高斯分量数；

p(x|μ_m,σ_m)表示第m个高斯分量，τ_m为p(x|μ_m,σ_m)的权重；

μ_m和σ_m分别表示第m个高斯分量的均值和标准差。

设影像总像素数为N，则τ_m满足：

$N=\underset{m=1}{\overset{M}{\Σ}}{\mathrm{\τ}}_m---\left(3\right)$

本具体实施中，使用局部窗口法获取高斯混合模型的初始模型参数，所述的模型参数包括高斯分量数M以及各高斯分量的权重τ_m、均值μ_m、标准差σ_m，并通过最大似然法(EM)对模型参数进行迭代计算，最终得到稳定的高斯混合模型。

本步骤的具体实施方式如下：

步骤2.1，设置局部窗口，使用局部最大值法获取灰度直方图的峰值点，峰值点数即高斯分量数M。若所获取峰值点数大于5，仅保留灰度值最大的5个峰值点，则高斯分量数M为5。记峰值点横坐标为P_i，i＝1,2,...M。

步骤2.2，使用最大类间方差法，获取相邻峰值点间的谷值点；同时，灰度直方图的最左和最右两端点也被标记为谷值点。记谷值点横坐标为V_j，j＝1,2,...M+1。谷值点V_k、V_k+1分别为位于峰值点P_k两侧，k＝1,2,...M。

步骤2.3，计算第m个高斯分量的权重、均值、标准差的初始值

${\mathrm{\μ}}_m^{\left(0\right)}=P_m---\left(4\right)$

${\mathrm{\σ}}_m^{2\left(0\right)}=\underset{x=V_m}{\overset{V_{m+1}}{\Σ}}\[h\left(x\right){(x-{\mathrm{\μ}}_m^{\left(0\right)})}^2\]---\left(5\right)$

${\mathrm{\τ}}_m^{\left(0\right)}=\underset{x=V_m}{\overset{V_{m+1}}{\Σ}}h\left(x\right)---\left(6\right)$

式(4)～(6)中，P_m表示第m个高斯分量峰值点横坐标；h(x)表示灰度直方图，即灰阶x对应的像素数。

步骤2.4，构建最大似然模型，迭代拟合模型参数，如下：

$R_m^{\left(t\right)}=\frac{\[{\mathrm{\τ}}_m^{\left(t\right)}p\left(x\right|{\mathrm{\μ}}_m^{\left(t\right)},{\mathrm{\σ}}_m^{\left(t\right)})\]}{\underset{m=1}{\overset{M}{\Σ}}\[{\mathrm{\τ}}_m^{\left(t\right)}p\left(x\right|{\mathrm{\μ}}_m^{\left(t\right)},{\mathrm{\σ}}_m^{\left(t\right)})\]}---\left(7\right)$

${\mathrm{\μ}}_m^{(t+1)}=\frac{\underset{x=X_{\min }}{\overset{X_{\max }}{\Σ}}\[xh\left(x\right)R_m^{\left(t\right)}\]}{\underset{x=X_{\min }}{\overset{X_{\max }}{\Σ}}\[h\left(x\right)R_m^{\left(t\right)}\]}---\left(8\right)$

${\mathrm{\σ}}_m^{2(t+1)}=\frac{\underset{x=X_{\min }}{\overset{X_{\max }}{\Σ}}\[(x-{\mathrm{\τ}}_m^{(t+1)}h(x)R_m^{\left(t\right)}\]}{\underset{x=X_{\min }}{\overset{X_{\max }}{\Σ}}\[h\left(x\right)R_m^{\left(t\right)}\]}---\left(9\right)$

${\mathrm{\τ}}_m^{(t+1)}=\underset{x=X_{\min }}{\overset{X_{\max }}{\Σ}}\[h\left(x\right)R_m^{\left(t\right)}\]---\left(10\right)$

式(7)～(10)中：

代表最大似然函数，t表示迭代次数。

当满足下述收敛条件，停止迭代：

${\mathrm{\&mu;}}_m^{(t+1)}-{\mathrm{\&mu;}}_m^{\left(t\right)}<\frac{{\mathrm{\&mu;}}_m^{\left(t\right)}}{100}---\left(11\right)$

${\mathrm{\&sigma;}}_m^{2(t+1)}-{\mathrm{\&sigma;}}_m^{2\left(t\right)}<\frac{{\mathrm{\&sigma;}}_m^{2\left(t\right)}}{100}---\left(12\right)$

${\mathrm{\&tau;}}_m^{(t+1)}-{\mathrm{\&tau;}}_m^{\left(t\right)}<\frac{{\mathrm{\&tau;}}_m^{\left(t\right)}}{100}---\left(13\right)$

各模型参数趋于收敛后，高斯混合模型趋于稳定。

步骤3，高斯混合模型特征分析。

在可见光及近红外波谱范围内，云层反射率较高，而多数地物反射率要低很多，表现在卫星图像上，即云像元和地物像元的辐射特性具有明显差异。定量地分析各分量之间的关系，主要包括步骤：

步骤3.1，灰度直方图有效灰阶范围[X_min,X_max]内，表征地物的高斯分量主要集中在灰度相对较低的区域，即灰度直方图的左侧；表征云层的高斯分量则主要集中在灰度直方图的右侧。因此，将有效灰阶范围[X_min,X_max]的中间点Y_ini作为云层和地物间的初始分类阈值：

$Y_{ini}=\frac{1}{2}(X_{max}+X_{\min })---\left(14\right)$

步骤3.2，表征地物和表征云层的高斯分量之间波峰相距较远，二者相对独立；而当表征不同地物的高斯分量之间波峰距离较近，二者高度相干。所以，将[μ_m-1.3σ_m,μ_m+1.3σ_m]看作第m个高斯分量的高频区间。由高斯概率密度分布函数可知，任意像素分布在此高频区间内的概率约为80％。

步骤3.3，理想的灰度阈值应位于高斯混合模型表征地物的高斯分量右侧、表征云层高斯分量左侧。基于此，将μ_m±Kσ_m标记为第m个高斯分量左右两侧的临界点。其中K值由影像类型确定。

步骤4，灰度阈值确定，见图2。

步骤4.1，找出权重最大的高斯分量，记为p(x|μ_λ,σ_λ)。若p(x|μ_λ,σ_λ)的均值μ_λ＞Y_ini，则认为该权重最大的高斯分量表征云，卫星影像含云量较大，以高斯分量p(x|μ_λ,σ_λ)的左临界点(μ_λ-Kσ_λ)作为灰度阈值；否则，认为该权重最大的高斯分量表征地物，卫星影像少云或无云，此时执行步骤4.2。

步骤4.2，判断表征地物的高斯分量p(x|μ_λ,σ_λ)与当前高斯分量p(x|μ_α,σ_α)的高频区间是否有交集，当前高斯分量p(x|μ_α,σ_α)初始化为高斯分量p(x|μ_λ,σ_λ)的右侧相邻高斯分量；若有交集，令p(x|μ_α,σ_α)右侧相邻高斯分量为当前高斯分量p(x|μ_α,σ_α)，重复步骤4.2；若无交集，以p(x|μ_α,σ_α)的右临界点μ_α+Kσ_α作为灰度阈值。

步骤5，基于灰度阈值分割的云区提取。

如果卫星影像是全色影像，系数K取2.5～3，得到灰度阈值Y；采用灰度阈值Y进行分割得到卫星影像的高亮区域，即被判定为云区域。

如果卫星影像是多光谱影像，系数K取2～2.5，分别计算红、绿、蓝三个波段的灰度阈值Y_R、Y_G、Y_B，将灰度值同时高于灰度阈值Y_R、Y_G、Y_B的像素判断为初始云区像素。

Claims (4)

1.一种基于高斯混合模型的卫星影像自动云检测方法，其特征是，包括：

S1统计卫星影像的原始灰度直方图，预处理灰度直方图以消除干扰；

S2对灰度直方图构建高斯混合模型，本步骤进一步包括：

2.1采用高斯混合模型表示灰度直方图函数；

2.2采用局部最大值法获得灰度直方图的峰值点，峰值点数即高斯分量数M；

2.2使用最大类间方差法获取相邻峰值点间的谷值点；同时，将灰度直方图的两端点也标记为谷值点，谷值点横坐标记为V_j，j＝1,2,...M+1；

2.3初始化高斯混合模型的模型参数，获得各高斯分量的权重、均值、标准差的初始值其中，P_m表示第m个高斯分量峰值点横坐标；h(x)表示灰度直方图中灰阶x对应的像素数；

2.4采用最大似然法对模型参数进行迭代拟合，获得高斯混合模型；

S3确定分割云区域的灰度阈值，本步骤进一步包括：

3.1根据预处理后灰度直方图的灰阶范围[X_min,X_max]确定云层和地物间的初始分类阈值

3.2设置第m个高斯分量的高频区间[μ_m-1.3σ_m,μ_m+1.3σ_m]，μ_m、σ_m分别为第m个高斯分量的均值和标准差，从子步骤2.4获得的高斯混合模型获得；

3.3找出权重最大的高斯分量，记为p(x|μ_λ,σ_λ)；若p(x|μ_λ,σ_λ)的均值μ_λ＞Y_ini，以p(x|μ_λ,σ_λ)的左临界点(μ_λ-Kσ_λ)作为灰度阈值；否则，执行子步骤3.4；

3.4判断p(x|μ_λ,σ_λ)与当前高斯分量p(x|μ_α,σ_α)的高频区间是否有交集，p(x|μ_α,σ_α)初始化为p(x|μ_λ,σ_λ)右侧相邻高斯分量；若有交集，令p(x|μ_α,σ_α)右侧相邻高斯分量为当前高斯分量p(x|μ_α,σ_α)，执行本子步骤；否则，以p(x|μ_α,σ_α)右临界点μ_α+Kσ_α为灰度阈值；

上述，若卫星影像为全色影像，系数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具体为：

构建最大似然模型采用最大似然函数对各高斯分量的权重、均值、标准差进行迭代，当满足收敛条件时，停止迭代；

其中，代表最大似然函数；分别为第t+1次迭代下第m个高斯分量的权重、均值、标准差；分别为第t次迭代下第m个高斯分量的权重、均值、标准差；表示第t次迭代下的第m个高斯分量。