王 琪

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

1 問題重述

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

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

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

2 模型假設(shè)

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

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

圖1 模板示意圖（mm）

圖2 10個位置示意圖（mm）

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

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

3.1 探測器單元之間的距離

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

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

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

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

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

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

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

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

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

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

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

3.2.3 模型求解

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

3.3 X射線的180個方向

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

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

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

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

4.1 X射線強度衰減與圖像重建的原理

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

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

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

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

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

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

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

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