基于運(yùn)動(dòng)擴(kuò)散模型的含零值降水序列頻率計(jì)算

李丹丹，宋松柏，李運(yùn)平

（西北農(nóng)林科技大學(xué) 水利與建筑工程學(xué)院，陜西 楊凌　712100）

1　研究背景

水文頻率計(jì)算是根據(jù)水文現(xiàn)象的統(tǒng)計(jì)特性，利用現(xiàn)有水文資料，分析水文變量設(shè)計(jì)值與出現(xiàn)頻率（或重現(xiàn)期）之間的定量關(guān)系的工作過程，對于水利工程規(guī)劃、設(shè)計(jì)、管理以及水資源配置等具有至關(guān)重要的作用［1］。降水設(shè)計(jì)值是進(jìn)行用水管理的主要依據(jù)［2］。通常情況下，水文頻率計(jì)算采用的水文序列不含有零值，但在干旱和半干旱地區(qū)的水文頻率分析中，經(jīng)常遇到包含若干個(gè)零值的水文序列，而零是變異范圍的下限［3-6］。這些水文序列可能是干旱地區(qū)的月降水量，水庫設(shè)計(jì)中的防洪庫容系列以及河道斷流等［7］。因此，常規(guī)的水文頻率計(jì)算方法無法進(jìn)行這類水文序列頻率計(jì)算［8］，而且以往研究中對于這類水文序列的頻率分析關(guān)注較少［9］。Jennings等［4］、Woo等［9］、Wang 等［5］提出了增加小水文值、忽略零值、全概率法以及應(yīng)用次序統(tǒng)計(jì)量構(gòu)造全序列分布等方法計(jì)算含零值水文序列的頻率。1980年代以來，我國一些學(xué)者也進(jìn)行了含零值水文序列的計(jì)算方法研究，如間接法、中值適線法和Ⅱ型乘法分布等方法［7，11-14］?！端姽こ淘O(shè)計(jì)洪水計(jì)算手冊》推薦使用頻率比例法和Ⅱ型乘法分布法［11，15-17］。張良等［18］利用這兩種方法對安固里河洪峰流量進(jìn)行了分析。2003年，Strupczewski等［8］提出了基于運(yùn)動(dòng)擴(kuò)散模型（Kinematic Diffusion Model，簡稱KD模型）的含零值序列頻率計(jì)算模型，這是一種新型的含零值水文序列頻率計(jì)算方法，該模型已成功地應(yīng)用于美國干旱地區(qū)洪水頻率計(jì)算，取得了較好的擬合效果。

目前，Strupczewski等推導(dǎo)了KD模型矩法、極大似然法參數(shù)計(jì)算公式。眾所周知，概率權(quán)重矩法是水文頻率分布參數(shù)計(jì)算精度較高的方法［19-21］。但是，KD模型目前沒有概率權(quán)重矩法參數(shù)計(jì)算公式，且這類模型尚未用于我國含零值水文序列頻率計(jì)算。本文從概率權(quán)重矩定義出發(fā)，應(yīng)用數(shù)學(xué)變換原理和數(shù)值計(jì)算原理，嚴(yán)格地推導(dǎo)了概率權(quán)重矩法推求KD模型參數(shù)計(jì)算公式，給出了相應(yīng)的數(shù)值計(jì)算方法，并與矩法、極大似然法估計(jì)KD模型參數(shù)進(jìn)行比較，為KD模型參數(shù)估計(jì)方法選擇提供依據(jù)。以陜西省6個(gè)測站2月份降水?dāng)?shù)據(jù)為例，驗(yàn)證概率權(quán)重矩法計(jì)算KD模型參數(shù)公式正確性，并與頻率比例法和II型乘法分布法進(jìn)行對比，研究KD模型在我國含零值降水序列頻率計(jì)算的普適性，以期為我國含零值水文序列頻率計(jì)算提供一種新的計(jì)算方法。

2　含零值水文序列概率模型及運(yùn)動(dòng)擴(kuò)散模型

2.1含零值水文序列概率模型基于概率論觀點(diǎn)，含零值水文序列中零值出現(xiàn)的概率不等于零，在概率密度函數(shù)圖中具有不連續(xù)性，因此，這類水文序列的概率密度函數(shù)可表示為［8］：

式中：β為零值事件發(fā)生的概率，即可按泊松過程進(jìn)行估計(jì)；為連續(xù)函數(shù)，且滿足R為矢量參數(shù)；δ(x)為Dirac delta函數(shù)，即

2.2運(yùn)動(dòng)擴(kuò)散模型根據(jù)完全線型圣維南方程可推得表示運(yùn)動(dòng)波的運(yùn)動(dòng)擴(kuò)散模型（KD模型），其脈沖響應(yīng)函數(shù)可表示為［8］：

式中：Pi(λ)服從泊松分布；服從Gamma分布；參數(shù)α，λ為關(guān)于河道形態(tài)和水流條件的函數(shù)；1(t)為單位階躍函數(shù)，

Strupczewski等［22］指出，對于有限河段而言，基于物理參數(shù)的馬斯京根模型等價(jià)于KD模型，式（2）可以表示為概念元素的網(wǎng)絡(luò)級聯(lián)，即線性水庫和水文中常用的線性河道。Einstein［23］指出式（2）所給函數(shù)可作為推移質(zhì)輸沙的確定性隨機(jī)模型。Eagleson［24］提出，在假設(shè)暴雨發(fā)生次數(shù)遵循泊松分布、暴雨深度服從指數(shù)分布的條件下，式（2）可以作為總降水量頻率計(jì)算的分布函數(shù)。因此，Strupczewski等［22］認(rèn)為式（2）可以用于含零值水文序列頻率計(jì)算。

3　運(yùn)動(dòng)擴(kuò)散模型概率密度函數(shù)推導(dǎo)

將KD模型脈沖響應(yīng)函數(shù)變量t替換為變量x，得到雙參數(shù)概率分布函數(shù)［8］：

式中：z為泊松分布隨機(jī)變量；x為Gamma分布隨機(jī)變量。

式（3）等號右側(cè)第一項(xiàng)可表示為：

根據(jù)第一類一階修正Bessel函數(shù)及x＞0時(shí)1(x)=1，第二項(xiàng)可表示為：

將式（4）—（5）代入等式（3）得：

式（6）即為KD模型的概率密度函數(shù)，簡稱KD-PDF。因此，只要參數(shù)α、λ已知，就可以通過KD模型求得零值概率和非零值事件發(fā)生概率。

4　運(yùn)動(dòng)擴(kuò)散模型參數(shù)估計(jì)及設(shè)計(jì)值推求

4.1參數(shù)估計(jì)

（1）矩法參數(shù)估計(jì)根據(jù)累計(jì)量與矩的關(guān)系，可得［8］：

式中：μ′1為一階原點(diǎn)矩；μ2為二階原點(diǎn)矩；Cv為變差系數(shù)；μ′1、μ2和Cv可用樣本估算值代替。

（2）極大似然法參數(shù)估計(jì)根據(jù)Strupczewski等［8］的推導(dǎo)，KD模型分布參數(shù)與極大似然估計(jì)值的關(guān)系為：

式中：n為序列總長度；為n階貝塞爾函數(shù)；xj為第j個(gè)樣本值。

顯然，式（10）為λ的函數(shù)，求解該非線性方程，即可獲得參數(shù)，代入式（9）可獲得參數(shù)。其初值確定可以通過零值概率推求。

（3）概率權(quán)重矩法參數(shù)估計(jì)根據(jù)概率權(quán)重矩的定義，KD模型的零階概率權(quán)重矩為：

KD模型一階概率權(quán)重矩為：

式中：

將式（6）和式（15）代入式（13）中有：

則M1=w1+w2，經(jīng)化簡可得：

令g=2s，當(dāng)s=0時(shí)，g=0，當(dāng)s=∞時(shí)，則

因此，KD模型概率權(quán)重矩為：

概率權(quán)重矩法與傳統(tǒng)矩法一樣，用樣本概率權(quán)重矩作為總體概率權(quán)重矩的估計(jì)量，而樣本概率權(quán)重矩與分布無關(guān)。因此KD模型參數(shù)與概率權(quán)重矩關(guān)系為：

式中：l為將樣本按由小到大順序排列的序號；為相應(yīng)于序號l的樣本值；其余參數(shù)意義同前。

4.2設(shè)計(jì)值推求對于給定的xp，其超越概率根據(jù)式（23）計(jì)算。

5　頻率比例法與Ⅱ型乘法分布法

5.1頻率比例法頻率比例法在進(jìn)行頻率分析時(shí)，只分析非零的部分序列，然后通過頻率轉(zhuǎn)化得到全序列的頻率曲線，又稱間接法［14］。

首先，將n2項(xiàng)非零資料視作一個(gè)獨(dú)立系列，按P-Ⅲ型分布（Cs=2Cv）進(jìn)行頻率分析計(jì)算，各點(diǎn)的經(jīng)驗(yàn)頻率為然后通過參數(shù)計(jì)算、調(diào)整和適線等常規(guī)方法，求得一條頻率曲線。但此線只能代表全部n項(xiàng)資料中一部分資料的分布情況，而任意一個(gè)樣本值xi的頻率，可以通過下列轉(zhuǎn)換求得：

式中：Pb相應(yīng)于全部系列的頻率；Pa相應(yīng)于n2項(xiàng)系列的頻率。點(diǎn)繪各xi點(diǎn)的Pb，可以連成一條新的頻率曲線，即為所求的頻率曲線。

5.2Ⅱ型乘法分布法Ⅱ型乘法分布法是從獨(dú)立隨機(jī)變量的概率乘法定理引伸出來的。對于標(biāo)準(zhǔn)化變量模比系數(shù)及均值頻率b=P(K≥1)，嚴(yán)格服從概率乘法定理的原始乘法分布即［14］：

考慮到充分利用均值頻率b的有益信息，以b值為分界把(0，1)頻率區(qū)間分為首尾兩部分，建立起Ⅱ型乘法分布，其概率密度函數(shù)為：

式中：a、c分別為首、尾部密度函數(shù)參數(shù)，b的變幅一般在0.2～0.8之間。

Ⅱ型乘法分布法中均值頻率b可以從樣本系列求出均值后，在經(jīng)驗(yàn)頻率曲線上查讀或內(nèi)插b=P(K≥1)。對于確定的b值，密度參數(shù)a和c可直接由式（28）估計(jì)。

在(0，1］區(qū)間內(nèi)的每一個(gè)尾部變量均可單獨(dú)提供參數(shù)c的信息，當(dāng)k=0時(shí)，用k=0.00001代替。在［1，∞)區(qū)間內(nèi)的每一個(gè)首部變量，均可以單獨(dú)提供參數(shù)a的信息。

6　模型應(yīng)用

以陜西省寶雞、涇陽、咸陽、乾縣、武功、興平6個(gè)測站2月份降水序列為研究對象，在矩法、極大似然法估計(jì)KD模型參數(shù)的基礎(chǔ)上，驗(yàn)證概率權(quán)重矩法在KD模型參數(shù)估計(jì)中的應(yīng)用效果，并與頻率比例法、Ⅱ型乘法分布模型進(jìn)行對比，采用AIC準(zhǔn)則、OLS準(zhǔn)則和殘差平方和（RSS）最小準(zhǔn)則進(jìn)行擬合優(yōu)度評價(jià)，研究KD模型的適用性。經(jīng)審查，資料滿足一致性要求，資料序列長度如表1所示。

表1　陜西省6測站2月份降水資料系列

6.1參數(shù)估計(jì)利用矩法和極大似然法以及概率權(quán)重矩法對KD模型參數(shù)α和λ進(jìn)行計(jì)算；在計(jì)算樣本均值、標(biāo)準(zhǔn)偏差以及變差系數(shù)基礎(chǔ)上，根據(jù)非零序列推求P-Ⅲ型分布參數(shù)α、β、γ；根據(jù)式（28）計(jì)算Ⅱ型乘法分布模型參數(shù)a、b、c。研究區(qū)6個(gè)測站分布參數(shù)見表2。

表2　陜西省6測站降水量分布3種計(jì)算方法參數(shù)估計(jì)結(jié)果

由表2可知，KD模型的3種參數(shù)估計(jì)方法所估計(jì)的參數(shù)值存在一定差異，但又有一定規(guī)律。3種參數(shù)估計(jì)方法中，極大似然法估計(jì)得到的α值最大，λ值最小，矩法估計(jì)得到的λ值最大，α值最小，概率權(quán)重矩法得到的得參數(shù)值介于二者之間?；赑-Ⅲ分布的頻率比例法在參數(shù)計(jì)算過程中，位置參數(shù)γ為負(fù)值。Phien等［25］提出，利用P-Ⅲ分布進(jìn)行降水量頻率計(jì)算時(shí)，若獲得γ為負(fù)值時(shí)，通常將γ設(shè)置為零，轉(zhuǎn)換為兩參數(shù)Gamma分布計(jì)算。經(jīng)計(jì)算，基于Gamma分布的頻率比例法參數(shù)見表3。

表3　基于Gamma分布的頻率比例法參數(shù)

6.2頻率曲線繪制根據(jù)參數(shù)估計(jì)結(jié)果，利用KD模型、頻率比例法以及Ⅱ型乘法分布法模型對6個(gè)測站含零值降水序列進(jìn)行頻率計(jì)算，并繪制頻率曲線，如圖1所示。

由圖1可以看出，KD模型對于含零值降水序列的擬合效果優(yōu)于頻率比例法與Ⅱ型乘法分布法。KD模型與Ⅱ型乘法分布法對于尾部擬合效果優(yōu)于頻率比例法，其主要原因是KD模型與Ⅱ型乘法分布法在頻率計(jì)算過程中將含零值水文序列作為一個(gè)整體，近似為連續(xù)分布進(jìn)行計(jì)算，而頻率比例法是將含零值水文序列分為零序列和非零序列兩部分，通過非零序列計(jì)算含零值水文序列。

6.3擬合優(yōu)度評價(jià)采用OLS準(zhǔn)則、AIC準(zhǔn)則和RSS最小準(zhǔn)則對KD模型、頻率比例法、Ⅱ型乘法分布法頻率曲線進(jìn)行擬合優(yōu)度評價(jià)，計(jì)算結(jié)果見表4。

圖1　陜西省6測站2月降水量頻率曲線

由表4可以看出，基于OLS準(zhǔn)則、AIC準(zhǔn)則和RSS最小準(zhǔn)則的曲線擬合優(yōu)度評價(jià)結(jié)果基本一致，且與頻率曲線圖所反映的結(jié)果相同，其中RSS結(jié)果明顯反映出KD模型擬合效果優(yōu)于頻率比例法和Ⅱ型乘法分布法。分析KD模型3種參數(shù)估計(jì)方法擬合優(yōu)度結(jié)果可知，矩法估計(jì)擬合效果最差，極大似然法估計(jì)擬合效果與概率權(quán)重矩法估計(jì)擬合效果相近，但概率權(quán)重矩法估計(jì)擬合效果略優(yōu)于極大似然法估計(jì)擬合效果。

7　結(jié)論

本文應(yīng)用概率權(quán)重矩法定義，嚴(yán)格地推導(dǎo)了KD模型含零值序列頻率計(jì)算參數(shù)估計(jì)公式。以陜西省6個(gè)測站2月份降水資料為例，選取3種頻率計(jì)算模型進(jìn)行含零值降水序列頻率分析，并進(jìn)行擬合優(yōu)度評價(jià)，探討KD模型的適用性，優(yōu)選適用于KD模型的參數(shù)估計(jì)方法，以期為含零值降水序列頻率計(jì)算提供理論依據(jù)。

表4　陜西省6測站2月份降水量的擬合優(yōu)度評價(jià)結(jié)果

（1）選用KD模型模型、頻率比例法、Ⅱ型乘法分布法對含零值降水序列進(jìn)行頻率計(jì)算，探討KD模型適用性，結(jié)果表明，KD模型可以用于含零值降水序列頻率計(jì)算，且擬合效果優(yōu)于頻率比例法及Ⅱ型乘法分布法。該方法具有一定的物理基礎(chǔ)，雖然包含一些復(fù)雜的數(shù)值計(jì)算，但是，借助于計(jì)算機(jī)編程，不難實(shí)現(xiàn)含零值降水序列的頻率計(jì)算。

（2）對比分析矩法、極大似然法、概率權(quán)重矩法估計(jì)KD模型參數(shù)，結(jié)果表明，概率權(quán)重矩法優(yōu)于矩法和極大似然法，文中推導(dǎo)的概率權(quán)重矩法計(jì)算KD模型參數(shù)公式是正確的。

（3）根據(jù)OLS準(zhǔn)則、AIC準(zhǔn)則和RSS最小準(zhǔn)則評價(jià)曲線擬合優(yōu)度，結(jié)果顯示，采用概率權(quán)重矩法參數(shù)估計(jì)的KD模型對含零值降水序列頻率曲線擬合度最優(yōu)。文中模型以期為含零值水文頻率計(jì)算提供一種新的計(jì)算方法。

參考文獻(xiàn)：

［1］黃振平，陳元芳.水文統(tǒng)計(jì)學(xué)［M］.北京：中國水利水電出版社，2011.

［2］馬曉曉，宋松柏.基于Copula函數(shù)的不完全降水序列頻率計(jì)算方法研究［J］.水力發(fā)電學(xué)報(bào)，2017，36（2）：9-17.

［3 ］ZHANG L，SINGH V P.Frequency analysis of flood damage［J］.Journal of Hydrologic Engineering，2005，10（2）：100-109.

［4］JENNINGS M E，BENSON M A.Frequency curves for annual flood series with some zero events or incomplete data［J］.Water Resources Research，1969，5（1）：276-280.

［5］WANG S X，SINGH V P.Frequency estimation for hydrological samples with zero values［J］.Journal of Water Resources Planning and Management，1995，121（1）：98-108.

［6］WEGLARCZYK S，STRUPCZEWSKI W G，SINGH V P.Three-parameter discontinuous distributions for hydrological samples with zero values［J］.Hydrological Processes，2005，19（15）：2899-2914.

［7］宋松柏，康艷.含零值水文序列頻率的計(jì)算原理與應(yīng)用［J］.水文，2012，32（5）：22-26.

［8］STRUPCZEWSKI W G，WEGLARCZYK S，SINGH V P.Impulse response of the kinematic diffusion model as a probability distribution of hydrologic samples with zero values［J］.Journal of Hydrology，2003，270（3/4）：328-351.

［9 ］WOO MingKo，WU Kai.Fitting annual floods with zero-flows［J］.Canadian Water Resources Journal，2013，14（2）：10-16.

［10］楊力行.Ⅱ型乘法分布在干旱資料含零系列統(tǒng)計(jì)分析中的應(yīng)用［J］.水文，1992（3）：21-24.

［11］李滿剛.具有零項(xiàng)系列的中值適線法［J］.山西水利科技，1998（1）：27-29.

［12］毛賽珠.具有零項(xiàng)系列的中值適線法［J］.水文，1985（6）：29-31.

［13］劉志強(qiáng)，王成雄，白霜 .頻率比例法在干旱地區(qū)含零系列分析中的應(yīng)用［J］.東北水利水電，1999（1）：31-35.

［14］水利部長江水利委員會(huì)水文局，水利部南京水文水資源研究所.水利水電工程設(shè)計(jì)洪水計(jì)算手冊［M］.北京：中國水利水電出版社，1995.

［15］楊力行，陳容.Ⅱ型乘法頻率曲線［J］.新疆農(nóng)業(yè)大學(xué)學(xué)報(bào)，1987（1）：48-58.

［16］劉大秀.Ⅱ型乘法頻率曲線的統(tǒng)計(jì)參數(shù)及特征［J］.新疆農(nóng)業(yè)大學(xué)學(xué)報(bào)，1987（1）：58-66.

［17］楊力行，陳容，劉大秀.Ⅱ型乘法頻率曲線應(yīng)用可行性研究［J］.新疆農(nóng)業(yè)大學(xué)學(xué)報(bào)，1988（1）：57-64.

［18］張良，劉鑫楊.含零序列水文頻率計(jì)算方法在張北水文站的應(yīng)用探討［J］.海河水利，2017（1）：36-38.

［19］宋德敦，丁晶.概率權(quán)重矩法及其在P-Ⅲ分布中的應(yīng)用［J］.水利學(xué)報(bào)，1988（3）：1-11.

［20］宋德敦.不連序系列統(tǒng)計(jì)參數(shù)計(jì)算的新方法——概率權(quán)重矩法［J］.水利學(xué)報(bào)，1989（9）：25-32.

［21］秦大庸，孫濟(jì)良.概率權(quán)重矩法在指數(shù)γ分布中的應(yīng)用［J］.水利學(xué)報(bào)，1989（11）：1-9.

［22］STRUPCZEWSKI W G，NAPIORKOWSKI J J，DOOGE J C Z.The distributed Muskingum model［J］.Journal of Hydrology，1989，111（1/4）：235-257.

［23］EINSTEIN H A.Formulas for the transportation of bed load［J］.Transactions of the American Society of Civil Engineers，1947，107（1）：561-597.

［24］EAGLESON P S.Climate，soil，and vegetation：2.The distribution of annual precipitation derived from observed storm sequences［J］.Water Resources Research，1978，14（5）：713-721.

［25］PHIEN H N，NGUYEN V T V.Derivation of the Pearson type（PT）III distribution by using the principle of maximum entropy（POME）［J］.Journal of Hydrology，1987，90（3/4）：351-355.