王月圓俞政濤趙興群*

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

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

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

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

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

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

1 骨密度測量算法

1.1 現(xiàn)有算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖1 a＝3、b＝3 時高能低能X 線透過值和骨密度值之間的關系

2 實驗驗證

2.1 實驗環(huán)境

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

2.2 實驗模體

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

圖2 BFP 模體實物圖

BFP 模體主要的4 個部分質量厚度如表1所示。

表1 BFP 模體各部分質量厚度

2.2 測試模體

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

表2 測試模體質量厚度

2.3 實驗過程與結果分析

實驗中，首先進行空拍，得到高能低能X 線的初始值，接著掃描標準BFP 模體，得到標準BFP 模體高能低能X 線透射圖像和對應的透射值。透射圖像如圖3 所示。

圖3 BFP 模體原始投影圖

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

圖4 圖像處理基本流程

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

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

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

表4 測試模體的高能低能X 線透射值

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

將測試模體的高能低能X 線透射值代入解析式，得到改進算法對質量厚度的擬合結果以及對應的相對誤差如表5 所示。

表5 測試模體的擬合結果

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

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

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

表7 BFP 模體擬合結果

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

3 結論

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

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