CN104931363B - 节理岩体变形模量测试方法 - Google Patents

Abstract

本发明涉及节理岩体变形模量测试技术，在建立黏弹性介质中自由界面入射波分离算法的基础上，提出了节理岩体变形模量的测试方法。本发明包含六个步骤：节理岩体变形模量测试；建立黏弹性介质中应力波在自由界面的反射系数计算公式；根据测点波形分离得到入射P波和SV波；球面发散的振幅补偿；计算节理岩体的弹性模量和黏性系数；计算节理岩体的变形模量。采用本发明测试获得不同峰值频率动载荷加载作用下岩体的变形模量，当峰值频率较小时，变形模量为静态变形模量；当峰值频率较高时，变形模量为动态变形模量，便于根据岩体工程所受的动载特性选择相应的变形模量进行相关的力学分析，具有操作简单和测试成本低的优点。

Description

节理岩体变形模量测试方法

技术领域

本发明涉及节理岩体变形模量测试技术，在建立黏弹性介质中自由界面入射波分离算法的基础上，提出了节理岩体变形模量的测试方法。

背景技术

岩体变形模量测试方法分为静态测试和动态测试两类，其中动态测试法包括共振法和波传播法。受岩体结构的影响，其力学参数具有多尺度特性，为更准确获取岩体力学参数，通常需要现场测试岩体力学参数。

现场静态测试岩体变形模量的方法主要有：平板千斤顶法，包括刚性和柔性承压板加载法；试洞法，包括水压加载和压力枕加载；钻孔法，包括全断面柔性加压钻孔膨胀计法、部分断面刚性或柔性加压钻孔弹模计法。

由于岩体属非均质、非线性天然材料，各种方法测定结果的差异很大。试洞法加压区大，代表性较好，但费用特高；平板法次之，费用较高，测点很少，且受表面松动影响，测试结果偏低；钻孔法较灵活，费用最低，特点是可利用地质勘探孔进行测试，覆盖面较大，测试深度较深，数据较多，通过统计分析可以获得最接近现场特性的数据。

共振法是基于结构的共振规律测试岩体力学参数，由于复杂结构的共振规律难以把握，因此通常将岩块加工成杆件，建立岩石杆件的共振方程，基于杆件的共振规律测试岩体力学参数，主要是测试岩石的弹性模量和阻尼系数。该方法具有测试成本低、测试结果能够反映介质的动力特性，该方法的不足之处是不便于大规模测试和节理发育的岩石难以加工成杆件等。波传播法是通过研究应力波在岩体中的传播规律并建立运动参数与岩体模量的关系，这种方法具有测试迅速、测试成本低，便于大范围、大规模测试等优点，但实测值大大高于静力法测值。由于设计需要静态值，所以许多研究人员致力于研究动、静两种测量结果的相关关系。动态测试测试结果大大高于静态测试结果，其主要原因是，岩体具有一定的粘滞性，在动力载荷作用下，岩体变形具有明显的率相关性，即不同的加载频率，岩体变形模量也不同。

发明内容

本发明的目的是提供一种节理岩体变形模量测试方法。

本发明的技术方案：一种节理岩体变形模量测试方法，包括以下步骤：

第一步：节理岩体变形模量测试

节理岩体变形模量测试包括两个方面的工作，

(1)现场选取岩块，将岩块加工成标准岩样测试密度，采用材料力学试验机测试岩样的泊松比；

(2)在测试区域内打一个钻孔，孔径为30～100mm，孔深范围5～40m，孔深由测试范围确定；在地表布置加速度或速度传感器，传感器个数大于3；孔内爆破产生冲击荷载，药量为40～100g，采集仪记录传感器的振动波形；

第二步：建立黏弹性介质中应力波在自由界面的反射系数计算公式

根据岩体的力学特性，选择合适的黏弹性本构关系描述岩体的变形特性，建立应力波在自由界面的反射系数计算公式；SV波入射至地表时，可能存在掠入射情况，为此该算法分两种情况，即SV波的入射角不超过临界入射角和SV波的入射角大于临界入射角；

(1)当SV波的入射角不超过临界入射角

当SV波的入射角不超过临界入射角时，应力波在自由界面的反射系数为：

式(1a)、(1b)、(1c)和(1d)中R_PP和R_PS分别为P波入射时反射P波和反射SV波的反射系数，R_SP和R_SS分别为SV波入射时反射P波和反射SV波的反射系数，和为岩体变形参数，当入射波是P波时，参数K、d_α和d_β的计算公式分别为：

式(2a)、(2b)和(2c)中ω为圆频率，θ为入射角度，ρ为岩体(岩样)密度，和为岩体变形参数；当入射波是SV波时，参数K、d_α和d_β的计算公式分别为：

式(3a)、(3b)和(3c)中ω为圆频率，θ为入射角度，ρ为岩体(岩样)密度，和为岩体变形参数；

(2)当SV波入射角大于临界入射角

当SV波的入射角大于临界入射角时，应力波在自由界面的反射系数R_PP与式(1a)相同，R_PS与式(1b)相同，R_SP和R_SS为：

R_SS＝-exp(i2γ) (4b)

式(4a)和(4b)中i为虚单位，θ为入射角度，参数s、γ和D的计算公式为：

式(5a)、(5b)和(5c)中ν为泊松比，θ为入射角度；

第三步：根据测点波形分离得到入射P波和SV波

在地表记录的波形不仅包含入射波和反射波，当入射角不为0时，应力波在地表发生波形转换；为使得到准确的入射波，需要进行波场分离，将第一步的实测波形代入下面式(6a)和式(6b)分离得到入射波P波和SV波；根据各测点起跳时差，计算岩体的P波和SV波波速，计算SV波的临界入射角，P波和SV波的分离算法为：

式(6a)和(6b)中P_j(t)和SV_j(t)分别为第j测点的入射P波和SV波，u_xj(t)和u_zj(t)分别为第j测点记录的x方向即沿测线方向和z方向即垂直测线且垂直地表方向的振动波形，B₁、B₂、C₁和C₂参数的计算公式为：

C₁＝(1+R_PP)sinθ-DR_PScosθ₁ (7c)

C₂＝(-1+R_PP)cosθ+DR_PSsinθ₁ (7d)

式(7a)、(7b)、(7c)和(7d)中R_PP和R_PS分别为P波入射时反射P波和反射SV波的反射系数，R_SP和R_SS分别为SV波入射时反射P波和反射SV波的反射系数，反射系数R_PP、R_PS、R_SP和R_SS由式(1a)、(1b)、(1c)、(1d)、(4a)和(4b)计算，θ₁为P波入射时产生反射SV波的反射角，θ₂为SV波入射产生反射P波的反射角，θ₁、θ₂和D的计算公式为：

θ₁＝arcsin(D^-1sinθ) (8a)

式(8a)、(8b)和(8c)中ν为泊松比，θ为入射角度；

第四步：球面发散的振幅补偿

应力波在粘弹性介质中传播时，振幅会发生明显衰减，振幅衰减包含两部分，即几何衰减和非弹性衰减；为准确计算节理岩体的力学参数，需进行几何衰减的振幅补偿，即将第三步分离得到各测点的入射P波乘以补偿系数；补偿系数ζ为：

式(9)中ζ_j为第j个测点的振幅补偿系数，下标j为测点编号，r_j为第j个测点与振源间的距离；

第五步：计算节理岩体的弹性模量和黏性系数

现对振幅补偿处理后的各测点P波进行带通滤波，带宽为100～400Hz，滤波后各测点的振动波形分解为一系列子波；读取所有测点的各阶子波最大振幅、峰值频率，并代入式(10)，采用最小二乘法计算岩体弹性模量和黏性系数；

A_j＝A₁exp[α(r₁-r_j)](j≥2) (10)

式(10)中A_j为第j个测点入射子波的最大振幅，A₁为第1个测点子波的最大振幅，r_j为第j个测点与振源间的距离，r₁为第1个测点与振源间的距离，α为子波振幅衰减系数，根据第二步中选择的黏弹性本构关系确定α的计算公式；

第六步：计算节理岩体的变形模量

根据第五步计算得到的弹性模量和黏性系数计算节理岩体的变形模量，计算公式为：

式中为节理岩体的变形模量，E_r为节理岩体的弹性模量，E_i为节理岩体的损耗模量；采用式(11)计算得到岩体变形模量便是本发明的测试结果。

本发明的有益效果：

(1)采用本发明测试获得不同峰值频率动载荷加载作用下岩体的变形模量，当峰值频率较小时，变形模量为静态变形模量；当峰值频率较高时，变形模量为动态变形模量，便于根据岩体工程所受的动载特性选择相应的变形模量进行相关的力学分析。

(2)本发明在充分考虑波形转换的基础上，提出了根据测点振动波形计算入射P波和SV波的算法，提高了测试结果的精度，在低频条件下，采用本发明的测试结果与钻孔弹模仪测试结果非常接近，误差约为3.8％，本发明测试结果是有效的。

具体实施方式

采用本发明在国内某露天开采的矿山进行了不公开测试，具体实施步骤如下。

第一步：节理岩体变形模量测试

节理岩体变形模量测试包括两个方面的工作，具体如下：

(1)现场选取岩块，将岩块加工成标准岩样，测试得到岩块密度为2520kg/m³，采用岩石力学试验系统(RMT-150C)测试岩样的泊松比为0.22。

(2)在某露天开采形成的临时台阶上进行节理岩体的变形模量动力测试试验，台阶宽度约为72m，用潜孔钻向下垂直钻孔，孔径为70mm，孔深为12.0m，在台阶地表上布置一条测线，沿测线布置5个测点，用石膏粉将加速度传感器粘结在测点上。传感器为二分量加速度传感器，其中x方向与测线方向一致，z方向与测线垂直。传感器参数为：电荷灵敏度范围为10.284～14.147pC/(m.s.²)，频响为0.2Hz～5.0kHz。在钻孔中安放50g乳化炸药，炸药与地表距离为10.0m，引爆炸药产生冲击荷载，采用INV306振动采集仪记录测点振动波形。

第二步：建立黏弹性介质中应力波在自由界面的反射系数计算公式

根据岩体的力学特性，采用Kelvin模型描述试验区域节理岩体的变形特性，建立应力波在自由界面的反射系数计算公式。SV波入射至地表时，可能存在掠入射情况，为此该算法分两种情况，即SV波的入射角不超过临界入射角和SV波的入射角大于临界入射角。

(1)当SV波的入射角不超过临界入射角

当SV波的入射角不超过临界入射角时，应力波在自由界面的反射系数为：

式(1a)、(1b)、(1c)和(1d)中R_PP和R_PS分别为P波入射时反射P波和SV波的反射系数，R_SP和R_SS分别为SV波入射时反射P波和SV波的反射系数，和为岩体变形参数，当入射波是P波时，参数K、d_α和d_β的计算公式分别为：

式(2a)、(2b)和(2c)中ω为圆频率，θ为入射角度，ρ为岩体(岩样)密度，和为岩体变形参数；当入射波是SV波时，参数K、d_α和d_β的计算公式分别为：

式(3a)、(3b)和(3c)中ω为圆频率，θ为入射角度，ρ为岩体(岩样)密度，和为岩体的岩体变形参数。

对于Kelvin模型，变形参数和为：

式(4a)和(4b)中ω为圆频率，i为虚单位，E_r为岩体弹性模量，η为岩体黏性系数，ν为泊松比。

(2)SV波入射角大于临界入射角

当SV波的入射角大于临界入射角时，应力波在自由界面的反射系数R_PP和R_PS与式(1a)相同，R_PS与式(1b)相同，R_SP和R_SS为：

R_SS＝-exp(i2γ) (5b)

式(5a)和(5b)中R_SP和R_SS分别为SV波入射时反射P波和SV波的反射系数，i为虚单位，参数s、γ和D的计算公式为：

式(6a)、(6b)和(6c)中ν为泊松比，θ为入射角度。

第三步：根据测点波形分离得到入射P波和SV波

在地表记录的波形不仅包含入射波和反射波，当入射角不为0时，应力波在地表发生波形转换。为使得到更准确的入射波，需要进行波场分离，将第一步的实测波形代入式(7a)和式(7b)分离得到入射波P波和SV波。根据各测点起跳时差，计算岩体的P波和SV波波速，计算SV波的临界入射角，提出P波和SV波的分离算法为：

式(7a)和(7b)中P_j(t)和SV_j(t)分别为第j测点的入射P波和SV波，u_xj(t)和u_zj(t)分别为第j测点记录的x方向(沿测线方向)和z方向(垂直测线且垂直地表的方向)的振动波形，B₁、B₂、C₁和C₂参数的计算公式为：

C₁＝(1+R_PP)sinθ-DR_PScosθ₁ (8c)

C₂＝(-1+R_PP)cosθ+DR_PSsinθ₁ (8d)

式(8a)、(8b)、(8c)和(8d)中R_PP和R_PS分别为P波入射时反射P波和反射SV波的反射系数，R_SP和R_SS分别为SV波入射时反射P波和反射SV波的反射系数，反射系数R_PP、R_PS、R_SP和R_SS由式(1a)、(1b)、(1c)、(1d)、(4a)和(4b)计算，θ₁为P波入射时产生反射SV波的反射角，θ₂为SV波入射产生反射P波的反射角，θ₁、θ₂和D的计算公式为：

θ₁＝arcsin(D^-1sinθ) (9a)

式(9a)、(9b)和(9c)中ν为泊松比，θ为入射角度。

第四步：球面发散的振幅补偿

面波在粘弹性介质中传播时，振幅会发生明显衰减，振幅衰减包含两部分，即几何衰减和非弹性衰减。为准确计算节理岩体的力学参数，需进行几何衰减的振幅补偿，即将第三步分离得到各测点的入射P波乘以补偿系数。补偿系数ζ为：

式(10)中ζ为补偿系数，下标j为测点编号，r为测点与振源间的距离。

第五步：计算节理岩体的弹性模量和黏性系数

现对振幅补偿处理后的各测点P波进行带通滤波，带宽为200Hz，滤波后各测点的振动波形分解为一系列子波。读取所有测点的各阶子波最大振幅、峰值频率，并代入式(10)，采用最小二乘法计算岩体弹性模量和黏性系数。

A_j＝A₁exp[α(r₁-r_j)](j≥2) (11)

式(11)中A_j为第j个测点入射子波的最大振幅，A₁为第1个测点子波的最大振幅，r_j为第j个测点与振源间的距离，r₁为第1个测点与振源间的距离，α为子波衰减系数，对于Kelvin模型，α的计算公式为：

式(12)中ω_f为子波的峰值圆频率，ρ为岩体(岩样)密度，E_r为岩体弹性模量，η为岩体黏性系数，ν为泊松比。

第六步：计算节理岩体的变形模量

根据第五步计算得到的弹性模量和黏性系数计算节理岩体的变形模量。岩体变形模量的计算公式为：

式(13)中为节理岩体的变形模量，E_r为岩体弹性模量，ω_f为子波的峰值圆频率，η为岩体黏性系数，对于Kelvin模型，损耗模量E_i＝ηω_f。

采用式(13)计算得到岩体变形模量便是本发明的测试结果，测试结果见表1，为验证本发明测试结果的有效性，采用由欧美大地仪器设备中国有限公司研制开发的Geo-ExpertsGY-90型钻孔弹模仪现场测试节理岩体的弹性模量，依据大循环压力变形曲线计算岩体弹性模量E＝(30-12)/(3.5605-1.5285)×178.3/1000＝1.58GPa。在低频条件下，采用本发明的测试结果与钻孔弹模仪测试结果非常接近，误差为3.8％，表明本发明测试结果是有效的。

表1基于不同峰值频率子波计算节理岩体的弹性模量

Claims (1)

1.一种节理岩体变形模量测试方法，包括以下步骤：

第一步：节理岩体变形模量测试

节理岩体变形模量测试包括两个方面的工作，

(1)现场选取岩块，将岩块加工成标准岩样测试密度，采用材料力学试验机测试岩样的泊松比；

(2)在测试区域内打一个钻孔，孔径为30～100mm，孔深范围5～40m，孔深由测试范围确定；在地表布置加速度或速度传感器，传感器个数大于3；孔内爆破产生冲击荷载，药量为40～100g，采集仪记录传感器的振动波形；

第二步：建立黏弹性介质中应力波在自由界面的反射系数计算公式

根据岩体的力学特性，选择合适的黏弹性本构关系描述岩体的变形特性，建立应力波在自由界面的反射系数计算公式；SV波入射至地表时，可能存在掠入射情况，为此该算法分两种情况，即SV波的入射角不超过临界入射角和SV波的入射角大于临界入射角；

(1)当SV波的入射角不超过临界入射角

当SV波的入射角不超过临界入射角时，应力波在自由界面的反射系数为：

$R_{PP}=\frac{4K^2\stackrel{\‾}{\mathrm{\μ}}{d}_{\mathrm{\β}}{d}_{\mathrm{\α}}-\[(\stackrel{\‾}{\mathrm{\λ}}+2\stackrel{\‾}{\mathrm{\μ}})d_{\mathrm{\α}}^2+\stackrel{\‾}{\mathrm{\λ}}{K}^2\]({d_{\mathrm{\β}}}^2-K^2)}{4K^2\stackrel{\‾}{\mathrm{\μ}}{d}_{\mathrm{\α}}{d}_{\mathrm{\β}}+\[(\stackrel{\‾}{\mathrm{\λ}}+2\stackrel{\‾}{\mathrm{\μ}})d_{\mathrm{\α}}^2+\stackrel{\‾}{\mathrm{\λ}}{K}^2\]}---\left(1a\right)$

$R_{PS}=\frac{-4{\mathrm{Kd}}_{\mathrm{\α}}\[(\stackrel{\‾}{\mathrm{\λ}}+2\stackrel{\‾}{\mathrm{\μ}})d_{\mathrm{\α}}^2+\stackrel{\‾}{\mathrm{\λ}}{K}^2\]}{4K^2\stackrel{\‾}{\mathrm{\μ}}{d}_{\mathrm{\α}}{d}_{\mathrm{\β}}+\[(\stackrel{\‾}{\mathrm{\λ}}+2\stackrel{\‾}{\mathrm{\μ}})d_{\mathrm{\α}}^2+\stackrel{\‾}{\mathrm{\λ}}{K}^2\](d_{\mathrm{\β}}^2-K^2)}---\left(1b\right)$

$R_{SP}=\frac{4\stackrel{\‾}{\mathrm{\μ}}{d}_{\mathrm{\β}}K(d_{\mathrm{\β}}^2-K^2)}{4K^2\stackrel{\‾}{\mathrm{\μ}}{d}_{\mathrm{\α}}{d}_{\mathrm{\β}}+\[(\stackrel{\‾}{\mathrm{\λ}}+2\stackrel{\‾}{\mathrm{\μ}})d_{\mathrm{\α}}^2+\stackrel{\‾}{\mathrm{\λ}}{K}^2\](d_{\mathrm{\β}}^2-K^2)}---\left(1c\right)$

$R_{SS}=\frac{4K^2\stackrel{\‾}{\mathrm{\μ}}{d}_{\mathrm{\β}}{d}_{\mathrm{\α}}-\[(\stackrel{\‾}{\mathrm{\λ}}+2\stackrel{\‾}{\mathrm{\μ}})d_{\mathrm{\α}}^2+\stackrel{\‾}{\mathrm{\λ}}{K}^2\](d_{\mathrm{\β}}^2-K^2)}{4K^2\stackrel{\‾}{\mathrm{\μ}}{d}_{\mathrm{\α}}{d}_{\mathrm{\β}}+\[(\stackrel{\‾}{\mathrm{\λ}}+2\stackrel{\‾}{\mathrm{\μ}})d_{\mathrm{\α}}^2+\stackrel{\‾}{\mathrm{\λ}}{K}^2\]}---\left(1d\right)$

式(1a)、(1b)、(1c)和(1d)中R_PP和R_PS分别为P波入射时反射P波和反射SV波的反射系数，R_SP和R_SS分别为SV波入射时反射P波和反射SV波的反射系数，和为岩体变形参数，当入射波是P波时，参数K、d_α和d_β的计算公式分别为：

$K=\mathrm{\ω}\sqrt{\frac{\mathrm{\ρ}}{\stackrel{\‾}{\mathrm{\λ}}+2\stackrel{\‾}{\mathrm{\μ}}}}sin\mathrm{\θ}---\left(2a\right)$

$d_{\mathrm{\α}}=\mathrm{\ω}\sqrt{\frac{\mathrm{\ρ}}{\stackrel{\‾}{\mathrm{\λ}}+2\stackrel{\‾}{\mathrm{\μ}}}}cos\mathrm{\θ}---\left(2b\right)$

$d_{\mathrm{\β}}=\mathrm{\ω}\sqrt{\frac{\mathrm{\ρ}}{\stackrel{\‾}{\mathrm{\λ}}+2\stackrel{\‾}{\mathrm{\μ}}}}\sqrt{\frac{\stackrel{\‾}{\mathrm{\λ}}+2\stackrel{\‾}{\mathrm{\μ}}}{\stackrel{\‾}{\mathrm{\μ}}}-{\sin }^2\mathrm{\θ}}---\left(2c\right)$

式(2a)、(2b)和(2c)中ω为圆频率，θ为入射角度，ρ为岩体(岩样)密度，和为岩体变形参数；当入射波是SV波时，参数K、d_α和d_β的计算公式分别为：

$K=\mathrm{\ω}\sqrt{\frac{\mathrm{\ρ}}{\stackrel{\‾}{\mathrm{\μ}}}}sin\mathrm{\θ}---\left(3a\right)$

$d_{\mathrm{\α}}=\mathrm{\ω}\sqrt{\frac{\mathrm{\ρ}}{\stackrel{\‾}{\mathrm{\μ}}}}\sqrt{\frac{\stackrel{\‾}{\mathrm{\μ}}}{\stackrel{\‾}{\mathrm{\λ}}+2\stackrel{\‾}{\mathrm{\μ}}}-{\sin }^2\mathrm{\θ}}---\left(3b\right)$

$d_{\mathrm{\β}}=\mathrm{\ω}\sqrt{\frac{\mathrm{\ρ}}{\stackrel{\‾}{\mathrm{\μ}}}}cos\mathrm{\θ}---\left(3c\right)$

式(3a)、(3b)和(3c)中ω为圆频率，θ为入射角度，ρ为岩体(岩样)密度，和为岩体变形参数；

(2)当SV波入射角大于临界入射角

当SV波的入射角大于临界入射角时，应力波在自由界面的反射系数R_PP与式(1a)相同，R_PS与式(1b)相同，R_SP和R_SS为：

$R_{SP}={\left\{\frac{{\tan }^2\mathrm{\θ}{\[D^2(1-s^2)-2\]}^2}{{\tan }^2\mathrm{\θ}{\[D^2(1-s^2)-2\]}^4+16s^2}\right\}}^{\frac{1}{2}}\exp \left(i\mathrm{\γ}\right)---\left(4a\right)$

R_SS＝-exp(i2γ) (4b)

式(4a)和(4b)中i为虚单位，θ为入射角度，参数s、γ和D的计算公式为：

$\mathrm{\γ}=arctan\left(\frac{4s}{tan\mathrm{\θ}{\[\frac{2-2v}{1-2v}(1-s^2)-2\]}^2}\right)---\left(5a\right)$

$s={\[1-\frac{1-2v}{(2-2v){\sin }^2\mathrm{\θ}}\]}^{\frac{1}{2}}---\left(5b\right)$

$D={\left(\frac{2-2v}{1-2v}\right)}^{0.5}---\left(5c\right)$

式(5a)、(5b)和(5c)中ν为泊松比，θ为入射角度；

第三步：根据测点波形分离得到入射P波和SV波

在地表记录的波形不仅包含入射波和反射波，当入射角不为0时，应力波在地表发生波形转换；为使得到准确的入射波，需要进行波场分离，将第一步的实测波形代入下面式(6a)和式(6b)分离得到入射波P波和SV波；根据各测点起跳时差，计算岩体的P波和SV波波速，计算SV波的临界入射角，P波和SV波的分离算法为：

$P_j\left(t\right)=-\frac{B_2{u}_{xj}\left(t\right)+B_1{u}_{zj}\left(t\right)}{B_1{C}_2-B_2{C}_1}---\left(6a\right)$

${\mathrm{SV}}_j\left(t\right)=\frac{C_2{u}_{xj}\left(t\right)-C_1{u}_{zj}\left(t\right)}{B_1{C}_2-B_2{C}_1}---\left(6b\right)$

式(6a)和(6b)中P_j(t)和SV_j(t)分别为第j测点的入射P波和SV波，u_xj(t)和u_zj(t)分别为第j测点记录的x方向即沿测线方向和z方向即垂直测线且垂直地表方向的振动波形，B₁、B₂、C₁和C₂参数的计算公式为：

$B_1=(1-R_{SS})cos\mathrm{\θ}+\frac{1}{D}{R}_{SP}{\mathrm{sin\θ}}_2---\left(7a\right)$

$B_2=(1+R_{SS})sin\mathrm{\θ}+\frac{1}{D}{R}_{SP}{\mathrm{cos\θ}}_2---\left(7b\right)$

C₁＝(1+R_PP)sinθ-DR_PScosθ₁ (7c)

C₂＝(-1+R_PP)cosθ+DR_PSsinθ₁ (7d)

式(7a)、(7b)、(7c)和(7d)中R_PP和R_PS分别为P波入射时反射P波和反射SV波的反射系数，R_SP和R_SS分别为SV波入射时反射P波和反射SV波的反射系数，反射系数R_PP、R_PS、R_SP和R_SS由式(1a)、(1b)、(1c)、(1d)、(4a)和(4b)计算，θ₁为P波入射时产生反射SV波的反射角，θ₂为SV波入射产生反射P波的反射角，θ₁、θ₂和D的计算公式为：

θ₁＝arcsin(D^-1sinθ) (8a)

${\mathrm{\θ}}_2=\left\{\begin{array}{cc}arcsin\left(Dsin\mathrm{\θ}\right)& \mathrm{\θ}\≤\arcsin D^{-1}\\ \frac{\mathrm{\π}}{2}& \mathrm{\θ}>\arcsin D^{-1}\end{array}\right.---\left(8b\right)$

$D={\left(\frac{2-2v}{1-2v}\right)}^{0.5}---\left(8c\right)$

式(8a)、(8b)和(8c)中ν为泊松比，θ为入射角度；

第四步：球面发散的振幅补偿

应力波在粘弹性介质中传播时，振幅会发生明显衰减，振幅衰减包含两部分，即几何衰减和非弹性衰减；为准确计算节理岩体的力学参数，需进行几何衰减的振幅补偿，即将第三步分离得到各测点的入射P波乘以补偿系数；补偿系数ζ为：

${\mathrm{\ζ}}_j=\frac{r_j}{r_1},(j\≥2)---\left(9\right)$

式(9)中ζ_j为第j个测点的振幅补偿系数，下标j为测点编号，r_j为第j个测点与振源间的距离；

第五步：计算节理岩体的弹性模量和黏性系数

现对振幅补偿处理后的各测点P波进行带通滤波，带宽为100～400Hz，滤波后各测点的振动波形分解为一系列子波；读取所有测点的各阶子波最大振幅、峰值频率，并代入式(10)，采用最小二乘法计算岩体弹性模量和黏性系数；

A_j＝A₁exp[α(r₁-r_j)](j≥2) (10)

式(10)中A_j为第j个测点入射子波的最大振幅，A₁为第1个测点子波的最大振幅，r_j为第j个测点与振源间的距离，r₁为第1个测点与振源间的距离，α为子波振幅衰减系数，根据第二步中选择的黏弹性本构关系确定α的计算公式；

第六步：计算节理岩体的变形模量

根据第五步计算得到的弹性模量和黏性系数计算节理岩体的变形模量，计算公式为：

$\stackrel{\‾}{E}=\sqrt{E_r^2+E_i^2}---\left(11\right)$

式中为节理岩体的变形模量，E_r为节理岩体的弹性模量，E_i为节理岩体的损耗模量。