2017数学建模A题

2017数学建模A题

CT系统参数标定及成像

摘要

本文根据CT系统的成像和获取物体信息的问题，使用主成分分析、解释结构模型、滤波反投影等方法，综合分析了CT系统的参数标定与其成像原理，分别构建了ISM解释结构模型、直接反射投影模型。然后使用Excel、MATLAB等软件整理数据和编写程序，详细的阐明了CT系统对未知结构体进行扫描成像的理论，并且根据其扫描数据反推出未知物体的基本信息，充分证明了所设计方案的合理性和可操作性。

针对问题一的要求，本文首先以CT系统基本指标、旋转方式以及附件一、二所给出的数据进行CT系统的位置及具体指标进行确认，其中，运用Excel工具对附件一、二进行数据处理进行辅助计算和判断。在此基础上建立主成分分析模型对各项指标进行主成分分析，确定其中特殊位置（吸收率最大和最小）的数值模块进而得出CT系统在正方形托盘中的位置，激光束单元之间的距离0.2857mm，以及CT系统旋转180次所转过的角度为180°。

针对问题二、三的要求，为了了解未知介质的位置、几何形状、吸收率。我们可以运用直接反射投影法，建立直接反射投影模型，即利用MALAB 工具对Excel表格数据进行导入，编写反投影对应的程序，得出其几何图形。并且结合问题一中的CT系统基本信息，我们通过衰变公式结合附件三、五的吸收率数据推出其在托盘中的位置，通过未知物体的位置从而给出10个规定位置的点的吸收率为：

测器上512个等距单元之间存在一定的距离，此等距单元间距会存在一定的误差，并且我们在对附件二的数据进行图表分析的时候发现其折线图存在着突变，即得出其稳定性和精度存在一定的缺陷，我们针对已有的问题，建立新的优化模型加以改进即垂直X轴为初始位置，以（50,50）为旋转中心进行断层扫描。

本文给出了此CT系统在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获得样品内部的结构信息的特点。并通过对几个未知样品的分析来验证此CT系统的

功能和作用。且在后续的论述中我们对此CT系统的不足进行了改进和更新，解决了此CT系统在标定精度和稳定性的缺陷，使得此系统更加的完善。

关键字：断层成像；吸收率；直接反射投影；PCA；Excel；MATLAB

§1问题的重述

一、背景知识

1.背景的介绍

随着医疗和工业的发展，CT技术在我们的生活中日益推广，在疾病更加复杂，设计更加深奥的情况下，人们的探索也越来越深刻，于是随之而诞生的一项新技术------CT成像诞生了，随着许多年来人们对CT成像的不断探索，我们所掌握的技术也越来越高级，对于许多未知的事物探索也越来越广泛，于是这项技术被应用于各个邻域，尤为显得突出的是医学领域和工业领域。在医学领域中，对病人全身的身体状况进行探索，对各种人眼难以判断和简单仪器难以探测的部位可以采用CT新型技术对此进行处理；在工业中，我们通过CT技术对一些位置物体的内部情况或者在不破坏其结构的情况下，对其进行断层成像，便于人们的技术研究。

2.问题的产生

在工业中，我们对一个未知物体进行CT扫描处理，会得到许多的数据信息，在这个数据信息分析的同时，我们要确定此未知物体的几何形状介质分布情况以及许多未知的信息。在许多CT系统中，大多数CT系统均采用的是光线穿透法，通过穿透率强的几种光线来对物体内部的结构进行判定，于是我们需要根据对许多组数据的处理和计算才能得出物体的基本信息。

3.已有的对策

我们通过查阅资料和阅读文献可以知道，在CT系统断层成像的过程中，我们可以认为的控制其扫描位置以及其扫描角度，并且对扫描过程中的能量吸收率有一个具体的数值计算。我们可以根据CT探测器的基本已知指标和所测得的数据进行统计分析即可得出结论，接下来进行具体分析。

二、题目重述

CT可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。一种典型的二维CT系统如图1所示，平行射入的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。

三、待解决的问题

1.在正方形托盘上放置两个均匀固体介质组成的标定模

板，模板的几何信息如图2所示，相应的数据见附件1，其中

每一点的数值反映了该点的吸收强度，这里称为“吸收率”。对

应于该模板的接收信息见附件2。请根据这一模板及其接收信

息，确定CT系统旋转中心在正方形托盘中的位置、探测器单

元之间的距离以及该CT系统使用的X射线的180个方向。

2.附件3是利用上述CT系统得到的某位置介质的接收信

息。利用（1）中得到的标定参数，确定该未知介质在正方形

托盘中的位置、几何形状和吸收率等信息。另外，请具体给出

图3所给的10个位置处的吸收率，相应的数据文件见附录4.

3.附件5是利用上述CT系统得到的另一个未知介质的接

收信息。利用（1）中得到的标定参数，给出该未知介质的相

关信息。另外，请具体给出图3所给的10个位置处的吸收率。

4.分析（1）中参数标定的精度和稳定性。在此基础上自

行设计新的模板、建立对应的标定模型，以改进标定精度和稳

定性，并说明理由。

图1. CT系统示意图图2. 模板示意图（单位：mm）图3. 10个位置示意图

§2问题的分析

一、 问题的总分析

表一

二、 对问题的具体分析

1. 对问题一的分析

本题要求我们确定CT 系统旋转中心在正方形托盘中的位置、探测器

单元之间的距离以及该CT 系统使用的X 射线的180个方向。我们从探测

器单元距离入手，通过Excel 数据分析得到图二右边物体的小球直径和光

线射数的比值即为单元距离。从而通过单元距离判定吸收率最大处的射入

光线的编号从而找到其入射的180个角度，在进而计算出其在托盘中的位

置。

CT 系统参考系数标定及成像

问题问题得出能量吸收MATLAB 得出图形形对探测器单元通过单元距离问题从精度和稳定进行多层阶分

2.对问题二、三的分析

两题要求我们根据问题一所得的CT系统的具体参数从而对放在正方形托盘上的未知物体的位置、几何形状、吸收率进行求解，并且还要给出规定的10处位置的吸收率。我们需根据问题一求出其能量变化方程，通过直接反射投影法用MATLAB导入Excel的数据从而得出物体的几何图形，通过几何图形和附录中的数据位置和值的大小通过求得的方程进行推导出物体的位置和其吸收率，以及给出的10个点的吸收率。

3.针对问题四的分析

我们考虑到在此给出的CT系统条件下，其对未知物体探索的过程中，会导致物体精度和稳定性出现偏差和波动，于是我们通过多层阶分析得出物体具体的参数指标，然后我们在此分析基础上对其进行修改建立一个新的模型。

§3模型的假设

1.假设我们收集整理的数据都是真实有效的。

2.假设除我们选取的因数外，其余影响因数的影响效果

可以忽略不计。

3.假设物体不掺杂其他介质。

4.假设有条件支持我们构建新的CT模型。

§4名词解释和符号说明

直接反射投影：断层平面中的某一点的密度值可以看做是这一平面内所有经过该点的射线的投影值之和。

PCA：主成分分析。通过正交变换将一组可能存在相关性的变量转换为一组线性不相关的变量，转换后的这组变量叫主成分。

§5模型建立与求解

一、对问题一的分析与求解

1.对问题一的分析

考虑到附件所给信息数据过多，于是我们通过对附件一和二数据进行筛选，在Excel中将空吸收率即“吸收率”为0的数据全部

用空格来替换掉，再对整体视图进行缩小，即可大致看出其几何

图形和扫描成像。根据题目要求，我们需要确定CT系统旋转中

心、探测器单元之间的距离以及X射线的180个方向。对比三个

提问，我们可根据缩小视图大致推断出其为中心旋转半周的工作

期，并且在工作期中，每一次都有单元射线入射进两个物体，于

是我们可运用主成分分析法对吸收强度最大和最小的情况进行分

析即可得出两探测器之间的单元距离，随后利用逆推法得出其初

始值，并计算出旋转中心。

2.对问题的求解

1）模型的建立

我们根据题目所给要求，确立所需求得的变量，CT系统旋转中心（X0，Y0），探测器单元之间的距离d，X射线的180个方向n （1,2,3……180）。

第一步：我们首先对探测器单元之间的距离进行求解。根据缩小视图进行分析可知，在旋转过程中均有射线射过右侧图形，并且每次射入的光束数量大致相等，结合附件一视图，可以推论出其为圆形，且每一个角度光束均穿过其直径。可得如下公式：

d=小球直径

入射光线的数量?1

第二步：我们确定了d值以后可以通过图2所给距离长度推出在最大值时探测器的对应旋转次数，从而通过顺时针旋转此次数即为初始位置，即可得出180组方向。

第三步：通过前两步确定的值算出探测器长度和相对正方形托盘的位置

以吸收率最大的光束对应的光束编号来推测中间光束的位置，然后顺时针旋转90°，两条中心光束的交点即为旋转中心。

2）模型的求解

先对附件一和附件二进行“吸收率”为0的值得替换处理，并进行视图放缩，可以得到如下两组视图：

图4 图5

通过图分析可知，其左侧图形大致为椭圆形图形，于是可以实行模型建立中第一步求间距d，在对图分析值，其代表右侧图形的视图显示，经过其内部的射线条数大致为固定值，我们进行数据分析得射线个数为众数

a=29

而由图2知，小球直径为8mm。所以可得：

d=8

28=2

7

mm

由此可以继续讨论第二步证明

我们通过Excel进行分析，找出吸收率最大的值，确定其位置为（ES，222），可知在最大吸收率处，旋转ES次后第222个单元射线射入椭圆长轴，即此时探测器接收板与正方形托盘下侧平行，即为S处。我们将ES转换为数字得

ES=149

即CT系统逆时针旋转149次到达了S处，我们通过数据转置大致得到下图

图6

图7

可以看出此为对称图形，即可以得到此CT系统大致旋转了180°，即

每一次大约旋转1°，即初始位置如图，第222条射线在左上角与椭圆形短

轴左侧成60°角处，即可确定180次旋转角度为从初始位置开始，逆时针

以每一度开始旋转180次。

我们再来探讨第三步，已知探测器由512个等距单元，可以求得探测

器的长度：

d 探=512?d =146.3mm

d 探>100mm

所以可知探测器必定超出正方形木板边缘，我们找到探测器板正中间射线

编号为512/2+1=257

Δd= (257-222 )* 27

=10mm

经过顺时针旋转90°可以得到探测板与正方形木板左侧平行位置，此时257

编号距托盘上侧距离为40mm ，我们以椭圆中心为原点建立直角坐标系，

可以得到中心旋转的位置为（-10,10）处。

图8

二、 对问题二，三的求解

1. 对问题二、三的分析

我们对问题二、三进行分析，两个提问大体上都要我们通过第一问求出来

的CT 系统的已知位置和基本标量对两个未知物体进行扫描判断出其几何

形状，托盘中所处的位置以及其吸收率。并且也额外给出10个已知坐标的

点求这些点的吸收率。通过对问题的对比，我们大致可以看出在图形的生

成和内部结构上存在较大的差异。接下来进行问题的求解。

2.对问题二、三的求解

1）模型的建立

在问题一的基础上，我们已知CT系统的旋转方式和基本指标。通过其扫描此未知介质可以得到其断层信息，在附件三、五中显示的是此两个未知物体的接收信息数据，我们通过MATLAB软件导入此两组数据以repmat 和floor编写程序生成二维几何图显示此两者未知物体的几何形状然后我们通过矩阵降维将此512*180的矩阵减低成256*256的正方形矩阵。通过吸收数值与空值来确定其在托盘中的位置。

2）模型的求解

首先对附件三进行分析，我们将其导入MATLAB中，通过直接反投影函数得出其初始基本形状如下：

图9

我们可以看出附件三给出的未知介质的基本图形大体上看为一个椭圆且在一边有空缺。我们通过对附件三进行降维处理，将此512*180的矩阵降阶成256*256的正方形矩阵，然后将其成比例的映射到正方形托盘上，按照做盘边缘建立二维坐标系，得到下图：

图10

根据如图所示的坐标刻度和其图形在坐标轴中占的位置可以描述出其在正方形托盘中的位置。

同理，对第三问进行处理。我们将附件五导入MATLAB中，通过直接反投影函数得出其初始基本形状如下：

图11

我们大致可以看出，其内部结构呈现网面状，初略表示为海绵体的内部结构。同样的，我们将之进行维度降维处理。然后将其成比例的映射到正方形托盘上，按照做盘边缘建立二维坐标系，得到下图：

图12

根据如图所示的坐标刻度和其图形在坐标轴中占的位置可以描述出其在正方形托盘中的位置。同样的，我们在对附件三和附件五进行矩阵维度处理的同时，我们将其变成256*256的矩阵，此换成的表格数据对应的值即为未知物体的吸收率。（见文件problem2.xls和problems3.xls）

另外的，题目给出了10个已知点的坐标位置，我们根据附录四和图三可以建立如下坐标图：

图13

我们根据衰变公式可以知道，放射性物质的原子核改变速度dN/dt 与当时的剩余的衰变原子核的总数N成正比,即

dN/dt=?kN

对原式进行积分

t=?1

?ln|N|+c

ln|N|?kc=?kt N=N0,t=0代入得c=1

k

?lnN0所以

ln(N 0)=?kt or N =N0?e ?kt N 为剩余的衰变原子核数，N 0为初始衰变原子核数，k 为衰变常数，t

为时间。若我们将物体均匀分为多段，各段物体的线性衰减系数分别为（μ1，μ2，…，μn ），其每段线段的长度分别为L 1，L 2，…L n ，则I =I 0e ?（μ1?L 1+μ2L 2+?+μn L n ）

上式可化为1N

J =-∑μj L j =Ln （I/I 0）

更一般地，如果物体在平面内都不均匀，衰减系数分布为f （x,y ）。在某一方向路径L 的射线变化强度为：

I =I 0e ?∫f (x,y )dl L

通过以上推导出的公式，我们根据每个点的对应坐标可以得到其吸收率如下图：

表三

三、对问题四的求解

1.对问题四的分析

在对上三问探讨的过程中，我们发现在进行成像还原的同时，会出现部分图像丢失的情况，于是我们在对其精度的要求上加以改进，以保证其在降维过程中不会出现物体投像超出储存单元的情况。并且我们在对此CT 系统的稳定性加以分析的同时，对附件二进行处理，对其所有平均值进行分析得到下列的折线图：

图14

由图可知，在一定的变换角度时，会出现较大的波动变化，这说明此种扫描稳定性存在一定的缺陷。我们需要寻找一种新型的扫描的CT系统来确保在人工可控的范围内进行图像成像，来优化此类算法。

2.对问题四的求解

1）模型的建立

我们根据图14看出，此CT系统扫描出来的图像的值十分不稳定，波动性较大，于是我们在此基础上不变的情况下我们改变探测板的位置，即探测板中心对应的照应于正方形托盘的中心，然后我们使此系统围绕正方形托盘中心进行旋转，可以提高此模型的精度和稳定性。

2）模型的求解

我们对附件二的数据进行分析，我们对其数据进行压缩，并假定第256个射线扫描点对应直射入正方形托盘的中心。即R256对应射于（50,50）的位置，然后我们以垂直X轴为初始位置，以（50,50）为旋转中心进行断

层扫描得到一组新的大数据，然后我们继续求此180组数据的平均值然后形成折线图得到如下的图形。

图15

可以看出，图15相比于图14，期稳定性得到了明显的改进，并且精度更加准确，能够更加线性的表示出数据的基本情况。可以得出，在此改进上可以使得此CT系统的基本性能加以提高，便于对应用上的更深入研究。

§6模型的评价及改进

一、模型的评价

（1）误差分析

1.收集的数据为真实数据，受到外界各方面的影响，存在一定

误差；

2.为了使结果更理想，对数据进行主观分析，进行了一些必要

地处理，带来一些误差。

（2）模型优点

1.在对问题的研究中，我们通过各种模型假设和模型方法对问题进行了

讨论和验证，结果不乏严谨性。

2.我们对基本模型进行了分析的同时给出了自己的观点，同时在此模型

不足上加以改进，建立了一个更加完整的模型，使得问题的探究更加深入，

能更好的与现实生活结合起来。

二、模型的改进

在现实生活中，未知物体的形状是多变的，并且在对物体进行扫描的同时，我们对物体构成的介质情况是不了解的，现实生活一个物体的物质构成是多元的，也就是说我们通过此种方法只能显示出其大致的形状，不能完全的反映出内部的结构构造。最重要的一点是，在生活中，物体的维度

是三维的，我们在对物体进行研究的时候要考虑到物体的立体特性，即在对其进行探索的时候要有一个全方位穿透的扫描。所以我们在假设条件允许我们构建新的模型的情况下，我们对CT系统的组成成分进行改进，具体参考问题四的解决。就对模型而言，我们通过直接反投影法对MATLAB 图形还原的时候，不能精确的反映具体部位的吸收情况，并且在对其生成图像的结果中，图形的分辨率和尺度位置都存在一定的不足。就此问题，我们需要进一步的研究，通过查阅大量资料，进行具体操作的详细学习，从而能够更加完善的对数据进行处理，进而做出更加满意的问题答复。

§7参考文献

[1] 闫鋲，李磊；CT图像重建算法；科学出版社

[2] 姜启源，谢金星，叶俊；数学模型；高等教育出版社；2003.8

[3] 恒盛杰资讯；Excel数据处理与处理经典；中国青年出版社；2007

[4] 杨桂元；数学建模[M]；上海财经大学出版社；2015

[5] 余胜威；MTALAB数学建模经典案例实战；清华大学出版社；2014

[6] 吴建国；数学建模案例精编；中国水利水电出版社；2005

[7] 周品，赵新芬；MATLAB数理统计分析；国防工业出版社；2009.4

[8] 周品；MATLAB数学计算与仿真应用；电子工业出版社；2013.1

[9] 马雪伶；Excel 2010实战技巧精粹；人民邮电出版社；2013

[10] 杨桂元，朱家明，数学建模竞赛优秀论文评析[M]．合肥：中国科学技术大学出版社，2013．

[11]直接反投影法https://www.360docs.net/doc/a814419876.html,/p-1844302741747.html

[12]图像重建原理算法基础

https://https://www.360docs.net/doc/a814419876.html,/view/88a85383e53a580216fcfebd.html

§8附录

程序1：

a=A;

b=imrotate(a,90,'crop')

P =b;

R=radon(P);

figure;imshow(R,[]);

figure;

imshow(P,[]);title('yuan图');

%反投影

P_3=zeros(l,l);

for i=1:l

for j=1:l

for k=1:180

theta=k/180*pi;

t=(j-l/2-0.5)*cos(theta)+(l/2+0.5-i)*sin(theta)+(H+ 1)/2;

t1=floor(t);

t2=floor(t+1);

P_3(i,j)=P_3(i,j)+(t2-t)*x(t1,k)+(t-t1)*x(t2,k);

end

end

end

P_3=pi/N*P_3;

figure;imshow(P_3,[]);title('滤波反投影法');

程序2：

%直接反投影法

l = pow2(nextpow2(size(a,1))-1);%重构图像的大小

P_2 = zeros(1,1);%用于存放重构后的图像

for i=1 : size(a,2)

tmp = imrotate( repmat(a(:,i),1,size(a,1)),i-1,'bilinear' );

tmp = tmp(floor(size(tmp,1)/2-l)+1:floor(size(tmp,1)/2+l),floor(s ize(tmp,2)/2-l)+1:floor(size(tmp,2)/2+l));

P_2=P_2+tmp;

end

P_2=P_2/size(a,2);

P_2=rot90(P_2);

imagesc(P_2);

axis on;