卡爾曼濾波解卷積算法的優(yōu)化問題*

閆曉瑾，何國榮，梅鐵民

（沈陽理工大學(xué)自動化與電氣工程學(xué)院，遼寧 沈陽 110159）

0 引言

信號解卷積是信號處理領(lǐng)域中一個重要的研究內(nèi)容，其廣泛應(yīng)用于圖像處理、語音信號處理、陣列信號處理、地震波信號處理、生物醫(yī)學(xué)、故障診斷等領(lǐng)域。在這些領(lǐng)域中，通常觀測信號是源信號與一線性系統(tǒng)沖激響應(yīng)的卷積，解卷積就是從觀測信號中恢復(fù)源信號或系統(tǒng)沖激響應(yīng)的過程。

已有的解卷積算法大致可分為兩類，一類是房間沖激響應(yīng)未知情況下的盲解卷積算法，另一類是房間沖激響應(yīng)已知的解卷積算法。盲解卷積是指在房間沖激響應(yīng)未知或不需要估計(jì)的情況下，由觀測信號直接估計(jì)解卷積濾波器。Bussgang 算法[1]是自適應(yīng)信號處理領(lǐng)域中最經(jīng)典的盲解卷積算法，該算法實(shí)現(xiàn)簡單且收斂速度快，但該算法需要源信號和解卷積噪聲的統(tǒng)計(jì)特性等先驗(yàn)知識。針 對SISO（Single Input Single Output）系 統(tǒng)，文獻(xiàn)[2]將觀測信號轉(zhuǎn)化為復(fù)基帶信號，然后利用最小均方算法（Least Mean Square,LMS）來實(shí)現(xiàn)盲解卷積，該算法適用于高斯信號和非高斯信號，對源信號的適用面較寬，但仍然存在收斂速度和計(jì)算復(fù)雜度之間的矛盾問題。在生物醫(yī)學(xué)中，心電圖信號、脈搏信號等可以被認(rèn)為是心跳信號與人體系統(tǒng)的卷積，通過對心電圖信號進(jìn)行解卷積處理即可獲得病人的心跳信號，文獻(xiàn)[3]中利用遞歸最小二乘算法（Recursive least squares,RLS）來求解卷積濾波器，利用定點(diǎn)卷積核補(bǔ)償（Fixed-point convolution kernel compensation,FP-CKC）求出RLS算法的期望信號，利用粒子群優(yōu)化算法（Particle Swarm optimization,PSO）求出逆濾波器的初始值,該算法可用于較高的噪聲水平。在語音去混響領(lǐng)域中，解卷積過程是其中重要的一部分。文獻(xiàn)[4]根據(jù)聲源語音信號的短時傅里葉變換系數(shù)的稀疏性，將非負(fù)矩陣分解（Non-negative Matrix Factorization,NMF）應(yīng)用到多通道線性預(yù)測（Multi-Channel Linear Prediction,MCLP）中來達(dá)到去混響的目的，該算法取得不錯的效果，不足之處在于明顯增加了計(jì)算量。在軸承故障診斷領(lǐng)域中，振動信號是由故障引起的周期性沖擊信號與機(jī)械部件的共振響應(yīng)卷積的結(jié)果，傳統(tǒng)的解卷積方法有最大相關(guān)峭度解卷積算法（Maximum Correlated Kurtosis Deconvolution,MCKD），該算法的解卷積信號的包絡(luò)譜中可以呈現(xiàn)出典型的周期性脈沖的特點(diǎn)[5]。文獻(xiàn)[6]提出最小熵解卷積算法（Mimum-Entropy Deconvolution，MED），該算法可適用于非最小相位系統(tǒng)，但要求源信號具有簡單稀疏的統(tǒng)計(jì)特性。除此之外，還有基于稀疏表示的盲解卷積算法[7]、子空間分解的盲解卷積方法[8]等。在語音去混響解卷積算法中，通常先估計(jì)出房間沖激響應(yīng)，然后再進(jìn)行解卷積。在SISO 系統(tǒng)中，文獻(xiàn)[9]利用復(fù)倒譜技術(shù)估計(jì)出房間沖激響應(yīng)，并在倒譜域內(nèi)求其逆濾波器，然后將其逆濾波器作為改進(jìn)的自然梯度算法的初始值來進(jìn)行迭代解卷積。MINT（the multiple-input/output inverse-filtering theorem）定理是Miyoshi 等人提出的在房間沖激響應(yīng)已知的情況下求逆濾波器的一種經(jīng)典方法[10]，文獻(xiàn)[11]和[12]利用觀測信號的自相關(guān)矩陣來獲得房間沖激響應(yīng)的先驗(yàn)知識，然后根據(jù)MINT 定理來求逆濾波器，但該算法要求源信號是平穩(wěn)的白噪聲信號。其實(shí)，利用MINT 定理求逆濾波器的方法在實(shí)際中并不可行，因?yàn)樵摲椒▽υ肼暦浅Ｃ舾?，針對此問題，文獻(xiàn)[13]提出在子帶最小二乘算法中使用正則化，該算法可以降低逆濾波器對房間沖激響應(yīng)的估計(jì)誤差和對觀測噪聲的敏感度。文獻(xiàn)[14]通過觀測信號的二階累積量估計(jì)出非最小相位系統(tǒng)的房間沖激響應(yīng)，對其卷積矩陣直接求逆，該方法在高階系統(tǒng)時的計(jì)算量將非常大。文獻(xiàn)[15]在估計(jì)出房間沖激響應(yīng)后，利用p-范數(shù)和窗函數(shù)來求逆濾波器，進(jìn)而實(shí)現(xiàn)解卷積。文獻(xiàn)[16]給出了在估計(jì)出房間沖激響應(yīng)的情況下不求逆濾波器的解卷積算法，利用卡爾曼濾波算法直接進(jìn)行解卷積，該算法具有較高的噪聲穩(wěn)定性,但其缺點(diǎn)是計(jì)算量大。

本文在卡爾曼濾波解卷積算法的基礎(chǔ)上，提出了一種降低計(jì)算量的優(yōu)化算法，該算法在保證算法穩(wěn)定性的同時，又能大幅度減少計(jì)算量。該算法不僅可以用于語音去混響，還可以應(yīng)用在其他領(lǐng)域。

1 標(biāo)準(zhǔn)的卡爾曼濾波解卷積算法

1.1 數(shù)學(xué)模型

在單輸入多輸出系統(tǒng)（Single Input Multiple Output,SIMO）中，觀測信號是聲源信號和房間沖激響應(yīng)的卷積，那么，觀測信號可以表示為：

其中，s(n)表示聲源信號，n為時間，hi表示第i路房間沖激響應(yīng)，xi(n)表示第i路麥克風(fēng)接收的觀測信號，vi(n)表示第i路觀測噪聲，L為房間沖激響應(yīng)的長度，N為麥克風(fēng)的數(shù)量。

利用（1）式建立卡爾曼濾波的狀態(tài)方程和觀測方程：

狀態(tài)方程：

測量方程：

其中S(n)為n時刻的狀態(tài)矢量，由聲源信號s(n)構(gòu)成狀態(tài)矢量S(n)=[s(n),s(n-1),…,s(n-L+1)]T；X(n)=[x1(n),x2(n),…,xN(n)]T為n時刻的觀測矢量；觀測矩陣由N路房間沖激響應(yīng)構(gòu)成，是一個N×L維矩陣（其中hi=[hi(0),hi(1),…,hi(L-1)]T）；u1(n)和u2(n)分別為均值為零的過程白噪聲和觀測白噪聲，它們的協(xié)方差矩陣分別為和（I為單位矩陣）；W(n+1,n)為L×L維狀態(tài)轉(zhuǎn)移矩陣，定義為：

并且，W(n,n+1)=WT(n+1,n)。

1.2 標(biāo)準(zhǔn)的卡爾曼濾波解卷積算法

已知房間沖激響應(yīng)（h1,h2,...hN），在一個解卷積塊中的標(biāo)準(zhǔn)卡爾曼濾波解卷積算法的迭代過程列于表1 中。

表1 標(biāo)準(zhǔn)卡爾曼濾波解卷積算法

其中G(n)為卡爾曼增益，P(n)為濾波的狀態(tài)矢量估計(jì)誤差的自相關(guān)矩陣，K(n+1,n)為狀態(tài)矢量一步預(yù)測誤差的自相關(guān)矩陣，α(n)為新息?？柭鼮V波器輸出的n時刻的聲源信號的最后一個分量。

標(biāo)準(zhǔn)卡爾曼濾波算法中五個變量迭代一次的計(jì)算量統(tǒng)計(jì)于表2 中。關(guān)于算法的計(jì)算量問題，本文只統(tǒng)計(jì)乘除法的運(yùn)算次數(shù)，不統(tǒng)計(jì)加減法。從表2中可以看出，標(biāo)準(zhǔn)卡爾曼濾波算法的計(jì)算量主要集中在G(n)、P(n)和K(n+1,n)上，相比之下，和α(n)的計(jì)算量可以忽略不計(jì)，所以要想減少計(jì)算量，需從G(n)、P(n)和K(n+1,n)入手。

表2 標(biāo)準(zhǔn)卡爾曼濾波算法的計(jì)算量

2 改進(jìn)的卡爾曼濾波解卷積算法

仔細(xì)研究標(biāo)準(zhǔn)卡爾曼濾波解卷積算法發(fā)現(xiàn)，G(n)、P(n)和K(n+1,n)三者的迭代計(jì)算在極短的時間內(nèi)即可達(dá)到收斂，并且G(n)只與H(n)有關(guān)，與觀測信號無關(guān)；而α(n)、除與H(n)有關(guān)外，還與當(dāng)前時刻的觀測信號X(n)有關(guān)。對于非時變系統(tǒng)而言，G(n)一旦收斂，就不再變化，因此可以把與G(n)計(jì)算有關(guān)的項(xiàng)提取出來單獨(dú)計(jì)算，待其收斂后再進(jìn)行α(n)和的計(jì)算，從而實(shí)現(xiàn)解卷積。因此，可以通過把標(biāo)準(zhǔn)卡爾曼濾波解卷積算法分解成兩個子循環(huán)的辦法來到達(dá)降低樣本平均計(jì)算量的目的。新算法的迭代步驟列于表3 中。

表3 改進(jìn)的卡爾曼濾波解卷積算法

通過改變標(biāo)準(zhǔn)算法的迭代次數(shù)和迭代順序，使其在保證穩(wěn)定性的同時，又能明顯減少計(jì)算量。在卡爾曼增益G(n)的計(jì)算中，Bg的選取要具體問題具體分析（例如，若算法要求計(jì)算量小，則Bg盡量??；若算法要求收斂性，則Bg可以大一些）。該算法通過犧牲收斂速度來降低算法的計(jì)算量，但犧牲的收斂速度有限。

在分塊解卷積算法中，認(rèn)為在一個塊內(nèi)系統(tǒng)沖激響應(yīng)是不變的，假設(shè)一個解卷積塊的長度為B，則從上述迭代過程中可以看出，改進(jìn)算法比標(biāo)準(zhǔn)算法減少的計(jì)算量為(B-Bg)(4L3+4L2N+2LN2+N3)。

3 仿真實(shí)驗(yàn)

在仿真實(shí)驗(yàn)中用到的數(shù)據(jù)是實(shí)測的房間沖激響應(yīng)（信道數(shù)N為4，長度為1300 點(diǎn)）和一段語音信號（長度為21000 點(diǎn)）。把已知的語音信號作為聲源信號s(n)，然后將聲源信號與實(shí)測的房間沖激響應(yīng)卷積（在不同的仿真實(shí)驗(yàn)中，房間沖激響應(yīng)的長度會有多不同）來獲得卷積信號。由于本文只討論算法計(jì)算量的問題，所以在仿真實(shí)驗(yàn)中不加入噪聲。

令解卷積信號估計(jì)誤差e(n)=se(n)-s(n)的均方值為E(n)；聲源信號的均方值為Q(n)。采用如下的滑動平均公式在線估計(jì)：

其中，滑動因子α=0.9；n是迭代次數(shù)；e(n)為算法輸出誤差；se(n)為算法的解卷積信號；s(n)為聲源信號。

定義解卷積信號的誤差-信號比作為評價算法性能的標(biāo)準(zhǔn)：

Esir曲線能夠反映出算法的收斂水平，單位是分貝（dB）。

3.1 卡爾曼增益的收斂速度

在第2 節(jié)中提到，對于非時變系統(tǒng)而言，G(n)一旦收斂，就不再變化，而卡爾曼增益G(n)的計(jì)算只與房間沖激響應(yīng)H(n)有關(guān)，所以本節(jié)的仿真內(nèi)容是確定卡爾曼增益的收斂速度。在本節(jié)中房間沖激響應(yīng)的長度分別取L=100，150，200（從已知的房間沖激響應(yīng)中任意截?。?。標(biāo)準(zhǔn)卡爾曼濾波算法中對于不同長度的房間沖激響應(yīng)時卡爾曼增益的收斂情況見圖1。由于卡爾曼增益在全局的收斂速度較快，所以圖1 中的卡爾曼增益的收斂情況為前400點(diǎn)的放大圖。從圖1 中可以看出，卡爾曼增益的收斂速度幾乎與房間沖激響應(yīng)的長度L相同。

3.2 改進(jìn)的卡爾曼濾波算法對解卷積結(jié)果的影響

在本節(jié)的仿真實(shí)驗(yàn)中，房間沖激響應(yīng)長度L=150，在上一節(jié)的仿真實(shí)驗(yàn)中發(fā)現(xiàn)卡爾曼增益的收斂速度幾乎與房間沖激響應(yīng)的長度L相同，所以本節(jié)中選取改進(jìn)算法的卡爾曼增益的迭代次數(shù)Bg=L+10=160。圖2（a）給出了仿真卷積信號（1）聲源信號（2）改進(jìn)算法的解卷積信號（3）和標(biāo)準(zhǔn)算法的解卷積信號（4）的波形比較。其中仿真卷積信號是第2 路的卷積信號。從圖2 中可以看出兩種算法都有較好的解卷積效果，改進(jìn)算法只是收斂速度慢一點(diǎn)，一旦收斂，二者的解卷積結(jié)果一致。

圖1 不同長度的房間沖激響應(yīng)對應(yīng)的卡爾曼增益的收斂情況（L=100，150，200）

圖2 聲源信號、卷積信號和兩種算法的解卷積信號的比較

3.3 改進(jìn)的卡爾曼濾波算法對收斂速度的影響

本節(jié)與上一節(jié)實(shí)驗(yàn)相同，房間沖激響應(yīng)長度L=150，B=5L=750，Bg=L+10=160。利用兩種算法的解卷積信號來繪制出Esir 曲線，如圖3 所示。

圖3 標(biāo)準(zhǔn)算法解卷積信號和改進(jìn)算法解卷積信號的Esir 曲線

從圖3 可以看出，雖然改進(jìn)算法的解卷積信號沒有標(biāo)準(zhǔn)算法解卷積信號那么高的信噪比，但是改進(jìn)算法的解卷積信號的Esir 也能夠快速達(dá)到50dB以下。從算法的運(yùn)行時間來看，在仿真實(shí)驗(yàn)所用的計(jì)算機(jī)上，本節(jié)實(shí)驗(yàn)的標(biāo)準(zhǔn)算法的運(yùn)行時間是5.92秒，改進(jìn)算法的運(yùn)行時間是1.56 秒，減少了73.6%的運(yùn)行時間，說明計(jì)算量減少了70%以上。

4 結(jié)語

本文提出了一種基于卡爾曼濾波解卷積的優(yōu)化算法，利用卡爾曼增益在短時間內(nèi)收斂的特點(diǎn)，將標(biāo)準(zhǔn)卡爾曼濾波算法用兩個子循環(huán)來完成，卡爾曼增益、狀態(tài)矢量一步預(yù)測誤差的自相關(guān)矩陣和狀態(tài)估計(jì)誤差的自相關(guān)矩陣的迭代作為一個子循環(huán)，狀態(tài)矢量和新息的迭代作為一個子循環(huán)。在仿真實(shí)驗(yàn)中可以看出，與標(biāo)準(zhǔn)卡爾曼濾波解卷積算法相比，本文提出的改進(jìn)算法犧牲的收斂速度有限，但降低的計(jì)算量是明顯的。