三峡工程运用后，荆江河段受上游来水的持续冲刷，床沙逐年粗化，只是粗化程度与粗化速率各异^[1]。同时，长江中游水沙条件多变，使得该河段年内床沙交换过程更为复杂，床面泥沙颗粒间的相对位置分布规律更具随机性，其内在机理有待进一步研究。河床粗化或床沙交换属于泥沙运动统计理论的主要内容，是研究非平衡输沙及河床演变规律的理论基础，也是河流动力学领域的重要课题之一^[2]。

目前，关于床沙交换及粗化过程的研究成果可大致分为两类：一是河床粗化，此部分内容基于床沙与推移质泥沙的交换，只考虑床沙起动或推移质落淤。Robert^[3]在Little等^[4]模型的基础上，定义粗化层厚度为E=Ad₉₅（A为常数），得到了计算粗化层级配的概率—体积方法，该计算方法可应用于水槽试验及天然河流中。许全喜^[5]、孙志林^[6]等提出考虑非均匀隐暴效应的泥沙起动概率，并建立基于起动概率的河床粗化层多步计算模式，该模式可同时计算河床冲刷深度；但是，许全喜等^[5]认为床沙在冲刷粗化过程中，活动层厚度E应随床沙组成的变化而变化，对于卵石夹沙河床，活动层厚度为E=d₉₅，E随着床沙粗化的进行不断增大；孙志林等^[6]提出的计算粗化层级配的多步模式未考虑床沙起悬及悬沙淤积，并假定床沙粗化过程中活动层厚度E保持不变，只改变床面活动层级配和高程。基于上述观点，作者认为在实际河流中，随着床沙粗化过程的进行，活动层厚度E应逐渐减小，并在粗化完成时趋于定值；同时，在计算床沙粗化级配变化时需要考虑冲淤过程的影响。二是床沙与推移质、悬移质泥沙之间的交换过程，此部分内容基于两态或三态转移概率，考虑了床沙、推移质或悬移质泥沙运动状态的相互转换，但关于三态转移概率应用方面的研究还很少。由于单颗粒泥沙起动概率是两态或三态转移概率的前提，因此一些学者对单颗粒泥沙进行受力分析，从不同角度得到均匀沙滚动概率、悬移概率^[7-9]；有些学者进一步考虑了水流脉动、床面泥沙颗粒间隐暴效应对泥沙起动的影响，建立非均匀沙滚动概率、跃移概率^[10-11]，并用水槽试验资料进行验证。之后，Sun等^[12]基于Komogorov微分方程及Markov定理，得到非均匀沙两态转移概率矩阵；首次提出状态占用概率，推导出非均匀沙分级输沙率公式。Yang等^[13]基于生物生灭过程及马尔可夫随机过程，建立床沙冲刷粗化计算模型，进而得到预测推移质输沙率与冲刷历时之间的分段函数关系式，与实测资料符合较好。另外，王涛^[14]、聂锐华^[15]等通过水槽试验，研究了卵石夹沙河床清水冲刷过程中粗化层破坏的临界条件及粗化层破坏过程中的最大推移质输沙率。值得注意的是，上述模型只考虑了床沙与推移质泥沙的相互交换，忽略了悬沙的存在，所以只能应用于清水冲刷试验。韩其为^[16]提出床面泥沙改变状态的6种基本概率，从理论角度比较系统的总结了床面泥沙4种运动状态（静止、滚动、跳跃、悬浮）之间的转移概率及交换强度，并由12种交换强度导出4种运动状态之间的9种输沙能力公式，与现有的推移质及悬移质输沙公式不同。

通过以上分析可知，应用概率模型计算天然河流床沙交换过程的研究很少。本文采用力学分析与概率统计理论相结合的方法，基于Markov三态转移概率矩阵，引入非均匀沙隐暴系数，建立卵石夹沙河床上床沙级配、推移质级配及悬沙级配的概率计算模型，该模型能够同时考虑前期水沙条件、床沙起悬及冲淤过程对床沙交换过程的影响。最后，采用本文的概率计算模型，计算了枝城断面卵石夹沙河床年际床沙粗化及年内床沙交换过程，并与2003—2009年部分实测资料进行了比较。

1 研究河段概况 1.1 宜昌—江口河段

如图1所示，宜昌—江口河段包括两部分：宜枝河段、枝江河段（枝城—江口河段）；宜枝河段位于三峡大坝下游约43.1 km处，上起宜昌，下迄枝城，全程长约59.0 km；枝江河段位于三峡大坝下游约102.4 km处，上起枝城，下迄江口。若根据地理位置划分，宜枝河段以虎牙滩为界又可分为宜昌河段、宜都河段，河长分别约为19.4、39.6 km；若根据河床组成划分，宜枝河段以宜都为界分为河段1（宜昌—宜都）、河段2（宜都—枝城），河长分别约为42.0、17.3 km。由于本研究适用条件为卵石夹沙河床，故以第2种划分结果作为依据。由图1可以看出，河段1沿程包括宜昌、胭脂坝、虎牙滩、红花套等地，河段2沿程包括宜都、白洋、枝城等地。其中，宜昌—宜都河段为卵石夹沙河床，河床上有卵砾石洲滩分布；宜都—枝城河段同样为卵石夹沙河床，但河床组成多为中细沙，相对河段1来说床沙粒径范围比较均匀，卵石层深埋在床面以下^[17-18]。另外，枝江河段（图1河段3所示）表层床沙组成多为卵砾石，为典型的卵石夹沙河床，也是本文重点研究河段。由于缺少宜昌等站水文资料，故本研究所用验证资料来自枝城站（图1河段2及河段3交界处，属于卵石夹沙河床），符合本文概率计算模型的适用条件。

1.2 近期河床冲刷粗化情况

在一定的水沙条件下，床沙、推沙（推移质泥沙）、悬沙（悬移质泥沙）相互交换，从而使河床发生冲淤变化。因此，为研究三峡工程下游宜枝河段河床冲淤过程，首先要了解长江中游水沙条件变化及床沙交换过程。为了进一步说明枝城断面的床沙交换及粗化过程，图2点绘了枝城站1994—2017年水量及沙量逐年变化过程^[19]。由图2可知，三峡工程运用后（2003—2017年），年均水量减少3.7%，但沙量减少近88.4%。

图3为三峡工程运用后，宜昌—江口河段荆3、荆4断面2002—2017年河床高程变化情况（枝城河段右岸有护岸工程）。如图2、3所示，三峡大坝建成后，大量泥沙在库前淤积，下泄水流含沙量处于严重不饱和状态，下泄清水持续冲刷河床，造成长江中游河段床沙逐年粗化，其中枝城河段粗化情况最为明显，且冲刷部位集中在枯水河槽^{[18, 20-21]}。

通过搜集整理，三峡工程下游枝城站2003—2017年实测年均床沙级配及中值粒径d₅₀资料如图4所示。

总体来说，受持续冲刷的影响，枝城段河床泥沙粒径逐年粗化，粒径范围变宽，并在2007年发生突变，发生突变的原因可能是在2006年12月到2007年1月时间段内，水流冲刷强度相对强烈，使得床沙迅速粗化。

由图4（a）可知，枝城站2003年初始床沙级配粒径范围为0.002～41.000 mm，中值粒径为d₅₀=0.23 mm；2017年河床粗化后床沙级配粒径范围为0.031～86.000 mm，中值粒径变为0.424 mm；由图4（b）可以看出，在2003—2017年间，该断面床沙中值粒径d₅₀总体呈上升趋势，并在2010年以后有一定波动；同时，2007—2012年粗化速率慢于2003—2007年粗化速率，但2012年以后（2012—2017年）粗化速率又明显增大。

2 Markov随机过程 2.1 三态转移概率

在一定的水流条件下，处于悬移状态的泥沙颗粒可与床面泥沙相互交换，由悬移状态3变为推移状态2、静止状态1，或者继续保持悬移状态3。同理，处于床面上的泥沙可能从静止状态变为推移或悬移状态；处于推移状态的推移质泥沙颗粒也可能变为静止状态或悬移状态。如图5、6所示，这一状态转移过程称为Markov随机过程。

在已有研究中，一般只考虑床沙与推移质或者床沙与悬移质之间的交换，即两态随机交换模型^{[7-8, 10-11, 22]}，且一般只适用于水槽试验、均匀沙及输沙平衡状态，在天然河流应用中存在局限性。因此，为了计算天然河流中的床沙交换过程，需要建立同时考虑床沙、推移质、悬移质转移概率的三态模型；若将悬移、跳跃视为同一运动状态（悬移），那么非均匀沙运动的三态转移概率矩阵 ${{\varepsilon}} _i^{(n)}$ 可以表示为：

${{\varepsilon}} _i^{(n)} = \left[\!\!\!\! {\begin{array}{*{20}{c}} {\varepsilon _{11i}^{(n)}}&{\varepsilon _{12i}^{(n)}}&{\varepsilon _{13i}^{(n)}}\\ {\varepsilon _{21i}^{(n)}}&{\varepsilon _{22i}^{(n)}}&{\varepsilon _{23i}^{(n)}}\\ {\varepsilon _{31i}^{(n)}}&{\varepsilon _{32i}^{(n)}}&{\varepsilon _{33i}^{(n)}} \end{array}}\!\!\!\! \right]$

(1)

式中：i为分组粒径； ${{\varepsilon}} _i^{(n)}$ 表示第n步状态转移矩阵； $\varepsilon _{xyi}^{(n)}$ 为状态转移概率（床沙交换），表示泥沙颗粒由x状态转移到y状态的概率（其中：x、y=1、2、3，“1”表示静止，“2”表示滚动或滑动，“3”表示悬移）。

Kuai等^[23]作出如下假设：1）由推移状态转移到悬移状态的概率 $\varepsilon _{23i}^{(n)}$ 与泥沙起悬概率 $P_{{\rm{S}}i}^{(n)}$ 近似相等，即 $\varepsilon _{23i}^{(n)} \approx P_{{\rm{S}}i}^{(n)}$ ；2）由悬移状态转移到静止状态的概率近似为0，即 $\varepsilon _{31i}^{(n)} \approx 0$ ；3）天然河流时间尺度 $\Delta t$ 假定为水深h与泥沙沉降速度w之比，即 $\Delta t = h/w$ 。从而得到，式（1）中各分组粒径i泥沙的三态转移概率矩阵 ${{\varepsilon}} _i^{(n)}$ ：

${{\varepsilon}} _i^{(n)} = \left[\!\!\!\! {\begin{array}{*{20}{c}} {1 - P_{\rm{T}i}^{{\rm{(}}n{\rm{)}}}}&{P_{\rm{T}i}^{{\rm{(}}n{\rm{)}}} - \varepsilon _{13i}^{(n)}}&{\varepsilon _{13i}^{(n)}}\\ {1 - \varepsilon _{12i}^{(n)}}&{\varepsilon _{12i}^{(n)} - P_{\rm{S}i}^{{\rm{(}}n{\rm{)}}}}&{P_{\rm{S}i}^{{\rm{(}}n{\rm{)}}}}\\ 0&{\beta \varepsilon _{21i}^{(n)}}&{1 - \beta \varepsilon _{21i}^{(n)}} \end{array}} \!\!\!\!\right]$

(2)

式中： $P_{{\rm{T}}i}^{(n)}$ 为总概率，且 $P_{{\rm{T}}i}^{(n)} = P_{{\rm{R}}i}^{(n)} + P_{{\rm{S}}i}^{(n)}$ ； $P_{{\rm{S}}i}^{(n)}$ 为起悬概率； $\;\beta $ 为概率修正系数，且 $\;\beta = {{\varepsilon _{32i}^{(n)}} / {\varepsilon _{21i}^{(n)}}}$ 。

在式（2）所述模型中，存在不足之处：1）假设 $\varepsilon _{23i}^{(n)} \approx P_{{\rm{S}}i}^{(n)}$ 、 $\varepsilon _{31i}^{(n)} \approx 0$ 不合理^[23]；2）天然河流时间尺度的选择，缺少天然河流实测资料验证；3）该模型中部分系数需要实测资料率定。

考虑到Kuai等^[23]提出的三态Markov非齐次离散模型存在以上不足，作者对三态Markov非齐次离散模型进行修正，从而得到基于隐暴效应的非均匀沙三态转移概率矩阵 ${{\varepsilon}} _{xy}^{(n)}$ 为：

${{\varepsilon}} _{xy}^{(n)} = \left[\!\!\!\! {\begin{array}{*{20}{c}} {1 - P_{{\rm{T}}i}^{(n)}}&{P_{{\rm{T}}i}^{(n)} - P_{{\rm{S}}i}^{(n)}}&{P_{{\rm{S}}i}^{(n)}}\\ {1 - P_{2i}^{(n)}}&{P_{{\rm{R}}i}^{(n)}P_{2i}^{(n)}}&{P_{2i}^{(n)}(1 - P_{{\rm{R}}i}^{(n)})}\\ {1 - P_{3i}^{(n)}}&{P_{3i}^{(n)}(1 - P_{{\rm{S}}i}^{(n)})}&{P_{{\rm{S}}i}^{(n)}P_{3i}^{(n)}} \end{array}} \!\!\!\!\right]$

(3)

式中， $P_{{\rm{R}}i}^{(n)}$ 为泥沙推移概率， $1 - P_{2i}^{(n)}$ 、 $1 - P_{3i}^{(n)}$ 分别为止滚概率、止悬概率。

在式（3）中，首先确定对角线转移概率 $\varepsilon _{11i}^{(n)}$ 、 $\varepsilon _{22i}^{(n)}$ 、 $\varepsilon _{33i}^{(n)}$ （保持原运动状态不变），第1列转移概率 $\varepsilon _{21i}^{(n)}$ 、 $\varepsilon _{31i}^{(n)}$ （由动到静）。然后根据各行概率之和为1，确定其他转移概率值。

2.2 起滚、起悬概率

对于非均匀沙来说，由于受到泥沙颗粒间隐暴效应的影响，可以认为泥沙起悬、起滚相互独立；根据韩其为等^[24]的研究，止动流速U₀一般比起动流速U小，即 ${U_0} = aU \approx 6.725{d_i}^{0.5}(a < 1)$ ；所以，止动概率 $1 - P_{2i}^{(n)}$ 可以表示为：

$1 - P_{2i}^{(n)} = \frac{1}{{\sqrt {2\text{π} } }}\int_{ - \frac{{{U_0}}}{{2{u_ * }}} - 2.7}^{\frac{{{U_0}}}{{2{u_ * }}} - 2.7} {{e^{ - \frac{{{t^2}}}{2}}}} {\rm{d}}t$

(4)

因此，泥沙的不起动概率（ $1 - P_{{\rm{R}}i}^{(n)}$ 、 $1 - P_{{\rm{S}}i}^{(n)}$ ）与止动概率 $\varepsilon _{21i}^{(n)}$ （止滚概率： $1 - P_{2i}^{(n)}$ ）、止悬概率 $\varepsilon _{31i}^{(n)}$ （ $1 - P_{3i}^{(n)}$ ）不等，即：

$\left\{\!\!\!\! {\begin{array}{*{20}{c}} {1 - P_{{\rm{R}}i}^{(n)} \ne 1 - P_{2i}^{(n)}} \\ {1 - P_{{\rm{S}}i}^{(n)} \ne 1 - P_{3i}^{(n)}} \end{array}} \right. \Rightarrow \left\{\!\!\!\! {\begin{array}{*{20}{c}} {P_{{\rm{R}}i}^{(n)} \ne P_{2i}^{(n)}} \\ {P_{{\rm{S}}i}^{(n)} \ne P_{3i}^{(n)}} \end{array}} \right.$

(5)

参考韩其为^[24]、胡春宏^[25]等的研究，止滚概率 $1 - P_{2i}^{(n)}$ 与止悬概率 $1 - P_{3i}^{(n)}$ 有如下关系：

$1 - P_{3i}^{(n)} = (1 - P_{{\rm{S}}i}^{(n)})(1 - P_{2i}^{(n)})$

(6)

将式（6）化简，则不止悬概率 ${P_{3i}}(l)$ 可以表示为：

$P_{3i}^{(n)} = P_{2i}^{(n)} + (1 - P_{2i}^{(n)})P_{{\rm{S}}i}^{(n)}$

(7)

根据Li等^[22]的研究，考虑隐暴效应的非均匀沙滚动、起悬概率分别为：

$\begin{aligned}[b] P_{{\rm{R}}i}^{(n)} &= \dfrac{1}{2}\left\{ {\frac{{0.21 - A}}{{\left| {0.21 - A} \right|}}\sqrt {1 - \exp \left[ { - {{\left( {\frac{{0.46}}{A} - 2.2} \right)}^2}} \right]}- } \right. \\ & \left. {\dfrac{{0.{\rm{135 - }}A}}{{\left| {0.{\rm{135 - }}A} \right|}}\sqrt {1 - \exp \left[ { - {{\left( {\frac{{0.295}}{A} - 2.2} \right)}^2}} \right]} } \right\} \end{aligned} $

(8)

$P_{{\rm{S}}i}^{(n)} = \frac{1}{2}\left\{ {1 - \frac{{0.21 - A}}{{\left| {0.21 - A} \right|}}\sqrt {1 - \exp \left[ { - {{\left( {\frac{{0.46}}{A} - 2.2} \right)}^2}} \right]} } \right\}$

(9)

式（8）～（9）中：A为综合系数， $A = \sqrt {{\xi _i}{\varTheta _i}{C_{\rm{L}}}} $ ； ${\xi _i}$ 为非均匀沙隐暴系数，且 ${\xi _i} = \sigma _{\rm{g}}^{0.25}{({d_i}/{d_{\rm{m}}})^{0.5}}$ ， ${\sigma _{\rm{g}}} = \sqrt {{{{d_{84.1}}} / {{d_{15.9}}}}} $ 为非均匀系数， ${d_{\rm{m}}}$ 为平均粒径； ${\varTheta _i}$ 为第i粒径组泥沙的水流强度参数； ${C_{\rm{L}}}$ 为上举力系数，一般取为0.4。

3 床沙交换概率计算模型 3.1 床沙级配概率计算模型

根据许全喜^[5]、孙志林^[6]等的相关研究，卵石夹沙河床的活动层沙量平衡方程为：

$EP_{{\rm{a}}i}^{(n - 1)} - \Delta {H^{(n)}}P_{{\rm{b}}i}^{(n)} + \Delta {H^{(n)}}{P_{0i}} = EP_{{\rm{a}}i}^{(n)}$

(10)

在式（10）中，未考虑床沙起悬及悬沙淤积对床沙级配变化的影响，床沙粗化过程中活动层厚度E保持不变，并认为活动层内每步被冲刷下移的沙量由等量的下层原始床沙来补充，不适用于天然河流床沙粗化过程的计算，所以对式（10）进行改进：

$\begin{aligned}[b] {E^{(n)}}P_{{\rm{a}}i}^{(n - 1)} &- \Delta H_1^{(n)}[P_{{\rm{b}}i}^{(n)} + P_{{\rm{s}}i}^{(n)}]+ \\ & \Delta H_2^{(n)}[P_{{\rm{b}}i}^{(n - 1)}\varepsilon _{21i}^{(n)} + P_{{\rm{s}}i}^{(n - 1)}\varepsilon _{31i}^{(n)}] = {E^{(n)}}P_{{\rm{a}}i}^{(n)} \\ \end{aligned} $

(11)

式（10）～（11）中：i为床沙粒级； ${E^{(n)}}$ 为第n步床沙交换活动层厚度，与床面泥沙的最大起动粒径 $d_{c\max }^{(n)}$ 有关； $P_{{\rm{a}}i}^{(n)}$ 、 $P_{{\rm{b}}i}^{(n)}$ 、 $P_{{\rm{s}}i}^{(n)}$ 分别为第n步床沙级配、推移质级配、悬移质级配，代表泥沙粒径小于d_i的沙重百分比； $\varepsilon _{13i}^{(n)}$ 为状态转移概率（床沙交换），表示泥沙颗粒由悬移状态转移到静止状态的概率（其中：“1”表示静止，“2”表示滚动或滑动，“3”表示悬移）； $\Delta {H^{(n)}}$ 为第n步床面总冲淤厚度； ${P_{0i}}$ 为初始床沙级配，且存在 $P_{{\rm{a}}i}^{(0)} = {P_{0i}}$ 。式（11）左边第2项 $\Delta H_1^{(n)}[P_{{\rm{b}}i}^{(n)} + $ $P_{{\rm{s}}i}^{(n)}]$ 为冲刷项，第3项 $\Delta H_2^{(n)}[P_{{\rm{b}}i}^{(n - 1)}\varepsilon _{21i}^{(n)} + P_{{\rm{s}}i}^{(n - 1)}\varepsilon _{31i}^{(n)}]$ 为淤积项。

3.2 推移质级配及悬移质级配

将床沙、推移质、悬沙视为整体（认为推移质落淤后剩余的推移质泥沙、床沙起滚部分泥沙、悬沙转移到推移质状态的泥沙共同组成第n步推移质级配），并假设泥沙总量为1，则第n步推移质级配 $P_{bi}^{(n)}$ 可以表示为：

$P_{{\rm{b}}i}^{(n)} = P_{{\rm{a}}i}^{(n - 1)}\varepsilon _{12i}^{(n)} + P_{{\rm{b}}i}^{(n - 1)}\varepsilon _{22i}^{(n)} + P_{{\rm{s}}i}^{(n{\rm{ - }}1)}\varepsilon _{32i}^{(n)}$

(12)

同理，可以得到第n步悬移质级配 $P_{{\rm{s}}i}^{(n)}$ ：

$P_{{\rm{s}}i}^{(n)} = P_{{\rm{a}}i}^{(n - 1)}\varepsilon _{13i}^{(n)} + P_{{\rm{b}}i}^{(n - 1)}\varepsilon _{23i}^{(n)} + P_{{\rm{s}}i}^{(n{\rm{ - }}1)}\varepsilon _{33i}^{(n)}$

(13)

式（12）～（13）中， $\varepsilon _{1yi}^{(n)}$ 为由静止状态转移到y状态的概率， $\varepsilon _{3yi}^{(n)}$ （其中：y=1、2、3）为由悬移状态转移到y状态的概率。

3.3 冲淤厚度

根据许全喜^[5]、孙志林^[6]等的研究，计算时段 $\Delta T$ 内冲刷厚度 $\Delta H_1^{(n)}$ 为：

$\Delta H_1^{(n)} = {T^{(n)}}E_1^{(n)}(\sum\limits_{i = 1}^{\max } {P_{{\rm{a}}i}^{(n - 1)}} \varepsilon _{13i}^{(n)} + \sum\limits_{i = 1}^{\max } {P_{{\rm{a}}i}^{(n - 1)}} \varepsilon _{12i}^{(n)})$

(14)

式中： $\varepsilon _{13i}^{(n)}$ 、 $\varepsilon _{12i}^{(n)}$ 分别为床沙由静止状态转移到悬移状态的概率、由静止状态转移到滚动（滑动）状态的概率； $E_1^{(n)}$ 为某一床沙交换过程中活动层的厚度，且 ${T^{(n)}}E_1^{(n)} = {E^{(n)}}$ ； ${T^{(n)}}$ 为相对计算步长（相对于一个交换过程），且 ${T^{(n)}} = \Delta T/\Delta t$ ， $\Delta T$ 为实际计算步长， $\Delta t$ 为完成一次床沙交换所需要的时间。

参考冲刷厚度的定义，计算时段 $\Delta T$ 内淤积厚度 $\Delta H_2^{(n)}$ 可以表示为：

$\Delta H_2^{(n)} = E_2^{(n)}(\sum\limits_{i = 1}^{\max } {P_{{\rm{s}}i}^{(n{\rm{ - }}1)}} \varepsilon _{31i}^{(n)} + \sum\limits_{i = 1}^{\max } {P_{{\rm{b}}i}^{(n{\rm{ - }}1)}} \varepsilon _{21i}^{(n)})$

(15)

式中： $\varepsilon _{31i}^{(n)}$ 为悬沙由悬移状态转移到静止状态的概率； $E_2^{(n)}$ 为淤积厚度参数，且 $E_2^{(n)} = \alpha {E^{(n)}}$ 。

从而得到，计算时段内床面总冲淤厚度 $\Delta {H^{(n)}}$ 为（假设冲刷为正，淤积为负）：

$\Delta {H^{(n)}} = \Delta H_1^{(n)} + \Delta H_2^{(n)}$

(16)

综上所述，联立式（11）、（14）～（16）可得，天然河流中床沙级配的计算公式为：

$\begin{aligned}[b] P_{{\rm{a}}i}^{(n)} = P_{{\rm{a}}i}^{(n - 1)} - {N_1}[P_{{\rm{b}}i}^{(n)} + P_{{\rm{s}}i}^{(n)}] + {N_2}[P_{{\rm{b}}i}^{(n - 1)}\varepsilon _{21i}^{(n)} + P_{{\rm{s}}i}^{(n - 1)}\varepsilon _{31i}^{(n)}]\\ \end{aligned} $

(17)

式中，

$N_1^{(n)} = \sum\limits_{i = 1}^{\max } {P_{{\rm{a}}i}^{(n - 1)}\varepsilon _{13i}^{(n)}} + \sum\limits_{i = 1}^{\max } {P_{{\rm{a}}i}^{(n - 1)}\varepsilon _{12i}^{(n)}} $

$N_2^{(n)} = \alpha (\sum\limits_{i = 1}^{\max } {P_{{\rm{s}}i}^{(n - 1)}\varepsilon _{31i}^{(n)}} + \sum\limits_{i = 1}^{\max } {P_{{\rm{b}}i}^{(n - 1)}\varepsilon _{21i}^{(n)}} )。$

4 枝城河段床沙交换计算及讨论 4.1 计算流程

1）利用枝城站2003—2009年实测流量过程、床沙级配、推移质级配、悬移质级配资料，将初始床沙级配 ${P_{0i}}$ ，第n–1步床沙级配 $P_{ai}^{(n - 1)}$ 及第n步水流条件（水深 ${h^{(n)}}$ 、流速 ${U^{(n)}}$ ）代入式（3），得第n步三态转移概率矩阵 ${{\varepsilon}} _{xy}^{(n)}$ 。

2）将第n–1步床沙级配 $P_{ai}^{(n - 1)}$ 、推移质级配 $P_{bi}^{(n - 1)}$ 及悬移质级配 $P_{si}^{(n - 1)}$ ，第n步转移概率矩阵 ${{\varepsilon}} _{xy}^{(n)}$ 分别代入式（12）、（13）得第n步推移质级配 $P_{bi}^{(n)}$ 、悬移质级配 $P_{si}^{(n)}$ ；代入式（14）、（15）得第n步冲刷厚度 $\Delta H_1^{(n)}$ 、淤积厚度 $\Delta H_2^{(n)}$ ，进而得到总冲淤厚度 $\Delta {H^{(n)}} = \Delta H_1^{(n)} + \Delta H_2^{(n)}$ 。

3）将1）和2）中计算的第n步转移概率矩阵 ${{\varepsilon}} _{xy}^{(n)}$ 、推移质级配 $P_{bi}^{(n)}$ 、悬移质级配 $P_{si}^{(n)}$ 、冲刷厚度 $\Delta H_1^{(n)}$ 及淤积厚度 $\Delta H_2^{(n)}$ 代入式（11），得到第n步床沙级配 $P_{ai}^{(n)}$ 。

4）将年内床沙交换计算结果取算术平均值，可得年际床沙级配及中值粒径变化过程。

4.2 年内床沙交换计算

三峡工程运用后，荆江河段枯水期缩短，中水期延长，最小流量增加^[26]，从而造成该河段年内水沙过程发生改变，年内床沙交换机理更为复杂，影响因素增加。在以上研究的基础上，作者采用枝城站2003—2009年实测流量过程、床沙级配、推移质级配、悬移质级配资料计算荆江河段年内床沙交换过程。

4.2.1 床沙交换过程分析

由前文叙述可知，长江中下游年内水沙条件多变，汛期一般集中在5—10月份，进而导致床沙交换过程比较复杂。为深入了解长江中游床沙交换过程，搜集整理了枝城站2003—2017年实测床沙资料，得到该断面年内月均床沙变化过程，如图7（a）、（b）、（c）所示。

通过统计分析，将枝城站年内月均床沙变化过程大致分为3类：1）平衡—粗化；2）平衡—细化—粗化；3）细化—粗化—细化—粗化。此外，2013年（细化—平衡—粗化—平衡—细化）及2015年（粗化—细化）变化过程如图7（d）所示。

由图7（a）、（b）、（c）可以看出，枝城站年内床沙变化复杂，但可总结出如下规律：

1）平衡—粗化：在汛期结束（9—10月）之前，床沙交换过程处于近似平衡状态，汛期之后随着水沙条件的变化逐渐粗化；

2）平衡—细化—粗化：在1—4月床沙交换处于近似平衡状态，然后床沙在4—5月迅速细化，5—12月完成粗化；

3）细化—粗化—细化—粗化：此过程比较复杂，各变化过程之间的时间节点不同，不作具体分析。

4.2.2 床沙交换过程计算

由于验证结果较为复杂，年内计算过程曲线较多，该断面选择具有代表性的2007年（平衡—细化—粗化）及2009年（细化—粗化—细化—粗化）床沙级配计算结果来说明本文概率模型的计算精度，如图8（a）、（b）所示。

由图8（a）可知，2007年床沙级配计算结果与前文分析一致，床沙交换过程近似为平衡—细化—粗化；其中，3月21日床沙级配与1月18日初始级配相近（近似平衡态），直至汛期开始之后（7月）完成细化，之后床沙逐渐粗化。

由图8（b）可知，2009年床沙交换过程计算结果体现为细化—粗化—细化—粗化，与前文分析一致；其中，1月15日为初始床沙，至汛期初第一次细化完成，之后床沙逐渐粗化，并在汛期末（10—11月）完成。

此外，为了说明本文概率模型年内床沙交换过程（床沙中值粒径d₅₀与时间t的变化关系），将枝城站2003—2009年年内床沙中值粒径d₅₀计算结果与实测资料进行对比，其年内床沙交换过程如图9所示。

由图9可见，2003—2009年，枝城断面年内床沙交换过程较为复杂，床沙中值粒径d₅₀变化幅度较大，包含了床沙交换过程的3种状态：近似平衡、粗化、细化，只是3种交换状态在每年发生的时段不同；除个别计算点外，本文模型计算结果与枝城站实测资料符合较好，变化趋势基本一致，能够应用于枝城河段卵石夹沙河床年内床沙交换过程的预报。此外，图9中个别数据点依然存在较大误差，可能是因为床沙取样时存在随机性，模型计算存在累计误差，或者来水来沙过程（江湖水沙交换等）变化复杂。

4.3 年际床沙粗化计算

结合本文第4.2节床沙交换过程的计算结果，运用枝城站2003—2009年实测流量过程、床沙级配、推移质级配、悬移质级配资料计算枝城断面年际床沙交换过程及变化规律，并将计算结果与实测资料进行对比，如图10所示（图10中：“P”代表“计算值”，“M”代表“实测值”）。

由图10可见，2003—2009年，该断面床沙发生不同程度的粗化，其中在2007年粗化程度最为明显，2010—2017年粗化速率慢于2003—2009年粗化速率；该模型能够较好的预测卵石夹沙河床年际床沙变化过程。另外，在计算年际床沙变化过程时，2003—2009年实测及计算床沙级配皆为该年日均床沙级配的算术平均值。

为了更直观的反映本文概率模型年际床沙交换过程的计算精度，将枝城站年际床沙中值粒径d₅₀计算结果与实测资料进行比较，验证结果如图11所示。

由图11可见：2003—2009年，枝城断面床沙逐年粗化，但粗化速率在2007年以后有变慢的趋势，但仍需后续补充实测资料进一步验证；此外，本文公式计算结果与实测资料符合较好，变化趋势一致，能够应用于三峡工程下游宜枝河段卵石夹沙河床年际床沙粗化过程的预报。

5 结　论

三峡工程运用后，下泄水流含沙量处于严重不饱和状态，下泄清水持续冲刷河床，造成长江中游河段床沙逐年粗化，对长江中下游河段水沙输移及河床演变规律产生重要影响。为研究宜枝河段、枝江河段卵石夹沙河床年际床沙粗化及年内床沙交换过程，基于Markov三态转移概率矩阵，引入非均匀沙隐暴系数，得到卵石夹沙河床上床沙级配、推移质级配及悬沙级配的概率计算模型，能够同时考虑前期水沙条件、床沙起悬及冲淤过程对床沙交换过程的影响。在此基础上，得到的主要结论如下：

1）三峡工程运用后，受持续冲刷的影响，宜枝河段床沙逐年粗化；2003—2017年，枝城站床沙逐年粗化（2007年粗化最为明显），泥沙颗粒级配分布曲线左移，2003年床沙粒径范围为0.002～41.000 mm，2017年变为0.031～86.000 mm；该断面床沙中值粒径总体呈上升趋势（中值粒径在2007年发生突变），并在2010年以后有一定波动，但波动幅度较小。同时，河床粗化层内并非全为粗颗粒泥沙，只是泥沙颗粒级配分布曲线发生了变化（分布曲线左移），床沙粒径范围逐渐变宽。

2）通过统计分析枝城站2003—2017年实测资料，以汛初、汛末为时间节点，将枝城站年内月均床沙变化过程大致分为三类（2013年及2015年除外）：“平衡—粗化”、“平衡—细化—粗化”以及“细化—粗化—细化—粗化”，并分析了各类床沙变化过程的特点。

3）通过三峡工程下游枝城站2003—2017年部分实测流量过程、床沙级配、推移质级配、悬移质级配资料验证可知：2003—2009年，枝城断面年内床沙交换过程较为复杂，床沙中值粒径d₅₀变化幅度较大，包含了床沙交换过程的三种状态：近似平衡、粗化、细化，只是三种交换状态在每年发生的时段不同；除个别计算点外，本文公式计算结果与枝城站实测资料符合较好，能够应用于宜枝河段卵石夹沙河床年际床沙粗化及年内床沙交换过程的预报，为今后研究三峡工程下游宜枝河段及枝江河段卵石夹沙河床非均匀沙运动及非平衡输沙机理提供理论基础。但是，由于缺少宜枝河段的详细实测资料（推移质、悬移质级配）验证，且长江中下游水沙条件多变，床沙交换过程复杂，其内在机理有待进一步研究。