浙江省丽水市亚热带森林景观格局对森林碳密度的影响

王剑武，季碧勇，王铮屹，朱程昊

（1.浙江省森林资源监测中心，浙江 杭州 310020；2.国家林业和草原局 华东调查规划院，浙江 杭州 310019）

气候变化是人类面临的全球性问题。到2050 年，为了将全球升温控制在1.5 ℃以内，必须增加约1 亿 hm2森林[1]。森林生态系统是全球陆地生态系统中最大的碳库[2]。森林碳密度是衡量森林生态系统固碳能力的重要指标，其影响机制越来越受到重视，已围绕树种组成[3]、气候驱动变化[4]、森林管护水平[5]、道路网络[6]和城市扩张[7]等自然与人为因素开展了富有成效的研究。

实际上，景观格局反映了自然与人为因素在不同时空尺度上作用的最终结果，并能够体现各生态过程在不同时空尺度上的相互作用关系[8]。森林景观格局反映了森林景观的组成、空间分布、数量、结构和功能等信息[9]，森林景观规模异质性的增加，改变了森林的结构、功能和过程[10]，关注森林景观异质性能更好地衡量碳循环等森林生态功能[11]。森林景观如何影响森林储量受到越来越多关注[12]，关于森林景观格局对森林碳储量的影响机制也进行了研究。森林景观格局优化可以促进森林碳储量提高[13]，扩大斑块面积从而提高森林景观连通性，促进次生林碳储量提升[14]，而森林破碎化则会减少森林边缘的碳储量[15]。张丹[16]研究了不同尺度城市森林景观格局与碳储量的相关性，并指出增加森林斑块面积和连结度，减小斑块的平均邻近距离与周长面积比均可增加城市森林碳储量。宋洁[17]分析了祁连山森林景观格局对森林碳储量的相关性，认为增加森林景观面积、景观形状复杂度、景观聚集度和连通性能够提升森林碳储量。吕海亮[18]研究了城市森林碳密度与景观格局关系，认为平均斑块面积与碳密度呈显著正相关，而景观形状指数则与其呈显著负相关。也有学者聚焦森林景观格局对森林碳密度的影响研究，通过耦合协调度模型探究碳密度与森林景观不同组成、分布等特征的关系[19]。然而，量化森林景观格局对森林碳密度影响的研究还不多，且以上研究虽然解析了森林景观格局对森林碳密度的相关性，但很少从空间上解释森林景观对森林碳密度的影响规律。实际上森林碳密度具有空间自相关性[20-21]，忽略地理因素的空间非平稳性，易导致结果有偏差或估计效率低下。

目前，地理加权回归模型(GWR)、空间计量模型、地理探测器等方法已被用于分析研究对象与各类驱动因素之间的空间关联性。其中，GWR 解释模型方差比其他模型更有效，且相比于普通最小二乘模型(OLS)，能较好地揭示解释变量的空间异质性，但GWR 模型是基于恒定带宽，忽略了解释变量尺度作用差异。为此，FOTHERINGHAM 等[22]提出了多尺度地理加权(MGWR)模型，解决了不同解释变量在不同尺度同一带宽的问题，即允许在多个空间尺度上建立自变量和因变量的关系模型。MGWR 模型已成功应用于生境质量[23]、生态效率[24]、森林火灾[25]等的驱动因素分析，但在森林景观研究中的应用仍处于探索阶段。

本研究选取了浙江省丽水市为研究区域，采用MGWR 模型探究亚热带森林景观格局对森林碳密度的影响机制和空间非平稳性，为研究区域森林碳密度时空特征提供框架，以期为制定国土绿化、森林质量提升政策，因地制宜开展森林保护与管理提供依据。

1 研究方法

1.1 研究区概况

丽水市位于浙江省西南部，27°25′～28°57′N 和118°41′～120°26′E，地貌以丘陵、中山为主，市域面积为1.73 万km2。丽水市属中亚热带季风气候，温暖湿润，雨量充沛，具有典型的山地气候。作为南方亚热带重要的集体林区，及浙江省森林资源最为丰富的设区市，截至2019 年，丽水市森林面积为142.14 万hm2，森林覆盖率为82.27%，森林蓄积量为0.96 亿m3，居浙江省前列。近年来，随着快速的城市化和社会发展，该地区土地利用发生了相当大的变化，这对森林提供生态系统服务能力产生了重要影响[26]。同时该地区国家公园和自然保护地建设，以及森林管理活动也对森林碳循环产生了重大影响[27]，增强了森林景观格局和森林碳密度的空间异质性。

1.2 数据来源

①土地利用矢量数据来源于浙江省森林资源年度变更成果。将研究区的土地利用分为乔木林地、竹林地、灌木林地、未成林造林地、苗圃地、迹地、宜林地和非林地等8 种类型。为更细致刻画森林景观特征，进一步把乔木林地细分为针叶林、阔叶林、针阔混交林等3 种类型。地类和树种分类按照DB33/T 640—2017《森林资源规划设计调查规程》[28]执行。②样地数据来源于丽水市716 个森林资源连续清查样地(以下简称连清样地)的现场调查，全市以公里格网4 km×6 km 间距布设样地，单个样地面积为0.08 hm2，形状为正方形。③数字高程模型来源于地理空间数据云(http://www.gscloud.cn/)，分辨率为30 m。④年降水量和年平均气温数据根据资源环境科学与数据中心下载的中国756 个气象站观测数据，通过普通Kriging 法进行空间插值形成栅格数据，分辨率为30 m。以上数据均取2012 和2019 年的数据。

1.3 测算森林碳密度

①利用2012 和2019 年连清样地实测的样木数据，将每株样木树种归类至表1 所列的树种类型。②使用树高-胸径曲线模型计算单株树(竹)高；③使用冠长模型计算单株树(竹)冠长；④使用单株立木(竹)生物量模型计算各样地内的单株活立木(竹)的地上和地下部分生物量，再通过与含碳率的乘积测算单株活立木(竹)的碳储量；⑤将样地内各株活立木(竹)的碳储量累加形成样地林分水平的碳储量，然后除以样地面积得到森林碳密度。⑥利用ArcMap 10.8 的球面半变异函数模型的普通Kriging 法，形成丽水市森林碳密度空间分布图。树高-胸径曲线模型、冠长模型、生物量模型见表1，含碳率见表2。

表1 树高-胸径曲线模型、冠长模型和生物量模型Table 1 Height-diameter curves model, crown length model and biomass model

表2 树种含碳率Table 2 Carbon contents of species

为验证碳密度空间分布结果精度可靠性，基于已有数据，先将空间分布结果由栅格转为矢量面要素的格网，并裁切提取出丽水市行政范围内的格网；再将每个格网的面积与其碳密度相乘并汇总加和得到丽水市森林碳储量；然后按照文献[29]所述的系统抽样统计方法，基于连清样地计算丽水市森林碳储量的估计中值、估计精度和估计区间；最后根据文献[30]，基于空间分布结果计算的森林碳储量如落入估计区间，则说明该结果是有精度保证的，森林碳储量是碳密度乘以丽水市面积得到，也可说明森林碳密度空间分布结果是有精度保证的。

1.4 量化森林景观格局指数

选择5 个常用且相关性较低的景观指标：香农多样性指数(SHDI)、最大斑块指数(LPI)、斑块密度(PD)、蔓延度指数(CONTAG)、景观形状指数(LSI)。各指标的含义和计算公式见文献[31]。另外，本研究选择了年平均降水量(PREP)、年平均气温(TEMP)、海拔(ELEV)作为控制变量。气温和降水量是影响森林固碳能力最常见的气候要素，其影响是双向的。气温会影响植被光合作用和净初级生产力，从而影响森林生长[32]，但高温引起的呼吸增加可能使碳储量随着温度升高而减少[33]。一定范围内的水资源会促进森林的生长，水分过多则会抑制土壤呼吸，影响植物生长[34]。使用Fragstats 4.2 软件，参照已有对市级尺度的景观研究[35-36]，选取5 km×5 km 为窗口大小，采用移动窗口法生成各景观指标栅格图，基于ArcMap 10.8 生成4 km×4 km 的格网(研究区共计1 204 个网格)，对每个格网的景观格局指标数值和控制变量数值分区统计。

1.5 空间自相关分析

地理学第一定律表明，在涉及具有地理空间的数据时，事物之间距离越近关联性越大[37]。采用Global Moran’sI和Local Moran’sI来描述森林碳密度的全域和局部聚类特征。计算公式可见文献[38]。Moran’sI介于-1 到1 之间，＞0 表示存在正相关性，＜0 表示存在负相关性，0 表明不存在空间自相关。使用GeoDa 1.12 绘制LISA (local indicators of spatial association)图，反映森林碳密度的空间格局，LISA 图包括4 个类别，“高—高”和“低—低”表示空间正相关，“高—低”和“低—高”表示空间负相关。Global Moran’sI也将用于量化MGWR 模型的残差空间自相关，如果残差存在显著的空间自相关，则表明模型中缺失了关键解释变量。

1.6 MGWR 模型

采用MGWR 模型探究森林景观格局对森林碳密度的多尺度影响。MGWR 模型通过不断寻找各解释变量的最优带宽并不断更新参数估计值，直至系数迭代收敛，从而形成最终回归结果[39]。该模型如下：

其中，yi是格网i的被解释变量；xij是格网i的第j个解释变量；βbwj(ui,vi)为格网i的第j个解释变量的局部回归系数，(ui,vi)为格网i的空间位置，bwj为第j个解释变量回归系数所使用的带宽bw；m为解释变量的个数；εi为误差项。本研究采用MGWR 2.2 进行模型构建，MGWR 形式选择Gaussian，空间Kernel 类型选择Adaptive Bisquare 函数，带宽搜寻方法采用Golden Section，模型最优化准则采用AICc 准则。

2 结果与分析

2.1 森林碳密度空间格局

2012—2019 年，研究区森林碳密度平均值从23.19 t·hm-2提高到31.96 t·hm-2，年均增量为1.25 t·hm-2。根据研究区森林资源连续清查结果，森林面积和森林蓄积量分别从2012 年的132.81 万hm2和0.65 亿m3增加到142.14 万hm2和0.96 亿m3，年均分别增加1.33 万hm2和0.04 亿m3，森林面积和森林蓄积量逐年上升促进了森林碳密度增加。森林碳密度空间分布格局存在显著差异，低森林碳密度主要在研究区的北部、西部和西南部的森林碳密度较高(图1)。2012 和2019 年研究区的Global Moran’I分别为0.889 和0.891 (P＜0.01)，表明森林碳密度较高或较低区域趋于聚集。由图2 发现：①“高—高”区域主要分布在西南部和中南部地区。西南部是瓯江源头，拥有钱江源-百山祖国家公园百山祖园区等多处自然保护地，得益于较少的人为干扰、公益林保护和森林生态保护修复，森林质量高。②“低—低”区域主要分布在北部，北部地区包含松谷平原和碧湖平原，这些地区耕地、建筑用地较多，受人类活动影响频繁。③“低—高”区域分布极少，说明出现周围森林碳密度高而本地森林碳密度低的空间结构概率较低。

图1 丽水市森林碳密度空间分布图Figure 1 Spatial distribution map of forest carbon density in Lishui City

图2 丽水市森林碳密度LISA 分布特征图Figure 2 Spatiotemporal characteristics of LISA map of forest carbon density in Lishui City

基于丽水市连清样地，系统抽样统计结果(表3)表明：2012 年全市森林碳储量估计中值为4 023.48 万t，估计区间为(3 729.05～4 317.90)万t，估计精度为92.68% (P＜0.05)。2019 年全市森林碳储量估计中值为5 541.19 万t，估计区间为(5 192.46～5 889.91)万t，估计精度93.71%为(P＜0.05)。2012 和2019 年基于空间分布结果计算的森林碳储量分别为4 023.08 万t和5 543.42 万t，均位于相应年份的估计区间内，说明森林碳密度空间分布结果是有精度保证的。

表3 分别基于连清样地、森林碳密度空间分布结果计算的丽水市森林碳储量Table 3 Forest carbon storage in Lishui City based on the continuous forest inventory and spatial distribution of forest carbon density, respectively

2.2 基于MGWR 模型的空间非平稳性分析

2.2.1 变量尺度效应 由表4 可见：2012 和2019 年的MGWR 模型的Adj-R2分别达0.898 和0.892。模型的带宽包含44、79、992 和1 204，由于研究区被划分为1 204 个网格，因此从行政边界角度看，带宽992 和1 204 的空间尺度可定义为市级，44 和79 可定义为县级。LSI 带宽为44，作用尺度较小，表明景观形状复杂度对森林碳密度的影响在空间上存在较大差异。CONTAG、LPI 和SHDI 的带宽都较大，说明森林景观连通性、面积优势度和景观多样性对森林碳密度影响的空间差异相对较小，也表明上述3 个影响因素回归系数空间异质性较不明显。PD 带宽从79 升至1 204，作用尺度在扩大，意味着景观破碎化程度对森林碳密度影响的空间差异在变小。

2.2.2 回归结果 由表4 可见：从景观格局指标的回归参数2 期均值的绝对值看，CONTAG 对森林碳密度影响最大，其次是SHDI 和PD，LPI 和LSI 相对影响较小。从图2 的局部R2看，研究区大部分区域拟合优度较高，模型可以揭示研究区2012 年森林碳密度总变化的17.8%～95.9%和2019 年的22.1%～96.6%。从各影响因素beta 系数的正负值占比看，除CONTAG 和SHDI 外，其余景观格局的参数估计值正负情况都存在，即对森林碳密度表现不同效应影响，说明MGWR 模型能将系数的局部特征呈现出来。

2.2.3 回归结果的空间非平稳性 由图3 可见：2012 和2019 年，CONTAG 对森林碳密度的负向影响均从研究区的西北向东南下降，但对森林碳密度的负向驱动力在研究期间逐渐减弱。这可能是因为西北部包含了松古平原，同时遂昌县境内分布有若干小盆地，斜坡平缓，茶叶等经济林种植面积广，应当通过加强农田林网建设和平原绿化增加片状森林面积，降低经济树种的景观连通度，提高区域森林碳密度。PD 系数在研究区范围内以负值为主，并且负值区域在逐渐增多，说明降低森林破碎化程度一定程度上可减少对森林碳密度的负向影响。LPI 系数以负值为主，意味着大部分区域需要减少人为干扰，维持植被自然生长状态。从2012 年到 2019 年，LPI 系数正值区域有所增加，在西北部和中西部等区域与森林碳密度呈正相关关系，说明这些地方通过适当增加森林抚育、补植造林等措施，能够对提高森林碳密度产生积极影响。从2012 年到2019 年，LSI 对森林碳密度的驱动力逐渐增强。说明在一定程度上，随着斑块形状的复杂度增加，森林碳密度也会随之增加。2019 年，LSI 对森林碳密度的正向影响由西向东下降，这表明相对于研究区东部，西部森林斑块复杂程度的增加对森林碳密度的正向影响更大。这可能是因为森林景观斑块形状复杂，能够促进斑块边缘的植被形成高效率的能量共生网络[40]，增强森林的边缘效应，而斑块边缘通常具有较高的初级生产力[41]。从2012 年至2019 年，SHDI 对森林碳密度的负向驱动力逐渐减弱，形成负向影响由西向东增加的格局，说明研究区的林分种类分布比较均匀。年平均降水量系数总体上呈现由东部大于西部转变为西部大于东部的趋势，但2 期正向与负向影响区域的数量差异较小。年平均气温对森林碳密度的正向影响区域逐渐增多，而海拔系数空间分布差异不大。2019 年，年平均降水量主要在西部和西南部与森林碳密度呈正相关关系，其余区域主要呈负相关关系。年平均气温系数空间分布与年平均降水量有一定相似，主要在西南部和西北部与森林碳密度呈正相关关系。研究区的绝大部分区域的海拔与森林碳密度呈正相关关系，并且中部区域的正相关关系更强。年平均降水量、年平均气温系数正值分布区域包含了钱江源-百山祖国家公园百山祖园区、九龙山国家级自然保护区、括苍山省级森林公园，充分利用丰富的水热资源，加强自然保护地内森林生态系统保护修复，有助于提高森林碳密度。

图3 MGWR 模型系数和局部R2 空间格局Figure 3 Spatial pattern of MGWR model coefficient and local R2

3 讨论

3.1 MGWR 模型能反映森林碳密度对森林景观格局响应的空间分异

本研究发现：森林碳密度具有显著的空间自相关性，考虑影响因素的空间尺度对模型的结果会产生重大影响。MGWR 模型拟合结果表明，斑块密度、最大斑块指数、年平均降水量、年平均气温和海拔等变量在不同区域的回归系数既有正值也有负值，因MGWR 模型使用自适应带宽，可以更精确地探索驱动因素对森林碳密度的作用机制和空间异质性影响，并明确各驱动因素的具体影响尺度。

3.2 森林景观格局驱动森林碳密度对森林管理的建议

森林景观格局对森林碳密度影响表现出不同大小的驱动力和多尺度空间特征，不同尺度需要关注的影响因素是不同的，因此林业管理者可以根据这些特征制定宏观与微观的政策。在市级层面，林业管理者需要加强农田林网和生态廊道建设，降低纯林的景观连通度，同时对生态保护修复进行全区域尺度的管理，将森林恢复和减少森林破碎化战略纳入空间规划中，减少山区森林生态系统中人为干扰和森林破碎化，以提高森林碳密度的整体水平。在县级尺度上，通过合理配置森林景观组成可以有效增强森林固碳量[42]，并发挥森林组成成分间相互作用效能，发挥森林边缘效应优势。

3.3 局限性和未来研究方向

区域森林碳密度计算是一项复杂的研究任务且有多种计算方法，如采用CASA 模型计算植被净初级生产力(NPP)，再将其转换为固碳量，或采用InVEST 模型计算固碳量。森林碳密度受各种因素的影响，本研究是在系统布设的固定样地计算结果基础上，采用空间插值形成区域森林碳密度，计算方法在某些方面可能仍需改进。由于缺少灌木林实测数据，本研究未测算样地内的灌木林碳储量。随着尺度增加，森林碳密度的影响因素也愈加复杂，由于景观格局具有尺度依赖性，MGWR 模型作为线性回归模型，对愈加复杂的影响因素解释能力需进一步探讨。此外，本研究使用边长4 km 正方形的网格尺度来创建景观格局指标，并探究其与森林碳密度的关系，这个尺度可能更适用于建模，对实际森林管理会有不便。在未来研究中可以以乡镇、村级行政单位为建模单元，进一步增强模型结果的实际应用性。鉴于森林碳密度的空间自相关性，可以考虑使用空间计量模型进一步探究森林碳密度对森林经营管理和经济社会发展的响应机制。

4 结论

本研究以浙江省丽水市为研究区，使用MGWR 模型探讨了森林景观格局对森林碳密度的影响。2012 和2019 年，森林碳密度的空间自相关度分别为0.889 和0.891，呈显著的空间自相关性。森林景观格局对森林碳密度的影响在空间上存在非平稳性，CONTAG 和SHDI 存在负向影响，PD、LPI、LSI 在空间上正向和负向影响并存。森林景观格局也表现出不同尺度效应，CONTAG、LPI、SHDI 的作用尺度较大，LSI 作用尺度较小。基于研究结果，建议进一步重视森林经营管理的空间决策支持，因地制宜实施自然演替和人工干预相结合的森林生态系统保护修复措施，促进森林生态系统碳汇功能提升。