基于多尺度CSRBFs的地層三維地質(zhì)建模

龐慶剛 車德福 賈慶仁2

（1.開灤集團有限責(zé)任公司錢家營礦業(yè)分公司，河北唐山063301；2.東北大學(xué)資源與土木工程學(xué)院，遼寧沈陽110819；3.東北大學(xué)深部金屬礦山安全開采教育部重點實驗室，遼寧沈陽110819）

礦山開采過程中，地層三維模型的建立對于儲量管理、災(zāi)害評估和礦山開采規(guī)劃有重要意義［1-3］。在礦山生命周期中，不同階段、來源的數(shù)據(jù)具有明顯的尺度差異［4］，如在勘查找礦時期，往往通過大范圍布設(shè)鉆孔獲取地層分布信息；進入到設(shè)計采掘流程后，通過井下工作來精確獲取局部范圍內(nèi)地層的位置和厚度等數(shù)據(jù)。因此，是否能夠融合多尺度數(shù)據(jù)、快速建立地層模型，成為評價建模方法優(yōu)劣的重要因素［5-6］。目前已有多種方法用于地層建模，如基于廣義三棱柱的礦體建模法［7］、基于層面插值的建模方法［8-10］、動態(tài)更新的地層建模方法［11-12］等。上述方法獲得地層信息需要人工干預(yù)，容易引入人為誤差。文獻［13］提出了經(jīng)驗貝葉斯克里金方法（Empirical Bayes Kriging，EBK），可通過構(gòu)造子集和模擬的過程來自動選擇參數(shù)，減少人工交互，但是該方法較為復(fù)雜、計算量大，在地層三維建模應(yīng)用中實現(xiàn)困難。

近年來，隱式建模在三維建模領(lǐng)域受到了關(guān)注，該方法可自動從樣本分布中重建三維標量場f(x,y,z)，并提取其零等值面（f=0）作為模型表面，通過多邊形網(wǎng)格化方法可轉(zhuǎn)化為光滑的網(wǎng)格表達。諸多隱式方法中，緊支撐徑向基隱式函數(shù)（Compactly Supported Radial Basis Function，CSRBF）作為一種精確插值器，可以保證插值結(jié)果嚴格遵守地質(zhì)數(shù)據(jù)約束［14］。Ohtake等［15］基于八叉樹數(shù)據(jù)結(jié)構(gòu)建立了多尺度的緊支撐徑向基函數(shù)方法，在快速求解CSRBF的同時，能夠修復(fù)數(shù)據(jù)缺失區(qū)域，這對于數(shù)據(jù)量大、尺度差異明顯的地層數(shù)據(jù)有較好的實際應(yīng)用價值。

近年來，多尺度CSRBFs插值方法已經(jīng)被用于實現(xiàn)DEM、溫度場等的插值［16-17］，但在基于多源、多尺度數(shù)據(jù)的地層建模方面鮮有成果問世。本研究將多尺度CSRBFs用于地層表面三維建模，從大規(guī)模非均勻數(shù)據(jù)集中重構(gòu)標量場隱式函數(shù)，提取地層表面模型并進行可視化，并給出了地層斷層重構(gòu)方法。以錢家營礦部分區(qū)域數(shù)據(jù)作為試驗數(shù)據(jù)，通過與克里金方法對比，驗證了多尺度CSRBFs方法從多源、多尺度數(shù)據(jù)集中重建地層模型的優(yōu)勢。

1 多尺度CSRBFs地層三維建模方法

1.1 地層建模

應(yīng)用CSRBFs方法進行隱式建模所采用的數(shù)據(jù)約束主要為位置坐標與該位置處曲面的法向量［18］。使用這些約束后，該方法最終獲得樣品分布標量場的隱式函數(shù)，再通過多邊形可視化方法提取其零等值面，并轉(zhuǎn)換為三角形網(wǎng)格后在計算機中進行繪制。本研究提出的地層多尺度數(shù)據(jù)建模流程主要由4個步驟組成（圖1）。

（1）樣品點獲取與單位法向量計算。首先從目標地層相關(guān)的各類數(shù)據(jù)中提取樣品點，包括鉆孔、地質(zhì)剖面圖和井下素描圖等，得到采樣點數(shù)據(jù)集；然后通過計算獲得其單位法向信息。

（2）通過八叉樹空間索引構(gòu)建層次點集。通過遞歸劃分建?？臻g，建立八叉樹模型。對于八叉樹的第k層有數(shù)據(jù)的節(jié)點，計算質(zhì)心處的坐標和法向量，作為該層樣品點集。

（3）通過多尺度CSRBFs方法在樣品點集上由粗到細逐層插值，計算地層表面隱式函數(shù)。并且在插值各層樣品點集過程中，將誤差作為插值目標向下一層傳遞，逐層構(gòu)建隱式函數(shù)。

（4）地層隱式曲面可視化。提取隱函數(shù)零等值面，通過 Bloomenthal方法［18］將表面轉(zhuǎn)換為 Delaunay TIN，作為目標地層表面模型，并在計算機中渲染三角形。

1.2 斷層建模

斷層造成地層不連續(xù)，容易引起礦井生產(chǎn)事故，對于礦山生產(chǎn)影響很大［19］。本研究在地層曲面格網(wǎng)模型基礎(chǔ)上，通過在CD-TIN上加入斷層與地層底板交界線（斷層交線）分割地層，實現(xiàn)斷層三維建模，主要步驟如下。

1.2.1 步驟1

在地層底板TIN上根據(jù)斷層數(shù)據(jù)插值斷層中心線，并添加斷層交線。由中心線上的斷層點計算斷層交線上的點位平面坐標(x',y')：

式中，x,y,z為已知斷層點的坐標值；h為所求點的高程（點位相對向上移動時，h=z+l/2，點位相對向下移動時，h=z-l/2，l為斷層點的落差值）；α1為斷層傾角；α2為斷層傾向。

1.2.2 步驟2

根據(jù)斷層傳播規(guī)律計算斷層影響域，標記出影響域內(nèi)的局部三角網(wǎng)。在連接斷層交線確定斷層影響域的過程中，有以下兩種情形。

（1）情形一。如圖2（a）所示，P0為斷層中心線與地層邊界的交點，P1和P2點是由上式計算出的斷層交線上的點，P3和P4點為地層邊界邊的兩個端點。刪除邊界邊P3P4，連接邊P1P3和邊P2P4，如圖 2（b）所示。另一不封閉端進行同樣處理，得到斷層影響域和上下邊界，如圖2（c）所示。

（2）情形二。如圖3（a）所示，P0點為尖滅點，P1、P2、P3點分別為P0點所在三角形的3個頂點。刪除ΔP1P2P3，連接P0與P1點、P0與P2點、P0與P3點，重構(gòu)3個新的三角形，如圖3（b）所示。刪除斷層影響域內(nèi)的三角形，確定影響域上下邊界，如圖3（c）所示。

1.2.3 步驟3

移動影響域內(nèi)三角形角點的位置，通過加入拓撲不連續(xù)性的方式，構(gòu)建斷層三維實體模型。正斷層對上影響域邊界與上斷層邊界線、下邊界與下斷層邊界線（圖4（a）加粗曲線），逆斷層對上邊界與下斷層邊界線、下邊界與上斷層邊界線（圖4（b）加粗曲線）包圍的區(qū)域分別進行局部三角剖分。如果當前斷層的斷層中心線與已建模斷層的斷層邊界線相交，則認為兩個斷層相交。在相交斷層有逆斷層時，必須首先判斷第一個斷層的正逆情況，因為逆斷層區(qū)域水平投影有重疊，然后對相交斷層進行逐一建模（圖4（c）、圖4（d）和圖4（e））。

2 地層建模結(jié)果

基于本研究方法開發(fā)了原型系統(tǒng)［20］，系統(tǒng)包括數(shù)據(jù)錄入管理、數(shù)據(jù)插值、數(shù)據(jù)可視化、斷層建模等功能模塊。運行該系統(tǒng)，由研究區(qū)地質(zhì)資料建立三維地層模型。

2.1 研究區(qū)概況

錢家營礦井田面積約88 km2，井田范圍內(nèi)有8個礦層，共有勘探鉆孔259個。本研究選擇某礦層10采區(qū)作為研究區(qū)（3 km×3 km）。該區(qū)域地層相關(guān)數(shù)據(jù)主要有勘探鉆孔數(shù)據(jù)、井下測量導(dǎo)線數(shù)據(jù)和開采素描圖數(shù)據(jù)。從研究區(qū)各類數(shù)據(jù)中提取采樣點，如從鉆孔測斜數(shù)據(jù)中獲取地層采樣點，從地質(zhì)素描圖中獲取礦體頂板和底板位置，通過插值獲取采樣點等，共獲得2 478個不同尺度的某礦層相關(guān)數(shù)據(jù)，數(shù)據(jù)尺度差異如表1所示。采樣點的法向量通過PCL庫構(gòu)建K-D樹和計算K鄰域得到（K=10）。

2.2 地層模型構(gòu)建

通過對建模區(qū)域數(shù)據(jù)建立八叉樹空間索引，將數(shù)據(jù)劃分為7層（圖5）。每一層數(shù)據(jù)在其所處的數(shù)據(jù)層子節(jié)點內(nèi)計算其質(zhì)心坐標與法向量，然后逐層進行插值。最終，根據(jù)點集層次逐層計算的隱式函數(shù)即為樣品點所處的標量空間。

使用Bloomenthal方法進行地層隱式函數(shù)可視化的結(jié)果如圖6所示。由圖6可知：由于2074E工作面數(shù)據(jù)密度最大（包含回采時觀測的素描數(shù)據(jù)），其表面起伏也最為明顯。

在地層可視化的基礎(chǔ)上加入斷層數(shù)據(jù)，建立了研究區(qū)內(nèi)含斷層的多個地層模型（編號分別為5th、7th、8th、9th和12th）。圖7為研究區(qū)部分區(qū)域的地層模型剖面。從剖面及等高線能看出多個地層被逆斷層切割，造成了地層表面的拓撲不連續(xù)。

3 插值方法比較與分析

為了探討多尺度CSRBFs方法在多尺度數(shù)據(jù)建模上的優(yōu)勢，在筆記本電腦（2.5 GHz Intel? Core（TM）i5-4200M，4GB RAM，Intel HD顯卡）上，基于相同的數(shù)據(jù)集，使用反距離加權(quán)（Inverse distance weighted，IDW）、普通克里金（Ordinary Kriging，OK）方法與多尺度CSRBFs方法（Multi-scale CSRBF，MsCS?RBF）構(gòu)建了地層表面模型，并對3種方法在建模插值精度和時間效率上的差異進行對比分析。試驗數(shù)據(jù)集包含從研究區(qū)獲得的2 478個點位坐標及其單位法向量，其中勘探鉆孔數(shù)據(jù)、開采素描圖數(shù)據(jù)密度差異較大，具體數(shù)據(jù)及密度如表1所示。

（1）試驗 1。采用十折交叉驗證［21］對比分析IDW、OK、MsCSRBF 3種方法的建模精度。將所有采樣點分為10個子集，將其中9個作為建模數(shù)據(jù)集，另外1個作為驗證集。插值處理建模數(shù)據(jù)集，統(tǒng)計在驗證數(shù)據(jù)集上的插值誤差，并將其作為插值精度衡量指標。其中，誤差定義為平面坐標（x，y）對應(yīng)的插值點和原始采樣點之間的高程偏差絕對值，反映了插值的準確程度，誤差標準差則用于衡量插值結(jié)果的穩(wěn)定性。3種方法評價結(jié)果如表2所示，MsCSRBF方法的插值誤差、誤差標準差都小于另外兩種方法。是因為MsCSRBF法由于額外加入了法向數(shù)據(jù)約束，且使用八叉樹分層方式，使每一層的數(shù)據(jù)都盡可能均勻，因此插值結(jié)果更加平滑。

（2）試驗2。在不同規(guī)模（分別為10、50、250、1 250、2 478個）點集上運行3種插值方法的耗時對比如圖8所示。OK方法的變異函數(shù)采用高斯模型進行擬合，當采樣點很多時，普遍采用鄰域法限制參與計算的點數(shù)目。本研究采用象限法確定鄰域，以目標位置為原點，將空間在XOY平面劃分為4個象限，從每個象限取距離最近的2個點，以此限制每次參與插值的點數(shù)不大于8個。由圖8可知：MsCSRBF法運行耗時明顯少于OK、IDW法。

4 結(jié)論

基于多尺度CSRBFs插值方法，提出了一種地層三維建模方法。該方法通過使用八叉樹、構(gòu)造分層點集并逐層插值的方式，能夠從非均勻數(shù)據(jù)中快速插值獲得地層三維表面TIN模型。在此基礎(chǔ)上，提出了作用于地層表面模型上的、基于TIN網(wǎng)重構(gòu)的斷層建模方法。與以往的研究相比，使用該方法建立地層三維模型的優(yōu)勢有：

（1）通過在不同數(shù)據(jù)集上的交叉驗證，多尺度CSRBFs方法的插值精度優(yōu)于IDW、OK方法。

（2）隨著數(shù)據(jù)集增大，多尺度CSRBFs方法對非均勻數(shù)據(jù)集的插值運算效率優(yōu)于IDW、OK方法。

（3）基于提出的斷層建模方法，可在地層TIN網(wǎng)上構(gòu)建正、逆斷層模型和復(fù)雜相交斷層模型。