随着我国BDS、俄罗斯GLONASS和欧洲Galileo等卫星定位系统迅速发展，基于全球导航卫星定位系统的位置服务已经从单一GPS模式发展到多模GNSS模式阶段^[1]。相对于单一的GPS，GNSS各个卫星定位系统在卫星定轨和信号质量等方面具有一定的差异性，导致了GNSS时间序列中的误差特性更加复杂^[2]，已经成为制约GNSS深入应用的瓶颈问题。特别是对于大坝变形、桥梁振动、边坡滑移和板块运动等长周期高精度的变形监测用户来说，GNSS时间序列的误差在时间轴线上还存在着一定积累的效应，造成的影响不可低估。

早期研究普遍认为GNSS时间序列的噪声主要包括白噪声和闪烁噪声^[3]。随着卫星导航定位技术的应用研究深入，人们逐渐认识到GNSS时间序列中还含有显著的随机漫步噪声^[4-5]。文献[6—7]分别采用谱指数估计方法，证明了GNSS时间序列中存在着以随机漫步特征为主的噪声。文献[8]利用极大似然估计法来探索坐标时间序列的噪声成分，也同样发现了部分GNSS时间序列中含有随机漫步噪声。文献[9]通过艾伦方差分析法，进一步证实了GNSS时间序列中随机漫步噪声的存在。文献[10]运用误差分析的方法，发现了美国中东部连续GNSS时间序列水平分量中含有较大随机漫步噪声。因此，有必要对GNSS时间序列中的随机漫步噪声开展系统的消噪研究，才能对它加以有效地控制^[11]。

时间序列的消噪方法一直都是GNSS研究的热点领域，已经涌现出卡尔曼滤波法、模极大值法和小波阈值法等一批经典的理论^[12]，其中小波阈值消噪法因简捷有效，最具有代表性^[13]。文献[14—15]较早地基于小波理论提出了硬阈值和软阈值两种消除白噪声的方法。针对硬阈值的消噪结果方差较大与软阈值有固定偏差的双重问题，文献[16]对阈值函数进行改进而提出了半软阈值法，较好地降低了消噪结果失真情况。21世纪以来，小波阈值消噪法有了新的进展。文献[17]基于小波理论与模糊数学相结合的方法，对半软阈值法的阈值计算过程进行优化，减少了阈值计算量，进一步扩大了半软阈值法的应用范围。文献[18]提出了一种基于二次曲线的半软阈值函数，大大缩减了阈值计算工作量。文献[19]对小于阈值部分采用缓变到零的消噪策略代替了直接置零方法，确保了有用信号在消噪过程中得到保留。上述半软阈值法的相关研究，其阈值函数主要以线性收缩策略为主，虽然对于以白噪声和闪烁噪声为主的GNSS时间序列具有较好的适应性^[20]，而对于具有量级小、频率低、敏感性差等复杂特征的随机漫步噪声，直接使用它们容易造成阈值计算量大，且存在着消噪效果不佳等问题。因此，有必要对现有半软阈值函数进行有效的改进，以提高对随机漫步噪声的适应性。

本文提出一种改进的半软阈值算法，重点解决阈值函数的连续性处理和平滑性的参数率定两个关键问题，形成了一套统一的半软阈值函数体系，提高常规半软阈值对随机漫步噪声的适应性。

1 改进半软阈值算法的基本原理 1.1 常规半软阈值法

半软阈值法是经典的一种小波阈值消噪方法，小波阈值消噪的过程是对分解后的小波系数进行阈值处理，然后将处理后的小波系数进行小波重构从而得到消噪后的信号。常规半软阈值法对小波系数进行分段处理，一般需要确定两个阈值，其表达式为

式中，ω_{i, j}为原始信号经小波分解后的小波系数; 为阈值函数处理后的小波系数；sgn (ω_{i, j})为符号函数; T₁、T₂为两个阈值(T₂>T₁)。

1.2 改进半软阈值法

1.2.1 阈值函数改进基本策略

针对阈值函数的改进，主要是将文献[18]中阈值函数二次曲线的优势和文献[19]中阈值函数缓变置零的处理策略进行有效结合，然后从阈值函数的连续性和平滑性两个角度进行处理以达到改进目的。

文献[18]在ω_{i, j}>T的函数部分采用二次曲线形式，符合小波系数中噪声的分布特点，消噪效果明显。该阈值函数其中ω_{i, j}>T的函数部分为

式中，T为阈值；a、b、c、K的计算方法为, K= max (|ω_{i, j}|)。该函数由于采用二次曲线，其平滑性较常规半软阈值有很强的优势，如图 1中改进后的实线曲线和改进前的虚线折线所示。但是，由于该阈值函数采用的是固定形式，无法依据GNSS时间序列中随机漫步噪声强度进行自我调节，直接应用效果并不好。因此，改进策略是引入一个可变参量，通过增加GNSS时间序列自动控制与调节功能，达到对随机漫步噪声的适应性。

文献[19]在|ω_{i, j}|≤T的函数部分没有将小波系数置零，而是缓变到零，这样有利于保留部分信号的小波系数。该阈值函数其中|ω_{i, j}|≤T的函数部分为

式中，T为阈值；N为调节因子(取正整数)。当N取值较小时，经阈值函数处理后的小波系数会偏大，考虑到随机漫步噪声的量级较小特点，往往会导致随机漫步噪声的消除不彻底。因此，在图 1中的|ω_{i, j}|≤T函数部分，改进策略是减去参量h(ω_{i, j})，达到降低该部分函数曲线波动幅度，从而实现对随机漫步噪声的有效控制。

综合上述两点改进策略，初步构建出改进后的阈值函数

式中，T为阈值。

1.2.2 阈值函数的连续性处理

由于改进的半软阈值法还是采用分段阈值函数，其在临界点处并不连续；考虑到随机漫步噪声一般都比较小，而且阈值函数连续时，对于经过阈值函数处理后得到的信号，可以更多地保留原始信号的连续的形态性质，最大限度地保留信号质量，所以有必要将分段函数连接成一个连续的函数曲线。这就要求分段阈值函数在彼此相邻的临界点处的函数值相等，也就是要求在临界点处的阈值函数值的差值为零。记两段函数在阈值处的差值为Δ，即

如果两段函数在阈值处连续，则在ω_{i, j}=T时，Δ_T=0。现在令Δ_T=0，得h(ω_{i, j})=T/(2N+1)；同理令Δ_-T=0，得到h(ω_{i, j})=-T/(2N+1)。由此可得，当|ω_{i, j}|≤T时，h(ω_{i, j})=ω_{i, j}/(2N+1)。

1.2.3 平滑性参数的确定

平滑性用于描述每一段函数曲线的重要特征。改进的半软阈值法不仅要求各个临界点是连续的，而且还要确保整个阈值函数的平滑性，因此需要将|ω_{i, j}|≥T和|ω_{i, j}|≤T的两部分函数中的参数进行统一。将式(2) 中[1-T²/[2T-K²]]的参数K与式(3) 中的参数N统一为一个可变因子m，形成一个完整的半软阈值函数算法体系。虽然式(2) 中的参数K与式(3) 中的参数N在各自阈值函数中分别是独立的变量，但是将它们统一为一个可变因子并不会影响其阈值函数，只是为了增加改进阈值函数的适应性，从而确定出最佳的消噪效果。由于式(3) 中参数N的取值为正整数，为了继承该阈值函数的优势，将可变因子的取值范围也确定为正整数。在具体计算过程中，采用消噪后信号的信噪比达到最大或者均方根误差达到最小作为评价标准，此时所对应的m值为最佳平滑性参数，如图 2所示。

综上所述，最终改进后的阈值函数表达式为

式中，m为可变因子(取正整数)。

1.3 评估方法

为了进一步验证本文所提出方法的可靠性，采用3种方法对改进算法的实际效果进行评估。第1种方法是利用最为常用的信噪比与均方根误差两个指标作为改进效果的评价标准。第2种方法采用时间序列的自相关系数进行评价，通过时间序列时域自相关计算结果可以判别出噪声的时间相关性，从而定性判断噪声成分^[21-22]。第3种方法是采取谱指数估计的分析方法，噪声序列谱指数K的分布范围一般在-3~1之间，当K=-2时，表现为随机漫步噪声；K=0时为白噪声，K=-1时为闪烁噪声，因此噪声序列的谱指数可以定量验证噪声成分^[23-24]

2 试验 2.1 数据

为了验证改进算法的消噪效果，采用文献[25]中GNSS日解算坐标时间序列模型来生成模拟仿真数据。其表达式为

式中，a为初始位置常数项；b为线性速率；c、d、e和f分别为年、半年周期项系数；t_i为时间；g_i为阶跃；H为阶梯函数；v_i为残差。从简化试验角度考虑，没有考虑半年周期项系数e和f，试验数据全部由Matlab平台生成，共1700个历元，如图 3所示。其中，图 3(a)中S0是GNSS时间序列，不含任何误差，坐标时间序列的参数取值依次为a=0、b=0.003 4和c =2。图 3(b)中S1是基于时间相关性特征来生成与S0同步的随机漫步噪声序列，并通过了相应谱指数指标的检验。图 3(c)中S2则是将上述GNSS时间序列S0和随机漫步噪声序列S1相叠加而成。

2.2 结果

2.2.1 常规半软阈值法的消噪结果

小波消噪过程中常用的阈值确定准则有混合准则(heursure)、无偏风险估计准则(rigrsure)、固定阈值准则(sqtwolog)和极大极小准则(minimaxi)。对于常规半软阈值法，采用无偏风险估计和固定阈值准则分别确定两个阈值 T₁和T₂，分别为0.674 3和3.872 1，其消噪结果如图 4中S3所示。

2.2.2 改进半软阈值的消噪结果

由于改进后的半软阈值只需确定一个阈值，本文采用产生均方误差较小的极大极小阈值选取准则来进行计算，其阈值T =2.371 7。然后，采用信噪比达到最大和均方根误差达到最小作为评价标准，分别计算m从1至50过程中两个指标的具体数值，整体上看两个评价指标结果较为一致，如图 5所示。当m=2时，信噪比最低，均方根误差最大，两者结果都表明此时消噪效果最差；当m=8时，信噪比最大，均方根误差最小，消噪效果最好；而且从8至50过程中，信噪比和均方根误差都保持稳定，消噪效果具有一定持续性。因此，可变因子取值为m=8。其实，这个结果同图 2可变因子的理论收敛过程是十分吻合的：当m取值较小时，函数曲线偏离原始信号比较远，从而导致信噪比与均方根误差偏差较大，当m取值较大时，函数曲线趋于稳定，所对应的信噪比与均方根误差趋于稳定。

此时，利用式(6) 进行消噪处理，得到改进半软阈值法的结果S4，如图 6所示。

3 分析与讨论 3.1 改进前后的信噪比/均方根误差变化分析

信噪比和均方根误差是最为常用的消噪效果评价指标。信噪比越大、均方根误差越小，则表明消噪后的信号中残留的噪声越少，消噪的效果越好；反之则效果越差。从图 4和图 6可以初步看出，改进后方法得到的消噪结果比改进前更加平滑，其信号特征保留更多。从表 1中可以进一步看出，改进半软阈值法相对于常规半软阈值法，信噪比从15.846提高到17.672 9，均方根误差从0.57下降到0.46，这都说明改进策略是有一定效果的。相对于现有成果中的白噪声或者闪烁噪声的消噪效果来说，虽然两个具体指标提升的绝对值并不大，但是考虑到随机漫步噪声自身数量级并不大这一实际情况，本文提出的改进算法效果是完全可以接受的。

3.2 改进前后消除噪声的形态特征讨论

为了从定性角度探索所消除噪声类型，采用时间域自相关函数分析其消除噪声的形态特征。分别利用上述常规半软阈值法的消噪结果S3和改进半软阈值法的消噪结果S4，同原始含噪信号S2进行求差运算处理，得到两种方法所对应消除的噪声序列分别为S5和S6，如图 7所示。从所消除噪声序列的波形对比角度来看，图 7的S6和图 3的S1并不是完全吻合的，在局部区间存在着一定的差异，这主要是噪声仿真的试验过程中对随机漫步噪声特征只能在整体上进行控制、而局部区间内并不能绝对控制所造成的。但是，从整体上看，它们在整体上还是具有类似的趋势。为了量化说明这个趋势，下面进一步开展了它们的自相关分析，结果如图 8所示。

图 8所示是利用改进前后方法消除的噪声序列和仿真随机漫步噪声序列所求得自相关系数。从图 8可知，改进前的常规半软阈值法所消除的噪声序列S5的自相关系数中，当延迟时间较小时( < 10 d)，噪声时间序列的自相关系数迅速下降至0；其后，随着延迟时间的不断增大，自相关系数函数始终在0附近波动，上述特征说明其中主要含有时域不相关的白噪声成分。改进后的半软阈值法所消除的噪声序列S6的自相关系数中，可以看出其主要含有与时域相关的有色噪声成分。进一步将改进后的方法消除的噪声S6与仿真的随机漫步噪声S1进行对比分析，可以看出在500 d以内这种相关性是非常显著的，这说明了S6跟S1在整体上具有类似的趋势。虽然随着延迟时间的加大，超过500 d以后(接近试验数据时间跨度的1/3)，相关性的确会减弱，主要是由于自相关分析方法所采用有限的时间序列造成的，并不说明数据整体的差异大。因此，可以初步确定该有色噪声的主要成分是随机漫步噪声。

3.3 改进前后消除噪声类型的定量识别

上述内容从定性的角度说明了改进方法的消噪效果，下面从定量角度来验证所消除噪声是否为随机漫步类型，为此分别计算改进前后消除噪声的谱指数大小。功率谱密度与频率的对数已经被证明服从反比的关系，当信号功率谱密度对数与对应频率对数的离散点回归直线斜率近似为-2时，就可以认为该信号是随机漫步噪声。图 9所示为改进前后方法消除的噪声序列和仿真随机漫步噪声序列的功率谱密度曲线以及计算拟合直线，其谱指数k依次为改进前-1.602 6、改进后-1.925 1和仿真的随机漫步噪声-1.965 1。对比来看，改进后方法所消除噪声，较改进前的更加接近于仿真的随机漫步噪声，这表明改进方法可以更好地消除随机漫步噪声。

4 结论

本文围绕GNSS时间序列中随机漫步噪声消除问题，提出了一种改进的半软阈值算法及其评估，通过验证分析，得到以下结论：

(1) 通过解决阈值函数的连续性处理和平滑性的参数确定两个关键问题，对常规半软阈值函数进行有效改进，在随机漫步消噪方面表现出优异性能；

(2) 从信噪比与均方根误差两个消噪指标、自相关系数和谱指数等多方评估结果来看，改进方法所消除掉的噪声类型属于随机漫步噪声；

(3) 通过引入的可变因子m ，形成了统一的针对随机漫步的半软阈值消噪方法体系，具有更广泛的潜在应用范围。

本文主要是从随机漫步噪声的特征角度入手提出相应的消除噪声方法，由于GNSS时间序列的噪声与其观测环境有相当大的关系，特别是随着GNSS系统更为深入应用的背景下，今后仍需要从噪声产生机理角度来研究GNSS时间序列消噪方法。