2017年高教社杯全国大学生数学建模竞赛题目-CT系统参数标定及成像

12018西交数模第一次模拟赛数学建模论文首页选题A队伍编号696班级学号姓名队长钱学森73班2173714071孔令辉队员1钱学森73班2171311098杨峥队员2电气642160400092李佳航2018年6月30日2摘要本题针对一种二维CT获取样品内部结构信息的工作方式及成像原理，意在通过借助已知结构样品进行参数标定，消除系统误差，而后对未知结构的样品进行成像，得出该未知介质的相关信息。问题一中，要求通过标定模板的相关，确定CT系统的旋转中心、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。本文利用标定模板几何参数，通过条数、以及探测器单元间距相等等信息，首先计算出探测器单元之间的距离为0.2778mm。根据X射线与标定模板的几何特性，以椭圆短轴所在直线为x轴，长轴所在直线为y轴建立直角坐标系，得旋转中心在所建立的坐标系中的坐标为7060.5,0396.9。最后通过近似确定X射线180个方向的旋转角度具有高度线性相关性，根据部分确定数据拟合出整体旋转角度。前十五组旋转角度为：29.649331.002731.557232.647633.679534.647835.648336.649037.647538.648639.648840.649141.648742.649743.6492问题二中，要求通过已求得的标定参数，确定未知介质在正方形托盘中的位置、几何形状和吸收率等信息，并具体给出所要求十个位置的吸收率。本文依据CT层析成像原理，利用逆拉东变换作出重构图像，并利用Excel中数据分布计算出原位置介质的相关性质。所求十个位置的吸收率为：序号12345吸收率0.01850.9701-0.00011.17031.0287序号678910吸收率1.44251.2815-0.0074-0.00250.0194问题三中，与问题二计算方法相似得十个位置的吸收率为：序号12345吸收率0.07362.53057.13270.03320.8015序号678910吸收率3.14875.57140.05308.13620.0272问题四中，要求对问题一中的标定模型进行改进以减小误差并增加稳定性，本文利用在问题一求解过程中遇到的问题进行思考，首先适当增大模板减小偶然误差，其次做出投射图像后应容易找到极值，并且图像应有一定的对称性；图像扫描后应尽量少地得出重复数据。关键词：CT层析成像Radon变换与逆变换吸收率MATLAB算法31.问题重述1.1问题背景CT(ComputedTomography)可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。1.2相关信息本题介绍了一种二维CT系统，平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。并且，为消除安装误差，需要对安装好的CT系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样品进行成像。题目附件中提供了标定模板的几何信息，接受信息，待测介质的接收信息以及图3所给位置的相应数据。1.3需要解决的问题1.3.1问题一在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息如图2所示，请根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。1.3.2问题二附件3是利用上述CT系统得到的某未知介质的接收信息。利用第一问中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。另外，请具体给出图3所给的10个位置处的吸收率。1.3.3问题三利用给出上述CT系统得到的另一个未知介质的接收信息。利用第一问中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出图3所给的10个位置处的吸收率。1.3.4问题四分析第一问中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。图1.CT系统示意图图2.模板示意图（单位：mm）图3.10个位置示意图2.问题分析42.1问题一分析结合图2和附件1表中数据，可以首先计算出CT系统探测器个数和模板长度度量的比值，运用程序1.1可以得出模板的几何形状如图2.1.1，可大致认为它是对称的。对于附件2，由于对180个方向尚无清晰地认识，首先用同样的方式做出数据分布图2.1.2，观察到图像比较平滑，因此认为按表格的顺序180个方向是相邻较密、不错位的。（图中有色区域表示该点有吸收率，蓝色部分表示吸收率大于100）图2.1.1附件1的数据分布图图2.1.2附件2的数据分布图在图2.1.2中易观查到有一红色条形图案，这是在不同的方向扫描到圆形时留下的，可以此为突破口首先求出探测器单元之间的距离。求出探测器单元之间的距离与增益比率之后，可以根据几何关系，自己设立坐标系并通过数学运算计算出旋转中心。最后在建立的坐标系内，将待求解的180个旋转方向转换成X射线的斜率进行数学运算。2.2问题二分析题目给出了未知介质的接收信息，要求出介质的相关信息，可以搜索相应的数学模型，将附件中给出的按照射线条数与旋转角度列成的表格，一一对应为相应坐标点的吸收情况，从而根据各坐标点的不同性质，还原回该未知介质的几何信息与吸收率等信息。2.3问题三分析问题三与问题二类似，可以大致看出数据分布更具有一般性，不容易描述出未知介质的相关信息，可以通过图形大致描绘出介质相关信息。2.4问题四分析问题四要求分析题目所给的二维CT系统，设计新的标定模板以提高原系统的精确度与稳定性。可以搜索相关资料，根据第一问的求解思路与求解过程，以规避求解过程中因标定模板自身性质而出现的误差为原则进行思路拓展，以设计高精度与稳定性的标定模板。3.模型假设假设在射线经过介质时能量只损失在介质中，及不考虑衍射等现象；假设附件中所给出数据是正确的、可以直接利用的；5假设旋转中心在相邻两条射线的中间直线上；4.模型的建立与求解4.1问题一模型的建立与求解4.1.1建立模型从题目中可以得出，由于x射线之间得间距相等，不管x射线怎么旋转，穿过托盘上圆的x射线条数应该是大致相等的。可将穿过圆形标定模板的X射线与模板建立模型示意图如下：图4.1.1.1X射线穿过圆形标定模板示意图如图4.1.1.1所示可以设穿过圆的在圆心两旁的射线被圆所截的距离分别为和。求解旋转中心时，以椭圆形标定模板短轴所在直线为X轴，长轴所在直线为Y轴建立直角坐标系，将托盘进行划分。由附件二的数据分布图的不对称性可知，旋转中心应在对称轴某一侧，大致确定旋转中心方位后，根据旋转中心两侧探测器个数不变且同一探测器接收的与距旋转中心的距离不变具体确定旋转中心的位置。求解X射线旋转角度时，设穿过椭圆的最边缘的射线到椭圆中心的距离为R，取椭圆中心为原点，模板的对称中心为x轴建立平面直角坐标系，X射线所在直线的斜率为k。4.1.2符号说明4.1.3模型求解求探测器单元间距离：d相邻探测器之间的距离N对应于圆形模板的射线条数或探测器个数μ标定模板的吸收率λ处理数据时的增益率φ圆形模板的直径6通过MATLAB编程求出计算出的条形带平均涉及探测器个数为8023.28N，沿直径方向上的平均吸收率1766.14；由图二圆形模板直径mm8，计算出探测器之间的平均距离：mm2778.0mm8023.288Nd探测系统平均增益率：7721.181766.14求旋转中心：分析图2.1.2可知，由于蓝色区域仅分布在后部分角度范围内，因此估计旋转中心的位置在对称轴的某一侧，红色区域分为两部分时表示该角度下由两部分射线分别照射经过椭圆形和圆形；为了讨论方便，现对正方形托盘做出如下划分：图4.1.3.1圆盘划分简单分析可知，如果旋转点在I区域，则沿180个方向照射后不会出现射线分成两部分的情况；在II或III或IV区域，当吸收率出现最大值的时候射线也被分成两部分，而不是像图1.2那样成为一部分，也排除；综合各因素可判断旋转中心应该在V区域。因为发射-接收系统逆时针旋转，且图2.1.2中两部分红色区域，间距缩小，说明旋转之后在垂直于发射-接收方向上二者的距离是缩短的，从而确定旋转中心在短轴的上半侧，先运用Excel对附件二第1列数据作图（即画出第一个方向上的扫描图像）如图4.1.3.2，发现此时得到两部分图像，说明两模板之间有一部分射线直接被探测器接收，直到第14列数据两部分图像结合在了一起，如图4.1.3.3。7图4.1.3.2第一列数据分布图4.1.3.3第14列数据分布结合托盘的几何特征，在垂直于对称轴方向上应该会出现最大的吸收率，利用MATLAB求得出现最大吸收率的方向为第151个方向，在此方向结合增益率得出的模板长度为80mm79.9884进一步验证了结果；平行于对称轴方向上最大的吸收率出现在沿对称轴的直线上，计算得出为第61个方向，此方向的模板长度8mm3037.9988，也验证了结果。根据极近似水平方向为第61方向、最大吸收率出现在235号探测器，极近似竖直方向为151方向、最大吸收率出现在223号探测器，设旋转中心到竖直轴、水平轴的距离分别为x,y。为求解还需要另一个方向的等量关系，选取椭圆和圆的一条外公切线的方向，经计算得出经椭圆与圆公切线所在直线的X射线为第372号射线。则在所建立坐标系内，切线方程为：8255.418148.0xy进而列出二元方程组：dyxydxy)235372(8148.018148.08225.41)223235(2解之得:7060.50396.9yx因此旋转中心在所建立的坐标系中的坐标为7060.5,0396.9。求180个旋转方向：由于x射线的发生装置是连续旋转的，所以在512条射线中穿过模板的射线条数应该是连续变化的，所以用MATLAB编写程序，计算每次旋穿过模板的x射线条数，并画出图像如图4.1.3.4。8图4.1.3.4穿过模板射线条数随旋转次数变化曲线图中横坐标为旋转次数，纵坐标为穿过模板的射线条数。从图中可以看出附表2中的数据是按照射线发射装置旋转的顺序依次给出的，而且可以看出，180次旋转后装置共旋转了180度，每次旋转的角度近乎相等。考虑到在前50次旋转中，穿过整个装置的x射线条数与穿过椭圆的条数相等。设穿过椭圆的最边缘的射线到椭圆中心的距离为R，取椭圆中心为原点，模板的对称中心为x轴建立平面直角坐标系，由于射线可以近似看成与椭圆和圆都相切，可以得到以下方程组：{𝑥𝑦𝑅𝑦𝑘𝑥𝑚{𝑥15𝑦401𝑦𝑘𝑥𝑚化简以上方程组可得：k√40−𝑅𝑅−15√40−𝑝𝑠𝑝𝑠−15利用MATLAB编程得到图像如图4.1.3.5所示9图4.1.3.5前50次旋转角度变化曲线用同样方式得出最后面25组数据的图像为：图4.1.3.6后25次旋转角度变化曲线从图像中可以看出每次旋转角度近似为1度。为了更加精确地计算旋转角度以便据此得到估算其他角度地依据，本文设计了另外一种算法。从附表二中可以看出，前15组数据中穿过椭圆的射线与穿过圆的射线没有交叉，所以数据中换算出来的最大吸收距离就是近似穿过椭圆中心的射线被椭圆所截的距离，弦长公式为：10{𝑥15𝑦401𝑦𝑘𝑥𝑚且这个距离由直线的斜率k唯一确定，已知k、m时可解得：l1200√𝑘115𝑘40√15𝑘−𝑚40所以用MATLAB编程计算得出了较为精确的前15次旋转角度的结果如下：29.649331.002731.557232.647633.679534.647835.648336.649037.647538.648639.