基于頻譜成像與最大似然估計的骨密度測量方法
【專利摘要】本發(fā)明屬于超聲醫(yī)學定量測量技術領域，公開了一種骨骼超聲系統(tǒng)中從背散射信號提取出骨小梁平均間距的方法。該方法將超聲系統(tǒng)中獲取的松質骨背散射信號作為分析對象，首先將多次測量獲得的背散射信號通過濾波、截斷、快速傅里葉變換等步驟后，組合成頻域的圖像；然后通過最大似然估計得到頻域中每一列的基頻估計；最后由線性回歸做擬合得到關于基頻的函數(shù)。由基頻值和已知的超聲傳播速度推出骨小梁平均間距。本發(fā)明首次提出了利用頻譜成像與最大似然估計的方法來研究平均骨小梁間距，相比于傳統(tǒng)的方法排除了單次測量的偶然性，同時對于強噪聲環(huán)境有較強的魯棒性。
【專利說明】
基于頻譜成像與最大似然估計的骨密度測量方法
技術領域
[0001] 本發(fā)明涉及人體骨密度測量問題，屬于超聲醫(yī)學定量測量技術領域。本發(fā)明結合 信號頻譜成像與最大似然估計，從統(tǒng)計的角度提出了一種新的從背散射信號中得到骨小梁 平均間距的方法。
【背景技術】
[0002] 隨著人口趨于老齡化，多種負面的問題接踵而來，骨質疏松癥及其并發(fā)骨折就是 其中之一。據(jù)不完全統(tǒng)計，目前全世界骨質疏松癥患者約有2億人，而我國的骨質疏松癥患 者也已經超過了 9000萬人，約占我國總人口的7%，其中大約有50%的老年婦女和20%的老 年男性患有骨質疏松癥，骨質疏松癥已經像心臟病、癌癥、糖尿病一樣嚴重威脅著中老年人 的健康。目前世界衛(wèi)生組織已將骨質疏松癥與糖尿病、心血管病共同列為危害中老年人健 康的三大殺手。鑒于我國即將全面進入老齡化社會，對于骨密度檢測的研究有重大的社會 價值和經濟價值。
[0003] 骨質可分為皮質骨和松質骨兩個部分，皮質骨是骨的外層結構，由位于骨膜下的 密質骨組成。松質骨為海綿狀，由許多針狀或者片狀的骨小梁相互交織而成，內部由骨髓填 充。骨質疏松癥患者的松質骨骨小梁數(shù)量明顯減少，骨小梁變薄，同時骨小梁間距明顯增 大。松質骨的超聲背散射信號中包含著骨小梁的微觀結構信息。
[0004] 松質骨中骨小梁成網狀排列，超聲波在通過時會向不同方向形成散射波，而其中 沿著射入路徑返回并被換能器接收到的信號為背散射信號。背散射信號本身產生的機理決 定了該信號的幅度較小，信噪比較低。這些因素都給平均骨小梁間距(MTBS)的估計帶來了 困難。當前從背散射信號中提取出骨小梁平均間距的方法主要有AR倒譜法、譜自相關算法 等。然而這些傳統(tǒng)的算法都還存在著缺陷：一方面只對單次測量的結果做分析，無法排除測 量過程中的偶然性；另一方面?zhèn)鹘y(tǒng)方法對于噪聲的敏感度較高，往往對于信號基頻出現(xiàn)錯 誤判斷。
[0005] 傳統(tǒng)方法對于信號的處理方式是將每一次測量和分析看成一個獨立的過程，然而 由于諸多因素的影響例如信號質量不理想、超聲聚焦位置不準確或者算法自身對噪聲敏 感，都會帶來錯誤的結果。采用頻譜圖和最大似然估計的方法，對多次的測量結果分組分 析，用統(tǒng)計的方法避免了偶然性帶來的錯誤，同時結果在頻譜圖中有更直觀的呈現(xiàn)。

【發(fā)明內容】

[0006] 由于松質骨中骨小梁在壓力方向上近似規(guī)則分布，它的背散射信號中包含周期信 號。本發(fā)明提供了一種基于頻譜圖和最大似然估計的方法獲取上述周期信號的頻率。該方 法包括信號提取、信號重組與頻譜成像、最大似然估計分析、結果擬合四個部分。
[0007] 超聲回波包含諸多干擾波，如皮質骨反射信號、隨機背散射信號、白噪聲等。超聲 在骨組織中傳播時，首先在皮質骨發(fā)生反射，再經過松質骨并發(fā)生散射，這些信號先后由超 聲換能器接收。我們感興趣的是松質骨背散射信號，因此可以在時間域中選取合適的時間 窗提取出松質骨背散射信號，并對提取出的信號做低通濾波。
[0008] 最大似然估計需要足夠的樣本才能得到可信的結果，因此對同一部位做多次測 量，將測量結果通過濾波、截斷、快速傅里葉變換得到松質骨的背散射頻域信號。將所得到 的頻域信號按照測量次序放入矩陣中，作為樣本空間。
[0009] 信號頻域的數(shù)值矩陣包含著基頻頻率，下面采用最大似然估計將基頻提取出來。 對于矩陣的各列信號可以由下式表示：
[0010]
[0011] 可以看到將信號表示為多次諧波線性疊加與噪聲v[t]之和。為了求取fn，將上式 表示為矩陣形式，如式7所示：
[0012] A(fn)xn+v[n] = s[n] (8)
[00123]其中A(fn)為N*2L的矩陣，第j行為：
[0014]
[0015] 假設v[n]為方差未知，期望為0的高斯分布的噪聲信號。構造關于參數(shù)向量χη的似 然函數(shù)：
[0016]
[0017]對似然函數(shù)f(Xn)取對數(shù)可得：
[0019] 由于高斯噪聲的方差未知，用以下運算求出估計量〇2:
[0018]
[0020]
[0021]
[0022] 不妨 <，將式（13)代入式（10)可得： 2 /=U
[0023] f(xn) = (2JiV2/N)-N/2exp(_N/2) (14)
[0024] 由式(14)可知，使得f(Xn)為極大值的相當于使V2最小，SP: 3
[0025]
[0026] 對于上式直接做最大似然估計，運算量較大。將= ,·?[/?]代入式（15)中， 得到關于fη的壓縮似然估計(compressed likelihood)表達式：
[0027]
[0028] 其中/= (7-)。令PA(fn) = A(f)使得式（16)為極小值的f"相當于使 |Α(Λ)?取極大值。為進一步減少運算量，將A(fn)作QR分解，即△(/") = (?/"~，取矩陣Q的 前2L列，則fnSI以由式（17)得到：
[0029]
(17)
[0030] 假設實際基頻值fQ，對矩陣的每N個樣本做最大似然估計得到的多個基頻值應服 從期望為f〇的正態(tài)分布。為了排除由各種因素引發(fā)的錯誤估計結果，首先將所得一系列基 頻值做平均，并求出每個基頻相對于平均值的偏差，將偏差值較大的數(shù)值取代以平均值;然 后在此基礎上做線性回歸，得到準確的基頻估計結果。
【附圖說明】
[0031] 圖1是本發(fā)明的流程示意框圖。
[0032]圖2是本發(fā)明采用的仿真信號。
[0033] 圖3是本發(fā)明最大似然估計算法框圖。
[0034] 圖4是本發(fā)明得到的頻譜灰度圖和算法結果。
【具體實施方式】
[0035] 下面結合附圖，以從一段仿真信號中提取平均骨小梁間距為例，介紹本發(fā)明的實 施過程。本發(fā)明的流程如圖1所示，包括四個部分:信號預處理、信號頻譜成像、最大似然估 計、線性回歸。對仿真信號做時域截斷，得到松質骨部分的超聲背散射信號。如圖2所示，采 用長度為1200的窗截取背散射信號，圖中紅色的框內部代表信號中被截取的部分。從矩陣 包含樣本中取出連續(xù)N個樣本做最大似然估計，流程如圖3所示，其中Nmx表示矩陣中最大樣 本數(shù)即矩陣的寬度。
[0036] 將所得到的基頻值取平均，設定閾值，將基頻值與平均值相差大于閾值的點去除， 并用平均值代替。最后，對基頻點做線性回歸，得到基頻在頻譜圖中的圖像。圖4是由仿真信 號得到的頻譜圖，寬表示測試的次數(shù)，高表示頻率值，每一個紅色的圓圈表示該樣本區(qū)間符 合式（17)的基頻點。由圖4可知，由于強噪聲的影響出現(xiàn)了錯誤估計，然而整體而言依然可 以看到正確的基頻分布位置。得到基頻估計值后，結合已知的超聲聲速，根據(jù)式（18)得到骨 小梁的平均間距，
[0037]
(18)
[0038] 式中c為超聲聲速，fo為基頻值。
【主權項】
1. 一種基于頻譜圖與最大似然估計的骨密度測量方法，其特征在于具體步驟為： (1) 對采集到的超聲回波信號采用窗函數(shù)截取松質骨部分并做低通濾波，得到松質骨 背散射信號； (2) 多次測量后，將信號做快速傅里葉變換并按照時間順序組合到信號矩陣中。將該矩 陣的成員值映射到[〇,255]灰度區(qū)間內并顯示成灰度圖像； (3) 觀測信號可表示為多次諧波線性疊加后與噪聲之和，諧波頻率為基頻的整數(shù)倍。將 基頻值和諧波系數(shù)看作待估計參數(shù)，頻譜圖像視為觀察到的總樣本空間。將頻譜圖像分割 為等大小的連續(xù)樣本區(qū)間，由最大似然估計得到信號圖像在各區(qū)間的基頻； (4) 排除奇異點，通過線性回歸得到所有基頻點的擬合函數(shù)。2. 根據(jù)權利要求1所述的方法，其特征是在于步驟（1)中，首先對信號做低通濾波，然后 用窗函數(shù)截取超聲背散射信號。3. 根據(jù)權利要求1所述的方法，其特征在于步驟(2)中，對于同一部位做多次測量。對每 次測量結果按照步驟(1)截取出松質骨背散射信號，并分別作快速傅里葉變換，然后將結果 按照時間先后放入信號矩陣中。將信號矩陣中的每一個值映射到[〇,255]的灰度區(qū)間內，得 到灰度圖像。4. 根據(jù)權利要求1所述的方法，其特征在于步驟(3)中，從信號矩陣內取出連續(xù)的N列的 樣本。骨小梁沿著壓力的方向呈近似規(guī)則分布，因此在背散射信號中包含著周期性的信號。 將觀測到的樣本表示為多次諧波的線性疊加與噪聲之和，如式(1):其中s為離散化后的信號;n表示諧波階數(shù)，階數(shù)越高對信號的細節(jié)復原越完整，然而要 求的數(shù)據(jù)量越多;an、bn代表諧波的系數(shù)，fn為基頻；。用L表示諧波的階數(shù)，上式的矩陣形式 可以表示為： A(fn)xn+v[ri] = s[n] (2) xn為諧波系數(shù)組成的向量，表示為： Xn- [Ell SL2 SL3 ... SLL bl \)2 ... t)L] (3) A(fn)為N*2L的矩陣，其中第j行為：將fn以及Xn看作待估計參數(shù)，對矩陣中的N個樣本做最大似然估計，公式如下：5.根據(jù)權利要求1所述的方法，其特征在于步驟(4)中，首先求得一系列基頻值的均值， 設定區(qū)間排除錯誤估計點；然后對于基頻做線性回歸，從而在頻譜圖上得到基頻的擬合函 數(shù)。該函數(shù)與頻率軸的交點即基頻值，最后由下式得出骨小梁平均間距：其中c為超聲在松質骨中傳播的速度，fo為基頻，MTBS為平均骨小梁間距。
【文檔編號】A61B8/00GK106037794SQ201610325585
【公開日】2016年10月26日
【申請日】2016年5月12日
【發(fā)明人】何愛軍, 唐彬, 唐國揚, 付思東, 陳自立, 孫琴, 湯冬冬, 李 昊
【申請人】南京大學