黃河流域重要斷面設計年徑流量計算研究

杜 懿，孟 越，陳昱楨，王大剛

（中山大學 地理科學與規(guī)劃學院，廣東 廣州 510275）

水文頻率分析的目的是為水利工程的規(guī)劃、設計和建設等提供科學依據(jù)，其首要工作就是對研究序列進行“三性”（可靠性、代表性和一致性）審查。 在過去相當長一段時間內，人類活動的范圍和強度有限，對自然環(huán)境的影響不大，一般經(jīng)整編后的水文資料序列均能滿足一致性要求。 但是，近年來，隨著全球氣候變化以及人類活動對自然界的大規(guī)模改造，水文時間序列的一致性要求不再得到普遍的滿足［1］。 為此，水文工作者進行了大量的相關研究，涌現(xiàn)出眾多針對非一致性水文頻率分析的新方法和新理論，可以簡單概括為間接法和直接法兩大類。

間接法主要基于“還原或還現(xiàn)”思想，即將原始的非一致性水文序列經(jīng)過一系列變換，使得序列重新滿足一致性要求后再進行頻率分析，有基于趨勢分析和基于跳躍分析等非一致性分析方法［2-4］；直接法則主要有分解合成法［5］、混合分布法［6］、時變矩法［7］、條件概率分布法［8］、廣義加法模型［9］等。 以上方法在我國各地區(qū)、各流域有著廣泛的應用，如成靜清等［10］將基于混合分布的非一致性頻率分析方法應用于陜北及關中地區(qū)的實測年徑流量序列處理中，研究表明理論頻率和經(jīng)驗分布擬合較好；胡義明等［11］提出了基于跳躍分析的非一致性水文頻率計算方法，并對金沙江流域某水文站資料進行了研究；謝平等［12］提出了一種基于小波分析的非一致性洪水頻率計算方法，并在西江梧州站做了實例驗證；劉丙軍等［13］運用時變矩方法研究了西北江三角洲地區(qū)的非一致性洪水頻率分析問題；李伶杰等［14］對條件概率模型在非一致性水文序列頻率計算中的應用進行了研究，該研究為變化環(huán)境下的水文頻率分析提供了新方法；張冬冬等［15］以GAMLSS模型為基礎，分析了大渡河流域降水頻率的非一致性特征，結果表明GAMLSS 模型可以很好地反映水文序列的非一致性特征。

1 研究對象和方法

以黃河流域上、中、下游的三大控制站（蘭州站、花園口站、利津站）的年徑流量時間序列為研究對象，在序列特性分析的基礎上，采用不同方法進行各站的設計年徑流量計算。 三大站從上到下分別為蘭州水文站、花園口水文站和利津水文站，其中蘭州水文站和花園口水文站具有1919—2018年共100 a 長度的年徑流量資料，利津水文站有1952—2018年共67 a 長度的年徑流量資料。

利用謝平等［16］提出的水文時間序列變異診斷系統(tǒng)對3 個站的年徑流量時間序列進行變異診斷（見圖1）。 如果研究序列存在變異特征，說明序列不滿足一致性要求，則需采用非一致性頻率計算方法來進行設計年徑流量計算；如果序列不存在顯著變異，說明序列保持著較好的一致性，則采用傳統(tǒng)水文頻率分析方法進行設計年徑流量計算。

圖1年徑流量時間序列的變異診斷

2 研究過程

2.1 水文時間序列一致性診斷

（1）初步診斷。 綜合利用過程線法和累計值法分別對蘭州站、花園口站和利津站的年徑流量序列進行一致性分析，這兩種方法均屬于直觀的定性分析，因原理簡單而被廣泛應用于水文時間序列一致性的初步診斷中。 具體診斷結果如圖2 所示。

圖2 三大控制站年徑流量序列一致性的初步診斷

由圖2可見，1919—2018年蘭州站和花園口站的年徑流量序列變化趨勢不明顯，而利津站1952—2018年的年徑流量序列呈現(xiàn)出較為明顯的下降趨勢；三大站的時序累計值結果顯示，蘭州站的年徑流量序列與45°斜率直線的擬合效果最好，而花園口站和利津站的擬合誤差則相對較大。 綜合分析，可以初步確定花園口站和利津站年徑流量序列的一致性可能遭到了破壞。

為了進一步準確地描述3 個水文站年徑流量序列的一致性變化情況，需要進行詳細診斷。

（2）詳細診斷。 一般地，可以認為造成水文時間序列發(fā)生變異的主要原因是序列出現(xiàn)了趨勢、跳躍和周期等成分。 由于目前我國絕大多數(shù)水文站的時間序列長度較短，較早建設的站也只有幾十年的觀測資料，因此周期成分一般不顯著，可不予考慮。 下面僅對趨勢性成分和跳躍性成分進行識別和檢驗。

分別采用相關系數(shù)法、線性趨勢回歸檢驗法、Spearman 秩次相關檢驗法以及Kendall 秩次相關檢驗法等［17］對蘭州、花園口和利津3 個站的年徑流量序列的趨勢性成分進行識別和檢驗，顯著性水平均取0.05，檢驗結果見表1。

表1 各站年徑流序列的趨勢性成分識別與檢驗結果

由表1 中的檢驗結果可以看出，蘭州站和花園口站近100 a 來的年徑流量序列變化的趨勢性不顯著；而利津站1952—2018年的年徑流量序列則呈現(xiàn)出較為顯著的趨勢性成分，且表現(xiàn)為下降趨勢。 以上結論與初步診斷結論一致，再次證明利津站的年徑流量序列已不再具備一致性，存在趨勢性變異，傳統(tǒng)的一致性水文頻率分析方法失效。

下面對蘭州站、花園口站和利津站的歷史年徑流量序列進行跳躍性成分的識別與檢驗。 雖然常用的水文時間序列突變檢驗方法眾多，但各方法的側重點不同，最終的檢測結果常有差異［18］。 本文主要采用檢驗結果較為一致的有序聚類分析、Lee-Heghinian 法以及Mann-Kendall 非參數(shù)檢驗法（M-K 法）對蘭州站、花園口站和利津站的年徑流量序列進行跳躍性變異診斷，結果如圖3～圖8 所示。

圖3 基于有序聚類分析和Lee-Heghinian 法的蘭州站年徑流量序列跳躍性診斷結果

當顯著性水平α取0.05時，蘭州站年徑流量序列有序聚類分析和Lee-Heghinian 法的統(tǒng)計量｜ T ｜=3.85＞tα／2=1.64，說明蘭州站1919—2018年年徑流量時間序列存在較為顯著的突變點，最大可能變異位置為1932年；觀察圖4可以看出，M-K 法檢測到了多個突變點，其中就包括了1932年。 綜合以上兩種方法的突變檢驗結果，確定1932年為蘭州站年徑流量序列的最大可能突變點。

圖4 基于M-K 法的蘭州站年徑流量序列跳躍性診斷結果

當顯著性水平取0.05 時，花園口站年徑流量序列有序聚類分析和Lee-Heghinian 法的統(tǒng)計量為4.21，說明花園口站1919—2018年年徑流量時間序列存在較為顯著的突變點，最大可能變異位置為1990年；觀察圖6可以看出，M-K 法檢測到的突變點也為1990年。綜合以上兩種方法的突變檢驗結果，確定1990年為花園口站年徑流量序列的最大可能突變點。

圖6 基于M-K 法的花園口站年徑流量序列跳躍性診斷結果

當顯著性水平α取0.05時，利津站年徑流量序列有序聚類分析和Lee-Heghinian 法的統(tǒng)計量｜ T ｜=5.04＞tα／2=1.64，說明利津站1952—2018年年徑流量時間序列存在較為顯著的突變點，最大可能變異位置為1985年和1990年；觀察圖8，可以看出M-K 法檢測到的突變點為1985年。 綜合以上兩種方法的突變檢驗結果，確定1985年為利津站年徑流量序列的最大可能突變點。

圖7 基于有序聚類分析和Lee-Heghinian 法的利津站年徑流量序列跳躍性診斷結果

圖8 基于M-K 法的利津站年徑流量序列跳躍性診斷結果

綜上，可以確定3 個水文站的歷史年徑流量序列均發(fā)生了變異，各序列不再滿足一致性要求。 其中：蘭州站和花園口站的年徑流量時間序列主要發(fā)生了跳躍性變異，而利津站的年徑流量時間序列既發(fā)生了趨勢性變異又發(fā)生了跳躍性變異。 針對各序列的不同變異特性，分別采用相應的非一致性水文頻率分析方法進行設計年徑流量計算。

2.2 水文頻率分析計算

2.2.1 蘭州站年徑流量序列頻率分析

針對蘭州站年徑流量序列的跳躍性變異特征，采用基于跳躍診斷的二次修正法進行非一致性水文頻率計算。 對于年徑流量序列xi，假設變異點為τ，首先根據(jù)變異點位置將原始序列劃分成前后兩個子序列，即（x1，x2，…，xτ） 和（xτ＋1，xτ＋2，…，xn） 。可以得到將變異點之前的序列還原到變異點之后狀態(tài)的表達式為

式中：Xt、St分別為修正前、后的變異點之前的子序列；Ex1、Ex2分別為原始序列變異點前、后兩個子序列的均值。

對還現(xiàn)處理后的新序列進行一致性檢驗，若仍不滿足一致性要求，則進行下一步處理。

如圖9，修正后的序列變異強度雖然有所弱化，但仍然存在一定程度的跳躍特征，以1989年最為顯著，故需要進一步處理以完全消除序列的跳躍特征。

圖9 一次修正后的基于有序聚類分析和Lee-Heghinian 法的蘭州站年徑流量序列跳躍性診斷結果

對一次修正后的序列繼續(xù)進行變異診斷，根據(jù)新的變異點位置將序列重新劃分成兩個子系列，均值仍然記為Ex1和Ex2，設Ex為平穩(wěn)狀態(tài)的振動中心，則有：

根據(jù)文獻［3］，權重γ的計算公式以及二次修正后的序列可表示如下：

從圖10可以看出，經(jīng)過連續(xù)兩次修正，蘭州站的年徑流量序列已經(jīng)不再表現(xiàn)出突變特征（UF、UB曲線無交點），說明二次修正后的年徑流量序列滿足一致性要求。

圖10 二次修正后的基于M-K 法的蘭州站年徑流量序列跳躍性診斷結果

對二次修正后的年徑流量時間序列進行水文頻率分析，以P-Ⅲ型曲線為分布線型，采用優(yōu)化適線法計算均值、變差系數(shù)和偏態(tài)系數(shù)等參數(shù)，蘭州站設計年徑流量計算結果見表2。

表2 蘭州站設計年徑流量計算結果

2.2.2 花園口站年徑流量序列頻率分析

由圖5、圖6 的分析結果可知，花園口站年徑流量序列的最大可能變異位置為1990年，則原序列被劃分成1919—1990年和1991—2018年兩個子序列。 兩段序列的長度均能滿足水文頻率分析的要求，故宜采用混合分布法計算花園口站的設計年徑流量。

圖5 基于有序聚類分析和Lee-Heghinian 法的花園口站年徑流量序列跳躍性診斷結果

混合分布法的主要思想是將非同分布的樣本序列視作若干個子分布混合而成，實際應用中為簡化計算，多將原始序列分解成兩個子序列，并將各子序列占原始序列長度的比例作為權重系數(shù)進行連接。可用以下公式進行描述：

式中：F（x） 為混合分布函數(shù)；F1（x） 和F2（x） 分別表示變異點前、后序列的分布函數(shù)；β為變異點前的子序列占原始序列長度的比例。

以P-Ⅲ型曲線為分布線型，采用優(yōu)化適線法計算P-Ⅲ型分布中的參數(shù)，分別對1919—1990年、1991—2018年兩個年徑流量子序列進行水文頻率分析，再將各子序列的設計值進行疊加，即可得到原序列的設計值。 花園口站設計年徑流量計算結果見表3。

表3 花園口站設計年徑流量計算結果

2.2.3 利津站年徑流量序列頻率分析

由前述分析可知，利津站歷史年徑流量序列既存在著趨勢性變異又存在著跳躍性變異。 為了更好地刻畫出該序列的趨勢性特征，本文采用分解合成法計算利津站的設計年徑流量。

先對原始年徑流量序列進行db3 小波分解，分解深度取6 級（d1～d6）。 如圖11 所示，s為原始序列，原始序列被分解成了7 個部分，其中a6為序列的趨勢項（單位為億m3）。可以看出，趨勢項序列呈現(xiàn)出較為明顯的下降趨勢，說明利津站1952—2018年年徑流量有逐年減少的趨勢，這一現(xiàn)象與前述分析結論也是一致的。

圖11 利津站年徑流量序列的db3 小波分解結果

一般把趨勢項看作是原始序列的確定性成分，而把趨勢項以外的其余部分視作隨機性成分，并假定隨機性成分仍保持著一致性。 對確定性成分進行函數(shù)趨勢擬合，根據(jù)得到的擬合方程計算未來年份的設計值（見圖12）；對隨機性成分的設計值采用傳統(tǒng)的一致性水文頻率分析方法進行計算（見表4）。

圖12 確定性成分的線性擬合結果

表4 利津站年徑流量隨機性成分序列設計值

將確定性成分的分析結果和隨機性成分的計算結果進行疊加合成，即可得到利津站各頻率的設計年徑流量。 為了體現(xiàn)出年徑流量序列隱含的趨勢特征，表5 給出了過去情景下（2010年）、現(xiàn)狀情景下（2020年）和未來情景下（2030年）的年徑流量設計值。

表5 利津站在不同情景水平年下的年徑流量設計值

3 結論

（1）經(jīng)初步診斷和詳細診斷，3 個水文站的年徑流量序列均不滿足一致性要求，無法直接進行水文頻率分析。 其中：蘭州站和花園口站均存在顯著的跳躍性變異，變異位置分別為1932年和1990年；利津站既存在趨勢性變異又存在跳躍性變異，整體上年徑流量序列呈現(xiàn)下降趨勢，最大可能突變位置為1985年。

（2）根據(jù)3 個水文站年徑流量序列的變異特點，分別采用不同的非一致性頻率計算方法，包括基于跳躍診斷的二次修正法、混合分布法和分解合成法，計算得到各水文站不同重現(xiàn)期的設計年徑流量。

（3）經(jīng)計算，蘭州站和花園口站的百年一遇設計年徑流量分別為498.51 億m3和852.14 億m3，利津站在3種不同情景下的百年一遇設計年徑流量分別為921.34 億、917.19 億、913.04 億m3。