浙江清凉峰自然保护区土壤肥力指标空间变异及其影响因素

姜霓雯,童根平,叶正钱,程樟峰,吕永强,傅伟军,*

1 浙江农林大学,亚热带森林培育国家重点实验室培育基地, 杭州 311300

2 浙江清凉峰国家级自然保护区管理局, 临安 311321

自然保护区是生物多样性保护的基石[1],在确保各类自然生态系统安全稳定、筑牢生态安全屏障、改善生态环境质量等方面发挥了重要的作用[2]。自然保护区不仅能够最大限度长期维护生物多样性,并且还为社会提供了一系列生态系统服务效益,例如提供清洁水、木材等[3]。其中,森林类型的自然保护区以其复杂的生态系统结构,在陆地生态系统中占据领先位置[4]。土壤层是植被生长的基础,森林生态功能的发挥和森林植物资源的保育都依赖于土壤资源,土壤资源的探明和保护也有利于区域内土壤动物以及微生物资源的保护[5—6]。土壤氮、磷、钾等肥力指标的空间异质性及其影响因素是近年来众多学者研究的热点之一,且多集中于对农田、经济林以及小流域的研究。如邹润彦等[7]利用地理加权回归模型,系统分析了环鄱阳湖区农田土壤有机碳影响因子的空间分布格局,并进行了制图研究；董佳琦等[8]探讨了香榧主产区林地土壤养分的分布规律及其主控因素；Wang等[9]揭示了太岳山亚高山森林流域土壤有机碳和全氮的空间变异规律。但是以往研究缺少对于南方丘陵区人为扰动较小的自然保护区内天然林地土壤肥力特征的报道,深入研究自然保护区内土壤肥力指标的空间变异情况能够精确了解天然林分特别是珍贵树种的生长条件,对森林生态系统的稳定和珍惜树种的保护起着重要作用。

清凉峰自然保护区地处华东皖浙丘陵区,属森林和野生动物类型保护区,1998年扩区并经国务院批准晋升为国家级自然保护区,不仅是我国一级重点保护野生动物华南梅花鹿的重要分布区域,也是我国经济发达的长三角地区保存完好的基因库。目前对清凉峰国家级自然保护区的研究主要侧重于森林植被和动物保护方面[10—11],对土壤养分等肥力特征缺乏全面、系统的调查研究。本研究旨在分析该区域土壤肥力指标的空间变异规律,并定量揭示地形、林分等因子对其空间变异的影响,以期为区域生态保护和森林管理提供科学参考。

1 研究区概况

清凉峰自然保护区位于浙江省杭州市临安区西北部(118°50′57″—119°13′23″E,30°00′42″—30°19′33″N),地处华东皖浙丘陵区,由千顷塘、龙塘山和顺溪坞3块区域组成,总面积为11252 hm2,主峰清凉峰,海拔1787.4 m,系浙西第一高峰,区域内地势高差悬殊,地形复杂多样。保护区位于中亚热带季风气候区北缘,年均温12.5℃,年降水量1862.2 —2331.9 mm之间。该区不同时代地层发育较为齐全,其中以侏罗系分布最广,主要包括中酸性火山岩、紫红色层状粉砂岩、粉晶灰岩、硅质页岩以及比较大面积的岩浆岩等。研究区内植物区系组成丰富,是国家极为珍贵的生物多样性宝库,主要珍稀特色植物种群包括华榛(Coryluschinensis)、台湾水青冈(Fagushayatae)、银缕梅(Parrotiasubaequalis)、香果树(Emmenopteryshenryi)、鹅掌楸(Liriodendronchinense)和巴山榧树(Torreyafargesii)等。本文研究区为千顷塘(图1),区域面积5690 hm2,是野生华南梅花鹿(Cervuspseudaxis)主要栖息地,在整个保护区地位举足轻重。研究区根据华南梅花鹿的分布程度分为核心区、缓冲区和实验区3部分,区域内土壤肥力指标的研究与梅花鹿的主要食物来源以及珍稀植物的生长情况息息相关。

图1 研究区样点分布及海拔信息

2 研究方法

2.1 样品采集与分析

研究区表层土壤样(0—20 cm)以地形图和植被类型分布特征为辅助信息,结合样点分布的均质性和科学性原则,采用1 km×1 km网格法进行室内布点,于2019年11月—2020年4月共实地采集土壤样品71个。将土壤样品及时送回实验室,风干、去杂、研磨过筛备用。采样时用GPS记录样点的实地坐标,以及海拔、植被种类等信息。

样品肥力指标测定:土壤有机质(SOM)用K2CrO7氧化还原滴定法,全磷(TP)用NaOH碱熔-钼锑抗比色法；全氮(TN)采用半微量凯氏法；全钾(TK)用NaOH熔融-火焰分光光度法；pH用pH酸度计测定；容重(bulk density)采用环刀法测定。

2.2 地形数据提取

研究区12.5 m×12.5 m数字高程模型(DEM)源于NASA官方网站(https://search.asf.alaska.edu/#/)；野外采样点实地坐标、高程数据在采样同时用GPS测定。参照国内外相关研究成果,本文利用ArcGIS 10.2空间分析工具直接提取基本地形数据,包括坡度(slope)、坡向(aspect)、曲率(curvature)、平面曲率(plancurvature,Ch)、剖面曲率(profilecurvature,Cv),其余衍生地形属性包括地形起伏度(amplitude of landforms,QFD)、地表粗糙度(roughness)、地形湿度指数(topographic wetness index,TWI),均结合空间分析工具、水文分析模块以及栅格计算器经复合计算获取。

2.3 地统计分析

半方差函数是地统计学中描述变量的定量工具,在土壤元素空间变异性研究中应用广泛,用于揭示区域化变量的随机性和结构性特征[12]。计算公式如下:

(1)

式中,γ(h)是间距为h时的半方差；Z(xi)和Z(xi+h) 分别是对应变量在xi和xi+h的实测值；N(h)是间距为h时的观测点对总数。

半方差分析的参数中,块金值(C0)表示由采样和检测分析误差引起的随机变异,基台值(C0+C)代表总空间变异程度,块基比(C0/(C0+C)表示随机变异的占比情况,用于衡量空间相关程度[13],变程(a)表示空间自相关的作用范围。

2.4 空间自相关计算

土壤属性在地理空间上与相邻区域的观测值存在的相互依赖性,即空间自相关,可用全局Moran′sI指数反应指标的空间自相关性大小,计算公式如下:

(2)

Moran′sI指数取值范围是-1—1,小于0表示负相关,等于0表示不相关,大于0表示正相关。Moran′sI的显著性水平采用下式检验:

(3)

式中,Z为检验Moran′sI指数的显著性的统计量；I为Moran′sI指数；E(I)为期望值；Var(I)为方差[14]。|Z|≥1.96、2.58分别是空间自相关显著和极显著的分界值。

2.5 数据处理

对肥力指标进行统计分析前,采用三倍标准差法(阈值法)剔除异常值[15],识别到的异常值使用新样本的最大值和最小值替代,保持原有样本量不变。应用SPSS 20.0软件进行土壤肥力指标的描述性统计、K-S法检验数据正态性,方差分析用于指标间的差异显著性检验、回归分析用以定量表达环境因素对指标空间变异的影响程度；半变异函数分析和相关模型参数的计算在GS+7.0中完成；利用Geoda进行空间自相关计算；地形因子的计算和插值图的绘制均借助ArcGIS 10.7软件完成。

3 结果与讨论

3.1 土壤肥力指标描述性统计分析

千顷塘土壤肥力指标元素含量描述性统计结果如表1所示,pH值为4.05—6.39,整体呈酸性；各指标含量均值对比全国第二次土壤普查属性分级标准[16],有机质、全氮、全钾均值含量分别为72.44 g/kg、1.91 g/kg、23.16 g/kg,总体含量较高,其中有机质含量处于极丰富水平(一级),全氮、全钾处于丰富水平(二级),全磷含量较低,处于中等水平(四级)。土壤各指标的变异系数介于10.38—48.54%之间,根据王政权等[17]划分标准,所有肥力指标均属于中等程度变异,其中pH变异程度较低。正态性检验结果显示pH、全氮符合正态分布,有机质、全磷、全钾经对数转换后符合正态分布,满足地统计学分析的要求。

表1 千顷塘土壤肥力指标描述性统计特征/(g/kg)

3.2 土壤肥力指标地统计分析

从图2和表2可以看出,pH、全磷、全氮、全钾的半方差理论模型均为指数模型,有机质为球状模型；R2分别为0.72,0.84,0.89,0.91和0.83,能较准确地反映千顷塘土壤肥力指标的空间变异结构特征；变程分别为0.03,0.99,0.19,0.02 km和7.17 km,说明不同元素的空间自相关尺度范围存在差异。块基比是用于研究区域化变量空间变异程度的重要参数[18],根据Gao等[19]划分标准,pH、全磷和全钾的块基比均处于0—25%之间,表现为强烈空间自相关,说明结构性因素起主要影响,例如地形因子、成土母质等。自然状态下,土壤中磷、钾两种元素主要来自于土壤母质的风化和枯落物的分解,其分布与土壤的淋溶特征有关,而地形因子较大程度上影响着土壤的淋溶作用[20]；有机质、全氮的块基比处于25—75%,均属于中等空间自相关,表明耕作活动、种植制度等随机因素存在一定的影响。研究区内植被垂直带分异明显,显著影响了土壤属性的异质性,也可能与复杂地势下多样性小生境的发育有关。

表2 土壤肥力指标半变异模型参数及全局Moran′s I指数

图2 千顷塘土壤肥力指标含量半变异函数分析

半方差函数的变程反映了有机质的空间变异尺度范围大于pH、全磷、全氮、全钾,表明环境因素在较大程度上控制着有机质的空间异质性,而其他指标元素的空间分布更加趋向于破碎化；通过块基比揭示pH、全磷、全钾的空间自相关性高于有机质和全氮,其参数分析结果可以为Kriging插值提供依据[21],进而对各类指标空间格局进行更准确的描述,但是无法对空间相关的显著性进行统计学检验[22]。全局Moran′sI指数能够较为直观的反映土壤属性在整个研究区域的聚集状态,采用随机条件下近似正态分布假设的标准差对数据进行标准化,置信区间双侧检验阈值为界限,得出有机质和全磷存在显著空间自相关(P<0.05,Z>1.96), 全氮具有极显著空间自相关性(P<0.01,Z>2.56),表明这3种指标呈聚集分布；而pH和全钾的空间自相关性不显著,表明其分布离散。从变程结果可以看出,pH和全钾变程很小,说明这两种指标在研究区内主控因素的空间连续性尺度较小,分布趋向于随机化,这与Moran′sI显著性检验结果较为一致。

3.3 土壤肥力指标空间分布格局

根据半方差分析结果参数,通过ArcGIS 10.7 软件进行普通Kriging插值绘制千顷塘各肥力指标含量的空间分布图(图3),综合来看,5种属性在空间上分布不均,斑块特征显著,表现出较为明显的空间异质性。对照土壤属性分级标准[16],对各级别指标含量比重进行统计(表3),千顷塘土壤肥力整体处于丰富水平,pH处于4.5—5.5的范围最大,占比达到70.42%,弱酸性pH 5.5—6.5的范围次之,研究区整体生态保存完好,土壤酸性主要来自于地层中的酸性土壤母岩,以火山岩和岩浆岩为主,具体包括流纹质斑岩、花岗闪长斑岩以及石英闪长岩等；有机质含量很高,含量>40 g/kg的点位占比87.32%,高值区分布在中偏西北部,整体呈现明显的“条带状”分布；全磷总体处于中等偏下水平,含量由西北向东南方向降低；全氮和全钾总体含量较为丰富,全氮的分布格局和全磷较为相似,全钾的高值区分布在西北部和东南部边缘地区。总体来看,研究区各类土壤肥力指标含量均较高,只有磷素含量偏低,可能与土壤酸性较强有关,相关研究显示,土壤pH<5.5的酸性环境下,铁铝等氧化物对磷素具有很强的固存能力[23]。各类土壤指标分布规律存在一定的相似性,具体体现在多种肥力指标在空间上呈现“两头高,中间低”的分布特点,高低值分布较为集中,斑块较大,分布较为连贯,表明该区域空间变异较小,空间相关性高。可见,千顷塘土壤肥力条件较好,对于森林生态系统的供给水平较高,研究区处于浙西海拔最高区域,日均气温较低,有机质矿化速度较慢,人为扰动小,不利于有机质分解,相比于其他地区更容易积累。由千顷塘的DEM图(图1)可以看出,千顷塘地势呈现由西北向中部地区升高,再向西南方向整体降低的趋势,与土壤肥力指标空间分布特点较为相似,相关研究也证实了地形因子导致的太阳辐射不均匀以及水分差异是造成土壤属性空间异质性的重要原因[24—25]。

表3 土壤属性分级比例统计表

图3 千顷塘土壤肥力指标空间分布图

3.4 土壤肥力指标空间异质性影响因子

初步研究显示土壤肥力指标主要受结构性因素影响,其中地形因子可通过影响水热条件以及成土过程的再分配来影响土壤属性[26]。研究区地形因子和土壤属性的相关分析结果显示(表4),有机质和全钾与海拔有极显著的相关关系(P<0.01),其中有机质与海拔呈正相关,表明有机质含量随着海拔的升高而增加；全钾与海拔呈负相关,其含量体现出随海拔升高而下降的趋势,这可能与钾的元素特性有关。海拔较高的低山丘陵区,受水流的侵蚀作用较强, 相比于土壤对氮磷元素较强的吸持力来说,钾元素易于淋失[27—28],部分元素随水流淋溶下渗,出现其含量与海拔呈负相关的特点。本文中坡度与全钾含量呈正相关,但未达到显著水平(P>0.05),这与俞月凤等[29],杨家慧等[30]研究结果存在差异,多数研究结果显示海拔和坡度与钾素具有一致的相关性,本研究结果可能与千顷塘的地势特点有关。千顷塘高海拔地区分布着大面积的火山天池,天池周边地势相对平坦,而陡坡险坡在低海拔处较为常见,由此出现全钾与坡度呈正相关的结果。

表4 土壤肥力指标与环境因子的Pearson相关系数

不同海拔带间湿度、温度差异会影响土壤元素的迁移、分解与累积[31],为了进一步研究两种元素随海拔变化的规律,将海拔带划分为500—800、800—1000、1000—1200、1200—1400 m这4个等级(下文分别以1、2、3、4级海拔代替),对有机质和全钾在不同海拔带间的差异性进行分析,结果如图4所示。有机质在4个海拔带上的含量分别为49.66 g/kg、65.18 g/kg、79.29 g/kg、108.34 g/kg,1、2级海拔间,2、3级海拔和4级海拔间均存在显著差异性(P<0.05),表明有机质随海拔梯度变化较快；全钾在4个海拔带上的含量分别为27.50 g/kg、26.71 g/kg、19.39 g/kg、15.86 g/kg,1、2、3级海拔和4级海拔间存在显著性差异(P<0.05),相比有机质,全钾变化较慢,变化趋势较为平缓,与影响全钾的主控因素空间连续性尺度较小有关。除海拔之外,其他地形因子与土壤指标之间的相关性未达到显著水平。

图4 不同海拔梯度下的土壤肥力指标特征

容重与有机质、全氮和全钾均存在极显著相关性(P<0.01),表明土壤容重在很大程度上影响着土壤属性特征。图5显示,有机质和全氮含量随着容重增加而减少,说明两者与容重呈反比关系,这与祁凯斌等[32]研究结果一致。有机质下降趋势较大,含量变化较快；全氮下降趋势较小,总体变化平缓；全钾含量先增加后降低,变化幅度和全氮类似。研究显示,容重与土壤贮水能力、结构和松紧度密切相关,容重越大,表示土壤紧实度越高,透气透水性较差；反之,容重越小,土壤中团粒结构相对较多,表示土壤疏松多孔,质地较好,提高了水分入渗效率,有利于植物细根的穿插和发育,进一步提高了土壤有机质回归能力[33]。

图5 不同土壤容重下的肥力指标特征

3.5 不同植被类型土壤肥力指标差异性

土壤肥力很大程度上也受到植物体养分的影响,最直接影响因素是通过凋落物归还土壤的形式实现的。实地调查发现千顷塘自然保护区内共有落叶阔叶林、针阔叶混交林、针叶林、竹林这4种林地类型,高海拔处分布着较大面积的草甸,以及部分人为开垦的农田。对千顷塘不同植被类型下土壤肥力指标含量进行统计(图6)可知,不同林分间有机质、全氮、全磷的差异性显著,以有机质为例,有机质含量均值依次为针叶林(99.56 g/kg)>针阔叶混交林(77.00 g/kg)>草甸(75.17 g/kg)>落叶阔叶林(61.77 g/kg)>农田(45.97 g/kg)>竹林(40.45 g/kg),其中针叶林中土壤有机质含量显著高于其他林分土壤,原因在于,相比于其他林分类型,针叶林普遍分布在高海拔区域,对于有机质而言,随着海拔高度增加,气温和蒸发量下降,湿度变大,成土过程中的生物化学作用减弱,有机质易于积累,而分布在山顶区域的草甸养分含量偏低的主要原因在于植被分布较为单一,其枯落物归还量较小,这与谢红花等[34]对云南乌蒙山、薛丽佳等[35]对武夷山的有机质含量研究结果一致。与中亚热带其他地区的土壤林分肥力指标相比,千顷塘的土壤肥力指标处于较高水平,如黄继育等[36]报道了浙江省安吉县毛竹林土壤有机质平均含量为36.33 g/kg；Zeng等[37]对中亚热带典型常绿阔叶林的土壤营养状况进行了研究,其中常绿阔叶林中土壤有机碳含量为41.96 g/kg。可见原生性较好的森林土壤的保肥效果较好,同时自然保护区的生态结构功能也发挥了一定优势。对于全氮、全磷、全钾元素含量的统计发现,全氮在农田土壤中的含量最高,显著高于落叶阔叶林和竹林土壤；全磷在农田中的含量也显著高于其他利用类型的土壤；全钾含量在各类土壤中虽然不存在显著差异性,但在农田土壤中也处于较高水平,由此推断千顷塘保护区内开垦的农田存在一定的施肥耕作措施,有一定的人为扰动迹象。

图6 不同植被类型下土壤肥力指标特征

3.6 土壤肥力指标与环境因子的回归分析

采用回归分析方法定量表达海拔、容重以及植被类型对千顷塘土壤肥力指标的影响程度(表5),由于环境因子对pH的影响均不显著,因此回归分析中不纳入该指标。本文对定性分类变量采用哑变量[38]进行赋值后采用逐步线性回归进行分析,为重点反映林分类型对土壤属性的影响,哑变量赋值时以农田为参考类别,分析结果根据德宾-沃森(D-W)值检验自变量的自相关性。

表5 不同因素对土壤肥力指标的回归分析结果

分析结果中德宾-沃森(D-W)值均接近2,表示自变量基本不存在的一阶自相关性[39]。海拔、容重和林分类型对有机质的影响均达到极显著水平,与有机质在空间上变程较大,其空间异质性受环境因素影响强烈的研究结果相一致。三者对有机质的空间变异的独立解释能力分别为21.3%、51.5%、17.0%,表明容重是该地区影响有机质空间变异性的最主要因素,林分类型的逐步回归分析结果中入选因子为针叶林,说明针叶林对有机质的贡献率最大。全磷的变异仅与林分类型的影响达到显著水平(P<0.05), 独立解释能力为18.6%,这是由于土壤中的磷大部分来源于基岩的风化,而表层土壤磷直接来源于凋落物中的磷,并通过植物的表聚作用在表层土壤积累[40]。容重对全氮、全钾的独立解释能力分别为23.6%,17.9%,均略高于其他因素的影响,说明容重对这两种元素的影响较大,但是其解释度并不高,可见研究区土壤肥力指标还受到其他因素的影响。周晓阳等[41]研究显示土壤养分的空间变异与土壤类型相关,且分类级别越低,反应能力越大；曹详会等[42]对土壤中微生物的活性、氮磷钾等养分元素的淋溶迁移以及植物本身的生产力进行研究,从而发现气候因素影响土壤中养分分解与输送的途径；吕圣桥等[43]发现土壤有机碳空间分布特征与土壤黏粒含量呈正相关。人类活动对土壤肥力的影响日渐增加,叶晶等[44]研究发现人为的土地整理措施会降低土壤有机碳的质量分数,影响土壤生态系统的稳定性。值得注意的是,千倾塘保护区内已经出现越来越多的农田开垦现象,长此以往,会破坏自然保护区内原生生态结构,对林区土壤肥力失衡势必造成一定的影响,严重的或将危害到野生华南梅花鹿的生境。此外,大型动物对草地的采食和践踏也直接影响着土壤理化特性[45],在今后的分析中,可以补充考虑这一因素,建立更加科学、完善的指标体系。

4 结论

(1)研究区内pH、有机质、全磷、全氮、全钾含量平均值分别为5.12、72.44、0.45、1.91、23.16 g/kg,均属于中等程度变异。半方差分析结果显示pH、全磷和全钾表现为强烈空间自相关,有机质和全氮表现为中等空间自相关,说明土壤肥力指标主要受结构性因素影响。有机质、全磷、全氮的空间自相关性达到显著性水平,Moran′sI分别为0.21、0.19、0.29,在空间上呈聚集分布。

(2)Kriging插值结果显示多种指标在空间上呈现“两头高,中间低”的分布特点,高低值分布较为集中,斑块较大,分布较为连贯,与千顷塘地势特点存在一定的相似性。

(3)有机质、全磷、全氮、全钾与环境因子海拔、容重和林分类型的相关性达到显著水平,回归分析显示容重相对于其他因素对有机质、全氮和全钾的解释度最高,是影响这几种指标空间变异的主控因素,但研究因子总体解释度较低说明肥力指标的空间变异还受其他因素影响,在今后研究中要充分考虑。