弹性理论主要研究物体受力与形变的关系, 其波动方程由本构方程、运动平衡微分方程和几何方程共同建立。在弹性各向同性介质中, 二维一阶速度-应力弹性波动方程为：

$\left\{ \begin{align} & {{\rho }_{j}}\frac{\partial u}{\partial t}=\frac{\partial {{\tau }_{xx}}}{\partial x}+\frac{\partial {{\tau }_{xz}}}{\partial z} \\ & {{\rho }_{j}}\frac{\partial v}{\partial t}=\frac{\partial {{\tau }_{xz}}}{\partial x}+\frac{\partial {{\tau }_{zz}}}{\partial z} \\ & \frac{\partial {{\tau }_{xx}}}{\partial t}=\left( \lambda +2\mu \right)\frac{\partial u}{\partial x}+\lambda \frac{\partial v}{\partial z} \\ & \frac{\partial {{\tau }_{zz}}}{\partial t}=\lambda \frac{\partial u}{\partial x}+\left( \lambda +2\mu \right)\frac{\partial v}{\partial z} \\ & \frac{\partial {{\tau }_{xz}}}{\partial t}=\mu \frac{\partial v}{\partial x}+\mu \frac{\partial u}{\partial z} \\ \end{align} \right.$ （1）

在时间上, 采用二阶精度的时间差分方程解的近似值, 在空间采用四阶差分方程解的近似值, 利用交错网格格式, 取

$\left\{ \begin{align} & x=i\Delta x \\ & z=j\Delta z \\ & t=k\Delta t \\ \end{align} \right.$ （2）

可推导出波动方程的差分格式：

$\left\{ \begin{align} & A_{i, j}^{k+{1}/{2}\; }=A_{i, j}^{k-{1}/{2}\; }+\frac{\Delta t}{{{\rho }_{i, j}}}\frac{1}{\Delta x}\sum\limits_{n=1}^{2}{C_{n}^{\left( 2 \right)}\left[ D_{i+{\left( 2n-1 \right)}/{2}\; , j}^{k}-D_{i-{\left( 2n-1 \right)}/{2}\; , j}^{k} \right]} \\ & +\frac{\Delta t}{{{\rho }_{i, j}}}\frac{1}{\Delta z}\sum\limits_{n=1}^{2}{C_{n}^{\left( 2 \right)}\left[ F_{i, j+{\left( 2n-1 \right)}/{2}\; }^{k}-F_{i, j-{\left( 2n-1 \right)}/{2}\; }^{k} \right]} \\ & B_{i, j}^{k+{1}/{2}\; }=B_{i, j}^{k-{1}/{2}\; }+\frac{\Delta t}{{{\rho }_{i, j}}}\frac{1}{\Delta x}\sum\limits_{n=1}^{2}{C_{n}^{\left( 2 \right)}\left[ F_{i+{\left( 2n-1 \right)}/{2}\; , j}^{k}-F_{i-{\left( 2n-1 \right)}/{2}\; , j}^{k} \right]} \\ & +\frac{\Delta t}{{{\rho }_{i, j}}}\frac{1}{\Delta z}\sum\limits_{n=1}^{2}{C_{n}^{\left( 2 \right)}\left[ E_{i, j+{\left( 2n-1 \right)}/{2}\; }^{k}-E_{i, j-{\left( 2n-1 \right)}/{2}\; }^{k} \right]} \\ & D_{i, j}^{k+1}=D_{i, j}^{k}+\left( \lambda +2\mu \right)\frac{\Delta t}{\Delta x}\sum\limits_{n=1}^{2}{C_{n}^{\left( 2 \right)}\left[ A_{i+{\left( 2n-1 \right)}/{2}\; , j}^{k}-A_{i-{\left( 2n-1 \right)}/{2}\; , j}^{k} \right]} \\ & +\lambda \frac{\Delta t}{\Delta z}\sum\limits_{n=1}^{2}{C_{n}^{\left( 2 \right)}\left[ B_{i, j+{\left( 2n-1 \right)}/{2}\; }^{k}-B_{i, j-{\left( 2n-1 \right)}/{2}\; }^{k} \right]} \\ & E_{i, j}^{k+1}=E_{i, j}^{k}+\lambda \frac{\Delta t}{\Delta x}\sum\limits_{n=1}^{2}{C_{n}^{\left( 2 \right)}\left[ A_{i+{\left( 2n-1 \right)}/{2}\; , j}^{k}-A_{i-{\left( 2n-1 \right)}/{2}\; , j}^{k} \right]} \\ & +\left( \lambda +2\mu \right)\frac{\Delta t}{\Delta z}\sum\limits_{n=1}^{2}{C_{n}^{\left( 2 \right)}\left[ B_{i, j+{\left( 2n-1 \right)}/{2}\; }^{k}-B_{i, j-{\left( 2n-1 \right)}/{2}\; }^{k} \right]} \\ & F_{i, j}^{k+1}=F_{i, j}^{k}+\mu \frac{\Delta t}{\Delta x}\sum\limits_{n=1}^{2}{C_{n}^{\left( 2 \right)}\left[ B_{i+{\left( 2n-1 \right)}/{2}\; , j}^{k}-B_{i-{\left( 2n-1 \right)}/{2}\; , j}^{k} \right]} \\ & +\mu \frac{\Delta t}{\Delta z}\sum\limits_{n=1}^{2}{C_{n}^{\left( 2 \right)}\left[ A_{i, j+{\left( 2n-1 \right)}/{2}\; }^{k}-A_{i, j-{\left( 2n-1 \right)}/{2}\; }^{k} \right]} \\ \end{align} \right.$ （3）

1.2.1 边界条件

在计算区域中左、右边界与实际条件基本相符, 按照常规反射边界来处理, 而上、下边界和4个边角点设置吸收条件。

1.2.2 初始条件

当t≤ 0时, 速度和位移为零, 即

$\left\{ \begin{align} & U\left( x, z, t \right)\left| _{t\le 0} \right.=0 \\ & V\left( x, z, t \right)\left| _{t\le 0} \right.=0 \\ & u\left( x, z, t \right)\left| _{t\le 0} \right.=0 \\ & v\left( x, z, t \right)\left| _{t\le 0} \right.=0 \\ \end{align} \right.$ （4）

1.2.3 振源条件

$U\left( x, z, 0 \right)=\left\{ \begin{matrix} {{R}_{\text{1}}}\left( {{x}_{0}}, {{z}_{0}}, 0 \right) & {} \\ 0 & {} \\ \end{matrix} \right.$ （5）

1.2.4收敛条件

分析误差传播和积累的情况, 是保证差分格式收敛性的基础。董良国等^[15]详细推导了交错网格高阶差分方程的稳定性条件。当介质为均匀各向同性时, 时间精度为二阶, 空间精度为四阶的稳定条件为：

$\frac{\Delta t}{\Delta x}\sqrt{v_{\text{p}}^{2}+v_{\text{s}}^{2}}\le {6}/{7}\; $ （6）

实验室超声波测试采用透射法进行测量。但由于页岩的非均质性及层理分布的不确定性, 给研究层理结构对超声波传播特性的影响带来了较大困难, 而且透射实验中岩心与探头接触不平整以及人为读取初始时间等也会造成误差。数值模拟基于被测岩心的纵向剖面建立模型, 以超声波探头的激发信号为振源, 基于波动理论推导出一阶位移-应力弹性波方程组, 并开展交错网格半程计算, 用Matlab编程实现超声波透射实验的模拟。

1.3.1波形对比验证

为了使模型与超声波透射实验相统一, 将探头频率设为250 kHz。其中入射波信号显示为初始波形, 当入射信号透过岩心到达接收端时, 实验室可获取透射信号, 显示为实验波形, 而数值计算的结果显示为模拟波形（见图1）。波形对比结果显示, 模拟波形与实验波形到达接收端的时间均为0.9× 10^-5 s左右, 电压均为-5~5 V, 二者范围相近（见图1）。

图1

	Figure Option ViewDownloadNew Window
	图1 标准页岩超声波透射实验与数值模拟结果对比图

将模拟波形和实验波形进行相关性计算, 将模拟波形信号设为X, 实验波形信号设为Y：

$r=\frac{\operatorname{cov}(X, Y)}{\sqrt{{{D}_{0}}(X)}\sqrt{{{D}_{0}}(Y)}}$ （7）

X和Y的标准差分别为：

$\begin{align} & {{D}_{0}}(X)={{E}_{0}}{{\left[ X-{{E}_{0}}(X) \right]}^{2}} \\ & {{D}_{0}}(Y)={{E}_{0}}{{\left[ Y-{{E}_{0}}(Y) \right]}^{2}} \\ \end{align}$ （8）

协方差：

$\operatorname{cov}(X, Y)={{E}_{0}}\left\{ \left[ X-{{E}_{0}}\left( X \right) \right]\left[ Y-{{E}_{0}}\left( Y \right) \right] \right\}$ （9）

通过相关性计算得到模拟波形与实验波形之间的相关系数为0.85, 吻合度较高, 所以该数值模拟方法可对超声波透射实验进行有效的模拟。

1.3.2数值计算验证

考虑到层理角度为页岩最显著的结构特征, 因此, 结合实际岩心建立类似不同层理角度的理想岩心模型, 在此基础上开展数值模拟计算。图2是渝东南地区下志留统龙马溪组露头岩心, 编号记为A₁、A₂、A₃、A₄。层理角度为层理方向与超声波传播方向的夹角, 其中A₁的层理角度为0° , A₂为30° , A₃为60° , A₄为90° （见图2）。为了进一步验证数值模拟方法的有效性, 本文还以A₁、A₂、A₃、A₄岩心为模板, 建立数值模型B₁、B₂、B₃、B₄（见图3）。在模型中对速度参数进行设置, 将超声波在层理中的速度值设置为1 900 m/s, 将其在骨架中的速度值设置为4 200 m/s, 该模型为理想的层理角度模型, 其角度与岩心A₁、A₂、A₃、A₄一一对应。在此基础上, 对超声波透射实验与数值模拟结果进行对比验证。

图2

	Figure Option ViewDownloadNew Window
	图2 渝东南地区下志留统龙马溪组露头岩心照片（据参考文献[9]修改）

图3

	Figure Option ViewDownloadNew Window
	图3 渝东南地区龙马溪组岩心数值模型（据参考文献[9]修改）

通过超声波透射实验完成对岩心A₁、A₂、A₃、A₄在不同探头频率条件下衰减系数的计算, 结果显示, 在同一层理角度下, 随着探头频率的增大, 衰减系数整体呈现出增大趋势（见图4a）。同时, 通过数值模拟也完成了对模型B₁、B₂、B₃、B₄在相同条件下的计算（见图4b）, 数值模拟计算得到的规律与超声波透射实验规律基本吻合^[16], 即在同一层理角度下, 衰减系数与探头频率均呈正相关。

图4

	Figure Option ViewDownloadNew Window
	图4 超声波衰减特性结果对比图（样品数4个）

上述数值计算是基于理想的层理角度模型。为了进一步验证上述数值模拟结果的可靠性, 利用数值图像处理技术提取出真实岩心的层理角度, 然后再开展不同层理角度的衰减系数数值计算。图5是渝东南地区下志留统龙马溪组露头岩心照片, 编号记为C₁、C₂、C₃、C₄, 其中C₁的层理角度为0° , C₂为30° , C₃为45° , C₄为90° （见图5）。

图5

	Figure Option ViewDownloadNew Window
	图5 渝东南地区下志留统龙马溪组露头岩心照片

以C₁、C₂、C₃、C₄岩心为研究对象, 对照片中的层理和基质骨架区域进行边缘提取、Hough直线检测^[17], 建立4个数值模拟模型, 编号分别记为D₁、D₂、D₃、D₄（见图6）, 其层理角度与C₁、C₂、C₃、C₄一一对应。由于数字图像是由矩形排列的像素点组成, 每个像素点是一个小正方形, 这些小正方形均可作为有限差分计算中的四边形单元, 故可将层理结构提取后的数字图像通过方形单元映射的方法可以很方便地转化成有限差分网格数据。首先, 每个单元的4个节点坐标通过图像实际尺寸与像素尺寸之间的比例关系, 转换为相应矢量空间的物理坐标, 整个图像就可转换为正方形有限差分网格, 此法即为方形单元映射法。然后根据每个像素点所属灰度值, 给每个单元赋予相应的参数, 从而把数字图像表征的层理结构引入到数值模型中。最后, 利用Matlab算法开展数值模拟计算, 并将计算结果与超声波透射实验结果对比。从衰减系数结果图来看, 数值模拟的衰减系数与层理角度呈正相关, 这与超声波透射实验中岩心的衰减系数变化规律基本一致（见图7）。

图6

	Figure Option ViewDownloadNew Window
	图6 渝东南地区下志留统龙马溪组岩心层理提取结果图

图7

	Figure Option ViewDownloadNew Window
	图7 实验结果与数值模拟结果图（样品数4个）

综上所述, 针对超声波透射实验成本高, 由于实验过程中岩心与探头接触不平整、人为数据读取等造成误差较大的现状^[18], 本文提出的页岩超声波透射实验的数值模拟方法与实际超声波透射实验结果吻合度较高, 因此, 可依托此数值模拟计算手段, 开展页岩层理结构声学特性研究。

页岩的显著特征是层理发育, 表征层理结构的参数主要包括层理厚度、层理密度和层理角度。为了开展与超声波透射实验条件相吻合的数值模拟计算, 本文将模型中的岩心柱塞高度设置为50mm, 直径为25mm, 页岩骨架中的纵波速度设为4200m/s, 横波速度设为2500m/s; 页岩层理中的纵波速度为1900m/s, 横波波速为1100m/s; 页岩骨架密度为2.6 g/cm³, 层理密度为2.3 条/cm²。提取实验室探头入射波作为初始波, 保持超声波传播方向不变, 将柱塞纵向剖面区域划分成长度为250mm、宽为125mm的网格, 空间网格间距为0.2mm, 时间步长为20ns。通过不断调整柱塞的层理模型参数来实现页岩结构的变化, 完成对超声波透射实验的模拟, 为探索层理结构对超声波特性的影响奠定基础。

2.1.1 层理厚度对超声波特性的影响

模型参数设置中, 我们将平直且与层面平行的层理称为水平层理, 定义平直且与层面垂直的层理为垂直层理。在层理角度和密度相同的情况下, 将垂直层理厚度设为0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0mm, 此组模型编号记为E₁— E₁₁。而另一组模型中水平层理厚度与垂直层理厚度设为相同, 编号记为F₁— F₁₁。两组模型中超声波传播方向均为由其设置的岩心柱塞底端向顶端传播。本文从中选取了层理厚度结构特征比较明显的E₂、E₁₀和F₂、F₁₀模型, 分析两组模型不同层理厚度对超声波特性的影响, 并通过其波场快照图进一步开展层理厚度对超声波特性影响的微观分析。其中, E₂的层理厚度为0.2 mm, E₁₀的层理厚度为1.8 mm, F₂层理厚度为0.2 mm, F₁₀层理厚度为1.8 mm。由波速数值计算结果可知, 随着层理厚度的增加, 无论水平层理还是垂直层理, 波速均呈线性递减, 水平层理的递减幅度更大（见图8a）。衰减系数计算结果显示, 随着层理厚度增加, 无论水平层理还是垂直层理, 衰减系数都呈增大趋势, 其中水平层理的衰减系数以幂函数规律递增, 而垂直层理以指数函数规律递增（见图8b）。超声波主频和主振幅计算结果显示, 随着层理厚度的变化, 水平层理中的超声波主频、垂直层理中的超声波主频与垂直层理中的超声波主振幅的变化无规律, 而水平层理中的超声波主振幅随层理厚度的增大整体呈下降趋势（见图8c）。

图8

	Figure Option ViewDownloadNew Window
	图8 不同层理厚度对超声波特性的影响（模型个数为22个）

前文预设模型中E₁₀垂直层理厚度比E₂的垂直层理厚度大, 超声波由岩心底端到顶端传播。波场快照图显示, 相同时刻下的超声波信号在E₁₀中传播所能到达的距离比在E₂中的更短（见图9a— 图9b）, 所以导致超声波波速在E₁₀中也较小（见图8a）。超声波在F₂与F₁₀中的水平层理传播也呈同样规律, 即随层理厚度增大, 在相同时刻下超声波信号在F₁₀中传播距离也更短, 说明超声波波速在F₁₀中也较小。图9中使用不同的声场幅值来表征超声波透射能量大小, 直观可见E₁₀的颜色总体较E₂更浅, 说明随着垂直层理厚度增大, E₁₀中透射超声波能量更小（见图9b）, 所以其衰减系数也较大（见图8b）。在模型F₂和F₁₀中的水平层理中也表现出与E₂、E₁₀同样的变化规律, 即随层理厚度增大, 透射超声波能量减小（见图9c— 图9d）, 衰减系数增大（见图8b）。

图9

	Figure Option ViewDownloadNew Window
	图9 层理厚度模型的波场快照

2.1.2层理密度对超声波特性的影响

在层理角度和厚度恒定的条件下, 通过调整层理数目来模拟层理密度变化, 即岩心剖面单位面积上层理的数目。将垂直层理厚度设为0.4mm, 层理数目设为0、1、2、4、5、6、10、12、15、18、20条, 对应的层理密度为0、0.08× 10⁴、0.16× 10⁴、0.32× 10⁴、0.40× 10⁴、0.48× 10⁴、0.80× 10⁴、0.96× 10⁴、1.20× 10⁴、1.44× 10⁴和1.60× 10⁴条/m², 此组模型编号记为G₁— G₁₁。另一组模型编号记为H₁— H₁₁, 为水平层理, 并将层理密度设置与上述一致。从中选取了能较好反映层理密度结构特征的G₃、G₇和H₃、H₇模型分析两组模型不同层理密度对超声波特性的影响, 并通过其波场快照图, 开展层理密度对超声波特性影响的微观分析。其中, G₃的层理密度为0.16× 10⁴条/m², G₇的层理密度为0.80× 10⁴条/m², H₃层理密度为0.16× 10⁴条/m², H₇层理密度为0.80× 10⁴条/m²。

由波速计算结果可知, 随着层理密度的增加, 无论水平层理还是垂直层理, 波速均呈线性递减, 垂直层理下降幅度相对较小（见图10a）。衰减系数数值计算结果显示, 衰减系数随着层理密度的增加呈线性递增（见图10b）。超声波主频和主振幅数值计算结果显示, 随着层理密度的增加, 垂直层理中的超声波主频和主振幅变化无规律; 但是在水平层理中, 随着层理密度的增加, 层理密度小于0.3× 10⁴条/m²时, 主频和主振幅呈下降趋势, 层理密度大于0.3× 10⁴条/m²后开始缓慢上升直到基本保持不变（见图10c）。

图10

	Figure Option ViewDownloadNew Window
	图10 不同层理密度对超声波特性的影响（模型个数为22个）

前文已对模型预设层理密度, 由于模型中G₇层理密度比G₃层理密度大, H₇层理密度比H₃层理密度大, 超声波由岩心底端到顶端传播, 在G₇中超声波通过层理与骨架的界面数目要比G₃多, 波场快照图显示, 相同时间下的超声波信号在G₇中传播的距离比在G₃中的更短（见图11）, 所以超声波速度也较小（见图10a）。超声波在H₇与H₃中的水平层理传播也呈同样规律, 即在相同时间下其超声波信号传播距离也更短, 说明超声波波速也较小。在衰减特性方面, 随着层理密度的增加, G₇中层理条数比G₃更多, 非均质性更强。由于层理与骨架之间的模量不等, 在骨架与层理之间产生一定的应变相位差, 从而在界面处产生内摩擦, 使机械能转为热能, 极大消耗了能量^[19]。图11中使用不同的声场幅值来表征超声波透射能量大小, 直观可见G₇颜色整体较G₃更浅（见图11）; 因此, G₇中的超声波透射能量要明显低于G₃, G₇对应的衰减系数更大（见图10b）。H₇与H₃类似, 也呈同样规律, 即随着层理密度的增大, H₇超声波透射能量更小, 说明衰减系数更大。

图11

	Figure Option ViewDownloadNew Window
	图11 层理密度模型的波场快照图

2.1.3层理角度对超声波特性的影响

在层理厚度和密度恒定的条件下, 将层理角度为0° 、10° 、20° 、30° 、40° 、50° 、60° 、70° 、80° 、90° 的模型编号记为I₁— I₁₀。文中选取了层理角度特征比较明显的I₃、I₈, 并通过其波场快照图进一步开展层理厚度对超声波特性影响的微观分析。其中, 模型I₃的层理角度为30° , I₈的层理角度为80° 。

数值计算结果显示, 随着层理角度的增加, 波速呈幂函数递减规律（见图12a）, 而衰减系数呈线性递增规律（见图12b）。

图12

	Figure Option ViewDownloadNew Window
	图12 不同层理角度对超声波特性的影响图（模型个数为11）

由于模型预设I₈的层理角度大于I₃的层理角度, 超声波由底端向顶端传播, 在层理密度相同的前提下, I₈的层理界面数目多于I₃, 波场快照图显示, 在相同时刻下I₈中超声波信号传播所能到达的距离比在I₃中的更短（见图13a— 图13b）, 对应的波速也就较小（见图12a）。在衰减特性方面, 由于I₈中的骨架与层理面之间的界面个数多于I₃, 导致阻碍超声波传播的路径也多于I₃, 图13中使用不同的声场幅值来表征超声波透射能量大小, 直观可见I₈的颜色总体较I₃更浅, 说明透射超声波能量更小, 所以衰减系数也较大, 对应的衰减系数更大（见图12b）。

图13

	Figure Option ViewDownloadNew Window
	图13 层理角度模型的波场快照图

在页岩数值模拟计算结果的基础上, 利用灰色系统关联法^[20]定量分析声学参数对层理结构的敏感性。首先, 将超声波速度、衰减系数等声学参数作为参考序列, 将层理角度、密度和厚度变化等层理结构变化参数分别作为比较序列：

${{Y}_{i}}=[{{y}_{i(1)}}, \ {{y}_{i(2)}}\ldots , \ {{y}_{i(n)}}]; \ {{X}_{i}}=[{{x}_{i(1)}}, \ {{x}_{i(2)}}\ldots , \ {{x}_{i(n)}}]$ （10）

对序列中的数据进行均值无量纲化处理, 即：

$\mathop{X}_{i(k)}=\frac{\mathop{x}_{i(k)}}{\sum\limits_{nk}^{In}{\mathop{x}_{i(k)}}}$, $\mathop{Y}_{i(k)}=\frac{\mathop{y}_{i(k)}}{\sum\limits_{nk}^{In}{\mathop{y}_{i(k)}}}$ （11）

求取序列的关联系数：

$\mathop{\xi }_{i(k)}=\frac{\underset{i}{\mathop{\min }}\, \underset{k}{\mathop{\min }}\, \left| \mathop{Y}_{i(k)}-\mathop{X}_{i(k)} \right|+0.5\underset{i}{\mathop{\max }}\, \underset{k}{\mathop{\max }}\, \left| \mathop{Y}_{i(k)}-\mathop{X}_{i(k)} \right|}{\left| \mathop{Y}_{i(k)}-\mathop{X}_{i(k)} \right|+0.5\underset{i}{\mathop{\max }}\, \underset{k}{\mathop{\max }}\, \left| \mathop{Y}_{i(k)}-\mathop{X}_{i(k)} \right|}$

（12）

获取整个数列的关联度：

$\mathop{r}_{i}=\frac{1}{n}\sum\limits_{k=1}^{n}{\mathop{\xi }_{i(k)}}$ （13）

声学参数中不仅速度对层理角度和层理密度的变化较为敏感（关联度值均为0.80以上）, 衰减系数对层理厚度的变化也较敏感（关联度值可达0.89）, 而主振幅、主频与层理结构相关联度则相对较小（关联度均在0.60以下）（见表1）。因此, 在进一步利用声学参数表征页岩层理结构的过程中, 除了选用常规敏感因子速度之外, 还需配合衰减系数来综合描述。

表1

表1

表1 层理结构参数和各个影响因子的关联度 层理结构 速度 关联度 衰减系数 关联度 主振幅 关联度 主频 关联度 层理厚度 0.76 0.89 0.52 0.58 层理密度 0.85 0.79 0.54 0.44 层理角度 0.81 0.67 0.58

表1 层理结构参数和各个影响因子的关联度

利用岩石波速建立动态力学参数模型是现场油气井开发中最常用的方法。其中各向同性模型和横向各向同性模型表达式如下所示：

各向同性岩石力学参数预测模型：

$\upsilon ={\left[ {1}/{2{{\left( {{{V}_{\text{p}}}}/{{{V}_{\text{s}}}}\; \right)}^{2}}-1}\; \right]}/{\left[ {{\left( {{{V}_{\text{p}}}}/{{{V}_{\text{s}}}}\; \right)}^{2}}-1 \right]}\; $ （14）

${{E}_{t}}={\left[ V_{\text{p}}^{2}\rho \left( 1+\upsilon \right)\left( 1-2\upsilon \right) \right]}/{\left( 1-\upsilon \right)}\; $ （15）

横向各向同性岩石力学参数预测模型^[21]：

${{\nu }_{\beta }}=\frac{{{\nu }_{2}}{{E}_{\beta }}}{{{E}_{2}}}-{{E}_{\beta }}\left( \frac{1}{{{E}_{1}}}+\frac{1}{{{E}_{2}}}+\frac{2{{\nu }_{2}}}{{{E}_{2}}}-\frac{1}{{{G}_{2}}} \right){{\sin }^{2}}\beta {{\cos }^{2}}\beta $（16）

$\frac{1}{{{E}_{\beta }}}=\left( \frac{1}{{{G}_{2}}}-\frac{2{{\nu }_{2}}}{{{E}_{2}}} \right){{\sin }^{2}}\beta {{\cos }^{2}}\beta +\frac{1}{{{E}_{1}}}{{\sin }^{4}}\beta +\frac{1}{{{E}_{2}}}{{\cos }^{4}}\beta $

（17）

通过数值模拟研究发现, 仅通过波速建立页岩的动态力学参数模型是不可行的, 还需配合对岩石各向异性特征敏感的衰减系数, 方能达到更好效果。利用实验测得的渝东南彭水地区井下龙马溪组页岩的波速和衰减系数, 基于麦夸特与全局优化算法^[22]在考虑量纲的基础上, 建立了该地区页岩的动态力学参数模型：

$\upsilon ={{p}_{1}}+{{p}_{\text{2}}}{{x}^{3}}+{{{p}_{3}}}/{\left( {{{V}_{\text{p}}}}/{{{V}_{\text{s}}}}\; \right)}\; $ （18）

${{E}_{t}}={{q}_{1}}+{{q}_{2}}{{V}_{\text{p}}}^{2}\rho +{{{q}_{3}}}/{\upsilon +}\; {{{q}_{\text{4}}}}/{{{\upsilon }^{2}}}\; $ （19）

采用渝东南ZY1井和YY1井的井下岩心开展基础物性实验获取岩心密度, 通过岩石力学实验获取弹性模量和泊松比, 利用超声波透射实验获取纵、横波速度（见表2）。在此基础上, 对本文模型进行验证。将表2中实验结果代入由公式（14）、（15）所建立的各向同性模型, 以及由公式（16）、（17）所建立的横向各向同性模型, 开展模型的对比分析。

表2

表2

表2 井下岩心实验结果 井名 层理 角度/(° ) 密度/ (g• cm-3) 垂向纵波速度/ (m• s-1) 垂向横波速度/ (m• s-1) 弹性 模量/ MPa 径向纵波速度/ (m• s-1) 径向横波速度/ (m• s-1) YY1 0 2.60 4 368 2 638 21 932 3 900 2 510 0 2.53 4 765 2 610 39 603 3 830 2 510 15 2.55 4 433 2 648 32 784 3 635 2 300 15 2.53 4 254 2 677 28 745 3 680 2 430 30 2.59 4 177 2 621 28 037 3 970 2 600 30 2.56 4 183 2 724 27 912 4 185 2 630 60 2.60 3 974 2 595 26 199 4 180 2 620 60 2.56 3 676 2 428 24 936 4 260 2 680 75 2.55 3 635 2 266 24 432 4 430 2 650 90 2.59 3 832 2 202 32 022 4 770 2 810 90 2.58 3 896 2 509 24 849 4 370 2 640 0 2.57 4 331 2 333 23 815 4 696 2 683 0 2.57 4 222 2 470 22 963 4 500 2 571 15 2.60 4 529 2 380 26 085 4 801 2 743 15 2.58 4 386 2 601 23 826 4 521 2 583 60 2.55 4 168 2 632 40 984 5 200 2 971 70 2.53 4 578 2 616 26 695 4 689 2 679 70 2.58 4 581 2 618 26 695 4 699 2 685 75 2.52 4 435 2 595 42 524 5 131 2 932 75 2.55 4 366 2 688 43 960 5 138 2 936 75 2.54 4 433 2 648 43 725 4 829 2 759 75 2.62 4 254 2 677 42 555 4 902 2 801 ZY1 80 2.60 4 634 2 729 26 198 4 788 2 736 80 2.53 4 410 2 748 26 022 4 773 2 727 80 2.55 4 646 2 563 26 796 4 756 2 718 80 2.53 4 711 2 773 27 783 4 431 2 532 90 2.52 4 598 2 658 44 516 4 707 2 690

表2 井下岩心实验结果

验证结果表明：利用本文模型计算得到的动、静态力学参数相关性整体较好, 其中在ZY1井中静态泊松比与动态泊松比相关系数可达0.83, 静态弹性模量与动态弹性模量相关系数可达0.84。在YY1井中静态泊松比与动态泊松比相关系数可达0.82, 静态弹性模量与动态弹性模量相关系数可达0.81（见图14）。本文的新模型模拟效果优于通常使用的各向同性模型和横向各向同性模型, 其中各向同性模型是指不同方向上物理性质相同的介质, 为均质模型; 横向各向同性模型是指垂直于轴线的任何方向上弹性相同, 为非均质模型。究其原因, 是因为各向同性模型只针对各向同性的岩石, 未考虑页岩的层理发育, 而横向各向同性模型将页岩等效为横向各向同性体处理, 模型结构过于理想化, 影响准确性, 而本文模型是基于表征页岩层理结构最为敏感的波速和衰减参数, 利用通用全局优化算法求得, 模型可较好地反映页岩层理结构的影响, 而且本文计算得到的岩石动态模量一般大于静态模量, 动态泊松比与静态泊松比基本为1︰1关系, 与Kuhlman等得到的关系基本一致^[23]。

图14

	Figure Option ViewDownloadNew Window
	图14 动态力学参数模型的对比验证图（R— 相关系数）

收集渝东南地区ZY2井的测井资料开展岩石力学动态参数剖面反演, 获取静态力学参数剖面结果。从图15中看出, 在反演深度范围内, 动态弹性模量数值普遍大于静态弹性模量数值, 动态泊松比与静态泊松比数值很接近, 基本为1︰1, 这与前文所述模型计算结果保持一致。同时, 取ZY2井龙马溪组不同深度的10块岩心开展岩石力学实验, 将实验结果与现场测井资料的反演数据进行验证。

图15

	Figure Option ViewDownloadNew Window
	图15 ZY2井下志留统龙马溪组岩石力学参数预测结果图

由表3可见, 10个不同深度反演得到的力学参数值与对应井下岩心的实验测试值基本一致。通过计算两者的相对误差可知, 泊松比的相对误差小于10%, 平均相对误差为6.94%; 弹性模量最大相对误差为21.25%, 平均误差为15.30%。可见, 利用测井资料反演的力学剖面能准确地反映地层力学参数真实情况, 本文建立的动态力学参数模型对层理发育的页岩地层有较好的适应性。

表3

表3

表3 现场静态力学参数验证误差分析表 深度/m 实验测量值 新模型反演计算 泊松比 弹性 模量/MPa 相对 误差值/% 泊松比 弹性 模量/MPa 相对 误差值/% 2539.74 0.26 26 356 8.86 0.28 30 542 15.88 2543.08 0.26 26 134 7.04 0.27 29 795 14.01 2547.20 0.25 25 135 9.52 0.27 29 588 17.72 2619.89 0.25 25 541 6.63 0.26 28 145 10.19 2627.20 0.24 23 839 8.48 0.26 26 797 12.40 2631.47 0.25 24 057 4.80 0.26 27 540 14.48 2660.43 0.25 23 036 6.58 0.26 27 933 21.25 2663.93 0.26 25 180 6.90 0.27 29 221 16.04 2668.81 0.25 25 289 7.17 0.27 28 609 13.12 2671.10 0.25 24 202 5.26 0.26 27 690 14.41

表3 现场静态力学参数验证误差分析表