CN115049568A - 一种基于超声信息熵图像与零差K分布α参数图像融合的定征生物组织的方法 - Google Patents

Abstract

一种基于超声信息熵图像与零差K分布α参数图像融合的定征生物组织的方法属于超声影像技术领域。此发明将超声信息熵图像和零差K分布α参数图像进行信息融合，并生成可以清晰显示组织变性边界的超声伪彩色图像。该融合方法识别生物组织变性区域的面积准确率达到87.68％。

Description

一种基于超声信息熵图像与零差K分布α参数图像融合的定征 生物组织的方法

技术领域

本发明属于超声影像技术领域。

背景技术

超声影像可以被广泛应用于产科、妇科、泌尿科等多个领域，为医护人员检查提供了参考信息。为获得受测组织的超声影像，探头向受测组织发射超声波，并经一定延时后接收到受测组织返回的超声回波。在对超声回波进行波束合成处理后得到若干超声回波数据，并可基于超声回波数据生成对应不同成像时刻的多帧超声影像，通过多帧超声影像的显示，医护人员可获得受测组织的实时影像。然而，在传统的B-mode图像中，生物组织变性部位边界模糊不易识别，使得检测人员不容易准确定位和分辨。Song等人在《Monitoringmicrowave ablation using ultrasound homodyned K imaging based on the noise-assisted correlation algorithm:An ex vivo study》中提出了将超声零差 K(homodyned-K)分布中的α参数成像应用于生物组织的变性监测领域，但由于噪声的存在以及不同超声系统和变性组织个体之间的差异较大，homodyned-K分布α参数图像在不同的个体中可能显示的效果不尽相同。Tsui等人在《Effects of Fatty Infiltration of theLiver on the Shannon Entropy of Ultrasound Backscattered Signals》中将超声信息熵成像首次用于组织定征，取得了显著的效果。

发明内容

1.此发明利用现有超声扫描设备，实现了在超声影像中对生物组织变性区域的识别与监测。

2.此发明将超声信息熵图像和零差K分布α参数图像进行信息融合，并生成可以清晰显示组织变性边界的超声伪彩色图像。

步骤1：如图1所示，由超声探头扫描目标区域后接收到的回波信号进行波束成形处理，即可得到原始的M行N列的超声背向散射信号矩阵。N列：生成的矩阵中有N条扫描线，M行：每条扫描线中包含M个采样点。两个相邻采样点之间的距离为Intaxi米，两个相邻扫描线之间的距离为Intlat米。

步骤2：对步骤1中得到的背向散射信号矩阵进行包络检测处理得到M行N列的超声包络信号矩阵。

步骤3：对步骤2中得到的M行N列包络信号矩阵，采用滑动窗口法来估算其大小为M行N列的homodyned-K分布α参数矩阵。其具体实施方式为：采用一个大小为M_wхN_w的矩形滑动框从包络矩阵的左上角第一个元素开始从左向右、从上到下以一个元素为步长逐行滑动遍历整个包络矩阵。每个矩形窗包含M_wхN_w个数据点，Mw＝<3хpulselength/Intlat>， Nw＝<3хpulselength/Intaxi>，其中pulselength为超声脉冲长度，<>表示向上取整，每滑动一次窗口位置，则对窗口内的包络信号的概率密度函数P建模为homodyned-K分布：

其中，A表示包络信号的振幅值；r是积分变量；J₀(·)是零阶第一类贝塞尔函数；α表示超声分辨率单元内有效散射子的个数；ε²表示相干散射信号能量；σ²表示弥漫散射信号的能量。

在本发明中采用Destrempes等人^[1]提出的XU算法来求解该分布中的α参数。XU法是基于信号强度的一阶矩、X统计量和U统计量的方法。U统计量和X统计量定义如式(2) 和(3)，I表示包络信号振幅的强度，I＝A²，包络信号的平均强度表示为

U＝E[logI]-logE[I] (2)

X＝E[IlogI]/E[I]-E[logI] (3)

XU法估算α参数可以通过求解如下的非线性方程组来获得：

X_HK和U_HK的定义，如式(5)、式(6)所示，其中下标HK表示homodyned-K分布:

式中，α，ε²和σ²的物理含义与(1)式相同，β＝ε²/(2σ²)为计算辅助变量，Q_α+1(·)表示α+1阶第二类修正贝塞尔函数，γ_E为欧拉常数在本发明中取值为0.5772156649，

表示双伽马函数，_pF_q(a₁,...,a_p；b₁,...,b_q；z)表示广义超几何级数。

该方程组有解的充分条件是：

由于(5)和(6)式中的双伽马函数和超几何级数在现实中有有限的数值，因此在实践中，α的最大值是有限的，将式(4)转换成式(8)进行实际的求解，其中：argmin是指使得argmin(U_HK-U)取值最小时的自变量取值；α_max为α的上限，本发明中取α_max＝59.5。

由于X_HK、U_HK是单调的函数，且U_HK是α的函数，X_HK是β的函数，因此可通过迭代的方法求解(8)得到参数α的估计值。

每滑动一次窗口，得到一个α估计值。当滑动窗遍历完整个包络信号矩阵后，得到一个大小为(M-M_w)х(N-N_w)的α参数矩阵，随后对其进行三样条插值，得到一个大小为MxN的α参数矩阵。

步骤4：对步骤3中得到的α参数矩阵进行伪彩色映射，得到homodyned-K分布的α参数图像。

步骤5：对步骤2中获得的MхN大小的包络矩阵，采用滑动窗口法估算其信息熵矩阵。具体实施方法为：采用一个大小为D_wхL_w的矩形滑动框从包络矩阵的左上角第一个元素开始从左向右、从上到下以一个元素为步长逐行滑动遍历整个包络矩阵。每个矩形窗包含D_wхL_w个数据点，Dw＝<pulselength/Intlat>，Lw＝<pulselength/Intaxi>，其中pulselength 为超声脉冲长度，<>表示向上取整。每滑动一次窗口位置，计算该窗口内的熵值。在此发明中，我们将每一次滑动窗口内的包络数据点重排为1列一维信号，使用公式(9)进行计算：

其中，g代表重排后的一维信号的幅度值，w(g)代表该信号的概率密度函数，g_max和g_min分别为g的最大值和最小值，H_c表示熵值的大小。在实践中，使用公式(9)的离散形式来计算获得信息熵。当滑动窗口遍历完整个包络矩阵后，得到一个大小为(M-D_w)х(N-L_w)的熵值矩阵，对此矩阵进行三样条插值，最终获得一个大小为MхN的熵值矩阵。对此矩阵进行伪彩色映射，得到超声信息熵图像。

步骤6：如图3所示，对步骤4得到的超声homodyned-K分布α参数图像和步骤5获得的超声信息熵图像进行图像融合，得到二者的融合图像。

其融合的具体实施方案为：首先，对步骤4得到的超声homodyned-K分布α参数图像和步骤5获得的超声信息熵图像分别进行二维离散小波分解。设f(x,y)表示一幅二维图像的像素值，x和y为像素点横纵坐标，定义二维母小波为ψ(x,y)，将母小波函数ψ(x,y)经伸缩和平移后得到二维小波序列：

式中：σ为伸缩因子，τ_x，τ_y为小波母函数分别在沿x轴和y轴两个维度上的平移，本发明中选择haar小波作为基函数。二维小波变换是在两个维度对图像进行处理，其小波变换为：

其中，

是

的共轭函数，公式为：

二维离散小波变换利用低通滤波器和高通滤波器对图像进行处理，将图像每一层分解为4个子带，分别为：水平和垂直方向低频子带LL、水平方向低频和垂直方向高频子带LH、水平方向高频和垂直方向低频子带HL、垂直和水平方向高频子带HH。图4显示了图像的2层分解示意图。图中的下标数字表示分解层数。

将步骤4中获得的α参数图像和步骤5中获得的超声信息熵图像转换为灰度图像后，使用Matlab软件中的’wavedec2’函数分别对得到的灰度图像进行二维离散小波分解，共计分解2层。如图3所示，每幅图像在分解后都得到高频部分和低频部分。

两幅图像分解后的低频部分反映了图像的主要结构信息，这里采用加权平均的融合规则进行融合处理。设A和B两幅灰度图像在(x，y)位置的像素值分别表示为f_A(x，y)和f_B(x， y)，则在(x，y)位置低频部分融合后的像素值I(x，y)为：

I(x,y)＝0.5*A(x,y)+0.5*B(x,y) (13)

两幅图像分解后的高频部分反映了图像的细节信息，为了使获得的图像具有良好的细节特征，采用区域特性量测的融合规则进行处理。其处理步骤如下：

1)计算每幅图像在相应分解层对应的区域能量，公式如下:

式中:E(x，y)为以(x，y)为中心的局部区域能量；h是一个3×3的模版，代表对应的权系数；s为图像值的逐点平方运算。h和s的计算公式为：

s(x,y)＝f(x,y)² (16)

2)计算两幅图像的区域匹配度M并设置匹配阈值T，公式为：

t(x,y)＝f_A(x,y)*f_B(x,y) (17)

式中：E_A(i，j)和E_B(i，j)为公式(14)计算的两幅图像的区域能量值。

当M≥T时，说明两幅图像的对应区域能量相接近，采用加权平均的融合方法；当M<T 时，说明两幅图像对应区域能量相差较大，则选取区域能量大的小波系数作为融合图像的系数，本发明中根据经验设置阈值T＝0.7。具体流程如图5所示：

最后，按照图4所示，对融合后的高频部分和低频部分进行小波重构，其重构公式如下：

其中，

为公式(12)所示的

的共轭函数，σ和τ的含义同(10)。该重构过程通过调用Matlab中的’waverec2’函数实现。

由此，可以获得了步骤4中得到的homodyned-K分布α参数图像与步骤5中得到的信息熵图像的融合图像。

步骤7：对步骤6中获得的融合图像分别沿每行和每列实行阶数为6的多项式拟合，并将拟合后的图像中显示出-6dB等高线作为识别凝固区的边界。

附图说明

图1超声B-mode成像原理和homodyned-K分布α参数成像原理图

图2超声信息熵成像计算流程示意图

图3基于小波变换的图像融合方法流程图

图4图像2层离散小波分解示意图

图5小波分解后高频部分融合算法流程图

图6实验仪器设备图

图7离体猪肝剖面图(金标准)

图8Bland-Altman图：融合信息熵图像与homodyned-K分布α参数图像

具体实施方式

本发明融合了超声信息熵图像和零差K分布α参数图像，实现了对超声影像中扫描的生物组织变性区域的识别和检测。此方法的实现不依赖于额外的设备，可以将此算法集成到现有超声扫描仪中，能够保证检测的实时性。此方法克服了单独使用B-mode图像监测过程中不能清晰显示生物组织变性区域边界的缺点。此方法克服了单独使用零差K(homodyned-K)分布参数成像监测微波消融过程中受噪声和个体差异影响的缺点。此方法可以将变性的生物组织的边界叠加显示在融合图像上，使得操作人员可以在实时观测并区分变性的生物组织和正常组织。

1.步骤3中对滑动窗口内的包络信号建模为homodyned-K统计分布，此分布模型也可以替换为Nakagami分布模型，通过求解Nakagami分布中的能反映散射子浓度的m 参数来替代homodyned-K分布中的α参数。

2.步骤3中对滑动窗口内的包络信号建模为homodyned-K统计分布，并求解该分布的α参数进行成像。可以将此α参数替换为求解该分布中另一个能够反映散射子结构的k参数进行成像。

3.步骤4中采用XU算法来求解homodyned-K分布中的α参数，可以将此求解算法替换为Hruska等人提出的RSK法^[2]。

实验验证：

为了对此发明方法的准确性进行验证，用微波消融了15例离体猪肝来模拟生物组织变性的情况，实验采用的超声设备为Terson T3000便携式超声扫描仪，对应的超声探头为中心频率为7.5MHz的线阵换能器，换能器脉冲长度0.7mm。使用的微波消融设备为南京康友医疗科技有限公司生产的KY-2000微波消融仪，如图6所示。离体猪肝的盛放于透明的亚克力材质盒，盒子中间有直径3mm空洞由于插入水冷式微波消融针。超声探头扫描平面垂直于消融针方向并覆盖针尖发热点，以此点作为消融中心。

消融开始后，利用步骤1采集得到超声背向散射信号，并以每秒2帧的速率进行保存。对所有15例离体猪肝采用80W功率加热1分。消融结束之后，沿超声扫描平面切开组织并拍照，如图7，以解剖面图像为金标准并将凝固坏死的组织变性区域轮廓视为椭圆形，测量实际的变性区域面积S_gold。

导出消融结束时刻由步骤1得到的那一帧背向散射信号，随后采用步骤2至7来获得超声信息熵图像和homodyned-K分布α参数图像的融合图像，通过融合图像可以检测识别的组织变性区域的面积。为了说明此方法识别凝固区的准确性，将步骤7中获得的-6dB 等高线内面积作为测量面积，与实际组织变性区域面积进行比较，利用公式(20)来计算面积识别准确率:

其中，ACC为组织变性区域面积检测精度；S_-6dB为-6dB等高线内计算的面积；S_gold为消融后离体猪肝剖面测量得到的变性区域金标准面积值。使用信息熵图像融合homodyned-K 分布α参数图像对15例由微波消融引起的的离体猪肝的组织变性区域面积估算与检测精度如表1所示。

表1融合信息熵图像与homodyned-K分布α参数图像组织变性区域检测结果

从表1可以看出，该融合方法识别生物组织变性区域的面积准确率达到87.68％。

15例数据都位于95％一致性界限区间，说明超声信息熵图像和零差K分布α参数图像融合估算的生物组织变性区域与实际变性区域具有一致性，证明本发明的融合方法可以更准确的在超声影像中检测生物组织变性区域。

1.DESTREMPES F,

J,CLOUTIER G.Estimation method of the homodynedK-distribution based on the mean intensity and two log-moments[J].SIAMJournal on Imaging Sciences,2013,6(3):1499-1530.

2.HRUSKA D P,OELZE M L.Improved parameter estimates based on thehomodyned K distribution[J].IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control,2009,56(11):2471-2481。

Claims (1)

1.一种基于超声信息熵图像与零差K分布α参数图像融合的定征生物组织的方法，其特征在于：

步骤1：由超声探头扫描目标区域后接收到的回波信号进行波束成形处理，即可得到原始的M行N列的超声背向散射信号矩阵；N列：生成的矩阵中有N条扫描线，M行：每条扫描线中包含M个采样点；两个相邻采样点之间的距离为Intaxi米，两个相邻扫描线之间的距离为Intlat米；

步骤2：对步骤1中得到的背向散射信号矩阵进行包络检测处理得到M行N列的超声包络信号矩阵；

步骤3：对步骤2中得到的M行N列包络信号矩阵，采用滑动窗口法来估算其大小为M行N列的homodyned-K分布α参数矩阵；

其具体实施方式为：采用一个大小为M_wхN_w的矩形滑动框从包络矩阵的左上角第一个元素开始从左向右、从上到下以一个元素为步长逐行滑动遍历整个包络矩阵；每个矩形窗包含M_wхN_w个数据点，Mw＝<3хpulselength/Intlat>，Nw＝<3хpulselength/Intaxi>，其中pulselength为超声脉冲长度，<>表示向上取整，每滑动一次窗口位置，则对窗口内的包络信号的概率密度函数P建模为homodyned-K分布：

其中，A表示包络信号的振幅值；r是积分变量；J₀(·)是零阶第一类贝塞尔函数；α表示超声分辨率单元内有效散射子的个数；ε²表示相干散射信号能量；σ²表示弥漫散射信号的能量；

XU法是基于信号强度的一阶矩、X统计量和U统计量的方法；U统计量和X统计量定义如式(2)和(3)，I表示包络信号振幅的强度，I＝A²，包络信号的平均强度表示为

U＝E[logI]-logE[I] (2)

X＝E[IlogI]/E[I]-E[logI] (3)

XU法估算α参数通过求解如下的非线性方程组来获得：

X_HK和U_HK的定义，如式(5)、式(6)所示，其中下标HK表示homodyned-K分布:

式中，α，ε²和σ²的物理含义与(1)式相同，β＝ε²/(2σ²)为计算辅助变量，Q_α+1(·)表示α+1阶第二类修正贝塞尔函数，γ_E为欧拉常数在取值为0.5772156649，

表示双伽马函数，_pF_q(a₁,...,a_p；b₁,...,b_q；z)表示广义超几何级数；

该方程组有解的充分条件是：

由于(5)和(6)式中的双伽马函数和超几何级数在现实中有有限的数值，因此在实践中，α的最大值是有限的，将式(4)转换成式(8)进行实际的求解，其中：argmin是指使得argmin(|U_HK-U|)取值最小时的自变量取值；α_max为α的上限，取α_max＝59.5；

由于X_HK、U_HK是单调的函数，且U_HK是α的函数，X_HK是β的函数，因此可通过迭代的方法求解(8)得到参数α的估计值；

每滑动一次窗口，得到一个α估计值；当滑动窗遍历完整个包络信号矩阵后，得到一个大小为(M-M_w)х(N-N_w)的α参数矩阵，随后对其进行三样条插值，得到一个大小为MxN的α参数矩阵；

步骤4：对步骤3中得到的α参数矩阵进行伪彩色映射，得到homodyned-K分布的α参数图像；

步骤5：对步骤2中获得的MхN大小的包络矩阵，采用滑动窗口法估算其信息熵矩阵；具体实施方法为：采用一个大小为D_wхL_w的矩形滑动框从包络矩阵的左上角第一个元素开始从左向右、从上到下以一个元素为步长逐行滑动遍历整个包络矩阵；每个矩形窗包含D_wхL_w个数据点，Dw＝<pulselength/Intlat>，Lw＝<pulselength/Intaxi>，其中pulselength为超声脉冲长度，<>表示向上取整；每滑动一次窗口位置，计算该窗口内的熵值；在此发明中，将每一次滑动窗口内的包络数据点重排为1列一维信号，使用公式(9)进行计算：

其中，g代表重排后的一维信号的幅度值，w(g)代表该信号的概率密度函数，g_max和g_min分别为g的最大值和最小值，H_c表示熵值的大小；在实践中，使用公式(9)的离散形式来计算获得信息熵；当滑动窗口遍历完整个包络矩阵后，得到一个大小为(M-D_w)х(N-L_w)的熵值矩阵，对此矩阵进行三样条插值，最终获得一个大小为MхN的熵值矩阵；对此矩阵进行伪彩色映射，得到超声信息熵图像；

步骤6：对步骤4得到的超声homodyned-K分布α参数图像和步骤5获得的超声信息熵图像进行图像融合，得到二者的融合图像；

其融合的具体实施方案为：首先，对步骤4得到的超声homodyned-K分布α参数图像和步骤5获得的超声信息熵图像分别进行二维离散小波分解；设f(x,y)表示一幅二维图像的像素值，x和y为像素点横纵坐标，定义二维母小波为ψ(x,y)，将母小波函数ψ(x,y)经伸缩和平移后得到二维小波序列：

式中：σ为伸缩因子，τ_x，τ_y为小波母函数分别在沿x轴和y轴两个维度上的平移，选择haar小波作为基函数；二维小波变换是在两个维度对图像进行处理，其小波变换为：

其中，

是

的共轭函数，公式为：

二维离散小波变换利用低通滤波器和高通滤波器对图像进行处理，将图像每一层分解为4个子带，分别为：水平和垂直方向低频子带LL、水平方向低频和垂直方向高频子带LH、水平方向高频和垂直方向低频子带HL、垂直和水平方向高频子带HH；

将步骤4中获得的α参数图像和步骤5中获得的超声信息熵图像转换为灰度图像后，使用Matlab软件中的’wavedec2’函数分别对得到的灰度图像进行二维离散小波分解，共计分解2层；每幅图像在分解后都得到高频部分和低频部分；

两幅图像分解后的低频部分反映了图像的主要结构信息，这里采用加权平均的融合规则进行融合处理；设A和B两幅灰度图像在(x，y)位置的像素值分别表示为f_A(x，y)和f_B(x，y)，则在(x，y)位置低频部分融合后的像素值I(x，y)为：

I(x,y)＝0.5*A(x,y)+0.5*B(x,y) (13)

两幅图像分解后的高频部分反映了图像的细节信息，采用区域特性量测的融合规则进行处理；其处理步骤如下：

1)计算每幅图像在相应分解层对应的区域能量，公式如下:

式中:E(x，y)为以(x，y)为中心的局部区域能量；h是一个3×3的模版，代表对应的权系数；s为图像值的逐点平方运算；h和s的计算公式为：

s(x,y)＝f(x,y)² (16)

2)计算两幅图像的区域匹配度M并设置匹配阈值T，公式为：

t(x,y)＝f_A(x,y)*f_B(x,y) (17)

式中：E_A(i，j)和E_B(i，j)为公式(14)计算的两幅图像的区域能量值；

当M≥T时，说明两幅图像的对应区域能量相接近，采用加权平均的融合方法；当M＜T时，说明两幅图像对应区域能量相差较大，则选取区域能量大的小波系数作为融合图像的系数，设置阈值T＝0.7；对融合后的高频部分和低频部分进行小波重构，其重构公式如下：

其中，

为公式(12)所示的

的共轭函数，σ和τ的含义同(10)；该重构过程通过调用Matlab中的’waverec2’函数实现；

由此，获得了步骤4中得到的homodyned-K分布α参数图像与步骤5中得到的信息熵图像的融合图像；

步骤7：对步骤6中获得的融合图像分别沿每行和每列实行阶数为6的多项式拟合，并将拟合后的图像中显示出-6dB等高线作为识别凝固区的边界。

Images