APP下载

岩体结构面模拟的统计方法改进研究

2019-07-25斌,黄

水力发电 2019年4期
关键词:密度估计蒙特卡洛直方图

张 斌,黄 俊

(1.中国华西工程设计建设有限公司,四川成都610031;2.地质灾害防治与地质环境保护国家重点实验室(成都理工大学),四川成都610059)

0 引 言

岩体是由岩块及结构面组成的非均质不连续体,而岩体中结构面的几何特征、组合特征较为复杂,使岩体结构表现出一定的复杂性[1-3]。目前,针对岩体结构面几何特性的研究,已经取得了许多成果,但是仍然存在些许不足。在结构面三维网络模拟过程中,需对采集的数据进行合理分组,然后对每组数据建立概率模型[4]。上述过程中,对结构面的分组没有定量标准,一般根据数据整体分布趋势进行分组,主要通过极点分布图或直方图[5- 8],而在直方图统计中,带宽的确定通常也是依经验而定,对同一组数据的特征,不同带宽的直方图的反映是不同的。因此,需要一种定量确定带宽的方法,提高直方图统计的可靠性。其次,对结构面要素的分布类型的判定过于依靠经验,卡方检验[9-10]是一种很好的方法,但其操作过于复杂,极大似然估计相对简便。此外,还是有一个无法回避的问题,即结构面不能完全服从目前已知分布,甚至相差很大。在非参数估计方法中,核密度估计方法[11]在样本不能较好服从已知分布的情况下,可以模拟得到与样本较一致的随机数。

鉴于此,本文针对目前结构面模拟中存在的问题,利用现场采集的结构面数据,使用定量方法确定结构面要素的分布类型及数据带宽,通过核密度估计方法提高结构面模拟精度,并针对传统蒙特卡洛方法,提出了新的工作思路,以简化工作流程。

1 研究区概况

本文采集数据源于我国西南某在建大型水电站坝基右岸岩体(见图1)。该水电站坝址位于金沙江下游,高原深谷地貌,呈不对称“V”字形峡谷。右岸为药山山脉西坡,山峰高程在3 000 m以上,主要为陡坡与缓坡相间的地形。金沙江由南向北流入坝区,枯水期水面宽50~115 m,水位约594 m,水深9~18 m不等。拱坝部位高程825 m处谷宽590~713 m。研究区域主要出露二叠系上统峨眉山组柱状节理玄武岩,采用测线法对坝基岩体进行结构面精测,得到结构面极点分布,见图2。

图1 右岸建基面岩体

图2 结构面极点分布

2 结构面模拟的蒙特卡洛方法

2.1 方法介绍

蒙特卡洛方法(Monte Carlo method)可以用来解决具有概率解释的任何问题,因其求解的问题本身具有内在的随机性,因此必须借助计算机的运算能力才能模拟这种随机的过程。在研究具有概率解释的问题时,对随机数抽样,使抽样结果符合问题的分布形式,从而对研究的问题进行模拟。由于需要大量随机数,蒙特卡洛方法必须依赖于计算机。尽管岩体结构面是随机分布的,但经验表明,结构面的空间分布在一定条件下依然符合一定的概率分布。因此,对精测数据使用蒙特卡洛方法,便可得到一定深度范围内的结构面特征数据。

2.2 工作流程

使用蒙特卡洛方法对结构面进行三维模拟,首先要判定数据所服从的分布类型,然后生成服从该分布的随机数[12]。由于岩体受到不同时期不同类型构造的影响,整体不易服从某一确定的分布,所以根据直方图起伏和集中趋势对数据进行分组,每组数据单独使用蒙特卡洛方法模拟。需要注意的是,每组数据在使用蒙特卡洛方法前后,其在数据总量中所占的比值应保持一致。

2.3 直方图统计的带宽优化

带宽在直方图统计中是一个自由参数,一般根据经验确定。直方图统计是一种对未知量的非参数估计。直方图中,带宽过大或过小,得到的直方图都不能反映问题的真实情况。带宽过大,造成直方图过于平滑,易忽略细节;带宽过小,放大了细节在直方图中的作用,得到的直方图会产生较多的尖刺,增加了异常的比重,也不能较好地反映真实情况。因此,带宽是否取值得当直接决定了直方图统计的质量,影响后续论述。为量化带宽的选取是否得当,需要一个误差函数来衡量,本文使用平均积分平方误差(MISE)作为带宽选取的误差,平均积分平方误差MISE表达式为[12-13]

(1)

式中,f(x)为未知的真实的密度分布;fh(x)为基于样品的以带宽h为独立变量的假定分布;E是样本的期望。

由于密度函数f(x)未知,无法直接得到方程的解,可通过计算机的多次计算,寻找该式的最小值,从而代替直接求解方程[14]。本文使用上述方法编写的程序对带宽进行优化,倾角优化后的最佳带宽为1.72,倾向为7.06,迹长为0.060 232 7,相对应的最佳区间数分别为52、51和50。

2.4 结构面要素的分布拟合优度检验

使用带宽优化之后的数据进行直方图统计(见图3),结果与极点分布图保持一致,能够均衡地反映结构面的优势产状和发育较弱的产状。经过带宽优化的直方图,降低了主观因素对直方图趋势的影响。依据统计结果对数据进行分组,根据数据的峰值和起伏状况,将倾角分为[0°, 50°)和[50°, 90°]2个组,倾向分为[0°, 180°)、[180°, 280°)和[280°, 360°]3个组,迹长不分组。

图3 结构面要素统计直方图

分布倾向0°~180°AD值P值极大似然比P倾向180°~280°AD值P值极大似然比P倾向280°~360°AD值P值极大似然比P正态5.634<0.005—2.238<0.005—0.7630.046—Box-Cox变换0.8490.028—1.58<0.005—0.3030.568—对数正态18.57<0.005—1.58<0.005—0.5510.152—三参数对数正态5.671∗01.559∗0.9030.34∗0.031指数42.048<0.003—109.343<0.003—39.783<0.003—双参数指数41.772<0.0100.04740.2<0.01005.999<0.0100Johnson变换——————0.2990.578—三参数对数Logistic2.396∗01.003∗0.528———

对结构面数据进行分组后,分别对每组数据进行参数估计,使用Minitab进行分布类型检验。本文使用Anderson-Darling检验(简称“AD检验”)进行极大似然估计(当能够进行似然比检验时显示极大似然比P,P值是一个概率,用来度量否定原假设的证据)。Anderson-Darling统计量(即AD值,表示假定分布与样本的差值的平方)用来衡量数据遵循特定分布的程度。对指定的样本和分布,分布越适合样本,AD值越小,其表达式为[15]

(2)

式中,n是样本数;f(x)是待检验的假定分布;fh(x)为样本的累积分布函数。本文使用相应的P值(当存在时)来检验数据是否服从假定分布。如果P值小于选定的阈值α(本文取0.05,但0.05只是一个常用的阈值,P值小于0.05并不是无效假设,只能说明其拒绝原假设的理由更充分),则拒绝数据来自该分布的假设。P值在数学上不存在时,Minitab软件就不显示AD检验的P值。对研究区结构面要素进行分布拟合优度检验,得到结构面要素的概率分布图(默认置信区间为95%)。对6组数据进行分布类型检验后的结果见表1~3。表中, *表示无法计算的缺失值。

表1 倾角分布类型检验结果

表3 迹长分布类型检验结果

综合AD值小、P值大的原则,确定原数据分布类型。其中,倾向0°~180°区间较符合Box-Cox变换后正态分布,倾向180°~280°区间较符合三参数对数正态分布,倾向280°~360°区间较符合Johnson 变换后正态分布;倾角0°~50°区间较符合Box-Cox变换后正态分布,倾角50°~90°区间较符合Box-Cox变换后正态分布;迹长服从Johnson变换后正态分布。各自对应的在95%置信区间的累积概率分布见图4~6。

图4 倾角累积概率分布

图5 倾向累积概率分布

图6 迹长累积概率分布

分布拟合优度检验结果表明,仅迹长在未经分组的情况下,对假定分布类型的接受程度较高。岩体作为非均质体存在,在各种地质营力的作用下,被结构面切割,在一定区域内形成几何相似的块体,因此在统计量没有大到足以代表整个岩体结构面时,精测数据应该是非连续数据。在特定地质营力作用下,形成的数量相对较多、且易于被观察和测量的部分结构面,即为优势结构面。优势结构面数量较大,易于被观测,因此优势结构面数据在精测数据占比较大。林德伯格-费勒中心极限定理亦表明,结构面迹长作为描述结构面出露程度的标量,其原始数据(未分组)应该服从正态分布。

2.5 模拟结果

本次模拟在Jupyter notebook中完成。Jupyterno-tebook是基于Python的交互式解释器。整个程序设计过程为:首先,程序初始化后导入数据,设置并初始化变量,对数据分组;然后,确定模拟的各组数据的总量,根据各组数据对应的变换方式分别对每组数据进行变换,使用NumPy工具库生成伪随机数,接着对生成的伪随机数反变换,得到模拟结果,见图7。

对比图3,使用蒙特卡洛方法模拟的结果,在集散趋势上与原数据保持一致。经过带宽优化后,为避免主观因素的影响,使用量化方法判断精测数据的分布类型。通过分布类型检验,获得了数据的一般分布类型,但由于细节的损失,模拟结果无法较好再现原数据的分布特征。

图7 模拟结果直方图

图8 原始数据与核密度估计结果对比

3 结构面模拟的核密度估计

使用蒙特卡洛方法进行结构面模拟时,需要对数据进行分布类型检验,而模拟结果的好坏取决于数据对假定分布类型的服从程度。对数据切割分组的不连续,可能降低模拟结果的准确性。为提高模拟结果与原始数据的相似度,使用核密度估计方法对结构面数据进行模拟。

3.1 核密度估计

核密度估计是估计随机变量的概率密度函数的非参数方法,这意味对随机变量的估计不再需要固定的参数,而且程序设计可以减少很多工作。在直方图统计中,对于区间(xk-1,xk]((xk-1,x]为其中某1个条带)中的所有样本都同等对待。如在[0,1]区间,有2个样本取值0+和1-,前者表达的信息是分布在0周围有密度,而后者应该在1周围有密度,两者存在不同的含义。因此,对于直方图统计,给定宽度h,样本Xi应该反应区间[Xi-h/2,Xi+h/2](该区间的分布函数为I(-h/2≤x≤h/2))的信息,定义[16]

Kh(x)=I(-h/2≤x≤h/2)/h

(3)

对应的估计fh(x)为

(4)

再进一步把宽度h看成参数,直接定义K(x)=I(|x|≤1)/2,若该区间的权重相等,则估计fh(x)为

(5)

(6)

Kh(· )=K(·/h)/h称为核函数,理论上任何非负单峰的概率密度函数都可以作为核函数。不同的核函数在对目标进行核密度估计时会产生不同估计模型,在每个估计区间体现核函数原本的性质。目前,常用的核函数是高斯(Gaussian)核函数,表示为K(x;h)∝exp[-x2/2h2]。

3.2 结构面模拟的核密度估计

使用核密度估计方法,依然要先对数据进行带宽优化。使用优化后的带宽,在Jupyter notebook中,利用Scikit-learn工具计算。为提高采样精度,Scikit-learn中的KDE模块引入了最近邻算法中的KDTree和BallTree,使采样模型更贴近原始数据。使用Gaussian内核的核密度估计方法对精测数据模拟。原始数据与核密度估计结果对比见图8。图8中,直方图为原始数据统计结果,曲线为核密度估计。核密度估计模拟结果见图9。

图9 核密度估计模拟结果

3.3 结果对比

从图8可知,高斯内核的核密度估计曲线较为平滑,其概率密度的变化均为渐变曲线。而直方图统计中,相邻区首尾相接处的数据是跃变的,而实际上在区间连接处数据应该是连续的。因此,使用核密度估计方法对数据进行统计,能够提高模拟的精度。核密度估计曲线与直方图统计数据增减趋势一致,但不完全贴近方块,为模拟结果预留了变化区间,表明核密度估计的结果受偶然因素的影响较小,可以减少过拟合的发生。在保证带宽为最佳的前提下,使用核密度估计既能够保留细部特征,又能够兼顾整体分布趋势。

从图9可知,对比原始数据和蒙特卡洛方法的结果,核密度估计结果更符合原始数据的集散趋势,并且对细部的反映更真实。观察蒙特卡洛模拟结果,其数据整体较为规律,数据的整体分布趋势符合公式生成数据的特点。而核密度估计方法模拟的结果更加贴近原数据方块。在整体趋势上,两者的结果较为相似,但在局部数据的呈现上,核密度估计更加符合原数据的特点,对细节的把握更充分。

在工作流程上,蒙特卡洛方法的工作量更大,实施步骤更多。核密度估计只需要进行带宽优化、直方图统计、核密度估计、随机数生成,最后成果检验,相对而言,工作量大大减少。从量化的角度看,蒙特卡洛方法在分组时存在主观判断,可能对模拟的成果有一定影响。因此,在结构面三维模拟中,使用核密度估计方法,有利于模拟精度的提高,同时减少工作量。

4 结 语

本文以西南某大型在建水电站建基面结构面精测数据为基础,对结构面要素进行统计分析,使用平均积分平方误差优化直方图统计带宽。引入核密度估计方法,对结构面模拟的统计方法进行改进,对比蒙特卡洛方法的流程和结果,得出以下结论:

(1)使用平均积分平方误差对直方图带宽进行优化后,直方图统计结果与极点分布图保持一致,能够照顾到结构面的优势产状和发育较弱的产状。经过带宽优化的直方图,能避免主观因素对直方图的影响,为结构面分布类型的确定和核密度估计做好铺垫。

(2)在蒙特卡洛模拟中,采用了极大似然估计,根据各组数据对应的AD值和P值,判定各组数据对假定分布类型的服从程度,从而确定分布类型。蒙特卡洛方法的结果表明,其把握数据的整体分布趋势较好,但细节丢失严重。

(3)在探索提升结构面模拟精度的过程中,引入核密度估计方法,结果表明,其结果能够较好地还原数据的整体分布特征。此方法省略了分布类型检验等步骤,忽略了主观因素对结果的影响,提高了模拟精度,减少了工作量。

猜你喜欢

密度估计蒙特卡洛直方图
面向鱼眼图像的人群密度估计
符合差分隐私的流数据统计直方图发布
基于MATLAB 的核密度估计研究
一种基于改进Unet的虾苗密度估计方法
基于自适应带宽核密度估计的载荷外推方法研究
基于FPGA的直方图均衡图像增强算法设计及实现
征服蒙特卡洛赛道
用直方图控制画面影调
基于蒙特卡洛法的车用蓄电池20h率实际容量测量不确定度评定
中考频数分布直方图题型展示