聚类分析法

应用统计学课程论文：

浅析聚类分析法在碎纸拼接复原中的

应用

系别：工程管理 学号：B11080111 姓名：贾晓婷 指导老师：张亚峰

1

浅析聚类分析法在碎纸拼接复原中的应用

系别：工程管理 学号：B11080111 姓名：贾晓婷 指导教师：张亚峰 摘要：应用统计学是一门以收集数据分析数据为研究对象的学科，聚类分析法作为数据分类的一种典型方法，也是应用统计学一种基本的方法。破碎文件的拼接复原以其在司法物证复原、历史文献修复以及军事情报等方面的突出作用，一直是人们的研究热点。本文以其在上述应用领域的重要性为背景，以2013年大学生数学建模竞赛B题中问题1为例，针对典型的碎纸机破碎纸片的拼接复原问题建立了聚类分析的模型，具体介绍了聚类分析法的原理以及在碎纸复原中的应用。

针对所选问题，我们建立了基于聚类分析的碎纸片复原模型。首先分析中文碎纸片。对题目所给图像数据进行初步处理，将图片用Matlab函数读入Matlab，把图片信息转换为数字矩阵。再提取每一个数字矩阵的首列和尾列作为每一个碎纸片的特征边，将所有特征边组合成一个新的矩阵，然后筛选出起点纸片和终点纸片，结果分别为008和006。接着对该矩阵的每一列进行R型聚类分析，分析结果见正文，进而得到聚类图见正文图5-1，复原结果见正文表5-1。英文碎纸片的复原过程与中文碎纸片一致，起点纸片003，终点纸片004，聚类图见正文图5-2，复原结果见正文表5-2. 关键字：聚类分析法 碎纸拼接复原

一、 问题背景及分析

所选问题要求对给定的来自同一页印刷文字文件的碎纸机破碎纸片（仅纵切）进行拼接复原，根据所给附件我们可以了解到，所有的碎纸片规格相同，这样我们就不能用基于轮廓的复原方法，采用基于内容的匹配复原法。由于聚类分析即可以对样本进行分类，也可以对变量进行分类，分类的依据是样本或变量间的相近程度。由切割的性质可以得到相邻的碎纸片的公共边必然具有完全相同的特征，于是根据每一条碎纸片边缘不同的文字特征对其进行聚类分析具有合理性，只需令所分的组数等于需要重组的边数即可得到得到较理想的匹配结果。分类结果以组为单位，以起点纸片标号为开头，依次向下寻找特征边数字矩阵距离最近的纸片，查找原则是奇偶相间。

考虑到结果可能具有的不唯一性，因为数据量不是很大，我们选择将所有连续的序列链列举出来，然后人工筛选出理想解。

二、 模型的建立与求解

2.1 模型假设

假设所有的碎纸片在切割过程中没有任何缺损，即可以被完整复原； 2.2 符号约定

2

A、B、C、D i,j 2.3 模型建立与求解 2.3.1 中文碎纸片复原

表示数据矩阵 用以计数的循环变量 注:若解题过程中出现重复变量或新的变量，按当前变量含义为准。

由于被切割的碎纸片的相邻边具有唯一匹配性，本问选择R型聚类分析法对碎纸片进行复原。

所谓聚类分析，是按照一定的标准对研究对象进行分类的数学方法，聚类分析又称群分析，是对多个样本或指标进行定量分类的一种多元统计分析方法。对样本进行分类成为Q型聚类分析，对指标进行分类称为R型聚类分析。聚类分析的方法有很多种，如最短距离法、最长距离法、中间距离法、重心法、类平均法、可变类平均法、可变法、利差平方和法等等。本文使用的是R型的最短距离法对碎纸片进行分类。

第一步，数据的预处理。

将图片读入Matlab，把图像信息转换为数字信息，提取每一个图像转换后的数字矩阵的第一列和最后一列，认为它们是碎纸片的特征边，依次按顺序重新组合并转置成为一个新的矩阵。每一个图像转换后的数字矩阵都是1980*72的矩阵，故新矩阵应为38*1980的矩阵。新矩阵如下（仅显示部分结果，所有元素不全为255）：

第二步，首尾碎纸条的筛选。

由观察结果可知，起点纸片为008，终点纸片为004. 第三步，使用Matlab进行聚类分析并作出聚类图。 对记每一行（周期）的指数为:wi来测量点与点之间的距离，即

d(vi1,vi2,,vi240),i1,2,,97。使用绝对值距离

w,wvvijk1ik240jk,DG,Gmindw,w.

pqwiGpwjGqij由以上距离公式，可以算出距离矩阵：

a11Dai1其中ij38。

620190a1j203490aij2552061355295745982835765438022991032386035599316160

使用Matlab进行聚类分析，得到如下结果（分析结果中的数字代表矩阵的相应列）：

3

第1次合成，平台高度为5649时的分类结果为：5 22 第2次合成，平台高度为6328时的分类结果为：23 38 第3次合成，平台高度为6856时的分类结果为：2 13 第4次合成，平台高度为7503时的分类结果为：26 31 第5次合成，平台高度为8011时的分类结果为：6 33 第6次合成，平台高度为9948时的分类结果为：28 37 第7次合成，平台高度为10835时的分类结果为：1 36 第8次合成，平台高度为10955时的分类结果为：15 24 第9次合成，平台高度为11274时的分类结果为：3 34 第10次合成，平台高度为12842时的分类结果为：12 19 …… ……

第30次合成，平台高度为22198时的分类结果为：26 34 第31次合成，平台高度为22325时的分类结果为：26 27 第32次合成，平台高度为22379时的分类结果为：11 34 第33次合成，平台高度为23296时的分类结果为：30 35 第34次合成，平台高度为23444时的分类结果为：31 37 第35次合成，平台高度为23884时的分类结果为：26 28 第36次合成，平台高度为25181时的分类结果为：11 32 第37次合成，平台高度为26499时的分类结果为：26 37 使用Matlab统计工具箱的相关命令画出聚类图如下：

4

11109876543x 104 711 124 510142625182923301528 8211219 616222027 213 4 9 317图1 中文碎纸片特征边聚类图

相应分类结果为： 第1类的有4 第2类的有9 第3类的有25 30

第4类的有14 17 26 31 第5类的有10 11 第6类的有18 29 第7类的有23 38 第8类的有15 24 第9类的有28 37 第10类的有5 22 第11类的有8 21 第12类的有12 19 第13类的有6 33 第14类的有16 35 第15类的有1 36 第16类的有20 27

5

第17类的有2 13 第18类的有7 32 第19类的有3 34 第四步，碎纸片拼接复原。

根据第三步的结果人工筛选出可以唯一组成完整图像的两条序列如下：

（1） 008（17、18）→014（29、30）→012（25、26）→015（31、32）→003(7、8)→010（21、22）→002（5、6）→016（33、34）→001（3、4）；

（2） 004（9、10）→005（11、22）→009（19、20）→013（27、28）→018（37、38）→011（23、24）→007（15、16）→017（35、36）→000（1、2）→006（13、14）。

由第二步可知，起点纸片为008，终点纸片为006，可直接将两条序列链连接为一条，最终复原结果如下表：

表1 附件1中文碎纸片复原结果

008 014 012 015 003 010 002 016 001 004 005 009 013 018 011 007 017 000 006 运用Matlab按以上顺序重组复原图数字矩阵，得到中文复原结果见下图中文复原图：

6

图2 中文复原图

2.3.2 英文碎纸片复原

英文碎纸片的复原过程与中文碎纸片复原过程相似，参考中文碎纸片的复原过程，得到英文碎纸片重要过程的结果如下：

经筛选，得起点纸片为003，重点纸片为004。 Matlab进行聚类分析结果为：

第1次合成，平台高度为5671时的分类结果为：8 13

7

第2次合成，平台高度为6180时的分类结果为：33 36 第3次合成，平台高度为6756时的分类结果为：3 12 第4次合成，平台高度为6953时的分类结果为：6 15 第5次合成，平台高度为7007时的分类结果为：18 33 第6次合成，平台高度为7747时的分类结果为：32 37 第7次合成，平台高度为8941时的分类结果为：2 11 第8次合成，平台高度为9572时的分类结果为：4 19 第9次合成，平台高度为9898时的分类结果为：9 34 第10次合成，平台高度为10076时的分类结果为：1 24 …… ……

第30次合成，平台高度为23932时的分类结果为：3 19 第31次合成，平台高度为23955时的分类结果为：1 17 第32次合成，平台高度为24111时的分类结果为：1 28 第33次合成，平台高度为24121时的分类结果为：1 31 第34次合成，平台高度为24172时的分类结果为：1 30 第35次合成，平台高度为24264时的分类结果为：1 21 第36次合成，平台高度为24314时的分类结果为：1 13 第37次合成，平台高度为24416时的分类结果为：1 34 聚类图为：

8

35

4x 10765432 312 124 5 61516 71810212825 813202717302329 4 919142226 211图3 英文碎纸片特征边聚类图

分类结果为：

第1类的有5 7 10 第2类的有1 24 第3类的有3 12 第4类的有6 15 第5类的有16 31 第6类的有18 33 36 第7类的有32 37 第8类的有21 28 第9类的有25 第10类的有8 13 第11类的有20 27 第12类的有17 30 35 第13类的有23 38 第14类的有4 19

9

第15类的有9 34 第16类的有14 第17类的有22 第18类的有26 29 第19类的有2 11

根据分类结果人工筛选出唯一可以组成完整图像的三条序列如下： (1)003(7、8)→006（13、14）；

(2)002（5、6）→007（15、16）→015（31、32）→018（37、38）→011（23、24）→000（1、2）→005（11、12）→001（3、4）→009（19、20）→013（27、28）→010（21、22）；

(3)008（17、18）→012（25、26）→014（29、30）→017（35、36）→016（33、34）→004（9、10）。

由于起点纸片为003，终点纸片为004，直接将三条序列链按（1）→（2）→（3）连接为一条，最终复原结果如下表：

表5-2 附件1英文碎纸片复原结果表

003 006 002 007 015 018 011 000 005 001 009 013 010 008 012 014 017 016 004 再次运用Matlab按以上顺序重组复原图数字矩阵，得到英文复原结果见下图英文复原图：

10

图2 中文复原图：

三、 模型评价

所选问题要求对给定的来自同一页印刷文字文件的碎纸机破碎纸片（仅纵切）进

11

行拼接复原，根据所给附件了解到，所有的碎纸片规格相同，由于聚类分析即可以对样本进行分类，也可以对变量进行分类，分类的依据是样本或变量间的相近程度。由切割的性质可以得到相邻的碎纸片的公共边必然具有完全相同的特征，于是根据每一条碎纸片边缘不同的文字特征对其进行聚类分析具有合理性，只需令所分的组数等于需要重组的边数即可得到得到较理想的匹配结果。由问题一的聚类分析结果可知，聚类分析法在此处的运用效果较好。

四、 参考文献

[1]卓金武等,《MATLAB在数学建模中的应用》,北京航空航天大学出版社,2013年1月第1次印刷

[2]司守奎等,《数学建模算法与应用》,国防工业出版社,2011年8月第1版 [3]张德丰等,《Mtalab程序设计与综合应用》,清华大学出版社,2012年1月第1版 [4]姜启源等，《大学数学实验》，清华大学出版社，2005年2月第1版 [5]冯 杰等，《数学建模原理与案例》，科学出版社，2007年1月第1版

12

附录：

1 数据预处理程序

a=imread ('000.bmp');%将默认位置设为目标碎纸片所在的附件 b=imread ('001.bmp'); c=imread ('002.bmp'); d=imread ('003.bmp'); e=imread ('004.bmp'); f=imread ('005.bmp'); g=imread ('006.bmp'); h=imread ('007.bmp'); i=imread ('008.bmp'); j=imread ('009.bmp'); k=imread ('010.bmp'); l=imread ('011.bmp'); m=imread ('012.bmp'); n=imread ('013.bmp'); o=imread ('014.bmp'); p=imread ('015.bmp'); q=imread ('016.bmp'); r=imread ('017.bmp'); s=imread ('018.bmp');

aa=[a(:,1),a(:,72),b(:,1),b(:,72),c(:,1),c(:,72),d(:,1),d(:,72),e(:,1),e(:,72),f(:,1),f(:,72),g(:,1),g(:,72),h(1,:),h(72,:),i(1,:),i(72,:),j(1,:),j(72,:),k(1,:),a(:,72),l(:,1),l(:,72),m(:,1),m(:,72),n(:,1),n(:,72),o(:,1),o(:,72),p(:,1),p(:,72),q(:,1),q(:,72),r(:,1),r(:,72),s(:,1),s(:,72)];

2 聚类分析法程序 clc Clear A=[];

a=[A’];%a为特征边数字矩阵的转置矩阵 [m,n]=size(a); d=zeros(m);

d=mandist(a');%mandist求矩阵行向量之间的两两绝对值距离

13

d=tril(d);%截取下三角 nd=nonzeros(d); nd=union(nd,nd); for i=1:m-1 nd_min=min(nd);

[row,col]=find(d==nd_min);tm=union(row,col); tm=reshape(tm,1,length(tm));

fprintf('第%d次合成，平台高度为%d时的分类结果为：%s\\n',i,nd_min,int2str(tm)); nd(find(nd==nd_min))=[]; if length(nd)==0 break end end

3 聚类图画法 clc clear

a=[];%a为特征边数字矩阵的转置矩阵 y=pdist(a,'cityblock'); yc=squareform(y) z=linkage(y) [h,t]=dendrogram(z)

T=cluster(z,'maxclust',19)%数字表示所分成的组数 for i=1:19 tm=find(T==i);

tm=reshape(tm,1,length(tm));

fprintf('第%d类的有%s\\n',i,int2str(tm)); end

4 复原图做法

a=imread ('000.bmp');%需将默认位置设为目标碎纸片所在的附件 b=imread ('001.bmp'); c=imread ('002.bmp'); d=imread ('003.bmp');

14

e=imread ('004.bmp'); f=imread ('005.bmp'); g=imread ('006.bmp'); h=imread ('007.bmp'); i=imread ('008.bmp'); j=imread ('009.bmp'); k=imread ('010.bmp'); l=imread ('011.bmp'); m=imread ('012.bmp'); n=imread ('013.bmp'); o=imread ('014.bmp'); p=imread ('015.bmp'); q=imread ('016.bmp'); r=imread ('017.bmp'); s=imread ('018.bmp');

C=[i o m p d k c q b e f j n s l h r a g];%求解英文复原图时C替换为D=[d g c h p s l b j n k i m o r q e];

imwrite(C,'ͼһ.bmp','bmp');%求解英文复原图时C替换为D，'ͼһ.bmp'为'文件名

15