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

李懷良 庹先國， 蔣 鑫

(1.核廢物與環境安全國防重點學科實驗室，綿陽 四川621010;2.地質災害防治與地質環境保護國家重點實驗室，成都 四川610059)

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

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

1 協同克里格法

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

式中，λk，i為第k 個區域化變量中的第i 個樣本點的權值，由最小化方差求出，即

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

該算法主要實現流程見圖1。

2 鉆孔條帶狀效應

西藏甲瑪銅多金屬礦最主要的數據有地面鉆孔數據、鉆孔化驗數據、小體重數據等。原始數據為截至2012 年甲瑪礦區勘探工程數據，包括58 個鉆孔共15 187 件樣品的數據。若某個鉆孔遇到斷層，鉆孔數據便會出現邊界截斷現象，如圖2 所示。

如圖3 所示，在進行協同克里格插值時首先需要確定一個搜索范圍，搜索范圍在三維空間內通常為一個橢球體的邊界截斷。圖1、圖2 所示的2 類數據可稱為鉆孔的有限條帶狀數據。

圖1 協同克立格法實現流程Fig.1 Implementation process of Co-Kriging interpolation algorithm

圖2 數據被斷層邊界截斷Fig.2 Drilling data truncated by boundary faults

圖3 數據被搜索鄰域截斷Fig.3 Drilling data truncated by search field

以Cu 和Ag 的協同區域化為例，對其中的一列有限的條帶狀數據進行協同克里格插值計算，其中變差函數以及2 種元素的交叉變差函數均采用球狀模型計算。其中Cu 塊金效應為0 的協同克里格權值如圖4 所示。

圖4 協同克里格插值引起的條帶狀效應Fig.4 Banned effect caused by Co-Kriging interpolation

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

3 權值校正

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

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

4 處理結果

4.1 數據預處理

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

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

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

表1 Cu 特異值處理結果Table 1 Results of singular value of Cu

采用權值校正協同克里格法估算資源量時，組合樣品長度設定為6 m。組合樣品中Cu 品位特征值的統計結果見表2。

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

4.2 協同克里格法估算儲量

4.2.1 變差函數及結構分析

對組合樣品進行統計分析，結果見表3。由表3可知，組合樣品中Cu 品位與Ag 品位相關性最大，因此以Ag 品位作為輔助變量對Cu 的品位進行估值。

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

表3 組合樣本統計分析Table 3 Composite samples analysis results

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

由表4 可知，2 種元素在各個方向上的變差函數參數相似，所有變差函數均為各向同性的。對Cu、Ag品位項的3 個軸方向上的變差函數進行結構套合，得到Cu、Ag 品位的變差函數模型

式中，h 為滯后距，m。

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

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

4.2.2 權值校正

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

圖8 權值校正結果Fig.8 Results of weight calibration

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

4.2.3 算法模型驗證

對甲瑪礦區化學樣品中Cu 品位估值情況進行交叉驗證，結果見圖9。協同克里格插值模型與權值校正協同克里格模型的交叉驗證結果見表5。

圖9 交叉驗證結果Fig.9 Results of interactive verification

表5 2 種插值模型交叉驗證結果Table 5 Interactive verification results of two interpolation models %

經交叉驗證，殘差值大體呈正態分布，各元素品位實際值與估計值的誤差趨于0，表明模型合理。同時，相對于普通協同克里格模型而言，在3 570 個組合樣品點上，權值校正協同克里格插值模型降低了估值的平均誤差，且其標準差相對較小，表明該模型儲量估算精度具有優勢。

4.2.4 儲量估算

針對三維塊段模型的儲量計算，首先考慮甲瑪礦區的勘探網度和變差函數的特征，確定單元塊段的尺寸等參數如表6 所示。

表6 塊段估值參數Table 6 Valuation parameters of block segment

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

圖10 搜索范圍Fig.10 Search field

5 結 語

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

［1］ 黎楓佶.西藏墨竹工卡縣甲瑪銅多金屬礦三維模型構建和資源量估算［D］.成都:成都理工大學，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 軟件的三維地質建模及儲量估算——以新疆巴里坤礦區某井田為例［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］ 程 勖. 地質統計學軟件開發與應用［D］. 長春:吉林大學，2009.

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

［4］ 唐 攀，唐菊興，唐曉倩，等. 傳統方法和地質統計學在礦產資源儲量分類中的對比分析［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］ 余 璨，李 峰，張達兵，等.基于DMINE 的紅龍廠銅礦床地質建模與儲量計算［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］.地球科學進展，2005，20(2):248-256.

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

［7］ 肖克炎，張曉華，王全明，等. 應用改進的克里格法分離重力區域異常與局部異常［J］.物探化探計算技術，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］.石油地球物理勘探，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 插值中條件數據點個數的選擇［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］ 汪 保，孫 秦.改進的Kriging 模型的可靠度計算［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］ 孫 卡，翁正平，李章林，等.排序掃描線算法在克里格估值中的應用［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］ 牛文杰，朱大培，陳其明. 滑動領域克里金插值法的改進［J］.計算機輔助設計與圖形學學報，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］ 李少華，王勇標，李 君，等. 克里金估值中條帶效應的校正［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］ 牛文杰，孟憲海，李吉剛.結合軟數據的同位置協同克里金估值新方法［J］.煤田地質與勘探，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］ 侯景儒，黃競先.地質統計學的理論與方法［M］.北京:地質出版社，1990.

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

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

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

［17］ 孫紅泉，康永尚，杜惠芝.實用地質統計學程序集［M］.北京:地質出版社，1997.

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

［18］ 易君君，汪 保，顏倩倩.基于新擬牛頓方程的優化算法設計及應用［J］.寧波工程學院學報，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］ 王小朋，劉翔峰.一種基于病態問題的修正牛頓法［J］.河南科技大學學報:自然科學版，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］ 張華軍，趙 金，羅 慧.基于擬牛頓法的同時擾動隨機逼近算法［J］.華中科技大學學報:自然科學版，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.