留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于高斯混合模型法的国产高分辨率卫星影像云检测

康一飞 潘励 孙明伟 陈奇 王越

康一飞, 潘励, 孙明伟, 陈奇, 王越. 基于高斯混合模型法的国产高分辨率卫星影像云检测[J]. 武汉大学学报 ● 信息科学版, 2017, 42(6): 782-788. doi: 10.13203/j.whugis20140875
引用本文: 康一飞, 潘励, 孙明伟, 陈奇, 王越. 基于高斯混合模型法的国产高分辨率卫星影像云检测[J]. 武汉大学学报 ● 信息科学版, 2017, 42(6): 782-788. doi: 10.13203/j.whugis20140875
KANG Yifei, PAN Li, SUN Mingwei, CHEN Qi, WANG Yue. Gaussian Mixture Model Based Cloud Detection for Chinese High Resolution Satellite Imagery[J]. Geomatics and Information Science of Wuhan University, 2017, 42(6): 782-788. doi: 10.13203/j.whugis20140875
Citation: KANG Yifei, PAN Li, SUN Mingwei, CHEN Qi, WANG Yue. Gaussian Mixture Model Based Cloud Detection for Chinese High Resolution Satellite Imagery[J]. Geomatics and Information Science of Wuhan University, 2017, 42(6): 782-788. doi: 10.13203/j.whugis20140875

基于高斯混合模型法的国产高分辨率卫星影像云检测

doi: 10.13203/j.whugis20140875
基金项目: 

国家973计划 2012CB719900

国家自然科学基金 41301519

详细信息
    作者简介:

    康一飞, 博士生, 主要从事数字摄影测量研究。2217707@163.com

    通讯作者: 潘励, 博士, 教授。panli@whu.edu.cn
  • 中图分类号: P236

Gaussian Mixture Model Based Cloud Detection for Chinese High Resolution Satellite Imagery

Funds: 

The Major State Basic Research Development National of China (973 Program) 2012CB719900

the National Natural Science Foundation of China 41301519

More Information
    Author Bio:

    KANG Yifei, PhD candidate, specializes in digital photogrammetry. E-mail:2217707@163.com

    Corresponding author: PAN Li, PhD, professor. E-mail:panli@whu.edu.cn
图(10) / 表(2)
计量
  • 文章访问数:  1246
  • HTML全文浏览量:  86
  • PDF下载量:  429
  • 被引次数: 0
出版历程
  • 收稿日期:  2015-10-22
  • 刊出日期:  2017-06-05

基于高斯混合模型法的国产高分辨率卫星影像云检测

doi: 10.13203/j.whugis20140875
    基金项目:

    国家973计划 2012CB719900

    国家自然科学基金 41301519

    作者简介:

    康一飞, 博士生, 主要从事数字摄影测量研究。2217707@163.com

    通讯作者: 潘励, 博士, 教授。panli@whu.edu.cn
  • 中图分类号: P236

摘要: 针对资源三号、高分一号等国产高分辨率卫星影像波段少、光谱范围受限的特点,提出一种通过高斯混合模型拟合影像灰度直方图从而自动获取灰度阈值的云检测算法。首先由影像灰度直方图自适应地获取高斯混合模型初始拟合参数,然后依据期望最大原则对初始参数进行调整,最后根据拟合模型中各高斯分量的分布特点,自动确定该波段影像中云与晴空之间的灰度阈值。实验表明,该算法不受限于卫星光谱范围,同时适用于含云和无云影像,检测精度较高,且不需要辅助信息和人工干预,可满足自动化生产的需要。

English Abstract

康一飞, 潘励, 孙明伟, 陈奇, 王越. 基于高斯混合模型法的国产高分辨率卫星影像云检测[J]. 武汉大学学报 ● 信息科学版, 2017, 42(6): 782-788. doi: 10.13203/j.whugis20140875
引用本文: 康一飞, 潘励, 孙明伟, 陈奇, 王越. 基于高斯混合模型法的国产高分辨率卫星影像云检测[J]. 武汉大学学报 ● 信息科学版, 2017, 42(6): 782-788. doi: 10.13203/j.whugis20140875
KANG Yifei, PAN Li, SUN Mingwei, CHEN Qi, WANG Yue. Gaussian Mixture Model Based Cloud Detection for Chinese High Resolution Satellite Imagery[J]. Geomatics and Information Science of Wuhan University, 2017, 42(6): 782-788. doi: 10.13203/j.whugis20140875
Citation: KANG Yifei, PAN Li, SUN Mingwei, CHEN Qi, WANG Yue. Gaussian Mixture Model Based Cloud Detection for Chinese High Resolution Satellite Imagery[J]. Geomatics and Information Science of Wuhan University, 2017, 42(6): 782-788. doi: 10.13203/j.whugis20140875
  • 以资源三号、高分一号为代表,国产高分辨率对地观测卫星影像广泛应用于资源调查、防灾减灾等领域。然而云层的存在不仅覆盖了地表信息,也给影像的配准、融合等处理造成诸多影响。

    当前较为常用的云检测算法主要有纹理分析法[1, 2]、同态滤波法[3, 4]、多光谱综合法[5, 6]等。然而现有算法并不能完全适用于国产高分辨率卫星影像,例如:① 纹理分析法基于统计方法提取云区和非云区的空间特性,可有效识别大片层云,但难以识别高分辨率卫星影像上常见的纹理较强的卷积云;② 同态滤波法不擅长处理厚云影像,且滤波的过程中会丢失一些有用信息,降低高分辨率卫星影像的质量;③ 多光谱综合法能够利用云和晴空反射率的差异,但它要求传感器配有多个热红外波段,探测波长需涵盖水或二氧化碳的吸收带[1];而国产高分辨率卫星影像一般仅含有蓝、绿、红、近红外共4个波段,光谱范围也一般为0.43~0.90 μm[7],多光谱综合法的优势难以体现。

    综合考虑上述情况,本文提出了一种基于高斯混合模型的云检测算法,并以资源三号、高分一号卫星影像为实验数据,分析算法的可行性和云检测精度。

    • 灰度直方图反映图像灰度分布的统计特征(即横坐标代表影像灰度,无单位)。由于地物本身辐射特征的复杂性以及传感器CCD在某些情况下的异常响应,直方图会受到噪声的干扰。因此在对直方图进行分析前,需要对其进行预处理。

      首先,舍弃位于直方图左右两端各占总数0.01%的像素,得到直方图的主体区间,设为[gmin, gmax];然后对该区间内的直方图用一维窗口 $\left[ {\frac{1}{{10}}, \frac{2}{{10}}, \frac{4}{{10}}, \frac{2}{{10}}, \frac{1}{{10}}} \right]$ 进行平滑,消除随机噪声的干扰。

      对于一幅理想的单波段含云影像,其直方图应该具有双峰,如图 1(a)所示,此时可用最大类间方差算法[8]确定双峰间的谷值,作为云与地物之间的阈值。对于含云量较少且云层较薄的影像,直方图在高亮区域并不会形成明显的波峰,而是呈平缓延伸状,如图 1(b)所示。无云影像则没有表示云层的波峰,如图 1(c)所示。

      图  1  资源三号卫星全色影像及其直方图

      Figure 1.  ZY-3 Panchromatic Image and Its Histogram

      进一步,对于一幅包含多种地物的卫星影像,其直方图分布可看作由多个简单分布混合形成的多峰形态[12],此时影像直方图函数h(x)可以用一个高斯混合模型GMM(x)来近似表示:

      $$ \left\{ \begin{array}{l} GMM\left( x \right) = \sum\limits_{m = 1}^M {\left[ {{\tau _m}G\left( {x\left| {{\theta _m}} \right.} \right)} \right]} \\ G\left( {x\left| {{\theta _m}} \right.} \right) = \frac{1}{{\sqrt {2\pi } {\sigma _m}}}\exp \left[ { - \frac{{{{\left( {x - {\mu _m}} \right)}^2}}}{{2\sigma _m^2}}} \right] \end{array} \right. $$ (1)

      式中,M为混合模型中高斯分量的数目;G(x|θm)为第m个分量对应的波形函数,m = 1, 2, …,Mμmσmτm分别为第m个高斯分量的均值、标准差和权重。设影像总像素数为N,则权重τm满足:

      $$ N = \sum\limits_{m = 1}^M {{\tau _m}} $$ (2)
    • 获取高斯混合模型的参数是直方图拟合的关键,可使用迭代方法建立最大似然方程,采用期望最大(expectation maximization, EM)算法[9-12]对参数进行估计。

      EM算法需要首先设置M以及各高斯分量的权重初值τm(0)、均值初值μm(0)、方差初值σm2(0)。一般情况下,M可根据经验给定一个较大的值[10],也可由贝叶斯信息准则得出[11];初始参数τm(0)μm(0)σm2(0)则由人工调查和目视解译[12]得到的样本集来确定。

      为了避免人工参与,本文设置一定大小的局部窗口,使用局部最大值法获取直方图峰值点;将过小的峰值点排除,剩余较大的峰值点的个数即为M;然后在相邻峰值之间,使用局部最小值法获取谷值点。

      Pm为第m个高斯分量峰值点横坐标,VmVm+1分别为第m个高斯分量左右两侧谷值点横坐标,则初始参数τm(0)μm(0)σm2(0)可分别表示为:

      $$ \left\{ \begin{array}{l} \tau _m^{\left( 0 \right)} = \sum\limits_{x \in L} {h\left( x \right)} \\ \mu _m^{\left( 0 \right)} = {P_m}\\ \sigma _m^{2\left( 0 \right)} = \sum\limits_{x \in L} {\left[ {h\left( x \right){{\left( {x - \mu _m^{\left( 0 \right)}} \right)}^2}} \right]/\tau _m^{\left( 0 \right)}} \end{array} \right. $$ (3)

      式中,L=[Vm, Vm+1]。

    • 根据上一步获取的初始值,计算样本n属于第m类的后验概率Qmn,其公式为:

      $$ {Q_{mn}} = \frac{{{\tau _m}}}{N}G\left( {x\left| {{\theta _m}} \right.} \right) $$ (4)

      标准化为Rmn

      $$ {R_{mn}} = \frac{{{Q_{mn}}}}{{{Q_n}}} = \frac{{{\tau _m}G\left( {x\left| {{\theta _m}} \right.} \right)}}{{\sum\limits_{m = 1}^M {\left[ {{\tau _m}G\left( {x\left| {{\theta _m}} \right.} \right)} \right]} }}] $$ (5)

      最大化式(5),得到新的参数值:

      $$ \left\{ {\begin{array}{*{20}{l}} {{{\tilde \tau }_m} = \sum\limits_{n = 1}^N {{R_{mn}}} }\\ {{{\tilde \mu }_m} = \frac{{\sum\limits_{n = 1}^N {\left( {{R_{mn}}{x_n}} \right)} }}{{N{{\tilde \tau }_m}}}}\\ {\tilde \sigma _m^2 = \frac{{\sum\limits_{n = 1}^N {\left[ {{R_{mn}}{{\left( {{x_n} - {{\tilde \mu }_m}} \right)}^2}} \right]} }}{{N{{\tilde \tau }_m}}}} \end{array}} \right. $$ (6)

      式中,xn为样本n的观测值;N为样本个数。针对数字影像计算中可将影像全部像素都选为样本的特殊情况,用xn代表像素n的灰度值,N代表总像素数,可看作直方图函数h(x)在区间L内的积分:

      $$ N = \sum\limits_{x \in L} {h\left( x \right)} $$ (7)

      式中,L=[gmin, gmax]。基于式(7),将式(6) 改写为:

      $$ \left\{ {\begin{array}{*{20}{l}} {{{\tilde \tau }_m} = \sum\limits_{x \in L} {\left[ {h\left( x \right){R_{mn}}} \right]} }\\ {{{\tilde \mu }_m} = \frac{{\sum\limits_{x \in L} {\left[ {h\left( x \right){R_{mn}}x} \right]} }}{{\sum\limits_{x \in L} {\left[ {h\left( x \right){R_{mn}}} \right]} }}}\\ {\tilde \sigma _m^2 = \frac{{\sum\limits_{x \in L} {\left[ {h\left( x \right){R_{mn}}{{\left( {x - {{\tilde \mu }_m}} \right)}^2}} \right]} }}{{\sum\limits_{x \in L} {\left[ {h\left( x \right){R_{mn}}} \right]} }}} \end{array}} \right. $$ (8)

      式中,L=[gmin, gmax]。

      迭代计算式(5)、式(8),直至 ${{\tilde \tau }_m}$ 、 ${{{\tilde \mu }_m}}$ 、 $\tilde \sigma _m^2$ 收敛,此时高斯混合模型趋于稳定,各参数得到确定。

    • 在可见光及近红外波谱范围内,云层反射率较高,而多数地物反射率要低很多,表现在卫星图像上,即云像元和晴空像元的辐射特性具有明显差异[13, 14]。观察分析大量国产高分辨率卫星影像的灰度直方图和高斯混合模型,可以得到两条统计规律:

      (1) 在影像直方图有效灰阶范围[gmin, gmax]内,表征地物的高斯分量主要集中在灰度相对较低的区域,即直方图的左侧,表征云层的高斯分量则主要集中在直方图的右侧。

      (2) 分别表征地物和云层的两个高斯分量之间,波峰相距较远且相对独立;表征不同地物的高斯分量之间波峰距离较近且高度相干。

      基于规律(1),将[gmin, gmax]区间的中间点Yini作为云层和地物之间的初始分类阈值:

      $$ {Y_{{\rm{ini}}}} = \frac{1}{2}\left( {{g_{\max }} + {g_{\min }}} \right) $$ (9)

      基于规律(2),将[(μm-k1σm), (μm+k1σm)]看作高斯分量G(x|θm)的高频区间,若两个高斯分量的高频区间存在交集,则说明两个波相干性较强,反之相干性较弱。其中方差系数k1的大小在一定程度上影响检测结果,当k1过大时,会造成薄云与地物高频区间相交,使算法对薄云不敏感;当k1过小时,又容易造成复杂实验区不同地物之间高频区间相互分离,使算法对高亮度地物产生误判。本文实验中,由经验设置k1=1.5;大量测试证明,此取值适用于资源三号、高分一号卫星全色和多光谱传感器影像。

    • 首先找出高斯混合模型中权重最大的高斯分量,记为G(x|θλ),若其均值μλ>Yini,则认为最大分量表征云,影像含云量较大;反之,则认为最大分量表征地物,影像少云或者无云。

      理想的阈值应位于高斯混合模型表征地物分量的右侧、表征云层分量的左侧。将μm±k2σm看作高斯分量G(x|θm)左右两侧的临界点,其中k2也是方差系数。若影像含云量较大,则以表征云的高斯分量G(x|θλ)的左临界点μλ-k2σλ作为阈值;若影像少云或者无云,则以表征地物的高斯分量G(x|θλ)为参照,通过判断相邻高斯分量的高频区间是否相交,将位于低亮度区域的数个相干性较强的高斯分量都看作地物。最后以地物类别中最右侧的高斯分量G(x|θα)的右临界点μα+k2σα作为阈值。阈值选取流程如图 2所示。

      图  2  阈值选取流程

      Figure 2.  Flowchart of Threshold Selection

      大量实验证实,对国产高分辨率卫星,如果是全色影像,取系数k2=3,得到阈值Y后即可通过阈值分割的方式获取云区;如果是多光谱影像,取系数k2=2.5,分别计算红、绿、蓝三个波段的阈值YRYGYB,当像素灰度同时满足三个阈值的限制时,方可判定为云。

    • 选取4幅国产高分辨率卫星影像作为实验数据。其中影像1、影像2为资源三号全色影像,分辨率为2.1 m,大小为24 530×24 575像素;影像3、影像4为高分一号多光谱影像,分辨率为8 m,大小为4 548×4 500像素。影像1中间区域被厚云覆盖,厚云边缘和影像左上角存在少量薄云,晴空部分为山地;影像2无云,主要地物为山地和城镇;影像3存在大量卷积云,晴空部分为河水和城镇;影像4无云,主要地物为农田和城镇。

      使用文献[15]中提到的带限定条件(对比实验需要部分人工干预)的Otsu算法进行对比实验,分别获取本文算法和对比算法计算得到的云检测阈值,根据阈值对原始影像做二值分割,目视评价检测效果。同时以人工解译勾选出的云区作为真值,利用混淆矩阵计算4幅影像检测结果的总体精度和Kappa系数,定量评价云检测精度。

      图 3~图 6可以看出,对于4组测试数据,使用本文算法计算的阈值,取得了良好的云检测效果,且适用于无云影像。由定量评价参数表 1表 2可看出,相对于带限定条件的Otsu算法,本文算法检测精度更高,总体精度在95%以上,Kappa系数在0.9以上。

      图  3  资源三号含云影像检测效果

      Figure 3.  Detection Results of ZY-3 Cloudy Images

      图  4  资源三号无云影像检测效果

      Figure 4.  Detection Results of ZY-3 Cloudless Images

      图  5  高分一号含云影像检测效果

      Figure 5.  Detection Results of GF-1 Cloudy Images

      图  6  高分一号无云影像检测效果

      Figure 6.  Detection Results of GF-1 Cloudless Images

      表 1  资源三号影像检测评价参数

      Table 1.  Parameters of ZY-3 Images

      参数 本文算法 Otsu算法
      含云 总体精度 0.978 541 0.946 691
      Kappa系数 0.918 724 0.732 822
      无云 总体精度 0.990 397 0.959 307
      Kappa系数 / /

      表 2  高分一号影像检测评价参数

      Table 2.  Parameters of GF-1 Images

      参数 本文算法 Otsu算法
      含云 总体精度 0.986 191 0.986 177
      Kappa系数 0.933 681 0.933 612
      无云 总体精度 0.991 470 0.960 120
      Kappa系数 / /

      图 7图 8图 9图 10中,各自的(a)、(b)、(c)子图均表示本文实验阈值选取的过程,其中绿色曲线表示高斯分量,黄色曲线表示最大高斯分量,蓝色曲线表示叠加后的高斯混合模型,初始阈值由红色虚线划出,精确阈值由红点标出。可以看出,本文算法获取的高斯混合模型均对直方图做出了较为精确的拟合,同时保证了阈值位于直方图中表征地物高斯分量的右侧。以上4图的(d)子图中的橙色点表示带限定条件的Otsu算法阈值。对于图 7(d)图 9(d),初次Otsu阈值可取得较好的分割效果,不需调整;对于图 8(d)图 10(d),需要以初次Otsu阈值(左侧标记点)为限定条件,对初始阈值右侧区间再次执行分割,得到带限定条件的Otsu阈值(右侧标记点)。

      图  7  资源三号含云影像直方图及阈值

      Figure 7.  Histogram and Threshold of ZY-3 Cloudy Images

      图  8  资源三号无云影像直方图及阈值

      Figure 8.  Histogram and Threshold of ZY-3 Cloudless Images

      图  9  高分一号含云影像直方图及阈值(红波段)

      Figure 9.  Histogram and Threshold of GF-1 Cloudy Images

      图  10  高分一号无云影像直方图及阈值(红波段)

      Figure 10.  Histogram and Threshold of GF-1 Cloudless Images

    • 针对国产高分辨率卫星影像的云检测,本文提出的基于高斯混合模型的阈值选取算法在一定程度上弥补了传统纹理分析法、同态滤波法、多光谱综合法的不足,可用于卫星影像的筛选与质量控制、云区自动提取、含云影像快速修补等工作。主要有:① 检测精度高,误判较少,同时适用于含云、无云影像;② 检测基于单波段影像直方图,不受限于卫星光谱范围,同时适用于全色、多光谱影像;③ 无需辅助信息和人工干预,计算速度快,可自动化处理。

      阈值法也有其固有的局限性,主要表现为容易对亮度较高的建筑物、裸地、积雪产生误判。当条件允许时,可将本文算法作为初始判别准则,使用复合决策的方式,取得更精准的云检测结果[15-18],例如使用形态学开操作消除对建筑物的误判,并平滑云区轮廓;使用归一化水指数、热红外传感器信息、卫星成像信息等来区分云和积雪。

参考文献 (18)

目录

    /

    返回文章
    返回