基于连清数据的广州市1979—2012年森林生态功能的变化

钱春花，赵邑晨，李明阳，李 陶

(1.南京林业大学林学院,江苏 南京 210037；2.苏州农业职业技术学院，江苏 苏州 215008；3.沈阳市勘察测绘研究院有限责任公司 辽宁 沈阳 110004)

森林生态功能是指森林生态系统及其生态过程所形成的有利于人类生存和发展的生态环境条件与效用，包含水源涵养功能、水土保持功能、气候调节功能、环境净化功能、生物多样性保护功能等[1-2]。森林作为陆地最大的生态系统，其生态功能的变化直接影响着人类生存环境的变化。国内外对生态功能评价的研究始于20世纪50年代，90年代后迅速增多[3-4]。我国对森林生态功能评价相关的研究多集中于森林生态系统对降水的分配与拦截效应、净化水质功能、保持水土功能、森林生产力、生物多样性保护等林分或生态系统尺度的单项生态功能方面[5-7]。国家、区域尺度的森林生态服务综合功能评价，通常以森林资源统计报表数据为主要信息源，研究内容缺少森林生态功能空间分析的内容，并且侧重于经济评价[8]。自20 世纪70 年代建立了森林资源连续清查(简称连清)体系以来，我国很多省份已经进行了8～9次森林资源清查。森林资源连续清查数据具有调查因子多、调查精度高、连续动态等优点。为避免以单一或少数因子进行生态功能评价带来的缺点，增强评价的综合性、权威性及可比性，2004年颁布的《国家森林资源连续清查技术规定》中提出了森林生态功能指数概念[9-10]。此后，这种基于连续数据的森林生态功能评价研究开始在区域生态功能评价中得到应用，但研究的时间跨度往往较小，同时在生态功能时空变化及其驱动因素分析方面较为薄弱。

根据第八次全国森林资源清查(2009—2013)结果，全国森林质量指数为0.60，总体处于中等质量水平，生产力低下、结构不合理、生态功能脆弱、生态稳定性差问题突出。为了提升森林质量、增强森林生态系统服务功能，2016年发布的《全国森林经营规划》(2016—2050年)提出，森林经营的目的是建立健康稳定、优质高效的森林生态系统，森林经营应该是在寻求生态合理性的基础上来追求经济效益和社会效益的合理性。由于森林经营活动对生态系统的影响及森林经营成效的体现是一个长期的过程，衡量森林经营是否可持续必须采用一个轮伐期以上的时间尺度。因此，长期森林经营规划的各种森林经营措施及林业方针政策的制定需要建立在对规划区长期、精准的森林生态功能综合评价基础之上。

广州地处经济发达的珠江三角洲，近40年快速的城市化及工业化进程、频繁的国际贸易等引发的林地占用、林地退化、外来有害生物入侵导致区域森林生态服务功能明显下降[11]。作为南方集体林区的一部分，林农是森林经营的主体，广州森林经营状况受林业政策法规变动的影响较大。在《广东省森林经营规划》(2016—2050年)中，广州属于珠三角平原丘陵生态风景林与江河防护林经营亚区，从而明确了森林生态功能的区域森林主要经营目标。在珠江三角洲，杉木、马尾松、常绿阔叶林等主要优势树种的综合轮伐期通常在30 a以上。因此，以广州市为研究区域，利用8期固定样地数据、采用生态功能指数的方法，进行1979—2012年长期森林生态功能时空变化驱动因素研究，以期为珠三角森林长期经营规划、森林经营措施及林业方针政策的制定提供科学参考。

1 材料与方法

1.1 研究区概况

广州(109°46′～117°21′E，20°08′～25°30′N)土地总面积7 434.4 km2。2019年的广州市常住人口为1 530.59万人，地区生产总值23 628.60亿元。广州市辖11区，本级统筹区7个，其中越秀、海珠、荔湾3区为老城区，番禺区、花都区、从化区、增城区为新4区。

广州地处亚热带沿海，气候属南亚热带季风海洋性气候，全年平均气温19.7～26.3 ℃，年降水量约为1 802.4 mm，平均相对湿度77%，具有温差较小、光照充足、温暖多雨及冬季短、夏季长等显著特征。广州属于丘陵地带，地势东北高、西南低，背山面海，北部是森林集中的丘陵山区，最高峰为北部从化市与龙门县交界处的天堂顶，海拔为1 210 m；东北部为中低山地，中部是丘陵盆地，南部为沿海冲积平原，是珠江三角洲的组成部分。

广州市2015年的林业用地面积29.2万hm2，其中有林地面积27.6万hm2，森林蓄积量达1 499.67万m3，森林覆盖率为42.03%。主要优势树种中常绿针叶树种有杉木(Cunninghamialanceolata)、马尾松(Pinusmassoniana)、湿地松(Pinuselliottii)，常绿阔叶树种有桉树(Eucalyptusrobusta)、黎蒴(Castanopsisfissa)、马占相思(Acaciamangium)、木荷(Schimareinw)。由于森林结构、龄级结构不合理，用材林单位面积蓄积量较低，沿江、沿河防护林质量不高。广州市生态效益低下的桉树纯林和马占相思纯林比例较大，单位面积乔木林蓄积量(55.68 m3) 低于全省平均水平(58.25 m3)。薇甘菊、松材线虫、松突圆蚧等外来林业有害生物对区域森林生态安全构成较大威胁。总体评价认为区域森林生态功能、森林景观质量不高。

1.2 数据来源与预处理

底图审图号：粤审(2018)121号。下同。图1 广州市森林资源连续清查固定样地空间分布Fig.1 Spatial distributions of fixed plots of continuous inventoryies of forest resources in Guangzhou

数据源主要有：①广州市1979—2012年8期148块森林资源连续清查固定样地数据，样地大小为0.067 hm2(图1)，样地间距为6 km×8 km，样地属性表中的郁闭度、植被覆盖度、自然度、群落结构、树种结构、林分蓄积和平均树高等因子与森林生态功能密切相关；②来源于中国科学院资源环境科学与数据中心(http://www.resdc.cn)的研究区2012年均气温和年降水量空间插值数据、2010年人口分布(简称人口密度，下同)与国内生产总值(简称人均GDP，下同)网格数据，数据的空间分辨率均为1 000 m；③源自美国国家海洋大气局国家环境信息中心(https://www.ngdc.noaa.gov)提供的2012年全球辐射定标夜间灯光数据(DMSP/OLS，简称灯光亮度，下同)，数据的空间分辨率为1 000 m；④来源于地理空间数据云平台(http://www.gscloud.cn)的2012年10月的Landsat ETM+遥感影像(轨道号122/44，云覆盖2.0%)、空间分辨率为30 m的数字高程模型(DEM)。

由于森林生物量、森林生产力等森林生态功能的高低与植被生长状况、地形条件、气候条件、社会经济状况密切相关[12]。预处理结束后，在ArcGIS平台上，利用空间分析工具箱中的Extract工具，提取每个固定样地的植被生长(植被覆盖度)、地形(海拔、坡度)、气候(年降水量、年平均温度)、社会经济(人口密度、人均GDP、灯光亮度)合计8个环境因子。其中，植被覆盖度是在计算出研究区归一化植被指数(NDVI)基础上，由二分法模型[13]生成。

1.3 生态功能评价指标构建

根据2004年的《国家森林资源连续清查技术规定》和GB/T 38590—2020《森林资源连续清查技术规程》，采用生态功能指数对研究区各个样地生态功能进行评价。从固定样地的60多个属性中，选择森林郁闭度、植被总盖度、森林自然度、群落结构、树种结构、枯枝落叶层厚度、森林蓄积和林分平均树高8个与森林生态功能密切相关的因子，按相对重要性通过加权平均构建森林生态功能指数来综合评定森林生态功能等级。各生态功能评价指标及其权重见表1。

表1 森林生态功能评价因子及类型划分标准Table 1 Evaluation factors and classification standards of forest ecological functions

在8个评价指标中，森林蓄积量、郁闭度、平均树高3个因子8期数据都有，并且属于定量因子，可以根据表1的标准直接划分等级。对于早期森林资源清查样地属性表中没有的评价指标，如自然度、群落结构、树种结构，根据样地中其他调查因子，按照2004年的《国家森林资源连续清查技术规定》和GB/T 38590—2020《森林资源连续清查技术规程》进行等级划分，然后再根据表1划分等级。如根据技术规定中的“森林自然度划分与赋值标准”(表2)，可以根据样地的地类、林种、树种、起源、林龄等其他因子判断出一个样地的森林自然度等级。

为便于计算分析，将类型Ⅰ、Ⅱ、Ⅲ分别赋值3、2、1。在进行指标综合之前，按照参考文献[14]的方法，分别大者为优(指标1、4、5、6、7)或小者为优(指标2、3、8)对各评价指标进行标准化。在此基础上，按下式计算综合得分：

(1)

式中：Y为森林生态功能指数，Xi为第i项评价因子类型的标准化得分值，Wi为第i项评价因子的权重。然后根据综合得分值评定每个样地的生态功能等级。

表2 森林自然度划分标准与赋值Table 2 Forest natural degree division standards and assignments

1.4 时间增减和空间格局变化分析

利用一元线性回归方程可以分析各个样地的生态功能时间变化趋势[15]，其计算公式为：

(2)

式中：θslope表示趋势斜率；n为研究时间段的个数(n=8)，Ci为第i年生态功能指数。θslope若为正值，则说明样地的生态功能在研究期间呈增加的趋势；θslope为0，表示研究期间样地生态功能没有发生变化；θslope为负值，则表示样地的生态功能在研究期间呈降低的趋势。

①分娩方式：统计两组产妇自然分娩率、难产率和剖宫产率。②产程进展：比较两组产妇第一产程、第二产程时间。③分娩结局：比较两组产后出血、新生儿窒息、早产、新生儿死亡发生率。其中，产后出血指胎儿娩出后2 h产妇阴道出血量>500 mL；新生儿窒息指新生儿出生后1 min Apgar评分<7分。

利用ArcGIS包含的空间统计(spatial analysis)工具，选择莫兰指数(Moran’sI，用于判断样地生态功能的空间自相关性)和高值点(热点)地理中心分析(central feature)来探讨广州市森林生态功能空间变化规律。采用ArcGIS 10.2空间统计工具箱中的聚集及特例分析工具(cluster and outlier analysis-anselin local moran’sI)，通过空间聚类，识别固定样地中的森林生态功能高值点(空间热点)[16]。空间热点地理中心分析的目的在于判断森林生态功能高值点样地的空间变化趋势。通过空间聚类，识别出固定样地中的森林生态功能高值点(空间热点)，采用ArcGIS空间统计工具箱中的Central Feature工具生成研究区1979—2012年的8期森林生态功能高值点地理中心，并与2012年研究区Landsat ETM假彩色合成叠加。莫兰指数是一种最常用的自相关系数，计算公式如下：

(3)

1.5 地理加权回归

地理加权回归(geographically weighted regression,GWR)模型加入了观测点的空间位置信息，即在每个观测点都进行一次独立的线性回归[17]，因而可以解释每个变量在各个区域影响的差异性，从而克服了传统线性回归的缺点。本研究采用GWR构建广州市森林生态功能空间差异驱动因素回归模型，计算公式为：

yi=β0(ui-vi)+∑βk(ui,vi)xik+εi。

(4)

式中：β为模型系数,(ui，vi)是第i个固定样地的空间位置坐标，βk(ui，vi)是连续函数在固定样地i的生态功能指数值，xik是变量xi在固定样地i的值，εi代表随机误差。以2012年为例，选择148块固定样地森林生态功能指数为被解释变量，选取人口密度、人均GDP、灯光亮度、年降水量、植被覆盖度、海拔、坡度共7项指标作为解释变量建立GWR模型进行研究区生态功能空间差异性驱动因素分析。

2 结果与分析

2.1 广州森林生态功能时间变化趋势分析

在计算148块样地8期生态功能指数的基础上，将生态功能指数值代入公式(2)，计算出每块样地1979—2012年的生态功能趋势斜率，以及斜率的平均值、标准差(σ)。按照均值-2σ、均值-1σ、均值+1σ、均值+2σ的方法设定4个阈值，将1979—2012年研究区148块固定样地的生态功能时间变化趋势分为5种类型：显著变高(A)、轻微变高(B)、基本稳定(C)、轻微变低(D)、显著变低(E)，并计算5种类型样地的占比(图2A)。

图2 1979—2012年广州市不同森林生态功能变化的类型样地比例及等级比例Fig.2 The proportion of different plots types of forest ecological function changes and function levels in Guangzhou from 1979 to 2012

从图2A可以看出，1979—2012年广州森林生态功能指数提高的样地数多于变低的样地数，森林生态功能整体向变好的方向发展。森林生态功能指数提高(轻微变高和显著变高)的样地数是变低(轻微变低和显著变低)样地数的1.82倍，显著变高样地数占研究区样地总数的19.59%，轻微变高样地数占比14.86%；轻微变低样地数占样地总数的9.47%，显著变低样地数仅占9.46%，森林生态功能指数基本稳定的样地数占46.62%。

在计算出研究区1979—2012年148块样地8期森林生态功能指数前提下，按照低(<0.3)、中等(0.3～0.6)、高(>0.6)的标准，统计各时段各个等级样地的比例变化(图2B)。从图2B可以看出，研究区1979—2012年森林生态功能指数等级为高、中等的样地数整体呈现曲折上升、平稳增加的趋势，森林生态功能指数等级为“低”的样地数呈阶梯状下降且总体减少趋势。

研究区森林生态功能呈现总体提升的良好态势，与20世纪80年代中后期开始实施的“造林灭荒、退耕还林、天然林保护、城镇绿化”[11]等生态工程有着密不可分的关系。从图2B可以看出，1988年，研究区生态功能高的样地比例出现较大幅度下降，这与1983年山林分山到户引发的过度砍伐造成森林资源质量下降有关；1992年研究区生态功能中、高等样地比例出现一定幅度上升，生态功能低的样地比例出现大幅度下降，这种变化源自于1985—1991年的广州市绿化达标、造林灭荒的生态工程建设成果。

2.2 广州森林生态功能空间变化分析

图3 1979—2012年广州市森林生态功能空间自相关性变化Fig.3 Changes in a spatial autocorrelation of forest ecological functions in Guangzhou from 1979to 2012

基于8期固定样地生态功能指数Moran’sI的研究区生态功能空间自相关随时间变化见图3。从图3可以看出，148块8期固定样地森林生态功能指数的Moran’sI值均>0，说明1979—2012年的研究区森林生态功能存在较为明显的空间正相关。1979—2012年的森林生态功能指数Moran’sI值变化呈现先增大再迅速下降，后又逐渐增大的趋势。森林生态功能指数Moran’sI的最大值为1988年的0.265，最小值为1992年的0.056。Moran’sI值的变化趋势揭示了研究区森林生态功能空间聚集性增大后急剧下降、又缓慢增强的复杂变化历程。这是由于1979年后，伴随着“林业三定”政策的实施，虽然研究区低海拔地区森林资源受到了严重破坏，但是广州市东北部高海拔、高坡度地区的森林由于交通不便得到了较好保护，表现为森林生态功能较高，故1979—1988年期间研究区森林生态功能指数空间自相关性呈明显增强趋势。1985年广东省采取“五年种树、十年绿化广东”的举措，使得研究区森林覆盖率较低的平原、丘陵地区森林面积增加，森林质量提高，平原、丘陵山区的森林生态功能差距持续减小，形成1988—1992年研究区森林生态功能的空间自相关性呈现较大幅度下降趋势。1993—1997年以开发区建设、房地产投资、股票和金融为代表的经济过热现象，导致研究区低海拔地区的森林被采伐、林地被蚕食的面积较大，森林生态功能空间自相关性再次呈现明显增强趋势。1997—2002年，随着经济治理整顿措施的落实，森林资源严重破坏的趋势得到遏制。2003年《中共中央国务院关于加快林业发展的决定》(www.moa.gov.cn/gk/zcfg/xzfg/200301/t20060123-541448.htm)确立了以生态建设为主的林业发展方向之后，研究区森林的干扰和破坏现象明显减少，森林生态功能的空间聚集性逐渐趋于稳定。

将1979—2012年的8期森林生态功能高值点地理中心与2012年研究区Landsat ETM假彩色合成叠加(图4)可以看出，1979—2012年广州市森林生态功能高值点地理分布中心主要集中在广州市中部偏北地区，结合中心局部放大图可以看出，研究区森林生态功能高值点地理分布中心呈现由南向北、由西向东逐渐偏移的变化趋势。

图4 1979—2012年广州市森林生态功能高值点地理中心的变动Fig.4 Changes in geographic centers of high-value points of forest ecological functions in Guang-zhou from 1979 to 2012

这种生态功能地理分布中心的变动趋势与研究区的地形条件差异、经济发展水平的高低、林业方针政策的变化存在密切关系。1978年后，研究区开始经历较为快速的城镇化、工业化进程。研究区的北部海拔高、坡度陡、交通不便，受人为干扰少，森林覆盖率高；中部、南部的丘陵盆地，森林覆盖率低，城镇化、工业化引发的林地被占用、外来生物入侵较为严重。因此，1979—1997年，随着时间的推移，研究区中部、南部的森林质量呈下降趋势，分布中心向北偏移；2000年后，随着退耕还林、天然林保护、美丽乡村建设的进行，研究区东部少林地区的森林覆盖率大幅度提高，森林质量呈现增加趋势；2000—2012年分布中心则呈现向东偏移的趋势。

2.3 广州生态功能空间差异性驱动因素分析

1979—2012年间广州市森林生态功能空间分布呈显著的空间正相关特征，表明采用普通线性回归最小二乘法对广州市森林生态功能空间差异性进行驱动因素分析，得到的结果有可能存在较大偏差。因此，采用地理加权回归方法，进行148块固定样地生态功能空间差异性驱动因素分析。

首先利用统计分析软件SPSS 23.0对8个驱动因素进行多重共线性分析，当容差≤0.1或方差膨胀因子(VIF，是容差的倒数)≥10时说明自变量间存在严重的共线性情况。年均气温与年降水量两因子的容差均≤0.1且VIF均≥10，证明二者存在共线性。因此，地理加权回归选取解释变量时删除年均气温因子，保留对森林生长较为敏感的年降水量因子。基于多重共线性分析、GWR模型计算结果发现，在输出的评价系数中，Cond最小值为16.675，最大为27.448，平均值为22.062，均小于30，表明模型模拟结果比较理想。模型的局部决定系数(R2)为0.649，说明模型可以解释64.9%的差异，拟合效果较好(表3)。由表3可知，海拔的局部平均回归系数均值最大，其次是植被覆盖度，第三是坡度，这3个环境变量与因变量均存在着较为显著的正相关关系。排在第4位的人均GDP对森林生态功能指数有负向影响。灯光亮度、人口密度、年降水量与因变量弱正相关。按照相关系数绝对值从大到小的顺序排列，影响研究区森林生态功能空间差异的前4个主要驱动因素为海拔>植被覆盖度>坡度>人均GDP。将7个环境因子分为植被、地形、气候、社会经济四大类，根据每类因子局部平均回归系数绝对值的均值排序，影响研究区森林生态功能的驱动因素按照重要性大小排序如下：植被生长>地形影响>社会经济状况>气候条件。

植被覆盖度的高低反映了林木生长、森林经营水平的高低。林木生长状况越好、森林经营水平越高，森林质量越高，森林生态服务功能越强。在区域尺度上，海拔越高、坡度越陡的地方，由于交通不便，过度砍伐、林地占用等人为不良干扰的频率、强度均较小，森林生态功能可以得到较好的保护。

表3 2012年广州市生态功能驱动因素地理加权回归系数表Table 3 Regression results of GWR for driving factors of ecological functions in 2012 in Gungzhou

在经济发达的珠三角地区，人均GDP的增加往往伴随着城镇化和工业化引发的林地占用、林地退化及过度砍伐，因此该因子与森林生态功能负相关；广州地处南亚热带，属于南方集体林区，较为充足的降水量可以促进林木生长；一定的人口密度与经济发展水平，有利于森林集约经营及森林质量提升。因此灯光亮度、人口密度、年降水量3个因子与森林生态功能存在比较弱的正相关关系。

3 讨 论

1979—2012年，引发广州森林生态功能时空变化的主要驱动因素是林业方针政策的变化、城镇化和工业化进程导致的人类干扰。1983年分山到户引发的过度砍伐、1993年经济过热引发的林地占用以及森林退化，是区域森林生态功能质量下降、空间格局变化的主要驱动因素。1985年发起的造林灭荒工程、2000年后的退耕还林及天然林保护与美丽乡村建设等生态工程，则较大程度上提升了区域森林生态功能，引起了森林生态功能高值点地理分布中心发生变化。总体而言，1978—2012年，研究区森林生态质量呈现上升趋势，森林生态功能的空间聚集性逐渐趋于稳定。研究结果与戴前石等[18]对湖南省森林生物量长期变化案例研究结果相似。

森林生态功能空间差异性GWR分析表明，在影响区域森林生态功能异质性的因素当中，植被生长>地形影响>社会经济状况>气候条件。由于地形因子难以改变，在短的时空尺度上，气候因子人为控制的难度较大。增大人力、物力、科技投入进行集约经营，通过精准提高森林质量，增强区域森林生态服务功能；与此同时，在发展经济的同时，加强林地保护力度，防止林地占用、林地退化，是增强区域森林生态服务功能的必由之路。

本研究的重点是广州市长期森林生态功能时空变化及驱动因素分析，因此生态功能指标及权重直接采用了2014年的《国家森林资源连续清查技术规定》和GB/T 38590—2020《森林资源连续清查技术规程》。虽然这些指标的选取及权重的确定经过了国家林业部门权威专家的评审和论证，但是否完全适用于研究区还需要进一步分析。此外，人均GDP、人口密度、灯光亮度、年均温度、年均降水量这些社会经济、气候因子空间分辨率较低，有的因子获取时间与样地调查时间不完全一致。由于固定样地间距较大、研究区空间范围较广，一定程度上抵消了对研究结果准确性的影响，但具体影响尚需进一步的量化分析。