2017国赛数模论文

论文题目摘要一、问题重述应用ng yong1构造ient in parabolic problemsJ. CT(Computed Tomography)可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。一种典型的二维CT系统如图1所示，平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。CT系统安装时往往存在误差，从而影响成像质量，因此需要对安装好的CT系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样品进行成像。问题一：在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息如图2所示，相应的数据文件见附件1，其中每一点的数值反映了该点的吸收强度，这里称为“吸收率”。对应于该模板的接收信息见附件2。请根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。问题二：附件3是利用上述CT系统得到的某未知介质的接收信息。利用问题一中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。另外，请具体给出图3所给的10个位置处的吸收率，相应的数据文件见附件4。问题三：附件5是利用上述CT系统得到的另一个未知介质的接收信息。利用问题一中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出图3所给的10个位置处的吸收率。问题四：分析问题一中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。 图1. CT系统示意图 图2. 模板示意图（单位：mm） 图3. 10个位置示意图二、问题假设1假设焦点是无限小的, 并且轨道是准确的。2假设射线光谱是单一的。3假设射线是没有宽度的, 且每次输出在能量上没有波动。4假设检测器也是完全没有大小的, 并且可以准确真实的测量投影的值。5假设采样几何为平行束, 并且这些投影无论是通道方向以及投影方向完全是连续的。6图像及投影测量是轴方向上是无限薄的。7假设图像最外侧与光线总是相切。三、符号说明四、问题分析4.1问题一的分析问题一要求根据附件一，二中所给信息求出CT系统的旋转中心在正方形托盘的位置，因此我们需要对所给图形进行几何分析，结合附件二中的吸收率找到旋转中心与所给图形的几何关系，通过几何证明和解析函数求得旋转中心的位置。其次，对附件二进行颜色处理，找到小圆形的吸收率的轨迹，通过小圆形所跨度的探测单元个数，和小圆形的直径长度，求出探测单元之间的距离。再对附件二图表的数据进行逐列分析与椭圆的弦长进行比较，进而求出CT系统使用的X射线的180个方向。4.2问题二的分析对附件二、三进行Radon变换，得出附件三中的被投射的物体，通过几何分析和比较的方法求取该未知介质在正方形托盘中的位置、几何形状和吸收率等信息，并具体给出图3所给的10个位置处的吸收率。4.3问题三的分析对附件五中的数据进行Radon变换，通过与附件二中数据的比较与分析求取该介质的相关信息，并具体给出图3所给的10个位置处的吸收率。4.4问题四的分析根据Radon变换建立的模型，在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性。五、问题一求解针对问题1，为研究问题方便，我们假设CT系统逆时针转动，由已知条件知道，旋转中心内与正方形托盘中心往往不重合，我们以表示旋转中心。512个探测器按如图所示标号，相邻两探测器之间距离。由关于CT系统工作原理可知，对于均质模板，X射线穿过的距离越长则该处吸收率越大，由此可以通过将附件2的Excel数据进行单元格规则修改，数值小的单元格为深色，数值大的单元格为浅色，结合CT系统逆时针转动，分析下图可有结论：（1） 图中较细曲线由X射线透过圆模板得来，另一部分图像由椭圆模板得来；（2） 图中LINE 1 处穿过圆模板的X射线对应的标号靠前，同时椭圆部分宽度最窄，此时接收器应位于第一、二象限内，X射线方向为；LINE 2处穿过圆模板的X射线位于椭圆中间位置，同时椭圆部分宽度最宽，结合CT系统逆时针转动，此时接收器应位于第一、四象限内，X射线方向为。（3） CT系统按固定角度旋转，可以认为相邻两列数据表示CT系统转动最小角度得来经Excel统计，结合以上分析，可得LINE 1 为第150组数据，LINE 2 为第60组数据。显然第180列数据表示X射线的最后状态，此时X射线方向为；第1列数据表示X射线的初始状态，此时X射线方向为。故CT系统共转动了，间隔，由到。然后确定CT系统相邻探测器单元之间的距离。模板中圆的几何参数题目已经给出，附件2中每个数据表示每个探测器的接收到的该点的吸收强度，可认为每个探测器接收一条发出的X射线。由Excel统计可知，X射线经过圆形模板的有29条，由于29个接收器宽度仅有一个模板圆直径大小，可做一下近似处理：认为模板圆两条平行切线出恰好穿过两条X射线，故两相邻探测器之间距离；同理，统计经过椭圆短轴周长的X射线有108条，得两探测器之间距离。为减小误差，取均得。示意图如下：为确定CT旋转中心在正方形托盘中的位置，以椭圆中心为原点，椭圆中心与圆中心连线为x轴建立坐标系，由已知条件，得椭圆方程与圆方程：椭圆方程： 圆形方程：这里我们取X射线发射端进行讨论，由Excel统计，对第150列数据取吸收率最大值，对应第223行，即第223行探测器此时位于y轴；对第60列数据取吸收率最大值，对应第236行，即第233行探测器此时位于x轴。两探测器之间距离。两探测器位置关系如图，作辅助线如图所示： 现证明CT系统的旋转中心必在直线CF上，直线CF方程为：。线段绕旋转中心旋转，则，且，同时显然线段AB上任意一点都满足垂直相等的关系。首先证明直线CF上任一点均能满足上述关系：在直线CF上取点M，连接BM，BM，过点M作DF、EF的垂线，分别交DF、EF于点G、H。若，只需证明，即。由于线段由旋转得来，故易证明，又，且、，得，又，因此。然后证明直线外一点N均不能满足上述关系：连接线段与，与直线CF相交于点，由上述证明可知，则，故与不重合，命题得证。这样还不能确定出旋转中心，因此需再找到这样一条直线与直线CF相交才能得出旋转中心。此时考虑能否再找出线段在其他旋转角度时的位置信息。联立圆锥曲线方程可得两圆锥曲线的外切线的方程： 或 同样的方法可以确定此时第101个探测器位于直线上，该探测器距第223和第235个探测器的距离 可得此时第223个探测器所在直线的方程为：第235个探测器所在直线的方程为：同理，与直线AC相交于点I，与y轴交于点J，直线IJ必过旋转中心P，直线IJ方程方程联立，可得旋转中心为：六、问题二求解6.1建立模型图像投影，就是说将图像在某一方向上做线性积分（或理解为累加求和）。如果将图像看成二维函数f(x, y)，则其投影就是在特定方向上的线性积分，比如f(x, y)在垂直方向上的线性积分就是其在x轴上的投影；f(x, y)在水平方向上的线积分就是其在y轴上的投影。引入Radon变换，Radon变换（拉东变换），就是将数字图像矩阵在某一指定角度射线方向上做投影变换。这就是说可以沿着任意角度theta来做Radon变换。附件二三就是根据图像投影原理所得到的数字图像矩阵，通过Matlab进行Radon变换的就会得到原图像。Radon变换原理N-1维超平面上的积分值是一个维函数的Radon变换，Radon逆变换相应的就是指从这些线积分值反求得到原始被积函数。当N=2时，在二维空间上的函数f(x,y)沿着任意投影射线方向的线积分就为其Radon变换。如图三所示，Radon变换表示为：Radon逆变换的公式为：投影射线的方程可以表示为：其中，代表投影射线的法线和x轴的夹角，r代表投影射线与坐标原点之间的距离。 图三设是X射线对介质在点的衰减系数，表示射线源的强度，即X射线进入介质的强度，简称入射强度，穿透介质后X射线的强度变为，即探测器接受的强度， 和是可以根据附件二中数据可以计算得到，由物理知识，即，其中L是X射线的投射路径。在扫描过程中,沿着任一条直线的线积分都可以通过探测获得，从所有这些积分中可以重构函数。图像重构的问题就是由测量出来的多个方向的投影值来求得物体内部多个部位的衰减系数当 X射线通过介质时，由于介质各个部分密度不同，所以衰减系数就不同于是图像重构所求出的值和断面的各部位应有值的对比，便可得到所求图形形状。 6.2模型求解利用 matlab 的radon变换函数，代入附件一数据对问题一模型进行求解可得介质的几何形状如下由Radon理论可知，经变换和反变换后得到的362*362矩阵的中心点即为旋转中心。程序中矩阵R即为该介质的吸收率信息。那么如何得知正方形托盘上某位置上介质的吸收率呢?矩阵R中的数据是如何与托盘任意位置相对应的呢?不难发现matlab程序得到的图像与托盘中的介质是同一个，并且二者不是在同一坐标系下进行表述的，下面将求解托盘与求解得到的图像的坐标变换关系以及长度对应关系，由此便可以得到托盘上任意位置的介质吸收率。首先将前面求解问题时建立的坐标系进行伽利略变换得到以旋转中心为原点，水平方向为x轴的坐标系，变换关系如下：然后再进行平面直角坐标系的旋转，得到坐标系，变换关系如下：式中，表示坐标系旋转的角度。在矩阵R中，我们可以通过比对接收信息即数值的大小，确定出椭圆与圆的中心所在的行数与列数（x1,y1），(x2,y2)，通过公式，求得角度。其中，圆圆心的行列数可以直接读取为（），椭圆中心的行列数可以由两切线切取椭圆求得。如图4.1读取两切线的列数，分别为，由椭圆对称性知椭圆中心横坐标为。同理可求纵坐标为。解得题目中已给出椭圆中心与圆中心的实际距离为45mm，而矩阵行列式中两点的距离为，由此可得矩阵R行列数距离与实际尺寸的对应关系为。将圆形中心点带入进行验证，模型在较小误差允许范围内是成立的。因为题目中给出问题1固体介质的吸收率为1，而经Radon变换后求得的是相对吸收率的分布矩阵。可以通过比较同一位置（椭圆）的吸收率，得到结论： 介质的吸收率=相对吸收率*2即按坐标变换查得对应的相对吸收率后乘以2得到最终结果。通过 matlab代入附件三数据对问题二模型进行求解可得图像由图像矩阵R可得灰色椭圆中心行列数为（106,210.5），通过进行逆变换可知该介质中心位于正方形托盘中坐标系的。对附件4中点进行标号，得：位置12345m10.0000 34.5000 43.5000 45.0000 48.5000 n18.0000 25.0000 33.0000 75.5000 55.5000 位置678910m50.0000 56.0000 65.5000 79.5000 98.5000 n75.5000 76.5000 37.0000 18.0000 43.5000 由上表变换到坐标系中位置如下表：位置12345x-40.0000 -15.5-6.5-5-1.5y-32-25-1725.55.5位置678910x0615.529.548.5y25.526.5-13-32-6.5再将上表进行伽利略变换得位置12345t-31.1025-6.60252.39753.89757.3975s-37.5591-30.5591-22.559119.9409-0.0591位置678910t8.897514.897524.397538.397557.3975s19.940920.9409-18.5591-37.5591-12.0591然后进行旋转变换，同时注意矩阵R行列数与实际尺寸的比例关系以及原点的行列数（181,181），从而得到上表位置在矩阵中的行列数：进而通过查找矩阵元素，可得这10个位置的吸收率如下表 七、问题三求解通过使用 matlab 的radon函数，代入附件五数据对问题三模型进行求解可得图像由以上分析知，对该未知介质10个位置进行坐标变换后，