阿拉尔垦区土壤盐渍化遥感监测及时空特征分析

代云豪， 管 瑶， 张钦凯， 孙峻杰， 贺兴宏,2

（1.塔里木大学水利与建筑工程学院，新疆 阿拉尔 843300；2.塔里木大学南疆岩土工程研究中心，新疆 阿拉尔 843000）

土壤盐渍化问题不仅影响区域农业经济发展，也会对该区域的生态环境造成严重威胁［1］。利用科学的技术手段及时准确掌握土壤盐渍化信息和预防次生盐渍土，对区域发展有着现实的科学意义［2-3］。遥感技术从多光谱到高光谱，再到微波遥感的发展，一直是国内外学者的研究热点，遥感技术凭借快速、及时、准确、高效的优势被广泛应用于土壤盐渍化研究［4-7］。国外早期以土壤盐分敏感波段研究为主，Dwivedi等［8］提出蓝波段、红波段、中红外波段合成影像获取盐渍土信息量最佳。21世纪后，土壤盐渍化研究不在仅限于敏感波段研究，通过构建相应参数或模型对土壤盐渍化进行定性、定量研究也成为当下热点。Sidike等［9］将归一化植被指数（Normalized difference vegetation index，NDVI）、土壤亮度指数等作为参数应用于土壤盐分反演模型以提高反演精度。El-Horiny［10］研究了电导率与光谱指数的关系性，表明盐度指标（SI-1和SI-2）与电导率的相关性最高，模型对研究区的预测准确率为77%。我国土壤盐渍化多发生在干旱、半干旱和滨海等区域，研究的热点区域主要位于河套灌区［11-12］、黄河三角洲［13-15］、新疆维吾尔自治区［16-19］等地，目前通常采用土壤光谱特征或是通过间接特征（植被长势、土壤要素、土地覆被类型等）对其进行研究［20］。研究中，多光谱遥感中光谱指数方便且易获取，成为土壤盐渍化研究热点之一，常以NDVI、修改型土壤调整植被指数（Modified soil-adjusted vegetation index，MSAVI）、盐分指数（Salinity index，SI）等建立二维特征空间和构建遥感盐分监测指数（Salinization detection index，SDI）模型，该方法反演研究区土壤盐分信息是土壤盐渍化遥感监测的一种先进方法［21-23］。卢晶等［24］采用MSAVI-SI 特征空间并构建改进型盐碱化监测指数（Modified salinization detection index，MSDI）模型反演河套灌区盐度，其精度可靠性较高；王飞等［25］采用NDVI-SI特征空间构建SDI 模型能准确反演塔里木南缘于田绿洲土壤地表盐量。其中SDI模型由NDVI 和SI 构建，而应用二维特征空间理论来发展盐渍化遥感监测的定量方法和指标有巨大的潜力，利用植被指数与SI 构建SDI 模型能有效对土壤盐渍化进行反演［26-28］。

阿拉尔垦区是受土壤盐渍化侵害的典型绿洲农田，土壤盐渍化分布较为广泛［29］。垦区作为新疆重要的粮棉产地和兵团向南发展的重要部署地，有着其重要的核心地位，但多年来局限于垦区棉田或部分试验田土壤盐渍化研究，系统研究整个阿拉尔垦区土壤盐渍化时空分布较少。因此，为进一步掌握垦区土壤盐渍化分布状况，提高土壤盐渍化精准治理与防治，本研究选取SI 和NDVI 建立二维特征空间并构建SDI 模型，利用该模型对阿拉尔垦区土壤盐渍化进行反演，反演结果可作为垦区土壤盐渍化治理的一种参考，对促进经济与社会发展具有一定的现实意义。

1 材料与方法

1.1 研究区概况

阿拉尔垦区位于新疆南部（80°30′～81°58′E，40°22′～40°57′N），属塔里木河冲积细土平原，气候条件干燥，生态环境脆弱，傍依塔里木河，地势由西北向东南倾斜。垦区辖七团、八团、九团、十团、十一团、十二团、十三团、十四团8个团场，拥有3个人工水库，包括多浪水库、胜利水库、上游水库。垦区总面积4197.58 km2，海拔高度997～1340 m，年降雨量稀少，地表蒸发强烈，年均气温4～20 ℃，年均降水量40.1～82.5 mm，年均蒸发量1876.6～2558.9 mm［30-31］。垦区种植农作物以棉花为主，林果则以红枣、核桃、香梨为主，农业灌水方式以冬春漫灌洗盐结合膜下滴灌为主（图1）。

1.2 数据来源

影像数据来源于地理空间数据云（http://www.gscloud.cn），筛选云量小于10%的阿拉尔垦区2011年8月20日Landsat 7 ETM+和2021年9月8日Landsat 8 OLI 影像数据，影像数据像元为30 m×30 m，影像行列号为146/32。

野外土壤样品主要采集0～10 cm 表层土壤，采集时间为2021年3月8—14日、2021年8月23日—9月3 日，其中：2 次土壤采样点分别为42 个和97 个，共计139 个（图1）。采样点随机均匀分布于垦区七团至十四团内，满足均匀布设原则。采样点地物类型为垦区典型种植农作物（棉花样点69 个、玉米样点2 个、辣椒样点6 个）、林果（枣树样点42 个、梨树样点8个、核桃树样点2个）和裸地（10个）。为减小单一样品的不确定性，每个样点根据三角形采样法［24］相距12～15 m，利用土钻采集3份样品各装入铝盒，利用GPS 在三角形中点处测定样点坐标，野外数据共采集土壤样品417份。

图1 阿拉尔垦区概况及采样点分布Fig.1 Overview and distribution of sampling points in Aral Reclamation Area

1.3 研究方法

1.3.1 数据处理原始遥感影像经几何校正、辐射定标、大气校正等预处理工作，裁剪得到研究区影像，其中2011 年影像数据条带丢失，为解决这一问题，采用ENVI 5.3软件中landsat_gapfill插件处理［32］。

野外土壤样品称重后放入105 ℃烘箱烘干8 h，将同一采样点的3份样品均匀混合磨碎成为1份样品，剔除植物根系、石块等杂物后过2 mm 筛。取土20 g，蒸馏水100 mL 充分混合搅动后静置过滤，制备得到的滤液用上海科佑DDS-11A 电导仪测定溶液电导率。

1.3.2 遥感光谱指数

（1）归一化植被指数（NDVI）

尽管植被覆盖会改变土壤的光谱响应，但植被生长状况在很大程度上受到盐分的影响和控制，可很好的反应土壤盐渍化程度。选取NDVI作为反应土壤盐渍化信息间接指标，其植被长势情况与土壤盐渍化程度成反比，土壤盐渍化程度越重，植被长势将受到抑制［33］。

式中：ρred、ρnir分别为landsat ETM+/OLI影像中的红光波段、近红外光波段。

（2）盐分指数（SI）

SI由蓝光波段和红光波段计算得到，Khan等［34］研究发现多光谱遥感中的蓝、红光波段对土壤盐分较为敏感，利用SI能较好反演盐分程度。

式中：ρblue为landsat ETM+/OLI影像中的蓝光波段。

1.3.3 掩膜处理阿拉尔垦区紧邻塔克拉玛干沙漠，垦区内有多浪水库、胜利水库、上游水库和塔里木河等水域，为避免水体和沙地对土壤盐渍化信息提取产生干扰影响，故利用SWIR 波段、NDVI 和湿度指数（Normalized difference moisture index,NDMI）对水体和无植被沙地进行掩膜并去除（图2）［24］。

图2 水体及无植被沙地提取Fig.2 Extraction of water bodies and unvegetated sandy land

式中：ρswir1为landsat ETM+/OLI 影像中的短波红外1波段。

1.3.4 归一化处理由于各指数的量纲不同，建立的模型需消除各指数之间的差异，故对各指数归一化处理［25］。

式中：SImax、SImin分别为盐分指数的最大值和最小值；NDVImax、NDVImin分别为归一化植被指数的最大值和最小值。

2 结果与分析

2.1 遥感盐分监测指数模型

以阿拉尔垦区2021 年9 月遥感影像数据计算的归一化SI 为坐标纵轴，NDVI 为坐标横轴建立SINDVI二维散点图［15］（图3），并计算二维散点的拟合度（表1）。由表1 可知，所有模型拟合度R2均大于0.8600，适合构建SDI 模型［25］，并能在宏观尺度上反演阿拉尔垦区土壤盐渍化信息。SDI模型计算公式如下：

表1 二维散点模型拟合性Tab.1 Fitting of two-dimensional scatter model

图3 SI-NDVI散点图Fig.3 SI-NDVI scatter diagram

2.2 阿拉尔垦区土壤盐渍化时空动态变化

2.2.1 土壤盐分遥感监测指数模型精度分析通过测定139 个野外土壤样点数据的电导率，并对SDI模型进行精度验证（图4）。野外实测土壤电导率与SDI模型之间的拟合度R2=0.7579，拟合效果较好，表明SDI模型可以对阿拉尔垦区较大范围的土壤盐渍化信息进行反演。

图4 实测电导值与SDI模型精度验证Fig.4 Verification of measured conductivity value and SDI model accuracy

为便捷、直观了解垦区土壤盐渍化程度，查阅南疆地区遥感盐渍土相关研究［23,35］，参考王亲遵等［36］土壤盐渍化等级划分与土壤盐渍化等级标准［37］所对应的土壤电导率值，利用SDI 模型与实测电导率构建的数学关系模型对阿拉尔垦区土壤盐渍化进行等级划分（表2）。

表2 土壤盐渍化等级标准划分Tab.2 Classification of soil salinization grade

2.2.2 阿拉尔垦区土壤盐渍化空间分布按标准电导率对应的SDI 模型对土壤盐渍化进行分级，得到阿拉尔垦区土壤盐渍化分布图（图5），计算分析得到2011—2021 年土壤盐渍化面积（表3）。2011 年阿拉尔垦区的非盐渍土面积达到972.98 km2，占垦区研究土壤面积的29.37%，主要分布在垦区西部和东部地区；轻度盐渍土面积251.09 km2，主要分布在西北、东北和西部地区，其分布呈零散化；中度盐渍土主要在南部塔里木河河床、水库附近；重度盐渍土和盐土主要分布在垦区东北地区和东南地区。2021年阿拉尔垦区非盐渍土区域分布明显扩大，面积达到1291.20 km2，其分布范围在垦区西部、北部和东部地区；轻度盐渍土主要分布在西部和北部地区；中度盐渍土主要在垦区中部地区和东北地区；重度盐渍土主要分布在垦区东北地区，东南地区重度盐渍土分布面积减小；盐土的分布同2011年大致相同。

图5 阿拉尔垦区土壤盐渍化分布Fig.5 Distribution of soil salinization in Aral Reclamation Area

由表3 可知，2011—2021 年阿拉尔垦区土壤盐渍化分布中非盐渍土、轻度盐渍土、中度盐渍土、重度盐渍土、盐土面积分别呈“上升、上升、下降、下降、上升”趋势，垦区非盐渍土面积增加318.22 km2，增幅达到32.70%，主要增加面积位于垦区北部和东部地区。轻度盐渍土增加0.80 km2，变化极小，主要呈零散化分散在垦区内。中度盐渍土和重度盐渍土分别减少76.84 km2、153.03 km2，面积明显减少区域主要在垦区西南、西北和东南地区。盐土面积增加了68.61 km2，增加区域位于垦区东北地区。

表3 2011—2021年阿拉尔垦区土壤盐渍化面积统计Tab.3 Statistics of soil salinization area in Aral Reclamation Area from 2011 to 2021 /km2

根据土壤实测数据可知，样品数量关系呈现为非盐渍土样品数量＞轻度、中度盐渍土样品数量＞重度盐渍土、盐土样品数量，其中非盐渍土的样品数量占总样品数量的72%。结合奥维地图分析2021年阿拉尔垦区土壤盐渍化空间格局分布，模型反演的非盐渍土以农田为主，轻度、中度盐渍土一部分零散分布于农田，一分部分布于塔河附近，重度盐渍土和盐土则多分布于沙漠附近。查询该区域相关土壤盐渍化研究，遥感解译与实际调查现状相符。因此，SDI模型能用于本研究区域土壤盐渍化研究。

2.2.3 阿拉尔垦区土壤盐渍化时空动态变化分析利用ArcGIS 10.3，生成2011—2021年阿拉尔垦区土壤盐渍化面积转移矩阵（表4）和土壤盐渍化动态变化图（图6）。由表4和图6可知，土壤盐渍化自2011—2021 年发生明显良好逆转的区域主要分布于垦区北部和南部等地区，其中非盐渍土主要由中度盐渍土转入155.89 km2，由重度盐渍土转入183.29 km2；轻度盐渍土主要由非盐渍土和中度盐渍土分别转入87.83 km2和60.40 km2，其分布呈零散化分布；中度盐渍土主要由非盐渍土和重度盐渍土转入58.94 km2和105.19 km2，其主要在垦区西部和东北方向。土壤盐渍化向严重趋势发展的区域主要位于垦区中部地区和东部地区，其中重度盐渍土主要由中度盐渍土和盐土分别转入70.42 km2和140.02 km2，盐土主要由重度盐渍土转入94.15 km2。

图6 2011—2021年阿拉尔垦区土壤盐渍化动态变化Fig.6 Dynamic changes of soil salinization in Aral Reclamation Area from 2011 to 2021

表4 2011—2021年阿拉尔垦区土壤盐渍化面积转移矩阵Tab.4 Transfer matrix of soil salinization area in Aral Reclamation Area from 2011 to 2021 /km2

3 讨论

3.1 模型选取

以NDVI 和SI 构建SDI 模型反演阿拉尔垦区土壤盐渍化精度为82.73%，与王飞等［25］采用此模型反演塔里木盆地南缘地区土壤盐分的精度基本一致，表明模型较适用于反演阿拉尔垦区土壤盐渍化。不同研究区自然条件、社会经济因素具有较大差异，南疆地区独特的气候条件和地位位置，致使SDI模型一定程度上反演内地土壤盐渍化精度有所不同［15］。SDI 模型较多，但是否适用于反演阿拉尔垦区土壤盐分，其适用性还需进一步的实验验证。总之，应充分考虑不同研究区地理位置等情况选择最合适的光谱指数构建SDI 模型，反演的土壤盐渍化信息也将会更加科学、合理。

3.2 植被指数与土壤盐分

植物生长状况与土壤盐分有着密切相关性，土壤盐分含量超过一定阈值便会影响植被生理参数，植被中所对应的光谱反射率便会发生改变，鉴于此，诸多学者认同植被指数能间接反映土壤盐渍化程度［36］。研究区计算NDVI值与土壤盐渍化分布有着较强的相似性，NDVI 值越高，SDI 模型反演土壤盐渍化等级越低，NDVI 与SI 构建二维特征空间散点图便可以直观观察植被与盐分关系成反比，随着土壤盐分增加，NDVI 值呈下降趋势，这与众多学者研究结果大致相同［13,15,24-25］。当然，植被指数种类繁多，NDVI 仅为其中最常用的一种，其余植被指数对垦区土壤盐分反演的精度还有待进一步验证。同时，南疆干旱区的植被生长期主要集中于6—10月，其余时间研究土壤盐分便需要考虑其他光谱指数或直接测定土壤光谱信息。

3.3 近10 a土壤盐渍化变化及驱动因素分析

为合理解释与辅助分析2011—2021 年阿拉尔垦区土壤盐渍化变化情况，由农业部规划设计研究院设施农业研究所温室数据共享平台（http://data.sheshiyuanyi.com/）获取阿拉尔垦区近10 a气候数据（图7），由兵团年鉴获取阿拉尔垦区近10 a 有效灌溉面积（图8）。

图7 近10 a阿拉尔垦区年均气温及年降水量Fig.7 Average annual temperature and annual precipitation in Aral Reclamation Area in recent 10 years

图8 近10 a阿拉尔垦区有效灌溉面积Fig.8 Effective irrigation area of Aral Reclamation Area in recent 10 years

分析阿拉尔垦区年均气温、年降水量和有效灌溉面积可知，年均气温波动范围在1～2 ℃之间，年降水量波动范围在20～100 mm，有效灌溉面积波动范围在0.5×105hm2。总体变化趋势为年均气温上升，年降水量降低，有效灌溉面积增加。结合SDI 模型分析的2011—2021 年非盐渍土、轻度盐渍土、中度盐渍土、重度盐渍土、盐土面积分别呈“上升、上升、下降、下降、上升”趋势，其中非盐渍土面积增长变化最为明显，盐土面积稍有增长。由上述分析得知，土壤盐渍化受自然因素和人为因素影响，阿拉尔垦区气候因素对土壤盐渍化变化影响较小，人为因素对土壤盐渍化变化影响较大，人为因素中灌溉面积则是造成垦区土壤盐渍化变化的主要驱动因素。

3.4 土壤盐渍化分布及其原因分析

土壤盐渍化分布受自然因素和人为因素影响［13］，阿拉尔垦区位于新疆南部干旱区，降雨量少、日照时间久、光照强、蒸发量大。由土壤盐渍化分布可知垦区农田以非盐渍土和轻度盐渍土为主，分布情况同吴家林［29］研究垦区0～30 cm棉田土壤盐渍化结果大致相同。垦区农田常年在春季或冬季采用漫灌洗盐，作物生长期常采用膜下滴灌控制水量和土壤盐分，随着种植年限增加，盐分快速向下淋洗，极大降低了农田表层土壤的盐分含量［38］。垦区农业灌水需求量大，塔里木河区域输入水量相对较少，灌溉后河床、草地等地区由于蒸发强烈，土壤表层水分快速蒸发，下层土壤水分携带着盐分向表层迁移聚集，造成表层盐分积累［39］，以中度盐渍土和重度盐渍土为主。土地荒漠化和紧邻塔克拉玛干沙漠区域附近，土壤因长期缺水，日照强烈，盐分多以盐斑集中呈现在土壤表层，此时表层土壤几乎无植被生长。当然，垦区受水文地质单元、植被类型、灌水水质、滴灌带质量等影响，农田局部区域水量不一致，土壤表层盐分含量有所差异，导致同一农田有不同程度的盐渍土。

此外，SDI 模型对阿拉尔垦区土壤盐渍化监测受Landsat 8遥感影像分辨率影响，只能对大范围进行反演，得到一个近似区域，对小范围土壤盐渍化精准反演，可考虑分辨率更高的遥感影像数据。

4 结论

（1）阿拉尔垦区农田以非盐渍土和轻度盐渍土为主，塔里木河河床、草地、林带等以中度盐渍土、重度盐渍土为主，土地荒漠化和沙漠附近以盐土为主。

（2）2011—2021 年，阿拉尔垦区土壤盐渍化等级变化中非盐渍土、轻度盐渍土、中度盐渍土、重度盐渍土、盐土面积分别呈“上升、上升、下降、下降、上升”趋势，总体可知近10 a阿拉尔垦区土壤盐渍化得到较好改善。

（3）阿拉尔垦区土壤盐渍化程度发生明显良好逆转区域主要位于垦区的北部和南部地区，非盐渍土主要由中度和重度盐渍土转入。垦区中部和东部地区土壤盐渍化程度加重，以重度盐渍土和盐土之间相互转化为主导。