基于GWR思茅松单木含碳量模型的构建及空间分布研究

刘胜刚余哲修欧光龙刘 畅

（1.西南林业大学林学院，云南昆明 650233；2.东营市河口区河庆路消防救援站，山东东营 257000；3.北京林业大学林学院，北京 100083）

森林植被在陆地生态系统中占据主要地位，其对全球碳循环起着至关重要的作用，以二氧化碳为主的温室气体的排放大幅上升，从而导致生态平衡遭到破坏，因此对森林碳储量的研究已成为国内外学者的热点[1]。

目前学者对森林碳储量的相关研究方法归纳起来主要有样地清查法、遥感估测法和模型估测法[2]。其中模型估测森林碳储量是目前较为流行且预估精度较高的方法之一[3-4]。有研究表明林木在生长期间会与周围林木、其他树种分布和森林干扰等多种因素相互作用相互影响，这种林木间的影响产生了空间效应[5]，即空间自相关性和空间异质性。而单木的生长同样会受到自身器官、邻近植物和土壤等环境因素的影响，存在空间效应。当前林木间空间效应的研究较多集中在区域大尺度森林生物量上[6-8]，在单木不同器官含碳量上的研究较为少见。传统的统计学假设林木数据相互独立与地理位置无关，因此采用传统统计学模型对森林碳储量的估测是有偏差的。为了解决空间效应对森林采样数据的影响，可以采用空间回归模型的方法，目前应用在森林碳储量估测的空间模型很多，比如空间误差模型、空间滞后模型、线性混合模型和地理加权回归模型等[9-12]。空间模型能够解决采样数据分布间存在的空间效应问题，但不同种类的空间模型适用条件各有不同。如空间滞后模型仅考虑了全局模型中没有考虑的空间自相关，而忽略了空间异质性的影响，该类模型虽然降低了自相关对模型的影响，提升了模型拟合效果和模型精度，但由于空间异质性存在，模型估测仍然存在一定偏差。Fotheringham等[12]为了更好的利用回归模型中的空间特性分析数据，在空间变异系数回归模型的基础上，提出地理加权回归模型，它能直观显示空间数据间非平稳性，因此该模型被认为是解决空间异质性问题时最有效的方法之一[13]，目前广泛应用于经济、生态和林业等领域。

思茅松（Pinuskesiyavar.langbianensis）具有树干端直高大、生长速度快、木材用途广、生态效益好和经济效益高等特点，是云南省重要的造林树种。目前对思茅松生物量及碳储量的相关研究已获得不少成果，但是对于单木不同器官的固碳能力研究相对较少。单木不同器官的固碳量是否存在空间效应，又分别对全树及样地的贡献力达到多少，都会影响森林经营决策和措施。综上所述，本研究以云南省澜沧县思茅松天然林为研究对象，对其空间效应进行分析的同时，采用测树学因子和环境因子构建单木树干、树枝、树叶、树根和全树含碳量的地理加权回归（GWR）和最小二乘回归模型（OLS），并分析模型优劣，这对及时准确地获取区域森林碳储量的分布信息以及帮助管理者更好的对森林进行经营管理具有重要意义。

1 研究区概况

研究区为云南省普洱市澜沧拉祜族自治县，该县位于普洱市西南部，介于99°29′～100°35′E，22°01′～23°16′N之间，总面积为8807 km2。澜沧县处于热带与亚热带交替区，因东临澜沧江而得名，其地形主要以山区为主，约占98.8%。海拔在600～1700 m；年平均气温17～22℃，光照时间长，全年光照时数长达2700 h；水资源年均总量为12270亿m³，水资源以降水补给为主要来源，年降雨量1052～3431 mm；土壤类型主要有砖红壤、赤红壤、红壤等。森林覆盖率约为74.7%，森林资源丰富，植物约为5600多种，主要树种类型为思茅松、水冬瓜（Alnus cremastogyne）、木荷（Schimasuperba）和各种栎类、西南桦（Betula alnoides）和竹类等[14]。

2 材料与方法

2.1 样木数据

本研究所使用的数据来源于2013年思茅松样地调查数据，在澜沧县内采用皆伐的方式共调查思茅松天然林标准木36株。每株样木分别记录样木位置（GPS坐标）、地形因子（海拔、坡度、坡向、坡 位 等）、林 分 因 子（树 龄、胸 径（DBH）、树高（H）、冠长、冠幅等）。

2.2 样木生物量测定及含碳量计算

单木生物量测定通过分器官分别测定得到，主干部分因无法全部称重因而采用材积密度法测定主干部分的生物量，伐倒木所测定的因子为长度、直径等因子，用于套算材积，分段称取鲜质量并取样；枝、叶使用分级标准枝法完成测定；主根生物量分段称重并取样，侧根全称质量并取样，同时记录主根根长及基径。将称重后的树干、树枝、树叶、树皮、树根及枯枝等样品烘干后用打磨机磨碎，用CN分析仪测量打磨后各样品的含碳率，再乘以对应的生物量计算单木各维度的含碳量[15-16]。单木各维度含碳量统计情况见表1。

表1 样木含碳量基本统计表Table 1 Basic statistics of carbon content of sample wood

2.3 空间效应分析

2.3.1 Moran'sI

定量描述空间自相关的方法主要有Moran'sI、Geary's C和Getis'G等[17-19]，本研究采用Moran'sI解释描述思茅松含碳量的空间相关性。全局Moran'sI定义公式如下：

式中：xi和xj为样点i 与j上的含碳量实际测定值，x¯为研究区域内所有含碳量实际测定值的平均值，n为样本数，wij(d)为根据样点i 与j之间距离计算的权重。

全局Moran's I数值在(-1, 1)之间，当Moran's I ＞0时，其相关性为正相关，Moran'sI数值＜0时，其相关性为负相关，其值为0时则不相关，表明空间数据随机分布。本研究采用投影坐标，使用Rookcase计算单木各器官的Moran'sI指数[20]。

2.3.2 半变异函数

半变异函数是描述空间中数据点的半变异值与点之间空间距离的函数，对半变异函数的图形描述可得到一个数据点与其相邻数据点的空间相关关系图，是描述区域化变量随机性和结构性特有的基本手段[21]。块金值表示随机部分引起的空间变异，如实验误差，取样间距较小；基台值反映了总的空间变异程度，即随机部分和固定部分空间变异；块金值与基台值的比值称之为基底效应，它反映了空间自相关的程度，该值越小空间自相关越高；变程（a）表示存在空间自相关的最大距离，当观测值之间的距离大于变程时，观测值的自相关程度消失。本研究采用ArcGIS地统计模块计算各器官含碳量的半变异函数。常用的半变异函数模型主要有圆形、球状、高斯、指数模型，经实验，球状模型效果较好。

式中：C0为块金值；C0+C为基台值；a为变程；C为偏基台值。

2.4 建模方法

2.4.1 最小二乘模型

最小二乘模型是一种全局模型，它以因变量的实测值y与估测值yi之间最小残差平方和（RSS）的拟合方法，模型内的常数及变量系数与地理位置无关，即不同的研究区间，变量系数不变，它体现不出各区域间的空间差异性。本研究采用GeoDa构建各器官含碳量的最小二乘模型，其定义公式如下：

式中：y为含碳量的估测值；p为 变量个数；βi为 变量拟合参数；ε为模型残差，ε ～N

2.4.2 地理加权回归模型

地理加权回归模型是一种局部模型，它是普通线性回归模型的扩展，将数据的地理位置嵌入到回归参数之中，它的公式定义如下：

式中：(ui,vi)为第i个采样点的坐标， βk(ui,vi)是第i个采样点上的第k个回归参数，是地理位置函数，

空间权函数和带宽的选取对地理加权回归模型的预估精度非常重要，常见的空间权函数类型有高斯函数（Gaussian）和截尾型函数法（bisquare）[22]。带宽是一个非负参数，用来表示权重与距离的关系，相对空间权函数而言，对地理加权回归模型的影响更大。带宽太大，会增加模型偏差；带宽太小，模型参数的估测方差会增大。常用的最优带宽的判定方法有最小信息指数（AIC）、贝叶斯信息准则（BIC）、交叉检验（CV）。本研究利用MGWR构建各器官含碳量的地理加权回归模型[23]，空间权函数选择固定距离bi-square，选择自动带宽，其评判方法为最小AIC信息指数。

2.5 模型评价

选用决定系数（R2）、Akaike信息指数（AIC）、均方根误差（RMSE）、模型有效性（ME）和模型残差系数（CRM）进行模型评价。R²越大，AIC、RMSE越小，模型拟合效果越好；ME越是接近于1，CRM越是接近于0，其模型估测精度越高，如果CRM的值大于0，表明估测值与实测值相比偏小；CRM的值如果小于0，则说明估测值与实测值相比偏大。

模型评价指标计算方法见式（5）～（9）。

式 中：yi为 实 测 值；为 估测 值；为 实 测 平均值；n为样本数；p为参数个数。

3 结果与分析

3.1 单木各器官含碳量的空间自相关

30～300 m范围内的单木树干、树皮、树枝、树叶、树根和全树含碳量的Moran'sI见图1。由图1可知，在滞后距离30 m范围内，树叶的全局Moran'sI为0.05，其余器官除树枝外的全局Moran'sI均大于0.6，表现出较强的正自相关；在250 m范围内，树干、树皮、树根和全树含碳量表现为正相关，树枝含碳量的空间自相关变化较大，在150 m范围内表现出正相关，随距离增加表现出负相关。树叶含碳量的Moran'sI始终随距离变化不显著，由于树叶部分含碳量相对其他器官含量较低引起的。综上所述，思茅松单木的树干、树皮、树枝、树根和全树含碳量均存在不同程度的空间自相关，且随距离增加，空间自相关性越小。

图1 思茅松各器官含碳量空间自相关关系Fig.1 Spatial autocorrelation of carbon content in each dimension of P.kesiya var.langbianensis

3.2 单木含碳量的空间异质性

30～300 m距离内各器官含碳量全向的块金值与基台值见图2。由图2可知，除树叶外，200 m内各器官含碳量均存在一定的空间异质性，其中树根含碳量的变异最大。说明树根的空间异质性受生长环境影响较大，如海拔、坡度和土壤成分。树叶的空间异质性不显著，是由于树叶含碳量相对其他器官来说碳含量较小，受地形和环境的影响不大。树干、树皮、树枝和全树含碳量在200 m内具有相同的变化，均是随着距离的增加，空间变异先增加后稳定，说明200 m范围内的单木具有相似的生长环境，其含碳量的变异情况差别不大，200 m是研究思茅松生态过程空间变异的拐点。

图2 单木各器官含碳量全向基底效应尺度变异Fig.2 The omnidirectional base effect scale variation of carbon content in individual wood organs

3.3 最小二乘模型构建

采用胸径、冠长、海拔、坡度和胸径平方乘以树高5个变量构建各器官含碳量的最小二乘模型，以R2和AIC评价模型，确定树干、树皮、树枝、树叶、树根和全树含碳量的模型形式。拟合情况见表2，除树叶决定系数为0.303，相对较差，其余维度模型拟合效果较好，树干、树皮、树根和全树含碳量的决定系数均在0.8以上。DBH2H除树叶和树枝含碳量的模型构建外，其余维度均未引入。

3.4 地理加权回归模型构建

基于最小二乘模型（OLS），构建各维度含碳量的地理加权回归模型（GWR），各维度回归参数见表3。对比表2～3，OLS的回归参数是固定的，地理加权回归模型会在每一个地理坐标上给出一个模型，所以每一个地理坐标都会有一套参数。表3为地理加权回归模型的局域系数的几种统计量，包括：最小值、第一四分位数值、中值、第三四分位数值及最大值。OLS的回归参数基本上接近GWR回归参数的中值。以思茅松全树含碳量OLS和GWR回归系数统计结果见表4，GWR变量的回归系数皆在OLS各变量回归系数正负1个标准差的区间，进一步说明模型数据间存在非平稳性，即思茅松含碳量存在空间异质性[24]。

表2 思茅松各维度含碳量OLS拟合结果Table 2 OLSfitting resultsof carbon content of P.kesiya var.langbianensis

表3 思茅松各维度含碳量GWR拟合结果Table 3 GWR fitting results of carbon content in each dimension of P.kesiya var.langbianensis

表4 全树含碳量OLS和GWR回归系数统计Table 4 Statistics of regression coefficients of OLSand GWR of the whole tree carbon content

续表 3

3.5 模型精度评价

各维度含碳量OLS和GWR精度评价结果见表5。由表5可知，整体上看GWR各器官的R²均在0.8以上，基本上都大于OLS，GWR模型的AIC信息指数和均方根误差都是小于OLS，GWR的ME值较OLS更接近于1，GWR的残差系数较OLS来说相对接近于0。以上表明GWR模型在本研究中优于OLS。

表5 各维度含碳量OLS和GWR精度评价比较Table 5 Comparison of accuracy evaluation of OLSand GWR of carbon content in various dimensions

续表 5

为进一步对比模型的预测效果，将各器官含碳量的OLS和GWR模型实测值和预测残差绘制成散点图，结果见图3。由图3可知，树干、树皮、树枝、树叶、树根和全树含碳量的GWR模型残差均在OLS模型残差范围内。综上所述，思茅松单木含碳量的GWR模型优于OLS，GWR解决了思茅松含碳量的空间异质性问题，显著提升模型估测精度，这与前人研究的结果一致[25]。

图3 各器官含碳量OLS和GWR残差对比Fig.3 Comparison of residual carbon content of OLSand GWR in various organs

4 结论与讨论

4.1 讨论

用不同滞后距离绘制的全局Moran'sI相关图表示思茅松各器官含碳量随距离的趋势变化。当距离为30 m时，思茅松各器官除树叶含碳量的Moran'sI均大于0.24，树干、树皮、树根和全树含碳量的Moran'sI均超过了0.61，说明思茅松各器官含碳量存在不同程度的空间自相关。总的来看，树叶含碳量的空间自相关不显著；其余器官含碳量的Moran'sI在250 m处趋于平稳，300 m处接近为0，且不再变化。说明300 m后，思茅松单木各维度含碳量空间自相关消失。

通过半变异函数对思茅松各维度含碳量空间异质性的定量分析和表达，确定思茅松不同尺度下的空间变异趋势，从而获得研究尺度的最优空间格局，也就是半变异函数的拐点。在本研究中，思茅松各器官含碳量均存在空间异质性，200 m是思茅松各器官含碳量空间变异的拐点；各器官含碳量在200 m之前的空间异质性是由自相关引起的，200 m之后的则是由随机因素造成的。

思茅松单木各器官含碳量均存在一定尺度内的空间自相关和空间异质性，这与思茅松各器官含碳量在生长过程中受其自身器官间相互作用和周围环境密切相关。其中树根的空间自相关性最大，树叶的最小，且随着距离的变化，其他各器官的空间自相关性均呈现逐渐降低至平稳的趋势，但树叶的含碳量随距离变化发的趋势并不明显；与此同时，从空间异质性的角度来看，所有器官的含碳量会随着距离的增加到拐点趋于平稳，即表现为随机分布，但是树叶则没有这样的变化。由此可知，树叶从空间相关性和空间异质性两方面来看，其空间变异并不大，这说明树冠的固碳量与环境因子并无关系，但是树叶本身的固碳量是最低的，而其他器官均表现出明显的空间效应，因此在构建思茅松含碳量模型时，采用全局模型，即在数据分析前就假定了数据间独立不相关，实际所得的结果是有偏的，对各器官含碳量进行估测必然导致一定偏差。

地理加权回归模型是一种局域模型，它基于不同地理位置间思茅松采样点与其周围的影响，采用权函数加权观测值，从而实现对每一个局域参数和模型估计。因此采用GWR能够获得不同地理位置上思茅松各器官含碳量与变量间的关系，降低了数据间空间关系的非平稳性，即消除了空间异质性，因此模型的估测精度和拟合效果得到了有效提升。本研究采用地理加权回归模型构建了各器官含碳量的空间回归模型，结果表明，思茅松含碳量空间分布的影响因素与测树因子和地形因子相关。胸径、冠长、海拔和坡度数据容易获取，相关研究表明胸径大的树木，它的含碳量相对高；海拔和坡度对生物多样性分布有着重要影响，为解释环境对空间效应的影响，引入海拔和坡度。为了突显空间模型的应用优势，本研究同时构建了最小二乘模型作为对比，思茅松树干、树皮、树根和全树含碳量的OLS模型决定系数在0.8以上，这是由于采样数据分布较为均匀，对于全局模型来说，数据越均匀，拟合效果越好。但由于空间异质性的存在，OLS模型的拟合效果和估测精度不及GWR高；GWR在估测森林碳储量时，带宽和权函数的选取对于模型估测精度的影响非常重要，因此GWR在森林碳储量的相关研究中，提升模型的精度如何选择相应的方法和验证指标，是今后探索的方向。

本研究表明，在思茅松单木各器官含碳量的研究中，空间效应不可忽视，Moran'sI和半变异函数能够有效的解释思茅松单木各器官含碳量的空间自相关和空间异质性。相对全局模型，GWR模型对于思茅松单木含碳量的研究具有较好的预测能力，因此在今后单木含碳量问题的相关研究上可以选择GWR进行估测研究。

4.2 结论

本研究采用Moran'sI描述了思茅松单木树干、皮、枝、叶、根和全树含碳量6个维度30～300 m间的空间自相关。从研究结果来看，思茅松单木各器官均存在空间自相关，说明单木各器官含碳量在一定距离内存在类似的变化，随着距离的增加，类似的变化减弱，空间自相关消失。

地统计半变异函数是描述空间异质性一种很好的方法。通过设置不同的步长大小计算不同尺度下树干、树皮、树枝、树叶、树根和全树含碳量的块金值和基台值。可以看出，随着步长的增加，块金值和基台值区域稳定，用块金值和基台值之比描述各维度含碳量的相对空间异质性。从研究结果来看，思茅松树干、树皮、树枝、树叶、树根和全树含碳量一定距离内存在空间异质性，这与思茅松的生长环境、结构性因素（如气温、降水、地形地貌、土壤成分）和非结构性因素（如自然、人为干扰）有关。

本研究采用GWR和GeoDa构建了思茅松单木树干、树皮、树枝、树叶、树根和全树含碳量的地理加权回归模型和最小二乘模型，研究结果表明，地理加权回归模型的预测精度较最小二乘模型要好，是由于OLS模型假设不同地理位置间的观测数据不变，没有考虑到空间自相关和空间异质性的影响，一定程度上存在局限性；GWR是一种局部模型，它不仅考虑了不同采样点的地理位置，同时也考虑了采样点受周围的影响。GWR可以得到不同空间位置上思茅松单木各器官含碳量与影响因子间的关系，解决了OLS模型无法解释的空间异质性，因而模型拟合效果、预估精度都得到了提升。本研究所用变量胸径、冠长、海拔和坡度获取方式比较简单，结合GWR模型可以为今后单木碳储量的估测提供更好的方法。