基于有序最小值与小波滤波的高光谱图像宽条带去除方法与流程

本发明属于信号与信息处理技术领域，具体涉及图像的有序滤波和小波变换多尺度分解滤波以及矩匹配滤波相结合的一种基于有序最小值与小波滤波的高光谱图像宽条带去除方法。

背景技术：

高光谱遥感技术诞生于20世纪80年度初，主要在紫外、可见光和红外波段中获得非常窄的波段图像，从而可获得反映地物化学材质特性的光谱曲线，使得许多原在多光谱图像中无法识别的特征变得非常容易实现。其中，短波红外波段的高光谱遥感图像是地物检测和识别的重要数据源，但是目前该波段的数据存在大量的随机噪声和条带噪声，特别是星载高光谱遥感图像，存在超宽的条带噪声，严重降低了高光谱成像质量，并影响了高光谱图像中地物光谱特征的提取与分析，给地物分类识别精度和有效探测带来了极大的挑战和困难。因此，对条带噪声的去除是高光谱遥感图像预处理的重要内容。

高光谱遥感在成像过程中会受到多噪声源的影响，而且引入的噪声还会与信号一起历经多种系统处理和转换，最终使得高光谱图像中的噪声分布相当复杂。同时，由高光谱成像机理产生的条带噪声也相当复杂，通常很难完全去除。因为这些条带噪声在不同波段的图像中所表现的形式是不完全相同，有些波段图像中会多一些，有些波段图像中会少一些，而且条带的宽度和散射强度也不固定。

条带噪声是由于高光谱成像传感器CCD阵列的响应函数不一致性造成的，而且每种传感器包含的CCD个数不同，品种不同，所以每种高光谱成像仪获得的遥感数据中的条带噪声也通常不相同。基于此，虽然许多学者针对高光谱图像中的条带噪声的成因和特点进行了大量的研究工作，提出了许多去条带噪声的相关理论和方法，但是不同方法基本上是针对特殊成像仪或图像，很难有通用性和普适性，也就是说都有一定的局限性和优点。最关键的是这些研究主要集中在原始遥感数据阶段，此阶段的噪声主要属于细条带噪声，而且有相应的周期和规律，比较好去除。此外，高光谱图像中还存在大量的宽带和超宽带条带噪声，用以往的方法很难去除宽条带噪声。

高光谱遥感图像中的条带噪声是由其成像机理产生的，而不是由外界因素造成的。虽然有关条带噪声滤除的方法有许多，例如，传统的直方图匹配法和矩匹配法及其大量的改进算法，空间域和变换域的滤波方法，基于变分理论和压缩感知理论的新方法等。但是，条带噪声是由CCD传感器之间的响应函数不一致导致的，而且不同的高光谱成像仪包含的传感器数量不同，使得现有的算法是为解决某类图像问题而提出，面对其它传感器图像，其效果并不理想。

目前，条带噪声的去除有两种策略，一是硬件去除，二是软件去除。由于硬件去除对传感器有较高的要求，实现难度比较大，所以通常采用软件的方法去除条带噪声。而单个的算法或方法很难获得理想的效果，例如，小波多尺度分解对于亮条带噪声处理的效果不理想，矩匹配要求地物分布比较均匀。

技术实现要素：

为了克服上述现有技术的不足，本发明的目的是提供一种基于有序最小值与小波滤波的高光谱图像宽条带去除方法，简称OWM算法，从高光谱成像机理和条带噪声产生机理以及条带噪声分布的特点出发，结合小波变换理论多尺度分解的特点和矩匹配的优势，该方法不仅能够有效去除暗条带、灰色条带和亮条带噪声，而且能对超宽条带进行滤除。本方法主要针对2级高光谱图像产品中的宽条带噪声，对细条带噪声其效果更好。在小波多尺度分解中，对于亮条带噪声，在各个分解尺度的分解系数中，其往往和图像中地物目标的几何细节信息混在一起，不易去除，小波分解只提供3个方向的分解系数，而像曲波分解和轮廓波分解能提供更多的方向的分解系数，反而不适合条带噪声的提取。经验模态分解方向性比较差，也不适合条带特征提取。条带噪声不同与随机噪声，它们产生的机理不同，表现的方式不同。因此，小波去除条带噪声的效果也很差。假设条带噪声是垂直分布的，用小波对含条带噪声的高光谱图像进行M尺度分解后，条带噪声主要存在于每个分解尺度上垂直子图像和低频近似子图像中，对这两个子图像进行判断，如果它们包含的信息几乎是条带噪声，把它们设置为零，如果不是，则对它们继续进行小波分解，直到包含的信息几乎是条带信息为止，然后进行系列的反变换，就初步完成了条带噪声的去除。对于小波分解理论，为便于处理和反变换，采用的是二维平稳小波，而不是二维离散小波。小波分解里面对于亮条带噪声，如果直接去除，会丢失许多几何细节信息，所以，通过最小值滤波，把亮条带转为暗条带，利用小波多尺度分解去除条带噪声充分利用了小波多尺度多方向性的特点和条带噪声具有方向性的优势相结合，进而有效去除，最后还要进行矩匹配滤波处理，而能获得好的去噪效果。

为了实现上述目的，本发明采用的技术方案是：

一种基于有序最小值与小波滤波的高光谱图像宽条带去除方法，包括以下步骤：

1)输入高光谱遥感图像；

2)对步骤1)输入高光谱图像进行对比度调整，首先获得输入图像灰度值的范围，然后把它调整到1～256范围内，调整对比度的目的是调整图像灰度值的动态范围，提高图像的辐射分辨率，同时亮条带更亮一些，暗条带更暗一些，有利于后续步骤对条带噪声的去除；

3)进行最小值滤波，对步骤2)调整过后的高光谱图像进行最小值滤波步骤如下：

高光谱图像中的条带噪声，有可能是明亮的条带，也有可能是暗淡的条带，还有可能是两者兼有，特别是明暗条带都存在时，比较难去除，在OWM算法中，选择最小值滤波，把亮的条带转变成较暗的条带，如果把暗条带转变成较亮条带，那么在小波变换处理中，这些条带噪声将和地物细节信息混合出现在低尺度的三个高频分解系数中，不利于条带噪声的去除；如果转成暗条带，被小波分解后，它们主要包含在每个分解尺度中的近似低频分解系数中。所以，OWM算法利用最小值的滤波原理把亮条带转变成暗条带，以图像中某点(x,y)为中心，设置一个滑动窗口，窗口大小为K×L，对窗口内的所有像素值与中心值I(x,y)进行比较，如果该值比中心像素值小，则进行替换，在本算法中，设置窗口参数K＝L＝3；

4)确定条带噪声的方向性；其具体步骤为：

判断输入图像中条带的方向性，首先假设的条带是垂直条带噪声，所以如果条带是垂直条带，则直接进入步骤6)；如果条带是横向条带，则进入步骤5)；

5)进行转置处理；

6)进行小波分解，去除条带噪声；其步骤为：

a、确定小波分解方式，选择二维平稳小波变换(SWT)对图像进行分解，而不是二维离散小波变换(DWT)，因为SWT变换是非下采样，分解后的所有系数子图像的大小跟原始图像的大小一样，这样有利于对分解系数进行处理，DWT是下2采样的，每次采样后系数子图像大小是上一尺度的一半，随着尺度的增加，系数子图像的大小越来越小，不便于对子图像进行处理；

b、确定小波基函数和分解尺度M，常用的小波函数即可，在OWM算法中选择dbN小波，小波分解的尺度不超过5，即M<＝5比较合适，分解尺度设为5，在实际应用过程中，根据条带噪声的出现情况对相应尺度进行处理；

c、确定需进行处理的系数子图像，小波变换对图像进行分解，每个分解尺度上可以获得四个子系数图像，包括一个近似低频图像和三个高频子图像，这三个高频子图像分别表示水平方向、垂直方向和对角方向的信息，由于处理的是垂直条带噪声，所以条带噪声信息主要存在于垂直系数子图像和近似系数子图像中，对这两个系数子图像进行处理即可；

7)对系数子图像进行SWT分解；其具体步骤为：

如果系数子图像中包含的信息除了条带外还有其它信息，就对该系数子图像进行小波分解，分解的尺度设为1或2即可，然后选择垂直子图像和低频近似子图像；

8)系数子图像进行归零处理；其具体步骤为：

把步骤7)中所获得的近似子图像和垂直子图像进行归零处理，即把该系数子图像的所有值设为零；

9)对处理过的系数子图像进行ISWT变换；

10)全部分解系数进行逆平稳小波变换(ISWT)；

11)进行矩匹配处理；其具体步骤为：

a、由于条带噪声出现的主要原因是各个CCD传感器响应函数的不一致，因此，矩匹配算法的思想是假设各CCD探元的响应函数为具有移不变性质的线性函数，令C_i为第i个CCD探元，则C_i的光谱响应函数可以用式(1)表示，即

Y_i＝k_iX+b_i+e_i(X) (1)

其中，Y_i为的输出值，也就是图像中像素的灰度值。X为该CCD记录的辐射值，k_i为响应函数的增益，b_i为偏移，e_i为随机噪声，如果图像的信噪比较高，那么式(1)中的随机噪声可以忽略，式(1)可以写为式(2)如下：

Y_i＝k_iX+b_i (2)

从式(2)可以看出，即使是同一辐射强度X，如果其增益k_i和偏移量b_i的取值不同，则会得到不同的结果，因而就产生了条带噪声，所以，根据条带噪声产生原理，如把不同的值归一化到相同的值，则可以有效地消除条带噪声，以某波段图像中的某CCD列作为参考列，利用公式将其它各CCD列的值校正到该参考CCD的辐射率上，就可实现条带噪声的消除，归一化匹配的数学模型如式(3)所示，即

其中，X和Y分别表示某波段图像的第i列像素校正前和校正后的灰度值，μ_r和σ_r分别是参考CCD列的均值和标准差，μ_i和σ_i分别为第i列的均值和标准差。参考值通常用整幅图像的均值和标准差代替参考列的均值和标准差；

步骤12：输出消除条带噪声的图像。

进一步，所述的步骤2)对高光谱遥感数据中的所有波段图像都进行处理，首先就某个波段的图像进行处理，获得该波段图像灰度值的最大值和最小值，确定新的灰度值范围，其范围一定要大于图像初始的灰度范围，实际上是一个拉伸的处理过程，然后，对所有波段的图像进行相同的处理，最后获得对比度增强后的高光谱图像数据。

进一步，所述的步骤3)对增强处理后的高光谱图像实施有序最小值滤波，以一个大小为K×L的滑动窗作为模板，在某个波段图像上进行移到，把滑动窗内的像素灰度值进行排序，选择最小的值作为模板中心像素的值，直到所有波段进行处理完毕，窗口大小一般根据条带宽带进行取值，如果条带噪声比较宽，K和L的取值相对大一点；如果条带较窄，K和L的值取相对小一些，通常情况下取7或9比较合适。

进一步，所述的步骤4)及步骤5)确定条带噪声的方向，并进行转置处理,，假设条带噪声是垂直条带，因此对于横向的条带噪声，先进行一个转置处理，变成垂直条带，进行转置处理之前，判断是否是横向还是纵向条带噪声。

进一步，所述的步骤6)对初步处理后的高光谱图像进行小波多尺度分解，再对各个分解尺度上的分解系数进行处理，达到消除条带噪声的目的，用小波对图像进行多尺度分解，图像被小波进行多尺度分解后，接着是对各个尺度上获得的分解系数进行处理，由于假设的条带噪声是垂直条带，因此除了最高的分解尺度外，其它每个分解层只要处理垂直高频系数子图像即可，判断垂直系数子图像中条带噪声的包含量，如果几乎是条带噪声，则直接把该系数值赋值为零，否则对该系数子图像继续进行小波分解，处理方式如前，直到某分解尺度上低频或垂直高频子图像几乎是条带噪声为止，处理完后，逐层进行小波逆变换，结果是经小波变换消条处理后的图像。

进一步，所述的步骤11)所述的进行矩匹配处理是对处理后的图像进行第二次消条处理。

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

本发明是一种新的、非常有效的条带噪声去除算法，特别是包含不同条带以及超宽条带的高光谱图像，其优势非常明显，具体体现在以下三个方面：

1)方法设计的优点：

设计上采用双重去条带噪声的思想和平稳小波分解的策略。由于单个的方法有较大的不足，要么去除条带噪声不好，要么普适性非常差。采用小波多尺度分解滤波处理和矩匹配处理来去除条带噪声，效果非常不错。利用二维平稳小波对图像进行处理，便于对小波分解系数的处理与重构。

2)数据处理的优点：

为了尽量去除条带噪声，对高光谱图像数据进行了对比度增强处理和最小值滤波处理。这样处理的目的或者结果是：把亮条带噪声彻底变成暗条带噪声，以便在后续的处理中小波多尺度变换能彻底消除条带噪声。

3)条带去除的优点：

用本方法高光谱遥感图像，不仅能有效消除不同宽度的条带噪声，例如，宽条带和超宽条带，而且能消除各类灰度条带噪声，如暗条带、灰条带和亮条带，使原来用一般方法无法消除的各种条带噪声能有效地被去除，表明该方法有较好的普适性。

附图说明

图1为本发明的方法流程示意图。

图2为实施例的第6、23、27和52波段图像。其中a为第6波段；b为第23波段；c为第27波段；d为第52波段。

图3为第6波段处理结果图，其中(A)是BLPF算法处理的结果，(B)是MMF算法获得的结果，(C)是WTZF算法获得的结果，(D)是OWM算法获得的实验结果。

图4为第23波段处理结果图，其中(A)是BLPF算法处理的结果，(B)是MMF算法获得的结果，(C)是WTZF算法获得的结果，(D)是OWM算法获得的实验结果。

图5为第27波段处理结果图，其中(A)是BLPF算法处理的结果，(B)是MMF算法获得的结果，(C)是WTZF算法获得的结果，(D)是OWM算法获得的实验结果。

图6为第52波段处理结果图，其中(A)是BLPF算法处理的结果，(B)是MMF算法获得的结果，(C)是WTZF算法获得的结果，(D)是OWM算法获得的实验结果。

具体实施方式

以下结合附图对本发明进一步叙述，但本发明不局限于以下实施例。

如图1所示，一种基于有序最小值与小波滤波的高光谱图像宽条带去除方法，其特征在于，包括以下步骤：

1)输入高光谱遥感图像；

2)对步骤1)输入高光谱图像进行对比度调整，首先获得输入图像灰度值的范围，然后把它调整到1～256范围内，调整对比度的目的是调整图像灰度值的动态范围，提高图像的辐射分辨率，同时亮条带更亮一些，暗条带更暗一些，有利于后续步骤对条带噪声的去除；

3)进行最小值滤波，对步骤2)调整过后的高光谱图像进行最小值滤波步骤如下：

高光谱图像中的条带噪声，有可能是明亮的条带，也有可能是暗淡的条带，还有可能是两者兼有，特别是明暗条带都存在时，比较难去除，在OWM算法中，选择最小值滤波，把亮的条带转变成较暗的条带，如果把暗条带转变成较亮条带，那么在小波变换处理中，这些条带噪声将和地物细节信息混合出现在低尺度的三个高频分解系数中，不利于条带噪声的去除；如果转成暗条带，被小波分解后，它们主要包含在每个分解尺度中的近似低频分解系数中。所以，OWM算法利用最小值的滤波原理把亮条带转变成暗条带，以图像中某点(x,y)为中心，设置一个滑动窗口，窗口大小为K×L，对窗口内的所有像素值与中心值I(x,y)进行比较，如果该值比中心像素值小，则进行替换，在本算法中，设置窗口参数K＝L＝3；

4)确定条带噪声的方向性；其具体步骤为：

判断输入图像中条带的方向性，首先假设的条带是垂直条带噪声，所以如果条带是垂直条带，则直接进入步骤6)；如果条带是横向条带，则进入步骤5)；

5)进行转置处理；

6)进行小波分解，去除条带噪声；其步骤为：

a、确定小波分解方式，选择二维平稳小波变换(SWT)对图像进行分解，而不是二维离散小波变换(DWT)，因为SWT变换是非下采样，分解后的所有系数子图像的大小跟原始图像的大小一样，这样有利于对分解系数进行处理，DWT是下2采样的，每次采样后系数子图像大小是上一尺度的一半，随着尺度的增加，系数子图像的大小越来越小，不便于对子图像进行处理；

b、确定小波基函数和分解尺度M，常用的小波函数即可，在OWM算法中选择dbN小波，小波分解的尺度不超过5，即M<＝5比较合适，分解尺度设为5，在实际应用过程中，根据条带噪声的出现情况对相应尺度进行处理；

c、确定需进行处理的系数子图像，小波变换对图像进行分解，每个分解尺度上可以获得四个子系数图像，包括一个近似低频图像和三个高频子图像，这三个高频子图像分别表示水平方向、垂直方向和对角方向的信息，由于处理的是垂直条带噪声，所以条带噪声信息主要存在于垂直系数子图像和近似系数子图像中，对这两个系数子图像进行处理即可；

7)对系数子图像进行SWT分解；其具体步骤为：

如果系数子图像中包含的信息除了条带外还有其它信息，就对该系数子图像进行小波分解，分解的尺度设为1或2即可，然后选择垂直子图像和低频近似子图像；

8)系数子图像进行归零处理；其具体步骤为：

把步骤7)中所获得的近似子图像和垂直子图像进行归零处理，即把该系数子图像的所有值设为零；

9)对处理过的系数子图像进行ISWT变换；

10)全部分解系数进行逆平稳小波变换(ISWT)；

11)进行矩匹配处理；其具体步骤为：

a、由于条带噪声出现的主要原因是各个CCD传感器响应函数的不一致，因此，矩匹配算法的思想是假设各CCD探元的响应函数为具有移不变性质的线性函数，令C_i为第i个CCD探元，则C_i的光谱响应函数可以用式(1)表示，即

Y_i＝k_iX+b_i+e_i(X) (1)

其中，Y_i为的输出值，也就是图像中像素的灰度值。X为该CCD记录的辐射值，k_i为响应函数的增益，b_i为偏移，e_i为随机噪声，如果图像的信噪比较高，那么式(1)中的随机噪声可以忽略，式(1)可以写为式(2)如下：

Y_i＝k_iX+b_i (2)

从式(2)可以看出，即使是同一辐射强度X，如果其增益k_i和偏移量b_i的取值不同，则会得到不同的结果，因而就产生了条带噪声，所以，根据条带噪声产生原理，如把不同的值归一化到相同的值，则可以有效地消除条带噪声，以某波段图像中的某CCD列作为参考列，利用公式将其它各CCD列的值校正到该参考CCD的辐射率上，就可实现条带噪声的消除，归一化匹配的数学模型如式(3)所示，即

其中，X和Y分别表示某波段图像的第i列像素校正前和校正后的灰度值，μ_r和σ_r分别是参考CCD列的均值和标准差，μ_i和σ_i分别为第i列的均值和标准差。参考值通常用整幅图像的均值和标准差代替参考列的均值和标准差；

步骤12：输出消除条带噪声的图像。

进一步，所述的步骤2)对高光谱遥感数据中的所有波段图像都进行处理，首先就某个波段的图像进行处理，获得该波段图像灰度值的最大值和最小值，确定新的灰度值范围，其范围一定要大于图像初始的灰度范围，实际上是一个拉伸的处理过程，然后，对所有波段的图像进行相同的处理，最后获得对比度增强后的高光谱图像数据。

进一步，所述的步骤3)对增强处理后的高光谱图像实施有序最小值滤波，以一个大小为K×L的滑动窗作为模板，在某个波段图像上进行移到，把滑动窗内的像素灰度值进行排序，选择最小的值作为模板中心像素的值，直到所有波段进行处理完毕，窗口大小一般根据条带宽带进行取值，如果条带噪声比较宽，K和L的取值相对大一点；如果条带较窄，K和L的值取相对小一些，通常情况下取7或9比较合适。

进一步，所述的步骤4)及步骤5)确定条带噪声的方向，并进行转置处理,，假设条带噪声是垂直条带，因此对于横向的条带噪声，先进行一个转置处理，变成垂直条带，进行转置处理之前，判断是否是横向还是纵向条带噪声。

进一步，所述的步骤6)对初步处理后的高光谱图像进行小波多尺度分解，再对各个分解尺度上的分解系数进行处理，达到消除条带噪声的目的，用小波对图像进行多尺度分解，图像被小波进行多尺度分解后，接着是对各个尺度上获得的分解系数进行处理，由于假设的条带噪声是垂直条带，因此除了最高的分解尺度外，其它每个分解层只要处理垂直高频系数子图像即可，判断垂直系数子图像中条带噪声的包含量，如果几乎是条带噪声，则直接把该系数值赋值为零，否则对该系数子图像继续进行小波分解，处理方式如前，直到某分解尺度上低频或垂直高频子图像几乎是条带噪声为止，处理完后，逐层进行小波逆变换，结果是经小波变换消条处理后的图像。

进一步，所述的步骤11)所述的进行矩匹配处理是对处理后的图像进行第二次消条处理。

实施例

为验证本文所提算法的有效性和可行性，进行了一系列实验，给出的一组实验数据如1所示。由天宫一号在2014年10月获取的高光谱图像数据，成像区域为陕西省某地，空间分辨率为20m。图2(a)-(d)分别表示第6、23、27和52波段图像，这四幅图像中的条带噪声具有代表性，即有明的又有暗的，而且宽度都超过一个像素。进行比较滤波的方法有巴特沃斯低通滤波(Butterworth Low-Pass Filter,BLPF)算法、矩匹配滤波(Moment Matching Filter,MMF)算法、小波变换归零滤波(Wavelet Transform Zero Filter,WTZF)算法以及OWM算法。

实验结果如图3、图4、图5和图6所示，其中(A)是BLPF算法处理的结果，(B)是MMF算法获得的结果，(C)是WTZF算法获得的结果，(D)是OWM算法获得的实验结果。

图2是图1(a)处理后的结果图。从图1(a)中可知，该图像中主要包含了一条较亮的条带噪声，而且相对而言比较窄。BLPF方法去条带效果不明显，虽然降低频率值可以加强条带噪声的消除，但是会造成图像显著模糊。MMF算法和WTZF算法均能对条带噪声有所消除，但非常不彻底。由于在小波变换域中，在去除条带噪声的同时也会丢掉一些细节信息，因此图像会有所模糊。图2表明OWM算法能消除较亮的条带噪声。

图1(b)中既包含较亮条带噪声，又包含较暗条带噪声，并且较暗的条带噪声的宽度比较宽，其实验结果如图3所示.BLPF算法和WTZF算法消除条带噪声的效果不太理想，如图3(A)和图3(C)所示。图3(B)和图3(D)表明MMF算法和OWM算法对条带噪声的消除有比较好的效果，但还不是十分满意。图1(b)是第27波段的图像，里面包含了一条较亮和一条较暗的两条条带噪声，而且条带噪声的灰度值与地物的辐射值差不多，即它们之间的差别不大。对于这类条带噪声，除了BLPF算法之外，其它三类算法均能在一定程度降低条带噪声的影响，结果如图4所示。

图1(d)中包含有多条条带噪声，而且整个图像地物的辐射值都比较低，条带噪声也主要属于暗条带。实验结果如图6所示，较满意的结果仍然是由MMF算法和OWM算法所获得，如图6(B)和图5(D)所示。BLPF算法和WTZF算法获得的结果不令人满意，因为条带噪声还是比较明显存在。