臺州市喬木林樹種結構變化模擬與預測

任典挺，張 軍，金 鑫，李 萍，王懿祥，葛宏立，陶吉興

（1.三門縣自然資源和規劃局，浙江 三門 317100；2.臺州市黃巖區林業局，浙江 臺州 318020；3.浙江農林大學 環境與資源學院，浙江 杭州 311300；4.浙江省森林資源監測中心，浙江 杭州 310020）

森林作為獨有的陸地生態系統，發揮著無可替代的作用與功能[1]。森林具有運動和變化的特征，隨著時間的推移，樹種組成會發生明顯的變化，從而引起整個森林組分的變化，這一現象稱為森林演替[2]。森林演替中森林組分的變化除了與木材的供應有關外，還與生產力、碳固存、動植物生境等有關，因此，對森林演替進行預測，在森林生態系統管理中尤為重要[3]。

長期以來生態學一直在深入探究和預測森林演替[4-5]。研究者使用高度參數化的模型從模擬葉片的光合作用、蒸騰作用等方面來解釋森林動態演替進程[6-7]。較多研究通過單木生長、光照、溫度、養分、水分、更替及死亡的林窗模型描述森林演替過程[2,8-9]。較大空間尺度有林分級混合模型、景觀級模型等，隨著模擬單位尺度的變大，林分水平、樣地水平的模型也得到更多的認識，并且伴隨著空間分辨率和時間間隔的增加，可以使用相對較少的參數模擬森林動態[10]。此外，也有馬爾科夫鏈的概率轉移理論考慮了外部擾動、空間異質性等因素描述森林轉移的過程[11]。

以單木為基本單元來描述森林是最為標準的，但不適用于大面積區域，此類模型往往伴隨著大量反映局部地塊特征的參數，能夠很好地模擬試驗區的森林動態變化，但是過度復雜的模型如果沒有大量的工作準備，對新的應用環境并不友好[2,12-13]。森林資源連續清查結果表明：因受松材線蟲Bursaphelenchusxylophilus病的危害，臺州市松類資源大幅減少，受經濟因素影響，杉類資源也在減少，同時闊葉類資源在大幅增加。但通過建立數學模型定量描述并預測這種變化的研究鮮見報道。現有模型過于依賴細節和大量變量，并不適用于大區域樹種結構變化情況的模擬[12-13]。鑒于此，本研究對臺州市森林樹種結構變化的特殊情況及其變化缺乏適用模型的情況進行了研究，以期為臺州市森林樹種結構調整提供科學依據。

1 研究地區與方法

1.1 研究區概況

臺州市(28°01′～29°20′N，120°17′～121°56′E)位于浙江省中部沿海，總面積約10 050 km2。該區屬中亞熱帶季風氣候，年均日照時數為1 800～2 037 h，全年降水日數為132～171 d，年均降水量為1 632 mm。夏、秋季受臺風影響比較大。地帶性植被為中亞熱帶常綠闊葉林，主要的闊葉喬木樹種有殼斗科Fagaceae 的甜櫧Castanopsiseyrei、山茶科Theaceae 的木荷Schimasuperba等，針葉樹種主要有馬尾松Pinusmassoniana、杉木Cunninghamialanceolata等。主要森林類型有常綠闊葉林、落葉闊葉林、針闊混交林、常綠落葉闊葉混交林、竹林等[14]，森林覆蓋率為67.0%。

1.2 研究數據

選用臺州市6 期(1994、1999、2004、2009、2014 和2019年)森林資源連續清查固定樣地和固定樣木的資料作為研究數據，各期固定樣地總數(含非林地、無林地等非喬木林樣地) 分別為379、390、388、388、388 和388 個。1994 和1999年2 期樣地數和后面各期相比略有差異，是因為前3 期每期有1/3 的樣地進行了移動。移動部分樣地的目的是為了檢驗固定樣地是否受到特殊對待，即樣地內部和外部是否具有相同的經營措施，結果表明沒有特殊對待的情況發生，即樣地固定后對總體的代表性不會改變。固定樣地為800 m2的正方形。森林資源連續清查體系是基于概率抽樣理論的機械抽樣，每5 a 復測1 次，是森林資源調查中最具權威性的數據，對于森林的動態變化分析及預測十分有利[15]，同時也便于更有效地傳播和利用現有數據進行可持續的森林生態系統經營和管理[16]。本研究的對象是喬木樹種，所以只選擇其中的喬木林樣地，根據地類選擇樣地不改變其隨機性和代表性。由于經濟林比較特殊，所以本研究不包括喬木經濟林。

由于樹種眾多，對每一樹種進行模擬是不現實的，因此需要將具有共同特性的樹種聚集為組進行模擬[17]。將樣木分為3 個樹種組：①松類包括馬尾松、黑松P.thunbergii、黃山松P.taiwanensis、濕地松P.elliottii等，其中馬尾松占絕對多數；②杉類包括杉木、柳杉Cryptomeriafortunei、水杉Metasequoia glyptostroboides、池杉Taxodiumascendens、紅豆杉Taxuswallichianavar.chinensis等，其中杉木占絕對多數；③闊葉類包括櫟屬Quercus、栲屬Castanopsis、楓香Liquidambar、木荷、樟科Lauraceae 等植物。

將喬木林樣地按優勢樹種組(樣地中蓄積比例最大的樹種組)統計，喬木林樣木胸徑最大值、均值、中位數和標準差等見表1。其中，最小值均為森林資源連續清查規定的5.0 cm，故沒列出。

表1 喬木林各期樹種組的胸徑統計Table 1 DBH statistics of tree species groups in different stages of arbor forest

1.3 數據預處理

將所有喬木林樣地上的樣木合在一起，統計3 個樹種組的總樣木數和總材積及其比例。材積數據來自森林資源連續清查數據庫，是根據一元材積表計算的。這些數據將用于建立自約束模型。統計各樹種組的進界木和采枯木株數(受松材線蟲病危害的林木在死亡后及時得到清理，所以調查時根據伐根難以區分采伐木和枯損木，所以本研究將兩者合在一起，簡稱采枯木)，用于解釋模型的預測結果。

1.4 自約束模型

設y1、y2和y3分別為松類、杉類和闊葉類樹種組的株數比例或材積比例，t為年份，可得：

式(1)中：λ1、λ2、γ1和γ2為模型參數。由于變量和參數自動滿足0＜y1，y2，y3＜1，y1+y2+y3=1，λ1，λ2＞0，γ1，γ2＜0，所以本研究將式(1)命名為非線性方程組自約束樹種結構模型。式(1)推導過程繁復，本研究略去推導過程的介紹。年份t的起點可以自由設定，本研究為了方便取起點，即1994年時的t取5，1999年時的t取10 等，當然也可以取實際數1994、1999 等。y1，y2，y3存在理論極限值y1∞，y2∞，y3∞(極限的點估計)，根據式(1)可得：

仿真時,將額定工頻50 Hz及采樣值輸入三角函數神經網絡進行訓練,即可快速獲得基波和各次諧波的諧波頻率、幅值及其相位,目標誤差eobj設為0.001。

1.5 極限值的區間估計

極限作為今后的發展趨勢，具有重要意義。設y∞=(y1∞,y2∞,y3∞)T為極限的點估計，各極限點估計的方差D為：

式(3)中： σ2λ1、σ2λ2、σλ1λ2分別為參數λ1、λ2的方差和它們之間的協方差。

式(1)有4 個獨立參數，而用于估計參數的數據只有6 組，自由度只有2，所以只能采用小樣本t分布方法。同時由于比例數不可能小于0、大于1，所以應該采用截尾t分布。設自由度為k的t分布概率密度函數為p(t,k)，則截尾t分布概率密度函數pa(t,k)為：

用式(4) 進行區間估計，具體操作利用Excel 的對應函數實現。設α為顯著性水平，tl(k)α/2、tr(k)α/2滿足以下條件：

則極限真值Y∞的1-α雙側區間估計為：

若tl(k)α、tr(k)α滿足以下條件：

2 結果與分析

2.1 數據預處理結果

2.1.1 各樹種組的樣木株數、材積及其比例 從表2和表3可以看出：株數比例和材積比例表現出了相同的規律，松類比例持續下降，闊葉類比例持續上升，杉類比例呈先上升后下降的變化趨勢。

表2 各樹種組樣木株數及其比例Table 2 Numbers of total sample trees and their proportions by species group

表3 各樹種組樣木材積及其比例Table 3 Volumes of total sample trees and their proportions by species group

2.1.2 進界木株數、采枯木株數和平均單株材積變化 從表4可見：松類和杉類的進界木株數在下降，闊葉類的進界木株數在上升，且總量遠大于松類和杉類之和。5 a 進界率和5 a 采枯率分別為監測周期內發生的進界木總株數和采枯木總株數與初期的活立木總株數的比值。從表5可見：松類平均單株材積單調上升，杉類早期有波動，后期也是明顯上升，闊葉類則是先下降后上升，樹種組的總平均值也整體呈上升趨勢。

表4 各樹種組的樣地進界木和采枯木株數變化Table 4 Numbers of ingrowth trees and numbers of felled and dead trees by species groups

表5 各樹種組的平均單株材積變化Table 5 Change of mean single tree volumes by species groups

2.2 建模及極限點估計結果

基于表2和表3數據，利用Stata 15.1 擬合自約束模型，結果見表6、圖1和圖2。從圖1和圖2可知2019年之后40 a 株數和材積的預測情況，株數比例模型和材積比例模型的σλ1λ2分別為0.239 4 和0.051 9。

圖1 各樹種組的株數比例模型預測圖Figure 1 Graph of predicted percentages of the numbers of trees by species group

圖2 各樹種組的材積比例模型預測圖Figure 2 Graph of predicted percentages of the volumes of trees by species group

表6 自約束模型的擬合結果Table 6 Fitting results of self-constrained models

2.3 極限區間估計結果

從表7可見：區間估計的可靠性均為95%，即α統一取0.05。各t值均為自由度為2 的標準化截尾t分布取值。雙側分布區間估計給出95%分布的下限和上限，例如松類極限值在1.50%～7.63%的概率為95%。單側分布區間估計下限表示極限值大于這個下限的概率為95%，例如松類株數比例的極限值大于2.10%的概率為95%。同理單側分布上限表示極限值小于這個上限的概率為95%，例如松類株數比例的極限值小于6.50%的概率為95%。

表7 區間估計結果Table 7 Results of interval estimation

3 討論

3.1 自約束模型的模擬和預測效果

本研究表明：松類的株數比例和材積比例均在持續單調下降，杉類均是先上升后下降，闊葉類均在快速單調增長。自約束模型很好地擬合了1994—2019年25 a 間臺州市喬木林3 大樹種組的株數比例和材積比例的變化情況。

松類、杉類、闊葉類株數比例的極限點估計分別為4.09%、7.02%、88.89%，材積比例為10.07%、11.94%、77.99%。按95% 的可靠性，松類、杉類、闊葉類株數比例極限分布雙側區間估計分別為1.50%～7.63%、2.67%～12.86%、79.83%～95.71%，材積比例分別為3.64%～18.90%、4.35%～22.31%、59.90%～91.76%。單側區間估計表明：松類、杉類、闊葉類株數比例極限分別大于2.10%、3.71%、82.72%，分別小于6.50%、10.99%、94.07%，材積比例分別大于5.13%、6.12%、65.52%，分別小于16.09%、19.01%、88.51%的可能性為95%。可見，闊葉樹將占絕對優勢，松類的比例將會是最低的，但永遠不會消失，將以低比例一直存在，這與已有研究結論相似[18]，杉類的比例比松類略高，也將以低比例一直存在。

本研究還表明：松類的單株材積在單調上升，尤其近些年上升速度更快。杉類整體也在上升，且也是近年上升更快。闊葉類是先下降后上升，但上升速度比松類和杉類低，尤其和松類相比低更多。闊葉類與松類單株材積的比值，從1994年的0.901 下降到2019年的0.425，闊葉類與杉木的比值從1994年的1.295 下降到2019年的0.626，中間有波動。單株材積的變化規律導致了松類和杉類的株數極限比例分別小于各自的材積比例，而闊葉類正好相反，株數比例大于材積比例。

平均材積的變化差異主要原因在于進界木株數的差異。本研究表明：從1999—2004年這個復查期開始，松類的進界木株數持續大幅減少，進界率也在下降，造成松類小徑階林木的占比減少，導致松類的平均單株材積上升，杉類也有類似的情況，而闊葉類進界木持續增加，導致其平均單株材積增長緩慢。

3.2 針減闊增的主要原因

本研究表明：以松類為優勢的樣地數從2004年的105 個下降到2019年的41 個，減少了60%；株數占比和材積占比分別從1994年的0.729 和0.785 單調下降至2019年的0.140 和0.246，下降速度很快。松類資源下降的主要原因是松材線蟲病的危害[19]。2019年，浙江省新增松材線蟲病發生面積為29.88 萬hm2，新增病枯死松樹達413.74 萬株[20]。從本研究可以看出，松類的采枯率高出闊葉類很多。隨著疫情的加重和疫區的擴大，正常的商業采伐和撫育采伐比例越來越小，而因疫情死亡的比例越來越高。因為松類樹種易受疫情危害，所以松類人工造林急劇下降。由于林分整體密度在增加，整體競爭逐步加劇，馬尾松等陽性樹種在正常的自然演替中處于相對劣勢地位，這也導致了隨著演替的進程馬尾松枯死現象的增多[21-24]。

杉類樹種的株數比例及材積比例在1994 至2004年間增加，其后緩慢減少。杉木優勢的樣地數由1994年的22 個增加至2004年的54 個，而同時期的松類樣地數變化并不大，說明這個時期臺州市造林以杉類為首選。此后樣地數略有下降。杉類中杉木占絕大部分，近年由于經濟效益的下降，杉木造林不多。但杉木具有萌芽的特性，樣地數減少的速度比松類低，進界率近年也比松類略高，但比闊葉類低很多。杉類由于不受松材線蟲病的危害，采枯率比松類低很多，但略高于闊葉類。

闊葉類占優勢的樣地數從1994年的7 個增加到2019年的95 個，株數比例和材積比例分別從1994年的0.052 和0.051 上升到2019年的0.612 和0.457。進界木株數絕對值從1994年的833 株上升到2019年的4 946 株，闊葉類進界率是松類的4.8 倍以上，對比杉類，除了1994—1999 是1.6 倍外，其余是3.3 倍以上。而最近3 個監測期的采枯率低于松類的2/5。闊葉林資源的增加主要有以下幾個方面：①為了公益林建設，大面積的封山育林和林種結構調整，使得闊葉幼林和部分灌木林發展成為闊葉林；②隨著近自然經營的興起，出現了大量的“針改闊”，通過針葉樹采伐或死亡后補植闊葉樹種或促進闊葉樹幼樹生長的形式調整林分的樹種組成[25]；③隨著人們對森林生態系統認識的提高，以及對病蟲害的考慮，有意識地增加了造林時對闊葉樹種的選擇[26]；④松材線蟲病危害導致了松林資源的損失，但同時為闊葉林資源的發展提供了空間；⑤隨著演替的不斷進行，密度增大，馬尾松等陽性先鋒樹種的競爭力下降，并逐漸被淘汰，闊葉樹種的優勢逐漸擴大，隨之達到樹種結構的平衡并趨于穩定[21-24]。

一般來說常綠闊葉林是中亞熱帶的頂級森林群落，具有最好的生態功能，但適量的松類和杉類樹種有利于生物多樣性的保持和生態功能的優化。在控制松材線蟲病的前提下，松類是很好的景觀樹種和用材樹種，杉類也是很好的用材樹種，尤其是大徑材培育的優先選擇樹種。本研究表明：松類和杉類的比例雖然在下降，但始終會以一定的比例存在而不會消亡。

雖然3 個樹種組的變化各異，但3 個樹種組的總株數和總材積均在穩定增長。喬木林的總樣地數從2004年開始基本穩定，但總株數和總材積卻在持續增長，說明林分的相對密度在持續增大。需引起重視的是闊葉類采枯木在快速增加，這其中采伐木的比例不大，主要是由于密度過大導致枯死增加，這說明需要適當加強闊葉林的撫育。

3.3 方法的合理性與不足

本研究是一個區域的喬木林樹種結構變化規律，需要大量的數據支撐。由于2 個原因沒有將數據分為建模和檢驗2 個部分：一是喬木林樣地不足200 個，如果分成2 個部分則會造成建模和檢驗的數據都不夠充分，得不到可靠的結果；二是6 期數據的喬木林樣地并不是固定不變的，前3 期每期樣地都有1/3 移動，后面雖然不再移動，但樣地的地類有變化，有的從喬木林變為其他地類，有的從其他地類變成喬木林，喬木林樣地不固定，喬木林總樣地數也在變化，這給數據建模和檢驗造成了麻煩。所以本研究采用基于建模結果進行可靠性理論分析的方法。

式(1)具有4 個獨立參數，而建模數據只有6 組，自由度只有2，所以本研究采用自由度為2 的小樣本t分布方法進行極限分布的區間估計。由于比例數據不可能小于0、大于1，所以又進行了截尾調整，即采用了截尾t分布模型。松類的株數比例極限分布95%的可能性在1.50%～7.63%，大于2.10%的可能性為95%，是具有實際參考價值的。

4 結論

非線性方程組自約束樹種結構模型適用于大區域樹種組結構的模擬和預測。它反映區域的整體變化情況，但不能反映局部細節。25 a 間臺州市喬木林資源總量在持續增加，但呈現“針減闊增”的現象：①受松材線蟲病危害以及自然演替中的不利因素影響松類比例持續降低；②受市場因素影響，杉類的比例在1999—2004年間達到最大值后開始下降，但幅度比松類小；③闊葉類的比例呈現持續、穩定增加趨勢，將占絕對優勢。松類和杉類的比例雖然在下降，但最終會穩定在一定的比例而不會消亡。這表明臺州市森林生態系統正朝著以闊葉樹種為主、結構更合理的植被群落方向正向演替。