常天龍

(中國鐵路設計集團有限公司，天津 300251)

隨著三維數(shù)值模擬逐漸得到推廣與應用，三維地質建模技術快速發(fā)展，在巖土工程和地質工程領域中得到了廣泛應用，但是現(xiàn)在的地質建模仍是主要依靠一些商用軟件如：Midas、ANSYS、ABAQUS等軟件中提供的建模模塊實現(xiàn)三維建模。人們起初對三維地質建模的要求是實現(xiàn)建立的三維地質模型可以真實反映出研究對象包括的地層及不良地質體，而隨著數(shù)值模擬的發(fā)展，為了提高數(shù)值計算結果的精確度，要求所建立的模型不僅要能真實反映研究對象的地質特征，而且模型還要具有較高質量的網(wǎng)格剖分度[1]。為了實現(xiàn)建模的效率與計算的準確性，往往為了突出研究的問題而采取了部分“有重點”的建模方法。

在三維數(shù)值模擬地質建模方面，許多學者提出了自己的方法，Jones提出了利用鉆孔與剖面結合的實體建模方法為MODFLOW提供計算網(wǎng)格[4]，王春祥等提出L-W拓撲模型建立三維地質模型并實現(xiàn)地質信息與有限元計算的結合[6]，王明華等通過研究層狀巖體的分布特征，提出基于松散模式的三維規(guī)則網(wǎng)格與FLAC3D之間的模型轉化方法[5]。李明超等結合數(shù)值模擬計算的特點提出地質建模與數(shù)值模擬耦合的方法來進行前處理的思路[7]。以上研究成果為三維地質建模服務于數(shù)值模擬計算提供了很多新思路、新方法，為未來通用程序的開發(fā)提供了不同的思路。

不同學者開發(fā)的計算程序多將滑面近似圓弧形去處理，對于某些非圓弧滑面邊坡往往很難得到準確的計算結果，本文結合前人的研究成果，提出了一種基于四棱柱網(wǎng)格單元的建模方法，結合Matlab矩陣計算的特點，實現(xiàn)數(shù)據(jù)的矩陣計算與存儲。依據(jù)地質體地層呈層分布特點，本文采用矩陣計算插分得到網(wǎng)格坐標，通過單元節(jié)點三維坐標間的運算求取四棱柱的體積、單元表面積等參數(shù)，從而建立一種基于四棱柱單元可以求取任意形狀滑面邊坡穩(wěn)定性的地質模型。

1 構建三維地質模型

1.1 網(wǎng)格節(jié)點坐標計算

通過現(xiàn)場測量得到邊坡表面地面點及勘探鉆孔坐標，將同一地層內的控制點三維坐標放入同一文本文件中作為插值點備用。Matlab中提供了四種插值函數(shù)1.nearest函數(shù)：最鄰近插值方法(nearest neighbor interpolation)；2.linear函數(shù)：線性插值(linear interpolation)；3.spline函數(shù)：三次樣條插值(cubic spline interpolation)；4. v5cubic函數(shù)：Matlab程序自帶插值計算函數(shù)。根據(jù)不同情況，可以根據(jù)地形起伏變化的大小選用不同的插值函數(shù)進行網(wǎng)格插值劃分，筆者在程序設計過程中對比不同的插值函數(shù)結果發(fā)現(xiàn)使用“v5cubic”函數(shù)進行插值可以得到比較平滑的曲面，減少表面凸點的出現(xiàn)，使得數(shù)據(jù)奇點更少，提高計算精度。使用自編的插值程序，各個地層按照自定義網(wǎng)格尺寸插值計算得出各個曲面的網(wǎng)格角點三維坐標。由于Matlab是矩陣運算程序，所有插值得到的網(wǎng)格角點坐標以矩陣的形勢保存，方便讀取修改(見圖1)。

圖1 插值得到的曲面網(wǎng)格

1.2 網(wǎng)格單元數(shù)據(jù)計算

在同一邊坡不同位置上滑面所在地層不盡相同，計算某滑面以上不同地層的總重量既為該滑面單元的滑體重量，對于滑面所在單元，滑面切割該單元為兩部分，滑面以上部分計入滑體重量，滑面以下需要剔除。針對不同的滑面與四棱柱單元相交情況可把單元分為三類12種單元。筆者按照滑面四個角點的相對位置分為：“四點同層”，“一點凸出”，“兩點凸出斜交”三種單元類型。(圖2中a為四點同層型，b為一點凸出型，c為兩點凸出斜交型)從圖中可以看出，假設滑面為H層，A、B、C層為地層面，H層平面矩形單元角點坐標沿順時針方向分別為H(i,j)、H(i+1,j)、H(i+1,j+1)、H(i,j+1)其他層角點坐標角標相同不予贅述，每個層面網(wǎng)格坐標構成規(guī)則的坐標矩陣。為了滿足不出現(xiàn)滑面跨層現(xiàn)象，網(wǎng)格尺寸需滿足：maxΔH﹤2h或2h/a>tanθ且2h/b>tanθ，其中a、b為網(wǎng)格長寬，h為分層層厚；θ為滑面最大傾角；為滑面相鄰兩個網(wǎng)格角點高程之差；只需要將網(wǎng)格劃分足夠細，避免出現(xiàn)網(wǎng)格跨度過大從而導致跨層網(wǎng)格出現(xiàn)即可。下面就對上述三類網(wǎng)格體積計算公式進行推導：

未被滑面切割的一般單元體積為：

V=0.25ab[∑A-∑B]

(1)

1.3 滑面所在單元體積計算方法

1.3.1 四點同層類網(wǎng)格

判斷條件為：B(i，j)

(2)

1.3.2 一點凸出類網(wǎng)格

判斷條件為：H(i,j)>B(i,j),H(i,j+1)

六棱柱A(i,j)A(i+1,j)A(i+1,j+1)A(i,j+1)H(i,j)H(i+1,j)H(i+1,j+1)H(i,j+1)體積為：

(3)

滑面H(i,j)H(i+1,j)H(i+1,j+1)H(i,j+1)面積為：

(4)

1.3.3 兩點凸出斜交型

判斷條件為：B(i,j)>H(i,j),B(i,j+1)>H(i,j+1),B(i+1,j)

滑面單元H(i,j)H(i+1,j)H(i+1,j+1)H(i,j+1)與單元面B(i,j)B(i+1,j)B(i+1,j+1)B(i,j+1)交點為O1、O2。

V=VAB-Vu+Vd

(5)

其中：VAB為滑面以上A、B單元體積，Vu為上三棱柱H(i+1,j)H(i+1,j+1)B(i+1,j)B(i+1,j+1)O1O2體積，Vd位下三棱柱H(i,j)H(i,j+1)B(i,j)B(i,j+1)O1O2體積

V=0.25ab[∑A-∑B]

(6)

(7)

(8)

a為四點同層型，b為一點凸出型，c為兩點凸出斜交型

以上給出了三類網(wǎng)格某種情況下的體積計算公式，其他同類型網(wǎng)格只需要調整各個角點的在公示中的順序即可得到相應的計算公式，在此不予以一一贅述。通過Matlab中循環(huán)判斷語句實現(xiàn)上文中的判斷計算過程，計算出各個網(wǎng)格的體積后乘以該單元密度即可得到單元重度。

1.4 滑面單元面積及傾角

在求解安全系數(shù)時還需要輸入每個滑面網(wǎng)格的表面積，滑面在XOY面上的投影即為a*b(a、b分別為設定的網(wǎng)格長寬)，通過求取每個滑面單元與水平面的夾角即可求得滑面的面積。面積和夾角公式為式9、10。

(9)

(10)

式中：α為滑面與水平面夾角

通過本節(jié)中各項計算公式，可以得到穩(wěn)定分析中所需的單元重量，滑面傾角，滑面面積等所需數(shù)據(jù)，從而實現(xiàn)三維地質建模，數(shù)據(jù)以矩陣形式存儲與調用。

2 結語

Matlab是基于矩陣計算的一種用于算法開發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計算的高級技術計算語言和交互式算法開發(fā)程序，本文介紹了一種基于Matlab自編程序實現(xiàn)地層網(wǎng)格化的方法，從而得到網(wǎng)格節(jié)點坐標，經過文中推導的一系列坐標運算公式，可以快速的實現(xiàn)網(wǎng)格單元特性計算，從而實現(xiàn)數(shù)值分析前處理三維地質建模的目的，本方法建立的三維地質模型不受滑面形狀所限，適用于任意形狀滑面邊坡計算，與邊坡實際地層較吻合。該算法具有較好的實用性和可拓展性，通過自編不同功能的算法可

以實現(xiàn)不同的計算目的。

[1]徐能雄. 適于數(shù)值模擬的三維工程地質建模方法[J]. 巖土工程學報.2009.31(11): 1710-1715.

[2]胡斌,張倬元,黃潤秋,等. FLAC3D 前處理程序的開發(fā)與仿真效果檢驗[J]. 巖石力學與工程學報.2002.21(9):1381-1397.

[3]張 煜, 溫國強, 王笑海. 三維體繪制技術在工程地質可視化中的應用[J]. 巖石力學與工程學報.2002.21(4): 563-567.

[4]LEMONAM, JONES. Building solid models from boreholes and user-defined cross-sections[J]. Computers & Geosciences, 2003.29(5): 547-555.

[5]王明華, 白 云. 層狀巖體三維可視化構模與數(shù)值模擬的集成研究[J]. 巖土力學.2005.26(7): 1123-1126.

[6]王純祥, 白世偉. 三維地層信息系統(tǒng)與有限元方法集成研究[J]. 巖石力學與工程學報.2004.23(21): 3695-3699.

[7]李明超, 鐘登華, 等. 基于三維地質模型的工程巖體結構精細數(shù)值建模[J]. 巖石力學與工程學報.2007.26(9): 1893-1898.

[8]Xing Z. Three-dimensional stability analysis of concave slopes in plan view[J]. Journal of the Geotechnical Engineering, 1988.114(6): 658-671.

[9]徐能雄, 武 雄, 汪小剛,等. 基于三維地質建模的復雜構造巖體六面體網(wǎng)格剖分方法[J]. 巖土工程學報.2006.28(8): 957-961 .

[10]谷天峰, 王家鼎, 王念秦. 呂梁機場黃土滑坡特征及其三維穩(wěn)定性分析[J]. 巖土力學.2013.