北京山區(qū)油松生長(zhǎng)模型研究

北京山區(qū)油松生長(zhǎng)模型研究

越來越多的森林工人關(guān)注林分生長(zhǎng)和收獲模型的兼容性。Cluttter首先提出林分的生長(zhǎng)量與收獲量應(yīng)該是一致的。Sullivan等和Pienaar等利用系統(tǒng)模型解決林分生長(zhǎng)與收獲模型的兼容性問題。在他們所建的模型中,利用限制性條件保證了林分生長(zhǎng)預(yù)測(cè)從A1到A3與從A1到A2、然后再從A2到A3預(yù)測(cè)的結(jié)果一樣,即階段無偏性(step-invariance)。然而正是由于這些限制性條件,降低了系統(tǒng)模型的預(yù)測(cè)精度。Ochi等提出利用林分的年生長(zhǎng)預(yù)測(cè)來克服這些限制性條件,雖然直觀上來看不是對(duì)林分生長(zhǎng)與收獲模型兼容性的研究,但是保證了林分生長(zhǎng)預(yù)測(cè)的階段無偏性。據(jù)查閱相關(guān)文獻(xiàn),國內(nèi)對(duì)階段無偏性這方面的研究很少。而且,目前森林資源監(jiān)測(cè)的發(fā)展趨勢(shì)正向監(jiān)測(cè)周期的年度化發(fā)展,連續(xù)定期清查的數(shù)據(jù)很難準(zhǔn)確的提供森林每年的生長(zhǎng)量和枯損量以及森林資源消長(zhǎng)動(dòng)態(tài)。為此,筆者利用北京山區(qū)油松Pinustabulaeformis的定期調(diào)查數(shù)據(jù),建立油松林分的年生長(zhǎng)預(yù)測(cè)模型,進(jìn)而達(dá)到以下兩個(gè)目的:一、為我國林分生長(zhǎng)預(yù)測(cè)階段無偏性的研究提供一種方法;二、提供林分每年的生長(zhǎng)變化情況,為資源監(jiān)測(cè)周期年度化的發(fā)展提供重要的理論依據(jù)。1樣地的復(fù)測(cè)計(jì)算本研究采用的數(shù)據(jù)來源于北京市林業(yè)調(diào)查設(shè)計(jì)院,其中油松的固定樣地有63個(gè),每個(gè)樣地面積為0.067hm2。樣地主要調(diào)查因子有:林木胸徑、方位角、林分年齡、林分優(yōu)勢(shì)平均高、郁閉度、水平距、坡向、坡位、坡度、海拔高度、土層厚度等因子。樣地每隔5年復(fù)測(cè)1次,本研究利用的數(shù)據(jù)是1986、1991、1996、2001年復(fù)測(cè)的一類清查數(shù)據(jù),利用這些復(fù)測(cè)數(shù)據(jù)可以組成間隔調(diào)查期5年的樣地有156塊,間隔10年的樣地有93塊,間隔15年的樣地有39塊的數(shù)據(jù)集。因此,建模樣地有110塊(間隔5年),檢驗(yàn)樣地共有178塊(間隔5年的樣地有46塊,間隔10年的樣地有93塊,間隔15年的樣地有39塊)。根據(jù)Ochi等和Yue等在研究中的檢驗(yàn)方法,本研究對(duì)3個(gè)間隔期進(jìn)行分開檢驗(yàn)。油松林分樣地分布情況及變量因子統(tǒng)計(jì)分別列于表1、表2。2林分健全at/at全林分生長(zhǎng)方程主要包括林分優(yōu)勢(shì)高、林分株樹密度、林分平均胸徑、林分?jǐn)嗝娣e、林分蓄積等方程。本研究利用林分平方平均胸徑代替林分算術(shù)平均胸徑,因?yàn)樵诹址纸y(tǒng)計(jì)中常用林分平方平均胸徑,而不是算術(shù)平均胸徑。年生長(zhǎng)模型在單木生長(zhǎng)模型中很常見,研究者們通過利用定期調(diào)查數(shù)據(jù),來估計(jì)年生長(zhǎng)模型的參數(shù)。Cao在對(duì)單木生長(zhǎng)模型進(jìn)行年生長(zhǎng)預(yù)測(cè)時(shí),提出了可變生長(zhǎng)率法,該方法考慮了單木每年生長(zhǎng)的變化,而不是按每年固定生長(zhǎng)率進(jìn)行的。張雄清等利用可變生長(zhǎng)率法對(duì)北京地區(qū)油松一類清查數(shù)據(jù)作了研究。那么,林分每年的生長(zhǎng)量也是有變化的,不應(yīng)該固定不變。因此,利用可變生長(zhǎng)率法,對(duì)林分優(yōu)勢(shì)高方程、林分株樹密度方程、林分平均胸徑方程、林分?jǐn)嗝娣e方程和林分蓄積方程的推導(dǎo)如下。(t+1)年時(shí):Hi+1=exp[(At/At+1)ln(Ht)+(1-At/At+1)(α1+α2Ht)];(1)Hi+1=exp{[(At/At+1)ln(Nt)+(1-At/At+1)(β1+β2At+β3/ln(Nt)]};(2)Dt+1=exp{(At/At+1)ln(Dt)+(1-At/At+1)[χ1+χ2At/ln(Nt)+χ3Dt]};(3)Bt+1=Bt+exp{(At/At+1)ln(Bt)+(1-At/At+1)[δ1+δ2Ht+δ3At/ln(Nt)+δ4ln(Bt)]};(4)Mt+1=Mt+exp{(At/At+1)ln(Mt)+(1-At/At+1)[γ1+γ2Ht+γ3/ln(Nt)+γ4ln(Bt)]}。(5)(t+q)年時(shí):Ht+q=exp[(At+q-1/At+q)ln(Ht+q-1)+(1-At+q-1/At+q)(α1+α2Ht+q-1)]};(6)Nt+q=exp{(At+q-1/At+q)ln(Nt+q-1)+(1-At+q-1/At+q)(β1+β2At+q-1+β3/ln(Nt+q-1)]};(7)Dt+q=exp{(At+q-1/At+q)ln(Dt+q-1)+(1-At+q-1/At+q)[χ1+χ2At+q-1/ln(Nt+q-1)+χ3Di+q-1]};(8)Bt+q=Bt+q-1+exp{(At+q-1/At+q)ln(Bt+q-1)+(1-At+q-1/At+q)[δ1+δ2Ht+q-1+δ3At+q-1/ln(Nt+q-1)+δ4ln(Bt+q-1]};(9)Mt+q=Mt+q-1+exp{(At+q-1/At+q)ln(Mt+q-1)+(1-At+q-1/At+q)[γ1+γ2Ht+q-1+γ3/ln(Nt+q-1)+γ4ln(Bt+q-1]}。(10)式(1)～(10)中:q為生長(zhǎng)期(調(diào)查間隔期,5a);At為t年時(shí)林分的平均年齡(a);Ht為t年時(shí)林分的優(yōu)勢(shì)木平均高(m);Nt為t年時(shí)林分的公頃株樹(株/hm2);Dt為t年時(shí)林分平方平均胸徑(cm);Bt為t年時(shí)林分?jǐn)嗝娣e(m2/hm2);Mt為t年時(shí)林分蓄積量(m3/hm2);α1、α2,β1、β3,χ1～χ3,δ1～δ5,γ1～γ4為待估參數(shù)。3參數(shù)估計(jì)與模型評(píng)價(jià)3.1林分平均胸徑和密度相關(guān)性全林整體模型是一個(gè)大的模型體系,該系統(tǒng)中的各模型從不同的側(cè)面對(duì)生長(zhǎng)過程進(jìn)行了描述,這些模型之間是相互關(guān)聯(lián)的。如:林分蓄積與林分平均高生長(zhǎng)、林分平均胸徑和林分密度息息相關(guān)(cov(εi,εj)≠0)。似乎不相關(guān)聯(lián)立方程組作為聯(lián)立方程組的一個(gè)特例,從系統(tǒng)方程上看,這些方程之間不同的因變量依賴的自變量不同或者部分相同,利用最小二乘法對(duì)這些模型進(jìn)行獨(dú)立估計(jì)時(shí),會(huì)產(chǎn)生系統(tǒng)偏差,而采用似乎不相關(guān)聯(lián)立估計(jì)則會(huì)解決系統(tǒng)偏差,保證參數(shù)估計(jì)的有效性。因此,本研究利用似乎不相關(guān)聯(lián)立估計(jì)全林分模型參數(shù),減少系統(tǒng)誤差所帶來的估計(jì)偏差,這樣就能夠提高參數(shù)估計(jì)的有效性和一致性。3.2立地因子及預(yù)測(cè)方法林分優(yōu)勢(shì)高模型、林分株樹密度模型、林分平均胸徑生長(zhǎng)模型、林分?jǐn)嗝娣e生長(zhǎng)模型、林分蓄積生長(zhǎng)模型等可以通過統(tǒng)計(jì)量平均偏差(EMD)、平均絕對(duì)偏差(EMAD)、均方根誤差(ERMSE)和決定系數(shù)(R2)進(jìn)行評(píng)價(jià)。它們的數(shù)學(xué)表達(dá)式分別為:EMD=∑(yi-y′i)/n。(11)EMAD=∑|yi-y′i|/n。(12)ERΜSE=√∑(yi-y′i)2n-p。(13)R2=1-∑(yi-y′i)2/∑(yi-ˉyi)2。(14)式(11)～(14)中:yi為實(shí)際值(林分優(yōu)勢(shì)高、林分株數(shù)密度、林分平均胸徑、林分?jǐn)嗝娣e和林分蓄積量等);y′i、ˉyi分別為它們的預(yù)測(cè)值和平均值;n為樣地個(gè)數(shù);p為參數(shù)個(gè)數(shù)。評(píng)價(jià)一個(gè)模型或者一種預(yù)測(cè)方法的優(yōu)劣,可以利用上述這幾個(gè)統(tǒng)計(jì)量來完成。平均偏差、平均絕對(duì)偏差和均方根誤差小,且決定系數(shù)大,則該模型或該預(yù)測(cè)方法為優(yōu)。本研究參數(shù)估計(jì)和模型評(píng)價(jià)均利用SAS軟件來完成。4林分樹密度模型的檢驗(yàn)表3列出了用林分優(yōu)勢(shì)高模型、林分株樹密度模型、林分平均胸徑生長(zhǎng)模型、林分?jǐn)嗝娣e生長(zhǎng)模型和林分蓄積生長(zhǎng)模型進(jìn)行獨(dú)立估計(jì)時(shí),生長(zhǎng)模型殘差的相關(guān)系數(shù)矩陣。根據(jù)表中數(shù)據(jù)可知,其中4個(gè)生長(zhǎng)模型的殘差在置信水平為0.01下相關(guān)性顯著,這也表明了林分平均胸徑、林分?jǐn)嗝娣e、林分蓄積和林分株樹密度這4個(gè)林分變量因子存在著一定的相關(guān)性,即:cov(εi,εj)≠0。因此,根據(jù)Borders的方法對(duì)這4個(gè)方程的參數(shù)進(jìn)行似乎不相關(guān)聯(lián)立估計(jì),這樣就能夠提高參數(shù)估計(jì)的有效性和一致性。各林分生長(zhǎng)模型的參數(shù)估計(jì)及標(biāo)準(zhǔn)誤見表4。由表4可知,各參數(shù)估計(jì)值比較穩(wěn)定,估計(jì)有效。表5列出了建模樣地及檢驗(yàn)樣地的模型統(tǒng)計(jì)量。從表5可知:在建模數(shù)據(jù)中,林分優(yōu)勢(shì)高模型的EMD為0.0977,EMAD為0.6099,ERMSE為0.8983,R2為0.8690;林分株樹密度模型的EMD為11.1066,EMAD為113.3270,ERMSE為185.0350,R2為0.8578;林分平均胸徑生長(zhǎng)模型的EMD為0.0812,EMAD為0.4920,ERMSE為0.6577,R2為0.9355;林分?jǐn)嗝娣e生長(zhǎng)模型的EMD為0.2836,EMAD為1.15485,ERMSE為1.6296,R2為0.9321;林分蓄積生長(zhǎng)模型的EMD為1.2365,EMAD為4.5492,ERMSE為6.5330,R2為0.9397。經(jīng)過Kolmogorov-Smirnov檢驗(yàn),林分優(yōu)勢(shì)高模型的殘差、林分株樹密度模型的殘差、林分平均胸徑生長(zhǎng)模型的殘差、林分?jǐn)嗝娣e生長(zhǎng)模型的殘差和林分蓄積生長(zhǎng)模型的殘差都服從正態(tài)分布,且分別對(duì)林分優(yōu)勢(shì)高預(yù)測(cè)值、林分株樹密度預(yù)測(cè)值、林分平均胸徑預(yù)測(cè)值、林分?jǐn)嗝娣e預(yù)測(cè)值和林分蓄積預(yù)測(cè)值都沒有明顯的估計(jì)偏差(見圖1)。表6列出了利用3類檢驗(yàn)數(shù)據(jù)所得的各林分生長(zhǎng)模型的評(píng)價(jià)統(tǒng)計(jì)量。q=5時(shí),林分株樹密度模型的EMD為36.0396,EMAD為39.2258,ERMSE為60.3986,R2為0.9810;q=10時(shí),林分株樹密度模型的EMD為52.1159,EMAD為197.1320,ERMSE為284.7610,R2為0.6492;q=15時(shí),林分株樹密度模型的EMD為96.0569,EMAD為238.4860,ERMSE為330.0300,R2為0.4891。比較這3個(gè)間隔期林分株樹密度模型EMD、EMAD、ERMSE和R2值的大小,可知,EMD(36.0396,q=5)<EMD(52.1159,q=10)<EMD(96.0569,q=15),EMAD(39.2258,q=5)<EMAD(197.1320,q=10)<EMAD(238.4820,q=15),ERMSE(60.3986,q=5)<ERMSE(284.7610,q=10)<ERMSE(330.0300,q=15),R2(0.9810,q=5)>R2(0.6492,q=10)>R2(0.4891,q=15)。很明顯,隨著預(yù)測(cè)時(shí)間變長(zhǎng),預(yù)測(cè)偏差越大。而且,在本研究中,主要是針對(duì)林分的年生長(zhǎng)預(yù)測(cè),建立林分年生長(zhǎng)預(yù)測(cè)模型,因此,預(yù)測(cè)間隔的時(shí)間越長(zhǎng),預(yù)測(cè)的偏差會(huì)更大。同樣,林分平均胸徑生長(zhǎng)模型、林分?jǐn)嗝娣e生長(zhǎng)模型和林分蓄積量生長(zhǎng)模型也一樣。5標(biāo)準(zhǔn)ce—結(jié)論林分年