CN104778668A - 基于可见光波段光谱统计特征的光学遥感图像薄云去除方法 - Google Patents

Abstract

基于可见光光谱统计特征的光学遥感图像薄云去除方法，其特征在于，包括如下步骤：步骤1：数据准备步骤；2：表观反射率转换；步骤3：大气校正；步骤4：区域分割；步骤5：统计特征参数提取；步骤6：利用去云辐射传输方程求解去云结果。本发明提出了一种基于可见光光谱统计特征的光学遥感图像薄云去除方法，从而使得光学遥感图像能够反映真实的地物光谱信息，提高多光谱光学遥感数据的使用质量、适用范围及应用能力。

Description

基于可见光波段光谱统计特征的光学遥感图像薄云去除方法

技术领域

本发明属于光学遥感图像去云技术领域，具体涉及一种基于可见光波段光谱统计特征的光学遥感图像薄云去除方法。

背景技术

由于大气对太阳辐射具有选择性的吸收及散射作用，使得光学遥感卫星的传感器接收到的地物信息包含了大气对其光谱的影响，减弱了遥感图像反映地物真实信息的能力，尤其在大气中存在云层时，将进一步影响光学遥感图像的质量，极大的限制了光学遥感图像的应用范围和实用价值。

随着大气校正技术的不断发展和完善，现在有许多较为成熟的大气校正算法和软件正在被广泛使用，包括：ACORN–AtmosphericCORrection Now(InSpec,2002)，ATCOR–the ATmosphericCORrection program(Thiemann and Hermann 2002)，ATREM–theATmospheric REMoval program[Center for the Study of Earth fromSpace(CSES),University of Colorado)],FLAASH–FastLine–of–sight Atmospheric Analysis of Spectral Hypercubes(Research Systems,Inc.,2003)。然而，由于云层的时空易变性及其下垫面的复杂性，如何去除云层对遥感图像的影响一直是研究的热点和难点。当云层较薄时，经由地物反射和散射的太阳辐射能够一定程度的穿透云层并由传感器接受并记录下来，这种部分穿透性使得复原受到薄云影响的地物信息成为可能。

本发明或本领域中所指薄云，是指虽然有云层阻挡，但卫星仍能透过云层获取云层以下的地表信息，只要云层不能全部覆盖地表信息，即可称为薄云。只有当云层较薄时，来自地物的辐射才能部分穿透云层到达光学传感器并被记录下来，这种部分穿透性是从单幅光学图像上恢复云层下的地物信息的必要条件。现有的薄云去除方法大致可分为两类：

一类是基于信号处理和图像复原的方法，例如同态滤波、IHS变换、小波变换等，具体实施方法参见[1]赵忠明,朱重光，1996，遥感图像中薄云的去除方法，环境遥感,11(3):195-199。[2]Souza,M.P.,E.B.Pereira,F.R.Martins,R.C.Chagas,and W.S.FreitasJr.,2003.The cloud cover fraction obtained from a ground CCDcamera and its effect on a radiative transfer model,EGS-AGU-EUG Joint Assembly,1:2601。[3]Guo,X.,2008.The studyof removal of thin cloud in CBERS-02CCD based on the wavelettransform,Science of Surveying and Mapping,33(S1).4415p.这类算法处理过程简便且在特定的区域取得了一定的效果，然而这些方法的物理意义不明确，同时在去除薄云影响的同时丢失了一定的地物信息，精度不能得到保证。

另一类是基于辐射传输模型的方法，例如S.Liang,H.Fang,andM.Chen,2001.“Atmospheric correction of Landsat ETM+landsurface imagery-part I:methods,”IEEE Transactions ongeoscience and remote sensing,39(11):2490-2498.当有条件获取研究区域的相关物理参数时，这类算法的精度较高。然而在实际应用中如何准确的获取相关实测参数限制了该类方法的适用范围。

Gao提出了一种简化的辐射传输方程：

${\mathrm{\ρ}}_i^{*}={\mathrm{\ρ}}_{\mathrm{ci}}+T_{\mathrm{ci}}*{\mathrm{\ρ}}_i$

其中为传感器记录的i波段的地表反射率；ρ_ci为薄云在i波段的反射率；T_ci薄云对i波段反射率的双向透过率(包含直射和散射)；ρ_i为i波段的地物表观反射率。根据Gao的研究，T_ci近乎于常量且通常大于0.9，ρ_ci则需要卫星搭载的传感器具有专门的薄云探测波段来进行计算。值得一提的是，现有的大部分光学卫星并没有相应的薄云探测波段，因此该方法的应用范围受到了很大限制。

发明内容

本发明的目的是为了解决上述已有的薄云去除方法的不足，提出了一种基于可见光光谱统计特征的光学遥感图像薄云去除方法，从而使得光学遥感图像能够反映真实的地物光谱信息，提高多光谱光学遥感数据的使用质量、适用范围及应用能力。

本发明所述基于可见光光谱统计特征的光学遥感图像薄云去除方法，包括如下步骤：

步骤1：数据准备

采集多光谱光学遥感图像数据及其头文件Hdr，它有I个可见光波段图层；从头文件Hdr中获得近红外波段反射率数据Nir；各波段修正因子和各波段偏移量太阳高度角θ_SE；其中i、I均为正整数且i＝1，2，3，...，I，i表示不同波段，I表示采集的波段总数；

步骤2：表观反射率转换

将步骤1中的多光谱光学数据代入公式中，分别求得I个波段的表观反射率其中i＝1，2，3，...，I；

步骤3：大气校正

将步骤2中各个波段的表观反射率矩阵进行大气校正，分别获得I个波段的地表反射率矩阵对每一下标i，在该步骤中计算的地表反射率矩阵中下标不同于的记做即j＝1，2，3，...，I且j≠i；

步骤4：区域分割

从步骤1中的多光谱光学遥感图像数据中选取任意一块无云的晴空子区域，记该区域内的蓝光波段反射率为ρ_b，该区域内的红光波段反射率为ρ_r；ρ_b、ρ_r均为矩阵形式，对ρ_b、ρ_r利用方程ρ_r＝t₁*ρ_b+t₂进行线性拟合，计算出t₁、t₂；通过公式计算得到HOT，其中为中蓝光波段在该晴空子区域上的子集，为中红光波段在该晴空子区域上的子集；

对HOT通过阈值分割方法得到二值图Mask₁；由计算得到NDWI，其中Nir为步骤1中得到的近红外波段反射率数据，为中的绿光波段反射率；将NDWI通过阈值分割方法划分得到二值图Mask₂；

步骤5：统计特征参数提取

将步骤3中的大气校正结果代入公式及公式分别得到晴空区域的反射率Clear_i和薄云覆盖的有水区域的反射率Cloud_i，其中符号“！”表示取反运算，i＝1，2，3，...，I；

在该步骤中计算的晴空区域的反射率中下标不同于Clear_i的记做Clear_j，即j＝1，2，3，...，I且j≠i；对Clear_i和Clear_j进行线性拟合，得到拟合线性方程为Clear_j＝k_ij*Clear_i+b_ij；对Cloud_i和Cloud_j进行线性拟合，对应的线性方程为Cloud_j＝m_ij*Cloud_i+n_ij；其中i分别取1，2，3，...，I，共得到I组(k_ij,b_ij,m_ij,n_ij)；

步骤6：利用去云辐射传输方程求解去云结果

将步骤5中提取的I组统计特征参数(k_ij,b_ij,m_ij,n_ij)及步骤3中的大气校正结果代入去云辐射传输方程

其中i分别取1，2，3，...，I；j＝1，2，3，...，I且j≠i；T_ci为薄云对i波段反射率的双向透过率；即可分别求解得到对应i波段的去云结果ρ_i，其中i＝1，2，3，...，I。

具体的，步骤4中所述蓝光波段反射率、红光波段反射率、绿光波段反射率分别为波长在0.450–0.515um、0.630–0.680um、0.525–0.600um范围内。

具体的，所述步骤4中利用阈值分割方法得到二值图Mask₁、Mask₂过程中的阈值均为零。

采用本发明所述基于可见光光谱统计特征的光学遥感图像薄云去除方法，本发明提出了一种基于可见光光谱统计特征的光学遥感图像薄云去除方法，从而使得光学遥感图像能够反映真实的地物光谱信息，提高多光谱光学遥感数据的使用质量、适用范围及应用能力。

附图说明

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

图2为采用Landsat-8多光谱遥感图像数据的可见光波段图像，(a)、(b)、(c)、(d)分别表示波长为波段1、波段2、波段3、波段4的图像；从上至下分别为薄云覆盖遥感图像数据及采用本发明所述方法处理后的图像。

图3为统计特征参数提取结果，反映可见光波段间统计特征关系的散点图及其拟合线；图3中(a)、(b)、(c)、(d)的横、纵坐标分别表示晴空区下波段1和波段2、波段2和波段3、波段3和波段1、波段4和波段3的反射率，晴空区下波段i的反射率即所述步骤5中所求的Clear_i。图3所示仅表示各点的线性拟合程度，各点的具体坐标数值不影响线性拟合度的判断。

图2及图3中波段1为0.433–0.453um波长范围，波段2为0.450–0.515um波长范围，波段3为0.525–0.600um波长范围，波段4为0.630–0.680um波长范围。

具体实施方式

下面结合附图，对本发明的具体实施方式作进一步的详细说明。

本发明所述基于可见光光谱统计特征的光学遥感图像薄云去除方法，包括如下步骤：

步骤1：数据准备

采集多光谱光学遥感图像数据及其头文件Hdr、近红外波段反射率数据Nir，以上数据通常可由卫星遥感数据中直接获取。多光谱光学遥感图像数据有I个可见光波段图层；各波段修正因子和各波段偏移量太阳高度角θ_SE可以从头文件Hdr中得到；其中i、I均为正整数且i＝1，2，3，...，I，i表示不同波段，I表示采集的波段总数；

步骤2：表观反射率转换

将步骤1中的多光谱光学数据代入公式中，分别求得I个波段的表观反射率矩阵其中i＝1，2，3，...，I；

步骤3：大气校正

将步骤2中各个波段的表观反射率矩阵进行大气校正，分别获得I个波段的地表反射率矩阵对每一下标i，在该步骤中计算得到与下标不同于地表反射率矩阵的记做关联矩阵即j＝1，2，3，...，I且j≠i；对每一有(I-1)个例如对于i＝1，则j＝2、3、4。

此处地表反射率是指地面反射辐射量与入射辐射量之比，表征地面对太阳辐射的吸收和反射能力。

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

$r_{\mathrm{\λ}}=\frac{L_{\mathrm{\λ}}-L_p}{{\mathrm{\μ}}_s{F}_0\mathrm{\τ}\left({\mathrm{\μ}}_s\right)\mathrm{\τ}\left({\mathrm{\μ}}_v\right)+(L_{\mathrm{\λ}}-L_p)S}$

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

步骤4：区域分割

从步骤1中的多光谱光学遥感图像数据中选取任意一块无云的晴空子区域，记该区域内的蓝光波段反射率为ρ_b，该区域内的红光波段反射率ρ_r；ρ_b、ρ_r均为矩阵形式，对ρ_b、ρ_r利用方程ρ_r＝t₁*ρ_b+t₂进行线性拟合，计算出t₁，t₂；由于ρ_b、ρ_r均为矩阵形式，在ρ_r＝t₁*ρ_b+t₂中，等号左右两边的形式也为矩阵，每一行或列对应相等，即实际存在多个等式，可以直接拟合出t₁，t₂，实际拟合时，等号左右两边的矩阵均可变形为列向量，再进行拟合。

通过公式计算得到HOT，其中为中蓝光波段在该晴空子区域上的子集，为中红光波段在该晴空子区域上的子集；

对HOT通过阈值分割方法得到二值图Mask₁；由计算得到NDWI，其中Nir为步骤1中得到的近红外波段反射率数据，为中的绿光波段反射率；将NDWI通过阈值分割方法划分得到二值图Mask₂；

特别的，当多光谱光学遥感图像数据不包含水体时，Mask₂为0矩阵。

二值图Mask₁由“0”和“1”两类元素组成，其中0代表晴空区域像素点的值，1代表薄云区域像素点的值；二值图Mask₂，二值图由“0”和“1”两类元素组成，其中0代表陆地区域像素点的值，1代表有水区域像素点的值。

阈值分割是一种常用图像分割方法，因其实现简单、计算量小、性能较稳定而成为图像分割中最基本和应用最广泛的分割技术。图像阈值分割的目的是要按照灰度级，对像素集合进行一个划分，得到的每个子集形成一个与现实景物相对应的区域，各个区域内部具有一致的属性，而相邻区域不具有这种一致属性。这样的划分可以通过从灰度级出发选取一个或多个阈值来实现。常用的特征包括：直接来自原始图像的灰度或彩色特征；由原始灰度或彩色值变换得到的特征。

例如设原始图像为f(x，y)，按照一定的准则f(x，y)中找到特征值T，将图像分割为两个部分A和B，若取A＝0，B＝1，即为图像二值化。根据HOT和NDWI的物理含义，优选的在步骤4中利用阈值分割方法得到二值图Mask₁、Mask₂过程中的阈值均为零。

步骤5：统计特征参数提取

将步骤3中的大气校正结果代入公式及公式分别得到晴空区域的反射率Clear_i和薄云覆盖的有水区域的反射率Cloud_i，其中符号“！”表示取反运算，i＝1，2，3，...，I；

在该步骤中计算的晴空区域的反射率中下标不同于Clear_i的记做Clear_j，即j＝1，2，3，...，I且j≠i；对Clear_i和Clear_j进行线性拟合，得到拟合线性方程为Clear_j＝k_ij*Clear_i+b_ij；对Cloud_i和Cloud_j进行线性拟合，对应的线性方程为Cloud_j＝m_ij*Cloud_i+n_ij；其中i分别取1，2，3，...，I，共得到I组(k_ij,b_ij,m_ij,n_ij)；

特别的，当多光谱光学遥感图像数据不包含水体时，m_ij＝1,n_ij＝0。

步骤6：利用去云辐射传输方程求解去云结果

将步骤5中提取的I组统计特征参数(k_ij,b_ij,m_ij,n_ij)及步骤3中的大气校正结果地表反射率矩阵关联矩阵代入去云辐射传输方程

其中i分别取1，2，3，...，I，j＝1，2，3，...，I且j≠i；T_ci为薄云对i波段反射率的双向透过率；即可分别求解得到对应i波段的去云结果后的表观反射率ρ_i，其中i＝1，2，3，...，I。

本发明所述可见光是电磁波谱中人眼可以感知的部分，一般人的眼睛可以感知的电磁波的波长在0.400～0.760um之间。可见光谱没有精确的范围，本发明中可见光波段特指波长范围在0.400～0.760um之间的电磁波。

对于步骤6中的去云辐射传输方程，本发明借鉴了GAO的研究，但本发明还基于以下两点假设：

1.可见光波段的晴空区域的地物表观反射率ρ_i与关联表观反射率ρ_j之间存在着不受地物类型影响的线性关系。

此处i分别取1，2，3，...，Ij＝1，2，3，...，I且j≠i；用公式描述为：

ρ_j＝k*ρ_i+b

上式中k，b分别为该线性方程系数

2.可见光波段的薄云反射率ρ_ci，关联薄云反射率ρ_cj线性相关，此处i分别取1，2，3，...，I；j＝1，2，3，...，I且j≠i；

用公式描述为：

${\mathrm{\ρ}}_{\mathrm{cj}}^{*}=m_{\mathrm{ij}}*{\mathrm{\ρ}}_{\mathrm{ci}}+n_{\mathrm{ij}}$

上式中m _ij ，n _ij 分别为该线性方程系数；

去云辐射传输方程的推导过程如下：

卫星接收到的i波段的反射率可以由薄云在该波段的反射率和地物的表观反射率线性表出，如公式(1)所示：

${\mathrm{\ρ}}_i^{*}={\mathrm{\ρ}}_{\mathrm{ci}}+T_{\mathrm{ci}}*{\mathrm{\ρ}}_i---\left(1\right)$

由假设1、2可知：

ρ_i＝k_ij*ρ_j+b_ij   (2)

ρ_ci＝m_ij*ρ_cj+n_ij   (3)

联立上述公式(1)、(2)、(3)可得i波段的表观反射率为：

${\mathrm{\ρ}}_{\mathrm{cj}}^{*}=m_{\mathrm{ij}}*{\mathrm{\ρ}}_{\mathrm{ci}}+n_{\mathrm{ij}}+k_{\mathrm{ij}}*T_{\mathrm{ci}}*{\mathrm{\ρ}}_j+T_{\mathrm{ci}}*b_{\mathrm{ij}}---\left(4\right)$

为消除公式(4)中的变量ρ_cj，引入j波段的辐射传输方程：

${\mathrm{\ρ}}_j^{*}={\mathrm{\ρ}}_{\mathrm{cj}}+T_{\mathrm{cj}}*{\mathrm{\ρ}}_j---\left(5\right)$

将公式(5)代入(4)中，则得到j波段的去云辐射传输方程：

${\mathrm{\ρ}}_j=\frac{{\mathrm{\ρ}}_i^{*}-m_{\mathrm{ij}}*{\mathrm{\ρ}}_j^{*}-b_{\mathrm{ij}}-n_{\mathrm{ij}}}{k_{\mathrm{ij}}*T_{\mathrm{ci}}-m_{\mathrm{ij}}*T_{\mathrm{cj}}}$

采用以上两个假设推导得到的去云辐射传输方程，对数据的要求简单，适用面广。适用于现有的绝大部分的光学卫星(搭载了可见光波段传感器的卫星，如Landsat、MODIS、SPOT等)图像数据，不需要其他来源图像或实测数据的辅助，克服了其他的单幅遥感图像去云方法数据要求苛刻、适用范围小等缺点。

该方法计算代价小，且可信度较高。与现有的薄云去除方法相比，如传统的同态滤波、IHS变换、小波变换等方法，不需要进行空频转换、小波变换及频域滤波等一系列时间复杂度较高的运算，在矫正云区图像数据的同时又很好的恢复了地物信息，可信度较高。

本发明中所公开的实施例描述的方法或算法的步骤可以直接用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储器(RAM)、内存、只读存储器(ROM)、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM、或技术领域内所公知的任意其它形式的存储介质中。

前文所述的为本发明的各个优选实施例，各个优选实施例中的优选实施方式如果不是明显自相矛盾或以某一优选实施方式为前提，各个优选实施方式都可以任意叠加组合使用，所述实施例以及实施例中的具体参数仅是为了清楚表述发明人的发明验证过程，并非用以限制本发明的专利保护范围，本发明的专利保护范围仍然以其权利要求书为准，凡是运用本发明的说明书及附图内容所作的等同结构变化，同理均应包含在本发明的保护范围内。

Claims (3)

1.基于可见光光谱统计特征的光学遥感图像薄云去除方法，其特征在于，包括如下步骤：

步骤1：数据准备

采集多光谱光学遥感图像数据及其头文件Hdr，它有I个可见光波段图层；从头文件Hdr中获得近红外波段反射率数据Nir；各波段修正因子和各波段偏移量太阳高度角θ_SE；其中i、I均为正整数且i＝1，2，3，...，I，i表示不同波段，I表示采集的波段总数；

步骤2：表观反射率转换

将步骤1中的多光谱光学数据代入公式中，分别求得I个波段的表观反射率其中i＝1，2，3，...，I；

步骤3：大气校正

将步骤2中各个波段的表观反射率矩阵进行大气校正，分别获得I个波段的地表反射率矩阵对每一下标i，在该步骤中计算的地表反射率矩阵中下标不同于的记做即j＝1，2，3，...，I且j≠i；

步骤4：区域分割

从步骤1中的多光谱光学遥感图像数据中选取任意一块无云的晴空子区域，记该区域内的蓝光波段反射率为ρ_b，该区域内的红光波段反射率为ρ_r；ρ_b、ρ_r均为矩阵形式，对ρ_b、ρ_r利用方程ρ_r＝t₁*ρ_b+t₂进行线性拟合，计算出t₁、t₂；通过公式计算得到HOT，其中为中蓝光波段在该晴空子区域上的子集，为中红光波段在该晴空子区域上的子集；

对HOT通过阈值分割方法得到二值图Mask₁；由计算得到NDWI，其中Nir为步骤1中得到的近红外波段反射率数据，为中的绿光波段反射率；将NDWI通过阈值分割方法划分得到二值图Mask₂；

步骤5：统计特征参数提取

将步骤3中的大气校正结果代入公式及公式分别得到晴空区域的反射率Clear_i和薄云覆盖的有水区域的反射率Cloud_i，其中符号“！”表示取反运算，i＝1，2，3，...，I；

在该步骤中计算的晴空区域的反射率中下标不同于Clear_i的记做Clear_j，即j＝1，2，3，...，I且j≠i；对Clear_i和Clear_j进行线性拟合，得到拟合线性方程为Clear_j＝k_ij*Clear_i+b_ij；对Cloud_i和Cloud_j进行线性拟合，对应的线性方程为Cloud_j＝m_ij*Cloud_i+n_ij；其中i分别取1，2，3，...，I，共得到I组(k_ij,b_ij,m_ij,n_ij)；

步骤6：利用去云辐射传输方程求解去云结果

将步骤5中提取的I组统计特征参数(k_ij,b_ij,m_ij,n_ij)及步骤3中的大气校正结果代入去云辐射传输方程

其中i分别取1，2，3，...，I；j＝1，2，3，...，I且j≠i；T_ci为薄云对i波段反射率的双向透过率；即可分别求解得到对应i波段的去云结果ρ_i，其中i＝1，2，3，...，I。

2.如权利要求1所述基于可见光光谱统计特征的光学遥感图像薄云去除方法，其特征在于，步骤4中所述蓝光波段反射率、红光波段反射率、绿光波段反射率分别为波长在0.450–0.515um、0.630–0.680um、0.525–0.600um范围内。

3.如权利要求1所述基于可见光光谱统计特征的光学遥感图像薄云去除方法，其特征在于，所述步骤4中利用阈值分割方法得到二值图Mask₁、Mask₂过程中的阈值均为零。