CT系统参数标定

1、CT系统参数标定问题研究摘要本文通过对体模系统和投影点坐标的分析，提出了较为充分合理的假设，并对数据 进行了拟合处理。先后建立了三个模型（理想状态下的模型，只考虑载物台倾斜的模型， 考虑载物台倾斜与探测器的俯仰角的模型）來求解CT系统参数标定的问题。在对问题一的分析中，只考虑在理想情况下：即载物台水平和探测器无偏转，并不 考虑机械误差的情况下各个参数之间的关系，以X光放射源为原点建立坐标系，沿水平 力向射向探测器的水V直线为X轴，竖点方向为Y轴，建立模型一，绘出侧视图如一所 示。运用几何知识作图，通过建立相似图形的等比例方程來研究确定儿个系统参数之间 的关系（X射线源的位置参数、载物中心线的位

2、置参数和探测器的位置参数）。同时求出 载物台半径等参数，最后讨论了引入像素、载物台偏角、观测器俯仰角和数据不确定度 对误差的影响。在分析问题二时，考虑到所给椭圆点集可能有偏差，于是首先运用拟合的方法对数 据行处理，得到更为确信的数据并由此求出修正后的“椭圆”方程（观察均方根误差, 得知此方程能符合条件，可以使用）。由于考虑多个误差因素将使模型计算变得极为复 杂，在此不同时考虑。假设探测器无偏转，只有载物台倾斜，建立模型二。易知倾斜 的载物台与安装没有偏差的水平观测器可以等效为水平的载物台与只纵向倾斜的观测 器。然后根据拟合后的数据求出两类椭圆长轴斜率，从而求出载物台倾斜角为0. 1146 &#

3、176;, 并依据此进一步修正数据。在模型二修正数据的基础上，考虑探测器的俯角5因素（易知S不影响模型二的求 解）建立模型三。由修正后的图知两“椭圆”长轴长等。由上方长轴长大于下方长轴长 知探测器向前倒。经过一定数据近似后，再根据等比例关系可求得探测器俯角5为 1.623。近一步推导可得坐标及光源到旋转台的旋转轴的距离为263.97mm,光源在成 像平面上的水平直射点，光源到探测器水半直射点的间距为688.57mm,。最后在问题三，对模型三的机械误差进行了分析计算。机械误差主要來源为像素、 载物台偏角，探测器的俯仰角，圆心距离的误差。并进行了灵敏度分析。关键词：坐标修正载物台偏角探测器俯仰角机

4、械误差一、问题的重述扫描物件围绕某一固定转轴旋转，每隔一定角度采集一张图像，然后根据采集的图 像采用3D图像重建算法即可将原始3D物件觅建出來。在实际中，由X射线源发出X射 线，经过扫描物体衰减后照射在探测器上，探测器根据接收到的光子数的计数实现比电 转换，从而形成灰度图像。一般平板探测器大小为3000X2000像素，每个像素为0. 127mm。在CT系统安装过程中往往存在机械误差，而这些误差对于物件重建的准确性往往 是致关重要的，实际中就需要对安装好的CT系统进行参数标定，系统参数主要包括X 射线源的位置、载物中心线的位置和探测器的位置参数等。对J：CT系统参数标定精度 要求是不够的（通常C

5、T重建过程中系统参数的机械误差不允许超过一个像素）。通过扫 描己知参数的体模，分析投影数据，估计系统参数的方法。针对工业CT要求解决如下 问题：（1）建立合适的坐标系，正确描述CT系统的各种参数和机械误差，并建模分析这 些参数的关系和可能的机械误差。（2）通常采用轴承钢球作为CT参数标定体模。因为它具有各向投影一致，边缘清 晰等优点，在里采用两个钢球实现CT的参数标定实验.两个钢球置于有机玻璃管中（X 射线容易透过有机玻璃，容易后期球心获取的图像处理），钢球直径为8mm±0. 0008mm, 两球心距离约为100mm±lmmo将有机玻璃管固定在旋转载物台上，载物台携带有机玻

6、璃 管以均匀速度旋转，与此同时X射线源发出的射线，探测器采集数据（旋转一周，等间 距采集180张）。两个钢球球心投影的轨迹是两个椭圆。试根据实验数据估计该CT系统 的参数值，给出CT系统的标定。（3）在参数标定过程中可能存在多种可能的机械误差，对参数标定可能的误差进 行分析。二、问题的分析问题一：本题误差主要來源于工业CT系统安装过程。除了装置本身的由于在装系统 时，可能会导致探测器扭转、前后俯仰，载物台的倾斜等，都会导致机械误差的存在。 我们在计算时会考虑这些因素,并建立相应的模型。所以，在标定时应该给出其机械误 差大小。可以根据体模的投影数据，利用儿何知识来求解主要包括X射线源的位置、载

7、物中心线的位置和探测器的位置并求出他们的关系。对丁载物台來说，由F载物台要进 行旋转，所以其倾斜方式只有一种。问题二：根据分析，载物台的倾斜可以等效为探测器投影绕一个点在竖直方向旋转 （图三所示）。由附录一中图分析可知，两个椭圆的长轴不是水平的，这是由于载物台 的倾斜造成的：上椭圆与下椭圆的长轴不一样长，这是由于探测器的前后俯仰造成的:“椭圆”图形左右不对称，是由于探测器的扭转造成的。我们分别建立模型将数据一步 一步进行处理，得到机械误差来源，再进行修正，从而给出准确的CT参数标定。问题三：本题中考虑的机械误差來源有两圆心距离误差，载物台倾斜，探测器前俯 后仰，利用问题二中所得数据，可求得总机

8、械误差。三、模型的假设与符号说明1. 模型的假设（1）考虑到此处光子的波动不明显，故忽略光的衍射现象，以光的直线传播进行处理；（2）忽略X光辐射源的体积和形状，将其视为质点：（3）一些题目中未说明但实际易得的参数，如载物台的半径等，假设可以测量得到， 此处不考虑具体方法及相关课差；（4）数据能客观反映实际情况：（5）载物台与有机玻璃管始终垂直；（6）两小球圆心视为质点；（7）在数据采集过程中，放射源点，载物台，探测器参数不再变化。2. 符号说明11： X光射线源到载物台中心垂线的距离；12： X光射线源到探测器的距离；13:探测器俯仰偏转后下方椭圆最下沿点到Y轴距离：14：探测器俯仰偏转后上方

9、椭圆最上沿点到Y轴距离； r：有机玻璃管和小球的轨迹圆半径； h：两个小球重心之间的距离； h1：椭圆间最近两点的距离：h,：椭圆间最远两点的距离：n1：上投影椭圆短轴长度；n2：上投影椭圆短轴长度； m：投影的竖直方向最大长度；载物台的倾斜角；%探测器偏转的俯仰角；tl：大椭圆的长轴长度；t2：小椭圆的长轴长度；tl2：两椭圆中心之间的距离；% :当偏差为一个像素时候，探测器的俯仰角度：At1：当偏差为一个像素时候的变化量；At3：当偏差为一个像素时候t?的变化量；At12：当偏差为一个像素时候站的变化量四、模型的建立和求解问题一：1、模型一的建立苗先考虑在理想情况下，即载物台水平和探测器无

10、偏转和不考虑机械误差的情况下 齐个参数之间的关系。建立坐标系时选定X光辐射源为坐标原点，沿水半方向射向探测 器的水平直线为X轴，竖直方向为Y轴，绘出侧视图如一所示。图二理想模型下的俯视图由平面儿何中相似图形成比例可以得出下列方程：h.li-i111 12hL + rm- (iij + n?)12r , km nm72肩9若考虑引入机械误差，则可能存在以下儿种机械误差：（1）载物台与地面有一定倾斜角，此角会使投影偏转。（2）探测器地而间存在一定俯仰角，此偏角的存在会使原本应是椭圆的投影产生一定 偏差：（3）x放射点与载物台轴线成面与探测器成角不成90度，会使原本应是椭圆的投影产 生一定偏差;（4

11、）像素不够大引起的误差，一般平板探测器大小为3000x2000,每个像素为0. 127inn若两个投影点之间的水平或竖直距离小于0. 127inn.在探测器上就会区分不出，例如以下 2组数据：上方球心坐标 （X, Y）： （1415, 1639）, （1375, 1639）；（5）题目中给出的数据（小球半径.球心距离等）不精确产生的误差；问歸一.1、處二的建立暂只考虑载物台的倾斜角Q的情况下，建立此模型。倾斜的裁物台与安装没有偏差的观测器可以等效为水半的载物台与只纵向倾斜的观测 器，如图二由丁数据可能有较大的误差，我们用mat lab对数据进行拟合。530002000图四：修正前的图形修正后1

12、800160014001200100080060040020050010003000150020002500图五：拟合修正后的图形°0#再由图 中所示方法求出长轴短轴大小及对应各个切点。於平行切线Z何的距离故 短为炷轴图六：椭圆的切线与椭I貝轴的关系再根据对应切点确定长轴斜率 tanx=(lxj-rxj)/ (lx-rx);电下方椭圆长轴斜率tans= (Is j-rs j ) / (ls-rs) ;E上方椭圆长轴斜率tan=(tans+tanx)/2;得 q)2 = arctaiiO.OO2根据相应切点求出圆心轨迹所成圆之圆心，水平射于探测器的点横纵坐标 xsO=(ls+rs)/2

13、;ysO=(lSj4-rsj)/2;%上方圆心轨迹圆圆心水平射于探测器的点横纵坐标xxO=(lx+rx)/2;yxO= (lxj + rxj )/2;%卜方圆心轨迹圆圆心水平射于探测器的点横纵坐标bili=(gs-ds)/(gs-ds)+(gx-dx);根据附件一中所给数据，计算出X放射点水平射于探测器的横纵坐标 yO=ds- (gx-ds) *bili;%x放射点水平射于探测器的纵坐标 xO=dsi- (gxi-dsi) *bili;Ex 放射 点水平 射于探 测7修正后坐标轴勺Y2修-5正前处标轴t Xo.Vo>#图七修正的原理图再根据对应切点确定长轴斜率h=abs(sqrt(tan

14、A2+l)一(tan*(x(i,j)-Is)-y(i,j)+lsj)/sqrt(tmnA2+l);$上方椭圆上任意一点至上方椭圆长轴的距离y(irj)=yO+(ysO-yO)*sqrt(tanA2+l)-(tan*(x(izj)-Is)-y(i,j)+lsj)/ sqrt (tanA2+l);咎禧正后的纵坐标x (i, j ) =xO+fuhao*hh;电修正后的纵横坐标修正后图形如图五。详细数据见附录一。2、模型三的建立我们可以建立在模型二的基础上考虑探测器的俯仰角度幻，如图八所示，由于载物台 的偏转角度Q与探测器俯仰角度幻不影响，所以可以把他们结合起来计算。考虑他们偏 转角的误差。由图八和

15、图可知，探测器向前偏转了匀，导致上椭圆和下椭圆大小不一样，根据几 何知识，可以求得方程(4)和方程(5)：于是，我们就可以求出在探测器前后俯仰时候的 偏转角5与体模投影点坐标之间的关系。找出偏转角s与其他参数之间的关系，进行修 正。得到准确的载物中心线的位置和探测器的位置参数，较楮确对CT进行标定。而tl与t2的值非常相近，即纟非常小。相似三角形对应边要小于对应1 cos幻，但由 与的非常小，相差不大，二者可约等于。'则由图三和图四可以得到下列方程如下：9llj cos%7-)h2 cos%9tlS非常小。则可以将v i3, 13的长度看成一样。即式子可改为:1、cos （fl h2

16、cos .（片tjlr&于是，我们就可以求出在探测器前后俯仰角3 =1.62% 由模型三的修正，我们可以得到准确的数据，就可以求得： h = 787 m = 1905 珂=347 n2 = 215 t】=2708 t3 = 2685 屯=1625.5 带入模型一的方程,由模型一，X射线源的位置在破点，即坐标为（00）,由射尬源到载物中心线的距离为1 = 2078.55（像素）=263.97inn】，即位置坐标为（0,263.97）,由射线源到探测器的距离为12 = 5421.85（像素）=68&57mm ,即位置坐标为（0,688.57）,我们还求出在载物台上物体绕载物台中心线

17、旋转的半径, 由于小球在载物台边缘，则旋转半径约等于载物台半径rr=468.56（像素）=59.5mm问题三：根据前面的描述，误差的来源主要有三种：像素、偏转角和球心距。对于CT系统 参数标定精度要求是其机械误差不超过一个像素，而由式（1）（2） （5）知球心距只与h有直接联系，与像素无直接联系,故此处主耍研究偏转角与投影像素Z间的关系。有式11q = At2怙血0=心土凹b tl+“I（5）可进一步推出它们的关系如下:(8)(9)(10)计算（8）得到>=2.009°：笑一0.386。； 计算（10）得却%=1764°；爲一的=0.141°；比较两次结果，

18、取最小值0.141%即要保证投影在探测板上的点位误差不超过一个 相位，探测板的俯仰角浮动范围最大不能超过0.141。五、模型的推广与改进方向考虑三维直角坐标系，能更精确的求出机械误差，具有更好的实际应用。考虑探测 器的扭转角对体模投影图像的影响。考虑多种误差同时影响O六、模型的优缺点优点：1、建立多个模型，模型的建立由易到难，比较有层次，能更清晰的分析问题，找 到它们之间的关系：2、对体模的投影点坐标进行数据的修正，计算的结果更加准确：3、考虑到的误差比较合理，有较大的实际意义；缺点：1、对探测器的扭转角考虑不全面，仅在理论上进行分析。2、由于像素造成的机械误差比较大，像素精度不高，此方法解决

19、这个问题效果欠 佳；七、参考文献1姜启源.数学模型Ml.北京：高等教育出版社，1999. 韩中庚.数学建模方法及其应用M.北京：高等教育出版社，2009.八、附录附录一：投影点坐标的修正 zhen=zeros (2000, 3000); for i=l:180for j=l:2zhen(y(i, j), x(i, .1)=1: endend%矩阵Z生成rs=1000;lx二1000;rx=1000;for i=l:3000for j=1:2000 ifzhenf j, i) =1 if i<lsif j>1400ls=i;lsj=j;endendif i<lxif j<

20、500lx=i;Ixj巧； end end if i>rs if j>1400 rs=i ; rsj=j;endendif i>rx if j<500 rx=i;rxj=j；endendendendendgs=1800; ds=1800;for i=l:3000 for j=l:2000 ifzhenfj, i) =1 if j>1000 if j<ds ds=j;dsi=i;endendend ifzhen； j, i)=l if j>gs ss=j;gsi=i;endendendend gx=200;dx=200； for i=l:300013fo

21、r j=1:2000 ifzhen! j, i) =1 if j<dxdx=j;dxi二i;endendifzhen ： j, i) =1if j<1000if j>gxgx二 j;gxi=i;endendendendendtanx二(lxj-rxj)/(lx-rx); tans=(lsj-rsj) /(ls-rs); tan=(tans+tanx) /2;%下方椭圆长轴斜率%上方椭圆长轴斜率xs0=(ls+rs) /2;ys0=(lsj+rsj)/2;标xxO二(1x+rx)/2;yxO=(lxj+rxj)/2;坐标bili = (gs-ds) /(gs-ds) + (gx-dx);yO=ds-(gx-dS)*b订i;%x放射点水半射于探测器的纵坐标 xO=dsi-(gxi-dsi)*bi 1 i :%x放射点水V射J：探测器的横坐标%上方圆心轨迹圆圆心水半射于探测器的点横纵坐%下方圆心轨迹関圆心水平射探测器的点横纵for j=l：2for i=l:180辻(y(i, j)>1500)h=(tan*(x(i, j) -xsO)-y(i, j)+ysO) /sqrt (tan 2+1);hh=sqrt (x(i, j)-xs0) "2+(y(i, j)-ys