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

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

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

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

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

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

1.1 研究區(qū)概況

臺州市(28°01′～29°20′N，120°17′～121°56′E)位于浙江省中部沿海，總面積約10 050 km2。該區(qū)屬中亞熱帶季風(fēng)氣候，年均日照時數(shù)為1 800～2 037 h，全年降水日數(shù)為132～171 d，年均降水量為1 632 mm。夏、秋季受臺風(fēng)影響比較大。地帶性植被為中亞熱帶常綠闊葉林，主要的闊葉喬木樹種有殼斗科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 的樣地進(jìn)行了移動。移動部分樣地的目的是為了檢驗固定樣地是否受到特殊對待，即樣地內(nèi)部和外部是否具有相同的經(jīng)營措施，結(jié)果表明沒有特殊對待的情況發(fā)生，即樣地固定后對總體的代表性不會改變。固定樣地為800 m2的正方形。森林資源連續(xù)清查體系是基于概率抽樣理論的機(jī)械抽樣，每5 a 復(fù)測1 次，是森林資源調(diào)查中最具權(quán)威性的數(shù)據(jù)，對于森林的動態(tài)變化分析及預(yù)測十分有利[15]，同時也便于更有效地傳播和利用現(xiàn)有數(shù)據(jù)進(jìn)行可持續(xù)的森林生態(tài)系統(tǒng)經(jīng)營和管理[16]。本研究的對象是喬木樹種，所以只選擇其中的喬木林樣地，根據(jù)地類選擇樣地不改變其隨機(jī)性和代表性。由于經(jīng)濟(jì)林比較特殊，所以本研究不包括喬木經(jīng)濟(jì)林。

由于樹種眾多，對每一樹種進(jìn)行模擬是不現(xiàn)實(shí)的，因此需要將具有共同特性的樹種聚集為組進(jì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ù)和標(biāo)準(zhǔn)差等見表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ù)預(yù)處理

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

1.4 自約束模型

設(shè)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)推導(dǎo)過程繁復(fù)，本研究略去推導(dǎo)過程的介紹。年份t的起點(diǎn)可以自由設(shè)定，本研究為了方便取起點(diǎn)，即1994年時的t取5，1999年時的t取10 等，當(dāng)然也可以取實(shí)際數(shù)1994、1999 等。y1，y2，y3存在理論極限值y1∞，y2∞，y3∞(極限的點(diǎn)估計)，根據(jù)式(1)可得：

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

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

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

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

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

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

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

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

2 結(jié)果與分析

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

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

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

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

圖2 各樹種組的材積比例模型預(yù)測圖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 的標(biāo)準(zhǔn)化截尾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 自約束模型的模擬和預(yù)測效果

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

松類、杉類、闊葉類株數(shù)比例的極限點(diǎn)估計分別為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)勢，松類的比例將會是最低的，但永遠(yuǎn)不會消失，將以低比例一直存在，這與已有研究結(jié)論相似[18]，杉類的比例比松類略高，也將以低比例一直存在。

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

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

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，新增病枯死松樹達(dá)413.74 萬株[20]。從本研究可以看出，松類的采枯率高出闊葉類很多。隨著疫情的加重和疫區(qū)的擴(kuò)大，正常的商業(yè)采伐和撫育采伐比例越來越小，而因疫情死亡的比例越來越高。因為松類樹種易受疫情危害，所以松類人工造林急劇下降。由于林分整體密度在增加，整體競爭逐步加劇，馬尾松等陽性樹種在正常的自然演替中處于相對劣勢地位，這也導(dǎo)致了隨著演替的進(jìn)程馬尾松枯死現(xiàn)象的增多[21-24]。

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

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

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é)果進(jìn)行可靠性理論分析的方法。

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

4 結(jié)論

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