CN104616253B - 一种利用独立成分分析技术的光学遥感图像薄云去除方法 - Google Patents

Abstract

本发明公开了一种利用独立成分分析技术的光学遥感图像薄云去除方法，它包括数据准备，大气校正，提取二值图，ICA变换，提取和处理云成分和ICA逆变换这些步骤。本发明利用光学遥感数据单幅图像自身数据，利用了独立成分分析技术提取云独立成分并处理，获得去云的结果。该方法摆脱了多数据或多传感器去云方法对数据要求苛刻等诸多束缚，一定程度上解决了其他单幅遥感影像去云方法中的不足。在多光谱光学遥感数据中，该方法去除了薄云对地物光谱特征的影响，提高了光学遥感数据的使用质量和图像的应用能力。本发明与现有的去薄云的方法相比，相对于单幅数据方法要进行时域变换、频域变换和频率域滤波等一系列复杂处理，克服了现有方法对数据要求苛刻、处理方法复杂等缺点，本发明所用的独立成分分析技术成熟，操作高效、简单，有很大的实际应用价值。

Description

一种利用独立成分分析技术的光学遥感图像薄云去除方法

技术领域

本发明属于光学遥感图像去云技术领域，具体涉及一种利用独立成分分析技术的光学遥感图像薄云去除方法。

背景技术

去云是光学遥感图像处理中一个重要的组成部分。光学遥感图像通过卫星传感器获取光谱信息，由于大气的存在，太阳能通过大气时会发生散射和吸收作用，这使卫星获取的遥感图像不能准确的捕捉地物的信息，从而影响遥感图像的质量。若大气中存在云层，将进一步影响光学遥感图像获取地物信息。

随着大气校正技术的发展和完善，大气校正算法和软件正被广泛的使用，包括：ACORN–Atmospheric CORrection Now(InSpec,2002),ATREM–the ATmospheric REMovalprogram[Center for the Study of Earth from Space(CSES),University ofColorado)],FLAASH–Fast Line–of–sight Atmospheric Analysis of SpectralHypercubes(Research Systems,Inc.,2003)等。然而由于云层的复杂特性和时空易变性，一定程度上阻碍了遥感图像去云技术的发展。本发明专注于有效的去薄云的方法。

目前，主要的去薄云的方法可以分为以下几种：

基于辐射传输特性的方法，该方法在辐射传输的过程中考虑薄云的影响，利用辐射传输过程中，不同光谱范围对薄云的反应特征，建立传输关系，从而得到去云结果，具体方法参见“Gao,B.C.,P.Yang,W.Han,R.Li,W.J.Wiscombe.2002.An algorithm usingvisible and1.38-μm channels to retrieve cirrus cloud reflectances fromaircraft and satellite data.IEEE Transactions on Geoscience and RemoteSensing,vol.40,No.8,pp.1659-1668.”此类方法要求有单独检测薄云的波段，对遥感光学数据要求高。

基于滤波的方法，该类方法主要是讲图像在频率域上去除低频成分，不可避免的去除了一些有用的信息，并且在选取截止频率时凭经验选择，不可靠。具体的方法参见“张波，季民河，沈琪.2011.基于小波变换的高分辨率快鸟遥感图像薄云去除[J].遥感信息,03:38–43.”

Haze Optimized Transformation(HOT)，在两个波段的散点图中(红光和蓝光波段为该方法推荐使用的波段)，HOT定义了一条“晴空线”，像素点到“晴空线”的距离反应了云或者尘埃的影响程度。具体方法参见“Y.Zhang,B.Guindon J.Cihlar.2002.An imagetransform to characterize and compensate for spatial variations in thin cloudcontamination of Landsat images,Remote Sensing of Environment,vol.82,No.2,pp.173-187.”该方法无法应用于具有高反射值或高亮度值的区域，如雪或者干沙地，也无法适用具有低反射值或低亮度值的区域，如开阔的水体表面。

基于缨帽变换的方法，该方法用于去云是缨帽变换的衍生和扩展，并且缨帽变换的第四分量被认为主要贡献是云和大气尘埃，具体方法参见“Richter,R.1996.Atmospheric correction of satellite data with haze removal including ahaze/clear transition region.Computer and Geosciences,22,675–681.”但是，值得注意的是设计该方法的初衷是用于评估土地和植被对地表辐射的贡献(如：亮度，绿度和湿度)。因此，用该方法的第四分量去去云值得商榷。

发明内容

本发明的目的是为了解决上述已有去薄云技术中存在的问题，提出了一种利用独立成分分析技术的光学遥感图像薄云去除方法，从而恢复在多光谱光学遥感数据中薄云覆盖区域处地物光谱信息，提高光学遥感图像的质量和应用能力。

为了方便描述本发明的内容，首先作以下定义：

定义1、大气校正

光学遥感的大气校正是去除遥感数据中的大气效应，获取地表反射率的过程。大气校正主要包括两部分：大气参数估计和地表反射率反演。对于水平均匀的大气和朗伯体地面，地表反射率r_λ是通过以下公式得到。

其中，L_λ为表观光谱辐亮度为，L_p是大气路径辐射，S是大气的半球反照率，F₀乘以π后是在大气顶部垂直于太阳光束入射的太阳能通量密度，τ(μ_s)和τ(μ_v)是太阳到地面和地表到传感器的总透过率。μ_s和μ_v是太阳角和观测角的余弦值。详见文献“定量遥感”，梁顺林等编著，科学出版社，2012。

定义2、大气校正软件模块

大气校正软件模块可以通过用户提供基础大气特征信息或者特定的大气吸收波段计算特定时间的大气散射和吸收特性以获取地表反射率。广泛用于大气校正的软件和算法包括：ACORN–Atmospheric CORrection Now(InSpec,2002)，ATCOR–the ATmosphericCORrection program(Thiemann and Hermann 2002)，ATREM–the ATmospheric REMovalprogram[Center for the Study of Earth from Space(CSES),University ofColorado)],FLAASH–Fast Line–of–sight Atmospheric Analysis of SpectralHypercubes(Research Systems,Inc.,2003)等。这些大气校正软件模块通常需要用户提供：

·遥感图像经纬度信息，

·遥感图像获取日期和时刻，

·遥感图像的海拔高度,

·卫星传感器的高度，

·大气模型(如：中纬度-夏季，中纬度-冬季，热带)，

·辐射校正的辐亮度数据(如：数据单位必须是W·cm^-2·sr^-1·μm^-1)，

·传感器特定的波段信息(如：波段的半高全宽)，

·遥感图像获取时的大气能见度等。

定义3、阈值分割

图像阈值化分割是一种传统的最常用的图像分割方法，因其实现简单、计算量小、性能较稳定而成为图像分割中最基本和应用最广泛的分割技术。基本原理是：通过设定不同的特征阈值，把图像象素点分为若干类。

常用的特征包括有直接来自原始图像的灰度或彩色特征或由原始灰度或彩色值变换得到的特征。

按照一定的准则找到特征值，将图像分割为两个部分，分割后的图像取：类别1＝0(黑)，类别2＝1(白)，即为我们通常所说的图像二值化。

定义4、卷云波段数据

卷云波段指波长在1370纳米左右，宽度在20纳米的红外波段。由于水对该波段的光谱能量有强吸收作用，几乎能量不能到达地面，所以能够被该波段检测到的信息都是被大气以上的物质反射的能量。基于这个特性，该波段主要用于检测平均高度在6千米以上的卷云的信息。

定义5、独立成分分析和ICA变换

给定随机变量的一组观测(x₁(t),x₂(t)...,x_n(t))，其向量表示为x，t为时间或者样本标号。假设它们由独立成分(s₁(t),s₂(t)...,s_n(t))组成，其向量表示为s，线性混合而产生：

或以向量的形式表示为：

x＝As

式中，A是一个未知矩阵。在仅能观测到x_i(t)的情况下，独立成分分析就要同时估计出矩阵A和x_i(t)。

独立成分分析寻找独立成分的方法可以大致分为以下几种：

极大非高斯性的估计方法，极大似然估计方法，极小互信息的估计方法，基于张量的估计方法等。

通过这些方法估计矩阵A和x_i(t)，并用公式s＝A^-1x计算独立成分s的过程称为ICA变换。

定义6、ICA逆变换

将独立成分s利用公式x＝As计算得到随机变量的一组观察x的过程称为ICA逆变换。

定义7、相关系数

相关系数是用以反映变量之间相关关系密切程度的统计指标。相关系数是按协差方法计算，其计算公式为：

其中，x_i和y_i表示两类样本第i个元素，i为正整数，1≤i≤n，n为样本个数；和表示两类样本的平均值。

本发明提供了一种利用独立成分分析技术的光学遥感图像薄云去除方法，该方法的步骤如下：

步骤1：数据准备

本发明提供的数据包括：多光谱光学遥感图像数据Q，它有n个图像，包括i个可见光波段图像B_v，j个近红外波段图像B_r，k个短波红外波段图像B_w和1个卷云波段图像B_c；其中i，j，k，n为正整数，且i+j+k+1＝n。

步骤2：大气校正

将步骤1中的多光谱光学遥感图像数据Q通过经典的大气校正软件模块进行传统的大气校正处理，获得大气校正结果Q_a。

步骤3：提取二值图矩阵(“0”代表无云，“1”代表有云)

将步骤1中的卷云波段图像B_c进行传统的阈值分割方法处理得到二值图D_c，图像由“0”和“1”两类元素组成，图像中0代表无云区域像素点的像素值，1代表有云区域像素点的像素值。

步骤4：ICA变换

将步骤2中得到的大气校正结果Q_a进行ICA变换得到转换结果IC，它有n个独立成分组成，分别是IC₁，IC₂，…,IC_n；

步骤5：提取云成分并处理

将步骤4得到的n个独立成分，分别与卷云波段图像B_c计算相关系数，得到n个相关系数；然后找出最大相关系数对应的独立成分，该独立成分即为云独立成分，记为IC_cloud。

利用步骤3得到的云二值图D_c和云独立成分IC_cloud，通过公式

M＝mean(IC_cloud*(！D_c))

计算得到云独立成分无云区域的平均值M，其中符号“mean(.)”为求平均值运算，“！”表示取反运算。最后将云独立成分IC_cloud的每一个像素值用平均值M替换得到新的云独立成分IC’_cloud。

将去云后的云独立成分IC’_cloud和其他n-1个独立成分重新组合形成新的独立成分IC’。

步骤6：ICA逆变换

将步骤5得到的IC’进行传统的ICA逆变换，得到n个新的多光谱遥感波段，即为去云结果。

本发明提供的一种利用独立成分分析技术的光学遥感图像薄云去除方法，该方法充分利用了单幅多光谱遥感图像数据的特征，与传统的多时相或者多传感器的去云方法相比，对数据要求相对较宽松，与传统的单幅多光谱遥感数据去云方法相比，处理手段更加简单，更容易实现。

附图说明

下面以Landsat-8遥感图像数据为例，用附图详细阐述本发明的实现过程，使本发明的目的和效果变得更加明显。

图1为本发明的流程示意图。

图2为Landsat-8多光谱遥感图像数据8个波段图像。(a)-(d)表示深蓝，蓝光，绿光和红光4个可见光波段图像；(e)表示近红外波段图像；(f)和(g)表示短波红外波段图像；(h)表示卷云波段图像。

图3为Landsat-8多光谱遥感图像数据进行大气校正处理后的结果。(a)-(d)表示4个可见光波段图像；(e)表示近红外波段图像；(f)和(g)表示短波红外波段图像；(h)表示卷云波段图像。

图4为卷云波段进行阈值分割方法得到的二值云图(0–无云，1表示有云)。

图5为大气校正后结果进行的ICA变换的结果。(a)-(h)分布表示IC₁-IC₈。

图6为卷云波段图像分布与IC₁-IC₈的相关系数。

图7为去云结果。(a)-(d)表示深蓝，蓝光，绿光和红光4个可见光波段图像；(e)表示近红外波段图像；(f)和(g)表示短波红外波段图像。

具体实施方法

为了对本发明的技术特征、目的和效果有更加清楚的理解，现说明本发明的具体实施方式。

参考图1处理流程图，本发明以Landsat-8遥感图像数据为例，具体阐述这种利用独立成分分析技术的光学遥感图像薄云去除方法。

步骤1：数据准备

本发明提供的数据包括：多光谱光学遥感图像数据Q，它有8个图像，包括4个可见光波段图像B_v，1个近红外波段图像B_r，2个短波红外波段图像B_w和1个卷云波段图像B_c。每个波段图像如图2所示。

步骤2：大气校正

将步骤1中的多光谱光学遥感图像数据Q通过经典的大气校正软件模块进行传统的大气校正处理，获得大气校正结果Q_a。每个波段数据如图3所示，(a)-(d)分别为深蓝，蓝光，绿光和红光4个可见光波段图像，(e)为近红外波段，(f),(g)为两个短波红外波段，(h)为卷云波段。

步骤3：提取二值图矩阵(“0”代表无云，“1”代表有云)

将步骤1中的卷云波段图像B_c进行传统的阈值分割方法处理得到二值图D_c，图像由“0”和“1”两类元素组成，图像中“0”代表无云区域像素点的像素值，“1”代表有云区域像素点的像素值。二值图如图4所示。

步骤4：ICA变换

将步骤2中得到的大气校正结果Q_a进行ICA变换得到转换结果IC，它有8个独立成分组成，分别是IC₁，IC₂，…,IC₈。8个独立成分如图5所示。

步骤5：提取云成分并处理

将步骤4得到的8个独立成分，分别与卷云波段图像B_c计算相关系数，得到8个相关系数，如图6中数据所示；然后找出最大相关系数对应的独立成分，本次实例中IC₁与卷云的相关系数最大，独立成分IC₁即为云独立成分，记为IC_cloud。

利用步骤3得到的云二值图D_c和云独立成分IC_cloud，通过公式

M＝mean(IC_cloud*(！Dc))

计算得到云独立成分无云区域的平均值M，其中符号“mean(.)”为求平均值运算，“！”表示取反运算。最后将云独立成分IC_cloud的每一个像素值用平均值M替换得到新的云独立成分IC’_cloud。

将去云后的云独立成分IC’_cloud和其他n-1个独立成分重新组合形成新的独立成分IC’。

步骤6：ICA逆变换

将步骤5得到的IC’进行传统的ICA逆变换，得到7个多光谱遥感图像去云的结果，如图7所示，(a)-(d)分别为深蓝，蓝光，绿光和红光4个可见光波段图像，(e)为近红外波段，(f),(g)为两个短波红外波段。

Claims (1)

1.一种利用独立成分分析技术的光学遥感图像薄云去除方法，其特征是它包括如下步骤：

步骤1：数据准备

提供的数据包括：多光谱光学遥感图像数据Q，它有n个图像，包括i个可见光波段图像B_v，j个近红外波段图像B_r，k个短波红外波段图像B_w和1个卷云波段图像B_c；其中i，j，k，n为正整数，且i+j+k+1＝n；

步骤2：大气校正

将步骤1中的多光谱光学遥感图像数据Q通过经典的大气校正软件模块进行传统的大气校正处理，获得大气校正结果Q_a；

步骤3：提取二值图矩阵(“0”代表无云，“1”代表有云)

将步骤1中的卷云波段图像B_c进行传统的阈值分割方法处理得到二值图D_c，图像由“0”和“1”两类元素组成，图像中0代表无云区域像素点的像素值，1代表有云区域像素点的像素值；

步骤4：ICA变换

将步骤2中得到的大气校正结果Q_a进行ICA变换得到转换结果IC，它有n个独立成分组成，分别是IC₁，IC₂，…,IC_n；

步骤5：提取云成分并处理

将步骤4得到的n个独立成分，分别与卷云波段图像B_c计算相关系数，得到n个相关系数，然后找出最大相关系数对应的独立成分，该独立成分即为云独立成分，记为IC_cloud；

利用步骤3得到的云二值图D_c和云独立成分IC_cloud，通过公式

M＝mean(IC_cloud*(！D_c))

计算得到云独立成分无云区域的平均值M，其中符号“mean(.)”为求平均值运算，“！”表示取反运算，最后将云独立成分IC_cloud的每一个像素值用平均值M替换得到新的云独立成分IC’_cloud；

将去云后的云独立成分IC’_cloud和其他n-1个独立成分重新组合形成新的独立成分IC’；

步骤6：ICA逆变换

将步骤5得到的IC’进行传统的ICA逆变换，得到n个新的多光谱遥感波段，即为去云结果。