✨ 长期致力于矿震、微震、球面波、点源与纵波、点涡与横波、采空区、正反演、震源定位、多层介质研究工作擅长数据搜集与处理、建模仿真、程序编写、仿真设计。✅ 专业定制毕设、代码✅如需沟通交流点击《获取方式》1含采空区岩土体的球面波散射解析模型与敲击试验验证基于弹性波动理论将采空区等效为圆柱形孔洞建立二维平面应变下的纵波入射散射场表达式。采用波函数展开法将入射波、散射波均用Hankel函数表示并在孔洞边界施加应力自由条件求解未知系数。推导出孔洞后方特定区域的加速度时程曲线发现存在指向波源的负加速度反转区。设计带孔有机玻璃板敲击试验使用加速度计阵列测量。在孔洞直径30mm、板厚8mm条件下孔后2倍孔径处监测到负向峰值达入射波峰值35%的反转信号与理论预测误差小于8%。利用有限元仿真分析采空区群的影响当空区间距小于直径时发生波导效应能量沿空区间柱体集中夹岩柱水平加速度可放大4倍。2多层水平介质中球面波前正演与震源深度几何平均定位法基于惠更斯原理和斯涅耳定律推导点源在N层水平介质中传播的走时曲线显式表达式。建立双曲线型波前曲面方程给出各层界面折射点的迭代求解方法:利用等时线从台站逆推至最后一层界面然后逐层向上计算。提出几何平均法确定震源深度:由两个不同偏移距台站的纵波初至时差构造关于深度的二次方程取正根作为深度估计。在二层混凝土试验中震源深度20cm几何平均法计算深度20.8cm误差4%。进一步结合发震时刻反演提出正反演联用法先假设震源水平坐标为零扫描深度得到时刻再反演水平偏移。该方法对倾斜界面同样适用通过坐标旋转将倾斜层转换成等效水平层。3任意观测系统下的波前正演定位算法与鲁棒性分析将发震时刻作为未知量加入非线性方程组利用各台站之间的到时差构造新的观测方程。采用信赖域反射算法求解最小二乘问题目标函数为理论走时与观测到时差的平方和。每次迭代中根据当前震源位置和发震时刻计算各个台站的理论走时并通过数值方法求解球面波在多层介质中的射线路径。为提高定位稳定性引入M估计器权重函数降低异常到时的影响采用Huber函数阈值设为1.5倍残差中位数。在模拟采空区扰动的环境下对10个随机台站进行定位测试传统Geiger法平均误差7.2m波前正演法降至2.3m。将该方法应用到某煤矿微震监测数据定位事件与采掘面位置吻合度提升39%。import numpy as np from scipy.optimize import least_squares class SphericalWaveInLayers: def __init__(self, layer_thick, layer_vp): self.nlayers len(layer_thick) self.h np.array(layer_thick) self.vp np.array(layer_vp) def travel_time(self, source_depth, offset): # ray tracing using Snells law if source_depth 0: source_depth 0 cum_h np.cumsum(self.h) if source_depth cum_h[0]: src_layer 0; z_in_layer source_depth else: src_layer np.searchsorted(cum_h, source_depth) z_in_layer source_depth - (cum_h[src_layer-1] if src_layer0 else 0) # simplified: constant ray parameter p p_guess 1.0/self.vp[0] def f(p): tt 0.0 x 0.0 for i in range(src_layer, self.nlayers): if isrc_layer: dz self.h[i] - z_in_layer else: dz self.h[i] if self.vp[i]*p 1: return 1e9 angle np.arcsin(self.vp[i]*p) tt dz/(self.vp[i]*np.cos(angle)) x dz*np.tan(angle) return x - offset from scipy.optimize import brentq try: p_opt brentq(f, 1e-6, 1.0/self.vp[-1]) tt 0.0 for i in range(src_layer, self.nlayers): if isrc_layer: dz self.h[i]-z_in_layer else: dz self.h[i] angle np.arcsin(self.vp[i]*p_opt) tt dz/(self.vp[i]*np.cos(angle)) return tt except: return np.inf def wavefront_inversion(self, stations_xy, arrivals, init_depth50): def residual(params): x0, y0, depth, t0 params err [] for (sx, sy), ta in zip(stations_xy, arrivals): offset np.sqrt((sx-x0)**2(sy-y0)**2) tt self.travel_time(depth, offset) if np.isfinite(tt): err.append(ta - (t0tt)) return err res least_squares(residual, [0,0,init_depth,0], bounds([-100,-100,0,-10],[100,100,200,10])) return res.x if __name__ __main__: layers [5.0, 10.0, 15.0] # thickness meters vp [1500, 2000, 2500] # m/s model SphericalWaveInLayers(layers, vp) tt model.travel_time(source_depth25, offset50) print(fTravel time for offset 50m: {tt:.4f}s) stations np.array([[10,10],[20,30],[40,0],[60,50]]) arrivals np.array([0.112, 0.135, 0.142, 0.167]) src model.wavefront_inversion(stations, arrivals, init_depth30) print(fInverted source: x{src[0]:.1f}m, y{src[1]:.1f}m, depth{src[2]:.1f}m, t0{src[3]:.4f}s)
采空与多层岩土体中球面波传播规律与正反演方法【附代码】
发布时间：2026/5/21 15:13:01
✨ 长期致力于矿震、微震、球面波、点源与纵波、点涡与横波、采空区、正反演、震源定位、多层介质研究工作擅长数据搜集与处理、建模仿真、程序编写、仿真设计。✅ 专业定制毕设、代码✅如需沟通交流点击《获取方式》1含采空区岩土体的球面波散射解析模型与敲击试验验证基于弹性波动理论将采空区等效为圆柱形孔洞建立二维平面应变下的纵波入射散射场表达式。采用波函数展开法将入射波、散射波均用Hankel函数表示并在孔洞边界施加应力自由条件求解未知系数。推导出孔洞后方特定区域的加速度时程曲线发现存在指向波源的负加速度反转区。设计带孔有机玻璃板敲击试验使用加速度计阵列测量。在孔洞直径30mm、板厚8mm条件下孔后2倍孔径处监测到负向峰值达入射波峰值35%的反转信号与理论预测误差小于8%。利用有限元仿真分析采空区群的影响当空区间距小于直径时发生波导效应能量沿空区间柱体集中夹岩柱水平加速度可放大4倍。2多层水平介质中球面波前正演与震源深度几何平均定位法基于惠更斯原理和斯涅耳定律推导点源在N层水平介质中传播的走时曲线显式表达式。建立双曲线型波前曲面方程给出各层界面折射点的迭代求解方法:利用等时线从台站逆推至最后一层界面然后逐层向上计算。提出几何平均法确定震源深度:由两个不同偏移距台站的纵波初至时差构造关于深度的二次方程取正根作为深度估计。在二层混凝土试验中震源深度20cm几何平均法计算深度20.8cm误差4%。进一步结合发震时刻反演提出正反演联用法先假设震源水平坐标为零扫描深度得到时刻再反演水平偏移。该方法对倾斜界面同样适用通过坐标旋转将倾斜层转换成等效水平层。3任意观测系统下的波前正演定位算法与鲁棒性分析将发震时刻作为未知量加入非线性方程组利用各台站之间的到时差构造新的观测方程。采用信赖域反射算法求解最小二乘问题目标函数为理论走时与观测到时差的平方和。每次迭代中根据当前震源位置和发震时刻计算各个台站的理论走时并通过数值方法求解球面波在多层介质中的射线路径。为提高定位稳定性引入M估计器权重函数降低异常到时的影响采用Huber函数阈值设为1.5倍残差中位数。在模拟采空区扰动的环境下对10个随机台站进行定位测试传统Geiger法平均误差7.2m波前正演法降至2.3m。将该方法应用到某煤矿微震监测数据定位事件与采掘面位置吻合度提升39%。import numpy as np from scipy.optimize import least_squares class SphericalWaveInLayers: def __init__(self, layer_thick, layer_vp): self.nlayers len(layer_thick) self.h np.array(layer_thick) self.vp np.array(layer_vp) def travel_time(self, source_depth, offset): # ray tracing using Snells law if source_depth 0: source_depth 0 cum_h np.cumsum(self.h) if source_depth cum_h[0]: src_layer 0; z_in_layer source_depth else: src_layer np.searchsorted(cum_h, source_depth) z_in_layer source_depth - (cum_h[src_layer-1] if src_layer0 else 0) # simplified: constant ray parameter p p_guess 1.0/self.vp[0] def f(p): tt 0.0 x 0.0 for i in range(src_layer, self.nlayers): if isrc_layer: dz self.h[i] - z_in_layer else: dz self.h[i] if self.vp[i]*p 1: return 1e9 angle np.arcsin(self.vp[i]*p) tt dz/(self.vp[i]*np.cos(angle)) x dz*np.tan(angle) return x - offset from scipy.optimize import brentq try: p_opt brentq(f, 1e-6, 1.0/self.vp[-1]) tt 0.0 for i in range(src_layer, self.nlayers): if isrc_layer: dz self.h[i]-z_in_layer else: dz self.h[i] angle np.arcsin(self.vp[i]*p_opt) tt dz/(self.vp[i]*np.cos(angle)) return tt except: return np.inf def wavefront_inversion(self, stations_xy, arrivals, init_depth50): def residual(params): x0, y0, depth, t0 params err [] for (sx, sy), ta in zip(stations_xy, arrivals): offset np.sqrt((sx-x0)**2(sy-y0)**2) tt self.travel_time(depth, offset) if np.isfinite(tt): err.append(ta - (t0tt)) return err res least_squares(residual, [0,0,init_depth,0], bounds([-100,-100,0,-10],[100,100,200,10])) return res.x if __name__ __main__: layers [5.0, 10.0, 15.0] # thickness meters vp [1500, 2000, 2500] # m/s model SphericalWaveInLayers(layers, vp) tt model.travel_time(source_depth25, offset50) print(fTravel time for offset 50m: {tt:.4f}s) stations np.array([[10,10],[20,30],[40,0],[60,50]]) arrivals np.array([0.112, 0.135, 0.142, 0.167]) src model.wavefront_inversion(stations, arrivals, init_depth30) print(fInverted source: x{src[0]:.1f}m, y{src[1]:.1f}m, depth{src[2]:.1f}m, t0{src[3]:.4f}s)
