張 斌，黃 俊

(1.中國華西工程設(shè)計(jì)建設(shè)有限公司，四川成都610031；2.地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室(成都理工大學(xué))，四川成都610059)

0 引 言

巖體是由巖塊及結(jié)構(gòu)面組成的非均質(zhì)不連續(xù)體，而巖體中結(jié)構(gòu)面的幾何特征、組合特征較為復(fù)雜，使巖體結(jié)構(gòu)表現(xiàn)出一定的復(fù)雜性[1-3]。目前，針對(duì)巖體結(jié)構(gòu)面幾何特性的研究，已經(jīng)取得了許多成果，但是仍然存在些許不足。在結(jié)構(gòu)面三維網(wǎng)絡(luò)模擬過程中，需對(duì)采集的數(shù)據(jù)進(jìn)行合理分組，然后對(duì)每組數(shù)據(jù)建立概率模型[4]。上述過程中，對(duì)結(jié)構(gòu)面的分組沒有定量標(biāo)準(zhǔn)，一般根據(jù)數(shù)據(jù)整體分布趨勢(shì)進(jìn)行分組，主要通過極點(diǎn)分布圖或直方圖[5- 8]，而在直方圖統(tǒng)計(jì)中，帶寬的確定通常也是依經(jīng)驗(yàn)而定，對(duì)同一組數(shù)據(jù)的特征，不同帶寬的直方圖的反映是不同的。因此，需要一種定量確定帶寬的方法，提高直方圖統(tǒng)計(jì)的可靠性。其次，對(duì)結(jié)構(gòu)面要素的分布類型的判定過于依靠經(jīng)驗(yàn)，卡方檢驗(yàn)[9-10]是一種很好的方法，但其操作過于復(fù)雜，極大似然估計(jì)相對(duì)簡便。此外，還是有一個(gè)無法回避的問題，即結(jié)構(gòu)面不能完全服從目前已知分布，甚至相差很大。在非參數(shù)估計(jì)方法中，核密度估計(jì)方法[11]在樣本不能較好服從已知分布的情況下，可以模擬得到與樣本較一致的隨機(jī)數(shù)。

鑒于此，本文針對(duì)目前結(jié)構(gòu)面模擬中存在的問題，利用現(xiàn)場(chǎng)采集的結(jié)構(gòu)面數(shù)據(jù)，使用定量方法確定結(jié)構(gòu)面要素的分布類型及數(shù)據(jù)帶寬，通過核密度估計(jì)方法提高結(jié)構(gòu)面模擬精度，并針對(duì)傳統(tǒng)蒙特卡洛方法，提出了新的工作思路，以簡化工作流程。

1 研究區(qū)概況

本文采集數(shù)據(jù)源于我國西南某在建大型水電站壩基右岸巖體(見圖1)。該水電站壩址位于金沙江下游，高原深谷地貌，呈不對(duì)稱“V”字形峽谷。右岸為藥山山脈西坡，山峰高程在3 000 m以上，主要為陡坡與緩坡相間的地形。金沙江由南向北流入壩區(qū)，枯水期水面寬50～115 m，水位約594 m，水深9～18 m不等。拱壩部位高程825 m處谷寬590～713 m。研究區(qū)域主要出露二疊系上統(tǒng)峨眉山組柱狀節(jié)理玄武巖，采用測(cè)線法對(duì)壩基巖體進(jìn)行結(jié)構(gòu)面精測(cè)，得到結(jié)構(gòu)面極點(diǎn)分布，見圖2。

圖1 右岸建基面巖體

圖2 結(jié)構(gòu)面極點(diǎn)分布

2 結(jié)構(gòu)面模擬的蒙特卡洛方法

2.1 方法介紹

蒙特卡洛方法(Monte Carlo method)可以用來解決具有概率解釋的任何問題，因其求解的問題本身具有內(nèi)在的隨機(jī)性，因此必須借助計(jì)算機(jī)的運(yùn)算能力才能模擬這種隨機(jī)的過程。在研究具有概率解釋的問題時(shí)，對(duì)隨機(jī)數(shù)抽樣，使抽樣結(jié)果符合問題的分布形式，從而對(duì)研究的問題進(jìn)行模擬。由于需要大量隨機(jī)數(shù)，蒙特卡洛方法必須依賴于計(jì)算機(jī)。盡管巖體結(jié)構(gòu)面是隨機(jī)分布的，但經(jīng)驗(yàn)表明，結(jié)構(gòu)面的空間分布在一定條件下依然符合一定的概率分布。因此，對(duì)精測(cè)數(shù)據(jù)使用蒙特卡洛方法，便可得到一定深度范圍內(nèi)的結(jié)構(gòu)面特征數(shù)據(jù)。

2.2 工作流程

使用蒙特卡洛方法對(duì)結(jié)構(gòu)面進(jìn)行三維模擬，首先要判定數(shù)據(jù)所服從的分布類型，然后生成服從該分布的隨機(jī)數(shù)[12]。由于巖體受到不同時(shí)期不同類型構(gòu)造的影響，整體不易服從某一確定的分布，所以根據(jù)直方圖起伏和集中趨勢(shì)對(duì)數(shù)據(jù)進(jìn)行分組，每組數(shù)據(jù)單獨(dú)使用蒙特卡洛方法模擬。需要注意的是，每組數(shù)據(jù)在使用蒙特卡洛方法前后，其在數(shù)據(jù)總量中所占的比值應(yīng)保持一致。

2.3 直方圖統(tǒng)計(jì)的帶寬優(yōu)化

帶寬在直方圖統(tǒng)計(jì)中是一個(gè)自由參數(shù)，一般根據(jù)經(jīng)驗(yàn)確定。直方圖統(tǒng)計(jì)是一種對(duì)未知量的非參數(shù)估計(jì)。直方圖中，帶寬過大或過小，得到的直方圖都不能反映問題的真實(shí)情況。帶寬過大，造成直方圖過于平滑，易忽略細(xì)節(jié)；帶寬過小，放大了細(xì)節(jié)在直方圖中的作用，得到的直方圖會(huì)產(chǎn)生較多的尖刺，增加了異常的比重，也不能較好地反映真實(shí)情況。因此，帶寬是否取值得當(dāng)直接決定了直方圖統(tǒng)計(jì)的質(zhì)量，影響后續(xù)論述。為量化帶寬的選取是否得當(dāng)，需要一個(gè)誤差函數(shù)來衡量，本文使用平均積分平方誤差(MISE)作為帶寬選取的誤差，平均積分平方誤差MISE表達(dá)式為[12-13]

(1)

式中，f(x)為未知的真實(shí)的密度分布；fh(x)為基于樣品的以帶寬h為獨(dú)立變量的假定分布；E是樣本的期望。

由于密度函數(shù)f(x)未知，無法直接得到方程的解，可通過計(jì)算機(jī)的多次計(jì)算，尋找該式的最小值，從而代替直接求解方程[14]。本文使用上述方法編寫的程序?qū)掃M(jìn)行優(yōu)化，傾角優(yōu)化后的最佳帶寬為1.72，傾向?yàn)?.06，跡長為0.060 232 7，相對(duì)應(yīng)的最佳區(qū)間數(shù)分別為52、51和50。

2.4 結(jié)構(gòu)面要素的分布擬合優(yōu)度檢驗(yàn)

使用帶寬優(yōu)化之后的數(shù)據(jù)進(jìn)行直方圖統(tǒng)計(jì)(見圖3)，結(jié)果與極點(diǎn)分布圖保持一致，能夠均衡地反映結(jié)構(gòu)面的優(yōu)勢(shì)產(chǎn)狀和發(fā)育較弱的產(chǎn)狀。經(jīng)過帶寬優(yōu)化的直方圖，降低了主觀因素對(duì)直方圖趨勢(shì)的影響。依據(jù)統(tǒng)計(jì)結(jié)果對(duì)數(shù)據(jù)進(jìn)行分組，根據(jù)數(shù)據(jù)的峰值和起伏狀況，將傾角分為[0°， 50°)和[50°， 90°]2個(gè)組，傾向分為[0°， 180°)、[180°， 280°)和[280°， 360°]3個(gè)組，跡長不分組。

圖3 結(jié)構(gòu)面要素統(tǒng)計(jì)直方圖

分布傾向0°～180°AD值P值極大似然比P傾向180°～280°AD值P值極大似然比P傾向280°～360°AD值P值極大似然比P正態(tài)5.634<0.005—2.238<0.005—0.7630.046—Box-Cox變換0.8490.028—1.58<0.005—0.3030.568—對(duì)數(shù)正態(tài)18.57<0.005—1.58<0.005—0.5510.152—三參數(shù)對(duì)數(shù)正態(tài)5.671?01.559?0.9030.34?0.031指數(shù)42.048<0.003—109.343<0.003—39.783<0.003—雙參數(shù)指數(shù)41.772<0.0100.04740.2<0.01005.999<0.0100Johnson變換——————0.2990.578—三參數(shù)對(duì)數(shù)Logistic2.396?01.003?0.528———

對(duì)結(jié)構(gòu)面數(shù)據(jù)進(jìn)行分組后，分別對(duì)每組數(shù)據(jù)進(jìn)行參數(shù)估計(jì)，使用Minitab進(jìn)行分布類型檢驗(yàn)。本文使用Anderson-Darling檢驗(yàn)(簡稱“AD檢驗(yàn)”)進(jìn)行極大似然估計(jì)(當(dāng)能夠進(jìn)行似然比檢驗(yàn)時(shí)顯示極大似然比P，P值是一個(gè)概率，用來度量否定原假設(shè)的證據(jù))。Anderson-Darling統(tǒng)計(jì)量(即AD值，表示假定分布與樣本的差值的平方)用來衡量數(shù)據(jù)遵循特定分布的程度。對(duì)指定的樣本和分布，分布越適合樣本，AD值越小，其表達(dá)式為[15]

(2)

式中，n是樣本數(shù)；f(x)是待檢驗(yàn)的假定分布；fh(x)為樣本的累積分布函數(shù)。本文使用相應(yīng)的P值(當(dāng)存在時(shí))來檢驗(yàn)數(shù)據(jù)是否服從假定分布。如果P值小于選定的閾值α(本文取0.05，但0.05只是一個(gè)常用的閾值，P值小于0.05并不是無效假設(shè)，只能說明其拒絕原假設(shè)的理由更充分)，則拒絕數(shù)據(jù)來自該分布的假設(shè)。P值在數(shù)學(xué)上不存在時(shí)，Minitab軟件就不顯示AD檢驗(yàn)的P值。對(duì)研究區(qū)結(jié)構(gòu)面要素進(jìn)行分布擬合優(yōu)度檢驗(yàn)，得到結(jié)構(gòu)面要素的概率分布圖(默認(rèn)置信區(qū)間為95%)。對(duì)6組數(shù)據(jù)進(jìn)行分布類型檢驗(yàn)后的結(jié)果見表1～3。表中， *表示無法計(jì)算的缺失值。

表1 傾角分布類型檢驗(yàn)結(jié)果

表3 跡長分布類型檢驗(yàn)結(jié)果

綜合AD值小、P值大的原則，確定原數(shù)據(jù)分布類型。其中，傾向0°～180°區(qū)間較符合Box-Cox變換后正態(tài)分布，傾向180°～280°區(qū)間較符合三參數(shù)對(duì)數(shù)正態(tài)分布，傾向280°～360°區(qū)間較符合Johnson 變換后正態(tài)分布；傾角0°～50°區(qū)間較符合Box-Cox變換后正態(tài)分布，傾角50°～90°區(qū)間較符合Box-Cox變換后正態(tài)分布；跡長服從Johnson變換后正態(tài)分布。各自對(duì)應(yīng)的在95%置信區(qū)間的累積概率分布見圖4～6。

圖4 傾角累積概率分布

圖5 傾向累積概率分布

圖6 跡長累積概率分布

分布擬合優(yōu)度檢驗(yàn)結(jié)果表明，僅跡長在未經(jīng)分組的情況下，對(duì)假定分布類型的接受程度較高。巖體作為非均質(zhì)體存在，在各種地質(zhì)營力的作用下，被結(jié)構(gòu)面切割，在一定區(qū)域內(nèi)形成幾何相似的塊體，因此在統(tǒng)計(jì)量沒有大到足以代表整個(gè)巖體結(jié)構(gòu)面時(shí)，精測(cè)數(shù)據(jù)應(yīng)該是非連續(xù)數(shù)據(jù)。在特定地質(zhì)營力作用下，形成的數(shù)量相對(duì)較多、且易于被觀察和測(cè)量的部分結(jié)構(gòu)面，即為優(yōu)勢(shì)結(jié)構(gòu)面。優(yōu)勢(shì)結(jié)構(gòu)面數(shù)量較大，易于被觀測(cè)，因此優(yōu)勢(shì)結(jié)構(gòu)面數(shù)據(jù)在精測(cè)數(shù)據(jù)占比較大。林德伯格-費(fèi)勒中心極限定理亦表明，結(jié)構(gòu)面跡長作為描述結(jié)構(gòu)面出露程度的標(biāo)量，其原始數(shù)據(jù)(未分組)應(yīng)該服從正態(tài)分布。

2.5 模擬結(jié)果

本次模擬在Jupyter notebook中完成。Jupyterno-tebook是基于Python的交互式解釋器。整個(gè)程序設(shè)計(jì)過程為：首先，程序初始化后導(dǎo)入數(shù)據(jù)，設(shè)置并初始化變量，對(duì)數(shù)據(jù)分組；然后，確定模擬的各組數(shù)據(jù)的總量，根據(jù)各組數(shù)據(jù)對(duì)應(yīng)的變換方式分別對(duì)每組數(shù)據(jù)進(jìn)行變換，使用NumPy工具庫生成偽隨機(jī)數(shù)，接著對(duì)生成的偽隨機(jī)數(shù)反變換，得到模擬結(jié)果，見圖7。

對(duì)比圖3，使用蒙特卡洛方法模擬的結(jié)果，在集散趨勢(shì)上與原數(shù)據(jù)保持一致。經(jīng)過帶寬優(yōu)化后，為避免主觀因素的影響，使用量化方法判斷精測(cè)數(shù)據(jù)的分布類型。通過分布類型檢驗(yàn)，獲得了數(shù)據(jù)的一般分布類型，但由于細(xì)節(jié)的損失，模擬結(jié)果無法較好再現(xiàn)原數(shù)據(jù)的分布特征。

圖7 模擬結(jié)果直方圖

圖8 原始數(shù)據(jù)與核密度估計(jì)結(jié)果對(duì)比

3 結(jié)構(gòu)面模擬的核密度估計(jì)

使用蒙特卡洛方法進(jìn)行結(jié)構(gòu)面模擬時(shí)，需要對(duì)數(shù)據(jù)進(jìn)行分布類型檢驗(yàn)，而模擬結(jié)果的好壞取決于數(shù)據(jù)對(duì)假定分布類型的服從程度。對(duì)數(shù)據(jù)切割分組的不連續(xù)，可能降低模擬結(jié)果的準(zhǔn)確性。為提高模擬結(jié)果與原始數(shù)據(jù)的相似度，使用核密度估計(jì)方法對(duì)結(jié)構(gòu)面數(shù)據(jù)進(jìn)行模擬。

3.1 核密度估計(jì)

核密度估計(jì)是估計(jì)隨機(jī)變量的概率密度函數(shù)的非參數(shù)方法，這意味對(duì)隨機(jī)變量的估計(jì)不再需要固定的參數(shù)，而且程序設(shè)計(jì)可以減少很多工作。在直方圖統(tǒng)計(jì)中，對(duì)于區(qū)間(xk-1，xk]((xk-1，x]為其中某1個(gè)條帶)中的所有樣本都同等對(duì)待。如在[0，1]區(qū)間，有2個(gè)樣本取值0+和1-，前者表達(dá)的信息是分布在0周圍有密度，而后者應(yīng)該在1周圍有密度，兩者存在不同的含義。因此，對(duì)于直方圖統(tǒng)計(jì)，給定寬度h，樣本Xi應(yīng)該反應(yīng)區(qū)間[Xi-h/2，Xi+h/2](該區(qū)間的分布函數(shù)為I(-h/2≤x≤h/2))的信息，定義[16]

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

(3)

對(duì)應(yīng)的估計(jì)fh(x)為

(4)

再進(jìn)一步把寬度h看成參數(shù)，直接定義K(x)=I(|x|≤1)/2，若該區(qū)間的權(quán)重相等，則估計(jì)fh(x)為

(5)

(6)

Kh(· )=K(·/h)/h稱為核函數(shù)，理論上任何非負(fù)單峰的概率密度函數(shù)都可以作為核函數(shù)。不同的核函數(shù)在對(duì)目標(biāo)進(jìn)行核密度估計(jì)時(shí)會(huì)產(chǎn)生不同估計(jì)模型，在每個(gè)估計(jì)區(qū)間體現(xiàn)核函數(shù)原本的性質(zhì)。目前，常用的核函數(shù)是高斯(Gaussian)核函數(shù)，表示為K(x;h)∝exp[-x2/2h2]。

3.2 結(jié)構(gòu)面模擬的核密度估計(jì)

使用核密度估計(jì)方法，依然要先對(duì)數(shù)據(jù)進(jìn)行帶寬優(yōu)化。使用優(yōu)化后的帶寬，在Jupyter notebook中，利用Scikit-learn工具計(jì)算。為提高采樣精度，Scikit-learn中的KDE模塊引入了最近鄰算法中的KDTree和BallTree，使采樣模型更貼近原始數(shù)據(jù)。使用Gaussian內(nèi)核的核密度估計(jì)方法對(duì)精測(cè)數(shù)據(jù)模擬。原始數(shù)據(jù)與核密度估計(jì)結(jié)果對(duì)比見圖8。圖8中，直方圖為原始數(shù)據(jù)統(tǒng)計(jì)結(jié)果，曲線為核密度估計(jì)。核密度估計(jì)模擬結(jié)果見圖9。

圖9 核密度估計(jì)模擬結(jié)果

3.3 結(jié)果對(duì)比

從圖8可知，高斯內(nèi)核的核密度估計(jì)曲線較為平滑，其概率密度的變化均為漸變曲線。而直方圖統(tǒng)計(jì)中，相鄰區(qū)首尾相接處的數(shù)據(jù)是躍變的，而實(shí)際上在區(qū)間連接處數(shù)據(jù)應(yīng)該是連續(xù)的。因此，使用核密度估計(jì)方法對(duì)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)，能夠提高模擬的精度。核密度估計(jì)曲線與直方圖統(tǒng)計(jì)數(shù)據(jù)增減趨勢(shì)一致，但不完全貼近方塊，為模擬結(jié)果預(yù)留了變化區(qū)間，表明核密度估計(jì)的結(jié)果受偶然因素的影響較小，可以減少過擬合的發(fā)生。在保證帶寬為最佳的前提下，使用核密度估計(jì)既能夠保留細(xì)部特征，又能夠兼顧整體分布趨勢(shì)。

從圖9可知，對(duì)比原始數(shù)據(jù)和蒙特卡洛方法的結(jié)果，核密度估計(jì)結(jié)果更符合原始數(shù)據(jù)的集散趨勢(shì)，并且對(duì)細(xì)部的反映更真實(shí)。觀察蒙特卡洛模擬結(jié)果，其數(shù)據(jù)整體較為規(guī)律，數(shù)據(jù)的整體分布趨勢(shì)符合公式生成數(shù)據(jù)的特點(diǎn)。而核密度估計(jì)方法模擬的結(jié)果更加貼近原數(shù)據(jù)方塊。在整體趨勢(shì)上，兩者的結(jié)果較為相似，但在局部數(shù)據(jù)的呈現(xiàn)上，核密度估計(jì)更加符合原數(shù)據(jù)的特點(diǎn)，對(duì)細(xì)節(jié)的把握更充分。

在工作流程上，蒙特卡洛方法的工作量更大，實(shí)施步驟更多。核密度估計(jì)只需要進(jìn)行帶寬優(yōu)化、直方圖統(tǒng)計(jì)、核密度估計(jì)、隨機(jī)數(shù)生成，最后成果檢驗(yàn)，相對(duì)而言，工作量大大減少。從量化的角度看，蒙特卡洛方法在分組時(shí)存在主觀判斷，可能對(duì)模擬的成果有一定影響。因此，在結(jié)構(gòu)面三維模擬中，使用核密度估計(jì)方法，有利于模擬精度的提高，同時(shí)減少工作量。

4 結(jié) 語

本文以西南某大型在建水電站建基面結(jié)構(gòu)面精測(cè)數(shù)據(jù)為基礎(chǔ)，對(duì)結(jié)構(gòu)面要素進(jìn)行統(tǒng)計(jì)分析，使用平均積分平方誤差優(yōu)化直方圖統(tǒng)計(jì)帶寬。引入核密度估計(jì)方法，對(duì)結(jié)構(gòu)面模擬的統(tǒng)計(jì)方法進(jìn)行改進(jìn)，對(duì)比蒙特卡洛方法的流程和結(jié)果，得出以下結(jié)論：

(1)使用平均積分平方誤差對(duì)直方圖帶寬進(jìn)行優(yōu)化后，直方圖統(tǒng)計(jì)結(jié)果與極點(diǎn)分布圖保持一致，能夠照顧到結(jié)構(gòu)面的優(yōu)勢(shì)產(chǎn)狀和發(fā)育較弱的產(chǎn)狀。經(jīng)過帶寬優(yōu)化的直方圖，能避免主觀因素對(duì)直方圖的影響，為結(jié)構(gòu)面分布類型的確定和核密度估計(jì)做好鋪墊。

(2)在蒙特卡洛模擬中，采用了極大似然估計(jì)，根據(jù)各組數(shù)據(jù)對(duì)應(yīng)的AD值和P值，判定各組數(shù)據(jù)對(duì)假定分布類型的服從程度，從而確定分布類型。蒙特卡洛方法的結(jié)果表明，其把握數(shù)據(jù)的整體分布趨勢(shì)較好，但細(xì)節(jié)丟失嚴(yán)重。

(3)在探索提升結(jié)構(gòu)面模擬精度的過程中，引入核密度估計(jì)方法，結(jié)果表明，其結(jié)果能夠較好地還原數(shù)據(jù)的整體分布特征。此方法省略了分布類型檢驗(yàn)等步驟，忽略了主觀因素對(duì)結(jié)果的影響，提高了模擬精度，減少了工作量。