基于栅格数据的局部奇异性分析迭代方法及其软件实现
2014-08-02，,3，
，,3，
(1.中国地质大学地质过程与矿产资源国家重点实验室，湖北武汉430074； 2.中国地质大学(武汉)资源学院，湖北武汉430074； 3.York大学地球空间科学系和地理系，加拿大多伦多M3J1P3)
基于栅格数据的局部奇异性分析迭代方法及其软件实现
陈志军1,2，成秋明1,2,3，陈建国1,2
(1.中国地质大学地质过程与矿产资源国家重点实验室，湖北武汉430074； 2.中国地质大学(武汉)资源学院，湖北武汉430074； 3.York大学地球空间科学系和地理系，加拿大多伦多M3J1P3)
近些年，弱缓化探异常识别已成为成矿预测和勘查评价中十分关键的问题。多重分形理论的局部奇异性分析因其崭新的思路、简便的方法、优良的效果而在弱缓异常识别中受到广泛关注。在深入剖析局部奇异性分析基本思想及计算方法的基础上，引入正规化尺度参数L，提出了迭代方法来改进局部奇异性指数的估值。给出了奇异性迭代算法并用C++编程实现，软件功能强劲，操作灵活，不仅适用于各向同性情形，还适用于各向异性尺度的奇异性计算和方向性奇异性计算。软件的动态链接库版本已在GeoDAS矿产资源定量预测专业软件中应用。
局部奇异性分析; 栅格数据; 多重分形; 迭代方法; 正规化尺度参数
0 引 言
分形与多重分形在勘查地球化学数据、勘查地球物理数据、矿石品位数据等地质勘查数据中的普遍存在性为成矿过程奇异性与矿产预测定量化的新理论与新方法研究提供了崭新的思路和新型的数学工具(成秋明, 2007; Agterberg, 2003)。不均匀性、自相似性(自放射性、广义自相似性)、奇异性等通常被用来描述地球化学场特征、矿床空间分布、热液成矿系统矿床品位-吨位分形模型的非线性性质。一种五参数的多重级联模型最近被提出，该模型对de Wijs’s锌数据等多重分形经典数据的非对称奇异谱进行了成功建模，为现实世界中地球化学场多重分形性的存在性提供了理论依据。分形维数、多重分形对称度、多重分形奇异谱等参数从非线性角度对描述成矿物质富集与贫化的统计规律发挥了重要作用(成秋明, 2000, 2007)，然而，对于矿产资源定量预测中的“定位”(矿产资源空间分布位置)预测功能则有所不足，需要将全局性统计转变为局部化模型。基于多重分形理论建立的局部奇异性分析模型，试图通过局部异常的有效识别来达到空间预测目的，该模型正逐渐在覆盖区勘查地球化学弱缓异常识别中表现出巨大的应用潜力，成为非线性矿产预测的有效方法之一(成秋明, 2012)。应用的需求也促进了局部奇异性分析自身模型的发展，出现了多种形式的扩展模型(陈志军, 2007)。主要从提高奇异性指数计算精度方面介绍局部奇异性分析迭代分析方法，设计算法并编程实现，同时指明准确计算奇异性指数需要注意的一些关键参数设置。
1 局部奇异性分析方法的基本思想
奇异性理论指出：奇异性是指在很小的时间-空间范围内具有巨大能量释放或巨量物质形成的现象(成秋明, 2006; Cheng et al, 2009)。假定在E维(拓扑)空间中有一由质点凝聚成分维数为D的分形(如矿物的富集)，L尺度上质点集合的质量M(L)与集合中质点的数目N(L)成正比，则该分形体的密度可表示为(陈颙等, 2005)：
(1)
由式(1)可知，当D 该模型的基本内容：记B(x,ε)为任意空间位置x，尺度大小为ε的盒子(可以是方形、矩形、或其他更复杂的形状)，在盒子B(x,ε)覆盖区域内的金属量为μ(B(x,ε))，则有： μ(B(x,ε))=c(x,ε)εα(x) (2) 为便于理解与制图，令： Δα(x) =E-α(x) (3) 设ε大小邻域上的平均密度为〈ρ(B(x,ε))〉，〈〉表示统计期望，式(2)于是写成密度形式，有： 〈ρ(B(x,ε))〉 =μ(B(x,ε)) /εE=c(x,ε)ε-Δα(x) (4) 式(4)中，E为研究对象的空间拓扑维数(E=1,2,3，分别对应一维、二维或三维问题)，待定系数α(x)称为空间位置x处的局部奇异性指数，它表征了奇异性的强弱程度，是与尺度无关的量，无量纲；若ε趋近于0且c(x,ε)存在，则此时的c值被认为是位置x处的“α(x)维分形密度”。若ε用最大尺度进行了归一化，则c(x,ε)值是对该最大尺度的滑动平均的逼近。对地球化学采样数据，测量所得数据通常表示的是浓度含量，可视为一种广义的“密度”数据。因此，利用式(4)可直接对化探测量原始数据建模，对结果解释更直观，故而在分析中式(4)比式(3)更常用。需要注意的是，式(3)也是非常有用的，特别是在奇异性指数的算法设计方面。式(3)在计算中可直接获得α值并且其拟合优度高于用式(4)直接获得的Δα。这是由于普通最小二乘(OLS)模型只考虑纵向误差，而α的均值水平趋近于E、Δα的均值水平常趋近于0的缘故。(-Δα)作为log-log坐标系中回归直线的斜率，若其结果近似为0，意味着自变量和应变量不存在相关关系；同样的数据用α来表示(此时近似为E)，自变量和应变量的相关关系将变得显著。无论用式(3)还是式(4)计算，获得的奇异性指数和局部系数是一致的，不同的是log-log坐标系下对应的回归方程具有不同的显著性水平。要解决此问题，对式(4)和式(3)在log-log坐标系下可选用标准加权最小二乘法(SWLS)，权重为1/(1+b2)(其中b为拟合直线斜率)，以观测点到回归直线的垂直距离平方和最小化来确定直线参数，当然该计算过程比普通最小二乘模型略显复杂。 奇异性指数α值(或Δα)有着优良的异常指示功能，在实际应用中常用来度量不同空间位置物质(或能量)较其周围相对富集或亏损的程度。α值与E偏差越大，说明局部富集或亏损越越强烈，也即奇异性越强。以二维化探异常研究为例，此时E取值为2，当α<2(即Δα>0)表示局部正异常，元素相对局部富集；α>2(即Δα<0)表示局部负异常，元素相对局部亏损；而α= 2(即Δα=0)表示无异常，属于背景。通常在全图范围内计算各个位置的α值，对α值(或Δα值)进行制图表达，从而来识别不同空间位置的局部化探异常，并进一步筛选出有用异常并在空间上圈定出来。没有奇异性的区域(α=E)所构成子集的分形维数接近于分形盒维数E，对于地球化学背景区域，它们的空间分布通常占全图的绝大部分，其统计分布规律通常符合正态分布或对数正态分布；而奇异性区域(α≠E)所构成子集的分形维数有很多个，可以由f(α) 利用奇异性指数识别化探异常的效果经常被用来与衬值异常进行比较。陈志军等(2009)探讨了α值和滑动衬值(简记为CV)之间的定量关系。当无标度区间上符合严格的幂律关系时，α的估值可任选2个尺度进行计算，有如下关系成立： (5) 式(5)中，la和lb分别表示较小和较大2个不同尺度的窗口；Δα和log(CV)相比较而非直接比较CV，在统计上更能区分辨别二者的异常能力。由于对数函数是典型的单调函数，log(CV)和CV具有相同的排序值，可以通过比较Δα和CV的排序值(或用百分位数)来判别异常识别能力(陈志军, 2007)。 局部奇异性指数的计算通常要考察多个尺度的趋势性变化，由统计分形关系来确定，在log-log双对数图中采用回归分析技术来估值，相对仅由2个窗口度量值所确定的滑动衬值来说，Δα值对局部异常强度的指示作用更稳健，更能代表局部范围内异常的总体水平。滑动衬值则常因2个窗口选择的随意性而易受波动。通常，局部奇异性指数相比滑动衬值异常识别能力更佳(陈志军等, 2009)。 2.1 栅格数据模型 从计算简便性考虑，二维情形下基于栅格数据模型的计算方案效率最高。在栅格形式表达中，虽然对于栅格的形状没有限制(如等三角形、等六边形等)，然而正方形网格的简单性、运算方便性、直观性使得这种形式的栅格数据模型最为普遍。这些形状一样的正方形网格作为基本单元(单元格或称像元)组成了形如数学中规则排列的矩阵。像元大小、空间位置及像元值构成了栅格数据的基本要素。像元大小表示数据的最小尺度，或称空间分辨率；像元的取值即为空间对象的测量值；每个像元有唯一的行索引和列索引，测量值的空间位置隐含于栅格行列位置，也即可根据行列号转换成相应的空间坐标。一个栅格数据集描述了某个专题内容在区域的位置、特性和空间上的相对位置。 众多GIS软件都支持栅格数据类型，并且有众多数据格式。表1展示了ArcGIS和Surfer 2种主流的栅格数据文本文件格式，由此可以直观地了解栅格数据的编码规则。在对这两类栅格数据文件进行相互转化中，需要注意ArcGIS ASC文件所定义的左下角位置是左下角栅格单元的左下边界的角点，而Surfer 6 Text Grid栅格数据所定义的X最小值和Y最小值位于栅格单元的中心位置，因此两者偏差半个像元大小(表1示例文件中该偏差值为50)。ArcGIS ASC文件文件头信息可自定义缺失数据标记值，如-99 999；而在Surfer 6 Text Grid文件中，缺失数据默认标记为1.701 41E+038。此外，两者数据阵列的存放方案是不一样的，具有上下对称的关系。 表1 ArcGIS和Surfer栅格数据文件格式示例 2.2 局部奇异性分析算法中的量纲问题 局部奇异性分析算法通常基于栅格数据的多窗口邻域统计来开展计算。利用栅格数据在邻域分析方面的简便性，可以快速获得各个空间位置上不同尺度的滑动平均值，进而在log-log坐标下进一步计算奇异性指数(空间维数E=2)。栅格数据模型具有统一的几何支撑，意味着像元分辨率范围内研究对象具有同质性。对于某个栅格数据来说，空间分辨率一旦确定，局部奇异性特征尺度区间的最小值也就随之确定。因此，获取的局部奇异性指数是在一定尺度区间上估计出来的，揭示的是统计分形规律，在现实中特征尺度不能无限趋近于0。 局部奇异性分析方法算法的一些关键计算技术如：空间覆盖系列盒子的各向同性、各向异性窗口的多种构造方式、空间加权处理、掩模矩阵构置、等效尺度的确定及log-log最小二乘线性拟合、数据边缘区扩边、多种测度参见文献(陈志军，2007)。为便于推导迭代算法，这里有必要分析一下尺度ε的处理问题。ε的量纲是否及如何影响计算结果？假定窗口由小到大依次递增序列为：lmin=l1 2.2.1ε=l具有长度量纲 此时，c值在理论上的量纲是[ρ][l]Δα，估计Δα值(或α值)和c值的回归方程可表示为： log〈ρ(x,li)〉 = log(c(x))-Δα(x)×log(li),i=1,2,3，…,m (6) 用α(x)表示，也即 log〈ρ(x,li)×li2〉=log(c(x))+α(x)×log(li)，i=1,2,3，…,m (7) 其中，Δα(x)和log(c(x))为待定系数。式(6)与(7)相比，式(6)在形式上更为简单，但是由于Δα(x)均值水平接近于0而使得对普通最小二乘(OLS)准则下回归方程显著性水平偏低，作者建议先用式(7)形式计算α(x)，再由式(3)转换成Δα(x)。此时，回归直线段斜率即代表α(x)，通常与栅格数据模型的维数2比较接近，然后由式(4)计算Δα(x)。当li=1(带长度单位)时，可以获得拟合值即为c(x)值。这里，需要注意到c(x,l1=1)、c(x,l1=lmin)、c(x,l1=lcellsize)三者可以是完全不同的，c(x,l1=lcellsize)才表示原始栅格分辨率意义下对原始测量值的估值。 2.2.2ε=l/L=N(l)/N(L) 消除量纲 引入正规化尺度参数L(为一定值)，通过比值l/L消除ε的量纲。在栅格数据模型中，可以直接用窗口边长所占像元个数的比值：N(l)/N(L)等效表达物理长度比值l/L，这里，N(·)操作符表示取长度覆盖的栅格像元数，允许取小数表示等效个数。对非方形窗口，借鉴面积与边长的关系，N(l)的大小可由窗口覆盖像元数的平方根来等效表达。此时，Δα(x)值与c值均无量纲，Δα(x)值大小与不消除量纲情形获得的结果相同，而c值则不同。式(3)可转变为： 〈ρ(x,l)〉=c(x,l)(l/L)-Δα(x) (8) 对x位置处奇异性指数及其系数的估计可由log(l/L)-log(ρ(x,l)×(l/L)2)图中回归直线来确定，对应回归方程可表示为： log〈ρ(x,li)×(li/L)2)〉=log(c(x))+α(x)×log(li/L)，i=1,2,3，…,m (9) L参数加入到模型中，增加了模型处理获得c值的灵活性，以下介绍的迭代方法与L参数有紧密联系。L参数的取值可以是多样化的，下面列举4种典型的情形。 (1)L=lcellsize(栅格数据空间分辨率)，此时，ε的大小就是栅格像元个数，也即ε=l/L=N(l)。对栅格数据计算窗口边长通过数栅格像元个数在实际中操作简单，仅需对计算窗口的行列号进行操作，不涉及具体空间距离计算，因此在奇异性指数计算中使用频繁，此种情形c(x)值为原始数据的非线性拟合值。 (2)L=lmax(尺度窗口序列的最大值)，这种取值方式在许多分形计算场合常被采用。log(c(x))值在式(9)中为log-log图中截距(此时li=L)，相当于在lmax×lmax窗口滑动平均对数值的拟合值。在log-log拟合中，当L取lmax时，log(c(x))值在横轴的最大值位置获得，从回归分析误差角度考虑，由于偏离拟合数据的重心而有较大的误差。 (3)L取值为尺度窗口序列{li} (i=1,2,3，…,m)的几何平均值： 或表示成等效窗口数： (9) 此种情形克服了式(2)中的不足，log-log回归分析中平移纵轴使之经过数据重心，从而使得c(x)估值的误差最小化。用式(9)确定的L进行奇异性分析时，用c(x)值拟合相应窗口大小的滑动平均值具有最佳的预测可信度。 (4)L取值小于栅格数据的分辨率lcellsize。例如L=lcellsize/ 3，此时c(x)估计值表示了x位置处(lcellsize/ 3)空间分辨率意义下的分形插值结果。 综上可见，ε的量纲对于Δα值(或α值)的估算没有影响，而c值的大小及其量纲有差异，它的大小取决于正规化尺度参数L的取值方式。某些作者在开展局部奇异性分析模型时，对尺度指标ε的取值张冠李戴，不加区分：在实际处理中采用了像元个数N(l)来进行计算，却声称窗口计算参数为空间距离某某千米变化到某某千米。尽管对于α(x)的估值不受影响，但严格来说这是不妥当的，忽视了c(x)值的数量水平、量纲及其应有的意义。 从模型计算的通用性角度考虑，笔者建议用形式上更为一般化的式(8)作为基于栅格数据模型的常规局部奇异性的算法。可以在连续区间[lmin,lmax]上内插预测及适度外推预测任意窗口大小的滑动平均结果，或者获得更小分辨率的分形插值结果。本研究所介绍的迭代方法是在L>lcellsize的情形下作出推导的，此时c(x)值与原始数据具有相同的量纲，且c(x)值代表了L尺度上的光滑成分。 3.1 迭代方法的基本思想 局部奇异性分析方法的特点在于采用局部化的α估值与空间拓扑维数E值的差异性来实现异常识别，于是合理估计α值(或Δα值)就成为一个关键问题。利用一系列方形窗口基于栅格模型的局部奇异性分析算法简单，运算快捷。然而，嵌套的这一系列尺度窗口一次性计算所得的α值是最好的吗？最优的α值应该符合怎样的数学条件？笔者注意到与α值所伴随的c值的特性有其重要作用。从式(4)可见，局部奇异性分析可将原始数据分离出两大因子Δα成分和c成分。如果各个位置处的Δα值完全刻画了奇异性的强度，那么由全体c值构成的c集就应成为一个非奇异性的成分；否则，若c集仍带有奇异性，那么对应的Δα集便与理想的奇异性指数不对应，c集还可以继续分离出新的Δα成分和c成分。Δα成分和c成分应该分别代表奇异性因子和非奇异性因子，c成分应足够光滑，没有奇异性，具有良好的连续性和可微性。在实际计算中，获取Δα成分的log-log拟合过程中总会存在剩余误差，它并没有完全保证c成分不再带有奇异性信息。设法去除c集中的奇异性成分，使之具有足够高的光滑性，而将奇异性信息都归属于Δα集，这成为笔者研究局部奇异性分析迭代方法来优化α估值的关键思路(Chen et al, 2005, 2007)。 3.2 迭代方法的数学表达 记α*(x)和c*(x)分别为期望获取的一种最优的局部奇异性指数和局部系数，仍有 ρ(x)=c*(x)εα*(x)-E=c*(x)ε-Δa*(x) (10) 成立。这里，ε=l/L且L>lcellsize，盒子B(x,ε) 内的平均密度简记为ρ(x)，ρ(x)在所有空间位置的取值构成的数据集合称为ρ集，c*(x)在所有空间位置的取值构成的数据集合称为c*集。注意到，c*集与ρ集量纲相同，c*集反映了比ρ集更大尺度上的滑动平均非线性拟合信息，c*集比ρ集具有更小的离散程度。c*集可看作是ρ集的一个粗尺度逼近，其细节(粗糙度)信息保存在Δα*图中，理论上的c*集应不再具有多重分形特征，这样ρ集的奇异性信息都应集中反映在α*集中。若对c*集再做奇异性分析，则其计算所得的各个α值都趋近于同一个常数，即空间拓扑维数E。 基于同尺度的局部系数图迭代方法可以不断校正奇异性指数，将具有奇异性的c(x)通过迭代不断演化成非奇异的(Chen et al, 2007)：首先将ρ(x)分解成c0(x)和α0(x)，被估计出来的c0(x)可以进一步估计其局部奇异性，即把它当作新的ρ(x)，再次进行分解得到新的c1(x)，(下标1表示进行第1次迭代数，下同)，可知c1(x)的离散程度比c0(x)小，比ρ(x)更小，不断重复这样的过程，可以预料，最终获取的局部系数图将趋近于一种平稳状态，最终将不再具有奇异性，此时它的局部奇异性指数α趋近于自身的欧氏空间维数E，即Δα趋近于0。迭代方法与常规的非迭代方法的区别在于前者考虑了局部系数c(x)的作用。 令c-1(x) =ρ(x)，于是有如下迭代形式： ck-1(x)=ck(x)ε-Δαk(x),k=0,1,2,... (11) c集不断视为新的输入数据集，再次进行同样的局部奇异性计算，如此循环，新产生的c集将变得越来越光滑，趋近于一种平稳状态……于是最终可获得所期望的c*集。迭代分解的示意见图1。 图1 局部奇异性分析迭代方法示意图 假定进行了n次迭代时达到理论上的迭代终止条件：αn→0时，研究获取包括第0次(非迭代)和n次迭代的结果，有式(12)成立。 对比式(10)和式(12)，于是可得最终局部奇异性指数及系数的估值，见式(13)。 (13) 从中可见，当n=0时，迭代方法退化为常规方法；n>1时，迭代方法获取的最终奇异性指数α*的结果中则多了一个附加校正项，Δα*的表达则更为简练，是历次迭代所产生的Δαk的累加。可见，迭代方法具有对常规方法结果进行校正的能力，且其形式简洁。奇异性迭代方法得到了学者的关注，α*被称为最终奇异性指数(Agterberg, 2012)。 3.3 迭代终止条件与边缘效应 理论上的迭代终止条件要求Δαn中处处取值为0，在实际计算中这一条件过于严格，可用均方根误差(RMSE)作为判别条件，而均方根误差在数值上与Δαn的标准方差等价，于是，可直接用s(Δαk) < eps(eps为精度要求，例如0.01)迭代约束条件。 当精度要求过高时，将导致迭代过程无法在期望时间内结束，因此与一般的迭代算法一样，可人为限制迭代最多次数N0，使之定能在有限时间内完成计算。当迭代计数器k>N0时，跳出循环，终止迭代。作者对地质勘查数据的大量计算表明，经过数次迭代(通常3～5次)就能大幅度地降低c集的不光滑度，获得较理想的奇异性指数。当N0取值为0时，迭代方法就退化为常规方法。 当多尺度窗口滑动到栅格数据边缘位置时，一些窗口部分位置因覆盖于数据范围之外而造成数据获取的不完整，从而产生边缘效应，窗口尺度越大且迭代次数越多，边缘效应会越来越严重。对此问题，可以采用傅里叶变换类似的手段来扩边，补齐数据范围之外的像元取值；还可以采用缺失数据处理手段，仅对有效数据作邻域统计量计算，获得等效的平均密度。 3.4 基于栅格数据的迭代算法设计 局部奇异性分析的迭代算法可设计如下： 设置多尺度窗口参数{li},i=1,2,3，…,m 设置L参数且L>lcellsize 设置迭代精度eps与最多次数N0 输入数据初始化为原始栅格数据 最终的局部奇异性指数Δα*所有元素初始化为0 for(intk= 0,i<=N0;k++) { //约束迭代次数不超过N0 for(inti= 0,i for(int j=0;j 计算[i][j]栅格焦元位置m个不同 大小窗口的邻域平均值 //** 双对数log-log回归分析获得奇异性 指数αk[i][j]和系数ck[i][j] } } Δαk=E-αk; //默认二维情形E=2；退化为方向性奇异性，则E=1 Δα*+=Δαk; //最终局部奇异性指数的累计计算 if(s(Δαk)< eps) break; //迭代达到精度要求，结束迭代 else输入数据赋值更新为计算所得的ck集 } c*=ck; //获得最终的局部系数 注：**步骤可以采用积分图快速算法。 利用C++编程语言在Windows平台下开发了RasterDataMining 1.0软件，实现了基于栅格数据模型的局部奇异性分析专业分析功能。以主流文本栅格格式ArcGIS ASC Grid、Surfer 6 Text Grid来读写文件，制图功能则可在ArcGIS或Surfer平台下完成。该软件从栅格数据局部奇异性分析的深度应用需求出发，具有如下功能(图2)。 图2 局部奇异性程序功能界面(图上带圈数字编号说明内容见正文功能介绍) 具体功能按数字编号依次介绍如下。 (1) 实现局部奇异性常规方法和迭代方法，并支持缺失数据计算。 (2) 支持各向同性及各向异性邻域形状并自动构建多尺度窗口。 在栅格模型下直接以像元空间分辨率为单位来表达尺度大小，尺度大小为窗口边长所占的像元数，正规化尺度参数L也采用此方案。邻域形状可选方形、圆形、矩形、椭圆，有宽度和高度2个维度。窗口数、最内窗口、半径增量、旋转角度几项关键参数来创建多尺度邻域形状，其中，最内窗口设置值约定取奇数，默认为1，也可以设置成3、5等其他奇数值，此时的最小尺度大于像元的空间分辨率。多尺度窗口序列有2种自动构建方式：线性等间距(Linear space)，形成等差数列；对数等间距(Log space)，形成等比数列。在自动生成的基础上，用户可以自由修改尺度序列。 程序支持各向同性情形的奇异性分析，通过设置非等宽高的窗口参数(还可包括旋转角度)来实现。开发了某个像元是否在1个椭圆内(带旋转角度)的快速判断算法。 程序支持方向性奇异性分析计算。程序由二维自动退化为一维情形，即E=1。选择矩形窗口(可设置旋转角度)，固定其中一个维度的窗口大小不变化。例如，在高度方向窗口数设置为1或者3等定值(长度单位为栅格空间分辨率)且尺度半径增量固定为0，于是奇异性指数沿着东西方向进行计算。若高度窗口数固定为1，相当于原始栅格数据被分成多条东西向的一维数据序列；若高度窗口固定为3，则等效于原始栅格数据作南北向3个窗口的滑动平均后的多条一维数据序列。此时的奇异性指数图可突出南北向的差异性。 (3) 支持自定义计算位置，包括2个方面。图2中3a：自定义感兴趣区域进行计算，有3种选择：数据格网所有位置、数据格网有效数据为准、用户文件自定义位置；图2中3b：指定某一空间位置进行奇异性探查与预测(图3)。 图3 局部奇异性分析特定计算位置奇异性探查与预测 (4) 支持批量计算，包括2个方面。图2中4a：对多套不同的多尺度窗口进行批量计算；图2中4b：对多个输入栅格数据进行批量计算。 (5) 输出结果多样化，满足进一步分析的需要。结果文件保存不局限于奇异性指数，还包括系数c集，拟合优度等多种指标，可以追踪历次迭代的详细结果，也可以只保留最终结果。 此外，程序还提供了友好的栅格数据管理功能(支持ArcGIS和Surfer栅格数据批量相互转换)以及地图代数常见分析功能。 为帮助用户优化参数设置，程序提供了单个位置的计算功能，用户可以获得该位置的感兴趣信息：数据子集提取、邻域窗口统计、双对数回归分析、奇异性模型计算结果、因变量预测。图3显示WX数据集在行列坐标(53,75)位置处的奇异性计算信息，该位置计算所得Δα值为0.67(对应α值为1.33)，表明具有较强的正奇异性，拟合优度系数达到95.6%，可以通过显著性检验。按照所建立的奇异性模型，当窗口大小为lcellsize/2和lcellsize/3时，预测结果分别为22.67，29.70，明显高于原始栅格像元值20.41，与该位置的正奇异性特征一致，这提供了一种非线性插值方法。对应窗口大小为3.5个像元单位时，也可以通过奇异性公式获得预测值，为6.20，该预测值为分数维窗口内的滑动平均估值。 在前述参数设置条件下，增加迭代计算约束条件：标准误差s(Δαk)≤0.01，最多迭代次为25次，对整个数据集计算全图范围的局部奇异性(图4)。程序运行结果表明，当迭代到第20次时，s(Δαk)=0.009 7，于是终止迭代。图4展示了20次迭代过程中计算结果空间模式和常见统计量的变化趋势。图4曲线展示了s(Δαk)随迭代次数k的变化，可见由s(Δα0)到s(Δα1)具有最大的变化性，迭代4～5次后，s(Δαk)的下降速度越来越缓慢，趋近于0。与之相对应，局部系数c图也变得越来越光滑。原始数据是一个典型的多重分形数据集，具有强烈的不均匀性，对原始数据制图，仅能分辨有限的若干高值(显性异常，形如“大型山峰”)，而更多局部邻域上相对高值(弱缓异常，形如“小型山峰”)难以识别，将数据转换成对数后在一定程度上达到了异常增强效果，可以识别较多的局部异常。通过局部奇异性常规分析所获得的Δα0图可见，原始数据还有更多空间位置具有局部弱缓异常存在。通过迭代方法，可以获得任意迭代次数后的最终奇异性指数，由图4下面所展示的3个结果来看，异常的总体面貌雷同，但清晰程度各有差别。实际上，欲获得较好的奇异性指数，结合s(Δαk)曲线的变化性，对该数据集迭代4～5次就可得到较优的结果。 此外，对于批量处理等功能，限于篇幅，在此不一一展示。综上所述，所研发的基于栅格数据的局部奇异性分析程序计算效率高，功能全面，实用性强。核心算法代码还编译生成了动态链接库，提供给第三方使用。软件GeoDAS就采用了本程序的动态链接库，将局部奇异性计算与可视化结合到1个环境中，促进了方法的推广并方便了用户的使用。 图4 局部奇异性迭代分析结果空间模式和常见统计量趋势(枚举了迭代次数为0, 1,10,20的c图和Δα图，图件由Surfer软件生成，图件渲染颜色仅表示该图自身取值的相对高低，不同图件相同颜色值所对应的数据值无对比性，常见统计量见左统计表) 在深入在剖析局部奇异性分析基本思想及计算方法的基础上，引入正规化尺度参数L，提出了迭代方法来改进局部奇异性指数的估值。迭代方法使得在计算中通常被忽视的局部系数c值的作用被挖掘并利用。迭代分析方法的成立是建立在正规化尺度参数L大于栅格空间分辨率的基础上，通过迭代，可以获得越来越光滑的c图，与此相伴，在Δα中将积累更多的奇异性信息，因此迭代分析方法对Δα估值有其良好的特性。不同的L参数取值，可以获得不同的c值，增加了奇异性模型的灵活性。基于建立的奇异性模型，可以开展非线性空间插值、滑动平均估值。所设计的软件功能强劲，为用户应用迭代分析新方法提供了高效便捷的计算工具。不仅支持各向同性情形，还适用于各向异性尺度的奇异性计算和方向性奇异性计算。软件的动态链接库版本也在GeoDAS软件中被推广应用，促进了方法的推广。 本研究提出的迭代方法的数学公式同样适用于矢量数据模型。矢量数据作为最主要的空间数据类型之一，由于矢量数据空间位置的不规则性，相比栅格数据模型，需要优化算法提高时间效率，这将是后续研究的方向之一。目前研究的迭代策略是全图幅同步进行的，迭代使用的邻域形状也是全区一样，这些方面也可能进行局部化而优化迭代分析的结果。迭代方法中局部系数c集在迭代过程中间接地获取更大范围上的数据信息，具有加权滑动平均的效应。随着迭代次数的增加，各空间位置的ck(x)值被磨光的速率可能是不一致的，全图幅同步迭代将导致一些位置被过度计算。迭代次数越多边缘效应越严重，这在解释数据时需要引起注意。迭代过程中使用固定的邻域形状、固定的尺度变化范围，简化了计算，但忽视了局部系数在迭代过程中自身所发生的变化。本程序提供的批处理功能方便了用户从多种计算方案中优选结果，但这还不够智能化，更好的方案是从空间数据自身特点着手来发现最合适的窗口形状及尺度变化范围，例如借助空间U统计量方法(Cheng et al, 1996, 2007)。这些问题、难题的解决将有助于基于多重分形的局部奇异性方法的进一步发展与完善。 在写作过程中，笔者与Agterberg博士进行了有益的探讨，在此表示由衷的感谢！ 成秋明.2000.多重分形理论和地球化学因素分布规律[J].地球科学：中国地质大学学报,25(3):311-318. 成秋明.2004.空间模式的广义自相似性分析和矿产资源评价[J].地球科学：中国地质大学学报,29(6):733-744. 陈颙,陈凌.2005.分形几何学[M].2版.北京:地震出版社. 成秋明.2006.非线性成矿预测理论:多重分形奇异性-广义自相似性-分形谱系模型与方法[J].地球科学：中国地质大学学报,31(3):337-348. 成秋明.2007.成矿过程奇异性与矿产预测定量化的新理论与新方法[J].地学前缘,14(5):42-53. 陈志军.2007.多重分形局部奇异性分析方法及其在矿产资源信息提取中的应用[D].武汉:中国地质大学(武汉). 陈志军,成秋明,陈建国.2009.利用样本排序方法比较化探异常识别模型的效果[J].地球科学：中国地质大学学报,34(2):353-364. 成秋明.2012.覆盖区矿产综合预测思路与方法[J].地球科学：中国地质大学学报，37(6):1109-1125. AGTERBERG F P. 2001. Multifractal simulation of geochemical map patterns[J]. Journal of China University of Geosciences, 12(1): 31-39. AGTERBERG F P. 2003. Past and future of mathematical Geology[J]. Journal of China University of Geosciences, 14(3): 191-198. AGTERBERG F P. 2012. Sampling and analysis of chemical element concentration distribution in rock units and orebodies[J]. Nonlin Processes Geophys, 19: 23-44. CHENG QIUMING, AGGTERBERG F P, BONHAM-CARTER G F. 1996.Spatial analysis method for geochemical anomaly separation[J]. Journal of Geochemical Exploration, 56(3): 183-195. CHENG QIUMING. 1999. Multifractality and spatial statistics[J]. Computers & Geosciences, 25(9): 949-961. CHEN ZHIJUN, CHENG QIUMING, CHEN JIANGUO. 2005. Significance of fractal measure in local singularity analysis of multifractal model[C]//CENG QIUMING, GRAEME BONHAM-CATER. Proceedings ofIAMG'05: GIS and Spatial Analysis. International Association for Mathematical Geology. Wuhan, China: China University of Geosciences Printing House, (1): 475-480. CHENG QIUMING. 2006. GIS-based fractal/multifractal anomaly analysis for modeling and prediction of mineralization and mineral deposits[C]//HARRIS J R. GIS Application in the Earth Sciences. St John’s, NL, Canada:Geological Association of Canada, Special book: 289-300.CHEN ZHIJUN, CHENG QIUMING. 2007. Generalized local singularity analysis and its application[J]. Journal of China University of Geosciences, 18(Special issue): 348-350. CHEN ZHIJUN, CHENG QIUMING, XIE SHUYUAN, et al. 2007. A novel iterative approach for mapping local singularities from geochemical data[J]. Nonlin Processes Geophys, 14: 317-324. CHENG QIUMING, AGTERBERG F P. 2009. Singularity analysis of ore-mineral and toxic trace elements in stream sediments[J]. Computers & Geosciences, 35(2): 234-244.EVERTZ C J G, MANDEBROT B B. 1992. Multifractal measures: Appendix B[M]//PEITGEN H O, JURGENS H, SAUPE D. Chaos and Fractals. New York: Springer Verlag, 922-953. Iterative approach to local singularity analysis and software implementation based on raster data CHENZhi-jun1,2,CHENGQiu-ming1,2,3,CHENJian-guo1,2 (1. State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Wuhan 430074, China; 2. Faculty of Earth Resources, China University of Geosciences (Wuhan), Wuhan 430074, China; 3. Department of Earth and Space Science, Department of Geography, York University, Toronto M3J1P3, Canada) Weak anomaly identification was a challenging problem in the metallogeny prognosis and mineral resource assessment in recent years. Great attentions were drawn to the local singularity analysis based on multifractal theory for its new research ideas, simple method and beneficial effects to anomaly recognition. The authors introduced the singularity principles and its computational methods. The authors took into account the regularization scale parameter L and brought forward the iterative approach for the singularity analysis to improve the estimation of the singularity index. The iterative algorithm and software were realized and implemented using C++ language. This iterative approach could work not only for the isotropic case, but also could apply to the anisotropic scaling case and directional singularity index estimation. The dynamic-link library version of this software had been applied in GeoDAS which was professional GIS software for quantitative prediction of mineral resources. Local singularity analysis; Raster data; Multifractals; Iterative approach; Regularization scale parameter 10.3969/j.issn.1674-3636.2014.04.613 2014-08-18；编辑：陆李萍 国家自然科学基金项目(41272361, 40802081),中国地质调查局工作项目(12120113089300),中央高校基本科研业务费专项资金资助项目(CUG120102) 陈志军(1978— )，男，副教授，博士，主要从事数学地质的教学与科研工作, E-mail: chenzhijunCS@163.com P628 ：A ：1674-3636(2014)04-0613-10
2 基于栅格数据的算法剖析
3 迭代方法及算法设计
4 软件功能实现及其效果
5 结 论
6 致 谢
