二维DOA估计方法、系统及计算机可读介质与流程

二维doa估计方法、系统及计算机可读介质
技术领域
1.本发明主要涉及天线阵列信号处理技术领域，具体地涉及一种二维doa估计方法、系统及计算机可读介质。


背景技术：

2.码分多址(code division multiple access，cdma)通信体制具有抗干扰、抗截获能力强，通信隐蔽，易于实现多址通信、语音激活、软容量和软切换等特点，在低轨卫星通信系统得到广泛应用。在cdma通信系统中，天线阵列接收的信号数往往远大于阵元数目，导致传统的波达方向(direction of arrival，doa)估计算法难以直接应用。即使在通信用户较少的情况下，由于扩频信号的信噪比很低，直接对阵元接收信号进行doa估计，效果也不甚理想。
3.针对上述问题，在通信用户扩频序列已知的情况下，通常使用码匹配滤波器对天线阵列的接收信号进行解扩预处理。由于扩频序列良好的自相关与互相关特性，对某一通信用户信号解扩后，其他通信用户信号可近似为噪声，解扩信号可看作该通信用户信号与噪声的合成信号，且信噪比得到提升。对于各通信用户解扩后的等效单源信号，可使用常规doa估计方法逐用户进行角度估计。常用的高分辨率doa估计算法如多信号分类(multiple signal classification，music)算法、基于旋转不变技术的信号参数估计(estimating signal parameters via rotational invariance techniques，esprit)算法等。
4.在对解扩信号逐用户进行doa估计时，由于常规doa估计算法需要谱峰搜索、特征分解等操作，在通信用户数量较多时，存在系统计算开销较为庞大的问题。寻找能够在cdma通信系统中具有稳定测向性能、低计算复杂度的doa估计算法，在理论与工程应用中都具有重要的研究价值。


技术实现要素：

5.本技术所要解决的技术问题是提供一种二维doa估计方法、系统及计算机可读介质，本技术的估计方法测向稳定，在保证测向精度的同时，计算复杂度大大降低。
6.本技术为解决上述技术问题而采用的技术方案是一种二维doa估计方法，该二维doa估计方法应用于cdma通信体制下，包括：将天线阵列的接收信号进行解扩预处理，获得解扩信号，将天线阵列划分为第一子阵a0、第二子阵a1和第三子阵a2，计算第一子阵a0解扩信号的自相关矩阵第一子阵a0解扩信号与第二子阵a1解扩信号的互相关矩阵第一子阵a0解扩信号与第三子阵a2解扩信号的互相关矩阵剔除自相关矩阵和互相关矩阵的主对角线元素与左下副对角线元素，并将自相关矩阵和互相关矩阵的剩余元素分别重构为列矢量形式，获得重构矢量和重构矢量以及，剔除自相关矩阵和互相关矩阵的主对角线元素与左下角((n-1)(n-2))阶矩阵的主对角线元素，并将自
相关矩阵和互相关矩阵的剩余元素分别重构为列矢量形式，获得重构矢量和重构矢量根据重构矢量和重构矢量计算估计矩阵根据重构矢量和重构矢量计算估计矩阵根据估计矩阵和估计矩阵使用元素平均法、相位平均法、最小二乘法中的任意一种方法计算相位估计和相位估计相位估计和相位估计包括解扩信号的doa信息，doa信息包括俯仰角和/或方位角；根据相位估计和相位估计分别计算天线阵列中x轴方向的方向余弦估计和y轴方向的方向余弦估计以及根据方向余弦估计和方向余弦估计计算俯仰角估计和/或方位角估计
7.在本技术的一实施例中，将天线阵列的接收信号进行解扩预处理，获得解扩信号的步骤包括：采用下面的公式计算解扩信号：
[0008][0009][0010]
其中，n表示第n个采样点，x(n)表示接收信号，k表示用户数量为k个，n(n)表示阵列噪声矢量，s表示用户的符号信息数，σk表示信号幅值，dk(i)∈{-1,1}(i＝0,1,...,s-1)表示第k个用户发送的符号信息，t表示符号周期，ak表示第k个用户信号对应的导向矢量；ck(n')∈{-1,1},n'＝0,1,...,g-1为分配给第k个用户的扩频序列，包括伪随机序列，g为扩频增益，即ck(n-it)表示ck(n')在第i个符号周期的延拓；
[0011]
yk(i)表示解扩信号，解扩信号具有等效单源特性，k表示第k个用户，k'表示第k个用户以外的其他用户，表示第k个用户的扩频序列与其他用户扩频序列的互相关函数；等效单源信号sk(i)＝gσ
kdk
(i)，等效噪声
[0012]
在本技术的一实施例中，天线阵列为n
×
n维的阵列，第一子阵a0、第二子阵a1和第三子阵a2均为(n-1)
×
(n-1)维的方阵，将天线阵列划分为第一子阵a0、第二子阵a1和第三子阵a2的步骤包括：将天线阵列剔除第n行阵元与第n列阵元获得第一子阵a0；将天线阵列剔除第n行阵元与第1列阵元获得第二子阵a1；将天线阵列剔除第1行阵元与第n列阵元获得第三子阵a2。
[0013]
在本技术的一实施例中，采用下面的公式分别计算自相关矩阵互相关矩阵和互相关矩阵
[0014]
[0015][0016][0017]
其中，s表示用户的符号信息数，x0(i)表示第一子阵a0解扩信号，表示第一子阵a0解扩信号的共轭转置，表示第二子阵a1解扩信号的共轭转置，表示第三子阵a2解扩信号的共轭转置。
[0018]
在本技术的一实施例中，根据重构矢量和重构矢量计算估计矩阵的步骤包括：采用下面的公式计算估计矩阵
[0019][0020]
其中，φ1为第二子阵a1的旋转矩阵，为全1矢量，[
·
]i表示矩阵第i行元素。
[0021]
在本技术的一实施例中，根据重构矢量和重构矢量计算估计矩阵的步骤包括：采用下面的公式计算估计矩阵
[0022][0023]
其中，φ2为第三子阵a2的旋转矩阵，为全1矢量，[
·
]i表示矩阵第i行元素。
[0024]
在本技术的一实施例中，根据估计矩阵使用元素平均法计算相位估计的步骤包括：采用下面的公式进行计算：
[0025][0026][0027]
其中，表示使用元素平均法得到的第二子阵a1的旋转矩阵估计值，n表示天线阵列的维度，表示使用元素平均法得到的相位估计，angle(
·
)函数用于求解相位角。
[0028]
在本技术的一实施例中，根据估计矩阵使用相位平均法计算相位估计的步骤包括：采用下面的公式进行计算：
[0029][0030]
其中，表示使用相位平均法得到的相位估计，n表示天线阵列的维度，angle(
·
)函数用于求解相位角。
[0031]
在本技术的一实施例中，根据估计矩阵使用最小二乘法计算相位估计的步骤包括：采用下面的公式进行计算：
[0032][0033][0034][0035]
其中，表示第二子阵a1的旋转矩阵的共轭，表示重构矢量的伪逆，表示重构矢量的共轭转置，表示使用最小二乘法得到的相位估计，angle(
·
)函数用于求解相位角。
[0036]
在本技术的一实施例中，根据相位估计和相位估计分别计算天线阵列中x轴方向的方向余弦估计和y轴方向的方向余弦估计的步骤包括：采用下面的公式计算方向余弦估计和方向余弦估计
[0037][0038]
其中，λ表示接收信号的波长，d表示天线阵列的阵元间距。
[0039]
在本技术的一实施例中，根据方向余弦估计和方向余弦估计计算俯仰角估计和/或方位角估计的步骤包括：采用下面的公式计算俯仰角估计和/或方位角估计
[0040][0041][0042]
本技术为解决上述技术问题还提出一种二维doa估计系统，包括：天线阵列，天线阵列用于接收信号；存储器，用于存储可由处理器执行的指令；处理器，用于执行指令以实现如上的二维doa估计方法。
[0043]
本技术为解决上述技术问题还提出一种存储有计算机程序代码的计算机可读介质，计算机程序代码在由处理器执行时实现如上的二维doa估计方法。
[0044]
本技术的技术方案通过对接收信号进行解扩预处理后，解扩信号具有等效单源特性，在通信用户数量较多时，进行逐用户的doa估计，可以更好地适应天线阵列的过载状态；将天线阵列划分为三个子阵，根据子阵间的相关关系，使用点除或元素平均法、相位平均法、最小二乘法中的任意一种方法代替标准esprit算法中系统计算开销较大的矩阵特征分解操作，在保证测向精度的同时，计算复杂度大大降低；根据天线阵列中的方向余弦估计计算俯仰角估计与方位角估计，测向稳定。
附图说明
[0045]
为让本技术的上述目的、特征和优点能更明显易懂，以下结合附图对本技术的具体实施方式作详细说明，其中：
[0046]
图1是本技术一实施例的二维doa估计方法的示例性流程图；
[0047]
图2是本技术一实施例的二维doa估计方法的天线阵列示意图；
[0048]
图3是本技术一实施例的二维doa估计方法的相关矩阵重构示意图；
[0049]
图4是本技术一实施例的二维doa估计方法仿真结果的俯仰角估计误差与输入snr关系曲线示意图；
[0050]
图5是本技术一实施例的二维doa估计方法仿真结果的方位角估计误差与输入snr关系曲线示意图；
[0051]
图6是本技术一实施例的二维doa估计方法仿真结果的俯仰角估计误差与用户信息符号数关系曲线示意图；
[0052]
图7是本技术一实施例的二维doa估计方法仿真结果的方位角估计误差与用户信息符号数关系曲线示意图；
[0053]
图8是本技术一实施例的二维doa估计方法仿真结果的俯仰角估计误差与用户数关系曲线示意图；
[0054]
图9是本技术一实施例的二维doa估计方法仿真结果的方位角估计误差与用户数关系曲线示意图；
[0055]
图10是本技术一实施例的二维doa估计系统的系统框图。
具体实施方式
[0056]
为让本技术的上述目的、特征和优点能更明显易懂，以下结合附图对本技术的具体实施方式作详细说明。
[0057]
在下面的描述中阐述了很多具体细节以便于充分理解本技术，但是本技术还可以采用其它不同于在此描述的其它方式来实施，因此本技术不受下面公开的具体实施例的限制。
[0058]
如本技术和权利要求书中所示，除非上下文明确提示例外情形，“一”、“一个”、“一种”和/或“该”等词并非特指单数，也可包括复数。一般说来，术语“包括”与“包含”仅提示包括已明确标识的步骤和元素，而这些步骤和元素不构成一个排它性的罗列，方法或者设备也可能包含其他的步骤或元素。
[0059]
本技术中使用了流程图用来说明根据本技术的实施例的系统所执行的操作。应当理解的是，前面或下面操作不一定按照顺序来精确地执行。相反，可以按照倒序或同时处理各种步骤。同时，或将其他操作添加到这些过程中，或从这些过程移除某一步或数步操作。
[0060]
本技术提出一种二维doa估计方法，该方法应用于cdma通信体制下。cdma是在无线通讯上使用的技术，cdma允许所有通信用户同时使用全部频带，且把其他通信用户发出讯号视为杂讯，完全不必考虑讯号碰撞问题。doa估计是一种波达方向定位技术，doa估计通过处理接收到的回波信号，获取目标的距离信息和方位信息。示例性地，本技术的二维doa估计方法可以应用在基于cdma通信体制下的低轨卫星通信系统，在低轨卫星通信系统中，通信用户例如是卫星系统里的卫星。本技术的二维doa估计方法也可以推广到其他使用cdma通信的系统中，本技术对通信用户不作限制。
[0061]
图1是本技术一实施例的二维doa估计方法的示例性流程图，参考图1所示，该实施例的二维doa估计方法包括以下步骤：
[0062]
步骤s110：将天线阵列的接收信号进行解扩预处理，获得解扩信号，将天线阵列划分为第一子阵a0、第二子阵a1和第三子阵a2，计算第一子阵a0解扩信号的自相关矩阵第一子阵a0解扩信号与第二子阵a1解扩信号的互相关矩阵第一子阵a0解扩信号与第三子阵a2解扩信号的互相关矩阵
[0063]
步骤s120：剔除自相关矩阵和互相关矩阵的主对角线元素与左下副对角线元素，并将自相关矩阵和互相关矩阵的剩余元素分别重构为列矢量形式，获得重构矢量和重构矢量以及，剔除自相关矩阵和互相关矩阵的主对角线元素与左下角((n-1)(n-2))阶矩阵的主对角线元素，并将自相关矩阵和互相关矩阵的剩余元素分别重构为列矢量形式，获得重构矢量和重构矢量
[0064]
步骤s130：根据重构矢量和重构矢量计算估计矩阵
[0065]
步骤s140：根据重构矢量和重构矢量计算估计矩阵
[0066]
步骤s150：根据估计矩阵和估计矩阵使用元素平均法、相位平均法、最小二乘法中的任意一种方法计算相位估计和相位估计相位估计和相位估计包括解扩信号的doa信息，doa信息包括俯仰角和/或方位角。
[0067]
步骤s160：根据相位估计和相位估计分别计算天线阵列中x轴方向的方向余弦估计和y轴方向的方向余弦估计
[0068]
步骤s170：根据方向余弦估计和方向余弦估计计算俯仰角估计和/或方位角估计。
[0069]
下面详细说明上述的步骤s110至步骤s170：
[0070]
在步骤s110中，将天线阵列的接收信号进行解扩预处理，获得解扩信号，将天线阵列划分为第一子阵a0、第二子阵a1和第三子阵a2，计算第一子阵a0解扩信号的自相关矩阵第一子阵a0解扩信号与第二子阵a1解扩信号的互相关矩阵第一子阵a0解扩信号与第三子阵a2解扩信号的互相关矩阵
[0071]
考虑同步直接序列码分多址(direct sequence-code division multiple access，ds-cdma)系统中k个用户的情况，假设第k个用户发送的符号信息为dk(i)∈{-1,1}(i＝0,1,...,s-1)，信号经扩频操作并以功率发射。并且假设信号在理想信道中传播，忽略路径衰减及多径效应的影响，在接收端使用m元天线阵列进行接收，并对接收信号进行chip级采样。chip级采样即对扩频序列码字符号进行采样，扩频序列的每一位码字称为一个chip。
[0072]
基于上述假设，接下来具体说明对天线阵列的接收信号进行解扩预处理，获得解扩信号的步骤以及前述三个子阵的划分方法。
[0073]
在一些实施例中，将天线阵列的接收信号进行解扩预处理，获得解扩信号的步骤包括：采用下面的公式(1)和(2)计算解扩信号：
[0074][0075][0076]
其中，在上述公式(1)中，n表示第n个采样点，x(n)表示接收信号，k表示用户数量为k个，n(n)表示阵列噪声矢量，s表示用户的符号信息数，σk表示信号幅值，dk(i)∈{-1,1}(i＝0,1,...,s-1)表示第k个用户发送的符号信息，t表示符号周期，ak表示第k个用户信号对应的导向矢量；ck(n')∈{-1,1},n'＝0,1,...,g-1为分配给第k个用户的扩频序列，通常为伪随机序列，g为扩频增益，即ck(n-it)表示ck(n')在第i个符号周期的延拓，相当于将ck(n')向后移动了it个采样点数，得到ck(n-it)。
[0077]
将接收信号x(n)输入到m个相关器组中，对于第k个用户信号，解扩后输出的第i个符号级信号记为指信号维度为m
×
1维。
[0078]
上述公式(2)中，yk(i)表示解扩信号，解扩信号具有等效单源特性，k表示第k个用户，k'表示第k个用户以外的其他用户，表示第k个用户的扩频序列与其他用户扩频序列的互相关函数；等效单源信号sk(i)＝gσ
kdk
(i)，等效噪声
[0079]
可以证明，当g
→
∞时，γ
kk'
(i)近似于高斯概率分布，第k个用户之外的信号通过码匹配滤波器后，可以近似当作噪声处理。因此，可将yk(i)看作单源信号sk(i)与噪声的组合，后续可使用常规doa估计算法对第k个用户信号的来波方向进行估计。其中，单源信号指的是使用第k个用户的扩频序列对含有多个用户信号的接收信号进行解扩后，解扩信号可近似看作只含有第k个用户信号，而其他用户信号近似当作噪声处理。
[0080]
通过对接收信号进行解扩预处理后，解扩信号具有等效单源特性，在通信用户数量较多时，进行逐用户的doa估计，可以更好地适应天线阵列的过载状态。
[0081]
将接收信号输入码匹配滤波器，得到k个如公式(2)所示的解扩信号，之后对每个用户依次进行doa估计。常用的高分辨率doa估计算法通常需要谱峰搜索、特征分解等操作，在用户数量较多时，系统开销较为庞大。因此本技术利用解扩后信号的等效单源特性，设计简化二维esprit算法，本技术的简化二维esprit算法将在后文详细说明。
[0082]
图2是本技术一实施例的二维doa估计方法的天线阵列示意图。
[0083]
在一些实施例中，参考图2所示，天线阵列200为n
×
n维的阵列，第一子阵a0、第二子阵a1和第三子阵a2均为(n-1)
×
(n-1)维的方阵，将天线阵列200划分为第一子阵a0、第二子阵a1和第三子阵a2的步骤包括：将天线阵列200剔除第n行阵元2001与第n列阵元2001获得第一子阵a0；将天线阵列200剔除第n行阵元2001与第1列阵元2001获得第二子阵a1；将天线阵列200剔除第1行阵元2001与第n列阵元2001获得第三子阵a2。
[0084]
设第k个用户信号的来波方向为由于解扩后信号可近似看作单源信号，为便于表述，将来波方向简记为类似地，分别将前文所述的等效单源信号sk(i)，信号幅值gσk和第k个用户信号对应的导向矢量ak也另记为单源信号s(i)，信号幅值σs和导向矢量a。参考图2所示，考虑阵元2001数为n
×
n的均匀平面天线阵列200，将该天线阵列200划分为三个子阵，即第一子阵a0、第二子阵a1和第三子阵a2。将x轴方向n个阵元2001对应的导向矢量用公式(3)表示：
[0085][0086]
将y轴方向对应的导向矢量用公式(4)表示：
[0087][0088]
在公式(3)和公式(4)中，j为虚数，d表示阵元间距，λ表示接收信号的波长。
[0089]
该平面阵列的导向矢量a可表示为a
x
与ay的克罗内克积(kronecker)积，用公式(5)表示：
[0090][0091]
对于第一子阵a0用公式(6)表示：
[0092][0093]
第一子阵a0对应的导向矢量可用公式(7)表示：
[0094][0095]
在公式(6)和公式(7)中，(0:n-2)中的0表示第0个元素，n-2表示第n-2个元素。
[0096]
类似地，对于第二子阵a1、第三子阵a2对应的导向矢量，分别可以用公式(8)、公式(9)表示：
[0097][0098][0099]
结合公式(7)、(8)、(9)，可以将第二子阵a1对应的导向矢量用公式(10)表示，将第三子阵a2对应的导向矢量用公式(11)表示：
[0100][0101][0102]
在公式(10)中，φ1为第二子阵a1的旋转矩阵，在公式(11)中，φ2为第三子阵a2的
旋转矩阵。在经典esprit算法中，旋转矩阵是维数等于入射信号数的对角矩阵，由于对角元素的相位包含信号的doa信息，因此可通过对旋转矩阵的相似矩阵进行特征分解来求取信号的doa信息。对于解扩信号，有效信源数为1，因此公式(10)、(11)中旋转矩阵φ1和φ2的维数也均为1，即旋转矩阵φ1和φ2是两个标量，可直接通过估计旋转矩阵的相位获取信号的doa信息。
[0103]
将旋转矩阵φ1和φ2的相位用公式(12)表示：
[0104][0105]
其中，angle(
·
)函数用于求解相位角。
[0106]
接下来具体描述使用点除或元素平均法、相位平均法、最小二乘法中的任意一种方法代替标准esprit算法中的特征分解操作，对旋转矩阵φ1和φ2的相位φ1和φ2进行估计。
[0107]
对旋转矩阵φ1的相位φ1进行估计、旋转矩阵φ2的相位φ2进行估计有如下文所述的推导过程。
[0108]
对于第一子阵a0解扩信号和第二子阵a1解扩信号可用公式(13)表示：
[0109][0110]
其中，x0(i)表示第一子阵a0解扩信号，x1(i)表示第二子阵a1解扩信号，s(i)表示等效单源信号，n0(i)表示第一子阵a0的第i个信号的等效噪声矢量，n1(i)表示第二子阵a1的第i个信号的等效噪声矢量。
[0111]
用下面的公式(14)计算第一子阵a0解扩信号x0(i)的自相关矩阵r
00
以及第一子阵a0解扩信号与第二子阵a1解扩信号的互相关矩阵r
01
：
[0112][0113]
其中，e[
·
]函数表示求均值，h表示矢量的共轭转置，表示第一子阵a0解扩信号的共轭转置，表示第二子阵a1解扩信号的共轭转置，表示等效单源信号s(i)的信号功率，表示第一子阵a0解扩信号的共轭转置，表示等效噪声功率，表示旋转矩阵φ1的共轭转置，i为(n-1)2阶单位矩阵；γ
01
可用下面的公式(15)表示：
[0114]
[0115][0116]
在公式(15)、(16)中，diag(
·
)函数用于构造一个对角矩阵，矩阵里的ο中元素均为0，公式(15)γ
01
的矩阵维度为(n-1)2×
(n-1)2维，公式(16)λ
n-1
的矩阵维度为(n-1)
×
(n-1)维，λ
n-1
是左下副对角线为1，其余元素为0的(n-1)阶矩阵；γ
01
是左下副对角线位置(n-2)个元素1和1个元素0交替出现，其余元素全为0的(n-1)2阶矩阵，即一个周期有n-1个元素，前n-2个是1，最后1个是0。
[0117]
用下面的公式(17)计算第一子阵a0解扩信号x0(i)的自相关矩阵r
00
以及第一子阵a0解扩信号与第三子阵a2解扩信号的互相关矩阵r
02
：
[0118][0119]
其中，e[
·
]函数表示求均值，h表示矢量的共轭转置，表示第一子阵a0解扩信号的共轭转置，表示第三子阵a2解扩信号的共轭转置，表示等效单源信号s(i)的信号功率，表示第一子阵a0解扩信号的共轭转置，表示等效噪声功率，表示旋转矩阵φ2的共轭转置，i为(n-1)2阶单位矩阵；γ
02
可用下面的公式(18)表示：
[0120][0121]
其中，γ
02
是(n-1)2阶矩阵，γ
02
左下角为((n-1)(n-2))阶单位矩阵，其余元素为0。
[0122]
从上述公式(14)中，分别减去后，用公式(19)表示：
[0123][0124]
由于旋转矩阵φ1为标量，根据公式(19)可以得到下面的公式(20)：
[0125][0126]
在公式(20)中，*表示共轭，表示旋转矩阵φ1的共轭，是一维的标量。
[0127]
根据公式(20)对矩阵进行点除，得到公式(21)：
[0128]
[φ1e]
i,j
＝[q
00
]
i,j
/[q
01
]
i,j
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(21)
[0129]
其中，φ1e为估计矩阵，e为(n-1)2阶的全1矩阵，[
·
]
i,j
表示矩阵第i行第j列元素。上述公式(21)表示：估计矩阵φ1e第i行第j列元素，由矩阵q
00
第i行第j列元素除以矩阵q
01
第i行第j列元素得到。可由矩阵φ1e获得对相位φ1的估计。由于计算q
00
和q
01
需要已知而对的估计通常需要特征分解自相关矩阵r
00
，为简化计算，避开特征分解操作，可利用自相关矩阵r
00
和互相关矩阵r
01
自身的结构特点，剔除指定位置元素得到重构矢量，对重构矢量进行点除操作。
[0130]
在实际应用中，自相关矩阵r
00
、互相关矩阵r
01
和互相关矩阵r
02
往往由各子阵解扩信号估计得到。
[0131]
在一些实施例中，采用下面的公式(22)、(23)、(24)分别计算自相关矩阵互相关矩阵和互相关矩阵
[0132][0133][0134][0135]
其中，s表示用户的符号信息数，x0(i)表示第一子阵a0解扩信号，表示第一子阵a0解扩信号的共轭转置，表示第二子阵a1解扩信号的共轭转置，表示第三子阵a2解扩信号的共轭转置。示例性地，是实际使用有限数量的接收数据样本点对自相关矩阵r
00
的估计，可表示为样本自相关矩阵。
[0136]
图3是本技术一实施例的二维doa估计方法的相关矩阵重构示意图。
[0137]
在步骤s120中，剔除自相关矩阵和互相关矩阵的主对角线元素与左下副对角线元素，并将自相关矩阵和互相关矩阵的剩余元素分别重构为列矢量形式，获得重构矢量和重构矢量以及，剔除自相关矩阵和互相关矩阵的主对角线元素与左下角((n-1)(n-2))阶矩阵的主对角线元素，并将自相关矩阵和互相关矩阵的剩余元素分别重构为列矢量形式，获得重构矢量和重构矢量
[0138]
参考图3所示，图3的五角星301表示需要剔除的元素，实心圆302表示需要保留的元素。
[0139]
由于计算q
00
和q
01
需要已知而对的估计通常需要特征分解自相关矩阵r
00
。为简化计算，避开特征分解操作，可利用自相关矩阵r
00
和互相关矩阵r
01
自身的结构特点，将自相关矩阵r
00
和互相关矩阵r
01
主对角线与左下副对角线的元素剔除，并将两矩阵的剩余元素分别重构为列矢量形式，得到重构矢量和重构矢量和为((n-1)
2-1)2×
1维矢量。类似的，将r
00
和r
02
的主对角线元素，以及左下角((n-1)(n-2))阶矩阵的主对角线元素剔除，将两个相关矩阵的剩余元素分别重构为列矢量形式，得到重构矢量和
[0140]
在步骤s130中，根据重构矢量和重构矢量计算估计矩阵的步骤包括：采用下面的公式(25)计算估计矩阵
[0141][0142]
其中，φ1为第二子阵a1的旋转矩阵，为((n-1)
2-1)2×
1维全1矢量，[
·
]i表示矩阵第i行元素。
[0143]
在步骤s140中，根据重构矢量和重构矢量计算估计矩阵的步骤包括：
采用下面的公式(26)计算估计矩阵
[0144][0145]
其中，φ2为第三子阵a2的旋转矩阵，为全1矢量，[
·
]i表示矩阵第i行元素。
[0146]
在步骤s150中，根据估计矩阵和估计矩阵使用元素平均法、相位平均法、最小二乘法中的任意一种方法计算相位估计和相位估计相位估计和相位估计包括解扩信号的doa信息，doa信息包括俯仰角和/或方位角。
[0147]
在一些实施例中，根据估计矩阵使用元素平均(element average，ea)法计算相位估计的步骤包括：采用下面的公式(27)、(28)进行计算：
[0148][0149][0150]
其中，表示使用元素平均法得到的第二子阵a1的旋转矩阵估计值，n表示天线阵列的维度，表示使用元素平均法得到的相位估计，angle(
·
)函数用于求解相位角。在公式(27)对中所有元素取均值，在公式(28)对取相位。
[0151]
类似的，根据估计矩阵使用元素平均法计算相位估计的步骤可以采用上述公式(27)、(28)进行计算，在此不再赘述。
[0152]
在一些实施例中，根据估计矩阵使用相位平均(phase average，pa)法计算相位估计的步骤包括：采用下面的公式(29)进行计算：
[0153][0154]
其中，表示使用相位平均法得到的相位估计，n表示天线阵列的维度，angle(
·
)函数用于求解相位角。
[0155]
类似的，根据估计矩阵使用相位平均法计算相位估计的步骤可以采用上述公式(29)进行计算，在此不再赘述。
[0156]
在一些实施例中，根据估计矩阵使用最小二乘(least squares，ls)法计算相位估计的步骤包括：采用下面的公式(30)、(31)、(32)进行计算：
[0157][0158][0159]
[0160]
其中，表示第二子阵a1的旋转矩阵的共轭，表示重构矢量的伪逆，表示重构矢量的共轭转置，表示使用最小二乘法得到的相位估计，angle(
·
)函数用于求解相位角。根据前述公式(20)可得到公式(30)，根据最小二乘法的思想对公式(30)求伪逆可得到公式(31)，由于只需要旋转矩阵φ1的相位信息，因此只需计算公式(31)中的部分即可，由此得到公式(32)。
[0161]
类似的，根据估计矩阵使用最小二乘法计算相位估计的步骤可以采用上述公式(30)、(31)、(32)进行计算，在此不再赘述。
[0162]
本技术将天线阵列划分为三个子阵，根据子阵间的相关关系，使用点除或元素平均法、相位平均法、最小二乘法中的任意一种方法代替标准esprit算法中系统计算开销较大的矩阵特征分解操作，在保证测向精度的同时，计算复杂度大大降低。
[0163]
接下来由相位估计和相位估计对用户来波方向进行估计。将方向余弦用公式(33)表示：
[0164][0165]
其中，u
x
表示天线阵列中x轴方向的方向余弦，uy表示天线阵列中y轴方向的方向余弦。
[0166]
在步骤s160中，根据相位估计和相位估计分别计算天线阵列中x轴方向的方向余弦估计和y轴方向的方向余弦估计的步骤包括：采用下面的公式(34)计算方向余弦估计和方向余弦估计
[0167][0168]
其中，λ表示接收信号的波长，d表示天线阵列的阵元间距。根据前述公式(12)和(33)可得公式(34)。
[0169]
在步骤s170中，根据方向余弦估计和方向余弦估计计算俯仰角估计和/或方位角估计的步骤包括：采用下面的公式计算俯仰角估计和/或方位角估计
[0170][0171][0172]
在实际应用中可根据需要计算俯仰角估计和方位角估计还可以只计算俯仰角估计或方位角估计本技术不做限制。
[0173]
接下来，进行计算复杂度分析。
[0174]
为方便描述，将各子阵阵元数简记为m
sub
＝(n-1)
×
(n-1)，使用元素平均法、相位
平均法、最小二乘法估计相位，将本技术的简化esprit算法分别记为s-esprit-ea、s-esprit-pa、以及s-esprit-ls算法。对于s-esprit-ea算法，先由和使用点除方法得到与再对各元素求平均，对所得均值取相位，其计算复杂度约为s-esprit-pa在得到与后，先对各元素取相位，再对相位值进行平均，计算复杂度增加为而s-esprit-ls算法直接使用重构矢量，在矢量乘法操作后取相位，计算复杂度约为根据参考文献，标准esprit算法的计算复杂度约为综上所示，本技术的二维doa估计方法相比标准esprit算法具有更低的计算复杂度，且相差约m
sub
的一个数量级。
[0175]
接下来，根据本技术二维doa估计方法的仿真结果进行说明。
[0176]
在仿真中，k个等功率用户信号入射到阵列，用户信息符号数为s，使用扩频增益g＝64的gold序列进行扩频。不失一般性，仿真结果均以第1个用户的doa估计为例，该用户信号入射方向的俯仰角和方位角分别为θ1＝30
°
,其余(k-1)个用户信号入射方向的俯仰角θk,k＝2,3,...,k服从区间[0
°
,90
°
]内的均匀分布，方位角服从区间[0
°
,180
°
]内的均匀分布。定义doa估计误差δα作为测向精度的评价指标，用公式(37)表示：
[0177][0178]
其中，nm为蒙特卡洛次数，后续仿真中取nm＝200，α(α代表θ或)为doa真值，为第i次的doa估计值。
[0179]
考虑阵元分布为n
×
n的矩形栅格阵列，阵元间距d＝λ/2。经码匹配滤波后，使用本技术二维doa估计方法，即简化esprit算法(s-esprit-ea，s-esprit-pa和s-esprit-ls)以及标准esprit算法进行二维doa估计，比较各算法在不同仿真条件下的测向性能。
[0180]
1)不同信噪比下简化esprit算法测向性能。
[0181]
图4是本技术一实施例的二维doa估计方法仿真结果的俯仰角估计误差与输入snr关系曲线示意图；图5是本技术一实施例的二维doa估计方法仿真结果的方位角估计误差与输入snr关系曲线示意图。
[0182]
在此，验证本技术的简化esprit算法在不同输入信噪比(signal-to-noise ratio，snr)下的测向性能。参考图4和图5所示，假设用户数k＝4，用户信息符号数s＝200，矩形栅格阵列行/列数n＝6。输入snr变化范围为[-25db,15db]。图4和图5分别为第1个用户信号入射方向的俯仰角θ1和方位角估计误差与输入snr的关系曲线。可以看出，本技术的简化esprit算法在降低了计算复杂度的同时，算法性能与标准esprit算法相差无几。在输入snr＞-10db时，各算法的测向误差曲线几乎完全重合。信噪比较低时，s-esprit-ea算法测向误差增大较为明显，而s-esprit-ls算法与标准esprit算法性能最为接近。
[0183]
2)不同用户信息符号数下简化esprit算法测向性能。
[0184]
图6是本技术一实施例的二维doa估计方法仿真结果的俯仰角估计误差与用户信息符号数关系曲线示意图；图7是本技术一实施例的二维doa估计方法仿真结果的方位角估
计误差与用户信息符号数关系曲线示意图。
[0185]
在此验证本技术简化esprit算法在不同用户信息符号数下的测向性能。参考图6和图7所示，假设用户数k＝4，矩形栅格阵列行/列数n＝6，输入snr固定为-10db,用户信息符号数s的变化范围为[50,300]。θ1和的估计误差与符号数s的关系曲线分别如图6、图7所示。可以看出，随着符号数的增加，本技术简化esprit算法与标准esprit算法的俯仰角与方位角测向误差均有所下降。由于上述算法都需要使用信号解扩后的协方差矩阵进行测向，所以对用户信息符号数变化的敏感程度也几乎一致。
[0186]
3)不同用户数下简化esprit算法测向性能。
[0187]
图8是本技术一实施例的二维doa估计方法仿真结果的俯仰角估计误差与用户数关系曲线示意图；图9是本技术一实施例的二维doa估计方法仿真结果的方位角估计误差与用户数关系曲线示意图。
[0188]
在此验证本技术简化esprit算法在不同用户数下的测向性能。参考图8和图9所示，假设矩形栅格阵列行/列数n＝4，用户信息符号数s＝200，输入snr固定为-10db,用户数分别取k＝1,5,15,45,60,64,65,70。θ1和的估计误差与用户数变化关系如图8、图9所示。由图可知，当用户数小于60时，本技术简化esprit算法与标准esprit算法测向性能受用户数变化的影响不明显，且俯仰角与方位角测向误差均在0.2
°
以内。当用户数大于64时，各算法的测向误差显著增大。对于阵元分布为4
×
4的均匀平面阵列，用户数在60左右时已明显超出阵列的自由度，由于各用户在码域接近于正交，对接收信号经码匹配滤波后进行逐用户的doa估计，可以更好地适应天线阵列的过载状态。
[0189]
本技术还包括一种二维doa估计系统，包括天线阵列、存储器和处理器。其中，该天线阵列用于接收信号；存储器用于存储可由处理器执行的指令；处理器用于执行该指令以实现前文所述的二维doa估计方法。
[0190]
图10是本技术一实施例的二维doa估计系统的系统框图。参考图10所示，该二维doa估计系统1000可包括内部通信总线1001、处理器1002、只读存储器(rom)1003、随机存取存储器(ram)1004以及通信端口1005。当应用在个人计算机上时，该二维doa估计系统1000还可以包括硬盘1006。内部通信总线1001可以实现该二维doa估计系统1000组件间的数据通信。处理器1002可以进行判断和发出提示。在一些实施例中，处理器1002可以由一个或多个处理器组成。通信端口1005可以实现该二维doa估计系统1000与外部的数据通信。在一些实施例中，该二维doa估计系统1000可以通过通信端口1005从网络发送和接受信息及数据。该二维doa估计系统1000还可以包括不同形式的程序储存单元以及数据储存单元，例如硬盘1006，只读存储器(rom)1003和随机存取存储器(ram)1004，能够存储计算机处理和/或通信使用的各种数据文件，以及处理器1002所执行的可能的程序指令。处理器执行这些指令以实现方法的主要部分。处理器处理的结果通过通信端口传给用户设备，在用户界面上显示。
[0191]
上述的二维doa估计方法可以实施为计算机程序，保存在硬盘1006中，并可加载到处理器1002中执行，以实施本技术的二维doa估计方法。
[0192]
本技术还包括一种存储有计算机程序代码的计算机可读介质，该计算机程序代码在由处理器执行时实现前文所述的二维doa估计方法。
[0193]
二维doa估计方法实施为计算机程序时，也可以存储在计算机可读存储介质中作
为制品。例如，计算机可读存储介质可以包括但不限于磁存储设备(例如，硬盘、软盘、磁条)、光盘(例如，压缩盘(cd)、数字多功能盘(dvd))、智能卡和闪存设备(例如，电可擦除可编程只读存储器(eprom)、卡、棒、键驱动)。此外，本文描述的各种存储介质能代表用于存储信息的一个或多个设备和/或其它机器可读介质。术语“机器可读介质”可以包括但不限于能存储、包含和/或承载代码和/或指令和/或数据的无线信道和各种其它介质(和/或存储介质)。
[0194]
应该理解，上文所描述的实施例仅是示意。本文描述的实施例可在硬件、软件、固件、中间件、微码或者其任意组合中实现。对于硬件实现，处理器可以在一个或者多个特定用途集成电路(asic)、数字信号处理器(dsp)、数字信号处理设备(dspd)、可编程逻辑器件(pld)、现场可编程门阵列(fpga)、处理器、控制器、微控制器、微处理器和/或设计为执行本文所述功能的其它电子单元或者其结合内实现。
[0195]
本技术的一些方面可以完全由硬件执行、可以完全由软件(包括固件、常驻软件、微码等)执行、也可以由硬件和软件组合执行。以上硬件或软件均可被称为“数据块”、“模块”、“引擎”、“单元”、“组件”或“系统”。处理器可以是一个或多个专用集成电路(asic)、数字信号处理器(dsp)、数字信号处理器件(dapd)、可编程逻辑器件(pld)、现场可编程门阵列(fpga)、处理器、控制器、微控制器、微处理器或者其组合。此外，本技术的各方面可能表现为位于一个或多个计算机可读介质中的计算机产品，该产品包括计算机可读程序编码。例如，计算机可读介质可包括，但不限于，磁性存储设备(例如，硬盘、软盘、磁带
……
)、光盘(例如，压缩盘cd、数字多功能盘dvd
……
)、智能卡以及闪存设备(例如，卡、棒、键驱动器
……
)。
[0196]
计算机可读介质可能包含一个内含有计算机程序编码的传播数据信号，例如在基带上或作为载波的一部分。该传播信号可能有多种表现形式，包括电磁形式、光形式等等、或合适的组合形式。计算机可读介质可以是除计算机可读存储介质之外的任何计算机可读介质，该介质可以通过连接至一个指令执行系统、装置或设备以实现通讯、传播或传输供使用的程序。位于计算机可读介质上的程序编码可以通过任何合适的介质进行传播，包括无线电、电缆、光纤电缆、射频信号、或类似介质、或任何上述介质的组合。
[0197]
上文已对基本概念做了描述，显然，对于本领域技术人员来说，上述申请披露仅仅作为示例，而并不构成对本技术的限定。虽然此处并没有明确说明，本领域技术人员可能会对本技术进行各种修改、改进和修正。该类修改、改进和修正在本技术中被建议，所以该类修改、改进、修正仍属于本技术示范实施例的精神和范围。
[0198]
同时，本技术使用了特定词语来描述本技术的实施例。如“一个实施例”、“一实施例”、和/或“一些实施例”意指与本技术至少一个实施例相关的某一特征、结构或特点。因此，应强调并注意的是，本说明书中在不同位置两次或多次提及的“一实施例”或“一个实施例”或“一替代性实施例”并不一定是指同一实施例。此外，本技术的一个或多个实施例中的某些特征、结构或特点可以进行适当的组合。
[0199]
一些实施例中使用了描述成分、属性数量的数字，应当理解的是，此类用于实施例描述的数字，在一些示例中使用了修饰词“大约”、“近似”或“大体上”来修饰。除非另外说明，“大约”、“近似”或“大体上”表明所述数字允许有
±
20％的变化。相应地，在一些实施例中，说明书和权利要求中使用的数值参数均为近似值，该近似值根据个别实施例所需特点
可以发生改变。在一些实施例中，数值参数应考虑规定的有效数位并采用一般位数保留的方法。尽管本技术一些实施例中用于确认其范围广度的数值域和参数为近似值，在具体实施例中，此类数值的设定在可行范围内尽可能精确。