次奈奎斯特采樣在超聲波成像中的應(yīng)用

林偉毅

（華南理工大學(xué)自動(dòng)化科學(xué)與工程學(xué)院，廣東廣州510640）

近年來(lái)，Donoho和Candes等人提出了一種采樣與壓縮同步進(jìn)行的理論——壓縮感知理論（Compressive Sensing，CS）[1]，與傳統(tǒng)的“先采樣，后壓縮”不同，CS理論是“邊采樣，邊壓縮”的方法。更重要的是，CS理論指出，信號(hào)的采樣頻率可以遠(yuǎn)遠(yuǎn)低于信號(hào)頻譜中最高頻率的2倍，而不影響信號(hào)的精確重建，只要信號(hào)是稀疏的，以低于奈奎斯特頻率對(duì)信號(hào)進(jìn)行采樣的方式，即次奈奎斯特采樣。將次奈奎斯特采樣應(yīng)用于超聲波成像中，可以有效的減少采樣點(diǎn)和采樣頻率，而這對(duì)于縮小超聲波成像系統(tǒng)的儀器體積以及減少工作過(guò)程中所消耗的電能，都有著積極的意義。

FRI[2]模型是傳統(tǒng)的采樣理論與壓縮感知相結(jié)合的采集信號(hào)的新方案，它最早由Vetterli等人提出來(lái)。FRI模型的提出是基于以下的設(shè)想的：在自然界中，很多信號(hào)都可以分解為一系列已知形狀的短脈沖（假設(shè)脈沖數(shù)為L(zhǎng)）。這樣的信號(hào)，可以由各個(gè)短脈沖信號(hào)的延遲時(shí)間和幅度進(jìn)行完備的表示。因此在特定的時(shí)間T內(nèi)，只有有限的稀疏度（2L），文獻(xiàn)[2]將稀疏度2L與時(shí)間T的比，定義為ROI（the Rate of Invovation），這樣的信號(hào)叫做FRI（Finite Rate of Innovation）信號(hào)。ROI是理論上，F(xiàn)RI信號(hào)能夠完美重建的最低采樣頻率，它遠(yuǎn)遠(yuǎn)低于奈奎斯特采樣頻率。超聲波成像系統(tǒng)的反射信號(hào)可以看成由一系列不同延遲時(shí)間和幅度的脈沖信號(hào)的疊加。因此，可以利用FRI模型中的采樣方法進(jìn)行信號(hào)采樣和重建。不幸的是，在之前的相關(guān)文獻(xiàn)中，提出的一些對(duì)FRI信號(hào)進(jìn)行采樣的方案，比較著名的像B-Spline濾波器，E-Spline濾波器[3]，Intergrators濾波器[4]以及Exponential濾波器[5]，對(duì)于帶噪聲的高ROI信號(hào)（L≥10）的重建效果并不理想。所以，尋找一種適合高噪聲、高ROI信號(hào)的次奈奎斯特采樣方案，對(duì)于次奈奎斯特采樣在超聲波成像中的應(yīng)用至關(guān)重要，本文的目的正在于此。

1 采樣框架

假設(shè)我們得到了如下的一個(gè)持續(xù)時(shí)間為（0，T）的有限長(zhǎng)FRI信號(hào)：

其中，h（t）是形狀已知的脈沖并且滿足關(guān)系式：

{tl，為未知的延遲時(shí)間和幅度，L為短脈沖的個(gè)數(shù)。假設(shè)我們得到了如下所示的M個(gè)傅里葉系數(shù)：

其中H（ω）為h（t）的連續(xù)傅里葉變換，κ是長(zhǎng)度為M連續(xù)整數(shù)取值區(qū)間。令H為M×M的對(duì)角矩陣，其中第k個(gè)對(duì)角元素為令V（t）為的矩陣，它的第（k，j）元素為其中向量t={t1，t2，…，tL}，令向量a={a1，a2，…，aL}，令向量x表示傅里葉系數(shù)向量。那么（2）式可以重新寫(xiě)成矩陣形式：

我們可以選擇區(qū)間，使得矩陣H是可逆的，定義y=H-1x，則有，

矩陣V是Vandermonde矩陣，因此只要滿足條件2L≤M，V則是列滿秩的。將方程組（5）進(jìn)行重寫(xiě)，則有

由式（3），（4），（5）可知，一旦得到了傅里葉系數(shù)向量x，并且保證2L≤M，以及ti≠tj（i≠j），就可以根據(jù)方程組（6）解出，本文采用Annihilating Filter Method[3]進(jìn)行信號(hào)重建，當(dāng)然也可以使用其它的譜分析[6]。

下面將討論一下如何在時(shí)間域上對(duì)FRI信號(hào)x（t）進(jìn)行采樣，以得到傅里葉系數(shù)向量x。

我們采用的是多通道混合采樣的方式，使得每個(gè)通道的采樣輸出都是傅里葉系數(shù)向量x的線性組合，以防止通道失效，傅里葉系數(shù)丟失的問(wèn)題。有限長(zhǎng)FRI信號(hào)采樣框架如圖1所示。

圖1 有限長(zhǎng)FRI信號(hào)采樣框架Fig.1 Sampling scheme of Finite stream of pulses

假設(shè)有P條通道，由圖1可知，第i條通道的載波信號(hào)為：

其中，每一條通道中的k∈κ不盡相同，設(shè)采樣向量為c={c1，c2，…，cP}，則有

令S為P×M的矩陣，它的第（i，k）元素為sik，那么（8）式可以改寫(xiě)為矩陣形式：

當(dāng)且僅當(dāng)S列滿秩的時(shí)候，即P≥M時(shí)，S左可逆，則有：

S由我們所選擇的載波信號(hào)si（t）所決定。

2 載波

我們采用方波pi（t）產(chǎn)生載波si（t），因?yàn)樵趯?shí)際電路中，方波的發(fā)生電路比較簡(jiǎn)單。方波的數(shù)學(xué)模型如下所示：

其中，ai[n]為長(zhǎng)度是N的系數(shù)向量（N≥M），它的元素只能取±1，p（t）為：

它的連續(xù)傅里葉變換為：

下面，我們將會(huì)介紹：以濾波的方式，從方波信號(hào)pi（t）得到實(shí)際載波信號(hào)si（t），以及推導(dǎo)出這種情況下，混合矩陣S的計(jì)算公式。

1）由方波信號(hào)pi（t）產(chǎn)生si（t）

由式（14）可以看出，pi（t）是一個(gè)周期函數(shù)。我們知道，周期函數(shù)可以由傅里葉系數(shù)如下表示：

其中，傅里葉系數(shù)由di[k]由以下公式計(jì)算得出：

可以看出，與式（7）相比，式（14）是無(wú)限相加的形式，因此，為了產(chǎn)生si（t），必須加一個(gè)濾波器g（t），對(duì)pi（t）進(jìn)行濾波

g（t）的連續(xù)傅里葉變換必須有如下限制：

那么，可得，

由式（7）和式（17）可知（18），下面推導(dǎo)S時(shí)將會(huì)用到：

2）混合矩陣S的計(jì)算公式

由式（11）、（15）計(jì)算可得：

結(jié)合（19）和（20），可得:

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

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

結(jié)合式（13）、（17），（23）可以改寫(xiě)為：

3 MATLAB仿真

1）利用超聲波工業(yè)成像系統(tǒng)采集一組一維原始信號(hào)，其中探頭的中心頻率為fc=3.423 5 MHz，系統(tǒng)的采樣頻率fs=40 MHz，成像深度為Rmax=0.16 m，在被測(cè)物體中，超聲波速C=1 540，那么可以計(jì)算出信號(hào)的持續(xù)時(shí)間為T=2×Rmax/C=2.08×10-4sec，因此可以計(jì)算出，若用傳統(tǒng)采樣方法，每次成像的數(shù)據(jù)點(diǎn)數(shù)目為SamplingNo_classic=T×fs=8 320，我們的仿真，只采用了一半數(shù)據(jù)，因此T=T/2=1.04×10-4sec，SamplingNo_classic=SamplingNo_classic/2=4 160；

2）根據(jù)超聲波工業(yè)探頭的特性，可知h（t）為高斯脈沖是比較合適的，其中高斯參數(shù)，σ=3×10-7

3）利用圖1所示的框架對(duì)原始信號(hào)進(jìn)行混合積分采樣，設(shè)L=4，P=N=M=2×L×oversampling+1，由于原始信號(hào)噪聲很大，因此設(shè)oversampling=5，則：P=N=M=4，κ={-20，-19，…，0，…，19，20}；

4）我們選擇方波信號(hào)作為載波，其中方波信號(hào)符合式（11）、（12）、方波的個(gè)數(shù)為41，周期為T，寬度為T/N，濾波器g（t）符合式（17），該濾波器為低通濾波器，帶寬為2π/T；

5）積分采樣后，得到向量c，根據(jù)方波信號(hào)寫(xiě)出A矩陣，根據(jù)式（22）、（24）得到S矩陣；

6）利用公式（10）計(jì)算出傅里葉系數(shù)向量x，利用公式，y=H-1x，求出向量y，從而得到方程組（6）；

7）利用Annihilating Filter Method求解方程（6），得到參數(shù)，重建效果，如圖2所示，其中信號(hào)大小和時(shí)間都?xì)w1了。

圖2 信號(hào)重建效果圖（虛線部分）Fig.2 Applying our method on real ultrasound imaging data.Results are shown vs.full demodulated signal

由圖2可知，延遲時(shí)間的重建十分準(zhǔn)確（其它的波峰波谷為噪聲），而幅值的重建則有偏差，尤其是第二個(gè)高斯脈沖。但這并不會(huì)影響到超聲波工業(yè)成像的效果，因?yàn)槌暡üI(yè)成像是為了找尋缺陷的位置，而延遲時(shí)間包含了位置的全部信息。比較傳統(tǒng)的香農(nóng)定理采樣與我們的采樣方法可知，我們方法的采樣頻率約是實(shí)際采樣頻率的0.024%（1/T/fs），采樣數(shù)據(jù)點(diǎn)數(shù)目約是傳統(tǒng)采樣方法的0.8%（41/4 160）。

4 結(jié)束語(yǔ)

文中提出了一種新型的多通道FRI信號(hào)的采集方法，詳細(xì)地介紹了它的信號(hào)采樣過(guò)程以及信號(hào)重建過(guò)程，豐富了次奈奎斯特采樣方式；為次奈奎斯特采樣在超聲波成像系統(tǒng)中的應(yīng)用邁出了堅(jiān)實(shí)的一步，并為以后的研究工作指明方向。

[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.