三峡坝前水位过程线特征参数的标定与统计特性

何金文,孙钧键,程 军

(三峡大学 水利与环境学院,湖北 宜昌 443002)

1 研究背景

三峡大坝是迄今为止世界上规模最大的水力发电工程，主要功能是防洪、发电与航运。水库汛末将库水位蓄至175 m，主汛期控制在145 m以腾出库容调蓄上游来水以护航下游行洪安全，洪峰消减可达40%[1]。库区每年30 m的水位消落与回升必然对岸坡的变形和稳定带来巨大影响。

三峡库区有数千处崩滑灾害，水位周期性变化，将导致老滑坡复活并诱发新滑坡[2]。Luo等[3]研究了三峡库区水位变化对藕塘古滑坡的变形特性及复活机制；尚敏等[4]认为三峡库水快速消落是影响白家包滑坡变形的主要因素；张富灵等[5]、卢书强等[6]分别研究了库区谭家湾滑坡、树坪滑坡变形随水位变化的关系；胡亚波等[7]认为库水消落将形成指向坡外的渗透压力；熊珅等[8]认为三峡库水消落速率的增大会导致八字门滑坡稳定性的减小。可见，库区岸坡的变形与稳定主要与库水变化与消落速率有关。目前三峡库区岸坡变形监测与稳定性计算绝大多数是根据运行调度图或坝前水位过程线进行的，而库区涉水岸坡的稳定性分析从确定性向非确定性发展是研究的必然趋势，需要用到水位变化时间与速率、汛期洪峰等参数的均值、标准差、分布类型及相关系数。

库水位是影响岸坡稳定的重要随机变量。吴世伟等[9]认为一般水库坝前年水位峰值的分布规律为正态分布或对数正态分布；丁晶等[10]、葛慧等[11]认为坝前年最高水位可由对数正态表征。然而，三峡水库调度使得坝前年最高水位在175 m附近，其随机特性无法用单一的年最高水位的统计特性来描述，而是受水库调度方式及上游来水影响的水位过程线的随机特性。

获取三峡水库多年水位过程线变化过程的统计特性是库区边坡变形预测及非确定性分析准确性的前提。首先根据三峡水库调度特点选取特征参数标定坝前水位过程线，然后运用数理统计方法获取特征参数的分布参数和相关系数、AIC准则识别特征参数的最优边缘分布类型、Bootstrap方法模拟小样本基础上特征参数分布类型及相关系数的统计不确定性。

2 坝前水位过程线特征参数的标定

2.1 水库调度方式

三峡水库按消落期、汛期、回升期、高水位期，综合考虑下游防洪、上游排沙、航道水深、汛末能否蓄满等情况进行库水蓄泄控制。图1给出了2008—2020年共13 a坝前年水位过程线及水库调度初步设计和优化调度方案。图1中，2008年、2009年、2010年水位过程线有其特殊性。2008年开始高水位试验性蓄水，2008和2009年汛后起蓄时间逐步前移，2010年是采用优化调度方案的第一年，之后调度方案不变。

图1 三峡水库调度方案及13条年水位过程线Fig.1 Scheduling of Three Gorges Reservoir and 13 WPLs

初步调度方案如下：11月份至次年1月份水库处于正常蓄水位175 m，1月1日库水开始消落，4月份调度曲线有一高程约160 m的平台，4月末之前水位≥155 m,以保证发电水头和库尾航道水深，5月份库水继续消落，6月上旬之前坝前水位降至145 m。汛期6—9月份，水库维持145 m运行以利于泥沙排出库外，遇大洪水时根据下游防洪需要拦蓄洪水，洪峰过后库水仍降至145 m运行，在泥沙较少的10月份开始蓄水到175 m，达到“蓄清排浑”的效果[12]。2008年水库开始高水位试验性蓄水， 2009年，针对长江上游水库群建设、三峡水库蓄水运用以来入库泥沙量减少、汛末来水量减少等新情况，三峡水库开始实施优化调度方案，将汛后蓄水时间由10月1日提前至9月15日，2010年进一步提前至9月10日[13]。2008年和2009年的汛后起蓄时间与2010年及之后的优化调度方式存在差异，为更好地反映最新调度方式导致的水位过程线的随机特性，采用2010年1月至2021年12月共11 a的水位过程线作为统计数据。

三峡水库调度方式决定了坝前年水位过程线的总体变化趋势，每年的水位过程线具有周期性变化特点。上游来水在时间和空间上的分布不均直接导致库水位年变化过程与运行调度方案之间存在差异。库水位变化在消落期与调度曲线较贴近，调度曲线在主汛期保持汛限水位运行，而实际水位过程线在汛期有洪峰出现。水库在汛期、汛末蓄水期的调度较频繁，水位波动变化较大[14-15]。因此，有必要拟定一些特征参数来标定考虑汛期洪峰和总体调度趋势的年水位过程线。

2.2 标定原理与步骤

2010—2020年的年水位过程线可以在优化调度过程线的基础上加汛期洪峰进行描述。每年1—10月份调度曲线包含6条直线段，汛期洪峰可用2条直线段描述，汛期洪峰把调度方案中的汛限段分为两截。因此，9条直线段便可描述每年1—10月份的水位过程线。

图2给出了库区2013年12月至2014年12月的水位过程线。9条直线段分别为水位消落期ab缓降段、bc消落平台段和cd陡降段，汛限期de峰前平台段、ef洪峰爬升段、fg洪峰下降段、gh峰后平台段，水位回升期hi陡升段和ij缓升段。相应的特征参数共14个，分别为ab缓降时长Tab和速率Kab、bc平台持时Tbc、cd陡降时长Tcd和速率Kcd、峰前和峰后平台的持时Tde和Tgh、洪峰水位Hf、ef洪升时长Tef、fg洪降时长Tfg、hi陡升时长Thi和水位变化速率Khi、ij缓升时长Tij与速率Kij。特征参数的标定便是用9条直线段拟合每年1—10月底的水位过程线，获取图2中点a—j的坐标并计算上述14个特征参数。

图2 水位过程线标定示意图Fig.2 Sketch diagram of calibration for WPL

从图1可以看出多年水位过程线库水缓降向陡降、陡升向缓升的过渡带在160.0～167.5 m，多年的消落期陡降段很贴近，汛限最低水位接近145.0 m且均<145.5 m。点d对应的时间可根据库水首次消落至[145.0,145.5] m的时间近似确定。库水每年10月27日左右蓄至175 m，汛期最大洪峰点f可直接从水位过程线上获取。特征参数的标定还需确定点b,c,e,f,g,h,i对应的时间，是典型的多参数寻优问题。可将年水位过程线S(t)分ad区间确定点b、c，df区间确定点e，fj区间确定点g、h、i，3个区间的目标函数分别为：

(1)

(2)

(3)

式中：Fad、Fdf、Ff j分别为3个区间的拟合误差；Kef和Kfg分别为直线段ef和fg的斜率;Ta—Tj分别为标定点a—j对应的时间；Hbc、Hde、Hgh为直线段bc、de、gh直线段的平均高程。根据以上标定原理，库区水位过程线特征参数的标定步骤如下。

步骤1：在水位过程线6—9月份内获取最大洪峰对应时间Tf和水位Hf， 5—7月份寻找库水首次消落至145.5 m对应的时刻Td。设定Te∈(Td,Tf)，计算直线de高程Hde、直线ef的斜率Kef和截距，按式(2)计算df区间的误差，根据误差最小原则确定Te和Hde。

步骤2：Ta取1月1日，在[Ta,Td]内寻找水位首次达到167.5 m和160.0 m对应的时刻T167.5和T160.0，设定Tb∈[T167.5-30,T167.5+30]，Tc∈[T160.0-30,T160.0+30],消落平台高程Hbc∈[160.0,167.5]。在Tb,Tc和Hbc所属区间均匀布点，对每个搜索点分别计算ab、cd直线段的斜率Kab、Kcd和截距。直线ab的截距根据点a计算，直线cd的截距根据点d计算。设置寻优约束条件Tb≤Tc，按式(1)计算消落期3条直线段的总误差，根据误差最小原则确定Tb、Tc和Hbc。

步骤3：Tj取10月27日，在水位过程线中找出点h对应的大致时间T′，设定Th∈[T′-5,T′+5],Tg∈[Tf,Th]，Ti∈[Th,Tj]，Hi∈[160.0,167.5]。在Th、Tg、Ti,Hi所属区间均匀布点，对每个搜索组合分别计算直线段fg、hi、ij的斜率Kfg、Khi、Kij和相应直线段的截距，计算 [Tg,Th]区间的水位均值Hgh，按式(3)计算fj区间的总误差，根据最小误差确定Th、Tg、Ti、Hi和Hgh。

步骤4：根据点a—j的坐标计算14个特征参数，例如：陡升持时Thi=Ti-Th，速度Khi=(Hi-Hh)/(Ti-Th)。

2.3 标定结果

表1给出了2010—2020年坝前水位过程线特征参数的标定结果，图3给出了实测水位和标定水位过程线分年度对比情况，同时在每年汛期谷底两侧标出了根据标定步骤1和步骤3获取的峰前平台和峰后平台水位值。11个年份的标定结果与水位实测过程线变化趋势一致。标定方法默认运行高水位为175 m，陡降与缓降的分界点在160.0～167.5 m之间，有效避开了2010年1—5月份的“异常”现象。2010年和2019年消落期没有平台出现，2017年消落期平台很窄，仅13 d，2014年和 2020年的bc平台持时在30 d以上。2010年、2012年和2020年汛期洪峰水位较高，分别为160.98、162.87、167.41 m，2015年汛期几乎没有洪峰出现。2018年和2020年汛后期水位未回落至汛限水位附近便进入蓄水阶段。从11条水位过程线的标定结果可知：峰前平台较明显，峰后平台最宽持时达44 d，最窄持时仅4 d，峰后平台最高水位为155.04 m，最低水位为146.10 m，均值为148.19 m。一定程度上反映了汛期来水的时间和强度上的不均匀性。

表1 水位变化特征参数标定结果及分布参数Table 1 Results of characteristic parameters of WPLs

图3 实测水位与标定水位过程线对比Fig.3 Comparison of WPLs between calibration and measurement

3 特征参数的分布参数

特征参数的随机性常用统计特性来表征，包括分布参数(均值、变异系数)、分布类型和相关性。

3.1 特征参数的均值与变异性

表1最后3行给出了14个特征参数的均值、标准差和变异系数。水位消落期，缓降时长和陡降时长均值分别为87.00 d和47.45 d，变异系数分别为0.24和0.23。消落平台持时均值为27.82 d，变异系数为0.64，这与有的年份没有平台、有的年份平台较宽的离散性大相吻合。缓降和陡降消落速率均值分别为14.53 cm/d和36.75 cm/d，变异系数分别为0.22和0.19，说明多年水位过程线在陡降段相对集中，在缓降段相对分散，与图1中水位过程线消落期松散程度一致。

汛限期，峰前和峰后平台的持时均值分别是29.09 d和19.00 d，汛期洪升和消落时长均值分别是18.45 d和13.73 d，变异系数分别为0.80和0.55。洪升时间越短，汛期洪峰水位也不会太高，符合一般规律。汛期4条直线段的持续时间的变异性较大，一定程度上说明汛期上游暴雨出现时间、持续时间和暴雨强度的分布不均及下游防洪需求时间的不均匀性。虽然汛期洪峰的标准差达到5.60 m，但是洪峰水位均值为156.95 m，两者相除后，变异系数值仅为0.04。

水位回升期，陡升段和缓升段的时长均值分别是22.18 d和32.27 d，库水回升速率分别是74.33 cm/d和37.88 cm/d，两个速率的变异系数基本一致。陡升时长的变异系数比缓升时长的变异系数大0.14，这与汛后期降雨时间的不确定性导致有些年份汛期洪峰还未完全消落便开始蓄水，从而压缩了陡升段宽度一致。

3.2 特征参数标定均值线与水位均值线

水位过程线的标定过程中并没有把峰前和峰后两平台水位作为特征参数，而是根据图3中峰前平台和峰后平台11 a水位数据计算其统计特性。采用14个特征参数均值和两平台高程的统计均值无法构造出一条闭环的水位过程线。按照表1中14个特征参数的均值及根据消落期缓降和陡降段、回升期陡升和缓升段的时长和速率乘积与175 m的差值计算得到的汛限期两平台高程绘制的水位过程线称“标定均值线”；根据汛期两平台高程的统计均值修正标定均值线获取的水位过程线称“平台统计均值线”。将水位过程线按公历日期累加求平均可获取“水位均值线”。图4给出了上述3条均值线及11条年水位过程线。3条均值线都能较好地反映坝前水位的变化过程，在回升期的缓升段和消落期十分接近，但是在峰后平台和陡升期差别较大。水位均值线的峰后平台高于标定均值线相应两平台高程，这是因为水位均值线在求平均时存在“削峰填谷”现象，在8月至9月上旬尤其严重。例如2020年8月下旬出现汛期最大洪峰167.41 m，而其它年份该时间段内水位绝大多数在150 m以下，甚至已回落至汛限水位附近，但是11 a数据在该时段平均后水位均值提高了至少1 m。同样的现象在2010年也出现过。平台统计均值线采用的是平台高程的统计均值，与水位均值线的平台高程基本一致，一定程度上也受“削峰填谷”的影响，而标定均值线不仅保留了水位调度曲线的主要特点，而且较好地避开了水位均值线直接平均产生的“削峰填谷”现象。

图4 11 a水位过程线、水位均值线、标定均值线对比Fig.4 Comparison among eleven WPLs, mean water level, and calibration mean results

此外，库区岸坡的变形与库水消落时长和速率相关[16]。汛期洪峰的消落对库岸变形同样重要，“削峰填谷”将直接压缩汛期洪峰的消落时间，如果采用水位均值线或平台统计均值线，将低估汛期库岸变形，高估汛期岸坡稳定性。因此，在库区涉水岸坡的库水周期性震荡影响下的变形预测或确定性分析中运用标定均值线更为合适。涉水岸坡的非确定性分析中除了需要分布参数外，还需要随机变量的分布类型。

4 特征参数的分布类型

4.1 AIC准则

采用AIC准则[17]识别随机变量的最优边缘分布类型，先假定特征参数可能的备选边缘分布，具有最小AIC值的边缘分布通常被认为是拟合实测数据概率分布的最优分布。AIC值计算公式为

(4)

式中：f(xi;p,q)为备选分布的概率密度;xi为样本观测值；k为备选分布参数个数。采用正态、对数正态、极值Ⅰ型和威布尔分布作为备选分布，其概率密度函数和分布参数p和q的换算参照文献[18]。

4.2 最优边缘分布的识别

表2左侧给出了14个特征参数的4种备选概率分布AIC值和据AIC准则识别出的最优边缘分布类型。运用对数正态拟合时，洪峰爬升时长、洪峰水位、洪峰消落时长和陡升时长的拟合效果最好。消落期缓降时长、平台持时和陡降速率、峰后平台持时、缓升时长，当备选分布类型为威布尔时AIC值最小。极值Ⅰ型是缓降速率、陡降时长、峰前平台持时、陡升速率和缓升速率，拟合频率直方图最优的分布类型。

表2 特征参数的分布类型识别结果Table 2 Results of distribution function of characteristic parameters

图5给出了基于11个样本计算的回升期和消落期特征参数的频率直方图和4种备选概率分布拟合曲线。一般地，对数正态和极值Ⅰ型能较好地拟合左偏型频率直方图，而对称型、右偏型频率直方图分别运用正态和威布尔分布拟合效果更好。

图5 消落期和回升期持续时间和库水变化速率直方图Fig.5 Histograms of duration and change rate in falling period and rising period

库水消落期，图5(a)和图5(d)中缓降时长和陡降速率的频率直方图为右偏型，威布尔拟合最佳，与表2左侧采用AIC值最小识别的结果一致；图5(b)和图5(c)中缓降速率和陡降时长的频率直方图为左偏型，陡降时长在41 d频率远大于其它区间，极值 Ⅰ 型在41 d的密度值大于对数正态分布，极值 Ⅰ 型是陡降时长的最优边缘分布；极值 Ⅰ 型在缓降速率直方图左右两侧较对数正态分布更为贴近，是缓降速率的最优分布。

库水回升期，图5(e)至图5(h)中缓升时长为右偏型，运用威布尔分布拟合较好；陡升时长、陡升速率和缓升速率的频率直方图为左偏型，对数正态分布和极值Ⅰ型分布拟合这3个特征参数的直方图时两类分布的概率密度曲线很接近，很难从图中辨认出最优的边缘分布类型。而通过表2中对数正态和极值Ⅰ型AIC值的大小可以较方便地确定相应特征参数的最优分布类型。

4.3 边缘分布类型的Bootstrap模拟

表2左侧识别的分布类型是基于11个样本，具有统计不确定性。Bootstrap 方法是模拟小样本统计不确定性较成熟的方法之一[19-20]。

表2右侧给出了Bootstrap 模拟20 000次的4种备选分布被识别为最优边缘函数的总次数，14个特征参数中没有一个备选分布能被100%选为最优边缘分布。除峰前平台时长和洪峰高程外，其它12个特征参数Bootstrap 模拟的最优边缘分布结果与原始样本AIC值识别结果一致。基于小样本AIC值确定的汛期洪峰的最优分布类型为对数正态，而Bootstrap 模拟中极值Ⅰ型和对数正态分布被选为最优分布的次数分别是6 757次和5 513次，极值Ⅰ型分布是考虑统计不确定性后的汛期洪峰的最优边缘分布。11个样本基础上峰前平台持时，极值Ⅰ型AIC值为82.62，略小于对数正态的82.66，因而认为极值Ⅰ型是最优边缘分布函数，而Bootstrap 模拟中对数正态被选为最优的次数比极值Ⅰ型多1 342次，对数正态是峰前平台持时Bootstrap 模拟的最优边缘分布。

5 特征参数间的相关系数

每年库区水位从175 m消落至145 m耗时约160 d，从145 m回升到175 m耗时50 d左右。不同的库水消落回升时长与速率组合之间必然存在相关性。

5.1 相关系数的计算方法

可靠度计算中常用的相关系数主要有Pearson线性相关系数ρ和Kendall秩相关系数τ。Pearson线性相关系数在相关非正态随机变量非线性变换前后值会发生变化，而Kendall秩相关系数非线性和线性单调变换后值保持不变。两者均可以根据标定出的特征参数进行计算，公式如下：

(5)

(6)

5.2 库区水位过程线特征参数的相关性

运用式(5)和式(6)分别计算表1中相邻2个特征参数之间的相关系数，其中相关性较大的6组列于表3。缓降和陡降段的持时和速率、陡升和缓升持时和速率均为负相关，即要消落或回升相同的水位差，持续时长越短，需要的库水上升和下降速率越大。汛期洪峰与洪峰下降天数成正比，即洪峰越高，洪水消去时间越长。消落期平台持时与缓降速率正相关，即缓降速率越大，从175 m消落到160 m左右的时间越短，160～165 m高程附近水位“逗留”的机动时间可长一些，上述现象符合库水调度的一般规律。

表3 特征参数间的相关系数Table 3 Coefficients of correlation among characteristic parameters

在给出的6组相关系数中库水位陡升时长和速率的负相关性最大，Pearson相关系数为-0.85，陡降时长和速率的Pearson相关系数次之，为-0.78；库水位缓降时长和速率、缓升时长与速率的负相关性最小，均为-0.67，但Kendall秩相关系数分别为-0.31和-0.70。消落期缓降速率和平台天数，汛期洪峰与洪峰下降天数的Pearson相关系数均为0.70，Kendall秩相关系数分别为0.46和0.60。

基于11个样本数据的特征参数相关系数之间同样存在统计不确定性。

5.3 相关性的Bootstrap模拟

在Bootstrap模拟中，计算20 000个子样本6组特征参数的相关系数，然后统计子样本相关系数的均值和变异系数，计算结果列于表3右侧。Bootstrap模拟相关系数均值结果与基于11个小样本计算得到的Pearson相关系数和Kendall秩相关系数基本一致，说明Bootstrap模拟继承了原始样本特征参数之间的相关特性。Bootstrap模拟的Pearson相关系数均值误差比Kendall秩相关系数均值误差小，子样本Pearson相关系数的变异系数小于Kendall秩相关系数的变异系数，变异系数绝大多数小于0.20。缓降时长与缓降速率之间相关性的误差是6组相关性中最大的，Pearson相关系数均值误差为2.90%，Kendall秩相关系数均值的误差达18.42%，Bootstrap模拟的变异系数也最大，分别为0.25和0.56。

6 结 论

本文对2010—2020年三峡坝前水位过程线进行特征参数的标定与统计，得到以下结论：

(1)根据三峡水库调度特点及考虑汛期洪峰，将每年1月1日至10月27日水位过程线分9条直线段进行拟合，消落期、汛限期、回升期分别用持续时长、高程、水位上升或下降速率等14个特征参数进行标定，标定均值过程线能有效避开多年水位平均线在汛限期和陡升期的削峰填谷现象，比水位均值线更适合作为库水周期性震荡影响下库岸边坡变形预测及稳定分析的典型水位过程线。

(2)基于11 a小样本特征参数的统计结果表明，库水消落期，缓降时长和水位消落速率分别服从威布尔和极值Ⅰ型分布，均值和变异系数分别为87.00 d、0.24和14.53 cm/d、0.22，Pearson相关系数为-0.67。陡降时长和水位下降速率分别服从极值Ⅰ型和威布尔分布，均值和变异系数分别为47.45 d、0.23和36.75 cm/d、0.19，Pearson相关系数为-0.78。库水回升期，陡升时长和水位上升速率分别服从对数正态和极值Ⅰ型分布，均值和变异系数分别为22.18 d、0.38和74.33 cm/d、0.29，Pearson相关系数为-0.85；缓升段时长和水位回升速率分别服从威布尔和极值Ⅰ型分布，均值和变异系数分别为32.27 d、0.24和37.88 cm/d、0.28，Pearson相关系数为-0.67。

(3)水位过程线特征参数分布类型和相关系数的统计不确定性模拟表明汛限期峰前平台持时和洪峰高程的分布类型应修正为对数正态和极值Ⅰ型分布；缓降时长与速度、缓降速度与消落期平台之间的Pearson相关系数均值和Kendall秩相关系数均值的变异系数较大。在库区岸坡稳定分析的风险分析中应引起重视。

(4)三峡库区范围很大，本文研究的是坝前水位过程线的统计特性，如何建立库区涉水滑坡所在地与坝前水位过程线特征参数之间的联系是值得研究的一个问题；如何考虑水位过程线特征参数分布类型及相关性的统计不确定性带来的库岸风险区间估计是今后需要研究的另一个问题。