1、CT 系统参数标定及成像摘要本文根据给定的 CT 标定板设计了一套平面 CT 标定的算法，具有较好的精度和稳定性，并建立了通过投影值反求吸收率的重构模型与算法，经过去噪后可以得到较好的吸收率图像。同时设计了一套新的 CT 标定板，通过充分利用几何特征来提高标定算法的性能。第一问中，首先通过一部分投影值与圆弦长的数据进行最小二乘拟合，确定了二者之间的线性关系，初步计算得增益系数 值为 1.7727。再进行整体的几何模型建模，建立坐标系，构造探测器接收信息与旋转中心、探测器位置，及探测器单元之间的距离的函数关系，基于最小二乘的思想进行非线性优化，为了增加模型的精度和稳定性，使用迭代优化的方法对参数
2、进行求解。求得旋转中心 R0的坐标为(-9.2696mm,6.2738mm)，传感器间距为 d=0.2768mm。根据所建坐标系，CT 系统 X 射线逆时针旋转的始值为-60.3465，终值 118.6439。第二问中，首先根据问题一中的标定数据，对附件 3 中的数据进行预处理，将其变换为旋转中心在正方形托盘正中心的数据。再分别建立连续、离散两种 CT 反投影重建模型。连续模型中，利用傅里叶中心切片定理，设计滤波反投影算法（FBP），先将投影数据进行傅里叶变换，滤波后逆傅里叶变换，将所得的值在反投影平面累加，实现吸收率图像重构；离散模型中使用代数迭代法（ART）对网格化的离散数据直接建模，用矩
3、阵迭代的方法对像素网格进行处理，实现图像的重构。最终我们基于两种不同的算法，分别对模型进行重建，给出了题中所要求的十个点处的吸收强度。第三问中，针对更加复杂的噪声，以及更复杂的数据，我们对噪声进行分析。用 SPSS软件对可能噪声值进行 Kolmogorov-Smirnov 分布检验，确定其分布，求出其噪声为均值0.1549 的均匀分布。在此基础上，使用同第二问的方法进行反变换、去噪，得到最终的重建图像，并给出十个指定点处的吸收强度。第四问中，我们对问题一中的标定算法进行稳定性和精度分析，并设计了一套新的标定模板。首先利用仿真的方法，通过设定所需标定参数，生成探测器接收信息数据，再通过数据反求其
4、标定参数。求得在无噪声的情况下，标定参数误差较小，标定精度较高；在数据有人工噪声的情况下，各个参数的相对误差在 1%以下，稳定性良好。同时，对于问题一中迭代算法选取的不同初值，模型均能求解出其标定参数，且误差较小，证明我们的模型对初值选取不敏感。最后，我们设计了一个具有显著几何特征的图形作为标定模板。通过充分利用几何特性，可以求得各个参数的初步估计值，以估计值作为初始值进行最小二乘的非线性优化，得到了良好的标定效果。关键词：最小二乘拟合；迭代优化；滤波反投影算法；代数迭代法；1获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众
5、号：数模加油站】国赛交流分享群：5444576571问题重述CT 系统是一种利用样品对射线的能量吸收特性对样品进行断层成像的技术，在不破坏样品的情况下获取样品内部的结构信息。问题中使用的是一种二维 CT 系统，探测器平面发出平行入射的 X 射线，探测器绕某固定的旋转中心逆时针旋转 180 次，可以获得 180组接收信息，每组信息有 512 个等距单元的数据。我们需要建立数学模型和算法，解决以下问题：1.根据题目给出的标定模板，以及标定模板的接收信息，对 CT 系统的旋转中心的位置、探测器单元之间的距离以及射线的 180 个方向进行确定。2.根据给出的某未知介质的接收信息，利用(1)中得到的标定
6、参数，对介质的位置、几何形状、吸收率进行确定。且给出指定位置的吸收率信息。3.根据给出的另一未知介质的接收信息，给出该未知介质的相关信息，以及指定未知的吸收率信息。4.对问题(1)中的参数进行精度和稳定性分析，并自行设计新模板，建立对应的标定模型，以改进参数标定。2问题分析第一问中，首先需要确定吸收率与衰减强度之间的数学关系。考虑使用一部分投影值与圆弦长的数据进行分析，确定二者之间的函数关系。再进行整体的几何模型建模，建立坐标系，根据旋转中心、探测器与坐标轴的夹角来建立关于弦长的投影模型。考虑到待优化参数较多，可迭代优化参数，再通过改进的最小二乘法，进行非线性优化，得到标定值。第二问，首先根据
7、问题 1 中求得的标定参数，对附件 3 中的数据进行预处理，使之旋转中心变换到正方形托盘中心。再进行 CT 反投影重建建模，考虑连续和离散两种建模求解方式。连续模型中利用傅里叶中心切片定理，设计滤波反投影应算法（FBP）进行重建。离散模型中使用代数迭代法（ART）对像素网格直接进行求解。再对吸收率图像进行降噪处理，考虑比较多种降噪算法表现，得到重建图像。第三问与第二问相比，数据加入了噪声，且图形更加复杂，对算法稳定性要求更高。首先进行去噪，用 SPSS 软件对可能噪声值进行 Kolmogorov-Smirnov 均匀分布检验，确定其分布，再按照同第二问的方法进行反变换、去噪，得到最终的重建图像
8、。第四问首先用模拟仿真的方法对第一问的标定模型进行检验，研究其精度；再考察其对于噪声和迭代初值的稳定性。最后考虑建立新模板并利用其几何特征估计参数并作为迭代初值，提升标定算法速度与精度，并进行稳定性验证。3基本假设1.载物托盘平面始终垂直于 CT 旋转轴，发射-接收系统平行于托盘平面2.X 射线不与样品发生衍射，散射，不考虑硬化效应2获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576573.模板加工精度足够高，不考虑模板几何形状的误差4问题一的求解在问题一中，我们需要根据已知几何信
9、息的标定模板对 CT 系统进行标定。首先我们的对标定模板进行建模。在这一问题中，我们需要得到的量有：旋转中心在正方形托盘中的位置、探测器单元之间的距离、以及 CT 系统使用的 X 射线的 180 个方向。4.1理论基础问题中，使用的 CT 系统的原理是使用 X 射线照射在样品上，物体吸收部分射线的能量，使得射线强度产生衰减，射线的衰减呈指数变化。3 那么对于长度为 l 的均匀同质物体，吸收率为，则探测器上对应位置得到的接收信息 D 应满足：D=f(l)(1)而对于平面内的不均匀介质，射线沿某一路径后得到的探测器的值应为：D=f(l)dl)(2)因此，在进一步分析之前，我们需要先确定经过增益等处
10、理之后得到的接收信息与其吸收量之间的关系。图 1:附件 1 几何形状图 2:附件 2 信息观察附件 2 中接收信息图可以发现，图像可以很明显的看成是两部分叠加而成。结合标定模板集合信息可知，这两部分分别为椭圆部分的投影与圆部分的投影。由于圆具有各方向投影均相等的良好性质，因此我们截取其中一部分数据进行试探性研究。选择第一个测量角度的属于圆的投影的非零数据，得到一组探测信息 D(i)。对于探测点 i，其发出的射线与圆心的距离为 di=id+d0，其中 d 为探测器上的射线间距，d0为探测器位置相3获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模
11、相关资料关注【公众号：数模加油站】国赛交流分享群：544457657对偏置。那么第 i 条射线与圆相交的弦长为 2R2 d2i，其中 R 为圆半径。我们推测探测器的接收数据与 l 成一次函数关系，即D(i)=2R2(id+d0)2+c(3)使用最小二乘法对曲线进行拟合，求出各个参数的值为：=1.7724,d=0.2768,d0=4.0688,c=0.0000图 3:拟合圆上的数据即是说对于标定模板，其吸收强度 D=l，对于参数，在后面的模型中我们会进一步进行求解，此处仅是为了寻找模型之间的关系进行初步的分析。4.2模型建立以正方形托盘的中心为坐标原点，椭圆中心与圆中心的连线方向为 x 轴，过坐
12、标原点垂直于 x 轴方向为 y 轴，建立平面直角坐标系。在这一坐标系中，设 CT 系统的旋转中心坐标为 R(x0,y0)，探测器平面与 x 轴的夹角为，探测器中心与旋转中心在探测器平面上的投影的距离为 d0。图 4:标准形态图 5:圆的弦长4获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576574.3理论公式由于这个问题中，标定模板是均匀介质，所有地方的吸收率都为 1，那么其探测器获得的数值即与穿过的长度有关。那么在这一问题中，我们转化为某一方向的直线穿过标定模板的长度的求解。标定
13、模板中，椭圆与圆的方程分别为：x2m2+y2n2=1,m=15,n=40;(x 45)2+y2=42为了对问题进行化简，我们设定一个坐标原点在探测器上的投影位于探测器中心的状态，此时对于探测器上与探测器中心相距为 d 的一条射线，可以对其在椭圆中交出的弦长进行求解。这条射线的方程可写出为：xcos+y sin=d与椭圆的方程进行联立求解，可解出该直线与椭圆相交的弦长：p=2mnr2 d2r2,这里r2=m2cos2+n2sin2对于圆上的部分，设圆心坐标为(G,0)，圆半径为 r0，其中 G=45,r0=4，那么容易求出其弦长表达式为：p1=2r20(Gcos d)2同时，由于在整个旋转过程中
14、，直线与椭圆或圆不一定有交点，那么这种情况下，我们对整个角度范围进行整合，可以得到总的投影长度为：pt=2mnmax(0,m2cos2+n2sin2 d2)m2cos2+n2sin2+2max(0,r20(Gcos d)2)考虑到实际情况中，探测器的中心与坐标原点有一定的偏移，那么在上式的基础上，我们需要对投影长度进行修正。易证，在探测器能探测到物体的前提下，沿垂直于探测器方向的探测器位置变化对投影长度没有影响。即在探测器位置发生变化时，只需要考虑沿探测器方向的位置变化。图 6:计算弦长5获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资
15、料关注【公众号：数模加油站】国赛交流分享群：544457657那么对于旋转中心坐标 R(x0,y0)，探测器中心与旋转中心在探测器平面上的投影的有向距离 d0，对于探测器上一条射线，其与探测器中心的距离为 d，那么容易求出，该点与坐标原点在探测器上的投影的距离 d的值为;d=x0cos+y0sin+d0+d即可转化到刚才求出的标准形式的公式上去。由于在实际数据中，探测器是有 512 个 X 射线组成的，对于第 i 个发射器，其与发射器中心的距离也可以简单的表示出来，为：d=(i 256.5)d这里距离的负号表示探测器中心与该点的连线沿探测器与 x 轴夹角 的反向。则d=x0cos+y0sin+
16、d0+(i 256.5)d标定模板的吸收强度处处相等且都为 1，为计算吸收强度与探测器探测到的数值的关系，我们引入增益系数，表示投影长度与最终结果的关系。即：V=pt综合上述各式，所以对于探测角度为 的探测器，其上第 i 条 X 射线探测所得的数值的计算公式为：V=2mnmax(0,m2cos2+n2sin2 (x0cos+y0sin+d0+(i 256.5)d)2)m2cos2+n2sin2+2max(0,r20(Gcos (x0cos+y0sin+d0+(i 256.5)d)2)问题即转化为，根据附件 2 给出的数据，对该公式中的参数进行求解。4.4模型求解为了确定模型中的各项参数信息，我
17、们需要对模型参数进行求解，使得理论值与实际测量值误差最小，我们使用均方根误差来作为衡量指标。即min=180j=1512i=1(V(i,j)D(i,j)2n 1这里 D(i,j)指附件 2 中，第 j 组探测数据中第 i 个探测器的测量值，V(i,j)指第 j 组探测数据中第 i 个探测器的计算值。使用 MATLAB 的非线性拟合函数工具进行模型求解，由于这是一个数据量较大的非线性优化模型，直接寻找到最优解比较困难，函数在进行求解时会从给定的初值开始迭代，初值的选取对计算结果有直接影响。因此，我们使用以下步骤对参数进行迭代求解，获取更准确的参数值：6获取更多数学建模相关资料关注【公众号：数模加
18、油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576571.选取参数 d0,x0,y0,d,i,i=1,.,180 的初值，这里参数的初值可以随意选取，我们设置为 d0=x0=y0=0,=1.7724（由前文粗计算得），使用 512 180 组数据，对所有参数进行拟合，获得第一次求解的参数结果。2.对第一次求解得到的 1i,i=1,.,180，使用局部加权回归散点平滑法进行平滑处理，接着将 i作为已知参数，使用 512 180 组数据，以第一次求解的结果作为初始点，对参数 d0,x0,y0,d 进行求解，得到第二次参数的求解结
19、果。3.使用第二次求解得到的参数 d0,x0,y0,d 作为已知参数，对于 180 组，每组 512 个数据，以第(2)步中得到的角度作为初值，分别对每组数据的角度参数进行求解，得到第(3)步的求解结果。4.使用第(2)步与第(3)步的结果作为初始值，再次使用 512 180 组数据，对所有参数进行拟合，得到最终结果。4.5求解结果使用上述算法对模型进行求解，得到各个参数如下：=1.7727,x0=9.2696,y0=6.2738,d0=0.0000,d=0.2768同时可以求出 180 个测量角度的数据，角度数据见附录中表9，部分数据见表1，由于我们求的是探测器的角位置，因此这里将所得数据转
20、换为题目要求的 X 射线的角度位置：表 1:X 射线角度信息角度序号1239091179180角度（度）-60.3465-58.9927-58.439228.638229.6437117.6509118.64394.6模型评价在问题(1)的模型中，我们使用均方根误差来进行衡量，第一次求解与最终结果的均方根误差如表2，从表中可以看出，我们的算法对原始数据的拟合程度较好，均方根误差很小，因此可以认为我们得到的参数信息是较为准确的，模型的可靠性较好。表 2:算法均方根误差算法步骤第一次求解最终结果RMSE4.41610.01485问题二的求解在问题 2 中，需要根据某未知介质的接受信息，利用问题 1
21、 中求出的参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率。第一题中分析了投影的过程，此问需7获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657要利用给定的投影数据进行逆过程。下文首先利用标定数据对投影数据进行了预处理，并分别用连续模型的 FBP（滤波反投影）算法和离散模型的 ART（代数迭代法）进行反投影重建，并比较了数种滤波降噪算法的效果。5.1数据预处理由问题一可知，由于安装过程的误差，CT 系统的旋转中心并不在正方形的中心，由前文推导公式可知，探测器上的第 i
22、 个传感器与坐标原点的沿探测器方向的距离为：d=x0cos+y0sin+d0+(i 256.5)d对于任意角度的探测器，我们都可以设置一个中心与原点重合的辅助探测器，将探测器上的数据转化到辅助探测器上，再对问题进行求解。即：d=(i 256.5)d=x0cos+y0sin+d0+(i 256.5)d则对于原始数据 y=data(i,j)，可以得到新的转化后的数据datapre(i,j)=data(i(i),j)=data(i(x0cos+y0sin+d0)/d,j)图 7:数据预处理图 8:连续模型5.2连续模型首先，我们把第一问中的投影问题用线积分的形式来描述。例如上图中，记灰色的封闭图形由
23、函数 f(x,y)描述，一条投影射线为 xcos+y sin=t，则 f(x,y)关于该射线的线积分为：P(t)=(,t)linef(x,y)ds8获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657则这里的 P(t)就相当于这个物体沿该射线方向上的投影值。接下来引入傅里叶中心切片定理(Fourier Slice Theorem)。1 它证明了若给定一组投影的数据 p，则可以通过二维傅里叶变换来估计原物体。定义线积分投影 P(t)的傅里叶变换为：S(w)=P(t)ej2wtdt原
24、二维图像的傅里叶变换定义为：F(u,v)=f(x,y)ej2(ux+vy)dt则根据中心切片定理，有：S(w)=F(wcos,wsin)实际情况中，我们考察的 f(x,y)均为离散的网格化。假设其中 x,y 的所属区间均为A/2,A/2。对于 F(u,v)的傅里叶逆变换可以写为 4：f(x,y)=1A2N/2m=N/2N/2n=N/2F(m/A,n/A)ej2(m/A)x+(n/A)y)傅里叶重心切片定理将投影的傅里叶变换和图像的傅里叶变换建立了联系，故可以通过充足数目方向上投影的傅里叶变换信息，来构建出原始函数的估计值。在实际应用过程中，我们使用滤波反投影法（FBP）来完成重构，其算法步骤如
25、下：遍历 K 个角度，属于 0-180 之间1.对 P(t)进行快速傅里叶变换，得到 S(w)2.将 S(w)乘以 2|w|/K3.将 Step2 中得到的值进行傅里叶逆变换，将得到的值累加得到图像平面。5.3离散模型：ART5.3.1ART由于 FBP（滤波反投影）算法具有的缺陷与限制，对于不完全的投影数据，可以使用迭代法进行重建，这类算法在求解时把图像离散化。ART（代数迭代法）是常用的一种算法，其求解速度相对较快，结果较优，并且在求解过程中可以使用先验信息对数据进行修正。首先引入投影矩阵的概念，设图片共有 J=n n 个像素，共有 I 条射线进行扫描，将射线看作宽为 的粗线，将图片像素值
26、集中在像素的中心，则有 I J 大小的投影矩阵R，其元素定义为rij=1,i 号射线通过 j 号像素中心0,其他9获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657记 x=x1,x2,.,xJT为 J 维图像矢量，xj表示图像上第 j 个像素点的吸收强度；p=p1,p2.pIT为 I 维投影矢量，pi表示第 i 条射线所经过的所有像素的总射线投影，则Rx=p左侧各分量称为伪射线和，右侧称为真射线和。要将图像的吸收强度解出来，需要对矩阵R 进行求逆。一般的，R 是大型稀疏矩阵，因
27、此直接求其逆矩阵有困难；而在实际应用中，误差的存在可能使得 R 奇异。故通常使用迭代法进行求解，下面介绍 ART（代数迭代法）。其思路是每次校正一条射线路径上的像素值，使得该条真伪射线和间误差减小。其迭代公式为：x(k+1)=x(k)+(k)qik rikx(k)|rik|2rik式中，k 是迭代次数，k=0,1.，ik=k(modI)+1；(k)称为松弛参数。可以证明，当0 1(k)2=50&d 0)=1|.min(Res(R2_start:R2_start+73 1,1)0)=1break;endendforR1_last=512:1:35if(min(Res(R1_last 35+1:R
28、1_last,1)0)=1|.min(Res(R1_last 37+1:R1_last,1)0)=1).&min(Res(R1_last 38+1:R1_last,1)0)=0.&min(Res(R1_last 35+1:R1_last+1,1)0)=0break;endendR1_c1=R1_last 18;R2_c1=R2_start+36;phi0=rad2deg(asin(R1_c1 R2_c1)*0.2768/60);%Second FrameforR2_start=1:512 71ifmin(Res(R2_start:R2_start+71 1,180)0)=1.|min(Res(
29、R2_start:R2_start+73 1,180)0)=1break;endendforR1_last=512:1:35if(min(Res(R1_last 35+1:R1_last,180)0)=1.|min(Res(R1_last 37+1:R1_last,180)0)=1).&min(Res(R1_last 38+1:R1_last,180)0).=0&min(Res(R1_last 35+1:R1_last+1,180)0)=028获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5
30、44457657break;endendR1_c2=R1_last 18;R2_c2=R2_start+36;R_p=(R1_c1+R2_c1+R1_c2+R2_c2)/4;d1=(R1_c1+R2_c2)/2;%Third FrameforR2_start=1:512 71ifmin(Res(R2_start:R2_start+71 1,90)0)=1.|min(Res(R2_start:R2_start+73 1,90)0)=1break;endendforR1_last=512:1:35if(min(Res(R1_last 35+1:R1_last,90)0)=1.|min(Res(R1
31、_last 37+1:R1_last,90)0)=1).&min(Res(R1_last 38+1:R1_last,90)0).=0&min(Res(R1_last 35+1:R1_last+1,90)0)=0break;endendR1_c3=R1_last 18;R2_c3=R2_start+36;O1=(R1_c1+R2_c1)/2;%O2=(R1_c2+R2_c2)/2;O3=(R1_c3+R2_c3)/2;R_p=(R1_c1+R2_c1+R1_c2+R2_c2)/4;R_r=0.2768*sqrt(O1 R_p)2+(O3 R_p)2);R_theta0=180+rad2deg(a
32、tan(O1 R_p)/(O3 R_p);29获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657R_theta=phi0 R_theta0;phi0=phi0;x0=R_r*cos(deg2rad(R_theta);y0=R_r*sin(deg2rad(R_theta);d0=256 R_p;x_init=x0,y0,d0,phi0:10:phi0+10*17;A=;b=;Aeq=;beq=;lb=;ub=;nonlcon=;options=optimoptions(fminc
33、on ,Display ,i t e r ,.Algorithm ,sqp ,MaxFunEvals ,15000,TolFun,1.000000e 10);%phiR1=5;R2=10;R3=5;d0=1;mu=1.7725;x0=0;y0=6;phi=deg2rad(30:149)+180*ones(1,180)+rand(1,180);Res=zeros(512,180);foriPhi=1:180for n=1:512theta=phi(iPhi);d=d0 0.2768*(256.5 n);Res(n,iPhi)=2*mu*(sqrt(max(0,R12.(tan(theta)*(3
34、0)x0*tan(theta)+y0 d/cos(theta)2./(1+tan(theta)2)+.sqrt(max(0,R22 (tan(theta)*(30)x0*tan(theta).+y0 d/cos(theta)2/(1+tan(theta)2);endend30获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657index=find(Res 0);Res(index)=Res(index)+0.01;%normrnd(0,0.1,length(index),1);save(Res.mat ,Res )31获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657
