建瓯市竹林分布与土壤有机质空间变异研究*

王 鑫 戴倩倩 徐梦洁 庄舜尧

(1.中国科学院南京土壤研究所，江苏 南京210008；2.南京农业大学 公共管理学院，江苏 南京210095)

我国是竹林资源最丰富的国家，竹林不仅具有巨大的经济效益和社会效益，而且其碳汇能力在应对气候变化中也扮演着举足轻重的作用[1-2]。土壤养分是竹林植被赖以发育的基础，土壤提供了植物生长必要的营养元素，是植被生长所需营养的重要来源[3-4]。对竹林土壤变异性的研究有助于改善竹林的综合效益，从而推动竹产业及上下游相关产业的发展[5]。目前许多学者对竹林土壤属性的空间分布特征以及竹林土壤养分的影响因素进行了相关研究，如高强伟等[6]对蜀南竹海景区毛竹Phyllostachys edulis 林土壤物理性质和空间异质性进行了分析，刘军等[7]在杭州市临安区和余杭区、苏木荣等[8]在佛山市以及沈德才等[9]在东莞市也开展了类似的研究。与此同时，随着遥感技术的快速发展，利用遥感技术提取森林信息成为获取林业基础数据的常用途径之一[10]，如许章华等[11]、官凤英等[12]利用遥感技术提取竹林分布信息，对福建省顺昌县竹林的空间分布特征和动态变化开展了研究，李阳光等[13]基于Landsat数据提取了浙江省的竹林信息并分析了其时空演变规律。因此为了准确、高效的获取竹林资源的变化信息，通过结合遥感技术和地理信息技术，丰富数据来源，及时更新林业信息，为林业资源动态监测和管理提供更加高效和多元化的管理措施，契合林业资源研究的前景[14]。

福建省建瓯市是全国重点林业县，2018 年森林覆盖率达79.6%，是竹林富集区之一，不少研究人员以此为研究区域开展了竹林林地土壤的研究。黄启堂等[15]分析了当地不同毛竹林林地土壤的理化性质，董蔚[16]总结了毛竹林土壤的相关属性的特点及其空间分布特征，林振清[17]着重研究了毛竹林土壤速效氮磷钾含量及其空间分布，Zhang等[18]关注土壤容重、砾石含量、阳离子交换量、土壤有机碳等属性的空间变异性，鲍思屹等[19]纳入时间维度，探讨了近40 a 来竹林土壤有机质的时空变异特征。这些研究通常在GIS 软件辅助下展开，使我们深入了解当地毛竹林土壤属性的空间分布及其伴随时间的动态变化，但是基于遥感影像的竹林土壤属性研究却很少，所使用的底图中竹林的分布信息往往相对滞后，因而影响了数据的即时性。与此同时，土壤有机质是土壤的重要组成部分，既能为植物生长提供各种营养成分，又在提升土壤肥力和维持土壤生态系统稳定方面发挥着重要作用[20-21]。因此本文以建瓯市为研究区域，基于遥感影像提取建瓯市竹林分布信息并形成研究底图，以乡(镇)为单位进行竹林土壤采样，对测定后的土壤数据进行分析，通过半方差分析和空间插值来揭示研究区竹林土壤有机质空间变异特征。研究能提供更新后的竹林分布专题图，研究结果也能为当地竹林产业的管理与规划提供参考。

1 研究区概况

建瓯市位于福建省北部，是福建省面积最大、闽北人口最多的县级市，素有“中国竹子之乡”的美誉。地处东经117°58′45″～118°57′11″，北纬26°38′54″～27°20′26″的东南沿海低山丘陵区，属亚热带海洋性季风气候，年平均气温19 ℃、降雨量1 600～1 800 mm，日照1 612 h，无霜期286 d，气候温和，雨量充沛，境内毛竹林面积、立竹量、竹材和鲜笋产量均居全国县(市)之首，竹产业已经成为当地的支柱产业[22]。

2 数据来源与分析方法

2.1 遥感数据

2.1.1 数据来源 本文选用的遥感影像来自于2017 年的Landsat-8 OLI 影像和同期30 m 分辨率的DEM 数据，下载自中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn/)。

2.1.2 数据处理 由于遥感影像的精度会受到如卫星姿态、卫星运行轨道和大气反射等因素的影响，因此在进行遥感解译前，需要对原始影像进行大气校正、几何精校正、图像裁剪和波段融合等预处理[23]，其中要用到的方法有临近点插值法[24]、直方图匹配法[25]、卷积算法[26]。

对于预处理后的遥感影像，需要筛选信息提取的最佳波段组合，本文使用最佳指数法[27](OIF)进行筛选，通过计算出影像数据各个波段间的标准差，构建相关关系矩阵。在相关关系矩阵的基础上，计算出三波段组合的OIF 值。OIF 值越大，说明相应组合获得的图像信息量就越大，组合效果趋于优化。

筛选出竹林提取的最佳波段组合后，再采用最大似然分类法(Maximum Likelihood Classification)[28]进行监督分类。最大似然分类法基于贝叶斯准则，假设遥感影像的每个波段数据都呈正态分布，通过训练样本计算出其平均值及方差、协方差等特征参数，从而得出总体的先验概率密度函数，然后将待分类像元的均值、协方差等特征参数带入判别函数计算其属于各类别的概率，最后将像元分到归属概率最大的类别中[29]。为了验证该方法的有效性，我们还采用了支持向量机法进行分类[30]。遥感数据预处理等程序利用ENVI 5.3 和ERADS 9.2 完成。

2.2 竹林土壤数据

2.2.1 采样与测定方法 2018 年4 月在研究区内以乡镇为单位，根据竹林面积大小共布设了209个竹林土壤采样点，研究人员在采样时利用手持GPS 确定采样预设点的位置，并导出其经纬度坐标，然后采用五点采样法在以预设点为中心的100 m×100 m 矩形4 个顶点和中心点处各采集一份0～20 cm 竹鞭分布层土样，混合均匀后带回实验室风干处理。土壤有机质采用低温外热重铬酸钾氧化-比色法测定[31]。

2.2.2 数据处理 (1)描述性统计分析 对研究区209 个样点的土壤属性数据进行异常值剔除，对剩余的205 个样点统计各属性的最大值、最小值、平均值、中位数和变异系数等指标，统计分析利用SPSS 22.0 完成。

(2)地统计学方法 对于服从近似正态分布的土壤属性数据进行半方差函数拟合，半方差函数包括球状 (Spherical) 、高斯 (Gaussian) 、指数 (Exponential) 和线性 (Linear) 等模型[32]，通常根据决定系数(R2)最大，残差(RSS)最小原则确定最佳拟合模型[33]，并获取块金值、基台值、变程等参数，以进一步描述区域化变量空间变异性的强弱、结构等；以此为基础，采用克里格方法对土壤属性进行空间插值，生成土壤属性的空间分布图[34]。我们还采用了一种确定性插值方法——反距离权重法（IDW），以插值点与样本点的距离为权重进行插值[35]。半方差函数及理论模型拟合利用地统计学软件GS+9.0，克里格插值在ArcGIS 10.2 下完成。

(3)精度检验 利用ArcGIS 10.2 软件从原始样点中随机抽取20%作为独立样本验证点，剩下的80%作为实验数据参与空间插值过程。通过计算平均误差(ME)、平均绝对误差(MAE)和均方根误差(RMSE) 来评价模型的预测精度，其中ME、MAE 和RMSE 越小，预测精度越高[36]。

表1 各波段组合OIF 指数Table 1 The OIF index table of band combination

3 结果与分析

3.1 竹林提取

在ERADS 9.2 软件支持下，选取同期Google Earth 影像(2 m)分辨率作为标准地图，结合实地野外调查，共选择了15 个地面点为控制点。采用临近点插值法对图像进行重采样操作，对重采样后影像进行几何精纠正操作，再对纠正后的单波段原始影像进行精度检验[37]。通过精度检验后，使用建瓯市行政边界依次裁剪三幅影像各波段文件，用直方图匹配方法对图像进行匀色处理。图像增强则采用卷积算法，以3×3 分析窗口进行，得到预处理后的影像以提取竹林分布信息。

B1-B7 波段中，将7 个波段进行三三组合，计算三波段组合的OIF 指数，结果如表1 所示，根据各组合的OIF 指数以及典型地物在各个波段的亮度均值，筛选出最佳波段组合方案为356、236 和367，结果如图1 所示。

建瓯市竹林分布情况较为复杂，竹林总面积近8.93×105hm2，遍布10 个镇、4 个乡和4 个街道。为了更高精度进行竹林的判读和提取，采用最大似然分类法。按照监督分类的一般步骤[38]，首先将研究区内地物划分为竹林、水体、非竹林和其他四类。通过直观观察和辅助验证点判定发现，在356波段组合处理下竹林为绿色明显区别于236 和367波段组合，因此以356 波段组合为竹林分类提取的波段组合，以356 波段组合图像作为分类基础图像，每个大类选择至少100 个以上训练区，合成分类地物色彩特征，构建地物分类模板。利用混淆矩阵对分类模板的精度进行检验，如表2 所示竹林解译的Kappa 系数为89.40%[39-40]，满足竹林信息提取要求。其Kappa 系数为88.78%，因此仍以最大似然分类法分类结果为基准，提取出竹林分布专题图(图2)。

图1 各波段组合Fig.1 Band combination map

表2 最大似然分类混淆矩阵Table2 The confusion matrix of maximum likelihood classification

图2 建瓯市竹林分布Fig.2 Bamboo forest distribution in Jian'ou city

3.2 竹林土壤有机质空间变异特征

3.2.1 样本描述性统计分析 由建瓯市毛竹林土壤属性描述性分析表(表3)可知，建瓯市竹林土壤有机质含量变幅为7.22～95.27 g/kg ，平均值为39.47 g/kg，最大值与最小值相差接近13 倍。变异系数(CV)可反映土壤特性的空间变异程度[41]。从表3 可以看出，有机质属于中等变异水平(CV 介于10%和100%之间)。

3.2.2 土壤有机质的趋势面特征 由于全局趋势的存在，会对半变异分析过程产生影响，因此在半变异分析过程中要对全局趋势进行剔除[42]。从建瓯市采样点土壤有机质数据的趋势分析(图3)可以看出，在东西方向上二阶趋势不明显，在南北方向呈“U”型变化。

3.2.3 土壤有机质的空间结构特征 数据的正态性检验结果表明，土壤有机质不符合正态分布，在对数据进行对数转换后，通过正态性检验，符合半方差分析的要求。半方差函数拟合的结果表明，土壤有机质的最优拟合模型是指数模型，其参数见表4。

块金值(C0)反映由实验误差和小于实验取样尺度引起的变异，如果块金值较大，则表明较小尺度上的某种过程不可忽视[43]。基台值(C0+C)表示系统属性的最大变异，块金值和基台值的比值(C0/C0+C)称为块基系数，它反映的是变量的空间相关程度[44]。块基系数小于25%说明变量的空间相关性强，土壤属性受到结构性因素(如气候、土壤母质、地形等)的影响较大，而受到随机因素(如施肥、灌溉、耕作制度等)的影响较小；系数在25%～75%之间说明变量具有中等的空间相关性，土壤属性受到随机因素和结构性因素的共同影响；系数大于75%，说明空间相关性很弱，土壤属性受到结构性因素的影响较大[45-46]。变程(A0)指区域化变量在空间上的最大相关距离[47]。由表4 可知研究区竹林土壤的有机质表现出中等的空间相关性，其空间相关性范围为24 990 m。

3.2.4 土壤有机质的空间分布特征 在ArcGIS 10.2 软件平台下，应用普通克里格法和反距离权重法进行插值，插值精度检验结果见表5。由表5可知，两种插值方法的平均误差和平均绝对误差虽然相近，但普通克里格法的均方根误差明显小于反距离权重法，表明普通克里格的插值结果更为精确。

以竹林分布图为掩膜对插值后的数据进行提取和重分类，如图4 展示了研究区内竹林土壤有机质的空间分布特征。插值后的竹林土壤中的有机质变幅为22.53～76.37 g/kg，较样本点数据变幅有所缩减。这是由于普通克里格插值的平滑效应以及将插值结果转换为栅格数据所致[48]。其中东北部的龙村乡和川石乡以及西南部的南雅镇竹林土壤有机质较高；中部的芝山街道、建安街道和瓯宁街道竹林土壤有机质含量较低。土壤有机质含量介于22.53～32.88 g/kg 的竹林面积最大，占比35.95%，土壤有机质含量介于60.97～76.37 g/kg 的竹林面积最小，仅占比7.11%。

表3 土壤有机质描述性分析Table 3 Descriptive statistics of SOM

图3 建瓯市竹林土壤有机质趋势分析Fig.3 Trend analysis of SOM in Jian'ou city

表4 半方差理论模型与参数Table 4 Semivariogram theoretical model andparameters

图4 建瓯市竹林土壤有机质空间分布Fig.4 Space distribution map of SOM in Jian'ou city

4 结论与讨论

4.1 本研究采用遥感影像解译的方式，以及时准确地更新竹林分布信息。遥感影像信息提取的方法众多，包括平行六面体分类法[49]、最小距离分类法[50]、马氏距离分类法[51]、神经网络分类法[52]、最大似然分类法[53]及支持向量机法[54]等。而研究人员[55]基本认同最大似然分类法的效果较好，在Kappa 系数和总体分类精度上往往优于其他分类方法的分类效果，因而应用较为广泛。如白宇兴[56]在提取西安市未央区土地利用分类信息时分别使用了最大似然分类法和支持向量机法进行了解译，前者的整体分类精度(97.80%)和Kappa 系 数(97.07%) 均 高 于 后 者(91.78%、88.81%)。本文在竹林信息提取时选用了最大似然分类法和支持向量机法进行解译，其中基于最大似然分类法的信息提取精度达到了89.4%，分类精度也高于支持向量机法，验证了这一观点。

就土壤属性而言，董蔚[16]测得的建瓯市毛竹林土壤有机质均值为39.8 g/kg，其空间分布以西南及东北区为高；林振清[17]的文献中土壤有机质含量均值为39.8 g/kg，呈从中部向西南及东北方向逐渐升高的趋势；Zhang 等[18]得到的有机质均值为40.8 g/kg，中部地区有机质含量最低。本研究测得的结果与他们较为相近，土壤有机质变幅为7.2～95.3 g/kg，均值为39.5 g/kg，其中西南部以及东北部的乡镇(街道)有机质较高，而中部的有机质含量较低。中间区域的海拔相对较低，通常海拔越高，有机质含量往往越高[29]。傅平[57]于2008 年的研究中指出建瓯市耕地土壤有机质含量平均值为28.5 g/kg，属中等偏上水平，其中最大值为85.9 g/kg，最小值为2 g/kg。本研究结果要明显高于该值，这与竹林土壤的特性有关。竹林土壤之上覆盖有森林凋落物，在降水、森林土壤动物和微生物等因素的作用下，转移到土壤中，因此林地土壤有机质含量普遍要高于农田。本研究测定的样本有机质及其变幅与已有文献基本一致，体现了数据的合理性，在此基础上通过插值得到的空间分布特征也呈现相同的特点，也体现了研究结果的合理性。

4.2 本文以建瓯市为研究区域，基于遥感影像提取建瓯市竹林分布信息并形成研究底图，在此基础上以乡(镇)为单位，进行竹林土壤采样，通过半方差分析和空间插值来揭示研究区竹林土壤有机质的空间变异特征。最佳指数法的计算结果表明，356 波段为竹林提取的最佳组合方式。在此基础上采用最大似然法得到的竹林解译信息精度较高。利用遥感影像来提取竹林信息，一方面可以获取及时更新的竹林分布专题图，另一方面也令后续的空间插值结果更为客观、准确。插值分析表明，研究区竹林土壤有机质在不同方向的空间分布存在差异。由于土壤有机质受到结构性和随机因素的共同影响，即自然条件和人为因素的综合影响，因此建瓯市竹林的管理和规划应关注人为因素的影响[58]。建议在今后的研究中纳入其他土壤属性并探讨其空间分布和变异特征，同时注重遥感数据的应用与更新，从而更好的推动当地竹林产业的发展。