-1-
基于 Ansys、FLAC3D和 Fortran的某抽水蓄能电站初
始应力场反演回归分析
岑成汉 1，张燎军 1，左晓霞 2
1.河海大学水利水电工程学院，江苏南京（210098）
2.浙江省水利水电勘测设计院，浙江杭州（310002）
摘 要：本文结合某抽水蓄能电站工程地质条件以及实测地应力资料，分析了研究区域初始
地应力的影响因素。在各影响因素单独作用下，采用 ansys、FLAC3D以及 Fortran三种计算
程序结合，对建立的地质概化计算模型进行计算；在实测点地应力值与计算得到的应力值之
间建立多元回归模型，通过多元回归分析，求出最优回归系数。通过实测点的计算应力值与
现场实测值的比较，两者在量值上很接近，表明经过回归得到的地应力场是合理的，从而获
得某抽水蓄能电站初始地应力场较为合理的分布规律，为电站开挖模拟及长期稳定性分析提
供了合理的三维初始地应力场。
关键词：FLAC3D，地应力场模拟，多元回归，逐步回归分析，反演
1．引言
某抽水蓄能电站主副厂房洞室长，其中主机段长，副厂房长，安
装场长；厂房顶拱高程，厂房总高度；吊车梁以下开挖宽度为，
吊车梁以上宽度为。地下厂房区受区域构造影响大，构造较发育，岩层产状扭曲变化
较大，上库进口至上水平段为N75~85°E NW∠35~45°，下平段为背斜核部，下平段至下游岩
层产状为N70~80°E SE∠15~35°，从上库进出水口至斜井段为一等斜褶皱段，其中上库进出
口至上平段为上库向斜的南翼，同时也是背斜的北翼，背斜轴与输水系统轴线近直交，背斜
轴部位于厂房顶部。
根据地勘资料，输水线路沿线揭露的断层共62条，主要位于尾水反坡段核上平段，与洞
线近直交，走向以NE方向为主，在此我们主要考虑对电站工程区影响较大的三个断层：F11，
F521，F113，受断层核褶皱作用，输水系统部位挤压破碎带和节理近平行主干断层发育，且
分布疏密不均，在岩质较软弱的泥质砂岩夹层内多形成挤压破碎带，在断层之间多组成节理
密集带发育。
由于地质构造复杂，利用ansys的强大建模功能建立合理的计算模型，利用FLAC3D和
Fortran程序，采用多元回归分析的方法对某抽水蓄能电站的三维初始地应力场进行了分析。
本文结合ansys，FLAC3D，Fortran对该电站地下厂房区进行分析，得到合理的应力场分布，
使其能满足某抽水蓄能电站开挖模拟和长期稳定性分析的需要。
2．地应力实测结果
该电站的实测地应力值采用水压致裂法进行测量，在地应力场回归分析中，以计算坐标
系xyz 中坐标应力分量为基本对象，而地应力实测值所在的坐标系和我们模型的计算坐标系
不一致，需要将实测地应力进行正确的坐标变换以得到计算坐标系下的应力分量。由于计算
坐标系和实测值所用的坐标系中z轴是一致的，所以只需要转换两个坐标系的x，y轴，于是
我们所需要的两个坐标系的之间的方向余弦为：
cos
cos
⎩⎨
⎧
=
=
β
α
M
L
（1）
-2-
式中L为计算坐标系中σi与x轴正方向之间的夹角，M为计算坐标系中σi与y轴正方向之
间的夹角
根据每组实测主应力量值及方位角、倾角，可由式（2）将主应力转换成计算坐标系下
的应力分量：
式中
'
xσ ， 'yσ ， 'zσ ， 'xyτ ， 'yzτ ， 'zxτ 为计算坐标系下的各个应力分量， xσ ， yσ ， zσ ，
xyτ ， yzτ ， zxτ 为实测坐标系下的各个应力分量。
经计算，实测的地应力的6个分量结果以及转换后的坐标应力分量如表1所示。
⎪⎪
⎪⎪
⎩
⎪⎪
⎪⎪
⎨
⎧
=
=
++−+++=
=
++=
++=
zxzx
yzyz
xyxyxyxxyyxy
zz
xyxyy
xyxyx
LMMLLMML
LMLM
LMML
ττ
ττ
τσστστστ
σσ
τσσσ
τσσσ
'
'
222'
'
22'
22'
)2()()(
2
2
（2）
表 1 地应力实测值（MPa）
1 2 3 4 测点
项目 实测值 转换值 实测值 转换值 实测值 转换值 实测值 转换值
σx
σy
σz
τxy
τyz
τzx
3．回归分析理论与方法
回归理论
本文根据地应力实测资料、山体的地形状况及岩体的地质构造、岩体的力学性质等因素，
建立大范围山体的计算区域. 由于形成初始地应力场的主要因素有岩体自重和地质构造运
动等，对于自重因素，计算时采用岩体实测密度获得自重应力场σ自，对于构造应力场σ构获
取比较困难，需考虑多种因素，本文通过选取足够大的区域在边界上加法向、切向分布荷载
来实现.
对确定的计算模型模型，若取z向为垂直向上，x、y分别取水平方向，则此构造场σ构可
以由下列几种因素组合而成：由垂直于x 轴的面加沿x 方向的法向分布荷载引起的σ构1（如
图1所示）；由垂直于y 轴的面加沿y 方向的法向分布荷载引起的σ构2（如图2所示）；由垂直
于x 轴的面加沿y 方向的切向分布荷载及由垂直于y 轴的面加沿x 方向的切向分布荷载引
起的σ构3（如图3所示）；由垂直于y 轴的面加沿z 方向的切向分布荷载引起的σ构4（如图4
所示）；由垂直于x 轴的面加沿z 方向的切向分布荷载引起的σ构5（如图4所示）[2]。
由于实际所在不同的面沿不同方向分布载荷的具体量值是未知的，可通过对计算域边界
加相应单位分布载荷，利用常规有限元法对各种因素单独作用下进行计算，即先求得σ构1，
σ构2，σ构3，σ构4，σ构5，再按照现实的可行性，根据地应力实测点的应力值，建立初始
-3-
地应力场与各影响因素之间的回归方程，由测得的地应力资料， 用统计分析方法，使残差
平方和达到最小，可求得各影响因素的唯一解，进而得到可供工程使用的三维地应力场。
图 1 x向法向分布荷载 图 2 y向法向分布荷载 图 3 水平切向分布荷载 图 4 竖直切向分布荷载
通过上述分析确定了初始地应力场计为σ的六种基本影响因素，设
kbbbbbb εσσσσσσσ ++++++= 665544332211 构构构构构构 （3）
式中：bi（i=1，2，……，6）为回归系数；εk是随机变量，对于有N个观测点，三维问
题有6N个观测地应力分量，则有：①误差εk的数学期望值全为零，即E（εk）=0，k =1，2，……，
6N；②各次观测误差服从正态分布；③各次观测值相互独立，并有相同的精度，即εk之间的
协方差可表示为：
6N1,2,3,hk,
hk
hk 0
),( 2 ），（ ……=⎩⎨
⎧
=
≠= σεε hkCOV （4）
协方差反映了各实测点误差（εk，εh）间的关联程度，如k=h，即为方差，如k≠h，由于
各观测值是相互独立的，所以为零。
回归分析就是根据6N组地应力的实测值（σk）和数学模型中的“实测值”（σk1，σk2，σk3，
σk4，σk5，σk6）得到回归方程的回归系数b1、b2、b3、b4、b5、b6的估计值，根据得到的估计
值，回归方程可以写成：
kk bbbbbb εσσσσσσσ ++++++= 665544332211 构构构构构构 （5）
取 'kσ 为σk的估计值
665544332211
'
构构构构构构 σσσσσσσ bbbbbbk +++++= （6）
根据最小二乘法原理，这时观测值与回归值的偏差的平方和Q达到最小值，即
Q＝∑ ′− 2)( kk σσ （7）
达到最小。
求解出回归系数的估计值，我们可以有很多方法，常见的就是最小二乘法，在此我们使
用逐步回归分析的方法来求解[1]。
逐步回归的思想就是：将变量的逐个引入，引入变量的条件是其偏回归平方和经检验是
显著的，同时每引入一个新的变量后，对已选入的变量要逐个检验，将不显著变量剔除，保
证最后所得的变量子集中所有变量都是显著的，这样经若干步便得“最优”的变量子集。
测点的插值
使用ansys和Flac3D可以计算出自重应力场以及各个构造应力场的值，由此，我们得到了
各个测点的在自重应力场以及构造应力场的值，利用以上提出的方法，我们利用回归分析的
方法和理论编写Fortran计算程序，对该电站地下厂房的初始地应力场进行了回归反演分析。
但是由于四个测点的位置并非都在模型的节点上，所以测点的应力值应通过插值而得，
-4-
可根据该点所在单元的结点上的应力值以及此点在该单元中的局部坐标值进行插值，插值方
法如下：找到该点所在的单元及在这个单元中的局部坐标，有了局部坐标值，再利用地应力
值通过形函数来插值求得该点的地应力值，公式如下：
∑
=
=
n
i
iiN
1
σσ
式中n为空间单元的结点个数，σ为空间单元结点上某一已知应力分量，Ni为应力插值形
函数，该插值过程可由Fortran程序完成。
当然，还可以在建模的时候就把实测点作为单元的结点，这在有限元软件ansys 中是可
以实现的，只需要把实测点设置成“硬点”，则在剖分单元的时候，这些“硬点”必定是单元的
结点，这是一种比较优越的方法。但是，对于一些重要的工程，随着工程的进行，往往会增
加实测点的数目，新增加的实测点的空间坐标是事先无法预知的。
4．厂区地应力场反演分析
模型建立
通过分析某抽水蓄能电站工程区的范围和该处工程地质与水文地质条件以及地应力实
测点的分布资料，确定计算区为：x，y轴的计算范围为 m× m；z轴竖直向上，从
-170m高程一直到自然边坡面。沿电站主厂房轴线方向取为x轴，y轴向垂直于主厂房轴线方
向指向上游。
利用ansys强大的建模功能，建立所需要计算的模型，同时考虑了对电站工程区影响较
大的三个断层：F11，F521，F113，因断层均有一定的厚度，所以采用实体单元进行模拟，
图5为在FLAC3D中的计算模型。
回归分析计算结果对比
通过对四个测点地应力值进行的回归分析计算，求得式子（5）中的各个回归系数分别
为：b1＝，b2＝，b3＝，b4＝，b5＝，b6＝，εk＝
，则该电站地下厂房的初始应力场的回归结果为：
654321 −−++++= kkkkkk σσσσσσσ
式中σ为回归的初始地应力场的值，σk1为自重应力场中的应力值，σk2为由垂直于x轴的
面加沿x方向的法向分布荷载产生的应力场中的应力值，σk3为由垂直于y轴的面加沿y方向的
法向分布荷载产生的应力场中的应力值，σk4为由垂直于x轴的面加沿y方向的切向分布荷载
及由垂直于y轴的面加沿x方向的切向分布荷载产生的应力场中的应力值，σk5为由垂直于y轴
的面加沿z方向的切向分布荷载产生的应力场中的应力值，σk6为由垂直于x轴的面加沿z方向
的切向分布荷载产生的应力场中的应力值。
图 5 计算模型 图 6 x＝0截面的 σz 等值线云图
-5-
图 7 x＝0截面的最小主应力等值线云图 图 8 x＝0截面的最大主应力等值线云图
回归分析中通过计算，可以得到复相关系数为，表明回归公式相关性较好，结果
表明回归应力值合理，且较接近于实测点的地应力值，测点地应力值和回归地应力及其比较
见表2，由表2可以看出回归值与实测值之间接近，计算值和实测值拟合的较好，图6为在
FLAC3D的模型的x＝0截面的σz等值线云图，图7、图8分别x＝0平面的最小和最大主应力等
值线云图，可以看出，应力场较符合规律性，因此，反演得到的初始地应力场可用于工作设
计中。
表 2 测点处地应力值及回归结果 （MPa）
构造场
测
点
应
力
分
量
自重场
σ 自 σ 构 1 σ 构 2 σ 构 3 σ 构 4 σ 构 5
回归值 实测值 误差
σx
σy
σz
τxy
τyz
Ι
τzx
σx
σy
σz
τxy
τyz
ΙΙ
τzx
σx
σy
σz
τxy
τyz
ΙΙΙ
τzx
σx
σy
σz
τxy
τyz
ΙV
τzx
5．结论
（1）本文简述了结合Flac3D、ansys、Fortran，根据地质勘测和实测地应力资料，反演求得
某抽水蓄能电站岩体初始地应力场，表3的比较分析结果表明，结合Flac3D、ansys、Fortran，
-6-
采用多元回归与三维地质模型的有限差分计算相结合的初始地应力反演回归方法能够得到
合理的应力场分布，反演回归应力值与实测应力值拟合较好，为电站开挖模拟及长期稳定性
分析提供了合理的三维初始地应力场。
（2）岩体的初始应力场主要是由重力以及地质构造运动所产生，由重力产生的自重应力场
比较容易确定，而由构造运动产生的应力场很难确定，从本文的地应力场回归效果来看，采
用三维有限元反演回归初始应力场较好地反映了山体地形、地貌和地质条件对初始应力场的
影响，能较好地与实测应力值吻合，是一种较为有效的反演回归方法。
参考文献
[1] 周纪芗. 实用回归分析方法[M] . 上海：上海科学技术出版社，1990.
[2] 邵国建. 岩体初始地应力场的反演回归分析. 水利水电科技进展. 2000（5）：36～38.
[3] 胡 斌，冯夏庭，黄小华，苏国韶，周 辉. 龙滩水电站左岸高边坡区初始地应力场反演回归分析. 岩石
力学与工程学报. 2005（22）：4055～4063.
[4] 博弈创作室.ANSYS 基础教程与实例详解.北京：中国水利水电出版社,2004.
[5] 刘波，韩彦辉. FLAC原理、实例与应用指南. 北京：人民交通出版社，2005.
REGRESSION ANALYSIS OF INITIAL GEOSTRESS
FIELD BASED ON ANSYS, FLAC3D, FORTRAN FOR A
PUMPED STORAGE POWER STATION
Cen Chenghan1, Zhang Liaojun1, Zuo Xiaoxia2
of Water Conservancy and Hydropower Engineering, Hohai Univ. Nanjing
Jiangsu（210098）
Design Institute of Water Conservancy & Hydro-electric Power,Hangzhou
Zhejiang (310002）
Abstract
Based on the measured data of geostress and engineering-geological conditions in the region of the a
pumped storage power station，possible factors affecting the initial geostress field are discussed in this
paper. Under individual effect of factors，the established simplified geological model is calculated by
ansys，FLAC3D and Fortran. Then，a multivariate regression model is established between the actually
measured and corresponding calculated results of geostresses. Through multiple regression analysis，the
optimum regression coefficient can be derived. Through the comparison between computed and
measured geostress values of measuring points，it has been found that the two are similar in values，
which suggests the accuracy and reasonability of the geostress field regression. Thereby， the
distribution tendency of initial geostress field of the a pumped storage power station region can be
obtained，providing a rational 3D initial geostress field for excavating simulation and analysis of
long-term stability of the Hongping pumped storage power station.
Keywords：FLAC3D；simulation of in-situ stresses；multi-factor regression；step-by-step regression
analysis；back analysis
作者简介：岑成汉，男，1980年生，广西来宾人，工学硕士，主要从事地下洞室稳定性的研
究。
