CN113109825A - 基于Radon变换的长骨相控超声信号定征与骨质评价系统 - Google Patents

基于Radon变换的长骨相控超声信号定征与骨质评价系统 Download PDF

Info

Publication number
CN113109825A
CN113109825A CN202110366930.9A CN202110366930A CN113109825A CN 113109825 A CN113109825 A CN 113109825A CN 202110366930 A CN202110366930 A CN 202110366930A CN 113109825 A CN113109825 A CN 113109825A
Authority
CN
China
Prior art keywords
wave
signal
fas
long bone
bone
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.)
Granted
Application number
CN202110366930.9A
Other languages
English (en)
Other versions
CN113109825B (zh
Inventor
许凯亮
凌非尧
陈琦
徐冬梅
杨春山
郑书航
他得安
李大为
王威琪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Shengai Medical Technology Co.,Ltd.
Original Assignee
Shanghai A & S Science Technology Development Co ltd
Fudan University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Shanghai A & S Science Technology Development Co ltd, Fudan University filed Critical Shanghai A & S Science Technology Development Co ltd
Priority to CN202110366930.9A priority Critical patent/CN113109825B/zh
Publication of CN113109825A publication Critical patent/CN113109825A/zh
Application granted granted Critical
Publication of CN113109825B publication Critical patent/CN113109825B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/88—Sonar systems specially adapted for specific applications
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/08—Detecting organic movements or changes, e.g. tumours, cysts, swellings
    • A61B8/0875—Detecting organic movements or changes, e.g. tumours, cysts, swellings for diagnosis of bone
    • 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/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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Biophysics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Veterinary Medicine (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Computer Hardware Design (AREA)
  • Rheumatology (AREA)
  • Orthopedic Medicine & Surgery (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Acoustics & Sound (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Computer Networks & Wireless Communication (AREA)

Abstract

本发明提供了一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，用于对待测长骨进行定征与评价，其特征在于，包括：信号采集模块利用超声探头采集到待测长骨对应的原始信号；降噪模块利用Radon正反变换对对原始信号进行处理得到降噪处理后信号；慢度截距图生成模块对降噪处理后信号进行Radon反变换得到慢度‑截距图；区域截取模块利用第一到达波与基阶反对称模态Lamb波的特征从慢度‑截距图中截取出FAS波区域与A0波区域；速度计算模块基于FAS波区域与A0波区域计算得到FAS波速度与A0波速度；长骨反演模块基于FAS波速度以及A0波速度对待测长骨进行反演得到长骨几何参数以及弹性参数从而进行评价得到评价结果。

Description

基于Radon变换的长骨相控超声信号定征与骨质评价系统
技术领域
本发明属于数据识别领域，具体涉及一种基于Radon变换的长骨相控超声信号定征与骨质评价系统。
背景技术
随着人口老龄化的加剧，骨质疏松症发病率仍有上升趋势。骨质状况的诊断和检测技术对于骨质疏松症防治具有重要意义。然而，目前我国骨质疏松症诊疗率仍有较大的城乡与地区差异，诊治率较低。原发性骨质疏松症诊疗指南(2017版)要点解读指出，当前我国骨质疏松症诊断率仅为2/3，接受有效治疗的不足1/4。因此，骨质状况评价方法的研究对于骨质疏松症的诊断与防治具有重要意义。
目前，临床上主要采用双能X射线技术(DXA)评价骨质状况。该技术可用于测量骨矿化密度和骨几何形态，而无法直接反映骨的力学性能。此外，DXA技术还有许多其他的缺点，如设备昂贵，体积笨重，有较强的电离辐射等，难以用于社区普查以及基层医疗。另一方面，超声是一种弹性波，介质的弹性系数和密度决定了声波的传播速度。因此，超声可用于骨质弹性评价。另外，超声设备体积小、价格低廉，无电离辐射，非常适于社区普查与基层医疗，具有很大的应用潜力。
目前，超声已被用于长骨皮质骨评价。长骨超声定征中常用的技术被称为轴向传播法，即通过测量沿着长骨长轴方向传播的超声波实现长骨皮质骨评价。在轴向传播的过程中，边界上的横波和纵波发生模式转换，最终会形成超声导波。因此，超声导波中携带着长骨中的弹性特征以及几何信息。在固体平板中，这种超声导波被称为Lamb波。Lamb波模型与动物长骨样本上测量的导波速度十分接近，因此轴向传播信号中的导波可以用Lamb波理论分析，进而为长骨的诊断提供帮助。
在轴向传输模式中，第一到达波(FAS)是非常重要的测量对象。FAS是长骨轴向传输信号中传播速度最快的信号成分，可以视作长骨中的纵波。目前已经发现，不同波长的FAS的速度可以反映骨的不同特性，如骨矿质化程度等。如果可以将FAS从多模式波中提取出来，就可以根据其速度反求出骨的弹性模量，从而对骨质状况做出评估。另外，基阶反对称模态Lamb波(A0波)是最慢的导波信号成分，是Lamb波理论中的最低阶反对称模式的解。A0波的速度可以较准确地反映长骨的皮质骨厚度以及弹性模量。然而，长骨中的Lamb波理论上有无穷多种模式的解，一种或者多种模式都可以用于长骨定征，如何从超声导波中准确地分离FAS成分以及A0波是一个难题。
发明内容
为解决上述问题，提供一种能够从超声导波中准确地分离FAS成分以及A0波的长骨相控超声信号定征与骨质评价系统，本发明采用了如下技术方案：
本发明提供了一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，用于对待测长骨进行定征与评价，其特征在于，包括：信号采集模块，利用超声探头从待测长骨中采集到原始信号；降噪模块，利用Radon正反变换对对原始信号进行处理得到降噪处理后信号；慢度截距图生成模块，对降噪处理后信号进行Radon反变换，得到Radon域上的慢度-截距图；区域截取模块，利用第一到达波的特征从慢度-截距图中截取出FAS波区域，利用基阶反对称模态Lamb波的特征从慢度-截距图中截取出A0波区域；速度计算模块，基于FAS波区域计算得到第一到达的传播速度作为FAS波速度，基于A0波区域计算得到基阶反对称模态Lamb波的传播速度作为A0波速度，并基于FAS波区域以及A0波区域计算得到实际频散曲线；长骨反演模块，基于FAS波速度以及A0波速度对待测长骨进行反演，从而得到该待测长骨对应的长骨几何参数以及弹性参数；以及骨质评价模块，基于长骨几何参数以及弹性参数，对待测长骨的骨质状况进行评价从而得到评价结果。
根据本发明提供的一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，还可以具有这样的技术特征，其中，超声探头包括探头壳体、至少一个发射阵元以及多个接收阵元，发射阵元与接收阵元分别设置在探头壳体的两端，发射阵元用于对待测长骨发射超声信号，接收阵元用于接收待测长骨因超声信号产生的信号并作为原始信号。
根据本发明提供的一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，还可以具有这样的技术特征，其中，FAS波速度通过如下步骤计算得到：步骤S1-1，对FAS波区域进行Radon正变换得到第一到达波在距离-时间域上的投影作为FAS波投影；步骤S1-2，对FAS波投影进行Hilbert变换得到第一到达波的包络信号作为第一包络信号；步骤S1-3，根据接收阵元与发射阵元之间的距离以及第一包络信号到达接收阵元的时间计算得到的速度作为FAS波速度，A0波速度通过如下步骤计算得到：步骤S2-1，对A0波区域进行Radon正变换得到基阶反对称模态Lamb波在距离-时间域上的投影作为A0波投影；步骤S2-2，对A0波投影进行Hilbert变换得到基阶反对称模态Lamb波的包络信号作为第二包络信号；步骤S2-3，根据接收阵元与发射阵元之间的距离以及第二包络信号到达接收阵元的时间计算得到的速度作为A0波速度。
根据本发明提供的一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，还可以具有这样的技术特征，其中，第一包络信号到达接收阵元的时间为第一包络信号达到最大值的时间，第二包络信号到达接收阵元的时间为第二包络信号达到最大值的时间。
根据本发明提供的一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，还可以具有这样的技术特征，其中，第一包络信号到达接收阵元的时间为第一包络信号达到最大值后幅度降为1/2的时间，第二包络信号到达接收阵元的时间为第二包络信号达到最大值后幅度降为1/2的时间。
根据本发明提供的一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，还可以具有这样的技术特征，其中，第一包络信号到达接收阵元的时间为第一包络信号达到最大值后幅度降为0的时间，第二包络信号到达接收阵元的时间为第二包络信号达到最大值后幅度降为0的时间。
根据本发明提供的一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，还可以具有这样的技术特征，其中，FAS波速度的计算方法为选取FAS波区域中幅度最大的点作为第一幅度最大点，将该第一幅度最大点对应的慢度值的倒数作为FAS波速度，A0波速度的计算方法为选取A0波区域中幅度最大的点作为第二幅度最大点，将该第二幅度最大点对应的慢度值的倒数作为A0波速度。
根据本发明提供的一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，还可以具有这样的技术特征，其中，FAS波速度的计算方法为选取FAS波区域中截距值最小的点作为第一截距最小点，将该第一截距最小点对应的慢度值的倒数作为FAS波速度，A0波速度的计算方法为选取A0波区域中截距值最小的点作为第二截距最小点，将该第二截距最小点对应的慢度值的倒数作为A0波速度。
根据本发明提供的一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，还可以具有这样的技术特征，其中，长骨反演模块包括：代价函数设计单元，根据实际频散曲线以及待测长骨的理论频散曲线涉及代价函数C：
Figure BDA0003007491300000051
式中，(i,j)为实际频散曲线上一点的坐标，(p,q)为理论频散曲线上一点的坐标，A(i,j)为实际频散曲线上点(i,j)的幅度；反演模型构建单元，选取参数范围作为预估范围，根据该预估范围构建反演模型；纵波速度计算单元，利用反演模型计算FAS波速度在不同厚度下的纵波速度；以及参数选定单元，根据代价函数C以及纵波速度计算得到不同参数下实际频散曲线与理论频散曲线之间的差值，将差值最小对应的理论频散曲线的参数作为目标参数，从目标参数中获取长骨几何参数以及弹性参数。
根据本发明提供的一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，还可以具有这样的技术特征，其中，降噪模块包括Radon反变换单元以及Radon正变换单元，Radon反变换单元对原始信号进行Radon反变换得到反变换后信号，Radon正变换单元对反变换后信号进行Radon正变换得到正变换后信号，并作为降噪处理后信号。
发明作用与效果
根据本发明的一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，由于降噪模块利用Radon正反变换对对原始信号进行处理得到降噪处理后信号，因此降噪处理后信号具有较高的信噪比。另外，由于区域截取模块基于第一到达波与基阶反对称模态Lamb波的特征从慢度-截距图中截取出FAS波区域与A0波区域，因此解决了傅里叶变换不能很好地提取出FAS和A0波信号的问题。除此之外，还由于速度计算模块基于FAS波区域与A0波区域计算得到FAS波速度与A0波速度，因此，便于后续根据FAS波速度与A0波速度精确反演得到待测长骨的几何参数以及弹性参数，更进一步地，提高了基于几何参数以及弹性参数的骨质评估结果的准确性。
附图说明
图1为本发明实施例的基于Radon变换的长骨相控超声信号定征与骨质评价系统的结构框图；
图2为本发明实施例的Radon变换对处理过程示意图；
图3为本发明实施例的降噪处理后信号处理过程示意图；
图4为本发明实施例的慢度-截距图中FAS截取过程示意图；
图5为本发明实施例的慢度-截距图中A0波截取过程示意图；
图6为本发明实施例的包络信号示意图；
图7为本发明实施例的反演结构示意图；
图8为本发明实施例的基于Radon变换的长骨相控超声信号定征与骨质评价系统工作过程的流程图。
具体实施方式
为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，以下结合实施例及附图对本发明的一种基于Radon变换的长骨相控超声信号定征与骨质评价系统作具体阐述。
<实施例>
本实例中，基于Radon变换的长骨相控超声信号定征与骨质评价系统对待测长骨进行定征与评价，其中，待测长骨选用的是1.8mm厚的骨板仿体，其横波速度为1.8km/s，纵波速度为4.0km/s。
图1为本发明实施例的基于Radon变换的长骨相控超声信号定征与骨质评价系统的结构框图。
如图1所示，基于Radon变换的长骨相控超声信号定征与骨质评价系统1包括信号采集模块11、降噪模块12、慢度截距图生成模块13、区域截取模块14、速度计算模块15、长骨反演模块16以及骨质评价模块17。
信号采集模块11利用超声探头从待测长骨中采集到原始信号。
其中，超声探头包括探头壳体、至少一个发射阵元以及多个接收阵元。
发射阵元与接收阵元分别设置在探头壳体的两端，从而有效地减小沿探头传播的超声信号对采集结果产生的影响。
发射阵元通过预定的宽带或窄带、编码或者非编码信号对待测长骨发射超声信号。发射阵元在发射超声信号时，每个发射阵元均为有效阵元，从而加大激励的强度，提高信噪比。
接收阵元用于接收待测长骨因超声信号产生的信号并作为原始信号。接收阵元在接收原始信号时，每个接收阵元均为有效阵元，从而接收沿待测长骨传播的信号。
本实施例中，超声探头型号为1MHz，探头上的总阵元总数为128个，其中，5个阵元为发射阵元，编号为1-5号阵元，设置在超声探头壳体的一端；64个阵元为接收阵元，编号为65-128号阵元。
在发射超声波时，5个激励阵元同时发射超声信号，其中心频率为1MHz。激励信号发射完成后，立即开始采集信号。
在接收超声波时，64个接受阵元均为有效接收阵元，从而接收到沿样本轴向传输的超声导波信号。
降噪模块12利用Radon变换对原始信号进行处理得到降噪处理后信号。
图2为本发明实施例的Radon变换对处理过程示意图。
其中，降噪模块12包括Radon反变换单元以及Radon正变换单元。
图2a为原始信号的距离-时间域图，图2b为反变换后信号的慢度-截距图，图2c为反变换后信号的距离-时间域图。
Radon反变换单元对原始信号(如图2中a部分所示)进行Radon反变换得到反变换后信号(如图2中b部分所示)。
从图2中a部分可以看出，由于环境噪声的干扰，导波信号与噪声信号相互混叠，因此原始信号的信噪比较低。同时，超声信号在骨板仿体上传播时还有多种模式波，因此原始信号是多种模式波混杂的结果。
其中，反变换后信号为原始信号在Radon域上的投影，从图2的b部分可以看出不同信号成分的截距与慢度都各有不同。
Radon正变换单元对反变换后信号进行Radon正变换得到正变换后信号，并作为降噪处理后信号(如图2中c部分所示)。
对比图2中a部分与图2中c部分，可以发现图2中c部分中的降噪处理后信号的信噪比有较大提升，可以更加清晰地分辨出沿骨板仿体传导的信号。
慢度截距图生成模块13对降噪处理后信号进行Radon反变换，得到Radon域上的慢度-截距图。
图3为本发明实施例的降噪处理后信号处理过程示意图。
图3中a-c分别为所有接收阵元对应的降噪处理后信号距离-时间域图、降噪处理后信号经过傅里叶变换后的频散曲线图以及降噪处理后信号经过Radon反变换后的变换结果对应的慢度-截距图。
图3中d-f分别为所有接收阵元对应的FAS波距离-时间域图、FAS波经过傅里叶变换后的频散曲线图以及FAS波经过Radon反变换后的变换结果对应的慢度-截距图。
图3中g-i分别为所有接收阵元对应的A0波距离-时间域图、A0波经过傅里叶变换后的频散曲线图以及A0波经过Radon反变换后的变换结果对应的慢度-截距图。
其中，图3g横坐标为距离(mm)，纵坐标为时间(μs)；图3h横坐标为频率(kHz)，纵坐标为波数(rad/mm)；图3i横坐标为慢度(s/km)，纵坐标为时间(μs)。
对比图3中a、d以及g，降噪处理后信号为一条较宽的带状信号，而FAS是整个降噪处理后信号中最上方的信号，A0波则混杂在带状信号中。
图3b与图3h中的细线为理论计算所得的A0波频散曲线。另外，从图3b，图3e以及图3h的对比结果可以看出，从频率域难以分离出不同信号成分。
区域截取模块14利用第一到达波的特征从慢度-截距图中截取出FAS波区域，利用基阶反对称模态Lamb波的特征从慢度-截距图中截取出A0波区域。
图4为本发明实施例的慢度-截距图中FAS截取过程示意图。
本实施例中，FAS的特征为沿骨板仿体传播的FAS大约为3.0-4.5km/s，因此，选择截取的FAS波区域为慢度范围0.22-0.33s/km、截距范围为12.0-15.5μs的区域(如图4中a部分圆圈部分所示)。
图4a为慢度-截距图，横坐标为慢度(s/km)，纵坐标为时间(μs)，从图4中a部分圆圈部分(即FAS波区域)可以看出，FAS位于整个慢度-截距图中各个信号区域的最上端，其余部分是超声信号传播过程中的其他模式导波。
为了显示FAS波区域截取前后的具体变化，本实施例中选取65、70、75、80以及85号接收阵元各自对应的降噪处理后信号进行FAS波区域截取前后的具体分析。
如图4中b部分所示，5根实线分别为65、70、75、80以及85号接收阵元对应的降噪处理后信号，虚线为FAS波，FAS波到达接收阵元的时间最短，对应的传播速度最快，但是信号幅度较小(即被虚线截断的靠前部分抖动幅度小)。
图4中c部分与d部分分别为其他模式导波(即FAS以外信号)对应的距离-时间域图以及FAS波区域(即FAS)对应的距离-时间域图，对比图5的b、c以及d部分，可以看出FAS从其他模式导波中完全分离了出来，同时对其它信号成分没有显著影响。
图5为本发明实施例的慢度-截距图中A0波截取过程示意图。
A0波的特征为沿骨板仿体传播的A0波的速度大约为1.0-1.8km/s，因此，选择截取的A0波区域为慢度范围为0.5-1.0s/km、截距范围为38.0-45.0μs的区域(如图5中a部分圆圈部分所示)。
图5a为慢度-截距图，横坐标为慢度(s/km)，纵坐标为时间(μs)，从图5中a部分圆圈部分(即A0波区域)可以看出，A0波位于整个慢度-截距图中各个信号区域的最右端，其余部分是超声信号传播过程中的其他模式导波。
为了展示A0波区域截取前后的具体变化，本实施例中也选取65、70、75、80以及85号接收阵元各自对应的降噪处理后信号进行A0波区域截取前后的具体分析。
图5b为降噪处理后信号(即降噪后信号)距离-时间图，横坐标为距离(mm)，纵坐标为时间(μs)。
如图5b部分所示，5根实线分别为65、70、75、80以及85号接收阵元对应的降噪处理后信号，在该降噪处理后信号中，A0波与其他模式导波混杂在一起。
图5中c部分与d部分分别为其他模式导波(即A0波以外信号)对应的距离-时间域图以及A0波区域(即A0波)对应的距离-时间域图，对比图5的b、c以及d部分，可以看出A0波从其他模式导波中完全分离了出来，同时对其它模式导波没有显著影响。
速度计算模块15基于FAS波区域计算得到第一到达的传播速度作为FAS波速度，基于A0波区域计算得到基阶反对称模态Lamb波的传播速度作为A0波速度，并基于FAS波区域以及A0波区域计算得到实际频散曲线。
图6为本发明实施例的包络信号示意图。
其中，FAS波速度通过如下步骤计算得到：
步骤S1-1，对FAS波区域进行Radon正变换得到第一到达波在距离-时间域上的投影作为FAS波投影。
步骤S1-2，对FAS波投影进行Hilbert变换得到第一到达波的包络信号作为第一包络信号。
为了展示包络信号，以图6中a部分作为展示，横坐标为时间，纵坐标为幅度，其中虚线为单个接收阵元对应的包络信号，实线为各种导波信号在时间-幅度域上的投影。
图6b部分为65、70、75、80以及85号接收阵元对应的第一包络信号(即FAS包络)，横坐标为距离(mm)(即发射阵元和接收阵元之间的距离)，纵坐标为时间(μs)。
步骤S1-3，根据接收阵元与发射阵元之间的距离以及第一包络信号到达接收阵元的时间计算得到的速度作为FAS波速度。
第一包络信号到达接收阵元的时间可以是第一包络信号达到最大值的时间(如图6a中虚线上圆形点对应的时刻)，也可以是第一包络信号达到最大值后幅度降为1/2的时间(如图6a中虚线上三角形点对应的时刻)，还可以是第一包络信号达到最大值后幅度降为0的时间(如图6a中虚线上正方形点对应的时刻)。
本实施例中，选择第一包络信号达到最大值的时间。
FAS波速度除了通过上述方式计算得到，还可以通过选取FAS波区域中幅度最大的点作为第一幅度最大点，将该第一幅度最大点对应的慢度值的倒数作为FAS波速度得到。还可以通过选取FAS波区域中截距值最小的点作为第一截距最小点，将该第一截距最小点对应的慢度值的倒数作为FAS波速度得到。
将图6b中5个第一包络信号中凸起点连成一条直线，该直线即为最大值连线，即FAS到达时刻的连线，计算得到该虚线的斜率为0.26s/km，将该斜率的倒数3.85km/s作为FAS波速度，在骨板仿体FAS的速度范围1.0-1.8km/s内。
A0波速度通过如下步骤计算得到：
步骤S2-1，对A0波区域进行Radon正变换得到基阶反对称模态Lamb波在距离-时间域上的投影作为A0波投影。
步骤S2-2，对A0波投影进行Hilbert变换得到基阶反对称模态Lamb波的包络信号作为第二包络信号。
图6c部分为65、70、75、80以及85号接收阵元对应的第二包络信号(即A0包括)，横坐标为距离(mm)，纵坐标为时间(μs)。
步骤S2-3，根据接收阵元与发射阵元之间的距离以及第二包络信号到达接收阵元的时间计算得到的速度作为A0波速度。
第二包络信号到达接收阵元的时间可以是第二包络信号达到最大值的时间，也可以是第二包络信号达到最大值后幅度降为1/2的时间，还可以是第二包络信号达到最大值后幅度降为0的时间。
本实施例中，选择第二包络信号达到最大值的时间。
A0波速度处理通过上述方法计算得到外，还可以通过选取A0波区域中幅度最大的点作为第二幅度最大点，将该第二幅度最大点对应的慢度值的倒数作为A0波速度得到。还可以通过选取A0波区域中截距值最小的点作为第二截距最小点，将该第二截距最小点对应的慢度值的倒数作为A0波速度得到。
将图6c中5个第二包络信号中凸起点连成一条直线，该直线即为最大值连线，即A0波到达时刻的连线，计算得到该虚线的斜率为0.63s/km，将该斜率的倒数1.50km/s作为A0波速度，在骨板仿体A0波的速度范围3.0-4.5km/s内。
长骨反演模块16基于FAS波速度以及A0波速度对待测长骨进行反演，从而得到该待测长骨对应的长骨几何参数以及弹性参数。
其中，长骨反演模块16包括代价函数设计单元、反演模型构建单元、纵波速度计算单元以及参数选定单元。
代价函数设计单元根据实际频散曲线以及待测长骨的理论频散曲线涉及代价函数C：
Figure BDA0003007491300000151
式中，(i,j)为实际频散曲线上一点的坐标，(p,q)为理论频散曲线上一点的坐标，A(i,j)为实际频散曲线上点(i,j)的幅度。
反演模型构建单元选取参数范围作为预估范围，根据该预估范围构建反演模型。
纵波速度计算单元利用反演模型计算FAS波速度在不同厚度下的纵波速度。
本实施例中，待测长骨的密度为1.85g/cm3，选取的反演参数为厚度以及横波速度，厚度区间为1.780-1.820mm，步长为0.001mm。横波速度区间为1.71-1.90km/s，步长为0.01km/s。
此时，纵波速度＝FAS速度/[0.25×(厚度×激励频率/4000)0.2+0.75]，可以应用在骨板厚度小于超声波长情境下。
图7为本发明实施例的反演结构示意图。
参数选定单元根据代价函数C以及纵波速度计算得到不同参数下实际频散曲线与理论频散曲线之间的差值，将差值最小对应的理论频散曲线的参数作为目标参数，从目标参数中获取长骨几何参数以及弹性参数。
图7a为不同参数下测量值与理论值的误差图，横轴为横波速度，纵轴为厚度，竖轴为代价函数值，从图7a可以看出，代价函数存在最小值(即图7a中三角形点)，该最小值对应的点的坐标即为最接近原始信号的反演结果。
图7b为最优厚度下不同横波速度大小的误差图，横坐标为横波速度，纵坐标为代价函数值，本实施例中，最优厚度为1.813mm，代价函数存在最小值(如图7b中三角形点)。
图7c为最优横波速度下不同厚度大小的误差图，横坐标为厚度，纵坐标为代价函数值，本实施例中，最优横波速度为1.79km/s，代价函数存在最小值(如图7c中三角形点)。
对比图7b与图7c，厚度为1.813mm且横波速度为1.79km/s时代价函数值最小，即反演结果最接近原始信号。
其中，几何参数包括但不限于厚度，弹性参数包括但不限于横波速度、纵波速度、弹性模量、泊松比、拉梅常数。
本实施例中，目标参数对应的骨板仿体的厚度为1.813mm，横波速度为1.79km/s，纵波速度为3.99km/s。而骨板仿体的实际值厚度为1.8mm，实际横波速度为1.8km/s，实际纵波速度为4.0km/s。对比上述目标参数与实际参数可以看出，目标参数十分接近实际参数，从而说明本发明的基于Radon变换的长骨相控超声信号定征与骨质评价系统1有较高的准确性。
另外，弹性模量根据骨板仿体的密度计算得到：弹性模量＝密度×速度2。
其中，骨板仿体的密度为1.85g/cm3，横波决定的弹性模量参数C55为5.93GPa，纵波决定的弹性模量参数C33为29.49GPa。
骨质评价模块17基于长骨几何参数以及弹性参数，对待测长骨的骨质状况进行评价从而得到评价结果。
图8为本发明实施例的基于Radon变换的长骨相控超声信号定征与骨质评价系统工作过程的流程图。
如图8所示，基于Radon变换的长骨相控超声信号定征与骨质评价系统1工作过程包括如下步骤：
步骤S3-1，信号采集模块11利用超声探头采集到待测长骨对应的原始信号，然后进入步骤S3-2；
步骤S3-2，降噪模块12利用Radon正反变换对对原始信号进行处理得到降噪处理后信号，然后进入步骤S3-3；
步骤S3-3，慢度截距图生成模块13对降噪处理后信号进行Radon反变换得到Radon域上的慢度-截距图，然后进入步骤S3-4；
步骤S3-4，区域截取模块14分别基于第一到达波与基阶反对称模态Lamb波的特征从慢度-截距图中截取出FAS波区域与A0波区域，然后进入步骤S3-5；
步骤S3-5，速度计算模块15分别基于FAS波区域与A0波区域计算得到FAS波速度与A0波速度以及实际频散曲线，然后进入步骤S3-6；
步骤S3-6，长骨反演模块16基于FAS波速度以及A0波速度对待测长骨进行反演得到长骨几何参数以及弹性参数，然后进入步骤S3-7；
步骤S3-7，骨质评价模块17基于长骨几何参数以及弹性参数，对待测长骨的骨质状况进行评价从而得到评价结果，然后进入结束状态。
实施例作用与效果
根据本实施例提供的基于Radon变换的长骨相控超声信号定征与骨质评价系统1，由于降噪模块利用Radon正反变换对对原始信号进行处理得到降噪处理后信号，因此降噪处理后信号具有较高的信噪比。另外，由于区域截取模块基于第一到达波与基阶反对称模态Lamb波的特征从慢度-截距图中截取出FAS波区域与A0波区域，因此解决了傅里叶变换不能很好地提取出FAS和A0波信号的问题。除此之外，还由于速度计算模块基于FAS波区域与A0波区域计算得到FAS波速度与A0波速度，因此，便于后续根据FAS波速度与A0波速度精确反演得到待测长骨的几何参数以及弹性参数，更进一步地，提高了基于几何参数以及弹性参数的骨质评估结果的准确性。
实施例中，由于发射阵元与接收阵元分别设置在探头壳体的两端，因此可以有效避免激励沿探头传播而不是沿骨板传播对采集信号产生的影响。
上述实施例仅用于举例说明本发明的具体实施方式，而本发明不限于上述实施例的描述范围。

Claims (10)

1.一种基于Radon变换的长骨相控超声信号定征与骨质评价系统，用于对待测长骨进行定征与评价，其特征在于，包括：
信号采集模块，利用超声探头从所述待测长骨中采集到原始信号；
降噪模块，利用Radon正反变换对对所述原始信号进行处理得到降噪处理后信号；
慢度截距图生成模块，对所述降噪处理后信号进行Radon反变换，得到Radon域上的慢度-截距图；
区域截取模块，利用第一到达波的特征从所述慢度-截距图中截取出FAS波区域，利用基阶反对称模态Lamb波的特征从所述慢度-截距图中截取出A0波区域；
速度计算模块，基于所述FAS波区域计算得到所述第一到达的传播速度作为FAS波速度，基于所述A0波区域计算得到所述基阶反对称模态Lamb波的传播速度作为A0波速度，并基于所述FAS波区域以及所述A0波区域计算得到实际频散曲线；
长骨反演模块，基于所述FAS波速度以及所述A0波速度对所述待测长骨进行反演，从而得到该待测长骨对应的长骨几何参数以及弹性参数；以及
骨质评价模块，基于所述长骨几何参数以及所述弹性参数，对所述待测长骨的骨质状况进行评价从而得到评价结果。
2.根据权利要求1所述的基于Radon变换的长骨相控超声信号定征与骨质评价系统，其特征在于：
其中，所述超声探头包括探头壳体、至少一个发射阵元以及多个接收阵元，
所述发射阵元与所述接收阵元分别设置在所述探头壳体的两端，
所述发射阵元用于对所述待测长骨发射超声信号，
所述接收阵元用于接收所述待测长骨因所述超声信号产生的信号并作为所述原始信号。
3.根据权利要求2所述的基于Radon变换的长骨相控超声信号定征与骨质评价系统，其特征在于：
其中，所述FAS波速度通过如下步骤计算得到：
步骤S1-1，对所述FAS波区域进行Radon正变换得到所述第一到达波在距离-时间域上的投影作为FAS波投影；
步骤S1-2，对所述FAS波投影进行Hilbert变换得到所述第一到达波的包络信号作为第一包络信号；
步骤S1-3，根据所述接收阵元与所述发射阵元之间的距离以及所述第一包络信号到达所述接收阵元的时间计算得到的速度作为所述FAS波速度，
所述A0波速度通过如下步骤计算得到：
步骤S2-1，对所述A0波区域进行Radon正变换得到所述基阶反对称模态Lamb波在距离-时间域上的投影作为A0波投影；
步骤S2-2，对所述A0波投影进行Hilbert变换得到所述基阶反对称模态Lamb波的包络信号作为第二包络信号；
步骤S2-3，根据所述接收阵元与所述发射阵元之间的距离以及所述第二包络信号到达所述接收阵元的时间计算得到的速度作为所述A0波速度。
4.根据权利要求3所述的基于Radon变换的长骨相控超声信号定征与骨质评价系统，其特征在于：
其中，所述第一包络信号到达所述接收阵元的时间为第一包络信号达到最大值的时间，
所述第二包络信号到达所述接收阵元的时间为第二包络信号达到最大值的时间。
5.根据权利要求3所述的基于Radon变换的长骨相控超声信号定征与骨质评价系统，其特征在于：
其中，所述第一包络信号到达所述接收阵元的时间为第一包络信号达到最大值后幅度降为1/2的时间，
所述第二包络信号到达所述接收阵元的时间为第二包络信号达到最大值后幅度降为1/2的时间。
6.根据权利要求3所述的基于Radon变换的长骨相控超声信号定征与骨质评价系统，其特征在于：
其中，所述第一包络信号到达所述接收阵元的时间为第一包络信号达到最大值后幅度降为0的时间，
所述第二包络信号到达所述接收阵元的时间为第二包络信号达到最大值后幅度降为0的时间。
7.根据权利要求1所述的基于Radon变换的长骨相控超声信号定征与骨质评价系统，其特征在于：
其中，所述FAS波速度的计算方法为选取所述FAS波区域中幅度最大的点作为第一幅度最大点，将该第一幅度最大点对应的慢度值的倒数作为所述FAS波速度，
所述A0波速度的计算方法为选取所述A0波区域中幅度最大的点作为第二幅度最大点，将该第二幅度最大点对应的慢度值的倒数作为所述A0波速度。
8.根据权利要求1所述的基于Radon变换的长骨相控超声信号定征与骨质评价系统，其特征在于：
其中，所述FAS波速度的计算方法为选取所述FAS波区域中截距值最小的点作为第一截距最小点，将该第一截距最小点对应的慢度值的倒数作为所述FAS波速度，
所述A0波速度的计算方法为选取所述A0波区域中截距值最小的点作为第二截距最小点，将该第二截距最小点对应的慢度值的倒数作为所述A0波速度。
9.根据权利要求1所述的基于Radon变换的长骨相控超声信号定征与骨质评价系统，其特征在于：
其中，所述长骨反演模块包括：
代价函数设计单元，根据所述实际频散曲线以及所述待测长骨的理论频散曲线涉及代价函数C：
Figure FDA0003007491290000051
式中，(i,j)为所述实际频散曲线上一点的坐标，(p,q)为所述理论频散曲线上一点的坐标，A(i,j)为所述实际频散曲线上点(i,j)的幅度；
反演模型构建单元，选取参数范围作为预估范围，根据该预估范围构建反演模型；
纵波速度计算单元，利用所述反演模型计算所述FAS波速度在不同厚度下的纵波速度；以及
参数选定单元，根据所述代价函数C以及所述纵波速度计算得到不同参数下所述实际频散曲线与所述理论频散曲线之间的差值，将差值最小对应的理论频散曲线的参数作为目标参数，从所述目标参数中获取所述长骨几何参数以及所述弹性参数。
10.根据权利要求1所述的基于Radon变换的长骨相控超声信号定征与骨质评价系统，其特征在于：
其中，所述降噪模块包括Radon反变换单元以及Radon正变换单元，
所述Radon反变换单元对所述原始信号进行Radon反变换得到反变换后信号，
所述Radon正变换单元对所述反变换后信号进行Radon正变换得到正变换后信号，并作为所述降噪处理后信号。
CN202110366930.9A 2021-04-06 2021-04-06 基于Radon变换的长骨相控超声信号定征与骨质评价系统 Active CN113109825B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110366930.9A CN113109825B (zh) 2021-04-06 2021-04-06 基于Radon变换的长骨相控超声信号定征与骨质评价系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110366930.9A CN113109825B (zh) 2021-04-06 2021-04-06 基于Radon变换的长骨相控超声信号定征与骨质评价系统

Publications (2)

Publication Number Publication Date
CN113109825A true CN113109825A (zh) 2021-07-13
CN113109825B CN113109825B (zh) 2022-06-14

Family

ID=76713959

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110366930.9A Active CN113109825B (zh) 2021-04-06 2021-04-06 基于Radon变换的长骨相控超声信号定征与骨质评价系统

Country Status (1)

Country Link
CN (1) CN113109825B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116098655A (zh) * 2023-04-11 2023-05-12 湖南工商大学 基于超声导波多重信号分类的骨骼参数检测装置和方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050004457A1 (en) * 2001-11-30 2005-01-06 Petro Moilanen Method and device for the non-invasive assessement of bones
US20090236145A1 (en) * 2008-03-20 2009-09-24 Schlumberger Technology Corporation Analysis refracted acoustic waves measured in a borehole
CN103906474A (zh) * 2011-11-01 2014-07-02 骨治医疗公司 利用电磁波的骨骼方法和安排
CN106597539A (zh) * 2016-12-28 2017-04-26 中国石油化工股份有限公司 针对黄土塬地区的曲波域Radon变换噪声压制方法
CN110068613A (zh) * 2019-04-09 2019-07-30 南京信息职业技术学院 一种结构导波响应群速度频散测试方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050004457A1 (en) * 2001-11-30 2005-01-06 Petro Moilanen Method and device for the non-invasive assessement of bones
US20090236145A1 (en) * 2008-03-20 2009-09-24 Schlumberger Technology Corporation Analysis refracted acoustic waves measured in a borehole
CN103906474A (zh) * 2011-11-01 2014-07-02 骨治医疗公司 利用电磁波的骨骼方法和安排
CN106597539A (zh) * 2016-12-28 2017-04-26 中国石油化工股份有限公司 针对黄土塬地区的曲波域Radon变换噪声压制方法
CN110068613A (zh) * 2019-04-09 2019-07-30 南京信息职业技术学院 一种结构导波响应群速度频散测试方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
KAILIANG XU ET AL.: "Multichannel processing for dispersion curves extraction of ultrasonic axial-transmission signals: Comparisons and case studies", 《ACOUSTICAL SOCIETY OF AMERICA》 *
刘丹 等: "超声导波定量评价长骨骨折状况", 《中国科学》 *
楼思佳 等: "时频域独立元分析方法实现超声兰姆波的多模式分离", 《声学学报》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116098655A (zh) * 2023-04-11 2023-05-12 湖南工商大学 基于超声导波多重信号分类的骨骼参数检测装置和方法

Also Published As

Publication number Publication date
CN113109825B (zh) 2022-06-14

Similar Documents

Publication Publication Date Title
Droin et al. Velocity dispersion of acoustic waves in cancellous bone
EP0123427B1 (en) Ultrasonic medium characterization
Nicholson et al. A comparison of time-domain and frequency-domain approaches to ultrasonic velocity measurement in trabecular bone
US7601120B2 (en) Method and device for the non-invasive assessment of bones
JPH026022B2 (zh)
EP2939598A1 (en) Ultrasonic diagnostic device and elasticity evaluation method
CN110261485B (zh) 一种超声波测量材料内部各处弹性模量及泊松比的方法
Moreau et al. Measuring the wavenumber of guided modes in waveguides with linearly varying thickness
TW200817659A (en) Ultrasonic stress measuring apparatus
Moilanen et al. Ultrasonically determined thickness of long cortical bones: two-dimensional simulations of in vitro experiments
CN105300856A (zh) 基于超声阻抗谱对颗粒浓度和尺寸的测量方法
Wear Cancellous bone analysis with modified least squares Prony’s method and chirp filter: Phantom experiments and simulation
CN113109825B (zh) 基于Radon变换的长骨相控超声信号定征与骨质评价系统
Kehlenbach et al. Identifying damage in plates by analyzing Lamb wave propagation characteristics
CN108852416B (zh) 一种剪切波传播速度的确定方法及装置
Shi et al. Resolution enhancement of ultrasonic imaging at oblique incidence by using WTFM based on FMC-AR
WO2019025510A1 (en) METHOD AND DEVICE FOR CHARACTERIZING A WAVEGUIDE
CN110045014B (zh) 基于贝叶斯学习的Lamb波频散消除方法及其系统
Laaboubi et al. Application of the reassignment time–frequency method on an acoustic signals backscattered by an air-filled circular cylindrical shell immersed in water
Wear Estimation of fast and slow wave properties in cancellous bone using Prony's method and curve fitting
EP2366997A1 (en) Method and device for determining the structural organization of an object with ultrasounds
JP2001337077A (ja) コンクリート構造物の剥離の非破壊検査方法
JP2002139478A (ja) 構造材料のクリープ損傷検出方法及び装置
Bauer et al. Bone sonometry: Reducing phase aberration to improve estimates of broadband ultrasonic attenuation
JP7302651B2 (ja) 超音波信号処理装置、超音波診断装置、超音波信号処理方法、およびプログラム

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
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20240321

Address after: Room 602, No. 10, Lane 999, Dangui Road, China (Shanghai) Pilot Free Trade Zone, Pudong New Area, Shanghai, March 2012 (property registration floor is 5th floor)

Patentee after: Shanghai Shengai Medical Technology Co.,Ltd.

Country or region after: Zhong Guo

Address before: 200433 No. 220, Handan Road, Shanghai, Yangpu District

Patentee before: FUDAN University

Country or region before: Zhong Guo

Patentee before: SHANGHAI A & S SCIENCE TECHNOLOGY DEVELOPMENT Co.,Ltd.