普通克里格法在草桃背鈾礦床資源量估算中的應用

馮巧云 馬 濤 鄧 宇

(江西省地質局有色地質大隊,江西 贛州 341000)

地質統計學是以區域化變量理論作為理論基礎,以變異函數作為主要工具,對既具有隨機性又具有結構性的變量進行研究的科學,其核心即克里格法,是一種無偏的最小估計方差的估值方法[1-5]。隨著應用領域不斷擴大,地質統計學理論更加成熟,各種克里格方法、技術不斷出現,如參數地質統計學中的普通克里格法、泛克里格法、析取克里格法,非參數地質統計學的指示克里格法、概率克里格法等,其中普通克里格法是其最基本、應用最為廣泛的估值方法[6-10]。

本研究借助MinePlan3D礦業軟件(原MineSight礦業軟件),以隱式建模方法對贛南草桃背鈾礦床進行三維地質建模,在充分研究該礦床地質特征及品位變化規律的基礎上,制訂地質統計學實施方案,以鈾邊界品位0.01%圈定礦化域,利用普通克里格法對鈾品位進行限制性插值,構建礦床資源模型,實現礦床邊際品位動態化選擇。根據變異函數擬合結果,揭示區域化變量鈾品位在垂深方向上比水平方向上礦化連續性更好,為該區后續深部隱伏找礦工作提供有益參考。

1 礦床地質特征

研究區位于武夷山成礦帶與南嶺成礦帶交匯處,大余-會昌Ⅱ級EW向構造與廣昌-周田Ⅱ級NE向構造的復合部位,會昌盆地東緣中基性火山巖體帶和富城巖體接合處(圖1)。區內出露地層主要為震旦-寒武系變余長石石英砂巖、千枚巖、板巖、碳質和硅質板巖,局部可見混合巖化。上白堊系紅色砂礫巖、砂巖、粉砂巖不整合在震旦-寒武系之上,中基性火山巖則產于紅層下部。富城巖體是區內主要巖體,橫跨贛閩二省,東西長約60 km,南北寬18～22 km,總面積為1 095 km2。北部侵入于震旦、寒武系,南部與混合巖呈侵入接觸,東部侵入于中侏羅系漳平統,西部被會昌上白堊系紅層所覆蓋,主要巖性為中粗粒似斑狀黑云母花崗巖、細粒似斑狀黑云母花崗巖,后期脈巖廣泛發育,常見花崗斑巖、石英斑巖、煌斑巖、石英脈等。構造以斷裂構造為主,以EW、NE向和環狀構造組成構造骨架,斷裂性質以壓扭性為主。

圖1 研究區區域地質特征Fig.1 Regional geological characteristics of the study area

根據礦體特征、構造特征、礦物共生組合、硫同位素測量、瀝青鈾礦全分析和絕對年齡測定結果,以及鈾礦脈之間的穿插關系,可將草桃背鈾礦床分為產于火山機構基底的細粒似斑狀黑云母花崗巖中巖漿期后熱液型(Ⅰ類型)、產于火山震碎帶的火山熱液型(Ⅱ類型)、產于火山震碎帶之紫色帶與綠色帶界面附近的淋積疊加型(Ⅲ類型),以及產于火山斜坡的火山熱液型(Ⅳ類型)4種鈾礦化類型。其中,Ⅱ類型鈾礦化及Ⅲ類型鈾礦化為本礦床主要礦化類型。礦體規模大小不一,礦體形態多為似層狀、透鏡狀,礦體產狀多變,受火山機構控制,主要礦體由40m×25m工程間距所控制,礦體產狀總體為EW向展布,但與火山口活動關系密切的Ⅱ類型和Ⅲ類型產狀受火山構造控制,Ⅰ類型產狀為近EW走向,傾向N,Ⅳ類型則主要受NE向構造控制,各礦化類型傾角變化較大,但傾角均較陡,一般為40°～88°,在走向和傾向上礦化都很不均勻,礦體品位變化較大。

2 三維地質建模與樣品統計分析

2.1 數據庫建立

本研究數據庫建立過程中收集了草桃背鈾礦床歷年來的勘探資料,包括鉆探、坑探、槽探等地質資料。收集整理了547個鉆孔工程的7 212個樣品數據,按照軟件對原始信息數據的要求,共建立了孔口坐標表、鉆孔測斜表、鉆孔化驗表與鉆孔巖性表4個數據文件。所有數據均經過細致全面的檢查校對,可作為礦床建模的基礎數據。

2.2 三維地質模型構建

本次三維地質建模構建了巖體模型及礦化域模型。巖體建模采用了隱示建模方法,其原理是基于隱函數對離散的三維點數據進行插值實現表面重建[3]。其中高斯徑向基函數(Radial Basis Function,RBF)是常用的隱函數,最早由HARDY引入[4],國際通用三維礦業軟件如Leapfrog、MinePlan3D、Micromine、Surpac等也采用該函數作為隱式建模的插值函數[11-15]。MinePlan3D軟件的隱式建模功能使用RBF函數生成表面較光滑自然的實體,提供了地表、風化、侵入、斷層、不整合、品位殼、巖脈及地層等預設選項以及相互之間的地質時代關系約束功能。草桃背礦床地質數據豐富,通過隱式建模,高效立體展示了火山機構的分布及各類巖體構造的接觸關系,彌補了以往只能用二維圖件展示地質體導致的認識不足,且建模效率與質量遠高于顯式建模,為周邊類似復雜地質條件的火山口(小富足、上富等)探索提供了有效工具及參考樣板。本研究通過隱式建模構建的安山巖實體模型如圖2所示。

礦化域模型是礦床成因及其在三維空間上變化特征的具體反映,通過礦化域模型不僅能在三維視角下展示礦化現象的幾何空間形態,同時也為后續品位模型插值提供了限制約束域。本研究采用剖面解譯法進行建模,在各個剖面對礦化邊界進行界定,按照幾何法的圈連和外推原則,以鈾邊界品位0.01%在剖面上圈定低品位礦化域,根據礦化域剖面在水平面的投影分布圈定各水平分段。

2.3 樣品統計特征分析

2.3.1 原始樣品分析統計

利用地質統計學方法進行礦床資源儲量評估時,有必要對樣品數據進行統計特征分析。只有在樣品數據的統計特征確定以后,方可有效選擇地質統計學插值方法[2]。根據建立的數據庫,對原始樣品鈾品位進行了統計分析,結果如圖3所示。通過圖3可以看出原始樣鈾品位服從偏正態分布,故本研究選用普通克里格法進行插值計算。

圖3 含特高品位的原始樣的鈾品位直方圖Fig.3 Uranium grade histogram of the original sample containing extremely high grade

2.3.2 特高品位處理

對于應用普通克里格法插值的原始數據需要對特高品位進行處理,以保證參數估值精度。目前,常用的特高品位判斷方法有置信區間判斷法、圖形判別法、公式判別法等[2],本研究采用圖形判別法確定特高品位。從原始樣鈾品位累計概率分布曲線圖(圖4)可以看出,在99.9%分位數上,其直線分布出現明顯的轉折拐點,根據統計學方法對特高品位的判斷與替換方法,該99.9%分位數上對應的品位值1.66%即為特高品位臨界值,故以1.66%品位值作為上限值代替原始樣鈾品位中的特高品位。

圖4 原始樣鈾品位累計概率分布曲線Fig.4 Cumulative probability distribution curve of original uranium grade

2.3.3 樣品組合與統計分析

地質統計學要求有效數據必須建立在固定長度的支撐上,因此,有必要對樣長不等的鉆孔、坑探品位數據按一定的長度進行組合計算并進行分析[1]。組合樣長度一般根據礦山開采方法、夾石剔除厚度、礦床類型及品位在垂深方向的變異程度等確定;通常要求一個采礦分段(臺階)高度內需有1個或若干個組合樣;組合樣的統計數字特征需與原始樣相近,不取過大的組合樣長度是為了降低樣品組合導致的品位平均化程度,并確保組合樣數量滿足地質統計學計算變異函數需有較多的樣本個體為基礎的要求。本研究建模選取1.25 m樣品組合長度進行統計分析。組合樣與原始樣的統計特征值比較見表1。由表1可知:無論是否處理了特高品位,二者的均值、標準差及變化系數等都相近。組合樣品位分布直方圖(圖5)分析表明:組合樣鈾品位服從正態分布,且與原始樣品位直方圖分布近似,因而取1.25 m作為組合樣長度較適宜。經過特高品位處理后,組合樣品位的變化系數明顯降低,說明數據的穩定性得到加強,為后續變異函數分析奠定了基礎。

表1 1.25 m組合樣鈾品位統計及其與原始樣的統計比較Table 1 Uranium grade statistics of 1.25 m composite sample and its statistical comparison with original sample

圖5 處理特高品位后的組合樣鈾品位直方圖Fig.5 Uranium grade histogram of combined sample after processing ultra-high grade

3 變異函數計算及結構分析

3.1 變異函數模型擬合與結構分析

礦區勘查控制工程基本網度間距為40 m(走向)×25 m(傾向),鉆孔垂深方向鈾品位變異大,因此在計算變異函數時,基本滯后距取目前最密工程間距的1/2即12.5 m,搜索角為15°,進行0°、15°、30°、…、165°共12個水平方向以及垂直方向(傾角分別設為0°、15°、30°、…、75°)的變異函數計算。

根據鈾的試驗變異函數等值線圖(圖6)分析,結合礦體產狀,選擇0°與90°這一組試驗變異函數進行擬合。以90°方向為例,得出變異函數及其擬合結果如圖7所示,各方向上的參數取值見表2。

表2 組合樣鈾品位球狀模型變異函數參數Table 2 Parameters of variogram of spherical model of uranium grade for composite samples

圖6 鈾的試驗變異函數等值線(水平方向)Fig.6 Isoline of test variogram of uranium(horizontal direction)

圖7 組合樣鈾品位試驗變異函數曲線(水平0°方位、傾角30°)Fig.7 Variable function curve of combined sample uranium grade test (horizontal 0° azimuth,dip 30°)

鈾的試驗變異函數結果擬合表明:鈾品位變異函數具有幾何異向性特點,表現為在水平0°/傾角30°、水平90°/傾角30°和垂深方向上具有相同的塊金常數值和基臺值,但變程不同。據上述鈾變異函數幾何異向性參數,礦業軟件可自動進行套合。

3.2 搜索橢球體參數確定

搜索橢球體用于克里格法插值對鄰域樣品搜索時,搜索半軸(半徑)一般不宜超過變異函數的變程[2]。根據變異函數擬合確定的變程等參數,估值的普通克里格法估計鄰域的半徑長度分別為長軸方向40 m,短軸方向25 m,垂直方向1.8 m;搜索橢球體全域扇區數8個,每個待估塊段參與克里格法估計的最少樣品數有2件,最多15件。

3.3 交叉驗證

試驗得到的變異函數模型有必要進一步進行理論檢驗。交叉驗證結果見表3。由表3可知:實際值與估計值的均值偏差近于0;且殘差分布圖服從正態分布,符合交叉驗證的2個理論判據,表明用該組變異函數參數對組合樣鈾品位進行估計是無偏的,滿足區域化變量的內蘊假設[2],即通過地質統計學方法進行品位插值建模具有合理性。

表3 交叉驗證結果Table 3 Cross-validation results

4 塊段模型及克里格估值

待估塊段尺寸確定需要綜合考慮礦床類型、礦體規模、夾石剔除厚度、探礦工程間距等因素。該礦主礦體厚度一般為2～6 m,傾向N,傾角30°左右,地下采礦作業不允許過度貧化,塊段尺寸在Z方向不宜超過3 m。由于本礦床勘查以40 m(走向)×25 m(傾向)工程間距求控制資源量,故本研究取該投影網度最小間距的1/2,即10 m×10 m作為待估塊段水平X、Y方向的尺寸。因此,模型塊段尺寸設置為10 m(東向)×10 m(北向)×2.5 m(垂直方向)。

設置好搜索橢球體參數,并在礦化域邊界限制的約束前提下對鈾品位進行普通克里格法估值,通過估值得到每個塊段的克里格品位及塊段的克里格方差值,模型塊段估值的三維效果如圖8所示。圖8顯示鈾礦化非常明顯地環繞火山口分布,表明鈾礦化嚴格受火山機構控制,且鈾礦體主要分布在火山口西部。

圖8 模型塊段估值的三維效果Fig.8 Three-dimensional effect of model segment estimation

5 資源量估算及驗證對比

品位插值后,可快速統計出不同邊際品位下的資源量。草桃背礦床模型以鈾克里格邊際品位0.047%估算礦床資源量,與原用地質塊段法(邊界品位0.03%,最低工業品位0.05%)估算結果進行了對比,發現兩者估算的礦石量相對誤差為0.34%,金屬氧化物量相對誤差為3%,平均品位相對誤差為3.2%,均在允許誤差范圍內。草桃背礦床模型資源量估算結果的鈾克里格邊際品位與礦石量的關系如圖9所示。由圖9可知:在礦石量a4與a5之間的靠近a4處,品位-噸位曲線出現明顯拐點,其對應的邊際品位低于0.05%,約0.047%。

圖9 草桃背礦床鈾礦塊邊界品位與礦石量關系Fig.9 Relationship between uranium block boundary grade and ore quantity in Caotaobei Deposit

對本次地質統計學法資源量估算結果進行全局驗證,首先分別計算統計出間隔2.5 m的每個水平(中段)的組合樣品平均值及該水平塊段模型克里格法估計平均值,分別繪制連線圖;然后在垂直方向上對比相同水平兩條連線的吻合程度,結果如圖10所示。

圖10 草桃背礦床鈾品位全局驗證曲線Fig.10 Global verification curves of uranium grade in Caotaobei Deposit

由圖10可知:組合樣鈾品位曲線與克里格法鈾品位曲線吻合度好,反映出克里格法估計結果可靠性好。

6 結 論

(1)本研究使用普通克里格法估算草桃背鈾礦床資源量,并與原用地質塊段法的估算結果進行了對比,兩種方法估算的礦石量相對誤差為0.34%,金屬氧化物量相對誤差為3%,平均品位相對誤差為3.2%,根據品位-噸位關系圖,提出了用鈾克里格品位0.047%作為邊際品位的單指標方案,為礦山今后適應市場經濟需要調整邊際品位奠定了基礎。

(2)采用隱式建模方法構建了該礦床巖體三維地質模型,彌補了以往僅有二維認識的不足,直觀地表征了火山機構內各巖石的分布形態,對紅盆底部的構造變異部位和紫色蝕變帶分布有了更清晰的了解,為進一步探索紅盆及開展火山巖型鈾礦深部找礦工作提供了地質依據。

(3)根據試驗變異函數擬合得到的鈾品位球狀模型理論變異函數,其水平方向主軸變程為63 m,次軸變程為56 m,鉆孔方向變程為75 m,即垂深方向變程比水平方向大,說明垂深方向上的礦化比水平方向上礦化連續性更好。