有序样品的最优分割的算法及其在MATLAB中的实现

1 / 7有序样品的最优分割算法及其在Matlab中的实现一、 有序样品聚类一一最优分割的概念地质数据中，有些样品有一定的排列顺序，如沿地层剖面采集的岩石标本， 由钻孔取得的岩芯样品，由测井曲线所得的数据，由岩体中心到围岩的蚀变剖面 的样品等，它们是有序地质变量，在对这些有序样品进行分类时，不能打乱样品 的前后次序。所以，一些不考虑样品排列顺序的数学处理方法， 对此并不适用。 有序样品的聚类分析就是对有序样品进行分段的统计方法。对n个有序样品进行分割，就可能有2n-1种划分方法，这每一种分法成为一种分割。在所有的这些分 割中，有一种分割使得各段内部之间差异性最小，而短语段之间差异性最大。这 种对n个样品分段并使组内离差平方和最小的分割方法，成为最优分割法。这类问题的提法如下：设有一批（N个）按一定顺序排列的样品，每个样品测得 p项指标，其原始 资料矩阵：X (pXN)X11X12?X1NX21X22?X1N?Xp1Xp2?XpN其中元素X表示第j个样品的第i个指标的观测值。现在要把此 N个样品 按顺序（不破坏序列的连续性）进行分割（分段或者分类）。其所有可能的分割 法共有C1N-1+C2N-1 + C 3N-1 + +CN-1 N-1 = 2N-1-1种。现在要求在所有分割中找出一种分割法，这种分割法使得各段内样品之间的差异最小，而各段之间的差异最大。各段内部差异最小，即各段内数值变化最小，段内数值变化可用变差或者极 差来表示，比如样品段 Xi、Xi+1、Xi+2、Xj :变差：dij=X i,ja=1dij表示样品段 Xi、Xi+八 Xi+2、Xj内样品间的差异情况，dij小表示段内 各样品之间数值比较接近，反之，dj大表示段内各样品数值之间的差异大。极差：dj二 =1 (maxi<3<jxa3- mini勻倍)对于单指标情况，则dij= (maxi <3<jx 3 - mg <3<jx 3)要各段内部的差异最小，即所分成各段变差的总和（即段内离差平方和，称 为总变差）为最小。总变差分解公式：S总=S段间+S段内S总=-m|=1j=1 (Xj|-x)2=m:|=1nij=1 (X j| -X|) + (X| -X)2=m:|=1j=1(Xj| -x 2 ,mX|)+|=1nl j=1(X| -X)2 + 2 mij=1 (Xij-X|)(X |-X)= m |=1j= (Xj| -2 .mX|)+|=1n| j=1(X| -x)2=S 段内 +S 段间其中：D = 2 m=i 善(Xij - X|)(X| - x)=2 mi (Xij - X|) j= (X| - x)=2 mi X| - x (n|X| - n|X|)=0s段内=m1善(Xj| - X|)2为段内离差平方和s段间=m=i善(X|-x)2为段间离差平方和所以S段间=S总-S段内对给定的N个样品，S总是个固定的量。若使段内离差平方和为最小，则段间 离差平方和必为最大。所以，使总变差（段内离差平方和）为最小的分割法就是 最优的分割法。二、有序样品聚类的最优分割意义最优分割在地质研究中是一个非常有用的手段，只要地质体的某些地球化学特征存在规律性的差异，采用最优分割的数学树立方法，就能按顺序在最理想的 地方进行分段。通过对地层中采集的具某些地球化学特征样品的最优分割，能在地层的划分对比中发挥重要的辅助作用； 在找矿过程中，该方法更显得天独厚的 优势，它能进行蚀变、矿化及矿体的准确分带；对岩浆岩相带划分及演化序列的 研究也十分有效。第四纪地层岩相和厚度变化很大，给第四纪地层划分和对比地 层带来了困难，但第四系地层是不同阶段和不同的环境条件下形成的， 如重矿物、 微量元素等地球化学特征的分布规律与一定的沉积阶段和沉积环境相对应，因此,最优分割法的数据处理，可能是第四纪地层划分和对比的有效方法。最优分割的计算步骤及其计算公式1. 数据正规化 设原始资料矩阵为X (pXN)将矩阵X中的元素Xj变换为:Xj - min 1 弓qx jZij =(i=1、max i 討 qXj- min 1 弓qXj而得矩阵X11X12?X1NX21X22?X1N? ?Xp1Xp2?XpN2、p；j=1、2、N)Z(p >N)=乙jMatlab 代码：function std= std1(vector)%对矩阵进行标准化%vector为待分割矩阵max1=max(vector);% 对歹U求最大值min1=min( vector);%对列求最小值a,b=size(vector); % 矩阵大小,a为行数,b为列数 for j=1:bstd(j)= (vector(j)-mi n1)/(max1-mi n1);end2. 计算极差(或变差)矩阵由上述极差(或变差)计算公式得到矩阵：d12 d13 ?d1Nd23d2ND =?dN- 1NMatlab 代码：fun ctio n D,a,b=ra nge1(vector)%D返回计算所得的极差矩阵a,b=size(vector); % 求矩阵大小,a为行数,b为列数 k=a;%当只计算单指标数据时，k=a=1for i=1:bfor j=i:bd(i,j)=max(vector(k,i:j)-mi n(vector(k,i:j);endendD=d;3. 进行最优二分割由D矩阵计算全部分两类的各种分割相应的总变差，即对每一个m ( m = N、N-1、2),求出相应的总变差Sm (2; j) (j = 1、2、m-1 )找出最小值，确定各子段的最优二分割点a 1 (m,即Sm (2; a (m) = mini弓奇-iSm(2;j)从而得出N个样品的最优二分割 Xi、x2、Xai(N) X al(N)+1、XNMatlab 代码：function S,alp=divi2(vector,n)%最优二分隔,S为最优二分割各段的分割点，%a1记录了二分割的序号d,a,b=ra nge1(vector);alp=ones(n-1,b);%al(i,j)表示前i个样品的第j次分割点S=zeros(b,b);for m=2:bfor j=1:m-1s(m,j)=d(1,j)+d(j+1,m);endS_temp(m,1)=mi n(s(m,1:m-1);for j=1:m-1if S_temp(m,1)=s(m,j);alp( n-1,m)=j;endend for t=1:mS(t,alp( n-1,t)=S_temp(t,1);endend4. 进行最优三分割对于 m=N N-1、4、3,由 S( a ; a )(j=2、3、4、m-1 )及 D 矩阵分别计算：Sm(3; a1(j),j)=S j(2; a1(j)+d j+1,N(m=N、N-1、4、3)然后求出最小值，即Sm (3; a (m), a (m)= min?討奇-1Sm(3; a j ,j)从而得到N个样品最优三分割 X1、X2、X a1( N) X al( N)+1、X a2 ( N) X a2 (N)+!、XN 5. 最优K分割完全类似4的办法，在最优三分割的基础上可以进行最优四分割，继而进行 五分割，以此类推，如果已经作出最优K-1分割，则可以产生最优 K分害叽Matlab 代码：function S,alp=divi(vector,n)%n为要分割的段数d,a,b=ra nge1(vector);alp=zeros(1,b);%al(i,j)表示前i个样品的第j次分割点for m=n:bfor j=n-1:m-1if n=2s(m,j)=d(1,j)+d(j+1,m);elseS,alp=divi(vector, n-1); s(m,j)=S(j,alp(j)+d(j+1,m);endendS=zeros(b,b);S_temp(m,1)=mi n( s(m, n-1:m-1);for j=1:m-1if S_temp(m,1)=s(m,j);alp(m)=j;endendfor t=1:mif alp(t)=0S(t,alp(t)=S_temp(t,1);endendendfunction array = sect(vector, n)%vector为样品矩阵，当直接对样品矩阵进行分割时调用该函数%array返回样品最优n分割的分割点号%n为要分割的段数a,b=size(vector);for num=n :-1:2S,alp=divi(vector, nu m);if num = narray( nu m-1)=alp(1,b);elsearray( nu m-1)=alp(array( nu m);endendfunction array=fsect(file name ,n)%filename为需要分割的样品数据的文件名，如d:temp.txt '%要对数据文件进行分割时调用该函数fid=fope n(filen ame,'r');A,cou nt=fsca nf(fid,'%f);vector=A'a,b=size(vector);for num=n :-1:2S,alp=divi(vector, nu m);if num = narray( nu m-1)=alp(1,b);elsearray( nu m-1)=alp(array( nu m);endend四、实例分析琼州海峡标志性钻孔古地磁样品的磁化率参数，是海底沉积物在地质历史上 各时期沉积产物受磁化作用强弱的重要参数，与沉积环境的磁场变化、氧化还原条件等因素有关，可据此对地层沉积的阶段性进行划分 以下是钻孔0-10m岩芯中取得的12个样品数据:ID磁化率ID磁化率ID磁化率16.055.798.326.066.3107.735.375.3117.744.084.71210.3则样品数据矩阵为：Vector=6.0 6.0 5.3 4.0 5.7 6.3 5.3 4.7 8.3 7.7 7.7 10.31、样品正规化>>v = stdl(Vector)v =0.3175 0.3175 0.20630 0.2698 0.3651 0.20630.1111 0.6825 0.5873 0.58731.00002、计算极差矩阵>>D,a,b=ra nge1(v);>>DD=U0CULL】D.31T5D.aiTBQ.3E&1a. 3E6Jou&h山 682 &0. E821叽 5825LOOOO00CL 3 LT5D. 31T60.3EB10i 30510U3S610bd&250. 69250, i6251, 00000000.20(13D.2E9A0.1 DEI0.3BE1OL 3051Dl B8550. 682EIL 0825LOOOO0(100DlSBSS0.3dElDL3GE1Dl ees E(L BBSS(L S8'25LOOOO01.1u000- 09520,168?(LOu &TH0<B7140, &J14山 996900000Qo. issrOl霑川0. 57140. 57H0. 5714o. ssa»00000-a0UD&D20. S7140. 57140. 5714o. eseaUUuDU00aD_ &?14U. £714U.射MOL HUH 9IJn0000QiD.贬0. roaz0,4J2T00000fla00D00-41270000a:aaDDD0.412?00a0aflaaD00D3、最优K分割设K=4>>array=sect(v,4)array =7811即样品段 Vector=6.0 6.0 5.3 4.0 5.7 6.3 5.3 4.7 8.3 7.7 7.7 10.3的最优四分割的分割点为8和11,即最优四分割方案为：6.0 6.0 5.3 4.0 5.7 6.3 5.3 4.7 8.3 7.7 7.7 10.3 五、 结论本文首先对有序样品的最优分割发原理及算法进行了深入分析，在深入理解原理的基础上，设计了 Matlab算法并独立完成了代码编写，使得这一较为复杂 的数学计算问题能够在计算机程序中快速的计算，以得到满意的答案。但同时，算法设计与程序代码中仍存在一些问题，现列如下：1、用以衡量段间和段内样品数据差异的指标因不同的应用环境而应采用不 同的指标，本文中作为示例仅以极差作为衡量指标。2、由于程序中使用了递归算法，使得当分割段数较多的情况下计算量迅速 增大，计算时间较长，算法有待改进。3、未能及时进行视窗设计，实现软件化，所以实用代码仅能在Matlab命令 窗里调用计算。7 / 7