基于GMM与EM的堆积层滑坡失稳模式概率判识方法及系统

基于gmm与em的堆积层滑坡失稳模式概率判识方法及系统
技术领域
1.本发明涉及堆积层滑坡失稳模式概率判识技术领域，特别涉及一种基于gmm与em的堆积层滑坡失稳模式概率判识方法及系统。


背景技术：

2.在滑坡失稳模式判识的现有技术方案中，随机有限元法得到了广泛的应用。随机有限元法，也称概率有限元法，是随机场理论与有限元方法相结合的产物，是在传统的有限元方法的基础上发展起来的随机数值分析方法。其中，非侵入式有限元法由于其易操作，流程相对简单且易于实现的特点，广受青睐；通过非侵入式随机有限元法对滑坡进行失稳模式判别的一般步骤为：1.建立滑坡确定有限元模型；2.建立滑坡土体参数的随机场；3.建立滑坡不确定模型；4.稳定性分析搜索最危险滑面；5.基于最危险滑面间相关系数，判别失稳模式。具体实现步骤为：
3.步骤1：通过现有有限元软件建立滑坡确定有限元模型，具体包括：
①
确定滑坡概化模型；
②
赋予材料属性；
③
设置求解分析步；
④
施加荷载和设置边界条件；
⑤
划分网格；
4.步骤2：建立滑坡土体参数随机场，主要包括：
①
选择所研究的土体参数作为随机变量，拟定其空间变异信息统计特征，包括均值μ、方差σ、变异系数cov、相关系数ρ、波动范围δ；
②
拟定随机场的离散方法，目前主要的随机场离散方法有：中心点离散法、局部平均法、随机场插值、随机场加权积分、随机场正交展开等。中心离散：是随机场最简单的一种离散方法。该方法用随机场在每个单元中心点的值来表征该随机场在每个单元的属性，因而随机场在每个单元内部都是常量，且等于它在各个单元中心的值。局部平均：用随机场在各离散单元内的局部平均值来表征单元体的特征。随机场的插值：将随机场在单元内的值用单元结点处值的插值函数来表示于是随机场的统计特性可由各单元结点处随机变量间的统计特性近似反映。随机场的加权积分法：在单元刚度矩阵的推导过程中采用随机场在单元高斯点上的加权积分以表征单元上的随机场。随机场的正交展开：将材料特性参数随机场进行一正交展开并由此推导出刚度矩阵的级数展开式，从而获得位移、应力的统计特性。
5.步骤3：建立滑坡不确定有限元模型：以步骤1建立的确定有限元模型网格数据为依据，将步骤2中获得的土体参数随机场数据，对应赋予给每一个模型，以此建立具有非均质材料属性的有限元模型，即不确定有限元模型；n个随机场，将得到n个不确定有限元模型。
6.步骤4：稳定性分析搜索最危险滑面：目前滑坡稳定性分析主要以传统的极限平衡法为主，包括bishop法、morgenster-price法、spencer法等，以及以强度折减法为代表的数值方法。当前，业内对于“滑坡失稳”状态的判断，主要以稳定系数fs的值来判断滑坡是否失稳或者破坏：fs＜1，滑坡失稳；fs＝1，滑坡处于极限平衡状态；fs＞1，滑坡处于稳定状态。每进行一次稳定性分析获得一个fs，且对应获得一个最危险滑面。
7.步骤5：基于最危险滑面间相关性，判别失稳模式：在岩土工程可靠度理论中，通常将滑坡看作是由多个失稳模式组成的联合系统，每一个最危险滑动面都看作一种失稳模
式，同时每一种失稳模式间具有明显的相关性，因此基于该特性研究不同滑面间安全系数的相关性来判别失稳模式，其判据通常以相关系数阈值来确定。
8.现有的滑坡失稳模式判识方法中，通常都以理想化(几何外形规整，利于计算)的滑坡模型进行分析计算，然而实际的工程研究中，滑坡的形态特征十分复杂，且滑坡类型诸多，传统的判识方法无法普遍适用于具体类型的滑坡，因此缺乏一定的针对性。在失稳模式的判别中，以不同滑面间安全系数的相关性作为依据，通过一个给定的相关系数阈值来判别失稳模式，然而，安全系数以1作为分界，当安全系数在1附近波动时，边坡的稳定状态并无根本性差异，但状态划分结果却完全相反。这种不严谨的方式缺乏一定的客观认识，同时也会影响到失稳模式判别的准确性。


技术实现要素：

9.为了解决上述现有技术方案客观存在的缺陷，本发明提供一种基于gmm与em的堆积层滑坡失稳模式概率判识方法及系统，以针对自然界大量存在的堆积层滑坡，在一定程度上弥补现有失稳模式判识方法，无法普遍适用于自然界某一种客观存在的滑坡类型的缺点；再者，引入机器学习的方法，以不确定有限元模型稳定性分析结果中的节点位移和变形体体量为观测数据，借助混合高斯模型(gmm)的期望最大化算法(em)，对堆积层滑坡进行失稳模式概率判识，此方法不单纯基于滑面的安全系数进行分析，因此现有技术的客观缺点在本发明中对失稳模式的判别准确度影响非常小。
10.根据本发明的其中一方面，本发明提供额了一种基于gmm与em的堆积层滑坡失稳模式概率判识方法，具体包括如下步骤：
11.s1、获取堆积层滑坡的计算参数，所述计算参数包括：堆积层滑坡模型的几何参数、堆积体的物理力学参数以及堆积体强度参数空间变异信息统计；
12.s2、利用s1中的计算参数，建立堆积层滑坡地质概化模型，以堆积层滑坡地质概化模型为基础，在有限元软件中建立堆积层滑坡确定有限元模型；
13.s3、对有限元软件进行二次开发，以提取步骤s2中堆积层滑坡确定有限元模型的各个单元的中心坐标；
14.s4、利用所述各个有限元模型的中心坐标，将黏聚力c和内摩擦角φ作为随机变量，进行相关随机场模拟，得到多个随机场；
15.s5、利用所述多个随机场，采用非侵入式随机有限元方法批量建立堆积层滑坡不确定有限元模型，并进行稳定性分析，获得每一个堆积层滑坡不确定有限元模型的节点位移；
16.s6、以节点位移作为观测值，通过k-means聚类分析，将滑坡模型整体聚为稳定体和变形体两类；每个观测值包含一个堆积层不确定有限元模型的所有节点位移；
17.s7：基于步骤s6得到的聚类结果，统计变形体体量，并基于统计结果，以变形体体量为观测对象，通过混合高斯模型的期望最大化算法进行失稳模式判识。
18.进一步地，在本发明的基于gmm与em的堆积层滑坡失稳模式概率判识方法中，堆积层滑坡模型的几何参数包括：边坡几何形状控制点的坐标；堆积体的物理力学参数包括：稳定基岩层的黏聚力cr、内摩擦角φr、杨氏模量er、泊松比νr和天然重度γr；堆积体强度参数空间变异信息统计包括：黏聚力的均值μc、标准差σc、变异系数covc，内摩擦角的均值μ
φ
、标
准差σ
φ
、变异系数cov
φ
，垂直波动范围δv、水平波动范围δh，黏聚力c和内摩擦角φ的相关系数ρ
c,φ
。
19.进一步地，在本发明的基于gmm与em的堆积层滑坡失稳模式概率判识方法中，步骤s3包括：提取有限元软件建模过程中生成的包含所有建模过程的命令文件，修改此命令文件，并将修改后的命令文件在有限元软件中运行，获得各个单元的中心坐标。
20.进一步地，在本发明的基于gmm与em的堆积层滑坡失稳模式概率判识方法中，步骤s4中，所述进行相关随机场模拟具体是指：假设堆积体关键强度参数粘聚力和内摩擦角均服从对数正态分布，基于随机场理论和所述中心坐标，使用基于乔列斯基分解的中心点法并采用拉丁超立方抽样进行随机场模拟；随机场模拟时，各个单元的杨氏模量、泊松比和天然重度相同。
21.进一步地，在本发明的基于gmm与em的堆积层滑坡失稳模式概率判识方法中，步骤s5具体包括：
22.s51、对所述随机场进行堆积体的物理力学参数的赋值；
23.s52、生成堆积层滑坡不确定有限元模型批计算脚本；
24.s53、调用有限元软件内核运行步骤s52中计算脚本，进行堆积层滑坡不确定有限元模型的稳定性分析批处理；
25.s54、批量提取稳定性分析批处理的处理结果中每一个堆积层滑坡不确定有限元模型的节点位移。
26.进一步地，在本发明的基于gmm与em的堆积层滑坡失稳模式概率判识方法中，步骤s7具体包括如下步骤：
27.s71：统计变形体体量x＝(x1,x2,
…
,xn)，n为观测数据总量，i＝1,2,3
…
n；xi为第i个观测数据点；设观测数据点x所属的失稳模式类别为z，其取值有m种(z
(1)
,z
(2)
,
…
,z
(m)
)；
28.s72：确定高斯混合模型的概率密度函数：
[0029][0030]
s73：初始化对数化后的均值方差权重系数λj；
[0031]
s74：使用初始化的均值方差权重系数λj，计算qi(z
(j)
)：
[0032][0033]
s75：使用上一步计算得到的qi(z
(j)
)值重新估计参数均值方差和权重系数λj：
[0034]
[0035][0036][0037]
s76：计算最大似然函数值：
[0038][0039]
s77：当步骤s74中3个参数均值方差及权重系数λj和步骤s76中计算得到的最大似然函数值趋于稳定时，迭代计算结束，否则利用步骤s74至s76进行迭代，直到计算收敛，得到未知参数λj,的最终解；
[0040]
其中，部分参数的含义如下:
[0041]
λj：第j种失稳模式所占权重；
[0042]
第j种失稳模式变形体体量经对数化后的标准差；
[0043]
第j种失稳模式变形体体量经对数化后的期望；
[0044]
xi：样本数据集合中第i个观测数据点；
[0045]
qi(z
(j)
)：xi对z的分布函数；
[0046]
nj：表示被分到第j类的观测数据点的有效数量。
[0047]
根据本发明的另一方面，本发明还提供了一种基于gmm与em的堆积层滑坡失稳模式概率判识系统，具备计算机存储介质，所述计算机存储介质存储有计算机程序，该计算机程序被处理器执行时实现上述任一项所述的基于gmm与em的堆积层滑坡失稳模式概率判识方法。
[0048]
与现有技术方案相比本发明的优势在于：
[0049]
现有滑坡失稳模式识别技术方案中，通常以理想化概化边坡作为技术方法的佐证算例，但在实际工程中，这种理想化的概化边坡模型并不存在，因此现有的大多数技术方案缺乏一定的普遍适用性，本发明提供一种基于gmm与em的堆积层滑坡失稳模式概率判识方法及系统，针对性的以自然界大量存在的堆积层滑坡为处理对象，在一定程度上弥补现有失稳模式判识方法，无法普遍适用于自然界某一种客观存在的滑坡类型的缺点；再者，引入机器学习的方法，以不确定有限元模型稳定性分析结果中的节点位移和变形体体量为观测数据，借助混合高斯模型(gmm)的期望最大化算法(em)，对堆积层滑坡进行失稳模式概率判识，获得每一种失稳模式的权重、期望及标准差。
附图说明
[0050]
图1为基于gmm与em的堆积层滑坡失稳模式概率判识方法一实施例的流程图。
[0051]
图2为堆积层滑坡地质概化模型示意图；
[0052]
图3为堆积层滑坡不确定有限元模型示意图；
[0053]
图4为堆积体关键强度参数相关随机场一次实现的图例；
[0054]
图5为变形体体量统计直方图；
[0055]
图6为实施算例计算所得失稳模式第一示意图；
[0056]
图7为实施算例计算所得失稳模式第二示意图。
具体实施方式
[0057]
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地描述。
[0058]
本发明提供了一种基于gmm与em的堆积层滑坡失稳模式概率判识方法及系统，本发明以堆积层滑坡为处理对象，将堆积体强度参数黏聚力和内摩擦角作为随机变量建立相关随机场，采用非侵入式随机有限元方法批量建立堆积层不确定有限元模型，进行堆积层不确定有限元模型的稳定性分析，获得每一个堆积层不确定有限元模型的节点位移；以节点位移为观测值，通过k-means聚类分析，将滑坡整体聚为稳定体与变形体两类；统计变形体体量，基于变形体体量统计结果，通过混合高斯模型(gmm)的期望最大化算法(em)，最后分离出各类失稳模式，获得每一类失稳模式的统计特征，从而进行堆积层滑坡失稳模式的概率判识。
[0059]
参考图1，图1为基于gmm与em的堆积层滑坡失稳模式概率判识方法一实施例的流程图。本实施的基于gmm与em的堆积层滑坡失稳模式概率判识方法具体包括如下步骤：
[0060]
s1、获取堆积层滑坡的计算参数，所述计算参数包括：堆积层滑坡模型的几何参数、堆积体的物理力学参数以及堆积体强度参数空间变异信息统计。
[0061]
堆积层滑坡模型的几何参数包括：边坡几何形状控制点的坐标；堆积体的物理力学参数包括：稳定基岩层的黏聚力cr、内摩擦角φr、杨氏模量er、泊松比νr和天然重度γr；堆积体强度参数空间变异信息统计包括：黏聚力的均值μc、标准差σc、变异系数covc，内摩擦角的均值μ
φ
、标准差σ
φ
、变异系数cov
φ
，垂直波动范围δv、水平波动范围δh，黏聚力c和内摩擦角φ的相关系数ρ
c,φ
。
[0062]
该实施算例中，确定该模型的具体计算参数包括：
[0063]
堆积层滑坡模型的几何参数，包括：模型的高度h＝400m,宽度l＝933m,边坡几何形状控制点坐标z1(933,0)，z2(933,50)，z3(902.2,50)，z4(875,75)，z5(850,98)，z6(735,149.7)，z7(668,170.8)，z8(658.6,172.8)，z9(605.2,178.7)，z10(588.4,183.7)，z11(572.3,191.7)，z12(523,214)，z13(325,288.5)，z14(48,397)，z15(41,400)，z16(38.5,400)，z17(0.,400)；堆积体与下层稳定基岩分界面控制点：y1(570,113)，y2(310,250.7)，y3(45.3,391)，y4(721.4,85.5)
[0064]
堆积体的物理力学参数：包括：堆积体粘聚力cs＝15kpa、内摩擦角φs＝28.5
°
、天然重度γs＝19.82kn/m3、杨氏模量es＝20.8mpa、泊松比νs＝0.35；稳定基岩层的黏聚力cr＝4800kpa、内摩擦角φr＝0.25
°
、杨氏模量er＝1200mpa、泊松比νr＝0.32和天然重度γr＝20.45kn/m3。
[0065]
堆积体关键强度参数空间变异信息统计：包括：堆积体黏聚力的均值μc＝15kpa、变异系数covc＝0.3，内摩擦角的均值μ
φ
＝28.5
°
、变异系数cov
φ
＝0.2，垂直波动范围δv＝40m、水平波动范围δh＝4.0m，黏聚力c和内摩擦角φ的相关系数ρ
c,φ
＝0。
[0066]
s2、利用s1中的计算参数，建立堆积层滑坡地质概化模型，如附图2所示。以堆积层
滑坡地质概化模型为基础，在有限元软件中建立堆积层滑坡确定有限元模型。
[0067]
步骤s2具体包括：
[0068]
根据堆积层滑坡模型的几何参数建立地质概化模型，并进一步对概化模型进行实体装配、赋予材料(物理力学参数赋值)、设置求解分析步、划分网格、建立任务，完成有限元模型构建，该些步骤均为现有技术，具体可参考相关的文件。
[0069]
s3、对有限元软件进行二次开发，以提取步骤s2中堆积层滑坡确定有限元模型的各个单元的中心坐标。
[0070]
具体为：提取有限元软件建模(堆积层滑坡确定有限元模型)过程中生成的包含所有建模过程的命令文件，修改此命令文件，将第三方库函数xlwt嵌入命令行中，将修改后的命令文件在有限元软件中运行，获得各个单元的中心坐标。
[0071]
s4、利用所述各个有限元模型的中心坐标，将黏聚力c和内摩擦角φ作为随机变量，进行相关随机场模拟，得到多个随机场。
[0072]
进行相关随机场模拟具体是指：假设堆积体关键强度参数粘聚力和内摩擦角均服从对数正态分布，基于随机场理论和所述中心坐标，使用基于乔列斯基分解的中心点法并采用拉丁超立方抽样进行随机场模拟。随机场模拟时，建立这个随机场需要堆积体关键强度参数粘聚力和内摩擦角这个两个参数各自的均值、变异系数、波动范围(水平和垂直)和两者的相关系数，以及杨氏模量、泊松比和天然重度，其中各个单元的杨氏模量、泊松比和天然重度相同。如附图4所示，为粘聚力c和内摩擦角φ相关标准随机场的一次典型实现。
[0073]
s5、利用所述多个随机场，采用非侵入式随机有限元方法批量建立堆积层不确定有限元模型，并进行稳定性分析，获得每一个堆积层不确定有限元模型的节点位移。步骤s5具体包括：
[0074]
s51、对所述随机场进行堆积体的物理力学参数的赋值；
[0075]
s52、生成堆积层滑坡不确定有限元模型批计算脚本；
[0076]
s53、调用有限元软件内核运行步骤s52中计算脚本，进行堆积层滑坡不确定有限元模型的稳定性分析批处理；
[0077]
s54、批量提取稳定性分析批处理的处理结果中每一个堆积层滑坡不确定有限元模型的节点位移。
[0078]
在粘聚力c和内摩擦角φ相关随机场模拟结果的基础上，再次通过修改有限元软件建模过程中生成的包含所有建模过程的脚本文件，脚本文件修改完成后，调用有限元软件内核运行该脚本，批量建立不确定有限元模型并进行稳定性分析，图3为一个积层不确定有限元模型。计算结束后，自动保存结果文件，通过结果文件提取所有模型节点位移并保存为nodalu.txt。
[0079]
s6、以节点位移作为观测值，通过k-means聚类分析，将滑坡模型整体聚为稳定体和变形体两类；每个观测值包含一个堆积层不确定有限元模型的所有节点位移。将节点位移作为观测值k＝(k1,k2,...,kn)，每个观测值x是m维实向量，表示具体节点的观测值，将n个观测值划分为2个集合(分别稳定体和变形体)，每个观测值属于具有最近均值的那一类，n为节点数。
[0080]
s7：基于步骤s6得到的聚类结果，统计变形体体量，并基于统计结果，以变形体体量为观测对象，通过混合高斯模型(gmm)的期望最大化算法(em)进行失稳模式判识。变形体
体量可以是变形区域的面积或者变形网格的总面积。每一个堆积层滑坡不确定有限元模型经计算后得到一个变形体体量，一个随机场产生一个堆积层滑坡不确定有限元模型，几个随机场就有几个堆积层滑坡不确定有限元模型。k-means聚类分析结束后得到每一个堆积层滑坡不确定有限元模型的变形体体量保存为mass.txt，统计所有堆积层滑坡不确定有限元模型的变形体体量保存为massall.txt用于失稳模式判识。
[0081]
通过步骤s6，基于不确定有限元模型稳定性分析后获得的节点位移，进行k-means聚类将不确定有限元模型整体分为变形体和稳定体，假设变形体体量样本(massall.txt)x＝(x1,x2,...,xn)，n为堆积层滑坡不确定有限元模型的个数。
[0082]
首先，对所求未知参数λj,的计算公式进行简单推导：
[0083]
先给出单一对数正态分布的概率密度函数：
[0084][0085]
式(1)中：x:变形体体量数据集合；
[0086]
x:变形体体量；
[0087]
μ
lnx
:x经对数化后的均值；
[0088]
σ
lnx
:x经对数化后的标准差；
[0089]
变形体体量的分布假定为多个对数正态分布的混合，设观测数据点x所属的失稳模式类别为z，其取值有m种(z
(1)
,z
(2)
,...,z
(m)
)，则混合模型的联合概率密度函数为：
[0090][0091]
其中xi对z的分布函数用qi(z)表示：
[0092][0093]
式中：z：失稳模式类别：
[0094]
m:失稳模式类别总数：
[0095]
λj：第j种失稳模式所占权重：
[0096]
第j种失稳模式变形体体量经对数化后的标准差；
[0097]
第j种失稳模式变形体体量经对数化后的期望；
[0098]
xi：样本数据集合x中第i个数据点。
[0099]
写出式(2)的对数似然函数为：
[0100][0101]
将的导数相对于混合分布中某类分布的均值设为零，得
到：
[0102][0103]
将式(5)两侧同时乘以结合式(3)化简得：
[0104][0105][0106]
式(7)中nj表示被分到第j类的样本数据点的有效数量。
[0107]
与式(5)同理，将的导数相对于混合分布中某类分布的方差设为零，得到：
[0108][0109]
利用拉格朗日乘法对权重系数λj最大化
[0110][0111]
将式(9)的导数相对于权重系数λi设为零，可得：
[0112][0113]
结合式(3)将式(10)两侧同时乘以λj得到η＝-n，因此式(10)可变换为：
[0114][0115]
完成上述公式推导后，可以得出步骤s7具体的实现步骤为：
[0116]
步骤1：统计变形体体量(massall.txt)x＝(x1,x2,
…
,xn)，绘制频率直方图如附图5所示。其中，n为观测数据总量(即堆积层滑坡不确定有限元模型个数)，i＝1,2,3
…
n；xi为第i个观测数据点；设观测数据点x(xi中的任意一个)所属的失稳模式类别为z，其取值有m种(z
(1)
,z
(2)
,
…
,z
(m)
)。
[0117]
步骤2：写出高斯混合模型(gmm)的概率密度函数：
[0118][0119]
下面执行期望最大化算法(em):
[0120]
步骤3：初始化对数化后的均值方差权重系数λj；
[0121]
步骤4：使用初始化的均值方差权重系数λj，计算qi(z
(j)
)：
[0122][0123]
步骤5：使用上一步计算得到的qi(z
(j)
)值重新估计参数均值方差和权重系数λj：
[0124][0125][0126][0127]
步骤6：计算式(11)的最大似然函数值：
[0128][0129]
步骤7：当步骤4中3个未知参数λj,和步骤5中计算得到的最大似然函数值趋于稳定时，迭代计算结束，否则重复步骤4,5,6直到计算收敛，得到未知参数λj,的最终解。
[0130]
公式中部分参数的含义:
[0131]
m:失稳模式类别总数；
[0132]
λj：第j种失稳模式所占权重；
[0133]
第j种失稳模式变形体体量经对数化后的标准差；
[0134]
第j种失稳模式变形体体量经对数化后的期望；
[0135]
xi：样本数据集合中第i个观测数据点；
[0136]
qi(z
(j)
)：xi对z的分布函数；
[0137]
nj：表示被分到第j类的观测数据点的有效数量。
[0138]
通过上述所有步骤，实施算例被划分为两类失稳模式即m＝2，λj,经式(18)，(19)变换后得到实施算例的两类失稳模式的λ,μ,σ统计于下表；并将两类失稳模式绘制于附图6和附图7。
[0139]
[0140][0141][0142]
根据本发明的另一方面，本发明还提供了一种基于gmm与em的堆积层滑坡失稳模式概率判识系统，具备计算机存储介质，所述计算机存储介质存储有计算机程序，该计算机程序被处理器执行时实现上述任一项所述的基于gmm与em的堆积层滑坡失稳模式概率判识方法。
[0143]
在不冲突的情况下，本文中上述实施例及实施例中的特征可以相互结合。
[0144]
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。