土壤风蚀造成的土地退化是全世界共同面临的重要环境问题，全球范围内约28%的土地正遭受土壤风蚀的危害^［1］。中国干旱与半干旱地区的土地退化和沙化现象尤为严重^［2-4］，受土壤风蚀影响的国土面积达到157.61万km²，占国土面积的16.42%^［5］。风蚀使土壤中富含营养元素的细颗粒流失，造成土地生产力下降、河流海洋污染，对农业生产和社会经济活动产生巨大的影响，小粒径土壤颗粒进入空气造成空气质量下降进而影响人体健康^［6-11］。土壤风蚀模型是科学表达土壤风蚀过程、进行土壤风蚀评价和土壤风蚀预报的重要手段，土壤风蚀模型的构建对土壤风蚀危害的科学评估以及土壤风蚀的有效防治具有重要意义。

风是土壤风蚀发生和发展的基本驱动力，无论基于经验统计还是基于过程的土壤风蚀模型，风力侵蚀因子都是土壤风蚀模型中最重要的输入因子^［12-13］，高时间分辨率和高空间分辨率风力数据对土壤风蚀模型中土壤风蚀模数计算具有重要意义^［14-15］。但在实际的风速观测中，因观测设备、观测条件等限制因素，难以获得较大研究范围内高时间分辨率的风速数据^［16-17］。中国气象局公布的中国范围内较长时间序列基本气象资料和产品中，以北京时间（Coordinated Universal Time+8，简称UTC+8）记录的逐日4次定时观测风速数据（UTC+8 02:00、08:00、14:00、20:00）和逐日2类风速数据（平均风速和最大风速）在风蚀研究中应用比较广泛^{［16，18］}。如Lin等^［19］和Du等^［20］用1日2类风速数据分别研究1982—2015年河西走廊的风蚀时空变化以及风蚀驱动力因子和2001—2014年中国北方的土壤风蚀量和粉尘释放量；王旭洋等^［21］和Guo等^［22］利用1日4次风速数据研究中国北方农牧交错带的土壤风蚀的时空变化。此外，逐小时风速、逐日8次风速数据（UTC+8 00:00、03:00、06:00、09:00、12:00、15:00、18:00、21:00）、日平均风速和月平均风速等多种类型的风速数据均被用于风蚀模型的预测研究^［23-26］。

为提高模型的模拟精度，通常采用WINDGEN、线性插值、降尺度等方法来提高风速数据的时间分辨率^{［14，27-30］}，但风速数据的类型不同、时间插值方法的差异都会对模型计算结果造成影响。Skidmore等^［31］和van Donk等^［32］发现在WEPS模型中可以使用1日24 h风速和1日4次风速两种时间分辨率的风速数据计算侵蚀风力密度（WPD），但使用1日24 h风速和1日4次的风速数据估算WPD产生的差异小于10%，这两类风速数据在估算土壤风蚀时造成的误差不大；随后Panebianco等^［17］发现在RWEQ模型中，使用小时尺度的风速数据获得的土壤损失量远低于以5 min为间隔获得的土壤损失量；Guo等^［30］发现以5、10、15、30、60 min等时间尺度统计平均风速计算的平均侵蚀风力密度（AWPD）和风因子（wf）均小于以1 min为统计间隔平均风速的计算结果，并且AWPD和wf随着风速数据统计间隔的增大而减小，不同时间分辨率风速会影响土壤风蚀模型的模拟效果。

空间插值是模型将有限数量观测点的观测值转换为观测指标空间连续分布曲面的必要过程，插值结果的准确程度受观测样点数量、分布、插值方法等多种因素影响。中国2 000多个基准气候站和基本气象站的气象数据在学术研究中应用广泛^［33］，在土壤风蚀研究领域使用风速等气象数据时，除不同学者使用相关数据的时间分辨率存在差异外，由于研究的时间尺度、研究区范围、位置等不同，还存在另外两方面的问题：一方面，不同学者选择的气象观测站类型不同；另一方面，不同学者使用的气象观测站数量存在差异^［20-26］。这些因素会导致土壤风蚀模型计算的区域风蚀模数的差异，进而影响区域土壤风蚀程度的客观评价和土壤风蚀预报的精度。

基于此，本文选取了中国北方农牧交错带130个基准气候站和基本气象站的气象数据，其中风速数据包括1日4次风速数据和1日2类风速数据等两种风速数据，通过设置不同气象站数量、不同风速数据类型及风速数据不同时间插值方法模拟了不同RWEQ模型的不同输入情景，分析该模型模拟结果对于输入参数的敏感性，研究在RWEQ模型中气象站的数量和风速数据类型对估算潜在风蚀模数和区域土壤风蚀状况的影响。

以处于西北干旱半干旱区和东部季风区过渡地带的中国北方农牧交错带为研究区（36°30′—46°42′N、106°16′—124°51′E），面积约为54.3万km²（图1）^［22］。研究区海拔49—3 133 m，地势西高东低；年平均气温为2.4—11.5 ℃；年降水量300—450 mm，季节分布不均，主要发生在夏秋两季。农牧交错带内春季大风和沙尘天气频繁，多年平均风速1.3—3.9 m·s^-1，最大风速可达到24 m·s^-1；农牧交错带自然植被以灌丛、草地和草甸为主，农田、林草地和沙地等多种地类交错分布，砂性土面积分布较为广泛，是中国北方土壤风力侵蚀较严重的地区。该区是京津地区沙尘天气的主要源地和传输路径，对该地区的风蚀程度进行科学评估，有助于京津及华北北部地区生态屏障的构建。

根据RWEQ模型输入数据的要求，收集了2000—2012年中国北方农牧交错带气象数据、土壤数据、植被盖度、数字高程模型（DEM）数据和土地利用类型数据（表1）。气象数据来源于中国气象科学数据共享服务网（http：//data.cma.cn），气象数据中包含了研究区逐日风速、风向、气温、降水、日照时数和湿度等；植被盖度数据来源于USGS网站（https：//www.usgs.gov）MODIS产品的D13A2、OD11A2和MOD09A1的遥感数据产品；土地利用类型数据主要来源于中国科学院资源环境科学数据中心（http：//www.resdc.cn）栅格数据；DEM （Digital Elevation Model）和土壤数据来源于寒区旱区科学数据中心（http：//westdc.westgis.ac.cn）中国1 km分辨率的高程模型数据集及世界土壤数据库中国土壤数据集。所有数据均采用Beijing1954坐标系下的Lambert投影。

修订版的风蚀方程模型（RWEQ）是由美国农业部开发的用于预测田间尺度的土壤风蚀预测模型^［34］。近年来，该模型将预测的范围从田间尺度扩展到区域尺度并得到广泛的应用^{［24-26，35］}。RWEQ模型充分考虑气候条件、地表植被覆盖、地表粗糙度、土壤可蚀性和结皮等因子，可表示为^［36］：

式中：Q_max表示风力的最大输沙能力（kg·m^-1）；WF表示气候因子（kg·m^-1）；EF表示土壤可蚀性因子（无量纲）；SCF表示土壤结皮因子（无量纲）；K′表示土壤粗糙度因子（无量纲）；COG表示结合残茬因子（无量纲）。

气候因子是气象要素对土壤风蚀的综合影响，可根据地面2 m高度风速计算：

式中：SW表示土壤湿度因子（无量纲）；SD表示雪盖因子（无量纲）；wf表示风因子（kg·m^-1·s^-1）；U₂ 表示2 m高度风速（m·s^-1）；U_t 表示2 m高度临界起沙风速（5 m·s^-1）；N_d表示观测风速的时间段，一般为15 d；N表示观测风速时间段内的观测次数；ρ表示空气密度（kg·m^-3）；g表示重力加速度（m·s^-2）。

式中：ET_P表示潜在蒸发量（mm·d^-1）；R表示测定时间内的降水量(mm)；I表示测定时间内的灌溉量(mm)；R_d表示测定时间内的降水日数(d)；N_d表示测定时间内的灌溉日数(d)；SR表示测定时间内的总太阳辐射量(MJ·m^-2·d^-1)；DT表示测定时间的平均气温（℃）；H₀表示大气上空的太阳辐射量（MJ·m^-2·d^-1）；S、S_L表示日照时数和日长(h)。

由于土壤结构和有机质的含量年际变化小，土壤中的黏土、有机质等含量随时间变化缓慢，因此土壤可蚀性因子和结皮因子在相当长的一段时间内不会产生变化^［37-38］。土壤可蚀性因子EF和土壤结皮因子SCF可由公式（7）、（8）表示：

式中：Sa表示土壤砂粒含量（%）；Si表示土壤粉粒含量（%）；Sa/Cl表示土壤砂粒和黏粒含量比（%）；Cl表示土壤黏粒含量（%）；OM表示土壤有机质含量（%）；CaCO₃表示土壤碳酸钙含量（%）。

结合残茬因子COG是平铺残茬SLR_f、直立残茬SLR_s和植被盖度SLR_c的乘积。然而研究区的作物残茬在作物收割后一般会使地表裸露^［22］。因此结合残茬因子可用植被盖度代替：

式中：cc表示植被盖度。

地表粗糙度因子^［24］：

式中：α表示地形坡度。

为研究气象站数量对于模型模拟结果的影响，根据研究区气象站的类型，本文分别选取研究区内130个（基准气候站和基本气象站）、47个（基准气候站）气象站的风速数据作为两类情形。1日4次风速数据和1日2类风速数据在土壤风蚀的研究中应用比较普遍，为对比不同类型风速数据对计算潜在风蚀模数的影响，针对于不同类型和数量的气象站，将风速数据分别做以下处理：1日4次风速数据用线性插值方法生成逐小时风速、1日2类风速数据用WINDGEN生成逐小时风速、1日4次风速数据直接参与模型计算，3类风速数据均采用克里金插值法，插值生成分辨率为1 km的栅格数据。不同类型风速数据与不同数量气象站组合为以下6种模拟情景（表2）。

模型评价包含模拟值与实测值的比较研究、敏感性分析和模型可靠性评价等多个方面的评价。敏感性分析是对模型输入参数与模型响应机制的研究，常用于确定模型主要影响因子^［39-40］。敏感性分析在土壤风蚀、水蚀模型研究中应用广泛，Hagen等^［41］用敏感性分析研究WEPS模型输入参数对模型的影响，Feng等^［40］进一步研究WEPS模型输入参数在土壤风蚀和PM₁₀粉尘释放中的影响。本研究利用敏感性分析法对比RWEQ模型计算的潜在风蚀模数对于不同情景下输入风速数据的响应见表3。敏感系数SS可表示为：

式中：I₂和I₁分别表示输入参数的最大值和最小值；O₂和O₁为对应参数的输出值；I₁₂和O₁₂分别是输入和输出参数的平均值。

风速是RWEQ模型中最敏感的输入参数（表3），土壤风蚀强度与大于临界侵蚀风速的风速数据变化趋势基本一致，风速数据的变化会影响风蚀模数的估算及风蚀强度的空间分布趋势，高时空分辨率的风速数据对区域风蚀模拟尤为重要；RWEQ模型对雪盖因子和植被盖度的变化较敏感，雪盖因子和植被盖度的增加均可有效抑制土壤风蚀；受气候条件的影响，雪盖因子在冬季能有效减缓风蚀，植被盖度在春季、秋季，尤其是夏季对风蚀的抑制效果最明显；降水量和土壤黏粒含量的增加在地表形成土壤结皮可减轻风蚀危害，其他土壤理化性质如砂粒、粉粒、有机质、CaCO₃含量以及地表粗糙度对风蚀的控制效果相比以上参数会减弱；而模型对气温、日照时数等参数最不敏感，这些气候因子在抑制或促进土壤风蚀发生过程中的作用并不明显。

依据中国土壤风蚀分类分级标准，将不同情景下利用RWEQ模型计算的中国北方农牧交错带潜在风蚀模数分为微度侵蚀（0—2 t·hm^-2·a^-1）、轻度侵蚀（2—25 t·hm^-2·a^-1）、中度侵蚀（25—50 t·hm^-2·a^-1）、强烈侵蚀（50—80 t·hm^-2·a^-1）、极强烈侵蚀（80—150 t·hm^-2·a^-1）和剧烈侵蚀（>150 t·hm^-2·a^-1）等6个不同等级^［42］。

在不同模拟情景下，研究区土壤风蚀状况均呈现出明显空间分异特征（图2）。总体上看，研究区风蚀状况以微度侵蚀和轻度侵蚀为主，土壤风蚀模数的高值区集中在毛乌素沙地、浑善达克沙地和科尔沁沙地；而毛乌素沙地以东、浑善达克沙地以南和科尔沁沙地东北方向的区域土壤风力侵蚀危害小，多以微度和轻度侵蚀为主。6种情景下中国北方农牧交错带2000—2012年的多年平均风蚀模数分别为27.33、19.34、22.12、15.73、39.22、28.72 t·hm^-2·a^-1（图3A）。在情景1、情景5和情景6模拟情况下，研究区整体风蚀状况为中度侵蚀，在情景2、情景3和情景4下模拟情况下则为轻度侵蚀。不同风蚀强度面积在所有模拟情景中也存在显著差异，如图3B所示，所有情景中微度侵蚀和轻度侵蚀面积约占研究区总面积的4/5，分别为78.28%、84.79%、80.31%、86.79%、73.62%和78.21%，情景2和情景4下计算的微度侵蚀和轻度侵蚀面积最大；强烈侵蚀、极强烈侵蚀和剧烈侵蚀情况在所有情景中比重较小，分别占总面积的11.63%、9.76%、10.50%、9.39%、19.19%和14.99%，情景5和情景6下强烈侵蚀、极强烈侵蚀和剧烈侵蚀面积最大。从估算的潜在风蚀模数和发生严重侵蚀的面积上看，使用130个站点风速数据模拟的情景1、情景3和情景5分别要大于使用47个站点风速数据模拟的情景2、情景4和情景6；直接使用1日4次未插值风速模拟的情景5和情景6大于其他4个使用线性插值及WINDGEN生成小时风速的模拟情景。

WINDGEN风速生成法被广泛地应用于RWEQ模型和WEPS模型等输入小时风速的土壤风蚀模型，因此，本文选取气象站分布密集的情景1作为标准情景，对比土壤风蚀状况随风速数据的分辨率和气象站数量的变化。就空间分布而言，土壤风蚀热点区集中在三大沙地及周围的农田和草地，但是相比情景1，在情景2中科尔沁沙地和浑善达克沙地的剧烈侵蚀和极强烈侵蚀面积面积减小，情景3和情景4剧烈侵蚀和极强烈侵蚀的面积进一步减小，情景4浑善达克沙地以极强烈侵蚀为主，在情景5和情景6中，浑善达克沙地和科尔沁沙地的剧烈侵蚀和极强烈侵蚀的面积扩大，毛乌素沙地土壤侵蚀危害增加，而且随气象站数量的增加，剧烈和极强烈土壤侵蚀面积增加。

RWEQ模型中潜在土壤风蚀模数以半月时间尺度计算，本研究将半月尺度计算结果按照季节尺度（3、4、5月划分为春季，6、7、8月划分为夏季，9、10、11月划分为秋季，12、1、2月划分为冬季）进行统计（图3C）。在土壤风蚀模数的季节变化上，北方农牧交错带土壤风蚀主要发生在春季，夏季土壤风力侵蚀危害最小；在6种模拟情景下，春季潜在土壤风蚀模数分别占全年的65.05%、70.96%、94.77%、70.79%、61.24%、59.25%，情景1、情景5和情景6下估算的春季潜在土壤风蚀模数占全年总风蚀模数比重较小，而情景3下估算的春季潜在土壤风蚀模数占全年风蚀模数的比例最高；冬季估算的风蚀模数仅次于春季，夏季风蚀模数最小，分别为0.60、0.42、0.51、0.36、1.12、0.81 t·hm^-2·a^-1。从土壤风蚀模数的年际变化看，6种情景中2000—2012年潜在土壤风蚀模数均呈明显波动下降趋势，研究时间段内北方农牧交错带土壤风蚀状况有明显的改善（图3D），但不同模拟情景下估算的年均潜在风蚀模数变率存在较大差异，除情景5，最大土壤风蚀模数出现在2000年，为74.44 t·hm^-2·a^-1，其他情景估算的最大土壤风蚀模数出现在2001年，分别为52.83、41.42、43.98、35.08、55.73 t·hm^-2·a^-1；最小的潜在土壤风蚀模数均出现在2011年，分别为8.82、7.27、7.03、5.68、14.46、12.04 t·hm^-2·a^-1。与标准情景1相比，情景2、情景3和情景4的年均土壤风蚀模数和各季节土壤风蚀模数均小于情景1，而情景5和情景6估算的土壤风蚀模数大于情景1。6种不同模拟情景估算的年均土壤风蚀模数和不同季节的土壤风蚀模数基本呈现情景5>情景6>情景1>情景3>情景2>情景4的关系。

风是土壤风蚀发生的基本驱动力，在RWEQ模型诸多输入参数中，风蚀模数的计算结果对风速数据的输入最敏感，高质量的风速数据能提高模型的模拟效果^{［29-30，41］}。在估算潜在风蚀量和区域土壤风蚀状况时，选取的气象站点的数量和风速数据的分辨率都会对模拟结果产生影响。在风速数据相同的条件下，130个气象站的气象数据估算的区域多年平均潜在土壤风蚀模数均大于47个气象站估算的潜在风蚀模数；不同分辨率的风速数据也会对潜在风蚀量和区域土壤风蚀状况的估算产生影响；随着风速数据分辨率的提高，RWEQ模型估算的潜在风蚀模数会提升。风速数据的时间插值方法也会对风蚀模拟的结果产生重要影响，根据本文设计的不同模拟情形，分别使用WINDGEN、线性插值将1日2类风速和1日4次风速数据生成小时风速并结合1日4次风速数据估算了2000—2012年的月均风因子（wf）和2 m高度大于临界风速（5 m·s^-1）出现频率（图4A）。6种模拟情景中，风因子和大于临界风速的风速数频率变化趋势一致，2000—2012年大于临界风速数在总风速数中所占百分比分别是5.77%、4.88%、4.89%、4.09%、6.44%和5.61%，估算的风因子分别为2.03、1.51、1.56、1.18、2.73 kg·m^-1·s^-1和2.06 kg·m^-1·s^-1。使用1日4次风速数据的情景5和情景6估算的结果明显大于使用WINDGEN生成小时风速的情景1和情景2，使用线性插值生成小时数据的情景3和情景4估算结果最小，且6种模拟情景的风因子和大于临界风速数频率与多年平均土壤风蚀模数的变化一致。

不同模拟情景下，2000—2012年各月风速大于临界风速数平均百分比的年际变化趋势见图4B，大于临界风速数的月均百分比与年均土壤风蚀模数的变化趋势一致；各模拟情景下，大于临界风速数的月均百分比2001年最大，2011年最小，且2001年到2003年有显著的下降；情景2、情景3和情景4的所占百分比明显低于情景1、情景5和情景6，以上结果与风速数据的时间插值方法相关。RWEQ模型和WEPS模型在计算中一般使用小时风速，WINDGEN风速生成法广泛地应用到RWEQ模型的模拟中，WINDGEN能较好地模拟每日小时风速数据的变化趋势，并在5—8 m·s^-1会较高地模拟风速^［14］；而线性插值拟合的小时风速会降低大于临界风速数出现的频率，使RWEQ模型估算的研究区潜在土壤风蚀模数偏小，也会降低模型模拟精度。

风速在RWEQ模型中是最敏感的输入参数，风速数据的变化会影响风蚀模数的估算及风蚀强度的空间分布。在估算区域潜在风蚀量时选取的气象站点的数量和风速数据类型都会对模拟结果产生影响。RWEQ模型估算的区域潜在风蚀模数会随气象站点的数量和风速分辨率的提升而增加；此外，风速数据的线性插值方法在RWEQ模型中应用效果不理想，与WEPS模型的WINDGEN风速生成方法相比，线性插值法使地面2 m高度大于临界起沙风速（5 m·s^-1）的风速频率降低，过低估算区域潜在土壤风蚀模数和土壤风蚀状况。此外，风速数据的空间插值方法对模型计算的潜在风蚀模数也存在影响，本文未做进一步讨论。

RWEQ模型是基于田间尺度构建的土壤风蚀模型，将该模型从田间尺度扩展到大范围的土壤风蚀研究，应对模型参数和各因子进行适用于区域研究状况的调整^［23-25］，并结合大量实测数据对模型进行修正。在风蚀模拟研究中，应尽可能输入高时空分辨率的模型参数，选择合理的时间和空间插值方法，以获取精度及时空分辨率较高的土壤风蚀模数数据，服务于区域土壤风蚀防治。