王建华，李 阳，梁树能，孙小飞
（1.中国科学院空天信息创新研究院，北京 100094；
2.北京航空航天大学，北京 100191；
3.自然资源部国土卫星遥感应用中心，北京 100048；
4.西南交通大学，成都 611756）
沙化是指在各种气候条件下，由于各种因素形成的地表呈现以沙（砾）物质为主要标志的土地退化，具有这种明显特征的退化土地为沙化土地[1]。土地沙化是我国生态环境和社会经济的重大威胁[2]，截至2014年，全国沙化土地总面积1 721 175 km2，占国土总面积的17.93%，具有明显沙化趋势的土地面积300 293 km2，分布在内蒙古、新疆、青海等12个省[3]。由于土地沙化区域自然环境恶劣，生态系统脆弱，资源利用问题突出，需要对土地沙化地区及时监测[4]，采取保护、治理以及防治手段，避免土地沙化程度加重。
遥感技术广泛应用于土地资源调查、生态环境监测、灾害预报与灾情评估等领域。自20世纪70年代起，遥感技术开始应用于土地沙化研究，目前多以遥感和地面调查相结合的技术方法进行土地沙化监测[5,6]，其观测范围大、速度快、精度高、具有“图谱合一”的特性，为监测土地沙化发生、发展状况，及时有效采取防治土地沙化措施提供了切实可靠的途径，在土地沙化监测领域发挥着无可取代的作用[6]。
随着遥感技术的发展，高光谱遥感为开展全球性和区域性的沙化研究提供了新的手段，土地沙化监测开始向定量化遥感方向发展[7-10]。高光谱遥感技术在光谱分辨率上有着巨大的优势，利用高光谱数据反演得到的地物反射光谱特征，能表征地球表面物体的分类、物质的成分、含量、存在状态、空间分布及动态变化[11，12]。
研究区行政区划隶属于内蒙古鄂尔多斯市伊金霍洛旗，位于黄河中游风沙砒砂岩区域，地处陕北黄土高原与库布齐沙漠接壤地带，为毛乌素沙漠与黄土丘陵沟壑区的过渡地带，地貌类型主要为黄土丘陵沟壑地貌，沟壑纵横，支离破碎，流域上游少部分为黄土盖沙区，下游沿黄河河谷一带是基岩沟谷丘陵区，流域内的黄土性土壤和风沙土两类，前者结构疏松，透水性强，是地表径流的有利形成条件[13]（图1）。
图1 研究区位置示意图Fig.1 The location of study area
研究区域为干旱、半干旱大陆性季风气候，结合了鄂尔多斯高原和黄土高原的气候特点，风大沙多，降水集中且强度大，多年平均降水量约为430 mm。在地域分布上，区域自东南部向西南部，雨量减少；
在时间分布上，降水年际变化大且年内分配不均，多以暴雨形式出现，汛期（6—9月）降雨量可占全年降雨量的80%，七、八两月降雨量占年雨量的54%。高强度暴雨是流域内径流、泥沙产生的主要原因，洪水期沙量占全年总沙量的96%以上，汛期沙量占全年沙量的99%以上，水土流失极为严重，占全流域面积的97%[14，15]。
本文基于高光谱卫星数据的光谱细分特征优势，选用植被覆盖度（Fraction Vegetation Coverage，FVC）、裸土指数（Bare Soil Index，BSI）、地表反照率（Albedo）、改进土壤调整植被指数（Modified Soil Adjusted Vegetation Index MSAVI）四种土地沙化综合监测指标，结合人机交互解译工作进行了区域土地沙化监测，基于人机交互解译土地沙化监测结果建立决策树分类模型，对全部研究区域的土地沙化程度进行分类，实现了土地沙化信息的识别和提取。
资源一号02D卫星是首颗基于资源环境观测需求的高光谱业务卫星，定位于中等分辨率、大幅宽观测和定量化遥感任务，卫星光谱获取能力突出，可应用于地矿勘探、土地监管过程中地物信息定量化反演[16]，满足对于土地沙化信息的识别和提取的需求（表1）。
表1 ZY-1 02D高光谱卫星数据源数据特征Tab.1 Characteristics of hyperspectral satellite data
本文选取了三景资源一号02D卫星影像，成像时间为2021年6月30日，为同一轨道成像数据，遥感影像成像质量高；
重点监测的土地沙化区域色调突出，植被和土地沙化区分度大，易于土地沙化信息的识别及提取。
经过辐射定标、大气校正、几何校正、拼接镶嵌等预处理过程，得到最终的遥感卫星地表反射率影像作为本文研究的基础数据，如图2所示，遥感影像图RGB合成方式为波段148（2 199.482 9 nm）、波段55（859.678 4 nm）、波段20（559.234 3 nm）。
图2 研究区域高光谱遥感影像Fig.2 The hyperspectral remote sensing image of the study area
2.1 土地沙化遥感监测指标
高光谱遥感具有精细的光谱特征刻画能力，能够实现对土地沙化的指示性指标的准确提取，而这些指示性指标是土地沙化遥感监测的关键。土地沙化过程受内外环境多种因素的影响，其外在表现多有不同，这反映在多种多样的土地沙化监测指标上。针对于土地沙化遥感监测的单一指标无法完全表征土地沙化情况，如植被指数、裸土指数等，最终的土地沙化监测结果具有一定程度的局限性。本文在植被指数、裸土指数指标的基础上引入地表反射率等其他土地沙化监测指标以提高土地沙化评价效果[17-22]。为了实现研究区内大尺度的土地沙化监测，准确的反映出土地沙化的特征，通过高光谱卫星影像提取、反演得到的植被覆盖度、裸土指数、地表反照率、改进土壤调整植被指数，将这四个指标用于土地沙化遥感监测。
（1）植被覆盖度、裸土指数
植被覆盖度和裸土指数是评价沙化程度的重要指标，植被覆盖度越低，裸土指数越高，土地沙化程度越高。
本文采用完全约束最小二乘法混合像元分解模型进行植被覆盖度、裸土指数的提取。遥感影像中混合像元普遍存在，该像元分解模型法考虑单一像元的多个组分对遥感传感器观测信息的贡献，并且认为是每个组分的线性宏观组合，则高光谱影像中的每一个像元光谱即可表示为土壤、植被等光谱的线性组合，利用该模型获得植被和土壤分布，进行植被覆盖度和裸土指数的估算[23]。
研究选取的端元为植被、裸地和水体。端元光谱的一般来源包括参考光谱库和高光谱影像。本文采用直接从影像中提取端元，一般包括纯净像元指数（Pixel Purity Index，PPI）、连续最大角凸锥（Sequential Maximum Angle Convex Cone，SMACC）等方法[24]。本文中采用的三景数据为同一轨道成像的高光谱卫星数据，因此本文借助PPI和N维可视化工具，针对研究区域镶嵌后的三景高光谱影像进行端元波谱采集，端元提取的流程和结果见图3。
图3 高光谱遥感影像端元提取流程Fig.3 The process of endmember extraction
基于高光谱遥感影像中提取的端元，经过混合像元分解，得到植被、裸地和水体的丰度，完成植被覆盖度、裸土指数的提取，提取结果如图4所示。
（2）地表反照率
地表反照率是土地沙化监测中地表温度、干燥度或湿度的指示因子。地表反照率的计算方法依据大气辐射传输模型建立，地表反照率可表示为[25]：
式中，ρB、ρR、ρNIR、ρSWIR2分别为影像的11波段（481.829 5 nm）、32波段（662.286 9 nm）、52波段（833.873 6 nm）和149波段（2216.3213nm）的地表反射率，提取结果见图5a。
（3）改进土壤调整植被指数
植被指数的变化能够表征土地生产力的变化，从而体现土地沙化程度的动态变化，为减小低密度植被下的表达误差，消除土壤的影响，采用改进土壤调整植被指数更准确的表达研究区域植被生长状况[26，27]。该指数基于土壤调整植被指数，对于低植被覆盖，MSAVI比SAVI具有更强的土壤变化表征能力，对植被覆盖的敏感性增加超过50%。该指数适用于植被相对稀疏，能够通过冠层看到土壤的地区，改进土壤调整植被指数可表示为：
式中，MSAVI为改进土壤调整植被指数，ρR和ρNIR分别为影像的32波段（662.2869 nm）和52波段（833.8736 nm）的地表反射率。改进土壤调整植被指数提取结果见图5b。
图5 地表反照率（Albedo）和改进土壤调整植被指数（OSAVI）提取结果Fig.5 The result of Albedo and OSAVI
2.2 人机交互解译
人工交互解译方法依赖于先验知识，土地沙化监测结果精度较高，能够为本文提供研究区域内的有效、可靠的土地沙化监测结果参考数据，为各个土地沙化程度监测指标的判别提供依据。参考“地质环境遥感监测技术要求（1/25万）”[28]，采用三级分类方式，沙质荒漠化程度按风积、风蚀地表形态占该地面积百分比、植被覆盖度及其综合地貌景观特征划分为轻度、中度、重度3个级别，土地沙化程度的划分（见表2）。
表2 人机交互解译土地沙化程度划分表Tab.2 The division of land desertification degree
根据遥感卫星数据特征并结合地面调查先验知识，确立不同沙化程度的遥感影像解译标志[29]，并针对研究区内土地沙化程度多样且存在较多极重度沙化情况的南部地区进行了土地沙化遥感解译工作（图6）。从图6中可看出区域总体分布特征从空间分布上由北向南，沙化程度加重；
自西向东，沙化程度减轻。
图6 土地沙化遥感解译分布图Fig.6 The distribution of land desertification
2.3 土地沙化信息识别及提取
根据第四次全国荒漠化和沙化监测技术规定将沙化区域按气候类型划分为亚湿润干旱区、半干旱区、干旱区、极干旱区，不同气候区内按照土地沙化程度划分为轻度、中度、重度、极重度4级。
以人机交互解译的土地沙化监测结果为参考数据，经过数字化处理后，将该参考影像和高光谱反演出的到植被覆盖度、裸土指数、地表反照率、改进土壤调整植被指数这四种土地沙化遥感监测指标的影像进行叠加分析，依据各个监测指标在参考影像中的不同土地沙化程度等级，得各个指标对应的土地沙化程度判别值（表3）。
表3 土地沙化监测指标Tab.3 The indicators of land desertification monitoring
将各指标进行累加运算，利用决策树分类方法根据土地沙化监测指标分别进行研究区域土地沙化程度的判定，土地沙化监测决策树结构见图7。
图7 土地沙化监测决策树结构Fig.7 The structure of decision tree
基于研究区域的高光谱遥感影像，分别提取植被覆盖度、裸土指数、地表反照率、改进土壤调整植被指数等监测指标，由决策树分类模型得到研究区土地沙化程度分布图，如图8所示。
图8 土地沙化监测结果Fig.8 The result of land desertification monitoring
采用土地沙化综合监测指标以及决策树分类的方式，得到的研究区域土地沙化遥感监测结果（表4）。
表4 研究区域土地沙化监测统计表Tab.4 The statistical table of land desertification in study area
除耕地、建设用地以及水体外，极重度沙化区面积为271.63 km2，占总面积的2.71%，主要分布在乌审旗乌兰陶勒盖镇、榆林市小壕兔乡和孟家湾乡；
重度沙化区面积为1 685.85 km2，占总面积的16.80%，主要分布在乌审旗、榆林市、神木县以及工作区北部的达拉特旗；
中度沙化区面积最大，为5462.70 km2，占总面积的54.44%，遍布整个研究区域；
轻度沙化区面积为2 457.77 km2，占总面积的24.49%，主要分布在伊金霍洛旗的阿勒腾锡热镇、扎萨克镇以及乌兰木伦镇；
非沙化区面积为156.45km2，占总面积的1.56%，主要分布在伊金霍洛旗的阿勒腾锡热镇、扎萨克镇以及乌兰木伦镇，地形上主要分布在沟谷地带。
采用已有的人机交互解译结果对现有的决策树分类沙化监测结果进行验证（图9）。在定性验证方面，现有结果的极重度沙化区和原有结果的重度沙化区套合情况较好，边界明显。
图9 人机交互解译结果与决策树分类结果套合情况Fig.9 The nested condition of result
在定量验证方面，采用接收者操作特征曲线（Receiver Operating Characteristic Curve，ROC）来验证现有结果。接收者操作特征曲线是一种对于灵敏度进行描述的功能图像[30]，被广泛应用于包括环境评价在内的各个领域。在ROC曲线分析中，曲线下面积(The Area Under the ROCCurve，AUC)与验证精度之间的定性和定量关系可分为差（0.5～0.6）、一般（0.6～0.7）、良好（0.7～0.8）、非常好（0.8～0.9）和优秀（0.9～1）。
ROC曲线如图10所示，决策树分类结果的验证精度AUC值为0.603，对验证结果进行分析，沙化程度匹配不一致的主要原因为：土地沙化程度分类级别在人机交互解译和决策树分类过程中存在差异；
采用土地沙化遥感综合检测指标和决策树分类模型的方法在对土地沙化的精细刻化上比人机交互解译方法表现更加优越。
图10 土地沙化监测结果ROC曲线Fig.10 The ROC of land desertification monitoring
本文基于高光谱卫星数据对研究区进行了土地沙化遥感解译、沙化程度识别提取和评估，人机交互解译方法能够提供真实可靠的研究区域土地沙化监测结果，是决策树分类模型建立的基础，同时能够为结果验证提供数据支撑。同时利用高光谱卫星数据精细光谱特征的优势，反演了植被覆盖度、裸土指数、地表反照率、改进土壤调整植被指数，将这四种指标作为土地沙化综合监测指标，将土地沙化程度分为极重度沙化、重度沙化、中度沙化、轻度沙化四个级别，基于人工交互解译的参考数据确定各个监测指标的土地沙化程度判别值，采用决策树分类方法完成土地沙化程度的判定，效果较好。
总体来看，研究区域土地沙化总体分布广泛，从空间分布上由北向南，沙化程度加重；
自西向东，沙化程度减轻，非沙化地区主要集中于沟谷地带，土地沙化分布和区域地形地貌，水文特征情况相匹配，显示出高光谱遥感技术可为黄河中游地区的土地沙化定量反演及监测研究提供行之有效的技术支撑。