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

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

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

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于內(nèi)蒙古烏蘭察布市興和縣蘇木山林場,地處陰山山脈東南段,地理坐標(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,年均無霜期110d。林場現(xiàn)有林地面積2.2萬hm2,主要喬木樹種有華北落葉松、白樺(Betulaplatyphylla)、樟子松(Pinussylvestrisvar.mongolica)、油松(Pinustabuliformis)等。華北落葉松人工林面積0.68萬hm2,海拔1 500～2 300m。林下灌木主要有繡線菊(Spiraeapubescens)、山刺玫(Rosadavurica)、懸鉤子(Rubussachalinensis)等,林下草本主要有菊科(Compositae)、毛茛科(Ranunculaceae)、禾本科(Gramineae)、莎草科(Cyperaceae)等。

1.2 數(shù)據(jù)來源

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

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

1.3 模型建立

1.3.1基礎(chǔ)模型選擇

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

h=a(1-e-bt)c

(E0)

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

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

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

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

1)選擇一個生長方程為基礎(chǔ)方程。

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

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

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

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

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

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

1)同樣選擇一個生長方程為基礎(chǔ)方程。

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

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

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

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

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

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

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

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

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

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

(1)

(2)

(3)

2 結(jié)果與分析

2.1 模型擬合結(jié)果

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

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

2.2 模型檢驗(yàn)

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

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

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

2.3 模型選擇

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

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

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

2.4 優(yōu)勢木樹高連年生長量

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

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

3 結(jié)論與討論

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

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

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

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