许云鹏，吴震宇，尹 川

(四川大学水利水电学院 水力学与山区河流开发保护国家重点实验室，四川 成都 610065)

1 概述

土质心墙堆石坝是我国西部高坝广泛采用的坝型之一，已建成的长河坝(240m)、糯扎渡坝(261.5m)、两河口坝(293m)等心墙堆石坝工程运行良好，在建的古水坝(305m)、双江口坝(314m)是超过300m的世界级水平大型心墙堆石坝工程[1]。防渗心墙的渗透稳定性是影响大坝安全的关键因素。库水位的变动会影响心墙堆石坝的渗透稳定性。一般来说，堆石坝心墙由浸润线分成饱和区和非饱和区，渗透系数在饱和区域是一个常数[2- 3]。库水位变化时，总存在一个非饱和区域，其渗透系数随孔隙水压力和给水度的变化而变化。但在通常的心墙堆石坝渗流数值模拟计算中，视渗透系数为常数，上述因素对计算的影响往往被忽略[4- 5]。实践证明，在水库水位变化情况下计算土坝渗流时，如果不考虑心墙非饱和特性因素，计算结果和实际情况相比会有很大的误差[2，6]。

非饱和渗透系数可根据土壤水分特征曲线(SWCC)进行推定[6]。土壤水分特征曲线成为非饱和渗流数值模拟的重要基础。直接测定土壤水分特征曲线成本较高且十分复杂，经验公式法预测SWCC逐渐应用开来。经验公式法是利用代数关系式来描述SWCC，代表性的模型主要有Brooks-Corey模型(BC模型)、van Genuchten模型(VG模型)、Campbell模型、Fredlund模型等[7]。在众多SWCC模型中，因为BC模型和VG模型适用性范围更广，准确度更高，故而相关研究中得到了广泛的应用[8]。同时BC模型和VG模型提供渗透系数函数模型，用于表示渗透系数随基质吸力的变化情况。Gardner模型[9]仅提供渗透系数函数，可与VG模型的SWCC结合，用于非饱和渗流分析，称为VG-G模型。

土-水特征曲线和渗透系数函数可以准确反映材料饱和度随孔隙水压力的变化情况，进而反映材料渗透系数随孔隙水压力的变化情况，直接影响非饱和非稳定渗流模型的正确性。因此基于渗流实测数据如何合理确定SWCC和渗透系数函数模型，是构建非饱和渗流模型、开展后续渗流分析研究的关键。模型选择不仅要考虑模型的复杂程度，也要考虑模型的预测效果，过于简单的模型不能充分描述数据的内在结构变化，造成预测偏差，模型复杂虽然能充分描述数据的内在关联，但过于复杂的模型也会导致预测的方差偏大，因此选择合适的模型就是要找到预测偏差与方差的平衡点。1974年日本学者Akaike在信息理论研究中，为平衡估计模型的复杂性和拟合数据优良性，提出了AIC准则[10]。1978年Schwarz使用蒙特卡洛法研究AIC准则时发现，当样本数量增加时，AIC准则计算精度逐渐下降，进而在AIC准则中引入了包含样本数量的惩罚项，提出BIC准则以满足大样本条件下的模型识别[11]。随着模型识别方法的逐渐完善，越来越多的国内外学者围绕模型识别开展了深入研究。汪建均[12]在识别非正态响应部分因子试验的显著性因子中，结合经验贝叶斯先验和变量后验，提出了基于广义线性模型的模型选择方法。段偲默[13]根据AIC准则和BIC准则从基于Copula函数的风光联合出力的静态模型和基于Copula函数的风光联合出力的动态模型中获取最优模型。本文以PBG砾石土心墙堆石坝为例，基于原观数据进行非饱和渗流参数反演然后基于AIC、BIC准则识别非饱和渗流模型。

2 非饱和渗流模型

达西定律被证明既适合饱和土体，也适合非饱和土体。与饱和土体相比，非饱和土体渗透系数不是一个常数，它随着土体饱和度或基质吸力的变化而变化。根据达西定律及流体运动连续性方程可导出二维非稳定饱和—非饱和渗流控制方程为：

(1)

(2)

式中，kx、ky—x、y方向的渗透系数；H—总水头；Q—边界流量；γw—水的重度；mw—比水重度，定义为土水特征曲线斜率；θw—土体含水率；ua—孔隙气压力；uw—孔隙水压力。

对于稳定流条件，基质吸力剖面的时间依赖性消失，等式(1)的右侧变为零。以下介绍了本工作中用于描述饱和和非饱和流动的主要本构模型和参数。Brooks和Corey基于土体含水量试验提出了有效饱和度和基质吸力之间关系的土水特征曲线模型：

(3)

式中，h—为吸力水头；hb—进气压力水头。

Van Genuchten提出了一种形式平滑、封闭的3参数土水特征曲线模型如下：

Se=[1+(-ah)n]-m

(4)

式中，a—约等于进气压力值的倒数，kPa-1；m、n—无量纲拟合参数，m=1-1/n。

渗透系数随基质吸力变化而变化。BC模型和VG模型分别给出了渗透系数函数关系方程如下：

(5)

(6)

式中，Ks—饱和土渗透系数。

除上述2种模型外，Gardner提出了幂函数形式的双参数渗透系数模型：

K(h)=Ks[1+(-αh)β]-1

(7)

式中，α—渗透系数开始下降的起点；β—渗透系数随基质吸力下降的斜率。

将Brooks-Corey表达土水特征曲线的公式(3)和渗透系数的公式(5)称为Brooks-Corey-Burdine模型，简称为Brooks-Corey模型(BC模型)；将van Genuchten表达土水特征曲线的公式(4)和渗透系数的公式(6)称为van Genuchten-Mualem模型，简称为van Genuchten模型(VG模型)。Gardner模型仅提供渗透系数函数，可与VG模型的土水特征曲线结合，用于非饱和渗流分析，称为VG-G模型。

3 模型识别方法

3.1 模型识别准则

最优模型准则包括似然函数最大化和模型未知参数最少化。似然函数值越大模型拟合效果越好，但易导致未知参数越来越多，模型过于复杂，因此最优模型应该综合拟合精度和参数个数最优化。AIC准则全称是赤池信息准则(Akaike Information Criterion)，主要思想是平衡模型的复杂度和数据的拟合精度，其表达式如下：

AIC=2k-2lnL

(8)

式中，k—参数个数，L—极大似然函数。

当样本容量较大时，AIC准则会选择参数较多的模型，造成过拟合现象。针对该问题，Schwartz提出BIC准则，全称是贝叶斯信息准则(Bayesian Information Criterion)，其表达式如下：

BIC=klnN-2lnL

(9)

式中，N—样本数量，其他参数意义同上式。

BIC准则考虑了样本数量，避免精度过高选择复杂模型。使用AIC、BIC准则进行模型识别时，值越小模型越优。

3.2 模型识别方法流程

首先选择待反演的堆石坝心墙非饱和区渗透系数、饱和区渗透系数、非饱和渗流模型特性参数作为待反演参数。根据渗透系数和非饱和模型特性参数取值区间，在材料正常的取值范围内上下波动参数进行正交试验，得到多组试验参数组合。根据坝体结构轮廓尺寸、材料分区和坝址地形地质条件等资料，构建二维心墙堆石坝有限元模型，输入实测库水位过程线，计算各组试验组合下非稳定非饱和渗流数值模拟，输出每组试验下的测点模拟渗压监测序列。以最低库水位为界，将心墙分为饱和区和非饱和区。根据上游水位变化，选取代表性水位对应时间节点作为渗压控制节点。以待反演参数为自变量，各组有限元模型控制节点下计算渗压值为因变量拟合响应面方程。基于最小二乘法原理，利用心墙实际测点渗压监测序列和拟合的响应面方程构建目标函数。将构造的多目标函数输入到NSGA-II算法中进行优化，输出Pareto最优解集，最终确定反演参数。反演参数代入有限元模型，计算得到心墙拟合渗压值，结合实测渗压值，采用AIC、BIC准则对每个模型进行评价，最终筛选出最优模型。模型识别方法流程如图1所示。

图1 模型识别方法

4 工程实例

4.1 有限元模型及参数

PBG大坝坝体断面尺寸和材料分区如图2所示，选取大坝渗流监测河床0+240断面，建立坝-地基非饱和渗流二维有限元模型，如图3所示。坝体最大坝高186m，坝顶宽度14m，上游坝坡比为1∶2和1∶2.25，下游坝坡比为1∶1.8。心墙顶宽4m，底宽96m，上、下游坡比均为1∶0.25；在心墙上、下游侧设置双层的反滤层，反滤层与堆石间设过渡层。坝基下有3层覆盖层，最大厚度75.4m，具有中等-强透水性。坝基设主、副两道各厚1.2m的混凝土防渗墙，墙中心间距14m，墙底嵌入基岩1.5m。主防渗墙底帷幕深度为50m；副墙在主墙上游侧，墙顶插入心墙内部10m，墙底帷幕深度为10m。其地基模拟范围为：上下游方向自坡脚延伸2倍坝高，深度方向延伸2倍坝高。离散后的坝-地基体系有限元模型含11555个单元，其中坝体5912个单元，防渗墙及墙底帷幕102个单元，坝基覆盖层①816个单元，覆盖层②1311个单元，基岩3414个单元。坝体坝基各类材料的渗透系数见表1。根据PBG工程心墙渗压计的布置和运行情况，选取0+240断面测点分析，测点布置如图4所示，其中正常在测测点有P12、P13、P15、P16、P19、P20、P21。

4.2 非饱和渗流参数反演

通过敏感性分析，结合上述提到的3种非饱和渗流模型，确定需要反演的参数有心墙上区渗透系数、心墙下区渗透系数、弱风化基岩渗透系数以及非饱和特性参数。输入2014年PBG工程运行时库水位历时过程线，每2个月的一天作为时间控制点，共6个控制节点。3种模型的参数分别进行正交试验设计，输入有限元中计算得到各组正交试验的渗压监测序列。采用不含交叉项的3次多项式，各时间控制节点试验渗压进行响应面方程拟合。各测点基于渗压实测值和响应面拟合值的误差平方和构建目标函数，采用多目标优化对多个目标函数求解，筛选得到合理反演参数。将反演参数代回有限元模型进行非饱和非稳定渗流计算，同时不考虑非饱和特性进行非稳定渗流计算。反分析误差分析结果见表2，拟合效果选取心墙上区测点P12和心墙下区测点P15、P16进行展示，如图5—7所示。

图2 PBG工程河床0+240断面图

图3 PBG0+240断面各分区有限元网格

表1 坝体及坝基各区材料渗透系数 单位：cm/s

图4 PBG堆石坝心墙特征点示意图

对比非饱和非稳定渗流与饱和非稳定渗流，非饱和特性对心墙下区渗压影响较小，2种渗流的心墙下区测点的平均绝对误差分别为1.18m(VG模型)和1.97m，但是对心墙非饱和区渗压影响较大，2种渗流的平均绝对误差分别为3.53m(VG-G模型)和8.75m，饱和非稳定渗流的最大绝对误差达到35.50m。因此非饱和渗流有限元模型不仅能反映饱和区渗流变化，也能反映非饱和区渗流变化，使心墙渗流的模拟更符合实际情况。同时在3种模型的对比中，对于心墙上区测点，VG模型、VG-G模型和BC模型的平均绝对误差分别为2.83、3.53、2.53m，对于心墙下区测点，VG模型、VG-G模型和BC模型的平均绝对误差分别为1.18、0.98、0.90m，因此BC模型对心墙饱和区和非饱和区的渗压变化的拟合情况更加准确。

4.3 模型识别

基于3种非饱和渗流模型计算的拟合渗压值和心墙实际渗压值，分别采用AIC和BIC准则对VG模型、BC模型、VG-G模型3种非饱和渗流模型进行识别，AIC和BIC准则计算结果见表3，如图8所示。VG模型、VG-G模型、BC模型的AIC均值分别为604.48、522.96、374.06，BIC均值分别为619.40、543.84、388.97，2种准则结果一致，均表明3种模型中BC模型最优，VG-G模型稍差，VG模型最差。

表2 反演结果误差分析 单位：m

图5 P12 VG、VG-G和BC模型反分析及非稳定渗压对比

图6 P15 VG、VG-G和BC模型反分析及非稳定渗压对比

单从非饱和区测点P12分析，VG-G模型的AIC值为689.09，大于BC模型的298.98，小于VG模型的787.32，表明对于非饱和区，BC模型最优，VG-G模型稍差，VG模型最差。饱和区测点模型识别结果表明，BC模型和VG-G模型的AIC和BIC值均接近，平均相差28.89，VG模型更劣于VG-G模型，平均相差73.17。总而言之，3种模型对心墙饱和区的拟合差异远小于心墙非饱和区，BC模型对于心墙渗压变化的拟合情况更加准确，更能整体反映大坝心墙渗压场的变化情况。

图7 P16 VG、VG-G和BC模型反分析及非稳定渗压对比

图8 非饱和渗流模型识别结果图

表3 非饱和渗流模型识别结果

5 结语

本文以PBG砾石土心墙堆石坝0+240断面为基础，以SWCC和渗透系数函数为基础，构建了3种非饱和非稳定渗流模型，结合实测渗压值和反演参数代入有限元计算得到模型拟合渗压值使用AIC、BIC准则对非饱和渗流模型进行最优模型识别。结果表明，3种非饱和渗流模型的拟合渗压值与实测值误差较小，可以反映心墙渗压变化；同时在这3种模型中，反演结果误差分析与AIC、BIC准则计算结果一致，均表明BC模型最为优越，对于心墙渗压变化拟合效果更加准确，更能整体反映堆石坝心墙渗压场的实际变化情况。