数学建模国赛一等奖论文VIP

A题：CT系统参数标定及成像,ParametercalibrationandimagingofCTsystem,报告人：陈淑珂队友：贾梦颖端启航,CONTENT,问题重述,RESTATEMENT,01,CT系统介绍,一种典型的二维CT系统如图1所示，平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。CT系统安装时往往存在误差，从而影响成像质量，因此需要对安装好的CT系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样品进行成像。,问题重述,RESTATEMENT,01,建立模型，解决以下问题,(1)在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息如图2所示，相应的数据文件见附件1，其中每一点的数值反映了该点的吸收强度，这里称为“吸收率”。对应于该模板的接收信息见附件2。请根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。(2)附件3是利用上述CT系统得到的某未知介质的接收信息。利用(1)中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。另外，请具体给出图3所给的10个位置处的吸收率，相应的数据文件见附件4。(3)附件5是利用上述CT系统得到的另一个未知介质的接收信息。利用(1)中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出图3所给的10个位置处的吸收率。(4)分析(1)中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。,图2：几何模版,图3：10个位置点,第一问,FIRSTQUESTION,02,Mathematica,Mathematica是一款科学计算软件，很好地结合了数值和符号计算引擎、图形系统、编程语言、文本系统、和与其他应用程序的高级连接。很多功能在相应领域内处于世界领先地位，它也是使用最广泛的数学软件之一。Mathematica的发布标志着现代科技计算的开始。Mathematica是世界上通用计算系统中最强大的系统。自从1988发布以来，它已经对如何在科技和其它领域运用计算机产生了深刻的影响。,第一问,FIRSTQUESTION,02,探测器单元间的距离,直径为8mm的圆遮挡了29个探测单元，于是可以计算出探测单元之间的距离为0.2759mm，我们还根据当探测器处于模版正上方和正左方时，椭圆的长轴和短轴的长度以及其投影在探测器单元上的个数，分别计算得到的距离为0.2759mm和0.2727mm。,旋转中心在正方形托盘中的位置,旋转中心在探测器上的投影位于第256个探测单元上。以探测器在正上方时最左边的点的下方512个单元的地方建立坐标系，记为坐标系1，那么旋转中心的坐标为(256，256)。,第一问,FIRSTQUESTION,02,旋转中心在正方形托盘中的位置,射线穿过的长度越长，衰减的能量越多，我们用数据画出了CT探测器接收到的衰减能量的分布图如下。我们假设探测器在正上方时左边为第1个单元，右边为第512个单元。当探测器位于模版正左方时（如右上图所示），圆中心与椭圆中心的连线和CT探测器垂直，所以检测到的数值最大的传感器对应的横坐标是235，即表示椭圆中心对应检测器单元序号为235。当探测器处于正下方，即两个物体成像的中心相距最远时（如右下图所示），小圆中心与椭圆中心的连线与CT射线垂直，所以得到了四组椭圆中心对应的探测器单元序数222、223、224、225，取平均值为223.5。,第一问,FIRSTQUESTION,02,旋转中心在正方形托盘中的位置,所以椭圆中心在坐标系1中坐标为（288.5,235），经过坐标变换将旋转中心坐标转换到坐标系2内，得出旋转中心坐标为（-32.5，21）。结合上一部分求出的探测器单元之间的距离，可以计算出旋转中心在距正方形托盘左边41.0345，距上边44.2069处。由图二中给出的各个长度数据我们可以看出旋转中心在椭圆的内部，右图一为任一角度时的旋转中心与椭圆的位置关系，右图二为旋转中心在坐标系2中的位置。,第一问,FIRSTQUESTION,02,CT系统使用的X射线的180个方向,建立几何模型，在选定的直角三角形中利用距离计算，运用反三角函数计算出探测器与水平线之间的角度，然后再转换成射线与水平线的夹角，也就是系统中射线的180个方向。要找出两个圆心在接收器上投影的位置，我们采用的方法是寻找极大值，通过圆心的射线一定比未通过圆心的射线衰减能量更多，所以通过圆心的射线在探测器上的接收数值一定是一个极大值。,第一问,FIRSTQUESTION,02,CT系统使用的X射线的180个方向,探测器能收到两个不同片区的信号，也即是椭圆和圆在某一角度范围内对X光的遮挡没有重叠作用。此时，我们已知椭圆中心和圆心之间的距离，并且可以计算出二者在探测器上对应的距离，由于探测器与X光垂直，所以在图1-5的三角形中，有,其中是探测器与水平线间的夹角，那么X光与水平线夹角为：,第一问,FIRSTQUESTION,02,CT系统使用的X射线的180个方向,探测器只能接收到一个片区的信号，也即是二者对X光的照射有重叠作用。所以我们将这一段单独拿出来，对数据进行筛选，选择那些被椭圆遮挡而未被圆遮挡的地方，进行曲线拟合，补齐被圆遮挡的部分，然后按照上面的几何模型选取点计算角度。由于函数值域为，所以超过的度数得到的结果为负，所以我们人为的做了些矫正，用减去得到的负值得到在内的角度。,第一问,FIRSTQUESTION,02,CT系统使用的X射线的180个方向,二三问,SECONDANDTHIRD,03,拉东变换,逆拉东变换,每个投影的变换对应物体变换中的一条线，我们根据投影得到物体的每条线的话，就能够重构整个物体。,滤波反投影,二三问,SECONDANDTHIRD,03,几何图形的重建,由附件三和附件五我们知道未知介质的接收信息，并且可以用拉动算法画出相对应的正弦图。标定模版、待测介质1、2的正弦图如下,二三问,SECONDANDTHIRD,03,我们利用上面得到的正弦图采用滤波反投影法反投影出模版原本的几何形状。采用Mathematica中逆拉东算法，将数据信息反投影，得到了初步的几何图形。,几何图形的重建,二三问,SECONDANDTHIRD,03,随后我们采用同样的方法制作了标定模版的逆拉东CT图来检验图形位置和角度的正确性。但是得到的图像与题中已给的图像信息不符，位置和角度都发生了偏移。我们把180组的每组数据选取接收信息的极大值（衰减能量最多的时候是射线经过椭圆圆心），将它对应到256的中心位置，做到对旋转中心角度的修正。所以以下是我们对图的进行矫正的过程。首先是旋转中心的矫正，如下图（左边是几何图形，右边是正弦图）：,几何图形的重建,二三问,SECONDANDTHIRD,03,旋转角度进行修正考虑到前面求出的旋转方向中初始位置的角度不是，而是，所以我们把超过的部分截取放到前面，进行旋转中心位置的调整。矫正之后我们绘制了下列图像：,几何图形的重建,二三问,SECONDANDTHIRD,03,修正处理后我们得到二三问中要求的几何图形（左二右三）：,几何图形的重建,二三问,SECONDANDTHIRD,03,对附件2中均匀标定介质接收信息进行逆拉东变换后得到的几何图像不同颜色区域对应不同的测量值，结果为一个的矩阵，记为，记附件一中吸收率的数据矩阵为，标定模板中吸收率为1的部分测量值分布图如下：,吸收率的计算,二三问,SECONDANDTHIRD,03,当时，显示该被测点处无介质；当时，计算对应的平均值为0.133。对附件三和附件五做同样的操作，我们分别得到了和的表格，于是我们如下比式：,吸收率的计算,然后导出了两个的Excel表格，就是待测物体1、2的吸收率的表格，,二三问,SECONDANDTHIRD,03,我们进行相应的坐标转换，由于题中给出的是点与边界之间的距离，所以我们要把长度坐标轴中点的坐标转换成导出的表格中的具体位置。,特定十个点的吸收率,二三问,SECONDANDTHIRD,03,我们进行相应的坐标转换，由于题中给出的是点与边界之间的距离，所以我们要把长度坐标轴中点的坐标转换成导出的表格中的具体位置。,特定十个点的吸收率,二三问,SECONDANDTHIRD,03,然后分别在待测介质1、2吸收率表格中找到相对应的点，就能得到对应点的吸收率。结果如下表,特定十个点的吸收率,第四问,FOURTHQUESTION,04,在计算180次旋转角度的时候，在附近的数据出现了一些偏差。由于CT系统是逆时针旋转，所以角度应该是逐渐增大的，但是在下表这个区域一些数值却变小了，我们在此方面的精确度还不够高。,精度的分析,第四问,FOURTHQUESTION,04,稳定性的分析,由上述模型知根据180个角度512个探测器单元检测到的射线衰减能量，可以得出待测介质的位置、几何形状与吸收率等信息；而当CT系统出现轻微故障，如探测单元或发射射线出现故障导致512个探测器单元接收信息不完整，或某时刻某角度受干扰等影响，导致某一角度探测器未能探测到射线的衰减能量，在这些情况下系统仍能较为准确地计算出各标定参数的能力，称为参数标定的稳定性。为了分析标定模板的稳定性，我们从附件二中随机去除一些数据（某角度某探测器单元若未能正确接收射线，则设该处接收信息为0），将在系统出现轻微故障时（假设由于干扰、故障等原因标定时1接收信息出错）的标定结果，与个数据完整时的结果进行对比，从而得到参数标定的稳定性。由于探测单元之间的距离根据诸多角度圆直径与椭圆长、短轴的投影距离计算，故探测单元间距不会受系统轻微故障影响；由于我们通过特殊点法对旋转中心进行标定，所以当故障数据并未出现在计算旋转中心所需的6个特殊点时，旋转中心的标定结果不会受影响；当系统轻微故障（随机选取1的接收信息为0）时，旋转中心的标定结果出现误差的概率大致为0.6%；关于180个角度的标定，我们假设出现故障的92个数据平均分布在180次探测中，每个角度的计算需要512个数据中的6个特殊点对应接收信息进行计算，则某个标定角度出现错误的概率大致为0.6%；,第四问,FOURTHQUESTION,04,新模版改进精度和稳定性,我们设计的新模板几何特性如下：,我们设计这样一个模版有如下原因：（1）3个圆形，相比于椭圆形，易于加工和批量生产高精度模板，提高精度。（2）3个不同大小的圆形，整体而言既不中心对称也不轴对称，在参数标定的时候能更好的确定探测器方位。相比原来只有两个特殊位置，3个圆在两两完全重叠的时候具有3种特殊状态，在确定CT系统旋转中心在正方形托盘中的位置时提高稳定性。（3）物体体积更小，减少投影相互重叠的概率，利于计算。（4）在计算探测器单元之间的距离时，3个圆的投影都可以用于计算，并且不依赖特殊位置，可以大幅提升可使用数据量，求平均值后更加稳定。（5）3个圆的吸收率不同，可以检测探测器和发射器接受是否正常，而原来只有一种吸收率是无法判断的。,第四问,FOURTHQUESTION,04,新模版改