CN105022072B

CN105022072B - 北斗卫星定位坐标连续时间序列的混合消噪方法

Info

Publication number: CN105022072B
Application number: CN201510257625.0A
Authority: CN
Inventors: 吴浩; 李奎; 张建华; 池秀文; 董元锋; 吴彩保; 殷亚
Original assignee: Wuhan University of Technology WUT
Current assignee: Wuhan University of Technology WUT
Priority date: 2015-05-19
Filing date: 2015-05-19
Publication date: 2017-12-29
Anticipated expiration: 2035-05-19
Also published as: CN105022072A

Abstract

本发明提供的一种基于小波系数的北斗卫星定位坐标连续时间序列的混合消噪方法，是在考虑到北斗卫星定位坐标连续时间序列中不仅含有白噪声而且含有闪烁噪声的基础之上提出的一种新的消噪方法。具体是：首先估计出与原始坐标连续时间序列(原始信号)中的真实白噪声具有相同方差的人工白噪声及其与原始信号在不同尺度下的小波系数；然后小波系数信息熵之间的关系提取出闪烁噪声的小波系数，从而将原始信号中闪烁噪声消除；最后运用小波软阈值处理进一步消除其中的白噪声。本发明可以有效的解决传统方法只考虑消除白噪声而导致的北斗卫星定位坐标连续时间序列存在的准确性和光滑性等方面的问题，提高了其在矿山高边坡变形监测中的可靠性。

Description

北斗卫星定位坐标连续时间序列的混合消噪方法

技术领域

本发明属于大地测量领域，特别是一种基于小波系数的北斗卫星定位坐标连续时间序列的混合消噪方法。在北斗卫星定位系统的矿山边坡位移变形监测过程中可用于同时消除北斗卫星定位坐标连续时间序列中的闪烁噪声和白噪声，提高北斗卫星定位坐标连续时间序列在边坡位移变形监测中的精度。

背景技术

近二十年来，卫星定位系统GPS(Global Positing System)已经广泛应用于自然灾害的监测，如地震、泥石流以及山体滑坡等。随着我国自主研发的北斗卫星定位系统逐渐成熟，已经完全具备了亚太区域的导航定位功能。但是，作为一种新型卫星导航定位系统，其定位误差特征和消噪方法仍然是空白。北斗卫星定位系统的监测结果的准确性，与卫星时钟、地形条件和大气延迟密切相关，它们都是导致北斗卫星定位系统的连续监测时间序列产生严重的误差的原因。针对这些因素的影响，目前主要从北斗卫星定位系统基线高精度解算的角度提出了一系列处理算法；但是，由于不同时间点的观测值之间缺乏一定的误差相关性，导致长时间北斗卫星定位系统的监测序列仍然残留显著的噪声，严重制约着利用北斗卫星定位系统进行灾害监测的可靠性和准确性。因此，研究如何准确的消除北斗卫星定位坐标连续时间序列中的噪声是非常必要的。

卫星定位系统的坐标连续时间序列的消噪方法研究，长期以来一直是卫星定位系统在变形监测领域的热点。Donoho等人在1994年提出了一种小波阈值消噪方法来消除卫星定位系统的连续坐标时间序列中存在的白噪声，从而大大的提高了卫星定位系统监测的可靠性。这充分的发挥了小波理论在消除北斗卫星定位坐标连续时间序列中消噪的巨大的潜力，具有重大的意义。基于Donoho的理论，许多研究人员也提出了一些新的方法来消除卫星定位系统所采集的信号的质量。如Gao 1997年运用阈值水平来使得重构后的信号尽可能不含有噪声。 Souza和Monico通过运用小波收缩的方法来消除北斗卫星定位坐标连续时间序列中相应位置的高频信号，从而达到消除噪声的目的。虽然他们在对北斗卫星定位坐标连续时间序列消噪方面都做出了巨大的贡献，但是他们几乎都忽略了信号光滑性的重要性。为了弥补这一缺陷，Han等人基于小波软阈值消噪理论提出了一种新的阈值算法来提高信号的光滑性。然而，北斗卫星定位坐标连续时间序列中的噪声和其它的地理现象一样也可以被定义为一个幂率的过程。Agnew(1992)已经指出了北斗卫星定位坐标连续时间序列中的噪声不仅含有独立且均匀分布的白噪声还含有非平稳变化且长期互相关的闪烁噪声。Williams(2004)也指出包含在北斗卫星定位坐标连续时间序列中的噪声可以被充分的描述为白噪声加闪烁噪声的模型。然而，闪烁噪声具有与白噪声完全不同的统计特性，所以运用处理白噪声的传统方法来消除闪烁噪声已经不再适用了。为了解决这一难题，2003年何凯等人利用香农熵理论将闪烁噪声与白噪声分离，从而达到消除闪烁噪声的目的。除此之外，Wornell和Oppenheim在1992年提出了在小波领域内运用极大似然估计来消除闪烁噪声的方法。Chen和Lin 1994年结合正交小波理论和维纳滤波器对北斗卫星定位坐标连续时间序列中的闪烁噪声进行了消除。但是，针对北斗卫星定位坐标连续时间序列中既含有白噪声又含有闪烁噪声这一特性，提出一个能够同时消除北斗卫星定位坐标连续时间序列中的白噪声和闪烁噪声的方法将会使北斗卫星定位系统的监测结果更为准确。

目前，据文献统计，虽然国内外已有众多学者分别针对北斗卫星定位坐标连续时间序列中的白噪声和闪烁噪声提出了很多的消噪方法，但是由于没有同时考虑到北斗卫星定位坐标连续时间序列中的白噪声和闪烁噪声对监测结果的影响，导致这些算法所得到的消噪结果还不是很准确。

发明内容

本发明所要解决的技术问题是：为弥补传统消噪方法将北斗卫星定位坐标连续时间序列中的噪声假设为纯白噪声进行消噪而忽略了闪烁噪声的重要性，无法准确的对北斗卫星定位坐标连续时间序列中的噪声进行消除，提出了一种基于小波系数的北斗卫星定位坐标连续时间序列的混合消噪方法，同时消除北斗卫星定位坐标连续时间序列中的白噪声和闪烁噪声，达到提高北斗卫星定位坐标连续时间序列分析可靠性。

本发明解决其技术问题的技术方案是：

本发明提供的北斗卫星定位坐标连续时间序列的混合消噪方法，是一种基于小波系数的北斗卫星定位坐标连续时间序列的混合消噪方法，该方法是：首先利用已估计的原始北斗卫星定位坐标连续时间序列即原始信号的信噪比来构造人工白噪声，然后通过小波分解分别计算出原始信号和人工白噪声在不同尺度下的小波系数香农熵；在此基础上，根据原始信号和人工白噪声在不同尺度下的小波系数香农熵之间的关系，计算出相应尺度下闪烁噪声的小波系数香农熵，得到不同尺度下闪烁噪声的小波系数；最后运用小波阈值消噪方法进一步消除其中的白噪声，得到既不含有闪烁噪声也不含有白噪声的北斗卫星定位坐标连续时间序列。

在构造人工白噪声的过程中，利用估计出原始信号的信噪比来得到真实白噪声的方差，从而运用matlab的rand函数构造出与真实白噪声具有相同方差的人工白噪声来代替真实白噪声。

在获取小波系数香农熵的过程中，先通过小波分解分别得到原始信号S1和人工白噪声在不同尺度下的小波系数，然后通过香农熵的定义由下述公式分别计算出原始信号S1和人工白噪声W在不同尺度下的小波系数香农熵，

式中：原始信号记为S1，为原始信号的小波系数模极大值，ψ(j,k)为基本小波，j为尺度因子，k为小波系数中模极大值的总数，代表为原始信号的小波系数香农熵。

在获取小波系数香农熵的基础上，根据原始信号和人工白噪声在不同尺度下的小波系数香农熵之间的关系，由下述公式来提取出属于闪烁噪声的小波系数，得到不同尺度下闪烁噪声的小波系数，以消除北斗卫星定位坐标连续时间序列中的闪烁噪声，

式中：和分别代表人工白噪声和闪烁噪声的小波系数香农熵，和分别代表人工白噪声和闪烁噪声的小波系数模极大值，T_NF代表门限值的大小。

所述T_NF门限值可以通过在不同尺度上与之间的比值变化来确定，当尺度a对应的比值发生急剧增大的情况时T_NF的值就为a。

本发明可以采用包括以下步骤的方法来消除北斗卫星定位坐标连续时间序列中的闪烁噪声：

(1)通过计算闪烁噪声与原始信号在不同尺度下小波系数香农熵的比例，计算出不含闪烁噪声的北斗卫星定位坐标连续时间序列S2在不同尺度上的小波系数，

式中：为闪烁噪声在不同尺度下的小波系数，代表闪烁噪声与原始信号在不同尺度下小波系数香农熵的比例；

(2)利用小波重构将不同尺度的小波系数进行重构得到不含闪烁噪声的北斗卫星定位坐标连续时间序列S2，

式中：S2为不含闪烁噪声的北斗卫星定位坐标连续时间序列。

本发明可以采用包括以下步骤的方法获取既不含有闪烁噪声也不含有白噪声的北斗卫星定位坐标连续时间序列S3：

(a)运用小波分解将不含闪烁噪声的北斗卫星定位坐标连续时间序列S2分解为不同尺度下的小波系数，将北斗卫星定位坐标连续时间序列由时间域转化到频率域；

(b)设置阈值，将含有白噪声的高频域去除，

式中：T_j代表给定的阈值，为信号S2在不同尺度下的小波系数，为信号S3在不同尺度下的小波系数，σ_j为尺度j下的小波系数的方差，N为所采集的样本点的总数；

(c)运用小波重构将不同尺度下的小波系数重构得到既不含有闪烁噪声也不含有白噪声的北斗卫星定位坐标连续时间序列S3。

本发明提供的方法与现有技术相比，具有以下的主要的优点：

(1)相比于传统的消噪方法而言，本方法不仅消除了具有独立且均匀分布的特点的白噪声，而且去除了具有非平稳变化且长期互相关的特点的闪烁噪声，使得消噪后的时间序列的光滑性得到明显的改善，信噪比较传统的消噪方法至少提高5倍。

(2)针对北斗卫星定位坐标连续时间序列中同时含有白噪声和闪烁噪声的特点，本方法采用在白噪声消除之前先利用香农熵理论将闪烁噪声从混合噪声中消除，这样极大的提高了小波阈值消噪方法对白噪声的去除效率。

附图说明

图1是本发明的方法流程图。

图2是原始信号S1在时间域和频率域下的示意图。

图3是消除闪烁噪声之后的信号S2在在时间域和频率域下的示意图。

图4是同时消除了闪烁噪声和白噪声之后的信号S3在在时间域和频率域下的示意图。

图5是传统消噪方法在消除噪声之后的信号S4在在时间域和频率域下的示意图。

图6是不同信号的信噪比(SNR)的对比结果图。

具体实施方式

下面结合实施例及附图对本发明作进一步说明，但并不局限于下面所述内容。

本发明针对传统消噪方法在对北斗卫星定位坐标连续时间序列进行噪声消除的过程中，将噪声简单的假设为纯白噪声进行消除而忽略了闪烁噪声的影像这一问题，提出了基于小波系数的北斗卫星定位坐标连续时间序列的混合消噪方法。该方法是：首先利用已估计的原始北斗卫星定位坐标连续时间序列(即原始信号S1，如图2所示)的信噪比来构造人工白噪声，然后通过小波分解分别得出原始信号和人工白噪声在不同尺度下的小波系数香农熵；在此基础上，根据信息熵理论，计算出在不同尺度上闪烁噪声的小波系数香农熵，利用闪烁噪声与原始信号的小波系数香农熵之间的关系将北斗卫星定位坐标连续时间序列中闪烁噪声消除得到信号S2(如图3所示)；最后运用小波的分解与重构消除其中的白噪声得到信号S3(如图4所示)。

本发明提供的上述方法包括以下的步骤：

一、数据准备

为了验证发明中所提出基于小波系数的北斗卫星定位坐标连续时间序列的混合消噪方法的消噪效果，我们选择了一组包含512个数据点的北斗卫星定位坐标连续时间序列进行实验。这组数据选择的是我国陕西金堆城钼矿在2012年在该矿区建立了一个北斗卫星定位连续监测网中的一个监测点(BB05)在2013年7月6日到2013年9月29号期间采集的北斗卫星定位坐标连续时间序列。该监测网共包含40个监测点，每个监测点每隔3个小时就会采集一次三维坐标(包括x、y和z方向)。由于系统误差、测量误差和偶然误差的存在，坐标时间序列中存在着许多噪声使得这些监测点所采集的监测点的位移时间序列不能够直接用于边坡的稳定性评价。图2a给出的是监测点BB05在2013年7月6日到2013年9月29号期间所采集的在z方向的位移变形数据(即原始信号S1)。显然，它包含着严重的误差以至于监测点的变形趋势无法确定。为了能够更加清楚的看到在效果过程中的信号质量的变化，我们运用傅立叶变换将信号由时间域转换到空间域。经过傅立叶变换之后，有用信号通常集中在低频段而噪声往往集中在高频部分。由图2b可以看出相对于低频部分(<＝50Hz)的振幅，高频部分(>＝50Hz)的振幅非常的明显。这说明了原始信号S1的质量严重受到噪声的影响。

二.北斗卫星定位坐标连续时间序列中闪烁噪声的消除

将上面准备好的原始信号S1和构造的人工白噪声进行小波分解得到不同尺度下的小波系数，根据信息熵理论，利用原始信号和人工白噪声的小波系数香农熵计算出闪烁噪声的小波系数香农熵。通过闪烁噪声与原始信号小波系数香农熵之间的比例关系得到闪烁噪声在不同尺度下的小波系数，从而得到消除闪烁噪声之后的信号S2，如图3所示。

具体包括以下步骤：

1.编写matlab代码：

基于小波系数的北斗卫星定位坐标连续时间序列的混合消噪方法，编写matlab代码。

北斗卫星定位坐标连续时间序列中闪烁噪声的消除主要以下步骤的计算方法得到：

(1)构造人工白噪声：

运用频率域估计方法来估计出原始信号的信噪比，再通过得到的信噪比可以计算出白噪声的方差，最后可以利用matlab中的rand函数来随机生成一个与真实白噪声具有相同方差的人工白噪声，其matlab代码如下：

White_noise＝(rand(1,L)-0.5)*sqrt(12*var)+mean

(2)计算闪烁噪声的小波系数香农熵：

对原始信号和人工白噪声分别进行小波分解，得到不同尺度小的小波系数，然后根据香农熵的定义计算出在不同尺度下原始信号和人工白噪声小波系数香农熵，具体的计算公式如下所示。其中，原始信号记为S1，为原始信号的小波系数模极大值，ψ(j,k)为基本小波，j为尺度因子，k为反应位移，代表为原始信号的小波系数香农熵

最后由原始信号和人工白噪声小波系数香农熵之间的关系计算出闪烁噪声在不同尺度下的小波系数香农熵。由下述公式来提取出属于闪烁噪声的小波系数，得到不同尺度下闪烁噪声的小波系数，以消除北斗卫星定位坐标连续时间序列中的闪烁噪声。其中，和分别代表人工白噪声和闪烁噪声的小波系数香农熵，和分别代表人工白噪声和闪烁噪声的小波系数模极大值，T_NF代表门限值的大小。

(a)小波分解的matlab代码如下：

load S1；

s＝S1；

％用db4小波函数对信号进行九尺度小波分解

[C,L]＝wavedec(s,9,'db4')；

figure(1)；

plot(s)；

title('原始信号')；

％提取尺度1的低频系数

cA1＝appcoef(C,L,'db4',1)；％用小波分解框架[C.L]计算1层低频系数的近似值，小波基为db4

％提取尺度2的低频系数

cA2＝appcoef(C,L,'db4',2)；％用小波分解框架[C.L]计算2层低频系数的近似值，小波基为db4

％提取尺度3的低频系数

cA3＝appcoef(C,L,'db4',3)；％用小波分解框架[C.L]计算3层低频系数的近似值，小波基为db4

％提取尺度4的低频系数

cA4＝appcoef(C,L,'db4',4)；％用小波分解框架[C.L]计算4层低频系数的近似值，小波基为db4

％提取尺度5的低频系数

cA5＝appcoef(C,L,'db4',5)；％用小波分解框架[C.L]计算5层低频系数的近似值，小波基为db4

％提取尺度6的低频系数

cA6＝appcoef(C,L,'db4',6)；％用小波分解框架[C.L]计算6层低频系数的近似值，小波基为db4

％提取尺度7的低频系数

cA7＝appcoef(C,L,'db4',7)；％用小波分解框架[C.L]计算7层低频系数的近似值，小波基为db4

％提取尺度8的低频系数

cA8＝appcoef(C,L,'db4',8)；％用小波分解框架[C.L]计算8层低频系数的近似值，小波基为db4

％提取尺度9的低频系数

cA9＝appcoef(C,L,'db4',9)；％用小波分解框架[C.L]计算9层低频系数的近似值，小波基为db4

figure(2)；

subplot(2,1,1)；

plot(cA1)；

title('尺度1的低频系数')；

subplot(2,1,2)；

plot(cA2)；

title('尺度2的低频系数')；

subplot(2,1,3)；

plot(cA3)；

title('尺度3的低频系数')；

subplot(2,1,4)；

plot(cA4)；

title('尺度4的低频系数')；

plot(cA1)；

title('尺度1的低频系数')；

subplot(2,1,5)；

plot(cA5)；

title('尺度5的低频系数')；

subplot(2,1,6)；

plot(cA6)；

title('尺度6的低频系数')；

subplot(2,1,7)；

plot(cA7)；

title('尺度7的低频系数')；

subplot(2,1,8)；

plot(cA8)；

title('尺度8的低频系数')；

subplot(2,1,9)；

plot(cA9)；

title('尺度9的低频系数')；

％提取尺度1的高频系数

cD1＝detcoef(C,L,1)；％用小波分解框架[C.L]计算1层高频系数的近似值，小波基为db4％提取尺度2的高频系数

cD2＝detcoef(C,L,2)；％用小波分解框架[C.L]计算2层高频系数的近似值，小波基为db4％提取尺度2的高频系数

cD3＝detcoef(C,L,3)；％用小波分解框架[C.L]计算3层高频系数的近似值，小波基为db4％提取尺度2的高频系数

cD4＝detcoef(C,L,4)；％用小波分解框架[C.L]计算4层高频系数的近似值，小波基为db4％提取尺度2的高频系数

cD5＝detcoef(C,L,5)；％用小波分解框架[C.L]计算5层高频系数的近似值，小波基为db4％提取尺度2的高频系数

cD6＝detcoef(C,L,6)；％用小波分解框架[C.L]计算6层高频系数的近似值，小波基为db4％提取尺度2的高频系数

cD7＝detcoef(C,L,7)；％用小波分解框架[C.L]计算7层高频系数的近似值，小波基为db4％提取尺度2的高频系数

cD8＝detcoef(C,L,8)；％用小波分解框架[C.L]计算8层高频系数的近似值，小波基为db4％提取尺度2的高频系数

cD9＝detcoef(C,L,9)；％用小波分解框架[C.L]计算9层高频系数的近似值，小波基为db4figure(3)；

subplot(3,1,1)；

plot(cD1)；

title('尺度1的高频系数')；

subplot(3,1,2)；

plot(cD2)；

title('尺度2的高频系数')；

subplot(3,1,3)；

plot(cD3)；

title('尺度3的高频系数')；

subplot(3,1,4)；

plot(cD4)；

title('尺度4的高频系数')；

subplot(3,1,5)；

plot(cD5)；

title('尺度5的高频系数')；

subplot(3,1,6)；

plot(cD6)；

title('尺度6的高频系数')；

subplot(3,1,7)；

plot(cD7)；

title('尺度7的高频系数')；

subplot(3,1,8)；

plot(cD8)；

title('尺度8的高频系数')；

subplot(3,1,9)；

plot(cD9)；

title('尺度9的高频系数')；

(b)计算小波系数香农熵的matlab代码如下：

(3)得到不含有闪烁噪声的北斗卫星定位坐标连续时间序列：

通过闪烁噪声和原始信号在不同尺度上的小波系数香农熵之比可以得到其小波系数之比，从而同构重构不好又闪烁噪声的小波系数就能得到不含有闪烁噪声的北斗卫星定位坐标连续时间序列S2。具体的计算公式如下所示。其中，为闪烁噪声在不同尺度下的小波系数，代表闪烁噪声与原始信号在不同尺度下小波系数香农熵的比例。

Matlab代码如下：

c1＝[cA8

cD8,(-0.5442)*cD7,(-0.5242)*cD6,(-0.4037)*cD5,(-0.1049)*cD4,(-0.2899)*cD3,(-0.0140)

*cD2,(-0.1189)*cD1]；

S2＝waverec(c1,L,'db4')；

figure(4)；

plot(S2)；

三.北斗卫星定位坐标连续时间序列中闪烁噪声的消除

本步骤针对步骤二中最终得到的不含有闪烁噪声的北斗卫星定位坐标连续时间序列S2，运用小波阈值消噪对其中的白噪声进行消除，最终得到既不含有闪烁噪声也不含有白噪声的北斗卫星定位坐标连续时间序列S3，具体的步骤如下所示。

(b)设置阈值，将含有白噪声的高频域去除。

(c)运用小波重构将不同尺度下的小波系数重构得到既不含有闪烁噪声也不含有白噪声的北斗卫星定位坐标连续时间序列S3。其中，T_j代表给定的阈值，为信号S2在不同尺度下的小波系数，为信号S3在不同尺度下的小波系数，σ_j为尺度j下的小波系数的方差，N为所采集的样本点的总数。

其具体的matlab代码如下：

load s2；

％用db4小波对信号s2进行8层分解并提取系数

[c,l]＝wavedec(s2,8,'db4')；

a8＝appcoef(c,l,'db4',8)；

d8＝detcoef(c,l,8)；

d7＝detcoef(c,l,7)；

d6＝detcoef(c,l,6)；

d5＝detcoef(c,l,5)；

d4＝detcoef(c,l,4)；

d3＝detcoef(c,l,3)；

d2＝detcoef(c,l,2)；

d1＝detcoef(c,l,1)；

thr＝1；

％进行阈值处理

ytsoftd8＝wthresh(d8,'s2',thr)；

ytsoftd7＝wthresh(d7,'s2',thr)；

ytsoftd6＝wthresh(d6,'s2',thr)；

ytsoftd5＝wthresh(d5,'s2',thr)；

ytsoftd4＝wthresh(d4,'s2',thr)；

ytsoftd3＝wthresh(d3,'s2',thr)；

ytsoftd2＝wthresh(d2,'s2',thr)；

ytsoftd1＝wthresh(d1,'s2',thr)；

c8＝[a3ytsoftd8ytsoftd7ytsoftd6ytsoftd5ytsoftd4ytsoftd3ytsoftd2ytsoftd1]；

s3＝waverec(c8,l,'db4')；

plot(s3)；

title('不含闪烁噪声和白噪声的北斗卫星定位坐标连续时间序列s3')；

本发明提供的上述方法可以有效解决传统消噪方法将北斗卫星定位坐标连续时间序列中的噪声假设为纯白噪声而忽略了闪烁噪声的重要性而导致的消噪结果在不精确的问题。通过同时考虑闪烁噪声和白噪声消除，诠释了闪烁噪声对消噪效果的影响，提高了北斗卫星定位坐标连续时间序列的精度。下面以具体实例说明：

1.以2012年在陕西省金堆城钼矿边坡上建立的40个边坡变形监测点所采集的北斗卫星定位坐标连续时间序列为基础数据，选择在地形变化最具有代表性的BB05监测点在2013年7月6日到2013年9月29号期间采集的在竖直方向的512个连续坐标时间序列作为实例研究的基础数据。

2.通过matlab代码对采集的监测数据进行粗差的提取，去除粗差之后的北斗卫星定位坐标连续时间序列(即原始信号S1)如图2a所示。

3.运用频率域估计方法估计出原始信号S1的信噪比，通过信噪比构造出与真实白噪声具有相同方差的人工白噪声。

4.利用小波分解分别得到原始信号和人工白噪声在不同尺度下的小波系数。

5.根据香农熵的定义分别计算出原始信号和人工白噪声在不同尺度下的小波系数香农熵，再利用原始信号和人工白噪声的小波系数香农熵之间的关系计算出在不同尺度下闪烁噪声的小波系数，从而通过小波重构得到不含有闪烁噪声的北斗卫星定位坐标连续时间序列S2，如图3中的上图所示。

6.将步骤5中得到的不含有闪烁噪声的北斗卫星定位坐标连续时间序列S2导入到matlab中，用db4小波对信号S2进行8层分解并提取小波系数，然后对其中的高频系数进行软阈值处理，最后通过小波重构来得到既不含有闪烁噪声也不含有白噪声的北斗卫星定位坐标连续时间序列S3，如图4中的上图所示。

消噪结果的对比分析：

功率谱密度是直观的描述北斗卫星定位坐标连续时间序列中有用信号和噪声的重要指标，能够清晰的反映北斗卫星定位坐标连续时间序列消噪过程中其中噪声的变化。因此，为了能够更加清楚的看到北斗卫星定位坐标连续时间序列的信号质量的变化，我们运用傅立叶变换将信号由时间域转换到空间域。经过傅立叶变换之后，有用信号通常集中在低频段(<50Hz)而噪声往往集中在高频部分(>＝50Hz)。

以2012年在陕西省金堆城钼矿北帮边坡上的监测点BB05在2013年7月6日到2013年9月29号期间采集的在竖直方向的512个北斗卫星定位坐标连续时间序列S1为例，采用基于小波系数的北斗卫星定位坐标连续时间序列的混合消噪方法来去除其中的白噪声和闪烁噪声，主要分为两个步骤：(1)基于香农熵理论在白噪声背景下删除闪烁噪声，得到不含有闪烁噪声的北斗卫星定位坐标连续时间序列(图3中的上图)。(2)再运用小波阈值消噪方法消除其中的白噪声，得到既不含有闪烁噪声也不含有白噪声的北斗卫星定位坐标连续时间序列S3(图4中的上图)。

通过图2a和图3中的上图的对比我们可以发现虽然信号S2在光滑性方面较信号S1没有很大的提高，但是从图3中的下图的傅立叶变换结果中可以看出信号S2在高频部分(>＝50Hz)的振幅较信号S1有明显的减小。这表明在消除原始信号中的闪烁噪声之后信号的质量有着明显的提高。图4中的上图是运用小波阈值消噪方法从信号S2中消除白噪声之后的结果。信号S3的光滑性较信号S2(图2a)有了明显的提高，可以明显的看出监测点的变形趋势。除此之外，图4中的下图与图2b中的高频部分(>＝50Hz)信号S3与信号S2的振幅大小也有着鲜明的对比。这说明了信号S3的质量较信号S2有了很大的提高。同时这一结果也有利的证明了北斗卫星定位坐标连续时间序列中的噪声主要是由白噪声和闪烁噪声组成。为了更加全面的证明消除闪烁噪声的必要性，我们选择了一个传统的消噪方法的消噪结果进行对比。因为在消除信号S2中的白噪声时我们选用了小波阈值消噪法，为了使对比结果更具说服力，在这一部分我们就选小波阈值消噪法进行消噪。图5中的上图是我们运用小波阈值消噪法消噪后得到的消噪结果信号S4。通过对比图4中的上图和图5中的上图我们可以发现传统消噪方法的消噪结果与本文所提出的消噪方法有着明显的不同。毫无疑问信号S4的光滑信号明显不如信号S3。在图5中的下图中给出的信号S4的傅立叶变化图我们也可以发现在低频部分振幅基本保持不变的情况下，信号S4在高频部分(>＝50Hz)振幅的下降明显不如信号S3。这一结果表明了信号S1中包含着白噪声和其它噪声的组合，而不是仅仅只含有白噪声。在此基础上我们可以得出在消除北斗卫星定位坐标连续时间序列中噪声的过程中考虑出噪声外其它类型噪声的存在是非常有必要的。

虽然在上述的讨论部分已经对不同信号的光滑性以及高频部分的振幅进行了分析，但是它不能非常直观的表达出不同信号的质量变化。为了能够进一步的了解在消噪过程中不同信号的质量变化，我们将运用不同信号的信噪比再一次进行比较。图6给出的是在消噪过程中出现的4个不同信号的信噪比的比较。从中可以看出信号S3的信噪比高达24.6319，是四个信号中信噪比最好的。这表明了本文中提出的算法所得到的消噪结果要比传统的消噪结果要好。同时我们可以发现信号S2的信噪比等于1.5472，仅仅是原始信号的2倍而已，然而信号S3的信噪比却是信号S2的近16倍。运用传统消噪方法得到的消噪结果信号S4仅仅是原始信号S1的6倍。造成这一现象是由于我们提前在消除白噪声之前提前将原始信号中的白噪声进行了消除，这就使得剩余在信号S2中的白噪声的特性更加的明显从而提高了随后的白噪声消除的效率。此外，在原始信息中的噪声都存在着相互的影响，也就是说一个噪声不仅会对原始信号中有用的信号产生影响，也会对其它类型的噪声产生潜在的影响。这一结果说明了提前对原始信号中的闪烁噪声的消除也有利于白噪声的去除。

Claims

1.一种北斗卫星定位坐标连续时间序列的混合消噪方法，其特征是一种基于小波系数的北斗卫星定位坐标连续时间序列的混合消噪方法，该方法是：首先利用已估计的原始北斗卫星定位坐标连续时间序列即原始信号的信噪比来构造人工白噪声，然后通过小波分解分别计算出原始信号和人工白噪声在不同尺度下的小波系数香农熵；在此基础上，根据原始信号和人工白噪声在不同尺度下的小波系数香农熵之间的关系，计算出相应尺度下闪烁噪声的小波系数香农熵，得到不同尺度下闪烁噪声的小波系数；最后运用小波阈值消噪方法进一步消除其中的白噪声，得到既不含有闪烁噪声也不含有白噪声的北斗卫星定位坐标连续时间序列；

<mrow> <msubsup> <mi>r</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>=</mo> <msup> <mn>2</mn> <mrow> <mo>-</mo> <mi>j</mi> <mo>/</mo> <mn>2</mn> </mrow> </msup> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <mi>S</mi> <mn>1</mn> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mi>&psi;</mi> <mrow> <mo>(</mo> <msup> <mn>2</mn> <mrow> <mo>-</mo> <mi>j</mi> </mrow> </msup> <mi>t</mi> <mo>-</mo> <mi>k</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>t</mi> </mrow>

<mrow> <msub> <mi>E</mi> <mi>n</mi> </msub> <mo>&lsqb;</mo> <msubsup> <mi>r</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>&rsqb;</mo> <mo>=</mo> <munder> <mo>&Sigma;</mo> <mi>k</mi> </munder> <msup> <mrow> <mo>(</mo> <msubsup> <mi>r</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> <msub> <mi>log</mi> <mn>2</mn> </msub> <msup> <mrow> <mo>(</mo> <msubsup> <mi>r</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow>

式中：原始信号记为S1，为原始信号的小波系数模极大值，ψ(j,k)为基本小波，j为尺度因子，k为小波系数中模极大值的总数，代表为原始信号的小波系数香农熵；

<mrow> <msub> <mi>E</mi> <mi>n</mi> </msub> <mo>&lsqb;</mo> <mover> <msubsup> <mi>x</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>&OverBar;</mo> </mover> <mo>&rsqb;</mo> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>E</mi> <mi>n</mi> </msub> <mo>&lsqb;</mo> <mover> <msubsup> <mi>r</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>&OverBar;</mo> </mover> <mo>&rsqb;</mo> <mo>-</mo> <msub> <mi>E</mi> <mi>n</mi> </msub> <mo>&lsqb;</mo> <mover> <msubsup> <mi>&omega;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>&OverBar;</mo> </mover> <mo>&rsqb;</mo> </mrow> </mtd> <mtd> <mrow> <mo>(</mo> <mi>j</mi> <mo><</mo> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mi>F</mi> </mrow> </msub> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>E</mi> <mi>n</mi> </msub> <mo>&lsqb;</mo> <mover> <msubsup> <mi>r</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>&OverBar;</mo> </mover> <mo>&rsqb;</mo> </mrow> </mtd> <mtd> <mrow> <mo>(</mo> <mi>j</mi> <mo>&GreaterEqual;</mo> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mi>F</mi> </mrow> </msub> <mo>)</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>

2.根据权利要求1所述的北斗卫星定位坐标连续时间序列的混合消噪方法，其特征是在构造人工白噪声的过程中，利用估计出原始信号的信噪比来得到真实白噪声的方差，从而运用matlab的rand函数构造出与真实白噪声具有相同方差的人工白噪声来代替真实白噪声。

3.根据权利要求1所述的北斗卫星定位坐标连续时间序列的混合消噪方法，其特征是公式中的T_NF门限值可以通过在不同尺度上与之间的比值变化来确定，当尺度a对应的比值发生急剧增大的情况时T_NF的值就为a。

4.根据权利要求1所述的北斗卫星定位坐标连续时间序列的混合消噪方法，其特征是采用包括以下步骤的方法消除北斗卫星定位坐标连续时间序列中的闪烁噪声：

<mrow> <msubsup> <mi>&omega;</mi> <mi>i</mi> <msup> <mi>j</mi> <mo>,</mo> </msup> </msubsup> <mo>=</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mover> <msub> <mi>R</mi> <mi>j</mi> </msub> <mo>&OverBar;</mo> </mover> <mo>)</mo> </mrow> <msubsup> <mi>r</mi> <mi>i</mi> <mi>j</mi> </msubsup> </mrow>

<mrow> <mover> <msub> <mi>R</mi> <mi>j</mi> </msub> <mo>&OverBar;</mo> </mover> <mo>=</mo> <msub> <mi>E</mi> <mi>n</mi> </msub> <mo>&lsqb;</mo> <mover> <msubsup> <mi>x</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>&OverBar;</mo> </mover> <mo>&rsqb;</mo> <mo>/</mo> <msub> <mi>E</mi> <mi>n</mi> </msub> <mo>&lsqb;</mo> <mover> <msubsup> <mi>r</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>&OverBar;</mo> </mover> <mo>&rsqb;</mo> </mrow>

式中：为闪烁噪声在不同尺度下的小波系数，为闪烁噪声与原始信号在不同尺度下小波系数香农熵的比例，i为反应位移；

<mrow> <mi>S</mi> <mn>2</mn> <mo>=</mo> <munder> <mo>&Sigma;</mo> <mi>j</mi> </munder> <munder> <mo>&Sigma;</mo> <mi>i</mi> </munder> <msup> <mn>2</mn> <mrow> <mo>-</mo> <mi>j</mi> <mo>/</mo> <mn>2</mn> </mrow> </msup> <msubsup> <mi>&omega;</mi> <mi>i</mi> <msup> <mi>j</mi> <mo>,</mo> </msup> </msubsup> <mi>&psi;</mi> <mrow> <mo>(</mo> <msup> <mn>2</mn> <mrow> <mo>-</mo> <mi>j</mi> </mrow> </msup> <mi>t</mi> <mo>-</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow>

5.根据权利要求1所述的北斗卫星定位坐标连续时间序列的混合消噪方法，其特征是采用包括以下步骤的方法获取既不含有闪烁噪声也不含有白噪声的北斗卫星定位坐标连续时间序列S3：

(b)设置阈值，将含有白噪声的高频域去除，

<mrow> <msubsup> <mi>&omega;</mi> <mi>i</mi> <msup> <mi>j</mi> <mrow> <mo>&prime;</mo> <mo>&prime;</mo> <mo>&prime;</mo> </mrow> </msup> </msubsup> <mo>=</mo> <msub> <mi>&eta;</mi> <mi>S</mi> </msub> <mrow> <mo>(</mo> <msubsup> <mi>&omega;</mi> <mi>i</mi> <msup> <mi>j</mi> <mrow> <mo>&prime;</mo> <mo>&prime;</mo> </mrow> </msup> </msubsup> <mo>,</mo> <msub> <mi>T</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msubsup> <mi>&omega;</mi> <mi>i</mi> <msup> <mi>j</mi> <mrow> <mo>&prime;</mo> <mo>&prime;</mo> </mrow> </msup> </msubsup> <mo>-</mo> <msub> <mi>T</mi> <mi>j</mi> </msub> </mrow> </mtd> <mtd> <mrow> <mo>(</mo> <msubsup> <mi>&omega;</mi> <mi>i</mi> <msup> <mi>j</mi> <mrow> <mo>&prime;</mo> <mo>&prime;</mo> </mrow> </msup> </msubsup> <mo>&GreaterEqual;</mo> <msub> <mi>T</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <mo>(</mo> <mo>|</mo> <msubsup> <mi>&omega;</mi> <mi>i</mi> <msup> <mi>j</mi> <mrow> <mo>&prime;</mo> <mo>&prime;</mo> </mrow> </msup> </msubsup> <mo>|</mo> <mo><</mo> <msub> <mi>T</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msubsup> <mi>&omega;</mi> <mi>i</mi> <msup> <mi>j</mi> <mrow> <mo>&prime;</mo> <mo>&prime;</mo> </mrow> </msup> </msubsup> <mo>+</mo> <msub> <mi>T</mi> <mi>j</mi> </msub> </mrow> </mtd> <mtd> <mrow> <mo>(</mo> <msubsup> <mi>&omega;</mi> <mi>i</mi> <msup> <mi>j</mi> <mrow> <mo>&prime;</mo> <mo>&prime;</mo> </mrow> </msup> </msubsup> <mo>&le;</mo> <mo>-</mo> <msub> <mi>T</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>

<mrow> <msub> <mi>T</mi> <mi>j</mi> </msub> <mo>=</mo> <msub> <mi>&sigma;</mi> <mi>j</mi> </msub> <msqrt> <mrow> <mn>2</mn> <mi>ln</mi> <mi> </mi> <mi>N</mi> </mrow> </msqrt> <mo>,</mo> </mrow>