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

摘要

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

主权项

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\lambda }}^2}|}_j\approx \frac{s\left({\lambda }_i\right)-2s\left({\lambda }_j\right)+s\left({\lambda }_k\right)}{{\left(\mathrm{\Delta \lambda }\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])/bw2 (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}}{\Sigma }}\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}}{\Sigma }}\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}{\Sigma }}{C}_i{y}_{j+i}}{N}----\left(10\right)$ 式中，y是原始光谱值，Y是平滑后光谱值；Ci是平滑窗口中第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}{\Sigma }}\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}{\Sigma }}\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中重构的高光谱图像文件。