基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件的方程改进方法与流程

1.本发明涉及基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件的方程改进方法，尤其涉及到时间限制以及有限采掘边界长度的电磁波强度限制的解决方法，属于隧道工程技术领域。


背景技术：

2.传统的电磁波波场中，瓦斯隧道掌子面空腔以及异构体的断层结构只要传统的运用矿井直流电法勘探、矿井瞬变电磁法勘探、探底雷达和地震勘探类方法。
3.在本发明中主要依托于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件进行研究，电磁波在地下岩层中传播时，由于各种岩石、矿电性的不同，它们所呈现出来的电磁波波场有一定的差异，且不同的电磁波波场对岩层的识别存在一定的限制，这就需要考虑到在探测过程中电磁波波场对巷道空腔的边界问题，通过改进有限差分方程，可以改进出有边界条件的巷道空腔电磁波探测方程。但是当电磁波在地下岩层传播时，如果存在着含水地段、陷落柱、断层、空洞或其它不均匀的地质构造，电磁波的波场能量就会被其吸收或者完全屏蔽，此时就会产生信号显著性异常，透视异常。
4.对于之前的有限差分相关算法，大多数学者都是基于其稳定性和精确性进行研究，少有考虑到精细网格细分改进差分方程的问题，所以本发明针对算法入手，对有限差分方程进行时间线性插值，以此实现减小误差的效果。此外，在重建的过程中，会出现电磁波之间不是平行束的情况，此时应考虑角度问题，本发明将运用对算法进行改进，以最大程度达到真值，实现对岩层的勘探实测目的。


技术实现要素：

5.本发明提供了基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件的方程改进方法，以用于解决矿井探测中使用的有效差分法中的网格化不够精细、时间长，不能稳定的解出有限采掘边界长度的电磁波波场强度问题。
6.本发明的技术方案是：基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件的方程改进方法，所述方法的具体步骤如下：
7.step1、根据探地雷达下正演的相关知识以及真实场景，运用非均匀网格技术，求解案例当前情况下的有限差分方程；
8.首先，根据背景知识，地质雷达主要通过发射天线m向探测的目标发射低频或高频电磁波，其通过高频电磁波在介质中的反射和折射来实现对地下介质的探测。在探测过程中，若探测到电性差异较大的界面时，会发生折射或投射等现象，之后被接天线n接收返回雷达主体，即地质雷达是通过雷达得到的电磁波信号分析波的形状、强度等特征，从而可以判断出地质体的结构。
9.本发明基于如图4所示的地质雷达探测图，对于隧道建设中出现的空洞进行质量
检测。电阻率、介电常数、磁导率等是影响雷达波在地下介质中传播的电性参数。如果在隧道中出现空洞，空洞内一般为空气，而空气(相对介电常数为1))的介电常数等电性参数与衬砌混凝土(相对介电常数6～10)有很大差异，电磁波在衬砌内传播时，会在空洞上下两个界面产生两次强反射，雷达剖面上会出现双曲线形态的强反射波，其同相轴与相邻道会发生错位。因此能够利用波形特征预计波场的分布情况来识别隧道中的空洞。
10.需要运用数值计算的方法对整个空间进行非均匀网格技术分割，将无限大场的分布空间截断为有限的计算空间，在分割的过程中主要采用的是引入合适的边界条件，使得电磁场扩散至边界时能没有明显的反射，运用数学公式得到每一个分割点的有限差分方程如下：
[0011][0012]
其中σ表示介质电导率；σa为异常体电导率；jm为磁流密度矢量；h
p
，hs为磁场强度；μ为介质磁导率。
[0013]
step2、引入加权平均值的概念，对网格中的各点周围的电导率进行面积加权平均，得到离散化后有限差分算法各个节点的电磁波场强的数学表达式；
[0014]
所述step2中包括的各个节点的电磁波场强的数学表达式如下：
[0015][0016]
其中，n表示当前时刻，n+1表示下一时刻，i+1就是节点(i,j,k)x轴正方向上多一个单位的数学表示，i-1就是节点(i,j,k)x轴负方向上多一个单位的数学表示，j+1就是节点(i,j,k)y轴正方向上多一个单位的数学表示，j-1就是节点(i,j,k)y轴负方向上多一个单位的数学表示，k+1就是节点(i,j,k)z轴正方向上多一个单位的数学表示，k-1就是节点(i,j,k)z轴负方向上多一个单位的数学表示，为可探测区域总场t＝(n+1)δt时的磁场内中心节点的场强扩散方程；是x，y，z方向的平均网度；为x轴方向的网度；为y轴方向的网度；为z轴方向的网度；为节点(i,j,k)周围电导率的面积加权平均值；i(t)为关断电流函数；s为所在节点处双侧封闭曲面面积，μ为介质磁导率；
[0017]
具体的操作步骤为：
[0018]
给有限差分方程的两边在δv内积分，得到
[0019][0020]
接下来，根据高斯公式对上述式子进行计算并整理，有
[0021][0022]
其中δv表示前后两个节点移动形成的体积差；s
x
、sy和sz分别为δv与各个方向轴正割的曲面正面；s-x
、s-y
和s-z
分别为δv与各个方向轴正割的曲面反面，σ表示介质电导率；σa为异常体电导率；jm为磁流密度矢量；h
p
，hs为磁场强度；接下来对各积分得到积分结果后，进行面积加权平均计算，得到：
[0023][0024]
为节点(i,j,k)周围电导率的面积加权平均，即为
[0025][0026]
其中i,j,k分别表示x,y,z方向上的坐标表示；σ
i,j,k
为节点(i,j,k)周围电导率的面积加权值；σ
i+1,j,k
为节点(i+1,j,k)周围电导率的面积加权值，同理，变换下标就求得不同节点的面积加权值；δxi为x轴上的网格距，δyj为y轴上的网格距，δzk为z轴上的网格距；
[0027]
最后运用时间线性插值算法，得到各个节点电磁波场强数学表达式。
[0028]
step3、监测此时矿井磁场中的各点电导率的稳定性，根据极值来判定监测的时间长短；
[0029]
作为本发明的进一步方案，所述step3中包括具体如下：运用局部细分网格的思路找到时间步长与网格之间的关系来得到趋于稳定监测的时间区域，用时间插值公式表示为：
[0030][0031]
其中，h为与全区域中网格不重合的局部网格内部强度，h为与局部网格重合的网格中以及全区域网格中的强度；v为时间步长；m为分的网格数量。
[0032]
用空间插值表示为：
[0033][0034]
[0035][0036]
局部细分网格有限差分算法的具体操作步骤如下：
[0037]
1)在包含细网格区的粗网格上利用步骤step 2中的式子(4)求得
[0038]
2)在局部细网格中给运用式子(5)迭代两次可以得到
[0039]
step4、运用中心差分近似法对step2中的结果计算差分网度以及最大网度值；
[0040]
作为本发明的进一步方案，所述step4中包括：运用欧拉法算出经过欧拉变换之后的差分方程为：
[0041][0042]
其中，令
[0043]
将时间导数运用中心差分近似有：经过差分近似后，得到新的差分方程为：
[0044][0045]
其中，r＝δt/μσδ2，此处的r为差分网度；为节点(i,j,k)周围电导率的面积加权平均；σ表示介质电导率，μ为介质磁导率；n表示当前时刻，n+1表示下一时刻；
[0046]
又因为欧拉法的稳定性条件是-1≤fe≤1，又因为在上述差分方程式子中fe≤1始终满足，故满上前半部分得到此处的即为最大网度。
[0047]
step5、根据巷道内满足的拉普拉斯方程，对整个矿井巷道空间求体积分，运用时间线性插值进行最后的整理，即可以得到矿井巷道内可控制边界的相应所需方程。
[0048]
作为本发明的进一步方案，所述step5具体包括如下：
[0049]
step5.1、根据laplace方程求积分，在矿井巷道中满足的laplace方程为
[0050]
运用有限差分法分割，对空腔巷道内部所示的δv内求解积分，得到
[0051][0052]
step5.2、时间线性插值；
[0053]
运用时间线性插值，运用欧拉法算出当前时刻经过欧拉变换之后的差分方程得到在矿井巷道内运用有限差分方程时相应所需方程为：
[0054][0055]
本发明的有益效果是：对有限差分方程进行改进，大多数学者都是基于其稳定性和精确性进行研究，少有考虑到精细网格细分改进差分方程的问题，所以本发明针对算法入手，对有限差分方程进行时间线性插值，以此实现减小误差的效果。此外，在重建的过程中，会出现电磁波之间不是平行束的情况，此时应考虑角度问题，本发明将运用对算法进行改进，以最大程度达到真值，实现对岩层的勘探实测。
附图说明
[0056]
图1为本发明的流程图；
[0057]
图2为本发明的物理场景图；
[0058]
图3为有限差分网格部分节点模型图；
[0059]
图4为地质雷达探测示意图。
具体实施方式
[0060]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下获得的所有其他实施例，都属于本发明保护的范围。需要说明的是，在不冲突的情况下，本技术中的实施例及实施例中的特征可以相互任意组合。
[0061]
在煤岩安全防控方面，对于基于瓦斯隧道的正演方法研究，只有极少数的隧道学者和专家研究电磁波波场情况和利用有限差分法方程实现边界控制的方式；在瓦斯隧道分层界面正演方面，巷道空腔识别和波场强度分布是重要的研究方面；通常瓦斯隧道中巷道空腔判断存在比较大的问题，网格化划分次数越高越精细，浪费时间越长。
[0062]
为了解决矿井探测中使用的有效差分法中的网格化不够精细、时间长，不能稳定的解出有限采掘边界长度的电磁波波场强度问题，本发明的实施例提供了基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件的方程改进方法。
[0063]
首先结合附图，对本发明的实施例进行说明。
[0064]
如图1-图4所示，基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件的方程改进方法，所述方法的具体步骤如下：
[0065]
step1、根据探地雷达下正演的相关知识以及真实场景，运用非均匀网格技术，求解案例当前情况下的有限差分方程；
[0066]
首先，根据背景知识，地质雷达主要通过发射天线m向探测的目标发射低频或高频电磁波，其通过高频电磁波在介质中的反射和折射来实现对地下介质的探测。在探测过程中，若探测到电性差异较大的界面时，会发生折射或投射等现象，之后被接天线n接收返回雷达主体，即地质雷达是通过雷达得到的电磁波信号分析波的形状、强度等特征，从而可以
判断出地质体的结构。
[0067]
本发明基于上述地质雷达探测图，对于隧道建设中出现的空洞进行质量检测。电阻率、介电常数、磁导率等是影响雷达波在地下介质中传播的电性参数。如果在隧道中出现空洞，空洞内一般为空气，而空气(相对介电常数为1))的介电常数等电性参数与衬砌混凝土(相对介电常数6～10)有很大差异，电磁波在衬砌内传播时，会在空洞上下两个界面产生两次强反射，雷达剖面上会出现双曲线形态的强反射波，其同相轴与相邻道会发生错位。因此能够利用波形特征预计波场的分布情况来识别隧道中的空洞。
[0068]
需要运用数值计算的方法对整个空间进行非均匀网格技术分割，将无限大场的分布空间截断为有限的计算空间，在分割的过程中主要采用的是引入合适的边界条件，使得电磁场扩散至边界时能没有明显的反射，运用数学公式得到每一个分割点的有限差分方程如下：
[0069][0070]
其中σ表示介质电导率；σa为异常体电导率；jm为磁流密度矢量；h
p
，hs为磁场强度；μ为介质磁导率。
[0071]
step2、引入加权平均值的概念，对网格中的各点周围的电导率进行面积加权平均，得到离散化后有限差分算法各个节点的电磁波场强的数学表达式；
[0072]
所述step2中包括的各个节点的电磁波场强的数学表达式如下：
[0073][0074]
其中，n表示当前时刻，n+1表示下一时刻，i+1就是节点(i,j,k)x轴正方向上多一个单位的数学表示，i-1就是节点(i,j,k)x轴负方向上多一个单位的数学表示，j+1就是节点(i,j,k)y轴正方向上多一个单位的数学表示，j-1就是节点(i,j,k)y轴负方向上多一个单位的数学表示，k+1就是节点(i,j,k)z轴正方向上多一个单位的数学表示，k-1就是节点(i,j,k)z轴负方向上多一个单位的数学表示，为可探测区域总场t＝(n+1)δt时的磁场内中心节点的场强扩散方程；是x，y，z方向的平均网度；为x轴方向的网度；为y轴方向的网度；为z轴方向的网度；为节点(i,j,k)周围电导率的面积加权平均值；i(t)为关断电流函数；s为所在节点处双侧封闭曲面面积，μ为介质磁导率；
[0075]
如图3所示，图中的网格单元，(i,j,k)等7个节点构成的平面可以将该网格单元划分成8个小单元，在每个小单元内都假设电导率都是均匀的，例如在单元oabcdefg内的电导率为σ
i+1,j+1,k
。
[0076]
具体的操作步骤为：
[0077]
给有限差分方程的两边在δv内积分，得到
[0078][0079]
接下来，根据高斯公式对上述式子进行计算并整理，有
[0080][0081]
其中δv表示前后两个节点移动形成的体积差；s
x
、sy和sz分别为δv与各个方向轴正割的曲面正面；s-x
、s-y
和s-z
分别为δv与各个方向轴正割的曲面反面，σ表示介质电导率；σa为异常体电导率；jm为磁流密度矢量；h
p
，hs为磁场强度；接下来对各积分得到积分结果后，进行面积加权平均计算，得到：
[0082][0083]
为节点(i,j,k)周围电导率的面积加权平均，即为
[0084][0085]
其中i,j,k分别表示x,y,z方向上的坐标表示；σ
i,j,k
为节点(i,j,k)周围电导率的面积加权值；σ
i+1,j,k
为节点(i+1,j,k)周围电导率的面积加权值，同理，变换下标就求得不同节点的面积加权值；δxi为x轴上的网格距，δyj为y轴上的网格距，δzk为z轴上的网格距；
[0086]
最后运用时间线性插值算法，得到各个节点电磁波场强数学表达式。
[0087]
step3、监测此时矿井磁场中的各点电导率的稳定性，根据极值来判定监测的时间长短；
[0088]
作为本发明的进一步方案，所述step3中包括具体如下：运用局部细分网格的思路找到时间步长与网格之间的关系来得到趋于稳定监测的时间区域，用时间插值公式表示为：
[0089][0090]
其中，h为与全区域中网格不重合的局部网格内部强度，h为与局部网格重合的网格中以及全区域网格中的强度；v为时间步
长；m为分的网格数量。
[0091]
用空间插值表示为：
[0092][0093][0094][0095]
局部细分网格有限差分算法的具体操作步骤如下：
[0096]
3)在包含细网格区的粗网格上利用步骤step 2中的式子(4)求得
[0097]
4)在局部细网格中给运用式子(5)迭代两次可以得到
[0098]
step4、运用中心差分近似法对step2中的结果计算差分网度以及最大网度值；
[0099]
作为本发明的进一步方案，所述step4中包括：运用欧拉法算出经过欧拉变换之后的差分方程为：
[0100][0101]
其中，令
[0102]
将时间导数运用中心差分近似有：经过差分近似后，得到新的差分方程为：
[0103][0104]
其中，r＝δt/μσδ2，此处的r为差分网度；为节点(i,j,k)周围电导率的面积加权平均；σ表示介质电导率，μ为介质磁导率；n表示当前时刻，n+1表示下一时刻；
[0105]
又因为欧拉法的稳定性条件是-1≤fe≤1，又因为在上述差分方程式子中fe≤1始终满足，故满上前半部分得到此处的即为最大网度。
[0106]
step5、根据巷道内满足的拉普拉斯方程，对整个矿井巷道空间求体积分，运用时间线性插值进行最后的整理，即可以得到矿井巷道内可控制边界的相应所需方程。
[0107]
作为本发明的进一步方案，所述step5具体包括如下：
[0108]
step5.1、根据laplace方程求积分，在矿井巷道中满足的laplace方程为
[0109]
运用有限差分法分割，对空腔巷道内部所示的δv内求解积分，得到
[0110][0111]
step5.2、时间线性插值；
[0112]
运用时间线性插值，运用欧拉法算出当前时刻经过欧拉变换之后的差分方程得到在矿井巷道内运用有限差分方程时相应所需方程为：
[0113][0114]
如图2所示，图中的含有空腔巷道的模型中，空腔均匀围岩介质的电阻率为ρ＝1ω
·
m，空腔内为不导电的空气。在整个运用有限差分方程来计算各点场强的过程中，均使用开发过程中掌子面空腔的场景来完成。此外，图3是对图2的网格化表示，可以更加明确的通过数学抽象表达的方式表示出不同节点的场强情况。
[0115]
上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

技术特征：
1.基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件的方程改进方法，其特征在于：所述方法的具体步骤如下：step1、根据探地雷达下正演的相关知识以及真实场景，运用非均匀网格技术，求解案例当前情况下的有限差分方程；step2、引入加权平均值的概念，对网格中的各点周围的电导率进行面积加权平均，得到离散化后有限差分算法各个节点的电磁波场强的数学表达式；step3、监测此时矿井磁场中的各点电导率的稳定性，根据极值来判定监测的时间长短；step4、运用中心差分近似法对step2中的结果计算差分网度以及最大网度值；step5、根据巷道内满足的拉普拉斯方程，对整个矿井巷道空间求体积分，运用时间线性插值进行最后的整理，即可以得到矿井巷道内可控制边界的相应所需方程。2.根据权利要求1所述的基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件的方程改进方法，其特征在于：所述step1中包括：需要运用数值计算的方法对整个空间进行非均匀网格技术分割，将无限大场的分布空间截断为有限的计算空间，在分割的过程中采用的是引入合适的边界条件，使得电磁场扩散至边界时能没有明显的反射，运用数学公式得到每一个分割点的有限差分方程如下：其中σ表示介质电导率；σ
a
为异常体电导率；j
m
为磁流密度矢量；h
p
，h
s
为磁场强度；μ为介质磁导率。3.根据权利要求1所述的基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件的方程改进方法，其特征在于：所述step2中包括的各个节点的电磁波场强的数学表达式如下：其中，n表示当前时刻，n+1表示下一时刻，i+1就是节点(i,j,k)x轴正方向上多一个单位的数学表示，i-1就是节点(i,j,k)x轴负方向上多一个单位的数学表示，j+1就是节点(i,j,k)y轴正方向上多一个单位的数学表示，j-1就是节点(i,j,k)y轴负方向上多一个单位的数学表示，k+1就是节点(i,j,k)z轴正方向上多一个单位的数学表示，k-1就是节点(i,j,k)z轴负方向上多一个单位的数学表示，为可探测区域总场t＝(n+1)δt时的磁场内中心节点的场强扩散方程；是x，y，z方向的平均网度；为x轴方向的网度；为y轴方向的网度；为z轴方向的网度；为节点(i,j,k)周围电导率的面
积加权平均值；i(t)为关断电流函数；s为所在节点处双侧封闭曲面面积，μ为介质磁导率；具体的操作步骤为：给有限差分方程的两边在δv内积分，得到接下来，根据高斯公式对上述式子进行计算并整理，有其中δv表示前后两个节点移动形成的体积差；s
x
、s
y
和s
z
分别为δv与各个方向轴正割的曲面正面；s-x
、s-y
和s-z
分别为δv与各个方向轴正割的曲面反面，σ表示介质电导率；σ
a
为异常体电导率；j
m
为磁流密度矢量；h
p
，h
s
为磁场强度；接下来对各积分得到积分结果后，进行面积加权平均计算，得到：行面积加权平均计算，得到：为节点9i,j,k)周围电导率的面积加权平均，即为其中i,j,k分别表示x,y,z方向上的坐标表示；σ
i,j,k
为节点(i,j,k)周围电导率的面积加权值；σ
i+1,j,k
为节点(i+1,j,k)周围电导率的面积加权值，同理，变换下标就求得不同节点的面积加权值；δx
i
为x轴上的网格距，δy
j
为y轴上的网格距，δz
k
为z轴上的网格距；最后运用时间线性插值算法，得到各个节点电磁波场强数学表达式。4.根据权利要求1所述的基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件的方程改进方法，其特征在于：所述step3中包括具体如下：运用局部细分网格的思路找到时间步长与网格之间的关系来得到趋于稳定监测的时间区域，用时间插值公式表示为：其中，h为与局部网格重合的网格中以及全区域网格中的强度；v为时间步长；m为分的网格数量。5.根据权利要求1所述的基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界
条件的方程改进方法，其特征在于：所述step4中包括：运用欧拉法算出经过欧拉变换之后的差分方程为：其中，令将时间导数运用中心差分近似有：经过差分近似后，得到新的差分方程为：其中，r＝δt/μσδ2，此处的r为差分网度；为节点(i,j,k)周围电导率的面积加权平均；σ表示介质电导率，μ为介质磁导率；n表示当前时刻，n+1表示下一时刻；又因为欧拉法的稳定性条件是-1≤f
e
≤1，又因为在上述差分方程式子中f
e
≤1始终满足，故满上前半部分得到此处的即为最大网度。6.根据权利要求3所述的基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件的方程改进方法，其特征在于：所述step5具体包括如下：step5.1、根据laplace方程求积分，在矿井巷道中满足的laplace方程为运用有限差分法分割，对空腔巷道内部所示的δv内求解积分，得到step5.2、时间线性插值；运用时间线性插值，运用欧拉法算出当前时刻经过欧拉变换之后的差分方程得到在矿井巷道内运用有限差分方程时相应所需方程为：

技术总结
本发明涉及基于有限差分法的电磁波波场中瓦斯隧道掌子面空腔及边界条件的方程改进方法，属隧道工程领域。本发明先根据探地雷达下正演的相关知识以及真实场景，运用非均匀网格技术，求解案例当前情况下的有限差分方程；其次引入加权平均值的概念，对网格中的各点周围的电导率进行面积加权平均，得到离散化后有限差分算法各个节点电导率值的抽象表达式；再次，监测此时矿井磁场中的各点电导率的稳定性，根据极值来判定监测的时间长短；然后，运用中心差分近似法对上方结果计算差分网度以及最大网度值；最后，根据巷道内满足的拉普拉斯方程对整个矿井巷道空间求体积分，运用时间线性插值进行最后的整理，即得到矿井巷道内可控制边界的相应所需方程。制边界的相应所需方程。制边界的相应所需方程。


技术研发人员：郑万波 董锦晓 谭春琳 冉丝蔓 朱榕 李磊 王耀 金正浩 钟佩玲
受保护的技术使用者：云南卫士盾科技有限公司
技术研发日：2023.05.07
技术公布日：2023/8/16