基于CWT的黑土有机质含量野外高光谱反演模型

林鹏达, 佟志军, 张继权, 赵云升, 李向前, 朱晓萌

(1.东北师范大学 环境学院, 长春130117; 2.东北师范大学 自然灾害研究所,长春130117; 3.东北师范大学 地理科学学院, 长春130024)

黑土是我国肥力最高的土壤之一，由于长期不合理的开发利用，黑土地土壤侵蚀及随之产生的土壤退化不仅危害农业生产，还带来了严重的生态环境问题。而土壤有机质是衡量土壤肥力的重要指标，也是土壤退化状态的重要指标。传统的土壤有机质含量的测定主要通过土壤野外取样化验分析来实现，精度虽高，但时效差、成本高，不能有效监测土壤的退化状况，快速准确测定土壤有机质含量已成为精准农业发展及土壤生态安全评价的必需。高光谱技术反演土壤有机质含量因其快速无损的优势得到了广泛的关注与应用，确定土壤有机质的光谱响应波段是研究的关键[1-2]。已有的研究多为室内特定条件下的土壤高光谱测量，在测试过程中存在光源功率、光源距离、照射角度等测试条件的不同，测试土样多为处理后的土壤，这在一定程度上增加了光谱数据的不确定性，也导致建立的土壤理化性质反演模型共享研究结果的困难[3-15]。重要的是在未考虑野外实际环境因素影响的情况下，室内光谱建立的反演模型可能无法直接应用于野外遥感反演。为此，基于野外光谱测量的土壤理化性质反演模型的构建十分必要。但野外条件下土壤光谱的测量会受到很多不确定的环境因素影响并引起光谱噪声，对此除了在野外光谱测量过程中规范操作，还需对采集到的高光谱数据进行相应的处理，前人研究多采用倒数、对数、微分、吸收峰深度等变换来消除土壤光谱噪声。目前，小波变换以其在信号去噪和数据压缩方面的独特优势在高光谱数据中得到应用，主要通过多尺度分解得到光谱信号中的近似特征和细节特征。小波变换中的离散小波变换(DWT)在分析高光谱数据时对输出参数解析存在困难，连续小波变换(CWT)在土壤有机质含量高光谱反演上的应用虽较少，但已有的研究中CWT处理效果较好[16-19]。因此，本研究对吉林省黑土进行采样、化验分析，在土壤野外光谱数据及其常规数学变换形式(原始光谱数据倒数的对数、一阶微分)的基础上，将连续小波变换应用于野外光谱去噪处理中，主要解决野外复杂条件下高光谱技术在土壤有机质反演方面存在的困难，为黑土有机质定量遥感反演提供支撑。

1　材料与方法

1.1　样本采集与制备

试验区选在公主岭市(43°11′—44°09′N，124°02′—125°18′E)、德惠市(44°02′—44°53′N，125°14′—126°24′E)，两市均位于吉林省黑土分布核心带上。土样采集于2016年10月1日和2日，采0—20 cm表层土，取土1 kg左右，在两市共采集35个土样(图1)，公主岭市18个土样，德惠市17个土样，采集的土壤属黑土土类中的黑土亚类，黄土质黑土土属，大部分呈中性或微酸性。在土样采集时，选择地势较为平坦、土壤裸露的耕地作为样区，土样采集前先进行土壤光谱数据测量。室内将采样点与东北典型黑土土壤分布数据进行叠加，将位于典型黑土带范围外的采样点去除，最终得到32个黑土土样。将所选定的样品自然风干、敲碎、去除杂物、磨细、过2 mm孔筛后，送由东北地理与农业生态研究所进行土壤有机质含量的测定。

图1土壤样点分布

1.2　光谱测量

野外地面土壤光谱数据利用美国Spectral Evolution便携式光谱仪测定，波长范围为350～2 500，350～1 000 nm光谱采样间隔为1.5 nm，光谱分辨率为3.5 nm；1 000～1 500 nm范围内光谱采样间隔为3.8 nm，光谱分辨率为10 nm；1 500～2 500 nm范围内光谱采样间隔为2.5 nm，光谱分辨率为7 nm。野外光谱测量时天气晴朗、无云、无风，光谱测定之前先用白板进行定标，采用8°视场角探头，距地物15 cm处垂直角度进行光谱采集；为了减少太阳高度对反射率的影响，测量时段范围为10：00—14：00。测试过程中，为消除随机噪声，每个土壤采样点随机选择临近的3个测试点，每个测试点光谱重复采集10次，把10组数据的平均值作为该测试点实际野外光谱反射率，3个测试点土壤光谱反射率算术平均后得到该样品的最终光谱反射率，并通过GPS进行定位。

1.3　数据处理

由于光谱仪器波谱两端的边缘波段存在较大的仪器噪声，因此，剔除350～399，2 451～2 500 nm波段范围的光谱数据；由于野外复杂环境条件，光谱存在某些突出的噪声波段，也进行剔除处理。每份土样仅选取400～1 880，1 940～2 450 nm的光谱数据用于后续处理分析。为寻找对土壤有机质敏感的光谱指标，除了对土壤样品的光谱反射率(R)分析外，还对光谱曲线作倒数的对数[lg(1/R)]、一阶微分(R′)的变换。

1.4　连续小波变换(CWT)

连续小波变换(continuous wavelet transformation，CWT)通过小波基函数将土壤野外高光谱数据分解成不同尺度上的一系列的小波系数[9,16-21]，公式如下：

式中：f(t)为土壤野外高光谱数据；t为光谱波段；ψa,b为小波基函数；a为尺度因子；b为平移因子。小波系数Wf(a,b)包含i和j两维，分别是分解尺度(i=1，2，…，m)和波段(j=1，2，…，n)组成的m×n矩阵，即CWT将一维光谱反射率数据转换为二维小波系数。

1.5　敏感波段的筛选

将土壤样本的有机质含量分别与原始反射率(R)、原始反射率倒数的对数[lg(1/R)]、原始反射率的一阶微分(R′)以及它们的CWT变换进行逐波段的相关分析，选取不同变换形式下显著性强相关的波段作为敏感波段；当相关性强的波段较为集中时，取其中最显著的波段为敏感波段，当相关性强的波段较为零散时，结合前人研究得到的黑土土壤有机质响应波段范围进行控制[13,22-24]。

1.6　模型的建立

多元逐步回归(Multiple stepwise regression，MSR)是对全部自变量逐一进行对因变量贡献的大小的比较，选择偏回归平方和显著的变量进入回归方程，当引入新的变量对因变量有显著影响时保留，否则剔除，通过这样的逐步过程得到回归方程。

偏最小二乘回归(Partial least square regression，PLSR)是一种基于因子分析的非参数回归方法，它囊括了主成分分析、相关分析和线性回归分析的优点，解决了自变量之间多重相关的问题，适合于自变量较多的光谱分析。

利用黑土土壤样本，将选取的敏感波段对应的原始光谱数据、不同变换形式数据以及经过CWT处理后的小波系数作为自变量，土壤有机质含量作为因变量，运用MATLAB软件建立土壤有机质含量高光谱多元逐步回归(MSR)及偏最小二乘回归(PLSR)反演模型。

1.7　模型的检验

模型稳定性采用决定系数R2来检验，R2越大，模型越稳定，可靠性越高；估算值与实测值比较采用均方根误差(root mean square error，RMSE)及1∶1线来进行评价，RMSE反映样点数据的极值效应与估值的灵敏度，RSME的值越小，表示预测结果越准确；1∶1线则在模型检验指标差异很小的情况下，通过实测值、估算值所构成的点偏离1∶1线的程度来评估模型的精度。

2　结果与分析

2.1　土壤有机质含量统计特征

选定32个样品(表1)随机分为两组，一组(21个)用于建模，另一组(11个)用于模型验证。32个土壤样品有机质含量最大为3.77%，最小为1.35%，平均为2.39%，标准差为0.52%，变异系数为0.22。

表1　土壤有机质含量描述性统计量　%

2.2　不同有机质含量的野外反射光谱特征

在野外条件下，对土壤光谱的影响因素较多，本文已将明显受干扰波段剔除(图2)。随着土壤有机质含量的增加，光谱反射率逐渐减小；光谱曲线在1 400，1 900，2 200 nm附近存在明显的水分吸收谷；土壤有机质对可见光、近红外波段的反射率具有较大的影响。土壤有机质在可见光波段存在较明显的吸收波段，主要是黑土颜色更深，并且含有较高有机质的原因；在近红外波段，光谱则主要受到NH，CH以及CO等基团的影响[25]。

2.3　基于常用光谱分析方法的相关性分析

将土壤有机质含量与光谱反射率(R)、光谱反射率倒数的对数[lg(1/R)]、光谱反射率一阶微分(R′)做相关性分析(图3)，土壤有机质含量与光谱反射率R2>0.3的波段主要有730～947，964～1 037 nm，取波段区间内777，799，826，913，1 008，1 016 nm作为反演模型的输入变量；土壤有机质含量与光谱反射率倒数的对数没有R2>0.3的波段，选取R2>0.15的波段区间，分别是546～1 410，1 440～1 820，2 150～2 358，1 845～1 880，1 940～1 960，2 393～2 430 nm间的单个波段，取799，843，879，970，1 008，1 102，1 207 nm作为反演模型的输入变量；土壤有机质含量与光谱反射率一阶微分R2>0.3的波段主要有440～637，675～790，1 624～1 664，800～900，1 700～1 760 nm间的单个波段，取波段区间内582，630，744，916，1 650，1 703 nm作为反演模型的输入变量。

图2不同有机质含量的土壤野外反射光谱曲线

2.4　基于连续小波变换的相关性分析

土壤高光谱曲线的吸收特征与Gaussian函数近似，因此，本研究选取gaus 4函数作为小波基函数，利用Matlab R2013对土样野外高光谱数据及其倒数的对数、一阶微分进行CWT变换，为了减少数据的冗余，将CWT的分解尺度设为21，22，23，…，29，即1—9个尺度。将变换后的9个尺度的小波系数与土壤有机质含量进行相关性分析，得到决定系数R2(图3—5)，其中浅色区域代表相关性强的区域。

经过R-CWT处理后，R2高于光谱反射率(R)的R2，最大提升约为0.3(图4)。利用R-CWT处理可以有效提取出光谱信息中的微妙信号，放大光谱吸收、反射特征中的局部频谱特征，能更好地捕捉土壤有机质中的特征信息。相关系数较高的波段区主要集中在410～580，610～760，830～1 280，1 750～1 870 nm；对比各分解尺度下的最大R2，发现分解尺度1，2明显较其他分解尺度小，因此，敏感波段的选取集中于3～9分解尺度，这7个分解尺度中R2在460，489，614，746，849，868，1 767 nm达到最大，其中，决定系数最高为0.6左右，这主要是由于黑土土壤有机质含量较高，黑土野外光谱曲线吸收特征较明显。

图3土壤有机质含量与R，lg(1/R)、R′的相关性分析

图4R小波系数与土壤有机质含量的相关系数

经过lg(1/R)-CWT处理后，R2同样得到了提高，最大提升约为0.3(图5)。利用lg(1/R)-CWT处理后，R2>0.3的波段区主要集中在400～430，460～500，710～760，810～1 300，1 760～1 840，2 030～2 200，2 400～2 450 nm；将各分解尺度下的R2进行降序排列并作对比，最终选取400，486，712，821，1 770，2 416 nm作为敏感波段，对应分解尺度为7，4，6，8，3，9，其中，决定系数R2最大为0.66。可见，经过CWT处理后能在一定程度上更明显地突出土壤有机质中隐藏的微弱信息。

图5lg(1/R)小波系数与土壤有机质含量的相关系数

经过R′-CWT处理后，R2变化不大且高值区域更窄(图6)。分析得出：R2>0.3的波段区主要集中在410～585，620～660，760～800，915～935，1 010～1 160，1 700～1 800 nm，根据图6的显示结合排序对比，剔除1，5，8，9这4个分解尺度数据；选取R2达到最大值的波段，依次是572，648，792，925，1 707 nm，对应分解尺度为6，7，3，2，4，其中，决定系数R2最大达到0.67。

图6R′小波系数与土壤有机质含量的相关系数

2.5　土壤有机质高光谱反演模型的优选

根据上述的分析，将R，lg(1/R)、R′得到的敏感波段对应的光谱数据以及它们对应的CWT变换后得到的小波系数作为自变量，土壤有机质含量为因变量，采用MSR，PLSR两种方法构建反演模型(表2)。由表2分析可知，R-PLSR，lg(1/R)-PLSR及R-MSR，lg(1/R)-MSR模型的建模集R2均较小，而R′-MSR 与R′-PLSR模型的建模精度均显著高于R及lg(1/R)的两种模型精度，R2最高提升了0.47，RMSE最高降低了0.18。经过CWT变换后建立的6个模型的建模集精度相较更高(R2≥0.75)，比常用光谱分析方法R2最大提升了0.59，RMSE最高降低了0.24。且在建模集中，PLSR模型明显优于对应的MSR模型。验证集的结果显示：R及lg(1/R)的两种模型的验证集R2均有所提高，而R′以及经CWT处理后的验证集R2却低于建模集，RMSE均明显增大，可能是受样本数量影响。与建模集结果相同的是：CWT变换比常规变换处理后的模型精度更高，且PLSR模型精度高于对应的MSR模型。总体上，利用CWT进行土壤高光谱数据预处理并结合PLSR的模型精度更高，能更好地筛选出对有机质敏感的信息。

表2　土壤有机质含量反演模型的建模集和验证集结果

将所有模型的估算值与实测值进行1∶1线分析(图7)。由图更清楚地可知，采用CWT对土壤高光谱数据及不同数学变换进行处理，经相关分析得到各分解层中R2最大的波段，并将其作为敏感波段，由敏感波段对应的小波系数建立的模型，其实测值和估算值样点大致分布在1∶1线附近，模型估算精度较高；而没有经过CWT处理过建立的模型效果稍差，且R，lg(1/R)、R′对应的两个模型1∶1线图中明显存在个别局外点，可能是该点在野外光谱测量时受到的外界干扰较大。在R，lg(1/R)、R′经CWT处理后建立的两个模型1∶1线图中，局外点更贴近1∶1线分布，表明选取与土壤高光谱吸收特征相似的小波基函数的连续小波变换比常规数学变换所得的结果效果更好。同时，也表明CWT处理可以在一定程度上降低外界环境条件对光谱反射曲线造成的干扰，提高光谱信噪比，进一步提高土壤有机质含量与光谱数据的相关性，精确敏感波段的选取，最终提高土壤有机质含量的反演精度。结合表2结果，经过CWT处理的R′数据更佳，PLSR因其同时考虑光谱矩阵和土壤有机质矩阵之间的相互关系较MSR更为适合。

图7　土壤有机质含量的实测值与估算值比较

3　结 论

本研究利用连续小波变换(CWT)方法对野外光谱反射率(R)、光谱反射率倒数的对数[lg(1/R)]、光谱反射率一阶微分(R′)进行不同尺度上连续小波的分解，将分解后的小波系数以及R，lg(1/R)，R′对应的敏感波段值与吉林省中部黑土土壤有机质质量分数进行了高光谱反演模型的构建，并对比了多元逐步回归(MSR)及偏最小二乘回归(PLSR)两种建模方法的反演精度。结果表明，连续小波变换作为一种信号数据处理方法，对野外土壤光谱数据噪声的去除效果更好，可以很好地检测出有效光谱信号，能够挖掘出光谱中的微弱的隐藏信息，对土壤有机质的野外高光谱反演有重要作用。R，lg(1/R)，R′间相比，R′建模集精度高于R，lg(1/R)，验证了导数光谱在一定程度上可提高土壤有机质和光谱反射率的相关性的结论，且PLSR建模精度比MSR更高；相比R，lg(1/R)，R′，它们经CWT处理后得到的小波系数与土壤有机质含量的R2均有较明显的提高，R-CWT，lg(1/R)-CWT的提高值在0.3左右；经CWT处理后的模型精度也更高，R-CWT，lg(1/R)-CWT，R′-CWT比对应的R，lg(1/R)，R′建模精度更高；经CWT处理后的两种模型的建模集中PLSR比MSR模型R2更高，验证集中PLSR与MSR模型整体精度不高，这可能是由于验证集样本数较少；由实测值样点和估算值样点在1∶1线的分布情况可得R′-CWT更集中，模型反演精度更高；总体而言，CWT处理后的模型精度更高，PLSR建模比MSR更适合土壤有机质高光谱的反演。今后，还需在吉林省黑土区增加土壤取样点以及野外光谱数据的测定，以提高模型反演精度。

参考文献：

[1]栾福明,张小雷,熊黑钢,等.基于不同模型的土壤有机质含量高光谱反演比较分析[J].光谱学与光谱分析,2013,33(1):196-200.

[2]纪文君,史舟,周清,等.几种不同类型土壤的VIS-NIR光谱特性及有机质响应波段[J].红外与毫米波学报,2012,31(3):277-282.

[3]Shi T, Cui L, Wang J, et al. Comparison of multivariate methods for estimating soil total nitrogen with visible/near-infrared spectroscopy[J]. Plant and Soil,2013,366(1/2):363-375.

[4]Gomez C, Le Bissonnais Y, Annabi M, et al. Laboratory Vis-NIR spectroscopy as an alternative method for estimating the soil aggregate stability indexes of Mediterranean soils[J]. Geoderma,2013,209:86-97.

[5]Nocita M, Kooistra L, Bachmann M, et al. Predictions of soil surface and topsoil organic carbon content through the use of laboratory and field spectroscopy in the Albany Thicket Biome of Eastern Cape Province of South Africa[J]. Geoderma, 2011,167:295-302.

[6]Xu L, Wang Q. Retrieval of soil water content in saline soils from emitted thermal infrared spectra using partial linear squares regression[J]. Remote Sensing, 2015,7(11):14646-14662.

[7]Xu C, Zeng W, Huang J, et al. Prediction of soil moisture content and soil salt concentration from hyperspectral laboratory and field data[J]. Remote Sensing,2016,8(1):42.

[8]Bao N, Wu L, Ye B, et al. Assessing soil organic matter of reclaimed soil from a large surface coal mine using a field spectroradiometer in laboratory[J]. Geoderma, 2017,288:47-55.

[9]Vohland M, Ludwig M, Harbich M, et al. Using variable selection and wavelets to exploit the full potential of visible-near infrared spectra for predicting soil properties[J]. Journal of near Infrared Spectroscopy, 2016,24(3):255-269.

[10]Gmur S, Vogt D, Zabowski D, et al. Hyperspectral analysis of soil nitrogen, carbon, carbonate, and organic matter using regression trees[J]. Sensors, 2012,12(8):10639-10658.

[11]薛利红,周鼎浩,李颖,等.不同利用方式下土壤有机质和全磷的可见近红外高光谱反演[J].土壤学报,2014,51(5):993-1002.

[12]武红旗,范燕敏,何晶,等.不同粒径土壤的反射光谱对荒漠土壤有机质含量的响应[J].草地学报,2014,22(2):266-270.

[13]刘焕军,吴炳方,赵春江,等.光谱分辨率对黑土有机质预测模型的影响[J].光谱学与光谱分析,2012,32(3):739-742.

[14]栾福明,熊黑钢,王芳,等.基于小波分析的土壤速效K含量高光谱反演[J].干旱区地理,2015(2):320-326.

[15]陈红艳,赵庚星,李希灿,等.基于小波变换的土壤有机质含量高光谱估测[J].应用生态学报,2011,22(11):2935-2942.

[16]廖钦洪,顾晓鹤,李存军,等.基于连续小波变换的潮土有机质含量高光谱估算[J].农业工程学报,2013,28(23):132-139.

[17]于雷,洪永胜,周勇,等.连续小波变换高光谱数据的土壤有机质含量反演模型构建[J].光谱学与光谱分析,2016,36(5):1428-1433.

[18]Wang G, Fang Q, Teng Y, et al. Determination of the factors governing soil erodibility using hyperspectral visible and near-infrared reflectance spectroscopy[J]. International Journal of Applied Earth Observation and Geoinformation, 2016,53:48-63.

[19]蒋金豹,何汝艳.基于连续小波变换的地下储存CO2泄漏高光谱遥感监测[J].煤炭学报,2015,40(9):2152-2158.

[20]Ullah S, Skidmore A K, Naeem M, et al. An accurate retrieval of leaf water content from mid to thermal infrared spectra using continuous wavelet analysis[J]. Science of the Total Environment, 2012,437:145-152.

[21]Cheng T, Rivard B, Sánchez-Azofeifa A G, et al. Predicting leaf gravimetric water content from foliar reflectance across a range of plant species using continuous wavelet analysis[J]. Journal of Plant Physiology, 2012,169(12):1134-1142.

[22]刘焕军,赵春江,王纪华,等.黑土典型区土壤有机质遥感反演[J].农业工程学报,2011,27(8):211-215.

[23]吴嵩.典型黑土区土壤有机质含量反演研究[D].长春：吉林大学,2016.

[24]武彦清,张柏,宋开山,等.松嫩平原土壤有机质含量高光谱反演研究[J].中国科学院大学学报,2011,28(2):187-194.

[25]Stenberg B, Rossel R A V, Mouazen A M, et al. Chapter five-visible and near infrared spectroscopy in soil science[J]. Advances in Agronomy, 2010,107:163-215.