土壤盐渍化作为荒漠化的一种表现形式, 会引起生态、环境、社会和经济等一系列问题, 这在干旱、半干旱区表现的尤为明显^[1], 严重阻碍着区域生态文明建设的推进；另一方面, 盐渍土作为重要的后备耕地资源, 迫于人口剧增、粮食不足、环境恶化和生态破坏等的压力, 开发利用盐渍土的局势也已经迫在眉睫^[2]。要达到治理与利用大范围盐渍土的目的, 必须及时获取有关盐渍土的可靠信息。因此, 土壤盐渍化成因、空间分布、变化规律、监测方法、提高监测精度和盐渍化预警能力等成为地理学、生态学、农学等学者们的研究热点。

在土壤盐渍化动态变化监测中, 具有尺度大、效率高和破坏小等特点的遥感技术近年来被广泛应用。其中, 拥有更精细光谱波段信息的高光谱遥感对于提高土壤盐渍化监测精度提供了有效途径, 所以国内外众多学者利用高光谱数据进行了大量土壤盐渍化监测理论与定量反演模型的研究。在1984年, Clark等^[3]已经初步开始利用高光谱反射率研究地物特征, 为后续研究地物高光谱特征及高光谱遥感反演提供了理论基础。Csillag等^[4-6]通过对不同盐分含量土壤在400—2500nm范围光谱进行测定, 分析归纳了不同盐分含量土壤高光谱敏感吸收特征波段, 指出大气对土壤盐渍化的高光谱诊断特征产生一定的影响, 土壤盐分含量的反演应集中于波谱形状, 而与高光谱吸收带参数关系较少。陈皓锐等^[7]测定内蒙古河套灌区沙壕渠灌域沙壤土样品的高光谱曲线并进行平滑处理、相关分析和去包络线处理, 分别采用偏最小二乘回归法和稳健估计法构建土壤电导率的估算模型, 为大面积快速获取含盐土壤电导率和盐渍化特征提供了参考。姚远^[8]基于EM38测得的盐渍土电导率数据和高光谱仪测得的盐渍土高光谱反射率数据, 对高光谱反射率进行11种光谱变换后与电导率数据作相关分析, 选取相关性最好的变换形式及响应波段计算5种盐分指数, 并筛选最优高光谱指数, 以此建立区域土壤盐渍化监测模型。彭杰等^[9]利用土壤样品高光谱数据以及室内测定的盐分和电导率数据, 分析耕作土壤盐分含量与电导率之间的关系, 比较高光谱信息对二者的敏感性以及高光谱反演模型的精度, 发现对于耕作土壤而言, 电导率与盐分含量之间没有显著的相关性, 因此不能利用电导率数据进行土壤盐渍化的高光谱遥感监测。

综上所述, 以往对于土壤盐渍化高光谱定量估算的研究是通过全波段反射率及其对应的数学变换在一维层面上选取单个敏感波段或多个敏感波段, 并利用多种回归方法建立预测模型, 模型精度有待进一步提高, 且基于两波段优化组合算法进行光谱指数的波段二维层面优化的应用研究相对较少。基于此, 本研究以新疆艾比湖流域为研究靶区, 进行野外土壤采样, 室内土壤高光谱采集及理化分析等工作, 尝试运用土壤原始高光谱反射率及对应的5种数学变换, 对简化光谱指数(nitrogen planar domain index, NPDI)进行波段优化计算, 分析基于不同形式变换光谱的优化光谱指数与盐渍土电导率的相关性, 通过变量重要性准则(variable importance in projection, VIP)筛选最优高光谱参数, 并利用偏最小二乘回归(partial least square regression, PLSR)分析法建立土壤电导率高光谱定量估算模型, 为土壤盐渍化动态监测及星载传感器等相关研究提供科学支持和应用参考。

1 研究区概况

地处新疆博尔塔拉境内的艾比湖流域, 区域年平均气温7.7℃、年均降水量102 mm, 潜在蒸散量达1447 mm, 典型土壤为灰漠土、灰棕漠土及风沙土等, 是典型的干旱半干旱地区盐渍化代表区^[10]。艾比湖作为新疆最大的咸水湖, 湖水主要依赖地表径流补给, 近十几年来, 由于地表径流被引入灌区、渗漏地下和消耗于地面蒸发与植物蒸腾的量增多, 艾比湖入不敷出, 湖泊面积严重萎缩, 干涸湖底面积不断增加, 湖滨荒漠化及周边区域土壤盐渍化程度不断加剧, 加之阿拉山口全年8级以上大风达164d, 常年侵蚀大面积裸露湖床及盐壳, 使艾比湖日渐成为中国西北部沙尘暴、盐尘暴的主要策源地之一, 严重影响天山北坡绿洲生态文明建设与可持续发展。

2 实验材料与数据采集 2.1 土壤样品采集与分析

土壤样品的采集时间为2016年10月15日至10月23日, 针对艾比湖流域典型自然盐渍化土壤, 设置5 m×5 m样方(图 1), 利用GPS记录每个样方位置, 采用5点混合法进行样品采集, 土壤样品采样深度为0—10 cm, 共计57个土壤样品。土壤样品带回实验室并进行自然风干、研磨后, 通过2 mm孔筛分为两部分, 分别用于室内高光谱数据采集及土壤电导率和含盐量分析。其中, 土壤样本的电导率和全盐量测定方法参照《土壤农业化学分析方法》, 土壤电导率在土水比例为1 : 5的土壤悬浊液中利用德国WTW(WissenschaftlichTechnischeWerkstätten)公司生产的inoLab^Ⓡ Cond 7310台式电导率测试仪进行测定, 全盐量采用水浴烘干法进行测定。

2.2 高光谱测定及预处理

本研究通过美国ASD(Analytical Spectral Devices)公司生产的FieldSpec3型高光谱仪测定土壤室内高光谱数据, 波段范围为350—2500 nm。高光谱波段在350—1000 nm区间的采样间隔为1.4 nm, 在1000—2500 nm区间的采样间隔为2 nm, 全波段范围重采样间隔为1 nm。将处理好的土壤样品分别装入黑色盛样皿(直径12 cm, 深1.8 cm)中, 对装满的土壤进行表面平滑处理。测定样品的光源为50 W卤素灯, 测量时距土壤样品表面为50 cm, 光源的天顶角为15°, 待测土壤样品表面与探头的距离为5 cm。每测定两个土壤样品后进行一次白板校正, 每个土壤样品重复测定5次, 取5条高光谱曲线的算术平均值作为该样品的实际高光谱反射率^[11]。

将ASD FieldSpec3高光谱仪测定的土壤高光谱数据通过ViewSpec Pro软件处理与导出, 为了降低噪声引起的影响, 去除信噪较低的边缘波段(350—399 nm及2401—2500 nm), 利用Savitzky-Golay滤波方法对57个土壤样品的高光谱数据(400—2400 nm)进行平滑去噪预处理。

2.3 研究方法

本研究利用Li等^[12]提出的三波段光谱指数(NPDI)归一化比值算法, 其计算公式为NPDI=(R_NIR/R_RED-1)/[(R_NIR-R_RED)/(R_NIR+R_RED)], 为进一步发挥高光谱数据的优势, 将光谱指数从三波段组合降级到两波段组合, 计算公式变为(R_i+R_l)/R_l, 其中i与l作为土壤高光谱反射率。基于此, 本文使用自主开发的《高光谱数据两波段组合软件V1.0》(No: 2018R11S177501), 将两波段优化算法应用于全波段范围进行光谱指数波段优化。通过该算法计算57个土壤样品的原始高光谱反射率R以及对应的5种光谱变换(倒数变换、对数变换、对数倒数变换、平方根变换、分数阶微分变换)在波段400—2400 nm之间进行所有可能两波段组合的优化光谱指数(NPDIs)。分数阶微分变换处理在《高光谱数据分数阶微分计算软件V1.0》(No: 2016SR006487)中完成。为了选取模型最佳自变量参数, 本文将变量投影重要性准则应用到其中, VIP值代表自变量对模型拟合的程度, 自变量对因变量的解释作用相同, 则自变量的VIP均接近于1。Wold^[13]建议VIP值小于1的自变量对因变量的贡献较小, 可以考虑剔除。因此, 本研究自变量VIP的阈值为大于等于1。基于PLSR建模方法的优点, VIP技术主要用于样本较少且几个自变量间相关性较强的情形, 在一定程度上弥补了传统线性回归的不足^[14]。

考虑到建模数据和验证数据需要能够充分反映本研究区土壤电导率的实际现状, 对57个土壤样本按照电导率降序排列, 等间隔抽取38个数据为建模集和19个数据为验证集, 分别用于模型的建立与精度验证。将通过显著性检验(P＜0.01)及VIP(≥1.0)的优化光谱指数(NPDIs)作为自变量, 利用PLSR法建立土壤电导率的高光谱估算模型。通过对比各模型的精度参数筛选出最优模型用以研究区土壤电导率的反演, 精度参数包括建模决定系数(R_Cal²)、建模均方根误差(RMSE_Cal)、预测决定系数(R_Pre²)、预测均方根误差(RMSE_Pre)、相对分析误差(RPD)、最小信息准则(AIC)。R_Cal²与R_Pre²判定模型的稳定性, 决定系数越接近于1时, 表明模型的稳定性越好；RMSE_Cal与RMSE_Pre用以表征模型的准确度, 均方根误差越接近于0说明模型的精度越高。RPD的阈值分为3类, 当(a < 1.4)时, 模型估算能力较低；在(1.4≤b < 2.0)范围时, 模型估算能力一般, 当(c≥2.0)时, 模型的定量估算能力表现较好。除此, AIC作为一种评价模型效率的指标, 可以权衡所估计模型的复杂度和此模型拟合数据的优良性, AIC的值越小, 表明模型可以用最少的自由参数最好地解释数据, 较好地避免过度拟合的情况^[15]。

3 分析与讨论 3.1 土壤电导率统计分析

由表 1可知, 建模集和验证集对应的土壤电导率最大值分别为55.70mS/cm和48.30mS/cm, 最小值分别为0.20mS/cm和0.07mS/cm, 均值分别为11.27mS/cm和11.18mS/cm, 变异系数分别为123.64%和117.50%；研究区所有采样点土壤电导率平均值为11.32mS/cm, 变异系数为122.70%, 介于建模集与验证集之间, 数据离散程度属于强变异性系数(C.V＞100%)。

3.2 土壤电导率与含盐量的相关性分析

表 2是57个土壤样品的电导率和含盐量数据, 从中可以看出, 含盐量的变异系数也大于100%, 表现为较强空间变异性。土壤电导率与含盐量的决定系数达到0.99(图 2), 表现为极显著相关关系, 因此在本研究中土壤电导率可以很好地指示土壤含盐量。

3.3 土壤高光谱曲线特征分析

基于上述分析, 本文利用盐渍土的分级标准^[16]对土壤样品高光谱反射率进行归类、求平均, 绘制不同电导率的土壤高光谱曲线(图 3), 以此大致分析土壤高光谱反射率对土壤电导率的响应。本研究区不同盐渍土高光谱反射率曲线的波动基本一致, 土壤样品的高光谱反射率与其电导率未呈现明显的正负相关关系。从图 3a中可知, 在可见光及近红外波段(400—1500 nm)范围内高光谱曲线呈平缓上升态势, 并在1400 nm、1950 nm、2350 nm左右处有3个比较明显的水分吸收谱段。总体而言, 重度盐渍土高光谱反射率偏低, 非盐渍土高光谱反射率偏高；除600—700 nm波段范围, 不同土壤电导率的高光谱反射率曲线较易区分, 且在1200—1600 nm波段区间, 4种盐渍土的高光谱曲线差异性最大(图 3)。

3.4 优化光谱指数与土壤电导率相关性分析

在波段400—2400 nm范围内, 通过两波段优化算法计算土壤样品的原始高光谱反射率R及其对应的5种光谱变换下所有可能两波段组合的优化光谱指数(NPDIs), 并与土壤电导率进行Pearson相关性分析, 在Matlab2012a软件中制作的相关性图如图 4所示, 其中分数阶微分方法以0.2为间隔进行0—2阶的10种光谱预处理, 经分析(|r| ≥0.870, ρ < 0.01)发现1.6阶时的相关性最高, |r|为0.888且自变量数量适中, 因此将1.6阶微分选定为本文的最优分数阶微分光谱变换方法。为了分析高光谱数据的5种数学变换对优化光谱指数的影响, 统计了不同数学变换下NPDIs与土壤电导率相关系数在(|r| ≥0.870, ρ < 0.01)范围内的数量, 统计结果见表 3。由图 4及表 3可知, 原数据和对数倒数、平方根、1.6阶微分变换后得出的优化光谱指数与土壤电导率呈负相关关系, 而对数、倒数变换为正相关关系；原数据及其对应的平方根、对数倒数、倒数和1.6阶微分变换后优化光谱指数与土壤电导率的相关性都较好, 相关性绝对值均大于0.8(图 4中红色部分), 且敏感指数(NPDI, |r| ≥0.870)的极值大小对应变换顺序为, 因此, 数学变换预处理对于提高优化光谱指数与土壤电导率的相关性具有一定的作用；对数和一阶微分变换却降低了相关性, 相关性绝对值最高仅达0.65与0.76, 这说明不是所有的光谱变换方式都有利于提高优化光谱指数与土壤电导率的敏感程度。从敏感波段组合的数量及范围方面来看, 以3种变换得出的NPDIs与土壤电导率的相关性绝对值超过0.870的个数呈增加趋势, 其中, 这种效果在对数倒数变换形式下最为明显, 而分数阶微分变换表现出相反的态势, 可以说不同的数学变换预处理在敏感波段组合的数量方面有不同的作用；不论原数据还是数学变换后的数据, 对土壤电导率最敏感的优化光谱指数均是在2040 nm和1880 nm中红外波段附近两波段组合得出的, 并未在不同程度盐渍土高光谱曲线差异性最大的波段区间, 这与李相等^[17]得出的土壤电导率高光谱特征波段结果相符。

3.5 土壤电导率PLSR估算模型建立及精度分析

偏最小二乘回归法是当前应用最为广泛的高光谱建模方法, 为PLSR法更深入地分析数据, 在建立研究区土壤电导率PLSR估算模型之前, 本文利用VIP技术进一步对自变量进行筛选, 图 5展示了VIP准则筛选过程与对应波段组合。由图 5可知, 基于原始高光谱反射率及其5种变换优化的NPDI在自变量选择的情况上基本一致, 筛选出的自变量数分别为17、17、28、21、23、15个, 较原有的敏感波段组合数, VIP准则筛选效果明显, 剔除了对模型贡献小及相对累赘的参数。

根据PLSR预测模型精度参数(表 4), 对建立的6种模型进行分析发现(另外6种没有采用VIP技术的模型在此仅作为对比, 不作深入分析), RPD大于2的模型有5个, 分别是基于原数据、对数倒数、倒数、平方根、1.6阶微分变换下优化的NPDI而建立的估算模型, 表明这些预测模型具有一定的定量反演能力, 模型精度R_pre²分别为0.79、0.79、0.82、0.80、0.84, RMSE_Pre分别为3.11、2.89、2.33、2.56、2.07 mS/cm, AIC分别为219.12、217.33、180.05、196.65、158.11。其中, 基于1.6阶微分变换预处理高光谱反射率并对光谱指数(NPDI)进行波段优化运算后, 通过VIP技术筛选的优化指数(x1=NPDI(2012, 1895 nm)、x2=NPDI(2020, 1893 nm)、x3=NPDI(2016, 1893 nm)、x4=NPDI(2020, 1892 nm)、x5=NPDI(2017, 1893 nm)、x6=NPDI(2019, 1893 nm)、x7=NPDI(2083, 1881 nm)、x8=NPDI(2078, 1882 nm)、x9=NPDI(2078, 1881 nm)、x10=NPDI(2083, 1882 nm)、x11=NPDI(2014, 1895 nm)、x12=NPDI(2084, 1881 nm)、x13=NPDI(2015, 1895 nm)、x14=NPDI(2018, 1893 nm)、x15=NPDI(2023, 1892 nm))为自变量建立的PLSR预测模型为最佳, 模型公式为：y=997.521-29.839×x1-34.239×x2-29.830×x3-31.939×x4-31.131×x5-32.134×x6-34.796×x7-35.302×x8-34.280×x9-35.723×x10-31.258×x11-34.454×x12-31.333×x13-30.969×x14-33.990×x15, 预测精度为R_pre²=0.84, RMSE_Pre=2.07mS/cm, RPD=2.94, AIC=158.11。而基于对数预处理高光谱数据优化的NPDI对土壤电导率预测能力较低, 预测精度R_pre²为0.67, RMSE_Pre为4.48 mS/cm, RPD为1.25, AIC为249.16。将1.6阶微分预测模型应用到全部采样点中, 得到对应土壤采样点的电导率预测值, 通过反距离权重插值法绘制研究区内土壤电导率的实测与预测图, 如图 6所示。预测值的插值平均绝对误差、均方根误差分别为0.98mS/cm、4.96mS/cm, 实测值插值图的为1.05、4.53mS/cm, 二者极值出现的地区以及整体的空间分布基本一致且差异较小。这说明本文使用的方法在土壤电导率的监测中具有一定的应用价值, 为快速获取盐渍土壤电导率提供了可能的途径。

3.6 讨论

盐渍土的形成与水盐运移有密切联系, 土壤水盐运移模型可以模拟区域土壤水盐运移过程, 预报土壤水盐动态变化, 对于改造和利用盐渍土具有重要的作用^[18]。表层土壤水盐是土壤水盐运移模型的重要边界条件参数, 准确的表层土壤水盐信息可以提高水盐运移模型的模拟与预测精度。通过高光谱遥感监测土壤电导率便可以及时高效地掌握表层土壤盐分的状况, 这是因为对于自然土壤而言, 尽管电导率和含盐量都可以反映出土壤盐渍化的程度, 但相关研究^{[17, 19]}表明土壤高光谱信息对土壤电导率的响应较含盐量敏感, 以土壤电导率替代含盐量进行土壤盐渍化高光谱估算研究是一种精度更高、速度更快的方法。

传统的高光谱处理方法对土壤电导率高光谱建模时, 敏感波段通常是在一维层面上以土壤电导率与高光谱反射率的相关性分析来确定, 相关性越高, 波段的敏感程度越高。对高光谱反射率的预处理目的就是提高土壤电导率与反射率之间的相关性, 并进一步提高预测模型的精度^[20-22]。高光谱丰富的波段信息为两波段优化算法提供了更多的可能组合, 海量光谱数据的两波段优化算法能充分提取与土壤电导率相关性最大的波段组合, 在复杂的高光谱参数中达到快速寻优的效果, 深度挖掘高光谱数据从而进一步提高土壤电导率的高光谱估算精度, 减少环境因素等对建模的影响^[23-28]。于是本文想探讨的主要问题就是对高光谱反射率的预处理是否同样能够有利于优化光谱指数更好地估算土壤电导率。

分数阶微分在阶数上对整数阶微分的概念进行了扩展, 相比整数阶微分, 具有记忆性、遗传性以及非局部性, 在系统控制与诊断、数字滤波、信号与图像处理等领域有着较为广泛的应用。在光谱分析领域, 近期的相关研究指出, 对于高光谱这类具有海量信息的高维数据源, 分数阶微分也能够很好地挖掘潜在信息, 弥补整数阶微分可能造成某些信息丢失的不足, 并极大地扩充光谱数据预处理的方法, 为高光谱研究提供一个全新的角度^[29-32]。但是, 分数阶微分在单波段高光谱预处理中存在的一个问题是：它虽然可以较好地增加敏感波段的数量, 却不能有效地提高相关性。因此, 本文将分数阶微分与传统光谱变换方法进行对比研究, 探究分数阶微分在优化光谱指数中的效用, 发现分数阶微分预处理方法可以有效地提高优化光谱指数相关性的极值。

为比较一维与二维高光谱数据处理效果的差异, 绘制了高光谱反射率及其不同数学变换后在一维层面上与土壤电导率进行相关性分析得出的最大相关系数图。由图 7可知, 传统光谱变换的相关性提升效果不明显, 平方根变换的相关性最高为0.633, 对数倒数、倒数变换的相关性明显降低, 最低为-0.566, 而在优化光谱指数中表现较差的对数变换没有太大的降低；在分数阶微分变换中, 多数都有提升效果, 最高的为一阶微分变换, 达到-0.673, 这与优化光谱指数的结果形成了鲜明的差异, 最低的为二阶微分, 仅达-0.436, 而在二维层面表现最好的1.6阶微分在一维中没有起到相同的作用；从相关性质方面来看, 原数据、对数和平方根变换呈正相关, 对数倒数、倒数变换呈负相关, 这与优化光谱指数的结果相异, 分数阶微分的相关性质以一阶微分为界, 0.2-0.8阶为正相关, 1.0—2.0阶为负相关。综上所述, 不管在一维层面上还是二维层面上, 合适的高光谱数据预处理方法都在一定程度上对提高相关性有所帮助, 分数阶微分总体上优于传统预处理方法, 而且优化光谱指数对比于传统的高光谱处理方法来说具有明显的优势^[33-34]。

与已有研究^[35-42]相比, 本文的研究特色在于：将高光谱数据分析从传统的一维层面上升至二维层面, 并结合分数阶微分预处理进行光谱指数波段优选, 用于建立土壤电导率估算模型, 以提高土壤电导率反演精度, 为土壤盐渍化相关研究提供一种新的思路和方法。而且本文得到的优化光谱参数可为快速准确寻求卫星传感器中监测干旱、半干旱地区土壤电导率的最佳波段提供依据, 此外, 波段的优化也可以为设计特定波段的主动传感器提供理论基础, 进一步减少高光谱海量数据处理的工作量, 为实现土壤盐分信息的高效监测服务。最后, 研究区虽然属于典型的干旱、半干旱区, 但是干旱、半干旱区乃至中国具有区域异质性^[43], 这就会不可避免地导致本研究确定的土壤电导率反演优化光谱参数仍有一定的地域局限性, 因此, 比较光谱指数优化算法在不同地区的最佳参数并找出普适性高光谱参数, 将是值得研究的方向。

4 结论

本研究以新疆艾比湖流域为研究靶区, 通过土壤高光谱数据对简化光谱指数(NPDI)进行波段优化运算, 分析二维层面上基于高光谱数据不同变换形式的优化光谱指数与土壤电导率的相关性, 通过相关系数显著性检验(P＜0.01, |r| ≥0.870, VIP≥1.0)筛选敏感高光谱参数, 并利用偏最小二乘回归分析法建立土壤电导率高光谱估算模型, 得出以下结论：

(1) 优化光谱指数与土壤电导率之间的相关性显著, 通过运用高光谱数据的5种数学变换形式发现, R、1/logR、1/R和FOD形式下的优化光谱指数均与土壤电导率呈现很好的相关性, 相关性最高达到±0.888。

(2) 原数据及其不同数学变换后发现, 主要位于2040 nm和1880 nm左右波段范围组合的简化光谱指数(NPDI)与土壤电导率之间相关性显著提高。

(3) 对土壤电导率高光谱PLSR建模发现, 基于1.6阶微分(FOD)变换预处理高光谱反射率并对简化光谱指数(NPDI)进行波段优化运算后, 通过VIP技术筛选的优化指数为自变量建立的PLSR预测模型为最佳, 预测精度为R_pre²=0.84, RMSE_Pre=2.07mS/cm, RPD=2.94, AIC=158.11。这为研究区土壤盐渍化动态监测提供了更精确的基础理论与方法。