CN111723477A - 一种基于动态层析扫描的管涌发展预测预警方法 - Google Patents

Abstract

本发明公开了一种基于动态层析扫描的管涌发展预测预警方法，包括如下步骤：解译含水层渐进侵蚀机理与渗流场动态响应规律；建立含水层渐进侵蚀流固耦合模型；布置监测井收集监测数据；建立含水层渐进侵蚀多源监测数据解译模型；计算当前渗流场；管涌流固耦合模型计算侵蚀情况；更新渗流场并迭代计算至步长要求。本发明方法结合了管涌流固耦合模型的特点以及水力层析扫描刻画非均质含水层的优势，对管涌通道动态发展进行预测，得到的结果相较于传统的监测方式，在准动态反应、实时反馈、精确定位等方面有较大的提升。

Description

一种基于动态层析扫描的管涌发展预测预警方法

技术领域

本发明涉及地层岩性、地质构造地球物理勘探领域，具体涉及 一种基于动态层析扫描的管涌发展预测预警方法。

背景技术

渗透破坏泛指堤坝内部或底部在渗流作用下颗粒流失并进一步 造成破坏的过程。世界上有46.1％的堤坝破坏都是由渗透破坏造成 的。根据最新的分类，渗透破坏可以分为集中渗漏侵蚀、接触侵蚀、 潜蚀、逆向侵蚀管涌等四种类型。逆向侵蚀管涌是造成堤防失稳破 坏的主要原因。历史上，我国的主要堤防，包括荆江大堤、北江大 堤、淮河大堤等都曾因管涌发生破坏。尤其在长江中下游，特殊的 二元堤基结构在洪峰过境时面临巨大的管涌风险。如今，我国已建 成各类堤防近30余万公里，因降雨分布不平均，每年仍有大量堤防遭受管涌威胁或因管涌发生破坏甚至溃决，管涌依然是严重威胁堤 防及两岸人民生命财产安全的重要问题。

探测管涌发生、监测管涌发展、预测管涌风险，是管涌防治的 关键，然而，由于管涌的隐蔽性与随机性，长期以来汛期管涌防范 依然依靠群众性巡查来实现。这不仅无法确保及时发现管涌，也无 法准确获取管涌通道发展的情况，防汛投入大但效果差。部分堤防已逐步试点安全监测预警系统，力求利用数据监测及物探手段，推 断管涌通道的置与规模，监测管涌发生发展。但在监测过程中，由 于缺少监测数据与管涌发展之间的定量关系，只能根据数据分布与 变化规律定性分析。逆向侵蚀管涌的物理本质是含水层的渐进式侵蚀。管涌口形成后，在渗流作用下堤基土体细颗粒发生移动并流失， 形成管涌通道并沿渗流相反的方向逐渐向上游发展，逐渐形成集中 管涌通道。颗粒侵蚀以及管涌通道的形成改变了含水层的水力特性， 进而引起渗流场重分布，并进一步影响管的发生发展。因此，逆向侵蚀管涌是一个流固耦合过程。因管涌物理过程的耦合特性、管涌 土体组成的多样性、管涌发展水力环境的复杂性等，逆向侵蚀管涌 的机理尚不清楚，管涌数学模型也需要进一步完善。因此，揭示管 涌发生发展机理、建立管涌预测预警模型、实现管涌自动在线识别， 对保障堤防两岸人民生命财产安全具有重要的意义。

发明内容

本发明的目的是提供一种基于动态层析扫描的管涌发展预测预 警方法。

为实现本发明目的而采用的技术方案是这样的，一种基于动态 层析扫描的管涌发展预测预警方法，包括以下步骤：

1)布置监测井收集监测数据：在监测区域布置若干监测井，若 干监测井呈网状分布，在每个监测井内设置监测点，监测点接收水 位、水温、流速流向的数据；

2)多源监测数据融合：通过贝叶斯理论将步骤1)中收集到的时 序响应数据进行融合，构建多源信息融合框架，建立含水层渐进侵 蚀多源监测数据的解译模型；

3)利用连续线性估计算法求解地下水运移方程：利用在前期勘 察资料中得到的先验信息，计算出场地内水力参数随机分布规律， 结合步骤2)中的多源信息融合框架，通过连续线性估计算法得到堤 基含水层渗流场的空间分布；

4)基于动态层析扫描的堤基逆向侵蚀管涌实时刻画：利用分区 界面上的物质交换关系，结合堤基管涌初始边值条件特点提出管涌 耦合模型的高效求解方法；通过流固耦合数学模型以及管涌数学模 型参数敏感度求解高效算法共同求解在步骤3)给出的渗流场中管涌 的发展情况，并计算颗粒流失情况；

5)通过步骤4)中的颗粒流失状况，更新水力参数场；

6)利用步骤5)中更新的水利参数场重复步骤4)和5)，直到计算 时长达到要求；

7)通过在线监测数据的实时传输与快速的在线求解算法，实现 监测数据的快速解译及管涌生成的实时在线识别。

进一步，步骤1)中所述监测井打穿监测区域的含水层，监测井 的内径为10cm，监测井沿其深度方向布设若干监测点，相邻监测点 的间距小于等于1m。

进一步，步骤2)包括以下步骤：

2.1)将步骤1)获得的渗透系数分布场T(x)和水头分布场H(x)进 行分解：

lnT(x)＝F+f(x) (1)

式(1)中：F为渗透系数分布场的均值，f(x)为均值F下的扰 动；

lnH(x)＝H+h(x) (2)

式(2)中：H为水头分布场的均值，h(x)为均值H下的扰动；

2.2)监测区域内任意一点的渗透系数的估计值为：

其中，λ_mi为渗透系数分布场对该点的贡献量，β_mj为水头分布场 对该点的贡献量，N_f为渗透系数分布场测量的总个数，N_h为水头分 布场测量的总个数。

进一步，步骤3)包括以下步骤：

3.1)根据前期勘探数据，获取含水层参数的均值、方差和相关尺 度；

3.2)将含水层参数的均值赋予分析区域内每一个网格内，结合边 界条件，建立监测区域三维正分析模型，利用监测区域三维正分析 模型，计算每一次刺激下每个探测位置在估计参数条件下的水头以 及流量；

3.3)利用连续线性估计算法对探测区域内待估点的含水层参数 值进行估计，未知参数的迭代计算公式为：

式中，u_c为待估的监测区域含水层参数向量；u_c ^(r+1)为参数向量 u_c在第r+1次的条件估计值，r＝0的参数来源于步骤3.2)正分析的结 果以及先验资料。d^*为每次刺激作用下的监测值，d^(r)为三维正分析 模型的模拟值，ω^(r)T表示权重系数，T为反演周期，权重系数矩阵ω 计算公式为：

[ε_dd+λdiag(ε_dd)]ω＝ε_du (5)

其中ε_dd是观测数据之间的协方差矩阵，而ε_du是观测数据与参 数之间的协方差矩阵。λ是Levenberg-Marquardt算法动态乘子，而 diag()运算符代表取对角阵。协方差由敏感度得到：

式中，J_du是观测数据对探测区域参数变化的敏感度矩阵，参数 协方差矩阵ε_uu在r＝0时由先验地质信息给出，随后每次迭代按照下 式逐步更新：

3.4)重复式(4)的迭代过程，直至得到的参数估计值条件下得 到的计算值与探测值之差小于等于0.1％，结束迭代；

3.5)通过达西定律计算渗流场在该时刻的空间分布：

式中，k_x、k_y、k_z分别为x、y、z方向上的渗透系数，v_x、v_y、 v_z分别为x、y、z方向上的渗透速度，

分别表示水头分 布场H对x、y、z的偏导。

进一步，步骤4)包括以下步骤：

4.1)通过伴随方程求解管涌流固耦合模型参数敏感度矩阵：

将敏感度方程乘以一个任意的可微参数

然后分别对体积和 时间积分，即：

式中，式中，K(x)为渗透系数，S_S是贮水系数，Q(x)是流量， Ω和T是空域和时域，α是参数向量，k＝1,…N,N是参数总数量；

通过格林恒等式将式(8)转变为：

通过边界条件以及状态方程，消去敏感度

4.2)通过计算得到的渗流场分布计算管涌正分析数学模型：

式(10a)表示流固耦合模型中的水流方程，式(10b)表示流 固耦合模型中的固态向的方程。其中，H是水头；S_S是贮水系数； C_m是土体侵蚀后运动颗粒的质量浓度；V是步骤3.5)中求得的渗流 场分布；I是单位时间内土体中颗粒的流失量；采用幂函数表示：

式中，α、β为常系数，τ_c为管涌临界条件；在非均质含水层中 渗透系数k、临界条件τ_c、系数α、β均具有空间变异性；

4.3)对管涌区域网格进行自适应划分，在管涌区使用加密网格， 在非管涌区使用稀疏网格。

本发明的技术效果是毋庸置疑的，本发明的层析扫描技术是结 合了管涌流固耦合模型的特点以及水力层析扫描刻画非均质含水层 的优势，并对管涌通道动态发展进行预测。得到的结果相较于传统 的监测方式，在动态反应、计算速度、精确定位等方面有较大的提 升。

附图说明

图1为本发明的工作流程图；

图2为求解域及网格划分示意图；

图3为隐伏渗漏通道随时间的发育情况。

具体实施方式

下面结合实施例对本发明作进一步说明，但不应该理解为本 发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思 想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换 和变更，均应包括在本发明的保护范围内。

实施例1：

参见图1，本实施例公开一种基于动态层析扫描的管涌发展预测 预警方法，包括以下步骤：

1)解译含水层渐进侵蚀机理与渗流场动态响应规律；解译在管 涌发展过程中含水层参数与水利要素的关系；解译含水层参数自相 关性和互相关性对数据时序响应的影响；

2)自动识别堤基含水层的分区情况，并利用物质输移过程中的 质量守恒、动量守恒以及能量守恒等手段解释不同分区内的水流与 颗粒运动规律的关系；

3)在监测区域布置监测井收集监测数据：在监测区域内布置监 测井，监测井上设置有监测点，用于接收水位、水温、流速流向的 数据；所述监测井打穿监测区域的含水层，监测井的内径为10cm， 监测井沿其深度方向布设若干监测点，相邻监测点的间距小于等于1m；

4)多源监测数据融合：通过贝叶斯理论将第3)步中收集到的各 种时序响应数据进行融合，构建多源信息融合框架，建立含水层渐 进侵蚀多源监测数据的解译模型；具体包括以下步骤：

4.1)假设在步骤1)中获得了渗透系数分布场T(x)和水头分布场 H(x)，对这两个物理量进行分解：

lnT(x)＝F+f(x) (1)

式(1)中：F为渗透系数分布场的均值，f(x)为均值F下的扰 动；

lnH(x)＝H+h(x) (2)

式(2)中：H为水头分布场的均值，h(x)为均值H下的扰动；

4.2)监测区域内任意一点的渗透系数的估计值为：

其中，λ_mi为渗透系数分布场对该点的贡献量，β_mj为水头分布场 对该点的贡献量，N_f为渗透系数分布场测量的总个数，N_h为水头分 布场测量的总个数，通过该式完成了二元信息的融合，当信息量增 加时，只需相应增加参数个数即可。

5)利用连续线性估计算法求解地下水运移方程：利用在前期勘 察资料中得到的先验信息，计算出场地内水力参数随机分布规律， 并结合第4)步中的多源信息融合框架通过连续线性估计算法得到堤 基含水层渗流场的空间分布；具体包括以下步骤：

5.1)根据前期勘探数据，获取含水层参数的均值、方差和相关尺 度；勘探数据包括但不局限于现场测量，文献资料等获取渠道，主 要用于估计各种物理参数的均值、方差和尺度，用于正分析模型中；

5.2)将含水层参数的均值赋予分析区域内每一个网格内，结合边 界条件，建立监测区域三维正分析模型，利用监测区域三维正分析 模型，计算每一次刺激下每个探测位置在估计参数条件下的水头以 及流量；

5.3)利用连续线性估计算法对探测区域内待估点的含水层参数 值进行估计，需要注意的是此处的待估点并不是监测点，连续线性 估计算法是通过监测点以及刺激点收集到的数据对待估点的值进行 估计，待估点可以是试验场地内的任意一点；未知参数的迭代计算 公式为：

式中，u_c为待估的监测区域含水层参数向量；u_c ^(r+1)为参数向量 u_c在第r+1次的条件估计值，r＝0的参数来源于步骤5.2)正分析的结 果以及先验资料，根据求解的问题不同，先验资料中包括水头、渗 透系数、电阻率等物理参数。d^*为每次刺激作用下的监测值，d^(r)为 三维正分析模型的模拟值，ω^(r)T表示权重系数，T为反演周期，权重 系数矩阵ω计算公式为：

[ε_dd+λdiag(ε_dd)]ω＝ε_du (5)

其中ε_dd是观测数据之间的协方差矩阵，而ε_du是观测数据与参 数之间的协方差矩阵。λ是Levenberg-Marquardt算法动态乘子，而 diag()运算符代表取对角阵。协方差由敏感度得到：

式中，J_du是观测数据对探测区域参数变化的敏感度矩阵，参数 协方差矩阵ε_uu在r＝0时由先验地质信息给出，随后每次迭代按照下 式逐步更新：

5.4)重复式(4)的迭代过程，直至得到的参数估计值条件下得 到的计算值与探测值之差小于等于0.1％，结束迭代，当计算结果不 能满足误差要求时，则根据预先设置的迭代步数来判断是否结束迭 代步骤；

5.5)通过达西定律计算渗流场在该时刻的空间分布：

式中，k_x、k_y、k_z分别为x、y、z方向上的渗透系数，v_x、v_y、 v_z分别为x、y、z方向上的渗透速度，

分别表示水头分 布场H对x、y、z的偏导。

6)基于动态层析扫描的堤基逆向侵蚀管涌实时刻画：利用分区 界面上的物质交换关系，结合堤基管涌初始边值条件特点提出管涌 耦合模型的高效求解方法；通过流固耦合数学模型以及管涌数学模 型参数敏感度求解高效算法共同求解在第5)步给出的渗流场中管涌 的发展情况，并计算颗粒流失情况；具体包括以下步骤：

6.1)通过伴随方程求解管涌流固耦合模型参数敏感度矩阵：

由于使用传统方法计算敏感度时要依次求解个单元的参数值会 浪费大量计算资源，伴随方程法的核心在于不去对每个单元求解敏 感度方程，即消掉敏感度的影响。为实现这一目的，将敏感度方程 乘以一个任意的可微参数φ，然后分别对体积和时间积分，即：

式中，K(x)为渗透系数，S_S是贮水系数，Q(x)是流量，Ω和 T是空域和时域，α是参数向量，k＝1,…N,N是参数总数量；

通过格林恒等式将式(8)转变为：

通过边界条件以及状态方程，消去敏感度φ；

6.2)通过计算得到的渗流场分布计算管涌正分析数学模型：

式(10a)表示流固耦合模型中的水流方程，式(10b)表示流 固耦合模型中的固态向的方程。其中，H是水头；S_S是贮水系数； C_m是土体侵蚀后运动颗粒的质量浓度；V是步骤3.5)中求得的渗流 场分布；I是单位时间内土体中颗粒的流失量；采用幂函数表示：

式中，α、β为常系数，τ_c为管涌临界条件；在非均质含水层中 渗透系数k、临界条件τ_c、系数α、β均具有空间变异性；

6.3)对管涌区域网格进行自适应划分，在管涌区使用加密网格， 在非管涌区使用稀疏网格。

7)通过第6)步中的颗粒流失状况，更新水力参数场；

8)利用步骤7)中更新的水利参数场重复步骤第6)和7)，直到计 算时长达到要求；

9)通过在线监测数据的实时传输与快速的在线求解算法，实现 监测数据的快速解译及管涌生成的实时在线识别。

实施例2：

本实施例主要步骤同实施例1，进一步，选取一块40m×60m的 含水层区域作为研究区域并对其进行数值建模，利用层析扫描方法 对含水层的渗透系数分布进行刻画。并假设在模型的左边界存在一 高水头区域，即为堤防中的水位边界，模拟一从左到右的渗流场。 在边界围成的求解域内进行单元网格划分，网格为正方形，尺度为 1m×1m，并共计布置了56个监测电极，并预设了6次放电-接收试 验方案，网格划分情况如图2所示。

规定了堤防内部构造的均值、方差与相关尺度，通过多源响应 数据对管涌发展过程进行模拟。将模拟步长设置为1天，每隔T＝1d 显示管涌通道的发展情况，图3代表了多次模拟过程中渗漏通道逐 渐发展的一个过程，分别展示了在第1天，第4天，第8天，第10 天的管涌通道发展情况以及渗流场的变化情况。

由图可知随着时间的发展，在左边高水头右边低水头情况下， 管涌通道逐渐形成，在管涌通道形成的过程中，其渗流场也在不断 的改变。在右侧的图中，展示了渗流场中流速流向以及水头的变化 趋势，可以看出左右两幅图是相互影响的，这也体现了管涌发展过 程中的流固耦合规律。

值得说明的是，动态层析扫描技术是结合了管涌流固耦合模型 的特点以及水力层析扫描刻画非均质含水层的优势，并对管涌通道 动态发展进行预测。得到的结果相较于传统的监测方式，在动态反 应、计算速度、精确定位等方面有较大的提升。

Claims (5)

1.一种基于动态层析扫描的管涌发展预测预警方法，其特征在于：包括以下步骤：

1)布置监测井收集监测数据：在监测区域布置若干监测井，若干监测井呈网状分布，在每个监测井内设置监测点，监测点接收水位、水温、流速流向的数据。

2)多源监测数据融合：通过贝叶斯理论将步骤1)中收集到的时序响应数据进行融合，构建多源信息融合框架，建立含水层渐进侵蚀多源监测数据的解译模型；

3)利用连续线性估计算法求解地下水运移方程：利用在前期勘察资料中得到的先验信息，计算出场地内水力参数随机分布规律，结合步骤2)中的多源信息融合框架，通过连续线性估计算法得到堤基含水层渗流场的空间分布；

4)基于动态层析扫描的堤基逆向侵蚀管涌实时刻画：利用分区界面上的物质交换关系，结合堤基管涌初始边值条件特点提出管涌耦合模型的高效求解方法；通过流固耦合数学模型以及管涌数学模型参数敏感度求解高效算法共同求解在步骤3)给出的渗流场中管涌的发展情况，并计算颗粒流失情况；

5)通过步骤4)中的颗粒流失状况，更新水力参数场；

6)利用步骤5)中更新的水利参数场重复步骤4)和5)，直到计算时长达到要求；

7)通过在线监测数据的实时传输与快速的在线求解算法，实现监测数据的快速解译及管涌生成的实时在线识别。

2.根据权利要求1所述的一种基于动态层析扫描的管涌发展预测预警技术，其特征在于：步骤1)中所述监测井打穿监测区域的含水层，监测井的内径为10cm，监测井沿其深度方向布设若干监测点，相邻监测点的间距小于等于1m。

3.根据权利要求1所述的一种基于动态层析扫描的管涌发展预测预警技术，其特征在于：步骤2)包括以下步骤：

2.1)将步骤1)获得的渗透系数分布场T(x)和水头分布场H(x)进行分解：

lnT(x)＝F+f(x) (1)

式(1)中：F为渗透系数分布场的均值，f(x)为均值F下的扰动；

lnH(x)＝H+h(x) (2)

式(2)中：H为水头分布场的均值，h(x)为均值H下的扰动；

2.2)监测区域内任意一点的渗透系数的估计值为：

其中，λ_mi为渗透系数分布场对该点的贡献量，β_mj为水头分布场对该点的贡献量，N_f为渗透系数分布场测量的总个数，N_h为水头分布场测量的总个数。

4.根据权利要求1所述的一种基于动态层析扫描的管涌发展预测预警技术，其特征在于：步骤3)包括以下步骤：

3.1)根据前期勘探数据，获取含水层参数的均值、方差和相关尺度；

3.2)将含水层参数的均值赋予分析区域内每一个网格内，结合边界条件，建立监测区域三维正分析模型，利用监测区域三维正分析模型，计算每一次刺激下每个探测位置在估计参数条件下的水头以及流量；

3.3)利用连续线性估计算法对探测区域内待估点的含水层参数值进行估计，未知参数的迭代计算公式为：

式中，u_c为待估的监测区域含水层参数向量；u_c ^(r+1)为参数向量u_c在第r+1次的条件估计值，r＝0的参数来源于步骤3.2)正分析的结果以及先验资料。d^*为每次刺激作用下的监测值，d^(r)为三维正分析模型的模拟值，ω^(r)T表示权重系数，T为反演周期，权重系数矩阵ω计算公式为：

[ε_dd+λdiag(ε_dd)]ω＝ε_du (5)

其中ε_dd是观测数据之间的协方差矩阵，而ε_du是观测数据与参数之间的协方差矩阵。λ是Levenberg-Marquardt算法动态乘子，而diag()运算符代表取对角阵。协方差由敏感度得到：

式中，J_du是观测数据对探测区域参数变化的敏感度矩阵，参数协方差矩阵ε_uu在r＝0时由先验地质信息给出，随后每次迭代按照下式逐步更新：

3.4)重复式(4)的迭代过程，直至得到的参数估计值条件下得到的计算值与探测值之差小于等于0.1％，结束迭代；

3.5)通过达西定律计算渗流场在该时刻的空间分布：

式中，k_x、k_y、k_z分别为x、y、z方向上的渗透系数，v_x、v_y、v_z分别为x、y、z方向上的渗透速度，

分别表示水头分布场H对x、y、z的偏导。

5.根据权利要求1所述的一种基于动态层析扫描的管涌发展预测预警技术，其特征在于：步骤4)包括以下步骤：

4.1)通过伴随方程求解管涌流固耦合模型参数敏感度矩阵：

将敏感度方程乘以一个任意的可微参数

然后分别对体积和时间积分，即：

式中，K(x)为渗透系数，S_S是贮水系数，Q(x)是流量，Ω和T是空域和时域，α是参数向量，k＝1,…N,N是参数总数量；

通过格林恒等式将式(8)转变为：

通过边界条件以及状态方程，消去敏感度

4.2)通过计算得到的渗流场分布计算管涌正分析数学模型：

式(10a)表示流固耦合模型中的水流方程，式(10b)表示流固耦合模型中的固态向的方程。其中，H是水头；S_S是贮水系数；C_m是土体侵蚀后运动颗粒的质量浓度；V是步骤3.5)中求得的渗流场分布；I是单位时间内土体中颗粒的流失量；采用幂函数表示：

式中，α、β为常系数，τ_c为管涌临界条件；在非均质含水层中渗透系数k、临界条件τ_c、系数α、β均具有空间变异性；

4.3)对管涌区域网格进行自适应划分，在管涌区使用加密网格，在非管涌区使用稀疏网格。

Images