王正輝,孫會霞,郭 慶
(1.防空兵學院,河南鄭州 450052;2.河南工業(yè)大學,河南鄭州 450052;3.河南農業(yè)大學,河南鄭州 450052)
我國是一個農業(yè)大國,也是世界上農業(yè)自然災害較為嚴重的國家之一。農業(yè)自然災害主要包括水災、旱災、風災、蟲災、凍災等。建國60年來,我國平均每年發(fā)生重大水災、旱災、風災等自然災害23.2次,全國農業(yè)每年平均有0.4億hm2農作物、2億多人口受災,經濟損失達數百億元。因此預測農業(yè)自然災害十分必要。
災害的準確預測是農業(yè)的可持續(xù)發(fā)展條件之一,是災害防御決策的重要環(huán)節(jié)。據我國1950-2007年受災面積統(tǒng)計數據,繪制受災面積的折線圖 (圖1),從圖中可見我國農業(yè)受災面積變化既有不斷增長的趨勢,又有周期變化的規(guī)律,還有隨機波動成分。
圖1 1950-2007年我國農業(yè)受災面積的變化
假設災害發(fā)生強度為D(t),從長時間看,災害發(fā)生強度的變化既包含有確定的一面,也包含有隨機的一面,確定的部分表現為某種趨勢T(t)和周期規(guī)律P(t),隨機部分可視為在有限的狀態(tài)之間轉移,表現為隨機波動成分S(t)。于是混合預測模型可表達為D(t)=T(t)+P(t)+S(t),t=1,2,3,…。
因此可用混合預測模型作為受災面積預測公式,建立動態(tài)模型并預測未來狀態(tài)。首先通過建立t=1,2,…,55(即1950-2004年)的混合預測模型,預測 t=56,57,58(即 2005,2006和2007年)的發(fā)生值,并和這3年的統(tǒng)計值作比較,通過對比來檢驗模型的準確度。
一般用一元非線性回歸來擬合趨勢項,一元線性回歸x(t)=a+bt為一元非線性回歸的基礎。為了能夠應用一般最小二乘法,需要把一般非線性模型x(t)=F(t,a,b)轉化為某種線性形式。轉換后的一般式為X(t)=A+BT。
其中變換式X(t)=Fx[x(t)],A=Fa[a],B=Fb[b],T=Ft[t]。例如模型 x(t)=a ebt,兩邊取對數,就可變換為線性形式,其中X(t)=ln(x(t)),A=ln(a),B=b,T=t。然后可按一元線性回歸模型,通過最小二乘法計算出其中參數A,B,T,最后進行代換,計算出 a,b,t,把參數代入原表達式,就可以建立趨勢項擬合模型[1]。
由混合預測模型可知,周期成分與隨機波動的迭加為原始受災面積數據D(t)與趨勢項預測值之差,從圖2中可以看出,繪制的數據折線圖顯示了該差值序列具有一定的周期性,考慮周期可能存在迭加,可采用三角函數P(t)=C+A sin(ω1+φ1)+B cos(ω2+φ2),對周期分量數據進行擬合。式中A,B為波動幅度,C為直流成分,ω,ω1,ω2為頻率,φ,φ1,φ2為初始相位。由于表達式是非線性的,所以無法利用線性回歸來估計參數模型[2]。
確定部分包含趨勢項T(t)和周期規(guī)律P(t),把趨勢項T(t)和周期項P(t)迭加在一起組成一個確定部分模型,計算結果為:
Q(t)=T(t)+P(t)=1 115.416t0.380+115.786+594.935sin(0.354t+4.368)-462.591cos(0.628t+2.706)。
圖2 周期成分和隨機波動迭加折線的變化
設有隨機過程{Xn,n∈T},若對于任意的整數n∈ T和任意的 i0,i1,…in+1∈ I,條件概率滿足P{Xn+1=in+1|X0=i0,X1=i1,…,Xn=in}=P{Xn+1=in+1|Xn=in},則稱{Xn,n∈T}為馬爾科夫鏈,簡稱馬氏鏈。馬爾科夫鏈是一種時間離散、狀態(tài)可數的無后效隨機過程,當前的狀態(tài)和狀態(tài)轉移概率矩陣決定著下一個狀態(tài)。利用馬爾科夫鏈預測狀態(tài)轉移趨勢的步驟可以分4個環(huán)節(jié),分別為:狀態(tài)劃分、狀態(tài)轉移矩陣的建立、初始狀態(tài)概率向量的建立、狀態(tài)轉移預測[3-4]。
根據我國農業(yè)自然災害受災面積的隨機波動情況,可以把隨機擾動程度劃分為5個狀態(tài) (表1)。根據預測波動值等于狀態(tài)取值區(qū)間的平均值和確定部分預測值的乘積,計算出每一年的波動預測值,從圖3中可以看出混合模型預測值已經擬合得非常好。
條件概率pij(n)=P{Xn+1=j|Xn=i}為馬爾科夫鏈{Xn,n∈T}在時刻n的一步轉移概率,其中i,j∈I。相當于隨機游動的質點在時刻n處于狀態(tài)i的條件下,下一步轉移到狀態(tài)j的概率,狀態(tài)轉移概率矩陣是對離散馬爾科夫過程狀態(tài)轉移的定量描述[5]。其一步轉移概率矩陣 P1= [pij]n×n的計算公式為 Pij=Nij/Ni,i,j=1,2,…,n。
表1 隨機波動的狀態(tài)劃分
圖3 受災面積和混合模型預測值折線的變化
式中,Ni為序列中出現第i狀態(tài)的次數,Nij為Ni個狀態(tài)i中前一年狀態(tài)為i,下一年狀態(tài)為j的次數。同理可以建立k步狀態(tài)轉移概率轉移矩陣Pk。且有如下的關系:
Pk=P1k及 Pk=P1Pk-1,k=1,2,… 。
根據我國農業(yè)自然災害隨機狀分量劃分的狀態(tài)(表1),可以得出1-4步的轉移矩陣:
設某一時刻狀態(tài)分布矩陣為At,則狀態(tài)預測模型:
At+k=AtPk,k=1,2,… 。
以此可計算出未來每一年各種災害狀態(tài)出現的概率,從而可確定隨機分量的值,2002,2003,2004年隨機分量的劃分狀態(tài)為3,5,2,即狀態(tài)分布向量分別為,A2002=[0 0 1 0 0],A2003=[0 0 0 0 1],A2004= [0 1 0 0 0],預測計算結果列于表2,其中波動預測值等于用來預測年份的確定部分和狀態(tài)區(qū)間平均值的乘積。
表2 波動值預測
把上述建立的3類模型迭加起來,可得到我國農業(yè)自然災害集成預測模型:
D(t)=115.786+1 115.416t0.380+594.935sin(0.35t+4.368)-462.591cos(0.628t+2.706)+波動預測值。
當 t=1,2,…,55 時,可 以 計 算 出 1950,1951,…,2004年的我國農業(yè)自然災害擬合值;當t=56,57,58時,可以計算出2005,2006,2007年的農業(yè)自然災害預測值。和2005,2006,2007年的實際值對比 (表 3),相對誤差分別為10.77%,3.03%,8.39%,2006年的預測結果較好,2008年的尚可,2005年的則不太理想。
表3 模型預測值與實際發(fā)生值的比較
我國是每年發(fā)生農業(yè)災害較多且受災面積較大的國家,對農業(yè)受災面積的準確預測是農業(yè)可持續(xù)發(fā)展條件之一。但農業(yè)受災面積變化受多種因素的影響,且諸多因素之間是一種多變量、強耦合、嚴重的非線性關系,這種關系具有動態(tài)性。傳統(tǒng)的預測模型多數基于最小二乘法的單一數學公式,因而傳統(tǒng)的一些預測方法精度不高。本研究模擬結果表明,用集成預測模型預測我國自然災害的方法優(yōu)于傳統(tǒng)的單一數學模型。
[1] 魏瑞江.河北省主要農作物農業(yè)氣象災害災損評估方法[J].中國農業(yè)氣象,2000,21(1):27-31.
[2] Park M W,Kim Y D.Asystematic procedure for setting parameters in simulated annealing algorithms [J].Computer&Operations Research,1998,25(3):207-217.
[3] 張穎.軟計算方法[M].北京:科學出版社,2002.
[4] 宮德吉,陳素華.農業(yè)氣象災害損失評估方法及其在產量預報中的應用[J].應用氣象學,1999,27(1):66-69.
[5] 趙素英,劉艷.應用灰色模型進行初、終霜災害預測[J].黑龍江八一農墾大學學報,17(3):31-38.