地球重力场是用于表征地球系统属性的基本物理场，确定地球重力场的精细结构及其随时间的变化规律，历来是物理大地测量学及相关地球学科的重点和热点问题．随着卫星观测技术的不断完善，特别是CHAMP，GRACE和GOCE卫星计划的成功实施，促使地球重力场模型精度大幅提升[1]．相较于重力卫星时代之前精度最高的EGM96模型，采用8 a的CHAMP卫星数据解算的地球重力场模型在60阶次之前精度提升约1个量级[2]，采用13 a的GRACE卫星观测数据解算的地球重力场模型在100阶次之前精度提升约2~3个量级[3]，采用4.5 a的GOCE卫星观测数据解算的地球重力场模型在100~200阶次精度提升约1个量级[4]．考虑到地球重力场在地球科学和国防安全中的重要地位，借鉴国际重力卫星的成功经验，我国正在推进新一代自主重力卫星．由于传统模型化的方法无法有效分离出非保守力对卫星轨道的影响，在重力卫星核心载荷研制过程中，须设计高精度加速度计以实时获取卫星在轨受到的非保守力．在国内，华中科技大学引力中心研制的静电悬浮加速度计已成功搭载天琴一号(TQ-1，Tianqin-1)卫星，其加速度计分辨率为5×10-12 m·s-2·Hz-1/2@0.1 Hz，在轨评估精度达到1×10-10 m·s-2·Hz-1/2@0.1 Hz和5×10-11 m·s-2·Hz-1/2@0.05 Hz[5]．TQ-1卫星是一颗国防科工局立项的技术试验卫星，其主要任务是对未来开展引力波空间探测所必须的六项共性关键技术进行在轨验证，具体包括：高精度惯性传感技术；无拖曳控制技术；微牛级连续可调微推进技术；高精度激光干涉测量技术；高精度质心控制技术；高精度温度控制技术．同时，TQ-1卫星搭载了高精度GNSS定位系统，可为卫星提供精密轨道．利用高精度GNSS定位系统和静电悬浮加速度计，理论上可通过高低卫星跟踪卫星技术直接获取地球重力场信息；因此，这里以最新获取的TQ-1观测数据为基础，评估我国高低卫星跟踪卫星技术获取地球重力场信息的能力．鉴于卫星重力技术的蓬勃发展，国内外诸多学者研制了多种解算方法，以获取高精度高分辨率的地球重力场模型．最常用的有动力法[6-9]、加速度法[10]、能量法[11]、短弧法[12-13]和天体力学法[14]．鉴于动力法可同时获取卫星精密轨道和地球重力场模型，国外官方机构通常采用该方法处理重力卫星数据，如美国德克萨斯大学空间研究中心(Center for Space Research，CSR)、美国宇航局喷气推进实验室(Jet Propulsion Laboratory，JPL)和德国地学研究中心(GeoForschungs Zentrum Potsdam，GFZ)．在国内，文献[15-19]也利用动力法获取了一系列地球重力场模型．考虑到动力法的诸多优点，本研究拟采用该方法实现TQ-1卫星精密定轨，并反演一组地球重力场模型．1 动力法基本原理及数据处理方法根据观测值和未知参数处理方式的不同，利用高低卫星跟踪卫星技术获取地球重力场的方法可分为一步法和两步法．一步法以原始跟踪数据为观测值，同时求解卫星精密轨道和地球重力场模型[14-15,20-21]；两步法首先获取卫星轨道，再利用卫星轨道解算地球重力场模型[18-19,22-26]．利用一步法解算地球重力场模型的代表机构有CSR，JPL和GFZ等，更多学者采用了两步法解算地球重力场模型，以ITSG-Grace2018，Tongji-Grace2018和HUST-Grace2020模型为代表的GRACE时变重力场模型序列，其精度优于CSR，JPL和GFZ等机构利用一步法最新解算的RL06模型．考虑到两步法相对高效并易于实现，采用两步法解算地球重力场模型．当利用高低卫星跟踪卫星数据反演地球重力场模型时，动力法以卫星轨道X(t)为观测值，通过建立卫星轨道摄动与初始状态向量误差ΔX(t0)、加速度计校准参数Δp和地球重力场位系数改正量Δx的线性关系，然后经过最小二乘直接法获取新的地球重力场模型．根据上述描述，动力法的观测方程可表示为ΔX(t)=X(t)-X¯(t)=∂X(t)∂X(t0)ΔX(t0)+∂X(t)∂pΔp+∂X(t)∂xΔx+εX, (1)式中：X(t)为原始轨道观测值，表示卫星在任意时刻的状态向量，包含了卫星的位置和速度信息；X¯(t)为积分轨道(即动力学轨道)，须要通过数值积分获取；ΔX(t)为轨道残差值，可用于评估动力学轨道精度．偏导数∂X(t)/∂X(t0)为状态转移矩阵，偏导数∂X(t)/∂p和∂X(t)/∂x为参数敏感矩阵，均可通过数值积分获取；εX为线性化误差．由动力学的基本原理可知，利用动力法反演地球重力场模型的步骤主要包括：以观测轨道X(t)为输入，计算出与观测轨道拟合度最高的动力学轨道X¯(t)，即动力学精密定轨；以轨道残差值ΔX(t)为输入，以状态转移矩阵和参数敏感矩阵为系数矩阵，利用最小二乘法反演地球重力场模型．1.1　动力学精密定轨以初始时刻的状态向量X(t0)(包括位置和速度)为初值，利用先验力模型计算出参考轨道X¯(t)．参考轨道的计算实际上是一个数值积分过程，计算中采用改进的Gauss-Jackson数值积分器[27]，积分器阶次选为14阶，积分步长选为5 s．先验力模型均采用最新发布的力模型(见表1)．10.13245/j.hust.220917.T001表1先验力模型摄动力及其说明摄动力力模型说明海洋潮汐EOT11a截断至120阶，含18个主潮波和234个次潮波N体摄动IERS2010太阳和月球的三体摄动及J2项的影响．采用DE421星历固体潮IERS2010含频率无关项和频率有关项固体极潮IERS2010计算了其对位系数C21和S21的影响海洋极潮Desai2002截断至60阶大气和海洋去混频AOD RL06截断至100阶，利用线性内插法计算其对各个历元观测值的影响相对论效应IERS2010仅考虑了Schwarzschild项的影响非保守力加速度计观测值每个弧段引入一组偏差因子除了上述保守力摄动模型，非保守力摄动的测量精度也与动力学定轨精度密切相关．为了获取高精度动力学轨道，选用高精度的加速度计观测值来量化卫星在轨所受的非保守力．由于受到加速度计内部噪声、卫星平台振动和观测环境温度等综合影响[28-29]，加速度计观测值不可避免地存在偏差；为了实现动力学精密定轨，加速度计的精密校准是一个关键环节．由于仅用30 h的数据参与计算，参考文献[30]的数据处理方法，在加速度计精密校准过程中未考虑尺度因子的影响，仅引入3个观测方向的偏差因子b．校准后的非保守力可表示为fcal=R[fobs+b]，(2)式中：fobs为加速度计观测值；R为星固系至惯性系的旋转矩阵，可由星象仪观测值获取．类似于GRACE卫星，TQ-1卫星的旋转矩阵R可由卫星姿态四元数(q1,q2,q3,q4)获取，即R=[R1,R2,R3]，(3)式中：R1=[q12-q22-q32+q42,2(q1q2+q3q4),2(q1q3-q2q4)];R2=[2(q1q2-q3q4),-q12+q22-q32+q42,2(q2q3+q1q4)];R3=[2(q1q3+q2q4),2(q2q3-q1q4),-q12-q22+q32+q42].基于上述摄动力模型和加速度计观测值，利用Gauss-Jackson数值积分器完成轨道积分获取参考轨道X¯(t)．利用卫星轨道观测值X(t)减去数值积分计算的参考轨道X¯(t)，得到式(1)中的轨道残差ΔX(t)；利用数值积分获取的状态转移矩阵∂X(t)/∂X(t0)和参数敏感矩阵∂X(t)/∂p构建设计矩阵，通过最小二乘平差法获取初始状态向量改正值ΔX(t0)和加速度计校准参数改正量Δp．利用上述求解的改正量修正初始状态向量和加速度计校准参数，迭代进行新一轮的轨道积分和未知数求解，直至轨道残差达到某一设定阈值．通过上述解算过程，可获取更新后的初始状态X(t0)和加速度计校准参数p，即实现了加速度计的精密校准和动力学精密定轨．1.2　地球重力场反演当利用式(1)反演地球重力场模型时，根据未知数处理方式的差异，可分为两种数据处理方式：保持动力学精密定轨获取的初始状态X(t0)和加速度计校准参数p不变，仅解算地球重力场位系数改正量Δx；以动力学精密定轨获取的初始状态X(t0)和加速度计校准参数p为参考值，同时解算三类未知数的改正量ΔX(t0)，Δp和Δx．文献[18，31]的研究结果表明：第一种处理方式依赖于先验重力场模型的精度，且存在隐性的正则化效应．鉴于此，选取第2种处理方式实现地球重力场反演．在实际数据处理中，考虑到数值积分误差和线性化误差的累积效应，选择固定弧长建立观测方程．将与弧段有关的未知参数(初始状态误差ΔX(t0)和加速度计校准参数Δp)记为局部参数xl，将地球重力场位系数改正量Δx记为全局参数xe，则每个弧段的观测方程可表示为l=[A,B][xl,xe]T+ε，(4)式中：l为观测向量；A和B为法方程系数值．根据式(4)可获取对应的法方程为NllNleNelNeexlxe=WlWe．(5)将式(5)的法方程进行约化，消除局部参数xl，则全局参数xe的解可表示为xe=(Nee-NelNll-1Nle)-1(We-WeNel-1Wl)．(6)利用解算得到的全局变量xe修正先验重力场模型，即可获取更新后的地球重力场模型．2 基于TQ-1数据的地球重力场反演2.1　反演地球重力场模型根据动力学基本原理，自主编写了反演地球重力场模型的动力学软件系统，包括动力学精密定轨和地球重力场反演两部分．其中，动力学精密定轨包括了数据批量预处理模块、Gauss-Jackson数值积分模块、摄动力精密建模模块和加速度计精密校准模块；地球重力场反演包括了观测方程构建模块、法方程构建及求解模块、反演模型精度评估模块．该软件已经成功应用于基于GRACE卫星数据的静态和时变重力场建模[3,18]．利用自主研制的数据处理软件，基于TQ-1观测数据实现了地球重力场反演．由于TQ-1卫星为引力波探测服务的实验卫星，搭载了多种载荷进行实验，主要任务是对未来开展引力波空间探测所必需的六项共性技术进行在轨验证，各项科学实验的观测任务和观测时间均有严格限制．TQ-1卫星搭载的加速度计也承担了多项科研任务，大部分工作时间用于完成高精度惯性传感在轨评估、无拖曳控制技术评估、微牛级连续可调微推进技术评估、高精度质心控制技术评估等工作．受到卫星能量限制，星载GNSS接收机、星载加速度计和星象仪仅在2020年8月7日11时34分42.983 s后同时开启，开展了约30 h的地球重力场应用试验；因此，仅选用该时段内的TQ-1卫星连续观测数据参与解算，具体选用时段为2020年8月7日12时0分—2020年8月8日18时0分，观测数据包括轨道数据、加速度计数据和姿态数据．由于搭载的加速度计可直接获取在轨非保守力，联合上述多种观测数据，即可实现基于高低卫星跟踪卫星技术的地球重力场反演．须要注意的是，由于观测仪器或观测环境的影响，观测数据不可避免地存在间断现象．对于不同类型的数据空白，处理方式各不相同：为保证轨道积分连续性，采用样条插值的方法内插获取间断区域的加速度计数据和星象仪数据；为减小插值误差的影响，轨道观测数据不进行内插，缺失部分数据不参与计算；为了保证数值积分的精度，除了轨道观测值保留为10 s采样率，其余观测值均重采样至5 s；为避免动力学轨道积分累积误差的影响，与GRACE卫星数据处理采用相同的处理策略，计算中积分弧长选为6 h．计算中首先以式(2)为基础，实现了加速度计精密校准，校准参数仅选取偏差因子．在完成动力学精密定轨后，依照1.2节的地球重力场反演方法，利用TQ-1观测数据反演地球重力场模型．由于仅采用30 h的观测数据参与反演，地球重力场截断阶次仅选为15阶．2.2　观测数据及校准结果分析首先利用地固系下的轨道观测数据计算了TQ-1卫星的星下点轨迹，结果如图1所示．将各观测历元的卫星状态向量转换为开普勒轨道根数，统计了开普勒轨道根数的平均值．TQ-1卫星的平均轨道高度为629.03 km，倾角为97.85°，离心率为0.001 3，其近圆近极地的轨道特性使得星下点轨迹几乎覆盖全球，两极存在大约7.85°的极空白区域；在飞行30 h的情况下，能够实现全球覆盖，但轨道分布较为稀疏，且主要呈现南北分布．10.13245/j.hust.220917.F001图1TQ-1卫星的飞行轨迹(地固系)其次，完成了加速度计的精密校准．以TQ-1卫星各个弧段第一个历元的位置和速度为初始状态，利用本文第1.1节提及的先验力模型计算各时刻的保守力，利用加速度计实际观测值获取非保守力，然后利用Gauss-Jackson数值积分器进行轨道积分，获取参考轨道X¯(t)．当计算非保守力时，加速度计校准因子仅包含偏差因子，初始值设为零；利用星象仪获取的姿态数据构建旋转矩阵，具体计算方法见式(3)，然后将星固系下的非保守力转换到惯性系下．完成轨道积分后，将TQ-1观测的观测轨道X(t)和计算的参考轨道X¯(t)相减，获取轨道残差ΔX(t)．仅以式(1)中的初始状态向量ΔX(t0)和加速度计校准参数Δp为未知数，以轨道残差为观测值ΔX(t)，通过式(4)构建观测方程，并利用最小二乘平差法获取加速度计校准因子．计算的加速度计偏差因子在三方向的值如表2所示．10.13245/j.hust.220917.T002表2TQ-1加速度计每个弧段各个方向的偏差因子弧段x/(10-5m∙s-2)y/(10-7m∙s-2)z/(10-6m∙s-2)1-1.944-4.4533.7772-1.957-2.4453.6813-1.987-2.7523.4554-1.964-5.6573.5065-1.973-2.7673.803为了便于分析加速度计校准因子的特性，图2给出了加速度计x方向的观测值．由图2可知：TQ-1加速度计观测值在x方向在10-5 m/s2量级，这主要与TQ-1卫星轨道特征有关，该卫星的平均高度为629.03 km，沿轨主要受到大气阻力和太阳光压的影响．相应地，TQ-1加速度计在该方向偏差因子的变化范围为-1.944×10-5~-1.973×10-5 m/s2，变化幅度仅为3×10-7 m/s2，表明卫星在x方向受到的各弧段非保守力变化差异较小．10.13245/j.hust.220917.F002图2加速度计x方向观测的非保守力对比表2和图2的结果也可知：在30 h的观测时段内，x方向的加速度计观测值逐步减小，对应的加速度计校准参数也逐步减小；进一步计算了5个弧段对应的加速度计观测平均值，分别为-1.965×10-5，-1.969×10-5，-1.973×10-5，-1.975×10-5和1.978×10-5 m/s2，与表2各弧段的加速度计校准偏差因子差异较小．值得注意的是：本研究在加速度计校准中的偏差因子初始值设置为零，尽管偏差因子初始值设置为零，最终解算的各个弧段之间的偏差因子之间差异较小，这一特征与GRACE实测数据处理结果得出的结论类似[3]．为了对比校准参数的先验值对反演结果的影响，进行了两组对比实验：采用弧段内加速度计观测值的平均值；先验值设为0．利用上述两种方式进行动力学精密定轨时，第一种方法仅需一次迭代即可收敛，而第二种方法需两次迭代，但两组对比实验获取的校准参数完全相同．基于TQ-1数据的校准结果均表明：最终解算的偏差因子和平均值之间差异较小．特别地，以零值为加速度计校准参数的先验值，通常需要两次迭代才能够使迭代收敛，而利用弧段内的加速度计平均值为校准参数的初始值，仅需一次迭代即可收敛．鉴于上述经验，在后续发展的重力卫星数据处理过程中，为提高解算效率，可将观测弧段的平均值作为加速度计偏差因子的初始值提供给科学用户．2.3　地球重力场反演结果分析受限于TQ-1卫星观测环境，仅能够采用30 h观测数据解算地球重力场模型；因此，首先利用30 h的模拟数据分析了其反演地球重力场模型的可行性．模拟计算中，为尽可能接近真实观测环境，以TQ-1和各弧段的第一历元观测值为初始状态向量，以GGM05s模型为“真实”静态重力场模型，利用表1中的力模型计算卫星受到的摄动力，利用Gauss-Jackson数值积分器计算“真实”观测轨道；为了接近TQ-1卫星轨道噪声水平，在TQ-1模拟轨道中加入与文献[32]评价结果一致的轨道噪声，即分别加入标准差为1.79，6.27和6.02 cm的观测噪声．值得注意的是，这里的轨道噪声不是白噪声，而是与飞行轨道特性相关的噪声(详见文献[32]图13)．利用模拟30 h的数据解算地球重力场模型，具体解算方法见1.2节．基于上述模拟策略可知：理论上反演模型越接近于“真实”模型GGM05s，则反演精度越高．将反演得到的地球重力场与GGM05s模型作差，计算反演模型的大地水准面阶次误差．由图3可知：利用30 h的TQ-1模拟轨道观测数据获取的地球重力场模型，在15阶之前的大地水准面阶次误差均小于地球重力场信号，即信噪比大于1．模拟计算结果表明：利用30 h的TQ-1数据可实现15阶地球重力场模型的反演．10.13245/j.hust.220917.F003图3模拟解算模型的大地水准面阶误差在完成上述模拟验证后，利用TQ-1实测数据反演了地球重力场模型．反演过程中，以2.1节中动力学精密定轨后的初始状态向量和加速度计校准参数为初值，将式(1)中的初始状态向量、加速度计校准参数和地球重力场位系数均设为未知数，将所有未知数同时求解，最终获取地球重力场模型TQ01s-20200807，具体解算方法见第1.2节．下面通过多种方式评估TQ-1数据反演地球重力场模型的精度．首先给出了TQ-1卫星重力场模型的内符合精度(即formal error)，结果如图4所示．采用最小二乘平差法解算的地球重力场模型，其内符合精度主要反映了卫星轨道特性对解算结果的影响．TQ-1卫星重力场模型的内符合误差在低阶项、带谐项和近带谐项部分略微偏大，而扇谐项部分偏小；内符合精度主要与轨道特性密切相关的：TQ-1卫星为近极地轨道，星下点轨迹全球覆盖，在赤道附近较为稀疏，从而导致仅利用30 h数据解算的重力场在2阶位系数反演结果较差；由于TQ-1卫星轨道倾角为97.85°，导致带谐项和近带谐项部分的内符合误差略微偏大．类似的结果也出现在GOCE卫星解中，其轨道倾角为96.7°，带谐项的位系数精度最差[33-34]．10.13245/j.hust.220917.F004图4TQ-1重力场模型的位系数内符合误差(色标单位：log10)其次，将最终解算的重力场模型与AIUB-CHAMP03S模型作差，获得位系数误差，以评估解算模型的外符合误差．AIUB-CHAMP03S是基于2002—2009年的CHAMP卫星轨道观测数据解算的重力场模型，可用于评估本研究中仅采用30 h观测数据反演的地球重力场模型．TQ-1卫星重力场模型的大地水准面阶次误差和位系数误差分别如图 5和图6所示．10.13245/j.hust.220917.F005图5TQ-1重力场模型的大地水准面阶误差10.13245/j.hust.220917.F006图6TQ-1重力场模型的位系数外符合误差(色标单位：log10)由图5可知：利用30 h的观测数据反演的TQ-1地球重力场模型在15阶以内的信噪比均大于1，可有效反演该频段内的地球重力场信息，且在15阶时的大地水准面阶次误差优于2 cm．由图6可知：由于赤道附近观测数据分布较为稀疏，2阶位系数精度最低，仅达到10-8量级；3阶以上的位系数精度较高，位系数误差均优于10-9量级，且有部分位系数的精度可达到10-10量级．最后，利用TQ-1卫星反演地球重力场模型计算了全球格网重力异常，结果分别如图7(a)所示；为了对比，同时给出了HUST-Grace2016s模型计算的全球格网重力异常，结果如图7(b)所示．类似地，图8给出了全球格网大地水准面的计算结果，并给出了与HUST-Grace2016s模型计算结果的差异．在计算过程中，HUST-Grace2016s模型均截断至160阶次．10.13245/j.hust.220917.F007图7TQ-1反演重力场(15阶)和HUST-Grace2016s模型(160阶)计算的全球格网重力异常(色标单位：mGal)10.13245/j.hust.220917.F008图8TQ-1反演重力场(15阶)和HUST-Grace2016s模型(160阶)计算的全球格网大地水准面及两者之差(色标单位：m)对比图7和图8结果可知：利用TQ-1卫星和HUST-Grace2016s模型计算的全球格网重力异常和全球格网大地水准面在空间分布特征上均较为一致；由于HUST-Grace2016s是采用13 a的GRACE卫星观测数据反演的重力场模型，且包含了高低卫星跟踪卫星和低低卫星跟踪卫星观测信息[3]，计算的全球格网重力异常在空间分辨率上显著优于TQ-1卫星结果．天琴二号卫星将在天琴一号卫星的基础上发展高精度激光测距技术，实现低低卫星跟踪卫星精确测量，星间测距设计指标优于GRACE卫星，预期将进一步提升我国重力卫星观测精度和重力场建模精度．3 结论基于TQ-1观测数据实现了地球重力场模型反演，该模型为首个采用我国自主研制的高低卫星跟踪卫星技术解算的地球重力场模型．利用30 h的观测数据实现了加速度计精密校准、动力学精密定轨和地球重力场模型反演，得到如下结论．a．利用模拟数据计算结果表明，利用30 h的高低卫星跟踪卫星技术可实现15阶次地球重力场的解算．在此基础上，利用实测TQ-1卫星数据反演的地球重力场在15阶次内的信噪比均大于1，表明利用实测数据可实现地球重力场反演，且大地水准面反演精度在15阶次处优于2 cm．b．由于仅采用30 h观测数据计算，且轨道倾角为97.85°，TQ-1卫星星下点轨迹极空白较大，利用TQ-1卫星反演的地球重力场在2阶项、带谐项和近带谐项位系数的精度较低，但整体位系数精度均优于10-8量级，3阶以上的位系数精度均优于10-9量级．c．利用TQ-1卫星反演的重力场计算了全球格网重力异常和全球格网大地水准面，该计算结果与HUST-Grace2016s模型在空间分布上较为一致．由于HUST-Grace2016s模型包含了13 a的GRACE观测数据，且为高低卫星跟踪卫星技术和低低卫星跟踪卫星技术联合解算模型，其空间分辨率显著优于TQ-1反演重力场模型．上述研究结果表明：我国已掌握利用高低卫星跟踪卫星技术获取地球重力场的能力，对我国今后发展自主重力卫星具有重要参考意义．由于TQ-1卫星并非专用重力卫星，轨道特性、测量时间和测量环境等因素制约了地球重力场模型反演精度．随着天琴二号卫星(TQ-2)的成功立项，利用TQ-2卫星的高精度星间激光测距技术，势必可获取更高精度和分辨率的地球重力场模型．