耦合馬爾科夫鏈和伽馬分布函數(shù)的降水模擬模型及其應(yīng)用

付 強，王 斌

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

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

1 逐日降水的隨機模擬模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

應(yīng)用哈爾濱48a的降水資料，統(tǒng)計出哈爾濱各月的P （W／D）、P （W／W）及f見圖1。由圖1可見，P （W／D）、P （W／W）與f的月份變化趨勢基本一致。統(tǒng)計分析結(jié)果表明，哈爾濱的P（W／D）、P （W／W）與f之間均存在很強的相關(guān)性，且P （W／D）與f之間的相關(guā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實測降水資料，求得哈爾濱12個月份的α、β及p見圖2。由圖2可見，α與p間的月份變化趨勢差別很大，而β與p各月的變化趨勢基本一致。統(tǒng)計分析結(jié)果表明，哈爾濱的α與p之間的線性關(guān)系很差，但β與p之間具有很強的線性關(guān)系，這也與文獻 ［2］的研究結(jié)論相同。β與p的回歸方程見式 （11）：

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

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

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

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

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

3 結(jié) 論

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

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

［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］王 斌，付 強，王 敏，等．幾種模擬逐日降水的分布函數(shù)比較分析 ［J］．數(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］王 斌，付 強，張金萍，等．逐日降水的隨機模擬模型及其在黑龍江省的應(yīng)用 ［J］．農(nóng)機化研究，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．