CN106037794A - 基于频谱成像与最大似然估计的骨密度测量方法 - Google Patents

Abstract

本发明属于超声医学定量测量技术领域，公开了一种骨骼超声系统中从背散射信号提取出骨小梁平均间距的方法。该方法将超声系统中获取的松质骨背散射信号作为分析对象，首先将多次测量获得的背散射信号通过滤波、截断、快速傅里叶变换等步骤后，组合成频域的图像；然后通过最大似然估计得到频域中每一列的基频估计；最后由线性回归做拟合得到关于基频的函数。由基频值和已知的超声传播速度推出骨小梁平均间距。本发明首次提出了利用频谱成像与最大似然估计的方法来研究平均骨小梁间距，相比于传统的方法排除了单次测量的偶然性，同时对于强噪声环境有较强的鲁棒性。

Description

基于频谱成像与最大似然估计的骨密度测量方法

技术领域

本发明涉及人体骨密度测量问题，属于超声医学定量测量技术领域。本发明结合信号频谱成像与最大似然估计，从统计的角度提出了一种新的从背散射信号中得到骨小梁平均间距的方法。

背景技术

随着人口趋于老龄化，多种负面的问题接踵而来，骨质疏松症及其并发骨折就是其中之一。据不完全统计，目前全世界骨质疏松症患者约有2亿人，而我国的骨质疏松症患者也已经超过了9000万人，约占我国总人口的7％，其中大约有50％的老年妇女和20％的老年男性患有骨质疏松症，骨质疏松症已经像心脏病、癌症、糖尿病一样严重威胁着中老年人的健康。目前世界卫生组织已将骨质疏松症与糖尿病、心血管病共同列为危害中老年人健康的三大杀手。鉴于我国即将全面进入老龄化社会，对于骨密度检测的研究有重大的社会价值和经济价值。

骨质可分为皮质骨和松质骨两个部分，皮质骨是骨的外层结构，由位于骨膜下的密质骨组成。松质骨为海绵状，由许多针状或者片状的骨小梁相互交织而成，内部由骨髓填充。骨质疏松症患者的松质骨骨小梁数量明显减少，骨小梁变薄，同时骨小梁间距明显增大。松质骨的超声背散射信号中包含着骨小梁的微观结构信息。

松质骨中骨小梁成网状排列，超声波在通过时会向不同方向形成散射波，而其中沿着射入路径返回并被换能器接收到的信号为背散射信号。背散射信号本身产生的机理决定了该信号的幅度较小，信噪比较低。这些因素都给平均骨小梁间距(MTBS)的估计带来了困难。当前从背散射信号中提取出骨小梁平均间距的方法主要有AR倒谱法、谱自相关算法等。然而这些传统的算法都还存在着缺陷：一方面只对单次测量的结果做分析，无法排除测量过程中的偶然性；另一方面传统方法对于噪声的敏感度较高，往往对于信号基频出现错误判断。

传统方法对于信号的处理方式是将每一次测量和分析看成一个独立的过程，然而由于诸多因素的影响例如信号质量不理想、超声聚焦位置不准确或者算法自身对噪声敏感，都会带来错误的结果。采用频谱图和最大似然估计的方法，对多次的测量结果分组分析，用统计的方法避免了偶然性带来的错误，同时结果在频谱图中有更直观的呈现。

发明内容

由于松质骨中骨小梁在压力方向上近似规则分布，它的背散射信号中包含周期信号。本发明提供了一种基于频谱图和最大似然估计的方法获取上述周期信号的频率。该方法包括信号提取、信号重组与频谱成像、最大似然估计分析、结果拟合四个部分。

超声回波包含诸多干扰波，如皮质骨反射信号、随机背散射信号、白噪声等。超声在骨组织中传播时，首先在皮质骨发生反射，再经过松质骨并发生散射，这些信号先后由超声换能器接收。我们感兴趣的是松质骨背散射信号，因此可以在时间域中选取合适的时间窗提取出松质骨背散射信号，并对提取出的信号做低通滤波。

最大似然估计需要足够的样本才能得到可信的结果，因此对同一部位做多次测量，将测量结果通过滤波、截断、快速傅里叶变换得到松质骨的背散射频域信号。将所得到的频域信号按照测量次序放入矩阵中，作为样本空间。

信号频域的数值矩阵包含着基频频率，下面采用最大似然估计将基频提取出来。对于矩阵的各列信号可以由下式表示：

$s\[t\]=\underset{\mathrm{\η}=1}{\overset{N_{\mathrm{\η}}}{\mathrm{\Σ}}}\left(a_{\mathrm{\η}}cos\right(2{\mathrm{\π\ηf}}_{n}t)+b_{\mathrm{\η}}sin(2{\mathrm{\π\ηf}}_{n}t\left)\right)+\mathrm{\ν}\[t\],t=n,n+1,\mathrm{....},n+N-1---\left(7\right)$

可以看到将信号表示为多次谐波线性叠加与噪声v[t]之和。为了求取f_n，将上式表示为矩阵形式，如式7所示：

A(f_n)x_n+v[n]＝s[n] (8)

其中A(f_n)为N*2L的矩阵，第j行为：

${\left[\begin{array}{c}cos\left(2{\mathrm{\πf}}_{n}j\right)\\ cos\left(2\mathrm{\π}2f_{n}j\right)\\ .\\ .\\ .\\ cos\left(2{\mathrm{\πLf}}_{n}j\right)\\ sin\left(2{\mathrm{\πf}}_{n}j\right)\\ \sin \left(2\mathrm{\π}2f_{n}j\right)\\ .\\ .\\ .\\ \sin \left(2{\mathrm{\πLf}}_{n}j\right)\end{array}\right]}^T---\left(9\right)$

假设v[n]为方差未知，期望为0的高斯分布的噪声信号。构造关于参数向量x_n的似然函数：

$f\left(x_n\right)={\left(2{\mathrm{\π\σ}}^2\right)}^{-N/2}\exp \[-\frac{1}{2{\mathrm{\σ}}^2}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{\left(A\right(f_n)x_n-s\[n\])}^2\]---\left(10\right)$

对似然函数f(x_n)取对数可得：

$\ln f\left(x_n\right)=-1/2ln2{\mathrm{\π\σ}}^2-\frac{1}{2{\mathrm{\σ}}^2}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{\left(A\right(f_n)x_n-s\[n\])}^2---\left(11\right)$

由于高斯噪声的方差未知，用以下运算求出估计量σ²：

$\frac{\∂}{\∂{\mathrm{\σ}}^2}\ln f\left(x_n\right)=0---\left(12\right)$

${\mathrm{\σ}}^2=\frac{1}{N}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{\left(A\right(f_n)x_n-s\[n\])}^2---\left(13\right)$

不妨令将式(13)代入式(10)可得：

f(x_n)＝(2πV²/N)^-N/2exp(-N/2) (14)

由式(14)可知，使得f(x_n)为极大值的相当于使V²最小，即：

$J(f_n,x_n)=min\left(\right|\left|A\right(f_n)x_n-s\[n\]|{|}_2^2)---\left(15\right)$

对于上式直接做最大似然估计，运算量较大。将代入式(15)中，得到关于f_n的压缩似然估计(compressed likelihood)表达式：

$min\left(J\right(f_n\left)\right)=min\left|\right|P_{A\left(f_n\right)}^{\⊥}s\[n\]|{|}_2^2---\left(16\right)$

其中令P_A(f_n)＝A(f)使得式(16)为极小值的f_n相当于使取极大值。为进一步减少运算量，将A(f_n)作QR分解，即取矩阵Q的前2L列，则f_n可以由式(17)得到：

$f_n^{\′}={\mathrm{argmax}}_{f_n}\left|\right|Q_{f_n}(:,1:2L)s\[n\]|{|}_2^2---\left(17\right)$

假设实际基频值f₀，对矩阵的每N个样本做最大似然估计得到的多个基频值应服从期望为f₀的正态分布。为了排除由各种因素引发的错误估计结果，首先将所得一系列基频值做平均，并求出每个基频相对于平均值的偏差，将偏差值较大的数值取代以平均值；然后在此基础上做线性回归，得到准确的基频估计结果。

附图说明

图1是本发明的流程示意框图。

图2是本发明采用的仿真信号。

图3是本发明最大似然估计算法框图。

图4是本发明得到的频谱灰度图和算法结果。

具体实施方式

下面结合附图，以从一段仿真信号中提取平均骨小梁间距为例，介绍本发明的实施过程。本发明的流程如图1所示，包括四个部分：信号预处理、信号频谱成像、最大似然估计、线性回归。对仿真信号做时域截断，得到松质骨部分的超声背散射信号。如图2所示，采用长度为1200的窗截取背散射信号，图中红色的框内部代表信号中被截取的部分。从矩阵包含样本中取出连续N个样本做最大似然估计，流程如图3所示，其中N_max表示矩阵中最大样本数即矩阵的宽度。

将所得到的基频值取平均，设定阈值，将基频值与平均值相差大于阈值的点去除，并用平均值代替。最后，对基频点做线性回归，得到基频在频谱图中的图像。图4是由仿真信号得到的频谱图，宽表示测试的次数，高表示频率值，每一个红色的圆圈表示该样本区间符合式(17)的基频点。由图4可知，由于强噪声的影响出现了错误估计，然而整体而言依然可以看到正确的基频分布位置。得到基频估计值后，结合已知的超声声速，根据式(18)得到骨小梁的平均间距，

$MTBS=\frac{c}{2*f_0}---\left(18\right)$

式中c为超声声速，f₀为基频值。

Claims (5)

1.一种基于频谱图与最大似然估计的骨密度测量方法，其特征在于具体步骤为：

(1)对采集到的超声回波信号采用窗函数截取松质骨部分并做低通滤波，得到松质骨背散射信号；

(2)多次测量后，将信号做快速傅里叶变换并按照时间顺序组合到信号矩阵中。将该矩阵的成员值映射到[0，255]灰度区间内并显示成灰度图像；

(3)观测信号可表示为多次谐波线性叠加后与噪声之和，谐波频率为基频的整数倍。将基频值和谐波系数看作待估计参数，频谱图像视为观察到的总样本空间。将频谱图像分割为等大小的连续样本区间，由最大似然估计得到信号图像在各区间的基频；

(4)排除奇异点，通过线性回归得到所有基频点的拟合函数。

2.根据权利要求1所述的方法，其特征是在于步骤(1)中，首先对信号做低通滤波，然后用窗函数截取超声背散射信号。

3.根据权利要求1所述的方法，其特征在于步骤(2)中，对于同一部位做多次测量。对每次测量结果按照步骤(1)截取出松质骨背散射信号，并分别作快速傅里叶变换，然后将结果按照时间先后放入信号矩阵中。将信号矩阵中的每一个值映射到[0，255]的灰度区间内，得到灰度图像。

4.根据权利要求1所述的方法，其特征在于步骤(3)中，从信号矩阵内取出连续的N列的样本。骨小梁沿着压力的方向呈近似规则分布，因此在背散射信号中包含着周期性的信号。将观测到的样本表示为多次谐波的线性叠加与噪声之和，如式(1)：

$s\[t\]=\underset{\mathrm{\η}=1}{\overset{N_{\mathrm{\η}}}{\Σ}}\left(a_{\mathrm{\η}}cos\right(2{\mathrm{\π\ηf}}_{n}t)+b_{\mathrm{\η}}sin(2{\mathrm{\π\ηf}}_{n}t\left)\right)+v\[t\],t=n,n+1,\mathrm{....},n+N-1---\left(1\right)$

其中s为离散化后的信号；η表示谐波阶数，阶数越高对信号的细节复原越完整，然而要求的数据量越多；a_η、b_η代表谐波的系数，f_n为基频；。用L表示谐波的阶数，上式的矩阵形式可以表示为：

A(f_n)x_n+v[n]＝s[n] (2)

x_n为谐波系数组成的向量，表示为：

x_n＝[a₁ a₂ a₃ ... a_L b₁ b₂ ... b_L]^T (3)

A(f_n)为N*2L的矩阵，其中第j行为：

${\left[\begin{array}{c}cos\left(2{\mathrm{\πf}}_{n}j\right)\\ cos\left(2\mathrm{\π}2f_{n}j\right)\\ .\\ .\\ .\\ cos\left(2{\mathrm{\πLf}}_{n}j\right)\\ sin\left(2{\mathrm{\πf}}_{n}j\right)\\ \sin \left(2\mathrm{\π}2f_{n}j\right)\\ .\\ .\\ .\\ \sin \left(2{\mathrm{\πLf}}_{n}j\right)\end{array}\right]}^T---\left(4\right)$

将f_n以及x_n看作待估计参数，对矩阵中的N个样本做最大似然估计，公式如下：

$J(f_n,x_n)=min\left({||A\left(f_n\right)x_n-s\[n\]||}_2^2\right)---\left(5\right)$

5.根据权利要求1所述的方法，其特征在于步骤(4)中，首先求得一系列基频值的均值，设定区间排除错误估计点；然后对于基频做线性回归，从而在频谱图上得到基频的拟合函数。该函数与频率轴的交点即基频值，最后由下式得出骨小梁平均间距：

$MTBS=\frac{c}{2*f_0}---\left(6\right)$

其中c为超声在松质骨中传播的速度，f₀为基频，MTBS为平均骨小梁间距。