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

1/7有序样品的最优分割算法及其在Matlab中的实现一、有序样品聚类——最优分割的概念地质数据中，有些样品有一定的排列顺序，如沿地层剖面采集的岩石标本，由钻孔取得的岩芯样品，由测井曲线所得的数据，由岩体中心到围岩的蚀变剖面的样品等，它们是有序地质变量，在对这些有序样品进行分类时，不能打乱样品的前后次序。所以，一些不考虑样品排列顺序的数学处理方法，对此并不适用。有序样品的聚类分析就是对有序样品进行分段的统计方法。对n个有序样品进行分割，就可能有2n-1种划分方法，这每一种分法成为一种分割。在所有的这些分割中，有一种分割使得各段内部之间差异性最小，而短语段之间差异性最大。这种对n个样品分段并使组内离差平方和最小的分割方法，成为最优分割法。这类问题的提法如下：设有一批（N个）按一定顺序排列的样品，每个样品测得p项指标，其原始资料矩阵:X（p×N）=x11x12⋯x1Nx21x22⋯x1N⋮⋮⋱⋮xp1xp2⋯xpN其中元素xij表示第j个样品的第i个指标的观测值。现在要把此N个样品按顺序（不破坏序列的连续性）进行分割（分段或者分类）。其所有可能的分割法共有C1N-1+C2N-1+C3N-1+…+CN-1N-1=2N-1-1种。现在要求在所有分割中找出一种分割法，这种分割法使得各段内样品之间的差异最小，而各段之间的差异最大。各段内部差异最小，即各段内数值变化最小，段内数值变化可用变差或者极差来表示，比如样品段｛xi、xi+1、xi+2、…、xj｝:变差：dij=[xα−xjα=i(i,j)]2xi,j=1j−i+1xαjα=1dij表示样品段｛xi、xi+1、xi+2、…、xj｝内样品间的差异情况，dij小表示段内各样品之间数值比较接近，反之，dij大表示段内各样品数值之间的差异大。极差：dij=(maxi≤β≤jxαβ−mini≤β≤jxαβ)pα=1对于单指标情况，则2/7dij=（maxi≤β≤jxβ−mini≤β≤jxβ）要各段内部的差异最小，即所分成各段变差的总和（即段内离差平方和，称为总变差）为最小。总变差分解公式：S总=S段间+S段内S总=ml=1(xjlnlj=1−x)2=ml=1[(xjlnlj=1−xl)+(xl−x)]2=ml=1(xjlnlj=1−xl)2+ml=1(xl−x)2nlj=1+2ml=1(xij−xl)(xl−x)nlj=1=ml=1(xjlnlj=1−xl)2+ml=1(xl−x)2nlj=1=S段内+S段间其中：D=2ml=1(xij−xl)(xl−x)nlj=1=2(xij−xl)ml=1(xl−x)nlj=1=2xl−x(nlxl−nlxl)ml=1=0S段内=ml=1(xjlnlj=1−xl)2为段内离差平方和S段间=ml=1(xl−x)2nlj=1为段间离差平方和所以S段间=S总-S段内对给定的N个样品，S总是个固定的量。若使段内离差平方和为最小，则段间离差平方和必为最大。所以，使总变差（段内离差平方和）为最小的分割法就是最优的分割法。二、有序样品聚类的最优分割意义最优分割在地质研究中是一个非常有用的手段，只要地质体的某些地球化学特征存在规律性的差异，采用最优分割的数学树立方法，就能按顺序在最理想的地方进行分段。通过对地层中采集的具某些地球化学特征样品的最优分割，能在地层的划分对比中发挥重要的辅助作用；在找矿过程中，该方法更显得天独厚的优势，它能进行蚀变、矿化及矿体的准确分带；对岩浆岩相带划分及演化序列的研究也十分有效。第四纪地层岩相和厚度变化很大，给第四纪地层划分和对比地层带来了困难，但第四系地层是不同阶段和不同的环境条件下形成的，如重矿物、微量元素等地球化学特征的分布规律与一定的沉积阶段和沉积环境相对应，因此，3/7最优分割法的数据处理，可能是第四纪地层划分和对比的有效方法。三、最优分割的计算步骤及其计算公式1.数据正规化设原始资料矩阵为X（p×N）=x11x12⋯x1Nx21x22⋯x1N⋮⋮⋱⋮xp1xp2⋯xpN将矩阵X中的元素xij变换为：zij=xij−min1≤j≤n{xij}max1≤j≤n{xij}−min1≤j≤n{xij}(i=1、2、…、p；j=1、2、…、N)而得矩阵Z(p×N)=[zij]Matlab代码：function[std]=std1(vector)%对矩阵进行标准化%vector为待分割矩阵max1=max(vector);%对列求最大值min1=min(vector);%对列求最小值[a,b]=size(vector);%矩阵大小,a为行数,b为列数forj=1:bstd(j)=(vector(j)-min1)/(max1-min1);end2.计算极差（或变差）矩阵由上述极差（或变差）计算公式得到矩阵：D=d12d13⋯d1Nd23…d2N⋱⋮dN−1NMatlab代码：function[D,a,b]=range1(vector)%D返回计算所得的极差矩阵[a,b]=size(vector);%求矩阵大小,a为行数,b为列数k=a;%当只计算单指标数据时，k=a=1fori=1:bforj=i:bd(i,j)=max(vector(k,i:j))-min(vector(k,i:j));endendD=d;3.进行最优二分割4/7由D矩阵计算全部分两类的各种分割相应的总变差，即对每一个m（m=N、N-1、…、2），求出相应的总变差Sm（2；j）（j=1、2、…、m-1）找出最小值，确定各子段的最优二分割点α1（m），即Sm（2；α1（m））=min1≤j≤m−1Sm(2;j)从而得出N个样品的最优二分割｛x1、x2、…、xα1（N）｝{xα1（N）+1、…、xN}Matlab代码：function[S,alp]=divi2(vector,n)%最优二分隔,S为最优二分割各段的分割点，%a1记录了二分割的序号[d,a,b]=range1(vector);alp=ones(n-1,b);%al(i,j)表示前i个样品的第j次分割点S=zeros(b,b);form=2:bforj=1:m-1s(m,j)=d(1,j)+d(j+1,m);endS_temp(m,1)=min(s(m,1:m-1));forj=1:m-1ifS_temp(m,1)==s(m,j);alp(n-1,m)=j;endendfort=1:mS(t,alp(n-1,t))=S_temp(t,1);endend4.进行最优三分割对于m=N、N-1、…、4、3，由Sj(α;α1(j))（j=2、3、4、…、m-1）及D矩阵分别计算：Sm(3;α1(j),j)=Sj(2;α1(j))+dj+1,N(m=N、N-1、…、4、3)然后求出最小值，即Sm（3；α1（m），α2（m））=min2≤j≤m−1Sm(3;α1j,j)从而得到N个样品最优三分割｛x1、x2、…、xα1（N）｝｛xα1（N）+1、…、xα2（N）｝｛xα2（N）+！、…、xN｝5.最优K分割5/7完全类似4的办法，在最优三分割的基础上可以进行最优四分割，继而进行五分割，以此类推，如果已经作出最优K-1分割，则可以产生最优K分割。Matlab代码：function[S,alp]=divi(vector,n)%n为要分割的段数[d,a,b]=range1(vector);alp=zeros(1,b);%al(i,j)表示前i个样品的第j次分割点form=n:bforj=n-1:m-1ifn==2s(m,j)=d(1,j)+d(j+1,m);else[S,alp]=divi(vector,n-1);s(m,j)=S(j,alp(j))+d(j+1,m);endendS=zeros(b,b);S_temp(m,1)=min(s(m,n-1:m-1));forj=1:m-1ifS_temp(m,1)==s(m,j);alp(m)=j;endendfort=1:mifalp(t)~=0S(t,alp(t))=S_temp(t,1);endendendfunction[array]=sect(vector,n)%vector为样品矩阵，当直接对样品矩阵进行分割时调用该函数%array返回样品最优n分割的分割点号%n为要分割的段数[a,b]=size(vector);fornum=n:-1:2[S,alp]=divi(vector,num);ifnum==narray(num-1)=alp(1,b);elsearray(num-1)=alp(array(num));endend6/7function[array]=fsect(filename,n)%filename为需要分割的样品数据的文件名，如’d:\temp.txt’%要对数据文件进行分割时调用该函数fid=fopen(filename,'r');[A,count]=fscanf(fid,'%f');vector=A';[a,b]=size(vector);fornum=n:-1:2[S,alp]=divi(vector,num);ifnum==narray(num-1)=alp(1,b);elsearray(num-1)=alp(array(num));endend四、实例分析琼州海峡标志性钻孔古地磁样品的磁化率参数，是海底沉积物在地质历史上各时期沉积产物受磁化作用强弱的重要参数，与沉积环境的磁场变化、氧化还原条件等因素有关，可据此对地层沉积的阶段性进行划分。以下是钻孔0-10m岩芯中取得的12个样品数据：ID磁化率ID磁化率ID磁化率16.055.798.326.066.3107.735.375.3117.744.084.71210.3则样品数据矩阵为：Vector=[6.06.05.34.05.76.35.34.78.37.77.710.3]1、样品正规化v=std1(Vector)v=0.31750.31750.206300.26980.36510.20630.11110.68250.58730.58731.00002、计算极差矩阵[D,a,b]=range1(v);DD=7/73、最优K分割设K=4[array]=sect(v,4)array=7811即样品段Vector=[6.06.05.34.05.76.35.34.78.37.77.710.3]的最优四分割的分割点为8和11，即最优四分割方案为：｛6.06.05.34.05.76.35.3｝｛4.7｝｛8.37.77.7｝｛10.3｝五、结论本文首先对有序样品的最优分割发原理及算法进行了深入分析，在深入理解原理的基础上，设计了Matlab算法并独立完成了代码编写，使得这一较为复杂的数学计算问题能够在计算机程序中快速的计算，以得到满意的答案。但同时，算法设计与程序代码中仍存在一些问题，现列如下：1、用以衡量段间和段内样品数据差异的指标因不同的应用环境而应采用不同的指标，本文中作为示例仅以极差作为衡量指标。2、由于程序中使用了递归算法，使得当分割段数较多的情况下计算量迅速增大，计算时间较长，算法有待改进。3、未能及时进行视窗设计，实现软件化，所以实用代码仅能在Matlab命令窗里调用计算。