锦屏复杂结构谷坡应力场反演模拟与特征分析

徐佩华，黄润秋，陈剑平，袁中凡

（1.吉林大学 建设工程学院，长春 130016；2.成都理工大学 地质灾害防治与地质环境保护国家重点实验室，成都 610059）

1 前 言

地应力是水利水电工程的基本荷载之一，对工程的布置和设计，尤其是对深挖工程岩体的变形和稳定具有重要影响[1]，因此，对河谷区岩体地应力场的研究是水利水电工程建设中普遍关注的一个重要研究课题。现场测试是目前工程中直接获得岩体地应力资料的惟一方法，但其存在诸如受经费限制、实测点数有限、受局部地质构造、地形、测试方法等的影响，测得数据很分散等缺点，很难从实测数据中获得整体地应力场特征。一些研究者提出了利用少量测量资料拟合出一定范围内的全场地应力分布的数值模拟方法，如地应力回归分析法[2]、地应力趋势分析法[3]、正交设计法[4]、应力试算法[5]。地应力回归分析法和趋势分析法的可靠性主要依赖于测点数量和精度，对地质条件的反映不足[4]；正交设计法是对描述河谷形成前小区域地应力场的状态参数K1、K2、T1、T2和地表形态参数D、D0，进行科学安排与分析多因素试验的一种方法，前提是实测点数据不受或少受局部构造、现今河谷形态等因素影响，能实际反映小区域地应力；应力试算法是在边界施加一定的水平应力或位移或速度，进行试算，当试算的结果与实测数据有较好吻合度时，认为此时的应力场分布即为实际斜坡的应力场。

河谷区地应力场是在区域地应力的基础上由于地表剥蚀、河流侵蚀等地质作用或受河谷地形影响而形成的局部地应力场[6]，若只是据边坡目前几何形状建立计算模型进行应力场数值模拟，则忽略了河床迅速下切过程中，侧向及垂向卸荷作用对边坡一定深度范围内应力场的影响[5]。这样反演出的谷坡应力场显然与实际会有较大的差别，尤其是当要研究与河谷下切过程相关的地质问题时，河谷下切过程的模拟必不可少。因此，人们在进行锦屏深裂缝及其有关变形与稳定性问题研究时，都对枢纽区河谷进行了河谷下切模拟[7－8]。

锦屏一级水电站为高地应力区，在水电可研、勘察及建设阶段都遇到了诸多由高地应力形成的独特现象，为此，人们运用ANSYS[9－10]、3D-σ[11]、FLAC3D[7]、UDEC[8]、3DEC[8]及其他有限元法[12]等多种手段进行了该区地应力场的反演计算。然而，以往的反演计算，没有进行河谷下切模拟，或者对谷坡结构的概化过于简单，对锦屏河谷应力场分布特征也没有进行系统地分析。因此，本文根据现场调查及前人研究成果，在分析坝址区工程地质条件的基础上，概化复杂结构谷坡，采用有限差分法模拟河谷演化发展，对照实测地应力数据，反演地应力场，并对其特征进行分析。

采用应力试算法进行反演，其结果取决于以下几个方面：（1）边界施加应力的量值和方式，包括采用何种边界，如应力边界、速度边界、位移边界；何种加载方式，如先计算重力场，然后施加均匀的水平向应力边界或是速度、位移边界，或是直接施加梯度边界。对称加载，还是一边加载，一边固定；一个方向加载还是两个甚至三个方向加载。（2）谷坡结构的概化，包括岩性、岩组结构、地质构造条件、风化界限等。（3）河谷发育条件概化，包括河谷类型、形态、发育时间、阶地高程、岸坡线演化等。本文对以上几个方面将给出详细介绍，并与前人的模拟过程进行对比，并根据计算结果分析锦屏现今河谷应力场特征。

2 工程概况与基本地质条件

锦屏一级水电站位于四川省西南部雅砻江上游大河湾水文站至手爬沟河段之间的下游段，地处四川省木里县和盐源县交界部位。该水电站设计坝高305 m，正常蓄水位为1 880 m，装机容量为3 200 MW，库容为77 亿 m3，是雅砻江流域规划拟建的最大库容水电站[11]。大地构造位置属扬子板块西缘，青藏高原东部边缘的川滇棱形断块中段东部的锦屏山断裂东侧。坝区主要结构为三滩向斜、普斯罗沟背斜及 4 组断层，左岸发育的断层主要有f5、f8、f2。发育的结构面可分为 4 组：层面N20º～35ºE/NW∠ 30º～45º，控制边坡稳定性的 2 组节理 ， 即 N50º～ 70ºE/SE∠50º～ 80º和 N20º～40ºE/SE∠60º～85º，以及 1 组不怎么发育的节理：N50º～70ºW/NE[13]。

坝区河流流向N25ºE，河道顺直狭窄，两岸山体雄厚，岸坡陡峻，为典型的高山峡谷地貌、深切V型河谷。左岸地层总体产状N15º～30ºE/NW∠30º～55º，反向坡，下陡上缓。由坝址从上游往下游1 850～1 900 m高程线以下为三叠系杂谷脑组段大理岩出露坡段，坡度为55º～70º。1 900 m高程以上为变质砂岩、粉沙质板岩出露段岸坡，坡度为40º～50º[8]。

3 模型概化

3.1 谷坡模型概化

本文选取与坝轴线相交的 II1勘探剖面为研究对象，图1为该勘探线的地质剖面图。由图可见，该剖面中的地质结构面除了前面介绍的 f5、f8、f2断层外，还有 f42-9，以及左右岸各有 1条黄斑岩脉，右岸有f13、f14断层。本文概化的模型涵盖了以上所有断层和黄斑岩脉，通过岩体质量的调整，概化了强卸荷线、中卸荷线及弱卸荷线，模型涵盖了左右岸。各岩体按图1中各岩性的岩体质量进行设置，岩体质量参数见表1，参数取值与文献[7]几乎一致。如图 2所示为概化后的 II1岸坡的模型图，X轴为剖面方向，即S62°E方向，由左岸指向右岸，宽2 833 m；Y轴为高程方向，底板高程为750 m，坡面最高高程为3 022 m，相对高差2 272 m。Z轴为河流方向，本次模型平行拉伸Z方向，构成假三维模型，进行三维应力计算与分析。

相比较而言，文献[10]中的模型采用对称边界简化了右岸，且左岸只有 f5断层及中卸荷线。文献[7]中的模型缺少f42-9、f14以及2条黄斑岩脉，只有1条卸荷线；文献[8]中的模型左岸只有f9、f2断层，右岸则有f13、f14断层和2条黄斑岩脉，模型中考虑了层面；文献[12]未见有关断层的描述，文献[7，10，12]模型的X方向与本文相同，而文献[8]模型中的X方向为正东，Z轴方向为正北，y轴正向向上。

图1 III勘探线地质剖面图Fig.1 The geological section III profile

3.2 河谷发育条件概化

图2 III勘探剖面的模型概化图Fig.2 The generalization numerical model of geological section III

分析已有的地质资料可知，该河段共有六级阶地，离锦屏最近的阶地高程为洼里阶地，因此参照洼里阶地的拔河高度，对 II1剖面模型进行阶地高程设计，分别为1 980、1 905、1 805、1 730、1 690、1 660 m，如图3所示。河谷发育过程，不仅仅对应着阶地高程，还有岸坡线的演化，这个演化完全是一个黑箱操作，因此只能进行一些推测。本文平移现今河谷岸坡线至 VI级阶地高程，作为初始岸坡线，下切过程中的岸坡线退化速度，与各阶地下切速度和经历时间成正比，最后设计结果见图 3。另外，下切过程中岸坡岩体的风化线演进，也是一个黑箱过程，本文设计的过程见表 2。表中，各数字代表着相应的组，共有57组，阿拉伯数字对应图3中的材料编号。各组位置见图 3，图中，阿拉伯数字为模型材料编号。

表2 模型下切过程及风化线演化过程表Table 2 The evolution process of bank slope frontier and belt of weathering in numerical model

相比较而言，文献[7，10，12]均未进行河谷下切模拟。文献[8]进行了河谷演化，先设计成平坦地质体，然后进行四级阶地演化，卸荷过程中岩体强度变化规律采用Mohr-Coulomb模型双线性应变软化遍布节理模型来实现。

4 应力场模拟

4.1 模拟方法

本文采用应力试算法，且前人研究[15－16]已得出结论，认为锦屏构造应力背景值在 7～12 MPa之间，这为本文的试算工作节省很多时间。采用 VI级阶地模型时，先用重力加速度算出重力场，再在边界上施加恰当的均布水平构造应力，然后进行河谷下切演化，得到现今构造应力场。由于控制河谷地区基本构造格局的构造运动发生在燕山期或更早，河流发育在晚第三纪才开始[17]，因此认为河谷发育期间，基本构造应力场保持不变，即整个下切过程中，施加边界的水平构造应力不变。另外，当河谷走向与区域地应力场中最大主压应力方位的锐角越大时，它对河谷地应力场分布的影响越大[6]。这种影响在应力试算法的数值模拟中如何体现？当模型边界与最大主压方向不一致时，在某一边界上施加的应力，可看作是水平构造应力的综合体现，但不是构造应力背景值。这里要提出的问题是，本文模型的X方向与区域主压方向一致，即与最大主应力方向一致，若只是在X方向施加构造应力，那么第二与第三主应力受构造应力的影响如何体现？通常的应力试算法，都只是在某一方向施加水平应力，如文献[7，10]均是在X边界上施加8 MPa的构造应力背景值。当笔者采用在 X边界上对称施加12 MPa（施加在VI级阶地模型的边界上，其他边界为固定边界）的构造背景应力值进行河谷下切演化时，发现此种边界条件下的现今河谷应力场只有最大主应力值与实测点应力数据相近，而第二、三主应力值均非常小。可见，当数值模型边界与区域主压方向一致时，只施加一个方向的边界应力，不足以很好地拟合真实的应力场。

文献[12]采用有限元回归分析法，经多次还原计算最终得到有限元模型中各节点的初始应力值。正如前述，该方法的准确性取决于测点数据的精度和个数，而且最后得到的是一个统一的回归方程。但是，对于谷坡应力场，不同位置，其应力场分布特征不一样，不可能用一个式子就能对其进行概括，因此认为该方法适合用于地形平坦，岩性较为均一的区域。文献[8]的初始模型为一平坦地面模型，通过一定的计算方法，得出各坐标方向上应力随深度变化的规律，见式(1)。

将该分布规律初始化到初始模型中，采用速度约束边界，对其进行河谷下切演化，得到现今河谷应力场。该文给出的初始应力分布规律与朱焕春等[18]提出的采用正交法得到的初始应力分布规律，见式（2）。

式中：T1、T2、K1、K2为常数；γ为岩石密度。

相类似，所不同的是式（2）中的竖向应力只是重力，不受构造应力的影响，而且水平向应力总是为最大、最小主应力，因此是一种潜在走滑型应力场。但竖直应力在多数情况下为最小主应力，在少数情况下为中间主应力，只有个别情况下为最大主应力[19]。而式(1)中的竖向应力，即Y方向应 力，不仅仅是重力，也受构造应力的影响，这与黄润秋等[5]提出的‘欧亚大陆板块垂直应力随深度增加梯度大于由岩石重度所确定的增加梯度’相一致。笔者认为文献[8]中初始应力的赋值更为合理。

综上，结合锦屏条件及本模型的构建情况，经多次试算后确定，在X方向对称施加12 MPa边界应力，Z方向对称施加8 MPa边界应力，Y方向底部采用固定边界。此种条件下得到的河谷现今应力场与实测值最为接近。

4.2 计算结果对比分析

II1勘探线前后有较多的地应力实测点，图 4为谷底钻孔 P122、P123、P124测点数据的最大水平应力分布图，各钻孔相对位置见图1。由图4可见，谷底应力包范围为河床下 80～145 m，约1 580～1 485 m 高程，最大水平应力为 22～37 MPa，且应力包偏向于右岸[20]。

图4 谷底测点数据最大水平应力分布图Fig.4 The scattergram on maximum horizontal stress with vertical depth of borehole at the bottom of valley

与钻孔 P122、P123、P124相对应位置处，截取了相应的数值模拟结果，即最大主应力随深度变化曲线。以25 MPa进行谷底“应力包”划定[16]，以内部应力变化斜率突变点为最深深度，各自的应力包范围分别为1 578～1 490 m，1 572～1 489 m，1 578～1 500 m，这与实测值非常一致。图 5为P123位置处的分布图及相应的测压力系数随深度变化曲线图。从图中曲线可以看出，模拟的应力包量值显然比实测值大，这是因为在河谷底部由于地形影响，导致构建数值模型时，产生一些畸形单元，这些单元会产生非常大的应力，可达上百MPa，如图6中蓝色的十字应力标，因此，FLAC3D软件在进行中间的插点计算时，这些畸形点周围的应力会显著增大，如图7的最大主应力云图，其最大值可达52 MPa。图6为各单元的应力标量图，为各个单元独立计算后显示的主应力，从图中可看出，排除极个别点后，应力包范围内的单元，最大主应力值范围为30～39 MPa，比实测值偏大10%～30%，属可接受范围。

图5 钻孔P123位置处的最大主应力与侧压力系数随深度变化曲线图Fig.5 Curves on the maximum principle stress and lateral pressure coefficient with vertical depth at the bottom of valley

本文所计算的河谷应力包与文献[8]计算的结力包，其最大主应力近水平，在25～35 MPa之间，应力包内及上部分布有应力值特别大的单元。

表3为左右岸实测点应力数据对比表，各编号位置见图中相应点所示。从表中可见，各测点数据与模拟数据，最大主应力和中间主应力十分接近，但是最小主应力差别较大。最小主应力之所以差别较大，笔者认为与泊松比的取值有关。计算过程中采用的是报告[14]中按岩体质量等级给出的参数，尤其是对位置较深部位的岩体，应按岩石参数进行取值计算，其结果或许更为合理。

5 锦屏现今河谷应力场特征

总体表现如下特征(在 FLAC3D软件中，以压为负，以拉为正，因此，在文中所有从该软件中得到的截图，即图6～8均是以压为负，以拉为正)。

（1）应力分异现象显著，且分异带宽度较宽，或者说是由坡表逐渐向内过渡，很难明确划分其范围（图6、8）。在该范围内，主应力方向发生偏转，越接近坡表，最大主应力越平行于坡面，最小主应力则越垂直于坡面。按照最大主应力与坡面的方向关系，将最大主应力与坡面基本平行的区域圈定为应力分异强烈带（见图8）。显然左右岸不对称，左岸受倾坡外的f42-9断层和黄斑岩脉影响，应力分异强烈带的深度显著加深。

（2）应力集中区主要表现为坡脚、谷底的1σ、τmax集中区，谷底应力包量值在30～39 MPa，略偏向于右侧（见图6、7），最大剪应力则达 20.7 MPa。

（3）应力分区（降低区、增高区、平稳区）现象较复杂，且左右岸情况不一致。模型中设置了1 620、1 700、1 750、1 800、1 870、1 950 m 高程的左右岸监测线，具体位置见图 1。从这些监测线的水平深度与最大主应力关系曲线图上看，左右岸应力分布不对称，左岸比右岸复杂。的水平深度与最大主应力关系曲线图上看，左右岸应力分布不对称，左岸比右岸复杂。

表3 实测点测量数据与模拟数据对比表Table 3 Comparison table between in-situ and simulation stress data

图6 锦屏现今地应力场主应力张量图(单位：10 MPa)Fig.6 The principle stress tensor map of Jinping current valley(unit: 10 MPa)

图7 锦屏现今地应力场最大主应力分布云图(单位：10 MPa)Fig.7 Nephogram of maximum principle stresses of Jinping current valley(unit: 10 MPa)

图8 现今地应力场主应力张量图(浅表放大部分)Fig.8 The enlarged principal stress tensor map at shallow

右岸1 750～1 950 m间的应力分布具有相似特征，如图9(a)为右岸1 800 m监测线上最大主应力与水平深度间的关系曲线，从图中可知，应力降低区与中卸荷线深度一致，水平深度约40 m。应力增高区深度约为280～610 m之间，最大主应力值达到最大，且变化平缓。在降低区与增高区之间有一个过渡区，区内，最大主应力随深度增加呈线性增大。应力增高区以内则为原岩应力区。1 750 m以下应力增高区的深度变浅，其宽度迅速减小并逐渐消失，如 1 700 m监测线上的应力增高区范围为225～300 m，过渡区内应力分布受f13、f14影响，波状起伏，而1 620 m监测线的最大主应力分布，受谷底应力包的影响，50～340 m的范围都可看作是应力增高区，过渡区几乎消失。左岸，在1 880 m高程（此处是f42-9与黄斑岩脉交汇处）以上，应力降低区的深度可划至 f42-9，以内则是应力增高区及逐渐过渡到原岩应力区，如图9(b)的1 950 m监测线的最大主应力与水平深度关系曲线所示。1 880 m以下至1 795 m段，f42-9在黄斑岩脉以外，该范围内的应力降低区深度划至中卸荷线，应力降低区与应力增高区之间有一个较宽的过渡带。

图9 最大主应力与水平深度关系曲线Fig.9 Relation curves of the maximum principal stress with horizontal depth

6 结 论

（1）进行河谷应力场反演模拟时，河谷模型概化的正确性是反演结果真实性的前提，而模型概化的细致程度则决定着，所模拟应力场能否体现研究区独特的应力场特征。本文所概化的模型，较前人所做的模拟都更为细致，认为本次模拟结果体现出了诸多独特的信息。

（2）采用试算边界应力法进行反演时，当边界与区域主压方向垂直时，只在该垂直边界施加水平应力不足以对第二、三主应力产生影响。建议在其他边界也施加适当的应力，或是改变边界方向。

（3）有着多级阶地的河谷，尤其是环青藏高原周边，下切迅速、地貌上呈V型谷的河谷，为更合理地反演应力场，对其进行河谷下切演化是很有必要的。

（4）河谷应力包的拟合度非常高，其他部位的实测数据则显示，第一、二主应力的拟合很好，但第三主应力稍微差些。

（5）左岸受倾坡外断层、黄斑岩脉影响，应力分布情况比右岸复杂，其应力分异强烈带、卸荷范围均比右岸深。应力分区现象，在不同高程体现出不同特征，左右岸不对称。

[1] 王涛, 周先前, 田树斌, 等.基于正交设计的河谷地应力场数值模拟方法及应用[J].岩土力学, 2003, 24(5):831－835.WANG Tao, ZHOU Xian-qian, TIAN Shu-bin, et al.Numerical simulation method for rock natural stress field of a valley and its application based on orthogonal experiments[J].Rock and Soil Mechanics, 2003, 24(5):831－835.

[2] 郭怀志, 马启超, 薛玺成, 等.岩体初始应力场的分析方法[J].岩土工程学报, 1983, 5(3): 64－75.GUo Huai-zhi, MA Qi-chao, XUE Xi-cheng, et al.The analytical method of the initial stress field for rock masses[J].Chinese Journal of Geotechnical Engineering, 1983, 5(3): 64－75.

[3] 张有天, 胡惠昌.地应力场的趋势分析[J].水利学报,1984, 15(4): 31－38.ZHANg You-tian, HU Hui-chang.Trend analysis of residual stress distribution in rock mass[J].Journal of Hydraulic Engineering.1984, 15(4): 31－38.

[4] 朱焕春, 赵海斌.河谷地应力场的数值模拟[J].水利学报, 1996, 27(5): 29－36.ZHU Huan-chun, ZHAO Hai-bin.Numerical simulation for rock stress field in a valley[J].Journal of Hydraulic Engineering, 1996, 27(5): 29－36.

[5] 黄润秋, 王士天, 张倬元.大型高边坡应力场数值模拟分析的理论与方法[J].成都地质学院学报.1990, (2):17－21.HUANG Run-qiu, WANG Shi-tian, ZHANG Zhuo-yuan.On the theory and methods of stress field numerical simulation of large high rock slopes[J].Journal of Chengdu Geology College, 1990, 17(2): 17－21.

[6] 朱焕春, 陶振宇, 黄德凡.河谷走向与河谷地应力分布[J].岩石力学与工程学报, 1995, 14(1): 17－24.ZHU Huan-chun, TAO Zheng-yu, HUANG De-fan.Influence of the trend on rock stress in valley[J].Journal of Rock Mechanics and Engineering, 1995, 14(1): 17－24.

[7] 伍法权, 祁生文, 宋胜伍, 等.复杂岩质高陡边坡变形与稳定性研究[M].北京: 科学出版社, 2008: 45－49.

[8] 荣冠, 朱焕春, 王思敬.锦屏一级水电站左岸边坡深部裂缝成因初探[J].岩石力学与工程学报, 2008, 27(增刊1): 2855－2863.RONG Guan, ZHU Huan-chun, WANG Si-jing.Primary research on mechanism of deep fractures formation in left bank of Jinping first stage hydropower station[J].Chinese Journal of Rock Mechanics and Engineering, 2008,27(Supp.1): 2855－2863.

[9] 祁生文, 伍法权, 兰恒星.锦屏一级水电站普斯罗沟左岸深部裂缝成因的工程地质分析[J].岩土工程学报,2002, 24(5): 596－599.QI Shen-wen, WU Fa-quan, LAN Heng-xing.Study on the mechanism of the deep fractures of the left abutment slope at the Jinping first stage hydropower station[J].Chinese Journal of Geotechnical Engineering, 2002,24(5): 596－599.

[10] 安关峰, 伍法权.锦屏水电站左坝肩岩体深卸荷带成因分析[J].岩土力学, 2003, 24(2): 300－303.AN Guan-feng, WU Fa-quan.Formation cause analysis of deep unloaded band of rock mass in left dam shoulder of Jinping hydropower station[J].Rock and Soil Mechanics,2003, 24(2): 300－303.

[11] 谢红强, 肖明砾, 何红达, 等.锦屏水电站坝区初始地应力场回归反演分析[J].长江科学院院报, 2008, 25(5):50－54.XIE Hong-qiang, XIAO Ming-li, HE Hong-da, et al.Regression analysis of initial geostress for dam region of Jinping hydropower station[J].Journal of Yangtze River Scientific Research Institute, 2008, 25(5): 50－54.

[12] 徐佩华, 陈剑平, 黄润秋, 等.雅砻江河谷下切三维数值模拟分析[J].吉林大学学报(地球科学版), 2003,33(2): 208－212.XU Pei-hua, CHEN Jian-ping, HUANG Run-qiu, et al.Three dimensions numerical simulation analysis of Yalongjiang valley[J].Journal of Jilin University(Earth Science Edition), 2003, 33(2): 208－212.

[13] QI S W, WU F Q, YAN F Z, et al.Mechanism of deep cracks in the left bank slope of Jinping first stage hydropower station[J].Engineering Geology, 2004, 73:129-144.

[14] 国家电力公司成都勘测设计研究院.雅砻江锦屏一级水电站可行性研究报告——近坝库岸及坝址区边坡稳定性研究[R].成都: 国家电力公司成都勘测设计研究院, 2003.

[15] 祁生文, 伍法权.高地应力地区河谷应力场特征[J].岩土力学, 2011, 32(5): 1460－1464.QI Sheng-wen, WU Fa-quan.Stress field characteristics of a valley slope in high geo-stress area[J].Rock and Soil Mechanics, 2011, 32(5): 1460－1464.

[16] 黄润秋.岩石高边坡发育的动力过程及其稳定性控制[J].岩石力学与工程学报, 2008, 27(8): 1525－1544.HUANG Run-qiu.Geodynamical process and stability control of high rock slope development[J].Chinese Journal of Rock Mechanics and Engineering, 2008,27(8): 1525－1544.

[17] 朱焕春, 陶振宇.地应力研究新进展[J].武汉水利电力大学学报, 1994, 27(5): 543－547.ZHUN Huan-chun, TAO Zhen-yu.New developments in study on rock stress[J].Journal of Wuhan University of Hydraulic and Electric Engineering, 1994, 27(5): 543－547.

[18] 朱焕春, 陶振宇.地形地貌与地应力分布的初步分析[J].水利水电技术, 1994, (1): 29－34.ZHUN Huan-chun, TAO Zhen-yu.A preliminary analysis on ground stress and topography morphology[J].Journal of Water Resources and Hydropower Engineering,1994, (1): 29－34.

[19] 张倬元, 王士天, 王兰生, 等.工程地质分析原理（第三版）[M].北京: 地质出版社, 2009.

[20] 谭成轩, 张鹏, 郑汉淮, 等.雅砻江锦屏一级水电站坝址区实测地应力与重大工程地质问题分析[J].工程地质学报, 2007, 16(2): 162－168.TAN Cheng-xuan, ZHANG Peng, Zhen Hanhuai, et al.An analysis on in-situ crustal stress measurements and major engineering geology issues at the dam site area of Jinping first stage hydropower station[J].Journal of Engineering Geology, 2007, 16(2): 162－168.