江浙滬地區(qū)降雨極值的半?yún)?shù)擬合

徐鵬程

(南京財(cái)經(jīng)大學(xué) 應(yīng)用數(shù)學(xué)學(xué)院, 江蘇 南京 210046)

0 引言

近年來,極端氣候事件出現(xiàn)的頻率和強(qiáng)度均有所增加,造成了嚴(yán)重的經(jīng)濟(jì)損失和人員傷亡,引起了人們廣泛的關(guān)注.我國是一個洪澇災(zāi)害頻發(fā)的國家.研究顯示：過去半個世紀(jì)以來,我國越來越多的地區(qū)降水有集中趨勢[1].氣候要素是一個隨機(jī)變量,極值就是這些隨機(jī)的變量的某種函數(shù),研究氣候變化可以從研究極值入手,利用統(tǒng)計(jì)極值理論方法研究降水,特別是極端降水變化規(guī)律,對洪水風(fēng)險(xiǎn)評估、區(qū)域防災(zāi)減災(zāi)、防洪設(shè)施工程設(shè)計(jì)等具有重要意義和參考價(jià)值.

極值理論和分布擬合方法廣泛被運(yùn)用于分析降雨、洪水等極值數(shù)據(jù)的分布規(guī)律等研究.目前研究極端降水的統(tǒng)計(jì)分布模型有對數(shù)正態(tài)分布、廣義極值分布(GEV)、威布爾(Weibull)分布、耿貝爾(Gumbel)分布、Gamma分布、廣義帕累托分布(GPD)、皮爾遜Ⅲ(Pearson-Ⅲ)分布等.如趙立將GEV和Gumbel分布應(yīng)用于擬合五華縣降水極值[2];王文琪等人利用GEV和GPD分布研究太湖流域湖西浙西區(qū)降水極值特性[3];杜曉陽等人用皮爾遜Ⅲ、對數(shù)正態(tài)、指數(shù)、耿貝爾等概率分布擬合分析廣州市單日降水量的極值特征[4];張玉虎等人探討了不同的分布函數(shù)對不同降雨極值序列的適用性[5];張延偉等人選用GEV和Gumbel分布研究新疆地區(qū)降水概率分布特征[6];李亞麗使用Gumbel分布擬合陜西短歷時(shí)降水百年極值[7];張萍應(yīng)用廣義極值分布、正態(tài)分布和Gamma分布擬合珠江流域極端降水[8]; 王紅利用Weibull分布擬合廣西一日最大降水量的概率分布[9].此外,還有很多學(xué)者對各地降雨的時(shí)空變化特征進(jìn)行了研究.自從Qin和Zhang建立半?yún)?shù)模型以來,半?yún)?shù)密度比模型就成為現(xiàn)代統(tǒng)計(jì)的研究熱點(diǎn)[10].例如,Folkianos建立了半?yún)?shù)的密度估計(jì)方法[11].Qin給出了來自不同樣本分布的混合比例的半?yún)?shù)統(tǒng)計(jì)推斷[12].Zhang利用最大似然方法[13],研究了半?yún)?shù)的分位數(shù)估計(jì)方法.

不同的站點(diǎn)數(shù)據(jù)通常需要用不同的概率分布進(jìn)行擬合[14],從而尋找到該站點(diǎn)的最優(yōu)擬合分布模型,而選擇多種概率分布函數(shù)不便于分析降雨極值的分布規(guī)律.為了用同一類分布函數(shù)擬合分析不同的站點(diǎn),本文由經(jīng)典logistic回歸模型構(gòu)造出一個多樣本半?yún)?shù)概率密度比模型,將此模型用在降雨極值擬合上研究和分析長序列降雨極值.

1 多樣本概率密度比模型及其參數(shù)估計(jì)

1.1 k樣本logistic模型

k樣本概率密度比模型可以由經(jīng)典logistic回歸模型導(dǎo)出.設(shè)X11,X12,…,X1n1是來自第1個總體的樣本,Xk1,Xk2,…,Xknk是來自第k個總體的樣本;k組樣本相互獨(dú)立令D=i代表第i個總體,i=1,2,…,k,對給定的測試數(shù)據(jù)X=x,多樣本的logistic回歸模型是:

(1)

1.2 k樣本概率密度比模型的密度函數(shù)

類似于Qin等的處理[10],由貝葉斯法則得到模型(1)等價(jià)于如下的半?yún)?shù)密度比模型:

Xk1,Xk2,…,Xknk～gk(x)

(2)

1.3 k樣本概率密度比模型的參數(shù)估計(jì)

首先引入一些記號,令{T1,T2,…,Tn}代表合并的樣本{X11,…,X1n1,X21,…,X2n2,…,Xk1,…,Xknk},且令n=n1+n2+…+nk.參數(shù)估計(jì)的方法采用Owen提出的經(jīng)驗(yàn)似然方法[16]進(jìn)行,樣本的經(jīng)驗(yàn)似然函數(shù)為:

記Pi=dGk(Ti),這里,Pi可以理解成概率分布在每個觀測值上賦予的概率,即Pi是概率的躍變且總和為1,故

其計(jì)算由R語言[17]實(shí)現(xiàn),這里j=1,2,…,k-1,l(αj,βj)是關(guān)于(αj,βj)的剖面對數(shù)似然函數(shù):

1.4 k樣本概率密度比模型的分布函數(shù)

1.5 分布擬合的Kolmogorov-Smirnov距離

K-S檢驗(yàn)的檢驗(yàn)方法是將樣本數(shù)據(jù)的累計(jì)頻數(shù)分布與特定理論分布比較,如果兩者的K-S距離很小,可推斷理論分布和經(jīng)驗(yàn)分布擬合很好.設(shè)y1,y2,…,yn是樣本量為n的樣本且已按升序進(jìn)行了排序,K-S檢驗(yàn)的原假設(shè)H1:經(jīng)驗(yàn)分布的累計(jì)頻率來自于理論分布;備擇假設(shè)H2:經(jīng)驗(yàn)分布的累計(jì)頻率不是來自于理論分布.用F(x)和Fn(x)分別表示理論分布和經(jīng)驗(yàn)分布的累積頻率,原假設(shè)H1為真時(shí),K-S檢驗(yàn)統(tǒng)計(jì)量定義為D=maxx∈R|Fn(x)-F(x)|.K-S距離表達(dá)式為:

(3)

其計(jì)算由R語言實(shí)現(xiàn).

2 降雨極值數(shù)據(jù)的半?yún)?shù)密度比模型擬合

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

降雨極值數(shù)據(jù)來自中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)提供的逐年日降雨極大值資料(單位:1 mm),時(shí)間為1951-2018年(站點(diǎn)的起始時(shí)間不完全一致),范圍為江浙滬地區(qū)10個氣象站點(diǎn)氣象站點(diǎn)包括徐州、贛榆、南京、東臺、寶山、徐家匯、杭州、定海、溫州、瑞安等,以上站點(diǎn)序號依次為1-10.樣本量最大為68,最小為59.

2.2 半?yún)?shù)密度比模型對降雨極值的參數(shù)擬合

選用極大似然估計(jì)方法進(jìn)行參數(shù)的擬合.第一步選擇r(x)=(x,x2)T,求出樣本的經(jīng)驗(yàn)似然,根據(jù)拉格朗日乘子法,解計(jì)分方程組得出參數(shù).第二步將由R語言計(jì)算所得的參數(shù)代入,從而得到擬合半?yún)?shù)概率密度比模型的分布函數(shù).第三步將擬合半?yún)?shù)概率密度比模型的分布函數(shù)與樣本累積頻率進(jìn)行比較,由軟件可繪制出擬合圖，圖1給出了部分站點(diǎn)半?yún)?shù)密度比模型的分布函數(shù)與樣本累積頻率.由圖中可以直接看出半?yún)?shù)概率密度比模型的分布函數(shù)和樣本累積頻率接近,擬合效果較好.

(a)站點(diǎn)1 (b)站點(diǎn)3

(c)站點(diǎn)6 (d)站點(diǎn)9

2.3 與其他分布法的K-S距離比較

廣義極值分布(GEV)多用于模擬歷史降水和徑流,估算不同重現(xiàn)期最大值;Gumbel分布被廣泛應(yīng)用于單站頻率計(jì)算和特征值分布規(guī)律模擬中;Gamma分布是氣候統(tǒng)計(jì)學(xué)中重要的偏態(tài)分布,能夠穩(wěn)定的描述降水量的分布;Weibull分布則廣泛應(yīng)用于對各種極值數(shù)據(jù)的擬合.當(dāng)?shù)玫桨雲(yún)?shù)概率密度比模型的分布函數(shù)后,由式(3)分別計(jì)算半?yún)?shù)概率密度比模型(以下簡稱Semi分布)與擬合GEV、Gumbel分布、Gamma分布、Weibull分布的K-S距離,如表1所示.

表1 Semi、GEV、Gumbel、Gamma、Weibull分布法的K-S距離

比較上述5種分布的K-S距離折線圖,如圖2所示.可以觀察到Gumbel分布、Gamma分布以及Weibull分布的K-S距離明顯大于Semi分布和GEV分布的K-S距離.可以看出Semi分布和GEV分布更適合所有站點(diǎn)的數(shù)據(jù)擬合，再將Semi分布的K-S距離和GEV的K-S距離進(jìn)行比較,觀察到10個站點(diǎn)中,有6個站點(diǎn)的Semi分布的K-S距離小于GEV的K-S距離.一般來說,Semi分布略優(yōu)于GEV分布.

用5表示最大距離,4次之,以此類推,1表示最小,這樣得到10個站點(diǎn)各分布的最小距離次序,如表2所示.比較上述5種分布的K-S距離次序,Semi分布有5次最小,GEV分布3次,Gumbel分布和Weibull分布各1次.Semi分布的K-S距離次序的平均數(shù)和方差分別為1.90和1.21,而GEV分布的K-S距離次序的平均數(shù)和方差分別為2.00和1.33,均高于Semi分布,這表明Semi分布比GEV分布更加穩(wěn)定.

表2 5種分布法的K-S距離排序

圖2 5種分布法的K-S距離折線圖

3 結(jié)論

由經(jīng)典Logistic回歸模型導(dǎo)出多樣本半?yún)?shù)概率密度比模型,采用極大似然方法估計(jì)多樣本半?yún)?shù)概率密度比模型的參數(shù),并擬合逐年日降雨量序列,利用R語言計(jì)算K-S距離,并與其他4種常見的用于擬合降雨量極值數(shù)據(jù)的分布比較,得出以下初步結(jié)論:

1)半?yún)?shù)概率密度比模型的分布與廣義極值分布、Gumbel分布、Gamma分布、Weibull分布擬合優(yōu)度度量表明：半?yún)?shù)概率密度比模型的分布與廣義極值分布對大部分站點(diǎn)擬合的K-S距離總體偏小,而Gumbel分布、Gamma分布和Weibull分布對于部分站點(diǎn)擬合的K-S距離偏大;半?yún)?shù)概率密度比模型的分布與廣義極值分布比較,半?yún)?shù)概率密度比模型的分布的K-S距離更小也更穩(wěn)定.所以半?yún)?shù)概率密度比模型可以適應(yīng)不同站點(diǎn)的降雨極值提高擬合的穩(wěn)定性.

2)半?yún)?shù)概率密度比模型能更好的模擬降雨極值的概率分布,通過選取適合的r(x)可以提高擬合準(zhǔn)確率,有助于分析降雨極值的分布規(guī)律,是研究氣候變化下極值事件統(tǒng)計(jì)特征的重要工具,為估算降雨極值的重現(xiàn)期提供了有效方法.