李 迪 胡乃聯(lián) 李國清 楊 樺

（金屬礦山高效開采與安全教育部重點實驗室）

基于Kriging與IDW融合算法的氧化帶DTM擬合方法*

李 迪 胡乃聯(lián) 李國清 楊 樺

（金屬礦山高效開采與安全教育部重點實驗室）

氧化礦帶取樣數(shù)據(jù)屬于廣域離散型數(shù)據(jù)，而傳統(tǒng)插值方法對于離散型數(shù)據(jù)不能很好地擬合，誤差較大。為此，基于Kriging法的區(qū)域化變量及距離冪次反比法（IDW），提出了將二者融合的新算法，并以國內某大型露天鉬礦為例進行研究，得出該算法簡單易行，且能大幅提高氧化帶數(shù)字地面模型（DTM）擬合精度的結論。

氧化帶 Kriging IDW DTM擬合

氧化礦是金屬礦床受氧化作用形成的氧化帶中的礦石。礦石受氧化作用后，其礦物組分、結構、構造均產生明顯變化，因而在加工利用時須采用不同于原生礦的方法與工藝流程。地質勘查時，需要通過系統(tǒng)的采樣和物相分析來確定氧化帶的界線及深度。當?shù)V床的氧化礦石占到總儲量的一定比例時，還必須進行選礦試驗，圈定范圍。因而氧化礦量的多少是選礦方法擇取的主要影響因素之一。為選擇合理的選礦方法與工藝流程，通常需要精確計算氧化礦量。從技術經(jīng)濟的角度看，對于氧化量占到一定比例的某些貴金屬礦，宜采用堆浸方式進行選礦處理，且氧化礦量決定了堆浸計量、規(guī)模等相關參數(shù)。因此，精確計算氧化礦量對于最大限度地開發(fā)利用礦產資源具有重要意義。

氧化帶是區(qū)分氧化礦與原生礦的重要依據(jù)。氧化帶的數(shù)字地面模型（DTM）精確程度對氧化礦量的計算具有重大影響，為此，迫切需要對建立精確氧化帶模型的方法進行深入研究。近年來，隨著可視化研究不斷深入，圖形仿真技術普及礦業(yè)領域，礦山數(shù)字化已成必然趨勢。在礦業(yè)軟件中，建立氧化帶DTM，能自動識別區(qū)分氧化礦與原生礦，對后續(xù)計算氧化礦量等相關操作必不可少。然而，現(xiàn)階段用于建立氧化帶DTM的方法較少，且結果并不理想，主要原因如下：其一，氧化率化驗數(shù)據(jù)屬廣域離散型數(shù)據(jù)，生成DTM過度生硬，擬合誤差大，需要進行插值計算；其二，傳統(tǒng)插值方法對于離散型數(shù)據(jù)擬合效果較差，導致插值精度不高；其三，氧化帶插值后與地表、原生礦帶DTM矛盾時如何處理需要進一步討論。為此，本研究提出Kriging與IDW融合算法，并以國內某大型露天鉬礦進行空間插值計算，借助3DMine軟件生成插值擬合模型。

1 IDW算法與K riging算法概述

空間點插值計算方法很多，其中常用的有IDW與Kriging插值算法［1］。

隨各行業(yè)數(shù)字化的發(fā)展，IDW（Inverse Distance Weighted）以其原理簡單、插值過程快速高效的特點，廣泛應用于各領域［2］。國內外諸多軟件引入IDW作為內部插值計算方法，如礦業(yè)軟件3Dmine、MicroMine、Surpac，遙感軟件ERmapper［3］，GIS系列軟件ArcGis、MapGis［4］及統(tǒng)計分析軟件等。

IDW是一種常用而簡便的空間插值方法，它以插值點與樣本點間的距離為權重進行加權平均，離插值點越近的樣本點賦予的權重越大.該方法對于未采樣點的估計值是未采樣點屬性真值的無偏估計。然而，IDW法單單只考慮了距離與待估值之間的關系，在本研究所要插值計算的氧化帶模型中并沒有考慮到空間地質因素，即氧化帶成因所造成的空間相關性，因而該方法具有一定的弊端。為提高估值精度，人們廣泛應用Kriging進行估值。

Kriging法是以礦化空間結構理論為基礎，區(qū)域化變量為核心，變異函數(shù)為工具的數(shù)學地質方法，是地質統(tǒng)計學中的核心方法。地質統(tǒng)計學發(fā)展近60 a，在國內外各學科領域有著廣泛的應用。在礦業(yè)領域，地質統(tǒng)計學不僅廣泛應用于資源儲量計算，在其他礦業(yè)活動，如異常評價［5］、找礦勘探［6］、礦體圈定［7］采礦設計［8］及地學科研［9］等方面也具有明顯的優(yōu)越性。

Kriging相比于IDW最突出的優(yōu)勢是區(qū)域化變量理論，區(qū)域化變量理論將地質因素加入到插值計算中，使得插值運算不再是純數(shù)學問題，而是與空間地質相關的地質統(tǒng)計學問題。本研究的插值對象是氧化率取樣點標高，即Z值，以Z（xu，xv，xw）作為區(qū)域化變量。區(qū)域化變量能夠反映地質變量的結構性及隨機性。當空間點x位置確定后，該點標高Z（x）為隨機變量，空間兩不同點x、x＋h處標高Z（x）、Z（x＋h）具有某種程度的自相關性，這種自相關性依賴于分隔兩點的距離h和地質構造、礦化規(guī)律，體現(xiàn)了氧化礦成礦過程的某種聯(lián)系性，即區(qū)域化變量的結構性。

Kriging法全面應用數(shù)理統(tǒng)計理論，遵循無偏估計和最優(yōu)估計原則，考慮了地質因素對差值計算的影響，因而估算精度相對較高。然而，Kriging的最大缺陷是其平滑效應［10］，這種平滑效應使得氧化面高程均化，嚴重影響插值精度。為提高模型擬合的精確度，還需進一步改進算法。

針對2種算法各有優(yōu)缺點的現(xiàn)狀，本研究考慮將2種算法融合。

2 基于K riging與IDW融合算法的數(shù)學模型

基于上述2種插值方法的優(yōu)缺點，本研究提出Kriging與IDW融合算法，將地質統(tǒng)計學中的核心優(yōu)勢——區(qū)域化變量理論引入IDW。本研究的插值對象雖為氧化帶標高，但能體現(xiàn)氧化面成礦規(guī)律的是其氧化率屬性，因而Kriging中的變異函數(shù)擬合對象是氧化率屬性。首先通過Kriging方法中針對氧化率屬性的變異函數(shù)擬合找到插值點間的空間地質關聯(lián)，Kriging中的空間地質關聯(lián)集中體現(xiàn)在插值權重等值面構成的搜索橢球體上，通過Kriging變異函數(shù)的擬合找到搜索橢球體定義參數(shù)，包括橢球體搜索半徑，橢球體傾角、方位角、傾伏角，以及主軸、半主軸、短軸之間的比率等。最后在IDW中對氧化面標高屬性的插值過程中根據(jù)Kriging求得的搜索橢球體作為插值中搜索樣品對的方式，繼而進行IDW插值計算。這種融合算法能夠有效避免Kriging的平滑效應與IDW中缺乏地質成因關聯(lián)的不足。

IDW法用于插值的基本公式為

式中，Z*（B）為待估點屬性值，這里指待估點氧化帶標高；Z（xi）為已知采樣點標高；λi為已知點權重。依據(jù)IDW法的基本思想，確定權重λi的方法為

式中，di為待估點與已知點間距離；k作為di的冪指數(shù)，其取值由具體情況確定，通?？梢匀?，2，3等整數(shù)。

Kriging方法中使用變異函數(shù)來描述區(qū)域化變量的空間結構性變化和隨機性變化，通過變異函數(shù)擬合尋找橢球體搜索參數(shù)。氧化面取樣點氧化率R（x）在點x、x＋h處增量方差的一半，即R（x）在該方向上的半變異函數(shù)。

在實踐中，把由有限實測樣品值構成的變異函數(shù)稱之為實驗變異函數(shù)，記為g*（h）。

式中，n（h）為相距h的樣品對數(shù)；R（xi）、R（xi＋h）分別為xi、xi＋h處的氧化率化驗值。

Kriging方程組求解思想是使估計方差最小［11］，經(jīng)過Kriging法求解尋得最優(yōu)擬合結果，得到相應的搜索橢球體基本參數(shù)，并將該系列參數(shù)帶入到IDW中，確定IDW中搜索樣品對的方式，最后由式（1）、式（2）、式（3）求得插值結果。

3 工程應用

3.1 工程背景

某露天鉬礦采礦權面積4.6 km2，開采深度640～0 m標高，為目前國內最大單一鉬礦床。根據(jù)地質勘探報告，得到如下氧化帶信息：該礦將勘探工作中氧化鉬含量與全鉬含量的比值（氧化率）≥10%的鉆孔樣品定義為氧化帶。礦體氧化帶最大深度51.93 m，平均16.82 m；氧化礦石最厚17.50 m，平均厚4.28 m；在松散蓋層以下，氧化帶平均氧化率最高76.45%，最低10.02%，平均37.42%；礦化巖石氧化率，一般由地表向深部逐漸降低。在地表，氧化率最高達100%，向深部最低為10.2%。局部處在后期斷裂構造帶上的氧化礦或礦化巖石，氧化率呈現(xiàn)高低相間分布，變化較大。簡言之，礦體氧化帶不甚發(fā)育，一般順地形起伏而變化，呈似層狀，但在后期斷裂構造帶附近，不僅氧化深度較大，且氧化率值變化無明顯規(guī)律。

3.2 建立氧化帶DTM

氧化帶上表面接近地表，區(qū)分氧化礦與原生礦的主要依據(jù)是氧化帶下表面（即下氧化面）。這里主要研究下氧化面DTM。本次研究共收集氧化率數(shù)據(jù)245個。首先將氧化率化驗數(shù)據(jù)導入數(shù)據(jù)庫，在常規(guī)數(shù)據(jù)庫基礎上新建oxidation表。oxidation表結構如表1。

表1 鉆孔氧化率數(shù)據(jù)oxidation結構

根據(jù)化驗點已知氧化率生成氧化面與原生礦面DTM，即分別由氧化率大于10%與小于10%的點生成，如圖1、圖2所示。

圖1 氧化面初步DTM

圖2 原生礦面DTM

氧化率大于10%的化驗取樣點被認為是氧化帶的組成部分，即插值計算中的已知點。針對已知點采用本研究提出的Kriging與IDW融合插值法。首先根據(jù)氧化率組合樣建立變異函數(shù)模型，繼而進行結構分析擬合，圖3顯示了主軸變異函數(shù)的擬合狀況，半主軸、短軸擬合過程類似，這里不再贅述。最后得到基臺值90.584 01，變程284.894，方位角139°，傾角10°，傾伏角0°，主軸與半主軸、短軸之比分別為1.27、1.00，橢球體如圖4所示。最后將得到的搜索橢球體應用到IDW中，作為IDW樣品對的搜索方式進行插值，由插值點與原樣品點共同參與生成下氧化面DTM，如圖5所示。

圖3 變異函數(shù)主軸擬合過程

3.3 模型修正

由于氧化面不能穿過地表，并且不能與化驗小于10%的樣品點矛盾，這里調入地表DTM與原生礦面DTM，對插值后的氧化帶模型進行修剪、校正，其過程如圖6、圖7所示。

圖4 K riging擬合的各向異性橢球體

圖5 插值后的實體模型

圖6 氧化面與地表相交處理

圖7 氧化面與原生礦面相交處理

3.4 結果分析

為進一步驗證融合算法的準確性，研究使用Kriging法對化驗點進行插值計算，并針對原化驗點將插值結果與融合算法進行對比（如圖8）：Kriging插值具有明顯的光滑效應，其結果均化明顯，對原取樣點變化趨勢的擬合效果較差，而融合算法所得插值結果更加服從原始數(shù)據(jù)，并較為準確地擬合出其趨勢效果。

圖8 K riging法與融合算法結果對比

圖9反映了插值前后氧化面的變化，經(jīng)對比發(fā)現(xiàn)插值計算得到的氧化帶DTM相較于插值前的模型更加平滑、精確，貼近實際，一定程度擴展了氧化面面積，能夠有效切割礦體，并由地表模型、原生礦面模型加以校正，增加模型擬合精度，為求取氧化礦石量提供了準確保障，對該礦山后續(xù)選礦方式及礦石堆存方式的擇取具指導性意義。

圖9 氧化面界面對比

4 結 論

（1）氧化帶化驗數(shù)據(jù)屬廣域離散型數(shù)據(jù)，生成氧化帶表面過渡生硬，其連接方式導致結果準確性差，面積有所限制。

（2）將Kriging方法與IDW法融合，針對廣域離散型數(shù)據(jù)進行插值計算，能夠有效避免Kriging的平滑效應與IDW中缺乏空間地質因素的不足，提高氧化帶表面擬合精度。

（3）該方法借助3DMine礦業(yè)軟件插值計算，操作相對便捷。

［1］張伊麗.基于CUDA的泛kriging算法的研究與設計［D］.北京：中國地質大學，2012.

［2］李章林，王 平，張夏林.距離冪次反比法的改進與應用［J］.金屬礦山，2008（4）：88-92.

［3］丁建華，肖克炎，婁德波，等.大比例尺三維礦產預測［J］.地質與勘探，2009（6）：729-734.

［4］羅周全，張 保，劉曉明，等.礦體品位和儲量統(tǒng)計分析的三維可視化方法［J］.有色金屬，2008（5）：23-27.

［5］Kringe D G.Some basic considerations in the application of geostatistics to valuation of ore in south African gold mines［J］.The South African Institute of Mining and Metallurgy，1976，76（9）：383-391.

［6］潘國成，李秀峰.地質統(tǒng)計學儲量估算方法在礦業(yè)中應用的若干實際問題（一）［J］.國外金屬礦山，1996（7）：7-12.

［7］潘國成，李秀峰.地質統(tǒng)計學儲量估算方法在礦業(yè)中應用的若干實際問題（二）［J］，國外金屬礦山，1996（8）：5-9.

［8］李曉暉，袁 峰，賈 蔡，等.基于地統(tǒng)計學插值方法的局部奇異性指數(shù)計算比較研究［J］.地理科學，2012（2）：136-142.

［9］李 健，吳順川，高永濤，等.基于Kriging與Closest Point融合算法的邊坡巖土層界面擬合［J］.北京科技大學學報，2012，34（5）：500-505.

Oxidation Belt DTM Fitting M ethod Based on K riging and IDW Nesting Algorithm

Li Di Hu Nailian Li Guoqing Yang Hua
（State Key Laboratory of High-Efficient Mining and Safety of Metal Mines，Ministry of Education）

Sampling data of oxidation belt iswide and discrete，and the traditional interpolationmethods are not fitting well with the discrete data.This resourch，based on the regionalized variable and variation function of Kriging along with Inverse DistanceWeighted（IDW），put forward a new algorithm nesting the twomethods.Based on investigation in a largemolybdenum deposit at domestic，the conclusion is that this algorithm is simple and can greatly improve the fitting accuracy of the digital terrain model（DTM）of oxidation belt.

Oxidation belt，Kriging，IDW，DTM fitting

2013-08-04）

*國家自然科學基金項目（編號：51104010），教育部基本科研業(yè)務費項目（編號：FRF－SD－12－001A），國家高技術研究發(fā)展計劃（863計劃）項目（編號：2012AA062201）。

李 迪（1988—），女，碩士研究生，100083北京市海淀區(qū)學院路30號。