最小二乘支持向量机-粒子群算法在地下厂房围岩参数反分析中的应用

杨继华， 齐三红， 郭卫新， 张党立

(黄河勘测规划设计有限公司， 河南 郑州 450003)

0 引言

目前，国内外水电站地下厂房呈现出大跨度(20～30 m)、高边墙(40～80 m)的技术特点，部分断面面积大于1 000 m2，为超大断面地下工程。水电站大型地下厂房结构复杂，其受岩性、断层、节理、岩脉、地下水及地应力等因素的影响远大于普通隧洞[1-2]。

围岩的力学参数及地应力场特征(如弹性模量、泊松比、侧压力系数等)是地下厂房支护设计和稳定性分析的重要依据。岩体作为一种地质体，其成因极为复杂，且一般经历过多次构造作用，具有不连续性、非均质性和各向异性的特点，在力学性质上表现出明显的非线性特征[3]。

目前，岩体力学参数的确定主要有解析法、试验法及理论与实测相结合的方法等3种[4]。解析法先引入一定的假定条件，但假设往往与实际情况会有所差别，因此计算的参数存在偏差的可能性较大； 试验可分为室内试验和现场试验，但试验结果只能反映试验点处的岩体力学参数，无法代表整个工程区域的岩体特性，且试验费用昂贵、时间较长； 理论与实测相结合的方法一般利用实测的数据，采用理论分析的方法确定岩体的力学参数，又称为反分析法，目前常用的反分析法主要是位移反分析法。位移是反映岩土体受力变形的重要物理量，是岩土体稳定性最直接和最宏观的表现，位移反分析法融合了最优化理论和数值计算等先进方法和计算技术，具有理论基础完善、计算精度高、可以进行非线性分析、时间和经济成本低等优点，能弥补传统研究方法的不足，具有很大的实用价值[5-8]。刘英棨等[9]提出以水平收敛与拱底沉降相结合的多位移反分析法，建立平面简化有限元模型，实时分析隧道支护结构受力状态； 江宗斌等[10]提出了基于位移-应力的多元信息联合知能反分析方法； 张研等[11]将粒子群优化算法与高斯过程机器学习方法相融合，结合FLAC3D数值计算程序，提出隧洞围岩位移优化反分析的粒子群-高斯过程智能协同优化方法； 文建华等[12]采用分层运算方法对地下洞室进行了黏弹性位移反分析。综上所述，目前研究多集中在隧洞等小断面洞室方面，对大断面地下厂房涉及较少。另外，在位移反分析中，多采用监测的位移数据，但实际上在洞室开挖后位移监测仪器安装之前，大部分位移已经发生且无法监测到。因此，采用位移增量进行反分析更能反映围岩的物理力学性质。

最小二乘支持向量机做为一种新的机器学习方法，在处理小样本、非线性问题上有独特的优势。粒子群算法作为一种新型仿生进化算法，有全局优化、收敛速度快等优点[13]。本文以正交设计、最小二乘支持向量机和粒子群算法等现代数学方法为基本手段，建立围岩力学参数分析方法，以多点位移计实测位移增量值为依据，以厄瓜多尔Coca-Codo Sinclair水电站(简称CCS水电站)大型地下厂房为工程背景，采用二维弹塑性有限元方法建立地质结构分析模型，对CCS水电站地下厂房区域岩体力学特性及地应力场特征进行反分析。

1 最小二乘支持向量机与粒子群算法位移反分析方法

1.1 目标函数

位移反分析一般是通过某种方法，使得一组待分析的参数及其相应的位移值逼近实测位移值。对于工程实际，通过逼近达到总体上的最优效果，因此目标函数可用式(1)表示：

(1)

式中：x=(x1,x2,…,xk,…,xm)，xk为待反分析的岩体参数，如弹性模量E、黏聚力c、内摩擦角φ等，m为参数的个数；fi(x)为第i个测点上位移计算值；ui为相应的实测位移值；n为测点数。

1.2 最小二乘支持向量机

最小二乘支持向量机采用不同的优化目标函数，并且用等式约束取代不等式约束，可用于回归或模式识别等问题的解决[14-16]。假定训练集为(xi,yi)，其中：i=1, 2,…,k；xi∈Rn为n维系统输入向量；yi∈R为相应的输出值。如果问题为非线性，可将非线性映射输入向量映射到高维空间，转化为类似线性问题加以解决。最小二乘支持向量机优化问题的最小化函数为：

(2)

式中：w为权向量；C为惩罚因子，即调节常数。

约束条件为：

yi-wφ(xi)=b+ξi。

(3)

对于最小二乘支持向量机问题，根据式(2)和式(3)定义Lagrange求解方程：

(4)

式中ai为Lagrange乘子。

参数a和b的最优值可以通过KKT条件获得：

(5)

由式(5)可得：

(6)

消去式(6)中的w和ξ，优化问题转化为求解如下方程：

(7)

(8)

核函数K(x,xi)为任意对称函数，其应满足Mercer条件，常用的核函数主要包括线性函数、径向基函数及多项式函数等。最小二乘支持向量机的建模精度通常受2个因素的影响，即核函数参数与惩罚因子C的取值。本文采用粒子群算法确定核函数与惩罚因子的取值。

1.3 粒子群算法

粒子群算法(particle swarm optimization)是一种迭代优化工具[17-18]。在粒子群算法中，把问题的解看作为搜索空间中的粒子。所有粒子都有被优化函数决定的适应值，同时所有粒子通过速度决定其运动的方向和距离，其他粒子追随当前的最优粒子在解空间中搜索。粒子群算法首先产生一组初始化的随机粒子，随后采用迭代的方法寻求最优解。在迭代过程中，粒子通过对2个极值的跟踪更新自己。其中一个是粒子在每次搜索中的最优解，称为个体极值Pbest，另外一个是粒子群全部粒子在每次搜索中的最优解，称为全局极值gbest。粒子群中第i个粒子在n维空间的位置可用xi=(xi1,xi2,…,xin)表示，其速度可用vi=(vi1,vi2,…,vin)表示，第i个粒子的个体极值可表示为Pbest=(Pi1,Pi2,…,Pin)，粒子群的全局极值表示为gbest=(g1,g2,…,gn)。在搜索到这2个极值后，用式(9)和(10)来更新粒子的速度和位置。

vi(k+1)=uvi(k)+c1rand1(Pbest-xi(k))+

c2rand2(gbest-xi(k));

(9)

xi(k+1)=xi(k)+vi(k+1)。

(10)

式(9)—(10)中：c1、c2为学习因子，其取值在(0，2)之间； rand1和rand2为随机数，取值在(0，1)之间；u表示动量系数，其值随迭代改变。

1.4 最小二乘支持向量机和粒子群算法的位移反分析步骤

1)首先根据工程实际资料，确定待反分析的岩体力学参数及其取值范围，采用正交设计的方法确定计算方案；

2)建立有限元模型对全部计算方案进行计算，得到各监测点对应的位移值，将待反分析的岩体力学参数作为输入参数xi，计算位移值作为输出参数yi，形成学习的样本；

3)初始化粒子群算法，主要包括粒子群的规模、各粒子的权重因子、计算迭代次数、产生的随机粒子群向量及各粒子向量对应的最小二乘支持向量机惩罚因子C与核函数的参数σ2，把学习样本集作为训练样本和检验样本，将各粒子的个体极值设为当前位置，代入最小二乘支持向量机进行训练并获取对应的预测位移值；

4)计算各粒子相应的预测值与真实值的平均相对误差，并将其作为各粒子的适应值。然后进行迭代计算，并更新各粒子的速度和位置参数，同时记忆个体与全局所对应的最优适应值，满足初始设置的最大迭代次数时计算终止，最后记忆最优的惩罚因子和核函数参数(C，σ2)；

5)将上一步粒子群算法计算得到的最优参数(C，σ2)代入最小二乘支持向量机模型，建立待反分析岩体力学参数与位移值之间的非线性映射关系；

6)利用待反分析岩体力学参数与位移值之间的非线性映射关系代替正分析中的有限元计算，将位移反分析的目标函数值作为每个粒子的适应值，利用粒子群算法计算与位移实测值最接近的待反分析岩体力学参数。

2 工程应用

2.1 工程概况[19]

厄瓜多尔CCS水电站位于纳波省和苏昆比奥斯省境内的Coca河下游，为引水式电站，发电厂房为地下式，总装机容量为1 500 MW。地下厂房洞室群由主厂房、主变室、施工排水洞及引水压力管道等组成。主变室和主厂房2大洞室平行布置，洞室走向为NW45°，洞室断面为圆拱直墙型，其中主厂房开挖尺寸为212.8 m(长)×27.5 m(宽)×46.8 m(高)，主变室开挖尺寸为192.0 m(长)×17.0 m (宽)×34.0 m(高)。

2.2 工程地质条件

CCS水电站地下厂房埋深为300～400 m，区内未见较大规模的断层，但局部可见小断层及节理、裂隙密集带。地层岩性主要为侏罗系-白垩系迷撒华林地层(J-Km)紫红色、青灰色及浅红色的火山凝灰岩、火山角砾岩及流纹岩等，岩体结构以整体块状为主，局部呈块状、次块状。厂房区地下水以基岩裂隙水为主。厂房区内受构造应力影响不明显，为中—低应力值水平。

2.3 地下厂房施工方案

主厂房和主变室施工开挖采用钻爆法，由于主厂房和主变室的断面大，无法采用全断面开挖，故采用分层开挖的方法，其中主厂房分6层开挖，主变室分4层开挖，采用开挖一层支护一层的方法，支护方式为锚杆、挂网及喷混凝土。

2.4 有限元模型建立

地下洞室的数值模拟分析主要有二维模拟分析和三维模拟分析，二维平面分析的优点是可快速建模、计算速度快，对于CCS地下厂房，其纵向尺寸远大于横向尺寸，故可采用二维平面模型。地下厂房洞室群中除了主厂房和主变室之外，还有母线洞、尾水洞等洞室，但是在地下厂房基本开挖完成后才施工的，且其尺寸远小于地下厂房尺寸，所以对地下厂房影响较小，因此在本文的计算中不作考虑。

采用有限元程序Phase 2，根据地质概化模型建立8#机组剖面数值模型，计算范围如下：X方向为垂直主厂房和主变室轴线方向，以8#机组中心向上游方向延伸175 m、向下游方向延伸200 m，以指向下游方向为正；Y方向为竖直方向，以主厂房底板向下延伸90 m、向上延伸至地表面，以向上为正。计算范围大于5倍洞室，根据弹性力学理论，能满足计算边界的要求。模型考虑了对主厂房和主变室围岩稳定性有影响的F10、F11、F13、F14等小规模断层及地下水位。在数值计算过程中，断层采用接触单元模拟，断层参数未参与反分析，主要是因为断层开挖揭露后，对断层带物质进行了现场及室内试验，得出了其基本参数，数值计算中采用的断层基本参数如下： 切向刚度，3 500 MPa/m； 法向刚度，9 000 MPa/m； 黏聚力，0.10 MPa； 内摩擦角，25°。厂房区8#机组有限元模型共划分了6 391个节点、11 588个单元，网格划分后的有限元计算模型如图1所示。模型边界条件如下： 两侧为X向位移约束，底部为固定约束，顶部为自由边界。

图1 8#机组剖面有限元模型

2.5 位移增量

在CCS水电站主厂房和主变室的施工开挖过程中，为监测围岩的变形，选取了多个断面布置多点位移计，以8#机组断面为例，在洞室的项拱、拱角和边墙位置共布置了12条多点位移计测线，如图2所示。由于多点位移计是在洞室开挖到相应高程后安装的，因此其所测的位移值是该处多点位移计安装后围岩发生的位移，而测前发生的位移已被损失掉[20-21]，关于测前位移量损失的大小目前尚无好的确定方法，为避免反分析的随意性，本文采用主厂房和主变室第2层开挖完成与第1层开挖完成后的顶拱中心的多点位移计实测位移增量进行反分析计算。多点位移计孔口实测BX1-29增量位移值为0.48 mm，BX1-34增量位移值为0.35 mm。

图2 主厂房和主变室多点位移计布置(单位： m)

Fig. 2 Layout of multi-point extensometers of main powerhouse and main converter room (unit: m)

2.6 围岩参数的取值范围

根据CCS水电站地下厂房勘察设计资料及地下厂房围岩的实际情况，待反分析的围岩参数取值范围如下： 弹性模量E∈[9，25]，GPa； 泊松比μ∈[0.20，0.28]； 黏聚力c∈[0.60，2.20]，MPa； 内摩擦角φ∈[40，52]，(°)； 侧压力系数kx∈[0.30，0.70]、kz∈[1.30，1.70]。

2.7 反分析样本设计

理论上，反分析参数的组合有无穷多个，如果对每组参数都进行计算，实施困难，且没有必要，因此本文采用正交试验的设计方法，选取一定有代表性的组合进行计算。正交试验方法是根据数理统计学与正交性原理，从大量试验点中选取合适的有代表性的点，再按照“正交表”安排试验组合方案，由于正交表具有“均衡分散性”和“整齐可比性”的构造原则，因此，按照此方法设计的试验次数少，并且能反映客观事物的变化规律。正交表是正交试验设计的关键，它必须满足以下2个条件： 每一列(因素)的不同水平在试验中出现的次数相同，以保证其均匀性； 任意2列(因素)的不同水平组合组成的数对在试验中出现的次数相同，以保证试验点分布的均匀性。具体到本文的计算方案组合，围岩参数共有6个，在每个参数的取值范围内均匀选取5个值，各因素的取值范围并不是随意选取的，其中岩石(体)的力学参数是根据现场和室内试验得到的范围值，而地应力的侧压力系数是根据前期的地应力测试得到的范围值，如kx的范围值为0.3～0.7，这里将其分为5个水平，即0.3、0.4、0.5、0.6、0.7,以0.1为等级，能满足精度的要求。由于最小二乘支持向量机在处理小样本数据具有优势，25个不同围岩参数组合能代表围岩情况，满足参数反分析的要求，不同计算方案下的参数组合及有限元计算所得的位移增量见表1。

2.8 反分析结果及检验

在有限元计算中，输入不同的围岩参数，即可得出相应的位移值，即围岩力学参数与位移增量之间已有明确的映射关系，但实际的位移监测值与位移计算值并不对应，所以实际上的围岩参数是未知的，因此需要通过建立围岩力学参数与位移增量之间的模糊映射关系，根据实际的监测位移值采用逼近的方法来反分析围岩参数。

本文首先根据反分析的参数样本，进行有限元计算，各参数及位移增量值计算结果见表1。再根据最小二乘支持向量机和粒子群算法位移反分析理论和步骤，采用MATLAB编制相应的计算程序。最后将表1中的数据作为学习样本，建立围岩力学参数与洞室顶拱位移增量之间的非线性映射关系，按照1.3节的方法，经反分析得到围岩的力学参数反分析值如下： 弹性模量E=14.15 GPa，泊松比μ=0.23，黏聚力c=1.46 MPa，内摩擦角φ=49°，侧压力系数kx=0.55、kz=1.45。为检验反分析所得围岩力学参数的合理性，利用上述参数进行正分析有限元计算，得到主厂房和主变室分层开挖各多点位移计测点位移增量的计算值，选取主厂房第Ⅵ层与第Ⅰ层开挖、主变室第4层与第1层开挖的计算值与实测位移值增量进行对比，结果见图3和图4。可以看出，计算值和实测值相对误差较小，最大误差不超过10%，如BX1-34测点的误差为6.5%，BX1-35测点的误差为2.4%。在实际工程中，受监测仪器本身精度、施工方法、邻近爆破等的影响，监测结果会有一定的误差，但结果误差在10%以内是可以接受的，这也说明了反分析围岩参数的合理性。

表1 围岩力学参数及有限元计算结果

图3 主厂房多点位移计位移增量计算值与实测值对比

Fig. 3 Comparison of displacement increment between calculation results and monitoring values of main powerhouse

图4 主变室多点位移计位移增量计算值与实测值对比

Fig. 4 Comparison of displacement increment between calculation results and monitoring values of main converter room

3 结论与讨论

1)最小二乘支持向量机可有效降低计算的复杂性，加快求解速度，增强抗干扰能力，具有良好的优化性能，在处理小样本、非线性问题上有较大的优势。粒子群算法作为一种新型仿生进化算法，有收敛速度快、参数较少易确定等优点。最小二乘支持向量机和粒子群算法的有机结合为岩土工程位移反分析提供了一种新的方法。

2)本文以正交设计和有限元计算为基础，采用最小二乘支持向量机和粒子群算法相结合的方法，建立CCS水电站地下厂房洞室顶拱实测位移与围岩力学参数之间的非线性映射关系，成功地反分析出围岩弹性模量、泊松比、黏聚力、内摩擦角和侧压力系数。

3)多点位移计一般都是在洞室开挖后进行安装，因此无法测量到测前洞室发生的位移，由于确定测前位移损失较为困难，本文采用主厂房和主变室第2层和第1层开挖所产生的位移增量进行反分析，简化了反分析过程，有效提高了围岩参数反分析的精度。

4)对CCS水电站地下厂房反分析所得围岩参数进行检验，有限元计算结果表明： 主厂房第Ⅳ层与第Ⅰ层开挖和主变室第4层与第1层开挖所产生的位移增量计算值与多点位移计实测值吻合较好，最大误差不超过10%，说明了反分析参数的合理性，同时也说明了采用基于最小二乘支持向量机和粒子群算法的位移反分析方法在工程上是可行的，且效果较为理想。

5)在目前的地下工程岩体参数反分析中采用的多是位移反分析，实际上地下洞室开挖后，岩体的应力、位移、塑性区等均有一定程度的改变，在以后的研究中，可考虑采用应力、位移及塑性区等多元信息的综合分析，以便相互验证，提高分析精度。