土壤中不同盐分离子对作物的毒害不同，土壤盐分离子的分异对作物的危害性往往高于总盐的危害^[1-2]。传统土壤盐分离子的测定需要耗费时间，且步骤繁琐，而高光谱遥感技术可以充分挖掘光谱信息和构建高精度的模型^[3]。土壤盐分离子的含量和特征在一定程度上决定了盐渍化的性状和程度，将土壤盐分离子和高光谱数据综合考虑，筛选其特征波段，可以明确研究区土壤盐分离子的特征光谱，提高土壤盐分离子的估算精度和实用性^[4-5]。郭鹏等^[6]采用混合波段和筛选波段的方法以提高土壤盐分离子估算模型的精度；张俊华等^[7]基于0 ~ 5 cm和0 ~ 20 cm对比分析土壤盐分离子在不同土层的估算模型。针对分数阶微分预处理高光谱数据的方法中，有采用Savitzky-Golay加一阶微分变换和以不同间隔的分数阶微分变换等不同方式^[8-9]。亚森江·喀哈尔等^[10]认为利用分数阶微分优化光谱指数对土壤电导率的特征波段敏感程度更强。高光谱数据常用的模型估算方法有神经网络模型、偏最小二乘、支持向量机等^[11-13]。袁玉芸等^[14]利用地理加权回归模型(GWR)分析于田绿洲土壤表层盐分高光谱数据估算精度优于OLS模型。GWR能够考虑到空间异质性，将采样点的坐标代入不同的系数，就能进行空间坐标分析，通过建立空间尺度上每个点的局部回归方程，能有效揭示土壤属性空间变异特征关系^[15]。空间关系约束的高光谱土壤盐分离子估算方面研究鲜有报道。针对干旱区湖滨绿洲土壤盐分离子的空间异质性，开展以博斯腾湖湖滨绿洲为研究区，分析土壤盐分离子与相关光谱参量间的关系，构建地理加权回归模型，将每个样本点加入空间坐标估算土壤盐分离子含量。

1 材料与方法 1.1 研究区概况

博斯腾湖湖滨绿洲行政区划隶属于新疆博湖县，位于新疆焉耆盆地东南部，地理位置介于86°15′ ~ 86°55′E，41°45′ ~ 42°10′N(图 1)。其是典型的人工绿洲和自然绿洲混合的山前湖泊绿洲，面积为1 360 km²。湖滨绿洲受到博斯腾湖水域的调节，年平均气温8.2 ~ 11.5 ℃，无霜期175.8 ~ 211.3 d，年蒸发量1 880.0 ~ 2 785.8 mm，年降水量47.7 ~ 68.1 mm，蒸降比高达40:1。研究区内土壤类型主要有沼泽地、草甸土、灌耕潮土、风沙土、盐土、棕漠土等，地下水埋深1.0 ~ 2.5 m，矿化度为0.1 ~ 10 g/L^[16-17]。研究区土地利用/覆盖类型主要为耕地与未利用地，土壤盐分以氯化钾钠盐和硫酸钾钠盐为主，主要土壤盐分离子有HCO₃^–、Cl^–、SO₄^2–、Ca²⁺、Mg²⁺、Na⁺+K⁺。

1.2 土壤样品采集及测定

根据研究区土壤类型的现状，结合研究区的土壤盐分状况，采用GPS定位技术，使采样点尽可能遍及研究区范围内的主要土地利用类型，考虑采样点的水文地貌条件、植被覆盖类型以及土壤盐渍化程度等因素，选取34个样点进行土壤采集(图 1)，每个样点3次重复。采样时间为2019年7月8日，采用分层随机取样的方法，每个样本点范围分0 ~ 10、10 ~ 20 cm采样土层进行，共计68份。土壤光谱测定采用ASD FieldSpec3便携式地物光谱仪，波段范围为350 ~ 2 500 nm，野外测量需要选择晴朗无风的天气，采集光谱前对光谱仪白板校正，且校正间隔控制在5 min^[18]，每个土壤样品重复测量20条光谱曲线，取其平均值作为此样本的光谱数据。对均值处理后的光谱曲线去除水分吸收带波段，对去除干扰波段的光谱曲线用Savitzky-Golay滤波方法进行平滑处理。

1.3 数据处理与研究方法 1.3.1 分数阶微分

分数阶微分是将整数阶微分的阶数扩展至任意阶。目前，在众多的分数阶微分定义形式中^[19-20]，Grünwald-Letnikov(G-L)定义的分数阶微分较为常用，主要是应用一元函数差分来实现，表达式为：

$ \begin{gathered} \frac{{{d^V}}}{{d{\lambda ^V}}}f\left( \lambda \right) = f\left( \lambda \right) + \left( { - V} \right)f\left( {\lambda - 1} \right) + \frac{{\left( { - V} \right)\left( { - V - 1} \right)}}{2} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;f\left( {\lambda - 2} \right) + \cdot \cdot \cdot + \frac{{\Gamma \left( { - V + 1} \right)}}{{n!\Gamma \left( { - V + n + 1} \right)}}f\left( {\lambda - n} \right) \\ \end{gathered} $

(1)

式中：V为阶数，$\Gamma $为Gamma函数，$\lambda $为对应的值，n为微分上下限之差。

1.3.2 光谱指数

根据研究区土壤盐分离子状况，结合前人的研究结果，光谱指数采用比值盐分指数(ratio salt indices，RSI)、差值光谱指数(difference spectral indices，DSI)和归一化光谱指数(normalized difference spectral indices，NDSI)^[21]。根据全光谱的第i波段和第j波段的所有组合计算二维相关系数，绘制形成光谱系数矩阵图。光谱指数计算公式：

$ {\rm{SI}} = {R_i}/{R_j} $

(2)

$ {\rm{DSI}} = {R_i} - R $

(3)

$ {\rm{NDSI}} = ({R_i} - {R_j})/({R_i} + {R_j}) $

(4)

式中：R_i和R_j是土壤盐分离子在第i波段和第j波段的土壤光谱反射率。

1.3.3 地理加权回归

Fotheringham等^[22]提出地理加权回归模型(geographically weighted regression，GWR)能够将空间坐标代入到回归参数中，采用的是局部光滑思想，这样能使局部加权的最小二乘法进行样点的估测，其中权重是回归点所在位置到其他观测点的位置之间距离的函数。该模型表达式^[23]为：

$ {y_i} = {\beta _{\rm{0}}}({u_i}, {v_i}) + \sum\limits_{i = 1}^k {{\beta _i}({u_i}, {v_i}){x_{ik}} + {\varepsilon _i}} \;i= 1, 2, 3, ···, \\ n $

(5)

式中：y_i为样点i的因变量；β₀为回归常数项；(u_i, v_i)为第i个采样点的坐标(如经纬度)；β_i(u_i, v_i)是i个采样点上的第k个回归参数，是地理位置的函数；ε_i为误差项；x_ik为第i个样点上的k个变量的观测值。采用带宽优化方法中易于与经典的线性回归模型进行比较的AIC信息准则，表达式为：

$ {\rm{AIC}} = 2nln\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\delta } + nln(2\pi ) + n\left[ {\frac{{n + tr(s)}}{{n - 2 - tr(s)}}} \right] $

(6)

式中：tr(s)是带宽b的函数；$\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\delta } $是随机误差项方差的极大似然估计，$\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\delta } $=RSS/n–tr(s)；当AIC值最小时，所选择的权函数的带宽最优^[24]。在GWR4.0软件支持下完成GWR模型回归过程，模型精度评价采用决定系数R²和均方根误差(RMSE)进行评价。

2 结果 2.1 主要土壤盐分离子含量描述性统计

由表 1可知，主要土壤盐分离子为HCO₃^–、Cl^–、SO₄^2–、Ca²⁺、Mg²⁺、Na⁺+K⁺。从均值看，研究区土壤的SO₄^2–、Na⁺+K⁺含量最高，分别为1.58和1.59 g/kg；阴离子中Cl^–次之；HCO₃^–含量最少，为0.11 g/kg。Ca²⁺、Mg²⁺含量分别占阳离子的10.64% 和10.95%，Cl^-占阴离子的15.06%。结合前人的研究结果，在后续建模中以3:1比例选取建模集与验证集^[25]。

2.2 一维光谱相关性特征波段选择

由图 2可知，对原始光谱反射率以0.2为间隔，在0 ~ 2阶进行微分处理，将分数阶微分变换下的光谱反射率与HCO₃^–、SO₄^2–、Ca²⁺、Mg²⁺、Na⁺+K⁺、Cl^–进行相关性分析，筛选出通过极显著性检验且相关性较高的阶数进行分析。1.6阶和2阶与HCO₃^–的相关性较小，通过极显著性检验的波段较少；1阶与Cl^–相关性在分数阶微分变化下最优，最大相关系数绝对值为0.57；SO₄^2–、Ca²⁺、Mg²⁺、Na⁺+K⁺、Cl^–的较大相关系数集中在正相关。波段筛选原则为：在样本数为68，通过显著性检验，即选取相关系数大于0.315的最优4条波段。

由表 2可知，对优选出的土壤盐分离子建模所需要的特征波段进行分析，HCO₃^–的特征波段主要集中在419 ~ 490 nm可见光紫、蓝光波段；Na⁺+K⁺的特征波段主要集中在468 ~ 874 nm、1 182 ~ 1 539 nm可见光绿、黄、橙、红光及近红外波段；Cl^–的特征波段主要集中在1 253 ~ 2 360 nm近红外波段；SO₄^2–的特征波段主要集中在468 ~ 830 nm可见光绿、黄、橙、红光波段，与近红外波段也有一定相关性；Ca²⁺的特征波段主要集中在1 080 ~ 2 316 nm近红外波段，Mg²⁺的特征波段主要集中在912 ~ 2 175 nm近红外波段。1 ~ 2阶的光谱反射率与土壤盐分离子的相关性较高，达到显著水平，从而筛选出较优的特征波段。分数阶微分变换下的光谱数据更能深度挖掘光谱中的潜在信息。

2.3 二维土壤指数特征波段选择

土壤盐分离子含量和光谱指数之间的相关性，能够利用二维光谱矩阵系数图进行可视化表达。通过建立HCO₃^–、SO₄^2–、Ca²⁺、Mg²⁺、Na⁺+K⁺、Cl^–分别与原光谱反射率的RSI、DSI、NDSI光谱矩阵系数图，从中选取相关性较高的光谱矩阵图进行分析，如图 3可知：同种离子在RSI、DSI、NDSI的光谱指数下二维光谱矩阵系数图具有相似性。其中Cl^–的相关性整体大于其他离子，特征波段主要集中在1 461 ~ 1 597 nm和1 710 ~ 1 749 nm，主要位于近红外波段；Mg²⁺在RSI和DSI的相关性与Cl^–的光谱矩阵图相似，但相关系数值比Cl^–低；HCO₃^–特征波段主要集中在1 550~ 1 760 nm；Na⁺+K⁺特征波段主要集中在1 742 ~ 2 395 nm；SO₄^2–、Ca²⁺、Mg²⁺的特征波段主要集中在1 977 ~ 2 340 nm的近红外波段，其中SO₄^2–的特征波段还集中在600 ~ 700 nm的可见光红光波段。

优选出的特征波段，用于构建模型的二维土壤指数见表 3，由于Mg²⁺、Ca²⁺在RSI、NDSI光谱指数下筛选的特征波段相同，分别为(1 989, 2 066)、(2 005, 1 982)，为避免冗余，选择RSI建模。

2.4 地理加权回归模型分析

一维相关优选特征波段和二维相关优选特征波段与土壤盐分离子都具有较高的相关性，为寻找较优的估算模型，分别以HCO₃^–、SO₄^2–、Ca²⁺、Mg²⁺、Na⁺+K⁺、Cl^–为因变量，优选的特征波段分别为自变量，构建地理加权回归模型。由表 4可知，由于一维相关和二维相关的RMSE相差较大，因此从一维和二维相关中各筛选出一个较优的建模集，再根据验证集的较优模型确定土壤盐分离子的最优模型。Cl^–估算模型中，1阶的R² = 0.74最大，RMSE = 0.03最小；DSI的R² = 0.82最大，RMSE = 0.11最小，均为建模集较优模型。SO₄^2–估算模型中，1.2阶的R² = 0.91最大，RMSE = 0.02最小；DSI的R² = 0.89最大，RMSE = 0.36最小，均为建模集较优模型。Na⁺+K⁺估算模型中，1.2阶和2阶的R² = 0.74最大，RMSE = 0.01最小；RSI的R² = 0.94最大，RMSE = 0.22最小，均为建模集较优模型。HCO₃^–估算模型中，2阶的R² = 0.23最大，RMSE = 0.01最小；NDSI的R² = 0.17最大，RMSE = 0.09最小，均为建模集较优模型。Ca²⁺估算模型中，1阶的R² = 0.75最大，RMSE = 0.03最小；DSI的R² = 0.74最大，RMSE = 0.18最小，均为建模集较优模型。Mg²⁺估算模型中，1.2阶的R² = 0.69最大，RMSE = 0.01最小；RSI的R² = 0.83最大，RMSE = 0.10最小，均为建模集较优模型。

由图 4可知估算模型的验证集，对于Cl^–而言，基于1阶的R² = 0.93，RMSE = 0.11；SO₄^2–基于DSI的R² = 0.77，RMSE = 0.32；Na⁺+K⁺基于RSI和DSI的R² = 0.74，RMSE = 0.19；HCO₃^–基于DSI的R² = 0.83，RMSE=0.05；Ca²⁺基于DSI的R² = 0.96，RMSE=0.08；Mg²⁺基于RSI的R² = 0.93，RMSE = 0.08，均为验证集较优模型。综合建模集和验证集来看，Cl^–的最优模型是1阶，其中1阶的AIC值为0.50；SO₄^2–的最优模型是1.2阶，1.2阶的AIC较小为122.87；Na⁺+K⁺的最优模型是RSI，AIC最小的是DSI，RSI比DSI的AIC高3.08；Ca²⁺的最优模型是DSI，其中1阶微分的AIC较小为1.79，DSI在光谱指数中AIC为14.61；Mg²⁺的最优模型是RSI，RSI的AIC为11.31，比最小的AIC高5.46。从整体来看，选择的最优模型，多是分数阶微分或光谱指数中AIC值较小的模型。Na⁺+K⁺和Mg²⁺在二维光谱指数构建的估算模型，建模集和验证集的决定系数均优于一维。

3 讨论

通过优选高光谱的特征波段，对优选的1种特征波段和优选的2种特征波段进行比较，优选的特征波段具有一定重合性，Na⁺+K⁺特征波段在光谱变换下集中在可见光绿、黄、橙、红光及近红外波段，2种波段选取中重合部分为近红外波段；SO₄^2–的特征波段在光谱变换下主要集中在可见光波段，2种波段选取重合部分也为可见光红光和近红外波段；Cl^–、Ca²⁺、Mg²⁺的特征波段在光谱变化下和2种波段选取中都集中在近红外波段；HCO₃^–在光谱变换及指数选取的特征波段分别集中在可见光紫光、蓝光波段和近红外波段。基于GWR模型对于研究区土壤盐分离子估算具有较优的结果。GWR模型对Na⁺+K⁺估算效果较好，这与宋梦洁和李新国^[26]在相同研究区下利用多元线性回归模型估算的分析结果相同。GWR模型能够对Na⁺+K⁺进行较好的预测，依据梁东等^[27]的研究结果，研究区的土壤盐分与Na⁺+K⁺呈极显著正相关，且Na⁺+K⁺在不同植被类型中含量都是最高。土壤盐分离子含量较高，可进行很好的预测，这或许也说明光谱反射特征机理是对主要离子组成的响应。土壤盐分离子含量、组成、结合方式的差异，在提取特征波段和光谱构建模型时难度较大。GWR模型提高土壤盐分离子的估算效果受到土壤盐分离子与各变量间相关关系的空间非稳定性程度^[28]。本研究的土壤盐分离子估算效果与代希君等^[1]在新疆南疆的研究结果基本相同，SO₄^2–、Ca²⁺、Mg²⁺、Cl^–的估算精度较高，对于土壤盐分离子中含量较低的HCO₃^–估算精度较低。

4 结论

1) Na⁺+K⁺的微分变换光谱特征波段主要集中在468 ~ 724 nm与1 182 ~ 1 539 nm，相关性值最大为0.48；二维土壤指数特征波段主要集中在1 742 ~ 2 395 nm；基于RSI的特征波段优选下地理加权回归模型对Na⁺+K⁺含量的估算效果较好，建模集R² = 0.94，RMSE = 0.22，验证集R² = 0.74，RMSE = 0.19。

2) SO₄^2–的微分变换特征波段主要集中在468 ~ 724 nm，其中在1阶、1.2阶、DSI的建模集和验证集精度R²值最大，R²值均大于0.71；Ca²⁺、Mg²⁺、Cl^–的特征波段均在近红外区域，且验证集精度高于建模集精度，其中Ca²⁺在光谱指数DSI的验证精度R²值最大为0.96，RMSE为0.08。

3) 基于地理加权回归模型的土壤盐分离子估算的最优模型R²值范围在0.74 ~ 0.96，最优模型的AIC值范围在0.50 ~ 122.87；土壤离子中含量相对较高的Na⁺+K⁺、SO₄^2–、Cl^–的地理加权回归模型精度高于含量较低的Ca²⁺、Mg²⁺。