地应力反演报告

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所示。