小波变换和局部插值融合的高光谱图像条带噪声消除方法与流程

本发明属于高光谱信号与信息处理技术领域，具体涉及小波变换和局部插值融合的高光谱图像条带噪声消除方法。

背景技术：

条带噪声普遍存在于高光谱遥感图像中，给高光谱图像的应用带来了极大挑战，有效去除条带噪声是高光谱图像很好应用的前提和基础。高光谱图像中条带噪声去除涉及到地物目标细节和边缘信息的保存，它们是一对矛盾的关系。高光谱图像中条带噪声的消除通常有两种思路，即硬件消除和软件消除。如果从硬件的角度进行消除，则需要改变现有传感器栅格阵列的成像机理，从根源消除条带噪声。但是，这种技术难度大，短时间内无法完成，而且工艺成本非常高。所以，在短时间内从硬件的角度彻底消除条带噪声的可能性不大，对于条带噪声的消除主要集中在软件上来实现。从高光谱遥感图像数据自身的角度进行处理来消除条带噪声的方法称为软件消除法。针对高光谱遥感图像数据的不同应用目的，不同的学者提出了许多处理方法，这些方法大概可以归纳成四种类型。第一种类型是传统的统计匹配方法，典型的代表是直方图匹配法和矩匹配法，这类方法容易实现，也不复杂，运算时间短，但对高光谱图像有较高的要求，要求高光谱图像中地物目标的分布比较均匀，否则处理的效果不佳，因此，对于地物分布差异大的高光谱遥感图像，这两种传统的方法已很难获得满意的消条效果，但是在它们的基础上，却有许多改进的方法被提出和应用；第二类是基于空间域的滤波方法，该类方法比较简单，容易实现，但是单个的空间域滤波方法在高光谱图像条带噪声滤除方面已显得无能为力，因为条带噪声去除效果不仅不太理想，而且残留噪声比较多，故单个滤波算法己经很少使用；第三类是基于变换域的滤波方法，这类的典型方法有傅里叶变换法和小波变换法等，由于傅里叶变换滤波法和其它的低通滤波一样，往往把条带噪声当作周期信号处理，这样很难找到一个合适的频率成分将条带噪声和信号进行完全分离，这类方法不仅不能彻底去除条带噪声，而且计算复杂，性能不稳定，同时还会造成信号损失，使图像细节模糊，小波变换理论用于图像滤波处理，主要有两种方式进行，即软阈值滤波和硬阈值滤波，这两种方式对于条带噪声效果也不理想，于是有学者提出利用小波分解系数归零的思路来进行处理，对于细条带噪声而言，效果还可以，但是对于宽条带噪声，滤波效果就不能让人满意了；第四种类型是基于一些新的理论和方法，如最近提出的变分理论、相位一致性、正规化低秩表示法等等，每种方法通常都有其针对性，在一定范围内都能去掉条带噪声，但是它们的普适性和通用性比较差，这是因为不同的高光谱图像具有不同的条带噪声特征，例如亮条带、灰色条带和暗条带，细条带、宽条带和超宽条带等，并且出现的位置也不是固定的，所以，如果要想有效消除噪声的影响，必须先对高光谱图像中条带噪声的产生机理和分布特征进行分析和研究，否则滤波效果很难达到预期目标。

通过对现有方法进行深入分析和研究，发现这些方法虽然能够去除高光谱遥感图像中的条带噪声，主要针对的细条带噪声和比较均匀的条带噪声，并且去除效果并不很理想，有时图像的清晰度不是很高，同时图像中的一些细节信息也经常被损失。例如，小波变换用于条带噪声的去除，不论是传统的阈值法还是分解系数归零法，都会模糊图像的一些几何细节信息，当然能达到去条带噪声的目的。对于宽带和超宽带以及分布随机的条带噪声，用传统的单个滤波方法，已经很难获得满意的结果。

技术实现要素：

为了克服上述现有技术的缺点，本发明的目的在于提供小波变换和局部插值融合的高光谱图像条带噪声消除方法，不仅能有效地消除一般的细条带噪声，而且对宽带和超宽带以及分布无规律的条带噪声均能起到消除作用，还能较好地保持图像中的一些细节和边缘几何信息，同时还能应用于不同的高光谱图像，有着较好的普适性。

为了达到上述目的，本发明采取的技术方案为：

小波变换和局部插值融合的高光谱图像条带噪声消除方法，包括以下步骤：

步骤1：输入高光谱图像数据立方体，高光谱图像数据立方体为经过预处理的2级产品数据，预处理包括随机噪声滤波、辐射校正和几何校正处理；

步骤2：对高光谱图像波段数据进行判断：高光谱图像数据在获取过程中受到各种因素的干扰，会使部分波段数据的质量受到损坏，损坏的波段图像就无需继续进行处理；

步骤3：剔除损坏的高光谱图像波段数据；

步骤4：读取正常的第k波段的高光谱图像数据：在剩下的高光谱图像数据中选择第k波段数据，分别进行步骤5和步骤6操作；

步骤5：用二维平稳小波进行高光谱图像数据的条带噪声去除；

步骤6：用局部插值进行高光谱图像数据的条带噪声滤除；

步骤7：对步骤5和步骤6的处理结果进行融合处理；

步骤8：判断处理的k波段是否为最后的波段，如果是最后波段，则结束；如果不是最后波段，则转入步骤4，即读取下一波段的图像。

所述的步骤2中高光谱图像波段数据判断的方法为：

采用列平均值的导数(DMC)和列平均值的梯度(GMC)的乘积δ(k)来进行判断；

δ(k)＝D(k)×G(k) (1)

列平均值的导数(DMC)和列平均值的梯度(GMC)的数学模型如下：

其中，k表示第k波段，D(k)和G(k)分别表示第k波段列平均值的导数(DMC)的导数平方和及列平均值的梯度(GMC)的平方和，μ_i(k)和分别表示第k波段图像的第i列探测单元的平均值和梯度值，δ(k)表示它们的乘积。

所述的步骤3中剔除损坏的高光谱图像波段数据的方法为：如果δ(k)值比预定的阈值ε大，即δ(k)＞ε，表明该波段的图像数据已受到损坏，直接删除该波段的数据。

所述的步骤5中用二维平稳小波进行条带噪声去除的具体步骤为：

步骤5.1：条带噪声类型转换：先把亮条带转变成暗条带；

步骤5.2：设置小波分解尺度数：二维平稳小波分解的尺度数设为5，即M＝5；

步骤5.3：选择小波函数：选择Daubechies(db1)函数作为小波分解的基函数，然后选择的k波段图像进行分解；

步骤5.4：对图像进行分解：利用步骤5.2和步骤5.3设置的参数，对输入的高光谱k波段图像进行小波分解；

步骤5.5：对小波分解后各尺度的分解系数进行处理：首先判断系数子图像是否还包含大量的条带噪声，如果系数子图像中包含的条带噪声较多，则转至步骤5.2，对该系数在图像继续进行分解，直到当某个分解系数子图像几乎是条带噪声时停止进行小波分解；然后对最后的分解层获得的系数值全部进行赋零值处理；最后逐系数逐尺度进行小波逆变换，最终得到各个分解尺度的小波系数；

步骤5.6：用步骤5.5得到的小波系数进行整个图像的小波逆变换处理，获得用小波变换处理条带噪声后的图像。

所述的步骤6中局部插值法去除条带噪声的具体步骤为：

步骤6.1：计算高光谱图像的列均值，并获得列均值曲线图；

步骤6.2：计算列均值的梯度，并获得列均值的梯度曲线图；

步骤6.3：通过梯度曲线图，获得波峰和波谷所在的位置，这就是条带噪声所在的位置；

步骤6.4：把条带噪声所在的列进行去除，然后进行插值处理，结果是获得了局部插值的图像；

步骤6.5：进行矩匹配处理：式(4)就是矩匹配的数学模型，

其中，k表示高光谱图像中第k波段图像，X_k(i,j)和Y_k(i,j)分别表示第k波段图像单个像素灰度值调整前后值，μ_kr和σ_kr分别表示整幅图像的均值和标准差，μ_kj和σ_kj分别为第j列的均值和标准差。

所述的步骤7中对步骤5和步骤6的处理结果进行融合处理的具体步骤为：

融合规则如下式：

I_f(k)＝a×I_wt(k)+b×I_li(k) (5)

式中的k表示第k波段，才表示融合后的图像，I_wt(k)表示小波变换去条带噪声后的图像，I_wt(k)表示局部区域插值去条带噪声后的图像，参数a和b分别为各自的权系数，它们的取值范围为[0,1]，并且a+b＝1。

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

(1)方法设计的优点：在处理方法上，首先设计了两次去除条带噪声，但是，这两次条带噪声的去除不是前后级联的关系，而是并列的关系。小波变换具有多尺度多方向的特点，而条带噪声有很好的方向性，因此小波变换能有效提取到条带噪声，可是同时也去掉了不少细节信息。通过寻找噪声所在的局部位置，然后进行去除，随后进行局部插值，这样保护了其它细节信息。再把两者的结果进行融合，达到了既去除条带又保护细节的目的。其次，用小波进行高光谱图像分解时，选择的是二维平稳小波，而不是二维离散小波。因为二维离散小波是下2采用分解，尺度每增加1，分解后的系数子图像变为上层的1/2，不便于对图像进行处理。采用二维平稳小波分解，每个尺度上的子图像，每个子图像分解后的子图像，都与原始图像保持相同大小，这样便于处理和小波逆变换。

(2)应用上的优点：条带噪声是高光谱成像机理产生，在硬件和工艺上很难消除，因此只能采用数据处理的方法进行去除。一般的方法均是针对不同的应用目标，适应性非常有限。本方法有较强的普适性，可用于不同传感器的高光谱图像条带噪声的消除，而且可同时去除细条带、宽条带、超宽条带，或者灰色条带、暗条带和亮条带。不仅能有效消除条带，还能较好地保持细节信息。

(3)数据处理的优点：高光谱图像是数据立方体，数据量非常大，为了提高处理速度，减少不必要的处理时间，本方法对数据进行了优化处理。首先，判断哪些数据是有用的，哪些数据已经受到损坏，对损坏的数据图像进行剔除，然后保留好的数据图像进入下一步处理，这样节省了对坏数据图像处理的时间。此外，为便于小波分解提取条带噪声，本方法对条带噪声类型进行了处理，有利于提高去条带的效果，例如亮条带变成暗条带。第三，最后对不同思路的处理结果进行了融合处理，使处理结果达到保持细节信息的同时有效滤除条带噪声。

附图说明

图1为本发明的流程图。

图2(a)为实施例原始的高光谱图像数据，图2(b)、图2(c)、图2(d)分别为D(k)、G(k)和δ(k)值的曲线图，横坐标为波段数，纵坐标为数值。

图3为实施例原始的高光谱图像，其中图a为6波段，图b为23波段，图c为27波段，图d为52波段。

图4为利用二维平稳小波变换对高光谱图像去除条带后的结果，其中图a为6波段，图b为23波段，图c为27波段，图d为52波段。

图5为图3(b)的梯度曲线图。

图6为局部插值法处理过程结果，其中图a为直接插值，图b为矩匹配。

图7为局部插值法消噪结果，其中图a为6波段，图b为23波段，图c为27波段，图d为52波段。

图8为矩匹配消噪结果，其中图a为6波段，图b为23波段，图c为27波段，图d为52波段。

图9为本发明方法消噪结果，其中图a为6波段，图b为23波段，图c为27波段，图d为52波段。

具体实施方式

下面结合附图和实施例对本发明作进一步阐述。

本实施例的图像截取大小为512×486，空间分辨率为20m，光谱分辨率为20nm，成像波长范围为1100-2500nm。

参照图1，小波变换和局部插值融合的高光谱图像条带噪声消除方法，包括以下步骤：

步骤1：输入高光谱图像数据立方体，高光谱图像数据立方体为经过预处理的2级产品数据，预处理包括随机噪声滤波、辐射校正和几何校正处理；所有的波段有75个，成像频谱为1100纳米到2500纳米，即短波红外成像区域；

步骤2：对高光谱图像波段数据进行判断：高光谱图像数据在获取过程中受到各种因素的干扰，会使部分波段数据的质量受到损坏，损坏的波段图像就无需继续进行处理；

采用列平均值的导数(DMC)和列平均值的梯度(GMC)的乘积δ(k)来进行判断；

δ(k)＝D(k)×G(k) (1)

列平均值的导数(DMC)和列平均值的梯度(GMC)的数学模型如下：

其中，k表示第k波段，D(k)和G(k)分别表示第k波段列平均值的导数(DMC)的导数平方和及列平均值的梯度(GMC)的平方和，μ_i(k)和分别表示第k波段图像的第i列探测单元的平均值和梯度值，δ(k)表示它们的乘积，如图2所示；

步骤3：剔除损坏的高光谱图像波段数据：如果δ(k)值比预定的阈值ε大，即δ(k)＞ε，表明该波段的图像数据已受到损坏，没有必要进入下一步操作，直接删除该波段的数据；

虽然D(k)和G(k)可以单独来判断高光谱图像的质量问题，从图2(b)和图2(c)中可知，波段1-4,19-22,38-42,67-75的数据有问题，但是在波段43-67之间时，它们很难判断，因为曲线的波动性很大。在图2(b)中可能会判断波段27有问题，实际上没问题，利用δ(k)可以较准确的判断哪些波段的数据有问题，波动性越大，数据越有问题，对于这些有问题的数据，没有必要进行下一步的处理。为了便于比较分析，选择了第6、23、27和52波段数据，分别如图3(a)、图3(b)、图3(c)、图3(d)所示，原因有两个：一是图像数据没有受到损坏；二是里面包含了不同类型的宽条带噪声，在图3中，为便于条带噪声的消除和视觉效果检测，对原始图像进行了转置保存，即把横向的条带噪声转置成垂直向条带噪声；

步骤4：读取正常的第k波段的高光谱图像数据：在剩下的高光谱图像数据中选择第k波段数据，分别进行步骤5和步骤6操作；

步骤5：用二维平稳小波进行高光谱图像数据的条带噪声去除：

步骤5.1：条带噪声类型转换：由于处理的图像是2级产品，已经过一系列的处理，如随机噪声滤波、辐射校正和几何校正等，当然也对细条带进行了处理，但是条带噪声仍然存在，特别是一些宽/粗的条带噪声以及超宽条带。这里选择的是平稳小波变换，利于条带噪声的去除处理，在小波分解中，亮条带噪声和几何细节信息很难区分提取，它们通常存在于高频子系数图像中。为了尽可能消除条带噪声，先把亮条带转变成暗条带，这样在小波分解后，条带噪声主要存在于高尺度分解层，而几何细节信息主要存在于低尺度分解系数中；

步骤5.2：设置小波分解尺度数：二维平稳小波分解的尺度数设为5，即M＝5；

步骤5.3：选择小波函数：选择Daubechies(db1)函数作为小波分解的基函数，然后选择的k波段图像进行分解；

步骤5.4：对图像进行分解：利用步骤5.2和步骤5.3设置的参数，对输入的高光谱k波段图像进行小波分解；

步骤5.5：对小波分解后各尺度的分解系数进行处理：首先判断系数子图像是否还包含大量的条带噪声，如果系数子图像中包含的条带噪声较多，则转至步骤5.2，对该系数在图像继续进行分解，直到当某个分解系数子图像几乎是条带噪声时停止进行小波分解；然后对最后的分解层获得的系数值全部进行赋零值处理；最后逐系数逐尺度进行小波逆变换，最终得到各个分解尺度的小波系数；

步骤5.6：用步骤5.5得到的小波系数进行整个图像的小波逆变换处理，获得用小波变换处理条带噪声后的图像；

如图4所示，图4中(a)-(d)对应图3(a)-(d)，从图4中可以看出，对于亮条带和灰度条带噪声，小波变换能有效去除，但是对于暗条带，小波变换不易去除。因为暗条带与图像的背景信息类似，去除了暗条带也去除了背景；

步骤6：用局部插值进行高光谱图像数据的条带噪声滤除；

步骤6.1：计算高光谱图像的列均值，并获得列均值曲线图；

步骤6.2：计算列均值的梯度，并获得列均值的梯度曲线图；

步骤6.3：通过梯度曲线图，获得波峰和波谷所在的位置，这就是条带噪声所在的位置；

步骤6.4：把条带噪声所在的列进行滤除，然后进行插值处理，结果是获得了局部插值的图像，因为只对局部区域进行处理，保留了绝大部分的几何细节信息不变；

步骤6.5：进行矩匹配处理：矩匹配方法法虽然是一种经典的条带噪声去除算法，但是单独使用时，要求地物类型分别均匀，故当地物分布非均匀时，滤波效果不太理想，如果要取得好的滤波效果，加大矩匹配算法的应用范围，需要不同的方法进行组合，式(4)就是矩匹配的数学模型，

其中，k表示高光谱图像中第k波段图像，X_k(i,j)和Y_k(i,j)分别表示第k波段图像单个像素灰度值调整前后值，μ_kr和σ_kr分别表示整幅图像的均值和标准差，μ_kj和σ_kj分别为第j列的均值和标准差；

以图3(b)为例进行说明，结果如图5所示，有5条明显的条带，而且条带的宽度都超过一个像素，从图5中可以清楚地看到条带噪声所在的位置。找到条带噪声位置后，接着是去除条带噪声和插值处理，结果如图6(a)所示，效果不是很好。对于窄条带，局部插值还可以，特别单个像素的条带，效果尤佳，随着条带宽度的增加，局部插值的去噪能力逐渐减弱，如图6(a)所示，在插值的基础上再进行矩匹配处理，效果得到较明显的改善，可以把超宽带噪声去除，如图6(b)所示。虽然效果还不是十分的理想，但是可以接受。如图7所示，图7(a)-(d)分别对应图3(a)-(d)的处理结果。从图7可以看到，局部插值法基本上能去除较宽的条带噪声，由于进行的是局部插值，所以细节信息也保护的比较好。

矩匹配方法实验结果如图8所示，从图8中可知，直接的矩匹配去条带噪声不彻底，留下的痕迹比较多，由于矩匹配的假设前提是地物分布类型相对均匀，对于分布不均匀的地物，其去噪效果不太好，因此，矩匹配很少单独使用；

步骤7：对步骤5和步骤6的处理结果进行融合处理；融合规则如下式：

I_f(k)＝a×I_wt(k)+b×I_li(k) (5)

式中的k表示第k波段，才表示融合后的图像，I_wt(k)表示小波变换去条带噪声后的图像，I_wt(k)表示局部区域插值去条带噪声后的图像，参数a和b分别为各自的权系数，它们的取值范围为[0,1]，并且a+b＝1，设置融合因子a＝b＝0.5；

步骤8：判断处理的k波段是否为最后的波段，如果是最后波段，则结束；如果不是最后波段，则转入步骤4，即读取下一波段的图像。

本发明方法的结果如图9所示，从图4、图7、图8和图9中可知，本发明方法的效果最好，达到了预期的目的。它不仅能有效去除条带噪声，还能保护更多的细节信息。