CN102254095A - 基于多维分形克里格方法的成矿异常提取方法 - Google Patents

Abstract

本发明公开了一种基于多维分形克里格方法的成矿异常提取方法，其特征是首先对待分析勘查地球化学数据集进行正态分布检验，对于服从正态分布的数据集计算实验变异函数值，并通过理论变异函数模型对实验变异函数值进行拟合，获得理论变异函数模型中的相关参数；在重建局部高值信息和突出弱异常的基础上更有效地基于多维分形克里格插值结果对异常下限值进行求取，从而获得更为可靠的成矿异常区域用于找矿勘探工作。

Description

基于多维分形克里格方法的成矿异常提取方法

技术领域

本发明涉及一种基于多维分形克里格方法的成矿异常提取方法，更具体地说是通过对勘查地球化学数据进行多维分形克里格方法分析，进而提取地球化学成矿异常区域，从而为准确圈定勘探靶区提供依据。

背景技术

勘查地球化学方法是矿产勘查的重要手段，通过对岩石、水系沉积物、水、土壤等原生和次生介质样品进行分析，可以有效地对区域地球化学背景和异常进行调查研究，从而为矿产资源勘查和环境评价等工作起到重要的指导和支持作用。在针对地球化学采样数据的地球化学场的异常识别工作中，为了研究数据的空间变化趋势以及异常范围，通常需要对采样点数据进行空间插值处理使其转化为连续的曲面。常见的插值方法包括反距离加权插值法(IDW)、克里格法(Kriging)、径向基函数法以及样条曲线法(Spline)等。这些方法的往往体现某种趋势或是某种平均估计结果，具有显著的低通平滑特点。虽然克里格插值方法通过变异函数可以表征一定尺度的局部变化性，但其作为滑动加权插值法的一种，其低通滤波的特性使其只适合于讨论均值周围的变化性，却难以反映高值异常信息，而这些局部高值异常信息却是资源预测以及找矿勘探最为关注的。

因此，如何刻画和重建局部高值变异信息已成为提高勘查地质学、勘查地球化学处理效果的重要途径之一，尤其是对于以矿产勘查为目的的地球化学场数据处理和异常分析，如何保持和突出与成矿有关的局部异常是数据处理和分析成败的关键；同时，在复杂的地质过程尤其是成矿作用过程中，地球化学场往往受到多期次地质过程的影响，造成异常场同多种不同的背景场相互叠加，弱异常常被掩盖在强变异性的背景中，因此如何有效地从复杂的地球化学场中分离出弱异常，这对于成矿信息提取和成矿预测也至关重要。

为了能同时表达勘查地球化学场的局部空间自相关性以及奇异性，Cheng(1999)基于非线性理论提出了多维分形插值方法。该方法通过线性插值或是地统计学变异函数来表达空间自相关性，同时利用多重分形模型对具有自相似性或统计自相似性的多重分形分布(multifractal distributions)的奇异性(α)进行刻画。多维分形插值方法能同时保持和突出场值的局部空间结构和奇异性信息，这对于保持和增强数据的局部结构信息、重建高值和提取弱异常信息都非常有利。多维分形插值方法的公式如下：

$Z\left(x_0\right)={\mathrm{\ϵ}}^{a\left(x_0\right)-2}\underset{\mathrm{\Ω}(x_0,\mathrm{\ϵ})}{\mathrm{\Σ}}\mathrm{\ω}\left(\right||x-x_0|\left|\right)Z\left(x\right)---\left(3\right)$

式中，ε表示分形估值尺度，α(x₀)表示x₀处的奇异性指数，ω为滑动平均加权的权系数。

多维分形插值是对待插值点的小邻域内的均值进行邻域奇异性校正的结果，而克里格方法虽然具有明显的低通平滑效应，但却是估计小邻域均值的最好方法之一，因此将多维分形方法与完善的克里格方法相结合，可以融合克里格和多维分形方法的共同优点。

由多维分形插值公式(3)可知，分形估值尺度ε可以决定奇异性因子对插值局部的强化和减弱的程度，从而对插值的效果产生直接的影响。用一个较大范围内的均值估计较小范围内的均值时，由于奇异性校正的因素，将不可避免的引入误差，这是由问题本身未知因素过多造成的。因此，具有不确定性的分形估值尺度ε直接影响到成矿异常的提取效果。

此外，目前地球化学勘查标准中的异常下限值计算方法(平均值

与2倍标准差σ之和法)是基于正态分布数据的计算方法，这种方法具有很大的局限性，对明显具有强变异性质且呈非正态分布的矿区土壤地球化学数据难以处理。由于多维分形方法对原数据进行了多维分形校正，因此更难采用传统方法对成矿异常区域进行识别。

发明内容

本发明是为避免上述现有技术所存在的不足之处，提供的一种基于多维分形克里格方法的成矿异常提取方法和流程。以期能够获取最优的多维分形克里格插值结果，在重建局部高值信息和突出弱异常的基础上更有效地基于多维分形克里格插值结果对异常下限值进行求取，从而获得更为可靠的成矿异常区域用于找矿勘探工作。

本发明解决技术问题采用如下技术方案：

本发明基于多维分形克里格方法的成矿异常提取方法的特点是按如下步骤进行：

(1)、对已经获得的勘查区内待分析勘查地球化学数据集进行正态分布检验，如果数据集服从正态分布则直接进入步骤(2)；如果数据集不服从正态分布，则采用剔除所述数据集中的异常值或对所述数据集进行正态变换的方式使所述数据集服从正态分布条件，在所述数据集服从正态分布条件后进入步骤(2)；

(2)、对于来自步骤(1)的服从正态分布的数据集，计算所述数据集的实验变异函数值，并通过理论变异函数模型对实验变异函数值进行拟合，获得理论变异函数模型中的相关参数，离散化的实验变异函数定义为：

$\mathrm{\γ}\left(h\right)=\frac{1}{2N\left(h\right)}\underset{i=1}{\overset{N\left(h\right)}{\mathrm{\Σ}}}{[Z\left(x_i\right)-Z(x_i+h)]}^2---\left(1\right)$

式(1)中，h为滞后距离，N(h)为距离等于h的点对数，Z(x_i)为点x_i处变量的真实值，Z(x_i+h)为与点x_i偏离h处变量的真实值；

(3)、对勘查区进行规则格网划分；结合理论变异函数及其拟合获取的理论变异函数模型中的相关参数，利用地统计学克里格插值方法对已服从正态分布的数据集进行克里格插值计算；如果所述勘察地球化学数据集经过正态变换，则须对插值结果进行逆变换处理；对克里格插值结果进行交叉验证，通过交叉验证结果对理论变异函数模型中的相关参数进行修正优化，得到最优的克里格插值结果以及所述数据集中任一数据点的克里格交叉验证结果；

(4)、结合所述待分析勘查地球化学数据集，对勘查区内格网的任一格网点进行奇异性指数计算；所述奇异性指数计算方法是：首先定义不同尺度r的窗体，r_min＝r₁＜r₂...＜r_n＝r_max，分别计算不同尺度r窗体下的元素浓度平均值C[A(r_i)]；r_i和C[A(r_i)]服从式(2)的线性关系：

Log C[A(r_i)]＝C+(α-2)Log(r)               (2)

由式(2)所述线性关系，对勘查区内的任一格网点，将不同尺度r下元素平均浓度C[A(r)]和不同尺度r在双对数坐标下进行线性拟合，由式(2)计算得到任一格网点的奇异性指数α值；

(5)、利用奇异性指数与克里格交叉验证结果对(3)式中的分形估值尺度ε进行优化选择，选取不同的分形估值尺度ε，运用多维分形方法对克里格交叉验证结果进行多维分形插值方法校正，所述多维分形方法公式为：

$Z\left(x_0\right)={\mathrm{\ϵ}}^{a\left(x_0\right)-2}\underset{\mathrm{\Ω}(x_0,\mathrm{\ϵ})}{\mathrm{\Σ}}\mathrm{\ω}\left(\right||x-x_0|\left|\right)Z\left(x\right)---\left(3\right)$

式(3)中，ε表示分形尺度，α(x₀)表示x₀处的奇异性指数，ω为滑动平均加权的权系数；

将不同分形估值尺度ε下的克里格交叉验证结果的多维分形校正结果同所述勘查数据集中的实际真实数据进行比较，采用绝对平均误差作为衡量分形估值尺度ε优劣的统计指标，以所述绝对平均误差最小的分形估计尺度ε做为最佳的分形估计尺度；所述绝对平均误差计算公式如下：

$E=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}|x_n-y_n|}{N}---\left(4\right)$

式(4)中，x_n为所述数据集中第n个数据点的克里格交叉验证结果的多维分形校正结果，y_n为所述勘查数据集中第n个数据点的实际真实数据，N为样品总数，E为绝对平均误差；

(6)、结合最优的分形估值尺度、克里格插值结果以及奇异性指数计算结果，利用(3)式所述的多维分形方法公式对勘查区内所有格网点进行计算，获得多维分形克里格插值结果；

(7)、利用步骤(6)获得的多维分形克里格插值结果绘制等值线图，在等值线图的基础上基于V-N元素含量-矿床数目累积频率模型对成矿异常区域进行提取和识别，所述V-N元素含量-矿床数目累积频率模型为：

$N\left(r\right)=\underset{r_a\≤r\≤r_b}{\mathrm{\Σ}}{t}_{\left(r\right)}---\left(5\right)$

式(5)中：r为元素含量值，r_a≤r≤r_b；r_a为能够包含勘查区内全部相关矿床所处空间位置的等值线含量值，r_b为不能包含勘查区内任何相关矿床所处空间位置的等值线含量值，∑t(r)为含量值为r的等值线所包含空间范围内的相关矿床数量；

应用(5)式需确定某一含量值r_a，该含量值等值线所包含的空间范围可以囊括勘查区内全部相关矿床所处的空间位置，从含量值为r_a的等值线开始，以等距递增，直到某一含量值等值线不能包含勘查区内任何相关矿床的空间位置为止；分别计算不同含量值等值线范围内囊括的相关矿床数量N(r)，并将含量值r与矿床数量N(r)绘制成散点图；之后采用最小二乘法对所述散点图进行分段线性拟合，得到若干交叉分界点；将各个交叉分界点所代表的含量值在等值线图中进行异常区域圈定，选择最具有异常指示效果的交叉分界点作为异常下限值；

(8)、将步骤(7)中获得的异常下限值在多维分形克里格插值结果等值线图中进行空间圈定，最终圈定获得的区域即为成矿异常区域。

与已有技术相比，本发明有益效果体现在：

1、通过本发明方法建立的基于多维分形克里格方法的成矿异常识别方法和流程，能够获取最优的多维分形克里格插值结果，并针对多维分形克里格插值结果更为有效的对异常下限进行计算并对异常区域进行提取和圈定。

2、本发明方法在步骤5中对奇异性指数和克里格交叉验证结果进行综合分析，采用绝对平均误差作为评价指标获取了最优的分形估值尺度，增加了多维分形克里格插值的确定性和准确性。

3、本发明方法在步骤7中应用V-N元素含量-矿床数目累积频率模型，较之已有技术方法，可以更为适用和准确的基于多维分形克里格插值结果进行异常区域的提取和圈定，获取更为可信的和潜在的成矿异常区域，为准确圈定勘探靶区提供更为可靠的依据。

附图说明

图1为本发明基于多维分形克里格方法的成矿异常提取方法流程图；

具体实施方式

本实施例基于多维分形克里格方法的成矿异常提取方法按如下过程进行：

(1)、对已经获得的勘查区内待分析勘查地球化学数据集进行正态分布检验，如果数据集服从正态分布则直接进入步骤(2)；如果数据集不服从正态分布，则需采用剔除数据集中的异常值或对数据集进行正态变换的方式使数据集服从正态分布条件，在数据集服从正态分布条件后进入步骤(2)；

正态变换采用对数变换、Box-Cox变换或Johnson变换；

(2)、对于来自步骤(1)的服从正态分布的数据集，计算其实验变异函数值，并通过理论变异函数模型对实验变异函数值进行拟合，获得理论变异函数模型中的相关参数，离散化的实验变异函数定义为：

$\mathrm{\γ}\left(h\right)=\frac{1}{2N\left(h\right)}\underset{i=1}{\overset{N\left(h\right)}{\mathrm{\Σ}}}{[Z\left(x_i\right)-Z(x_i+h)]}^2---\left(1\right)$

式(1)中，h为滞后距离，N(h)为距离等于h的点对数，Z(x_i)为点x_i处变量的真实值，Z(x_i+h)为与点x_i偏离h处变量的真实值；

理论变异函数模型包括：球状模型、指数模型、高斯模型、幂函数模型、线性模型和纯块金模型，理论变异函数具有不同的形状特征，需针对实验变异函数的形状特征进行选择，常用于土壤地球化学数据分析的理论变异函数有球状模型、线性模型和纯块金模型，具体形式如下：

球状模型为：

$\mathrm{\γ}\left(h\right)=\left\{\begin{array}{cc}{C}_0+C_1[1.5(h/a)-0.5{(h/a)}^3],& 0\≤h\≥a\\ C_0+C_1,& h>a\end{array}\right.$

线性模型为：

$\mathrm{\γ}\left(h\right)=\left\{\begin{array}{cc}{C}_0+C_1\left[h\right],& 0\≤h\≥a\\ C_{1}a,& h>a\end{array}\right.$

纯块金模型为：

$\mathrm{\γ}\left(h\right)=\left\{\begin{array}{cc}0,& h=0\\ C_0,& h>a\end{array}\right.$

式中C₀为块金值；C₀+C₁为基台值；a为变差距离或相关尺度；

(3)、对勘查区进行规则格网划分；结合理论变异函数及其拟合获取的理论变异函数模型中的相关参数，利用地统计学克里格插值方法对已服从正态分布的数据集进行克里格插值计算；如果勘察地球化学数据集经过正态变换，则须对插值结果进行逆变换处理；对克里格插值结果进行交叉验证，通过交叉验证结果对理论变异函数模型中的相关参数进行修正优化，得到最优的克里格插值结果以及数据集中任一数据点的克里格交叉验证结果；

交叉验证方法是指在克里格插值的过程中，在每一有实际真实样品数值的地方，将此实际真实样品数值暂时去除，用其他剩余的观察值以选定好的变异函数模型和参数采用克里格法对每一个观察点进行估计，并得到一系列估计值；通过对实际值和估计值进行比较可以对理论变异函数模型以及克里格插值的精度做出评价，继而对理论变异函数模型的相关参数进行修正优化，最后得到具有最佳精度的克里格插值结果。

(4)、结合待分析勘查地球化学数据集，对勘查区内格网的任一格网点进行奇异性指数计算；奇异性指数计算方法是：首先定义不同尺度r的窗体，r_min＝r₁＜r₂...＜r_n＝r_max，分别计算不同尺度r窗体下的元素浓度平均值C[A(r_i)]；r_i和C[A(r_i)]服从式(2)的线性关系：

Log C[A(r_i)]＝C+(α-2)Log(r)                        (2)

由式(2)线性关系，对勘查区内的任一格网点，将不同尺度r下元素平均浓度C[A(r)]和不同尺度r在双对数坐标下进行线性拟合，由式(2)计算得到任一格网点的奇异性指数α值；

(5)、利用奇异性指数与克里格交叉验证结果对(3)式中的分形估值尺度ε进行优化选择，选取不同的分形估值尺度ε，运用多维分形方法对克里格交叉验证结果进行多维分形插值方法校正，多维分形方法公式为：

$Z\left(x_0\right)={\mathrm{\ϵ}}^{a\left(x_0\right)-2}\underset{\mathrm{\Ω}(x_0,\mathrm{\ϵ})}{\mathrm{\Σ}}\mathrm{\ω}\left(\right||x-x_0|\left|\right)Z\left(x\right)---\left(3\right)$

式(3)中，ε表示分形尺度，α(x₀)表示x₀处的奇异性指数，ω为滑动平均加权的权系数。

如果数据点x₀位于背景场和非奇异场时，奇异性指数α(x₀)＝2，通过多维分形方法计算的结果与通常的加权平均插值结果相同；而当数据点处于含量富集地段而且局部场具有奇异性时，奇异性指数α(x₀)＜2，通过多维分形方法计算的结果将高于加权平均插值结果；相反，当处于含量贫化地段时，奇异性指数α(x₀)＞2，通过多维分形方法计算的结果将低于加权平均插值结果。

将不同分形估值尺度ε下的克里格交叉验证结果的多维分形校正结果同勘查数据集中的实际真实数据进行比较，采用绝对平均误差作为衡量分形估值尺度ε优劣的统计指标，以绝对平均误差最小的分形估计尺度ε做为最佳的分形估计尺度；绝对平均误差计算公式如下：

$E=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}|x_n-y_n|}{N}---\left(4\right)$

式(4)中，x_n为数据集中第n个数据点的克里格交叉验证结果的多维分形校正结果，y_n为数据集中第n个数据点的真实数据，N为样品总数，E为绝对平均误差；

(6)、结合最优的分形估值尺度、克里格插值结果以及奇异性指数计算结果，利用(3)式的多维分形方法公式对勘查区内所有格网点进行计算，获得多维分形克里格插值结果；

(7)、利用步骤(6)获得的多维分形克里格插值结果绘制等值线图，在等值线图的基础上基于V-N元素含量-矿床数目累积频率模型对成矿异常区域进行提取和识别，V-N元素含量-矿床数目累积频率模型为：

$N\left(r\right)=\underset{r_a\≤r\≤r_b}{\mathrm{\Σ}}{t}_{\left(r\right)}---\left(5\right)$

式(5)中：r为元素含量值，r_a≤r≤r_b；r_a为能够包含勘查区内全部相关矿床所处空间位置的等值线含量值，r_b为不能包含勘查区内任何相关矿床所处空间位置的等值线含量值，∑t(r)为含量值为r的等值线所包含空间范围内的相关矿床数量；

应用式(5)需确定某一含量值r_a，该含量值等值线所包含的空间范围可以囊括勘查区内全部相关矿床所处的空间位置，从含量值为r_a的等值线开始，以等距递增，直到某一含量值等值线不能包含勘查区内任何相关矿床的空间位置为止；分别计算不同含量值等值线范围内囊括的相关矿床数量N(r)，并将含量值r与矿床数量N(r)绘制成散点图；之后采用最小二乘法对散点图进行分段线性拟合，得到若干交叉分界点；将各个交叉分界点所代表的含量值在等值线图中进行异常区域圈定，选择最具有异常指示效果的交叉分界点作为异常下限值；

(8)、将步骤(7)中获得的异常下限值在多维分形克里格插值结果等值线图中进行空间圈定，最终圈定获得的区域即为成矿异常区域。

Claims (1)

1.基于多维分形克里格方法的成矿异常提取方法，其特征是按如下步骤进行：

(1)、对已经获得的勘查区内待分析勘查地球化学数据集进行正态分布检验，如果数据集服从正态分布则直接进入步骤(2)；如果数据集不服从正态分布，则采用剔除所述数据集中的异常值或对所述数据集进行正态变换的方式使所述数据集服从正态分布条件，在所述数据集服从正态分布条件后进入步骤(2)；

(2)、对于来自步骤(1)的服从正态分布的数据集，计算所述数据集的实验变异函数值，并通过理论变异函数模型对实验变异函数值进行拟合，获得理论变异函数模型中的相关参数，离散化的实验变异函数定义为：

$\mathrm{\γ}\left(h\right)=\frac{1}{2N\left(h\right)}\underset{i=1}{\overset{N\left(h\right)}{\mathrm{\Σ}}}{[Z\left(x_i\right)-Z(x_i+h)]}^2---\left(1\right)$

式(1)中，h为滞后距离，N(h)为距离等于h的点对数，Z(x_i)为点x_i处变量的真实值，Z(x_i+h)为与点x_i偏离h处变量的真实值；

(3)、对勘查区进行规则格网划分；结合理论变异函数及其拟合获取的理论变异函数模型中的相关参数，利用地统计学克里格插值方法对已服从正态分布的数据集进行克里格插值计算；如果所述勘察地球化学数据集经过正态变换，则须对插值结果进行逆变换处理；对克里格插值结果进行交叉验证，通过交叉验证结果对理论变异函数模型中的相关参数进行修正优化，得到最优的克里格插值结果以及所述数据集中任一数据点的克里格交叉验证结果；

(4)、结合所述待分析勘查地球化学数据集，对勘查区内格网的任一格网点进行奇异性指数计算；所述奇异性指数计算方法是：首先定义不同尺度r的窗体，r_min＝r₁＜r₂...＜r_n＝r_max，分别计算不同尺度r窗体下的元素浓度平均值C[A(r_i)]；r_i和C[A(r_i)]服从式(2)的线性关系：

Log C[A(r_i)]＝C+(α-2)Log(r)                       (2)

由式(2)所述线性关系，对勘查区内的任一格网点，将不同尺度r下元素平均浓度C[A(r)]和不同尺度r在双对数坐标下进行线性拟合，由式(2)计算得到任一格网点的奇异性指数α值；

(5)、利用奇异性指数与克里格交叉验证结果对(3)式中的分形估值尺度ε进行优化选择，选取不同的分形估值尺度ε，运用多维分形方法对克里格交叉验证结果进行多维分形插值方法校正，所述多维分形方法公式为：

$Z\left(x_0\right)={\mathrm{\ϵ}}^{a\left(x_0\right)-2}\underset{\mathrm{\Ω}(x_0,\mathrm{\ϵ})}{\mathrm{\Σ}}\mathrm{\ω}\left(\right||x-x_0|\left|\right)Z\left(x\right)---\left(3\right)$

式(3)中，ε表示分形尺度，α(x₀)表示x₀处的奇异性指数，ω为滑动平均加权的权系数；

将不同分形估值尺度ε下的克里格交叉验证结果的多维分形校正结果同所述勘查数据集中的实际真实数据进行比较，采用绝对平均误差作为衡量分形估值尺度ε优劣的统计指标，以所述绝对平均误差最小的分形估计尺度ε做为最佳的分形估计尺度；所述绝对平均误差计算公式如下：

$E=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}|x_n-y_n|}{N}---\left(4\right)$

式(4)中，x_n为所述数据集中第n个数据点的克里格交叉验证结果的多维分形校正结果，y_n为所述勘查数据集中第n个数据点的实际真实数据，N为样品总数，E为绝对平均误差；

(6)、结合最优的分形估值尺度、克里格插值结果以及奇异性指数计算结果，利用(3)式所述的多维分形方法公式对勘查区内所有格网点进行计算，获得多维分形克里格插值结果；

(7)、利用步骤(6)获得的多维分形克里格插值结果绘制等值线图，在等值线图的基础上基于V-N元素含量-矿床数目累积频率模型对成矿异常区域进行提取和识别，所述V-N元素含量-矿床数目累积频率模型为：

$N\left(r\right)=\underset{r_a\≤r\≤r_b}{\mathrm{\Σ}}{t}_{\left(r\right)}---\left(5\right)$

式(5)中：r为元素含量值，r_a≤r≤r_b；r_a为能够包含勘查区内全部相关矿床所处空间位置的等值线含量值，r_b为不能包含勘查区内任何相关矿床所处空间位置的等值线含量值，∑t(r)为含量值为r的等值线所包含空间范围内的相关矿床数量；

应用(5)式需确定某一含量值r_a，该含量值等值线所包含的空间范围可以囊括勘查区内全部相关矿床所处的空间位置，从含量值为r_a的等值线开始，以等距递增，直到某一含量值等值线不能包含勘查区内任何相关矿床的空间位置为止；分别计算不同含量值等值线范围内囊括的相关矿床数量N(r)，并将含量值r与矿床数量N(r)绘制成散点图；之后采用最小二乘法对所述散点图进行分段线性拟合，得到若干交叉分界点；将各个交叉分界点所代表的含量值在等值线图中进行异常区域圈定，选择最具有异常指示效果的交叉分界点作为异常下限值；

(8)、将步骤(7)中获得的异常下限值在多维分形克里格插值结果等值线图中进行空间圈定，最终圈定获得的区域即为成矿异常区域。

Images