数模国赛2017A题原创优秀论文

三、模型假设
1.假设CT光源的旋转中心在探测器的中垂线上。

2.假设X光不会发生衍射等其他影响吸收强度的现象。

四、符号说明
五、模型建立与求解
1.问题一
1.1．建立坐标系
椭圆方程较为复杂，为方便分析，选择在椭圆中心建立直角坐标系，可得模板椭圆和圆的方程为：
1.2. 增益的确定
1.2.1 的模型
查阅资料可知X光吸收强度与其穿过的介质长度和密度有关，令模板的密度函数为，可得
由于椭圆和圆模板均为均匀介质，可认为为常数，可得
可知X光吸收强度和其穿过的介质长度呈正比，令增益，即可得
1.2.2 的计算
选取中非0数据最多的六列数据，可以有效减小系统误差。

取每一列数据数值最大的几个值，其表示椭圆短轴和圆直径吸收衰减后的X射线能量经增益处理的量值，取六个方向平均值，对应为38；同理选取中非0数据最少的六列数据，此时探测器位于平行于x 轴的位置，两段不为0 数据中的最大值分别表示椭圆长半轴和圆直径吸收衰减后的射线能量增益后的量值，取三个方向平均值分别得,
对应的,为80 和8。

对这三组数据用excel进行最小二乘法拟合，得到μ=1.7713。

过程如图所示：
1.3 探测器间距离确定
通过附件2，可知中每一列非0数据的个数，即为X光源截得相应弦长，对应的探测器的个数。

则当探测器平行于y轴时，探测器的个数最多；平行于x轴时，探测器的个数最少。

将附件2数据，用Matlab可视化,如图
可确定在，有最少个数探测器;,有最多个
数探测器。

得到当时，之间，有
个探测器；当时，之间，有个探
测器。

最终可算出
取均值得
1.4 旋转中心的确定
当时，设第行, 使得取到最大值；当时，设第行, 使得取到最大值,。

显然当时，其X射线路径通过原点。

其截得模板的长度分别为椭圆长轴和短轴。

有1.3图像可知：
将在这两个位置将椭圆中心即坐标系原点与旋转中心之间的探测器单元数目差值分别确定，找到模板和探测器系统的相对位置，代入d 值，分别求得纵坐标和横坐标。

则
1.5 180次旋转角度的确定
1.5.1 180次旋转角度的模型
在发射系统，绕(,)(-9.2233,6.0182)旋转过程中，可用表示出发射系统所对应的直线的方程
该直线截的模板长度之和有3种情形。

用3种情形下,分别对应的,的关系
发射系统与探测器平行且长度相等，则探测器之间的距离之和，即为发射系统与所截得模板的长度。

n为对应的探测器的数量，即为每列中非0数据的个数
1.5.2 180次旋转角度的求解
根据上述模型，采用excel统计相关数据，Matlab程序求解出180个角度，Matlab程序见附表
2.问题二
2.1 未知介质几何形状的确定
2.1.1 未知形状介质几何形状的模型
由1.2可知，吸收强度与介质密度函数呈现线积分关系
经查阅相关资料，可知此关系满足radon变换。

用直线到原点的距离，直线倾斜角为参数，直线方程为
用狄拉克函数对积分区域进行限制
则可用
显示出X光将介质每一条直线的信息投影到空间内。

空间内，每一个点的数值大小反应吸收率的强弱，反应了介质的几何分布。

2.1.2 未知形状介质几何形状的计算
对于空间内的每一个点，可用radon逆变换求得其在原像空间中的直线信息
可用Matlab的iradon函数进行操作，Matlab代码见附表。

几何形状如图所示。

其中颜色越浅的部分，吸收率越强，吸收率可视图如下
2.2 未知介质位置确定
由iradon函数得到像素矩阵，其对应在托盘上，可
得平均每个像素所占长度
以正方形托盘中心为坐标原点作图
2.3 10个点的吸收率的确定
观察到根据附件2的模板接收器信息，iradon函数返回值与附件1模板真实信息之间，存在一个比值m，定义为
用excel进行最小二乘法拟合得到m = 2.0751，即接收信息为原始图像的真实“吸收率”为函数返回值的2倍。

通过附件4，确定10个点的坐标，查询iradon
函数的返回矩阵后，得到吸收率，如下图
3.问题三
3.1 未知介质相关信息的确定
3.1.1.未知介质的几何分布
(1)直接反投影法的局限性
与问题二相同，使用radon逆变换即可求得未知介质的几何分布。

通过附件5的数据，用Matlab直接反投影，发现图像具有大量阴影，如下图。

程序见附表。

直接反投影法把每次测得的投影数据“原路”反投影到投影线的各个像素上。

即指定投影线上所有各点的值等于所测得的投影值。

如图，被探查介质只是在原点位置上的一个致密点。

在一个角度上测量到其投影值时，就把这个值赋给投影线上的所有点。

于是，从不同角度进行反投影后的重建图像是由以原点为中心的一系列辐射线。

显然，原点位置上的分布密度最高；愈往四周，密度愈低。

直接反投影只是粗略地把图像恢复出来，但除了密度最大的中心点，四周出现了逐渐变浅的晕状阴影。

(2)采用中值滤波及剔除低值后---反投影法
1.中值滤波原理简述：
图像滤波,即在尽量保留图像细节特征的条件下对目标图像的噪声进行抑制，是图像处理中不可缺少的操作，其处理效果的好坏将直接影响到后续图像处理和分析的有效性和可靠性。

中值滤波将局部像素矩阵的中心像素改为全部点的像素的中值。

经过中值滤波后，反投影图像，如下图所示，其实现代码，见附表
可发现，图像上依然有较多“毛刺”，使得图像依然不够清晰。

经判断是由X 光系统扫描到吸收率较小的杂质而产生的噪音，故通过程序，筛除这部分杂
音，得到较为清晰的图像，实现代码见附表
同2.2，以矩阵中心为坐标原点，建立坐标系
3.1.3 10个点吸收率的确定
个点的坐标，比对数据矩阵，可得10个点的吸收
同2.3，先确定10
率
4.问题四
4.1 误差的确定
(1)X射线旋转方向引起的误差
1.精度分析：
在该CT 系统中，模板或待检测物质固定，X 射线绕着旋转中心旋转180 次，由于系统本身在每次旋转时存在误差，累加CT 系统每次旋转对像素点的
处理所带来的误差，我们假设在每次旋转时产生的误差为，那么，在180次旋转中带来的总误差为180，这大大增加了重建图像时的误差。

2.稳定性分析:
随着使用时间的增加，CT 系统会出现一些老化、磨损，测量精度不够等问题，进而影响系统的稳定性。

(2)探测器单元间距引起误差
1.精度分析：
CT系统实际上只对每一个像素点进行操作，而非连续性操作操。

探测器间距
为d=0.2767mm，而X量级，可以认为X 射线没有宽度。

如果假设间距 d 足够小（几乎可以认为每个探测器单元是连续变化的），那么所得的接收信息可以认为是包含了待检测物质的所有信息。

而在实际中，由于间距d 与X 射线量级之间的差别，产生了系统误差，这样的误差会让重建原始图像时，所得图像与真实图像之间存在一定差距。

2.稳定性分析：
随着旋转次数的增加，CT 系统的探测器单元在旋转时，由于上述的探测器
单元间的微小偏差，会使得系统稳定性下降。

(3)旋转中心引起的误差
1.精度分析：
在问题一中，以椭圆中心为坐标原点建立直角坐标系，当X 射线分别与x轴平行和与y轴平行时，利用探测器中心与模板椭圆中心所产生的差值来确定旋转中心。

长轴对应表格的每一列中非0 的数据个数最多的那N列。

而由于CT 系统在成像的过程中是对像素点进行操作，与实际距离存在一定偏差，即观察附录数据，取了非0 数据个数最多的3列(理论上只有 1 或 2 列)，这使得求解旋转中心时对数据的选择产生误差，进而导致重建图像时的误差。

(2)稳定性分析：
由于所求得的旋转中心存在误差，且X 射线平行入射，使得每次计算的旋转
角度都不相同，从而导致稳定性降低。

4.2 创建新模板
4.2.1 新模板示意图
4.2.2新模板的评价
与椭圆模板相比，新模板有以下优点：
(1)旋转中心：几何中心与旋转中心重合越高，重建效果越好。

选择图像几何中
心为旋转中心。

有效避免了信息的偶然性造成的误差。

同时，在最大的N列，该模板较好的减少了这几列非0数据的个数，从而提高了旋转中心的精度。

(2)探测器单元间距：椭圆模板在求间距时采用长轴对应接收数据最大值。

而在
该模板上采用“三边平分法”，即每个底边对应的探测器个数求平均。

从而避免了单一数据造成的偶然性，进而提高精度。

(3)旋转角度：椭圆模板在求旋转角度时，利用椭圆从与坐标系水平到与坐标
系垂直，求出所旋转的次数，进而求出每次旋转的角度。

而该在求旋转角度时，“最长三法”，意味着在探测单元经过最长边时，经过正三级形其中某一边长，此时数据显示为最长，即0 数据个数最少，当探测单元转过60时，经过相邻的一边，此时数据显示再次最长，两次之间旋转的次数便于得到，正三角形每个角的角度均为60，这样便求得每次旋转角。

为了增大精度，减小误差。