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

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

（金屬礦山高效開(kāi)采與安全教育部重點(diǎn)實(shí)驗(yàn)室）

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

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

（金屬礦山高效開(kāi)采與安全教育部重點(diǎn)實(shí)驗(yàn)室）

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

氧化帶 Kriging IDW DTM擬合

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

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

1 IDW算法與K riging算法概述

空間點(diǎn)插值計(jì)算方法很多，其中常用的有IDW與Kriging插值算法［1］。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 工程應(yīng)用

3.1 工程背景

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

3.2 建立氧化帶DTM

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

表1 鉆孔氧化率數(shù)據(jù)oxidation結(jié)構(gòu)

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

圖1 氧化面初步DTM

圖2 原生礦面DTM

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

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

3.3 模型修正

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

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

圖5 插值后的實(shí)體模型

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

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

3.4 結(jié)果分析

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

圖8 K riging法與融合算法結(jié)果對(duì)比

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

圖9 氧化面界面對(duì)比

4 結(jié) 論

（1）氧化帶化驗(yàn)數(shù)據(jù)屬?gòu)V域離散型數(shù)據(jù)，生成氧化帶表面過(guò)渡生硬，其連接方式導(dǎo)致結(jié)果準(zhǔn)確性差，面積有所限制。

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

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

［1］張伊麗.基于CUDA的泛kriging算法的研究與設(shè)計(jì)［D］.北京：中國(guó)地質(zhì)大學(xué)，2012.

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

［3］丁建華，肖克炎，婁德波，等.大比例尺三維礦產(chǎn)預(yù)測(cè)［J］.地質(zhì)與勘探，2009（6）：729-734.

［4］羅周全，張 保，劉曉明，等.礦體品位和儲(chǔ)量統(tǒng)計(jì)分析的三維可視化方法［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］潘國(guó)成，李秀峰.地質(zhì)統(tǒng)計(jì)學(xué)儲(chǔ)量估算方法在礦業(yè)中應(yīng)用的若干實(shí)際問(wèn)題（一）［J］.國(guó)外金屬礦山，1996（7）：7-12.

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

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

［9］李 健，吳順川，高永濤，等.基于Kriging與Closest Point融合算法的邊坡巖土層界面擬合［J］.北京科技大學(xué)學(xué)報(bào)，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）

*國(guó)家自然科學(xué)基金項(xiàng)目（編號(hào)：51104010），教育部基本科研業(yè)務(wù)費(fèi)項(xiàng)目（編號(hào)：FRF－SD－12－001A），國(guó)家高技術(shù)研究發(fā)展計(jì)劃（863計(jì)劃）項(xiàng)目（編號(hào)：2012AA062201）。

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