臺州市喬木林樹種結(jié)構(gòu)變化模擬與預測

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

（1.三門縣自然資源和規(guī)劃局，浙江 三門 317100；2.臺州市黃巖區(qū)林業(yè)局，浙江 臺州 318020；3.浙江農(nóng)林大學 環(huán)境與資源學院，浙江 杭州 311300；4.浙江省森林資源監(jiān)測中心，浙江 杭州 310020）

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

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

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

1 研究地區(qū)與方法

1.1 研究區(qū)概況

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

1.2 研究數(shù)據(jù)

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

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

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

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

1.3 數(shù)據(jù)預處理

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

1.4 自約束模型

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

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

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

1.5 極限值的區(qū)間估計

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

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

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

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

則極限真值Y∞的1-α雙側(cè)區(qū)間估計為：

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

2 結(jié)果與分析

2.1 數(shù)據(jù)預處理結(jié)果

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

表2 各樹種組樣木株數(shù)及其比例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 進界木株數(shù)、采枯木株數(shù)和平均單株材積變化 從表4可見：松類和杉類的進界木株數(shù)在下降，闊葉類的進界木株數(shù)在上升，且總量遠大于松類和杉類之和。5 a 進界率和5 a 采枯率分別為監(jiān)測周期內(nèi)發(fā)生的進界木總株數(shù)和采枯木總株數(shù)與初期的活立木總株數(shù)的比值。從表5可見：松類平均單株材積單調(diào)上升，杉類早期有波動，后期也是明顯上升，闊葉類則是先下降后上升，樹種組的總平均值也整體呈上升趨勢。

表4 各樹種組的樣地進界木和采枯木株數(shù)變化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 建模及極限點估計結(jié)果

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

圖1 各樹種組的株數(shù)比例模型預測圖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 自約束模型的擬合結(jié)果Table 6 Fitting results of self-constrained models

2.3 極限區(qū)間估計結(jié)果

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

表7 區(qū)間估計結(jié)果Table 7 Results of interval estimation

3 討論

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

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

松類、杉類、闊葉類株數(shù)比例的極限點估計分別為4.09%、7.02%、88.89%，材積比例為10.07%、11.94%、77.99%。按95% 的可靠性，松類、杉類、闊葉類株數(shù)比例極限分布雙側(cè)區(qū)間估計分別為1.50%～7.63%、2.67%～12.86%、79.83%～95.71%，材積比例分別為3.64%～18.90%、4.35%～22.31%、59.90%～91.76%。單側(cè)區(qū)間估計表明：松類、杉類、闊葉類株數(shù)比例極限分別大于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%?？梢?，闊葉樹將占絕對優(yōu)勢，松類的比例將會是最低的，但永遠不會消失，將以低比例一直存在，這與已有研究結(jié)論相似[18]，杉類的比例比松類略高，也將以低比例一直存在。

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

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

3.2 針減闊增的主要原因

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

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

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

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

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

3.3 方法的合理性與不足

本研究是一個區(qū)域的喬木林樹種結(jié)構(gòu)變化規(guī)律，需要大量的數(shù)據(jù)支撐。由于2 個原因沒有將數(shù)據(jù)分為建模和檢驗2 個部分：一是喬木林樣地不足200 個，如果分成2 個部分則會造成建模和檢驗的數(shù)據(jù)都不夠充分，得不到可靠的結(jié)果；二是6 期數(shù)據(jù)的喬木林樣地并不是固定不變的，前3 期每期樣地都有1/3 移動，后面雖然不再移動，但樣地的地類有變化，有的從喬木林變?yōu)槠渌仡悾械膹钠渌仡愖兂蓡棠玖?，喬木林樣地不固定，喬木林總樣地?shù)也在變化，這給數(shù)據(jù)建模和檢驗造成了麻煩。所以本研究采用基于建模結(jié)果進行可靠性理論分析的方法。

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

4 結(jié)論

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