西藏甲瑪銅多金屬礦床儲量的協同克立格估值

蔣 鑫，庹先國，，柳炳琦，李懷良

（1．成都理工大學地質災害防治與地質環境保護國家重點實驗室，成都610059；2．西南科技大學核廢物與環境安全國防重點學科實驗室，綿陽621010）

西藏甲瑪銅多金屬礦床儲量的協同克立格估值

蔣 鑫1，庹先國1，2，柳炳琦1，李懷良2

（1．成都理工大學地質災害防治與地質環境保護國家重點實驗室，成都610059；2．西南科技大學核廢物與環境安全國防重點學科實驗室，綿陽621010）

西藏甲瑪銅多金屬礦礦區元素分布復雜，傳統的地質統計學方法對其進行儲量估算時忽略了多金屬的相互影響，因此為了反映出不同金屬元素的空間變異情況，采用協同克里格法對該礦區的金屬元素儲量進行估算。這里首先介紹了協同克里格算法的相關理論和相關技術，然后以此為基礎，對不同方向上的空間變差函數進行結構套合的優化，并對協同克里格方程組進行降維處理。最后以2012年甲瑪礦區勘探工程的數據為例，以Cu為主區域化變量，以Ag為協同區域化變量，計算了各自的實驗變差函數和交差實驗變差函數，分別進行協同克里格法插值和普通克里格法插值。交叉驗證結果表明，協同克里格估值的標準差為0．6477，在儲量計算上面精度更高，并能廣泛應用于西藏甲瑪銅多金屬礦的地質屬性、儲量估算等空間數據建模。

儲量估算；協同克里格；空間變異；結構套合

0 引言

西藏墨竹工卡縣甲瑪夕卡巖－角巖型銅多金屬礦床是目前我國同類礦床礦石品質最好、成礦元素多，探明的資源量已經達到超大型規模的礦床［1］。目前對于該地區的儲量估算已經采用的是傳統斷面法［4］、SD法［2］和地質統計學方法［3］，其中地質統計學方法不但能給出結果的估算精度而且是無偏且估計方差最優的，這些優勢恰好是傳統方法所缺少

的［5］。

西藏甲瑪銅多金屬礦的區域內地層比較復雜，地層曲面變化比較大，甚至可能存在地層曲面相交和錯斷的情況，并且多種金屬之間會出現協同區域化［11－14］的現象，即多種金屬的儲量分布會相互影響。在這樣一種情況下，就需要我們解決如何尋找一種合適的方法消除多金屬相關關系的影響。這里采用協同克里格法對該礦區的空間變異進行分析，并且對礦區的銅元素儲量估算進行實例驗證。驗證結果發現，協同克里格法的估算結果更符合實際情況，對深化礦業和石油行業開發地質隨機建模理論和有效開采剩余礦產石油資源，具有重要的理論意義和現實意義。

1 協同克里格算法原理

在西藏甲瑪銅多金屬礦的觀測數據中，常遇到多個變量協同區域化的現象，如銀、銅、鉛、鋅等幾個變量的協同區域化，當某一個變量的取樣量不足以獲得所需精度的估計量，而其他變量卻有較充足的取樣量時，研究這個變量與其他變量間的空間相關關系，借助其他變量的樣品信息用協同克里格法就可以提高對這個變量的估計精度。設有K個二階平穩的變量構成協同區域化變量｛Zk（x），k＝1，2，…，K｝，其數學期望為：

設k0為k＝1、2、…、K中某一特定值，即設Zk0（x）為我們要估計的主變量，設要估計的是中心點在x0處的承載Vk0（x0）上Zk0（x）的平均值：

有效數據集是｛Zak，ak＝1，2，…，nk｝，k＝1、2、…、K，而Zak是確定在承載Vak上的平均值，即：

要得到待估樣本點的無偏最優估計值，必須滿足最小方差無偏估計的要求：

為使估計方差達到最小，則要求權重系數要滿足協同克立格方程組，即：

其中：?ak＝1、2、…、nk；k＝1、2、…、K，這個方程組的未知數是個權系數λak及K個拉格朗日參數μk；共有＋K個未知數，同時也有＋K個線性方程構成這個方程組。

2 協同克里格算法的改進與實現

2．1 交叉變差函數的建立

如式（1）所示的二階平穩的協同區域化變量，其互協方差函數即交叉變差函數的定義為：

對交叉變差函數和變差函數進行離散計算時，需要用球狀模型擬合理論變差函數。球狀函數模型的一般公式為：

其中：C0為塊金效應，它表示h很小時兩個樣品點之間品位的變化；a為變程，當h＞a時，任意兩樣品點就不在具有相關性。以球狀函數為基礎，通過計算可以得出交叉變差函數的理論模型。

2．2 求解協同克里格方程組

將式（4）轉化成矩陣形式為：

普通克里格方程組只要在二階平穩條件下，總是可以用協方差函數表示，又可以用變差函數表示。而協同克里格方程組不同，二階平穩條件下，只有當Ck′k（h）＝Ckk＇（h），從而λk′k（h）＝Ck′k（0）－Ck′k（h）成立時，才可用交叉協方差函數表示。

2．3 優化協同克里格方程組

這里主要從兩個方面來優化［6－10］協同克立格算法，①變差函數；②從求解協同克立格方程組。因為在實際問題中，各個方向的變差函數是不一樣的，就需要研究不同方向上的結構套合問題，主要是通過計算不同方向的變差函數值后，再用球狀模型來擬合實驗變差函數，最后進行結構套合。改進后變差函數流程圖如圖1所示。

在求解協同克立格方程組系數時，協同克立格方程組的維數與樣品數據個數有直接關系，也就是說樣品數據越多，方程組維數越大。這就需要對協同克立格方程組進行降維處理或者減少協同克立格方程組的數目，協同克立格方程組的矩陣形式為式（8）。

針對研究數據樣本，樣本數據構形相同，協同克立格矩陣［K］相同，所以只需計算一次［K］－1，假設每次估值計算時，待估塊段大小相等，幾何特性相等，那么［M］（M表示協方差函數向量）相同。這時，就只需要求解一次協同克立格方程組，所得到的權重系數［λ］＝［K］－1［M］，就適合所有的待估塊，就能進行比較快速的計算。權值校正的協同克里格算法流程如圖2所示。

圖1 變差函數流程圖Fig．1 Flow chart of variogram

圖2 協同克立格插值算法流程圖Fig．2 Flow chart of Cokriging

圖3 地表及鉆孔點三維模型圖Fig．3 3Dmodel diagram of surface and drill

3 應用實例

作者以西藏自治區墨竹工卡縣甲碼銅多金屬礦床儲量估算為例。原始數據為截止2012年甲瑪礦區勘探工程的數據。包括58個鉆孔，共計15 187件化學樣。每個鉆孔均有完整的測斜數據共568個。礦床地表模型及鉆孔點的三維模型如圖3所示。

3．1 數據預處理

對礦床儲量的協同克里格法估值，首先需要對鉆孔數據進行特異值處理與樣品組合處理，以Cu元素的計算進行實例說明。首先根據Cu的品位數據，在正態概率紙上作累積頻率曲線圖（圖4）。

圖4表明Cu品位的累積頻率符合正態分布，并且存在少量異常偏大的離散點，我們才認為它們是真正的特異值。對這種異常的特高品位點用上截品位代替，取第一個發生偏離的品位值為上截品位，將所有高于上截品位的品位數據全部用上截品位值代替。代替前、后的統計分析見表1。原始數據樣品長度各不相同，需要進行樣品組合，組合后的Cu元素特征值的統計結果見表2。

圖4 Cu元素累積頻率曲線圖Fig．4 Cumulative frequency curves of Cu

表1 Cu元素特異值處理結果Tab．1 The results of specific values of Cu

3．2 協同克里格法估算儲量

3．2．1 變差函數及結構分析

對組合樣本進行統計分析，分析結果如表3所示，組合樣品中Cu的品位數據與Ag的品位數據相關程度最大，所以以Ag元素作為輔助變量對Cu元素的品位進行估值。

進行協同克里格法插值的第一步是得到Cu品位和Ag品位的變差函數，以及兩者的交叉變差函數。把所有鉆孔點兩兩組合成點對，計算出其距離值，再以球狀模型擬合成理論實驗變差函數。Cu品位項和Ag品位項的三個變差函數軸的方向定位及其計算結果見表4。

表3 組合樣本統計分析表Tab．3 Combined sample analysis table

表4 Cu和Ag元素變差函數參數表Tab．4 Variogram Parameter of Cu and Ag

圖5 變差函數圖Fig．5 Diagram of variogram

兩種元素品位的變差函數圖見圖5。

從圖5可以看出，兩種元素在各個方向上的變差函數比較相似，因此所有變差函數都是各向同性的。對Cu品位項和Ag品位項的三個軸方向上的變差函數進行結構套和，套和結果如下：

Cu品位的變差函數模型為：

Ag品位的變差函數模型為：

以這兩種品位為基礎計算它們的交叉變差函數，同樣交叉變差函數的變化趨勢滿足球狀模型，所以用球狀模型擬合計算出交叉變差函數的基臺值為0．58，拱高為0．42，變程為55m，其變差函數圖如圖6所示。

圖6 交叉變差函數圖Fig．6 Diagram of interactive variogram

3．2．2 算法模型驗證

為了檢驗協同克里格插值的合理性，采用交叉驗證方法對鉆孔點依次求解。其基本原理是從信息樣集合中臨時抽取一個點作為已知值，利用協同克里格插值模型去估算這個已知值，然后對計算結果和真實值的差值進行統計分析。以Cu元素的品位估值情況進行交叉驗證為例，驗證結果見圖7。

同理用普通克里格插值模型和已知樣品去估算已知值，與協同克里格模型相比結果見表5。

表5 協同克里格驗證結果表Tab．5 Result of Cokriging

從圖7可以看出，Cu元素品位的殘差值大部分在3以內，與實際鉆孔值相當接近，表明模型確定合理，能用于下一步的儲量計算。與普通克里格模型相比發現，在3 570個組合樣品點上，協同克里格插值模型減少了估值的平均誤差，并且普通克里格插值模型的標準差總是高于協同克里格插值模型的標準差，這表明協同克里格模型在儲量計算上面精度更高。

3．2．3 儲量計算

針對三維塊段模型的儲量計算，運用變差函數的特征確定單元塊段的尺寸參數（表6）。

圖7 交叉驗證結果圖Fig．7 Result of interactive verification

表6 塊段估值參數表Tab．6 Valuation parameter block segment table

由表6可以看出，甲瑪礦區總共劃分出153 600個單元塊段，對每個單元塊段的Cu品位進行估值時，需要確定對于該估值點的搜索范圍。由于Cu品位的變差函數模型以及與輔助變量之間的交叉變差函數模型都呈球狀模型變化趨勢，而球狀模型的特點是變差函數值隨著距離變化到變程時，變差函數值就穩定為基臺值，所以變程大小直接影響待估單元塊段的搜索區域。進行協同克里格估值時，搜索范圍定義為具有一個長軸和兩個短軸的三維橢球體，綜合Cu品位球狀模型參數和礦區鉆孔分布間隔的因素，搜索的三維橢球體設長軸為200m，次軸為200m，短軸為30m，在整個橢球體內部的單元塊段都要參與估值計算（圖8）。

通過計算，2012年西藏甲瑪礦區的勘探區域共探獲銅金屬量為229×104T（只考慮銀元素為輔助變量）。

圖8 搜索范圍圖Fig．8 Diagram of search field

4 結語

對于西藏自治區墨竹工卡縣甲碼銅多金屬礦床儲量估算，當用傳統普通克里格方法進行插值計算時，沒有考慮不同金屬元素的相關關系，插值精度不高。基于此，作者采用協同克里格法，并對協同克里格法進行優化，能快速地計算西藏甲碼銅多金屬礦的資源儲量。作者首先對鉆孔數據資料進行了統計分析和預處理，然后計算Cu元素和Ag元素的變差函數及其交叉變差函數以構造協同克里格方程組，通過求解方程組得出權值系數。最后對甲瑪礦區的Cu元素儲量進行了估算，得到Cu金屬量為229× 104T。同時通過方差驗證，與普通克里格法的驗證結果進行對比，結果發現，協同克里格法因為反映了不同金屬元素品位間的空間變異情況，可靠性更高，能為礦山開采設計提供科學、可靠的依據。

［1］黎楓佶．西藏墨竹工卡縣甲瑪銅多金屬礦三維模型構建和資源量估算［D］．成都：成都理工大學，2010．

LI F J．Constructing 3Dmodels and estimating reserves of Jiama polymetallic copper deposit，MoZhugongka county，Tibet［D］．Chengdu：Chengdu University of Technology，2010．（In Chinese）

［2］柳炳利．基于SD法的固體礦產資源研究［D］．成都：成都理工大學，2009．

LIU B L．Research of solid mineral resources based on the SD method［D］．Chengdu：Chengdu University of Technology，2009．（In Chinese）

［3］程勖．地質統計學軟件開發與應用［D］．長春：吉林大學，2009．

CHENG X．Geostatisties sottware development and application［D］．Changchun：Jilin University，2009．（In Chinese）

［4］施俊法，唐金榮，周平，等．世界找礦模型與礦產勘查［M］．北京：地質出版社，2010．

SHI J F，TANG J R，ZHOU P，et al．The world model of prospecting and exploration of mineral resources［M］．Beijing：geological publishing house，2010．（In Chinese）

［5］楊學善，秦德先，鄧明國．地學空間數據庫建設快捷方法研究［J］．煤田地質與勘探，2005，33（5）：16－19．

YANG X SH，QIN D X，DENG M G．Shortcut methods of spatial database constru－ction for geosciences［J］．Coal Geology ＆Exploration，2005，33（5）：16－19．（In Chinese）

［6］李慶謀．多維分形克里格方法［J］．地球科學進展，2005，20（2）：248－256．

LI Q M．Multi－fractal Krige interpolation method［J］．Advances in Earth Science，2005，20（2）：248－256．（In Chinese）

［7］肖克炎，張曉華，王全明，等．應用改進的克里格法分離重力區域異常與局部異常［J］．物探化探計算技術，1995，17（2）：19－25．

XIAO K Y，ZHANG X H，WANG Q M，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．（In Chinese）

［8］高美娟，朱慶忠，張淑華．利用貝葉斯－克里金估計技術進行儲層參數預測［J］．石油地球物理勘探，1999，34（4）：390－397．

GAO M J，ZHU Q ZH，ZHANG SH H．Reservoir parameter prediction using Bayes－Kriging estimation technique［J］．OGP，1999，34（4）：390－397．（In Chinese）

［9］汪保，孫秦．改進的Kriging模型的可靠度計算［J］．計算機仿真，2011，28（2）：113－116．

WANG B，SUN Q．Structural reliability computation based on Kriging model［J］．Computer Simulation，2011，28（2）：113－116．（In Chinese）

［10］牛文杰，孟憲海，李吉剛．結合軟數據的同位置協同克里金估值新方法［J］．煤田地質與勘探，2011，39（2）：13－17．

NIU W J，MENG X H，LI J G．A new estimation method of collocated CoKing combined with soft data［J］．Coal Geology ＆Exploration，2011，39（2）：13－17．（In Chinese）

［11］孫紅泉．地質統計學及其應用［M］．北京：中國礦業大學出版社，1990．

SUN H Q．Geostatistics and its application［M］．Beijing：China University of Mining Press，1990．（In Chinese）

［12］侯景儒，黃競先．地質統計學的理論與方法［M］．北京：地質出版社，1990．

HOU J R，HUANG J X．The theory and method of Geostatistics［M］．Beijing：geological publishing house，1990．（In Chinese）

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

SUN R D．the Spatial Variability Theory and Application［M］．Beijing：Science Press，2005．（In Chinese）

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

SUN H Q，KANG Y SH，DU H ZH．The utility of geostatistical procedures set［M］．Beijin：geological publishing house，1997．（In Chinese）

Estimation of reserves based on Cokriging in Jiama polymetallic copper deposit，Tibet

JIANG Xin1，TUO Xian-guo1，2，LIU Bing－qi1，LI Huai-liang2
（1．State Key Laboratory of Geohazard Prevention and Geoenvironment Protection，Chengdu University of Technology，Chengdu 610059，China；2．Key Laboratory for Radioactive Waste and Environmental Security，Southwest University of Science and Technology，Mianyang 621010，China）

The distribution of mining area element is very complex in Jiama polymetallic copper deposit，Tibet．Traditional geostatistical method ignored the phenomenon of mutual influence more metal on the reserves estimation．In order to reflect different metal elements of space mutation，metal reserves of the mine area is estimated based on the cokriging method in this paper．The paper first introduces the relevant theory and technology of the Cokriging algorithm，then optimizes spatial variation function in different directions to structure of nested，and decreases the dimension of cokriging equations．Finally，by taking an example of the data of exploration engineering for Jiama mining area in 2012，with Cu as the main regionalized variables and with Ag as collaborative regionalized variables，this paper calculate the experimental variogram and the interactive experimental variogram，respectively，Cokriging interpolation and the ordinary kriging interpolation．Cross validation results point out that the standard deviation of Cokriging is 0．6477，and the Cokriging is higher in the calculation of reserves above precision．The Cokriging can be widely used in the geological properties，reserve estimation and spatial data modeling of Jiama polymetallic copper deposit，Tibet．

estimating reserves；Cokriging；spatial variation；structure of nested

P 632

A

10．3969／j．issn．1001-1749．2015．03．16

1001-1749（2015）03-0372-07

2014-07-01 改回日期：2014-08-13

國家自然科學基金重大科研儀器設備研制專項（41227802）；國家杰出青年基金（41025015）

蔣鑫（1991－），男，碩士，研究方向為計算機應用，現主要從事勘探設備軟件開發，E-mail：jiangxin＿cdut＠163．com。