李晓晖
, 袁峰, 贾蔡, 张明明, 周涛发
合肥工业大学资源与环境工程学院, 安徽 合肥 230009
LI Xiao-hui, YUAN Feng, JIA Cai, ZHANG Ming-ming, ZHOU Tao-fa
中图分类号: P628
文献标识码: A
文章编号: 1000-0690(2012)02-0136-07
通讯作者:
收稿日期: 2011-03-10
修回日期: 2011-06-11
网络出版日期: 2012-02-20
版权声明: 2012 《地理科学》编辑部 本文是开放获取期刊文献，在以下情况下可以自由使用：学术研究、学术交流、科研教学等，但不允许用于商业目的.
基金资助:
作者简介:
作者简介：李晓晖（1986-）,男,河北定州市人,博士研究生,主要从事多维分形及地质体三维建模及预测研究。E-mail:lxhlixiaohui@163.com
展开
摘要
以铜陵矿集区土壤Pb元素为例,研究稀疏采样条件下地统计学克里格方法,序贯高斯模拟方法对奇异性指数计算的影响。研究结果表明,序贯高斯模拟方法强调了短距离范围内的空间不确定性,弥补了克里格方法平滑效应的不足,对于精细重建土壤元素的空间分布特征具有更好的效果。对于稀疏采样的数据集,较之原始数据和克里格方法,基于序贯高斯模拟方法获取的奇异性指数能够更精细的刻画局部空间结构,更好的应用于土壤地球化学异常的提取和识别。
关键词：
Abstract
Using data of Pb concentrations within soil samples from the Tongling Mining District area as an example, Kriging and Sequential Gaussian Simulation were used to determine the local singularity exponent within a dataset with low spatial density sampling. Within this highly variable dataset, when Sequential Gaussian Simulation was used ,the calculated maximum concentration of Pb and the standard deviation and coefficient of variation were similar to those of the raw dataset, whereas when the method of Kriging was applied the same statistics were significantly lower than those within the raw dataset. This means that the Sequential Gaussian Simulation method can be used to interpolate soil geochemical data without significant smoothing of the dataset, enabling the highlighting of spatial variations over even short distances and potentially being a better method of interpolation prior to interpretation than Kriging, which may leading to a loss of resolution. Furthermore, the method of Kriging can also compress the range of the singularity exponent calculated by this method, whereas when the Sequential Gaussian Simulation method was applied to the data, the range and contrast between individual data points were significantly increased, improving the user’s ability to identify anomalies within the dataset. Semi-variograms constructed from data interpolated using the Sequential Gaussian Simulation method are more homogenous and correspond well with the semi-variograms constructed from the raw data, especially between shorter distances, when compared to those constructed from Kriged data. For datasets with low spatial density, for example those generated by regional governmental or reconnaissance soil sampling, when compared with the raw and Kriged data, the Singularity Exponent calculated using Sequential Gaussian Simulation is significantly better at locating spatial structures and highlighting significant anomalies. Within the Tongling mining district this is demonstrated by the good correlation between areas highlighted during Sequential Gaussian Simulation and areas with known Pb-Zn mineralisation; this is not necessarily the case when using the raw soil Pb concentrations and when using data that has been Kriged. This indicates that, for the datasets with low spatial density sampling, calculation of a singularity exponent based on Sequential Gaussian Simulation could produce significantly improved results, and therefore improved interpretation, than using the data of raw or Kriged during the identification of anomalies within soil geochemical data.
Keywords：
土壤地球化学异常的提取和识别是勘查地球化学和环境地球化学工作重要目标[1, 2]。异常信息往往与背景信息相叠加,或是弱异常掩盖在强变异性的背景中,使得常规统计方法识别异常信息的效果并不显著,无法很好的应用于区域地球化学的异常识别[3,4]。多维分形理论的提出为地球化学数据的异常识别和处理提供了强有力的工具[1,3~6]。其中,基于多维分形理论的局部奇异性指数从多维分形的观点解释了尺度独立性,利用局部自相似性来度量奇异性,表明了空间变量随尺度变化的统计特征[4]。通过对这种统计特征进行研究,可以有效的应用于低缓、微弱、难识别异常信息的获取,并可以从背景中区分复杂的叠加异常[7]。近年来,针对局部奇异性指数的研究得到了不断的发展,基于空间U统计量[8]的广义局部奇异性指数计算方法[9]和局部奇异性分析迭代方法[10]均为局部奇异性指数计算的优化提供了新的途径。
在局部奇异性指数计算过程中,当勘查地球化学数据密度较为稀疏时,基于原始样本计算得到的奇异性指数很难展现局部奇异性的精细空间结构,通常需要对原始数据进行插值加密预处理。插值方法中,以地统计学为理论基础的克里格插值法由于其无偏最优特性[11],被广泛应用于地质、环境、生态等行业[12]。但克里格方法基于滑动加权平均原理,以变异函数为基础对局部范围进行加权平均,因而不可避免的对数据造成平滑效应[13]。而局部奇异性指数恰恰需要突出数据的差异和异常,克里格方法的平滑效应将显著降低数据集的高值和方差从而无法重建原始信息中的高值异常信息[14]。序贯高斯模拟方法同样属于地统计学范畴,与克里格方法相比,其理论上更适于定量刻画某一属性的非均质性和不确定性[15],得到的模拟值不但能保持与样本数据的数学期望、方差和分布函数相同,还能保持与协方差函数或变异函数的一致性[16],这在一定程度上弥补了克里格方法的平滑效应。
因此,本文以铜陵矿集区土壤Pb元素为实例,研究在稀疏采样条件下,序贯高斯模拟方法对奇异性指数计算的影响,并与原始样本、克里格方法进行对比,以期为更有效的计算奇异性指数提供理论和实践依据。
本文采用的勘查地球化学数据为安徽省铜陵矿集区土壤Pb元素含量数据。数据来源于安徽省地质调查院“安徽省江淮流域多目标区域地球化学调查”资料,土壤样品为较稳定地块的表层土壤,按间距为2 km稀疏网格采样（图 1）,样本数共计204件。
普通克里格法是一种基于样本数据的空间相关和分布性质的最优线性无偏估计（BLUE）方法[17]。变异函数作为度量数据空间相关和变异程度的工具,可为克里格插值提供空间变异性的定量化支持,离散化的变异函数可以由下式定义[18]：
式（1）中,h为样本点xi和xi+h的距离,N(h)为距离为h的点对数。
序贯高斯模拟假设一个高斯随机域,基于变异函数特征,对每一个待模拟的空间位置利用克里格方法求构建高斯分布函数,之后随机在其中抽取一个数值作为该点的模拟值[19]。序贯高斯模拟在一定程度上可以弥补克里格方法的平滑效应,得出的模拟值不但能保持与样本数据的数学期望、方差和分布函数的一致性,而且还能保持协方差函数或变异函数的一致性,同时在各实测点处的模拟值与该点的实测值一致。
序贯高斯模拟算法包括以下步骤[20]：① 对研究区内规则网格随机建立一条经过每一个节点的路径;② 对于每一个网格节点,基于变异函数,利用克里格方程组求取邻域内原始数据和已模拟数据的均值和方差;③ 利用克里格均值和方差建立高斯条件累计分布函数;④ 从高斯条件累计分布函数中获取一个随机值,作为该网格节点的模拟值,并将这个模拟值引入条件数据集参与之后随机路径上其他网格节点的模拟计算。
局部奇异性指数α基于多维分形理论[21],与传统统计学度量面积密度或非奇异性数据不同,奇异性指数从多维分形的角度刻画了场的分布特征,度量了异常的局部标度性和奇异性[22]。局部奇异性指数可以通过局部的自相似性来表征[13]：
ρ(ε)∝εα-2 (2)
式（2）中,ε代表尺度范围,ρ(ε)为尺度范围内的平均密度,α即为奇异性指数。
当奇异性指数α接近2时,表明该地段为正常的未受到成矿或外源影响太大的地段,这时元素密度基本不随分布范围缩小而变化;当奇异性指数α<2时,指示该地段受到成矿或外源影响而造成元素富集,这时元素密度随分布范围缩小而增大;当奇异性指数α> 2 时,指示该地段受到成矿或外源影响而造成元素贫化,这时元素密度随分布范围缩小而减小[23]。
局部奇异性指数求取常采用基于窗口方法[1, 3, 23]。对于待计算网格中任一个节点,定义一系列以网格节点为中心的不同尺度窗口,εmin=ε1<ε2<…<εn=εmax,之后通过计算每一个窗口内元素含量的平均密度值并与相应窗口尺度建立线性关系,求取奇异性指数,如下式：
log[ρ(ε)]=C+(α-2) log (ε) （3）
式（3）中,α-2即为该线性关系的斜率,可以通过最小二乘法求得。
由表1可见,Pb元素含量数据的最大值是最小值的83.85倍,变异系数达到152.92%,表明研究区内土壤Pb元素具有很强的空间异质性[24,25]。
由于普通克里格法和序贯高斯模拟方法均以变异函数为基础,数据的非正态分布将会影响到数据的平稳性,同时也会使实验变异函数产生比例效应,即提高基台和块金值,增大估计误差,使变差值点的波动变大,甚至会掩盖其固有的结构特征 [26,27]。同时,克里格法的线性最优的特性建立在数据集呈正态分布且符合二阶平稳假设或准二阶平稳假设之上[17],序贯高斯模拟方法则更需以数据场呈高斯分布作为模拟的前提[22]。
为了检验数据的正态分布特征,本文采用Kolmogorov-Smirnov（K-S）方法对原始数据进行正态检验[28]。本文设置信度α=0.05,若检验得到的p<0.05,则否定原假设,断定总体呈非正态分布。K-S正态验证方法显示数据明显不符合正态分布（p<0.01）。因此,本文对原数据进行了正态变换使数据满足正态分布要求（p>0.05）,从而更好的符合地统计学克里格及高斯序贯模拟方法的理论前提（表1）。
表1 土壤Pb元素含量
Table 1 Summary statistics of soil Pb concentrations
| 平均值 | 中位数 | 最小值 | 最大值 | 标准方差 | 变异系数 (%) | 偏度 | 峰度 | p1 | p2 | |
|---|---|---|---|---|---|---|---|---|---|---|
| Pb (μg/g) | 82.26 | 106 | 13 | 1090 | 125.79 | 152.92 | 4.19 | 23.72 | <0.01 | >0.15 |
图2可见,应用原始数据直接计算得到的变异函数趋势混乱、跃动明显,半变异函数先逐渐增加到一稳定值后又表现出降低的趋势,呈现不平稳的特性。而通过正态变换的数据变异函数则具有良好的空间平稳结构,呈现明显的球状理论变异函数特征,因此本文采用球状理论变异函数模型对实验变异函数进行拟合（表2）。
表2 变异函数参数
Table 2 Parameters of semi-variogram
| 块金值 | 基台值 | 变程 | 块金值/基台值 | 决定系数(r2) | 剩余平方和(RSS) |
|---|---|---|---|---|---|
| 0.217 | 1.203 | 20.32 | 0.18 | 0.994 | 5.150E-03 |
本文的克里格插值与序贯高斯模拟结果均采用GS+软件计算获得,数据的空间分布等值线图则通过ArcGIS9.2软件进行绘制。为保持数据的局部平稳性和对比的条件一致性,在克里格插值和序贯高斯模拟过程中统一设定最大搜索点数为12,插值网格间距为0.2 m;插值结果采用ArcGIS 软件中的分位点（Quantile）法进行分级,该方法可以保证各级包含有相同数量网格单元,从而更显著表达插值结果（图3）。
图3 土壤Pb元素空间分布克里格插值(a)、序贯高斯模拟(b)结果 （μg/g）
Fig.3 Spatial distribution map of soil Pb concentrations generated by Kriging and Sequential Gaussian Simulation estimation
由图3可见,克里格插值方法得到的空间分布图呈现相对平滑的空间分布特征,土壤Pb元素富集的区域虽然能与铅锌矿床的位置大体吻合,但范围较大,只能提供一定范围的富集趋势。而序贯高斯模拟的结果总体虽然相对破碎凌乱,但却显示出大量的细节信息,在一定程度上强调小范围内的空间不确定性。对比克里格插值和序贯高斯模拟结果的基本统计信息可见（表3）,克里格插值结果的最大值、标准方差以及变异系数较之原始数据均有很大程度降低,说明克里格插值具有非常强的平滑效应,忽略了部分高值变异特征和信息;而序贯高斯模拟方法的结果却呈现与原始数据一致的特征,从统计量的角度来看,序贯高斯模拟的各项统计量同原始数据更为一致。相比于克里格方法,序贯高斯模拟方法对于重建Pb元素的空间分布更有优势。
表3 估值结果基本统计信息
Table 3 Summary statistics of estimation results
| 平均值（μg/g） | 中位数（μg/g） | 最小值（μg/g） | 最大值（μg/g） | 标准方差 | 变异系数 (%) | 偏度 | 峰度 | |
|---|---|---|---|---|---|---|---|---|
| 原始数据 | 82.26 | 106 | 13 | 1090 | 125.79 | 152.92 | 4.19 | 23.72 |
| 克里格插值 | 51.85 | 29.81 | 19 | 531 | 51.02 | 98.39 | 3.32 | 17.62 |
| 序贯高斯模拟 | 79.64 | 37.35 | 13 | 1090 | 126.15 | 158.4 | 4.73 | 27.70 |
本文利用变异函数对克里格插值和序贯高斯模拟结果的空间结构特征进行分析（图4）。可见,所有数据都具有较好的结构特征,变程差异不大;但是克里格插值结果由于平滑作用导致方差和变异性变小,块金值和基台值都有了很大程度的下降,一定程度上改变了原始数据的空间结构;而序贯高斯模拟则与原始数据的变异函数曲线有着较好的一致性,特别是在分离距离较小的范围内,两者空间变异特征基本一致,较之克里格方法,序贯高斯模拟结果在空间结构上与原始数据更为相似。
本文根据样品的采集间距,分别设置边长ε=2,4,6,…,12 km的正方形计算不同尺度范围内Pb元素含量的平均密度ρ(ε),并在双对数图上利用最小二乘法求取每一点的奇异性指数α（图5）。
可见,由于采样间距较为稀疏,根据原始数据计算得到的奇异性指数分布图中出现大量规则的斑块状区域,斑块内奇异性指数大体相等,无法很好表征奇异性指数空间分布的精细结构。而基于克里格插值结果的奇异性指数对这种“斑块状”特征有较好改善,富集区域较为平滑,但由于克里格方法的平滑效应,其奇异性指数的数值范围（1.51~2.42）较原始数据（1.09~3.08）有大幅缩小,克里格方法不但对于原始高值数据有平滑效果,同时也降低了奇异性指数对比程度,这对于异常信息的识别和提取非常不利。基于序贯高斯模拟数据的奇异性指数的数值范围（1.39~2.61）较克里格方法有较大改善,提高了奇异性指数对比程度,并能够更精细刻画局部奇异性指数空间分布特征。同时,由于奇异性指数表征的是局部场的分布特征和空间自相似性,因此一定程度降低了序贯高斯模拟的破碎凌乱程度,对序贯高斯模拟结果也起到一定的改善作用。
为了定量分析比较基于上述3种数据集获取的奇异性指数的异常识别效果,本文分别对奇异性指数（α<2）富集范围内的已知矿床数量进行了统计（表4）。可见较之原始数据和克里格方法,基于序贯高斯模拟数据获得的奇异性指数富集区域（奇异性指数α<2）能识别出更多数量的已知矿床,最高富集等级范围对已知矿床的识别也更加显著。这均表明序贯高斯模拟可以更好的重建数据的高值异常特征,弥补克里格插值的平滑效应,在稀疏数据集条件下可以更好的进行异常信息识别,具有良好适用性。
表4 奇异性指数富集范围内已知矿床数统计
Table 4 Statistics of known deposits within enrichment area of singularity exponent
| α | 矿床数 | α（最高富集等级） | 矿床数 | |
|---|---|---|---|---|
| 原始数据 | <2 | 16 | <1.52 | 6 |
| 克里格插值 | <2 | 14 | <1.74 | 4 |
| 序贯高斯模拟 | <2 | 18 | <1.69 | 7 |
1) 对具有强变异性质的数据集,克里格插值结果的最大值、标准方差和变异系数较原始数据均大幅降低;而序贯高斯模拟结果的基本统计信息与原始数据更为一致,其强调短距离范围内的空间不确定性,可以弥补克里格方法平滑效应的不足,对于精细重建土壤元素的空间分布特征具有更好的效果。
2) 序贯高斯模拟与原始数据的变异函数曲线具有较好的一致性,特别是在分离距离较小的范围内,两者空间变异特征基本一致,较之克里格方法,序贯高斯模拟结果在空间结构上与原始数据更为相似。
3) 克里格方法不但对原始数据有平滑效果,也压缩了基于其计算得到的奇异性指数的数据范围,而基于序贯高斯模拟数据的奇异性指数较之克里格方法有较大改善,扩大数值范围、提高对比程度。对稀疏采样数据集,较之原始数据和克里格方法,基于序贯高斯模拟数据获取的奇异性指数可以更精细的刻画局部空间结构,突出局部异常。铜陵矿集区实例显示,序贯高斯模拟数据得到的土壤Pb元素富集区域与已知铅锌矿床的分布更加吻合,具有更好的异常识别效果。
The authors have declared that no competing interests exist.
| [1] |
Application of singularity mapping technique to identify local anomalies using stream sediment geochemical data, a case study from Gangdese,Tibet,western China [J]. |
| [2] |
Techniques for analysis and visualization of lithogeochemical data with applications to the Swayze Greenstone Belt,Ontario [J]. |
| [3] |
Mapping singularities with stream sediment geochemical data for prediction of undiscovered mineral deposits in Gejiu, Yunnan Province,China [J]. |
| [4] |
Multifractal IDW interpolation and fractal filtering method in environmental studies:an application on regional stream sediments of (Italy),Campania region [J]. |
| [5] |
Multifractal modeling and spatial statistics [J].
|
| [6] |
The Determination of Geochemical Abnormity of Soil Cd in Tongling Mining Area Using Multifractal Method [J]. |
| [7] |
成矿过程奇异性与矿产预测定量化的新理论与新方法 [J]. |
| [8] |
A spatial analysis method for geochemical anomaly separation [J]. |
| [9] |
多重分形局部奇异性分析方法及其在矿产资源信息提取中的应用[D] . |
| [10] |
A novel iterative approach for mapping local singularities from geochemical data [J]. |
| [11] |
Principles of geostatistics [J]. |
| [12] |
季节性冻土区不同时期土壤剖面水分空间变异特征研究 [J]. |
| [13] |
多重分形与地质统计学方法用于勘查地球化学异常空间结构和奇异性分析 [J].( |
| [14] |
多维分形克里格方法 [J]. |
| [15] |
农田土壤表层饱和导水率的条件模拟 [J]. |
| [16] |
序贯条件模拟方法研究及应用[D] . |
| [17] |
|
| [18] |
福建省经济空间增长变异特征及驱动机制 [J]. |
| [19] |
Characterization of soil lead by comparing sequential Gaussian simulation,simulated annealing simulation and kriging methods [J]. |
| [20] |
Using robust kriging and sequential Gaussian simulation to delineate the copper- and lead-contaminated areas of a rapidly industrialized city in Yangtze River Delta,China [J]. |
| [21] |
Multifractality and spatial statistics [J]. |
| [22] |
奇异性理论在个旧锡铜矿产资源预测中的应用:成矿弱信息提取和复合信息分解 [J].( |
| [23] |
非线性成矿预测理论:多重分形奇异性—广义自相似性—分形谱系模型与方法 [J].( |
| [24] |
区域尺度土壤养分空间变异及其影响因素研究 [J]. |
| [25] |
等 .淮海经济区县域经济差异变动的空间分析 [J]. |
| [26] |
区域景观破碎化的多尺度空间变异研究——以贵州省乌江流域为例 [J]. |
| [27] |
松嫩平原西部土壤盐碱化空间变异与微地形关系研究 [J]. |
| [28] |
Using sequential indicator simulation to assess the uncertainty of delineating heavy-metal contaminated soils [J]. |
/
| 〈 |
|
〉 |
