CN112932539B - 一种经颅超声估计血管血流速度场的装置及方法 - Google Patents
一种经颅超声估计血管血流速度场的装置及方法 Download PDFInfo
- Publication number
- CN112932539B CN112932539B CN202110092110.5A CN202110092110A CN112932539B CN 112932539 B CN112932539 B CN 112932539B CN 202110092110 A CN202110092110 A CN 202110092110A CN 112932539 B CN112932539 B CN 112932539B
- Authority
- CN
- China
- Prior art keywords
- blood flow
- ultrasonic
- blood vessel
- blood
- flow velocity
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000017531 blood circulation Effects 0.000 title claims abstract description 118
- 238000000034 method Methods 0.000 title claims abstract description 44
- 238000002604 ultrasonography Methods 0.000 title claims abstract description 28
- 230000002792 vascular Effects 0.000 title claims abstract description 17
- 210000004204 blood vessel Anatomy 0.000 claims abstract description 150
- 238000009826 distribution Methods 0.000 claims abstract description 41
- 238000001914 filtration Methods 0.000 claims abstract description 37
- 238000003384 imaging method Methods 0.000 claims abstract description 37
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 36
- 238000010586 diagram Methods 0.000 claims abstract description 25
- 230000002490 cerebral effect Effects 0.000 claims description 23
- 238000004458 analytical method Methods 0.000 claims description 16
- 238000004364 calculation method Methods 0.000 claims description 16
- 238000000354 decomposition reaction Methods 0.000 claims description 11
- 210000004556 brain Anatomy 0.000 claims description 8
- 239000008280 blood Substances 0.000 claims description 7
- 230000015572 biosynthetic process Effects 0.000 claims description 6
- 210000004369 blood Anatomy 0.000 claims description 6
- 238000003786 synthesis reaction Methods 0.000 claims description 6
- 239000002245 particle Substances 0.000 claims description 3
- 238000002608 intravascular ultrasound Methods 0.000 claims description 2
- 238000005457 optimization Methods 0.000 claims 1
- 230000003313 weakening effect Effects 0.000 abstract 1
- 238000000520 microinjection Methods 0.000 description 13
- 238000012545 processing Methods 0.000 description 9
- 210000003625 skull Anatomy 0.000 description 8
- HRPVXLWXLXDGHG-UHFFFAOYSA-N Acrylamide Chemical compound NC(=O)C=C HRPVXLWXLXDGHG-UHFFFAOYSA-N 0.000 description 7
- 238000002474 experimental method Methods 0.000 description 7
- 239000012530 fluid Substances 0.000 description 7
- 238000006073 displacement reaction Methods 0.000 description 6
- 230000006870 function Effects 0.000 description 6
- 239000007788 liquid Substances 0.000 description 6
- 238000012360 testing method Methods 0.000 description 6
- 210000001519 tissue Anatomy 0.000 description 6
- 239000011159 matrix material Substances 0.000 description 5
- 238000005259 measurement Methods 0.000 description 4
- 238000010146 3D printing Methods 0.000 description 3
- 238000005314 correlation function Methods 0.000 description 3
- 125000004122 cyclic group Chemical group 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 239000002504 physiological saline solution Substances 0.000 description 3
- 239000004626 polylactic acid Substances 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 238000012285 ultrasound imaging Methods 0.000 description 3
- 238000000827 velocimetry Methods 0.000 description 3
- 238000005452 bending Methods 0.000 description 2
- 230000003727 cerebral blood flow Effects 0.000 description 2
- 238000005520 cutting process Methods 0.000 description 2
- 201000010099 disease Diseases 0.000 description 2
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 2
- 230000014509 gene expression Effects 0.000 description 2
- 230000033001 locomotion Effects 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 238000000691 measurement method Methods 0.000 description 2
- 239000000693 micelle Substances 0.000 description 2
- 238000011084 recovery Methods 0.000 description 2
- 238000003756 stirring Methods 0.000 description 2
- 230000001502 supplementing effect Effects 0.000 description 2
- 230000002123 temporal effect Effects 0.000 description 2
- 239000004698 Polyethylene Substances 0.000 description 1
- 230000003187 abdominal effect Effects 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 210000001367 artery Anatomy 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 210000001715 carotid artery Anatomy 0.000 description 1
- 238000002585 cerebral angiography Methods 0.000 description 1
- 230000019771 cognition Effects 0.000 description 1
- 238000013329 compounding Methods 0.000 description 1
- 238000007428 craniotomy Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000002592 echocardiography Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 238000007917 intracranial administration Methods 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 231100000915 pathological change Toxicity 0.000 description 1
- 230000036285 pathological change Effects 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 239000004033 plastic Substances 0.000 description 1
- 229920003023 plastic Polymers 0.000 description 1
- -1 polyethylene Polymers 0.000 description 1
- 229920000573 polyethylene Polymers 0.000 description 1
- 230000002035 prolonged effect Effects 0.000 description 1
- 239000000700 radioactive tracer Substances 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000013341 scale-up Methods 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
- 238000010998 test method Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
- 230000017105 transposition Effects 0.000 description 1
- 238000011282 treatment Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/06—Measuring blood flow
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/5207—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of raw data to produce diagnostic data, e.g. for generating an image
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/5215—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of medical diagnostic data
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Heart & Thoracic Surgery (AREA)
- Molecular Biology (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Physics & Mathematics (AREA)
- Biomedical Technology (AREA)
- Veterinary Medicine (AREA)
- Medical Informatics (AREA)
- Biophysics (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Hematology (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
Abstract
本发明公开了一种经颅超声估计血管血流速度场的装置及方法。装置主要包括低频超声换能器、可编程超声平台和计算机，方法重点在于使用低频超声换能器发射附有线性调频信号的平面波进行成像；对射频数据进行时空滤波，以及造影微泡的定位、追踪得到血管结构图，并确认血管直径及位置坐标；以及将传统图像互相关算法与快速傅里叶变换算法结合应用于速度分布估计中。本发明能够有效地克服经颅超声衰减造成的回波信号强度减弱，并在低频超声换能器成像分辨率低的情况下，较精准地确定血管位置区域，特别适用于经颅超声小血管血流速度分布的估计。
Description
技术领域
本发明属于声学测量技术领域，涉及一种利用超声成像估计血管血流速度场的装置以及对应的血流速度场估计方法。
背景技术
脑小血管病对人的身体、认知和精神方面都有较大的损害，因此，在对脑小血管病的基础研究及诊疗中有必要了解和掌握反映脑小血管生理、病理变化的参量及数据。血管血流速度场是指血管内径向切面上血流速度的空间分布，流体力学中将血管腔近似于圆管，血液在血管内的流动符合圆管内层流流动的情况，血管血流速度场呈现抛物线形状，通过研究血管血流速度场可以进一步得到血管内的各项血流动力学参数。
基于超声散斑的图像测速方法可以用于血管血流速度场的估计，其测量对象均为浅表组织内的血管(以颈动脉、腹动脉等大血管和小血管为主)，采用高频超声换能器发射聚焦超声波并采集造影微泡在血液中流动的回波信号，得到分辨率较高的血管区域内的超声图像，进而利用连续两帧血管图像进行血流速度场的估计。然而经颅超声成像受到颅骨遮挡，虽然通过在血管内释放造影微泡可以增强信号，但受半波长限制和穿颅后声波衰减的影响，无法获得较高分辨率的超声图像，不能有效分辨血管区域，因此，不适合应用于估计经颅亚毫米级别血管(管径0.1～1mm的血管)以及微小血管(管径小于0.1mm的血管)等脑小血管的血流速度场。另外，传统超声成像扫描方式是线扫描，即超声换能器按照阵元排序逐个发射信号，并逐个接受回波信号，这种扫描方式会导致阵元收到的信号并不是同一时刻的，存在帧率不足的问题，影响血流信号的获得。
目前利用超声识别脑小血管并完成其血流速度场的估计较困难的原因不仅在于超声穿透颅骨时会发生明显衰减，还在于在得到的超声图像中无法获取精确有效的血管位置，导致无法在血管内运用超声散斑的图像测速方法估计血流速度分布。中国专利CN107753062A“基于马尔科夫链蒙特卡洛多目标追踪的经颅超声脑血管造影超分辨率成像方法”、中国专利CN107714091A“经颅低频超声线性调频脉冲逆转微泡成像方法”中提到的成像方法有效地得到了脑血管结构图，但脑血管结构图只是由微泡追踪点连接形成的线条图，不包含物理信息，并不能作为超声散斑图像测速方法估计血流速度分布的对象。
发明内容
针对利用超声成像无创估计脑小血管血流速度场的问题，本发明提供一种经颅超声估计血管血流速度场的装置及方法。
为达到上述目的，本发明采用了以下技术方案：
一种经颅超声估计血管血流速度场的方法，该估计血管血流速度场的方法包括以下步骤：
1)利用超声扫描颅脑内部区域，并采集超声回波信号；
2)通过对超声回波信号进行时空滤波、成像及图像定位和追踪，得到超声成像区域内较精确的血管结构图；
3)选取血管结构图中的血管，并推断该血管的位置坐标信息；
4)根据血管的位置坐标信息从时空滤波后通过成像形成的超声序列图像(该时空滤波后的图像是只有血流回波信号的图像)中提取血管内血流的图像，并对提取的相邻帧血管内血流的图像采用划分多级判读窗口进行迭代的快速傅里叶变换相关算法进行分析，得到血管血流速度梯度分布。
优选的，所述步骤1)具体包括以下步骤：采用可编程超声平台控制超声换能器向颅内发射低频线性调频平面波，并由可编程超声平台采集所述超声换能器接收到的回波信号，得到原始射频数据，对原始射频数据进行处理(例如，参照B模式图像的生成)，生成原始超声序列图像，对该超声序列图像未做任何后续处理。
优选的，所述步骤2)具体包括以下步骤：对原始射频数据进行波束合成及奇异值分解(时空滤波，即时间空间滤波采用的主要方法是奇异值分解)，对经过波束合成及奇异值分解的原始射频数据(只保留血流回波信号)进行处理，生成时空滤波后的超声序列图像(图像的血管区域并不精确，是得到血管位置坐标信息后要映射回去的对象，经过映射后才能得到较为准确的血管区域)，对该超声序列图像采用高斯峰值拟合定位得到对应图像中造影微泡的位置，并采用马尔科夫链蒙特卡洛算法对造影微泡的流动进行追踪。
优选的，所述步骤3)具体包括以下步骤：采用测量血管结构图曲线半高宽的方法确定血管直径及位置坐标。
优选的，所述血管选自直径≤1mm的脑小血管(例如，直径0.1～1mm的亚毫米级的血管等脑小血管)。
优选的，所述步骤4)中，对步骤2)中利用经过奇异值分解保留的血流回波信号生成的超声序列图像(即时空滤波后的超声序列图像)进行血管内超声粒子图像信息的截取、扩大及插值，得到血管内血流的图像信息。
优选的，所述估计血管血流速度场的方法还包括以下步骤：采用传统(例如，参考超声散斑的图像测速方法)互相关算法对步骤2)中利用经过奇异值分解保留的血流回波信号生成的超声序列图像(即时空滤波后的超声序列图像)进行速度分布分析，得到有效峰值速度，将有效峰值速度与步骤4)得到的血流速度梯度分布复合，从而实现对血管血流速度场的精准估计。
一种经颅超声估计血管血流速度场的装置，该估计血管血流速度场的装置包括超声换能器、可编程超声平台和计算机；
所述超声换能器用于在可编程超声平台控制下通过发射超声波(例如，上述低频线性调频平面波)扫描颅脑内部区域并接收超声回波信号(回波信号进一步由可编程超声平台采集，并输入计算机)；
所述计算机包括经颅超声估计血管血流速度场的软件系统(用于执行上述经颅超声估计血管血流速度场的方法)，该软件系统包括血管标定模块及血流速度梯度分布计算模块；
所述血管标定模块用于对超声回波信号进行时空滤波、成像及图像定位和追踪，以及推断血管的位置坐标信息；
所述血流速度梯度分布计算模块用于根据血管的位置坐标信息从时空滤波后通过成像形成的超声序列图像(即时空滤波后的超声序列图像)中提取血管内血流的图像，并对提取的相邻帧血管内血流的图像采用划分多级判读窗口进行迭代的快速傅里叶变换相关算法进行血流速度梯度分布分析。
优选的，所述软件系统还包括血流速度梯度分布优化模块，该模块用于采用互相关算法对时空滤波后通过成像形成的超声序列图像(即时空滤波后的超声序列图像)进行速度分布分析，以及用于将采用互相关算法分析得到的有效峰值速度与血流速度梯度分布计算模块计算得到的血流速度梯度分布复合，得到更为精确的血管血流速度场。
优选的，所述超声换能器选自低频相控阵或线阵超声换能器。
优选的，所述超声换能器采用阵元齐发齐收的平面波扫描方式，有助于提高帧率，帧率的提高对流动微泡的回波信号的捕捉能力也提高，更有助于细化微泡流动过程中的瞬时变化。
优选的，所述估计血管血流速度场的装置还包括颅脑小血管仿血液流动系统，用于模拟脑小血管内血液的流动，具体是将无壁血管仿体固定在颅骨模型中，以及将无壁血管仿体分别与微量注射泵及液体回收仓连通，并完成血液流动的模拟。
优选的，所述无壁血管仿体由中空腔体的具有一定弯曲弧度的壁面形成。
本发明的有益效果体现在：
本发明可以在血管位置不明的情况下，通过依据时空滤波及图像定位和追踪对血管位置坐标的标定并映射到时空滤波后的图像中，从而依据图像相应血管区域的血流图像完成对血管血流速度场的估计。相较于直接针对B超图像进行处理的方法而言，能够实现颅内亚毫米等级别脑小血管的血流速度场的估计，且估计结果精准度较高。
进一步的，本发明通过将采用计算互相关函数的算法分析得到有效峰值速度与划分多级判读窗口进行迭代的快速傅里叶变换相关算法分析得到的血流速度梯度分布复合，可以解决血流速度梯度分布中峰值流速与理论值误差大的问题。
进一步的，本发明中超声换能器发射低频线性调频平面波，从而有效地解决经颅超声衰减后导致成像回波信号强度降低的情况，能够有效地保留血流回波信号中的微小微泡信号。同时可以利用平面波提升数据采集的帧率(在追踪以及相邻帧间的互相关计算时，都要求得到的同一帧图像信号是同一时刻的)。
附图说明
图1是经颅超声小血管血流速度场测试装置的示意图。
图2是经颅超声小血管血流速度场测试方法的流程图。
图3是超声低频线性调频平面波的波形示意图。
图4是测量得到的仿体中血管血流速度向量图。
图5是仿体中血管血流速度场测试实验计算值与模拟值之间的吻合度示意图，其中，抛物线为模拟值，离散数据为通过5次实验得到的血流速度梯度计算值的均值与标准差。
图6是利用循环移位法及血管结构图曲线半高宽法估计血管直径及位置坐标的示意图；其中：(A)循环移位前后血管中心点分布及血管位置，(B)结构图曲线半高宽估计直径。
具体实施方式
下面结合附图和实施例对本发明作进一步详细说明。所述实施例仅用于解释本发明，而非对本发明保护范围的限制。
参见图1，本发明提出的经颅超声小血管血流速度场测试装置主要包括微量注射泵(微量注射泵由微量注射控制器8和微量注射执行器7组成，微量注射执行器7由受控于微量注射控制器8的执行单元和安装在该执行单元上的注射器组成)、低频超声换能器3、换能器三维升降支架2和可编程超声平台1。可编程超声平台1和计算机构成一体，是采集和处理超声回波信号的基本设备(可通过编程实现相应地超声功能)，可编程超声平台1与低频超声换能器3连接在一起。在实验中，低频超声换能器3通过耦合剂与3D打印聚乳酸(PLA)颅骨4紧密贴合，换能器三维升降支架2用于固定和调节低频超声换能器3的位置。所述微量注射控制器8、微量注射执行器7、3D打印聚乳酸(PLA)颅骨4、置于3D打印聚乳酸(PLA)颅骨4内的丙烯酰胺脑小血管仿体5(进口端与微量注射执行器7相连通)，以及与该丙烯酰胺脑小血管仿体5出口端相连通的液体回收仓6共同组成了颅脑小血管仿血液流动系统。通过设置微量注射泵的参数来调节仿血液的液体(例如，生理盐水)在丙烯酰胺脑小血管仿体5中的流动速度，其调节范围符合实际的血液流动情况，从而可以模拟在血液中释放造影微泡作为示踪剂的过程。
上述丙烯酰胺脑小血管仿体5(简称为仿体血管)采用内部呈中空状且具有一定弯曲弧度的丙烯酰胺块体，使其中空部分所形成的腔体(直径0.2mm～0.5mm)更加接近于一段具有进出口的真实血管，确保仿体血管内的超声回波信号被更好的接收，避免采用具有较硬管壁(即声阻抗与仿体组织差异较大的管材料)的仿体血管(例如，聚乙烯塑料管)时存在的强管壁反射对血流回波信号干扰(声阻抗差越大，声压反射系数越大，那么管壁的强反射信号对管内血流回波信号产生干扰)。
上述液体在通过微量注射泵注入颅脑小血管仿血液流动系统前，采用磁力搅拌器搅拌其中混有的造影微泡，搅拌的过程使得液体中的造影微泡分布均一，注射进仿体血管后造影微泡在液体中的跟随性较好、流动均匀，成像后的图像可以用于速度分布的估计。
本发明还提出了一种估计脑小血管血流速度场的方法，方法的重点在于使用低频超声换能器3发射附有线性调频信号的平面波进行成像以及将图像互相关算法与快速傅里叶变换相关算法结合应用于脑小血管血流速度场精准估计中。
本发明采用上述经颅超声小血管血流速度场测试装置，进行了仿体血管血流速度场的测试实验，通过实验对上述估计脑小血管血流速度场的方法进行了验证。实验中，首先通过在丙烯酰胺脑小血管仿体5中注入含有一定量造影微泡的生理盐水模拟血液流动和示踪粒子的释放，并采用可编程超声平台1采集相应超声射频数据，进一步对其运用时空滤波去除颅脑仿体(颅脑仿体是将颅骨4和血管仿体5配置在一起而形成的，血管仿体5的具体位置为位于颞窗内侧)组织回波等干扰信号，最大程度地保留血流回波信号，然后选取连续帧图像(时空滤波后的图像，只保留了血流回波信号)按照以下两种方式完成血流速度分布的分析计算：(1)对整体图像划分ROI区域(考虑并去除入口效应)后利用互相关算法完成血流速度分布分析，得到成像平面内的有效峰值速度；(2)对时空滤波后的图像进行造影微泡的定位、追踪，得到仿体血管的结构图，进而完成对血管中血流图像的精确提取，通过对该图像采用划分多级判读窗口进行迭代的快速傅里叶变换相关算法进行分析，从而较精准地计算得到血管内的血流速度梯度分布，结合以上两种方式的分析计算结果实现对血管血流速度场的更精准估计。
参见图2，测试实验的具体步骤如下：
步骤1将颅脑仿体固定，调整换能器三维升降支架2，使得低频超声换能器3经过颞窗向颅内有效发射超声波和接收超声回波信号，并且使得仿体血管呈现的成像平面最大可能地处于中心平面(中心平面即为仿体血管的中轴线部分)，通过微量注射泵向丙烯酰胺脑小血管仿体5内推注含有造影微泡的生理盐水，并达到稳定流动状态。
步骤2在步骤1的基础上，运用可编程超声平台1，执行已编译的程序控制低频超声换能器3发射低频线性调频平面波(低频线性调频平面波低中心频率的特性在穿透颅骨时能够有效减少声衰减、线性调频信号能够提高回波信号的可检测性以及获得超快平面波成像的高帧率特点)，在设定的帧频和帧数范围内(一般情况帧频为100Hz～1000Hz，帧数为1000～2000)，可编程超声平台1采集低频超声换能器3接收的序列超声射频数据并进行保存。
结合线性调频信号的发射波形(其频带范围与低频超声换能器一致，图3)表示为：
其中，0≤t≤T，w(t)为线性调频信号x(t)的窗函数，f0为线性调频信号x(t)的中心频率，B为线性调频的带宽，T为脉冲持续时间。
线性调频信号的脉冲持续时间越长，发射的能量越大，作为示踪剂的造影微泡散射回波的能量越大，成像的信噪比越高。但脉冲持续时间过长会造成浅表成像盲区以及超声换能器温度过高，因此，在实际的超声射频数据采集的过程中，在考虑各方制约参数的情况下，尽可能大的延长线性调频信号的脉冲持续时间，有利于提升成像的分辨率。
所述低频超声换能器3的类型可以是线阵或相控阵，其发射方式为平面波，即换能器所有阵元同时发射及接收信号。低频超声换能器3的频率可以是1MHz～4MHz，符合经颅超声的成像频率范围。
步骤3对采集的超声射频数据进行波束合成后运用奇异值分解保留所采集数据中的血流回波信号的特征，并有效地去除杂波和颅脑仿体组织回波信号，形成连续多帧图像，即时空滤波后的图像，然后运用高斯峰值拟合定位获得各帧时空滤波后的图像中造影微泡的位置，并使用马尔科夫链蒙特卡洛算法完成微泡流动的追踪，获得仿体血管的结构图，根据仿体血管结构图的曲线半高宽计算确定仿体血管直径(血管结构图可以看成由很多图像信息组成，将其用曲线的方式表达出来，并取其半高宽来标定血管直径)及对应位置(追踪出来的血管在图像中的坐标点)，从而实现对血管位置分布的精准定位(指在成像平面即中心平面内定位位置分布，由射频数据处理得到的超声图像都是成像平面内的)。
所述奇异值分解从物理的角度将波束合成得到的图像信息用更简而小的几个矩阵的相乘来表示，通过这几个小矩阵来描述超声图像中的特征信息，并从中提取颅脑仿体中血液流动的特征。在这种奇异值分解中，颅脑仿体组织位移因其较高地时空相干性使得大量空间像素在时间尺度上呈现了一致性，即其特征主要通过第一奇异值和奇异向量来描述，而血流回波信号的时空相干性较低，因此奇异值相对要低的很多。采用奇异值分解进行的时空滤波表示为：
Sf＝SVIfV*＝UΔfV*
其中，Sf是滤波后的数据集；If是矩阵滤波器，即第一个对角线元素为零的矩阵，对应于去除颅脑仿体组织运动的奇异值的阈值对角矩阵，S是输入的数据集，U和V是正交矩阵，*表示共轭转置，Δf是对角矩阵。
所述高斯峰值拟合定位通过构造椭圆高斯核K沿着时空滤波后的图像的宽度x和深度z进行图形拟合，然后，选取二维高斯包络的峰值表示在成像平面中的各个造影微泡的中心位置。
所述高斯核K定义为：
其中，Δx是横向像素间距；Δy是纵向像素间距；T表示矩阵转置；
α是常数，由
决定；B是常数，由
决定；σx(y)和σz(y)是
这个对角阵中的变量，可通过利用y方向的单位间距进行定位，具体表示为：
其中，σ0x和σ0z是超声成像设备和超声换能器的性能(σ0x对应于σ0x，0z前面的部分，σ0z对应于σ0x，0z后面的部分)，根据图像数据进行经验估计得到；y和d也是取决于超声成像设备和超声换能器的性能。
在使用马尔科夫链蒙特卡洛算法完成对微泡流动的追踪时，首先将从所有帧图像(指各帧时空滤波后的图像)得到的定位点坐标(指通过定位得到的造影微泡的中心坐标)集Q分解为一些列子集M，其中M＝0，1，...，m(m表示追踪轨迹的数量)，便可得到数据关联结果ω。将虚警结果联合为τ0，而其它被认为可以构成轨迹的正确造影微泡坐标则被连接成轨迹τ1至τm。
其中，P(ω|Q)表示建模数据关联的后验概率，Ω表示总的数据关联集，ωmax表示最大化数据关联集。
上述ωmax是概率最大化情况下众多目标位置关联形成的各个轨迹集，这些轨迹集最后就构成了仿体血管的结构图。
参见图6，根据追踪得到的仿体血管的结构图，选取一段血管区域(即选取一定长度的血管)作为处理区，在该区域标定血管的上下腔壁(首先，标定是指依据结构图中的微泡轨迹信息，判断血管所处的位置。因为这些轨迹信息构成的结构图是有一定宽度的，因微泡是在血管内流动的，这里的宽度就是血管的直径，而标定就是判断结构图中轨迹的上下界面；另外，选取一段血管区域作为处理区是因为后续要测定血流分布的也是这段区域，这个区域可以理解为选取一定长度的血管)，并据此确定血管的中心位置的坐标，将其中一个中心位置作为标定点，然后使用循环移位法将其他的中心位置的坐标移动至同一纵坐标高度，然后对图像信号幅值做归一化处理得到曲线，并通过测量曲线的半高宽估计血管的尺寸(指血管直径)及血管在图像中的位置坐标信息。
步骤4利用互相关算法分析时空滤波后的图像(指序列图像中的每一帧经过时空滤波后的图像)，对整体图像划分ROI区域，划分的原则是考虑了流体流动的入口段效应。其次，将入口和出口段截除后，在最大程度地保留有效的血流回波信号的情况下，对划分ROI区域后的整体图像采用互相关算法进行分析，计算出ROI区域流体流动的情况，得到成像平面的ROI区域内血流的有效峰值速度。
假设所估算段的血流符合层流的条件，根据层流入口段的长度Le和血管管径d以及雷诺数Re之间的比值可确定入口效应存在的距离：
其中，雷诺数由流体的流速、密度、粘性系数和血管的当量直径共同决定。
本发明利用散斑图像分析计算流体速度分布中采用的互相关算法，即选取连续的两帧图像，通过划分矩形判读窗口(在ROI内划分判读窗口)，将两帧图像中对应位置的判读窗口进行互相关运算，得到互相关函数的相关峰相对判读窗口(连续两帧图像中同一位置划分的判读窗口)中心的距离和方向，即为判读窗口所代表的流体微团的位移的大小和方向。具体的判读窗口之间的互相关函数R(m，n)表示为：
其中，连续两帧图像的判读窗口大小均为M×N；a(m，n)和b(m，n)为判读窗口内图像灰度函数(a、b分别对应两个图像各自的判读窗口)；-M≤m≤M，-N≤n≤N。
以上采用的互相关算法计算的对象是整幅时空滤波后的图像，其也具有速度分布，但是这个速度分布是针对血管存在的大概区域的(这里没有对血管精确的位置进行定位)，但是通过与以下步骤5得到的血管内的血流速度分布对比发现，互相关算法得到速度中的中心峰值速度(即速度分布中位于管中心线的最大速度)与步骤5中得到的峰值速度相比，更为接近于理论模拟值，分析认为这是互相关算法固有的计算精度高的特点导致的，因此，这里认为互相关算法计算得到的中心峰值速度比步骤5得到的中心峰值速度更有效，称前者为有效峰值速度。
步骤5根据步骤3得到的血管直径及位置信息，对时空滤波后的图像中的每一帧图像，依据对应的血管位置坐标信息完成截取并扩大到原比例，并采用双三次插值的方法对图像信息进行还原补充，提高其信噪比，然后利用划分三级窗口迭代的快速傅里叶变换相关算法分析得到血管血流速度分布(图4)。
上述扩大到原比例的操作实例如下：
对一根横向血管进行处理，血管长度方向按照原尺寸截取，直径方向按照位置坐标截取相应的宽度，并做图像扩大，扩大至原来时空滤波的图像的纵向尺寸，扩大过程中缺少的像素通过插值补足；扩大并插值后的图像像素间的距离为一个超声波波长时符合信号分析中的香农采样定理。
所述三级窗口迭代的快速傅里叶变换相关算法，用于对截取并插值放大后得到的血管区域内的图像进行分析，在迭代的过程中需要逐级减小判读窗口尺寸，直至判读窗口的尺寸满足空间分辨率和测量精度的要求，从而精确地估计得到血管直径范围内血流速度梯度分布，这里采用的判读窗口逐级递减的方式需要满足第一级判读窗口尺寸的四分之一大于连续两帧图像中微泡团的最大位移(最大位移若超过窗口的四分之一，窗口之间的相关性就会降低，计算就会产生较大误差)。
所述三级窗口迭代的快速傅里叶变换相关算法的数学表达式具体描述为：
b(m，n)＝a(m，n)*s(m，n)
其中，b(m，n)表示t0+Δt时刻判读窗口的图像灰度函数，a(m，n)表示t0时刻判读窗口的图像灰度函数，s(m，n)表示位移函数，*表示卷积；
利用卷积定理，得：
B(U，V)＝A(U，V)*S(U，V)
其中，B(U，V)表示b(m，n)的傅里叶变换，A(U，V)表示a(m，n)的傅里叶变换，S(U，V)表示s(m，n)的傅里叶变换，*表示卷积；
令A*(U，V)为A(U，V)的复共轭，则有：
其中，φ(U，V)＝A*(U，V)B(U，V)，只要求得φ(U，V)的反变换
就可以得到
的峰值位置，然后识别并拟合相关系数峰值，确定流体微团的运动速度(检测到峰值的位置，则该位置离开判读窗口中心的距离为微泡团在该窗口的位移)。
步骤6针对步骤4和步骤5得到的血流速度分布数据，分别选取其优势部分进行复合，即得到更为精确的血管血流速度场的估计结果。
从流体速度分布的分析计算结果上看，互相关算法较快速傅里叶变换算法的计算精度更高，且实际测算结果显示互相关算法的中心峰值速度计算结果与理论值有更好的吻合度。但在亚毫米级血管(0.1～1mm)内完成计算时，单一窗口无法有效地识别像素间的变化，而采用尺寸逐级递减的判读窗口进行快速傅里叶变换，可以较为细致、精确地测算得到血管内部血流速度梯度分布。因此，将互相关算法得到的中心峰值速度(有效峰值速度)替换该血流速度梯度分布中的中心峰值速度，即得到了结合两种方法(步骤4和步骤5)的血管血流速度场的估计结果(计算值)。
实验结果表明，两种方法(步骤4和步骤5)结合后所得计算值与以下模拟值契合较好(图5)，实现了对比验证。
理论模拟值可通过成像平面内任意一点(-d，y)处流速的理论表达式得到：
v(-d，y)＝(1-(d2+y2)/a2)v0
d表示扫描平面与中心平面之间的距离(成像平面就在中心平面，因此d＝0)；v0表示血管中心线最大速度；a表示为血管管径的二分之一；y为成像平面上径向上分布的点。
总之，本发明为脑小血管血流速度场的估计提供一种可行的方法，能够在经颅超声成像结果分辨率低、血管边界无法直接辨别的情况下，实现血管血流速度场的精确估计。
Claims (6)
1.一种经颅超声估计血管血流速度场的方法，其特征在于：该估计血管血流速度场的方法包括以下步骤：
1)利用超声扫描颅脑内部区域，并采集超声回波信号；
2)通过对超声回波信号进行时空滤波、成像及图像定位和追踪，得到超声成像区域的血管结构图；
所述步骤2)具体包括以下步骤：对采集的超声回波信号通过波束合成及奇异值分解进行时空滤波，时空滤波后只保留血流回波信号，然后生成时空滤波后的超声序列图像，对该超声序列图像采用高斯峰值拟合定位得到对应图像中造影微泡的位置，并对造影微泡的流动进行追踪；
3)选取血管结构图中的血管，并推断该血管的位置坐标信息；
所述步骤3)具体包括以下步骤：采用测量血管结构图曲线半高宽的方法确定血管直径及位置坐标；
4)根据血管的直径及位置坐标从时空滤波后通过成像形成的超声序列图像中提取血管内血流的图像，并对提取的相邻帧血管内血流的图像采用划分多级判读窗口进行迭代的快速傅里叶变换相关算法进行分析，得到血管血流速度梯度分布；
所述估计血管血流速度场的方法还包括以下步骤：采用互相关算法对步骤2)中利用时空滤波保留的血流回波信号生成的超声序列图像进行速度分布分析，得到有效峰值速度，将有效峰值速度与步骤4)得到的血流速度梯度分布复合。
2.根据权利要求1所述一种经颅超声估计血管血流速度场的方法，其特征在于：所述步骤1)具体包括以下步骤：采用可编程超声平台控制超声换能器向颅内发射低频线性调频平面波，并由可编程超声平台采集所述超声换能器接收到的回波信号。
3.根据权利要求1所述一种经颅超声估计血管血流速度场的方法，其特征在于：所述血管选自直径≤1mm的脑小血管。
4.根据权利要求1所述一种经颅超声估计血管血流速度场的方法，其特征在于：所述步骤4)中，对步骤2)中利用时空滤波保留的血流回波信号生成的超声序列图像进行血管内超声粒子图像信息的截取、扩大及插值，得到血管内血流的图像信息。
5.一种经颅超声估计血管血流速度场的装置，其特征在于：该估计血管血流速度场的装置包括超声换能器(3)、可编程超声平台(1)和计算机；
所述超声换能器(3)用于在可编程超声平台(1)控制下通过发射超声波扫描颅脑内部区域并接收超声回波信号；
所述计算机包括经颅超声估计血管血流速度场的软件系统，该软件系统包括血管标定模块及血流速度梯度分布计算模块；
所述血管标定模块用于对超声回波信号进行时空滤波、成像、图像定位和追踪后采用测量血管结构图曲线半高宽的方法确定血管直径及位置坐标；其中，对采集的超声回波信号通过波束合成及奇异值分解进行时空滤波，时空滤波后只保留血流回波信号，然后生成时空滤波后的超声序列图像，对该超声序列图像采用高斯峰值拟合定位得到对应图像中造影微泡的位置，并对造影微泡的流动进行追踪；
所述血流速度梯度分布计算模块用于根据血管的直径及位置坐标从时空滤波后通过成像形成的超声序列图像中提取血管内血流的图像，并对提取的相邻帧血管内血流的图像采用划分多级判读窗口进行迭代的快速傅里叶变换相关算法进行血流速度梯度分布分析；
所述软件系统还包括血流速度梯度分布优化模块，该模块用于采用互相关算法对时空滤波后通过成像形成的超声序列图像进行速度分布分析，以及用于将采用互相关算法分析得到的有效峰值速度与采用划分多级判读窗口进行迭代的快速傅里叶变换相关算法计算得到的血流速度梯度分布复合。
6.根据权利要求5所述一种经颅超声估计血管血流速度场的装置，其特征在于：所述估计血管血流速度场的装置还包括颅脑小血管仿血液流动系统，用于模拟脑小血管内血液的流动。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110092110.5A CN112932539B (zh) | 2021-01-23 | 2021-01-23 | 一种经颅超声估计血管血流速度场的装置及方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110092110.5A CN112932539B (zh) | 2021-01-23 | 2021-01-23 | 一种经颅超声估计血管血流速度场的装置及方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112932539A CN112932539A (zh) | 2021-06-11 |
| CN112932539B true CN112932539B (zh) | 2023-05-09 |
Family
ID=76236126
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110092110.5A Active CN112932539B (zh) | 2021-01-23 | 2021-01-23 | 一种经颅超声估计血管血流速度场的装置及方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112932539B (zh) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113724562B (zh) * | 2021-07-20 | 2023-04-14 | 西安交通大学 | 一种用于经颅超声扫查的仿体颅脑模型及其制备方法 |
| CN115770063A (zh) * | 2021-09-08 | 2023-03-10 | 深圳开立生物医疗科技股份有限公司 | 一种血管成像方法、装置、电子设备及可读存储介质 |
| CN114098817A (zh) * | 2021-11-18 | 2022-03-01 | 西安建筑科技大学 | 一种高帧率超声图像血管壁运动细节追踪方法、系统、设备和可读存储介质 |
| CN114366163B (zh) * | 2022-01-11 | 2023-08-25 | 深圳市德力凯医疗设备股份有限公司 | 基于快速扫描的脑血流数据采集方法、系统及智能终端 |
| CN114419016B (zh) * | 2022-01-25 | 2022-11-04 | 华润武钢总医院 | 用于甲状腺结节良恶性鉴别诊断的超声超分辨方法及系统 |
| CN114533122B (zh) * | 2022-03-11 | 2024-08-23 | 清华大学 | 一种超声微血流成像的信号处理方法和系统 |
| CN114848013B (zh) * | 2022-04-29 | 2024-11-26 | 西安交通大学 | 一种利用超声微泡次谐波测量颅内绝对血压分布图的方法 |
| CN117562577B (zh) * | 2023-12-14 | 2024-04-26 | 华润武钢总医院 | 一种基于奇异值分解滤波的超声矢量流速成像方法及系统 |
| CN118351079B (zh) * | 2024-04-18 | 2025-03-28 | 山东探微医疗技术有限公司 | 毛细血管光滑程度判定方法及装置 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2018222724A1 (en) * | 2017-05-31 | 2018-12-06 | Mayo Foundation For Medical Education And Research | Methods for super-resolution ultrasound imaging of microvessels |
| CN107753062B (zh) * | 2017-11-27 | 2021-03-16 | 西安交通大学 | 基于马尔科夫链蒙特卡洛多目标追踪的经颅超声脑血管造影超分辨率成像方法 |
| CN107714091B (zh) * | 2017-11-27 | 2019-12-20 | 西安交通大学 | 经颅低频超声线性调频脉冲逆转微泡成像方法 |
| CN108836392B (zh) * | 2018-03-30 | 2021-06-22 | 中国科学院深圳先进技术研究院 | 基于超声rf信号的超声成像方法、装置、设备及存储介质 |
| CN113180734A (zh) * | 2018-12-27 | 2021-07-30 | 深圳迈瑞生物医疗电子股份有限公司 | 一种超声血流成像方法及系统 |
-
2021
- 2021-01-23 CN CN202110092110.5A patent/CN112932539B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN112932539A (zh) | 2021-06-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112932539B (zh) | 一种经颅超声估计血管血流速度场的装置及方法 | |
| Harput et al. | 3-D super-resolution ultrasound imaging with a 2-D sparse array | |
| CN110381845B (zh) | 具有用于导出成像数据和组织信息的神经网络的超声成像系统 | |
| Ricci et al. | Real-time blood velocity vector measurement over a 2-D region | |
| Voorneveld et al. | 4-D echo-particle image velocimetry in a left ventricular phantom | |
| US8992426B2 (en) | Feedback in medical ultrasound imaging for high intensity focused ultrasound | |
| KR102006035B1 (ko) | 초음파 촬영을 이용한 관 특성기술 | |
| US20140147013A1 (en) | Direct echo particle image velocimetry flow vector mapping on ultrasound dicom images | |
| US10332250B2 (en) | Three-dimensional cavitation quantitative imaging method for microsecond-resolution cavitation spatial-temporal distribution | |
| US20080015440A1 (en) | Echo particle image velocity (EPIV) and echo particle tracking velocimetry (EPTV) system and method | |
| KR20180013956A (ko) | 단일 추적 위치 전단파 탄성 이미징을 위한 방법, 시스템 및 컴퓨터 프로그램 제품 | |
| Foroozan et al. | Microbubble localization for three-dimensional superresolution ultrasound imaging using curve fitting and deconvolution methods | |
| Wei et al. | High frame rate volumetric imaging of microbubbles using a sparse array and spatial coherence beamforming | |
| US20220386996A1 (en) | Ultrasonic shearwave imaging with patient-adaptive shearwave generation | |
| Salles et al. | Full 3-D transverse oscillations: a method for tissue motion estimation | |
| Liang et al. | Velocity field estimation in transcranial small vessel using super-resolution ultrasound imaging velocimetry | |
| Caudoux et al. | Curved toroidal row column addressed transducer for 3D ultrafast ultrasound imaging | |
| Wei et al. | High-frame-rate volumetric porcine renal vasculature imaging | |
| CN103340620B (zh) | 一种管壁应力相位角的测量方法和系统 | |
| EP3755230B1 (en) | Method and apparatus for simultaneous 4d ultrafast blood flow and tissue doppler imaging of the heart and retrieving quantification parameters | |
| CN112823295B (zh) | 检测流动不稳定性的方法和装置 | |
| Wang et al. | Plane-wave ultrasound imaging based on compressive sensing with low memory occupation | |
| Zhou et al. | Transcranial ultrafast ultrasound Doppler imaging: A phantom study | |
| Jansen et al. | Increasing abdominal aortic aneurysm curvature visibility using 3D dual probe bistatic ultrasound imaging combined with probe translation | |
| Azhari et al. | Ultrasound Imaging |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |