邢俊景 曲智林

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

?

基于混合效應模型的大興安嶺地被可燃物含水率模型1)

邢俊景 曲智林

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

基于黑龍江省大興安嶺林區(qū)南甕河生態(tài)站的落葉松(Larixgmelinii)林、蒙古櫟(Quercusmongolica)林、落葉松白樺(Betulaplatyphylla)混交林3種典型林分的288組可燃物含水率數(shù)據(jù)，選擇基于平衡含水率的可燃物含水率實時變化模型為基礎模型，采用非線性混合效應(NLME)模型方法，以林分因子作為隨機效應，建立具有混合效應的可燃物含水率的實時變化預測模型，并通過給殘差方差增加權重的方法解決異方差性問題。結果表明:考慮隨機效應和異方差結構的可燃物含水率實時變化NLME預測模型的擬合效果(MAE=0.716 7，MRE=0.026 6)優(yōu)于不含隨機效應的可燃物含水率實時變化預測模型(MAE=0.815 6，MRE=0.031 2)；其中以常數(shù)加冪函數(shù)作為異方差結構的模型精度最高(AIC=547.72，BIC=581.29，-2LL=527.72)且明顯優(yōu)于未給殘差方差增加權重的可燃物含水率實時變化NLME預測模型(AIC=961.65，BIC=988.50，-2LL=945.65)。利用獨立樣本數(shù)據(jù)對模型進行檢驗，檢驗結果表明，對于可燃物含水率實時變化的預測，考慮隨機效應和異方差結構的NLME模型的檢驗精度(MAE=0.495 8，MRE=0.034 2)比利用最小二乘法擬合的多元非線性回歸模型(MAE=0.588 5，MRE=0.588 5)有所提高,說明基于混合效應模型的可燃物含水率實時變化模型可以很好地描述區(qū)域尺度上不同林分類型的可燃物含水率實時變化規(guī)律。

大興安嶺；可燃物含水率；非線性混合效應模型；異方差

Daxing’an Mountains; Fuel moisture content; Nonlinear mixed effects models (NLMEMs); Heteroscedasticity

森林可燃物是林火發(fā)生的物質基礎和首要條件[1]，其含水率多少是決定森林火災發(fā)生的重要指標之一,研究森林可燃物含水率變化的規(guī)律,有利于人類了解森林火災發(fā)生的規(guī)律。目前，可燃物含水率預測模型的方法主要包括4種：基于平衡含水率的方法、氣象要素回歸法、遙感估測法、基于過程模型的方法。在這4種方法中，基于平衡含水率的預測方法應用最廣，在美國、加拿大等國家的森林火險等級系統(tǒng)中已廣泛應用，且預測效果較好[2]。近年來，許多學者針對不同林分類型對可燃物含水率變化進行研究。金森等[3-4]針對江西、云南典型林分建立可燃物含水率預測模型。楊博文等[5]分析了帽兒山兩林分氣溫與地表可燃物溫度差異及對可燃物含水率預測的影響。周振超等[6]對哈爾濱城市林業(yè)示范基地典型林分地表死可燃物含水率與氣象因子的關系進行了研究。胡海清等[7]預測了大興安嶺典型林分地表死可燃物含水率動態(tài)變化。上述研究都是在林分尺度上對可燃物含水率進行研究，林分類型的不同決定了可燃物含水率變化規(guī)律不同，為了便于分析不同林分類型對地被可燃物含水率變化的影響，建立含有林分因子的可燃物含水率模型是非常重要的。

混合效應模型是指模型中部分參數(shù)或全部參數(shù)由固定效應和隨機參數(shù)的混合組成,該模型為分析分組數(shù)據(jù)提供了靈活和強有力的工具。因此，利用混合效應模型來分析不同林分因子的森林可燃物含水率變化規(guī)律是可行的。

基于黑龍江省大興安嶺林區(qū)南甕河生態(tài)站采集的落葉松林、蒙古櫟林、落葉松白樺混交林3種典型林分的地被可燃物含水率和相應的氣象因子數(shù)據(jù)，采用非線性混合效應(NLME)模型的方法，考慮林分因子的隨機效應和數(shù)據(jù)間的異方差性，建立可燃物含水率實時變化模型，并利用獨立樣本數(shù)據(jù)對模型的精度進行檢驗，以期建立具有混合效應的可燃物含水率預測模型，進而使區(qū)域尺度上不同林分類型的可燃物含水率模型具有相容性，為區(qū)域尺度上的林火預測提供理論依據(jù)。

1 研究區(qū)概況

研究區(qū)位于黑龍江省大興安嶺南甕河生態(tài)站(125°7′55″～125°50′5″E，51°5′7″～51°39′24″N)。該區(qū)為大興安嶺支脈，伊勒呼里山南坡，屬低山丘陵地貌，地形起伏不大，地勢為北高南低，西高東低，海拔高一般為500～800 m。氣候屬寒溫帶大陸性季風氣候，冬季受西伯利亞寒流的影響，異常寒冷，冬季漫長，長達9個月，年平均氣溫-3 ℃，極端最低溫度-48 ℃。相反，溫暖季節(jié)甚短，夏季最長不超過1個月，極端最高氣溫36 ℃。年降水量500 mm左右，80%以上皆集中于溫暖季節(jié)(7—8月)。而5—6月常有明顯旱象，形成云霧少，日照強，溫度低的氣候特點，致使林木草原火險增多。主要林分為落葉松(Larixgmelinii)、山楊(Populusdavidiana)、白樺(Betulaplatyphylla)、蒙古櫟(Quercusmongolica)、落葉松白樺混交林、山楊白樺混交林等。

2 研究方法

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

在南甕河生態(tài)站分別在落葉松林、蒙古櫟林、落葉松白樺混交林內設置樣地。在樣地內設置數(shù)據(jù)收集點，觀測儀器每小時自動收錄地被可燃物含水率，以及生態(tài)站實時觀測氣象數(shù)據(jù)，主要包括氣溫、空氣相對濕度、風速、降雨量等。本文所使用的數(shù)據(jù)為2015年6月11日00時—2015年6月14日23時的觀測數(shù)據(jù)，其中2015年6月11日00時—2015年6月13日23時的觀測數(shù)據(jù)數(shù)據(jù)作為建模數(shù)據(jù)，2015年6月14日00時—2015年6月14日23時的觀測數(shù)據(jù)作為檢驗數(shù)據(jù)，用于模型檢驗，統(tǒng)計信息見表1。數(shù)據(jù)處理利用R和SPSS 18.0軟件完成。

表1 建模數(shù)據(jù)和檢驗數(shù)據(jù)統(tǒng)計信息

注：Mt、Tt、Ht分別表示t時刻的可燃物含水率、氣溫和空氣相對濕度。

2.2 混合效應模型

降水對可燃物含水率的變化有明顯的影響，而降水量與可燃物含水率之間的關系比較復雜，本文建立的模型和使用的方法均基于無降水的情況。

由于可燃物含水率的瞬時變化率主要受可燃物的含水率與其平衡含水率的影響[8]。因此可燃物含水率變化模型為：

dMt/dt=k(Et-Mt)。

(1)

式中：Mt為可燃物t時刻的含水率；Et為可燃物t時刻的平衡含水率。

將方程(1)離散化，可得到可燃物含水率實時變化預測模型(1 h內變化)：

Mt+1=(1-k)Mt+k·g(Tt，Ht)。

(2)

式中：Tt為t時刻的溫度；Ht為t時刻的空氣相對濕度；Et=g(Tt，Ht)為環(huán)境因子對平衡含水率的響應函數(shù)。

Lindstrom et al[9]定義NLME模型為：

yij=f(φij，vij)+εij，i=1、…、N，j=1、…、n。

(3)

式中：yij=Mij+1為第i個林分j+1時刻的可燃物含水率；yij為第i個林分j+1時刻的可燃物含水率所對應的影響因子，包括第i個林分j時刻的可燃物含水率(Mij)、氣溫(Tij)和空氣相對濕度(Hij)；N為林分因子的個數(shù)，取值為3；ni為第i個林分可燃物含水率數(shù)據(jù)的組數(shù)，取值為72；εij為第i個林分j時刻觀測數(shù)據(jù)的誤差項值；φij是出現(xiàn)在非線性函數(shù)f中，由固定效應和隨機效應組成的形式參數(shù)向量，非線性函數(shù)f的表達式：

f=(1-k)Mij+k·g(Tij，Hij)。

(4)

2.3 形式參數(shù)構造與誤差結構

根據(jù)Pinheiro et al[10]在基礎模型中選擇不同的固定效應參數(shù)組合添加隨機效應組成混合參數(shù)進行模擬，選擇模型收斂及模擬精度最高的混合效應模型作為該水平的最優(yōu)模型，即采用赤池信息量準則(AIC)、貝葉斯信息準則(BIC)、負二倍的對數(shù)似然值(-2LL)等模型的評價指標，對模型的精度進行比較分析，這3個模型的評價指標值越小，其模擬精度越高。

為避免模型參數(shù)過多，本文利用似然比檢驗[11]：

LRT=2lg(L1/L2)=2(lgL1-lgL2)。

(5)

其中：L1為復雜模型的最大似然函數(shù)值；L2為簡單模型的最大似然函數(shù)值，LRT近似服從λ分布，且其自由度為k1-k2。給定顯著性水平α=0.05，當LRT≥λα(k1-k2)拒絕原假設，說明復雜模型與簡單模型差異顯著，選擇隨機效應參數(shù)較多的復雜模型；反之，接受原假設選擇隨機效應參數(shù)較少的簡單模型。

協(xié)方差矩陣R反映了數(shù)據(jù)間的異方差性，本研究試圖通過給殘差方差增加權重的方法消除數(shù)據(jù)間的異方差性。從自變量為當前時刻的含水率的指數(shù)函數(shù)、冪函數(shù)和常數(shù)加冪函數(shù)3個候選模型中由AIC、BIC、-2LL和似然比檢驗確定一個效果最好的殘差方差模型[12]。其結構表達式為：

(6)

Var(εij)=σ2exp(2δMij)；

(7)

(8)

(9)

其中：σ2為剩余方差；Ri為i林分樣地數(shù)據(jù)的協(xié)方差矩陣；為Gi為ni×ni維的對角矩陣，反映混合效應模型的隨機效應異方差性；Ii為ni×ni維的方差矩陣，描述了隨機效應的自相關性，假定Ii為ni×ni維的單位矩陣；Mij為第i林分樣地內第j時刻的可燃物含水率；δ、δ1、δ2為待估參數(shù)。

2.4 模型誤差分析

本研究利用含水率平均絕對誤差(MAE)、含水率平均相對誤差(MRE)這兩個指標對非線性多元回歸模型和NLME模型進行精度評價，MAE和MRE的值越接近0，說明模型的精度越高[13]。

(10)

(11)

3 結果與分析

3.1 基礎模型的確定

選擇四種不同的環(huán)境影響因子對可燃物平衡含水率的響應方程，采用最小二乘法擬合可燃物含水率實時變化模型，結果見表2。

表2 基于不同平衡含水率方程的可燃物含水率模型擬合結果比較

注：*表示一般顯著；** 表示顯著；*** 表示非常顯著。

從表2中可以看出模型Ⅰ的平均絕對誤差(MAE)和平均相對誤差(MRE)最小，從評價指標的大小方面考慮應該選擇模型Ⅰ為基礎模型。但是從系數(shù)的顯著性來看，模型Ⅰ的參數(shù)c的顯著水平大于0.1，顯著性不強，因此，選擇每一個參數(shù)顯著性水平都小于0.02，且平均絕對誤差(MAE)和平均相對誤差(MRE)相對最小的模型Ⅳ作為基礎模型：

Mij+1=(1-k)Mij+k(a+bHij+cTijHij)。

(12)

3.2 隨機效應參數(shù)的確定

以林分因子作為隨機效應，考慮不同隨機效應參數(shù)的組合，利用R語言的nlme功能，基于模型12建立NLME模型，并對模型的擬合精度進行比較。不同隨機效應參數(shù)組合的模型擬合精度見表3。

對于可燃物含水率模型，只有8種NLME模型收斂，其中模型12-1、12-2、12-3、12-4的AIC、BIC值完全相等，模型12-5、12-6、12-7、12-8的AIC、BIC值完全相等，并且各個模型差別不大。參數(shù)k代表水分的擴散系數(shù)，參數(shù)a、b、c與平衡含水率相關，從實際意義考慮林分因子對這4個參數(shù)均有影響，可以添加隨機效應，本文選擇以k和b同時作為混合參數(shù)的模型(12-5)建立可燃物含水率的NLME模型：

Mij+1=(1-(k+uik))Mij+(k+uik)(a+(b+uib)Hij+cTijHij)。

(13)

式中：uik為模型中參數(shù)k考慮林分因子的隨機效應參數(shù)，uib為模型中參數(shù)b考慮林分因子的隨機效應參數(shù)。

表3 基于不同隨機效應參數(shù)組合的可燃物含水率模型擬合精度比較

模型編號隨機效應參數(shù)AICBIC-2LL12-1k957.65977.79945.6512-2a957.65977.79945.6512-3b957.65977.79945.6512-4c957.65977.79945.6512-5k、b961.65988.50945.6512-6k、c961.65988.50945.6512-7a、c961.65988.50945.6512-8b、c961.65988.50945.65

3.3 異方差結構的NLME模型

借助模型的擬合值-殘差分布圖可以直觀地判斷出數(shù)據(jù)間是否存在異方差性[14]，圖1a、圖1b分別是傳統(tǒng)非線性多元回歸模型與考慮隨機效應和異方差結構的混合效應模型的可燃物含水率擬合值-殘差分布結果。從圖1a可以看出，采用最小二乘法擬合的可燃物含水率的殘差是隨著擬合值的增大而增大的，表明數(shù)據(jù)間存在異方差。為解決異方差問題，用冪函數(shù)、指數(shù)函數(shù)和常數(shù)加冪函數(shù)作為異方差結構對模型13進行改進，結果見表4。從表4中可知，AIC、BIC、-2LL等評價指標均明顯降低，說明考慮異方差結構的NLME模型明顯提高了模型的預估能力，且與不考慮異方差結構的NLME模型相比具有顯著差異，其中以常數(shù)加冪函數(shù)作為異方差結構的NLME模型(模型13-3)效果最好。利用建模數(shù)據(jù)對模型13-3進行擬合，最終得到可燃物含水率的NLME模型：

Mij+1=(1-(0.071+uik))Mij+(0.071+uik)(13.866+(0.182+uib)Hij-0.007TijHij)。

(14)

圖1b為模型14的擬合值-殘差分布圖，可以看出該模型能夠在一定程度上消除數(shù)據(jù)間的異方差。

a.傳統(tǒng)最小二乘法模擬的殘差分布 b.混合效應模擬的殘差分布

圖1 傳統(tǒng)非線性多元回歸模型和混合效應模型模擬結果殘差分布

3.4 模型檢驗

基于建模數(shù)據(jù)和檢驗數(shù)據(jù)分別計算模型12和模型14的MAE和MRE指標(見表5)。從表5可知，可燃物含水率的NLME模型(模型14)的擬合精度和檢驗精度均優(yōu)于基礎模型(模型12)。

表5 基礎模型與混合效應模型模擬結果比較

3.5 不同林分可燃物含水率的差異性

與一般的非線性回歸模型不同，混合效應模型考慮了不同林分類型可燃物含水率變化的差異，把誤差分解成兩部分，一部分為因林分因子不同而產生的隨機效應，另一部分是隨機誤差。由于混合效應模型中各參數(shù)具有明確的生物學意義，因此，根據(jù)各自隨機效應取值大小可反映出相應的隨機效應因子對可燃物含水率大小的影響程度[14]。由于模型參數(shù)中，固定效應參數(shù)恒取一個值，在相同氣象條件下(即溫度與相對濕度相同)對于同一個初始可燃物含水率，隨機效應參數(shù)uik越小，可燃物含水率的變化越??；而隨機效應參數(shù)uib越大，平衡含水率越大，可燃物含水率的變化越大。

表6 不同林分的隨機作用參數(shù)

表6是不同林分隨機作用參數(shù)值，由于這些隨機效應是由林分因子的差異產生的，因此它們的值可以反映在相同的氣象條件下當初始可燃物含水率相同時，經(jīng)過1小時后不同林分可燃物含水率的差異性。不同林分對水分擴散因子的影響程度由大到小的順序為：落葉松林、落葉松白樺混交林、蒙古櫟林；不同林分對平衡含水率的影響程度由大到小的順序為：蒙古櫟林、落葉松白樺混交林、落葉松林。

4 結論

1)利用單水平NLME模型方法，基于多元非線性回歸模型建立了大興安嶺典型林分地表可燃物含水率實時變化預測模型，并考慮了不同隨機效應參數(shù)的組合，用AIC、BIC、和-2LL及LRT對模型的擬合精度進行檢驗。結果表明：選擇以k和b同時作為混合參數(shù)建立的單水平NLME模型，從實際生物意義角度考慮是合理的，充分反映出林分因子的影響。

2)在模型建立過程中，通過殘差方差增加權重的方法，解決數(shù)據(jù)間的異方差問題，與傳統(tǒng)的多元線性回歸模型相比，考慮隨機效應和異方差結構的混合效應模型，可以很好地反映總體的可燃物含水率變化趨勢，還可以反映不同林分對可燃物含水率變化的差異，同時消除數(shù)據(jù)間的異方差性，進而提高模型的預估精度。

3)選擇大興安嶺林區(qū)3個典型林分的地被可燃物含水率的數(shù)據(jù)來建立的NLME模型，分析表明林分類型的不同對可燃物含水率變化的影響是不同的。由于數(shù)據(jù)的原因，對于大興安嶺林區(qū)其他林分類型本文沒有進行研究，因此，本文建立的混合效應模型只適用于文中提到的3種典型林分。

4)具有混合效應的可燃物含水率模型，是反映林分尺度上的可燃物含水率實時變化規(guī)律，如果去掉混合效應，則可以描述區(qū)域尺度上的可燃物含水率實時變化規(guī)律。如果具有混合效應的模型與各林分分別建模差別較小，并且各林分單獨建模的數(shù)據(jù)較少時，也可通過建立具有混合效應模型來實現(xiàn)。

[1] 田甜，邸雪穎．森林地表可燃物含水率變化機理及影響因子研究概述[J]．森林工程，2013，29(2):21-25.

[2] MATTHEWS S， MCCAW W L， NEAL J E， et al． Testing a process-based fine fuel moisture model in two forest types[J]. Canadian Journal of Forest Research，2006，37(1):23-35.

[3] 金森，劉萬龍.江西南昌典型林分地表死可燃物含水率預測：方法優(yōu)選與適用性分析[J].中南林業(yè)科技大學學報，2014，34(11)：1-8.

[4] 金森，劉周勇.昆明典型地表死可燃物含水率預測模型的研究[J].中南林業(yè)科技大學學報，2014，34(12)：7-15.

[5] 楊博文，陳鵬宇，金森.帽兒山兩林分氣溫與地表可燃物溫度差異及對可燃物含水率預測的影響[J].林業(yè)科學，2015，51(7)：91-98.

[6] 周振超，劉安婧，趙磊森，等.典型林分地表死可燃物含水率與氣象因子的關系：以哈爾濱城市林業(yè)示范基地典型可燃物為例[J].林業(yè)科技情報，2016，48(2)：4-11.

[7] 胡海清，陸昕，孫龍，等.大興安嶺典型林分地表死可燃物含水率動態(tài)變化及預測模型[J].應用生態(tài)學報，2016，27(7)：2212-2224.

[8] NELSON R M. Prediction of diurnal change in 10hour fuel stick moisture content[J]. Canadian Journal of Forest Research，2000，30(7):1071-1087.

[9] LINDSTROM M J， BATES D M． Nonlinear mixed effects models for repeated measures data[J]． Biometrics，1990,46(3):673-687.

[10] PINHEIRO J C， BATES D M． Mixed effects models in SandS-Plus[M]． New York: Spring-Verlag，2000.

[11] FANG Z， BAILEY R L． Nonlinear mixed effects modeling for slash pine dominant height growth following intensive silvicultural treatments[J]． Forest Science，2001,47( 3):287-300.

[12] 符利勇，孫華.基于混合效應模型的杉木單木冠幅預測模型[J].林業(yè)科學，2013，49(8)：65-74.

[13] 陸昕，胡海清，孫龍，等.大興安嶺地表細小死可燃物含水率預測模型[J].東北林業(yè)大學學報，2016，44(7)：84-90.

[14] 陳東升，孫曉梅，李鳳日.基于混合模型的落葉松樹高生長模型[J].東北林業(yè)大學學報，2013，41(10)：60-64.

邢俊景，女，1989年11月生，東北林業(yè)大學理學院，碩士研究生，E-mail：973063328@qq.com。

曲智林，東北林業(yè)大學理學院，教授。E-mail：q_zhilin@nefu.edu.cn。

2016年12月20日。

S762.1

1)林業(yè)公益性行業(yè)科研專項(201404402)。

責任編輯：王廣建。

Ground Surface Fuel Moisture Content by Mixed Effects Models in Daxing’an Mountains//Xing Junjing, Qu Zhilin(Northeast Forestry University, Harbin 150040, P. R. China)//Journal of Northeast Forestry University，2017,45(3)：58-62.

A real-time fuel moisture prediction model with the mixed effects was built with nonlinear mixed effect model (NLMEM) and stand factors as the stochastic effect. The fuel moisture model was with 288 sets of data related to the fuel moisture which were collected from three typical forests includingLarixgmelinii,QuercusmongolicaFischer and mixture ofLarixgmeliniiandBetulaplatyphylla. This model can solve the problems of heteroscedasticity by adding more weights to the residual variance. The results show that the fitting effect of the real-Time NLME prediction model (MAE=0.716 7,MRE=0.026 6) with the random effects and heteroscedasticity are better than that of the model (MAE=0.815 6,MRE=0.031 2) without considering the random effects. Moreover, the accuracy of the model (AIC=547.72, BIC=581.29, -2LL=527.72) using a constant plus power function as the heteroscedastic structure is the highest, and it is obvious superior to the real-time NLME prediction model (AIC=961.65, BIC=988.50, -2LL=945.65) without adding more weights to the residual variance. The results from the model test by using independent samples show that the test precision of the NLME model (MAE=0.495 8，MRE=0.034 2) with the random effects, and heteroscedasticity structure has some improvements in comparison to the multivariate nonlinear regression model (MAE=0.588 5,MRE=0.588 5) fitted by the least square method. The real-time fuel moisture content prediction model based on the mixed effects method can well describe the change laws in fuel moisture content for different types of forest in the regional scale of interest.