一种高心墙堆石坝变形参数反演方法与流程

本发明涉及高心墙堆石坝监测技术领域，特别涉及一种高心墙堆石坝变形参数反演方法。

背景技术：

堆石坝是用块石、砂砾石等作为主体材料，经抛填或碾压而成的土石坝。坝体由作为支承的堆石和防渗体及他们之间的过渡层、反滤层组成。土质心墙可以与坝基内垂直和水平防渗体系相连接，因此，由土质心墙作防渗体的堆石坝，可以在深厚覆盖层上修建，是高、中坝常用的坝型。

在我国西部水电开发中，在建和已建的很多高心墙堆石坝所处的地形、地质环境复杂，坝体与地基载荷大、应力水平高，大坝与基础变形机理复杂，坝体开裂风险高，调控难度大。对于正在运行的心墙堆石坝，由于坝体和坝基材料的蠕变特性，大坝的变形往往历时多年才会达到稳定。实际工程中虽然布置了变形监测点，能够真实反映大坝的变形特性，然而监测点数量有限，只能捕捉到大坝局部变形特征。因此如何全面反映大坝整体变形特性，并准确地预测大坝变形发展，从而提前做好变形调控措施至关重要，一般采用数值模拟技术来实现这个目的。

数值模拟是依靠电子计算机，结合有限元软件，通过数值计算和图像显示的方法，实现对工程问题进行研究的一种方法。在此项目中，通过数值模拟技术可以全面反映大坝整体变形特性，并预测大坝变形发展规律，从而更好地调控大坝运行。然而在数值模拟过程中，选取的高心墙堆石坝数值模拟参数是否合适直接关系到模拟结果的准确性，而这些参数一般是通过室内试验或者工程类比的方法获得，室内试验由于尺寸效应以及样本与现场坝体的真实受力状态差别较大，并不能代表坝体结构的真实参数。因此，通过反演分析技术获得坝体结构的真实参数是一种有效途径。

目前，高心墙堆石坝变形参数的主要反演分析方法是以试验参数或类比值为基准值，在变形参数范围内改变参数取值，然后代入有限元计算，并将计算得到的结果与监测值作比较，计算结果与监测值越接近，取的参数越接近真实参数，也就是反演参数值。每取一组参数，就要进行一次高心墙堆石坝变形有限元计算，因此，该方法需要大量的有限元计算，效率较低。另外，根据高心墙堆石坝的运行特点，坝体结构产生的总变形是瞬时变形和蠕变变形叠加的结果，而目前反演分析有限元计算中，往往只计算瞬时变形，把瞬时变形当成总变形与监测值作比较进行反演参数优化，从而导致反演得到的瞬时变形参数偏大。

技术实现要素：

本发明的目的在于：提供了一种高心墙堆石坝变形参数反演方法，解决了传统高心墙堆石坝变形参数反演分析方法中的两点弊端：一是大量的有限元计算造成的反演效率低的问题；二是反演分析有限元计算中，未计算蠕变变形，把瞬时变形当成总变形与监测值作比较进行反演参数优化，从而导致反演得到的瞬时变形参数偏大的问题。

本发明采用的技术方案如下：

一种高心墙堆石坝变形参数反演方法，包括以下步骤：

s1：结合心墙堆石坝的地质剖面图和变形监测布置图，建立心墙堆石坝三维有限元模型，并根据沉降仪埋设位置，在有限元模型中相应位置布置节点；

s2：进行参数敏感性分析，即从若干材料力学参数中选取对变形敏感的参数，即得到待反演参数；

s3：运用均匀设计方法设计待反演参数样本；

s4：将步骤s3中得到的待反演参数样本输入步骤s1中建立的三维有限元模型中，进行有限元计算，得到变形值样本；

s5：将步骤s4中得到的变形值样本进行标准化处理，得到标准化变形值，建立标准化变形值与待反演参数样本之间不含交叉项的响应面函数，并根据最小二乘法原理拟合响应面函数的系数；

s6：结合心墙堆石坝原型观测数据，利用响应面函数，选取多个不同测点，建立目标优化函数；

s7：引入模拟退火算法对步骤s6中建立的目标优化函数进行迭代优化，得到满足心墙堆石坝变形趋势的最佳变形参数集；

s8：对变形参数集进行检验，即将反演参数代入有限元模型计算，将计算结果与原型观测数据对比，从而评价反分析参数是否合理。

传统反演方法中，在目标函数的构建上，多选用最小二乘法、极大似然法、贝叶斯法等，然而这些方法普遍存在对实测数据敏感，反演得到的参数容易偏离真实值，或者需要大量的有限元正分析计算，工作量大的问题。本技术方案在步骤s5中，采用了响应面法，来构建了响应面函数，能够经过少量的有限元计算，较准确地把特征量和待反演参数间的复杂关系用较简单的表达式表达出来，并且响应面函数的物理意义明确，设计灵活，能够避免了有限元正分析计算工作量大的缺点，计算精度较高；同时，经过标准化处理，避免了由于坝体不同位置的变形方向不同而引起变形抵消，或者变形量级不同而导致变形较小的测点对目标函数的贡献小等问题。

而在目标函数的优化上，传统反演方法多选用单纯形法和复合形法等直接搜索法以及人工神经网络训练法、粒子群优化算法等智能化反分析技术，然而这些方法普遍存在计算工作量大、耗时长、收敛速度慢等弊端，以及容易陷入局部极小值，解的稳定性差等问题。本技术方案在步骤s7中，引入模拟退火算法来对上述目标函数进行迭代优化，以一种时变且最终趋于零的概率转移到下一个新解，保证了能够搜索到空间内的每一个点，从而使得算法具有全局最优。

进一步地，所述步骤s1与步骤s2之间还包括：

s9：根据心墙堆石坝的变形特征确定反演时段。

其中，所述反演时段包括第一反演时段和第二反演时段，其中第一反演时段选取蓄水间歇期作为蠕变变形的反演时段，该时段外部荷载不变，符合蠕变变形的受力特点；第二反演时段将第一反演时段中得到的反演蠕变参数作为已知参数，选取整个填筑期和运行期作为瞬时变形的反演时段。

其中，所述第一反演时段中采用修正沈珠江模型；所述第二反演时段中采用邓肯-张模型。

修正沈珠江模型，如公式(1-1)所示。

ε(t)＝ε0+εf(1-e^-αt)(1-1)

其中ε0、εf分别为瞬时变形和最终蠕变变形，最终蠕变变形分解为最终体积蠕变量和最终剪切蠕变量，如公式(1-2)、公式(1-3)所示。

其中σ3d为体积蠕变量随应力水平增加时减小或增加的界限值，该蠕变模型的试验参数为α,b,c,d,m1,m2,m3，总的蠕变量用prandtl-reuss假设，总的应变率张量可表示为公式(1-4)。

邓肯-张模型，如公式(2-1)所示。

dσij＝3λdεmδij+2gtdεij(2-1)

公式(2-1)中λ,gt为拉梅常数和剪切模量，分别按公式(2-2)、公式(2-3)计算。

λ＝3bt(3bt-et)/(9bt-et)(2-2)

gt＝et/2(1+νt)(2-3)

公式(2-2)、公式(2-3)中et,bt为切线模量和体积模量，按公式(2-4)、公式(2-5)计算。

et＝kepa(σ3/pa)ⁿ[1-rfsl]²(2-4)

bt＝kbpa(σ3/pa)^m(2-5)

公式(2-4)中sl为应力水平，按公式(2-6)计算。

其中，是摩擦角，考虑了围压对摩擦角的修正，按公式(2-7)计算。

进一步地，所述步骤s3中，待反演参数样本的取值范围在基准值的基础上，上下浮动20％。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

1.本发明在反演过程前，选取了合适的反演时段，第一反演时段选取蓄水间歇期作为蠕变变形的反演时段，该时段外部荷载不变，符合蠕变变形的受力特点；第二反演时段将第一反演时段中得到的反演蠕变参数作为已知参数，选取整个填筑期和运行期作为瞬时变形的反演时段。

2.本发明在合适的反演时段的基础上，采用了响应面法构建了目标优化函数，能够经过少量的有限元计算，较准确地把特征量和待反演参数间的复杂关系用较简单的表达式表达出来，为坝体提供了更加准确的物理力学参数；同时，目标优化函数经过标准化处理，避免了由于坝体不同位置的变形方向不同而引起变形抵消，或者变形量级不同而导致变形较小的测点对目标函数的贡献小等问题。

3.本发明还引入模拟退火算法来对目标优化函数进行迭代优化，以一种时变且最终趋于零的概率转移到下一个新解，保证了能够搜索到空间内的每一个点，从而使得算法具有全局最优，对坝体的蠕变变形参数和瞬时变形参数进行了高效、准确的反演分析。

4.综合来说，本发明一种高心墙堆石坝变形参数反演方法，选取合适的反演时段，采用了响应面法构建了目标优化函数，引入模拟退火算法来对目标优化函数进行迭代优化，能够用于高心墙堆石坝变形参数反演分析模拟，为更加准确的进行数值模拟提供技术支撑。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1是某一高心墙堆石坝的剖面图；

图2是某一高心墙堆石坝的三维有限元模型；

图3是本发明所述反演技术的流程图；

图4是用反演参数和室内试验参数计算的位移与原型观测位移的对比图；

图5是模拟退火算法流程图；

图中，1-汛期限制水位、2-正常蓄水位、3-死水位、4-堆石料、5-上游堆石区、6-混凝土防渗墙、7-基岩、8-覆盖层、9-弃渣压重区、10-下游主堆石区、11-下游次堆石区、12-过渡区、13-砾石土心墙、14-反滤层、15-位移计v1-v13。

具体实施方式

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

需要说明的是，术语“第一”和“第二”等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

下面结合图1至图5对本发明作详细说明。

实施例1

一种高心墙堆石坝变形参数反演方法，包括以下步骤：

s1：结合心墙堆石坝的地质剖面图和变形监测布置图，建立心墙堆石坝三维有限元模型，并根据沉降仪埋设位置，在有限元模型中相应位置布置节点；

s2：进行参数敏感性分析，即从若干材料力学参数中选取对变形敏感的参数，即得到待反演参数；

s3：运用均匀设计方法设计待反演参数样本；

s4：将步骤s3中得到的待反演参数样本输入步骤s1中建立的三维有限元模型中，进行有限元计算，得到变形值样本；

s5：将步骤s4中得到的变形值样本进行标准化处理，得到标准化变形值，建立标准化变形值与待反演参数样本之间不含交叉项的响应面函数，并根据最小二乘法原理拟合响应面函数的系数；

s6：结合心墙堆石坝原型观测数据，利用响应面函数，选取多个不同测点，建立目标优化函数；

s7：引入模拟退火算法对步骤s6中建立的目标优化函数进行迭代优化，得到满足心墙堆石坝变形趋势的最佳变形参数集；

s8：对变形参数集进行检验，即将反演参数代入有限元模型计算，将计算结果与原型观测数据对比，从而评价反分析参数是否合理。

在步骤s5中，采用了响应面法，来构建了响应面函数，能够经过少量的有限元计算，较准确地把特征量和待反演参数间的复杂关系用较简单的表达式表达出来，并且响应面函数的物理意义明确，设计灵活，能够避免了有限元正分析计算工作量大的缺点，计算精度较高；同时，经过标准化处理，避免了由于坝体不同位置的变形方向不同而引起变形抵消，或者变形量级不同而导致变形较小的测点对目标函数的贡献小等问题。

在步骤s7中，引入模拟退火算法来对上述目标函数进行迭代优化，以一种时变且最终趋于零的概率转移到下一个新解，保证了能够搜索到空间内的每一个点，从而使得算法具有全局最优。

实施例2

一种高心墙堆石坝变形参数反演方法，包括以下步骤：

s1：结合心墙堆石坝的地质剖面图和变形监测布置图，建立心墙堆石坝三维有限元模型，并根据沉降仪埋设位置，在有限元模型中相应位置布置节点；

s9：根据心墙堆石坝的变形特征确定反演时段；

s2：进行参数敏感性分析，即从若干材料力学参数中选取对变形敏感的参数，即得到待反演参数；

s3：运用均匀设计方法设计待反演参数样本；

s4：将步骤s3中得到的待反演参数样本输入步骤s1中建立的三维有限元模型中，进行有限元计算，得到变形值样本；

s5：将步骤s4中得到的变形值样本进行标准化处理，得到标准化变形值，建立标准化变形值与待反演参数样本之间不含交叉项的响应面函数，并根据最小二乘法原理拟合响应面函数的系数；

s6：结合心墙堆石坝原型观测数据，利用响应面函数，选取多个不同测点，建立目标优化函数；

s7：引入模拟退火算法对步骤s6中建立的目标优化函数进行迭代优化，得到满足心墙堆石坝变形趋势的最佳变形参数集；

s8：对变形参数集进行检验，即将反演参数代入有限元模型计算，将计算结果与原型观测数据对比，从而评价反分析参数是否合理。

在步骤s9中，选取了合适的反演时段，所述反演时段包括第一反演时段和第二反演时段，其中第一反演时段选取蓄水间歇期作为蠕变变形的反演时段，该时段外部荷载不便，符合蠕变变形的受力特点；第二反演时段将第一反演时段中得到的反演蠕变参数作为已知参数，选取整个填筑期和运行期作为瞬时变形的反演时段。

实施例3

本实施例为实施例2的补充说明。

所述第一反演时段中采用修正沈珠江模型。

修正沈珠江模型，如公式(1-1)所示。

ε(t)＝ε0+εf(1-e^-αt)(1-1)

其中ε0、εf分别为瞬时变形和最终蠕变变形，最终蠕变变形分解为最终体积蠕变量和最终剪切蠕变量，如公式(1-2)、公式(1-3)所示。

其中σ3d为体积蠕变量随应力水平增加时减小或增加的界限值，该蠕变模型的试验参数为α,b,c,d,m1,m2,m3，总的蠕变量用prandtl-reuss假设，总的应变率张量可表示为公式(1-4)。

实施例4

本实施例为实施例2的补充说明。

所述第二反演时段中采用邓肯-张模型。

邓肯-张模型，如公式(2-1)所示。

dσij＝3λdεmδij+2gtdεij(2-1)

公式(2-1)中λ,gt为拉梅常数和剪切模量，分别按公式(2-2)、公式(2-3)计算。

λ＝3bt(3bt-et)/(9bt-et)(2-2)

gt＝et/2(1+νt)(2-3)

公式(2-2)、公式(2-3)中et,bt为切线模量和体积模量，按公式(2-4)、公式(2-5)计算。

et＝kepa(σ3/pa)ⁿ[1-rfsl]²(2-4)

bt＝kbpa(σ3/pa)^m(2-5)

公式(2-4)中sl为应力水平，按公式(2-6)计算。

其中，是摩擦角，考虑了围压对摩擦角的修正，按公式(2-7)计算。

实施例5

结合具体案例说明本发明一种高心墙堆石坝变形参数反演分析方法，包括以下步骤：

s1：某一高心墙堆石坝的剖面图和心墙内部变形监测点如图1所示，坝顶高程856m，三维有限元模型如图2所示，根据沉降仪埋设位置，在有限元模型中相应位置布置了节点。

s2：该高心墙堆石坝施工填筑过程基本连续，最大的施工间歇期为1个月，发生在2008年8月6日至2008年9月6日；蓄水主要分2步，首先蓄水至死水位高程790.00m，稳定近5个月后，蓄水至正常蓄水位850.00m，此后，大坝在死水位和正常蓄水位周期循环作用下运行。因此，选取稳定水位期2009年12月20日至2010年5月20日作为流变蠕变参数反分析时段，将得到的流变蠕变变形参数用于全过程的变形分析，在此过程中对填筑料的邓肯-张参数进行反分析，反演流程如图3所示：

s3：心墙堆石坝因材料分区的影响，只反演对变形影响较大的材料分区，即主堆石区、次堆石区、心墙料和覆盖层。通过敏感性分析得到对变形影响较大的蠕变模型参数有a,b,c,d，覆盖层参数与主堆石一致，且b＝c，因此待反演的蠕变模型参数共有9个；同理根据敏感性分析对变形敏感的邓肯-张模型参数，各参数的敏感顺序依次为ke,n,kb,m,rf，其中心墙和覆盖层的参数对变形不敏感，因此待反演的邓肯-张模型参数总共有22个。

s4：由于该工程未作蠕变参数的室内试验，因此参考类似工程的蠕变参数作为基准值；邓肯-张模型参数反演以表1中室内试验参数作为样本基准值。用于正分析的参数样本利用正交设计上下浮动范围为基准值的20％。

表1

s5：将上述中的每组用于正分析的参数样本分别代入到有限元中进行计算，得到位移δ。

s6：再将用于正分析的参数样本和对应的有限元计算的位移δ带入蠕变变形参数的响应面方程和瞬时变形参数的响应面方程中，如公式(3-1)和公式(3-2)所示。

并且采用最小二乘法原理拟合响应面函数，得到a19＝[0.013,0.016,...,0.014,0.019]^t，a22＝[0.312,2.210,...,-37.016,-3.198]^t，其余bki的结果、cki的结果如表2至表5所示。

表2

表3

表4

表5

s7：利用步骤s6中的响应面方程，选取多个不同测点，构建目标优化函数，本模拟中蠕变模型参数选取了22个测点，邓肯-张模型参数反演选取了20个测点，构建的目标优化函数如公式(3-3)、公式(3-4)所示。

s8：引入模拟退火算法，对步骤s7中的响应面方程进行优化，最终得到的最优参数，就是反演得到的高心墙堆石坝的变形参数，结果如表6，以及表1中反演分析部分所示。

表6

s9：将步骤s8中反演得到的参数带入有限元计算，得到的心墙内部的沉降与实测值进行对比，结果如图4所示。从与实测位移值的对比来看，本发明提出的变形参数反演方法能够高效、准确地反演得到高心墙堆石坝的物理力学参数。

以上所述，仅为本发明的优选实施方式，但本发明的保护范围并不局限于此，任何熟悉本领域的技术人员在本发明所揭露的技术范围内，可不经过创造性劳动想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书所限定的保护范围为准。