CN103549978B - 一种瞬时弹性成像快速时移估计方法 - Google Patents

Abstract

本发明公开了一种瞬时弹性成像快速时移估计方法，该方法包括以下步骤：接收不同时刻的N个超声回波信号s₁，…，s_N；根据预先设定的窗长和重叠率，将所述N个超声回波信号s₁，…，s_N分为若干个信号片段，所述N个超声回波信号s₁，…，s_N上相同深度的信号片段分别表示为s_1k，…，s_Nk，计算所述信号片段s_2k相对所述超声回波信号s₁的时移d_2k，并以所述时移d_2k作为先验值计算所述信号片段s_3k相对所述超声回波信号s₁的时移d_3k，以此类推，直到以d_(n-1)k作为先验值计算所述信号片段s_nk和所述超声回波信号s₁的时移d_nk。本发明可以有效地缩小时移估计的范围，提高计算速度。

Description

一种瞬时弹性成像快速时移估计方法

技术领域

本发明涉及快速时移估计方法，尤其涉及一种瞬时弹性成像的快速时移估计方法。

背景技术

瞬时弹性成像是一种实用的间接测量组织弹性的方法，该方法的原理是根据剪切波在组织里面的传播速度与组织剪切模量的理论关系，利用测得的剪切波的传播速度来计算组织的弹性系数。其通常的做法是采用外振动器在组织表面产生50-200Hz的低频脉冲振动，低频振动在组织中会产生剪切波，在软组织里面剪切波的传播速度只有几到几十米每秒，而在组织成像里面通常所用的超声信号是一种压缩波，其波速约为1540米每秒。因此可以用超声波跟踪低频振动产生的剪切波在组织中的传播过程。通过对超声M模式信号的分析，来确定低频振动产生的剪切波在组织中的传播速度，从而计算组织的弹性系数。

在上述瞬时弹性成像方法中，需要从剪切波引起的组织应变图上得到剪切波的传播速度，因此分析不同时刻的超声回波信号，来确定剪切波传播造成的组织偏移是瞬时弹性成像信号处理过程中的关键步骤之一。对于每次低频脉冲振动，发射接收N次超声波，根据这N次超声回波信号进行位移/应变估计的步骤为：1)确定一定的窗长度和重叠率；2)按照确定的窗长和重叠率，确定每一对对应射频信号片段的时移，即相对位移；3)利用得到的位移计算组织的应变分布。

在第二步中，确定对应射频信号片段的时移时，通常是以第一根射频信号为基准，即确定射频信号片段与第一根射频信号上的对应信号片段间的时移。有两种实现方式，一种是直接确定射频信号片段与第一根射频信号上的对应信号片段间的时移；另一种是先确定相邻两根射频信号对应信号片段的时移，再通过累加得到相对于第一根射频信号对应信号片段的时移。由于组织的位移是连续的，所以第二种方式搜索时移的范围比较小，计算量小，但是由于受信号采样频率的限制，会造成误差累积，影响精度。虽然可以通过插值来解决这一问题，但是通过插值来提高采样频率会增加计算量。而第一种方式没有误差累积问题，但是时移的搜索范围较第二种方式大，计算量也随之增大。

发明内容

本发明的目的是，针对上述直接计算射频信号片段相对于第一根射频信号的对应信号片段的时移计算量大的问题，提供一种用于瞬时弹性成像快速时移估计方法，利用组织运动连续性，缩小搜索范围，降低计算量，避免误差累积问题或插值带来的计算量增大问题。

为了实现上述目的，本发明提供了一种瞬时弹性成像快速时移估计方法，该方法包括以下步骤：

接收不同时刻的N个超声回波信号s₁,…,s_N；根据预先设定的窗长和重叠率，将所述N个超声回波信号s₁,…,s_N分为若干个信号片段，所述N个超声回波信号s₁,…,s_N上相同深度的信号片段分别表示为s_1k,…,s_Nk，所述s_nk＝s_n(k:M1+k-1)，其中信号片段s_nk为超声回波信号s_n的第k至第M1+k-1个采样点，n＝1,…,N，M1为窗长；计算所述信号片段s_2k相对所述超声回波信号s₁的时移d_2k，并以所述时移d_2k作为先验值计算所述信号片段s_3k相对所述超声回波信号s₁的时移d_3k，以此类推，直到以d_(n-1)k作为先验值计算所述信号片段s_nk和所述超声回波信号s₁的时移d_nk，其中n＝2,…,N。

进一步地，所述计算所述信号片段s_2k相对所述超声回波信号s₁的时移d_2k步骤包括：

分别计算所述信号片段s_2k与所述超声回波信号s₁(k+i:M+k+i-1)的绝对差和(SumofAbsoluteDifferences，简称SAD)值sad(i)，i＝-L,…,L，其中L是一个正整数，为时移估计的搜索范围，即

$sad\left(i\right)=\underset{q=k}{\overset{M1+k-1}{\Σ}}|s_2\left(q\right)-s_1(q+i)|,i=-L,...,L.$

计算d_2k＝arg(min(sad(i))，其中i＝-L,…,L，即找出sad(i)的最小值出现的位置d，所述信号片段s_2k相对s₁的时移d_2k＝d。

进一步地，所述以d_(n-1)k作为先验值计算所述信号片段s_nk和所述超声回波信号s₁的时移d_nk步骤包括：

分别计算s_nk与s₁(k+i:M+k+i-1)的SAD值sad(i)，i＝d_(n-1)k-L,…,d_(n-1)k+L，即

$sad\left(i\right)=\underset{q=k}{\overset{M1+k-1}{\Σ}}|s_n\left(q\right)-s_1(q+i)|,i=d_{(n-1)k}-L,...,d_{(n-1)k}+L.$

计算d_nk＝arg(min(sad(i))，其中i＝d_(n-1)k-L,…,d_(n-1)k+L，即找出sad(i)的最小值出现的位置d，所述信号片段s_nk相对s₁的时移d_nk＝d。

本发明可以有效地缩小时移估计的范围，提高计算速度。另外，因为直接估计了每次超声回波信号与第一次超声回波信号的时移，避免了误差累积。无需对信号进行插值，避免了计算量的增加。

附图说明

图1是本发明实施例提供的一种瞬时弹性成像快速时移估计方法流程图。

具体实施方式

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

图1为本发明实施例提供的一种瞬时弹性成像快速时移估计方法流程图。如图1所示，本发明实施例将接收不同时刻的N个超声回波信号s₁,…,s_N，根据预先设定的窗长和重叠率将N个超声回波信号s₁,…,s_N划分为若干信号片段，设信号s₁,…,s_N上相同深度的信号片段分别表示为s_1k,…,s_Nk，长度都为m1(即窗长为M1)，且s_nk＝s_n(k:M1+k-1)，即信号片段s_nk为信号s_n的第k至第M1+k-1个采样点，n＝1,…,N。

本发明实施例首先计算信号片段s_2k相对超声回波信号s₁的时移d_2k，并以时移d_2k作为先验值计算信号片段s_3k相对超声回波信号s₁的时移d_3k，以此类推，直到以d_(n-1)k作为先验值计算所述信号片段s_nk和所述超声回波信号s₁的时移d_nk，其中n＝2,…,N。

优选地，本发明实施例采用绝对差和(SumofAbsoluteDifferences，简称SAD)方法进行时移计算。

以下通过SAD方法计算信号片段s_2k相对超声回波信号s₁的时移d_2k，其实现步骤包括：

1)分别计算信号片段s_2k与超声回波信号s₁(k+i:M+k+i-1)的SAD值sad(i)，i＝-L,…,L，其中L是一个正整数，为时移估计的搜索范围，即

$sad\left(i\right)=\underset{q=k}{\overset{M1+k-1}{\Σ}}|s_2\left(q\right)-s_1(q+i)|,i=-L,...,L.$

2)计算d_2k＝arg(min(sad(i))，其中i＝-l,…,L，即找出sad(i)的最小值出现的位置d，则信号段s_2k相对s₁的时移d_2k＝d。

以下以d_(n-1)k作为先验值，用SAD方法估计信号段s_nk和超声回波信号s₁的时移d_nk，其实现步骤包括：

1)分别计算s_nk与s₁(k+i:M+k+i-1)的SAD值sad(i)，i＝d_(n-1)k-L,…,d_(n-1)k+L，即

$sad\left(i\right)=\underset{q=k}{\overset{M1+k-1}{\Σ}}|s_n\left(q\right)-s_1(q+i)|,i=d_{(n-1)k}-L,...,d_{(n-1)k}+L.$

2)计算d_nk＝arg(min(sad(i)),其中i＝d_(n-1)k-L,…,d_9n-10k+L，即找出sad(i)的最小值出现的位置d，则信号段s_nk相对s₁的时移d_nk＝d。

本发明实施例利用组织位移的连续性特征，以前一次时移估计值为先验值，每次估计时移的需要计算2L+1次SAD值，而如果不利用这个先验值，则搜索时移时需要计算SAD值的次数一定要远大于2L+1才能得到正确的时移估计。本发明实施例可以有效地缩小时移估计的范围，提高计算速度。另外，因为直接估计了每次超声回波信号与第一次超声回波信号的时移，避免了误差累积。无需对信号进行插值，避免了计算量的增加。

显而易见，在不偏离本发明的真实精神和范围的前提下，在此描述的本发明可以有许多变化。因此，所有对于本领域技术人员来说显而易见的改变，都应包括在本权利要求书所涵盖的范围之内。本发明所要求保护的范围仅由所述的权利要求书进行限定。

Claims (1)

1.一种瞬时弹性成像快速时移估计方法，其特征在于，包括以下步骤：

接收不同时刻的N个超声回波信号s₁、…、s_N；

根据预先设定的窗长和重叠率，将所述N个超声回波信号s₁、…、s_N分为若干个信号片段，所述N个超声回波信号s₁、…、s_N上相同深度的信号片段分别表示为s_1k、…、s_Nk，所述s_nk＝s_n(k:M1+k-1)，其中信号片段s_nk为超声回波信号信号s_n的第k至第M1+k-1个采样点，n＝1、…、N，M1为窗长；

计算所述信号片段s_2k相对所述超声回波信号s₁的时移d_2k，并以所述时移d_2k作为先验值计算所述信号片段s_3k相对所述超声回波信号s₁的时移d_3k，以此类推，直到以d_(n-1)k作为先验值计算所述信号片段s_nk和所述超声回波信号s₁的时移d_nk，其中n＝2、…、N；

所述计算所述信号片段s_2k相对所述超声回波信号s₁的时移d_2k步骤包括：

分别计算信号片段s_2k与超声回波信号s₁(k+i:M+k+i-1)的绝对差和SAD值sad(i)，i＝-L、…、L，其中L是一个正整数，为时移估计的搜索范围，即

$sad\left(i\right)=\underset{q=k}{\overset{M1+k-1}{\Σ}}|s_2\left(q\right)-s_1(q+i)|,i=-L,...,L$

计算d_2k＝arg(min(sad(i))，其中i＝-L、…、L，即找出sad(i)的最小值出现的位置d，所述信号片段s_2k相对s₁的时移d_2k＝d；

所述以d_(n-1)k作为先验值计算所述信号片段s_nk和所述超声回波信号s₁的时移d_nk步骤包括：

分别计算s_nk与s₁(k+i:M+k+i-1)的SAD值sad(i)，i＝d_(n-1)k-L、…、d_(n-1)k+L，即

$sad\left(i\right)=\underset{q=k}{\overset{M1+k-1}{\Σ}}|s_n\left(q\right)-s_1(q+i)|,i=d_{(n-1)k}-L,...,d_{(n-1)k}+L$

计算d_nk＝arg(min(sad(i))，其中i＝d_(n-1)k-L、…、d_(n-1)k+L，即找出sad(i)的最小值出现的位置d，所述信号片段s_nk相对s₁的时移d_nk＝d。