徂徕山—莲花山地区环境生态指数时空变化与驱动因素分析

马 超，蔡盼丽

(1.河南理工大学测绘与国土信息工程学院，焦作 454000；2.矿山空间信息国家测绘与地理信息局重点实验室，焦作 454000；3.中国科学院大学资源与环境学院，北京 100190)

0 引言

20世纪90年代，国际地圈—生物圈计划(international geosphere-biosphere program，IGBP)和全球环境变化中的人文领域计划(international human dimensions programme，IHDP)联合提出“土地利用和土地覆被变化(land use and cover change，LUCC)”这一核心研究领域[1-2]，使土地利用变化研究成为目前全球变化研究的前沿和热点。区域变化过程(process)、格局(pattern)与驱动力(driver)的对比研究被认为是揭示全球与区域尺度土地利用/覆盖时空变化规律的有效方法[3-5]。21世纪以来，随着3S技术的迅速发展，我国LUCC时空变化研究增温增速，如刘纪远等[6]依据5 a土地利用变化区域分异的显著特征，从自然地理、经济发展与国家宏观政策因素综合考虑，揭示了21世纪初土地利用变化的空间格局与驱动因素；佟光臣等[7]采用1986年、1996年、2002年和2013年的4期土地利用/覆被数据从土地利用类型结构、变化速度、转化关系和景观格局等方面分析了南京市近30 a间土地利用/覆被和景观格局的时空变化及其驱动因素。但LUCC驱动力的研究仍显不足，主要体现在：①研究目标较为单一，主要研究集中在经济发展较快、快速城镇化区域[8]，如京津冀地区[9]、三江平原[10]、长江三角洲地区[11-12]、黄河三角洲区域[13]和新疆乌鲁木齐[14-15]等地区；②研究方法有局限性，对于土地利用变化与驱动因素的研究方法集中在主成分分析法、转化矩阵等，受主观因素的影响大，缺乏客观性对比；③研究内容较为有限，通常从面积变化反映土地利用变化，考虑经济、政策、市场等人为因素及自然因素对土地利用变化的影响不充分。

山东省徂徕山—莲花山地区外受泰岳山脉环绕，内为两座中型工业化城市及所属人口密度极高的乡村所包围，受到自然隔绝与人类活动的双重影响，人地关系紧张，生态格局演化剧烈，生态结构退化严重。该区域的生态自然演变与人类改造行为交织，自然生境与人类活动时空交错，人地关系促成了内陆自然保护区典型的岛屿型生境，而目前对于该区域的研究仅限于文化考古[16-17]、动植物分布[18-20]、地质构造[21-22]和旅游开发[23-24]等。土地利用是生态环境变化的重要组成部分，该区域土地利用与驱动因素分析研究，将对自然保护、旅游开发和民生改善具有现实意义。

本文利用(1982—2016年)35 a间的遥感数据和(1982—2012年)31 a的气象数据，通过一系列预处理，获得该区域的裸地、水体、建筑用地、耕地、林地信息，时间序列归一化植被指数(normalized difference vegetation index，NDVI)和归一化建筑指数(normalized difference barren index，NDBI)，分析区域内建筑用地、耕地、林地的几何变化与NDVI和NDBI的时空变化特征，区域内的气候变化以及“气候—土地利用—生态指数”之间的相关性，辨识该区域引起景观格局与土地利用变化的驱动因素，为内陆封闭生境的生态环境变化及可持续发展提供参考。

1 研究区概况与数据源

1.1 研究区概况

山东省徂徕山—莲花山地区，属暖温带季风气候区，年均气温11～14 ℃，年均降水量在550～950 mm之间。徂徕山(图1)位于山东省泰安市东南20 km，E117°13′50″，N36°02′24.8″；核心区占地253.368 km2，主峰太平顶海拔1 028 m。莲花山位于山东泰沂山脉中段，E117°39′49″，N36°01′44″；南北兼跨新泰、莱芜2市，核心区占地354.398 km2，主峰海拔999 m。两山相对高度虽然不大，但是坡度较陡，坡度大部分在25°～35°之间，个别地方在45°以上，核心区不适于耕种。

图1 徂徕山—莲花山地区所处的地理位置Fig.1 The geographical location of the Culai-Lianhua area

1.2 数据源

研究采用Landsat 3/5 MSS/TM及Landsat 8 OLI(1982—2016年)数据。根据研究区8月份的植被指数达到最大值的特点，以7月27日—10月03日间、云量<10%作为下载阈值，每隔2～3 a下载一期，共获得16期遥感影像，空间分辨率30 m(表1)。各期数据均进行了土地利用分类。

表1 遥感影像列表(1982—2016年)Tab.1 List of remote sensing images (1982—2016)

气象数据来自国家气象中心网站。选取济南、泰山、淄川、沂源和泰安5个站点的31 a(1982—2012 年)数据进行了克里金插值。数字高程模型(digital elevation model，DEM)数据采用空间分辨率为30 m的ASTER GDEM。

2 研究方法

2.1 流域分析

徂徕山—莲花山地区(图2)由牟汶河、柴汶河(汇成大汶河)和苗庄河围限，利用ASTER GDEM进行流域分析，获得该区域的分水岭、汇水线和等高线，上述3条河流的汇水线勾勒出一个完整的研究区边界。

图2 徂徕山—莲花山区域范围图Fig.2 The total area of the Culai-Lianhua basin

2.2 NDVI计算

经辐射定标、大气校正和波段运算后获得研究区1988—2016年的NDVI值(1982年NDVI数据异常，未采用)。对其进行密度分级(图3)，获得该区域时间序列上的植被覆盖空间变化趋势图。

(a)1982-09-15 (b)1988-08-04 (c)1993-09-09 (d)1998-09-01

(e)2003-09-15 (f)2009-08-30 (g)2014-07-27 (h)2016-09-02

图3 研究区时间序列NDVI密度分级图

Fig.3TimeseriesNDVIdensitysegmentationmapping

2.3 人类活动环境要素获取

在对原始影像进行预处理后，采用IGBP GLCC V2.0全球土地覆盖分类体系对徂徕山—莲花山地区进行监督分类，将研究区分为裸地、湖泊、建筑用地、耕地、林地和未分类6种土地利用类型(图4)。

(a)1982-09-15 (b)1985-07-27 (c)1987-10-05 (d)1988-08-04

(e)1991-09-30 (f)1993-09-09 (g)1998-09-01 (h)2001-09-09

(i)2003-09-09 (j)2004-10-03 (k)2006-08-22 (l)2009-08-30

(m)2011-09-21 (n)2013-09-16 (o)2014-07-27 (p)2016-09-02

图4 1982—2016年研究区土地利用状况图

Fig.4TimeseriesofLUCCthematicmappingduring1982—2016

可以看出建筑用地面积呈上升趋势，耕地面积呈先上升后下降趋势，林地呈先减少后增长的趋势。据此可以对气候要素(降水、气温)和人类活动要素(NDBI、建筑用地、耕地和林地面积)以及生态指标(NDVI值)进行相关性分析，探讨研究区域生态环境变化的相关因素及主要驱动力。

2.4 气候变化分析

在ArcGIS的支持下，对研究区40 km范围内的5个站点31 a的年平均降水量和气温进行克里金插值[25]，获得徂徕山—莲花山地区的平均气温和平均降水变化趋势图(图5)。可以看到：

(a)气温变化趋势 (b)降水变化趋势

图5 1982—2012年徂徕山—莲花山气温和降水变化趋势

Fig.5Thetrendoftemperatureandprecipitationinthestudyareaduring1982—2012

1)研究区年均气温如图5(a)所示，介于10.6～12.7 ℃之间，总体呈上升趋势。31 a间出现2个高值，即12.7 ℃(1998年)和13.2 ℃(2007年)；出现3个相对低值，即10.9℃(1984年)和10.6℃(1993和2011年)。对其进行滑动平均，气温变化可以分为3个阶段：第一阶段(1982—1993年间)气温变化起伏较小，平均温度为11.52 ℃(图5(a)中的t1值)；第二阶段(1994—2003年间)气温波动增强，平均气温为11.98 ℃(图5(a)中的t2值)；第三阶段(2004—2012年)气温变化最为剧烈，最高气温13.2 ℃(2007年)和最低气温10.6℃(2011年)，平均温度为11.74 ℃(图5(a)中的t3值)。总体来说，气温变化剧烈，整体呈上升趋势，与全球变化的事实吻合[26]。

2)研究区年均降水趋势如图5(b)所示，在波动中上升，以10 a为一周期进行比较，有3期相对低值，即1989年的424.5 mm、1992年的606.752 mm和2002年的504.1 mm；有4期相对高值，即1990年(1 138.5 mm)、1994年(1 022.5 mm)、2003年(1 094.1 mm)和2004年(1 153.4 mm)。对其进行滑动平均，降水变化也分为3个阶段：1982—1991年间平均值为744.3 mm(图5(b)中的p1值)；1992—2001年间平均值为788.0 mm(图5(b)中的p2值)，2002—2012年间平均值为833.4 mm(图5(b)中的p3值)。

当然，气候的变化由多种因素所引起，并非简单的线性变化。在2007年第四次和2013年第五次IPCC评估报告[27]中指出：人类活动是近50年影响气候变化的主要因素。在全球气温变化总体上升的趋势下，区域的气候变化有不同的变化。1982—2012的31 a间，研究区气温在起伏波动中呈上升趋势，增速为+0.046 ℃/10 a；降水量呈缓慢上升趋势，速度为+43.0 mm/10 a。

3 结果与分析

3.1 人文因素变化分析

通过对监督分类结果进行统计，获取了1982—2016年间以5 a为间隔的耕地、林地和建筑用地面积变化信息，如图6所示。

(a)耕地面积变化 (b)林地面积变化 (c)建筑用地面积变化

图6 1982—2016年间徂徕山—莲花山土地利用变化

Fig.6TrendoftheLUCCinCulai-Lianhuabasinduring1982—2016

由图6可知：①耕地面积(图6(a))先增后降。除2016年外，整体呈波动性变化，无明显的突变现象。造成这种状况的原因较多，影像时相差异会产生误差，另外与国家提倡保护环境、退耕还林等土地政策不无关系。②林地面积(图6(b))呈先降后增的变化趋势，在2004年达到最低水平，2005年之后林地面积呈快速上升趋势，特别是2016年林地面积有较大跃升。主要原因之一为针对环境资源环境压力逐渐加大，国务院提出落实科学发展观加强环境保护的政策，决心扭转“重经济，轻环境”的倾向，坚持“环境优先，保护优先”。林地面积的变化与耕地面积变化的趋势相反，呈显著的负相关关系。③建筑用地面积(图6(c))上升趋势显著。除1998年短暂下降外，建筑用地面积均持续上升，2003年面积是1982年的2倍，社会经济发展，城镇化建设的加快是建筑用地面积的增加的一个主要因素。这一发现与采用NDBI的研究结果具有一致性。整体来看，建筑用地面积增长显著，耕地与林地面积此消彼长，可以发现人类活动导致的土地利用变化显著。

为验证上述结论的可靠性，与国家基础地理信息中心的最新研究成果GlobeLand30(glc30.tianditu.com)进行了对比，如图7所示。该区域2010年相对于2000年耕地面积减少了1.6%，林地面积增加了1.9%，居民地面积增加了5.8%；虽然成像的季节不同，仍与监督分类获得的2001年、2011年结果趋势一致。相比于2001年，2011年耕地减少了9.3%，林地增加了15.7%，居民地及建筑用地增加了20.5%。

(a)Global Land 30_2000 (b)Global Land 30_2010

图7 2000，2010年徂徕山-莲花山地区GlobeLand30产品

Fig.7GlobeLand30productsinCulaiLianhuaareain2000and2010

NDBI可准确反映建筑用地信息，其数值越大建筑用地比例越高，建筑密度越高。对所获得的数据进行一系列处理得到了区域NDBI均值的变化信息。为了更准确地进行分析比较，本文分别对核心区、核心区周围地段进行了统计分析，结果如图8所示。

(a)核心区NDBI均值变化 (b)核心区周围NDBI均值变化

图8 1988—2016年徂徕山—莲花山NDBI均值变化

Fig.8VariationtrendofNDBIinCulai-Lianhuabasinduring1988—2016

在核心区内(图8(a))，海拔较高、坡度较大，NDBI均值波动较小，整体呈下降趋势。有3期相对高值，为1998年(-0.183)，2003年(-0.151)，2016年(-0.192)，1期相对低值2014年(-0.254)，在29 a(1988—2016年)内核心区的NDBI均值降速为-0.014/10 a。在核心区周围地段(图8(b))，地势较为平坦，NDBI均值变化幅度较大，整体呈上升趋势。峰值出现在2003年(-0.144)，2016年(-0.181)，谷值出现在1993年(-0.203)，2014年(-0.254)，NDBI均值增速为0.003 5/10 a。

3.2 植被NDVI年均值变化

植被指数是反映生态环境变化的一个重要指标，为对研究区进行更为准确的分析比较，也将研究区分为核心区、核心区周围地段2个分区，分别提取2个区域的NDVI均值变化，如图9所示。

(a)核心区NDVI均值变化 (b)核心区周围NDVI均值变化

图9 1988—2016年徂徕山—莲花山NDVI均值变化趋势

Fig.9VariationtrendofNDVIinCulai-Lianhuabasinduring1988—2016

核心区内(图9(a))地势较高，坡度较大，不适于耕种，NDVI均值(0.61～0.69)处于较高水平。在1993年(0.661)、2009年(0.684)出现极大值，核心区NDVI整体呈上升趋势，增速为0.006 9/10 a；核心区周围地段(图9(b))人类活动频繁，1988—2016年NDVI呈下降趋势，2003年出现的低值(0.547)对应于NDBI的高值，可能与人类活动有关，该区域NDVI均值呈下降趋势，下降速度为0.018/10 a。

3.3 NDVI生态剖面分析

在NDVI影像中，沿南北两个方向分别对徂徕山、莲花山做剖面，得到徂徕山—莲花山区域时间序列NDVI剖面图，如图10—11所示。穿过徂徕山区域的NDVI南北剖面(图10)变化较为剧烈。结合遥感影像，剖面两端为耕地，且中间分布着较多居民地，NDVI值变化较大；中部NDVI值较高区域(224～542像元)为林地，整体呈耕地—林地—耕地的结构变化，林地—耕地之间为生态交错区。对莲花山也做南北剖面，获得了NDVI变化图(图11)。该图中，NDVI剖面两侧为耕地；剖面中部区域(244～628像元，684～830像元)为林地；在187～243像元和831～879像元处NDVI值有较大变化，因该区域存在裸地和建筑用地，原有的水域遭到了分割破坏；629～683像元处，植被恢复成效显著。

图10 1988—2016年间徂徕山NDVI影像南北剖面图Fig.10 North-south section of NDVI in Culai basin during 1988—2016

图11 1988—2016年莲花山NDVI影像南北剖面图Fig.11 North-south section of NDVI in Lianhua basin during 1988—2016

对比时间序列徂徕山—莲花山两条NDVI南北剖面，林区部分内部结构稳定性较好，耕地及建筑用地部分稳定性差，交错区变化较为剧烈。土地利用类型较少，植被覆盖度较高，耕地分布在林地边缘，建筑用地分布在耕地边缘，形成了耕地包围林地，居民地与耕地交错分布的空间格局。

3.4 NDVI结构变化分析

根据研究区地表覆被类型，将NDVI影像分割为6个层级，小于0的为水体，0～0.2之间为裸地，0.2～0.4之间为稀疏草地，0.4～0.6之间为草地及湿地，0.6～0.8之间为灌丛及林地，0.8～1.0之间为优质林地，结果如图12所示。将密度分割统计数据进行百分比制图(图12(a))，可以看到：整个研究区灌丛及林地(0.6

(a)总体状况 (b)核心区

图12 1988—2016年徂徕山—莲花山NDVI密度分割百分比

Fig.12ThepercentageofNDVIdensitydivisioninCulai-Lianhuabasinduring1988—2016

4 讨论与结论

4.1 讨论

徂徕山—莲花山地区受自然环境演变和人类活动的双重影响，其生态环境的变化可通过遥感植被指数及景观格局的变化做出评价。

1)自然环境演变是漫长的，影响该区1982—2012年的31 a间自然演变的主要因素是气候变化。徂徕山—莲花山地区表现为暖湿化特征；气温呈波动上升趋势，增速为+0.046 ℃/10 a，降水量呈缓慢上升趋势，速度为+43.0 mm/10 a。年均气温、降水变化起伏较大，与耕地、林地、建筑用地面积变化无显著的对应关系。

2)遥感指数能较好地反映该区域的生态环境的变化。1988—2016年的29 a间核心区和核心区周围地段NDVI均值在2003年为相对低值，在2009年达到相对高值。同时2003年NDBI为相对高值，在2014年为相对低值。29 a间研究区核心区NDVI均值整体上升，NDBI均值下降，核心区周围NDVI均值下降，NDBI均值上升，NDVI与NDBI呈显著负相关关系，人类活动影响着区域植被指数的变化。

3)21世纪，随着经济的发展，城镇化建设加强，保护环境、可持续发展战略的实施，林地和建筑用地面积显著增加，人类活动对生态环境的影响程度加深，建筑用地面积的迅速增加，增加了耕地转化为建筑用地和林地的速率，灌丛、林地和草地向稀疏草地转变，人类活动已触发了生态格局的转变。

4)把1982—2016 年的35 a间获得的人文与自然要素：NDVI均值、NDBI均值导入SPSSⓒ，两两做相关性分析，可得到Pearson相关系数阵如图13所示。图中**表示在0.01水平(双侧)上显著相关；C表示核心区(core area)；A表示核心区周围地段(around the core area)。

图13 人文与自然要素相关性矩阵Fig.13 The correlation matrix between humanity and natural elements

PCNDVI-CNDBI=-0.946(在0.01水平(双侧)上显著相关，样本数n=7，缺1982年数据(下同))；PANDVI-ANDBI=-0.955(在0.01水平(双侧)上显著相关，n=7)；PCNDBI-ANDVI=-0.930(在0.01水平(双侧)上显著相关，n=7)；PCNDVI-ANDBI=0.988(在0.01水平(双侧)上显著相关，n=7)；PCNDBI-ANDBI=0.953(在0.01水平(双侧)上显著相关，n=7)；PCNDVI-ANDVI=0.988(在0.01水平(双侧)上显著相关，n=7)。另外，P林地面积—耕地面积=-0.787(在0.05(双侧)上显著相关，n=8，含1982年数据)。

分析研究区NDVI与NDBI变化得知，核心区NDVI与核心区NDBI呈负相关；核心区周围NDVI与核心区周围NDBI呈负相关，整个研究区出现明显的生态界限。

分析时间序列的林地面积变化的影响因素可以看出，在1982—2016年的35 a间，经济迅速发展，人类活动加强，旅游开发行业发展，退耕还林政策的实施，林地面积受耕地面积变化的影响，两者呈明显负相关，与气候因素不显著相关，土地政策大于气候变化对LUCC的影响，人类活动已成为生态环境变化的主要影响因素。

4.2 结论

1)徂徕山—莲花山地区降水量对全球气候变化响应敏感，1982—2012年的31 a间平均气温增速为+0.046 ℃/10 a，降水量呈缓慢上升趋势，速度为+43.0 mm/10 a。该区域NDVI均值、林地面积与气候变化因素不显著相关，即土地覆被、土地利用面积变化与气候变化因素不显著相关。

2)1982—2016年的35 a间徂徕山—莲花山地区人类活动显著增加，建筑用地面积增加迅速，增速为6.253 km2/a、受退耕还林政策的影响，耕地面积下降，林地面积增加，部分耕地转为建筑用地和林地，且林地面积与耕地面积呈明显负相关关系，人类活动已成为影响土地利用变化的主要因素。

3)徂徕山-莲花山核心区与核心区周围出现了明显的生态界限，1988—2016年的29 a间核心区NDVI增速为0.006 9/10 a，NDBI均值降速-0.014/10 a；核心区周围NDVI降速为-0.018/10 a，NDBI均值增速为0.003 5/10 a。形成了建筑用地、耕地包围林地，林地中分散着耕地和居民地的格局。

4)随着徂徕山—莲花山地区被开辟为风景区，旅游业发展迅速，游客的增多带来了经济的发展的同时，生态环境受到了显著影响；虽然退耕还林政策改善了交错区的生态状况，仍形成了明显的生态界限，不利于该区生态环境的可持续发展。

5)本文的不足之处在于数据不充分造成了某些指标波动性较大；由于云量的影响，部分年份的NDVI没能计算；研究区相同时相的数据难以获得，部分数据成像时间差距较大，如1988和2013年2期数据获取时间差了53 d，对植被覆盖度的统计影响较显著。

志谢：本文系对蔡盼丽本科毕业论文的提炼，该毕业论文曾获评河南理工大学校级优秀论文(指导教师：马超，mac@hpu.edu.cn)。感谢张文豪同学(同济大学2019级研究生)对本文的英文摘要提供的润色。