次奈奎斯特采样在超声波成像中的应用

林伟毅

（华南理工大学自动化科学与工程学院，广东广州510640）

近年来，Donoho和Candes等人提出了一种采样与压缩同步进行的理论——压缩感知理论（Compressive Sensing，CS）[1]，与传统的“先采样，后压缩”不同，CS理论是“边采样，边压缩”的方法。更重要的是，CS理论指出，信号的采样频率可以远远低于信号频谱中最高频率的2倍，而不影响信号的精确重建，只要信号是稀疏的，以低于奈奎斯特频率对信号进行采样的方式，即次奈奎斯特采样。将次奈奎斯特采样应用于超声波成像中，可以有效的减少采样点和采样频率，而这对于缩小超声波成像系统的仪器体积以及减少工作过程中所消耗的电能，都有着积极的意义。

FRI[2]模型是传统的采样理论与压缩感知相结合的采集信号的新方案，它最早由Vetterli等人提出来。FRI模型的提出是基于以下的设想的：在自然界中，很多信号都可以分解为一系列已知形状的短脉冲（假设脉冲数为L）。这样的信号，可以由各个短脉冲信号的延迟时间和幅度进行完备的表示。因此在特定的时间T内，只有有限的稀疏度（2L），文献[2]将稀疏度2L与时间T的比，定义为ROI（the Rate of Invovation），这样的信号叫做FRI（Finite Rate of Innovation）信号。ROI是理论上，FRI信号能够完美重建的最低采样频率，它远远低于奈奎斯特采样频率。超声波成像系统的反射信号可以看成由一系列不同延迟时间和幅度的脉冲信号的叠加。因此，可以利用FRI模型中的采样方法进行信号采样和重建。不幸的是，在之前的相关文献中，提出的一些对FRI信号进行采样的方案，比较著名的像B-Spline滤波器，E-Spline滤波器[3]，Intergrators滤波器[4]以及Exponential滤波器[5]，对于带噪声的高ROI信号（L≥10）的重建效果并不理想。所以，寻找一种适合高噪声、高ROI信号的次奈奎斯特采样方案，对于次奈奎斯特采样在超声波成像中的应用至关重要，本文的目的正在于此。

1 采样框架

假设我们得到了如下的一个持续时间为（0，T）的有限长FRI信号：

其中，h（t）是形状已知的脉冲并且满足关系式：

{tl，为未知的延迟时间和幅度，L为短脉冲的个数。假设我们得到了如下所示的M个傅里叶系数：

其中H（ω）为h（t）的连续傅里叶变换，κ是长度为M连续整数取值区间。令H为M×M的对角矩阵，其中第k个对角元素为令V（t）为的矩阵，它的第（k，j）元素为其中向量t={t1，t2，…，tL}，令向量a={a1，a2，…，aL}，令向量x表示傅里叶系数向量。那么（2）式可以重新写成矩阵形式：

我们可以选择区间，使得矩阵H是可逆的，定义y=H-1x，则有，

矩阵V是Vandermonde矩阵，因此只要满足条件2L≤M，V则是列满秩的。将方程组（5）进行重写，则有

由式（3），（4），（5）可知，一旦得到了傅里叶系数向量x，并且保证2L≤M，以及ti≠tj（i≠j），就可以根据方程组（6）解出，本文采用Annihilating Filter Method[3]进行信号重建，当然也可以使用其它的谱分析[6]。

下面将讨论一下如何在时间域上对FRI信号x（t）进行采样，以得到傅里叶系数向量x。

我们采用的是多通道混合采样的方式，使得每个通道的采样输出都是傅里叶系数向量x的线性组合，以防止通道失效，傅里叶系数丢失的问题。有限长FRI信号采样框架如图1所示。

图1 有限长FRI信号采样框架Fig.1 Sampling scheme of Finite stream of pulses

假设有P条通道，由图1可知，第i条通道的载波信号为：

其中，每一条通道中的k∈κ不尽相同，设采样向量为c={c1，c2，…，cP}，则有

令S为P×M的矩阵，它的第（i，k）元素为sik，那么（8）式可以改写为矩阵形式：

当且仅当S列满秩的时候，即P≥M时，S左可逆，则有：

S由我们所选择的载波信号si（t）所决定。

2 载波

我们采用方波pi（t）产生载波si（t），因为在实际电路中，方波的发生电路比较简单。方波的数学模型如下所示：

其中，ai[n]为长度是N的系数向量（N≥M），它的元素只能取±1，p（t）为：

它的连续傅里叶变换为：

下面，我们将会介绍：以滤波的方式，从方波信号pi（t）得到实际载波信号si（t），以及推导出这种情况下，混合矩阵S的计算公式。

1）由方波信号pi（t）产生si（t）

由式（14）可以看出，pi（t）是一个周期函数。我们知道，周期函数可以由傅里叶系数如下表示：

其中，傅里叶系数由di[k]由以下公式计算得出：

可以看出，与式（7）相比，式（14）是无限相加的形式，因此，为了产生si（t），必须加一个滤波器g（t），对pi（t）进行滤波

g（t）的连续傅里叶变换必须有如下限制：

那么，可得，

由式（7）和式（17）可知（18），下面推导S时将会用到：

2）混合矩阵S的计算公式

由式（11）、（15）计算可得：

结合（19）和（20），可得:

其中k′=-（k-■M/2」）。那么混合矩阵S可以表示为：

A为P×N的矩阵，第（i，n）元素Ain=ai[n]；W为N×M的矩

结合式（13）、（17），（23）可以改写为：

3 MATLAB仿真

1）利用超声波工业成像系统采集一组一维原始信号，其中探头的中心频率为fc=3.423 5 MHz，系统的采样频率fs=40 MHz，成像深度为Rmax=0.16 m，在被测物体中，超声波速C=1 540，那么可以计算出信号的持续时间为T=2×Rmax/C=2.08×10-4sec，因此可以计算出，若用传统采样方法，每次成像的数据点数目为SamplingNo_classic=T×fs=8 320，我们的仿真，只采用了一半数据，因此T=T/2=1.04×10-4sec，SamplingNo_classic=SamplingNo_classic/2=4 160；

2）根据超声波工业探头的特性，可知h（t）为高斯脉冲是比较合适的，其中高斯参数，σ=3×10-7

3）利用图1所示的框架对原始信号进行混合积分采样，设L=4，P=N=M=2×L×oversampling+1，由于原始信号噪声很大，因此设oversampling=5，则：P=N=M=4，κ={-20，-19，…，0，…，19，20}；

4）我们选择方波信号作为载波，其中方波信号符合式（11）、（12）、方波的个数为41，周期为T，宽度为T/N，滤波器g（t）符合式（17），该滤波器为低通滤波器，带宽为2π/T；

5）积分采样后，得到向量c，根据方波信号写出A矩阵，根据式（22）、（24）得到S矩阵；

6）利用公式（10）计算出傅里叶系数向量x，利用公式，y=H-1x，求出向量y，从而得到方程组（6）；

7）利用Annihilating Filter Method求解方程（6），得到参数，重建效果，如图2所示，其中信号大小和时间都归1了。

图2 信号重建效果图（虚线部分）Fig.2 Applying our method on real ultrasound imaging data.Results are shown vs.full demodulated signal

由图2可知，延迟时间的重建十分准确（其它的波峰波谷为噪声），而幅值的重建则有偏差，尤其是第二个高斯脉冲。但这并不会影响到超声波工业成像的效果，因为超声波工业成像是为了找寻缺陷的位置，而延迟时间包含了位置的全部信息。比较传统的香农定理采样与我们的采样方法可知，我们方法的采样频率约是实际采样频率的0.024%（1/T/fs），采样数据点数目约是传统采样方法的0.8%（41/4 160）。

4 结束语

文中提出了一种新型的多通道FRI信号的采集方法，详细地介绍了它的信号采样过程以及信号重建过程，丰富了次奈奎斯特采样方式；为次奈奎斯特采样在超声波成像系统中的应用迈出了坚实的一步，并为以后的研究工作指明方向。

[1] Candes E，Wakin M.An introduction to compressive sampling[J].IEEE Signal Process Magazine，2008，25（2）:21-30.

[2] Vetterli M，Marziliano P，Blu T.Sampling signals with?nite rate of innovation[J].IEEE Transactions on Signal Processing，2002，50（6）:1417-1428.

[3] Dragotti P L，Vetterli M，Blu T.Sampling moments and reconstructing signals of fnite rate of innovation:Shannon meets strang-fix[J].IEEE Transactions on Signal Processing，2007，55（5）:1741-1757.

[4] Kusuma J，Goyal V.Multichannel sampling of parametric signals with a successive approximation property[J].IEEE Int.Conf.Image Process，2006，50(6）:1265-1268.

[5] Olkkonen H，Olkkonen J.Measurement and reconstruction of impulse train by parallel exponential filters[J].IEEE Signal Process.Lett.，2008，15（1）:241-244.

[6] Stoica P，Moses R.Introduction to Spectral Analysis[M].3rd ed.Englewood Cliffs，NJ:Prentice-Hall，1997.