土壤风蚀是土壤颗粒在风力作用下发生蠕移、跃移、悬移的动态过程^［1］，是干旱、半干旱区耕地土壤养分流失的重要原因^［2］。作物收割后留在地表的残茬是土壤风蚀影响因子——粗糙干扰因子的重要组成部分^［3］，其对农田土壤风蚀影响体现在两方面，其一可以减小风力及跃移颗粒对地表的剪切力，其二可减小地表裸露的面积，从而抑制农田土壤风蚀现象发生^［4-5］。在风洞试验中，同等风速条件下，翻耕地土壤风蚀现象较为明显，而免耕地几乎不发生土壤风蚀^［6］，翻耕地风蚀量是免耕地风蚀量的3~8倍^［7］。在野外观测中，农田地表输沙通量随作物残茬覆盖度增加呈指数形式下降^［8］，当作物残茬覆盖度达到35%及以上，土壤风蚀状况会明显降低，作物残茬覆盖的保护性耕作农田风蚀速率较传统耕作降低20%~70%^［9］。因此，获得准确作物残茬覆盖度是提高农田土壤风蚀量计算精度的关键，可为农田风蚀状况评估和保护性耕作措施实施提供必要数据支撑。

遥感为作物残茬覆盖度研究提供了新的手段，当前土壤风蚀研究中，普遍使用表征绿色植被特征的归一化植被指数（NDVI）计算残茬覆盖度，并作为植被覆盖因子输入到土壤风蚀模型中计算风蚀量，但由于风蚀季的作物残茬普遍缺乏叶绿素，NDVI对其表征能力不足，从而影响残茬覆盖度的计算精度，进而导致风蚀量计算结果出现偏差。作物残茬和土壤在400~1 100 nm波长光谱特征较为相似，需结合实测光谱数据进行进一步分析，寻找作物残茬和土壤在特定波段处光谱曲线的差异从而区分两者^［10］，作物残茬的荧光含量在440 nm波段处高于土壤^［11］，其和土壤的荧光含量相对差值在不同含水率条件下保持一致^［12］，结合荧光值检测方法利用400 nm处波段可对90%以上短轮作周期作物残茬进行有效辨识^［13］。枯黄植被的短波红外反射特性在区分作物残茬与土壤方面显示出较好效果，作物残茬中广泛存在的木质素、半纤维素和纤维素成分在2 100 nm附近的吸收峰表现出与土壤反射率上的差异^［14-15］。基于这些光谱特征开发了归一化差异耕作指数（Normalized Difference Tillage Index，NDTI）^［16］、归一化差分作物指数（Normalized Differential Senescent Vegetation Index，NDSVI）^［17］、归一化差值残茬指数（Normalized Difference Residue Index，NDRI）^［18］等两波段光谱指数，并建立了这些指数与残茬覆盖度的关系，在光谱指数基础上引入雷达数据构建综合光谱指数，提升了作物残茬覆盖度估算精度^［19］。两波段光谱指数形式简单，但易受背景噪声干扰，面对复杂特征时易饱和，导致灵敏度下降^［20-21］，综合光谱指数能更好地捕捉作物残茬光谱信息^［22-23］，但其形式复杂且需要不同来源遥感数据支撑，影响了使用便捷性。本研究以分析残茬光谱特征及其与覆盖度相关关系为切入点，通过第三波段引入补偿机制，构建抗干扰能力强、可抑制饱和效应的作物残茬光谱指数，以提高作物残茬覆盖度计算精度^［24-25］，克服NDVI指数在计算农田土壤风蚀状况评估及预报的精度方面的不足。

科尔沁沙地位于中国北方典型半干旱脆弱生态区，横跨内蒙古自治区东部、吉林省西部及辽宁省西北部，是距京津地区最近的天然风沙源，也是科尔沁、浑善达克沙地歼灭战主战场。受大陆性季风气候与地表覆被特征共同影响，该区域在春季易发生严重土壤风蚀现象，其土壤风蚀强度居中国北方风蚀区前列，剧烈的风蚀过程不仅导致农田表层土壤有机质下降、农业减产，也对周边地区生态环境造成严重威胁。历史上，过度开垦、放牧等不合理生产活动对该地区生态环境造成了极大的危害，加剧了该地区的土壤风蚀程度，在“三北”工程、京津风沙源治理工程等一系列重点生态工程的治理下，该地区发生了“人进沙退”“林进沙退”的实质性改变，但生态状况仍然存在脆弱性、不稳定性和反复性等问题。本研究通过分析科尔沁沙地的作物残茬实测光谱特征，构建反映作物残茬覆盖度的光谱指数，结合遥感影像建立风蚀季作物残茬覆盖度计算模型，提高研究区作物残茬覆盖度及农田土壤风蚀评估精度，助力科尔沁沙地歼灭战。

本研究以地处科尔沁沙地南缘的奈曼旗为研究区，奈曼旗（42°14′—43°32′N、120°12′—121°35′E）位于内蒙古自治区通辽市西南部。该地区属北温带大陆性季风干旱气候，年均气温7.79 ℃，年平均降水量354.3 mm，降水分配不均匀，集中在夏季（6—8月），年蒸发量1 982.8~2 081.8 mm，年均风速为3.5 m·s^-1，春季风速最大，且多大风，8级以上大风日数年均21次，土壤类型主要为风沙土。研究区主要农作物类型为玉米，是通辽市黄玉米生产地带和内蒙古自治区玉米种植示范区，《奈曼旗2022年国民经济和社会发展统计公报》公布的数据表明，奈曼旗玉米种植面积占耕地的87%，其余为小面积分布的小麦、高粱、薯类等作物，耕作制度为一年一熟，5月上中旬播种，10月中下旬收割，农田地表每年约有6个月时间缺少植被覆盖，风蚀季农田地表多为作物残茬覆盖，因此，对研究区风蚀季玉米残茬覆盖状况的准确识别与计算将有助于科学评估研究区农田土壤风蚀状况。

于2023年3月25日至4月5日开展野外调查及观测，根据实地调查的玉米残茬覆盖度差异，共布设12个200 m×200 m样区（图1A），每个样区嵌套3个10 m×10 m样方，每个样方内随机分布3个1 m×1 m样点，共采集98个样点数据。在晴朗无风天气下，利用ASD Fieldspec 3的25°视场角探头在垂直距残茬1.2 m处测定样点光谱信息，每次测定前均进行白板标定（图1B），通过变换探头位置连续采集15条光谱曲线，同时，在样点附近选择无残茬覆盖区域采集20~30条土壤光谱曲线。使用数码相机垂直向下拍摄样点照片，使用ENVI对照片进行监督分类以区分玉米残茬和裸土，通过玉米残茬所占照片比例计算玉米残茬覆盖度（图1C）。

使用ASD Fieldspec3配套软件处理野外实测玉米作物残茬高光谱数据。首先，剔除受噪声和大气水汽影响较大的波段范围（1 350~1 430、1 795~1 970、2 400~2 500 nm）；其次，取每个样方采集的多条光谱曲线平均值作为该样方光谱曲线；最后，使用Savitzky-Golay滤波对光谱反射率数据进行平滑处理（图2）。

遥感数据来自Sentinel-2A，该数据具有13个波段（B₂、B₃、B₄、B₈空间分辨率为10 m，B₅、B₆、B₇、B_8a、B₁₁、B₁₂空间分辨率为20 m，B₁、B₉、B₁₀空间分辨率为60 m），时间分辨率为5 d。在GEE平台调用研究区L2A级“COPERNICUS/S2_SR”数据集，数据已经过几何精校正和大气校正处理，筛选云量少于20%的影像，利用平台提供的去云函数去除其余云量，利用重采样工具将数据空间分辨率统一为10 m，生成研究区2018—2022年风蚀季（10月至次年5月）逐月数据产品。

将实测玉米残茬的高光谱反射率进行多种组合运算，构建不同玉米残茬光谱指数（表1），将每个光谱指数与对应的玉米残茬覆盖度进行相关性分析，并计算相关系数（R²），绘制不同指数与残茬覆盖度相关性矩阵图，选取相关性最高的波段范围，为残茬覆盖度计算模型提供建模依据。

根据不同波长组合构建的玉米残茬光谱指数与残茬覆盖度的相关性矩阵图，选取二者相关系数较高的波段组合，对比构成对应光谱指数的高光谱数据特征波长与Sentinel-2A多光谱数据各波段中心波长，确定用于计算玉米残茬覆盖度的Sentinel-2A对应波段，构建Sentinel-2A多光谱数据与残茬覆盖度的关系，通过决定系数（R² ）和均方根误差（RSME）确定研究区玉米残茬覆盖度的最优遥感计算模型，并计算研究区2018—2022年风蚀季逐月玉米残茬覆盖度。

土壤风蚀量受区域风力侵蚀因子、土壤抗蚀因子及粗糙干扰因子等多种因素的影响，其计算方法与评价指标也存在多种量化方式。本研究采用修正土壤风蚀方程（RWEQ）计算研究区的土壤风蚀量^［31-32］：

式中：SL为实际土壤风蚀量（kg·m^-2）；s为关键地块长度（m）；z为下风向距离（m）；Q_max表示风力的最大输沙量（kg·m^-1）；WF表示气候因子（kg·m^-1）；EF为土壤可蚀性因子；SCF为土壤结皮因子；K′为土壤粗糙度因子；COG为结合残茬因子。本研究风蚀量计算涉及的数据来自美国国家海洋和大气管理局、美国国家气候数据中心、中国寒区旱区数据资源中心、地理空间数据云。

利用本研究构建模型计算的玉米残茬覆盖度（FVC）计算结合残茬因子：

式中：COG表征作物残茬覆盖下的土壤流失比率。为评估本文构建的玉米残茬覆盖度计算方法对研究区土壤风蚀量计算的影响，另外，基于其他研究中普遍采用的NDVI像元二分模型计算的作物残茬覆盖度（FVC_NDVI ）计算COG，并输入RWEQ计算土壤风蚀量，对比两种不同方法计算的玉米残茬覆盖度导致的风蚀量差异。

根据实地调查状况，将样点玉米残茬覆盖度划分为5%~10%、10%~20%、20%~30%、30%~40%、40%~50% 5个等级，计算各覆盖度等级实测高光谱数据的平均反射率，并与土壤光谱曲线平均值进行对比。整体上，玉米残茬光谱反射率随波长增加呈先增加后减小趋势，反射率在350~1 350 nm波长范围随覆盖度增加而上升，在1 431~1 750 nm波长范围增幅缓慢；在1 971~2 130、2 125~2 200 nm波长范围内呈现双反射峰，主要是由于木质素和纤维素的特征光谱吸收作用（图2A）。土壤光谱反射率平均值变化与残茬较相似，但在350~1 350 nm波长范围内，土壤光谱反射率低于玉米残茬，在1 971~2 400 nm波长范围内，土壤光谱反射率高于玉米残茬。不同覆盖度玉米残茬及土壤反射率一阶微分显示，不同覆盖度玉米残茬一阶微分光谱在1 000、1 450 nm处及1 100~1 200、1 600~1 800、2 000~2 400 nm波长范围呈现显著零点穿越特征，可作为玉米残茬覆盖度计算的特征波段（图2B）。

在基于实测高光谱数据构建的不同指数与玉米残茬覆盖度相关性二维矩阵图中，基于不同光谱指数计算的玉米残茬覆盖度相关性具有相似的分布趋势，由在1 400~1 600、680~1 000、1 650~1 800、750~900 nm波长范围内反射率构成的NDI₁、RI、NDI₂与玉米残茬覆盖度的相关性较高（图3A、3B、3C）；在立体三维矩阵图中，由在600~1 000、800~1 000、1 000~1 200、2 000~2 200、2 100~2 400 nm波长范围内反射率构成的DFI与玉米残茬覆盖度的相关性较高（图3D），由在350~500、500~800、500~1 000、800~1 000、2 000~2 200、2 100~2 400、2 200~2 400 nm波长范围内反射率构成的CAI与玉米残茬覆盖度的相关性较高（图3E）。整体上看，NDI₁、RI、NDI₂等光谱指数与玉米残茬覆盖度的相关系数可达0.58，DFI和CAI等光谱指数与玉米残茬覆盖度的相关系数可达0.63，后者引入的第三波段，可有效提升玉米残茬覆盖度的计算精度。根据相关性分析结果，分别选择相关性较高的波段组合并与多波段遥感数据进行对比，用于玉米残茬覆盖度遥感指数构建。

基于与玉米残茬覆盖度相关性较高光谱指数的构成特征波长，结合高光谱数据与Sentinel-2A多光谱数据各波段中心波长的对应关系，选择Sentinel-2A对应波段代替高光谱数据重新构建各光谱指数，所选波段分别为短波红外特征波段B₁₀、B₁₁、B₁₂，近红外波段B₈、B₆，可见光波段B₄、B₁，基于Sentinel-2A光谱指数与玉米残茬覆盖度线性回归分析结果见表2和图4。

由B₁₂、B₆和B₁波段组合计算的CAI指数与利用样点照片实测玉米残茬覆盖度回归分析结果的决定系数R² 最高，为0.77，均方根误差为5.77%，较NDI₁、RI、NDI₂等两波段光谱指数对玉米残茬覆盖度的计算精度平均可提升约48%，RMSE平均可降低约30%，可相对准确实现对玉米残茬覆盖度的计算。故基于CAI指数进行研究区玉米残茬覆盖度和风蚀量计算。

研究区2018—2022年风蚀季平均玉米残茬覆盖度约为19.3%，总体呈现低覆盖度状况。从季节上看，研究区玉米残茬覆盖度秋季最高，冬季次之，在土壤风蚀主要发生的春季最低。从月际变化上看，研究区玉米残茬覆盖度10—12月总体降低，次年1—3月缓慢升高，3—5月又呈快速降低趋势（图5）。主要原因为研究区在9—10月进行农作物收获，作物收获后的大部分残茬残留在农田地表，随着时间推移，受到大风、降雪覆盖等天气影响，玉米残茬覆盖度总体降低，而在次年1—3月地表积雪融化，玉米残茬覆盖度缓慢升高，受当地耕作制度影响，研究区4—5月进行农田翻耕、播种，原残茬覆盖的农田变为了裸露的状态。

研究区2018—2022年风蚀季玉米残茬覆盖度存在明显空间分异特征（图6）。研究区玉米残茬覆盖度主要在0~20%，占研究区玉米种植面积的69%；其次为20%~40%，占研究区玉米种植面积的28%，覆盖度在40%以上地区面积较小。从整体上看，研究区南部黄花塔拉苏木、义隆永镇、新镇和北部东明镇风蚀季玉米残茬覆盖度较高，中部固日班花苏木、白音塔拉苏木和西部的苇莲苏乡等地区风蚀季玉米残茬覆盖度较低。研究区地处科尔沁沙地南缘，中部及东西两侧区域主要为沙地，农田多集中分布在人口集中的南部地区，玉米残茬覆盖度的空间分布与研究区南、北部耕地密集、中部耕地稀疏的分布情况一致。

基于玉米残茬覆盖度计算的研究区2018—2022年多年平均风蚀季土壤风蚀量约为6.702 t·hm^-2，月均土壤风蚀量为0~5 t·hm^-2。根据中国土壤侵蚀分类分级标准，将研究区土壤风蚀程度按月尺度分为微度侵蚀（0~0.2 t·hm^-2）、轻度侵蚀（0.2~2.5 t·hm^-2）、中度侵蚀（2.5~5.0 t·hm^-2）、强烈侵蚀（5.0~8.0 t·hm^-2）、极强烈侵蚀（8.0~15.0 t·hm^-2）、剧烈侵蚀（>15.0 t·hm^-2）6个不同等级，研究区春季农田风力侵蚀为微度侵蚀和中度侵蚀，秋季农田风力侵蚀为微度侵蚀。具体表现为春季风蚀量较大，约占风蚀季风蚀量的68.1%，其中，4月风蚀量最大，秋季10—11月、次年冬季1—3月和5月风蚀量较小，整体上与研究区月风速分布特征一致（图7）。玉米残茬在土壤风蚀抑制中起到了重要作用，秋季作物收割后，残茬覆盖度较高，增加了地表粗糙程度，加之冬季1—2月积雪覆盖面积较大，两者都能够有效抵制土壤风蚀现象，因此，秋季11月、次年冬季1—3月土壤风蚀量相对较小。随着残茬衰败和腐烂，其覆盖度逐渐降低，伴随着春耕，更多农田处于裸露状态，残茬覆盖度在春季降至最低值，同时，受风力环境影响，土壤风蚀量在4月达到最大值。

研究区2018—2022年风蚀季土壤风蚀量空间分布与玉米残茬覆盖度空间分布呈现相反特征，玉米残茬覆盖度较低的区域，土壤风蚀量较高，反之，则土壤风蚀量较低，不同月份土壤风蚀量空间分布特征未显示明显差异（图8）。整体上，研究区以微度侵蚀和轻度侵蚀为主，其面积约占研究区耕地面积的54%，集中分布于南部与北部的耕地，强烈侵蚀与极强烈侵蚀区域面积较小，约占研究区耕地面积的18%，主要位于中部和北部边缘。受研究区农田分布特征影响，北部耕地面积较小，中部及两侧主要为裸露沙地，地表粗糙度较小，加剧了土壤风蚀现象。

本研究基于玉米残茬实测高光谱特征，利用3个波段构建了CAI计算玉米残茬覆盖度（FVC_CAI），CAI通过引入可捕捉2 100 nm附近木质素与纤维素吸收特征的短波红外波段，克服二波段指数的不足，有效区分玉米残茬与裸土背景，提升玉米残茬覆盖度的计算精度。此外，使用在相关研究中广泛使用的基于NDVI计算植被覆盖度的方法计算研究区的玉米残茬覆盖度（FVC_NDVI），并与本研究构建方法计算的FVC_CAI进行对比，以进一步评估本研究构建的玉米残茬覆盖度计算方法对玉米残茬覆盖度计算精度的影响。

基于两种方法计算的玉米残茬覆盖度呈现显著差异（图9），整体上，研究区月均FVC_CAI约为19.3%，而月均FVC_NDVI约为13.1%，FVC_CAI月均值普遍高于FVC_NDVI，前者约为后者的1.47倍，反映绿色植被光谱特征的NDVI无法有效表达缺失叶绿素的玉米残茬，明显低估了玉米残茬覆盖度。从时间角度来看，基于两种方法计算的玉米残茬覆盖度呈现较一致的月际变化特征，玉米残茬覆盖度均在10月达到最大值，与玉米收割后残茬集中残留在农田的实际情况相符^［33］，12月到次年1月玉米残茬覆盖度呈现降低趋势。3—4月，气温回升，覆盖在农田地表的冰雪融化，残茬逐渐出露，FVC_CAI出现缓慢增加趋势，但FVC_NDVI未表现出明显的变化趋势，主要原因为NDVI无法识别地表较小的残茬覆盖变化。4—5月，由于农田翻耕及播种，玉米残茬被掩埋或移除，导致覆盖度下降，但5月FVC_NDVI呈现出显著增加趋势，主要原因为该时期绿色植被返青及作物出苗，地表绿色植被覆盖度增大。上述结果表明，CAI能够更精准地反映风蚀季作物残茬覆盖的动态变化，提升土壤风蚀量计算精度^［34］，但在光合植被与非光合植被同时存在的时期，应开发新的方法反映地表的实际覆盖状况。

基于两种指数计算的玉米残茬覆盖度分别输入RWEQ，计算的研究区2018—2022年风蚀季月风蚀量变化存在相似特征，但整体存在显著差异（图10）。总体上，采用NDVI计算的月土壤风蚀量均大于采用CAI计算的月土壤风蚀量，采用NDVI计算的年土壤风蚀量约是采用CAI计算的年土壤风蚀量的0.55倍，二者在4、10、11月的差异较大，其中4月的差异最大，表明如果仅是用NDVI计算研究区玉米残茬覆盖度，将在这些月份普遍高估土壤风蚀量，而在12、1、2、3、5月，两种残茬覆盖度计算方法对土壤风蚀量计算影响较小。

对2018—2022年奈曼旗风蚀季采用NDVI计算的月土壤风蚀量和采用CAI计算的月土壤风蚀量进行差值计算（图11）。采用NDVI进行计算，在90%以上区域系统性高估了土壤风蚀量，采用CAI进行计算仅在研究区中部及北部变化等零星区域轻微高估了土壤风蚀量。采用NDVI计算的月土壤风蚀量和采用CAI计算的月土壤风蚀量的差异在4月最大，相差1.69 t·hm^-2，10、11月次之，二者差异较大的区域仍集中在南部及东北部的耕地集中区域。根据徐斌等^［34］、陈悠等^［35］在奈曼旗实地风沙观测结果，1987—2002年研究区4月平均风蚀速率达到5~6 t·hm^-2。采用NDVI进行计算高估了研究区土壤风蚀量，而采用CAI计算结果与实地观测结果较为接近，表明以FVC_CAI 作为植被覆盖因子输入风蚀模型计算的土壤风蚀量能够更好地反映风蚀季农田风蚀状况。

根据奈曼旗玉米残茬的实测高光谱特征，基于高光谱反射率（2 190、740、443 nm）波段组合反演的纤维素吸收指数（CAI）能够克服二波段指数的不足，有效区分玉米残茬与裸土背景，与玉米残茬覆盖度（FVC）相关性最强，对应Sentinel-2A的B₁₂、B₆、B₁波段构建玉米残茬覆盖度计算方程为FVC=1.344CAI-0.07，R²=0.77，RSME为5.77%。研究区2018—2022年风蚀季玉米残茬覆盖度约为19.3%，4、5月玉米残茬覆盖度最低，基于CAI较基于NDVI计算的玉米残茬覆盖度能更准确反映玉米残茬覆盖的真实状况，基于CAI计算的年均玉米残茬覆盖度是基于NDVI计算的年均玉米残茬覆盖度的1.47倍，后者明显低估了风蚀季玉米残茬覆盖度。研究区多年平均风蚀量约为6.702 t·hm^-2，基于CAI计算的农田土壤风蚀量更接近风蚀观测结果，基于NDVI计算的风蚀季农田土壤风蚀量明显高估。研究区农田风蚀量春季较大，4、5月风蚀量约占风蚀季风蚀量的68.1%，南部与北部风蚀量以微度侵蚀、轻度侵蚀为主，中部与北部边缘以强烈侵蚀、极强烈侵蚀为主。在助力“科尔沁沙地歼灭战”过程中，应积极探索、推进作物残茬覆盖、免（少）耕播种等保护性耕作措施的推广，充分发挥作物残茬在农田土壤风蚀防治中的作用。