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

張毅 顧鳳岐

(東北林業大學，哈爾濱，150040)

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

1 研究區概況

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

2 材料與方法

2.1 試樣解析

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

圖1 樹高生長與樹齡關系的數據模擬

2.2 模型構建

2.2.1 傳統樹高理論生長模型

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

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

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

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

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

Richards生長方程的一般形式為：y=A(1-e-rt)b，A、r、b>0。A為樹木生長的最大值參數，A=ymax；r為生長速率參數；b為同化作用冪指數m有關的參數，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為零均值平穩過程；εt為白噪聲過程。第一個模型為確定性趨勢模型，第二個模型為自回歸(AR)模型。

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

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

3 結果與分析

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

3.1 傳統理論生長方程擬合模型

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

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

3.2 殘差自回歸模型

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

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

圖2 樹齡-樹高殘差自相關圖和偏自相關圖

表4 樹齡-樹高殘差自相關模型結果

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

圖3 樹齡-樹高殘差Q統計量檢驗

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

3.3 模型擬合

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

表5 2種模型擬合結果

圖4 樹齡-樹高模型擬合結果

3.4 模型預測

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

表6 2種模型預測結果

4 結論

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