图5 剪切力分布图
图6 二元变形网格图
图7 铜锣峡构造长兴组岩石破坏接近程度等值线图
岩地层的破坏接近程度值的相对大小，可将研 究区划分为不同的岩石破坏区或裂缝发育区。 1）在断层和断层附近地区，其破坏接近程度值很 高，属地层破碎区。 2）在铜锣峡主构造区，其破坏接近程度值大于 1.248，岩石破坏程度比较高，属破裂发育区。 3）在铜锣峡构造西北边平缓带及东南边部分高陡 地带, 破坏接近程度值在1.142～1.248之间，为破 裂临界发育区。 4）在东南边部分高陡地带, 破坏接近程度值在 0.99～1.142之间，为欠发育区。 5）小于0.99为破裂不发育区。
2
F
W W x y
2 2
W G 1 y

2
W W 1 y x
2
W xy
2
W W 1 y x W W 1 y x
g ( x, h ) 1 Z ( x) Z ( x h)2 2
根据样品点计算的变差函数叫做实验变差函教， 其计算公式为：
1 N (h) Z ( xi ) Z ( xi h)2 r (h) 2 N (h) i 1
i
式中：xi—为第个观测点的坐标； Z(i)、Z(xi+h)—分别为xi及xi+h两点处的 观测值； h—为两观测点间的距离; N(h)—为相距h数据对数目； r(h)—为实验变差函数的值。
前处理主要有作图、识别分析域、约束条件设定、 荷载条件的设定、网格生成、材料参数的设定等功 能，作图窗口用于定义分析对象的形状、材料范围、 挖掘断面等几何信息。识别分析域的功能用于区分 不同材料（地层），并指定各区域材料。约束条件 设定是设定所选择边界约束，载荷条件的设定用于 设定所选择边界或点的载荷。网格生成根据所设定 的网格分割数自动生成有限元计算网格。材料参数 设定通过填写参数设定对话框来完成材料参数的设 定。通过前处理，得到有限元分析所需要的所有数 据后，便可进行分析计算。
应力场那样用实测点来进行拟合模拟，因为古构造 应力场，特别是比较久远的构造应力场，现在还没 有办法对个别点的数值大小和方向进行实测。目前 探索使用的Kaiser效应方法是一个新的方法，对于 近期的地应力的测量比较准确，而对于年代久远的 地应力的测量结果尚在探索中。所以古构造应力场 只能根据现今的构造行迹来进行应变场的模拟，这 种模拟需要有正确的构造发展及演化的认识结果， 并且模拟结果仅是相对值，不能当实际应力值使用。 3、有限元软件简介 二维计算软件主要分为前处理、计算和后处理三 大功能块，现将其功能简单介绍如下：
表1 岩体物理力学参数表
材料名 灰岩 断层 大断层 高点 弹性模量 45000 34000 32000 34000 泊松比 0.26 0.23 0.22 0.22 残余内聚 力 4 2.4 2 2.6 岩体抗拉强 度 45 29 28 39
断层高陡带
陡坡 缓坡 构造发育带
35000
42000 40000 38000
① 趋势面拟合法
趋势面拟合可以构造曲面整体拟合，也可以局部内插加密拟 合，此次研究采用后一种方法。通常局部内插三次趋势拟合度一 般达到90%以上。其趋势面方程为：
f Ax3 By3 Cx 2 y Dxy2 Exy Fx2 Gy 2 Hx Iy J
将构造面抽象为一个数学曲面 W W ( x, y)
图3 最大主应力值分布图
图4 最小主应力分布图
图5为最小主应力分布图，其最小主应力方位基 本上沿北东—南西向展布，一般在16—26 Mpa之间 变化。在断层端部形成应力集中区，而在铜锣峡西 侧平缓带应力值较小。 剪应力分布特征十分明显，如图6所示，在最大 主应力值比较高而最小主应力值比较小的地方会出 现剪应力集中的情况，由最大主应力值分布图和最 小主应力分布图可看出，在铜锣峡构造带的南高点 及铜1、6、7井地区，最大、最小应力值相差较大， 因而在剪应力分布图上表现为剪应力集中带。另外， 在断层处往往出现剪应力集中，这有利于岩石产生 破裂。
K LN M 2 EG F 2
它在点（X，Y）处的三维最大主曲率为：
1 R1, 2 H H 2 K
LN M 2 其中：K EG F 2
EN 2 FM GL H 2( EG F 2 )
W E 1 x
2W L x 2 M
区域现今应力场的总体规律，可以在调查断层的新 构造活动特征、震源机制解和地应力实测的基础上 得出初步认识。但要定量的反映区域应力场，找出 应力集中部位，则需要通过地应力实测和数值模拟 来实现。所谓区域应力场的有限元反演，就是采用 有限单元法根据已有的已知地应力实测点和震源机 制解来推求整个区域的地应力场。其方法是首先根 据区域地质调查结果，建立研究区的地质力学模型； 然后通过不断改变边界力作用方式和大小量值（包 括大小和方向）与已有地应力实测结果和地震震源 机制（最大主应力大小和方向）达到最佳拟合。由 此即可得出反映研究区现今应力—形变场的真实情 况。然而，古构造应力场不能像现代区域
（2）计算模型范围 计算模型是根据新解释的《四川盆地铜锣峡构 造下三叠统飞仙关组底界地震反射构造图》为基础 建立的，模型的边界北以84TLX—D323，南以 86TLX—D285测线为界，东以主测线85TLX— D668及86TLX—D619为界，西以砂坪场两岔一线 为界组成一个矩形模型。 （3）地质结构模型 根据《四川盆地铜锣峡构造下三叠统飞仙关组 底界地震反射构造图》及其它地质资料，抽象出用 于实际计算的地质结构模型，如图所示。
0.235
0.255 0.245 0.24
2.8
3.5 3.4 3
31
40 36 33
（4）模型的离散化 在整个模型区域内采用8节点四边形单元和六节 点三边形单元将整个计算模型离散为9432个节点， 3081个单元，经离散后的计算模型如图2所示。
图2 计算模型图
（5）计算结果及分析 在上述地质模型的基础上，利用日本软脑公司的二 维有限元计算软件2D—δ进行数值分析（计算所得的应 力场遵循弹性力学约定，即张为正、压为负），可得到 模拟计算区的最大、最小主应力值等色谱图，最大剪应 力等值色谱图。 从图3可知，在晚二叠纪时期，研究区持续受到来 自东南—西北向的挤压力作用。图4为最大主应力值分 布图，表明研究区内最大主应力主要分布在西北—东南 方向，在各构造高点上出现最大主应力值，最大值为 115.022Mpa。在断层内部，由于应力释放，形成低应 力区，但在断层周围的岩体边缘却有应力集中，因此， 往往在断层附近，主应力值变化较快。
式中：
H 12 (2Y0 H 3 ) X 2 X0 2X
H 12 H 32 2 X 2 Y0 2(H 1 H 3 )
H 3 H 2 H 3 H 1 H 1 H 2
裂缝孔隙度及渗透率计算
裂缝孔隙度有多种计算方法，其中Murray的计算 方法应用较多。该方法是以岩层受拉张力产生弯曲裂 开的一个单元为基础 ，利用几何形态来计算该单元的 裂缝孔隙度。其裂缝孔隙度为：
差分法的曲率计算与趋势面中的曲率计算相同。
③ 样条函数法
这是一种曲线拟合的方法，美、英称为Spline函 数。通常以三次样条函数来模拟复杂的曲线轨迹，其 计算公式为：
u ui 1 u ui 1 S (u ) xi 1F0 h xi F1 h i i mi 1G0 u ui 1 m1Gi u ui 1 hi h h i i
曲率半径计算公式：
R 1 d 2 z / dx 2
④ 三点圆弧法
数据网格化后，在每一行上依次提取等距三点的 海拨高度，在构造剖面上便构成一个圆弧。根据点距 及深度即可计算圆心坐标，进而求出三点的平均曲率R。
ห้องสมุดไป่ตู้0 X
H2
L2
2X
H3
H1
L1
L
利用圆心坐标计算曲率半径的公式为：
R 1 X 02 Y02
用多种不同的材料表示出了该区的岩性（生物灰 岩），以及主要构造、断层、高点、凹陷等特征， 涉及到了铜锣峡、统景场、明月场、天口场、两岔 等构造，还有东北等方向的断层，较准确的反映了 实际地层的情况。 物理力学参数计算采用弹性结构模型，用基于 Mohr—Coulomb原理破坏接近度来判断岩石材料 的破坏状态，计算所涉及的物理力学参数主要包括 弹性模量、泊松比、残余内聚力及岩体抗拉强度， 其确定原则主要依据实验测试结果和工程地质类比 方法确定，各参数选取结果见表1。
后处理主要用于查看及分析结果。后处理为用 户提供了丰富的功能，主要有位移和应力的色谱图、 等值线图、断面曲线、轴力曲线、各分量的数值表。 4、铜锣峡构造应力场二维有限元计算 （1）计算模型的建立 据有关研究，川东铜锣峡地区长兴组地层主要 受海西期构造应力的作用，以东南—西北的构造挤 压力为主，从而形成了现今北东—南西的构造面貌， 为了模拟古地理构造应力场，在二维有限元构造应 力作用的模拟中，采取了地质结构模型和离散化网 格加不同的载荷约束边界来实现。
2
2
2W N y 2
2
式中： R1,2—曲率半径； K—总曲率； H—平均曲率。
② 曲面差分法
地质统计学认为，在空间观测到的一个地质参数 往往不能用一个简单的确定性的函数来表示，它们是 随机的。但这些值又显示出一定的空间连续性，相互 之间又有一定的相关性，对这既具有空间上的相关性 又具有随机性的变量称为区域化变量。 变差函数是指区域化变量Z在与两点处的增量的方 差之半，即区域化变量在相距为h的任意两点处的平方 差值的一半。也称理论变差函数，其公式为：
（6）裂缝预测
表2 铜锣峡地区裂缝预测的K值标准