1 引言

世界时(UT1)是以地球自转运动为基础的时间计量系统,目前被广泛应用于大地测量和天文导航等相关领域,它是实现天球与地球坐标系转换的必要联系参数之一。UT1反映了地球自转角速度的变化,它是地球定向参数(EOP)之一, UT1是EOP中变化最大、最快且最不易精准预测的量,每1 s内的变化都是不可忽略的^[1-2]。现阶段,人们可以通过国际地球自转服务(IERS)获取定期发布的UT1计算结果,其获取的UT1计算结果的时间分辨率为1 d。随着现代空间导航和深空探测等技术的不断发展,对UT1的及时性要求也越来越高,有的应用甚至需要实时的UT1计算结果,而通过IERS获取的UT1参数有时无法满足应用需求。

与地球固连的大型光纤陀螺仪(FOG)可以实现对地球自转运动的实时监测^[3-4],利用陀螺仪实时输出的地球自转角速度测量值可以用于解算UT1参数,该测量方法可以有效提高UT1参数获取的时间分辨率。本文从光纤陀螺仪的基本原理出发,研究一种利用光纤陀螺仪测量数据解算UT1的方法,并结合实测数据对计算方法进行验证,并与IERS提供的UT1参数结果进行对比,论证基于大型光纤陀螺仪测量解算高时间分辨率UT1参数的可行性。

2 基本原理

2.1 光纤陀螺仪的观测方程

利用光纤陀螺仪测量地球自转角速度的基本原理是基于Sagnac效应,它是法国物理学家Sagnac在20世纪初研究以太与地球之间的相互关系时发现的。光纤陀螺仪产生的Sagnac效应描述的是两束相同的光在同一光纤环路中沿相反方向传播,当光纤环路所在平面内存在旋转角速度时,两束光经过光纤环路传播后产生的相位差值与旋转角速度之间的转换关系^[5-6]可表示为

Δϕ=Kn·Ω,K=8πSλc,(1)

式中:Δϕ为Sagnac相位差;K为光纤陀螺仪的标度系数,其值与光纤环路所包围的面积成正比,与光源的平均波长λ和光速c之积成反比;n为光纤环平面的法向量,反映光纤环平面法线的瞬时指向;Ω为旋转角速度矢量。(1)式为理想状态下的Sagnac效应表达式,在实际测量过程中,需要考虑测量系统存在的噪声及各项误差,可以将(1)式改写为包含地球自转参数、陀螺仪定向变化以及陀螺仪仪器噪声的关系式^[7-8],即

Δϕ=Kn·Ω+Δϕinstr=Kn·(Ω0+ΔΩ)+Δϕinstr,(2)

式中: Δϕ_instr表示光纤陀螺仪的固有噪声引起的等效相位差变化量,该项噪声主要来源于背向散射、温度变化以及数据采集过程中光电探测器产生的噪声等;Ω₀为地球的平均自转角速度矢量;ΔΩ为地球自转角速度矢量的变化量,该变化包括自转轴的变化和自转角速率的变化。在大地坐标系下,地球的瞬时自转角速度矢量可以表示为

Ω=Ω0(m1,m2,1+m3)T,(3)

式中:Ω₀为地球平均自转角速率;m₁,m₂表示地球自转轴的瞬时位置,m₃为地球自转角速率的瞬时变化。

对于水平安放于地面的光纤陀螺仪,其光纤环平面的法向量可以表示为

n=(cosαcosβ,cosαsinβ,sinα)。(4)

图1所示为观测方程展开过程中涉及的各方向角的示意图,Ω为地球瞬时自转角速率(时间t的函数),(α,β)为向量n在大地坐标系下的两个指向角,其值与测站所在的大地经纬度(α₀,β₀)和陀螺仪测量平台在局部地壳形变的作用下产生的实际倾斜角度有关^[9-10]。将(3)式、(4)式代入(2)式中,则Sagnac相位差可表示为

Δϕ=K'[sinα0+m1cosα0cosβ0+m2cosα0sinβ0+m3sinα0+dTnscosα0]+Δϕinstr,(5)

式中:K'=(8πS/λc)Ω₀=KΩ₀,Ω₀=7.2921150×10^-5 rad/s;dT_ns表示n在南北方向产生的倾斜变化。(5)式即为光纤陀螺仪的观测方程。

图 1. 观测方程展开过程中各方向夹角示意图

Fig. 1. Schematic of each direction angle in expansion of observation equation

下载图片 查看所有图片

2.2 光纤陀螺仪的定向误差改正

在测量地球自转的过程中,光纤环平面法向量n发生偏移会导致光纤陀螺仪测量产生附加相位差,从而使输出的地球自转角速度值存在误差。通常情况下,法向量n的偏移是由测站地壳几何形变引起的,其主要原因是测站周围地理环境变化引起的地壳几何形变以及潮汐形变等^[11-12]。因此,测量过程中需要利用多台水平仪实时监测平台的倾斜变化,以便于对光纤陀螺仪进行定向改正。

利用水平仪可以完整测量出地壳几何形变引起的倾斜变化,但是除了测量到地壳几何形变外,水平仪还可以同时测量到地方铅垂线的变化,该变化主要由潮汐效应引起。测量过程中,光纤陀螺仪只对地壳的几何形变敏感,对测站铅垂线的变化不敏感。依据固体潮理论模型,测站在固体潮影响下产生的铅垂线变化可以表示为^[13-14]

TVD=-(1+k2-l2)∂W2(m)ag0∂θ,(6)

式中:k₂,l₂为2阶勒夫数; W2(m)为2阶m(m=0,1,2)级外部引潮力位;a为地球赤道平均半径;g₀为赤道平均重力加速度。

通过以上分析可知,单纯利用水平仪测量数据对光纤陀螺仪进行定向误差改正在消除地壳几何形变产生的影响时,会引入铅垂线变化的影响。因此,可以先利用水平仪测量结果对光纤陀螺仪进行倾斜改正,再通过(6)式利用固体潮理论模型模拟铅垂线变化结果对水平仪改正后结果进行修正,以补偿引入的铅垂线变化的影响,最终完成对光纤陀螺仪的定向误差改正。

2.3 角速率与UT1的转换关系

UT1是以地球自转运动为基础的时间,它与地球自转速率直接相关^[15-16],具体可以表示为

tUT1(t)=tUT1(t0)+1Ω0∫t0tΩ(t)dt,(7)

式中:t_UT1(t)表示协调世界时(UTC)的t时刻对应的UT1值;t与t₀分别表示时刻;Ω(t)为地球瞬时自转角速率;Ω₀表示地球平均自转角速率,所对应的转动周期为1恒星日^[17]。

通常,UT1在UTC(或TAI)时刻的值以改正的形式(UT1-UTC)给出,下文中以δ_UT1表示(UT1-UTC),设地球自转角速率变化量ΔΩ(t)=Ω(t)-Ω₀,则(7)式可改写为

δUT1(t)=δUT1(t0)+1Ω0∫t0tΔΩ(t)dt,(8)

式中:δ_UT1(t)表示在UTC的t时刻对应的UT1的改正值。光纤陀螺仪的优点是可以实时监测并提供地球自转角速度的测量值,利用其输出的角速度测量值可以计算出地球自转角速率的变化量ΔΩ(t),对t₀~t时间段内的ΔΩ(t)积分就可以得到δ_UT1的变化量ξ为

ξ=δUT1(t)-δUT1(t0)=1Ω0∫t0tΔΩ(t)dt。(9)

利用光纤陀螺仪测量数据理论上可以实现每秒计算并提供UT1参数的计算结果,该方法可以提高UT1获取的时间分辨率,同时该方法也可以依据用户的不同需求提供任意时刻的UT1参数结果,每次的计算结果以对应时间段内δ_UT1的变化量ξ的形式给出。

3 测量数据来源与处理流程

依托国家授时中心的光纤陀螺仪实验平台,利用光纤陀螺仪输出的地球自转角速度测量值验证UT1的解算方法。实验所用光纤陀螺仪光纤环直径为40 cm,光纤总长为30 km,所用光源的平均波长为1550 nm;陀螺仪采用双偏振双端口设计^[18],该设计既可有效抑制光纤过长而导致的温度梯度效应,也可以最大程度地减小光纤环平面的面积;利用多匝光纤环绕的方式实现放大Sagnac效应的效果。光纤陀螺仪经过初期的调试,现已稳定运行。光纤陀螺仪水平放置在实验室的恒温、恒湿的密闭环境中,尽可能地避免周围环境对陀螺仪的干扰,测量平台上安放有两台高精度电子水平仪,可以实时监测平台的倾斜变化情况。图2所示为2018-11-08UTC00:00:00—2018-11-13UTC00:00:00时间段的陀螺仪实测数据,该数据采样频率为10 Hz,在采集时间段内连续不间断测量。光纤陀螺仪输出的测量数据以旋转角速率(单位:rad/s)的形式给出。结合实测数据,依据(5)式的光纤陀螺仪观测方程,设计了如图3所示的UT1数据解算流程,由于目前光纤陀螺仪还无法完全测量出地球瞬时自转轴的变化,所以在接下来的分析中暂时忽略该部分的影响。

图 2. 大型光纤陀螺仪测量的地球自转角速率的原始数据

Fig. 2. Raw data of Earth’s rotation angular rate measured from large-scale FOG

下载图片 查看所有图片

图 3. 基于大型光纤陀螺仪的UT1数据解算流程

Fig. 3. Flow chart of UT1 calculation process by using data measured from large-scale FOG

下载图片 查看所有图片

4 数据分析与讨论

对光纤陀螺仪原始测量数据进行处理,首先需要消除测量数据中存在的若干数据异常值以及噪声,异常值并不是地球自转变化的真实反映,并且该值会对后续的数据处理过程产生影响。采用PauTa准则剔除数据中的异常值^[19],通过计算可知该段数据的标准偏差σ=2.07093×10^-8 rad/s,将实际测量值与平均值之差大于3σ的值视为异常值予以剔除,并用该段数据的算术平均值替换异常值。消除数据异常值后需要进一步消除测量数据中的噪声,利用光纤陀螺仪测量的地球真实自转角速度信息通常淹没在背景噪声中,该背景噪声通常具有白噪声的统计特性。利用Allan方差分析法可以分析光纤陀螺仪测量数据中的噪声特性^[20-21],如图4所示。在Allan方差分析曲线中,光纤陀螺仪的噪声随着平均时间的增加而逐渐减小,而平均时间越短,表示计算得到地球自转角速率的时间分辨率越高。对于实验所用光纤陀螺仪,取平均时间为5 min可以使该测量数据中的噪声降到最低。因此,采用对5 min内的光纤陀螺仪测量数据取算术平均值的方式来消除原始数据中的背景噪声,消噪后的结果如图5所示。所选取的平均时间5 min也是利用光纤陀螺仪的输出数据解算所得的UT1参数的时间分辨率,即通过光纤陀螺仪测量可以每隔5 min进行一次UT1参数的解算。

图 4. 光纤陀螺仪测量数据的Allan方差分析曲线

Fig. 4. Curve of Allan variance of data measured from FOG

下载图片 查看所有图片

图 5. 平均消噪后的数据结果(平均时间为5 min)

Fig. 5. Data after average denoising (average time is 5 min)

下载图片 查看所有图片

校正定向误差的目的是消除光纤陀螺仪因定向变化而产生的误差数据。本研究首先通过水平仪测量数据对陀螺仪进行倾斜改正,再利用固体潮理论模型消除因固体潮引起的地方铅垂线变化,最终实现定向误差的校正。定向误差校正结果如图6所示,其中由于地壳形变产生的倾斜变化量可以通过(5)式转换为等价的角速率变化。

图 6. 光纤陀螺仪测量数据的定向改正

Fig. 6. Orientation correction of data measured from FOG

下载图片 查看所有图片

消除定向误差后,利用(9)式可以计算得到5 min内的ξ_{5 min},进而通过(8)式可以实现在测量时间段内每隔5 min计算一个δ_UT1参数,将计算结果与IERS提供的δ_UT1结果进行对比,如图7所示,两段数据的起始时刻均为11月8日0时(UTC)。

图 7. 利用大型光纤陀螺仪测量数据解算的δUT1结果与IERS提供的数据对比

Fig. 7. Comparison between δUT1 calculated by data measured from large-scale FOG and δUT1 from IERS

下载图片 查看所有图片

通过对比可以看出,利用目前的光纤陀螺仪输出数据计算的δ_UT1值围绕着真实的δ_UT1上下波动,在一天内波动变化范围主要集中在±20 ms内。鉴于目前实验所用陀螺仪在5 min平均时间内测量角速度的灵敏度可达到3.15×10^-10 rad/s(对应产生的一天内δ_UT1的计算误差介于±500 ms之间),实际计算的波动范围并没有超出理论的误差范围,说明目前的计算结果已经基本达到了实验光纤陀螺仪的测量极限。

5 结论

通过上述实测数据分析,初步实现了利用大型光纤陀螺仪测量的地球自转角速度数据解算UT1的目标。利用目前的大型光纤陀螺仪的测量数据可以使UT1的获取时间间隔缩短为5 min,说明基于大型光纤陀螺仪可以提高UT1参数获取的时间分辨率,理论上可以实现实时提供UT1参数。但是,目前受限于光纤陀螺仪自身测量极限,计算结果与IERS提供的UT1参数结果还具有一定的差距,但是随着光纤陀螺仪技术的不断发展,测量灵敏度的不断提高,当灵敏度提高到7×10^-12 rad/s时,一天内通过光纤陀螺仪测量解算的UT1参数的波动范围可限制在±1 ms以内。该方法有潜力成为一种新型的UT1测量方法,通过大型光纤陀螺仪测量UT1也可以对目前现有的UT1测量方法(天顶筒、VLBI等)起到良好的互补作用。