基于快速平行因子分解的聲矢量傳感器陣二維DOA 估計(jì)

陳未央，徐 樂(lè)，張小飛

（南京航空航天大學(xué)電子信息工程學(xué)院，南京211106）

陣列信號(hào)處理中的接收信號(hào)波達(dá)方向（Direc?tion of arrival，DOA）估計(jì)，近年來(lái)獲得了廣泛的關(guān)注和迅速的發(fā)展［1?3］，而在水聲信號(hào)處理中，新型的聲矢量傳感器也越來(lái)越多地被應(yīng)用到信號(hào)參數(shù)估計(jì)中。聲矢量傳感器一般由聲壓傳感器和振速傳感器組成，相比于傳統(tǒng)的信號(hào)接收傳感器，聲矢量傳感器可以同時(shí)測(cè)量空間中某點(diǎn)聲場(chǎng)的標(biāo)量和矢量，即同步共點(diǎn)測(cè)量聲場(chǎng)中一點(diǎn)處的聲壓和質(zhì)點(diǎn)振速若干正交分量，具備較強(qiáng)的抗各向同性噪聲干擾的能力［4］，這些特點(diǎn)使得聲矢量傳感器陣列對(duì)于信號(hào)參數(shù)估計(jì)具有較大優(yōu)勢(shì)。

近年來(lái)，學(xué)者們已相繼提出了許多關(guān)于聲矢量傳 感 器 陣 列DOA 估 計(jì) 的 算 法［5?11］，包 括 多 重 分 類（Multiple signal classification，MUSIC）算法［5?6］、借助旋轉(zhuǎn)不變性進(jìn)行信號(hào)參數(shù)估計(jì)（Estimation of signal parameters via rotational invariance tech?nique，ESPRIT）算 法［7?8］、傳 播 算 子（Propagator method，PM）［9］，以及最大似然（Maximum likeli?hood，ML）算法［10］等。其中，基于子空間類的方法，如MUSIC 算法和ESPRIT 算法，通常要對(duì)接收信號(hào)協(xié)方差矩陣進(jìn)行特征值分解，其計(jì)算復(fù)雜度較高。而PM 算法利用傳播算子矩陣獲得信號(hào)的參數(shù)估計(jì)，避免了接收信號(hào)協(xié)方差矩陣分解過(guò)程，其復(fù)雜度較低，但在低信噪比的情況下，該算法參數(shù)估計(jì)性能較差。

平行因子（Parallel factor，PARAFAC）模型最早于生理學(xué)中提出，用于對(duì)多維數(shù)據(jù)分析［12?13］，近年來(lái)，該技術(shù)被成功地運(yùn)用到了信號(hào)處理領(lǐng)域。文獻(xiàn)［14］利用平行因子模型中的PARAFAC 分解方法，能夠得到聲矢量傳感器陣列二維DOA 估計(jì)。然而傳統(tǒng)PARAFAC 算法直接利用隨機(jī)矩陣進(jìn)行初始化，如采用高斯隨機(jī)矩陣初始化，后進(jìn)行循環(huán)迭代更新，直到分解收斂。但該過(guò)程所需收斂次數(shù)較多，迭代分解過(guò)程收斂較慢，導(dǎo)致算法計(jì)算復(fù)雜度很高，特別是在陣列陣元數(shù)較大或信號(hào)采樣快拍數(shù)較大的情況下。

為了降低平行因子方法的復(fù)雜度，本文提出了一種適合任意矢量傳感器陣列結(jié)構(gòu)的基于快速PARAFAC 分解的二維DOA 估計(jì)算法。該算法首先將接收信號(hào)構(gòu)建為平行因子模型，然后在數(shù)據(jù)域?qū)υ撃Ｐ椭械膮?shù)矩陣進(jìn)行初始估計(jì)，最后利用三線性交替最小二乘（Trilinear alternating least square，TALS）算法［15?16］得到和信號(hào)一一匹配的仰角和方位角估計(jì)。所提算法基于減少迭代分解的次數(shù)、加快收斂速度的思想，對(duì)傳統(tǒng)平行因子方法改進(jìn)，對(duì)所構(gòu)建的平行因子模型中的參數(shù)矩陣進(jìn)行了初始化估計(jì)，大大加快了PARAFAC 分解的迭代過(guò)程。與普通的PARAFAC 算法相比，該算法復(fù)雜度較低。同時(shí)，其二維DOA 估計(jì)性能非常接近于PARAFAC 算法，同時(shí)優(yōu)于ESPRIT 算法和PM 算法。

1 數(shù)據(jù)模型

假設(shè)有K 個(gè)遠(yuǎn)場(chǎng)信號(hào)被陣列所接收，其中信源數(shù)K 已知，陣列由M 個(gè)任意分布在三維空間中的聲矢量傳感器組成，其結(jié)構(gòu)如圖1 所示。

圖1 信號(hào)接收陣列結(jié)構(gòu)圖Fig.1 Structure of signal receiving array

假設(shè)陣列所接收的信號(hào)均為非相干的窄帶信號(hào)，所含噪聲均為加性高斯白噪聲，且與信號(hào)獨(dú)立［14］。第k 個(gè)信號(hào)的波達(dá)方向?yàn)?φk，θk)，其中φk為信號(hào)方位角，θk為信號(hào)仰角。則陣列中第m 個(gè)陣元在t 時(shí)刻的接收信號(hào)可以表示為［14］

式中：sk(t)為第k 個(gè)信源的傳輸信號(hào)；τm為相較于參考陣元，信號(hào)到達(dá)第m 個(gè)傳感器的時(shí)延；nm(t)為接收噪聲；hk為第k 個(gè)信源的方位矢量，表示信號(hào)在各個(gè)方向上的正交分量，可寫為［4］

整個(gè)陣列在t 時(shí)刻的接收信號(hào)可以寫為［14］

式中：“°”表示矩陣的Khatri?Rao 積；A∈CM×K為整個(gè)陣列的方向矩陣［14］；s(t)=[s1(t)，s2(t)，…，sK(t)]T為信源矢量；H =[h1，h2，…，hK]∈C4×K。則陣列J個(gè)快拍的接收信號(hào)可以表示為

式中：Dm(·)表示取矩陣的第m 行組成的對(duì)角陣，SJ×K=[s(t1)，…，s(tJ)]T為 信 源 矩 陣；Nx=[n(t1)，…，n(tJ)]為J 個(gè) 快 拍 的 信 號(hào) 噪 聲；Xm=HDm(A)ST+Nxm(m=1，2，…，M )，其 中Nxm為切片噪聲。Xm可以表示成平行因子三線性模型［15］，即有

式中：am，k和hp，k分別為矩陣A 和H 的第（m，k）和（p，k）元素，sj，k為S 的第（j，k）元素，nm，p，j為噪聲切片Nxm的 第（p，j）元 素；m=1，2，…，M；p=1，2，3，4；j=1，2，…，J。該結(jié)構(gòu)還可以寫成另外兩種表達(dá)形式［14］

式中Ny和Nz為變換后相應(yīng)的噪聲。

2 二維DOA 估計(jì)算法

普通的PARAFAC算法通常直接利用三線性交替最小二乘算法對(duì)該模型進(jìn)行分解。但該迭代分解過(guò)程收斂較慢，導(dǎo)致PARAFAC 算法計(jì)算復(fù)雜度較高。為了加快接收信號(hào)模型分解收斂速度，本文算法先在數(shù)據(jù)域?qū)δＰ椭械膮?shù)矩陣進(jìn)行初估計(jì)，然后再利用PARAFAC分解得到信號(hào)的二維DOA估計(jì)。

2.1 參數(shù)矩陣初估計(jì)

式（4）中的接收信號(hào)，經(jīng)過(guò)初等行變換可以表示為

式中：N′為變換后相應(yīng)的噪聲；G ∈C4M×4M為行變換矩陣，表示為

定義B=H °A，對(duì)矩陣B 進(jìn)行分塊

式中：B1∈CK×K；B2∈C(4M-K)×K。則存在一個(gè)矩陣P ∈CK×(4M-K)滿足［9］

接收信號(hào)的協(xié)方差矩陣為

對(duì)R 進(jìn)行分塊得到

式中：R1∈C4M×K；R2∈C4M×(4M-K)。在無(wú)噪聲情況下R1和R2滿 足

可以通過(guò)求解下式的最小值來(lái)獲得矩陣P 的估計(jì)［9］

式中‖ · ‖F(xiàn)表 示矩陣 的Frobenius 范數(shù)，通過(guò) 式（16）可以求得P 的估計(jì)為

得到P?的估計(jì)后，可以得到矩陣P?c=[ IK×K，P?]T。對(duì)P?c進(jìn)行分塊

式 中P1、P2、P3和P4均 為M ×K 的 矩 陣。則 根 據(jù)PcB1=B 以及B=H °A，有

根據(jù)式（19，20）可以得到P2B1=P1B1D2(H )，然后可得

式中：(·)?表示矩陣廣義逆，對(duì)P?1P2進(jìn)行特征值分解，其特征值即為D2(H )的估計(jì)值D2(H?)。同時(shí)，可 以 得 到B1的 估 計(jì) 值B?1，進(jìn) 而 可 得B 的 估 計(jì)B?=P?cB?1。將B?進(jìn)行分塊可以得到

式 中B1、B2、B3和B4均 為M ×K 的 矩 陣。由 式（22）即可得到矩陣A 的初始估計(jì)A?ini=B1。同時(shí)可以得到B3=B1D3(H )，B4=B1D4(H )，然后可以得到D3(H )和D4(H )的估計(jì)為

至此，已經(jīng)獲得了矩陣A 的估計(jì)A?ini。同時(shí)，D2(H?)、D3(H )和D4(H )的對(duì)角元素即分別為矩陣H 第2、3、4 行的估計(jì)，由于H 第1 行全部由元素1 組成，因此最終也可獲得H 的初始估計(jì)H?ini。

2.2 三線性分解與二維DOA 估計(jì)

三線性交替最小二乘算法是平行因子模型中常用的分解方法，該算法每次更新模型中的一個(gè)矩陣的估計(jì)，直到收斂。

通過(guò)式（4）可以得到最小二乘擬合為［14］

由式（25）可以得到S 的估計(jì)為

式中A?和H?為上一次迭代所得到的A 和H 的估計(jì)值。

類似地，式（6）的最小二乘擬合為

由式（27）可以得到A 的估計(jì)為

式中H?和S?為上一次的迭代所得到的H 和S 的估計(jì)值。

同理，通過(guò)式（7）可以得到H?的估計(jì)為

式中S?和A?為上一次迭代所得到的S 和A 的估計(jì)值。

定理［16］考 慮 平 行 因 子 模 型 Xm=HDm(A)ST+Nxm(m=1，2，…，M )， 其 中 ，H ∈C4×K、A∈CM×K、S ∈CJ×K，3 個(gè) 矩 陣 的k 秩［15］分別為kH、kA和kS。如果

則S、A 和H 對(duì)于列交換和尺度變換是唯一的。

最終可以獲得矩陣H 的估計(jì)為

式中：Π 為列模糊矩陣，Λ 為尺度因子，W 為估計(jì)誤差。列模糊問(wèn)題對(duì)角度的估計(jì)沒(méi)有影響，尺度模糊則可以通過(guò)歸一化方法解決。

對(duì)H?進(jìn) 行 歸 一 化 后，設(shè) 其 第k 列 為h?k，定 義rk=h?k(2)+jh?k(3)，則 可 以 得 到 接 收 信 號(hào) 的 仰 角和方位角的估計(jì)分別為

式中：arcsin(·)為反正弦函數(shù)，angle(·)表示取元素相角；k=1，2，…，K。上述過(guò)程中，由于仰角和方位角從矩陣H?同一列中得到，因此所獲得的二維角度已自動(dòng)配對(duì)。此外，由于H?中只包含信號(hào)的角度信息，與陣元位置無(wú)關(guān)，因此本算法能夠應(yīng)用于任意結(jié)構(gòu)的聲矢量陣列。

2.3 算法步驟及復(fù)雜度、優(yōu)點(diǎn)分析

上述已經(jīng)給出了基于快速PARAFAC 分解的二維DOA 估計(jì)算法，其主要步驟總結(jié)如下：

（1）計(jì)算接收信號(hào)協(xié)方差矩陣，并根據(jù)式（14）對(duì)協(xié)方差矩陣進(jìn)行分塊，然后通過(guò)式（17～24）獲得矩陣A 和H 的初始估計(jì)；

（2）根據(jù)式（26，28，29），重復(fù)地更新矩陣S、A 和H 的估計(jì)值，直到收斂；

（3）根據(jù)所獲得的H 的估計(jì)值，通過(guò)式（32，33）得到接收信號(hào)的方位角和仰角估計(jì)。

對(duì)于本文所提算法，參數(shù)矩陣初估計(jì)過(guò)程計(jì)算復(fù) 雜 度 為O {16M2(J+K )+14K2M +4K3}，單次迭代分解過(guò)程復(fù)雜度為O {3K3+K2(J+M+4)+6K2(4M+JM+4J )}，因此所提算法總的復(fù)雜度為O {16M2(J+K )+14K2M+4K3+n(3K3+K2(J+M+4)+6K2(4M+JM+4J ))}。其 中n為分解迭代次數(shù)。圖2 和圖3 分別為本文所提算法與PARAFAC 算法的收斂速度對(duì)比圖以及在不同的快拍數(shù)下計(jì)算復(fù)雜度對(duì)比圖。從圖2 和圖3中可以看出，所提算法加快了迭代分解過(guò)程，計(jì)算復(fù)雜度大大降低。

圖2 迭代分解收斂速度對(duì)比Fig.2 Comparision of decompisition convergence rate

圖3 計(jì)算復(fù)雜度對(duì)比Fig.3 Comparision of computational complexity

本文所提算法有如下優(yōu)點(diǎn)：

（1）本文所提算法能夠應(yīng)用于任意結(jié)構(gòu)的聲矢量陣列，同時(shí)能得到配對(duì)的二維DOA 估計(jì)；

（2）本文所提算法借助于矩陣的初估計(jì)過(guò)程，加快了分解收斂速度，降低了算法的復(fù)雜度；

（3）本文算法角度估計(jì)性能非常接近于PARA?FAC 算法，同時(shí)優(yōu)于ESPRIT 算法和PM 算法。

3 仿真結(jié)果

根據(jù)前文描述的算法原理及應(yīng)用場(chǎng)景，本節(jié)通過(guò)MATLAB 仿真驗(yàn)證本文所提算法的有效性。在仿真中，假設(shè)有3 個(gè)不相干的遠(yuǎn)場(chǎng)信源發(fā)射的信號(hào)，被聲矢量傳感器陣列所接收，信號(hào)波達(dá)方向分別 為 (φ1，θ1)=(10°，15°)，(φ2，θ2)=(20°，25°) 和(φ3，θ3)=(30°，35°)。M、J 分別為陣元數(shù)和快拍數(shù)。仿真中，采用L 次蒙特卡洛仿真來(lái)評(píng)估算法的角度估計(jì)性能，定義角度估計(jì)的求根均方誤差（Root mean squares error，RMSE）為

式中：φk、θk分別為第k 個(gè)信源的準(zhǔn)確方位角和仰角；和分 別 為 第l 次 仿 真 中φk和θk的 估 計(jì)值；L 為仿真次數(shù)，仿真中L 取1 000。

圖4 為本文所提出的算法在信噪比（Signal noise ratio，SNR）為15 dB 下仿真100 次的估計(jì)結(jié)果分布圖。仿真中陣元數(shù)M=15、快拍數(shù)J=200。從圖4 可以看出本算法可以有效地用于聲矢量傳感器陣列的二維DOA 估計(jì)。圖5 為本文算法與PARAFAC 算 法［14］、ESPRIT 算 法［8］以 及PM 算法［9］角度估計(jì)性能對(duì)比仿真結(jié)果。仿真中M=10，J=200。從圖中可以看出，所提算法角度估計(jì)性能非常接近于PARAFAC 算法，同時(shí)優(yōu)于ESPRIT算法和PM 算法。

圖6 給出了本文算法在不同陣元數(shù)下的仿真結(jié)果，仿真中快拍數(shù)J=200。從圖6 可以看出，該算法性能隨著陣元數(shù)M 的增大而提高，其原因是陣元數(shù)的增多提高了空間分集增益。圖7 為本文算法在不同快拍數(shù)下的仿真結(jié)果，仿真中陣元數(shù)M=14。從圖7 可知，該算法角度估計(jì)誤差隨著快拍數(shù)的增大而減小，其原因?yàn)榭炫臄?shù)增大提高了時(shí)間分集增益。

圖4 本文所提算法角度估計(jì)結(jié)果Fig.4 Angle estimation of the proposed algorithm

圖5 不同算法的角度估計(jì)性能對(duì)比Fig.5 Comparison of angle estimation performance of dif?ferent algorithms

圖6 本文算法在不同陣元數(shù)下的角度估計(jì)性能Fig.6 Angle estimation performance under different M of the proposed algorithm

圖7 在不同快拍數(shù)下的角度估計(jì)性能Fig.7 Angle estimation performance under different J of the proposed algorithm

4 結(jié) 論

針對(duì)任意結(jié)構(gòu)的聲矢量傳感器陣列中的二維DOA 估計(jì)問(wèn)題，本文提出的快速PARAFAC 分解算法，能夠以較低的計(jì)算復(fù)雜度得到接收信號(hào)匹配的二維DOA 估計(jì)。仿真表明，該算法的二維角度估計(jì)性能接近于PARAFAC 算法，同時(shí)優(yōu)于ES?PRIT 算 法 和PM 算 法。