基于地震记录和SSI方法的高土石坝模态识别

王茂华; 迟世春; 周雄雄

doi:10.11779/CJGE202107013

基于地震记录和SSI方法的高土石坝模态识别 English Version

王茂华^{1, 2,},
迟世春^{1, 2, ,},
周雄雄³

1.
大连理工大学海岸与近海工程国家重点实验室,辽宁　大连　116024
2.
大连理工大学建设工程学部水利工程学院工程抗震研究所,辽宁　大连　116024
3.
西北农林科技大学水利与建筑工程学院,陕西　杨凌　712100

基金项目:

国家重点研发计划项目 2016YFB0201001

详细信息

作者简介:
王茂华(1994— ),男,博士研究生,主要从事高土石坝数值分析方面的研究工作。E-mail：wmhgryx1@163.com

通讯作者:
迟世春, E-mail：schchi@dlut.edu.cn

中图分类号: TU435
计量
- 文章访问数: 338
- HTML全文浏览量: 29
- PDF下载量: 104
出版历程
- 收稿日期: 2020-06-01
- 网络出版日期: 2022-12-02
- 刊出日期: 2021-06-30

Modal identification of high earth-rock dams based on seismic records and SSI method

1.
State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China
2.
Institute of Earthquake Engineering, School of Hydraulic Engineering, Faculty of Infrastructure Engineering, Dalian University of Technology, Dalian 116024, China
3.
College of Water Resources and Architectural Engineering, Northwest A & F University, Yangling 712100, China

摘要

摘要: 采用实测地震记录来识别高土石坝的模态参数是进行大坝动力特性研究的一种可行的方式。根据地震中坝体监测点记录的数据构建协方差驱动的随机子空间模型(stochastic subspace identification, SSI),通过构造多个不同维度的Hankel矩阵来形成稳定图,可以很好地剔除虚假模态,形成较清晰的稳定轴,采用谱系聚类的筛选方法,不仅可以有效地识别坝体的固有频率和阻尼比,而且实现了物理模态参数的自动识别,规避了人为选择引入的误差。通过两个数值算例验证了提出的随机子空间法应用于地震数据分析的有效性和精度,并将其应用于糯扎渡心墙堆石坝的模态参数识别,获得了合理的坝体固有频率和阻尼比,表明该方法可以较好地识别出坝体的模态参数,具有很好的工程应用,可以为高土石坝抗震安全有限元分析提供一定的依据和参考。
- 土石坝 /
- 模态参数 /
- 随机子空间 /
- Hankel矩阵 /
- 稳定图 /
- 谱系聚类
Abstract: Using the measured seismic records to identify the modal parameters of high earth-rock dams is a feasible way to study the dynamic characteristics of dams. Constructing the covariance-driven stochastic subspace identification (SSI) based on the data recorded at the monitoring points of a dam during the earthquake and the stable graph by establishing multiple Hankel matrices with different dimensions can effectively eliminate false modes. The hierarchical clustering can effectively identify the natural frequency and damping ratio of the dam, and realize the automatic identification of physical modal parameters, avoiding the errors introduced by human selection. The effectiveness and accuracy of the proposed SSI method for seismic data analysis are verified by numerical examples. It is applied to the identification of modal parameters of Nuozhadu core-wall rockfill dam, and the reasonable frequency and damping ratio are obtained, indicating that the method can identify the modal parameters of the dam well and has a good engineering application, and it can provide a certain basis and reference for the finite element analysis of seismic safety of high earth-rock dams.
- earth-rock dam /
- modal parameter /
- stochastic subspace identification /
- Hankel matrix /
- stabilization diagram /
- hierarchical clustering

HTML全文

0. 引言

在基坑工程中，地下水回灌常被用于控制或消除基坑降水对周围环境的影响。地下水回灌可以分为定流量回灌和定水头回灌，其中定水头回灌常用于低渗透性含水层的水位抬升^[1-3]。最早的定水头井流模型是由Jacob等^[4]建立的。该模型描述了水平无限延伸的均质承压含水层中完整井定降深条件下的径向水流运动。在此基础上，国内外众多学者考虑含水层越流、表皮效应、非完整井、含水层类型等因素的影响，对Jacob等模型进行了改进并推导了相应的解析解或半解析解^[1-3]。由于地下水回灌通常被视为地下水开采的逆过程，上述单井定降深井流模型可同样为定水头回灌条件下的含水层渗流响应分析提供参考。然而，由于相反的水流流动方向改变了井周的物理条件，不同于抽汲地下水在井周形成的渗透增高地带，回灌地下水携带的细颗粒会随着水流速度下降逐渐沉积而在井周形成一个渗透衰减地带^[5]。与此同时，回灌水中溶解的气泡、回灌水与地下水或含水层介质间的化学反应以及滤管周围微生物的生长繁殖均可能进一步加剧井周含水层渗透性衰减。上述物理、生物与化学过程导致的渗透衰减使得基于恒定渗透系数的抽水井流模型不再适用于分析回灌井流问题。

为考虑回灌过程中渗透衰减效应，何满潮等^[6]通过对地热水回灌的测试数据分析，建立了考虑井周物理与化学阻塞效应的渗透系数衰减模型。随后，何满潮等^[5]将该服从指数衰减的渗透模型引入经典的Theis模型^[7]，推导了考虑渗透系数动态衰减的单井回灌公式，并分析了渗透衰减对于回灌井周渗流场的影响。基于上述单井回灌公式，闫峭等^[8]利用改进遗传算法研究了考虑渗透系数衰减的水文地质参数反演问题。上述研究均只考虑了定流量回灌，并未涉及定水头回灌。对于由回灌压力控制的回灌井流问题，李炯等^[9]将何满潮等^[6]提出的渗透衰减方程引入至Jacob等^[4]模型，并利用Weber变换得到了模型的解析解，但该模型未考虑回灌井非完整性的影响。此外，上述研究均假设整个含水层的渗透系数存在同步衰减，而现实中回灌堵塞通常仅存在于回灌井邻近的柱状包围带^[5]，堵塞对主含水层渗透性的影响非常轻微。

为分析定水头非完整井回灌过程中井周堵塞效应对承压含水层渗流场的影响，本文选用指数形式的渗透方程代替恒定渗透系数反映堵塞作用下的井周含水层渗透衰减行为，建立了考虑井周堵塞效应的非完整井定水头回灌井流模型。通过变量代换、Laplace变换和有限余弦Fourier变换及其逆变换，推导了井流模型在Laplace空间的解，并利用Stehfest数值逆变换方法得到了实时空间内的半解析解，并进一步分析了回灌堵塞诱发的渗透衰减对定水头非完整井回灌渗流场的影响规律。研究结果可以为定水头回灌的工程设计、参数反演、堵塞治理提供理论参考。

1. 问题的提出与求解

1.1 渗透衰减模型

由于物理、生物以及化学因素引起的堵塞在回灌井周围并不是均匀分布的，距离回灌井越近，地层堵塞越严重，渗透系数下降越显著，使得回灌堵塞影响下的井周含水层渗透系数兼有时间和空间变异性。然而，当前试验条件下准确获得径向流动过程中的堵塞时空发展规律及其导致的渗透时空变异规律是非常困难的。同时，考虑时空变异的渗透衰减方程对于井流模型的解析求解亦提出了更高的要求。为便于求解与分析，本文仅考虑渗透系数的时变效应。需要指出的是，这里的渗透系数表征的是整个井周堵塞区而非某一给定位置的水力传导性能，是对堵塞区实际渗透系数在空间上的平均唯象近似。

基于上述分析，为考虑回灌堵塞对含水层导水特性的影响，本文假设井周含水层在水平和竖直方向上的渗透系数均服从Bianchi等^[10]基于室内水平砂柱试验建立的反映堵塞作用下渗透衰减问题的概括性方程。该方程具有如下形式：

${K_{\text{r}}}(t) = {K_{{\text{r}},\infty }} + ({K_{{\text{r}},0}} - {K_{{\text{r}},\infty }})\exp ( - {\lambda _0} \cdot t) \;\;\text{，}$

(1a)

${K_{\text{z}}}(t) = {K_{{\text{z}},\infty }} + ({K_{{\text{z}},0}} - {K_{{\text{z}},\infty }})\exp ( - {\lambda _0} \cdot t) \;\;。$

(1b)

式中：K_r和K_z分别为含水层水平和竖直方向上的渗透系数；“0”和“∞”分别表示渗透系数的初始与渐近值； ${\lambda _0}$ 为渗透系数随时间的衰减指数。当K_r = K_z且K_∞ = 0时，上述渗透衰减模型可退化至何满潮等^[6]提出的渗透衰减模型。此外，为便于数学分析，本文假设含水层水平和竖向渗透系数具有相同的衰减速率，即相同的 ${\lambda _0}$ ，且两者之比的初始值和稳定值保持不变，即K_{r, 0}/K_{z, 0} = K_{r, ∞}/K_{z, ∞} = ${\alpha ^2}$ 。此时，在整个渗透衰减过程中K_r(t)与K_z(t)的比值亦为定值 ${\alpha ^2}$ 。

1.2 数学模型

本文选用的非完整井定水头回灌模型位于图 1所示的含水层系统中。该含水层系统由一个均质等厚、侧向无限延伸的水平各向同性承压含水层及其上覆和下伏隔水层组成，是一种工程实践中常见的含水层系统，便于开展相关的野外现场测试和室内模型试验，从而有利于验证理论模型的适用性和准确性。同时，该含水层系统也被众多经典井流模型选用，如Yang等^[3]建立的考虑表皮效应的两区井流模型。鉴于此，本文采用了与Yang等模型几乎完全相同的含水层假设。但不同于Yang等假定井周表皮层渗透系数随时间恒定，本文模型须考虑回灌堵塞造成的井周含水层渗透衰减效应。对于井周含水层回灌堵塞区，参考Yang等^[3]对表皮层的处理方式，本文假设井周堵塞区的径向延伸r_c沿深度保持不变，呈圆柱状分布，如图 1所示。在该柱状包围带内，含水层具有良好的均质性，不同位置的渗透系数受堵塞效应影响而随时间同步衰减。

图 1 承压含水层中非完整井定水头回灌示意

Figure 1. Schematic graph of constant-head injection in confined aquifers using a partially penetrating well

下载: 全尺寸图片幻灯片

基于上述假设，以回灌井轴线与含水层底板交点为原点、回灌井轴线竖直向上为z轴、含水层底板水平径向延伸为r轴建立柱坐标系。分别用式（1a）和式（1b）代替Yang等^[3]模型中水平与竖直方向上的常渗透系数，则考虑堵塞效应的非完整井定水头回灌模型的控制方程可以描述为

${K_{\text{r}}}(t)\left[ {\frac{{{\partial ^2}s}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial s}}{{\partial r}}} \right] + {K_{\text{z}}}(t)\frac{{{\partial ^2}s}}{{\partial {z^2}}} = \frac{S}{B}\frac{{\partial s}}{{\partial t}}\;\;。$

(2)

式中：s为含水层的水头抬升；r和z分别为到回灌井轴线的径向距离和到含水层底板的垂向距离；t为回灌时间；S为含水层储水系数；B为含水层厚度。

回灌开始前，井周含水层初始水头变化为零：

$s(r,z,t) = 0,{\text{ }}t = 0\;\;。$

(3)

回灌过程中，井筒内的平均压力水头恒为s_w：

$\frac{1}{{l - d}}{\left. {\int_d^l {s(r,z,t)} {\text{d}}z} \right|_{r = {r_{\text{w}}}}} = {s_{\text{w}}}{\text{ (}}d \leqslant z \leqslant l) \;\;。$

(4)

式中：r_w为回灌井半径；l和d分别为回灌井过滤器顶端与底端的竖向坐标。

在回灌井内，假设定水头作用下回灌流量Q(t)沿滤管均匀分布，则有

${\left. { - 2{\text{π }}r(l - d){K_{\text{r}}}(t){{\partial s(r,z,t)} \mathord{\left/ {\vphantom {{\partial s(r,z,t)} {\partial r}}} \right. } {\partial r}}} \right|_{r = {r_{\text{w}}}}}{\text{ = }}Q(t){\text{ (}}d \leqslant z \leqslant l) \;\;。$

(5)

回灌过程中，沿井壁套筒流量为零，含水层与下伏和上覆弱透水层的交界面完全隔水，则有

${\left. {{{\partial s(r,z,t)} \mathord{\left/ {\vphantom {{\partial s(r,z,t)} {\partial r}}} \right. } {\partial r}}} \right|_{r = {r_{\text{w}}}}} = {\text{0 (}}z \leqslant d,{\text{ }}z \geqslant l) \;\;\text{，}$

(6)

${{\partial s(r,z,t)} \mathord{\left/ {\vphantom {{\partial s(r,z,t)} {\partial z}}} \right. } {\partial z}} = 0\;\;\;{\text{ (}}z = 0)\;\;\text{，}$

(7)

${{\partial s(r,z,t)} \mathord{\left/ {\vphantom {{\partial s(r,z,t)} {\partial z}}} \right. } {\partial z}} = 0\;\;\;{\text{ (}}z = B)\;\;。$

(8)

回灌过程中，井周堵塞区与主含水层的交界面上的水头抬升如下式所示：

$s\left(r_{\mathrm{c}}, z, t\right)=s_{\mathrm{c}}(z, t) \quad(0 \leqslant z \leqslant B) \;\;。$

(9)

式中，s_c(z, t)为回灌流量等于Q(t)且不考虑井周堵塞效应时求得的r = r_c处的水头抬升，其表达式与Feng等^[1]推导的变流量条件下承压含水层中水头变化公式相同。式（9）实际上忽略了井周堵塞对主含水层的渗流场的影响。根据Yeh等^[11]对考虑表皮效应的两区井流模型的分析结果可知，试井周围的低渗透区域仅会在试验初期对主含水层渗流场有轻微影响。随着井周渗流场在很短时间内达到拟稳态，该影响仅会在试验初期带来非常微小的误差，是可以接受的。

1.3 问题的求解

式（2）~（9）构成了考虑井周堵塞效应的非完整井定水头回灌井流模型的定解问题。为便于求解，首先对定解问题作如下变量代换：

$r^{\prime} = {r \mathord{\left/ {\vphantom {r \alpha }} \right. } \alpha } \;\;\text{，}$

(10)

$t^{\prime} = {K_{{\text{z}},\infty }}t/{K_{{\text{z}},0}} + ({K_{{\text{z}},0}} - {K_{{\text{z}},\infty }})(1 - {{\text{e}}^{ - {\lambda _0}t}})/({\lambda _0}{K_{{\text{z}},0}}) \;\;\text{，}$

(11)

式中， $\alpha = \sqrt {{{{K_{{\text{r,0}}}}} \mathord{\left/ {\vphantom {{{K_{{\text{r,0}}}}} {{K_{{\text{z,0}}}}}}} \right. } {{K_{{\text{z,0}}}}}}} = \sqrt {{{{K_{{\text{r,}}\infty }}} \mathord{\left/ {\vphantom {{{K_{{\text{r,}}\infty }}} {{K_{{\text{z,}}\infty }}}}} \right. } {{K_{{\text{z,}}\infty }}}}}$ 。经过坐标变换之后，含水层的水头抬升s(r, z, t)用s^'(r^', z, t^')表示，则定解问题转换为

$\frac{{{\partial ^2}s^{\prime}}}{{\partial {{r^{\prime}}^2}}} + \frac{1}{{r^{\prime}}}\frac{{\partial s^{\prime}}}{{\partial r^{\prime}}} + \frac{{{\partial ^2}s^{\prime}}}{{\partial {z^2}}} = \frac{S}{{B{K_{{\text{z}},0}}}}\frac{{\partial s^{\prime}}}{{\partial t^{\prime}}} \;\;\text{，}$

(12)

$s^{\prime}\left( {r^{\prime},z,t^{\prime}} \right) = 0,{\text{ }}t^{\prime} = 0\;\;\text{，}$

(13)

$\frac{1}{{l - d}}{\left. {\int_d^l {s^{\prime}(r^{\prime},z,t^{\prime})} {\text{d}}z} \right|_{r^{\prime} = {{r}_{\text{w}}^{\prime}}}} = {s_{\text{w}}}\;\;\;{\text{ (}}d \leqslant z \leqslant l) \;\;\text{，}$

(14)

${\left. {\frac{{\partial s^{\prime}(r^{\prime},z,t^{\prime})}}{{\partial r^{\prime}}}} \right|_{r^{\prime} = {{r}_{\text{w}}^{\prime}}}}{\text{ = }} - \frac{{{{Q^{\prime}}_1}(t^{\prime})}}{{2{\text{π }}{{r}_{\text{w}}^{\prime}}(l - d){K_{{\text{r,0}}}}}}\;\;\;{\text{ (}}d \leqslant z \leqslant l) \;\;\text{，}$

(15)

${\left. {{{\partial s^{\prime}(r^{\prime},z,t^{\prime})} \mathord{\left/ {\vphantom {{\partial s^{\prime}(r^{\prime},z,t^{\prime})} {\partial r^{\prime}}}} \right. } {\partial r^{\prime}}}} \right|_{r^{\prime} = {{r}_{\text{w}}^{\prime}}}} = 0\;\;\;{\text{ (}}z \leqslant d,{\text{ }}z \geqslant l) \;\;\text{，}$

(16)

${{\partial s^{\prime}(r^{\prime},z,t^{\prime})} \mathord{\left/ {\vphantom {{\partial s^{\prime}(r^{\prime},z,t^{\prime})} {\partial z}}} \right. } {\partial z}} = 0\;\;\;{\text{ (}}z = 0)\;\;\text{，}$

(17)

${{\partial s^{\prime}{\text{(}}r^{\prime},z,t^{\prime}{\text{)}}} \mathord{\left/ {\vphantom {{\partial s^{\prime}{\text{(}}r^{\prime},z,t^{\prime}{\text{)}}} {\partial z}}} \right. } {\partial z}} = 0\;\;\;{\text{ (}}z = B)\;\;\text{，}$

(18)

$s^{\prime}({r^{\prime}_{\text{c}}},z,t^{\prime}){\text{ = }}{s^{\prime}_{\text{c}}}(z,t^{\prime})\;\;\;{\text{ (}}0 \leqslant z \leqslant B) \;\;。$

(19)

其中，式（15）中 ${Q^{\prime}_1}\left( {t^{\prime}} \right){\text{ = }}{{Q^{\prime}\left( {t^{\prime}} \right){K_{{\text{r,0}}}}} \mathord{\left/ {\vphantom {{Q^{\prime}\left( {t^{\prime}} \right){K_{{\text{r,0}}}}} {{{K^{\prime}}_{\text{r}}}\left( {t^{\prime}} \right)}}} \right. } {{{K^{\prime}}_{\text{r}}}\left( {t^{\prime}} \right)}}$ ，Q^'(t^')和K_r^'(t^')以及式（19）中的 ${s^{\prime}_{\text{c}}}(z,t^{\prime})$ 分别为回灌流量Q(t)、渗透系数方程K_r(t)以及堵塞区外部水头抬升 ${s_{\text{c}}}(z,t)$ 作关于时间的变量代换之后得到的表达式。由于根据式（11）获得t^'关于t的表达式十分困难，使得上述表达式的具体形式同样难以确定。根据 ${s^{\prime}_{\text{c}}}(z,t^{\prime}) = {s_{\text{c}}}(z,t(t^{\prime}))$ ， ${s^{\prime}_{\text{c}}}(z,t^{\prime})$ 可表示为 ${s_{\text{c}}}(z,t^{\prime})$ 与 $\Delta {s_{\text{c}}}(z,t^{\prime})$ 之和，其中 $\Delta {s_{\text{c}}}(z,t^{\prime})$ 的表达式通过非线性拟合求得。考虑到后续求解需要对 $\Delta {s_{\text{c}}}(z,t^{\prime})$ 分别进行关于z和t^'的积分变换，因此 $\Delta {s_{\text{c}}}(z,t^{\prime})$ 的拟合公式宜采用如下可分离变量的形式：

$\Delta {s_{\text{c}}}(z,t^{\prime}) = \Delta {s_{\text{c}}}(z) \cdot \Delta {s_{\text{c}}}(t^{\prime}) \;\;。$

(20a)

同时，根据Laplace变换和有限Fourier变换的定义，当原函数为指数函数时，可以直接求得积分变换后的象函数。因此， $\Delta {s_{\text{c}}}(z)$ 和 $\Delta {s_{\text{c}}}(z)$ 宜采用如下形式的分段指数函数进行拟合：

$\Delta {s_{\text{c}}}(t^{\prime}) = \left[ {U(t^{\prime} - {{t^{\prime}}_{i + 1}}) - U(t^{\prime} - {{t^{\prime}}_i})} \right]\left[ {{A_i} + {A_{i + 1}}\exp ({\alpha _i}t^{\prime})} \right] \;\;\text{，}$

(20b)

$\Delta {s_{\text{c}}}(z) = \left[ {U(z - {z_{j + 1}}) - U(z - {z_j})} \right]\left[ {{C_j} + {C_{j + 1}}\exp ({\varepsilon _j}z)} \right] \;\;。$

(20c)

式中：A_i， ${\alpha _i}$ ，C_j及ε_j为根据 ${s_{\text{c}}}\left[ {z,t(t^{\prime})} \right]$ 与 ${s_{\text{c}}}(z,t^{\prime})$ 的差值序列确定的拟合参数；U(·)为Heaviside跃阶函数，当t < ${t_i}$ 时，U $(t - {t_i})$ = 0，当t ≥ ${t_i}$ 时，U $(t - {t_i})$ =1， ${t^{\prime}_i}$ 和 ${z_j}$ 为分段点。

对定解问题作关于时间 $t^{\prime}$ 的Laplace变换和关于z的有限Fourier余弦变换可得

$\frac{\partial^2 \hat{\bar{s}^{\prime}}}{\partial r^{\prime 2}}+\frac{1}{r^{\prime}} \frac{\partial \hat{\bar{s}^{\prime}}}{\partial r^{\prime}}-\xi(n, p)^2 \hat{\bar{s}^{\prime}}=0\;\;\;,$

(21)

$\frac{1}{{l - d}}\int_d^l {\bar s^{\prime}({{r}_{\text{w}}^{\prime}},z,p)} {\text{d}}z = \frac{{{s_{\text{w}}}}}{p}{\text{ (}}d \leqslant z \leqslant l)\;\;\text{，}$

(22)

$\left.K_{\mathrm{r}, 0} \frac{\partial \hat{\bar{s}^{\prime}}\left(r^{\prime}, n, p\right)}{\partial r^{\prime}}\right|_{r^{\prime}=r_{\mathrm{w}}^{\prime}}=\frac{-\bar{Q}_1^{\prime}(p)}{2 \pi r_{\mathrm{w}}^{\prime}(l-d)} \cdot A(n)\;\;\text{，}$

(23)

$\hat{\bar{s}^{\prime}}({r^{\prime}_{\text{c}}},n,p){\text{ = }}{\hat{\bar{s}^{\prime}}_{\text{c}}}(n,p) + \Delta {\hat{\bar{s}^{\prime}}_{\text{c}}}(n,p) \;\;。$

(24)

式中： $\hat{\bar{s}^{\prime}}(r^{\prime},n,p)$ 和 $\Delta {\hat{\bar{s}^{\prime}}_{\text{c}}}(n,p)$ 分别为水头抬升 $s(r^{\prime},n,p)$ 和水头差 $\Delta {s_{\text{c}}}(n,p)$ 在Laplace-Fourier空间的表达式，符号“-”和“^”分别表示Laplace和Hankel空间的项，p和n分别为Laplace变量和Fourier变量， $\xi {({r}^{\prime },n,p)}^{2}={w}_{n}^{2}+Sp/(B{K}_{\text{z},0})，A(n)=[\mathrm{sin}({w}_{n}l)-$ $\sin ({w_n}d)]/{w_n}$ ， ${w_n} = n{\text{π /}}B$ 。式（21）为零阶修正Bessel方程，通解为

$\begin{array}{l} \hat{\bar{s}}(r^{\prime},n,p) = {C_1}(n,p){{\text{I}}_0}\left[ {{\xi _0}(n,p)r^{\prime}} \right] + {C_2}(n,p){{\text{K}}_0} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {{\xi _0}(n,p)r^{\prime}} \right] \;\;。 \end{array}$

(25a)

式中：待定系数C₁(n, p)和C₂(n, p)根据边界条件式（23）和（24）确定表达式如下所示：

${C}_{1}(n,p)=\frac{\widehat{\overline{{s}}^{\prime}}({{r}}_{\text{c}}^{\prime },n,p){\text{K}}_{1}\left[{\xi }_{0}(n,p){{r}}_{\text{w}}^{\prime }\right]-\frac{{\overline{{Q}}}_{1}^{\prime}(p)A(n){\text{K}}_{0}\left[{\xi }_{0}(n,p){{r}}_{\text{c}}^{\prime }\right]}{2\text{π}{{r}}_{\text{w}}^{\prime }(l-d){K}_{\text{r,0}}{\xi }_{0}(n,p)}}{{\text{I}}_{\text{1}}\left[{\xi }_{0}(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{K}}_{0}\left[{\xi }_{0}(n,p){{r}}_{\text{c}}^{\prime }\right]+{\text{K}}_{1}\left[{\xi }_{0}(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{I}}_{0}\left[{\xi }_{0}(n,p){{r}^{\prime }}_{}\right]}\text{，}$

(25b)

${C}_{2}(n,p)=\frac{\widehat{\overline{{s}}^{\prime}}({{r}}_{\text{c}}^{\prime },n,p){\text{I}}_{1}\left[{\xi }_{0}(n,p){{r}}_{\text{w}}^{\prime }\right]+\frac{{\overline{{Q}}}_{1}^{\prime}(p)A(n){\text{I}}_{0}\left[{\xi }_{0}(n,p){{r}}_{\text{c}}^{\prime }\right]}{2\text{π}{{r}}_{\text{w}}^{\prime }(l-d){K}_{\text{r,0}}{\xi }_{0}(n,p)}}{{\text{I}}_{1}\left[{\xi }_{0}(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{K}}_{0}\left[{\xi }_{0}(n,p){{r}}_{\text{c}}^{\prime }\right]+{\text{K}}_{1}\left[{\xi }_{0}(n,p){{r}}_{\text{w}}^{\prime }\right]{I}_{0}\left[{\xi }_{0}(n,p){{r}}_{\text{c}}^{\prime }\right]}。$

(25c)

式中：I_j(·) (j = 0, 1)和K_j(·) (j = 0, 1)分别为第一类和第二类j阶虚宗量Bessel函数。对式（25）作有限Fourier余弦逆变换，则定解问题在Laplace空间的解可写作

$\bar s^{\prime}(r,z,p) = \frac{1}{B}\hat{\bar{s}^{\prime}}(r,0,p) + \frac{2}{B}\sum\limits_{n = 1}^\infty {\hat{\bar{s}^{\prime}}(r,n,p)} \cdot \cos ({w_n}z)\;\;。$

(26)

式中， ${\bar Q^{\prime}_1}(p)$ 的表达式仍然未知，可通过将式（25）和式（26）代入边界条件式（22）确定其表达式。需要注意的是，在计算 ${\bar Q^{\prime}_1}(p)$ 时还需首先确定受堵塞区外部边界的水头抬升 ${s_{\text{c}}}(z,t)$ ， ${s_{\text{c}}}(z,t)$ 又取决于通过过滤器的回灌流量Q(t)，而Q(t)的表达式仍然未知。为此，可以首先指定Q(t)一个初始值并计算该初始值条件下的 ${s_{\text{c}}}(z,t)$ ，进一步根据式（26）和 ${\bar Q^{\prime}_1}(p)$ 的表达式对Q(t)进行更新，不断循环迭代即可最终确定Q(t)的近似表达式。在迭代过程中，Q(t)可用下式近似：

$Q(t){\text{ = }}{Q_\infty }(t) + \left[ {{Q_0}(t) - {Q_\infty }(t)} \right] \cdot \sum {\left[ {{B_i}\exp ( - {\lambda _i}t)} \right]} 。$

(27)

式中，Q₀(t)和Q_∞(t)分别为井周含水层渗透系数为K_{r, 0}和K_{r, ∞}时的回灌流量，其表达式与Yang等^[3]推导的承压含水层中非完整定降深两区模型的抽水流量解相同，B_i和 ${\lambda _i}$ 为待定拟合参数，式（27）在i = 1且B_i = 1， ${\lambda _1}$ = ${\lambda _0}$ 时的计算结果可作为计算Q(t)的迭代初始值。

1.4 回灌初期与后期的近似解

在回灌伊始堵塞效应尚不明显，井周含水层与主含水层的导水特性无明显差异。此时，描述定水头回灌条件下井周含水层渗流场的定解问题与同等条件下主含水层渗流场的定解问题完全相同。将定解问题（2）~（9）中的K_r(t)和K_z(t)分别用K_{r, 0}和K_{z, 0}替换，则式（11）可以简化为 $t^{\prime} = t$ 。同时，令r_c → ∞，无穷远处的水头抬升可以忽略不计，即s_c(z, t) = 0。保持其他条件不变，得到新的定解问题。重复前述变量代换以及积分变换过程，得到新定解问题在Laplace空间的解为

$\begin{array}{l} \bar s(r^{\prime},z,p) = \frac{1}{B}\frac{{\bar q^{\prime}(p)A(0){{\text{K}}_0}\left[ {{\xi _0}(0,p)r^{\prime}} \right]}}{{{\xi _0}(0,p){K_{{\text{r,0}}}}{{\text{K}}_1}\left[ {{\xi _0}(0,p){{r}_{\text{w}}^{\prime}}} \right]}} + \\ \;\;\;\;\;\;\;\;\;\;\frac{2}{B}\sum\limits_{n = 1}^\infty {\frac{{\bar q^{\prime}(p){{\text{K}}_0}\left[ {{\xi _0}(n,p)r^{\prime}} \right]A(n)}}{{{\xi _0}(n,p){K_{{\text{r,0}}}}{{\text{K}}_1}\left[ {{\xi _0}(n,p){{r}_{\text{w}}^{\prime}}} \right]}}} \cos ({w_n}z) \;\;。 \end{array}$

(28a)

式中， $\bar q^{\prime}(p) = {{\bar Q^{\prime}(p)} \mathord{\left/ {\vphantom {{\bar Q^{\prime}(p)} {\left[ {2{\text{π }}{{r}_{\text{w}}^{\prime}}(l - d)} \right]}}} \right. } {\left[ {2{\text{π }}{{r}_{\text{w}}^{\prime}}(l - d)} \right]}}$ ，可根据边界条件式（22）确定，其表达式为

$\begin{aligned} \bar{q}(p)= & \frac{s_{\mathrm{w}}}{p}\left\{\frac{1}{B} \frac{\mathrm{K}_0\left[\xi_0(0, p) r_{\mathrm{w}}^{\prime}\right]}{\xi_0(0, p) K_{\mathrm{r}, 0} \mathrm{~K}_1\left[\xi_0(0, p) r_{\mathrm{w}}^{\prime}\right]}+\right. \\ & \left.\frac{2}{B(l-d)} \sum\limits_{n=1}^{\infty} \frac{\mathrm{K}_0\left[\xi_0(n, p) r_{\mathrm{w}}^{\prime}\right][A(n)]^2}{\xi_0(n, p) K_{\mathrm{r}, 0} \mathrm{K}_1\left[\xi_0(n, p) r_{\mathrm{w}}^{\prime}\right]}\right\}^{-1}\;\;。 \end{aligned}$

(28b)

式（28）即为非完整井定水头回灌初期水头抬升在Laplace空间的解，也是同等回灌条件下不考虑堵塞效应的承压含水层水头抬升在Laplace空间的解，且与Feng等^[1]推导的变流量条件下承压含水层中水头降深公式相同。特殊地，该解在 $r^{\prime}{\text{ = }}{r^{\prime}_{\text{c}}}$ 处的表达式即为式（9）右端项 ${s^{\prime}_{\text{c}}}(z,t^{\prime})$ 在Laplace空间的表达式。

当回灌试验进行足够长时间，井周堵塞影响下的含水层导水特性早已稳定。井周含水层水平与竖向渗透系数分别等于K_{r, ∞}和K_{z, ∞}。将式（2）~（9）中的K_r(t)和K_z(t)分别用K_{r, ∞}和K_{z, ∞}替换，得到新的定解问题。根据式（10）进行关于r的坐标变换并重复前述Laplace变换和Fourier变换，得到新定解问题在Laplace空间的解，其待定系数C₁(n, p)和C₂(n, p)表示如下：

${C}_{1}(n,p)=\frac{\widehat{\overline{{s}}^{\prime}}({{r}}_{\text{c}}^{\prime },n,p){\text{K}}_{1}\left[{\xi }_{\infty }(n,p){{r}}_{\text{w}}^{\prime }\right]-\frac{\overline{{Q}}^{\prime}(p)A(n){\text{K}}_{0}\left[{\xi }_{\infty }(n,p){{r}}_{\text{c}}^{\prime }\right]}{2\text{π}{{r}}_{\text{w}}^{\prime }(l-d){K}_{\text{r,}\infty }{\xi }_{\infty }(n,p)}}{{\text{I}}_{1}\left[{\xi }_{\infty }(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{K}}_{0}\left[{\xi }_{\infty }(n,p){{r}}_{\text{c}}^{\prime }\right]+{\text{K}}_{1}\left[{\xi }_{\infty }(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{I}}_{0}\left[{\xi }_{\infty }(n,p){{r}}_{\text{c}}^{\prime }\right]}\text{，}$

(29a)

${C}_{2}(n,p)=\frac{\widehat{\overline{{s}}^{\prime}}({{r}}_{\text{c}}^{\prime },n,p){\text{I}}_{1}\left[{\xi }_{\infty }(n,p){{r}}_{\text{w}}^{\prime }\right]+\frac{\overline{{Q}}^{\prime}(p)A\text{(}n{\text{)I}}_{0}\left[{\xi }_{\infty }(n,p){{r}}_{\text{c}}^{\prime }\right]}{2\text{π}{{r}}_{\text{w}}^{\prime }(l-d){K}_{\text{r,}\infty }{\xi }_{\infty }(n,p)}}{{\text{I}}_{1}\left[{\xi }_{\infty }(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{K}}_{0}\left[{\xi }_{\infty }(n,p){{r}}_{\text{c}}^{\prime }\right]+{\text{K}}_{1}\left[{\xi }_{\infty }(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{I}}_{0}\left[{\xi }_{\infty }(n,p){{r}}_{\text{c}}^{\prime }\right]}。$

(29b)

式（25a），（26），（29a），（29b）即为非完整井定水头回灌后期水头抬升在Laplace空间的解，其中 $\hat{\bar{s}^{\prime}}({r^{\prime}_{\text{c}}},n,p)$ 的表达式与式（28a），（28b）相同， ${\xi _\infty }{(n,p)^2} = w_n^2 +$ $Sp/(B{K_{{\text{z}},\infty }})$ 。此时， $\hat{\bar{s}^{\prime}}({r^{\prime}_{\text{c}}},n,p)$ 表达式中的 $\bar Q^{\prime}(p)$ 仍然未知，可联立式（25a）~（25c），（26），（28a），（28b），（29a），（29b）确定其表达式。

由于上述Laplace空间的解形式复杂，无法直接利用解析方法求得其在实时空间的表达式。本文在对比多种Laplace变换数值反演方法后，选用兼具计算效率和精度的Stehfest^[12]数值反演方法计算实时空间的水头抬升。需要指出的是，Stehfest算法中求和级数的项数N对反演结果具有显著影响，试算结果表明，当8≤N≤16时不同N值对应的反演结果差别不大，因此本文在进行Laplace变换的数值反演假设N = 16。

2. 结果验证与分析

目前由于缺少相关测试数据，本节首先通过数值计算方法对本文推导的半解析解进行对比验证，随后重点分析非完整井定水头回灌过程中堵塞诱发的渗透衰减对含水层渗流响应的影响规律。本节用于的验证和分析的非完整井定水头回灌井流模型如图 1所示。如无特殊声明，本节在进行模型计算时选取的试井参数和含水层参数分别如表 1，2所示。

表 1 井身结构参数

Table 1. Parameters of well structure 单位: m

试井类型	回灌压力	井径	滤管长度	滤管顶端高度	滤管底端高度
回灌井	200	0.1	10	15	5

下载: 导出CSV

| 显示表格

表 2 含水层水文地质参数

Table 2. Hydrological parameters of recharged aquifer

堵塞区厚度/m	初始渗透系数/(m·h^-1)		渐近渗透系数/(m·h^-1)		渗透衰减速率/h^-1	含水层厚度/m	储水系数
堵塞区厚度/m	水平	竖向	水平	竖向	渗透衰减速率/h^-1	含水层厚度/m	储水系数
1.0	0.06	0.006	0.03	0.003	0.05	20	0.004

下载: 导出CSV

| 显示表格

2.1 解的验证

当l = B和d = 0时，本文建立的定解问题退化为完整井定水头回灌井流模型。此时，含水层中的渗流为径向一维流动，回灌引起的水头抬升与深度无关。这时如果对Laplace域内堵塞区外部水头抬升 $\bar s^{\prime}({r^{\prime}_{\text{c}}},p)$ 作有限Fourier余弦变化，则可得到

$\hat{\bar{s}^{\prime}}\left(r_{\mathrm{c}}^{\prime}, n, p\right)=\left\{\begin{array}{ll} B \cdot \bar{s}^{\prime}\left(r_{\mathrm{c}}^{\prime}, p\right) & (n=0) \\ 0 & (n \neq 0) \end{array}\right. \;\;。$

(30)

此外，将l = B，d = 0代入A(n)化简得到

$A(n)=\left\{\begin{array}{ll} B & n=0 \\ 0 & n \neq 0 \end{array}\right. \;\;。$

(31)

将式（30），（31）代入式（25），（26）即得到完整井定水头回灌井流模型在Laplace空间的半解析解。此外，完整井定水头回灌井流模型还可通过MATLAB的偏微分方程求解函数pdepe进行数值求解。该函数可以求解如下形式的偏微分定解问题：

$c\left(x,t,u,\frac{\partial u}{\partial x}\right)\frac{\partial u}{\partial t}=\frac{1}{x}\frac{\partial }{\partial x}\left(xf\left(x,t,u,\frac{\partial u}{\partial x}\right)\right)+s\left(x,t,u,\frac{\partial u}{\partial x}\right)\text{，}$

(32a)

$u(x,0) = {u_0}(x) \;\;\text{，}$

(32b)

$p(x,t,u) + q(x,t)f(x,t,u,{{\partial u} \mathord{\left/ {\vphantom {{\partial u} {\partial x}}} \right. } {\partial x}}) = 0 \;\;。$

(32c)

其中，式（32a）~（32c）分别为定解问题的控制方程、初值条件和边值条件，x和t分别为空间变量和时间变量，c(x, t, u, ∂u/∂t)为控制方程非定常项的系数矩阵，f(x, t, u, ∂u/∂t)和s(x, t, u, ∂u/∂t)分别为控制方程的流项和源项。将完整井定水头回灌井流模型按照式（32a）~（34c）进行改写，并编写MATLAB计算程序，得到定解问题的数值解。图 2展示了不同位置的水头抬升的半解析计算值与数值计算值的对比结果。对比结果表明，利用Stehfest算法反演得到的实时空间的解在N = 16时与相应数值解符合良好，证明本文推导半解析解的准确性以及Stehfest算法在解决考虑变渗透系数的定水头非稳定流问题的适用性。

图 2 本文半解析解与数值解的比较

Figure 2. Comparison between semi-analytical and numerical solutions

下载: 全尺寸图片幻灯片

2.2 回灌堵塞程度的影响

图 3所示为不同堵塞程度，即不同K_{r, ∞}值，对井周堵塞区内不同位置（z = 10 m，r = 0.2，0.5，1.0 m）处水头抬升的影响，其中图 3（a）为水头抬升s时程曲线，图 3（b）为考虑渗透衰减（K_{r, ∞} < 0.06 m/h）与不考虑渗透衰减（K_{r, ∞} = 0.06 m/h）情形下的水头差Δs时程曲线。在回灌初期，堵塞效应尚不明显，任一给定位置处不同K_{r, ∞}值下的s随时间以相同速度逐渐增大。在该阶段不同位置处考虑堵塞效应与否的水头差Δs基本为零。随着回灌的进行，K_{r, ∞} = 0.06 m/h时的s继续随时间逐渐增大，而堵塞效应下K_{r, ∞} < 0.06 m/h时的s则呈现一定下降趋势，且K_{r, ∞}越小，r越大，下降趋势越显著，Δs越大。至回灌后期，堵塞引起的渗透衰减趋于稳定，此时不同K_{r, ∞}下的s均随时间缓慢增大，Δs在达到峰值后基本保持不变或轻微下降。

图 3 不同K_{r, ∞}值下堵塞区内水头抬升s和水头差

$\Delta s$ 变化规律

Figure 3. Development of s and

$\Delta s$ under different values of K_{r, ∞}

下载: 全尺寸图片幻灯片

图 4所示为不同K_{r, ∞}值对含水层水头抬升分布的影响，其中图 4（a）为不同时刻z = 10 m深度处的水头抬升s的径向分布曲线。从图 4中可以看出，任一给定时刻t的水头抬升随径向距离r逐渐减小，且K_{r, ∞}越小，s下降趋势越显著，s的大小及分布范围也越小。同时，水头抬升在沿径向距离下降的过程中，在堵塞区域与主含水层交界面处存在明显的转折，K_{r, ∞}越小，t越大，该转折越明显。对考虑渗透衰减与否两种情形下所得水头抬升分布曲线作差，计算结果如图 4（b）所示。可以发现，水头差分布沿径向距离先迅速增大，在井周含水层与主含水层交界面达到峰值，且K_{r, ∞}越小，t越大，水头差峰值越大，之后随着r的增大逐渐下降至零。

图 4 不同K_{r, ∞}值下含水层水头抬升s和水头差

$\Delta s$ 分布规律

Figure 4. Distribution of s and

$\Delta s$ under various values of K_{r, ∞}

下载: 全尺寸图片幻灯片

图 5所示为不同K_{r, ∞}对定水头条件下试验井回灌流量的影响，其中为回灌流量Q时程曲线，为考虑渗透衰减与否所得回灌流量差 $\Delta Q$ 时程曲线。从图 5中可以看出，堵塞诱发的渗透衰减效应显著降低了回灌流量。在回灌初期，渗透衰减效应尚不明显，不同K_{r, ∞}下回灌流量相差非常微小；随着渗透衰减效应的逐渐显现，K_{r, ∞}越小，回灌流量衰减越迅速，且考虑渗透衰减与否引起的流量差随着时间增加明显；当回灌时间足够长，渗透系数不再继续衰减，含水层中的渗流趋于拟稳态，此时不同K_{r, ∞}下的Q随着时间非常缓慢地下降，而 $\Delta Q$ 也在非常缓慢地下降。

图 5 不同K_{r, ∞}值下回灌流量Q和流量差

$\Delta Q$ 变化规律

Figure 5. Development of Q and

$\Delta Q$ under different values of K_{r, ∞}

下载: 全尺寸图片幻灯片

2.3 回灌堵塞速率的影响

所示为不同堵塞速率，即不同 $\lambda$ 值，对z = 10 m深度处井周堵塞区内含水层水头抬升的影响，其中图 6（a）为不同径向距离（r = 0.2，0.5，1.0 m）处水头抬升s时程曲线。从中可以看出，任一给定位置处的s在发展过程均存在显著的下降阶段，且 $\lambda$ 越大，r越大，下降阶段发生越早，下降程度越显著。在回灌初期和后期，不同的 $\lambda$ 值对含水层水头抬升几无影响，而在回灌中期，较大的 $\lambda$ 值可降低水头抬升。将 $\lambda$ ≠ 0.06/h时与 $\lambda$ = 0.06/h时的s作差，计算结果如所示。可以发现，相对于 $\lambda$ = 0.06/h时的水头抬升，越小的 $\lambda$ 引起的水头差Δs随时间增加越迅速、所能达到峰值越大、达到峰值所需时间越长。同时，Δs峰值大小和达到峰值所需的时间分别与r成正比和反比关系。此后，随着不同 $\lambda$ 下的渗透系数趋于相同的K_{r, ∞}，由于 $\lambda$ 不同引起的Δs随时间逐渐衰减至零。

图 6 不同

$\lambda$ 值下堵塞区内水头抬升s和水头差

$\Delta s$ 变化规律

Figure 6. Development of s and

$\Delta s$ under different values of

$\lambda$

下载: 全尺寸图片幻灯片

所示为井周含水层不同 $\lambda$ 值对不同时刻含水层水头抬升s分布的影响，其中为s的径向分布曲线，为因 $\lambda$ 不同而引起的水头抬升差Δs的径向分布曲线。从图中可以看出，任一给定时刻t的水头抬升随径向距离r逐渐减小。在回灌初期（t = 10 h）和后期（t = 1000 h），不同 $\lambda$ 下的s分布曲线接近重合在一起，此时的Δs分布几乎等于零。而在回灌中期（t = 50 h）， $\lambda$ 越小，任意r处的s越大，相应的s分布范围也越大。此时， $\lambda$ 不同引起的Δs随r先增大并在井周含水层与主含水层交界面达到最大值，且 $\lambda$ 越小，Δs随r增速越大、所能达到的峰值越大。此后，随着r继续增大，Δs开始缓慢下降并逐渐减少至零。

图 7 不同

$\lambda$ 值下含水层水头抬升s和水头差

$\Delta s$ 分布规律

Figure 7. Distribution of s and

$\Delta s$ under various values of

$\lambda$

下载: 全尺寸图片幻灯片

所示为不同 $\lambda$ 对定水头条件下试验井回灌流量的影响，其中为不同 $\lambda$ 情形下的回灌流量时程曲线。从图中可以看出， $\lambda$ 的大小对于试验初期和后期的回灌量几乎没有影响，但在试验中期，较大的 $\lambda$ 可一定程度上降低回灌流量。将 $\lambda$ ≠ 0.06/h时与 $\lambda$ = 0.06/h的回灌流量时程曲线作差，计算结果如所示。因 $\lambda$ 不同而引起的回灌流量差首先随着时间增大， $\lambda$ 越小，流量差的增加速度越快，流量差幅值越大，达到峰值所需的时间越短。此后，随着试验的进行，不同 $\lambda$ 计算得到的流量差逐渐减小为零。

图 8 不同

$\lambda$ 值下回灌流量Q和流量差

$\Delta Q$ 变化规律

Figure 8. Development of Q and

$\Delta Q$ under different values of

$\lambda$

下载: 全尺寸图片幻灯片

3. 分析与讨论

本文的参数分析表明，堵塞程度同时改变了含水层水头抬升速度以及含水层的最终渗流状态，而堵塞速率则主要改变了含水层的渗流响应速度，并不影响含水层最初和最终的渗流状态。当井周含水层发生堵塞时，含水层渗透性下降增加了回灌水从试井流入含水层的难度，导致回灌流量的下降，进而引起整个含水层水头抬升不同程度的降低。同时，由于井筒定水头边界的约束作用，这种堵塞引起的水头抬升下降在井筒附近非常微小。随着远离井筒，水头抬升下降逐渐增大并在堵塞区与主含水层交界面达到最大值。此后，随着径向距离和过水断面的增大，井周堵塞导致的水头差逐渐减小并最终降为零。

上文的参数分析同时还表明，当考虑堵塞引起的渗透系数衰减时，含水层中水头抬升呈现先快速上升，继而明显下降，随后又继续缓慢上升的发展过程。这一现象与Wen等^[13]报道的抽水流量呈指数形式衰减时引起的水头变化过程非常类似。结合Wen等^[13]的研究，可以推断这种水头抬升过程中的下降趋势是由于渗透系数衰减导致恒定压力下回灌流量下降引起的。当堵塞越严重、回灌流量衰减越显著时，这种水头抬升的下降趋势越明显。由于该特征与渗透性恒定含水层中进行定水头回灌引起的承压水头持续抬升存在显著差异，这一特征可以作为判别井周含水层是否发生堵塞以及堵塞程度的判断标准。

除了水头抬升过程中出现的下降现象，井周含水层渗透性衰减还导致了水头抬升分布在井周堵塞区与主含水层交界面存在一定的折点，且该折点随着堵塞的发展愈加明显。分布曲线折点的存在主要是由于井周含水层与主含水层存在渗透差异造成的，这与均质含水层中的光滑分布曲线形成明显区别。因此，这一特征可以用于判别堵塞范围r_c的大小。需要指出的是，由于本文建立的井流模型假设r_c为定值，上述r_c的确定方法仅适用于试井周围堵塞分布在水平和竖直方向上相对均匀的情形，而对于实际工程中可能存在的堵塞分布不均匀的问题，相关模型的建立以及堵塞范围的判定还有待进一步研究。

4. 结论

本文以井周含水层为研究对象，建立了考虑井周堵塞效应的非完整井定水头回灌井流模型。通过变量代换、Laplace变换和有限余弦Fourier变换，推导了井流模型在Laplace-Fourier空间的解，利用有限余弦Fourier逆变换及基于Stehfest法的Laplace数值逆变换得到了井流模型在实时空间的解，并通过与数值解对比证明了本文求解方法的正确性。在此基础上，本文分析了非完整井定水头回灌条件下堵塞效应对渗流场的影响规律，得到以下3点结论。

（1）较小的K_{r, ∞}会减小回灌过程中水头抬升和回灌流量，并降低含水层的拟稳态渗流响应，而较大的 $\lambda$ 仅降低了回灌过程中水头抬升和回灌流量，并不影响含水层初始和达到拟稳态时的渗流状态。

（2）受定水头回灌边界的影响，因K_{r, ∞}和λ不同而引起的水头差沿径向距离先增大后减小，在井周堵塞区与主含水层交界面达到最大值，较小的K_{r, ∞}或较大的λ均会引起水头差分布与时程曲线的增大。

（3）回灌堵塞使得水头抬升分布在井周堵塞区与主含水层交界面上存在明显的折点，造成水头在抬升过程中存在明显的下降阶段，上述特征在堵塞越严重时越显著，可为回灌堵塞的评估与预测提供参考。

图 1 双自由度系统模型

Figure 1. Two degrees of freedom system model

下载: 全尺寸图片幻灯片

图 2 EI Centro地震波

Figure 2. EI Centro seismic waves

下载: 全尺寸图片幻灯片

图 3 原始稳定图(i=60)

Figure 3. Original stable graph

下载: 全尺寸图片幻灯片

图 4 剔除虚假极模态后稳定图

Figure 4. Stable graph after removing false polar modes

下载: 全尺寸图片幻灯片

图 5 频率–阻尼比簇分析图

Figure 5. Cluster analysis chart of frequency-damping ratio

下载: 全尺寸图片幻灯片

图 6 堆石坝网格图

Figure 6. FEM mesh of rockfill dam

下载: 全尺寸图片幻灯片

图 7 原始稳定图(i=100)

Figure 7. Original stable graph (i=100)

下载: 全尺寸图片幻灯片

图 8 剔除虚假极模态后稳定图

Figure 8. Stable graph after removing false polar modes

下载: 全尺寸图片幻灯片

图 9 频率–阻尼比簇分析图

Figure 9. Cluster analysis chart of frequency-damping ratio

下载: 全尺寸图片幻灯片

图 10 坝体典型断面及强震仪布置图

Figure 10. Typical cross-section with material zones in Nuozhadu core-wall rockfill dam

下载: 全尺寸图片幻灯片

图 11 坝体T2、T6、T7、T8测站的顺河向加速度记录

Figure 11. Downstream acceleration records of dam body at monitoring points T2, T6, T7 and T8

下载: 全尺寸图片幻灯片

图 12 糯扎渡大坝原始稳定图

Figure 12. Original stable graph of Nuozhadu dam

下载: 全尺寸图片幻灯片

图 13 糯扎渡大坝剔除虚假极模态后稳定

Figure 13. Stable graph after removing false polar modes

下载: 全尺寸图片幻灯片

图 14 糯扎渡频率–阻尼比簇分析图

Figure 14. Cluster analysis chart of frequency-damping ratio

下载: 全尺寸图片幻灯片

图 15 糯扎渡大坝三维有限元模型

Figure 15. 3D FEM of Nuozhadu dam

下载: 全尺寸图片幻灯片

表 1 双自由度模态参数识别结果

Table 1 Identified model parameters of two degrees of freedom system model

振型	频率/Hz		阻尼比
振型	理论值	计算值	理论值	计算值
1阶	4.653	4.614	0.107	0.105
2阶	12.458	12.856	0.234	0.246

下载: 导出CSV

表 2 均质堆石坝模态参数识别结果

Table 2 Identified model parameters of homogeneous rockfill dam

振型	频率/Hz		阻尼比
振型	计算值	FEM	计算值	FEM
1阶	1.11	1.09	0.023	0.020
2阶	1.75	1.76	0.016	0.020
3阶	1.93	1.96	0.022	0.020
4阶	2.20	2.24	0.015	0.020

下载: 导出CSV

表 3 糯扎渡大坝模态参数识别结果

Table 3 Identified model parameters of Nuozhadu dam

振型	频率/Hz			阻尼比
振型	计算值	FEM1	FEM2	计算值	FEM
1阶	1.15	1.10	1.22	0.051	0.04
2阶	1.35	1.27	1.42	0.034	0.04
3阶	1.83	1.54	1.73	0.027	0.04
4阶	2.37	1.74	1.96	0.043	0.04
注：FEM1为相互作用有限元模型,FEM2为刚性边界模型。

下载: 导出CSV

参考文献(32)

[1]	李红军, 朱凯斌, 赵剑明, 等. 基于设定地震场地相关反应谱的高土石坝抗震安全评价[J]. 岩土工程学报, 2019, 41(5): 934-941. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201905021.htm LI Hong-jun, ZHU Kai-bin, ZHAO Jian-ming, et al. Safety evaluation of high rock-fill dams subjected to earthquakes based on the site-related response spectra[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(5): 934-941. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201905021.htm
[2]	朱亚林, 孔宪京, 邹德高, 等. 河谷地形对高土石坝动力反应特性影响的分析[J]. 岩土工程学报, 2012, 34(9): 1590-1597. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201209008.htm ZHU Ya-lin, KONG Xian-jing, ZOU De-gao, et al. Effect of valley topography on dynamic response properties of high earth-rockfill dams[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(9): 1590-1597. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201209008.htm
[3]	CHENG L, ZHENG D J. The identification of a dam's modal parameters under random support excitation based on the Hankel matrix joint approximate diagonalization technique[J]. Mechanical Systems and Signal Processing, 2014, 42(1/2): 42-57.
[4]	王茂华, 迟世春, 刘振平. 考虑相互作用影响的堆石料动力参数反演[J]. 岩土工程学报, 2019, 41(10): 1967-1976. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201910029.htm WANG Mao-hua, CHI Shi-chun, LIU Zhen-ping. Back analysis of dynamic parameters of rock-fill materials considering interaction effects[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(10): 1967-1976. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201910029.htm
[5]	ALVES S W, HALL J F. System identification of a concrete arch dam and calibration of its finite element model[J]. Earthquake Engineering and Structural Dynamics, 2006, 35(11): 1321-1337. doi: 10.1002/eqe.575
[6]	LOH C H, WU T C. System identification of Fei-Tsui arch dam from forced vibration and seismic response data[J]. Journal of Earthquake Engineering, 2000, 4(4): 511-537.
[7]	LOH C H, WU T S. Identification of Fei-Tsui arch dam from both ambient and seismic response data[J]. Soil Dynamics and Earthquake Engineering, 1996, 15(7): 465-483. doi: 10.1016/0267-7261(96)00016-4
[8]	YANG J, JIN F, WANG J T, et al. System identification and modal analysis of an arch dam based on earthquake response records[J]. Soil Dynamics and Earthquake Engineering, 2017, 92: 109-121. doi: 10.1016/j.soildyn.2016.09.039
[9]	李帅, 潘坚文, 罗广衡, 等. 溪洛渡拱坝模态参数识别[J]. 水力发电学报, 2020, 39(3): 86-95. https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB202003009.htm LI Shuai, PAN Jian-wen, LUO Guang-heng, et al. Modal parameter identification for Xiluodu arch dam[J]. Journal of Hydroelectric Engineering, 2020, 39(3): 86-95. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB202003009.htm
[10]	谯雯, 罗佩, 刘国明. 基于自然激励技术和HHT变换的重力坝模态分析[J]. 水利学报, 2014, 45(8): 958-966. doi: 10.13243/j.cnki.slxb.2014.08.009 QIAO Wen, LUO Pei, LIU Guo-ming. Modal parameter identification of gravity dam based on natural excitation technique and Hilbert-Huang transform[J]. Journal of Hydraulic Engineering, 2014, 45(8): 958-966. (in Chinese) doi: 10.13243/j.cnki.slxb.2014.08.009
[11]	何蕴龙, 刘俊林, 熊堃. 汶川地震冶勒大坝动力响应规律分析[J]. 四川大学学报(工程科学版), 2009, 41(3): 157-164. https://www.cnki.com.cn/Article/CJFDTOTAL-SCLH200903026.htm HE Yun-long, LIU Jun-lin, XIONG Kun. Seismic response of Yele dam during Wenchuan earthquake[J]. Journal of Sichuan University (Engineering Science Edition), 2009, 41(3): 157-164. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SCLH200903026.htm
[12]	孔宪京, 周扬, 邹德高, 等. 汶川地震余震记录及紫坪铺面板堆石坝余震反应研究[J]. 岩土工程学报, 2011, 33(5): 673-678. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201105003.htm KONG Xian-jing, ZHOU Yang, ZOU De-gao, et al. Aftershock records of Wenchuan Earthquake and seismic response of Zipingpu Concrete Face Rock-fill Dam[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(5): 673-678. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201105003.htm
[13]	苗君, 何蕴龙, 曹学兴, 等. 芦山地震冶勒大坝强震监测资料分析[J]. 岩土力学, 2015, 36(1): 225-232, 256. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201501031.htm MIAO Jun, HE Yun-long, CAO Xue-xing, et al. Analysis of strong motion seismograph data at rockfill Yele dam during Lushan earthquake[J]. Rock and Soil Mechanics, 2015, 36(1): 225-232, 256. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201501031.htm
[14]	刘振平, 迟世春, 赵显波, 等. 鲤鱼潭大坝坝料动力参数反演[J]. 岩土工程学报, 2015, 37(4): 761-768. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201504029.htm LIU Zhen-ping, CHI Shi-chun, ZHAO Xian-bo, et al. Back analysis of dynamic parameters of Liyutan dam materials[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(4): 761-768. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201504029.htm
[15]	王茂华, 迟世春, 相彪, 等. 弱震情况下高土石坝坝料动力参数反演分析[J]. 岩土工程学报, 2020, 42(2): 289-298. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202002014.htm WANG Mao-hua, CHI Shi-chun, XIANG Biao, et al. Back analysis of dynamic parameters of high earth-rock dam materials under weak earthquakes[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(2): 289-298. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202002014.htm
[16]	李平, 薄景山, 李孝波, 等. 安宁河河谷及邛海地区土层场地对地震动的放大作用[J]. 岩土工程学报, 2016, 38(2): 362-369. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201602026.htm LI Ping, BO Jing-shan, LI Xiao-bo, et al. Amplification effect of soil sites on ground motion in Anning River valley and Qionghai Lake area[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(2): 362-369. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201602026.htm
[17]	章国稳. 环境激励下结构模态参数自动识别与算法优化[D]. 重庆: 重庆大学, 2012. ZHANG Guo-wen. Modal Parameter Automatic Identification for Structures Under Ambient Excitation and Algorithm Optimization[D]. Chongqing: Chongqing University, 2012. (in Chinese)
[18]	REYNDERS E, ROECK G D. Reference-based combined deterministic-stochastic subspace identification for experimental and operational modal analysis[J]. Mechanical Systems and Signal Processing, 2008, 22(3): 617-637.
[19]	ZHANG L M, WANG T, TAMURA Y. A frequency- spatial domain decomposition (FSDD) method for operational modal analysis[J]. Mechanical Systems and Signal Processing, 2010, 24(5): 1227-1239.
[20]	张永祥, 刘心, 褚志刚, 等. 基于随机子空间法的模态参数自动提取[J]. 机械工程学报, 2018, 54(9): 187-194. https://www.cnki.com.cn/Article/CJFDTOTAL-JXXB201809023.htm ZHANG Yong-xiang, LIU Xin, CHU Zhi-gang, et al. Autonomous modal parameter extraction based on stochastic subspace identification[J]. Journal of Mechanical Engineering, 2018, 54(9): 187-194. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JXXB201809023.htm
[21]	汤宝平, 章国稳, 陈卓. 基于谱系聚类的随机子空间模态参数自动识别[J]. 振动与冲击, 2012, 31(10): 92-96. https://www.cnki.com.cn/Article/CJFDTOTAL-ZDCJ201210021.htm TANG Bao-ping, ZHANG Guo-wen, CHEN Zhuo. Automatic stochastic subspace identification of modal parameters based on hierarchical clustering method[J]. Journal of Vibration and Shock, 2012, 31(10): 92-96. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-ZDCJ201210021.htm
[22]	MAGALHÃES F, CUNHA Á, CAETANO E. Online automatic identification of the modal parameters of a long span arch bridge[J]. Mechanical Systems and Signal Processing, 2009, 23(2): 316-329.
[23]	PEETERS B, DE ROECK G. Reference-based stochastic subspace identification for output-only modal analysis[J]. Mechanical Systems and Signal Processing, 1999, 13(6): 855-878.
[24]	BAKIR P G. Automation of the stabilization diagrams for subspace based system identification[J]. Expert Systems with Applications, 2011, 38(12): 14390-14397.
[25]	HONG A L, UBERTINI F, BETTI R. New stochastic subspace approach for system identification and its application to long-span bridges[J]. Journal of Engineering Mechanics, 2013, 139(6): 724-736.
[26]	大崎顺彦. 振动理论[M]. 谢礼立,译. 北京: 地震出版社, 1990. AKIHIKO O. Vibration Theory[M]. XIE Li-yi, trans. Beijing: Seismological Press, 1990. (in Chinese)
[27]	徐斌, 邹德高, 孔宪京, 等. 高土石坝坝坡地震稳定分析研究[J]. 岩土工程学报, 2012, 34(1): 139-144. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201201015.htm XU Bin, ZOU De-gao, KONG Xian-jing, et al. Seismic stability of slopes of high rockfiU dams[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(1): 139-144. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201201015.htm
[28]	孔宪京, 娄树莲, 邹德高, 等. 筑坝堆石料的等效动剪切模量与等效阻尼比[J]. 水利学报, 2001, 32(8): 20-25. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB200108003.htm KONG Xian-jing, LOU Shu-lian, ZOU De-gao, et al. The equivalent dynamic shear modulus and equivalent damping ratio of the rockfill for dam[J]. Journal of Hydraulic Engineering, 2001, 32(8): 20-25. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB200108003.htm
[29]	中国水电顾问集团昆明勘测设计研究院. 200 m级以上高心墙堆石坝坝料特性及坝料设计准则研究[R]. 昆明: 中国水电顾问集团昆明勘测设计研究院, 2010. Hydro China Kunming Engineering Corporation. Study on Dam Material Characteristics and Design Criteria of Dam Materials for High Core Wall Rockfill Dam Above 200 m[R]. Kunming: Hydro China Kunming Engineering Corporation, 2010. (in Chinese)
[30]	楚金旺, 朱晟, 黄亚梅. 基于实际震害的土石坝永久变形估算[J]. 中国水利水电科学研究院学报, 2017, 15(6): 409-417. https://www.cnki.com.cn/Article/CJFDTOTAL-ZGSX201706001.htm CHU Jin-wang, ZHU Sheng, HUANG Ya-mei. Estimation of earthquake-induced permanent deformation for earth dam based on seismic damage[J]. Journal of China Institute of Water Resources and Hydropower Research, 2017, 15(6): 409-417. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-ZGSX201706001.htm
[31]	毛雯娟. 强震区高土石坝三维动力反应分析[D]. 大连: 大连理工大学, 2008. MAO Wen-juan. 3-D Dynamic Analysis of High Rockfill Dam in Meizoseismal Area[D]. Dalian: Dalian University of Technology, 2008. (in Chinese)
[32]	杨玉生, 刘小生, 刘启旺, 等. 双江口心墙堆石坝动力特性的振动台模型试验研究[J]. 水力发电学报, 2011, 30(1): 114-119. https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB201101020.htm YANG Yu-sheng, LIU Xiao-sheng, LIU Qi-wang, et al. Shaking table model tests on dynamic characteristics of Shuangjiangkou high earth-rockfill dam[J]. Journal of Hydroelectric Engineering, 2011, 30(1): 114-119. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB201101020.htm