土壤是人类赖以生存的重要资源, 是植物所需营养物质的重要供给源^[1]。了解土壤特性不仅能避免土壤养分失衡, 还有助于保持土壤肥力水平从而提高作物品质和产量^[2]。然而近几十年来, 不科学的经营管理方式导致我国许多地区的林地土壤养分失衡, 土壤退化的趋势日益明显^[3-4]。而土壤养分失衡不仅会导致作物生长受阻产量下降, 还会造成环境污染^[5]。因此, 了解土壤养分现状及其影响因素能为科学土壤管理提供理论依据。

香榧(Torreya grandis cv. Merrillii)为榧树(Torreya grandis Fortune ex Lindley)中的人工栽植品种, 是我国特有的珍贵经济树种, 其果实具有较高的食用、药用和经济价值^[6]。长期食用可降低患心、脑血管疾病的风险, 同时还可以明目健脑^[7]。成年香榧林每年每亩产值超过1万元, 其栽培效益居经济林之首^[8]。由于其较高的经济价值, 使得现阶段香榧产业发展迅速。然而快速的发展也带来了弊端, 当前香榧林地土壤经营中存在着许多问题^[7]。储开江等^[9]发现导致香榧枝条黄化现象的原因不是病虫害, 而是不适应的生长环境所产生不良反应。戴文圣等^[6]也报道了大量的施肥, 导致香榧产量和质量显著下降。目前有关香榧的研究主要集中于繁殖栽培、生态习性以及病虫害防治等方面, 关于区域性土壤养分状况及影响因素的研究较少^[10-12]。

地统计学和克里格插值法在土壤养分空间异质性研究中得到了广泛的应用。许多学者将GIS和地统计学方法相结合, 系统揭示了特定区域土壤养分的空间分布特征, 该领域的研究已成为土壤学和农业生态学的热点之一^[13-15], 例如Chen等^[5]采用克里格法对亚热带平原粮食生产区土壤性状和肥力进行研究, 以此进行针对性指导施肥。但大部分研究集中于农作物和蔬菜土壤, 缺乏对经济林土壤的研究。

以浙江香榧主产区土壤为研究对象, 采用土壤综合肥力评价法计算研究区土壤养分水平, 利用地统计学结合GIS技术绘制研究区养分指标分布图, 探明土壤养分空间分布特征及其影响因素, 为香榧林地土壤可持续发展提供数据支持及对策。

1 研究区概况

本研究以浙江省香榧主产区为研究区域, 选取了位于浙江会稽山区的诸暨、嵊州、柯桥及东阳四市(区)。香榧性喜温暖潮湿、山谷纵横、溪流迂回交叉的生态环境, 主要分布在海拔200—800 m的丘陵山地^[8]。该研究区属于亚热带季风气候, 年平均气温14—17℃, 年平均降水量1100—1700 mm, 年平均日照时数1900 h, 无霜期207—240 d^[10]。成土母岩主要包括流纹岩、凝灰岩和流纹质凝灰岩等, 土壤类型以红壤、黄壤和黄红壤为主^[6]。近年来浙江省香榧种植面积约为1.33万hm², 年产量800—1500 t, 占全国总产量的94%^[16], 其中研究区所在诸暨、嵊州、柯桥、东阳等会稽山脉香榧种植面积占全省的65%, 产量占全省的80%^[17-18]。

2 研究方法 2.1 样品采集和实验分析

本研究综合考虑了样点分布的均匀性及代表性, 在香榧林地主产区以1.0 km×1.0 km网格布设土壤采样点, 并进行准确定位。结合香榧主产区实际种植和分布状况于2019年5月, 共采集121个土壤样本。在选定的样点上, 按“梅花”形状布点, 在半径为10 m的范围内, 采集5个子样点的表层(0—20 cm)土壤样品, 混合均匀组成一个土壤样品, 重量约为1 kg。同时记录香榧林地采样点的立地条件、农户经营管理和香榧产量等信息。

将土壤样品带回实验室并进行自然风干, 捡去石块和动、植物残体等, 研磨并通过2 mm的筛子。从中取出一部分, 用玛瑙研钵研磨过0.149 mm的筛子, 并密封在聚乙烯袋中, 贴好标签, 保存待分析。

土壤pH采用土水比为1∶2.5的悬浊液测定；土壤有机质选取重铬酸钾氧化法测定；土壤碱解氮、有效磷和速效钾分别采用碱解扩散法、Olsen法和乙酸铵浸提-火焰光度法测定^[19]。

2.2 空间自相关分析

空间自相关分析现今已被广泛的应到地理学的研究中^[20]。Moran′s I指数能系统的揭示研究变量的空间自相关性, 通常包括全局和局部Moran′s I^[21]。全局Moran′s I用来描述整个研究区域的空间自相关性, 并通过单一值反映区域空间变量的相似性, 其取值范围是-1到1, 大于0表示正相关, 反之则表示负相关, 等于0则表示不相关^[22]。而局部Moran′s I主要是用来计算每个特定位置的空间自相关程度, 可以识别出局部空间的聚类并进行异常值的分析, 弥补全局性分析的不足^[23]。全局和局部Moran′s I的公式如下：

全局Moran′s I：

(1)

局部Moran′s I：

(2)

式中, z是变量Z的平均值；Z_i和Z_j分别是变量在空间i和j处的数值(i≠j)；σ²是变量Z的方差；W_ij是采样点间的距离权重。当局部Moran′s I > 0时则表示研究区的目标采样点与其附近的观测值具有一定的相似性, 即空间集聚区, 主要包括高值集聚区和低值集聚区两部分；当局部Moran′s I < 0时则表示空间区域异常, 主要包括高-低值异常和低-高值异常。本研究用全局Moran′s I值来呈现研究区土壤养分的空间自相关水平, 并采用LISA (Local indicators of spatial association) 局部自相关分布图来揭示局部Moran′s I的自相关水平^[24]。

2.3 地统计分析

地统计学是一种基于区域化变量理论, 以半变异函数为研究工具的空间分析方法, 是一门研究空间变异与结构的自然现象的方法^[22]。

克里格插值法是通过区域化变量的原始数据和半方差函数的结构性, 对未采样点区域化变量值进行的无偏最佳估算的一种广义线性回归的方法^[5]。半方差函数(semi-variance)是地统计学的基础, 被广泛的应用于定量描述土壤变量空间结构的变异性^[24]。其公式如下：

(3)

式中：γ(h)为样点空间间隔距离为h时的半方差函数；Zx_i和Zx_i+h分别是变量Z(x)在x_i和x_i+h位置的实测值；Nh是样点空间间隔距离为h时的所有观测样点的成对数目。在实际应用中, 当半方差函数γ(h)随着空间间隔距离h的增加而增加时, 从非零值达到基本稳定的常数, 这个常数被称为基台值(C₀+C)；在h=0时的半方差函数值称为块金值(C₀), 这可能是因采样误差和小于采样尺度的随机因素引起的^[15]。指数模型、球状模型、高斯模型和线性模型是常用的变异函数理论模型^[22]。

2.4 土壤肥力评价方法

土壤肥力一词由来已久, 几个世纪以来一直被用来指土壤支持作物生长的能力。建立完整的土壤评价结构有以下三个步骤：(1)评价指标选择；(2)计算指标隶属度值和权重；(3)计算综合肥力。

2.4.1 评价指标选择

基于前人的研究成果, 本研究选取了能够最大限度代表土壤肥力质量的重要指标, 包括pH、有机质、碱解氮、有效磷和速效钾。

2.4.2 指标隶属度值和权重

不同肥力指标实测值的量纲各异, 为消除其影响需将各肥力指标标准化, 本研究采用隶属度函数方法对指标进行归一化处理, 将每个指标转化为0.1到1的无量纲值^[25]。因植物的效应曲线不同, 将隶属度函数分成S型和抛物线型两种, 并将曲线型函数转换成折线型函数。结合前人研究以及林地土壤的肥力特征^[26], 各指标在折线型函数中折线点的取值见表 1。

本研究中仅有pH一种指标属于抛物线型, 其余四种指标均属于S型函数, 其公式如下：

(4)

(5)

式中：W_i是指标隶属度；x是指标的测量值；x₁、x₂、x₃和x₄是指标的转折点。

土壤评价指标的权重通过因子分析法确定。各评价指标的公因子方差所占的比例为权重值^[26]。

2.4.3 土壤综合肥力评价指数

在对各指标进行评价后, 需将单因素评价结果转换为由各指标所构成的综合土壤肥力评价结果。本文采用加法合成的方法, 将各指标的隶属度值进行加权求和, 以计算土壤综合肥力评价指数^[5], 其公式如下：

(6)

式中：IFI为土壤综合肥力评价指数；W_i为第i项土壤肥力指标的隶属度值；N_i为第i项土壤肥力指标的权重值。

2.5 数据处理

本文采用Kolmogorov-Smirnov(K-S)检验数据集的正态性, K-S_p> 0.05表示该数据符合正态分布, 需对不符合正态的数据进行对数转换。采用Pearson相关分析法研究土壤各项养分指标间的相关性。

采用SPSS 20.0对数据进行描述性统计、正态分布检验和多元统计分析。利用Geoda软件来分析和识别空间聚类和异常值以及计算空间自相关系数；采用GS+ 9.0软件进行半方差分析；采用R3.3.3进行相关分析；所有土壤理化指标空间分布图均使用ArcGIS 10.2软件绘制。

3 结果与讨论 3.1 香榧林地土壤养分的描述性分析

如表 2所示, 浙江香榧主产区土壤pH值的范围为3.58—6.81, 平均值为4.91。戴文圣等^[27]研究表明pH值会影响香榧的生长发育, 过度酸黏的土壤会导致产量降低, 品质变差, pH值为5.2时最佳。本研究区土壤酸度的变幅虽较小, 但过度酸化不利于提升香榧的产量和品质, 应及时引起当地林农和有关部门的重视^[7]。研究区土壤有机质平均值为30.60 g/kg, 碱解氮、有效磷和速效钾的平均值分别为136.77、15.02、153.42 mg/kg, 参照《浙江林业土壤》土壤养分分级标准, 研究区土壤有机质、碱解氮、有效磷和速效钾含量偏高^[28]。根据实地调查显示, 大多数地区林农大量施用复合肥(N∶P₂O₅∶K₂O=15∶15∶15/17∶17∶17)和除草剂。大量施肥虽使土壤养分含量总体提高, 但过高的氮和磷会影响香榧产量和质量, 导致落叶落果^[28]。土壤钾含量过高对植物生长虽影响较小, 但会造成资源浪费和环境污染等一系列问题。在浙江省其他经济林种植园(山核桃、竹子和茶叶)中, 不科学的肥料配比也导致了土壤养分失衡^[29-31]。

变异系数(CV)值可以用来描述研究变量的变异度, 能更好的反映离散程度。根据王政权等^[32]的研究报道, 当变异系数 < 10%时, 为弱变异; 当10%≤变异系数≤30%时, 为中度变异; 变异系数>30%时, 为高度变异。如表 2所示, 仅pH属于中度变异, 其余4种养分元素均表现为高度变异, 表明香榧林地土壤肥力之间的差异较大。

3.2 空间聚类和空间异常值分析

土壤pH、有机质、碱解氮、有效磷和速效钾的全局Moran′s I值均大于0, 分别为0.14、0.24、0.20、0.29、0.23, 呈现了显著正的空间自相关性(P < 0.05), 说明在整个研究区域土壤pH、有机质、碱解氮、有效磷和速效钾的目标值和它们临近的采样点具有一定的相似性。

局部Moran′s I的空间分布特征表明了土壤pH的高值集聚区主要分布在诸暨市, 低值集聚区主要分布在柯桥区(图 2)。土壤有机质、碱解氮、有效磷和速效钾的空间自相关分布较为相似, 高值集聚区主要分布在柯桥区和嵊州市, 低值集聚区主要分布在诸暨市和东阳市。不同地区的施肥方式存在差异且没有规范统一的科学管理措施, 是产生空间格局异质性的主要原因^[33-35]。而大多数高-低空间异常值主要分布在低值集聚区的附近区域, 与之相反, 低-高空间异常值主要分布在高值集聚区的附近地区。

3.3 香榧林地土壤养分的空间变异结构特征

采用地统计学方法对浙江香榧主产区土壤养分进行空间结构和变异特征分析, 进行半方差函数拟合, 并根据最大R²选取最佳拟合模型^[15]。如图 3和表 3可知, 土壤pH符合高斯模型；土壤有机质和碱解氮符合指数模型；土壤有效磷和速效钾符合球状模型。

块基比又称为基台效应, 即块金值与基台值的比值, 可以衡量空间依赖的程度。Cambardella等^[36]将块基比比值分类为25%、25%—75%和75%, 分别表示强、中、弱的空间相关性。变程也是半方差函数中的一个重要的指标, 它主要反映了在一定区域尺度下空间自相关性的作用和影响范围^[20]。从表 3可知, 土壤pH、有机碳、碱解氮、有效磷和速效钾的块基比均在25%—75%之间, 表示中等变异程度的空间相关性。土壤养分的分布是自然因素和人为活动共同作用的结果。自然因素会增强土壤养分变量的空间相关性, 而人为活动将削弱其空间相关性并向均质化方向发展^{[34, 37]}。根据块基比(表 3)可知, 这些变量受到自然因素(地形、土壤类型、母质、气候等)和人为经营(施肥、除草等)的共同影响。五种变量的变程均较小, 其中pH的变程最小为3.29 km, 有机质、碱解氮、有效磷和速效钾的变程较为相似, 分别为8.52、5.84、8.82、6.49 km, 表明各变量与人类活动密切相关^[22]。

3.4 香榧林地土壤养分的空间分布格局

土壤养分状况需要借助土壤养分水平分级标准来进行评价。目前对香榧林地土壤的研究较少, 迄今还没有相应的土壤养分分级标准来衡量香榧土壤养分状况。本研究基于浙江林业常用的养分等级划分标准将土壤pH、有机质、碱解氮、有效磷和速效钾划分为四级^[28]。通过普通克里格插值法, 绘制香榧主产区土壤养分的空间分布图。如图 4所示, 土壤pH低值区主要分布在柯桥区和嵊州市, 高值区主要分布在诸暨市和东阳市。总体来说, 大部分地区的土壤pH值较低, 柯桥区和嵊州市的土壤酸化情况尤为严重。香榧林地土壤pH值会受到成土母岩以及人为施肥等因素的共同影响, 尤其是长期滥用化学肥料会导致土壤酸化日趋严重^[7]。马闪闪等^[31]报道了采用生石灰、土壤调理剂有利于改良临安山核桃土壤酸化。研究区土壤有机质、碱解氮、有效磷和速效钾的空间分布格局具有一定的相似性, 高值区主要分布在柯桥区和嵊州市, 少数分布在诸暨市东部, 低值区主要分布在诸暨市和东阳市。诸暨东部为早期的香榧种植区, 长期集约经营造成了土壤养分的富集^[38]。土壤养分受气候和地形等自然条件以及人为因素的共同作用^[34]。本研究中4个区域的自然气候和地理条件相似, 地形均以丘陵为主, 因此经营管理措施的不同是导致本研究区土壤养分差异的主要因素。研究区土壤养分空间分布特征与上述局部Moran′s I所揭示的空间分布特征一致。

3.5 香榧主产区土壤肥力指标相关性分析

相关性分析是揭示土壤养分元素和环境变量关系的有效方法。如图 5所示, 海拔仅与土壤有机质呈显著的正相关(P < 0.01), 而其他地形变量(坡度和坡向)与土壤养分元素的相关性较弱。这一结果表明, 地形变量不是影响土壤pH和养分的主要因素, 这与Dai等^[15]的研究结果一致。从图 5可知, 林龄与土壤pH呈显著负相关, 与有效磷呈极显著正相关(P < 0.01), 与有机质、碱解氮和速效钾呈显著正相关(P < 0.05)。表明林龄对土壤pH、有机质和碱解氮含量具有显著影响。

本研究区所有的土壤养分元素两两间均呈现出显著的相关性(P < 0.01)。其中土壤pH与有机质、碱解氮、有效磷和速效钾均呈现显著的负相关关系, 说明在一定程度上, 土壤酸化会影响土壤养分的供给, 过量施肥使林地土壤养分含量提高的同时也会导致土壤酸化^[8]。土壤有机质与pH、碱解氮、有效磷和速效钾均呈现显著的正相关关系, 其中与碱解氮的相关系数(0.747)最高, 表明土壤有机质含量与土壤的供氮能力密切相关。张建杰等^[39]进一步研究了土壤氮素和有机质的空间变异规律, 表明土壤氮素绝大部分来自于有机质。而土壤碱解氮、有效磷和速效钾的分析结果具有相似性, 与pH的相关系数相对较低, 与其余元素的相关系数较高且差异不大。这一结果可能是由于在人为经营条件下, 肥料是调控土壤氮磷钾三要素有效态水平的主要因素^[15]。

3.5.1 林龄的影响

不同林龄条件下各地区香榧林地土壤养分含量的差异性如图 6所示。根据各地区林龄实际情况, 将林龄分为 < 10年、10—30年、30—100年和>100年。由图 6可知, 土壤pH随着林龄的增长有降低的趋势, 且 < 10年和>100年的土壤pH差异显著。表明随着香榧林地种植年限的增长, 土壤呈现出酸化的趋势^[8]。对于土壤肥力质量而言, 土壤的进一步酸化不利于实现可持续发展的目标。随着林龄的增长土壤有机质、碱解氮和有效磷含量也在逐渐增加, 且幼林(< 10年)均与古榧林(>100年)差异显著。林地经营年限的增长, 立地环境基本稳定, 有机质的积累速率逐渐大于分解速率, 再加上有机肥的施用等经营措施, 使土壤有机质和碱解氮含量逐渐增长^[7]；磷素在土壤中移动性较弱, 所以随着林地经营年限的增长有效磷在土壤中得到了富集^[40]。土壤速效钾含量与其他养分含量变化规律不同, 没有随经营年限的增长得到明显的累积。钾素在土壤中的移动性较强, 易受水土流失、淋溶流失、径流流失的影响^[7]。不同年份气候状况的差异导致不同的径流和淋溶效应, 使各林龄下土壤速效钾含量具有差异。

3.6 香榧主产区土壤综合肥力评价

土壤综合肥力评价法有助于各地区农业管理者、相关研究人员和决策者更加深入的了解区域农业生产管理系统的土壤质量状况。本研究选取pH、有机质、碱解氮、有效磷和速效钾5个指标, 通过土壤综合肥力指数法(IFI)计算土壤肥力。采用隶属度函数并根据每个评价指标的阈值将变量转化为0.1—1.0之间的值(表 1)。通过因子分析法得出各评价指标的公因子方差值和权重值(表 4)。各变量被分配到的权重值差异不大, pH和有机质分配到的权重较高(0.215和0.209), 碱解氮、有效磷和速效钾的权重分别是0.185、0.194和0.197。

根据土壤综合肥力评价方法(IFI)计算研究区土壤肥力指数, 并通过克里格插值法绘制出土壤肥力空间分布图(图 7)。土壤肥力指数(IFI)>0.47属于中、高质量^[5], 本研究区有61%地区的土壤肥力属于中、高质量, 有39%地区的土壤肥力水平相对较低。如图 7所示, 高肥力土壤主要分布在柯桥区和嵊州市, 少数分布在诸暨市东部, 低肥力土壤主要分布在诸暨市和东阳市。表明研究区香榧林地土壤肥力水平总体较高, 大多数已达到肥沃水平, 少数地区土壤综合肥力指数较低。根据实地调查, 近年来各地区林农大量施用复合肥(N∶P₂O₅∶K₂O=15∶15∶15/17∶17∶17)以达到提高产量的目的。虽使香榧林地土壤养分得到了明显提升, 但一味的追求高产量, 而不注重成本和肥效不仅会使肥料浪费环境污染, 还会导致果实质量和产量下降^{[34, 41-42]}。土壤肥力水平不仅取决于土壤养分和作物吸收能力, 还受各因子协调程度的影响。因此, 应根据香榧林地土壤肥力的实际状况, 采用测土配方施肥制定施肥结构和用量, 以满足实际生产需要, 又不污染环境, 最终实现香榧产业的可持续发展。

4 结论

香榧林地土壤有机质、碱解氮、有效磷和速效钾低值区主要分布在诸暨市和东阳市, 高值区主要分布在柯桥区和嵊州市, 而pH的空间分布格局则与之相反。土壤综合肥力指数表明柯桥区和嵊州市土壤较为肥沃, 诸暨市和东阳市土壤综合肥力指数较低。香榧主产区土壤酸化以及养分失衡现象较为严重, 建议采用生石灰、土壤调理剂等进行酸化土壤改良, 在施肥过程中采用“稳氮降磷控钾”等方式, 以调节香榧林地土壤养分, 满足香榧不同生长阶段的养分需求。