簡(jiǎn)介面向地質(zhì)建模與GIS領(lǐng)域的工程師和研究者一套圍繞巖體結(jié)構(gòu)面自動(dòng)提取與產(chǎn)狀計(jì)算的代碼包覆蓋點(diǎn)云預(yù)處理濾波、降噪、邊緣檢測(cè)、連通性分析以及TIN構(gòu)建、密度分析等關(guān)鍵技術(shù)環(huán)節(jié)可根據(jù)點(diǎn)云密度自動(dòng)識(shí)別結(jié)構(gòu)面并輸出傾角、走向等產(chǎn)狀參數(shù)適用于礦產(chǎn)勘查、隧道掘進(jìn)、地質(zhì)災(zāi)害評(píng)估等場(chǎng)景。壓縮包共713個(gè)文件以710個(gè)txt數(shù)據(jù)文件含點(diǎn)云法向量、三角網(wǎng)、重采樣、密度、距離等中間結(jié)果為主另附1個(gè)Python核心腳本和2個(gè)zip數(shù)據(jù)包整體29.39MB目錄與文件名對(duì)應(yīng)清晰便于按流程逐模塊對(duì)照學(xué)習(xí)。已有838人學(xué)習(xí)下載既能幫助初學(xué)者理解從點(diǎn)云到結(jié)構(gòu)面產(chǎn)狀的完整算法邏輯也可供從業(yè)者直接修改腳本并應(yīng)用于自有數(shù)據(jù)省去大量手工標(biāo)注和參數(shù)試錯(cuò)的時(shí)間提升分析效率與結(jié)果一致性尤其適合需要快速驗(yàn)證地質(zhì)結(jié)構(gòu)面識(shí)別方案的項(xiàng)目前期階段。1. 巖體結(jié)構(gòu)面自動(dòng)提取為什么三維點(diǎn)云不能直接算產(chǎn)狀三維激光掃描儀能在幾分鐘內(nèi)把一面高邊坡變成上百萬(wàn)個(gè)三維點(diǎn)但這些點(diǎn)并不區(qū)分哪個(gè)點(diǎn)屬于哪個(gè)結(jié)構(gòu)面。點(diǎn)云里沒(méi)有“面”的實(shí)體、沒(méi)有邊界、沒(méi)有ID——所謂結(jié)構(gòu)面只是取自一個(gè)斷裂面上的一簇大致共面的點(diǎn)。人工在CloudCompare里框選節(jié)理面再擬合產(chǎn)狀單個(gè)面通常要花十分鐘以上一個(gè)邊坡幾十個(gè)結(jié)構(gòu)面人工作業(yè)需要一到兩天。而這個(gè)Python處理管線(xiàn)Automatic extraction of discontinuities.py做的事就是自動(dòng)從點(diǎn)云中找出屬于同一個(gè)不連續(xù)面的點(diǎn)分組并計(jì)算產(chǎn)狀。它依賴(lài)的關(guān)鍵數(shù)據(jù)包括法向量、TIN三角網(wǎng)、重采樣點(diǎn)云、密度場(chǎng)和距離場(chǎng)。適合做邊坡勘察、隧道超前地質(zhì)預(yù)報(bào)、礦山邊坡穩(wěn)定性分析的巖土工程師。2. 結(jié)構(gòu)面提取數(shù)據(jù)管線(xiàn)法向量、TIN與密度場(chǎng)的配合邏輯2.1 點(diǎn)云到TIN為什么選三角網(wǎng)而不是體素柵格原始點(diǎn)云沒(méi)有拓?fù)潢P(guān)系點(diǎn)與點(diǎn)之間的鄰接無(wú)法直接定義。要計(jì)算法向量、曲率或者做區(qū)域生長(zhǎng)必須先建立鄰域關(guān)系。常見(jiàn)做法是用kd-tree做近鄰查詢(xún)但在巖體表面這種幾何復(fù)雜場(chǎng)景下單純用k近鄰會(huì)把互不相連的懸空點(diǎn)拉進(jìn)鄰域尤其是在結(jié)構(gòu)面邊緣和陡坎位置。TIN不規(guī)則三角網(wǎng)用Delaunay三角化把點(diǎn)云連接成連續(xù)的三角形網(wǎng)格每個(gè)頂點(diǎn)可以直接通過(guò)三角形的邊找到真實(shí)鄰接頂點(diǎn)這種拓?fù)潢P(guān)系更接近地質(zhì)表面的實(shí)際連通性。為什么不選體素柵格體素化把空間離散成規(guī)則小立方體邊界棱線(xiàn)在柵格化后會(huì)被磨圓而且體素大小需要反復(fù)調(diào)。TIN完全基于原始點(diǎn)坐標(biāo)結(jié)構(gòu)面邊界保持精度更高。腳本輸出的data_triangle.txt保存的正是這個(gè)TIN的頂點(diǎn)索引每一行三個(gè)數(shù)字對(duì)應(yīng)一個(gè)三角形的三個(gè)頂點(diǎn)編號(hào)編號(hào)指向重采樣后的點(diǎn)云行號(hào)。拿到這個(gè)文件后先檢查三角形是否覆蓋了所有有效點(diǎn)有些情況下邊緣區(qū)域會(huì)出現(xiàn)懸空三角形需要在后續(xù)處理中排除。2.2 data_normals.txt 與法向量估算PCA方法與鄰域大小法向量是結(jié)構(gòu)面提取的核心依據(jù)。對(duì)每個(gè)點(diǎn)收集它在TIN上的鄰接頂點(diǎn)構(gòu)成局部點(diǎn)集對(duì)點(diǎn)集做PCA主成分分析最小特征值對(duì)應(yīng)的特征向量就是該點(diǎn)的法向量。data_normals.txt里每一行是歸一化后的三維法向量分量nx, ny, nz歸一化意味著nx2ny2nz21。拿到這個(gè)文件后我一般先檢查法向量是否朝向同一個(gè)半空間比如全部翻轉(zhuǎn)成z分量非負(fù)避免后面聚類(lèi)時(shí)出現(xiàn)方向歧義。鄰域大小直接決定法向量質(zhì)量。鄰域太大法向量被跨越結(jié)構(gòu)面邊界的點(diǎn)污染鄰域太小噪聲占主導(dǎo)。實(shí)際處理時(shí)鄰域半徑取平均點(diǎn)間距的2到3倍或者直接用固定k值。點(diǎn)云密度均勻時(shí)k取20到30效果都不錯(cuò)密度波動(dòng)大就按半徑搜索而不是按k。腳本里這個(gè)參數(shù)通常叫radius或k_neighbors建議先跑一組對(duì)比實(shí)驗(yàn)再定不同巖性的表面粗糙度對(duì)最優(yōu)鄰域大小影響很明顯。2.3 密度場(chǎng)與距離場(chǎng)邊界處的隱性約束data_density.txt記錄每個(gè)點(diǎn)附近的點(diǎn)云密度單位通常是點(diǎn)/平方米。密度場(chǎng)有兩個(gè)作用一是識(shí)別噪聲區(qū)掃描時(shí)被灌木遮擋形成的碎點(diǎn)密度極低二是檢測(cè)結(jié)構(gòu)面邊界因?yàn)榻Y(jié)構(gòu)面交界處往往會(huì)因?yàn)檎趽醍a(chǎn)生密度突變。區(qū)域生長(zhǎng)時(shí)如果兩個(gè)鄰接點(diǎn)的密度比值超過(guò)1.5到2倍就應(yīng)該停止生長(zhǎng)避免跨過(guò)邊界。data_distance.txt記錄每個(gè)點(diǎn)到某個(gè)局部擬合平面的垂直距離。這個(gè)距離場(chǎng)在處理緩傾結(jié)構(gòu)和弧形巖面時(shí)非常關(guān)鍵——單一法向量約束下弧形面上的點(diǎn)會(huì)被錯(cuò)誤歸為同一個(gè)面但距離場(chǎng)會(huì)在彎曲處產(chǎn)生明顯抬升把弧形拆成若干近似平面。下表匯總幾個(gè)文件在管線(xiàn)中的角色數(shù)據(jù)文件內(nèi)容在提取管線(xiàn)中的作用data_resample.txt重采樣后的點(diǎn)坐標(biāo)統(tǒng)一點(diǎn)間距降低計(jì)算量保證鄰域搜索一致性data_triangle.txtTIN三角形頂點(diǎn)索引提供頂點(diǎn)級(jí)拓?fù)溧徑雨P(guān)系供區(qū)域生長(zhǎng)遍歷data_normals.txt每個(gè)點(diǎn)的歸一化法向量聚類(lèi)與生長(zhǎng)的核心判定依據(jù)data_density.txt鄰域點(diǎn)密度值識(shí)別低密度噪聲與邊界突變data_distance.txt點(diǎn)到擬合平面的垂直距離拆分弧形結(jié)構(gòu)約束面片合并整套管線(xiàn)順序是重采樣 → 構(gòu)建TIN → 估算法向量 → 計(jì)算密度場(chǎng)和距離場(chǎng) → 自動(dòng)提取結(jié)構(gòu)面。前四個(gè)文件都是中間產(chǎn)物最后一步Automatic extraction of discontinuities.py讀入這些文件輸出結(jié)構(gòu)面分組與產(chǎn)狀。實(shí)際跑數(shù)據(jù)時(shí)重采樣這一步最容易被跳過(guò)導(dǎo)致后面的鄰域參數(shù)在不同區(qū)域完全不可比。數(shù)據(jù)文件與code的配合方式直接決定了提取效果的上限。3. 自動(dòng)提取算法實(shí)現(xiàn)Python區(qū)域生長(zhǎng)與面片分割3.1 讀入TIN并構(gòu)建鄰接表自動(dòng)提取的第一步不是聚類(lèi)而是把TIN的拓?fù)潢P(guān)系轉(zhuǎn)換成可以快速遍歷的數(shù)據(jù)結(jié)構(gòu)。data_triangle.txt里有M個(gè)三角形每行三個(gè)頂點(diǎn)索引需要先把它轉(zhuǎn)成每個(gè)頂點(diǎn)對(duì)應(yīng)的鄰接頂點(diǎn)列表。用Python實(shí)現(xiàn)這一步很直接def build_adjacency(n_points, triangle_file): adj [[] for _ in range(n_points)] with open(triangle_file, r) as f: for line in f: parts line.split() if len(parts) 3: continue i, j, k (int(p) for p in parts[:3]) adj[i] [j, k] adj[j] [i, k] adj[k] [i, j] # 去重避免重復(fù)鄰接關(guān)系拖慢生長(zhǎng) return [list(set(nb)) for nb in adj]這段代碼把每個(gè)三角形展開(kāi)成三條邊再把邊兩端的頂點(diǎn)互加為鄰居。set去重是因?yàn)樵诿芗蔷W(wǎng)里一個(gè)頂點(diǎn)的鄰居可能超過(guò)20個(gè)重復(fù)索引會(huì)在后續(xù)遍歷中產(chǎn)生大量無(wú)效訪(fǎng)問(wèn)。鄰接表構(gòu)建完成后區(qū)域生長(zhǎng)每次訪(fǎng)問(wèn)頂點(diǎn)時(shí)只需查這個(gè)列表時(shí)間復(fù)雜度從O(M×k)降到O(N)。3.2 區(qū)域生長(zhǎng)從種子點(diǎn)開(kāi)始的同向面片擴(kuò)張有了鄰接表和法向量剩下的核心問(wèn)題是哪些點(diǎn)屬于同一個(gè)結(jié)構(gòu)面最常用的方法是區(qū)域生長(zhǎng)。先選一個(gè)種子點(diǎn)從種子點(diǎn)出發(fā)把法向量夾角小于閾值的鄰接點(diǎn)并入當(dāng)前面片再以這些新點(diǎn)繼續(xù)向外擴(kuò)展直到?jīng)]有滿(mǎn)足條件的鄰居為止。def region_growing(adj, normals, angle_threshold_deg25.0, min_points50): n len(normals) labels np.full(n, -1, dtypeint) cos_thr np.cos(np.deg2rad(angle_threshold_deg)) label 0 for seed in range(n): if labels[seed] ! -1: continue labels[seed] label stack [seed] while stack: p stack.pop() for q in adj[p]: if labels[q] ! -1: continue # 用點(diǎn)積判斷法向量夾角是否在閾值內(nèi) if np.dot(normals[p], normals[q]) cos_thr: labels[q] label stack.append(q) if np.sum(labels label) min_points: labels[labels label] -1 # 丟棄過(guò)小面片 label - 1 label 1 return labels這里用當(dāng)前點(diǎn)p的法向量與鄰居q的法向量做點(diǎn)積而不是用種子點(diǎn)的法向量做全局比較。這樣生長(zhǎng)路徑可以沿結(jié)構(gòu)面自然彎曲延展適合巖體表面不是絕對(duì)平面的情況。angle_threshold_deg的物理含義是相鄰微面的最大夾角偏差工程上20到30度是常見(jiàn)區(qū)間角度越小分割越碎角度越大越容易合并多個(gè)結(jié)構(gòu)面。min_points用于過(guò)濾孤立碎面這些通常是噪點(diǎn)或掃描碎片。3.3 用密度場(chǎng)和距離場(chǎng)修正邊界只靠法向量的區(qū)域生長(zhǎng)有兩個(gè)典型失敗模式一是相鄰結(jié)構(gòu)面產(chǎn)狀接近法向量夾角只有四五度生長(zhǎng)會(huì)順著接縫跨過(guò)去二是弧形巖面被整體歸成一個(gè)面產(chǎn)狀卻在空間上連續(xù)變化。這兩個(gè)問(wèn)題需要密度場(chǎng)和距離場(chǎng)兜底。生長(zhǎng)過(guò)程中增加兩個(gè)約束條件。第一如果兩個(gè)鄰接點(diǎn)的密度比值超過(guò)density_ratio經(jīng)驗(yàn)值1.5到2.0即使法向量夾角滿(mǎn)足閾值也停止生長(zhǎng)。第二記錄每個(gè)點(diǎn)相對(duì)當(dāng)前結(jié)構(gòu)面擬合平面的距離距離超過(guò)max_distance通常取平均點(diǎn)間距的1到2倍的點(diǎn)不能并入。實(shí)現(xiàn)上密度檢查在入棧前判斷距離檢查在面片完成生長(zhǎng)后做一次離群點(diǎn)剔除兩輪串行可以讓邊界更干凈。3.4 運(yùn)行腳本與參數(shù)一覽資源里的主腳本可以直接從命令行跑數(shù)據(jù)文件作為參數(shù)傳入python Automatic_extraction_of_discontinuities.py \ --points data_resample.txt \ --triangles data_triangle.txt \ --normals data_normals.txt \ --density data_density.txt \ --distance data_distance.txt \ --angle-threshold 25 \ --min-points 100 \ --density-ratio 1.8 \ --output discontinuity_set.txt參數(shù)含義如下表參數(shù)取值建議影響--angle-threshold20–30角度越大合并越多越小越碎--min-points50–200過(guò)濾小面片噪聲--density-ratio1.5–2.0控制邊界是否跨越密度突變區(qū)--max-distance平均點(diǎn)間距×1~2控制點(diǎn)到擬合面的最大垂直距離跑完之后會(huì)輸出每個(gè)結(jié)構(gòu)面的編號(hào)、包含點(diǎn)數(shù)、擬合法向量與產(chǎn)狀。第一次跑建議用默認(rèn)參數(shù)先出結(jié)果再根據(jù)輸出結(jié)構(gòu)面數(shù)量反推閾值方向而不是一開(kāi)始就追求一次到位。4. 產(chǎn)狀計(jì)算與精度驗(yàn)證SVD擬合、傾向傾角與人工對(duì)比4.1 SVD擬合平面從點(diǎn)集到位姿一個(gè)結(jié)構(gòu)面的產(chǎn)狀本質(zhì)上是擬合平面法向量的問(wèn)題。把結(jié)構(gòu)面內(nèi)的所有點(diǎn)收集起來(lái)做中心化處理后進(jìn)行SVD分解最小奇異值對(duì)應(yīng)的右奇異向量就是平面法向量。def fit_plane_normal(points): centroid np.mean(points, axis0) centered points - centroid _, _, vt np.linalg.svd(centered, full_matricesFalse) normal vt[-1] # 最小奇異值對(duì)應(yīng)的右奇異向量 if normal[2] 0: # 統(tǒng)一向上 normal -normal return normal, centroidSVD比直接求協(xié)方差矩陣特征分解更穩(wěn)定尤其當(dāng)點(diǎn)集接近退化比如點(diǎn)分布近似一條線(xiàn)時(shí)協(xié)方差矩陣可能接近奇異SVD的數(shù)值行為更穩(wěn)健。如果提取階段把邊界噪聲點(diǎn)也包含了進(jìn)來(lái)SVD擬合時(shí)這些離群點(diǎn)會(huì)拉偏法向量。我一般會(huì)在擬合前用RANSAC迭代剔除離群點(diǎn)或者用距離場(chǎng)文件判斷哪些點(diǎn)是結(jié)構(gòu)面內(nèi)部的可靠點(diǎn)。4.2 從法向量到產(chǎn)狀傾向與傾角的換算地質(zhì)上產(chǎn)狀用傾向Dip Direction和傾角Dip Angle描述也有工程習(xí)慣用走向Strike加傾角。從法向量求產(chǎn)狀沒(méi)有歧義只需要一個(gè)坐標(biāo)變換def normal_to_orientation(normal): n normal / np.linalg.norm(normal) dip_angle np.degrees(np.arccos(np.clip(n[2], -1.0, 1.0))) dip_dir np.degrees(np.arctan2(n[0], n[1])) if dip_dir 0: dip_dir 360.0 strike (dip_dir 90.0) % 360.0 return dip_angle, dip_dir, strike這里傾角是法向量與豎直向上方向的夾角傾向是法向量水平投影的方位角從北方向順時(shí)針計(jì)算?？匆粋€(gè)實(shí)際數(shù)字如果法向量為(0.34, 0.47, 0.81)傾角約36度傾向約36度走向約126度。拿到這些參數(shù)后可以對(duì)照野外羅盤(pán)記錄驗(yàn)證。需要注意不同軟件對(duì)走向的定義存在差異DIPS和CloudCompare里展示的走向可能基于右手規(guī)則建議導(dǎo)出時(shí)確認(rèn)參考系。4.3 與人工測(cè)量數(shù)據(jù)的對(duì)比驗(yàn)證驗(yàn)證自動(dòng)提取結(jié)果常見(jiàn)做法是在CloudCompare里手動(dòng)框選同一組結(jié)構(gòu)面擬合平面后讀取產(chǎn)狀再與自動(dòng)結(jié)果對(duì)比。統(tǒng)計(jì)指標(biāo)一般看三點(diǎn)對(duì)比項(xiàng)可接受誤差傾角差值小于5度傾向差值小于10度結(jié)構(gòu)面數(shù)量匹配率大于80%數(shù)量匹配率低時(shí)先判斷是不是過(guò)分割導(dǎo)致同一個(gè)面被切成幾片。如果是把a(bǔ)ngle_threshold調(diào)大、min_points調(diào)大如果欠分割導(dǎo)致多個(gè)面并成一個(gè)則調(diào)小角度閾值。這個(gè)驗(yàn)證過(guò)程同時(shí)也是標(biāo)定參數(shù)的過(guò)程每個(gè)場(chǎng)地因?yàn)閹r石完整性和風(fēng)化程度不同最優(yōu)閾值會(huì)略有差異同一套參數(shù)換一個(gè)工地后重新標(biāo)定是常態(tài)。5. 參數(shù)調(diào)優(yōu)與巖體結(jié)構(gòu)面提取的邊界問(wèn)題5.1 關(guān)鍵參數(shù)的調(diào)整順序?qū)嶋H工程里參數(shù)調(diào)整有先后順序。先固定重采樣間距和鄰域半徑保證法向量穩(wěn)定再調(diào)angle_threshold觀(guān)察結(jié)構(gòu)面數(shù)量變化最后用min_points和density_ratio清理小碎面。如果結(jié)果偏碎優(yōu)先增大min_points而不是增大角度閾值因?yàn)榻嵌乳撝底兇髸?huì)讓真正獨(dú)立的相鄰結(jié)構(gòu)面合并。處理多站拼接點(diǎn)云時(shí)先按站點(diǎn)分別提取再合并結(jié)果比直接處理全部點(diǎn)云更容易控制誤差。5.2 植被遮擋與掃描盲區(qū)的處理植被是點(diǎn)云結(jié)構(gòu)面提取最大的干擾源。低矮灌木和樹(shù)枝會(huì)在巖面上方形成一層碎點(diǎn)密度比巖面低但法向量雜亂。預(yù)處理階段用密度閾值過(guò)濾低密度點(diǎn)能去掉大部分植被點(diǎn)。掃描盲區(qū)是另一個(gè)問(wèn)題結(jié)構(gòu)面被遮擋后只剩部分點(diǎn)擬合出的產(chǎn)狀往往偏向可見(jiàn)部分這時(shí)不能憑點(diǎn)數(shù)判斷可靠性要看每個(gè)結(jié)構(gòu)面的點(diǎn)云覆蓋范圍是否均勻必要時(shí)補(bǔ)測(cè)。5.3 把結(jié)果接入DIPS與CloudCompare自動(dòng)提取的輸出通常是帶產(chǎn)狀和點(diǎn)集的面片列表。接入DIPS分析時(shí)直接導(dǎo)入傾向和傾角兩列即可DIPS會(huì)生成極點(diǎn)圖和赤平投影。若要在CloudCompare里可視化面片可以按結(jié)構(gòu)面編號(hào)把點(diǎn)云分組導(dǎo)出為多個(gè)文本文件再以標(biāo)量場(chǎng)形式著色分段。最后給一個(gè)實(shí)用檢查看結(jié)構(gòu)面邊界與TIN三角形的關(guān)系。如果邊界位置出現(xiàn)大量細(xì)長(zhǎng)三角形說(shuō)明該處點(diǎn)云密度不夠或存在重疊掃描產(chǎn)狀結(jié)果可信度要打折。這個(gè)檢查只需統(tǒng)計(jì)輸出面片內(nèi)的平均三角形邊長(zhǎng)與整體點(diǎn)云平均邊長(zhǎng)的比值比值超過(guò)2時(shí)優(yōu)先補(bǔ)掃描或者調(diào)低該區(qū)域的權(quán)重。本文還有配套的精品資源點(diǎn)擊獲取