王大洋,杜 懿,王大剛,龍鎧豪,安 程
(中山大學地理科學與規劃學院,廣東 廣州 510275)
汛期是水庫發揮防洪功能的主要時期;因此,汛期的科學劃定是保障水庫更好發揮防洪調度功能的前提[1-2]。以往學者的研究多為直接對收集的水文資料進行整體分析,篩選出反映汛期特性的要素,如降水、蒸發、徑流等,然后直接對總體數據進行汛期研究劃分[3- 8]。這種方法雖然較為直接,但未能考慮水文要素的周期性變化規律特性。籠統地對全部要素進行分析可能會掩蓋部分因周期性變化產生的現象和趨勢,不能全面反映汛期各分期的變化過程。因此,本文在分析水文要素周期性變化規律的基礎上,以周期長度為時間尺度,借鑒滑窗理論[9],研究水庫汛期分期問題,旨在為現有的汛期分期研究方法提供新的思路和參考。
研究首先采用小波分析方法對水文要素進行周期性分析,從而確定水文要素變化的主周期;其次,以譜聚類為基礎,對長系列水文要素進行分析,確定非汛期與汛期的界限;之后,以主周期為時間尺度,參考滑窗理論,研究主周期長度時間尺度下,汛期各分期的變化規律;最后,通過各個分期的起止時間,確定最終的汛期各分期。
小波分析(Wavelet Analysis)是以Fourier變換為基礎,通過利用一簇小波函數系列來近似逼近某一信號序列。其常用于分析水文氣象要素的多時間演變特征[10-12]。在周期性研究中,Morlet小波函數應用較為廣泛,其小波函數和連續變換表達式為
(1)
(2)
式中,Wf(a,b)為小波變換函數;a、b分別為尺度因子和平移因子。小波方差值可以反映水文要素序列能量波動隨時間尺度的變化狀態。因此,小波方差圖中極大值點對應的時間尺度為水文序列變化的主周期。小波方差表達式為
(3)
譜聚類(Spectral Clustering)算法來源于譜圖劃分理論。該算法的主要思想是把所有聚類的對象看作是某個維度空間中的點,這些點用邊線進行連接。邊線上具有權重,距離較遠的兩個點之間的邊線權重較小,距離較近的兩點之間的邊線權重較大,通過對所有數據點進行切割,使得切割后不同的子圖間邊線權重盡可能小,而子圖內的邊線權重盡可能大,從而達到聚類的目的。目前,譜聚類算法在計算機視覺、機器學習、圖像識別等領域具有較為廣泛的用途[13]。
如果需要將圖G劃分成兩個子圖A和B,劃分的原則為:子圖內的邊線權重最大化,而各子圖間的邊線權重最小化,則可以建立目標函數Ncut。有
(4)
(5)
(6)
式中,cut(A,B)表示子圖A和B的邊線,也稱作“邊切集”;此外,若能實現目標函數的最小化,則分割后的結果自然能滿足不同類樣本間的相似度最小,同時也能滿足同一類樣本相識度最大。如果需要同時劃分成多個子圖,則上式可以調整為

(7)
集對分析(Set Pair Analysis)是一種將中國哲學中的對立統一和普遍聯系進行有機統一,從全新的角度分析變量之間不確定關系的方法。該方法的核心思想是對不確定系統的兩個有聯系的集合構建集對,通過對集對的特性進行統一性、差異性、對立性分析,建立集對的同、異、反聯系度,從而定量分析變量之間的不確定關系[14-16]。
若有集合X和Y,其均有n項屬性表征其特性,則有X(a1,a2,…,an)和Y(b1,b2,…,bn),則可以構建集對H(X,Y),那么表征集對H(X,Y)聯系度的表達式為
(8)
式中,S為同一性個數;F為差異性個數;P為對立性個數,三者總和為n;I為差異度系數,其取值在(-1,1)區間,根據集合反映的具體現象確定,但有時僅作為表示差異的符號;J為對立系數,取值為-1;若令a=S/n,b=F/n,c=P/n,則有μX~Y=a+bI+cJ;式中,a,b,c分別為集對H(X,Y)的同一度,差異度和對立度。
對于差異度,其可以進行進一步挖掘擴展,如令bI=b1I1+b2I2+…,則可得到更多元的聯系度
μX~Y=a+b1I1+b2I2+…+bk-2Ik-2+cJ
(9)
式中,a+b1+b2+…+bk-2+c=1;b1,b2,…,bk-2為差異度分量,表示分量的不同程度,如輕度、中度、重度等;I為不同等級的差異性系數。
研究以珠江水系紅水河上游龍灘水庫為對象,龍灘水庫位于廣西壯族自治區天峨縣境內,處中國西南部桂、黔、滇三省交界地帶,兼具發電、防洪、灌溉等多項功能。壩址以上流域面積為9.85萬km2,占據整個紅水河流域面積的71%,有效庫容為205.3億m3,具有多年調節能力。流域內多年平均降水量為760~1 860 mm,多集中在汛期。其中,汛期降水占據全年降水總量的85%以上。
研究的水文要素采用龍灘水庫下游天峨水文站1960年~2012年逐日徑流量數據。考慮到汛期洪水在時段上的延續性,以“候”為時間單位[17](古代氣候學中,一候是五天,三候為一節氣,一年二十四節氣,對應七十二候)。通過對53 a的逐日徑流量數據進行整理分析,得到影響汛期的5個因子,分別為多年平均候流量、多年平均候內最大1日洪量、多年平均候內最大3日洪量、多年平均候內最大3日洪量變差系數Cv、候內年最大洪峰出現次數,共同組成分期指標體系。
徑流周期分析是對龍灘水庫53 a的年徑流序列進行分析,尋找徑流周期性變化的規律,確定年徑流變化的主周期,以此作為滑動窗口的時間尺度。通過采用Morlet小波分析方法,對龍灘水庫年徑流序列變換分析,得到徑流Morlet小波實部圖和小波方差圖(見圖1、2)。

圖1 年徑流Morlet小波實部

圖2 年徑流Morlet小波方差
從圖1、2可知:①年徑流以22~26 a為尺度的周期震蕩信號最為明顯,20~25 a為尺度的小波實部圖中曲線高度閉合,整個值域內的特征明顯;②在小波方差圖中,共出現了4個極大值點,分別對應5、10、15 a和24 a。其中,24 a對應的小波方差峰值最大,表明震蕩能量最為劇烈,其余3個峰值相差不大。因此,年徑流的第一主周期為24 a,以此時間尺度作為汛期分期滑動窗口的“窗口寬度”。
水文要素的特性在汛期和非汛期時域內表現特征是不同的,汛期要素相較非汛期會發生明顯的變化。對于水庫流域而言,為了能夠從容地應對汛期洪水災害,實現適應性調度,一般采用年度長序列進行汛期和非汛期的劃分,汛期起止時間相對固定。汛期和非汛期的水文要素應當歸屬為兩個不同的類別。為此,利用譜聚類的獨特優勢,對指標體系進行聚類分析,進而確定汛期和非汛期的時間界限。
首先,對分期指標體系進行歸一化處理,消除指標之間因數據量級差異對分類結果產生的影響;其次,以一候為一個樣本,每個樣本由5個指標組成,對應5維空間上的一個向量。之后,分別計算5維空間中相鄰兩個向量之間的邊線權重,并構建目標函數。最后,通過對目標函數進行優化分析,使其達到最小值來確定樣本的歸屬類別。指標體系中共有72候,將此72候進行分割后,對于同一類的樣本用相同編號標識,汛期類別序號為2,非汛期為1,見圖3。

圖3 譜聚類劃分結果
從譜聚類劃分結果可知,72個候可以歸為兩類,標識非汛期的候序號為{1~21,63~72},其對應的時間段分別為1月1日~4月15日,11月11日~12月31日;相應地,標識汛期的候序號為{22~62},其對應的時間段為4月16日~11月10日。從汛期時間結果來看,相比傳統的簡單以月為單位的劃分情況(常以4月~10月為汛期),以候為單位則更加精確和科學。從汛期的起止時間和延續長度分析,汛期的起始時間延后了15 d,結束時間延后了10 d,汛期總體持續時間約為7個月,和傳統汛期持續時間長度基本保持一致。
通過譜聚類分析,確定了汛期的起止時間,即第22候~第26候。通過小波變換周期分析,確定年徑流序列的主周期為24 a。以兩者計算結果為基礎,引入滑窗理論,從1960年開始,以24 a年為尺度,分別將53 a的候序列進行滑動分組,如第1組為1960年~1983年,第2組為1961年~1984年,依次類推,共計30組,逐組進行汛期時間段(第22~62候)分期指標體系構建,然后以集對分析方法對汛期各個分期進行劃定。
針對每一組,計算過程大致可分為3個步驟:
(1)建立劃分標準集合Y={Y1,Y2},其中Yi的形狀與分期指標體系一致(汛期有41個候,共5個指標,為41×5的矩陣),且令Y1中所有元素均為1,則有Y1={1,1,…,1},令Y2中所有元素均為0,則有Y2={0,0,…,0};本研究中Y1標識主汛期,Y2標識非主汛期,依據時間確定為前汛期或者后汛期。
(2)建立分期指標體系符號集合M,M的形狀與分期指標體系一致,其中元素的確定方法為,以各項指標的均值為標準,對分期指標體系中的各個元素逐項進行判別歸類,當元素大于均值時,符號集合中的與該元素相同位置標記為1,反之標記為0。由此可確定所有分期指標體系中的符號元素類型,也就是指標體系符號集合。

表1 1960年~1983年汛期各個候聯系度及所屬類
(3)構建集對H=(Mi,Yj),并計算集對的聯系度。其中,Mi為分期指標體系符號集合,Yj為分期標準集合。通過將集合Mi中的元素逐一與Y1,Y2進行對照比較,統計符號相同的個數,記作S;統計符號相差為1的個數,記作F。根據下式可以得到2元聯系度。式中,因F為差異度系數,根據經驗取值為0,通過計算符號集合M的聯系度,可以判斷每候的歸類。即
(10)
式中,S1,S2分別為集合Mi中的元素逐一分別與Y1,Y2進行對照比較時符號相同的個數;F1,F2為相差為1的個數。
通過計算可得每組汛期的41個候的聯系度,以第一組1960年~1983年為例,其聯系度見表1。
從表1劃分結果可知,候序號{22~31}和{49~62}的所屬類均為Y2,候序號{32~48}所屬類為Y1,因此該組可以將汛期劃分為前汛期、主汛期和后汛期共3期。前汛期時域為4月16日~6月5日;主汛期時域為6月6日~8月31日;后汛期時域為9月1日~11月10日。
通過時間尺度進行滑動,可得到每組汛期的分期情況,為了更直觀地表現劃分結果,以候序列為橫坐標,以起始年份為縱坐標(如第1組年份為1960年,代表范圍1960年~1983年),對每組汛期分期情況進行統計,結果如圖4所示。圖中方框中○表示所屬類為Y2的候,方框中●表示所屬類為Y1的候。

圖4 各組汛期劃分結果匯總
從30個組的整體情況看,每組的前汛期結束時間點變化不大,基本在第30、31和32候結束。對于主汛期,情況則較為復雜,存在類似兩個Y1中間夾雜一個Y2的情況,如{○,●,○}。為避免單個不同類別摻雜影響劃分的情況,規定只有當Y2前后面均連續出現2次及以上的Y1時,如{●,●,○,●,●},此Y2方可轉化為Y1,其余情況均以Y2首次出現的時間為主汛期截止時間。因此,按照此法重新進行修正后得到圖5。

圖5 各組汛期劃分結果修正后匯總
從圖5可知,隨著年份的向后滑動,主汛期的時間范圍有變大的趨勢,從第1組的17個候(約102 d),變化到第21組的23個候(約132 d),整體而言,主汛期結束的時間點多為第48~58候。對前汛期和主汛期結束時間進行統計分析,可得到表2。

表2 前汛期和主汛期結束候序號統計比例
由表2可知,前汛期在第30候結束時的比例最大,為73.33%,平均候序號為30.37(約為第30候),其對應時時間為4月30日;主汛期在第53候結束的比例最大,為43.33%,平均候序號為51.06(約為第51候),其對應的時間為9月15日。因此可以確定龍灘水庫的汛期各分期起止時間,前汛期為4月16日~5月31日、主汛期為6月1日~9月15日、后汛期為9月16日~11月10日。
結合之前有關學者對龍灘水庫汛期分期的相關研究,將結果進行匯總比較(見表3)。

表3 不同劃分方法研究結果對比
由表3可知,史文海[18]采用分形方法對龍灘水庫進行劃分時,劃定汛期起訖時間的依據是樣本序列中第一場暴雨和最后一場暴雨出現的時間,如此進行劃定略顯粗糙。且該法在進行分期時段劃分時是以旬(10 d)為單位,不夠精確。模糊集合分析法[19]在進行汛期和非汛期的劃分時則是直接采用珠江全流域的劃分標準,以4月開始10月結束,未能考慮和捕捉到各子流域在汛期時域上存在的微小差異性。考慮水文周期的劃分方法兼顧考慮了劃分時間的精準性和子流域間汛期時域存在的差異性情況,先對汛期時域進行劃定,在此基礎上考慮水文要素的周期性變化規律,進行汛期各分期的劃分,從而使得劃分結果更科學可信。
(1)本研究首次將水文要素的周期性變化規律引入到水庫汛期分期的研究中,以水文要素的周期為時間尺度,借鑒滑窗理論,對時間尺度內的汛期進行劃分。對于汛期各子分期劃分時,采用能同時反映要素的宏觀性和微觀聯系的集對分析方法。最后,對劃分結果進行統計分析,并與之前有關研究進行比較,從而確定龍灘水庫的汛期分期情況。
(2)科學的汛期分期可以為更好地計算各分期的設計洪水,提出適應性調度規則,制定科學的水庫管理方案提供有力支撐。