殘差自回歸模型在人工林紅松樹高生長規(guī)律預(yù)測中的應(yīng)用1)

張毅 顧鳳岐

(東北林業(yè)大學(xué)，哈爾濱，150040)

樹高是森林調(diào)查中最重要的因子之一，與其它因子相比，樹高測量比較困難，測量誤差通常也較大。林業(yè)上比較常用的測量方法，除了直接測量法[1]之外，還經(jīng)常使用模型預(yù)測法[2]進(jìn)行預(yù)測，與直接測量相比，用模型預(yù)測的方式更為節(jié)省成本?，F(xiàn)在常用的樹高預(yù)測模型多是采用樹高-胸徑方程的混合效應(yīng)模型[3-8]，是以胸徑為自變量的預(yù)測樹高方法。而以樹齡為自變量的樹高預(yù)測方法多使用傳統(tǒng)的理論生長模型[9-12]，所以本研究以我國東北林區(qū)紅松人工林為研究對象，依據(jù)實測的樹高數(shù)據(jù)，將時間序列分析方法中的殘差自回歸模型與傳統(tǒng)樹木理論生長方程進(jìn)行組合，對比組合后的模型和單獨使用傳統(tǒng)樹木理論生長方程后的結(jié)果，以擬合程度和預(yù)測精度作為評判標(biāo)準(zhǔn)，選出合適的樹高生長模型，可以為紅松的生長狀況和運(yùn)營提供參考。

1 研究區(qū)概況

研究地區(qū)位于吉林省通化縣三棚林場。三棚林場(125°17′～126°25′E，41°20′～42°7′N)于1958年建場，地處長白山腹地，是紅松自然分布的原生地，屬于北溫帶大陸性季風(fēng)氣候。冬季嚴(yán)寒漫長，夏季炎熱短促，氣溫分布在-38.7～36.1 ℃，積雪厚度可達(dá)50 cm，降水量可達(dá)1 120 mm，海拔310～1 220 m。森林資源經(jīng)營總面積3 198 hm2，其中林地面積為3 050 hm2，紅松干果林總面積為1 119.6 hm2，森林總蓄積為28.6萬m3。

2 材料與方法

2.1 試樣解析

2018年10月，在通化縣三棚林場中隨機(jī)選取生長正常、無病蟲害、不斷梢的平均木1株，進(jìn)行樹干解析。將解析木伐倒后，按中央斷面積區(qū)分方法將樹干按1 m區(qū)分，在伐根處、每個區(qū)分段及1.3 m處截取解析圓盤，厚度為5 cm。每個圓盤截取完畢，立即裝入塑料袋中，帶回實驗室。將各圓盤刨光后用掃描儀掃描后將圖像輸入計算機(jī)，使用數(shù)字年輪分析系統(tǒng)(WinDENDRO TM V6.5)測量南北東西4個垂直方向上的各圓盤年輪數(shù)和各年輪寬度。然后，采用樹干解析計算程序，計算解析木的樹齡(74 a)、胸徑(35.1 cm)、樹高(22.57 m)、樹干材積(0.956 09 m3)，每個齡階的樹高生長量數(shù)據(jù)模擬見圖1。

圖1 樹高生長與樹齡關(guān)系的數(shù)據(jù)模擬

2.2 模型構(gòu)建

2.2.1 傳統(tǒng)樹高理論生長模型

在樹木生長模型研究中，根據(jù)生物學(xué)特性作出某種假設(shè)，建立關(guān)于樹木總生長曲線的微分方程或微積分方程，求解后并代入其初始條件或邊界條件，從而獲得該微分方程的特解，這類生長方程叫理論方程。理論方程具有邏輯性強(qiáng)、適用性大、參數(shù)可獨立進(jìn)行檢驗、對未來的生長趨勢進(jìn)行預(yù)測等特點。在林業(yè)上，傳統(tǒng)的一些理論生長方程，包括約翰遜-舒馬赫(Johnson-Schumacher)生長方程、邏輯斯蒂(Logistic)生長方程、米切利希(Mitscherlich)生長方程、坎派茲(Gompertz)生長方程、理查德(Richards)生長方程等[13]，這些方程都可以很好地擬合樹木生長趨勢。所以選擇這5種理論生長方程作為擬合樹高數(shù)據(jù)的生長方程。

Johnson-Schumacher生長方程的一般形式為：y=Ae(-k/(t+a))，A、k、a>0。A為樹木生長的最大值參數(shù)，A=ymax；k、a為方程參數(shù)，且k=1/r。

Logistic生長方程的一般形式為：y=A/(1+ce-rt)，A、c、r>0。A為樹木生長的最大值參數(shù)，A=ymax；c為與初始值有關(guān)的參數(shù)；r為潛在的最大生長率。

Mitscherlich生長方程的一般形式為：y=A(1-e-rt)，A、r>0。A為樹木生長的最大值參數(shù)，A=ymax；r為生長速率參數(shù)。

Gompertz生長方程的一般形式為：y=Aexp(-ce-rt)，A、c、r>0。A為樹木生長的最大值參數(shù)，A=ymax；c為與初始值有關(guān)的參數(shù)；r為潛在的最大生長率。

Richards生長方程的一般形式為：y=A(1-e-rt)b，A、r、b>0。A為樹木生長的最大值參數(shù)，A=ymax；r為生長速率參數(shù)；b為同化作用冪指數(shù)m有關(guān)的參數(shù)，b=1/(1-m)。

2.2.2 殘差自回歸模型

殘差自回歸模型的一般形式[14]為：yt=m(t)+ωt、ωt=α1ωt-1+α2ωt-2+…+αpωt-p+εt，t={1，2,…，T}，T為樣本量；ωt為零均值平穩(wěn)過程；εt為白噪聲過程。第一個模型為確定性趨勢模型，第二個模型為自回歸(AR)模型。

殘差自回歸模型一般形式中的確定性趨勢，實際上是數(shù)據(jù)增長的趨勢，而樹木生長一般呈“S”形增長，所以可以使用傳統(tǒng)樹木理論生長模型替代確定性趨勢模型。生長數(shù)據(jù)一般都具備自相關(guān)性和異方差性，當(dāng)對數(shù)據(jù)建模擬合時，這兩種特性會干擾模型的擬合程度，降低預(yù)測精度。而自回歸模型可以將數(shù)據(jù)中存在的自相關(guān)性去除，合適的自回歸模型在去除數(shù)據(jù)自相關(guān)性的同時，也會削弱數(shù)據(jù)異方差性。

殘差自回歸模型由傳統(tǒng)樹木生長模型和自回歸模型組合而成，所以同時具備兩種模型的優(yōu)點：既可以沿用傳統(tǒng)理論生長方程關(guān)于樹木生長的生物學(xué)性質(zhì)，也可以對生長方程中未提取完的信息進(jìn)行二次提取，同時去除數(shù)據(jù)自相關(guān)性和削弱異方差性，從而提高擬合程度和預(yù)測精度。

3 結(jié)果與分析

將74組數(shù)據(jù)分成擬合組和預(yù)測組，擬合組使用前71組數(shù)據(jù)建立模型并進(jìn)行擬合程度比較，預(yù)測組使用后3組數(shù)據(jù)進(jìn)行預(yù)測精度檢驗。

3.1 傳統(tǒng)理論生長方程擬合模型

通過R軟件代入擬合組數(shù)據(jù)，非線性最小二乘估計迭代參數(shù)得到5種理論生長方程參數(shù)的顯著性檢驗(見表1)、5種樹木理論生長方程的最終形式和擬合效果(見表2)。由表2可見，5種理論方程中最適合擬合數(shù)據(jù)的理論方程為Gompertz生長方程：y=23.93exp(-3.671e-0.051 19t)。

表1 5種理論生長方程參數(shù)的顯著性檢驗

3.2 殘差自回歸模型

將Gompertz生長方程與原數(shù)據(jù)產(chǎn)生的殘差使用Eviews軟件建立自回歸模型。①保證殘差序列是平穩(wěn)過程，對殘差序列先進(jìn)行一階差分，再對差分后的序列進(jìn)行單位根檢驗(擴(kuò)張的迪基-福勒檢驗(ADF檢驗))，得到檢驗的t統(tǒng)計量值(見表3)。由表3可見，殘差序列單位根檢驗值小于顯著性水平在0.01下的t統(tǒng)計量值，認(rèn)為序列是平穩(wěn)的。②保證殘差序列存在自相關(guān)性，對坎派茲(Gompertz)生長模型的殘差進(jìn)行殘差Q統(tǒng)計量檢驗，比較不同滯后階數(shù)下的殘差序列自相關(guān)性，發(fā)現(xiàn)滯后1～10階下的Q統(tǒng)計量檢驗顯著，p<2.2×10-16。結(jié)果顯示，Gompertz模型的殘差序列存在自相關(guān)性，可以建立自回歸模型。

從自相關(guān)圖和偏自相關(guān)圖(見圖2)中可看出，殘差序列的自相關(guān)拖尾、偏自相關(guān)拖尾，說明模型應(yīng)建立自回歸移動平均(ARMA)模型；綜合考慮模型的赤池信息準(zhǔn)則(AIC)值和施瓦茲準(zhǔn)則(SC)值最小原則，建立模型(見表4)。

圖2 樹齡-樹高殘差自相關(guān)圖和偏自相關(guān)圖

表4 樹齡-樹高殘差自相關(guān)模型結(jié)果

由表4可見，模型滿足參數(shù)顯著性。自回歸系數(shù)多項式根的倒數(shù)為0.63、0.49和移動平均系數(shù)多項式根的倒數(shù)為0.99、0.99，均滿足小于1，說明模型滿足平穩(wěn)可逆性。數(shù)據(jù)的純隨機(jī)性見圖3，由圖3可見，殘差Q統(tǒng)計量對應(yīng)的概率全部滿足p>0.05，所以模型滿足參數(shù)檢驗要求，自回歸模型選用ARMA(2,(2))。檢驗殘差序列的異方差性，使用懷特(white)檢驗，懷特檢驗的統(tǒng)計量nR2=1.953，對應(yīng)的概率p=0.923 9，認(rèn)為模型最后無異方差。

圖3 樹齡-樹高殘差Q統(tǒng)計量檢驗

在模型建立的過程中，雖然對殘差建立的不是單純的AR模型，但由于所有的ARMA(p,q)過程都可以通過在等式兩端除以移動平均系數(shù)多項式得到AR(∞)過程[6]，因此綜合得到最終的殘差自回歸模型為：

3.3 模型擬合

對上述2種模型擬合的樹高數(shù)據(jù)情況進(jìn)行匯總(見表5)，2種模型擬合后的數(shù)據(jù)見圖4。

表5 2種模型擬合結(jié)果

圖4 樹齡-樹高模型擬合結(jié)果

3.4 模型預(yù)測

分別使用2種模型對預(yù)測組數(shù)據(jù)進(jìn)行預(yù)測(見表6)，經(jīng)計算得到2種模型的預(yù)測均方誤差分別是0.186 4、0.072 1。

表6 2種模型預(yù)測結(jié)果

4 結(jié)論

在對樹高生長數(shù)據(jù)的擬合上，Gompertz生長模型和殘差自回歸模型的擬合優(yōu)度(R2)都很接近1，均方誤差(MSE)都很接近0，但是相比之下殘差自回歸模型的擬合優(yōu)度更高，均方誤差更小，說明殘差自回歸模型的擬合程度更好；在對樹高生長數(shù)據(jù)的預(yù)測上，殘差自回歸模型比Gompertz生長模型的均方誤差更小，說明殘差自回歸模型的預(yù)測精度更高。綜合模型擬合和預(yù)測結(jié)果，殘差自回歸模型比傳統(tǒng)樹木生長模型更適合描述人工林紅松的樹高生長過程。