基于CWT-sCARS的东北旱作农田土壤有机质高光谱反演

勾宇轩 赵云泽 李 勇 卓志清 曹 梦 黄元仿

(1.中国农业大学土地科学与技术学院， 北京 100193； 2.自然资源部农用地质量与监控重点实验室， 北京 100193；3.浙江大学农业遥感与信息技术应用研究所， 杭州 310058；4.中国建筑一局(集团)有限公司生态园林分公司， 北京 100071)

0 引言

土壤有机质(Soil organic matter, SOM)作为反映土壤肥力状况和退化程度的重要指标，其含量的动态变化受到研究者的广泛关注[1-2]。传统SOM含量的测定主要通过湿烧法、干烧法、重铬酸钾容量法等化学分析法来实现，操作复杂、周期长、成本高，且为点状信息数据，不适合快速、大面积测定，难以满足土壤退化动态监测和精准防治的要求。高光谱技术凭借快速、简便、无污染、不破坏等优势，在土壤属性预测研究中得到了广泛应用。有机质含量不同的土壤在可见光和近红外光谱区域存在独特的光谱特征，为SOM含量的快速测定开辟了新的途径[3-4]。近年来，国内外学者利用光谱技术进行SOM的预测取得了较好的效果[5-8]。

然而，由于受到土壤样本表面、光谱测试环境、光谱高频随机噪声、杂散光等干扰及光谱数据的信息冗余、多重共线性、吸收峰重叠等现象影响，模型的性能和精度易受到严重影响。前人研究多采用倒数、对数、微分、吸收峰深度等传统变换来消除土壤光谱噪声，通过无信息变量消除(Uninformative variables elimination，UVE)、连续投影算法(Successive projections algorithm，SPA)、竞争性自适应重加权采样(Competitive adaptive reweighted sampling，CARS)等方法滤除无效变量或冗余变量[9-10]，优选出SOM的敏感波段，以降低模型复杂度，提升模型的预测精度和稳定性。如何对高光谱数据进行处理，消除干扰，滤除冗余、共线信息，筛选敏感波段，已成为当前研究的热点之一。

连续小波变换(Continuous wavelet transform，CWT)通过多尺度分解，得到光谱信号中的近似特征和细节信息，在高光谱数据的噪声去除和数据压缩上具有独特优势[11-13]。CARS算法在敏感波段的筛选上具有很好的性能，但只考虑了变量回归系数，并不能总是反映变量重要性的真实信息。ZHENG等[14]构建了稳定性竞争自适应重加权采样(Stability competitive adaptive reweighted sampling, sCARS)算法克服了这一缺点。刘国海等[15]提出了基于sCARS策略的特征波长筛选方法，与UVE、CARS等现有方法相比，能够有效地增强光谱特征波段选择的准确性和稳定性。何东健等[16]比较了sCARS、SPA、CARS和RF(随机蛙跳算法)4种波长选择方法对建模效果的影响，发现sCARS的结果最佳。李冠稳等[17]使用sCARS算法结合随机森林模型，实现了SOM含量的快速、无损、精准估测。

东北地区土层深厚，土壤性状良好，是我国重要的商品粮基地，但是由于长期的高强度利用，加之侵蚀作用的影响，其SOM含量显著下降，建立可靠的SOM含量动态监测方法，对东北土壤退化防治和耕地质量提升有重要的现实意义。本文探究东北旱作农田典型土壤类型的光谱反射特性，采用CWT对高光谱数据进行分析，利用sCARS算法筛选敏感波段，提取有益信息，提高SOM预测模型的精度和稳定性，为建立可靠的土壤退化监测网络提供有力的支撑，进一步推动高光谱技术在SOM含量预测中的应用。

1 材料与方法

1.1 土壤样品采集

在黑龙江、吉林、辽宁3省地形坡度小于5°、1 km2网格旱地占耕地面积40%以上的区域内，以15 km×15 km网格布点，结合土壤类型进行分层抽样，综合考虑种植体系、分布面积和集中程度等因素，共确定118个采样点(图1)，其中黑土采样点37个，黑钙土采样点33个，潮土采样点28个，棕壤采样点20个。土壤样品采集时间为2017年5—6月，采样深度为0～20 cm。采集的土壤样品自然风干、研磨、过2 mm筛后分为两份，一份用于土壤有机质含量测定，一份用于土壤光谱测试。

图1 研究区及采样点地理位置Fig.1 Study area and sampling point

1.2 土壤有机质含量测定

采用重铬酸钾-外加热法测定样品的SOM含量，统计特征见表1。东北旱作农田SOM含量均值为26.70 g/kg，整体处于较高水平，数据离散程度一般，属于中等变异强度(15%～100%)。黑土SOM含量最高，且变异系数最小，离散程度最低，潮土和棕壤SOM含量最低。

表1 土壤有机质含量的统计特征Tab.1 Statistical characteristics of SOM content

1.3 土壤光谱测试

在暗室中使用ASD FieldSpec 4型便携式地物光谱仪进行土壤光谱数据测试。将土壤样本放置于深1.5 cm、半径5 cm的玻璃盛样皿内，选择50 W卤素灯为光源，保持15°顶角，距离土壤样品表面30 cm，采用8°视场角，使探头垂直位于土壤样本正上方15 cm处，使用前先进行白板校正，仪器稳定后再进行实验。每个土样采集10条光谱曲线，算术平均后得到该土样的光谱反射率。仪器的光谱波长为350～2 500 nm，重采样间隔为1 nm，共输出2 151个波段。

1.4 光谱预处理

剔除波段350～399 nm和2 401～2 500 nm，以消除边缘噪声的干扰。通过Savitzky-Golay算法(多项式阶数为2，点数为11)平滑400～2 400 nm的土壤光谱曲线，并采用倒数对数、一阶微分、连续统去除(Continuum removal，CR)和CWT 4种形式对平滑后的光谱曲线进行处理。

CWT源于傅里叶算法，可同时在频率和时间开展分析，从信号中提取有效信息，其通过小波基函数将一维的土壤高光谱数据分解成包含波段和不同分解尺度的二维小波系数，计算式为

(1)

(2)

式中Wf(a,b)——小波系数，包含i、j两维，分别是分解尺度(i=1,2,…,m)和波段(j=1,2,…,n)组成的m×n矩阵

f(t)——土壤高光谱反射率

t——光谱波段(400～2 400 nm)

ψa,b——小波基函数

a——尺度因子

b——平移因子

1.5 敏感波段筛选

sCARS算法首先计算各波长变量的稳定性值，然后利用自适应重加权采样技术(Adaptive reweighted sampling，APS)和指数衰减函数(Exponentially decreasing function，EDF)筛选出回归系数绝对值大且稳定性高的变量，多次循环迭代此过程。最终通过十折交互检验对每次循环后所得的变量子集进行检验，优选出交互验证均方根误差(RMSECV)最小的变量子集[18]。

1.6 模型建立与验证

SOM反演模型构建采用偏最小二乘回归(Partial least square regression，PLSR)，其有利于解决自变量之间多重相关的现象，能够避免光谱建模中的过拟合问题。选取决定系数(Determination coefficients，R2)和均方根误差(Root mean squared error，RMSE)从模型的稳定性和预测能力两方面对模型的精度进行验证。

光谱数据的平滑、PLSR建模分析在Unscrambler X(CAMO Inc., 挪威)软件中实现，CR处理通过ENVI 5.3(ESRI Inc., 美国)软件实现，CWT处理和sCARS波段筛选在Matlab R2018a(The Mathworks Inc., 美国)环境下完成。

2 结果与讨论

2.1 土壤光谱特征分析

按照SOM含量将所有土壤样品分为6个等级，计算各等级所有土壤光谱反射率平均值，得到不同有机质含量下土壤的光谱反射率曲线(图2)。研究区土壤光谱反射率随着SOM含量的增加而减小，不同有机质含量等级的土壤光谱反射率曲线整体变化趋势相近。整体来看，光谱曲线表现平缓，光谱反射率随波长而增大，在400～900 nm波段区域反射率增加较快，斜率较大，而在近红外的900～2 400 nm波段反射率增加较为平缓。

图2 不同有机质含量的土壤光谱反射率曲线Fig.2 Reflectance spectra curves of soils with different SOM contents

由不同土壤类型的光谱反射率曲线变化特征可知(图3)，黑土和黑钙土的光谱反射率较低，这是由于黑土、黑钙土SOM含量均值较高，颜色较深，土壤质地粘细，其光谱曲线属于富含有机质类型[19]。一般情况下，随着SOM含量的增加，光谱反射率减小。潮土和棕壤虽然SOM含量均值相近，但与棕壤相比，潮土的光谱反射率相对较低，可能受到氧化铁含量等其他因素的影响[20]。整体来看，不同土壤类型的光谱特征吸收带出现的位置总体上一致，只是深度有所不同。在可见光波段，主要受土壤发色团和有机质本身颜色的影响，存在较宽的吸收波段；在近红外波段，则主要受到N—H、C—H和C—O等基团的影响。在波长1 400、1 900、2 200 nm处土壤光谱曲线有3个明显的水分吸收谷，这可能是黏土矿物中所含的水分子和羟基的吸收带[21]。

图3 不同类型土壤的光谱反射率曲线Fig.3 Reflectance spectra of different soils

2.2 传统光谱变换的相关性分析

3种传统光谱变换对相关性的提升效果并不理想(图4)。原始光谱与SOM含量的决定系数最大为0.41，倒数对数变换后，土壤光谱与SOM含量的相关性无明显变化；经过连续统去除处理，相关性有所降低，决定系数最大仅为0.36；一阶微分处理后相关性有一定的提升，决定系数最大可达0.46，提高了0.05，这是由于微分技术对分离重叠样品、消除平行噪声的干扰有一定的效果。原始光谱和倒数对数R2>0.4的波段主要分布在670～890 nm，并在794 nm处达到最大值；一阶微分R2>0.4的波段主要分布在530～630 nm和1 920～1 950 nm，并在594 nm和1 930 nm处达到最大值。

图4 不同光谱变换形式的相关性分析Fig.4 Correlation analysis of different spectral transformation forms

2.3 连续小波变换的相关性分析

由于土壤光谱曲线特征与Gaussian函数近似，选取Gaussian4函数作为小波基函数[22]，CWT分解尺度设定为21、22、…、210，即分解尺度1～10。将变换后的小波系数与SOM含量进行相关性分析，得到不同类型土壤小波系数与SOM含量的决定系数矩阵(图5)，图中红色区域代表相关性强的区城，蓝色区域代表相关性弱的区城。

由图5可以看出，有效光谱信息主要集中于中高尺度分解。在分解尺度为1～3时，信息相对均一，部分光谱细节特征消失，有效信息比较匮乏；随着分解尺度增大，光谱对SOM的敏感性也逐步增大，CWT处理后土壤光谱与SOM含量的相关性明显高于原始光谱，这表明土壤的部分有效光谱信息被隐藏，利用连续小波分析，适度增大分解尺度可有效挖掘其内的隐含信息[23]。

CWT处理后的潮土光谱与SOM含量的相关性最强，决定系数最高可达0.69，而黑土较差。将决定系数最高的5%区域作为SOM的敏感区域，黑土的敏感区域主要分布在563～595 nm、563～612 nm、788～987 nm、2 084～2 126 nm、425～468 nm、755～902 nm、857～947 nm、1 282～1 630 nm，对应的尺度分别为第5、6、6、7、8、8、9、9尺度；黑钙土的敏感区域主要分布在1 384～1 425 nm、2 319～2 359 nm、479～532 nm、872～912 nm、1 416～1 469 nm、2 166～2 207 nm、2 257～2 301 nm、1 498～1 703 nm、1 125～1 265 nm，对应的尺度分别为第4、5、6、6、7、7、7、9、10尺度；潮土的敏感区域主要分布在2 117～2 173 nm、2 252～2 282 nm、1 548～1 609 nm、431～542 nm、1 035～1 167 nm、1 823～2 090 nm、400～537 nm，对应的尺度分别为第7、7、8、9、9、9、10尺度；棕壤的敏感区域主要分布在419～525 nm、607～737 nm、1 621～1 694 nm、1 751～2 125 nm、417～507 nm、717～829 nm，对应的尺度分别为第7、7、7、7、8、8尺度。整体而言，不同类型土壤的SOM敏感波段基本上仍保持在相同的位置，并未发生大幅变动，且主要分布在近红外范围内，可能是因为东北旱作农田SOM主要源于农作物残体，其由糖类化合物、含氮化合物、纤维素等组成，这些成分中C—H键、C—O键、N—O键、N—H键等的光谱响应波段位于近红外区域[24]。

2.4 土壤有机质敏感波段筛选

以黑钙土原始光谱为例，采用sCARS算法筛选SOM敏感波段的过程见图6，可以看出，随着迭代次数的增加，所保留的敏感波段数量逐渐减少，RMSECV呈现先减小再增大的变化趋势，当运行次数为73次时，RMSECV最小(4.86 g/kg)。这是由于在1～73次变量筛选运行过程中，不断剔除与SOM含量相关性较弱，对建模结果影响不大的波段，使RMSECV减小；而在73次之后，开始剔除与SOM含量相关性强的波段，导致RMSECV增大。当运行次数为73次时，得到的敏感波段子集最佳，共筛选了13个敏感波段。图7为筛选后13个黑钙土敏感波段的分布情况，其主要分布在小于1 200 nm和大于2 200 nm的波段。不同处理筛选后敏感波段的数量(表2)仅占初始波段的0.07%～2.85%，有效去除了冗余的信息变量，提取了与SOM相关的重要特征信息变量，减少了参与建模的变量数量，降低了模型的复杂度。

图6 sCARS敏感波段筛选过程Fig.6 Screening process of sensitive bands by sCARS method

图7 sCARS筛选后黑钙土敏感波段分布Fig.7 Distribution of sensitive band of chernozem by sCARS method

表2 土壤有机质含量PLSR反演模型Tab.2 PLSR prediction model of soil organic matter content

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

将筛选后的敏感波段作为反演模型的输入变量，采用PLSR方法构建SOM含量反演模型，进行内部交叉验证(表2)。总体来看，各类型土壤的有机质反演效果均较理想，仅黑土原始光谱(R)、倒数对数(lg(1/R))和连续统去除(CR)模型验证R2低于0.50。黑土、黑钙土、潮土和棕壤最大R2分别达到0.83、0.88、0.93和0.93，潮土和棕壤的模型精度及稳定性略优于黑土和黑钙土。

比较分析4种处理的模型效果，一阶微分(R′)和连续小波变换(CWT)模型的效果明显优于倒数对数和连续统去除，与相关性分析结果一致。不同类型土壤建模集最佳模型均为CWT模型，验证集黑土、潮土的最佳模型为CWT模型，黑钙土、棕壤一阶微分模型验证效果最好。与原始光谱相比，倒数对数处理后，仅潮土的R2提升了0.04，RMSE下降了0.46 g/kg，其他类型土壤的模型精度和稳定性都略有降低；连续统去除处理后，不同类型土壤的R2提升了0.01～0.06，模型效果得到了一定程度提高，但效果不明显。而一阶微分和连续小波变换处理后，不同土壤类型的模型预测能力和稳定性均有较大提高，一阶微分处理后，建模集、验证集的R2最大提升了0.11和0.16，RMSE最大降低了1.33、1.44 g/kg；连续小波变换处理后，建模集、验证集的R2最大提升了0.13、0.28，RMSE最大降低了2.48、2.40 g/kg。这是因为连续小波技术不仅能有效地挖掘土壤光谱内隐含的有效信息，还能够压制噪声干扰的负面影响，提升光谱信息的信噪比，从而提升对SOM的估测精度与稳定性，而一阶微分处理虽然有效消除了基线和低频噪声的干扰，但同时也引入了高频噪声。

3 结论

(1)东北旱作农田黑土和黑钙土有机质含量丰富，光谱反射率较低；潮土和棕壤虽然有机质含量均值相近，但受氧化铁含量等其他因素的影响，潮土的光谱反射率低于棕壤。

(2)一阶微分对于分离重叠样品、消除平行噪声干扰有一定效果；CWT不仅可以抑制背景和噪声的干扰，还能够挖掘土壤光谱内隐含的有效信息，提升光谱信息的信噪比，可以有效提升土壤光谱与有机质含量的相关性。

(3)sCARS算法能够提取与SOM相关的重要特征信息变量，去除冗余、重叠的光谱信息，提高建模效率。

(4)CWT构建的SOM含量反演模型预测精度和稳定性最佳，黑土、黑钙土、潮土和棕壤R2分别达到0.83、0.88、0.93和0.93，CWT技术结合sCARS算法，能够较好地反演东北旱作农田典型土壤类型的SOM含量，为SOM含量的高光谱预测提供了新途径。