龐慶剛 車德福 賈慶仁2

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

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

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

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

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

1.1 地層建模

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

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

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

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

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

1.2 斷層建模

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

1.2.1 步驟1

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

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

1.2.2 步驟2

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

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

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

1.2.3 步驟3

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

2 地層建模結(jié)果

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

2.1 研究區(qū)概況

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

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

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

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

在地層可視化的基礎(chǔ)上加入斷層數(shù)據(jù)，建立了研究區(qū)內(nèi)含斷層的多個(gè)地層模型（編號分別為5th、7th、8th、9th和12th）。圖7為研究區(qū)部分區(qū)域的地層模型剖面。從剖面及等高線能看出多個(gè)地層被逆斷層切割，造成了地層表面的拓?fù)洳贿B續(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ìn)行對比分析。試驗(yàn)數(shù)據(jù)集包含從研究區(qū)獲得的2 478個(gè)點(diǎn)位坐標(biāo)及其單位法向量，其中勘探鉆孔數(shù)據(jù)、開采素描圖數(shù)據(jù)密度差異較大，具體數(shù)據(jù)及密度如表1所示。

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

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

4 結(jié)論

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

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

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

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