基于SPEI和信息擴(kuò)散加權(quán)一階局域法的旱澇特征分析

段春青，陳曉楠，許新勇，靳燕國(guó)

(1.北京市水務(wù)局 政務(wù)服務(wù)中心，北京 100071;2.中國(guó)南水北調(diào)集團(tuán)中線有限公司，北京 100038;3.華北水利水電大學(xué) 水利學(xué)院，河南 鄭州 450045）

旱澇是全球最常見(jiàn)的經(jīng)濟(jì)損失最嚴(yán)重的自然災(zāi)害。據(jù)統(tǒng)計(jì),每年旱澇災(zāi)害導(dǎo)致的經(jīng)濟(jì)損失占全球各類自然災(zāi)害造成損失的43%[1-2］。我國(guó)旱澇頻繁發(fā)生,且范圍廣、持續(xù)時(shí)間長(zhǎng),給生產(chǎn)、生活特別是農(nóng)業(yè)生產(chǎn)造成巨大損失。

旱澇災(zāi)害的時(shí)空演化特征分析是當(dāng)前災(zāi)害學(xué)領(lǐng)域研究熱點(diǎn),國(guó)內(nèi)外學(xué)者已取得較豐富的成果,并建立了一些科學(xué)實(shí)用的旱澇程度的量化評(píng)估指標(biāo),目前常見(jiàn)的有降水距平指標(biāo)[3］、降水溫度均一化指標(biāo)[4］、帕爾默干旱指數(shù)(palmer drought severity index,PDSI)[5］、標(biāo)準(zhǔn)化降雨指數(shù)(standard precipitation index,SPI)[6］、標(biāo)準(zhǔn)化降水蒸散指數(shù)(standard precipitation evapotranspiration index,SPEI)[7-8］等。

SPEI是Vicente-Serrano等于2010年提出的。該指數(shù)在SPI的基礎(chǔ)上,結(jié)合了SPI和PDSI的優(yōu)點(diǎn),加入潛在蒸散因素,具備計(jì)算簡(jiǎn)單和多時(shí)間尺度的特點(diǎn)。SPEI一般只需月平均氣溫和降水資料即可計(jì)算區(qū)域蒸散量和水分盈虧量,并以標(biāo)準(zhǔn)化指數(shù)呈現(xiàn),可用于不同時(shí)間、不同區(qū)域的旱澇變化特征比較。近些年,我國(guó)學(xué)者基于SPEI指數(shù)進(jìn)行了大量的旱澇評(píng)價(jià)和規(guī)律分析,取得了豐碩成果,同時(shí)驗(yàn)證了該指數(shù)在國(guó)內(nèi)具有很好的適用性[9-11］。

在旱澇時(shí)空分布特征和演變規(guī)律研究中,常見(jiàn)的分析方法有Mann-Kendall(M-K)趨勢(shì)檢驗(yàn)[12］、滑動(dòng)平均模型[13］、logistic函數(shù)[14］、多項(xiàng)式回歸[15］等。這些方法一般需要較多的樣本數(shù)據(jù)作為支撐,而實(shí)際應(yīng)用中所獲得的數(shù)據(jù)不多,一般只有幾十年資料,如何在較少樣本條件下得到較好的非線性回歸結(jié)果是重點(diǎn)研究問(wèn)題之一。信息擴(kuò)散技術(shù)是近些年興起的針對(duì)不完備樣本的一種有效處理手段,是利用適當(dāng)?shù)臄U(kuò)散函數(shù)實(shí)現(xiàn)對(duì)小樣本數(shù)據(jù)的集值化模糊數(shù)學(xué)方法,通過(guò)把單值樣本點(diǎn)擴(kuò)散到各個(gè)控制點(diǎn)上來(lái)獲取更多的信息[16-17］。此外,隨著混沌理論及相關(guān)時(shí)間序列分析技術(shù)的不斷發(fā)展,許多研究表明水文、氣象等系統(tǒng)具有一定的混沌特征,通過(guò)對(duì)時(shí)間序列重構(gòu)相空間,可挖掘其中豐富的動(dòng)力學(xué)信息[18-19］。傳統(tǒng)的各種混沌時(shí)間序列預(yù)測(cè)方法需基于大量數(shù)據(jù),對(duì)于樣本數(shù)據(jù)量較少情況,鑒于信息擴(kuò)散技術(shù)處理小樣本數(shù)據(jù)的優(yōu)勢(shì),可考慮在回歸中采用信息擴(kuò)散技術(shù)進(jìn)行數(shù)據(jù)處理,建立新的基于信息擴(kuò)散近似推理的混沌時(shí)間序列分析法。

1 研究方法

1.1 信息擴(kuò)散技術(shù)

1.1.1 信息擴(kuò)散頻率分布

信息擴(kuò)散頻率是在小樣本條件下,通過(guò)構(gòu)建離散論域的控制點(diǎn),利用信息擴(kuò)散函數(shù)將樣本點(diǎn)攜帶的信息分配到控制點(diǎn)上,得到模糊集合,在此基礎(chǔ)上分析頻率的分布,主要步驟如下[17］。

(1)設(shè)樣本序列為X＝(x1,x2,…,xl),l為樣本數(shù)量。根據(jù)樣本序列中的最大值和最小值,構(gòu)建離散論域U＝｛u1,u2…,uc｝,c為論域U的離散點(diǎn)數(shù)量,對(duì)每個(gè)樣本點(diǎn)xi通過(guò)下面的信息擴(kuò)散函數(shù),將其信息分配到論域中的每個(gè)控制點(diǎn)上:

h為信息擴(kuò)散系數(shù),可根據(jù)樣本長(zhǎng)度l和樣本中的a、b來(lái)計(jì)算。

式中:a、b分別為樣本中的最小值、最大值。

(2)對(duì)樣本點(diǎn)xi的信息分配結(jié)果進(jìn)行歸一化處理:

(3)針對(duì)控制點(diǎn)uj,根據(jù)所有樣本點(diǎn)在其上的信息分配,計(jì)算分配的信息總量:

(4)計(jì)算每個(gè)控制點(diǎn)的頻率值:

(5)計(jì)算各控制點(diǎn)的超越頻率值:

通過(guò)計(jì)算兩個(gè)控制點(diǎn)的超越頻率的差,即可計(jì)算出控制點(diǎn)區(qū)間內(nèi)的頻率分布。

1.1.2 信息擴(kuò)散回歸分析

信息擴(kuò)散回歸分析通過(guò)信息擴(kuò)散技術(shù)建立變量間的模糊關(guān)系,能較好地直觀反映出函數(shù)的變化趨勢(shì),得到光滑的函數(shù)曲線,主要步驟如下[18］。

(1)設(shè)自變量和因變量組成的樣本序列為[(x1,y1),(x2,y2),…,(xl,yl)］,分別根據(jù)輸入和對(duì)應(yīng)輸出樣本的數(shù)據(jù)分布情況,確定輸入、輸出的離散論域U＝｛u1,u2,…,us｝,V＝｛v1,v2,…,vr｝,并通過(guò)信息擴(kuò)散將樣本(xi,yi)轉(zhuǎn)換為模糊集合:

式中:hx、hy分別為輸入樣本和輸出樣本的信息擴(kuò)散系數(shù),根據(jù)式(2)來(lái)計(jì)算。

(2)對(duì)每組樣本(xi,yi)得到的模糊集Ai和Bi建立模糊關(guān)系Ri:

(3)對(duì)給定的輸入數(shù)據(jù)x0,通過(guò)信息擴(kuò)散回歸模型估算輸出數(shù)據(jù)y0。首先,將輸入數(shù)據(jù)進(jìn)行集值化處理:

式中:d＝｜u1-u2｜。

然后,根據(jù)得到的模糊集和模糊關(guān)系進(jìn)行模糊推理計(jì)算:

假設(shè)v′滿足

則v′是由第i組樣本(xi,yi)推算得出y0的一個(gè)估計(jì)值,記 作y＾i,并將其對(duì)應(yīng)的隸屬度看作估計(jì)值的權(quán)重wi。

(4)按照上述步驟,分別利用每組樣本進(jìn)行計(jì)算,將得到l個(gè)y0的估計(jì)值和相應(yīng)的權(quán)重,加權(quán)平均得到最終的估計(jì)值:

1.2 信息擴(kuò)散混沌時(shí)間序列分析

混沌是確定性系統(tǒng)中出現(xiàn)貌似隨機(jī)的現(xiàn)象,在自然界中普遍存在?；煦缣匦允菍?duì)初始值的敏感性,雖長(zhǎng)期不可預(yù)測(cè),但短期預(yù)測(cè)可行。目前,混沌時(shí)間序列的分析已經(jīng)廣泛應(yīng)用于各個(gè)領(lǐng)域中,主要預(yù)測(cè)方法有全域法、局域法、加權(quán)零階局域法、加權(quán)一階局域法等。其中,加權(quán)一階局域法最為常見(jiàn),在對(duì)時(shí)間序列數(shù)據(jù)重構(gòu)相空間后,采用線性局部擬合相空間中的樣本,一般需要較多的數(shù)據(jù)。當(dāng)時(shí)間序列樣本較少時(shí),重構(gòu)相空間后高維向量樣本更少,可考慮采用擅長(zhǎng)小樣本處理的信息擴(kuò)散回歸技術(shù)進(jìn)行處理,建立信息擴(kuò)散加權(quán)一階局域法,主要步驟如下[19］。

(1)設(shè)時(shí)間序列(xn｜n＝1,2,…,N),對(duì)其進(jìn)行快速傅里葉變換,計(jì)算平均軌道周期P。

(2)采用C-C方法[21］計(jì)算時(shí)間序列嵌入維數(shù)m和時(shí)間延遲τ,重構(gòu)相空間得到高維樣本序列Yi＝(xi,xi＋τ,…,xi＋(m-1)τ)(i＝1,2,…,M,其中M＝N-(m-1))。

(3)按照限制短暫分離的原則,要求｜i-j｜＞P。對(duì)每個(gè)點(diǎn)Yi尋找最鄰近點(diǎn)Yj,設(shè)最鄰近的“點(diǎn)對(duì)”初始距離為di(0)。

(4)對(duì)相空間中的每個(gè)最鄰近的“點(diǎn)對(duì)”,進(jìn)行k步演進(jìn)后,兩者之間的距離為di(k)。

(5)對(duì)于每個(gè)k,計(jì)算所有的ln[di(k)］的平均值y(k),由下式計(jì)算:

式中:y(k)為計(jì)算得到的序列；g為非0的di(k)的數(shù)量；T為采樣周期,本文取T＝1；K為實(shí)際中最多演進(jìn)的步長(zhǎng)。

(6)對(duì)序列y(k)進(jìn)行一元線性回歸,直線的斜率即為最大李雅普諾夫(Lyapunov)指數(shù)λ1。當(dāng)λ1＞0時(shí),認(rèn)為系統(tǒng)具有混沌特征,λ1越大混沌特征越明顯。若時(shí)間序列具有混沌特性,采用信息擴(kuò)散加權(quán)一階局域法進(jìn)行預(yù)測(cè)。

(7)尋找鄰近點(diǎn)。計(jì)算YM的q個(gè)最鄰近點(diǎn)Yk(k＝a1,a2,…,aq)。設(shè)Yk到Y(jié)M的距離為di(i＝1,2,…,q),dm是di中的最小值。定義Yk的權(quán)重Wk為:

(8)信息擴(kuò)散加權(quán)預(yù)測(cè)。利用信息擴(kuò)散回歸法根據(jù)YM中分量xN估算xN＋1,計(jì)算過(guò)程中將式(14)得到的權(quán)重Wk與式(12)權(quán)重相乘,并進(jìn)行歸一化處理后作為最終信息擴(kuò)散回歸中權(quán)重系數(shù),推算xN＋1的估計(jì)值。

1.3 標(biāo)準(zhǔn)化降水蒸散指數(shù)

根據(jù)區(qū)域的降水和氣溫?cái)?shù)據(jù)量化計(jì)算標(biāo)準(zhǔn)化降水蒸散指數(shù)(SPEI),計(jì)算主要步驟如下[2］。

(1)根據(jù)氣溫?cái)?shù)據(jù)及研究區(qū)位置,應(yīng)用Thornthwaite法計(jì)算月潛在騰發(fā)量ETi:

其中:

式中:ETi為月潛在騰發(fā)量；Ti為月平均氣溫；K為由區(qū)域緯度得出的平均晝長(zhǎng)修正系數(shù)；I為年熱量指數(shù)；m為指數(shù),由I計(jì)算得出。

(2)根據(jù)降水和計(jì)算的潛在騰發(fā)量分析水分盈虧量:

式中,Di為月水分盈虧量,正值意味著水分盈余,負(fù)值則表示水分虧缺,0代表水分恰好平衡；Pi為月降水量。

(3)基于log-logistic概率分布計(jì)算SPEI:

式中:Γ()為Gamma函數(shù)；N為序列數(shù)；s＝0,1,2。

設(shè)P＝1-F(x),當(dāng)P≤0.5時(shí):

當(dāng)P＞0.5時(shí):

其中,c0＝2.5155,c1＝0.8029,c2＝0.0103；d1＝1.4328,d2＝0.1893,d3＝0.0013。

SPEI的旱澇分級(jí)標(biāo)準(zhǔn)為:特旱(-∞,-2］,重旱(-2,-1.5］,中旱(-1.5,-1］,輕旱(-1,-0.5］,正常(-0.5,0.5］,輕澇(0.5,1］,中澇(1,1.5］,重澇(1.5,2］,特澇(2,∞］。

1.4 信息擴(kuò)散和混沌理論的旱澇分析

根據(jù)研究區(qū)域月尺度下的月均降水量和月均氣溫歷史數(shù)據(jù),利用SPEI指數(shù)計(jì)算得出月旱澇程度；利用旱澇時(shí)間序列數(shù)據(jù),基于信息擴(kuò)散頻率分布模型,得出不同旱澇程度分布規(guī)律；對(duì)旱澇時(shí)間序列進(jìn)行信息擴(kuò)散回歸分析,得出旱澇變化趨勢(shì)曲線,研究其隨時(shí)間交替演變的特征；應(yīng)用信息擴(kuò)散加權(quán)一階局域法預(yù)測(cè)旱澇變化趨勢(shì),驗(yàn)證模型的有效性。

2 應(yīng)用實(shí)例

以陜西省西安市為典型研究區(qū),分析旱澇演變規(guī)律。西安市是陜西省的政治、經(jīng)濟(jì)、文化中心,地處渭河流域中部,位于北緯33°39′—34°44′、東經(jīng)107°40′—109°49′。西安四季分明,屬于溫帶半干旱、半濕潤(rùn)大陸性季風(fēng)氣候區(qū),多年年均降水量約740 mm,年內(nèi)分布不均,年際變化大；多年平均氣溫約13℃；主要農(nóng)作物為冬小麥、夏玉米等。本研究選擇9月份為計(jì)算分析時(shí)段,該時(shí)段為玉米的成熟期,旱澇影響較大。研究區(qū)1951—2015年的月均降水量和月均氣溫?cái)?shù)據(jù)通過(guò)中國(guó)氣象數(shù)據(jù)網(wǎng)獲得。

2.1 旱澇等級(jí)頻率分布分析

根據(jù)SPEI指數(shù)計(jì)算西安市1951—2015年每年9月份的旱澇程度序列見(jiàn)圖1。

圖1 SPEI旱澇程度序列

以圖1中SPEI旱澇程度序列為基礎(chǔ),利用信息擴(kuò)散技術(shù)分析不同等級(jí)頻率分布。本文針對(duì)樣本數(shù)據(jù)的實(shí)際分布,選取離散論域控制點(diǎn)101個(gè),控制點(diǎn)旱澇程度最小值u1＝-2.026、最大值u101＝2.122,控制點(diǎn)等距離0.04148；其范圍包含全部樣本點(diǎn)。計(jì)算得出每個(gè)控制點(diǎn)上的超越頻率,并根據(jù)SPEI等級(jí)的劃分標(biāo)準(zhǔn)計(jì)算出不同旱澇等級(jí)分布頻率(見(jiàn)圖2)。

圖2 旱澇等級(jí)分布頻率

西安市9月份旱澇程度分布規(guī)律為:正常情況的頻率最大,為39.6%；特旱和特澇的頻率相對(duì)最小,分別為0.9%和1.0%；輕旱至輕澇的頻率達(dá)到63.0%；中旱及以上程度的干旱頻率為19.5%；中澇及以上程度的漬澇頻率為17.5%。

與傳統(tǒng)頻率計(jì)算方法進(jìn)行對(duì)比。利用傳統(tǒng)頻率計(jì)算方法計(jì)算得出不同旱澇等級(jí)分布頻率為:特旱1.5%、重旱6.2%、中旱13.8%、輕旱9.2%、正常40.0%、輕澇10.8%、中澇7.7%、重澇9.2%、特澇1.5%。信息擴(kuò)散計(jì)算結(jié)果與傳統(tǒng)方法整體上相差不大。但是,假設(shè)增加1個(gè)旱澇程度為-2(特旱等級(jí))的年份,數(shù)據(jù)樣本增加66組,則用傳統(tǒng)方法計(jì)算得出特旱等級(jí)的頻率為3.0%,相比原65組樣本計(jì)算結(jié)果增大了1倍。而利用信息擴(kuò)散方法計(jì)算結(jié)果為1.7%,是因?yàn)闃颖緮?shù)量較少,當(dāng)發(fā)生個(gè)別變化時(shí),傳統(tǒng)方法計(jì)算結(jié)果很不穩(wěn)定,但信息擴(kuò)散技術(shù)可充分利用每個(gè)樣本位置信息,得出較穩(wěn)定的結(jié)果,體現(xiàn)出小樣本條件下信息擴(kuò)散技術(shù)分析頻率的優(yōu)越性。

2.2 旱澇演化趨勢(shì)分析

利用信息擴(kuò)散回歸分析方法對(duì)旱澇程度數(shù)據(jù)進(jìn)行處理。將年份序號(hào)作為自變量,相應(yīng)的計(jì)算時(shí)段旱澇程度作為因變量,進(jìn)行關(guān)系擬合。在計(jì)算中,兩個(gè)變量所在論域的控制點(diǎn)個(gè)數(shù)均選取101個(gè),且等間隔分布,最大控制點(diǎn)和最小控制點(diǎn)均分別選取樣本的最大值、最小值外延0.001。信息擴(kuò)散近似推理完成擬合后,處理前后對(duì)比情況如圖3所示。圖3中原旱澇程度序列呈現(xiàn)鋸齒狀,逐年變化波動(dòng)大,不易直接看出演變規(guī)律。經(jīng)過(guò)信息處理后得到光滑曲線,能夠較為清晰地反映出演變的趨勢(shì)和波動(dòng)起伏的特征。

圖3 信息擴(kuò)散旱澇趨勢(shì)

由圖3可以看出,旱澇呈現(xiàn)交替變化規(guī)律:1951—1962年呈偏旱趨勢(shì),1963—1986年呈現(xiàn)偏澇趨勢(shì),1987—2002年又呈現(xiàn)偏旱趨勢(shì),2003年后再次呈現(xiàn)偏澇趨勢(shì)。1951—2015年的整體趨勢(shì)線基本處于水平狀態(tài),旱澇圍繞正常狀態(tài)上下波動(dòng),趨勢(shì)線未有明顯的上升或下降趨勢(shì)。

2.3 旱澇趨勢(shì)信息擴(kuò)散混沌預(yù)測(cè)

利用上述信息擴(kuò)散回歸得到的旱澇程度趨勢(shì)序列數(shù)據(jù),計(jì)算時(shí)間序列的Lyapunov指數(shù)最大值為0.1492、大于0,表明序列具有一定的混沌特性。利用本文提出的信息擴(kuò)散加權(quán)一階局域法對(duì)旱澇趨勢(shì)進(jìn)行預(yù)測(cè)分析,結(jié)果如圖4所示。

圖4 旱澇趨勢(shì)混沌預(yù)測(cè)

以1951—1984年實(shí)際趨勢(shì)數(shù)據(jù)為基礎(chǔ)預(yù)測(cè)1985年旱澇程度,之后再以1951—1985年實(shí)際趨勢(shì)數(shù)據(jù)為基礎(chǔ)預(yù)測(cè)1986年旱澇程度,以此遞推。對(duì)比1985—2015年實(shí)際趨勢(shì)數(shù)據(jù)和逐年預(yù)測(cè)的情況,可以看到利用信息擴(kuò)散混沌預(yù)測(cè)的結(jié)果與實(shí)際的旱澇變化過(guò)程吻合較好,預(yù)測(cè)曲線與實(shí)訓(xùn)曲線比較接近,曲線轉(zhuǎn)折點(diǎn)也較好地吻合,表明利用信息擴(kuò)散混沌分析方法對(duì)旱澇的趨勢(shì)預(yù)測(cè)可獲得較好效果,在實(shí)際防災(zāi)、減災(zāi)中能提供一定的參考。

利用信息擴(kuò)散一階加權(quán)法預(yù)測(cè)平均誤差為0.0863,而傳統(tǒng)線性回歸一階加權(quán)局域法預(yù)測(cè)平均誤差為0.1063,改進(jìn)后的方法預(yù)測(cè)的精度更高。

3 結(jié)論

根據(jù)西安市歷年9月份的氣象資料,采用SPEI指數(shù)對(duì)旱澇程度進(jìn)行量化計(jì)算。在此基礎(chǔ)上,利用小樣本數(shù)據(jù)分析的信息擴(kuò)散技術(shù)對(duì)旱澇等級(jí)的分布頻率進(jìn)行計(jì)算,得出不同旱澇等級(jí)對(duì)應(yīng)的發(fā)生頻率。利用信息擴(kuò)散回歸方法對(duì)旱澇程度序列進(jìn)行擬合,得到光滑的趨勢(shì)曲線,較清晰地反映出旱澇交替變化規(guī)律及整體變化趨勢(shì)。結(jié)合信息擴(kuò)散和傳統(tǒng)局域混沌時(shí)間序列分析方法,建立了信息擴(kuò)散加權(quán)一階局域法預(yù)測(cè)新模型,利用信息擴(kuò)散近似推理實(shí)現(xiàn)相空間里的高維向量擬合,并對(duì)旱澇趨勢(shì)進(jìn)行預(yù)測(cè),取得了較好的效果。研究表明,借助SPEI指數(shù)對(duì)旱澇量化評(píng)估,以及信息擴(kuò)散技術(shù)和混沌分析技術(shù)能夠較好地實(shí)現(xiàn)區(qū)域旱澇演變規(guī)律的分析,為當(dāng)?shù)胤罎晨购堤峁┘夹g(shù)支撐。