D0:10.13374/.j.issn1001053x.2013.01.001 第35卷第1期 北京科技大学学报 Vol.35 No.1 2013年1月 Journal of University of Science and Technology Beijing Jan.2013 基于偏最小二乘法的地应力场拟合 陈章华)小，陈磊》，纪洪广2) 1)北京科技大学数理学院.北京1000832)北京科技大学土木与环境工程学院，北京100083 d通信作者，E-mail:chenzhanghua Gustb.edu.cn 摘要采用偏最小二乘回归分析方法分析矿区地应力场分布.以乌鞘岭隧道岭脊地段的有限测点的地应力测量结果为 样本，使用商业有限元软件ANSYS对岭脊地段区域的地应力场进行数值模拟.利用地应力实测数据与计算结果进行偏 最小二乘回归计算，得到拟合的地应力场数据.与传统最小二乘分析方法相比，本文的方法通过主成分分析，可以有效 解决多因变量回归计算时由于响应变量之间较高的相关性而导致的回归计算误差，并可以有效地反映变量间的关系，为 地应力场的建立提供有效的数据支持， 关键词偏最小二乘：地应力场：相关性：回归分析 分类号TD32 Geo-stress field fitting based on the partial least square method CHEN Zhang-hua),CHEN Lei 1,JI Hong-guang 2) 1)School of Mathematics and Physics,University of Science and Technology Beijing,Beijing 100083.China 2)School of Civil and Environmental Engineering,University of Science and Technology Beijing,Beijing 100083,China Corresponding author.E-mail:chenzhanghuaustb.edu.cn ABSTRACT Geo-stress distribution in a mining area was analyzed by means of a partial least squares regression method (PLSRM).Based on data measured in Wushaoling tunnel ridge section,a numerical simulation of the initial geo-stress field around the area was performed using the commercial software ANSYS.Partial least squares regression analysis was carried out through actual measurement data and calculation results of geo-stress,and then fitting data was ready.Compared with the traditional least squares method,due to the advantage of principal component analysis in applying PLSRM,the present method decreases the regression error caused by the correlation between response variables.It can more effectively reflect the features of discontinuous geo-stress distribution and thus provides effective data support for establishing the geo-stress field. KEY WORDS partial least squares;geo-stress field:correlation:regression analysis 随着统计方法的发展及工程建设的需要，分析参数取值的问题，提出了使用偏最小二乘(PLS)回 地应力场的方法越来越受到重视和应用.许多学者 归进行数据处理的理论，解决了变量间多重共线性 都对地应力场的分析提出具有建设性的理论方法.的问题，提出了一种适用于地应力场分析的数据处 蔡美峰等川就对玲珑金矿地应力场测量结果进行 理方法.在此基础上，聂善文针对上述方法耦合 分析并对地应力场进行建模.通过现场实测获得矿岩土参数本身的模糊性和随机性，对偏最小二乘算 区5个水平、12个测点的三维地应力状态.根据实法进行了进一步的完善，丰富了地应力场数据处理 测结果分析了矿区地应力场分布规律及其与地质构的方法.地应力场的分析在实际应用中也包括对矿 造的关系，建立了矿区地应力场分布模型，提高了 区的稳定性分析，马宇4针对金川主力矿山二矿区 测量结果的可靠性和精度.同时，在用数值方法对开展研究工作，通过对二矿区深部地质构造、岩石 地应力场进行计算的时候，还涉及到力学参数的求力学特性、构造应力场的现场调查和测试，以及室 解.马莎等2]针对建立岩土地应力场中关键的力学 内试验和数字仿真模拟，研究了矿区应力分布，在 收稿日期：201203-31 基金项目：国家高技术研究发展计划资助项目(2008AA062104)
第 卷 第 期 年 月 北 京 科 技 大 学 学 报 基于偏最小二乘法的地应力场拟合 陈章华 门, 陈 磊 , 纪洪广 北京科技大学数理学院, 北京 。。 北京科技大学土木与环境工程学院, 北京 困 通信作者 , 一, 一 川 、 · 摘 要 采 用偏 最小 二乘 回归分析方 法分析 矿区地应 力场分 布 以乌鞘 岭隧道 岭脊地 段的 有限测 点的地应 力测 量结 果为 样本, 使用商业有限元软件 对岭脊地段区域的地应力场进行数值模拟 利用地应力实测数据与计算结果进行偏 最小二乘回归计算, 得到拟合的地应力场数据 与传统最小二乘分析方法相比, 本文的方法通过主成分分析, 可 以有效 解决多因变量回归计算时由于响应变量之间较高的相关性而导致的回归计算误差, 并可以有效地反映变量间的关系, 为 地应 力场 的建 立提供 有效 的数据 支持 关键词 偏最小二乘 地应力场 相关性 回归分析 分类号 一 `刀百万 乳 了夕一人二 只 , `,万百 五 乞 , 万 万 夕一夕。 , `夕 场 , 、 , 、 , 郡 , , 困 一, 一 回 刀 仆 一 一 加 二、一 认、 朋 罗 、 , 、 、 、 ℃ 、 一 、 一 、 , , , , 一 , , 一 , , 品 、 旋 一 、 一 一 飞 、 一 、 、 随着统计方法的发展及工程建设的需要 , 分析 地应力场的方法越来越受到重视和应用 许多学者 都对地应 力场 的分析提 出具有建 设性的理论方法 蔡美峰等 川 就对玲珑金矿地应力场测量结果进行 分析并对地应力场进行建模 通过现场 实测获得矿 区 个水平 、 个测点的三维地应力状态 根据实 测结果分析了矿区地应力场分布规律及其与地质构 造的关系, 建立了矿区地应力场分布模型 , 提高了 测量结果的可靠性和精度 同时, 在用数值方法对 地应力场进行计算的时候, 还涉及到力学参数的求 解 马莎等图针对建立岩土地应力场中关键的力学 收稿 日期 。一 一 基金项目 国家高技术研究发展计划资助项 目 参数取值的问题 , 提出了使用偏最小二乘 回 归进行数据处理的理论 , 解决了变量间多重共线性 的问题, 提 出了一种适用于地应力场分析的数据处 理方法 在此基础上, 聂善文 网 针对上述方法祸合 岩土参数本身的模糊性和随机性, 对偏最小二乘算 法进行了进一步的完善, 丰富了地应力场数据处理 的方法 地应力场的分析在实际应用 中也包括对矿 区的稳定性分析 , 马宇叫针对金川主力矿 山二矿区 开展研究工作 , 通过对二矿区深部地质构造 、岩石 力学特性 、构造应力场的现场调查和测试 , 以及室 内试验和数字仿真模拟 , 研究了矿区应力分布 , 在 DOI :10.13374/j .issn1001 -053x.2013.01.001
.2 北京科技大学学报 第35卷 此基础上结合地下实际工程，分析了重点工程的稳 以利用有限元软件，在建模过程中对网格的划分 定性.现场地应力测量地点与工程所在位置，往往 方式、单元的类型以及计算方法进行不同的选择 有一定距离，因而测点处的地应力常常不能代表工以提高精度，但是得到的应力值只是有限元模拟 程所在位置的地应力状态.特别是在地形起伏及构 计算得到的数值，过于依赖建模过程中使用到的 造复杂的地区尤其如此.为此，王薇等问通过三维 岩土参数以及边界条件.此外，袁海平等川实现 地应力场的反演模拟，由测点处的地应力测量结果了基于工程有限测点实测应力的地应力场的自动反 推知工程所在地点一青藏铁路羊八井隧道轴线 演，发展和完善了地应力场反演方法.此种方法 上的地应力状态，为隧道的设计提供了地应力依据。依靠数学中的多元线性回归模型，反演得出应力 本文主要从偏最小二乘算法以及更高的地应 场：但是多元回归时，若自变量间相关性较高就 力拟合要求入手，通过回归计算，比较实测值、拟 会对结果造成无法克服的影响，可能导致在建立 合值与传统多元线性回归方法得到的拟合值之间的 自变量与因变量的拟合函数时，各因素系数的计 数据，为初始地应力场的建立提供较好的数据支持， 算产生误差.目前，经典的初始地应力场的反演 从而得出相应的结论 方法一般是基于有限元模拟计算的结果，再通过 多元线性回归，对计算值与实测值进行回归分析， 1偏最小二乘算法提出的背景 得到关于各个工况的拟合值，最后得到初始地应 仅依据工程中的实测结果往往无法表述整体 力场的估算分布，从而达到反演的目的.张勇慧 应力分布规律.目前普遍使用的回归反演方法是根 等⑧就利用此方法对大岗山水电站地下厂房区进 据地应力实测资料及地质构造条件，通过数值计算 行了三维地应力场反演 求得地应力最优回归系数，从而反演出相应的初始 总的来说，地层内部的应力由各方面因素共同 地应力场 作用影响，单纯依靠某种数学方法计算反演得到初 反演分析是岩土力学中的一个重要分支，因为 始地应力场显然不科学.下面具体介绍一般的建立 在工程实践中，材料参数的取值往往依赖于反演分 初始地应力场的步骤： 析.因此，根据现场量测信息（如地应力，地质材料 (1)根据现场地质地形勘测资料，利用有限元 弹性模量、泊松比和密度)来反演地应力数据以及 程序建立计算区域三维地质模型： 采用偏最小二乘回归方法对岩体参数进行分析，其 (2)根据地质力学分析，把影响初始地应力场 实用性为解决岩土工程计算参数的选取开辟了独特 形成的主要因素作为待定因素，对每种因素可利用 途径.采用偏最小二乘回归方法，可以解决样本点 数值计算得到已知量测点的应力值，然后在每种因 个数少、自变量间相关性强等问题，从而得到拟合 素已知量测点实测地应力值和计算应力值之间建立 值，用来建立自变量与因变量之间的关系.实际工 多元回归方程： 程中，利用地应力的计算值与实测值建立多元线性 (3)利用最小二乘法，根据残差平方和最小原 回归的模型是比较常见的方法.传统方法中直接利 则求得多元回归方程中各自变量系数的最优解，从 用计算值进行回归时，若不同要素的相关性太大， 而获得区域的初始地应力场. 会使回归计算时重复考虑同一种影响因素，从而使 回归函数的计算产生偏差.本文的方法是在计算应 3计算模型 力值的基础上，通过偏最小二乘算法，先得到独立 3.1数学模型 性较强的潜在变量再进行回归拟合，排除了传统方 偏最小二乘回归是一种新型的多元统计数据 法中自变量间相关性对结果造成的影响，得到的拟 分析方法，其优点在于回归建模过程中采用了信 合值更加可靠 息综合与筛选技术，不是直接考虑因变量集合与自 2 初始地应力场反演 变量集合的回归建模，而是在变量系统中提取若干 对系统具有最佳解释能力的新综合变量（又称为主 反演初始地应力场的作用是可以了解地应力 成分)，特别对于解释变量多重相关性及个数大于观 场的分布，可用于考察隧道围岩稳定性等问题.陈 察个体数的情况，是一种具有较好发展前景的新型 香梅等©根据地应力实测数据和有限元数学模型， 数据分析方法0-1刂.偏最小二乘法简单地说就是 提出了区域构造应变回归分析方法，反演了隧道围降维的最小二乘法.首先通过主成分分析的方法 岩初始地应力场.该方法虽然计算精度较高，可 提取变量中代表最大信息量的主成分，实现降维的
北 京 科 技 大 学 学 报 第 卷 此基础上结合地下实际工程 , 分析了重点工程的稳 定性 现场地应力测量地点与工程所在位置 , 往往 有一定距离, 因而测点处的地应力常常不能代表工 程所在位置的地应力状态 特别是在地形起伏及构 造复杂的地区尤其如此 为此 , 王薇等 通过三维 地应力场的反演模拟 , 由测点处的地应力测量结果 推知工程所在地 点 — 青藏铁路羊八井隧道轴线 上的地应力状态, 为隧道 的设计提供了地应力依据 本文主要 从偏最小二乘算法 以及更 高的地应 力拟合要求入手 , 通过回归计算, 比较实测值 、拟 合值与传统多元线性回归方法得到的拟合值之间的 数据, 为初始地应力场的建立提供较好的数据支持 , 从而得出相应的结论 偏最小二乘算法提 出的背景 仅依据工程 中的实测结果往 往无法 表述整 体 应力分布规律 目前普遍使用 的回归反演方法是根 据地应力实测 资料及地质构造条件, 通过数值计算 求得地应力最优回归系数 , 从而反演出相应的初始 地应力场 反演分析是岩土力学中的一个重要分支, 因为 在工程实践 中, 材料参数的取值往往依赖于反演分 析 因此 , 根据现场量测信息 如地应力, 地质材料 弹性模量 、泊松 比和密度 来反演地应力数据 以及 采用偏最小二乘回归方法对岩体参数进行分析 , 其 实用性为解决岩土工程计算参数的选取开辟 了独特 途径 采用偏最小二乘回归方法 , 可以解决样本点 个数少 、 自变量间相关性强等 问题 , 从而得到拟合 值 , 用来建立 自变量与因变量之间的关系 实际工 程中, 利用地应力的计算值与实测值建立多元线性 回归的模型是 比较常见的方法 传统方法中直接利 用计算值进行 回归时 , 若不同要素的相关性太大, 会使回归计算时重复考虑 同一种影响因素 , 从而使 回归函数的计算产生偏差 本文的方法是在计算应 力值的基础上 , 通过偏最小二乘算法 , 先得到独立 性较强的潜在变量再进行回归拟合 , 排除了传统方 法中 自变量间相关性对结果造成的影响, 得到的拟 合值更加可靠 初始 地应 力场 反 演 反演初始地应 力场的作用是可 以了解地 应力 场的分布 , 可用于考察隧道围岩稳定性等问题 陈 香梅等 根据地应力实测数据和有限元数学模型, 提 出了区域构造应变回归分析方法, 反演 了隧道围 岩初始地应力场, 该方法虽然计算精度较 高, 可 以利用有 限元软件 , 在 建模过程 中对 网格 的划 分 方式 、 单元的类 型 以及计 算方法进行不 同的选 择 以提 高精度 , 但是得到 的应力值 只是有 限元模拟 计算得到的数值 , 过于依赖建模过程 中使 用到的 岩土参数 以及边界条件 此外 , 袁海平等 实现 了基于工程有 限测点实测应力的地应力场的自动反 演 , 发展和完善 了地应 力场 反演 方法 此种方法 依靠数学中的多元线性 回归模型 , 反演得 出应 力 场 但是多元回归时 , 若 自变量 间相关 性较高就 会对结果造成无 法克服 的影响 , 可能导致在建立 自变量与因变量 的拟合 函数 时, 各 因素 系数的计 算产生误 差 目前 , 经典 的初始 地应力场 的反演 方法一般 是基于有 限元模拟计算 的结果, 再通过 多元线性回归, 对计算值与实测值 进行 回归分析 , 得到关于各个工况的拟合值 , 最后得到初始地应 力场 的估算分布 , 从而达到反演 的 目的 张勇慧 等 图 就利用此方法对大岗山水 电站地下厂房区进 行 了三维地应力场反演 总的来说 , 地层 内部的应力由各方面因素共同 作用影响 , 单纯依靠某种数学方法计算反演得到初 始地应力场显然不科 学 下面具体介绍一般的建立 初始地应力场的步骤 根据现场地质地形勘测资料 , 利用有限元 程序建立计算区域三维地质模型 根据地质力学分析 , 把影响初始地应力场 形成的主要因素作为待定因素, 对每种因素可利用 数值计算得到己知量测点的应力值 , 然后在每种因 素己知量测点实测地应力值和计算应力值之间建立 多元回归方程 利用最小二乘法 , 根据残差平方和最小原 则求得多元回归方程中各 自变量系数的最优解, 从 而获得 区域 的初始地应力场 计算模型 数学模型 偏最小二乘 回归是 一种新型 的多元统计数据 分析方法网, 其优点在于回归建模过程 中采用 了信 息综合与筛选技术, 不是直接考虑因变量集合与 自 变量集合的回归建模, 而是在变量系统 中提取若干 对系统具有最佳解释能力的新综合变量 又称 为主 成分 , 特别对于解释变量多重相关性及个数大于观 察个体数的情况, 是一种具有较好发展前景的新型 数据分析方法 `。一` 偏最小二乘法简单地说就是 降维的最小二乘法 首先通过主成分分析的方法 , 提取变量中代表最大信息量 的主成分 , 实现降维的
第1期 陈章华等：基于偏最小二乘法的地应力场拟合 3· 目的：然后用偏最小二乘回归方程，计算回归参数， 上施加法向线性分布荷载（如图1所示）作为荷载 得到残差矩阵，继续迭代运算，直到检验出的主成 条件来实现.通过计算不同荷载条件下的应力值并 分足够代表原变量：最后求得主成分与因变量的关 与实测应力值进行偏最小二乘法回归分析，最后得 系模型，继而得出回归模型 到拟合关系与拟合值，为最终回归初始地应力场提 3.2计算步骤 供有效的数据支持.这就是本文中具体的回归拟合 具体的计算步骤如下 步骤.与一般的传统多元线性回归方法不同的是， (1)对数据进行标准化处理 本文提出利用偏最小二乘法的数学方法处理数据、 至于此方法的优势在之后的结果分析中进行阐述 e (en -ei)/Si (1) 4.2确定自变量与因变量 ,表示标准化之后的值，e表示原始值，E是样 以乌鞘岭隧道工程某些特定测点的应力回归 本的平均值，S,是样本的标准差.对数据进行标准 求取为例进行分析.根据地质资料，在计算区域内 化处理的目的是使样本点集合的重心与坐标原点重 考虑的主要岩性有砂岩夹砾岩及泥岩、砂岩夹页层 合，压缩处理可以消除由量纲不同而引起的虚假变 及薄层煤、板岩夹千枚岩、安山岩等，主要考虑的 异信息，使分析结果更加合理.为方便起见，仍记 断层带中不同的岩体力学参数如表1所示. 标准化处理的矩阵为X. 表1岩体的力学参数表 (2)计算标准化数据矩阵X的协方差矩阵V, Table 1 Mechanical parameters of the rock 这时V又是X的相关系数矩阵. 介质序号围岩级别弹性模量/GPa泊松比密度/(gcm-3) (3)求V的前m个特征值，以及对应的特征 1 IV 1.80 0.35 2.50 向量，要求它们是标准正交的 2 1.70 0.45 2.40 (4)求第h主成分F,有 IV 1.80 0.35 2.50 V 1.70 0.45 2.40 F=Xah=∑anjt (2) 5 IV 1.80 0.35 2.50 式中，a是主轴ah的第j个分量.所以，主成分 本次实测值的选取来自探井5#.对5#井来 F是原变量x1,x2.·,xp的线性组合，组合系数 说，所处的地质层属于上表中第一个，即序号1的 恰好为a,从这个角度，又可以说F是一个新的 介质.所以首先模拟5#井所在的1000m×1400 综合变量. m×1500m地质层模型，利用ANSYS软件，使用 (5)具体地，在得到了第一主成分以后，偏最小 SOLD185单元，模拟出包含约10万个单元的模 二乘算法要求通过多次的运算迭代，找出最能代表 型.图1表示具有边界线性分布侧压的网格图，图 自变量信息的所有成分 2表示计算得到的第三主应力云图 简单地说，每次通过回归求得残差矩阵，就能 然后可以根据5#井的位置，在有限元模型中， 得到下一个成分.最后根据有效性检验，把满足条 得出实测点位置的计算应力值，如表2所示.根据 件的所有成分提取出来，以所有成分为新的自变量， 反演地应力场的方法，将表2中的计算应力值作为 因变量不变，建立它们之间的关系模型.这样就能 自变量，待之后的偏最小二乘回归计算.同时，在各 根据所得的模型，与原始实测数据进行比较. 实测点处，已经测得了最大主应力、最小主应力等 数据，如表3所示.将这些数据提取出，作为因变 4实例计算分析 量，以便之后的计算 4.1分析策略 把计算应力值作为自变量，实测值作为因变量 在地矿工程中，为获得地下岩土材料和应力分 进行回归.由于在使用ANSYS建模，形成如图1 布数据，必须进行地下取样和地下埋点实测.测得 的地质层时，需要输入表1中1号介质层的各项数 数据包括岩土层的弹性模量、泊松比、特定点的主据作为建模的参数，因此此模型中已经包含影响地 应力、主方向等.反演初始地应力场的过程就是利应力的材料力学常数，包括弹性模量E、泊松比“ 用这些数值尽可能准确地得到整个地质层的应力分 等.回归系数用于调整有限元计算所施加的边界位 布情况.根据力学知识可知，弹性模量、泊松比、密 移、边界线性分布载荷集度、自重场的重度等参数 度等数据是计算三维地应力场不可缺少的条件.在 在得到了计算数据时，可以看出自变量间存在 本研究中，选取岩体自重、地质构造位移和在边界 着多重相关性，在采用回归方法的时候，要选择可
第 期 陈章 华等 基 于偏最 小二 乘法 的地应 力场 拟合 目的 然后用偏最小二乘回归方程 , 计算 回归参数, 得到残差矩 阵, 继续迭代运算, 直到检验出的主成 分足够代表原变量 最后求得主成分与因变量的关 系模型, 继而得 出回归模型 计算步骤 具体的计算步骤如下 对数据进行标准化处理, 布, 。`了一引 匀 瓦, 表示标准化之后的值 , 。, 、表 示原始值, 弓 是样 本的平均值 , 是样本的标准差 对数据进行标准 化处理的 目的是使样本点集合的重心与坐标原点重 合 , 压缩处理可 以消除 由量纲不 同而引起的虚假变 异信息 , 使分析结果更加合理 为方便起见, 仍记 标准化处理 的矩阵为 计算标准化数据矩阵 的协方差矩阵 , 这时 又是 的相关系数矩阵 求 的前 。 个特征值, 以及对应的特 征 向量 , 要求它们 是标准 正交的 似 求第 主成分 月, , 有 上施加法向线性分布荷载 如 图 所示 作为荷载 条件来实现 通过计算不同荷载条件下 的应力值并 与实测应力值进行偏最小二乘法回归分析 , 最后得 到拟合关系与拟合值 , 为最终回归初始地应力场提 供有效的数据支持 这就是本文 中具体 的回归拟合 步骤 与一般的传统多元线性回归方法不 同的是 , 本文提 出利用偏最小二乘法的数学方法处理数据 , 至于此方法的优势在之后的结果分析 中进行阐述 确定 自变量与因变量 以乌鞘 岭隧道工程某些特 定测 点的应力 回归 求取为例进行分析 根据地质 资料 , 在计算区域 内 考虑的主要岩性有砂岩夹砾岩及泥岩 、砂岩夹页层 及薄层煤 、板岩夹千枚岩 、安山岩等, 主要考虑的 断层带 中不同的岩体力学参数如表 所示 表 岩体的力学参数表 恤 人 飞 介质序号 围岩级别 弹性模量 泊松 比 密度 。一丁勺 几一 ,, 一艺 、 , · 生 〕 刁 式中, 勺、 是主轴 的第 , 个分量 所以, 主成分 凡, 是原变量 二, ·… 二。的线性组合, 组合系数 恰好为 嘶 从这个角度 , 又可以说 几 是一个新的 综合变量 具体地 , 在得到了第一主成分 以后 , 偏最小 二乘算法要求通过多次的运算迭代 , 找出最能代表 自变量信息的所有成分 简单地说, 每次通过回归求得残差矩阵, 就能 得到下一个 成分 最后根据有效性检验 , 把满足条 件的所有成分提取 出来 , 以所有成分为新的自变量 , 因变量不变, 建立它们之间的关系模型 这样就能 根据所得的模型, 与原始实测数据进行 比较 实例 计 算分 析 分析策略 在地矿 工程中 , 为获得地下岩土材料和应力分 布数据, 必须进行地下取样和地下埋点实测 测得 数据包括岩土层的弹性模量 、泊松 比 、特定点的主 应力 、主方 向等 反演初始地应力场 的过程就是利 用这些数值尽可能准确地得到整个地质层的应力分 布情况 根据力学知识可知, 弹性模量 、泊松 比 、密 度等数据是计算三维地应力场不可缺少的条件 在 本研究中, 选取岩体 自重 、地质构造位移和在边界 本 次实测值的选取来 自探井 对 井来 说 , 所处的地质层属于上表中第一个 , 即序号 的 介质 所以首先模拟 井所在的 。 地质层模型, 利用 软件, 使用 单元 , 模拟 出包含约 万个单元的模 型 图 表示具有边界线性分布侧压的网格图, 图 表示计算得到 的第三主应力云图 然后可 以根据 井的位置 , 在有 限元模型 中, 得 出实测点位置 的计算应力值 , 如表 所示 根据 反演地应力场的方法, 将表 中的计算应力值作为 自变量, 待之后的偏最小二乘回归计算 同时 , 在各 实测点处, 已经测得了最大主应力 、最小主应力等 数据 , 如表 所示 将这些数据提取 出, 作为因变 量, 以便之后的计算 把计算应力值作为 自变量 , 实测值作为因变量 进行回归 由于在使用 建模, 形成如 图 的地质层时 , 需要输入表 中 号介质层的各项数 据作为建模的参数, 因此此模 型中己经包含影响地 应力的材料力学常数 , 包括弹性模量 、泊松 比 户 等 回归系数用于调整有限元计算所施加的边界位 移 、边界线性分布载荷集度 、自重场的重度等参数 在得到了计算数据时, 可 以看出 自变量间存在 着多重相关性 , 在采用回归方法的时候, 要选择可
北京科技大学学报 第35卷 以减少相关性干扰的方法.否则，就会跟传统的多 系，同时得到回归的拟合值，用以分析比较.从而 元回归方法一样，对拟合结果产生干扰。为此，本 为反演得到初始地应力场提供较为准确的数据.将 文提出使用偏最小二乘法，利用主成分分析和有效 自变量数据构成的矩阵，记为E:将因变量数据构 性检验结合的方法，得出自变量和因变量之间的关 成的矩阵，记为F. N一 分15：32 0.755E+07 0.118E+07 0.992E+07 -0.187E+8 -0.274E+08 -0.319E+07 -0.555E+07 -0.143E+08-0.230E+08 0.318E+08 图1三维有限元模拟后的地应力场示意图（单位：Pa) Fig.1 Schematic diagram of the geo-stress field after three-dimensional finite element simulation (unit:Pa) AN 1B29212 09:1n:b2 0.317E+08 -0.209E+08 -0.101E+08 710996 .115E+0 .263E+0 0.155E+08 -0.469E+07 0.613E+07 0.170E+08 图2三维有限元模拟后的地层应力云图（单位：Pa) Fig.2 Nephogram of geo-stress after three-dimensional finite element simulation (unit:Pa)
第1期 陈章华等：基于偏最小二乘法的地应力场拟合 .5 表25#井各测点计算应力值 过程序计算，最后得到的应力拟合值如表4所示. Table 2 Calculated stress values of measuring points in Well 5妆 开始 X向应力/ Y向应力/ 2向应力/ XY向切应力/ MPa MPa MPa MPa -17.2 -2.55 -2.59 -2.26 输入自变量与因变量 -17.5 -2.62 -2.51 -2.38 矩阵 -17.5 -2.62 -2.51 -2.38 -17.1 -2.60 -2.25 -2.62 -17.7 -2.67 -2.35 -2.60 标准化操作 -17.7 -2.67 -2.35 -2.60 -17.8 -2.75 -2.10 -2.78 -17.9 -2.72 -1.71 -2.78 求残差矩阵.并 求协方差矩阵 -15.5 -2.64 -1.23 -2.58 代替原矩阵 表35#井各测点实测应力值 得第一主成分 Table 3 Measured stress values of measuring points in Well 不 满 5# 足 水平面内最大/ 水平面内最小/ 自重应力/ 主应力MPa 主应力MPa MPa 有效性检验 -4.50 -2.50 -13.98 -11.50 -6.00 -14.07 -9.59 -6.09 -14.11 警 -23.00 -12.00 -14.15 -17.00 -8.50 -14.17 -21.00 -10.50 -14.27 得所有主成分 -16.00 -8.00 -14.32 -11.69 -6.69 -14.38 -19.31 -11.31 -14.70 主成分与因变量关系 4.3应用偏最小二乘回归得到拟合值 (1)首先对自变量E和因变量F进行标准化 自变量与因变量关系 处理，以方便它们在同一量纲上进行运算操作，得 到E和F 因变量拟合值 (2)通过协方差矩阵的计算，求出最大特征值 及其所对应的单位特征向量，继而求出主成分 图3主要程序流程图 (3)根据最小二乘迭代和普遍适用且简便的有 Fig.3 Program flow chart 效性判别法（特征根大于1法），筛选出所有主成分 表4偏最小二乘拟合值 主成分t是根据E和w的矩阵相乘而来，其中w Table 4 Stress fitting by PLS 是指特征向量组成的矩阵 最大主应力/MPa 最小主应力/MPa 自重应力/MPa (4)之后，利用最小二乘回归原理，得到主成分 实测值拟合值 实测值拟合值 实测值拟合值 与因变量的关系.此时，可求出因变量的拟合值以 -4.50 -13.51 -2.50 -7.23 -13.98 -13.99 及因变量与自变量的关系式，最后将标准化之后的 -11.50 -10.86 -6.00 -5.90 -14.07 -14.08 -9.59 -10.86 -6.09 -5.90 -14.11 -14.08 因变量的值，进行标准化反向计算，得到实际的因 -23.00 -24.79 -12.00 -12.58 -14.15-14.17 变量拟合值 -17.00 -14.79 -8.50 -7.68 -14.17 -14.18 由于偏最小二乘算法中涉及的原始数据和算 -21.00 -14.79 -10.50 -7.68 -14.27 -14.18 法的使用要求较高的精度，因此使用计算机编程的 -16.00 -13.90 -8.00 -7.32 -14.32 -14.36 -11.69 -10.82 -6.69 -6.09 -14.38-14.44 算法可以又准又快地得出拟合结果.考虑到较多的 -19.31 -19.25 -11.31 -11.22 -14.70-14.67 涉及矩阵的运算，运用MATLAB编程计算可以快 速方便地得出结果.主要的程序流程如图3所示.通 在求出偏最小二乘拟合值的同时，可以拟合出
第 期 陈章华等 基于偏最小二乘法的地应力场拟合 表 井各测点计算应力值 过程序计算 , 最后得到的应力拟合值如表 所示 认 、 三 向应力 向应力 向应力 一 一 一 一 一 一 ` 一 一 一 一 向切应力 一 一 一 了 一 一 一 一 一 月 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 开始 输人 自变量 与因变量 矩 阵 标准化操作 表 井各测点实测应力值 入 、 水平面内最大 水平面 内最小 自重应力 主应力 主应力 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 召 一 一 求残差矩阵 并 求协方差矩阵 代替原矩阵 得第一主成分 应用偏最小二乘回归得到拟合值 首先对 自变量 和 因变量 进行标准化 处理 , 以方便它们在 同一量纲上进行运算操作, 得 到 和 通过协方差矩阵的计算 , 求出最大特征值 及其所对应 的单位特征向量, 继而求 出主成分 根据最小二乘迭代和普遍适用且简便 的有 效性判别法 特 征根大于 法 , 筛选 出所有主成分 主成分 是根据 和 叨 的矩阵相乘而来 , 其 中 叨 是指特征向量组成的矩阵 之后, 利用最小二乘回归原理, 得到主成分 与 因变量的关系 此时, 可求出因变量的拟合值 以 及因变量与 自变量的关系式 , 最后将标准化之后的 因变量 的值, 进行标准化反 向计算, 得到实际的因 变量拟合值 由于偏 最小二乘算法 中涉及 的原始数据和算 法的使用要求较高的精度 , 因此使用计算机编程 的 算法可 以又准 又快地得出拟合结果 考虑到较多的 涉及矩阵的运算 , 运用 ' , 编程计算可 以快 速方便地得出结果 主要的程序流程如图 所示 通 得所有主成分 主成分与因变量关系 自变量 与因变量关 系 因变量拟合值 图 主要程序流程图 表 偏最小二乘拟合值 妞 最大主应力 实测值 拟合值 一 一 一 一 一 一 名 一 一 一 一 一 一 一 一 一 习 一 一 一 一 最小主应力 实测值 拟合值 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 ` 一 一 习 一 一 , 自重应力 实测值 拟合值 一 一 一 一 一 一 , 一 一 一 一 一 一 一 一 一 一 一 一 乃 在求出偏最小二乘拟合值的同时, 可以拟合出
北京科技大学学报 第35卷 自变量与因变量之间的解析表达式： -13.0 -13.2 一◆一实测值 Y=XRx+Co (3) -13.4 一偏最小二乘值 -13.6 多元回归值 其中X与Y分别为两种变量的矩阵，Rx为系数 -13.8 矩阵，C为常数项 -14.0 -7.5779 -3.8794 -0.0699 念 -14.2 -80.4147 -37.2359 1.2116 -14.4 Rx= 12.1325 5.0290 -0.3268 -14.6 51.8505 23.7156 0.0128 -14.8 5 6 C为9行3列矩阵，每一行的值相同，每一列的 待测点 值依次为200.3128、102.2844和12.9189 图5自重应力拟合值与实测值比较 此时，根据三类主应力与方位角、倾角的拟合 Fig.5 Comparison between fitted gravity stress and mea- 分布值，就可得到相应的x、,和z方向应力与各 sured gravity stress 切向应力的拟合值，也就是初始地应力场的建立 4.4传统多元线性回归与PLS拟合比较 5结论 传统的地应力场回归方法是多元线性回归的 方法.下面，根据赵德安等在文献12]中得到的 (1)本文通过有限元数值模拟、实测应力并结 通过多元线性回归计算的数据，即图4和图5中 合偏最小二乘回归方法对乌鞘岭隧道区域岭脊地段 所示的多元回归值，与使用偏最小二乘方法得到的 地应力场的分布特点进行了拓展研究，得到了较为 数据进行比较.为了使比较更为直观，把待测点的 合理的应力拟合值.相对于多元线性回归方法，本 实测值、偏最小二乘拟合值和多元线性回归拟合值 文使用的方法在回归准确性方面有所改进， 用折线图的形式表现出来.为了不造成冗余，只对 (②)反演初始地应力场时，首先利用有限元方 最大主应力与自重应力两类数据进行比较，分别如 法对指定地质层建模，模型中采用该地质层的岩土 图4和图5所示.由两图可以看出，多元线性回归 材料常数并施加反映构造运动的位移边界条件、线 得到的拟合值，与实测值相比差距较大，而偏最小 性分布的边界应力和自重.由此得到计算应力值， 二乘拟合回归得到的值更为准确地接近实测值.图 以便进行拟合回归 4中最大主应力的偏最小二乘波动曲线很好地复合 (③)通过地下实测应力，获得测试样本点数据 了实际的数值走向，这是因为它去除了自变量间相 是反演初始地应力场非常重要的工作内容.通常由 关性的影响，使结果比多元回归更加准确.在图5 于成本限制，实测样本点数相对很少，测试得到的 的自重应力比较中，由于提高了数据的精度再进行 有效数据非常有限.反演初始地应力场所涉及变量 比较，可以看出，偏最小二乘依然很好地拟合了实 多，且一些变量之间存在多重线性相关性.在这种 际的结果，而多元回归与之相比就要差上不少 情况下，运用偏最小二乘回归方法进行回归分析非 常有效.传统的多元线性回归方法实质上是对回归 →一实测值 ·一偏最小二乘值 区域内的数据做了较为平滑的拟合，但不能在反映 ·一多元回归值 局部区域内，由于材料分布不同而引起的应力响应 突变.偏最小二乘法不仅考虑了影响应力分布的材 料常数、边界作用、自重影响以及实测结果，同时 20 还检测并最大可能的剔除了自变量间的相关性，可 以很好地反映应力响应突变的规律，得到的拟合值 与实测数据相比，减少了误差，数据走向也符合实 456 际情况. 待测点 (4)本文所用的自变量是利用有限元软件AN 图4最大主应力拟合值与实测值比较 SYS计算所得的与地下实测应力点位置相对应的计 Fig.4 Comparison between fitted maximum principal stress 算点的正应力和切应力.自变量可形成n行4列的 and measured maximum principal stress X矩阵，n表示样本点个数.当xy面与地面平行
北 京 科 技 大 学 学 报 第 卷 自变量与因变量之间的解析表达式 十 今 实测值 偏最小二乘值 一“ 多元回归值 皿逻州只侧山`川 其 中 与 分别为两种变量的矩阵, 矩阵, 为常数项 为系数 一 专 一 , 一 一 一 一 一 一 一 一 卫 一 一 石 一 名 一 一 一 滩 一 一 为 行 列矩阵, 每一行的值相 同, 每一列 的 值依次为 、 和 此时, 根据三类主应力与方位角 、倾角的拟合 分布值 , 就可得到相应的 、 和 方向应力与各 切 向应力的拟合值, 也就是初始地应力场的建立 传统多元线性回归与 拟合比较 传 统的地应力场 回归方法 是多元线性回归的 方法 下面 , 根据赵德安等在文献 中得到的 通过多元线性 回归计算 的数据 , 即图 和 图 中 所示的多元 回归值 , 与使用偏最小二乘方法得到的 数据进行 比较 为 了使 比较更为直观 , 把待测点的 实测值 、偏最小二乘拟合值和多元线性回归拟合值 用折线 图的形式表现出来 为 了不造成冗 余, 只对 最大主应力与 自重应力两类数据进行 比较 , 分别如 图 和图 所示 ` 由两图可以看出, 多元线性回归 得到的拟合值, 与实测值相比差距较大, 而偏最小 二乘拟合回归得到的值更为准确地接近实测值 图 中最大主应力的偏最小二乘波动 曲线很好地复合 了实际的数值走 向, 这是因为它去除了自变量 间相 关性的影 响, 使结果 比多元 回归更加准确 在 图 的 自重应力比较中, 由于提高了数据的精度再进行 比较 , 可 以看 出, 偏最小二乘依然很好地拟合 了实 际的结果, 而多元回归与之相比就要差上不少 待测点 图 自重应力拟合值与实测值比较 · 一 一 实测值 偏最小二乘值 门曰目口氏︸ 书一 只侧州拭唱岌川 待测点 图 最大主应力拟合值与实测值 比较 结论 本文通过有 限元数值模拟 、 实测应力并结 合偏最小二乘回归方法对乌鞘岭隧道 区域岭脊地段 地应力场的分布特点进行 了拓 展研究, 得到了较为 合理的应力拟合值 相对于多元线性 回归方法 , 本 文使用 的方法在回归准确性方面有所改进 反演初始地应力场 时, 首先利用有限元方 法对指定地质层建模, 模型中采用该地质层的岩土 材料常数并施加反映构造运动 的位移边界条件 、线 性分布 的边界应力和 自重 由此得到计算应力值, 以便进行拟合回归 通过地下实测应力 , 获得测试样本点数据 是反演初始地应力场非常重要的工作 内容 通常 由 于成本限制, 实测样本点数相对很少 , 测试得到的 有效数据非常有限 反演初始地应力场所涉及变量 多, 且一些变量之间存在多重线性相关性 在这种 情况下, 运用偏最小二乘回归方法进行回归分析非 常有效 传统的多元线性 回归方法实质上是对 回归 区域内的数据做了较为平滑的拟合 , 但不能在反映 局部区域内, 由于材料分布不 同而 引起的应力响应 突变 偏最小二乘法不仅考虑 了影响应力分布 的材 料常数 、边界作用 、 自重影响 以及 实测结果 , 同时 还检测并最大可能的剔除了 自变量间的相关性, 可 以很好地反映应力响应突变 的规律, 得到的拟合值 与实测数据相比, 减少 了误差 , 数据走 向也符合实 际情况 本文所用 的自变量是利用有 限元软件 计算所得的与地下实测应力点位置相对应的计 算点的正应力和切应力 自变量可形成 行 列的 矩阵, 表示样本点个数 当 夕面与地面平行
第1期 陈章华等：基于偏最小二乘法的地应力场拟合 7· 时，4表示采用影响权重比较大的x向正应力、y向 dimensional crustal stress field and its application in tun- 正应力、之向正应力以及xy切应力共4个应力分 nel design.Acta Geosci Sin,2004,25(5):71 量.若采用边界条件的数值不同，仅影响如表2所 (王薇，王连捷，乔子江，等.三维地应力场的有限元模拟及 示的作为自变量的计算值.利用偏最小二乘回归分 其在隧道设计中的应用，地球学报，2004,25(5)：71) 析所得的拟合值不会发生太大改变，以此拟合值去 [6]Chen X M,Ren Y Q,Xu Z Y,et al.The inversion of 校正指定点的应力值，因此最终拟合的结果必然趋 initial geostress field of tunnel surrounding rock based on 近实际分布规律 FEM model.Chin J Eng Geophys,2011,8(6):767 (陈香梅，任永强，徐智勇，等.基于有限元数学模型的隧道 (5)由于反演初始地应力场要求能够获得地质 围岩初始地应力场反演.工程地球物理学报，2011,8(6)： 层内指定点的应力，当采用偏最小二乘回归分析时， 767) 首先根据有限元软件得到场中此指定点的计算应力 [7 Yuan H P.Zhao K,Wang J L,et al.On the method of 值，然后采用偏最小二乘回归所得的拟合系数可以 multivariate linear regression for geo-stress field.Nonfer- 确定指定点的应力拟合值，从而确定在此点处工程 rous Met Sci Eng,2011,2(1):59 计算所需要的各种应力值. (袁海平，赵科，王继伦，等.多元线性回归地应力场反演方 法研究.有色金属科学与工程，2011,2(1)：59) (8]Zhang Y H,Wei Q.Sheng Q.et al.Three dimensional 参考文献 back analysis of geostress field in underground powerhouse zone of Dagangshan hydropower station.Rock Soil Mech, [1]Cai M F,Wang S H.Results and their analysis of in-side 2011,32(5:1523 earth stress measurement in Linglong gold mine.China (张勇慧，魏侍，盛谦，等。大岗山水电站地下厂房区三维地 Mm,2000,9(5):46 (蔡美峰，王双红.玲珑金矿地应力场测量结果及其分析 应力场反演分析，岩土力学，2011,32(5)：1523) 中国矿业，2000,9(5)：46) [9 Wang H W.Partial Least Squares Regression Method and 2]Ma S,Wu L F,Zhao M X,et al.The value analysis Its Application.Beijing:National Defense Industry Press 1999 of the rock deformation module according to the partial least-squares regression.J North China Inst Water Con- (王惠文.偏最小二乘回归方法及其应用.北京：国防工业 servancy Hydroelectr Power,2005,26(3):51 出版社，1999) (马莎，吴林峰，赵明献，等.偏最小二乘回归在岩体力学 [10]Wold S.Trygg J,Berglund A,et al.Some recent develop- 参数取值中的应用.华北水利水电学院学报，2005,26(3)： ments in PLS modeling.Chemom Intell Lab Syst,2001 51) 58(2):131 [3]Nie S W.Geotechnical mechanics parameters based on [11]Lu Q H,Ren K L,Zhou F X.Achieve non-linear regres- partial least squares regression method and random-fuzzy sion analysis based on partial least squares.Gansu Sci theory.J Guilin Univ Technol,2010,30(1):80 Technol,,2005,21(11):146 (聂善文，基于偏最小二乘法和随机模糊理论的岩土力学 (鲁庆华，任康乐，周风玺，基于偏最小二乘法实现非线性 参数研究.桂林理工大学学报，2010,30(1)：80) 回归分析.甘肃科技，2005,21(11：146) [4]Ma Y.Research on Stress Field and Engineering Stability (12]Zhao D A,Li G L,Chen Z M,et al.Three-dimensional of the Second Mine Block in the Jinchuan Mine [Disser- FE regression analysis of multivariate geostress field of tation].Beijing:Chinese Academy of Geological Sciences, Wushaoling tunnel.Chin J Rock Mech Eng,2009, 2006:18 28(Suppl1):2687 (马宇.金川二矿区地应力场研究与地下工程稳定性学位 (赵德安，李国良，陈志敏，等，乌鞘岭隧道三维地应力场多 论文].北京：中国地质科学院，2006：18) 元有限元回归拓展分析.岩石力学与工程学报，2009,28（增 (5]Wang W.Wang L J,Qiao Z J,et al.Modeling of three f刊1)：2687)
第 期 陈章华等 基于偏最小二乘法的地应力场拟合 时, 表示采用影响权重 比较大的 向正应力 、军向 正应 力 、 向正应力 以及 切应力共 个应力分 量 若采用边 界条件 的数值不同, 仅影响如表 所 示的作为 自变量的计算值 利用偏最小二乘回归分 析所得的拟合值不会发生太大改变, 以此拟合值去 校正指定点的应力值 , 因此最终拟合的结果必然趋 近实际分布规律 由于反演初始地应力场要求能够获得地质 层 内指定点的应力, 当采用偏最小二乘回归分析时, 首先根据有限元软件得到场中此指定点的计算应力 值 , 然后采用偏最小二乘回归所得的拟合系数可以 确定指定点的应 力拟合值 , 从而确定在此点处工程 计算所需要的各种应力值 参 考 文 献 , 一 ” 盯 二, , 蔡美峰, 王双红 玲珑金矿地应力场测量结果及其分析 中国矿业 , , , , 一 此 乞几 玩 艺 不 。艺已 几 二 二, 御 二 二叽 , 马莎, 吴林峰, 赵 明献, 等 偏最小二乘回归在岩体力学 参数取值中的应用 华北水利水电学院学报, , 一 。乞二吮 乞。 , , 聂善文 基于偏最小二乘法和随机模糊理论的岩土力学 参数研究 桂林理工大学学报, , 化 ” 叱 八 几夕 几 代几夕 乞, 从 二 无 二 六。 二 、 械 。 【 」 〕 , 马宇 金川二矿区地应力场研究与地下工程稳定性 【学位 论文」北京 中国地质科学院, 、 , 、 , , , 【 〔 〕 」 一 · 云 、 乞二, , 王薇 王连捷, 乔子江, 等 三维地应力场的有限元模拟及 其在隧道设计中的应用 地球学报, , , , , 入 ,丁` · 几 几夕 刀万, , 陈香梅, 任永强, 徐智勇, 等 基于有限元数学模型的隧道 围岩初始地应力场反演 工程地球物理学报, , , , 、从 】 , , ,, ” 已。 阳七 对 亡 肠 夕, , 袁海平, 赵科, 王继伦, 等 多元线性回归地应力场反演方 法研究 有色金属科学与工程, , , 、叭 , 无 , , 张勇慧 魏倩, 盛谦, 等 大 岗山水电站地下厂房区三维地 应力场反演分析 岩土力学, 、 对 亡 七 化 爬 ” 艺 尹尹 , , 王惠文 偏最小二乘回归方法及其应用 北京 国防工业 出版社, 、认 , 肠 , , · 吼 , , , , , 一 几、 二 , , 鲁庆华, 任康乐, 周凤玺 基于偏最小二乘法实现非线性 回归分析 甘肃科技, , 一 , , , 一 倪 确 夕, , 赵德安, 李国良, 陈志敏, 等 乌鞘岭隧道三维地应力场多 元有限元回归拓展分析 岩石力学与工程学报, , 增 刊