多重分形技术与主成份分析

多重分形技术与主成份分析陈建国张生元中国地质大学（武汉）[e-mail]一、多重分形及其S-A异常分解技术分形理论多重分形理论多重分形滤波技术及其软件实现1.地球化学场的分解观察的地球化学场T(x,y)一般被认为是各种大尺度的地质过程和小尺度的成矿过程形成的，这两种不同的过程可能导致在同一区域上的同一元素的富集具有截然不同的分布，分别用B(x,y)和A(x,y)表示。T(x,y)=B(x,y)+A(x,y)其中，T(x,y)表示地球化学元素在位置（x,y）的观察值；B(x,y)代表了由区域地质事件所引起的部分，通常被称为背景值；A(x,y)代表了局部地质事件，包括成矿过程等所产生的地球化学元素的富积，被称为异常值。常用分解的方法有：(1)基于阈值的分解方法：设T0为给定的阈值，则A=T-T0，B=T0。这种方法的关键就是怎样确定阈值T0。阈值的确定方法通常有：（a）均值+n倍均方差（b）分形C-A法（c）证据权法(2)趋势面插值和剩余异常分解：该方法将把T分解为趋势面B和剩余异常A＝T-B。(3)频率滤波方法：是应用傅立叶变换将地球化学场变换到频率域，在频率域构造适当的频率滤波器，对特定频率的信号进行过滤，在对滤波后的信号应用傅立叶逆变换将其变回到空间域的方法。(4）地质统计学泛克立格法：漂移(背景）与涨落（剩余异常）(5）分形滤波技术S-A：分形滤波技术是将频率滤波和C－A方法相结合，通过在频率空间应用C－A方法构造异常滤波器和背景滤波器的能谱密度滤波方法。(6)衬度异常2分形理论分形通常是指一种不规则的具有自相似特征的几何体。1960年，曼德尔布罗特在研究棉价变化的长期性态时，发现了价格在大小尺度间的对称性。在对尼罗河水位和英国海岸线的数学分析中，发现类似规律。他总结自然界中很多现象从标度变换角度表现出的对称性。他将这类集合称作自相似集。他认为，欧氏测度不能刻划这类集的本质，转向维数的研究，发现维数是尺度变换下的不变量，主张用维数来刻划这类集合。1983年，曼德尔布罗特的新著《自然界的分形几何》出版，将分形定义为局部以某种方式与整体相似的集，重新讨论盒维数，正式将这一自相似性命名为分形(fractal)。分形的定义分形的实质是由所谓幂指数函数来控制的.对于分形几何体来说，对它们的某种度量M(δ)和度量尺度之间的关系服从幂指数形式：M(δ)∝δE-D这里∝表示“呈比例”，E、D和E-D分别表示拓扑维数(topologicaldimension)、分形维数(fractaldimension)和剩余分形维数(co-dimension)。幂指数函数的特点之一是其具有尺度不变性，即改变度量尺度δ不影响幂指数函数的类型。该函数形式完全由分形维数所确定。这一性质就决定了分形具有自相似性或统计自相似性。设F是Rn中任一非空有界子集，记N（F，δ）表示覆盖F所需的直径为δ的闭盒子的最小数目，则F的计盒维数D（分形维Fractaldimensiom）的计算公式为：分形维数的确定:计盒维数（分形维）)/1ln()),(ln(limdim0FNFDB3多重分形理论多重分形通常所描述的是定义在某一面积(二维)或体积(三维)中的一种测度。如果这种测度具有空间自相似性或统计自相似性，那么这种测度叫作多重分形。通过这种测度的奇异性可将所定义的区域分解成这样一系列空间上镶嵌的子区域，每一子区域均构成单个分形，这样形成的分形除具有分形维数外，还具有各自度量的奇异性和一系列的分形维数。一些复杂的过程，如，重复叠加过程、湍流和布朗运动已被认为服从多重分形规律。多重分形模型分析技术(1)计盒方法计盒方法是将图像的像素进行组合形成互不相交盒子的方法。例如，ε0表示像素的尺度，ε表示盒子的尺度，若ε=3ε0，则每个盒子由9个像素组成。由计盒方法所产生的盒子总数与盒子的尺度和空间的大小（L）成比例，即其中，D是空间的拓扑维数。DLN)/()((2)矩方法定义多重分形系统的配分函数为即将它定义为测度μi(ε)的q阶矩，q可以从－∞到＋∞，N(ε)表示集合S中μi(ε)≠0的尺度为ε的盒子总数其中，τ(q)称为质量指数这样我们共引入了奇异性指数、多重分形谱和质量指数。作为描述同一物理对象的三个参数，它们之间的联系由统计物理中的勒让德变换来表示，有如下的关系式：)(1)()(Niqiq)log(loglim)log()(loglim)()(100Niqiqq)()()()()()()(qdqqdqqqqfdqqdq配分函数与盒子尺度的log－log图质量指数与指数q关系图奇异指数分布图多重分形谱分布图几个重要的分形参数计盒维数：当q＝0时，τ(0)为计盒维数；信息维数：当q＝1时，τ(1)称为信息维数；相关维数：当q＝2时，τ(2)称为相关维数；多重分形度：)0()1(2)2()1(''能谱空间上多重分形与A（S)的关系：分形滤波技术成秋明等（2000）提出了分形滤波（S－A）技术。它在傅立叶能谱空间上度量物化探异常所对应的各向异性的广义自相似性，通过识别不同的广义自相似性并借助设计适当的分形滤波器将能谱密度进行滤波，进而利用傅立叶逆变换对物化探异常和背景进行分解。这样所圈定的物化探异常不仅具有形式的多样性（比如，不同的异常强度、不同的范围、处于不同的背景等），而且它们在频率域具有与背景场表现不同的自相似性。这种自相似性可以由关系表示为其中S是能谱密度，A(S)是在能谱密度空间上能谱密度大于S的面积，β为分形模型的指数系数SSA)(对S与A(S)同时取对数，则显然S和A(S)之间存在线性关系。若将S与绘制在双对数图上，则S与A(S)所服从指数关系，在图上表现为线性关系。其斜率与指数系数相对应。通常在ln-ln图上，S的不同取值区间对应于不同的线性关系。不同的直线段代表了不同的分形关系，两条直线的交点所对应的横坐标值（谱能量密度值）可被用来确定分形滤波器的阈值。能谱密度S与累积面积双对数图,直线段由最小二乘法拟合(1)分形滤波器的构造三种类型的分形滤波器可被构造：低通、高通和带通谱能量密度滤波器。两条线段相交，取交点横坐标S0作为阈值，在S0的两边的两条线段具有不同的斜率，反映了满足不同的分形规律，定义两类滤波器:00S)S(0S)S(1)(AG00S)S(0S)S(1)(BGGA是一个高频低能谱密度滤波器，而则GB是低频高能谱量密度滤波器，通常GA被称为异常滤波器，GB被称为背景滤波器。在log(S)-log(A(S))图上取某直线段的两端为S1和S2阈值，构造能谱带通滤波器0S)S(S1)(21其它CG该滤波器将滤掉能谱密度小于S1或者大于S2的能谱成份，只保留在区间〔S1,S2〕的能谱成份，所以，可以GC看作是具有特定分形特征的分形滤波器。应用Fourier逆变换，把在频率域滤波后的结果变回到空间域，B=F-1(F(T)GB)A=F-1(F(T)GA)C=F-1(F(T)GC)这里F和F-1分别表示Fourier变换和Fourier逆变换分形滤波程序流程图多重分形模块菜单和弹出菜单(1)S与A(S)双对数散点图及其设置首先对数据进行傅立叶变化，得到其能谱密度，对能谱密度进行分组，并计算在各组中的频数和累积频数，绘制S与A(S)的双对数图。通常，绘制能谱密度在整个范围的双对数图效果并不好，需要重新设定能谱密度分组。能谱密度范围的设置可从菜单项“能谱密度分组”直接设置，也可以通过直接移动两端边界线，通过弹出菜单“重新计算分组频数”。(2)直线拟合由于S和A(S)之间满足分形关系，所以，在Ln-Ln图上表现为线性关系，这种关系包含在Ln(S)和Ln(A(S))散点图中，为了体现这种线性关系，有必要进行线性拟合。线性拟合可通过自动拟合和手动拟合来实现。自动拟合：通过观察Ln－Ln图确定Ln(S)和Ln(A(S))之间可能存在的直线段数，从菜单中设置直线段数，系统将按照散点图的斜率的变化自动拟合，自动拟合的效果未必是最好的，一部还需要通过观察，并不断手动调整段点的位置。手动拟合：观察Ln－Ln图，在散点图斜率明显改变的位置用鼠标手动插入分界点，也可以将不合适的点删除，或者将一个断点移动到合适的位置上，直线拟合由系统自动完成。设置分形滤波器当设置好断点后，系统将自动进行最小二乘拟合，并绘制拟合直线。(3)分形滤波结果输出首先要进行输出设置，默认设置是把滤波结果作为一个新变量加入到原数据文件中，也可以设置输出到指定的新文件中。然后根据要提取信息的内容选择不同的分形滤波器。分形滤波器的选择是分形虑波技术的关键，一般选择多种滤波器进行比较是必要的。有三种滤波器可供选择，异常滤波器、背景滤波器和分形特征滤波器。选择滤波器后，系统自动把分形滤波的结果保存到指定的文件中保存并显示。用最小二乘法模拟的右边的两条线段的交点确立的阈值lnS0=9.05。SS0通常代表异常，SS0代表背景(Cheng,XuandGrunsky,2000)。通常，出于圈定和评价靶区的目的，我们只对异常感兴趣。通过逆Fourier变换，并取lnS0=9.05作为阈值绘制Cu异常图（图2.11）。该图表明，大多数已知Cu矿床分布在Cu异常区，同时提供了一些新的异常区（找矿远景区）分形滤波应用实例Cu原始含量图分形滤波应用实例Cu的LnA(S)-lnS图，使用最小二乘法模拟三条具有不同斜率的线段，并分别获取三个临界点：LnS0=9.05,LnS1=7.3,andLnS2=6.5。S-A法获取的Cu异常取lnS0=9.05作为阈值绘制Cu异常图。该图表明，大多数已知Cu矿床分布在Cu异常区，同时提供了一些新的异常区（找矿远景区）二、主成份分析1.主成份分析的基本原理2.高阶主成份分析3.加权主成份分析4.主成份分析计算机实现主成份分析的基本原理主成份分析的目的是从多个原始变量中取若干线性组合，能尽可能多地保留原始变量中的信息。设X1，…，Xp是p个变量，考虑其线性变换。pppppppppppXlXlXlXlYXlXlXlXlYXlXlXlXlY2211'2222112'221221111'11显然：jijiiiillYYllYVar''),cov()(i,j=1,2,…,p如果要用Y1尽可能多地保留原始变量X的信息，经典的办法是使的方差尽可能大，这需要对线性变换的系数a1加以限制，一般要求它是单位向量，即l1’l1=1。其它的各Yi也希望尽可能多地保留X的信息，但前面的Y1,..,Yi-1已保留的信息就不再保留，即要求cov(Yi,，Yj)=0(j=1,…,i-1)，同时要求li’l=1，在这样的条件下使var（Yi)最大。设协方差阵的特征值为，相应的单位特征向量分别为a1,a2,…,ap(当特征根有重根时单位特征向量不唯一）。这时的第i个主成分为:021piiii＝λYXaY)var(p,1,2,i'且记：A＝（a1,a2,…,ap)，Y=(Y1,Y2,…,Yp)，则A是正交矩阵，Y＝AX，的主对角元素是其中ii11,pipiiii被称为因子载荷，（的相关系数是与变量主成份),(,1),(,),(p,1,2,ki,/)Y:XY121iikikikpkikpikikiiikiiXYXYXYaX(1)对样本数据的标准化(2)计算相关矩阵(3)求特征值和特征向量(4)求主成分(5)确定的主成份个数主成份分析的计算过程加权主成份分析方法为了反映空间上各点数据在主成份分析中的不同作用，可以考虑加权主成份分析方法，设n为样本观察数，（w1,w2,…,wn)对于的权重，假设数据已经标准化，定义加权相关系数矩阵为：pjiwxxwr