CN102279971B - 基于变分方法的中分辨率成像光谱仪条带噪声去除方法 - Google Patents

Abstract

本发明公开了一种基于变分方法的中分辨率成像光谱仪(MODIS)条带噪声去除方法，首先对MODIS数据按照分辨率进行数据分离；其次将获得的探测器子图像分为两类，一类是不含条带噪声或条带噪声不明显的子图像

一类是含有明显条带噪声的子图像

接着对含有明显条带噪声的子图像

利用变分方法进行条带噪声去除；最后将所有去除条带后探测器子图像

和不含条带噪声子图像

按照编号顺序进行合成，得到去条带后图像。本发明能够实现中分辨率成像光谱仪数据的探测器间条带和随机噪声条带的并行去除，避免了去除一种条带噪声而使另外一种条带噪声得到增强的情况，并且可在有效解决条带噪声影响的同时最大限度的保留原数据的物理特性。

Description

基于变分方法的中分辨率成像光谱仪条带噪声去除方法

技术领域

本发明属于遥感图像处理领域，特别是一种基于变分方法的中分辨率成像光谱仪条带噪声去除方法。

背景技术

中分辨率成像光谱仪(MODIS)是一种先进的空间对地观测传感器。该传感器具有36个波段，采用多元并扫方式进行扫描，其中250m分辨率40元并扫，500m分辨率20元并扫，1km分辨率10元并扫。多元并扫技术在提高数据采集和扫描效率的同时，由于多个探测器之间很难达到完全一致的响应水平，加上空间复杂的电磁环境和器件本身的误差导致了条带噪声的产生。而条带噪声的存在不仅严重的影响了波谱特征增强及空间特征增强和提取的效果，对于进行大气校正进而定量反演水体等光学特性来说也是十分不利的(因为条带效应引入了处理中的噪声，如若不进行条带的消除工作，进行大气校正的结果就会退化)。因此，必须进行条带去除，这对于遥感影像数据后续应用的有效性至关重要。

目前，用于去除遥感影像条带噪声的方法主要分为两大类：一类是通过傅立叶变换、小波变换等频率域滤波的方法通过滤波算子去除周期性噪声的频率成分，然后反变换回空域获得去噪后图像。该方法存在的一个很大问题就是不易选择正确的频率成分，而将处于条带频率的非条带的有用信息去除掉。另一类是针对图像本身的灰度值特征进行归一化和匹配的方法，典型的有直方图匹配和矩匹配方法。在实际应用中，通常每个通道选择一个确定的探测器作为参考，其他探测器以此为标准进行调整，从而实现图像条带的消除。

MODIS数据中存在三种条带噪声：镜边条带、探测器间条带和随机噪声条带。为了去除这三种条带噪声，通常的做法是针对三种条带噪声特性分别进行去除。PreesanRakwatin在“stripe noise reduction in MODIS data by combining histogram matching withfacet filter”一文中提出先利用直方图匹配方法去除探测器间条带，再利用迭代加权最小二乘小面滤波法去除随机噪声条带。但是这种方法在去除探测器间条带的同时可能增强随机噪声，使后续的去噪处理难以奏效，同时，遥感图像的重要细节部分会受到模糊。

发明内容

本发明所解决的技术问题在于提供一种基于变分方法的中分辨率成像光谱仪条带噪声去除方法，从而能够实现探测器间条带和随机噪声条带的并行有效去除，并最大限度地保持数据自身的物理特性。

实现本发明目的的技术解决方案为：一种基于变分方法的中分辨率成像光谱仪条带噪声去除方法，包括以下步骤：

步骤1、对MODIS数据按照分辨率进行数据分离；

步骤2、将步骤1获得的探测器子图像分为两类，一类是不含条带噪声或条带噪声不明显的子图像

一类是含有明显条带噪声的子图像

下标n₁和n₂为探测器子图像编号，且n₁和n₂编号不重复；

步骤3、对含有明显条带噪声的子图像利用变分方法进行条带噪声去除；

步骤4、将所有去除条带后探测器子图像

和不含条带噪声子图像

按照编号顺序进行合成，得到去条带后图像。

本发明与现有技术相比，其显著优点为：本发明针对MODIS数据含有三种条带噪声的特点，能够实现探测器间条带和随机噪声条带的并行去除，避免了去除一种条带噪声而使另外一种条带噪声得到增强的情况，并且可在有效解决条带噪声影响的同时最大限度的保留原数据的物理特性。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是本发明的基于变分方法的中分辨率成像光谱仪条带噪声去除方法流程图。

图2是本发明的实施例所用数据。

图3是分离探测器子图像，其中(a)-(j)分别对应1-10号探测器子图像。

图4是整幅图像去条带结果图，其中(a)为具体实例原图(b)为基于传统直方图匹配方法去条带结果图像(c)为基于变分方法去条带结果图像。

具体实施方式

本发明的一种基于变分方法的中分辨率成像光谱仪条带噪声去除方法，包括以下步骤：

步骤1、对MODIS数据按照分辨率进行数据分离；对MODIS数据按照分辨率进行数据分离时，由于MODIS的250m分辨率数据为40元并扫，500m分辨率数据为20元并扫，1km分辨率数据为10元并扫，同时考虑到镜像扫描，将250m分辨率数据看作2组40个探测器共80个探测器，500m分辨率数据看作2组20个探测器共40个探测器，1km分辨率数据看作2组10个探测器共20个探测器；设定数据中的第一行由一号探测器生成，第二行由二号探测器生成，以此类推，得到对应数量的探测器子图像，编号从1开始。

步骤2、将步骤1获得的探测器子图像分为两类，一类是不含条带噪声或条带噪声不明显的子图像

一类是含有明显条带噪声的子图像

下标n₁和n₂为探测器子图像编号，且n₁和n₂编号不重复；

步骤3、对含有明显条带噪声的子图像

利用变分方法进行条带噪声去除；利用变分方法去除含有明显条带噪声的子图像

中的条带噪声具体包括以下步骤：

步骤31、从含条带噪声子图像集合中选取一幅子图像，设为

步骤32、对于子图像

利用相关性计算公式，从

中计算得到与其相关性最强的相邻子图像

作为参考探测器，相关性计算公式如下：

$H_{n_1}=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}[X_{n_1}(i,j)-{\mathrm{\μ}}_{{\mathrm{Xn}}_1}][Y_{K_1}(i,j)-{\mathrm{\μ}}_{Y_{K_1}}]}{\sqrt{\left\{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{[X_{n_1}(i,j)-{\mathrm{\μ}}_{{\mathrm{Xn}}_1}]}^2\right\}\left\{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{[Y_{K_1}(i,j)-{\mathrm{\μ}}_{Y_{K_1}}]}^2\right\}}}$

式中：H为两相邻子图像的相关系数，

表示探测器子图像中坐标(i，j)处的值，

表示探测器子图像

中坐标(i，j)处的值，M和N分别为子图像的高和宽，μ表示子图像的灰度均值；

步骤33、利用基于变分的直方图匹配法去除探测器间条带噪声：

$\frac{\∂Y_{K_1}}{\∂t}=[M\×N-H\left(Y_{K_1}\right)]-A[(v,w):Y_{K_1}(v,w,t)\≥X_{K_2}(x,y,t)]$

式中，

是输入的第K₁号含条带噪声子图像，

是

最相关的相邻探测器子图像，H(·)是

的累计直方图，A(·)表示面积函数，即

中满足

的像素个数；

步骤34、利用正则化各向异性扩散方程消除随机条带噪声：

$\frac{\∂Y_{K_1}}{\∂t}=\mathrm{div}[g\left(\right|\▿Y_{K_1\mathrm{\σ}}\left|\right)\▿Y_{K_1}]$

其中，div为散度算子，

表示梯度，g(·)是选择的边缘函数，

G_σ表示方差为σ的高斯函数，^*表示卷积；

步骤35、利用梯度下降法得到步骤33和步骤34的离散数值求解公式，公式为：

$Y_{K_1}^{n+1}=Y_{K_1}^n+\mathrm{\Δt}\left\{\right[M\×N-H\left(Y_{K_1}^n\right)]-A[(v,w)|Y_{K_1}^n(v,w,t)\≥X_{K_2}^n(x,y,t)]$

$+\mathrm{\λ}\mathrm{div}[g\left(\right|\▿Y_{K_1\mathrm{\σ}}^n\left|\right)\▿Y_{K_1}^n]\}$

式中，上标n和n+1分别表示在n和n+1时刻的图像值，Δt为时间迭代参数，参数λ用来平衡直方图匹配法和各向异性扩散去噪方法对图像的贡献；

步骤36、利用步骤35中的公式对进行求解，在求解过程中的每次迭代后，判断是否满足迭代终止条件，如果满足则终止迭代，并执行步骤37，否则继续执行本步骤直至满足迭代终止条件；所述迭代终止条件为：1)前后两次迭代图像之差的平方和除以当前图像平方和是否小于预设阈值ε，或2)是否达到迭代次数，所述迭代次数根据需要进行预设；

步骤37、判断含条带噪声子图像集合

是否处理完毕，如果处理完毕则结束，否则选取下一幅子图像，从步骤32开始继续处理，对于每个

得到对应去除条带后探测器子图像

步骤4、将所有去除条带后探测器子图像

和不含条带噪声子图像按照编号顺序进行合成，得到去条带后图像。

下面结合实施例对本发明做进一步详细的描述：

本实施例的流程如图1所示，所用数据为2005年4月13日Terra MODIS的1KM分辨率的波段27数据，采用10元并扫，图像大小为2000×200，如图2所示。具体步骤为：

第一步，对MODIS数据按照其空间分辨率进行数据分离。为了确定对应的探测器子图像，设原始影像的第一行数据由MODIS传感器的第一个探测器扫描产生，第二行数据由第二个探测器扫描产生，依此类推。

由于MODIS使用了双边扫描镜，为了去除镜边条带噪声，处理的时候需要将MODIS数据看作由20个探测器组成，即每10个探测器对应一个扫描镜边。从而可将原始影像数据划分为20幅子图像(如图3(a)-(j)所示，分别对应1-10号探测器子图像)，每幅子图大小为100×200，编号从1到20。

第二步，将20个探测器子图像分为含条带的子图像和不含条带的子图像，如表1所示。

表1探测器划分结果

第三步，对含条带噪声的子图像利用变分方法进行条带噪声去除。

(1)从含条带噪声子图像集合中选取第一幅子图像，1号子图像。

(2)利用相关性计算公式找出与选取条带子图像最相关的相邻参考子图像。

(3)利用离散数值求解公式对选取条带子图像进行迭代计算。在本实例中选取迭代参数Δt＝0.02，λ＝0.2。

(4)判断是否满足迭代终止条件。在本实例中选取ε＝10^-3，迭代次数为20。

(5)判断含条带噪声子图像集合是否处理完毕，如果是则转第四步，否则选取下一幅条带子图像，从(2)开始处理，直到将所有含有条带的探测器子图像处理完毕。

第四步，将去除条带后的各子图像及不含条带的子图像按照原顺序编号进行组合，得到去条带后图像，如图4(c)所示。

为了对比去条带效果，同时给出原图(如图4(a)所示)和基于传统直方图匹配方法去条带结果图像(如图4(b)所示)。

由结果图像和原图像比较可知，提出的方法可以实现探测器间条带和随机噪声条带的并行去除，避免了去除一种条带噪声而使另外一种条带噪声得到增强的情况，并且可在有效解决条带噪声影响的同时最大限度的保留原数据的物理特性。

Claims (1)

1.一种基于变分方法的中分辨率成像光谱仪条带噪声去除方法，其特征在于，包括以下步骤：

步骤1、对MODIS数据按照分辨率进行数据分离，由于MODIS的250m分辨率数据为40元并扫，500m分辨率数据为20元并扫，1km分辨率数据为10元并扫，同时考虑到镜像扫描，将250m分辨率数据看作2组40个探测器共80个探测器，500m分辨率数据看作2组20个探测器共40个探测器，1km分辨率数据看作2组10个探测器共20个探测器；设定数据中的第一行由一号探测器生成，第二行由二号探测器生成，以此类推，得到对应数量的探测器子图像，编号从1开始；

步骤2、将步骤1获得的探测器子图像分为两类，一类是不含条带噪声或条带噪声不明显的子图像

一类是含有明显条带噪声的子图像下标n₁和n₂为探测器子图像编号，且n₁和n₂编号不重复；

步骤3、对含有明显条带噪声的子图像

利用变分方法进行条带噪声去除，具体包括以下步骤：

步骤31、从含条带噪声子图像集合

中选取一幅子图像，设为

步骤32、对于子图像利用相关性计算公式，从

中计算得到与其相关性最强的相邻子图像

作为参考探测器，相关性计算公式如下：

$H_{n_1}=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}[X_{n_1}(i,j)-{\mathrm{\μ}}_{{\mathrm{Xn}}_1}][Y_{K_1}(i,j)-{\mathrm{\μ}}_{Y_{K_1}}]}{\sqrt{\left\{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\right[X_{n_1}(i,j)-{\mathrm{\μ}}_{{\mathrm{Xn}}_1}{]}^2\left\}\right\{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{[Y_{K_1}(i,j)-{\mathrm{\μ}}_{K_1}]}^2\}}}$

式中：H_n1为两相邻子图像的相关系数，

表示探测器子图像

中坐标(i，j)处的值，

表示子图像

中坐标(i，j)处的值，M和N分别为子图像的高和宽，μ表示子图像的灰度均值；

步骤33、利用基于变分的直方图匹配法去除探测器间条带噪声：

$\frac{{\∂Y}_{K_1}}{\∂t}=[M\×N-H\left(Y_{K_1}\right)]-A[(v,w):Y_{K_1}(v,w,t)\≥X_{K_2}(x,y,t)]$

式中，

是输入的第K₁号含条带噪声子图像，

是

相关性最强的相邻子图像，H(·)是的累计直方图，A(·)表示面积函数，即中满足 $Y_{K_1}(v,w,t)\≥X_{K_2}(x,y,t)$ 的像素个数；

步骤34、利用正则化各向异性扩散方程消除随机条带噪声：

$\frac{\∂Y_{K_1}}{\∂t}=\mathrm{div}[g\left(\right|{\▿Y}_{K_1\mathrm{\σ}}\left|\right)\▿Y_{K_1}]$

其中，div为散度算子，

表示梯度，g(·)是选择的边缘函数，

G_σ表示方差为σ的高斯函数，^*表示卷积；

步骤35、利用梯度下降法得到步骤33和步骤34的离散数值求解公式，公式为：

$Y_{K_1}^{n+1}=Y_{K_1}^n+\mathrm{\Δt}\left\{\right[M\×N-H\left(Y_{K_1}^n\right)]-A[(v,w):Y_{K_1}^n(v,w,t)\≥X_{K_2}^n(x,y,t)]$

$+\mathrm{\λ}\mathrm{div}[g\left(\right|{\▿Y}_{K_1\mathrm{\σ}}^n\left|\right)\▿Y_{K_1}^n]$

式中，上标n和n+1分别表示在n和n+1时刻的图像值，Δt为时间迭代参数，参数λ用来平衡直方图匹配法和各向异性扩散去噪方法对图像的贡献；

步骤36、利用步骤35中的公式对进行求解，在求解过程中的每次迭代后，判断是否满足迭代终止条件，如果满足则终止迭代，并执行步骤37，否则继续执行本步骤直至满足迭代终止条件；所述迭代终止条件为：1)前后两次迭代图像之差的平方和除以当前图像平方和是否小于预设阈值ε，或2)是否达到迭代次数，所述迭代次数根据需要进行预设；

步骤37、判断含条带噪声子图像集合是否处理完毕，如果处理完毕则结束，否则选取下一幅子图像，从步骤32开始继续处理，对于每个

得到对应去除条带后探测器子图像

步骤4、将所有去除条带后探测器子图像

和不含条带噪声子图像

按照编号顺序进行合成，得到去条带后图像。

Images