Dynamic changes and development of water retention in Liandu District

莲都区地处浙江省丽水市瓯江上中游，位于28°06′~28°44′N，119°32′~120°08′E，总面积为1 502 km²。地处括苍山、洞宫山、仙霞岭3条山脉之间。境内地形可分为河谷平原、丘陵、山地3种，其中，低丘和高丘占全区总面积的57.0%，低山、中山面积占全区总面积的30.2%，平均海拔为381 m。气候属中亚热带季风气候类型，雨水充沛，温暖湿润，四季分明，具有明显的山地立体气候特征^[18]。

利用InVEST模型对水源涵养功能进行评价。该模型主要根据水量平衡原理，通过降水、植物蒸腾、地表蒸发、根系深度和土壤深度等参数计算产水量。在产水量的基础上，考虑土壤厚度、渗透性、地形等因素的影响，利用地形指数、流速系数和土壤饱和导水率对产水量进行修正。根据式（1）计算获得水源涵养量：

式（1）中：R为水源涵养量（mm）；V_e为流速系数；K为土壤饱和导水率（cm·d^-1），根据式（2）计算；I_T为地形指数，无量纲，根据式（3）计算；Y为产水量，根据式（4）计算。

土壤饱和导水率基于Cosby土壤饱和导水率传递函数计算。式（2）中，S和C分别代表砂粒和黏粒的含量（%）。

式（3）中：D_A为集水区栅格数量；S_D为土壤深度；P_S为百分比坡度。

式（4）中：Y_jx为年产水量；P_x为年均降水量；$A_{xj}^{\rm{ET}} $为土地利用类型j栅格上单元x的年均蒸散量，由式（5）计算。

式（5）中：R_xj为土地利用类型j上栅格单元x的干燥指数，无量纲，为潜在蒸发量与降水量的比值，根据式（6）计算；${{\omega }_{x}} $为植物年需水量和降水量的比值，是描述自然气候和土壤性质的参数，根据式（8）计算。

式（6）中：k为作物系数，是作物蒸散量E^T与潜在蒸散量$E_{\rm{0}}^{\rm{T}} $的比值；$ E_{\rm{0}}^{\rm{T}}$是指假设平坦地面被特定矮秆绿色植物全部遮蔽，同时土壤保持充分湿润情况下的蒸散量，采用式（7）计算。

式（7）中：$ E_{\rm{0}}^{\rm{T}}$为潜在蒸发量（mm·d^-1）；R_A为太阳顶层辐射（MJ·m^-2·d^-1）；T_avg为日最高气温均值和日最低气温均值的平均值（℃）；T_D为日最高气温均值和日最低气温均值的差值（℃）；P为月均值降水量。

式（8）中：Z为Zhang系数，表征降水季节性特征的一个常数，其值为1~10，降水主要集中在冬季时，其值接近于10，降水主要集中与夏季或季节分布均匀时，其值接近于1。$A_{x}^{\rm{WC}} $为植被有效可利用水。

元胞自动机（cellular automata, CA）是由Von NEUMANN等在20世纪40年代最早提出的一种时间、空间和状态都离散的网格动力学模型。马尔可夫链（Markov）是预测事件发生概率的一种数学方法，是基于随机理论形成的一种具备“无后效性”的特殊随机运动过程，即t+1时刻的状态只与t时刻的状态有关^[19]。CA-Markov模型是综合了Markov链对长时间序列的预测和CA元胞自动机基于空间关系和规则动力学模拟的优点，能更加准确地从时间和空间上模拟土地利用类型变化的情况。

式（9）中：S_t，S_t+1为t，t+1时刻土地利用系统状态；P_ij为状态转移矩阵。式（10）中，S′_t, S′_t+1为t，t+1时刻元胞状态；U_d为d维的元胞空间及网格单元；N_t为t时刻邻居的状态组合；f为转换规则。

根据模型所需参数，对各参数进行获取及本地化处理。以1994年5月12日、2001年4月13日Landsat 5遥感影像、2008年7月5日Landsat 5遥感影像以及2015年10月13日Landsat 8 OLI影像为数据源，参考“GLC2000土地覆被分类系统”将土地覆被类型^[20]主要分为裸地、农田、城镇、水体、针叶林、灌丛、阔叶林、竹林和针阔混交林等9类，采用随机森林法进行土地覆被类型分类。降水等气候数据根据莲都区及周边地区气象站记录的各年份数据进行插值获得，2030年气候数据来自于国家气候中心参与政府间气候变化专门委员会（IPCC）时根据北京气候中心模型模拟所得^[21-22]；土壤深度参考浙江省森林资源二类清查数据和1:50万浙江省土壤数据；根系深度及流速系数参数根据研究区前期外业调查和周边地区已有研究，进行本地化处理获得，尽量保证参数符合实际情况；土壤饱和水率^[23]参考1:50万浙江省土壤数据，利用COSBY等^[24]的土壤传递函数计算获得；植物可利用水根据外业调查所得田间持水量数据，参考周文佐^[25]的研究结果进行处理获得；地形数据来源于日本地球遥感数据分析中心提供的30 m分辨率全球数字高程模型数据。

本研究利用InVEST模型计算水源涵养功能时均以30 m × 30 m分辨率的栅格图为参数来源，同样也是以该分辨率栅格作为最小单位对结果进行输出计算。输出结果为水源涵养深度（mm）。通过单位换算汇总得出研究区水源涵养总量（m³）。以乡镇矢量边界对研究区各地区水源涵养量进行统计，从而分析研究区水源涵养功能的数量及空间变化特征。

对遥感影像进行解译后，由于地物光谱特征存在“异物同谱”和“同物异谱”现象，因此，2015年土地覆被类型是利用高精度影像和野外实地调查数据对研究区土地覆被类型分类进行精度验证。对于1994，2001和2008年土地覆被类型分类无法进行野外实地调查验证的问题，根据土地覆被类型具有一定规律性的特点，结合相关年份的统计年鉴、森林资源二类小班调查数据和当时的土地利用现状图等辅助方法加以验证。验证结果显示，研究区1994，2001，2008和2015年土地覆被类型遥感解译精度分别为76.3%，78.6%，81.4%和84.3%，达到了研究所需要求。

从表 1可以看出：2015年研究区针叶林的面积最大，达6.08×10⁴ hm²，其次是针阔混交林和阔叶林，分别为2.55×10⁴ hm²和2.16×10⁴ hm²，裸地、水体和灌木林的面积较小。1994-2015年，城镇、阔叶林、竹林和针阔混交林面积均有不同程度增加，其中，针阔混交林增幅最大，其次分别是城镇、阔叶林和竹林。裸地、农田、水体和针叶林面积均有不同程度的下降，其中，针叶林和裸地下降幅度最大，其次分别是农田和水体面积。

从图 1研究区土地覆被空间分布状况可以看出：1994-2001年，研究区中北部地区裸地面积显著减少，被以针叶林为主的林地所替代；竹林、阔叶林和针阔混交林面积增长明显，主要集中在研究区南部，峰源乡和莲都林场等地。2001-2008年，裸地面积持续减少，竹林、阔叶林和针阔混交林面积持续增加，在北部和东部部分地区增长明显；同时在此期间，城镇面积增长迅速，主要在万象街道和白云街道原有城镇基础上呈现辐射型扩张，并在南明山街道地区发展形成新城区。2008-2015年，城镇面积持续在研究区中部南明山街道等地快速扩张。同时，研究区北部和南部地区针叶林面积下降，阔叶林和针阔混交林面积提高。1994-2015年，莲都区北部地区裸地被以针叶林为主的林地替代，南部地区针叶林逐渐被竹林、阔叶林和针阔混交林所替代，土地覆被类型空间分布格局变化较大。此外，21 a间莲都区中部城镇扩张明显，人为活动增加，对生态环境干扰加剧。

Figure 1.  Spatial characterisitic of land cover from 1994 to 2015

通过查询黄渡、小白岩、上显滩3个水文站和五里亭水库的水文记录，获得1994，2001，2008和2015年研究区的径流深度数据。根据气象数据查询可知：研究区1994，2001，2008和2015年降水量分别为1 453.0，1 460.0，1 147.0和1 674.0 mm。根据相关研究^[26]可知：蒸散量约占实际降水量的30%。通过研究区降水量、蒸散量和径流深度可以计算水源涵养深度（V），公式为V=P-E-C。其中：P为降水量，E为蒸散量，C为径流深度。将相关参数代入计算，得出研究区各年份各水文记录点的水源涵养深度值，并将实测结果与InVEST模型所得估算结果进行检验，结果如图 2所示。由图 2可知：估算值对实测值的回归方程为y=0.9x+36.17，R²为0.83，表明方程拟合程度较高，InVEST模型能够较为真实地反映研究区水源涵养功能。

Figure 2.  Accuracy test of InVEST model

1994-2015年，研究区水源涵养功能提高明显。从表 2可以看出：1994年研究区水源涵养总量为1.11×10⁸ m³，水源涵养深度为69.98 mm；2001年研究区水源涵养总量为1.57×10⁸ m³，水源涵养深度为98.40 mm；2008年研究区水源涵养总量达1.23×10⁸ m³，水源涵养深度为82.36 mm；2015年研究区水源涵养总量达1.83×10⁸ m³，水源涵养深度为114.01 mm。21 a间研究区水源涵养功能提高了62.92%。研究区水源涵养功能提高并非匀速，不同时期变化情况有所不同。1994-2001年水源涵养总量提高了0.46×10⁸ m³，水源涵养深度提高了28.42 mm，年均水源涵养功能增长5.10%，2001-2008年水源涵养总量减少了0.34×10⁸ m³，水源涵养深度降低了16.04 mm，年均水源涵养功能减少2.33%; 2008-2015年，水源涵养总量增长了0.60×10⁸ m³，水源涵养深度提高了31.65 mm，年平均增长率为5.49%。但从研究区各年度水源涵养深度占降水量比例来看，1994年为4.82%，2001年为6.47%，2008年为7.18%，一直表现增长趋势，但2015年由于降水量过高，水源涵养深度占降水量6.81%，稍有下降。从表 2可以看出：21 a间研究区各乡镇水源涵养功能提高程度差异明显，其中，水源涵养总量增长明显的乡镇为雅溪镇和黄村乡，1994-2015年分别增长了9.72×10⁶ m³和9.68×10⁶ m³，分别提高了75.44%和71.42%。水源涵养总量增长较少的是万象街道，增长0.02×10⁴ m³，仅提高了5.71%。水源涵养深度提高最为明显的依然是雅溪镇和黄村乡，分别提高了82.70 mm和92.73 mm，提高比例分别为75.12%和70.94%。万象街道的水源涵养深度提高最少，仅1.10 mm，提升比例仅为4.07%。从表 3可以看出：研究区水源涵养功能变化呈现先增后减再增，总体提高的趋势。但各乡镇地区水源涵养功能前中后期变化特征并不一致，其中万象街道、莲都林场和大港头镇后期年平均增长速率高于前期；白云街道、紫金街道、仙渡乡、丽新畲族乡、峰源乡和南明山街道前后期年平均增长速率相近；其余各乡镇前期年平均增长速率高于后期。

从图 3可以看出：研究区1994，2001，2008，2015年的水源涵养功能空间差异明显。1994年水源涵养深度为69.98 mm，仅东北部黄村乡、北部雅溪镇和南部峰源乡部分地区水源涵养功能较高。2001年水源涵养深度为98.40 mm，南部峰源乡、莲都林场和东北部黄村乡水源涵养功能好于其他地区，个别地区水源涵养深度达到300 mm以上，中部地区水源涵养功能较低；水源涵养深度较大的区域面积明显增大。2008年研究区水源涵养深度为82.36 mm，水源涵养功能空间格局整体与1994年相似，但略好于1994年。2015年水源涵养深度为114.01 mm，水源涵养功能较高的地区为南部峰源乡和莲都林场，其次是东北部黄村乡和雅溪镇，水源涵养功能较低的区域主要分布在中部白云街道、万象街道、南明山街道和碧湖镇等地。水源涵养深度较大的区域面积明显进一步增大。1994-2001年，研究区水源涵养深度仅中部城镇和农田呈现减少趋势，其他地区均有不同程度的增长，大部分地区增长量为20~40 mm，部分地区增长量达40 mm以上。2001-2008年，研究区水源涵养功能下降明显，其中北部和中部地区降幅较大。2008-2015年，研究区水源涵养功能增长迅速，其中中部和北部地区水源涵养功能变化最为剧烈，南部地区水源涵养功能已处于较高水准，因此增长并不明显。1994-2015年，研究区各地区水源涵养功能均有不同程度的增长，其中南部峰源乡和莲都林场尤为明显，其次是莲都区北部黄村乡、雅溪镇、岩泉街道和太平乡等地，中部地区虽有一定程度的提高，但部分地区水源涵养功能仍处于较低的水平。对比土地覆被类型变化和水源涵养功能可以看出，森林植被面积的提高能有效改善区域水源涵养功能状况。不过研究区2008年由于整体降水量较低，导致当年水源涵养量并不理想。

Figure 3.  Spatial distribution characteristics of water retention in Liandu District

选取2015年研究区水源涵养结果及其相关的12项因子（降水量、潜在蒸发量、植物可利用水、土壤深度、流速系数、土壤饱和导水率、地形指数、数字高程数据、实际蒸散量、森林面积比例、裸地城镇面积比例和阔叶、针阔混交林面积比例），基于Arc GIS平台，使用Arc GIS流域分析算法，将研究区划分为122个小流域，并以小流域为单位统计水源涵养量及各因子平均值。利用SPSS 19.0对统计数据进行逐步回归分析，保留显著因子、剔除非显著因子，方程最后得到R²为0.894，表明方程拟合程度较好，可信程度较高。由表 4可知：土壤饱和导水率对研究区水源涵养功能的正面影响最大，阔叶、针阔叶混交林面积比例和降水量的正面影响次之，并且土壤饱和导水率会随着植被的改善和群落的演替而提高^[27-28]，因此，森林结构的改善将进一步改善水源涵养；裸地、城镇面积比例有较高的负面影响，实际蒸散量的负面影响程度相对较小，但具有较高的负相关性。由此可以看出：影响研究区水源涵养功能的因素主要可分为两大类：气候因素和下垫面特征，其中下垫面特征对水源涵养功能的影响程度更高。下垫面特征因素表明当区域内土壤属性改善，阔叶林、针阔混交林面积增加，裸地及城镇面积减少时，相应的会增强区域内的水源涵养功能。气候因素表明当区域内获得更多的降水，并同时减少蒸散的情况下，水源涵养功能将随着提高。以上因素是相辅相成，共同对水源涵养功能发挥作用的。研究区1994-2015年水源涵养功能提高的原因，主要是21 a间研究区森林面积的增加和结构的改善。

基于1994和2001年土地覆被类型分类结果数据，利用CA-Markov模型对莲都区2015年的土地覆被类型进行模拟，利用模拟结果与2015年土地覆被实际类型进行对比，检验其精度，得到Kappa系数为0.72，为中等一致性，模拟结果可信。基于2001和2015年土地覆被类型分类结果，利用CA-Markov模型对2030年莲都区土地覆被类型进行模拟预测（表 5）。如图 4所示：研究区2030年土地覆被类型整体与2015年差异不大。城镇、阔叶林、竹林和针阔混交林均有不同程度的增长，其中针阔混交林增幅最大，面积和比例分别增长了3.40×10⁴ hm²和2.27%，城镇、阔叶林和竹林的面积比例分别增长了1.10%，1.81%和1.49%。针叶林、灌木林、裸地、农田和水体呈现下降趋势，其中针叶林降幅最大，面积和比例分别下降了9.57×10⁴ hm²和6.37%；裸地、农田、水体和灌木林略有下降，降幅较小。从空间上来看，2015-2030年，研究区北部变化相对明显，原有针叶林逐渐被阔叶林、针阔混交林和竹林替代，变化部分主要发生于斑块交接区域。研究区城镇在原有基础上稳定发展，在瓯江两岸逐渐扩张并相互连接，同时碧湖镇的村落也逐步发展增大。虽然研究区城镇发展分布广泛，但由于研究区城镇周边公益林面积较大，限制了城镇的向外快速扩张。

Figure 4.  Prediction of land cover and water retention of Liandu District in 2030

预测结果显示，研究区2030年水源涵养总量为1.81×10⁸ m³，水源涵养深度为112.90 mm，与2015年相比整体差异较小，略有下降。2030年研究区水源涵养能力最好为峰源乡和黄村乡，其水源涵养深度分别为228.27和223.29 mm，其次为莲都林场、市白云山林场、岩泉街道、雅溪镇、太平乡、紫金街道和老竹畲族乡，水源涵养深度均在100 mm以上，其余乡镇地区水源涵养深度均在100 mm以下。黄村乡和雅溪镇水源涵养总量较高，分别为2.32×10⁷ m³和2.23×107 m³。根据模拟结果来看，莲都区在未来增强水源涵养功能的任务仍十分艰巨。重点应该关注北部和西部地区，这2个地区有较高的水源涵养潜力。建议对其林型结构进行逐步改善，提高当地的生物多样性水平，增强生态系统稳定性。南部地区为水源涵养最高的地区，应该继续维护其目前的森林建设水平，保持现状。中部地区是以城镇分布为主，难以实现较好的水源涵养能力。因此，应该注重城镇周边山体植被修复和城市内部绿地建设，有效利用各类绿地，提高单位面积的水源涵养能力。

1994-2015年，莲都区土地覆被类型变化明显，森林面积和结构改善显著。针阔混交林、阔叶林和竹林面积分别增加了1.42×10⁴，1.06×10⁴和0.84×10⁴ hm²，裸地面积减少了1.91×10⁴ hm²。森林面积的增加有效地改善生态系统环境，有利于莲都区水源涵养功能建设和区域可持续发展。

21 a间，莲都区水源涵养功能提高了69.92%，其中前期的增速要高于后期；从空间分布上来看，莲都区水源涵养功能南部地区好于北部地区，中部地区最差，其中雅溪镇和黄村乡的水源涵养功能最好，21 a间水源涵养深度分别提高了82.70和92.73 mm，万象街道的水源涵养功能最差，21 a间仅提高了1.10 mm。

2030年模拟结果显示，莲都区水源涵养功能与2015年差异不大，水源涵养深度和总量分别为112.90 mm和1.81×10⁸ hm³，西部和北部地区水源涵养提升潜力较大。

气候因素和下垫面特征是影响莲都区水源涵养功能的2个主要因素。下垫面因素中，土壤饱和导水率和不同土地覆被类型的面积比例对水源涵养功能有重要影响；在气候因素变化难以控制情况下，唯一措施就是通过改善植被来提高区域水源涵养功能，保障区域水资源安全。