CN102819827A - 一种基于灰度分割的自适应矩匹配条带噪声去除方法 - Google Patents

Abstract

一种基于灰度分割的自适应矩匹配条带噪声去除方法，包括获取原始含噪影像的亮度值曲线；通过采用不同大小的处理窗口对亮度值曲线进行中值滤波去除基线，获取处理单元中心；根据两个相邻处理单元中心的间距进行调整，确定处理单元；对单个处理单元采用标准矩匹配的处理方式进行影像处理。该方法针对数据量大的遥感全色影像，基于基本影像处理理论，提供了一种有效的影像噪声去除方法。该方法能够有效地复原影像、改善影像质量，算法效率也很高。

Description

一种基于灰度分割的自适应矩匹配条带噪声去除方法

技术领域

本发明属于遥感影像处理领域，涉及一种基于灰度分割的自适应矩匹配条带噪声去除方法。

背景技术

条带噪声是影响光学卫星影像成像质量的一个重要因素，抑制或去除条带噪声是卫星地面预处理进行辐射处理的基本环节之一。光学卫星的成像系统由于内外部因素的干扰,如：CCD器件（电荷耦合器件）的性能随时间发生变化、大气干扰等，在进行影像均一化相对辐射校正之后，影像中仍然会残留有条带噪声，这类噪声的存在，极大的降低了影像的清晰度，为影像的后续判读处理增加了难度，因此必须予以剔除。

通过对影像中残留的条带噪声进行分析发现该类条带噪声存在有以下特点：1.噪声出现的空间位置是随机的。2.噪声与周围地物之间高度非线性相关。

目前常用的去条带噪声的方法归纳起来可以分为两类：一类是针对影像空间域特征提出的去噪方法；另一类是将空间域和频率域相结合，采用适当的滤波算子去除条带噪声的方法。其中空间域去噪的典型算法有Wegener提出的基于改进直方图匹配的多载荷成像的条带噪声去除法，该方法为了得到期望直方图，需要统计大量的有代表性的卫星影像，工作量巨大,不利于单景影像实时处理，而且影像中条带噪声去除的效果依赖于样本影像的选取，因此该算法的处理效果不稳定。于是，Gadallahd等人提出了基于矩匹配的多载荷影像的条带噪声去除法，该方法要求影像灰度均一，且条带噪声影像与其相邻影像之间存在较高的相关性，这类算法对于处理局部的地物单一的卫星影像具有很好的效果，但是处理整幅地物丰富的卫星影像，效果不佳。在此基础上，刘正军等人提出了基于均值补偿的矩匹配方法，该方法利用高光谱影像的多光谱信息来纠正矩匹配算法产生的灰度畸变，该方法很好的处理了标准矩匹配算法纠正影像时造成的探元间的灰度畸变，但是并没有纠正标准矩匹配算法造成的探元内的灰度畸变，且该方法对于单波段的全色影像并不适用。空间域和频率域相结合的去噪算法目前主要是利用小波变换的时频特点，通过对影像进行小波变换，研究噪声的小波系数变换规律，从而提炼出条带噪声的成分并对其进行剔出，但小波变换计算量大，并且会对光谱信息造成较大损失，对于整幅卫星影像的处理，计算速度和光谱保留度并不理想，同时由于影像中条带噪声空间位置是随机的，因此完全定位该类条带噪声在影像中的位置，可实现性尚存在质疑。

发明内容

本发明所要解决的问题是，针对全色影像中存在的随机条带噪声，提供一种有效的影像噪声去除方法。

本发明的技术方案为一种基于灰度分割的自适应矩匹配条带噪声去除方法，包括以下步骤：

步骤1，获取原始含噪影像的亮度值曲线；

步骤2，通过采用不同大小的处理窗口对步骤1所得亮度值曲线进行中值滤波去除基线，获取处理单元中心；

步骤3，根据两个相邻处理单元中心的间距进行调整，确定处理单元；

步骤4，对单个处理单元采用标准矩匹配的处理方式进行影像处理。

而且，设原始含噪影像的大小为M×N，(i,j)表示原始含噪影像中第i行第j列的像素，其中i的取值范围为0,1…M-1，j的取值范围为0,1…N-1，步骤1中获取原始含噪影像的亮度值曲线是根据如下公式计算得到，

$S\left(k\right)=\left\{\underset{i=0,j=0}{\overset{i=M-1,j=N-1}{\mathrm{\Σ}}}{I}_E(i,j)\right|I_E(i,j)=k\}$

其中，k表示对应的灰度级，I_E(i,j)表示影像中位于原始含噪影像(i,j)处像素的灰度值。

而且，步骤2包括如下子步骤，

步骤2.1，设中值滤波的处理窗口大小记为W_F，初始化设置W_F的值，设迭代次数Count=0；

步骤2.2，执行以下计算公式获取去基线的亮度值曲线S(k)_f；

S(k)_f＝S(k)-med(S(k-t)|t∈W)

其中，W＝{t|-(W_F-1)/2≤t≤(W_F-1)/2}，med是指中值滤波处理；

步骤2.3，对去基线的亮度值曲线S(k)_f采用均值μ(k)_f和方差σ(k)_f进行归一化处理，处理公式如下，

$S{\left(k\right)}_n=\frac{S{\left(k\right)}_f-\mathrm{\μ}{\left(k\right)}_f}{\mathrm{\σ}{\left(k\right)}_f}$

其中，S(k)_n就是均一化后的结果；

步骤2.4，根据预设的阈值T选择候选处理中心，选择公式如下，

$P{\left(k\right)}_c=\left\{\begin{array}{cc}1& \left|S{\left(k\right)}_n\right|>T\\ 0& \end{array}\right.$

当去基线的亮度值曲线S(k)_f上某点的候选值P(k)_c为1时，该点是候选处理中心，否则不是；

步骤2.5，令Count=Count+1，判断Count是否达到预设次数，否则进入步骤2.6，是则减小W_F的值，返回步骤2.2迭代执行中值滤波；

步骤2.6，对迭代所得的所有候选处理中心进行过滤，

首先计算判断函数值，计算公式如下，

$P\left(k\right)=\left\{\begin{array}{cc}1& P{\left(k\right)}_c>\mathrm{Count}/2\\ 0& \end{array}\right.$

将某候选处理中心的判断函数值P(k)带入到步骤1所得亮度值曲线S(k)中，若对应的像素个数大于预设阈值I，则确定为处理单元中心，否则去除。

而且，预设阈值I采用以下公式计算，

$I=\frac{1}{2^m-1}$

其中，m表示影像量化的位数。

而且，步骤3中对处理单元中心进行调整，实现方式为：

当两个相邻处理单元中心的间距小于预设阈值h时，则将两个相邻处理单元中心合并为一个处理单元中心，合并后的处理单元中心的值为这两个相邻处理单元中心的均值；

当两个相邻处理单元中心的间距大于或等于预设阈值h时，两个相邻处理单元中心之间分裂产生新的处理单元中心，分裂后新的处理单元中心的值为原始两相邻处理单元中心的平均值。

而且，预设阈值h采用原始含噪影像的随机噪声灰度跃变最大值的4~5倍。

本发明只通过全色影像本身的信息来抑制标准矩匹配算法在去噪过程中产生的灰度畸变。该方法能够有效地复原影像、改善影像质量，运行效率也很高。

附图说明

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

具体实施方式

根据影像处理的基本理论可以知道，影像直方图的不同峰值代表了影像中包含的不同地物类型，标准矩匹配方法虽然在处理地物丰富的影像时会造成灰度畸变，即在影像直方图具有多个峰值或者直方图峰值曲线跨度较大的情况下，采用标准矩匹配方法会造成灰度畸变现象，但是对于灰度范围跨度较小的单个峰值的曲线分布，标准矩匹配方法在对影像进行噪声去除处理时，并不会造成影像中产生灰度畸变现象。因此将单一探元内部的成像数据根据影像直方图分布的差异(即不同的灰度等级)进行区分，然后分别求取不同灰度区域成像数据的方差和均值，最后对不同的灰度区域的成像数据，分别进行标准矩匹配处理，这样就可以解决标准矩匹配由于均值改变导致的探元间的灰度畸变以及由于方差改变导致的探元内的灰度畸变。

以下结合附图和实施例详细说明本发明技术方案。

基于灰度分割的自适应矩匹配条带噪声去除方法，其原理公式和处理步骤如下所示：

$y=\frac{\underset{m=0}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\σ}}_{\mathrm{mr}}}{\underset{m=0}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\σ}}_{\mathrm{mi}}}x+\underset{m=0}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{mr}}-\underset{m=0}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{mi}}\frac{{\mathrm{\σ}}_{\mathrm{mr}}}{{\mathrm{\σ}}_{\mathrm{mi}}}---\left(1\right)$

式(1)中，x，y分别为影像第i个扫描列各像素校正前、后的灰度值；σ_mr，μ_mr为灰度等级为m时，无噪声参考列的方差和均值；σ_mi，μ_mi为灰度等级为m时第i个扫描列的方差和均值。m的取值为0,1…M，灰度等级共有M+1个。对于每个灰度等级的成像数据，并非所有的探元的数据都需要参与相应的校正计算，只有当卫星载荷对应的探元在扫描成像时包含有对应灰度等级的成像数据，才参与相应的校正计算。

本发明技术方案可由本领域技术人员采用计算机软件技术实现自动运行流程，关键在于分割处理，即处理单元中心的自适应获取。实施例规定灰度等级划分的原则如下：

1)处理单元的左边取本单元中心和左边相邻中心的均值为界限，处理单元的右边取本单元中心和右边相邻中心的均值为界限。

2)灰度级跨度范围应小于某一阈值，以免处理单元的灰度跨度较大造成标准矩匹配处理时产生灰度畸变现象；

3)依据影像直方图，处理单元中心对应的灰度值所包含的像素个数应该在单元内处最大值，且要大于一个域值以保证处理单元中大部分是地物反射成分而不是噪声成分，从而达到去噪的目的而并非噪声增强处理；

4)处理单元的划分应该自动收敛，使得处理无需人工干预。

获取处理单元的关键在于准确定位直方图的峰值点，为了提高直方图峰值点获取的准确性，本发明实施例采用多层滤波的方式，并设置限定条件，最后对于灰度跨度范围依据噪声的特点进行限定，从而在去除条带噪声的同时抑制标准矩匹配算法带来的灰度畸变现象。

通过对国产卫星全色影像中残留的条带噪声进行分析发现，随机条带噪声和其四周影像的能量相关性较高，因此随机条带噪声也可以根据成像数据灰度等级区分方式，被分配到不同的灰度等级之中，根据这个原理，本发明实施例的处理步骤如下所示，如图1：

步骤1，获取原始含噪影像的亮度值曲线。

实施例中，假设一幅影像的大小为M×N，(i,j)表示影像中第i行第j列的像素，其中i的取值范围为0,1…M-1，j的取值范围为0,1…N-1。那么亮度值分别的曲线根据如下公式计算得到：

$S\left(k\right)=\left\{\underset{i=0,j=0}{\overset{i=M-1,j=N-1}{\mathrm{\Σ}}}{I}_E(i,j)\right|I_E(i,j)=k\}---\left(2\right)$

式(2)中，k表示对应的灰度级，I_E(i,j)表示影像中位于图像(i,j)处像素的灰度值。亮度值曲线S(k)通常由两部分组成：变化较平缓的基线和一些峰值，这些变换较平缓的基线说明了影像包含地物的丰富程度；而一些峰值则说明了影像中一些反映纹理特征的地物类型,由于随机条带噪声和地物类型存在高非线性相关性，因此随机条带噪声的成分随机分布于这些峰值周围，它们共同组成影像中地表地物反射率。

步骤2，通过采用不同大小的处理窗口对步骤1所得亮度值曲线进行中值滤波去除基线，获取处理单元中心。

依据步骤1所得亮度值曲线S(k)的特点可知，获取处理单元中心首先要进行亮度值曲线中的基线成分去除处理。本实施例采用中值滤波的方法来求取基线，事实证明它能有效的抑制峰值的影响.为了去掉亮度值曲线中的基线成分，使用一个大小为W_F的窗口作为中值滤波的处理窗口，输出的去基线的亮度值曲线S(k)_f为：

S(k)_f＝S(k)-med(S(k-t)|t∈W)                                (3)

式(3)中集合W＝{t|-(W_F-1)/2≤t≤(W_F-1)/2}，med是指中值滤波处理。

将去基线的亮度值曲线S(k)_f采用均值μ(k)_f和方差σ(k)_f进行归一化处理，处理公式如式(4)所示：

$S{\left(k\right)}_n=\frac{S{\left(k\right)}_f-\mathrm{\μ}{\left(k\right)}_f}{\mathrm{\σ}{\left(k\right)}_f}---\left(4\right)$

式(4)中S(k)_n就是均一化后的结果，其中μ(k)_f可以假设为0，因为S(k)_f经过了去基线处理。这里设定一个阈值T，那么选取的处理中心的候选值为：

$P{\left(k\right)}_c=\left\{\begin{array}{cc}1& \left|S{\left(k\right)}_n\right|>T\\ 0& \end{array}\right.---\left(5\right)$

这里T的选取比较重要，如果T选择大了，那么一些峰值点将会被排除在外，如果T选择小了，那么一些非峰值点将会被囊括进来，为了确保T的合理性，本发明建议选取最大|S(k)_n|的一半作为阈值。

上述步骤中中值滤波的窗口大小如果选择过大，那么一些峰值点将被划分到基线中去，如果窗口大小选择太小，那么一些非峰值点会被囊括到最后的结果当中，为了保证选取峰值点的最优结果，本算法中将窗口大小采用逐步递减的方式来获取多组去除基线的亮度值曲线，即获取多组处理中心的候选值，其中窗口大小W_F=21,20,19….5,4,3。每一次迭代时，若第k个灰度级是峰值点，那么P(k)_c就加1，本算法的迭代次数记为Count，其中本实施例中Count=21-3+1=19。因此，自动获取处理单元中心的迭代过程可以为，进行中值滤波，若窗口还未达到最小，继续减小窗口返回进行种植滤波，直到迭代完毕，每次迭代都可以得到可能作为处理单元中心的峰值点。

迭代结束后对获取的所有峰值点进行判定，判断其是否为处理中心，其判定函数公式如式(6)所示：

$P\left(k\right)=\left\{\begin{array}{cc}1& P{\left(k\right)}_c>\mathrm{Count}/2\\ 0& \end{array}\right.---\left(6\right)$

将求取的峰值点函数P(k)中的峰值带入到第一步中获取的亮度值曲线当中，若其对应的像素个数大于某一个域值I，那么保留这个处理中心，否则当作伪中心点予以去除，其中域值I的求取公式如下所示：

$I=\frac{1}{2^m-1}---\left(7\right)$

式(7)中，m表示原始含噪影像灰度值量化的位数。

具体实施时，可以采用窗口是否达到最小或迭代次数等作为迭代结束条件。实施例采用以下流程实现本步骤：

步骤2.1，设中值滤波的处理窗口大小记为W_F，初始化设置W_F的值（实施例设W_F=21），设迭代次数Count=0；

步骤2.2，获取去基线的亮度值曲线S(k)_f；

步骤2.3，对去基线的亮度值曲线S(k)_f采用均值μ(k)_f和方差σ(k)_f进行归一化处理，得到均一化后的结果S(k)_n；

步骤2.4，根据预设的阈值T选择候选处理中心，当去基线的亮度值曲线S(k)_f上某点的候选值P(k)_c为1时，该点是候选处理中心，否则不是；

步骤2.5，令Count=Count+1，判断Count是否达到预设次数（实施例设为19），否则进入步骤2.6，是则减小W_F的值（实施例设W_F=W_F-1），返回步骤2.2迭代执行中值滤波；

步骤2.6，对迭代所得的所有候选处理中心进行过滤，包括首先计算判断函数值P(k)，再将某候选处理中心的判断函数值P(k)带入到步骤1所得亮度值曲线S(k)中，若对应的像素个数大于预设阈值I，则确定为处理单元中心，否则去除。

步骤3，根据两个相邻处理单元中心的间距进行调整，确定处理单元。

实施例求取所有相邻处理单元中心之间的距离，通过比较相邻处理单元中心的距离与域值h的大小，对相邻处理单元进行合并和增加处理，获取最后的处理单元。实际处理过程中，根据国产卫星随机噪声特点，一般设定域值h的值为图像中随机噪声灰度跃变最大值的4~5倍。在处理过程中，当相邻处理单元中心的距离大于h时，为了确保处理单元灰度跨度的要求，需要对处理单元中心进行增加处理，即将相邻处理单元中心值取平均来获取新的处理中心，进而获取新的处理单元；当相邻处理单元中心的距离小于h时，可以将两个相邻处理单元进行合并作为新的处理单元，其中心值为两个相邻处理单元中心的均值。

步骤4，对单个处理单元采用标准矩匹配的处理方式进行影像处理。

根据现有标准矩匹配处理技术，实施例中标准矩匹配处理的公式如式(8)所示：

$y=\frac{{\mathrm{\σ}}_r}{{\mathrm{\σ}}_i}x+{\mathrm{\μ}}_r-{\mathrm{\μ}}_i\frac{{\mathrm{\σ}}_r}{{\mathrm{\σ}}_i}---\left(8\right)$

式(8)中，x，y分别为影像第i个扫描列中各像素校正前、后的灰度值；σ_r，μ_r为无噪声参考列的方差和均值；σ_r，μ_r为第i个扫描列的方差和均值。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (6)

1.一种基于灰度分割的自适应矩匹配条带噪声去除方法，其特征在于，包括以下步骤：

步骤1，获取原始含噪影像的亮度值曲线；

步骤2，通过采用不同大小的处理窗口对步骤1所得亮度值曲线进行中值滤波去除基线，获取处理单元中心；

步骤3，根据两个相邻处理单元中心的间距进行调整，确定处理单元；

步骤4，对单个处理单元采用标准矩匹配的处理方式进行影像处理。

2.如权利要求1所述基于灰度分割的自适应矩匹配条带噪声去除方法，其特征在于：设原始含噪影像的大小为M×N，(i,j)表示原始含噪影像中第i行第j列的像素，其中i的取值范围为0,1…M-1，j的取值范围为0,1…N-1，步骤1中获取原始含噪影像的亮度值曲线是根据如下公式计算得到，

其中，k表示对应的灰度级，I_E(i,j)表示影像中位于原始含噪影像(i,j)处像素的灰度值。

3.如权利要求2所述基于灰度分割的自适应矩匹配条带噪声去除方法，其特征在于：步骤2包括如下子步骤，

步骤2.1，设中值滤波的处理窗口大小记为W_F，初始化设置W_F的值，设迭代次数Count=0；

步骤2.2，执行以下计算公式获取去基线的亮度值曲线S(k)_f；

S(k)_f=S(k)-med(S(k-t)|t∈W)

其中，W={t|-(W_F-1)/2≤t≤(W_F-1)/2}，med是指中值滤波处理；

步骤2.3，对去基线的亮度值曲线S(k)_f采用均值μ(k)_f和方差σ(k)_f进行归一化处理，处理公式如下，

其中，S(k)_n就是均一化后的结果； 

步骤2.4，根据预设的阈值T选择候选处理中心，选择公式如下，

当去基线的亮度值曲线S(k)_f上某点的候选值P(k)_c为1时，该点是候选处理中心，否则不是；

步骤2.5，令Count=Count+1，判断Count是否达到预设次数，否则进入步骤2.6，是则减小W_F的值，返回步骤2.2迭代执行中值滤波；

步骤2.6，对迭代所得的所有候选处理中心进行过滤，

首先计算判断函数值，计算公式如下，

将某候选处理中心的判断函数值P(k)带入到步骤1所得亮度值曲线S(k)中，若对应的像素个数大于预设阈值I，则确定为处理单元中心，否则去除。

4.如权利要求3所述基于灰度分割的自适应矩匹配条带噪声去除方法，其特征在于：预设阈值I采用以下公式计算，

其中，m表示影像量化的位数。

5.如权利要求1或2或3或4所述基于灰度分割的自适应矩匹配条带噪声去除方法，其特征在于：步骤3中对处理单元中心进行调整，实现方式为，

当两个相邻处理单元中心的间距小于预设阈值h时，则将两个相邻处理单元中心合并为一个处理单元中心，合并后的处理单元中心的值为这两个相邻处理单元中心的均值；

当两个相邻处理单元中心的间距大于或等于预设阈值h时，两个相邻处理单元中心之间增加一个新的处理单元中心，新的处理单元中心的值为原始两相邻处理单元中心的平均值。

6.如权利要求5所述基于灰度分割的自适应矩匹配条带噪声去除方法，其特征在于：预设阈值h采用原始含噪影像的随机噪声灰度跃变最大值的4～5倍。 

Images