一種條帶狀銅多金屬礦儲量估算插值方法

李懷良 庹先國， 蔣 鑫

(1.核廢物與環(huán)境安全國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，綿陽 四川621010;2.地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，成都 四川610059)

西藏墨竹工卡縣甲瑪矽卡巖－角巖型銅多金屬礦床是我國目前同類礦床中礦石品質(zhì)最好、成礦元素最多、探明的資源量已經(jīng)達(dá)到超大型規(guī)模的礦床［1］。目前，一般采用傳統(tǒng)斷面法、SD 法和普通克里格法［2-3］估算該地區(qū)的礦床儲量。傳統(tǒng)斷面法雖然具有一些不足，但經(jīng)過幾十年的發(fā)展，在礦產(chǎn)資源儲量估算方面已日趨完善，只要合理選取地質(zhì)參數(shù)，不但能夠簡化計(jì)算，而且能夠達(dá)到較高的精度［4-5］。近年來，SD 法在我國發(fā)展也十分迅速，北京恩地科技發(fā)展有限責(zé)任公司開發(fā)的SD 法礦產(chǎn)資源儲量計(jì)算系統(tǒng)已通過國土資源部組織的評審鑒定，該系統(tǒng)可廣泛應(yīng)用于西藏甲瑪銅多金屬礦的礦床勘查至開采完成各個階段的儲量計(jì)算。普通克里格法的優(yōu)勢也較為明顯，不但能夠給出儲量估算的精度，而且是無偏且估計(jì)方差最優(yōu)的［6-9］。

西藏甲瑪銅多金屬礦區(qū)內(nèi)地層復(fù)雜，地層曲面變化較大，甚至可能存在地層曲面相交和錯斷的現(xiàn)象，因此鉆孔數(shù)據(jù)常呈條帶狀分布［10-14］，并且多種金屬間會出現(xiàn)協(xié)同區(qū)域化以及多種金屬儲量分布相互影響的現(xiàn)象，因此有必要消除多種金屬的相關(guān)性以及條帶狀效應(yīng)造成的鉆孔兩端數(shù)據(jù)權(quán)值偏大的影響［15-17］。為此，采用權(quán)值校正［18］的協(xié)同克里格法對協(xié)同克里格法插值得到的權(quán)值進(jìn)行校正，在估值方差基本不變的前提下使得協(xié)同克里格法估值的權(quán)值更為合理。

1 協(xié)同克里格法

采用協(xié)同克里格法［6-10，19-20］進(jìn)行儲量估算時(shí)，對由K 個變量構(gòu)成的協(xié)同區(qū)域化變量{Zk(x)}(k = 1，2，3，…，K)，有相鄰m 個位于位置ui(i = 1，2，3，…，m)處的數(shù)據(jù)點(diǎn)，其中心點(diǎn)的品位數(shù)據(jù)設(shè)為Z(uk，i)(i= 1，2，3，…，m;k = 1，2，3，…，K)，在待估點(diǎn)Z(u0)處的協(xié)同克里格的估計(jì)量為

式中，λk，i為第k 個區(qū)域化變量中的第i 個樣本點(diǎn)的權(quán)值，由最小化方差求出，即

設(shè)樣品數(shù)據(jù)點(diǎn)ui(i = 1，2，3，…，m)與待估點(diǎn)u0的距離為di，對于式(1)中的權(quán)值，如果di＜dj，則λk，i＞λk，j(i，j = 1，2，3，…，m)，且有

該算法主要實(shí)現(xiàn)流程見圖1。

2 鉆孔條帶狀效應(yīng)

西藏甲瑪銅多金屬礦最主要的數(shù)據(jù)有地面鉆孔數(shù)據(jù)、鉆孔化驗(yàn)數(shù)據(jù)、小體重?cái)?shù)據(jù)等。原始數(shù)據(jù)為截至2012 年甲瑪?shù)V區(qū)勘探工程數(shù)據(jù)，包括58 個鉆孔共15 187 件樣品的數(shù)據(jù)。若某個鉆孔遇到斷層，鉆孔數(shù)據(jù)便會出現(xiàn)邊界截?cái)喱F(xiàn)象，如圖2 所示。

如圖3 所示，在進(jìn)行協(xié)同克里格插值時(shí)首先需要確定一個搜索范圍，搜索范圍在三維空間內(nèi)通常為一個橢球體的邊界截?cái)唷D1、圖2 所示的2 類數(shù)據(jù)可稱為鉆孔的有限條帶狀數(shù)據(jù)。

圖1 協(xié)同克立格法實(shí)現(xiàn)流程Fig.1 Implementation process of Co-Kriging interpolation algorithm

圖2 數(shù)據(jù)被斷層邊界截?cái)郌ig.2 Drilling data truncated by boundary faults

圖3 數(shù)據(jù)被搜索鄰域截?cái)郌ig.3 Drilling data truncated by search field

以Cu 和Ag 的協(xié)同區(qū)域化為例，對其中的一列有限的條帶狀數(shù)據(jù)進(jìn)行協(xié)同克里格插值計(jì)算，其中變差函數(shù)以及2 種元素的交叉變差函數(shù)均采用球狀模型計(jì)算。其中Cu 塊金效應(yīng)為0 的協(xié)同克里格權(quán)值如圖4 所示。

圖4 協(xié)同克里格插值引起的條帶狀效應(yīng)Fig.4 Banned effect caused by Co-Kriging interpolation

圖4 顯示的權(quán)值體現(xiàn)了采用協(xié)同克里格法估值時(shí)所產(chǎn)生的條帶狀效應(yīng)，對于單個鉆孔數(shù)據(jù)的估值，與待估點(diǎn)距離最遠(yuǎn)的鉆孔數(shù)據(jù)兩端樣品點(diǎn)的權(quán)值偏大，反之則普遍較小。根據(jù)克里格估值法中的屏蔽效應(yīng)［13-16］，該圖中計(jì)算的權(quán)值是不符合實(shí)際情況的，特別當(dāng)鉆孔數(shù)據(jù)分布于斷層中或搜索范圍截出不平穩(wěn)的鉆孔數(shù)據(jù)時(shí)，便會出現(xiàn)端點(diǎn)數(shù)據(jù)處的權(quán)值明顯大于或小于中間數(shù)據(jù)權(quán)值的現(xiàn)象，因此有必要對其進(jìn)行校正。

3 權(quán)值校正

對權(quán)值系數(shù)進(jìn)行約束需滿足式(3)中的最小化問題約束條件。采用擬牛頓法［11］對權(quán)值進(jìn)行校正，流程見圖5。

圖5 擬牛頓法校正權(quán)值流程Fig.5 Flow of weight correction by Quasi-Newton method

4 處理結(jié)果

4.1 數(shù)據(jù)預(yù)處理

采用權(quán)值校正協(xié)同克里格法估計(jì)礦床儲量時(shí)，應(yīng)首先選取區(qū)域化變量。區(qū)域化變量必須具備隨機(jī)性和結(jié)構(gòu)性的雙重性質(zhì)，由于西藏墨竹工卡縣甲瑪夕卡巖－角巖型銅多金屬礦床為Cu、Mo、Au、Ag、Pb、Zn 等復(fù)合多金屬礦床，并有大量的鉆孔化驗(yàn)數(shù)據(jù)，因此可選取金屬元素的品位作為區(qū)域化變量，本研究采用Cu 品位進(jìn)行分析。在礦產(chǎn)儲量和品位計(jì)算之前，對鉆孔數(shù)據(jù)進(jìn)行預(yù)處理，包括特異值處理與樣品組合處理。根據(jù)Cu 品位數(shù)據(jù)，在正態(tài)概率紙上繪制累積頻率曲線圖見圖6。

圖6 Cu 品位累積頻率曲線Fig.6 Cumulative frequency curve of Cu grade

由圖6 可知，Cu 品位累積頻率曲線基本符合正態(tài)分布。在直線上部存在少量離散點(diǎn)，且偏離直線，該類離散點(diǎn)為真正的特異值。此次估算采用頻率曲線法將所有高于上截品位的數(shù)據(jù)全部用上截品位值代替。取1.5% 為上截品位值，礦區(qū)Cu 品位高于1.5%的共有782 件樣品。上截品位代替前、后Cu 特異值統(tǒng)計(jì)分析結(jié)果見表1。

表1 Cu 特異值處理結(jié)果Table 1 Results of singular value of Cu

采用權(quán)值校正協(xié)同克里格法估算資源量時(shí)，組合樣品長度設(shè)定為6 m。組合樣品中Cu 品位特征值的統(tǒng)計(jì)結(jié)果見表2。

表2 組合樣品中Cu 品位特征值Table 2 Statistical results of Cu grade characteristic values of composite samples

4.2 協(xié)同克里格法估算儲量

4.2.1 變差函數(shù)及結(jié)構(gòu)分析

對組合樣品進(jìn)行統(tǒng)計(jì)分析，結(jié)果見表3。由表3可知，組合樣品中Cu 品位與Ag 品位相關(guān)性最大，因此以Ag 品位作為輔助變量對Cu 的品位進(jìn)行估值。

權(quán)值校正協(xié)同克里格法插值的第一步是得到Cu、Ag品位的變差函數(shù)及兩者的交叉變差函數(shù)。首先將所有鉆孔點(diǎn)兩兩組合成點(diǎn)對，計(jì)算出其距離值;然后以球狀模型擬合成理論試驗(yàn)變差函數(shù)曲線。Cu、Ag 品位項(xiàng)的3 個變差函數(shù)軸的方向定位及其計(jì)算結(jié)果見表4。

表3 組合樣本統(tǒng)計(jì)分析Table 3 Composite samples analysis results

表4 Cu 和Ag 變差函數(shù)參數(shù)Table 4 Variation function Parameters of Cu and Ag

由表4 可知，2 種元素在各個方向上的變差函數(shù)參數(shù)相似，所有變差函數(shù)均為各向同性的。對Cu、Ag品位項(xiàng)的3 個軸方向上的變差函數(shù)進(jìn)行結(jié)構(gòu)套合，得到Cu、Ag 品位的變差函數(shù)模型

式中，h 為滯后距，m。

以上述2 種元素品位為基礎(chǔ)，用球狀模型擬合計(jì)算出交叉變差函數(shù)的基臺值為0.58 m，拱高為0.42 m，變程為55 m，其交叉變差函數(shù)曲線見圖7。

圖7 交叉變差函數(shù)曲線Fig.7 Cross-variogram function curve

4.2.2 權(quán)值校正

選取鉆孔編號為2012jm－8576 的條帶狀數(shù)據(jù)為測試數(shù)據(jù)，選取鉆孔中間位置為待估位置，求解協(xié)同克里格方程組得到的權(quán)重系數(shù)，按圖5 所示的擬牛頓法進(jìn)行權(quán)值校正，單個鉆孔校正后的權(quán)值見圖8。

圖8 權(quán)值校正結(jié)果Fig.8 Results of weight calibration

由圖8 可知，待估樣品點(diǎn)位于197#樣品點(diǎn)附近，距離待估點(diǎn)附近的權(quán)值偏大，反之則權(quán)值越小，符合實(shí)際情況。

4.2.3 算法模型驗(yàn)證

對甲瑪?shù)V區(qū)化學(xué)樣品中Cu 品位估值情況進(jìn)行交叉驗(yàn)證，結(jié)果見圖9。協(xié)同克里格插值模型與權(quán)值校正協(xié)同克里格模型的交叉驗(yàn)證結(jié)果見表5。

圖9 交叉驗(yàn)證結(jié)果Fig.9 Results of interactive verification

表5 2 種插值模型交叉驗(yàn)證結(jié)果Table 5 Interactive verification results of two interpolation models %

經(jīng)交叉驗(yàn)證，殘差值大體呈正態(tài)分布，各元素品位實(shí)際值與估計(jì)值的誤差趨于0，表明模型合理。同時(shí)，相對于普通協(xié)同克里格模型而言，在3 570 個組合樣品點(diǎn)上，權(quán)值校正協(xié)同克里格插值模型降低了估值的平均誤差，且其標(biāo)準(zhǔn)差相對較小，表明該模型儲量估算精度具有優(yōu)勢。

4.2.4 儲量估算

針對三維塊段模型的儲量計(jì)算，首先考慮甲瑪?shù)V區(qū)的勘探網(wǎng)度和變差函數(shù)的特征，確定單元塊段的尺寸等參數(shù)如表6 所示。

表6 塊段估值參數(shù)Table 6 Valuation parameters of block segment

甲瑪?shù)V區(qū)總共劃分出153 600 個單元塊段，對每個單元塊段的Cu 品位進(jìn)行估值時(shí)，應(yīng)首先確定待估點(diǎn)的搜索范圍。由于Cu 品位的變差函數(shù)模型以及與輔助變量之間的交叉變差函數(shù)模型均呈球狀模型變化趨勢，因此采用權(quán)值校正協(xié)同克里格法進(jìn)行估值時(shí)，搜索范圍可定義為如圖10 所示三維橢球體，綜合考慮Cu 品位球狀模型參數(shù)和礦區(qū)鉆孔分布間隔等因素，分別取三維橢球體長軸為200 m，次軸為200 m，短軸為30 m。經(jīng)計(jì)算，甲瑪?shù)V區(qū)共探獲夕卡巖型礦體中銅金屬量約209 萬t(含伴生銅金屬量約28 萬t)，基本符合實(shí)際情況。

圖10 搜索范圍Fig.10 Search field

5 結(jié) 語

依據(jù)西藏甲碼銅多金屬礦的鉆孔數(shù)據(jù)資料，在對原始鉆孔數(shù)據(jù)進(jìn)行預(yù)處理的基礎(chǔ)上，通過計(jì)算Cu、Ag的變差函數(shù)及其交叉變差函數(shù)，構(gòu)造協(xié)同克里格方程組，求解得到權(quán)值。采用擬牛頓法對不符合實(shí)際情況的權(quán)值進(jìn)行校正，并對西藏甲瑪?shù)V區(qū)的Cu 儲量進(jìn)行了估算，估算結(jié)果與實(shí)際情況基本一致，可為礦山開采設(shè)計(jì)提供依據(jù)。

［1］ 黎楓佶.西藏墨竹工卡縣甲瑪銅多金屬礦三維模型構(gòu)建和資源量估算［D］.成都:成都理工大學(xué)，2010.

Li Fengji.Constructing 3D Models and Estimating Reserves of Jiama Copper Polymetallic Deposit，Mozhugongka County，Tibet［D］.Chengdu:Chengdu University of Technology，2010.

［2］ 余牛奔，齊文濤，王立歡，等.基于3DMine 軟件的三維地質(zhì)建模及儲量估算——以新疆巴里坤礦區(qū)某井田為例［J］.金屬礦山，2015(3):138-141.

Yu Niuben，Qi Wentao，Wang Lihuan，et al. Three-dimensional geological modeling and reserve estimation based on 3DMine software:taking a well field of Balikun mining area，Xinjiang as an example［J］.Metal Mine，2015(3):138-141.

［3］ 程 勖. 地質(zhì)統(tǒng)計(jì)學(xué)軟件開發(fā)與應(yīng)用［D］. 長春:吉林大學(xué)，2009.

Cheng Xu. Development and Application of Geostatisties Software［D］.Changchun:Jilin University，2009.

［4］ 唐 攀，唐菊興，唐曉倩，等. 傳統(tǒng)方法和地質(zhì)統(tǒng)計(jì)學(xué)在礦產(chǎn)資源儲量分類中的對比分析［J］.金屬礦山，2013(11):107-109.

Tang Pan，Tang Juxing，Tang Xiaoqian，et al. The comparative research of traditional method and geostatistics method in mineral resource/reserve classification［J］.Metal Mine，2013(11):107-109.

［5］ 余 璨，李 峰，張達(dá)兵，等.基于DMINE 的紅龍廠銅礦床地質(zhì)建模與儲量計(jì)算［J］.金屬礦山，2015(2):108-111.

Yu Can，Li Feng，Zhang Dabing，et al.Geological modeling and calculation of the reserves of Honglongchang copper deposit based on DIMINE［J］.Metal Mine，2015(2):108-111.

［6］ 李慶謀.多維分形克里格方法［J］.地球科學(xué)進(jìn)展，2005，20(2):248-256.

Li Qingmou. Multi-fractal Krige interpolation method［J］. Advances in Earth Science，2005，20(2):248-256.

［7］ 肖克炎，張曉華，王全明，等. 應(yīng)用改進(jìn)的克里格法分離重力區(qū)域異常與局部異常［J］.物探化探計(jì)算技術(shù)，1995，17(2):19-25.

Xiao Keyan，Zhang Xiaohua，Wang Quanming et al.Separation of regional field and local field by improved Kriging method［J］.Computing Techniques for Geophysical and Geochemical Exploration，1995，17(2):19-25.

［8］ 高美娟，朱慶忠，張淑華. 利用貝葉斯－克里金估計(jì)技術(shù)進(jìn)行儲層參數(shù)預(yù)測［J］.石油地球物理勘探，1999，34(4):390-397.

Gao Meijuan，Zhu Qingzhong，Zhang Shuhua. Reservoir parameter prediction using Bayes-Kriging estimation technique［J］. Oil Geophyical Prospecting，1999，34(4):390-397.

［9］ 李 君，李少華，毛 平，等. Kriging 插值中條件數(shù)據(jù)點(diǎn)個數(shù)的選擇［J］.斷塊油氣田，2010，17(3):277-279.

Li Jun，Li Shaohua，Mao Ping et al. Choosing the number of conditional data in Kriging interpolation［J］.Fault-Block Oil ＆ Gas Feild 2010，17(3):277-279.

［10］ 汪 保，孫 秦.改進(jìn)的Kriging 模型的可靠度計(jì)算［J］.計(jì)算機(jī)仿真，2011，28(2):113-116.

Wang Bao，Sun Qin. Structural reliability computation based on Kriging model［J］.Computer Simulation，2011，28(2):113-116.

［11］ 孫 卡，翁正平，李章林，等.排序掃描線算法在克里格估值中的應(yīng)用［J］.金屬礦山，2009(7):73-76.

Sun Ka，Weng Zhengping，Li Zhanglin，et al.Application of sequential scanning algorithm in Kriging prediction［J］.Metal Mine，2009(7):73-76.

［12］ 牛文杰，朱大培，陳其明. 滑動領(lǐng)域克里金插值法的改進(jìn)［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2001，13(8):752-756.

Niu Wenjie，Zhu Dapei，Chen Qiming. Improvement of moving neighborhood Kriging method［J］. Journal of Computer-Aided Design ＆ Compute Graphics，2001，13(8):752-756.

［13］ 李少華，王勇標(biāo)，李 君，等. 克里金估值中條帶效應(yīng)的校正［J］.石油地球物理勘探，2012，47(6):979-983.

Li Shaohua，Wang Yongbiao，Li Jun，et al. Correcting the string effect of Kriging［J］. Oil Geophyical Prospecting，2012，47(6):979-983.

［14］ 牛文杰，孟憲海，李吉剛.結(jié)合軟數(shù)據(jù)的同位置協(xié)同克里金估值新方法［J］.煤田地質(zhì)與勘探，2011，39(2):13-17.

Niu Wenjie，Meng Xianhai，Li Jigang. A new estimation method of collocated Co-Kring combined with soft data［J］. Coal Geology ＆Exploration，2011，39(2):13-17.

［15］ 侯景儒，黃競先.地質(zhì)統(tǒng)計(jì)學(xué)的理論與方法［M］.北京:地質(zhì)出版社，1990.

Hou Jingru，Huang Jingxian. Theory and Method of Geostatistics［M］.Beijing:Geological Publishing House，1990.

［16］ 孫仁鐸.空間變異理論及其應(yīng)用［M］.北京:科學(xué)出版社，2005.

Sun Renduo.Spatial Variability Theory and Application［M］. Beijing:Science Press，2005.

［17］ 孫紅泉，康永尚，杜惠芝.實(shí)用地質(zhì)統(tǒng)計(jì)學(xué)程序集［M］.北京:地質(zhì)出版社，1997.

Sun Hongquan，Kang Yongshang，Du Huizhi. Utility of Geostatistical Procedures Set［M］. Beijing:Geological Publishing House，1997.

［18］ 易君君，汪 保，顏倩倩.基于新擬牛頓方程的優(yōu)化算法設(shè)計(jì)及應(yīng)用［J］.寧波工程學(xué)院學(xué)報(bào)，2015，27(1):12-18.

Yi Junjun，Wang Bao，Yan Qianqian.Design and application of optimization algorithm based on new Quasi-Newton equation［J］.Journal of Ningbo University of Technology，2015，27(1):12-18.

［19］ 王小朋，劉翔峰.一種基于病態(tài)問題的修正牛頓法［J］.河南科技大學(xué)學(xué)報(bào):自然科學(xué)版，2015，36(1):86-91.

Wang Xiaopeng，Liu Xiangfeng. A modified Newton method based on ill-conditioning problem［J］.Journal of Henan University of Science and Technology:Natural Science，2015，36(1):86-91.

［20］ 張華軍，趙 金，羅 慧.基于擬牛頓法的同時(shí)擾動隨機(jī)逼近算法［J］.華中科技大學(xué)學(xué)報(bào):自然科學(xué)版，2014，42(9):1-4.

Zhang Huajun，Zhao Jin，Luo Hui. Simultaneous perturbation stochastic approximation algorithm based on Quasi-Newton method［J］.Journal of Huazhong University of Science ＆ Technology:Natural Science Edition，2014，42(9):1-4.