SoS 濾波器及其在工業(yè)超聲波成像中的應(yīng)用*

戴光智 韓國(guó)強(qiáng) 林偉毅 歐陽顯躍

(1.華南理工大學(xué) 計(jì)算機(jī)科學(xué)與工程學(xué)院，廣東 廣州 510006;2.華南理工大學(xué) 自動(dòng)化科學(xué)與工程學(xué)院，廣東 廣州 510640;3.湖南大學(xué) 嵌入式與網(wǎng)絡(luò)計(jì)算湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410082)

Donoho 和Candès 等［1-4］提出的采樣與壓縮同步進(jìn)行的壓縮感知(CS)理論與傳統(tǒng)的“先采樣后壓縮”不同，它采用“邊采樣邊壓縮”的方法.CS 理論指出，只要信號(hào)是稀疏的，那么信號(hào)的采樣頻率可以遠(yuǎn)遠(yuǎn)低于信號(hào)頻譜中最高頻率的2 倍而不影響信號(hào)的精確重建.將CS 理論應(yīng)用于工業(yè)超聲波成像系統(tǒng)中，可以有效地減少采樣點(diǎn)和采樣頻率.目前，國(guó)內(nèi)外將CS 應(yīng)用于超聲領(lǐng)域的研究［5-13］主要集中在醫(yī)用超聲方面，在工業(yè)超聲方面的研究很少.有限新息率(FRI)模型是傳統(tǒng)采樣理論與CS 相結(jié)合的信號(hào)采集方法［14］.具有FRI 性質(zhì)的信號(hào)可以由各個(gè)短脈沖信號(hào)的延時(shí)和幅度進(jìn)行完備的表示，因此可以有效地減少采樣數(shù)據(jù)和采樣頻率.工業(yè)超聲波成像系統(tǒng)的反射信號(hào)可以看成一系列不同延時(shí)和幅度的高斯脈沖信號(hào)的疊加，因此，工業(yè)超聲波成像系統(tǒng)的反射信號(hào)具有FRI 性質(zhì)，可以利用FRI 模型中的采樣方法進(jìn)行信號(hào)采樣和重建.但現(xiàn)有的FRI 信號(hào)采集方法對(duì)高頻信號(hào)進(jìn)行采樣時(shí)很不穩(wěn)定.通常，工業(yè)超聲波成像系統(tǒng)的反射信號(hào)頻率很高，導(dǎo)致現(xiàn)有的FRI 模型并不能很好地應(yīng)用于工業(yè)超聲波成像系統(tǒng)中.因此，文中研究能夠采集高頻信號(hào)的FRI 信號(hào)采集方法，對(duì)CS 理論在工業(yè)超聲波成像系統(tǒng)中的應(yīng)用具有重大的意義.

為有效地減少采樣數(shù)據(jù)和采樣頻率，文中提出了一種適用于工業(yè)超聲信號(hào)處理的新型FRI 采集方法——辛格函數(shù)之和(SoS)濾波器法，推導(dǎo)了這種采集方法應(yīng)具備的條件，給出了SoS 濾波器的構(gòu)造方式，并利用SoS 濾波器對(duì)一維超聲波成像原始數(shù)據(jù)進(jìn)行處理，最后通過仿真實(shí)驗(yàn)驗(yàn)證了該方法的可行性.

1 壓縮感知理論和FRI 模型

1.1 壓縮感知理論

CS 理論指出，當(dāng)信號(hào)滿足稀疏或可壓縮性時(shí)，可以采集少量觀測(cè)信號(hào)，然后通過求解一個(gè)優(yōu)化問題準(zhǔn)確重建原始信號(hào).由信號(hào)理論可知，時(shí)域信號(hào)X(XRn×1)能夠用一組基Ψ(Ψ=(ψ1，ψ2，…，ψn))線性表示，即

式中，α 為n 維系數(shù)向量，α=(α1，α2，…，αn).如果式(1)中系數(shù)αi僅有k 個(gè)非0，則稱信號(hào)X 是稀疏的，稱Ψ 為信號(hào)X 的稀疏基或稀疏字典.如果對(duì)X的測(cè)量也是在時(shí)域上進(jìn)行的，則有

其中Θ(Θ=ΦΨ)為m×n 階矩陣.由于式(3)中的α 是k 稀疏的，且k ＜m ＜n，故α 可以通過式(4)的最優(yōu)化問題獲得，即

可以通過基于l1范數(shù)的最優(yōu)化方法(如貝葉斯方法)求解目標(biāo)信號(hào)α，但其前提是測(cè)量矩陣必須滿足RIP 性質(zhì).

1.2 FRI 模型

假設(shè)得到如下一個(gè)周期為 的無限周期信號(hào):

式中，h(t)為形狀和周期已知的脈沖，tl和al(l=1，2，…，L)分別為未知的延時(shí)和幅度，L 為一個(gè)周期內(nèi)短脈沖的個(gè)數(shù).假設(shè)一個(gè)周期內(nèi)得到M 個(gè)連續(xù)的傅里葉系數(shù)，則x(t)的傅里葉變換為

選擇M個(gè)連續(xù)的整數(shù)，使得矩陣H 是可逆的，定義y=H-1x ，則有

式中，V為Vandermonde 矩陣，因此只要L ≤M，則V 是列滿秩的.

將方程組(8)進(jìn)行重寫，有

由式(7)-(9)可知，一旦得到傅里葉系數(shù)向量x，并保證2L≤M，就可以根據(jù)方程組(8)解出tl和al(l= 1，2，…，L).文中采用零化濾波重建算法［15］求解方程組(8)，當(dāng)然也可以使用其他的方法，如矩陣束算法［16］.

2 FRI 濾波器組及SoS 濾波器

2.1 FRI 濾波器組

傅里葉系數(shù)的采集是稀疏重建成像的關(guān)鍵.在時(shí)間域上對(duì)FRI 信號(hào)x(t)進(jìn)行采樣(采樣函數(shù)為s*())，可以得到M(2L≤M)個(gè)傅里葉系數(shù)，能實(shí)現(xiàn)此功能的方法的集合稱為FRI 濾波器組.假設(shè)待采樣信號(hào)為周期信號(hào)，其采樣框架如圖1 所示.

圖1 單通道采樣框架Fig.1 Framework of single channel sampling

由圖1 可知，對(duì)于FRI 信號(hào)x(t)及采樣周期T，實(shí)際采集到的向量為c= (c［1］，c［2］，…，c［N］)，其中

式中:n=1，2，…，N;S*(ω)為S(ω)的伴隨矩陣的元素，S(ω)為s(t)的連續(xù)傅里葉變換，ω=2 k/ .實(shí)際上，只要選擇

可見，對(duì)濾波后得到的采樣向量c 進(jìn)行離散傅里葉變換DFT{c}，再乘以S-1，就能得到傅里葉系數(shù)x.

在文獻(xiàn)［5］的FRI 濾波器模型中，s*(- t)=Bsin(-cBt)，B=M/ ，2L ≤M ≤N，并且s(t)是一個(gè)帶寬為B 的理想低通濾波器:

對(duì)于一個(gè)有限信號(hào)，采用FRI 濾波器的采樣模型定義如下:

假設(shè)h(t)滿足條件

經(jīng)過推算可以驗(yàn)證，當(dāng)R ≤ 時(shí)，有以下公式成立:

其中，s3p(t)=s(t - )+ s(t)+ s(t + ).

2.2 SoS 濾波器

由于理想低通濾波器帶寬是有限的，而脈沖信號(hào)在時(shí)域上是無限的，因此文獻(xiàn)［5］的FRI 濾波器具有局限性，它不能擴(kuò)展到無限脈沖信號(hào).通常，工業(yè)超聲回波信號(hào)的頻率很高，但超聲探頭及轉(zhuǎn)換硬件性能的限制使所獲得的超聲信號(hào)帶寬遠(yuǎn)遠(yuǎn)低于實(shí)際的信號(hào)帶寬，故存在較大的信息丟失.為此，基于FRI 濾波器組，在文獻(xiàn)［6］的基礎(chǔ)之上，文中提出了帶寬不受限制的SoS 濾波器:

式中，bk≠0，bk=，

將它轉(zhuǎn)化到時(shí)域上，有

同樣，有限信號(hào)采用SoS 濾波器采樣模型，定義為

選擇正確的參數(shù)bk(kκ)，能夠增強(qiáng)濾波器的抗噪性能，文中取長(zhǎng)度為M 的漢明窗系數(shù):

可知，bk=，因此濾波后的采樣系數(shù)c［n］必定為實(shí)數(shù).

3 實(shí)驗(yàn)及結(jié)果分析

3.1 Matlab 仿真結(jié)果分析

SoS 濾波法的實(shí)現(xiàn)步驟如下:

(1)輸入包含L 個(gè)高斯脈沖的信號(hào)x(t)，高斯脈沖為

其中σ=7 ×10-3，周期= 1 s，時(shí)間延時(shí)和幅度是隨機(jī)選擇的;

(2)根據(jù)實(shí)際的超聲檢測(cè)，取脈沖個(gè)數(shù)L=5，κ={- L，- L +1，…，L -1，L}，M== 11 ;

(3)采用漢明窗系數(shù)滿足式(23)(計(jì)算結(jié)果見圖2)的g(t)對(duì)輸入信號(hào)進(jìn)行濾波;

圖2 bk 的選取Fig.2 Selection of bk

根據(jù)式(7)-(9)對(duì)圖3 所示信號(hào)進(jìn)行濾波，SoS濾波器的實(shí)部如圖4 所示，采樣結(jié)果如圖5 所示，原始信號(hào)的估計(jì)延時(shí)和幅度如圖6 所示，所得到的重建信號(hào)如圖7 所示，由圖可見重建效果很好.

圖3 包含不同周期的原始信號(hào)Fig.3 Original signals containing different periods

圖4 SoS 濾波器的實(shí)部Fig.4 Real component of SoS filter

圖5 濾波輸出信號(hào)中采樣的21 個(gè)值Fig.5 Twenty-one sampling values of filter output signal

圖6 原始信號(hào)的估計(jì)延時(shí)和幅度Fig.6 Estimation time delay and amplitude of original signal

圖7 仿真重建信號(hào)與原始信號(hào)對(duì)比Fig.7 Comparison of reconstruction signal with original signal in simulation

3.2 實(shí)測(cè)結(jié)果分析

(1)利用工業(yè)超聲波成像系統(tǒng)采集一組一維原始信號(hào)(如圖8 所示)，其中單一探頭的中心頻率fc=3.4235 MHz，系統(tǒng)的采樣頻率fs= 40 MHz，成像深度Rmax= 0.16 m，在被測(cè)物體中，超聲波速c=1540 m/s，可以計(jì)算出信號(hào)的持續(xù)時(shí)間 '=2Rmax/c=2.08×10-4s;

(2)根據(jù)超聲波工業(yè)探頭的特性，將h(t)定為高斯脈沖是比較合適的，選擇高斯參數(shù)σ=3 ×10-7s，因σ ≤ '，故可以使用濾波器(21)對(duì)采樣數(shù)據(jù)進(jìn)行濾波，結(jié)果見圖9;

(3)取L=5，為了增強(qiáng)抗噪性能，設(shè)M=2 ×5 ×2+1=21，N=M=21，κ= {- 10，- 9，…，0，…，9，10}，采樣周期T= '/N=2.08 ×10-4/21 s=9.91 μs;

(4)利用式(20)-(22)構(gòu)造濾波器，計(jì)算采樣向量c;

(5)利用式(14)計(jì)算傅里葉系數(shù)向量x;

(6)利用y=H-1x 求出向量y，從而得到方程組(9);

(7)利用零化濾波重建算法求解方程(9)，得到估計(jì)的延時(shí)和幅度如圖10 所示，重建信號(hào)見圖11.

圖8 實(shí)測(cè)超聲信號(hào)Fig.8 Actually measured ultrasonic signals

圖9 濾波輸出信號(hào)的采樣值Fig.9 Sampling values of filter output signal

圖10 原始超聲信號(hào)的估計(jì)延時(shí)和幅度Fig.10 Estimation time delay and amplitude of original ultrasonic signal

圖11 實(shí)測(cè)原始超聲信號(hào)和重建信號(hào)比較Fig.11 Comparison of original ultrasonic signal with reconstruction signal in actual measurment

由重建效果可以看出，延時(shí)系數(shù)是精確的.在超聲波成像系統(tǒng)中，延遲系數(shù)代表著缺陷的位置，因此，只要延時(shí)系數(shù)精確，就能準(zhǔn)確地確定缺陷的位置.

4 結(jié)語

根據(jù)工業(yè)超聲波信號(hào)頻率較高的特點(diǎn)，基于FRI 濾波器組，文中提出了一種新的適用于工業(yè)超聲信號(hào)處理的FRI 信號(hào)采集方法——SoS 濾波器法，Matlab 仿真和實(shí)測(cè)結(jié)果表明該方法是有效的，從而驗(yàn)證了CS 理論在工業(yè)超聲波檢測(cè)系統(tǒng)中應(yīng)用的可行性.文中就CS 理論在工業(yè)超聲波檢測(cè)系統(tǒng)中的應(yīng)用做了有益的探討和實(shí)測(cè)數(shù)據(jù)仿真，進(jìn)一步豐富了CS 理論.

［1］Donoho D L.Compressed sensing［J］.IEEE Transactions on Information Theory，2006，52(4):1289-1306.

［2］Candès E J，Romberg J，Tao T.Robust uncertainty principles:exact signal reconstruction from highly incomplete frequency information［J］.IEEE Transactions on Information Theory，2006，52(2):489-509.

［3］Candès E J.Compressive sampling［C］∥Proceedings of the International Congress of Mathematics.Madrid:European Mathematical Society，2006:1433-1452.

［4］Candès E J，Wakin M B.An introduction to compressive sampling［J］.IEEE Signal Processing Magazine，2008，25(2):21-30.

［5］Tur R，Eldar Y C，F(xiàn)riedman Z.Innovation rate sampling of pulse streams with application to ultrasound imaging［J］.IEEE Transactions on Signal Processing，2011，59(4):1827-1842.

［6］Wagner N，Eldar Y C，F(xiàn)euer A，et al.Compressed beamforming in ultrasound imaging［J］.IEEE Transactions on Signal Processing，2012，60(9):4643-4657.

［7］Chernyakova T，Eldar Y C.Fourier domain beamforming:the path to compressed ultrasound imaging［J/OL］.Arxiv，2013:11(2013-07-24).http:∥arxiv.org/abs/1307.6345.

［8］Wagner N，Eldar Y C，F(xiàn)euer A，et al.Compressed beamforming with applications to ultrasound imaging ［C］∥Proceedings of 2012 IEEE International Conference on Acoustics，Speech and Signal Processing.Kyoto:IEEE，2012:3641-3644.

［9］Friboulet D，Liebgott H，Prost R.Compressive sensing for raw RF signals reconstruction in ultrasound ［C］∥Proceedings of 2010 IEEE Ultrasonics Symposium.San Diego:IEEE，2010:367-370.

［10］Quinsac C，Basarab A，Kouamé D.Frequency domain compressive sampling for ultrasound imaging［J/OL］.Advances in Acoustics and Vibration，2012，2012:16(2012-04-05).http:∥dx.doi.org/10.1155/2012/231317.

［11］Quinsac C，Basarab A，Kouame D，et al.3D compressed sensing ultrasound imaging［C］∥Proceedings of 2010 IEEE Ultrasonics Symposium.San Diego:IEEE，2010:363-366.

［12］Quinsac C，Basarab A，Girault J，et al.Compressed sensing of ultrasound images:sampling of spatial and frequency domains ［C］∥Proceedings of 2010 IEEE Workshop on Signal Processing Systems.San Diego:IEEE，2010:231-236.

［13］Shi Guangming，Chen Chongyu，Lin Jie，et al.Narrowband ultrasonic detection with high range resolution:separating echoes via compressed sensing and singular value decomposition［J］.IEEE Transactions on Ultrasonics，F(xiàn)erroelectrics and Frequency Control，2012，59(10):2237-2253.

［14］Vetterli M，Marziliano P，Blu T.Sampling signals with finite rate of innovation［J］.IEEE Transactions on Signal Processing，2002，50(6):1417-1428.

［15］Stoica P，Moses R.Introduction to spectral analysis［M］.Englewood Cliffs:Prentice-Hall，2000.

［16］Sarkar T K，Pereira O.Using the matrix pencil method to estimate the parameters of a sum of complex exponentials［J］.IEEE Antennas and Propagation Magazine，1995，37(1):48-55.