， ，，

(1.南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國家重點實驗室, 南京 210029； 2.中國科學(xué)院南京地理與湖泊研究所 中國科學(xué)院流域地理學(xué)重點實驗室，南京 210008)

基于核密度估計的AM-MCMC算法在徑流模擬中的應(yīng)用

童坤1，2，劉恒1，耿雷華1，徐澎波1

(1.南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國家重點實驗室, 南京 210029； 2.中國科學(xué)院南京地理與湖泊研究所 中國科學(xué)院流域地理學(xué)重點實驗室，南京 210008)

無資料或資料稀缺地區(qū)的徑流概率模擬， 是目前水文研究難點問題之一。 基于此, 利用Kernal核密度估計法估算出流量的月徑流概率密度函數(shù)， 采用基于自適應(yīng)采樣算法(Adaptive Metropolis algorithm， AM)的馬爾可夫鏈蒙特卡羅(Markov Chain Monte Carlo,MCMC)模擬方法求解， 最后給出月徑流量的模擬預(yù)測。 實例表明基于Kernel核密度估計的AM-MCMC算法模型計算結(jié)果精度較高， 有良好的應(yīng)用價值， 可在資料較少地區(qū)推廣使用。

徑流模擬；概率分布；核密度估計； AM-MCMC算法；羅岙水庫

1 研究背景

水文變量概率密度函數(shù)復(fù)雜且未知，某一指定概率分布與真實分布存在著差異。目前對于水文概率密度估計主要有參數(shù)估計和非參數(shù)估計2種類型：參數(shù)估計法是密度函數(shù)結(jié)構(gòu)，此時的密度估計就是傳統(tǒng)的參數(shù)估計問題(如P-Ⅲ型分布已知求cv,cs)；非參數(shù)估計是密度函數(shù)未知，僅從既有的樣本出發(fā)得出密度函數(shù)的表達式。

非參數(shù)估計方法主要有核密度法、直方圖法、樣條函數(shù)法以及混合概率密度法等，其中核密度估計是較常用和有效的方法。本文提出基于Kernel核密度估計的AM算法模型，通過Kernel估計法，估計出壩址月徑流量的密度分布，然后利用AM算法對其進行模擬。使用核函數(shù)方法進行初始概率密度估計，推導(dǎo)出條件概率密度估計，以分布權(quán)重為抽樣概率計算出均值、方差，并由此計算出模擬值。基于AM算法的馬爾可夫鏈蒙特卡羅(MCMC)模擬不再依賴于參數(shù)的推薦分布也不受參數(shù)的先驗區(qū)間的影響,能夠保證算法的遍歷性和較高的抽樣效率[1]。與過去單純的馬爾可夫鏈蒙特卡羅方法用于計算月徑流量過程相比，該方法精度更高，可用于進一步推廣應(yīng)用。

2 核密度估計法

核密度估計法(Kernel density estimation)又名Parzen窗(Parzen window)，由Rosenblatt[2](1956)和Parzen[3](1962)提出，在概率論中被用來估計未知隨機變量的密度函數(shù)。核密度估計的概念是直方圖概念的推廣。Scott等[4](1996)證明當(dāng)平均平移直方圖(Averaged Shifted Histogram,ASH)的直方數(shù)量趨于無窮時，直方圖便趨于Kernel核密度。1994年，Ruppert等[5]利用數(shù)據(jù)集密度函數(shù)聚類算法對其進行修正。1996年我國學(xué)者崔恒建等[6]引入核密度估計法對直徑分布進行擬合，2000年開始應(yīng)用到我國各個研究領(lǐng)域[7-9];唐林俊等[10](2006)在單變量核密度估計的基礎(chǔ)上建立了風(fēng)險價值預(yù)測的預(yù)測模型；王文圣等[11](2001)首次將單變量核密度估計用于徑流隨機模擬。核密度估計法由于其從樣本出發(fā)研究現(xiàn)有數(shù)據(jù)的分布特征而受到高度重視和應(yīng)用。相對于其他方法，核密度估計法主要優(yōu)點就是對數(shù)據(jù)分布不附加任何假定，不需要了解數(shù)據(jù)分布的先驗知識。

設(shè)X1,X2,…，Xn是密度函數(shù)為f(x)的未知總體的獨立同分布的隨機變量，Kernel核密度估計為

(1)

(2)

多元函數(shù)核密度估計的是通過將一維核函數(shù)連乘得到乘積核函數(shù)，利用乘積核函數(shù)進行估計。乘積核函數(shù)表達式為

(3)

式中:xij是第j個分量的第i個觀測值;hd為維度d下的窗寬。

本文只涉及到一個變量，因此使用式(1)進行計算。

從Kernel密度估計的表達式可以看出，該方法是通過將核函數(shù)K分別置于以樣本Xi為中心處，然后計算此時所有點的函數(shù)值，最后對其進行平均。常用的K(t)有正態(tài)函數(shù)、三角核函數(shù)、雙權(quán)重核函數(shù)、三權(quán)重核函數(shù)、Epanechnikov核函數(shù)等，如表1所示。

表1 典型核函數(shù)表達式Table 1 Typical expressions of kernel function

3 窗寬的確定(LSCV法)

目前，有關(guān)窗寬選擇的方法可分為以下幾種：交錯鑒定方法、懲罰函數(shù)法、插入法以及對比方法[12-15]。常用的是插入法，即把未知函數(shù)的估計插入到漸近公式里以選擇最佳窗寬。本文采用固定窗寬法，固定窗寬就是在每一個擬合點取等窗寬，基于最小平方差(LSCV)的思想，根據(jù)積分均方誤差(MISE)最小，求出最優(yōu)窗寬。MISE為

(4)

(5)

(6)

(7)

(8)

因此，如需MISE最小，則AMISE達到最小即可，對AMISE求一階導(dǎo)數(shù)并令其等于0，求得最優(yōu)窗寬h*為

(9)

4 AM-MCMC抽樣算法

4.1 馬爾可夫鏈蒙特卡羅方法(MCMC)

馬爾可夫鏈蒙特卡羅方法(MCMC)通過構(gòu)造馬爾可夫鏈，使得計算出的馬爾可夫鏈穩(wěn)定分布為要求的目標(biāo)分布，也就是說希望通過利用馬爾可夫鏈取得的樣本就可以直接當(dāng)成是目標(biāo)分布中產(chǎn)生的樣本來利用，其核心思想是建立轉(zhuǎn)換函數(shù)使得無論初始值取何值，最后馬爾可夫鏈都會收斂到目標(biāo)分布。馬爾可夫鏈蒙特卡羅方法作為一種隨機模擬方法，早已被應(yīng)用到物理、天文、氣象、通信等各個領(lǐng)域，MCMC的關(guān)鍵是如何選擇推薦分布轉(zhuǎn)移密度使采樣更加有效，常用的方法有Gibbs采樣[16]、Metropolis-Hasting取樣[17]和自適應(yīng)取樣[18](Adaptive metropolis-hasting method，AM算法)。這3種抽樣方法中，只有AM算法不依賴于事先的推薦分布且收斂速度較快，因此本文采用AM算法。

4.2 收斂判斷準(zhǔn)則

(10)

其中

(11)

(12)

5 誤差判別

選取相對誤差與均方差作為計算結(jié)果的判別函數(shù)，其中，相對誤差公式為

(13)

均方根誤差RMS(Root Mean Squared Residual)公式為

(14)

因為RMS計算公式中沒有考慮擬合來水量變化幅度對模型精度的影響，因此，引入另外一個更加準(zhǔn)確的判別參數(shù):標(biāo)準(zhǔn)化殘差均方根Normalized RMS，其計算公式為

(15)

式中:(Xobs)max為觀測值中的極大值;(Xobs)min為觀測值中的極小值。

6 案例應(yīng)用

羅岙水庫位于健跳港上游，見圖1。壩址以上集水面積12.3 km2，多年平均入庫水量約1 100萬m3。該水庫為1969年動工，1979年完工的以灌溉為主，兼顧防洪、供水的小型水庫，由六敖鎮(zhèn)鎮(zhèn)政府管理。2005年，該水庫被列為浙江省“千庫保安”建設(shè)計劃和臺州市重點工程，開展加高加固。羅岙水庫工程樞紐建筑物由攔河壩、溢洪道、輸水隧洞等組成。除險加固以前，水庫總庫容478萬 m3，正常蓄水位6.9 m，相應(yīng)庫容339萬 m3，死水位1.50 m，相應(yīng)庫容12萬 m3，興利庫容327萬 m3。

圖1 羅岙水庫及雨量、流量站分布示意圖Fig.1 Map of Luo’ao reservoir and distribution of rainfall and discharge gauging stations

根據(jù)海游站1952—2006年共55 a逐月降水量與蒸發(fā)量資料，對羅岙水庫壩址長系列來水月徑流量進行模擬分析。壩址長系列歷年逐月平均徑流量見圖2。

圖2 羅岙水庫壩址逐月平均徑流量變化過程曲線Fig.2 Monthly average flow at the dam site of Luo’ao reservoir

對1956—2006年平均月徑流量進行Kernel估計，假設(shè)預(yù)先給定的概率密度分布是正態(tài)分布，則核函數(shù)K(t)選擇正態(tài)分布函數(shù)，在每個樣本處利用核函數(shù)進行該處的密度估計，然后樣本間進行差值得到連續(xù)的估計密度曲線。圖3是抽樣樣本均值與方差軌線圖，可以看出所取隨機數(shù)4 000個以后基本趨于穩(wěn)定。

圖3 均值方差軌跡線Fig.3 Trajectory of mean value and variance

利用得出的密度分布，結(jié)合AM算法對其進行抽樣，設(shè)定抽樣次數(shù)50 000次。根據(jù)LSCV法，通過多次擬合計算，最終確定窗寬h=0.01，頻率直方圖及理論概率密度曲線如圖4，可見樣本的頻率直方圖與理論密度曲線擬合很好。

圖4 概率密度分布Fig.4 Distribution of probability density

圖5 殘差分布Fig.5 Distribution of residual

對水庫壩址的實際月徑流量進行模擬分析，如圖6所示，根據(jù)模擬結(jié)果計算得出RMS=0.48%,Normalized RMS=0.002%，說明所采集的樣本合理，可以用來對水庫徑流量進行模擬計算。

圖6 壩址月徑流量擬合曲線Fig.6 Curve fitting of monthly runoff at dam site

7 結(jié) 語

本文推薦的核密度估計法不需要事先了解模型的一系列參數(shù)，且基于AM算法的MCMC抽樣模型無需預(yù)先給定推薦分布，不依賴于過去的經(jīng)驗，遍歷性好、適應(yīng)性強，該方法對于資料較少的地區(qū)可以進行很好的應(yīng)用。從案例的實際應(yīng)用中可以看出，基于核密度估計的AM-MCMC算法具有較高的計算精度，利用該模型計算出來的長系列水庫月徑流量與實測值差異不大，標(biāo)準(zhǔn)化殘差均方根僅為0.002%，可以用于實際應(yīng)用中估計水庫的水文風(fēng)險。不過本文對于核密度估計中的窗寬采用固定窗寬的辦法，如何選擇最優(yōu)窗寬優(yōu)化核密度估計，提高模擬的精度是下一步要研究的對象。

[1] 邢貞相, 芮孝芳, 崔海燕,等. 基于AM-MCMC算法的貝葉斯概率洪水預(yù)報模型[J]. 水利學(xué)報, 2007, 38(12):1500-1506.

[2] ROSENBLATT M. Remarks on Some Nonparametric Estimates of a Density Function[J]. Annals of Mathematical Statistics, 1956, 27(3): 832-837.

[3] PARZEN E. On Estimation of a Probability Density Function and Mode[J]. Annals of Mathematical Statistics, 1962, 33(3): 1065-1076.

[4] SCOTT D W, WHITTAKER G. Multivariate Applications of the ASH in Regression[J]. Communications in Statistics—Theory and Methods, 1996, 25(11): 2521-2530.

[5] RUPPERT D, CLINE D B. Bias Reduction in Kernel Density Estimation by Smoothed Empirical Transformations[J]. The Annals of Statistics, 1994:185-210, doi: 10.1214/aos/1176325365.

[6] 崔恒建, 王雪峰. 核密度估計及其在直徑分布研究中的應(yīng)用[J]. 北京林業(yè)大學(xué)學(xué)報, 1996,18(2):67-72.

[7] 劉 銳, 胡偉平, 王紅亮, 等. 基于核密度估計的廣佛都市區(qū)路網(wǎng)演變分析[J]. 地理科學(xué), 2011,31(1):81-86.

[8] 趙 淵, 沈智健, 周念成, 等. 基于序貫仿真和非參數(shù)核密度估計的大電網(wǎng)可靠性評估[J]. 電力系統(tǒng)自動化, 2008,32(6):14-19.

[9] 凌建國, 劉爾琦, 梁海燕, 等. 基于核密度估計的紅外目標(biāo)提取方法[J]. 紅外與毫米波學(xué)報, 2006,25(6):434-438.

[10] 唐林俊, 楊 虎, 張洪陽. 核密度估計在預(yù)測風(fēng)險價值中的應(yīng)用[J]. 數(shù)學(xué)的實踐與認(rèn)識, 2006,35(10):29-35.

[11] 王文圣, 丁 晶. 單變量核密度估計模型及其在徑流隨機模擬中的應(yīng)用[J]. 水科學(xué)進展, 2001,12(3):367-372.

[12] JONES M C, MARRON J S, SHEATHER S J. A Brief Survey of Bandwidth Selection for Density Estimation[J]. Journal of the American Statistical Association, 1996, 91(433): 401-407.

[13] SAIN S R, BAGGERLY K A, SCOTT D W. Cross-validation of Multivariate Densities[J]. Journal of the American Statistical Association, 1994, 89(427): 807-817.

[14] BOLANCé C, GUILLEN M, NIELSEN J P. Kernel Density Estimation of Actuarial Loss Functions[J]. Insurance: Mathematics and Economics, 2003, 32(1): 19-36.

[15] FUKUNAGA K, HOSTETLER L. The Estimation of the Gradient of a Density Function, with Applications in Pattern Recognition[J]. IEEE Transactions on Information Theory, 1975, 21(1): 32-40.

[16] SMITH A F, ROBERTS G O. Bayesian Computation via the Gibbs Sampler and Related Markov Chain Monte Carlo Methods[J]. Journal of the Royal Statistical Society. Series B (Methodological), 1993:3-23, doi: 10.2307/2346063.

[17] CHIB S, GREENBERG E. Understanding the Metropolis-Hastings Algorithm[J]. The American Statistician, 1995, 49(4): 327-335.

[18] HAARIO H, SAKSMAN E, TAMMINEN J. An Adaptive Metropolis Algorithm[J]. Bernoulli, 2001, 7(2): 223-242.

[19] GELMAN A, RUBIN D B. Inference from Iterative Simulation Using Multiple Sequences[J]. Statistical Science, 1992, 7(4): 457-472.

AM-MCMC Algorithm for Runoff Simulation ModelBased on Kernel Density Estimation

TONG Kun1,2, LIU Heng1, GENG Lei-hua1, XU Peng-bo1

(1.State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Nanjing Hydraulic Research Institute, Nanjing 210029, China; 2.Key Laboratory of Watershed Geographic Science, Nanjing Institute of Geography and Limnology, Chinese Academy of Sciences, Nanjing 210008, China)

The simulation of runoff probability in an area in lack of runoff data is a difficulty in hydrological research. In this article, we try to establish the probability density function of monthly runoff flow by adopting kernal density estimation method, and give the solution by Markov Chain Monte Carlo (MCMC) simulation method based on Adaptive Metropolis (AM) algorithm. Case study shows that the AM-MCMC algorithm model based on kernel density estimation is of high accuracy and good application value. It can be used in areas in lack of data.

runoff simulation; probability distribution; kernel density estimation; AM-MCMC algorithm; Luo’ao Reservoir

2016-08-17;

2016-09-29

南京水利科學(xué)研究院院基金項目(Y516011);水利部公益性項目(201201020)

童 坤(1986-)，女，江蘇高郵人，工程師，博士，主要從事水資源配置方面的研究工作。E-mail:tongkun0502@hotmail.com

10.11988/ckyyb.20160843

TV214

A

1001-5485(2018)01-0036-04

(編輯：王 慰)