付 強(qiáng)，王 斌

（東北農(nóng)業(yè)大學(xué) 水利與建筑學(xué)院，哈爾濱 150030）

由于降水現(xiàn)象具有隨機(jī)性特點(diǎn)，在水文水利計(jì)算等科研實(shí)踐工作中，一般應(yīng)用數(shù)理統(tǒng)計(jì)方法研究降水現(xiàn)象，應(yīng)用這類方法的前提是降水資料的樣本系列較長(zhǎng)。然而，我國(guó)大部分地區(qū)難以收集到長(zhǎng)系列的降水?dāng)?shù)據(jù)，其資料系列一般只有幾十年，且多以逐日數(shù)據(jù)為主。近些年來，各國(guó)學(xué)者均對(duì)降水的模型化和模擬化產(chǎn)生了濃厚的興趣，很通用的一種方法是采用馬爾科夫鏈和某種分布函數(shù)相結(jié)合的模型模擬逐日的降水過程，即用馬爾科夫鏈描述降水日的發(fā)生，再用選定的某種分布函數(shù)產(chǎn)生降水日的降水量［1－4］。本文應(yīng)用哈爾濱降水資料驗(yàn)證這種延展降水系列的方法，以期為評(píng)估我國(guó)高寒區(qū)降水對(duì)農(nóng)業(yè)的影響、旱作區(qū)雨水的高效利用、制定合理的灌溉制度等方面提供科學(xué)依據(jù)。

1 逐日降水的隨機(jī)模擬模型

1．1 馬爾科夫鏈的兩狀態(tài)轉(zhuǎn)移概率

大量研究表明，在一定的環(huán)境條件下，用一階馬爾科夫鏈描述降水的發(fā)生是簡(jiǎn)單有效的，這種模型不僅能夠保持降水的很多重要特性，并且適用性也很好［5－8］。一階馬爾科夫鏈的兩狀態(tài)分別是降水日 （wet day）和非降水日 （dry day），狀態(tài)的轉(zhuǎn)換可以用兩個(gè)轉(zhuǎn)移概率來描述，即由降水日到降水日的轉(zhuǎn)移概率P （W／W）和由非降水日到降水日的轉(zhuǎn)移概率P （W／D）。Shu Geng等對(duì)7個(gè)站點(diǎn)降水?dāng)?shù)據(jù)進(jìn)行了統(tǒng)計(jì)分析，結(jié)果表明，在每個(gè)站點(diǎn)，各月份P （W／D）、P （W／W）均與同月降水日出現(xiàn)的頻率f間存在線性關(guān)系，尤其P （W／D）與f之間具有很強(qiáng)的相關(guān)性，可用下式表示［2］：

式中f為某月降水日出現(xiàn)頻率的多年平均值，可根據(jù)實(shí)測(cè)降水資料推求；a、b為回歸系數(shù)。

Shu Geng等發(fā)現(xiàn)利用式 （1）可以解釋各站點(diǎn)不同時(shí)間轉(zhuǎn)移概率總變異的96．5%，因此建議采用式 （3）、式 （4）推求轉(zhuǎn)移概率［2］：

1．2 伽瑪分布及其參數(shù)

如果連續(xù)型隨機(jī)變量X的密度函數(shù)為式 （5），則稱X服從伽瑪 （Gamma）分布：

式中α為伽瑪分布的形狀參數(shù)；β為伽瑪分布的尺度參數(shù)。

設(shè)X為表示日降水量的隨機(jī)變量，當(dāng)X服從伽瑪分布時(shí)，根據(jù)伽瑪分布函數(shù)的統(tǒng)計(jì)特性，參數(shù)α和β也可根據(jù)實(shí)測(cè)降水資料求得。Shu Geng等在分析各站點(diǎn)降水?dāng)?shù)據(jù)時(shí)還發(fā)現(xiàn)，參數(shù)β與多年平均各月降水日的降水量之間也存在很強(qiáng)的線性關(guān)系，應(yīng)用7個(gè)站點(diǎn)的平均結(jié)果得到了式 （6）、式 （7），在不同地區(qū)和不同的月份間，利用該式也可以解釋96．5%總變異［2］：

式中p為某月降水日降水量的多年平均值。

1．3 逐日降水序列的產(chǎn)生

當(dāng)轉(zhuǎn)移概率和伽瑪分布參數(shù)確定以后，利用計(jì)算機(jī)產(chǎn)生 ［0，1］區(qū)間上均勻分布的隨機(jī)數(shù)，并將這個(gè)隨機(jī)數(shù)與P （W／D）和P （W／W）相比較，估計(jì)該日是否產(chǎn)生降水。對(duì)降水日降水量的模擬就是求出指定P所對(duì)應(yīng)的隨機(jī)變量X的取值xP，即求出的xP應(yīng)滿足F （X＞xP）＝P，亦即：

當(dāng)為降水日時(shí)，利用式 （8）通過數(shù)值積分即可求得降水日的降水量。

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

利用f和p估計(jì)轉(zhuǎn)移概率和伽瑪分布參數(shù)，不僅簡(jiǎn)化了逐日降水模擬模型，而且擴(kuò)大了該模型的適用范圍。本文應(yīng)用哈爾濱1961～2008年的逐日降水?dāng)?shù)據(jù)，檢驗(yàn)該方法在我國(guó)高寒地區(qū)應(yīng)用的有效性。由于直接積分式 （8）非常繁雜，因此應(yīng)用Matlab軟件的伽瑪分布函數(shù) （GAMRND）實(shí)現(xiàn)式（8）的計(jì)算過程，即將α、β兩個(gè)參數(shù)提供給GAMRND函數(shù)，令該函數(shù)產(chǎn)生服從伽瑪分布的隨機(jī)數(shù)列［3－4］。

2．1 轉(zhuǎn)移概率和伽瑪分布參數(shù)的確定

應(yīng)用哈爾濱48a的降水資料，統(tǒng)計(jì)出哈爾濱各月的P （W／D）、P （W／W）及f見圖1。由圖1可見，P （W／D）、P （W／W）與f的月份變化趨勢(shì)基本一致。統(tǒng)計(jì)分析結(jié)果表明，哈爾濱的P（W／D）、P （W／W）與f之間均存在很強(qiáng)的相關(guān)性，且P （W／D）與f之間的相關(guān)性更強(qiáng)，這與文獻(xiàn) ［2］的研究結(jié)論相同，P （W／D）、P （W／W）與f的回歸方程見式 （9）、式 （10）：

圖1 哈爾濱降水的轉(zhuǎn)移概率和降水日發(fā)生頻率Fig．1 Precipitation transition probability and occurrence frequency per day in Harbin

根據(jù)哈爾濱48a實(shí)測(cè)降水資料，求得哈爾濱12個(gè)月份的α、β及p見圖2。由圖2可見，α與p間的月份變化趨勢(shì)差別很大，而β與p各月的變化趨勢(shì)基本一致。統(tǒng)計(jì)分析結(jié)果表明，哈爾濱的α與p之間的線性關(guān)系很差，但β與p之間具有很強(qiáng)的線性關(guān)系，這也與文獻(xiàn) ［2］的研究結(jié)論相同。β與p的回歸方程見式 （11）：

對(duì)比式 （9）、式 （10）與式 （3）、式 （4），以及式 （11）與式 （6）可以看出，應(yīng)用哈爾濱實(shí)測(cè)降水?dāng)?shù)據(jù)求得的回歸模型與文獻(xiàn) ［2］中建議的經(jīng)驗(yàn)公式存在較大差別，但與文獻(xiàn) ［4］的研究結(jié)果相近。因此，即使在外國(guó)能適應(yīng)很大環(huán)境變化范圍的模型，卻不一定適合我國(guó)的部分地區(qū)。

圖2 哈爾濱降水的伽瑪分布參數(shù)和降水日降水量Fig．2 Precipitation Gamma distribution parameters and daily precipitation in Harbin

2．2 逐日降水過程模擬及分析

以50、100、200、500、1 000a為模擬年限 （對(duì)于每種年限的模擬，均假設(shè)從一個(gè)非降水日開始），采取兩種方案模擬哈爾濱的逐日降水過程，兩種方案依次為：①統(tǒng)計(jì)分析實(shí)測(cè)降水?dāng)?shù)據(jù)直接確定P （W／D）、P （W／W）、α、β；②根據(jù)統(tǒng)計(jì)實(shí)測(cè)降水?dāng)?shù)據(jù)得到的f、p，再利用式 （9）、式 （10）、式 （11）、式（7）間接推求出P （W／D）、P （W／W）、α、β。模擬的月份平均結(jié)果見表1、表2，精度分析 （各月份的模擬數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)的相關(guān)系數(shù)r）見表3。

表1 哈爾濱降水量模擬值Table 1 Precipitation simulation values in Harbin ／mm

表2 哈爾濱降水天數(shù)模擬值Table 2 Simulation value of precipitation day in Harbin ／d

表3 模擬值與觀測(cè)值的相關(guān)系數(shù)Table 3 The correlation coefficient between simulation values and observed values

由表1～表3可見：對(duì)于不同模擬年限的降水系列，模型模擬的精度很高，不同模擬年限間的模擬結(jié)果間并無顯著差別。因此，根據(jù)實(shí)際需要，可以應(yīng)用上述方法，選擇一定的模擬年限模擬哈爾濱的逐日降水過程。對(duì)于哈爾濱周邊地區(qū)，當(dāng)降水資料不全時(shí)，如果能收集到該地區(qū)每月降水天數(shù)和降水量的多年平均值，即可確定f、p，利用式 （9）、式 （10）、式 （11）、式 （7）等即可模擬該地區(qū)的逐日降水過程。

3 結(jié) 論

哈爾濱的P （W／D）和P （W／W）與f之間，以及β與p之間均存在很強(qiáng)的線性關(guān)系，這與國(guó)外的研究結(jié)論相同，但建立的回歸模型系數(shù)與國(guó)外地區(qū)的經(jīng)驗(yàn)公式差別較大，因此，在引進(jìn)國(guó)外模型模擬逐日降水過程時(shí)，模型及其參數(shù)的有效性需要實(shí)地?cái)?shù)據(jù)驗(yàn)證。

利用哈爾濱實(shí)測(cè)降水資料確定P （W／D）、P（W／W）、α、β后，應(yīng)用一階馬爾科夫鏈和伽瑪分布函數(shù)相結(jié)合的隨機(jī)模型可以模擬哈爾濱的逐日降水過程，模型的模擬精度很高；本文建立的3個(gè)回歸模型的R2＞0．97，應(yīng)用此3個(gè)模型，也可為資料不全的哈爾濱周邊地區(qū)模擬逐日降水過程時(shí)提供參考。

［1］Stern R．D．，Coe R．The use of rainfall models in agricultural planning ［J］．Agric．Meteorogy，1982，26：35－50．

［2］Shu Geng，F(xiàn)rits W．T．Penning de Vries，Iwan Supit．A simple method for generating daily rainfall data ［J］．Agricultural and Forest Meteorology，1986，36：363－376．

［3］王 斌，付 強(qiáng)，王 敏，等．幾種模擬逐日降水的分布函數(shù)比較分析 ［J］．?dāng)?shù)學(xué)的實(shí)踐與認(rèn)識(shí)，2011，41 （9）：128－133．Wang Bin，F(xiàn)u Qiang，Wang Min，et al．Comparative analysis on three distribution functions simulating daily rainfall ［J］．Mathematics in Practice and Theory，2011，41 （9）：128－133．

［4］王 斌，付 強(qiáng)，張金萍，等．逐日降水的隨機(jī)模擬模型及其在黑龍江省的應(yīng)用 ［J］．農(nóng)機(jī)化研究，2011，（9）：65－69．Wang Bin，F(xiàn)u Qiang，Zhang Jinping，et al．Stochastic model of generating daily rainfall data and its application in Heilongjiang Province ［J］．Journal of Agricultural Mechanization Research，2011，（9）：65－69．

［5］Ahmed J．，Van Bavel，C．H．M．，et al．Optimization of crop irrigation strategy under a stochastic weather regime：a simulation study ［J］．Water Resour．Res．，1976，12：1 241－1 247．

［6］Delleur J．W．，Kavvas M．L．Stochastic models for monthly rainfall forecasting and synthetic generation［J］．J．AppL MeteoroL，1978，17：1 528－1 536．

［7］Nicks A．D．，Harp J．F．Stochastic generation of temperature and solar radiation data ［J］．J．Hydrol．，1980，48：1－17．

［8］Bruhn J．A．，F(xiàn)ry W．E．，F(xiàn)ick G．W．Simulation of daily weather data using theoretical probability distributions［J］．J．Appl．Meteorol．，1980，19：1 029－1 036．