CN103513294B - 一种低低星星跟踪卫星重力场测量性能解析计算方法 - Google Patents

Abstract

一种低低星星跟踪卫星重力场测量性能解析计算方法，获取低低星星跟踪重力卫星系统参数；计算重力卫星载荷测量误差对地球引力非球形摄动位功率谱的影响，进而得到反演重力场模型的位系数阶误差方差；与Kaula准则给出的位系数阶方差比较；计算反演重力场模型的大地水准面阶误差及其累积误差、重力异常阶误差及其累积误差；将计算得到的重力场测量有效阶数、大地水准面阶误差及其累积误差、重力异常阶误差及其累积误差汇总，即为低低星星跟踪重力场测量性能。本发明可以快速、定量评估重力场测量效果，获取重力卫星系统参数对重力场测量性能的影响规律，避免了卫星重力场测量数值模拟所带来的计算时间长、无法获取系统参数影响规律等缺陷。

Description

一种低低星星跟踪卫星重力场测量性能解析计算方法

技术领域

本发明属于卫星重力场测量技术领域，涉及解析计算低低星星跟踪卫星重力场测量性能，用于低低星星跟踪重力场测量卫星系统的参数设计。

背景技术

重力场是地球的基本物理场，在地球科学研究、国土资源勘探和地质灾害预报等方面具有重要应用，历来是大地测量学研究的核心问题。随着航天技术的发展，卫星重力场测量以其全球高覆盖率、全天候、不受地缘政治和地理环境影响等独特优势，受到了越来越多的重视，在理论研究和工程实践上均取得了长足发展，已成为获取全球重力场模型的最有效手段[1,2]。

根据观测数据的不同，卫星重力场测量可以分为轨道摄动、低低星星跟踪和重力梯度三种原理[3]。其中，轨道摄动原理适宜于低阶重力场测量，它的主要观测数据是卫星摄动轨道；低低星星跟踪原理适宜于中高阶重力场测量，它的主要观测数据是两个低轨卫星之间的距离及其变化率；重力梯度原理适宜于高阶重力场测量，它的主要观测数据是重力梯度值。已成功实施或正在研制的重力卫星均采用了以上测量原理或其组合，如CHAMP卫星利用轨道摄动原理恢复低阶重力场，GRACE、GRACEFollow-on、NGGM卫星同时利用了轨道摄动和低低星星跟踪原理恢复中高阶重力场，GOCE卫星分别利用轨道摄动和重力梯度原理恢复低阶和高阶重力场。虽然这三种原理适宜于不同的重力场测量频段，但是重力场测量有效阶数和精度最终要取决于重力卫星的载荷指标。在针对中高阶重力场测量的低低星星跟踪和重力梯度方式下，以目前的载荷性能指标分析可知，低低星星跟踪测量完全可以达到甚至超过重力梯度的测量水平。为此，2007年在荷兰召开的“未来重力卫星测量”专题研讨会上决定，目前国际重力卫星继续采用低低星星跟踪重力场测量方式，同时考虑提高改善系统参数以提高重力场测量性能。

针对低低星星跟踪重力场测量系统轨道参数和载荷指标设计，传统上主要采用数值模拟法分析系统参数对重力场测量性能的影响，进而确定系统设计参数。但是，重力场测量数值模拟对计算机性能要求非常高，计算时间非常长，不利于分析系统参数对重力场测量的影响规律，不便于进行系统参数的优化设计。为克服这一缺陷，本发明从能量守恒原理出发，建立了分析低低星星跟踪重力场测量性能的解析方法，可以快速获取重力场测量有效阶数、大地水准面误差和重力异常误差等重力场测量性能，确定系统参数对重力场测量的影响规律，对低低星星跟踪重力卫星系统参数设计具有重要指导意义。

引用的文献

[1]宁津生.跟踪世界发展动态致力地球重力场研究.武汉大学学报信息科学版,2001,26(6):471-474.

[2]宁津生.卫星重力探测技术与地球重力场研究.大地测量与地球动力学,2002,22(1):1-5.

[3]谷振丰,刘红卫,王兆魁,张育林.基于引力位系数相对权重的卫星重力场测量分析.地球物理学进展,2013,28(1):17-23.

发明内容

本发明的目的在于，针对低低星星跟踪重力场测量卫星系统，建立了重力场测量有效阶数、大地水准面误差和重力异常误差等重力场测量性能与卫星轨道高度、星间距离、星间距离变化率测量精度、定轨精度、非引力干扰、测量数据采样间隔、任务周期等卫星系统参数之间的解析关系，可以快速计算低低星星跟踪重力场测量性能，分析系统参数对重力场测量性能的影响规律，指导低低星星跟踪重力场测量卫星系统参数优化设计。

为了实现上述目的，本发明采用的技术方案如下：

一种低低星星跟踪卫星重力场测量性能解析计算方法，包括步骤如下：

步骤1：获取低低星星跟踪重力卫星系统参数；

步骤2：根据低低星星跟踪重力卫星系统参数，计算重力卫星载荷测量误差对地球引力非球形摄动位功率谱的影响，进而得到反演重力场模型的位系数阶误差方差；

步骤3：将得到的反演重力场模型位系数阶误差方差与Kaula准则给出的位系数阶方差比较；随着重力场模型阶数的增加，位系数阶误差方差逐渐增加，而位系数阶方差则逐渐减小，当阶误差方差等于阶方差时，认为达到重力场测量的最高有效阶数；

步骤4：根据反演重力场模型的阶误差方差，计算反演重力场模型的大地水准面阶误差及其累积误差、重力异常阶误差及其累积误差；

步骤5：将计算得到的重力场测量有效阶数、大地水准面阶误差及其累积误差、重力异常阶误差及其累积误差汇总，即为低低星星跟踪重力场测量性能。

优选的，步骤1中的所述低低星星跟踪重力卫星系统参数，包括但不限于重力卫星系统轨道参数和所述重力卫星系统载荷指标；

所述卫星重力场测量性能指标包括重力场反演的最高有效阶数N_max、n阶对应的大地水准面阶误差△_n、n阶对应的大地水准面累积误差△、n阶对应的重力异常误差△g_n和n阶对应的重力异常累积误差△g中的一种或几种；

所述重力卫星系统轨道参数，包括重力卫星轨道高度h和两星的地心矢量夹角θ₀；

所述重力卫星系统载荷指标，包括星间距离变化率测量误差、卫星定轨位置误差(△r)_m、非引力干扰△F、星间距离变化率数据采样间隔卫星轨道位置数据采样间隔(△t)_△r、非引力干扰数据间隔(△t)_△F和重力场测量任务寿命T。

优选的，步骤2具体包括：

建立关于低低星星跟踪重力场测量位系数阶误差方差δσ_n ²满足的解析关系式：

$\begin{array}{c}{{\mathrm{\δ\σ}}_n}^2=\frac{1}{2n+1}\underset{k=0}{\overset{n}{\Σ}}\[{\left(\mathrm{\δ}{\stackrel{\‾}{C}}_{nk}\right)}^2+{\left(\mathrm{\δ}{\stackrel{\‾}{S}}_{nk}\right)}^2\]=\frac{1}{\underset{k=0}{\overset{n}{\Σ}}\[B_1(r_0,n,k,{\mathrm{\θ}}_0)+B_2(r_0,n,k,{\mathrm{\θ}}_0)+B_3(r_0,n,k,{\mathrm{\θ}}_0)\]}\×\\ \frac{1}{2n+1}\frac{2(n+1){r_0}^{2n}}{{\mathrm{\π\μ}}^2{{\mathrm{Ta}}_e}^{2n-2}}{\[-D\sqrt{(\frac{{\left(\mathrm{\Δ}F\right)}^2{T_{arc}}^2}{36}{\left(\mathrm{\Δ}t\right)}_{\mathrm{\Δ}F}+{{\left(\mathrm{\Δ}r\right)}_m}^2{\left(\mathrm{\Δ}t\right)}_{\mathrm{\Δ}r})f_{\mathrm{\δ}r}\left(n\right)}+\sqrt{\frac{\mathrm{\μ}}{r_0}}\sqrt{(\frac{{\left(\mathrm{\Δ}F\right)}^2{T_{arc}}^2}{4}{\left(\mathrm{\Δ}t\right)}_{\mathrm{\Δ}F}+{{\left(\mathrm{\Δ}\stackrel{\·}{\mathrm{\ρ}}\right)}_m}^2{\left(\mathrm{\Δ}t\right)}_{\mathrm{\Δ}\stackrel{\·}{\mathrm{\ρ}}})f_{\mathrm{\δ}\stackrel{\·}{\mathrm{\ρ}}}\left(n\right)}\]}^2\end{array}$

其中，

$\left\{\begin{array}{c}{B}_1(r_0,n,k,{\mathrm{\θ}}_0)={A_n}^2(n-1,k,{\mathrm{\θ}}_0)+\left(\frac{a_e}{r_0}\right)\left|A_n(n-1,k,{\mathrm{\θ}}_0)A_n(n,k,{\mathrm{\θ}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^2\left|A_n(n-1,k,{\mathrm{\θ}}_0)A_n(n+1,k,{\mathrm{\θ}}_0)\right|\\ B_2(r_0,n,k,{\mathrm{\θ}}_0)={\left(\frac{a_e}{r_0}\right)}^2{A_n}^2(n,k,{\mathrm{\θ}}_0)+\left(\frac{a_e}{r_0}\right)\left|A_n(n-1,k,{\mathrm{\θ}}_0)A_n(n,k,{\mathrm{\θ}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^3\left|A_n(n,k,{\mathrm{\θ}}_0)A_n(n+1,k,{\mathrm{\θ}}_0)\right|\\ B_3(r_0,n,k,{\mathrm{\θ}}_0)={\left(\frac{a_e}{r_0}\right)}^4{A_n}^2(n+1,k,{\mathrm{\θ}}_0)+{\left(\frac{a_e}{r_0}\right)}^2\left|A_n(n-1,k,{\mathrm{\θ}}_0)A_n(n+1,k,{\mathrm{\θ}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^3\left|A_n(n,k,{\mathrm{\θ}}_0)A_n(n+1,k,{\mathrm{\θ}}_0)\right|\end{array}\right.$

$A_n(l,k,{\mathrm{\θ}}_0)={\mathrm{\δ}}_k\underset{0}{\overset{\mathrm{\π}}{\∫}}\[{\stackrel{\‾}{P}}_{lk}\left(\cos \right(\mathrm{\θ}+{\mathrm{\θ}}_0\left)\right)-{\stackrel{\‾}{P}}_{lk}\left(\cos \mathrm{\θ}\right)\]{\stackrel{\‾}{P}}_{nk}\left(\cos \mathrm{\θ}\right)\sin \mathrm{\θ}d\mathrm{\θ}$

${\mathrm{\δ}}_k=\left\{\begin{array}{cc}2,& k=0\\ 1,& k\≠0\end{array}\right.$

$D=K(\frac{cos\frac{{\mathrm{\θ}}_0}{2}}{2}-\frac{4}{3})\frac{{\mathrm{\μsin\θ}}_0}{{r_0}^4}$

$f_{\mathrm{\δ}r}\left(n\right)=\underset{k=0}{\overset{n}{\Σ}}\underset{0}{\overset{\mathrm{\π}}{\∫}}{\[{\stackrel{\‾}{P}}_{nk}\left(cos\mathrm{\θ}\right)\]}^2{\cos }^2(\mathrm{\θ}+\mathrm{\ξ}){\sin }^2\mathrm{\θ}d\mathrm{\θ}=\frac{3}{8}(2n+1)\mathrm{\π}$

$f_{\mathrm{\δ}\stackrel{\·}{\mathrm{\ρ}}}\left(n\right)=(n+\frac{1}{2})\mathrm{\π}$

r₀＝a_e+h

δσ_n ²是反演重力场模型位系数的阶误差方差，r₀是卫星的地心距，θ₀是两星地心矢量的夹角，a_e是地球平均半径，μ是万有引力常数和地球质量的乘积，h是卫星的轨道高度，T_arc是积分弧长，△F是非引力干扰，(△t)_△F是非引力干扰数据间隔，(△r)_m是卫星定轨的位置误差，(△t)_△r是卫星轨道数据的采样间隔，是星间距离变化率测量误差，是星间距离变化率数据采样间隔，T是重力场测量的任务寿命。

优选的，系数K和相位ξ的取值分别为：

$K=1.3476\×{10}^{11}{m}^2,\mathrm{\ξ}=-\frac{\mathrm{\π}}{2}$

优选的，对于非引力干扰数据间隔(△t)_△F，如果低低星星跟踪重力卫星系统对非引力干扰进行测量，那么(△t)_△F指非引力干扰的测量间隔；如果低低星星跟踪重力场测量系统对非引力干扰进行抑制，那么(△t)_△F指非引力干扰的抑制间隔。

优选的，步骤3具体包括建立的关于重力场测量最高有效阶数N_max满足解析关系式：

${{\mathrm{\δ\σ}}_{N_{max}}}^2=\frac{1}{2N_{max}+1}\frac{1.6\×{10}^{-10}}{{N_{max}}^3}$

优选的，步骤4具体包括：

建立的关于n阶对应的大地水准面阶误差为：

${\mathrm{\Δ}}_n=a_e\sqrt{(2n+1){{\mathrm{\δ\σ}}_n}^2}$

和/或

建立的关于n阶对应的大地水准面累积误差为：

$\mathrm{\Δ}=\sqrt{\underset{k=2}{\overset{n}{\Σ}}{\left({\mathrm{\Δ}}_k\right)}^2}$

和/或

建立的关于n阶对应的重力异常阶误差为：

${\mathrm{\Δg}}_n=\frac{\mathrm{\μ}}{{a_e}^2}(n-1)\sqrt{(2n+1){{\mathrm{\δ\σ}}_n}^2}$

和/或

建立的关于n阶对应的重力异常累积误差为：

$\mathrm{\Δ}g=\sqrt{\underset{k=2}{\overset{n}{\Σ}}{\left({\mathrm{\Δg}}_k\right)}^2}$

本发明“一种低低星星跟踪卫星重力场测量性能解析计算方法”包括四部分，分别是获取低低星星跟踪重力卫星系统参数、计算重力场测量的有效阶数、确定大地水准面阶误差及其累积误差、确定重力异常阶误差及其累积误差。其中，获取低低星星跟踪重力卫星系统参数，就是得到低低星星跟踪重力卫星系统的轨道高度、星间距离、星间距离变化率测量精度、定轨精度、非引力干扰、测量数据采样间隔、任务寿命等参数，为重力场测量性能计算提供输入参数；计算重力场测量的有效阶数，就是根据所得到的低低星星跟踪重力卫星系统参数，利用重力场测量性能解析分析公式计算反演重力场模型位系数的阶误差方差，然后将阶误差方差与Kaula准则给出的位系数阶方差比较，当阶误差方差等于阶方差时认为达到重力场测量的最高有效阶数；确定大地水准面阶误差及其累积误差，就是根据上一步得到的反演重力场模型的位系数阶误差方差，利用大地水准面计算公式得到有效阶数范围内的大地水准面阶误差及其累积误差；确定重力异常阶误差及其累积误差，就是根据反演重力场模型的位系数阶误差方差，利用重力异常计算公式得到有效阶数范围内的重力异常阶误差及其累积误差。

本发明“一种低低星星跟踪卫星重力场测量性能解析计算方法”实施步骤如下：

步骤1：获取低低星星跟踪重力卫星系统参数，包括轨道高度、星间距离、星间距离变化率测量精度、定轨精度、非引力干扰、测量数据采样间隔、任务寿命等；

步骤2：根据低低星星跟踪重力卫星系统参数，计算重力卫星载荷测量误差对地球引力非球形摄动位功率谱的影响，进而得到反演重力场模型的位系数阶误差方差；

步骤3：将得到的反演重力场模型位系数阶误差方差与Kaula准则给出的位系数阶方差比较。随着重力场模型阶数的增加，位系数阶误差方差逐渐增加，而位系数阶方差则逐渐减小，当阶误差方差等于阶方差时，认为达到重力场测量的最高有效阶数；

步骤4：根据反演重力场模型的阶误差方差，计算反演重力场模型的大地水准面阶误差及其累积误差、重力异常阶误差及其累积误差；

步骤5：将计算得到的重力场测量有效阶数、大地水准面阶误差及其累积误差、重力异常阶误差及其累积误差汇总，即为低低星星跟踪重力场测量性能。

本发明的优点在于：

建立了低低星星跟踪重力卫星测量性能与系统参数之间的解析关系，可以快速、定量评估重力场测量效果，获取重力卫星系统参数对重力场测量性能的影响规律，避免了卫星重力场测量数值模拟所带来的计算时间长、无法获取系统参数影响规律等缺陷。

附图说明

图1低低星星跟踪重力场测量卫星编队；

图2星间距离变化率幅值与r₀和θ₀的关系；

图3GRACE系统参数下的重力场测量位系数阶误差方差曲线；

图4GRACE重力场测量的大地水准面阶误差及其累积误差；

图5GRACE重力场测量的重力异常阶误差及其累积误差。

附表说明

表1不同轨道高度和星间地心矢量夹角下的星间距离变化率幅值(m/s)；

表2低低星星跟踪重力场测量中的物理参数；

表3GRACE卫星系统参数。

具体实施方式

为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

如图1-图5所示的一种低低星星跟踪卫星重力场测量性能解析计算方法，包括步骤如下：

步骤1：获取低低星星跟踪重力卫星系统参数；

步骤2：根据低低星星跟踪重力卫星系统参数，计算重力卫星载荷测量误差对地球引力非球形摄动位功率谱的影响，进而得到反演重力场模型的位系数阶误差方差；

步骤3：将得到的反演重力场模型位系数阶误差方差与Kaula准则给出的位系数阶方差比较；随着重力场模型阶数的增加，位系数阶误差方差逐渐增加，而位系数阶方差则逐渐减小，当阶误差方差等于阶方差时，认为达到重力场测量的最高有效阶数；

步骤4：根据反演重力场模型的阶误差方差，计算反演重力场模型的大地水准面阶误差及其累积误差、重力异常阶误差及其累积误差；

步骤5：将计算得到的重力场测量有效阶数、大地水准面阶误差及其累积误差、重力异常阶误差及其累积误差汇总，即为低低星星跟踪重力场测量性能。

在更加优选的是实施例中：

步骤1中的所述低低星星跟踪重力卫星系统参数，包括但不限于重力卫星系统轨道参数和所述重力卫星系统载荷指标；

所述卫星重力场测量性能指标包括重力场反演的最高有效阶数N_max、n阶对应的大地水准面阶误差△_n、n阶对应的大地水准面累积误差△、n阶对应的重力异常误差△g_n和n阶对应的重力异常累积误差△g中的一种或几种；

所述重力卫星系统轨道参数，包括重力卫星轨道高度h和两星的地心矢量夹角θ₀；

所述重力卫星系统载荷指标，包括星间距离变化率测量误差卫星定轨位置误差(△r)_m、非引力干扰△F、星间距离变化率数据采样间隔卫星轨道位置数据采样间隔(△t)_△r、非引力干扰数据间隔(△t)_△F和重力场测量任务寿命T。

步骤2具体包括：

建立关于低低星星跟踪重力场测量位系数阶误差方差δσ_n ²满足的解析关系式：

$\begin{array}{c}{{\mathrm{\δ\σ}}_n}^2=\frac{1}{2n+1}\underset{k=0}{\overset{n}{\Σ}}\[{\left(\mathrm{\δ}{\stackrel{\‾}{C}}_{nk}\right)}^2+{\left(\mathrm{\δ}{\stackrel{\‾}{S}}_{nk}\right)}^2\]=\frac{1}{\underset{k=0}{\overset{n}{\Σ}}\[B_1(r_0,n,k,{\mathrm{\θ}}_0)+B_2(r_0,n,k,{\mathrm{\θ}}_0)+B_3(r_0,n,k,{\mathrm{\θ}}_0)\]}\×\\ \frac{1}{2n+1}\frac{2(n+1){r_0}^{2n}}{{\mathrm{\π\μ}}^2{{\mathrm{Ta}}_e}^{2n-2}}{\[-D\sqrt{(\frac{{\left(\mathrm{\Δ}F\right)}^2{T_{arc}}^2}{36}{\left(\mathrm{\Δ}t\right)}_{\mathrm{\Δ}F}+{{\left(\mathrm{\Δ}r\right)}_m}^2{\left(\mathrm{\Δ}t\right)}_{\mathrm{\Δ}r})f_{\mathrm{\δ}r}\left(n\right)}+\sqrt{\frac{\mathrm{\μ}}{r_0}}\sqrt{(\frac{{\left(\mathrm{\Δ}F\right)}^2{T_{arc}}^2}{4}{\left(\mathrm{\Δ}t\right)}_{\mathrm{\Δ}F}+{{\left(\mathrm{\Δ}\stackrel{\·}{\mathrm{\ρ}}\right)}_m}^2{\left(\mathrm{\Δ}t\right)}_{\mathrm{\Δ}\stackrel{\·}{\mathrm{\ρ}}})f_{\mathrm{\δ}\stackrel{\·}{\mathrm{\ρ}}}\left(n\right)}\]}^2\end{array}$

其中，

$\left\{\begin{array}{c}{B}_1(r_0,n,k,{\mathrm{\θ}}_0)={A_n}^2(n-1,k,{\mathrm{\θ}}_0)+\left(\frac{a_e}{r_0}\right)\left|A_n(n-1,k,{\mathrm{\θ}}_0)A_n(n,k,{\mathrm{\θ}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^2\left|A_n(n-1,k,{\mathrm{\θ}}_0)A_n(n+1,k,{\mathrm{\θ}}_0)\right|\\ B_2(r_0,n,k,{\mathrm{\θ}}_0)={\left(\frac{a_e}{r_0}\right)}^2{A_n}^2(n,k,{\mathrm{\θ}}_0)+\left(\frac{a_e}{r_0}\right)\left|A_n(n-1,k,{\mathrm{\θ}}_0)A_n(n,k,{\mathrm{\θ}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^3\left|A_n(n,k,{\mathrm{\θ}}_0)A_n(n+1,k,{\mathrm{\θ}}_0)\right|\\ B_3(r_0,n,k,{\mathrm{\θ}}_0)={\left(\frac{a_e}{r_0}\right)}^4{A_n}^2(n+1,k,{\mathrm{\θ}}_0)+{\left(\frac{a_e}{r_0}\right)}^2\left|A_n(n-1,k,{\mathrm{\θ}}_0)A_n(n+1,k,{\mathrm{\θ}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^3\left|A_n(n,k,{\mathrm{\θ}}_0)A_n(n+1,k,{\mathrm{\θ}}_0)\right|\end{array}\right.$

$A_n(l,k,{\mathrm{\θ}}_0)={\mathrm{\δ}}_k\underset{0}{\overset{\mathrm{\π}}{\∫}}\[{\stackrel{\‾}{P}}_{lk}\left(\cos \right(\mathrm{\θ}+{\mathrm{\θ}}_0\left)\right)-{\stackrel{\‾}{P}}_{lk}\left(\cos \mathrm{\θ}\right)\]{\stackrel{\‾}{P}}_{nk}\left(\cos \mathrm{\θ}\right)\sin \mathrm{\θ}d\mathrm{\θ}$

${\mathrm{\δ}}_k=\left\{\begin{array}{cc}2,& k=0\\ 1,& k\≠0\end{array}\right.$

$D=K(\frac{cos\frac{{\mathrm{\θ}}_0}{2}}{2}-\frac{4}{3})\frac{{\mathrm{\μsin\θ}}_0}{{r_0}^4}$

$f_{\mathrm{\δ}r}\left(n\right)=\underset{k=0}{\overset{n}{\Σ}}\underset{0}{\overset{\mathrm{\π}}{\∫}}{\[{\stackrel{\‾}{P}}_{nk}\left(cos\mathrm{\θ}\right)\]}^2{\cos }^2(\mathrm{\θ}+\mathrm{\ξ}){\sin }^2\mathrm{\θ}d\mathrm{\θ}=\frac{3}{8}(2n+1)\mathrm{\π}$

$f_{\mathrm{\δ}\stackrel{\·}{\mathrm{\ρ}}}\left(n\right)=(n+\frac{1}{2})\mathrm{\π}$

r₀＝a_e+h

δσ_n ²是反演重力场模型位系数的阶误差方差，r₀是卫星的地心距，θ₀是两星地心矢量的夹角，a_e是地球平均半径，μ是万有引力常数和地球质量的乘积，h是卫星的轨道高度，T_arc是积分弧长，△F是非引力干扰，(△t)_△F是非引力干扰数据间隔，(△r)_m是卫星定轨的位置误差，(△t)_△r是卫星轨道数据的采样间隔，是星间距离变化率测量误差，(△t)_△ρ&是星间距离变化率数据采样间隔，T是重力场测量的任务寿命。

系数K和相位ξ的取值分别为：

$K=1.3476\×{10}^{11}{m}^2,\mathrm{\ξ}=-\frac{\mathrm{\π}}{2}$

步骤3具体包括建立的关于重力场测量最高有效阶数N_max满足解析关系式：

${{\mathrm{\δ\σ}}_{N_{max}}}^2=\frac{1}{2N_{max}+1}\frac{1.6\×{10}^{-10}}{{N_{max}}^3}$

步骤4具体包括：

建立的关于n阶对应的大地水准面阶误差为：

${\mathrm{\Δ}}_n=a_e\sqrt{(2n+1){{\mathrm{\δ\σ}}_n}^2}$

和/或

建立的关于n阶对应的大地水准面累积误差为：

$\mathrm{\Δ}=\sqrt{\underset{k=2}{\overset{n}{\Σ}}{\left({\mathrm{\Δ}}_k\right)}^2}$

和/或

建立的关于n阶对应的重力异常阶误差为：

${\mathrm{\Δg}}_n=\frac{\mathrm{\μ}}{{a_e}^2}(n-1)\sqrt{(2n+1){{\mathrm{\δ\σ}}_n}^2}$

和/或

建立的关于n阶对应的重力异常累积误差为：

$\mathrm{\Δ}g=\sqrt{\underset{k=2}{\overset{n}{\Σ}}{\left({\mathrm{\Δg}}_k\right)}^2}$

在更加优选的实施例中：

对于非引力干扰数据间隔(△t)_△F，如果低低星星跟踪重力卫星系统对非引力干扰进行测量，那么(△t)_△F指非引力干扰的测量间隔；如果低低星星跟踪重力场测量系统对非引力干扰进行抑制，那么(△t)_△F指非引力干扰的抑制间隔。

下面举例详细说明：

低低星星跟踪重力卫星系统由两颗低轨卫星沿迹向形成跟飞编队，通过测量两星之间的距离及其变化率来恢复地球重力场。为使测量数据精度尽可能一致，同时便于卫星姿态和轨道控制，通常选取卫星轨道为圆轨道或近圆轨道。为满足全球覆盖测量要求，通常选取卫星轨道为极轨道或近极轨道。为此，在建立低低星星跟踪重力场测量性能解析分析方法的过程中，假设卫星轨道偏心率为0，轨道倾角为90°。

1低低星星跟踪重力场测量的能量守恒方程

设组成低低星星跟踪重力场测量卫星编队的两个卫星分别为A和B，如图1所示。其中，B为引导星，A为跟踪星。卫星A和卫星B之间存在如下能量守恒关系

$\begin{array}{c}{V}_A-V_B\\ =\frac{1}{2}({{\stackrel{\·}{r}}_A}^2-{{\stackrel{\·}{r}}_B}^2)-\mathrm{\ω}\[{(r_x{\stackrel{\·}{r}}_y-r_y{\stackrel{\·}{r}}_x)}_A-{(r_x{\stackrel{\·}{r}}_y-r_y{\stackrel{\·}{r}}_x)}_B\]\underset{t_0}{\overset{t}{\∫}}(a_A\·{\stackrel{\·}{r}}_A-a_B\·{\stackrel{\·}{r}}_B)dt-V_{t,AB}-E_{0,AB}\end{array}---\left(1\right)$

在(1)式中，左边项为卫星A和卫星B的引力位差。已知引力位的球谐展开式为

$V(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{\μ}}{r}\[1+\underset{n=2}{\overset{\mathrm{\∞}}{\Σ}}\underset{k=0}{\overset{n}{\Σ}}{\left(\frac{a_e}{r}\right)}^n({\stackrel{\‾}{C}}_{nk}\cos k\mathrm{\λ}+{\stackrel{\‾}{S}}_{nk}\sin k\mathrm{\λ}){\stackrel{\‾}{P}}_{nk}\left(\cos \mathrm{\θ}\right)\]---\left(2\right)$

其中，(r,θ,λ)是地球固连坐标系中的球坐标，μ是地球引力常数，a_e是地球半径，和分别是n阶k次引力位系数的余弦项和正弦项，是完全规格化的缔合勒让德多项式。引力位V(r,θ,λ)可以分为中心引力位和非球形摄动位R(r,θ,λ)两部分

$V(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{\μ}}{r}+R(r,\mathrm{\θ},\mathrm{\λ})---\left(3\right)$

$R(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{\μ}}{r}\underset{n=2}{\overset{\mathrm{\∞}}{\Σ}}\underset{k=0}{\overset{n}{\Σ}}{\left(\frac{a_e}{r}\right)}^n({\stackrel{\‾}{C}}_{nk}\cos k\mathrm{\λ}+{\stackrel{\‾}{S}}_{nk}\sin k\mathrm{\λ}){\stackrel{\‾}{P}}_{nk}\left(cos\mathrm{\θ}\right)---\left(4\right)$

从而，(1)式的左边项可以表示为

$\begin{array}{c}{V}_A(r,\mathrm{\θ},\mathrm{\λ})-V_B(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{\μ}}{r_A}+R_A(r,\mathrm{\θ},\mathrm{\λ})-\frac{\mathrm{\μ}}{r_B}-R_B(r,\mathrm{\θ},\mathrm{\λ})\\ =(\frac{\mathrm{\μ}}{r_A}-\frac{\mathrm{\μ}}{r_B})+\[R_A(r,\mathrm{\θ},\mathrm{\λ})-R_B(r,\mathrm{\θ},\mathrm{\λ})\]\end{array}---\left(5\right)$

在(1)式中，等式右边第一项是卫星A和B的动能差。设e是卫星A指向卫星B的单位向量

$e=\frac{r_{AB}}{\left|r_{AB}\right|}---\left(6\right)$

卫星A和B的动能差可以表示为

$\frac{1}{2}({{\stackrel{\·}{r}}_A}^2-{{\stackrel{\·}{r}}_B}^2)=-\frac{1}{2}({\stackrel{\·}{r}}_A+{\stackrel{\·}{r}}_B)\·{\stackrel{\·}{r}}_{AB}=-\frac{1}{2}({\stackrel{\·}{r}}_A+{\stackrel{\·}{r}}_B)\·\{({\stackrel{\·}{r}}_{AB}\·e)e+\[{\stackrel{\·}{r}}_{AB}-({\stackrel{\·}{r}}_{AB}\·e)e\]\}---\left(7\right)$

考虑到卫星A和B的距离很近，它们的平均速度近似沿两星连线方向，即

$\frac{1}{2}({\stackrel{\·}{r}}_A+{\stackrel{\·}{r}}_B)\·\[({\stackrel{\·}{r}}_{AB}\·e)e\]>>\frac{1}{2}({\stackrel{\·}{r}}_A+{\stackrel{\·}{r}}_B)\·\[{\stackrel{\·}{r}}_{AB}-({\stackrel{\·}{r}}_{AB}\·e)e\]---\left(8\right)$

所以，(8)式可以近似表示为

$\frac{1}{2}({{\stackrel{\·}{r}}_A}^2-{{\stackrel{\·}{r}}_B}^2)\≈-\frac{1}{2}({\stackrel{\·}{r}}_A+{\stackrel{\·}{r}}_B)\·\[({\stackrel{\·}{r}}_{AB}\·e)e\]=-\frac{1}{2}({\stackrel{\·}{r}}_A+{\stackrel{\·}{r}}_B)\·\left(\stackrel{\·}{\mathrm{\ρ}}e\right)---\left(9\right)$

其中，是卫星A和B之间的距离变化率，等于在单位向量e上的投影。由于假设卫星A和B的轨道为圆轨道，同时考虑到A和B的地心距基本相等，设其平均值为r₀，则(9)式可以进一步表示为

$\frac{1}{2}({{\stackrel{\·}{r}}_A}^2-{{\stackrel{\·}{r}}_B}^2)\≈-\sqrt{\frac{\mathrm{\μ}}{r_0}}\stackrel{\·}{\mathrm{\ρ}}---\left(10\right)$

在(1)式中，等式右端第二项是由于地球自转引起的卫星A和B能量差，数值计算表明该项比动能差小4个数量级，在解析分析中可以忽略不计。(1)式等式右端第三项是非引力干扰引起的卫星A和B能量差。在重力场测量中，非引力干扰通过时间累积引起卫星轨道偏移纯引力轨道，同时引起星间距离变化率偏移纯引力作用下的星间距离变化率。其中，纯引力轨道偏移量与定轨误差共同决定了对纯引力轨道的获取能力，星间距离变化率偏移量和测量误差共同决定了对纯引力星间距离变化率的获取能力。这样非引力干扰的影响可以在分析纯引力轨道误差和纯引力星间距离变化率误差时引入。(1)式等式右端第四项是三体引力、潮汐摄动等引起的卫星A和B能量差。三体引力、潮汐摄动等具有高精度模型，可以满足静态重力场测量要求，对重力场测量的影响可忽略。(1)式右端第五项是积分常数，对重力场测量无影响。

将(5)和(10)式代入(1)式，得到

$R_A(r,\mathrm{\θ},\mathrm{\λ})-R_B(r,\mathrm{\θ},\mathrm{\λ})=-\sqrt{\frac{\mathrm{\μ}}{r_0}}\stackrel{\·}{\mathrm{\ρ}}-(\frac{\mathrm{\μ}}{r_A}-\frac{\mathrm{\μ}}{r_B})-E_{0,AB}---\left(11\right)$

对上式变分，得到

$\begin{array}{c}{\mathrm{\δR}}_A-{\mathrm{\δR}}_B-\frac{1}{2}\sqrt{\frac{\mathrm{\μ}}{{r_0}^3}}\stackrel{\·}{\mathrm{\ρ}}\left(\mathrm{\δ}r\right)-\sqrt{\frac{\mathrm{\μ}}{r_0}}\mathrm{\δ}\stackrel{\·}{\mathrm{\ρ}}+(\frac{\mathrm{\μ}}{{r_A}^2}-\frac{\mathrm{\μ}}{{r_B}^2})\left(\mathrm{\δ}r\right)\\ \≈\frac{1}{2}\sqrt{\frac{\mathrm{\μ}}{{r_0}^3}}\stackrel{\·}{\mathrm{\ρ}}\left(\mathrm{\δ}r\right)-\sqrt{\frac{\mathrm{\μ}}{r_0}}\mathrm{\δ}\stackrel{\·}{\mathrm{\ρ}}+\frac{2\mathrm{\μ}}{{r_0}^3}(r_B-r_A)\left(\mathrm{\δ}r\right)\end{array}---\left(12\right)$

其中，δr为纯引力轨道的位置误差，包括定轨误差和非引力干扰引起的纯引力轨道偏移量；为纯引力轨道星间距离变化率误差，包括星间测距误差和非引力干扰引起的星间距离变化率偏移量。下面建立纯引力轨道位置误差δr和纯引力轨道星间距离变化率误差之间的关系。

2两星地心距之差与星间距离变化率之间的几何关系

已知卫星A和B的星间距离变化率为

$\stackrel{\&CenterDot;}{\mathrm{\&rho;}}={\stackrel{\&CenterDot;}{r}}_{AB}\&CenterDot;e=\left|{\stackrel{\&CenterDot;}{r}}_{AB}\right|cos<{\stackrel{\&CenterDot;}{r}}_{AB},e>---\left(13\right)$

采用极坐标推导两星地心距差r_B-r_A与星间距离变化率之间的几何关系。在由r_A和r_B决定的平面内，如图1所示，设直线L₂平分r_A和r_B的夹角，直线L₁垂直于L₂。选取地心为极点，直线L₁为极轴，逆时针方向为角度正方向，则卫星A和B的位置矢量可以表示为

$r_A=r_A{e}^{i\frac{\mathrm{\&pi;}-{\mathrm{\&theta;}}_0}{2}}---\left(14\right)$

$r_B=r_B{e}^{i\frac{\mathrm{\&pi;}+{\mathrm{\&theta;}}_0}{2}}---\left(15\right)$

其中，θ₀是卫星A和B的地心矢量夹角。卫星A和B的速度矢量可以表示为

${\stackrel{\&CenterDot;}{r}}_A={\stackrel{\&CenterDot;}{r}}_A{e}^{i\frac{2\mathrm{\&pi;}-{\mathrm{\&theta;}}_0}{2}}---\left(16\right)$

${\stackrel{\&CenterDot;}{r}}_B={\stackrel{\&CenterDot;}{r}}_B{e}^{i\frac{2\mathrm{\&pi;}+{\mathrm{\&theta;}}_0}{2}}---\left(17\right)$

于是，可以得到卫星A和B的位置矢量差和速度矢量差为

$r_{AB}=r_B-r_A=r_B{e}^{i\frac{\mathrm{\&pi;}+{\mathrm{\&theta;}}_0}{2}}-r_A{e}^{i\frac{\mathrm{\&pi;}-{\mathrm{\&theta;}}_0}{2}}---\left(18\right)$

${\stackrel{\&CenterDot;}{r}}_{AB}={\stackrel{\&CenterDot;}{r}}_B-{\stackrel{\&CenterDot;}{r}}_A={\stackrel{\&CenterDot;}{r}}_B{e}^{i\frac{2\mathrm{\&pi;}+{\mathrm{\&theta;}}_0}{2}}-{\stackrel{\&CenterDot;}{r}}_A{e}^{i\frac{2\mathrm{\&pi;}-{\mathrm{\&theta;}}_0}{2}}---\left(19\right)$

在地球非球形引力摄动下，卫星A和B的地心距会产生周期性变化。设在时刻t卫星A的地心距为r₀，卫星B的地心距为r₀+△r₀。下面计算矢量e的幅角，将(18)式代入(6)式，得到

$e=\frac{r_{AB}}{\left|r_{AB}\right|}=\frac{(r_0+{\mathrm{\&Delta;r}}_0)e^{i\frac{\mathrm{\&pi;}+{\mathrm{\&theta;}}_0}{2}}-r_0{e}^{i\frac{\mathrm{\&pi;}-{\mathrm{\&theta;}}_0}{2}}}{\mathrm{\&rho;}}---\left(20\right)$

其中，ρ为卫星A和B的星间距离。整理(20)式，得到

$e=\frac{-(2r_0+{\mathrm{\&Delta;r}}_0)sin\frac{{\mathrm{\&theta;}}_0}{2}+{\mathrm{i\&Delta;r}}_{0}cos\frac{{\mathrm{\&theta;}}_0}{2}}{\mathrm{\&rho;}}---\left(21\right)$

从而得到e的幅角为

$\arg \left(e\right)=\mathrm{\&pi;}+arctan(-\frac{{\mathrm{\&Delta;r}}_0}{(2r_0+{\mathrm{\&Delta;r}}_0)tan\frac{{\mathrm{\&theta;}}_0}{2}})---\left(22\right)$

下面计算的幅角。在圆轨道假设条件下，卫星A和B的速度大小可近似为

${\stackrel{\&CenterDot;}{r}}_A=\sqrt{\frac{\mathrm{\&mu;}}{r_0}}---\left(23\right)$

${\stackrel{\&CenterDot;}{r}}_B=\sqrt{\frac{\mathrm{\&mu;}}{r_0+{\mathrm{\&Delta;r}}_0}}---\left(24\right)$

将(23)和(24)式代入(19)式，得到

${\stackrel{\&CenterDot;}{r}}_{AB}={\stackrel{\&CenterDot;}{r}}_B-{\stackrel{\&CenterDot;}{r}}_A=\sqrt{\frac{\mathrm{\&mu;}}{r_0+{\mathrm{\&Delta;r}}_0}}{e}^{i\frac{2\mathrm{\&pi;}+{\mathrm{\&theta;}}_0}{2}}-\sqrt{\frac{\mathrm{\&mu;}}{r_0}}{e}^{i\frac{2\mathrm{\&pi;}-{\mathrm{\&theta;}}_0}{2}}---\left(25\right)$

整理(25)式，得到

${\stackrel{\&CenterDot;}{r}}_{AB}=(\sqrt{\frac{\mathrm{\&mu;}}{r_0}}-\sqrt{\frac{\mathrm{\&mu;}}{r_0+{\mathrm{\&Delta;r}}_0}})cos\frac{{\mathrm{\&theta;}}_0}{2}-i(\sqrt{\frac{\mathrm{\&mu;}}{r_0+{\mathrm{\&Delta;r}}_0}}+\sqrt{\frac{\mathrm{\&mu;}}{r_0}})sin\frac{{\mathrm{\&theta;}}_0}{2}---\left(26\right)$

从而，得到虚部和实部的比为

$\frac{\mathrm{Im}\left({\stackrel{\&CenterDot;}{r}}_{AB}\right)}{\mathrm{Re}\left({\stackrel{\&CenterDot;}{r}}_{AB}\right)}=-\frac{(\sqrt{\frac{\mathrm{\&mu;}}{r_0+{\mathrm{\&Delta;r}}_0}}+\sqrt{\frac{\mathrm{\&mu;}}{r_0}})sin\frac{{\mathrm{\&theta;}}_0}{2}}{(\sqrt{\frac{\mathrm{\&mu;}}{r_0}}-\sqrt{\frac{\mathrm{\&mu;}}{r_0+{\mathrm{\&Delta;r}}_0}})\cos \frac{{\mathrm{\&theta;}}_0}{2}}\&ap;-\frac{(4r_0+{\mathrm{\&Delta;r}}_0)tan\frac{{\mathrm{\&theta;}}_0}{2}}{{\mathrm{\&Delta;r}}_0}---\left(27\right)$

从而，得到的幅角为

$\arg \left({\stackrel{\&CenterDot;}{r}}_{AB}\right)=arctan\left(\frac{{\mathrm{\&Delta;r}}_0}{(4r_0+{\mathrm{\&Delta;r}}_0)tan\frac{{\mathrm{\&theta;}}_0}{2}}\right)-\frac{\mathrm{\&pi;}}{2}---\left(28\right)$

由(22)和(28)式，得到e和之间的夹角为

$\begin{array}{c}<{\stackrel{\&CenterDot;}{r}}_{AB},e>=\arg \left({\stackrel{\&CenterDot;}{r}}_{AB}\right)-\arg \left(e\right)\\ =\arctan \left(\frac{{\mathrm{\&Delta;r}}_0}{(4r_0+{\mathrm{\&Delta;r}}_0)\tan \frac{{\mathrm{\&theta;}}_0}{2}}\right)+\arctan \left(\frac{{\mathrm{\&Delta;r}}_0}{(2r_0+{\mathrm{\&Delta;r}}_0)\tan \frac{{\mathrm{\&theta;}}_0}{2}}\right)-\frac{3\mathrm{\&pi;}}{2}\end{array}---\left(29\right)$

于是，

$cos<{\stackrel{\&CenterDot;}{r}}_{AB},e>=cos\&lsqb;arctan\left(\frac{{\mathrm{\&Delta;r}}_0}{(4r_0+{\mathrm{\&Delta;r}}_0)tan\frac{{\mathrm{\&theta;}}_0}{2}}\right)+arctan\left(\frac{{\mathrm{\&Delta;r}}_0}{(2r_0+{\mathrm{\&Delta;r}}_0)tan\frac{{\mathrm{\&theta;}}_0}{2}}\right)-\frac{3\mathrm{\&pi;}}{2}\&rsqb;---\left(30\right)$

化简，得到

$cos<{\stackrel{\&CenterDot;}{r}}_{AB},e>\&ap;-\frac{3{\mathrm{\&Delta;r}}_0}{4r_0\tan \frac{{\mathrm{\&theta;}}_0}{2}}---\left(31\right)$

由(13)和(31)式，得到

$\left|{\stackrel{\&CenterDot;}{r}}_{AB}\right|=\frac{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}{cos<{\stackrel{\&CenterDot;}{r}}_{AB},e>}=-\frac{4r_0\left(tan\frac{{\mathrm{\&theta;}}_0}{2}\right)\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}{3{\mathrm{\&Delta;r}}_0}---\left(32\right)$

由(26)式，得到

$|{\stackrel{\&CenterDot;}{r}}_{AB}{|}^2={(\sqrt{\frac{\mathrm{\&mu;}}{r_0}}-\sqrt{\frac{\mathrm{\&mu;}}{r_0+{\mathrm{\&Delta;r}}_0}})}^2{\cos }^2\frac{{\mathrm{\&theta;}}_0}{2}+{(\sqrt{\frac{\mathrm{\&mu;}}{r_0+{\mathrm{\&Delta;r}}_0}}+\sqrt{\frac{\mathrm{\&mu;}}{r_0}})}^2{\sin }^2\frac{{\mathrm{\&theta;}}_0}{2}---\left(33\right)$

由(32)和(33)式，得到两星地心距之差△r₀和星间距离变化率之间的关系为

$r_B-r_A={\mathrm{\&Delta;r}}_0=-\frac{2r_0}{3\cos \frac{{\mathrm{\&theta;}}_0}{2}}\sqrt{\frac{r_0}{\mathrm{\&mu;}}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}---\left(34\right)$

将(34)式代入(12)式，得到低低星星跟踪重力场测量的能量守恒关系满足

${\mathrm{\&delta;R}}_A-{\mathrm{\&delta;R}}_B=(\frac{1}{2}-\frac{4}{3cos\frac{{\mathrm{\&theta;}}_0}{2}})\sqrt{\frac{\mathrm{\&mu;}}{{r_0}^3}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\mathrm{\&delta;}r-\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}---\left(35\right)$

为了计算(35)式右端项的功率谱密度，需要得到随时间的变化规律。

3星间距离变化率的数学表示

由于地球引力的摄动影响，卫星轨道半长轴会出现周期性变化，导致星间距离变化率也出现周期性变化。其中，J₂项摄动是主要摄动项，其影响可以分为长期项和短周期项。J₂摄动的长期项对轨道半长轴无影响，短周期项对半长轴的影响为

$a_s\left(t\right)=\frac{3}{2}\frac{J_2}{a}\{\frac{2}{3}(1-\frac{3}{2}{\sin }^{2}i)\&lsqb;{\left(\frac{a}{r}\right)}^3-{(1-e^2)}^{-3/2}\&rsqb;+{\sin }^{2}i{\left(\frac{a}{r}\right)}^3\cos 2(f+\mathrm{\&omega;})\}---\left(36\right)$

其中，a是半长轴，i是轨道倾角，r是卫星地心距，e是轨道偏心率，f是真近点角，ω是近地点角距。对于低低星星跟踪重力卫星而言，有

从而，(36)式可以近似表示为

$r_s\left(t\right)=\frac{3}{2}\frac{J_2}{r}cos(2\mathrm{\&theta;}+\mathrm{\&xi;})---\left(38\right)$

其中，θ是余纬，ξ是初始相位角。从而，低低星星跟踪的两个卫星之间的瞬时地心距之差为

${\mathrm{\&Delta;r}}_0=\frac{3}{2}\frac{J_2}{r}\&lsqb;cos(2\mathrm{\&theta;}+2{\mathrm{\&theta;}}_0+\mathrm{\&xi;})-cos(2\mathrm{\&theta;}+\mathrm{\&xi;})\&rsqb;---\left(39\right)$

可知，两星地心距之差变化的幅值正比于其中θ₀是两个卫星地心矢量的夹角。考虑到(34)式，可知

$Am\left(\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\right)\&Proportional;\frac{cos\frac{{\mathrm{\&theta;}}_0}{2}}{{r_0}^2}\sqrt{\frac{\mathrm{\&mu;}}{r_0}}{\mathrm{sin\&theta;}}_0---\left(40\right)$

其中，Am()表示幅值。星间距离变化率的变化周期应当与轨道周期相同，从而

$\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\left(t\right)=K\frac{cos\frac{{\mathrm{\&theta;}}_0}{2}}{{r_0}^2}\sqrt{\frac{\mathrm{\&mu;}}{r_0}}{\mathrm{sin\&theta;}}_{0}cos(nt+\mathrm{\&xi;})---\left(41\right)$

其中，K是待定系数，n是轨道角速度，ξ＝-π/2是初始相位。考虑到卫星轨道为极轨道，所以有

$\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\left(\mathrm{\&theta;}\right)=K\frac{cos\frac{{\mathrm{\&theta;}}_0}{2}}{{r_0}^2}\sqrt{\frac{\mathrm{\&mu;}}{r_0}}{\mathrm{sin\&theta;}}_{0}cos(\mathrm{\&theta;}+\mathrm{\&xi;})---\left(42\right)$

为了得到系数K，在不同的轨道高度和星间距离下，进行纯引力轨道积分，根据积分轨道得到星间距离变化率随时间振荡的幅值，如表1所示。

表1

	0.2°	0.4°	0.6°	0.8°	1.0°
						200km	8.52384717×10-2	1.70054302×10-1	2.54878559×10-1	3.39375389×10-1	4.24052520×10-1
250km	8.39648169×10-2	1.66955167×10-1	2.49980172×10-1	3.32888860×10-1	4.15999957×10-1
						300km	8.22614068×10-2	1.63879328×10-1	2.45575028×10-1	3.26890250×10-1	4.08533929×10-1
350km	8.07749317×10-2	1.60715522×10-1	2.40566235×10-1	3.20445000×10-1	4.00451832×10-1
						400km	7.90730919×10-2	1.57561129×10-1	2.35965947×10-1	3.14446876×10-1	3.92894479×10-1
	1.2°	1.4°	1.6°	1.8°	2.0°
						200km	5.08769221×10-1	5.93621975×10-1	6.78081254×10-1	7.62879389×10-1	8.47098944×10-1
250km	4.99688820×10-1	5.82412568×10-1	6.65497074×10-1	7.48424917×10-1	8.31446183×10-1
						300km	4.90153499×10-1	5.71442505×10-1	6.53303618×10-1	7.34316987×10-1	8.15750088×10-1
350km	4.80516090×10-1	5.60716101×10-1	6.40541874×10-1	7.20322030×10-1	8.00195000×10-1
						400km	4.71103131×10-1	5.49539036×10-1	6.28048672×10-1	7.06380954×10-1	7.84346464×10-1

以为横坐标，以表1中的星间距离变化率幅值为纵坐标，得到图2所示直线，说明(42)式给出的星间距离变化率表达式是合理的，拟合得到系数K为

K＝1.3476×10¹¹m²(43)

将(42)式代入(35)式，得到低低星星跟踪的能量变分方程为

${\mathrm{\&delta;R}}_A-{\mathrm{\&delta;R}}_B=K(\frac{cos\frac{{\mathrm{\&theta;}}_0}{2}}{2}-\frac{4}{3})\frac{{\mathrm{\&mu;sin\&theta;}}_0}{{r_0}^4}cos(\mathrm{\&theta;}+\mathrm{\&xi;})\mathrm{\&delta;}r-\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}---\left(44\right)$

下面将基于(44)式，建立低低星星跟踪重力场测量的阶误差方差，进而得到重力场测量的有效阶数、大地水准面误差和重力异常误差等重力场测量性能。

4低低星星跟踪重力场测量的阶误差方差

由(44)式，得到第n阶上的阶误差方差关系为

${P_n}^2({\mathrm{\&delta;R}}_A-{\mathrm{\&delta;R}}_B)={P_n}^2\&lsqb;K(\frac{cos\frac{{\mathrm{\&theta;}}_0}{2}}{2}-\frac{4}{3})\frac{{\mathrm{\&mu;sin\&theta;}}_0}{{r_0}^4}cos(\mathrm{\&theta;}+\mathrm{\&xi;})\mathrm{\&delta;}r-\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\&rsqb;---\left(45\right)$

对于(45)式，由功率谱定义得到

${P_n}^2({\mathrm{\&delta;R}}_A-{\mathrm{\&delta;R}}_B)=\underset{k=-n}{\overset{n}{\&Sigma;}}{\&lsqb;\frac{1}{4\mathrm{\&pi;}}\underset{0}{\overset{2\mathrm{\&pi;}}{\&Integral;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}({\mathrm{\&delta;R}}_A-{\mathrm{\&delta;R}}_B){\stackrel{\&OverBar;}{Y}}_{nk}(\mathrm{\&theta;},\mathrm{\&lambda;})sin\mathrm{\&theta;}d\mathrm{\&theta;}d\mathrm{\&lambda;}\&rsqb;}^2---\left(46\right)$

其中，

${\stackrel{\&OverBar;}{Y}}_{nk}(\mathrm{\&theta;},\mathrm{\&lambda;})={\stackrel{\&OverBar;}{P}}_{n\left|k\right|}\left(cos\mathrm{\&theta;}\right)Q_k\left(\mathrm{\&lambda;}\right)---\left(47\right)$

$Q_k\left(\mathrm{\&lambda;}\right)=\left\{\begin{array}{cc}\cos k\mathrm{\&lambda;}& k\&GreaterEqual;0\\ \sin \left|k\right|\mathrm{\&lambda;}& k<0\end{array}\right.---\left(48\right)$

由于低低星星跟踪的两个卫星沿均极轨道运行，并且在同一个轨道面内，所以它们的非球形摄动引力位可以分别表示为

$R_A(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{\&mu;}}{r_0}\underset{n=2}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\underset{k=0}{\overset{n}{\&Sigma;}}{\left(\frac{a_e}{r_0}\right)}^n({\stackrel{\&OverBar;}{C}}_{nk}\cos k\mathrm{\&lambda;}+{\stackrel{\&OverBar;}{S}}_{nk}\sin k\mathrm{\&lambda;}){\stackrel{\&OverBar;}{P}}_{nk}\left(cos\right(\mathrm{\&theta;}+{\mathrm{\&theta;}}_0\left)\right)---\left(49\right)$

$R_B(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{\&mu;}}{r_0}\underset{n=2}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\underset{k=0}{\overset{n}{\&Sigma;}}{\left(\frac{a_e}{r_0}\right)}^n({\stackrel{\&OverBar;}{C}}_{nk}\cos k\mathrm{\&lambda;}+{\stackrel{\&OverBar;}{S}}_{nk}\sin k\mathrm{\&lambda;}){\stackrel{\&OverBar;}{P}}_{nk}\left(cos\mathrm{\&theta;}\right)---\left(50\right)$

将(49)和(50)式代入(46)式，得到

$\begin{array}{c}{P_n}^2({\mathrm{\&delta;R}}_A-{\mathrm{\&delta;R}}_B)\\ =\underset{k=-n}{\overset{n}{\&Sigma;}}{\&lsqb;\frac{1}{4\mathrm{\&pi;}}\underset{0}{\overset{2\mathrm{\&pi;}}{\&Integral;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\left(\frac{\mathrm{\&mu;}}{r_0}\underset{l=2}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\underset{m=0}{\overset{l}{\&Sigma;}}{\left(\frac{a_e}{r_0}\right)}^l\right(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{lm}\cos m\mathrm{\&lambda;}+\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{lm}\sin m\mathrm{\&lambda;})\&lsqb;{\stackrel{\&OverBar;}{P}}_{lm}(\cos \left(\mathrm{\&theta;}+{\mathrm{\&theta;}}_0\right))-{\stackrel{\&OverBar;}{P}}_{lm}(\cos \mathrm{\&theta;})\&rsqb;){\stackrel{\&OverBar;}{Y}}_{nk}(\mathrm{\&theta;},\mathrm{\&lambda;})\sin \mathrm{\&theta;}d\mathrm{\&theta;}d\mathrm{\&lambda;}\&rsqb;}^2\\ \&ap;\frac{1}{16}\frac{{\mathrm{\&mu;}}^2}{{r_0}^2}\underset{k=0}{\overset{n}{\&Sigma;}}{\&lsqb;\underset{l=1,n,n+1}{\&Sigma;}\mathrm{\&delta;}\left(k\right)\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{lk}\right){\left(\frac{a_e}{r_0}\right)}^l\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\&lsqb;{\stackrel{\&OverBar;}{P}}_{lk}\left(\cos \right(\mathrm{\&theta;}+{\mathrm{\&theta;}}_0\left)\right)-{\stackrel{\&OverBar;}{P}}_{lk}\left(\cos \mathrm{\&theta;}\right)\&rsqb;{\stackrel{\&OverBar;}{P}}_{nk}\left(\cos \mathrm{\&theta;}\right)\sin \mathrm{\&theta;}d\mathrm{\&theta;}\&rsqb;}^2\\ +\frac{1}{16}\frac{{\mathrm{\&mu;}}^2}{{r_0}^2}\underset{k=1}{\overset{n}{\&Sigma;}}{\&lsqb;\underset{l=n-1,n,n+1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{lk}\right){\left(\frac{a_e}{r_0}\right)}^l\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\&lsqb;{\stackrel{\&OverBar;}{P}}_{lk}\left(\cos \right(\mathrm{\&theta;}+{\mathrm{\&theta;}}_0\left)\right)-{\stackrel{\&OverBar;}{P}}_{lk}\left(\cos \mathrm{\&theta;}\right)\&rsqb;{\stackrel{\&OverBar;}{P}}_{nk}\left(\cos \mathrm{\&theta;}\right)\sin \mathrm{\&theta;}d\mathrm{\&theta;}\&rsqb;}^2\end{array}---\left(51\right)$

其中，

${\mathrm{\&delta;}}_k=\left\{\begin{array}{cc}2,& k=0\\ 1,& k\&NotEqual;0\end{array}\right.---\left(52\right)$

设

$A_n(l,k,{\mathrm{\&theta;}}_0)={\mathrm{\&delta;}}_k\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\&lsqb;{\stackrel{\&OverBar;}{P}}_{lk}\left(cos\right(\mathrm{\&theta;}+{\mathrm{\&theta;}}_0\left)\right)-{\stackrel{\&OverBar;}{P}}_{lk}\left(cos\mathrm{\&theta;}\right)\&rsqb;{\stackrel{\&OverBar;}{P}}_{nk}\left(cos\mathrm{\&theta;}\right)sin\mathrm{\&theta;}d\mathrm{\&theta;}---\left(53\right)$

则(51)式变为

$\begin{array}{c}{P_n}^2({\mathrm{\&delta;R}}_A-{\mathrm{\&delta;R}}_B)\\ =\frac{1}{16}\frac{{\mathrm{\&mu;}}^2}{{r_0}^2}{\left(\frac{a_e}{r_0}\right)}^{2(n-1)}\underset{k=0}{\overset{n}{\&Sigma;}}{\&lsqb;\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{n-1,k}\right)A_n(n-1,k,{\mathrm{\&theta;}}_0)+\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{nk}\right)\left(\frac{a_e}{r_0}\right)A_n(n,k,{\mathrm{\&theta;}}_0)+\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{n+1,k}\right){\left(\frac{a_e}{r_0}\right)}^2{A}_n(n+1,k,{\mathrm{\&theta;}}_0)\&rsqb;}^2\\ +\frac{1}{16}\frac{{\mathrm{\&mu;}}^2}{{r_0}^2}{\left(\frac{a_e}{r_0}\right)}^{2(n-1)}\underset{k=1}{\overset{n}{\&Sigma;}}{\&lsqb;\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{n-1,k}\right)A_n(n-1,k,{\mathrm{\&theta;}}_0)+\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{nk}\right)\left(\frac{a_e}{r_0}\right)A_n(n,k,{\mathrm{\&theta;}}_0)+\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{n+1,k}\right){\left(\frac{a_e}{r_0}\right)}^2{A}_n(n+1,k,{\mathrm{\&theta;}}_0)\&rsqb;}^2\end{array}---\left(54\right)$

下面计算(54)式中第一个中括号中的表达式，将平方项展开并利用不等式关系，可得

$\begin{array}{c}\underset{k=0}{\overset{n}{\&Sigma;}}{\&lsqb;\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{n-1,k}\right)A_n(n-1,k,{\mathrm{\&theta;}}_0)+\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{nk}\right)+\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{nk}\right)\left(\frac{a_e}{r_0}\right)A_n(n,k,{\mathrm{\&theta;}}_0)+\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{n+1,k}\right){\left(\frac{a_e}{r_0}\right)}^2{A}_n(n+1,k,{\mathrm{\&theta;}}_0)\&rsqb;}^2\\ \&le;\underset{k=0}{\overset{n}{\&Sigma;}}\left[\begin{array}{c}{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{n-1,k}\right)}^2\&lsqb;{A_n}^2(n-1,k,{\mathrm{\&theta;}}_0)+\left(\frac{a_e}{r_0}\right)\left|A_n(n-1,k,{\mathrm{\&theta;}}_0)A_n(n,k,{\mathrm{\&theta;}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^2\left|A_n(n-1,k,{\mathrm{\&theta;}}_0)A_n(n+1,k,{\mathrm{\&theta;}}_0)\right|\&rsqb;\\ +{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{nk}\right)}^2\&lsqb;\left(\frac{a_e}{r_0}\right){A_n}^2(n,k,{\mathrm{\&theta;}}_0)+\left(\frac{a_e}{r_0}\right)\left|A_n(n-1,k,{\mathrm{\&theta;}}_0)A_n(n,k,{\mathrm{\&theta;}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^3\left|A_n(n,k,{\mathrm{\&theta;}}_0)A_n(n+1,k,{\mathrm{\&theta;}}_0)\right|\&rsqb;\\ {\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{n+1,k}\right)}^2\&lsqb;{\left(\frac{a_e}{r_0}\right)}^4{A_n}^2(n+1,k,{\mathrm{\&theta;}}_0)+{\left(\frac{a_e}{r_0}\right)}^2\left|A_n(n-1,k,{\mathrm{\&theta;}}_0)A_n(n+1,k,{\mathrm{\&theta;}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^3\left|A_n(n,k,{\mathrm{\&theta;}}_0)A_n(n+1,k,{\mathrm{\&theta;}}_0)\right|\&rsqb;\end{array}\right]\end{array}---\left(55\right)$

设

$\left\{\begin{array}{c}{B}_1(r_0,n,k,{\mathrm{\&theta;}}_0)={A_n}^2(n-1,k,{\mathrm{\&theta;}}_0)+\left(\frac{a_e}{r_0}\right)\left|A_n(n-1,k,{\mathrm{\&theta;}}_0)A_n(n,k,{\mathrm{\&theta;}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^2\left|A_n(n-1,k,{\mathrm{\&theta;}}_0)A_n(n+1,k,{\mathrm{\&theta;}}_0)\right|\\ B_2(r_0,n,k,{\mathrm{\&theta;}}_0)={\left(\frac{a_e}{r_0}\right)}^2{A_n}^2(n,k,{\mathrm{\&theta;}}_0)+\left(\frac{a_e}{r_0}\right)\left|A_n(n-1,k,{\mathrm{\&theta;}}_0)A_n(n,k,{\mathrm{\&theta;}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^3\left|A_n(n,k,{\mathrm{\&theta;}}_0)A_n(n+1,k,{\mathrm{\&theta;}}_0)\right|\\ B_3(r_0,n,k,{\mathrm{\&theta;}}_0)={\left(\frac{a_e}{r_0}\right)}^4{A_n}^2(n+1,k,{\mathrm{\&theta;}}_0)+{\left(\frac{a_e}{r_0}\right)}^2\left|A_n(n-1,k,{\mathrm{\&theta;}}_0)A_n(n+1,k,{\mathrm{\&theta;}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^3\left|A_n(n,k,{\mathrm{\&theta;}}_0)A_n(n+1,k,{\mathrm{\&theta;}}_0)\right|\end{array}\right.---\left(56\right)$

则(55)式可以表示为

$\begin{array}{c}\underset{k=0}{\overset{n}{\&Sigma;}}{\&lsqb;\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{n-1,k}\right)A_n(n-1,k,{\mathrm{\&theta;}}_0)+\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{nk}\right)\left(\frac{a_e}{r_0}\right)A_n(n,k,{\mathrm{\&theta;}}_0)+\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{n+1,k}\right){\left(\frac{a_e}{r_0}\right)}^2{A}_n(n+1,k,{\mathrm{\&theta;}}_0)\&rsqb;}^2\\ \&le;\underset{k=0}{\overset{n}{\&Sigma;}}\&lsqb;{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{n-1,k}\right)}^2{B}_1(r_1,n,k,{\mathrm{\&theta;}}_0)+{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{nk}\right)}^2{B}_2(r_0,n,k,{\mathrm{\&theta;}}_0)+{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{n+1,k}\right)}^2{B}_3(r_0,n,k,{\mathrm{\&theta;}}_0)\&rsqb;\\ \&ap;\underset{k=0}{\overset{n}{\&Sigma;}}{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{nk}\right)}^2\&CenterDot;\frac{1}{n+1}\underset{k=0}{\overset{n}{\&Sigma;}}\&lsqb;B_1(r_0,n,k,{\mathrm{\&theta;}}_0)+B_2(r_0,n,k,{\mathrm{\&theta;}}_0)+B_3(r_0,n,k,{\mathrm{\&theta;}}_0)\&rsqb;\end{array}---\left(57\right)$

同理，得到(54)式中第2个中括号的表达式为

$\begin{array}{c}\underset{k=1}{\overset{n}{\&Sigma;}}{\&lsqb;\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{n-1,k}\right)A_n(n-1,k,{\mathrm{\&theta;}}_0)+\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{nk}\right)\left(\frac{a_e}{r_0}\right)A_n(n,k,{\mathrm{\&theta;}}_0)+\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{n+1,k}\right){\left(\frac{a_e}{r_0}\right)}^2{A}_n(n+1,k,{\mathrm{\&theta;}}_0)\&rsqb;}^2\\ \&ap;\underset{k=1}{\overset{n}{\&Sigma;}}{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{nk}\right)}^2\&CenterDot;\frac{1}{n+1}\underset{k=0}{\overset{n}{\&Sigma;}}\&lsqb;B_1(r_0,n,k,{\mathrm{\&theta;}}_0)+B_2(r_0,n,k,{\mathrm{\&theta;}}_0)+B_3(r_0,n,k,{\mathrm{\&theta;}}_0)\&rsqb;\end{array}---\left(58\right)$

将(57)和(58)式代入(54)式，得到

$\begin{array}{c}{P_n}^2({\mathrm{\&delta;R}}_A-{\mathrm{\&delta;R}}_B)\\ =\frac{1}{16(n+1)}\frac{{\mathrm{\&mu;}}^2}{{r_0}^2}{\left(\frac{a_e}{r_0}\right)}^{2(n-1)}\underset{k=0}{\overset{n}{\&Sigma;}}\&lsqb;B_1(r_0,n,k,{\mathrm{\&theta;}}_0)+B_2(r_0,n,k,{\mathrm{\&theta;}}_0)+B_3(r_0,n,k,{\mathrm{\&theta;}}_0)\&rsqb;\&CenterDot;\underset{k=0}{\overset{n}{\&Sigma;}}\&lsqb;{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{nk}\right)}^2+{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{nk}\right)}^2\&rsqb;\end{array}---\left(59\right)$

下面推导(45)式的右端部分。设

$D=K(\frac{cos\frac{{\mathrm{\&theta;}}_0}{2}}{2}-\frac{4}{3})\frac{{\mathrm{\&mu;sin\&theta;}}_0}{{r_0}^4}---\left(60\right)$

根据功率谱密度的定义，得到(45)式右端部分为

$\begin{array}{c}{P_n}^2\&lsqb;D\cos (\mathrm{\&theta;}+\mathrm{\&xi;})\mathrm{\&delta;}r-\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\&rsqb;\\ =\underset{k=-n}{\overset{n}{\&Sigma;}}{\left\{\frac{D}{4\mathrm{\&pi;}}\underset{0}{\overset{2\mathrm{\&pi;}}{\&Integral;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\cos (\mathrm{\&theta;}+\mathrm{\&xi;})\mathrm{\&delta;}r{\stackrel{\&OverBar;}{Y}}_{nk}(\mathrm{\&theta;},\mathrm{\&lambda;})\sin \mathrm{\&theta;}d\mathrm{\&theta;}d\mathrm{\&lambda;}\right\}}^2+\underset{k=-n}{\overset{n}{\&Sigma;}}{\left\{\frac{1}{4\mathrm{\&pi;}}\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\underset{0}{\overset{2\mathrm{\&pi;}}{\&Integral;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}{\stackrel{\&OverBar;}{Y}}_{nk}(\mathrm{\&theta;},\mathrm{\&lambda;})\sin \mathrm{\&theta;}d\mathrm{\&theta;}d\mathrm{\&lambda;}\right\}}^2\\ -\underset{k=-n}{\overset{n}{\&Sigma;}}\{2\&CenterDot;\frac{D}{4\mathrm{\&pi;}}\underset{0}{\overset{2\mathrm{\&pi;}}{\&Integral;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\cos (\mathrm{\&theta;}+\mathrm{\&xi;})\mathrm{\&delta;}r{\stackrel{\&OverBar;}{Y}}_{nk}(\mathrm{\&theta;},\mathrm{\&lambda;})\sin \mathrm{\&theta;}d\mathrm{\&theta;}d\mathrm{\&lambda;}\&CenterDot;\frac{1}{4\mathrm{\&pi;}}\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\underset{0}{\overset{2\mathrm{\&pi;}}{\&Integral;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}{\stackrel{\&OverBar;}{Y}}_{nk}(\mathrm{\&theta;},\mathrm{\&lambda;})\sin \mathrm{\&theta;}d\mathrm{\&theta;}d\mathrm{\&lambda;}\}\end{array}---\left(61\right)$

(61)式中的各项分别为

$\begin{array}{c}\underset{k=-n}{\overset{n}{\&Sigma;}}{\left\{\frac{D}{4\mathrm{\&pi;}}\underset{0}{\overset{2\mathrm{\&pi;}}{\&Integral;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\cos (\mathrm{\&theta;}+\mathrm{\&xi;})\mathrm{\&delta;}r{\stackrel{\&OverBar;}{Y}}_{nk}(\mathrm{\&theta;},\mathrm{\&lambda;})\sin \mathrm{\&theta;}d\mathrm{\&theta;}d\mathrm{\&lambda;}\right\}}^2\\ ={\left(\frac{D}{4\mathrm{\&pi;}}\right)}^2\frac{2{\mathrm{\&pi;S}}_{\mathrm{\&delta;}r}}{T}\underset{k=0}{\overset{n}{\&Sigma;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}{\&lsqb;{\stackrel{\&OverBar;}{P}}_{nk}\left(\cos \mathrm{\&theta;}\right)\&rsqb;}^2{\cos }^2(\mathrm{\&theta;}+\mathrm{\&xi;}){\sin }^2\mathrm{\&theta;}d\mathrm{\&theta;}\end{array}---\left(62\right)$

$\begin{array}{c}\underset{k=-n}{\overset{n}{\&Sigma;}}{\left\{\frac{1}{4\mathrm{\&pi;}}\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\underset{0}{\overset{2\mathrm{\&pi;}}{\&Integral;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}{\stackrel{\&OverBar;}{Y}}_{nk}(\mathrm{\&theta;},\mathrm{\&lambda;})\sin \mathrm{\&theta;}d\mathrm{\&theta;}d\mathrm{\&lambda;}\right\}}^2\\ ={\left(\frac{1}{4\mathrm{\&pi;}}\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\right)}^2\frac{2{\mathrm{\&pi;S}}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}}{T}\underset{k=0}{\overset{n}{\&Sigma;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}{\&lsqb;{\stackrel{\&OverBar;}{P}}_{nk}\left(\cos \mathrm{\&theta;}\right)\&rsqb;}^2{\sin }^2\mathrm{\&theta;}d\mathrm{\&theta;}\end{array}---\left(63\right)$

其中，S_δr是纯引力轨道位置误差δr的功率谱密度，是纯引力作用下的星间距离变化率误差的功率谱密度，T是重力场测量总时间。为便于表示，在(62)和(63)式中设

$f_{\mathrm{\&delta;}r}\left(n\right)=\underset{k=0}{\overset{n}{\&Sigma;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}{\&lsqb;{\stackrel{\&OverBar;}{P}}_{nk}\left(cos\mathrm{\&theta;}\right)\&rsqb;}^2{\cos }^2(\mathrm{\&theta;}+\mathrm{\&xi;}){\sin }^2\mathrm{\&theta;}d\mathrm{\&theta;}---\left(64\right)$

$f_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}\left(n\right)=\underset{k=0}{\overset{n}{\&Sigma;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}{\&lsqb;{\stackrel{\&OverBar;}{P}}_{nk}\left(cos\mathrm{\&theta;}\right)\&rsqb;}^2{\sin }^2\mathrm{\&theta;}d\mathrm{\&theta;}---\left(65\right)$

由(61)～(65)式，得到

${P_n}^2\&lsqb;Dcos(\mathrm{\&theta;}+\mathrm{\&xi;})\mathrm{\&delta;}r-\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\&rsqb;\&le;{\&lsqb;\frac{-D}{4\mathrm{\&pi;}}\sqrt{\frac{2{\mathrm{\&pi;S}}_{\mathrm{\&delta;}r}}{T}{f}_{\mathrm{\&delta;}r}\left(n\right)}+\frac{1}{4\mathrm{\&pi;}}\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\sqrt{\frac{2{\mathrm{\&pi;S}}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}}{T}{f}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}\left(n\right)}\&rsqb;}^2---\left(66\right)$

已知纯引力轨道位置误差δr和纯引力星间距离变化率误差的方差与功率谱密度的关系为

${{\mathrm{\&sigma;}}_{\mathrm{\&delta;}r}}^2=\underset{-f_{\mathrm{\&delta;}r,\max }}{\overset{f_{\mathrm{\&delta;}r,\max }}{\&Integral;}}{S}_{\mathrm{\&delta;}r}\left(f\right)df=2S_{\mathrm{\&delta;}r}{f}_{\mathrm{\&delta;}r,max}---\left(67\right)$

${{\mathrm{\&sigma;}}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}}^2=\underset{-f_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}},max}}{\overset{f_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}},max}}{\&Integral;}}{S}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}\left(f\right)df=2S_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}{f}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}},max}---\left(68\right)$

其中，最高频率f_δr,max和通常取为Nyquist频率

$f_{\mathrm{\&delta;}r,max}=\frac{1}{2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&delta;}r}}---\left(69\right)$

$f_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}},max}=\frac{1}{2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}}---\left(70\right)$

其中，(△t)_δr和分别是纯引力轨道位置误差和纯引力星间距离变化率的数据采样间隔。将(67)～(70)式代入(66)式，得到

$\begin{array}{c}{P_n}^2\&lsqb;D\cos (\mathrm{\&theta;}+\mathrm{\&xi;})\mathrm{\&delta;}r-\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\&rsqb;\\ \&le;{\&lsqb;\frac{-D}{4\mathrm{\&pi;}}\sqrt{\frac{2{{\mathrm{\&pi;\&sigma;}}_{\mathrm{\&delta;}r}}^2\left(\mathrm{\&Delta;}t\right)}{T}{f}_{\mathrm{\&delta;}r}\left(n\right)}+\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\sqrt{\frac{2{{\mathrm{\&pi;\&sigma;}}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}}^2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}}{T}{f}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}\left(n\right)}\&rsqb;}^2\\ =\frac{1}{8\mathrm{\&pi;}T}{\&lsqb;-D\sqrt{{{\mathrm{\&sigma;}}_{\mathrm{\&delta;}r}}^2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&delta;}r}{f}_{\mathrm{\&delta;}r}\left(n\right)}+\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\sqrt{{{\mathrm{\&sigma;}}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}}^2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}\left(n\right)}\&rsqb;}^2\end{array}---\left(71\right)$

由(45)、(59)和(71)式，得到低低星星跟踪重力场测量的阶误差方差为

$\begin{array}{c}{{\mathrm{\&delta;\&sigma;}}_n}^2=\frac{1}{2n+1}\underset{k=0}{\overset{n}{\&Sigma;}}\&lsqb;{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{nk}\right)}^2+{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{nk}\right)}^2\&rsqb;=\frac{1}{\underset{k=0}{\overset{n}{\&Sigma;}}\&lsqb;B_1(r_0,n,k,{\mathrm{\&theta;}}_0)+B_2(r_0,n,k,{\mathrm{\&theta;}}_0)+B_3(r_0,n,k,{\mathrm{\&theta;}}_0)\&rsqb;}\&times;\\ \frac{2(n+1)}{2n+1}\frac{{r_2}^{2n}}{{\mathrm{\&pi;\&mu;}}^2{{\mathrm{Ta}}_e}^{2n-2}}{\&lsqb;-D+\sqrt{{{\mathrm{\&sigma;}}_{\mathrm{\&delta;}r}}^2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&delta;}r}{f}_{{\mathrm{\&delta;}}_r}\left(n\right)}+\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\sqrt{{{\mathrm{\&sigma;}}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}}^2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}{f}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}\left(n\right)}\&rsqb;}^2\end{array}---\left(72\right)$

在(72)式中，纯引力轨道位置误差σ_δr ²包括纯引力轨道定轨误差(△r)_m ²和非引力干扰引起的纯引力轨道偏移(△r)_△F ²。同样，纯引力作用下的星间距离变化率误差包括星间距离变化率测量误差和非引力干扰引起的星间距离变化率偏移即

σ_δr ²(△t)_δr＝(△r)_△F ²(△t)_△F+(△r)_m ²(△t)_△r(73)

${{\mathrm{\&sigma;}}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}}^2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}={{\left(\mathrm{\&Delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\right)}_{\mathrm{\&Delta;}F}}^2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&Delta;}F}+{{\left(\mathrm{\&Delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\right)}_m}^2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&Delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}---\left(74\right)$

其中，(△t)_△F是非引力干扰数据间隔，(△t)_△r是卫星轨道数据的采样间隔，是星间距离变化率数据的采样间隔。需要说明的是，如果低低星星跟踪重力卫星系统对非引力干扰进行测量，那么(△t)_△F指非引力干扰的测量间隔；如果低低星星跟踪重力场测量系统对非引力干扰进行抑制，那么(△t)_△F指非引力干扰的抑制间隔。

非引力干扰δF引起纯引力轨道位置累积误差和星间距离变化率累积误差，其最大值对应匀加速直线运动条件下的累积误差，其平均累积误差为

${\left(\mathrm{\&Delta;}r\right)}_{\mathrm{\&Delta;}F}=\frac{1}{T_{arc}}\underset{0}{\overset{T_{arc}}{\&Integral;}}\frac{1}{2}\left(\mathrm{\&Delta;}F\right)t^{2}dt=\frac{\left(\mathrm{\&Delta;}F\right){T_{arc}}^2}{6}---\left(75\right)$

${\left(\mathrm{\&Delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\right)}_{\mathrm{\&Delta;}F}=\frac{1}{T_{arc}}\underset{0}{\overset{T_{arc}}{\&Integral;}}\left(\mathrm{\&Delta;}F\right)tdt=\frac{\left(\mathrm{\&Delta;}F\right)T_{arc}}{2}---\left(76\right)$

其中，T_arc为重力场反演中的积分弧长。将式(73)～(76)代入(72)式，得到阶误差方差为

$\begin{array}{c}{{\mathrm{\&delta;\&sigma;}}_n}^2=\frac{1}{2n+1}\underset{k=0}{\overset{n}{\&Sigma;}}\&lsqb;{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{nk}\right)}^2+{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{nk}\right)}^2\&rsqb;=\frac{1}{\underset{k=0}{\overset{n}{\&Sigma;}}\&lsqb;B_1(r_0,n,k,{\mathrm{\&theta;}}_0)+B_2(r_0,n,k,{\mathrm{\&theta;}}_0)+B_3(r_0,n,k,{\mathrm{\&theta;}}_0)\&rsqb;}\&times;\\ \frac{1}{2n+1}\frac{2(n+1){r_0}^{2n}}{{\mathrm{\&pi;\&mu;}}^2{{\mathrm{Ta}}_e}^{2n-2}}{\&lsqb;-D\sqrt{(\frac{{\left(\mathrm{\&Delta;}F\right)}^2{T_{arc}}^2}{36}{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&Delta;}F}+{{\left(\mathrm{\&Delta;}r\right)}_m}^2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&Delta;}r})f_{\mathrm{\&delta;}r}\left(n\right)}+\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\sqrt{(\frac{{\left(\mathrm{\&Delta;}F\right)}^2{T_{arc}}^2}{4}{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&Delta;}F}+{{\left(\mathrm{\&Delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\right)}_m}^2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&Delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}})f_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}\left(n\right)}\&rsqb;}^2\end{array}---\left(77\right)$

其中，

$\left\{\begin{array}{c}{B}_1(r_0,n,k,{\mathrm{\&theta;}}_0)={A_n}^2(n-1,k,{\mathrm{\&theta;}}_0)+\left(\frac{a_e}{r_0}\right)\left|A_n(n-1,k,{\mathrm{\&theta;}}_0)A_n(n,k,{\mathrm{\&theta;}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^2\left|A_n(n-1,k,{\mathrm{\&theta;}}_0)A_n(n+1,k,{\mathrm{\&theta;}}_0)\right|\\ B_2(r_0,n,k,{\mathrm{\&theta;}}_0)={\left(\frac{a_e}{r_0}\right)}^2{A_n}^2(n,k,{\mathrm{\&theta;}}_0)+\left(\frac{a_e}{r_0}\right)\left|A_n(n-1,k,{\mathrm{\&theta;}}_0)A_n(n,k,{\mathrm{\&theta;}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^3\left|A_n(n,k,{\mathrm{\&theta;}}_0)A_n(n+1,k,{\mathrm{\&theta;}}_0)\right|\\ B_3(r_0,n,k,{\mathrm{\&theta;}}_0)={\left(\frac{a_e}{r_0}\right)}^4{A_n}^2(n+1,k,{\mathrm{\&theta;}}_0)+{\left(\frac{a_e}{r_0}\right)}^2\left|A_n(n-1,k,{\mathrm{\&theta;}}_0)A_n(n+1,k,{\mathrm{\&theta;}}_0)\right|+{\left(\frac{a_e}{r_0}\right)}^3\left|A_n(n,k,{\mathrm{\&theta;}}_0)A_n(n+1,k,{\mathrm{\&theta;}}_0)\right|\end{array}\right.---\left(78\right)$

$A_n(l,k,{\mathrm{\&theta;}}_0)={\mathrm{\&delta;}}_k\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\&lsqb;{\stackrel{\&OverBar;}{P}}_{lk}\left(cos\right(\mathrm{\&theta;}+{\mathrm{\&theta;}}_0\left)\right)-{\stackrel{\&OverBar;}{P}}_{lk}\left(cos\mathrm{\&theta;}\right)\&rsqb;{\stackrel{\&OverBar;}{P}}_{nk}\left(cos\mathrm{\&theta;}\right)sin\mathrm{\&theta;}d\mathrm{\&theta;}---\left(79\right)$

${\mathrm{\&delta;}}_k=\{\begin{array}{cc}2,& k=0\\ 1,& k\&NotEqual;0\end{array}---\left(80\right)$

$D=K(\frac{cos\frac{{\mathrm{\&theta;}}_0}{2}}{2}-\frac{4}{3})\frac{{\mathrm{\&mu;sin\&theta;}}_0}{{r_0}^4}<0---\left(81\right)$

$\begin{array}{c}{f}_{\mathrm{\&delta;}r}\left(n\right)=\underset{k=0}{\overset{n}{\&Sigma;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}{\&lsqb;{\stackrel{\&OverBar;}{P}}_{nk}\left(\cos \mathrm{\&theta;}\right)\&rsqb;}^2{\cos }^2(\mathrm{\&theta;}+\mathrm{\&xi;}){\sin }^2\mathrm{\&theta;}d\mathrm{\&theta;}=2.3562n+1.1781\\ =\frac{3}{8}(2n+1)\mathrm{\&pi;}\end{array}---\left(82\right)$

$f_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}\left(n\right)=\underset{k=0}{\overset{n}{\&Sigma;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}{\&lsqb;{\stackrel{\&OverBar;}{P}}_{nk}\left(cos\mathrm{\&theta;}\right)\&rsqb;}^2{\sin }^2\mathrm{\&theta;}d\mathrm{\&theta;}=(n+\frac{1}{2})\mathrm{\&pi;}---\left(83\right)$

r₀＝a_e+h(84)

$K=1.3476\&times;{10}^{11}{m}^2,\mathrm{\&xi;}=-\frac{\mathrm{\&pi;}}{2}---\left(85\right)$

上述公式中的物理参数说明如表2所示。

表2

5低低星星跟踪重力场测量性能计算

由阶误差方差可以确定重力场测量的有效阶数、大地水准面误差和重力异常误差。根据Kaula准则，地球重力场模型的阶方差为

${\mathrm{\&sigma;}}_n^2=\frac{1}{2n+1}\underset{k=0}{\overset{n}{\&Sigma;}}({{\stackrel{\&OverBar;}{C}}_{nk}}^2+{{\stackrel{\&OverBar;}{S}}_{nk}}^2)\&ap;\frac{1}{2n+1}\frac{1.6\&times;{10}^{-10}}{n^3}---\left(86\right)$

阶误差方差δσ_n ²是重力场模型阶数n的增函数，阶方差σ_n ²是n的减函数。随着n的增加，当δσ_n ²等于σ_n ²时，认为达到重力场反演的最高有效阶数N_max，即

${{\mathrm{\&delta;\&sigma;}}_{N_{max}}}^2=\frac{1}{2N_{max}+1}\frac{1.6\&times;{10}^{-10}}{{N_{max}}^3}---\left(87\right)$

由(87)式得到n阶大地水准面的阶误差为

$n=a_e\sqrt{(2n+1){{\mathrm{\&delta;\&sigma;}}_n}^2}---\left(88\right)$

进而得到n阶大地水准面的累积误差为

$\mathrm{\&Delta;}=\sqrt{\underset{k=2}{\overset{n}{\&Sigma;}}{\left({\mathrm{\&Delta;}}_k\right)}^2}---\left(89\right)$

由(87)式得到n阶重力异常的阶误差为

${\mathrm{\&Delta;g}}_n=\frac{\mathrm{\&mu;}}{{a_e}^2}(n-1)\sqrt{(2n+1){{\mathrm{\&delta;\&sigma;}}_n}^2}---\left(90\right)$

进而得到n阶重力异常的累积误差为

$\mathrm{\&Delta;}g=\sqrt{\underset{k=2}{\overset{n}{\&Sigma;}}{\left({\mathrm{\&Delta;g}}_k\right)}^2}---\left(91\right)$

6低低星星跟踪重力场测量性能解析分析方法的验证

由于GRACE重力卫星采用了低低星星跟踪重力场测量原理，所以可以根据GRACE卫星的系统参数，利用本发明“一种低低星星跟踪卫星重力场测量性能解析计算方法”计算GRACE重力场测量性能，然后与GRACE重力场测量数值模拟结果对比，可以验证本发明的正确性。

已知GRACE重力卫星系统参数如表3所示，利用本发明给出的低低星星跟踪重力场测量性能解析分析方法计算得到重力场测量的阶误差方差、大地水准面误差和重力异常误差，如图3～5所示。可知，基于本发明得到的GRACE重力场测量有效阶数为146，对应的大地水准面累积误差是11.73cm，重力异常累积误差是2.51mGal，与GRACE重力场测量性能数值模拟结果基本吻合，从而验证了本发明“适用于低低星星跟踪原理的卫星重力场测量性能解析计算方法”的正确性。

表3

以上通过具体的和优选的实施例详细的描述了本发明，但本领域技术人员应该明白，本发明并不局限于以上所述实施例，凡在本发明的精神和原则之内，所作的任何修改、等同替换等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种低低星星跟踪卫星重力场测量性能解析计算方法，其特征在于,包括步骤如下：

步骤1：获取低低星星跟踪重力卫星系统参数；

步骤2：根据低低星星跟踪重力卫星系统参数，计算重力卫星载荷测量误差对地球引力非球形摄动位功率谱的影响，进而得到反演重力场模型的位系数阶误差方差；

步骤3：将得到的反演重力场模型位系数阶误差方差与Kaula准则给出的位系数阶方差比较；随着重力场模型阶数的增加，位系数阶误差方差逐渐增加，而位系数阶方差则逐渐减小，当阶误差方差等于阶方差时，认为达到重力场测量的最高有效阶数；

步骤4：根据反演重力场模型的阶误差方差，计算反演重力场模型的大地水准面阶误差及其累积误差、重力异常阶误差及其累积误差；

步骤5：将计算得到的重力场测量有效阶数、大地水准面阶误差及其累积误差、重力异常阶误差及其累积误差汇总，即为卫星重力场测量性能指标；

步骤1中的所述低低星星跟踪重力卫星系统参数，包括但不限于重力卫星系统轨道参数和所述重力卫星系统载荷指标；

所述卫星重力场测量性能指标包括重力场反演的最高有效阶数N_max、n阶对应的大地水准面阶误差Δ_n、n阶对应的大地水准面累积误差Δ、n阶对应的重力异常误差Δg_n和n阶对应的重力异常累积误差Δg中的一种或几种；

所述重力卫星系统轨道参数，包括重力卫星轨道高度h和两星的地心矢量夹角θ₀；

所述重力卫星系统载荷指标，包括星间距离变化率测量误差卫星定轨位置误差(Δr)_m、非引力干扰ΔF、星间距离变化率数据采样间隔卫星轨道位置数据采样间隔(Δt)_Δr、非引力干扰数据间隔(Δt)_ΔF和重力场测量任务寿命T；

步骤2具体包括：

建立关于低低星星跟踪重力场测量位系数阶误差方差δσ_n ²满足的解析关系式：

$\begin{array}{c}{{\mathrm{\&delta;\&sigma;}}_n}^2=\frac{1}{2n+1}\underset{k=0}{\overset{n}{\mathrm{\&Sigma;}}}\&lsqb;{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{nk}\right)}^2+{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{nk}\right)}^2\&rsqb;=\frac{1}{\underset{k=0}{\overset{n}{\mathrm{\&Sigma;}}}\&lsqb;B_1\left(r_0,n,k,{\mathrm{\&theta;}}_0\right)+B_2\left(r_0,n,k,{\mathrm{\&theta;}}_0\right)+B_3\left(r_0,n,k,{\mathrm{\&theta;}}_0\right)\&rsqb;}\&times;\\ \frac{1}{2n+1}\frac{2(n+1){r_0}^{2n}}{{\mathrm{\&pi;\&mu;}}^2{{\mathrm{Ta}}_e}^{2n-2}}{\&lsqb;-D\sqrt{\left(\frac{{\left(\mathrm{\&Delta;}F\right)}^2{T_{arc}}^4}{36}{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&Delta;}F}+{{\left(\mathrm{\&Delta;}r\right)}_m}^2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&Delta;}r}\right)f_{\mathrm{\&delta;}r}\left(n\right)}+\sqrt{\frac{\mathrm{\&mu;}}{r_0}}\sqrt{\left(\frac{{\left(\mathrm{\&Delta;}F\right)}^2{T_{arc}}^2}{4}{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&Delta;}F}+{{\left(\mathrm{\&Delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\right)}_m}^2{\left(\mathrm{\&Delta;}t\right)}_{\mathrm{\&Delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}\right)f_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}\left(n\right)}\&rsqb;}^2\end{array}$

其中，

$\left\{\begin{array}{c}{B}_1\left(r_0,n,k,{\mathrm{\&theta;}}_0\right)={A_n}^2\left(n-1,k,{\mathrm{\&theta;}}_0\right)+\left(\frac{a_e}{r_0}\right)\left|A_n\left(n-1,k,{\mathrm{\&theta;}}_0\right)A_n\left(n,k,{\mathrm{\&theta;}}_0\right)\right|+{\left(\frac{a_e}{r_0}\right)}^2\left|A_n\left(n-1,k,{\mathrm{\&theta;}}_0\right)A_n\left(n+1,k,{\mathrm{\&theta;}}_0\right)\right|\\ B_2\left(r_0,n,k,{\mathrm{\&theta;}}_0\right)={\left(\frac{a_e}{r_0}\right)}^2{A_n}^2\left(n,k,{\mathrm{\&theta;}}_0\right)+\left(\frac{a_e}{r_0}\right)\left|A_n\left(n-1,k,{\mathrm{\&theta;}}_0\right)A_n\left(n,k,{\mathrm{\&theta;}}_0\right)\right|+{\left(\frac{a_e}{r_0}\right)}^3\left|A_n\left(n,k,{\mathrm{\&theta;}}_0\right)A_n\left(n+1,k,{\mathrm{\&theta;}}_0\right)\right|\\ B_3\left(r_0,n,k,{\mathrm{\&theta;}}_0\right)={\left(\frac{a_e}{r_0}\right)}^4{A_n}^2\left(n+1,k,{\mathrm{\&theta;}}_0\right)+{\left(\frac{a_e}{r_0}\right)}^2\left|A_n\left(n-1,k,{\mathrm{\&theta;}}_0\right)A_n\left(n+1,k,{\mathrm{\&theta;}}_0\right)\right|+{\left(\frac{a_e}{r_0}\right)}^3\left|A_n\left(n,k,{\mathrm{\&theta;}}_0\right)A_n\left(n+1,k,{\mathrm{\&theta;}}_0\right)\right|\end{array}\right.$

$A_n(l,k,{\mathrm{\&theta;}}_0)={\mathrm{\&delta;}}_k\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\&lsqb;{\stackrel{\&OverBar;}{P}}_{lk}\left(cos\left(\mathrm{\&theta;}+{\mathrm{\&theta;}}_0\right)\right)-{\stackrel{\&OverBar;}{P}}_{lk}\left(cos\mathrm{\&theta;}\right)\&rsqb;{\stackrel{\&OverBar;}{P}}_{nk}\left(cos\mathrm{\&theta;}\right)sin\mathrm{\&theta;}d\mathrm{\&theta;}$

${\mathrm{\&delta;}}_k=\left\{\begin{array}{cc}2,& k=0\\ 1,& k\&NotEqual;0\end{array}\right.$

$D=K(\frac{cos\frac{{\mathrm{\&theta;}}_0}{2}}{2}-\frac{4}{3})\frac{{\mathrm{\&mu;sin\&theta;}}_0}{{r_0}^4}$

$f_{\mathrm{\&delta;}r}\left(n\right)=\underset{k=0}{\overset{n}{\&Sigma;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}{\&lsqb;{\stackrel{\&OverBar;}{P}}_{nk}\left(cos\mathrm{\&theta;}\right)\&rsqb;}^2{\cos }^2(\mathrm{\&theta;}+\mathrm{\&xi;}){\sin }^2\mathrm{\&theta;}d\mathrm{\&theta;}=\frac{3}{8}(2n+1)\mathrm{\&pi;}$

$f_{\mathrm{\&delta;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}\left(n\right)=(n+\frac{1}{2})\mathrm{\&pi;}$

r₀＝a_e+h

δσ_n ²是反演重力场模型位系数的阶误差方差，r₀是卫星的地心距，θ₀是两星地心矢量的夹角，a_e是地球平均半径，μ是万有引力常数和地球质量的乘积，h是卫星的轨道高度，T_arc是积分弧长，ΔF是非引力干扰，(Δt)_ΔF是非引力干扰数据间隔，(Δr)_m是卫星定轨的位置误差，(Δt)_Δr是卫星轨道数据的采样间隔，是星间距离变化率测量误差，是星间距离变化率数据采样间隔，T是重力场测量的任务寿命；是引力位系数余弦项的误差，是引力位系数正弦项的误差，是完全规格化的缔合勒让德多项式，K是系数，ξ是相位。

2.根据权利要求1所述的低低星星跟踪卫星重力场测量性能解析计算方法，其特征在于，系数K和相位ξ的取值分别为：

$K=1.3476\&times;{10}^{11}{m}^2,\mathrm{\&xi;}=-\frac{\mathrm{\&pi;}}{2}.$

3.根据权利要求1所述的一种低低星星跟踪卫星重力场测量性能解析计算方法，其特征在于，对于非引力干扰数据间隔(Δt)_ΔF，如果低低星星跟踪重力卫星系统对非引力干扰进行测量，那么(Δt)_ΔF指非引力干扰的测量间隔；如果低低星星跟踪重力场测量系统对非引力干扰进行抑制，那么(Δt)_ΔF指非引力干扰的抑制间隔。

4.根据权利要求1所述的一种低低星星跟踪卫星重力场测量性能解析计算方法，其特征在于，步骤3具体包括建立的关于重力场测量最高有效阶数N_max满足解析关系式：

${{\mathrm{\&delta;\&sigma;}}_{N_{\max }}}^2=\frac{1}{2N_{max}+1}\frac{1.6\&times;{10}^{-10}}{{N_{\max }}^3}.$

5.根据权利要求1所述的一种低低星星跟踪卫星重力场测量性能解析计算方法，其特征在于，步骤4具体包括：

建立的关于n阶对应的大地水准面阶误差为：

${\mathrm{\&Delta;}}_n=a_e\sqrt{(2n+1){{\mathrm{\&delta;\&sigma;}}_n}^2}$

和/或

建立的关于n阶对应的大地水准面累积误差为：

$\mathrm{\&Delta;}=\sqrt{\underset{k=2}{\overset{n}{\&Sigma;}}{\left({\mathrm{\&Delta;}}_k\right)}^2}$

和/或

建立的关于n阶对应的重力异常阶误差为：

${\mathrm{\&Delta;g}}_n=\frac{\mathrm{\&mu;}}{{a_e}^2}(n-1)\sqrt{(2n+1){{\mathrm{\&delta;\&sigma;}}_n}^2}$

和/或

建立的关于n阶对应的重力异常累积误差为：

$\mathrm{\&Delta;}g=\sqrt{\underset{k=2}{\overset{n}{\&Sigma;}}{\left({\mathrm{\&Delta;g}}_k\right)}^2}.$