基于热源示踪与水力联合层析反演的致灾构造超前预报方法

1.本发明涉及地下工程领域致灾构造超前预报技术领域，尤其是一种基于热源示踪与水力联合层析反演的致灾构造超前预报方法。


背景技术：

2.随着水利水电工程与铁路公路工程的建设重点向复杂地形地质条件下的西部山区转移，地下工程建设呈现出埋深大、地质复杂、灾害频发的特点，工程施工中易遭遇突水、突泥等地质灾害。岩体中发育的不良地质体具有十分显著的致灾作用，包括溶洞、岩溶管道、暗河等岩溶不良地质体、天然和人工洞穴、断层破碎带、强风化带等特殊岩土体等。其中，含导水断层及充水溶洞等构造在水头压力和施工扰动的共同作用下，隔水岩层或阻水结构可能发生失稳，进而诱发突水突泥灾害，成为造成突水突泥灾害的主要原因。地下工程施工需要提前准确探明工程前方赋存的富水致灾构造，防患于未然。
3.应用地球物理领域的传统勘探手段是开展灾害源探测常采用的手段之一，其所获取的电导率、磁导率等参数，与地层水力特性间的转换过程较为复杂且模糊，同时相关转换常与具体地层条件有关，难以为具体工程提供普适性参考。结合水文地质学领域最新研究，针对地下工程开展示踪试验及抽注水试验等传统水文地质试验所获取的数据信息，包含了许多描述致灾构造富水性的非冗余信息，结合水头降深和示踪走时数据开展联合反演，为致灾构造超前探测研究提供了新思路。近年来，温度作为一种低廉、便捷且无污染的新型示踪剂开始得到关注，双封隔分层成井技术、分布式光纤温度传感器的迅速发展也为高精度、低成本的原位实时测量创造了可能，为水头、温度示踪数据监测提供了技术支持。在此背景下，亟待结合热源示踪与水力信息研究一种适用于致灾构造的联合反演超前预报方法。


技术实现要素：

4.本发明所要解决的技术问题在于，提供一种基于热源示踪与水力联合层析反演的致灾构造超前预报方法，用于获取水文地质参数分布及探查富水致灾构造位置，为工程安全施工提供保障。
5.为解决上述技术问题，本发明提供一种基于热源示踪与水力联合层析反演的致灾构造超前预报方法，包括如下步骤：
6.(1)基于具体工程研究区的地质资料、物探及钻探资料，针对富水构造破碎带、岩溶发育区段及短中距离预报中的物探异常区域，结合超前地质探测施作的超前钻孔长度、深度及倾角情况，圈定掌子面前方监测区域；随后，结合地下工程超前钻探布置位置、数量及隧道、巷道洞径，确定热源示踪试验及分层补源水力试验监测方案，包括监测组数及每组数量、监测间距；
7.(2)根据已确定的监测区域及试验监测方案，实施热源示踪试验及抽注水水力试验，监测记录示踪穿透曲线及水头演化曲线；
8.(3)开展曲线平滑消噪滤波处理，提取示踪旅行时及水头降深数据；
9.(4)考虑到温度传播过程存在热扩散效应及弥散效应，基于温度穿透时间曲线的一阶导数变化最速时刻确定延迟因子，实现示踪旅行时纠偏。为规避早期旅行时选取的不确定性，针对纠偏后的旅行时数据集，进一步开展早期旅行时诊断分析；
10.(5)基于旅行时数据集和联合迭代重建算法，开展迭代实现介质特征重构，选取最佳重建旅行时及估计参数分布；
11.(6)结合基于旅行时得到的估计水文地质参数分布，确定含水层参数的初始估计均值、方差及相关尺度等，作为进一步开展水力层析分析的先验信息。
12.(7)结合同时连续线性估计算法，开展水力层析反演分析，逐次迭代直至满足阈值设定，随后计算确定渗透系数、储水系数及导压系数等水文地质参数分布，结合高渗区圈定富水致灾构造位置及分布。
13.优选的，步骤(1)中，取掌子面前方最长钻孔的1.2-1.5倍深度范围为监测区域，为30米至100米范围，如工程所需探测精度较高且研究区岩溶发育程度高、富水性强，采用的试验监测间距为0.5-1倍洞径，且最大采用2-3倍洞径作为监测间距；如超前钻孔已揭示潜在富水致灾构造，则应结合潜在致灾构造区域处3倍洞径范围内的多个钻孔，并尽可能加密监测点，缩小监测间距为0.3-0.5倍钻孔最小间距，重新试验开展精细探测以确定致灾构造规模及特征。
14.优选的，步骤(2)中，热源示踪试验采用单点投放、多点监测的形式，热水水源投放宜选取在掌子面中心钻孔的监测点内实施，所用热源温度应高于环境背景温度40-60摄氏度，且注入时长为监测总时长的一半后，距离注入端口最近的监测端口温度宜上升20摄氏度以上；热源示踪试验停止注入后，在保证前序试验监测孔温度恢复至初始背景温度后方可进行后续试验，温度变化低于5％视为恢复。
15.优选的，步骤(2)中，水力试验采用单点或多点抽注，同时多点持续监测，且开展实验不少于三次；当监测点数少于十个时，开展抽注水实验不少于五次；当监测点数多于十个，则适当减少试验次数；试验早期采用短监测时间间隔，具体间隔短于抽注水试验总时间的1％，尽可能多的收集早期水头数据，并在水头达到稳定后停止试验，选取设计试验总时间的四分之一时长作为抽注水早期阶段。
16.优选的，步骤(3)中，开展曲线平滑消噪滤波处理，提取示踪旅行时及水头降深数据具体包括如下步骤：
17.(31)研究设计窗口加权随方差变化的平滑低通滤波函数，开展预消噪处理，进一步基于常规卡尔曼滤波与强跟踪卡尔曼滤波算法的自适应卡尔曼滤波方法，依据噪声动态变化，实时调整状态预报误差的协方差矩阵及增益矩阵，开展深度滤波消噪；
18.(32)针对平滑低通滤波，就温度示踪穿透曲线采用高斯滤波函数进行消噪，根据窗口内数据方差与整体数据集总方差的相对大小确定高斯函数方差，当窗口内数据方差大于数据集总方差时，高斯函数方差选取为1.5，反之取为0.5，进而结合高斯函数方差，确定窗口内各数据的权重向量，根据权重及实测温度数据获取预消噪温度曲线结果；
19.(33)针对平滑低通滤波，就水头曲线采用五阶或六阶多项式函数进行拟合消噪，且所用多项式阶次需小于窗口长度，进而对给定高阶多项式开展最小二乘拟合，逐步移动平均得到预消噪水头曲线；
20.(34)假设系统噪声和观测噪声互不相关，结合滤波本身开展系统噪声方差及量测
噪声方差的估计和修正，引入时变衰减因子l
k+1
依据噪声动态变化，修正常规卡尔曼滤波中的预测协方差矩阵u
k+1,k
，实时调整状态预报误差的协方差矩阵及增益矩阵；
[0021][0021][0022]
其中，trace[]为求矩阵迹，ik为单位矩阵，为状态转移矩阵，uk为滤波协方差矩阵，v
k-1
为残差协方差矩阵，zk为新息序列，rk为系统噪声协方差矩阵，qk为量测噪声协方差矩阵；
[0023]
(35)针对消噪后的温度示踪穿透曲线，开展蒙特卡洛分析提取峰值旅行时数据；针对消噪后的水头曲线提取降深数据，将开始试验后四分之一总时长视为抽水试验早期，停止抽水后至水头完全恢复视为抽水试验晚期，其余时间为中期；结合早、中、晚期数据，在每组水头曲线选取不少于五个降深数据，其中在早、晚期时间内至少各选取两个数据。
[0024]
优选的，步骤(4)中，确定温度函数导数峰值时间tb与峰值旅行时td的定量关系，进一步得到所述的转换因子为
[0025][0026][0027]
tb＝δtdꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0028]
上式中，和分别为热源温度信号的一阶、二阶时间导数，tb和td分别为导数峰值时间及峰值时间，δ为转换因子，c为对流换热系数，d
l
为热扩散率；
[0029]
针对纠偏后的旅行时数据集，开展蒙特卡洛分析确定特定不同旅行时均值及其标准差，包括t-5％、t-25％、t-50％、t-75％及t-100％作为旅行时数据集，依据多组数据集依次开展迭代。
[0030]
优选的，步骤(5)中，将热源温度示踪运移过程视为一系列溶质浓度的波前扩展形式，假设初始慢度参数场为均质分布，根据研究域规模及反演精度，建立模型域并进行网格划分；基于旅行时数据集和联合迭代重建算法依次开展迭代进行介质特征重构，结合不同单元射线路径长度，确定各个单元慢度校正增量，对所有单元进行逐次加权校正，以匹配所测射线投影，逐次迭代直至旅行时差异小于预定阈值，再将慢度场转换为渗透参数分布，具体转换为：
[0031][0032]
其中，r为热阻系数，φ为孔隙率，i为局部水力梯度，k为渗透系数；结合步骤(6)中
得到的多组旅行时数据集，在迭代次数相同的情况下，开展参数敏感性分析，通过对比观测值与重建值误差，选取最佳重建旅行时得到的参数分布。
[0033]
优选的，步骤(5)中，在迭代次数相同的情况下，开展参数敏感性分析，通过对比观测值与重建值误差，选取最佳重建旅行时得到的参数分布。
[0034]
优选的，步骤(6)中，根据示踪旅行时迭代获得的参数分布进行水文地质参数统计参数的预估，包括均值及方差，随后将模型域的参数均设定为满足高斯分布的非均质场，或采用示踪旅行时迭代获得的参数分布作为初始场，结合钻孔取芯及物探手段，估计参数场统计特性；同时，结合钻孔取芯情况得到的地层情况，确定相关尺度。
[0035]
优选的，步骤(7)中，基于同时连续线性估计算法，执行水头数据的地质统计反演分析，结合初始参数分布，求解条件有效水头场，确定水文地质参数均值、方差，以及其参数与水头响应的无条件协方差及交叉协方差、相关贡献系数信息；
[0036]
随后，确定第一次水文地质参数估计，得到首次迭代的参数扰动及各组抽水的观测水头与模拟水头差异，更新参数条件均值及其协方差；
[0037]
基于首次迭代后的参数分布更新条件有效水头场，更新水头扰动及相关协方差，进一步计算确定水头扰动在不同时刻下不同抽水过程，对各个单元位置处各参数的贡献度系数，重复迭代改善前一次参数估计，直至更新后的参数条件协方差或水头误差满足阈值设定。
[0038]
本发明的有益效果为：本方法实施基础为地下工程中施作的超前地质钻孔，无需重新重复布孔作业，对不同地质情况下的地下工程具有普适性；同时采用的试验监测布置方案紧密结合超前地质钻探的钻孔长度、深度、倾角、位置及数量条件，并与隧道(洞)、巷道洞径有关，并可根据钻探揭露的前方含导水断层及充水溶洞等地质构造位置及分布情况动态调整具体监测方案；采用的热源示踪法因其监测便捷、无污染具备良好的应用前景，并且采用平滑低通滤波与自适应卡尔曼滤波相结合的滤波方式，能够实现状态预报误差的协方差矩阵及增益矩阵的动态调整，具有跟踪突变状态能力，有效避免滤波发散；通过开展水文地质数据的联合层析反演，实现了充分融合利用示踪试验及抽注水试验等水文地质试验获取的常规数据信息，显著降低了反问题的多解性，解译了多源水文地质数据中包含的致灾构造富水性的非冗余信息，为提前准确探明工程前方赋存的致灾构造提供理论依据。
附图说明
[0039]
图1为本发明地下工程热源示踪及水力试验监测布置整体示意图。
[0040]
图2为本发明地下工程热源示踪及水力试验监测布置分解示意图。
具体实施方式
[0041]
一种基于热源示踪与水力联合层析反演的致灾构造超前预报方法，包括如下步骤：
[0042]
(1)基于具体工程研究区的地质资料、物探及钻探资料，针对富水构造破碎带、岩溶发育区段及短中距离预报中的物探异常区域，结合超前地质探测施作的超前钻孔长度、深度及倾角情况，圈定掌子面前方最长钻孔的1.2-1.5倍深度范围为监测区域，常为30至100米范围。随后，结合地下工程超前钻探布置位置、数量及隧道(洞)、巷道洞径等，确定热
源示踪试验及分层补源水力试验监测方案，包括监测组数及每组数量、监测间距。
[0043]
在隧道、隧洞及巷道等地下工程中，结合超前钻孔施工设计，将监测点布置于工程掌子面前方三维空间，图1和图2为示例方案之一，可结合超前钻孔的数量确定监测组数，根据隧道、隧洞及巷道等的洞径尺寸确定每组的具体监测点数，如工程所需探测精度较高且研究区岩溶发育程度高、富水性强，可采用的试验监测间距为0.5-1倍洞径，否则最大可采用2-3倍洞径作为监测间距。另外，如超前钻孔已揭示潜在富水致灾构造，则应结合潜在致灾构造区域附近3倍洞径范围内的多个钻孔，尽可能加密监测点，其监测间距可考虑采用0.3-0.5倍钻孔最小间距，重新试验开展精细探测，以确定致灾构造规模及特征；
[0044]
(2)根据已确定的监测区域及试验监测方案，实施热源示踪试验及抽注水水力试验，监测记录示踪穿透曲线及水头曲线。
[0045]
热源示踪试验可采用单点投放、多点监测的形式，热水水源投放宜选取在掌子面中心钻孔的监测点内，所用热源温度应高于环境背景温度40-60摄氏度，且注入时长为监测总时长的一半后，距离注入端口最近的监测端口温度宜上升20摄氏度以上。停止注入后，在保证前序试验后监测孔温度恢复至初始背景温度后方可进行后续试验，一般温度变化低于5％视为恢复。
[0046]
水力试验采用单点或多点抽注，同时多点持续监测，当监测点数少于十个时，开展抽注水实验不少于五次；当监测点数多于十个，则可适当减少试验次数。由于试验早期水头变化快，在试验早期应采用短监测时间间隔，具体间隔应短于抽注水试验总时间的1％，尽可能多的收集早期水头数据，并在水头达到稳定后停止试验，可考虑选取设计试验总时间的四分之一时长作为抽注水早期阶段。
[0047]
(3)开展曲线平滑消噪滤波处理，提取示踪旅行时及水头降深等数据。研究设计窗口加权随方差变化的平滑低通滤波函数，开展预消噪处理，随后结合自适应卡尔曼滤波方法开展深度滤波消噪；
[0048]
针对平滑低通滤波，就温度示踪穿透曲线可考虑采用高斯滤波函数进行消噪，为了在保证消除噪音的同时，尽可能多的保留曲线特征，可根据窗口内数据方差与整体数据集总方差的相对大小确定高斯函数方差，当窗口内数据方差大于总方差时，高斯函数方差选取为1.5，反之取为0.5。进而结合高斯函数方差，确定窗口内各数据的权重向量，根据权重及实测温度数据获取预消噪结果；另一方面，就水头曲线可考虑采用五阶或六阶多项式函数进行拟合消噪，且所用多项式阶次需小于窗口长度，根据窗口内数据方差与整体数据集总方差的占比确定加权系数，进而对给定高阶多项式开展最小二乘拟合，逐步移动平均得到预消噪水头曲线。
[0049]
针对自适应卡尔曼滤波，假设系统噪声和观测噪声互不相关，结合滤波本身开展系统噪声方差及量测噪声方差的估计和修正，引入时变衰减因子l
k+1
依据噪声动态变化，修正常规卡尔曼滤波中的预测协方差矩阵u
k+1,k
，实时调整状态预报误差的协方差矩阵及增益矩阵，保证增益矩阵能够依据模型的不确定性与外界噪声的变化来自动进行调整，令滤波器具有跟踪突变状态的能力，避免滤波发散。
[0050]
[0051][0052]
uk为滤波协方差矩阵，为状态转移矩阵，qk为量测噪声的协方差矩阵，i为单位矩阵，trace[]为求矩阵迹，zk为新息序列，rk为系统噪声协方差矩阵，v
k-1
为残差协方差阵。
[0053]
针对消噪后的温度示踪穿透曲线，进行蒙特卡洛分析提取峰值旅行时数据。针对消噪后的水头曲线提取降深数据，一般将开始试验后四分之一总时长视为抽水试验早期，停止抽水后至水头完全恢复视为抽水试验晚期，其余时间为中期。结合早、中、晚期数据，在每组水头曲线选取不少于五个降深数据，其中在早、晚期时间内至少各选取两个数据。
[0054]
(4)考虑到温度传播过程中的热扩散效应及弥散效应，基于温度穿透时间曲线的一阶导数变化最速时刻确定延迟因子，实现示踪旅行时纠偏。为规避早期旅行时选取的不确定性，针对纠偏后的旅行时数据集，进一步开展早期旅行时诊断分析。
[0055]
根据热对流扩散控制方程，确定温度函数导数峰值时间与峰值旅行时的定量关系，进一步可得到所述的转换因子为
[0056][0057][0058]
tb＝δtdꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0059]
上式中，和分别为热源温度信号的一阶、二阶时间导数，tb和td分别为导数峰值时间及峰值时间，δ为转换因子，c为对流换热系数，d
l
为热扩散率。
[0060]
针对纠偏后的旅行时数据集，开展蒙特卡洛分析确定特定不同旅行时均值及其标准差，包括t-5％、t-25％、t-50％、t-75％及t-100％作为旅行时数据集，依据多组数据集依次开展迭代。
[0061]
(5)基于旅行时数据集和联合迭代重建算法，开展迭代进行介质特征重构，选取最佳重建旅行时及其估计参数分布。
[0062]
具体做法为，将热源温度示踪运移过程视为一系列溶质浓度的波前扩展形式，假设初始慢度参数场为均质分布，根据研究域规模及反演精度，建立模型域并进行网格划分。
[0063]
随后，基于温度旅行时数据集和联合迭代重建算法依次开展迭代进行介质特征重构，结合不同单元射线路径长度，确定各个单元慢度校正增量，对所有单元进行逐次加权校正，以匹配所测射线投影，逐次迭代直至旅行时差异小于预定阈值，再将慢度场转换为渗透参数分布，具体转换可参考：
[0064][0065]
其中，r为热阻系数，φ为孔隙率，i为局部水力梯度，k为渗透系数。在迭代次数相同的情况下，开展参数敏感性分析，通过对比观测值与重建值误差，选取最佳重建旅行时得到的参数分布。
[0066]
(6)结合基于旅行时得到的估计水文地质参数分布，确定含水层参数的初始估计均值、方差及相关尺度等信息。
[0067]
在一些实施案例中，根据示踪旅行时迭代获得的参数分布进行水文地质参数统计参数的预估，包括均值及方差等，随后将模型域的参数均设定为满足高斯分布的非均质场；或者可采用示踪旅行时迭代获得的参数分布作为初始场，并结合钻孔取芯及物探等手段，估计参数场统计特性。此外，确定相关尺度应充分考虑钻孔取芯情况得到的地层情况。
[0068]
(7)结合同时连续线性估计算法，开展水文地质参数反演分析，逐次迭代直至满足阈值设定，随后计算确定渗透系数、储水系数及导压系数等参数分布，结合高渗区圈定富水致灾构造位置及分布。
[0069]
具体的，基于同时连续线性估计算法，执行水头数据的地质统计反演分析，结合初始参数分布，求解条件有效水头场，确定水文地质参数均值、方差，以及其参数与水头响应的无条件协方差及交叉协方差、相关贡献系数等信息。
[0070]
随后，确定第一次水文地质参数估计，得到首次迭代的参数扰动及各组抽水的观测水头与模拟水头差异，更新参数条件均值及其协方差。基于首次迭代后的参数分布更新条件有效水头场，更新水头扰动及相关协方差，进一步计算确定水头扰动在不同时刻下不同抽水过程，对各个单元位置处各参数的贡献度系数，重复迭代改善前一次参数估计，直至更新后的参数条件协方差或水头误差满足阈值设定。