一种深埋隧道围岩介电特性反演与松动圈范围识别方法与流程

1.本发明涉及深埋隧道及地下工程领域，尤其是针对围岩松动圈识别的一种高精确的研究方法，具体为一种基于探地雷达的深埋隧道围岩介电特性反演与松动圈范围识别方法。


背景技术：

2.进入21世纪以来，我国隧道及地下工程得到了前所未有的发展，目前我国已经是世界上隧道及地下工程规模最大、数量最多、修建技术发展最快的国家。截至2020年底，我国投入运营的铁路隧道共16798座，总长约1963万米。从1980年至2020年的40年间，中国共建成隧道12412座，总长约17621km(占中国铁路隧道总长度的90％)。特别是近15年来，中国铁路隧道发展极为迅速，共建成铁路隧道9270座，总长约15321km(占中国铁路隧道总长度的78％)，由此可见，我国已经进入了隧道工程建设的高峰期，与此同时，深埋隧道逐渐增加，如：峨汉公路隧道最大埋深1944m；川藏铁路隧道最大埋深2600m。此外，在水电隧洞领域、矿山井巷领域以及能源储库领域也有大批深部工程建设，未来，我国深部地下工程建设将区域常态。
3.深埋隧道工程赋存地质环境极端复杂，高地应力、高地温、高渗透压等问题突出，与浅部地层相比，地质灾害发生频率更高、致灾机理更复杂。深埋隧道爆破开挖时不可避免的会对围岩产生损伤扰动，导致岩石原本稳定的应力状态受到不平衡影响，进入不稳定状态，岩石应力由三向应力状态转变成二向应力状态，一定范围内围岩会发生破裂，形成一个逐步趋于稳定的破碎带，此范围内围岩状态处于屈服状态，称为塑性松动圈，形成松动圈。
4.根据经验，围岩松动圈范围的厚度，是判断围岩是否稳定的依据，也是支护的关键依据。复杂条件下深埋隧道围岩松动圈范围更加难以确定。常用的测试方法中：声波法仪器操作简单但抗干扰性较差；多点位移法操作简单、测点少但周期长，测量精度低；地震波法测试精度高但安装困难，测试成本高；电阻率法布点方便、测试简单但对仪器精度要求高。


技术实现要素：

5.发明目的：
6.本发明提出一种深埋隧道围岩介电特性反演与松动圈范围识别方法，其目的在于解决常规方法对深埋隧道爆破后松动圈范围测定遇到的测量精度不高、测量数据不可靠以及操作复杂等问题。基于地质雷达正演模拟与围岩介电特性反演理论，为松动圈识别提供了一种新方法。
7.本发明采用的技术方案为：
8.一种深埋隧道围岩介电特性反演与松动圈范围识别方法，其特征在于：
9.步骤一：在深埋隧道开挖前，进行雷达实地探测，确定区域的大小、介质类型以及地下工况，并输出探测数据和雷达扫描图；雷达实地探测后进行正演模拟，确定区域介电常数ε
前
、电导率σ
前
和磁导率μ
前
，并得到区域的扫描图，通过开挖前波速模型求出岩石原始状态
下雷达波速v
前
；
10.步骤二：在深埋隧道开挖后，进行雷达实地探测和正演模拟输出松动圈材料的磁导率μ，通过介电特性反演输出介电常数ε和电导率σ，将松动圈材料的磁导率μ、介电常数ε代入开挖后波速模型求出深埋隧道开挖后波速v；
11.步骤三：根据波速v计算松动圈厚度；通过波速v与雷达波速v
前
比值的平方作为对围岩松动圈的完整程度进行定量分析，根据岩体完整性指数划分情况表得出松动圈等级。
12.进一步的，步骤一中开挖前波速模型为：
[0013][0014]
式中：ε
前
为隧道开挖前材料的介电常数；ε0为真空的介电常数，参考值ε0＝8.854
×
10-9
f/m；μ
前
为隧道开挖前材料的磁导率；μ0为真空的磁导率，参考值μ0＝4π
×
10-7
h/m；c为光速，取3
×
108m/s；v
前
为隧道开挖前电磁波在介质中的波速。
[0015]
进一步的，步骤二中开挖后波速模型为：
[0016][0017]
式中：ε为松动圈材料的介电常数；μ为松动圈材料的磁导率；ε0为真空的介电常数，参考值ε0＝8.854
×
10-9
f/m；μ0为真空的磁导率，参考值μ0＝4π
×
10-7
h/m；c为光速，取3
×
108m/s；v为松动圈的电磁波在介质中的波速。
[0018]
进一步的，步骤二中介电特性反演步骤为：
[0019]
2.3.1：通过正演模拟得到数据与实地探测数据进行匹配，基于数据残差建立目标函数；
[0020]
2.3.2：根据直线传播理论选择初始模型；
[0021]
2.3.3：使用初始模型计算正演波场d(ε,σ)；
[0022]
2.3.4：采用局部优化算法，计算迭代方向
[0023]
2.3.5：计算迭代梯度；
[0024]
2.3.6：计算迭代步长ζ
ε
及ζ
σ
；
[0025]
2.3.7：更新模型参数；
[0026]
2.3.8：设ε＝ε
upd
，σ＝σ
upd
重复步骤2.3.3至步骤2.3.7，直至收敛或达到指定的迭代次数，输出介电常数ε和电导率σ。
[0027]
进一步的，步骤三中松动圈厚度求解公式为：
[0028]
h＝v
×
t
[0029]
式中：h为松动圈厚度；t为雷达波单程走时。
[0030]
有益效果：
[0031]
本发明基于雷达探测结合松动圈形成原理，提出一种计算松动圈厚度的方法，判断出松动圈围岩等级，根据介电特性对开挖后的松动圈进行了精细刻画，解决了松动圈测量困难、测量精度不高和数据不可靠等的问题。
附图说明
[0032]
图1是流程总图；
[0033]
图2是探地雷达数据处理流程图；
[0034]
图3是探地雷达探测方式中心点法；
[0035]
图4是正演模拟流程图；
[0036]
图5是裂缝正演灰度图；
[0037]
图6是空洞正演灰度图；
[0038]
图7是介电特性反演流程图；
[0039]
图8是松动圈形状图；
[0040]
图9是围岩等级图。
具体实施方式
[0041]
以下结合说明书附图更详细的说明本发明。
[0042]
本发明在保证支护结构合理的情况下，采用应用范围广、精度高、分辨率高、抗干扰能力强的探地雷达法。探地雷达向地下发射高频电磁波，当发射的高频电磁波遇到不同介电特性的介质时会发生反射、折射以及衍射，反射的电磁波被接收天线所接收，经过信号处理和数据分析形成雷达波的反射图像，根据反射波的波形特征和振幅等情况可以推测出地下介质的分布情况及属性。但只能根据实际雷达剖面粗略地推断出地下具有电性差异的目标层界面、空间位置分布以及目标体形态、尺寸等信息，却无法提供对介质介电特性(介电常数、电导率)的定量估计和对近地表结构的精细刻画。
[0043]
为了弥补探地雷达自身不足，本发明将探地雷达与介电特性反演结合。反演利用了波场的运动学和动力学特征(振幅和相位)，具有更高的反演分辨率和建模精度，提供了一个可信的波速分布，具有强非线性的映射特点。反演以正演模拟为基础，通过正演模拟得到数据与实地探测数据进行匹配，基于数据残差建立目标函数，利用数据反演地下介质地球物理参数的方法，求出精确的介电常数和电导率，达到对地下复杂结构精确刻画的目的。
[0044]
本发明中，在深埋隧道开挖后，进行探地雷达和正演模拟；通过正演模拟得到数据与实地探测的数据进行匹配，基于数据残差建立目标函数，利用数据反演出介电特性(介电常数、磁导率)，介电特性可以更好的反映出地质状况，根据电磁波在介质中的传播规律，计算出电磁波的传播速度。首次提出松动圈的计算方法：即根据波速与单程走时求出松动圈厚度；根据隧道开挖前后波速比的平方，确定出岩石破坏程度，根据岩石完整性划分情况表，确定出松动圈等级。基于松动圈厚度和松动圈等级对围岩进行更合理的支护设计。
[0045]
如图1所示，本发明采用的技术方案分为三个步骤：
[0046]
步骤一：在深埋隧道开挖前，进行雷达实地探测，确定区域的大小、介质类型以及地下工况，并输出探测数据和雷达扫描图；雷达实地探测后进行正演模拟，确定区域介电常数ε
前
、电导率σ
前
和磁导率μ
前
，并得到区域的扫描图，通过开挖前波速模型求出岩石原始状态下雷达波速v
前
；
[0047]
在深埋隧道开挖前，进行雷达实地探测进行如下操作：
[0048]
1.1.1：雷达探测前，需确认探地雷达探测深度、所用波源等自身功能条件。
[0049]
1.1.2：雷达探测，确认地质条件，结合探地雷达选择合适的雷达波；根据电磁波的
波长，选择合适的探测深度。
[0050]
1.1.3：确定探测方式：采用现有的中心点法(图2)，保持发射和接收天线的中心点位置保持不变的情况下，逐渐变大两个天线离中心点的距离，从而记录整个测线上发射雷达波数据。
[0051]
1.1.4：确定发射天线、接收天线位置，记录收发时间。
[0052]
探地雷达经发射天线接到控制单元发出的指令后会产生特定频率的电信号，电信号经发射天线转换为电磁波信号并发射出来，而接收天线则负责接收目标反射的电磁波信号，并转换成电信号进行存储。反射信号中包含了介电特性差异界面的许多重要信息，如目标的位置、尺寸、性质等等。
[0053]
1.1.5：因探测过程中，受到环境影响以及探地雷达自身缺陷，会对原始数据进行干扰，造成雷达图像不清晰，需对探测所得原始数据进行处理(如图3)，剔除无用信息。
[0054]
数据处理包括：数据合并、废道剔除、道内均衡、背景消除。
[0055]
①
数据合并和废道剔除：采取一种简易处理方法，就是直接用与废道相邻的数据代替废道数据；
[0056]
②
道内均衡：将不同道的能量强的波按照相应的比例进行压缩，而相应的按对应比例将弱波的进行放大，以使得两种波的振幅在特定范畴内进行波动，简言之，即把各道数据在对应时间内的振幅和一个确定的值进行相乘即可；
[0057]
③
背景消除：首先对需要处理部分数据求平均值，这些平均值代表的就是有规则的水平干扰信号，再把每道数据减去这个均值背景干扰信号，这样水平干扰信号就会被减弱或者去除，使得没有规律的有用信息的反射信号得到了相对加强。
[0058]
经过对原始数据处理后得到较为清晰准确的探测数据和扫描图，记录保存。
[0059]
实地探测后，进行正演模拟，计算雷达波速v
前
。
[0060]
正演模拟前，根据探测实况，确定探测深度、波源、探测层数以及介质层厚度，根据实地探测大小设置正演模型。
[0061]
根据探测介质，查找介质的介电常数、电导率、相对渗透率和磁损耗，并设置收发天线位置，确定时窗。
[0062]
如图4，正演模拟步骤如下：
[0063]
1.2.1：基于时域有限差分法，并且确定正演模型所需的吸收边界、解的稳定性和激励源等参数。
[0064]
所述正演模拟基于时域有限差分法(fdtd)，以有限差分和微分形式的maxwell(麦克斯韦)旋度方程为基础，对电磁场进行时域数值分析。其表达式：
[0065][0066]
电磁场的本构方程如下：
[0067]
d＝εe,b＝μh,j＝σe,jm＝ρh
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0068]
式中：e为电场强度(v/m)；h为磁场强度(a/m)；d为电通量密度(c/m)；b为磁通量密度(wb/m2)；jm为磁流密度(v/m2)；j为电流密度(a/m2)；ε为介电常熟(f/m)；μ为磁导率(h/m)，文中磁导率被认为是不变的；ρ为磁损耗(ω/m)；t为时间(s)。
[0069]
1.2.2：确定吸收边界、解的稳定性和激励源，保证正演模拟边界处的波不受有限空间影响，并保证解的计算收敛。完全匹配层吸收边界公式、解的稳定性限制公式如下：
[0070]
完全匹配层吸收边界进行设置，计算公式：
[0071][0072]
σ＝ρ＝0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0073]
解的稳定性公式：
[0074][0075]
1.2.3：通过.in文件编写命令，建立仿真模型，在gprmax仿真软件运行，确定区域介电常数、电导率和磁导率等的正演数据，并得到区域的扫描图。
[0076]
仿真模型包括六种命令：
[0077]
基本命令：包括模型大小、网格和时窗；
[0078]
一般命令：对基本命令进一步说明；
[0079]
介质命令：引入不同材料，确定材料参数(即介电常数、电导率、相对渗透率和磁损耗)；
[0080]
对象构造：构建具有不同本构参数的几何形状；
[0081]
激励和输出命令：模型中放置源点和输出点；
[0082]
pml命令：吸收边界条件的高级定制和优化。
[0083]
正演结束后，使用matlab查看数据，并记录。
[0084]
正演结束，为求出开挖前波速v
前
则将正演的扫面图，与实地雷达探测的扫描图进行对比，两者基本一致时，使用不变的磁导率μ
前
和正演的介电常数ε求波速v
前
，开挖前波速模型如下：
[0085][0086]
式中：ε
前
为隧道开挖前材料的介电常数；ε0为真空的介电常数，参考值ε0＝8.854
×
10-9
f/m；μ
前
为隧道开挖前材料的磁导率；μ0为真空的磁导率，参考值μ0＝4π
×
10-7
h/m；c为光速，取3
×
108m/s；v
前
为隧道开挖前电磁波在介质中的波速(m/s)。
[0087]
步骤二：在深埋隧道开挖后，进行雷达实地探测和正演模拟输出松动圈材料的磁导率μ，通过介电特性反演输出介电常数ε和电导率σ；由于正演过程中认为磁导率是不变的，则磁导率无需计算。将松动圈材料的磁导率μ、介电常数ε代入开挖后波速模型求出深埋隧道开挖后波速v；
[0088]
隧道开挖常用爆破的方式，爆炸冲击波的威力很大，巷道的采掘过程不是一次爆破能完成的，数次爆破后，对岩体造成很大的损伤破坏。
[0089]
研究多次爆破振动对巷道围岩累计损伤特性更能揭示岩体的破坏过程和动态失稳。从宏观和微观上进行累计损伤弱化效应的演变。
[0090]
所述宏观方面：一次爆破不会对中远区岩体造成破坏，但是多次爆破后中远区岩体表面也会出现裂纹，部分严重区域的岩体性能己完全破坏甚至有岩爆发生。
[0091]
所述微观方面：岩体本身内部会有或多或少的裂纹，每一次的爆破都会对内部裂纹造成不同程度的影响，爆炸应力波作用多次后，小裂纹逐渐扩展成大裂纹，裂纹之间也会相互贯通，最终形成大的主裂缝。
[0092]
巷道爆破开挖时会对围岩产生损伤扰动，引起应力重分布形成松动圈。
[0093]
松动圈在爆破作用下，其发展分为四个阶段，分别为应力调整阶段、松动圈产生阶段、松动圈快速发展阶段和松动圈稳定阶段。
[0094]
松动圈探测：探测区域与开挖前保持一致，保证探测中雷达采样是连续记录。
[0095]
2.1：隧道开挖后松动圈进行雷达探测，同步骤一中的雷达实地探测步骤，输出探测数据。
[0096]
2.2：隧道开挖后进行正演模拟，其步骤同步骤一中正演模拟步骤。正演模拟直观反映出在裂缝和空洞的工况下雷达波的反射状态；模型正演：尺寸：5m
×
3m矩形，激励源：800mhz，三层介质：黏土、砂岩和岩石，在砂岩层出现裂缝和空洞，其裂缝正演灰度图如图5，其空洞正演灰度图如图6。建立裂缝和空洞的工况，正演灰度图可以更直观的反映出雷达波的反射状态，同时也可以根据反射特征更好的判断地质状态。
[0097]
2.3：建立目标函数进行介电特性反演(流程如图7)，介电特性反演步骤如下：
[0098]
2.3.1：介电特性反演主要是寻找介电常数ε与电导率σ的空间分布。通过正演模拟得到数据与实地探测数据进行匹配，基于数据残差建立目标函数，目标函数表达式：
[0099][0100]
其中，d
obs
为实地探测的数据；d(ε,σ)是从正演得到数据记录；t为转置符号；s，n，τ分别表示震源、接收点和观测时间。该式在源s，接收点n以及观测时间τ三个方向上求和。
[0101]
2.3.2：根据直线传播理论选择初始模型；
[0102]
初始模型理论：本质上是数据拟合，即根据模拟数据与实地探测数据之间的误差，去寻找与探测数据最匹配的模型，拟合出和探测数据最符合的模型参数。由于细节也可通过反演进行拟合，则该方法提高了准确性和分辨率。
[0103]
2.3.3：使用初始模型计算正演波场d(ε,σ)；
[0104]
2.3.4：采用局部优化算法，计算迭代方向
[0105]
2.3.5：计算迭代梯度；
[0106]
方程(7)的关键是求解电场，由于电场是关于介电常数和电导率的方程，确定梯度，需要对这些介电特性的电场进行线性化处理，得到扰动目标函数的一阶近似，公式如下：
[0107][0108]
对公式(8)进行展开，并线性化处理，引入接受点位置xn，得到每一个发射源、每一
个接收位置、每一个接收时刻的全尺寸空间的介电常数和电导率的梯度：得到梯度公式如下：
[0109][0110]
式中：t
s,n,τ
为多点源，即每一个接收点位置都有一个源，在时间上为逆时传播。
[0111]
计算目标函数的每个发射源、接收点和观测时间的各自贡献，然后累加这些梯度项，来估计总梯度
[0112]
2.3.6：计算迭代步长ζ
ε
及ζ
σ
；
[0113]
迭代步长中介电常数和电导率遵循以下公式：
[0114][0115][0116]
计算迭代步长ζ
ε
和ζ
σ
，为了使梯度方向函数最小化，引入扰动因子η
ε
和η
σ
，扰动因子的选取，即要保证足够小又保证扰动模型处于线性范围内。公式为：
[0117][0118][0119]
2.3.7：更新模型参数。
[0120]
采用共轭梯度算法，更新模型参数，方程式如下：
[0121]mω+1
＝m+ζδm
ω
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0122]
式中：m
ω
和m
ω+1
分别为第ω次和第ω+1次迭代的模型参数，ζ为更新步长，δm
ω
为模型的更新方向。
[0123]
更新方向δm
ω
，采用共轭梯度算法，表达式如下：
[0124]
δm
ω
＝-n(m
ω
)+β
ω
δm
ω-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0125]
式中：δm
ω
与δm
ω-1
互为共轭；n：表示梯度；β
ω
可通过fletcher-reeves法求取，其表达式：
[0126][0127]
共轭梯度算法具有收敛快、所需存储量小，不需要引入任何外来参数等优点。利用当前次迭代和上一次迭代所求取的梯度，构造一组共轭方向作为新的更新方向，并沿这组方向进行搜索，最终求取目标函数的极小点，从而实现快速迭代求取最优解的目标。共轭梯度算法步骤如下：
[0128]
a：设置初始模型m0和终止迭代精度0≤ε≤1，计算初始迭代梯度且令
迭代次数ω＝0。
[0129]
b：当||n
ω
||≤ε，迭代运算停止，输出最终结果m
ω
；否则，继续迭代更新。
[0130]
c：计算更新方向δm
ω
：
[0131][0132]
所述β
ω
由(式16)确定。
[0133]
d：寻找一个沿负梯度方向的最小目标函数l，寻找迭代步长ζ。
[0134]
e：求取更新后的模型参数：m
ω+1
＝m
ω
+ζδm
ω
，并求取当前迭代的梯度
[0135]
f：令ω＝ω+1，转回第二步。
[0136]
在matlab中需编写梯度、优化步长、共轭梯度算法迭代的程序，实现介电特性反演。
[0137]
2.3.8：设ε＝ε
upd
，σ＝σ
upd
重复步骤2.3.3至步骤2.3.7，直至收敛或达到指定的迭代次数，输出介电常数和电导率。
[0138]
反演结束后，输出介电常数ε和电导率σ，与正演中不变的磁导率μ求深埋隧道开挖后的波速v，开挖后波速模型为：
[0139][0140]
式中：ε为松动圈材料的介电常数；μ为松动圈材料的磁导率；ε0为真空的介电常数，参考值ε0＝8.854
×
10-9
f/m；μ0为真空的磁导率，参考值μ0＝4π
×
10-7
h/m；c为光速，取3
×
108m/s；v为松动圈的电磁波在介质中的波速(m/s)。
[0141]
步骤三：确定松动圈范围：提出一种松动圈计算方法，根据波速v计算松动圈厚度；通过波速v与雷达波速v
前
比值的平方作为对围岩松动圈的完整程度进行定量分析，根据岩体完整性指数划分情况表得出松动圈等级。
[0142]
深埋隧道爆破开挖后，引起应力重分布形成松动圈。
[0143]
松动圈形状最终呈圆角菱形，第一阶段沿垂直岩层方向发展，其形状近似椭圆；第二阶段沿平行于岩层方向快速发展，这是由于松动圈形成后围岩的应力调整所致。松动圈形状整体呈现出正交各向异性特点，以垂直层理方向大于平行层理方向为特征的圆角菱形(图8)。
[0144]
围岩松动圈是在围岩应力重新分布过程中，在围岩应力超过围岩强度后形成的，其力学特性表现为应力和强度均降低，宏观上岩体微裂隙增多，变形增大，岩体内聚力、摩擦角降低，其波速减小。
[0145]
3.1.1：根据松动圈波速v计算松动圈厚度，其公式如下：
[0146]
h＝v
×
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(19)
[0147]
式中：h为松动圈厚度(m)；t为雷达波单程走时(s)。
[0148]
单程走时，根据隧道开挖后实地探测记录的收发天线之间总时间的一半。经过介电特性反演获得更准确的波速与实地探测时间确定出松动圈厚度，对松动圈的判断，将传统的速度传感器测试的峰值作为松动圈判断的依据转化为把松动圈进行量化，更直观，更
具体的体现松动圈的厚度，为支护方案提供更准确的依据。
[0149]
3.1.2：确定松动圈等级：
[0150]
以波速v与雷达波速v
前
比值的平方作为对围岩松动圈的完整程度进行定量分析，根据岩体完整性指数划分情况表得出松动圈等级。
[0151]
围岩松动圈本质是一个“松弛破碎带”，由于岩体开挖应力场中附加应力超过了围岩强度，造成围岩破碎。结合松动圈形成原理，引入完整性系数指标来确定围岩松动圈的范围。
[0152]
所述松动圈本质是“松弛破碎带”，可近似的认为松动圈完整性系数等于破碎岩石完整性系数，其公式：
[0153][0154]
式中：rv为破碎岩石完整性系数；lv为松动圈完整性系数。
[0155]
根据波速在介质中的传播规律，通过反演出更准确的介电常数和电导率，提出一种计算松动圈的方法：根据波速与单程走时求出松动圈厚度；根据隧道开挖前后波速比的平方，确定出岩石破坏程度，根据岩石完整性划分情况表，确定出松动圈等级。通过公式计算，得出松动圈厚度数值，更直观，更准确的确定松动圈范围，根据岩石完整性系数，确定松动圈等级，根据介电特性可以更精确的刻画松动圈状况，在弥补松动圈探测不准，支护不合理缺点的同时为松动圈支护方案提供理论依据。
[0156]
岩体完整性指数划分情况表：
[0157]
围岩等级iiiiiiivv完整系数＞0.750.75～0.550.55～0.350.35～0.15＜0.15
[0158]
根据岩体完整性指数划分情况表确定松动圈等级(图9)，根据松动圈等级进行对应的围岩加固措施和初期支护设计。