尚振華，徐必根

(長(zhǎng)沙礦山研究院金屬礦山安全技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖南長(zhǎng)沙410012)

基于蒙特卡羅模擬的巖石點(diǎn)載荷試驗(yàn)數(shù)據(jù)處理*

尚振華，徐必根

(長(zhǎng)沙礦山研究院金屬礦山安全技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖南長(zhǎng)沙410012)

結(jié)合某鐵礦采空區(qū)研究的巖石力學(xué)試驗(yàn)，以蒙特卡羅模擬處理點(diǎn)荷載試驗(yàn)結(jié)果，利用概率論與統(tǒng)計(jì)理論，結(jié)合Excel插件@Risk對(duì)巖石點(diǎn)載荷強(qiáng)度指標(biāo)進(jìn)行分析及最終取值。該方法采用計(jì)算機(jī)多次模擬，較傳統(tǒng)的處理方式易分析、易理解，同時(shí)提高了點(diǎn)載荷強(qiáng)度的精度與可信度，另外受問(wèn)題條件限制的影響較小，可很好的提高可行性。

巖石強(qiáng)度;點(diǎn)載荷試驗(yàn);試驗(yàn)數(shù)據(jù);蒙特卡羅;概率分布;@Risk

點(diǎn)荷載試驗(yàn)由于具有現(xiàn)場(chǎng)操作簡(jiǎn)單、取樣要求不高和試驗(yàn)周期短等特點(diǎn)，在巖土工程中得到越來(lái)越多的應(yīng)用［1］。一般情況下，試驗(yàn)數(shù)據(jù)按傳統(tǒng)方式處理后將其結(jié)果作為準(zhǔn)確值的一個(gè)估計(jì)，但該估計(jì)值與實(shí)際值之間將會(huì)產(chǎn)生一個(gè)估計(jì)誤差，而這個(gè)誤差將會(huì)或多或少地帶來(lái)一定的損失［2］。另外，在試驗(yàn)過(guò)程中，由于試驗(yàn)條件的限制使得操作無(wú)法嚴(yán)格按照設(shè)計(jì)執(zhí)行，以致一些外部因素使試驗(yàn)結(jié)果出現(xiàn)了隨機(jī)性，此外，巖體的空間變異性和不同的隸屬關(guān)系使得工程巖組的分類具有了不確定性。

1 傳統(tǒng)方式的數(shù)據(jù)處理

在現(xiàn)實(shí)中，一方面，并不能保證試驗(yàn)過(guò)程完全按照期望的要求進(jìn)行，另一方面，過(guò)于保守或者過(guò)于樂(lè)觀的估計(jì)結(jié)果所包含的錯(cuò)誤綜合在一起，常常會(huì)導(dǎo)致所得結(jié)果與真實(shí)結(jié)果有很大差異，于是根據(jù)該“預(yù)期”方式做出的試驗(yàn)值可能是錯(cuò)誤的。

傳統(tǒng)試驗(yàn)結(jié)果處理方法是通過(guò)綜合分析模型變量的單“點(diǎn)”估計(jì)來(lái)預(yù)測(cè)一種結(jié)果，采用標(biāo)準(zhǔn)Excel電子表格進(jìn)行簡(jiǎn)單的整理、計(jì)算，但是該計(jì)算的變量只能估計(jì)，因?yàn)闊o(wú)法確切得知實(shí)際將產(chǎn)生的值。

某鐵礦采空區(qū)研究中要進(jìn)行巖石力學(xué)實(shí)驗(yàn)，但由于該鐵礦圍巖之一花崗斑巖風(fēng)化較為嚴(yán)重，較難取得完整的巖芯試樣，于是在該礦+130m中段、+100m中段和+70m中段分別采集了部分試樣，在長(zhǎng)沙礦山研究院巖石力學(xué)實(shí)驗(yàn)室利用A125型攜帶式點(diǎn)載荷儀進(jìn)行點(diǎn)載荷強(qiáng)度試驗(yàn)，按照傳統(tǒng)的ISRM法及回歸分析計(jì)算得出花崗斑巖的點(diǎn)載荷強(qiáng)度指標(biāo)Is(50)=2.70，抗壓強(qiáng)度σc=64.80MPa。

2 蒙特卡羅分析的基本原理

蒙特卡羅(MonteCarlo)方法，又稱計(jì)算機(jī)隨機(jī)模擬方法，是一種基于“隨機(jī)數(shù)”的計(jì)算方法。這一方法源于美國(guó)在第二次世界大戰(zhàn)研制原子彈的“曼哈頓計(jì)劃”，該計(jì)劃的主持人之一數(shù)學(xué)家馮·諾伊曼用馳名世界的賭城——摩納哥的MonteCarlo來(lái)命名這種方法［3］。

它是用一系列隨機(jī)數(shù)來(lái)近似解決問(wèn)題的一種方法，是通過(guò)尋找一個(gè)概率統(tǒng)計(jì)的相似體并用試驗(yàn)取樣過(guò)程來(lái)獲得該相似體的近似解的處理數(shù)學(xué)問(wèn)題的一種手段。運(yùn)用該近似方法所獲得的問(wèn)題的解更接近于試驗(yàn)結(jié)果，而不是經(jīng)典數(shù)值計(jì)算結(jié)果。

它的基本思想是為了求解數(shù)學(xué)、物理、工程技術(shù)以及生產(chǎn)管理等方面的問(wèn)題，首先建立一個(gè)概率模型或隨機(jī)過(guò)程，使它的參數(shù)等于問(wèn)題的解;然后通過(guò)對(duì)模型或過(guò)程的觀察或抽樣試驗(yàn)來(lái)計(jì)算所求參數(shù)的統(tǒng)計(jì)特征;最后給出所求解的近似值，而求解的確定精確度可用估計(jì)值的標(biāo)準(zhǔn)誤差來(lái)表示［4］。

假設(shè)所要求的量x是隨機(jī)變量ε的數(shù)學(xué)期望E(ε)，那么近似確定x的方法是對(duì)ε進(jìn)行N次重復(fù)抽樣，產(chǎn)生相互獨(dú)立的ε值的序列ε1、ε2、…、εn，并計(jì)算其算術(shù)平均值:

根據(jù)柯爾莫哥羅夫(Komogorov)的加強(qiáng)大數(shù)定律有:

3 基于蒙特卡羅模擬的數(shù)據(jù)處理

3.1 試驗(yàn)數(shù)據(jù)的分布擬合

Excel插件@Risk采用蒙特卡羅模擬，對(duì)采集到的試驗(yàn)數(shù)據(jù)(巖心直徑的平方、試樣破壞時(shí)的總荷載P)應(yīng)用MicrosoftExcel的插件@Risk進(jìn)行分布擬合，利用蒙特卡羅抽樣得出二者分別基本符合三角分布(Triangledistribution)和威布爾分布(Weibull distribution)，見圖1、圖2。

圖1 巖心直徑的平方的分布擬合

圖2 試樣破壞時(shí)的總荷載P的分布擬合

同樣，對(duì)整理的數(shù)據(jù)(標(biāo)準(zhǔn)點(diǎn)載荷強(qiáng)度Is(50)、抗壓強(qiáng)度σc)應(yīng)用相同方法得出二者均基本符合正態(tài)分布(Normaldistribution)，見圖3、圖4。

圖3 標(biāo)準(zhǔn)點(diǎn)載荷強(qiáng)度Is(50)的分布擬合

圖4 抗壓強(qiáng)度σc的分布擬合

3.2 數(shù)學(xué)模型的建立

(1)定義輸入變量及相關(guān)。將分布擬合得出的花崗斑巖的各個(gè)統(tǒng)計(jì)量(、P、Is(50)、σc)分別寫入Excel不同的單元格，即對(duì)這些統(tǒng)計(jì)量的不確定值指定概率分布函數(shù)，該分布包含了每個(gè)參數(shù)的具體情況，如Is(50)服從正態(tài)分布RiskNormal(2.78261，0.33965)。利用傳統(tǒng)的Excel功能對(duì)和P進(jìn)行數(shù)值擬合，得出二者相關(guān)系數(shù)為0.89。然后通過(guò)@Risk的“定義相關(guān)”功能，以散點(diǎn)矩陣的方式輸入并以散點(diǎn)圖的方式顯示二者之間的相關(guān)性。

(2)定義輸出變量。利用@Risk的分布函數(shù)RiskOutput()輸出重點(diǎn)關(guān)注的點(diǎn)載荷的強(qiáng)度指標(biāo)Is(50)和σc。其中，Is(50)的輸出為RiskOutput("Is(50)(MPa)")+的輸出為RiskOutput("σc(MPa)")+24×Is(50)。這里添加的RiskOutput函數(shù)中的“+”，表示為根據(jù)輸出函數(shù)需要來(lái)選擇合適的單元格作為模擬輸出，并將輸出分別命名為“Is(50)(MPa)”和“σc(MPa)”。

3.3 模擬計(jì)算及結(jié)果分析

該模型利用@Risk對(duì)采集到的試驗(yàn)數(shù)據(jù)應(yīng)用基于最優(yōu)化的蒙特卡羅法進(jìn)行分析處理，并在Excel中以概率分布的形式用圖形顯示可能的多種結(jié)果及其概率，這允許在存在不確定因素的情況下做出最好的點(diǎn)載荷試驗(yàn)處理結(jié)果。

此次點(diǎn)載荷試驗(yàn)的數(shù)據(jù)處理結(jié)果在蒙特卡羅模擬的基礎(chǔ)上執(zhí)行5個(gè)循環(huán)，每個(gè)循環(huán)迭代1000次，計(jì)算結(jié)果見圖5，計(jì)算結(jié)果的分布擬合見圖6。

圖5 σc的蒙特卡羅模擬結(jié)果

圖6 σc模擬結(jié)果的分布擬合

該計(jì)算結(jié)果基本符合對(duì)數(shù)正態(tài)分布(Logarithmicnormaldistribution)。假設(shè)給定顯著性水平n=0.05，即置信水平為95%時(shí)，由蒙特卡羅方法得出的置信區(qū)間為［45.7，95.0］，最終σc的取值為期望值67.60MPa。較傳統(tǒng)方式處理結(jié)果提高了2.80 MPa，約4.32%，也進(jìn)一步驗(yàn)證了傳統(tǒng)的“點(diǎn)”估計(jì)在該次試驗(yàn)結(jié)果的處理過(guò)程中略偏保守。

通過(guò)繪制“龍卷風(fēng)”圖(見圖7)，很直觀地得出對(duì)σc的取值影響較大的參數(shù)系數(shù)分別為Is和，影響系數(shù)分別為0.9和0.3，其中π。這也進(jìn)一步驗(yàn)證了點(diǎn)載荷試驗(yàn)中巖石試樣的形狀對(duì)最終的試驗(yàn)結(jié)果有較大的影響。

運(yùn)用@Risk對(duì)所取得的試驗(yàn)數(shù)據(jù)進(jìn)行多種組合的模擬計(jì)算，得到的模擬結(jié)果一方面提高了試驗(yàn)數(shù)據(jù)處理的精度，另一方面可以以更直觀的角度顯示多種試驗(yàn)組合帶來(lái)的不同結(jié)果。在該模擬中通過(guò)對(duì)模擬結(jié)果的查看，至少可以得到1000組試驗(yàn)數(shù)據(jù)及組合，這對(duì)進(jìn)一步了解和分析巖石點(diǎn)載荷試驗(yàn)過(guò)程中的各個(gè)環(huán)境變量對(duì)其最終結(jié)果的取值有一定的幫助。

圖7 對(duì)σc回歸分析的“龍卷風(fēng)”

4 結(jié)論

點(diǎn)載荷試驗(yàn)是巖體工程中巖體分級(jí)(RMR，CSMR等)研究和強(qiáng)度參數(shù)的一個(gè)十分重要的指標(biāo)，但由于傳統(tǒng)方法以“點(diǎn)”估計(jì)的方法對(duì)試驗(yàn)結(jié)果處理，常常導(dǎo)致實(shí)際情況的結(jié)果與估計(jì)的結(jié)果有很大差異。借助@Risk基于蒙特卡羅法對(duì)后期試驗(yàn)數(shù)據(jù)分析處理，并針對(duì)現(xiàn)實(shí)情況中存在的不確定因素加以組合、模擬，能以更加直觀、準(zhǔn)確的方式顯示所有可能“預(yù)期”的結(jié)果。

@Risk基于蒙特卡羅模擬以概率統(tǒng)計(jì)理論為基礎(chǔ)，以隨機(jī)抽樣(隨機(jī)變量的抽樣)為手段，在很多方面有重要的應(yīng)用。它在巖石點(diǎn)載荷試驗(yàn)數(shù)據(jù)處理方面的優(yōu)點(diǎn)表現(xiàn)在3個(gè)方面:方法和程序的結(jié)構(gòu)簡(jiǎn)單，易分析、易理解;收斂的概率性和收斂速度與問(wèn)題的維數(shù)無(wú)關(guān)，很好地避免了維數(shù)問(wèn)題;受問(wèn)題條件限制的影響較小，很好地提高可行性。

［1］李茂蘭，鐘光宙.巖石點(diǎn)荷載試驗(yàn)及其應(yīng)用［M］.西安:西安交通大學(xué)出版社，1994.

［2］譚文輝，喬蘭，李治平.巖石點(diǎn)荷載強(qiáng)度指標(biāo)取值的模糊隨機(jī)方法［J］.金屬礦山，2001，(10).

［3］王坤.MonteCarlo方法及其簡(jiǎn)單應(yīng)用［J］.科技信息，2010，(14).

［4］徐鐘濟(jì).蒙特卡羅方法［M］.上海:上?？茖W(xué)技術(shù)出版社，1985:5－6.

863計(jì)劃重點(diǎn)課題(2008AA062104).

2011－09－08)

尚振華(1985－)，男，河南安陽(yáng)人，碩士研究生，主要從事礦山巖石力學(xué)及地壓監(jiān)控技術(shù)研究，Email:szh207@gmail.com。