尖滅分叉型煤層三維地質建模研究:基于地質統計學和符號矩陣的方法

王艷美, 劉星

(安徽理工大學地球與環境學院, 淮南 232001)

三維地質建模(3D geology modeling),是運用計算機技術,在虛擬三維環境下將空間信息管理、地質解譯、空間分析和預測、地學統計、實體內容分析及圖形可視化等工具結合起來,用于地質研究的一門技術[1]。為實現煤炭精準開采智慧化,礦井透明地質條件是重要基礎[2]。煤礦智能化開采是中國煤炭行業升級轉型的必由之路[3],也是煤炭行業未來發展的大趨勢,三維地質建模及其可視化技術促進了地質資料信息的管理和共享,在煤炭行業的應用非常廣泛且具有極其重要的意義[4]。針對三維地質建模相關研究的各種空間建模技術和方法也不斷涌現,根據建模核心過程是否由數學模型控制,可以將其分為顯式建模和隱式建模兩類[5]。顯示建模方法以人機交互為基礎,通?；趯用婺Ｐ蚚6]來構建,該方法針對區域范圍小、構造簡單的近平行沉積地層應用效果良好,而針對區域范圍大、地層交錯復雜、在不同剖面上產狀變化大的地質要素進行建模時,需要大量的時間和高要求的專業素養,實際應用程度較低,且所建模型通常與實際地質構造出入較大,模型棱角尖銳,在多z值現象(多逆斷層)的地質體中人機都無法對相應地質界線進行連接,無法獲得實際的三維地質體,難以進行動態調整且不同的人圈定的地質體差異較大[7]。隱式建模主要核心為空間插值和模擬技術,空間插值的分類方式較多[8-9],常見的空間插值方法有克里金插值、距離反比插值和離散光滑插值等[10]。畢林等[11]提出了一種基于輪廓線和距離場的礦體隱式建模方法,并利用支持向量機和Poisson曲面進行了改進以提高這種方法的自動化程度。鄒艷紅等[12]提出一套基于隱函數曲面的地質體隱式建模方法,可以快速構建礦體三維模型。Zhang等[13]、唐丙寅等[14]提出一種新的基于不規則三角網與角點網格(triangulate irregular network and corner-point grid,TIN-CPG)的空間數據模型,并基于克里格插值方法實現了三維地層結構與屬性模型的一體化展示。郭甲騰等[15]利用全局徑向基函數構造的Coons曲面擬合地質特征,并利用布爾操作對模型進行裁剪優化,實現了三維模型的精細構建與系統集成。師素珍等[16]利用地震反演技術進行煤層的分叉合并,證明波阻抗反演對于煤層分叉合并及厚度變化具有較高的精度。空間插值方法針對地層平緩、地質構造簡單的地質體建模適用度較高,但對數據量的要求較高,若數據量稀疏容易出現地層穿插現象,且針對復雜地質體的三維地質模型構建過程復雜,難以達到理想狀態。因此,不能只單一使用空間插值方法,需要結合其他方法來解決上述問題。

鑒于此,提出一種基于地質統計學和符號矩陣的三維地質建模方法,以適用于大面積區域,不受煤層尖滅或分叉等復雜現象的限制,能夠精確定位煤層的位置,并且不限制區域內煤層的層數,使煤層之間的連通性更好。研究成果為煤炭資源精準透明開采設計提供依據。

1 方法

1.1 符號矩陣

符號矩陣是由4個符號學要素組成的顯示矩形圖式。它是一種不同于數值矩陣的特殊符號表達式[17],其結果可在軟件中建模。符號矩陣的計算在許多學科的理論研究與科技計算中都占有重要地位,其一般是二元對立要素。一些簡單的符號矩陣計算可手工完成,然而實際的科研工作中往往會產生大量冗繁的符號矩陣或表達式,對此手工計算幾乎是不可能的,需要一種不同于數值軟件的符號計算工具,對這種符號形式的矩陣進行處理并得到正確結果。主要利用MATLAB軟件進行計算,首先對研究區未知位置進行符號距離函數標記,然后再將符號距離函數以{-1, 1}的符號矩陣的形式來標記煤層和非煤層位置,從而有效標識煤層分布情況。

1.2 泛克里金插值方法簡介

克里金插值方法是法國G.Matheron教授以南非礦業工程師D.G.Krige(克里格)名字命名的一種對空間分布數據求最優、線性、無偏內插估計量的方法,是地質統計學的重要組成部分,也是地質統計學的核心[18]。由于沉積環境等情況的復雜性,區域化變量在實際應用中很難滿足二階平穩假設,因此選用泛克里金插值方法對待插點煤厚進行預測。在假設非平穩區域化變量Z(x)的變異函數存在,并已知Z(x)的漂移形式為多項式的條件下,討論各種線性無偏最優估計的問題,稱為泛克里金法。區域變量Z(x)在待插值點x0處的真實值Z(x0)是未知的,因此需要用x0影響范圍內的測量值的加權平均來估計待插值點的值。

(1)

式(1)中:Z*(x0)為待插值點x0處的煤厚估計值;Z(xi)為已知鉆孔的實際測量值;Wi為對應的權重。

求Z*(x0),需要先求權重,其權重不僅取決于待插點與已知點之間的距離h,還與預測區域內的地質環境和變異結構密切相關。變異函數與x無關,只與分隔兩點之間的距離h有關,定義為

(2)

式(2)中:h為滯后距;E為方差。

球狀模型、指數模型、高斯模型是常用的擬合函數模型[19]。

(1)球狀模型的表達式為

(3)

式(3)中:C0為塊金效應;C為總基臺值;a為變程。

(2)指數模型的表達式為

(4)

(3)高斯模型的表達式為

(5)

泛克里金矩陣表示為

A×W=B

(6)

式(6)中:

其中第n個鉆孔到未知待預測第p個鉆孔的距離,p為待測區域中未知待預測位置的個數;Xp和Yp為泛克里金估值的位置;a1和a2為待求的偏差參數;Xi、Yi為已知點的坐標,i=1,2,…,n。

求出權重之后,則可利用式(1)求出待插值點處的煤厚。

1.3 反距離加權插值法

反距離加權插值法(inverse distance weighting IDW)是常用的空間插值方法之一,由美國國家氣象局提出,是以待插點與已知點之間的距離為權重的插值方法,離待插點越近的已知點被賦予的權重越大,其權重值與距離成反比,其一般表達式為

(7)

(8)

(9)

p越大時較近的已知點將被賦予更高的權重,p越小時權重則更均勻地分配給各已知點。當p=1時稱為距離反比法,是一種常用且簡便的空間插值方法;當p=2時稱為距離平方反比法,選取p=2作為冪指數,來求待插值點的距離函數估值[20]。

2 建模流程

主要利用地質統計學方法和符號矩陣方法對于尖滅、分叉煤層進行三維地質建模,研究的關鍵是通過已知鉆孔數據,利用符號矩陣方法按照一定精度離散煤層和非煤層,通過泛克里金法和反距離加權法對煤層位置進行預測,建立三維地質模型。其流程如圖1所示。

圖1 模型流程圖Fig.1 Model flow chart

2.1 數據準備

通過區域中已知鉆孔的實際測量值,求出它的聚煤量和距離函數值,其計算方式如圖2所示,聚煤量是指在鉆孔處與煤層相交的總厚度,其計算公式為

a、b為鉆孔編號;在鉆孔a中,該鉆孔下有一層煤,T為鉆孔a的聚煤量,t1為鉆孔a與煤層相交的煤厚,此時T=t1; 在鉆孔b中,該鉆孔下有四層煤,T為鉆孔b的聚煤量,t1為鉆孔b與第一層煤層相交的煤厚,t2為鉆孔b與第二層煤層相交的煤厚,t3為鉆孔b與第三層煤層相交的煤厚,t4為鉆孔b與第四層煤層相交的煤厚,此時T=t1+ t2+ t3+t4圖2 數據預處理示意圖Fig.2 Data preprocessing diagram

(10)

式(10)中:T(a)為在鉆孔a處的聚煤量;n(a)為與鉆孔a相交的煤層層數;ti(a)為與鉆孔a中第i層相交的煤厚。

距離函數與到兩域界面的距離有關,兩域分別是煤層和非煤層。垂直方向上的柵格點的距離函數為其域到另一個域中最近柵格點的距離,本文定義煤層內距離函數為負值,非煤層內距離函數為正值,則假設兩域界面的距離函數為0。每個柵格點是二元對立的,則定義煤層內i(u)=1,非煤層內i(u)=-1。在柵格點u處,確定離它最近的柵格點u′,使得u-u′最小,即為已知點的鉆孔在柵格點u處的距離函數F(u),可表示為

(11)

在兩域界面處的距離函數為0,非煤層內的距離函數為正,其值在遠離界面時增大;煤層內的距離函數為負,其絕對值在遠離界面時增大。

2.2 煤厚計算

利用泛克里金法求待插點的煤厚。通過已知鉆孔的聚煤量,確定區域化理論變差模型,即可求出變異函數,進而代入式(6),求出權重值,該權重值可用于在每個待插點進行估計,最后利用式(1),求出區域內待插點的煤厚。

2.3 距離函數估值計算

對三維區域進行柵格化,為三維網格中的每個柵格點進行距離函數估值。指定一個垂直搜索窗口來限制使用哪些柵格點對當前柵格點進行估計,并且不限制水平搜索窗口,直接使用煤厚計算求出的權重值,將分配給已知鉆孔的泛克里金權重值應用于該鉆孔中位于垂直搜索窗口內的每個柵格點?？梢灾苯邮褂迷摍嘀刂?也可以利用反距離加權法對泛克里金權重值進行加權,使其在地質環境上更接近被估計柵格點的數據,無論如何權重值之和都為初始分配給已知鉆孔的泛克里金權重值。其距離估值的計算公式為

E*(u)=

(12)

式(12)中:E*(u)為待求柵格點的距離估值,F(ui,j)為第i個鉆孔垂直搜索窗口中第j個柵格點的距離函數;n為已知鉆孔個數;m為垂直搜索窗口中柵格點的個數;d(ui,j)為待求柵格點到第i個鉆孔垂直搜索窗口中第j個柵格點的距離。

2.4 根據距離函數估值匹配煤厚

通過計算三維網格中每個柵格點的距離函數可以確定煤層的位置。已知三維網格各列中每個柵格點的深度和距離函數,通過這兩個矢量即可確定煤層頂底板。為此,需要考慮距離函數閾值,距離函數估值位于距離函數閾值的相反側,煤層與非煤層的邊界位于閾值的位置?？紤]距離函數閾值為0,其示例如圖3(a)所示,煤頂在4～4.5 m之間,其距離函數估值為0.35和-0.14;煤底在8～8.5 m,其距離函數估值為-0.27和0.23,則煤厚估計值是4 m。為了使煤厚精確度更高,需對煤頂和煤底進行精細化,如圖3(b)所示,初始垂直分辨率為0.5 m,煤厚估計值是4 m,令垂直分辨率為0.05 m,煤頂在4.35～4.4 m,其距離函數估值為0.01和-0.04;煤底在8.25～8.3 m,其距離函數估值為-0.02和0.03,則煤厚估計值是3.9 m,與泛克里金法預測的煤厚更加相近,但距離函數閾值為0 時的煤厚估值與泛克里金預測的煤厚并不總是相近,為此考慮距離函數閾值為0.5 m,其示例如圖4所示,煤頂在3.8～3.85 m,其距離函數估值為0.55和0.49;煤底在8.75～8.8 m,其距離函數估值為0.49和0.54,則煤厚估計值是4.95 m。由此可知,距離函數閾值與煤厚估計值成正比,增加距離函數閾值煤厚估計值也增加,減小距離函數閾值煤厚估計值也相應減小。其在多煤層中也同樣適用。將距離函數閾值應用于距離函數估值這一列中可以識別不同的煤層和不同的煤層厚度,其示例如圖5所示,低距離函數閾值識別

圖3 距離函數閾值為0Fig.3 Distance function threshold is 0

圖4 距離函數閾值為0.5Fig.4 Distance function threshold is 0.5

圖5 距離函數閾值圖Fig.5 Distance function threshold diagram

出兩個煤層,高距離函數閾值識別出一個煤層,聚煤量也在增加。

2.5 構建三維地質模型

可以繪制三維地質模型的軟件有很多,如RMS、GOCAD、EVS、ArcGIS軟件等,選擇GOCAD作為建模軟件[21],首先在該軟件上對三維地質模型進行柵格化,生成S-Grid網格模型[22],其為M×N×K個子單元的規則拓撲結構模型,其每個單元體有8個柵格點,其柵格點的坐標可以單獨儲存,將其以Object的形式導出,對符號矩陣進行計算,其計算公式為

(13)

式(13)中:Classi,j,k為第(i,j,k)個子單元的符號矩陣值,i=1,2,…,M;j=1,2,…,N;k=1,2,…,K,其構建出K個M×N的符號矩陣,將符號矩陣中的值對properties進行重新賦值,重新導入即可看到煤層分布情況。

3 應用實例

3.1 研究區地質背景

潘二煤礦隸屬于安徽省淮南市潘集區,礦區范圍為116°45′00″E～116°56′15″E,32°45′00″N～32° 50′00″N,其主要含煤地層為石炭系、二疊系,石炭系煤層發育差,無工業價值,二疊系煤層主要為山西組、下、上石河子組[23]。由于研究區范圍廣,數據量龐大,選取潘集煤礦西北部進行三維地質建模,主要研究6、7、8號煤層群在地下的分布狀況,其主要煤層為7-1煤層,厚0～4.95 m,可采指數為0.9,其煤層由于環境變化、后期侵蝕等出現分叉、尖滅情況,煤層結構較為復雜,含暗煤、亮煤,煤層頂底板均以泥巖為主。研究區共有21個鉆孔勘探點,其二維分布如圖6所示。

圖6 研究區鉆孔分布圖Fig.6 Borehole distribution map of the study area

3.2 模型匹配

已知其研究區內鉆孔的聚煤量及距離函數,將該區域內的已知鉆孔數據集稱之為區域化變量,對其進行模擬,求得變異函數如圖7所示,其理論模型為高斯模型, 塊金效應C0=0.001,總基臺值C=0.954,變程a=2 251.666。利用泛克里金法求出待插點的煤厚,并結合反距離加權法求待插點處的各柵格點的距離函數估值,選取距離函數閾值為0,垂直窗口的大小為5時,其光滑度更高,與泛克里金法求出的煤厚值最為匹配。

圖7 變異函數圖Fig.7 Variation function diagram

3.3 三維地質建模構建

根據研究區邊界范圍建立SGrid模型,利用符號矩陣方法,對該研究區內鉆孔數據進行預處理,按照一定精度離散煤層和非煤層,進行區域空間建模,通過泛克里金法確定二維網格煤厚,反距離加權法確定三維網格中的距離函數估值,局部窗口對煤層進行精準定位,確定與煤層厚度相匹配的距離函數截距,對各六面體單元屬性進行預測及結果可視化,對其進行剖切顯示,得到研究區三維地質模型如圖8所示,由于該區域不含復雜地質結構,整個區域煤層按一定趨勢平緩變化,其煤層分布比較集中,且煤層連通性也比較好,眾多薄煤層也可以精準定位,并在適當位置分叉、尖滅,所建三維模型符合實際規律。如圖9(a)所示,鉆孔14-2與鉆孔12-4附近有尖滅現象,在圖9(b)中鉆孔L4-1和鉆孔13-1中煤層附近有分叉現象,其過渡自然,符合地質規律。

圖8 研究區三維地質模型柵欄圖(z方向拉伸5倍)Fig.8 Three-dimensional geological model fence map of the study area(stretch 5 times in z direction)

圖9 三維地質模型中煤層尖滅、分叉剖面與實際鉆孔媒段對比Fig.9 Comparison between the coal seam pinch-out and bifurcation sections in the three-dimensional geological model and the actual borehole media section is drawn

創建了煤層剖面示意圖,水平250 m,垂直20 m,包括兩個相距150 m的鉆孔,每個鉆孔與4 m厚的煤層相交。圖10說明改變垂直窗口個數的效果,在每種情況下,兩個鉆孔相交接縫之間的垂直偏移和煤層位置都是相同的。初始煤層呈斷裂狀態,隨著垂直窗口大小的增加,接縫的連接變得平滑。結果表明:垂直窗口較小時,煤層會出現不連續現象,實際應用時,選擇垂直窗口為7時滿足建模要求,圖11展示了改變煤層位置的效果,在每種情況下,垂直窗口大小保持不變,煤層連通性隨著兩鉆孔相交煤層的位置的距離增加而變化。初始煤層呈平滑狀態,接縫似乎一直在拉伸,在兩鉆孔相交煤層的位置為10 m時,煤層開始斷裂,直到最后以約6 m的間距“斷裂”。兩鉆孔相交煤層距離在10 m之內,煤層不會出現斷裂現象,在超過10 m時出現斷裂,與地質規律一致。

圖10 垂直窗口變化后對煤層連通性影響Fig.10 Influence of vertical window change on coal seam connectivity

圖11 煤層位置變化后對煤層連通性影響Fig.11 Influence of coal seam connectivity after the change of coal seam position

3.4 模型檢驗

在研究區內隨機抽取3個參與建模鉆孔,將原始鉆孔數據與模型預測數據對比,其結果如表1所示。由表1可知,二者具有良好的一致性,精度約為0.3 m,滿足開采要求[24]。

表1 原始鉆孔數據與模型預測數據精度對比Table 1 Comparison of accuracy between original borehole data and model prediction data

4 結論

(1)以鉆孔數據為基礎,提出一種基于地質統計學和符號矩陣的三維地質建模方法,可以快速建立分叉、尖滅連接項,精度可靠,誤差小于0.3 m,滿足開采要求。

(2)本文方法使用符號矩陣方法來構建煤層幾何模型。利用符號矩陣方法計算出的距離函數被局部截斷匹配正確的煤層厚度,接縫自然向外擴張并分叉,局部窗口對對距離函數模型進行精細化,以實現高精度的單個煤層的標高,且跟傳統方法相比,本文方法減少人工交互,自動化程度提高,避免層面建模需要人工經驗,先解譯再連接的難題。

(3)實踐表明,通過在潘集煤礦中使用該方法,成功地預測出7號煤在地層下的分布情況,符合地質規律,該方法中的算法運行速度和精確定位煤層和非煤層邊界的能力使其成為一個強大的建模工具,它可以模擬地下的煤層分布,使地質體三維可視化,為煤炭資源精準透明開采設計提供依據。