基于多维分形模型与指示克里格方法的地球化学异常识别研究

李 晓 晖，袁 峰，贾 蔡，张 明 明，周 涛 发

基于多维分形模型与指示克里格方法的地球化学异常识别研究

李 晓 晖，袁 峰，贾 蔡，张 明 明，周 涛 发

（合肥工业大学资源与环境工程学院，安徽 合肥 230009）

指示克里格方法是一种不依赖于分布假设条件的非参数估值方法，对于异常值和偏态分布都具有良好的稳健能力，因此可以作为地球化学异常研究的理想工具。阈值作为指示克里格方法最重要的参数，对于插值结果具有很大影响。该文从多维分形理论出发，利用浓度－面积（C-A）模型计算地球化学异常下限值，并将其作为阈值参与指示克里格法插值计算。为了解决多维分形测度时的不确定性，利用Voronoi图的唯一性对C-A模型进行了改进，并通过安徽省某铜矿区大比例尺化探Cu元素地球化学数据对上述方法进行了实例研究。结果显示，较之反距离加权插值法、普通克里格插值法，指示克里格法获取的最高累计概率范围与已知矿体的空间吻合程度更高，具有更好的地球化学异常识别能力。对于地球化学数据空间变异性强烈的地区，指示克里格方法在稳定变异函数形状和高值信息重建能力方面具有很大的优势。

指示克里格;多维分形;浓度－面积模型;Voronoi图;地球化学

0 引言

在针对勘查地球化学数据进行异常识别的研究中，为了获取地球化学场的空间变化趋势以及异常范围，通常需要对采样点数据进行空间插值分析和处理［1］。以地统计学为基础的普通克里格插值法由于充分度量了空间变异信息以及具有无偏最优的特点，被广泛应用于地质、环境、生态等领域［2－7］。然而，地球化学数据往往包含有部分异常数据，常具有偏斜的非正态分布特征，虽然正态变换和异常值剔除等方法可以增强变异函数和普通克里格方法的稳健性［8－10］，但同时会大幅增加克里格方法本身具有的平滑效应［11］，这对于基于地球化学数据的异常识别极其不利。

指示克里格方法是一种不依赖于分布假设条件的非参数估值方法，对于具有异常值和偏态分布的数据都具有良好的稳健能力［12］。阈值作为其最重要的参数，对于指示克里格插值的结果具有很大影响。通过设定表征异常下限的阈值，利用指示克里格方法可以获取小于或大于该阈值的累计频率［3］，从而从空间频率的角度为地球化学数据的异常识别提供科学依据。

多维分形理论的提出为异常下限值的获取提供了有利工具。由于受到多期次的地质过程或人为活动的影响，地球化学场往往具有多维分形的特征。通过对地球化学场自相似性的研究，可以准确把握地球化学场的多维分形特征，从而实现地球化学场背景与异常的分离［13－15］。近年来多种用于提取异常下限值的分形和多维分形模型相继被提出，其中以浓度－面积模型（Concentration-Area Model）［16］的应用最为广泛，其他如周长－面积模型（Perimeter-Area Model）［17］、含 量 -距 离 模 型 （Concentration-Distance Model）［18］、含 量 -样 品 数 量 模 型（Concentration-Number Model）［19］也 有 一 定 的 应用。然而，目前基于异常下限值的应用研究多局限于等值线图 中 的 异 常 区 域 圈 定［15，20，21］，对 于 其 他 方面的应用未进行更深入的研究。

本文将上述浓度－面积多维分形模型与指示克里格方法结合，不仅能获得更为合理的表征地球化学异常下限的阈值，还可以定量计算高于该阈值的概率分布信息。为了验证上述方法的实际效果，以安徽省某铜矿区Cu元素地球化学数据进行实例验证，并与其他插值方法进行对比。

1 指示克里格方法

指示克里格法属于非线性克里格方法范畴，可用于估计某一位置超过指定阈值zk的累计频率［3］。与普通克里格法相比，它不严格依赖于空间现象的平稳性假设，也不要求区域化变量服从某种分布［22］;并且由于其在插值前根据一定的阈值将数据转换为指示变量，所以对于异常值和偏态分布都具有良好的稳健功能，指示变量的变换公式［3］为：

式中：zk为设定的阈值。如果样品数据大于等于zk，则赋予指示变量 0;如果小于zk，则赋予指示变量1。

与普通克里格方法相似，对于任何一个待估点x0，z（x）≤zk的概率可以通过对邻域内指示变量进行线性相加获得［12］。普通指示克里格估值公式为：

式中：μ是拉格朗日算子，γi（xα－xβ;zk）是第α个和第β个样品点的指示变量的变异函数值;γi（xα－x0;zk）是待估点x0与第α个样品点的指示变量的变异函数值。

指示克里格方法采用的实验变异函数计算方法同普通克里格方法基本相同，不同的是指示克里格方法采用指示变量计算实验变异函数，公式［12］为：

式中：N代表滞后距离区间h内的样品对数。上述公式计算得到的实验变异函数值γ（h）还需通过拟合求解理论变异函数的参数，从而参与指示克里格方程组的计算［23］。

2 多维分形模型

为了研究地球化学数据的空间域和频率域特征，许多分形和多维分形模型被提出［16－19，24］，这些模型均假设地球化学数据呈现分形和多维分形分布并服从power-law函数关系。Power-law函数关系的最重要特性即尺度的独立性［16，25］：

式中：测度M和尺度δ构建了一个函数关系，E是拓扑维数，D是分形维数。D或分形余维（E-D）可以通过建立M与δ的双对数图由拟合直线的斜率估计得到。

测度M反映了分形的尺度独立性，可用于计算分形维数。其中以面积作为测度的浓度－面积（CA）模型［16，26］应用最为广泛，公式如下：

式中：测度A（≥c）是指包围在高于或等于c的等值线内的面积，其也可以通过度量栅格数据得到。由于地球化学采样点往往不遵循格网分布，且分布稀疏，局部区域分布的采样点密度也不同，因此上述模型需要通过内插处理获取等值线和栅格数据。然而插值结果具有显著的不确定性，反距离加权插值法、克里格插值法等方法产生的结果显著不同，特别是对于具有异常数据的强变异性数据［14，27］。

为了消除上述插值方法带来的影响，谢淑云等改进C-A方法，提出了面积校正累计频率法［27，28］，首先对研究区域进行网格划分，对各网格内的样品数据进行综合计算，然后在双对数坐标下建立C-N（C＞c）曲线获取异常下限值。面积校正累计频率法虽然消除了样品点分布不均的影响，但在网格划分中，对于平均网格密度d的设定仍具有主观因素，过小的d常导致部分网格中没有样品，过大的d则导致数据产生“平滑”效应。

鉴于上述问题，本文采用Voronoi图对上述CA模型进行改进。Voronoi图由一组连接两邻点直线的垂直平分线组成的连续多边形组成，N个在平面上有区别的点按照最邻近原则划分平面，每个点与其最近邻区域相关联［29］。由于Voronoi图具有唯一的性质［30］，因此可以用于解决上述构建多维分形测度时的不确定性，这对于采样密度不均或空间分布形状不规则的数据集在理论上更具优势。具体步骤为：1）根据样品点数据的空间分布获取研究区域边界，并根据1/2倍基本采样间距对区域边界进行扩展;2）根据样品点空间分布和研究区边界建立Voronoi图;3）计算各个Voronoi多边形的面积并对应于每个样品点数值;4）选定一组在对数坐标下的含量值c＝｛ci｝（i＝1，2，…，n），n代表分类数，统计所有平均含量值C＞c的Voronoi多边形的累计面积A，并在双对数坐标下绘制C-A（C＞c）曲线;5）根据曲线的趋势进行分段线性拟合，拟合直线的交点即可作为需要获取的异常下限值。

3 实例应用及讨论

3.1 实例数据和基本统计特征

本文实例数据为安徽省某铜矿区大比例尺土壤化探Cu元素数据，采样点分布如图1所示。土壤样品按照网格采样，采样间距为200 m，共计241件;由于部分区域难以取样，因此存在采样点分布不规则或缺失现象。Cu元素数据的基本统计结果见表1，研究区域内Cu元素的变异系数高达381.99%，具有非常强的空间变异性质，偏度和峰度值则显示出数据呈非正态分布，且明显正偏，具有含量值很大的异常值。

图1 研究区域及采样点分布Fig.1 Study area and the distribution of sample points

表1 Cu元素含量基本统计量 /Table 1 The statistical results of Cu concentration mg kg

3.2 多维分形模型异常下限值提取

本文应用Voronoi图进行面积的测度计算，并通过C-A模型求取Cu元素的异常下限值。其中Voronoi图通过ArcGIS软件计算得到，异常下限值计算和拟合通过Matlab软件进行分析处理。

首先定义研究区域的边界（图1），边界的获取对应于采样点的空间分布形状，向外扩展1/2倍采样间距，以保证边缘和中心样品对结果作用的一致性。基于样品和研究区边界计算得到Voronoi图（图2）。

图2 采样点Voronoi图Fig.2 Voronoi diagram of sample points

计算Voronoi图中每个多边形的面积，在双对数坐标下建立Cu元素的C-A（浓度－Voronoi多边形累计面积）多维分形模型（图3）。从图3可见，离散的数值点具有两个明显的似线性段，第一似线性段可以认为是背景值条件下的低值波动［27］，而第二似线性段可以认为是异常环境下的分形特征，也可以认为是反映不同期次地质或人为活动的地球化学场。通过对两条似线性段采用最小二乘方法进行拟合可以获得两条拟合直线，通过对两条拟合直线交点的横坐标C求取反对数即可获得待求的异常下限值，拟合直线方程及求取的异常下限值见表2。

图3 C-A多重分形模型双对数图Fig.3 Log-Log plots of C-A multifractal model

表2 拟合直线方程及异常下限值Table 2 Equation of fitted line and the threshold

3.3 空间变异分析

变异函数分析是了解数据空间变异性质的有效工具，也是克里格插值的必要条件。由于本文采用的数据分布形状极不规则，且Voronoi图显示数据并无明显的各向异性，因此仅基于各向同性条件对变异函数进行计算。通过ArcGIS10软件分别计算了普通克里格方法和指示克里格方法下的实验变异函数，并根据实验变异函数的形态特征，应用带有块金效应的球状变异函数模型对其进行拟合。其中，指示克里格方法下的变异函数采用C-A方法获取的异常下限值作为阈值参数。由于本研究侧重高值异常信息，因此对于大于所定义阈值的样品数据赋予指示变量1，如果小于定义阈值则赋予指示变量0。计算得到的变异函数和拟合结果如图4所示，理论变异函数参数见表3。

表3 变异函数参数Table 3 Parameters of the variograms

图4 实验变异函数与拟合变异函数模型Fig.4 Experimental variograms with fitted variogram models：Ordinary Kriging and Indicator Kriging

从图4a可见，由于普通克里格方法采用的变异函数由原始数据直接计算得到，同时原始数据具有很强的变异性质且含有很高的异常数值，因此实验变异函数值的趋势混乱且跃动明显，难以拟合出符合实验变异函数趋势的理论变异函数模型。块金值和基台值的比值大于0.8则说明区域内随机波动占据了空间变异性质的绝大部分，因此很难反映数据的空间变异特征。图4b中基于指示变量的实验变异函数则显示出良好的空间变异结构，变化趋势稳定且符合球状理论变异函数结构，较之普通克里格方法能够更好地描述空间数据的相关性和结构性。

3.4 指示克里格插值

为了对插值结果进行对比，本文同时对数据进行反距离加权插值及采用平均值作为阈值的指示克里格插值计算。插值的过程统一设定最大搜索点数为12，指示克里格插值的概率结果采用等距（Equal Intervals）方法进行分级，而克里格和反距离加权插值结果则采用几何间距（Geometric Intervals）方法进行分级，以保证其获得良好的对比程度，插值结果见图5。

图5 插值结果等值线Fig.5 Contour maps of interpolation results

由插值结果可见，反距离加权插值和普通克里格插值方法均无法以最高值等值线区域识别出已知矿体（图5a、图5b）。普通克里格插值方法由于更注重表征数据的空间分布趋势，因此较之反距离加权插值法对于数据具有更显著的平滑效应，具体表现为在图5b中高值区域的缺失以及中高值区域在大范围内的均一分布。而基于Voronoi图的C-A方法获取阈值的指示克里格方法则具有很好的地球化学异常识别能力和高值信息重建能力，插值结果（图5c）中最高值等值线区域（0.9～1.0）不但能够很好的与已知矿体的空间位置吻合，而且较之反距离加权插值法、普通克里格插值方法能够更好的突出反映原始数据的高值异常信息。图5d是采用数据的平均值作为阈值的指示克里格插值结果，可见，由于研究区域内存在大量高值异常数据，因此高累计频率区域仅位于研究区的南部，而矿体位置则表现为中低的累计频率分布。传统的异常下限值通常利用平均值与1倍或者2倍标准差之和的方法求取［21，31］，其远大于平均值，因而更难以将其作为指示克里格方法的阈值对矿体位置进行有效识别。因此，本文采用的基于Voronoi图的C-A方法，可以更为有效地求取指示克里格方法的阈值参数，并参与指示克里格插值计算，更好地服务于针对具有多维分形性质的非平稳数据的异常识别研究。

4 结论

Voronoi图具有唯一性，可用于解决多维分形模型计算测度时的不确定性。基于Voronoi图的C－A方法可以更好的应用于地球化学数据的多维分形特征分析及异常下限的提取，可作为指示克里格方法阈值参数求取的有效方法。实例研究显示，对于地球化学数据空间变异性强烈的地区，较之反距离加权插值法和普通克里格插值方法，利用上述方法获取阈值的指示克里格方法具有更好的地球化学异常识别能力和高值信息重建能力，所得结果的最高累计频率值范围与已知矿体具有更好的空间吻合程度。对于高偏斜、强变异且具有异常值的数据集，较之平均值作为指示克里格方法的阈值参数，基于上述C-A方法获取阈值的指示克里格方法可以更为有效地进行地球化学异常识别研究。

［1］成秋明.多重分形与地质统计学方法用于勘查地球化学异常空间结构和奇异性分析［J］.中国地质大学学报（地球科学），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］张仁铎.空间变异理论及应用［M］.北京：科学出版社，2005.

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

［5］孙洪泉.地质统计学及其应用［M］.徐州：中国矿业大学出版杜，1990.1－252.

［6］王政权.地统计学及在生态学中的应用［M］.北京：科学出版社，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］李蒙文，战明国，赵财胜，等.稳健估计方法在内蒙古新忽热地区水系沉积物测量异常评价中的应用［J］.矿床地质，2006，25（1）：27－35.

［9］张长波，吴龙华，骆永明.稳健变异函数在土壤污染物来源识别中的应用：以某重金属污染场地为例［J］.环境科学，2008，29（3）：804－808.

［10］李晓晖，袁峰，白晓宇，等.典型矿区非正态分布土壤元素数据的正态变换方法对比研究［J］.地理与地理信息科学，2010，26（6）：102－105.

［11］李庆谋.多维分形克里格方法［J］.地球科学进展，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］成秋明，张生元，左仁广，等.多重分形滤波方法和地球化学信息提取技术研究与进展［J］.地学前缘，2009，16（2）：185－198.

［14］文战久，高星，姚振兴.基于“元素含量－面积”模型方法的地球化学场的多重分形模式分析［J］.地球科学进展，2007，22（6）：598－604.

［15］白晓宇，袁峰，周涛发，等.多重分形方法识别铜陵矿区土壤中Cd的地球化学异常［J］.矿物岩石地球化学通报，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］孙忠军.矿产勘查中化探异常下限的多重分形计算方法［J］.物探化探计算技术，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］李章林，张夏林，翁正平.指示克里格法在矿体储量计算方面的研究与应用［J］.矿业快报，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］成秋明.多维分形理论和地球化学元素分布规律［J］.中国地质大学学报（地球科学），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］谢淑云，鲍征宇.地球化学场的连续多重分形模式［J］.地球化学，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图的矿产储量新算法［J］.地理与地理信息科学，2004，20（1）：54－59.

［30］周小平，周瑞忠.基于Voronoi图的新型几何插值及其与传统代数插值方法的比较［J］.岩石力学与工程学报，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－hui，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

新世纪优秀人才支持计划项目（NCET－10－0324）;安徽省科技攻关计划项目（08010302200）;安徽省公益性地质（科技）工作项目（2009－13）;安徽省优秀青年科技基金项目（08040106907、04045063）

李晓晖（1986－），男，博士研究生，主要从事多维分形及地质体三维建模预测研究。E－mail：lxhlixiaohui＠163.com