基于连续小波变换的土壤有机质含量高光谱反演

陈昊宇, 杨光, 韩雪莹, 刘昕, 刘峰, 王宁

(内蒙古农业大学沙漠治理学院， 内蒙古自治区风沙物理与防沙治沙工程重点实验室， 呼和浩特 010010)

精准农业作为目前农业发展的主要方向，是一种基于信息和知识管理的现代化生产系统，主要是通过3S(GPS、GIS和RS)技术与现代农业相结合，最大限度地提高农业生产力。所以快速、无损、精确地获取土壤中水分、养分的空间分布成为了实现精准农业的关键环节，近年来，光谱分析在土壤化学分析领域得到了迅猛发展，为实现土壤养分的快速诊断提供了新思路[1]。有机质是土壤养分供应能力和肥力的重要指标之一，在全球碳循环中发挥着重要作用。因此，快速准确地估测土壤有机质含量对于发展精准农业具有重要意义[2]。

传统的土壤有机质测定方法虽然精度比较高，但周期较长、成本较高，只能达到瞬测量，很难进行长时间大面积测量。高光谱遥感具有波段多、波段窄、信息丰富和实时高效等特点，为快速测量土壤有机质含量提供了一种新的方法和手段[3]。

国内外已经有大量研究表明，通过对光谱数据进行不同的数学变换(主要通过对光谱进行倒数、对数、微分、平方根、吸收峰深度、包络线去除等方法)可以有效提高光谱数据与土壤有机质含量之间的相关系数，有效筛选出光谱信息中的敏感波段[4]。现在各学者主要将研究重心放到了模型建立上[5]，普遍运用的线性模型有多元逐步回归与偏最小二乘回归[6]；常见的非线性模型包括BP神经网络[7]、支持向量机[8]、决策树[9]等，而且随着非线性模型算法的逐步改良与完善，在土壤有机质含量估算中已经成为不可取代的一部分。随着小波算法的改进与发展，最初仅运用于植物叶绿素、冠层成分含量预测中[10-11]，目前已成为土壤养分预测的热点问题[12-13]，连续小波变换是目前被广泛应用的一种方法。王祥浩[14]选择土地裸露地区为样区，利用神经网络算法对光谱连续小波变换、一阶导数、对光谱的平均值处理、光谱背景及深度4种方法建模，模型结果表明，小波变换方法得到的神经网络模型精度最高；包青岭等[15]选择渭干河-库车河三角洲具有代表性的干旱区绿洲为研究区，对光谱进行8层分解，结果表明小波变换不同分解层，从低频到高频范围内与土壤有机质含量的相关性呈现先减后增的趋势，结合随机森岭模型可以对干旱区土壤有机质含量进行有效的估算；王延仓等[16]以北京东部区潮土为例，对不同梯度重采样的光谱进行连续小波变换后，利用偏最小二乘法建立模型，结果表明连续小波分析算法可深入挖掘土壤光谱内的有益信息，提升对有机质含量的估测能力，与土壤高光谱反射率相比，经连续小波技术处理后，模型精度得到了有效的提升；叶红云等[17]同样针对干旱区土壤，通过对两种常用光谱变换R′、Ln(1/R)进行连续小波变换建立偏最小二乘模型，结果表明连续小波变换不会因人类干扰程度的提高而使模型精度大幅度降低，更加适用于干旱区有机质含量的预测；林鹏达等[18]通过解决黑土有机质高光谱野外反演的困难，同样证明了连续小波变换可有效提升模型精度。小波技术在土壤有机质高光谱反演研究中逐渐趋于成熟，但目前学者的研究多数都在同一土壤类型下或同一区域内，对于不同土壤类型及土地利用下土壤有机质高光谱反演是否存在影响的研究目前并不多。本文研究区内土壤类型主要包括3类：沙壤土、栗钙土、盐碱土，且部分区域土壤盐渍化程度严重，导致土壤养分空间分布上存在较大差异，取样表层土地利用类型主要包括：耕地、林地、草地、盐渍地、荒地。

通过对原始光谱(R)、原始光谱倒数(1/R)、原始光谱对数(LnR)以及原始光谱一阶微分(R′)4种不同情况进行连续小波变换，利用BP神经网络以及支持向量机2种模型，探究了不同土壤类型与不同土地利用类型下是否会对土壤有机质高光谱反演模型产生影响，小波变换前后土壤有机质反演模型的精度，旨为区域土壤有机质含量监测及实现精准农业提供理论与技术支持。

1 材料与方法

1.1 研究区概况

托克托县隶属于内蒙古自治区呼和浩特市，位于自治区中部、大青山南麓、黄河上中游分界处北岸的土默川平原上(图1)。地理坐标东经111°2′30″—111°32′21″、北纬40°5′55″—40°35′15″，总面积1 409.67 km2，平均海拔1 117 m，属于温带大陆性干旱气候，年均气温7.3 ℃，年均降雨362 mm。托克托县耕地总面积达400 km2，其中古城镇、新营子镇和五申镇的耕地较多，占全县耕地面积的60%以上[19]，主要作物包括小麦、玉米、莜麦。工农业及生产生活用水主要来源于大黑河和黄河水资源，整个地形以大黑河为轴，呈现由丘陵向平原过渡的趋势，地势为东南高、西北和西南低。东南向西北土壤类型依次为栗钙土、砂壤石灰性冲积土、盐渍化石灰性冲积土[20]，土壤类型的不同导致土壤养分存在差异性分布。植被类型从西向东依次为草甸草原、干草原和退化灌丛草原分布。以Landsat8OLI影像为基础数据源，运用人工目视解译与BP神经网络分类法得到托克托县2019年7月份土地利用数据，其中耕地面积最大为730.12 km2，占51.79%；林草地338.7 km2，占24.02%；盐碱地141.1 km2，占10.00%。详细土地利用空间分布见图1。

图1 土样采集点及土地利用空间分布

1.2 研究方法

1.2.1土样采集与处理土壤样本点均匀地分布在托克托县境内，采集方法为五点采样法，采集深度为0—20 cm，共采集120个点。采集的土样置于通风干燥室内进行自然风干、研磨，过10目筛，进行土壤光谱测定；过100目筛，采用重铬酸钾外加热法进行土壤有机质含量测定。

1.2.2光谱测量及光谱处理土壤光谱于暗室内测量，采用SVC HR-1024(北京东方佳气科技有限公司)便携式光谱仪，光谱范围在350～2 500 nm。在350～1 000 nm波段之间光谱分辨率≤3.5 nm；在1 000～1 850 nm波段之间，光谱分辨率≤9.5 nm；在1 850～2 500 nm波段之间，光谱分辨率≤6.5 nm。光源采用与太阳光接近的50 W卤素灯，将土壤样品放入深2 cm、宽10 cm的黑色器皿内，用直尺将土壤表面刮平，探头距离土样10 cm，光源距离土壤表面30 cm，天顶角为15°。测量前用白板进行标定，每个土样采集5条光谱作为该土样的光谱数据。

由于受噪音与仪器暗电流的的影响，导致光谱数据混入噪音等信息，因此删除350～399 nm和2 400～2 500 nm的波段，采用五点平滑法对光谱进行平滑处理，并将光谱重采样至5 nm，同时对原始光谱(R)进行一阶微分(R′)、倒数(1/R)、对数(LnR)等传统数学变换。

1.2.3连续小波变换采用连续小波变换，并用Mexh小波母函数对原始光谱、原始光谱的倒数、对数、一阶微分进行10层小波变换，生成一系列小波系数。

(1)

式中，a为伸缩因子，b为平移因子，λ为土壤高光谱数据的波段数。

(2)

式中，f(λ)为土壤光谱反射率，小波系数Wf(a，b)包含二维，分别为波长(350～2 500)与分解尺度(1，2，3…10)， 故小波系数行为尺度数，列为波长数的矩阵[16]。

1.2.4模型及精度验证采用BP神经网络与支持向量机模型(support vector machine，SVM)建立土壤有机质预测模型，支持向量机采用线性核函数，相对于径向基函数(radial basis function, RBF)来说计算高效,不易过拟合。BP神经网络的迭代次数设置为1 000，学习率0.01，训练的均方根误差(root mean square error，RMSE)小于0.001。

依据相关系数筛选的特征波段以及小波系数作为自变量，土壤有机质含量为因变量，分别建立模型，模型精度采用决定系数(R2)、均方根误差(RMSE)、相对分析误差(relative percent deviation，RPD)以及1∶1线共同评价。R2表征模型的稳定性，越接近于1模型越稳定，拟合程度越好。均方根误差(RMSE)用来检验模型的预报能力，RMSE越小则表明模型的估测能力越好。RPD是样本的标准差与RMSE的比值，RPD<1.4时，模型无法对样品进行预测；1.4≤RPD<2时，模型效果一般，可以用来对样品进行粗略评估；RPD≥2时，模型具有极好的预测能力。1∶1线表示实测值与预测值构成的点偏离y=x线的程度[21]。

2 结果与分析

2.1 土壤有机质含量统计分析

建模样品集、不同土地利用方式、不同土壤类型下土壤有机质含量描述性统计见表1。本研采样点内土地利用方式主要包括林地、草地、耕地、盐渍地，土壤有机质在草地内均值含量最大(0.80%)，其次为林地(0.72%)、耕地(0.67%)、盐渍地有机质含量最低(0.63%)；土壤有机质含量最大值位于耕地(1.28%)，最小值位于林地(0.19%)。采样点内主要土壤类型为栗钙土、沙壤土、盐碱土，沙壤土有机质含量最高(0.77%)，其次为盐碱土(0.68%)和栗钙土(0.67%)，土壤有机质含量最大值位于沙壤土内(1.28%)，最小值位于盐碱土内(0.19%)。

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

2.2 土壤反射光谱特征

对R、1/R、LnR、R′进行小波变换，变换结果如图2所示，R、1/R、LnR光谱曲线较为平滑，分解曲线随波峰波谷变化.R′其光谱曲线并不规则存在较多波峰波谷，分解小波系数与前三者不同。R、1/R、LnR、R′分解后，小波系数均随分解尺度的增加而增加，同时可以看出，由Mexh小波母函数进行的连续小波变换，对于光谱波峰与波谷有较高的敏感性，对于放大、挖掘光谱信息有着显著的作用。

图2 连续小波变换光谱特性

2.3 相关性分析

2.3.1不同导数变换光谱与土壤有机质含量相关性土壤有机质含量与光谱相关性曲线及敏感波段见图3。R与土壤有机质含量呈负相关关系(相关系数r=-0.463)，主要集中于735～780 nm处波段；1/R与土壤有机质的相关性则与R相反，呈正相关关系(r=0.462)，集中于600～800 nm与1 800～2 200 nm处波段；LnR的相关性曲线图与R相关性曲线类似，总体呈现负相关关系，相关系数(r=-0.465)，主要集中于745～795 nm处的波段；R′相关性在500 nm(r=-0.589)与1 400 nm(r=-0.411)处为负相关，在8 00 nm(r=0.408)与1 380 nm(r=0.412)处为正相关，相关系数曲线变换趋势与前三者不同，呈无规律变化。

图3 土壤光谱相关性曲线及敏感波段

2.3.2不同分解尺度小波系数与土壤有机质含量的相关性图4为不同光谱变换方式经过连续小波变换后与土壤有机质含量的相关系数矩阵图，其中红色代表相关性高的区域，蓝色代表相关性低的区域。R在800～1 000、1 400～1 600 nm处相关性明显增加，在500、800、2 200 nm波段处相关系数达到最大值(r=0.667)；1/R在800～1 200 nm处相关系数达到最大值(r=0.552)，在2 400～2 500 nm处相系数达到0.4，受噪音和仪器本身的影响，此波段的相关系数不进行相关性参考；LnR在分解尺度1下相关性较低，在2～10尺度下，相关性出现最大值(r=0.664)；R′相关性主要集中在500～900、1 200～1 600、2 100～2 300 nm处。筛选的敏感波段与尺度如表2所示。有效的光谱信息主要存在于低分解尺度，随分解尺度的增加呈递减趋势，相关性最大值较未处理前分别增加了0.204、0.09、0.199、0.252，对于挖掘潜在光谱信息有着重要意义。

图4 土壤有机质与小波系数相关性

表2 筛选的敏感波段

2.4 土壤有机质高光谱模型建立

2.4.1BP神经网络预测模型采用BP神经网络构建反演模型，结果如表3所示。未进行连续小波变换处理的模型中，BP-R与BP-R′效果较好，R2分别为0.69和0.73，RPD为1.45与1.53，模型能粗略估算土壤有机含量，BP-LnR与BP-1/R样本外预测能力较差，同时RPD未达到1.4以上，不能对土壤有机质未能进行有效预测；连续小波变换处理之后的模型，仅BP-CWT-1/R模型RPD未达到预测水平，其余3种模型R2与RPD较未处理前均有所增加，RMSE均减少，其中BP-CWT-LnR模型预测效果较好，RPD达到2.12可以有效地对土壤有机质进行预测。将BP-CWT处理的4个模型的实测值与预测值进行1∶1线分析。由图5可知，除BP-CWT-1/R模型外，其余模型的实测值与预测值样点基本分布在1∶1线附近，BP-CWT-LnR效果较为明显，且估算精度高，可较好地进行土壤有机质含量的估算。

表3 土壤有机质BP神经网络估测模型结果

图5 BP-CWT模型土壤实测值与预测值对比

2.4.2支持向量机预测模型SVM构建反演模型，结果如表4所示。未经过连续小波处理的光谱特征波段未能较好地对土壤有机质进行预测反演，经过CWT后模型SVM-CWT-R与SVM-CWT-R′预测结果较之前有较大的提升，R2分别达到了0.50与0.56，二者RPD均达到1.4以上，可以粗略地对土壤有机质进行预测。同时根据图6，SVM-CWT模型进行1：1线分析，二者实测值与预测值分布情况在4种模型下较好，虽然模型SVM-CWTLnR分布同样较为集中，但其样本外预测情况较差(RPD=1.38)，综合考虑不对其进行土壤有机质预测。结合表3和表4的结果分析，连续小波变换能够有效地提升模型精度与模型泛化能力，对于光谱信息挖掘有着重要意义，BP神经网络与支持向量机对CWT-R与CWT-R′都能够提升R2减少RMSE，可对土壤有机质做出较好的预测。虽然BP神经网络与支持向量机在处理非线性回归问题中有较强的能力，但本身模型中存在不稳定性，对模型的环境设置同样要求较高，所以未能对所有数据集进行良好的预测。

表4 土壤有机质支持向量机估测模型结果

图6 利用SVM-CWT模型土壤实测值与预测值的对比

3 讨论

本研究采用连续小波变换对光谱进行处理，用BP神经网络与支持向量机(SVM)两种模型对土壤有机质含量进行反演预测。未经过连续小波变换前，R、1/R、LnR、R′与土壤有机质的相关系系数最大值分别为-0.463、0.462、-0.465、0.589，可以看出，R′与土壤有机质的相关系数最高，与吴倩等[22]、张新乐等[23]的研究结果相同；经过连续小波变换后，CWT-R、CWT-1/R、CWT-LnR、CWT-R′相关系数最大值分别为0.667、0.552、0.664、0.662，较之前分别增加了0.20、0.09、0.19、0.07。王延仓等[1]、于雷等[4]、叶红云等[17]等同样证明连续小波变换可有效提高与土壤有机质含量的相关系数。不同分解尺度对于光谱数据的深度挖掘有着重要意义，本研究只利用Mexh小波母函数进行处理，未对其他函数进行考虑，分解层数同样是根据前人经验所得[4，10]，小波技术的研究与发展仍然有很大的探索空间。

相对于两种模型来看，未进行连续小波处理的支持向量机模型中，只有SVM-R′模型R2最高达到0.43，其余三者均未到达0.4。综合多种模型评价方法，由于其RPD未达到1.4以上，无法对土壤有机质含量进行预测。经过连续小波处理后，各模型的R2有明显提高，其中SVM-CWT-R与SVM-CWT-R′模型效果较好，R2分别提高了0.29、0.13，RPD达到1.62与1.53实现了对土壤有机质有效的预测，但预测结果较BP神经网络较低。在BP神经网络预测模型中，未进行连续小波变换前，BP-R与BP-R′预测效果较好，R2达到0.69与0.73，RPD为1.45与1.53；进行连续小波处理后，除SVM-CWT-1/R模型未到达预测效果，其余3种模型预测结果较之前均有明显改善，可实现对土壤有机质较好的预测，其中BP-CWT-LnR预测模型效果最佳R2达到0.76，RPD达到2.12。根据1：1线分析图也可看出，其实测值与预测值分布较为集中，于雷等[4]、叶红云等[17]、林鹏达等[18]同样通过连续小波变换有效提升了模型的精度与泛化能力。

针对土壤有机质高光谱反演研究中，姚聪[24]对耕层土壤通过BP神经网络与支持向量机模型，反演精度R2分别为0.42与0.67；叶红云等[17]采用连续小波变换对干旱区土壤有机质反演，模型精度R2=0.75、EMSE=0.71；谢文[25]在森林土壤有机质反演研究中，BP神经网络模型R2=0.78、EMSE=0.77，支持向量机模型R2=0.87、EMSE=0.76。本研究对耕地、林草地、盐碱地、栗钙土、沙壤土、盐渍土等不同土地利用类型与土壤类型进行综合反演，最佳反演模型为BP-CWTLnR，R2=0.76、EMSE=0.15、RPD=2.12，与前人研究的结果基本相符，证明通过连续小波变换处理，不同土壤类型与土地利用类型未对土壤反演模型精度产生影响。所以采用连续小波变换进行光谱数据挖掘，采用BP-CWT-LnR神经网络建立反演模型，可对不同土地利用于土壤类型条件下土壤有机质高光谱反演提供一定的理论支持与应用价值。