超声回波信号模型

简介

超声回波信号模型(Echo Model, EM)是对超声回波信号经验性的建模。选择适当的模型，可以用少量参数表示一个回波信号。虽然是经验性的，部分模型参数仍与具有实际意义的物理参数相关联，例如与飞行时间(Time of Flight, ToF)一致的时延参数。超声回波信号建模、模型参数估计方法以及模型参数的物理解释等研究从上世纪末至今一直较为活跃。

时域模型的发展

多数应用中，超声时域信号直接由模数转换器(Analog-to-Digital Converter, ADC)转换而来。时域模型即是对超声应用中最容易获取的时域信号的建模。
R. Demirli相继提出了高斯回波模型和不对称高斯回波模型，郭纲提出了双指数回波模型，卢振坤提出了混合指数回波模型。

分类

时域模型

从信号角度看，超声回波的时域信号是一种时域有限信号，意味着回波信号仅出现在时域部分区域。截取一段含有单个超声回波的信号，可以使用包络和振荡项的乘积来近似表示该信号。这就是时域回波模型的基本思想：
$$ s(t)=env(t)osc(t) $$
接下来分别从包络和振荡项两方面介绍常见的时域回波模型。

包络(Envelope)

高斯(Gaussian)
高斯包络是最基本、最常见的包络形式，表达式为：
$$ env_\mathrm{g}(t)=\exp(-\frac{t^2}{\alpha^2}) $$
式中$\alpha$常称为带宽系数。

grid();
axis([0, 0], [0, 1]);
axis([0, 0], [1, 0]);

a = slider([-2, 1.2], [2, 1.2], [0.01, 1, 2]) <<name: 'α'>>;
f = function (x) {return exp(-x*x/V(a)/V(a));};
g = functiongraph(f) <<withLabel: true, name: 'env_g(t)'>>;

$board.setView([-5, 1.5, 5, -0.3]);

不对称高斯(Asymmetric Gaussian)
不对称高斯包络是在高斯包络基础上，使其对称轴左右两侧的$\alpha$不同的一种包络。最直接的实现可以为：
$$
env_{\mathrm{agd1}}(t)=\exp(-\left[\frac{t^2}{\alpha_l^2}u(-t)+\frac{t^2}{\alpha_r^2}u(t)\right])
,\ u(t)=\begin{cases}1,&t>0 \\ 0,&t\le0 \end{cases}
$$
或
$$
env_{\mathrm{agd2}}(t)=\exp(-\frac{t^2}{\alpha^2}(1+\mathrm{sgn}(t)r))
,\ r\in(-1,1)
,\ \mathrm{sgn}(t)=\begin{cases}1,&t>0 \\ 0,&t=0 \\ -1,&t<0 \end{cases} \quad
$$
这类直接(direct)实现的缺陷在于包络在$t=0$时不可导。因此采用平滑化的方法($\tanh(\kappa t)$是对$\mathrm{sgn}(t)$很好的一种平滑化近似)，得到：
$$
env_{\mathrm{ag}}(t)=\exp(-\frac{t^2}{\alpha^2}(1+\tanh(\kappa t)r))
,\ r\in(-1,1)
,\ \kappa > 0
$$
式中$r$用于控制不对称程度，$\kappa$用于控制在$t=0$附近的平滑程度。先来观察这个平滑化的$\mathrm{sgn}(t)$：
```
grid();
axis([0, 0], [0, 1]);
axis([0, 0], [1, 0]);

k = slider([-4, 3.2], [4, 3.2], [0.1, 1, 10]) <<name: 'κ'>>;
r = slider([-4, 2.5], [4, 2.5], [-0.99, 0.2, 0.99]) <<name: 'r'>>;
f1 = function (x) {return tanh(V(k)*x);};
f2 = function (x) {return 1+V(r)*tanh(V(k)*x);};
g1 = functiongraph(f1) <<withLabel: true, name: 'f(t)=tanh(κt)'>>;
g2 = functiongraph(f2) <<withLabel: true, name: 'f(t)=1+tanh(κt)r'>>;

$board.setView([-10, 3.5, 10, -1.5]);
```
整个包络的形状则如下图所示：
```
grid();
axis([0, 0], [0, 1]);
axis([0, 0], [1, 0]);

a = slider([-4, 1.6], [4, 1.6], [0.2, 3, 10]) <<name: 'α'>>;
k = slider([-4, 1.35], [4, 1.35], [0.1, 1, 2]) <<name: 'κ'>>;
r = slider([-4, 1.1], [4, 1.1], [-0.99, -0.5, 0.99]) <<name: 'r'>>;
f = function (x) {return exp(-x*x/V(a)/V(a)*(1+V(r)*tanh(V(k)*x)));};
g = functiongraph(f) <<withLabel: true, name: 'env_{ag}(t)'>>;

$board.setView([-10, 1.8, 10, -0.3]);
```
混合指数(Mixed Exponential)
混合指数包络用于描述上升快下降慢的包络，其形式为：
$$ env_{\mathrm{me}}(t)=t^m\exp(-\frac{t}{T}),\ m\in\mathbb{N}_+,\ T>0 $$
可以看出，与高斯包络的最大区别是控制变化率的项由指数中的平方项，变为了外部的幂函数乘数项。求导数后不难得到$t=mT$时，函数有极大值$ env_{me}(mT)=m^mT^me^{-m}$。为了方便可视化，考虑将包络归一化：
$$ env_{\mathrm{men}}(t)=\left(\frac{t}{mT}\right)^m\exp(m-\frac{t}{T}),\ m\in\mathbb{N}_+,\ T>0 $$
```
grid();
axis([0, 0], [0, 1]);
axis([0, 0], [1, 0]);

m = slider([2, 1.6], [14, 1.6], [1, 2, 8]) <<name: 'm', snapWidth: 1>>;
T = slider([2, 1.35], [14, 1.35], [0.1, 1.4, 2.5]) <<name: 'T'>>;

f = function (x) {return (x/V(m)/V(T))^V(m)*exp(V(m)-x/V(T));};
g = functiongraph(f, 0, 30) <<withLabel: true, name: 'env_{men}(t)'>>;

$board.setView([-1, 1.8, 20, -0.3]);
```
双指数(Dual Exponential)
双指数包络同样用于描述上升快下降慢的包络，其形式为：
$$ env_{\mathrm{de}}(t)=e^{-\alpha t}-e^{-\beta t},\ \alpha>0,\ \beta>0 $$
在郭纲的论文中，比较了混合指数包络和双指数包络，认为双指数包络可以实现更快的上升速度，更适用于文中的电子白板定位信号的处理。类似的，可以求得该包络的极值点：
$$t_M=\frac{\ln\alpha-\ln\beta}{\alpha-\beta}$$
以及极值：
$$
\begin{align}
env_{\mathrm{de}}(t_M)&=\exp(-\alpha\frac{\ln\alpha-\ln\beta}{\alpha-\beta}) - \exp(-\beta\frac{\ln\alpha-\ln\beta}{\alpha-\beta}) \\
&=(\frac{\alpha}{\beta})^{-\alpha/(\alpha-\beta)}-(\frac{\alpha}{\beta})^{-\beta/(\alpha-\beta)}
\end{align}
$$
所以同样为了方便可视化，可以归一化如下：
$$
env_{\mathrm{den}}(t)=env_{\mathrm{de}}(t)/env_{\mathrm{de}}(t_M)
$$
```
grid();
axis([0, 0], [0, 1]);
axis([0, 0], [1, 0]);

a = slider([2, 1.6], [7, 1.6], [0.2, 3, 10]) <<name: 'α'>>;
b = slider([2, 1.35], [7, 1.35], [0.2, 5, 10]) <<name: 'β'>>;
pole = function () {return ln(V(a)/V(b))/(V(a)-V(b));};
f = function (x) {return exp(-x*V(a))-exp(-x*V(b));};
fn = function (x) {return f(x)/f(pole());};
g = functiongraph(fn, 0, 30) <<withLabel: true, name: 'env_{den}(t)'>>;

$board.setView([-1, 1.8, 10, -0.3]);
```

小结
以上共总结了4种用于回波信号模型的包络，其实除了这些包络之外，还有其他形式的包络，例如这篇文献中的拉盖尔函数包络。上述的4种包络中，双指数包络和混合指数包络与高斯型包络(包括对称与非对称)有很大差别，主要体现在极点上。由于极点与包络参数直接相关，双指数包络和混合指数包络的适用性并不如高斯型包络广泛。

振荡项(Oscillator)

正弦(Sinudoidal)
因为和常用的超声激励信号相同，正弦形式的振荡项是很自然就想到的一种振荡形式：
$$
osc_{\sin}(t)=\cos(2\pi f_ct+\phi)
$$
式中共有两个参数，$f_c$为信号频率，$\phi$为初始相位。
线性调频(Chirp)
定义瞬时频率(Instaneous Frequency, IF)如下：
$$
f_I(t)=\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}
$$
在正弦振荡项中，总相位$\Phi(t)=2\pi f_ct+\phi$，故正弦信号的瞬时频率为：
$$
f_{I,\sin}(t)=\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=2\pi f_c
$$
为时间无关的常量。如果信号频率随时间线性变化，即得Chirp振荡项，有：
$$
f_{I,\mathrm{chirp}}(t)=\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=2\pi f_c+kt
$$
故可得：
$$
\Phi_\mathrm{chirp}(t) = \int (f_c+\hat{k}t)dt = 2\pi f_ct + kt^2 + \phi
$$
所以Chirp振荡项可以表达如下：
$$
osc_\mathrm{chirp}(t) = \cos(\Phi_\mathrm{chirp}(t)) = \cos(2\pi f_ct + kt^2 + \phi)
$$
对称线性调频(Symmetric Chirp, SChirp)¹
在Chirp中，信号的瞬时频率总是单调的(递增或递减)，如果希望信号瞬时频率以某个中心点为对称轴，两侧递增或递减时，可以设计SChirp如下:
$$
\begin{align}
f_{I,\mathrm{schirp}}(t) &= \frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=2\pi f_c+k|t| \\
\Phi_\mathrm{schirp}(t) &= \int (f_c+\hat{k}|t|)dt = 2\pi f_ct + kt^2\mathrm{sgn}(t) + \phi \\
osc_\mathrm{schirp}(t) &= \cos(2\pi f_ct + kt^2\mathrm{sgn}(t) + \phi)
\end{align}
$$
```
grid();
axis([0, 0], [0, 1]);
axis([0, 0], [1, 0]);

k = slider([0.5, 1.6], [3, 1.6], [-1.5, 0, 1.5]) <<name: 'k'>>;
fc = slider([-4.5, 1.6], [-2.5, 1.6], [0, 2, 4]) <<name: 'f_c'>>;
p = slider([-4.5, 1.3], [-2.5, 1.3], [0, 0, 2*PI]) <<name: 'φ'>>;
f = function (x) {return cos(2*PI*V(fc)*x+V(k)*x*x*SGN(x)+V(p));};
g = functiongraph(f);

$board.setView([-5, 1.8, 5, -1.3]);
```
需要注意，$k<0$时，SChirp的瞬时频率会在降为0后递增。

完整的回波模型

包络项和振荡项任意相乘即可得到完整的回波模型。例如最基本，也是最常用的高斯回波模型：
$$
s_\mathrm{gauss}(t)=\beta env_\mathrm{g}(t)osc_\mathrm{sin}(t)=\beta \exp(-\frac{t^2}{\alpha^2})\cos(2\pi f_ct + \phi)
$$
在基于模型(Model Based)的超声信号处理中，使用延时项来匹配实际回波信号，可以表示为：
$$
s_\mathrm{GEM}(\theta;t)=s_\mathrm{gauss}(t-\tau)=\beta \exp(-\frac{(t-\tau)^2}{\alpha^2})\cos(2\pi f_c(t-\tau) + \phi),\ \theta=[\alpha,\beta,f_c,\tau,\phi]
$$

grid();
axis([0, 0], [0, 1]);
axis([0, 0], [1, 0]);

a = slider([-4.5, 5], [-2, 5], [0.05, 1, 4]) <<name: 'α'>>;
b = slider([-4.5, 4.5], [-2, 4.5], [0.1, 2, 4]) <<name: 'β'>>;
t = slider([0.5, 5], [3, 5], [-2, 0, 2]) <<name: 'τ'>>;
fc = slider([0.5, 4.5], [3, 4.5], [0.1, 2, 5]) <<name: 'f_c'>>;
p = slider([0.5, 4], [3, 4], [0, 0, 2*PI]) <<name: 'φ'>>;
f = function (x) {xm=x-V(t); return V(b)*exp(-xm*xm/V(a)/V(a))*cos(2*PI*V(fc)*xm+V(p));};
g = functiongraph(f);

$board.setView([-5, 5.5, 5, -4]);

也可以形成较为复杂的模型，例如AGCM:
$$
\begin{align}
s_\mathrm{AGCM}(\theta;t)&=\beta env_\mathrm{ag}(t-\tau)osc_\mathrm{chirp}(t-\tau) \\
&=\beta\exp(-\frac{t^2}{\alpha^2}[1+\tanh(\kappa (t-\tau))r])\cos(2\pi f_ct + k(t-\tau)^2 + \phi) \\
\theta&=[\alpha,\beta,f_c,\tau,\phi,\kappa,r,k],\ r\in(-1,1)
\end{align}
$$

grid();
axis([0, 0], [0, 1]);
axis([0, 0], [1, 0]);

a = slider([-4.5, 5.5], [-2, 5.5], [0.05, 1, 4]) <<name: 'α'>>;
b = slider([-4.5, 5], [-2, 5], [0.1, 2, 4]) <<name: 'β'>>;
t = slider([0.5, 5.5], [3, 5.5], [-2, 0, 2]) <<name: 'τ'>>;
fc = slider([0.5, 5], [3, 5], [0.1, 2, 5]) <<name: 'f_c'>>;
p = slider([0.5, 4.5], [3, 4.5], [0, 0, 2*PI]) <<name: 'φ'>>;
kappa = slider([0.5, 4], [3, 4], [0.1, 1, 2]) <<name: 'κ'>>;
k = slider([-4.5, 4.5], [-2, 4.5], [-2, 0, 2]) <<name: 'k'>>;
r = slider([-4.5, 4], [-2, 4], [-0.4, 0, 0.4]) <<name: 'r'>>;
f = function (x) {
    xm=x-V(t); 
    return V(b)*exp(-xm*xm/V(a)/V(a)*(1+tanh(V(kappa)*xm*V(r))))*cos(2*PI*V(fc)*xm+V(k)*xm*xm+V(p));
};
g = functiongraph(f);

$board.setView([-5, 6, 5, -4]);

个人想法，暂无文献支持。 ↩︎