一种基于B超信号评估呼吸和心动周期的方法与流程

本发明涉及医学图像处理领域，更具体地说，涉及一种基于b超信号评估呼吸和心动周期的方法。

背景技术：

呼吸和心动作为高等动物的重要生命体征，其准确监控对临床诊疗起着重要作用。现有方法中，对呼吸的监测一般基于生物体整体运动特征，对心动的监测一般基于心电信号或局部血管脉动。由于检测原理的不同，两者在临床上都是单独监控的，这就不可避免地增加了诊疗过程的复杂性和成本。

由于其非侵入性和非破坏性的特点，b型超声成像系统已被广泛应用于临床医学诊断。在对目标对象无需侵入性外科手术的情况下，b超可实时地提供目标对象的高分辨率图像。因此，b超已成为临床诊疗不可缺少的重要工具之一。考虑到呼吸和心动都会以引起生物组织、特别是胸、腹腔内的组织位移，b超成像系统有潜力同时实现呼吸和心动的测量与监控，从而减少统一诊疗场景中的设备数量需求，降低诊疗过程中发生误判的可能性。

中国专利申请号：cn201480005221.6“消除医学图像中由生理功能引起的运动影响”通过识别超声图像强反射结构位置和亮度随时间的变化来判断生理运动周期，所提取的信息限于时域信息，与本发明中结合相关系数的计算与线性调频z变换提取频域信息来获得生理运动周期有明显的不同。

中国专利申请号：cn201680053911.8“针对使用超声的四维计算机断层摄影成像的呼吸运动补偿”使用四维计算机断层摄影(4dct)成像采集计算机断层摄影数据，然后将所采数据与超声数据同步，从而确定呼吸周期。本发明只需采集超声数据进行二维互相关，且可同步计算呼吸周期和心动周期，与之有明显不同。

中国专利申请号：cn201680075482.4“用于呼吸监视的超声方法和装置”通过采集隔膜附近的超声信号，利用第一和第二回声信号之间的差异来监测呼吸。这与本发明采集单个b超探头信号，通过频谱特性来计算呼吸周期以及确定心动周期有明显的不同。

中国专利申请号：cn201710616876.2“生成超声图像的方法、超声系统和存储介质”利用心电图(edg)信号来确定生物的心动周期，与本发明中基于b超信号同步获取呼吸周期和心动周期有明显不同。

中国专利申请号：cn201711248270.4“用于在心回波描记术中进行相位确定的心率辅助”使用无线心率传感器采集患者的脉搏信号，根据脉搏信号来采集超声图像，利用脉搏信号确定超声图像的相位。这与本发明直接利用b超图像确定生理运动周期有明显的差别。

中国专利申请号：cn201810443553.2“心脏超声设备及其快速选择心动周期时相对应图像的方法”针对心脏超声图像进行处理，自动选择心动周期，且需要预设心动周期并与实时扫查的图像相关联，与本发明采集b超信号，通过频谱分析来计算生理运动周期有明显不同，且本发明还涉及呼吸周期的确定。

技术实现要素：

1.要解决的技术问题

针对现有技术中存在的现有的常规呼吸和心跳监控相互独立进行，实施方式复杂的问题，本发明提供了一种基于b超信号评估呼吸和心动周期的方法，它可以实现常规呼吸和心跳监控，也可用于基于b超的生理门控应用，还可以用于b超图像处理、超声测温、超声弹性成像等问题中的生理运动补偿，对呼吸和心动周期的判断具有准确性高、易实施等特点。

2.技术方案

本发明的目的通过以下技术方案实现。

本发明的一种基于b超信号评估呼吸和心动周期的方法，包括以下步骤：

步骤一、在生物组织热消融过程中采集若干b超图像信号；

步骤二、根据所得b超图像序列，计算获得一系列互相关系数-时间曲线；

得到参考帧对应的互相关系数-时间曲线c1(t)；互相关系数计算方式为：

其中，*表示共轭运算，设biref(x,y)为参考帧信号，bitar(x,y)为目标帧信号，(i0,j0)为roi区域最左下角的像素点的坐标，对于复值b超图像，该运算产生对应的复共轭图像；对于实值b超图像，该运算不使原图像发生变化；

步骤三、根据所得b超图像序列，评估呼吸和心动周期；具体方法如下：

3.1对步骤二中所得互相关系数-时间曲线c1(t)，利用线性调频z变换计算其频谱f1(f)；线性调频z变换计算如下：

其中：其中m，n为正整数，a0，w0为常数，θ0为抽样点的相角，φ0为两抽样点之间的角度差，x(n)为原始信号；选取频谱中幅值最大的峰值所处频率为呼吸频率的可能值fres,1；

3.2选择当前参考帧之后、时间间隔不小于0.1秒的第k帧为参考帧，按照步骤二和步骤3.1，重新计算互相关系数-时间曲线ck(t)、对应的频谱fk(f)和呼吸频率的可能值fres,k；如此重复，计算不少于x可能的呼吸频率值；

3.3对步骤3.2中得到的不少于x个可能的呼吸频率值，进行高斯分布拟合，取拟合所得高斯函数分布中心为最终确定的呼吸频率值fres；则呼吸周期tres确定为1/fres；

3.4对步骤3.2中所得不少于x个可能的呼吸频率值，计算对应的呼吸周期，通过阈值判断对应的参考帧为合格或不合格，超过阈值，则记对应的参考帧为不合格参考帧；剩余参考帧记为合格参考帧；

3.5如步骤3.4中所得合格参考帧少于x个，则重复3.2-3.4中的步骤，使合格参考帧不少于x个；

3.6对任意合格参考帧m，在其对应的频谱fm(f)中，将频率fres,m后的第一个幅度谱峰所处频率记为fm，将频率fm后的第一个幅度谱峰所处频率记为fr，依据判断条件进行判定：

0.5(1-r)(fr+fres,m)≤fm≤0.5(1+r)(fr+fres,m)

其中：r为判断系数，r不大于0.1，r的初始值设置为[0,0.1]之间的任意值；

3.7如上述判定成立，则将fr记为合格参考帧m对应的可能心动频率之一；随后，将步骤3.6中fr后的第一个幅度谱峰所处频率记为fr，重新按步骤3.6中的条件进行判定；如此重复；

3.8将fm后的第一个幅度谱峰所处频率记为fm后，将fm后的第一个幅度谱峰所处频率记为fr，按照步骤3.6-3.7进行判定；如此重复；对合格参考帧m对应的所有可能心动频率，取对应的fm处幅值最大的情形，记其中的fr为参考帧m对应的心动频率；

3.9对步骤3.5中所得全部合格参考帧，重复步骤3.6-3.8，得到所有参考帧对应的心动频率；对所有心动频率值进行高斯分布拟合，取拟合所得高斯函数分布中心为待选心动频率值，对应的高斯分布尺度参数为σ；

3.10将步骤3.6中的调节判定系数的值，使待选心动频率值分布最集中；选择使3.9中的高斯分布尺度参数最小的判断系数r，将对应的待选心动频率值记录为最终确定的心动频率值fhb；则心动周期thb确定为1/fhb。

更进一步的，步骤一，具体步骤包括：

1.1根据需求确定b超探头的类型、波束控制方式和几何形状；

1.2利用所选探头结合b超仪主机对活体进行监控，采集按时间等间隔排列的b超图像序列。

更进一步的，步骤二具体包括：

2.1对采集的每帧b超图像，将像素点在横向和纵向的坐标分别记为x和y；对所有帧，确定b超图像的矩形roi区域，roi为感兴趣区域，roi区域横向和纵向像素点数分别记为i和j；roi区域为整幅b超图像或b超图像中的局部；选取较大的roi区域有利于获得更为准确的评估结果，选取较小的roi区域有利于加快评估运算的速度；

2.2将采集得到的所有b超图像按时序排列并编号；将第i帧的编号除以ffps获得各帧对应的时间坐标，i＝1～n；将时间变量记为t，单位为秒；将频率变量记为f，单位为赫兹；选取第1帧图像作为参考帧信号；

2.3以所有b超图像帧分别作为目标帧信号，逐一与参考帧计算互相关系数；以目标帧的时间坐标为横坐标，以互相关系数为纵坐标，得到参考帧对应的互相关系数-时间曲线c1(t)；互相关系数计算方式为：

其中，*表示共轭运算，设biref(x,y)为参考帧信号，bitar(x,y)为目标帧信号，(i0,j0)为roi区域最左下角的像素点的坐标，对于复值b超图像，该运算产生对应的复共轭图像；对于实值b超图像，该运算不使原图像发生变化。

更进一步的，步骤1.1选择的b超探头类型为凸阵探头、线阵探头或相控阵探头；波束控制方式为线扫、相控阵、机械扇扫或面阵；波束几何形状为弧形、圆形或矩形。

更进一步的，所述步骤1.2中，b超监控区域包括生物活体的腹部或胸部。

更进一步的，步骤1.2中采集的b超图像序列的采集帧频ffps不低于5帧/秒，采集时长不短于4秒；采集的b超图像序列中每帧图像由探头接收的原始超声回波信号经波束合成后形成或原始超声回波经正交检波、降低采样率后，再经波束合成形成；每帧图像为复值图像或实值图像；所有图像帧的格式相同。

更进一步的，步骤3.3中高斯分布拟合，公式为，a是曲线尖峰的高度，b是尖峰中心的坐标，σ称为标准方差。

更进一步的，步骤3.4中，阈值选择为：可能呼吸频率与确定呼吸频率相比误差超过后者15％所对应的参考帧记为不合格参考帧，对应的互相关系数-时间曲线的频谱不用于计算心动周期。

更进一步的，步骤3.2、3.3、3.4、3.5中“x”范围为5-10。

更进一步的，步骤3.10调节判定系数的值使待选心动频率值分布最集中具体方法为，判断系数r以不超过0.05的间隔步进，每步进一次均重复步骤3.6-3.9。

3.有益效果

相比于现有技术，本发明的优点在于：

(1)本方案可根据b超图像实时评估呼吸和心动周期，因而可以增强临床诊断设备的集成度，有效地减少所需的诊断或监控设备的数量；

(2)基于b超图像序列，可以同步地输出心动周期和呼吸周期；更重要的是，本方法实际上判断了生物组织由于呼吸和心动导致的位移规律，因此可以用于开发对超声图像进行运动补偿相关的图像处理技术；

(3)在基于生理门控技术的各类b超应用中，如超声测温、超声弹性成像，根据本方法获得的呼吸和心动周期可以直接反馈到b超主机，从而控制超声信号的发射和接收，而不必额外增加任何其它设备，成本低。

附图说明

图1为本发明方法步骤流程图；

图2为本发明效果验证方式示意图；

图3为一组实验的互相关系数-时间图，以及相应的频谱幅度曲线；

图4为一组实验中呼吸和心动周期评估效果分布图；

图中标号说明：

1-超声诊断仪；2-成像探头；3-生物活体；4-心率监控仪；5-呼吸机；6-计算机。

具体实施方式

下面结合说明书附图和具体的实施例，对本发明作详细描述。

实施例1

本发明涉及一种基于b超信号评估呼吸和心动周期的方法，具体而言是利用b超信号，通过互相关系数计算、信号处理与心率评估公式相结合，测量生物体呼吸和心动周期的方法；属于医学图像处理、医学超声成像领域。本发明与已有的基于b超图像计算生理运动周期的方法相比，采用了序列图像的互相关计算，及将这种互相关信号转换至频域，或再次进行自相关后转换至频域分析的方法，可同步获得呼吸和心动周期。

步骤一、在生物组织热消融过程中采集若干b超图像信号；

具体步骤包括：

1.1根据需求确定b超探头的类型、波束控制方式和几何形状；选择的b超探头类型为凸阵探头、线阵探头或相控阵探头；波束控制方式为线扫、相控阵、机械扇扫或面阵；波束几何形状为弧形、圆形或矩形。

1.2利用所选探头结合b超仪主机对活体进行监控，b超监控区域包括生物活体的腹部或胸部。采集按时间等间隔排列的b超图像序列。采集的b超图像序列的采集帧频ffps不低于5帧/秒，采集时长不短于4秒；采集的b超图像序列中每帧图像由探头接收的原始超声回波信号经波束合成后形成或原始超声回波经正交检波、降低采样率后，再经波束合成形成；每帧图像为复值图像或实值图像；所有图像帧的格式相同。

步骤二、根据所得b超图像序列，计算获得一系列互相关系数-时间曲线；

2.1对采集的每帧b超图像，将像素点在横向和纵向的坐标分别记为x和y；对所有帧，确定b超图像的矩形roi区域，roi为感兴趣区域，roi区域横向和纵向像素点数分别记为i和j；roi区域为整幅b超图像或b超图像中的局部；选取较大的roi区域有利于获得更为准确的评估结果，选取较小的roi区域有利于加快评估运算的速度；

2.2将采集得到的所有b超图像按时序排列并编号；将第i帧的编号除以ffps获得各帧对应的时间坐标，i＝1～n；将时间变量记为t，单位为秒；将频率变量记为f，单位为赫兹；选取第1帧图像作为参考帧信号；

2.3以所有b超图像帧分别作为目标帧信号，逐一与参考帧计算互相关系数；以目标帧的时间坐标为横坐标，以互相关系数为纵坐标，得到参考帧对应的互相关系数-时间曲线c1(t)；互相关系数计算方式为：

其中，*表示共轭运算，设biref(x,y)为参考帧信号，bitar(x,y)为目标帧信号，(i0,j0)为roi区域最左下角的像素点的坐标，对于复值b超图像，该运算产生对应的复共轭图像；对于实值b超图像，该运算不使原图像发生变化。

步骤三、根据所得b超图像序列，评估呼吸和心动周期；具体方法如下：

3.1对步骤2.3中所得互相关系数-时间曲线c1(t)，利用线性调频z变换计算其频谱f1(f)；线性调频z变换计算如下：

其中：其中m，n为正整数，a0，w0为常数，θ0为抽样点的相角，φ0为两抽样点之间的角度差，x(n)为原始信号；选取频谱中幅值最大的峰值所处频率为呼吸频率的可能值fres,1；可能存在的方法还包括快速傅里叶变换：其中x(n)为原始信号；n为偶数。

3.2选择当前参考帧之后、时间间隔不小于0.1秒的第k帧为参考帧，按照步骤2.3和步骤3.1，重新计算互相关系数-时间曲线ck(t)、对应的频谱fk(f)和呼吸频率的可能值fres,k；如此重复，计算不少于x可能的呼吸频率值；

3.3对步骤3.2中得到的不少于x个可能的呼吸频率值，进行高斯分布拟合，取拟合所得高斯函数分布中心为最终确定的呼吸频率值fres；则呼吸周期tres确定为1/fres；步骤3.3中高斯分布拟合，公式为，a是曲线尖峰的高度，b是尖峰中心的坐标，σ称为标准方差。分布拟合函数也不仅限于高斯函数，

如还可以为，三角窗函数:

a是曲线尖峰的高度，b是尖峰中心的坐标，n为正整数。

汉宁窗函数:

a是曲线尖峰的高度，b是尖峰中心的坐标，；n为正整数。

海明窗函数:

a是曲线尖峰的高度，b是尖峰中心的坐标，n为正整数。

凯泽窗函数:

i0是第一类零阶贝塞尔函数,β是调整窗函数性能的参数,a是曲线尖峰的高度，b是尖峰中心的坐标,n为正整数等。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的方式及实施例，均应属于本发明的保护范围。

3.4对步骤3.2中所得不少于x个可能的呼吸频率值，计算对应的呼吸周期，通过阈值判断对应的参考帧为合格或不合格，超过阈值，则记对应的参考帧为不合格参考帧；剩余参考帧记为合格参考帧；阈值选择为：可能呼吸频率与确定呼吸频率相比误差超过后者15％所对应的参考帧记为不合格参考帧，对应的互相关系数-时间曲线的频谱不用于计算心动周期。

3.5如步骤3.4中所得合格参考帧少于x个，则重复3.2-3.4中的步骤，使合格参考帧不少于x个；

3.6对任意合格参考帧m，在其对应的频谱fm(f)中，将频率fres,m后的第一个幅度谱峰所处频率记为fm，将频率fm后的第一个幅度谱峰所处频率记为fr，依据判断条件进行判定：

0.5(1-r)(fr+fres,m)≤fm≤0.5(1+r)(fr+fres,m)

其中：r为判断系数，r不大于0.1，r的初始值设置为[0,0.1]之间的任意值；

3.7如上述判定成立，则将fr记为合格参考帧m对应的可能心动频率之一；随后，将步骤3.6中fr后的第一个幅度谱峰所处频率记为fr，重新按步骤3.6中的条件进行判定；如此重复；

3.8将fm后的第一个幅度谱峰所处频率记为fm后，将fm后的第一个幅度谱峰所处频率记为fr，按照步骤3.6-3.7进行判定；如此重复；对合格参考帧m对应的所有可能心动频率，取对应的fm处幅值最大的情形，记其中的fr为参考帧m对应的心动频率；

3.9对步骤3.5中所得全部合格参考帧，重复步骤3.6-3.8，得到所有参考帧对应的心动频率；对所有心动频率值进行高斯分布拟合，取拟合所得高斯函数分布中心为待选心动频率值，对应的高斯分布尺度参数为σ；

3.10将步骤3.6中的调节判定系数的值，使待选心动频率值分布最集中，判断系数r以不超过0.05的间隔步进，每步进一次均重复步骤3.6-3.9；选择使3.9中的高斯分布尺度参数最小的判断系数r，将对应的待选心动频率值记录为最终确定的心动频率值fhb；则心动周期thb确定为1/fhb。

上述步骤3.2、3.3、3.4、3.5中“x”范围为5-10。

具体的实施例如下，

实施例1

如图1、图2、图3和图4所示，一种评估活体生理运动周期，即呼吸周期和心动周期的超声方法，该方法用b型超声成像仪进行成像并收集其原始回波信号。基于b超时序排列所有超声图像，并利用互相关算法计算互相关系数-时间曲线，然后用线性调频z变换算法计算互相关系数的频谱，选取最大峰的对应的位置作为呼吸频率。以呼吸频率为基准，依据判断条件识别心率，进而可以得到活体的呼吸周期和心动周期。

以生物活体为例，该方法具体实施步骤如下：

利用超声诊断仪1的成像探头2对活体猪3进行b型扫描成像，按时间顺序连续输出未经处理的原始超声信号。其中超声诊断仪1，心率监控仪4由计算机6控制及收集数据，呼吸机5单独控制。

步骤一、对猪腹部进行b型超声扫描成像，按时序连续输出原始回波信号。

步骤二、选取初始时刻的一帧b超图像作为参考帧。

将信号中其他帧分别作为目标帧，逐一与参考帧进行相关系数的计算，得到一系列随时间变化的互相关系数，即时间-互相关系数曲线相关系数计算公式如下：如图3(a)所示，

其中biref为参考帧图像信号，bitar为目标帧图像信号；i和j分别为roi区域的横向和纵向的采样点数，(i0,j0)是roi区域的左下角点。这里roi区域即为整幅图像。

3.1、上述步骤中得到的互相关系数-时间曲线，利用线性调频z变换变换计算其频谱；线性调频z变换计算公式：

其中：其中m，n为正整数，a0，w0为常数，q0为抽样点的相角，φ0为两抽样点之间的角度差，其频谱如图3(b)所示，并选取幅度谱中最大峰值，即图3(b)中第一个峰，对应的频率为呼吸频率的可能值。

3.2、以0.5s为时间间隔，多次改变初始帧的位置，按照步骤二-步骤3.1重新计算呼吸频率，得到如图4(b)所示的呼吸频率统计图。

3.3、对所得不少于5个可能呼吸频率值，进行高斯分布拟合，a是曲线尖峰的高度，b是尖峰中心的坐标，σ称为标准方差，对图4(b)中的呼吸数据，得到分布中心值，即得呼吸频率。

3.4、对可能呼吸频率值偏离确定呼吸频率值超过15％的情形，去除对应参考帧的互相关曲线频谱，不用于计算心动频率。

3.5、步骤3.4中所得合格参考帧足够，不需要重复步骤；

3.6、在互相关系数-时间曲线的频谱中，以呼吸频率fres，其后第一个峰值频率fm、fm后的第一个峰值频率fr，依据判断条件进行判别，若判断条件成立，则记录fr为可能的心动频率。判断条件：

0.5(1-r)(fr+fres)≤fm≤0.5(1+r)(fr+fres)

初始条件下，r可设置为0.03。

3.7、在互相关系数-时间曲线的幅度谱中，将fr向右平移至下一个峰值，重复步骤3.6；如此循环；每次判定成立，则记录fr为可能的心动频率。

将fm向右平移至其后一个峰值，重复步骤3.7，如此循环。

3.8、在计算得到的所有可能心动频率中，选取fm处幅值最大的情形，将对应的fr记为当前参考帧对应的心动频率。

3.9、依次选取其它合格参考帧，重复步骤3.6-3.8，获得每个合格参考帧对应的心动频率。对所有心动频率值进行高斯分布拟合，获得分布宽度因子σ，将分布中心记为心动频率待选值。

3.10、将判断系数r在[0,0.1]之间以0.01的步长进行步进，重复步骤3.6-3.9，选取使最小的r值作为判断系数，对应的心动频率待选值为最终确定的心动频率，计算的心动频率值如图4(a)所示。

根据得到的呼吸频率和心动频率分别计算呼吸周期和心动周期。

实施案例：基于上述方案和实施方法，本发明以活体猪为例进行应用。

基于一套自研的超声诊断仪进行超声数据采集，生物活体为猪。超声诊断探头采用128阵元、中心频率为3.5兆赫兹的微凸探头，超声主机和心率仪由计算机控制，计算及同时负责数据收集和存储。深度睡眠条件下，呼吸机设定的猪的呼吸周期为5秒，经心率仪评估，深度睡眠状态下猪的心动周期为1秒左右。

利用b型超声诊断仪进行成像，每次数据采集时长12秒，帧频为38帧/秒，采样频率为40兆赫兹，每幅超声图像的尺寸为128像素×2048像素。在步骤(4)中，采用roi区域为整个超声图像(128像素×2048像素)；对b超图像信号进行正交变换，得到复值图像信号。经步骤(13)计算，判断系数r确定为0.04。

经应用本发明提出的评估方法，可得到互相关系数与时间图，相应的幅度谱图，以及呼吸和心动周期分布图。通过高斯分布拟合，确定呼吸周期为4.93秒，心动周期为0.96秒。

以上示意性地对本发明创造及其实施方式进行了描述，该描述没有限制性，在不背离本发明的精神或者基本特征的情况下，能够以其他的具体形式实现本发明。附图中所示的也只是本发明创造的实施方式之一，实际的结构并不局限于此，权利要求中的任何附图标记不应限制所涉及的权利要求。所以，如果本领域的普通技术人员受其启示，在不脱离本创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均应属于本专利的保护范围。此外，“包括”一词不排除其他元件或步骤，在元件前的“一个”一词不排除包括“多个”该元件。产品权利要求中陈述的多个元件也可以由一个元件通过软件或者硬件来实现。第一，第二等词语用来表示名称，而并不表示任何特定的顺序。