基于自由面上環路對應節理無限切割的關鍵塊體搜索

石 露，李小春，白 冰

（中國科學院武漢巖土力學研究所 巖土力學與工程國家重點實驗室，武漢 430071 ）

1 引 言

自然界的巖體經歷了多次反復的地質作用，經受過變形破壞，是由一定巖性的巖塊和結構面組成的地質體[1]。結構面的存在是巖體的力學性質不同于完整巖石材料的根本原因。在組成巖體的兩大基本單元中，巖塊和結構面所處的地位不同，從強度上來比較，結構面強度遠小于巖塊的強度，為巖體中的薄弱部位，常常是控制巖體強度和變形的關鍵因素。從結構上來講，巖塊的形態、規模等特征由結構面所控制，完全取決于結構面的組合、交切與圍限的情況。因此，在巖體的結構中，通常結構面占據著更為重要的位置[2]。

關鍵塊體理論是由 Goodman和石根華[3]于1982年創立，具有完善的理論體系，包括有限性定理和可移動性定理，目前在分析地下硐室、邊坡和基礎穩定性中得到廣泛的應用。塊體理論創立之初即尋找工程在自重條件下的不穩定的關鍵塊體，并假定結構面無限大。鑒于工程中，除層面或斷層在研究區域內可視為無限平面，多數遍布節理均是有限長的，由此后續發展了基于三維隨機網絡的巖體關鍵塊體搜索[4－6]，可對塊體的形態、大小及空間分布特征等進行統計分析。目前，塊體理論已發展成全空間內的塊體識別[7－10]，不僅可用于關鍵塊體分析，也可用于三維巖體不連續變形分析。

雖然關鍵塊體理論在巖質邊坡、硐室的穩定性評價中起到重要的作用，但目前尚無通用軟件、程序實現困難以及節理眾多造成的程序運行緩慢，嚴重阻礙該方法的進一步普及。

基于區域三維隨機節理網絡的關鍵塊體搜索方法主要有3種：

（1）基于組成自由面上閉合環路節理的圍交判斷的關鍵塊體搜索[2]；

（2）基于有向性和封閉性原理，根據節理上的閉合回路進行全空間內的塊體拓撲搜索[8－9]；

（3）基于節理區域無限切割以及子塊體拼裝的全空間塊體搜索[10]。

上述方法中，方法（1）完全基于關鍵塊體識別出發，認為關鍵塊體的圍閉節理必然在自由面上形成閉合環路，而后根據組成閉合環路節理的交切判斷是否形成塊體，但塊體表面較多時，節理的空間交切非常復雜，且不能適應復雜的研究區域；方法（2）和（3）均能夠勝任全空間的塊體識別，且算法本身具有非常好的健壯性，能夠研究復雜的區域以及實現凹面體的識別，而方法（3）完全由節理切割凸形子塊體，非常易于編程實現。但全空間上的塊體識別無疑需要龐大的計算量。此外，全空間搜索算法也不易根據關鍵塊體搜索特點、工程特點對搜索過程進行適當簡化，避免極少數偶然出現的塊體所帶來巨大的計算量。

本文首先采用方法（1）的閉合環路搜索，然后根據方法（3）利用環路組成節理及相關后緣節理進行空間無限切割再拼裝來識別該環路是否對應關鍵塊體，方法針對性強，能夠很好地適應人為規定的判別條件，如只搜索楔形體以及有無后緣切割等，且能夠實現凹面體關鍵塊體的搜索，此外也繼承了方法（3）編程實現簡單的特點。

2 主要實現過程

首先對本文的關鍵塊體搜索方法的基本步驟介紹如下：

（1）將研究區域劃分成凸形區域組合，如圖 1所示的多平臺邊坡和弧頂型硐室，這樣節理的切割完全轉化成平面切割凸形面和體的計算幾何問題，避免凹形區域切割帶來的麻煩。

（2）根據現場實測統計獲得的結構面產狀、長度與間距的統計特征，采用計算機偽隨機生成三維結構面網絡，并假定結構面為圓盤狀；確定的結構面如 III級以上的結構面可根據其實際空間形態添置進去。此外，三維結構面中心在區域中的生成是采用文獻[2]提出的歸一化方法，可以避免常規均勻分布帶來的節理間距不呈指數分布的情況。

（3）將一個環路的組成節理作為無限平面切割研究區域，判斷該環路對應的空間多面體的每個面是否真正存在，如不存在再進行相關的后緣節理切割。如果是非凸的環路，或者研究區域非凸還需要進行塊體的拼裝，如圖2所示的一個多平臺邊坡，其中一塊體跨平臺，切割后需要進行拼裝得到最終塊體。

（4）最后，對于第（3）步找到的塊體，計算出體積與滑動方向。對于單面滑動和雙面滑動計算其安全系數[11]，如果安全系數小于1則此塊體即為關鍵塊體，而對于墜落式塊體則直接判定為關鍵塊體。

具體的流程如圖3所示。接下來將對關鍵塊體實現過程中的閉合環路搜索和節理無限切割這兩個關鍵環節進行詳細闡述。

圖1 研究區域的凸形分解Fig.1 Convex decomposition of study region

圖2 非凸塊體的拼裝Fig.2 Assembling of non-convex blocks

圖3 關鍵塊體識別流程Fig.3 Flowchart of key block identification

3 閉合環路的搜索

正如前文所述，任一個可滑動的塊體，肯定對應著自由面上跡線圍成的一個閉合環路，且閉合環路的形狀也基本上決定了關鍵塊體的形狀。另外，根據實際工程情況，如若只考慮重力作用，可根據環路最低點以及與該點共棱的塊體內部點的高低位置關系，首先排除一些不能滑出的塊體，也可以非常方便地根據工程經驗對搜索過程作必要的簡化，如工程上經常搜索的楔形體對應的環路為三角形。因此，搜索塊體之前則需要找到自由面上的閉合環路。一般情況下，對于面積較小的環路可不做工程考慮。下面則重點介紹閉合環路的篩選過程。

3.1 無效節理的剔除

對于在自由面上可能構成閉合環路的跡線，必須與兩個或兩個以上的跡線包括自由邊界（如邊坡平臺的眉線）相交，這是一個必要條件，否則該跡線肯定不能構成閉合回路，顯然應被剔除在該閉合環路的第1類節理集合中，該過程即為文獻[6]所述的“砍樹”（tree cut）。上述過程是個循環過程，如圖4所示。圖4(a)為自由面的初始跡線，可以看出與跡線1和2相交的跡線數目分別為0和1，不具備構成閉合環路的必要條件，因此，應當被剔除，而當跡線2被剔除后，與跡線3相交的跡線數目則降為了1條，同樣也需要剔除，依次循環下去，直至所有的跡線都有兩條或兩條以上相交的跡線而終止，最終如圖圖4(b)所示。以上過程，前提是必須已知各跡線間的相交關系。

圖4 無效跡線剔除示意圖Fig.4 Schematic diagram of rejecting invalid traces

3.2 閉合環路的搜索

由于節理數目眾多，不可能無限制地生成復雜的閉合環路，在搜索過程中構成閉合環路的跡線數目需要根據工程經驗以及實際情況設一個上限，以降低問題的復雜性。事實上，根據大量的工程實踐，邊坡的關鍵塊體多以四面體的楔形體為主，即三角形的閉合環路，而洞室圍巖的可動塊體多為四面體、五面體或六面體，有時也會有七面體，但多于七面的可能性較小[2]。

對于邊數較少的閉合環路搜索比較容易實現，例如與一跡線相交的兩條跡線如果也相交的話，則可形成封閉的三角形環路，而如果此兩條跡線有一個共同的相交跡線，則可形成一個封閉的四邊形回路。然而，隨著環路邊數的增加，以上這種尋找閉合環路的編程思路則較難實現。為此，需要一個通用的環路搜索方法。采用邏輯矩陣（元素只能存儲0和1，分別代表不相交與相交）表示跡線相交的信息，是一種圖的存儲結構，如圖5所示的跡線相交圖，其相交信息A見式（1）。

圖5 跡線示意圖Fig.5 Schematic diagram of traces

為了能夠順利進行閉合環路搜索，還需設計一個堆棧管理數據進出，其基本的規則就是越晚得到的數據越早輸出[12]，堆棧的目的則是存儲邏輯矩陣的元素的位置信息，該堆棧的兩個基本函數為壓入（PUSH）和彈出（POP）函數，另外還有重復檢測函數以防止搜索過程中原路返回的現象發生。環路搜索的程序流程見圖6，圖中R表示堆棧中存儲行信息的向量，N為堆棧中數據的組數，從廣義上說，此屬于圖遍歷方法中的深度優先搜索方法。假定搜索閉合環路邊數上限為m，依次從一對角元素(i,i)開始，將該對角元素位置壓入堆棧中形成初始堆棧，而后進行條件循環。在每次循環中彈出一個位置數據 (j0,k0)，如果 j0=i且 A(j0,k0)=1那么堆棧中行位置元素加上i即形成一個環路。接下來，找到該位置同列中的下一個非0元素位置 (j1,k0)，將這一數據壓入到堆棧中，判斷 j1和堆棧中數據的重復性，如果不重復且堆棧個數據未數達到m-1，壓入(i,j1)，如果達到m-1，則不壓入(i,j1)，但如果A(i,j1)= 1，同樣也形成一個閉合環路。之后，繼續進行重復彈棧操作。當棧中元素為空時，壓入對角元素(i+1,i+1)，如此直至對角元素循環結束。

圖5所示的跡線相交圖，經過上述的算法可以得到的閉合環路有[2,4,1]、[2,4,5,3,1]、[2,5,3,1]、[2,5,4,1]、[3,5,2,4,1]、[3,5,4,1]和[4,5,2]這 7 個結果，這其中[2,5,4,1]這種交叉型的互連不能夠成閉合環路需要對其排除。一般情況下，研究塊體的閉合環路為凸性。否則造成其安全系數難以求解，即使如圖2(d)所示的非凸塊體，其自由面上的閉合環路仍為凸的。因此，根據環路的對應交點是否均在每個環路節理的同側可排除掉交叉型和非凸性的回路，即排除了[2,5,4,1]和[3,5,2,4,1]，從而剩下 5個閉合回路。

圖6 跡線閉合回路搜索流程圖Fig.6 Flowchart of searching closed loop

4 基于節理無限切割的關鍵塊體識別

已知了自由面上的閉合環路，則以閉合環路為單位進行關鍵塊體搜索。

4.1 空間多面體的數據結構

如圖2所示，對于一個閉合回路，確定其可能形成的最大塊體，在該過程中必須有規則的數據鏈表和拓撲意義上的數據結構。一個完整的多面體結構為一樹狀的結構（如圖7所示），其層次為體-面-環-棱-點，除點含有坐標表示幾何信息外，其余均為拓撲信息。

圖7 多面體的邏輯表示Fig.7 Logic representation of polyhedron

多面體的存儲采用類半邊數據結構[13]，但不使用“半邊”層，而是在環到邊的拓撲信息中采用一個邏輯數記錄該邊棱的方向（如圖8所示），共有5個雙向的鏈表來依次表示體、面、環、邊以及點這5層幾何或拓撲信息，對于任一層上的某節點，均能通過指針訪問其上下層拓撲關系，也可以訪問該層內其前一個與后一個指針，如采用FORTRAN95定義環（LOOP）結構體如下

TYPE :: loop_type

integer :: id

type(includeE_type),pointer :: includeE

type(face_type),pointer :: belongF_1 !面指針

type(face_type),pointer :: belongF_2 !面指針

type(loop_type),pointer :: prev

type(loop_type),pointer :: next

END TYPE loop_type

這其中有

TYPE :: includeE_type

type(edge_type),pointer :: edge !邊指針

logical :: dir !邊方向

type(includeE_type),pointer :: prev

type(includeE_type),pointer :: next

END TYPE includeE_type

可以看出，一個環節點即含有了環到邊的拓撲信息，亦含有環屬于面的拓撲信息，其他各層節點亦是如此。另外，includeF、includeL和includeE分別表示體、面和環的組成，亦是一個雙向鏈表，includeF節點包含著指向面的指針和面的外法向信息；includeL節點則是指向環的指針，且該鏈表的頭節點定為“外部”邊界，其余節點均為面上的“孔”；includeE的節點如上FORTRAN代碼定義所示，包括指向邊的指針和判別邊方向的邏輯變量。

圖8 多面體的層次結構與數據結構Fig.8 Hierarchy and data structure of polyhedron

有了上述完整地表征多面體幾何和拓撲信息的數據結構，則很方便進行下面的切割與拼裝的計算幾何操作。

4.2 塊體的切割與拼裝

規定空間內平面的正法向為：當平面法向與 z軸夾角小于5°時，該法向與 z軸正方向夾角為銳角；否則，當其與x軸夾角小于45°，那么該法向與x軸正方向夾角為銳角，如仍不成立，則以與y軸正方向夾角為銳角的方向為正。已知空間平面P方程為

式中：(a,b,c)為該平面的單位正法向向量；d為系數。對于任一點 v(X,Y,Z)，則點與平面的位置關系L為

式中：ε為個小正數。L=0、1和-1分別表示在面上、面的正側和負側。

若邊、平面以及多面體的所有頂點和平面的位置關系中既存在1亦存在-1，則很容易判斷其與該平面相交。圖9為多面體與節理的相交示意圖，由式（3）可以判斷出節理B與C所在的平面肯定與多面體不相交。此外，對于節理D與多面體是沒有任何關系的，則應避免節理切割該塊體，但當塊體數量比較龐大時，如果嚴格地判斷塊體和節理是否交切的自身計算量亦會非常大，因此，需要較為經濟的方法判斷兩者是否呈切割關系。文獻[10]通過節理與多面體的范圍箱是否有重疊來決定兩者是否進行相交計算。范圍箱指的是一個長、寬、高分別平行于空間坐標軸的長方體，使節理或多面體恰好落在該長方體內。對于多面體其范圍箱可根據所有頂點的坐標最大、最小值求出；而對于節理則需根據其在各面上的投影求得。有了范圍箱后，只需進行大小判斷即可斷定范圍箱是否有重疊，從而決定是否進行平面與多面體的切割計算。有了該步判斷，圖9的節理D則不會進行切割計算。

圖9 多面體的與節理交切關系示意圖Fig.9 Schematic diagram of cutting relationship between polyhedron and joints

下面分步驟介紹平面切割塊體：

（1）平面切割邊。該子程序將返回一個頂點指針。當邊與平面相交時，在頂點鏈表的尾部增加一個節點，并更新邊到環以及點到邊的拓撲信息，將返回指針指向該新頂點；若平面恰與邊交與邊的終端點，無需更新拓撲信息，直接返回終端點指針；若以上都不滿足則返回空指針。

（2）平面切割環。該子程序返回一個邊的指針，當環的某一邊在平面上時，直接返回該邊的指針。圖10為一平面切割環E12→E23→E34→E41，對環所有的邊進行第（1）步操作后，環變成E15→E52→E23→E36→E64→E41，并根據第（1）步返回的兩個頂點指針在邊鏈表尾部生成一個新邊 E56（或E65）。接下來，配置兩個環指針，對環所有的邊進行循環，并以一邏輯變量為指示，若為真將該邊原封不動地放入環一，反之則放入環二，若當前邊的終點為5或6的一點時，將對應E56或E65放入當前環的邊鏈表中，并對邏輯變量取反，根據該循環可得到兩個環 E15→E56→E64→E41和 E52→E23→E36→E65，用兩個新環中一個覆蓋原環，另一個添加到環鏈的尾部，最后再將新邊E56作為指針返回。

圖10 平面切割環Fig.10 Loops cutting by plane

（3）平面切割體。該程序返回一個面的指針。圖 11是一體被平面切割示例，首先是體的每個面（環）進行第（2）步操作后，根據返回的所有邊指針，生成一個新的環和面；配置兩個體指針，體一和體二分面位于平面的正負方向，將所有的面按照與平面的位置關系分別放置至兩個體中，并保持外法向不變，再將新面添加到兩個體中，其外法向分別對應著該平面的負法向和正法向；最后將平面正方向的體替換原體，而將負方向上的體放置在體鏈表的尾部。

圖11 平面切割體Fig.11 Solids cutting by plane

有了切割后的子塊體，可以將節理恢復成圓盤狀，從而判斷子塊體是否有效，如該子塊體含有研究區域邊界面，則同樣也是無效的。

首先，由形成自由環路的節理切割研究區域，含有閉合環路的塊體則稱為主塊體。若主塊體是有效的，則無需進行后緣切割；否則，則采用與閉合環路相交的第2類節理（與第1類節理相交的節理）對主塊體進行切割，通常為了減少復雜程度，不考慮第3類節理（與第2類節理相交的節理）。有了一系列將主塊體切割后的子塊體，則可進行子塊體的拼裝。

文獻[10]給出了通用的合并方法，拼裝的基礎是共面刪除邊、共邊刪除點的拓撲信息改變過程，在此不做過多贅述。主塊體被切割后的子塊體進行拼裝后，若不含有研究區域邊界，則此塊體則為該閉合環路對應的塊體。

5 算 例

本節以某掛幫礦開采的礦柱為例搜索一頂柱的關鍵塊體。

該區域的節理層面信息如表1所示，此外區域的層面產狀為345°∠73°，層面的間距為1.21 m。頂柱的尺寸為15 m×40 m×12 m，頂柱的幾何模型及三維節理網絡如圖12所示。由于頂柱的塊體基本為墜落式，因此，不需要根據摩擦角和凝聚力計算其安全系數。設H為塊體的頂點距自由面的最大距離，根據上述關鍵塊體搜索過程獲得H＞0.2 m關鍵塊體如圖13所示。

該頂柱的關鍵塊體以四面體（106個）為主，五面體（31個）次之，亦有少量的六面體（20個），無七面體，主要是節理1、2、4和層面相交生成四面體與五面體，節理3近水平從而對塊體起到后緣切割作用。頂柱的關鍵塊體最大體積為2.739 m3，開采過程可表現出小規模的塌方。

表1 節理產狀統計Table1 Statistical parameters of joint occurrence

圖12 頂柱的三維節理網絡Fig.12 3D joint network in the region of top pillar

圖13 頂柱關鍵塊體的分布圖Fig.13 Distribution of key blocks in top pillar

6 結 論

本文針對關鍵塊體搜索的特點，將研究區域分解成凸子區域，搜索自由面上的閉合環路，然后利用環路組成節理及后緣相交節理進行無限切割和拼裝來識別該環路是否對應關鍵塊體，得到如下結論：

（1）將自由面節理跡線相交信息用圖表示，利用堆棧進行圖的遍歷操作，可以適應任意邊數目的閉合環路搜索。

（2）采用類半邊數據結構存儲空間多面體，除點鏈含有幾何信息外，其余均為拓撲信息，拓撲結構明確，在節理平面切割邊、環和面時能夠自由的更新切割后的拓撲關系。

（3）結合自由面上的閉合環路搜索與節理平面無限切割的關鍵塊體搜索方法，避免了采用環路節理圍閉在節理較多時復雜的判斷以及全空間塊體搜索塊體龐大的計算量，方法針對性強，較好地適應人為規定的判別條件，可使關鍵塊體搜索達到效率與效果的平衡。

[1]蔡美峰,何滿潮,劉東燕. 巖石力學與工程[M]. 北京:科學出版社,2002.

[2]賈洪彪,唐輝明,劉佑榮,等. 巖體結構面三維網絡模擬與工程應用[M]. 北京: 科學出版社,2008.

[3]GOODMAN R E,SHI G H. Geology and rock slope stability—Application of the key block concept for rock slopes[C]//Proceedings of 3rdInternational Conference on Stability in Surface Mining. New York: AIME/SME,1982.

[4]SHI G H,GOODMAN R E. The key blocks of unrolled joint traces in developed maps of tunnel walls[J].International Journal for Numerical Analytical Methods in Geomechanics,1989,13(2): 131－158.

[5]鄔愛清,任放,郭玉. 節理巖體開挖面上塊體隨機分布及錨固方式研究[J]. 長江科學院院報,1991,8(4): 27－34.WU Ai-qing,REN Fang,GUO Yu. A study of block random distribution and anchor design for joint rock excavation[J]. Journal of Yangtze River Scientific Research Institute,1991,8(4): 27－34.

[6]張奇華. 塊體理論的應用基礎研究與軟件開發[D]. 武漢: 武漢大學,2004.

[7]LU Jun. Systematic identification of polyhedral rock blocks with arbitrary joints and faults[J]. Computers and Geotechnics,2002,29(1): 49－72.

[8]SHI Gen-hua. Producing joint polygons,cutting joint blocks and finding key blocks from general free surfaces[J]. Chinese Journal of Rock Mechanics and Engineering,2006,25(11): 2161－2170.石根華. 一般自由面上的多面節理生成、節理塊切割與關鍵塊搜尋方法[J]. 巖石力學與工程學報,2006,25(11): 2161－2170.

[9]張奇華,鄔愛清. 邊坡及洞室巖體的全空間塊體拓撲搜索研究[J]. 巖石力學與工程學報,2008,27(10): 2072－2078.ZHANG Qi-hua,WU Ai-qing. Study of spatial block topological identification of slope and cavern rock mass[J]. Chinese Journal of Rock Mechanics and Engineering,2008,27(10): 2072－2078.

[10]YU Q,OHNISHI Y,XUE G,et al. A generalized procedure to identify three-dimensional rock blocks around complex excavations[J]. International Journal for Numerical and Analytical Methods in Geomechanics,2009,33(3): 355－375.

[11]于青春,薛果夫,陳德基. 裂隙巖體一般塊體理論[M].北京: 中國水利水電出版社,2007.

[12]劉大有,唐海鷹,孫舒楊,等. 數據結構[M]. 北京: 高等教育出版社,2001.

[13]伏玉琛,周洞汝. 計算機圖形學——原理、方法與應用[M]. 武漢: 華中科技大學出版社,2003.