一種雙能X 線骨密度儀的測(cè)量算法研究

王月圓俞政濤趙興群*

(1.東南大學(xué)生物科學(xué)與醫(yī)學(xué)工程學(xué)院，江蘇南京 210096；2.南京澳思泰生物科技有限公司，江蘇南京 210046)

骨密度值是衡量骨質(zhì)量的重要指標(biāo)，對(duì)評(píng)價(jià)骨質(zhì)流失、診斷骨質(zhì)疏松并給予及時(shí)治療對(duì)治療骨質(zhì)疏松起著至關(guān)重要的作用[1]。目前中國50 歲以上的骨質(zhì)疏松人數(shù)預(yù)計(jì)有6 940 萬，骨質(zhì)疏松性骨折帶來高致死率和高致殘率。其中，髖部骨折是骨質(zhì)疏松骨折中后果最為嚴(yán)重的一種，骨折后一年患者的死亡率高達(dá)20%。全球每年因骨質(zhì)疏松導(dǎo)致髖部骨折的患者人數(shù)約166 萬人，而隨著人口老齡化日益突出，2050 年全球髖部骨折發(fā)生率將增長4倍，達(dá)到每年約700 萬人[2]。

人口老齡化問題，使得對(duì)骨密度精確測(cè)量的需求急劇增加。雙能X 線骨密度儀使用具有2 種不同能量的X 線，可在人體軟組織和骨骼中得到兩種不同能量的透過性，通過對(duì)其吸收特性和各種能的透過率的比較，求得人體骨密度。和傳統(tǒng)的單能和超聲骨密度儀相比，雙能X 線骨密度儀利用了人體骨骼和軟組織的吸收系數(shù)不同的特點(diǎn)，提高了骨密度測(cè)量結(jié)果的準(zhǔn)確性，是目前骨密度測(cè)定的金標(biāo)準(zhǔn)[3]。雙能X 線骨密度儀得到高能低能X 線圖像后，現(xiàn)有計(jì)算骨密度值的算法主要分為2 大類:理論公式法和擬合法。

理論公式計(jì)算法最早由Mazess[4]提出，被廣泛用于骨質(zhì)疏松的診斷以及后續(xù)治療，根據(jù)骨骼和軟組織對(duì)不同能量X 射線的衰減差異聯(lián)立方程組，求得人體骨密度值。但理論公式法在實(shí)際使用過程中由于高壓偏移、電流不穩(wěn)定、溫度變化等原因，在實(shí)際應(yīng)用中能譜分布測(cè)量也是非常困難的，因此存在較大誤差。針對(duì)這些問題，武伯歌[5]提出將廣譜的X 射線近似當(dāng)成單色光來處理，選擇沒有骨骼只有軟組織的部分進(jìn)行測(cè)量，進(jìn)行軟組織消除。這種方法在普通醫(yī)用X 線機(jī)獲得了較好的效果。張峰等[6]提出一個(gè)系統(tǒng)，這個(gè)系統(tǒng)產(chǎn)生的X 射線的質(zhì)量明顯高于由市電升壓、整流得到的高壓所產(chǎn)生的X射線，使得低能光子所占比例小，可被認(rèn)為是單色射線。陳敏聰[7]提出利用硬件修正的方法，對(duì)X 射線進(jìn)行預(yù)硬化，使其具有近似單能光子的衰減規(guī)律。盡管上述優(yōu)化方法一定程度上提高理論公式法在實(shí)際測(cè)量中的準(zhǔn)確性，但由于實(shí)際測(cè)量環(huán)境復(fù)雜，誤差不可避免，因此擬合法在實(shí)際測(cè)量過程中更為常用。

擬合法目前最常用的是Cardinal 等[8]提出的二次三次圓錐曲面擬合法，用二次三次圓錐曲面模擬骨密度和高能低能X 線之間的關(guān)系，該方法和傳統(tǒng)的多項(xiàng)式擬合法相比更符合X 線衰減規(guī)律，所需擬合數(shù)據(jù)較少，但也存在準(zhǔn)確度欠缺，擬合公式復(fù)雜，受環(huán)境影響較大等問題。

針對(duì)這些問題，在對(duì)現(xiàn)有骨密度測(cè)量算法進(jìn)行深入研究的基礎(chǔ)上，提出了一種基于對(duì)數(shù)擬合的改進(jìn)算法，引入高能低能圖像背景灰度值，從而減小不同測(cè)量環(huán)境帶來的誤差，提高骨密度測(cè)量的準(zhǔn)確度。通過實(shí)驗(yàn)，和現(xiàn)有算法進(jìn)行比較，驗(yàn)證了該算法的穩(wěn)定性和準(zhǔn)確度。

1 骨密度測(cè)量算法

1.1 現(xiàn)有算法

現(xiàn)有的針對(duì)雙能X 線骨密度儀的骨密度測(cè)量算法主要分為2 大類:理論公式法和擬合法。

理論公式計(jì)算法假設(shè)高能和低能的X 射線均為理想單能射線，則其穿過組織之前的高低能射線強(qiáng)度IOH和IOL與穿過介質(zhì)之后的射線強(qiáng)度IH和IL滿足朗伯比爾衰減規(guī)律:

上述方程組中，ubh、ubl分別表示骨組織在高低能下的質(zhì)量衰減系數(shù)；ush、usl表示軟組織在高低能下的質(zhì)量衰減系數(shù)，單位是cm2/g；Ms、Mb表示軟組織和骨骼的面密度，單位是g/cm2。從定義中不難發(fā)現(xiàn)，Mb即為要求的骨密度值。方程組變形可以得到:

式中的值Rs、Rb表示對(duì)于某種介質(zhì)的高低能量線性衰減系數(shù)之比Rs＝usl/ush和Rb＝ubl/ubh，是介質(zhì)與光子發(fā)生作用的光電-康普頓比。當(dāng)介質(zhì)已知時(shí)，R值可以通過查表得到。這種骨密度計(jì)算方法雖然符合X 射線衰減規(guī)律，但是，由于實(shí)際得到的高能低能X射線能譜復(fù)雜，且會(huì)受到高壓偏移、電流不穩(wěn)定、溫度變化等影響，在實(shí)際工程中，無法得到較為準(zhǔn)確的高能低能X 線衰減系數(shù)等參數(shù)，基于理論公式的骨密度計(jì)算方法無法得到更廣泛的應(yīng)用。

為了解決理論公式法穩(wěn)定性和準(zhǔn)確性的問題，提出了擬合法，擬定擬合公式，并通過大量實(shí)驗(yàn)，確定擬合公式中的系數(shù)，來模擬骨密度值和高能低能X 線透過值之間的規(guī)律。已有的使用比較廣的是多項(xiàng)式擬合法和二次三次圓錐曲面擬合法。

多項(xiàng)式擬合法是雙能X 線骨密度測(cè)量中比較常用的一種擬合方法，用多項(xiàng)式方程模擬骨密度值和高能低能X 線的之間的關(guān)系。x、y分別代表高能、低能X 線的透過值，f(x，y)表示骨密度值，擬合關(guān)系可以用式(5)表示:

用最小二乘法擬合求解出這個(gè)pij的方程組，可以得到多項(xiàng)式擬合的解析式，從而得到骨密度值和高能低能X 線透過值之間的關(guān)系。

二次三次圓錐曲面擬合用二次或三次圓錐曲面公式代替?zhèn)鹘y(tǒng)多項(xiàng)式曲面，以改進(jìn)傳統(tǒng)多項(xiàng)式曲面擬合對(duì)數(shù)據(jù)量要求較多的缺點(diǎn)，同時(shí)提高計(jì)算準(zhǔn)確度。骨密度值F滿足圓錐曲面方程:

式中:A、B、C為高能、低能X 線的透過值x、y的二次或三次泰勒展開。根據(jù)求根公式可以得到式(6)的解。用迭代非線性最小二乘擬合可以求得擬合表達(dá)式中的系數(shù)，從而得到骨密度值和高能低能X 線透過值之間的關(guān)系。

1.2 基于對(duì)數(shù)擬合的改進(jìn)算法

上述2 種擬合算法，多項(xiàng)式擬合法被廣泛地應(yīng)用在各個(gè)領(lǐng)域，但擬合公式并不能反映X射線的衰減規(guī)律，因此，準(zhǔn)確度和魯棒性存在一定的問題，達(dá)到同樣精度所需的數(shù)據(jù)量也相對(duì)較多；圓錐曲面擬合法針對(duì)多項(xiàng)式擬合法這些問題進(jìn)行了改進(jìn)，但是擬合公式比較復(fù)雜，使公式有理化的過程也會(huì)帶來一定的誤差，影響測(cè)量準(zhǔn)確度。針對(duì)這2 種擬合算法存在的問題，同時(shí)結(jié)合骨密度計(jì)算的理論公式，提出了一種基于對(duì)數(shù)擬合的改進(jìn)算法。

根據(jù)式(3)的骨密度理論計(jì)算公式，骨密度值Mb和高能X 線出射光強(qiáng)測(cè)量值與高能X線入射光強(qiáng)測(cè)量值的對(duì)數(shù)比、低能X 線出射光強(qiáng)測(cè)量值與低能X線入射強(qiáng)度測(cè)量值的對(duì)數(shù)比有關(guān)，而其中的系數(shù)Rs、ubl、ubh均與高能低能X 射線源的能量值有關(guān)，理論公式計(jì)算時(shí)，這些參數(shù)均通過查表獲得，但考慮到實(shí)際使用的X 射線源并不是理想的單能X 射線，通過查表無法獲得準(zhǔn)確的參數(shù)值，從而在計(jì)算時(shí)帶來極大的誤差?；趯?duì)數(shù)擬合的改進(jìn)算法就是針對(duì)這一問題進(jìn)行改進(jìn)，式(3)可改寫為:

式中:Rs＝usl/ush，可以看出和前的系數(shù)只與高能低能X 射線對(duì)骨頭和軟組織的吸收系數(shù)相關(guān)，無論X 射線是否為理想的單能，系數(shù)最終都是2 個(gè)常量，只隨X 射線源能量的變化而變化。因此，可以將系數(shù)看成一個(gè)整體，得到以下擬合公式:

用x、y分別代表高能、低能X線的透過值，f(x，y)表示骨密度值，擬合關(guān)系可以表示為:

其中穿過組織之前的高低能射線強(qiáng)度IOH和IOL在X 射線源確定的情況下，仍會(huì)隨著測(cè)量環(huán)境變化，例如測(cè)量時(shí)的電壓變化而變化。因此每次更換測(cè)量環(huán)境，均需對(duì)IOH、IOL進(jìn)行校正，以提高骨密度測(cè)量準(zhǔn)確度，故不可將IOH、IOL看作不可變系數(shù)進(jìn)行合并。根據(jù)擬合關(guān)系式進(jìn)行線性最小二乘擬合，假設(shè)實(shí)驗(yàn)得到m(m＞2)組高能低能X 線透過值和骨密度值數(shù)據(jù)，將和看作整體參數(shù)X1和X2，則有超定方程組:

引入殘差平方和函數(shù)S:

通過對(duì)S(β)進(jìn)行微分求最小值，可以得到:

如果矩陣XTX非奇異，則有唯一解:

即可得到系數(shù)a、b的最優(yōu)解，得到解析式，也就是骨密度值和高能低能X 線透過值之間的關(guān)系。圖1 以系數(shù)a、b為3，IOH為3 000，IOL為10 000 為例，展示了基于對(duì)數(shù)擬合的改進(jìn)算法中，骨密度值和高能低能X 線透過值之間的關(guān)系。

圖1 a＝3、b＝3 時(shí)高能低能X 線透過值和骨密度值之間的關(guān)系

2 實(shí)驗(yàn)驗(yàn)證

2.1 實(shí)驗(yàn)環(huán)境

實(shí)驗(yàn)中用雙能X 線骨密度儀采集圖像，其中X射線球管采用的是萬東電子的XD51-6、20/125 型號(hào)球管，該球管的最高工作電壓為125 kV，焦點(diǎn)大小可以達(dá)到0.5 mm。雙能X 射線產(chǎn)生選擇的是K邊緣過濾法，通過對(duì)X 射線球管施加80 kV 的高壓，使其產(chǎn)生能量為100 keV 的連續(xù)X 線譜，然后使用稀土元素過濾射線，得到能量峰在40 keV 和70 keV的高低能射線。

2.2 實(shí)驗(yàn)?zāi)ｓw

實(shí)驗(yàn)中作為擬合的標(biāo)準(zhǔn)數(shù)據(jù)，使用的是CIRS公司生產(chǎn)的專門為DEXA 系統(tǒng)設(shè)計(jì)的BFP(Bona Fide Phantom)脊椎模體。該模體長22 cm 寬19 cm高15 cm，覆蓋了從0.7 g/cm2到1.5 g/cm2的常見骨密度范圍。材料為鈣羥基磷灰石，該模體模擬脊柱正位，弧度邊緣對(duì)邊緣檢測(cè)算法要求高。這種成熟的模體用于校準(zhǔn)機(jī)器骨密度測(cè)量的精度和測(cè)試算法效果。

圖2 BFP 模體實(shí)物圖

BFP 模體主要的4 個(gè)部分質(zhì)量厚度如表1所示。

表1 BFP 模體各部分質(zhì)量厚度

2.2 測(cè)試模體

由于鋁和有機(jī)玻璃在40keV 和70keV 附近的衰減系數(shù)和人體骨骼和軟組織的衰減系數(shù)相近，因此用不同厚度的鋁片和有機(jī)玻璃片進(jìn)行組合模擬人體骨密度范圍，作為測(cè)試模體[7]，質(zhì)量厚度0.8 g/cm2到2.0 g/cm2。同時(shí)，為了驗(yàn)證擬合公式是否可以有效去除軟組織對(duì)骨密度測(cè)量的影響，選擇鋁片數(shù)量相同時(shí)，增減有機(jī)玻璃片的數(shù)量，觀察計(jì)算骨密度計(jì)算結(jié)果是否有明顯變化。實(shí)驗(yàn)中用到的不同組合的質(zhì)量厚度如表2 所示。

表2 測(cè)試模體質(zhì)量厚度

2.3 實(shí)驗(yàn)過程與結(jié)果分析

實(shí)驗(yàn)中，首先進(jìn)行空拍，得到高能低能X 線的初始值，接著掃描標(biāo)準(zhǔn)BFP 模體，得到標(biāo)準(zhǔn)BFP 模體高能低能X 線透射圖像和對(duì)應(yīng)的透射值。透射圖像如圖3 所示。

圖3 BFP 模體原始投影圖

獲得高能低能X 線圖像后，接著對(duì)圖像進(jìn)行分割，得到感興趣的區(qū)域。再對(duì)分割得到的圖像進(jìn)行必要的濾波處理，得到較為理想的高能低能X 圖像。最后對(duì)各個(gè)感興趣區(qū)域的灰度求平均值，得到各部分對(duì)應(yīng)的透射值。具體處理流程如圖4 所示。

圖4 圖像處理基本流程

表3 為實(shí)際得到的各部分高能低能X 線透射值。

表3 BFP 模體的高能低能X 線透射值

接著掃描測(cè)試模體，得到不同層數(shù)對(duì)應(yīng)的高能低能X 線透射圖，同樣對(duì)圖像進(jìn)行處理，最終得到的不同層數(shù)對(duì)應(yīng)的透射值如表4 所示。

表4 測(cè)試模體的高能低能X 線透射值

然后用BFP 模體作為擬合數(shù)據(jù)，測(cè)試模體作為測(cè)試數(shù)據(jù)進(jìn)行擬合。在實(shí)際計(jì)算過程中，可以看出，當(dāng)標(biāo)準(zhǔn)數(shù)據(jù)較少時(shí)，多項(xiàng)式擬合和圓錐曲面擬合的擬合關(guān)系式的次數(shù)僅能達(dá)到一次，而這2 種方法擬合的準(zhǔn)確度與擬合關(guān)系式的次數(shù)相關(guān)性較高，導(dǎo)致只用一次擬合誤差較大，最大相對(duì)誤差超過骨密度儀國家標(biāo)準(zhǔn)誤差范圍10%，因此不做過多展示?；趯?duì)數(shù)擬合的改進(jìn)算法得到的解析式中，對(duì)應(yīng)式(9)中的參數(shù)a、b分別為3.696 2 和3.153 8。

將測(cè)試模體的高能低能X 線透射值代入解析式，得到改進(jìn)算法對(duì)質(zhì)量厚度的擬合結(jié)果以及對(duì)應(yīng)的相對(duì)誤差如表5 所示。

表5 測(cè)試模體的擬合結(jié)果

鋁片數(shù)量相同時(shí)，對(duì)數(shù)擬合的計(jì)算結(jié)果的最大誤差如表6 所示。

表6 鋁片數(shù)量相同時(shí)的最大測(cè)量誤差

最后用測(cè)試模體作為擬合數(shù)據(jù)，BFP 模體作為測(cè)試數(shù)據(jù)進(jìn)行交叉驗(yàn)證?；趯?duì)數(shù)擬合的改進(jìn)算法得到的解析式中，對(duì)應(yīng)式(9)中的參數(shù)a、b分別為3.502 2 和2.999 7。將BFP 模體的高能低能X 線透射值代入解析式，改進(jìn)算法對(duì)質(zhì)量厚度的擬合結(jié)果以及對(duì)應(yīng)的相對(duì)誤差如表7 所示。

表7 BFP 模體擬合結(jié)果

從實(shí)驗(yàn)結(jié)果可以看出，對(duì)數(shù)擬合的測(cè)量結(jié)果和理論值的平均相對(duì)誤差2.53%，最大相對(duì)誤差小于8%，符合國家標(biāo)準(zhǔn)誤差范圍，測(cè)量結(jié)果較為準(zhǔn)確。同時(shí)在鋁片數(shù)量相同，而有機(jī)玻璃數(shù)量不同的情況下，測(cè)量的誤差在2%以內(nèi)，比較好地排除了有機(jī)玻璃對(duì)測(cè)量結(jié)果產(chǎn)生的影響。

3 結(jié)論

雙能X 線骨密度儀測(cè)量骨密度的算法主要分為理論公式法和擬合法。理論公式法在實(shí)際使用過程中存在較大誤差，因此擬合法在實(shí)際測(cè)量過程中更為常用。常見的擬合方法有多項(xiàng)式擬合、二次三次圓錐曲面擬合2 種。二次三次圓錐曲面擬合，在多項(xiàng)式擬合的基礎(chǔ)上進(jìn)行改進(jìn)，減少擬合所需數(shù)據(jù)量的同時(shí)提高了準(zhǔn)確度。

在對(duì)現(xiàn)有骨密度算法進(jìn)行深入研究的基礎(chǔ)上，提出了基于對(duì)數(shù)擬合的改進(jìn)算法，擬合公式簡單，進(jìn)一步減少了擬合對(duì)數(shù)據(jù)量的要求，同時(shí)也更加符合骨密度值與X 線透射值之間的理論關(guān)系。實(shí)際實(shí)驗(yàn)過程中，測(cè)量平均相對(duì)誤差2.53%，最大相對(duì)誤差小于8%，算法準(zhǔn)確度較高。鋁片數(shù)量相同，有機(jī)玻璃數(shù)量不同的情況下，測(cè)量的誤差小于2%，比較好地排除了有機(jī)玻璃對(duì)測(cè)量結(jié)果產(chǎn)生的影響，可以推測(cè)在實(shí)際測(cè)量過程中可以較好地排除軟組織對(duì)骨密度測(cè)量的影響，進(jìn)一步證明了該算法的準(zhǔn)確度。