基于多維分形模型與指示克里格方法的地球化學(xué)異常識別研究

李 曉 暉，袁 峰，賈 蔡，張 明 明，周 濤 發(fā)

基于多維分形模型與指示克里格方法的地球化學(xué)異常識別研究

李 曉 暉，袁 峰，賈 蔡，張 明 明，周 濤 發(fā)

（合肥工業(yè)大學(xué)資源與環(huán)境工程學(xué)院，安徽 合肥 230009）

指示克里格方法是一種不依賴于分布假設(shè)條件的非參數(shù)估值方法，對于異常值和偏態(tài)分布都具有良好的穩(wěn)健能力，因此可以作為地球化學(xué)異常研究的理想工具。閾值作為指示克里格方法最重要的參數(shù)，對于插值結(jié)果具有很大影響。該文從多維分形理論出發(fā)，利用濃度－面積（C-A）模型計算地球化學(xué)異常下限值，并將其作為閾值參與指示克里格法插值計算。為了解決多維分形測度時的不確定性，利用Voronoi圖的唯一性對C-A模型進(jìn)行了改進(jìn)，并通過安徽省某銅礦區(qū)大比例尺化探Cu元素地球化學(xué)數(shù)據(jù)對上述方法進(jìn)行了實例研究。結(jié)果顯示，較之反距離加權(quán)插值法、普通克里格插值法，指示克里格法獲取的最高累計概率范圍與已知礦體的空間吻合程度更高，具有更好的地球化學(xué)異常識別能力。對于地球化學(xué)數(shù)據(jù)空間變異性強(qiáng)烈的地區(qū)，指示克里格方法在穩(wěn)定變異函數(shù)形狀和高值信息重建能力方面具有很大的優(yōu)勢。

指示克里格;多維分形;濃度－面積模型;Voronoi圖;地球化學(xué)

0 引言

在針對勘查地球化學(xué)數(shù)據(jù)進(jìn)行異常識別的研究中，為了獲取地球化學(xué)場的空間變化趨勢以及異常范圍，通常需要對采樣點數(shù)據(jù)進(jìn)行空間插值分析和處理［1］。以地統(tǒng)計學(xué)為基礎(chǔ)的普通克里格插值法由于充分度量了空間變異信息以及具有無偏最優(yōu)的特點，被廣泛應(yīng)用于地質(zhì)、環(huán)境、生態(tài)等領(lǐng)域［2－7］。然而，地球化學(xué)數(shù)據(jù)往往包含有部分異常數(shù)據(jù)，常具有偏斜的非正態(tài)分布特征，雖然正態(tài)變換和異常值剔除等方法可以增強(qiáng)變異函數(shù)和普通克里格方法的穩(wěn)健性［8－10］，但同時會大幅增加克里格方法本身具有的平滑效應(yīng)［11］，這對于基于地球化學(xué)數(shù)據(jù)的異常識別極其不利。

指示克里格方法是一種不依賴于分布假設(shè)條件的非參數(shù)估值方法，對于具有異常值和偏態(tài)分布的數(shù)據(jù)都具有良好的穩(wěn)健能力［12］。閾值作為其最重要的參數(shù)，對于指示克里格插值的結(jié)果具有很大影響。通過設(shè)定表征異常下限的閾值，利用指示克里格方法可以獲取小于或大于該閾值的累計頻率［3］，從而從空間頻率的角度為地球化學(xué)數(shù)據(jù)的異常識別提供科學(xué)依據(jù)。

多維分形理論的提出為異常下限值的獲取提供了有利工具。由于受到多期次的地質(zhì)過程或人為活動的影響，地球化學(xué)場往往具有多維分形的特征。通過對地球化學(xué)場自相似性的研究，可以準(zhǔn)確把握地球化學(xué)場的多維分形特征，從而實現(xiàn)地球化學(xué)場背景與異常的分離［13－15］。近年來多種用于提取異常下限值的分形和多維分形模型相繼被提出，其中以濃度－面積模型（Concentration-Area Model）［16］的應(yīng)用最為廣泛，其他如周長－面積模型（Perimeter-Area Model）［17］、含 量 -距 離 模 型 （Concentration-Distance Model）［18］、含 量 -樣 品 數(shù) 量 模 型（Concentration-Number Model）［19］也 有 一 定 的 應(yīng)用。然而，目前基于異常下限值的應(yīng)用研究多局限于等值線圖 中 的 異 常 區(qū) 域 圈 定［15，20，21］，對 于 其 他 方面的應(yīng)用未進(jìn)行更深入的研究。

本文將上述濃度－面積多維分形模型與指示克里格方法結(jié)合，不僅能獲得更為合理的表征地球化學(xué)異常下限的閾值，還可以定量計算高于該閾值的概率分布信息。為了驗證上述方法的實際效果，以安徽省某銅礦區(qū)Cu元素地球化學(xué)數(shù)據(jù)進(jìn)行實例驗證，并與其他插值方法進(jìn)行對比。

1 指示克里格方法

指示克里格法屬于非線性克里格方法范疇，可用于估計某一位置超過指定閾值zk的累計頻率［3］。與普通克里格法相比，它不嚴(yán)格依賴于空間現(xiàn)象的平穩(wěn)性假設(shè)，也不要求區(qū)域化變量服從某種分布［22］;并且由于其在插值前根據(jù)一定的閾值將數(shù)據(jù)轉(zhuǎn)換為指示變量，所以對于異常值和偏態(tài)分布都具有良好的穩(wěn)健功能，指示變量的變換公式［3］為：

式中：zk為設(shè)定的閾值。如果樣品數(shù)據(jù)大于等于zk，則賦予指示變量 0;如果小于zk，則賦予指示變量1。

與普通克里格方法相似，對于任何一個待估點x0，z（x）≤zk的概率可以通過對鄰域內(nèi)指示變量進(jìn)行線性相加獲得［12］。普通指示克里格估值公式為：

式中：μ是拉格朗日算子，γi（xα－xβ;zk）是第α個和第β個樣品點的指示變量的變異函數(shù)值;γi（xα－x0;zk）是待估點x0與第α個樣品點的指示變量的變異函數(shù)值。

指示克里格方法采用的實驗變異函數(shù)計算方法同普通克里格方法基本相同，不同的是指示克里格方法采用指示變量計算實驗變異函數(shù)，公式［12］為：

式中：N代表滯后距離區(qū)間h內(nèi)的樣品對數(shù)。上述公式計算得到的實驗變異函數(shù)值γ（h）還需通過擬合求解理論變異函數(shù)的參數(shù)，從而參與指示克里格方程組的計算［23］。

2 多維分形模型

為了研究地球化學(xué)數(shù)據(jù)的空間域和頻率域特征，許多分形和多維分形模型被提出［16－19，24］，這些模型均假設(shè)地球化學(xué)數(shù)據(jù)呈現(xiàn)分形和多維分形分布并服從power-law函數(shù)關(guān)系。Power-law函數(shù)關(guān)系的最重要特性即尺度的獨立性［16，25］：

式中：測度M和尺度δ構(gòu)建了一個函數(shù)關(guān)系，E是拓?fù)渚S數(shù)，D是分形維數(shù)。D或分形余維（E-D）可以通過建立M與δ的雙對數(shù)圖由擬合直線的斜率估計得到。

測度M反映了分形的尺度獨立性，可用于計算分形維數(shù)。其中以面積作為測度的濃度－面積（CA）模型［16，26］應(yīng)用最為廣泛，公式如下：

式中：測度A（≥c）是指包圍在高于或等于c的等值線內(nèi)的面積，其也可以通過度量柵格數(shù)據(jù)得到。由于地球化學(xué)采樣點往往不遵循格網(wǎng)分布，且分布稀疏，局部區(qū)域分布的采樣點密度也不同，因此上述模型需要通過內(nèi)插處理獲取等值線和柵格數(shù)據(jù)。然而插值結(jié)果具有顯著的不確定性，反距離加權(quán)插值法、克里格插值法等方法產(chǎn)生的結(jié)果顯著不同，特別是對于具有異常數(shù)據(jù)的強(qiáng)變異性數(shù)據(jù)［14，27］。

為了消除上述插值方法帶來的影響，謝淑云等改進(jìn)C-A方法，提出了面積校正累計頻率法［27，28］，首先對研究區(qū)域進(jìn)行網(wǎng)格劃分，對各網(wǎng)格內(nèi)的樣品數(shù)據(jù)進(jìn)行綜合計算，然后在雙對數(shù)坐標(biāo)下建立C-N（C＞c）曲線獲取異常下限值。面積校正累計頻率法雖然消除了樣品點分布不均的影響，但在網(wǎng)格劃分中，對于平均網(wǎng)格密度d的設(shè)定仍具有主觀因素，過小的d常導(dǎo)致部分網(wǎng)格中沒有樣品，過大的d則導(dǎo)致數(shù)據(jù)產(chǎn)生“平滑”效應(yīng)。

鑒于上述問題，本文采用Voronoi圖對上述CA模型進(jìn)行改進(jìn)。Voronoi圖由一組連接兩鄰點直線的垂直平分線組成的連續(xù)多邊形組成，N個在平面上有區(qū)別的點按照最鄰近原則劃分平面，每個點與其最近鄰區(qū)域相關(guān)聯(lián)［29］。由于Voronoi圖具有唯一的性質(zhì)［30］，因此可以用于解決上述構(gòu)建多維分形測度時的不確定性，這對于采樣密度不均或空間分布形狀不規(guī)則的數(shù)據(jù)集在理論上更具優(yōu)勢。具體步驟為：1）根據(jù)樣品點數(shù)據(jù)的空間分布獲取研究區(qū)域邊界，并根據(jù)1/2倍基本采樣間距對區(qū)域邊界進(jìn)行擴(kuò)展;2）根據(jù)樣品點空間分布和研究區(qū)邊界建立Voronoi圖;3）計算各個Voronoi多邊形的面積并對應(yīng)于每個樣品點數(shù)值;4）選定一組在對數(shù)坐標(biāo)下的含量值c＝｛ci｝（i＝1，2，…，n），n代表分類數(shù)，統(tǒng)計所有平均含量值C＞c的Voronoi多邊形的累計面積A，并在雙對數(shù)坐標(biāo)下繪制C-A（C＞c）曲線;5）根據(jù)曲線的趨勢進(jìn)行分段線性擬合，擬合直線的交點即可作為需要獲取的異常下限值。

3 實例應(yīng)用及討論

3.1 實例數(shù)據(jù)和基本統(tǒng)計特征

本文實例數(shù)據(jù)為安徽省某銅礦區(qū)大比例尺土壤化探Cu元素數(shù)據(jù)，采樣點分布如圖1所示。土壤樣品按照網(wǎng)格采樣，采樣間距為200 m，共計241件;由于部分區(qū)域難以取樣，因此存在采樣點分布不規(guī)則或缺失現(xiàn)象。Cu元素數(shù)據(jù)的基本統(tǒng)計結(jié)果見表1，研究區(qū)域內(nèi)Cu元素的變異系數(shù)高達(dá)381.99%，具有非常強(qiáng)的空間變異性質(zhì)，偏度和峰度值則顯示出數(shù)據(jù)呈非正態(tài)分布，且明顯正偏，具有含量值很大的異常值。

圖1 研究區(qū)域及采樣點分布Fig.1 Study area and the distribution of sample points

表1 Cu元素含量基本統(tǒng)計量 /Table 1 The statistical results of Cu concentration mg kg

3.2 多維分形模型異常下限值提取

本文應(yīng)用Voronoi圖進(jìn)行面積的測度計算，并通過C-A模型求取Cu元素的異常下限值。其中Voronoi圖通過ArcGIS軟件計算得到，異常下限值計算和擬合通過Matlab軟件進(jìn)行分析處理。

首先定義研究區(qū)域的邊界（圖1），邊界的獲取對應(yīng)于采樣點的空間分布形狀，向外擴(kuò)展1/2倍采樣間距，以保證邊緣和中心樣品對結(jié)果作用的一致性?；跇悠泛脱芯繀^(qū)邊界計算得到Voronoi圖（圖2）。

圖2 采樣點Voronoi圖Fig.2 Voronoi diagram of sample points

計算Voronoi圖中每個多邊形的面積，在雙對數(shù)坐標(biāo)下建立Cu元素的C-A（濃度－Voronoi多邊形累計面積）多維分形模型（圖3）。從圖3可見，離散的數(shù)值點具有兩個明顯的似線性段，第一似線性段可以認(rèn)為是背景值條件下的低值波動［27］，而第二似線性段可以認(rèn)為是異常環(huán)境下的分形特征，也可以認(rèn)為是反映不同期次地質(zhì)或人為活動的地球化學(xué)場。通過對兩條似線性段采用最小二乘方法進(jìn)行擬合可以獲得兩條擬合直線，通過對兩條擬合直線交點的橫坐標(biāo)C求取反對數(shù)即可獲得待求的異常下限值，擬合直線方程及求取的異常下限值見表2。

圖3 C-A多重分形模型雙對數(shù)圖Fig.3 Log-Log plots of C-A multifractal model

表2 擬合直線方程及異常下限值Table 2 Equation of fitted line and the threshold

3.3 空間變異分析

變異函數(shù)分析是了解數(shù)據(jù)空間變異性質(zhì)的有效工具，也是克里格插值的必要條件。由于本文采用的數(shù)據(jù)分布形狀極不規(guī)則，且Voronoi圖顯示數(shù)據(jù)并無明顯的各向異性，因此僅基于各向同性條件對變異函數(shù)進(jìn)行計算。通過ArcGIS10軟件分別計算了普通克里格方法和指示克里格方法下的實驗變異函數(shù)，并根據(jù)實驗變異函數(shù)的形態(tài)特征，應(yīng)用帶有塊金效應(yīng)的球狀變異函數(shù)模型對其進(jìn)行擬合。其中，指示克里格方法下的變異函數(shù)采用C-A方法獲取的異常下限值作為閾值參數(shù)。由于本研究側(cè)重高值異常信息，因此對于大于所定義閾值的樣品數(shù)據(jù)賦予指示變量1，如果小于定義閾值則賦予指示變量0。計算得到的變異函數(shù)和擬合結(jié)果如圖4所示，理論變異函數(shù)參數(shù)見表3。

表3 變異函數(shù)參數(shù)Table 3 Parameters of the variograms

圖4 實驗變異函數(shù)與擬合變異函數(shù)模型Fig.4 Experimental variograms with fitted variogram models：Ordinary Kriging and Indicator Kriging

從圖4a可見，由于普通克里格方法采用的變異函數(shù)由原始數(shù)據(jù)直接計算得到，同時原始數(shù)據(jù)具有很強(qiáng)的變異性質(zhì)且含有很高的異常數(shù)值，因此實驗變異函數(shù)值的趨勢混亂且躍動明顯，難以擬合出符合實驗變異函數(shù)趨勢的理論變異函數(shù)模型。塊金值和基臺值的比值大于0.8則說明區(qū)域內(nèi)隨機(jī)波動占據(jù)了空間變異性質(zhì)的絕大部分，因此很難反映數(shù)據(jù)的空間變異特征。圖4b中基于指示變量的實驗變異函數(shù)則顯示出良好的空間變異結(jié)構(gòu)，變化趨勢穩(wěn)定且符合球狀理論變異函數(shù)結(jié)構(gòu)，較之普通克里格方法能夠更好地描述空間數(shù)據(jù)的相關(guān)性和結(jié)構(gòu)性。

3.4 指示克里格插值

為了對插值結(jié)果進(jìn)行對比，本文同時對數(shù)據(jù)進(jìn)行反距離加權(quán)插值及采用平均值作為閾值的指示克里格插值計算。插值的過程統(tǒng)一設(shè)定最大搜索點數(shù)為12，指示克里格插值的概率結(jié)果采用等距（Equal Intervals）方法進(jìn)行分級，而克里格和反距離加權(quán)插值結(jié)果則采用幾何間距（Geometric Intervals）方法進(jìn)行分級，以保證其獲得良好的對比程度，插值結(jié)果見圖5。

圖5 插值結(jié)果等值線Fig.5 Contour maps of interpolation results

由插值結(jié)果可見，反距離加權(quán)插值和普通克里格插值方法均無法以最高值等值線區(qū)域識別出已知礦體（圖5a、圖5b）。普通克里格插值方法由于更注重表征數(shù)據(jù)的空間分布趨勢，因此較之反距離加權(quán)插值法對于數(shù)據(jù)具有更顯著的平滑效應(yīng)，具體表現(xiàn)為在圖5b中高值區(qū)域的缺失以及中高值區(qū)域在大范圍內(nèi)的均一分布。而基于Voronoi圖的C-A方法獲取閾值的指示克里格方法則具有很好的地球化學(xué)異常識別能力和高值信息重建能力，插值結(jié)果（圖5c）中最高值等值線區(qū)域（0.9～1.0）不但能夠很好的與已知礦體的空間位置吻合，而且較之反距離加權(quán)插值法、普通克里格插值方法能夠更好的突出反映原始數(shù)據(jù)的高值異常信息。圖5d是采用數(shù)據(jù)的平均值作為閾值的指示克里格插值結(jié)果，可見，由于研究區(qū)域內(nèi)存在大量高值異常數(shù)據(jù)，因此高累計頻率區(qū)域僅位于研究區(qū)的南部，而礦體位置則表現(xiàn)為中低的累計頻率分布。傳統(tǒng)的異常下限值通常利用平均值與1倍或者2倍標(biāo)準(zhǔn)差之和的方法求取［21，31］，其遠(yuǎn)大于平均值，因而更難以將其作為指示克里格方法的閾值對礦體位置進(jìn)行有效識別。因此，本文采用的基于Voronoi圖的C-A方法，可以更為有效地求取指示克里格方法的閾值參數(shù)，并參與指示克里格插值計算，更好地服務(wù)于針對具有多維分形性質(zhì)的非平穩(wěn)數(shù)據(jù)的異常識別研究。

4 結(jié)論

Voronoi圖具有唯一性，可用于解決多維分形模型計算測度時的不確定性?；赩oronoi圖的C－A方法可以更好的應(yīng)用于地球化學(xué)數(shù)據(jù)的多維分形特征分析及異常下限的提取，可作為指示克里格方法閾值參數(shù)求取的有效方法。實例研究顯示，對于地球化學(xué)數(shù)據(jù)空間變異性強(qiáng)烈的地區(qū)，較之反距離加權(quán)插值法和普通克里格插值方法，利用上述方法獲取閾值的指示克里格方法具有更好的地球化學(xué)異常識別能力和高值信息重建能力，所得結(jié)果的最高累計頻率值范圍與已知礦體具有更好的空間吻合程度。對于高偏斜、強(qiáng)變異且具有異常值的數(shù)據(jù)集，較之平均值作為指示克里格方法的閾值參數(shù)，基于上述C-A方法獲取閾值的指示克里格方法可以更為有效地進(jìn)行地球化學(xué)異常識別研究。

［1］成秋明.多重分形與地質(zhì)統(tǒng)計學(xué)方法用于勘查地球化學(xué)異?？臻g結(jié)構(gòu)和奇異性分析［J］.中國地質(zhì)大學(xué)學(xué)報（地球科學(xué)），2001，26（2）：161－164.

［2］KRIGE D.A statistical approach to some basic mine evaluation problems on the Witwateround［J］.Chim.Min.Soc.South-Africa，1951，52：119－139.

［3］張仁鐸.空間變異理論及應(yīng)用［M］.北京：科學(xué)出版社，2005.

［4］JOURNEL A，HUIJBREGTS C.Mining Geostatistics［M］.New York：Academic San Diego，1978.1－600.

［5］孫洪泉.地質(zhì)統(tǒng)計學(xué)及其應(yīng)用［M］.徐州：中國礦業(yè)大學(xué)出版杜，1990.1－252.

［6］王政權(quán).地統(tǒng)計學(xué)及在生態(tài)學(xué)中的應(yīng)用［M］.北京：科學(xué)出版社，1999.1－195.

［7］MCGRATH D，ZHANG C，CARTON O T.Geostatistical analyses and hazard assessment on soil lead in Silvermines area，Ireland［J］.Environmental Pollution，2004，127（2）：239－248.

［8］李蒙文，戰(zhàn)明國，趙財勝，等.穩(wěn)健估計方法在內(nèi)蒙古新忽熱地區(qū)水系沉積物測量異常評價中的應(yīng)用［J］.礦床地質(zhì)，2006，25（1）：27－35.

［9］張長波，吳龍華，駱永明.穩(wěn)健變異函數(shù)在土壤污染物來源識別中的應(yīng)用：以某重金屬污染場地為例［J］.環(huán)境科學(xué)，2008，29（3）：804－808.

［10］李曉暉，袁峰，白曉宇，等.典型礦區(qū)非正態(tài)分布土壤元素數(shù)據(jù)的正態(tài)變換方法對比研究［J］.地理與地理信息科學(xué)，2010，26（6）：102－105.

［11］李慶謀.多維分形克里格方法［J］.地球科學(xué)進(jìn)展，2005，20（2）：248－255.

［12］LIN Y P，CHANG T，SHIH C，et al.Factorial and indicator Kriging methods using a geographic information system to delineate spatial variation and pollution sources of soil heavy metals［J］.Environmental Geology，2002，42（4）：900－909.

［13］成秋明，張生元，左仁廣，等.多重分形濾波方法和地球化學(xué)信息提取技術(shù)研究與進(jìn)展［J］.地學(xué)前緣，2009，16（2）：185－198.

［14］文戰(zhàn)久，高星，姚振興.基于“元素含量－面積”模型方法的地球化學(xué)場的多重分形模式分析［J］.地球科學(xué)進(jìn)展，2007，22（6）：598－604.

［15］白曉宇，袁峰，周濤發(fā)，等.多重分形方法識別銅陵礦區(qū)土壤中Cd的地球化學(xué)異常［J］.礦物巖石地球化學(xué)通報，2008，27（3）：306－310.

［16］CHENG Q.Spatial and scaling modelling for geochemical anomaly separation［J］.Geochemical Exploration，1999，65（3）：175－194.

［17］CHENG Q.The perimeter-area fractal model and its application to geology［J］.Mathematical Geology，1995，27（1）：69－82.

［18］LI C，MA T，SHI J.Application of a fractal method relating concentrations and distances for separation of geochemical anomalies from background［J］.Geochemical Exploration，2003，77（2－3）：167－175.

［19］MAO Z，PENG S，LAI J，et al.Fractal study of geochemical prospecting data in south area of Fenghuanshan copper deposit，Tongling Anhui［J］.Earth Sciences and Environment，2004，26（4）：11－14.

［20］孫忠軍.礦產(chǎn)勘查中化探異常下限的多重分形計算方法［J］.物探化探計算技術(shù)，2007，29（1）：54－57.

［21］BAI J，PORWAL A，HART C，et al.Mapping geochemical singularity using multifractal analysis：Application to anomaly definition on stream sediments data from Funin Sheet，Yunnan，China［J］.Geochemical Exploration，2010，104（1－2）：1－11.

［22］李章林，張夏林，翁正平.指示克里格法在礦體儲量計算方面的研究與應(yīng)用［J］.礦業(yè)快報，2008，24（1）：11－15.

［23］WEBSTER R，OLIVER M A.Geostatistics for Environmental Scientists［M］.Chichester：Wiley，2001.

［24］CHENG Q，XU Y，GRUNSKY E.Multifractal power spectrum-area method for geochemical anomaly separation［J］.Natural Resources Research，2000，9（1）：43－51.

［25］成秋明.多維分形理論和地球化學(xué)元素分布規(guī)律［J］.中國地質(zhì)大學(xué)學(xué)報（地球科學(xué)），2000，25（3）：311－317.

［26］CHENG Q，AGTERBERG F，BALLANTYNE S.The separation of geochemical anomalies from background by fractal methods［J］.Geochemical Exploration，1994，51（2）：109－130.

［27］謝淑云，鮑征宇.地球化學(xué)場的連續(xù)多重分形模式［J］.地球化學(xué)，2002，31（2）：191－200.

［28］XIE S Y，BAO Z Y.Fractal and multifractal properties of geochemical fields［J］.Mathem Atical Geology，2004，36（7）：847－864.

［29］吳立新，郝海森，殷作如.基于鉆孔點集Voronoi圖的礦產(chǎn)儲量新算法［J］.地理與地理信息科學(xué)，2004，20（1）：54－59.

［30］周小平，周瑞忠.基于Voronoi圖的新型幾何插值及其與傳統(tǒng)代數(shù)插值方法的比較［J］.巖石力學(xué)與工程學(xué)報，2005，24（1）：133－138.

［31］HARRIS J，WILKINSON L，GRUNSKY E，et al.Techniques for analysis and visualization of lithogeochemical data with applications to the Swayze Greenstone Belt，Ontario［J］.Geochemical Exploration，1999，67（1－3）：301－334.

Study on Anomaly Recognition from Geochemical Data Based on Multifractal Model and Indicator Kriging Method

LI Xiao－h(huán)ui，YUAN Feng，JIA Cai，ZHANG Ming－ming，ZHOU Tao－fa
（SchoolofResourcesandEnvironmentalEngineering，HefeiUniversityofTechnology，Hefei230009，China）

Indicator Kriging is a nonparametric estimate method which does not rely on the assumption of distribution.It has strongly robust capability for outliers and skewed distribution，so it can be used for the anomaly recognition of geochemical data very well.The threshold is the most important parameter of indicator Kriging，which can influence the results heavily.Based on multifractal theory，this paper calculates the threshold of geochemical data by using the Concentration-Area（C-A）model，and then uses the threshold for the indicator Kriging interploation method.This paper uses the uniqueness of Voronoi diagram to slove the uncertainty in calculating the multifractal measure of the C-A model，and gives a case study by using a macrocale Cu concentration geochemical data of a copper mine in Anhui Province.The results show that，compared with inverse distance weighted method and ordinary Kriging method，the highest rank of cumulative probability which obtained by indicator Kriging method can more effectively highlight the geochemical anomalies which associated with the known mineralization.Furthermore，for the strongly spatially variable data，indicator Kriging method has more advantages in reconstruction of highly information and maintaining the robustness of variogram.

indicator kriging;multifractal;concentration-area model;Voronoi diagram;geochemistry

P595

A

1672－0504（2011）06－0023－05

2011－06- 05;

2011－08－01

新世紀(jì)優(yōu)秀人才支持計劃項目（NCET－10－0324）;安徽省科技攻關(guān)計劃項目（08010302200）;安徽省公益性地質(zhì)（科技）工作項目（2009－13）;安徽省優(yōu)秀青年科技基金項目（08040106907、04045063）

李曉暉（1986－），男，博士研究生，主要從事多維分形及地質(zhì)體三維建模預(yù)測研究。E－mail：lxhlixiaohui＠163.com