浙江清涼峰自然保護區土壤肥力指標空間變異及其影響因素

姜霓雯,童根平,葉正錢,程樟峰,呂永強,傅偉軍,*

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%,高值區分布在中偏西北部,整體呈現明顯的“條帶狀”分布；全磷總體處于中等偏下水平,含量由西北向東南方向降低；全氮和全鉀總體含量較為豐富,全氮的分布格局和全磷較為相似,全鉀的高值區分布在西北部和東南部邊緣地區?？傮w來看,研究區各類土壤肥力指標含量均較高,只有磷素含量偏低,可能與土壤酸性較強有關,相關研究顯示,土壤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)有機質、全磷、全氮、全鉀與環境因子海拔、容重和林分類型的相關性達到顯著水平,回歸分析顯示容重相對于其他因素對有機質、全氮和全鉀的解釋度最高,是影響這幾種指標空間變異的主控因素,但研究因子總體解釋度較低說明肥力指標的空間變異還受其他因素影響,在今后研究中要充分考慮。