CN104375116A - 一种基于无线传感器阵列的到达方向检测方法 - Google Patents

Abstract

本发明公开了一种基于无线传感器阵列的到达方向检测方法，属于雷达技术领域。本发明通过所求导向矢量构造函数，获取上述函数的多个零点，通过最大零点反解方位角，导向矢量的构造考虑了二维角度对应的所有阵元接收数据的范德蒙结构，克服了ESPRIT方法只通过单个特征值确定角度的弊端，提高了估计的准确性。

Description

一种基于无线传感器阵列的到达方向检测方法

技术领域

本发明涉及雷达技术领域，特别涉及一种基于无线传感器阵列的到达方向检测方法。

背景技术

在雷达领域中，确定电磁波信源的到达方向一直是研究的重要课题。

在现有的技术中，常用的有最大似然和MUSIC(Multiple SignalClassification，多信号分类)和ESPRIT(Estimating Signal Parameters viaRotational Invariance Techniques，借助旋转不变技术估算信号参数)方法，其中ESPRIT通过计算闭式解的方法，就可以得到信源的方位角和俯仰角两个重要参数，从而完成对DOA的估计，不需像最大似然和MUSIC方法那样要经过对谱峰进行搜索，可以显著降低相关数据的计算量和存储量。

在实现本发明的过程中，发明人发现现有技术至少存在以下问题：

作为其中较为优秀的方法，ESPRIT方法在运算时需要分维计算和参数配对，例如在进行一维分量重叠时会出现错误配对的情况，此时就会出现计算错误，导致无法对电磁波信源的到达方向进行准确估计。

发明内容

为了解决现有技术的问题，本发明提供了一种基于无线传感器阵列的到达方向检测方法，所述方法包括：

步骤一，部署信号接收阵列，所述信号接收阵列中包含有线性排列的M+1个传感器，其中相邻的两个传感器之间的间距d小于波长λ的二分之一；

步骤二，通过所述信号接收阵列接收到信源的第一信号值x₁(t)＝A₁s(t)+n₁(t)、第二信号值x₂(t)＝A₂s(t)+n₁(t)＝A₁Φs(t)+n₁(t)，将所述第一信号值x₁(t)、第二信号值x₂(t)分别进行相关运算，得到相关矩阵

$R_{11}=E\left(x_1{x}_1^H\right)=A_1{R}_s{A}_1^H+{\mathrm{\σ}}^2{I}_M,$

$R_{21}=E\left(x_2{x}_1^H\right)=A_1\mathrm{\Φ}{R}_s{A}_1^H,$

其中，R_s为所述信源的自相关矩阵，σ²为噪声平均功率，对所述噪声平均功率进行估计，得到去噪后的矩阵

C₁₁＝R₁₁-σ²I_M，

C₂₁＝R₂₁。

步骤三，获取双重旋转去噪矩阵 $C={\left[\begin{array}{cc}{C}_{11}^T& C_{21}^T\end{array}\right]}^T,$ 将所述双重旋转去噪矩阵进行奇异值分解，得到分解后的矩阵C_svd＝U₁Λ₁V₁ ^H，将所述分解后的矩阵中U₁进行分块处理，得到

$C_{\mathrm{svd}}=\left[\begin{array}{cc}\begin{array}{c}{U}_4\\ U_5\end{array}& U_3\end{array}\right]{\mathrm{\Λ}}_1{V_1}^H,$

在所述双重旋转去噪矩阵中代入x₁(t)＝A₁s(t)+n₁(t)、x₂(t)＝A₁Φs(t)+n₁(t)，有 $C={\left[\begin{array}{cc}{C}_{11}^T& C_{21}^T\end{array}\right]}^T={\left[\begin{array}{cc}{A}_1^T& {\left(A_1\mathrm{\Λ}\right)}^T\end{array}\right]}^T,$ 结合

$C_{\mathrm{svd}}=\left[\begin{array}{cc}\begin{array}{c}{U}_4\\ U_5\end{array}& U_3\end{array}\right]{\mathrm{\Λ}}_1{V_1}^H$

即存在 ${\left[\begin{array}{cc}{A}_1^T& {\left(A_1\mathrm{\Λ}\right)}^T\end{array}\right]}^T$ 和 ${\left[\begin{array}{cc}{\left(U_4\right)}^T& {\left(U_5\right)}^T\end{array}\right]}^T$ 张成相同的列空间，有A₁＝U₄F和A₁Λ＝U₅F，进一步有(U₄)^-1U₅＝FΛF^-1，对(U₄)^-1U₅＝FΛF^-1进行特征值分解，获取导向矢量

步骤四，选取所述导向矢量中第m行p列元素为a_mp，构造函数ρ(μ)如下：

$\mathrm{\ρ}\left(\mathrm{\μ}\right)=\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{a}_{m-1,p}{a}_{\mathrm{mp}}^{*}{e}^{\mathrm{\β}(m-1)\mathrm{\μ}},$

其中，μ＝sinθ_p，对所述函数ρ(μ)的模值进行平方，整理后有

${\left|\mathrm{\ρ}\left(\mathrm{\μ}\right)\right|}^2=\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{a}_{m-1,p}{a}_{\mathrm{mp}}^{*}{e}^{\mathrm{\β}(m-1)\mathrm{\μ}}\underset{n=2}{\overset{M}{\mathrm{\Σ}}}{a}_{n-1,p}^{*}{a}_{\mathrm{np}}{e}^{-\mathrm{\β}(n-1)\mathrm{\μ}}=\underset{m=2}{\overset{M}{\mathrm{\Σ}}}\underset{n=2}{\overset{M}{\mathrm{\Σ}}}{a}_{m-1,p}{a}_{\mathrm{mp}}^{*}{a}_{n-1,p}^{*}{a}_{\mathrm{np}}{e}^{\mathrm{\β}(m-n)\mathrm{\μ}},$

在所述|ρ(μ)|²中对μ求导，得到

$\frac{\∂{\left|\mathrm{\ρ}\left(\mathrm{\μ}\right)\right|}^2}{\∂\mathrm{\μ}}=\mathrm{\β}\underset{m=2}{\overset{M}{\mathrm{\Σ}}}\underset{n=2}{\overset{M}{\mathrm{\Σ}}}(m-n)a_{m-1,p}{a}_{\mathrm{mp}}^{*}{a}_{n-1,p}^{*}{a}_{\mathrm{np}}{\left(e^{\mathrm{\β\μ}}\right)}^{(m-n)},---\left(1\right)$

令e^βμ＝δ，并令公式(1)为0，有

$\underset{m=2}{\overset{M}{\mathrm{\Σ}}}\underset{n=2}{\overset{M}{\mathrm{\Σ}}}(m-n)a_{m-1,p}{a}_{\mathrm{mp}}^{*}{a}_{n-1,p}^{*}{a}_{\mathrm{np}}{\left(\mathrm{\δ}\right)}^{(m-n)}=0,---\left(2\right)$

获取公式(2)的所有零点中最大的零点δ_p，从所述δ_p反解出天线的方位角的估计值

${\widehat{\mathrm{\θ}}}_p=\frac{180}{\mathrm{\π}}\arcsin \left[\frac{\mathrm{\λ}}{2\mathrm{\πd}}\mathrm{angle}\left({\mathrm{\δ}}_p\right)\right],$

其中，angle(δ_p)为取复角主值。

本发明提供的技术方案带来的有益效果是：

通过所求导向矢量构造函数，获取上述函数的多个零点，通过最大零点反解方位角，导向矢量的构造考虑了二维角度对应的所有阵元接收数据的范德蒙结构，克服了ESPRIT方法只通过单个特征值确定角度的弊端，提高了估计的准确性。

附图说明

为了更清楚地说明本发明的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明提供的一种基于无线传感器阵列的到达方向检测方法的流程示意图；

图2是本发明提供的接收信号的模型示意图；

图3是本发明提供的分辨成功概率比较示意图；

图4是本发明提供的角度估计的均方根误差随信噪比变化的对比示意图；

图5是本发明提供的角度估计的均方根误差随快拍数变化的对比示意图；

图6是本发明提供的角度估计的均方根误差随振元数变化的对比示意图。

具体实施方式

为使本发明的结构和优点更加清楚，下面将结合附图对本发明的结构作进一步地描述。

实施例一

本发明提供了一种基于无线传感器阵列的到达方向检测方法，如图1所示，该方法包括：

步骤一，部署信号接收阵列，信号接收阵列中包含有线性排列的M+1个传感器，其中相邻的两个传感器之间的间距d小于波长λ的二分之一。

在实施中，天线由M+1个均匀线性传感器组成,它由两个具有相同结构的子阵组成，两个子阵列均含M个天线阵元，具体结构如图2所示，第1个子阵包含阵元1到M，第2个子阵包含阵元2到M+1，M的取值为从1开始的自然数。

本步骤使用线性传感器对信源发射的信号进行接收，用于对信源方位角这个一维参数进行估计。

步骤二，通过信号接收阵列接收到信源的第一信号值x₁(t)＝A₁s(t)+n₁(t)、第二信号值x₂(t)＝A₂s(t)+n₁(t)＝A₁Φs(t)+n₁(t)，将第一信号值x₁(t)、第二信号值x₂(t)分别进行相关运算，得到相关矩阵

$R_{11}=E\left(x_1{x}_1^H\right)=A_1{R}_s{A}_1^H+{\mathrm{\σ}}^2{I}_M,$

$R_{21}=E\left(x_2{x}_1^H\right)=A_1\mathrm{\Φ}{R}_s{A}_1^H,$

其中，R_s为信源的自相关矩阵，σ²为噪声平均功率，对噪声平均功率进行估计，得到去噪后的矩阵

C₁₁＝R₁₁-σ²I_M，

C₂₁＝R₂₁。

在实施中，信源中设置有多个信号发射目标，设其中第p个信号发射目标的方位角为θ_p。在第t(t＝1,2,…,T)次采样下信源所有p个信号发射目标在第m个阵元上的总接收信号为

$x_m\left(t\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{s}_p\left(t\right)e^{j\frac{2\mathrm{\π}}{\mathrm{\λ}}(m-1)d\sin \left({\mathrm{\θ}}_p\right)}+n_m\left(t\right),$

其中s_p为第p个信源的复包络，n_m为第m个阵元上的噪声。设阵元1到M为第1个子阵，阵元2到M+1为第2个子阵，则第1个子阵接收的数据向量为

x₁(t)＝[x₁(t),x₂(t),…,x_M(t)]^T，

＝A₁s(t)+n₁(t)

其中A₁＝[a₁,a₂,…,a_P]，a_p为列向量，其第i个元素为exp(j2π(i-1)θ_p/λ)，s(t)＝[s₁(t),s₂(t),…,s_P(t)]^T，n₁(t)＝[n₁(t),n₂(t),…,n_M(t)]^T，第2个子阵接收的数据向量为

x₂(t)＝A₂s(t)+n₁(t)＝A₁Φs(t)+n₁(t)，

其中Φ为对角矩阵， $\mathrm{\Φ}=\mathrm{diag}[e^{{\mathrm{j\ω}}_1},e^{{\mathrm{j\ω}}_2}...,e^{{\mathrm{j\ω}}_P}],$ ${\mathrm{\ω}}_p=\frac{2\mathrm{\π}d\sin {\mathrm{\θ}}_p}{\mathrm{\λ}},$ 噪声为n₂(t)＝[n_M+1(t),n_M+2(t),…,n_2M(t)]^T。将第一信号值x₁(t)、第二信号值x₂(t)分别进行相关运算，得到相关矩阵

$R_{11}=E\left(x_1{x}_1^H\right)=A_1{R}_s{A}_1^H+{\mathrm{\σ}}^2{I}_M,$

$R_{21}=E\left(x_2{x}_1^H\right)=A_1\mathrm{\Φ}{R}_s{A}_1^H,$

其中，R_s为信源的自相关矩阵，σ²为噪声平均功率，对R₁₁进行特征值分解，根据分解后M个特征值中P个较小的特征值的平均值对噪声平均功率进行估计，从R₁₁、R₂₁中减去得到的噪声估计，从而得到去噪后的矩阵为

C₁₁＝R₁₁-σ²I_M，

C₂₁＝R₂₁。

本步骤中，对信号接收阵列接收到的信号值进行相关运算、去噪等一系列处理，最终得到去噪后的信号矩阵，上述处理去除了其中的噪声等无关变量的影响，使得最终对信源方位的估计更加准确。

步骤三，获取双重旋转去噪矩阵 $C={\left[\begin{array}{cc}{C}_{11}^T& C_{21}^T\end{array}\right]}^T,$ 将双重旋转去噪矩阵进行奇异值分解，得到分解后的矩阵C_svd＝U₁Λ₁V₁ ^H，将分解后的矩阵中U₁进行分块处理，得到

$C_{\mathrm{svd}}=\left[\begin{array}{cc}\begin{array}{c}{U}_4\\ U_5\end{array}& U_3\end{array}\right]{\mathrm{\Λ}}_1{V_1}^H,$

在双重旋转去噪矩阵中代入x₁(t)＝A₁s(t)+n₁(t)、x₂(t)＝A₁Φs(t)+n₁(t)，有 $C={\left[\begin{array}{cc}{C}_{11}^T& C_{21}^T\end{array}\right]}^T={\left[\begin{array}{cc}{A}_1^T& {\left(A_1\mathrm{\Λ}\right)}^T\end{array}\right]}^T,$ 结合

$C_{\mathrm{svd}}=\left[\begin{array}{cc}\begin{array}{c}{U}_4\\ U_5\end{array}& U_3\end{array}\right]{\mathrm{\Λ}}_1{V_1}^H$

即存在 ${\left[\begin{array}{cc}{A}_1^T& {\left(A_1\mathrm{\Λ}\right)}^T\end{array}\right]}^T$ 和 ${\left[\begin{array}{cc}{\left(U_4\right)}^T& {\left(U_5\right)}^T\end{array}\right]}^T$ 张成相同的列空间，有A₁＝U₄F和A₁Λ＝U₅F，进一步有(U₄)^-1U₅＝FΛF^-1，对(U₄)^-1U₅＝FΛF^-1进行特征值分解，获取导向矢量

在实施中，将之前获取到的去噪矩阵C₁₁、C₂₁分别进行转置，得到对应的转置矩阵将两个转置矩阵进行矩阵乘的运算后，将得到的矩阵再次进行转置，得到转置后的矩阵C，接着将矩阵C进行奇异值分解，得到分解后的矩阵C_svd＝U₁Λ₁V₁ ^H，进一步的，将其中的部分矩阵U₁进行分块，具体的：

(1)在列的层级上，分为具有相同列数量的两部分，即 $\left[\begin{array}{c}{U}_4\\ U_5\end{array}\right]$ 和U₃两部分。

(2)在行的层级上，将 $\left[\begin{array}{c}{U}_4\\ U_5\end{array}\right]$ 分为具有相同行数量的两部分，即U₄和U₅。在分块完毕后，在矩阵C中代入x₁(t)＝A₁s(t)+n₁(t)、x₂(t)＝A₁Φs(t)+n₁(t)，即有 $C={\left[\begin{array}{cc}{C}_{11}^T& C_{21}^T\end{array}\right]}^T={\left[\begin{array}{cc}{A}_1^T& {\left(A_1\mathrm{\Λ}\right)}^T\end{array}\right]}^T,$ 与按奇异值分解后的矩阵

$C_{\mathrm{svd}}=\left[\begin{array}{cc}\begin{array}{c}{U}_4\\ U_5\end{array}& U_3\end{array}\right]{\mathrm{\Λ}}_1{V_1}^H$

对比即可发现，由于矩阵C_svd与C的关系，可以得知 ${\left[\begin{array}{cc}{A}_1^T& {\left(A_1\mathrm{\Λ}\right)}^T\end{array}\right]}^T$ 和 ${\left[\begin{array}{cc}{\left(U_4\right)}^T& {\left(U_5\right)}^T\end{array}\right]}^T$ 中所有的列进行线性组合而形成的空间相同，即 ${\left[\begin{array}{cc}{A}_1^T& {\left(A_1\mathrm{\Λ}\right)}^T\end{array}\right]}^T$ 和 ${\left[\begin{array}{cc}{\left(U_4\right)}^T& {\left(U_5\right)}^T\end{array}\right]}^T$ 可以张成相同的列空间，因此有A₁＝U₄F和A₁Λ＝U₅F，更进一步有(U₄)^-1U₅＝FΛF^-1，对(U₄)^-1U₅＝FΛF^-1进行特征值分解，获取导向矢量

在本步骤中，通过对进行两次转置后得到的矩阵C进行特征值分解，进而对分解后矩阵进行分块处理，最终得到由分块后的局部矩阵构成的导向矢量从而进一步的提高了对信源方位进行估计的准确性。

步骤四，选取导向矢量中第m行p列元素为a_mp，构造函数ρ(μ)如下：

$\mathrm{\ρ}\left(\mathrm{\μ}\right)=\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{a}_{m-1,p}{a}_{\mathrm{mp}}^{*}{e}^{\mathrm{\β}(m-1)\mathrm{\μ}},$

其中，μ＝sinθ_p，对函数ρ(μ)的模值进行平方，整理后有

${\left|\mathrm{\ρ}\left(\mathrm{\μ}\right)\right|}^2=\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{a}_{m-1,p}{a}_{\mathrm{mp}}^{*}{e}^{\mathrm{\β}(m-1)\mathrm{\μ}}\underset{n=2}{\overset{M}{\mathrm{\Σ}}}{a}_{n-1,p}^{*}{a}_{\mathrm{np}}{e}^{-\mathrm{\β}(n-1)\mathrm{\μ}}=\underset{m=2}{\overset{M}{\mathrm{\Σ}}}\underset{n=2}{\overset{M}{\mathrm{\Σ}}}{a}_{m-1,p}{a}_{\mathrm{mp}}^{*}{a}_{n-1,p}^{*}{a}_{\mathrm{np}}{e}^{\mathrm{\β}(m-n)\mathrm{\μ}},$ 在|ρ(μ)|²中对μ求导，得到

$\frac{\∂{\left|\mathrm{\ρ}\left(\mathrm{\μ}\right)\right|}^2}{\∂\mathrm{\μ}}=\mathrm{\β}\underset{m=2}{\overset{M}{\mathrm{\Σ}}}\underset{n=2}{\overset{M}{\mathrm{\Σ}}}(m-n)a_{m-1,p}{a}_{\mathrm{mp}}^{*}{a}_{n-1,p}^{*}{a}_{\mathrm{np}}{\left(e^{\mathrm{\β\μ}}\right)}^{(m-n)},---\left(1\right)$

令e^βμ＝δ，并令公式(1)为0，有

$\underset{m=2}{\overset{M}{\mathrm{\Σ}}}\underset{n=2}{\overset{M}{\mathrm{\Σ}}}(m-n)a_{m-1,p}{a}_{\mathrm{mp}}^{*}{a}_{n-1,p}^{*}{a}_{\mathrm{np}}{\left(\mathrm{\δ}\right)}^{(m-n)}=0,---\left(2\right)$

获取公式(2)的所有零点中最大的零点δ_p，从δ_p反解出天线的方位角的估计值

${\widehat{\mathrm{\θ}}}_p=\frac{180}{\mathrm{\π}}\arcsin \left[\frac{\mathrm{\λ}}{2\mathrm{\πd}}\mathrm{angle}\left({\mathrm{\δ}}_p\right)\right],$

其中，angle(δ_p)为取复角主值。

在实施中，构造一个与导向矢量估计矩阵中的元素相关的函数ρ(μ)，通过对ρ(μ)一系列的变换、求导后，得到最大零点，从该零点反解出的角度方

程如 ${\widehat{\mathrm{\θ}}}_p=\frac{180}{\mathrm{\π}}\arcsin \left[\frac{\mathrm{\λ}}{2\mathrm{\πd}}\mathrm{angle}\left({\mathrm{\δ}}_p\right)\right]$

所示，其中，angle(δ_p)为复角主值，是已知的，因此，通过上述运算，可以得到方位角的估计值。

下面为了验证本方法的准确性，建立步骤一的信号接收阵列，具体的为M＝13个间距半波长的等距线阵，接收远场p＝3个窄带信号θ＝(0°,5°,10°)，各阵元噪声为零均值白复高斯噪声，以下实验均为200次独立实验平均结果。

实验1：信噪比为-11dB到2dB，本发明算法和ESPRIT算法的分辨成功概率如图3所示。分辨成功定义为估计值与真实值的误差不超过参考角度处3dB波束宽度的一半，算法的分辨门限一般认为是分辨成功概率大于90％所对应的信噪比，可见本发明算法的分辨门限为-8dB，优于ESPRIT算法0dB的分辨门限。

实验2：信噪比为-5dB到20dB，设n＝1,2,…,N为实验次数，试验中角度估计的均方根误差RMSE定义为

$10\mathrm{lo}{g}_{10}\left(\sqrt{\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\left[\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{({\widehat{\mathrm{\θ}}}_p\left(n\right)-{\mathrm{\θ}}_p)}^2\right]}\right)$

图4所示为快拍数固定在300次，不同信噪比下ESPRIT算法和本发明算法及理论最优值的仿真曲线。由图可见，两种算法在低信噪比下性能差别较大，且两种算法离最优性能随信噪比的降低而扩大，这是由于两种算法都是建立矩阵特征分解基础上的特征结构类算法，该类算法在低信噪比下，信号空间和噪声空间的分辨误差明显增大，尤其在信号和噪声等功率的0dB以下，性能显著变差，随着信噪比的增加，在0dB到5dB之间，两种算法性能迅速接近并趋近于最优值，在信噪比大于5dB后，算法性能差距趋于稳定。

实验3：信噪比固定为5dB，快拍数为100到1000，图5所示为不同快拍下ESPRIT算法和本发明算法及最优值的仿真曲线。由图可见随着快拍数的增加，两种算法对到达方向估计的均方根误差逐渐变小。

实验4：信噪比固定为5dB，快拍数为300，图6所示为传感器阵元数变化下ESPRIT算法和本发明算法及最优值的仿真曲线。由图可见随着振元数的增加，两种方法角度估计更加精确，均方根误差不断减小。

本实施例中提出的，通过所求导向矢量构造函数，获取上述函数的多个零点，通过最大零点反解方位角，导向矢量的构造考虑了二维角度对应的所有阵元接收数据的范德蒙结构，克服了ESPRIT方法只通过单个特征值确定角度的弊端，提高了估计的准确性。

以上所述仅为本发明的实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于无线传感器阵列的到达方向检测方法，其特征在于，所述方法包括：

步骤一，部署信号接收阵列，所述信号接收阵列中包含有线性排列的M+1个传感器，其中相邻的两个传感器之间的间距d小于波长λ的二分之一；

步骤二，通过所述信号接收阵列接收到信源的第一信号值x₁(t)＝A₁s(t)+n₁(t)、第二信号值x₂(t)＝A₂s(t)+n₁(t)＝A₁Φs(t)+n₁(t)，将所述第一信号值x₁(t)、第二信号值x₂(t)分别进行相关运算，得到相关矩阵

$R_{11}=E\left(x_1{x}_1^H\right)=A_1{R}_s{A}_1^H+{\mathrm{\σ}}^2{I}_M,$

$R_{21}=E\left(x_2{x}_1^H\right)=A_1\mathrm{\Φ}{R}_s{A}_1^H,$

其中，R_s为所述信源的自相关矩阵，σ²为噪声平均功率，对所述噪声平均功率进行估计，得到去噪后的矩阵

C₁₁＝R₁₁-σ²I_M，

C₂₁＝R₂₁；

步骤三，获取双重旋转去噪矩阵 $C={\left[\begin{array}{cc}{C}_{11}^T& C_{21}^T\end{array}\right]}^T,$ 将所述双重旋转去噪矩阵进行奇异值分解，得到分解后的矩阵C_svd＝U₁Λ₁V₁ ^H，将所述分解后的矩阵中U₁进行分块处理，得到

$C_{\mathrm{svd}}=\left[\begin{array}{c}{U}_4\\ U_5\end{array}\begin{array}{c}{U}_3\end{array}\right]{\mathrm{\Λ}}_1{V_1}^H,$

在所述双重旋转去噪矩阵中代入x₁(t)＝A₁s(t)+n₁(t)、x₂(t)＝A₁Φs(t)+n₁(t)，有 $C={\left[\begin{array}{cc}{C}_{11}^T& C_{21}^T\end{array}\right]}^T={\left[\begin{array}{cc}{A}_1^T& {\left(A_1\mathrm{\Λ}\right)}^T\end{array}\right]}^T,$ 结合

$C_{\mathrm{svd}}=\left[\begin{array}{c}{U}_4\\ U_5\end{array}\begin{array}{c}{U}_3\end{array}\right]{\mathrm{\Λ}}_1{V_1}^H$

即存在 ${\left[\begin{array}{cc}{A}_1^T& {\left(A_1\mathrm{\Λ}\right)}^T\end{array}\right]}^T$ 和[(U₄)^T (U₅)^T]^T张成相同的列空间，有A₁＝U₄F和A₁Λ＝U₅F，进一步有(U₄)^-1U₅＝FΛF^-1，对(U₄)^-1U₅＝FΛF^-1进行特征值分解，获取导向矢量

步骤四，选取所述导向矢量中第m行p列元素为a_mp，构造函数ρ(μ)如下：

$\mathrm{\ρ}\left(\mathrm{\μ}\right)=\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{a}_{m-1,p}{a}_{\mathrm{mp}}^{*}{e}^{\mathrm{\β}(m-1)\mathrm{\μ}},$

其中，μ＝sinθ_p，对所述函数ρ(μ)的模值进行平方，整理后有 ${\left|\mathrm{\ρ}\left(\mathrm{\μ}\right)\right|}^2=\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{a}_{m-1,p}{a}_{\mathrm{mp}}^{*}{e}^{\mathrm{\β}(m-1)\mathrm{\μ}}\underset{n=2}{\overset{M}{\mathrm{\Σ}}}{a}_{n-1,p}^{*}{a}_{\mathrm{np}}{e}^{-\mathrm{\β}(n-1)\mathrm{\μ}}=\underset{m=2}{\overset{M}{\mathrm{\Σ}}}\underset{n=2}{\overset{M}{\mathrm{\Σ}}}{a}_{m-1,p}{a}_{\mathrm{mp}}^{*}{a}_{n-1,p}^{*}{a}_{\mathrm{np}}{e}^{\mathrm{\β}(m-n)\mathrm{\μ}},$ 在所述|ρ(μ)²中对μ求导，得到

$\frac{\∂{\left|\mathrm{\ρ}\left(\mathrm{\μ}\right)\right|}^2}{\∂\mathrm{\μ}}=\mathrm{\β}\underset{m=2}{\overset{M}{\mathrm{\Σ}}}\underset{n=2}{\overset{M}{\mathrm{\Σ}}}(m-n)a_{m-1,p}{a}_{\mathrm{mp}}^{*}{a}_{n-1,p}^{*}{a}_{\mathrm{np}}{\left(e^{\mathrm{\β\μ}}\right)}^{(m-n)},---\left(1\right)$

令e^βμ＝δ，并令公式(1)为0，有

$\underset{m=2}{\overset{M}{\mathrm{\Σ}}}\underset{n=2}{\overset{M}{\mathrm{\Σ}}}(m-n)a_{m-1,p}{a}_{\mathrm{mp}}^{*}{a}_{n-1,p}^{*}{a}_{\mathrm{np}}{\left(\mathrm{\δ}\right)}^{(m-n)}=0,---\left(2\right)$

获取公式(2)的所有零点中最大的零点δ_p，从所述δ_p反解出天线的方位角的估计值

${\widehat{\mathrm{\θ}}}_p=\frac{180}{\mathrm{\π}}\arcsin \left[\frac{\mathrm{\λ}}{2\mathrm{\πd}}\mathrm{angle}\left({\mathrm{\δ}}_p\right)\right],$

其中，angle(δ_p)为取复角主值。