0 引言

干旱是指水分收支或供求不平衡所造成的水分短缺现象^[1].据统计，气象灾害引起的损失占各类自然灾害的85%，而干旱又占其损失的50%^[2].独特的地理环境与气候特点导致我国易受干旱影响，且易造成损失^[3].干旱指数作为可以表征干旱发生的重要指标，能够将干旱的严重性进行等级区分和量化.目前为止，已经发展了多种干旱监测指数，如基于气象站点数据计算的标准化降水指数（Standardized Precipitation Index，SPI）、标准化降水蒸散指数（Standardized Precipitation Evapotranspiration Index，SPEI）和帕默尔干旱指数（Palmer Drought Severity Index，PDSI）^[4-6]等，但由于站点数量和位置分布不均，干旱监测缺乏空间连续性.随着遥感技术的发展，不同类型传感器可获取大范围近实时的降水、植被和土壤湿度等信息，为区域干旱监测提供了方法，提高了干旱监测的空间连续性.一些遥感干旱指数，例如温度状况指数（Temperature Condition Index，TCI）、降水状况指数（Precipitation Condition Index，PCI）、植被状况指数（Vegetation Condition Index，VCI）、温度植被干旱指数（Temperature Vegetation Dryness Index，TVDI）和条件温度植被指数（Vegetation Temperature Condition Index，TVCI）^[7-9]等广泛应用于大范围干旱监测，但遥感数据往往因其时间分辨率低影响干旱监测效果.

由于干旱与降水、地形、潜在蒸散发和土壤湿度等多种变量有着密切的关系，导致干旱成因复杂.虽然单一影响因子或干旱监测指数在有关行业和领域有着较好的应用，但会出现对于干旱因子的影响考虑不足或受限于应用条件，往往无法充分反映区域内的综合干旱发生状况.为此，许多学者利用多种数据及多种方法，开展了综合干旱监测研究^[10-12].目前，综合干旱指数主要通过权重组合、联合分布和机器学习^[13]等3类方法进行构建.权重组合法由于指数的选择以及权重的分配限制了指数的应用.联合分布方法在需要综合的干旱因子或指数较多时，不同因子或指数之间关系的复杂性会加大多元分布估计的难度，影响监测效果.机器学习方法因其具有处理非线性问题的能力，能够对多种干旱因子进行信息挖掘，在综合干旱监测上被广泛应用.然而，由于众多不同来源的遥感数据精度不同，单纯依靠遥感数据构建的综合干旱监测指数对于空间连续的干旱监测还存在不足.考虑到综合干旱监测在日常农业生产和抢险救灾等方面的重要性，需要更合适的方法和数据以构建日尺度的综合干旱监测模型.

中国气象局陆面数据同化系统（CMA Land Data Assimilation System，CLDAS）是我国自主建立的陆面数据同化系统，其通过融合多种气象站点观测、卫星遥感观测和再分析资料，输出多种高质量的气象及陆面数据，为日尺度气象监测提供重要的数据支撑^[14-16].然而，当前缺少将CLDAS数据应用于综合干旱监测的研究，尚不清楚CLDAS应用于综合干旱监测的效果.因此，本文通过引入2015—2018年连续日尺度CLDAS数据，再辅以相应的中等分辨率成像光谱仪（Moderate-resolution Imaging Spectroradiometer，MODIS）数据，通过梯度提升机（Gradient Boosting Machine，GBM）方法构建日尺度综合干旱监测模型，统筹考虑与干旱发生的多种相关因子，以获得区域日尺度下干旱发生的整体状况，通过评估模型对于日尺度综合干旱的监测效果，以期获得更加精细合理的干旱监测结果.

1 研究区与研究数据

1.1 研究区

选取华北地区为研究区（图1）^[17]，经纬度范围为110°~120°E，34°~43°N，包括山西省、河北省、北京市、天津市以及内蒙古自治区中部部分地区.研究区南部临近黄河，北部为内蒙古高原，东部有海河平原，西部为吕梁山脉，地势总体西高东低，下垫面情况复杂.华北地区是我国重要的农业主产区和粮食生产基地，主要种植冬小麦和玉米，研究区域大部分属暖温带半湿润半干旱季风气候区，年降雨量大约为400~500 mm，雨季主要集中在夏季，雨热同期，干湿期分隔明显.在全球气候变化的背景下，该区域年际降水分配不均，气候多变，干旱灾害事件频繁发生^[18].

1.2 研究数据

1.2.1 MODIS数据

反射率数据来自于美国国家航空航天局（https://ladsweb.modaps.eosdis.nasa.gov）的Terra MODIS表面反射率产品MOD09GA，该产品为表面波谱反射估计，数据经过大气和气溶胶的校正，空间分辨率为500 m，时间分辨率为1 d，提供1~7波段每日栅格化的数据.本文所用到的波段为红波段（620~670 nm）和近红外波段（841~876 nm）反射率数据，研究时段为2015—2018年.将地表反射率数据进行投影一致性转换、研究区的矢量边界裁剪等预处理后，通过双线性内插法重采样至0.062 5°× 0.062 5°分辨率，使其与CLDAS空间分辨率一致.

1.2.2 CLDAS数据

CLDAS大气驱动场产品V2.0覆盖亚洲区域（60°~160°E，0°~65°N），空间分辨率为0.062 5°×0.062 5°，时间分辨率为逐小时，包括2 m气温、2 m比湿、10 m风速、地面气压、降水和短波辐射等6个要素的等经纬网格融合分析产品.CLDAS土壤体积含水量分析产品V2.0覆盖区域和时空分辨率与大气驱动场产品一致，垂直分为5层，分别为[0 cm，5 cm]、[0 cm，10 cm）、[10 cm，40 cm）、[40 cm，100 cm）、[100 cm，200 cm].CLDAS地表温度数据分析产品V2.0与CLDAS大气驱动和土壤体积含水量产品覆盖区域与时空分辨率相同^[15，19].本研究使用2015—2018年大气驱动产品中的地面气压、短波辐射、2 m气温、降水要素，以及0~10 cm土壤体积含水量和地表温度数据分析产品，并将逐小时数据通过求平均或者加和合成为日尺度数据.

1.2.3 气象站点数据

气象站资料来源于国家气象信息中心2015—2018年中国地面气象站逐小时的气温和降水数据（http://data.cma.cn/），选取研究区域内158个气象站点和68个农业气象站点，将逐小时降水加和、逐小时的气温求取平均值，分别得到研究区逐日的降水数据和气温数据.将通过气象站点数据计算得出的CI监测数据，75%作为训练集，25%作为测试集，将其作为综合干旱监测模型的建模和验证输入数据.

1.2.4 其他辅助数据

1）DEM数据.数字高程模型（Digital Elevation Model，DEM）为美国“奋进号”航天飞机的雷达地形测绘（Shuttle Radar Topography Mission，SRTM）获得的DEM数据，来源于中国科学院资源环境科学数据中心（http://www.resdc.cn），范围为60°N~60°S.采用基于4.1版本的DEM数据，利用新的插值算法得到，空间分辨率为1 km.将高程数据进行投影一致性转换、研究区的矢量边界裁剪等预处理后，通过双线性内插法重采样至0.062 5°×0.062 5°分辨率，使其与CLDAS空间分辨率一致.

2）土壤水文参数数据.土壤水分亏缺指数（Soil Water Deficit Index，SWDI）的计算需要利用土壤水文参数数据得出.本文所用的土壤水文参数数据来源于北京师范大学发布的面向陆面过程模型的土壤水文数据集^[20]，分辨率为30角秒.利用此数据的-33 kPa和-1 500 kPa土壤水势数据，作为田间持水量数据（Field Capacity，FC）和土壤凋萎湿度（Wilting Point，WP）.土壤数据垂直剖面分为7层，分别为[0 m，0.045 m）、[0.045 m，0.091 m）、[0.091 m，0.166m）、[0.166 m，0.289 m）、[0.289 m，0.493 m）、[0.493 m，0.829 m）、[0.829 m，1.383 m]，最深为1.383 m.研究首先利用双线性内插法将其重采样至0.062 5°×0.062 5°，然后再插值到研究需要的0~10 cm土层深度，以田间持水量（FC）数据为例，具体计算方法如下：

式中:C_F为重新插值后（0~10 cm）土壤深度的田间持水量（cm³·cm^-3）; C_F，i为0~10 cm内第i层田间持水量平均（cm³·cm^-3）; Z_i为第i层田间持水量的土层深度（cm）; Z为所需的田间持水量土层深度即10 cm.

2 研究方法

2.1 综合干旱监测模型的构建

干旱受到降水、地表温度（或气温）、土壤、地表植被状况和地形等关键因素的综合影响.干旱主要是由于大气水分环流异常，导致降水时空分布不均，从而出现降水不足造成的.本研究利用CLDAS降水数据（PRE）反映研究区的降水异常状况，同时通过降水与蒸散的差值（AWD）反映地表水分的亏缺状况.土壤湿度则会受到降水不足和土壤性质的影响，土壤水分亏缺引起干旱.本研究选用CLDAS土壤湿度（SM）反映土壤水分不足的状况^[21]，同时利用土壤水分亏缺指数（SWDI）反映不同土壤类型下土壤湿度对干旱的响应.当干旱发生时，植物受到水分胁迫影响其生长过程，通过CLDAS地表温度（Land Surface Temperature，LST）反映地表温度对植被生长的制约情况，植被供水指数（Vegetation Supplication Water Index，VSWI）反映植被生长受到水分胁迫的严重程度.海拔（DEM）以及坡度（slope）所导致的下垫面状况差异对于区域干旱的发展同样有重要的影响^[22].为此，利用梯度提升机（GBM），将能够综合反映降水量与潜在蒸散发的综合气象干旱指数值（CI）作为因变量，能够反映干旱状况的多种因子和指数作为自变量，构建以CI=f（PRE，SM，LST，AWD，SWDI，VSWI，DEM，slope）为形式的综合干旱监测模型，技术路线如图2所示.

2.1.1 干旱指数

1）大气水分亏缺指数（Atmospheric Water Deficit，AWD）

降水与蒸散发都是影响干旱发生的重要因子，尤其是降水量的减少，往往会直接引起气象干旱的发生.大气水分亏缺指数（AWD）^[23]为降水与潜在蒸散发的差值，选用其作为构建综合干旱指数的因子来反映大气的水分状况，具体计算方法^[24-25]如下：

式中:i_AWD，i为某天的大气水分亏缺量; P_i为降雨量; m_PET，i为潜在蒸散量.当AWD＜0时，说明大气出现了水分亏缺的状况.

2）土壤水分亏缺指数（Soil Water Deficit Index，SWDI）

土壤干湿状况是反映农业干旱发生的重要指标，而土壤水分亏缺指数（SWDI）^[26]能够表征土壤水分缺失的严重程度，选取其作为综合干旱指数构建的一部分，来表征土壤状况，具体计算方法如下：

式中:θ是土壤含水量; θ_FC表示田间持水量; θ_AWC表示土壤有效含水量; θ_WP表示土壤凋萎湿度; θ_AWC是θ_FC与θ_WP的差值，即θ_AWC=θ_FC-θ_WP.

3）植被供水指数（Vegetation Supplication Water Index，VSWI）

植被供水指数（VSWI）通过植被生长受到制约的状况反映干旱程度，将其作为构建综合干旱指数的因子^[27]，来表征干旱发生与植被变化的响应关系，VSWI值越小，说明干旱发生状况越严重，VSWI值越大，说明干旱发生程度轻.具体计算方法如下：

式中:T_LS为地表温度; i_NDV为归一化植被指数.日尺度NDVI的计算使用MOD09GA产品的红波段（620~670 nm）和近红外波段（841~876 nm），计算公式如下：

式中:R_NI表示近红外波段反射率; R表示红波段反射率.

4）综合气象干旱指数（CI）

综合气象干旱指数（CI）广泛应用于气象干旱监测，在我国具有广泛的适用性.综合气象干旱指数是通过近30天和近60天的标准化降水指数（SPI），以及近30天的相对湿润度指数（MI）综合计算得到的，其具体计算公式^[28]如下：

式中：a为近30天标准化降水系数，取0.4；b为近60天标准化降水系数，取0.4；c为近30天相对湿润系数，取0.8；I_SP，30为近30天的标准化降水指数值；I_SP，60为近60天的标准化降水指数值；I_M，30为近30天的相对湿润度指数.CI指数的划分标准^[28]如表1所示.

2.1.2 梯度提升机（GBM）

在机器学习中，提升方法（Boosting）是一种通过组合一群复杂程度和训练成本程度低、不容易过拟合的弱分类器（Weak Learner），建立N个模型，并尝试在每次分类中都将上一次分错的数据权重提高一点再进行分类，来获得一个强分类器（Strong Learner）的方法^[21，29].梯度提升机（GBM）属于Boosting的一种，是主要用于解决分类和回归问题的机器学习算法，通过集成多个弱学习器，每个弱学习拟合先前累加模型的损失函数的负梯度，使加上该弱学习器后的累积模型损失往负梯度的方向减少，通过优化损失函数来生成最终强预测模型.梯度提升机机器学习算法对输入变量的相关性没有严格要求，其每轮训练目标是使预测值与真实值的平方误差最小.在利用梯度提升机进行模型构建时，考虑到可能出现的过拟合现象，所以在进行模型构建时采用交叉验证的方法来尽可能避免模型出现过拟合.

本研究首先将基于CLDAS与MODIS数据日尺度的干旱因子和指数值提取至站点并进行归一化后，与基于站点的日尺度CI指数共同作为输入数据，建立梯度提升机回归模型.之后，将基于CLDAS和MODIS数据的干旱因子和指数输入建立的回归模型中，生成对研究区内综合干旱监测模型的预测值，用于干旱监测.

2.2 模型的验证及相关性分析

参照综合气象干旱指数（CI）的旱情等级划分标准，将模型预测得到的干旱监测指数值分为特旱、重旱、中旱、轻旱和无旱5个等级，据此，将CI指数与模型监测的干旱等级进行对比，计算总体一致率，具体计算方法如下：

为了了解构建的综合干旱监测模型的监测效果，研究将其与标准化降水蒸散指数（Standardized Precipitation Evapotranspiration Index，SPEI）和土壤相对湿度（Relative Soil Moisture，RSM）进行相关性分析，具体计算方法如下：

式中:x和y分别为进行相关性分析的两个因子; R_xy为因子相关系数; x_i和y_i分别为两个因子的第i个值; $\stackrel{-}{x}$和$\stackrel{-}{y}$分别为两个因子的平均值.

3 结果与分析

3.1 综合干旱监测模型验证

表2所示75%训练集综合干旱监测值与实测综合气象干旱指数（CI）间的决定系数（R²）达到0.945，25%测试集综合干旱监测值与实测综合气象干旱指数（CI）的R²为0.655，两者有着较好的相关性，测试集的均方根误差（RMSE）和平均绝对误差（MAE）分别为0.082、0.062，误差较小，说明所构建的综合干旱监测模型具有较高的精度.

综合干旱监测模型与CI指数监测的干旱等级一致率比较结果表明（表3），模型监测值与CI值的等级一致率达到65.1%~89.4%，7、8、9月等级一致率较高，在82%以上，其余月份等级一致率也都在65%以上，总体上两者监测结果表现出较好的一致率.

3.2 与标准化降水蒸散指数（SPEI）的相关性

标准化降水蒸散指数（SPEI）是重要的气象干旱监测指数，可采用它来验证构建的综合干旱监测模型对气象干旱的监测能力^[30].2015—2018年日尺度SPEI指数与综合干旱监测模型监测结果的相关分析表明，两者具有较好的相关性（图3），各月相关系数在0.634~0.735之间，且各个月份间相关系数相差较小，其中4月和7月相对较高，分别为0.735和0.734，显著性检验p值均小于0.01，相关性结果都达到极显著水平，说明综合干旱监测模型能够较好地反映气象干旱的发生，具有较好的气象干旱监测能力.

3.3 与土壤相对湿度（RSM）指标的相关性

土壤相对湿度是土壤含水量与田间持水量的百分比，是反映土壤水分平衡的重要指标，通常用来反映干旱对于农业的影响^[31].研究分析了综合干旱监测模型与0~10 cm深度土壤相对湿度的相关性（图4），4月和5月的相关系数分别为0.574和0.51，这可能与春季华北地区春小麦播种时农田灌溉，对土壤湿度状况产生影响有关，其余月份相关系数均大于0.6，显著性检验p值均小于0.01，相关性结果都达到极显著水平，模型能够较好地反映土壤相对湿度变化对于干旱的发生的影响.

3.4 典型干旱情况监测分析

华北地区春季升温较快，当蒸发量较大，且尚未进入雨季时，极易导致春旱发生.我国2015—2018年均有春旱发生，尤其2016年北方冬麦区春旱^[32]，3月中旬，山西、河北等省部分地区开始发生旱情，4月上旬，旱情快速发展，作物因旱大面积受损.本文选取4月作为研究月份，对综合干旱监测模型与利用克里金进行CI插值的监测结果进行了比较，计算出每月发生干旱的天数.

图5a结果显示，2015年4月华北地区干旱天数较多，如河北省西北部、山西中部和北部，以及内蒙古中部地区，干旱发生都超过了15 d，而山西南部、河北南部干旱天数相对较少，低于10 d.2016年4月华北地区干旱发生如图5b所示，除山西南部部分地区外，华北大部干旱发生天数大于15 d.如图5c所示，2017华北地区4月干旱发生较多地区主要在河北中北部以及内蒙古中部地区，干旱发生天数约在20 d以上，而河北南部、山西中南部干旱发生相对较少.如图5d所示，2018年4月干旱发生较多的地区主要集中在内蒙古中部、山西大部、河北少部分地区，天数约在15 d以上.

此外，相较于模型监测结果，CI指数通过插值由区域尺度转化为空间尺度时，较少考虑地表温度、地表土壤或植被等下垫面分布状况对干旱发生的影响.CI监测结果（图5f）显示山西中部地区2016年4月干旱发生小于10 d，与综合干旱监测模型结果相比较低（图5b），可能是由于此处的吕梁山脉黄土广布同时植被覆盖度较低，产生的降水难以下渗或被截留，水土流失严重，所以干旱发生状况较为严重，CI监测结果低估.2018年4月模型监测结果（图5d）显示河北西北部与中南部干旱状况明显，这是由于河北地形的差异，河北西北部山区高原与谷地错落分布，以及林地广布使得下垫面状况较为复杂，而中南大部为平原地区，地形平坦且多为粮食种植区.但由于CI监测对于下垫面差异的因素考虑较少（图5h），所以CI监测结果显示河北2018年4月干旱天数差异较少，未能充分反映河北省西北部和中南部干旱分布差异的细节.

本文使用克里金插值将由站点计算的CI值从区域尺度转换为空间尺度（图5e—h）的结果表明，虽然站点插值结果基本反映了干旱发展过程，然而由于站点数量有限同时站点分布不均，结果无法有效反映部分干旱状况，所以其结果在没有气象站点分布的地方缺乏可信度.例如内蒙古中部地区2015— 2018年存在较严重的干旱状况，干旱发生天数在20 d以上，但CI站点监测结果则显示此处干旱状况偏轻，干旱严重程度低于CI监测结果，这是由于插值方法为单纯的数学过程，无法顾及到站点数量较少或位置分布不均所导致的监测结果不准确的问题，从而影响监测效果.

4 结论与讨论

本文通过综合考虑能够影响或反映干旱的影响因子，充分利用MODIS卫星遥感资料和CLDAS数据的优势，选取相应降水、地表温度、土壤湿度、大气水分亏缺指数、植被供水指数、土壤水分亏缺指数、高程和坡度等监测指标作为自变量，以综合气象干旱指数（CI）作为因变量，利用梯度提升机（GBM）方法，构建了综合干旱日尺度监测模型，并进行了验证和相关性分析，且通过2015—2018年华北地区的干旱发生情况，对模型的适用性进行分析，结果表明：

1）综合干旱监测模型监测值与CI值的训练集的相关决定系数（R²）为0.945，与测试集R²为0.665，有较好的相关性，训练集均方根误差（RMSE）与平均绝对误差（MAE）分别为0.033与0.021，测试集RMSE与MAE分别为0.082与0.062，误差较小.生长季（4—10月）模型监测结果与CI干旱监测等级一致率均在65%以上，说明模型预测有较高的精度.

2）2015—2018年生长季（4—10月）逐天的综合干旱监测模型监测结果与基于站点的标准化降水蒸散指数（SPEI）和土壤相对湿度（RSM）相关系数也大多在0.6以上，呈现较好的相关性.说明模型不仅对气象干旱具有一定的监测能力，同时对于土壤相对湿度的变化及区域农业干旱也有一定的监测能力.

3）典型干旱情况监测对比结果表明，与站点CI监测方法相比，综合干旱监测模型能够较为准确地识别出干旱情况.同时，模型监测结果能够综合考虑各种干旱影响因素，如下垫面地形地势以及植被生长状况等，较好地表征综合干旱状况且不受站点数量和位置分布的影响.

此外，本文综合干旱监测模型仍存在一定的局限性，需要在今后的研究工作中改进：

1）本文主要利用机器学习中梯度提升机（GBM）方法从MODIS数据以及CLDAS数据中提取干旱发生的相关信息，探索利用多源数据进行区域的综合干旱监测方法，但对于更多其他构建综合干旱监测模型的方法，还需进一步研究.另外，本文提出的综合干旱监测模型，在构建时仍依赖研究区站点观测数据，这会影响模型的建立和模型的应用和推广.因此，如何不受站点观测数据限制建立综合干旱监测模型，仍需要进一步研究.

2）本文从大气、土壤和植被多方面综合考虑多种影响干旱因子，进行综合干旱监测模型的构建，并利用CLDAS数据，对日尺度综合干旱监测进行了初步探索，并取得了较好的效果.但是由于干旱成因复杂、致旱因子复杂多样，如植被类型、地表反照率或人类活动等都会对干旱发生产生影响，另外，由于所选择的降水、植被等自变量的前期状况也会对长历时或现时干旱的发生产生影响，而本文并未考虑到降水、植被等变量影响的滞后效应对于长历时或现时干旱的影响.因此，紧密结合区域特性，更精细地探究不同致旱因子在干旱发生中的作用，引入更长时间序列的CLDAS精确数据，辅以更加优越精确的干旱监测模型，例如处理延迟序列具有优势的长短期记忆网络（LSTM）模型，将滞后效应的影响纳入分析，进行综合干旱监测的研究，这些仍是需要进一步探讨的重要问题.

参考文献