生态系统服务的概念为人从生态系统中获得的各种惠益^[1]，在联合国千年生态系统评估中生态系统服务被分为供给、调节、支持、文化4类，且随着人类活动的扩张加剧，24种主要生态系统服务类型中的15种已发生了较大程度的退化^[2]。应用模型法对于生态系统服务进行量化评价与空间化展示，可以明确研究区域的生态现状与亟需提升的区域分布，提升城市可持续发展的能力^[3]。常用模型包括InVEST模型、SWAT模型以及ARIES模型等，如InVEST模型中的产水量模块与土壤保持模块在自然社会发展现状不同的区域进行了广泛的应用，进行水供给服务与土壤保持服务的研究^[4-6]；SWAT模型这一分布式的水文模型被应用于水量调节服务的研究中^[7]；ARIES模型可通过生态系统服务流动的时空动态分析，应用子模块SPAN(Service Path Attribution Network)对“源”“汇”和“使用人”进行分析^[8-9]。生态系统服务簇是共同出现较为频繁的生态系统服务集^[10]，可通过聚类的方法实现区域的划定。相关研究集中于不同类型生态系统服务的选取^[11-12]、多年间生态系统服务簇的转移特征^[13]以及不同生态系统服务簇内各生态系统服务类型的驱动因子等^[12]，针对生态系统服务簇的研究有利于以不同尺度单元进行研究区域的管理与规划。

汉江流域是南水北调中线工程的水源地，同时也是重要的生态屏障区，该区域良好的水资源环境现状是实现区域可持续发展的重要基础。自20世纪以来，汉江流域的生态环境都是国家及省市各级部门建设与改善重点，退耕还林还草、林地保护以及水土保持等生态建设工程提升了汉江流域的植被覆盖度，扩大了水系面积，从整体上改善了水循环的效率^[14]，但是随着气候变化以及人类活动的加剧，汉江流域降雨量出现了汛期与非汛期不平衡的现状，水资源总量有下降的趋势^[15]。其中汉江流域中下游区域城镇化的发展以及灌溉用水需求的增加较为显著，水生态环境发展不平衡的问题日益突出^[14]。汉江流域水生态及生态系统服务的研究主要集中在以下方面：针对生态系统服务的量化与演变，胡砚霞^[14]、陈泽怡^[16]以水生态相关的水源涵养与水质净化为重点，进行多年间生态系统服务时空变化的分析，明确了各类生态系统服务的变化特征；针对生态系统服务间的权衡与协同，高艳丽^[17]、王鹏涛等^[18]对土壤保持服务、产水服务、食物供给服务以及固碳服务等进行了权衡与协同分析，明确了生态系统服务权衡区域与协同区域在汉江流域的分布，并为区域发展与生态保护策略提供依据；针对不同环境变量对汉江流域水环境的影响，森林植被^[19]、土地利用^[20]以及气象条件等^[21]因素都被考虑在内，研究证明森林植被、降雨对水涵养量以及水质均存在积极的影响，而不同的土地利用类型对水质的影响程度不同；针对不同发展情景下的水环境发展前景，Hong等^[22]针对设置的现状、短期规划以及长期规划进行分析，明确了汉江中下游水供给服务易损性以及脆弱性较高的区域。

随着国土空间规划的实施，在湖北省国土空间规划(2021—2035年)的生态修复重大工程中，汉江流域也被纳入到水生态综合修复与水土保持综合修复的建设中，生态系统服务功能中的水源涵养服务功能与土壤保持服务功能可在量化层面为该区域的研究提供数据支撑。其中前者为生态系统提供的调节径流、拦蓄降雨以及影响降雨量的能力，实质为生态系统对于降雨的存储能力；后者为生态系统提供的减小土壤侵蚀流失与保持储存泥沙的能力。目前汉江流域缺乏以子流域尺度将生态系统服务应用于空间布局与规划的研究，在发展不平衡的汉江中下游区域开展相关生态系统服务的研究对于明确水源涵养的重点区域、实现水资源的可持续发展有重要意义。本研究首先应用InVEST模型中的产水量模块与土壤保持模块进行多年水源涵养量以及土壤保持强度的研究；其次基于多年降雨量、水源涵养量与土壤保持强度的特征应用聚类的方法在子流域尺度下进行不同生态分区的划定，明确在气候变化背景下的生态分区转移特征；在此基础上选取环境变量，并明确不同生态分区各环境变量与生态系统服务的关联特征。本研究将为未来汉江流域中下游区域的水生态环境可持续发展提供科学依据，在发展的基础上实现生态环境的保护。

1 材料与方法 1.1 研究区概况

汉江流域中下游(30°4′—32°50′N，110°16′—114°17′E)位于湖北省北部，毗邻河南省，连通十堰、神农架、襄阳、荆门等10个行政区，占地面积3.68×10⁴ km²。该区域地貌以平原、山地、丘陵为主，包含中部地区最大的平原江汉平原以及部分神农架林地的山地，地势较高区域位于研究区域西侧。汉江是南水北调中线工程的水源地，水资源相对比较丰富，气候类型以亚热带季风气候为主，湖北省水利厅数据表明2000—2020年该区域降雨量约为721 mm，且呈上升的趋势。综合考虑降雨量变化的周期性特征与土地利用变化的线性特征，本研究考虑到2010年后数据较为完善，选取2010，2015和2020年与整体降雨量变化特征较为一致的3 a作为研究的时间节点。

1.2 数据来源

本研究的基础数据包括流域、气象、DEM、土壤、土地利用等。其中流域数据来源于国家冰川冻土沙漠科学数据中心(http://www.ncdc.ac.cn/portal/)，为全国1∶25万流域数据集；气象数据来源于国家气象科学数据中心(http://data.cma.cn)，包括研究区域周边50个气象站点2000，2015和2020年3 a的降雨、气温、湿度、风速、日照时长月值数据；DEM数据来源于地理空间数据云(http://www.gscloud.cn/search)，为ASTER GDEM 30 m分辨率的高程数据；土壤数据来源于国家青藏高原科学数据中心(http://www.tpdc.ac.cn/zh-hans/)，为基于HWSD(v1.2)的土壤数据集；土地利用数据来源于地理国情监测云平台，包括2010，2015和2020年这3 a的30 m栅格土地利用数据。

1.3 研究方法 1.3.1 水源涵养服务

水源涵养服务以产水量(Y)、地形指数(TI)、土壤饱和导水率(K)、径流流速系数(V)为基础，综合计算得到水源涵养量^[23-24]。其中InVEST模型中的产水量模块被广泛应用于产水量的计算。水源涵养量的计算方法为：

$ Y(x)=\left[1-\frac{\operatorname{AET}(x)}{P(x)}\right] \cdot P(x) $

(1)

式中：AET(x)表示栅格单元x的年实际蒸散量(mm)，通过Penman-Monteith方程计算^[25]；P(x)表示栅格单元x的年降水量(mm)，蒸散量与降水量通过克里金插值法得到。

$ \frac{\operatorname{AET}(x)}{P(x)}=1+\frac{\operatorname{PET}(x)}{P(x)}-\left\{1+\left[\frac{\operatorname{PET}(x)}{P(x)}\right]^\omega\right\}^{1 / \omega} $

(2)

式中：AET(x)/P(x)即Budyko曲线的近似值^[26]；PET(x)表示潜在蒸发量(mm)；ω为经验参数。

$ \operatorname{PET}(x)=K_c\left(l_x\right) \cdot \operatorname{ET}_0(x) $

(3)

式中：K_c(l_x)为栅格单元x中不同土地覆被类型的植被蒸散系数；ET₀为参考作物蒸散系数。

$ \omega(x)=Z \frac{A W C(x)}{P(x)}+1.25 $

(4)

$ \mathrm{TI}=\log _2\left(\frac{D}{S} \cdot P\right) $

(5)

$ \begin{aligned} W C(x)= & \min \left(\frac{1,249}{V}\right) \cdot \min \left(1,0.9 \times \frac{T I}{3}\right) \cdot \\ & \min \left(1, \frac{K}{300}\right) \cdot Y(x) \end{aligned} $

(6)

式中：Z为经验常数；AWC(x)为栅格单元x的土壤有效含水量(mm)，由土壤估算模型得到^[27]；TI为地形指数，通过区域汇水量D、土层深度S(mm)以及坡度比P计算得到^[28]；V为流速系数^[28]；K为土壤饱和导水率(cm/d)。

1.3.2 土壤保持服务

InVEST模型中，土壤保持模块在水土流失方程(USLE)的基础上考虑到地表植被对上游泥沙的拦截，弥补了USLE模型的不足，并实现精细的空间化表达^[24]。模型的主要算法为：

$ \begin{aligned} \operatorname{SEDRET}(x)= & R(x) \cdot K(x) \cdot L S(x) \cdot \\ & {[1-C(x) \cdot P(x)]+\operatorname{SEDR}(x) } \end{aligned} $

(7)

$ \begin{aligned} \operatorname{SEDR}(\mathrm{x})= & \operatorname{SE}(x) \sum_\limits{y=1}^{x-1} \operatorname{USLE}(y) \\ & \prod_\limits{z=y+1}^{x-1}[1-\operatorname{SE}(z)] \end{aligned} $

(8)

$ \begin{aligned} \operatorname{USLE}(x)= & R(x) \cdot K(x) \cdot LS(x) \cdot \\ & C(x) \cdot P(x) \end{aligned} $

(9)

式中：SEDRET(x)为栅格x的土壤保持量(t)；R(x)为栅格x的降水侵蚀力因子，通过Wischmeier公式计算得到^[29]；K(x)为栅格x的土壤可蚀性因子，通过Williams提出的公式计算得到^[29]；LS(x)为栅格x的坡度坡长因子；C(x)，P(x)分别为栅格x的植被与经营管理因子和水土保持措施因子^[28-29]；SEDR(x)为栅格x的泥沙持留量(t)；SE(x)为栅格x的泥沙持留效率；USLE(x)，USLE(y)分别为栅格x与其上坡栅格y的实际土壤侵蚀量(t)。

1.3.3 生态分区划定

基于研究区域的自然现状与发展需求，选取降雨特征以及与其显著相关的水源涵养功能与土壤保持功能作为生态分区划定的基础。k均值聚类算法是一种结合迭代概念的聚类分析方法，具有简洁性与高效率的特征^[30]。生态分区的定义运用SPSS 27中的k均值聚类工具，分析对象为2010，2015和2020年这3 a各子流域降雨、水源涵养量与土壤保持强度值。所得到的不同生态分区在降雨、生态系统服务功能以及土地利用类型的构成方面有着不同的特征，同时多年间同一子流域会发生生态分区的转移。

1.3.4 环境变量与生态系统服务的关联特征分析

进行关联特征分析时，现有研究常应用研究区域内随机点或样方的方法进行分析对象的选取，存在一定的偶然性^[31]。为保证研究对象可代表研究区域的综合特征，本文将所选分析对象进行10 km×10 km的网格化处理后进行分析。选取响应变量为水源涵养量与土壤保持强度的年间变化值；解释变量综合考虑自然因素以及人为因素，选取年间变化较为明显的降雨变化量以及各土地利用类型面积变化量，探索降雨及土地利用类型变化对生态分区变化的驱动作用。分别应用Spearman分析与冗余分析的方法，对解释变量解释单一响应变量以及响应变量集的能力进行分析。

2 结果与分析 2.1 降雨空间格局与变化特征

图 1展示了2010，2015和2020年降雨空间格局以及相邻两年之间降雨的变化特征。整体来看，3 a降雨均呈现自北向南逐渐增加的趋势。2015年与2010年降雨相比，在研究区域北侧出现了较大面积的减小，幅度达到120 mm；降雨增加区域多位于研究区域的南侧，幅度多为70 mm以上。2020与2015年降雨相比实现了全域的增加，且增加幅度较大多为320 mm以上，整体呈现东侧相比西侧增加较多的趋势。

2.2 ES空间格局与变化特征 2.2.1 水源涵养量

图 2展示了汉江中下游水源涵养量的空间分布特征，从整体来看，水源涵养量高值位于研究区域西侧以及东侧的山地地区，中部平原地区较低，整体呈现“东西高，中部低；南高北低”的分布格局，其中2010与2015年的分布格局较为接近。3个年份中，水源涵养量在90 mm以上的区域分别占研究区域面积的57.00%，56.82%，82.53%。2015与2010年相比水源涵养量存在明显的增加与降低区域，且变化特征与降雨量的增加减少特征相统一，所占面积分别为55.04%与44.96%。2020年与2015年相比水源涵养量在93.54%以上区域均实现了提升，增加较高(＞200 mm)的区域位于山地东侧边缘与中部区域，占总面积的16.15%。

2.2.2 土壤保持强度

图 3展示了汉江中下游土壤保持强度的空间分布特征，从整体来看与水源涵养量的分布格局较为类似，两者均呈现出中部较低，两侧较高的分布格局。其中2020年为土壤保持强度高值(＞600 t/hm²)分布最广的区域，占研究区域面积达到14.27%。2015年与2010年相比，土壤保持强度减少较多区域位于西侧山地，减少幅度在500 t/hm²以上；在除西侧山地以外区域变化较小，减少量在0~200 t/hm²之间；2020年与2015年相比，土壤保持强度在绝大面积区域增加，最大增加幅度在500 t/hm²以上，占研究区域总面积的10.47%。

2.3 生态分区划定及特征 2.3.1 生态分区空间分布及变化特征

通过k均值聚类的方法在子流域尺度上进行不同降雨量、水源涵养量以及土壤保持强度的生态分区划定，最终划定为5类不同的区域，特征如图 4—5所示。其中第1类区域为土壤保持强度最大的区域，第5类区域为降雨量最大的区域，第2类区域为水源涵养量最大的区域。在3个年份研究区域中比例最大的区域分别为第3类区域、第3类区域、第2类区域，面积比例分别达到了70.54%，72.92%，45.53%。相较于2010与2015年的主导区域类型第3类区域，2020年的主导区域类型第2类区域在降雨、水源涵养以及土壤保持上均有增加。

生态分区转移特征如图 6所示。当某一子流域由第a类区域转变为第b类区域时，将其定义为a→b区域，如1→3区域、1→1区域。2015年与2010年相比，生态分区保持不变的占10个子流域，面积比例为63.21%；2015年与2020年相比，保持不变的占3个子流域，所占面积比例为15.47%。

针对不同生态分区转移特征区域进行土地利用类型变化特征的分析，各区域土地利用类型保持不变的面积均大于83%；各土地利用类型中，农田与林地保持不变的面积之和比例均较高，在60%~92%之间。其中2010—2015年1→1区域中，土地利用保持不变的面积比例最大，达到96.23%；而同样作为生态分区保持不变的区域，2015—2020年5→5区域中土地利用保持不变的面积比例最小，为83.38%。

2.3.2 生态分区影响因子变化特征

经过对降雨量以及水源涵养量、土壤保持强度年间变化差值的分析，可得到变化量较大的区域均位于研究区域北侧(图 7)。3类因子中水源涵养量与降雨量在变化量与变化的正负特征上均较为一致；而土壤保持强度变化特征与降雨量变化特征一致性较差。

2.4 生态系统服务与环境变量关联特征 2.4.1 生态系统服务变化与环境变量变化的相关性分析

对研究区域整体进行的相关性分析表明，降雨量变化与水源涵养量以及土壤保持强度变化均呈强正相关关系，而在土地利用类型中呈现显著相关的为农田、林地、水体及建设用地(表 1)。其中，农田面积变化与水源涵养量及土壤保持强度变化均呈正相关，林地面积变化仅与水源涵养量变化呈弱正相关，水体面积变化仅与水源涵养量变化呈弱负相关，建设用地面积变化与水源涵养量及土壤保持强度变化均呈负相关。

2.4.2 生态系统服务变化与环境变量变化的冗余分析

对汉江流域中下游2010—2015年生态分区发生变化区域的响应变量进行去趋势分析，第一轴长度均小于3适于使用冗余分析，结果如图 8所示。

由图 8可以看出，在不同的区域中降雨量的变化(a)均对模型的贡献率较大。而在土地利用面积变化方面：在1→3区域水体面积变化(e)对模型的贡献率较大；在1→4区域草地面积变化(d)、水体面积变化(e)对模型的贡献率较大；在2→3区域，未利用地面积变化(g)对模型的贡献率较大；在3→4区域，农田面积变化(b)与未利用地面积变化(g)对模型的贡献率较大；在3→5区域，农田面积变化(b)、林地面积变化(c)、水体面积变化(e)对模型的贡献率较大。

在冗余分析的基础上对降雨量变化与土地利用变化两组解释变量进行方差分解(图 9)，明确两组变量对响应变量的解释程度，结果表明2→3区域与3→4区域两组变量可共同解释的响应变量变化程度较大，分别为68.1%与42.5%；2→3区域降雨量变化解释部分最大，为63.9%；1→4区域土地利用变化解释部分最大，为24.6%。

对2015—2020年生态分区发生变化区域的响应变量进行去趋势分析，第一轴长度均小于3适于使用冗余分析，结果如图 10所示。结果表明，在不同的区域中降雨量的变化(a)均对模型的贡献率较大。而在土地利用面积变化方面：在3→1区域水体面积变化(e)、建设用地面积变化(f)对模型的贡献率较大；在3→2区域林地面积变化(c)、建设用地面积变化(f)对模型的贡献率较大；在4→1区域农田面积变化(b)对模型的贡献率较大；在3→5，4→2，4→3区域草地面积变化(d)对模型的贡献率较大。

在冗余分析的基础上对降雨量变化与土地利用变化两组解释变量进行方差分解(图 11)，明确两组变量对响应变量的解释程度，结果表明：3→1区域以及4→1区域两组变量可共同解释的响应变量变化程度较大，分别为42.6%与69.4%；3→5区域以及4→3区域两组变量可共同解释的响应变量变化程度较小，分别为0.6%与0%；4→1区域降雨量变化解释部分最大，为69.4%；3→1区域土地利用变化单独解释部分最大，为41.6%。

3 讨论与结论 3.1 讨论

本研究在综合前人对汉江流域关键生态系统服务筛选与量化研究的基础上，以子流域为尺度进行了汉江流域中下游各生态分区空间分布及转移特征、自然及人为因素对生态分区转变的驱动能力研究。研究结果表明，水源涵养量与土壤保持强度的空间分布及变化特征与汉江流域相关研究一致^{[14, 16-17, 32]}，证明InVEST模型在该流域有良好的适应性与评估能力。另外，解释变量与响应变量的关联特征分析表明降雨量、农田面积以及建设用地面积的变化对生态系统服务变化影响较大，同时不同区域降雨量变化与土地利用变化对生态分区转变的贡献率不同，方差分解结果表明，位于研究区域西北侧的子流域更容易受到土地利用的变化导致生态分区转变。通过该方法可以明确环境变量对生态分区变化的影响程度，避免出现环境变量相关性较高但影响程度较低的情况。针对土地利用变化解释部分较大的区域，未来应加强对于关键土地利用类型的管控：严格保证耕地红线，控制建设用地的扩张趋势，并通过林地的培育减小斑块破碎的现状。针对降雨量变化解释部分较大的区域，应通过绿色及灰色基础设施建设的方式提升生态系统服务的稳定性，如通过提升水系连通性、地表水引蓄以及海绵设施的建设等增强水源涵养能力，通过生物、工程、蓄水保土措施等增强土壤保持能力；另外应加强对未来降雨量变化的预测与监测，避免因剧烈气象变化导致的生态系统服务供需失衡。目前的研究仍具有一定的局限性：从成因机制来看，南水北调中线工程自正式通水以来已累计向北输水3.00×10¹⁰ m³，造成下泄到下游的水量减小，也对生态系统服务以及生态分区的演变发生影响，而现研究无法明确其影响程度；在研究年份的选择方面，本文仅以5 a为间隔进行选择，应提高研究年份数量以减少年降雨量周期性变化的影响；在环境变量选取方面，本文仅选取了降雨变化量与各土地利用类型面积变化量进行分析，缺乏对于更多环境变量的筛选。后续将以本研究为基础增加研究年份与环境变量类型，并通过未来不同发展模式下生态系统服务预测结果优化改进生态规划与保护策略，提升现实意义。

3.2 结论

(1) 2010，2015和2020年生态分区面积比例最大的类型分别为第3，第3和第2类区域，分别占70.54%，72.92%，45.53%。在生态分区不同转变特征区域内，土地利用保持不变的区域面积比例均大于83%；而生态分区影响因子变化量较大的区域均位于北侧及西北侧。

(2) 解释变量与响应变量的相关性分析表明，降雨量变化与水源涵养量、土壤保持强度的变化均呈现显著性较高的正相关；而在土地利用类型中，农田与建设用地面积变化分别与两生态系统服务变化量呈现显著性较高的正相关与负相关。

(3) 在降雨量与土地利用变化对生态分区转变的贡献率方面，土地利用变化贡献率较高的区域面积比例较大。其中土地利用变化贡献率较高的区域均位于研究区域西侧较高高程的区域；降雨量变化贡献率高的区域位于研究区域中部及南部。