Dynamic adaptability of digital fringe processing for the laser interference absolute gravimeter and its shipboard applications evaluation
-
摘要: 激光干涉绝对重力仪干涉信号处理算法的动态适应性研究是进行动态绝对重力测量的基本前提。本文基于希尔伯特变换、直接正交法、过零点法三种瞬时相位处理算法，提出了基于背景振动物理量的合理动态约束条件，即反向的振动速度不能超过落体的下落速度。构建的单频振动信号仿真实验表明：在满足约束条件时，采样频率为60 MHz，希尔伯特变换算法的精度优于10−13 m/s2，过零点算法的精度优于10−9 m/s2，直接正交算法的精度为（−7.9±2.0）×10−8 m/s2。基于海浪模拟平台的实测试验表明：满足约束条件时，这三类瞬时相位处理算法均适用于动态环境，并获得了标准差为4.6×10−5 m/s2的绝对重力测量值。更进一步，基于本文动态适应性结论对系泊和海面船载动态环境进行评估，结果表明：测量船在3级以下海况可以进行10−5 m/s2量级绝对重力的测量；在3—4级海况下，需根据振动信号对测量结果进行筛选以获得10−5 m/s2量级绝对重力测量结果；4级以上海况则不适宜进行动态绝对重力测量。Abstract: The dynamic adaptability of digital fringe processing is a fundamental prerequisite for dynamic absolute gravity measurements. We have developed a dynamic restriction based on physical parameters of the background vibration according to the instantaneous phase methods, i.e., the reverse vibration velocity cannot exceed the falling velocity of the falling drop. This method involves the Hilbert transform, the direct quadrature, and the zero-crossing method. We have conducted simulation experiments with different vibration intensities and frequencies. The results show that this dynamic restriction is reasonable. The precision of the Hilbert transform method is better than 10−13 m/s2, the precision of the zero-crossing method is better than 10−9 m/s2, and the precision of the direct quadrature method is approximately （−7.9±2.0）×10−8 m/s2 at a 60 MHz sampling rate when the dynamic restriction is satisfied. Furthermore, we conducted absolute gravity measurements on a wave simulator, and obtained results with a standard deviation of 4.6×10−5 m/s2 . This verifies that the instantaneous phase methods apply to the dynamic environment when the restriction is satisfied. Based on the dynamic adaptability conclusion, we further evaluated the mooring system and the sea surface shipboard dynamic environment on the vessel. The evaluation revealed that the measurement vessel can conduct absolute gravity measurements with an accuracy of 10−5 m/s2 under sea conditions of grade 3 or below. Under sea conditions of grade 3−4, the measurement results need to be screened based on vibration signals to obtain absolute gravity measurements at the 10−5 m/s2 level. However, conducting dynamic absolute gravity measurements is not recommended under sea conditions above grade 4.
-
引言
地球重力场的精确测量可涉及多个领域的研究，比如在大地测量学中，重力等势面被用作高度参考，进而描述大陆和海洋地形；在地球物理学中，重力数据提供了关于地下质量分布及变化的信息，用于研究地质构造、反映火山和地震运动、监测冰川融化和水体变化、勘探石油矿藏等；在惯性导航中，高精度的导航算法也依赖于精确的重力场模型（Sandwell et al，2014；刘敏等，2017a；孙和平等，2021）。
静态条件下的重力测量技术相对成熟，陆基绝对重力测量的精度可以优于5×10−8 m/s2 （高景龙，1993；Niebauer et al，1995；吴琼，2011；胡华等，2012）。动态条件下，即在船舶、飞机等运动载体上进行重力测量，精度为10−5 m/s2量级，具有测量覆盖效率高、适应性广等特点（Verdun et al，2002；Sokolov et al，2016；刘敏等，2017b）。然而，目前在运动载体上的重力测量均是相对重力测量，绝对重力测量与相对重力测量相比有以下优势：① 不需要在测量前后前往重力基点进行标定，可提高测量效率；② 因为每一次测量都相对独立（Bidel et al，2018），可以克服仪器零漂和误差积累问题。因此，开展动态条件下的绝对重力测量技术研究具有重要的实用意义（Brown et al，2001；Merlet et al，2010；Baumann et al，2012；吴彬等，2020；Bidel et al，2020；吴燕雄等，2022）。
然而，进行陆基绝对重力测量需考虑使用的落体自由下落干涉信号相位提取算法是否可用，即，激光干涉原理的绝对重力测量仪器有较大的振动直接作用在干涉测量的参考点上时，在静态条件下普遍使用的干涉信号相位提取算法（Svitlov et al，2012；Svitlov，Araya，2014），是否适用于较大振动，适用于何种程度的振动，是动态条件下干涉信号处理算法拓展应用的重要依据，也是动态绝对重力测量仪器设计的重要约束。
本文拟针对动态绝对重力测量环境，分析激光干涉绝对重力仪干涉信号瞬时相位提取算法的效果并评估其可适应的垂向振动的动态上限。在能够获取干涉信号高速数字采样信号的前提下，把干涉信号处理转换成瞬时相位求解，从可测量的背景振动物理量角度，提出了激光干涉绝对重力测量中干涉信号瞬时相位求解算法的动态约束条件。通过构建单频振动信号模拟实验，模拟在不同强度、不同频率振动信号下的绝对重力测量，以验证约束条件，并给出三种干涉信号瞬时相位提取算法的误差对比。通过海浪模拟平台试验验证动态环境下瞬时相位算法的适用性，应用动态约束条件对船载系泊和海面船载动态环境的振动进行评估，给出不同海况下，船载绝对重力测量作业的建议，以期为船载等动态环境下进行激光干涉绝对重力测量的数据处理提供理论支撑。
1. 基于瞬时相位的干涉信号处理算法
1.1 瞬时相位
激光干涉绝对重力仪的测量原理是基于变型设计的迈克尔逊干涉仪，完成真空环境中自由下落落体（角棱镜）相对于参考点（参考棱镜）时间位移的测量，通过最小二乘拟合求解重力加速度绝对值（滕云田等，2013）。具体是稳频激光器发射的激光束经过分光镜获得参考光束和测量光束，参考光束的光程固定，测量光束经自由下落的角棱镜、参考棱镜反射后，与测量光束形成干涉信号，被光电传感器接收。光电传感器的输出U可以描述为
$$ U＝{U}_{0}\mathrm{cos}\left(2{\pi }\frac{2\Delta {\textit{z}}}{\lambda }\right) \text{，} $$ (1) 式中：U0是干涉信号最大光强时对应的输出电压；λ为激光波长；∆z是角棱镜的下落位移，2∆z是角棱镜的下落位移引起光程2倍的变化。角棱镜在真空腔中自由下落5—20 cm，忽略重力梯度的影响，下落位移∆z是一小段抛物线和参考棱镜振动叠加的函数
$$ \mathrm{\Delta }{\textit{z}} （ t ） ＝{{\textit{z}}}_{0}＋{v}_{0}t＋\frac{1}{2}g{t}^{2}＋\mathrm{\delta }{\textit{z}} （ t ） \text{，} $$ (2) 式中：z0，v0和g分别为起始位置、初速度和重力加速度，δz（t）表示测量环境振动等引起的参考棱镜振动。若把位移随时间的变化调制到了三角函数的相位上，即把式（2）带入式（1），可得到一个变频的干涉信号
$$ U （ t ） ＝{U}_{0}\mathrm{cos}\varnothing （ t ） \text{，} $$ (3) 其中相位$\varnothing （ t ） $为：
$$ \varnothing （ t ） ＝\frac{4\pi }{\lambda \left[{{\textit{z}}}_{0}＋{v}_{0}t＋\dfrac{1}{2}g{t}^{2}＋\mathrm{\delta }{\textit{z}} （ t ） \right]} {\text{.}} $$ (4) 因此，干涉信号处理的目标就是解算出（ti，$\varnothing （ t_i ） $）序列。
干涉信号的相位是随时间变化的信号，可以把干涉信号的处理转化为瞬时相位（频率）求解问题。相关文献中多称为“瞬时频率”问题（Huang et al，2009），即把瞬时相位求解出来后，再对时间求导，通过$\omega （ t ） ＝{\mathrm{d}\varnothing （ t ） }/{\mathrm{d}t}$ 取得。而在干涉信号处理的应用场景中，只需要获得瞬时相位，不需要再对时间求导，为描述统一，在本文中统一称“瞬时相位”。
1.2 瞬时相位求解方法
本文分别选取了三种比较推荐的瞬时相位算法：希尔伯特变换、直接正交法和过零点法来求解干涉信号的瞬时相位。过零点法和希尔伯特变换并不是首次用于绝对重力测量的干涉信号处理，本文旨在瞬时相位理论框架下，对算法的动态适应性进行分析和评估。
三种算法的原理（Huang et al，2009）为：
1） 希尔伯特变换（Hilbert transform，缩写为HT）。对于变量x（t），其希尔伯特变换定义为：
$$ H ［ x （ t ） ］ ＝\hat{x} （ t ） ＝\frac{1}{\pi }{\int }_{-\infty }^{＋\infty }\frac{x （ \tau ） }{t－\tau } \text{，} $$ (5) 式中，$\tau $为积分变量。
实信号及其希尔伯特变换构成唯一确定的解析信号：
$$ {\textit{z}} （ t ） ＝x （ t ） ＋{\rm{i}}\hat{x} （ t ） ＝A （ t ） {{\rm{e}}}^{-{\rm{i}}\theta （ t ） } \text{，} $$ (6) 式中，A（t）为解析信号的幅值，θ（t）为解析信号的相位。解析信号的相位也被认为是信号的瞬时相位，可以通过式（7） 和相位解缠算法计算获得，即
$$ \theta （ t ） ＝\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{\hat{x} （ t ） }{x （ t ） } \text{，} $$ (7) 2） 直接正交法（direct quadrature，缩写为DQ）。对于式（3）的纯调频（相位）信号， 其正值正交信号为
$$ {U}_{0}\mathrm{sin} ［ \varnothing （ t ） ］ ＝\sqrt{1－{ ［ {U}_{0}\mathrm{c}\mathrm{o}\mathrm{s}\varnothing^2 （ t ） ］ }}＝\sqrt{1－{ U^{2} （ t ） }} \text{，} $$ (8) 这时可以有两种方法计算相位：一种是反余弦，但这种方法在局部极值点处误差较大；另一种是通过四象限反正切和相位解缠算法进行相位展开，即
$$ \varnothing （ t ） ＝\mathrm{sign}\cdot \mathrm{arctan}\frac{\sqrt{1-{U^{2} （ t ） }}}{U （ t ） } \text{，} $$ (9) 式中，sign是信号正交项的符号，可以通过三角函数的变化方向判断。在实际使用中，由于极值点附近的采样点稀疏，会造成波形畸变，通过三点或五点中值滤波器可以得到一定的改善。
3） 过零点法（zero crossing，缩写为ZC）。过零点算法是计算局部瞬时相位的最基础的一种方法，在高速数字采样设备普及前，通过硬件检测电路就可以实现。干涉信号的过零时刻，就是式（1）为零的时刻，此时瞬时相位为kπ，即
$$ \mathrm{\varnothing } （ t ） ＝2{\pi }\frac{2\Delta {\textit{z}} （ t ） }{\lambda }＝k\pi \text{，}\qquad k＝1\text{，}2，\cdots ，n ，$$ (10) 则可以解得
$$ \Delta {\textit{z}} （ t ） ＝\frac{k}{4}\lambda \text{，}\qquad k＝1\text{，}2，\cdots ， n {\text{.}}$$ (11) 1.3 瞬时相位求解算法的动态约束条件
Huang等（1998）在前人研究的基础上，提出瞬时相位对单组份信号才有意义。对于非线性非稳态过程，多组份信号要经过经验模态分解（empirical mode decomposition， 缩写为EMD）得到本征模态函数（intrinsic mode function， 缩写为IMF）。IMF需要满足两个条件：① 局部极大值以及局部极小值的数目之和必须与零交越点的数目相等或最多只能差1，也就是说一个极值后面必需马上接一个零交越点；② 在任何时间点，局部最大值所定义的上包络线与局部极小值所定义的下包络线，取平均需接近于零。
若使HT得到的解析信号瞬时相位有意义，不仅需要是单组份信号，还要满足Bedrosian定理（Bedrosian，1963），即包络和载波的傅里叶谱不重叠：
$$ H ［ a （ t ） \mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\varnothing } （ t ） ］ ＝a （ t ） H ［ \mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\varnothing } （ t ） ］ . $$ (12) 对于式（3）的干涉信号，幅值是常数，所以满足Bedrosian定理。
此外，还需要满足Nuttall定理（Nuttall，Bedrosian，1966）：对任何给定的函数$x （ t ） ＝a （ t ） \mathrm{cos}\mathrm{\varnothing } （ t ） $，对于任意a（t）和$ \mathrm{\varnothing } （ t ） $，如果x（t）的HT由xh（t）给出，x（t）的正交信号为xq（t），那么HT和正交信号相同，即${H}\left\{a （ t ） \mathrm{c}\mathrm{o}\mathrm{s} ［ \varnothing （ t ） ］ \right\}＝\mathrm{sin} ［ \varnothing （ t ） ］ $的充要条件是
$$ E＝{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{ ［ xh （ t ） －xq （ t ） ］ }^{2}\mathrm{d}t＝2{\int }_{-\mathrm{\infty }}^{{\omega }_{0}}{F}_{\mathrm{q}} （ \omega ） \mathrm{d}\omega ＝0 \text{，} $$ (13) 式中：Fq（ω）是信号的正交频谱，${F}_{\mathrm{q}} （ \omega ） ＝F （ \omega ） ＋\mathrm{i}{\displaystyle\int }_{-\infty }^{＋\infty }a （ t ） \mathrm{sin}\varnothing （ t ） {\mathrm{e}}^{-{\rm{i}}\omega t}\mathrm{d}t$；F（ω）是信号的频谱。Nuttall定理揭示了当信号正交频谱不够对称分布在零线上下时，HT和正交变换之间是有误差的。
以上研究说明瞬时相位算法是有约束条件的，并不是所有的干涉信号都可以通过瞬时相位算法求解出瞬时相位，如果不能满足约束条件，可能就无法获得相位的瞬时值。因此判断干涉信号的瞬时相位算法是否有效，对于绝对重力数据处理至关重要。
而IMF条件和Nuttall定理都是对信号形态的约束，即需要先得到干涉信号，再判断信号是否满足条件。而实际中更需要的是可以预判动态环境是否满足约束条件，再开展测量或采取针对性的处理。因此本文尝试将可测量的振动物理量，作为预判动态环境下瞬时相位算法有效的约束条件。
对于某些信号，希尔伯特算法求解出的瞬时频率可能为负值，而负频率无实际物理意义。激光干涉绝对重力仪干涉信号的瞬时相位由式（4）获得，对相位求导获得瞬时频率
$$ f （ t ） ＝\frac{1}{2\mathrm{\pi }}\frac{\mathrm{d}}{\mathrm{d}t} ［ \varnothing （ t ） ］ \text{，} $$ (14) 对振动位移求导，得到振动速度vvib（t）。因此有
$$ f （ t ） ＝\frac{2}{\lambda } ［ gt＋{v}_{0}＋{v}_{\mathrm{v}\mathrm{i}\mathrm{b}} （ t ） ］ {\text{.}} $$ (15) 只有瞬时频率非负才有物理意义，1.2节的瞬时相位求解算法才可以应用，即
$$ -{v}_{\mathrm{v}\mathrm{i}\mathrm{b}} （ t ） ＜ gt＋{v}_{0} {\text{.}}$$ (16) 因此，得到瞬时相位存在的约束条件为：反向的振动速度不能超过落体的下落速度。同一时刻干涉信号、瞬时频率的计算值与瞬时相位求解误差的对应情况如图1所示。当瞬时频率的计算值存在负值时，对应的干涉信号产生了明显畸变，且局部极值点多于过零点两个以上，信号没有均匀分布在x轴上下，故无法同时满足IMF条件和Nuttall定理。如果HT变换准确，那么式（6）的解析信号的虚部就应该是实部的正交信号，其模值恰好是信号幅值U0，所以瞬时相位求解的误差情况可以用HT变换得到的解析信号模值与信号幅值的差反应，即
![]()
图 1 干涉信号、瞬时相位求解误差和瞬时频率的计算值的对应情况Figure 1. The fringe signal and the phase’s instantaneous value error compared with the calculated instantaneous frequency$$ \mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}＝\frac{\left|{\textit{z}} （ t ） \right|－{U}_{0}}{{U}_{0}} \text{，} $$ (17) 由图1还可以看出，HT变换的误差在畸变附近区间也明显增大。
如果把振动简化为一个单一频率的信号进行分析，即式（4）中的振动波动量
$$ {\delta }{\textit{z}} （ t ） ＝\varepsilon \mathrm{cos} （ \omega t＋\theta ） \text{，} $$ (18) 式中，ε，ω和θ依次为幅值、频率和初始相位。则振动速度为
$$ {v}_{\mathrm{v}\mathrm{i}\mathrm{b}} （ t ） ＝\frac{\mathrm{d}\delta {\textit{z}} （ t ） }{\mathrm{d}t}＝-\varepsilon \omega \mathrm{sin} （ \omega t＋\theta ） \text{，} $$ (19) 代入约束条件式（16），可得
$$ \varepsilon \omega \mathrm{sin} （ \omega t＋\theta ） {\text{≤}} gt＋{v}_{0} \text{，} $$ (20) 把式（20）进行适当简化可知，当满足约束条件
$$ \varepsilon \omega {\text{≤}} {v}_{0} \text{，} $$ (21) 则瞬时频率为正，可以用2.2节的方法求解干涉信号的瞬时相位。εω可以认为是速度信号的幅值，进一步εω2是加速度信号的幅值。图2a−c给出了振动信号速度幅值满足约束条件时，瞬时频率恒正，通过HT方法得到的瞬时相位高精度地逼近模拟输入的瞬时相位（固定相位差）；图2d−f显示出不满足约束条件时，瞬时频率存在负值，HT方法得到的瞬时相位与模拟输入的瞬时相位有明显的偏差。
![]()
图 2 满足（左）和不满足（右）εω≤v0条件的振动信号对瞬时相位求解的影响（a） 振动信号；（b） 计算出的干涉信号的瞬时频率；（c） 通过HT方法求解的瞬时相位与干涉信号的实际相位对比Figure 2. The effect of whether the vibration signal satisfies the restrictions（a）Vibration signal ；（b） Instantaneous frequency of the calculated interference signal；（c） Comparison of the instantaneous phase solved by the HT method with the actual phase of the interference signal1.4 瞬时相位算法误差分析
在满足动态约束条件的情况下，DQ方法的误差主要源于采样点稀疏，以致极值点未被采样到。如图3a所示，通过DQ方法得到的正交信号与仿真信号直接变换的所得理想正交信号的偏差同干涉信号的对应关系可以看出，DQ法计算的正交信号在峰值处的误差较大。图3b给出了单次下落中，用DQ算法解算出的位移与真值之间的误差，且误差随时间增大而增大。因为在下落中，落体速度越来越快，对应干涉信号的瞬时频率不断升高，同一个完整干涉条纹的采样点数目就不断减少，导致了瞬时相位求解误差的不断增大。
![]()
图 3 瞬时相位算法的误差分析（a） DQ算法正交信号和仿真信号的正交信号之差与干涉信号的对比关系；（b—d）分别为单次下落中，DQ （b），HT （c）和ZC （d）算法位移与仿真信号的位移之差随时间的变化Figure 3. Error analysis of the instantaneous phase methods（a） the difference between the quadrature signal obtained by the DQ method and the ideal quadrature signal comparing with the fringe signal；（b）−（d） The deviation of the DQ， HT， ZC calculated displacement from the actual valueHT方法的误差主要来自吉布斯现象造成的端点效应，如图3c所示，信号两端的误差明显高于中间信号。因此，在HT方法中，需对信号进行适当截取，舍弃掉误差最大的两端各1%的数据。
ZC方法的误差主要来自于对过零点的估计，如图3d所示，误差随下落时间增大而增大。与DQ法一样，干涉信号的频率增加，相对的采样点越稀疏，对过零点的估计误差也随之增大。
提高DQ和ZC算法精度最直接有效的方式是提高采样频率。如图4所示，以不同的采样频率进行实验，统计符合动态约束条件的实验结果，可以看出，DQ法和ZC法的误差，随着采样频率的增大而显著减小，再一次验证了DQ法和ZC法的误差与采样点的稀疏程度有关。但HT法截取掉两端的信号后，提高采样频率对误差的影响并不明显。
![]()
图 4 HT （a），DQ （b）和ZC （c）算法在不同采样频率下g的误差Figure 4. g errors of HT （a），DQ （b） and ZC （c） methods at different sampling rates2. 验证实验
2.1 动态约束条件仿真验证实验
为了验证动态约束条件 ［ 式（21） ］ ，构造不同频率ω和幅值ε的单频振动信号以进行仿真实验，同时对比了干涉信号瞬时相位处理算法在满足和不满足约束条件时的精度。
第一步，构造干涉信号，标准g0值取9.801 103 43 m/s2，落体下落的测量时间段为30—130 ms，采样频率为60 MHz，自由落体下落位移叠加单频振动信号为
$$ V （ t ） ＝P*\mathrm{cos} （ 2\pi {f}_{j}t ） ，$$ (22) 以模拟重力仪的机械干扰和主要环境振动干扰。船上200 Hz以上的干扰能量较小，εω不会超过约束条件；频率低于0.01 Hz，幅值超过50 m才可能超过约束条件，因此构造船载可能会出现超过约束条件的振动频率范围，$ {f}_{j} $ 取0.01—200 Hz。通过调节比例因子P，把振动信号带入式（1），（2），生成干涉信号。
第二步，分别通过HT，DQ和ZC方法，求解瞬时相位 ［ ti，$\varnothing $（ti） ］ 序列，并根据式（1）把瞬时相位转化成落体的下落位移 ［ ti，∆z（ti） ］ 轨迹。
第三步，补偿振动，为突出动态条件对瞬时相位提取算法的约束，计算结果中需要根据已知的振动信号进行振动干扰修正，即从包含振动的位移轨迹里，减去加入的振动信号，即式（22）。因为垂向振动本身就会引起测量误差，而在仿真实验中，振动干扰信号已知，补偿掉振动以后，存在的误差就是瞬时相位算法本身引入的解算误差。 对于HT和DQ方法，时间坐标是对齐的，因此可以从下落位移中直接减去振动干扰V（ti）；对于ZC方法，由于下落位移轨迹是等距离间隔流，而振动是等时间间隔流，需先对振动数据进行样条插值，待时间对齐后再做减法。
第四步，对去补偿振动以后的位移轨迹进行最小二乘法多项式拟合，得到最优g。比较g与g0得到三种瞬时相位算法的引入误差。
单频振动干扰信号实验中，振动信号的幅值、频率以及算法误差的关系如图5所示。根据实验参数，落体速度是从0.29 m/s变化到1.3 m/s。图5a表示幅值-频率平面根据归一化的εω量级p进行划分， 即
![]()
图 5 单频振动信号实验（a） 幅值和频率的乘积εω的归一化量级；（b—d） HT，DQ和ZC方法得到的重力值误差与振动信号的εω的对应关系Figure 5. Single frequency vibration signal experiments（a） εω，the product of amplitude and frequency of the vibration signal；（b−d） Errors of HT，DQ and ZC methods related to the εω of the vibration signal$$ p＝\mathrm{l}\mathrm{g}\left(\frac{\varepsilon \omega }{{v}_{0}} \right)\text{，} $$ (23) 式中，v0＝0.29 m/s是落体测量时的初始速度。当p为负值，即$ \varepsilon \omega ＜ {v}_{0} $，确定满足约束条件；p略大于0，即$ \varepsilon $ω≥v0，因为三角函数的存在，也有可能满足约束条件，可以被看作是过渡带；过渡带以上，p明显大于0时，εω不再满足约束条件，HT，DQ，ZC三种瞬时相位算法的精度都迅速降低。不同振动强度和频率下，HT，DQ和ZC算法得到的重力加速度g值误差如图5b−c所示，当εω没有超过约束条件时，同种瞬时相位算法本身的精度基本一致，可见瞬时相位算法适用于叠加不同振动信号的求解。精度对比结果列于表1。HT方法的误差优于10−13 m/s2；ZC方法的误差优于10−9 m/s2；DQ方法的误差为（−7.9±2.0）×10−8 m/s2。因为在激光绝对重力测量中其它因素引起的误差也在10−8 m/s2量级，因此小于10−9 m/s2的误差可以忽略不计。HT和ZC算法的误差都可以忽略不计，HT的误差更小。当振动信号速度幅值εω明显超过约束条件时，HT，DQ和ZC的误差可达1 m/s2以上，瞬时相位算法不再适用。
表 1 不同方法的瞬时相位误差Table 1. Errors in the instantaneous phase methods方法 g－g0 误差/ m·s−2 εω≤v0 εω$\gg $v0 HT （−5.9±3.3）×10−14 绝对值＞1 DQ （−7.9±2.0）×10−8 绝对值＞1 ZC （2.4±1.0）×10−10 绝对值＞1 2.2 海浪模拟平台实测试验
为了进一步验证满足动态约束条件下，瞬时相位算法在动态环境下的适用性，我们开展了海浪模拟平台上的实测试验，如图6a所示。设计的船载绝对重力测量系统由绝对重力测量部分、惯性稳定平台、加速度计和GPS组成。其中惯性稳定平台主要抑制纵摇和横摇带来的角度变化干扰。垂直方向上，GPS信号采样率低于10 Hz，用于测量仪器海拔高度的变化，补偿重力梯度变化。加速度计采集参考棱镜的垂向振动，与绝对重力测量系统进行同步采集，主要用于垂向振动的补偿（吴燕雄等，2022）。但干涉信号的瞬时相位准确求解是各种补偿修正处理的前提。
![]()
图 6 海浪模拟平台试验装置（a）和模拟平台运动状态（b）（周期为17 s）Figure 6. Sea wave simulator devices （a） and wave simulator movement （b） （the period is 17 s）图6b给出了海浪模拟平台三个方向的运动加速度，信号周期17 s，垂向加速度略小于0.01 m/s2，模拟1—2级平静海况。估算出的速度信号幅值εω＜0.03 m/s，满足1.3节讨论的动态适用条件，所以瞬时相位算法适用。
因HT本身的算法精度最高，同时数据点与加速度计采样点一致，易于实现振动抑制算法，因此基于HT算法进行了振动抑制处理，结果如图7所示。每分钟测量1次，连续进行了约1.5小时的测量，用3δ准则排除外点后得到80个测量结果（肖凡，2012），测量的处理结果满足正态分布。以每组10个分组，其组标准差为4.6×10−5 m/s2，有效验证了满足动态约束条件时，瞬时相位算法可以适用于动态环境下的绝对重力测量。
![]()
图 7 周期为17 s的海浪模拟平台试验结果（a） 80个测量结果的直方图分布；（b） 以每组10个进行的分组测量的结果与均值间的偏差Figure 7. Results of 17 s period wave simulator experiment（a） Histogram distribution of 80 measurements，satisfies normal distribution；（b） Deviation of group measurements from the mean value for 10 units/group3. 船载动态环境下的应用仿真
由于有动态约束条件的限制，因此在开展动态环境下的激光绝对重力测量前，有必要对测量环境的振动参数进行评估。实际的振动干扰信号可以看成是若干频率信号的叠加，通过考查幅值和频率乘积较大的频率分量，来预估瞬时相位算法的有效性。
3.1 船载系泊状态
为了开展船载绝对重力测量，在测量船上安装了加速度计，采集了系泊船上甲板的背景振动，振动主要来自于海浪的影响和发动机的振动，如图8所示。评估振动信号，频率主要集中在0.1—1 Hz和50—100 Hz频率段，加速度幅值εω2不超过10−2 m/s2，那么速度的幅值εω小于0.1 m/s。根据约束条件预估，在船载系泊条件下瞬时相位算法是适用的。
![]()
图 8 船系泊时船体垂向加速度（a）及加速度信号功率谱（b）Figure 8. The vertical acceleration （a） of the vessel during mooring，and the power spectrum （b） of the acceleration signal与单一频率振动实验类似的，构造仿真实验验证。选取连续50分钟的振动数据，每分钟进行2次模拟下落测量，共计100次下落。分别用HT，DQ和ZC方法，并基于已知的振动信号进行振动误差补偿，得到的重力加速度g值与真值g0的偏差和100个结果的标准差列于表2。船载系泊振动信号的实验验证了预判的结果，在未超过振动约束条件下，三种瞬时相位算法均有效。这也是激光绝对重力仪能够进行船载系泊状态下测量的前提条件，其中HT方法的算法精度最高，其引入误差为（8.2×10−15±3.1×10−13） m/s2，可以忽略。
表 2 船系泊状态下g值的偏差和标准差Table 2. The bias of the average g and the standard deviations of g at vessel mooring方法 g与g0 偏差/（m·s−2） 标准差/（m·s−2） HT −8.2×10−15 3.1×10−13 DQ −8.4×10−8 2.7×10−7 ZC 1.5×10−10 1.2×10−7 3.2 海面船载动态环境
通过船载的惯性传感器记录到了测量船在海面航行时的垂荡（垂直运动），分别选取1—2级海况（时段1）和3—4级海况（时段2）进行分析。海浪频率主要在0.1—1 Hz，根据动态约束条件和图5的实验结论，时段1垂荡幅值小于0.5 m，满足动态约束条件；时段2垂荡最大幅值超过1.5 m，有超过动态约束条件的情况，因此预判时段2在部分时刻瞬相位方法不适用。海浪的波动周期远大于落体下落的测量时间时，在海浪波动的局部，仍可满足约束条件的，部分时刻瞬时相位算法仍适用。
![]()
图 9 船航行在时段1 （a）和时段2 （b）的垂荡位移Figure 9. Raw ship heave data on sailing at the period 1 （a） and period 2 （b）与系泊状态实验相似，我们用不同海况下的船载真实振动数据构造绝对重力仪干涉信号，并模拟每分钟2次下落进行测量，连续50分钟，共计100次下落。分别用HT，DQ和ZC方法，基于同样已知的振动信号对测量结果进行振动误差修正，得到的重力加速度g值与真值g0的偏差及100个结果的标准差如表3所示。
表 3 船航行状态下g值的偏差和标准差（单位：10−8 m·s−2）Table 3. The bias of the average g and the standard deviations of g on ship sailing state （unit：10−8 m·s−2）HT方法 DQ方法 ZC方法 g－g0 误差 标准差 g－g0误差 标准差 g－g0 误差 标准差 时段1 −9.0×10−7 3.5×10−6 −8.4 0.24 1.1×10−2 9.4×10−2 时段2（全） ＞ 108 ＞ 108 ＞ 108 时段2（筛选） −1.5×10−6 4.2×10−6 −8.4 0.27 1.8×10−2 0.18 因为时段2存在不满足动态约束条件的时刻，所以有些瞬时相位的解算是可靠的，有些结果误差较大，如图10所示。根据同步振动速度的大小进行筛选，排除掉不满足动态约束条件的18次测量结果后，采用HT，DQ和ZC方法的最终结果也列于表3。从表3可以看出，在海况比较平稳（如1—2级海况）的情况下，可以进行船载绝对重力的连续测量，而在海况波动较大（如3—4级海况）下，则必须根据同步记录的振动信号的速度大小进行筛选，选取满足振动约束条件的时刻，以获得可靠测量结果。如果海况波动加重，可靠测量值的比例就会进一步降低，所以4级以上海况不适合进行船载绝对重力测量作业。
![]()
图 10 3—4级海况下连续100次测量的HT，DQ和ZC瞬时相位算法误差及船的垂荡速度Figure 10. Errors of HT，DQ and ZC methods for 100 consecutive measurements at level 3−4 sea state and the vertical velocity of the vessel4. 结论
在动态环境下开展激光干涉原理的绝对重力测量时，必须考虑是否满足干涉信号瞬时相位算法的动态约束条件。本文通过理论分析，给出了基于背景振动物理量的定量动态约束条件形式，可用于预判动态环境能否进行绝对重力测量。仿真实验计算结果表明，此种形式的动态约束条件是合理可行的，在满足动态约束条件的振动下，HT和ZC算法本身的误差均优于0.1×10−8 m/s2，适用于动态环境下绝对重力测量的数据处理。
在海浪模拟平台的实测试验中，基于课题组开发的海洋绝对重力测量系统，利用约束条件开展了动态环境下的绝对重力实验性测试，结果表明，瞬时相位处理算法适用于动态环境，重力测量值标准差为4.6×10−5 m/s2。
结合船载环境的动态特性，本文分别进行了系泊和海面船载动态环境下，瞬时相位算法适用性评价，预判在系泊状态和1—2级海况下，背景振动满足动态约束条件，瞬时相位算法适用；在3—4级海况下，背景振动不满足动态约束条件，仿真实验结果也显示，部分测量结果不可靠，需要根据同步振动速度的大小对测量结果进行筛选；并且不建议在4级以上的海况进行船载绝对重力测量作业。
-
![]()
图 1 干涉信号、瞬时相位求解误差和瞬时频率的计算值的对应情况
Figure 1. The fringe signal and the phase’s instantaneous value error compared with the calculated instantaneous frequency
![]()
图 2 满足（左）和不满足（右）εω≤v0条件的振动信号对瞬时相位求解的影响
（a） 振动信号；（b） 计算出的干涉信号的瞬时频率；（c） 通过HT方法求解的瞬时相位与干涉信号的实际相位对比
Figure 2. The effect of whether the vibration signal satisfies the restrictions
（a）Vibration signal ；（b） Instantaneous frequency of the calculated interference signal；（c） Comparison of the instantaneous phase solved by the HT method with the actual phase of the interference signal
![]()
图 3 瞬时相位算法的误差分析
（a） DQ算法正交信号和仿真信号的正交信号之差与干涉信号的对比关系；（b—d）分别为单次下落中，DQ （b），HT （c）和ZC （d）算法位移与仿真信号的位移之差随时间的变化
Figure 3. Error analysis of the instantaneous phase methods
（a） the difference between the quadrature signal obtained by the DQ method and the ideal quadrature signal comparing with the fringe signal；（b）−（d） The deviation of the DQ， HT， ZC calculated displacement from the actual value
![]()
图 4 HT （a），DQ （b）和ZC （c）算法在不同采样频率下g的误差
Figure 4. g errors of HT （a），DQ （b） and ZC （c） methods at different sampling rates
![]()
图 5 单频振动信号实验
（a） 幅值和频率的乘积εω的归一化量级；（b—d） HT，DQ和ZC方法得到的重力值误差与振动信号的εω的对应关系
Figure 5. Single frequency vibration signal experiments
（a） εω，the product of amplitude and frequency of the vibration signal；（b−d） Errors of HT，DQ and ZC methods related to the εω of the vibration signal
![]()
图 6 海浪模拟平台试验装置（a）和模拟平台运动状态（b）（周期为17 s）
Figure 6. Sea wave simulator devices （a） and wave simulator movement （b） （the period is 17 s）
![]()
图 7 周期为17 s的海浪模拟平台试验结果
（a） 80个测量结果的直方图分布；（b） 以每组10个进行的分组测量的结果与均值间的偏差
Figure 7. Results of 17 s period wave simulator experiment
（a） Histogram distribution of 80 measurements，satisfies normal distribution；（b） Deviation of group measurements from the mean value for 10 units/group
![]()
图 8 船系泊时船体垂向加速度（a）及加速度信号功率谱（b）
Figure 8. The vertical acceleration （a） of the vessel during mooring，and the power spectrum （b） of the acceleration signal
![]()
图 9 船航行在时段1 （a）和时段2 （b）的垂荡位移
Figure 9. Raw ship heave data on sailing at the period 1 （a） and period 2 （b）
![]()
图 10 3—4级海况下连续100次测量的HT，DQ和ZC瞬时相位算法误差及船的垂荡速度
Figure 10. Errors of HT，DQ and ZC methods for 100 consecutive measurements at level 3−4 sea state and the vertical velocity of the vessel
表 1 不同方法的瞬时相位误差
Table 1 Errors in the instantaneous phase methods
方法 g－g0 误差/ m·s−2 εω≤v0 εω$\gg $v0 HT （−5.9±3.3）×10−14 绝对值＞1 DQ （−7.9±2.0）×10−8 绝对值＞1 ZC （2.4±1.0）×10−10 绝对值＞1 
下载: 导出CSV
表 2 船系泊状态下g值的偏差和标准差
Table 2 The bias of the average g and the standard deviations of g at vessel mooring
方法 g与g0 偏差/（m·s−2） 标准差/（m·s−2） HT −8.2×10−15 3.1×10−13 DQ −8.4×10−8 2.7×10−7 ZC 1.5×10−10 1.2×10−7
下载: 导出CSV
表 3 船航行状态下g值的偏差和标准差（单位：10−8 m·s−2）
Table 3 The bias of the average g and the standard deviations of g on ship sailing state （unit：10−8 m·s−2）
HT方法 DQ方法 ZC方法 g－g0 误差 标准差 g－g0误差 标准差 g－g0 误差 标准差 时段1 −9.0×10−7 3.5×10−6 −8.4 0.24 1.1×10−2 9.4×10−2 时段2（全） ＞ 108 ＞ 108 ＞ 108 时段2（筛选） −1.5×10−6 4.2×10−6 −8.4 0.27 1.8×10−2 0.18
下载: 导出CSV
-
高景龙. 1993. NIM-3型新的轻小高精度可移式绝对重力仪[J]. 测绘学报,22(3):223–229. Gao J L. 1993. A new type NIM-3 transportable absolute gravimeter of small size and light weight[J]. Acta Geodaetica et Cartographica Sinica,22(3):223–229 (in Chinese).
胡华,伍康,申磊,李刚,王力军. 2012. 新型高精度绝对重力仪[J]. 物理学报,61(9):099101. Hu H,Wu K,Shen L,Li G,Wang L J. 2012. A new high precision absolute gravimeter[J]. Acta Physica Sinica,61(9):099101 (in Chinese). doi: 10.7498/aps.61.099101
刘敏,黄谟涛,欧阳永忠,邓凯亮,翟国君,陆秀平,吴太旗,陈欣. 2017a. 海空重力测量及应用技术研究进展与展望(一):目的意义与技术体系[J]. 海洋测绘,37(2):1–5. Liu M,Huang M T,Ouyang Y Z,Deng K L,Zhai G J,Lu X P,Wu T Q,Chen X. 2017a. Development and prospect of air-sea gravity survey and its applications,part Ⅰ:Objective,significance and technical system[J]. Hydrographic Surveying and Charting,37(2):1–5 (in Chinese).
刘敏,黄谟涛,欧阳永忠,邓凯亮,翟国君,陆秀平,吴太旗,陈欣. 2017b. 海空重力测量及应用技术研究进展与展望(二):传感器与测量规划设计技术[J]. 海洋测绘,37(3):1–11. Liu M,Huang M T,Ouyang Y Z,Deng K L,Zhai G J,Lu X P,Wu T Q,Chen X. 2017b. Development and prospect of air-sea gravity survey and its applications,part Ⅱ:Sensor,plan and design of survey[J]. Hydrographic Surveying and Charting,37(3):1–11 (in Chinese).
孙和平,孙文科,申文斌,申重阳,祝意青,付广裕,吴书清,崔小明,陈晓东. 2021. 地球重力场及其地学应用研究进展:2020中国地球科学联合学术年会专题综述[J]. 地球科学进展,36(5):445–460. Sun H P,Sun W K,Shen W B,Shen C Y,Zhu Y Q,Fu G Y,Wu S Q,Cui X M,Chen X D. 2021. Research progress of Earth's gravity field and its application in geosciences:A summary of Annual Meeting of Chinese Geoscience Union in 2020[J]. Advances in Earth Science,36(5):445–460 (in Chinese).
滕云田,吴琼,郭有光,张兵,张涛. 2013. 基于激光干涉的新型高精度绝对重力仪[J]. 地球物理学进展,28(4):2141–2147. doi: 10.6038/pg20130459 Teng Y T,Wu Q,Guo Y G,Zhang B,Zhang T. 2013. New type of high-precision absolute gravimeter base on laser interference[J]. Progress in Geophysics,28(4):2141–2147 (in Chinese).
吴彬,周寅,程冰,朱栋,王凯楠,朱欣欣,陈佩军,翁堪兴,杨秋海,林佳宏,张凯军,王河林,林强. 2020. 基于原子重力仪的车载静态绝对重力测量[J]. 物理学报,69(6):060302. doi: 10.7498/aps.69.20191765 Wu B,Zhou Y,Cheng B,Zhu D,Wang K N,Zhu X X,Chen P J,Weng K X,Yang Q H,Lin J H,Zhang K J,Wang H L,Lin Q. 2020. Static measurement of absolute gravity in truck based on atomic gravimeter[J]. Acta Physica Sinica,69(6):060302 (in Chinese).
吴琼. 2011. 高精度绝对重力仪关键技术研究[D]. 北京: 中国地震局地球物理研究所: 82–92. Wu Q. 2011. The Study of the Key Technology in High-Precision Absolute Gravimeter[D]. Beijing: Institute of Geophysics, China Earthquake Administration: 82–92 （in Chinese）.
吴燕雄,滕云田,吴琼,徐行,张兵. 2022. 船载绝对重力仪测量系统的误差修正模型及不确定度分析[J]. 武汉大学学报•信息科学版,47(4):492–500. Wu Y X,Teng Y T,Wu Q,Xu X,Zhang B. 2022. Error correction model and uncertainty analysis of the shipborne absolute gravity measurement system[J]. Geomatics and Information Science of Wuhan University,47(4):492–500 (in Chinese).
肖凡. 2012. FG5高精度绝对重力测量影响因素研究[D]. 郑州: 中国人民解放军信息工程大学: 25–28. Xiao F. 2012. Research on the Effect Factors for the High Precision Absolute Gravity Surveying Using FG5[D]. Zhengzhou: PLA Information Engineering University: 25–28 （in Chinese）.
Baumann H,Klingelé E E,Marson I. 2012. Absolute airborne gravimetry:A feasibility study[J]. Geophys Prospect,60(2):361–372. doi: 10.1111/j.1365-2478.2011.00987.x
Bedrosian E. 1963. A product theorem for Hilbert transforms[J]. Proc IEEE,51(5):868–869.
Bidel Y,Zahzam N,Blanchard C,Bonnin A,Cadoret M,Bresson A,Rouxel D,Lequentrec-Lalancette M F. 2018. Absolute marine gravimetry with matter-wave interferometry[J]. Nat Commun,9(1):627. doi: 10.1038/s41467-018-03040-2
Bidel Y,Zahzam N,Bresson A,Blanchard C,Cadoret M,Olesen A V,Forsberg R. 2020. Absolute airborne gravimetry with a cold atom sensor[J]. J Geod,94:20. doi: 10.1007/s00190-020-01350-2
Brown J M, Niebauer T M, Klingele E. 2001. Towards a dynamic absolute gravity system[C]//Gravity, Geoid and Geodynamics 2000. Banff, Alberta, Canada: Springer: 223–228.
Huang N E,Shen Z,Long S R,Wu M C,Shih H H,Zheng Q A,Yen N C,Tung C C,Liu H H. 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proc Roy Soc A:Math Phys Eng Sci,454(1971):903–995. doi: 10.1098/rspa.1998.0193
Huang N E,Wu Z H,Long S R,Arnold K C,Chen X Y,Blank K. 2009. On instantaneous frequency[J]. Adv Adapt Data Anal,1(2):177–229. doi: 10.1142/S1793536909000096
Merlet S, Le Gouët J, Bodart Q, Clairon A, Landragin A, Dos Santos F P, Rouchon P. 2010. Vibration rejection on atomic gravimeter signal using a seismometer[C]//Gravity, Geoid and Earth Observation. Chania, Crete, Greece: Springer: 115–121.
Niebauer T M,Sasagawa G S,Faller J E,Hilt R,Klopping F. 1995. A new generation of absolute gravimeters[J]. Metrologia,32(3):159–180. doi: 10.1088/0026-1394/32/3/004
Nuttall A H,Bedrosian E. 1966. On the quadrature approximation to the Hilbert transform of modulated signals[J]. Proc IEEE,54(10):1458–1459. doi: 10.1109/PROC.1966.5138
Sandwell D T,Müller R D,Smith W H F,Garcia E,Francis R. 2014. New global marine gravity model from CryoSat-2 and Jason-1 reveals buried tectonic structure[J]. Science,346(6205):65–67. doi: 10.1126/science.1258213
Sokolov A V,Krasnov A A,Konovalov A B. 2016. Measurements of the acceleration of gravity on board of various kinds of aircraft[J]. Meas Tech,59(6):565–570. doi: 10.1007/s11018-016-1009-y
Svitlov S,Rothleitner C,Wang L J. 2012. Accuracy assessment of the two-sample zero-crossing detection in a sinusoidal signal[J]. Metrologia,49(4):413–424. doi: 10.1088/0026-1394/49/4/413
Svitlov S,Araya A. 2014. Homodyne interferometry with quadrature fringe detection for absolute gravimeter[J]. Appl Opt,53(16):3548–3555. doi: 10.1364/AO.53.003548
Verdun J,Bayer R,Klingelé E E,Cocard M,Geiger A,Halliday M E. 2002. Airborne gravity measurements over mountainous areas by using a LaCoste & Romberg air-sea gravity meter[J]. Geophysics,67(3):807–816. doi: 10.1190/1.1484525
计量
- 文章访问数: 125
- HTML全文浏览量: 67
- PDF下载量: 36
