CN114994593A - 一种基于测向设备分析信号相干关系的方法 - Google Patents

Abstract

本发明公开了一种基于测向设备分析信号相干关系的方法，包括利用阵列接收的数据构建协方差矩阵R，计算协方差矩阵R的奇异值的个数J；估计各个来波方向θ₁,θ₂,…,θ_J及其对应的相干组响应矩阵P₁,P₂,…,P_J；定义实际阵列流形矩阵：A_k+1(θ)＝[P₁，P₂，…，P_k，a(θ)]，其中，k＝1,2,…,J‑1；对θ从0°到360°进行搜索，计算任一来波方向θ_k+1的相干组响应矩阵P_k+1，得到相干组响应矩阵P₂,P₃,…,P_J；对于每个相干组，利用其相干组矩阵P_j，构建局部协方差矩阵R_j；利用R_j调用极大似然算法获得全部信号的来波方向以及该相干组内信号的个数。本发明能够明确地指出信号之间的相干关系，实现对同频信号的相干组标记，且不需改变现有测向设备的硬件。

Description

一种基于测向设备分析信号相干关系的方法

技术领域

本发明涉及无线电监测技术领域，具体的说，是一种基于测向设备分析信号相干关系的方法。

背景技术

相干组标记指的是，在利用测向算法对多个同频信号进行测向时，除给出各个信号对应的来波方向以外，还能明确指出其中哪些信号是相干信号或非相干信号。现实应用中，某些空间谱测向算法(如极大似然估计算法)仅能给出各个信号的来波方向，但无法获知各个信号之前的相干关系。相干信号通常对应于同一个物理发射源。因此，在现实的监测测向工作中，相干组标记技术可以帮助监测人员明确，哪些信号来自于同一个物理发射源。目前能够实现相干组标记的监测测向设备较少，主要的原因有：1)实际的监测设备的天线大多是均匀圆阵，不具备平移不变性，不可以直接进行空间平滑。且由于实际的监测设备一般都并非是理想阵列，采用将圆阵转化为虚拟线阵再进行空间平滑的方式，通常测向效果不佳。因此，可以直接用于对相干信号进行分辨的算法种类有限。2)若采用极大似然估计算法，理论上而言，可以估计出各方向发射源对应的协方差矩阵，然后通过其非对角线元素的复数模值来判定各个来波方向之间的相干关系。但是，由于实际监测设备一般并非理想阵列，采用这种方式进行计算，误差太大，经过实测，并无明显的相干组标记的效果。

发明内容

本发明的目的在于提供一种基于测向设备分析信号相干关系的方法，用于解决现有技术中测向设备无法有效地识别各个来波方向信号的相干关系的问题。

本发明通过下述技术方案解决上述问题：

一种基于测向设备分析信号相干关系的方法，包括：

步骤S100、利用阵列接收的数据构建协方差矩阵R，计算协方差矩阵R的奇异值的个数，作为相干组的个数J；

步骤S200、估计各个来波方向θ₁,θ₂,…,θ_J，以及计算各个来波方向对应的相干组响应矩阵P₁,P₂,…,P_J；

步骤S300、对于每个相干组，利用其相干组矩阵P_j，构建局部协方差矩阵R_j：

R_j＝P_jP_j ^H

其中，j＝1,2,…,J，利用局部协方差矩阵R_j，调用极大似然算法获得全部信号的来波方向以及该相干组内信号的个数。

所述步骤S200具体包括：

步骤S210、对来波方向θ从0°到360°进行搜索，得到第一个相干组的主信号来波方向θ₁：θ₁＝argmaxTrace{R·a(θ)(a^H(θ))^-1a^H(θ)}，其中a(θ)表示该阵列在来波方向为θ时的阵列流形矩阵；

步骤S220、利用协方差矩阵R和θ₁计算第一个相干组的主信号来波方向θ₁的相干组响应矩阵P₁；

步骤S230、定义实际阵列流形矩阵：A_k+1(θ)＝[P₁，P₂，…，P_k，a(θ)]，其中，k＝1,2,…,J-1；

对θ从0°到360°进行搜索，对于任一来波方向θ_k+1有：

θ_k+1＝argmaxTrace{R·A_k+1(θ)(A_k+1 ^H(θ)A_k+1(θ))^-1A_k+1 ^H(θ)}

对于每个k，计算来波方向θ_k+1的相干组响应矩阵P_k+1，得到相干组响应矩阵P₂,P₃,…,P_J。

步骤S220和步骤S230中相干组响应矩阵的计算方法包括：

步骤A、对协方差矩阵R进行奇异值分解R＝USU^H，得到过程参数U和过程参数S：

步骤B、利用公式Z＝US^1/2计算Z；

步骤C、计算矩阵(Z^HZ)^-1(Z^Ha(θ)a^H(θ)Z)的最大特征值对应的特征向量p_i；

步骤D、利用P＝Zp计算P_i。

在靶图上，标记出所有的信号的来波方向，用相同的颜色或样式来表示同一个相干组，进行结果展示。

本发明与现有技术相比，具有以下优点及有益效果：

(1)本发明在同频信号当中同时存在相干信号和非相干信号的情况下，能够明确地指出信号之间的相干关系，所获取的同一个相干组内的各个信号都是相干信号，而不同的相干组内的信号不是相干信号。

(2)本发明并不改变现有的空间谱测向设备的硬件，在现有设备的基础上，即可实现对同频信号的相干组标记。

附图说明

图1为本发明的流程图；

图2为三个同频信号(其中两个信号相干)情况下的实验结果，其中，(a)为协方差矩阵的奇异值从大到小排列图，(b)为实验结果靶图；

图3为三个同频信号(不相干)情况下的实验结果，其中，(a)为协方差矩阵的奇异值从大到小排列图，(b)为实验结果靶图；

图4为两个同频相干信号情况下的实验结果，其中，(a)为协方差矩阵的奇异值从大到小排列图，(b)为实验结果靶图；

图5为四个同频信号(其中两个信号相干，另外两个信号相干)情况下的实验结果，其中，(a)为协方差矩阵的奇异值从大到小排列图，(b)为实验结果靶图。

具体实施方式

下面结合实施例对本发明作进一步地详细说明，但本发明的实施方式不限于此。

实施例：

结合附图1所示，一种基于测向设备分析信号相干关系的方法，包括：

步骤S100、利用阵列接收到的数据，构建协方差矩阵R，并计算其较大的奇异值的个数，作为相干组的个数J；

步骤S200、估计各个来波方向θ₁,θ₂,…,θ_J，以及各个来波方向对应的相干组响应矩阵P₁,P₂,…,P_J；

本发明采用相干组响应矩阵P₁,P₂,…,P_J进行相干组内的测向的思路的推导如下：

阵列响应矩阵是指当信号源为1且来波方向为θ时，各个阵列接收到的信号的复数表示所构成的复数矩阵A(θ)，当信号源为s(t)，各阵列所接收的信号所构成的矩阵为A(θ)s(t)，假设存在L个相干信号s₁(t),s₂(t),…,s_K(t)，分别对应来波方向θ₁,θ₂,…,θ_J，则一定存在复数标量σ₂,σ₃,…,σ_L，使得s₂(t)＝σ₂s₁(t)，s₃(t)＝σ₃s₁(t)，…，s_L(t)＝σ_Ls₁(t)。各阵元接收的信号构成矩阵X：

X＝A(θ₁)s₁(t)+A(θ₂)s₂(t)+...+A(θ_L)s_L(t)＝(A(θ₁)+σ₂A(θ₂)+...+σ_LA(θ_L))s₁(t)

定义相干组响应矩阵P：

P＝A(θ₁)+σ₂A(θ₂)+...+σ_LA(θ_L)

相干组响应矩阵的物理意义是：在一组彼此相干、且各信号的来波方向确定的信号存在的前提下，且第1个信号定义为1，此时由各阵元所接收到的信号所构成的矩阵称为相干组响应矩阵。

则有：X＝Ps₁(t)

此时阵列的协方差矩阵R为：

R＝E{XX^H}＝(PP^H)·E{s₁(t)s₁ ^H(t)}

由于(PP^H)的秩为1，且E{s₁(t)s₁ ^H(t)}为标量，故R的秩为1。

同理，当存在J个相干组，第j个相干组有K_j个信号时，表示为s_j,1(t)，s_j,2(t)，…，

分别对应来波方向θ_j,1，θ_j,2，…

和上文同理，一定存在复数标量σ_j,2，σ_j,3，..

使得：s_j,2(t)＝σ_j,2s_j,1(t)，s_j,3(t)＝σ_j,3s_j,1(t)，…，

各阵元所接收到的信号所构成的矩阵X为：

定义第j个相干组所对应的相干组响应矩阵P_j为：

得到：

阵列的协方差矩阵为：

由于当i≠j时，E{s_i(t)s_j ^H(t)}＝0，因此有：

由矩阵理由可知：

由于[P₁ P₂ ... P_J]的秩为J，故R的秩为J。由此可以得出：协方差矩阵的秩等于相干组的个数。如果能够利用对协方差矩阵R的分解，分别计算出P₁，P₂，…，P_J，就可以分别利用各个相干组的响应矩阵进行相干组内的测向。

以下为计算相干组响应矩阵的具体方法：

步骤S210、对θ从0°到360°进行搜索，得到第一个相干组的主信号来波方向θ₁：θ₁＝argmaxTrace{R·a(θ)(a^H(θ))^-1a^H(θ)}，其中a(θ)表示该阵列在来波方向为θ时的阵列流形矩阵；

步骤S220、利用协方差矩阵R和θ₁计算第一个相干组的主信号来波方向θ₁的相干组响应矩阵P₁；

本发明采用来波方向估计相干组响应矩阵的理论推导如下：

结合上文的推导过程，令

定义W＝[λ₁P₁ λ₂P₂ ... λ_JP_J]，可知，R＝WW^H，由于满足公式的W矩阵并不唯一，不能直接通过对R的分解求得矩阵W，对于任何一个矩阵Z，只要它满足R＝ZZ^H(注：Z和W并不一定相等)，对于W任何一列w_i＝λ_iP_i，易知w_i一定位于Z所张成的空间内。故一定存在一个J×1矩阵p_i，使得P_i＝Zp_i。Z可以通过以下步骤求得:由于R是Hermite矩阵，故对R进行奇异值分解可求得U和S，使得R＝USU^H，则Z＝US^2/1。已知来波方向θ_i，可以通过以下的准则求得P_i：

约束：P_i位于Z所张成的空间中。

a(θ_i)表示该阵列在来波方向为θ_i时的阵列流形矩阵。例如，对于N元均匀圆阵有：

其中，f为频率，r为阵列半径，c为光速。幂项开头的j为虚数单位。

上式的物理意义在于，使得P_i位于信号子空间内，且与a(θ_i)的矢量夹角尽可能地最小。

将P_i＝Zp_i带入，上述的优化准则可以写为，

根据矩阵论可知，对于Hermite矩阵B，最大化问题

的解为：x是B的最大特征值所对应的特征向量。

其中C为满秩Hermite矩阵。

根据矩阵论，必存在满秩方阵E，使得C＝E^HE，则可进行如下变换：

y＝Ex A-2)

则该问题可转化为

y的最优解y_opt为矩阵E^-HBE^-1的最大特征值λ_max所对应的特征向量，故有：

E^-HBE^-1y_opt＝λ_maxy_opt

根据公式(A-2)，x的最优解x_opt满足y_opt＝Ex_opt，即x_opt＝E^-1y_opt，则有：

E^-HBx_opt＝λ_maxEx_opt，

则，E^-1E^-HBx_opt＝λ_maxx_opt，即C^-1Bx_opt＝λ_maxx_opt，故x_opt是矩阵C^-1B的特征向量，且对应特征值为λ_max。

上文的推导步骤皆可逆，故可以看出，C^-1B和E^-HBE^-1拥有相同的特征值，故λ_max也是C^-1B的最大特征值。因此，

的最优解是：x为矩阵C^-1B的最大特征值所对应的特征向量。

因此，上面公式

可与式(A-1)通过下式联系起来：

C＝Z^HZ，

D＝Z^Ha(θ_i)a^H(θ_i)Z，

故有：

C^-1D＝(Z^HZ)^-1(Z^Ha(θ_i)a^H(θ_i)Z)

故矩阵(Z^HZ)^-1(Z^Ha(θ_i)a^H(θ_i)Z)的最大特征值对应的特征向量即为所求的p_i，按照算式P_i＝Zp_i即可求得第i个相干组所对应的相干组响应矩阵P_i。

由上述推导过程可知，已知协方差矩阵R和来波方向θ_i，估计该来波方向所对应的相干组响应矩阵P_i的步骤总结如下：

第1)步，对R进行奇异值分解R＝USU^H，获得S和U；

第2)步，利用Z＝US^1/2公式计算Z。

第3)步，计算矩阵(Z^HZ)^-1(Z^Ha(θ_i)a^H(θ_i)Z)的最大特征值所对应的特征向量p_i；

第4)步，利用P_i＝Zp_i计算P_i；

因此，步骤S220中，仅需要将θ₁代入矩阵(Z^HZ)^-1(Z^Ha(θ_i)a^H(θ_i)Z)、结合R＝USU^H、Z＝US^1/2和P_i＝Zp_i则可以计算出P₁；

计算出P₁和θ₁后，在对θ₂进行估计的时候，已经将相干组响应矩阵P₁作为已知，并在对θ₂进行估计的时候，考虑了θ₁对应的相干组的影响。P₁可以看作是第一个相干组所对应的实际阵列流形。

步骤S230、定义实际阵列流形矩阵：A_k+1(θ)＝[P₁，P₂，…，P_k，a(θ)]，其中，k＝1,2,…,J-1；对θ从0°到360°进行搜索，对于任一来波方向θ_k+1有：

θ_k+1＝argmaxTrace{R·A_k+1(θ)(A_k+1 ^H(θ)A_k+1(θ))^-1A_k+1 ^H(θ)}

对于每个k，采用第1)步-第4)步的方法，计算来波方向θ_k+1的相干组响应矩阵P_k+1，得到相干组响应矩阵P₂,P₃,…,P_J；

在每个相干组内，该算法仅能估计出最强的信号。这是因为，每估计出一个来波方向以后，该算法计算出了其对应的相干组响应矩阵，参与对下一个来波方向估计的算法过程，因此，在对下一个来波方向进行估计的算法过程中，已经排除了已经估计出的所有来波方向所在的相干组的全部信息。

步骤S300、对于每个相干组，利用其相干组响应矩阵P_j，构建局部协方差矩阵R_j；

步骤S400、对于每个相干组，利用其相干组矩阵P_j，构建局部协方差矩阵R_j：

R_j＝P_jP_j ^H

其中，j＝1,2,…,J，利用局部协方差矩阵R_j，调用极大似然算法获得全部信号的来波方向以及该相干组内信号的个数。

利用局部协方差矩阵R_j，调用极大似然算法能够得到全部信号的来波方向以及该相干组内信号个数的推到过程如下：

假设有N个阵元，G个信号。G个信号对应的来波方向为θ₁，θ₂，…，θ_G，X表示各阵列接收到的信号所构成的N×1的矩阵，A为各个来波方向对应的方向矢量所构成的N×G矩阵，与θ₁，θ₂，…，θ_G有关。S为各个发射源的信号构成的G×1矩阵，n为各个阵元的噪声构成的N×1矩阵。则有：

X＝AS+n

此时对S的极大似然估计是

S＝(A^HA)^-1A^HX

拟合残差为

E||X-AS||²

＝E||(I-(A^HA)^-1A^H)X||²

＝Trace{E((I-(A^HA)^-1A^H)X)((I-(A^HA)^-1A^H)X)^H}

＝Trace{(I-(A^HA)^-1A^H)E(XX^H)(I-(A^HA)^-1A^H)^H}

＝Trace{(I-(A^HA)^-1A^H)R(I-(A^HA)^-1A^H)^H}

＝Trace{(I-A(A^HA)^-1A^H)R}

测向的优化问题是对拟合残差的最小化，即，

等效为：

在估计出θ₁，θ₂，…，θ_G以后，即可求得A，然后可以按如下方式估计出信号的幅值：

G个信号相互之间的协方差矩阵R_S为：

R_S＝ESS^H

＝(A^HA)^-1A^HE{XX^H}A(A^HA)^-1

＝(A^HA)^-1A^HRA(A^HA)^-1

R_S的对角线元素的平方根，即为估计出的各个信号源的幅值。判断信号个数，可以通过对信号源幅值的估计来进行。在未知信号源个数的情况下，可以依次增大预设的信号个数G的值，求得各个信号源的幅值，按大小排列。真实存在的信号，幅值估计较大，而不存在的信号，幅值估计明显较小，以此来判断真实信号的个数。

步骤S500、结果展示：在靶图上，标记出所有的信号的来波方向，用相同的颜色或样式来表示同一个相干组。

采用9阵元的实际测向设备同步采集各阵元的IQ数据，用于分析和验证本方案的效果。用同一台信号源所发射的信号，经过功分器以后，按照两个或多个不同的方向用天线进行发射，则这些天线发射的信号为相干信号。以下给出四组实验数据：

实验一

如图2所示，三个同频信号，其中两个信号相干，三个同频信号的参考来波方向为102°，256°，319°，其中方向为256°和319°的信号相干。其中，图2中(a)为R矩阵的奇异值的排列图。从图中可以看出，共有两个较大的奇异值，因此可以判定相干组为2。采用本发明的方法，计算出各个相干组最大的信号对应来波方向为θ₁＝101°，θ₂＝256°；其对应的相干组响应矩阵分别为：

P₁＝[0.3388-0.2263i -0.4433+0.0242i 0.3638-0.3125i 0.2486-0.2772i -0.1255-0.0810i 0.2737-0.1631i -0.0049+0.1444i 0.1202+0.2763i 0.0570+0.1715i]；

P₂＝[0.5151-0.1574i -0.0604+0.0539i 0.0785+0.1068i -0.0600+0.2828i0.2889+0.1840i 0.0137-0.2926i 0.1556+0.2397i 0.3524-0.1771i 0.1469-0.3744i]；

利用P₁计算出该相干组内信号的来波方向为101°，利用P₂计算出该相干组内信号的来波方向为256°和319°。相干组分辨的结果如图2中(b)所示。

实验二

如图3所示，三个同频信号，不相干。三个同频信号的参考来波方向为103°，258°，317°。图3中(a)为R矩阵的奇异值的排列图，从图中可以看出，共有3个较大的奇异值，因此可以判定相干组为3。采用本发明的方法，计算出各个相干组最大的信号对应来波方向为θ₁＝103°，θ₂＝259°；θ₃＝317°，其对应的相干组响应矩阵分别为：

P₁＝[0.1615+0.3570i 0.0787-0.4114i 0.1916+0.5109i 0.1088+0.2954i0.1487-0.1661i 0.0380+0.2833i -0.0161-0.1420i -0.2074+0.0048i -0.2800+0.0302i]；

P₂＝[-0.3473-0.3381i 0.1155+0.0760i 0.2034-0.1521i 0.1577+0.1621i -0.1673-0.2498i -0.2237+0.1757i -0.2073-0.3779i -0.3117-0.3123i -0.0894+0.2657i]；

P₃＝[0.3028-0.1706i -0.1892+0.2822i 0.1724-0.2321i -0.2260+0.1032i -0.1883+0.2274i 0.1509-0.2700i -0.2349+0.0502i 0.3475-0.0529i 0.4848+0.1389i]；

利用P₁计算出该相干组内信号的来波方向为103°，利用P₂计算出该相干组内信号的来波方向为259°，利用P₃计算出该相干组内信号的来波方向为317°，相干组分辨的结果如图3中(b)所示。

实验三

如图4所示，两个同频信号，相干，两个同频相干信号的参考来波方向为256°，321°。图4中(a)所示为R矩阵的奇异值的排列图。从图中可以看出，共有1个较大的奇异值，因此可以判定相干组为1。采用本发明的方法，计算出各个相干组最大的信号对应来波方向为θ₁＝256°，其对应的相干组响应矩阵分别为：

P₁＝[-0.4616-0.0000i 0.2202-0.1409i -0.0295+0.1124i 0.1827+0.0311i -0.0096-0.2797i -0.2181+0.3147i -0.0077-0.2878i -0.4413-0.1463i -0.3678+0.1078i]；

利用P₁计算出该相干组内信号的来波方向为256°，321°。相干组分辨的结果如图4中(b)所示。

实验四

如图5所示，4个同频信号，其中两个相干，另外两个也相干。4个同频信号的参考来波方向为76°，140°，256°，320°。其中76°和140°信号相干，256°和320°信号相干。图5中(a)为R矩阵的奇异值的排列图。从图中可以看出，共有两个较大的奇异值，因此可以判定相干组为2。采用本发明的方法，计算出各个相干组最大的信号对应来波方向为θ₁＝75°，θ₂＝255°；可以计算出P₁和P₂。利用P₁计算出该相干组内信号的来波方向为75°和140°。利用P₂计算出该相干组内信号的来波方向为255°和320°。相干组分辨的结果如图5中(b)所示。

由上述可以看出，本发明能够有效地识别各个来波方向信号的相干关系，能够实现对同频信号的相干组标记。

尽管这里参照本发明的解释性实施例对本发明进行了描述，上述实施例仅为本发明较佳的实施方式，本发明的实施方式并不受上述实施例的限制，应该理解，本领域技术人员可以设计出很多其他的修改和实施方式，这些修改和实施方式将落在本申请公开的原则范围和精神之内。

Claims (3)

1.一种基于测向设备分析信号相干关系的方法，其特征在于，包括：

步骤S100、利用阵列接收的数据构建协方差矩阵R，计算协方差矩阵R的奇异值的个数，作为相干组的个数J；

步骤S200、估计各个来波方向θ₁,θ₂,…,θ_J，以及计算各个来波方向对应的相干组响应矩阵P₁,P₂,…,P_J；

步骤S300、对于每个相干组，利用其相干组矩阵P_j，构建局部协方差矩阵R_j：

R_j＝P_jP_j ^H

其中，j＝1,2,…,J，利用局部协方差矩阵R_j，调用极大似然算法获得全部信号的来波方向以及该相干组内信号的个数。

2.根据权利要求1所述的一种基于测向设备分析信号相干关系的方法，其特征在于，所述步骤S200具体包括：

步骤S210、对来波方向θ从0°到360°进行搜索，得到第一个相干组的主信号来波方向θ₁：θ₁＝argmaxTrace{R·a(θ)(a^H(θ))^-1a^H(θ)}，其中a(θ)表示该阵列在来波方向为θ时的阵列流形矩阵；

步骤S220、利用协方差矩阵R和θ₁计算第一个相干组的主信号来波方向θ₁的相干组响应矩阵P₁；

步骤S230、定义实际阵列流形矩阵：A_k+1(θ)＝[P₁，P₂，…，P_k，a(θ)]，其中，k＝1,2,…,J-1；

对θ从0°到360°进行搜索，对于任一来波方向θ_k+1有：

θ_k+1＝argmaxTrace{R·A_k+1(θ)(A_k+1 ^H(θ)A_k+1(θ))^-1A_k+1 ^H(θ)}

对于每个k，计算来波方向θ_k+1的相干组响应矩阵P_k+1，得到相干组响应矩阵P₂,P₃,…,P_J。

3.根据权利要求2所述的一种基于测向设备分析信号相干关系的方法，其特征在于，步骤S220和步骤S230中相干组响应矩阵的计算方法包括：

步骤A、对协方差矩阵R进行奇异值分解R＝USU^H，得到过程参数U和过程参数S：

步骤B、利用公式Z＝US^1/2计算Z；

步骤C、计算矩阵(Z^HZ)^-1(Z^Ha(θ)a^H(θ)Z)的最大特征值对应的特征向量p_i；

步骤D、利用P＝Zp计算P_i。

Images