王志波,季 蒙,李永樂(lè),李銀祥,馬世明,張海東

(1.內(nèi)蒙古自治區(qū)林業(yè)科學(xué)研究院,呼和浩特 010010；2.興和縣蘇木山林場(chǎng),內(nèi)蒙古 興和 013650；3.內(nèi)蒙古烏蘭察布市林業(yè)科學(xué)研究所,內(nèi)蒙古 集寧 012000)

森林立地是指氣候、地貌、土壤及其它生物等對(duì)森林生長(zhǎng)發(fā)育具有重要作用的環(huán)境因子總體[1]?？茖W(xué)地評(píng)價(jià)立地質(zhì)量,準(zhǔn)確掌握環(huán)境因子對(duì)森林生產(chǎn)力的影響程度,對(duì)宜林性和林木生長(zhǎng)量進(jìn)行合理預(yù)測(cè),是開展?fàn)I造林作業(yè)設(shè)計(jì)、生產(chǎn)潛力評(píng)估、森林撫育管理等工作的基礎(chǔ)[2-3]。評(píng)價(jià)立地質(zhì)量的方法有多種,其中,地位指數(shù)法由于簡(jiǎn)便明了,使用最為普遍[4-5],是最有效、最客觀的一種方法[6]。

地位指數(shù)是指林分在基準(zhǔn)林齡時(shí)優(yōu)勢(shì)木平均高度值,地位指數(shù)模型的構(gòu)建主要采用兩種方法,一是固定基準(zhǔn)林齡的靜態(tài)方程,另一種是可變基準(zhǔn)林齡的動(dòng)態(tài)方程[7-8]。靜態(tài)方程在推導(dǎo)優(yōu)勢(shì)木樹高生長(zhǎng)曲線時(shí),需預(yù)先確定基準(zhǔn)林齡,并利用基準(zhǔn)林齡時(shí)的樹高(地位指數(shù))推算其它時(shí)刻的樹高；而動(dòng)態(tài)方程的基準(zhǔn)林齡是可變的,在推導(dǎo)優(yōu)勢(shì)木樹高生長(zhǎng)曲線時(shí),選擇任何一個(gè)林齡都可以得到一個(gè)相同的地位指數(shù)曲線[9]。因此,動(dòng)態(tài)方程構(gòu)建的地位指數(shù)模型比固定基準(zhǔn)林齡的靜態(tài)方程更加精確,更具有優(yōu)勢(shì)。

自1963年Clutter[10]首次將差分生長(zhǎng)模型(projection growth models)應(yīng)用于美國(guó)火炬松生長(zhǎng)與收獲預(yù)估后,差分模型被廣泛應(yīng)用于模擬林分生長(zhǎng)過(guò)程[11-12]。Bailey等[13]于1974年利用代數(shù)差分法(algebraic difference approach,ADA),通過(guò)指定生長(zhǎng)方程中一個(gè)與形狀有關(guān)的參數(shù)為自由參數(shù),構(gòu)建差分方程,并建立了新西蘭輻射松人工林動(dòng)態(tài)地位指數(shù)模型。使用ADA法推導(dǎo)的差分地位指數(shù)模型,由于只設(shè)置一個(gè)自由參數(shù),構(gòu)建的地位指數(shù)曲線簇僅能滿足多水平漸進(jìn)極值或多形性中的一個(gè)特性[14]。針對(duì)這一問(wèn)題,Cieszewski等[15]擴(kuò)展了ADA法,提出廣義代數(shù)差分法(generalized algebraic difference approach,GADA),并進(jìn)行了詳細(xì)的研究和應(yīng)用[16-18]。此法在推導(dǎo)時(shí)設(shè)置多個(gè)自由參數(shù),能夠構(gòu)建具有多條水平漸近線和多形性兩個(gè)特性的位指數(shù)曲線簇[12,19],因而受到了廣泛關(guān)注。

Richards方程[20]是在Von Bertalanffy生長(zhǎng)理論方程的基礎(chǔ)上擴(kuò)展而來(lái),具有很強(qiáng)的靈活性和廣泛的適應(yīng)性,可對(duì)單木或林分的生長(zhǎng)規(guī)律進(jìn)行很好的擬合與預(yù)測(cè),是近代林業(yè)應(yīng)用最為廣泛的一類生長(zhǎng)曲線方程[21-22]。基于Richards生長(zhǎng)方程,童佳鳴等[23]構(gòu)建了方竹枝葉生長(zhǎng)動(dòng)態(tài)模型；董晨等[24]構(gòu)建了杉木胸徑和樹高生長(zhǎng)模型；范闊等[25]運(yùn)用改進(jìn)的Richards方程構(gòu)建了馬尾松人工林地位指數(shù)模型；李斌成等[26]、牛亦龍等[19]分別利用ADA法、GADA法對(duì)基礎(chǔ)理論生長(zhǎng)方程進(jìn)行推導(dǎo),并構(gòu)建杉木、長(zhǎng)白落葉松人工林地位指數(shù)模型,其中,Richards方程推導(dǎo)的差分地位指數(shù)模型擬合效果最佳,能夠更為精確地描述優(yōu)勢(shì)木樹高生長(zhǎng)進(jìn)程。

華北落葉松(Larixprincipis-rupprechtii)為松科落葉松屬的落葉針葉喬木,在內(nèi)蒙古主要分布在陰山、燕山和大興安嶺等地區(qū),是重要的山地造林樹種。該樹種具有耐寒、生長(zhǎng)速度快,材質(zhì)優(yōu)良等特點(diǎn),是重要的用材林樹種,在水土保持、涵養(yǎng)水源、調(diào)節(jié)氣候等方面也發(fā)揮了重要作用。建立華北落葉松人工林地位指數(shù)模型,科學(xué)的評(píng)估林分立地質(zhì)量,對(duì)人工林經(jīng)營(yíng)管理和提高林地生產(chǎn)力等具有重要作用。本研究以陰山山地蘇木山林場(chǎng)華北落葉松人工林為研究對(duì)象,基于Richards生長(zhǎng)方程,分別采用ADA法、GADA法構(gòu)建差分地位指數(shù)模型,并對(duì)各模型擬合情況進(jìn)行全面具體的分析,篩選出適合該林區(qū)生產(chǎn)需求地位指數(shù)模型,用于指導(dǎo)實(shí)踐,為華北落葉松造林、分類經(jīng)營(yíng)及林分質(zhì)量的精準(zhǔn)提升提供參考和依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于內(nèi)蒙古烏蘭察布市興和縣蘇木山林場(chǎng),地處陰山山脈東南段,地理坐標(biāo)為東經(jīng)113°38′45″～114°05′05″,北緯40°26′58″～40°42′23″,屬中溫帶大陸性季風(fēng)半干旱氣候區(qū),年均氣溫4.2℃,年均降水量397mm,年均蒸發(fā)量2 060mm,年均風(fēng)速3.8m/s,年均無(wú)霜期110d。林場(chǎng)現(xiàn)有林地面積2.2萬(wàn)hm2,主要喬木樹種有華北落葉松、白樺(Betulaplatyphylla)、樟子松(Pinussylvestrisvar.mongolica)、油松(Pinustabuliformis)等。華北落葉松人工林面積0.68萬(wàn)hm2,海拔1 500～2 300m。林下灌木主要有繡線菊(Spiraeapubescens)、山刺玫(Rosadavurica)、懸鉤子(Rubussachalinensis)等,林下草本主要有菊科(Compositae)、毛茛科(Ranunculaceae)、禾本科(Gramineae)、莎草科(Cyperaceae)等。

1.2 數(shù)據(jù)來(lái)源

本研究所使用的數(shù)據(jù)為2018—2019年調(diào)查的43株華北落葉松優(yōu)勢(shì)木解析木數(shù)據(jù)。解析木選取時(shí),根據(jù)研究區(qū)華北落葉松人工林分布現(xiàn)狀,分別在不同立地條件、不同林齡、林相完整的林分內(nèi)設(shè)置標(biāo)準(zhǔn)地,規(guī)格為20m×30m。設(shè)置完成后,在記錄紙上畫出標(biāo)準(zhǔn)地縮略圖,記錄地理位置、地形地貌、坡度、坡向、坡位、海拔以及造林年度、林分密度等因子。對(duì)胸徑大于5cm的樹木每木檢尺,測(cè)量胸徑、樹高、枝下高、冠幅等生長(zhǎng)指標(biāo)；根據(jù)測(cè)得的數(shù)據(jù),每塊標(biāo)準(zhǔn)地挑選出1株優(yōu)勢(shì)木進(jìn)行樹干解析。表1為解析木基本信息,將43株解析木數(shù)據(jù)隨機(jī)分為兩部分,32株用于模型擬合,11株作為獨(dú)立數(shù)據(jù)用于模型檢驗(yàn)。

表1 優(yōu)勢(shì)木基本信息

1.3 模型建立

1.3.1基礎(chǔ)模型選擇

本研究選擇Richards方程作為推導(dǎo)差分地位指數(shù)模型的基礎(chǔ)方程。表達(dá)式為:

h=a(1-e-bt)c

(E0)

式中:h為林分優(yōu)勢(shì)木樹高(m),t為林齡(a)；a為生長(zhǎng)極限參數(shù)(或漸近線參數(shù)),表示生長(zhǎng)因子的極值,是一個(gè)時(shí)間序列內(nèi)的生物量上限；b為生長(zhǎng)率參數(shù)(或尺度參數(shù)),影響曲線拐點(diǎn)的位置,與生長(zhǎng)速率密切相關(guān),起著調(diào)節(jié)生長(zhǎng)時(shí)間的作用,此值越大,生長(zhǎng)時(shí)間便越短；c為形狀參數(shù),決定曲線的拐點(diǎn)位置和形狀[19-21,25];e為自然常數(shù)。

1.3.2差分方程推導(dǎo)

1.3.2.1一般代數(shù)差分方程(ADA)

利用ADA法推導(dǎo)差分地位指數(shù)模型通常包括以下幾個(gè)步驟[9]:

1)選擇一個(gè)生長(zhǎng)方程為基礎(chǔ)方程。

2)設(shè)定方程的一個(gè)參數(shù)與不可觀測(cè)的立地變量X有關(guān),即參數(shù)=X。

3)假定初始時(shí)刻值為(t1,h1),將其代入基礎(chǔ)方程中便得到初始方程。

4)從初始方程中解出X的值,再將其解代入基礎(chǔ)方程中去,這樣我們便可以得到一個(gè)差分方程。

這種方法所推導(dǎo)出來(lái)的差分方程根據(jù)所選擇自由參數(shù)不同,可得到單形可變漸近線或多形單漸近線兩種曲線簇。本研究根據(jù)該法推導(dǎo)得到的差分方程如表 2所示。

1.3.2.2廣義代數(shù)差分方程(GADA)

利用GADA法推導(dǎo)差分地位指數(shù)模型通常包括以下步驟[12,19]:

1)同樣選擇一個(gè)生長(zhǎng)方程為基礎(chǔ)方程。

2)指定方程中的兩個(gè)或多個(gè)參數(shù)為自由參數(shù)。

3)提出一個(gè)與立地質(zhì)量有關(guān)的變量X0,并假設(shè)自由參數(shù)與變量成各種數(shù)量關(guān)系(如,線性、反函數(shù)、二次方、指數(shù)等)。

4)將上述函數(shù)關(guān)系代入基礎(chǔ)方程,并解出(t1,h1)時(shí)刻X0的表達(dá)式。其中,t1為指定林齡,h1為指定林齡下的優(yōu)勢(shì)木樹高。當(dāng)t1為基準(zhǔn)林齡時(shí),h1即為地位指數(shù)。

5)將解出的X0表達(dá)式代入基礎(chǔ)方程中,即得到了通過(guò)GADA法推導(dǎo)的差分地位指數(shù)方程。

這種方法構(gòu)建的地位指數(shù)模型可以滿足多條水平漸近線和多形性兩個(gè)屬性。

參考曹元帥等[7]、趙磊等[14]、牛亦龍等[19]廣義代數(shù)差分方程推導(dǎo)過(guò)程中參數(shù)設(shè)定方法,本研究將方程參數(shù)a和c指定為自由參數(shù)進(jìn)行推導(dǎo),結(jié)果如表3所示。

表2 一般代數(shù)差分方程

表3 廣義代數(shù)差分方程

1.4 模型擬合與檢驗(yàn)

使用SAS 9.0軟件的NLIN過(guò)程進(jìn)行模型參數(shù)的估計(jì)和擬合統(tǒng)計(jì)量的計(jì)算。檢驗(yàn)?zāi)Ｐ偷暮脡闹饕罁?jù)兩點(diǎn):一是模型及參數(shù)的生物學(xué)涵義,二是統(tǒng)計(jì)指標(biāo)表征的模型實(shí)際擬合效果。本研究采用回歸分析中常用的決定系數(shù)(R2)、均方根誤差(RMSE)和平均絕對(duì)誤差(MAE)3 個(gè)統(tǒng)計(jì)量指標(biāo)來(lái)評(píng)價(jià)模型的擬合質(zhì)量。利用擬合樣本計(jì)算決定系數(shù)(R2)和均方根誤差(RMSE)；決定系數(shù)越接近于1,均方根誤差越小,說(shuō)明模型越好。利用獨(dú)立檢驗(yàn)樣本數(shù)據(jù)計(jì)算平均絕對(duì)誤差(MAE),數(shù)值越小說(shuō)明模型預(yù)測(cè)效果越好。各指標(biāo)計(jì)算公式如下:

(1)

(2)

(3)

2 結(jié)果與分析

2.1 模型擬合結(jié)果

使用華北落葉松人工林優(yōu)勢(shì)木建模數(shù)據(jù)(表1)對(duì)Richards基礎(chǔ)生長(zhǎng)方程及各差分方程(表2、表3)進(jìn)行擬合,結(jié)果如表4所示。6個(gè)模型所有參數(shù)的估計(jì)值均收斂較好,達(dá)到極顯著水平(P<0.01),全部通過(guò)檢驗(yàn)。模型E0的決定系數(shù)(R2)為0.891 5,均方根誤差(RMSE)為0.772 6,平均絕對(duì)誤差(MAE)為1.735 3；利用ADA法和GADA法推導(dǎo)的差分模型(E1～E5),決定系數(shù)(R2)均在0.96以上,均方根誤差(RMSE)在0.86～0.96之間,平均絕對(duì)誤差(MAE)在0.42以下。擬合結(jié)果表明,無(wú)論是一般代數(shù)差分方程還是廣義代數(shù)差分方程,對(duì)于建模數(shù)據(jù)的擬合效果均優(yōu)于基礎(chǔ)方程,具有更好的預(yù)測(cè)能力,差分方程明顯提高了模型的擬合精度。從擬合統(tǒng)計(jì)指標(biāo)計(jì)算結(jié)果來(lái)看,模型E2最優(yōu),其次為模型E4,再次為模型E3,E5和E1,模型E0最差,初步篩選出模型E2,E4和E3為備選模型進(jìn)行檢驗(yàn),并與基礎(chǔ)模型對(duì)比分析。

表4 華北落葉松人工林差分地位指數(shù)模型的參數(shù)估計(jì)值及相關(guān)統(tǒng)計(jì)指標(biāo)

2.2 模型檢驗(yàn)

模型E2,E4,E3和E0的殘差隨優(yōu)勢(shì)木樹高預(yù)測(cè)值變化情況如圖1所示。4個(gè)模型的殘差都隨機(jī)分布在y=0兩側(cè),模型E2,E4沒(méi)有明顯的變化趨勢(shì),符合等方差性假設(shè)；模型E3對(duì)較高優(yōu)勢(shì)木樹高的預(yù)測(cè)偏大；模型E0殘差隨優(yōu)勢(shì)木樹高預(yù)測(cè)值的變化呈現(xiàn)非線性關(guān)系。模型E0的殘差波動(dòng)范圍在-4～4m之間,對(duì)較低的優(yōu)勢(shì)木樹高擬合效果較好,殘差絕對(duì)值隨著林齡的增大呈現(xiàn)增大趨勢(shì)。差分地位指數(shù)模型殘差波動(dòng)范圍在-1.5～1.5m之間,明顯優(yōu)于基礎(chǔ)模型,殘差絕對(duì)值隨林齡的增大呈現(xiàn)先減小后增大的趨勢(shì)。模型E2,E4在預(yù)測(cè)值10m時(shí),殘差最集中,對(duì)中等的優(yōu)勢(shì)木樹高擬合效果最好,對(duì)較高的優(yōu)勢(shì)木樹高擬合效果相對(duì)較差。

圖1 模型E2,E4,E3和E0殘差圖

圖2 差分地位指數(shù)模型 E2,E4及E3所繪制的地位指數(shù)曲線簇

2.3 模型選擇

理想的地位指數(shù)模型滿足3個(gè)條件,即S形生長(zhǎng)曲線、多形性和多條水平漸近線[7],利用備選的3個(gè)較優(yōu)模型(E2,E4,E3),根據(jù)解析木數(shù)據(jù),樹高生長(zhǎng)量在25a時(shí)開始趨于穩(wěn)定,因此以25a為基準(zhǔn)林齡,地位指數(shù)范圍從7～15m,地位指數(shù)級(jí)距為2m繪制地位指數(shù)曲線簇,結(jié)果如圖2所示。模型E2具備S形生長(zhǎng)曲線和多形性,但趨向于同一水平漸近線,對(duì)低地位指數(shù)(7)的擬合效果不好,曲線形狀較為平直,不符合樹高隨林齡的變化規(guī)律。模型E3具有多形性,但不滿足S 形生長(zhǎng)曲線和多條水平漸近線的要求,對(duì)低地位指數(shù)(7,9)的低林齡段(0～10a)樹高擬合也偏小。模型E4具備S 形生長(zhǎng)曲線、多形性和多條水平漸近線特征,是理想的地位指數(shù)模型,其表達(dá)式為:

h2=eX0(1-e-0.0313t2)(4.6150-1.002X0)

式中:h2為待預(yù)測(cè)林齡t2時(shí)的優(yōu)勢(shì)木樹高；h1是已知實(shí)測(cè)林齡t1時(shí)的樹高,當(dāng)t1為基準(zhǔn)林齡(25a)時(shí),h1即為地位指數(shù)(7,9,11,13,15)。

2.4 優(yōu)勢(shì)木樹高連年生長(zhǎng)量

利用篩選的最優(yōu)模型(E4)計(jì)算不同地位指數(shù)華北落葉松優(yōu)勢(shì)木樹高連年生長(zhǎng)量,結(jié)果如圖3所示。連年生長(zhǎng)量隨著林齡的增長(zhǎng)迅速達(dá)到最大值,然后緩慢下降。地位指數(shù)7,9,11,13,15的連年生長(zhǎng)量達(dá)到最大時(shí)的林齡為19,15,13,10,6a,對(duì)應(yīng)值分別為0.33,0.41,0.49,0.58,0.69m。優(yōu)勢(shì)木樹高連年生長(zhǎng)量達(dá)到最大值時(shí)的點(diǎn)即為地位指數(shù)曲線的拐點(diǎn)。

圖3 不同地位指數(shù)優(yōu)勢(shì)木樹高連年生長(zhǎng)量

3 結(jié)論與討論

本研究以Richards基礎(chǔ)生長(zhǎng)方程為原型,推導(dǎo)了5個(gè)與基準(zhǔn)林齡無(wú)關(guān)的差分地位指數(shù)模型,利用內(nèi)蒙古興和縣蘇木山林場(chǎng)43株解析木數(shù)據(jù)對(duì)模型參數(shù)進(jìn)行擬合,各模型擬合精度均較高。從各模型的參數(shù)估計(jì)值和擬合統(tǒng)計(jì)指標(biāo)計(jì)算結(jié)果來(lái)看,差分模型的擬合效果均優(yōu)于基礎(chǔ)模型,這與李斌成等[26]對(duì)杉木差分地位指數(shù)模型研究,段愛國(guó)等[27]對(duì)杉木多形地位指數(shù)模型研究結(jié)果相一致。ADA法推導(dǎo)的3個(gè)差分地位指數(shù)模型中(E1,E2,E3),以b或c為自由參數(shù)的模型擬合結(jié)果要好于以a為自由參數(shù)的模型,這與相聰偉[9]的研究結(jié)果相同。GADA法推導(dǎo)的2個(gè)差分模型中,假設(shè)自由參數(shù)c與變量X0成線性關(guān)系更符合華北落葉松優(yōu)勢(shì)木樹高生長(zhǎng)規(guī)律。

以25a為基準(zhǔn)林齡,地位指數(shù)級(jí)距為2m,地位指數(shù)范圍為7～15m,繪制地位指數(shù)曲線簇。以b或c為自由參數(shù)的模型(E2,E3)所繪制的地位指數(shù)曲線簇中,優(yōu)勢(shì)木樹高生長(zhǎng)曲線為多形的,能夠很好地區(qū)分開不同立地質(zhì)量?jī)?yōu)勢(shì)木樹高生長(zhǎng)過(guò)程,但是這類曲線簇中的曲線由于模型擁有一個(gè)共同的漸近線參數(shù)a,所有曲線都會(huì)趨向于同一條上漸近線,表現(xiàn)出不同立地質(zhì)量?jī)?yōu)勢(shì)木樹高值最終是相同的,這與華北落葉松實(shí)際生長(zhǎng)過(guò)程中立地質(zhì)量越好其優(yōu)勢(shì)木樹高漸近值就越大相違背。倪成才等[12]也認(rèn)為一般代數(shù)差分法構(gòu)建的差分模型僅指定一個(gè)參數(shù)為自由參數(shù),對(duì)林分生長(zhǎng)過(guò)程差異的解釋能力是有限的。利用最優(yōu)模型(E4)繪制地位指數(shù)曲線簇,并計(jì)算不同地位指數(shù)華北落葉松優(yōu)勢(shì)木樹高連年生長(zhǎng)量,結(jié)果表明,地位指數(shù)越大,優(yōu)勢(shì)木樹高及其連年生長(zhǎng)量的極值也越大,曲線的拐點(diǎn)也越靠近y軸。這說(shuō)明在較好的立地條件下,華北落葉松優(yōu)勢(shì)木樹高生長(zhǎng)量在前期快速增長(zhǎng),并在較短的時(shí)間內(nèi)達(dá)到最大值,到達(dá)拐點(diǎn)位置后,生長(zhǎng)量下降的速度也越快。牛亦龍等[19]認(rèn)為,這是由于立地質(zhì)量較差時(shí),林木生長(zhǎng)的環(huán)境條件是影響優(yōu)勢(shì)木樹高生長(zhǎng)的主要因素,而立地條件較好時(shí),林木自身的生物學(xué)特性是影響優(yōu)勢(shì)木樹高生長(zhǎng)的主要因素。

通過(guò)統(tǒng)計(jì)分析、殘差分析和地位指數(shù)曲線簇比較,篩選出采用GADA法,以自由參數(shù)a=eX0,c=c1+c2X0推導(dǎo)的模型E4為最佳模型,雖然從擬合統(tǒng)計(jì)指標(biāo)計(jì)算結(jié)果上看,它不是最好的,但利用該模型繪制的地位指數(shù)曲線簇符合S 形生長(zhǎng)曲線、多形性和多條水平漸近線的特性,是理想的地位指數(shù)模型。利用GADA 法,合理地選擇自由參數(shù)與變量的數(shù)量關(guān)系,構(gòu)建動(dòng)態(tài)地位指數(shù)模型,可以滿足生物學(xué)和統(tǒng)計(jì)學(xué)兩方面的要求[7],能夠科學(xué)地對(duì)研究區(qū)華北落葉松人工林立地質(zhì)量進(jìn)行預(yù)測(cè)和評(píng)價(jià)。

本研究?jī)H基于Richards方程進(jìn)行差分方程推導(dǎo),而且廣義差分方程也僅設(shè)置了兩種自由參數(shù)的假設(shè),利用其它生長(zhǎng)方程及不同自由參數(shù)數(shù)量關(guān)系的假設(shè),所得到的模型擬合效果如何,還有待進(jìn)一步研究。