钱塘江河口围垦回淤过程预测探讨
摘要：在河口进行围垦，由于边界条件的改变和纳潮量的变化，一般都会引起河床的冲淤变化，对淤积过程的预测是一个非常重要且困难的问题。本文根据定床潮流的结果，利用河床变形方程求得围垦后的初始淤积速率，利用最小能耗原理求得极限淤积平衡状态，而由灰色理论求得整个淤积过程，经动床物理模型的印证和实测资料的检验以及工程实施后的过程的验证，均表明本文提出的预测方法是可行的，具有一定的可靠性和预测精度。
关键词：河床变形 回淤过程预测 灰色模型
1 前言
钱塘江河口潮强流急，宽浅多变，根据若干年来河口治理实践，对钱塘江河口主要采用趁潮淤围涂、固定江道、以围代坝等治江与围涂相结合的措施来进行河口整治。围垦之后，由于边界条件的改变，导致了水沙运行环境的改变，涨落潮流特性随之变化，净输沙格局也相应变化。一般地，围垦后，槽蓄量减小，涨、落潮动力减小，潮流所挟带的沙量也随之减小，与之相适应，河床过水断面将减小，产生回淤现象，回淤的幅度和时间过程对围垦工程而言是一个非常重要且困难的问题。由于泥沙问题的复杂性和不成熟性，采用泥沙数学模型进行长历时的计算受到诸多影响因素的制约，计算精度也就受到影响；采用物理模型代价又很大。因此，如何估算泥沙回淤过程是一个很有价值和值得研究的问题。本文从半理论半经验的途径出发，在定床潮流计算的基础上，利用河床变形方程、最小能耗原理、灰色模型对此问题进行了探讨，与物理模型结论和实测资料分析结论相一致，并且与工程实施后的回淤过程相吻合。
2 回淤机理概述
河流系统是一个开放的非线性系统，它与周围环境存在着物质和能量的交换，系统内部的自组织、自我调整作用是其演变发展的内因，外部环境的改变则是其变化的外因，外因要通过内因而起作用。围垦工程的实施，改变了来水来沙条件和水沙运行环境，必将引起河床的自我调整。决定河床形态变化趋势和程度的自我调节机制是受河流体系能量分配和耗散原理的制约，能量分配的调整使体系内能量分配具有最大或然率。这要求体系在调整过程中力求使熵达到最大，因而使各个可能出现状态的或然率相等。根据最小能耗原理，当流域的来水来沙条件发生变化时，河流将主要通过调整能量耗散率VJ来影响这种变化，而水流能量耗散率的沿程调整将由河宽B和水深H的沿程调整均等地承担，从而有
VJ=V·V2·n2/H4/3=q3n2/H13/3=min | (1) |
由式(1)可知，为适应边界条件的变化，河流系统中反应最为敏感、最易调整的因子是水深H，其次是单宽流量q，再次是糙率n，而变化相对缓慢、滞后的是min，它与施加于系统的约束有关。在围垦后，边界条件发生改变，在河流系统的自我调整作用下，各水力要素均会相应地发生变化以适应变化了的环境条件，由(1)式易得
H2/H1=m(q2/q1)9/13 | (2) |
m=(n2/n1)6/13(min1/min2)3/13 | (3) |
式中”1”代表工程前，”2”代表工程后，H为水深，q为单宽流量。
该形式已在钱塘江河口和其它航道整治工程中得到了广泛的应用，本文也将据此来推求河床调整平衡后的可能冲淤厚度。钱塘江河口以悬移质造床为主，在围垦前后，其糙率变化不大，即n1≈n2;并且局部的围垦尚不足以引起河口系统的能耗最小值的大幅度相对变化，因此也可假设(min1/min2)3/13≈1，为此近似地取m≈1。根据(2)式，围垦后纳潮量减小，即q2减小，则H2应减少，河床将发生淤积，其回淤关系将满足式(2)，利用定床水流计算所得的水力要素和(2)式可求得围垦后可能的平衡淤积厚度。
另一方面，对于潮汐河口，由于潮波变形，涨、落潮不对称，对输水输沙的作用是不一样的。在某时刻的单宽输沙率可写成
| (4) |
式中
,
,
分别为潮周期的垂线平均流速、含沙量和水深，V′，S′，H′分别为某一时刻平均流速、含沙量和水深与垂线平均流速、含沙量和水深的偏差。以潮周期为计算单元时，有
∑V′=∑S′=∑H′=0 | (5) |
故潮周期平均单宽输沙率 为 |
|
(T1) (T2) (T3) (T4) (T5) | (6) |
式中T1项中，，恒为正，正负两可，因此，T1的正负(方向)仅取决于涨、落潮流的强度对比。T1表征了优势流对净输沙的贡献。T2项反映了stokes漂流效应对净输沙的贡献，在河口地区该项为负值。T3的符号依V′S′的正负而定，后者既取决于流速过程线的含沙量过程线和相位关系，还与涨、落潮掀沙效应的差异有关。涨潮含沙量大于落潮含沙量时，T3为负，反之为正。故T3反映了涨、落潮流挟沙强度对净输沙的贡献。T4为含沙量与水深变化的相关项，因研究区段处在动力平衡带，AKV-的绝对值通常较小，故T4的绝对值也较小，T5项依赖于V，S，H的相关性，其绝对值与前三项相比小得多。
由方程(4)可知，潮流对净输沙的贡献取决于涨、落潮流特性的相对变化。围垦之后，由于流速、潮量的减小，涨、落潮动力减小，潮流所挟带的沙量也随之减小，涨、落潮优势流和挟沙强度的相对变化决定了围垦之后河口的净输沙格局的变化和演变发展趋势。工程前后，河道的冲淤主要取决于优势流与挟沙强度的相对变化，为此可用输水输沙的相对变化来判断可能发生的冲淤情况。比如，对涨潮冲刷槽而言，当涨潮输沙量大于落潮输沙量时，将发生冲刷，可用优势沙表示
G=(VHSt)f/(VHSt)e∝(Vf3tf)/(Ve3te)S∝V2/H | (7) |
式中下标f表示涨潮，下标e表示落潮，S为含沙量。
同样，定义优势流为涨潮潮量与落潮潮量之比，即优势流为涨、落潮相对输水量之比
F=(Qf/Qe) | (8) |
显然，当?G2>G1,F2>F1时，河床将发生冲刷，当G2<G1,F2<F1时，河床将发生淤积。由优势流和优势沙的概念，可以利用定床计算的水力要素来判断围垦后河床发生冲、淤的部位及可能发生的演变趋势。根据上述概念所得到的冲、淤平面分布与工程实施后出现的分布情况相一致。
3 回淤预测公式一般形式
根据河床变形方程[1]可建立如下冲淤模式
| (9) |
式中 k为沉降机率(在钱塘江河口可取1～2)；ω为泥沙沉降速度(4～5×10-4m/s);s为含沙量；s*为挟沙力；γ′为泥沙干容重。
在围垦之前(天然情况下)，可认为河床处于相对冲淤平衡状态，即
，有s1=s1*，在潮汐河口，有
s*=α(v2/gh) | (10) |
围垦后初期，假定来沙不变，即s2=s1，由于围垦引起水流条件的改变，必然引起挟沙能力的相应变化，造成河床的冲淤变化，其冲淤强度可由(9)式的差分而求得，即
p=kw(s2-s2*)/γ′ | (11) |
代入挟沙关系(10)可得 | |
p=kws1/γ′[(1-(v2/v1)2(h1/h2)] | (12) |
上式中，当p→0时，可认为得到冲淤极限平衡状态，即h2→H2，从而可求得平衡水深H2 | |
H2=(q2/q1)aH1 | (13) |
式中α=2/3,q1=v1h1,q2=v2h2，水力要素均由定床潮流求得。上式在钱塘江河口中已达到了广泛的应用，它与前文根据最小能耗原理所得的平衡关系式(2)相接近，但该式隐含了工程前后挟沙力不变的假设，从而也说明了最小能耗原理更具有一般性，可用于一般的河床演变分析。当围垦工程实施后，由于河床变形相对滞后，可假设围垦后初期含沙量不变，从而可利用定床计算所得的水力要素，由式(12)可求得围垦后初期的淤积速率p。
对于潮汐河口，正如前述，由于涨、落潮对泥沙的输移作用是不一样的，为了利用上述由冲积河流所得的公式(12)，可以将公式中对应的水力要素取全潮平均进行计算。当取半潮平均时，可以分别得到半潮淤积速率p1(涨潮)和p2(落潮)，则淤积速率取二者相比较大者，p=max(p1,p2)，或根据当地的冲淤特性来选取，对涨潮占优者取p1，对落潮占优取p2，对涨、落持平者，取p=(p1+p2)/2作为初期淤积速率。
4 回淤过程的灰色预测
上文已求得围垦后初期的可能淤积速率和平衡淤积厚度。由于河床的自我调整作用，随着淤积的，水流要素也会受到影响而发生变化，从而进一步地影响河床的回淤速率。从原则上讲，可由(12)式进行迭代求解整个回淤过程，但由于在迭代过程中，水沙要素的变化不易确定，为进一步求解带来了一定的困难。在此，采用灰色预测方法来预估回淤过程。由于泥沙系统是一个本征性的灰色系统，可以采用GM(1，1)[2]模型来进行预测，利用上文已求得的初始淤积速率和平衡淤积厚度，可以求得回淤过程和淤积速率的变化过程。
对GM(1.1)模型有
| (14) |
其中,
为待定系数，z(1)为累积淤积厚度。
对灰色模型(14)而言，有解析解
| (15) |
对上式求导，可得淤积速率满足 | |
| (16) |
为方便起见，以初始地形为基准，即取z0=0。显然，当淤积得到平衡时，有 | |
；在淤积初期有
，根据上述条件和前文已求得的平衡淤积厚度和初始淤积速率，可以求得待定系数5 应用实例*
根据上述思想，在钱塘江河口尖山河湾南股槽整治研究中，建立了平面二维潮流模型进行数值模拟，利用该定床模型所求得的水流参数和实测含沙量资料，对围垦后的河床变形和演变趋势进行了分析预测。
*宋立松等. 钱塘江尖山河湾南股槽整治数学模型——整治研究之三.浙江省河口海岸研究所(内部报告)，1996.12。
5.1 定床潮流模型
针对钱塘江尖山河湾特点，建立了平面二维潮流模型，其控制方程为连续方程：
| (17) |
运动方程： |
|
| (18) |
| (19) |
式中z为潮位，即水面到达一基准面的距离；U，V分别为x,y方向上的垂线平均流速分量；h为水深；g为重力加速度；f为柯氏力参数(f=2ωesinφ,φ为纬度，ωe为地球自转速度)；Cz为谢才系数；εx,εy分别为水流在x,y方向的涡动扩散系数；τx,τy分别为x,y方向的风应力分量。 上述方程组的初始条件
|
|
边界条件 | |
水边界: | z(x,y,t)=z*(x,y,t)“*”表示已知值， |
陆边界: |
|
法线方向流速为零； | |
在南股槽整治研究中，选用了无网格的差分方法对控制方程进行离散，即对式(18)，(19)应用不规则网格差分离散，对离散方程进行了数值求解，求得围垦前后的有关水力要素(U1,V1,h1,U2,V2,h2),从而可得到：
,q1=v1h1, 
,q2=v2h2,有了以上参数，即可应用本文方法进行围垦后回淤过程的预测。
5.2 回淤预测 从优势流(F)与优势沙(G)的相对比较来看，二者所反映的趋势是基本一致的，淤积基本上以东进闸为中心(图1)，处于淤积环境中，个别地方反而略有冲刷，下游较远处(临海浦以下)基本上可维持平衡。利用定床水流的水力结果(v2,v1,h2,h1,q1=v1h1,q2=v2h2)和实测资料分析所得的泥沙参数(s1,ω)，由式(2)和式(9)求得平衡淤积厚度和初始淤积速率如表1。 |
| ||||||||||||||||||||||||||||||||||||||||
根据表1，由(15)和(16)求得东近闸附近的回淤过程满足如下方程
对近期方案(围垦6万亩)有
| (20) |
| (21) |
对中期方案有(围垦8.6万亩：即近期6万亩及中期2.6万亩) | |
| (22) |
| (23) |
由上述方程可知，近期方案实施后的淤积平衡厚度约为4.44m，初始淤积速率为4.1cm/d，得到平衡所需时间约为300～500天左右；中期方案实施后的淤积平衡厚度为8.39m，初始淤积速率约为5.09cm/d，得到平衡的时间约为500～700天左右(如图2所示)，可见其淤积速度是较快的。
|
图2 近期、中期方案实施后东进闸附近的回淤过程 |
Predicted aggradation process near Dongjing barrage |
6 结果评价与讨论
根据70年代的实测资料分析*，尖山河湾的时段冲淤量(划分冬半年、夏半年)约0.3～0.6×108m3。从尖山河湾多次汊道发生淤积的实测资料表明(1969年7.5洪水北沙串通，南股槽淤积和1973年，1985年尖山一带主槽由南靠北，南岸汊道的淤积)，当其中的一股汊道出现淤积环境时，其淤积速度是很快的，一般约0.5～1年。当床面高程为-1～0m时，淤积率2～3cm/d。南股槽内水深较大，因而估计围垦后其淤积速率前期会略大于上述实测淤积率，因此，按上述模式所求得的淤积强度略大是可以接受的，并且由上述方法所得到的冲淤部位，淤积平衡厚度、回淤时间等都与动床物理模型**的结论相一致，也与实测资料分析结果相符合，1996～1997年近期方案实施后出现的淤积面貌及其过程均与本文预测吻合良好(图2)，从而也说明上述预测的河床变形有一定的价值和可信度。而且利用所得到的灰色模型，还可以从数学的角度[3，4]对河床演变进行敏感性、稳定性等各种动态行为的分析，为进一步揭示河口泥沙运动提供了一种途径。
*余祈文等，钱塘江尖山河湾南股槽的形成及河床演变特性研究，浙江省河口海岸研究所(内部报告)，1992.12。
**熊绍隆等，钱塘江尖山河湾南股槽整治动床模型——整治研究之二，浙江省河口海岸研究所(内部报告)，1996.12。
当然，由于河口泥沙问题的复杂性和特殊性，要精确地预测河床长历时的变形是很困难的，目前尚无有效手段。采用泥沙数学模型进行长历时的计算受到诸多影响因素的制约，计算可靠性及精度也就受到影响；采用物理模型费时又费钱，代价很大；定床潮流模型只能反映水力变化，不能反映河床变形；最小能耗原理或河相关系只能反映河床的平衡状态，难以揭示其发展变化过程；灰色模型可以描述过程，但其动态行为对初值、终值有一定的依赖性，也不能完全确定河床演变过程。本文从半理论半经验的途径出发，根据定床潮流计算的结果，利用河床变形方程求得围垦后的初始淤积速率，利用最小能耗原理求得极限淤积平衡状态，而由灰色理论求得整个淤积过程，充分利用不同方法的长处对此问题进行了初步探讨，鉴于河床长历时演变的复杂性，本文方法尚需进一步的完善和实践的检验。
1 谢鉴衡(主编).河流模拟。电力出版社，1988.
2 邓聚龙.灰色预测与决策.华中工学院出版社，1986.
3 邓聚龙.灰色控制系统.华中工学院出版社，1985.
4 邓聚龙.多维灰色规划.华中理工大学出版社，1989.