王 瑤,盧蕓瀟,劉 淼,2*
(1.伊犁師范大學數學與統計學院,新疆伊寧 835000;2.伊犁師范大學應用數學研究所,新疆伊寧 835000)
新疆位于我國西北邊疆,是我國荒漠面積最大、生態環境較為脆弱的地區;同時也是我國森林資源較少的省區之一,并且林地面積小、結構不合理、分布不均勻.火災干擾在新疆的生態系統中起著主導作用,改變了森林演替、生物地球化學循環和碳存儲.
尋找統計模型來描述這一地理區域火災和燃燒面積的年數趨勢是森林科學家非常感興趣的問題.事實上,對這些模型的詳細分析可能會揭示火災發生模式的新方面,并就火災機理提出重要的想法.有幾種統計方法可以根據火災數據揭示前兆地震活動.曼達拉茲和Ye發現了泊松概率模型,很好地描述了森林野火燃燒的過程[1-2].本文提出簡單的泊松模型,作為模擬野火發生的隨機性的一種實用方法;還利用復合泊松模型對大火災造成的損害進行了火災的總燒毀面積的建模;并通過二重復合泊松過程對新疆森林火災造成的經濟損失進行描述,希望所得結果對相關部門今后掌握新疆森林火災的變化規律提供一定的借鑒與參考.
1.1.1 時間序列
時間序列是指某一現象的統計指標值按時間順序所排列的數據,時間序列分析方法用于確定時間序列的周期、趨勢和概率結構.基于該結構建立的模型通常用于模擬隱含在時間序列中的過程,從而對未來進行預測,用于解決實際問題.
1.1.2ARIMA模型簡介
ARIMA(p,d,q)(簡稱ARIMA模型)全稱為差分自回歸移動平均模型(Autoregressive Integrated Moving Average Model),是由博克斯(Box)和詹金斯(Jenkins)于20世紀70年代初提出的一種著名的時間序列預測方法,所以又稱為Box-Jenkins模型、博克斯-詹金斯法.ARIMA模型是指將非平穩時間序列轉化為平穩時間序列,然后將因變量僅對它的滯后值,以及隨機誤差項的現值和滯后值進行回歸所建立的模型[3].
ARIMA模型是自回歸模型(AR模型)和移動平均模型(MA模型)的混合模型.其中,p為自回歸項;q為移動平均項;d為時間序列變成平穩時間序列時所做的差分次數.本文所使用的是非季節性的ARIMA模型,即
或簡寫為
其中,B為后移算子;?d為向后差分算子;φ(B)=1-φ1B-φ2B2-...-φpBp為自回歸算子;θ(B)=1-θ1B-θ2B2-...-θqBq為移動平均算子[4].
1.2.1 新疆森林燒毀面積的數據初步判斷
本文就表1 所示的2003~2019 年新疆森林火災年度燒毀面積[5],建立新疆森林火災燒毀面積的ARIMA(p,d,q)模型.
表1 2003~2019新疆森林火災數據
利用SPSS軟件對新疆2003~2019年的森林火災燒毀面積的數據進行分析.
圖1為2003~2019年新疆森林火災燒毀面積的時間序列圖.
圖1 2003~2019年新疆森林火災燒毀面積序列圖
從圖1中可以看出,2003~2019年新疆森林火災燒毀面積的序列圖是一個非平穩時間序列.本文對數據作二階差分后得到一個平穩的時間序列,見圖2.由圖2可以看出對原始序列作二階差分能使該時間序列基本達到平穩狀態.
圖2 二階差分后的時間序列圖
1.2.2 模型的識別與定階
利用ACF圖、PACF圖與貝葉斯信息準則(Bayesian Information Criterions,BIC)確定模型的階數,采用最大似然估計或最小二乘的方法對識別階段提供的模型進行參數估計,并假設檢驗,用以判斷模型是否恰當[6].利用統計軟件繪制自相關圖及偏相關圖(見圖3、圖4),進行階數的初步判斷.
圖3 二階差分后的自相關圖
圖4 二階差分后的偏相關圖
對該時間序列取不同參數值來比較BIC值(BIC值越小越好),以及比較平穩的R2值和MAE值,反復擬合取(p,d,q)=(2,2,0)時,時間序列模型的BIC值達到最小值12.113,R2=0.611,模型擬合效果比較顯著.
1.2.3 模型的參數估計
取(p,d,q)=(2,2,0)時,模型的參數估計結果如表2所示.由表2可以看出,回歸系數的p值都小于0.05,接近于0,于是得到結論:t檢驗效果顯著.
表2 參數估計值
1.2.4 模型的顯著性檢驗
模型的顯著性檢驗是對序列的原始數據與擬合數據的誤差序列進行相關檢驗,看其是否與實際相吻合,是否能夠很好地反映實際,對時間序列模型的檢驗即就是看殘差序列是否為白噪聲序列[7].本文檢驗是否為白噪聲序列的方法選擇觀測殘差序列的ACF圖和PACF圖進行檢驗,見圖5.
圖5 ARIMA模型殘差的相關函數圖
由圖5可知,殘差的ACF和PACF都是平穩的,兩者的數值都趨向于0,但都不等于0,并且殘差序列數值之間沒有相關性.因此殘差序列近似為白噪聲序列.于是得到的模型ARIMA(2,2,0)能夠較好地用來擬合新疆2003~2019年森林火災燒毀面積的時間序列.
1.2.5ARIMA(2,2,0)模型的預測
利用建立的模型預測新疆森林火災燒毀面積,得到預測值與實際值的對比結果,如圖6:
圖6 新疆森林火災燃燒面積預測值與實際值對比圖
對比圖6中的原始序列和實際序列,模型的預測是合適的.使用該模型對未來兩年(2020~2021年)新疆森林火災燒毀面積進行預測,結果如表3 所示.由表3 預測得到2020 年新疆森林火災燒毀的面積約為37.57 hm2;預測2021年新疆森林火災所致的燒毀面積約為43.82 hm2.
表3 新疆2020~2021年森林火災燒毀面積預測值
定義1[8]稱計數過程{N(t),t≥0},為具有參數λ>0的泊松過程,如果{N(t),t≥0} 滿足下列條件:
(1)N(0)=0;
(2)N(t)是獨立平穩增量過程;
(3)P{N(h)=1}=λh+o(h);
(4)P{N(h)≥2}=o(h).
定義2[9]設{N(t),t≥0} 是參數為λ的泊松過程,若其中{Xn,n=1,2,...} 為獨立同分布隨機變量序列且與{N(t),t≥0} 獨立,則稱{Y(t),t≥0} 為復合泊松過程.
定義3[10]若{Y(t),t≥0} 是一復合泊松過程,{ηn,n=1,2,...} 是獨立同分布的非負整數值隨機變量序列,且{ηn} 與{Y(t)}相互獨立.令,t≥0,則稱{Z(t),t≥0} 為二重復合泊松過程.
定理1若{N(t),t≥0} 是參數為λ的泊松過程,且E(X12)<∞,則復合泊松過程,t≥0的基本數字特征如下:
(1)特征函數為gY(t)(u)=exp{λt[gX(u)-1]},其中gX(u)是隨機變量X1的特征函數,λ表示速率;
(2)均值函數為E[Y(t)]=λtE(X1);
(3)方差函數為D[Y(t)]=λtE(X12).
由于{N(t)≥0,t≥0} 與{Xn,n=1,2,...},是相互獨立的,則
(2)根據特征函數與矩的關系得到數學期望為
(3)均方值為
于是Y(t)的概率生成函數為
若s=0,gY(t)(s)的k階導數如下:
當m=1,x=0,1 且n=1,2,...時,P(Xn=0)=p0,P(Xn=1)=1-p0,此時易見Y(t)是參數為λ(1-p0)t的泊松過程.
當m=2時,Y(t)的概率生成函數為
當m>2時,利用上述同樣的方法進行微分及代數運算,得到m>2的概率
通過歸納得
2.3.1 模型的建立
從《新疆統計年鑒》中收集到2003~2019年有關新疆森林火災的數據如圖7.
圖7 2003~2019年新疆森林火災發生次數
利用圖7 中的數據進行分析.令{N(t),t≥0} 為某一給定季節內新疆地區森林火災發生次數的齊次泊松過程,則相應的森林火災的燒毀面積Xn,n=1,2,...為獨立同分布隨機變量,且與{N(t),t≥0} 獨立,則在[0,t]時間段內,新疆森林火災燒毀的總面積可看作一個復合泊松過程Y(t)=X1+X2+...+XN(t)=當每次新疆地區森林火災燃燒面積相互獨立時,即可確定新疆森林火災燒毀總面積的概率密度函數.在t時間段內新疆森林火災造成的經濟損失可看成一個二重復合泊松過程.當每一單位面積燒毀造成的經濟損失相互獨立時,便可得到新疆森林火災造成總經濟損失的數學模型.
2.3.2 模型檢驗
本文利用Kolmogoroν‐Smirnoν檢驗統計量(ρ=0.095),得到{N(t),t≥0} 是一個齊次泊松過程,λ=36.88(年),這意味著每年新疆森林火災預期達到36.88次.從圖8和Kolmogoroν-Smirnoν檢驗中得到了新疆森林火災燃燒面積的頻率分布,隨機變量Xn,n=1,2,...為對數正態分布.根據Spearman-ρ檢驗(Spearman-ρ=0.204;p=0.433),滿足獨立性假設條件.因此得到是一個復合泊松過程.利用隨機變量Xn與N(t)的分布特性同樣可得新疆森林火災燒毀總面積的復合泊松過程的概率函數.
圖8 新疆森林年度燃燒面積圖
通過Kolmogoroν-Smirnoν檢驗得到隨機變量ηn,n=1,2,...為一個正態分布.根據Spearman-ρ檢驗(Spearman-ρ=0.074;p=0.859),我們接受獨立性假設.因此我們認為是一個二重復合泊松過程.
本文首先介紹了一種時間序列模型,用于對新疆未來兩年森林火災的燃燒面積進行預測;同時提出了一種復合泊松模型,用于模擬給定季節內新疆森林火災的年度燃燒面積和火災次數的分布.利用新疆2003~2019年森林火災的歷史數據來擬合模型,最后得到新疆森林火災燃燒總面積的復合泊松過程的數字模型和新疆森林火災所造成經濟損失的二重復合泊松過程的數學模型,從而說明建立復合泊松模型,可以很好地對新疆森林火災燃燒總面積和燃燒次數進行概率擬合.本研究中得到的新疆森林火災的分布規律對于未來相關部門防控和管理新疆森林火災系統有一定的參考價值.