CN114325570A - 一种高精度噪声自适应测向方法 - Google Patents

Abstract

本发明涉及阵列信号处理领域，公开了一种高精度噪声自适应测向方法，包括建立基于均匀线阵的噪声自适应阵列测向模型，得到阵列接收数据矢量；对阵列接收数据进行去冲击平滑预处理：利用冲击噪声的分布特性进行噪声类型判定，对冲击噪声背景的阵列接收数据采用部分中位数滤波进行高斯化平滑预处理；计算阵列接收数据协方差矩阵，并通过矩阵重构平滑噪声；然后，对阵列接收数据协方差矩阵进行特征分解获得信号子空间与噪声子空间，利用子空间投影技术进一步消除信号子空间中的噪声成分；最后，采用ESPRIT测向算法进行测向。本方法同时适应冲击噪声和高斯噪声的测向，实现速度快，且在特征指数小于0.5的强冲击噪声背景下仍具有较高的测向精度。

Description

一种高精度噪声自适应测向方法

技术领域

本发明涉及阵列信号处理领域，具体涉及一种高精度噪声自适应测向方法。

背景技术

经典的DOA估计算法，大都在假定阵元噪声为空时独立高斯白噪声的前提下提出的。然而，在实际的应用环境中，有许多具有冲击特性的非高斯噪声，如水下噪声，低频大气噪声和各种人造噪声等。现有技术中针对高斯噪声和冲击噪声提出了一系列的性能优良的测向算法，但是这些算法只能适用于高斯噪声或者冲击噪声中的一种，不具备噪声自适应特性。

为解决上述问题，使测向算法能够应对复杂多变的实际噪声环境，现有技术提出了在针对任意阵列相干测向提出动态随机加权低阶协方差-空间平滑-MUSIC测向方法，其能够同时适用于高斯噪声和冲击噪声，估计性能良好。然而该方法采用谱峰搜索，计算量繁重，且在强冲击噪声背景时的测向精度低。

因此现在急需一种噪声自适应测向方法，能够同时适应冲击噪声和高斯噪声背景的测向，且在强冲击噪声背景下仍具有较高的测向精度。

发明内容

本发明意在提供一种高精度噪声自适应测向方法，用来解决现有的噪声自适应方法在强冲击噪声背景下精度低的技术问题。

本发明提供的基础方案为：一种高精度噪声自适应测向方法，包括以下步骤：

步骤一：建立基于均匀线阵的噪声自适应阵列测向模型，得到阵列接收数据矢量；

步骤二：对阵列接收数据进行去冲击平滑预处理，根据冲击噪声分布特性判定噪声背景：

若噪声背景判定为高斯噪声，则直接执行步骤三；若判定为冲击噪声，则需先采用部分中位数滤波方式对阵列接收数据进行去冲击平滑预处理，再执行步骤三；

步骤三：计算阵列接收数据的协方差矩阵并重构，再次平滑噪声成分；

步骤四：利用子空间投影技术进一步消除噪声成分并测向：

对重构后的数据协方差矩阵进行特征值分解，得到与信号相对应的信号子空间和与噪声相对应的噪声子空间；利用子空间投影技术进一步消除噪声成分，得到处理后的信号子空间；

根据ESPRIT算法，使用经子空间投影去噪后的信号子空间完成目标信号的DOA估计。

本发明的工作原理及优点在于：

本发明一种高精度噪声自适应测向方法，首先，建立基于均匀线阵的噪声自适应阵列测向模型，得到阵列接收数据矢量；其次，对阵列接收数据进行去冲击平滑预处理，利用冲击噪声的分布特性进行噪声类型判定；然后，计算阵列接收数据协方差矩阵，并通过矩阵重构再次平滑噪声；再然后，对阵列接收数据协方差矩阵进行特征分解获得信号子空间与噪声子空间，利用子空间投影技术进一步消除信号子空间中的噪声成分，提高测向算法估计精度；最后，采用经典的ESPRIT测向算法进行测向。

有益效果：现有技术大多基于分数低阶矩计算或者多维优化搜索，处理较为复杂，计算量大，且现有技术很少有能够同时适应冲击噪声和高斯噪声背景的测向；同时，现有的自适应处理方法普遍存在强冲击噪声背景下性能严重恶化甚至失效的情况。而相比于现有技术，本方法对阵列接收数据进行去冲击平滑预处理，并利用噪声的分布特性对冲击噪声进行识别，并通过矩阵重构和子空间投影技术进行去冲击平滑处理，进一步消除了噪声成分的影响，弥补了现有技术计算复杂且对强冲击噪声适应性差的缺陷，从而降低了计算复杂度，使处理过程更简便，不仅对强冲击噪声适应性优良，且在极恶劣的强冲击噪声环境仍然具有良好的测向性能。

进一步，步骤一中，基于均匀线阵的噪声自适应阵列测向模型的建立方法包括：

设均匀线阵的阵元个数为M，阵元间距为d，N个不相关的信号入射到该均匀线阵，入射角度为θ_i，i＝1，2，…，N；阵列接收噪声为高斯噪声或者冲击噪声，且信号与噪声之间相互独立；

阵列在时刻t的阵列接收信号数据可以表示为X(t)＝AS(t)+N(t)；其中，X(t)＝[x₁(t)，…，x_M(t)]^T，x_j(t)表示t时刻第j(j＝1，2，…，M)个阵元位置处的接收数据，[*]^T表示矩阵向量的转置运算；A＝[a(θ₁)，…，a(θ_N)]，为阵列流型，

表示第i个信号的阵列导向矢量，且

d为阵元间距，λ为信号波长；S(t)＝[s₁(t)，…，s_N(t)]^T，N(t)＝[n₁(t)，…，n_M(t)]^T，分别表示t时刻的信号数据矩阵和噪声数据矩阵；

N(t)为高斯噪声或冲击噪声，使用SαS分布对冲击噪声进行建模，其中α为特征指数。

有益效果：采用本方案中的模型满足ESPRIT算法的使用条件，使得最终进行测向时的计算更简单。

进一步，在步骤二中，判定噪声背景的方法包括：

步骤21：求各阵元接收数据x_j的二阶矩：γ_j＝x_jx_j ^H/L，L表示数据快拍数，则各阵元接收数据协方差表示为：γ＝[γ₁，γ₂，…，γ_M]；

步骤22：计算γ中每个元素与其均值的绝对差值并进行归一化：γ_norm＝abs[γ-mean(γ)]，其中，mean(*)表示求数组均值，abs(*)表示求绝对值；

步骤23：根据步骤22的计算结果进行判定：若存在γ_norm大于0.1的数值则判定该阵列接收数据噪声背景为冲击噪声，反之则判定为高斯噪声。

有益效果：冲击噪声的分布具有显著的非均匀特性，其只在少数位置产生很大的冲击数值；因此，在冲击噪声背景下，各阵元的接收数据协方差将存在较大的差异，特征指数α越小，噪声冲击性越强烈，各阵元数据协方差之间存在的差异越大，利用该分布特性能够清楚地识别冲击噪声和高斯噪声，便于后续对冲击噪声进行平滑预处理，将冲击噪声高斯化，既满足了噪声的自适应需求，又便于后续选择针对性的算法处理，降低后续处理难度。

进一步，在步骤二中，部分中位数滤波方式包括：

设平滑窗口长度为奇数K，分别对M个阵元接收数据矢量进行部分中位数滤波平滑处理，第j(j＝1，2，…，M)个阵元的接收数据矢量x_j在t时刻的数据进行部分中位数滤波处理后可表示为：

其中，amp(*)表示求复数的幅度运算，A_jt＝med(x_j(t))，表示以K为窗口长度对处数据进行幅度值中位数运算的结果。经上述去冲击平滑处理能够有效平滑x_j(t)中的冲击噪声成分，从而使得平滑处理后的阵列接收数据的二阶矩阵具有测向可用性。经过冲击噪声平滑预处理后的阵列接收矩阵可表示为：

有益效果：根据冲击噪声的随机分布特性和冲击特性，门限处理或中值滤波这类噪声平滑方法均可对此类噪声进行平滑处理；然而门限处理方法只能对于一个特定门限设定方案，在不同冲击强度时估计性能存在较大差异，需要频繁设置不同的门限，增大处理难度；而本方案采用部分中位数滤波方式来对阵列接收数据进行平滑处理，相比于门限处理更加简单；且相比于均值滤波等方式能够使冲击噪声高斯化处理的效果更好，进而使最终的DOA估计精度更高。

进一步，步骤三中，再次平滑噪声包括以下步骤：

步骤31：计算经过冲击噪声自适应处理所得的阵列接收数据，得到数据协方差矩阵

步骤32：根据步骤31得到的数据协方差矩阵的结构特性进行矩阵重构：

其中，

有益效果：在经过噪声背景判定后，针对不同的噪声再次进行平滑处理，使平滑处理的效果更好，从而降低噪声对测向精度的影响。

进一步，步骤四中，对重构后的阵列数据协方差矩阵

进行特征值分解的公式为：

其中，

为数据协方差矩阵，U_S为信号子空间，U_N为与噪声相对应的噪声子空间。

有益效果：提取信号子空间，满足使用子空间类测向算法条件。

进一步，步骤四中，利用子空间投影技术进一步消除噪声包括：

步骤41：构造噪声子空间的正交补空间：

步骤42：计算不含噪声成分的信号子空间的投影矩阵：

步骤43：得出剔除噪声成分的信号子空间

有益效果：进一步消除了信号子空间中的噪声，使得最终的测向精度更高。

附图说明

图1为本发明一种高精度噪声自适应测向方法实施例一的流程图；

图2为广义信噪比为10dB特征指数为0.2时本发明方法对4个独立信号的测向结果；

图3为广义信噪比为10dB特征指数为1.0时本发明方法和DSW-SS-MUSIC方法对4个独立信号的测向结果；

图4为广义信噪比为10dB特征指数为1.6时本发明方法和DSW-SS-MUSIC方法对4个独立信号的测向结果；

图5为广义信噪比为10dB高斯噪声时本发明方法和DSW-SS-MUSIC方法对4个独立信号的测向结果；

图6为广义信噪比为20dB时本发明方法在不同特征指数(噪声强度)时的估计成功概率；

图7为广义信噪比为20dB时本发明方法不同特征指数(噪声强度)时的估计均方根误差。

具体实施方式

下面通过具体实施方式进一步详细的说明：

实施例一：

实施例基本如图1所示，一种高精度噪声自适应测向方法，包括以下步骤：

步骤一：建立基于均匀线阵的噪声自适应阵列测向模型：

步骤11：设均匀线阵的阵元个数为M，阵元间距为d，N个不相关的信号入射到该均匀线阵，入射角度为θ_i，i＝1，2，…，N；阵列接收噪声为高斯噪声或者冲击噪声，且信号与噪声之间相互独立；

阵列在时刻t的阵列接收信号数据可以表示为X(t)＝AS(t)+N(t)，其中，且X(t)＝[x₁(t)，…，x_M(t)]^T，x^j(t)表示t时刻第j(j＝1，2，…，M)个阵元位置处的接收数据，[*]^T表示矩阵向量的转置运算；A＝[a(θ₁)，…，a(θ_N)]，为阵列流型，

表示第i个信号的阵列导向矢量，且

d为阵元间距，λ为信号波长；S(t)＝[s₁(t)，…，s_N(t)]^T，N(t)＝[n₁(t)，…，n_M(t)]^T，分别表示i时刻的信号数据矩阵和噪声数据矩阵；

步骤12：N(t)为高斯噪声或冲击噪声，当噪声背景为冲击噪声时，使用特征值为0～2之间的SαS分布对N(t)进行建模。

步骤二：利用噪声自适应处理模块处理阵列接收数据，根据冲击噪声分布特性判定噪声背景：

步骤21：求各阵元接收数据x_j的二阶矩，γ_j＝x_jx_j ^H/L，L表示数据快拍数，则各阵元接收数据协方差可表示为：γ＝[γ₁，γ₂，…，Y_M]；

步骤22：计算γ中每个元素与其均值的绝对差值并进行归一化：γ_norm＝abs[γ-mean(γ)]，其中，mean(*)表示求数组均值，abs(*)表示求绝对值；

步骤23：判定噪声背景，若存在γ_norm大于0.1的数值则判定该阵列接收数据噪声背景为冲击噪声，反之则判定为高斯噪声；

步骤24：针对不同噪声背景选择平滑预处理方法，若噪声背景判定为高斯噪声，则直接执行步骤三；若判定为冲击噪声，则需先采用部分中位数滤波方式对阵列接收数据进行去冲击平滑预处理，再执行步骤三；

步骤24中，部分中位数滤波方式具体内容如下：

设平滑窗口长度为奇数K，分别对M个阵元接收数据矢量进行部分中位数滤波平滑处理，第j(j＝1，2，…，M)个阵元的接收数据矢量x_j在t时刻的数据进行部分中位数滤波处理后可表示为：

其中，amp(*)表示求复数的幅度运算，A_jt＝med(x_j(t))，表示以K为窗口长度对处数据进行幅度值中位数运算的结果。经上述去冲击平滑处理能够有效平滑x_j(t)中的冲击噪声成分，从而使得平滑处理后的阵列接收数据的二阶矩阵具有测向可用性。经过冲击噪声平滑预处理后的阵列接收矩阵可表示为：

步骤三：计算阵列接收数据的协方差矩阵并重构：

步骤31：利用经过冲击噪声自适应处理所得的阵列接收数据计算数据协方差矩阵，得到：

步骤32：利用阵列数据协方差矩阵的结构特性进行矩阵重构，再次平滑噪声成分，得到：

其中，

步骤四：利用子空间投影技术进一步消除噪声成分并测向：

步骤41：对重构后的阵列数据协方差矩阵

进行特征值分解，

得到与信号相对应的信号子空间U_S和与噪声相对应的噪声子空间U_N；

步骤42：利用子空间投影技术进一步消除噪声成分包括以下步骤：

步骤421：构造噪声子空间的正交补空间：

步骤422：不含噪声成分的信号子空间的投影矩阵可表示为：

步骤423：剔除噪声成分的信号子空间

可表示为：

步骤43：根据经子空间投影去噪后的信号子空间

按照经典的ESPRIT算法原理完成目标信号的DOA估计。

由于ESPRIT算法的优点在于计算简单，实现速度快，实时响应能力更强，相比于ESPRIT算法，MUSIC算法的在估计偏差、估计方差、成功概率等方面性能几乎全面领先，即使是经典ESPRIT算法的改进算法也很难在性能和精度上超过MUSIC算法，因此现有技术在研究提高DOA估计精度的算法时通常都是选择MUSIC算法或者MUSIC算法的改进算法，而很少有人在研究精度时采用ESPRIT算法；即使有人在研究如何提高精度时选择ESPRIT算法，也通常会选择改进后的LS-RSPRIT或TLS-RSPRIT，或者其他改进算法，不会选择经典的ESPRIT算法。然而虽然MUSIC算法有以上的优势，但是MUSIC算法通常需要进行多维搜索，在计算量和实时性方面都有欠缺。而在本方法中，通过噪声类型判定和部分中位数滤波处理，使得冲击噪声背景下的阵列接收数据也能够具有有效的二阶矩阵，并结合矩阵重构再次平滑噪声，然后通过子空间投影技术进一步消除信号子空间中的噪声成分，极大地消除的噪声成分对测向精度的恶劣影响，相比于其他的自适应算法，本方法在强冲击噪声下的测向精度优势尤为突出；同时，由于本方法通过噪声背景判定后针对性的选择最佳方法，结合矩阵重构和子空间投影技术，经过两次噪声平滑后再消除，极大的降低了噪声影响，从而能够得到更高的精度，弥补了经典ESPRIT算法因阵列孔径损失带来性能下降的缺陷，因此使经典ESPRIT算法得以应用，有效扩展了经典ESPRIT测向算法的噪声适应性，进而在噪声自适应且提高DOA估计精度的基础上进一步简化了计算，加快了实现速度，能够同时兼顾精度和实时性。

下面通过图2-图7证明本方案的有效性：

其中，由于DSW-SS-MUSIC方法同样能够同时适应高斯噪声、弱冲击噪声和强冲击噪声的测向，因此选用DSW-SS-MUSIC方法在相同条件下进行对比，阵列结构等仿真建模条件参照本方法进行。本实施例中，噪声自适应测向方法的参数设置为：阵元个数M＝8，阵元间距d＝λ/2，当DOA估计结果与实际角度误差小于1°时认为估计成功。

图2为广义信噪比为10dB时，特征指数为0.2的强冲击噪声背景下对4个独立信号进行30次独立估计的测向结果，信号入射角度分别为40°、60°、80°和100°。可见，估计值和真实值之间存在一定偏差，但偏差较小，基本都能在1°之内，因此本方案对强冲击噪声有良好的测向能力。

图3-图5为广义信噪比为10dB时，特征指数分别为1.0和1.6，以及高斯噪声时本发明方法和DSW-SS-MUSIC方法对4个独立信号进行30次独立估计的向结果，信号入射角度分别为40°、60°、80°和100°。从图中可以看出，本发明方法估计值和真实值之间几乎完美重合，故所设计的测向方法对弱冲击噪声和高斯噪声具有优良的测向性能。由图3可知，在特征指数α为1时，冲击噪声较强，DSW-SS-MUSIC方法估计误差增大，估计性能存在一定程度的恶化；且随着特征指数α越来越小，即冲击噪声越来越强，DSW-SS-MUSIC方法的估计误差也越来越大；其他能够同时适应高斯噪声和冲击噪声的测向方法虽然未在此处作出对比，但普遍存在强冲击噪声背景(特征指数α＜0.5)下性能严重恶化甚至失效的情况。

图6-图7为广义信噪比为20dB时，特征指数α从0.1-2变化时，本发明的高精度噪声自适应测向方法对2个独立信号(65°，80°)进行200次独立估计后统计所得的估计成功概率和估计均方根误差。

从这两幅图也再次验证，本方法估计成功概率高，精度优良，对于强冲击噪声、弱冲击噪声和高斯噪声都具有很好的测向性能。而相比于本方法，其他方法大多数根本不能适应同时处理冲击噪声和高斯噪声，而少数能同时处理两种噪声的方法在强冲击噪声下精度也达不到这么高，即使有其他方法能达到类似的精度，也不能保证测向的实时性。

以上所述的仅是本发明的实施例，方案中公知的具体结构及特性等常识在此未作过多描述，所属领域普通技术人员知晓申请日或者优先权日之前发明所属技术领域所有的普通技术知识，能够获知该领域中所有的现有技术，并且具有应用该日期之前常规实验手段的能力，所属领域普通技术人员可以在本申请给出的启示下，结合自身能力完善并实施本方案，一些典型的公知结构或者公知方法不应当成为所属领域普通技术人员实施本申请的障碍。应当指出，对于本领域的技术人员来说，在不脱离本发明结构的前提下，还可以作出若干变形和改进，这些也应该视为本发明的保护范围，这些都不会影响本发明实施的效果和专利的实用性。本申请要求的保护范围应当以其权利要求的内容为准，说明书中的具体实施方式等记载可以用于解释权利要求的内容。

Claims (7)

1.一种高精度噪声自适应测向方法，其特征在于，包括以下步骤：

步骤一：建立基于均匀线阵的噪声自适应阵列测向模型，得到阵列接收数据矢量；

步骤二：根据冲击噪声分布特性判定噪声背景，并根据判定结果对阵列接收数据进行去冲击平滑处理：

若噪声背景判定为高斯噪声，则直接执行步骤三；若判定为冲击噪声，则需先采用部分中位数滤波方式对阵列接收数据进行去冲击平滑预处理，再执行步骤三；

步骤三：计算阵列接收数据的协方差矩阵并重构，再次平滑噪声成分；

步骤四：利用子空间投影技术进一步消除噪声成分并测向：

对重构后的数据协方差矩阵进行特征值分解，得到与信号相对应的信号子空间和与噪声相对应的噪声子空间；

利用子空间投影技术进一步消除噪声成分，得到处理后的信号子空间；

根据ESPRIT算法，使用经子空间投影去噪后的信号子空间完成目标信号的DOA估计。

2.根据权利要求1所述的一种高精度噪声自适应测向方法，其特征在于，步骤一中，基于均匀线阵的噪声自适应阵列测向模型的建立方法包括：

设均匀线阵的阵元个数为M，阵元间距为d，N个不相关的信号入射到该均匀线阵，入射角度为θ_i，i＝1，2，…，N；阵列接收噪声为高斯噪声或者冲击噪声，且信号与噪声之间相互独立；

阵列在时刻t的阵列接收信号数据可以表示为X(t)＝AS(t)+N(t)；其中，X(t)＝[x₁(t)，…，x_M(t)]^T，x_j(t)表示t时刻第j(j＝1，2，…，M)个阵元位置处的接收数据，[*]^T表示矩阵向量的转置运算；A＝[a(θ₁)，…，a(θ_N)]，为阵列流型，

表示第i个信号的阵列导向矢量，且

d为阵元间距，λ为信号波长；S(t)＝[s₁(t)，…，s_N(t)]^T，N(t)＝[n₁(t)，…，n_M(t)]^T，分别表示t时刻的信号数据矩阵和噪声数据矩阵；

N(t)为高斯噪声或冲击噪声，使用SαS分布对冲击噪声进行建模，其中α为特征指数。

3.根据权利要求1所述的一种高精度噪声自适应测向方法，其特征在于，在步骤二中，判定噪声背景的方法包括：

步骤21：求各阵元接收数据x_j的二阶矩：γ_j＝x_jx_j ^H/L，L表示数据快拍数，则各阵元接收数据协方差表示为：γ＝[γ₁，γ₂，…，γ_M]；

步骤22：计算γ中每个元素与其均值的绝对差值并进行归一化：γ_norm＝abs[γ-mean(γ)]，其中，mean(*)表示求数组均值，abs(*)表示求绝对值；

步骤23：根据步骤22的计算结果进行判定：若存在γ_norm大于0.1的数值则判定该阵列接收数据噪声背景为冲击噪声，反之则判定为高斯噪声。

4.根据权利要求1所述的一种高精度噪声自适应测向方法，其特征在于，在步骤二中，部分中位数滤波方式包括：

设平滑窗口长度为奇数K，分别对M个阵元接收数据矢量进行部分中位数滤波平滑处理，第j(j＝1，2，…，M)个阵元的接收数据矢量x_j在t时刻的数据进行部分中位数滤波处理后可表示为：

其中，amp(*)表示求复数的幅度运算，A_jt＝med(x_j(t))，表示以K为窗口长度对处数据进行幅度值中位数运算的结果；上述去冲击平滑处理能够有效平滑x_j(t)中的冲击噪声成分，从而使得平滑处理后的阵列接收数据的二阶矩阵具有测向可用性；过冲击噪声平滑预处理后的阵列接收矩阵可表示为：

5.根据权利要求1所述的一种高精度噪声自适应测向方法，其特征在于，步骤三中，再次平滑噪声包括以下步骤：

步骤31：计算经过冲击噪声自适应处理所得的阵列接收数据，得到数据协方差矩阵

步骤32：根据步骤31得到的数据协方差矩阵的结构特性进行矩阵重构：

其中，

6.根据权利要求1所述的一种高精度噪声自适应测向方法，其特征在于，步骤四中，对重构后的阵列数据协方差矩阵

进行特征值分解的公式为：

其中，

为数据协方差矩阵，U_S为信号子空间，U_N为与噪声相对应的噪声子空间。

7.根据权利要求6所述的一种高精度噪声自适应测向方法，其特征在于，步骤四中，利用子空间投影技术进一步消除噪声包括：

步骤421：构造噪声子空间的正交补空间：

步骤422：计算不含噪声成分的信号子空间的投影矩阵：

步骤423：得出剔除噪声成分的信号子空间

Images