青海省不同生境下青海云杉胸徑生長模型研究

馬 浩，陳科屹，徐干君，黨禹杰，何友均，王建軍＊

(1. 國家林業和草原局西北調查規劃院，陜西西安 710048；2. 中國林業科學研究院林業科技信息研究所 北京 100091)

單木生長模型以單木為基本單位，通過數學模型來模擬林木的生長過程，揭示林木生長規律，并且能夠提供較詳細的林分結構特征及動態變化信息，對于森林經營管理和保護等方面有十分重要的實踐價值[1-3]。例如，李慧婷[2]采用5 種理論生長方程建立了青海海東地區華北落葉松的樹干材積生長模型。與傳統的非線性回歸方法相比，混合效應模型能夠有效提高模型的估計精度[3-8]。李春明和唐守正等[4]建立了考慮樣地的隨機效應、觀測數據的時間序列相關性及間伐強度的落葉松云冷杉林斷面積生長混合模型，結果表明混合效應模型顯著地提高了林分斷面積的預估精度。

單木生長模型在區域水平上的應用具有一定的局限性，其影響因素主要有樹木立地條件、人為經營措施、空間結構多樣性等[9-11]。張雷等[12]研究發現，隨海拔的升高，祁連山北坡青海云杉林的結構出現明顯變化，平均胸徑增加，平均樹高呈“單峰”變化。張中惠等[13]探究寧夏六盤山華北落葉松單木樹高生長對立地因子的響應，發現影響最大的立地因子是海拔，其次是坡向、坡度。海拔、坡向與坡位在地形上改變了溫度、濕度、光照等資源的分配，進而影響植被分布與林木生長發育[14]。因此，研究林木在不同環境下的生長特點，不僅對科學描述林木生長動態十分重要，還能為森林多種功能與效益的發揮提供基礎參考。

青海云杉（Picea crassifoliaKom.）是青海地區主要的森林更新樹種和造林樹種，在涵養水源、水土保持、改善環境等方面的作用顯著。國內外學者對該樹種的生理生態特征、樹輪生長特性等方面進行了廣泛研究[15-17]。但是，關于青海云杉單木生長模型的相關研究較少[2]。因此，本研究以青海云杉為研究對象，分析不同起源、不同坡位下單木胸徑生長的規律，構建單木胸徑生長模型和考慮起源和坡位影響的混合效應模型，旨在為今后有效保護、合理經營青海云杉林提供經驗模型和有益參考。

1 研究區概況

青海位于我國西北部，地理坐標為31°36'～39°19' N，89°35'～103°04' E，屬于高原大陸性氣候，冬季寒冷干燥，夏季短暫且涼爽濕潤，日照時間長，太陽輻射強。年平均氣溫-5.1～9.0℃，年均降水量15～750 mm，降水分布由東南向西北逐漸減少[18]。該地區的地形、地貌復雜多樣，地勢總體呈西北高東南低，東部多山，西部為高原和盆地，河流湖泊分布多，平均海拔3 000 m 以上[18-19]。主要的植被類型有常綠針葉林、落葉闊葉林、高寒灌叢、高寒草甸等，主要建群種有青海云杉、祁連圓柏（Juniperus przewalskiiKom.）、青楊（Populus cathayanaRehd.）、油松（Pinus tabuliformisCarr.） 等， 主要灌木有沙棘（Hippophae rhamnoidesLinn.）、檸條（Caragana korshinskiiKom.）、錦雞兒屬（CaraganaFabr.）等。

2 材料與方法

2.1 數據收集與處理

采樣點位于青海省3 個縣：尖扎縣、祁連縣、大通回族土族自治縣（簡稱“大通縣”）。采樣點選取受人為干擾少、遠離道路及防火道、環境較為穩定的地帶。在不同起源、不同立地條件下的典型林分中選取標準木作為樣木。其中，以坡位劃分，在不同海拔、坡向、坡度內進行采樣，保證數據的區域代表性和均勻性。于2022 年7 月采集青海云杉樣木共172 株。其中，天然起源的青海云杉74 株（上坡19 株，中坡9 株，下坡22 株，平坡24株），人工起源的青海云杉88 株（上坡29 株，中坡15 株，下坡28 株，平坡16 株）。在每1 株樣木的胸徑高（1.3 m）的位置，用生長錐分別在南、北2 個方向上鉆取樹芯，將樹芯樣本保存并分類編碼，記錄樣木胸徑（DBH）、樹種、起源、坡位等信息。

將樣本帶回實驗室后，按樹輪年代法的標準方法進行風干、固定、切割/打磨等預處理；使用LinTab 年輪分析儀對樹輪寬度測定，測定精度為0.01 mm。年輪條測量獲得的測量值為某一方向上的半徑生長量，對測量值的平均值乘以2，得到該年齡下的胸徑生長量。計算胸徑平均生長量（MAI）和胸徑連年生長量（CAI），齡階間距為5 a。樣本信息見表1。

表1 樣本信息 Table 1 Basic information of samples

此外，在典型天然林中選取4 株標準木進行樹干解析，記錄樹干長、胸徑以及圓盤的直徑與年齡，以確定林木地徑生長到胸徑位置時的年限。解析木的胸徑范圍為25.9～34.8 cm，樹齡范圍為69～135 a，樹高范圍為16.1～26.5 m，到達胸徑位置的年齡范圍為21～26 a。根據解析木數據分析，青海云杉地徑生長到胸徑位置時的年限為24 a。樹齡為樹芯的年輪數加上林木生長到胸徑位置時的年限。

2.2 研究方法

2.2.1 生長模型建立 結合樹木的生長規律及前人的研究成果，本研究選取理查德（Richards）方程、邏輯斯蒂（Logistic）方程、考爾夫（Korf）方程、坎派茲（Gompert）方程和單分子式（Mitscherlich）方程5 種樹木生長理論方程對胸徑生長量與年齡的關系曲線進行擬合[2,20]。利用ForStat 2.2 統計軟件擬合生長模型。

2.2.2 非線性混合效應模型建立 非線性混合效應模型的建立考慮了回歸函數依賴于固定和隨機效應的非線性關系[21]。多水平非線性混合效應模型的形式（以兩水平模型為例）如下：

式中，yijk是第i個第一水平中的第j個第二水平內的第k次觀察值，m、mi分別是第一水平、第二水平的分組數量，nij是第i個第一水平中的第j個第二水平內的觀測次數，f是含有參數向量φijk和協變量向量vijk的非線性函數，Aijk是設計矩陣，β是（p×1）維固定效應向量，Bi,jkbi、Bijkbij分別是第一水平、第二水平的隨機效應設計矩陣，bi、bij分別是第一水平、第二水平的隨機參數向量，D1、D2分別是第一水平、第二水平的隨機參數方差-協方差矩陣，bi和bij不相關，εijk是服從正態分布的誤差項，σ2是方差，Rij是第i個第一水平中的第j個第二水平內的方差-協方差矩陣。

對于參數效應的確定，本研究將所有不同隨機效應參數組合的模型都進行擬合。模型的擬合優度指標為赤池信息準則（AIC）、貝葉斯信息準則（BIC）和對數似然值（LogLik），一般采用最小AIC和BIC，以及最大Loglik值的標準確定具有最優參數組合的模型。隨機效應參數的方差-協方差矩陣D反映了隨機效應在個體之間的差異性。采用常見的廣義正定矩陣作為隨機效應參數的方差-協方差矩陣。方差-協方差矩陣Rij主要用于解決數據中存在的自相關和異方差問題[22]。由于胸徑生長量數據不是重復觀測數據，因此不考慮數據的自相關問題。通過加權回歸方法來消除異方差問題，常用的異方差函數有常數加冪函數、冪函數和指數函數。利用R 語言的nlme 包進行混合效應模型的參數估計。

2.2.3 模型評價 為了對模型的擬合結果進行評價和比較，采用的指標有決定系數（R2）和均方根誤差（RMSE）。采用全部數據計算模型估計值的精度指標，選擇總體相對誤差（TRE）、平均系統誤差（MSE）、平均預估誤差（MPE）和平均百分標準誤差（MPSE）4 項誤差指標對模型進行綜合評價。根據R2較大，RMSE、TRE、MSE、MPE和MPSE較小原則選取最優模型。

3 結果與分析

3.1 胸徑生長過程

以林木生長到胸高位置的年齡為起始年限（24 a），5 a 為齡階，分析不同起源、坡位生長下青海云杉的胸徑生長過程（圖1）。CAI與MAI均呈現“升—降—平緩變化”的趨勢，MAI大于CAI。其中，天然林內青海云杉的快速生長期為29—44 a，CAI與MAI均在0.40 cm 以上；29 a時CAI與MAI相交且有最大值，分別為0.509 cm、0.515 cm。人工林內青海云杉的速生期在29—39 a，CAI與MAI均大于0.40 cm，34 a 時CAI與MAI有最大值，分別為0.473 cm、0.434 cm；CAI與MAI出現多個相交點，MAI在39 a 快速下降后變化不大，CAI在樹齡164—259 a 波動較大，在244 a 出現第2 個峰值0.302 cm。

圖1 不同起源與坡位青海云杉單木胸徑生長曲線Fig. 1 Individual-tree DBH growth for Picea crassifolia at different origins and slope positions

天然林不同坡位中，青海云杉在上坡的CAI與MAI在29—84 a 均保持0.35 cm 以上的較高值，為生長速生期，并且出現2 個相交點；82 a之后CAI先降后升。在中坡的青海云杉速生期為29—49 a，其中CAI與MAI在29—39 a 均大于0.50 cm，29—94 a 之間波動幅度較大，出現3 處相交點。在下坡的青海云杉的CAI與MAI變化平緩，速生期為29—64 a，24—104 a 之間兩者出現3 個相交點。在平坡的青海云杉CAI與MAI逐年下降，速生期為29—54 a，其中CAI與MAI在29—44 a 均在0.50 cm 以上。

人工林不同坡位中，青海云杉MAI與CAI多次相交。在上坡的青海云杉CAI與MAI先降再平緩變化，生長速生期為29—39 a，CAI與MAI均在0.35 cm 以上；229—289 a 之間CAI的波動幅度較大。在中坡的青海云杉CAI與MAI變化平緩，但兩者均小于0.30 cm。在下坡的青海云杉CAI與MAI逐年下降后平緩變化；速生期為29—39 a；CAI波動幅度較大，159—204 a 時期CAI先升后降。在平坡的青海云杉呈現先升后降的變化趨勢，速生期在29—54 a，出現1 個相交點。

3.2 胸徑生長模型的擬合與評價

3.2.1 胸徑生長模型的擬合與評價 采用5 種生長模型擬合不同起源青海云杉單木胸徑生長量與年齡關系曲線，模型擬合效果良好（表2，圖2）。天然林單木生長模型中，Gompertz 模型的擬合結果表現最好，R2為0.915；對于人工林，Korf 模型的擬合效果最優，R2為0.946。2 個模型的TRE在 ± 1%以內，MSE在 ± 2%以內，MPE均在3%以內，MPSE在35%以內。

圖2 不同起源單木胸徑生長擬合曲線對比Fig. 2 The fitting curves for Individual-tree DBH growth models at different origins

表2 不同起源最優單木胸徑生長模型參數及評價指標Table 2 The parameter estimates and evaluation indices of optimal individual-tree DBH growth models at different origins

不同坡位青海云杉單木胸徑生長模型的擬合結果如表3、圖3～4 所示。上坡青海云杉單木胸徑生長擬合效果最優的生長模型是Richards 模型，中坡的最優生長模型是Logistic 模型，下坡的最優生長模型是Gompertz 模型，平坡的最優生長模型是Korf 模型。其R2分別為0.982、0.915、0.913、0.996，TRE和MSE在 ± 2%以內，MPE在6%范圍內，MPSE大部分在35%范圍內。人工林不同坡位下，上坡、中坡、下坡和平坡的最優模型分別是Korf 模型、Richards 模型、Korf 模型、Logistic模型。其R2分別為0.940、0.973、0.980、0.955，TRE和MSE在 ± 2%以內，MPE在5%范圍內，MPSE在35%范圍內。

圖3 不同坡位青海云杉天然林單木胸徑生長擬合曲線對比Fig. 3 The fitting curves for Individual-tree DBH growth model for Picea crassifolia natural forests at different slope positions

表3 不同坡位青海云杉最優單木胸徑生長模型參數及評價指標 Table 3 The parameter estimates and evaluation indices of optimal individual-tree DBH growth models at different slope positions

3.2.2 混合效應模型的擬合與評價 根據生長模型擬合結果，本研究分別以Gompertz 模型、Korf 模型和Logistic 模型作為基礎模型，采用全部數據，在基礎模型中考慮起源、坡位的混合效應，對所有不同隨機參數組合進行擬合和檢驗，各基礎模型的最優模型結果見表4。由表可知，混合效應模型的檢驗結果均優于基礎模型，考慮混合效應顯著提高了模型的擬合效果。Gompertz 模型、Korf模型、Logistic 模型的最優隨機參數組合分別為(a、c)、(b)、(b、c)。以Korf 模型為基礎模型的混合效應模型具有最小的AIC、BIC和最大的LogLik，但是Korf 模型在考慮異方差函數過程中不收斂。因此，選擇擬合結果其次的Gompertz 模型（AIC=12 700.520，BIC= 12 757.610，LogLik= -6 340.259），并進行異方差校正。由表4 可知，3 種異方差函數均能改善模型的擬合效果。其中，添加冪函數的混合效應模型的擬合結果最優，AIC為12 394.490，BIC為12 457.300，LogLik為-6 186.247。

表4 混合效應模型的模擬結果 Table 4 Fitting results of mixed-effect model

青海云杉胸徑生長基礎模型和混合效應模型的參數估計以及擬合結果如表5 所示。從表中可以看出，最優混合效應模型具有更好的擬合精度，其R2為0.702，較基礎模型提高了32.7%；RMSE為4.14，較基礎模型下降了20.8%。十折交叉驗證法檢驗結果顯示，混合效應模型的檢驗指標均優于基礎模型?；旌闲Ｐ蚑RE、MSE、MPE、MPSE的降幅分別為89.3%、83.5%、20.6%、15.1%。

3.2.3 模型的預測 利用最優模型對青海云杉胸徑的生長擬合見圖5。由圖可以看出，青海云杉胸徑總生長量隨年齡逐漸增大，其中天然林上坡和人工林下坡的單木胸徑生長趨勢較好。天然林整體、天然林中坡和天然林下坡的單木連年生長量CAI先增大后減小，平均生長量MAI逐年下降；天然林上坡與平坡的單木CAI和MAI均呈現逐年下降趨勢。除了人工林平坡，人工林的單木CAI和MAI均逐年下降。

圖5 青海云杉單木胸徑生長擬合曲線Fig. 5 Individual-tree DBH growth curve for Picea crassifolia

4 討論

總體來看，青海云杉生長到胸高位置后（24 a），單木胸徑生長量隨著年齡的增加而逐漸減少（圖1）。青海云杉生長緩慢，樹齡20 a 左右的林木樹高約1.4 m[23-24]。溫度、降水、光、空間等是該地區林木生長的限制因素[25-27]。隨著胸徑生長量的增大，林木個體在林分中能夠獲取的可利用資源變少，養分吸收的有效性降低，使得林木在后期生長緩慢。不同起源對比，天然起源的青海云杉在樹齡29—44 a 保持相對較高的單木胸徑生長量，CAI和MAI均在0.40 cm 以上，比人工林的持續時間長（29—39 a）；隨后CAI和MAI的變化趨勢平緩（圖1）。在西北地區已有類似的研究結果[28-30]。王學福[29]研究發現祁連山青海云杉胸徑連年生長量高峰期在40—50 a，平均生長量高峰期在90—110 a。楊文娟[30]指出青海云杉單木胸徑的最大連年生長量出現在40 a 時，最大平均生長量出現在70 a 時，分別為0.30 cm、0.22 cm。不同坡位條件下，胸徑的生長規律具有差異（圖1）。坡位上的差異主要表現在土壤肥力和土壤水分方面，其中低坡位能夠最有效地提供林木生長所需的養分，進而促進林木的徑向生長[31-33]。蘇妮爾等[32]指出下坡位的土壤養分含量最高，是大徑材紅皮云杉培育的最優坡位。

單木胸徑生長模型的擬合結果顯示，不同生境壓力下林木生長過程具有不同的理論生長方程，但是大部分都符合傳統的“S”形生長曲線（除了Mitscherlich 模型）。各最優胸徑生長模型的R2均在0.913 以上，TRE和MSE均在 ± 2%以內，MPE大多在5%以內，MPSE在35%以內（表2～3）。這與同一地區生長模型研究得到的結果接近。馬克西等[34]構建了青海省4 個主要樹種組的胸徑生長率模型，其中模型的R2均在0.93 以上，MPE均在0.3%之內；胸徑生長量預估的MPE均在2%之內，MPSE均在60%左右。

在基礎模型上考慮起源和坡位對胸徑生長量的影響，建立了青海云杉單木胸徑生長量混合效應模型，結果顯示混合模型的擬合效果均優于基礎模型（表4～5）。在模型檢驗中，混合模型仍然呈現較高的擬合精度，其TRE、MSE、MPE和MPSE較基礎模型分別下降89.3%、83.5%、20.6%、15.1%。這與大多數學者的研究結論一致[4,6,8,21,26]。此外，模型擬合的精度還可能受到海拔、溫度、降水、干旱等因素的影響[12,17,35]。因此，未來氣候變化下的森林管理與保護應充分考慮更多的地形、氣候差異對林木生長的影響。

5 結論

不同生境壓力下的青海云杉胸徑生長有不同的變化規律。單木胸徑生長量隨著年齡的增加而下降，其中天然林的速生期較人工林的持續時間長，不同坡位條件下胸徑生長規律具有差異。青海云杉單木胸徑混合效應模型的擬合效果優于基礎模型。因此，在未來的森林經營及保護中，這些模型適用于青海省不同生境壓力下青海云杉的林木胸徑生長量估算。