天山雪嶺云杉林土壤生態化學計量特征的空間分布模擬

熱依法提?艾尼瓦爾 許仲林 常亞鵬 李路

摘要：土壤碳（C）、氮（N）和磷（P）含量以及它們之間的生態化學計量特征在養分限制、群落動態、養分利用效率和生物地球化學循環中起著重要作用，為了深入探討陸地生態系統的物質循環，使用多元線性回歸（MLR）、逐步回歸（STR）、嶺回歸（RDR）、套索回歸（LSR）方法估算天山雪嶺云杉林土壤C、N、P含量以及其化學計量比的空間分布。結果表明：（1）雪嶺云杉林土壤C、N、P含量和C ∶N、C ∶P、N ∶P 的均值分別為（55.76±25.24）、（4.63±2.11）、（0.71±0.21） g/kg和（12.24±3.73）、（79.95±31.89）、（28.16±3.18）;（2）隨經度的增加，雪嶺云杉林內土壤C、P含量，C ∶P、C ∶N遞減，N 含量、N ∶P遞增;（3）對雪嶺云杉林土壤N、P含量及其化學計量比空間格局的模擬中MLR、STR的表現優于RDR、LSR;（4）相對于原始變量，使用主成分作為預測因子的估算并未提高預測結果的精度。建議在未來的工作中更多地考慮自變量（如生物、非生物、人為因素等）以提高空間模擬的精度。

關鍵詞：C;N;P;化學計量比;空間分布;元素含量;雪嶺云杉;天山山脈

中圖分類號：S714?? 文獻標志碼： A

文章編號：1002-1302（2021）17-0208-08

收稿日期：2021-01-19

基金項目：新疆維吾爾自治區高校科研計劃（編號：XJEDU2020Y008）。

作者簡介：熱依法提·艾尼瓦爾（1995—），男，新疆克拉瑪依人，碩士研究生，主要從事生態化學計量學相關研究。E-mail：rayifat@163.com。

通信作者：許仲林，博士，教授，主要從事GIS及環境建模方面的研究工作。E-mail：zlxu@xju.edu.cn。

生態化學計量學側重于研究養分比例的變異性及其影響因素，已被成功應用于養分限制[1-2]、種群動態[3]、微生物養分狀況[4]、共生關系[5]、陸地及水生生態系統的養分利用效率[6]、生物地球化學循環[7-8]等領域。土壤是全球生物地球化學循環的重要組成部分，土壤碳（C）、氮（N）和磷（P）的空間格局及其化學計量特征是生物地球化學循環模型的關鍵輸入變量，也是研究陸地生態系統對全球環境變化反饋的重要參考[9]。目前，相關研究主要集中在土壤C、N、P的垂直模式（土壤深度梯度、海拔梯度）以及植物生長、凋落物分解、土壤風化過程中的生態化學計量特征變化[9-11]，而在生態化學計量的空間分布方面，仍須深入考察。

雪嶺云杉（Picea schrenkiana）林是天山山脈中典型的亞高山森林（海拔1 600～2 800 m），分布范圍從烏茲別克斯坦天山西部延伸到中國新疆天山東部，該物種分布范圍的廣域性和獨特性為研究中亞干旱區山地森林生態系統土壤C、N、P含量及其化學計量比的空間分布提供了絕佳對象。目前，對雪嶺云杉林土壤生態化學計量特征的研究較少，如陳曦等研究結果顯示，受年平均降水量（MAP）限制，雪嶺云杉林土壤養分含量與其他地區相比相對較低[12];Dai等考察了雪嶺云杉林線和土壤養分特征的空間變異性，研究表明天山中部地區的C儲量相對較高，而西部地區的N和P總量相對較高[13]。可以看出，對雪嶺云杉林生態系統生態化學計量特征的考察，仍然相對缺乏。基于此，本研究在實地采樣和和實驗室分析的基礎上，構建土壤C、N、P含量及其化學計量比等因變量與海拔、年平均溫度、年平均降水量等自變量之間的回歸模型，并在此基礎上應用模型估計雪嶺云杉林土壤及其化學計量比特征的空間分布特征。以期為估計山區復雜地形條件下土壤生態化學計量特征的空間分布提供借鑒。

1 材料與方法

1.1 研究區

研究區位于天山山脈的中部和東部區域（42°35′～44°20′ N，80°14′～88°07′E，圖1），屬典型大陸性氣候，冬季寒冷干燥、夏季溫暖潮濕。由于獨特的地理位置和地形特征，研究區溫度日較差高于周邊地區，隨海拔升高，年平均溫度從13.3 ℃降低到-7.3 ℃，年平均降水量從小于100 mm增加到500 mm以上[14]。一般情況下，森林中的土壤在11月下旬開始凍結，在第2年4月初開始解凍。夏、冬季2季相對較長，春、秋2季較短。天山的植被類型在海拔梯度上由低至高分別為草原、森林草原、亞高山灌叢草甸、永久冰雪[14]。雪嶺云杉林是海拔 1 600～2 700 m之間森林草原帶的純林。在森林邊緣生長的灌木物種有黑果栒子（Cotoneaster melanocarpus）、異果小檗（Berberis heteropoda）、密刺薔薇（Rosa spinosissima）、金絲桃葉繡線菊（Spiraea hypericifolia）、新疆方枝柏（Juniperus pseudosabina）、白皮錦雞兒（Caragana leucophloea）、剛毛忍冬（Lonicera hispida）等。林下草本植物包括圓葉老鸛草（Geranium rotundifolium）、天山羽衣草（Alchemilla tianschanica）、羊角芹（Aegopodium podagraria）[15]。

1.2 樣本采集與處理

本研究采樣工作在2012—2018年進行，共選擇173塊樣地進行土壤取樣（圖1），位置信息由全球定位系統（GPS）記錄。在每塊樣地中，隨機選擇5個采樣點，每個采樣點分別采集深度為0～30、30～50、50～80 cm 的土壤，5個樣點同一深度的土壤樣本經混合獲得均值樣本。將土壤樣本帶回實驗室經標準化處理獲得待分析土樣，后分別使用重鉻酸鉀氧化-外加熱法[16]、凱氏定氮法[17]、高氯酸消化法[18]分別測定土壤C、N、P含量。

1.3 回歸建模

本研究利用的氣候數據為WorldClim生物氣候數據集[19]，該數據基于全球超過40 000個氣象站的觀測數據，使用薄板樣條插值算法[19]得到19個溫度和降水變量，已在不同類型的空間建模研究中得到廣泛應用[20]。本研究首先對數據集中的變量進行重采樣（原始空間分辨率為1 km×1 km，將其重采樣得到30m×30m的分辨率）以獲得適用于

山區的氣候變量。除了WorldClim數據集中的變量外，本研究建模中使用的變量還包括海拔、土壤有效含水量（簡稱AWC，體積分數）、土壤容重（簡稱SBD，g/cm3）、陽離子交換量（簡稱CEC，cmol/kg）、沙粒含量（簡稱SDC，質量百分比）、粉粒含量（簡稱CRF，質量百分比）、黏粒含量（簡稱CLC，質量百分比）、土壤pH值，這些變量提取自SoilGrids數據集[21]。考慮到使用全部26個變量可能增加模型的復雜性，因此本研究采用主成分分析（PCA）實現降維處理。

雪嶺云杉林的實際分布來自于中國西部環境與生態科學數據中心。首先，在地理信息系統平臺中將原始矢量數據轉換為柵格數據，然后，將其重采樣為30 m×30 m分辨率的數據。本研究使用的回歸模型包括多元線性回歸（MLR）、逐步回歸（STR）、嶺回歸（RDR）、套索回歸（LSR）。MLR構造了因變量和自變量之間的回歸關系：Y=β0+β1X1+β2X2+…+βkXk+ε，其中，Y為因變量（C、N、P含量及其化學計量比），Xk為前述26個環境變量（或它們的主成分），βk是自變量Xk的回歸系數，ε是誤差。STR通過引入或排除解釋變量構建回歸模型，判斷變量的選擇標準是否提高了模型的擬合精度。嶺回歸為最小二乘法的補充，建模過程匯中增加了正則化選項。與嶺回歸類似，套索回歸也增加了正則化選項并對回歸系數的絕對值之和進行了界定[22]。使用SPSS 22.0確定多元線性回歸、逐步回歸、嶺回歸的回歸系數，套索回歸的系數由R語言的glmnet包獲得。對以上模型性能的評價，使用平均絕對誤差（MAE）、均方根誤差（RMSE）、調整的決定系數（Adjusted R2）實現。建模過程中，如果回歸模型能夠通過0.05水平的顯著性檢驗，則認為模型性能是顯著的。

2 結果與分析

2.1 取樣點土壤C、N、P含量及其生態化學計量特征

雪嶺云杉采樣點C、N、P含量分別為3.12～139.94、1.78～16.23、0.37～1.46 g/kg，平均值分別為（55.76±25.24）、（4.63±2.11）、（0.71±0.21） g/kg;C ∶N、C ∶P、N ∶P分別為 1.33～31.15、2.26～172.70、1.99～28.16，平均值分別為（12.24±3.73）、（79.95±31.89）、（28.16±3.18）。

2.2 預測變量的主成分

根據表1可知，前4個主成分占所有變量變異性的85.93%。其中，成分1中溫度變量的載荷較高，主要描述與溫度相關的預測因子;成分2在降水變量上載荷較高，主要表征了降水量的變異性;成分3在土壤沙粒含量、粉粒含量和黏粒含量上的載荷較高，主要描述了土壤的粒度特性。

2.3 C、N、P化學計量比和環境變量之間的模型性能和相關性

如表2所示，C、N、P含量空間分布模型的F值（1.93～19.04）表明，4種模型具有較高精度。相對而言，模型對C、N、P化學計量比的模擬精度較低。此外，F值和調整的R2表明，MLR模型的表現優于其他模型。就變量選擇而言，利用初始28個變量建模得到的回歸模型（調整的R2為 0.05～0.51）優于主成分回歸模型（調整的R2為0.02～0.33）。

如表3所示，Pearson相關系數矩陣顯示土壤C含量與最暖月的最高溫度、氣溫年較差、最暖季平均溫度、季節性降水變化具有相對較高的相關性;N含量與最暖月最高溫度、最濕季平均溫度之間的相關性較高;平均氣溫日較差是唯一與C ∶N的相關系數絕對值>0.2的變量。海拔、最暖月最高溫度、溫度年較差、最濕季平均溫度、最暖季平均溫度和季節性降水變化與C ∶P具有相對較高的相關系數。N ∶P與海拔、年平均溫度、最暖月最高溫度、最濕季平均溫度和最暖季平均溫度相關性較高。

2.4 土壤C、N、P含量及其化學計量比的空間格局

根據圖2可知，研究區模擬的C、N、P含量范圍分別為45.12～75.97、3.16～6.40、0.45～1.40 g/kg。與實際值的變異性相比，模擬值的變化范圍較小。由圖3可知，預測的C ∶N、C ∶P、N ∶P分別為8.64～14.03、59.29～112.43、5.86～9.28。

與中部、西南部相比，研究區東部和西北部雪嶺云杉林土壤的C含量較低。N含量呈現出與C含量相同的趨勢：研究區東部區域相對較低，而西部和中部則相對較高;而P含量的空間分布趨勢從研究區西部到中部遞增，從中部到東部呈現遞減趨勢。C ∶N在研究區西部變異程度較高，在整個研究區呈現出自西向東的遞增趨勢;C ∶P在西部較高，中部和東部較低;N ∶P呈現出與C ∶P相似的趨勢（圖3）。

2.5 實際值和模擬值間的比較

如圖4可知，實測C含量和模擬C含量間的相關系數較高（0.68，P<0.01，圖4-a）。相對而言，P含量的實測值和模擬值之間的相關系數較低（0.38，P<0.01，圖4-c），N含量的實測值和模擬值之間的相關系數（0.61，P<0.01，圖4-b）介于二者之間。C ∶N實測值和模擬值之間的相關系數為0.64（P<0.01）;C ∶P、N ∶P的實測值和模擬值之間的相關系數分別為0.39、0.44，相對較低，表明模型對C ∶N的模擬精度較高，而對C ∶P和N ∶P的模擬效果較差（圖4-e和4-f ）。

3 結論與討論

3.1 雪嶺云杉林中C、N、P含量及其化學計量比的空間格局

對土壤生態化學計量特征空間格局的準確認識有助于深入了解區域生物和生態過程的機制。本研究發現，C、P含量和N ∶P從研究區西部到東部遞增，而N含量、C ∶N和C ∶P沿西-東方向遞減（圖2、圖3）。 這與Dai等土壤有機碳自西向東遞增，而總氮量自西向東遞減的研究結果相似;該研究還指出，雪嶺云杉林土壤養分含量的規律性是由非生物變量的變化引起的，包括溫度、降水、蒸發能力、太陽輻射[13]。本研究結果也表明，溫度、降水等土壤理化性質的非生物因子與雪嶺云杉林土壤生態化學計量特征之間關系顯著;另外，不同的土壤生態化學計量特征環境變量的貢獻是不一致的，表明溫度、降水以及土壤理化性質在土壤生態化學計量特征的變異性方面具有不同的作用。關于雪嶺云杉林土壤生態化學計量特征的空間變異性方面，Cui等也指出與研究區西部相比，東部土壤C和P含量相對較低[22-23]。綜上所述，前人研究與本研究均表明，雪嶺云杉林內土壤C、N、P含量及其化學計量比在非生物因子的作用下，呈現出空間上的連續變化趨勢。

3.2 回歸模型的可靠性

生態化學計量特征在養分限制[1-2]、養分利用效率[6]、種群動態[3]、共生關系[5]以及區域和生物地球化學循環[7-8]中起著重要作用。因此，確定生態系統養分元素及其生態化學計量特征與環境變量之間的關系，并建立模型對生態化學計量特征的空間分布進行準確預測具有重要的科學意義。目前，學者們多采用相關分析和線性回歸分析方法進行研究[24-26]。但是，相關方法均是基于樣點數據的相關和回歸方法，旨在討論生態化學計量特征對環境因子的響應關系，并未給出生態化學計量特征的空間分布。盡管有學者嘗試使用地統計學和空間插值方法獲得了化學計量特征的空間分布[27-28]，但空間插值方法在地形復雜的地區（例如山區），可能受制于樣本量、環境變量的空間分辨率等因素而得到精確度有限的插值結果。最近，有學者提出生態化學計量分布模型（StDMs），并應用該模型模擬養分元素的空間分布和地表植物對該空間分布的響應，結果表明StDMs具有較高精度[29]。與StDMs類似，本研究所采用的方法也適用于在地形復雜區域給出較為可靠的生態化學計量特征的空間分布。

3.3 預測變量選擇

本研究中土壤pH值、陽離子交換量、土壤有效含水量等變量與化學計量特征之間的相關系數相對較低，表明氣候變量可能是大尺度土壤生態化學計量特征變異性的控制因素，這也在其他研究中得到了證實[30-31]。理論上，溫度和降水主要決定植物的功能類型，并控制陸地生態系統中植物的生物量積累[32]。此外，氣候變量也會對微生物活動產生影響，這對于有機物的分解至關重要[33-34]。目前，生態模型中最常用的氣候變量主要包括溫度和降水的平均值（例如年平均溫度和年平均降水量），但本研究表明，季節性溫度變化（例如溫度和降水的年較差）以及極值（例如最暖月最高溫度、最濕季降水量等）等因子，與土壤生態化學計量特征間也可能存在一定的相關關系，因此也是相關模型需要重點考慮的因子。值得注意的是，氣候變量之間可能是線性相關關系[19，35]，因此，需要使用合適的方法來找到具有最低共線性的數據集。本研究使用逐步回歸和主成分分析對共線性進行了處理，研究結果表明，就變量選擇而言，逐步回歸的性能優于主成分分析。

3.4 回歸方法的優點和局限

本研究發現回歸方法具有以下優勢：（1）使用該模型可以可靠地量化養分含量及其生態化學計量比對自變量的響應;（2）可以借助地理信息系統技術實現對生態化學計量特征空間分布的模擬;（3）模擬結果具有較高的精度。需要指出的是，其他模型在類似的應用中也具有潛力。例如，人工神經網絡（ANN）能夠檢測相關變量和自變量之間的復雜非線性關系[36];多元自適應回歸樣條（MARS）可以擬合變量之間的相互作用[37];分類回歸樹（增強的回歸樹、隨機森林等）也是分析復雜生態數據具有價值的方法[38];遺傳算法（GA）是一種全局優化方法，可用于解決復雜的優化問題[39]。因此，在后續工作中，綜合考慮多種模型有助于提高生態化學計量特征空間分布模擬結果的精度。

參考文獻：

[1]Feller I C，Mckee K L，Whigham D F，et al. Nitrogen vs. phosphorus limitation across an ecotonal gradient in a mangrove forest[J]. Biogeochemistry，2003，62（2）：145-175.

[2]Hgberg P，Nsholm T，Franklin O，et al. Tamm review：on the nature of the nitrogen limitation to plant growth in Fennoscandian boreal forests[J]. Forest Ecology and Management，2017，403：161-185.

[3]Johnson M T，Agrawal A A. Plant genotype and environment interact to shape a diverse arthropod community on evening primrose （Oenothera biennis）[J]. Ecology，2005，86（4）：874-885.

[4]Hill B H，Elonen C M，Seifert L R，et al. Microbial enzyme stoichiometry and nutrient limitation in US streams and rivers[J]. Ecological Indicators，2012，18（4）：540-551.

[5]Mariotte P，Canarini A，Dijkstra F A. Stoichiometric N：P flexibility and mycorrhizal symbiosis favour plant resistance against drought[J]. Journal of Ecology，2017，105（4）：958-967.

[6]He W M，Yu F H，Zhang L L. Physiological integration impacts nutrient use and stoichiometry in three clonal plants under heterogeneous habitats[J]. Ecological Research，2010，25（5）：967-972.

[7]Midgley M G，Phillips R P. Resource stoichiometry and the biogeochemical consequences of nitrogen deposition in a mixed deciduous forest[J]. Ecology，2016，97（12）：3369-3378.

[8]Schmidt S K，Porazinska D，Concienne B L，et al. Biogeochemical stoichiometry reveals P and N limitation across the post-glacial landscape of Denali National Park，Alaska[J]. Ecosystems，2016，19（7）：1164-1177.

[9]Yang Y H，Fang J Y，Guo D L，et al. Vertical patterns of soil carbon，nitrogen and carbon：nitrogen stoichiometry in Tibetan grasslands[J]. Biogeosciences Discuss，2010，7（1）：1-24.

[10]Jobbágy E G，Jackson R B. The vertical distribution of soil organic carbon and its relation to climate and vegetation[J]. Ecological Applications，2000，10（2）：423-436.

[11]王 棣，耿增超，佘 雕，等. 秦嶺典型林分土壤有機碳儲量及碳氮垂直分布[J]. 生態學報，2015，35（16）：5421-5429.

[12]陳 曦，許文強，羅格平，等. 天山北坡不同環境條件下雪嶺云杉（Picea schrenkiana）林限土壤屬性[J]. 生態學報，2008，28（1）：53-61.

[13]Dai L，Li Y Z，Luo G P，et al. The spatial variation of alpine timberlines and their biogeographical characteristics in the northern Tianshan Mountains of China[J]. Environmental Earth Sciences，2013，68（1）：129-137.

[14]Li L，Chang Y P，Li X F，et al. Carbon sequestration potential of cropland reforestation on the northern slope of the Tianshan Mountains[J]. Canadian Journal of Soil Science，2016，96（4）：461-471.

[15]王慧杰，常順利，張毓濤，等. 天山雪嶺云杉森林群落的密度制約效應[J]. 生物多樣性，2016，24（3）：252-261.

[16]Yeomans J C，Bremner J M.A rapid and precise method for routine determination of organic carbon in soil[J]. Communications in Soil Science & Plant Analysis，1988，19（13）：1467-1476.

[17]Bremner J M，Tabatabai M A. Use of an ammonia electrode for determination of ammonium in Kjeldahl analysis of soils[J]. Communications in Soil Science & Plant Analysis，1972，3（2）：159-165.

[18]Sherman M. Colorimetric determination of phosphorus in soils. Provision for eliminating the interference of arsenic[J]. Industrial & Engineering Chemistry Analytical Edition，1942，14（2）：182-185.

[19]Hijmans R J，Cameron S E，Parra J L，et al. Very high resolution interpolated climate surfaces for global land areas[J]. International Journal of Climatology，2010，25（15）：1965-1978.

[20]Fick S E，Hijmans R J. WorldClim 2：new 1‐km spatial resolution climate surfaces for global land areas[J]. International Journal of Climatology，2017，37（12）：4302-4315.

[21]Hengl T，Mendes D J J，Heuvelink G B M，et al. SoilGrids250m：global gridded soil information based on machine learning[J]. PloS One，2017，12（2）：e0169748.

[22]Cui W，Li Z，Chang Z，et al. Soils in Xinjiang[M]. Beijing：Science Press，1996.

[23]Guisan A，JrT C E，Hastie T. Generalized linear and generalized additive models in studies of species distributions：setting the scene[J]. Ecological Modelling，2002，157（2/3）：89-100.

[24]gren G I，Wetterstedt JM，Billberger M F K. Nutrient limitation on terrestrial plant growth-modeling the interaction between nitrogen and phosphorus[J]. New Phytologist，2012，194（4）：953-960.

[25]Feng D F，Bao W K，Pang X Y. Consistent profile pattern and spatial variation of soil C/N/P stoichiometric ratios in the subalpine forests[J]. Journal of Soils & Sediments，2017，17：1-12.

[26]Zinke P J，Stangenberger A G. Elemental storage of forest soil from local to global scales[J]. Forest Ecology & Management，2000，138（1/2/3）：159-165.

[27]Smith D B，Cannon W F，Woodruff L G，et al. Geochemical and mineralogical maps for soils of the conterminous United States[R]. Reston：U.S. Geological Survey，2014.

[28]Wang Y Q，Zhang X C，Huang C Q. Spatial variability of soil total nitrogen and soil total phosphorus under different land uses in a small watershed on the Loess Plateau，China[J]. Geoderma，2009，150（1/2）：141-149.

[29]Leroux S J，Wal E V，Wiersma Y F，et al. Stoichiometric distribution models：ecological stoichiometry at the landscape extent[J]. Ecology Letters，2017，20（12）：1495-1506.

[30]Elser J J，Sterner R W，Gorokhova E，et al. Biological stoichiometry from genes to ecosystems[J]. Ecology Letters，2000，3（6）：540-550.

[31]Tian H，Chen G S，Zhang C，et al. Pattern and variation of C ∶N ∶P ratios in Chinas soils：a synthesis of observational data[J]. Biogeochemistry，2010，98（1/2/3）：139-151.

[32]de Long J R，Sundqvist M K，Gundale M J，et al. Effects of elevation and nitrogen and phosphorus fertilization on plant defence compounds in subarctic tundra heath vegetation[J]. Functional Ecology，2016，30（2）：314-325.

[33]Cleveland C C，Liptzin D. C ∶N ∶P stoichiometry in soil：is there a “Redfield ratio” for the microbial biomass？[J]. Biogeochemistry，2007，85：235-252.

[34]Delgado-Baquerizo M，Reich P B，Khachane A N，et al. It is elemental：soil nutrient stoichiometry drives bacterial diversity[J]. Environmental Microbiology，2017，19：1176-1188.

[35]Kriticos D J，Webber B L，Leriche A，et al. CliMond：global high-resolution historical and future scenario climate surfaces for bioclimatic modelling[J]. Methods in Ecology & Evolution，2012，3：53-64.

[36]Martiny A C，Vrugt J A，Primeau F W，et al. Regional variation in the particulate organic carbon to nitrogen ratio in the surface ocean[J]. Global Biogeochemical Cycles，2013，27（3）：723-731.

[37]Lek S，Guégan J F. Artificial neural networks as a tool in ecological modelling，an introduction[J]. Ecological Modelling，1999，120（2/3）：65-73.

[38]DeAth G，Fabricius K E. Classification and regression trees：a powerful yet simple technique for the analysis of complex ecological data[J]. Ecology，2000，81（11）：3178-3192.

[39]Hamblin S. On the practical usage of genetic algorithms in ecology and evolution[J]. Methods in Ecology and Evolution，2013，4（2）：184-194.