绝对重力仪振动特性研究*

胡若，吴书清，牟宗磊，冯金扬，牛学犇

（1.中国计量科学研究院 北京，100029）

（2.山东科技大学电气与自动化工程学院 青岛，266590）

引言

高精度绝对重力场信息在地球科学研究、资源勘探、导航与制导等方面有着重要的意义，目前主要依靠绝对重力仪进行当地重力场信息的直接测量［1-3］。自由落体式绝对重力仪通过测量落体自由下落时相对于参考棱镜的位移和时间间隔，数据拟合求解得到测点的重力加速度［4］。在测量绝对重力值的过程中，影响重力测量精度的地面振动主要分为环境噪声和仪器噪声两种。其中，环境噪声的产生与测量时间和测量地点的坐标、地基以及周围环境密切相关，可以通过频率进行区分：①人类活动引起的噪声（频率一般大于1 Hz）；②地脉动噪声（频率一般为0.1～1 Hz）；③大气运动噪声（频率一般小于0.1 Hz）［5］。仪器噪声主要来自于落体上升、下降时电机以及传动机构产生的振动［6］。如何进一步提高绝对重力仪的测量精度、减小振动因素干扰已经成为该领域的一个重要研究方向［7-8］。

绝对重力仪的振动信号是极为复杂的非稳态时变信号，传统的时域或频域分析方法很难精细化提取出振动的特征信息，不利于振动特性分析以及后期的振动补偿和高精度减振装置的设计研发。Huang等［9］提出EMD 是一种自适应分解方法，可以将非平稳信号自适应分解为一系列具有不同物理意义的IMF，从时频域对信号进行分析。EMD 方法虽然具有诸多优点，在很多领域得到了广泛应用，但存在的一个主要问题是模态混叠［10］。为了解决模态混叠问题，总体经验模态分解（ensemble empirical mode decomposition，简称EEMD）和补充的总体经验模态分解（complete ensemble empirical mode decomposition，简 称CEEMD）被相应提出［11-12］。EEMD 通过向原始信号中多次加入不同的白噪声进行EMD 分解，将多次分解的结果进行平均运算后得到最终的IMF。CEEMD 算法主要是通过向原始信号中成对添加两个相反的白噪声信号再进行EMD 分解，CEEMD 在保证分解结果与EEMD 相近的情况下，减小了由原始信号噪声部分造成的重构误差，但是计算量大，如果添加白噪声幅值和迭代次数不合适，会产生很多伪分量，需要对IMF 进行重新组合或者进一步处理［13］。文献［14］提出了改进集总平均经验模态分解（modified ensemble empirical mode decomposition，简 称MEEMD）解决EMD 分解过程中的模态混叠和CEEMD 计算量过大等问题，基于排列熵的随机性检测，剔除分解过程中的异常信号，但是分解过程中仍然可能存在虚假分量，需要后续处理。HHT 是一种线性变换，它强调信号的局部性质，在处理非线性、非平稳性信号时具有更加直观的物理意义［15］。

基于以上问题，笔者重点开展了绝对重力仪振动特性的相关研究，并对MEEMD 方法进行了改进，提出了一种基于改进MEEMD-HHT 的绝对重力仪振动特性分析方法，对实测重力仪振动数据进行自适应分解后，从时频瞬时特征和能量谱特征等方面对振动数据进行了分析研究。

1 NIM-3A 型绝对重力仪测量原理

基于自由落体原理的NIM-3A 型绝对重力仪如图1 所示。

图1 NIM-3A 型绝对重力仪Fig.1 NIM-3A type absolute gravimeter

测量过程中，落体角锥棱镜在真空环境中自由下落，通过激光干涉仪测量其在下落过程中的距离s，用高精度时间间隔计数器测量对应的下落时间间隔t，通过数据拟合算法计算实测的距离si和时间间隔ti（i=1，2，…，n），再通过固体潮、大气压、振动和极移等参数综合修正得到所测点的绝对重力值［16］。其中，落体自由下落的距离通过激光干涉条纹数N乘以激光半波长λ2 计算得到［17］。通过检测干涉条纹过零点的时间，可得到干涉条纹每个过零点对应的时间-位移坐标信息（si，ti）。在实际测量过程中，重力场的梯度和光速有限性会对测量结果产生一定的影响，考虑到落体自由下落距离比较短，因此设定重力场梯度值在整个下落过程中为常量。对于光速有限性导致的条纹生成时间滞后于落体下落时间的问题，进行时间修正，则自由下落运动轨迹可表示为

其中：=t0-(si-s0)/c为延迟时间；γ为测量处的重力梯度值；c为光速。

更进一步，式（1）可以写成如下矩阵形式

式（2）可以简化为

利用最小二乘拟合算法可以得出

其中：S为不同时刻落体位移组成的向量矩阵；H为位移对应的时间间隔矩阵；为利用最小二乘法得到的最佳估计。

通过以上几式可以看出，在绝对重力值测量时，除去干涉带提取信息算法的误差外，测量点的背景振动、落体伺服控制系统的振动是误差的主要来源，影响测量精度。因此，有效采集振动信息，对振动特性进行深入研究，了解振动各频段的能量分布情况，采取合适的隔振措施，对测量进行合理的振动补偿是提高绝对重力值测量精度的有效途径。

2 改进的MEEMD-HHT 算法

2.1 改进的MEEMD 算法

基于文献［18］中的研究，笔者所提的改进MEEMD 算法自适应分解原始信号S（t）步骤如下。

1）在原始信号S（t）中成对添加两组绝对值相等的正负白噪声信号ωn+(t)和ωn-(t)，分别对其进行集总平均次数相等的EEMD 自适应分解可得

2）将分解得到的两组IMF 分量均值运算后，可极大地消除白噪声在信号中的残差

3）分解得到的di(t)分量并不一定是标准的IMF 分量，可能是存在模态混叠问题的异常信号，根据排列熵值［19］和互相关系数判定信号是异常信号还是平稳信号。

4）如果判断di(t)是异常信号，需要对这组分量继续进行EMD 自适应分解，剔除虚假信息以得到标准的IMF 分量，并消除模态混叠问题。以第1个IMF 分量进行EMD 自适应分解为例

其中：c1(t)为d1(t)经过EMD 自适应分解得到的第1 个标准IMF 分量；q1(t)为残差分量；ck(t)为分解得到的第k个标准IMF 分量；qk(t)为对应的残差分量；k=2，3，…，m。

5）将得到的IMF 分量按高频到低频依次排列，则最终的改进MEEMD 自适应分解可以表示为

其中：r(t)为改进MEEMD 自适应分解得到的最终残差分量。

通过改进MEEMD 自适应分解得到的各IMF分量，可以剔除绝对重力仪振动信号中的虚假分量，得到振动信号在不同频段的时频分布规律。

2.2 基于HHT 的能量谱特征参数提取

对EMD 自适应分解得到的每个IMF 分量进行HHT 后，原始信号可以表示为

式（11）中省略了残差R，把信号X（t）用幅值与频率皆为时间函数t的n个分量之和来表示。同理，信号X（t）用Fourier 级数可以表示为

其中：ai（t），ωi（t）为时间变量t的函数。

信号幅值在三维空间中既可以用时间-瞬时频率的函数来表示，也可以用时间-频率平面上的等高线来表示，时间-频率平面上的幅值分布即Hilbert 时频谱［20］，其表达式为

H（ω，t）对时间t积分运算后得到Hilbert 边际谱

在频域范围内对H（ω，t）的平方积分运算后得到Hilbert 瞬时能量谱［21］为

瞬时能量可以表述信号能量随时间的变化规律，而在时域内H（ω，t）的平方对时间t积分运算得到Hilbert 能量谱

通过Hilbert 能量谱可以求出绝对重力仪振动信号各个频段在整个时间长度内所积累的能量，可以得到不同频段的能量分布规律。

3 算法仿真分析与应用

为了验证所提算法的有效性和先进性，既接近真实的信号特性又不失一般性，笔者通过多种特性信号组合成的模拟信号来对CEEMD 和改进的MEEMD 两种算法的性能进行对比分析。采用的模拟信号X1～X3及其组合而成的合成信号如图2 所示。其中：X1是频率为30 Hz 的正弦周期信号，用来模拟信号中的周期主频部分；X2是频率为8 Hz 的调幅信号，用来模拟信号中不同幅值、不同频段的有效信息；X3为间歇信号，用来模拟实测信号中的突变部分，合成信号Signal=X1+X2+X3。

图2 模拟信号及其合成信号Fig.2 The simulation signal and its synthesis signal

添加30 对正负白噪声，噪声幅值为原始信号标准差（standard deviation，简称SD）的0.2 倍，集总平均次数为200 次，两种算法自适应分解结果如图3、图4 所示。

通过图3、图4 的对比分析可以看出：

图3 改进的MEEMD 自适应分解结果Fig.3 The decomposition results of improved MEEMD

图4 CEEMD 自适应分解结果Fig.4 The decomposition results of CEEMD

1）由于间歇信号的影响，CEEMD 分解得到的高频IMF1，IMF2，IMF3分量中含有明显的模态混叠现象和端点效应，同时分解过程中出现了很多虚假分量，高频IMF 部分未能把间歇信号很好地剔除出去，因此在分解含有间歇信号、脉冲信号这样的瞬态信号过程中，CEEMD 分解很难得到正确的IMF分量。

2）在采用同样噪声幅值的白噪声信号和同样的集总平均次数下，改进的MEEMD 剔除了信号中的间歇性噪声信号，准确得到了有效的IMF，抑制了模态混叠、端点效应和噪声残留的问题，更适用于含有间歇类噪声的非线性、非平稳性瞬态类信号的自适应分解。

4 实测信号分析与应用

为了验证改进的MEEMD-HHT 算法在绝对重力仪振动特性分析中的实际应用效果，用本研究方法对实测的振动数据进行了分析研究。实验在中国计量科学院昌平基地重力比对实验室进行，实验平台采用基于激光干涉法的NIM-3A 型绝对重力仪，主要由激光干涉仪、自由落体装置、数据采集装置、拾振器和测控终端组成，如图5 所示。其中，拾振器位于绝对重力仪正下方，与激光干涉仪直接相连，与重力仪主机部分相对独立。

图5 硬件实验平台Fig.5 The hardware experiment platform

实验过程中激光干涉条纹信号和振动数据同步采集，其中：激光干涉条纹由PXI-5124 采集，带宽为150 MHz，采样频率为200 MHz，AD 量化位数为12-Bit；振动数据由VSE-311M 拾振器同步获取，采样频率为500 Hz，带宽为0.03～70 Hz，分辨率为5×10-8g。根据信号特性和采样定理的相关分析，选用的实验平台满足系统测试需求。

相同实验条件下，共进行了6 组实验，每组实验做20 次自由落体运动。实验结束后，随机抽取了一次实验的激光干涉条纹信号和对应的振动信号如图6 所示，其中图6（a）为激光干涉条纹信号，图6（b）为重力仪振动信号。

图6 激光干涉条纹信号和振动信号Fig.6 The laser interference fringes and vibration signal

在绝对重力仪的测试过程中，激光干涉条纹信息和振动数据同步采集，对应交点时刻落体上的时间-位移坐标受到的振动是线性相关的，本次实验对应重力加速度g的拟合残差如图7 所示。

图7 拟合残差Fig.7 The fitting residual error

从图6、图7 的分析可以看出，测量过程中的振动噪声导致重力加速度g的拟合过程中存在残差，幅值范围在±3 nm 以内。

为了更好地研究环境噪声和仪器噪声引入的振动误差问题，在重力仪不工作的状态下，通过拾震器采集了3 组实验室环境下的地面振动，采样频率为500 Hz，采样时间为4.5 s，结合重力仪工作时的振动数据，进行了对比分析。其中，重力仪不工作时的3 组地面振动数据如图8（a）所示，随机抽取重力仪工作时的3 组振动数据如图8（b）所示。

通过图8（a，b）两幅图的对比分析可以看出：重力仪不工作时，拾振器采集的3 组振动数据具有一定的周期性和规律性，振动幅值较小，变化平稳；重力仪工作时振动信号的前半段（电机未启动、落体未运动时），振动较平缓；1.26 s 左右电机启动、落体开始运动时，振动明显加强。

图8 不同振动信号的对比分析Fig.8 Comparative analysis of different vibration signals

鉴于排列熵［22］算法可以验证时间序列的复杂程度，在此引入排列熵来验证以上6 组振动信号的复杂性。经过计算得到重力仪不工作时3 组振动数据的排列熵分别为0.667 3，0.673 6 和0.660 3，而重力仪工作时3 组振动数据的排列熵分别为0.851 6，0.844 8 和0.854 8，可见由于仪器振动的加入，排列熵值明显增大。综上，测量过程中的仪器噪声加剧了振动对于测量精度的影响，在进行隔振系统设计时需要关注仪器的振动噪声和周围环境的振动噪声。

分别采用CEEMD 和改进的MEEMD 对图6 选取的振动信号进行了自适应分解，限于篇幅，只列出MEEMD 的分解结果如图9 所示。其中，CEEMD和MEEMD 中添加的噪声幅值和对数分别为0.2和100。

图9 改进MEEMD 自适应分解实测振动信号结果Fig.9 The results of vibration signal decomposed by improved MEEMD

对自适应分解得到的各IMF 进行HHT 变换后，所得希尔伯特-黄谱如图10 所示，其中右边的彩条代表信号的能量分布情况，颜色越深，能量越大，反之则越小。

图10 希尔伯特-黄谱Fig.10 The Hilbert-Huang spectrum

通过图9、图10 的分析可以看出：改进的MEEMD-HHT 对实测信号自适应分解后得到6 个IMF（IMF1～IMF6）和1 个余项R7，分解的分量个数较为合理，信号频率较低，幅值较小。从希尔伯特-黄谱时频特性分析可以看出，信号主要集中在70 Hz 以内，大部分为低频信号。

鉴于振动特性与各频段的能量分布密切相关，振动在某频段内的能量越大，对自振频率在相同频段内仪器的振动影响越强。为了更加直观地观测所测振动信号各频段内能量的高低分布情况，在此将频率划分为多个频段后，用能量直方图展现各频段内的能量分布，如图11 所示。

图11 各频段的能量分布Fig.11 Energy distribution in each frequency band

为了研究能量分布的规律性，对以上6 组实验的数据进行归一化处理后，从每组实验中随机抽取1 次实验中的振动数据，用改进的MEEMD-HHT 算法自适应分解后，提取各IMF 的能量谱如图12 所示。

通过图11、图12 的分析可以看出：重力仪的振动特性具有较好地复现性和一致性，能量谱可以很好地展现出各个频段的能量分布；振动噪声主要由低频段的地脉动噪声（0.1～5 Hz）和部分高频段的自振噪声组成，幅值较小；振动噪声的能量主要集中在低频段（IMF4～IMF6），衰减速度快，高频段的能量较小，而衰减速度相对较慢。因此，设计在超低频段具有更好隔振效果的隔振装置以克服地面振动的影响，研究系统自振产生的原因并消除系统自振的影响，是提高绝对重力仪测量精度的有效手段。

图12 各IMF 的能量谱Fig.12 Energy spectrum of the IMFs

5 结束语

在自由落体式激光干涉绝对重力仪的测量过程中，参考棱镜在背景振动和系统自振的作用下呈现复杂的振动特性，影响绝对重力仪的测量精度。笔者在借鉴现有技术的基础上，对绝对重力仪的振动特性进行了分析研究，提出了一种基于改进MEEMD-HHT 的绝对重力仪振动特性分析方法。通过模拟仿真信号验证了该方法在自适应分解瞬态类、非平稳信号方面的优势，有效抑制模态混叠的同时还能剔除无效的IMF 分量，信号重构误差为10-15量级，实现了分解结果的准确性。采用改进的MEEMD 算法对实测绝对重力仪振动数据进行自适应分解，选取有效IMF 后，进行HHT 变换，得到了各频段振动信号的时频特征和能量谱特征，0.1～20 Hz 低频段能量占比95.73%，从能量的角度对振动特性进行了分析研究。本研究成果对进一步提高绝对重力仪的测量精度，以及更高精度的隔振系统的设计研发具有一定的指导意义。