Spatial distribution for leaf area index and canopy closure of Phyllostachys edulis stand using Cokriging

研究区为浙江省安吉县（图 1）。安吉县地处浙江省西北部，30°23′~30°53′N，119°14′~119°53′E，东邻湖州市吴兴区、德清县；南接杭州市余杭区、临安市；西与安徽省宁国市、广德县交界；北连湖州市长兴县。安吉县气候宜人，属亚热带海洋性季风气候，多年平均气温为15.5 ℃。安吉县多年平均降水量为1 400 mm，属亚热带东部常绿阔叶林亚区、中亚热带常绿阔叶林北部亚地带。安吉县竹类资源十分丰富，其中毛竹林为5.527×10⁴ hm²，占有林地面积的45％左右。

Figure 1.  Research area and field plots

外业样地调查于2008年8下旬到9月初完成，共调查了55个30 m×30 m不同经营水平的毛竹纯林样地，其中集约经营11个，中等经营水平21个，粗放经营23个。调查内容包括样地经纬度、海拔、坡度、坡向、胸径、年龄、样地内毛竹株数、郁闭度等，其中，样地株数最少的为46号样地（153株），最大的为15号地（16 033株）。另外，样地毛竹林年龄结构主要以1~3度为主，4~5度也有少量分布。研究区毛竹林分布及调查样地如图 1所示^[25]，可见样地基本涵盖了安吉县毛竹林分布区域。在森林资源调查中，目测法是郁闭度最常用且简便迅速的方法，2003年国家林业局《森林资源规划设计调查主要技术规定》中指出，有林地小班森林郁闭度可以通过目测法获取^[26]，因此，本研究样地郁闭度是在安吉县林业局技术人员指导下采用目测法获取的。样地毛竹叶面积指数采用式（1）计算^[27]。

式（1）中：$ \bar{D}$表示样地内毛竹平均胸径，N表示样地内毛竹总株数。

克里格法是利用区域化变量即原始数据和变异函数的结构性，对未采样点的区域化变量取值进行线性无偏最优估计的一种方法^{[8, 15]}，因此，变异函数及其理论模型是克里格估算的基础。变异函数定义为间距为h的2个区域化变量及差的方差的一半，如式（2）所示：

式（2）中：N（h）为距离相隔为矢量h的所有点对的个数。变异函数理论模型一般包括球状模型、指数模型、高斯模型、线性模型等，对半方差函数进行理论模拟，并在模型优选的基础上通过模型参数如块金值、基台值、变程以及结构比分析区域化变量的空间分布特征^[28]。通常选择决定系数（R²）大，残差小的理论模型作为变异函数拟合模型^[29]。

克里格包括简单克里格、普通克里格、趋势克里格、协同克里格、因子克里格、块状克里格等，其中协同克里格法（Cokriging）是用1个或多个次要变量（协变量）对感兴趣变量（主变量）进行插值估算。协变量与主变量都有相关关系，并且假设变量之间的相关关系能用于提高主要预测值的精度，估算公式如下：

式（3）中：Z_cok（u₀）为待估点u₀处的估测值；Z_a（u_i）和Z_b（u_i）分别是主变量和协变量在u_i和u_j处的实测值；λ_a和λ_b分别为主变量和协变量的权重且$ \sum{{{\lambda }_{{{a}_{i}}}}=1}, \sum{{{\lambda }_{{{b}_{j}}}}=0};$；n和m分别为参与估测的主变量和协变量样本个数。

首先，设置主变量和协变量。研究表明：叶面积指数和林分密度存在指数关系，与冠层郁闭度也有密切的相关性^[30]，因此，根据协同克里格算法，在对叶面积指数进行空间估算时，叶面积指数作为主变量，冠层郁闭度作为协变量，而对冠层郁闭度进行空间估算时，冠层郁闭度作为主变量，叶面积指数作为协变量。其次，主变量和协变量设置好之后，利用ArcGIS地统计模块，对55个样地数据统计分析、变异函数计算、理论变异函数的最优拟合及检验，并在此基础上进行克里格插值，实现整个研究区叶面积指数和冠层郁闭度估算；最后，利用研究区毛竹林分布图（图 1），对估算结果进行裁剪，从而得到研究区毛竹林叶面积指数和冠层郁闭度空间分布。

克里格插值考虑了空间距离对估算结果的影响即离采样点越近，克里格估计误差就越小^[9]，但没有考虑地形条件。山区地形复杂，在野外布设样地时难免会遇到地形起伏，所以选择相对平坦的样地对结果进行检验是否合理^[27]。本研究从样地中筛选出30个相对平坦的样地，对叶面积指数和冠层郁闭度估算结果进行验证，验证的指标包括实测值与估计值之间的相关性R²，均方根误差（E_RMS，root mean square error），平均标准误差（E_AS，average mean error），标准均方根误差（S_RMS，root-mean-square standardized）^[29]，均方根值（S_RM，root mean square），其中均方根值是样本平方平均值再开方，平均标准误差是指标准误差的平均值，均方根误差均方根误差计算公式如下：

式（4）中：X_{obs, i}和X_{model, i}分别为第i个样本的真实值和预测值，n为样本个数。

另外，也可以通过均方根误差减少的百分数（R_RMSE）和平均标准误差减少的百分数（R_ASE）来表示预测精度的提高程度^[31]。R_RMSE和R_ASE的计算公式（5）（6）如下：

式（5）和式（6）中：E_RMSk和E_ASk为普通克里格预测的均方根误差和平均标准误差；E_RMScok和E_ASk为协同克里格预测的均方根误差和平均标准误差。

变异函数结构分析一般要求区域化变量符合正态分布，否则会产生比例效应^[32]。通过对毛竹林叶面积指数和冠层郁闭度统计参数（表 1）及直方图（图 2~3）分析表明：样地毛竹林叶面积指数和冠层郁闭度的平均值为与中值接近，存在轻微负偏，变异系数较小，分别为13.66%和14.08%。

Figure 2.  Histogram of canopy closure

Figure 3.  Histogram of leaf area index

P—P概率图（P—P probability plots）可以检验样数据的正态性，即根据变量分布累积比和正态分布累积比生成的图形进行正态性判断，如果数据是正态分布，被检验数据成一条直线^[33-34]。图 4a和4b分别是叶面积指数和冠层郁闭度的正态P—P图，由图可见叶面积指数和冠层郁闭度的正态P—P图基本成一条直线，且检验表明，在0.05置信区间内，叶面积指数和冠层郁闭度均服从正态分布。

Figure 4.  Normal P—P plot of leaf area index (a) and canopy closure (b)

根据决定系数最大和残差值最小的原则，得到研究区叶面积指数和冠层郁闭度最优模型理论变异函数均为球状模型，模型相关参数如表 2。

2中C₀为块金效应，表示区域化变量的随机方差；C为区域化变量的结构方差；C₀＋C为变异函数基台，为区域化变量总方差；C/（C₀+C）为结构比，是区域化变量空间相关性程度的指标，大于75.0%，说明具有强烈的空间相关性，在25.0%~75.0%，为中等的空间相关性，小于25.0%，说明系统空间相关性很弱^[35]；变程是区域化变量影响范围的大小，在变程以内，区域化变量是空间自相关的。

由表 2可见：叶面积指数变异函数的结构比为83.6%，冠层郁闭度变异函数的结构比为78.4%，均大于77.0％，说明两者变程范围内存在具有强烈的空间自相关性；另外，叶面积指数和冠层郁闭度间自相关范围也远远大于最小抽样间距（360 m），说明本次采样密度是合理的^[11]，因此，在变异函数模型的基础上可以对叶面积指数和冠层郁闭度进行空间估算。

协同克里格估算结果与毛竹林分布信息叠加裁剪得到的研究区毛竹林叶面积指数空间分布图，如图 5a所示。由于安吉县北部几个乡镇很少有毛竹林分布，没有布设样地，因此这些区域没有预测值。根据图 5a可看出：研究区内叶面积指数分布整体上呈西南到东北逐渐递减空间分布格局，西南地区几个乡镇叶面积指数值最高，最大值可达10.731，北部几个乡镇的叶面积指数值较低。叶面积指数协同克里格估算值与实测值之间的相关关系如图 5b所示。由图 5可见：两者具有较好的线性关系，R²为0.635 1，均方根误差较小为0.541 6（表 3），说明基于协同克里格插值得到的叶面积指数结果较好。

Figure 5.  Spatial distribution (a) and correlation between the predicted and measured (b) of leaf area index based on Cokriging

图 6a是研究区毛竹林冠层郁闭度空间分布估算结果。图 6a表明：研究区内冠层郁闭度分布整体上呈西南到东北逐渐递减趋势，中南部地区几个乡镇冠层郁闭度值最高，最大值可达0.870，北部几个乡镇的冠层郁闭度值较低。这一空间分布格局与叶面积指数具有相似的特点，也印证了先前关于两者具有密切关系的描述及可以相互作为协变量进行估算的合理性。图 7b是冠层郁闭度协同克里格估算值与实测值之间的相关关系。由图可见，两者之间的决定系数R²为0.428 5，相关系数为0.654 6，尽管相关性低于叶面积指数，但均方根误差也较小（表 4），说明估计结果在一定程度也是可靠的。

Figure 6.  Spatial distribution (a) and correlation between the predicted and measured (b) of canopy closure based on Cokriging

Figure 7.  Relationships between predicted and measured of leaf area index and canopy clisure based on ordinary Kriging

图 7是普通克里格插值得到的叶面积指数和冠层郁闭度预测值与实测值之间的相关关系，从图 5~6对比可以看出，基于协同克里格插值方法得到的研究区的叶面积指数和冠层郁闭度预测值和实测值的相关性明显高于普通克里格方法。表 3~4详细比较了2种方法估算的精度。通过对比发现：协同克里格估测叶面积指数和冠层郁闭度的均方根误差、均方根、平均标准误差、标准均方根误差等误差标均小于普通克里格法，其中叶面积指数均方根误差减少2.00%，平均标准误差减少0.18%，且预测值和实测值的相关系数提高了1.94%，而冠层郁闭度均方根误差减少1.90%，平均标准误差减少1.30%，预测值和实测值的相关系数提高了4.82%。以上分析说明，相对于普通克里格法而言，协同克里格估算得到叶面积指数和冠层郁闭度估计算精度高，更能反映：安吉县毛竹林叶面积指数和冠层郁闭度实际空间分布情况。

在地面调查的基础上，利用协同克里格插值法对研究毛竹林叶面积指数、冠层郁闭度进行空间分布估算研究，并与普通克里格插值法进行比较。研究表明：基于协同克里格插值得到的毛竹林叶面积指数和冠层郁闭度预测值与实测值之间的决定系数R²分别为0.635 1和0.428 5，说明叶面积指数、冠层郁闭度空间估算结果在一定程度上是可靠的，且协同克里格法预测精度高于普通克里格法。当然，受条件所限，本研究叶面积指数和冠层郁闭度预测与实测之间的相关性还有待进一步提高，将遥感技术和地统计分析方法结合进行综合估算，可能会改善两者预测精度。另一方面，地统计空间估算除协同克里格法和普通克里格法外，还包括趋势克里格、漂移克里格、因子克里格、块状克里格等多种方法，需要根据实际情况进行选择或分析比较，当区域化变量分布存在某种趋势时，可以考虑采用漂移克里格法，而当变量实测样本数量较少时，可以考虑采用协同克里格法^[8]。因此，本研究根据样本的实际情况选择协同克里格法进行叶面积指数、冠层郁闭度空间分布估算。另外，安吉县毛竹林面积占有林地面积的45%左右，对毛竹这个单一物种，本研究样地基本涵盖了毛竹林主要分布区域，且样本叶面积指数和冠层郁闭度呈正态分布的统计结果（图 4）也在一定程度上反映了样本的代表性和合理性，为空间插值提供了基础。地统计学是一种空间统计学，它在对样本进行统计分析时，与传统统计学最大的区别就是考虑了区域化变量的空间分布特征值，因此，可以在变异函数模拟的基础上进行空间分布估算及其变异分析，并将“点上”信息向“面上”信息尺度转换^{[8, 15]}。研究表明：球状变异函数模型能够揭示安吉县毛竹林叶面积指数、冠层郁闭度空间变异特征，并在此基础上通过协同克里格法实现了叶面积指数和冠层郁闭度样地信息（图 1）转换为空间分布信息（图 7a和图 8a）。分析表明，叶面积指数和冠层郁闭度在空间上具有强烈的空间自相关性，且其空间分布基本呈现出从西南到东北逐渐递减的格局。对比安吉县森林资源规划设计调查成果报告^[36-37]和相关文献[11, 25]发现：这一结果在很大程度上符合安吉县毛竹林立竹度高、生物量大等特征，也反映了安吉县西南部为毛竹产业示范区域、毛竹林经营集约化程度高的实际情况，说明采用协同克里格法和地面调查数据估算安吉县毛竹林叶面积指数和冠层郁闭度空间分布在一定程度上是可行的。

区域化变量的空间变异特征是土壤、气候、施肥、耕作方式、种植制度、管理措施、经营水平等众多内外因素影响的综合反映。浙江省安吉县为著名的“中国竹乡”，竹林资源经营管理和产业发展在该县林业发展中占据重要的位置。本研究关于毛竹林叶面积指数和郁闭度空间分布特征的研究结果，可为安吉县毛竹林经营管理提供参考价值。