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所示.通 在求出偏最小二乘拟合值的同时，可以拟合出
第 期 陈章华等 基于偏最小二乘法的地应力场拟合 表 井各测点计算应力值 过程序计算 , 最后得到的应力拟合值如表 所示 认 、 三 向应力 向应力 向应力 一 一 一 一 一 一 ` 一 一 一 一 向切应力 一 一 一 了 一 一 一 一 一 月 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 开始 输人 自变量 与因变量 矩 阵 标准化操作 表 井各测点实测应力值 入 、 水平面内最大 水平面 内最小 自重应力 主应力 主应力 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 召 一 一 求残差矩阵 并 求协方差矩阵 代替原矩阵 得第一主成分 应用偏最小二乘回归得到拟合值 首先对 自变量 和 因变量 进行标准化 处理 , 以方便它们在 同一量纲上进行运算操作, 得 到 和 通过协方差矩阵的计算 , 求出最大特征值 及其所对应 的单位特征向量, 继而求 出主成分 根据最小二乘迭代和普遍适用且简便 的有 效性判别法 特 征根大于 法 , 筛选 出所有主成分 主成分 是根据 和 叨 的矩阵相乘而来 , 其 中 叨 是指特征向量组成的矩阵 之后, 利用最小二乘回归原理, 得到主成分 与 因变量的关系 此时, 可求出因变量的拟合值 以 及因变量与 自变量的关系式 , 最后将标准化之后的 因变量 的值, 进行标准化反 向计算, 得到实际的因 变量拟合值 由于偏 最小二乘算法 中涉及 的原始数据和算 法的使用要求较高的精度 , 因此使用计算机编程 的 算法可 以又准 又快地得出拟合结果 考虑到较多的 涉及矩阵的运算 , 运用 ' , 编程计算可 以快 速方便地得出结果 主要的程序流程如图 所示 通 得所有主成分 主成分与因变量关系 自变量 与因变量关 系 因变量拟合值 图 主要程序流程图 表 偏最小二乘拟合值 妞 最大主应力 实测值 拟合值 一 一 一 一 一 一 名 一 一 一 一 一 一 一 一 一 习 一 一 一 一 最小主应力 实测值 拟合值 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 ` 一 一 习 一 一 , 自重应力 实测值 拟合值 一 一 一 一 一 一 , 一 一 一 一 一 一 一 一 一 一 一 一 乃 在求出偏最小二乘拟合值的同时, 可以拟合出