一种基于Radon变换的矿山微震层析成像定位算法

一种基于radon变换的矿山微震层析成像定位算法
技术领域
1.本发明涉及一种基于radon变换的矿山微震层析成像定位算法，属于矿山微震监测技术领域。


背景技术：

2.微震监测技术是矿山安全监测的利器，也是采矿工艺优化的重要手段。随着应用领域的扩展和应用要求的提升，微震定位精度亟待提升。


技术实现要素：

3.本发明的目的在于提供一种基于radon变换的矿山微震层析成像定位算法，以radon变换为基础，从傅里叶变换理论出发，推导建立了地震层析成像技术中的投影和反投影计算公式，求取了矿山微地震记录的慢度谱，提出了层析扫描参数的选择准则，建立了利用慢度谱中的能量极值识别微地震事件及其波型的准则和矿山微震自动定位的判据，实现了对微震事件的精准化定位计算。
4.本发明的目的是通过以下技术方案实现：
5.一种基于radon变换的矿山微震层析成像定位算法，包括如下步骤：
6.步骤s1，结合高精度微震监测系统的监测结果，对采动诱发的矿山微地震信号进行采集和精细化处理，将四维地震记录转换为球坐标系时空域二位记录；
7.步骤s2，根据地震层析成像原理和fourier变换理论，推导建立了地震层析成像技术中的投影和反投影计算公式；
8.步骤s3，将微震记录投影到球坐标时空域中，获得震源表达式，求取了矿山微地震记录的慢度谱；
9.步骤s4，提出了层析扫描参数的选择准则；
10.步骤s5，建立了利用慢度谱中的能量极值识别微地震事件及其波型的准则，提出了矿山微震自动定位的判据，实现了矿山微地震的自动定位；
11.步骤s6，提出了矿山微震震源参数计算模型，计算得出了微震震源有效物理参数。
12.进一步来说，步骤s1中将大地坐标系时空域(r,tr)表示的四维地震记录g(xr,yr,zr,tr)转换成球坐标系时空域(r,t)的二维地震波记录f(r,t)，公式如下：
[0013][0014]
其中，xs、ys、zs：表示大地坐标系中的震源位置；ts：表示微地震发震时刻；xr、yr、zr：表示大地坐标系中的检波器位置；tr：表示地震波到达rr(xr,yr,zr)处的时刻；g(xr,yr,zr,tr)：表示大地坐标时空域中的地震波记录；r，t：表示球坐标时空域中的地震波传播距离、地震波旅行时，在球坐标系中，震源位置r＝0、激发时刻t＝0，f(r,t)表示球坐标时空域
中的地震波记录。
[0015]
进一步来说，步骤s2中地震层析成像技术中的投影和反投影计算公式：
[0016][0017][0018][0019]
式中，检波器接收到的地震信号为一个二维函数f(r,t)，r表示震源位置到检波器位置之间的距离，t表示接收时间长度。
[0020]
进一步来说，步骤s3中微地震记录在球坐标时空域中的投影表达式如下：
[0021][0022]
其中：i＝1,2,3...n：表示接收道序号；j＝1,2,3...m：表示可能的震源位置序号；k＝1,2,3...k：表示扫描参数值的序号；t表示截距；pk表示斜率；g(x
ri
,y
ri
,z
ri
,t
ri
)表示微地震记录。
[0023]
进一步来说，步骤s4中扫描参数pk的采样间隔
△
p应满足下列表达式：
[0024][0025]
其中：f
max
是所有微地震接收道中信号的最大频率；r
max
是震源位置与各检波器之间距离的最大值。
[0026]
进一步来说，步骤s5中震源自动定位计算步骤：
[0027]
第一步，分析并选定矿区内事故易发区域作为风险区域，风险区域的范围可以设定为rz(x
min
,x
max
,y
min
,y
max
,z
min
,z
max
)；
[0028]
第二步，对参与定位计算的微震记录进行预处理，保证参与层析成像投影的微地震记录，是经过自动识别的单个微地震事件的数据记录；
[0029]
第三步，在风险区内，利用表达式对单个微地震事件的数据记录进行层析成像投影，获得该微地震事件的慢度谱；
[0030]
第四步，将慢度谱中最大的能量极值所对应的空间坐标值作为自动定位的判据。
[0031]
进一步来说，最大能量值的计算公式为bs(ps,ts)＝max[bj(pk,tj)]，其中：rs(xs,ys,zs)：最大能量值对应的空间坐标值；ts：最大能量值对应的时间坐标值；ps：最大能量值对应的慢度坐标值；bs(ps,t)：最大能量值对应的慢度时间信号。
[0032]
进一步来说，步骤s6中微震震源有效物理参数计算步骤如下：
[0033]
第一步，首先将一个微地震事件在各检波器上的观测信号进行层析成像投影，获得能量最大值bs(ps,0)；
[0034]
第二步，利用震源层析成像自动定位的最大能量准则定位判据，确定最大能量值对应的慢度时间信号bs(ps,t),表达式为：bs(ps,t)＝∫f(r,t+psr)dr；
[0035]
第三步，利用震源位置与检波器位置的间距计算几何扩散补偿因子α(r)；其中，震源函数的计算公式：f1(0,t)＝α(r)bs(ps,t)；为几何扩散能量补偿因子；
[0036]
第四步，由ps道信号与几何扩散补偿因子α(r)相乘，可以直接求取该地震事件的震源函数f1(0,t)。
[0037]
进一步来说，步骤s6中微震震源破裂半径计算步骤如下：
[0038]
第一步，利用ω-2
模型计算中小地震的震源位移谱，可以表示为：其中ω(f)表示ω-2
模型的震源谱；ω(0)表示震源谱的零频极限值；fc表示震源谱的拐角频率；
[0039]
第二步，利用零频极限值ω(0)和fc拐角频率，可以计算该微地震事件的地震矩、震源破裂半径和应力降。其中震源破裂半径的计算公式为其中：v＝1/ps是微地震波传播速度；
[0040]
第三步，提出将ω-2
模型应用于矿山微地震事件的研究，建立微震震源振动位移函数的振幅谱u(f)计算模型f(t)表示震源振动速度函数；f(f)表示经傅里叶变换后可得微地震震源振动速度的振幅谱
[0041]
第四步，利用最小二乘法求解震源谱ω(f)和微震震源谱u(f)具有最小残差率，进而求取微震事件的破裂半径
[0042]
采用上述技术方案，具有以下有益效果：本方法相较传统定位算法，弱化了对初至到时的强烈依赖，初步实现了矿山微震事件的全自动定位计算，为矿山煤岩动力灾害的精准预警、防治提供了技术支撑。此外，所求取的微震震源函数及其参数(震源有效辐射能量、微地震零频极限值、拐角频率、震源破裂半径等)，为研究矿山岩体应变过程甚至破裂失稳问题提供了重要依据。
附图说明
[0043]
为了更清楚地说明本发明的实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是示例性的，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图引伸获得其它的实施附图。
[0044]
本说明书所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，故不具技术上的实质意义，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功
效及所能达成的目的下，均应仍落在本发明所揭示的技术内容得能涵盖的范围内。
[0045]
图1为本发明提供的流程示意图。
[0046]
图2为本发明提供的层析扫描参数间隔选择示意图。
[0047]
图3为本发明提供的理论模型和观测系统示意图。
[0048]
图4为本发明提供的理论模型中的微地震震源函数示意图。
[0049]
图5为本发明提供的理论模型和观测系统(震源在风险区内)示意图。
[0050]
图6为本发明提供的微地震理论记录(震源在风险区内)示意图。
[0051]
图7为本发明提供的层析成像自动定位结果图(震源在风险区内)示意图。
[0052]
图8为本发明提供的层析成像流程示意图。
[0053]
图9为本发明提供的采区微地震事件平面分布图。
[0054]
图10为本发明提供的工作面震源事件水平投影图。
[0055]
图11为本发明提供的工作面震源走向投影图。
[0056]
图12为本发明提供的工作面震源倾向投影图。
具体实施方式
[0057]
下面结合附图对本发明的具体实施方式作进一步说明。在此需要说明的是，对于这些实施方式的说明用于帮助理解本发明，但并不构成对本发明的限定。此外，下面所描述的本发明各个实施方式中所涉及的技术特征只要彼此之间未构成冲突就可以相互组合。
[0058]
实施例
[0059]
参见图1所示，一种基于radon变换的矿山微震层析成像定位算法，包括如下步骤：
[0060]
步骤s1，结合高精度微震监测系统的监测结果，对采动诱发的矿山微地震信号进行采集和精细化处理，将四维地震记录转换为球坐标系时空域二位记录；
[0061]
步骤s2，根据地震层析成像原理和fourier变换理论，推导建立了地震层析成像技术中的投影和反投影计算公式；
[0062]
步骤s3，将微震记录投影到球坐标时空域中，获得震源表达式，求取了矿山微地震记录的慢度谱；所述慢度谱为在球坐标时空域中微地震记录的投影结果；
[0063]
步骤s4，提出了层析扫描参数的选择准则；
[0064]
步骤s5，建立了利用慢度谱中的能量极值识别微地震事件及其波型的准则，提出了矿山微震自动定位的判据，实现了矿山微地震的自动定位；
[0065]
步骤s6，提出了矿山微震震源参数计算模型，计算得出了微震震源有效辐射能量、微地震零频极限值、拐角频率、震源破裂半径等震源物理参数。
[0066]
具体来说，在步骤s1中，将四维地震记录转换为球坐标系时空域二位记录的方式为：通过下列坐标变换公式，可以将大地坐标系时空域(r,tr)表示的四维地震记录g(xr,yr,zr,tr)转换成球坐标系时空域(r,t)的二维地震波记录f(r,t)：
[0067][0068]
其中，xs、ys、zs：表示大地坐标系中的震源位置；ts：表示微地震发震时刻；xr、yr、
zr：表示大地坐标系中的检波器位置；tr：表示地震波到达rr(xr,yr,zr)处的时刻；g(xr,yr,zr,tr)：表示大地坐标时空域中的地震波记录；r，t：表示球坐标时空域中的地震波传播距离、地震波旅行时，在球坐标系中，震源位置r＝0、激发时刻t＝0；f(r,t)：表示球坐标时空域中的地震波记录。
[0069]
步骤s2中地震层析成像技术中的投影和反投影计算公式如下：
[0070][0071][0072][0073]
式中，假设检波器接收到的地震信号为一个二维函数f(r,t)，这里r表示震源位置(坐标原点)到检波器位置之间的距离，t表示接收时间长度。当检波器接收到的地震信号为一个二维函数f(r,t)时，沿r-t平面上任一条直线l:t＝t+pr上的积分，称为f(r,t)的radon正变换。参数p是直线l的斜率，其物理意义是地震波传播速度的倒数，即慢度。
[0074]
步骤s3利用地震层析成像技术，微地震记录在球坐标时空域中的投影表达式如下：
[0075][0076]
在实际投影过程中，由于震源位置是未知的，无法求取震源到检波器之间的距离，因此必须首先假定一个震源位置为r
sj
(x
sj
,y
sj
,z
sj
)，发震时间和传播速度也是未知的，选择具有不同截距t和不同斜率pk(常常称为扫描参数)的直线，并沿着直线对各接收道的微地震记录g(x
ri
,y
ri
,z
ri
,t
ri
)进行叠加，最后获得一个相对于该震源位置的投影，其离散表达式为：
[0077][0078]
其中，i＝1,2,3...n：表示接收道序号；j＝1,2,3...m：表示可能的震源位置序号；k＝1,2,3...k：表示扫描参数值的序号。
[0079]
步骤s4中层析扫描参数的选择准则为：根据数字信号处理与分析原理，对连续型函数实行离散化处理，必须满足离散化采样定理，否则就会产生假频现象。扫描参数pk(其中k＝1,2,3...k)的间隔
△
p＝p
k-p
k-1
选定，也必须满足离散化采样定理。利用扫描参数对所有接收道上的微地震信号进行扫描(以为r＝0基准)时，对于任何一道上的微地震波，应当在一个周期t内至少达到两次扫描(如图2所示)。
[0080]
扫描参数pk的采样间隔
△
p应满足下列表达式：
[0081][0082]
其中：f
max
是所有微地震接收道中信号的最大频率；r
max
是震源位置与各检波器之间距离的最大值。
[0083]
步骤s5中震源层析成像的能量极值：
[0084]
微地震记录在球坐标时空域中的投影，是沿着具有不同截距t和不同斜率pk的直线，对各接收道的微地震记录g(x
ri
,y
ri
,z
ri
,t
ri
)进行叠加完成的。这条直线的表达式为：
[0085]
tr＝t+pr
ꢀꢀꢀꢀꢀꢀ
(8)
[0086]
由前文讨论可知：为震源与检波器之间的间距，现令：扫描参数代入上述直线方程(4-20)，当t＝0时：
[0087][0088]
根据地震波传播理论可知：表达式(9)恰好是均匀介质中透射波在球坐标系中的时距方程，震源位置和发震时间位于球坐标时空域(r,t)中的原点，其中v是传播速度，由于p是传播速度v的倒数，故称扫描参数p为慢度。
[0089]
由于在球坐标系中，震源位置r＝0、激发时刻t＝0，所以当沿着穿过原点的直线(时距方程)对地震信号的振幅的绝对值进行叠加时，可以获得一个极大值。由于地震信号的振幅与能量成平方关系，因此将极大值简称为能量极值。分析能量极值，可以提取能量极值对应的坐标参数：
[0090]
·rs
：能量极值对应的空间坐标值向量(xs,ys,zs)；
[0091]
·
ts：能量极值对应的时间坐标值；
[0092]
·
ps：能量极值对应的慢度坐标值；
[0093]
·
b(ps,t)：能量极值对应的慢度时间信号。
[0094]
如果是对连续监测的微地震记录进行层析成像投影，那么，在获得的慢度谱上，可能存在多个能量极值的情况，分析这些能量极值对应的坐标参数的异同情况，可以进一步解释微地震事件及其波型特征，如表1所示。
[0095]
表1能量极值与微地震事件及其波型的对应关系
[0096][0097]
利用表1所示的能量极值与微地震事件及其波型的对应关系，通过解释慢度谱上的能量极值，可以容易和快速地读取其对应的坐标参数值。对能量极值的坐标参数值进行不同组合和分析，达到分析微地震震源参数和传播介质结构参数的目的。
[0098]
步骤s5中震源自动定位计算的步骤为：
[0099]
第一步，分析并选定矿区内事故易发区域作为风险区域，风险区域的范围可以设定为rz(x
min
,x
max
,y
min
,y
max
,z
min
,z
max
)；
[0100]
第二步，对参与定位计算的微震记录进行预处理，保证参与层析成像投影的微地震记录，是经过自动识别的单个微地震事件的数据记录。
[0101]
第三步，在风险区内，利用表达式(6)，对单个微地震事件的数据记录进行层析成像投影，获得该微地震事件的慢度谱，根据表1所示的能量极值与微地震事件及其波型的对应关系。
[0102]
最大能量值的计算公式如下：
[0103]bs
(ps,ts)＝max[bj(pk,tj)]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0104]
如果风险区的范围设置足以包含所关注的微地震震源位置，那么，利用上述自动定位判据，可以自动提取最大能量值对应的震源参数如下：
[0105]
·
最大能量值对应的空间坐标值：rs(xs,ys,zs)；
[0106]
·
最大能量值对应的时间坐标值：ts；
[0107]
·
最大能量值对应的慢度坐标值：ps；
[0108]
·
最大能量值对应的慢度时间信号：bs(ps,t)。
[0109]
第四步，将慢度谱中最大的能量极值所对应的空间坐标值作为自动定位的判据，如表2所示。
[0110]
表2最大能量准则的定位判据
[0111][0112]
步骤s6中求解特征频率等震源参数的方法：
[0113]
第一步，首先将一个微地震事件在各检波器上的观测信号进行层析成像投影，获得能量最大值bs(ps,0)；
[0114]
第二步，利用震源层析成像自动定位的最大能量准则定位判据，确定最大能量值对应的慢度时间信号(称为ps道信号)bs(ps,t)，实际上bs(ps,t)可以用下列表达式表示：
[0115]bs
(ps,t)＝∫f(r,t+psr)dr
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0116]
第三步，利用震源位置与检波器位置的间距计算几何扩散补偿因子α(r)。其中，震源函数的计算公式：
[0117]
f1(0,t)＝α(r)bs(ps,t)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0118]
其中：为几何扩散能量补偿因子。
[0119]
第四步，由ps道信号与几何扩散补偿因子α(r)相乘，可以直接求取该地震事件的震源函数f1(0,t)。
[0120]
步骤s6中求解震源破裂半径的方法：
[0121]
第一步，利用ω-2
模型计算中小地震的震源位移谱，可以表示为：
[0122][0123]
其中，ω(f)表示ω-2
模型的震源谱；ω(0)表示震源谱的零频极限值；fc表示震源谱的拐角频率。
[0124]
第二步，利用零频极限值ω(0)和fc拐角频率，可以计算该微地震事件的地震矩、震源破裂半径和应力降。其中震源破裂半径的计算公式为：
[0125][0126]
其中：v＝1/ps是微地震波传播速度。
[0127]
第三步，提出将ω-2
模型应用于矿山微地震事件的研究，建立微震震源振动位移函数的振幅谱u(f)计算模型。
[0128]
假设矿山微地震监测到的震源振动速度函数(震源函数)为f(t)，经傅里叶变换后可得微地震震源振动速度的振幅谱为f(f)，则微震震源振动位移函数的振幅谱为：
[0129][0130]
第四步，利用最小二乘法求解震源谱ω(f)和微震震源谱u(f)具有最小残差率，进
而求取微震事件的破裂半径。
[0131]
把ω(0)和fc作为独立变量，利用最小二乘法，使ω-2
模型震源谱ω(f)和微地震震源谱u(f)具有最小残差率：
[0132][0133]
这样可以确定该微地震事件的零频极限值ω(0)和拐角频率fc，进而可以计算出该微震事件震源破裂半径。
[0134]
案例估算
[0135]
示例：
[0136]
步骤1.1，对采动影响区域进行初步估算，确立相应的采区范围，建立高精度微震监测系统体系。
[0137]
通过理论模型数据和实测微地震信号的试算，描述微地震震源层析成像方法的应用，同时验证该方法的正确性和实用性。
[0138]
第一步，确立采区影响范围。东西向从3000m到3860m、南北向从5000m到5460m、地下深度从-600m到-900m；介质速度为2500m/s，如3所示。
[0139]
第二步，优化设计微震观测系统。在该采区内的采矿巷道中埋置了10个检波器(其中道号i＝1,2,3,
…
10)，各检波器在采区内的空间位置，如图2中的蓝色圆圈所示，其中r
ri
(x
ri
,y
ri
,z
ri
)为检波器位置坐标。
[0140]
第三步，确立仪器监测参数。为了讨论方便，以微地震监测仪器启动时刻作为计时系统，t＝0时开始记录微地震信号，假设仪器的采样间隔
△
t＝0.0005s，采样点数为1024。
[0141]
步骤1.2，为简化计算，在采区影响范围内划分出一个极易发生事故的风险区，并将风险区进一步划分成网格体。同时计算检波器到网格节点之间的距离，最后计算出层析扫描参数。
[0142]
第一步，确立风险区范围。假设事故易发的风险区域为：东西向从3000m到3500m、南北向从5000m到5460m、地下深度从-630m到-770m，如图3中所示的风险区域。
[0143]
第二步，进行网格体划分。将风险区域按步长为10m划分成三维网格体，网格体的节点假设为可能的震源位置，其东西向e、南北向n和深度方向h的坐标值如下：
[0144][0145]
第三步，计算检波器与网格节点间的距离。根据各检波器位置坐标r
ri
(x
ri
,y
ri
,z
ri
)与风险区域内三维网格体节点位置坐标(x
l
,ym,zn)，可以计算出检波器到网格节点之间的距离r
ilmn
，其计算公式如下：
[0146][0147]
由此可知：检波器与网格节点之间的最大距离r
max
＝max(r
lmn
)。
[0148]
第四步，确立层析成像扫描参数。地震波内的传播速度是难以获取的，只能根据地质知识和经验给出一个速度范围。假设地震波传播速度范围在v
min
＝1000m/s到v
max
＝
3000m/s之间，于是扫描参数在p
min
＝1/v
max
＝0.0003s/m到p
max
＝1/v
min
＝0.001s/m之间。另外，微地震信号的频率在几赫兹到1000多赫兹之间，根据扫描参数的选择准则，可以选定一个扫描参数间隔
△
p
lmn
，于是扫描参数为：
[0149]
p
lmn_k
＝p
min
+(k-1)
△
p
lmn
(k＝1、2、3...k
kmn
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(19)
[0150]
其中：k
lmn
＝(p
max-p
min
)/
△
p
lmn
为采样点数。
[0151]
步骤1.3，层析成像自动定位计算。基于上一步中理论分析，利用层析成像自动定位算法获取矿震频繁发生的位置。
[0152]
第一步，在进行层析成像自动定位运算之间，设计两种情况的震源，一种情况是震源位置位于风险区内，另一种情况是震源位置位于风险区外。两种情况除了震源位置不同外，震源函数是相同的，均可以用下列公式表达：
[0153]
f(t)＝ae-βt
sin(2πft+φ)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(20)
[0154]
其中：振幅a＝100、衰减系数β＝0.8、频率f＝80hz和相位φ＝0。
[0155]
假设，微地震震源函数的震动延续时间为0.05s，采样间隔
△
t＝0.0005s、采样点数为100，如图4所示。
[0156]
在采区内，rs(xs,ys,zs)＝(3200,5200,-700)位置处，当时间tf＝0.05s时，有一个微地震事件发生，如图5中圆点所示。
[0157]
第二步，在各检波器接收到的微地震信号f(i,t)(如图6所示)，可以由下列公式表达：
[0158][0159]
其中：i＝1,2,3...10表示接收道号；f0(i,t)＝0表示没有接收到微地震信号时的道记录；为震源到各检波器之间的距离；表示微地震信号的延续时间。
[0160]
图6中所示的横轴是时间轴(仪器启动时刻为零点)，左纵轴是以站点编号排列的记录道号，右纵轴表示各站点位置相对于震源位置之间的距离值ri，由于各道ri值不同，所以各道接收到的地震波到时t
ri
也不同，ri值小的道先于ri值大的道接收到地震波。图6中所示的各道接收到的地震波振幅值做了道间平衡处理。
[0161]
第三步，计算微地震信号的慢度时间域，并求解最大能量值及其对应坐标参数。为了确定微地震震源位置rs(xs,ys,zs)，首先将各检波器接收到的微地震信号f(i,t)投影到慢度时间域中：
[0162][0163]
其中：j＝1,2,3...1024；k＝1,2,3...k
kmn
。
[0164]
第四步，根据最大能量准则的定位判据，即可获得该微震事件的定位结果。具体而言，求出能量最大值max(b
lmn
(pk,j))＝0.58539及其对应坐标参数，如图7所示。
[0165]
图7(a)表示能量最大值max(b
lmn
(pk,j))＝0.58539处对应的慢度谱，从图中可见，能量最大值对应的坐标为(0.053,2500)，其中2500表示速度v＝2500m/s，0.053表示微地震事件的震源时间ts，能量最大值对应的震源时间坐标ts＝0.053s与微地震发震时间tf＝
0.05s相差0.003s，产生这个差别的原因是：震源函数峰值处时刻(peak time)与发震时刻(trigger time)存在着0.003s的时差。由于层析成像的结果是对各道微地震记录的振幅值叠加，其极值是各道记录中大振幅值叠加的结果，因此，层析成像能量最大值对应的时间ts是微地震事件震动强度最大的时刻，而不是微地震事件启动时刻，也就是前文所属的成核时间。
[0166]
图7(b)、(c)和(d)分别代表能量最大值处对应的东西向、南北向和深度方向的慢度谱，从图中可见，能量最大值对应的坐标分别为(3200,2500)、(5200,2500)和(-700,2500)，由此可见，微地震事件的震源位置为rs(xs,ys,zs)＝(3200,5200,-700)，地震波传播速度v＝2500m/s。
[0167]
实施例
[0168]
具体操作流程如图8所示：
[0169]
步骤1，初始化模块。
[0170]
每当启动一个新的监测项目时，必须启动并运行“初始化模块”：初始化模块包括风险区坐标和速度范围定义、观测系统数据输入。如表3和表4所示。
[0171]
根据表3和表4所示的输入数据，初始化模块将自动计算出下列数据：
[0172]
第一步，风险区的3d网格节点数据矩阵grids_3c＝[x,y,z]；
[0173]
第二步，风险区网格节点与检波器间距的数据矩阵radius＝[行号＝网格节点号，列号＝检波器编号]；
[0174]
第三步，基于监测工区网格节点到检波器之间，微地震波旅行时差的数据矩阵
△
t＝[g
×v×
r]；
[0175]
第四步，显示风险区的图和基于监测工区的观测系统图。
[0176]
表3风险区坐标参数和速度模型参数
[0177][0178]
表4观测系统参数
[0179][0180]
步骤2，震源层析成像模块。
[0181]
在层析成像参数初始化的基础上，每当采集到一个微地震事件时，系统自动调用
并运行震源层析成像模块，并依次完成如下功能：
[0182]
第一步，自动读取微地震事件信号数据：包括该微地震事件的微地震记录数据s＝[t
×
r]和微地震记录数据相关的参数，道数；采样间隔；记录长度等；
[0183]
第二步，调用监测工区网格节点到检波器之间，微地震波旅行时差的数据矩阵
△
t＝[g
×v×
r]；
[0184]
第三步，计算微地震记录数据层析成像投影成像矩阵tomo＝[g
×v×
t]，自动判定并求出该微地震事件的震源位置、传播速度、发震时间和震源函数；
[0185]
第四步，自动求取震源函数的能量、零频极限值、拐角频率和震源破裂半径；
[0186]
第五步，自动保存震源位置、传播速度、发震时间、能量、零频极限值、拐角频率和震源破裂半径进入数据库，将震源函数保存为数据文件。根据需要还可以将上述震源参数显示出来。
[0187]
步骤3，震源投影展示模块。
[0188]
步骤3.1，加载矿区底图。
[0189]
加载矿区底图用于加载所要分析的采区底图，为该采区微地震事件的显示、分析和解释作准备。
[0190]
步骤3.2，震源水平投影。
[0191]
如图9所示，震源水平投影是将微地震事件显示在采区开拓平面图中，从图中可以分析该时间段内和能量范围内的微地震事件平面分布特征。
[0192]
步骤3.3，计算微震事件统计分析结果，并投影到平面图上，其中x轴为走向坐标，y轴为倾向坐标。
[0193]
震源在空间域的投影是将微地震事件在工作面平面图、走向剖面图、倾向剖面图中进行投影，这三个投影图构成微地震事件的空间分布。主要包括下列三个模块
[0194]
第一步，计算并制作水平投影图。如10所示，水平投影图是将微地震事件投影到工作面平面图上，通过微地震水平投影图上的微地震事件位置，可以在平面上分析解释微地震活动特征。
[0195]
第二步，计算并制作走向投影图。如图11所示，走向投影图是将微地震事件投影到工作面走向剖面图上，通过走向投影图找到微地震事件的走向位置及层位，可以在平面上分析解释微地震活动特征，从而判断上覆岩层的应变高度。
[0196]
第三步，计算并制作倾向投影图。如图12所示，倾向投影图是将微地震事件投影到工作面倾向剖面图上，通过倾向投影图找到微地震事件的倾向位置，可以在平面上分析解释微地震活动特征，便于找到微地震事件距离煤层顶底板的距离及分布规律。
[0197]
以上结合附图对本发明的实施方式作了详细说明，但本发明不限于所描述的实施方式。对于本领域的技术人员而言，在不脱离本发明原理和精神的情况下，对这些实施方式进行多种变化、修改、替换和变型，仍落入本发明的保护范围内。

技术特征：
1.一种基于radon变换的矿山微震层析成像定位算法，其特征在于：包括如下步骤：步骤s1，结合高精度微震监测系统的监测结果，对采动诱发的矿山微地震信号进行采集和精细化处理，将四维地震记录转换为球坐标系时空域二位记录；步骤s2，根据地震层析成像原理和fourier变换理论，推导建立了地震层析成像技术中的投影和反投影计算公式；步骤s3，将微震记录投影到球坐标时空域中，获得震源表达式，求取了矿山微地震记录的慢度谱；步骤s4，提出了层析扫描参数的选择准则；步骤s5，建立了利用慢度谱中的能量极值识别微地震事件及其波型的准则，提出了矿山微震自动定位的判据，实现了矿山微地震的自动定位；步骤s6，提出了矿山微震震源参数计算模型，计算得出了微震震源有效物理参数。2.根据权利要求1所述的基于radon变换的矿山微震层析成像定位算法，其特征在于：步骤s1中将大地坐标系时空域(r,t
r
)表示的四维地震记录g(x
r
,y
r
,z
r
,t
r
)转换成球坐标系时空域(r,t)的二维地震波记录f(r,t)，公式如下：其中，x
s
、y
s
、z
s
：表示大地坐标系中的震源位置；t
s
：表示微地震发震时刻；x
r
、y
r
、z
r
：表示大地坐标系中的检波器位置；t
r
：表示地震波到达r
r
(x
r
,y
r
,z
r
)处的时刻；g(x
r
,y
r
,z
r
,t
r
)：表示大地坐标时空域中的地震波记录；r，t：表示球坐标时空域中的地震波传播距离、地震波旅行时，在球坐标系中，震源位置r＝0、激发时刻t＝0，f(r,t)表示球坐标时空域中的地震波记录。3.根据权利要求1所述的基于radon变换的矿山微震层析成像定位算法，其特征在于：步骤s2中地震层析成像技术中的投影和反投影计算公式：步骤s2中地震层析成像技术中的投影和反投影计算公式：步骤s2中地震层析成像技术中的投影和反投影计算公式：式中，检波器接收到的地震信号为一个二维函数f(r,t)，r表示震源位置到检波器位置之间的距离，t表示接收时间长度。4.根据权利要求1所述的基于radon变换的矿山微震层析成像定位算法，其特征在于：步骤s3中微地震记录在球坐标时空域中的投影表达式如下：
其中：i＝1,2,3...n：表示接收道序号；j＝1,2,3...m：表示可能的震源位置序号；k＝1,2,3...k：表示扫描参数值的序号；t表示截距；p
k
表示斜率；g(x
ri
,y
ri
,z
ri
,t
ri
)表示微地震记录。5.根据权利要求1所述的基于radon变换的矿山微震层析成像定位算法，其特征在于：步骤s4中扫描参数p
k
的采样间隔
△
p应满足下列表达式：其中：f
max
是所有微地震接收道中信号的最大频率；r
max
是震源位置与各检波器之间距离的最大值。6.根据权利要求1所述的基于radon变换的矿山微震层析成像定位算法，其特征在于：步骤s5中震源自动定位计算步骤：第一步，分析并选定矿区内事故易发区域作为风险区域，风险区域的范围可以设定为rz(x
min
,x
max
,y
min
,y
max
,z
min
,z
max
)；第二步，对参与定位计算的微震记录进行预处理，保证参与层析成像投影的微地震记录，是经过自动识别的单个微地震事件的数据记录；第三步，在风险区内，利用表达式对单个微地震事件的数据记录进行层析成像投影，获得该微地震事件的慢度谱；第四步，将慢度谱中最大的能量极值所对应的空间坐标值作为自动定位的判据。7.根据权利要求6所述的基于radon变换的矿山微震层析成像定位算法，其特征在于：最大能量值的计算公式为b
s
(p
s
,t
s
)＝max[b
j
(p
k
,t
j
)]，其中：r
s
(x
s
,y
s
,z
s
)：最大能量值对应的空间坐标值；t
s
：最大能量值对应的时间坐标值；p
s
：最大能量值对应的慢度坐标值；b
s
(p
s
,t)：最大能量值对应的慢度时间信号。8.根据权利要求1所述的基于radon变换的矿山微震层析成像定位算法，其特征在于：步骤s6中微震震源有效物理参数计算步骤如下：第一步，首先将一个微地震事件在各检波器上的观测信号进行层析成像投影，获得能量最大值b
s
(p
s
,0)；第二步，利用震源层析成像自动定位的最大能量准则定位判据，确定最大能量值对应的慢度时间信号b
s
(p
s
,t),表达式为：b
s
(p
s
,t)＝∫f(r,t+p
s
r)dr；第三步，利用震源位置与检波器位置的间距计算几何扩散补偿因子α(r)；其中，震源函数的计算公式：f1(0,t)＝α(r)b
s
(p
s
,t)；为几何扩散能量补偿因子；第四步，由p
s
道信号与几何扩散补偿因子α(r)相乘，可以直接求取该地震事件的震源函数f1(0,t)。9.根据权利要求8所述的基于radon变换的矿山微震层析成像定位算法，其特征在于：步骤s6中微震震源破裂半径计算步骤如下：
第一步，利用ω-2
模型计算中小地震的震源位移谱，可以表示为：其中ω(f)表示ω-2
模型的震源谱；ω(0)表示震源谱的零频极限值；f
c
表示震源谱的拐角频率；第二步，利用零频极限值ω(0)和f
c
拐角频率，可以计算该微地震事件的地震矩、震源破裂半径和应力降。其中震源破裂半径的计算公式为其中：v＝1/p
s
是微地震波传播速度；第三步，提出将ω-2
模型应用于矿山微地震事件的研究，建立微震震源振动位移函数的振幅谱u(f)计算模型f(t)表示震源振动速度函数；f(f)表示经傅里叶变换后可得微地震震源振动速度的振幅谱第四步，利用最小二乘法求解震源谱ω(f)和微震震源谱u(f)具有最小残差率，进而求取微震事件的破裂半径

技术总结
本发明公开了一种基于Radon变换的矿山微震层析成像定位算法，包括步骤：S1，将四维地震记录转换为球坐标系时空域二位记录；S2，建立地震层析成像技术中的投影和反投影计算公式；S3，将微震记录投影到球坐标时空域中，获得震源表达式，求取矿山微地震记录的慢度谱；S4，提出了层析扫描参数的选择准则；S5，利用慢度谱中的能量极值识别微地震事件及其波型的准则，提出了矿山微震自动定位的判据，实现了矿山微地震的自动定位；S6，提出了矿山微震震源参数计算模型，计算得出了微震震源有效物理参数。相较传统定位算法，弱化了对初至到时的强烈依赖，初步实现了矿山微震事件的全自动定位计算，为矿山煤岩动力灾害的精准预警、防治提供了技术支撑。了技术支撑。了技术支撑。


技术研发人员：朱权洁 缪华祥 刘晓云 隋龙琨 杨磊 朱斯陶 陈学习
受保护的技术使用者：华北科技学院（中国煤矿安全技术培训中心）
技术研发日：2023.03.23
技术公布日：2023/8/24