技术领域
本发明涉及医学超声成像技术领域，特别涉及一种定量剪切波弹性成像方法。
背景技术
近年来，肿瘤疾病已经成为危害我国人民身体健康的主要杀手。肿瘤疾病的防治关键在于早期发现。临床上，人体组织的力学特征变化往往是肿瘤疾病最为重要的早期预警信号，如乳腺、甲状腺等肿瘤病变。随着肿瘤的生长，病变组织与正常组织相比，黏弹性发生较大变化，如乳腺癌与周围正常组织大约相差90倍，纤维化组织、非浸润性癌变和浸润性癌变组织的剪切弹性系数亦存在较大差异，因此，获得人体组织与粘弹性相关的参数信息(例如形变位移，剪切波速等)在医学诊断领域尤其对乳腺癌、甲状腺癌和前列腺癌的早期检测具有重要意义。
医学超声弹性成像是一种通过组织运动显示组织弹性模量或软硬程度的超声成像方式，传统的挤压式弹性成像方式需要医生利用探头挤压检测部位使其受力运动，从而得到组织弹性信息，这种方式的不足在于：医生需要以适当的操作挤压检测部位，主观性强；由于需要医生挤压，不同医生操作或同一医生不同时期的操作的结果很难直接进行对比，检测结果难以复现，往往只能作为定性的检测结果，无法获得定量的检测信息，用以疾病的跟踪和术后观察。
目前能进行定量分析的超声弹性成像技术，主要是基于声辐射力弹性成像(Acoustic Radiation Force Imaging，ARFI)，ARFI利用医学超声功率范围内的聚焦超声激励脉冲在生物黏弹性组织聚焦区域产生声辐射力，使组织发生形变，然后利用检测脉冲的回波信号通过相关性的时延计算方法在不同时间点检测组织的形变位移情况，将其成像后定性地反映组织的黏弹性特性。ARFI弹性成像克服了传统徒手弹性成像无法从体外对深部组织有效施压以及因操作者的使用习惯不同造成可重复性较差的缺点。但不足之处在于图像显示的是组织相对的形变位移差异，而无法完全定量地估算组织的弹性模量。基于ARFI技术衍生的超声剪切波成像(SWI，shear wave imaging)，通过调整脉冲激励方式，根据水平方向多个点的位移-时间曲线估算横向传播的剪切波波速，然后利用剪切波传播特性与生物组织弹性特征之间的内在联系，最终定量重构组织的弹性模量以形成二维的图像。
但是，现有的定量弹性成像SWI需要使用超声前端波束合成前的原始射频信号数据，需要基于平面波发射和接收，平面波技术每秒需要N*c/2z次计算，其中N为并行接收声束数目，c为声速，z为扫描深度。而传统波束合成方法的计算量则要小得多，仅为M/P分之一，其中M为发射次数，P为并行数，例如M等于100，P等于2时，平面波方法的计算量是传统ARFI方法的至少五十倍(有时甚至高达数百倍)；除了计算量的巨大差异外，基于平面波技术的接收端软件波束合成还需要存储大量原始回波信号，大幅度提高了数据存储和传输硬件成本，不便于与医学超声系统集成和推广应用。
发明内容
本发明的目的在于克服现有定量弹性成像技术需要采用原始射频信号(不经过波束合成)带来的计算量大、所需硬件设计困难和成本高昂的问题，提供一种对设备要求低、计算量小的超声定量弹性成像方法及系统。
为了实现上述发明目的，本发明提供了以下技术方案：
一种定量剪切波弹性成像方法，包含如下步骤：
进行ARFI检测，获取目标区域中各个位置应变ε，其中包含指定参考位置应变εref；
进行剪切波检测，获取目标区域中指定参考位置单点剪切波波速cref；
根据所述单点剪切波波速计算各个位置剪切波波速c；
根据各个位置应变计剪切波波速获取各个位置定量剪切波弹性图像E。
进一步的，所述目标区域中各位置应变其中，N为梯度滑动计算窗口大小，zi表示各位置的具体深度坐标，表示梯度滑动计算窗口内各个位置坐标的均值，yi表示梯度滑动计算窗口内当前点位移，表示梯度滑动计算窗口内各个位置位移的均值。
进一步的，所述单点剪切波波速根据指定测量点的横向距离、轴向距离，及各个指定测量点达到最大横向位移所需要的时间采用二维线性拟合法计算得出。
进一步的，所述单点剪切波波速cref通过公式X＝Aβ得出，其中，式中，xi表示指定测量点的横向距离，ti表示指定测量点达到最大位移所需时间，zi表示指定测量点的轴向距离，β1即为所求的指定测量点单点剪切波速cref，β0和β2为利用上式计算单点剪切波速cref过程中同时计算得出的回归参数，本方法中不使用这两个参数。
进一步的，所述各个位置剪切波波速其中，在聚焦区域，γ取值1；在聚焦区域之外的轴向区域，其中，z为轴向距离，n和σ为表示激励声场的系统常数，取值范围均为0～10，n和σ两个参数的取值可以相同或不同。
进一步的，所述聚焦区域DOF＝8(f#)2λ，其中，λ表示波长，f#为孔径控制参数，其为在0～5之间的实数，可见聚焦区域通常由设备参数决定，其大小范围只与设备参数相关，与目标区域大小无关，目标区域可能大于聚焦区域，也可能等于聚焦区域，或者小于聚焦区域。
进一步的，所述定量剪切波弹性图像E＝3ρc2，其中，ρ为被测区域介质密度。
以上方法可以运行在如下的超声定量弹性成像系统中，该系统包括，
剪切波脉冲发生器，用于发射剪切波激励脉冲及检测脉冲信号；
ARFI脉冲发生器，用于发射ARFI激励脉冲及检测脉冲信号；
探头，用于接收剪切波回波信号或ARFI回波信号；
控制装置，用于控制剪切波脉冲发生器、ARFI脉冲发生器的开关，接收剪切波回波信号和ARFI回波信号，并，根据ARFI回波信号计算目标区域各个位置的应变ε，其中包含指定参考位置应变εref；根据剪切波回波信号计算指定参考位置单点剪切波波速cref；根据单点剪切波波速计算各个位置剪切波波速c；根据各个位置应变计剪切波波速c及应变εref获取各个位置定量剪切波弹性图像E；
显示装置，用于将所述定量剪切波弹性图像E显示。
所述探头包括波束合成器，其用于将ARFI回波信号或剪切波回波信号经过波束合成，合成为回波射频信号。
所述控制装置包括剪切波信号处理模块、单点剪切波速计算模块；
所述剪切波信号处理模块自探头接收剪切波回波射频信号，并经过解调、位移估算得出被测组织目标区域各位置或者仅获取指定参考位置的剪切波位移-时间数据；
所述单点剪切波速计算模块根据剪切波位移-时间数据计算指定参考位置剪切波波速cref，其通过公式X＝Aβ得出，其中，式中，xi表示指定测量点(指定测量点包含指定参考位置)的横向距离，ti表示指定测量点达到最大位移所需时间，zi表示指定测量点的轴向距离，β1即为所求的指定测量点单点剪切波速cref；
所述控制装置还包括ARFI信号处理模块、应变计算模块；
所述ARFI信号处理模块用于自探头接收ARFI回波射频信号，并经过解调、位移估算得出被测组织目标区域各位置的ARFI位移-时间数据；
所述应变计算模块根据所述ARFI位移-时间数据计算目标区域各位置应变其中，N为梯度滑动计算窗口大小，zi表示各位置具体深度坐标，表示梯度滑动计算窗口内各个位置坐标的均值，yi表示梯度滑动计算窗口内当前点位移，表示梯度滑动计算窗口内各个位置位移的均值，ε中包含指定参考位置的应变εref。
所述控制装置还包括弹性模量计算模块；
所述控制装置还包括波速计算模块及弹性模量计算模块；
所述波速计算模块用于根据指定参考位置单点剪切波波速计算各个位置剪切波波速其中，在聚焦区域，γ取值1；在聚焦区域之外的轴向区域，其中，z为轴向距离，n和σ为表示激励声场的系统常数，取值范围均为0～10，n和σ两个参数的取值可以相同或不同；
所述弹性图像计算模块用于根据c及应变εref获取各个位置定量剪切波弹性图像E＝3ρc2，其中，ρ为被测区域介质密度。
所述控制装置控制所述剪切波脉冲发生器、ARFI脉冲发生器交替发射。
所述剪切波脉冲发生器及所述ARFI脉冲发生器为同一电路实现，其由所述控制装置控制发送剪切波脉冲或ARFI脉冲；，
所述剪切波脉冲发生器及所述ARFI脉冲发生器为各自独立电路实现，其由所述控制装置控制启动、关闭时间。
所述剪切波脉冲发生器和/或ARFI脉冲发生器共有24路以上。
所述剪切波脉冲发生器或ARFI脉冲发生器包括依次串接的驱动放大器、脉冲发生器、过流过压保护电路、Tx/Rx开关，以及自Tx/Rx开关反馈的限幅放大器、数模转换电路。
与现有技术相比，本发明的有益效果：
本发明提供的超声定量弹性成像方法基于滑动窗口线性拟合的应变和使用二维线性拟合的剪切波波速检测算法，其结果鲁棒性更好、抗噪声能力更强，结果更可靠。同时在不额外增加超声前端存储与传输模块负荷的情况下，实现了超声全域定量弹性成像，大幅地减低了超声定量弹性成像系统的设计难度和设备成本。
本发明方法运行的系统中设置有有效电路切换和保护模块，让系统满足ARFI与剪切波不同发射时序要求的同时，有效保护人体和系统电路，并最大限度的进行声激励以达到更好检测效果。
附图说明：
图1为本发明提供的超声定量弹性成像方法流程图。
图2是本发明提供的超声定量弹性成像系统框图。
图3是本发明中剪切波脉冲发生器和/或所述ARFI脉冲发生器构成示例的结构框图。
图4是本发明中本发明中声激励使用的发射脉冲序列示意图。
图5是剪切波波速计算使用的脉冲发射方式示意图。
具体实施方式
下面结合附图及具体实施例对本发明作进一步的详细描述。但不应将此理解为本发明上述主题的范围仅限于以下的实施例，凡基于本发明内容所实现的技术均属于本发明的范围。
实施例1：如图1所示，本发明提供一种定量剪切波弹性成像方法，包含如下步骤：
S100：进行ARFI检测，获取目标区域中各个位置应变ε，其中包含指定参考位置应变εref；
S200：进行剪切波检测，获取目标区域中指定参考位置单点剪切波波速cref；
步骤S100和S200的执行顺序并无特殊要求，如可以先执行S100，再执行S200，也可以先执行S200，再执行S100.
S300：根据所述单点剪切波波速计算各个位置剪切波波速c；
S400：根据各个位置应变计剪切波波速获取各个位置定量剪切波弹性图像E。
具体的，所述目标区域中各位置应变其中，N为梯度滑动计算窗口大小，zi表示各位置具体深度坐标，表示梯度滑动计算窗口内各个位置坐标的均值，yi表示梯度滑动计算窗口内当前点位移，表示梯度滑动计算窗口内各个位置位移的均值。
所述单点剪切波波速根据指定测量点的横向距离、轴向距离，及各个指定测量点达到最大横向位移所需要的时间采用二维线性拟合法计算得出。
所述单点剪切波波速cref通过公式X＝Aβ得出，其中，式中，通常我们在轴向深度上预设I个深度，每个深度上预设J个指定位置(点)，共，m个指定位置(点)，m＝I*J；我们将每个位置(点)在剪切波作用下达到最大位移的时间看做是剪切波传播到该位置(点)的时间，通过剪切波位移-时间数据可以得到剪切波传递到该位置(点)的时间，即ti；同时，zi表示指定测量位置(点)在轴向深度的轴向距离(表征深度位置)，xi表示指定测量位置(点)的横向距离(表征横向位置)，即，当各个指定位置(点)选定时，X中各个元素值以及zi值便已确定，我们只需从剪切波位移-时间数据中找出各个指定位置(点)达到最大位移所需时间ti即可完成计算；公式中，β1即为所求的指定测量点单点剪切波速cref；同时，应注意的是，β0和β2为利用上式计算单点剪切波速cref过程中同时计算得出的回归参数，本方法中并不使用这两个参数，其作用和含义对本发明没有实质意义。
所述各个位置剪切波波速其中，在聚焦区域γ取值1；在聚焦区域之外的轴向区域(当目标区域大于聚焦区域时，需用到此处)，其中，z为轴向距离，n和σ为表示激励声场的系统常数，取值范围均为0～10，n和σ两个参数的取值可以相同或不同。
所述定量剪切波弹性图像E＝3ρc2，其中，ρ为被测区域介质密度。
本实施例提供的方法可以运行在如图2、图3所示的超声定量弹性成像系统中，该系统包括，
剪切波脉冲发生器，用于发射剪切波激励脉冲及检测脉冲信号；
ARFI脉冲发生器，用于发射ARFI激励脉冲及检测脉冲信号；
探头，用于接收剪切波回波信号或ARFI回波信号；
控制装置1，用于控制剪切波脉冲发生器、ARFI脉冲发生器的开关，接收剪切波回波信号和ARFI回波信号，并，根据ARFI回波信号计算目标区域各个位置的应变ε，其中包含指定参考位置应变εref；根据剪切波回波信号计算指定参考位置单点剪切波波速cref；根据单点剪切波波速计算各个位置剪切波波速c；根据各个位置应变计剪切波波速c及应变εref获取各个位置定量剪切波弹性图像E。
显示装置，用于将所述定量剪切波弹性图像E显示。具体的，在显示之前，还包括对控制装置计算得出的弹性图像E数据进行扫描转换，平滑处理等一系列的处理过程，最后将该图像通过显示器展现给用户。
所述探头包括波束合成器2，其用于将ARFI回波信号或剪切波回波信号经过波束合成，合成为回波射频信号。
所述控制装置1包括剪切波信号处理模块13、单点剪切波速计算模块14；
所述剪切波信号处理模块13自探头接收剪切波回波射频信号，并经过解调成为IQ数据，进一步对该IQ数据进行位移估算得出被测组织目标区域各位置的剪切波位移-时间数据；
所述单点剪切波速计算模块14根据剪切波位移-时间数据计算指定参考位置剪切波波速cref，其通过公式X＝Aβ得出，其中，式中，xi表示指定测量点的横向距离，ti表示指定测量点达到最大位移所需时间，zi表示指定测量点的轴向距离，β1即为所求的指定测量点单点剪切波速cref；
所述控制装置1还包括ARFI信号处理模块11、应变计算模块12；
所述ARFI信号处理模块11用于自探头接收ARFI回波射频信号，并经过解调形成IQ数据，并进一步通过对IQ数据进行位移估算得出被测组织目标区域各位置的ARFI位移-时间数据；
所述应变计算模块12根据所述ARFI位移-时间数据计算目标区域各位置应变其中，N为梯度滑动计算窗口大小，zi表示各位置(此处的各位置通常也是从目标区域中选取的)具体深度坐标，表示梯度滑动计算窗口内各个位置坐标的均值，yi表示梯度滑动计算窗口内当前点位移，表示梯度滑动计算窗口内各个位置位移的均值。
所述控制装置1还包括波速计算模块15及弹性模量计算模块16；
所述波速计算模块15用于根据指定参考位置单点剪切波波速计算各个位置剪切波波速其中，在聚焦区域，γ取值1；在聚焦区域之外的轴向区域，其中，z为轴向距离，n和σ为表示激励声场的系统常数，取值范围均为0～10，n和σ两个参数的取值可以相同或不同；
所述弹性图像计算模块16用于根据c及应变εref获取各个位置定量剪切波弹性图像E＝3ρc2，其中，ρ为被测区域介质密度。
所述控制装置1控制所述剪切波脉冲发生器、ARFI脉冲发生器交替发射。
具体的，如图3所述，一些实施例中，控制装置1硬件部分可由PC、控制电路板和FPGA芯片组组成，另外一些实施例中，控制装置也可以是完全集成在一起的集成电路构成，上述的剪切波信号处理模块、单点剪切波速计算模块、ARFI信号处理模块、应变计算模块、弹性模量计算模块按功能分布在上述硬件装置中。
所述剪切波脉冲发生器及所述ARFI脉冲发生器为同一电路实现，其由所述控制装置控制发送剪切波脉冲或ARFI脉冲；或，
所述剪切波脉冲发生器及所述ARFI脉冲发生器为各自独立电路实现，其由所述控制装置控制启动、关闭时间。
具体的，在所述剪切波脉冲发生器及所述ARFI脉冲发生器为同一电路实现的实施例中，如图3所示，发生器包括依次串接的驱动放大器100、脉冲发生器101、过流过压保护电路102、Tx/Rx开关103，以及自Tx/Rx开关反馈的限幅放大器104、数模转换电路105。过流过压保护电路102可以让系统满足ARFI与剪切波不同发射时序要求的同时，有效保护人体和系统电路，并最大限度的进行声激励以达到更好检测效果。
所述剪切波脉冲发生器和/或ARFI脉冲发生器共有24路以上，如48路、64路、128路、256路等。
使用时，控制装置根据用户设定的电压、脉冲长度和相位等参数信息控制脉冲发生器(剪切波脉冲发生器和/或所述ARFI脉冲发生器)产生激励和检测脉冲信号，在T/R开关的控制下发射电路产生的超声波进入生物组织，接收电路接收回波信号。系统在发射或接收信号时，通过设置孔径控制参数f#及聚焦深度z，并利用公式来控制活动阵元(脉冲发生器)的个数，公式中，D是探头孔径大小，即用户通过设置孔径控制参数f#及聚焦深度z来控制活动探头孔径的大小来实现对活动阵元个数的控制，用户可根据情况在最大活动阵元数(如64路)内选择实际需要的活动阵元数量，例如，系统如果包括64路脉冲发生器，则系统可以实现的最大探头孔径大小即为64个阵元所分布的面积，孔径控制参数f#通常为大于0小于或等于5的实数。因为激励组织产生形变的力是由探头根据用户设定的选项发射的短时长脉冲波在聚焦区域形成的，与操作者的操作手法无关，因此可以认为在有效成像区域内力的大小是均匀的，一致的。聚焦区域的深度(Depth of Focus)可以描述为：DOF＝8(f#)2λ，λ表示波长；例如，对于5MHz的激励脉冲，如果孔径控制参数f#＝2，生物组织中的声速近似等于1540m/s，那么DOF的值接近1cm。ARFI弹性成像采用的是逐线方式扫描，每一个水平位置都会根据相同的系统参数发射激励脉冲使得组织产生微小形变，因此可以认为在焦点附近1cm范围内的声辐射力可以看作是均匀的，一致的，这个区域我们称为聚焦区域，如上所述，在聚焦区域中，剪切波波速计算公式中的γ取值为1。
发射脉冲序列示意如图4所示，根据ROI窗口深度位置，信号中心频率大小、PRF(pulse repeat frequency，脉冲重复频率)、激励电压以及脉冲长度等参数控制发射脉冲序列。序列的头部是1个或多个检测脉冲(Detect beam,较高电压(例如80V，一般与B模式相同)，短脉冲(例如2个周期))作为组织位移信息的参考信号；接下来是多组激励脉冲/检测脉冲对，激励脉冲(Push beam,较低电压(例如20～40V)，长脉冲(100～250个周期))用于在组织内焦点附近的局部区域激励产生微小位移。
检测脉冲(高电压的短脉冲)用于跟踪声辐射力载荷在加载的过程中组织的形变情况；之后是一系列的检测脉冲跟踪声辐射力载荷消失后组织的形变情况。激励脉冲的作用是使焦点位置处的局部组织在安全条件下能够振动起来，长回波信号由于其空间分辨率差不能用于位移估算，所以探头阵元在发射激励脉冲时关闭信号变迹，所有活动阵元的发射信号幅值相同，增强信号携带的能量以便在安全条件下产生更大的声辐射力。用于激励的长脉冲波频率和用于检测的短脉冲波中心频率可选择相同或不同，优选方案是中心频率不同，这样可以为了避免激励脉冲对检测脉冲回波信号的干扰，同时在信号处理时易于区分激励和检测的回波信号。
另外发射脉冲序列中的长脉冲波数目是可变的，整个序列的波数目也是可变的，这里我们将一次激励整个序列的波数目定义为取样容积数目(ensemble size)，例如16，24，32。长脉冲波数目(激励波)+短脉冲波数目(参考和检测波)＝取样容积数目，梯度滑动计算窗口N的大小为0～取样容积数之间的任意值，具体其为采用自适应算法自动选择，此处不再赘述。
一组剪切波发射脉冲包含若干个不同空间位置如图4所示的发射脉冲序列，其中激励脉冲空间位置固定，而不同发射序列的检测脉冲空间位置不同；如图5所示，D为检测波，P为激励波，假如有一组剪切波发射脉冲采用5个间隔为1mm的发射脉冲序列，这5个序列的激励波P1、P2、P3、P4、P5均在同一空间位置，而5个发射脉冲序列中的检测波D1、D2、D3、D4、D5在相同或不同空间位置，同时，至少一个检测脉冲的空间位置和激励脉冲位置相同，如图中的D1。
ARFI脉冲发生器使用的发射脉冲序列与剪切波脉冲发生器相同，但一帧ARFI成像扫描过程中激励波保持与检测波为同一空间位置。