阵列信号处理笔记-波达方向DOA-子空间方法

前言

波达方向-Direction Of Arrival是研究波束形成的重要课题，引用之前的老图，DOA要估算出来的就是两个角度：俯仰角$\phi$和方位角$\theta$，而如果是全向麦克风组成的线性阵列（z轴），那么方位角就可以省略了，只研究俯仰角就可以了，所以很多算法简化假设条件，然而实际中无法省略，不过可以通过阵列的摆放，设计的估算角度范围小一些，算法也会容易一些。

经典的DOA估算方法有波束形成测向方法（base），Capon最小功率估计器，ML极大似然估计器，MUSIC 多重信号分类方法，ESPRIT旋转不变量信号参数估计方法，随着时代进步，开源代码提供了丰富晚上的功能，本文跟随pyroomacoustic实现的方法，在参考文献中找出来这些神作，学习记录，以备后用。要研究这些算法，我们必须先了解一下子空间的基础知识。

子空间和子空间数据模型

说句题外话，我已经记不起来大学里学没学过这个概念了，不过在科幻电影主导荧幕的21世纪，提到多维空间的子空间是神秘而且令人兴奋的。只是线代里的子空间定义非常滴扼要，以至于看完了概念仍然一头雾水，尽管如此还是现学现卖了抄写一下公式吧。

子空间

【from 百度百科】假设设W是域P上的线性空间V的一个非空子集合，若对于V中的加法及域P与V的纯量乘法构成域P上的一个线性空间，则称W为V的线性子空间（或向量子空间），或简称子空间。
换句话讲，V的非空子集W是子空间的充分必要条件是：
（1）子集合W的任意两个向量α与β之和α+β仍是W中的向量；
（2）域P的任一数k与子集合W的任意一个向量α的积kα仍是W中的向量。
用直观的概念来理解，就是三维坐标中，任意由两个轴组成的平面即是这三维空间的二维子空间，所以多维空间中，子空间也有很多性质，都可以用三维空间来联想一下，方便记忆和理解。

子空间数据模型

【9】将阵列信号看成矩阵$A$的参数估计问题，假设：
$$x(t)=\mathbf{A}(\mathbf{\eta })s(t)+n(t)$$
这里先重新声明一下变量的意义和所属集合：
$x(t)\in {\mathbb{C}}^M$，可以理解为有M个传感器变换的矩阵，在$t$时刻采集到的信号；
$s(t)\in {\mathbb{C}}^k$，是远场输入的k个信号，$n(t)\in {\mathbb{C}}^M$是加性噪声；
$\mathbf{A}\in {\mathbb{C}}^{M\times k}$是转换矩阵，假设在观测窗内具有时不变特性，而$\mathbf{\eta }\in {\mathbb{R}}^q$是决定这个矩阵的参数向量，如果是线阵理想的情况，就是入射角$\theta$。
另外分析此类问题，还是离不开假设：

1. $\mathbf{A}(\mathbf{\eta })$在观测周期时不变；
2. 模型在 $\mathbf{A}$和$s$上满足双线性（bilinear）；
3. 加性噪声；

基于以上所有的假设和定义，那么当$k<M$，或者说入射信号数量少于观测传感器数量的情况，那么信号$x(t)\in {\mathbb{C}}^M$能够被压缩到${\mathbb{C}}^M$一个$k$维子空间${\mathbb{C}}^k$，这个子空间被称为“信号子空间”。而噪声被假设为在所有维度上都有能量分布的，所以，此时阵列信号公式可以理解为满秩（full-rank）噪声数据模型所污染的低秩（low-rank）信号，这就是基于子空间的数据模型定义。
这个模型的几何意义，可以用三维坐标来说明，下图引自知乎上的博文：

如果$x(t)\in {\mathbb{C}}^3$被变换到子空间${\mathbb{C}}^2$子空间来分析，那么弥散在整个三维空间的噪声大部分被降维规避了，分布于${\mathbb{C}}^2$上的噪声能量势必远小于整个三维空间。

子空间模型的参数估计

本段基本翻译理解【9】，基于子空间模型的参数估计大致分三步：

1.参数建模

\quad 建模一个参数$\mathbf{\eta }$的估计向量$$\mathbf{\eta }=[{\mathbf{\eta }}_1,...,{\mathbf{\eta }}_k]$$由（列）向量组成（矩阵）$$\mathbf{A}(\mathbf{\eta })=[\mathbf{a}({\mathbf{\eta }}_1),...,\mathbf{a}({\mathbf{\eta }}_{\mathbf{k}})]$$ $\mathbf{A}(\mathbf{\eta })$中的列向量$\mathbf{a}({\mathbf{\eta }}_{\mathbf{k}})$表示为第$k$个信号源的响应，线阵的方向列向量$$a({\theta }_k)=[1,e^{-j2\pi dsin\theta k/\lambda 0},..e^{-j2\pi (M-1)dsin\theta k/\lambda 0}]$$就是在线阵DOA估计中的一个实例。而一般情况下，${\mathbf{\eta }}_{\mathbf{k}}$本身也是多维向量，有可能包含坐标信息，信号载波频率等等，从算法上只要这个向量的维度（假设$p$）小于阵元数量$M$，那么$\mathbf{a}({\mathbf{\eta }}_{\mathbf{k}})$能够在${\mathbb{C}}^M$空间中，以${\mathbf{\eta }}_{\mathbf{k}}$作为变量，绘制一个$p$维曲面，这个曲面被叫做阵列流形。数学上定义如下为 A = { a ( η k ) : η ∈ Θ } \mathscr{A}=\{\boldsymbol{a(\eta_k}): \boldsymbol{\eta} \in \boldsymbol{\Theta}\} A={a(ηk​):η∈Θ}这里$\mathbf{\Theta }$是参数向量所有可能取值的集合。很显然线阵的DOA情况$p=1$，${\mathbf{\eta }}_{\mathbf{k}}={\theta }_k$，那么$\mathcal{A}$是一个一维的曲线，这个曲线向绳子一样贯穿整个${\mathbb{C}}^M$空间，在【9】中有一个图示如下：

从几何上来看，窄带DOA估计问题可以理解为流形和信号子空间的交汇的地方（懵懂中），这就是我们期望的估计点。

一般情况下，如果$\mathcal{A}$是确定的（unambiguous），并且阵元的个数$M$多于信源的个数$k$,那么矩阵$\mathbf{A}(\mathbf{\eta })$将是一个k维满秩的矩阵。集合阵元$X$和信源$S$的信号向量可以定义：$$X=\mathbf{A}(\mathbf{\eta })S$$按照子空间的理论，这个瞬间观察值（snapshot） 是矩阵$\mathbf{A}(\mathbf{\eta })$列的线性组合，或者说是每次观测都被限定在${\mathbb{C}}^M$的$k$维子空间，这个子空间是被$\mathbf{A}(\mathbf{\eta })$的$k$列所定义的。而如果信源$S$本身是满秩的（$k$），那么$\mathbf{A}(\mathbf{\eta })$的生成子空间将等价于$\mathbf{A}$的生成子空间$span(\mathbf{A})=span(\mathbf{A}(\mathbf{\eta }))$，那么对阵元的观察向量$X$将填满这个低秩子空间，为后面的矩阵参数估计提供了数学假设和基础。在加入了零均值加性噪声之后的阵元信号$$X=\mathbf{A}(\mathbf{\eta })S+N$$

噪声就像鬼魂一样无处不在，看不见摸不着。但如鬼最好画，噪声也被各路大神描绘的惟妙惟肖，在这里将噪声假设成复稳态循环高斯随机过程，采样点之间完全不相关，以及空间协方差矩阵被定义为 σ 2 Σ = ξ { n ( t ) n ∗ ( t ) } \sigma^2\Sigma=\boldsymbol{\xi}\{n(t)n^*(t)\} σ2Σ=ξ{n(t)n∗(t)} ，这样来表示噪声的数学期望，$\Sigma$是归一化的矩阵即行列式$det(\Sigma )=1$，${\sigma }^2$表示的噪声能量。更大胆的假设是空间白噪声，那么用单位阵$I$直接替代$\Sigma$。
\quad 开完了噪声的玩笑，但引入噪声后的实际信号向量会变得满秩，就是无法用低秩序列来直接表征了。如果说不考虑噪声的情况，利用$N\ge k$个数据向量所定义的信号子空间，寻找阵列流形和子空间的焦点作为解决方案。当噪声进来以后，就是要从数据中估计出来信号子空间，进而得到$\eta$，让生成的阵列流形最佳适配这个估计。实值上，带噪信号的处理效率是子空间方法的关键问题。接下来我们看看带噪的子空间估计（SVD方法此处被提到），阵列协方差矩阵${\mathbf{R}}_{XX}$可以被写成$${\mathbf{R}}_{XX}=\mathbf{E}[X{X}^{\ast }]=\mathbf{A}(\mathbf{\eta }){\mathbf{R}}_{SS}{\mathbf{A}}^{\ast }(\mathbf{\eta })+{\sigma }^2\mathbf{I}$$这里${\mathbf{R}}_{SS}=\mathbf{E}[S{S}^{\ast }]$是信源协方差矩阵，这个矩阵是无法测量的，但从前面的推导，可以看出$\mathbf{A}(\mathbf{\eta }){\mathbf{R}}_{SS}{\mathbf{A}}^{\ast }(\mathbf{\eta })$这个矩阵有一个低秩的特点，因为信源的个数也是茫然的，所以假设这是一个秩为$d^{\prime }<M$的矩阵，通过特征值分解（eigendecomposition）得出$${\mathbf{R}}_{XX}=\sum _{k=1}^M{\lambda }_k{e}_k{e}_k^{\ast }$$
${\lambda }_1>{\lambda }_2>.....{\lambda }_M$为特征值， { e k } \{e_k\} {ek​}是相对应的特征向量。根据$${\mathbf{R}}_{XX}=\mathbf{A}(\mathbf{\eta }){\mathbf{R}}_{SS}{\mathbf{A}}^{\ast }(\mathbf{\eta })+{\sigma }^2\mathbf{I}$$公式，${\mathbf{R}}_{XX}$是有一个秩为$d^{\prime }$的矩阵再加上一个定标的单位阵，这样会让你容易的推算出最后的$M-d^{\prime }$个特征值是非常小的，即${\lambda }_{d^{\prime }+1}={\lambda }_{d^{\prime }+2}=.....{\lambda }_M={\sigma }^2$，由此，定义$${\mathbf{E}}_{\mathbf{S}}=[e_1,e_2,...e_{d^{\prime }}]$$为信号子空间，和定义$${\mathbf{E}}_{\mathbf{N}}=[e_{d^{\prime }+1},e_{d^{\prime }+2},...e_M]$$作为噪声子空间，这两个空间是正教的（orthogonal complement）。而信号子空间包含于$\mathbf{A}(\mathbf{\eta })$的生成空间。到此，基于子空间估计的建模算是推导完毕。

信号采样

在测量空间，采样获取信号子空间的估计值$\hat{\mathbf{A}}$；实际中我们是估计${\hat{\mathbf{R}}}_{XX}$，那么基于N次采样的估计公式为$${\hat{\mathbf{R}}}_{XX}=\frac{1}{N}\sum _{t=1}^{N}x(t)x^{\ast }(t)$$经过特征分解后的$({\hat{\lambda }}_k,{\hat{e}}_k),1\le k\le M$是${\hat{\mathbf{R}}}_{XX}$的特征对，其中对最大的特征值进行保存，形成了信号子空间的持续估计，而信号子空间的维度也是基于特征值分布得到的。经典的方法有似然比（ikelihood ratio），MDL(minimum description length)和AIC(akaike information criterion)。到此可以看出计算矩阵的特征值分解将是最大的负担。

参数估计

估计一组参数$\hat{\mathbf{\eta }}$，使得 $\mathbf{A}(\hat{\mathbf{\eta }})$能够从某个层面最好的匹配$\hat{\mathbf{A}}$；但这个$\hat{\mathbf{\eta }}$如何与$\hat{\mathbf{A}}$建立联系，最终得到一直想要却还没出现的DOA呢？

DOA 估计问题

了解了子空间的概念和一般的分析方法，这些方法可以用在估算信号的方向、波长、频率、幅度，甚至利用多普勒频移估算信号的速度。而用在DOA上，是如何建模的呢？这种方法叫做子空间拟合（subspace fitting）。【9】子空间拟合办法就是一个公式$$\left[\hat{\mathbf{A}},\hat{\mathbf{T}}\right]=\arg \underset{A,T}{\min }||\mathbf{M}-\mathbf{A}(\mathbf{\eta })\mathbf{T}|{|}_F^2$$相比于之前的一堆公式，多出来的$M$是一个$m\times q$矩阵，这个矩阵要从采样数据中重建。$m\times p$的矩阵是利用$\eta$参数化的，联想一下线阵的角向量。$\mathbf{T}$是一个任意$p\times q$矩阵，其中一个办法是替换一个$$\hat{\mathbf{T}}={\mathbf{A}}^{†}\mathbf{M}$$，那个奇怪的符号表示conjugate transpose，又称埃尔米特共轭。用这个替换重新获得一个等价的问题 A ^ = arg ⁡ max ⁡ A T r { P A M M ∗ } , P A = A ( A ∗ A ) − 1 A ∗ \hat \mathbf A=\arg \max\limits_{A}Tr\{P_AMM^*\}, P_A=A(A^*A)^{-1}A^* A^=argAmax​Tr{PA​MM∗},PA​=A(A∗A)−1A∗$P_A$叫做投影矩阵，投影到$\mathbf{A}$的列空间。
\quad 矩阵$M$的选择和包含$\mathbf{A}$的集合就成了解决子空间拟合的基本问题，在【9】中有一个概括的表格，比较了各种子空间拟合方法：
M & A { A ∈ A d } { A ∈ A } { A ∈ ξ } M M ∗ = R ^ X X D e t − M L B e a m f o r m i n g M L − E S P R I T M = E ^ S M D − M U S I C M U S I C T L S − E S P R I T M = E ^ S W o p t 1 / 2 W S F W e i g h t e d   M U S I C W e i g h t e d   E S P R I T \begin{array}{c|lcr} {M\&A} & \{A\in\mathscr{A}^d\} & \{A\in\mathscr{A}\} & \{A\in\boldsymbol{\xi}\}\\ \hline MM^*=\boldsymbol{\hat R}_{XX} & Det-ML & Beamforming & ML-ESPRIT \\ M=\hat E_S & MD-MUSIC & MUSIC & TLS-ESPRIT \\ M= \hat E_SW_{opt}^{1/2}& WSF & Weighted\ MUSIC & Weighted\ ESPRIT \end{array} M&AMM∗=R^XX​M=E^S​M=E^S​Wopt1/2​​{A∈Ad}Det−MLMD−MUSICWSF​{A∈A}BeamformingMUSICWeighted MUSIC​{A∈ξ}ML−ESPRITTLS−ESPRITWeighted ESPRIT​​

Deterministic maximum likelihood方法

这种方法不是基于统计模型设计的，这样数据阵$X_N$最大化似然函数被认为是这样一个等价问题 min ⁡ η , S N T r { ( X N − A ( η ) S N ) ∗ ( X N − A ( η ) S N ) } = min ⁡ η , S N ∣ ∣ ( X N − A ( η ) S N ) ∣ ∣ F 2 \min\limits_{\eta,S_N}Tr\{(X_N-\mathbf{A}(\boldsymbol{\eta}) \mathbf S_N)^*(X_N-\mathbf{A}(\boldsymbol{\eta}) \mathbf S_N)\}=\min\limits_{\eta,S_N}||(X_N-\mathbf{A}(\boldsymbol{\eta}) \mathbf S_N)||_F^2 η,SN​min​Tr{(XN​−A(η)SN​)∗(XN​−A(η)SN​)}=η,SN​min​∣∣(XN​−A(η)SN​)∣∣F2​一般认为$S_N$是确定的，只是未知，所以假设${\hat{\mathbf{S}}}_N={\mathbf{A}}^{†}{\mathbf{X}}_N$。代入上式替换${\mathbf{S}}_N$，经过推导（我寄几是推不出来）最后得出优化问题的公式 η ^ = arg ⁡ max ⁡ η T r { P A ( η ) R ^ X X } \hat \mathbf {\eta}=\arg \max\limits_{\eta}Tr\{ \boldsymbol{P_{A(\eta)}}\boldsymbol{\hat R}_{XX}\} η^​=argηmax​Tr{PA(η)​R^XX​}, ${\hat{\mathbf{R}}}_{XX}=samplecovariancematrix$。这个方法是上面表格的第一个方法，计算负担很重，实际用的已经不多。

Beam forming方法

这里是指传统的’delay-and-sum method’，核心思想就是寻找加权输出的最大功率角度，阵列会扫描所有可能的方向已获得最优解。从公式的角度是取建模归一化的输出功率$P({\eta }_j)$，利用列向量$a({\eta }_j)$以及矩阵方法求解 P ( η j ) = 1 N ∑ t = 1 N ∣ a ∗ ( η j ) x ( t ) ∣ 2 ∣ a ( η j ) ∣ 2 = a ∗ ( η j ) R ^ X X a ( η j ) a ∗ ( η j ) a ( η j ) = T r { P a ( η j ) R ^ X X } P(\eta_j)=\frac{1}{N}\sum_{t=1}^N\frac{|a^*(\eta_j)x(t)|^2}{|a(\eta_j)|^2}=\frac{a^*(\eta_j)\hat R_{XX}a(\eta_j)}{a^*(\eta_j)a(\eta_j)}=Tr\{ \boldsymbol{P_{a(\eta_j)}}\boldsymbol{\hat R}_{XX}\} P(ηj​)=N1​t=1∑N​∣a(ηj​)∣2∣a∗(ηj​)x(t)∣2​=a∗(ηj​)a(ηj​)a∗(ηj​)R^XX​a(ηj​)​=Tr{Pa(ηj​)​R^XX​}，$P_{a({\eta }_j)}=a(a^{\ast }a{)}^{-1}{a}^{\ast }$，最大功率方向就是DOA的方向。比较DML方法，如果信源只有一个，两种方法是等价的；而分离度比较好的声源，两者的估计也相当。而要分辨靠得很近的声源辐射，阵列孔径的音素不得不考虑。

MUSIC和ESPRIT

Multiple Signal Classification多信号分类方法是最经典的空间谱估计方法，也是在加性噪声环境下成功的构建数据模型的第一个高精度算法。而Estimation of Signal Parameters via Rotational Invariance Techniques 则对MUSIC的精度和限制进行了改进，有人戏称ESPRIT是MUSIC的儿子【8】，两者都是以子空间（subspace）技术为基础的方法，认为信号和噪声是在不同的子空间。pyroomacoustic没有采用ESPRIT，猜想是因为MUSIC本身的稳定性和鲁棒性更好一些。

MUSIC

“音乐”方法是1981年Ralph Otto Schmidt在斯坦福的博士论文《A Signal Subspace Approach to Multiple Emmitter Location and Spectral Estimation》提出来的方法，距今已经有40年了。他在Multiple emitter location and signal parameter estimation文章里有算法的经典推导过程。

SRP-PHAT

Steered Response Power – Phase Transform

CSSM

Coherent Signal Subspace Method

WAVES

Weighted Average of Signal Subspaces

TOPS

Test of Orthogonality of Projected Subspaces

FRIDA

Finite Rate of Innovation Direction of Arrival

术语和约定-notational conventions

ARMA

Autoregressive moving average model 自回归滑动平均模型， 是研究时间序列的重要方法，由自回归模型（简称AR模型）与移动平均模型（简称MA模型）为基础“混合”构成。

FSD

fast subspace decomposition

low-rank && full-rank

azimuth and elevation angles

方位角和仰角

wavefronts：planar or curvature

complex：in-phase and quadrature

coherent & incoherent

covariance协方差

协方差（Covariance）在概率论和统计学中用于衡量两个变量的总体误差。而方差是协方差的一种特殊情况，即当两个变量是相同的情况。协方差表示的是两个变量的总体的误差，这与只表示一个变量误差的方差不同。 如果两个变量的变化趋势一致，也就是说如果其中一个大于自身的期望值，另外一个也大于自身的期望值，那么两个变量之间的协方差就是正值。 如果两个变量的变化趋势相反，即其中一个大于自身的期望值，另外一个却小于自身的期望值，那么两个变量之间的协方差就是负值。（摘自百度百科）

协方差矩阵

统计学与概率论中，协方差矩阵的每个元素是各个向量元素之间的协方差，是从标量随机变量到高维度随机向量的自然推广。

特征向量、特征值

Tr. 矩阵的迹

Maximum likelihood estimation

似然和概率是鸡和蛋的问题，在给定条件$\theta$下计算事件出现的概率可以表达为$$P(x|\theta )$$，而根据事件发生的概率（结果）反推出条件就是似然的过程$$\mathcal{L}(\theta |x)$$，所以最大似然估计是利用已知的样本的结果，在使用某个模型的基础上，反推最有可能导致这样结果的模型参数值。很显然求解DOA的过程也可以利用最大似然估计的办法。

Singular Value Decomposition（奇异值分解）

参考文献

【1】.R. Schmidt, Multiple emitter location and signal parameter estimation, IEEE Trans. Antennas Propag., Vol. 34, Num. 3, pp 276–280, 1986
【2】.J. H. DiBiase, A high-accuracy, low-latency technique for talker localization in reverberant environments using microphone arrays, PHD Thesis, Brown University, 2000
【3】.H. Wang, M. Kaveh, Coherent signal-subspace processing for the detection and estimation of angles of arrival of multiple wide-band sources, IEEE Trans. Acoust., Speech, Signal Process., Vol. 33, Num. 4, pp 823–831, 1985
【4】.E. D. di Claudio, R. Parisi, WAVES: Weighted average of signal subspaces for robust wideband direction finding, IEEE Trans. Signal Process., Vol. 49, Num. 10, 2179–2191, 2001
【5】.Y. Yeo-Sun, L. M. Kaplan, J. H. McClellan, TOPS: New DOA estimator for wideband signals, IEEE Trans. Signal Process., Vol. 54, Num 6., pp 1977–1989, 2006
【6】.H. Pan, R. Scheibler, E. Bezzam, I. Dokmanic, and M. Vetterli, FRIDA: FRI-based DOA estimation for arbitrary array layouts, Proc. ICASSP, pp 3186-3190, 2017
【7】.ESPRIT-Estimation of Signal Parameters Via Rotational Invariance Techniques
【8】.ESPRIT Algorithm- The Son of MUSIC
【9】Subspace Methods for Directions-of-Arrival, A. Paulraj, B. Ottersten, R. Roy, A. Swindlehurst, G. Xu and T. Kailath
【10】Performance Analysis of MUSIC and ESPRIT DOA Estimation Algorithms for Adaptive Array Smart Antenna in MobileCommunication
【11】Direction-of-Arrival Estimation Methods: APerformance-Complexity Tradeoff Perspective
【】