研究表明, 在干旱区由于盐渍化和次生盐渍化所造成的大面积土地被迫撂荒和植被的第一性生产力下降, 严重降低了土地生产力, 限制着农业生产的发展^[1]。同时由于为短期利益的农业活动而忽视了长期的后果是加重盐渍化的重要原因^[2]。现在传统的方法来监测盐渍化表现出多方面的不足, 而融合了遥感RS、地理信息系统GIS的技术方法却显示出很多优势^[3]。国外学者利用遥感提取盐渍土信息主要有以下方法：波段选择、主成分分析、K-T变换、IHS变换、图像比值运算, 以及最大相似性分类、神经网络、决策树、表面特征分解、模糊分类以及多源数据集成建模等信息提取技术^[4]。国内学者对土壤盐渍化状况遥感监测的研究多从静态角度, 盐渍化信息提取的方法主要有目视解译、基于影像光谱特征的数字图像处理技术、遥感信息地理信息综合分类、基于地学知识发现的识别与分类、人工神经网络分类器 (RBF或BP神经网络模型) 等^[4-5]。国内的一些学者还运用主要由各种植被指数、盐分指数或碱性指数等来构建特征空间来提取盐渍化信息^[6-9]。很多研究中都用到了归一化差值植被指数 (NDVI) 来作为参数, 但NDVI值与土壤电导率之间存在的相关性可用于区分盐渍土和非盐渍土, 但无法区分盐土和碱土, 同时由于不同类型的植被生长在不同的盐渍化水平下, NDVI表现出很大的不确切性^{[2, 4]}。因此单纯的考虑一种表征植被状态的指数并不可取。综合国内外已有的研究可以发现存在以下一些可以改进的地方：

(1) 多数对于监测盐渍化选择的参数是基于光谱特征^[10]、母质成因参数^[11-12]、地貌特征和数字高程模型 (DEM)^[13]。但是对于植被的信息大多是用如NDVI、修改型土壤调整植被指数 (MSAVI) 和转换型土壤调整植被指数 (TSAVI) 等相关的指数。这些指数虽然各有优势, 但解决的也只能是较理想下的土壤背景和植被覆盖等两大地类信息, 不能用在每一种自然环境状态下^[2]。同时, 很多研究只是以盐生植被和非盐生植被来介绍植被在不同盐分下的长势情况, 但较少对盐生植被与非盐生植被进行区别研究 (即划分出盐生植被和非盐生植被, 然后选择不同的方法或设置不同的参数进行监测, 最终把不同监测的结果进行汇总, 得到整个研究区盐渍化空间分布信息), 这很可能使得最后提取的盐渍化信息精度有所降低。根据相关研究^[14-15], 这种植被类型的划分更能有效的表征盐渍化程度的分布情况, 因此本文也通过探讨划分不同的植被类型进行分区来反演盐渍化信息。

(2) 很多研究对于盐渍化信息的提取只是针对某一特定时期进行的, 选择的影像数据及其他辅助数据都围绕着采样时间进行, 其结果是影像等辅助数据选择只是一期的的结果, 但盐渍化在一定的区域内是多种自然原因和人类原因长期共同作用的结果^[5], 而MODIS长时间序列的植被指数数据满足长期观测的条件。

(3) 不同的植被与土壤类型, NDVI等相关植被指数、盐分指数的变化^[16]、管理措施、虫害、疾病、气候和土壤属性^[2]等会使得这些参数与盐分信息呈现出非线性的变化, 因此有必要探讨运用人工神经网络强大而复杂的非线性拟合能力来反演盐渍化信息的时空分布情况。

(4) 反演出的植被指数与含盐量的相关性很低, 是因为宽区间的光谱波段不能捕捉由于盐分的变化而出现的生物物理化学上的差异^[17], 因此需要寻求一种能平衡时间分辨率、空间分辨率和波谱分辨率的数据源, 而MODIS数据同时具备这些优点^[15], 是一种良好的遥感数据源。

(5) 研究表明当土壤足够干燥时, 能够用遥感数据来进行较高精度的反演^[14], 因此这限制了遥感数据主要适用于区域中降水稀少的干季, 而不适用于降水集中的湿季^[6]。同时地表上一般均有植被覆盖, 特别是在农耕区, 因此监测作物的长势情况是一种潜力巨大, 能够代替传统方式的方法, 因为植被长势差的地区一般表征着重盐渍化程度集中的区域^[14]。这种方法在一些地方已经有成功运用的案例, 特别是在完全没有植被存在的情况下^[18-19]。

(6) 在先前的研究中, MODIS连续数据和盐分的相关性很可能被低估, MODIS连续时间序列的数据并没有被充分利用, 特别是在物候信息方面^[15]。已有的报道中虽有用物候参数来进行研究, 但研究较少, 主要考虑到的只有植被类因素^[14-15], 精度难免有所影响, 极有可能存在一定的不确定性, 因此本文通过加入物候信息来探讨反演盐渍化的信息。

综上, 对于植被监测盐渍化信息的能力还没有充分挖掘。一般情况下植被受到盐渍化的影响会有着较低的光合作用, 同时使得可见光反射率增高, 近红外反射率降低^[20-21]。盐生植被一般具有较高的耐盐性, 因此是区分盐渍化地区和非盐渍化地区的理想指标^[22]。MODIS连续序列的数据平衡了光谱分辨率、时间分辨率和空间分辨率的优势^[15], 可以用来进行大尺度的监测盐渍化信息。因此本文以渭干河-库车河三角洲绿洲为研究区, 基于MOD13A1-NDVI数据来提取2013年植被多种物候信息, 综合物候和非物候影响因子来反演研究区的盐渍化信息。

1 研究区概况

渭干河—库车河三角洲绿洲 (以下简称渭库绿洲) 位于新疆南部的塔里木盆地中北部, 属大陆性暖温带干旱气候^[7]。根据对实地采样区域确定的边界范围为: 82°10′—83°50′E；41°06′—41°40′N (图 1)。属于大陆性温带干旱气候。多年平均气温11.6℃, 多年平均降水量为52mm, 蒸降比为54:1^[23]。土壤主要为潮土和草甸土^[7]。绿洲自然植被以盐生植被为主, 如柽柳 (Tamarix)、盐节木 (Halocnemum strabliaceum)、盐爪爪 (Kalidium foliatum)、花花柴 (Karelinia caspica)、芦苇 (Phragmites communis)、苦豆子 (Sophoraio pecuroides)、白刺 (Nitrarria sibirica)、黑刺 (Lycium ruthenicum)、骆驼刺 (Aihagi sparsifolia)、盐穗木 (Halostachys caspica)、胡杨 (Populuseu phratica) 等肉质多汁泌盐植物和盐生灌丛, 主要呈片状分布在绿洲外围^[24]。农田植被主要由小麦、玉米、棉花等种植作物构成, 主要分布在绿洲内部水资源丰富、种植环境良好、排灌渠设施完善的种植区^[24]。渭库绿洲盐类沉积规模大, 土壤普遍盐渍化, 盐渍化总面积超过研究的50%, 重度盐渍化 (5—10 g/kg) 面积达到30%, 因此绿洲的发展受到严重威胁^{[23, 25]}。

2 研究方法与原理

综合物候参数、植被指数、盐渍化指数、干旱指数、特征空间、影像反射率、地形参数等22种参数来反演研究区盐渍化的空间分布信息。这些参数分别代表了物候特征、植被长势情况、盐渍化分布特征、干旱情况和地形变化特征等。物候参数的定义见表 1^{[15, 26]}, 非物候参数的定义见表 2。

2.1 电导率数据的采集

选择具有代表性的土壤表层单元, 测量单元的位置、数量、土壤质地和盐分状况和植被类型以及土地覆盖情况等进行选取 (图 1)。对每个测量单元进行表层 (0—10cm) 土壤样品的采集。然后带回实验室中经过自然风干, 磨碎, 过0.5mm孔径筛、制作1:5比例的浸提液来测定土壤电导率和含盐量, 采样时间为2013年8月, 剔除不能提取物候参数对应像元的采样点, 最终选出87个采样点来进行分析 (图 1)。

2.2 地形参数的获取

研究表明地形参数中高程、坡度对盐渍化的聚集有重大影响^[13]。地形的高低直接决定了地下水与地表水的流向, 这使得盐分会随着水流流向而发生不同程度的聚集, 地势较高的区域一般盐渍化程度较低, 地势较低的区域盐渍化程度一般较高。Moore等^[34]研究发现土壤盐分和地形湿度指数有很强的相关性。地形湿度指数具有明确的物理意义:具有相同地形湿度指数的不同空间位置对降雨的水文响应相同, 地形湿度指数频率分布相同的不同流域具有水文相似性。因此其频度分布及空间分布在流域的土壤分布及分布式水文模型等研究中具有重要的意义^[33]。同时地形湿度指数表征了地表水文特征的湿度程度, 被用来对盐渍化程度进行分类^[34]。所有这些参数都反演于500m分辨率的DEM数据, 具体计算公式参见表 2。

2.3 地表覆盖数据的获取

获取渭库绿洲2009年500m分辨率的土地利用产品MODIS-MCD12Q1的土地利用数据, 从中国科学院资源环境科学数据中心 (http://www.resdc.cn) 下载2010年土地利用数据。由于中科院的土地利用数据分辨率是1000 m, 因此2013年使用的是2009年分辨较高的土地利用数据。

根据本研究目的和根据MODIS产品的土地覆盖划分的类型和文献^[15]的分类, 分为农田植被 (代表不同作物类型的植被, 主要是非盐生植被)、盐生植被 (主要是稀疏植被, 分布于为绿洲外围和荒漠区) 和其他地类 (代表了混合了农用地和稀疏植被的区域, 主要包括稠密灌丛、稀疏灌丛、稀疏草原、草地、永久草地、自然植被等地类) 等3大类。但由于其他地类中的样点主要位于绿洲荒漠交错带区域, 受人类活动影响较大，同时在这个漫长的交错带上, 其他地类分布的样点显得特别分散, 如果单独进行反演极有可能会使得结果有较大的不确定性, 因此其他地类中反演结果用不分区时的结果替代。

2.4 温度植被干旱指数 (TVDI)

目前, 基于植被指数 (VI) 和地表温度 (Ts) 的二维特征空间获取温度植被干旱指数 (TVDI) 综合了两个参数特有的生理生态意义, 能有效地减小植被覆盖度对干旱监测的影响, 准确性更高, 实用性更强；且特征空间所需数据较少, 获取方便, 计算得到的TVDI的物理意义明确, 应用较为广泛^[35]。因此本文以TVDI来表征地表的土壤水分, 因为土壤水分和盐分之间的关系直接关系到利用遥感光谱信息监测土壤盐渍化的精度^[36]。很多研究表明TVDI对于不同地区的旱情监测或土壤湿度时空差异方面取得了显著的效果^[37]。同时刘立文^[38]等研究发现引入DEM数据对Ts做地形校正, 减少了地形起伏对能量平衡的影响, 所用模型为C校正模型。但是尽管C校正模型易于使用, 但在森林覆盖的山区, 其模型的几何关系与植被向地生长的特性不相符^[39]。为此, 一些研究者提出了适用于森林山区遥感影像地形校正的SCS模型^[39]。由于SCS与Cosine模型一样存在过度校正的情况, 所以依照C模型的建立方法, 经验参数C被引入SCS形成了SCS+C模型, 取得了较好的校正效果^[40]。因此本文采用基于奈曼最优分配法^[41]的SCS+C模型来校正Ts, 最后得到TVDI (表 2)。

2.5 盐渍化监测模型SDI的建立

综合分析归一化差值植被指数 (NDVI)、盐分指数 (SI) 二者之间的关系, 在此基础之上提出NDVI-SI特征空间概念, 并构建土壤盐渍化遥感监测指数模型 (SDI)。研究结果表明土壤表层含盐量与SDI相关性较高, 对于干旱区盐渍地信息提取具有重要意义^[28]。因此根据采样日期选择MODIS-NDVI数据和计算盐分指数, 构造二维特征空间NDVI-SI, 最后得到盐渍化监测模型SDI, 把SDI作为神经网络模型的一个输入变量。

图 2显示的是特征空间SDI模型的示意图, AEB是特征空间轨迹曲线。空间中距离点D (1, 0) 越远, 表示盐渍化程度越严重, 越近表示盐渍化程度越轻^[28]。根据两点之间的距离公式就可以得到从C点到D点的距离L:

(1)

从而建立盐渍化监测模型 (SDI), 表达式为:

(2)

2.6 植被指数与盐分指数的计算

植被是盐渍化程度的一种重要指示剂, 因此选择能够较好地反应土壤盐渍化程度的4种光谱指数 (表 2):归一化差值植被指数 (NDVI)、修改型土壤调整植被指数 (MSAVI)、转换型土壤调整植被指数 (TSAVI) 和盐分指数 (SI)。经过统计, NDVI和MSAVI、TSAVI的相关性较低, 相关系数均在0.2左右；MSAVI和TSAVI具有高度的相关性, 相关系数接近于1；SI与MSAVI和TSAVI的相关性较低, 均为-0.16；SI与NDVI的相关性较高, 相关系数为-0.84；SDI和SI、NDVI相关性较高, 相关系数均大于0.85(与NDVI是负相关), 与TSAVI和MSAVI相关性较低, 相关系数均为-0.21。综上分析, 尽管部分指数之间具有高度的相关性, 但为了精确的捕捉盐渍化信息细节上的变化, 把这5(包括SDI) 个指标均作为反演的输入变量。

2.7 遥感反射率数据

研究表明地表反射率对于监测地表盐渍化的效果显著^[42]。因此选择蓝波段、红波段和近红外波段反射率作为参数来表征盐渍化的信息。经统计可得蓝波段和红波段的相关系数达到了0.98, 其他波段之间的相关系数的绝对值都小于0.2, 但蓝波段和红波段在不同土壤水分含量的盐渍地中的响应并不一致, 同时作为反演的一个因子能更好的捕捉盐渍化变化的细节信息, 因此这3个波段反射率全部用来作为反演输入的变量。

2.8 植被物候参数的计算

尽管MOD13A1-NDVI数据经过了大气和气溶胶的校正, 但是厚重云雾仍然存在, 因此在做平滑处理之前, 需要根据MODIS质量控制产品 (QA) 对连续时间的序列数据赋予不同的权重^{[15, 26, 43]}。因此通过比较常用的Savitzky-Golay、Asymmetric Gaussian和Double logistic model等3种方法的滤波效果, 本文最终选择用Savitzky-Golay滤波法来对影像数据做处理, 得到重建后的NDVI曲线。拟合公式为:

(3)

权重C_j=1/(2n+1), 数据值y_i被位置t_i的平均值所取代, 这是从二次多项式最小二乘法的2n+1个点的滑动平均窗口中获得。同时加入MODIS产品质量控制影像文件, 赋予像元不同的权重后可以有效的降低噪声的影响。滤波结果如图 3所示, 图中蓝线 (具有不平滑特征的曲线) 表示原始的序列曲线, 红线 (较为平滑的曲线) 代表平滑去噪后的曲线。

采用动态阈值法从NDVI数据集中提取渭库绿洲的物候参数。主要是因为与传统的固定阈值法比较, 动态阈值法与每个像元的NDVI季节变化幅度紧密相关, 根据研究条件的不同, 可以动态地确定阈值, 并消除不同土壤背景值和植被类型的影响^[44]。根据相关者学者^[45-47]对中国及中亚不同地区的研究, 本文将提取生长季开始日期 (SOS) 和生长季结束日期 (EOS) 的系数都设置为30%。用儒略历计算法 (1月1日为一年中第1天) 进行物候期统计分析^[48]。

在经过了平滑和参数设置后, 用TIMESAT软件计算季节性参数包括SOS、EOS、生长季长度 (LEN) 和生长季丰度 (AMP) 等10个参数。表 1中阐述了不同指标所代表的生态学意义^{[15, 26]}。

2.9 BP-ANN模型的预测研究

由于光谱反射率和土壤参数间的转换函数是非线性的, 因此人工神经网络 (Artificial neural network, ANN) 已被广泛用于获取土壤参数。目前利用人工神经网络的研究多集中于对土壤质地的土壤特性研究^[49]。同时盐渍土光谱数据的参数反演是将多维空间向低维空间的映射过程, 这种映射关系是复杂的、非线性的。BP (Backpropagation) 神经网络是智能计算技术的重要分支, 具有快速、有效地处理多维、非线性数据的能力。将BP神经网络用于反演盐渍土含盐量能大大提高反演精度, 推动盐渍土遥感的定量化发展^[4]。所以结合前人的研究, 本研究采用3层BP神经网络, 来预测研究区盐渍化空间分布信息, 见图 4。

神经网络是一种黑箱结构, 多层神经网络能以任意精度逼近给定样本的分布和结构。训练过程无需人工操作, 并且自动获得每个神经元的权重^[50]。它具有运算效率高、自学习能力强、适应面宽等优点^[51]。利用BP神经网络主要分为模型校正和模拟, 主要思路为:建立三层神经网络, 包括输入层, 隐藏层和输出层, 利用随机采样的数据对其进行训练, 在利用训练好的网络对整个研究区的数据进行模拟与预测 (图 4)。

输入层将信号输出到隐藏层, 隐藏层第j个神经元所收到的信号为^[50]:

(4)

式中, k为输出信号对应的模拟单元, t为模拟时间, net_j(k, t) 为隐藏层第j个神经元所接收到的信号, w_{i, j}为为连接输入层和隐藏层之间的参数, 即权重值。

隐藏层对该函数的响应为:

(5)

所以输出层第l个神经元所接收到的信号为:

(6)

式中, PV_l(k, t) 为单元n在模拟时间t时刻第l个神经元所接收到的信号; w_{j, l}为隐藏层和输出层之间的参数, 即权重值。

将每个像元中的22个变量作为输入层数据, 输出层为预测的电导率数据, 由于隐藏层的神经元的数目至少为2n/3个^[52](其中n为输入层神经元的数目), 因此本研究中隐藏层的神经元数目为14个。

2.10 验证方法

选用均方根误差 (RMSE)、平均误差 (ME) 和可决系数 (R²)(实测值与预测值之间线性拟合方程的可决系数) 来评价不同方法预测的精度^{[6, 13]}。

(7)

(8)

(9)

式中, P_i代表预测值, O_i代表实测值, O_ave代表实测平均值, P_ave代表预测平均值, n表示验证点的数量。RMSE越小, R²越高, ME绝对值越接近于0, 精度越高, 反之越低。

3 实验结果和分析 3.1 不同参数的提取结果

依照上述描述的步骤计算出的各非物候参数见图 5。从图 5中可以得到以下结论: ① 从NDVI分布图可以得出, 渭库的绿洲内部的植被覆盖程度明显高于绿洲荒漠交错带和荒漠地区。② 从TVDI分布图中可以得出绿洲内部的土壤含水量明显高于绿洲外围, 从绿洲内部—绿洲外围—荒漠区依次递减。③ 从坡度分布图 (SLOPE) 中可以得到研究区坡度和高程的总体呈北高南低、西高东低的趋势, 但总体差异不大, 是明显的平原地区。因此水流的总体趋势是从北往南和从西向东, 研究区东部是盐渍化的“重灾区”^[53]。④ 地形湿度指数空间分布信息显示出TWI较低的地区多位于地势较高的西北部地区, 较高值多位于东部和南部地区。⑤ MSAVI和TSAVI的分布趋势与NDVI基本相同。⑥ 分布在绿洲内部农用地区域的SI数值明显低于荒漠区域。⑦ SDI空间分布结果表明绿洲农用地的SDI值明显低于非绿洲地区, 在绿洲的大部分区域SDI小于0.4, 非绿洲区域SDI集中于大于0.7, 这在一定程度上表明在绿洲的盐渍化程度普遍比绿洲荒漠交错带和荒漠区低。

3.2 物候参数的提取结果

为了更好的挖掘不同物候参数在不同植被类型下提取盐渍化信息的能力, 用野外实测的土壤含盐量数据与这10种物候参数做最佳线性或非线性方程的拟合, 拟合结果见表 3。

表 3中分为了不分植被类型和分植被类型的两种情况。基于研究目的, 只对最重要的代表非盐生植被的农田区实测样点数据和代表盐生植被的绿洲外围与荒漠地区的实测样点数据进行了拟合。可得以下结论:① 在真实的田间测量情况下, 可决系数R²的变化范围为0.06—0.48, 因此不分地类的拟合结果精度较低, 不足以来解释盐渍化信息的细节信息。② 划分地类后的农用地R²有了明显的提高, 大多在0.55以上, 最大值到达了0.92, 但是对于稀疏植被区R²较低, 因此单纯的用某个拟合曲线不能很好的反应盐渍化的空间分布特征。③ 分地类和不分地类的最佳拟合方程大多是非线性的。④ 从农用地最佳拟合曲线可以得出拟合方程多为单调递减函数, 即随着含盐量的增加, 其对应物候参数的值在逐渐降低, 最明显的就是LSI、RD、LD、PV这几个参数, 这些参数都代表了光合作用的情况信息, 说明农用地植被对于盐分情况的变化极为敏感, 盐分越重对于农作物的生产量的影响越大。⑤ 对于稀疏植被, 情况更为复杂。大多是呈下开方向的二次多项式方程, 即先随着盐分的增加对应物候参数的值也在增加, 到达一个最大之后又随着盐分的增加而减小, 如EOS、AMP、BASE、PV、LD、RD和SSI。这就说明盐生植被的确具有一定耐盐性, 但是当盐分含量到达一定程度之后也会对其生长造成较大的影响。⑥ 通过以上5点的分析, 可以得出尽管通过深挖植被的信息可以较好的反演盐渍化信息, 但至少对于盐生植被区域只用植被的物候信息反演盐渍化信息的效果并不理想, 需要借用其他方面的参数来参与反演。

通过以上6点的分析可以得出植被信息只是反演盐渍化的一种重要指示剂, 但不能适用于所有植被类型, 用最佳曲线拟合的方式来反演盐渍化信息不具有普适性。基于以上的分析, 本研究运用BP-ANN人工神经网络来结合物候参数、地形信息、影像数据、野外实测数据等22个参数来反演盐渍化信息。

因为从表 3的结果来看指示光合作用的LSI和SSI指数与盐分含量的相关性较高, 这表明生物累积量的信号能更好反映盐渍化信息。可能主要是因为对于不同类型的植被对于盐渍化程度的响应不同, 同种植被对于不同盐含量的响应也不同, 不同植被对于同一水平的盐分含量响应也会不同, 但是这一切都会原原本本的记录在植被的体内, 最突出的就是它的长势情况, 最终就会导致生物累积量的不同, 所以LSI和SSI的表征盐渍化的能力要优于其他几个物候参数, 因此图 6中只显示LSI和SSI分布图。

从图 6中的LSI空间分布可以得出位于绿洲农用地植被区域的累积生物量明显大于绿洲荒漠交错带与荒漠地区。农用地的生物量值最高, 值域在4.9—9.2之间；绿洲荒漠过渡带的生物量值较低, 值域在1.1—4.9之间；荒漠地区生物量值最低, 也是盐渍化程度最严重的区域, 值域在0—1.1之间。图 6中SSI分布的情况与LSI类似, 不同的是对于同一区域的累积生物量的值要明显低于LSI的值, 如在农田区的SSI值大部分处于4.9—6.4之间, 而LSI值大部分在6.4—9.2之间。SSI值在绿洲和荒漠地区低于4.9的区域面积明显大于LSI低值区。

3.3 不分区BP-ANN神经网络盐渍化信息反演

把提取出来的10种物候参数 (SOS、EOS、AMP、LEN、BASEV、PV、LD、RD、LSI、SSI) 和12种非物候参数 (NDVI、SI、SDI、TVDI、MSAVI、TSAVI、DEM、坡度 (SLOPE)、TWI、蓝/红/近红外波段反射率) 作为BP-ANN的输入因子。由于含盐量和电导率之间有着极好的相关性, 可以用电导率表征含盐量^[13], 因此把野外实测的电导率 (1:5比例的浸提液电导率) 作为输出结果。隐藏层的神经元数目设置为14个。

具体过程如下:

① 采样利用野外实际采样的87个采样点, 其中50个样点作为训练样本, 37个样点作为测试样本。

② 预测2013年盐渍化分布情况利用Matlab创建BP神经网络, 利用训练好的网络来对测试数据和逐像元数据进行预测, 精度验证结果见表 4, 预测结果见图 7。

3.4 分区BP-ANN神经网络盐渍化信息反演

分别把位于农用地、稀疏植被区和其他地类的野外采样点提取出来。用位于不同地类的采样点来训练神经网络, 分别得到分区之后的结果图, 然后用各自地类的矢量边界分别裁剪对应的结果图, 最后在ArcGIS软件中进行镶嵌处理, 得到分区之后的盐渍化信息分布图。其中由于其他地类区的反演结果用不分区的结果代替。

具体过程如下:

① 采样按照不同地物类型采集样本, 然后按照如上的设置提取训练样本和测试样本 (实测数据中60%用于训练, 40%用于测试)。

② 预测2013年盐渍化分布情况利用Matlab创建BP神经网络, 利用训练好的网络来对测试数据和逐像元数据进行预测, 得到不同分区的预测结果, 结果见图 8, 镶嵌后的结果如图 9所示, 精度验证结果见表 4。

同时只用非盐渍化的12个参数来反演盐渍化分布情况, 精度比较结果见表 4。

从表 4中可以得出当无物候参数参与反演时可得：不分区和分区的总体差别不显著, 几乎一致。RMSE从11.55增加为12.30, ME绝对值从0.88增加为2.46, R²从0.68下降为0.67, 从这几个方面可知分区后的总体精度略有下降, 但是分区后农用地部分的反演精度较高, RMSE从不分区时的11.55陡降为农用地部分的4.21, R²也从0.68上升为0.73；对于稀疏植被, 分区后的精度反而明显降低, RMSE从11.55上升为14.22, R²从0.68下降为0.54, 可见对于稀疏植被区, 单纯的用物候参数进行反演效果不理想, 需要探讨用其他方法以及加入其他参数后来进行反演。当有物候参数参与反演时可得：分区后的精度明显要高于不分区时的精度, 如RMSE从不分区时的9.55减少到分区后的7.73, ME从0.42减少为0.04, R²从0.79升高到0.8787, 这个0.8787并没有位于0.79和0.8765(农用地) 之间, 主要原因是镶嵌时的其他地类中的采样点拟合精度高于不分区时的整体精度, R²达到了0.90, 所以分区后镶嵌的结果的R²更高。因此总体而言, 按照植被类型分区后的预测精度明显高于不分区时的精度。这可能主要是因为不分区时神经网络训练的过程中是对所有样点进行全局最优拟合, 在误差反馈的过程中权重的分配根据所有样点而变化, 但是如地理坐标投影一样, 局部区域的误差会变大, 因此神经网络训练好的结果难免会受到其他样点的影响, 精度较低。相反, 在按照不同植被类型进行分区后, 对每种植被区的样本进行专门训练就不会受其他植被类型样点的干扰, 拟合出的方程与匹配的权重更加符合局部区域的实际情况, 所以最终的预测精度较高。对比国内相关学者对于该区域的研究^{[6-7, 24]}, 发现其RMSE普遍较低, 实测值与预测值的相关性一般都在0.84以上, 但是对比国外相关学者^[13]的研究, 发现其RMSE普遍较高, 一般数值在20以上, 本文结果介于二者之间, 但是从相关性和可决系数R²上看, 本文反演结果的精度较高。主要原因可能是采样时间与研究的时段不同和研究方法上的差异等。但是国外的研究一般都是基于大数据量的样点统计出来的结果, 更加具有统计学意义, 因此以后应该对研究区多采集样点, 来更好的对比国内外研究的精度。

从图 9中可以得到以下结论：①图 9中有部分无值区, 这主要是因为对于无植被覆盖或分布极其稀疏的区域, 不能有效的提取物候参数所致。但是通过出野外的实际观测, 这些无值区一般都是重盐渍化区域。② 绿洲内部农用地区域主要为非盐渍化地和轻盐渍化地, 占有效像元的比例达到53.42%, 主要位于中海拔冲积扇平原、中海拔冲积洪积扇平原、低海拔冲积扇平原的地貌类型区。③ 位于研究区绿洲东部和东南部盐渍化程度最为严重, 集中于大于16ms/m。整个研究区的重度盐渍化地比例为32.87%, 主要位于低海拔固定草灌丛沙堆的地貌区域中。④ 研究区中度盐渍化比例为13.71%, 主要零星分布于绿洲内部, 较为集中的分布在研究区东北部和绿洲最南端, 所处的地貌类型主要是低海拔固定草灌丛沙堆。从以上分析可以得出盐渍化程度受高程和地貌类型的影响显著, 高程较低的地区盐渍化程度普遍较重, 较高的区域盐渍化程度普遍较低。因此对研究区以主要地貌类型进行分区后, 用BP-ANN预测后的的精度分析见表 5, 结果见图 10。

从表 5可得：① 当物候参数参与反演时, 按地貌类型进行分区后的预测精度有了明显的提高, RMSE从9.55减小为8.18, 大于按植被类型分区时的7.73；R²从0.79增加为0.85, 小于按植被类型分区时的0.88；ME绝对值从0.42减少为0.11, 大于按植被类型分区时的0.04, 偏离真实值的幅度较大。② 在无物候参数参与反演时, 按地貌类型分区后的预测精度略微下降, RMSE从11.55增加为13.52, R²从0.68减少为0.60, ME绝对值从0.88增加为3.78。从图 10中可以得出：盐渍化与地貌类型呈高度对应的聚集分布。在较高海拔区的中海拔冲积扇平原和中海拔冲积洪积扇平原的盐渍化程度较低, 集中于0—5.50ms/m。处于稀疏植被区和荒漠区的低海拔半固定草灌丛沙堆、低海拔固定草灌丛沙堆地貌类型的盐渍化程度较重, 主要位于研究区东北部、南部和西部绿洲边缘区, 电导率普遍大于16.49ms/m。低海拔冲积扇平原是研究区主要的地貌类型, 盐渍化普遍较重, 集中于5.50—16.49ms/m, 但位于研究区南部边的部分区域盐渍化较轻, 集中于0—5.50ms/m。位于研究区南部和塔里木河北岸交界的低海拔冲积河漫滩盐渍化较重, 但集中于小于34ms/m。综上分析可以得出基于地貌类型分区后的整体预测精度比根据植被分区低, 但显著高于不分区时的精度, 在以后的分析中以植被分区来作为依据更为适宜。

4 结论与讨论

通过反演得到10种物候参数和计算得到植被指数、盐渍化指数、地形参数等12种非物候参数作为BP-ANN模型的输入因子, 把在实验室实测的电导率数据作为输出数值, 最后反演得到盐渍化信息空间分布图, 得到以下结论:

(1) 盐渍化的形成受多种因素的影响, 与物候参数大多呈非线性关系, 不能单纯的以某拟合公式来进行表达, 需借助人工神经网络超强的非线性拟合能力来反演盐渍化信息。

(2) 代表生物累积量的LSI和SSI最能显著的表征盐渍化的变化情况。

(3) 通过深挖植被多方面的信息, 在融入物候参数后的盐渍化反演精度显著提高, 可决系数从0.68增加为0.79, 但是仅依据植被方面的信息不足以来精确反演盐渍化信息, 需要加入地形、影像数据和土壤水分空间分布情况等方面的信息。

(4) 按照植被类型分区后的盐渍化信息预测精度进一步提升, R²达到了0.88。极可能是因为通过植被分区研究后能够用BP-ANN模型拟合得到对应最佳局部的拟合方程, 因此十分有必要按照植被类型进行分区来反演盐渍化信息。

(5) 按照地貌类型分区后的预测精度相比于不分区时有显著的增加, RMSE从9.55减小为8.18, R²从0.79增加为0.85, 但综合预测精度比按照植被类型时偏低。

(6) 研究区非盐渍化和轻度盐渍化、中度盐渍化和重度盐渍化地比例分别为53.42%、13.71%和32.87%, 盐渍化地普遍存在, 研究区东部和东南部盐渍化最为严重。高程和地貌类型对盐渍化分布影响显著。

综上分析可得出融合物候特征和按照植被类型或地貌类型分区后盐渍化的反演精度显著提高。但本文的采样点较少, 只有87个样点, 在建立BP-ANN神经网络模型时数据量略显不足, 但是渭库绿洲本身绿洲-荒漠形态的界限鲜明, 在绿洲农用地区、绿洲荒漠交错带和荒漠对应区域的同质性较强, 因此野外实测的样点虽较少, 但是均匀分布在了这3个大区域中, 所以文中分析的结果基本符合实际情况。因此本文侧重点在于介绍一种结合物候与非物候参数来进行反演盐渍化信息的方法, 以及按照植被类型和地貌类型分区进行预测盐渍化信息的思维模式。物候参数代表盐渍化长期受影响的因素, 非物候参数主要代表短期内的影响因素, 通过探讨结合长期和短期内的影响因子来反演盐渍化信息, 最后得出这种结合方法的预测精度较高, 具有较高的应用潜力。以后应该要针对研究区多采集样点, 以更好的挖掘物候参数在预测盐渍化信息空间分布的能力, 并进一步探讨如同时考虑按植被和地貌分区后的预测精度, 以及按照植被和地貌类型分区后对于预测结果的影响机理。