基于水文分区的海河流域洪水演变特性分析

陈 旭，韩瑞光

（海河水利委员会水文局，天津 300170）

海河流域位于东经112°～120°、北纬35°～43°，东临渤海，南界黄河，西靠云中、太岳山，北倚蒙古高原，流域总面积31.8 万km2，呈扇形分布，具有河系复杂、水系众多、过渡带短、源短流急等特点。流域洪水年际变化很大，一般洪水与特大洪水在量级上相差悬殊。流域内河系众多，流域面积小，分流入海，降雨落区预报难度大。加之地下水超采严重，下垫面变化剧烈，加大了洪水预报的难度。国内学者针对海河流域洪水特性分析已开展了大量的研究工作[1-3]，这些研究多集中在某一子流域或某一场暴雨洪水关系，且研究方法单一，研究成果对海河流域洪水演变特性分析缺乏系统性和全面性。因此，本文以海河流域为研究对象，采用多种统计方法从趋势、突变和周期等方面开展流域洪水演变特性分析，从中长期的角度研究洪水的演化规律，为流域洪水定性评估提供参考。

1 研究方法

1.1 主成分分析法和聚类分析法

主成分分析法（Principal Component Analysis，简称PCA）能够保留研究要素的关键信息，同时达到简化数据和降维的目的[4]。

K均值（K-means）聚类法通过方差分析确定最优分类数，根据各个样本到聚类中心的距离，将样本进行样本划分，然后在此计算聚类中心，通过对比前后两次聚类中心的变化大小，确定算法是否达到最优解[5]。

1.2 趋势分析法

为增加趋势分析结果的可靠性，本文同时采用滑动平均法[6]以及非参数检验方法Mann-Kendall（以下简称M-K）法和Spearman 秩次相关检验法[7-9]对海河流域洪水要素的变化趋势进行分析。

研究表明，序列间的相关性将增加显著性趋势检验概率[10]，因此在进行趋势检验之前，需要对数据作消除相关性处理。本文采用趋势自由预白化（TFPW）方法来消除序列间的相关性，在此基础上采用非参数研究法对海河流域洪水要素的趋势特性进行分析。

1.3 突变分析法

受人类活动影响，水文要素自然变化情况常常表现出突发性，因此对水文时间序列进行突变分析有助于掌握其变化趋势。常用的突变分析方法主要有M-K法[11]、佩蒂特法[12，13]、勒帕热法[14]、滑动t检验法[15]等，各方法在原理上具有一定的相似性。本文主要采用滑动t检验法、佩蒂特法、勒帕热法和BG（Bernaola-Galvan）启发式分割法[16]对流域洪水序列进行突变检验。

1.4 Morlet小波法

小波变换是基于傅里叶（Fourier）变换的一种时间-频率（时间-尺度）分析方法，其原理与傅里叶分析类似，但小波分析用的是多种能衰减的小波基通过伸缩和平移来表示某个信号。由于小波分析比傅里叶分析多了一个时域，所以其在处理非平稳水文时间信号中有出色表现[17，18]。因此，本文采用小波理论对洪水序列进行周期分析。

2 海河流域水文分区

2.1 水文分区影响因素

下垫面要素和气候要素是影响流域产汇流变化的两大主导因素（见图1）。下垫面因素主要包括：①地形特征，主要包括平均坡度、最大高程差和坡度＜1%面积比例；②土地利用；③土壤；④植被覆盖。气候因素主要包括降水和蒸发。降水是最主要、最直接的洪水来源，而蒸发则是一种水分损失。

影响水文分区的下垫面和气候要素提取结果，如图1所示。

图1 海河流域影响水文分区的下垫面及气候要素空间分布

2.2 水文分区划分

在DEM 的基础上通过流向计算、洼地提取、阈值计算、填洼及流域计算等过程将海河流域划分为425 个子流域，在此基础上，对影响洪水变化的下垫面和气候要素按子流域进行分析提取，得到海河流域不同子流域的下垫面和气候因子。

为了消除原始因子的重叠信息，需要将提取的各子流域下垫面和气候因子进行主成分分析。首先计算原始变量累计方差贡献并选取贡献率大于0.92时的特征根个数作为主成分个数，本文最终提取5个主成分。然后利用四次方最大法对因子荷载矩阵进行正交变换来简化因子载荷阵的结构，使因子便于解释和命名。表1给出了下垫面和气候因子的主成分荷载矩阵，可以看出第一主成分刻画了坡度、最大高程差、耕地面积比例及林地面积比例5个变量，第二主成分刻画了坡度＜1%面积比例和土壤砂含量2个变量，第三主成分刻画了植被覆盖度和降水2个变量，第四主成分刻画了草地面积比例，第五主成分刻画了蒸发。

表1 下垫面和气候因子主成分荷载矩阵

将主成分分析计算得到的5个主成分作为K均值聚类分析的变量。综合考虑研究的实际需要及最小类内方差原则，海河流域425个子流域被划分为5类，即5个水文类型分区（见图2），分别为位于东南沿海的第一水文分区、北部和西部的第二水文分区、东南部和西南部的第三水文分区、西北部和海河中部、南部平原区的第四水文分区、中东部的第五水文分区。

图2 海河流域水文类型分区

3 海河流域洪水演变特性分析

根据水文分区结果分别选取分区内典型水文站点，对洪水要素进行特性分析，其中第二、三、四和五分区选定的站点分别为王快水库、大黑汀水库、册田水库和张坊。第一水文分区位于东南沿海，覆盖面积较小，区域内没有具有长系列洪水资料的典型水文站点，故本文未对第一水文分区进行研究。

3.1 趋势特性分析

本文同时采用滑动平均法、M-K 法和Spearman秩次相关检验法对海河流域不同水文分区所选典型水文站的洪峰序列和最大3 d 洪量序列进行趋势分析，得到各水文分区洪峰、3 d洪量序列的变化趋势，如图3所示，详见表2—4。

表2 各水文分区洪峰、3 d洪量趋势变化分析结果

表3 各水文分区洪峰、3 d洪量自相关系数

表4 各水文分区洪峰、3 d洪量变化趋势检验

根据海河流域各水文分区洪峰、3 d 洪量的5 a滑动平均过程线（见图3）和线性倾向估计（见表2）可知：各水文站洪峰、3 d 洪量的变化趋势表现出较强的一致性，趋势变化基本相同；各水文站洪峰、3 d洪量的倾向值均小于零，即海河流域洪峰、3 d 洪量整体表现为下降趋势，其中第三水文分区洪峰、3 d洪量变化率最大，分别为-61.072 m3/s/a、-0.068亿m3/a，第四水文分区最小，分别为-20.957m3/s/a、-0.012亿m3/a；第二水文分区洪峰、3 d 洪量在1959—1961、1967—1972、1979—1986 和 1992—2012 年时间段表现为下降趋势，在1955—1959、1961—1963、1972—1979和1986—1992 年时间段表现为上升趋势；第三水文分区洪峰、3 d 洪量在1962—1971 年时间段波动下降，1955—1962年时间段波动上升，1979—2011年时间段表现为下降趋势，1971—1979 年时间段表现为上升趋势；第四水文分区洪峰、3 d 洪量在1967—1972 年时间段波动下降，1958—1967 和1972—1981 年时间段波动上升，1955—1958、1981—1989和1995—2011 年时间段表现为下降趋势，1989—1995 年时间段表现为上升趋势；第五水文分区洪峰、3 d 洪量在1956—1961 年时间段波动下降，1963—1966 年时间段表现为稳定波动状态，1966—1971 年时间段表现为下降趋势，1961—1963 年时间段表现为上升趋势，1971 年之后洪峰表现为稳定波动状态，洪量在1971—1977 年时间段表现为上升趋势、1977—1986 和 2000—2011 年时间段表现为下降趋势、1988—2000 年时间段波动上升。

图3 各水文站洪峰、3 d洪量径流变化及其5 a滑动平均过程

为了定量分析该研究区域洪峰、3 d洪量的变化规律，本文选用M-K 法和Spearman秩次相关检验法进一步对研究区域的洪水要素趋势变化特性进行分析。运用M-K 法对序列进行趋势检验前，首先需要计算洪水要素序列的自相关系数，计算结果详见表3。由表3 可以看出，除第三水文分区洪水要素外，其余水文分区洪水要素序列的自相关系数均大于0.1，因此在进行M-K法检验前需要对其进行预白热化处理。利用M-K 法和Spearman 秩次相关检验法对各水文分区洪峰、3 d 洪量序列进行趋势性检验，检验结果详见表4。由表4 可以看出，2 种统计方法一致检测到4 个水文分区洪峰、3 d 洪量序列均呈现显著的下降趋势，其中第四水文分区下降最为显著、第五水文分区次之。分析各站的Kendall 倾斜度β可知，4 个水文分区洪峰、3 d 洪量序列的 Kendall 倾斜度均为负值，这与M-K 法和Spearman秩次相关检验法统计检验的结果吻合。第三水文分区洪峰、3 d洪量减少速率最快，分别以28.293 m3/s/a、0.032 亿m3/a的速率递减，最慢的为第五水文分区，这与滑动平均法略有差别，主要是因为M-K 法利用中值来计算变化率，而线性倾向估计法直接计算序列的斜率得到倾向值。

3.2 突变特性分析

本文同时运用滑动t 检验法、佩蒂特法、勒帕热法和BG 启发式分割法对海河流域不同水文分区所选典型水文站的洪峰和最大3 d 洪量序列的突变特性进行分析，由于篇幅限制，本文仅给出第二水文分区王快水库的计算结果，如图4所示。

图4 第二水文分区洪峰、3 d洪量序列突变检验结果

由图4 可知，第二水文分区洪峰可能的突变点为 1964、1979、1985、1990、1992、1994、1996 和 1997年，其中1964、1996 和1997 年（延迟变异点）的突变点由特大值或者较大值造成，佩蒂特法和BG 启发式分割法同时诊断1990 年为突变点，佩蒂特法和勒帕热法同时诊断1992 年为突变点，因此综合考虑最终确定第二水文分区洪峰最可能的突变点为1990和1992 年；第二水文分区最大3 d 洪量可能的突变点为1964、1972、1974、1979、1982、1989 和2000 年，其中1964 年的突变点由特大值或者较大值造成，另外，王快水库洪峰序列的研究年限为1951—2012年，从统计上讲靠近序列末端2000 年作为突变点不可靠，滑动t 检验法和佩蒂特法同时诊断1979 年为突变点，因此综合考虑最终确定第二水文分区最大3 d洪量最可能的突变点为1979年。

第三—五水文分区各方法突变点诊断结果，详见表5。最终，确定第三水文分区大黑汀水库洪峰最可能的突变点为1979 年，最大3 d 洪量最可能的突变点为1979、1987 和1990 年；第四水文分区册田水库洪峰最可能的突变点为1967、1983 和1999 年，最大3 d 洪量最可能的突变点为1982 和1983 年；第五水文分区张坊站洪峰最可能的突变点为1964、1979 和1982 年，最大3 d 洪量最可能的突变点为1964、1979和1982年。

表5 海河流域4个水文分区典型站点洪峰、3 d洪量序列突变检验结果

3.3 周期特性分析

采用Morlet小波法对海河流域洪水要素的周期特性进行分析，计算得到小波系数实部等值线、小波系数模等值线、小波系数模方等值线以及小波方差，由于篇幅限制，本文仅给出第二水文分区王快水库洪峰、洪量Morlet小波法计算结果，如图5—6所示。

图5 第二水文分区洪峰小波分析

图6 第二水文分区3 d洪量小波分析

由第二水文分区小波系数实部等值线图5（a）可以看出，洪峰演变过程中存在着5～9、9～18和1～4 a的3 类尺度的周期变化规律，其中在5～9 a尺度上出现丰枯交替的14次震荡，且时段初至90年代中期表现较为稳定；在9～18 a 尺度上存在准7次震荡，且直到2012 年等值线未闭合，说明当前正处于洪峰偏枯期，未来一段时间洪峰将继续偏小；且这2个时间尺度在整个分析时段表现较为稳定，具有全域性；在1～4 a 尺度上出现丰枯交替的23 次震荡。图5（b）—（c）反映了第二水文分区洪峰在不同时间尺度上的周期震荡能量的强弱分布情况，在5～9 a尺度上周期震荡能量最强，但具有明显的局部特性（50 年代中期—60 年代末）；在9～18 a 尺度上周期震荡能量在1965年后减弱，且在1959年左右存在一个明显的震荡中心；在1～4 a 尺度上周期震荡能量最弱，在1963年左右存在一个明显的震荡中心。图5（d）的小波方差图中有3 个较为明显的峰值，依次对应着6、12和3 a 的尺度，分别为第二水文分区洪峰的第一、第二和第三主周期，这3 个周期波动控制着第二水文分区洪峰在整个时间域内的变化特征。

由图6（a）—（d）可知，第二水文分区王快水库洪量存在着 4～8、10～19 和 26～32 a 的 3 类尺度的周期变化规律。其中，洪量的第一、二主周期与洪峰的第一、二主周期表现出相似的演变特性。第三主周期在26～32 a 尺度上表现出丰枯交替的3.5 次震荡，且直到2012 年等值线未闭合，说明当前正处于洪量偏丰期，未来一段时间洪量将继续偏大。

第三—五水文分区Morlet 小波法周期分析结果，详见表6。

表6 海河流域4个水文分区洪水要素小波分析结果

第三水文分区大黑汀水库存在着17～32、4～10、10～16 和 1～4 a 的 4 类尺度的周期变化规律。第一、二、三和四主周期分别为26、6、12 和3 a，其中在17～32 和10～16 a 尺度上分别出现丰枯交替的3.5 次和7.5 次震荡，且在整个分析时段表现较为稳定，具有全域性；在 4～10 和 1～4 a 尺度上分别存在准 15.5 次和准24.5次震荡，且均在50年代中期—70年代末表现较为稳定；在4～10 a尺度上周期震荡能量最强，但具有明显的局部特性（50 年代中期—60 年代末），在1961 年左右存在一个明显的震荡中心；在17～32 a尺度上周期震荡能量在70 年代末后有所减弱；在10～16 a尺度上周期震荡能量较弱，无明显的震荡中心；在1～4 a 尺度上周期震荡能量在60 年代初到中期最强，在1963 年左右存在一个明显的震荡中心。第三水文分区大黑汀水库洪量存在着22～32、9～18和5～9 a的3 类尺度的周期变化规律，其中第一主周期（26 a）与洪峰的第一主周期表现出相似的演变特性，周期震荡能量最强，但具有明显的局部特性（70 年代中期—20世纪末）；第二主周期在9～18 a尺度上表现丰枯交替的6.5 次震荡，且直到2012 年等值线未闭合，说明当前正处于洪量偏丰期，未来一段时间洪量将继续偏大，其周期震荡能量在70年代中期—80年代中期较强，且在1979年左右存在一个明显的震荡中心；第三主周期在5～9 a尺度上表现丰枯交替的12.5次震荡，且70年代中期—2012年表现较为稳定，其周期震荡能量在80 年代最强，之后稳定震荡。

第四水文分区册田水库洪峰存在着19～29 和4～11 a 的2 类尺度的周期变化规律。第一、二主周期分别为23、8 a，其中在19～29 a尺度上出现丰枯交替的4.5次次震荡，其周期震荡能量具有明显的局部性，在70 年代末后明显减弱；在4～11 a 尺度上存在准12 次震荡，且在50 年代中期—70 年代末表现较为稳定，其周期震荡能量较强，但具有明显的局部特性（50年代中期—60年代末），在1953年左右存在一个明显的震荡中心。第四水文分区册田水库洪量与洪峰的演变特性表现出高度的一致性。

第五水文分区张坊洪峰存在着9～18、4～8和23～30 a 的3 类尺度的周期变化规律。第一、二、三主周期分别为12、6、26 a，其中在9～18 a尺度上出现丰枯交替的7.5 次次震荡，80 年代以前表现较为稳定，其周期震荡能量最强，但具有明显的局部特性（50年代中期—60 年代末），在1962 年左右存在一个明显的震荡中心；在4～8 a尺度上存在准15次震荡，且在50年代中期—70年代末表现较为稳定，其周期震荡能量具有明显的局部性，在70 年代末后明显减弱；在 23～30 a 尺度上存在准 3.5 次震荡，且在 70 年代中期—20 世纪末表现较为稳定，其周期震荡能量较弱。第五水文分区张坊洪量存在着9～18、4～8 和25～32 a 的3 类尺度的周期变化规律，其中第一、二主周期与洪峰的第一、二主周期表现出相似的演变特性；第三主周期在25～32 a尺度上表现丰枯交替的3 次震荡，且在整个分析时段表现较为稳定，具有全域性，直到2012 年等值线未闭合，说明当前正处于洪量偏丰期，未来一段时间洪量将继续偏大。

根据上述分析，海河流域4 个不同水文分区洪水均表现出 6～8 和 12～14 a 的周期性。截至 2012年，在长周期上洪水等值线均未封闭，正处于偏丰阶段，表明在长时间尺度上海河流域洪水将以较低能量继续偏丰；在短时间尺度上洪水等值线也未封闭，正处于偏枯阶段，表明海河流域洪水将以较高能量继续偏枯，按照最强震荡尺度6 和12 a 的周期特征推算，洪水偏枯状态将持续到2016—2022年。

4 结论

为系统掌握海河流域洪水演变特性，本文在水文分区划分的基础上，采用多种统计学方法对流域洪水要素的趋势、突变和周期特性进行研究，主要结论如下。

（1）基于下垫面和气候要素，海河流域共被划分为5 个水文分区，分别为位于东南沿海的第一水文分区、北部和西部的第二水文分区、东南部和西南部的第三水文分区、西北部和海河中部、南部平原区的第四水文分区、中东部的第五水文分区。

（2）4个水文分区洪峰、3 d洪量序列均呈现显著的下降趋势，其中第四水文分区下降最为显著、第三水文分区下降速率最快。

（3）第二水文分区洪峰的变异点为1990和1992年，最大3 d洪量的变异点为1979年；第三水文分区洪峰的变异点为1979 年，最大3 d 洪量的变异点为1979、1987 和1990 年；第四水文分区洪峰的变异点为 1967、1983 和 1999 年，最大 3 d 洪量的变异点为1982 和1983 年；第五水文分区洪峰的变异点为1964、1979 和 1982 年，最大 3 d 洪量的变异点为1964、1979和1982年。

（4）4个水文分区洪水要素均具有6～8和12～14 a的周期性，在短时间尺度上洪水等值线未封闭，正处于偏枯阶段，按照最强震荡尺度6 和12 a 的周期特征推算，洪水偏枯状态将持续到2016—2022年。