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

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

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

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

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

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

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

1 3D數(shù)字化結(jié)構(gòu)面跡線模型

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

將RTK-GPS量測數(shù)據(jù)導(dǎo)出后經(jīng)過必要整理，導(dǎo)入任何三維建模軟件，經(jīng)上述過程即可建立3D數(shù)字化結(jié)構(gòu)面模型，如圖1所示。

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

圖1 3D數(shù)字化結(jié)構(gòu)面模型

2 基于地形化數(shù)字測窗的跡線面密度估算結(jié)果

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

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

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

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

(1)

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

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

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

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

2.2 矩形測窗估算裂隙密度

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

(2)

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

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

數(shù)字為測窗布置編號

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

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

2.3 測區(qū)結(jié)構(gòu)面分組及基本參數(shù)

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

表1 結(jié)構(gòu)面分組及參數(shù)

2.4 面密度估算結(jié)果及分析

4組結(jié)構(gòu)面對應(yīng)于圖2測窗布置情況的面密度的估算結(jié)果如圖3所示。

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

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

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

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

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

3 分塊體密度優(yōu)化及3D結(jié)構(gòu)面網(wǎng)絡(luò)模型

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

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

當(dāng)半徑服從負(fù)指數(shù)分布時，面密度向體密度轉(zhuǎn)化計算公式[8-9]為

(3)

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

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

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

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

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

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

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

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

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

圖7 按優(yōu)化體密度生成的3D結(jié)構(gòu)面離散網(wǎng)絡(luò)模型

4 結(jié)論

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

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

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

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