巖體結構面面密度統計及體密度優(yōu)化

周愛蓮，章楊松*，李曉昭，石杏喜

(1.南京理工大學土木工程系，南京 210094；2.南京大學地球科學與工程學院，南京 210093)

巖體裂隙特征是工程巖體質量評價和穩(wěn)定性分析及裂隙巖體滲流分析的重要基礎。密度是表征巖體裂隙特征的重要參數之一。傳統方法是在現場布置人工測線或測窗采用手工量測記錄裂隙特征，通過有關估算公式得出一定產狀露頭面上的線密度或面密度。用于3D裂隙網絡建模的體密度不能直接通過測量獲取，一般是通過線密度或面密度計算并經檢驗修正來確定。

測窗法常用在露頭面的比較分析[1]，但矩形測窗測量結果可能受測窗長短邊方向偏差的影響。Mauldon[2]和Zhang等[3]提出了用有限大小圓形測窗估算真實跡長分布(平均值和標準差)和面密度的方法。然而，在現場布置測窗進行統計，還存在著諸多困難。一些現代量測技術，如RTK-GPS(real time kinematic-global positioning system)技術[4]和數字攝影測量[5]等明顯地改進了測量巖體結構面特性的能力，能方便獲取準確的結構面網絡，為統計密度提供了可能。郭亮等[6]利用RTK-GPS測量技術獲取甘肅北山結構面信息，對參數數字化統計后建立模型，為結構面網絡模擬提供了新的思路。趙佳斌等[7]和Xu等[8]根據數字攝影測量技術得到的三維點云數據，提出了一種結構面自動識別方法。測窗法[9]已直接用于地面或航空照片和衛(wèi)星影像上裂隙的分析。

由于地形往往為非平面狀的，在數字化模型上用平面狀的測窗統計，仍存在一定誤差。Sturzenegger等[10]提出了一種稱為“擬合地形或地形化”的圓形測窗。該方法適應了巖石斜坡起伏不平的變化。用隨地形變化的測窗允許考慮在平面狀取樣窗兩邊一定厚度范圍的帶內的不連續(xù)面的跡線。這種方法被認為是當天然或人工露頭不是完全平整，并且如果僅考慮與平面狀的測窗交切的不連續(xù)面作為結果而產生測量偏差的一種更符合實際的選擇。

在建立裂隙網絡三維模型時，結構面的體密度(確定了模擬區(qū)域內結構面數量)是重要參數之一[11]。由于現有轉化計算公式的適用條件和問題的復雜性，采用計算公式[12-13]確定的體密度只能作為初始值，需要進行校核并不斷調整，而當模擬區(qū)域較大時，在整個區(qū)域采用同一體密度是不符合實際的。王述紅等[14]采用了根據結構面發(fā)育程度劃分分區(qū)進行結構面三維網絡模擬，并采用擴大模擬區(qū)域來考慮結構面中心在周圍區(qū)域對模擬區(qū)的影響，同時采用結構面平均間距和跡長進行動態(tài)校核的方法。

基于此，現提出了基于差分進化算法對模型分區(qū)體密度進行優(yōu)化的結構面三維網絡建模方法，并研究了結構面跡線3D數字化模型及數字化測窗建立方法，以期建立與實際更加符合的結構面網絡模型，為工程建設提供更加可靠的依據。

1 3D數字化結構面跡線模型

跡線圖能夠記錄描述結構面交切模式，終止關系和其他空間集聚特性。但采用傳統方法獲取結構面空間跡線圖不易且費時，RTK-GPS技術[15]是采用載波相位動態(tài)實時差分(real time kinematic,RTK)技術與全球定位系統(global positioning system,GPS)技術相結合，測得的是結構面控制點的空間坐標，可達cm級精度，將屬于同一條結構面的點連起來即可方便得到較準確的結構面空間跡線圖。同時利用所測得的點云，經過三角形網格化，及必要的插值加密，可以得到近似的露頭面模型。為便于后面的應用，還可在三角形網格化及插值的基礎上將露頭面用矩形擬合(實際為曲面矩形)，即數字化面模型(digital surface model, DSM)。將空間跡線與擬合的露頭面疊加在一起，得到完整的數字化結構面跡線模型(digital discontinuity trace model, DDTM)。

將RTK-GPS量測數據導出后經過必要整理，導入任何三維建模軟件，經上述過程即可建立3D數字化結構面模型，如圖1所示。

圖1為甘肅北山笈笈槽的一個測區(qū)，其近似為一矩形, 其平面尺寸(XY平面內)面積近似為50×110=5 500 m2，X軸正方向指向東(E),Y軸正方向指向北(N)，Z軸為高程方向，向上高程增加。采用模糊聚類方法將裂隙分成4組。從裂隙的分布可看出，該區(qū)第二組(走向北東)結構面發(fā)育,跡線相對較長，而且在一定程度上對其他各組裂隙有限制作用。其他三組結構面跡線均較短,且分布不均一。從圖1中可看出，數字化結構面跡線與數字化地面貼合較好，在該模型上布置測窗進行統計分析是可行的。

圖1 3D數字化結構面模型

2 基于地形化數字測窗的跡線面密度估算結果

在數字化地形模型及跡線模型和相應數據庫建立后，可用于進一步分析研究，例如測窗取樣和測量。本研究中在3D模型上同時使用正方形和圓形測窗進行跡線密度估算。為比較測窗大小對統計結果的影響，也采用了幾種不同半徑的圓形測窗。

2.1 圓形測窗估算裂隙跡線密度

圓形測窗以跡長分布—自由端點建立估算式，可以避免長度偏差，刪節(jié)偏差。這種方法用到與圓相交的跡線數(N)，兩端點都被測窗刪剪(交切)的跡線數(N0)和兩端點都在窗內觀察到(包容)的跡線數(N2)。

Mauldon[2]進一步發(fā)展了這一方法，用圓形測線和測窗來估算跡線參數。這種方法基于測窗內跡線端點計數m和測窗面積進行計算：在圓形測窗內的跡線端點數, 即兩端點都在窗內觀察到的跡線端點數為2, 一端在測窗內的跡線端點數為1, 跡線密度(λa)的計算公式為

(1)

式(1)中：λa為跡線密度；r為測窗半徑；m為測窗內跡線端點數。

由于在現場很難布置理想的平面狀測窗，因此首先選擇相對比較平整的露頭面，再采用適當的方法在該露頭上圈出測窗范圍，取樣時實際上將該測窗范圍內在露頭上出露的某組裂隙跡線都考慮在內，而并非是真正采用平面狀的測窗來進行取樣。采用數字化模型時，盡管可布置平面狀測窗，但是數字化露頭面本身也并非完全平整的，跡線在空間是非直線狀的，所以平面狀測窗無法與露頭面或測窗范圍內所有結構面跡線吻合，盡管通過調整模型視向并將跡線投影到測窗上進行統計，可以近似得到結果，但這樣會產生較大誤差，并且不方便自動計算。

“擬合地形或地形化”的圓形測窗計算時考慮在平面狀取樣窗兩邊一定厚度范圍的帶內不連續(xù)面的跡線[10]。為了應用最原始的適合平面狀圓形測窗而提出的公式, 假定在地形化的測窗內終止和交切的跡線能夠投影到一局部圓盤上, 該圓盤的面積和周長是等效于包含地形化測窗的圓柱的底面。

根據以上思想，反過來只要采用垂直于局部地形的圓柱與地形面交切，即得到地形化圓形測窗，如圖2所示。圖2中顯示第4組跡線, 其中正方形測窗以顏色區(qū)分，圓形測窗僅顯示半徑為4 m的情況，圖2中測窗并非平面狀，而是隨地形變化的。同時由于式(1)不涉及測窗內跡線長度，只與跡線和圓形測窗的幾種幾何關系有關，所以如果將整個3D模型和“地形化”測窗投影到合適的平面內時(緩地形投影到水平面，陡露頭則投影到某一擬合平面)，以上關系(即m計數)保持不變，這就給統計分析帶來方便。

2.2 矩形測窗估算裂隙密度

目前基于矩形窗口估算面密度的計算公式有多種，代表性的有：一種是Kulatilake等[12]提出的矩形測窗面密度估算方法，該公式計算復雜，需要獲取研究區(qū)域跡長的概率密度函數。另一種是Mauldon[2]估算方法，其估算面密度方法基于凸多邊形，圓形矩形都適用。Mauldon[2]給出的矩形測窗跡線中點面密度估算公式為

(2)

式(2)中：λa為跡線中點面密度；N1為測窗內一端可見跡線的條數；N2為未刪節(jié)(兩端可見)跡線的條數；W為矩形測窗的寬度；H為該測窗的高度。

在3D數字模型上“地形化”圓形測窗的實現，需要相對繁雜的幾何操作或運算才能完成，而地形化矩形測窗則容易實現，可以直接采用多個擬合地面的小矩形曲面組合得到，如圖2所示。圖2中不同矩形測窗以顏色加以區(qū)分，如果將數字化模型視向調整到合適的平面內，則容易統計出各矩形測窗內的N1和N2。式(2)中的W和H則用每個矩形測窗的4個角點坐標計算得到。為了與圓形測窗對應，采用近似正方形的矩形測窗，這一方面使得同一位置的正方形測窗和圓形測窗覆蓋的裂隙跡線能最大程度上相似；另一方面，正方形測窗也能在一定程度上減小一般矩形測窗長短邊不同方向布置造成的方向偏差。

數字為測窗布置編號

在估算面密度時，關于單個測窗尺寸的確定，需考慮兩個方面因素：①由于采用Mauldon方法估算，為盡可能地減少N0型(兩端點都被測窗刪剪的跡線數)跡線對估算結果的影響，需選取較大測窗以使得N0類跡線數目盡可能少；②單個測窗面積越小，就越能精確描述研究區(qū)域內局部面密度情況。因此，測窗尺寸的選擇需平衡這兩方面因素。本例布置在研究區(qū)域各近似正方形測窗的尺寸相同，面積均為44.3 m2，圓形測窗分別取半徑3 m、3.75 m與正方形測(窗等面積)和4 m三種。

測窗產狀是計算和校核體密度所必需的參數之一，矩形測窗的產狀也更容易得到。在數字化模型上自動批量提取地形化正方形測窗的4個角點的空間坐標，采用這四點坐標用最小二乘法擬合得到局部測窗的產狀。

2.3 測區(qū)結構面分組及基本參數

考慮到用測窗法估算跡長影響因素較多[16]，對低密度裂隙巖體，由于不易布置滿足有關條件的合適大小的測窗，所以經常難以得到較確切的跡長估計數據，而本文研究所在場地露頭條件優(yōu)越，采用RTK-GPS量測方法能得到較完整的全跡長數據，所以裂隙跡長可以根據測量跡線的端點計算，則單條裂隙跡長的計算直接求笛卡爾坐標系下的首尾兩點距離即可。采用模糊聚類方法，得到裂隙分組結果有關參數如表1所示。

表1 結構面分組及參數

2.4 面密度估算結果及分析

4組結構面對應于圖2測窗布置情況的面密度的估算結果如圖3所示。

根據估算過程和結果可知，估算的面密度大小除與該測窗內實際裂隙疏密有關，實際上還與測窗范圍大小有關。范留明等[17]在采用密度劃分均質區(qū)時，認為如果步長(即范圍)選擇過大，則會丟失較多密度細節(jié)，甚至損失一些重要信息。實際上大范圍的密度是該范圍內的平均值，范圍越小，越能反映局部結構面的發(fā)育情況。如圖3所示，采用半徑為3 m的圓形測窗，各組結構面曲線的波峰明顯高于同位置半徑4 m和3.75 m圓形測窗及正方形測窗的對應面密度值，而波谷則是小范圍測窗得到的，密度更小(密度為零的測窗除外)，所以測窗范圍偏小時，從整體而言，存在過高(跡線分布密度的位置)或過低估計(密度低時)面密度。不同于跡長估計，密度是一個相對值，一旦確定好合適的測窗大小，在后續(xù)建模時也采用相應的子區(qū)大小即可。

表2 正方形測窗和與等面積的圓形測窗面密度平均值的比較

圖3 不同測窗面密度統計分析(第1組)

可見相對誤差都很小。其他各測區(qū)這兩種測窗的統計結果表明，無論是對應位置的兩種測窗密度單值還是按平均密度值考察，均具有與本例相似的結果。可以認為按正方形測窗可以得到圓形測窗幾乎相同的面密度結果。而正方形測窗也有優(yōu)點，如可以完全連續(xù)覆蓋整個測區(qū)，局部測窗的產狀更方便確定。而一旦正方形測窗確定以后，按其中心位置和產狀即可確定圓形數字化測窗參數，并可在同一位置布置不同半徑的圓形測窗。可同時采用兩種測窗相互補充、相互印證。

從計算結果還可發(fā)現，部分測窗面密度為零，這是相對低密度巖體的一種空穴效應。然而實際上如果隨著測窗范圍增大，由于平均化作用，則會掩蓋這種空穴現象。由此還可以推測，如果在現場統計量測工作中，僅選擇結構面出露好的露頭布置測窗量測統計，會存在過高估計面密度的可能性。

3 分塊體密度優(yōu)化及3D結構面網絡模型

用于3D結構面網絡建模的體密度是不能直接測定的，一般根據一定假設條件下推導的轉化公式計算得到。由計算參數直接建立的3D結構面網絡模型的面密度與對應的實測值往往有一定的差異，需要校核和修正才能提高模型精度。以往的研究多采用單一切面(一個產狀的測窗)和單一參數(平均面密度)進行校核。本研究嘗試采用差分進化算法對3D結構面網絡模型進行分塊體密度優(yōu)化。

差分進化算法(differential evolution，DE)是Storn等[18]在1997年為求解切比雪夫多項式首次提出。DE算法進化流程主要分為：變異、交叉和選擇三個步驟。該算法相對簡單，易于實現，具有很好的收斂性，是一種有效的全局優(yōu)化技術。基于差分進化分塊體密度優(yōu)化，其基本思想是尋找各分塊體密度，使與露頭面對應產狀的模型上各個測窗結構面的面密度與實測面密度之差(采用絕對值)小于預先設定的某一比較小的限制值。

當半徑服從負指數分布時，面密度向體密度轉化計算公式[8-9]為

(3)

式(3)中：E(λv)為結構面體密度期望值；E(λa)為該組結構面面密度期望值；E(D)為該組結構面直徑的期望值；v為結構面平均走向與測窗所在面的夾角。當半徑服從其他分布形式時不能采用式(3)，本研究區(qū)結構面跡長近似服從對數正態(tài)分布，假設半徑分布與跡長分布具有相同的形式，所以不能采用式(3)。而且即使分布滿足條件，采用式(3)計算的體密度仍需按實測面密度進行校核和修正。

對于大范圍分塊組合形成的三維結構面網絡模型，在進行面密度校核及體密度調整過程中估算面密度時，需考慮相鄰分塊之間的相互影響，當分塊數量多時, 其過程比較繁瑣復雜，人工調整不易得到最佳結果。而要使模型中多個測窗與實際測窗面密度差異最小，可歸為多目標優(yōu)化問題。

為在三維結構面建模過程中考慮多個連續(xù)測窗的實測面密度校核，分塊體密度優(yōu)化并生成結構面模型的步驟如圖4所示。其中幾點強調如下：①首先對模型區(qū)按面密度統計分塊相同的XY坐標進行模擬分區(qū)，暫時假設不考慮豎直方向的體密度變化；②各分塊初始體密度按式(3)計算，各分塊體密度范圍取相同值，下限0，上限在所有分塊初始體密度最大值基礎上適當放大；③在模型每個分塊位置以和實測相同產狀及大小的測窗(平面狀正方形和圓形測窗)進行檢驗計算，將模型面密度和實測面密度進行對比，對比時采用的是正方形測窗和等面積圓形測窗所得面密度的平均值；④采用差分進化算法不斷調整各分塊結構面體密度，直到模型各測窗面密度和實測面密度之間的誤差滿足事先設定的目標值，取|λi模型-λi實測|為目標函數, 本例設定目標值為0.001個/m2。

圖4 考慮體密度優(yōu)化的結構面3D模型建模流程

由于結構面面跡線網絡不僅與其數量還與跡長等有關，故將體密度差分進化優(yōu)化求解嵌入三維隨機結構面生成程序, 與模型其他裂隙參數一起產生，這樣得出的結構面模型數據更合理。所有過程編程由計算機完成。

第1組各分塊初始體密度及優(yōu)化體密度的對比，如圖5所示。結果顯示最終采用的體密度與按式(3)計算的體密度在某些測窗存在很大差異，而且各分塊不能采用統一的修正系數來進行修正。第1組模型中各測窗的面密度及其與對應實測面密度的誤差，如圖6所示。計算結果表明, 4組結構面優(yōu)化后所有測窗中的面密度最大誤差為0.000 919個/m2, 小于設定的0.001個/m2, 對應測窗的面密度誤差相對值為1.016 7%。可見經過優(yōu)化得到的模型體密度精度高。對于大量連續(xù)測窗相互影響的復雜情況，如果靠人工進行校核和調整既費時，而且?guī)缀醪豢赡芡瑫r達到這樣高的精度。

圖5 初始體密度及優(yōu)化體密度的對比(第1組)

對分塊體密度優(yōu)化后，按圖4流程生成的結構面數據建立的三維結構面網絡模型，如圖7所示。這樣建立的模型能反映實測面密度隨位置的變化(圖6)，特別是局部空穴效應(體密度為零)的存在。顯然這種分塊組合模型與整個模型區(qū)內按平均體密度建立的模型相比較，在后續(xù)的巖體體穩(wěn)定性及滲透性評價等方面的結果無疑均會有所不同。

圖6 優(yōu)化體密度3D模型的面密度值及其誤差(第1組)

圖7 按優(yōu)化體密度生成的3D結構面離散網絡模型

4 結論

(1)采用地形化圓形和正方形測窗，并與跡線端點計數算法相結合，減小了由于應用平面狀測窗在起伏不平的3D模型上統計計算的困難和可能的偏差，但同一測窗內地形起伏不宜過大。結合數字化方法，確定結構面面密度的精度和效率明顯高于傳統現場量測統計方法。其方法也適用于攝影測量或三維激光掃描等數字化方法獲取的結構面的密度統計。

(2)測窗法統計確定結構面面密度時關于測窗大小的選取不同于估算跡長時測窗大小的選取。在大的測窗范圍統計時是對測窗內的密度進行了平均化。建議在滿足N0型跡線出現情況的條件下，采用較小的測窗尺寸統計面密度，這樣能更好地反映局部面密度情況。在后續(xù)3D結構面網絡建模應采用同樣的分塊尺寸。

(3)普遍認為圓形測窗優(yōu)于矩形測窗，因為圓形測窗消除了方向偏差。研究表明, 地形化正方形測窗統計計算的面密度與等面積圓形測窗統計的結果非常接近，且與圓形測窗相比，正方形測窗可以完全連續(xù)覆蓋整個測區(qū)，沒有空隙。因此正方形測窗更加合理，可用正方形測窗代替圓形測窗。在數字化模型上由于可實現自動計算，可同時采用這兩種測窗，充分利用兩者各自的優(yōu)點，相互補充及校核。

(4)將巖體分區(qū)與體密度優(yōu)化方法相結合，通過編程實現對體密度的自動優(yōu)化。將差分進化算法優(yōu)化體密度的程序嵌入3D結構面模型建模程序，同時實現了自動對大量連續(xù)分布測窗的面密度校核和對分區(qū)體密度進行調整優(yōu)化，并同步生成與優(yōu)化體密度相對應的結構面模型數據。這樣生成的模型中結構面的空間分布更符合實際情況，為后續(xù)巖體穩(wěn)定性及滲流分析奠定了基礎。