不同降噪方式下基于高分五号影像的土壤有机质反演

刘焕军，鲍依临，孟祥添，崔 杨，张艾明，刘云超，王丹丹

（1. 赤峰学院资源环境与建筑工程学院，赤峰 024000；2. 东北农业大学公共管理与法学院，哈尔滨 150030；3. 中国科学院东北地理与农业生态研究所，长春 130012）

0 引 言

土壤有机质（Soil Organic Matter，SOM）含量是土壤肥力和土壤安全的重要指标，实现快速、精确的SOM空间制图可以为精准农业和区域土地可持续发展提供监测依据[1-2]。高光谱卫星系统被设计成通过电磁光谱的可见光-近红外区域，将表面光学性质分离成数百个光谱分辨率小于10 nm的光谱[3]，根据其在400～2 500 nm波长范围内的反射率，可以为预测土壤的物理、化学和生物学特性提供解决方案[4-5]。然而，卫星传感器的周期性偏移，载荷元器件间的电磁干扰以及波段间相互作用产生的噪声降低了图像的质量，甚至掩盖数字图像中真正的辐射信息[6-7]；同时，相比于植被、道路、建成区等地物，裸土的反射率低，光谱更易受到噪声影响，因此，对高光谱卫星数据的校正及噪声过滤是获取精准土壤光谱反射率的前提与保证[8]。

高光谱图像具有高分辨率和光谱间高相关性的特性，基于低秩的降噪方法成为研究热点之一[9]。鉴于噪声同时存在于空间域和光谱域，孟令军[10]提出基于高光谱图像像元光谱间相关性的光谱域降噪算法，分类精度提高了近8%。然而，其研究数据为合成高光谱图像，包含Cuprite Nevada的部分数据和Indian Pines，这种数据类型并不常见，其结果在指导类似研究时具有一定的限制。黄冬梅等[11]提出基于多云协同的遥感图像降噪方案，将遥感图像共享给多个云服务器，云服务器根据改进的非局部均值降噪算法在密文域中完成图像降噪。胡立栓[12]基于信息量和类间可分性的评价标准，分别提出了最优指数因子法与二进制粒子群算法，设计了基于该算法的特征选择和分类器参数同步寻优策略。证实了降噪算法在降低处理成本、提升分类精度的有效性。Meng等[13]以反射率一阶微分为基础，结合离散小波算法对高分五号（GF-5）高光谱影像进行降噪处理，使土壤有机碳的预测精度提升了约10%。现有的降噪方法可分为空间域、光谱域和空谱结合 3种。基于信号频率成分的奇异值分解方法隶属于光谱域降噪方法；离散小波变换可以分离出高频噪声信息并保留光谱中的基本成分，属于空谱结合的降噪方式[14]。中值滤波是基于排序统计理论的非线性信号处理技术，属于典型的空间域降噪方法，尤其适用于对椒盐噪声的处理[15-16]。

上述降噪方法在图像及信号处理中得到了广泛的应用，然而，大部分研究仅应用于低维（图片、音频、多光谱影像）数据类型中。GF-5高光谱卫星影像是国内光谱分辨率最高的卫星，也是国际上首次实现对大气和陆地进行综合观测的全谱段高光谱卫星，它由数以百计的波段信息组成，高维的空间信息中夹杂着较多的噪声；此外，同类像素间的光谱反射率具有较强的一致性，导致同类地物之间的光谱特征相似。由于黑龙江省西部明水县位于黑土带中心处，SOM含量空间差异明显，更具有代表性，因此本文选取明水县为研究区，利用多种降噪方式对GF-5影像进行降噪处理，评估不同方式下反射率与SOM含量的相关性，提取二维光谱指数，运用降维方法对输入量进行筛选并利用随机森林算法构建SOM反演模型，以提升 SOM预测精度，并实现高精度的SOM土壤制图，以期为实时定量监测土壤肥力变化提供依据。

1 研究区概况及数据来源

1.1 研究区概况

明水县（124°18′～125°21′E，46°44′～47°29′N）位于黑龙江省西南部，松嫩平原西北部，地处高纬度地带，平均海拔249.2 m，地势由中部向东西两侧逐渐降低，属中温带亚湿润气候，年平均气温约2.9 ℃，年平均降水量约480 mm。明水县主要包括黑土、黑钙土、风沙土，其中黑土和黑钙土的面积占总面积 60%以上。黑中含有较多腐殖质，SOM含量高；风沙土中黏粒含量较少，保水保肥能力较差，SOM含量较低。县域总面积为2 400 km2，其中耕地面积为1 087 km2，草原面积为467 km2，林地面积为360 km2。图1为明水县GF-5高光谱影像及采样点空间位置分布。

图1 研究区位置及样点布置Fig.1 Location of study fields and sampling point distribution

1.2 数据获取与处理

1.2.1 土壤数据获取

2016年5月，在研究区范围内0～20 cm深度选取采样点共38个（图1），利用GPS记录采样点的空间坐标。主要选取原则为沿公路两侧，兼顾土壤类型及物理性状选取开阔平坦的非均质性区域作为样地，样本间隔约10 km。每个样本是1个复合样本，由5～6个子样本组成，子样本在30 m×30 m大小的网格范围内随机选取，以确保SOM含量能够反映该区域平均水平。将获取的土壤样本置于布制土壤袋中存放，带回实验室进行风干处理，并去除结石，杂草根和其他杂质，研磨、干燥筛分至≤2 mm，采用重铬酸钾法测定样品的 SOM 含量。分析38个采样点 SOM 含量数据，结果表明，SOM的质量分数介于3.45%～5.53%之间，平均值为4.44%，样本标准差为4.28%，变异系数9.6。

1.2.2 高光谱卫星影像数据获取及预处理

研究区耕地范围内采用传统播种方式，作物生长期一般为一年一熟制，从 4月中旬陆续进行播种，此时土壤直接暴露于表面，地表无作物植被及秸秆残留，无积雪覆盖，这一时期被定义为“裸土窗口期”[17]。获取裸土时期（2019年 4月 24日）GF-5影像（http://data.cresda.com:90/#/home）。影像获取日期与土壤样本采集日期接近，SOM含量不会发生明显变化。GF-5影像数据的空间分辨率为30 m，地面覆盖宽度60 km，包含 330个波段，其中在可见光近红外（Visible and Near-infrared，VNIR）范围内共包含150个波段，光谱分辨率为4.28 nm；在短波红外（Shortwave Infrared，SWIR）范围内共包含180个波段，光谱分辨率为8.42 nm。由于传感器在400～430、900～1 050、2 450～2 500 nm波谱范围内信噪比较低，且在1 350～1 451、1 771～1 982 nm范围内受到大气中水蒸气吸收的影响，造成光谱数据不连续。因此，本文选取430～900、1 050～1 350、1 451～1 771、1 982～2 450 nm作为本次研究的光谱区间。使用ENVI 5.1辐射校准工具和GF-5数据头文件中的增益和偏差比对影像数据进行辐射校准，利用FLAASH大气校正模型对遥感影像进行大气校正。

2 研究方法

2.1 降噪方法

2.1.1 奇异值分解

奇异值分解（Singular Value Decomposition，SVD）可以被看作一个能够反映矩阵信息的矩阵代表值。利用奇异值分解的方法确定光谱反射信号的Hankle矩阵，对矩阵进行奇异值分解，构造逼近矩阵对含噪干涉信号进行降噪处理，奇异值越大时，它代表的信息越多，这一原理与主成分分析相似。假设Hm是m×z阶矩阵，其元素全部属于域Hm，则存在一个分解使得Hm=U∑V＊，式中U是m×m阶酉矩阵；V＊即V的共轭转置，是z×z阶酉矩阵；Σ是半正定m×z阶对角矩阵。当信号中存在噪声时，Hm为满秩矩阵，使该分解的奇异值之间存在如下关系：

式（1）中σk是SVD获得的第k个奇异值，为噪声信号中的分界点。因此，k个奇异值包含信号能量，随后的奇异值对应于信号中的噪声分量[18]，奇异值越大所包含的信息越多，取k前面若干个最大的奇异值和对应的左右奇异向量，可以基本还原出数据本身，因此可用于数据降维与压缩。该降噪方式在Matlab2016a中实现。

2.1.2 离散小波分解

离散小波分解（Discrete Wavelet Transform，DWT）是以 2的整数次幂为底对连续小波的尺度和平移参数采样，包含了信号分解与信号重构 2个过程。依据信号的长度和小波基长度将原始信号分解为高频信号和低频信号，其中高频信号主要由噪声信息构成，低频信号就是小波系数，其后的每一步均分解上一层的低频信号，迭代进行直至分解的最大尺度。在每个分解层j，建立 1个第j层的高频信号（Dj），因此原始高光谱信号（S）就等于最终分解层的低频信号（Aj）与所有分解层Dj之和，即

小波重构过程是将每层的低频信息进行重构，从而得到各层的特征光谱。常用于对高频噪声的过滤，且对于光谱曲线中的“毛刺”部分的平滑效果显著[13]。本研究在 Matlab2016a中测试了大量的母函数，最终选择‘db4’小波母函数。

2.1.3 中值滤波

中值滤波（Median Filtering，MF）是基于排序统计理论的一种广泛应用的平滑噪声和保持图像边缘的非线性滤波方法，其基本原理为：用像素点邻域灰度值的中值代替该像素点的灰度值，使周围的像素值接近真实值从而消除孤立的噪声点。对于给定的图像I(i,j)，中值滤波过程如式（3）所示：

式中(u,v)∈(1,2,…,h)×(1,2,…,l)，Imf(u,v)为处理后得到的图像；h和l是图像的宽度和高度，r和s代表窗口位移大小，W是方形窗口内的坐标集。该中值滤波法在取出脉冲噪声、椒盐噪声的同时能保留图像的边缘细节，在光学测量条纹图像的相位分析处理方法中有特殊作用，是经典的图像平滑噪声的方法。本研究在Matlab2016a中测试了不同窗口滤波器，最终选择 5×5方形窗的滤波器对影像进行降噪处理。

2.2 光谱指数的构建

光谱指数是不同波段反射率的线性或非线性组合形式，光谱指数的构建有助于提高遥感数据预测SOM含量的精度[19]。本研究探讨差值指数（Difference Index，DI）、比值指数（Ratio Index，RI）、归一化指数（Normalization Index，NDI）[20-21]与SOM之间的关系，计算过程如下：

式中Ro和Rp为波段o和p获取的光谱反射率。

在900～1 050 nm受到传感器自身缺陷的影响、且在1 350～1 451、1 771～1 982 nm受到大气中水蒸气吸收的影响，光谱数据不连续，因此，SOM估计的最佳光谱波段从 430～900、1 050～1 350、1 451～1 771、1 982～2 450 nm共计237个波段内选择。波段选择在Matlab2016a中实现。

2.3 主成分分析

主成分分析（Principal Component Analysis，PCA）的实质是基于数据统计特征的多维正交线性变化分析，将原有的高度相关变量的数据集转换为一组不相关变量，各变量间相互独立，前几个主成分信息可以最大限度地反映原始多个变量的信息[22]。在本研究中，保留了总光谱方法解释量＞85%的主成分信息作为本研究的输入量。这一过程是在Matlab2016a软件中通过princomp函数实现的。

2.4 RF预测模型

RF（Random Forest，RF）预测模型是一种集成模型，预测变量集在每个分裂点中都是随机限制的[23]，避免了计算中树与树之间的相关性问题，从而实现高精度的预测。RF的性能取决于以下3个参数的设置：生成树的数量（ntree）、叶片最小数量（nodesize）和每个节点处用于分割节点的预测变量数（mtry）。Ramón等[24]使用更多的树的数量（ntree）可以获得更稳定的估计变量重要性的结果，因此使用ntree =500。mtry默认值是预测因子总数的1/3，随后，使用迭代方法以最高预测精度确定最佳nodesize值。

2.5 SOM预测模型精度验证

研究采用建模集与验证集之比为 2:1的方式建立模型，建模集26个，验证集12个，校准和验证样本在研究区域内空间分散。为了检验不同降噪方法的效果，模型的性能通过精度度量决定系数（R2）、均方根误差（Root-Mean-Square-Error，RMSE）、测定值标准偏差与标准预测误差的比值（Residual Predictive Deviation，RPD）进行评估。R2越大，RMSE越小，证明预测精度越高；当 1.5＜RPD＜2时表明模型可以对样品含量进行粗略估测，当 RPD＞2时表明模型具有较好的定量预测能力[25]。

3 结果与分析

3.1 DWT最优分解尺度的确定

图2为9个DWT分解层（L1～L9）的光谱反射率曲线。由L1～L4可见，DWT降低了光谱曲线中的噪声，使曲线变得更加平滑，特别是在2 200～2 450 nm之间较为明显。整体来看，L1～L4可以较好地保留原始光谱曲线的特征。在L5～L9中，DWT消除了光谱曲线中大部分的细节信息，使得不同SOM含量的土壤除反射率高低的差异外，无其他明显差异，因此，过高的分解层数不利于光谱细节的保留以及SOM预测精度的提升。原始光谱曲线中，在2 200～2 450 nm处存在“小毛刺”形状，这是由于原始反射率噪声传递导致的，DWT可以将这些高频信号进一步去除，减弱噪声传递现象。将整个波段范围反射率与SOM做相关分析，由于波段信息丰富，表1列出了相关系数最高的3个波段，这3个波段主要集中在600～700 nm之间。

图2 不同离散小波分解（DWT）分解层数下不同SOM含量土壤反射光谱曲线Fig.2 Reflectance spectral curves of soils with different SOM contents under different Discrete Wavelet Transform (DWT) decomposition layers

表1 不同DWT分解层光谱反射率与SOM相关性最强的前3名Table 1 Top three based on correlation between spectral reflectance of different DWT decomposition layers and SOM

由表1可知，L1～L3均可提升光谱曲线与SOM含量的相关性，且当分解尺度为1层时，反射率与SOM含量的相关性最高。同时，这些相关性高的波段出现的位置与原始反射率（给出值）最接近，因此，将高光谱卫星影像进行1层DWT时，可以实现较高SOM预测精度，下文将选择1层DWT进行SOM反演。

3.2 不同降噪方式下土壤光谱曲线与SOM含量相关性

土壤反射率随SOM含量的增加而降低（图3），且不同 SOM 含量的土壤反射光谱曲线形状大体一致。与OR相比，不同的降噪方法均有效降低了光谱曲线中的噪声，使曲线变得更加平滑。

图3 不同降噪方法下的土壤反射光谱曲线Fig.3 Soil reflection spectral curves under different noise reduction methods

在600～1 300 nm范围内，反射率随着波长增加而增加，光谱曲线的形状差异比较明显，其中，在 900～1 050 nm间出现1个反射峰，这个反射峰是由于2个传感器连接处的噪声较大造成的。在1 450～1 700 nm范围内，反射率随着波长增加呈现先降低，后升高的趋势，且能够较好地保留光谱细节；在1 950～2 450 nm内，反射率随着波长增加先升高后降低。对比来看，SVD及DWT均对光谱曲线的不同位置进行了平滑，有效地实现了对高光谱卫星影像的降噪；然而，通过 MF降噪下的光谱曲线可知，SOM 质量分数处于 3.5%～4.0%与＞4.0%～4.5 %的光谱曲线不符合反射率随SOM含量增加而降低的规律，此外，该降噪方式缩减了不同SOM含量间光谱曲线的差异，不利于SOM含量的预测。

降噪方法可有效地提高反射率与 SOM 之间的相关性，如SVD及DWT（图4）。在800 nm附近，不同降噪条件的相关性由高到低依次为DWT、SVD、OR、MF，与 SOM 含量的相关系数依次为 0.625、0.557、0.581、0.481；在1 000～1 750 nm内，不同降噪方法的差别进一步扩大；在2 000～2 450 nm内，SVD、OR与SOM之间的相关系数大致相等。

图4 不同降噪方法下光谱反射率与SOM含量的相关性Fig.4 Correlation between spectral reflectance and SOM content under different noise reduction methods

3.3 光谱指数的选择

图5展示了不同降噪处理方式下，所选取的光谱指数（DI，RI，NDI）与SOM之间的相关关系。将与SOM相关性最高的波段组合列于表2，整体而言，与OR相比，基于DWT降噪方式下的光谱指数的相关性有所提升，SVD方式下仅RI相关性有所提升，MF处理下的相关性均低于OR反射率的相关性，不具备SOM反演的优势。

图5 土壤有机质含量与二维光谱指数相关图Fig.5 Correlation diagram of SOM content and 2-D spectral index

表2 土壤有机质含量与二维光谱指数相关性Table 2 Correlation between SOM and 2-D spectral index

3.4 不同降噪方式下的SOM预测结果

每种降噪方式包含237波段信息及3个二维光谱指数，共计240个输入变量。利用PCA方法提取贡献率大于 85%的信息进行预测，各种降噪处理下的主成分信息见表3。由表3可知，不同降噪处理下，累计贡献率均高于95%。可见，不同降噪处理后利用PCA进行波段降维，仍能够保留约95%的原始信息。

表3 基于主成分分析法的输入量统计表Table 3 Statistics of inputs based on principal component analysis method

建模集SOM预测精度R2为0.62～0.80，RMSE为1.80%～2.35%，RPD为1.63～2.23，其中以DWT降噪方式下模型模拟精度最高（图6）。验证集SOM的预测精度R2约为0.49～0.69，RMSE为2.26%～2.65%，RPD为1.41～1.80。结合 1:1线比较 SOM 的预测值与实测值，当以OR和SVD作为输入量时，SOM预测值分布比较集中；以MF作为输入量时，无论SOM含量处于较低或较高时，其预测值与实测值相差较大，说明该方式未能对数据进行有效的降噪。而使用DWT降噪方法时，其SOM预测值拟合曲线与1:1线更为接近，模拟精度更高，误差更小，且DWT方法与其他降噪方式相比具有更好的降噪效果，能够实现较高的SOM反演精度。

3.5 基于GF-5高光谱卫星数据反演图

明水县土壤分布以黑土为主，土壤肥力较高。中部黑钙土形成于低地貌环境，受低平原地下水控制，有利于钙质成分聚集，而形成钙积层，表层土被有充分的成壤条件，使腐殖质能够保存并得以逐渐加厚。西部地区地势较低，以草甸土为主，草甸土通过土壤侵蚀等作用沉积了较多的土壤养分（图7）。整体而言，中部漫川漫岗地形，SOM含量高且相对均一；东南部侵蚀严重，影响了SOM含量的积累。

图6 基于RF预测模型的SOM预测结果Fig.6 Results of SOM inversion based on RF model

图7 基于GF-5的RF模型明水县SOM分布Fig.7 Distribution of GF-5-based SOM predicted by RF model in Mingshui County

4 讨 论

高光谱卫星影像中包含丰富的波段信息，但由于受多种因素影响，影像中包含着大量噪声信息，不利于进行SOM的反演研究。因此，本文选取常见的降噪方法对高光谱卫星数据进行处理，同时，判断降噪后的影像信息是否可以提升SOM预测的精度。本文从空间域、光谱域和空谱结合 3种降噪方法出发，分别选取不同降噪原理的方式进行SOM预测。SVD方法通过保留影像中权重大的矩阵信息，实现光谱数据的降维；DWT是基于不同的分解尺度对高频信息进行过滤[26]；MF则采用象元平均化的方式以“抹平”噪声。相比于土壤本身的反射率，影像中的噪声更偏重于高频噪声。SVD和DWT降噪方法可以消除光谱数据中的冗余信息，增大不同SOM含量之间光谱曲线的差异，有助于实现高精度土壤属性的预测。对比来看，DWT将光谱曲线进行不同尺度的分解，得到不同尺度下土壤的光谱特征，使光谱与SOM含量的相关性得以提升。相比之下，MF这种象元平均化的方式未能有效的消除土壤光谱中的噪声，且降噪后的光谱与原始反射率相比，其与SOM含量的相关性有所降低。

利用高光谱遥感进行SOM预测研究的数据降维主要通过两方面实现：输入量和模型[27-28]。偏最小二乘回归模型适合处理自变量大于样点数的数据，适用于服从于正态分布的数据集，因此不具备普适性。这项研究中，每种降噪方式下均包含大量的波段信息，致使数据冗余，不利于模型的计算及分析。因此，本文采用PCA方法，将高维信息投射到更小维度的子空间中，保留原始数据大部分的信息，以缩减计算机运行时间，提升SOM预测精度[29]。本文利用不同的降噪方法对GF-5高光谱卫星数据进行降噪处理，为降低高光谱卫星影像中的噪声提供思路，然而，研究未能完全识别并消除影像中存在的噪声，在今后的研究中，应该与室内实测数据结合，充分校正数据中潜在的误差信息，将其用作开发其他领域（例如植被和水体）的参考方法。

5 结 论

考虑到高光谱卫星数据中的噪声不利于获取真实的地物反射光谱曲线，本文基于空间域、光谱域和空谱结合3种降噪方法，筛选用于GF-5高光谱数据的最佳降噪方式，主要结论如下：

1）奇异值分解（Singular Value Decomposition，SVD），离散小波变换（Discrete Wavelet Transform，DWT）能够有效降低高光谱数据中的噪声，与未降噪处理的光谱曲线相比，这 2种方法有效地减少了曲线中的“毛刺”部分，使光谱曲线更加平滑。

2）光谱指数可以提升波段信息与土壤有机质（Soil Organic Matter，SOM）的相关性，不同降噪方法对于SOM预测的潜力由高到低依次为 DWT、SVD、原始反射率（Original Reflectance，OR）、中值滤波（Median Filtering，MF），与 SOM 含量的相关性由高到低依次为 0.625、0.557、0.581、0.481。

3）降噪后，采用组成分分析降维，基于随机森林构建SOM模型。比较各降噪方式，最佳降噪方法为DWT，其验证集R2为0.69，RMSE为2.26%，RPD为1.80。基于DWT降噪方式，采用随机森林模型预测的明水县SOM分布与实际相符，进一步证实了降噪方法应用于处理高光谱卫星数据的可行性，在未来的研究中，建议使用室内实测光谱曲线作为高光谱卫星影像方法的补充，以更好地模拟不同光谱区间的噪声特征，拓展高光谱卫星在其他领域中的应用。