王 琪

（江西財(cái)經(jīng)大學(xué)經(jīng)濟(jì)學(xué)院， 江西 南昌 330000）

1 問(wèn)題重述

CT（Computed Tomography），即計(jì)算機(jī)斷層掃描，利用樣品對(duì)射線能量的吸收特性對(duì)樣品進(jìn)行斷層成像，由此獲取樣品內(nèi)部的結(jié)構(gòu)信息。平行入射的X射線垂直于探測(cè)器平面，每個(gè)探測(cè)器單元看成一個(gè)接收點(diǎn)，且等距排列。X射線的發(fā)射器和探測(cè)器相對(duì)固定不變，整個(gè)發(fā)射—接收系統(tǒng)繞某旋轉(zhuǎn)中心逆時(shí)針旋轉(zhuǎn)180次。對(duì)于每一個(gè)X射線方向，在探測(cè)器上測(cè)量經(jīng)位置固定不動(dòng)的二維待檢測(cè)介質(zhì)吸收衰減后的射線能量，并經(jīng)過(guò)增益等處理后得到180組接收信息。為減小誤差，需要對(duì)安裝好的CT系統(tǒng)進(jìn)行參數(shù)標(biāo)定，并據(jù)此對(duì)未知結(jié)構(gòu)的樣品進(jìn)行成像[1-2]。

問(wèn)題1：在正方形托盤上放置兩個(gè)均勻固體介質(zhì)組成的標(biāo)定模板，模板的幾何信息如圖1所示，其中每點(diǎn)的數(shù)值反映該點(diǎn)的吸收強(qiáng)度，即“吸收率”。要求根據(jù)模板及其接收信息，確定CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的位置、探測(cè)器單元之間的距離以及該CT系統(tǒng)使用的X射線的180個(gè)方向。

問(wèn)題2：給出利用上述CT系統(tǒng)得到的某未知介質(zhì)的接收信息。利用問(wèn)題1中得到的標(biāo)定參數(shù)，確定該未知介質(zhì)在正方形托盤中的位置、幾何形狀和吸收率等信息。另外，給出圖2所給的10個(gè)位置處的吸收率。

2 模型假設(shè)

1）假設(shè)將每個(gè)探測(cè)器視為一個(gè)接收點(diǎn)，其自身的寬度忽略不計(jì)。

2）假設(shè)X射線垂直于橢圓短軸透射樣品時(shí)在探測(cè)器上得到的投影數(shù)據(jù)最大。

圖1 模板示意圖（mm）

圖2 10個(gè)位置示意圖（mm）

3）假設(shè)X射線垂直于橢圓長(zhǎng)軸透射樣品時(shí)在探測(cè)器上得到的投影面積最大。

3 給定模板下的CT系統(tǒng)參數(shù)設(shè)定模型

3.1 探測(cè)器單元之間的距離

本文假設(shè)每個(gè)探測(cè)器單元看成一個(gè)接收點(diǎn)，且512個(gè)探測(cè)器等距排列。設(shè)所給模板中的圓形介質(zhì)直徑為d，接收到經(jīng)過(guò)圓形介質(zhì)的X射線的探測(cè)器個(gè)數(shù)為k，探測(cè)器單元之間的距離為l，則有：

根據(jù)180個(gè)不同方向上，512個(gè)探測(cè)器對(duì)模板的接收信息，選取第1個(gè)方向的數(shù)據(jù)為例，做出不同探測(cè)器得到的吸收率的圖像。易知，圖像第一段波動(dòng)曲線即為X射線經(jīng)過(guò)橢圓形介質(zhì)得出的吸收率，第二段波動(dòng)曲線即為X射線經(jīng)過(guò)圓形介質(zhì)得出的吸收率。由此，接收到經(jīng)過(guò)圓形介質(zhì)的X射線的探測(cè)器個(gè)數(shù)即為第二段波動(dòng)曲線在軸上的截距。由計(jì)算可讀取得出，接收到經(jīng)過(guò)圓形介質(zhì)的X射線的探測(cè)器為第402至第430個(gè)，共29個(gè)探測(cè)器。又由圖1模板示意圖可知，圓形介質(zhì)的直徑為8 mm。因此，探測(cè)器單元之間距離為0.275 8 mm.

3.2 CT系統(tǒng)旋轉(zhuǎn)中心在托盤中的位置

3.2.1 對(duì)于旋轉(zhuǎn)中心在托盤中心附近的證明

觀察所給每列數(shù)據(jù)，存在兩種情況：只有一個(gè)連續(xù)的正數(shù)群，被0值隔開(kāi)的兩個(gè)連續(xù)的正數(shù)群，可知無(wú)論是分開(kāi)還是聚合的正數(shù)群，在最中間第256和第257個(gè)探測(cè)器處的數(shù)值M附近存在距離ε個(gè)探測(cè)器數(shù)據(jù)L，有：

令橢圓中心與旋轉(zhuǎn)中心距離為d'，已知探測(cè)器間距d，有：

由此可得，旋轉(zhuǎn)中心位于正方形托盤中心的附近。

3.2.2 旋轉(zhuǎn)中心在托盤中的位置模型

以橢圓介質(zhì)中心為原點(diǎn)建立平面直角坐標(biāo)系，水平于正方形托盤的方向?yàn)閤軸且向右為正方向，垂直方向?yàn)閥軸且向上為正方向。為方便計(jì)算，本文選取X射線分別垂直于x軸與y軸的兩種情況，分別計(jì)算旋轉(zhuǎn)中心的橫坐標(biāo)與縱坐標(biāo)。

當(dāng)X射線垂直于y軸，即探測(cè)器旋轉(zhuǎn)于坐標(biāo)系左側(cè)時(shí)，此時(shí)可根據(jù)計(jì)算誤差值得出旋轉(zhuǎn)中心縱坐標(biāo)。由所給數(shù)據(jù)可得非零數(shù)值最多的一列，記為Cmax，表示在該方向上讀取到經(jīng)過(guò)介質(zhì)的X射線的探測(cè)器數(shù)量最多。設(shè)該列上最左端探測(cè)器和最右端探測(cè)器分別為第a個(gè)和第b個(gè)探測(cè)器。取兩端探測(cè)器平均值，根據(jù)平均值所在位置與原點(diǎn)差異即可得出旋轉(zhuǎn)中心縱坐標(biāo)y0，即：

當(dāng)X射線垂直于x軸，即探測(cè)器旋轉(zhuǎn)于坐標(biāo)系正上方時(shí)，此時(shí)可根據(jù)計(jì)算誤差值得出旋轉(zhuǎn)中心橫坐標(biāo)。由所給數(shù)據(jù)可得非零數(shù)值最少的一列，記為Cmin，表示在該方向上讀取到經(jīng)過(guò)介質(zhì)的X射線的探測(cè)器數(shù)量最少。設(shè)該列上最左端探測(cè)器和最右端探測(cè)器分別為第e個(gè)和第f個(gè)探測(cè)器。取兩端探測(cè)器平均值，根據(jù)平均值所在位置與原點(diǎn)差異即可得出旋轉(zhuǎn)中心橫坐標(biāo)x0，即：

3.2.3 模型求解

用MATLAB編程，算法如下：首先，找出非零數(shù)值最多的一列Cmax及最少的一列Cmin；其次，分別找出以上兩列的最左端探測(cè)器和最右端探測(cè)器；最后將數(shù)值代入式（4）（5）中求解。求得結(jié)果為Cmax=58，其最左端探測(cè)器和最右端探測(cè)器分別為第92個(gè)和第169個(gè)探測(cè)器；Cmin=150，其最左端探測(cè)器和最右端探測(cè)器分別為第276個(gè)和第380個(gè)探測(cè)器。旋轉(zhuǎn)中心坐標(biāo)為（-9.239 3，5.516 0）。

3.3 X射線的180個(gè)方向

對(duì)于X射線的180個(gè)方向，可依據(jù)特殊方向及每次的旋轉(zhuǎn)角度得出。當(dāng)X射線從橢圓介質(zhì)的短軸方向射入時(shí)，對(duì)應(yīng)數(shù)據(jù)中該列數(shù)據(jù)具有最高的對(duì)稱性，用MATLAB編程搜索可得為第61列，可知第61個(gè)方向上不同探測(cè)器得到的吸收率。同理，當(dāng)X射線從橢圓介質(zhì)的長(zhǎng)軸方向射入時(shí)，X射線穿過(guò)物體的厚度最大，物體對(duì)射線的吸收量最大，即吸收率最大。由所給數(shù)據(jù)可知吸收率最大的一列為第151列，可知第151個(gè)方向上不同探測(cè)器得到的吸收率。

從第61次旋轉(zhuǎn)到第151次，旋轉(zhuǎn)90°，則每次旋轉(zhuǎn)角度為1°。又第61次旋轉(zhuǎn)后X射線垂直y軸，由題意知從初始位置逆時(shí)針旋轉(zhuǎn)至第61個(gè)方向，中間轉(zhuǎn)過(guò)角度為61°。則初始入射位置與y軸夾角為29°。X射線180個(gè)方向結(jié)果見(jiàn)圖3所示。

圖3 X射線180個(gè)方向?qū)?yīng)結(jié)果圖

4 基于Radon變換重建的吸收率求解模型

4.1 X射線強(qiáng)度衰減與圖像重建的原理

當(dāng)X射線穿過(guò)均勻材料的介質(zhì)時(shí)，根據(jù)Lambert Beers定理，射線強(qiáng)度按指數(shù)規(guī)律衰減，滿足公式：

其中，I為穿過(guò)介質(zhì)后射線的強(qiáng)度，I0為入射前的射線強(qiáng)度，μ為衰減系數(shù)，即介質(zhì)對(duì)射線的吸收系數(shù)，x為介質(zhì)的厚度。

當(dāng)X射線穿過(guò)一個(gè)含有不同吸收系數(shù)材料的非均勻物體時(shí)，則有：

其中，μ（x，y）是沿L的線積分，即隨路徑L變化的衰減系數(shù)的函數(shù)。將函數(shù)μ（x，y）在平面上沿直線L的線積分定義為：

上式稱為 μ（x，y）的 Radon 變換，μ（x，y）關(guān)于某直線的Radon變換就是μ（x，y）沿該直線的一維投影。Radon變換的逆變換則對(duì)應(yīng)圖像二維函數(shù)的重建。

4.2 未知介質(zhì)的幾何形狀和吸收率

利用變換的逆變換，可編寫(xiě)程序得出問(wèn)題2中所給數(shù)據(jù)的未知介質(zhì)平面圖，未知介質(zhì)的幾何形狀為橢圓形，橢圓右側(cè)有兩個(gè)小橢圓形小洞。

用MATLAB可找出大橢圓形最左端、最右端、最上端及最下端的4個(gè)坐標(biāo)點(diǎn)分別為（6.06,44.90），（87.6，49.86），（48.48，25.62），（47.38，69.42）。將題目所給10個(gè)位置處轉(zhuǎn)化為橢圓形圖中對(duì)應(yīng)點(diǎn)的坐標(biāo)，分別找出對(duì)應(yīng)的10個(gè)吸收率，從左至右分別為0，0，0.485 1，0，0.497 7，0，0，0，0，0。