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

1、多重分形技术与主成份分析,陈建国 张生元 中国地质大学（武汉）,一、多重分形及其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)代表

2、了局部地质事件，包括成矿过程等所产生的地球化学元素的富积，被称为异常值。,常用分解的方法有： (1)基于阈值的分解方法：设T0为给定的阈值，则A=T-T0，B=T0。这种方法的关键就是怎样确定阈值T0。 阈值的确定方法通常有： （a） 均值+n倍均方差 （b） 分形C-A法 （c） 证据权法 (2)趋势面插值和剩余异常分解：该方法将把T分解为趋势面B和剩余异常AT-B。 (3)频率滤波方法： 是应用傅立叶变换将地球化学场变换到频率域，在频率域构造适当的频率滤波器，对特定频率的信号进行过滤，在对滤波后的信号应用傅立叶逆变换将其变回到空间域的方法。 (4）地质统计学泛克立格法：漂移(背景）与涨落（

3、剩余异常） (5）分形滤波技术S-A：分形滤波技术是将频率滤波和CA方法相结合，通过在频率空间应用CA方法构造异常滤波器和背景滤波器的能谱密度滤波方法。 (6)衬度异常,2 分形理论,分形通常是指一种不规则的具有自相似特征的几何体。 1960 年，曼德尔布罗特在研究棉价变化的长期性态时，发现了价格在大小尺度间的对称性。在对尼罗河水位和英国海岸线的数学分析中，发现类似规律。他总结自然界中很多现象从标度变换角度表现出的对称性。他将这类集合称作自相似集。他认为，欧氏测度不能刻划这类集的本质，转向维数的研究，发现维数是尺度变换下的不变量，主张用维数来刻划这类集合。 1983年，曼德尔布罗特的新著自然界

4、的分形几何出版， 将分形定义为局部以某种方式与整体相似的集，重新讨论盒维数，正式将这一自相似性命名为分形(fractal)。,分形的定义,分形的实质是由所谓幂指数函数来控制的.对于分形几何体来说，对它们的某种度量M()和度量尺度之间的关系服从幂指数形式： M() E-D 这里表示“呈比例”，E、D和E-D分别表示拓扑维数(topological dimension)、分形维数(fractal dimension)和剩余分形维数(co-dimension)。 幂指数函数的特点之一是其具有尺度不变性，即改变度量尺度不影响幂指数函数的类型。该函数形式完全由分形维数所确定。这一性质就决定了分形具有自相

5、似性或统计自相似性。,设F是Rn中任一非空有界子集，记N（F，）表示覆盖F所需的直径为的闭盒子的最小数目， 则F的计盒维数D （分形维Fractal dimensiom）的计算公式为：,分形维数的确定:计盒维数（分形维）,3 多重分形理论,多重分形通常所描述的是定义在某一面积(二维)或体积(三维)中的一种测度。如果这种测度具有空间自相似性或统计自相似性，那么这种测度叫作多重分形。 通过这种测度的奇异性可将所定义的区域分解成这样一系列空间上镶嵌的子区域，每一子区域均构成单个分形，这样形成的分形除具有分形维数外，还具有各自度量的奇异性和一系列的分形维数。一些复杂的过程，如，重复叠加过程、湍流和布朗

6、运动已被认为服从多重分形规律。,多重分形模型分析技术,(1) 计盒方法 计盒方法是将图像的像素进行组合形成互不相交盒子的方法。例如， 0表示像素的尺度，表示盒子的尺度，若=3 0，则每个盒子由9个像素组成。由计盒方法所产生的盒子总数与盒子的尺度和空间的大小（L）成比例，即 其中，D是空间的拓扑维数。,(2)矩方法 定义多重分形系统的配分函数为 即将它定义为测度i()的q阶矩，q可以从到，N()表示集合S中i() 0的尺度为的盒子总数 其中，(q)称为质量指数 这样我们共引入了奇异性指数、多重分形谱和质量指数。作为描述同一物理对象的三个参数，它们之间的联系由统计物理中的勒让德变换来表示，有如下的

7、关系式：,配分函数与盒子尺度的loglog图,质量指数与指数q关系图,奇异指数分布图,多重分形谱分布图,几个重要的分形参数 计盒维数：当q0时， (0)为计盒维数； 信息维数：当q1时， (1)称为信息维数； 相关维数：当q2时， (2)称为相关维数； 多重分形度：,能谱空间上多重分形与A（S)的关系： 分形滤波技术,成秋明等（2000）提出了分形滤波（SA）技术。它在傅立叶能谱空间上度量物化探异常所对应的各向异性的广义自相似性，通过识别不同的广义自相似性并借助设计适当的分形滤波器将能谱密度进行滤波，进而利用傅立叶逆变换对物化探异常和背景进行分解。这样所圈定的物化探异常不仅具有形式的多样性（比

8、如，不同的异常强度、不同的范围、处于不同的背景等），而且它们在频率域具有与背景场表现不同的自相似性。这种自相似性可以由关系表示为 其中S是能谱密度，A(S)是在能谱密度空间上能谱密度大于S的面积，为分形模型的指数系数,对S与A(S)同时取对数，则显然S和A(S)之间存在线性关系。若将S与绘制在双对数图上，则S与A(S)所服从指数关系，在图上表现为线性关系。其斜率与指数系数相对应。 通常在ln-ln图上，S的不同取值区间对应于不同的线性关系。不同的直线段代表了不同的分形关系，两条直线的交点所对应的横坐标值（谱能量密度值）可被用来确定分形滤波器的阈值。,能谱密度 S与累积面积双对数图,直线段由最小

9、二乘法拟合,(1) 分形滤波器的构造 三种类型的分形滤波器可被构造：低通、高通和带通谱能量密度滤波器。 两条线段相交，取交点横坐标S0作为阈值，在S0的两边的两条线段具有不同的斜率，反映了满足不同的分形规律，定义两类滤波器 :,GA是一个高频低能谱密度滤波器，而则GB是低频高能谱量密度滤波器，通常GA被称为异常滤波器，GB被称为背景滤波器。,在log(S)-log(A(S)图上取某直线段的两端为S1和S2阈值，构造能谱带通滤波器,该滤波器将滤掉能谱密度小于S1或者大于S2的能谱成份，只保留在区间S1,S2的能谱成份，所以，可以GC看作是具有特定分形特征的分形滤波器。,应用Fourier 逆变换

10、，把在频率域滤波后的结果变回到空间域， 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)的双对数图。 通常，绘制能谱密度在整个范围的双对数图效果并不好，需要重新设定能谱密度分组。 能谱密度范围的设置可从菜单项“能谱密度分组”直接设置，也可以通过直接移动两端边界线，通过弹出菜单“重新计算分组频

11、数” 。,(2)直线拟合 由于S和A(S)之间满足分形关系，所以，在Ln-Ln图上表现为线性关系，这种关系包含在Ln(S)和Ln(A(S)散点图中，为了体现这种线性关系，有必要进行线性拟合。线性拟合可通过自动拟合和手动拟合来实现。 自动拟合：通过观察LnLn图确定Ln(S)和Ln(A(S)之间可能存在的直线段数，从菜单中设置直线段数，系统将按照散点图的斜率的变化自动拟合，自动拟合的效果未必是最好的，一部还需要通过观察，并不断手动调整段点的位置。 手动拟合：观察LnLn图，在散点图斜率明显改变的位置用鼠标手动插入分界点，也可以将不合适的点删除，或者将一个断点移动到合适的位置上，直线拟合由系统自动

12、完成。设置分形滤波器 当设置好断点后，系统将自动进行最小二乘拟合，并绘制拟合直线。,(3)分形滤波结果输出 首先要进行输出设置，默认设置是把滤波结果作为一个新变量加入到原数据文件中，也可以设置输出到指定的新文件中。 然后根据要提取信息的内容选择不同的分形滤波器。分形滤波器的选择是分形虑波技术的关键，一般选择多种滤波器进行比较是必要的。有三种滤波器可供选择，异常滤波器、背景滤波器和分形特征滤波器。 选择滤波器后，系统自动把分形滤波的结果保存到指定的文件中保存并显示。,用最小二乘法模拟的右边的两条线段的交点确立的阈值lnS0=9.05。SS0 代表背景(Cheng, Xu and Grunsky,

13、 2000)。通常，出于圈定和评价靶区的目的，我们只对异常感兴趣。通过逆Fourier 变换，并取lnS0=9.05作为阈值绘制Cu异常图（图2.11）。该图表明，大多数已知Cu矿床分布在Cu异常区，同时提供了一些新的异常区（找矿远景区）,分形滤波应用实例,Cu 原始含量图,分形滤波应用实例,Cu的LnA(S)-lnS图， 使用最小二乘法模拟三条具有不同斜率的线段，并分别获取三个临界点： LnS0=9.05, LnS1=7.3, and LnS2=6.5。,S-A法获取的Cu异常,取lnS0=9.05作为阈值绘制Cu异常图。该图表明，大多数已知Cu矿床分布在Cu异常区，同时提供了一些新的异常区

14、（找矿远景区）,二、主成份分析,主成份分析的基本原理 高阶主成份分析 加权主成份分析 主成份分析计算机实现,主成份分析的基本原理,主成份分析的目的是从多个原始变量中取若干线性组合，能尽可能多地保留原始变量中的信息。 设 X1，Xp是p个变量，考虑其线性变换。,显然：,i, j=1,2,p,如果要用Y1 尽可能多地保留原始变量X 的信息，经典的办法是使 的方差尽可能大，这需要对线性变换的系数a1加以限制，一般要求它是单位向量，即 l1l1=1。其它的各Yi也希望尽可能多地保留X 的信息，但前面的Y1,.,Yi-1 已保留的信息就不再保留，即要求cov(Yi,，Yj)=0 (j=1,i-1)，同时

15、要求lil=1，在这样的条件下使var（Yi )最大。设协方差阵 的特征值为 ，相应的单位特征向量分别为 a1,a2,ap (当特征根有重根时单位特征向量不唯一）。这时 的第i个主成分为:,记：A（a1,a2,ap)，Y=(Y1,Y2,Yp)，则A是正交矩阵，YAX，,(1)对样本数据的标准化 (2)计算相关矩阵 (3) 求特征值和特征向量 (4) 求主成分 (5)确定的主成份个数,主成份分析的计算过程,加权主成份分析方法,为了反映空间上各点数据在主成份分析中的不同作用，可以考虑加权主成份分析方法，设n为样本观察数，（w1,w2,wn)对于的权重，假设数据已经标准化，定义加权相关系数矩阵为：,

16、高阶相关系数,为了在主因子中充分体现高偏差值的作用，可以考虑高阶相关系数矩阵作为计算主因子的数据矩阵,主成份分析模块的实现,本模块是通过一个空间主成份分析向导逐步完成的，当给选定数据文件后，将在向导的引导下逐步完成空间主成份分析，具有非常好的用户界面和互动性。程序流程图如图1所示,输入MORPAS GRID数据文件,选择参加主成份分析的变量,设置计算模 式,确定主成份的个数,确定旋转模式,输出设置,结果输出,图1 空间主成份分析程序流程图,1.变量选择,模块将首先显示所有的变量，用户可根据需要选择进行主成份分析的变量。由于地球化学数据经常服从对数正态分布，所以对选择的变量，可以设置是否取对数，

17、如图2所示,图2 变量选择对话框,2.计算模式选择,进行主成份分析时，由于变量之间的量纲经常不一样，需要对数据进行标准化处理。数据是否标准化对应两种计算模式：相关系数模式和协方差模式。当选择相关系数模式时，就意味对数据进行标准化处理，这是本模块的默认模式。 原始数据的采集经常受到多种因素的影响，有些时候在某些位置上无法获得观察数据，这种情况在空间数据中表现为数据缺失，在进行主成份分析时需要确定这些缺失值的在空间数据中的表达方法以及对这些值缺失值的处理规则。如图3所示。在本模块中，缺失值的处理采取两种方法，一种是去掉该样品，或者用其平均值代替,图3 计算模式设置对话框,3.主成份个数的确定,在本

18、模块中主成份个数的确定对用户是透明的，在对话框中将显示因子陡坡图及和因子累积贡献图，如图4所示。确定主成份的个数可以通过观察该图按照下列方法确定。 由因子累积贡献率来确定，一般要求因子累积贡献率大于70%。 通过观察因子陡坡图，直接确定主成份的个数。 特征值大于1。,图4 确定主因子个数对话框,4.主轴旋转,因子载荷给出了观测变量和主成份之间相关程度的大小，这意味着在某一主成份上的载荷大的变量对该主成份的影响较大，因子的实际意义在很大程度上取决于这些变量。但是，实际中往往会出现所有变量在一个主成份上的负载都比较大的情形，这为因子的解释带来了困难。因子旋转（rotation of factors）为因子解释提供了便利。因子旋转的目的是使某些变量在某个主成份上的载荷较大，而在其它因子上的载荷则显著的低。在本模块中仅给出2种正交旋转方法，最大方差旋转和四次方最大旋转，如。图5所示。是否旋转或者选择什么方式旋转可通过观