CN109615677B - 一种基于低采样率b超图像计算热应变分布的方法 - Google Patents
一种基于低采样率b超图像计算热应变分布的方法 Download PDFInfo
- Publication number
- CN109615677B CN109615677B CN201910112033.8A CN201910112033A CN109615677B CN 109615677 B CN109615677 B CN 109615677B CN 201910112033 A CN201910112033 A CN 201910112033A CN 109615677 B CN109615677 B CN 109615677B
- Authority
- CN
- China
- Prior art keywords
- ultrasonic
- signal
- image
- sequence
- frame
- 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
- 238000000034 method Methods 0.000 title claims abstract description 94
- 238000005070 sampling Methods 0.000 title claims abstract description 68
- 239000000523 sample Substances 0.000 claims abstract description 44
- 238000001914 filtration Methods 0.000 claims abstract description 37
- 230000008569 process Effects 0.000 claims abstract description 6
- 238000001514 detection method Methods 0.000 claims abstract description 4
- 230000009467 reduction Effects 0.000 claims abstract description 4
- 238000002604 ultrasonography Methods 0.000 claims description 35
- 238000012545 processing Methods 0.000 claims description 27
- 238000002679 ablation Methods 0.000 claims description 19
- 230000004044 response Effects 0.000 claims description 19
- 230000003595 spectral effect Effects 0.000 claims description 19
- 238000006073 displacement reaction Methods 0.000 claims description 10
- 238000001228 spectrum Methods 0.000 claims description 7
- 230000001502 supplementing effect Effects 0.000 claims description 7
- 239000013589 supplement Substances 0.000 claims description 6
- 238000003780 insertion Methods 0.000 claims description 5
- 230000015572 biosynthetic process Effects 0.000 claims description 4
- 230000037431 insertion Effects 0.000 claims description 4
- 238000012544 monitoring process Methods 0.000 claims description 4
- 238000007674 radiofrequency ablation Methods 0.000 claims description 4
- 238000003786 synthesis reaction Methods 0.000 claims description 4
- 238000000608 laser ablation Methods 0.000 claims description 3
- 238000003672 processing method Methods 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 238000004891 communication Methods 0.000 claims description 2
- 238000013461 design Methods 0.000 claims description 2
- 230000009466 transformation Effects 0.000 claims description 2
- 238000012285 ultrasound imaging Methods 0.000 claims 1
- 238000004364 calculation method Methods 0.000 abstract description 20
- 238000003384 imaging method Methods 0.000 abstract description 8
- 238000000605 extraction Methods 0.000 abstract 2
- 210000001519 tissue Anatomy 0.000 description 24
- 238000010438 heat treatment Methods 0.000 description 13
- 238000009529 body temperature measurement Methods 0.000 description 10
- 230000008859 change Effects 0.000 description 9
- 230000008901 benefit Effects 0.000 description 4
- 230000001186 cumulative effect Effects 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 210000000577 adipose tissue Anatomy 0.000 description 3
- 210000005228 liver tissue Anatomy 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 206010028980 Neoplasm Diseases 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000009191 jumping Effects 0.000 description 2
- 238000002560 therapeutic procedure Methods 0.000 description 2
- 238000010317 ablation therapy Methods 0.000 description 1
- 238000001816 cooling Methods 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000008506 pathogenesis Effects 0.000 description 1
- 230000002285 radioactive effect Effects 0.000 description 1
- 230000001105 regulatory effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000001356 surgical procedure Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/48—Diagnostic techniques
- A61B8/485—Diagnostic techniques involving measuring strain or elastic properties
-
- 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
- A61B8/5223—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of medical diagnostic data for extracting a diagnostic or physiological parameter from medical diagnostic data
-
- G06T5/70—
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10016—Video; Image sequence
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10132—Ultrasound image
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Molecular Biology (AREA)
- Heart & Thoracic Surgery (AREA)
- Veterinary Medicine (AREA)
- Public Health (AREA)
- Biophysics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Animal Behavior & Ethology (AREA)
- Surgery (AREA)
- General Physics & Mathematics (AREA)
- Physiology (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
- Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
Abstract
本发明公开了一种基于低采样率B超声图像计算热应变分布的方法，属于超声成像领域。本发明基于B超输出的图像信号，具体针对由超声波探头输出的RF(射频)信号经正交检波、降低采样率而获得的信号；对每条纵向超声线依次进行数据点间隔补零、低通滤波、群延迟补偿和抽取，从而增加B超图像的纵向采样率；对每条横向超声线依次进行一次低通滤波、间隔补零、二次低通滤波、群延迟补偿和抽取，从而增加B超图像的横向采样率。对随时间变化的超声图像序列，本发明将每帧B超图像进行如上处理，对相邻两帧图像的对应区域分别空间分割后进行二维幅度变迹，再利用相关性计算和差分滤波，获得热应变分布图像，最终得到随时间变化的热应变图像序列。
Description
技术领域
本发明涉及一种基于低采样率B超声图像计算热应变分布的方法，属于医学图像处理、医学超声成像领域。
背景技术
随着医学科技的发展和对肿瘤发病机制的深入研究，当前临床上已可通过外科手术、放射线和化学药物等手段治疗部分肿瘤；但是，经此类治疗后的患者生活质量通常较差。近年来，为在提高治愈率的同时改善患者的生活质量，微创介入治疗(包括射频消融、微波消融、超声聚焦消融等)以其微创/无创、安全的优点受到临床医学界的广泛重视。这类技术利用不同的热源对生物组织进行加热，消融局部病变组织。基于上述原理的热消融治疗对于热剂量的控制有着严格的要求。因此，对被加热区域的组织进行温度监控非常重要。
超声波测温技术是近年来快速发展的一种新型测温技术。超声测温的基本原理一般又分为三类：(1)基于组织受热后声速的变化；通常测温误差较大，且容易出现测温点误判，近年来已经较少使用；(2)基于组织背散射能量的变化；可以相对准确地判断温度变化区域，但测温精度不高；(3)基于组织热应变分布的变化；该方法具有较强的精确测温潜力，且可以准确地判断温度变化区域；其原理是：生物组织在温度变化时热胀冷缩，因而在热疗时的升温或降温过程中产生出持续的热应变分布；将热应变分布与不同组织的性质(包括声速、热容量、热扩散系数、热弛豫时间等)相结合，即可能获得准确的组织温度分布。
组织热应变分布可通过B超图像序列进行计算。但商用B超成像仪限于硬件传输速率和系统带宽，其输出数据一般由原始射频信号(采样率通常40MHz左右)经波束合成、正交解调和降采样率后形成；这种数据因具有较低的采样率(采样率通常降至5MHz以下)，难以直接通过其获得准确的热应变分布图像。以图像上的一个被加热点为例，经低分辨率B超图像序列计算获得的热应变-时间曲线通常具有严重的“台阶状”，因而难以反映组织温度变化的实际过程，用于组织测温会导致严重的误差。这在该技术的发展中是一个亟待解决的问题。
发明内容
本发明的目的在于提供一种基于低采样率数字B超图像，精确计算生物组织内热应变二维分布的方法。
为实现上述目的，本发明采取的技术方案是：一种基于低采样率B超声图像计算热应变分布的方法，该方法包括以下步骤：
步骤一、在生物组织热消融过程中采集超声图像信号，利用所选探头结合B超仪主机对治疗区域进行监控，采集并输出按时间排列的数字超声图像序列；
步骤二、对低采样率B超图像序列进行纵向插值，获得高采样率图像序列；
步骤三、对纵向插值后的图像序列进行横向插值，使每帧具有更多的超声线；
步骤四、根据插值后的B超图像序列，计算热应变实时分布图像。
进一步的，所述步骤一的具体步骤包括：
(1)根据应用需求，确定B超探头的类型、波束控制方式、几何形状和阵元数目；
(2)对生物组织进行局部热消融；
(3)利用所选探头结合B超仪主机对治疗区域进行监控，采集并输出按时间排列的数字超声图像序列；该序列中每帧图像由探头接收的原始超声回波信号经正交检波、降低采样率、波束合成后形成。
更进一步的，所述B型超声成像仪所使用的超声探头的型号与工作频率根据应用需求选用，超声探头类型包括凸阵探头、线阵探头、相控阵探头；波束控制方式包括线扫、相控阵、机械扇扫和面阵；几何形状包括弧形、圆形和矩形；阵元的数目为单阵元或多阵元。
更进一步的，所述对生物组织进行局部热消融的方法包括微波消融、射频消融、超声消融、激光消融和红外线消融。
进一步的，所述步骤二中的纵向插值可以采取两种有效的信号插值方法——“时域补零后低通滤波”法和“频域插零后逆傅里叶变换”法，应用两种不同的方法，其中，
所述“时域补零后低通滤波”法的具体步骤包括：
(1)设共有N帧按时间顺序排列的B超图像，每帧含有L条横向排列的超声线；每条超声线含有M个纵向分布的实数数据点；选取任意一帧B超图像；
(2)对任意一条超声线上的所有数据点，按照其距离B超探头的远近依次用自然数从1开始标记序号，如第m个数据点则标记为m；将每条超声线上的数据点按标记值的奇偶依次分为两组，以奇数组为实部、偶数组为虚部构建复数信号；
(3)对所得复信号进行快速傅里叶变换，根据所得谱图获得数字截止频率ωc1，使0～ωc1范围内频谱能量不小于总频谱能量；能量由所有谱线幅度的平方和求得；
(4)根据纵向采样率所需提升的有理数倍数I1，将I1表示为I1＝P1/Q1，其中P1和Q1均为正整数；设计一个具有线性相位性质的有限长单位抽样响应低通数字滤波器，其通带截止频率为ωc1/P1，群延迟为正整数τ1；
(5)对(2)中所得复信号，在每个数据点后补充P1个0，在信号末尾再补充τ1个0；利用(4)中所得滤波器对该信号进行滤波，对结果信号去除其最前面τ1个数据点；
(6)对(5)中所得信号，选择序号为1～Q1的任意数据点为保留数据点，随后每隔Q1-1保留一个数据点，其余数据点均舍弃；至此，得到(2)中超声线的纵向插值信号；
(7)对(1)中选中的图像，按(2)-(6)对所有L条纵向超声线进行处理，所得的L条插值信号线按原顺序进行组合，得到纵向插值后的一帧高采样率B超图像；
(8)对(1)中的所有N帧B超图像，按(2)-(7)的步骤对所有帧进行处理，得到包含N帧的纵向插值图像序列，每帧含有L条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点，M1为接近M×I1/2的正整数；
所述“频域插零后逆傅里叶变换”法的具体步骤包括：
(Ⅰ)设共有N帧按时间顺序排列的B超图像，每帧含有L条横向排列的超声线；每条超声线含有M个纵向分布的实数数据点；选取任意一帧B超图像；
(Ⅱ)对任意一条超声线上的所有数据点，按照其距离B超探头的远近依次用自然数从1开始标记序号，如第m个数据点则标记为m；将每条超声线上的数据点按标记值的奇偶依次分为两组，以奇数组为实部、偶数组为虚部构建复数信号；
(Ⅲ)对所得复信号进行快速傅里叶变换，得到点数为T0的频谱序列(T0＝M/2)，按照变换前复信号的先后顺序，从1开始标记序号，如第k个数据点则标记为k，k取值范围为1～T0，该点对应的值表示为X1(k)；
(Ⅳ)根据纵向采样率所需提升的有理数倍数I1，将I1表示为I1＝P1/Q1，其中P1和Q1均为正整数；
(Ⅴ)对(Ⅲ)中所得的频谱序列在频域插零。令长度为T1的空序列Y1(n)，其中T1＝T0*P1，n取值1～T1，Y1(n)表示对应n点的值。当T0为奇数，在n∈[1,(T0+1)/2]时，Y1(n)＝P1*X1(n)。在n∈[(T0+1)/2+1,T1-(T0-1)/2]，Y1(n)＝0。在n∈[T1-(T0-1)/2+1,T1]，Y1(n)＝P1*X1(n-T1+T0)；当T0为偶数，在n∈[1,T0/2]时，Y1(n)＝P1*X1(n)。n＝T0/2+1时，Y1(n)＝0.5*P1*X1(n)。在n∈[T0/2+2,T1-T0/2]，Y1(n)＝0。在n＝T1-T0/2+1时，Y1(n)＝0.5*P1*X1(T0/2+1)。在n∈[T1-T0/2+2,T1]，Y1(n)＝P1*X1(n-T1+T0)。
(Ⅵ)对(Ⅴ)中得到的频域插零后的Y1(n)序列，进行傅里叶逆变换，对所得信号，选择序号为1～Q1的任意数据点为保留数据点，随后每隔Q1-1保留一个数据点，其余数据点均舍弃；至此，得到(Ⅱ)中超声线的纵向插值信号；
(Ⅶ)对(Ⅰ)中选中的图像，按(Ⅱ)-(Ⅵ)对所有L条纵向超声线进行处理，所得的L条插值信号线按原顺序进行组合，得到纵向插值后的一帧高采样率B超图像；
(Ⅷ)对(Ⅰ)中的所有N帧B超图像，按(Ⅱ)-(Ⅶ)的步骤对所有帧进行处理，得到包含N帧的纵向插值图像序列，每帧含有L条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点，M1为接近M×I1/2的正整数；
更进一步的，所述“时域补零后低通滤波”法中纵向插值所用低通滤波器具有线性相位、有限长单位抽样响应、通带最大衰减不超过3dB、阻带最小衰减不小于20dB、截止频率内信号能量不低于总能量80％、群延迟为正整数的特点。
进一步的，所述步骤三中的横向插值可以采取两种有效的信号插值方法——“时域补零后低通滤波”法和“频域插零后逆傅里叶变换”法，应用两种不同的方法，其中，
所述“时域补零后低通滤波”法的具体步骤包括：
(1)共有N帧按时间顺序排列的纵向插值后B超图像，每帧含有L条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点；选取任意一帧B超图像；
(2)将每条超声线的第k个复数数据点按照超声线的排列顺序进行组合，得到第k条横向超声信号；每帧图像共包含M1条横向超声信号，每条长度为L；
(3)对第k条横向超声信号按照从左到右依次用自然数从1开始标记序号，如第r个数据点则标记为r；对该横向信号进行快速傅里叶变换，根据谱图获得数字截止频率ωc2，使0～ωc2范围内频谱能量不小于总频谱能量；
(4)设计一个具有线性相位性质的有限长单位抽样响应一次低通数字滤波器，通带截止频率为ωc2，群延迟为正整数τ2；
(5)对(3)中选中的横向超声信号，在其末尾补充τ2个0后，利用(4)中所得滤波器对其进行滤波；将滤波所得信号最前面τ2个数据点舍弃；得到第k条防混叠横向超声信号；
(6)根据横向采样率所需提升的有理数倍数I2，将I2表示为I2＝P2/Q2，其中P2和Q2均为正整数；设计一个具有线性相位性质的有限长单位抽样响应二次低通数字滤波器，其通带截止频率为ωc2/P2，群延迟为正整数τ3；
(7)对(5)中所得复信号，在每个数据点后补充P2个0，在信号末尾再补充τ3个0；利用(6)中所得滤波器对该信号进行滤波，对结果信号去除其最前面τ3个数据点；
(8)对(7)中所得信号，选择序号为1～Q2的任意数据点为保留数据点，随后每隔Q2-1保留一个数据点，其余数据点均舍弃；至此得到(3)中横向超声信号的插值结果；
(9)对(1)中选中的图像，按(2)-(8)对所有M1条横向超声线进行处理，所得的M1条插值信号线按原顺序排列，得到纵、横向均插值后的一帧高采样率B超图像；
(10)对(1)中的所有N帧B超图像，按(2)-(9)的步骤对所有帧进行处理，得到包含N帧的图像序列，每帧含有L1条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点，L1为接近L×I2的正整数。
所述“频域插零后逆傅里叶变换”法的具体步骤包括：
(Ⅰ)共有N帧按时间顺序排列的纵向插值后B超图像，每帧含有L条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点；选取任意一帧B超图像；
(Ⅱ)将每条超声线的第k个复数数据点按照超声线的排列顺序进行组合，得到第k条横向超声信号；每帧图像共包含M1条横向超声信号，每条长度为L；
(Ⅲ)对第k条横向超声信号按照从左到右依次用自然数从1开始标记序号，如第r个数据点则标记为r；对该横向信号进行快速傅里叶变换，根据谱图获得数字截止频率ωc2，使0～ωc2范围内频谱能量不小于总频谱能量；
(Ⅳ)设计一个具有线性相位性质的有限长单位抽样响应一次低通数字滤波器，通带截止频率为ωc2，群延迟为正整数τ2；
(Ⅴ)对(Ⅲ)中选中的横向超声信号，在其末尾补充τ2个0后，利用(4)中所得滤波器对其进行滤波；将滤波所得信号最前面τ2个数据点舍弃；得到第k条防混叠横向超声信号；
(Ⅵ)对(Ⅴ)所得横向超声信号进行快速傅里叶变换，得到点数为L的频谱序列，按照变换前复信号的先后顺序，从1开始标记序号，如第k个数据点则标记为k，k取值范围为1～L，该点对应的数值表示为X2(k)；
(Ⅶ)根据横向采样率所需提升的有理数倍数I2，将I2表示为I2＝P2/Q2，其中P2和Q2均为正整数；
(Ⅷ)对(Ⅵ)中所得的频谱序列在频域插零。令长度为T2的空序列Y2(n)，其中T2＝L*P2，n取值1～T2，Y2(n)表示对应n点的值。当L为奇数，在n∈[1,(L+1)/2]时，Y2(n)＝P2*X2(n)。在n∈[(L+1)/2+1,T2-(L-1)/2]，Y2(n)＝0。在n∈[T2-(L-1)/2+1,T2]，Y2(n)＝P2*X2(n-T2+L)；当L为偶数，在n∈[1,L/2]时，Y2(n)＝P2*X2(n)。n＝L/2+1时，Y2(n)＝0.5*P2*X2(n)。在n∈[L/2+2,T2-L/2]，Y2(n)＝0。在n＝T2-L/2+1时，Y2(n)＝0.5*P2*X2(L/2+1)。在n∈[T2-L/2+2,T2]，Y2(n)＝P2*X2(n-T2+L)。
(Ⅸ)对(Ⅷ)中得到的在频域插零后的Y2(n)序列，进行傅里叶逆变换，对所得信号，选择序号为1～Q2的任意数据点为保留数据点，随后每隔Q2-1保留一个数据点，其余数据点均舍弃；至此，得到(Ⅱ)中超声线的横向插值信号；
(Ⅹ)对(Ⅰ)中选中的图像，按(Ⅱ)-(Ⅸ)对所有M1条横向超声线进行处理，所得的M1条插值信号线按原顺序排列，得到纵、横向均插值后的一帧高采样率B超图像；
(Ⅺ)对(Ⅰ)中的所有N帧B超图像，按(Ⅱ)-(Ⅹ)的步骤对所有帧进行处理，得到包含N帧的图像序列，每帧含有L1条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点，L1为接近L×I2的正整数。
更进一步的，所述“时域补零后低通滤波”法中横向插值所用一次低通滤波器具有线性相位、有限长单位抽样响应、通带最大衰减不超过3dB、阻带最小衰减不小于20dB、截止频率内信号能量不低于总能量80％、群延迟为正整数的特点；二次低通滤波器具有线性相位、有限长单位抽样响应、通带最大衰减不超过3dB、阻带最小衰减不小于20dB、群延迟为正整数的特点。
进一步的，所述步骤四中计算热应变图像时，由每条纵向信号线获得热应变信号通过滑动平均、提前补零、数字差分滤波、群延迟补偿四个步骤进行。
进一步的，所述步骤四的具体步骤包括：
(1)选择任意一帧图像为参考帧，其后续一帧图像为目标帧；定义图像中的横向坐标为i，纵向坐标为j，i为整数，取值范围1～L1，j为整数，取值范围1～M1；选择L2为5～4L1/5之间的任意奇数自然数，M2为5～4M1/5之间的任意奇数自然数，L3为3～4L2/5之间的任意奇数自然数，M3为3～4M2/5之间的任意奇数自然数；其中L1为B超图像经插值后、横向分布的超声信号线的条数，M1为B超图像经插值后，纵向分布的超声信号线的数据点数；
(2)在目标帧中选择一个点，其坐标为i＝l，j＝m，其中l为正整数，其取值范围为(L2+1)/2～L1-(L2-1)/2，m为正整数，其取值范围为(M2+1)/2～M1-(M2-1)/2；以该点为中心选择i＝l-(L3-1)/2～l+(L3-1)/2、j＝m-(M3-1)/2～m+(M3-1)/2的矩形区域，将该区域内的二维超声信号记为Uo(l,m)；
(3)在参考帧中选择一个点，其坐标为i＝l+p，j＝m+q，其中p为整数，其取值范围为-(L2-L3)/2～(L2-L3)/2，q为整数，其取值范围为-(M2-M3)/2～(M2-M3)/2；以该点为中心选择i＝l+p-(L3-1)/2～l+p+(L3-1)/2、j＝m+q-(M3-1)/2～m+q+(M3-1)/2的矩形区域，将该区域内的二维超声信号取共轭，记为Ur(l,m；p,q)；
(4)设计一个中心位于i＝i0，j＝j0、横向尺寸为L3、纵向尺寸为M3的二维幅度变迹窗函数，可选用如下设计形式：
其中，
(5)将(4)中窗函数的中心(i0,j0)设置于(l,m)，使其与二维信号Uo(l,m)相乘，得到目标帧中以(l,m)为中心的二维变迹加权信号Vo(l,m)；将(4)中窗函数的中心(i0,j0)设置于(l+p,m+p)，使其与二维信号Uo(l,m)相乘，得到参考帧中以(l+p,m+p)为中心的二维变迹加权信号Vr(l,m；p,q)；
(6)对(2)中选定的每一个l和m的具体值，改变p和q的取值，使参数组合(p,q)遍历二维取值空间[-(L2-L3)/2～(L2-L3)/2,-(M2-M3)/2～(M2-M3)/2]，寻找使如下相关函数取值最大的(p,q)组合(plm,qlm)：
(7)plm和qlm即分别为位于i＝l、j＝m的组织由于热膨胀导致的横、纵向位移；对(2)中选中的(l,m)参数组合进行变化，使其遍历二维空间
[(L2+1)/2～L1-(L2-1)/2,(M2+1)/2～M1-(M2-1)/2]，可获得plm和qlm的二维分布图象Dp(l,m)和Dq(l,m)；
(8)取(7)中获得的纵向位移图像Dq(l,m)；其包含L1-L2+1条横向分布的信号线，每条信号线在纵向的长度为M1-M2+1；对其每一条信号线应用K点滑动平均，K为任意不小于3的奇数，经处理后纵向位移图像的纵向长度减小为M1-M2-K+2；
(9)设计一个具有线性相位性质的有限长单位抽样响应数字差分滤波器，其最大波纹不超过3dB，群延迟为正整数τ4；对(8)中的任意一条信号线，在其末尾补充τ4个0，再利用所得差分滤波器对其滤波，将所得结果最前方τ4个数据点舍弃；至此，即得到该条信号线对应的热应变信号；
(10)对(8)中的纵向位移图像，针对每条信号线均按(9)进行处理，可获得当前图像帧对应的组织热应变的二维分布图像；
(11)将(1)中选择的目标帧设定为新的参考帧，其后续一帧作为新的目标帧，重复(2)-(10)中的处理方法；依此类推，直至所有超声图像处理完毕，即可得到N-1个按时间排列的热应变二维分布图像。
本发明与已有的基于B超图像计算生物组织应变的方法相比，可应用较低采样率的B超图像即获得较高精度的应变分布图像。本发明所采用的B超图像信号的所有数据点为实数，且经过降采样率处理；实际可进行处理的信号不局限于此。可能存在的形式还有降采样率后的复数图像(则可省略步骤二的(2)/(Ⅱ)中构建复信号的过程)、未降采样率的超声射频(RF)信号(则可通过本方法进一步提升采样率)等。步骤二、三中设计的两种不同的信号插值方法，经验证有相接近的插值效果，应用于计算热应变分布也有相接近的结果。步骤四中设计的二维幅度变迹函数也不仅限于所列出的形式，还可能存在二维三角形变迹函数、二维凯泽窗变迹函数、以迭代优化形式设计的二维幅度变迹窗函数等。此外，对低分辨率B超图像进行插值时，也可只进行纵向插值或只进行横向插值，但所得热应变计算结果精度会有所下降。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的方式及实施例，均应属于本发明的保护范围。
综上，本发明具有如下显著优点：
(1)已有利用低分辨率B超图像计算热应变分布的方法中，所得热应变随时间的变化呈现“阶梯状”，误差通常较大。本发明提出的方法可有效地克服该问题，提高热应变分布的计算精度；
(2)已有利用高分辨率B超图像计算热应变分布的方法中，超声图像的采集对B超硬件的传输速率较高，要求的信号带宽也较大。本发明提供的方法可有效地克服该问题，降低高精度热应变测量对B超系统硬件的需求；
(3)已有利用B超图像计算生物组织应变分布的方法中，参考帧与目标帧之间进行相关性计算前均未通过幅度变迹处理，不容易突出两个相关区域中心部位的贡献，导致相关性运算结果中的噪声较大，热应变计算误差也同样较大。本发明通过二维幅度变迹处理，强化了两个相关区域中心部位的贡献，可降低热应变计算误差。
附图说明
图1为基于低采样率B超图像计算热应变分布的流程；
图2为生物组织加热过程中采集B超图像和温度数据的示意图；其中：图2a为数据采集系统，附图标记为：1-超声诊断仪，2-成像探头，3-生物组织，4-加热设备，5-加热探针；图2b和图2c为一种新型结构的超声波聚焦探头的平面图和立体图；
图3a为本发明实施例1中，第1帧未经插值的低采样率B超复值图像的幅度；
图3b为本发明实施例1中，低采样率B超图像计算所得热应变图像序列的第130帧；
图3c为本发明实施例1中，图3a经横向、纵向2倍插值后所得复值图像的幅度；
图3d为本发明实施例1中，B超图像经插值后计算所得热应变图像序列的第130帧；
图3e为本发明实施例1中，热电偶所在空间点累积热应变随时间变化的曲线，比较了低采样率B超图像计算结果和B超图像纵向、横向均插值2倍后的计算结果；两者均对热电偶测得温度曲线进行最小方差拟合；
图4a为本发明实施例2中，低采样率B超图像计算所得热应变图像序列的第125帧；
图4b为本发明实施例2中，B超图像经插值后计算所得热应变图像序列的第125帧；
图4c为本发明实施例2中，热电偶所在空间点累积热应变随时间变化的曲线，比较了低采样率B超图像计算结果和B超图像纵向、横向均插值2.5倍后的计算结果；两者均对热电偶测得温度曲线进行最小方差拟合。
具体实施方式
下面结合附图和具体实施例对本发明进行详细说明。
本发明的基于低采样率B超图像计算热应变分布的方法，包括以下步骤：
(1)根据应用需求，确定B超探头的类型、波束控制方式、几何形状和阵元数目；
(2)对生物组织进行局部热消融，通过机械或电子的方式，将加热源与B超仪的成像区域调整至同一范围，确保B超成像区域覆盖加热区域；
(3)利用所选探头结合B超仪主机对治疗区域进行监控，采集并输出按时间排列的数字超声图像序列。设共有N帧按时间顺序排列的B超图像，每帧含有L条横向排列的超声线；每条超声线含有M个纵向分布的实数数据点；选取任意一帧B超图像；
(4)对任意一条超声线上的所有数据点，将每条超声线上的数据点按位置的奇偶依次分为两组，以奇数组为实部、偶数组为虚部构建复数信号；
(5)对所得复信号进行快速傅里叶变换，根据所得谱图获得数字截止频率ωc1，使0～ωc1范围内频谱能量不小于总频谱能量的80％；能量由所有谱线幅度的平方和求得；
(6)根据纵向采样率所需提升的有理数倍数I1，将I1表示为I1＝P1/Q1，其中P1和Q1均为正整数；设计一个低通数字滤波器，其通带截止频率为ωc1/P1，群延迟为正整数τ1；对(4)中所得复信号，在每个数据点后补充P1个0，在信号末尾再补充τ1个0；利用所得滤波器对该信号进行滤波，对结果信号去除其最前面τ1个数据点；
(7)对(6)中所得信号，选择序号为1～Q1的任意数据点为保留数据点，随后每隔Q1-1保留一个数据点，其余数据点均舍弃，得到(4)中超声线的纵向插值信号；
(8)对(3)中选中的图像，按(4)-(7)对所有L条纵向超声线进行处理，所得的L条插值信号线按原顺序进行组合，得到纵向插值后的一帧高采样率B超图像；
(9)对(3)中的所有N帧B超图像，按(4)-(8)的步骤对所有帧进行处理，得到包含N帧的纵向插值图像序列，每帧含有L条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点，M1为接近M×I1/2的正整数；选取其中任意一帧图像；
(10)对第k条横向超声信号按照从左到右依次用自然数从1开始标记序号，如第r个数据点则标记为r；对该横向信号进行快速傅里叶变换，根据谱图获得数字截止频率ωc2，使0～ωc2范围内频谱能量不小于总频谱能量的80％；
(11)设计一个低通数字滤波器，通带截止频率为ωc2，群延迟为正整数τ2；对(10)中选中的横向超声信号，在其末尾补充τ2个0后，利用所得滤波器对其进行滤波；将滤波所得信号最前面τ2个数据点舍弃；
(12)根据横向采样率所需提升的有理数倍数I2，将I2表示为I2＝P2/Q2，其中P2和Q2均为正整数；设计一个低通数字滤波器，其通带截止频率为ωc2/P2，群延迟为正整数τ3；对(11)中所得复信号，在每个数据点后补充P2个0，在信号末尾再补充τ3个0；利用所得滤波器对该信号进行滤波，对结果信号去除其最前面τ3个数据点；
(13)对(12)中所得信号，选择序号为1～Q2的任意数据点为保留数据点，随后每隔Q2-1保留一个数据点，其余数据点均舍弃；
(14)对(9)中选中的图像，按(10)-(13)对所有横向超声线进行处理后按原顺序排列，得到纵、横向均插值后的一帧高采样率B超图像；
(15)对(9)中的所有N帧B超图像，按(10)-(14)的步骤对所有帧进行处理，得到包含N帧的图像序列，每帧含有L1条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点，L1为接近L×I2的正整数；
(16)对(15)中获得的图像序列，选择任意一帧图像为参考帧，其后续一帧图像为目标帧；定义图像中的横向坐标为i，纵向坐标为j。i为整数，取值范围1～L1，j为整数，取值范围1～M1；选择奇数自然数L2、M2、L3、M3；
(17)在目标帧中选择一个点，其坐标为i＝l，j＝m，其中l为正整数；以该点为中心选择一个矩形区域，将该区域内的二维超声信号记为Uo(l,m)；
(18)在参考帧中选择一个点，其坐标为i＝l+p，j＝m+q；以该点为中心选择一个矩形区域，将该区域内的二维超声信号取共轭，记为Ur(l,m；p,q)；
(19)设计一个中心位于i＝i0，j＝j0、横向尺寸为L3、纵向尺寸为M3的二维幅度变迹窗函数。将窗函数中心(i0,j0)设置于(l,m)，使其与Uo(l,m)相乘，得到目标帧中以(l,m)为中心的二维信号Vo(l,m)；将窗函数的中心(i0,j0)设置于(l+p,m+p)，使其与Uo(l,m)相乘，得到参考帧中以(l+p,m+p)为中心的二维信号Vr(l,m；p,q)；
(20)对每一组(l,m)，改变p和q的取值，使参数组合(p,q)遍历二维取值空间，寻找使Vo(l,m)和Vr(l,m；p,q)相关性最大的(p,q)组合(plm,qlm)；对(l,m)参数组合进行变化，使其遍历二维空间，可获得plm和qlm的二维分布图象Dp(l,m)和Dq(l,m)；
(21)取纵向位移图像Dq(l,m)；其包含L1-L2+1条横向分布的信号线，每条信号线在纵向的长度为M1-M2+1；对其每一条信号线应用K点滑动平均；经处理后纵向位移图像的纵向长度减小为M1-M2-K+2；
(22)设计一个数字差分滤波器，群延迟为正整数τ4；对(21)中的任意一条信号线，在其末尾补充τ4个0，再利用所得差分滤波器对其滤波，将所得结果最前方τ4个数据点舍弃；针对每条信号线同样处理，可获得当前图像帧对应的组织热应变的二维分布图像；
(23)将(16)中选择的目标帧设定为新的参考帧，其后续一帧作为新的目标帧，重复(17)-(22)中的处理方法；依此类推，直至所有超声图像处理完毕，即可得到N-1个按时间排列的热应变二维分布图像。
作为上述方法的优选，所述步骤(1)中，可用超声探头类型包括凸阵探头、线阵探头、相控阵探头；可用波束控制方式包括线扫、相控阵、机械扇扫和面阵；可用几何形状包括弧形、圆形和矩形；可用阵元的数目包括单阵元或多阵元。
作为上述方法的优选，所述步骤(2)中，热消融方法包括微波消融、射频消融、超声消融、激光消融和红外线消融。
作为上述方法的优选，所述步骤(3)中，采集的超声图像序列中每帧图像由探头接收的原始超声回波信号经正交检波、降低采样率、波束合成后形成。
作为上述方法的优选，所述步骤(6)和(11)中，滤波器均具有低通、线性相位、有限长单位抽样响应、群延迟为正整数、通带最大衰减不大于3dB、阻带最小衰减不小于20dB、通带内能量不小于总能量80％的性质。
作为上述方法的优选，所述步骤(12)中，滤波器均具有低通、线性相位、有限长单位抽样响应、群延迟为正整数、通带最大衰减不大于3dB、阻带最小衰减不小于20dB的性质。
作为上述方法的优选，所述步骤(16)中，奇数自然数L2的取值范围为5～4L1/5，M2的取值范围为5～4M1/5，L3的取值范围为3～4L2/5，M3的取值范围为3～4M2/5。
作为上述方法的优选，所述步骤(17)中，l为正整数，取值范围为(L2+1)/2～L1-(L2-1)/2，m为正整数，取值范围为(M2+1)/2～M1-(M2-1)/2；所选矩形窗的横向范围为l-(L3-1)/2～l+(L3-1)/2，纵向范围为m-(M3-1)/2～m+(M3-1)/2。
作为上述方法的优选，所述步骤(18)中，p为整数，取值范围为-(L2-L3)/2～(L2-L3)/2，q为整数，取值范围为-(M2-M3)/2～(M2-M3)/2；所选矩形窗的横向范围为l+p-(L3-1)/2～l+p+(L3-1)/2，纵向范围为m+q-(M3-1)/2～m+q+(M3-1)/2。
作为上述方法的优选，所述步骤(19)中，二维幅度变迹函数可为如下形式：
其中，
作为上述方法的优选，，所述步骤(20)中，变量组合(p,q)的遍历空间为矩形区域[-(L2-L3)/2～(L2-L3)/2,-(M2-M3)/2～(M2-M3)/2]；参数组合(l,m)的遍历空间为矩形区域[(L2+1)/2～L1-(L2-1)/2,(M2+1)/2～M1-(M2-1)/2]。
作为上述方法的优选，，所述步骤(20)中，二维信号和的相关性按下式计算：
作为上述方法的优选，，所述步骤(21)中，K为任意不小于3的奇数。
作为上述方法的优选，所述步骤(22)中，数字差分器具有线性相位、有限长单位抽样响应、最大波纹不超过3dB的性质。
下面应用上述方法结合具体实例进行介绍。
实施例1
计算猪脂肪组织中的热应变，B超图像纵向、横向插值倍数均为2倍。具体步骤为：
步骤一、对一块猪脂肪组织进行微波加热，利用商用B超采集图像数据，将热电偶探针插入脂肪组织，使其前端正好位于超声成像平面上，如图2所示。微波加热功率1瓦，B超输出图像数据横向128个点，纵向512个点，帧率为10帧/秒，从微波加热开始时刻起，共采集13秒B超和热电偶数据。
步骤二、不对低分辨率B超图像进行任何插值处理，即在前述步骤(4)构建复信号完成后，直接跳转到步骤(16)进行热应变计算。第一帧构建完毕的B超幅度图像如图3a所示。此时每帧复信号B超图像大小为：L1＝128、M1＝256。取参数L2＝7、M2＝41、L3＝3、M3＝31；不应用二维幅度变迹窗函数；步骤(21)中，对位移数据进行滑动平均时选择K＝3。所得热应变图像的尺寸为横122点、纵213点。所得130帧热应变图像中第130帧如图3b所示。
步骤三、对低分辨率B超图像在纵向和横向均进行2倍插值处理。第一帧插值完毕的B超幅度图像如图3c所示。插值后每帧复信号B超图像大小为：L1＝256、M1＝512。取参数L2＝13、M2＝81、L3＝5、M3＝61；二维幅度变迹窗函数的形式选第一种；步骤(21)中，对位移数据进行滑动平均时选择K＝5。所得热应变图像的尺寸为横244点、纵427点。所得130帧热应变图像中第130帧如图3d所示。
步骤四、取步骤二所得每帧热应变图像中横纵坐标为(79，67)的数据点，该位置对应热电偶所在空间点，将总计130个热应变值按所在图像帧的顺序排列后，计算其累加序列，即得到该空间点对应的累积热应变序列，与热电偶测得的温度曲线具有相关性。为此，对累积热应变序列×C(C为拟合常数)与热电偶测得温度曲线进行最小方差拟合。根据拟合方差的大小，判断热应变计算的准确性。方差越小，则准确性越好。结果如图3e所示。
步骤五、取步骤三所得每帧热应变图像中横纵坐标为(158，134)的点，该位置也对应热电偶所在的空间点，因此将所有热应变图像中总计130个热应变值按所在图像帧的顺序后，计算其累加序列，即得到该空间点对应的累积热应变序列。对累积热应变序列×C(C为拟合常数)与热电偶测得温度曲线进行最小方差拟合。结果也呈现在图3e中。
由图3a和图3c比较可见，两幅图像的形态基本一致。这表明，本发明提供的插值方法仅提升了图像的采样率，并未改变图像本身的空间特征。
由图3b和图3d比较可见，本发明提供的热应变计算方法一方面提高了热应变图像的采样率，另一方面也使热应变集中区域的图像更加清晰。
由图3e中三条曲线比较可见，本发明提供的热应变计算方法所得结果相比于未经插值、无幅度变迹处理的方法，其结果更接近热电偶测量数据，误差更小，有效地改善了曲线呈现“台阶状”的问题，有利于提升超声测温的精度。
实施例2
计算猪肝组织中的热应变，B超图像纵向、横向插值倍数均为2.5倍。具体步骤为：
步骤一、对一块猪肝组织进行微波加热，利用商用B超采集图像数据，将热电偶探针插入猪肝组织，使其前端正好位于超声成像平面上，如图2所示。微波加热功率0.5瓦，B超输出图像数据横向128个点，纵向512个点，帧率为10帧/秒，从微波加热开始时刻起，共采集12.5秒B超和热电偶数据。
步骤二、不对低分辨率B超图像进行任何插值处理，即在前述步骤(4)构建复信号完成后，直接跳转到步骤(16)进行热应变计算。此时每帧复信号B超图像大小为：L1＝128、M1＝256。取参数L2＝7、M2＝41、L3＝3、M3＝31；不应用二维幅度变迹窗函数；步骤(21)中，对位移数据进行滑动平均时选择K＝3。所得热应变图像的尺寸为横向122点、纵向213点。所得125帧热应变图像中第125帧如图4a所示。
步骤三、对低分辨率B超图像在纵向和横向均进行2.5倍插值处理。插值后每帧复信号B超图像大小为：L1＝320、M1＝640。取参数L2＝17、M2＝101、L3＝9、M3＝81；二维幅度变迹窗函数的形式选第二种；步骤(21)中，对位移数据进行滑动平均时选择K＝7。所得热应变图像的尺寸为横向304点、纵向532点。所得125帧热应变图像中第125帧如图4b所示。
步骤四、取步骤二所得每帧热应变图像中横纵坐标为(72，88)的数据点，该位置对应热电偶所在空间点，将总计125个热应变值按所在图像帧的顺序排列得到一个序列，计算该序列的累加序列，即得到该空间点对应的累积热应变序列。该序列应与热电偶测得的温度曲线具有相关性。为此，对累积热应变序列×C(C为拟合常数)与热电偶测得温度曲线进行最小方差拟合。结果如图4c所示。
步骤五、取步骤三所得每帧热应变图像中横纵坐标为(180，220)的点，该位置也对应热电偶所在的空间点，将所有热应变图像中总计125个热应变值按所在帧的顺序排列得到一个序列后，计算其累加序列，即得到该空间点对应的累积热应变序列。对累积热应变序列×C(C为拟合常数)与热电偶测得温度曲线进行最小方差拟合。结果也呈现在图4c中。
由图4a和图4b比较可见，本发明提供的热应变计算方法一方面提高了热应变图像的采样率，另一方面也使热应变集中区域的图像更加清晰。
由图4c中三条曲线比较可见，本发明提供的热应变计算方法所得结果相比于未经插值、无幅度变迹处理的方法，其结果更接近热电偶测量数据，误差更小，有效地改善了曲线呈现“台阶状”的问题，有利于提升超声测温的精度。
以上显示和描述了本发明的基本原理、主要特征和优点。本领域的普通技术人员应该了解，上述实施例不以任何形式限制本发明的保护范围，凡采用等同替换等方式所获得的技术方案，均落于本发明的保护范围内。
本发明未涉及部分均与现有技术相同或可采用现有技术加以实现。
Claims (8)
1.一种基于低采样率B超声图像计算热应变分布的方法，其特征在于，该方法包括以下步骤：
步骤一、在生物组织热消融过程中采集超声图像信号，利用所选探头结合B超仪主机对治疗区域进行监控，采集并输出按时间排列的数字超声图像序列；
步骤二、对低采样率B超图像序列进行纵向插值，获得高采样率图像序列；纵向插值采取以下两种有效的信号插值方法的任一种，两种方法分别为“时域补零后低通滤波”法和“频域插零后逆傅里叶变换”法；其中，
所述“时域补零后低通滤波”法的具体步骤包括：
(1)设共有N帧按时间顺序排列的B超图像，每帧含有L条横向排列的超声线；每条超声线含有M个纵向分布的实数数据点；选取任意一帧B超图像；
(2)对任意一条超声线上的所有数据点，按照其距离B超探头的远近依次用自然数从1开始标记序号，第m个数据点则标记为m；将每条超声线上的数据点按标记值的奇偶依次分为两组，以奇数组为实部、偶数组为虚部构建复数信号；
(3)对所得复数信号进行快速傅里叶变换，根据所得谱图获得数字截止频率ωc1，使0～ωc1范围内频谱能量不小于总频谱能量；能量由所有谱线幅度的平方和求得；
(4)根据纵向采样率所需提升的有理数倍数I1，将I1表示为I1＝P1/Q1，其中P1和Q1均为正整数；设计一个具有线性相位性质的有限长单位抽样响应低通数字滤波器，其通带截止频率为ωc1/P1，群延迟为正整数τ1；
(5)对(2)中所得复数信号，在每个数据点后补充P1个0，在信号末尾再补充τ1个0；利用(4)中所得滤波器对该信号进行滤波，对结果信号去除其最前面τ1个数据点；
(6)对(5)中所得信号，选择序号为1～Q1的任意数据点为保留数据点，随后每隔Q1-1保留一个数据点，其余数据点均舍弃；至此，得到(2)中超声线的纵向插值信号；
(7)对(1)中选中的图像，按(2)-(6)对所有L条纵向超声线进行处理，所得的L条插值信号线按原顺序进行组合，得到纵向插值后的一帧高采样率B超图像；
(8)对(1)中的所有N帧B超图像，按(2)-(7)的步骤对所有帧进行处理，得到包含N帧的纵向插值图像序列，每帧含有L条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点，M1为接近M×I1/2的正整数；
所述“频域插零后逆傅里叶变换”法的具体步骤包括：
(Ⅰ)设共有N帧按时间顺序排列的B超图像，每帧含有L条横向排列的超声线；每条超声线含有M个纵向分布的实数数据点；选取任意一帧B超图像；
(Ⅱ)对任意一条超声线上的所有数据点，按照其距离B超探头的远近依次用自然数从1开始标记序号，第m个数据点则标记为m；将每条超声线上的数据点按标记值的奇偶依次分为两组，以奇数组为实部、偶数组为虚部构建复数信号；
(Ⅲ)对所得复数信号进行快速傅里叶变换，得到点数为T0的频谱序列，T0＝M/2，按照变换前复数信号的先后顺序，从1开始标记序号，第k个数据点则标记为k，k取值范围为1～T0，该点对应的值表示为X1(k)；
(Ⅳ)根据纵向采样率所需提升的有理数倍数I1，将I1表示为I1＝P1/Q1，其中P1和Q1均为正整数；
(Ⅴ)对(Ⅲ)中所得的频谱序列在频域插零，令长度为T1的空序列Y1(n)，其中T1＝T0*P1，n取值1～T1，Y1(n)表示对应n点的值；当T0为奇数，在n∈[1,(T0+1)/2]时，Y1(n)＝P1*X1(n)，在n∈[(T0+1)/2+1,T1-(T0-1)/2]，Y1(n)＝0，在n∈[T1-(T0-1)/2+1,T1]，Y1(n)＝P1*X1(n-T1+T0)；当T0为偶数，在n∈[1,T0/2]时，Y1(n)＝P1*X1(n)；n＝T0/2+1时，Y1(n)＝0.5*P1*X1(n)，在n∈[T0/2+2,T1-T0/2]，Y1(n)＝0；在n＝T1-T0/2+1时，Y1(n)＝0.5*P1*X1(T0/2+1)，在n∈[T1-T0/2+2,T1]，Y1(n)＝P1*X1(n-T1+T0)；
(Ⅵ)对(Ⅴ)中得到的频域插零后的Y1(n)序列，进行傅里叶逆变换，对所得信号，选择序号为1～Q1的任意数据点为保留数据点，随后每隔Q1-1保留一个数据点，其余数据点均舍弃；至此，得到(Ⅱ)中超声线的纵向插值信号；
(Ⅶ)对(Ⅰ)中选中的图像，按(Ⅱ)-(Ⅵ)对所有L条纵向超声线进行处理，所得的L条插值信号线按原顺序进行组合，得到纵向插值后的一帧高采样率B超图像；
(Ⅷ)对(Ⅰ)中的所有N帧B超图像，按(Ⅱ)-(Ⅶ)的步骤对所有帧进行处理，得到包含N帧的纵向插值图像序列，每帧含有L条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点，M1为接近M×I1/2的正整数；
步骤三、对纵向插值后的图像序列进行横向插值，使每帧具有更多的超声线；横向插值采取以下两种有效的信号插值方法的任一种，两种方法分别为“时域补零后低通滤波”法和“频域插零后逆傅里叶变换”法；其中，
所述“时域补零后低通滤波”法的具体步骤包括：
(1)共有N帧按时间顺序排列的纵向插值后B超图像，每帧含有L条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点；选取任意一帧B超图像；
(2)将每条超声线的第k个复数数据点按照超声线的排列顺序进行组合，得到第k条横向超声信号；每帧图像共包含M1条横向超声信号，每条长度为L；
(3)对第k条横向超声信号按照从左到右依次用自然数从1开始标记序号，第r个数据点则标记为r；对该横向信号进行快速傅里叶变换，根据谱图获得数字截止频率ωc2，使0～ωc2范围内频谱能量不小于总频谱能量；
(4)设计一个具有线性相位性质的有限长单位抽样响应一次低通数字滤波器，通带截止频率为ωc2，群延迟为正整数τ2；
(5)对(3)中选中的横向超声信号，在其末尾补充τ2个0后，利用(4)中所得滤波器对其进行滤波；将滤波所得信号最前面τ2个数据点舍弃；得到第k条防混叠横向超声信号；
(6)根据横向采样率所需提升的有理数倍数I2，将I2表示为I2＝P2/Q2，其中P2和Q2均为正整数；设计一个具有线性相位性质的有限长单位抽样响应二次低通数字滤波器，其通带截止频率为ωc2/P2，群延迟为正整数τ3；
(7)对(5)中所得复数信号，在每个数据点后补充P2个0，在信号末尾再补充τ3个0；利用(6)中所得滤波器对该信号进行滤波，对结果信号去除其最前面τ3个数据点；
(8)对(7)中所得信号，选择序号为1～Q2的任意数据点为保留数据点，随后每隔Q2-1保留一个数据点，其余数据点均舍弃；至此得到(3)中横向超声信号的插值结果；
(9)对(1)中选中的图像，按(2)-(8)对所有M1条横向超声线进行处理，所得的M1条插值信号线按原顺序排列，得到纵、横向均插值后的一帧高采样率B超图像；
(10)对(1)中的所有N帧B超图像，按(2)-(9)的步骤对所有帧进行处理，得到包含N帧的图像序列，每帧含有L1条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点，L1为接近L×I2的正整数；
所述“频域插零后逆傅里叶变换”法的具体步骤包括：
(Ⅰ)共有N帧按时间顺序排列的纵向插值后B超图像，每帧含有L条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点；选取任意一帧B超图像；
(Ⅱ)将每条超声线的第k个复数数据点按照超声线的排列顺序进行组合，得到第k条横向超声信号；每帧图像共包含M1条横向超声信号，每条长度为L；
(Ⅲ)对第k条横向超声信号按照从左到右依次用自然数从1开始标记序号，第r个数据点则标记为r；对该横向信号进行快速傅里叶变换，根据谱图获得数字截止频率ωc2，使0～ωc2范围内频谱能量不小于总频谱能量；
(Ⅳ)设计一个具有线性相位性质的有限长单位抽样响应一次低通数字滤波器，通带截止频率为ωc2，群延迟为正整数τ2；
(Ⅴ)对(Ⅲ)中选中的横向超声信号，在其末尾补充τ2个0后，利用(4)中所得滤波器对其进行滤波；将滤波所得信号最前面τ2个数据点舍弃；得到第k条防混叠横向超声信号；
(Ⅵ)对(Ⅴ)所得横向超声信号进行快速傅里叶变换，得到点数为L的频谱序列，按照变换前复数信号的顺序，从1开始标记序号，第k个数据点则标记为k，k取值范围为1～L，该点对应的数值表示为X2(k)；
(Ⅶ)根据横向采样率所需提升的有理数倍数I2，将I2表示为I2＝P2/Q2，其中P2和Q2均为正整数；
(Ⅷ)对(Ⅵ)中所得的频谱序列在频域插零；令长度为T2的空序列Y2(n)，其中T2＝L*P2，n取值1～T2，Y2(n)表示对应n点的值；当L为奇数，在n∈[1,(L+1)/2]时，Y2(n)＝P2*X2(n)，在n∈[(L+1)/2+1,T2-(L-1)/2]，Y2(n)＝0，在n∈[T2-(L-1)/2+1,T2]，Y2(n)＝P2*X2(n-T2+L)；当L为偶数，在n∈[1,L/2]时，Y2(n)＝P2*X2(n)，n＝L/2+1时，Y2(n)＝0.5*P2*X2(n)，在n∈
[L/2+2,T2-L/2]，Y2(n)＝0；在n＝T2-L/2+1时，Y2(n)＝0.5*P2*X2(L/2+1)，在n∈[T2-L/2+2,T2]，Y2(n)＝P2*X2(n-T2+L)；
(Ⅸ)对(Ⅷ)中得到的频域插零后的Y2(n)序列，进行傅里叶逆变换，对所得信号，选择序号为1～Q2的任意数据点为保留数据点，随后每隔Q2-1保留一个数据点，其余数据点均舍弃；至此，得到(Ⅱ)中超声线的横向插值信号；
(Ⅹ)对(Ⅰ)中选中的图像，按(Ⅱ)-(Ⅸ)对所有M1条横向超声线进行处理，所得的M1条插值信号线按原顺序排列，得到纵、横向均插值后的一帧高采样率B超图像；
(Ⅺ)对(Ⅰ)中的所有N帧B超图像，按(Ⅱ)-(Ⅹ)的步骤对所有帧进行处理，得到包含N帧的图像序列，每帧含有L1条横向排列的超声线；每条超声线含有M1个纵向分布的复数数据点，L1为接近L×I2的正整数；
步骤四、根据插值后的B超图像序列，计算热应变实时分布图像。
2.根据权利要求1所述的一种基于低采样率B超声图像计算热应变分布的方法，其特征在于，所述步骤一的具体步骤包括：
(1)根据应用需求，确定B超探头的类型、波束控制方式、几何形状和阵元数目；
(2)对生物组织进行局部热消融；
(3)利用所选探头结合B超仪主机对治疗区域进行监控，采集并输出按时间排列的数字超声图像序列；该序列中每帧图像由探头接收的原始超声回波信号经正交检波、降低采样率、波束合成后形成。
3.根据权利要求1或2所述的一种基于低分辨率B超图像计算热应变分布的方法，其特征在于，所述B超仪主机所使用的超声探头的型号与工作频率根据应用需求选用，超声探头类型包括凸阵探头、线阵探头、相控阵探头；波束控制方式包括线扫、相控阵、机械扇扫和面阵；几何形状包括弧形、圆形和矩形；阵元的数目为单阵元或多阵元。
4.根据权利要求2所述的一种基于低分辨率B超图像计算热应变分布的方法，其特征在于，所述对生物组织进行局部热消融的方法包括微波消融、射频消融、超声消融、激光消融和红外线消融。
5.根据权利要求1所述的一种基于低采样率B超声图像计算热应变分布的方法，其特征在于，在步骤二的“时域补零后低通滤波”插值法中，纵向插值所用低通滤波器具有线性相位、有限长单位抽样响应、通带最大衰减不超过3dB、阻带最小衰减不小于20dB、截止频率内信号能量不低于总能量80％、群延迟为正整数的特点。
6.根据权利要求1所述的一种基于低采样率B超声图像计算热应变分布的方法，其特征在于，在步骤三的“时域补零后低通滤波”插值法中，横向插值所用一次低通滤波器具有线性相位、有限长单位抽样响应、通带最大衰减不超过3dB、阻带最小衰减不小于20dB、截止频率内信号能量不低于总能量80％、群延迟为正整数的特点；二次低通滤波器具有线性相位、有限长单位抽样响应、通带最大衰减不超过3dB、阻带最小衰减不小于20dB、群延迟为正整数的特点。
7.根据权利要求1所述的一种基于低分辨率B超图像计算热应变分布的方法，其特征在于，所述步骤四中计算热应变图像时，由每条纵向信号线获得热应变信号通过滑动平均、提前补零、数字差分滤波、群延迟补偿四个步骤进行。
8.根据权利要求1所述的一种基于低采样率B超声图像计算热应变分布的方法，其特征在于，所述步骤四的具体步骤包括：
(1)选择任意一帧图像为参考帧，其后续一帧图像为目标帧；定义图像中的横向坐标为i，纵向坐标为j，i为整数，取值范围1～L1，j为整数，取值范围1～M1；选择L2为5～4L1/5之间的任意奇数自然数，M2为5～4M1/5之间的任意奇数自然数，L3为3～4L2/5之间的任意奇数自然数，M3为3～4M2/5之间的任意奇数自然数；其中L1为B超图像经插值后、横向分布的超声信号线的条数，M1为B超图像经插值后，纵向分布的超声信号线的数据点数；
(2)在目标帧中选择一个点，其坐标为i＝l，j＝m，其中l为正整数，其取值范围为(L2+1)/2～L1-(L2-1)/2，m为正整数，其取值范围为(M2+1)/2～M1-(M2-1)/2；以该点为中心选择i＝l-(L3-1)/2～l+(L3-1)/2、j＝m-(M3-1)/2～m+(M3-1)/2的矩形区域，将该区域内的二维超声信号记为Uo(l,m)；
(3)在参考帧中选择一个点，其坐标为i＝l+p，j＝m+q，其中p为整数，其取值范围为-(L2-L3)/2～(L2-L3)/2，q为整数，其取值范围为-(M2-M3)/2～(M2-M3)/2；以该点为中心选择i＝l+p-(L3-1)/2～l+p+(L3-1)/2、j＝m+q-(M3-1)/2～m+q+(M3-1)/2的矩形区域，将该区域内的二维超声信号取共轭，记为Ur(l,m；p,q)；
(4)设计一个中心位于i＝i0，j＝j0、横向尺寸为L3、纵向尺寸为M3的二维幅度变迹窗函数，可选用如下设计形式：
其中，
(5)将(4)中窗函数的中心(i0,j0)设置于(l,m)，使其与二维信号Uo(l,m)相乘，得到目标帧中以(l,m)为中心的二维变迹加权信号Vo(l,m)；将(4)中窗函数的中心(i0,j0)设置于(l+p,m+p)，使其与二维信号Uo(l,m)相乘，得到参考帧中以(l+p,m+p)为中心的二维变迹加权信号Vr(l,m；p,q)；
(6)对(2)中选定的每一个l和m的具体值，改变p和q的取值，使参数组合(p,q)遍历二维取值空间[-(L2-L3)/2～(L2-L3)/2,-(M2-M3)/2～(M2-M3)/2]，寻找使如下相关函数取值最大的(p,q)组合(plm,qlm)：
(7)plm和qlm即分别为位于i＝l、j＝m的组织由于热膨胀导致的横、纵向位移；对(2)中选中的(l,m)参数组合进行变化，使其遍历二维空间
[(L2+1)/2～L1-(L2-1)/2,(M2+1)/2～M1-(M2-1)/2]，可获得plm和qlm的二维分布图象Dp(l,m)和Dq(l,m)；
(8)取(7)中获得的纵向位移图像Dq(l,m)；其包含L1-L2+1条横向分布的信号线，每条信号线在纵向的长度为M1-M2+1；对其每一条信号线应用K点滑动平均，K为任意不小于3的奇数，经处理后纵向位移图像的纵向长度减小为M1-M2-K+2；
(9)设计一个具有线性相位性质的有限长单位抽样响应数字差分滤波器，其最大波纹不超过3dB，群延迟为正整数τ4；对(8)中的任意一条信号线，在其末尾补充τ4个0，再利用所得差分滤波器对其滤波，将所得结果最前方τ4个数据点舍弃；至此，即得到该条信号线对应的热应变信号；
(10)对(8)中的纵向位移图像，针对每条信号线均按(9)进行处理，可获得当前图像帧对应的组织热应变的二维分布图像；
(11)将(1)中选择的目标帧设定为新的参考帧，其后续一帧作为新的目标帧，重复(2)-(10)中的处理方法；依此类推，直至所有超声图像处理完毕，即可得到N-1个按时间排列的热应变二维分布图像。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910112033.8A CN109615677B (zh) | 2019-02-13 | 2019-02-13 | 一种基于低采样率b超图像计算热应变分布的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910112033.8A CN109615677B (zh) | 2019-02-13 | 2019-02-13 | 一种基于低采样率b超图像计算热应变分布的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109615677A CN109615677A (zh) | 2019-04-12 |
| CN109615677B true CN109615677B (zh) | 2023-05-12 |
Family
ID=66018870
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910112033.8A Active CN109615677B (zh) | 2019-02-13 | 2019-02-13 | 一种基于低采样率b超图像计算热应变分布的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109615677B (zh) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104188689A (zh) * | 2013-10-25 | 2014-12-10 | 中国科学院深圳先进技术研究院 | 基于超声回波射频信号的组织位移估算方法和系统 |
| CN107569256A (zh) * | 2017-09-25 | 2018-01-12 | 南京广慈医疗科技有限公司 | 基于热膨胀和门控算法测量生物组织温度变化的超声方法 |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8403850B2 (en) * | 2008-03-25 | 2013-03-26 | Wisconsin Alumni Research Foundation | Rapid two/three-dimensional sector strain imaging |
| US9610063B2 (en) * | 2010-03-26 | 2017-04-04 | The Johns Hopkins University | Methods and apparatus for ultrasound strain imaging |
| US9320491B2 (en) * | 2011-04-18 | 2016-04-26 | The Trustees Of Columbia University In The City Of New York | Ultrasound devices methods and systems |
| CN102961166A (zh) * | 2011-08-31 | 2013-03-13 | 通用电气公司 | 用于检测和跟踪针的方法 |
| CN108784736B (zh) * | 2018-05-23 | 2020-02-14 | 成都信息工程大学 | 一种二维迭代的超声弹性成像应变估计方法 |
| CN109188409A (zh) * | 2018-10-24 | 2019-01-11 | 重庆大学 | 一种基于Chirp码的正交稀疏字典设计方法 |
-
2019
- 2019-02-13 CN CN201910112033.8A patent/CN109615677B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104188689A (zh) * | 2013-10-25 | 2014-12-10 | 中国科学院深圳先进技术研究院 | 基于超声回波射频信号的组织位移估算方法和系统 |
| CN107569256A (zh) * | 2017-09-25 | 2018-01-12 | 南京广慈医疗科技有限公司 | 基于热膨胀和门控算法测量生物组织温度变化的超声方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109615677A (zh) | 2019-04-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Porter et al. | Three-dimensional registration and fusion of ultrasound and MRI using major vessels as fiducial markers | |
| Pesavento et al. | New real-time strain imaging concepts using diagnostic ultrasound | |
| EP1757955A1 (en) | Apparatus and method for processing an ultrasound image | |
| CN107569256B (zh) | 基于热膨胀和门控算法测量生物组织温度变化的超声方法 | |
| US20160249882A1 (en) | MRI Compatible 3-D Intracardiac Echography Catheter and System | |
| EP2148622A2 (en) | Improved imaging system | |
| CN106821500B (zh) | 一种用于微创手术导航系统 | |
| CN105520752B (zh) | 波束成形设备、超声成像设备以及波束成形方法 | |
| Mirzaei et al. | 3D normalized cross-correlation for estimation of the displacement field in ultrasound elastography | |
| WO2015124069A1 (zh) | 基于rf数据超声成像处理方法及系统 | |
| JP3833597B2 (ja) | 超音波撮像装置及び超音波撮像方法 | |
| CN109615677B (zh) | 一种基于低采样率b超图像计算热应变分布的方法 | |
| Pesavento et al. | System for real-time elastography | |
| Ingle et al. | Three-dimensional sheaf of ultrasound planes reconstruction (SOUPR) of ablated volumes | |
| CN111272305B (zh) | 一种基于非线性热膨胀评估温度变化的超声方法及系统 | |
| Rivaz et al. | Tracked regularized ultrasound elastography for targeting breast radiotherapy | |
| Krueger et al. | Ultrasonic strain imaging of the female breast using phase root seeking and three-dimensional" optical flow" | |
| EP1623245A2 (en) | Mri probe design and tracking, and efficient mri reconstruction and deblurring | |
| Rohling et al. | Issues in 3-D free-hand medical ultrasound imaging | |
| Grube et al. | Ultrasound shear wave velocity estimation in a small field of view via spatio-temporal deep learning | |
| Yin et al. | Thermal strain imaging in vivo using interpolated IQ-images | |
| Lorenz et al. | Three-dimensional strain imaging and related strain artifacts using an ultrasonic 3D abdominal probe | |
| CN110931130B (zh) | 一种基于b超信号评估呼吸和心动周期的方法 | |
| Tierney et al. | Perfusion imaging with non-contrast ultrasound | |
| US20240058580A1 (en) | Magnetic resonance imaging guided active injection needle for radiation-oncology brachytherapy |
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 | ||
| CB03 | Change of inventor or designer information | ||
| CB03 | Change of inventor or designer information |
Inventor after: Kong Xiangqing Inventor after: Guo Xiasheng Inventor after: Yin Chuhao Inventor after: Zhang Dong Inventor after: Xue Honghui Inventor before: Guo Xiasheng Inventor before: Yin Chuhao Inventor before: Zhang Dong Inventor before: Kong Xiangqing Inventor before: Xue Honghui |
|
| GR01 | Patent grant | ||
| GR01 | Patent grant |