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

李 迪 胡乃聯 李國清 楊 樺

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

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

李 迪 胡乃聯 李國清 楊 樺

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

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

氧化帶 Kriging IDW DTM擬合

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

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

1 IDW算法與K riging算法概述

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

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

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

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

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

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

針對2種算法各有優缺點的現狀，本研究考慮將2種算法融合。

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

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

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

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

式中，di為待估點與已知點間距離；k作為di的冪指數，其取值由具體情況確定，通常可以取1，2，3等整數。

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

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

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

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

3 工程應用

3.1 工程背景

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

3.2 建立氧化帶DTM

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

表1 鉆孔氧化率數據oxidation結構

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

圖1 氧化面初步DTM

圖2 原生礦面DTM

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

圖3 變異函數主軸擬合過程

3.3 模型修正

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

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

圖5 插值后的實體模型

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

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

3.4 結果分析

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

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

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

圖9 氧化面界面對比

4 結 論

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

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

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

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

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

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

［4］羅周全，張 保，劉曉明，等.礦體品位和儲量統計分析的三維可視化方法［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］潘國成，李秀峰.地質統計學儲量估算方法在礦業中應用的若干實際問題（一）［J］.國外金屬礦山，1996（7）：7-12.

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

［8］李曉暉，袁 峰，賈 蔡，等.基于地統計學插值方法的局部奇異性指數計算比較研究［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），教育部基本科研業務費項目（編號：FRF－SD－12－001A），國家高技術研究發展計劃（863計劃）項目（編號：2012AA062201）。

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