邵雅琪,王春丽,肖 玲,王美林,姜群鸥,*

1 北京林业大学水土保持学院, 北京 100083 2 北京林业大学水土保持与荒漠化防治教育部重点实验室, 北京 100083

国土空间是国家政治、经济、文化活动的场所和载体[1],然而由于区域自然和社会经济等条件的差异以及人们对国土资源利用的目的和方式的不同,国土开发秩序混乱,资源环境问题日益严重,山水林田湖草生态系统遭到破坏,生态退化明显。习近平总书记在党的十九大报告中提出了山水林田湖草系统治理[2],构建“山水林田湖草生命共同体”的基本理念,指出我国生态文明建设的方向。在这种背景下,以山水林田湖草空间格局为基本单元,开展区域山水林田湖草空间格局时空演变特征及其生态过程研究,势必将为缓解国土资源开发与生态环境保护间的矛盾提供科学参考。

目前,以山水林田湖草空间格局为基本单元的研究相对较少,相关研究多是从区域尺度分析土地利用类型转换及其功能变化规律[3-5]。而我国相关研究集中在近几年,内容主要是内涵特征、机制体制、管理路径和启示作用等方面,例如学者专家结合山水林田湖草生命共同体,对其自然资源用途管制实施路径进行探索并提出建议[6-8]；又或采用系统聚类分析方法对重庆三峡库区山水林田湖生命共同体的健康状况、格局进行诊断[9]。尽管我国已有部分以山水林田湖草生态系统为对象的研究成果,但大部分仍停留在政策层面,落实到山水林田湖草空间格局的指标定量化分析还较为少见。因此,本文将在现有研究基础上,对妫水河流域山水林田湖草的空间格局及其演变特征进行定量化分析,阐明山水林田湖草空间格局的自相关关系,揭示不同格局下的生态变化过程。这势必会为区域生态环境修复,构建“山水林田湖草生命共同体”提供重要思路。

妫水河是官厅水库的入库水系,也是延庆区重要的水域空间,2019年世园会以及2022年冬奥会的举办对妫水河流域的生态环境和水源保障提出了更高的要求。但由于社会经济发展,妫水河流域山水林田湖草各生态要素正在遭受不同程度的破坏,因此,解析流域内山水林田湖草空间格局演变特征和生态变化过程不仅对流域提升生态环境质量具有重要指导意义,也是政府进行合理开发,同时推进生态修复的主要依据之一。本文将对妫水河流域进行区域划分,结合空间自相关分析,系统地探究低山丘陵、山前台地和山区的山水林田湖草空间格局演变特征及其生态过程,并提出相应生态修复建议。

1 研究区概况

妫水河流域(图1)位于北京市西北部延庆区,属于海河流域永定河水系,是北京西北重要的生态屏障。流域总面积约926 km2,地处115°40′48″—116°16′12″E,40°13′12″—40°40′12″N。流域地势东北高,西南低,东、南、北三面环山,西面为官厅水库。妫水河流域属大陆性季风气候区,地处北温带,是半干早与半湿润的过渡地带。多年平均气温为10 ℃,多年平均降水量为390 mm,降水多集中在6—8月[10]。流域内以林地和灌溉作物为主,林地主要沿研究区边界分布,灌溉作物则主要分布在中部地区。妫水河流量主要来源于降水和地下水补给,由于近几年干旱少雨,加之研究区内城市经济发展迅速,生活用水增加,因此流量急剧下降,面临着严重的水资源短缺等生态问题。

图1 妫水河流域地理位置图Fig.1 Location of the Guishui River Basin

2数据与方法2.1研究数据

(1)土地利用数据

本研究所采用的土地利用空间分布数据是基于购买的2008、2013和2017年3期空间分辨率为2.5 m的SPOT遥感影像数据,采用人机交互解译法所得到,解译精度达96.8%。土地利用类型根据研究内容需要,分为耕地、林地、园地、草地、水域、建设用地和裸地。

(2)数字高程模型数据

数字高程模型(Digital Elevation Model，DEM)数据来源于国家基础地理信息中心,空间分辨率为30 m,基于DEM数据进行地形分析处理,提取绝对高程以及坡度等地形因子。

(3)遥感影像数据

根据遥感影像的可获取性和植被生长的季相,本研究所采用的遥感数据为2008、2013以及2017年的3期Landsat遥感影像,数据来源于中国科学院地理空间数据云(http://www.gscloud.cn/),空间分辨率为30 m。为了保证图像光谱值的原始性,对图像进行辐射定标、大气校正等预处理,并通过掩膜提取得到研究区范围影像。

2.2 研究方法

2.2.1山地区域划分

高程和倾斜度是实现山地区域划分的关键性指标[11]。本研究采用绝对高程和坡度作为描述高程和倾斜度特征的指标,结合妫水河流域的地形地貌特征以及李炳元等对中国陆地基本地貌类型的划分指标[12-13],将妫水河流域地貌类型划分为低山丘陵、山前台地及山区(图2)。具体划分指标如表1所示。

表1 妫水河流域地貌类型划分

图2 妫水河流域地貌类型划分结果Fig.2 Geomorphic type classification of the Guishui River Basin

2.2.2空间自相关

山水林田湖草空间格局演变作为地理过程现象,通过探究其空间自相关特征,可以解析某一位置的要素观测值与相邻位置上观测值的关联程度[14-15]。一般可分为全局空间自相关和局部空间自相关分析[16]。

(1)全局空间自相关

全局空间自相关可以描述地理要素属性值在整个区域的空间特征,一般采用全局莫兰指数(Global Moran′sI)进行表征。Global Moran′sI取值范围为[-1,1],指数为正表示空间正相关性,其值越大,空间相关性越明显；指数为负表示空间负相关性,其值越小,空间差异越大；指数为零表示空间分布呈随机性。计算公式如下：

(2)局部空间自相关

为了确切展示地理要素集聚或异常发生的具体空间位置,通常采用空间关联局域指标(Local Indicators of Spatial Association,LISA)进行进一步分析。LISA指数被定义为：

式中,S2是观测单元xq的方差；m为观测单元总数目；wpq是要素q和p的权重,邻接取值为1,否则为0。若Ip>0,表示该区域单元周围相似值(高值或低值)在空间上的集聚；反之,Ip<0,表示非相似值在空间上的集聚。

2.2.3植被覆盖度

随着遥感技术的发展,基于遥感可实现植被覆盖度大范围长时间序列动态变化监测,进而形成了多种植被覆盖度估算方法,目前最常用的就是基于归一化植被指数(Normalized Difference Vegetation Index, NDVI)的像元二分模型[17]。NDVI计算公式为：

式中,ρR为红光波段的反射率,ρNIR为近红外波段的反射率。

基于NDVI与植被覆盖度之间的可靠相关性,建立基于NDVI的植被覆盖度估算模型：

VFC=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

其中,NDVIsoil为完全是裸土或无植被覆盖区域的NDVI值,NDVIveg地及山区的空间关联局部指则代表完全被植被所覆盖的像元的NDVI值,即纯植被像元的NDVI值。计算公式分别为：

NDVIsoil=(VFCmax×NDVImin-VFCmin×NDVImax)/(VFCmax-VFCmin)

NDVIveg=((1-VFCmin)×NDVImax-(1-VFCmax)×NDVImin)/(VFCmax-VFCmin)

3 结果与分析

3.1 山水林田湖草空间格局时空演变特征

基于对山地区域划分的结果,采用空间动态模型,分析了2008年、2013年和2017年3期妫水河流域山水林田湖草空间格局特征(图3)。2008—2017年低山丘陵内分布着大面积的山田,但其面积存在明显的持续缩减,减少面积约为26.6 km2,一方面可能是城镇建设占用耕地所致,另一方面可能是向山林的转换。山林和山草面积变化趋势较为相似,2008—2013年为面积急速增长期,2013—2017年是平缓增长期,其面积分别增加16.4、20.3 km2。同时,山裸地的面积也有所增加,这可能主要是由于冬奥会、世园会诸多工程的修建,短时间内在一定程度上破坏了生态环境。山水面积呈现缩减态势。

2008—2017年山前台地区山田和山林的面积基数相对较大。山田主要分布在流域东北部以及南侧山地前缘,研究时段内山田面积减少10.2 km2,呈持续减少态势。由于地势原因,山前台地中的山林则多是与山区山林相间分布于流域南北两侧,2008—2017年山林面积增加24.8 km2。山草格局面积变化趋势与山田一致,减少9.72 km2,山水面积经历了先减少后增加的过程,整体减少0.4 km2。而山裸地面积则持续增加,这在一定程度上反映出山前台地区域生态环境发生退化。

山区内成片分布着大量山林,但其面积在研究时段内减少3.2 km2,原因可能在于近年来人类在山区的活动日益频繁,对山林植被造成了破坏,使其发生退化。2008—2017年山区内的山草沿着北部山区底缘向南移动,2013年以后山区山草主要出现在西庄科、张山营一带。除此以外,2013年后流域西北部出现了明显的山裸地。山田和山水年际间变化并不显著,可能在于其面积基数较小,使得格局面积年际变化程度也较小。

图3 2008—2017年妫水河流域山水林田湖草空间格局Fig.3 The spatial pattern of Mountain-Water-Forest-Cropland-Lake-Grass System in the Guishui River Basin during 2008—2017 a. 低山丘陵,b. 山前台地,c. 山区

3.2 山水林田湖草空间自相关特征

3.2.1山水林田湖草空间格局的全局自相关

本研究基于标准化后的山水林田湖草面积百分比数据,利用全局莫兰指数分析妫水河流域低山丘陵、山前台地和山区内山水林田湖草空间格局自相关的演变特征(表2)。对于整体而言,除2008年低山丘陵中的山草和山裸地全局空间自相关未通过显著性检验以外,其余山水林田湖草格局均通过P<0.01的显著性水平检验,且Moran′sI均为正值,表明不同格局的空间分布不是随机的,均在空间上表现出显著的集聚性特征,呈现出空间正自相关关系。在低山丘陵区域,山田的空间集聚程度呈持续下降趋势,2008—2017年Moran′sI减少0.07,而山林的空间集聚程度却不断上升,研究时段内Moran′sI增加了0.27。2008—2017年山水格局的空间集聚程度经历了先急速下降后趋于平稳的过程,2013年后稳定在0.22左右。对于山前台地,尽管2008—2013年山田空间集聚性不断下降,但2013年后其Moran′sI稳定在0.59,依旧表现出较强的空间集聚特征。而山林、山草的空间集聚性不断上升,2017年分别达到0.66、0.37。此外,山水和山裸地的空间集聚性均表现出先减小后增大的特征。在地势较高的山区,成片分布着大面积的山林,因此山林表现出较强的空间集聚性,但2008—2017年这种集聚程度却不断下降,原因可能在于近年来人类活动对山区整体生态环境扰动加剧,使得山林受到破坏,完整性不断降低。2008—2017年山草的空间集聚性经历了先增大后减小的过程,但整体波动幅度较小。而由于近年来山区的施工建设,山裸地空间集聚程度有较大程度的升高。除此以外,山区的山田和山水格局面积基数较小,空间集聚程度维持在相对较低的水平。

表2 山水林田湖草全局空间自相关显著性检验

Z是标准差的倍数,P表示概率,Z与P相关联,Z<-1.96或Z>1.96时P<0.05,即置信区间大于95%

3.2.2山水林田湖草空间格局的局部自相关

为了进一步揭示山水林田湖草空间格局面积百分比的高值和低值空间集聚状态,了解局部的空间差异性,基于面积百分比数据和空间权重分别获取了低山丘陵、山前台地及山区的空间关联局部指标(LISA)空间分布图(图4—图6)。

低山丘陵区域内山田的面积相对较广,呈现显著的HH、LL集聚趋势,分别主要分布在流域中西部和东北部,2008—2017年山田HH集聚区数量减少15.6%,而LL集聚区却增加了19.2%。2008年山林主要呈现LL集聚趋势,分布在官厅水库周围及流域东北部,中部有零星的HH集聚区和LH异常区分布。2013年以后HH和LL集聚区呈块状显著增加,可能是由于近年来田间造林所致。与山林相似,2008年山草空间格局分布区域随机化,2013年后官厅水库周围出现了HH集聚区域,东桑园、马坊以及榆林堡南部的山草呈现较显著的LL集聚趋势。2008—2017年山水格局HH、LL集聚区呈分散趋势,这可能是由于近年来气候干旱、用水量增多,使得部分河道断流。除此以外,流域内山裸地HL、LH异常区域较为显著。

图4 低山丘陵区域山水林田湖草空间关联局部指标分布图(P<0.05)Fig.4 LISA distribution of Mountain-Water-Forest-Cropland-Lake-Grass System in the low hilly area (P<0.05)

图5 山前台地山水林田湖草空间关联局部指标分布图(P<0.05)Fig.5 LISA distribution of Mountain-Water-Forest-Cropland-Lake-Grass System in the piedmont platform (P<0.05)

对于山前台地区域,山田呈现显著的HH、LL集聚趋势,HH集聚区主要分布在流域南部和东北部的山前平坦区,地势较为低洼,水源相对充足,山田面积占比较大且连片分布,而LL集聚区则主要出现在流域的北部。2008—2013年LL集聚区面积呈收缩态势,面积减少10.8%,2013后LL集聚区面积基本维持不变；而2008—2017年尽管流域HH集聚区位置发生一定变化,但整体数量基本保持稳定。山林与山田类似,呈HH、LL集聚趋势,但其HH集聚区主要分布在流域南北两侧,而LL集聚区则几乎由西南向东北贯穿整个流域。2008—2017年山林的HH集聚区面积增加8.2%,LL集聚区面积却减少4.5%。山草主要呈现LL集聚趋势,2008—2017年其LL集聚区位置不断向流域北侧边缘蔓延覆盖,表明2008年以后流域北侧山地出现了较多数量的低密度山草地。山水的HH集聚区与LH集聚区相间分布,说明山水格局的稳定性较低。2008年以后流域北侧山裸地面积增加,LL集聚区呈条带状分布于流域北侧。

在山区,山林呈现显著的HH、LL集聚分布,流域北部的山区山林面积较大且呈块状分布,而LL集聚区则主要分布在流域东南侧山区,2008—2017年HH、LL集聚区面积均减少约11.8%。2008年山草主要在流域西北部呈现LL集聚趋势,但2013年以后山草逐渐转为HH集聚,但HH集聚面积相较于2008年的LL集聚区面积有明显缩小,集中在西庄科、水峪及西五里营。2008年山裸地主要呈现LH异常分布,而2013年以后出现呈块状的HH集聚区,这可能在于近年来随着2019年世园会和2022年冬奥会场馆及配套基础设施建设项目的施工,不得不从山区开凿管道,导致山区地表植被破坏。除此以外,山田、山水的HH、LL集聚区域面积均相对较小,这可能是由于在山区这3种格局的面积基数较小所致。

图6 山区山水林田湖草空间关联局部指标分布图(P<0.05)Fig.6 LISA distribution of Mountain-Water-Forest-Cropland-Lake-Grass System in the mountainous area (P<0.05)

3.3 山水林田湖草生态变化过程

3.3.1妫水河流域生态特征空间格局演变分析

图7 2008—2017年妫水河流域不同地貌类型下的植被覆盖度Fig.7 Vegetation coverage of different geomorphic types in the Guishui River Basin during 2008—2017

植被覆盖度是反映生态系统和地表植被分布特征的重要参数,是衡量区域生态环境改变的重要指标。本研究估算了2008年、2013年、2017年妫水河流域低山丘陵、山前台地、山区的植被覆盖度,用以解析其空间格局的生态过程,从而为科学合理地规划生态环境保护方案奠定基础(图7)。2008—2017年低山丘陵的生态环境质量退化情况明显加剧,植被覆盖度指数降低了21.35%,可见在妫水河流域低山丘陵的耕地和建设用地聚集范围内,由于人为干扰程度的加剧,整体生态格局的的异质性下降、破碎化增加,致使生态格局的稳定性和安全性相较于山区和山前台地明显下降。对于占妫水河流域研究区面积62.74%的山前台地,2008—2017年平均植被覆盖度呈现先减小后增加的趋势,但整体上平均植被覆盖度上升了1.79%。其中高值区主要分布在研究区边缘以及西北地区,低值区主要分布在研究区中部及东北地区,可见自妫水河流域的生态建设随着经济社会和观念的发展而不断深入,生态修复情况取得了一定成效,但分布在妫水河流域山前台地的东桑园、小泥河、古城、下垙等一些乡镇的生态环境退化较为严重,应引起政府的高度重视。相较于山前台地和低山丘陵,山区的生态环境质量相对较好,其中2008年的植被覆盖度为0.86。但从2008—2013年生态环境质量也呈现了急剧退化趋势,平均植被覆盖度减少了10.93%,随着时间的推移,生态环境质量有所改善,2017年山区的平均植被覆盖度为0.82。其中高值区域主要分布在南部的里炮地区和北部杨树河以及龙庆峡等地区,低值区主要分布在陆家河和佛峪口左支的局部北部地区。因此,在未来的妫水河小流域的山区治理过程中,应充分考虑其地理特点和发展规划,有针对性的治理,从而建设生态、协调并可持续发展的小流域。

3.3.2山水林田湖草空间格局的生态演变特征

结合植被覆盖度,本研究进一步定量分析了不同山水林田湖草空间格局的生态环境质量变化特征(图8)。在低山丘陵,山林、山草整体的植被覆盖度较高,但2008—2017年均有不同程度的下降,可见山林、山草空间格局的生态稳定性较差,生态退化较为严重,极易受人类活动和经济发展的影响。其中,山林2008年植被覆盖度最高,但随着时间的变化,植被覆盖度呈现明显的下降趋势；而山草随时间呈先增加后减少的趋势,2008和2017年植被覆盖度差异较小。妫水河流域水面深度较浅,有大量生长的植物,将部分水体掩盖,因此,山水中也存在一定的植被覆盖度,但2008—2017年山水的植被覆盖度降幅明显,2017年山水的植被覆盖度仅为0.23,其生态问题较为严峻。2008—2017年山田的植被覆盖度呈现稳步上升的趋势,2017年植被覆盖度达0.40。

图8 2008—2017年妫水河流域山水林田湖草空间格局下的植被覆盖度Fig.8 Vegetation coverage of Mountain-Water-Forest-Cropland-Lake-Grass System in the Guishui River Basin during 2008—2017

就山前台地而言,山林、山草的生态退化与改善趋势大致相同,植被覆盖度呈现先降低后升高趋势,从2008—2013年,山林、山草植被覆盖度减少了11.63%、25.39%,这是由于林地、草地面积的减少直接导致了其植被覆盖度的降低；随着时间的推移,其生态质量分别有所改善,但由于此类景观类型处于人类活动干扰强烈的耕地和城镇建设用地之间,生态恢复程度较慢。山田空间格局下植被覆盖度呈现持续上升趋势,2008—2017年植被覆盖度增加了32.15%。山水的植被覆盖度呈现先增高后降低的趋势,表明自2013年后山水空间格局下生态发生退化,严重影响了流域周边生产生活活动。

相较山前台地和低山丘陵区域,山区的平均植被覆盖度较高且不同年份间不同组合类型的植被覆盖度差异较小,可见山区的生态环境较为稳定,且退化程度较弱。其中,山区格局中山林、山草、植被覆盖度差异较相似,2008年的植被覆盖度较高,分别达到了0.86、0.88,但随着时间的流逝,山林、山草格局的生态均有不同程度的退化。山田植被覆盖度从2008—2013年急速下降然后趋于平缓,且有一定的恢复。

4 结论与讨论

本研究基于土地利用数据和高程数据等,分析了妫水河流域山水林田湖草空间格局,采用植被覆盖度,解析了不同空间格局下的生态变化过程。尽管低山丘陵和山前台地区域中的山田在2008—2017年间由于田间林地建设空间集聚程度下降,但其植被覆盖度呈现明显上升趋势。因此,建议在妫水河流域内的不同田块之间间作林地,形成植被过滤带。在林地植被的拦截作用下,可有效阻拦农田面源污染物随地表径流的运移,起到一定的污染物净化效果；同时,相对高大的林地植被也可以起到天然防风屏障的作用,减少沙尘暴。

基于本文对山林格局及植被覆盖度年际变化的分析可知,妫水河流域山区山林、山草复合格局相对较少,因此建议在对山区原始植被进行封禁的同时,加强对人为管护的山林地区草地植被保护,形成“山林-山草”复合格局。在林下种植一些具有经济效益的草本植物,涵养水源,修复土壤,以期在提升景观的同时也为当地居民带来一定的经济收益。另外,在低山丘陵,河流沿岸有较大数量的“山草”格局分布,因此我们建议可在该区域内种植水生植物,水位线以上种植观赏乔灌木,树下间作花卉和药材,同样可形成“山水-山林-山草”复合格局,打造妫水河沿线生态河滨带[18]。

妫水河流域近年来由于地区经济开发建设,用水量增加[19],导致河流断流现象严重,在空间上表现为集聚程度逐渐降低,同时,延庆区旅游业的不断发展也加重了城市绿化保育用水的负荷,低山丘陵、山前台地以及山区的山水格局下的植被覆盖度均维持在相对较低的状态,这在一定程度上表明山水格局的生态问题日益加重[20],因此,山水格局的生态修复工作仍是山水林田湖草综合治理的关键所在。本文针对当前妫水河流域低山丘陵、山前台地及山地的山水格局及其生态现状提出以下生态修复建议：以防治并重,治管结合的原则对低山丘陵和山前台地区域进行治理,开展农林、水利、水土保持相结合的流域综合治理工程,加强河道治理措施等,增加生态基流,恢复断流河道。针对山区范围内的山水区域妥善处理地面径流,防治风力侵蚀、水力侵蚀,减少人为干扰,保护山区水脉。