2017数学建模A题二等奖论文

1CT系统参数标定及成像摘要目前，如何通过CT机的探测器已知的投影信息还原出物体几何形状以及内部结构，探测已知模块得到的数据推导CT机的参数以及校正，是研究精密仪器的重要方面。本文通过建立数学模型，给出了通过分析表格数据来解决上述所有问题的方法。对于问题一，以正方形左下角为原点建立平面直角坐标系对方向角进行几何抽象，通过介质上任意一点在转动过程中在光屏上的投影点的运动轨迹都是在旋转中心的投影点的上下波动这一几何特性，引入三角函数、直流分量和交流分量并与几何方法结合，建立了已知模板接收信息、吸收率、CT参数、余弦轨迹的傅里叶展开之间关系的综合模型。101coswsinwnyaanbn，tan(tan)kkyxbag求解模型时，将附件的探测器接收数据和带入，利用光路方向角与椭圆中心和圆心之间的距离的几何关系求出180个方向角度，旋转中心坐标（50.2775,56.1050）,相邻探测器间距离为0.2775mm。对于问题二和问题三，本文根据中心切片定理，建立起投影物体的密度与其二维傅里叶变换后的函数的关系模型。在每一个确定角度下，分别使用直接反投影法和滤波反投影法(,)|(cos,sin)|()fxyfRRGR()(,)[(,)|]GFfxyF将变换后的函数作傅立叶逆变换就可以得到原物体的密度函数，并进行插值填充由MATLAB生成吸收率矩阵。问题二的介质几何中心为（43.1641,47.4609），形状大致是椭圆形的嵌套（见图5-9）。问题二的10个点吸收率（附件4坐标点自上到下）分别为：0、0、2.0053、1.9714、1.9613、1.9644、0、0、0、0。图片见本文4.2.2模型求解部分。问题三的10个点吸收率（附件4坐标点自上到下）分别为：0、0、0.6479、0、0、4.7242、、0、10.5331、0、0。图片轮廓形状类似于组织细胞图，图片见本文4.2.3模型求解部分。对于问题四，问题1中的相邻探测器间距离参数，是影响成像精度的重要参数，距离越小精确度越高。旋转中心在托盘上的位置和探测器的个数，是影响成像稳定性的重要参数，增加单位长度的探测器个数，并且调整旋转中心和探测器中心相对位置，使得探测器每转一个角度，接收到信息探测器尽可能在中部（避免下次转过一个角度时，介质的投影信息溢出光屏）。探测器绕旋转中心扫过的方向角度尽量稍微大于180°，使信息采集的更加完备，从而增加了CT工作的稳定性。关键词：余弦轨迹滤波反投影法中心切片定理,傅里叶变换傅里叶级数分析2一、问题重述1.1CT(ComputedTomography)可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。一种典型的二维CT系统如图1所示，平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。CT系统安装时往往存在误差，从而影响成像质量，因此需要对安装好的CT系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样品进行成像。1.2问题的提出1.在正方形托盘上放置两个均匀固体介质组成的标定模板，其中每一点的数值反映了该点的吸收强度，这里称为“吸收率”。请根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。2.利用(1)中得到的标定参数，确定附件3中的未知介质在正方形托盘中的位置、几何形状和吸收率等信息。另外，请具体给出图3所给的10个位置处的吸收率，相应的数据文件见附件4。3.附件5是利用上述CT系统得到的另一个未知介质的接收信息。利用(1)中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出图3所给的10个位置处的吸收率。4.分析(1)中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。二、问题分析2.1问题一分析以正方形托盘左下角为坐标原点建立平面直角坐标系。讨论光线方向角和探测器与x轴夹角的关系，根据投影方向角、旋转中心的坐标、旋转中心到最近的探测器的标号之间的关系，建立光路方程。光路方程被模板的几何图形所截长度正比于所给的吸收率。经过旋转，圆心的运动轨迹随角度的变化是余弦图并且其直流分量k即为旋转中心所对应的探测器的编号，对其曲线拟合和傅3里叶级数分析求出k值。联立两次K探测器接收的光线的直线方程，交点就是旋转中心，投影数据分析求探测器间距、180个方向。2.2问题二分析问题二是一个已知角度下的投影，还原原函数密度的问题。根据中心切片定理：函数密度在某一方向上的投影的一维傅里叶变换函数是原函数密度的二维傅里叶变换函数在同一角度且过原点直线上的值。据此定理结合附件数据求解变换后的傅里叶函数在每一个角度下的频域函数，再取傅立叶逆变换。2.3问题三分析问题三与问题二分析和求解类似，并且都是建立在第一问的已知180次测量角度上的，由于上述处理方法非特殊取值的方法，因此该方法具有普遍适用型，此处不再赘述。2.4问题四的分析模板对稳定性的影响在于模板在正方形托盘位置的分布，由几何概型可知，模板分布越集中于托盘中部，被探测器检测到的概率越大，稳定性越高。因此本文设计出的模板模型是圆心位于正方形托盘重心处的圆。三、模型假设1.忽略空气对X射线的吸收率。2.不考虑光的折射和反射。3.假设所测量的数据均为正常测量的情况，不考虑测量过程的抖动因素。4.假设CT的滤波器是理想的。四、符号说明θ:光线的方向角，光线矢量与y轴正方向夹角,逆时针旋转为增大的正方向β:探测器和x轴正方向夹角i:探测器编（当探测器旋转到平行于正方形上方时，从左向右以0―512编号）K:距离旋转中心最近的探测器编号Q:吸收率d:椭圆中心和圆心之间对应多个探测器的距离l:相邻探测器间距离4y1,y2:包络线位置函数max1y、max2y：椭圆中心和圆心对应探测器的编号n：测量次数（1≤n≤180）w：傅里叶展开函数的基频率（a,b）：旋转中心坐标1f、2f：圆的8mm直径对应28个间距、29个间距的探测器，出现频率(,)fxy：未知物体的密度函数()GR：密度函数在方向上的投影函数()G：密度函数在方向上的投影函数一维傅里叶变换函数(,)F：(,)fxy的二维傅里叶变换函数(,)Wfxy：重建后的密度函数()Hw：滤波器的频域表达式(,)WHF：(,)F滤波后的结果()n：每一次测量角度的弧度制12tt、:存储相邻两次旋转的稳定因子矩阵tl：程序中的矩阵的长度N：值为512，代表512个探测器五、模型的建立与求解5.1.1问题1的模型建立建立平面直角坐标系求解θ和经过K探测器的光路直线方程（分三类讨论）（1）090o5图5-1（090o旋转过程图）由几何关系可知：β=θcosdlβarccosdl经过K探测器的光路直线方程：tan(tan)kkyxbag（2）90180oo图5-2（90180oo旋转过程图）同理可得：90oβ6dsin=lβ=arcsin90dlo经过K探测器的光路直线方程：tan(90)[tan(90)]kkyxbaoog（3）180270oo图5-3（180270oo旋转过程图）同理可得180oβdsin=lβ=arcsin180dlomax1max2||dyylg经过K探测器的光路直线方程：tan(180)[tan(180)]kkyxbaoog求K由几何关系可知，对于单一的几何图形来说，在任一确定角度上，通过椭圆中心的弦和圆的圆心的弦是最长的，即对应的探测器接收数据最大。将附件2的数7据，对每一列数据使用牛顿迭代法找到极大值max1y、max2y极大值对应的探测器编号即为穿过圆心的光线对应的探测器。横轴代表180次测量的转次（暂未知角度），纵轴代表一次测量时，每一个的传感器测量状态（状态：传感器有无吸收率，1≤i≤512）。图5-4圆与椭圆双旋转图曲线拟合方法：由于两条包络是围绕同一旋转中心得到的，所以它们的相位相等，将两条包络线进行取平均得到y。图5-5椭圆和圆连线中点在接收器上的投影位置图8假使探测器不动，模板相对探测器旋转，圆心绕旋转中心的轨迹，分时投影到探测器上的曲线是余弦函数曲线形状，旋转后的结果图如下：（旋转动图见附件）。5-6几何中心点围绕旋转中心旋转产生的余弦图像对它利用傅里叶级数展开可以得到直流分量和余弦分量，直流分量的大小就是K。101coswsinwnyaanbn进行一级傅立叶级数展开：1101coswsinwyaab1120cos2wsin2wyaab1130cos3wsin3wyaab……101coswsinwnyaanbn其中y1、y2、y3已知，a0、a1、b1是未知。用上面求出的θ值对应的每一次位置，计算出每一次经过K探测器的光路直线方程，求两两直线的焦点，求出旋转中心。由平面几何分析可得，每一次经过K探测器的光路直线方程为：tan(tan)kkyxbag联立任意两个角度下的经过K探测器的光路直线方程，计算出多组两直线的交点，利用统计求取均值，所得就是旋转中心点。由于圆在探测器上投影的永远是直径，为定值，因此在每一个角度确定的情况下，通过圆直径所对应的探测器的间距个数，求出l12882928lff95.1.2问题1模型的求解根据数据按列迭代求取极大值，从而确定过椭圆中心和圆心每一个θ下对应的d的大小，最后利用4.1.2中三个情况下的θ数学表达式，求出每一个角度，得到下表（见附件“180个”.xls）：n12345678910角度θ29.498431.574431.574432.892734.165734.786835.398437.180537.758438.8923n11121314151617181920角度θ39.999140.542841.612943.177243.688644.697346.177652.217852.661253.5405n21222324252627282930角度θ53.976354.840955.269856.120956.963657.382058.625559.445459.852761.0649n31323334353637383940角度θ61.865262.263163.448264.621165.396866.167667.315468.075368.831169.9577n41424344454647484950角度θ70.704371.447472.556273.658274.389475.118076.206477.289678.009279.0851n51525354555657585960角度θ80.157081.225581.936082.999684.060789.957989.957989.957989.957989.9579n61626364656667686970角度θ89.957989.957989.957989.957989.957989.957995.586196.292897.354698.0640n71727374757677787980角度θ99.4865100.1999101.2731101.9908103.0709103.7936104.8820105.6106106.7084107.4438n81828384858687888990角度θ108.5526109.6686110.