貝葉斯變點(diǎn)分析模型在水文時(shí)間序列變點(diǎn)分析中的應(yīng)用

貝葉斯變點(diǎn)分析模型在水文時(shí)間序列變點(diǎn)分析中的應(yīng)用

在氣候變化和流域下部人類活動(dòng)的嚴(yán)重影響下，我們必須確保水文時(shí)間序列是否發(fā)生變化。水文時(shí)間序列變異分析主要包括趨勢(shì)分析。跳躍分析是指對(duì)水文時(shí)間序列變化點(diǎn)的評(píng)估和驗(yàn)證。這是水文統(tǒng)計(jì)分析的重要組成部分，在分析和預(yù)測(cè)異常現(xiàn)象方面發(fā)揮著重要作用。目前，在水文中常用的可變點(diǎn)分析方法包括mann-kdandall方法、有序聚類法、最小離散法、最大概率法和貝葉斯推理法。與其他分散分析相比，波葉斯推理法具有優(yōu)點(diǎn)。對(duì)于可能發(fā)生的可變點(diǎn)的概率分布，可以提供概率分布的可能性。因此，在水文序列變量分析中得到廣泛應(yīng)用。隨著貝葉斯變點(diǎn)模型要求時(shí)間序列按照正態(tài)分布，因此首先需要檢查時(shí)間序列的正態(tài)性。對(duì)于不服從正態(tài)分布的時(shí)間序列，可以通過(guò)統(tǒng)計(jì)轉(zhuǎn)換生成符合正態(tài)分布的時(shí)間序列，如豆包c(diǎn)ox轉(zhuǎn)換。此外，考慮到時(shí)間序列觀測(cè)樣本的隨機(jī)性，應(yīng)對(duì)采用k時(shí)間序列特征的最小概率后驗(yàn)概率。為了驗(yàn)證k時(shí)間序列中確定的離散點(diǎn)，有必要根據(jù)最大概率隨時(shí)間序列k確定的變異點(diǎn)進(jìn)行深入的可靠性分析。時(shí)間序列前面、后面的兩個(gè)部分之間是否在時(shí)序列的隨機(jī)波動(dòng)范圍內(nèi)，以及時(shí)間序列中是否存在變形。1方法1.1變點(diǎn)概率分布給定水文時(shí)間序列x,假設(shè)變點(diǎn)的任何可能位置為k.變點(diǎn)將整個(gè)序列分割為兩部分,這兩部分的某些統(tǒng)計(jì)特性(如:均值、方差等)會(huì)有明顯的不同.在變點(diǎn)之前(包括變點(diǎn))的時(shí)間序列稱為前部分,在變點(diǎn)之后(不包括變點(diǎn))的時(shí)間序列稱為后部分.采用如下符號(hào)表示整個(gè)水文序列和前、后兩部分:變點(diǎn)分析的任務(wù)就是檢驗(yàn)變點(diǎn)最有可能發(fā)的時(shí)間位置和所研究的統(tǒng)計(jì)特征值的變化幅度.對(duì)于某一固定時(shí)刻,假設(shè)水文序列的觀測(cè)值的發(fā)生服從某一概率分布.考慮到數(shù)學(xué)推導(dǎo)的復(fù)雜性以及假設(shè)的合理性,這個(gè)分布可以采用正態(tài)分布.如果產(chǎn)生水文序列的物理機(jī)制在某時(shí)刻發(fā)生了突變,那么在該時(shí)刻(即變點(diǎn))的前后兩部分,觀測(cè)值所服從的正態(tài)分布的統(tǒng)計(jì)參數(shù)將會(huì)不再相同.一般假定變點(diǎn)發(fā)生位置的先驗(yàn)分布為均勻分布,在此假定下,變點(diǎn)發(fā)生位置的后驗(yàn)分布密度函數(shù)最終可以簡(jiǎn)化為式中:p(k|X,μa,μb)表示在給定的μa、μb條件下,觀測(cè)資料X={Xk,Xk+1}在k位置發(fā)生變異的后驗(yàn)概率密度函數(shù);p(X|k,μa,μb)表示在給定μa、μb的情況下觀測(cè)資料X={Xk,Xk+1}的聯(lián)合分布函數(shù).μa和μb則假設(shè)它們服從正態(tài)分布.變點(diǎn)發(fā)生位置的后驗(yàn)概率期望值為式(3)可由蒙特卡洛隨機(jī)抽樣法估計(jì).1.2采用正態(tài)法測(cè)試和非正態(tài)分布變換方法Box-Cox冪分布簇取決于使用一個(gè)單個(gè)參數(shù)λ,λ可以通過(guò)大似然法來(lái)估計(jì).1.3時(shí)間序列的相關(guān)性分析由于水文時(shí)間序列的復(fù)雜性和隨機(jī)性,其分析結(jié)果有很大的偶然性.實(shí)際上沒有變異點(diǎn)或者不存在趨勢(shì)的時(shí)間序列由于自相關(guān)性、互相關(guān)性和隨機(jī)波動(dòng)等影響,可能導(dǎo)致分析結(jié)果存在變異點(diǎn)或者存在趨勢(shì),使得分析結(jié)果與實(shí)際情況有偏差.為了估計(jì)時(shí)間序列隨機(jī)性對(duì)分析結(jié)果的影響,一般采用再抽樣方法分析結(jié)果的可靠性[9,10,11,12,13,14,15,16,17,18].基本思想是:首先計(jì)算實(shí)測(cè)時(shí)間序列中用于分析時(shí)間序列變異的統(tǒng)計(jì)量T的值,記為T0;然后對(duì)實(shí)測(cè)時(shí)間序列隨機(jī)再抽樣N次,得到N個(gè)再抽樣樣本,計(jì)算N個(gè)再抽樣樣本中統(tǒng)計(jì)量T的值,記為Ti(1≤i≤N).將Ti值從小到大排序,計(jì)算其不及概率,得到統(tǒng)計(jì)量T的隨機(jī)再抽樣分布.由此可以查得T0值在統(tǒng)計(jì)量T的隨機(jī)再抽樣分布中的顯著性水平α,α表示Ti值大于等于T0值的頻率.若α小于給定的顯著性水平pint,則實(shí)測(cè)時(shí)間序列的T0值在隨機(jī)再抽樣時(shí)間序列的Ti值中是顯著的,即實(shí)測(cè)時(shí)間序列的變異程度相對(duì)于隨機(jī)再抽樣時(shí)間序列的變異程度是顯著的,超出了時(shí)間序列隨機(jī)波動(dòng)范圍,實(shí)測(cè)時(shí)間序列確實(shí)存在變異.水文時(shí)間序列分析中的再抽樣方法主要采用自助法(bootstrap),自助法是一種有返還的抽樣方法.對(duì)于實(shí)測(cè)時(shí)間序列{x1,x2,…,xn},使用隨機(jī)方法,等概率、獨(dú)立、有放回地從集合中抽取n個(gè)數(shù)據(jù),組成一個(gè)新的有順序的時(shí)間序列{x1i,x2i,…,xni}(1≤i≤N),稱為一個(gè)自助樣本,照此重復(fù)抽取得到N個(gè)自助樣本.在自助法中,每個(gè)數(shù)據(jù)可能出現(xiàn)多次,當(dāng)自助法再抽樣1000次以上時(shí),就可以達(dá)到基本穩(wěn)定、一致的結(jié)果.考慮到時(shí)間序列的隨機(jī)性,也需要對(duì)實(shí)測(cè)時(shí)間序列中貝葉斯變點(diǎn)模型確定的變點(diǎn)位置k(對(duì)應(yīng)著最大后驗(yàn)概率pmax的時(shí)間位置)進(jìn)行深入的可靠性分析,特別是針對(duì)最大后驗(yàn)概率pmax相對(duì)較小的時(shí)間序列.如果對(duì)實(shí)測(cè)時(shí)間序列隨機(jī)再抽樣N次得到的N個(gè)自助樣本時(shí)間序列中位置k可能為變點(diǎn)的后驗(yàn)概率p(k)i(1≤i≤N),大部分都能達(dá)到甚至超過(guò)實(shí)測(cè)時(shí)間序列在位置k可能為變點(diǎn)的后驗(yàn)概率pmax,即N個(gè)隨機(jī)再抽樣時(shí)間序列在位置k前后兩部分的差異程度大部分都能達(dá)到甚至超過(guò)實(shí)測(cè)時(shí)間序列位置k前后兩部分的差異程度,那么此時(shí)實(shí)測(cè)時(shí)間序列以最大后驗(yàn)概率pmax對(duì)應(yīng)的位置k作為變點(diǎn)是沒有統(tǒng)計(jì)學(xué)意義的,實(shí)測(cè)時(shí)間序列位置k前后兩部分的差異不顯著,在時(shí)間序列的隨機(jī)波動(dòng)范圍內(nèi).對(duì)貝葉斯變點(diǎn)分析模型中確定的變點(diǎn)(對(duì)應(yīng)著最大后驗(yàn)概率的時(shí)間位置)可靠性分析步驟如下:首先通過(guò)貝葉斯變點(diǎn)模型分析得到實(shí)測(cè)時(shí)間序列{x1,x2,…,xn}的最大后驗(yàn)概率pmax,以及所對(duì)應(yīng)的時(shí)間位置k(1≤k<n);然后對(duì)實(shí)測(cè)時(shí)間序列再抽樣N次,本文中取N=1000,得到N個(gè)自助樣本;計(jì)算每個(gè)自助樣本時(shí)間序列在位置k可能為變點(diǎn)的后驗(yàn)概率p(k)i;對(duì)N個(gè)p(k)i按升序排列,計(jì)為p(k)[r],按下面公式計(jì)算p(k)[r]對(duì)應(yīng)的不及頻率:式中,r代表排序后的秩,由此可以得到隨機(jī)再抽樣時(shí)間序列位置k可能為變點(diǎn)的后驗(yàn)概率分布.據(jù)此可以查得pmax在隨機(jī)再抽樣時(shí)間序列位置k可能為變點(diǎn)的后驗(yàn)概率分布中的顯著性水平α,α表示隨機(jī)再抽樣時(shí)間序列在位置k可能為變點(diǎn)的后驗(yàn)概率p(k)k大于等于pmax的頻率.如果α小于給定的顯著性水平pint,則表示pmax在p(k)k中是顯著的,即實(shí)測(cè)時(shí)間序列位置k前后兩部分的差異程度相對(duì)于隨機(jī)再抽樣時(shí)間序列中位置k前后兩部分的差異程度是顯著的.實(shí)測(cè)時(shí)間序列中位置k前后兩部分的差異程度超出了時(shí)間序列的隨機(jī)波動(dòng)范圍,k位置確實(shí)發(fā)生了變異,因此以位置k為變點(diǎn)是可靠的.2使用示例2.1年度流量最大化時(shí)間序列本文所分析的某水文站1942—2008年年徑流最大值時(shí)間序列曲線如圖1所示.其中年徑流最大值均值Qmax=10485m3/s.2.2box-cox轉(zhuǎn)換分析考慮到貝葉斯變點(diǎn)模型預(yù)測(cè)變異點(diǎn)要求時(shí)間序列服從正態(tài)分布,因此檢驗(yàn)了實(shí)測(cè)年徑流最大值時(shí)間序列的正態(tài)性,正態(tài)概率圖如圖2所示.由圖2可以看出實(shí)測(cè)年徑流最大值時(shí)間序列的散點(diǎn)圖與直線有一定的偏差,為了進(jìn)一步檢驗(yàn)時(shí)間序列的正態(tài)性,采用W檢驗(yàn)和D檢驗(yàn),結(jié)果如表1所示.由表1可以看出針對(duì)實(shí)測(cè)年徑流最大值時(shí)間序列,W檢驗(yàn)和D檢驗(yàn)的顯著性水平都遠(yuǎn)小于5%,因此在5%的顯著性水平下拒絕原假設(shè),認(rèn)為實(shí)測(cè)年徑流最大值序列不服從正態(tài)分布.因此對(duì)實(shí)測(cè)年徑流最大值時(shí)間序列進(jìn)行Box-Cox轉(zhuǎn)換,使轉(zhuǎn)換后的時(shí)間序列服從正態(tài)分布,根據(jù)殘差最小平方和準(zhǔn)則計(jì)算得到λ=0.342.將λ值代入公式(6),得到實(shí)測(cè)年徑流最大值時(shí)間序列Box-Cox轉(zhuǎn)換后時(shí)間序列(以下簡(jiǎn)稱Box-Cox轉(zhuǎn)換后時(shí)間序列),其正態(tài)概率圖如圖3所示.由圖3可以看出Box-Cox轉(zhuǎn)換后時(shí)間序列的散點(diǎn)圖基本上接近直線,為了進(jìn)一步檢驗(yàn)序列的正態(tài)性,采用W檢驗(yàn)和D檢驗(yàn),結(jié)果如表2所示.由表2可以看出針對(duì)Box-Cox轉(zhuǎn)換后時(shí)間序列,W檢驗(yàn)和D檢驗(yàn)的顯著性水平均大于5%,因此在5%的顯著性水平下不拒絕原假設(shè),認(rèn)為BoxCox轉(zhuǎn)換后時(shí)間序列服從正態(tài)分布.2.3變點(diǎn)分析計(jì)算針對(duì)實(shí)測(cè)年徑流最大值時(shí)間序列Box-Cox轉(zhuǎn)換后時(shí)間序列,建立貝葉斯變點(diǎn)分析模型分析實(shí)測(cè)年徑流最大值時(shí)間序列的跳躍點(diǎn).因?yàn)锽ox-Cox轉(zhuǎn)換僅改變時(shí)間序列的參數(shù)分布情況,并不會(huì)影響時(shí)間序列的變異性,故可以將Box-Cox轉(zhuǎn)換后時(shí)間序列的變點(diǎn)分析結(jié)果對(duì)應(yīng)到實(shí)測(cè)時(shí)間序列的變點(diǎn)分析.經(jīng)過(guò)分析計(jì)算,變點(diǎn)發(fā)生位置的先驗(yàn)和后驗(yàn)概率分布如圖4所示.由圖4可知,k=1989所對(duì)應(yīng)的后驗(yàn)概率最大,為0.35(pmax=0.35).1942—1989年,實(shí)測(cè)年徑流最大值時(shí)間序列的均值為13272m3/s,而2001—2008年實(shí)測(cè)年徑流最大值時(shí)間序列的均值為3445m3/s,均值比1989年前下降了9827m3/s,降幅74%,在5%的顯著性水平下年徑流最大值流量序列的均值變化(減少)是顯著的.經(jīng)計(jì)算,該結(jié)論與兩類有序聚類分析法分析年徑流最大值時(shí)間序列跳躍點(diǎn)的結(jié)果一致.3隨機(jī)再抽樣時(shí)間序列及相關(guān)概率分布為了檢驗(yàn)分析結(jié)果的可靠性,采用自助法再抽樣方法,對(duì)實(shí)測(cè)年徑流最大值時(shí)間序列Box-Cox轉(zhuǎn)換后時(shí)間序列進(jìn)行再抽樣1000次(N=1000).計(jì)算1000個(gè)自助樣本中1989年可能為變點(diǎn)的后驗(yàn)概率p(1989)i(1≤i≤1000),將p(1989)i按升序排列,采用公式(7)計(jì)算不及概率,得到隨機(jī)再抽樣時(shí)間序列1989年可能為變點(diǎn)的后驗(yàn)概率分布曲線,如圖5所示.從圖5可以看出,后驗(yàn)概率0.12所對(duì)應(yīng)的顯著性水平α幾乎接近0;因此得出實(shí)測(cè)年徑流最大值時(shí)間序列中1989年可能為變點(diǎn)的后驗(yàn)概率,0.35在隨機(jī)再抽樣時(shí)間序列1989年可能為變點(diǎn)的后驗(yàn)概率分布中的顯著性水平幾乎為0,小于顯著性水平pint=5%,即實(shí)測(cè)年徑流最大值時(shí)間序列1989年前后兩部分的差異程度相對(duì)于隨機(jī)再抽樣時(shí)間序列中的1989年前后兩部分的差異程度是顯著的,超出了時(shí)間序列的隨機(jī)波動(dòng)范圍,可見1989年確實(shí)發(fā)生了變異,因此實(shí)測(cè)年徑流最大值時(shí)間序列中以1989年作為變點(diǎn)是可靠的.4時(shí)間序列隨機(jī)性檢驗(yàn)與正態(tài)分布本文通過(guò)對(duì)實(shí)際觀測(cè)時(shí)間序列數(shù)據(jù)處理,建立了用于時(shí)間序列變點(diǎn)分析的貝葉斯模型,推求得到了變點(diǎn)位置的后驗(yàn)概率分布,確定變點(diǎn)的最可能位置,并對(duì)結(jié)果進(jìn)行了可靠性檢驗(yàn),本例中以1989年作為變點(diǎn)是可靠的.考慮到時(shí)間序列的隨機(jī)性,建議在今后工作中采用貝葉斯變點(diǎn)分析模型確定變異點(diǎn)時(shí),采用自助再抽樣方法對(duì)結(jié)果進(jìn)行可靠性分析.由于上述貝葉斯變點(diǎn)分析模型中假設(shè)