1 引言

土壤有机质主要由各种动植物的残体、微生物体及其分解和合成的各种有机质组成,虽然其仅占土壤总量的很小一部分,但在土壤肥力上的作用却十分重要。不同土壤类型和区域由于气候以及人类干扰程度的不同,使得有机质含量分解程度有所不同,其光谱反射特性也存在差异。

灰钙土是北疆典型的土壤类型之一,其有机质质量分数大多<2%,研究表明有机质含量较低会大幅减弱其对光谱反射率的吸收^[1-2]。为提高土壤有机质含量的预测精度,研究者多采用微分、倒数、对数、去包络线等变换对土壤光谱进行研究并取得了较好的效果^[3-7]。连续小波变换(CWT)近几十年来在高光谱遥感领域以其在信号去噪和数据压缩方面独特的优势为地物光谱信息提取提供了新途径,已有研究表明其比常用的光谱数学变换预测效果更优^[8-14],但目前对于干旱区土壤有机质的连续小波变换研究尚未涉及。

生态系统的干扰已经成为一种广泛存在的现象,其中人类干扰形式主要有对森林和草原植被的砍伐和开垦、狩猎、捕捞等。在我国干旱区、半干旱区,由于人类放牧、农业活动等的干扰活动,人类将大面积的荒地开垦为耕地、林地。人类活动对土壤原有生态平衡造成的破坏和改变,使土壤水分、有机质、盐分等指标发生变化^[15-17]。但目前很少有针对干旱区不同人类干扰程度下的土壤有机质进行预测的研究。因此,如何利用光谱准确估算其含量是尚需解决的问题。

本文以新疆北部阜康荒地为研究对象,从人类干扰程度出发,结合植被特征、土地利用方式等指标,将研究区划分为轻度干扰区、中度干扰区和重度干扰区3种类型,以CWT为研究手段估算其土壤有机质含量,研究人类不同程度的干扰对土壤有机质含量的影响及其预测精度,以期为干旱区土壤资源合理规划以及农业产业结构的调整与布局提供借鉴和参考。

2 材料与方法

2.1 研究区概况

研究区位于新疆北部的阜康市(87°44′E~88°46′E,43°29′N~45°45′N)。地势南高北低,平均海拔高度为452 m。气候为典型的温带大陆性气候,夏季高温,冬季严寒,光能资源充足,热量资源丰富,降水稀少且空间分布不均,年降水量为164 mm,年蒸发潜力约为2000 mm ^[18]。

实验将研究区分为I、II、III区。I区为轻度干扰区,因距离人类居住地较远,并未开发利用。该区人类活动很少,偶有人类进入,基本保持其原有风貌;II区为中度干扰区,处于新疆建设兵团102团附近,大部分为废弃耕地,人类对其有一定程度的干扰; III区为重度干扰区,位于研究区的西南部,目前人类对其开发利用程度较大,主要的开发利用方式为林地(包括人工梭梭林地、榆树林地、育苗地),尚处于开发利用的初期,如表1所示。

参照余作岳等^[19]提出的人类干扰分类体系,本文所研究的I区长期以来的主导干扰因子是人类偶尔进入该区,并未对该区原有植被类型造成干扰,即几乎未引起土壤中组分的变化,属于无效干扰;而II区、III区经过人类开垦、植树等,植被类型发生改变,正常演替序列被打破,土壤中的盐分、有机质等组分也随之发生变化,即这两种干扰改变了生态系统的正常动态发展,属于有效干扰。因此,以轻度干扰区的各项指标作为对照,比较分析中度干扰和重度干扰区土壤有机质含量的变化等信息。

2.2 实地调查与采样

研究于2017年10月进行了土壤样本数据采集。在Ⅰ区由南向北布设5条间距500~700 m的东西向采样线,每条采样线布设6个采样点;Ⅱ区相同方向布设5条间距400~600 m的采样线,每条采样线上选择5~7个具有代表性的采样点;Ⅲ区面积较小,布设6条间距200~400 m的采样线,每条采样线5个采样点。上述每个区均有30个采样点,由此构成一个由90个采样点组成的空间网格,并分别利用全球定位系统(GPS)对其进行定位。采集的野外数据主要包括:

1) 野外土壤高光谱数据的采集

光谱测量采用美国ASD公司生产的便携式FieldSpec 3Hi-Res光谱仪,波段值为350~2500 nm,其中350~1000 nm 光谱采样间隔为1.4 nm,光谱分辨率为3 nm;1001~2500 nm范围内光谱采样间隔为1.1 nm,光谱分辨率为8 nm。由于光谱测量采用野外实测,为了减少太阳高度对光谱采集造成的不利影响,实验均选在当地时间11:00~15:00、晴朗少云、无风的天气进行。每次采集光谱前对光谱仪进行白板校正以去除暗电流的影响。采用25°视场角探头,且距采样地面15 cm处垂直角度,在每个采样点按梅花桩对表层土壤原始光谱进行采集,每个位置重复测量10次,得到的50条光谱曲线的平均值即为该采样点的实测光谱值,共测定90个采样点的光谱曲线。野外测量时观察每个采样点的实测光谱值,对出现异常的光谱曲线给予删除并重新测定,以便获得更加准确的实测光谱反射率。

2) 土壤有机质含量的采集

土样的采集与高光谱数据的测定同步进行,选择地势平坦、能代表样点周围区域特征的地点作为采样单元,每个采样点分别从周边1 m范围内不同的地点挖取3个0~10 cm的土坑采样,并混合均匀放入密封袋中,封口标记编号,用手持GPS定位后记录下经纬度。将采集的样品带回实验室经自然风干、去除杂质后,研磨过 1 mm筛,送至中国科学院新疆生态与地理研究所,由专业人士采用重铬酸钾容量—外加热法测定^[20]有机质含量。

2.3 高光谱数据处理

由于光谱仪器波谱两端的边缘波段存在较大的仪器噪声,需要将350~399 nm和2451~2500 nm数据剔除,运用OriginPro 9.0对光谱反射率进行Savitaky-Golay(9点)平滑去噪处理,再通过观察研究区90个采样点的土壤光谱曲线,找出野外实测光谱因受水分吸收带和大气的影响而在1350~1450 nm和1830~1950 nm附近出现的异常,并剔除在1400、1900 nm附近所产生的较大噪声。

为了进一步减少外界其他干扰因素对原始光谱反射率的影响,查找土壤有机质的敏感波段,除了对土壤样品的光谱反射率R进行分析外,还对光谱曲线做了一阶微分R′及lg(1/R)的变换。

2.4 连续小波变换

CWT是一种线性变换方法,是将原始光谱信号在连续波段上进行分解,分解后的系数与原始的光谱波段一一对应^[21]。其变换公式为

Wf(a,b)≤fψa,b≥∫-∞+∞ftψa,b(t)dt,(1)ψa,b=1aψt-ba,(2)

式中:f(t)为研究区野外实测高光谱数据;t为光谱波段;ψ_a,b为小波基函数;a为尺度因子,b为平移因子;W_f(a,b)为小波系数。

2.5 模型的建立与检验

模型的精度检验主要通过判定系数R²、均方根误差R_RMSE和相对分析误差R_RPD检验实测值和估测值的拟合效果。R²越大,R_RMSE越小,则模型精度越高^[22]。Viscarra等^[23]对模型R_RPD做了如下分类:R_RPD<1.0,预测能力极差;1.0<R_RPD<1.4,预测能力较差;1.4<R_RPD<1.8,预测能力一般;1.8<R_RPD<2.0,预测能力较好;2.0<R_RPD<2.5,预测能力极好。

3 结果与分析

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

选用光谱理化值共生距离(SPXY)算法^[24]分别计算Ⅰ、Ⅱ、Ⅲ区各个样品有机质含量之间的欧氏距离,按照3∶2的比例划分为建模集和检验集。研究区土壤有机质含量总体较低,且随着干扰强度的减少,经营方式逐渐回归自然状态,有机质含量逐渐增加,变异系数(CV)逐渐减少,如表2所示。其中Ⅰ区原始植被状态保存完好,主要分布梭梭、柽柳、猪毛菜、盐爪爪等,植被覆盖度≥30%,平均有机质含量为10.257 g/kg,明显高于其他区,这是因为受人类活动干扰的影响,土地利用方式多样,主要有人工林、枸杞育苗地以及其他农业土地利用类型,使土壤有机质空间异质性与其他区相比得到加强,CV也相应地增加。此外,人类干扰活动加速了土壤中水分的散失,加上干旱区降水较少,使得盐分集聚土壤表面,不利于植物的生长,导致有机质含量较低。

3.2 人类不同干扰程度下的野外反射光谱特征

由于对土壤光谱外部影响因素较多,可以采用野外测定的光谱数据进行分析,并剔除明显受到干扰的波段(图1)。Ⅰ、Ⅱ和Ⅲ区的土壤光谱曲线形态基本相似,且随着人类干扰程度的增加,土壤光谱反射率值随之升高。光谱曲线在可见光范围反射率较低,随着波长的增加呈快速增长趋势,在近红外波段范围内趋于平缓,并在1400、1900 nm附近存在较为明显的水分吸收谷。由于在近红外波段光谱会受到NH、CH等基团的影响^[25],3条光谱曲线的吸收深度和吸收面积均存在差异。

图 1. 人类不同程度干扰下的土壤光谱曲线

Fig. 1. Soil spectral curves under different human disturbance degrees

下载图片 查看所有图片

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

分别对土壤野外原始光谱反射率进行一阶微分和倒数对数处理,并进行相关性分析。Ⅰ、Ⅱ和Ⅲ区原始光谱与土壤有机质含量的相关性均未通过0.01显著水平检验,且通过0.05显著水平的波段也较少(图2)。这是因为本研究采用的是野外实测光谱,无法避免植被、土壤质地等外界因素的影响,使得光谱信息受到一定的干扰。而且原始光谱反射率在可见光区差异小,因光照引起的乘性因素会产生较大的噪声^[26]。总体来看,Ⅰ、Ⅱ、Ⅲ区原始光谱反射率与有机质含量相关性较低,因此,如何将较为微弱的光谱信息从高光谱曲线中提取出来,成为干旱区土壤有机质含量准确估测的关键。

图 2. 土壤有机质含量与光谱反射率及其变换的相关性分析。(a) I区;(b) Ⅱ区;(c) Ⅲ区

Fig. 2. Correlation analysis among soil organic matter content, spectral reflectance, and its transformation. (a) Zone Ⅰ; (b) zone Ⅱ; (c) zone Ⅲ

下载图片 查看所有图片

通过对原始光谱反射率进行微分R′和倒数对数lg(1/R)数学变换后,Ⅰ、Ⅱ、Ⅲ区通过0.05水平检验的波段增多,且有小部分波段通过了0.01显著性水平。这是因为经过光谱微分变换可以平缓背景干扰或消除基线漂移,有助于吸收特征的提取;而光谱反射率经过倒数对数变换之后,不仅能够增强可见光部分的光谱微弱差异,还可以减少因光照、地形等外界条件引起的随机因素的干扰^[27]。以0.05显著水平以及相关性最大的原则,分别选择其中相关系数较大的5个敏感波段作为反演模型的输入变量,见表3。

3.4 基于CWT的相关性分析

由于土壤高光谱曲线的吸收特征与高斯函数近似,因此,选取Gaus4函数作为小波基函数,利用Matlab R2017b对不同人类干扰下的原始光谱反射率进行CWT。为了减少数据冗余,将CWT尺度设为2¹,2²,2³,…,2¹⁰,同时将生成的小波系数与土壤有机质含量进行相关性分析,如图3,图4,图5所示。

3.4.1 轻度干扰区CWT相关性分析

采用CWT对Ⅰ区原始光谱反射率进行变换分解,生成的小波系数与土壤有机质含量的决定系数较高的区域主要出现在420~440、460~490、740~770、1635~1670和2110~2120 nm(图3)。将各分解尺度下的决定系数进行降序排列并对比分析,最终选择432、480、761、1645、2115 nm作为敏感波长,其对应的尺度分别为3、4、9、8和2。其中决定系数最高为0.48,比原始光谱反射率R最大提升约0.33。这主要是因为通过CWT处理可以放大光谱吸收,能够有效提取光谱中的微妙信号,更好地捕捉土壤中的有机质光谱特征。

图 3. 轻度干扰区R小波系数与土壤有机质含量的相关系数图

Fig. 3. Correlation scalogram between R wavelet coefficient and soil organic matter content in mild disturbance zone

下载图片 查看所有图片

3.4.2 中度干扰区CWT相关性分析

将Ⅱ区原始光谱反射率进行CWT分解后,小波系数与土壤有机质含量的R²同样得到了提高(图4)。决定系数最大可达0.45,与原始光谱反射率和土壤有机质的决定系数相比提高约0.28,其高值区主要集中在450~480、580~600、690~710、1670~1690和2320~2340 nm,结合排序对比分析得出最大值分别出现在465、584、709、1682、2334 nm处,对应的尺度分别为5、3、6、1、4。

图 4. 中度干扰区R小波系数与土壤有机质的相关系数图

Fig. 4. Correlation scalogram between R wavelet coefficient and soil organic matter in moderate disturbance zone

下载图片 查看所有图片

3.4.3 重度干扰区CWT相关性分析

Ⅲ区原始光谱反射率通过CWT结合土壤有机质含量数据得到二者的二维决定系数图(图5),其决定系数较高的区域分解尺度主要分布在1、4、5、6和7,当分解尺度为5时,最大的决定系数R²为0.40。敏感波段主要集中在585~610、480~520、2100~2130、1310~1320和1790~1810 nm,分别在509、552、2110、1315、1798 nm达到最大值。

图 5. 重度干扰区R小波系数与土壤有机质含量相关系数图

Fig. 5. Correlation scalogram between R wavelet coefficient and soil organic matter content in severe disturbance zone

下载图片 查看所有图片

结果表明,针对野外采集的实测高光谱反射率,在有机质含量较低且受人类不同干扰程度下,通过CWT可以提高土壤高光谱中反映的土壤有机质微弱信息,有助于土壤有机质含量的准确预测。

3.5 土壤有机质高光谱估测模型的优选

基于上述分析,将R、R′、lg(1/R)选取的敏感波长以及CWT选取的小波系数为自变量,土壤有机质含量为因变量,构建偏最小二乘模型(表4)。对R进行R′、lg(1/R)、CWT变换后,在Ⅰ、Ⅱ、Ⅲ区无论是建模集的R²、R_RMSE,还是验证集的R²、R_RMSE均有不同程度的改善。以Ⅰ区CWT建模效果为例,R建模集和预测集的R²、R_RMSE分别为0.372、1.350和0.335、2.391,而CWT对应的数据分别为0.752、0.848和0.717、1.132。因此,对R进行R′、lg(1/R)、CWT变换处理,对于提高模型的预测精度以及模型的稳定性具有重要作用。综合比较3个区的4种光谱指标的估算精度,发现以CWT为自变量建模所得的建模集和验证集的R²和R_RPD均高于其他光谱指标的对应值,而R_RMSE则均低于其他光谱指标的对应值。以原始光谱反射率R进行建模时,Ⅰ、Ⅱ、Ⅲ区的R_RPD值均<1.4,预测能力较差;以R′、lg(1/R)为自变量所建模型的R_RPD值均未达到2,只能对有机质含量进行粗略的估计;而以CWT为自变量所建模型得到的Ⅰ、Ⅱ、Ⅲ区的R_RPD均>2,分别为2.150、2.090、2.013,表明其预测能力很好,可以对有机质含量进行精确的估测。此外,通过对比3个区CWT与原始光谱反射率R建模集和验证集的R²来看,Ⅰ、Ⅱ、Ⅲ区的CWT变换与R相比,建模集和验证集的R²分别提高了0.380、0.382,0.432、0.400,0.407、0.389,表明人类干扰活动的强弱并不会对CWT估算产生较大影响,该方法具有较好的适用性。

综合比较建模效果可以看出,Ⅰ区>Ⅱ区>Ⅲ区,这是因为Ⅰ区受人类活动干扰较小,大部分区域保持原始风貌,土壤有机质含量变异系数低,有利于有机质含量的预测;而Ⅱ、Ⅲ区由于不同地块受人类活动干扰的强度不同,使有机质分布变异性较大,预测精度受到影响。

利用检验样本将研究区R-CWT预处理的3个模型的预测值和实测值绘制散点图,如图6所示。研究区有机质含量实测值和预测值分布较为集中,表明原始光谱经过CWT后能够更好地筛选出土壤有机质敏感信息,利用其建立的土壤有机质含量高光谱估测模型具有较高的精度,更适应用于干旱区土壤有机质含量的预测。

图 6. 土壤有机质预测值与实测值散点图。(a) Ⅰ区;(b) Ⅱ区;(c) Ⅲ区

Fig. 6. Scatter plots of predicted and measured values of soil organic matter. (a) Zone Ⅰ; (b) zone Ⅱ; (c) zone Ⅲ

下载图片 查看所有图片

3.6 模型预测效果的评价

为了更加直观地观察模型预测效果,对研究区实测值与预测值进行Kriging插值分析(图7)。3个区土壤有机质含量预测值和实测值分布整体上具有较好的一致性。其中,Ⅰ区土壤有机质低值主要位于西北部,高值位于东南方向;Ⅱ区有机质含量从西南向东北方向递减,有机质含量高值主要分布在西南,低值主要位于西北和东北方向;Ⅲ区低值主要位于南部,北部有机质含量较高。研究表明,土壤有机质含量预测值与实测值的空间分布高度吻合,反映出CWT方法能够最大程度地提取高光谱中的有效信息,有利于提高干旱区土壤有机质含量的估算精度。

图 7. 土壤有机质含量实测值与预测值的Kriging插值图。(a) Ⅰ区;(b) Ⅱ区;(c) Ⅲ区

Fig. 7. Kriging interpolation plots of measured and predicted values of soil organic matter content. (a) Zone Ⅰ; (b) zone Ⅱ; (c) zone Ⅲ

下载图片 查看所有图片

4 结论

利用CWT对新疆阜康市干旱区人类不同干扰程度下的土壤有机质含量进行了估算,并与野外高光谱反射率R、光谱反射率一阶微分R′、光谱反射率倒数的对数lg(1/R)进行对比。主要得出以下结论:随着人类干扰程度的加强,有机质含量的空间变异性也逐渐增加,依次为Ⅰ、Ⅱ、Ⅲ区;Ⅰ、Ⅱ、Ⅲ区的R与土壤有机质含量相关性较低,均未通过0.05水平的检验,其中Ⅰ、Ⅱ区的R与土壤有机质含量呈负相关,Ⅲ区的则呈正相关;相比R、R′、lg(1/R),Ⅰ、Ⅱ、Ⅲ区原始光谱反射率经过CWT处理后得到的小波系数与土壤有机质含量的相关系数提高幅度最大,与R相比Ⅰ、Ⅱ、Ⅲ区CWT决定系数R²分别最大提高0.33、0.28、0.37;3个区利用CWT所建模型的预测值和实测值均具有较好的一致性,Ⅰ、Ⅱ、Ⅲ区R²分别为0.717、0.689、0.630,R_RPD分别为2.150、2.090、2.013,相对于原始高光谱反射率R所建模型的R²分别提高了0.382、0.400、0.389,且由预测能力较差上升为预测能力极好,均能很好地预测土壤有机质含量。