CN1710445A - 航空高光谱遥感图像光谱域噪声自检测与去除方法 - Google Patents

Abstract

一种涉及航空高光谱遥感图像光谱域噪声自检测与去除方法，尤指一种针对国产高光谱遥感器的高光谱图像进行处理，可以检测并去除高光谱图像光谱域中噪声，属工程科学技术中的遥感技术应用领域。该方法是通过一种适用于光谱噪声评定的变量、判定适用于光谱噪声评定的条件及找到适用于高光谱图像光谱域滤波器的方法，实现检测区分高光谱图像光谱域噪声或细微光谱特征，对光谱域较大噪声去除；主要解决如何找到适用于高光谱图像光谱域滤波器的方法系统软件及有关硬件等技术问题。本发明的优点是：能有效的去除光谱中存在噪声，并能保留光谱原有的大部分光谱特征，是地物实测光谱和高光谱图像预处理的一种有效手段。

Description

航空高光谱遥感图像光谱域噪声自检测与去除方法

技术领域

本发明涉及一种航空高光谱遥感图像光谱域噪声自检测与去除方法，该方法针对国产高光谱遥感器的高光谱图像进行处理，可以检测并去除高光谱图像光谱域中噪声，属工程科学技术中的遥感技术应用领域。

背景技术

在遥感图像数字处理中，图像噪声的去除是一项较重要而基础的工作。传统的遥感图像去噪与滤波总是从空间域对遥感图像进行噪声检测与滤波以消除噪声。高光谱图像具有空间图像和地物光谱两个方面的信息，因而噪声对高光谱图像的影响最终也会表现为空间域噪声和光谱域噪声两方面。传统的多光谱遥感空间噪声去除方法把高光谱每个波段图像看作是孤立的单幅图像，往往在空间域滤波后，光谱域中噪声并没有得到去除，有时反而会增强光谱域中噪声。反射率转换是高光谱图像预处理的重要步骤之一。

目前高光谱图像反射率转换的途径主要有三种：基于地面同步实测光谱的经验线形法、基于大气辐射传输模型的转换方法和同时使用大气辐射传输模型和地面实测光谱的混合方法，我国目前使用的一般是基于地面同步实测光谱的经验线性法。然而由于高光谱图像成像时自然光照明条件影响、地面地形影响、混合像元问题等原因引入各种噪声，表现为转换后的反射率曲线呈锯齿状。这种噪声存在于光谱域中，是需要检测并予以去除的。

Boardman，J W.注意到ATREM得到的光谱表现出一种波段间的系统误差。这种误差表现为锯齿状噪声，使得光谱看起来很不自然。他认为这些连贯误差是增益误差、仪器校正误差和大气与太阳模型影响的累积。在此基础上他提出了所谓经验平面场优化反射率变换EFFORT(the Empirical Flat Field Optimal ReflectanceTransformation)方法。经验平面场优化反射率变换EFFORT方法通过高光谱图像本身的统计找到一种轻微的改正(这种改正的增益系数接近1，偏移系数接近0)，这种改正可以使误差减小进而提高反射率的精度，使得反射率曲线有了一个视觉上的改进。EFFORT也可以使用与图像中对应地物的实测光谱参与改正，以保存一些光谱中的明显变化(比如植被红边)。

John R.Jensen也提到过这个问题，他认为反射率转换后光谱中会存在相当大的噪声，并采用一种线性校正方法称为单光谱增强SSE(Single Spectrum Enhancement)的方法(在美国高光谱“即时大气校正”ACORN软件中包含了此种方法)。该方法使用了同一地面目标(应为研究区域的主要地物类型)的转换后反射率曲线与实测光谱反射率曲线，并将转换后光谱拟合到实测地物光谱上。

常用的曲线平滑法包括均值滤波平滑法、加权平均法与最小二乘平滑法，这些方法也可以用来去除光谱域噪声。但如果使用目前常用的一维空间曲线平滑方法进行光谱噪声去除，往往会在去除噪声的同时，也将有用的细微光谱吸收等光谱特征同时去除了。而经验平面场优化反射率变换EFFORT和单光谱增强SSE方法需要一定的地面实测资料，且处理的高光谱图像区域中地表覆盖类型不能过于复杂，其应用受到一定限制。

中国于80年代中后期开始着手发展自己的高光谱图像分析系统，90年代初，中科院上海技术物理所完成了模块式航空成像光谱仪MAIS新型模块化航空成像光谱仪实用化系统；在“九五”期间，中科院上海技术物理所又相继研制成功128波段实用型模块化成像光谱仪系统(OMIS)、244波段推帚式高光谱成像仪(PHI)。这些成果标志着我国在高光谱系统领域取得了长足进展。然而由于高光谱系统本身的限制，高光谱图像在光谱域仍然会存在着较大噪声。为促进高光谱技术的进一步应用，有必要针对高光谱数据的特点，研制出一种用于光谱域噪声自检测与去除的方法，用以提高数据的质量。

发明内容

为了克服上述不足之处，本发明的主要目的旨在提供一种适用于光谱噪声评定的变量、判定适用于光谱噪声评定的条件及找到适用于高光谱图像光谱域滤波器的方法，实现检测区分高光谱图像光谱域噪声或细微光谱特征，对光谱域较大噪声去除，对光谱域的细微特征或较小噪声不予处理或进行更小程度处理的航空高光谱遥感图像光谱域噪声自检测与去除方法。

本发明要解决的技术问题是：要解决如何找到一种适用于光谱噪声评定的变量、判定适用于光谱噪声评定的条件问题；要解决如何找到适用于高光谱图像光谱域滤波器的方法系统软件及有关硬件等技术问题。

本发明解决其技术问题所采用的技术方案是：航空高光谱遥感图像光谱域噪声自检测与去除方法是通过一种适用于光谱噪声评定的变量、判定适用于光谱噪声评定的条件及找到适用于高光谱图像光谱域滤波器的方法，实现检测区分高光谱图像光谱域噪声或细微光谱特征，对光谱域较大噪声去除，对光谱域的细微特征或较小噪声不予处理或进行更小程度处理，其具体工作步骤是：

步骤1：高光谱遥感图像

高光谱遥感图像为标准格式的机载航空高光谱遥感反射率图像；

步骤2：图像光谱反射率的提取

a).高光谱遥感图像模块的输出信号传送到图像光谱反射率曲线模块；

b).打开高光谱图像，读取高光谱图像的基本信息

其中：设置波段之间宽度为bw纳米，图像存储类型为整型数int、浮点数float或其他，像元值存储字节长度nc，文件头字节数no，图像波段数nb，图像的宽度上像元数wp和高度上像元数hp；

c).计算反射率值，建立二维数组

由图像基本信息可以得到图像第i个像元的第j个波段的反射率值，从文件头偏移值ro由下式计算：

ro＝no+(wp*hp)*(j-1)*nc+i-1                   (1)

建立二维数组rb[][]，rb维数分别为wp*hp和nb，从p＝1到p＝wp*hp，b＝1到b＝nb依次读取每个像元每个波段的反射率值，并存储在rb[p][b]中；

d).二维数组rb[][]为提取出的图像光谱反射率曲线；

步骤3：图像光谱反射率曲线

图像光谱反射率曲线为步骤2中提取出的数组rb[][]，图像光谱反射率曲线模块的输出信号一路传送到一次滤波后反射率曲线模块的输入端，另一路传送到二次滤波后反射率曲线模块的输入端；

步骤4：反射率二阶导数的计算

a).图像光谱反射率曲线模块的输出信号传送到反射率曲线二阶导数模块；

b).用高光谱图像光谱的二阶导数作为光谱噪声评定的变量其中：光谱导数技术包括：

对反射光谱进行数学模拟和计算不同阶数的导数值；

确定光谱弯曲点及最大最小反射率的波长位置；

光谱导数处理强调曲线的变化和压缩均值影响；

c).计算光谱的二阶导数

光谱的二阶导数计算公式如下：

$\frac{d^{2}s}{d{\mathrm{\λ}}^2}{|}_j\≈\frac{s\left({\mathrm{\λ}}_i\right)-2s\left({\mathrm{\λ}}_j\right)+s\left({\mathrm{\λ}}_k\right)}{{\left(\mathrm{\Δ\λ}\right)}^2}$

(2)

式中：Δλ＝λ_k-λ_j＝λ_j-λ_i，λ_k＞λ_j＞λ_i；

λ_i为高光谱图像第i个波段的波长，Δλ为高光谱图像的波段之间宽度bw，s(λ_i)为像元反射率光谱曲线第i个波段的反射率值；

d).具体计算二阶导数值

对图像光谱反射率曲线数组rb[][]，对第p个像元其光谱二阶导数计算时，第b个波段的二阶导数值为：

(rb[p][b-1]-2*rb[p][b]+rb[p][b+1])/bw²        (3)

从p＝1到p＝wp*hp，b＝2到b＝nb-1依次进行光谱二阶导数计算，并将计算结果保存在数组sd[][]中，维数分别为wp*hp和nb；

在数组sd[p][b]中，当b＝1和b＝nb时，数组元素值为0；

e).数组sd[][]为反射率曲线二阶导数；

步骤5：反射率曲线二阶导数

反射率曲线二阶导数为步骤4中提取出的数组sd[][]；

步骤6：高光谱图像光谱域噪声自检测及反射率曲线噪声判定系数的计算

a).反射率曲线二阶导数模块的输出信号传送到反射率曲线噪声判定系数模块；

b).用高光谱图像光谱二阶导数的标准差作为光谱噪声的判定条件

其中：光谱的低阶导数处理对噪声影响敏感性较低，高阶微分对噪声影响敏感度高；

一阶导数反映了反射率曲线的斜率，二阶导数反映了噪声的实际分布情况；

c).用光谱二阶导数对光谱曲线进行噪声影响程度检测

1).进行光谱二阶导数计算

将二阶导数值s″(λ_i)与其平均值μ之差和给定的某一阈值作比较；

2).判定该波段噪声的大小；

3).确定该波段滤波平滑窗的大小；

4).判断不等式

对某一波段i而言，

如果不等式

|s″(λ_i)-μ|＞σ                             (4)

成立的话，那么该波段反射率值认为是有较大噪声存在，反射率曲线噪声判定系数为1，表明该波段噪声较大，式中：σ是光谱二阶导数的标准差；

如果不等式

|s″(λ_i)-μ|≤σ                             (5)

成立的话，那么该波段反射率值认为存在噪声较小，反射率曲线噪声判定系数为0；

d).具体计算二阶导数的平均值

1).对反射率曲线二阶导数数组sd[][]，第p个像元的反射率曲线二阶导数的平均值msd为：

$\mathrm{msd}=\left(\underset{b=1}{\overset{\mathrm{nb}}{\mathrm{\Σ}}}\mathrm{sd}\right[p\left]\right[b\left]\right)/\mathrm{nb}$ (6)

2).第p个像元的反射率曲线二阶导数的标准差ssd为：

$\mathrm{ssd}=\sqrt{\left[\underset{b=1}{\overset{\mathrm{nb}}{\mathrm{\Σ}}}\left(\mathrm{sd}\right[p\left]\right[b]-\mathrm{msd})\right]/\mathrm{nb}}$ (7)

3).建立二维数组de[][]

其维数分别为wp*hp和nb；

4).判断不等式

对第p个像元的第b个波段而言，

如果不等式：

|sd[p][b]-msd|＞ssd                           (8)

成立的话，de[p][b]＝1；

如果不等式：

|sd[p][b]-msd|≤ssd                           (9)

成立的话，de[p][b]＝0；从p＝1到p＝wp*hp，b＝1到b＝nb依次进行反射率曲线噪声判定系数计算，并将计算结果保存在数组de[][]中；

e).数组de[][]为反射率曲线噪声判定系数；

步骤7：反射率曲线噪声判定系数

反射率曲线噪声判定系数为步骤6中计算得到的数组de[][]；

步骤8：噪声的去除及图像光谱反射率曲线噪声的第一次滤波

a).反射率曲线噪声判定系数模块的输出信号传送到一次滤波后反射率曲线模块；

b).用赛维特斯基-高勒Savitzky-Golay平滑滤波器；

c).用简化的最小二乘拟合卷积方法对曲线进行平滑处理；

d).计算平滑后曲线各阶导数

1).简化后的最小二乘卷积方程式如下：

$Y_j=\frac{\underset{i=-m}{\overset{m}{\mathrm{\Σ}}}{C}_i{y}_{j+i}}{N}$ (10)

式中，y是原始光谱值，Y是平滑后光谱值；C_i是平滑窗口中第i个光谱值的系数，N是卷积中点值个数，j是沿原始数据纵坐标数据列的计算点下标；

2).平滑滤波法计算的卷积点为25点，并计算到光谱的第六阶导数；

3).各阶导数平滑系数的计算公式

其零阶导数二次或三次多项式拟合的任意点数平滑窗系数计算公式如下：

$p_s^{\left(0\right)}=\frac{3({3m}^2+3m-1-{5s}^2)}{(2m+3)(2m+1)(2m-1)},s=-m,...,0,...m$ (12)

式中，m为平滑窗的宽度的一半；

e).检测

在应用光谱学中光谱处理方法之前，对光谱中误差严重程度进行检测；

1).对噪声严重的地方进行高点数的平滑；

2).对噪声较小部分应用较小平滑窗进行平滑；

f).具体计算第一次滤波后图像反射率曲线

1).建立一维数组sg1[11]；

2).用来存储11点赛维特斯基-高勒滤波器系数；

3).计算sg1[]的值

按如下公式：

$\mathrm{sg}1[6+m]=\frac{3(3*5^2+3*5-1-5*m^2)}{(2*5+3)(2*5+1)(2*5-1)},m=-5,...0,...5$ (13)

计算出sg1[]的值；

4).建立二维数组rf1[][]

其维数分别为wp*hp和nb；

5)判断

如果第p个像元的第b个波段的反射率曲线噪声判定系数de[p][b]＝1，那么对第p个像元的第b个波段反射率值进行光谱域滤波，按如下公式：

$\mathrm{rf}1\left[p\right]\left[b\right]=\frac{\underset{i=-5}{\overset{5}{\mathrm{\Σ}}}\mathrm{sg}1[6+i]*\mathrm{rb}\left[p\right][b+i]}{11}$ (14)

式中5＜b＜nb-5；

如果第p个像元的第b个波段的反射率曲线噪声判定系数de[p][b]＝0，按如下公式：

rf1[p][b]＝rb[p][b]                           (15)

对图像中每个像元逐波段进行计算；

g).二维数组rf1[][]为第一次滤波后，去除了较大噪声的图像反射率曲线；

步骤9：一次滤波后反射率曲线

一次滤波后反射率曲线为步骤8中计算出的二维数组rf[][]；

步骤10：噪声的去除及图像光谱反射率曲线噪声的第二次滤波(10)

a).一次滤波后反射率曲线模块的输出信号传送到二次滤波后反射率曲线模块；

b).具体计算第二次滤波后图像反射率曲线

1).建立一维数组sg2[5]；

2).用来存储5点赛维特斯基-高勒滤波器系数；

3).计算sg2[]的值

按如下公式：

$\mathrm{sg}1[3+m]=\frac{3(3*2^2+3*2-1-2*m^2)}{(2*2+3)(2*2+1)(2*2-1)},m=-2,...0,...2$ (16)

计算出sg2[]的值；

4).建立二维数组rf2[][]

其维数分别为wp*hp和nb；

5)判断

对第p个像元的第b个波段反射率值进行光谱域滤波，

按如下公式：

$\mathrm{rf}2\left[p\right]\left[b\right]=\frac{\underset{i=-5}{\overset{5}{\mathrm{\Σ}}}\mathrm{sg}1[3+i]*\mathrm{rf}1\left[p\right][b+i]}{5}$ (17)

式中2＜b＜nb-1；

当b＝1，2，nb-1，nb时，按如下公式：

rf2[p][b]＝rb[p][b]                           (18)

对图像中每个像元逐波段进行计算；

c).二维数组rf2[][]为第二次滤波后的图像反射率曲线；

步骤11：二次滤波后反射率曲线

二次滤波后反射率曲线为步骤10中得到的二维数组rf2[][]；

步骤12：由滤波后反射率曲线重构高光谱图像

a).二次滤波后反射率曲线模块的输出信号传送到光谱域滤波后高光谱遥感图像模块；

b).建立新文件

按航空高光谱图像的格式标准，高光谱图像文件大小为(no+wp*hp*nb)*nc个字节，建立新文件；

c).写入新文件

将步骤1中读取出的高光谱图像文件头的no*nc个字节，写入新文件中；

d).写入

对p＝1，…，wp*hp，依次写入rf2[p][1]，然后再从p＝1，…，wp*hp，依次写入rf2[p][2]，直到从p＝1，…，wp*hp，依次写入rf2[p][nb]；

e).重构完成

高光谱图像文件重构完成；

步骤13：光谱域滤波后高光谱遥感图像

光谱域滤波后高光谱遥感图像为步骤12中重构的高光谱图像文件。

本发明的有益效果是：本方法能有效的去除光谱中存在噪声，并能保留光谱原有的大部分光谱特征，是地物实测光谱和高光谱图像预处理的一种有效手段；按本方法对高光谱图像进行光谱域噪声的自检测与滤波，可以在去除高光谱图像光谱域噪声的同时去除因光谱噪声引起的空间域噪声。

附图说明

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

附图1为本发明总体流程示意图；

附图2为本发明原始图像中植被反射率光谱示意图；

附图3为本发明用本方法对噪声自检测与去除后图像中植被反射率光谱示意图；

附图4为本发明含有噪声的图像反射率光谱示意图；

附图5为本发明五点赛维特斯基-高勒Savitzky-Golay滤波后光谱曲线示意图；

附图6为本发明十一点赛维特斯基-高勒Savitzky-Golay滤波后光谱曲线示意图；

附图7为用本发明方法处理后光谱曲线示意图；

附图8为本发明实测反射率光谱示意图；

附图9为本发明经验线性法得到的图像反射率光谱示意图；

附图10为本发明最小噪声变换MNF去噪后光谱示意图；

附图11为本发明经验平面场优化反射率变换EFFORT方法处理后光谱示意图；

附图12为本发明光谱均值滤波后光谱示意图；

附图13为本发明用本方法处理后光谱示意图；

附图中标号说明：

1-高光谱遥感图像；

2-图像光谱反射率的提取；

3-图像光谱反射率曲线；

4-反射率二阶导数的计算；

5-反射率曲线二阶导数；

6-噪声判定系数的计算；

7-反射率曲线噪声判定系数；

8-反射率曲线噪声的第一次滤波；

9-一次滤波后反射率曲线；

10-噪声的第二次滤波；

11-二次滤波后反射率曲线；

12-重构高光谱图像；

13-光谱域滤波后高光谱遥感图像；

具体实施方式：

请参阅附图1所示，本发明解决的最重要问题是找到方法，能检测区分高光谱图像光谱域噪声或细微光谱特征，对光谱域较大噪声去除，对光谱域的细微特征或较小噪声不予处理或进行更小程度处理。本技术分三个内容：1)找到一种适用于光谱噪声评定的变量；2)判定适用于光谱噪声评定的条件；3)找到适用于高光谱图像光谱域的滤波器。

本发明通过一种适用于光谱噪声评定的变量、判定适用于光谱噪声评定的条件及找到适用于高光谱图像光谱域滤波器的方法，实现检测区分高光谱图像光谱域噪声或细微光谱特征，对光谱域较大噪声去除，对光谱域的细微特征或较小噪声不予处理或进行更小程度处理，其具体工作步骤是：

步骤1：高光谱遥感图像(1)

高光谱遥感图像(1)为标准格式的机载航空高光谱遥感反射率图像(可参考中科院上海技术物理研究所的标准)；

步骤2：图像光谱反射率的提取(2)

高光谱遥感图像(1)模块的输出信号传送到图像光谱反射率曲线(3)模块；

打开高光谱图像，读取高光谱图像的基本信息：波段之间宽度bw(单位为纳米)，图像存储类型(整型数(int)、浮点数(float)或其他)，像元值存储字节长度nc(整型数为2个字节，浮点数为4个字节)，文件头字节数no，图像波段数nb，图像的宽度上像元数wp和高度上像元数hp。由图像基本信息可以得到图像第i个像元的第j个波段的反射率值在从文件头偏移值ro由下式计算：

ro＝no+(wp*hp)*(j-1)*nc+i-1                   (1)

建立二维数组rb[][]，rb维数分别为wp*hp和nb。从p＝1到p＝wp*hp，b＝1到b＝nb依次读取每个像元每个波段的反射率值，并存储在rb[p][b]中，二维数组rb[][]即为提取出的图像光谱反射率曲线。

步骤3：图像光谱反射率曲线(3)

图像光谱反射率曲线(3)即为步骤2中提取出的数组rb[][]，图像光谱反射率曲线(3)模块的输出信号一路传送到一次滤波后反射率曲线(9)模块的输入端，另一路传送到二次滤波后反射率曲线(11)模块的输入端。

步骤4：反射率二阶导数的计算(4)

图像光谱反射率曲线(3)模块的输出信号传送到反射率曲线二阶导数(5)模块；

本技术采用高光谱图像光谱的二阶导数作为光谱噪声评定的变量。光谱导数技术包括对反射光谱进行数学模拟和计算不同阶数的导数值以迅速地确定光谱弯曲点及最大最小反射率的波长位置。光谱导数处理强调曲线的变化和压缩均值影响。光谱导数技术对光谱信噪比非常敏感。光谱的二阶导数计算公式如下：

$\frac{d^{2}s}{d{\mathrm{\λ}}^2}{|}_j\≈\frac{s\left({\mathrm{\λ}}_i\right)-2s\left({\mathrm{\λ}}_j\right)+s\left({\mathrm{\λ}}_k\right)}{{\left(\mathrm{\Δ\λ}\right)}^2}$ (2)

其中Δλ＝λ_k-λ_j＝λ_j-λ_i，λ_k＞λ_j＞λ_i。λ_i为高光谱图像第i个波段的波长，Δλ为高光谱图像的波段之间宽度bw(见步骤2)，s(λ_i)为像元反射率光谱曲线第i个波段的反射率值。

具体计算步骤为：

对图像光谱反射率曲线数组rb[][]，对第p个像元其光谱二阶导数计算时，第b个波段的二阶导数值为：

(rb[p][b-1]-2*rb[p][b]+rb[p][b+1])/bw²        (3)

从p＝1到p＝wp*hp，b＝2到b＝nb-1依次进行光谱二阶导数计算，并将计算结果保存在数组sd[][]中，维数分别为wp*hp和nb。注意，由于公式(3)的限制，在数组sd[p][b]中，当b＝1和b＝nb时，数组元素值为0。数组sd[][]即为反射率曲线二阶导数。

步骤5：反射率曲线二阶导数(5)

反射率曲线二阶导数即为步骤4中提取出的数组sd[][]。

步骤6：高光谱图像光谱域噪声自检测-反射率曲线噪声判定系数的计算(6)

反射率曲线二阶导数(5)模块的输出信号传送到反射率曲线噪声判定系数(7)模块；

本技术采用高光谱图像光谱二阶导数的标准差作为光谱噪声的判定条件。研究表明，光谱的低阶导数处理对噪声影响敏感性较低，而高阶微分对噪声影响敏感度高。但一阶导数也同时反映了反射率曲线的斜率，二阶导数则很好的反映了噪声的实际分布情况。因此可以考虑用光谱二阶导数对光谱曲线进行噪声影响程度检测。首先进行光谱二阶导数计算，将二阶导数值s″(λ_i)的与其平均值μ之差和给定的某一阈值作比较，以判定该波段噪声的大小，从而确定该波段滤波平滑窗的大小。计算结果表明，对某一波段i而言，如果不等式

|s″(λ_i)-μ|＞σ                             (4)

成立的话，那么该波段反射率值认为是有较大噪声存在，反射率曲线噪声判定系数为1，表明该波段噪声较大，这里的σ是光谱二阶导数的标准差。如果不等式

|s″(λ_i)-μ|≤σ                             (5)

成立的话，那么该波段反射率值认为存在噪声较小，反射率曲线噪声判定系数为0。

具体计算步骤如下：

对反射率曲线二阶导数数组sd[][]，第p个像元的反射率曲线二阶导数的平均值msd为：

$\mathrm{msd}=\left(\underset{b=1}{\overset{\mathrm{nb}}{\mathrm{\Σ}}}\mathrm{sd}\right[p\left]\right[b\left]\right)/\mathrm{nb}$ (6)

第p个像元的反射率曲线二阶导数的标准差ssd为：

$\mathrm{ssd}=\sqrt{\left[\underset{b=1}{\overset{\mathrm{nb}}{\mathrm{\Σ}}}\left(\mathrm{sd}\right[p\left]\right[b]-\mathrm{msd})\right]/\mathrm{nb}}$ (7)

建立二维数组de[][]，其维数分别为wp*hp和nb。对第p个像元的第b个波段而言，如果不等式：

|sd[p][b]-msd|＞ssd                           (8)

成立的话，de[p][b]＝1。如果不等式：

|sd[p][b]-msd|≤ssd                           (9)

成立的话，de[p][b]＝0。从p＝1到p＝wp*hp，b＝1到b＝nb依次进行反射率曲线噪声判定系数计算，并将计算结果保存在数组de[][]中，数组de[][]即为反射率曲线噪声判定系数。

步骤7：反射率曲线噪声判定系数(7)

反射率曲线噪声判定系数(7)即为步骤6中计算得到的数组de[][]。

步骤8：噪声的去除-图像光谱反射率曲线噪声的第一次滤波(8)

反射率曲线噪声判定系数(7)模块的输出信号传送到一次滤波后反射率曲线(9)模块；

本技术采用赛维特斯基-高勒(Savitzky-Golay)平滑滤波器。使用简化的最小二乘拟合卷积方法对曲线进行平滑处理并可计算平滑后曲线各阶导数。简化后的最小二乘卷积一般方程式如下：

$Y_j=\frac{\underset{i=-m}{\overset{m}{\mathrm{\Σ}}}{C}_i{y}_{j+i}}{N}$ (10)

式中，y是原始光谱值，Y是平滑后光谱值。C_i是平滑窗口中第i个光谱值的系数，N是卷积中点值个数，j是沿原始数据纵坐标数据列的计算点下标。赛维特斯基-高勒(Savitzky-Golay)平滑滤波法可以计算的卷积点最高为25点，并可计算到光谱的第六阶导数。Madden改正了原赛维特斯基-高勒(Savitzky-Golay)系数的一些错误，给出了各阶导数平滑系数的计算公式。其零阶导数二次或三次多项式拟合的任意点数平滑窗系数计算公式如下：

$p_s^{\left(0\right)}=\frac{3({3m}^2+3m-1-{5s}^2)}{(2m+3)(2m+1)(2m-1)},s=-m,...,0,...m$ (12)

式中，m为平滑窗的宽度的一半。使用Madden提供的方程可以计算零到六阶导数的最小二乘拟合卷积系数值，从而计算光谱的零到六阶平滑后的导数值。由于卷积公式本身限制，光谱两端各m个点值没有进行计算。因此，得到平滑后的光谱要比原来的光谱要短。

赛维特斯基-高勒(Savitzky-Golay)平滑滤波因为可保留光谱的一些细微特征(如光谱吸收峰)，因此在光谱学中应用非常广泛，但由于高光谱遥感数据与光谱学中光谱的差别，光谱学中应用的方法不一定能直接用于高光谱遥感数据的分析中。光谱学中的数据是在严格的实验室控制条件下进行采集的，可近似认为数据的误差呈一定分布。而高光谱数据在成像时由于受自然光照明条件、地面地形起伏、大气衰减和镜头的工作状态等影响，数据的误差呈现不确定性。在应用光谱学中光谱处理方法之前需要对光谱中误差严重程度进行检测，对噪声严重的地方进行更高点数的平滑，而对噪声较小部分应用较小平滑窗进行平滑。

具体计算步骤如下：

建立一维数组sg1[11]用来存储11点赛维特斯基-高勒滤波器系数。按如下公式：

$\mathrm{sg}1[6+m]=\frac{3(3*5^2+3*5-1-5*m^2)}{(2*5+3)(2*5+1)(2*5-1)},m=-5,...0,...5$ (13)

计算出sg1[]的值。建立二维数组rf1[][]，其维数分别为wp*hp和nb。如果第p个像元的第b个波段的反射率曲线噪声判定系数de[p][b]＝1，那么对第p个像元的第b个波段反射率值进行光谱域滤波，按如下公式：

$\mathrm{rf}1\left[p\right]\left[b\right]=\frac{\underset{i=-5}{\overset{5}{\mathrm{\Σ}}}\mathrm{sg}1[6+i]*\mathrm{rb}\left[p\right][b+i]}{11}$ (14)

注意，式中5＜b＜nb-5。如果第p个像元的第b个波段的反射率曲线噪声判定系数de[p][b]＝0，按如下公式：

rf1[p][b]＝rb[p][b]                           (15)

对图像中每个像元逐波段进行计算。得到的二维数组rf1[][]即为第一次滤波后去除了较大噪声的图像反射率曲线。

步骤9：一次滤波后反射率曲线(9)

一次滤波后反射率曲线即为步骤8中计算出的二维数组rf[][]。

步骤10：噪声的去除-图像光谱反射率曲线噪声的第二次滤波(10)

一次滤波后反射率曲线(9)模块的输出信号传送到二次滤波后反射率曲线(11)模块；

具体计算步骤如下：

建立一维数组sg2[5]用来存储5点赛维特斯基-高勒滤波器系数。按如下公式：

$\mathrm{sg}1[3+m]=\frac{3(3*2^2+3*2-1-2*m^2)}{(2*2+3)(2*2+1)(2*2-1)},m=-2,...0,...2$ (16)

计算出sg2[]的值。建立二维数组rf2[][]，其维数分别为wp*hp和nb。对第p个像元的第b个波段反射率值进行光谱域滤波，按如下公式：

$\mathrm{rf}2\left[p\right]\left[b\right]=\frac{\underset{i=-5}{\overset{5}{\mathrm{\Σ}}}\mathrm{sg}1[3+i]*\mathrm{rf}1\left[p\right][b+i]}{5}$ (17)

注意，式中2＜b＜nb-1。当b＝1，2，nb-1，nb时，按如下公式：

rf2[p][b]＝rb[p][b]                           (18)

对图像中每个像元逐波段进行计算。得到的二维数组rf2[][]即为第二次滤波后的图像反射率曲线。

步骤11：二次滤波后反射率曲线(11)

二次滤波后反射率曲线(11)即为步骤10中得到的二维数组rf2[][]。

步骤12：由滤波后反射率曲线重构高光谱图像(12)

二次滤波后反射率曲线(11)模块的输出信号传送到光谱域滤波后高光谱遥感图像模块；

按航空高光谱图像的格式标准，高光谱图像文件大小为(no+wp*hp*nb)*nc个字节，建立新文件。将步骤1中读取出的高光谱图像文件头的no*nc个字节，写入新文件中。然后对p＝1，…，wp*hp，依次写入rf2[p][1]，然后再从p＝1，…，wp*hp，依次写入rf2[p][2]，直到从p＝1，…，wp*hp，依次写入rf2[p][nb]。高光谱图像文件重构完成。

步骤13：光谱域滤波后高光谱遥感图像(13)

光谱域滤波后高光谱遥感图像(13)即为步骤12中重构的高光谱图像文件。

本发明的具体实施例：

按照本发明的技术流程图各步骤，对上海2003年世博会区域PHI高光谱图像进行光谱域噪声的自检测与去除，处理效果比较好。

请参阅附图2所示，为原始图像中植被反射率光谱；

请参阅附图3所示，为用本方法对噪声自检测与去除后图像中植被反射率光谱。从图3中可以看出，本专利方法DSGF对高光谱图像最大的作用首先体现在对光谱平滑上。未经光谱域噪声去除的图像光谱曲线由于受各种噪声影响而呈现锯齿状(图2)，在400nm-500nm与750nm-980nm之间波段尤其明显。经本方法处理后，图像光谱在这两个光谱区间明显平滑许多(图3)，同时从图3还可以看到，光谱曲线仍保持了一些如吸收峰等细微的光谱特征。

请参阅附图4所示，为含有噪声的图像反射率光谱；

请参阅附图5所示，为五点赛维特斯基-高勒(Savitzky-Golay)滤波后光谱曲线图；

请参阅附图6所示，为十一点赛维特斯基-高勒(Savitzky-Golay)滤波后光谱曲线图；

请参阅附图7所示，为本发明方法处理后光谱曲线图。赛维特斯基-高勒(Savitzky-Golay)滤波的效果由参与滤波的点数决定，点数越高，平滑程度也越高，但初始光谱的特征就越难以保留。从图4中可以看出，初始反射率光谱曲线的两个矩形框区域呈现不同程度的锯齿状，其中尤其以近950nm附近波段为据。这两个区域的光谱值表明该波段区域含有噪声而且噪声的程度不一。从图5中看到，5点滤波使得噪声情况有所缓解，同时也保持了大部分的光谱细微特征，但在噪声较严重区域仍然呈锯齿状；从图6中看到，11点滤波后的光谱非常平滑，然而一些正常的光谱特征也被掩盖；从图7中看到，本发明方法处理后的光谱曲线在噪声严重的区域能较彻底的清除噪声，同时也能保持光谱原有的细微特征，这证明本发明方法是一种较好的高光谱图像光谱域噪声自检测与去除方法。

为了与高光谱其他一些噪声去除方法进行比较，我们还使用了多种方法对高光谱图像进行处理。

请参阅附图8所示，为实测反射率光谱；

请参阅附图9所示，为经验线性法得到的图像反射率光谱；

请参阅附图10所示，为最小噪声变换MNF去噪后光谱；

请参阅附图11所示，为经验平面场优化反射率变换EFFORT方法处理后光谱；

请参阅附图12所示，为光谱均值滤波后光谱；

请参阅附图13所示，为本发明方法处理后光谱。

从图9可以看出，与图8的图像光谱相比，利用经验线性法得到的图像反射率包含有较大噪声，如果在高光谱遥感应用中直接利用这种包含较大噪声的反射率值的话，无疑会使得应用结果也呈现噪声值。因此，在进一步高光谱应用处理步骤之前，对图像进行光谱域的噪声去除是有必要的。图10为最小噪声变换MNF后，对最小噪声变换MNF变换结果中波段进行噪声去除，然后进行MNF逆变换后得到的图像光谱，最小噪声变换MNF方法应用正交变化进行噪声分离，总体上说来，对噪声的去除效果还是比较不错的，但位于950nm波段附近的反射率值仍然有锯齿状噪声存在。

图11中可以看到经验平面场优化反射率变换EFFORT方法对光谱处理的效果不是很理想。图12和图13表明，光谱均值滤波和本专利方法DSGF方法滤波对光谱的处理效果相对其他方法而言要好。而光谱均值滤波后的光谱过于平滑，虽然在视觉上与实测光谱非常的相符，但由于高光谱的多种噪声影响，加上高光谱的波段宽度比地物光谱仪的波段宽度要大，因而图像光谱与实测光谱应该是有着较大的区别的，这些区别应该表现在光谱的一些细微特征上，光谱均值滤波虽然保留了光谱的整体趋势，但是光谱的细微特征不能辨别。

为了进一步比较各种处理方法处理后图像光谱与原图像光谱和实测光谱的相似程度，对上述六种光谱两两之间底光谱欧几里德距离和光谱角距离进行计算，得出了如表1和表2的结果。表1四种方法处理后的光谱与实测光谱和原图像光谱的欧几里德距离

	图十光谱	图十一光谱	图十二光谱	图十三光谱
					图八光谱	0.2136	O.2414	0.2033	0.1956
图九光谱	0.0973	0.0431	0.1258	0.1126

(接下页)

表2四种方法处理后的光谱与实测光谱和原图像光谱的光谱角距离

	图十光谱	图十光谱	图十一光谱	图十二光谱
					图八光谱	0.0835	0.0970	0.0773	0.0769
图九光谱	0.0425	0.0188	0.0549	0.0489

由欧几里德距离与光谱角距离的定义可知，距离值越小，两光谱越相似。从表1、表2中可以看出，图十三中所示用本发明方法处理后光谱与均值滤波光谱相比，图十三所示光谱与实测光谱的距离和与图像光谱的距离都要小于均值滤波光谱与实测光谱和图像光谱的距离。这说明，均值滤波后光谱过于平滑，而本发明方法处理后光谱保留了光谱的大部分特征，在进行光谱匹配时最具有优势。

Claims (1)

1、一种航空高光谱遥感图像光谱域噪声自检测与去除方法，其特征在于：该方法通过一种适用于光谱噪声评定的变量、判定适用于光谱噪声评定的条件及找到适用于高光谱图像光谱域滤波器的方法，实现检测区分高光谱图像光谱域噪声或细微光谱特征，对光谱域较大噪声去除，对光谱域的细微特征或较小噪声不予处理或进行更小程度处理，其具体工作步骤是：

步骤1：高光谱遥感图像(1)

高光谱遥感图像(1)为标准格式的机载航空高光谱遥感反射率图像；

步骤2：图像光谱反射率的提取(2)

a).高光谱遥感图像(1)模块的输出信号传送到图像光谱反射率曲线(3)模块；

b).打开高光谱图像，读取高光谱图像的基本信息

其中：设置波段之间宽度为bw纳米，图像存储类型为整型数int、浮点数float或其他，像元值存储字节长度nc，文件头字节数no，图像波段数nb，图像的宽度上像元数wp和高度上像元数bp；

c).计算反射率值，建立二维数组

由图像基本信息可以得到图像第i个像元的第j个波段的反射率值，从文件头偏移值ro由下式计算：

ro＝no+(wp*hp)*(j-1)*nc+i-1                    (1)

建立二维数组rb[][]，rb维数分别为wp*hp和nb，从p＝1到p＝wp*hp，b＝1到b＝nb依次读取每个像元每个波段的反射率值，并存储在rb[p][b]中；

d).二维数组rb[][]为提取出的图像光谱反射率曲线；

步骤3：图像光谱反射率曲线(3)

图像光谱反射率曲线(3)为步骤2中提取出的数组rb[][]，图像光谱反射率曲线(3)模块的输出信号一路传送到一次滤波后反射率曲线(9)模块的输入端，另一路传送到二次滤波后反射率曲线(11)模块的输入端；

步骤4：反射率二阶导数的计算(4)

a).图像光谱反射率曲线(3)模块的输出信号传送到反射率曲线二阶导数(5)模块；

b).用高光谱图像光谱的二阶导数作为光谱噪声评定的变量

其中：光谱导数技术包括：

对反射光谱进行数学模拟和计算不同阶数的导数值；

确定光谱弯曲点及最大最小反射率的波长位置；

光谱导数处理强调曲线的变化和压缩均值影响；

c).计算光谱的二阶导数

光谱的二阶导数计算公式如下：

${\frac{d^{2}s}{{\mathrm{d\λ}}^2}|}_j\≈\frac{s\left({\mathrm{\λ}}_i\right)-2s\left({\mathrm{\λ}}_j\right)+s\left({\mathrm{\λ}}_k\right)}{{\left(\mathrm{\Δ\λ}\right)}^2}---\left(2\right)$

式中：Δλ＝λ_k-λ_j＝λ_j-λ_i，λ_k＞λ_j＞λ_i；

λi为高光谱图像第i个波段的波长，Δλ为高光谱图像的波段之间宽度bw，s(λ_i)为像元反射率光谱曲线第i个波段的反射率值；

d).具体计算二阶导数值

对图像光谱反射率曲线数组rb[][]，对第p个像元其光谱二阶导数计算时，第b个波段的二阶导数值为：

(rb[p][b-1]-2*rb[p][b]+rb[p][b+1])/bw²            (3)

从p＝1到p＝wp*hp，b＝2到b＝nb-1依次进行光谱二阶导数计算，并将计算结果保存在数组sd[][]中，维数分别为wp*hp和nb；

在数组sd[p][b]中，当b＝1和b＝nb时，数组元素值为0；

e).数组sd[][]为反射率曲线二阶导数；

步骤5：反射率曲线二阶导数(5)

反射率曲线二阶导数(5)为步骤4中提取出的数组sd[][]；

步骤6：高光谱图像光谱域噪声自检测及反射率曲线噪声判定系数的计算(6)

a).反射率曲线二阶导数(5)模块的输出信号传送到反射率曲线噪声判定系数(7)模块；

b).用高光谱图像光谱二阶导数的标准差作为光谱噪声的判定条件

其中：光谱的低阶导数处理对噪声影响敏感性较低，高阶微分对噪声影响敏感度高；

一阶导数反映了反射率曲线的斜率，二阶导数反映了噪声的实际分布情况；

c).用光谱二阶导数对光谱曲线进行噪声影响程度检测

1).进行光谱二阶导数计算

将二阶导数值s″(λ_i)与其平均值μ之差和给定的某一阈值作比较；

2).判定该波段噪声的大小；

3).确定该波段滤波平滑窗的大小；

4).判断不等式

对某一波段i而言，

如果不等式

|s″(λ_i)-μ|＞σ                            (4)

成立的话，那么该波段反射率值认为是有较大噪声存在，反射率曲线噪声判定系数为1，表明该波段噪声较大，式中：σ是光谱二阶导数的标准差；

如果不等式

|s″(λ_i)-μ|≤σ                           (5)

成立的话，那么该波段反射率值认为存在噪声较小，反射率曲线噪声判定系数为0；

d).具体计算二阶导数的平均值

1).对反射率曲线二阶导数数组sd[][]，第p个像元的反射率曲线二阶导数的平均值msd为：

$\mathrm{msd}=\left(\underset{b=1}{\overset{\mathrm{nb}}{\mathrm{\Σ}}}\mathrm{sd}\right[p\left]\right[b\left]\right)/\mathrm{nb}---\left(6\right)$

2).第p个像元的反射率曲线二阶导数的标准差ssd为：

                 -3-

$\mathrm{ssd}=\sqrt{\left[\underset{b=1}{\overset{\mathrm{nb}}{\mathrm{\Σ}}}\left(\mathrm{sd}\right[p\left]\right[b]-\mathrm{msd})\right]/\mathrm{nb}}---\left(7\right)$

3).建立二维数组de[][]

其维数分别为wp*hp和nb；

4).判断不等式

对第p个像元的第b个波段而言，

如果不等式：

|sd[p][b]-msd|＞ssd                         (8)

成立的话，de[p][b]＝1；

如果不等式：

|sd[p][b]-msd|≤ssd                         (9)

成立的话，de[p][b]＝0；从p＝1到p＝wp*hp，b＝1到b＝nb依次进行反射率曲线噪声判定系数计算，并将计算结果保存在数组de[][]中；

e).数组de[][]为反射率曲线噪声判定系数；

步骤7：反射率曲线噪声判定系数(7)

反射率曲线噪声判定系数(7)为步骤6中计算得到的数组de[][]；

步骤8：噪声的去除及图像光谱反射率曲线噪声的第一次滤波(8)

a).反射率曲线噪声判定系数(7)模块的输出信号传送到一次滤波后反射率曲线(9)模块；

b).用赛维特斯基-高勒Savitzky-Golay平滑滤波器；

c).用简化的最小二乘拟合卷积方法对曲线进行平滑处理；

d).计算平滑后曲线各阶导数

1).简化后的最小二乘卷积方程式如下：

$Y_j=\frac{\underset{i=-m}{\overset{m}{\mathrm{\Σ}}}{C}_i{y}_{j+i}}{N}----\left(10\right)$

式中，y是原始光谱值，Y是平滑后光谱值；C_i是平滑窗口中第i个光谱值的系数，N是卷积中点值个数，j是沿原始数据纵坐标数据列的计算点下标；

2).平滑滤波法计算的卷积点为25点，并计算到光谱的第六阶导数；

3).各阶导数平滑系数的计算公式

其零阶导数二次或三次多项式拟合的任意点数平滑窗系数计算公式如下：

$p_s^{\left(0\right)}=\frac{3({3m}^2+3m-1-{5s}^2)}{(2m+3)(2m+1)(2m-1)},s=-m,...0,...m---\left(12\right)$

式中，m为平滑窗的宽度的一半；

e).检测

在应用光谱学中光谱处理方法之前，对光谱中误差严重程度进行检测；

1).对噪声严重的地方进行高点数的平滑；

2).对噪声较小部分应用较小平滑窗进行平滑；

f).具体计算第一次滤波后图像反射率曲线

1).建立一维数组sg1[11]；

2).用来存储11点赛维特斯基-高勒滤波器系数；

3).计算sg1[]的值

按如下公式：

$\mathrm{sg}1[6+m]=\frac{3(3*5^2+3*5-1-5*m^2)}{(2*5+3)(2*5+1)(2*5-1)}m=-5,...0,...5---\left(13\right)$

计算出sg1[]的值；

4).建立二维数组rff1[][]

其维数分别为wp*hp和nb；

5)判断

如果第p个像元的第b个波段的反射率曲线噪声判定系数de[p][b]＝1，那么对第p个像元的第b个波段反射率值进行光谱域滤波，按如下公式：

$\mathrm{rf}1\left[p\right]\left[b\right]=\frac{\underset{i=-5}{\overset{5}{\mathrm{\Σ}}}\mathrm{sg}1[6+i]*\mathrm{rb}\left[p\right][b+i]}{11}---\left(14\right)$

式中5＜b＜nb-5；

如果第p个像元的第b个波段的反射率曲线噪声判定系数de[p][b]＝0，按如下公式：

rf1[p][b]＝rb[p][b]                        (15)

对图像中每个像元逐波段进行计算；

g).二维数组rf1[][]为第一次滤波后，去除了较大噪声的图像反射率曲线；

步骤9：一次滤波后反射率曲线(9)

一次滤波后反射率曲线(9)为步骤8中计算出的二维数组rf[][]；

步骤10：噪声的去除及图像光谱反射率曲线噪声的第二次滤波(10)

a).一次滤波后反射率曲线(9)模块的输出信号传送到二次滤波后反射率曲线(11)模块；

b).具体计算第二次滤波后图像反射率曲线

1).建立一维数组sg2[5]；

2).用来存储5点赛维特斯基-高勒滤波器系数；

3).计算sg2[]的值

按如下公式：

$\mathrm{sg}1[3+m]=\frac{3(3*2^2+3*2-1-2*m^2)}{(2*2+3)(2*2+1)(2*2-1)}m=-2,...0,...2---\left(16\right)$

计算出sg2[]的值；

4).建立二维数组rf2[][]

其维数分别为wp*hp和nb；

5)判断

对第p个像元的第b个波段反射率值进行光谱域滤波，

按如下公式：

$\mathrm{rf}2\left[p\right]\left[b\right]=\frac{\underset{i=-5}{\overset{5}{\mathrm{\Σ}}}\mathrm{sg}1[3+i]*\mathrm{rf}1\left[p\right][b+i]}{5}---\left(17\right)$

式中2＜b＜nb-1；

当b＝1，2，nb-1，nb时，按如下公式：

rf2[p][b]＝rb[p][b]                                (18)

对图像中每个像元逐波段进行计算；

c).二维数组rf2[][]为第二次滤波后的图像反射率曲线；

步骤11：二次滤波后反射率曲线(11)

二次滤波后反射率曲线(11)为步骤10中得到的二维数组rf2[][]；

步骤12：由滤波后反射率曲线重构高光谱图像(12)

a).二次滤波后反射率曲线(11)模块的输出信号传送到光谱域滤波后高光谱遥感图像模块；

b).建立新文件

按航空高光谱图像的格式标准，高光谱图像文件大小为(no+wp*hp*nb)*nc个字节，建立新文件；

c).写入新文件

将步骤1中读取出的高光谱图像文件头的no*nc个字节，写入新文件中；

d).写入

对p＝1，...，wp*hp，依次写入rf2[p][1]，然后再从p＝1，...，wp*hp，依次写入rf2[p][2]，直到从p＝1，...，wp*hp，依次写入rf2[p][nb]；

e).重构完成

高光谱图像文件重构完成；

步骤13：光谱域滤波后高光谱遥感图像(13)

光谱域滤波后高光谱遥感图像(13)为步骤12中重构的高光谱图像文件。

Images