巖體結構面模擬的統計方法改進研究

張 斌，黃 俊

(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)在探索提升結構面模擬精度的過程中，引入核密度估計方法，結果表明，其結果能夠較好地還原數據的整體分布特征。此方法省略了分布類型檢驗等步驟，忽略了主觀因素對結果的影響，提高了模擬精度，減少了工作量。