CN104820216B - 基于阵列响应旋转不变性的多径信号波达角估计方法 - Google Patents

Abstract

本发明公开了一种基于阵列响应旋转不变性的多径信号波达角估计方法，用于解决现有多径信号波达角估计方法实用性差的技术问题。技术方案是首先在两个不同子阵的接收信号中分别计算观测响应矩阵，即对每一个子阵的天线振子的输出信号进行解调和低通滤波后，得到基带接收信号，再将基带接收信号与伪随机序列做滑动相关，得到两个子阵的观测冲激响应，通过对两个子阵的观测冲激响应进行离散化，抽取出多径，得到两个子阵的观测冲激响应矢量，分别将两个子阵的所有天线振子的观测冲激响应矢量整合得到两个子阵的观测响应矩阵；最后利用两个子阵接收信号协方差矩阵的旋转不变性估计多径信号的波达角。本发明方法简单，运算结果精度高，因此实用性强。

Description

基于阵列响应旋转不变性的多径信号波达角估计方法

技术领域

本发明涉及一种多径信号波达角估计方法，特别涉及一种基于阵列响应旋转不变性的多径信号波达角估计方法。

背景技术

文献“ESPRIT算法估计性能分析，《信息技术》，2011年第3期，100页至102页”中提到的ESPRIT(Estimation of Signal Parameters Via Rotational InvarianceTechniques)方法是一种典型的信号波达角估计方法。该方法将阵列分为两个完全相同的子阵，利用两个子阵接收信号协方差矩阵的旋转不变性进行波达角估计。但是该方法存在如下问题，该方法要求多径数量小于天线阵列的振子个数。而实际信道中的多径数量很大，考虑硬件成本，难以实现更大规模的接收天线阵列，因此ESPRIT方法难以推广到实际应用。

发明内容

为了克服现有多径信号波达角估计方法实用性差的不足，本发明提供一种基于阵列响应旋转不变性的多径信号波达角估计方法。该方法首先在两个不同子阵的接收信号中分别计算观测响应矩阵，即对每一个子阵的天线振子的输出信号进行解调和低通滤波后，得到基带接收信号，再将基带接收信号与伪随机序列做滑动相关，得到两个子阵的观测冲激响应，通过对两个子阵的观测冲激响应进行离散化，抽取出多径，得到两个子阵的观测冲激响应矢量，分别将两个子阵的所有天线振子的观测冲激响应矢量整合得到两个子阵的观测响应矩阵；最后利用两个子阵接收信号协方差矩阵的旋转不变性估计多径信号的波达角。本发明方法简单、运算复杂度低、运算结果精度高，适用于对大量多径信号波达角进行估计，突破了现有ESPRIT方法对于天线振元的限制，实用性强。

本发明解决其技术问题所采用的技术方案是：一种基于阵列响应旋转不变性的多径信号波达角估计方法，其特点是包括以下步骤：

(A)定义信号。

选择长度为X的伪随机序列作为基带探测信号a(t)，其表达式为

其中，X表示伪随机序列的长度，t表示时间，表示宽度为T_b的矩形脉冲信号。K个PN序列组成一个探测帧u(t)，其表达式为

其中，T_p＝XT_b。探测帧u(t)是基本探测信号，该探测帧调制后经天线发射出去。

对于包含了L条多径的传播环境，其空间信道冲激响应模型h(t)表示如下

其中，L为传播环境中的多径条数，表示第l条多径的信道复响应，是一个复常数，τ_l是第l条多径的时延值，δ(t)表示冲激函数。

接收天线阵列是一维线形天线阵列，包含M+1个天线振子，其中，M＞1。这M+1个天线振子等距离排列，振子间距表示为d，并且天线振子的方向图都相同。从几何上将天线阵列分为完全相同的两个子阵：子阵Z_α和子阵Z_β。其中，子阵Z_α由原天线阵列中的第1～M号天线振子组成，子阵Z_β由原天线阵列中第2～M+1号天线振子组成。两个子阵都包括有M个天线振子，因此把两个子阵中的天线振子都重新编号为1～M。子阵Z_α上的第m个天线振子和子阵Z_β上的第m个天线振子的距离表示为Δ。

对于子阵中的第m个天线振子，则其在的来波方向θ上的复响应s(θ)为

其中，e表示自然底数，j表示虚数，表示射频信号的波长，θ表示来波方向。

对于子阵Z_α上的第m个天线振子，按照式(5)生成射频接收信号

其中，N′_m(t)是服从高斯分布的白噪声，θ_l表示第l条径的入射角度，和τ_l为多径的复响应和时延，u′(t)为调制后的探测帧。对子阵Z_α上的所有天线振子，都生成相应的接收信号。

对于子阵Z_β上的第m个天线振子，按照式(6)生成射频接收信号

其中，N′_m(t)是服从高斯分布的白噪声，θ_l表示第l条径的入射角度，和τ_l为多径的复响应和时延，u′(t)为调制后的探测帧。对于子阵Z_β上的所有天线振子，都生成相应的接收信号。

子阵Z_α上的第m个天线振子输出的接收信号经解调后表示为

其中，N_m(t)是第m个天线振子所接收到的噪声信号，θ_l是第l条径入射信号的波达角，s_m(θ_l)是第m个天线振子对波达角为θ_l的入射信号的复响应。

子阵Z_β上的第m个天线振子输出的接收信号经解调后表示为

其中，e是自然底数，j表示虚数。由于子阵Z_β上的第m个天线振子相对于子阵Z_α上的第m个天线振子存在距离Δ，这种距离导致了来波方向上的波程差。进而对于入射角为θ_l的多径信号，子阵Z_β的接收信号相对于子阵Z_α的接收信号出现了相移ξ(θ_l)。ξ(θ_l)的表达式为

其中，表示射频信号的波长。

(B)计算观测冲激响应矩阵。

对于子阵Z_α和子阵Z_β上的第m个天线振子，计算观测冲激响应矩阵分为如下四个步骤：

(B-1)对天线振子的输出信号进行解调、低通滤波后，分别得到基带接收信号和

(B-2)分别将接收探测信号和与标准的本地伪随机序列a(t)做滑动相关，即可得到观测冲激响应和表达式分别如下

其中，表示扩频增益，N′_m(t)是与本地伪随机序列相关后的噪声信号。

(B-3)分别将观测冲激响应和离散化。选择值最大的L个峰值点，分别表示为和将组成观测冲激响应向量将组成观测冲激响应向量 和的数学表达式分别为

其中，表示噪声矢量。

(B-4)对于子阵Z_α和子阵Z_β上第1个至第M个接收天线振子的接收信号都执行以上三个步骤，分别得到它们的观测冲激响应矢量，表示为和将这些观测冲激响应矢量按照如下形式组织成观测冲激响应矩阵H_α和H_β

H_α是M行L列的矩阵，H_β是M行L列的矩阵，它们的数学表达式为

从式(16)和式(17)中看到，观测冲激响应矩阵的第m行代表了第m个振子的观测冲激响应向量，第l列表示第l条径在整个阵列上的观测冲激响应，称为阵列冲激响应向量。

(C)估计波达角。

空间信道包含L条多径，一次估计其中的p条，1≤p≤M，具体的步骤如下。

(C-1)从观测冲激响应矩阵H_α中选择p列，将这p列编号为p₁,…,p_p，并构成矩阵 的数学表达式为

其中，表示噪声向量，表示噪声矩阵，A是导向矢量矩阵，定义如下

矩阵其中diag[]表示构造对角矩阵。

再从观测冲激响应矩阵H_β中同样的位置选择p列，同样将这p列编号为p₁,…,p_p，并构成矩阵 的表达式如下

其中，矩阵和矩阵均是M×p矩阵。如果所选择的p列所对应的p条径的入射信号角度互不相同，则证明矩阵和矩阵的秩均为p。

求取矩阵和矩阵的秩，如果结果均等于p，则说明所选择的p条径的入射信号波达角互不相同，则继续执行(C-2)步。如果结果不等于p，则说明所选择的p条径中至少有两条径对应的入射信号的波达角是相同的。此时将p的值减去1，并重新执行本子步骤。证明当p最终减至1时，矩阵和矩阵的秩一定等于1，且p＝1时本方法仍然成立。因此本方法必成立。

(C-2)分别按照式(18)和式(19)求取和的协方差矩阵和

其中，和分别表示和的共轭转置矩阵。矩阵和都是M×M矩阵。证明，矩阵和的秩均是p。

(C-3)对矩阵进行特征分解，并将特征值按照从大到小的顺序排列，表示为求取这些特征值对应的特征向量，表示为按照式(20)构建矩阵U_α。

同样地对矩阵进行特征分解，并将特征值按照从大到小的顺序排列，表示为求取这些特征值对应的特征向量，表示为按照式(21)构建矩阵U_β。

矩阵U_α和U_β均是M×p的矩阵。证明，矩阵U_α和U_β的秩均为p。

(C-4)按照式(22)求取矩阵Ψ。

Ψ＝(U_α)⁺U_β (25)

其中，()⁺表示求矩阵的广义逆。矩阵Ψ是p×p的矩阵。

对矩阵Ψ求取特征值，分别表示为证明即是矩阵Φ的对角元素。根据即可求出所选择的p列所对应的p条径的信号波达角。求取的方法见式(26)。

其中，argsin()表示求反正弦，angle()表示求复数的辐角。完成p条径的信号波达角估计。

(C-5)转至(C-1)步骤，再选择p条径进行估计，直至所有的L条径都完成估计。

本发明的有益效果是：该方法首先在两个不同子阵的接收信号中分别计算观测响应矩阵，即对每一个子阵的天线振子的输出信号进行解调和低通滤波后，得到基带接收信号，再将基带接收信号与伪随机序列做滑动相关，得到两个子阵的观测冲激响应，通过对两个子阵的观测冲激响应进行离散化，抽取出多径，得到两个子阵的观测冲激响应矢量，分别将两个子阵的所有天线振子的观测冲激响应矢量整合得到两个子阵的观测响应矩阵；最后利用两个子阵接收信号协方差矩阵的旋转不变性估计多径信号的波达角。由于本发明方法中，一个估计子过程只处理一部分多径信号的波达角估计，大量的多径信号可以分不同的估计子过程进行处理。因此可估计的多径波达角数量不受接收天线阵列规模的限制。由于本发明方法不直接使用采样信号，而是根据采样信号处理得到的观测冲激响应矩阵进行估计，因此涉及的数据量较小。此外，本发明方法不涉及迭代或搜索，因此运算量较小。故本发明方法简单、运算复杂度低、运算结果精度高，适用于对大量多径信号波达角进行估计，突破了现有ESPRIT方法对于天线振元的限制，实用性强。

以下结合具体实施方式详细说明本发明。

具体实施方式

本发明基于阵列响应旋转不变性的多径信号波达角估计方法具体步骤如下：

(A)定义信号；(B)计算观测冲激响应矩阵；(C)估计波达角。

(A)定义信号。

选择长度为X的伪随机序列作为基带探测信号a(t)，其表达式为

其中，X表示伪随机序列的长度，t表示时间，表示宽度为T_b的矩形脉冲信号。K个PN序列组成一个探测帧u(t)，其表达式为

其中T_p＝XT_b。探测帧u(t)是本方法中的基本探测信号，该探测帧调制后经天线发射出去。

对于包含了L条多径的传播环境，其空间信道冲激响应模型h(t)表示如下

其中，L为传播环境中的多径条数，表示第l条多径的信道复响应，是一个复常数，τ_l是第l条多径的时延值，δ(t)表示冲激函数。

接收天线阵列是一维线形天线阵列，包含M+1个天线振子(M＞1)。这M+1个天线振子等距离排列，振子间距表示为d，并且天线振子的方向图都相同。从几何上将天线阵列分为完全相同的两个子阵：子阵Z_α和子阵Z_β。其中，子阵Z_α由原天线阵列中的第1～M号天线振子组成，子阵Z_β由原天线阵列中第2～M+1号天线振子组成。两个子阵都包括有M个天线振子，因此把两个子阵中的天线振子都重新编号为1～M。子阵Z_α上的第m个天线振子和子阵Z_β上的第m个天线振子的距离表示为Δ。

对于子阵中的第m个天线振子，则其在的来波方向θ上的复响应s(θ)为

其中e表示自然底数，j表示虚数，表示射频信号的波长，θ表示来波方向。

对于子阵Z_α上的第m个天线振子，按照式(5)生成射频接收信号

其中N′_m(t)是服从高斯分布的白噪声，θ_l表示第l条径的入射角度，和τ_l为多径的复响应和时延，u′(t)为调制后的探测帧。对子阵Z_α上的所有天线振子，都生成相应的接收信号。

对于子阵Z_β上的第m个天线振子，按照式(6)生成射频接收信号

其中N′_m(t)是服从高斯分布的白噪声，θ_l表示第l条径的入射角度，和τ_l为多径的复响应和时延，u′(t)为调制后的探测帧。对于子阵Z_β上的所有天线振子，都生成相应的接收信号。

子阵Z_α上的第m个天线振子输出的接收信号经解调后表示为

其中，N_m(t)是第m个天线振子所接收到的噪声信号，θ_l是第l条径入射信号的波达角，s_m(θ_l)是第m个天线振子对波达角为θ_l的入射信号的复响应。

子阵Z_β上的第m个天线振子输出的接收信号经解调后表示为

其中，e是自然底数，j表示虚数。由于子阵Z_β上的第m个天线振子相对于子阵Z_α上的第m个天线振子存在距离Δ，这种距离导致了来波方向上的波程差。进而对于入射角为θ_l的多径信号，子阵Z_β的接收信号相对于子阵Z_α的接收信号出现了相移ξ(θ_l)。ξ(θ_l)的表达式为

其中，表示射频信号的波长。

(B)计算观测冲激响应矩阵。

对于子阵Z_α和子阵Z_β上的第m个天线振子，计算观测冲激响应矩阵分为如下四个步骤：

(B-1)对天线振子的输出信号进行解调、低通滤波后，分别得到基带接收信号和

(B-2)分别将接收探测信号和与标准的本地伪随机序列a(t)做滑动相关，即可得到观测冲激响应和表达式分别如下

其中，表示扩频增益，N′_m(t)是与本地伪随机序列相关后的噪声信号。

(B-3)分别将观测冲激响应和离散化。选择值最大的L个峰值点，分别表示为和将组成观测冲激响应向量将组成观测冲激响应向量 和的数学表达式分别为

其中，表示噪声矢量。

(B-4)对于子阵Z_α和子阵Z_β上第1个至第M个接收天线振子的接收信号都执行以上三个步骤，分别得到它们的观测冲激响应矢量，表示为和将这些观测冲激响应矢量按照如下形式组织成观测冲激响应矩阵H_α和H_β

H_α是M行L列的矩阵，H_β是M行L列的矩阵，它们的数学表达式为

从式(16)和式(17)中看到，观测冲激响应矩阵的第m行代表了第m个振子的观测冲激响应向量，第l列表示第l条径在整个阵列上的观测冲激响应，称为阵列冲激响应向量。

(C)估计波达角。

空间信道包含L条多径，本方法一次估计其中的p条(1≤p≤M)，具体的步骤如下(其余的多径在第二次、第三次…..中估计出，如C-5所述)。

(C-1)从观测冲激响应矩阵H_α中选择p列(1≤p≤M)，将这p列编号为p₁,…,p_p，并构成矩阵 的数学表达式为

其中，表示噪声向量，表示噪声矩阵，A是导向矢量矩阵，定义如下

矩阵其中diag[]表示构造对角矩阵。

再从观测冲激响应矩阵H_β中同样的位置选择p列，同样将这p列编号为p₁,…,p_p，并构成矩阵 的表达式如下

其中矩阵和矩阵均是M×p矩阵。如果所选择的p列所对应的p条径的入射信号角度互不相同，则证明矩阵和矩阵的秩均为p。

求取矩阵和矩阵的秩，如果结果均等于p，则说明所选择的p条径的入射信号波达角互不相同，则继续执行(C-2)步。如果结果不等于p，则说明所选择的p条径中至少有两条径对应的入射信号的波达角是相同的。此时将p的值减去1，并重新执行本子步骤。证明当p最终减至1时，矩阵和矩阵的秩一定等于1，且p＝1时本方法仍然成立。因此本方法必成立。

(C-2)分别按照式(18)和式(19)求取和的协方差矩阵和

其中，和分别表示和的共轭转置矩阵。矩阵和都是M×M矩阵。证明，矩阵和的秩均是p。

(C-3)对矩阵进行特征分解，并将特征值按照从大到小的顺序排列，表示为求取这些特征值对应的特征向量，表示为按照式(20)构建矩阵U_α。

同样地对矩阵进行特征分解，并将特征值按照从大到小的顺序排列，表示为求取这些特征值对应的特征向量，表示为按照式(21)构建矩阵U_β。

矩阵U_α和U_β均是M×p的矩阵。证明，矩阵U_α和U_β的秩均为p。

(C-4)按照式(22)求取矩阵Ψ。

Ψ＝(U_α)⁺U_β (25)

其中，()⁺表示求矩阵的广义逆。矩阵Ψ是p×p的矩阵。

对矩阵Ψ求取特征值，分别表示为证明即是矩阵Φ的对角元素。根据即可求出所选择的p列所对应的p条径的信号波达角。求取的方法见式(26)。

其中，argsin()表示求反正弦，angle()表示求复数的辐角。这样就完成了p条径的信号波达角估计。

(C-5)转至(C-1)步骤，再选择p条径进行估计，直至所有的L条径都完成估计。以上(C-1)至(C-4)步骤称为一个估计子过程。需要说明的是，不同估计子过程的p值不一定相同。

以下介绍一种具体的仿真实施方式。

(a)定义信号。

按照以下步骤定义信号。

(a-1)生成探测帧信号。使用长度为1023的m序列作为伪随机序列，基带探测信号a(t)的码速率为100兆比特/秒，也即式(1)中的T_b＝10ns，其中ns表示纳秒。一个探测帧u(t)由两个伪随机序列连接组成，也即式(2)中K＝2。探测帧通过BPSK调制，载波频率为2.5GHz。调制后的探测帧表示为u′(t)。

(a-2)生成多径信息。假设环境中包含20条多径。这20条多径的时延、复响应和入射信号角度对发送端和接收端是未知的。为进行仿真，将这20条径的时延和复响应按照表1进行设置。

表1不同多径的时延和复响应

序号	复响应	时延(ns)		序号	复响应	时延(ns)
							1	ξ	100		11	0.80ξ	200
2	0.98ξ	110		12	0.78ξ	210
							3	0.96ξ	120		13	0.76ξ	220
4	0.94ξ	130		14	0.74ξ	230
							5	0.92ξ	140		15	0.72ξ	240
6	0.90ξ	150		16	0.70ξ	250
							7	0.88ξ	160		17	0.68ξ	260
8	0.86ξ	170		18	0.66ξ	270

9	0.84ξ	180		19	0.64ξ	280
							10	0.82ξ	190		20	0.62ξ	290

其中，ξ是一个复常数，可自由设置。时延的单位是纳秒。这20条多径的入射角可随机生成，角度分辨率为1度，取值范围为-90～90度。

(a-3)生成接收天线阵列。接收天线是线形天线阵列，包含6枚天线振子，相邻天线振子的距离表示为d，d等于射频信号波长的一半，也即6厘米。则子阵Z_α由1～5号天线振子组成，子阵Z_β由2～6号天线振子组成。将子阵Z_α和子阵Z_β上的天线振子都重新编号为1～5号。则两个子阵上相同编号振子的距离为Δ，易知本实施例中Δ＝d＝6厘米。

(a-4)生成子阵Z_α和子阵Z_β在不同波达方向上的导向矢量。由于两个子阵的振子和阵列结构相同，因此导向矢量也相同。对于子阵中的第m个天线振子，则其在的来波方向θ上的复响应s(θ)为

其中e表示自然底数，j表示虚数，表示射频信号的波长，在本例中为12厘米(0.12米)。对所有20条径的入射信号，都按照式(4)计算出子阵中1～5号天线振子在其来波θ₁,…，θ₂₀上的复响应s₁(θ₁),…,s₁(θ₂₀),s₂(θ₁),…,s₅(θ₂₀)。

(a-5)生成子阵相对相移。所谓子阵相对相移，是指对于入射角为θ_l的第l条径入射信号，子阵Z_β上的振子相对于子阵Z_α上同一个编号振子所输出信号的相移ξ(θ_l)。按照式(9)计算出所有径的入射信号的子阵相对相移ξ(θ₁),…,ξ(θ₂₀)。

(a-6)生成接收信号。对于子阵Z_α上的第m个天线振子，按照式(5)生成射频接收信号

其中N′_m(t)是服从高斯分布的白噪声，信噪比设置为0dB。和τ_l即为表1中所示的径的复响应和时延。对子阵Z_α上的所有5枚天线振子，都生成相应的接收信号。

对于子阵Z_β上的第m个天线振子，按照式(6)生成射频接收信号

对于子阵Z_β上的所有5枚天线振子，都生成相应的接收信号。

(b)计算观测冲激响应矩阵。

按照如下步骤计算观测冲激响应矩阵。

(b-1)设置天线振子序号m为1。

(b-2)分别对射频接收信号和进行BPSK解调、低通滤波(滤波器带宽100MHz)，得到基带探测帧和

(b-3)将基带探测帧和基带探测信号a(t)做滑动相关，得到观测冲激响应求得的最大值，表示为设置阀值Thr_α，其值为

再将基带探测帧和基带探测信号a(t)做滑动相关，得到观测冲激响应求得的最大值，表示为设置阀值Thr_β，其值为

(b-4)从t＝0开始，找出满足的20个峰值点，这20个峰值点的值即为离散化的观测冲激响应，将这20个点的值按照式(12)的形式组成一个行向量，即为观测冲激响应矢量

再从t＝0开始，找出满足的20个峰值点，这20个峰值点的值即为离散化的观测冲激响应，将这20个点的值按照式(13)的形式组成一个行向量，即为观测冲激响应矢量

(b-5)天线振子序号m加1，转回到(b-2)步骤执行，直至所有5个天线振子的观测冲激响应矢量都求取完成

(b-6)完成以上子步骤后，按照式(14)所示把子阵Z_α上振子的观测冲激响应矢量组合成为一个观测冲激响应矩阵，表示为H_α。再按照式(15)所示把子阵Z_β上振子的观测冲激响应矢量组合成为一个观测冲激响应矩阵，表示为H_β。本例中H_α和H_β均是5×20的矩阵。

(c)估计波达角。

按照如下步骤进行多径信号波达角的估计。

(c-1)以η表示估计子过程的序号，设置η为1。设置一次估计的多径信号波达角数量p设置为3。

(c-2)取出观测冲激响应矩阵H_α的第(η-1)×p+1至第η×p列。如果η×p大于20，则取观测冲激响应矩阵H_α的第(η-1)×p+1至第20列，并更新p的值为20-(η-1)×p。

将这p列编号为1～p。按照式(18)构成矩阵再取出观测冲激响应矩阵H_β的第((η-1)×p+1)至第(η×p)列，将这p列同样编号为1～p。按照式(20)构成矩阵求取矩阵的秩，表示为R_α，再求取矩阵的秩，表示为R_β。如果R_α或R_β不等于p，则将p的值减1，重新执行本子步骤。

(c-3)分别按照式(21)和式(22)计算矩阵和的协方差矩阵和本例中，矩阵和均是5×5的矩阵。

(c-4)对矩阵进行特征分解，并将特征值按照从大到小的顺序排列，表示为再求特征值对应的特征向量，表示为按照式(23)构建矩阵U_α。

同样地对矩阵进行特征分解，并将特征值按照从大到小的顺序排列，表示为再求特征值对应的特征向量，表示为按照式(24)构建矩阵U_β

(c-5)按照式(25)计算矩阵Ψ，求取矩阵Ψ的特征值，表示为按照式(26)分别求出这p条径的信号波达角，分别表示为并保存。

(c-6)估计子过程序号η的值加1，设置p的值为3，并转回至(c-2)步骤执行，直至所有径的信号波达角都完成估计。

将所有估计子过程输出的信号波达角估计值组成集合，即为所求20条径的入射信号波达角的估计值集合。

Claims (1)

1.一种基于阵列响应旋转不变性的多径信号波达角估计方法，其特征在于包括以下步骤：

(A)定义信号；

选择长度为X的伪随机序列作为基带探测信号a(t)，其表达式为

$a\left(t\right)=\underset{n=0}{\overset{X-1}{\Σ}}{b}_n{\mathrm{rect}}_{T_b}(t-{\mathrm{nT}}_b),b_n\∈\{+1,-1\}---\left(1\right)$

其中，X表示伪随机序列的长度，t表示时间，表示宽度为T_b的矩形脉冲信号；K个PN序列组成一个探测帧u(t)，其表达式为

$u\left(t\right)=\underset{k=0}{\overset{K-1}{\Σ}}a(t-{\mathrm{kT}}_p)---\left(2\right)$

其中，T_p＝XT_b；探测帧u(t)是基本探测信号，该探测帧调制后经天线发射出去；

对于包含了L条多径的传播环境，其空间信道冲激响应模型h(t)表示如下

$h\left(t\right)=\underset{l=1}{\overset{L}{\Σ}}{h}_l^{CH}\mathrm{\δ}(t-{\mathrm{\τ}}_l)---\left(3\right)$

其中，L为传播环境中的多径条数，L＝20；表示第l条多径的信道复响应，是一个复常数，τ_l是第l条多径的时延值，δ(t)表示冲激函数；

接收天线阵列是一维线形天线阵列，包含M+1个天线振子，其中，M>1；这M+1个天线振子等距离排列，振子间距表示为d，并且天线振子的方向图都相同；从几何上将天线阵列分为完全相同的两个子阵：子阵Z_α和子阵Z_β；其中，子阵Z_α由原天线阵列中的第1～M号天线振子组成，子阵Z_β由原天线阵列中第2～M+1号天线振子组成；两个子阵都包括有M个天线振子，因此把两个子阵中的天线振子都重新编号为1～M；子阵Z_α上的第m个天线振子和子阵Z_β上的第m个天线振子的距离表示为Δ；

对于子阵中的第m个天线振子，则其在的来波方向θ上的复响应s(θ)为

其中，e表示自然底数，j表示虚数，表示射频信号的波长，θ表示来波方向；

对于子阵Z_α上的第m个天线振子，按照式(5)生成射频接收信号

$y_{{\mathrm{\α}}_m}^{\′}\left(t\right)=\underset{l=1}{\overset{20}{\Σ}}{h}_l^{CH}{s}_m\left({\mathrm{\θ}}_l\right)u^{\′}(t-{\mathrm{\τ}}_l)+N_m^{\′}\left(t\right)---\left(5\right)$ 其中，N′_m(t)是服从高斯分布的白噪声，θ_l表示第l条径的入射角度，和τ_l为多径的复响应和时延，u′(t)为调制后的探测帧；对子阵Z_α上的所有天线振子，都生成相应的接收信号；

对于子阵Z_β上的第m个天线振子，按照式(6)生成射频接收信号

$y_{{\mathrm{\β}}_m}^{\′}\left(t\right)=\underset{l=1}{\overset{20}{\Σ}}{h}_l^{CH}{s}_m\left({\mathrm{\θ}}_l\right)e^{j\mathrm{\ξ}\left({\mathrm{\θ}}_l\right)}{u}^{\′}(t-{\mathrm{\τ}}_l)+N_m^{\′}\left(t\right)---\left(6\right)$

其中，N′_m(t)是服从高斯分布的白噪声，θ_l表示第l条径的入射角度，和τ_l为多径的复响应和时延，u′(t)为调制后的探测帧；对于子阵Z_β上的所有天线振子，都生成相应的接收信号；

子阵Z_α上的第m个天线振子输出的接收信号经解调后表示为

$y_{{\mathrm{\α}}_m}\left(t\right)=\underset{l=1}{\overset{L}{\Σ}}{h}_l^{CH}{s}_m\left({\mathrm{\θ}}_l\right)u(t-{\mathrm{\τ}}_l)+N_m\left(t\right)---\left(7\right)$

其中，N_m(t)是第m个天线振子所接收到的噪声信号，θ_l是第l条径入射信号的波达角，s_m(θ_l)是第m个天线振子对波达角为θ_l的入射信号的复响应；

子阵Z_β上的第m个天线振子输出的接收信号经解调后表示为

$y_{{\mathrm{\β}}_m}\left(t\right)=\underset{l=1}{\overset{L}{\Σ}}{h}_l^{CH}{s}_m\left({\mathrm{\θ}}_l\right)e^{j\mathrm{\ξ}\left({\mathrm{\θ}}_l\right)}u(t-{\mathrm{\τ}}_l)+N_m\left(t\right)---\left(8\right)$

其中，e是自然底数，j表示虚数；由于子阵Z_β上的第m个天线振子相对于子阵Z_α上的第m个天线振子存在距离Δ，这种距离导致了来波方向上的波程差；进而对于入射角为θ_l的多径信号，子阵Z_β的接收信号相对于子阵Z_α的接收信号出现了相移ξ(θ_l)；ξ(θ_l)的表达式为

其中，表示射频信号的波长；

(B)计算观测冲激响应矩阵；

对于子阵Z_α和子阵Z_β上的第m个天线振子，计算观测冲激响应矩阵分为如下四个步骤：

(B-1)对天线振子的输出信号进行解调、低通滤波后，分别得到基带接收信号和

(B-2)分别将接收探测信号和与标准的本地伪随机序列a(t)做滑动相关，即可得到观测冲激响应和表达式分别如下

${\tilde{h}}_{{\mathrm{\α}}_m}\left(t\right)=\underset{l=1}{\overset{L}{\Σ}}{\mathrm{\ρh}}_l^{CH}{s}_m\left({\mathrm{\θ}}_l\right)\mathrm{\δ}(t-{\mathrm{\τ}}_l)+N_m^{\′}\left(t\right)---\left(10\right)$

${\tilde{h}}_{{\mathrm{\β}}_m}\left(t\right)=\underset{l=1}{\overset{L}{\Σ}}{\mathrm{\ρh}}_l^{CH}{s}_m\left({\mathrm{\θ}}_l\right)e^{j\mathrm{\ξ}\left({\mathrm{\θ}}_l\right)}\mathrm{\δ}(t-{\mathrm{\τ}}_l)+N_m^{\′}\left(t\right)---\left(11\right)$

其中，表示扩频增益，N′_m(t)是与本地伪随机序列相关后的噪声信号；

(B-3)分别将观测冲激响应和离散化；选择值最大的L个峰值点，分别表示为和将组成观测冲激响应向量将组成观测冲激响应向量 和的数学表达式分别为

$\begin{array}{c}{\tilde{a}}_{{\mathrm{\α}}_m}=\[\begin{array}{ccc}{\tilde{h}}_{{\mathrm{\α}}_m}\left({\mathrm{\τ}}_1\right)& \mathrm{...}& {\tilde{h}}_{{\mathrm{\α}}_m}\left({\mathrm{\τ}}_L\right)\end{array}\]\\ =\mathrm{\ρ}\[\begin{array}{ccc}{h}_1^{CH}{s}_m\left({\mathrm{\θ}}_1\right)& \mathrm{...}& h_L^{CH}{s}_m\left({\mathrm{\θ}}_L\right)\end{array}\]+{\tilde{N}}_m^{\′}\end{array}---\left(12\right)$

$\begin{array}{c}{\tilde{a}}_{{\mathrm{\β}}_m}=\[\begin{array}{ccc}{\tilde{h}}_{{\mathrm{\β}}_m}\left({\mathrm{\τ}}_1\right)& \mathrm{...}& {\tilde{h}}_{{\mathrm{\β}}_m}\left({\mathrm{\τ}}_L\right)\end{array}\]\\ =\mathrm{\ρ}\[\begin{array}{ccc}{h}_1^{CH}{s}_m\left({\mathrm{\θ}}_1\right)e^{j\mathrm{\ξ}\left({\mathrm{\θ}}_1\right)}& \mathrm{...}& h_L^{CH}{s}_m\left({\mathrm{\θ}}_L\right)e^{j\mathrm{\ξ}\left({\mathrm{\θ}}_L\right)}\end{array}\]+{\tilde{N}}_m^{\′}\end{array}---\left(13\right)$

其中，表示噪声矢量；

(B-4)对于子阵Z_α和子阵Z_β上第1个至第M个接收天线振子的接收信号都执行以上三个步骤，分别得到它们的观测冲激响应矢量，表示为和将这些观测冲激响应矢量按照如下形式组织成观测冲激响应矩阵H_α和H_β

$H_{\mathrm{\α}}=\left[\begin{array}{c}{\tilde{a}}_{{\mathrm{\α}}_1}\\ .\\ .\\ .\\ {\tilde{a}}_{{\mathrm{\α}}_M}\end{array}\right]---\left(14\right)$

$H_{\mathrm{\β}}=\left[\begin{array}{c}{\tilde{a}}_{{\mathrm{\β}}_1}\\ .\\ .\\ .\\ {\tilde{a}}_{{\mathrm{\β}}_M}\end{array}\right]---\left(15\right)$

H_α是M行L列的矩阵，H_β是M行L列的矩阵，它们的数学表达式为

从式(16)和式(17)中看到，观测冲激响应矩阵的第m行代表了第m个振子的观测冲激响应向量，第l列表示第l条径在整个阵列上的观测冲激响应，称为阵列冲激响应向量；

(C)估计波达角；

空间信道包含L条多径，一次估计其中的p条，1≤p≤M，具体的步骤如下；

(C-1)从观测冲激响应矩阵H_α中选择p列，将这p列编号为p₁,…,p_p，并构成矩阵 的数学表达式为

其中，表示噪声向量，表示噪声矩阵，A是导向矢量矩阵，定义如下

矩阵其中diag[]表示构造对角矩阵；

再从观测冲激响应矩阵H_β中同样的位置选择p列，同样将这p列编号为p₁,…,p_p，并构成矩阵 的表达式如下

其中，矩阵和矩阵均是M×p矩阵；如果所选择的p列所对应的p条径的入射信号角度互不相同，则证明矩阵和矩阵的秩均为p；

求取矩阵和矩阵的秩，如果结果均等于p，则说明所选择的p条径的入射信号波达角互不相同，则继续执行(C-2)步；如果结果不等于p，则说明所选择的p条径中至少有两条径对应的入射信号的波达角是相同的；此时将p的值减去1，并重新执行本子步骤；证明当p最终减至1时，矩阵和矩阵的秩一定等于1，且p＝1时本方法仍然成立；因此本方法必成立；

(C-2)分别按照式(21)和式(22)求取和的协方差矩阵和

${\tilde{C}}_{\mathrm{\α}}={\hat{H}}_{\mathrm{\α}}{\hat{H}}_{\mathrm{\α}}^H---\left(21\right)$

${\tilde{C}}_{\mathrm{\β}}={\hat{H}}_{\mathrm{\β}}{\hat{H}}_{\mathrm{\β}}^H---\left(22\right)$

其中，和分别表示和的共轭转置矩阵；矩阵和都是M×M矩阵；证明，矩阵和的秩均是p；

(C-3)对矩阵进行特征分解，并将特征值按照从大到小的顺序排列，表示为求取这些特征值对应的特征向量，表示为按照式(23)构建矩阵U_α；

$U_{\mathrm{\α}}=\left[\begin{array}{ccc}{u}_{{\mathrm{\α}}_1}& \mathrm{...}& u_{{\mathrm{\α}}_p}\end{array}\right]---\left(23\right)$

同样地对矩阵进行特征分解，并将特征值按照从大到小的顺序排列，表示为求取这些特征值对应的特征向量，表示为按照式(24)构建矩阵U_β；

$U_{\mathrm{\β}}=\left[\begin{array}{ccc}{u}_{{\mathrm{\β}}_1}& \mathrm{...}& u_{{\mathrm{\β}}_p}\end{array}\right]---\left(24\right)$

矩阵U_α和U_β均是M×p的矩阵；证明，矩阵U_α和U_β的秩均为p；

(C-4)按照式(25)求取矩阵Ψ；

Ψ＝(U_α)⁺U_β (25)

其中，()⁺表示求矩阵的广义逆；矩阵Ψ是p×p的矩阵；

对矩阵Ψ求取特征值，分别表示为证明即是矩阵Φ的对角元素；根据即可求出所选择的p列所对应的p条径的信号波达角；求取的方法见式(26)；

其中，argsin()表示求反正弦，angle()表示求复数的辐角；完成p条径的信号波达角估计；

(C-5)转至(C-1)步骤，再选择p条径进行估计，直至所有的L条径都完成估计。