CN103135091A - 一种波达方向估计系统中的自适应冲击噪声消除方法 - Google Patents

Abstract

本发明涉及阵列信号处理领域。所要解决的技术问题是提供一种波达方向估计系统中的自适应冲击噪声消除方法，能利用修正经验模式分解及相应手段，无需任何先验信息，对含有冲击噪声的信号进行自适应的噪声消除。本发明对阵列天线各阵元获得的信号采用修正的经验模式分解进行处理，得到多个信号成分；对各个信号成分的统计特性进行判断，找出其中某一个或多个具有冲击噪声特性的成分；从原信号中剔除冲击噪声成分，将经过噪声抑制的信号输入波达方向估计模块进行处理。本发明对阵列天线的各个阵元所获得的信号进行噪声抑制预处理，使得波达方向估计技术可以在背景噪声为冲击噪声的环境中进行使用并保持较高的估计精度和稳定性。

Description

一种波达方向估计系统中的自适应冲击噪声消除方法

技术领域

本发明涉及阵列信号处理领域。 

技术背景

波达方向估计是阵列信号处理中非常重要的研究方向，在雷达、声纳以及通信领域有着广泛的应用。现有的波达方向估计通常假设随机信号和噪声是服从高斯分布的，基于该假设，包括MUSIC和ESPRIT等经典波达方向估计方法均利用阵列接收数据的二阶以上的统计量获得理想的估计结果。大量的试验数据表明，实际应用所遇到的随机信号和噪声并不都是服从高斯分布的，很多情况下是具有拖尾概率密度函数的冲击噪声。由于用于表示冲击噪声的对称α稳定分布不存在二阶及以上矩，基于高斯噪声假设的波达方向估计方法性能出现严重退化。为解决这一问题，包括FLOM-MUSIC和ROC-MUSIC在内的在冲击噪声背景下具有高稳健性的波达方向估计方法相继被提出。但是这类方法需要冲击噪声分布特性的先验信息或是估计值，这在实际应用中难以做到；而且这类方法均基于MUSIC算法，需要进行时间代价较大的谱峰搜索。能够适用于处理效率较高的ESPRIT方法在冲击噪声背景下进行稳健波达方向估计的方法暂未提出。 

发明内容

本发明所要解决的技术问题是为克服现有技术的不足，提供一种波达方向估计系统中的自适应冲击噪声消除方法，能利用修正经验模式分解及相应手段，无需任何先验信息，对含有冲击噪声的信号进行自适应的噪声消除。 

为解决上述技术问题，本发明的技术解决方案是： 

一种波达方向估计系统中的自适应冲击噪声消除方法，其特征在于： 

所述方法，包括以下步骤： 

步骤一：对阵列天线各阵元获得的信号采用修正的经验模式分解进行处理，得到多个信号成分； 

设第i个阵元接收到的信号为s(t)具体操作如下： 

步骤1.1找到信号s(t)的局域极值点； 

步骤1.2计算s(t)的局域均值 

$\stackrel{\‾}{s}\left(t\right)=\frac{{\∫}_{t-\frac{T}{2}}^{t+\frac{T}{2}}\mathrm{\ω}(\mathrm{\τ}-t)x\left(t\right)\mathrm{d\τ}}{{\∫}_{t-\frac{T}{2}}^{t+\frac{T}{2}}\mathrm{\ω}(\mathrm{\τ}-t)\mathrm{d\τ}}---\left(1\right)$

其中ω(t)为窗函数，设为 

ω(t)＝0.42-0.5cos(2πt/T)+0.08cos(4πt/T)，0≤t≤T    (2) 

T为窗函数的宽度，由s(t)局域极值的间距确定，设极值点分别出现在时刻(t₁，t₂，…t_N)，取 $T=\frac{3}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(t_{i+1}-t_i);$

步骤1.3从信号s(t)中减去局域均值 

步骤1.4将步骤1.3得到的h(t)作为s(t)重复代入步骤1.1～步骤1.3进行计算；直到计算所得的h(t)符合以下标准： 

(1)整个数据段里，极值点的数量(包括极大值和极小值)和过零点的数量最多相差1； 

(2)局域均值 

约为0； 

称h(t)为一阶内秉模式，记为IMF₁； 

步骤1.5将得到的IMF₁从原信号s(t)减去就可以得到残余量r(t)：r(t)＝s(t)-h(t)； 

步骤1.6判断r(t)是否为趋势项；设r(t)的上、下包络曲线分别为e_max和e_min，通过e_max和e_min可以获得参考曲线 

$\mathrm{\δ}\left(t\right)=\frac{|e_{\max }+e_{\min }|}{e_{\max }-e_{\min }}---\left(3\right)$

如果符合以下条件： 

(1)δ(t)＜θ₁的时刻的个数与整个信号时间长度之比不小于1-α；θ₁＝0.05，α＝0.05； 

(2)δ(t)＜θ₂；θ₂＝10θ₁； 

则认为r(t)为趋势项，分解操作结束； 

如果残余量r(t)并不符合上述条件，则将其视作一个新的信号并将其做如下替换，r(t)→s(t)，并重复以上的各个步骤来提取出其他的各阶IMF_l，(i＝1，2，…，N)； 

最终，原信号s(t)被分解为： 

$s\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{IMF}}_i+r\left(t\right)---\left(4\right)$

其中r(t)为残余量； 

步骤二：对各个信号成分的统计特性进行判断，找出其中某一个或多个具有冲击噪声特性的成分； 

对分解得到的各个子成分即各阶IMF是否具有冲击噪声特性，通过该子成分的信息熵H和峭度值K来进行判断： 

设一时间序列x_i，i＝1，2，…，N的概率分布p＝p_i，0≤p_i≤1，i＝1，2，…，N，其信息熵H可表示为： 

$H=-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{p}_i\ln p_i---\left(5\right)$

作为噪声成分的信号应满足H＞3； 

信号的峭度值K可以通过其四阶累积量和二阶累积量进行估计： 

$K=\frac{E{[s\left(t\right)-\mathrm{\μ}]}^4}{{\mathrm{\σ}}^4}---\left(6\right)$

其中μ和σ分别为s(t)的均值和标准差，E[□]表示期望计算；冲击噪声的峭度值应满足K＞10； 

如果对以一个IMF上述两类判据均达到要求，则认为该IMF具有冲击噪声特性； 

步骤三：从原信号s(t)减去冲击噪声成分； 

通过步骤二可以对原信号s(t)分解得到的N个IMF的性质进行判断，并将具有冲击噪声的成分去除，得到不含冲击噪声的信号s′(t)， 

$s^{\′}\left(t\right)=s\left(t\right)-\underset{i\∈M}{\mathrm{\Σ}}{\mathrm{IMF}}_i---\left(7\right)$

其中M为表现冲击噪声特性的IMF的阶数的集合。 

本发明可带来以下有益效果： 

本发明对阵列天线的各个阵元所获得的信号进行噪声抑制预处理，使得波达方向估计技术可以在背景噪声为冲击噪声的环境中进行使用并保持较高的估计精度和稳定性。该方法不仅不需要对背景冲击噪声的分布参数进行先验性的估计或分析，而且对后端的波达方向估计操作没有任何技术性限制，可以与MUSIC和ESPRIT等各类波达方向估计技术相结合。因此本发明提出的冲击噪声抑制技术，使基于ESPRIT方法的高效波达方向估计技术可以在冲击噪声的环境中使用，保证了波达方向估计该技术的实时运行能力。 

附图说明

图1自适应冲击噪声消除方法原理图 

图2修正经验模式分解流程图 

图3采用降噪处理前、后DOA估计结果比较图 

图4不同广义信噪比条件下采用降噪处理前、后DOA估计误差比较图 

具体实施方式

下面结合附图和实施例对本发明作进一步详细说明。 

如图1所示： 

本发明具体实现方案如下： 

步骤一：对阵列天线各阵元获得的信号采用修正的经验模式分解进行处理，得到多个信号成分； 

步骤二：对各个信号成分的统计特性进行判断，找出其中某一个或多个具有冲击噪声特性的成分； 

步骤三：从原信号中剔除冲击噪声成分，将经过噪声抑制的信号输入波达方向估计模块进行处理。 

具体实施如下： 

步骤一：应用修正经验模式分解对阵列天线各个阵元接收到的信号进行分解处理，如图2所示。 

设第i个阵元接收到的信号为s(t)，具体操作如下： 

1找到信号s(t)的局域极值点； 

2计算s(t)的局域均值 

$\stackrel{\‾}{s}\left(t\right)=\frac{{\∫}_{t-\frac{T}{2}}^{t+\frac{T}{2}}\mathrm{\ω}(\mathrm{\τ}-t)x\left(t\right)\mathrm{d\τ}}{{\∫}_{t-\frac{T}{2}}^{t+\frac{T}{2}}\mathrm{\ω}(\mathrm{\τ}-t)\mathrm{d\τ}}---\left(1\right)$

其中ω(t)为窗函数，设为 

ω(t)＝0.42-0.5cos(2πt/T)+0.08cos(4πt/T)，0≤t≤T    (2) 

T为窗函数的宽度，由 

局域极值的间距确定，设极值点分别出现在时刻(t₁，t₂，…t_N)，取 $T=\frac{3}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(t_{i+1}-t_i).$

3从信号s(t)中减去局域均值 

4将得到的h(t)作为s(t)重复代入上述步骤进行计算，直到计算所得的h(t)符合以下标准： 

(1)整个数据段里，极值点的数量(包括极大值和极小值)和过零点的数量最多相差1； 

(2)局域均值 

约为0。 

称h(t)为一阶内秉模式(Intrinsic Mode Function，IMF)，记为IMF₁。 

5将得到的IMF₁从原信号s(t)减去就可以得到残余量r(t)：r(t)＝s(t)-h(t)； 

6判断r(t)是否为趋势项。设r(t)的上、下包络曲线分别为e_max和e_min，通过e_max和e_min可以获得参考曲线 

$\mathrm{\δ}\left(t\right)=\frac{|e_{\max }+e_{\min }|}{e_{\max }-e_{\min }}---\left(3\right)$

如果符合以下条件： 

(1)δ(t)＜θ₁的时刻的个数与整个信号时间长度之比不小于1-α(实际操作建议取θ₁＝0.05，α＝0.05)； 

(2)δ(t)＜θ₂(实际操作建议取θ₂＝10θ₁)。 

则认为r(t)为趋势项，分解操作结束。 

如果残余量r(t)并不符合上述条件，则将其视作一个新的信号并将其做如下替换，r(t)→s(t)，并重复以上的各个步骤来提取出其他的各阶IMF_i，(i＝1，2，…，N)。 

最终，原信号s(t)被分解为： 

$s\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{IMF}}_i+r\left(t\right)---\left(4\right)$

其中r(t)为残余量。 

步骤二：对各个信号成分的统计特性进行判断。由步骤一可以获得通过各个阵元接收到的信号s(t)分解得到的一系列IMF，本步骤的目的是判断出其中某个或哪些IMF是原信号中的冲击噪声成分。 

判据1：信息熵 

信息熵可以作为信号不确定性的测度，设一时间序列x_i，i＝1，2，…，N的概率分布p＝p_i，0≤p_i≤1，i＝1，2，…，N，其信息熵H可表示为： 

$H=-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{p}_i\ln p_i---\left(5\right)$

信号越是“随机”，越是具有噪声特性，其熵就越大；反之其熵越小。因此可以利用信息熵来判断一个时间序列是噪声成分还是规则信号。实际操作中，要求作为噪声成分的信号满足H＞3。 

判据2：峭度 

冲击噪声一般符合对称α稳定分布模型，其概率密度函数较普通噪声(如高斯噪声)具有明显的尖峰，且具有较大的拖尾。通过峭度值可以有效的将具有此类概率密度分布的信号检测出来。信号的峭度值K可以通过其四阶累积量和二阶累积量进行估计： 

$K=\frac{E{[s\left(t\right)-\mathrm{\μ}]}^4}{{\mathrm{\σ}}^4}---\left(6\right)$

其中μ和σ分别为s(t)的均值和标准差，E[·]表示期望计算。高斯噪声的峭度约为3，冲击噪声的峭度值应远大于该值，实际操作中选取判断标准K＞10。 

如果对以一个IMF上述两类判据均达到要求，则认为该IMF具有冲击噪声特性。 

步骤三：从原信号s(t)减去冲击噪声成分。 

通过步骤二可以对原信号s(t)分解得到的N个IMF的性质进行判断，并将具有冲击噪声的成分去除，得到不含冲击噪声的信号s′(t)， 

$s^{\′}\left(t\right)=s\left(t\right)-\underset{i\∈M}{\mathrm{\Σ}}{\mathrm{IMF}}_i---\left(7\right)$

其中M为表现冲击噪声特性的IMF的阶数的集合。 

经过实测，经该方法处理原信号中的冲击噪声成分能够得到有效的抑制，波达方向估计结果具有很高的精度和稳定度。 

如图3所示，采用降噪处理前(星号所示)、后(加号所示)DOA估计结果比较，所估计的两个信号的方位信息(方位角α和俯仰角θ)分别为(40°10°)和(60°30°)，可以看出，采用降噪处理的估计结果较为准确而未采用降噪处理的估计结果误差非常大。 

如图4所示，不同广义信噪比条件下采用降噪处理前(由菱形表示)、后(由圆圈表示)DOA估计误差比较，可以看出经过降噪处理后估计误差得到了有效降低，尤其在低广义信噪比的情况下(GSNR＜-10dB)可有效降低估计误差。 

Claims (1)

1.一种波达方向估计系统中的自适应冲击噪声消除方法，其特征在于：

所述方法，包括以下步骤：

步骤一：对阵列天线各阵元获得的信号采用修正的经验模式分解进行处理，得到多个信号成分；

设第i个阵元接收到的信号为s(t)，具体操作如下：

步骤1.1找到信号s(t)的局域极值点；

步骤1.2计算s(t)的局域均值

$\stackrel{\‾}{s}\left(t\right)=\frac{{\∫}_{t-\frac{T}{2}}^{t+\frac{T}{2}}\mathrm{\ω}(\mathrm{\τ}-t)x\left(t\right)\mathrm{d\τ}}{{\∫}_{t-\frac{T}{2}}^{t+\frac{T}{2}}\mathrm{\ω}(\mathrm{\τ}-t)\mathrm{d\τ}}---\left(1\right)$

其中ω(t)为窗函数，设为

ω(t)＝0.42-0.5cos(2πt/T)+0.08cos(4πt/T)，0≤t≤T    (2)

T为窗函数的宽度，由s(t)局域极值的间距确定，设极值点分别出现在时刻(t₁，t₂，…t_N)，取 $T=\frac{3}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(t_{i+1}-t_i);$

步骤1.3从信号s(t)中减去局域均值

步骤1.4将步骤1.3得到的h(t)作为s(t)重复代入步骤1.1～步骤1.3进行计算；直到计算所得的h(t)符合以下标准：

(1)整个数据段里，极值点的数量(包括极大值和极小值)和过零点的数量最多相差1；

(2)局域均值

约为0；

称h(t)为一阶内秉模式，记为IMF₁；

步骤1.5将得到的IMF₁从原信号s(t)减去就可以得到残余量r(t)：r(t)＝s(t)-h(t)；

步骤1.6判断r(t)是否为趋势项；设r(t)的上、下包络曲线分别为e_max和e_min，通过e_max和e_min可以获得参考曲线

$\mathrm{\δ}\left(t\right)=\frac{|e_{\max }+e_{\min }|}{e_{\max }-e_{\min }}---\left(3\right)$

如果符合以下条件：

(1)δ(t)＜θ₁的时刻的个数与整个信号时间长度之比不小于1-α；θ₁＝0.05，α＝0.05；

(2)δ(t)＜θ₂；θ₂＝10θ₁；

则认为r(t)为趋势项，分解操作结束；

如果残余量r(t)并不符合上述条件，则将其视作一个新的信号并将其做如下替换，r(t)→s(t)，并重复以上的各个步骤来提取出其他的各阶IMF_i，(i＝1，2，…，N)；

最终，原信号s(t)被分解为：

$s\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{IMF}}_i+r\left(t\right)---\left(4\right)$

其中r(t)为残余量；

步骤二：对各个信号成分的统计特性进行判断，找出其中某一个或多个具有冲击噪声特性的成分；

对分解得到的各个子成分即各阶IMF是否具有冲击噪声特性，通过该子成分的信息熵H和峭度值K来进行判断：

设一时间序列x_l，i＝1，2，…，N的概率分布p＝p_i，0≤p_i≤1，i＝1，2，…，N，其信息熵H可表示为：

$H=-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{p}_i\ln p_i---\left(5\right)$

作为噪声成分的信号应满足H＞3；

信号的峭度值K可以通过其四阶累积量和二阶累积量进行估计：

$K=\frac{E{[s\left(t\right)-\mathrm{\μ}]}^4}{{\mathrm{\σ}}^4}---\left(6\right)$

其中μ和σ分别为s(t)的均值和标准差，E[·]表示期望计算；冲击噪声的峭度值应满足K＞10；

如果对以一个IMF上述两类判据均达到要求，则认为该IMF具有冲击噪声特性；

步骤三：从原信号s(t)减去冲击噪声成分；

通过步骤二可以对原信号s(t)分解得到的N个IMF的性质进行判断，并将具有冲击噪声的成分去除，得到不含冲击噪声的信号s′(t)，

$s^{\′}\left(t\right)=s\left(t\right)-\underset{i\∈M}{\mathrm{\Σ}}{\mathrm{IMF}}_i---\left(7\right)$

其中M为表现冲击噪声特性的IMF的阶数的集合。

Images