基于经验模态分解的超声背散射零差K成像方法与流程

本发明属于医学信号处理领域，特别涉及一种经验模态分解的超声背散射零差k成像方法，建立基于经验模态分解即emd的信号分解技术，并发展以零差k模型为基础的有效散射子个数以及相干散射参数成像，用于生物组织的超声组织定征，提供相关微结构信息，以此提供具体散射结构物理意义及参数数值，并提供可视化的参数成像。

背景技术：

超声成像具有实时性好、无电离辐射、价格低廉等优点，已经成为成像检测技术中主流的检测技术手段。在现如今的超声检测技术中，由超声探头发射的超声波如果遇到直径小于入射波长的界面会产生散射，其能量向各个方向辐射。其中背离探头的散射波即背向散射波。由此，我们可以将生物组织(如肝脏、乳腺等)结构视为许多微小散射子所构成的声学结构模型，散射子的相关信息直接反映了生物组织微结构的相关信息，所以生物组织的超声组织定征可以通过测量散射子相关信息所完成。当超声波入射至生物组织中，入射波与散射子交互作用下，探头所接受到的信号，是个别散射子所贡献背散射信号于空间上的叠加。生物组织在发生病变的情况下，其组织本身的微结构将会发生变化，由此会影响到所接收背散射信号的波形特征，基于背散射信号的随即性质，通过分析其概率分布模式，可以归纳出背散射信号统计分布与微小散射子的关联性。因生物体组织散射模式的多样化，后来衍生出许多利用数学上的统计分布，来描述和模式化散射统计，以此来进行生物组织微结构的研究，进行超声组织的定征。

近年来，随着超声检测技术的不断发展，以及相关统计模型应用的不断完善，国内外学者提出了许多基于统计模型的超声组织定征的方法，其主要包括瑞利模型、声学结构定量asq、nakagami模型、复合nakagami、零差k模型等。然而在现有超声组织定征方法中，不论asq技术还是nakagami统计参数成像技术，两者对于散射结构的物理意义都无法提供太多具体信息，由于nakagami参数m在组织定征中极容易出现参数饱和的情况，导致其估算的微观散射子结构不够准确，存在局限性，并且在生物组织的背散射信号中，不光存在目标生物组织的背散射信号，还包含其他噪声信号，这对后续的超声生物组织的定征以及微结构的测量产生很大影响。本发明提出一种基于经验模态分解的超声背散射零差k成像生物组织定征新方法，通过定量有效散射子个数以及相干散射参数，可用于生物组织的超声组织定征，提供微结构信息，并提供可视化的参数成像。由此满足在研究领域对生物组织微结构上的定量研究，具有重要的研究价值和应用前景。

技术实现要素：

本发明旨在提出一种基于经验模态分解的超声背散射零差k成像生物组织定征新方法。本发明由两个部分组成，第一部分是将超声仪器采集到的原始射频信号进行经验模态分解，由此来获取第一个和第二个本征模态函数imf1和imf2，其主要目的是来应对非目标组织所产生的噪声背散射信号对生物体组织微结构定征准确性的影响，提升参数估算的准确性。第二部分是本发明的核心部分，即基于零差k模型的参数估算以及参数可视化成像，首先在参数估算之前会生成b模式超声图像，再通过超声图像来勾画我们要进行参数估算的感兴趣区域，在感兴趣区域内应用滑动窗口技术来进行参数估算，参数估算主要涉及两个估算方法：(1)基于信噪比、偏度、峰度的rsk参数估算方法；(2)另一个是基于x统计和u统计的xu参数估算方法，通过估算得到有效散射子个数μ以及相干散射信号与弥漫散射信号的比值k，其参数可以生成参数图像，以上参数可以进行生物组织结构信息的定征，提供微结构数据，参数成像可以提供可视化信息。

具体技术内容包括以下步骤：

一种基于经验模态分解的超声背散射零差k成像方法，包括以下步骤：

步骤1.b模式超声成像及感兴趣区域获取，具体包括：

步骤1.1.对于原始超声射频数据i，获取其数据的列数为l，其数据的行数为d；

步骤1.2.利用希尔伯特变换解调原始超声射频数据i以构建包络图像，并且使用动态范围为40db进行对数压缩获得包络数据is来构建b模式超声图像b；

步骤1.3.对图像b进行多边形选取，对选取的多边形转换成二值化矩阵的掩模，其中区域内为1，区域外为0，由此便得到了感兴趣区域r；

步骤2.利用经验模态分解对原始射频数据进行信号分解，获取本征模态函数即imf，具体包括：

步骤2.1.进行初始化，令r0＝i,i＝1，其中r0表示还未分解的残余分量，此时即为原始超声射频数据i，其中i为经验模态分解次数的计数；

步骤2.2.令hj-1＝ri-1,j＝1其中hj-1为每一步分解后减去本证模态函数的新序列，此时即为上一步的残余分量ri-1，其中j为判定是否为本征模态函数次数的计数；

步骤2.3.找出hj-1的局部极值点，对其极大值点和极小值点进行三次样条函数插值，形成上下包络线，计算上下包络线的平均值mj-1，hj＝hj-1-mj-1；

步骤2.4.hj若不存在负的局部极大值和正的局部极小值，hj是imf函数，则第i个本证模态函数imfi＝hj，否则，j加一，从步骤2.3重新开始；

步骤2.5.分解后残余分量ri＝ri-1-imfi，如果ri极值点个数仍多于2个，则i加一，从步骤2.2重新开始，否则，分解结束，由此得到本征模态函数imfi。

步骤3.对分解后的数据进行零差k模型参数估算和参数成像，具体包括：

步骤3.1.采用滑动窗口技术作为构造超声参数图像的方法，使用imfi包络信号，利用窗边长度wsl为脉冲长度1～9倍的正方形窗口获取信号数据，在此窗口内进行零差k模型参数估算，窗口以距离增量在图像数据的整个范围内移动，分别以横向和纵向移动，窗口每个增量包含对应于窗口重叠率wor的像素数，窗口移动增量由单个窗口轴向像素数和横向像素数决定，首先确定窗口的相关参数，其中窗口轴向像素数inte＝wsl*p/t，窗口横向像素数inte1＝wsl*p/scan_step，窗口横向移动距离lateral＝l-inte1，窗口横向移动步长l1＝intel1/wor，窗口轴向移动距离axial＝d-inte，窗口轴向移动步长a1＝intel/wor，其中p为脉冲长度、t为轴向点距离、scan_step为扫描步长；

步骤3.2.计算窗口横向移动格数x＝lateral/l1，计算窗口轴向移动格数y＝axial/a1；

步骤3.3.此时信号数据被窗口划分成m*n的矩阵，其中m＝1、2、3…x，n＝1、2、3…y，窗口先在第1行横向移动，在每个窗口内利用rsk法和xu法分别计算得到其对应的参数μ11、μ12、μ13…μ1n，依次类推，直到窗口在第m行完成移动为止，由此获得参数矩阵μmn和kmn，即为有效散射子个数参数矩阵μ和相干散射信号与弥漫散射信号的比值参数矩阵k；

步骤3.4.将参数μ和k利用颜色映射形成参数图像rsk-k、rsk-μ、xu-k和xu-μ,其中rsk-k为利用rsk法所得到的k参数图像，rsk-μ为利用rsk法所得到的μ参数图像，xu-k为利用xu法所得到的k参数图像，xu-μ为利用xu法所得到的μ参数图像；

步骤3.5.利用参数图像rsk-k、rsk-μ、xu-k和xu-μ，计算感兴趣区域r内参数μ和k的平均值，用于生物组织的超声组织定征，提供相关微结构信息，以此提供具体散射结构物理意义及参数数值，并提供可视化的参数成像。

有益效果：

1.本发明方法采用了新的零差k模型参数成像方法，并且使用滑动窗口技术，为参数估算提供了更为便利的方法，为生物组织的超声组织定征开辟了新的研究渠道。

2.本发明采用基于经验模态分解的信号分解技术，可以有效避免非目标生物组织的噪声散射信号对定征的影响。

3.本发明算法可以将评价结果进行可视化展示，不光实现了生物组织相关微结构的量化分析，也提供了可视化分析。

附图说明

图1：本发明方法的流程图；

图2：通过包络信号所生成的b模式超声图像；

图3：本发明方法中零差k超声背散射统计参数成像流程(基于rsk和xu参数估算方法)；

具体实施方式

结合附图和实际例子对提取过程进行具体的描述。所使用数据为使用terason机器来收集的，采集部位是肝右叶，影像深度depth是8cm，探头中心频率fc是3mhz，采样率fs是12mhz，声波速度v是1540，窗口重叠率wor是0.5。本实施例中的生物组织以肝脏组织为例。总共采集了43例肝脏数据。下面分步进行介绍：

1.使用原始超声射频数据i，将原始射频数据以二维数组的形式读取，并获取其数据的列数为l，其数据的行数为d，利用希尔伯特变换解调i以构建包络信号，并且使用动态范围为40db进行对数压缩获得包络数据is来构建b模式图像b，图像如图2所示。根据医学经验，对图像b进行多边形选取(选取肝实质部分)，对选取的多边形转换成矩阵(二值图)的掩模，由此便得到了肝实质的感兴趣区域，同时保存感兴趣区域r为.mat文件以便后续使用。

2.使用经验模态分解算法对原始超声射频数据i进行信号分解，以便获取前两个本征模态函数imf1和imf2。具体实施步骤如下：

(1).进行初始化，令r0＝i,i＝1；

(2).得到第i个imf

(a)h0＝ri-1,j＝1；

(b)找出hj-1的局部极值点；

(c)对其极大值点和极小值点进行三次样条函数插值，形成上下包络线；

(d)计算上下包络线的平均值mj-1；

(e)hj＝hj-1-mj-1；

(f)若hj是imf函数，则imfi＝hj，否则，j＝j+1，转到(b)；

(3).ri＝ri-1-imfi；

(4).如果ri极值点个数仍多于2个，则从步骤(2)重新开始，否则，分解结束，由此得到本征模态函数imfi。算法最后可得

3.在基于零差k统计分布的基础上，利用rsk和xu参数估算方法进行有效散射子成像和相干弥漫比成像。参数成像流程图详见图3。使用未压缩的包络信号imf，采用基于滑动窗口算法进行。具体步骤如下：

(a)首先确定滑动窗口大小，计算窗口轴向像素数inte＝round(wsl*p/y)，计算窗口横向像素数inte1＝round(wsl*p/scan_step)

(b)计算窗口横向和轴向移动距离，横向移动距离lateral＝l-inte1，窗口轴向移动格数axial＝d-inte；

(c)计算窗口横向和轴向移动的步长，窗口横向移动步长l1＝intel1/wor，窗口轴向移动步长a1＝intel/wor，窗口横向移动格数x＝lateral/l1，计算窗口轴向移动格数y＝axial/a1；

(d)此时信号数据被窗口划分成m*n的矩阵，其中m＝1、2、3…x，n＝1、2、3…y，窗口先在第1行横向移动，在每个窗口内利用rsk法和xu法分别计算得到其对应的参数μ11、μ12、μ13…μ1n，依次类推，直到窗口在第m行完成移动为止，由此获得参数矩阵μmn和kmn，即为有效散射子个数参数举证μ和相干散射信号与弥漫散射信号的比值参数矩阵k

其中rsk算法计算流程如下：

rsk法采用包络振幅a的信噪比(snr)、偏度(skewness)和峰度(kurtosis)来估算零差k分布模型参数。信噪比r、偏度s和峰度k与a的任意阶矩的关系如式

式中，v为正实数，e表示数学期望。

零差k分布的任意阶矩可表示为：

e[a^v]＝(2σ²/u)j(k,μ,v)(4)

式中，j(k,u,v)定义如式

式中，pfq(a1,...,ap；b1,...,bq；c)表示广义超几何级数，csc(·)是余割函数，η定义如下：

使用线性插值计算j(k,u,v)：

同理，也可得到sv和kv关于k和μ的函数。通过对比r、s、k的估计值与零差k分布模型预测的r、s、k的理论值，即通过在(k,μ)二维参数空间搜寻等高线(levelcurve)，v取不同值时，等高线的交点即为k和μ的最优估计。

其中xu算法流程如下：

定义参数β为相干散射能量s²与漫散射能量2σ²α的比值β＝s²/(2σ²α)，则xu统计法根据包络振幅的强度i的均值和x统计、u统计的值来估算α和β。u统计和x统计如式

u＝e[logi]-loge[i](9)

x＝e[ilogi]/e[i]-e[logi](10)

则模型参数的估算可转化为以下非线性方程组的求解：

式中，uhk和xhk的定义如下：

式中，γe为欧拉常数，ψ(·)表示双伽马函数。根据μ和β估算k：

(e)最后集合矩阵数据会得到整体的有效散射子μ和相干弥漫比成像k，对其参数进行颜色映射，最终可得到参数图像rsk-k、rsk-μ、xu-k和xu-μ。

(f)利用参数图像rsk-k、rsk-μ、xu-k和xu-μ，计算感兴趣区域r内参数μ和k的平均值，用于肝脏的超声组织定征，因为肝脏组织中散射子的变化反应了肝脏组织微结构的信息，即参数μ和k，并且此参数可以进行可视化参数成像，因此实现基于经验模态分解的超声背散射零差k成像方法。