基于边界元方法的次级断裂信息挖掘试验研究_陆自清.pdf

第 48 卷 第 5 期 煤田地质与勘探 Vol. 48 No.5 2020 年 10 月 COAL GEOLOGY 2. Xi’an Research Institute Co. Ltd., China Coal Technology and Engineering Group Corp., Xi’an 710077, China Abstract In order to study the development of secondary faults and fracture, a boundary element BEM combined with seismic interpretation of main faults and geo-mechanics is discussed. Firstly, with the Anderson’s theory of faulting, the classification of secondary fault tectonic evolution events is improved under the combination of fault and stress; secondly, the stress tensor order reduction and linear superposition principle are applied to sim- plify the stress and displacement balance equation, and the Monte Carlo is used to inverse the paleo-stress field. After all, the equilibrium equation of stress, strain and discontinuous displacement is solved to obtain the current stress distribution of stratum and to analyze the state of secondary fault and fracture, fracture opening or fracture relative density. Based on the geological data of the experimental area, the boundary element geo-mechanical model is established to simulate the secondary faults. The distribution of the present stress field and the ination of the secondary faults are extracted. The experimental results show that the minimum hori- zontal principal stress is mainly in the direction of NW, and the cross fault zones are widely developed between the main faults, and the fault development zones have strong connectivity, which provides the reference data for the geological engineering design and hidden disaster assessment. Keywords BEM; secondary fault and fracture; stress field; Monte Carlo; fracture 地层中广泛发育的、与主断裂系统密切相关的 小尺度断裂，能显著提高流体气相、液相的导通 能力[1]，将这些与主断层密切相关的次生伴生断裂 定义为次级断裂[2]。 次级断裂的主要研究方法为地震信号处理与识 别、地质力学正演模拟两大类[3-5]。地震信号处理与 212 煤田地质与勘探 第 48 卷 识别方法，如全方位地震资料各向异性叠前反演[6]， 叠后地震信号拓频处理[7]，曲率体、相干体、方差 体等单属性、共偏移距向量片各向异性分析、多属 性复合分析计算[8-10]，以及混沌、断裂分形、蚂蚁 追踪、深度学习多层神经网络，卷积图形学，计算 机视觉模式识别等[11-14]，这些方法是对地震信号、 图像的处理与解释， 受地震资料采集与处理的影响 较大，易漏失真实的断裂响应，或者将处理流程中 不合理的噪声带入到结果中。地质力学正演模拟 研究包括有限元、有限差分、边界元等，如丁中 一等[15]提出的二元与破裂率法，联合裂缝发育程度 与岩石应力应变或者主应力、剪应力等参数建立数 学模型，利用库伦破裂准则定性、半定量预测裂缝 分布区；季宗镇[16]应用多参数预测法建立地应力与 裂缝参数关系式，预测裂缝的密度、开度信息等。 边界元法Boundary Element ， 简称 BEM 是在经典的积分方程基础上发展而来的，它与有限 元类似，但由于边界元法只在边界上划分，将数值 问题降维，具有数据输入少、计算速度快的优点。 断裂是地质演化时期古应力作用下不连续位移的结 果，几何形态上复杂多变，使用有限元、有限差分 方法，地质体网格剖分复杂，容易导致断裂描述简 单化，计算效率偏低的问题。边界元法的三角网格 剖分则能很好地表示弯曲的复杂不连续界面，方便 处理无界区域问题， 基本解满足无穷远处边界条件， 在无穷远处边界上的积分恒等于零，所以边界元法 特别适用于古应力条件下构造演化时的几何体破 裂与位移有关问题的研究。经典边界元地质力学 工具包 Poly3D通过自由多边形组合提高了弯曲断 面描述的精度，直接求解格林函数方程获取应力、 应变与位移，该技术在研究断裂方面取得了良好 的效果[17-23]。但受制于原始 Poly3D 程序本身的局 限性，上述案例中断层多边形数量上限较小，断层 多边形与应力作用模式设定单一，研究目标基本是 尺度较大的大地构造和相对定性的区域性宏观断裂 系统，缺少对较小规模的次级断裂研究。 次级断裂的发育程度和宽度受主断层性质、断 距大小及断层附近的岩石力学性质控制。从地质成 因出发，构造演化过程中产生的次级断裂，与主断 层构成形态及古应力场密切相关。岩石断裂的本质 是地层中聚集能量的耗散和应力场的重新分布，应 力应变逐渐形成次级断裂带，这是断层附近应力集 中分布的结果，断层形成后的应力场解除促使地层 内部应力场重新分布，重新达到平衡后形成现今应 力场。可以利用地震解释资料解释主断层，结合古 应力场反演，获取边界元模拟所需的多边形和应力 边界条件，使用边界元方法研究与地层主断层相关 的次级断裂系统。 1 次级断裂信息挖掘思路 边界元方法挖掘次级断裂信息正是基于断层发 育的地质成因，整合了岩体力学、地质力学、断裂力 学理论与方法，基本思路是利用地震资料解释成果， 建立边界元地质模型， 研究应力场约束条件， 预测地 下空间特定目的层裂缝分布范围和密度等属性。 如图 1 所示，通过地震解释主断层结果构建边界元几何模 型， 古应力场反演求解边界条件， 使用边界元模拟方 法，从中提取现今应力场分布与次级断裂属性。 图 1 次级断裂信息挖掘流程 Fig.1 The workflow of secondary fault and fracture ination mining 第 5 期 陆自清 基于边界元方法的次级断裂信息挖掘试验研究 213 2 边界元方法 岩体的断裂体系相当复杂，在绝大多数条件下 只能观察到有限的露头或者钻井岩心，岩体内部未 知的区域仍是绝大多数，在主断层周围存在规模不 等的次生伴生小断层、裂缝、裂纹等，是主断裂在 应力作用下的构造演化的结果[24-25]。 裂缝模拟是构造演化的数字化再现，分解为 3 个 主要步骤首先在地震解释主断层的基础上，依据 安德森断层分类系统和应力的组合模式进行分类， 配置构造演化独立事件；然后在露头、钻井等观测 信息的控制下，反演古应力场方位与相对大小；最 后综合古应力场边界条件与边界元多边形，使用边 界元法模拟构造演化过程，分析平衡状态下的地层 现今应力分布，分析衍生的次级断裂状况。 2.1 构造演化事件分解 构造演化事件中多边形与应力状态的组合多 变，直接求解应力与多边形构造形变会导致串行计 算过程冗长、低效。正逆断层与平滑断层应力状态 不同，以及应力与断层多边形的组合模式是构造演 化事件分类的基础。为了简化构造演化事件，需将 按照主断层几何形态与应力的组合关系进行分解。 安德森等学者分析了形成断层的应力状态[26]， 认为形成断层的三轴应力中的一个主应力趋于垂直 水平面， 并以此为依据提出了形成正断层图 2a、 平移断层图 2b、逆断层图 2c的 3 种应力状 态图 2， 矢量箭头方向与长短代表应力分量 σ1、 σ2、 σ3的方向与大小，其中 σ1σ2σ3，分别为最大主应 力、中间主应力、最小主应力。 图 2 安德森断层分类与应力状态 Fig.2 Stress states of Anderson’s theory of faulting 在安德森断层应力状态的基础上，根据断面与应 力相对的空间位置关系，不同断层应力组合特征下的 断裂演化事件，主要形成节理、缝合线、剪切断裂这 3 类次级断裂， 据此将主要构造演化事件分为 9 种图 3。 通过 9 类基本构造演化事件的区分，减少模拟计 算量，使得有限的模拟次数能够涵盖主要的应力与断 裂组合配置。 图 3 断层应力组合模式下 3 种裂缝形成结果 Fig.3 Three kinds of cracks under the combined stress mode of fault 2.2 古应力反演 基于应力张量降阶、线性叠加原理，使用蒙特 卡洛方法估计古应力场分布。 2.2.1 应力张量降阶 使用张量表述应力状态，如下式，应力边界条件 R σ由应力分量矩阵组成[27-28]。 R xxxyxz yxyyyz zxzyzz σσσ σσσ σσσ       σ 1 由于刚体连续介质的柯西应力张量是对称张量， 即满足 ijji ，式1即可简化为 R xxxyxz yyyz zz σσσ σσ σ      σ 2 安德森断层理论中[29]，设定唯一垂向主应力，x、 y 平面垂向分量为 0，式2可以简化为 R 0 0 xxxy yy zz σσ σ σ      σ 3 纵向静应力 zz σ不影响 R σ，式3可以转换为 R 0 0 0 xxzzxy yyzz σσσ σσ      - -σ 4 如式4中所示，将 R σ简化为 3 个未知分 量 xxzz σσ-、  yyzz σσ-、 xy σ的组合，因 R σ线性 独立，可以对 R σ进行分解 T R σR σR 5 214 煤田地质与勘探 第 48 卷 式5中， cossin sincos       R 6 1 2 0 0       σ 7 式6为二维可逆转矩阵，其中为最大水平主应 力方位角；式7为应力特征值矩阵，其中 1 ， 2 分 别对应最大、最小特征值。设定应力比φ，令 21 σφσ，将其代入式5 T R 10 0 θθ φ     1 σRR 8 当1固定为常量，则有  T R 10 , 0 θθ θ φ φ    σRR 9 据式9可知，根据已知点测量值 θ φ,即可表述 应力边界条件 R σ， 将应力张量的表达形式从三维转化 到二维。 2.2.2 线性叠加原理 构建独立断裂模拟方程的基础是基于线性叠加原 理的弹性位移方法，叠加原理避免了断层演化过程中 力与位移平衡方程的复杂性，应力与位移的线性组合、 分解是本文所用模拟方法的基本前提，使得边界元程 序能多线程并行模拟各分量的地质力学演化事件。 线性独立的常规表达式为 1 1223 3 fa fa fa f 10 式中 123 ,a a a,为实参数向量组； 123 ,ff f,为线性 独立的张量。根据格林函数，应变μ、应力 σ 、位 移 u之间的关系，存在如下平衡条件 123 123 123 123 123 123 aaa a σa σa σ a ua ua u        μ σ u 11 据式11，应变、应力、位移均可以表示为参 数向量组与 3 线性分量的组合，确定应力即需要估 计参数向量组，计算应力分量。 2.2.3 应力反演 在所有条件中，古应力是求解最为困难的未知 条件，给定主断层几何结构的多边形后，古应力场 求解即是求解应力、位移、形变之间的关系。存在 唯一的三元实数向量 α 组， 式11中应力线性叠加形 式简式为 σ A ασ 12 式12等式两边项各乘以 1 σ A-可得 1 σ  αA σ 13 式13表示给定 1 σ - A与应力 σ 后，则参数向量组 α 可解。由式9，使用蒙特卡洛方法搜索满足最小 成本函数的二维变量 θ φ, [30]， 再由式13获得实参 数向量组 α ，获得 3 个线性独立的应力分量 13 i σi≤ ≤，通过上述过程即完成了古应力反演。 2.3 边界元模拟 经典边界元 Poly3D 程序于 1993 年由 A. L. Thomas 使用 C 语言开发[18]，早期运行在 Sun 公司 的 Sun SPARC 处理器上，基于 Unix 运行环境，正 演模拟过程中，对所有的构造演化事件均采用线性 反演、元素迭代法直接求解，单线程串行计算过程 中，数值模拟效率不高，处理的多边形数量有限。 后来移植到 Linux 系统中，通过增加 MPI 补丁库的 方式提高并行性能，但由于正反演求解方法的局限 性，构造演化事件区分度不够，并行部分占比较低， 大部分运算过程依然是串行处理，全局效率提高幅 度较为有限。 通过地质演化事件分类，基于独立构造演化事 件进行线性分解与叠加，在此基础上对应力张量进 行 降 阶 ， 降 低 运 算 量 ， 同 时 结 合 QThread 、 QtConcurrent 等计算机线程控制 API 重新构建边界 元程序，使得边界元模拟方法适配现代计算机多核 多线程调度机制， 高效地模拟复杂的构造演化事件。 3 现场测试 3.1 区域概况 试验区块位于杨柳矿区童亭背斜的东北端，经历 多期构造运动，地质构造极为复杂。地层走向由西往 东为 NENWSN 向，变化较大，地层倾角较小， 一般为 510，断层发育，且断层相互切割剧烈。 H8 地层地震解释成果如图 4 所示，H8 地层为 图 4 试验区 H8 地层地震解释 Fig.4 Seismic interpretation in H8 ation of the experimental area 第 5 期 陆自清 基于边界元方法的次级断裂信息挖掘试验研究 215 一单斜构造，其走向为 NNE，倾向为 SEE，区域性 NE 向大断层将地层切割成 NE 向条状地堑、 地垒式 构造；受上述构造控制，区内断裂构造发育特征表 现为以 NE、NNE 向断裂构造发育为主；断层以 高角度正断层为主。F1F5 为地震解释主断层， W1W4 为钻井。 3.2 边界元试验 地质力学模拟过程中，古应力条件有一定的不 确定性，需要通过应力状态筛选，确定应力参数。 在井孔处提取不同应力状态模拟的结果，与井孔实 际统计数据进行比较，当各井孔统计数据与模拟结 果基本吻合时，以此时的应力条件作为断裂系统构 造演化的基准古应力边界条件。以 W1 井为例，如 图 5 所示，应力状态 1 的模拟结果与 W1 井点处断 裂统计数据基本吻合。 确定应力状态后，结合主断层多边形，完成 边界元模拟， 提取计算结果并形成 H8 地层次级断 裂相对密度平面图图 6，结果显示 H8 地层中广 泛发育与主断层交叉的次级断裂带，主断层之间 通过断裂发育带连通，其中 F1 与 F2 断层中部及 西北部连通性强，F3、F4、F5 之间的断裂发育连 通更好。 图 5 现今应力参数筛选 Fig.5 Screen of the stress parameters 图 6 试验区次级断裂平面分布 Fig.6 Plane distribution of secondary faults and fracture 3.3 结果分析 边界元数值模拟结果中获得的地层水平主应力 场分布及相对裂缝密度、开度等地层属性，可以应 用于水平井设计施工、煤矿防治水等相关工程。 图 7 为试验区 H8 地层最小水平主应力矢量图， 线段表示水平应力的大小与方位， 以近北西向为主。 煤层气储层的应力敏感性较强，分析应力场可以提 高开发煤层气的效率，煤层气压裂开发中，钻井方 位角的选择需要考虑最小水平主应力的方向，与最 小水平主应力方向垂直的近北东向为水平井钻井优 势压裂方位。 图 8 为裂缝平均开度预测。高含水砂岩到煤层 顶底板的流体运移通道不仅仅取决于主断层，次生 伴生断裂系统的流体导通能力也很强，同样是隐蔽 灾害的主控因素之一。由于裂缝发育导致地质结构 稳定性不强，裂缝发育区需要增强灾害防控与地下 水治理的力度。一般来说，裂缝开度与流体导通能 力成正比，以裂缝平均开度为研究对象模拟预测流 体优势运移通道，地层流体优势运移通道如图 8 所 示，暖色区域连通性较好，冷色反之。 图 7 试验区最小水平主应力矢量图 Fig.7 Vector diagram of minimum horizontal principal stress in test area 4 结 论 a. 边界元方法下的应力场数据模拟，在主断层 216 煤田地质与勘探 第 48 卷 图 8 试验区次级断裂平均开度 Fig.8 Average open degree of the secondary faults and fractures in the test area 解释的基础上，通过地质力学、边界元方法补充断 裂解释结果，能够实现地层应力状态、次生伴生断 裂的定性预测。 b. 应力场数值模拟方法一定程度上受制于观 测点的揭露状况，除孔中裂缝分析外，露头观测也 是一个重要的基础，将来如能应用激光断层摄影、 数字露头技术，多尺度、多分辨率的应力分析会更 好地支撑边界元数值模拟。 c. 目前准确定量表述次级断裂有一定的难度， 基础条件允许的情况下，如能通过岩石物理试验建 立本地区力学特征模型，将裂缝发育特征同地层弹 性参数和各向异性参数联系起来， 结合边界元方法， 定量描述不同级别断裂发育状况，对于工程实践会 有更积极的现实意义。 请听作者语音介绍创新技术成果 等信息，欢迎与作者进行交流 参考文献References OSID 码 [1] 李超峰. 黄陇煤田综放采煤顶板导水裂缝带高度发育特征[J]. 煤田地质与勘探，2019，472129–136. LI Chaofeng. Characteristics of height of water flowing fractured zone caused during fully-mechanized caving mining in Huan- glong coalfield[J]. Coal Geology Exploration，2019，472 129–136. [2] 刘源， KONIETZKY H. 基于离散元方法对走滑拉分盆地演化 及次级断裂扩展过程研究[J]. 地质力学学报，2019，255 840–852. LIU Yuan，KONIETZKY H. Particle-based modeling of crack propagation during pull-apart basin development[J]. Journal of Geomechanics，2019，255840–852. [3] 刘敬寿，丁文龙，肖子亢，等. 储层裂缝综合表征与预测研究 进展[J]. 地球物理学进展，2019，3462283–2300. LIU Jingshou，DING Wenlong，XIAO Zikang，et al. Advances in comprehensive characterization and prediction of reservoir fractures[J]. Progress in Geophysicsin Chinese，2019，346 2283–2300. [4] 刘志伟. 地震波方程时间域边界元法数值模拟[D]. 长春吉 林大学，2004. LIU Zhiwei. Value imitation of seismic wave equation using time-domain boundary element [D]. ChangchunJilin University，2004. [5] 管西竹，符力耘，陶毅，等. 复杂地表边界元–体积元波动方 程数值模拟[J]. 地球物理学报，2011，5492357–2367. GUAN Xizhu，FU Liyun，TAO Yi，et al. Boundary-volume integral equation numerical modeling for complex near sur- face[J]. Chinese Journal of Geophysics， 2011， 549 2357–2367. [6] 杨平，李海银，胡蕾，等. 提高裂缝预测精度的解释性处理技 术及其应用[J]. 石油物探，2015，546681–689. YANG Ping， LI Haiyin， HU Lei， et al. Interpretative processing techniques and their applications in improving fracture prediction accuracy[J]. Geophysical Prospecting for Petroleum，2015， 546681–689. [7] 王飞，程礼军，刘俊峰，等. 叠后地震属性识别页岩气储层裂 缝研究及应用[J]. 煤田地质与勘探，2015，435113–116. WANG Fei，CHENG Lijun，LIU Junfeng，et al. Research and application of post-stack seismic attributes in recognizing shale gas reservoir fracture[J]. Coal Geology Exploration，2015， 435113–116. [8] 狄贵东，孙赞东，庞雄奇，等. 应力场模拟约束下的碳酸盐岩 裂缝综合预测以塔中地区 ZG8 井区为例[J]. 石油物探， 2016，551150–156. DI Guidong，SUN Zandong，PANG Xiongqi，et al. Compre- hensive fracture prediction technology constrained by stress field simulation A case study from ZG8 area of central Tarim basin[J]. Geophysical Prospecting for Petroleum，2016，551150–156. [9] 田中英， 张茂省， 毋远召. 综合物探在高家湾滑坡群前缘裂缝 探查中的应用[J]. 工程地球物理学报，2019，166822–828. TIAN Zhongying，ZHANG Maosheng，WU Yuanzhao. Appli- cation of integrated geophysical exploration in the cracks detec- tion of the front edge in Gaojiawan[J]. Chinese Journal of Engi- neering Geophysics，2019，166822–828. [10] 张震，张新涛，徐春强，等. 渤海海域 428 潜山构造演化及其 对油气成藏的控制[J]. 东北石油大学学报， 2019， 434 69–77. ZHANG Zhen， ZHANG Xintao， XU Chunqiang， et al. Tectonic evolution and its controlling on hydrocarbon accumulation of 428 Buried Hill in Bohai Sea[J]. Journal of Northeast Petroleum University，2019，43469–77. [11] 樊晓东，李忠权，陈国飞，等. 井震联合断层识别技术在南贝 尔凹陷东次凹中的应用[J]. 成都理工大学学报自然科学版， 2018，455640–648. FAN Xiaodong， LI Zhongquan， CHEN Guofei， et al. Application of joint well-seismic technique to fault interpretation in Southern 第 5 期 陆自清 基于边界元方法的次级断裂信息挖掘试验研究 217 Beier Sag，Hailaer basin，China[J]. Journal of Chengdu Uni- versity of TechnologyScience Technology Edition，2018， 455640–648. [12] 庄益明， 宋利虎， 刘镜竹. 蚂蚁追踪技术在三维地震精细解释 中的应用以淮北祁南煤矿 82 采区为例[J]. 煤田地质与勘 探，2018，462173–176. ZHUANG Yiming，SONG Lihu，LIU Jingzhu. Application of ant tracking technology in 3D seismic fine interpretation of faultsA case study on mining district No.82 in Huaibei Qinan coal mine[J]. Coal Geology Exploration，2018，462 173–176. [13] 吴正阳，莫修文，柳建华，等. 裂缝性储层分级评价中的卷积 神经网络算法研究与应用[J]. 石油物探， 2018， 574 618–626. WU Zhengyang， MO Xiuwen， LIU Jianhua， et al. Convolutional neural network algorithm for classific action uation of frac- tured reservoirs[J]. Geophysical Prospecting for Petroleum， 2018，574618–626. [14] 周长松，邹胜章，夏日元，等. 太原盆地地裂缝发育规律及防 控措施[J]. 煤田地质与勘探，2016，44273–78. ZHOU Changsong，ZOU Shengzhang，XIA Riyuan，et al. Development regularity and control measures of ground fissures in Taiyuan basin[J]. Coal Geology Exploration，2016，442 73–78. [15] 丁中一， 钱祥麟， 霍红， 等. 构造裂缝定量预测的一种新方法 二元法[J]. 石油与天然气地质，1998，1913–5. DING Zhongyi，QIAN Xianglin，HUO Hong，et al. A new for quantitative prediction of tectonic fracturesTwo factor [J]. Oil Gas Geology，1998，1913–5. [16] 季宗镇. 卞闵杨地区阜宁组储层裂缝定量研究[D]. 青岛 中 国石油大学，2010. JI Zongzhen. Quantitative prediction of reservoir fracture of Funing Group in Bian-Min-Yang region[D]. QingdaoChina University of Petroleum，2010. [17] JEYAKUMARAN M，RUDNICKI J W，KEER L M. Modeling slip zones with triangular dislocation elements[J]. Bulletin of the Seismological Society of America，1992，8252153–2169. [18] THOMAS A L. Poly3DA three-dimensional，polygonal ele- ment，displacement discontinuity boundary element computer program with applications to fractures， faults， and cavities in the earth’s crust[D]. CaliforniaStanford University，1993. [19] POLLARD D， MAERTEN F， MAERTEN L， et al. Improved 3D modeling of complex fault geometries using Poly3D，an elastic boundary element code[C]//AGU Fall Meeting Abstracts，2001. [20] MEADE B J. Algorithms for the calculation of exact displace- ments， strains， and stresses for triangular dislocation elements in a uni elastic half space[J]. Computers Geosciences， 2007， 3381064–1075. [21] 龚发雄，单业华，林舸，等. 松辽盆地中部后裂谷早期正断层 系形成机制[J]. 地球科学中国地质大学学报，2008，334 547–554. GONG Faxiong，SHAN Yehua，LIN Ge，et al. Mechanism of early post-rift normal faults in the central Songliao basin， northeastern China[J]. Earth ScienceJournal of China Univer- sity of Geosciences，2008，334547–554. [22] WU K，OLSON J E. A simplified three-dimensional displace- ment discontinuity for multiple fracture simulations[J]. International Journal of Fracture，2015，1932191–204. [23] SWYER M， DAVATZES N， CLADOUHOS T， et al. Washington play fairway analysis-poly 3D Matlab fault modeling scripts with data to create permeability potential models[R]. DOE Geo- thermal Data Repository；AltaRock Energy Inc，2017. [24] GHOLAMI R， MOJOLAGBE J， MENSHOV A， et al. H-matrix arithmetic for fast direct and iterative of moment solution of surface-volume-surface EFIE for 3-D radiation problems[J]. Progress in Electromagnetics Research B，2018，82189–210. [25] PLATEAUX R，JOONNEKINDT J P，MAERTEN F M，et al. Maps of natural fracture reactivation likelihoodA comparison between homogeneous and heterogeneous stress fields[C]//Third EAGE Workshop on Naturally Fractured Reservoirs. European Association of Geoscientists Engineers，201811–7. [26] 朱志澄，曾佐勋，樊光明. 构造地质学[M]. 北京中国地质 大学出版，2008. ZHU Zhicheng，ZENG Zuoxun，FAN Guangming. Structural geology[M]. BeijingChina University of Geosciences Press， 2008. [27] MAERTEN L， GILLESPIE P， DANIEL J M. Three-