快速评估溃决洪水是堰塞湖应急处置中的一项重要内容，当前的堰塞坝溃坝模拟软件普遍存在结果分散、计算不稳定等问题，其中，堰塞体冲刷参数是影响模拟结果的关键因素，合理地选取冲刷模型参数是准确计算溃决洪水的前提。通常情况下，冲刷模型参数依据试验确定，但在堰塞湖应急处置条件下，此方法用时较长，难以满足应急需求，故需要提出一种快速的冲刷模型参数确定方法。

当前应用最为广泛的冲刷模型是过剪应力模型和Wilson模型。过剪应力模型是由Partheniades^[1]最先提出，表达式为：

$ \dot{\textit{z}}={\textit{z}}_{0}{\bigg(\frac{\tau }{{\tau }_{{\rm{c}}}}-1\bigg)}^{\alpha } $

(1)

式中：τ_c、z₀和α为过剪应力模型的冲刷参数，其中，τ_c为起动切应力，z₀为冲刷系数，α为无量纲指数；τ为水流施加给土体的切应力， $ \dot{\textit{z}}$ 为土体冲刷速率。

钟启明等^[2]采用过剪应力模型模拟了每层堰塞体材料的冲蚀过程，表明该模型能较好地揭示堰塞体受水流冲刷的作用机制。将过剪应力模型的冲刷参数τ_c、z₀和α放入模型中可较好地体现堰塞体材料的冲刷机制。Shields^[3]通过大量水槽试验提出了Shields图，用以表征土体中值粒径和起动切应力的关系，至今仍被广泛使用，但Shields图中黏性土颗粒粒径与起动切应力之间没有明显的相关性。宗全利等^[4]通过对荆江段河岸黏性土体起动条件研究，得到了黏性土起动切应力与土体液性指数之间的关系和冲刷系数z₀与起动切应力之间的关系。洪大林等^[5]对中运河、长江长兴岛、淮河入海水道等河道原状土样开展起动试验，表明对于不同的河流，由于土质结构、形成机制的不同，尽管原状土样具有相同的物理特性，但起动条件有较大差别；对于同一条河流，起动条件与土体物理力学指标有明显的规律性。

另一个当前广泛应用的模型是Wilson模型。Wilson^[6-7]基于土体水下受力规律和牛顿第二定律推导出的冲刷速率模型为：

$ \dot{\textit{z}}=\frac{{b}_{0}\sqrt{\tau }}{{\rho }_{{\rm{b}}}}\left[1-\mathrm{exp}\left\{-\mathrm{exp}\left(3-\frac{{b}_{1}}{\tau }\right)\right\}\right] $

(2)

式中，b₀、b₁为Wilson模型的冲刷参数，ρ_b为土体堆积密度。

Al–Madhhachi等^[8-9]对两种黏性土进行了12组冲刷试验，提出随着含水量的增加冲刷参数b₀与冲刷参数z₀都有下降趋势；并研究了渗透力对黏性土冲刷性能的影响，指出随着渗透力的增加，冲刷参数b₀变大而冲刷参数b₁减小。Criswell^[10]对平均粒径分别为0.45、0.60、1.30、1.90 mm的砂土进行了冲刷试验，建立了冲刷参数b₁–τ_c、b₀–b₁和z₀–τ_c之间的关系式。

Knapen等^[11]总结了过剪应力模型中冲刷参数与土性参数的关系，并研究了冲刷参数的范围；同时，Knapen等提到不同学者得到的关系不一致，甚至有些相反的结论，比如：Alberts^[12]提出冲刷系数与黏粒含量呈正相关关系，而Rapp^[13]和Allen^[14]等提出冲刷系数与黏粒含量呈负相关关系。因此，需要进一步对模型中冲刷参数τ_c、z₀和α的取值范围和变化规律进行深入的研究。

作者利用收集的冲刷试验数据，总结得到了各类土体冲刷参数的范围，并通过回归分析和相关性分析得出了冲刷参数与土性参数之间的关系式。在堰塞湖应急处置时，可依据现场确定的堰塞体物理性质，应用所提出的关系式，较为快速、准确地确定冲刷模型参数，为快速评估溃决洪水提供依据。

1 数据来源和分析方法 1.1 数据来源

堰塞体的材料类型较为广泛，因此，在收集已有试验资料时，尽可能覆盖堰塞体的级配。通过技术报告、国内外期刊^[15-23]及作者开展的室内试验^[24]，总结出了279组冲刷试验数据的数据库。将数据分成粗、细颗粒土两类，其中，粗颗粒土（粒径大于0.075 mm的颗粒超过全重50%）84组，细颗粒土（粒径大于0.075 mm的颗粒不超过全重50%）195组。图1为数据库中包含的常规土性参数频率直方图，包括塑性指数（PI）、重度（ $\gamma $ ）、含水量（w）、细颗粒含量（P_f）、黏粒含量（P_c）、抗剪强度（τ_f）、平均粒径（D₅₀）、不均匀系数（C_u）和曲率系数（C_c）。

将起动切应力定义为当冲刷速率小到0.1 mm/h时所对应的切应力^[25]，冲刷参数z₀等于切应力达到两倍起动切应力时所对应的冲刷速率。若冲刷速率曲线上存在切应力等于两倍的起动切应力的点，则该点所对应的冲刷速率的值为z₀；若冲刷速率曲线不存在该点，则对切应力小于和大于两倍起动切应力的点采用线性插值的方法得到两倍起动切应力点所对应的冲刷速率。在起动切应力τ_c和冲刷参数z₀确定之后，采用拟合曲线工具得到指数系数α。通过使用拟合工具也可以类似地得到冲刷参数b₀、b₁。图2和3分别列举了过剪应力模型和Wilson 模型的冲刷参数频率直方图。

1.2 分析方法

在进行回归分析之前，需要了解冲刷参数与常规土性参数的关系。应用相关系数r衡量两者之间的相关性，r的范围从–1到+1，负值代表负相关，正值代表正相关，0代表没有相关性。

$ r=\frac{{ \displaystyle\sum ({X}_{i}-\overline{X})({Y}_{i}-\overline{Y})}}{\sqrt{{\displaystyle \sum {({X}_{i}-\overline{X})}^{2}}}\sqrt{{ \displaystyle\sum {({Y}_{i}-\overline{Y})}^{2}}}} $

(3)

式中，X_i和Y_i为自变量， $ \overline{X}$ 和 $ \overline{Y}$ 为平均值。为避免偶然因素，利用表1中的概率P值检测结果的显著性，P可通过统计性软件求得，若小于显著性水平（0.05），则认为结果在统计意义上显著。

通过回归分析得到冲刷参数与常规土性参数的回归方程。通过4个评价指标，即R²、MSE、F_value/F_statistic和CV衡量回归方程的准确性。其中，R²为决定系数，是回归线对各观测点拟合紧密程度的测度，数学定义为：

$ {R}^{2}=1-\frac{{ \displaystyle\sum {\left({y}_{i}-{f}_{i}\right)}^{2}}}{{ \displaystyle\sum {\left({y}_{i}-\overline{y}\right)}^{2}}} $

(4)

式中，y_i为试验值， $ \overline{y}$ 为试验点的平均值，f_i为回归方程的预测值。R²的范围在0～1之间，越接近1，说明拟合效果越好，通常认为R²超过0.6便可接受^[22]。均方误差MSE为另一个判断拟合优度的参数，反映了预测值与试验值之间方差的平均值，表达式如下：

$ MS\!E=\frac{1}{n}{{ \sum \left({y}_{i}-{f}_{i}\right)}}^{2} $

(5)

式中，y_i为试验值，f_i为预测值。MSE的值越小，说明模型的预测精确越高。

R²和MSE并不能完全代表回归方程的统计学意义。比如：R²过于依赖数据点的数量，当预测一个参数时，如果仅有一个数据点，那么任何回归模型所得到的R²均为1。为避免这种情况，需引进另外一个参数F_value：

$ {F}_{\rm{value}}=\frac{S\!S\!E(reg)}{S\!S\!E(res)}\times \left(\frac{D}{V}-1\right)=\frac{{ \displaystyle\sum ({f}_{i}-\overline{y}{)}^{2}}}{{ \displaystyle\sum {({y}_{i}-{f}_{i})}^{2}}}\times \left(\frac{D}{V}-1\right) $

(6)

式中，D为数据点个数，V为自变量的个数，SSE（reg）为回归平方和，SSE（res）为残差平方和。从式（6）可以看出，如果数据点的个数远大于自变量的个数，残差平方和较小时，F_value将会变得很大。如果F_value远大于F检验的临界值，可以判断回归方程的可信度较高。F临界值可通过查找F分布表得到，这里显著性水平取0.05。

过拟合是在工程领域进行统计学分析时常遇到的问题，比如，用全部的数据拟合得到一个等式，就会发生过拟合现象。过拟合会导致增加一个数据或减少一个数据，公式的拟合优度会发生很大改变。为得到泛化能力强的等式，需对候选模型进行K折交叉检验^[26]。K折交叉检验的定义为将数据集随机分成K份容量相同的子集，依次选第k份子集作为测试集，其余的K–1份作为训练集，最后对K份泛化误差（CV）进行平均，便得到了K折交叉验证的CV值，如式（7）所示：

$ C{V}_{k}=\frac{1}{n}{ \sum ({y}_{i}}-{f}_{i}^{-k}{)}^{2}, CV=\frac{1}{K}{ \sum C{V}_{k}} $

(7)

式中：CV_k为用除第k份以外的样本作为训练集进行训练，之后用第k份对训练结果进行预测得到的泛化误差；f_i^–k表示去掉第k份，在其余K–1份拟合得到的预测值。本文采用10折交叉验证(K=10)。10折交叉验证选择模型的方法是在候选模型中选择使得泛化误差（CV）最小的模型。

所建立的回归方程应使R²接近1，MSE较小，F_value/F_statistic远远大于1，CV值较小。通过以下3个步骤筛选得到较好模型：1）综合考虑R²和MSE，在所有候选模型中挑选出拟合优度较好的几组模型。2）为了得到可信度较高的模型，从第一步筛选得到的候选模型中挑出F_value/F_statistic较大的2组模型。3）为了得到泛化能力较好的模型，对这2组模型进行10折交叉检查，通过比较CV值，选出最小的模型，从而得到最终回归方程。

2 冲刷参数关联性分析 2.1 细颗粒土的相关分析

细颗粒土冲刷参数与常规土性参数之间的相关性如表1所示，若P值小于0.05，则回归结果显著。细颗粒土的冲刷参数b₀与塑性指数和黏粒含量呈显著负相关，与平均粒径呈显著正相关；细颗粒土的冲刷参数b₁与含水量和细颗粒含量呈正相关；细颗粒土冲刷参数z₀与塑性指数、黏粒含量呈显著负相关，与平均粒径呈显著正相关；细颗粒土的起动切应力τ_c与含水量、细颗粒含量呈显著正相关，与不均匀系数呈显著负相关。黏聚力是描述土体抗冲刷特性的重要参数，受到土体的组成、土体的密实度和含水量等因素影响^[10]。颗粒粒径降低、细颗粒含量增加和黏粒含量增加，会导致颗粒的比表面积增加，在静电力和范德华力的作用下，使土体黏聚力增大。塑性指数越大，表明土颗粒越细，颗粒的比表面积增加，从而使土体黏聚力增加。土体含水量小于最优含水量时，干密度会随着含水量的增加而增加，干密度反映了土体的密实度，其增加会引起土体黏聚力的增加。

2.2 粗颗粒土的相关分析

粗颗粒土冲刷参数与常规土性参数之间的相关性如表2所示。

粗颗粒土的冲刷参数b₀与重度、细颗粒含量、不均匀系数和黏粒含量呈显著负相关；冲刷参数b₁与黏粒含量呈显著正相关；粗颗粒土的冲刷参数z₀与重度、抗剪强度、细颗粒含量、不均匀系数呈显著负相关；粗颗粒土的起动切应力τ_c与抗剪强度、黏粒含量呈显著正相关；冲刷参数α与不均匀系数C_u呈显著负相关。土体重度的增加、细颗粒含量和黏粒含量的增加，同样引起土体比表面积的增加，在静电力和范德华力的作用下，土体黏聚力增大。随着不均匀系数的增加，土颗粒之间嵌固更为密实，从而引起黏聚力增加^[27]。

通过对粗颗粒土和细颗粒土的冲刷参数（z₀、τ_c、b₀、b₁）与土性参数的相关性分析可以得出：细颗粒土的冲刷参数b₀和z₀随着塑性指数、黏粒含量的减小而增大，随着平均粒径的增加而增加；细颗粒土的冲刷参数b₁和τ_c随着含水量、细颗粒含量的增加而增加。粗颗粒土的冲刷参数b₀、z₀随着重度、抗剪强度、细颗粒含量、不均匀系数的增加而减小。

2.3 细颗粒土的回归分析

通过分析细颗粒土和粗颗粒土的冲刷参数z₀、τ_c、b₀、b₁与土性参数的相关性可以发现，冲刷参数与单个土性参数的相关性并不理想，故本文选择多元幂函数模型进行回归分析，表达式为：

$ Y=A\times {\left({P}_{1}\right)}^{{\beta }_{1}}\times {\left({P}_{2}\right)}^{{\beta }_{2}}\times {\left({P}_{3}\right)}^{{\beta }_{3}}\times \cdot \cdot \cdot \times {\left({P}_{n}\right)}^{{\beta }_{n}} $

(8)

式中：Y为因变量；P₁、P₂、P₃为所选择的土性参数，包括PI、 $\gamma $ 、w、P_f、τ_f、D₅₀、P_c、C_u和C_c；A和 ${{\;\beta }_{i}} $ 为常数。基于第1.2节的分析方法，根据R²和MSE筛选得出拟合优度较好的几组模型。在确定冲刷参数b₁表达式时，发现候选模型的R²较小，无法满足拟合模型所需要的拟合优度，Wilson通过推导公式得出冲刷参数b₁与起动切应力存在一定的关系，因此本文通过研究冲刷参数b₁与起动切应力关系得到回归方程。表3～7分别列举了细颗粒土冲刷参数τ_c、z₀、α、b₀和b₁拟合优度较好模型的表达式、数据点个数、R²、MSE、F_value/F_statistic和CV值。然后，通过比较F_value/F_statistic值，在拟合优度较好模型中选出两组可信度较高的模型。最终对这两组模型进行10折交叉验证，通过比较CV值，得到冲刷参数的回归方程。

2.4 粗颗粒土的回归分析

根据Briaud等^[15]的研究，当粗颗粒土的平均粒径大于0.3 mm时，起动切应力与平均粒径呈线型关系。因此，本文在研究粗颗粒土起动切应力时，将数据分成D₅₀<0.3 mm和D₅₀>0.3 mm两部分。研究D₅₀>0.3 mm数据起动切应力的定量关系式时，土性参数只需选择平均粒径预计最佳的回归方程。研究D₅₀<0.3 mm数据的起动切应力定量关系式的方法与研究细颗粒土相同。

根据R²和MSE在候选模型中筛选得到拟合优度较好的模型，表8～12分别列举了粗颗粒土冲刷参数τ_c、z₀、α、b₀和b₁拟合优度较好模型的表达式、数据点个数、R²、MSE、F_value/F_statistic和CV值；由F_value/F_statistic和CV值，得到冲刷参数的回归方程。

图4和5分别为通过试验数据得到的冲刷参数与通过表达式预测得到的冲刷参数散点图。细颗粒土的起动应力τ_c随着重度、抗剪强度和黏粒含量的增加而增加，这与工程实际是一致的；冲刷参数z₀随着重度、含水量和抗剪强度的增加而增加；冲刷参数α随着含水量、细颗粒含量、黏粒含量的增加和抗剪强度的减小而增加；冲刷参数b₀随着含水量和重度的增加而增加。

对于粗粒土，当平均粒径小于0.3 mm时，起动切应力随着不均匀系数的减小而增加，当平均粒径大于0.3 mm时，起动切应力随着平均粒径的增加而增加；冲刷参数z₀随着含水量和曲率系数的增加而增加；冲刷参数b₀随着塑性指数的减小和平均粒径的增加而增加；冲刷参数b₁随着抗剪强度的增加而增加。对于细颗粒土，黏聚力是描述土体抗冲蚀特性的重要指标，重度增加、平均粒径降低、细颗粒含量增加、塑性指数增加和黏粒含量增加，会导致颗粒的比表面积增加，在静电力和范德华力的作用下，土体的黏聚力增大，从而使细颗粒土的抗冲蚀性能增强。对于粗颗粒土，土体的冲蚀特性与颗粒的重量有着密切的关系，颗粒粒径的增加，导致土颗粒重量的增加，使粗颗粒土的抗冲蚀特性增加。

2.5 冲刷参数的范围

粗、细颗粒土的过剪应力模型和Wilson模型的冲刷参数范围分别如图6和7所示。细颗粒土的冲刷参数z₀、指数系数α和起动切应力的范围分别为0.89～5.77 mm/h、0.57～1.16、0.5～5.2 Pa；粗颗粒土的冲刷参数z₀、指数系数α和起动切应力的范围分别为1.5～28.0 mm/h、0.77～1.74、0.2～1.9 Pa。细颗粒土的冲刷参数b₀和b₁的范围分别为1.03～7.50 gm^–1s^–1N^–0.5、4.02～45.89 Pa；粗颗粒土的冲刷参数b₀和b₁的范围分别为10.85～642.59 gm^–1s^–1N^–0.5、4.02～45.89 Pa。

为更细致地了解不同土料冲刷参数的范围，将细颗粒土细分为黏土和粉土（细粒组含量≥50%且粗颗粒含量≤25%；若PI ≥0.73（w_L–20）且PI ≥7，为黏土；若PI<0.73（w_L–20）且PI<4，为粉土；w_L为土体液限），所收集的数据包括黏土167组、粉土28组。将粗颗粒土分为砂类土（砾类组含量≤50%）和砾类土（2 mm≤d≤60 mm，砾类组含量>50%）。如图8所示：黏土的冲刷参数z₀、指数系数α和起动切应力的范围分别为0.82～3.88 mm/h、0.91～1.24、0.77～5.05 Pa；粉土的冲刷参数z₀、指数系数α和起动切应力的范围分别为1.62～21.46 mm/h、0.59～1.64、0.23～1.92 Pa；砂类土冲刷参数z₀、指数系数α和起动切应力的范围分别为1.50～26.01 mm/h、0.76～1.44、0.21～1.11 Pa；砾类土冲刷参数z₀、指数系数α和起动切应力的范围分别为2.3～120.9 mm/h、0.83～2.03、6.50～10.46 Pa。如图9所示：黏土的冲刷参数b₀和b₁的范围分别为1.93～9.10 gm^–1s^–1N^–0.5、13.79～70.85 Pa；粉土冲刷参数b₀和b₁的范围分别为2.63～113.48 gm^–1s^–1N^–0.5、2.59～44.83 Pa；砂类土冲刷参数b₀和b₁的范围分别为9.29～155.08 gm^–1s^–1N^–0.5、3.85～37.58 Pa，砾类土的冲刷参数b₀和b₁的范围分别为2.62～274.12 gm^–1s^–1N^–0.5、70.85～105.76 Pa。

为得到不同类型土的抗冲刷性能，需要综合考虑冲刷参数z₀和起动切应力或冲刷参数b₀和b₁的大小。图10（a）为通过综合考虑冲刷参数z₀和起动切应力的平均值得到的不同类型土的冲刷特性，其中，黏土、粉土、砂类土和砾类土的冲刷参数z₀平均值分别为8.27、20.34、28.74、87.80 mm/h，黏土、粉土、砂类土和砾类土的起动切应力平均值分别为4.73、2.44、1.85、7.92 Pa。由图10（a）可知：τ_c（砾类土）>τ_c（黏土）>τ_c（粉土）>τ_c（砂类土），z₀（砾类土）>z₀（砂类土）>z₀（粉土）>z₀（黏土）。因此，黏土的抗冲刷性能较好，其次是粉土，砂类土抗冲刷性能较差；砾类土的起动切应力最大，冲刷参数z₀也最大。图10（b）为通过综合考虑冲刷参数b₀和b₁平均值得到的不同类型土的冲刷特性，其中，黏土、粉土、砂类土和砾类土的冲刷参数b₀平均值分别为14.85、92.46、375.98、391.9 gm^–1s^–1N^–0.5，黏土、粉土、砂类土和砾类土的冲刷参数b₁平均值分别为53.42、49.65、28.81、89.49 Pa。由图10（b）可知：b₀（砾类土）≈b₀（砂类土）>b₀（粉土）>b₀（黏土），b₁（砾类土）>b₁（黏土）>b₁（粉土）>b₁（砂类土）；同样得到黏土的抗冲刷性能较好，粉土中等，砂类土抗冲刷性能较差。砾类土和砂类土的冲刷参数b₀相同，但砾类土的冲刷参数b₁大于砂类土的冲刷参数b₁，故砾类土的抗冲刷性能比砂类土要好。

3 白格堰塞体冲刷参数

在应急条件下，土体的颗粒级配可以快速得出，而受事故现场条件的限制，无法将第2节数据库中的土性参数完全检测出来，第2.3和2.4节给出的细、粗颗粒土冲刷参数公式的应用也受到限制，因此需要根据现场条件，依据第1.1节中的数据库及第1.2节的分析方法，得到冲刷参数与有限的土性参数的回归方程，从而快速确定冲刷参数的大致范围。

2018 年 10 月 11 日和 11 月 3 日，在西藏自治区江达县波罗乡白格村与四川省白玉县绒盖乡则巴村交界处金沙江西藏岸（右岸）先后两次发生大规模高位滑坡，堵塞金沙江，形成堰塞湖^[28]。图11为白格堰塞体土样的颗粒级配曲线，平均粒径D₅₀为5 mm，不均匀系数C_u为46.73%，曲率系数C_c为1.43%，细颗粒含量P_f为3.2%。由于相比过剪应力模型，Wilson模型的物理机理更加明确，故应用Wilson模型计算白格堰塞体的冲刷速率曲线。将冲蚀参数b₀与土性参数D₅₀、C_u、C_c和P_f的所有候选模型罗列出来，依据第1.2节的分析方法，得到冲刷参数b₀与土体的不均匀系数和曲率系数的关系式：

$ {b}_{0}=\frac{3\;162}{{{C}^{0.54}_{\rm{u}}}\times {C}^{0.23}_{{\rm{c}}}{}}$

(9)

预测白格堰塞湖坝址处土样的冲刷参数b₀为403.8 gm^–1s^–1N^–0.5，第2.5节给出了粗颗粒冲刷参数b₁的参考范围为4.02～45.89 Pa，本节b₁取50 Pa。陈祖煜等^[29]对白格堰塞体材料进行了冲蚀特性试验，本文预测的白格堰塞湖冲刷速率曲线与实测试验结果对比如图12所示。可以看出，与实测的堰塞体材料的冲蚀速率比较，预测所得的数值小于实测值，可能是由于实际堰塞体的密实度较低，比较松散，所得到的冲蚀速率较大。

4 结　论

快速确定堰塞体冲刷参数是堰塞湖应急处置中溃决洪水评估的基础。本文基于收集的279组冲刷试验数据，应用相关分析和回归分析方法分析了冲刷参数与土性参数的相关性，给出了黏土、粉土、砂类土和砾类土冲刷参数的大致范围，得到了粗、细颗粒土拟合优度（R²和MSE）、可信度（F_value/F_statistic）和泛化能力（CV）较好的冲刷参数回归方程，为快速确定堰塞体冲刷参数奠定了基础。

研究结果表明：细颗粒土的冲刷参数b₀和z₀随着塑性指数、黏粒含量的减小而增大，随着平均粒径的增加而增加；冲刷参数b₁和τ_c随着含水量、细颗粒含量的增加而增加。粗颗粒土的冲刷参数b₀、z₀随着重度、抗剪强度、细颗粒含量、不均匀系数的增加而减小。

基于本文统计分析结果，以白格堰塞体为例，对其冲刷参数和冲刷速率进行了快速分析。依据本文建立的冲刷参数数据库及统计关系，可通过堰塞体的颗粒级配，快速确定堰塞冲刷参数，为应急处置条件下快速评估溃决洪水风险提供依据。