CN109239651B - 互质面阵下的二维doa跟踪方法 - Google Patents

互质面阵下的二维doa跟踪方法 Download PDF

Info

Publication number
CN109239651B
CN109239651B CN201810815417.1A CN201810815417A CN109239651B CN 109239651 B CN109239651 B CN 109239651B CN 201810815417 A CN201810815417 A CN 201810815417A CN 109239651 B CN109239651 B CN 109239651B
Authority
CN
China
Prior art keywords
array
matrix
angle
doa
sub
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.)
Active
Application number
CN201810815417.1A
Other languages
English (en)
Other versions
CN109239651A (zh
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.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
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 Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN201810815417.1A priority Critical patent/CN109239651B/zh
Publication of CN109239651A publication Critical patent/CN109239651A/zh
Application granted granted Critical
Publication of CN109239651B publication Critical patent/CN109239651B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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/74—Multi-channel systems specially adapted for direction-finding, i.e. having a single antenna system capable of giving simultaneous indications of the directions of different signals
    • 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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D30/00—Reducing energy consumption in communication networks
    • Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种互质面阵下的二维DOA跟踪方法，利用PASTd方法实时更新接收数据的信号子空间，然后基于互质面阵下的2D‑ESPRIT算法得到精确的DOA估计，并且通过联合子面阵间的角度关系来消除角度模糊问题，从而实现互质面阵下的低复杂度的DOA实时跟踪。本发明无需分解特征值和构造接收数据的协方差矩阵就能实现2D‑DOA信息的实时跟踪，复杂度低；对比互质阵下基于PAST进行DOA跟踪的方法，本发明具有更好的跟踪性能。

Description

互质面阵下的二维DOA跟踪方法
技术领域
本发明涉及阵列信号处理技术领域，尤其是一种互质面阵下的二维DOA跟踪方法。
背景技术
相对于传统面阵，互质面阵的主要思想是通过联合两个满足互质关系的子面阵进行空间谱估计，不仅能有效解决测向模糊问题，还可以获得较大的阵列孔径，显著提高波束形成和DOA估计中的自由度，因此得到广泛应用。传统的经典超分辨率DOA估计算法需要对接收信号的协方差矩阵进行特征值分解或者奇异值分解，计算量太大，难以满足DOA实时跟踪的需求，因此研究DOA跟踪算法变得很有实际意义。针对这个问题，国内外学者对其中的难点提出了诸多方法。其中比较有代表性的算法是子空间跟踪算法，包括PAST、PASTd、双迭代最小二乘(Bi-Iterative Least Squares，Bi-ILS)算法、双边迭代奇异值分解(Bi-iteration Singular Value Decomposition，Bi-SVD)算法等。本发明将互质面阵与PASTd算法相结合，提出互质面阵下的DOA实时跟踪算法，实现DOA角度的自动配对与跟踪，由于无需进行协方差矩阵的构造和特征值分解，算法复杂度低。
发明内容
本发明所要解决的技术问题在于，提供一种互质面阵下的二维DOA跟踪方法，能够实现互质面阵下的任意动态目标的低复杂度的DOA实时跟踪。
为解决上述技术问题，本发明提供一种互质面阵下的二维DOA跟踪方法，包括如下步骤：
(1)通过联合互质面阵中的两个子面阵构建阵列接收信号的数据模型；
(2)利用PASTd方法实时更新接收数据的信号子空间；
(3)基于阵列接收数据的信号子空间，利用互质面阵下的二维ESPRIT算法得到DOA估计，并通过联合阵列子阵间的角度关系实现解模糊。
优选的，步骤(1)具体为：
阵元总数为M2+N2-1的互质面阵，由阵元数分别为N×N和M×M的两个均匀子面阵构成，且两个子阵的阵元间距分别为Mλ/2和Nλ/2；对于子面阵i，其在X轴和Y轴上阵元的方向矩阵可以分别表示为
Figure BDA0001740197600000011
Figure BDA0001740197600000021
其中，方向向量/>
Figure BDA0001740197600000022
Figure BDA0001740197600000023
di为阵元间距；因此，子面阵i的接收信号xi(t)可以表示为
Figure BDA0001740197600000024
其中，s(t)为信源矢量，ni(t)为子面阵i的噪声，
Figure BDA0001740197600000025
和о分别表示Kronecker积和Khatri-Rao积。
优选的，步骤(2)具体为：
(21)选择适当的初始值λk(0)，W(0)；
(22)对每一个t＝1,2,…,J(J为快拍数)，使得x1(t)＝X(t)；
(23)对每一个k＝1,2,…,K(K为信源数)，分别更新以下变量：阵列接收数
Figure BDA0001740197600000026
特征值/>
Figure BDA0001740197600000027
特征向量/>
Figure BDA0001740197600000028
Figure BDA0001740197600000029
以及阵列接收数据xk+1(t)＝xk(t)-Wk(t)yk(t)；
(24)当步骤(23)中的k＝K后，使得t＝t+1，再次从步骤(22)开始计算。
优选的，步骤(3)具体为：
在互质面阵中，对于子面阵i，按下式构造Ai1，Ai2
Figure BDA00017401976000000215
Figure BDA00017401976000000216
其中，Mi＝M或Mi＝N，Ai2＝Ai1Φiy，旋转因子
Figure BDA00017401976000000210
将信号子空间Eis分解成Eix＝Es(1:Mi(Mi-1),:)，
Figure BDA00017401976000000211
其中Eix为Eis的1到Mi(Mi-1)行，Eiy为Eis的Mi+1到/>
Figure BDA00017401976000000212
行，Eix、Eiy也可表达为Eix＝Ai1Ti，Eiy＝Ai1ΦiyTi，矩阵/>
Figure BDA00017401976000000213
满秩，进而有
Figure BDA00017401976000000214
其中，
Figure BDA0001740197600000031
则Φiy的对角线的元素就是Ψi的特征值，由最小二乘法则可以得到Ψi的估计
Figure BDA0001740197600000032
将
Figure BDA0001740197600000033
特征值分解可以获得/>
Figure BDA0001740197600000034
然后根据/>
Figure BDA0001740197600000035
的特征矢量，就能得到T的估计结果/>
Figure BDA0001740197600000036
在不考虑噪声的条件下有
Figure BDA0001740197600000037
Πi表示一个置换矩阵；因为
Figure BDA0001740197600000038
和/>
Figure BDA0001740197600000039
的特征值是一致的，所以将/>
Figure BDA00017401976000000310
特征值分解就可以得到/>
Figure BDA00017401976000000311
Figure BDA00017401976000000312
的估计为
Figure BDA00017401976000000313
其中
Figure BDA00017401976000000314
为/>
Figure BDA00017401976000000315
的第k个特征值；
同样对矩阵Eis重构，可得
Figure BDA00017401976000000316
类似地可以将E'is分解成矩阵E'ix＝E'is(1:Mi(Mi-1),:)和矩阵/>
Figure BDA00017401976000000317
构造矩阵
Figure BDA00017401976000000318
那么，/>
Figure BDA00017401976000000319
在不考虑噪声的条件下有/>
Figure BDA00017401976000000320
其中Πi表示一个置换矩阵，类似于uik，可以得到/>
Figure BDA00017401976000000321
的估计值
Figure BDA00017401976000000322
其中
Figure BDA00017401976000000323
表示的是矩阵(E'ix)+E'iy第k个特征值；
通过联合两个子阵间的角度关系，可以消除互质面阵中角度估计的模糊值，即两个子阵相同的估计值即为真实角度，其余值为模糊值。
根据uik和vik的表达式可知，uik和vik的估计值有相同的列模糊，因此估计角度可以自动配对，角度配对完成后，根据以下公式即可得到信源仰角和方位角的估计值
Figure BDA00017401976000000324
Figure BDA00017401976000000325
本发明的有益效果为：本发明无需分解特征值和构造接收数据的协方差矩阵就能实现2D-DOA信息的实时跟踪，复杂度低；对比互质阵下基于PAST进行DOA跟踪的方法，本发明具有更好的跟踪性能。
附图说明
图1为本发明的阵列结构示意图。
图2为本发明的PASTd算法流程示意图。
图3为本发明实现的信源方位角的60s跟踪图。
图4为本发明实现的信源俯仰角的60s跟踪图。
图5为本发明实现的信源DOA的60s跟踪图。
图6为本发明在不同阵元数下的跟踪性能对比示意图。
图7为本发明与PAST算法在不同信噪比下的跟踪性能对比示意图。
具体实施方式
一、数据模型
互质面阵的阵列结构如图1所示，与传统均匀面阵相比，互质面阵可以分成两个子阵，N表示第一个子阵在X轴和Y轴方向的阵元数，M表示第二个子阵在X轴和Y轴方向的阵元数。子阵一是由阵元数为N×N的均匀面阵构成，子阵二由阵元数为M×M的均匀面阵构成。子阵阵元在原点位置重合，故互质面阵的阵元总数为M2+N2-1。其中子阵一的阵元间距为d1＝Mλ/2，子阵二的阵元间距为d2＝Nλ/2，各阵元的位置可表示为如下集合
Ds＝{(md1,nd1)|0≤m,n≤N-1}∪{(pd2,qd2)|0≤p,q≤M-1}
假设空间有K个非相干信源入射到上述互质面阵，θk,
Figure BDA0001740197600000041
分别代表入射信源的俯仰角和方位角。考虑面阵中的子面阵i，其在X轴和Y轴上阵元的方向矩阵可以分别表示为
Figure BDA0001740197600000042
其中，方向向量/>
Figure BDA0001740197600000043
Figure BDA0001740197600000044
因此，子面阵i的接收信号xi(t)可以表示为
Figure BDA0001740197600000051
其中，s(t)为信源矢量，ni(t)为子面阵i的噪声。
Figure BDA0001740197600000052
和/>
Figure BDA0001740197600000053
分别表示Kronecker积和Khatri-Rao积。
二、利用PASTd方法实时更新信号子空间
定义一个无约束代价函数
J(W)＝E{||x-WWHx||2}
＝E{(x-WWHx)H(x-WWHx)}
＝E{xHx}-2E{xHWWHx}+E{xHWWHWWHx}
可以发现
Figure BDA0001740197600000054
E{xHWWHx}＝tr(E{WHxxHW})＝tr(WHCW)
E{xHWWHWWHx}＝tr(E{WHWxxHWHW})＝tr(WHCWWHW)
其中C表示接收数据x的自相关矩阵。假设W秩为N，则J(W)可以表示为
J(W)＝tr(C)-2tr(WHCW)+tr(WHCWWHW)
由Bin Yang提出的定理可以知道当目标函数J(W)取极小值时，W的列空间等价于信号子空间。
在现实情况中，为了更新得到t时刻的子空间W(t)，需要利用t-1时刻的子空间W(t-1)以及t时刻的阵元数据x(t)。在此我们选择梯度下降法，选取下降梯度为
Figure BDA0001740197600000055
所以
Figure BDA0001740197600000056
其中，μ>0表示需要适当选择的步长值。将C(t)＝x(t)xH(t)，y(t)＝WH(t-1)x(t)代入上式，得到
W(t)＝W(t-1)-μ[2x(t)yH(t)-x(t)yH(t)
×WH(t-1)W(t-1)-W(t-1)WH(t-1)y(t)yH(t)]
由于J(W)取得极小值时，WH(t-1)W(t-1)＝I，可得
W(t)＝W(t-1)+μ[x(t)-W(t-1)y(t)]yH(t)
由于W(t)跟踪时变子空间的能力较差，算法收敛慢。为了解决这个问题可以定义一个新的指数加权函数
Figure BDA0001740197600000061
其中0<β≤1，表示遗忘因子，主要是保证在非稳定环境下将过去的数据降低权重，以保证跟踪的稳定性。当β＝1时对应常见的滑动窗口。
Figure BDA0001740197600000062
进一步，可以认为
y(i)＝WH(i-1)x(i)≈WH(t)x(i)
因此得到已经修正的目标函数
Figure BDA0001740197600000063
当目标函数全局最小时，W(t)可以用自相关矩阵Cyy(t)和互相关矩阵Cxy(t)来表示，minJ(W(t))最优解是Wiener滤波器，即
Figure BDA0001740197600000064
其中，Cyy(t)和Cxy(t)更新公式为：
Figure BDA0001740197600000065
Figure BDA0001740197600000066
PASTd方法主要利用紧缩技术通过一定顺序来估计主分量，首先在目标数K为1时通过PAST方法估计出当前最主要的特征向量Wk(t)，再把当前t时刻的数据xk(t)减去它在特征向量Wk(t)上的投影得到xk+1(t)，然后重复上述步骤计算出Wk+1(t),Wk+2(t),…
利用PASTd方法求解信源信号子空间的具体步骤如下：
1)选择适当的初始值λk(0)，W(0)；
2)对每一个t＝1,2,…,J(J为快拍数)，使得x1(t)＝X(t)；
3)对每一个k＝1,2,…,K(K为信源数)，分别更新一下变量：阵列接收数
Figure BDA0001740197600000071
特征值/>
Figure BDA0001740197600000072
特征向量
Figure BDA0001740197600000073
以及阵列接收数据xk+1(t)＝xk(t)-Wk(t)yk(t)；
4)当步骤(3)中的k＝K后，使得t＝t+1，再次从步骤(2)开始计算；
PASTd的算法流程图如图2所示，算法最后一步通过xk(t)减去C(t)的第k个特征向量Wk(t)达到算法紧缩。
三、利用2D-ESPRIT得到DOA信息
在互质面阵中，对于子面阵i，按下式构造Ai1，Ai2
Figure BDA0001740197600000074
Figure BDA0001740197600000075
其中，Mi＝M或Mi＝N，Ai2＝Ai1Φiy，旋转因子
Figure BDA0001740197600000076
将由PASTd得到的信号子空间Eis分解成Eix＝Es(1:Mi(Mi-1),:)，
Figure BDA0001740197600000077
其中Eix为Eis的1到Mi(Mi-1)行，Eiy为Eis的Mi+1到/>
Figure BDA0001740197600000078
行，Eix、Eiy也可表达为Eix＝Ai1Ti，Eiy＝Ai1ΦiyTi。矩阵/>
Figure BDA0001740197600000079
满秩。进而有
Figure BDA00017401976000000710
其中，
Figure BDA00017401976000000711
则Φiy的对角线的元素就是Ψi的特征值。由最小二乘法则可以得到Ψi的估计
Figure BDA0001740197600000081
将
Figure BDA0001740197600000082
特征值分解可以获得/>
Figure BDA0001740197600000083
然后根据/>
Figure BDA0001740197600000084
的特征矢量，就能得到T的估计结果/>
Figure BDA0001740197600000085
在不考虑噪声的条件下有
Figure BDA0001740197600000086
Πi表示一个置换矩阵。因为
Figure BDA0001740197600000087
和/>
Figure BDA0001740197600000088
的特征值是一致的，所以将/>
Figure BDA0001740197600000089
特征值分解就可以得到/>
Figure BDA00017401976000000810
Figure BDA00017401976000000811
的估计为
Figure BDA00017401976000000812
其中
Figure BDA00017401976000000813
为/>
Figure BDA00017401976000000814
的第k个特征值。
同样对矩阵Eis重构，可得
Figure BDA00017401976000000815
类似地可以将E'is分解成矩阵E'ix＝E'is(1:Mi(Mi-1),:)和矩阵/>
Figure BDA00017401976000000816
构造矩阵
Figure BDA00017401976000000817
那么，/>
Figure BDA00017401976000000818
在不考虑噪声的条件下有/>
Figure BDA00017401976000000819
其中Πi表示一个置换矩阵。类似于uik，可以得到/>
Figure BDA00017401976000000820
的估计值
Figure BDA00017401976000000821
其中
Figure BDA00017401976000000822
表示的是矩阵(E'ix)+E'iy第k个特征值。
由于互质面阵中阵元间距大于半波长，因此会在非满秩面阵中会产生角度模糊。假设有空间角度为
Figure BDA00017401976000000823
的非相干信源入射到互质面阵，考虑面阵中的子面阵i，子面阵i在X轴和Y轴方向上相邻阵元的相位差与阵元间距之间的关系可分别表示为
Figure BDA00017401976000000824
Figure BDA00017401976000000825
其中kx和ky都为整数，取值范围分别为
Figure BDA00017401976000000826
Figure BDA00017401976000000827
同时kx和ky还需满足/>
Figure BDA00017401976000000828
当/>
Figure BDA00017401976000000829
和/>
Figure BDA00017401976000000830
确定时，存在一个或多个空间角度满足上式，当d＝Mλ/2(M为大于1的整数)时，kx和ky的取值分别为M和M/2个，即使考虑到角度配对问题，仍有多个角度估计值满足上式，这些角度中只有一个为真实角度，其他均为模糊值。大量文献已证明，互质面阵的角度模糊消除可以通过联合两个子阵间的角度关系解决，两个子阵相同的估计值即为真实角度，其余值为模糊值。
根据uik和vik的表达式可知，uik和vik的估计值有相同的列模糊，因此估计角度可以自动配对。角度配对完成后，根据以下公式即可得到信源仰角和方位角的估计值
Figure BDA0001740197600000091
Figure BDA0001740197600000092
下面利用MATLAB仿真对本发明的算法性能进行分析。其中，采用求根均方误差(Root Mean Square Error，RMSE)来评估算法DOA估计性能，RMSE定义如下
Figure BDA0001740197600000093
其中,K表示空间中的信源个数，L表示Monte Carlo试验次数。
Figure BDA0001740197600000094
和/>
Figure BDA0001740197600000095
分别表示第q次Monte Carlo试验时第k个信源仰角θk和方位角/>
Figure BDA0001740197600000096
的估计值，θk,q和/>
Figure BDA0001740197600000097
分别表示其精确值。
N和M分别表示互质面阵中的两个子阵在X轴和Y轴方向的阵元数，阵列阵元总数为M2+N2-1。假设信源是以直线或曲线运动，方位角和仰角在时间T内线性变化，且整个跟踪过程中信源数保持不变。跟踪时间T设为60s。每隔1s跟踪一次，1s内快拍数J设为200。遗忘因子β设为0.97。PASTd算法的初始参数设置如下：λk＝1(k＝1,…,K)，W(0)＝[IK,0]T，IK为K×K的单位阵。
图3-5是SNR为20dB时基于PASTd算法的DOA跟踪结果，阵列参数M＝4，N＝5。从图3-5可以看出该算法能够有效的对信源的方位角和仰角进行跟踪。图5显示了两个信源的方位角和仰角在一个图中的变化轨迹，可看出本发明能够较精确的进行DOA跟踪。
图6是本发明在不同阵元数下的DOA跟踪性能比较图。在快拍数J为300的情况下，固定参数M为5，N＝[3,4,6]，可以看出随着阵元数的增加，本发明的跟踪性能变得越来越好。
图7显示了本发明与PAST算法在不同信噪比下的跟踪性能对比图，PAST算法同样基于图1的互质面阵模型。PASTd算法随着信噪比提高性能逐步改善，而在信噪比升高后性能增长减缓，同时可以看出本发明性能优于PAST算法。本发明在高信噪比时仍可以较精确的实现信源二维DOA跟踪。
本发明的复杂度分析：在图1所示的阵列模型下，PASTd算法每次更新的复杂度为O(4(M2+N2)K+2K)，在互质面阵中得到信号子空间的情况下利用2D-ESPRIT算法进行二维DOA估计的复杂度为O(4K2M(M-1)+4K2N(N-1)+16K3)。由此可见，本发明的复杂度较低，利于实时跟踪DOA信息。
符号说明：小写(大写)粗体字来表示向量(矩阵)。(·)T、(·)H分别表示矩阵或向量的转置、共轭转置。E(·)是统计期望。
Figure BDA0001740197600000101
和/>
Figure BDA0001740197600000102
分别表示Kronecker积和Khatri-Rao积。diag(·)代表使用向量的元素作为对角元素的对角矩阵。angle(·)表示取相角。triu(·)表示取矩阵上三角元素。/>
Figure BDA0001740197600000103
表示对值x的估计。

Claims (3)

1.互质面阵下的二维DOA跟踪方法，其特征在于，包括如下步骤：
(1)通过联合互质面阵中的两个子面阵构建阵列接收信号的数据模型；
(2)利用PASTd方法实时更新接收数据的信号子空间；
(3)基于阵列接收数据的信号子空间，利用互质面阵下的二维ESPRIT算法得到DOA估计，并通过联合阵列子阵间的角度关系实现解模糊；具体为：
在互质面阵中，对于子面阵i，按下式构造Ai1，Ai2
Figure FDA0004183919680000011
Figure FDA0004183919680000012
其中，Mi＝M或Mi＝N，Ai2＝Ai1Φiy，旋转因子
Figure FDA0004183919680000013
将信号子空间Eis分解成Eix＝Eis(1:Mi(Mi-1),:)，
Figure FDA0004183919680000014
其中Eix为Eis的1到Mi(Mi-1)行，Eiy为Eis的Mi+1到/>
Figure FDA0004183919680000015
行，Eix、Eiy表达为Eix＝Ai1Ti，Eiy＝Ai1ΦiyTi，矩阵/>
Figure FDA0004183919680000016
满秩，进而有
Eiy＝EixTi -1ΦiyTi＝EixΨi
其中，Ψi＝Ti -1ΦiyTi，则Φiy的对角线的元素就是Ψi的特征值，由最小二乘法则得到Ψi的估计
Figure FDA0004183919680000017
将
Figure FDA0004183919680000018
特征值分解获得/>
Figure FDA0004183919680000019
然后根据/>
Figure FDA00041839196800000110
的特征矢量，就能得到T的估计结果/>
Figure FDA00041839196800000111
在不考虑噪声的条件下有
Figure FDA00041839196800000112
Πi表示一个置换矩阵；因为
Figure FDA00041839196800000113
和/>
Figure FDA00041839196800000114
的特征值是一致的，所以将/>
Figure FDA00041839196800000115
特征值分解就得到
Figure FDA00041839196800000116
的估计为
Figure FDA00041839196800000117
其中
Figure FDA00041839196800000118
为/>
Figure FDA00041839196800000119
的第k个特征值，θk和/>
Figure FDA00041839196800000120
分别为第k个信号的俯仰角和方位角，k＝1,2,…,K，K为信源数；
同样对矩阵Eis重构，得
Figure FDA0004183919680000021
类似地将E'is分解成矩阵E'ix＝E'is(1:Mi(Mi-1),:)和矩阵/>
Figure FDA0004183919680000022
构造矩阵
Figure FDA0004183919680000023
那么，/>
Figure FDA0004183919680000024
在不考虑噪声的条件下有/>
Figure FDA0004183919680000025
其中Πi表示一个置换矩阵，类似于uik，得到/>
Figure FDA0004183919680000026
的估计值
Figure FDA0004183919680000027
其中
Figure FDA0004183919680000028
表示的是矩阵(E'ix)+E'iy第k个特征值；
通过联合两个子阵间的角度关系，消除互质面阵中角度估计的模糊值，即两个子阵相同的估计值即为真实角度，其余值为模糊值；
根据uik和vik的表达式可知，uik和vik的估计值有相同的列模糊，因此估计角度自动配对，角度配对完成后，根据以下公式得到信源仰角和方位角的估计值
Figure FDA0004183919680000029
Figure FDA00041839196800000210
2.如权利要求1所述的互质面阵下的二维DOA跟踪方法，其特征在于，步骤(1)具体为：
阵元总数为M2+N2-1的互质面阵，由阵元数分别为N×N和M×M的两个均匀子面阵构成，且两个子阵的阵元间距分别为Mλ/2和Nλ/2；对于子面阵i，其在X轴和Y轴上阵元的方向矩阵分别表示为
Figure FDA00041839196800000211
Figure FDA00041839196800000212
其中，方向向量/>
Figure FDA00041839196800000213
Figure FDA00041839196800000214
di为阵元间距，因此，子面阵i的接收信号xi(t)表示为
Figure FDA00041839196800000215
其中，s(t)为信源矢量，ni(t)为子面阵i的噪声，
Figure FDA0004183919680000031
和/>
Figure FDA0004183919680000036
分别表示Kronecker积和Khatri-Rao积，Aix为子面阵i在X轴上阵元的方向矩阵，对于子面阵1而言，Mi就等于其方向轴上的阵元个数N，而对于子面阵2，Mi就等于方向轴上的阵元个数M。
3.如权利要求2所述的互质面阵下的二维DOA跟踪方法，其特征在于，步骤(2)具体为：
(21)选择适当的初始值λk(0)，W(0)；
(22)对每一个t＝1,2,…,J，J为快拍数，使得x1(t)＝X(t)；
(23)对每一个k＝1,2,...,K，K为信源数，分别更新以下变量：阵列接收数
Figure FDA0004183919680000032
特征值/>
Figure FDA0004183919680000033
特征向量/>
Figure FDA0004183919680000034
Figure FDA0004183919680000035
以及阵列接收数据xk+1(t)＝xk(t)-Wk(t)yk(t)；
(24)当步骤(23)中的k＝K后，使得t＝t+1，再次从步骤(22)开始计算。
CN201810815417.1A 2018-07-24 2018-07-24 互质面阵下的二维doa跟踪方法 Active CN109239651B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810815417.1A CN109239651B (zh) 2018-07-24 2018-07-24 互质面阵下的二维doa跟踪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810815417.1A CN109239651B (zh) 2018-07-24 2018-07-24 互质面阵下的二维doa跟踪方法

Publications (2)

Publication Number Publication Date
CN109239651A CN109239651A (zh) 2019-01-18
CN109239651B true CN109239651B (zh) 2023-06-20

Family

ID=65072985

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810815417.1A Active CN109239651B (zh) 2018-07-24 2018-07-24 互质面阵下的二维doa跟踪方法

Country Status (1)

Country Link
CN (1) CN109239651B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110673085B (zh) * 2019-09-25 2023-06-13 南京航空航天大学 一种均匀面阵下基于快速收敛平行因子的相干信源测向方法
CN110736959B (zh) * 2019-10-25 2021-07-09 北京理工大学 一种基于和差协同阵构建的平面互质阵列设计方法
CN111239679B (zh) * 2020-02-12 2022-04-08 南京航空航天大学 一种用于互质面阵下相干信源doa估计的方法
CN111580040A (zh) * 2020-03-29 2020-08-25 重庆邮电大学 双基地展开互质阵列mimo雷达dod和doa降维估计方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106785486A (zh) * 2017-01-09 2017-05-31 南京航空航天大学 一种广义互质面阵天线结构及角度估计方法
CN107102292A (zh) * 2017-06-19 2017-08-29 哈尔滨工业大学 一种基于贝叶斯方法的目标方位跟踪方法
CN108120967A (zh) * 2017-11-30 2018-06-05 山东农业大学 一种平面阵列doa估计方法及设备

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106785486A (zh) * 2017-01-09 2017-05-31 南京航空航天大学 一种广义互质面阵天线结构及角度估计方法
CN107102292A (zh) * 2017-06-19 2017-08-29 哈尔滨工业大学 一种基于贝叶斯方法的目标方位跟踪方法
CN108120967A (zh) * 2017-11-30 2018-06-05 山东农业大学 一种平面阵列doa估计方法及设备

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
DOD and DOA Tracking Algorithm for Bistatic MIMO Radar Using PASTd without Additional Angles Pairing;Hailang Wu 等;《2012 IEEE fifth International Conference on Advanced Computational Intelligence(ICACI)》;20130218;全文 *
互质阵中空间谱估计研究进展;张小飞 等;《南京航空航天大学学报》;20171031;第49卷(第5期);全文 *
基于PASTd的圆阵ESPRIT算法;梁炎夏 等;《系统仿真技术》;20100731;第6卷(第3期);全文 *
改进型PASTd双基地MIMO雷达相干目标角度跟踪算法;张正言 等;《火力与指挥控制》;20160430;第41卷(第4期);全文 *
运动目标DOA自适应跟踪算法;赵汇强 等;《火力与指挥控制》;20120630;第37卷(第6期);全文 *

Also Published As

Publication number Publication date
CN109239651A (zh) 2019-01-18

Similar Documents

Publication Publication Date Title
CN109239651B (zh) 互质面阵下的二维doa跟踪方法
CN106772226B (zh) 基于压缩感知时间调制阵列的doa估计方法
CN105403856B (zh) 基于嵌套式最小冗余阵列的波达方向估计方法
CN110208735B (zh) 一种基于稀疏贝叶斯学习的相干信号doa估计方法
CN105824002B (zh) 基于嵌套式子阵阵列的波达方向估计方法
CN110113085B (zh) 一种基于协方差矩阵重构的波束形成方法及系统
CN107436421B (zh) 一种稀疏贝叶斯学习框架下混合信号doa估计方法
CN108344967B (zh) 基于互质面阵的二维波达方向快速估计方法
CN109116293B (zh) 一种基于离格稀疏贝叶斯的波达方向估计方法
CN107544051A (zh) 嵌套阵列基于k‑r子空间的波达方向估计方法
CN112698264A (zh) 增广互质阵列脉冲噪声环境下相干信源的doa估计方法
CN108872926A (zh) 一种基于凸优化的幅相误差校正及doa估计方法
Moses et al. An autoregressive formulation for SAR backprojection imaging
CN112731278B (zh) 一种部分极化信号的角度与极化参数欠定联合估计方法
CN108761380B (zh) 一种用于提高精度的目标波达方向估计方法
CN112099007A (zh) 适用于非理想天线方向图的方位向多通道sar模糊抑制方法
Qi et al. Time-frequency DOA estimation of chirp signals based on multi-subarray
CN113567913B (zh) 基于迭代重加权可降维的二维平面doa估计方法
CN115236584A (zh) 基于深度学习的米波雷达低仰角估计方法
CN108614235B (zh) 一种多鸽群信息交互的单快拍测向方法
CN112763972B (zh) 基于稀疏表示的双平行线阵二维doa估计方法及计算设备
CN113189538A (zh) 一种基于互质稀疏排列的三元阵列及其空间谱估计方法
CN110579737B (zh) 一种杂波环境中基于稀疏阵列的mimo雷达宽带doa计算方法
CN116699511A (zh) 一种多频点信号波达方向估计方法、系统、设备及介质
CN106844886B (zh) 基于主分量分析的目标波达方向获取方法

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
GR01 Patent grant
GR01 Patent grant