CN115308680A - 二维doa估计方法、系统及计算机可读介质 - Google Patents

二维doa估计方法、系统及计算机可读介质 Download PDF

Info

Publication number
CN115308680A
CN115308680A CN202210932982.2A CN202210932982A CN115308680A CN 115308680 A CN115308680 A CN 115308680A CN 202210932982 A CN202210932982 A CN 202210932982A CN 115308680 A CN115308680 A CN 115308680A
Authority
CN
China
Prior art keywords
array
matrix
estimation
sub
phase
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210932982.2A
Other languages
English (en)
Inventor
杨佳硕
梁广
赵璐璐
马菁涛
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Engineering Center for Microsatellites
Innovation Academy for Microsatellites of CAS
Original Assignee
Shanghai Engineering Center for Microsatellites
Innovation Academy for Microsatellites of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Engineering Center for Microsatellites, Innovation Academy for Microsatellites of CAS filed Critical Shanghai Engineering Center for Microsatellites
Priority to CN202210932982.2A priority Critical patent/CN115308680A/zh
Publication of CN115308680A publication Critical patent/CN115308680A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/02—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using radio waves
    • G01S3/14—Systems for determining direction or deviation from predetermined direction

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Radio Transmission System (AREA)

Abstract

本发明涉及一种二维DOA估计方法、系统及计算机可读介质。该估计方法包括：根据第一子阵A0、第二子阵A1和第三子阵A2的解扩信号，计算自相关矩阵
Figure DDA0003782236260000011
互相关矩阵
Figure DDA0003782236260000012
和
Figure DDA0003782236260000013
剔除
Figure DDA0003782236260000014
和
Figure DDA0003782236260000015
中指定元素，获得重构矢量
Figure DDA0003782236260000016
和
Figure DDA00037822362600000125
根据
Figure DDA0003782236260000017
和
Figure DDA0003782236260000018
计算估计矩阵
Figure DDA0003782236260000019
根据
Figure DDA00037822362600000110
和
Figure DDA00037822362600000111
计算估计矩阵
Figure DDA00037822362600000112
根据
Figure DDA00037822362600000113
和
Figure DDA00037822362600000114
使用元素平均法、相位平均法、最小二乘法中的任意一种方法计算相位估计
Figure DDA00037822362600000115
和
Figure DDA00037822362600000116
根据
Figure DDA00037822362600000117
和
Figure DDA00037822362600000118
分别计算方向余弦估计
Figure DDA00037822362600000119
和
Figure DDA00037822362600000120
根据
Figure DDA00037822362600000121
和
Figure DDA00037822362600000122
计算俯仰角估计
Figure DDA00037822362600000123
和/或方位角估计
Figure DDA00037822362600000124
本申请利用解扩信号的等效单源特性，使用点除或最小二乘方法代替标准ESPRIT算法中的矩阵特征分解操作，在保证测向精度的同时，计算复杂度大大降低。

Description

二维DOA估计方法、系统及计算机可读介质
技术领域
本发明主要涉及天线阵列信号处理技术领域，具体地涉及一种二维DOA估计方法、系统及计算机可读介质。
背景技术
码分多址(Code Division Multiple Access，CDMA)通信体制具有抗干扰、抗截获能力强，通信隐蔽，易于实现多址通信、语音激活、软容量和软切换等特点，在低轨卫星通信系统得到广泛应用。在CDMA通信系统中，天线阵列接收的信号数往往远大于阵元数目，导致传统的波达方向(Direction Of Arrival，DOA)估计算法难以直接应用。即使在通信用户较少的情况下，由于扩频信号的信噪比很低，直接对阵元接收信号进行DOA估计，效果也不甚理想。
针对上述问题，在通信用户扩频序列已知的情况下，通常使用码匹配滤波器对天线阵列的接收信号进行解扩预处理。由于扩频序列良好的自相关与互相关特性，对某一通信用户信号解扩后，其他通信用户信号可近似为噪声，解扩信号可看作该通信用户信号与噪声的合成信号，且信噪比得到提升。对于各通信用户解扩后的等效单源信号，可使用常规DOA估计方法逐用户进行角度估计。常用的高分辨率DOA估计算法如多信号分类(MultipleSignal Classification，MUSIC)算法、基于旋转不变技术的信号参数估计(EstimatingSignal Parameters via Rotational Invariance Techniques，ESPRIT)算法等。
在对解扩信号逐用户进行DOA估计时，由于常规DOA估计算法需要谱峰搜索、特征分解等操作，在通信用户数量较多时，存在系统计算开销较为庞大的问题。寻找能够在CDMA通信系统中具有稳定测向性能、低计算复杂度的DOA估计算法，在理论与工程应用中都具有重要的研究价值。
发明内容
本申请所要解决的技术问题是提供一种二维DOA估计方法、系统及计算机可读介质，本申请的估计方法测向稳定，在保证测向精度的同时，计算复杂度大大降低。
本申请为解决上述技术问题而采用的技术方案是一种二维DOA估计方法，该二维DOA估计方法应用于CDMA通信体制下，包括：将天线阵列的接收信号进行解扩预处理，获得解扩信号，将天线阵列划分为第一子阵A0、第二子阵A1和第三子阵A2，计算第一子阵A0解扩信号的自相关矩阵
Figure BDA0003782236240000021
第一子阵A0解扩信号与第二子阵A1解扩信号的互相关矩阵
Figure BDA0003782236240000022
第一子阵A0解扩信号与第三子阵A2解扩信号的互相关矩阵
Figure BDA0003782236240000023
剔除自相关矩阵
Figure BDA0003782236240000024
和互相关矩阵
Figure BDA0003782236240000025
的主对角线元素与左下副对角线元素，并将自相关矩阵
Figure BDA0003782236240000026
和互相关矩阵
Figure BDA0003782236240000027
的剩余元素分别重构为列矢量形式，获得重构矢量
Figure BDA0003782236240000028
和重构矢量
Figure BDA0003782236240000029
以及，剔除自相关矩阵
Figure BDA00037822362400000210
和互相关矩阵
Figure BDA00037822362400000211
的主对角线元素与左下角((N-1)(N-2))阶矩阵的主对角线元素，并将自相关矩阵
Figure BDA00037822362400000212
和互相关矩阵
Figure BDA00037822362400000213
的剩余元素分别重构为列矢量形式，获得重构矢量
Figure BDA00037822362400000214
和重构矢量
Figure BDA00037822362400000215
根据重构矢量
Figure BDA00037822362400000216
和重构矢量
Figure BDA00037822362400000217
计算估计矩阵
Figure BDA00037822362400000218
根据重构矢量
Figure BDA00037822362400000219
和重构矢量
Figure BDA00037822362400000220
计算估计矩阵
Figure BDA00037822362400000221
根据估计矩阵
Figure BDA00037822362400000222
和估计矩阵
Figure BDA00037822362400000223
使用元素平均法、相位平均法、最小二乘法中的任意一种方法计算相位估计
Figure BDA00037822362400000224
和相位估计
Figure BDA00037822362400000225
相位估计
Figure BDA00037822362400000226
和相位估计
Figure BDA00037822362400000227
包括解扩信号的DOA信息，DOA信息包括俯仰角和/或方位角；根据相位估计
Figure BDA00037822362400000228
和相位估计
Figure BDA00037822362400000229
分别计算天线阵列中x轴方向的方向余弦估计
Figure BDA00037822362400000230
和y轴方向的方向余弦估计
Figure BDA00037822362400000231
以及根据方向余弦估计
Figure BDA00037822362400000232
和方向余弦估计
Figure BDA00037822362400000233
计算俯仰角估计
Figure BDA00037822362400000234
和/或方位角估计
Figure BDA00037822362400000235
在本申请的一实施例中，将天线阵列的接收信号进行解扩预处理，获得解扩信号的步骤包括：采用下面的公式计算解扩信号：
Figure BDA00037822362400000236
Figure BDA00037822362400000237
其中，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个符号周期的延拓；
yk(i)表示解扩信号，解扩信号具有等效单源特性，k表示第k个用户，k'表示第k个用户以外的其他用户，
Figure BDA0003782236240000031
表示第k个用户的扩频序列与其他用户扩频序列的互相关函数；等效单源信号sk(i)＝Gσkdk(i)，等效噪声
Figure BDA0003782236240000032
在本申请的一实施例中，天线阵列为N×N维的阵列，第一子阵A0、第二子阵A1和第三子阵A2均为(N-1)×(N-1)维的方阵，将天线阵列划分为第一子阵A0、第二子阵A1和第三子阵A2的步骤包括：将天线阵列剔除第N行阵元与第N列阵元获得第一子阵A0；将天线阵列剔除第N行阵元与第1列阵元获得第二子阵A1；将天线阵列剔除第1行阵元与第N列阵元获得第三子阵A2。
在本申请的一实施例中，采用下面的公式分别计算自相关矩阵
Figure BDA0003782236240000033
互相关矩阵
Figure BDA0003782236240000034
和互相关矩阵
Figure BDA0003782236240000035
Figure BDA0003782236240000036
Figure BDA0003782236240000037
Figure BDA0003782236240000038
其中，S表示用户的符号信息数，x0(i)表示第一子阵A0解扩信号，
Figure BDA0003782236240000039
表示第一子阵A0解扩信号的共轭转置，
Figure BDA00037822362400000310
表示第二子阵A1解扩信号的共轭转置，
Figure BDA00037822362400000311
表示第三子阵A2解扩信号的共轭转置。
在本申请的一实施例中，根据重构矢量
Figure BDA00037822362400000312
和重构矢量
Figure BDA00037822362400000313
计算估计矩阵
Figure BDA00037822362400000314
的步骤包括：采用下面的公式计算估计矩阵
Figure BDA00037822362400000315
Figure BDA00037822362400000316
其中，Φ1为第二子阵A1的旋转矩阵，
Figure BDA00037822362400000317
为全1矢量，[·]i表示矩阵第i行元素。
在本申请的一实施例中，根据重构矢量
Figure BDA0003782236240000041
和重构矢量
Figure BDA0003782236240000042
计算估计矩阵
Figure BDA0003782236240000043
的步骤包括：采用下面的公式计算估计矩阵
Figure BDA0003782236240000044
Figure BDA0003782236240000045
其中，Φ2为第三子阵A2的旋转矩阵，
Figure BDA0003782236240000046
为全1矢量，[·]i表示矩阵第i行元素。
在本申请的一实施例中，根据估计矩阵
Figure BDA0003782236240000047
使用元素平均法计算相位估计
Figure BDA0003782236240000048
的步骤包括：采用下面的公式进行计算：
Figure BDA0003782236240000049
Figure BDA00037822362400000410
其中，
Figure BDA00037822362400000411
表示使用元素平均法得到的第二子阵A1的旋转矩阵估计值，N表示天线阵列的维度，
Figure BDA00037822362400000412
表示使用元素平均法得到的相位估计，angle(·)函数用于求解相位角。
在本申请的一实施例中，根据估计矩阵
Figure BDA00037822362400000413
使用相位平均法计算相位估计
Figure BDA00037822362400000414
的步骤包括：采用下面的公式进行计算：
Figure BDA00037822362400000415
其中，
Figure BDA00037822362400000416
表示使用相位平均法得到的相位估计，N表示天线阵列的维度，angle(·)函数用于求解相位角。
在本申请的一实施例中，根据估计矩阵
Figure BDA00037822362400000417
使用最小二乘法计算相位估计
Figure BDA00037822362400000418
的步骤包括：采用下面的公式进行计算：
Figure BDA00037822362400000419
Figure BDA00037822362400000420
Figure BDA00037822362400000421
其中，
Figure BDA00037822362400000422
表示第二子阵A1的旋转矩阵的共轭，
Figure BDA00037822362400000423
表示重构矢量
Figure BDA00037822362400000424
的伪逆，
Figure BDA00037822362400000425
表示重构矢量
Figure BDA00037822362400000426
的共轭转置，
Figure BDA00037822362400000427
表示使用最小二乘法得到的相位估计，angle(·)函数用于求解相位角。
在本申请的一实施例中，根据相位估计
Figure BDA00037822362400000428
和相位估计
Figure BDA00037822362400000429
分别计算天线阵列中x轴方向的方向余弦估计
Figure BDA00037822362400000430
和y轴方向的方向余弦估计
Figure BDA00037822362400000431
的步骤包括：采用下面的公式计算方向余弦估计
Figure BDA0003782236240000051
和方向余弦估计
Figure BDA0003782236240000052
Figure BDA0003782236240000053
其中，λ表示接收信号的波长，d表示天线阵列的阵元间距。
在本申请的一实施例中，根据方向余弦估计
Figure BDA0003782236240000054
和方向余弦估计
Figure BDA0003782236240000055
计算俯仰角估计
Figure BDA0003782236240000056
和/或方位角估计
Figure BDA0003782236240000057
的步骤包括：采用下面的公式计算俯仰角估计
Figure BDA0003782236240000058
和/或方位角估计
Figure BDA0003782236240000059
Figure BDA00037822362400000510
Figure BDA00037822362400000511
本申请为解决上述技术问题还提出一种二维DOA估计系统，包括：天线阵列，天线阵列用于接收信号；存储器，用于存储可由处理器执行的指令；处理器，用于执行指令以实现如上的二维DOA估计方法。
本申请为解决上述技术问题还提出一种存储有计算机程序代码的计算机可读介质，计算机程序代码在由处理器执行时实现如上的二维DOA估计方法。
本申请的技术方案通过对接收信号进行解扩预处理后，解扩信号具有等效单源特性，在通信用户数量较多时，进行逐用户的DOA估计，可以更好地适应天线阵列的过载状态；将天线阵列划分为三个子阵，根据子阵间的相关关系，使用点除或元素平均法、相位平均法、最小二乘法中的任意一种方法代替标准ESPRIT算法中系统计算开销较大的矩阵特征分解操作，在保证测向精度的同时，计算复杂度大大降低；根据天线阵列中的方向余弦估计计算俯仰角估计与方位角估计，测向稳定。
附图说明
为让本申请的上述目的、特征和优点能更明显易懂，以下结合附图对本申请的具体实施方式作详细说明，其中：
图1是本申请一实施例的二维DOA估计方法的示例性流程图；
图2是本申请一实施例的二维DOA估计方法的天线阵列示意图；
图3是本申请一实施例的二维DOA估计方法的相关矩阵重构示意图；
图4是本申请一实施例的二维DOA估计方法仿真结果的俯仰角估计误差与输入SNR关系曲线示意图；
图5是本申请一实施例的二维DOA估计方法仿真结果的方位角估计误差与输入SNR关系曲线示意图；
图6是本申请一实施例的二维DOA估计方法仿真结果的俯仰角估计误差与用户信息符号数关系曲线示意图；
图7是本申请一实施例的二维DOA估计方法仿真结果的方位角估计误差与用户信息符号数关系曲线示意图；
图8是本申请一实施例的二维DOA估计方法仿真结果的俯仰角估计误差与用户数关系曲线示意图；
图9是本申请一实施例的二维DOA估计方法仿真结果的方位角估计误差与用户数关系曲线示意图；
图10是本申请一实施例的二维DOA估计系统的系统框图。
具体实施方式
为让本申请的上述目的、特征和优点能更明显易懂，以下结合附图对本申请的具体实施方式作详细说明。
在下面的描述中阐述了很多具体细节以便于充分理解本申请，但是本申请还可以采用其它不同于在此描述的其它方式来实施，因此本申请不受下面公开的具体实施例的限制。
如本申请和权利要求书中所示，除非上下文明确提示例外情形，“一”、“一个”、“一种”和/或“该”等词并非特指单数，也可包括复数。一般说来，术语“包括”与“包含”仅提示包括已明确标识的步骤和元素，而这些步骤和元素不构成一个排它性的罗列，方法或者设备也可能包含其他的步骤或元素。
本申请中使用了流程图用来说明根据本申请的实施例的系统所执行的操作。应当理解的是，前面或下面操作不一定按照顺序来精确地执行。相反，可以按照倒序或同时处理各种步骤。同时，或将其他操作添加到这些过程中，或从这些过程移除某一步或数步操作。
本申请提出一种二维DOA估计方法，该方法应用于CDMA通信体制下。CDMA是在无线通讯上使用的技术，CDMA允许所有通信用户同时使用全部频带，且把其他通信用户发出讯号视为杂讯，完全不必考虑讯号碰撞问题。DOA估计是一种波达方向定位技术，DOA估计通过处理接收到的回波信号，获取目标的距离信息和方位信息。示例性地，本申请的二维DOA估计方法可以应用在基于CDMA通信体制下的低轨卫星通信系统，在低轨卫星通信系统中，通信用户例如是卫星系统里的卫星。本申请的二维DOA估计方法也可以推广到其他使用CDMA通信的系统中，本申请对通信用户不作限制。
图1是本申请一实施例的二维DOA估计方法的示例性流程图，参考图1所示，该实施例的二维DOA估计方法包括以下步骤：
步骤S110：将天线阵列的接收信号进行解扩预处理，获得解扩信号，将天线阵列划分为第一子阵A0、第二子阵A1和第三子阵A2，计算第一子阵A0解扩信号的自相关矩阵
Figure BDA0003782236240000071
第一子阵A0解扩信号与第二子阵A1解扩信号的互相关矩阵
Figure BDA0003782236240000072
第一子阵A0解扩信号与第三子阵A2解扩信号的互相关矩阵
Figure BDA0003782236240000073
步骤S120：剔除自相关矩阵
Figure BDA0003782236240000074
和互相关矩阵
Figure BDA0003782236240000075
的主对角线元素与左下副对角线元素，并将自相关矩阵
Figure BDA0003782236240000076
和互相关矩阵
Figure BDA0003782236240000077
的剩余元素分别重构为列矢量形式，获得重构矢量
Figure BDA0003782236240000078
和重构矢量
Figure BDA0003782236240000079
以及，剔除自相关矩阵
Figure BDA00037822362400000710
和互相关矩阵
Figure BDA00037822362400000711
的主对角线元素与左下角((N-1)(N-2))阶矩阵的主对角线元素，并将自相关矩阵
Figure BDA00037822362400000712
和互相关矩阵
Figure BDA00037822362400000713
的剩余元素分别重构为列矢量形式，获得重构矢量
Figure BDA00037822362400000714
和重构矢量
Figure BDA00037822362400000715
步骤S130：根据重构矢量
Figure BDA00037822362400000716
和重构矢量
Figure BDA00037822362400000717
计算估计矩阵
Figure BDA00037822362400000718
步骤S140：根据重构矢量
Figure BDA00037822362400000719
和重构矢量
Figure BDA00037822362400000720
计算估计矩阵
Figure BDA00037822362400000721
步骤S150：根据估计矩阵
Figure BDA00037822362400000722
和估计矩阵
Figure BDA00037822362400000723
使用元素平均法、相位平均法、最小二乘法中的任意一种方法计算相位估计
Figure BDA00037822362400000724
和相位估计
Figure BDA00037822362400000725
相位估计
Figure BDA00037822362400000726
和相位估计
Figure BDA00037822362400000727
包括解扩信号的DOA信息，DOA信息包括俯仰角和/或方位角。
步骤S160：根据相位估计
Figure BDA00037822362400000728
和相位估计
Figure BDA00037822362400000729
分别计算天线阵列中x轴方向的方向余弦估计
Figure BDA00037822362400000730
和y轴方向的方向余弦估计
Figure BDA00037822362400000731
步骤S170：根据方向余弦估计
Figure BDA00037822362400000732
和方向余弦估计
Figure BDA00037822362400000733
计算俯仰角估计
Figure BDA00037822362400000734
和/或方位角估计。
下面详细说明上述的步骤S110至步骤S170：
在步骤S110中，将天线阵列的接收信号进行解扩预处理，获得解扩信号，将天线阵列划分为第一子阵A0、第二子阵A1和第三子阵A2，计算第一子阵A0解扩信号的自相关矩阵
Figure BDA0003782236240000081
第一子阵A0解扩信号与第二子阵A1解扩信号的互相关矩阵
Figure BDA0003782236240000082
第一子阵A0解扩信号与第三子阵A2解扩信号的互相关矩阵
Figure BDA0003782236240000083
考虑同步直接序列码分多址(Direct Sequence-Code Division MultipleAccess，DS-CDMA)系统中K个用户的情况，假设第k个用户发送的符号信息为dk(i)∈{-1,1}(i＝0,1,...,S-1)，信号经扩频操作并以功率
Figure BDA0003782236240000084
发射。并且假设信号在理想信道中传播，忽略路径衰减及多径效应的影响，在接收端使用M元天线阵列进行接收，并对接收信号进行chip级采样。chip级采样即对扩频序列码字符号进行采样，扩频序列的每一位码字称为一个chip。
基于上述假设，接下来具体说明对天线阵列的接收信号进行解扩预处理，获得解扩信号的步骤以及前述三个子阵的划分方法。
在一些实施例中，将天线阵列的接收信号进行解扩预处理，获得解扩信号的步骤包括：采用下面的公式(1)和(2)计算解扩信号：
Figure BDA0003782236240000085
Figure BDA0003782236240000086
其中，在上述公式(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)。
将接收信号x(n)输入到M个相关器组中，对于第k个用户信号，解扩后输出的第i个符号级信号记为
Figure BDA0003782236240000087
指信号维度为M×1维。
上述公式(2)中，yk(i)表示解扩信号，解扩信号具有等效单源特性，k表示第k个用户，k'表示第k个用户以外的其他用户，
Figure BDA0003782236240000091
表示第k个用户的扩频序列与其他用户扩频序列的互相关函数；等效单源信号sk(i)＝Gσkdk(i)，等效噪声
Figure BDA0003782236240000092
可以证明，当G→∞时，γkk'(i)近似于高斯概率分布，第k个用户之外的信号通过码匹配滤波器后，可以近似当作噪声处理。因此，可将yk(i)看作单源信号sk(i)与噪声
Figure BDA0003782236240000093
的组合，后续可使用常规DOA估计算法对第k个用户信号的来波方向进行估计。其中，单源信号指的是使用第k个用户的扩频序列对含有多个用户信号的接收信号进行解扩后，解扩信号可近似看作只含有第k个用户信号，而其他用户信号近似当作噪声处理。
通过对接收信号进行解扩预处理后，解扩信号具有等效单源特性，在通信用户数量较多时，进行逐用户的DOA估计，可以更好地适应天线阵列的过载状态。
将接收信号输入码匹配滤波器，得到K个如公式(2)所示的解扩信号，之后对每个用户依次进行DOA估计。常用的高分辨率DOA估计算法通常需要谱峰搜索、特征分解等操作，在用户数量较多时，系统开销较为庞大。因此本申请利用解扩后信号的等效单源特性，设计简化二维ESPRIT算法，本申请的简化二维ESPRIT算法将在后文详细说明。
图2是本申请一实施例的二维DOA估计方法的天线阵列示意图。
在一些实施例中，参考图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。
设第k个用户信号的来波方向为
Figure BDA0003782236240000094
由于解扩后信号可近似看作单源信号，为便于表述，将来波方向简记为
Figure BDA0003782236240000095
类似地，分别将前文所述的等效单源信号sk(i)，信号幅值Gσk和第k个用户信号对应的导向矢量ak也另记为单源信号s(i)，信号幅值σs和导向矢量A。参考图2所示，考虑阵元2001数为N×N的均匀平面天线阵列200，将该天线阵列200划分为三个子阵，即第一子阵A0、第二子阵A1和第三子阵A2。将x轴方向N个阵元2001对应的导向矢量用公式(3)表示：
Figure BDA0003782236240000101
将y轴方向对应的导向矢量用公式(4)表示：
Figure BDA0003782236240000102
在公式(3)和公式(4)中，j为虚数，d表示阵元间距，λ表示接收信号的波长。
该平面阵列的导向矢量A可表示为Ax与Ay的克罗内克积(Kronecker)积，用公式(5)表示：
Figure BDA0003782236240000103
对于第一子阵A0用公式(6)表示：
Figure BDA0003782236240000104
第一子阵A0对应的导向矢量可用公式(7)表示：
Figure BDA0003782236240000105
在公式(6)和公式(7)中，(0:N-2)中的0表示第0个元素，N-2表示第N-2个元素。
类似地，对于第二子阵A1、第三子阵A2对应的导向矢量，分别可以用公式(8)、公式(9)表示：
Figure BDA0003782236240000106
Figure BDA0003782236240000107
结合公式(7)、(8)、(9)，可以将第二子阵A1对应的导向矢量用公式(10)表示，将第三子阵A2对应的导向矢量用公式(11)表示：
Figure BDA0003782236240000108
Figure BDA0003782236240000109
在公式(10)中，Φ1为第二子阵A1的旋转矩阵，在公式(11)中，Φ2为第三子阵A2的旋转矩阵。在经典ESPRIT算法中，旋转矩阵是维数等于入射信号数的对角矩阵，由于对角元素的相位包含信号的DOA信息，因此可通过对旋转矩阵的相似矩阵进行特征分解来求取信号的DOA信息。对于解扩信号，有效信源数为1，因此公式(10)、(11)中旋转矩阵Φ1和Φ2的维数也均为1，即旋转矩阵Φ1和Φ2是两个标量，可直接通过估计旋转矩阵的相位获取信号的DOA信息。
将旋转矩阵Φ1和Φ2的相位用公式(12)表示：
Figure BDA0003782236240000111
其中，angle(·)函数用于求解相位角。
接下来具体描述使用点除或元素平均法、相位平均法、最小二乘法中的任意一种方法代替标准ESPRIT算法中的特征分解操作，对旋转矩阵Φ1和Φ2的相位φ1和φ2进行估计。
对旋转矩阵Φ1的相位φ1进行估计、旋转矩阵Φ2的相位φ2进行估计有如下文所述的推导过程。
对于第一子阵A0解扩信号和第二子阵A1解扩信号可用公式(13)表示：
Figure BDA0003782236240000112
其中，x0(i)表示第一子阵A0解扩信号，x1(i)表示第二子阵A1解扩信号，s(i)表示等效单源信号，n0(i)表示第一子阵A0的第i个信号的等效噪声矢量，n1(i)表示第二子阵A1的第i个信号的等效噪声矢量。
用下面的公式(14)计算第一子阵A0解扩信号x0(i)的自相关矩阵R00以及第一子阵A0解扩信号与第二子阵A1解扩信号的互相关矩阵R01：
Figure BDA0003782236240000113
其中，E[·]函数表示求均值，H表示矢量的共轭转置，
Figure BDA0003782236240000114
表示第一子阵A0解扩信号的共轭转置，
Figure BDA0003782236240000115
表示第二子阵A1解扩信号的共轭转置，
Figure BDA0003782236240000116
表示等效单源信号s(i)的信号功率，
Figure BDA0003782236240000117
表示第一子阵A0解扩信号的共轭转置，
Figure BDA0003782236240000118
表示等效噪声功率，
Figure BDA0003782236240000121
表示旋转矩阵Φ1的共轭转置，I为(N-1)2阶单位矩阵；Γ01可用下面的公式(15)表示：
Figure BDA0003782236240000122
Figure BDA0003782236240000123
在公式(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。
用下面的公式(17)计算第一子阵A0解扩信号x0(i)的自相关矩阵R00以及第一子阵A0解扩信号与第三子阵A2解扩信号的互相关矩阵R02：
Figure BDA0003782236240000124
其中，E[·]函数表示求均值，H表示矢量的共轭转置，
Figure BDA0003782236240000125
表示第一子阵A0解扩信号的共轭转置，
Figure BDA0003782236240000126
表示第三子阵A2解扩信号的共轭转置，
Figure BDA0003782236240000127
表示等效单源信号s(i)的信号功率，
Figure BDA0003782236240000128
表示第一子阵A0解扩信号的共轭转置，
Figure BDA0003782236240000129
表示等效噪声功率，
Figure BDA00037822362400001210
表示旋转矩阵Φ2的共轭转置，I为(N-1)2阶单位矩阵；Γ02可用下面的公式(18)表示：
Figure BDA00037822362400001211
其中，Γ02是(N-1)2阶矩阵，Γ02左下角为((N-1)(N-2))阶单位矩阵，其余元素为0。
从上述公式(14)中，分别减去
Figure BDA00037822362400001212
后，用公式(19)表示：
Figure BDA00037822362400001213
由于旋转矩阵Φ1为标量，根据公式(19)可以得到下面的公式(20)：
Figure BDA0003782236240000131
在公式(20)中，*表示共轭，
Figure BDA0003782236240000132
表示旋转矩阵Φ1的共轭，
Figure BDA0003782236240000133
是一维的标量。
根据公式(20)对矩阵进行点除，得到公式(21)：
[Φ1E]i,j＝[Q00]i,j/[Q01]i,j (21)
其中，Φ1E为估计矩阵，E为(N-1)2阶的全1矩阵，[·]i,j表示矩阵第i行第j列元素。上述公式(21)表示：估计矩阵Φ1E第i行第j列元素，由矩阵Q00第i行第j列元素除以矩阵Q01第i行第j列元素得到。可由矩阵Φ1E获得对相位φ1的估计。由于计算Q00和Q01需要已知
Figure BDA0003782236240000134
而对
Figure BDA0003782236240000135
的估计通常需要特征分解自相关矩阵R00，为简化计算，避开特征分解操作，可利用自相关矩阵R00和互相关矩阵R01自身的结构特点，剔除指定位置元素得到重构矢量，对重构矢量进行点除操作。
在实际应用中，自相关矩阵R00、互相关矩阵R01和互相关矩阵R02往往由各子阵解扩信号估计得到。
在一些实施例中，采用下面的公式(22)、(23)、(24)分别计算自相关矩阵
Figure BDA0003782236240000136
互相关矩阵
Figure BDA0003782236240000137
和互相关矩阵
Figure BDA0003782236240000138
Figure BDA0003782236240000139
Figure BDA00037822362400001310
Figure BDA00037822362400001311
其中，S表示用户的符号信息数，x0(i)表示第一子阵A0解扩信号，
Figure BDA00037822362400001312
表示第一子阵A0解扩信号的共轭转置，
Figure BDA00037822362400001313
表示第二子阵A1解扩信号的共轭转置，
Figure BDA00037822362400001314
表示第三子阵A2解扩信号的共轭转置。示例性地，
Figure BDA00037822362400001315
是实际使用有限数量的接收数据样本点对自相关矩阵R00的估计，可表示为样本自相关矩阵。
图3是本申请一实施例的二维DOA估计方法的相关矩阵重构示意图。
在步骤S120中，剔除自相关矩阵
Figure BDA00037822362400001316
和互相关矩阵
Figure BDA00037822362400001317
的主对角线元素与左下副对角线元素，并将自相关矩阵
Figure BDA00037822362400001318
和互相关矩阵
Figure BDA00037822362400001319
的剩余元素分别重构为列矢量形式，获得重构矢量
Figure BDA00037822362400001320
和重构矢量
Figure BDA00037822362400001321
以及，剔除自相关矩阵
Figure BDA00037822362400001322
和互相关矩阵
Figure BDA0003782236240000141
的主对角线元素与左下角((N-1)(N-2))阶矩阵的主对角线元素，并将自相关矩阵
Figure BDA0003782236240000142
和互相关矩阵
Figure BDA0003782236240000143
的剩余元素分别重构为列矢量形式，获得重构矢量
Figure BDA0003782236240000144
和重构矢量
Figure BDA0003782236240000145
参考图3所示，图3的五角星301表示需要剔除的元素，实心圆302表示需要保留的元素。
由于计算Q00和Q01需要已知
Figure BDA0003782236240000146
而对
Figure BDA0003782236240000147
的估计通常需要特征分解自相关矩阵R00。为简化计算，避开特征分解操作，可利用自相关矩阵R00和互相关矩阵R01自身的结构特点，将自相关矩阵R00和互相关矩阵R01主对角线与左下副对角线的元素剔除，并将两矩阵的剩余元素分别重构为列矢量形式，得到重构矢量
Figure BDA0003782236240000148
和
Figure BDA0003782236240000149
重构矢量
Figure BDA00037822362400001410
和
Figure BDA00037822362400001411
为((N-1)2-1)2×1维矢量。类似的，将R00和R02的主对角线元素，以及左下角((N-1)(N-2))阶矩阵的主对角线元素剔除，将两个相关矩阵的剩余元素分别重构为列矢量形式，得到重构矢量
Figure BDA00037822362400001412
和
Figure BDA00037822362400001413
在步骤S130中，根据重构矢量
Figure BDA00037822362400001414
和重构矢量
Figure BDA00037822362400001415
计算估计矩阵
Figure BDA00037822362400001416
的步骤包括：采用下面的公式(25)计算估计矩阵
Figure BDA00037822362400001417
Figure BDA00037822362400001418
其中，Φ1为第二子阵A1的旋转矩阵，
Figure BDA00037822362400001419
为((N-1)2-1)2×1维全1矢量，[·]i表示矩阵第i行元素。
在步骤S140中，根据重构矢量
Figure BDA00037822362400001420
和重构矢量
Figure BDA00037822362400001421
计算估计矩阵
Figure BDA00037822362400001422
的步骤包括：采用下面的公式(26)计算估计矩阵
Figure BDA00037822362400001423
Figure BDA00037822362400001424
其中，Φ2为第三子阵A2的旋转矩阵，
Figure BDA00037822362400001425
为全1矢量，[·]i表示矩阵第i行元素。
在步骤S150中，根据估计矩阵
Figure BDA00037822362400001426
和估计矩阵
Figure BDA00037822362400001427
使用元素平均法、相位平均法、最小二乘法中的任意一种方法计算相位估计
Figure BDA00037822362400001428
和相位估计
Figure BDA00037822362400001429
相位估计
Figure BDA00037822362400001430
和相位估计
Figure BDA00037822362400001431
包括解扩信号的DOA信息，DOA信息包括俯仰角和/或方位角。
在一些实施例中，根据估计矩阵
Figure BDA00037822362400001432
使用元素平均(Element average，EA)法计算相位估计
Figure BDA00037822362400001433
的步骤包括：采用下面的公式(27)、(28)进行计算：
Figure BDA0003782236240000151
Figure BDA0003782236240000152
其中，
Figure BDA0003782236240000153
表示使用元素平均法得到的第二子阵A1的旋转矩阵估计值，N表示天线阵列的维度，
Figure BDA0003782236240000154
表示使用元素平均法得到的相位估计，angle(·)函数用于求解相位角。在公式(27)对
Figure BDA0003782236240000155
中所有元素取均值，在公式(28)对
Figure BDA0003782236240000156
取相位。
类似的，根据估计矩阵
Figure BDA0003782236240000157
使用元素平均法计算相位估计
Figure BDA0003782236240000158
的步骤可以采用上述公式(27)、(28)进行计算，在此不再赘述。
在一些实施例中，根据估计矩阵
Figure BDA0003782236240000159
使用相位平均(Phase Average，PA)法计算相位估计
Figure BDA00037822362400001510
的步骤包括：采用下面的公式(29)进行计算：
Figure BDA00037822362400001511
其中，
Figure BDA00037822362400001512
表示使用相位平均法得到的相位估计，N表示天线阵列的维度，angle(·)函数用于求解相位角。
类似的，根据估计矩阵
Figure BDA00037822362400001513
使用相位平均法计算相位估计
Figure BDA00037822362400001514
的步骤可以采用上述公式(29)进行计算，在此不再赘述。
在一些实施例中，根据估计矩阵
Figure BDA00037822362400001515
使用最小二乘(Least Squares，LS)法计算相位估计
Figure BDA00037822362400001516
的步骤包括：采用下面的公式(30)、(31)、(32)进行计算：
Figure BDA00037822362400001517
Figure BDA00037822362400001518
Figure BDA00037822362400001519
其中，
Figure BDA00037822362400001520
表示第二子阵A1的旋转矩阵的共轭，
Figure BDA00037822362400001521
表示重构矢量
Figure BDA00037822362400001522
的伪逆，
Figure BDA00037822362400001523
表示重构矢量
Figure BDA00037822362400001524
的共轭转置，
Figure BDA00037822362400001525
表示使用最小二乘法得到的相位估计，angle(·)函数用于求解相位角。根据前述公式(20)可得到公式(30)，根据最小二乘法的思想对公式(30)求伪逆可得到公式(31)，由于只需要旋转矩阵Φ1的相位信息，因此只需计算公式(31)中的
Figure BDA00037822362400001526
部分即可，由此得到公式(32)。
类似的，根据估计矩阵
Figure BDA00037822362400001527
使用最小二乘法计算相位估计
Figure BDA00037822362400001528
的步骤可以采用上述公式(30)、(31)、(32)进行计算，在此不再赘述。
本申请将天线阵列划分为三个子阵，根据子阵间的相关关系，使用点除或元素平均法、相位平均法、最小二乘法中的任意一种方法代替标准ESPRIT算法中系统计算开销较大的矩阵特征分解操作，在保证测向精度的同时，计算复杂度大大降低。
接下来由相位估计
Figure BDA0003782236240000161
和相位估计
Figure BDA0003782236240000162
对用户来波方向进行估计。将方向余弦用公式(33)表示：
Figure BDA0003782236240000163
其中，ux表示天线阵列中x轴方向的方向余弦，uy表示天线阵列中y轴方向的方向余弦。
在步骤S160中，根据相位估计
Figure BDA0003782236240000164
和相位估计
Figure BDA0003782236240000165
分别计算天线阵列中x轴方向的方向余弦估计
Figure BDA0003782236240000166
和y轴方向的方向余弦估计
Figure BDA0003782236240000167
的步骤包括：采用下面的公式(34)计算方向余弦估计
Figure BDA0003782236240000168
和方向余弦估计
Figure BDA0003782236240000169
Figure BDA00037822362400001610
其中，λ表示接收信号的波长，d表示天线阵列的阵元间距。根据前述公式(12)和(33)可得公式(34)。
在步骤S170中，根据方向余弦估计
Figure BDA00037822362400001611
和方向余弦估计
Figure BDA00037822362400001612
计算俯仰角估计
Figure BDA00037822362400001613
和/或方位角估计
Figure BDA00037822362400001614
的步骤包括：采用下面的公式计算俯仰角估计
Figure BDA00037822362400001615
和/或方位角估计
Figure BDA00037822362400001616
Figure BDA00037822362400001617
Figure BDA00037822362400001618
在实际应用中可根据需要计算俯仰角估计
Figure BDA00037822362400001619
和方位角估计
Figure BDA00037822362400001620
还可以只计算俯仰角估计
Figure BDA00037822362400001621
或方位角估计
Figure BDA00037822362400001622
本申请不做限制。
接下来，进行计算复杂度分析。
为方便描述，将各子阵阵元数简记为Msub＝(N-1)×(N-1)，使用元素平均法、相位平均法、最小二乘法估计相位，将本申请的简化ESPRIT算法分别记为S-ESPRIT-EA、S-ESPRIT-PA、以及S-ESPRIT-LS算法。对于S-ESPRIT-EA算法，先由
Figure BDA0003782236240000171
和
Figure BDA0003782236240000172
使用点除方法得到
Figure BDA0003782236240000173
与
Figure BDA0003782236240000174
再对各元素求平均，对所得均值取相位，其计算复杂度约为
Figure BDA0003782236240000175
S-ESPRIT-PA在得到
Figure BDA0003782236240000176
与
Figure BDA0003782236240000177
后，先对各元素取相位，再对相位值进行平均，计算复杂度增加为
Figure BDA0003782236240000178
而S-ESPRIT-LS算法直接使用重构矢量，在矢量乘法操作后取相位，计算复杂度约为
Figure BDA0003782236240000179
根据参考文献，标准ESPRIT算法的计算复杂度约为
Figure BDA00037822362400001710
综上所示，本申请的二维DOA估计方法相比标准ESPRIT算法具有更低的计算复杂度，且相差约Msub的一个数量级。
接下来，根据本申请二维DOA估计方法的仿真结果进行说明。
在仿真中，K个等功率用户信号入射到阵列，用户信息符号数为S，使用扩频增益G＝64的Gold序列进行扩频。不失一般性，仿真结果均以第1个用户的DOA估计为例，该用户信号入射方向的俯仰角和方位角分别为θ1＝30°,
Figure BDA00037822362400001711
其余(K-1)个用户信号入射方向的俯仰角θk,k＝2,3,...,K服从区间[0°,90°]内的均匀分布，方位角
Figure BDA00037822362400001712
服从区间[0°,180°]内的均匀分布。定义DOA估计误差Δα作为测向精度的评价指标，用公式(37)表示：
Figure BDA00037822362400001713
其中，NM为蒙特卡洛次数，后续仿真中取NM＝200，α(α代表θ或
Figure BDA00037822362400001714
)为DOA真值，
Figure BDA00037822362400001715
为第i次的DOA估计值。
考虑阵元分布为N×N的矩形栅格阵列，阵元间距d＝λ/2。经码匹配滤波后，使用本申请二维DOA估计方法，即简化ESPRIT算法(S-ESPRIT-EA，S-ESPRIT-PA和S-ESPRIT-LS)以及标准ESPRIT算法进行二维DOA估计，比较各算法在不同仿真条件下的测向性能。
1)不同信噪比下简化ESPRIT算法测向性能。
图4是本申请一实施例的二维DOA估计方法仿真结果的俯仰角估计误差与输入SNR关系曲线示意图；图5是本申请一实施例的二维DOA估计方法仿真结果的方位角估计误差与输入SNR关系曲线示意图。
在此，验证本申请的简化ESPRIT算法在不同输入信噪比(Signal-to-NoiseRatio，SNR)下的测向性能。参考图4和图5所示，假设用户数K＝4，用户信息符号数S＝200，矩形栅格阵列行/列数N＝6。输入SNR变化范围为[-25dB,15dB]。图4和图5分别为第1个用户信号入射方向的俯仰角θ1和方位角
Figure BDA0003782236240000181
估计误差与输入SNR的关系曲线。可以看出，本申请的简化ESPRIT算法在降低了计算复杂度的同时，算法性能与标准ESPRIT算法相差无几。在输入SNR＞-10dB时，各算法的测向误差曲线几乎完全重合。信噪比较低时，S-ESPRIT-EA算法测向误差增大较为明显，而S-ESPRIT-LS算法与标准ESPRIT算法性能最为接近。
2)不同用户信息符号数下简化ESPRIT算法测向性能。
图6是本申请一实施例的二维DOA估计方法仿真结果的俯仰角估计误差与用户信息符号数关系曲线示意图；图7是本申请一实施例的二维DOA估计方法仿真结果的方位角估计误差与用户信息符号数关系曲线示意图。
在此验证本申请简化ESPRIT算法在不同用户信息符号数下的测向性能。参考图6和图7所示，假设用户数K＝4，矩形栅格阵列行/列数N＝6，输入SNR固定为-10dB,用户信息符号数S的变化范围为[50,300]。θ1和
Figure BDA0003782236240000182
的估计误差与符号数S的关系曲线分别如图6、图7所示。可以看出，随着符号数的增加，本申请简化ESPRIT算法与标准ESPRIT算法的俯仰角与方位角测向误差均有所下降。由于上述算法都需要使用信号解扩后的协方差矩阵进行测向，所以对用户信息符号数变化的敏感程度也几乎一致。
3)不同用户数下简化ESPRIT算法测向性能。
图8是本申请一实施例的二维DOA估计方法仿真结果的俯仰角估计误差与用户数关系曲线示意图；图9是本申请一实施例的二维DOA估计方法仿真结果的方位角估计误差与用户数关系曲线示意图。
在此验证本申请简化ESPRIT算法在不同用户数下的测向性能。参考图8和图9所示，假设矩形栅格阵列行/列数N＝4，用户信息符号数S＝200，输入SNR固定为-10dB,用户数分别取K＝1,5,15,45,60,64,65,70。θ1和
Figure BDA0003782236240000183
的估计误差与用户数变化关系如图8、图9所示。由图可知，当用户数小于60时，本申请简化ESPRIT算法与标准ESPRIT算法测向性能受用户数变化的影响不明显，且俯仰角与方位角测向误差均在0.2°以内。当用户数大于64时，各算法的测向误差显著增大。对于阵元分布为4×4的均匀平面阵列，用户数在60左右时已明显超出阵列的自由度，由于各用户在码域接近于正交，对接收信号经码匹配滤波后进行逐用户的DOA估计，可以更好地适应天线阵列的过载状态。
本申请还包括一种二维DOA估计系统，包括天线阵列、存储器和处理器。其中，该天线阵列用于接收信号；存储器用于存储可由处理器执行的指令；处理器用于执行该指令以实现前文所述的二维DOA估计方法。
图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所执行的可能的程序指令。处理器执行这些指令以实现方法的主要部分。处理器处理的结果通过通信端口传给用户设备，在用户界面上显示。
上述的二维DOA估计方法可以实施为计算机程序，保存在硬盘1006中，并可加载到处理器1002中执行，以实施本申请的二维DOA估计方法。
本申请还包括一种存储有计算机程序代码的计算机可读介质，该计算机程序代码在由处理器执行时实现前文所述的二维DOA估计方法。
二维DOA估计方法实施为计算机程序时，也可以存储在计算机可读存储介质中作为制品。例如，计算机可读存储介质可以包括但不限于磁存储设备(例如，硬盘、软盘、磁条)、光盘(例如，压缩盘(CD)、数字多功能盘(DVD))、智能卡和闪存设备(例如，电可擦除可编程只读存储器(EPROM)、卡、棒、键驱动)。此外，本文描述的各种存储介质能代表用于存储信息的一个或多个设备和/或其它机器可读介质。术语“机器可读介质”可以包括但不限于能存储、包含和/或承载代码和/或指令和/或数据的无线信道和各种其它介质(和/或存储介质)。
应该理解，上文所描述的实施例仅是示意。本文描述的实施例可在硬件、软件、固件、中间件、微码或者其任意组合中实现。对于硬件实现，处理器可以在一个或者多个特定用途集成电路(ASIC)、数字信号处理器(DSP)、数字信号处理设备(DSPD)、可编程逻辑器件(PLD)、现场可编程门阵列(FPGA)、处理器、控制器、微控制器、微处理器和/或设计为执行本文所述功能的其它电子单元或者其结合内实现。
本申请的一些方面可以完全由硬件执行、可以完全由软件(包括固件、常驻软件、微码等)执行、也可以由硬件和软件组合执行。以上硬件或软件均可被称为“数据块”、“模块”、“引擎”、“单元”、“组件”或“系统”。处理器可以是一个或多个专用集成电路(ASIC)、数字信号处理器(DSP)、数字信号处理器件(DAPD)、可编程逻辑器件(PLD)、现场可编程门阵列(FPGA)、处理器、控制器、微控制器、微处理器或者其组合。此外，本申请的各方面可能表现为位于一个或多个计算机可读介质中的计算机产品，该产品包括计算机可读程序编码。例如，计算机可读介质可包括，但不限于，磁性存储设备(例如，硬盘、软盘、磁带……)、光盘(例如，压缩盘CD、数字多功能盘DVD……)、智能卡以及闪存设备(例如，卡、棒、键驱动器……)。
计算机可读介质可能包含一个内含有计算机程序编码的传播数据信号，例如在基带上或作为载波的一部分。该传播信号可能有多种表现形式，包括电磁形式、光形式等等、或合适的组合形式。计算机可读介质可以是除计算机可读存储介质之外的任何计算机可读介质，该介质可以通过连接至一个指令执行系统、装置或设备以实现通讯、传播或传输供使用的程序。位于计算机可读介质上的程序编码可以通过任何合适的介质进行传播，包括无线电、电缆、光纤电缆、射频信号、或类似介质、或任何上述介质的组合。
上文已对基本概念做了描述，显然，对于本领域技术人员来说，上述申请披露仅仅作为示例，而并不构成对本申请的限定。虽然此处并没有明确说明，本领域技术人员可能会对本申请进行各种修改、改进和修正。该类修改、改进和修正在本申请中被建议，所以该类修改、改进、修正仍属于本申请示范实施例的精神和范围。
同时，本申请使用了特定词语来描述本申请的实施例。如“一个实施例”、“一实施例”、和/或“一些实施例”意指与本申请至少一个实施例相关的某一特征、结构或特点。因此，应强调并注意的是，本说明书中在不同位置两次或多次提及的“一实施例”或“一个实施例”或“一替代性实施例”并不一定是指同一实施例。此外，本申请的一个或多个实施例中的某些特征、结构或特点可以进行适当的组合。
一些实施例中使用了描述成分、属性数量的数字，应当理解的是，此类用于实施例描述的数字，在一些示例中使用了修饰词“大约”、“近似”或“大体上”来修饰。除非另外说明，“大约”、“近似”或“大体上”表明所述数字允许有±20％的变化。相应地，在一些实施例中，说明书和权利要求中使用的数值参数均为近似值，该近似值根据个别实施例所需特点可以发生改变。在一些实施例中，数值参数应考虑规定的有效数位并采用一般位数保留的方法。尽管本申请一些实施例中用于确认其范围广度的数值域和参数为近似值，在具体实施例中，此类数值的设定在可行范围内尽可能精确。

Claims (13)

1.一种二维DOA估计方法，所述二维DOA估计方法应用于CDMA通信体制下，其特征在于，包括：
将天线阵列的接收信号进行解扩预处理，获得解扩信号，将所述天线阵列划分为第一子阵A0、第二子阵A1和第三子阵A2，计算所述第一子阵A0解扩信号的自相关矩阵
Figure FDA0003782236230000011
所述第一子阵A0解扩信号与所述第二子阵A1解扩信号的互相关矩阵
Figure FDA0003782236230000012
所述第一子阵A0解扩信号与所述第三子阵A2解扩信号的互相关矩阵
Figure FDA0003782236230000013
剔除所述自相关矩阵
Figure FDA0003782236230000014
和所述互相关矩阵
Figure FDA0003782236230000015
的主对角线元素与左下副对角线元素，并将所述自相关矩阵
Figure FDA0003782236230000016
和所述互相关矩阵
Figure FDA0003782236230000017
的剩余元素分别重构为列矢量形式，获得重构矢量
Figure FDA0003782236230000018
和重构矢量
Figure FDA0003782236230000019
以及，剔除所述自相关矩阵
Figure FDA00037822362300000110
和所述互相关矩阵
Figure FDA00037822362300000111
的主对角线元素与左下角((N-1)(N-2))阶矩阵的主对角线元素，并将所述自相关矩阵
Figure FDA00037822362300000112
和所述互相关矩阵
Figure FDA00037822362300000113
的剩余元素分别重构为列矢量形式，获得重构矢量
Figure FDA00037822362300000114
和重构矢量
Figure FDA00037822362300000115
根据所述重构矢量
Figure FDA00037822362300000116
和所述重构矢量
Figure FDA00037822362300000117
计算估计矩阵
Figure FDA00037822362300000118
根据所述重构矢量
Figure FDA00037822362300000119
和所述重构矢量
Figure FDA00037822362300000120
计算估计矩阵
Figure FDA00037822362300000121
根据所述估计矩阵
Figure FDA00037822362300000122
和所述估计矩阵
Figure FDA00037822362300000123
使用元素平均法、相位平均法、最小二乘法中的任意一种方法计算相位估计
Figure FDA00037822362300000124
和相位估计
Figure FDA00037822362300000125
所述相位估计
Figure FDA00037822362300000126
和所述相位估计
Figure FDA00037822362300000127
包括所述解扩信号的DOA信息，所述DOA信息包括俯仰角和/或方位角；
根据所述相位估计
Figure FDA00037822362300000128
和所述相位估计
Figure FDA00037822362300000129
分别计算所述天线阵列中x轴方向的方向余弦估计
Figure FDA00037822362300000130
和y轴方向的方向余弦估计
Figure FDA00037822362300000131
以及
根据所述方向余弦估计
Figure FDA00037822362300000132
和所述方向余弦估计
Figure FDA00037822362300000133
计算所述俯仰角估计
Figure FDA00037822362300000134
和/或所述方位角估计
Figure FDA00037822362300000135
2.如权利要求1所述的二维DOA估计方法，其特征在于，将天线阵列的接收信号进行解扩预处理，获得解扩信号的步骤包括：采用下面的公式计算所述解扩信号：
Figure FDA00037822362300000136
Figure FDA0003782236230000021
其中，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个符号周期的延拓；
yk(i)表示所述解扩信号，所述解扩信号具有等效单源特性，k表示第k个用户，k'表示第k个用户以外的其他用户，
Figure FDA0003782236230000022
表示第k个用户的扩频序列与其他用户扩频序列的互相关函数；等效单源信号sk(i)＝Gσkdk(i)，等效噪声
Figure FDA0003782236230000023
3.如权利要求1所述的二维DOA估计方法，其特征在于，所述天线阵列为N×N维的阵列，所述第一子阵A0、所述第二子阵A1和所述第三子阵A2均为(N-1)×(N-1)维的方阵，将所述天线阵列划分为第一子阵A0、第二子阵A1和第三子阵A2的步骤包括：将所述天线阵列剔除第N行阵元与第N列阵元获得所述第一子阵A0；将所述天线阵列剔除第N行阵元与第1列阵元获得所述第二子阵A1；将所述天线阵列剔除第1行阵元与第N列阵元获得所述第三子阵A2。
4.如权利要求1所述的二维DOA估计方法，其特征在于，采用下面的公式分别计算所述自相关矩阵
Figure FDA0003782236230000024
所述互相关矩阵
Figure FDA0003782236230000025
和所述互相关矩阵
Figure FDA0003782236230000026
Figure FDA0003782236230000027
Figure FDA0003782236230000028
Figure FDA0003782236230000029
其中，S表示用户的符号信息数，x0(i)表示所述第一子阵A0解扩信号，
Figure FDA00037822362300000210
表示所述第一子阵A0解扩信号的共轭转置，
Figure FDA00037822362300000211
表示所述第二子阵A1解扩信号的共轭转置，
Figure FDA0003782236230000031
表示所述第三子阵A2解扩信号的共轭转置。
5.如权利要求1所述的二维DOA估计方法，其特征在于，根据所述重构矢量
Figure FDA0003782236230000032
和所述重构矢量
Figure FDA0003782236230000033
计算估计矩阵
Figure FDA0003782236230000034
的步骤包括：采用下面的公式计算所述估计矩阵
Figure FDA0003782236230000035
Figure FDA0003782236230000036
其中，Φ1为所述第二子阵A1的旋转矩阵，
Figure FDA0003782236230000037
为全1矢量，[·]i表示矩阵第i行元素。
6.如权利要求1所述的二维DOA估计方法，其特征在于，根据所述重构矢量
Figure FDA0003782236230000038
和所述重构矢量
Figure FDA0003782236230000039
计算估计矩阵
Figure FDA00037822362300000310
的步骤包括：采用下面的公式计算所述估计矩阵
Figure FDA00037822362300000311
Figure FDA00037822362300000312
其中，Φ2为所述第三子阵A2的旋转矩阵，
Figure FDA00037822362300000313
为全1矢量，[·]i表示矩阵第i行元素。
7.如权利要求5所述的二维DOA估计方法，其特征在于，根据所述估计矩阵
Figure FDA00037822362300000314
使用元素平均法计算相位估计
Figure FDA00037822362300000315
的步骤包括：采用下面的公式进行计算：
Figure FDA00037822362300000316
Figure FDA00037822362300000317
其中，
Figure FDA00037822362300000318
表示使用元素平均法得到的所述第二子阵A1的旋转矩阵估计值，N表示所述天线阵列的维度，
Figure FDA00037822362300000319
表示使用元素平均法得到的相位估计，angle(·)函数用于求解相位角。
8.如权利要求5所述的二维DOA估计方法，其特征在于，根据所述估计矩阵
Figure FDA00037822362300000320
使用相位平均法计算相位估计
Figure FDA00037822362300000321
的步骤包括：采用下面的公式进行计算：
Figure FDA00037822362300000322
其中，
Figure FDA00037822362300000323
表示使用相位平均法得到的相位估计，N表示所述天线阵列的维度，angle(·)函数用于求解相位角。
9.如权利要求5所述的二维DOA估计方法，其特征在于，根据所述估计矩阵
Figure FDA0003782236230000041
使用最小二乘法计算相位估计
Figure FDA0003782236230000042
的步骤包括：采用下面的公式进行计算：
Figure FDA0003782236230000043
Figure FDA0003782236230000044
Figure FDA0003782236230000045
其中，
Figure FDA0003782236230000046
表示所述第二子阵A1的旋转矩阵的共轭，
Figure FDA0003782236230000047
表示所述重构矢量
Figure FDA0003782236230000048
的伪逆，
Figure FDA0003782236230000049
表示所述重构矢量
Figure FDA00037822362300000410
的共轭转置，
Figure FDA00037822362300000411
表示使用最小二乘法得到的相位估计，angle(·)函数用于求解相位角。
10.如权利要求1所述的二维DOA估计方法，其特征在于，根据所述相位估计
Figure FDA00037822362300000412
和所述相位估计
Figure FDA00037822362300000413
分别计算所述天线阵列中x轴方向的方向余弦估计
Figure FDA00037822362300000414
和y轴方向的方向余弦估计
Figure FDA00037822362300000415
的步骤包括：采用下面的公式计算所述方向余弦估计
Figure FDA00037822362300000416
和所述方向余弦估计
Figure FDA00037822362300000417
Figure FDA00037822362300000418
其中，λ表示所述接收信号的波长，d表示所述天线阵列的阵元间距。
11.如权利要求10所述的方法，其特征在于，根据所述方向余弦估计
Figure FDA00037822362300000419
和所述方向余弦估计
Figure FDA00037822362300000420
计算所述俯仰角估计
Figure FDA00037822362300000421
和/或所述方位角估计
Figure FDA00037822362300000422
的步骤包括：采用下面的公式计算所述俯仰角估计
Figure FDA00037822362300000423
和/或所述方位角估计
Figure FDA00037822362300000424
Figure FDA00037822362300000425
Figure FDA00037822362300000426
12.一种二维DOA估计系统，其特征在于，包括：
天线阵列，所述天线阵列用于接收信号；
存储器，用于存储可由处理器执行的指令；
处理器，用于执行所述指令以实现如1-11任一项所述的二维DOA估计方法。
13.一种存储有计算机程序代码的计算机可读介质，其特征在于，所述计算机程序代码在由处理器执行时实现如1-11任一项所述的二维DOA估计方法。
CN202210932982.2A 2022-08-04 2022-08-04 二维doa估计方法、系统及计算机可读介质 Pending CN115308680A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210932982.2A CN115308680A (zh) 2022-08-04 2022-08-04 二维doa估计方法、系统及计算机可读介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210932982.2A CN115308680A (zh) 2022-08-04 2022-08-04 二维doa估计方法、系统及计算机可读介质

Publications (1)

Publication Number Publication Date
CN115308680A true CN115308680A (zh) 2022-11-08

Family

ID=83858335

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210932982.2A Pending CN115308680A (zh) 2022-08-04 2022-08-04 二维doa估计方法、系统及计算机可读介质

Country Status (1)

Country Link
CN (1) CN115308680A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115291160A (zh) * 2022-08-04 2022-11-04 中国科学院微小卫星创新研究院 二维doa估计方法、系统及计算机可读介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4965732A (en) * 1985-11-06 1990-10-23 The Board Of Trustees Of The Leland Stanford Junior University Methods and arrangements for signal reception and parameter estimation
CN108872930A (zh) * 2018-08-28 2018-11-23 天津大学 扩展孔径二维联合对角化doa估计方法
CN109116295A (zh) * 2018-08-31 2019-01-01 上海微小卫星工程中心 基于相控阵选取基线的无源测向算法
CN113341371A (zh) * 2021-05-31 2021-09-03 电子科技大学 一种基于l阵和二维esprit算法的doa估计方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4965732A (en) * 1985-11-06 1990-10-23 The Board Of Trustees Of The Leland Stanford Junior University Methods and arrangements for signal reception and parameter estimation
CN108872930A (zh) * 2018-08-28 2018-11-23 天津大学 扩展孔径二维联合对角化doa估计方法
CN109116295A (zh) * 2018-08-31 2019-01-01 上海微小卫星工程中心 基于相控阵选取基线的无源测向算法
CN113341371A (zh) * 2021-05-31 2021-09-03 电子科技大学 一种基于l阵和二维esprit算法的doa估计方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
王伟伟;梁广真;: "基于均分阵元的空域平滑技术研究", 现代电子技术, no. 09, 1 May 2013 (2013-05-01) *
许培培: "一种改进的卫星空时DOA矩阵算法", 计算机工程, 20 October 2010 (2010-10-20) *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115291160A (zh) * 2022-08-04 2022-11-04 中国科学院微小卫星创新研究院 二维doa估计方法、系统及计算机可读介质

Similar Documents

Publication Publication Date Title
EP1253434B1 (en) Method for estimating a direction of arrival
Zhang et al. Estimation of two-dimensional direction-of-arrival for uncorrelated and coherent signals with low complexity
Dakulagi et al. Efficient coherent direction-of-arrival estimation and realization using digital signal processor
Zhai et al. DOA estimation of noncircular signals for unfolded coprime linear array: Identifiability, DOF and algorithm (May 2018)
CN115308680A (zh) 二维doa估计方法、系统及计算机可读介质
CN114884841A (zh) 基于高阶统计和非均匀阵列的欠定参数联合估计方法
CN113296049A (zh) 互质阵列脉冲环境下非圆信号的共轭增广doa估计方法
Qi et al. Time-frequency DOA estimation of chirp signals based on multi-subarray
CN109946663B (zh) 一种线性复杂度的Massive MIMO目标空间方位估计方法和装置
CN111679123B (zh) 适用于多模调制体制的符号边沿与频率估计方法与系统
Qi et al. An improved multiple-Toeplitz matrices reconstruction algorithm for DOA estimation of coherent signals
CN113238184B (zh) 一种基于非圆信号的二维doa估计方法
Peng et al. A novel weighted spatial smoothing DOA estimation algorithm for coherent signals
CN113391266B (zh) 基于非圆多嵌套阵降维子空间数据融合的直接定位方法
CN113406562B (zh) 北斗和超宽带系统中一种toa与doa联合估计降维方法
CN115291160A (zh) 二维doa估计方法、系统及计算机可读介质
CN115545067A (zh) 一种基于改进music算法的多径时延估计方法
CN114488064A (zh) 一种距离-速度联合估计方法
CN109683128B (zh) 冲击噪声环境下的单快拍测向方法
Zhai et al. DOA estimation of noncircular signals for coprime linear array via polynomial root-finding technique
Jiang et al. The DOA estimation for a mixture of uncorrelated and coherent sources via decomposing the coprime array
Xu et al. Two‐Dimensional Direction‐of‐Arrival Estimation of Noncircular Signals in Coprime Planar Array with High Degree of Freedom
Varum et al. Detect and Pointing Algorithms Performance for a 2D Adaptive Antenna Array
Tayem et al. Capon root-MUSIC-like direction of arrival estimation based on real data
Mismar et al. DOA estimation by controlling the nulls of the antenna array factor

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination