有序樣品的最優(yōu)分割的算法及其在MATLAB中的實(shí)現(xiàn)

1、 /7有序樣品的最優(yōu)分割算法及其在Matlab中的實(shí)現(xiàn)一、有序樣品聚類最優(yōu)分割的概念地質(zhì)數(shù)據(jù)中，有些樣品有一定的排列順序，如沿地層剖面采集的巖石標(biāo)本，由鉆孔取得的巖芯樣品，由測(cè)井曲線所得的數(shù)據(jù)，由巖體中心到圍巖的蝕變剖面的樣品等，它們是有序地質(zhì)變量，在對(duì)這些有序樣品進(jìn)行分類時(shí)，不能打亂樣品的前后次序。所以，一些不考慮樣品排列順序的數(shù)學(xué)處理方法，對(duì)此并不適用。有序樣品的聚類分析就是對(duì)有序樣品進(jìn)行分段的統(tǒng)計(jì)方法。對(duì)個(gè)有序樣品進(jìn)行分割，就可能有2-種1劃分方法，這每一種分法成為一種分割。在所有的這些分割中，有一種分割使得各段內(nèi)部之間差異性最小，而短語(yǔ)段之間差異性最大。這種對(duì)個(gè)樣品分段并使組內(nèi)離差平方

2、和最小的分割方法，成為最優(yōu)分割法。這類問(wèn)題的提法如下：設(shè)有一批（個(gè)）按一定順序排列的樣品，每個(gè)樣品測(cè)得項(xiàng)指標(biāo)，其原始資料矩陣x11xi2x1NxX(pxN)=21x22xiNxp1xp2xpN其中元素表示第個(gè)樣品的第個(gè)指標(biāo)的觀測(cè)值。現(xiàn)在要把此個(gè)樣品按順序（不破壞序列的連續(xù)性）進(jìn)行分割（分段或者分類）。其所有可能的分割法共有Cl+C2+C3+CN-1=2N-1-1N-1N-1N-1N-1種?，F(xiàn)在要求在所有分割中找出一種分割法，這種分割法使得各段內(nèi)樣品之間的差異最小，而各段之間的差異最大。各段內(nèi)部差異最小，即各段內(nèi)數(shù)值變化最小，段內(nèi)數(shù)值變化可用變差或者極差來(lái)表示，比如樣品段、變差：dij=j.xa

3、x（i,j）2ija=ia1jx1，j=r-T+Txaa=1表示樣品段、各樣品之間數(shù)值比較接近，反之,表示樣品段、各樣品之間數(shù)值比較接近，反之,極差：dij=a=l（maxixaBm叫Bj%）對(duì)于單指標(biāo)情況，則dijdij=(maximiniBj要各段內(nèi)部的差異最小，即所分成各段變差的總和（即段內(nèi)離差平方和，稱為總變差）為最小。所以總變差分解公式：所以S總=m=ij=l1(Xjl-X)2m=l=1jn=l1(Xjl-Xl)+(Xl-X)2ml=1jn=l1(Xjl-Xl)2+lm=1jn=l1(Xl-X)2+2lm=1jn=l1(Xijml=1jn=l1(Xjl-Xl)2+lm=1jn=l1(

4、Xl-X)2S嚴(yán)段間+S段內(nèi)X)(XX)其中：段內(nèi)段間D2ml=1jn=l1(Xij-Xl)(Xl-X)2m(Xl=1ij-Xl)jn=l1(Xl-X)2mXl=1Xl-X(nlXl-nlXl)0S段內(nèi)ml=1jn=l1(Xjl-Xl)2為段內(nèi)離差平方和S段間ml=1jn=l1(Xl-X)2為段間離差平方和=S段內(nèi)+S段間S=S-S段間總段內(nèi)對(duì)給定的個(gè)樣品，總是個(gè)固定的量。若使段內(nèi)離差平方和為最小，則段間離差平方和必為最大。所以，使總變差（段內(nèi)離差平方和）為最小的分割法就是最優(yōu)的分割法。二、有序樣品聚類的最優(yōu)分割意義最優(yōu)分割在地質(zhì)研究中是一個(gè)非常有用的手段，只要地質(zhì)體的某些地球化學(xué)特征存在規(guī)律

5、性的差異，采用最優(yōu)分割的數(shù)學(xué)樹立方法，就能按順序在最理想的地方進(jìn)行分段。通過(guò)對(duì)地層中采集的具某些地球化學(xué)特征樣品的最優(yōu)分割，能在地層的劃分對(duì)比中發(fā)揮重要的輔助作用；在找礦過(guò)程中，該方法更顯得天獨(dú)厚的優(yōu)勢(shì)，它能進(jìn)行蝕變、礦化及礦體的準(zhǔn)確分帶；對(duì)巖漿巖相帶劃分及演化序列的研究也十分有效。第四紀(jì)地層巖相和厚度變化很大，給第四紀(jì)地層劃分和對(duì)比地層帶來(lái)了困難，但第四系地層是不同階段和不同的環(huán)境條件下形成的，如重礦物、微量元素等地球化學(xué)特征的分布規(guī)律與一定的沉積階段和沉積環(huán)境相對(duì)應(yīng)，因此，最優(yōu)分割法的數(shù)據(jù)處理，可能是第四紀(jì)地層劃分和對(duì)比的有效方法。三、最優(yōu)分割的計(jì)算步驟及其計(jì)算公式數(shù)據(jù)正規(guī)化設(shè)原始資料矩陣

6、為X(X(pxN)x11X12X1NX21X22X1Nxp1Xp2XpN將矩陣中的元素將矩陣中的元素Xjj變換為：Zj=勺_詢（i=1、2、p；j=1、2、N）jmax1jnxijmin1jnxij而得矩陣(pxN)(pxN)Matlab代碼:functionstd=std1(vector)%對(duì)矩陣進(jìn)行標(biāo)準(zhǔn)化Matlab代碼:functionstd=std1(vector)%對(duì)矩陣進(jìn)行標(biāo)準(zhǔn)化%vector為待分割矩陣maX1=maX(vector);%對(duì)列求最大值min1=min(vector);%對(duì)列求最小值a,b=size(vector);%矩陣大小,a為行數(shù),b為列數(shù)forj=1:bst

7、d(j)=(vector(j)-min1)/(maX1-min1);end計(jì)算極差（或變差）矩陣由上述極差（或變差）計(jì)算公式得到矩陣：12dd12dd1233dd12NNdN1NMatlab代碼：functionD,a,b=range1（vector）%D返回計(jì)算所得的極差矩陣a,b=size（vector）;%求矩陣大小,a為行數(shù),b為列數(shù)k=a;%當(dāng)只計(jì)算單指標(biāo)數(shù)據(jù)時(shí)，k=a=1fori=1:bforj=i:bd(i,j)=maX(vector(k,i:j)-min(vector(k,i:j);endendD=d;3.進(jìn)行最優(yōu)二分割由矩陣計(jì)算全部分兩類的各種分割相應(yīng)的總變差，即對(duì)每一個(gè)(-

8、、)，求出相應(yīng)的總變差S(2；j)(j=1、2、m-1)m找出最小值，確定各子段的最優(yōu)二分割點(diǎn)a(),即S(2；a1(m)=min1.S(2；j)m11jm1m從而得出個(gè)樣品的最優(yōu)二分割x_、XMXM八、x)12a1(N)a1(N)+1NMatlab代碼：functionS,alp=divi2(vector,n)%最優(yōu)二分隔,S為最優(yōu)二分割各段的分割點(diǎn)，%a1記錄了二分割的序號(hào)d,a,b=range1(vector);alp=ones(n-1,b);%al(i,j)表示前i個(gè)樣品的第j次分割點(diǎn)S=zeros(b,b);form=2:bforj=1:m-1s(m,j)=d(1,j)+d(j+1,

9、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.進(jìn)行最優(yōu)三分割對(duì)于、，由aa(、)及矩陣分別計(jì)算：Sm(3;a1(j),j)=Sj(2;a1(j)+dj+1,Nm1(j)j1(j)j+1,N(m=N、N-1、4、3)然后求出最小值，即Sm(3；a1(m)，a2(m)=min2jm1Sm(3；1j，j)從而得到N個(gè)樣品最優(yōu)三分割x、X“zx八、X,x、12a1(N)a1(N)+1a2(N)a2(

10、N)+！N最優(yōu)分割完全類似4的辦法，在最優(yōu)三分割的基礎(chǔ)上可以進(jìn)行最優(yōu)四分割，繼而進(jìn)行五分割，以此類推，如果已經(jīng)作出最優(yōu)分割，則可以產(chǎn)生最優(yōu)分割。Matlab代碼：functionS,alp=divi(vector,n)%門為要分割的段數(shù)d,a,b=range1(vector);alp=zeros(1,b);%al(i,j)表示前i個(gè)樣品的第j次分割點(diǎn)form=n:bforj=n-1:m-1ifn=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_t

11、emp(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);endendendfunctionarray=sect(vector,n)%vector為樣品矩陣，當(dāng)直接對(duì)樣品矩陣進(jìn)行分割時(shí)調(diào)用該函數(shù)%array返回樣品最優(yōu)n分割的分割點(diǎn)號(hào)%門為要分割的段數(shù)a,b=size(vector);fornum=n:-1:2S,alp=divi(vector,num);ifnum=narray(num-1)=alp(1,b);elsearra

12、y(num-1)=alp(array(num);endendfunctionarray=fsect(filename,n)%filename為需要分割的樣品數(shù)據(jù)的文件名，如d:temp.txt%要對(duì)數(shù)據(jù)文件進(jìn)行分割時(shí)調(diào)用該函數(shù)fid=fopen(filename,r);A,count=fscanf(fid,%f);vector=A;a,b=size(vector);fornum=n:-1:2S,alp=divi(vector,num);ifnum=narray(num-1)=alp(1,b);elsearray(num-1)=alp(array(num);endend四、實(shí)例分析瓊州海峽標(biāo)志性

13、鉆孔古地磁樣品的磁化率參數(shù)，是海底沉積物在地質(zhì)歷史上各時(shí)期沉積產(chǎn)物受磁化作用強(qiáng)弱的重要參數(shù)，與沉積環(huán)境的磁場(chǎng)變化、氧化還原條件等因素有關(guān)，可據(jù)此對(duì)地層沉積的階段性進(jìn)行劃分。以下是鉆孔巖芯中取得的個(gè)樣品數(shù)據(jù):磁化率磁化率磁化率則樣品數(shù)據(jù)矩陣為:、樣品正規(guī)化、計(jì)算極差矩陣00D.J1UD.3175D-31V60.36510.36510.36510.63250.6S250.6S25i.ooao00D.1LD.31T5D-31750-36510.3651IX3651i：L6325i：L6825i：L6825i,ooai：i000D.2D030.26930.36510.30510.30510.6825f

14、l.6825fl.A825i.ooao00000.2693fl.36510.3fi510.3fi510.6825fl.6825fl.開(kāi)25.0000000000.05520.15870.26400,ET140.E7140,67140.88390000|：00.158?0.25400.5T140.57140-5714|：|.3S39U0UUQaa0BD952CL57140-57140.57140.88390000000Q0.0.5714fl.5714fl.88990n00000000.D9520,096E0.4J2T0000|：aa00DD0.412T0000|：aa00DD0.41270000r03最優(yōu)分割設(shè)即樣品段的最優(yōu)四分割的分割點(diǎn)為8和11，即最優(yōu)四分割方案為：0.5五、結(jié)論本文首先對(duì)有序樣品的最優(yōu)分割發(fā)原理及算法進(jìn)行了深入分析，在深入理解原理的基礎(chǔ)上