CN104820972B - 一种基于在轨分类统计的红外影像me噪声去除方法 - Google Patents

Abstract

一种基于在轨分类统计的红外影像ME噪声去除方法，利用卫星影像自身的辐射信息，使用ME检测模板遍历整幅影像检测ME噪声，对检测到的ME噪声使用迭代的方法进行去除，弥补了系统脉冲响应函数不足的问题。同时考虑到红外影像在不同辐亮度下ME噪声表现形式不同的特点，该方法引入分类的思想，对不同的辐亮度下ME噪声使用不同的方法进行去除。最后，使用自适应矩匹配对处理后的图像进行残余条纹噪声去除，进一步提高红外影像的辐射质量。

Description

一种基于在轨分类统计的红外影像ME噪声去除方法

技术领域

本发明属于航天遥感技术领域，涉及一种遥感红外影像的噪声处理方法。

背景技术

记忆效应(ME，Memory Effect)是光机扫描成像仪所成红外遥感影像的一种条带噪声，例如LANDSAT卫星上的TM的热红外影像、CBERS01卫星上的IRMSS影像。ME噪声主要表现为一景图像中某些扫描条带与其上下扫描条带之间存在整体的灰度差，灰度差可能会随着扫描角度不同而不同，在沿着扫描方向的扫描条带中，景物有明显亮暗突变的位置尤为明显，例如云或海岸线的边界位置，这会严重影响影像的目视效果。

目前可查到的ME噪声去除算法为LANDSAT卫星上的TM热红外影像ME校正方法，该方法的基本原理为：使用系统的脉冲响应函数构建复原滤波器，使用复原滤波器对图像进行复原从而去除ME噪声。使用该方法的前提为确定精确的遥感器系统脉冲响应函数，但是精确的系统脉冲响应函数比较难测得，这对于很多普通的遥感影像用户来说非常困难，所以，此类方法的ME噪声去除很难实现。

发明内容

本发明解决的技术问题是：克服现有ME噪声去除方法过于依赖精确的遥感器系统脉冲响应函数的问题，提供了一种基于在轨分类统计的红外影像ME噪声去除方法，根据卫星影像自身的辐射信息，使用ME检测模板遍历整幅影像检测ME噪声，对检测到的ME噪声使用迭代的方法进行分类去除，可以在去除红外影像ME噪声的同时保证影像辐射信息的有效性。

本发明的技术解决方案是：一种基于在轨分类统计的红外影像ME噪声去除方法，包括如下步骤：

(1)针对卫星在轨获取的红外原始地物影像，依据影像各像素点的DN值使用K均值算法进行分类，得到K个辐亮度分类区间C₁,C₂,...,C_K；

(2)建立ME检测模板，并利用ME检测模板对原始影像进行遍历，进行匹配检测ME噪声；

所述的ME检测模板其中E为N×M的对角线元素为1其余元素为0的矩阵，N为获取原始地物影像的探测器的像元数，M为1到10的之间的任一正整数；

所述匹配检测ME噪声的方法为：将原始地物影像分成若干个大小为3N×M的子图像块，计算ME检测模板G与各子图像块的相关系数r，当r大于阈值的时候认为该图像块存在ME噪声，

其中，为ME检测模板G的DN均值，是图像块所包括像素点的DN均值，Y(i,j)为子图像块第i行第j列的像素点的DN值，G(i,j)为模板第i行第j列的像素点的DN值，i＝1,2,3,…,3N，j＝1,2,3,…,M；

(3)将步骤(2)得到的存在ME噪声的图像块Y进行分块，使得Y₀，Y₁，Y₂为N×M的矩阵，然后根据存在ME噪声的图像块Y所在的辐亮度分类区间分别进行ME噪声去除，去除方法为：

如果图像块那么去除ME噪声后对应的图像块其中±与C_i对应的ME_i的正负相反；

然后判断图像块Y′是否同时满足如下两个最优条件：1)2)|y′₁-y′₀|>0，|y′₁-y′₂|>0；如果满足则进入下一步；否则令Y＝Y′，并再次进行ME噪声去除，直至最后一次得到的图像块Y′同时满足两个最优条件后进入下一步；

式中b₁为子影像Y₀最后一行与子影像Y₁第一行对应各像素点的DN值差，b₂为子影像Y₁最后一行与子影像Y₂第一行对应各像素点的DN值差，b′₁为子影像Y₀′最后一行与子影像Y₁′第一行对应各像素点的DN值差，b₂′为子影像Y₁′最后一行与子影像Y₂第一行对应各像素点的DN值差，为b₁的DN均值，为b₂的DN均值，为b′₁的DN均值，为b′₂的DN均值，y₀，y₁，y₂分别为Y₀，Y₁，Y₂的DN均值，y′₀，y′₁，y′₂分别为Y₀′，Y₁′，Y₂′的DN均值；

(4)对去除ME噪声后的影像，采用自适应矩匹配滤波的方法去除影像中残余的条纹噪声，得到最终的红外影像。

所述的K均值算法中辐亮度分类区间的个数K由原始地物影像的ME噪声表现形式决定，满足关系式ME_i＝DN_i-DN′_i，式中，为ME噪声去除前各像素点的DN值，i＝1,2,…,K，DN′_i为ME噪声去除后各像素点的DN值，ME_i为正或负并且相邻的区间段正负号不同。

所述的K均值算法中，K个辐亮度分类区间C₁,C₂,...,C_K的初始聚类中心分别为其中N_bit为影像的量化位数。

本发明与现有技术相比的优点在于：

(1)本发明方法利用卫星影像自身的辐射信息，使用ME检测模板遍历整幅影像检测ME噪声，对检测到的ME噪声使用迭代的方法进行去除，可以有效解决现有方法过于依赖系统脉冲响应函数的问题；

(2)本发明方法充分考虑了红外影像在不同辐亮度下ME噪声表现形式不同的特点，利用红外影像对地物进行分类处理，对不同的辐亮度下ME噪声使用不同的方法进行去除，可以在去除红外影像ME噪声的同时保证影像辐射信息的有效性；

(3)本发明方法根据ME噪声的主要特征生成检测模板，对整个图像遍历检测，确保了ME噪声检测的准确率；

(4)本发明方法最终使用自适应矩匹配去除残余的条纹噪声，进一步提高了红外影像的辐射质量。

附图说明

图1为本发明方法的流程框图；

图2为本发明K均值算法处理示意图。

具体实施方式

如图1所示，为本发明方法的流程框图，主要步骤如下：

步骤1：使用K均值算法，对卫星在轨影像获取的地物依据灰度差异进行分类，从而获取辐亮度分类区间。

K均值聚类算法是一类划分聚类算法，属于硬聚类，最终使目标函数达到最小值。根据初始聚类中心对数据集中的数据对象进行分类并重新计算聚类中心和数据对象分类，迭代过程的结束标志着聚类准则函数已达收敛。

假设带有ME噪声的原始影像为X，H和L为原始影像X的长和宽。以影像的灰度量化信息作为样本集X(x₁,x₂,...,x_n)，现要将n(n＝H×L)个样本x_i划分为K个类别C₁,C₂,...,C_K，K的取值是由遥感图像的ME噪声表现形式来决定的，ME噪声在这K个不同亮度区间段表现形式不同，如式1所示。

ME_i＝DN_i-DN′_i (1)

式中，为ME噪声去除前的DN值，i＝1,2,…,K，DN′_i为ME噪声去除后的DN值，ME_i为正或负，并且相邻的区间段正负号不同。

取距离函数为欧几里得距离函数d(V_i,V_j)＝||V_i-V_j||，其中V_i、V_j为向量表示，K均值算法具体实现描述如下：

i.样本聚类中心初始化。

选择K个固定的样本作为初始聚类中心m₁,m₂,...,m_K，此处为了方便计算，假设影像为N_bit位量化，可以选择影像量化DN值区间的K个中心点作为对应的K个初始聚类中心：

ii.K均值算法类聚。

如果d(x_j,m_p)≤d(x_j,m_i)，1≤p≤K，i＝1,2,...,K，则分配x_j到第p类，j＝1,2,3...,n，具体如图2所示，m₁,m₂,...,m_K为K个样本的聚类中心值，针对第j个像元x_j，j＝1,2,3...,n，计算其到每个聚类中心的距离：d(x_j,m_i)，i＝1,2,...,K，如果x_j到第p类的距离最小，即d(x_j,m_p)≤d(x_j,m_i)，1≤p≤K，i＝1,2,...,K，那么x_j分配到第p个聚类中，即x_j∈C_p。

iii.重新计算每个聚类的中心。

式中，N_i是第i聚类C_i中的样本数目。

重复步骤ii和iii直到m_i不再变化。

步骤2：使用检测模板对整幅影像进行遍历，检测ME噪声。

i.ME检测模板生成。

由于ME噪声主要表现为某些扫描条带与其上下扫描条带之间存在整体的灰度差，所以，根据成像探测器像元数生成ME检测模板模板大小为3N×M，其中N为探测器像元数，M可调(为1到10的正整数)，E为对角线元素为1其余为0的N×M矩阵。

ii.使用模板遍历整幅影像进行匹配检测ME噪声。

模板匹配的方法为：使用模板G与待处理影像中相同大小的子影像Y计算其相关系数，当r大于阈值r₀(r₀一般取0.7，也可根据实际情况略做调整)的时候认为此处存在ME噪声，计算公式如下：

其中，为是检测模板G的DN均值(各元素值相加再除以元素总个数)，是子影像Y的DN均值，Y(i,j)为子影像Y第i行第j列像元的DN值，G(i,j)为模板第i行第j列位置元素的DN值(0或者1)，i＝1,2,3,…,3N，j＝1,2,3,…,M。

步骤3：对上步中检测出的ME噪声，根据其所属地物辐亮度区间进行相应的去除。

i.计算ME噪声所在位置的子影像Y的均值。

假设上步检测到的ME噪声所在位置的子影像为Y₀，Y₁，Y₂为N×M的矩阵，然后计算Y，Y₀，Y₁，Y₂的DN均值y，y₀，y₁，y₂；

ii.判断子影像Y属于哪个聚类区间，然后分别进行相应的ME噪声去除处理。

如果y∈C_i，那么根据此样本集C_i的ME特性进行如下处理：其中，±与此样本集内ME_i的正负相反。

iii.判断处理后的影像Y′是否同时满足最优条件，如果满足条件则停止运算，否则，令Y＝Y′，然后继续步骤ii的处理。

考虑到相邻扫描条带的灰度连续性，判断处理后的影像Y′是否同时满足最优条件1)：条件2)：|y′₁-y′₀|>0，|y′₁-y′₂|>0；其中b₁为子影像Y₀最后一行与子影像Y₁第一行各对应像元的DN值差，b₂为子影像Y₁最后一行与子影像Y₂第一行各对应像元的DN值差，b′₁为子影像Y₀′最后一行与子影像Y₁′第一行各对应像元的DN值差，b₂′为子影像Y₁′最后一行与子影像Y₂第一行各对应像元的DN值差，为b₁的DN均值，为b₂的DN均值，为b′₁的DN均值，为b′₂的DN均值，y₀，y₁，y₂分别为Y₀，Y₁，Y₂的DN均值，y′₀，y′₁，y′₂分别为Y₀′，Y₁′，Y₂′的DN均值。

步骤4：对上步得到的结果进行自适应矩匹配滤波去除残余的条纹噪声。

自适应矩匹配滤波去除残余的条纹噪声时，主要包括获取亮度值曲线图、去基线处理、获取处理单元中心、去除伪中心点、获取处理单元、对每个处理单元进行矩匹配处理等过程。具体的方法可参见文章：张炳先，王密，潘俊，“采用灰度分割的自适应矩匹配条带噪声去除方法”，武汉大学学报.信息科学版，2012，37(12):1465。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (2)

1.一种基于在轨分类统计的红外影像ME噪声去除方法，其特征在于包括如下步骤：

(1)针对卫星在轨获取的红外原始地物影像，依据影像各像素点的DN值使用K均值算法进行分类，得到K个辐亮度分类区间C₁,C₂,...,C_K；所述的K均值算法中辐亮度分类区间的个数K由原始地物影像的ME噪声表现形式决定，满足关系式ME_i＝DN_i-DN′_i，式中，为ME噪声去除前各像素点的DN值，i＝1,2,…,K，DN′_i为ME噪声去除后各像素点的DN值，ME_i为正或负并且相邻的区间段正负号不同；

(2)建立ME检测模板，并利用ME检测模板对原始影像进行遍历，进行匹配检测ME噪声；

所述的ME检测模板其中E为N×M的对角线元素为1其余元素为0的矩阵，N为获取原始地物影像的探测器的像元数，M为1到10的之间的任一正整数；

所述匹配检测ME噪声的方法为：将原始地物影像分成若干个大小为3N×M的子图像块，计算ME检测模板G与各子图像块的相关系数r，当r大于阈值的时候认为该图像块存在ME噪声，

$r=\frac{\underset{i=1}{\overset{3N}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\[Y(i,j)-\stackrel{\‾}{Y}\]\×\[G(i,j)-\stackrel{\‾}{G}\]}{\sqrt{\underset{i=1}{\overset{3N}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\[Y(i,j)-\stackrel{\‾}{Y}\]}^2\underset{i=1}{\overset{3N}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\[G(i,j)-\stackrel{\‾}{G}\]}^2}}$

其中，为ME检测模板G的DN均值，是图像块所包括像素点的DN均值，Y(i,j)为子图像块第i行第j列的像素点的DN值，G(i,j)为模板第i行第j列的像素点的DN值，i＝1,2,3,…,3N，j＝1,2,3,…,M；

(3)将步骤(2)得到的存在ME噪声的图像块Y进行分块，使得Y₀，Y₁，Y₂为N×M的矩阵，然后根据存在ME噪声的图像块Y所在的辐亮度分类区间分别进行ME噪声去除，去除方法为：

如果图像块那么去除ME噪声后对应的图像块其中±与C_i对应的ME_i的正负相反；

然后判断图像块Y′是否同时满足如下两个最优条件：1)2)|y′₁-y′₀|>0，|y′₁-y′₂|>0；如果满足则进入下一步；否则令Y＝Y′，并再次进行ME噪声去除，直至最后一次得到的图像块Y′同时满足两个最优条件后进入下一步；

式中b₁为子影像Y₀最后一行与子影像Y₁第一行对应各像素点的DN值差，b₂为子影像Y₁最后一行与子影像Y₂第一行对应各像素点的DN值差，b′₁为子影像Y₀′最后一行与子影像Y₁′第一行对应各像素点的DN值差，b₂′为子影像Y₁′最后一行与子影像Y₂第一行对应各像素点的DN值差，为b₁的DN均值，为b₂的DN均值，为b′₁的DN均值，为b′₂的DN均值，y₀，y₁，y₂分别为Y₀，Y₁，Y₂的DN均值，y′₀，y′₁，y′₂分别为Y₀′，Y₁′，Y₂′的DN均值；

(4)对去除ME噪声后的影像，采用自适应矩匹配滤波的方法去除影像中残余的条纹噪声，得到最终的红外影像。

2.根据权利要求1所述的一种基于在轨分类统计的红外影像ME噪声去除方法，其特征在于：所述的K均值算法中，K个辐亮度分类区间C₁,C₂,...,C_K的初始聚类中心分别为其中N_bit为影像的量化位数。