基于方差中值滤波的极/超低频信道大气噪声抑制方法与流程

本发明涉及无线电长波通信技术领域，具体涉及一种基于方差中值滤波的极/超低频信道大气噪声抑制方法。

背景技术：

极/超低频通信(3hz～300hz)的信号衰减率小，能够穿透地层和海水到达一定的深度，因而在矿井应急通信与定位、水下无人航行器导航等领域获得了一定的应用前景。影响极/超低频通信性能的主要因素是接收端信噪比。通常，极/超低频通信的电磁干扰主要来源于信道大气噪声和本地人为干扰。信道大气噪声由自然界雷电活动引发。我国南海区域雷暴活动频繁，尤其是在夏季。在实际通信应用中，如遇本地发生雷电，大气噪声干扰表现为一簇强烈脉冲，这些脉冲会严重恶化接收信噪比，造成错报甚至漏报，使得可靠通信无法建立。极/超低频信道大气噪声呈非高斯性，传统的线性最优滤波理论难以取得理想的抑制效果。本地人为干扰是指由人类活动引发的电磁干扰，在极/超低频通信中主要表现为工频谐波干扰，这些干扰广泛存在于电力线和电力电子设备附近。上述干扰通过多种方式耦合至接收天线，给极/超低频通信带来严重影响，因此干扰抑制是极/超低频通信走向应用必须解决的一个关键问题。

针对极/超低频信道大气噪声的处理，人们已经提出了多种方法，按其思路可分为三类：1)基于噪声模型参数估计的最优接收处理，在已知大气噪声的统计模型(一阶概率密度函数)的前提下，通过比较似然比与一个固定的门限值，完成接收判决；2)非线性处理(次优接收机)，如判决器(hardlimiter)、限幅器(clipper)和凿孔器(holepuncher)，直接对接收序列的幅值进行压缩，除判决器外，其他两种非线性处理均需要依据噪声统计模型参数来确定最优触发电平；3)中值滤波处理，以滑动窗中样点的中值替换中心样点值来实现脉冲的平滑处理。上述1)和2)本质上都属于一种参数化处理方式，该类方法依赖于大气噪声的统计模型假设，然而目前获得应用的大气噪声统计模型有多种，如field-lewinstein模型(经验模型)、middleton的classa和classb模型(统计物理模型)等，每种模型均有不同的参数设置，这就为模型的选用带来困难。同时，为获取模型参数，需要完成建模、验模以及参数估计等一系列操作，运算复杂度极高。3)所述中值滤波是一种非参数化的脉冲处理方法，该技术无需脉冲噪声统计模型假设，运算复杂度低，以其良好的脉冲抑制能力和波形细节保护能力，在处理图像脉冲噪声如椒盐噪声和随机值噪声等方面获得广泛应用。然而，中值滤波存在两个缺陷：1)、脉冲抑制性能依赖于噪声的脉冲性(impulsiveness)的大小；2)、对于高密度脉冲噪声，中值滤波需要在脉冲抑制和波形细节损失两个方面进行折衷。

为了最大限度抑制带外干扰以降低噪声对信号检测的影响，工程应用中的极/超低频接收机应属于窄带接收机。在该条件下，大气噪声带宽大于接收机带宽，噪声将在接收机上产生显著的暂态响应，表现为脉冲时域波形被展宽，高频部分被滤除，脉冲性降低，接收序列呈现连续的非孤立脉冲点，局部脉冲密度显著增大，此时，常用的中值滤波技术往往难以取得令人满意的脉冲抑制效果。

技术实现要素：

针对上述现有方法在极/超低频信道大气噪声抑制方面存在的不足，本发明的目的在于提出一种基于方差中值滤波的极/超低频信道大气噪声抑制方法，以有效抑制大气噪声，提高极/超低频通信接收信噪比。

本发明的技术方案为：一种基于方差中值滤波的极/超低频信道大气噪声抑制方法，包括：

接收系统采用多组参考传感器与一个主传感器本地正交布设的方案，构建多参考磁矢量传感系统；

利用所述多参考磁矢量传感系统多通道同步接收极/超低频信号，得到磁感应输出样点矩阵x；

按照通道样点时序，对所述磁感应输出样点矩阵x进行分组，估计各分组的方差，得到方差矩阵d；

基于大数判决准则，对方差矩阵d进行脉冲检测、第一级自适应中值滤波和第二级自适应中值滤波；

根据所述第二级自适应中值滤波的输出和所述第一级自适应中值滤波的输入进行序列重建，得到输出矩阵y；

对所述输出矩阵y进行工频谐波干扰和残存大气噪声抑制。

进一步的，所述多参考磁矢量传感系统构建过程中，每组所述参考传感器包含两个正交的传感器，每组参考传感器均与所述主传感器构成三轴磁矢量传感器结构，所述主传感器放置于电屏蔽铝壳腔内，多组所述参考传感器沿所述主传感器磁芯轴线方向，等间隔均匀固定于铝壳外表面。

进一步的，所述按照通道样点时序，对所述磁感应输出样点矩阵x进行分组，估计各分组的方差，得到方差矩阵d包括：

通道u第i组的分组序列vu,i＝[xu((i-1)n+1),xu((i-1)n+2),...,xu(in)]^t；

分组矩阵

利用公式对通道u第i组的分组序列vu,i的方差进行估算，得到方差矩阵

其中，i＝1，2……i，u＝0,1，2……2u，所述i表示组数，所述u表示共有u组参考传感器，所述u表示第u个通道，i表示第i组，n表示每个分组内的样点个数，di为方差矩阵d的列向量。

进一步的，所述基于大数判决准则，对方差矩阵d进行脉冲检测包括：构造中心元素为du,i的滑动窗fu,i＝[du,i-m,...,du,i,...,du,i+m]，根据滑动窗fu,i中最大值(fu,i)max与du,i+t的大小关系判断所述方差矩阵d中各元素是否为脉冲元素；所述基于大数判决准则，仅当di中至少有u+1个元素被判决为脉冲元素时，才判定di为脉冲向量，所述m用于表征窗口长度l(l＝2m+1)，t为判决门限值。

进一步的，所述对方差矩阵d进行第一级自适应中值滤波包括：

各通道同步采用相同参数设置的自适应中值滤波算法，对于非脉冲向量不作滤波；

对于脉冲向量仅当以各元素为中心元素构建的通道滑动窗的中值mn均判定为非脉冲时，才用中值mn替换脉冲矢量作为中值滤波输出，否则，增大滑动窗长，再次进行脉冲检测。

进一步的，所述基于大数判决准则，对方差矩阵d进行第二级自适应中值滤波包括：对所述脉冲检测过程中检测到的脉冲向量及其邻近向量的第一级滤波输出结果，进行基于大数判决准则的脉冲检测和第二级自适应中值滤波，其他向量的第一级滤波输出结果保持不变，得到二级自适应中值滤波输出方差矩阵

进一步的，根据所述第二级自适应中值滤波的输出和所述第一级自适应中值滤波的输入进行序列重建，得到输出矩阵y包括：

比较第二级自适应中值滤波输出方差矩阵与第一级自适应中值滤波输入方差矩阵d的元素取值，确定比例系数

根据比例系数τu,i对通道u第i组的分组序列vu,i进行线性幅度压缩，按照分组的时序，构建幅值压缩后的通道样点矩阵即输出矩阵：

进一步的，所述输出矩阵y的工频谐波干扰采用级联二阶iir陷波滤波器进行压制，残存大气噪声采用多通道指数加权递归最小二乘(rls)算法的自适应噪声抵消技术进行压制。

本发明的有益效果：本发明对样点矩阵进行分组，对分组的方差进行自适应中值滤波，引入大数判决准则辅助判定脉冲状态，依据自适应中值滤波前后方差矩阵元素的大小变化确定比例系数，按比例压缩通道样点的幅值，从而达到抑制雷电脉冲的目的。对于工频谐波干扰和残存的大气噪声，分别采用级联二阶iir陷波滤波器和多通道指数加权递归最小二乘算法(rls)的自适应噪声抵消器进行抑制。

本发明与现有技术相比，具有以下优势：

第一、不依赖于噪声统计模型假设，操作步骤简单、方便、高效。

第二、以样点分组的方差作为自适应中值滤波器的输入，增强了脉冲性，改进了传统的中值滤波在局部高密度脉冲噪声条件下性能低下的缺陷。

第三、引入大数判决准则，提高了脉冲样点判决的准确率。

附图说明

图1为本发明涉及的信号处理流程示意图；

图2为本发明涉及的基于大数判决的自适应中值滤波算法流程示意图(不含第二级自适应中值滤波)；

图3为采用本发明涉及的多参考磁矢量传感系统，主传感器输出信号波形图；其中，大气噪声由大气噪声模拟器产生，频率范围为1hz～50khz；主传感器谐振频率为85hz，带宽为7hz；工频谐波干扰由室内电力线产生，其磁场分布未知；通信信号为85hz单音信号，其发射天线采用空心线圈，线圈环面法线与主传感器磁芯平行；

图4为采用本发明涉及的多参考磁矢量传感系统，主通道输出信号功率谱图，其中fft样点数为16384；

图5为采用本发明涉及的多参考磁矢量传感系统，经过步骤(1)处理后，主通道序列的分组的方差示意图，其中，组长n为25；

图6为采用本发明涉及的基于方差的中值滤波方法，经过第二级自适应中值滤波后，主通道信号波形图；

图7为采用本发明涉及的基于方差的中值滤波方法，经过残存干扰抑制处理后，主通道信噪功率比得到有效提升的示意图，其中，通信信号功率谱密度高于底噪6db。

具体实施方式

以下结合附图和具体实施例对本发明作进一步的详细说明：

图1为本发明涉及的信号处理流程示意图，各步骤详细说明如下：

第一步：接收系统采用多组参考传感器与一个主传感器本地正交布设的方案，构建多参考磁矢量传感系统。其中，主传感器放置于电屏蔽铝壳腔内，u组参考传感器沿主传感器磁芯轴线方向，以一定的空间间隔均匀固定于铝壳外表面，每组参考传感器包含2个正交的传感器，与主传感器构成三轴磁矢量传感器结构；传感器采用小型化、窄带、高灵敏度感应式磁场传感器，调整匹配电容使主传感器谐振于通信信号中心频率。

第二步：多通道同步采样。2u+1路传感器同时接收，经多通道ad转换器实现同步采样，得到多路磁感应输出数据矩阵：x＝[x0,x1,...,x2u]^t，其中x0表示主传感器输出序列，x1,...,x2u表示2u路参考传感器的输出序列，每个序列共k个样点，即xu＝[xu(1),xu(2),...,xu(k)]^t,那么x可表示为：

其中，k表示采样点数，其取值依据实际存储电路的可接受容量确定。

图3所示为采用本专利所述的多参考磁传感系统，主通道接收时域波形，其中，大气噪声由大气噪声模拟器产生，通信信号为85hz单音信号，工频谐波干扰来源于室内电力线和电器设备。

图4所示为主通道功率谱密度，fft点数为16384.

第三步：按照时序，对磁感应输出序列x进行分组，并估计分组的方差。设每分组包含n(n≤k)个样点，以第u通道序列xu为例，分组后的序列表示为：

即

vu,i＝[xu((i-1)n+1),xu((i-1)n+2),...,xu(in)^t(3)

其中，vu表示通道u的分组序列，vu,i表示第i组，i＝1,2,...,i，表示向下取整运算。那么，分组矩阵v表示为：

然后，估计序列样点局部方差。分组vu,i的方差利用下式进行估算：

由此可得方差矩阵为：

其中，di＝[d0,i,d1,i,...,d2u,i]^t表示d的列向量，i＝1,2,...,i.图5所示为主通道分组方差[d0,1,d0,2,...,d0,655]的分布图，考虑到超低频频段工频谐波干扰功率主要来源于50hz的基波、二次谐波、三次谐波以及五次谐波，每个分组中应包含整数个谐波干扰周期，当采样频率为1.25khz时，取n＝25。

第四步：基于大数判决准则，对方差矩阵d进行自适应中值滤波。

首先，检测d中的“脉冲元素”，即脉冲检测阶段。构造中心元素为du,i的滑动窗fu,i＝[du,i-m,...,du,i,...,du,i+m]，窗长l＝2m+1，1≤m≤(i-1)/2，初始m＝3，脉冲元素判决方法由下式表示：

其中，(fu,i)max表示滑动窗fu,i中的最大值，t为判决门限值，通常取t＝0。su,i用于表征du,i的状态，su,i＝0表示du,i为非脉冲元素，su,i＝1表示du,i为脉冲元素。然后，利用大数判决准则，当仅至少u+1条通道(对应于d中u+1行)均检测到脉冲元素时，才判定对应的列矢量是一个“脉冲矢量”。即对于第i分组，如果：

则判定di为脉冲矢量，否则，判定di为非脉冲矢量。

然后，对脉冲矢量进行自适应中值滤波。设滑动窗fu,i的中值为mu,i，按照下式判定mu,i的状态：

其中，wu,i＝0表示mu,i为脉冲元素，wu,i＝1表示mu,i为非脉冲元素。设脉冲矢量记为如果所有通道滑动窗中值mn＝[m0,n,m1,n,...,m2u,n]^t元素均为非脉冲，即

则用mn替换作为中值滤波输出。否则，增大滑动窗长度，m＝m+1，返回脉冲检测阶段。对于非脉冲的矢量，中值滤波输出与输入相同。

上述脉冲检测和自适应中值滤波过程的流程图如图2所示。

最后，对脉冲矢量及其邻近矢量进行第二级自适应中值滤波。窄带接收的极/超低频信号，其大气噪声呈现出被展宽的连续的脉冲波形，表现为非孤立脉冲序列。上述一级滤波过程往往难以充分压制非孤立脉冲点，此时，采用第二级自适应中值滤波过程，对脉冲矢量及其邻近矢量进行第二级自适应中值滤波。设上述第一级中值滤波输出为检测到的脉冲矢量为那么构造新的方差矩阵d′n：

其中，l为正整数且l≤(i-1)/2。以d′n作为中值滤波的输入，重复执行上述脉冲检测和自适应中值滤波处理，得到二级自适应中值滤波输出方差矩阵

第五步：序列重建。比较第二级自适应中值滤波输出方差矩阵与第一级滤波输入方差矩阵d的元素取值，按如下方式确定比例系数τu,i：

由比例系数τu,i对组vu,i进行线性幅度压缩，得到幅值压缩后的分组矩阵：

从而实现滤波输出序列的重建。

重建后的主通道样点序列如图6所示，其中，第二级中值滤波参数l＝10。

第六步：工频谐波干扰和残存大气噪声的抑制。

工频谐波干扰采用iir二阶陷波滤波器进行压制，待陷频率设置为工频谐波频点，所有通道采用相同的滤波参数设置，以确保陷波输出同步。采用级联滤波器形式，第n个陷波器的z域传递函数表达式为：

其中，fn＝50·nhz表示待陷频率，n＝1,2,...,6，fs表示采样频率，参数μ用以控制陷波增益和带宽。第u通道陷波滤波输出序列yu,n(k)为:

yu,n(k)＝xu(k)-2cos(2πfn/fs)xu(k-1)+xu(k-2)+

2(1-μ)cos(2πfn/fs)yu,n(k-1)-(1-μ)²yu,n(k-2)(16)

其中，当p≤0时，yu,n(p)＝0,xu(p)＝0.

采用多通道指数加权递归最小二乘(rls)算法的自适应噪声抵消(anc)技术，对残存的多通道相关的大气噪声进行压制。

经第六步处理后的输出序列的功率谱如图7所示，可以看出，85hz单音信号功率谱密度高于底噪6db，与图4相比，经本专利所述的技术方案处理后，信噪比得到有效提高。

以上所述，仅为本发明的具体实施方式，应当指出，任何熟悉本领域的技术人员在本发明所揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。