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

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

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

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

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

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

0 引言

在針對(duì)勘查地球化學(xué)數(shù)據(jù)進(jìn)行異常識(shí)別的研究中，為了獲取地球化學(xué)場(chǎng)的空間變化趨勢(shì)以及異常范圍，通常需要對(duì)采樣點(diǎn)數(shù)據(jù)進(jìn)行空間插值分析和處理［1］。以地統(tǒng)計(jì)學(xué)為基礎(chǔ)的普通克里格插值法由于充分度量了空間變異信息以及具有無(wú)偏最優(yōu)的特點(diǎn)，被廣泛應(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］，但同時(shí)會(huì)大幅增加克里格方法本身具有的平滑效應(yīng)［11］，這對(duì)于基于地球化學(xué)數(shù)據(jù)的異常識(shí)別極其不利。

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

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

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

1 指示克里格方法

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

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

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

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

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

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

2 多維分形模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 采樣點(diǎn)Voronoi圖Fig.2 Voronoi diagram of sample points

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

圖3 C-A多重分形模型雙對(duì)數(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ù)并無(wú)明顯的各向異性，因此僅基于各向同性條件對(duì)變異函數(shù)進(jìn)行計(jì)算。通過(guò)ArcGIS10軟件分別計(jì)算了普通克里格方法和指示克里格方法下的實(shí)驗(yàn)變異函數(shù)，并根據(jù)實(shí)驗(yàn)變異函數(shù)的形態(tài)特征，應(yīng)用帶有塊金效應(yīng)的球狀變異函數(shù)模型對(duì)其進(jìn)行擬合。其中，指示克里格方法下的變異函數(shù)采用C-A方法獲取的異常下限值作為閾值參數(shù)。由于本研究側(cè)重高值異常信息，因此對(duì)于大于所定義閾值的樣品數(shù)據(jù)賦予指示變量1，如果小于定義閾值則賦予指示變量0。計(jì)算得到的變異函數(shù)和擬合結(jié)果如圖4所示，理論變異函數(shù)參數(shù)見表3。

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

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

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

3.4 指示克里格插值

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

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

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

4 結(jié)論

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

［1］成秋明.多重分形與地質(zhì)統(tǒng)計(jì)學(xué)方法用于勘查地球化學(xué)異?？臻g結(jié)構(gòu)和奇異性分析［J］.中國(guó)地質(zhì)大學(xué)學(xué)報(bào)（地球科學(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)計(jì)學(xué)及其應(yīng)用［M］.徐州：中國(guó)礦業(yè)大學(xué)出版杜，1990.1－252.

［6］王政權(quán).地統(tǒng)計(jì)學(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)明國(guó)，趙財(cái)勝，等.穩(wěn)健估計(jì)方法在內(nèi)蒙古新忽熱地區(qū)水系沉積物測(cè)量異常評(píng)價(jià)中的應(yīng)用［J］.礦床地質(zhì)，2006，25（1）：27－35.

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

［10］李曉暉，袁峰，白曉宇，等.典型礦區(qū)非正態(tài)分布土壤元素?cái)?shù)據(jù)的正態(tài)變換方法對(duì)比研究［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é)場(chǎng)的多重分形模式分析［J］.地球科學(xué)進(jìn)展，2007，22（6）：598－604.

［15］白曉宇，袁峰，周濤發(fā)，等.多重分形方法識(shí)別銅陵礦區(qū)土壤中Cd的地球化學(xué)異常［J］.礦物巖石地球化學(xué)通報(bào)，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ì)算方法［J］.物探化探計(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］李章林，張夏林，翁正平.指示克里格法在礦體儲(chǔ)量計(jì)算方面的研究與應(yīng)用［J］.礦業(yè)快報(bào)，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］.中國(guó)地質(zhì)大學(xué)學(xué)報(bào)（地球科學(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é)場(chǎng)的連續(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］吳立新，郝海森，殷作如.基于鉆孔點(diǎn)集Voronoi圖的礦產(chǎn)儲(chǔ)量新算法［J］.地理與地理信息科學(xué)，2004，20（1）：54－59.

［30］周小平，周瑞忠.基于Voronoi圖的新型幾何插值及其與傳統(tǒng)代數(shù)插值方法的比較［J］.巖石力學(xué)與工程學(xué)報(bào)，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)秀人才支持計(jì)劃項(xiàng)目（NCET－10－0324）;安徽省科技攻關(guān)計(jì)劃項(xiàng)目（08010302200）;安徽省公益性地質(zhì)（科技）工作項(xiàng)目（2009－13）;安徽省優(yōu)秀青年科技基金項(xiàng)目（08040106907、04045063）

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