地应力反演报告

3.2. 三维地应力场的多元线性回归原理
根据多元回归法原理，将地应力回归计算值k ∧
σ作为因变量，把有限元计算求得的自重应力场和构造应力场相应于实测点的应力计算值i
k σ作为自变量，则
回归方程的形式为
i
k
n
i i k L σσ∑=∧
=1
式中：k 为观测点的序号；k ∧
σ 为第k 观测点的回归计算值；i
L 为相应于自
变量的多元回归系数；k ∧
σ和 i
k σ为相应应力分量实测值和计算值的单列矩阵，n
为工况数。

假定有m 个观测点，则最小二乘法的残差平方和为
2
16
11
*)
(∑∑∑===-=m
k j n
i i jk
i jk
L S σσ
残 (2)
式中：
*jk
σ为k 观测点j 应力分量的观测值，
i jk
σ为i 工况下k 观测点j 应力分量的
有限元计算值。

根据最小二乘法原理，使得
残
S 为最小值的方程式为
∑∑∑∑∑∑∑∑∑∑∑∑∑∑∑∑∑∑===================
m
k j n jk jk m
k j jk
jk
m k j jk jk n
m
k j n jk m k j n jk
jk
m k j jk
m k j n jk jk m k j jk jk
m k j jk
L L L 16
1*16
1
2*
16
11*2
1
161
2
16
12161
2
2
16
1116
12116
1
2
1
)()
()
(σσσσσσ
σσσ
σ
σ
σσ
σ
σ
称
对
(3)
解此方程式，得n 个待定回归系数
T
n L L L L ),...,,(21= ，则计算域内任一点P
的回归初始应力，可由该点各工况有限元计算值迭加而得
∑==n
i i jp
i jp L 1
σσ
(4)
式中j=1，2，…,6对应六个初始应力分量。

最后，可以通过复相关系数r来衡量回归计算的效果，通过偏相关系数Vi将不显著因素剔除，以防止其余因素的退化。

3.3. 区域地应力反演计算
本次厂房区地应力场回归采用国际通用的大型三维有限元计算软件FLAC3D 与回归分析程序接口进行计算，最后基于应力叠加原理合成得到整个研究区域地应力场。

3.3.1. 三维计算模型
三维初始地应力场计算模型的坐标原点在大地坐标(X地-O-Y地)的位置为：X0 = 4617259.8013，Y0 = 39462944.75，Z0=1100m。

沿X轴和Y轴的计算范围分别为906m和952m，竖直方向从海拔254m到地表。

三个坐标的方位分别为：X轴S向，Y轴E向，Z轴与大地坐标的高程方向（垂直向上）重合，厂房轴线方向与X轴平行。

图3-1三维计算模型图3-2 地应力反演区域
地应力反演区域和三维数值模型分别如图3-1、3-2所示，红色的即为厂房的计算模型。

模型中主要考虑f347、f348、f350、f363、f369、f370、f371、f375八条大断层，各断层在局部坐标系的位置及产状如表1所示。

对计算区域共划分了578174个单元和99257个节点。

地下厂房区域围岩以Ⅲ类微风化岩石为主，岩体力学参数建议值见表3-2。

厂区岩性基本为花岗岩。

表3-1 影响地下洞室稳定的结构面描述
断层号走向倾向倾角宽度（m）
f347 18 NW 72 0.1~1.2
f348 30 NW 50 0.3~1
f350 298 SW 54 0.7~6
f363 10 SE 75 1.5~1.8
f369 50 SE 32 0.05~0.35
f370 300 NE 62 0.1~0.3
f371 60 SE 85 0.1~0.5
f375 305 SW 45 0.4~0.6
表3-2岩体力学参数
类别干容重(kN/m3) 弹性模量(GPa) 泊松比地层Ⅲ26.1 20 0.23
断层Ⅳ21 0.8 0.4
3.3.2. 实测地应力分析
在本次地应力回归计算域内，共有4个钻孔，为ZK305、ZK307、ZK309、ZK330(如图1所示)。

其中ZK309位于地表山体压力管道位置，为垂直孔；ZK305、ZK307、ZK309位于PD1探硐内，均为垂直孔，钻孔位置如表3-3。

丰宁水电站采用水压致裂法进行地应力测试。

在本次计算中，四个钻孔共有26个测点，选取其中规律性较好的16个测点进行优化计算，其他测点数值因离散性较大未考虑在拟合计算中，其中16个测点的地应力位置和数值如表3-4。

表3-3 钻孔位置
X(m) Y(m) 孔口高程(m) 钻孔
ZK305 4617033.8900 39463805.7700 1042.74
ZK307 4616913.8189 39463459.4229 1061.00
ZK309 4616821.8338 39463464.2161 1061.00
ZK330 4616730.8287 39463469.1173 1061.00
表3-4多元回归分析中采用的实测地应力位置及数值
测点深度测点高程σHσhσV 编号
ZK305（1）39.10~39.70 1021.9 12.98 7.18 8.10
ZK305（2）46.90~42.52 1014.1 13.46 8.06 8.31
ZK305（3）52.90~53.50 1008.1 14.12 8.32 8.46
ZK305（4）105.20~105.80 955.8 15.13 9.32 9.85
ZK307（1）34.00~34.60 1027 11.53 7.53 6.85
ZK307（2）48.50~49.10 1012.5 11.28 7.48 7.24
ZK307（3）52.90~53.50 1008.1 12.02 8.02 7.35
ZK307（4）93.20~93.80 967.8 13.01 8.71 8.42
ZK309（1）284.10~207.65 1136.3 7.53 5.34 7.53
ZK309（2）309.50~310.55 1110.9 8.8 6.1 8.20
ZK309（3）338.40~339.45 1082 13.57 8.88 8.97
ZK309（4）354.00~355.05 1066.4 13.74 9.54 9.38
ZK330（1）39.20~39.80 1022.2 13.18 7.03 6.78
ZK330（2）51.00~51.60 1010.4 13.43 6.52 7.09
ZK330（3）62.00~62.60 999.4 17.61 10.61 7.38
ZK330（4）67.70~68.30 993.7 16.86 9.46 7.53
3.3.3. 多元回归分析
(1)计算工况
首先将实测地应力数据的主应力与主应力方向转换到计算坐标下的6个应力分量值。

然后确定地应力场形成的拟合荷载构造模式为：①自重力、②在x方向的挤压构造力、③在y方向的挤压构造力、④在xy方向剪切构造力、⑤在yx向剪切构造。

其计算模型工况与边界如表3-5所示。

表3-5多种地应力因素计算工况及计算条件
(2)系数回归
通过对上述工况的计算表明,自重＋X向梯形构造应力+Y向梯形构造应力+XY 剪切构造应力+YX剪切构造应力，这一种工况相关系数最大，复相关系数r=0.9264121，说明在这种工况下回归较好。

回归系数1l=-1592291.75，2l= -1.192249，3l=-0.255401，4l=-0.272275, 5l=4.293182，6l=0.04408。

则测点的实测值与回归值的比较如表3-6、表3-7。

表3-6 测点的应力分量实测值与回归值对比表
编号类别x
σ
(mpa)
y
σ
(mpa)
z
σ
(mpa)
xy
τ
(mpa)
yz
τ
(mpa)
zx
τ
(mpa)
ZK305（1）ZK305（2）实测应力7.919 12.24081 8.1 -1.93416 0 0
回归应力7.9803 9.9873 7.6619 -3.6397 0.40203 -1.2521
ZK305（2）实测应力8.329 13.19069 8.31 -1.17552 0 0 回归应力8.1188 10.071 7.6767 -3.543 0.4221 -1.2775
ZK305（3）实测应力8.609 13.83061 8.46 -1.26286 0 0
回归应力8.2459 10.224 8.0323 -3.5498 0.2952 -1.2585
ZK305（4）实测应力9.226 15.13394 9.85 -0.3944 0 0
回归应力9.7474 11.92 10.281 -3.184 0.14675 -2.0668
ZK307（1）实测应力7.729 11.33112 6.85 -0.86951 0 0
回归应力7.4963 10.831 7.7164 -3.8351 0.20666 -1.4138
ZK307（2）实测应力7.669 11.09119 7.24 -0.82577 0 0
回归应力7.6137 11.196 8.0294 -3.8105 0.28268 -1.3946
ZK307（3）实测应力8.219 11.82131 7.35 -0.86912 0 0
回归应力7.7541 11.242 8.1603 -3.7722 0.24997 -1.4041
ZK307（4）实测应力8.924 12.79644 8.42 -0.93423 0 0 回归应力8.2682 12.24 9.209 -3.7359 0.28259 -1.387
ZK309（1）实测应力 5.593 7.27738 7.53 -0.69961 0 0
回归应力7.2427 8.6495 6.0432 -2.7768 -0.47793 -1.6954
ZK309（2）实测应力 6.412 8.4883 8.2 -0.86286 0 0
回归应力7.8362 10.255 6.7788 -2.6535 -0.41373 -1.6022
ZK309（3）实测应力9.422 13.02754 8.97 -1.49999 0 0
回归应力8.2373 10.161 7.5854 -2.2602 -0.6861 -1.8026
ZK309（4）实测应力10.025 13.255 9.38 -1.34219 0 0
回归应力8.2373 10.161 7.5854 -2.2602 -0.6861 -1.8026
ZK330（1）实测应力7.148 13.062 6.78 0.84367 0 0
回归应力 6.9217 11.941 7.8349 3.5305 0.26098 -1.4639
ZK330（2）实测应力 6.653 13.2973 7.09 0.94832 0 0
回归应力7.0843 12.22 8.1554 3.494 0.26473 -1.4808
ZK330（3）实测应力11.036 17.08446 7.38 -1.8447 0 0
回归应力7.4083 12.665 8.5161 -3.5194 0.27943 -1.4747
ZK330（4）实测应力10.016 16.3036 7.53 -1.9514 0 0
回归应力7.4636 12.747 8.7635 -3.5251 0.25262 -0.50794
表3-7测点的主应力实测值与回归值对比表
16个测点对六个应力分量实测值和回归值的对比图分别如图3-3～3-6所示。

总体上来说回归值和实测值还是基本吻合的。

图3-3 zk305测点的应力对比图
图3-4 zk307测点的应力对比图
图3-5 zk309测点的应力对比
图3-6 zk330测点的应力对比图
3.4. 地应力特征分析
受重力场、水平构造运动、断层等影响，从主应力实测值和回归值对比表格中可以看出最大主应力接近水平方向，最小主应力接近垂直方向，图3-7～3-12给出了厂房区域三个面的主应力分布等势线图，从图中可以看出：
(1)最大主应力等值线分布较均匀，最小主应力等值线在浅处相对较密。

从各面的主应力等值线分布情况可以看出：主应力值基本上都是从上到下逐渐增大,等值线在浅层受地形变化影响较大，说明反演回归应力场充分反映了地形、地貌的影响。

(2)在断层经过区域，受断层挤压剪切（数值计算中仅反映了挤压趋势）影响，该区域初始地应力表现为应力变化剧烈。

由于断层的存在及其自身变形模量较小，导致断层处岩体应力变化梯度较大，应力量值却较小。

(3) 由数值模型得到的厂房区最大主应力范围为13～18Mpa，而最小主应力范围为4.2~7.9Mpa。

与实测地应力相比，误差较小。

最大主应力方向与实测值相差约10°左右，对稳定性计算结果影响不大。

(4) 在稳定性计算中，应用侧压系数进行计算，选取厂房区中心点计算侧压系数，在此点最大主应力为15.7Mpa，中间主应力为9.6Mpa，最小主应力为6.1Mpa。

(5) 随着洞室埋深增大，主要洞室围岩的应力水平递增，从目前初始应力场的分布来看，三大洞室各处围岩的初始地应力水平以尾水事故闸门室的区域相对较高，主变洞次之，主厂房相对较低。

图3-7 主厂房轴线上最大主应力等值线图
图3-8 主厂房轴线最小主应力等值线图
图3-9 主变洞轴线最大主应力等值线图
图3-10 主变洞最小主应力等值线图
图3-11 尾水事故闸门室轴线最大主应力等值线
图3-12尾水事故闸门室轴线最小主应力等值线
3.5. 结论
本次厂房区地应力场回归采用国际通用的大型三维有限元计算软件FLAC3D 与回归分析程序接口进行计算，最后基于应力叠加原理合成得到整个研究区域地应力场。

获得如下认识：
(1) 运用逐步回归方法，“挑选”出对地应力场起显著影响作用的基本因素和经优选后的实测数据，使回归效果更加逼近实际应力场，同时尽量保证应力场解的唯一性。

通过对地应力的计算值和实测值对比以及对侧压比规律的综合分析，表明得到的初始应力场能够较好地反映地形、地貌和地质构造的影响，其结果具有较好的可靠性，可作为厂房设计、施工和稳定性评价的参考依据。

(2) 对丰宁蓄能电站初始应力场影响显著的因数有岩体自重、厂房轴向的水平挤压构造力以及垂直厂房轴向的水平挤压构造力，从计算结果可以看出σH钻孔截面内大次应力最大，水压致裂法和应力解除法的测量结果均表明最大水平主应力在y方向要占相对主导的地位，该工程区y方向的水平构造作用力要大于垂直向
重力作用。

但x方向则稍小于垂直方向的重力性地应力。

(3) 由数值模型得到的厂房区最大主应力范围在13～18Mpa之间，而最小主应力范围为4.2~7.9Mpa。

(4) 从较直观的角度来说，在包括厂房及主变洞的中心区域，在x方向的侧压系数Kx约0.9，而在y方向Ky则为1.3左右。