2.1　明渠输水水力过渡过程模拟

2.1.1　复杂输水工程的常态输水过程

2.1.1.1　明渠水力计算基本方程

明渠输水系统恒定流计算和非恒定流计算的基本方程为圣维南方程。明渠水流通常指侧面和底面为固壁而上表面与大气接触的水流，如渠道、河道以及横断面未充满的管道中的水流等。当明渠中水流流速和水深等水力要素不随时间变化时，称为恒定流;当明渠中水流的流速和水深等随时间变化时，称为非恒定流。明渠非恒定流计算已经成为计算水力学的一个重要组成部分，虽然其基本理论和技术已经成熟，但仍有诸多问题尚需进一步研究，如临界流问题、堰闸孔流边界的过渡问题和明满流交替问题等。

(1)恒定流基本方程。恒定流是输水渠道运行中非常重要的一种水流流态。运行调度的实质就是通过调整分水闸分水流量、泵站抽水流量使输水系统从初始输水流量下的恒定流状态平稳过渡到目标流量下的恒定流状态，恒定流计算模型应能计算不同分水流量、各泵不同抽水组合情况下明渠恒定流状态下的水面线，即渠道稳定输水时的水流状态。在稳定流情况下，圣维南方程中所有关于t的微分项均为零，见式(2.1-1):

式中:x为空间坐标，m;Q为明渠输水流量，m3/s;q为单位长度渠道上的侧向入流量，m2/s;A为断面过流面积，m2;h为明渠断面水深，m;g为重力加速度，m/s2;Sf、S0分别为摩阻坡度和渠道底坡。其中，Sf采用曼宁公式近似:

式中:K为流量模数。

(2)非恒定流基本方程。非恒定流计算采用一维圣维南方程组，方程组由连续性方程和动量方程组成:

式中:t为时间坐标;其余参数意义同恒定流基本方程。

(3)倒虹吸基本方程。对于长度较短的倒虹吸，可以对水击波速进行计算，从而得出明渠水波在倒虹吸中的传播时间。水击波速的计算公式为

式中:c为倒虹吸中的水击波速;K为水的体积弹性系数;σ为水在20℃下的密度;E为倒虹吸壁的弹性系数;D为倒虹吸的等效半径;e为倒虹吸的厚度。其中，K=2.39×109;σ=1×103;E在本书中取混凝土的弹性系数，为30×109;e取0.5m。计算得出:c=953.47m/s，水波在倒虹吸中的传播时间约为0.31s。由于明渠计算时间步长Δt远远大于水波在倒虹吸中的传播时间，因此倒虹吸可当作水头损失处理。在计算有压非恒定流时，两节点间的水位差为进口渐变段h1、进口闸室段h2、出口渐变段h3、出口闸室段h4的局部损失以及洞身沿程水头损失之和。

1)进口渐变段:

式中:v1为进口渐变段的进口流速，可以采用上游明渠中的流速来代替;v为进口渐变段的出口流速，可以采用倒虹吸直管段的流速，在已知倒虹吸过流量的条件下，可以采用恒定流公式求得。

2)出口渐变段:

式中:v2为出口渐变段的出口流速，可以采用下游明渠中的流速来代替;v为出口渐变段的进口流速，可以采用洞身的流速，在已知洞身过流量的条件下，可以采用均匀流公式求得。

3)进口闸室段:

4)出口闸室段:

式中:v为闸室内的流速，可近似为倒虹吸直管段的流速，在已知倒虹吸过流量的条件下，可以采用均匀流公式求得。

结合均匀流方程，满足以下连接条件:

式中:L为倒虹吸的长度。

2.1.1.2　方程离散

2.1.1.2.1　Preissmann格式

圣维南方程组的数值解法主要有有限差分法、特征线法、有限体积法和有限元法等。通过对各种解法进行比较，Preissmann格式具有收敛性好、稳定性和精度高的特点，采用的Preissmann格式为

时间方向离散:

连续方程空间方向离散:

动量方程空间方向离散:

式中:θ为连续方程空间方向离散权重;ψ为时间方向权重;ϕ为动量方程空间方向离散权重。

连续方程离散为

动量方程离散为

对Preissmann格式的稳定性研究表明，当水流是缓流(Fr<1)时，因两个根λ1>0和λ2<0，有Cr=λ1>0和Cr=λ2<0。

因Cr>0和Cr<0同时存在，所以两个离散方程不会同时满足无条件的稳定状态。Cr<0时是有条件的，其条件为

如果权重参数选择不合理，则在计算过程中会出现波动现象，从而可能会出现不收敛，因此在使用Preissmann格式时，应该对两个方程采用不同的权重，这样可以保证计算的无条件稳定。

2.1.1.2.2　离散后方程的线性化

离散后的连续方程和动量方程是非线性的，这样就需循环迭代求解离散后的两个方程。目前有直接求解h、Q和求解d h、d Q两种求解离散非线性方程的方法，这里采用后者。在循环迭代的过程中，用到下列关系式:

式中:上标*代表上一个循环步的变量值;ΔA、Δh、ΔQ分别为过流面积、水深和流量的增量;B为水面宽。

将式(2.1-18)代入式(2.1-15)，有

其中

为了对动量方程进行线性化，要用到以下关系式:

将式(2.1-20)代入式(2.1-16)，有

其中

同理，对于倒虹吸基本方程，可以对式(2.1-9)和式(2.1-10)线性化处理得

其中

把式(2.1-19)和式(2.1-21)化简，其系数为

a2j=0

2.1.1.3　线性方程的求解算法

式(2.1-19)和式(2.1-21)组成的系数矩阵再加上内外边界条件后，可以用高效率的双扫描法求解。

为了方便用矩阵形式表示，将式(2.1-19)和式(2.1-21)改写成下述形式:

式中各微增量的系数与式中的系数一一对应，i=1，2，…，n。

渠道线性化方程加上内外边界条件后的矩阵形式为

式中:A为五对角带状矩阵;X为未知增量Δh、ΔQ所组成的列向量。

A、X、D内部矩阵形式可表示为

假设系数b1≠0，如上游为水位边界条件，则可以采用消元法把式(2.1-26)变换为

其中

式中元素由下述递推公式计算:

线性方程式(2.1-27)中的矩阵B是对角线上元素全部为0的上三角矩阵，可以采用回代过程递推求解:

当系数b1=0时，取列向量X=[ΔQ1 Δh1 … ΔQ2n+2 Δh2n+2]，与此同时调整矩阵A每一列上元素的位置，因而上述求解方法仍然适用。

由此可见，离散方程求解由一个消元过程和一个回代过程组成。在已知上一时刻t0明渠各断面流量和水深的条件下，采用下述过程求解t0+Δt时刻的明渠各断面流量和水深。

(1)假设t0+Δt时刻明渠各断面的流量和水深的初始值。

(2)计算矩阵A和列向量D的非零元素。

(3)推求线性方程式(2.1-27)中的X。

(4)判断|xi|≤ε，i=1，2，…，2n+2。如果所有条件都成立，则用求解出的流量和水深作为t0+Δt时刻的解。如果有任何一个条件不成立，则用hi+Δhi和Qi+ΔQi分别代替hi和Qi，重复步骤(2)～(4)。

2.1.2　冰期输水过程水力学模拟

对于大型渠道，根治冰害最适宜的方法是通过对渠道水流的调节、管理来控制冰花的生成及冰盖的形成。渠道中的泵站、闸或坝等建筑物可被看作壅水建筑物，通过水流调节促使积聚冰盖首先在建筑物处形成并逐渐发展增长。密云水库调蓄工程可以看作串联沿程分流(十三陵水库分水口、李史山小中河分水口)的渠系，控制方程以时间和流动距离为自变量，以水深、断面平均流速、断面平均水温、断面平均水中冰浓度及浮冰浓度作为基本因变量，根据质量守恒、动量守恒和能量守恒原理，建立明渠冰期基本控制方程组。

明渠冰期水力学模拟模型包括水流模型、热力模型和冰冻模型。水流模型用于计算渠道中的流场和其他水力要素;热力模型用于计算水体热交换、水温分布和降温过程;冰冻模型用于模拟冰冻的产生、输运发展和消融过程。这3个过程是相互联系、相互制约的:水力条件影响热力交换和冰冻过程，热力条件决定冰冻过程，冰冻条件又反过来影响水力条件和热力条件。

水流模型用于计算渠道中的流场和其他水力要素，既能解敞流期的流量与水位过程，也能解在有冰盖情况下的流量与水位过程，同时该模型适用于棱柱形渠道和非棱柱形渠道。模型采用有限差分的计算方法求解。

热力模型包括热交换模型、水流温度场模型。

(1)热交换模型。冬季河槽的温度状态取决于水(冰)体与周围环境之间不停发生的热交换。这一热交换过程决定了渠道水流初期的冰情和以后的冰情发展。

(2)水流温度场模型。水温是冰情研究中的重要组成部分，水面初冰是水体表面薄层的水温降低到0℃或0℃以下的结果;水面封冰是水体向水表面传热与水表面向大气散热互相平衡的结果，而上述这两种热量的数值都与水温有关;水面终冰是冰厚融消到零的结果，而且冰厚融消的计算中也需要知道水温，因此水温是冰情发生、发展和消失过程中起着重要作用的物理量，其计算十分重要。本书采用一维时均水流非稳定温度的数学模型并应用有限差分方法进行数值计算，从而得到水温随时间、空间的变化规律。

冰冻模型包括浮冰、水中冰、岸冰计算模型，结冰期连续冰盖形成的动态数学模型和冰盖厚度增长与消融的动态模型。

2.1.2.1　水流模型

2.1.2.1.1　冰盖下水面线数学模型式中:A为水流的过流面积;x、t分别为距离和时间;B为水面宽度;h为水深;S0为底坡;hi为冰盖厚度;ρi为冰的密度;ρ为水的密度;v为平均流速;S为摩阻坡度，是作用于单位重量液体上的阻力;为由于非棱柱形渠道

由分层流理论可得到以水深和流速表示的基本方程组。

连续性方程:

运动性方程:

f展宽面而增加的面积，对于棱柱形渠道，该项为零。

联解方程组并使其符合初始条件和边界条件，就可得出冰盖下水流的流速和水深随流程和时间的变化关系。

2.1.2.1.2　冰期水流流动阻力的确定

渠道中形成冰盖，这时的阻力应由两部分组成:水流与床面的摩擦切应力和水流与冰盖下表面的摩擦切应力，可用曼宁公式来表示式(2.1-30)中的阻力项，即

式中:R为封冻时的水力半径;nc为复合糙率系数，反映渠道糙率nb和冰盖下表面糙率ni的综合糙率。

结冰渠道的糙率是反映渠道阻力大小的主要参数，特别是对于大型渠道，冰盖阻力所产生的影响比起冰盖厚度的影响要大得多。

2.1.2.2　热力模型

2.1.2.2.1　水流温度场计算

水温是冰情研究中重要的组成部分，水面初冰是水体表面薄层的水温降低到0℃的结果;水面封冻是水体向水表面传热和水表面向大气散热相平衡的结果，而上述这两种热量的数值均与水温有关;水面终冰是冰厚消融到零的结果，而在冰厚消融的计算中也需要知道水温。因此，水温在冰情的发生、发展和消失的过程中一直起着重要的作用，它的计算是十分重要的。

对于输水水流温度场计算，一般采用一维计算，方程为

式中:Cp为水的比热;ρ为水的密度;T为水温;Q为流量;Ex为综合扩散系数，也称为混合系数;A为过流断面面积;B为水面宽度;∑S为水流与周围环境的热交换通量。

给定初始条件和边界条件，求解式(2.1-32)可得到水流温度场。

2.1.2.2.2　热交换过程

冬季渠道的温度状态取决于水体与周围环境不停发生的热交换。影响水流热平衡的因素很多，但从其作用的性质可分为两种，即增热因素和失热因素。一般以使水体增热为正。各水流热平衡因素常用单位水面上、单位时间内热的交换量来度量，单位为W/m2。

渠道中水体热交换包括水与大气热交换S1、水与渠底热交换S2、水与冰盖热交换S3。其中，S1在水体的热量与周围环境介质的交换中起着主导作用。在严寒的冬季，无论是夜间还是白天，水体的吸入热量总是小于实际向外散失的热量，一般为负值，故

(1)水与大气的热交换。单位水面与大气的热交换S1，与气象条件有着密切的关系，以水表面受热为正，可由式(2.1-34)表示:

式中:Ss为太阳短波净辐射通量;Sa为大气长波辐射通量;Se为水表面蒸发(或凝结)热通量，表面受热为正;Sc为对热流损失，表面受热为正。

1)太阳短波净辐射通量Ss。太阳能以短波辐射的形式通过地球大气层直接传到水面，称为短波太阳辐射(或称为太阳总辐射)。这是太阳辐射的一部分，不能被大气所吸收，表现为可见光。该项为表面受热，符号为正。短波太阳辐射包括直射Ss1和射散Ss2两部分。

投射到水面上的总辐射并不能被水体全面吸收，有一部分要被反射回去，称为反射辐射Ss3。因此，水表面所吸收的，或者说太阳短波净辐射通量Ss为

太阳的直射和散射一般来讲与当地的天文辐射量、大气透明程度和天气晴朗程度有关。从我国许多气象台的实测资料来看，这部分能量与相对日照数的关系较密切。

2)大气长波辐射通量Sa。长波辐射，或称为有效辐射，是指渠水上面的大气辐射。其值与地理纬度关系不大，而与气温、云量、湿度有较密切的关系。长波辐射由水(冰)面发出的长波辐射和水(冰)表面吸收的净大气热辐射两部分组成:

式中:Sa1为由水面发出的水体长波福射，为表面失热;Sa2为水表面吸收的净大气长波福射，该项表面受热，符号为正。

3)水表面蒸发(或凝结)热通量Se。

其中

式中:L为潜热;γ为水的容重;E为单位水面的体积蒸发率;v为地面2m处风速;es为水表面温度时的饱和水汽压;ea为空气中的实际水汽压。

4)对流热损失Sc。水体与大气之间存在温差，可通过热传导进行能量交换，以水面受热为正。对流热损失Sc主要取决于水体和空气两者之间的温差及紊动交换，是水温与气温的差值和风的函数。采用Rhnsha-Donchenko公式计算:

当有冰盖覆盖时，考虑冰盖影响，引入系数C:

(2)水流和渠底之间的热交换。水流和渠底的热变换，一般包括两部分:渠水与地下水之间的热交换和水与底部土壤之间的热传导。

1)渠水与地下水之间的热交换。对渠水与地下水之间的热交换，目前还没有办法进行处理。但当地下水水位较深时，这一热交换是很小的，从收集到的资料来看，均没有考虑河水与地下水之间的热交换。

2)水与底部土壤之间的热传导。从水体热平衡诸因素的影响程度来看，由土壤进入水中的热量是不大的，故常常可忽略不计，特别是明流时，水和渠底热交换的作用很小。但当渠道封冻后，则需考虑水与底部土壤之间的热传导:

式中:a为水流层与土壤的热交换系数;Tw为水流底层的温度;Tb为土壤下层的温度。

Tb可按一维热传导方程求解:

式中:Kb为土壤的热传导系数。

(3)水与冰盖之间的热交换。当水体表面形成冰盖后，水与大气之间的热交换变为水与冰盖的热交换。水流与冰盖之间的热交换，对冰盖形成、冰盖厚度及冰盖的消融有很大影响。采用Ashton于1982年提出的水和冰盖线性热交换的计算公式，即

其中

式中:Tb为冰盖底部温度;hwi为冰与水的热交换系数;v为平均流速;h为水深。

2.1.2.3　冰冻模型

2.1.2.3.1　浮冰、水中冰、岸冰的形成及演变模型

(1)浮冰、水中冰的初始形成。浮冰的形成决定于水面温度和紊动强度。Matousek在1984年通过对Ohio河的现场观测，发现浮冰和水中冰的形成与水面温度、水深平均温度、当地流速、风速、水面宽度等因素有关，并根据观测结果提出了以下浮冰和水中冰的判别法则。

1)Tws≥0℃，不出现冰情。

2)Twd≥0℃:如果<Tws<0℃和vb≥vz，将出现表面浮冰;如果将出现大块表面浮冰。

3)Tws≤0℃，将出现水中冰。

注:为临界水面温度，Tws为水面温度，Twd为水深平均温度。

(2)浮冰、水中冰的沿程分布。当水温达到0℃左右时，如继续失去热量，冰就开始形成。此时称为冰期起始时刻，水体中出现浮冰和水内冰。

浮冰浓度分布可由以下一维对流-扩散方程表示为

水中冰浓度分布为

其中

式中:Cs为表面流冰浓度;Cc为水中冰浓度;ρi为冰的密度;Li为结冰潜热;α为水内冰转变为浮冰的比例系数;ui为浮力作用所产生的冰粒上浮速度;vz为由于水流紊动作用所引起的垂向速度分量。

(3)静状冰(岸冰)的形成与发展。静状冰是指生成于渠道缓流区的岸冰，是由冰晶在过冷却表层形成且停滞在水面以较为缓慢的方式生成。其生成及消融受热力因素影响较大。由于它一旦形成后，就会减小敞露水面，因此需要考虑岸冰的影响。本书中，为了一致性，采用Matousek准则来模拟岸冰的形成。

岸冰形成后，由于水面浮冰(流冰)积聚，岸冰还可能沿横向发展，其横向增长取决于浮冰块与岸冰接触时的稳定性。Michel等通过对圣安纳河道中的冰情观测，提出了以下岸冰增长的计算公式:

式中:vc为浮冰黏附于岸冰的最大允许流速;ΔW为给定时段内岸冰宽度增长量;∑S为给定时段内、单位面积上的热交换(热损失);N为流冰的面密度。

vc值的确定与水流条件和流冰情况有密切的关系，取决于浮冰块与岸冰接触时稳定情况。浮冰块与岸冰接触时稳定与否则取决于作用在两者的合力。

2.1.2.3.2　冰盖动态模型

(1)结冰期连续冰盖形成的动态数学模型。初冬，各种冰正在形成并沿河往下游输移，其中大部分漂浮到水面。随着气温的不断下降，渠道产冰量逐渐增多，流冰密度增大。当密集的流冰在急弯、浅滩、建筑物、断面束窄处(包括岸冰阻截)等渠段发生卡堵时，自上游不断向下输运的冰花即从此平铺上溯，于是该渠段形成封冻冰盖，且不断向上游发展。初始冰盖厚度随形成条件而异。当流冰出现时，由于表层流冰堵塞，或者由于天然或人工障碍(如闸门、拦冰索或冰桥)，在适当的流动条件下，就可能出现初始冰盖。初始冰盖一旦形成，便将通过积聚上游的来冰逐渐向上游推进。其推进速度取决于上游来冰情况、冰盖体前缘厚度及水流条件等。当流速缓慢时将形成光滑的冰盖;而在流速较大的渠段，由于来冰在冰盖前缘断面推挤、下潜和附着，使得形成的冰盖犬牙交错，起伏不平。

冰盖推进方式取决于上游来冰情况和水流条件，渠道中冰盖向上游发展的过程中，有冰盖并置推进(平封，冰块拼接)、水力增厚推进(窄封)、机械增厚推进(即紧冰型推进，宽封)3种基本形式。

1)冰盖并置推进(平封)。在流速较为缓慢的渠段，由于水流作用较小，浮冰块通常只能形成简单堆积。冰盖前缘厚度与上游输移下来的冰花厚度相同，冰盖前缘以等同于单个浮冰块的厚度，以并置积聚的形式向上游推进。平封仅发生在流速较低的情形。

初始冰盖厚度取决于前缘来流推力。根据Pariset等于1966年提出的拟静力冰盖理论，确定初始冰盖厚度。

满足冰盖并置推进的水流条件，采用Pariset和Hausser于1961年提出的向下旋转并潜没的临界弗劳德数，或称为第一临界弗劳德数Fr1判别:

其中

ej=ep+(1-ep)ec

式中:为冰块的形状系数，在0.66～1.3之间取值;ti为冰块厚度(冰盖平衡的初始厚度);li为冰块长度;h为水深;v为上游前缘流速;ρi为冰块的密度;ρ为水的密度;ej为堆积冰盖体的总空隙率;ec为单个流冰块的空隙率(约为0.4);ep为流冰块间的堆积体空隙率。

当冰盖前缘的水流弗劳德数小于第一临界弗劳德数Fr1时，冰盖以聚集形式向上游推进，积聚冰盖厚度存在某一极限值。

如果水流弗劳德数超过第一临界弗劳德数Fr1，则冰盖前缘将停止向上游发展，上游来冰将潜入冰盖下，并堆积在冰盖底部，形成冰塞。随着冰塞体厚度不断增加，上游水位壅高，同时，流冰潜入冰塞下所需的弗劳德数增大。当其厚度达到一定值后，上游流冰不再下潜，此时，冰层前缘继续向上游推进。

2)水力增厚推进(窄封)。当冰盖前缘的水流弗劳德数超过第一临界弗劳德数Fr1时，单一冰块厚度的并置推进将不可能维持。这时冰盖将以水力增厚的方式向前推进。上游下来的流冰将会在冰盖前缘附近下潜，潜入冰盖底部并在附近堆积，形成大于单个流冰块厚度的冰盖层，导致冰盖厚度增加、前缘水位壅高、底部水流分离。

在时段Δt内，长度为Δx的河段上的冰盖层厚度变化为

浮冰块潜入冰盖底部所需的水流流速随冰盖厚度的增加而增加。因此，当冰盖前缘来流量一定时，随着冰盖前缘厚度的逐渐增加，水流弗劳德数逐渐下降。当小于第一临界弗劳德数Fr1时，冰盖前缘继续向前推进。

3)紧冰型推进(宽封)。在宽式渠道中，作用于冰盖的纵向作用力可能会超过岸边阻力，使冰盖破坏，然后又插堵堆积，冰盖增厚，直到达到平衡为止。平衡条件为

其中

式中:f为冰盖纵向力;ti为水流作用于冰盖下方的剪应力;τg为沿冰盖方向的自重分量;τai为沿冰盖的风应力;ρa为空气的密度;va为风速;θa为风速与水流下游方向的夹角;Cai为与表面粗糙度有关的阻力系数(约为1.55×10-3);μ1为岸边阻力系数;k1为侧向推挤系数(0.342);τc为岸边阻力的凝聚力分量;B为河宽。

4)初始冰盖空隙率。初始冰盖的空隙率随冰盖的形成方式和表层流冰的性质而异。对于由厚度为hi的固结冰壳和冰壳下方厚度为hf的淤积冰絮所组成的流冰来说，空隙率为

流冰之间的空隙率，当最大表面浓度为Cmax时，ec=1-Cmax，单层冰盖的整体空隙率为

对于多层冰盖，由于较大驱动力的挤压作用，空隙率略小于单层冰盖，但目前尚无足够资料来加以确定，形式上可以表示为

式中:αp为压缩因子(<1)。

5)冰盖下输冰和潜冰塞演变。潜冰塞是水内屑冰或上游来冰在冰盖下面的堆积体，它常常形成于较大流速敞流段下游冰盖的下表面上。在此模型中，采用输冰能力的理论来确定冰盖下潜冰的堆积。

根据现场观测和室内试验，Shen和Wang(1995)提出了冰盖下输冰能力公式:

其中

式中:ϕ为无量纲输冰能力;Θ为无量纲水流强度;vWi为冰粒上浮速度;qi为单宽体积输冰流量;τi为潜冰塞下表面上的剪应力;v*，i为潜冰塞下表面上的剪切速度;ΘC为无量纲临街剪切力，ΘC=0.041;F为下沉速度系数，对于球形冰粒，其值近似等于1.0;dn为冰粒的标称粒径。

由面冰疏密度计算结果可以确定冰盖推进前缘断面处进入封冻渠段的单宽面冰流量。

冰盖底部的冰流量连续方程为

其中

式中:tf为潜冰塞厚度;qi为冰盖下表面上的单宽冰流量;eu为潜冰塞空隙率;为与悬移冰层之间的净交换冰流量。

当冰盖下表面上的挟冰超过其挟冰能力时，屑冰将在冰盖下表面上堆积;对于冲蚀情形，从堆积冰体上冲蚀下来的冰块将加入下面冰流动。堆积和冲蚀均受可供水量的限制。

在时间步长Δt内，长度为Δx的渠段里，水内屑冰堆积或冲蚀引起的潜冰厚度的平均变化可由式(2.1-59)计算:

式中:上标u、d分别代表上、下游断面。

(2)冰盖厚度增长与消融的动态模型。渠道中形成连续冰盖后，水体通过冰盖传导与大气交换热量，此时水体失热速度明显减缓，水体失热主要体现在冰厚的增长，也就是说冰盖厚度将随着冰和大气的热交换程度而变化。此外，由于水流失去了过冷却的条件，所以封冻渠段不再产生水内冰。

将渠道作为由大气、冰、水和河床组成的一个完整封闭系统，根据系统热交换过程理论，考虑详细的能量交换过程研制冰层厚度变化模型。

冰盖内一维热传导方程为

式中:T为冰盖内的温度;ρi为冰的密度;Ci为冰的比热;Ki为冰的热传导系数;A(z，t)为冰盖内的热源。

冰盖上表面的边界条件为

式中:∑S为气与冰表面的热通量损失;Li为冰的潜热;hi为冰盖厚度。

冰盖下表面的边界条件为

式中:Swi为从水到冰的热通量。

式(2.1-61)和式(2.1-62)构成了非线性边值问题，由差分方法求解。

2.1.3　渠道渗漏损失计算技术

渠道渗漏损失是指在渠道输配水过程中，通过渠床表面渗入渠床土体而损失的水量。正确估算渠道的渗漏损失量是计算渠道水利用系数的关键。由于影响渠道渗漏的因素甚多，目前尚无通用的渠道渗漏损失计算公式。一些经验公式是基于渠道现场试验及统计分析得出的，如戴维斯和威尔逊总结的衬砌渠道计算公式，印度采用的对土质渠床损失水量计算的经验公式等。国内广泛采用的是考斯加可夫公式，一些学者还为提高计算精度对公式进行了改进。

计算渠道渗漏损失的经验公式一般是在对实测资料进行统计分析的基础上，选取某些主要影响因素作为变量进行推导或率定的。影响渠道渗漏的主要因素可分为以下几类。

(1)渠床土壤。不同土壤的渗透性不同，土壤越黏重，土粒间孔隙就越小，水流渗漏的阻力越大、水力传导度越小，因此重质土壤的渠床渗漏量小于轻质土壤。渠道渗漏量公式中常用渗透系数反映不同土壤的渗透性。

(2)断面形式及水力特性。主要包括湿周、水深、流速和流量等水力参数。湿周越大，渗漏面积越大，渗流量越大;水深与渠床土壤水势梯度正相关，因而与渗漏通量正相关。同时渠道水深也与湿周正相关，进而影响渗漏量。流速对渠道渗漏的影响主要体现在两个方面:当流量一定时，流速越大，渠道过水断面越小，水流渗漏面积也越小，渗流量随之变小;流速越大，相同水量的水流停留时间越短，则渗漏损失越小。流量对渠道渗漏也有较大影响，对于相同断面形式的渠道而言，渠道流量越大，水深和渠道渗漏面积越大，因此渠道渗漏量越大。实际上这些因素彼此关联，相互影响，实际中难以进行明确区分，均属于渠道的断面形式及水力特性，因此，不同的公式所采用的参数有所不同。

(3)地下水埋深。渠道在输水过程中的渗漏可分为两个不同的阶段，即自由渗漏和顶托渗漏。自由渗漏的损失量主要取决于渠道的渠床土壤性质和渠道断面尺寸。顶托渗漏的损失量基本上是由地下水埋深决定的，受渠道断面尺寸影响较小。当渠道地下水水位较高时，渠道水受地下水顶托作用较大，减少了渠道水下渗的水力梯度和储水空间，对渠道水的下渗起到抑制作用，理论计算时，采用地下水水位与渠内水位的差值为变量估算渗漏量，但由于地下水水位的时空不稳定性，资料难以获取，实际计算时，先计算自由渗漏，对有地下水水位顶托的情况再采用折减系数修正。

(4)衬砌条件。不同衬砌条件的防渗效果不同，与普通土渠相比，一般采用黏土夯实、混凝土衬砌和塑料薄膜衬砌的渠道渗漏损失量能分别减少45%、70%和80%左右。理论计算时，需要测试衬砌材料的厚度、渗透系数和土壤的分层厚度、渗透系数，并结合地下水埋深和渠道的断面形式等计算渗漏量。因所需的资料往往难以获取，故在实际计算时，常常先根据渠床土质计算无衬砌时的渗漏量，再乘以衬砌折减系数得到衬砌条件下的渗漏量。

京密引水渠反向输水渠道在输水过程中渠道水量损失主要包括蒸发、充水期水量损失、渗漏损失等。但是渗漏损失为主要水量损失方式。京密引水渠反向输水路线主要由明渠、渐变段、倒虹吸等组成，明渠段较长，水量损失大，倒虹吸尺度小，水量损失小。因而在计算渠道水量损失时有以下5个假设条件。

(1)水量损失均集中在渗漏损失，忽略水量蒸发损失。

(2)水量损失全部在明渠段，倒虹吸不考虑水量损失。

(3)水量损失与流量没有直接关系，只与水位高低和渠道过流断面大小有关。

(4)渠段水量损失是沿程均匀分布的，即将水量损失值平均到各计算段内。

(5)明渠段渗漏损失符合达西定律，即输水渠道渠底高程高于地下外水位。

通过对实测资料进行统计分析，明渠段渗漏损失符合达西定律，公式为

其中

式中:v为过流断面渗透流速，m/s;K为土壤渗透系数，cm/s;J为渗透坡降。

京密引水渠道过流断面示意如图2.1-1所示，主体过流断面为梯形断面，上面部分为直墙段，边坡为m，底宽为b，水深为h。

通过工程设计资料可知，京密引水渠道为预制板衬砌，梯形断面部分为现浇混凝土衬砌，直墙段为浆砌石衬砌。为了简化处理，将梯形断面部分和直墙段作为一个整体分析，假设整体综合渗透系数为K。假设渠道渠底高于地下水水位，渠道渠底和渠道边坡混凝土衬砌厚度δ1=20cm，渗透坡降直墙段浆砌石衬砌厚度δ2=15cm，渗透坡降

对于长度为1m的渠道，底面渗透流量为

水深为x处，边坡上单位面积渗透流量为

单侧边坡渗透流量为

对于相邻泵站间的渠段，假设有n个计算断面，桩号依次为S1，S2，…，Sn，渗透流量依次为Q1，Q2，…，Qn，则总渗透流量为

断面总渗透流量为假设各计算渠段综合渗透系数相同，即

Ki-1=Ki=K

令

则

式中:K为渗透系数，cm/s;X为根据渠段水位、几何参数计算值，m3/cm。

密云水库调蓄工程团城湖—怀柔水库段共有6座泵站，分为5大段，由于各个渠段的衬砌修复水平不一致，渠道长度、泵站运行水位也不一致，因而实际水量损失计算时需要针对不同的渠段分别考虑这些因素对输水损失的影响，5个渠段的系数分别为α1，α2，…，α5，各渠段水量损失公式为

水量损失较大时取α>1，较小时取0<α<1，具体取值根据实测数据进行率定。

2.1.4　明渠输水实时校正技术

由于实测数据误差或模型参数误差，常导致水动力学模型计算结果与真值差异过大，需对模型进行实时校正。水动力学模型的校正包含状态量校正和参数在线估计。状态量校正需要根据渠道若干个断面出现的水位、流量“新息”(即当前时刻预报值与实测值之差)全面合理地对各断面的计算水位、流量做校正;参数估计需要实现糙率的自动调整。在实现了河段渠段状态量校正和参数校正后，校正量需要通过河道节点，以某种途径扩散到整个输水渠系，即实现整个输水渠系的校正。实时校正技术首推卡尔曼滤波技术，将若干个断面的校正分量合理全面地扩散到整个河段是它的长处，且是最小方差估计;可以用它来预报和校正糙率，以此形成河道水流和糙率的联解方法。

2.1.4.1　卡尔曼滤波法

2.1.4.1.1　概述

信号处理的实际问题，常常是要解决在噪声中提取信号的问题，因此，需要寻找一种所谓有最佳线性过滤特性的滤波器。这种滤波器当信号与噪声同时输入时，在输出端能将信号尽可能精确地重现出来，而噪声却受到最大抑制。卡尔曼滤波就是用来解决这样一类从噪声中提取信号问题的一种过滤(或滤波)方法。卡尔曼滤波理论采用时域法以状态方程为其数学工具，以多变量控制、最优控制与估计以及自适应控制为主要内容。其适用于线性随机理论，提出了平稳与非平稳、线性与非线性、集中与分布和多输入、输出系统相当广泛领域内的一种现代动态系统分析技术。

卡尔曼滤波不需要过去全部的观察数据，它只是根据前一个估计值和最近一个观察数据来估计信号的当前值。它是用状态方程和递推方法进行估计的，而它的解是以估计值(常常是状态变量的估计值)的形式给出的。卡尔曼滤波的信号模型是从状态方程和量测方程得到的。卡尔曼首先提出用一个状态方程和一个量测方程来完整地描述线性动态过程。卡尔曼滤波以系统状态空间模型为分析对象，用一个具有随机初始状态的向量来描述状态随时间的变化规律，称为状态方程。可以假定，这个方程受到某些随机干扰的影响以及描述模型不准确的干扰等，这些干扰称为动态噪声。量测变量与状态变量之间，可以假定有某些函数依赖关系，并可用方程描述，即量测方程。同时，还存在随机量测误差，这种误差称为量测噪声。对于离散时间系统，如果用{Xk}、{Zk}分别表示n维状态变量序列和m维量测变量序列，假设该动态系统为线性的，则常常将动态方程和量测方程一般化为如下形式。

状态方程:

量测方程:

式中:Φk为在k-1时刻的状态条件下k时刻的线性状态转移矩阵;Hk为k时刻的线性量测矩阵;X、Y分别为状态变量和量测变量;Wk为动态噪声;Vk为测量噪声。

Wk和Vk是均值为零的白噪声序列，且两者互不相关。

卡尔曼滤波运用现代随机估计理论给出了系统状态的无偏最小均方误差的递推估计值。如果估计的是系统状态当前时刻的值，称为“滤波”，如果是系统状态将来的值，称为“预报”。估计理论已经证明:无偏最小方差估计是线性估计中性能最佳的估计。卡尔曼滤波即是一种最佳线性估计方法。而其递推公式既可以得到滤波估计值，又可以得到误差的方差阵，即可以完成自身的误差分析。

卡尔曼滤波技术中最常用的为“正规卡尔曼滤波技术”，它对模型、观测噪声的统计特性作了一些假定。为了解决实际使用中存在的各种问题，卡尔曼滤波还包含了多种滤波处理技术，如描述非线性系统的扩展卡尔曼滤波方法，为解决有色噪声而产生的状态变量扩维法和量测求差法，为克服滤波发散而采用的渐消记忆滤波、限定记忆滤波、自适应滤波方法等。

2.1.4.1.2　卡尔曼滤波技术的优点及应用难点

卡尔曼滤波是线性无偏最小方差估计，主要有以下几个优点。

(1)可以应用于任何线性随机系统，校正能力全面而合理。

(2)可以综合处理模型误差与水文测验误差。

(3)只要满足卡尔曼滤波的条件，就可以得到最优估计。但做到这点必须与水动力学模型有机地结合起来，构建合适的状态空间方程。

(4)可以充分利用流域内的实测信息来校正状态量。对于那些黑箱子模型，它们仅仅使用出口断面流量来进行实时校正，而对于流域内部许多宝贵的实时水文信息却弃之不用。仅靠一个出口流量，信息量太少，又丧失了宝贵的预见期，这是导致许多模型实时校正效果不佳的主要原因。

卡尔曼滤波推理严密，这是其优点，但也意味着其使用条件较多，在实际应用上存在以下难点。

(1)模型噪声和量测噪声不为白噪声。例如，模型噪声往往难以满足许多滤波算法要求的白噪声这一约束性条件，模型参数率定的不准确将导致模型噪声序列的相关性，即Cov(ωk，ωj)(k≠j)不等于零，甚至E{ωk}不等于0。解决办法主要有:增加测报系统内水文站数目;采用时变线性系统;采用扩大状态维数等办法，把有色噪声看作是扩大维数后的状态向量的一部分，从而把有色噪声当白噪声处理;提高水文测验精度等。从理论上看，似乎考虑有色噪声是必需的，但在有色噪声情形下的线性滤波中，仍然需要引进新的参数，通过数值实验来确定其大小，其本质也是一个率定的过程，对于水文预报来说，不如仔细率定模型参数来得比较实际。目前看来，在解决这个问题方面，增加水文站数目效果是最好的。

(2)模型噪声和量测噪声难以统计，且不为常数。很多自适应滤波算法在运算多步后，都会进入稳态，需要进一步研究解决在线估计模型噪声和量测噪声的方法。

(3)滤波的发散问题。由于数学模型不能正确反映实际物理过程，对模型噪声和量测噪声的统计特性了解不够，使得两者方差取值不合适，加上计算机有限字长使得计算过程中的舍入误差积累过大，导致计算出来的协方差矩阵不正定，甚至不对称，最后使滤波出现反常现象，即尽管协方差矩阵的计算值看上去不断变小，但状态估计值的实际误差却不断增加，以致得到完全不符合实际的结果。出现发散的原因千差万别，误差补偿技术名目繁多。在水文系统中以衰减记忆法和增益直接加权法用得较多。它们都体现了误差补偿的基本思想:当系统模型不够好的时候，尽量削减模型的影响，强化对实测资料的依赖。但衰减记忆因子和人工干预增益调节因子的选取都要通过实验来解决。

2.1.4.1.3　扩展卡尔曼滤波基本原理

众所周知，河网水力学方程为非线性方程，而卡尔曼滤波法是在假设系统是线性的，噪声是白色的、高斯型的条件下的一种递推资料处理方法。因此，必须将非线性问题线性化以方便线性滤波。扩展卡尔曼滤波法是一种可以描述非线性系统的方法，为更好地处理水动力学仿真、预报提供了途径。

一般地，非线性系统的状态方程和量测方程可以表示为

假设在测量时刻k以前，已经得到滤波值。现将状态方程中的Φk围绕进行泰勒级数展开取其线性项，并用代替Γk(Xk-1，k-1)，可得状态方程的近似表达式:

式中:称为切线性状态算子。

同理，把量测方程中的Hk围绕状态预测值进行泰勒级数展开获取线性项，可得量测方程的近似表达式:

式中:称为切线性量测算子。

于是，非线性系统可线性化成如下的近似模型:

根据卡尔曼滤波递推算法的推导过程，最终可得非线性系统的扩展卡尔曼滤波递推算法:

由于线性化以后略去了高阶项，为了提高计算精度，可采用迭代扩展卡尔曼滤波，即每进行一次滤波后，将最新的滤波值视为预测值，并保持预报误差方差阵和增益矩阵不变，重复滤波过程对状态预测值进行多次校正。

2.1.4.2　水动力学模型实时校正

2.1.4.2.1　水动力学模型实时校正方法

水动力学模型的实时校正在提高预报精度和校核实测水情数据方面有着重要的实用价值。水动力学模型的校正包含状态量校正和参数在线估计。状态量校正需要根据河道若干个断面出现的水位、流量“新息”全面合理地对各断面的计算水位、流量做校正;参数估计需要实现渠道糙率、倒虹吸损失系数的自动调整。在实现了河段状态量校正和参数校正后，校正量需要通过河道节点，以某种途径扩散到整个河网，即实现整个河网系统的校正。此途径要求尽可能地降低算法所需要的内存空间。

实时校正技术首推卡尔曼滤波技术，将若干个断面的校正分量合理全面地扩散到整个河段是它的长处，且是最小方差估计。对水动力学模型进行实时校正，必须注意校正方法可能引起的水动力学模型不稳定。将卡尔曼滤波技术应用于水动力学模型的校正，由于滤波模型对实际系统概括的近似性破坏了水量平衡关系，可能导致水动力学模型自身的不稳定。所以，卡尔曼滤波和水动力学模型联合应用于实时校正需要攻克的难点问题是如何在选取稳定性好的差分格式的基础上，构建适当的状态空间方程，使得推导滤波公式的数学模型与实际物理系统的状态尽可能相符，防止可能引起的发散。此外，避免滤波算法的发散，正确估计滤波参数也是极为重要的。

2.1.4.2.2　河网水情数据同化的卡尔曼滤波算法(把圣维南方程组嵌入即可)

河网水力学方程一般可表示为

式中:X为状态变量(包括流量、水位等);t为时刻。

将河网水力学方程沿空间和时间两个方向离散，并考虑白噪声的影响，可得非线性的差分方程:

另外，河网量测方程一般为线性方程:

合并式(2.1-79)和式(2.1-80)，可得河网非线性动态系统的状态方程和量测方程:

将式(2.1-81)所示的状态方程线性化处理后有

最终，可得河网水情数据同化的扩展卡尔曼滤波递推算法如下。

(1)预测计算:

(2)滤波技术:

(3)增益矩阵:

(4)预报误差方差阵:

(5)滤波误差方差阵:

式中:为切线性状态算子;一般取单位矩阵I;状态转移矩阵Φk通过河网水力学方程的差分形式确定;观测矩阵Hk可根据测站监测的对象确定。

2.1.4.2.3　噪声协方差矩阵的估计

应用卡尔曼滤波的难点之一是如何估计模型噪声协方差矩阵Qk和量测噪声协方差矩阵Rk。

正规卡尔曼滤波公式有以下几个假定:

这里要求Qk和Rk是已知的。由于真值是未知的，Qk和Rk不能准确统计，而增益矩阵Kk的估计完全依赖于Qk和Rk，因此需要对Qk和Rk进行近似的估计。

梯级泵站输水系统中由于水位与流量数量级差异较大，一般而言，水位测量相对准确，一般采用水位作为状态变量。此时，量测噪声协方差矩阵和模型噪声协方差矩阵可以用以下方法计算。

量测噪声和测量仪器的误差有密切关系，通过试验分析，量测噪声协方差矩阵的元素值r的选择范围在0.0001～0.0025之间，则量测噪声协方差矩阵为

对于模型噪声协方差矩阵，其反映的是模型计算值与真值的偏差程度。虽然模型误差是一个非平稳序列，但是考虑到水文系统是一个时间渐变系统及模型误差无法避免的相关性，相邻时刻间的误差统计特性具有一定程度的相似性，故可以把实测值看作真值，在模拟计算中，动态选择一个较小的样本，计算模型噪声。

则模型噪声协方差矩阵为

2.1.4.3　实时校正技术在梯级泵站输水系统中的应用

对于工程运行管理人员来讲，更多的是关心如何对梯级泵站输水系统的水位、流量等状态量进行修正。然而在实际运行调度过程中发现，河道糙率、倒虹吸损失系数、泵站特性曲线等均会对系统的仿真模拟和预测产生较大影响。

(1)糙率校正。输水系统的河道糙率具有季节性变化的特点。春夏季节，输水河道内部水草滋生，会导致河道糙率增大;秋季，河道两岸树木落叶进入输水河道，也会导致河道计算糙率的变化;冬季，由于冰的加入，又会形成新的糙率。因此，梯级泵站输水系统的河道糙率应当分季节率定。同时，由于梯级泵站输水系统由泵站工程区分成若干个渠段，各渠段由于内部建筑物布置情况不同，导致一维模拟过程中综合糙率的不同，因此输水渠道的糙率还应该分渠段率定。

渠段糙率校正有两种途径:第一种可以采用实时校正之后的河网水情数据(水位、流量)作为目标值，通过优化算法或者人工试错法对糙率进行调整优化;第二种可以将糙率作为状态变量，对糙率进行实时校正，考虑糙率本身无法测量得到，故这里需要通过水位与糙率之间的函数关系，将量测变量由水位的函数表示，此时系统的状态方程和量测方程可以表示成如下形式:

式中:N由各相邻渠段之间的糙率n的平方组成;f(Zk)为水位的函数，水位与糙率最为紧密相关。

(2)倒虹吸损失系数率定。倒虹吸是梯级泵站输水系统中的常用渠系建筑物，由其引起的水力损失较大，水动力模型的仿真模拟结果对倒虹吸水力损失系数较为敏感。而实际计算过程中，倒虹吸的水力损失系数一般由查表得到，选择系数并不能反映其真实情况，故需要对其进行率定。

倒虹吸的水力损失系数较多，采用敏感性分析方法确定其中较为敏感的局部水力损失系数进行率定。率定所采用的参考值为数据同化后的水情数据，即对倒虹吸所属渠段的水位、流量过程进行校正，用校正结果率定倒虹吸水力损失系数。

(3)泵站特性曲线校正。梯级泵站输水系统中，泵站内部使用的抽水装置特性曲线一般是由模型试验数据通过相似定律转换到原型上来的，需要用实测数据对曲线进行校正。但是在实际工程中，由于测量仪器设备的误差、读数误差等导致实测数据与其真值有一定偏差。因此，可采用2.1.4.2节中提到的水情数据同化的方法，将实测的状态量进行修正，进而采用修正的状态量校核泵站特性曲线。