原子钟数据的改进经验模态分解降噪

惠 恬，赵高长，苏 军，龚莹莹

(西安科技大学 理学院，陕西 西安 710600)

0 引言

精准的时间尺度需要根据参与计算的原子钟模型及其频率稳定度来设计[1-2]，由于铯原子钟自身存在的噪声会对其输出信号的频率稳定度和精度造成影响，因此在计算时间尺度之前，对铯原子钟数据进行消噪处理以提高其频率稳定度，具有重要意义。

目前，可以通过滤波的方法对原子钟的噪声进行处理。例如，文献[3]应用小波方法、Vondrak方法和卡尔曼(Kalman)滤波等7种方法，对原子钟数据进行消噪处理，均达到了比较好的降噪效果，其中小波方法的效果最好。在实际应用中，Kalman滤波[4-5]和小波方法[6-7]都有一定的缺陷。Kalman滤波的过程过于复杂，小波分析则需要选择一个合适的小波函数，分解层数还要选择合适的阈值。经验模态分解(empirical mode decomposition，EMD)方法是一种新的自适应信号处理技术，采用该方法对原子钟数据进行降噪处理，可以消除随机噪声的影响，提高原子钟的短期稳定度，获得更加精确的时间尺度，同时避免了Kalman滤波和小波阈值方法的缺点。文献[8]采用EMD方法，对高性能5071A原子钟信号进行了去噪和频率预测。文献[9-10]在文献[8]的基础上，采用EMD方法对铯原子钟的白色调频噪声进行去噪处理，并应用经过处理后的原子钟数据进行时间尺度计算，研究了EMD方法对合成频率稳定度的影响。然而，传统的EMD方法将信号分解成多个固有模函数(intrinsic mode function，IMF)分量后，将含有噪声的高频IMF分量舍弃，再将剩余的低频IMF分量与残差合成作为平滑后的原子钟数据。实际上，每一个IMF分量中既有噪声也有有用信号，简单地舍弃一个或多个IMF分量，会导致相应的有用信号也同时被舍弃，造成信号的严重失真，影响降噪效果。

因此，针对EMD方法分解后的噪声分量筛选及处理的问题，本文在传统EMD方法的基础上，结合小波阈值降噪法对其进行改进。针对原子钟频率数据非线性、非平稳的特点，利用窗口做出划分，对每一个窗口进行降噪处理。分别从时域和频域对该方法的降噪效果做出分析，并与传统的EMD方法及常见的小波阈值降噪方法做出比较，以期提升原子钟数据降噪效果。

1 原子钟数据的模型

铯原子钟的核心部件是原子频标，在实际运行的过程中，原子钟受到内部噪声及测量噪声的影响，导致原子钟的输出信号存在相位偏差和频率偏差[11-12]，降低了原子钟的频率稳定度。相位数据和频率数据都可以用来描述原子钟的频率稳定度，且相位数据转化为频率数据，可以通过计算相位数据的一次差分除以取样时间获得，即：

(1)

其中：yi为频率数据；xi和xi+1分别为i时刻和i+1时刻对应的相位数据；τ0为取样时间间隔。

通常情况下，原子钟的相位相对变化量x(t)和相对瞬时频率偏差y(t)，都是由确定性变化部分和随机变化部分组成，表达式分别见式(2)和式(3)：

(2)

y(t)=y0+Dt+εy(t)，

(3)

其中：x0为初始相位偏差；y0为初始频率偏差；D为线性频漂；εx(t)和εy(t)为随机变化量，即随机噪声部分。关于原子钟的随机变化部分，目前国际上普遍认为原子钟的随机噪声模型，可以看作是5种独立的噪声线性叠加而得到的符合幂律谱模型的经验公式[13-14]。

铯原子钟数据具有非线性、非平稳的特点。设含噪声的铯原子钟频率数据为X(t)，可以表达为：

X(t)=f(t)+e(t)，

(4)

其中：X(t)为含有噪声的原子钟数据；f(t)为有用数据；e(t)为噪声数据。为了提高结果的置信度，充分利用铯原子钟数据，采取重叠ALLAN方差对原子钟的频率稳定度进行评估。

2 降噪方法

EMD方法是在Hilbert-Huang变换(Hilbert-Huang transform，HHT)的基础上，将信号分解至不同频段，是一种循环迭代的算法[15-16]。EMD方法将原始信号分为若干个IMF和一个余项之和，不同频率的IMF分量按频率由高到低的顺序排列，每一个IMF可以看作是一个函数且满足以下条件：

(Ⅰ)在所有的数值序列中，最大点数和过零点数必须相等或最多相差1；

(Ⅱ)在任意点上，最大极值点和最小极值点的包络的平均值为0。

EMD方法分解与原始信号都在同一尺度的时域中，保持了信号的非平稳性，具有简单、高效、自适应的特点。在频域上，经EMD分解出的IMF分量的频率几乎按2的负幂次方的形式递减[17]，有用信号比较平稳，对应IMF的低频分量，而含有噪声的信号波动较大，对应IMF的高频分量。因此，可以通过频谱分析找出高频IMF分量，利用小波阈值对高频IMF分量进行二次处理后再重构，得到降噪信号。利用改进EMD方法对原子钟频率数据降噪，具体步骤如下：

(Ⅱ)找出子窗口内数据Xi所有的局部极大值点和极小值点，利用3次样条插值分别连接所有的极大值点和极小值点，作为上下包络线u1(t)和v1(t)，计算上下包络线的包络均值m1(t)曲线，

(5)

(Ⅲ)计算窗口内数据Xi和包络均值m1(t)的差作为第一分量h1(t)，

h1(t)=X(t)-m1(t)。

(6)

(Ⅳ)若该分量没有满足IMF分量的两个条件，则将该分量看作是窗口内的数据，重复第Ⅱ步，得到h11(t),h12(t),…，h1k(t)，直到h1k(t)=h1(k-1)(t)-m1k(t)满足IMF分量的两个条件，称h1k(t)为一阶IMF分量，是该窗口内数据Xi中最高频率分量，记为：

c1(t)=h1k(t)。

(7)

(Ⅴ)将c1(t)从初始信号中分离，得到第1个剩余分量r1(t)，

r1(t)=X(t)-c1(t)。

(8)

(Ⅵ)由于r1(t)中仍存在较长周期的原始信号分量，因此重复上述步骤，得到：

(9)

(Ⅶ)直到rn(t)变成单调序列或不再满足IMF调节时结束[18]，得到经过EMD分解后的窗口内的数据，可以表示为：

(10)

(Ⅷ)对IMF分量做出频谱分析，将噪声IMF分量筛选出来。

(Ⅸ)利用软阈值函数对噪声IMF分量进行小波阈值降噪，软阈值函数为：

(11)

(Ⅹ)将降噪处理后的IMF分量与信号分量进行信号重构，得到降噪后的窗口内数据为：

(12)

(Ⅺ)在每个窗口内重复以上步骤，得到降噪后的数据，

(13)

为了验证改进的EMD方法的有效性，分别从信噪比(signal-to-noise ratio,SNR)、均方根误差(root mean square error,RMSE)两方面进行评价[19-20]。信噪比的数值越大，说明降噪效果越好。均方根误差的数值越小，说明降噪效果越好。

3 实证分析

在实际计算中，选取一个月的铯原子钟数据，采样间隔τ0=1 h，在应用铯原子钟频率数据前，已对原始数据做出预处理。

采用EMD方法对铯原子钟频率数据进行分解，如图1所示，将铯原子钟频率数据分解为8个IMF分量和1个残差。对IMF分量做出频谱分析，如图2所示。由图2可以看出：通过EMD分解后，各个IMF分量的频率由高到低依次递减，IMF1中包含了大部分噪声。根据信号与IMF分量的对应关系，可以看出IMF1～IMF5为高频分量即噪声分量，其余为低频分量即信号分量。

图1 EMD分解后的IMF分量

在小波阈值降噪中，阈值选取和阈值函数的选择十分重要，不同的选取方式对去噪效果有很大影响。常见的阈值选取方式有基于无偏估计原理的自适应阈值(rigrsure)、固定阈值(sqtwolog)、启发式阈值(heursure)和极大极小阈值(minimax)。分别根据4种不同的阈值选取方式对高频IMF分量做处理，由于处理结果存在的差异不是非常明显，需要利用信噪比及均方根误差两个评价指标来判定降噪结果，见表1。根据降噪的评价指标可知：选用极大极小阈值进行处理的结果最好。

表1 不同阈值选取方式的降噪结果

表2 不同分组的降噪效果对比

根据表2，将铯原子钟频率数据分为5组分别做处理，去噪前后对比如图3所示，蓝色线条对应原始数据，红色线条对应降噪数据。由图3可以看出：改进EMD降噪后，噪声得到了一定的抑制且保留了原始数据的主要图像特征。从频域上分析，分别做出降噪前后的频谱，如图4所示，其中，图4a为降噪前的原始信号频谱，图4b为降噪后的信号频谱。由图4可以看出：经过降噪处理后，波动明显减小，说明噪声能够得到一定压制。同时从时域上分析，为了提高结果的置信度，充分利用铯原子钟数据，采取重叠ALLAN方差对原子钟的频率稳定度进行评估。计算了降噪前后铯原子钟频率数据对应不同时间间隔的重叠ALLAN方差，如图5所示，蓝色线条代表数据降噪前对应不同时间间隔的重叠ALLAN方差，红色线条为数据降噪后对应不同时间间隔的重叠ALLAN方差。由图5可以看出：相同的时间间隔，通过改进EMD方法降噪的数据频率稳定度小于原始数据的频率稳定度。综上所述，改进后的方法能够有效提高其频率稳定度。

图3 原始数据与降噪数据

(a) 原始信号频谱

分别应用改进后的EMD方法、传统的EMD方法及小波阈值方法，对同一组铯原子钟数据做降噪处理，结果如图6所示。图6a为通过传统EMD方法降噪的对比结果图，原子钟频差数据经过EMD分解后，将高频的IMF分量舍弃，剩余的低频IMF分量重构作为降噪后的信号。由图6a可以看出：简单地去除一个或多个IMF分量的传统EMD方法，会导致去除的IMF分量中的有用信息也被去掉，该方法较为粗糙，只能保持原始数据的大概趋势。图6b和图6c分别是小波阈值算法及改进EMD方法的降噪结果的对比图。由图6b和图6c可以看出：两种方法都可以对噪声产生一定的抑制作用，但是效果的差异并不是很明显。利用信噪比及均方根误差两个评价指标来判定降噪结果，见表3。由表3可以发现：改进后的EMD方法降噪的效果优于传统的EMD方法及小波阈值降噪，说明该降噪方法可以有效降低原子钟频率数据中的噪声。

图6 不同方法降噪结果对比

表3 不同方法降噪的结果对比

4 结论

(1)传统EMD方法分解后的噪声数据筛选会造成部分有用信号的浪费，进而导致原数据的失真，在传统EMD方法的基础上，结合小波阈值对噪声数据进行二次降噪，可以极大程度地保留噪声数据中含有的少数有用数据，避免了有用数据的浪费及原始数据失真的情况。

(2)利用改进EMD方法对铯原子钟频差数据进行降噪，既保留了传统EMD方法的优点，又避免了Kalman滤波计算的复杂性和小波基函数的选择，具有较好的自适应性，适用于具有非线性、非平稳特点的铯原子钟数据，减少了噪声对时间尺度计算的影响，进一步得到精确的时间尺度，具有很强的实用价值。

(3)对原始数据进行加窗处理，效果比直接对数据进行处理的效果好，窗口数越多，得到的效果越好，但当数据分组超过一定数量时，由于每一组的数据量过少得到的效果反而很差。对噪声数据进行二次降噪时，选用极大极小值的信噪比为3.453 9，均方根误差为7.582 4e-14，均优于其他3种阈值选择方式，说明选用极大极小值的效果更好。

(4)改进EMD方法与传统EMD方法和小波阈值方法相比较，信噪比从0.676 1和3.321 8提高到3.652 3，均方根误差从1.044 0e-13和7.698 6e-14降低到7.411 1e-14，说明该方法可以抑制原子钟频差数据中的噪声，进一步提高了原子钟的频率稳定度。