基于啞變量的濕地松林分斷面積生長模型

（福建農(nóng)林大學 林學院，福建 福州 350002）

濕地松Pinus elliottii原生于古巴、美國東南部等地[1]，其適應性強，已成為我國南方重要造林樹種之一[2-4]。濕地松木材質(zhì)量較好[5-6]且松脂產(chǎn)量較高，具有不錯的經(jīng)濟效益，但采脂對林木生長有抑制作用，其對胸徑生長影響比對樹高生長影響更明顯[7-8]，胸徑是外業(yè)調(diào)查的主要測樹因子，且為測算林分斷面積的重要因子，其測量簡單并具有較好的穩(wěn)定性與可預測性。斷面積是估算木材材積的重要依據(jù)，建立采脂與未采脂林分兼用的斷面積模型可作為濕地松林分規(guī)劃經(jīng)營模式的參考依據(jù)。斷面積模型是林分生長模型的核心[9]和森林可持續(xù)經(jīng)營的重要工具[10]，國內(nèi)外學者對生長模型進行了許多研究，其類型主要包括單木模型和全林分模型[11]。單木模型研究對象為單株立木不考慮其實際年齡，只需確定調(diào)查的間隔年數(shù)，且調(diào)查間隔期越長產(chǎn)生的誤差越大，該方法適用于預估單株林木的生長狀況，用于預測全林分時誤差會累加導致估測精度降低[12-13]。而全林分模型不適用于預估單木生長，但對于林分水平預測具有較好的效果[14]。

如果考慮采脂與未采脂對濕地松林分斷面積的影響，則對于這兩種不同經(jīng)營措施的林分需分別建立林分斷面積模型，在實際工作中會耗費更多人力物力，且可能會出現(xiàn)模型不相容的問題，采用啞變量模型為兼顧模型擬合效果和減小工作量提供了途徑[15]。朱光玉等[16]在研究櫟類天然林林分斷面積生長模型時引入了立地質(zhì)量等級為啞變量，李忠國等[17]在日本落葉松生長模型研究引入?yún)^(qū)域特征為啞變量，王亞楠等[18]在華山松樹高生長模型研究中引入不同種源為啞變量，Zeng等[19]在生物量的研究引入了森林起源為啞變量，王金池等[20]在云南松林分蓄積量生長模型的研究引入間伐與未間伐指標為啞變量，華偉平等[21]在黃山松地位級指數(shù)模型的研制引入立地質(zhì)量等級為啞變量，賈煒瑋等[22]利用啞變量方法構建考慮不同地域的林分碳儲量模型，華偉平等[23]以黃山松為研究對象，構建了兼容性林分生物量模型及全林分收獲模型。本研究將引入采脂與未采脂林分這兩種經(jīng)營措施作為啞變量，構建濕地松林分斷面積啞變量模型，以達到提升模型精度的目的。

本研究根據(jù)濕地松樣地數(shù)據(jù)，參考前人的研究[16-23]采用了4 類基礎模型，比較所選模型均方差（RMSE）、殘差平方和（SSE）、決定系數(shù)（R2）、模型精度（v）和模型準確度（P）等評價指標，得到最優(yōu)基礎模型，在此基礎上，將采脂林分與未采脂林分作為啞變量引入最優(yōu)基礎方程不同參數(shù)中分別建模，運用含熵權值的TOPSIS 法選出最優(yōu)啞變量模型，為測算采脂與未采脂林分濕地松木材材積與林分經(jīng)營模式提供參考。

1 研究區(qū)概況

研究區(qū)位于福州市，位于福建省東部，閩江下游與沿海地區(qū)。全市陸地總面積約為19 968 km2，其地貌為典型的河口盆地，海拔在600～1 000 m之間，地處25°15′～26°39′N，118°08′～120°31′E。該市為亞熱帶季風氣候，溫暖濕潤，年均溫度20～25℃，無霜期326 d，年均降水量900～2 100 mm，陽光充足，年均日照數(shù)1 700～1 980 h，濕地松是該區(qū)域重要造林樹種之一。

2 研究方法

2.1 數(shù)據(jù)采集和整理

在福州市濕地松人工林中設置臨時樣地217塊，記錄樣地內(nèi)林木平均胸徑、平均樹高、優(yōu)勢木平均高、平均年齡、郁閉度和海拔等因子。其中采脂林分樣地128 塊，未采脂林分樣地89 塊。對217 塊樣地進行隨機抽樣，從中隨機抽取172個樣本作為建模數(shù)據(jù)（其中采脂樣本92 個，未采脂樣本80 個），其余45 個樣本作為檢驗數(shù)據(jù)（其中采脂樣本26 個，未采脂樣本19 個），樣地的基本情況見表1。

表1 樣地基本情況Table 1 The conditions of plots

2.2 林分密度指數(shù)

林分密度指數(shù)是將現(xiàn)實林分株數(shù)換算為標準平均直徑時所具有的單位面積株數(shù)[24]，其計算公式如下：

式（1）中：ISD是林分密度指數(shù)；N是現(xiàn)實林分的公頃株數(shù)；D0是標準平均直徑；D是現(xiàn)實林分平均直徑；b為自然稀疏率；參考前人的研究標準平均直徑取10 cm，自然稀疏率取-1.605[23]。

2.3 模型的選擇

參考前人的研究[16-23]選用了理查德方程、邏輯斯蒂模型、Mitscherlich 模型和Schumacher 模型。所選模型具有較好的適應性及解釋性，且形式和計算較為簡單?；谶@4 個模型引入優(yōu)勢木平均高和年齡因子，選取林分密度指數(shù)作為林分密度指標用于擬合林分斷面積模型，模型形式如下：

式（2）～（5）中：G為林分公頃斷面積；H為林分優(yōu)勢木平均高；t為林分年齡；ISD為林分密度指數(shù)；b1、b2、b3、b4、b5、b6為待求參數(shù)。

2.4 模型評價指標

采用均方差（RMSE）、殘差平方和（SSE）、決定系數(shù)（R2）、模型精度（v）和模型準確度（P）等模型評價指標。其中高優(yōu)指標為R2和P數(shù)值越大越好，低優(yōu)指標為RMSE、SSE 和v數(shù)值越小越好，計算公式如下：

式（10）中：yi為林分實測值；為模型預估值；為模型預估值平均值；為實測值平均值；n為樣本數(shù)量；t0.05為置信水平為95%時t分布值；p為模型中參數(shù)的個數(shù)。

2.5 啞變量的設置

啞變量（虛擬變量）是定性變量，通常取值為0 或1，將定性因子做（0，1）化的展開，變量δ（x，i）表示為：δ（x，i）=0 或1，當x為第i等級時=1，否者=0，因此稱變量δ（x，i）為啞變量。本研究依據(jù)林分不同經(jīng)營類型使用采脂與未采脂作為啞變量，將用采脂與未采脂定性代碼1 或0 表示，第i種林分經(jīng)營類型表示為Si，定性數(shù)據(jù)Si轉(zhuǎn)化為（0，1）形式：Si=1 或0，當x為第i等級時為1，否者為0。式中：i=1,2；S1,S2分別是采脂與未采脂的定性代碼。

2.6 含熵權值的TOPSIS 法

本研究采用含熵權值的TOPSIS法，在計算模型評價指標權重時運用熵值法取代一般的主觀權重法，避免了人為主觀性的影響[25-26]。其計算步驟如下：

1）根據(jù)模型指標求解結果構建m個評價對象、n個評價指標的判斷矩陣，公式如下：

2）將判斷矩陣歸一化處理，得到無量綱化后的矩陣B，求解出的指標值范圍為0 至1，其值越大代表指標越優(yōu)，最優(yōu)解的值為1，最差解值為0。高優(yōu)指標（效益型指標）公式如下：

低優(yōu)指標（成本型指標）公式如下：

3）計算熵權值Wj

3 結果與分析

3.1 斷面積基礎模型擬合結果

使用172 塊樣地的建模數(shù)據(jù)，利用R 軟件的遺傳算法對基礎模型進行參數(shù)求解，并計算相關的評價指標，結果見表2。

從表2可知：4 個基礎模型的R2值均高于0.9，模型精度均大于95%，說明基礎模型的擬合效果較好。模型決定系數(shù)由大到小的排序為模型5＞模型3＞模型2＞模型4，精度由大到小的排序為模型5＞模型3＞模型2＞模型4。其中模型5的R2和P值分別為0.983 31 與98.20%高于其他模型且RMSE、SSE 和v等低優(yōu)指標值為1.92、275.7、98.20%、0.049 05 均低于其他模型，說明Schumacher 模型可能更適合用于模擬福州市濕地松人工林斷面積增長。模型5 的各項評價指標均為最優(yōu)，可直接表明其擬合效果最好，無需進一步的篩選。因此，將模型5 作為最優(yōu)基礎模型用于構建啞變量模型。

表2 模型參數(shù)及評價指標Table 2 Parameters and evaluation indexs of models

3.2 啞變量模型及TOPSIS 結果

依據(jù)基礎模型擬合結果，在模型5 中引入啞變量，經(jīng)過嘗試發(fā)現(xiàn)，在模型不同參數(shù)引入啞變量時模型各項評價指標存在差異，因此在模型的不同參數(shù)中加入啞變量，其模型參數(shù)求解值見表3，評價指標值見表4。

表3 啞變量模型參數(shù)?Table 3 Dumb variable model parameter

表4 啞變量模型評價指標Table 4 Dumb variable model evaluation index

模型擬合效果與預估精度直接相關，選用的5個評價指標從不同方面體現(xiàn)模型的優(yōu)度。在模型5參數(shù)b3中加入啞變量后所得模型的決定系數(shù)最高達到0.998 96，于模型5 參數(shù)b5中引入啞變量后所得到模型精度最高為99.383%，且RMSE、SSE和v等低優(yōu)指標最低，不同指標最優(yōu)值分別來自不同模型，無法直接篩選最優(yōu)啞變量模型。因此，采用含熵權值的TOPSIS 法對不同模型評價指標進行綜合分析，各評價指標相對最優(yōu)解距離見表5，各評價指標熵權值見表6。

表5 各評價指標相對最優(yōu)解距離Table 5 Relative optimal solution distance for each evaluation index

表6 各評價指標熵權值Table 6 Entropy weight of each evaluation index

在模型5 不同參數(shù)中加入啞變量R2值與其最優(yōu)解距離由大到小的排序為b2＞b1＞b4＞b6＞b5＞b3，其中b3參數(shù)加入啞變量時R2值與最優(yōu)解距離最小為0，表明其為R2最優(yōu)解，b2參數(shù)加入啞變量時其R2值與其最優(yōu)解距離最大為0.036 100 9，表明其為R2最差解。指標R2、RMSE、SSE、P、v熵權值分別為0.190 002 311、0.202 248 703、0.202 519 782、0.202 444 871、0.202 784 333，由大到小的排序為SSE＞RMSE＞P＞v＞R2。依據(jù)啞變量模型各項評價指標相對最優(yōu)解的距離及其對應的熵權值計算得到距最優(yōu)解總距離。其數(shù)值越小代表模型優(yōu)度越高。在模型5 的b3參數(shù)中引入啞變量各指標距離最優(yōu)解距離總和最小為0.000 016 7，因此該模型為最優(yōu)啞變量模型，其表達式如下：

式（17）中：S1、S2分別代表采脂林分和未采脂林分。

3.3 模型擬合效果檢驗

使用檢驗數(shù)據(jù)比較最優(yōu)啞變量模型對采脂林分和未采脂林分的擬合效果，結果見表7。

表7 采脂與未采脂林分啞變量模型擬合效果Table 7 Fitting effect of dummy variable model for harvested and non-harvested forest

由表7可知，最優(yōu)啞變量模型對采脂與未采脂林分的R2值均超過了0.99，模型精度均大于99.5%，對采脂林分RMSE、SSE、v值分別為0.197 5、1.443 7、0.000 3，對未采脂林分RMSE、SSE、v值分別為0.566 9、12.213 7、0.002 2，說明該模型對采脂與未采脂林分斷面積擬合效果均較好，且模型對采脂林分斷面積擬合精度與模型準確性更好，模型對未采脂林斷面積擬合的相關系數(shù)更高。使用檢驗數(shù)據(jù)對最優(yōu)啞變量模型與最優(yōu)基礎模型擬合效果進行比較，結果見表8。

由表8可知，最優(yōu)啞變量模型較最優(yōu)基礎模型R2值于P值均有提升，其他低優(yōu)指標均下降，表明模型擬合效果在引入啞變量后得到了提升。圖1為使用檢驗數(shù)據(jù)得到的散點圖，分別為林分斷面積實測值與最優(yōu)啞變量模型預測值散點圖，及斷面積實測值與最優(yōu)啞變量模型殘差散點圖。

表8 最優(yōu)啞變量模型與最優(yōu)基礎模型擬合效果比較Table 8 Comparison of the effects of dummy variables and optimal foundation models

圖1中左側(cè)部分是根據(jù)檢驗數(shù)據(jù)中的斷面積實測值和最優(yōu)啞變量模型預測值建立的線性回歸方程，其R2值達到0.999 2，該回歸方程常數(shù)項較為接近0，相關系數(shù)值較為接近于1。右側(cè)部分是最優(yōu)啞變量模型殘差分布圖，從圖1中可看出殘差基本均勻分布在橫軸兩側(cè)。根據(jù)上述結果，表明該模型的擬合效果較好，可作為估測采脂與未采脂不同經(jīng)營措施下濕地松林分斷面積的參考。

圖1 最優(yōu)啞變量模型散點圖像Fig.1 Optimal dummy variable model scatter plot

4 結論與討論

利用樣地調(diào)查數(shù)據(jù)，以理查德方程、邏輯斯蒂模型、Mitscherlich 模型和Schumacher 模型為基礎模型建立福州市濕地松林分斷面積生長模型，結果表明Schumacher 模型為最優(yōu)基礎模型，其決定系數(shù)為0.985 4，模型精度為98.8%。經(jīng)過嘗試，在模型中不同參數(shù)加入啞變量時，模型各項評價指標存在差異，因此分別在b1、b2、b3、b4、b5、b6等參數(shù)中加入啞變量，其中在b3參數(shù)中引入啞變量后得到的模型擬合效果最好，其決定系數(shù)和模型精度較最優(yōu)基礎模型均提高，分別達到了0.999 1 與99.4%，且其他低優(yōu)指標較最優(yōu)都有下降，表明引入啞變量后模型擬合效果得到了提升。在參數(shù)b5中引入啞變量模型擬合效果略低于b3，但其模型精度和準確性最好分別為99.383%和0.005 95，且均方差、殘差平方和最小分別為0.663 746 579 和33.0。b3為與林分密度指數(shù)相關的參數(shù)，采脂與未采脂對樹木胸徑生長的影響較大，胸徑直接影響到林分的密度指數(shù)，因此，在b3參數(shù)引入啞變量較為合理，且模型擬合效果達到最優(yōu)，b4參數(shù)為樹高相關的參數(shù)，采脂與未采脂對樹木樹高生長有一定的影響但小于對胸徑的影響，所以在該參數(shù)引入啞變量的擬合效果較好但略低于b3。

采用含熵權值的TOPSIS 法對不同啞變量模型選優(yōu)，選用的各項評價指標熵權值由其本身數(shù)值決定，該方法具有較好的客觀性[24]，評價過程避免了人為主觀確定各個指標權重而影響綜合評價結果。本研究分別于最優(yōu)基礎模型參數(shù)b1、b2、b3、b4、b5、b6中引入啞變量得到不同啞變量模型，其評價指標R2、RMSE、SSE、P、v熵權值均不同，表明不同指標在啞變量模型中重要性存在差異，且各指標最優(yōu)解來自不同啞變量模型，因此有必要采用含熵權值的TOPSIS 法對啞變量模型進行綜合評價。

本研究考慮采脂與未采脂對濕地松林分斷面積生長的影響，從理論角度考慮，將采脂與未采脂的經(jīng)營措施作為啞變量加入到模型構建中，能規(guī)避由于經(jīng)營措施不同而出現(xiàn)模型預測有偏和不同經(jīng)營措施單獨建模不相容的情況。從實際工作的角度考慮，模型中引入啞變量能避免重復建模問題從而減小工作量，且模型使用起來更加便利。模型的選用不僅要考慮準確性還需要顧及實際操作簡易性，本研究選取了外業(yè)調(diào)查比較容易獲得的因子，以胸徑和樹高為主要因子加入了其他因子來提升模型擬合效果。但本研究還未構建相關的蓄積量模型，與已有的蓄積量模型可能會出現(xiàn)不相容的問題。今后的研究中可以進一步采集數(shù)據(jù)，使用度量誤差的方法求解模型，構建具有相容性的斷面積與蓄積量模型。