CT系統(tǒng)參數標定數學模型

CT系統(tǒng)參數標定數學模型摘要本文通過對體模系統(tǒng)和投影點坐標的分析，提出了較為充分合理的假設，并對數據進行了擬合處理。先后建立了兩個模型（理想狀態(tài)下的模型，只考慮載物臺和探測器等一些因素誤差的模型）來求解CT系統(tǒng)參數標定的問題。在對問題一的分析中，只考慮在理想情況下：即載物臺水平和探測器無偏轉，并不考慮機械誤差的情況下各個參數之間的關系，我們以探測器左下角為（0，0），右方向為x軸正方向，上方向為y軸正方向，垂直于xoy平面的走直線為z軸建立如圖4所示的空間直角坐標系。運用幾何知識作圖，通過建立相似圖形的等比例方程來研究確定幾個系統(tǒng)參數之間的關系（X射線源的位置參數、載物中心線的位置參數和探測器的位置參數）。最后討論了引入載物臺中心軸傾斜、載物轉臺不均勻轉動和數據不確定度對誤差的影響。在分析問題二時，考慮到所給橢圓點集可能有偏差，于是首先運用MATLAB軟件對數據進行處理擬合，得到更為確信的數據并由此求出修正后的“橢圓”方程。從擬合出來的橢圓方程，我們可以得知橢圓的中心坐標還有其大小，并因此看出我們擬合出來的橢圓并不和我們的理想模型相符合。所以我們想到了有機械誤差存在。由于同時考慮多個誤差因素將使模型計算變得極為復雜，對此，我們假設探測器無偏轉，只有載物臺傾斜，建立模型二，并對探測器上的兩個橢圓進行旋轉修正。易知修正前后兩個橢圓的大小并沒有改變，只是相對位置發(fā)生了變化，發(fā)射源A在探測器上的投影坐標也沒有改變。對此，我們以A點為旋轉中心對其進行旋轉修正并根據擬合后的數據求出兩個類橢圓中心連線的斜率。從而求出載物臺傾斜角ａ=arctan0.00210563最后在問題三，對模型三的機械誤差進行了分析計算。機械誤差主要來源為探測器上數據像素誤差、載物臺的轉軸偏斜，圓心距離的誤差，并對其分析了原因。關鍵詞CT系統(tǒng)載物臺轉軸偏斜數據修正擬合機械誤差參數標定問題重述CT自發(fā)明以來，被公認為自倫琴發(fā)現X射線后是放射領域的最重要的發(fā)明之一。工業(yè)CT一般采用掃描物件旋轉和面探測器的結構，工作示意圖如圖1示：圖SEQ圖\*ARABIC1工業(yè)CT系統(tǒng)工作原理示意圖掃描物件圍繞某一固定轉軸旋轉，每隔一定角度采集一張圖像，然后根據采集的圖像采用3D圖像重建算法即可將原始3D物件重建出來。在實際中，由X射線源發(fā)出X射線，經過掃描物體衰減后照射在探測器上，探測器根據接收到的光子數的計數實現光電轉換，從而形成灰度圖像。一般平板探測器大小為3000×2000像素，每個像素為0.127mm。在CT系統(tǒng)安裝過程中往往存在機械誤差，而這些誤差對于物件重建的準確性往往是致關重要的，實際中就需要對安裝好的CT系統(tǒng)進行參數標定，系統(tǒng)參數主要包括X射線源的位置、載物中心線的位置和探測器的位置參數等。傳統(tǒng)的方法是采用尺子測量法，對于CT系統(tǒng)參數標定精度要求是不夠的（通常CT重建過程中系統(tǒng)參數的機械誤差不允許超過一個像素）。目前，一般通過實驗的方法實現參數標定。所謂實驗的方法標定參數就是通過掃描已知參數的體模，分析投影數據，估計系統(tǒng)參數的方法。針對工業(yè)CT要求解決如下問題：（1）建立合適的坐標系，正確描述CT系統(tǒng)的各種參數和機械誤差，并建模分析這些參數的關系和可能的機械誤差。（2）通常采用軸承鋼球作為CT參數標定體模。因為它具有各向投影一致，邊緣清晰等優(yōu)點，在里采用兩個鋼球實現CT的參數標定實驗（如圖3所示）。兩個鋼球置于有機玻璃管中（X射線容易透過有機玻璃，容易后期球心獲取的圖像處理），鋼球直徑為8mm±0.0008mm，兩球心距離約為100mm±1mm。將有機玻璃管固定在旋轉載物臺上，載物臺攜帶有機玻璃管以均勻速度旋轉，與此同時X射線源發(fā)出射線，探測器采集數據（旋轉一周，等間距采集180張）。附件1給出了根據180張采集圖像提取出的球心投影坐標。兩個鋼球球心投影的軌跡是兩個橢圓（如圖4所示）。試根據實驗數據估計該CT系統(tǒng)的參數值，即給出CT系統(tǒng)的標定。（3）在參數標定過程中可能存在多種可能的機械誤差，試就你的參數標定可能的誤差進行分析。圖SEQ圖\*ARABIC2雙球體模型示意圖圖SEQ圖\*ARABIC3雙球投影軌跡（雙球軌跡坐標見附錄1）問題分析題目要求我們求出工業(yè)CT系統(tǒng)的各項參數，這些參數包括以下六個：=1\*GB3①放射源到轉盤的轉軸的水平距離L1；=2\*GB3②轉盤的轉軸到探測器平面的水平距離L2；=3\*GB3③放射源到探測器平面的水平距離L3；=4\*GB3④放射源到xoz平面的垂直距離Ya；（平面xoz具體位置見問題（1）解答中的坐標系）=5\*GB3⑤放射源與yoz平面的水平距離Xa；=6\*GB3⑥轉軸到y(tǒng)oz平面的水平距離d。如圖所示：為了使參數更容易確定，我們選取探測器作為參照物，這樣各個參數就有了參考的標準。問題一要求我們求出系統(tǒng)各個參數之間的關系，還有根據我們所列的關系式進行誤差分析。本題誤差主要來源于工業(yè)CT系統(tǒng)安裝過程。除了裝置本身的由于在裝系統(tǒng)時，可能會導致物臺中心軸傾斜、載物轉臺不均勻轉動，探測器前后俯仰等，都會導致機械誤差的存在。我們在計算時會考慮這些因素,并建立相應的模型。所以，在標定時應該給出其機械誤差大小?？梢愿鶕w模的投影數據，利用幾何知識來求解主要包括X射線源的位置、載物中心線的位置和探測器的位置并求出他們的關系。問題二則以題目給出的360個坐標為基礎，運用matlab軟件來擬合兩個橢圓，考慮到這些點的坐標本身就有一點的誤差，經過一定的修改和篩選后，進行多次擬合進而求出橢圓的方程。求出橢圓的方程后，再利用在問題（1）中各個參數的關系，就可以求出各個參數的值。根據分析，我們知道了擬合出來的兩個橢圓的中心連線并不垂直于X軸，我們可以利用繞點旋轉的方法對其進行修正。問題三中的誤差分析，本題中考慮的機械誤差來源有兩圓心距離誤差，載物臺的中心轉軸傾斜，這些根據問題二中的數據，我們已經對它進行了修正。模型的假設與符號說明模型假設假設光子波動不明顯，即忽略光的衍射現象；將兩個鋼球的球心看作一個質點；假設實驗所收集的數據能客觀反映實際情況；玻璃管在隨著轉盤轉動的過程中，始終垂直固定于轉盤；載物臺在轉動過程中，玻璃管與兩個鋼球始終保持相對靜止。符號說明及名詞解釋符號符號含義數據類型L放射源到轉軸的水平距離待求參數L轉軸到探測器平面的水平距離待求參數L放射源到探測器平面的水平距離待求參數X放射源到y(tǒng)oz平面的距離待求參數Y放射源到xoz平面的垂直距離待求參數d轉軸到y(tǒng)oz平面的距離待求參數r有機玻璃管和小球的軌跡圓半徑中間值a上橢圓的短軸中間值a下橢圓的短軸中間值b兩橢圓長軸的平均值中間值m兩橢圓最高點間的距離中間值k兩橢圓最低點間的距離中間值t放射源在探測器上的投影與上橢圓最低點的距離中間值h兩鋼球球心間的距離已知量模型建立及問題求解問題一建立合適的坐標系，正確描述CT系統(tǒng)的各種參數和機械誤差，并建模分析這些參數的關系和可能的機械誤差。將探測器所在的平面定為參照物，結合題目已給出的二維坐標系xoy，在探測器平面左下方原點（0，0）引一條Z軸，形成三維坐標系，如下圖示：依照題意，我們首先得求出載物臺的轉盤半徑r，為了求出轉盤的半徑，我們先作了這樣的一個實驗，我們以探測器左下角為（0，0），右方向為x軸正方向，上方向為y軸正方向，垂直于xoy平面的走直線為z軸建立如圖1所示的空間直角坐標系。，再讓發(fā)射源與玻璃管中的下邊小球處于同一個小平面上進行實驗。在實驗過程中，我們能在探測器上得到玻璃管中上面小球所形成的一條直線為2b，而下面小球的投影則形成一個類橢圓，我們可以通過探測器上的數據擬合得到橢圓的方程，從而得到我們想要橢圓的長軸長和短軸長還有中心坐標，進而可以計算出有機玻璃管和小球的軌跡圓半徑r。圖4為了簡單分析，我們分別畫出圖1的正視圖和俯視圖如下：俯視圖正視圖根據正視圖和俯視圖，我們利用相似三角形和勾股定理得到了以下關系式：L1+L1-rrb=聯(lián)立（1）（2）（3）可解得：r=b*h2*m*h1根據以上的式子，我們可以根據實驗得出的數據求出r的值，然而對于相同的這個裝置來說，以下的實驗中有機玻璃管和小球的軌跡圓半徑r是不會改變的，還有因為其俯視圖中的2b的值是發(fā)射源垂直于xoy平面發(fā)射出來的投影的長度，而在其他情況下是求不出那個值，為了計算方便，在以下的運算中，我們將有機玻璃管和小球的軌跡圓半徑r作為一個常數a進行計算，即令r=a（a則為由上面（4）求出的數值）接著，我們按照題目的要求保持發(fā)射源A與載物臺和探測器的水平距離不變，即相對位置不變，而在豎直方向上移動發(fā)射源A，使其與兩個小球都不在同一個水平面進行以下實驗，其三維圖如下：為了方便分析，我們作出它的正視圖如下：根據其正視圖，利用相似三角形定理和勾股定理可以得到以下關系式：L1+rL1-ryB-r=a（其中a為一個常數）（8）根據CT系統(tǒng)的工作原理，即令掃描物件圍繞某一固定轉軸旋轉，并每隔一定角度采集一張圖像，然后根據采集的圖像采用3D圖像重建算法將原始3D物件重建出來。再根據重建出來的圖像來標定我們的系統(tǒng)參數，由此可知，凡是在重建圖像所需要的數據誤差都會影響我們的參數標定。因此，根據上方我們所列出來的關系式可以看出此過程中存在的誤差有以下這幾種：載物臺沒有正常勻速轉動，引起重建圖像不連續(xù)掃描物件的放置發(fā)生傾斜或在運動過程產生抖動，引起重建圖像不規(guī)則載物臺轉動角度發(fā)生偏離或載物臺的轉軸發(fā)生傾斜，引起重建圖像的數據不好處理實驗中兩個鋼球的半徑大小有一定的誤差，使得有機玻璃管和小球的軌跡圓半徑r存在誤差從而影響參數標定。問題二根據照題目給出的數據，我們用MATLAB軟件進行擬合，繪制出上下兩個橢圓的圖形（算法詳見附錄2），由圖形反映擬合情況，并解出橢圓一般方程的系數p1、p2數據圖形擬合上方橢圓的中心坐標為（1421.7199，1812.5604，0）長半軸長為1351.2899，短半軸長為172.8553（單位為像素）下方橢圓的中心坐標為（1427.2216，188.4994，0）長半軸長為1342.6651，短半軸長為上方橢圓的方程為：(x-1421.7199)下方橢圓的方程為：(x-1427.2216)由此，我們可以得出B,D,C,E各點的縱坐標分別為y將其坐標代入（5）（6）（7）（8）式可得：LLLya=811.83225，x即我們可以得出發(fā)射源A的坐標為（1423.832，811.83225，11.58469ａ）綜上所述，我們可以對給定CT系統(tǒng)參數進行標定（單位為像素）：X射線源的位置與載物臺中心線的距離：LX射線源的位置與探測器的距離：L若以坐標和方程的形式給出，則表示為：發(fā)射源A的坐標為（1423.832，811.83225，11.58469ａ）載物臺中心線的方程x=探測器的位置：x-o-y平面模型的進一步分析與修正在上述的分析過程中，我們發(fā)現經過擬合出來的兩個橢圓的中心連線并不垂直于X軸，與實際中我們預測分析的有所不同。由于影響系統(tǒng)的機械誤差有很多，在此，我們只認為這是因為載物臺在轉動過程中所繞的轉軸傾斜產生偏角所導致的，與其它因素無關。所以，我們對其作了修正如下（左為修正前，右為修正后）：由問題二我們已經求出了發(fā)射源A的坐標為（1423.832，811.83225，11.58469），由于此過程中只在xoy平面上修正，所以修正后發(fā)射源的坐標是不會改變的。對此，我們以A點在xoy平面上的投影A’為中心，線段OBOC繞著中心A’旋轉到線段OBOC垂直于X軸，此時，上下兩個橢圓的中心坐標分別為OB（1423.832，1812.56265，0上方橢圓的方程：(x-下方橢圓的方程：(x-此時我們還可以求出修正的角度ａ，其修正的圖形在xoy平面的角度ａ如下圖所示：以上方橢圓的中心坐標為標準，則修正前的角度ａ有如下關系：tan代入數據可以求出：tanａ即修正的角度的大小為：ａ=arctan經過了這樣的修正，我們可以更加準確的利用探測器上的數據進行擬合圖像，使得重建的圖像更加具有準確性，我們能夠更好的進行標定CT系統(tǒng)的參數。問題三除了問題二中所提到的載物臺的轉軸偏斜引起誤差之外，還存在以下兩種誤差可能引起參數標定不準確。首先，如果載物臺沒有正常勻速轉動，那么掃描物件的旋轉速度會受到影響，探測器所采集到的具體數據可能會發(fā)生斷點或重點，如下圖的斷點情況示意圖：其次，掃描物件的放置發(fā)生傾斜或在運動過程產生抖動會導致探測器所采集到的圖像數據發(fā)生偏差，影響了Matlab的擬合效果以及相關數據的準確性。如圖12所示，很明顯地，投影得到的圖形發(fā)生了扭曲，形成機械誤差，影響到了參數的標定。模型評價與改進優(yōu)點：對體模的投影點坐標進行數據的修正，使得計算的結果更加準確；考慮到的誤差比較合理，有較大的實際意義運用MATLAB軟件進行數據運算和圖表擬合，使得計算更加簡單準確。缺點：機械誤差存在的因素考慮還不夠全面，且僅限于理論上的分析，有待實踐過程中的進一步檢驗與改進在數據誤差處理時僅用理論知識進行處理，實際操作有些難度。參考文獻[1]百度文庫/[2]中國知網/kns50/index.aspx[3]謝照鴻，范正森等人編著，《數學建模技術》，中國水利水電出版社，2003年9月[4]周永正，詹棠森，方成鴻，邱望仁，《數學建?！?，同濟大學出版社，2011年1月附錄附錄1雙球模型投影在探測器上的坐標上方球心坐標下方球心坐標x軸y軸x軸y軸2762179227372142754178627272172745178127152202735177527012242723176926872272710176426712312695175826532342679175326352372661174726152402643174225942432623173725722462602173225492492580172725252522557172225002552532171724742572507171324472602481170824192622454170423902652425170023602672396169623302692367169222992712336168822672742305168422352752273168122022772240167721682792207167421342812173167120992832138166820642842103166520282852068166319922872032166019552881995165819182891959165518812901921165318432911884165118052921846165017672931807164817282941769164616902941730164516512951691164416122951652164315732961613164215332961573164114942961534164014552961494164014152961454164013762961415163913372961375163912972961336164012582951296164012192951256164011802941217164111412931178164111032931139164210642921100164310262911061164498829010231646951289985164791428894716498772869101651841285873165280528383616547692828001657734280764165970027872916616662776941664632275660166760027362616705682715931673536268561167650626652916794762644991683447261468168641825943916903912564101694364254383169833925135617023142483301706291245305171126824228117152472392591720227236237172520723321617301902301971735173226179174015822316217451442191461750131216132175612021211917611102091071767102205971773952028917789019882178486194761790841907217968318770180285183691808881797018149217573182199171771827107167841833117164921839129160102184514315611318521581521271858176148143186419514516018702171411791877240137201188326513422418892911302491895320127276190035112330419063831203351912417117367191745311440119234901104381928529108475193357010551419386121025551942656100598194770197642195174795687195579593734195984491782196389389832196694487882197099686933197210498598619751103841039197711578310931979121182114819811267811203198313228112591984137881131519851434811371198514898114271986154581148419861600821540198516558315961985171083165219841764841707198318178617611981187087181519801922891868197819739119211975202392197219732071952022197021199720721967216599212019632210102216619602253104221219562295107225519522336110229719482374113233919432412116237819382447119241519342480123245019292512126248419232542130251619182570133254719132596137257519072621140260219012643144262618962664148264918902683152267018842699156268918782714159270618722727163272218662739167273518592748171274618532756175275618472761179276418412765183277018352768186277418282768190277618222767194277718162765198277618102760201277318042754205276917982747209注：探測器左下角為（0，0），右方向為x軸正方向，上方向為y軸正方向附錄2用matlab軟件求橢圓方程的程序和輸出結果：求上橢圓方程（先將上橢圓的180個點的坐標導進，記作x1，y1，用矩陣表示）F=@(p,x)p(1)*x(:,1).^2+p(2)*x(:,2).^2+p(3)*x(:,1)+p(4)*x(:,2)+p(5);p0=[0.0050.0050.0050.0050.005];x=[x1,y1]warningoffp=nlinfit(x,zeros(size(x,1),1),F,p0)plot(x(:,1),x(:,2),'ro');holdon;xmin=min(x(:,1));xmax=max(x(:,1));ymin=min(x(:,2));ymax=max(x(:,2));ezplot(@(x,y)F(p,[x,y]),[-1+xmin,1+xmax,-1+ymin,1+ymax]);formatlongpA=sqrt(p(3)^2/(4*p(1)^2)+p(4)^2/(4*p(1)*p(2))-p(5)/p(1))B=sqrt(p(3)^2/(4*p(1)*p(2))+p(4)^2/(4*p(2)^2)-p(5)/p(2))x0=p(3)/(2*p(1))y0=p(4)/(2*p(2))C=A^2D=B^2輸出結果：p=1.0e-022*0.00000.0000-0.0000-0.00010.1249p=1.0e-022*0.000000000621500.00000003798134-0.00000176718735-0.000137686937760.12490432520128A=1.351289908892398e+003B=1.728