基于双星能量插值原理的卫星重力反演方法

摘要

本发明涉及卫星大地测量学、地球物理学、空间科学等交叉技术领域，提出一种基于双星能量插值原理的卫星重力反演方法。通过将星载K波段测量仪的高精度星间距离插值引入双星相对轨道动能项，进而建立双星能量插值观测方程，旨在精确和快速反演全球重力场。该方法卫星重力反演精度较高，观测方程物理含义明确，利于重力卫星系统误差分析，易于感测中高频重力场信号，计算机性能要求较低。由于可实质性提高120阶GRACE地球重力场的感测精度，因此本发明提出的双星能量插值法是建立高精度和高空间分辨率的新一代全球重力场模型的有效方法。

主权项

1.一种基于双星能量插值原理的卫星重力反演方法，包含下列步骤：步骤一：对卫星观测数据进行预处理，具体包括：1.1)采集星载K波段测量仪得到的星间距离数据ρ₁₂：基于t检验准则即罗曼诺夫斯基准则，剔除星间距离数据中存在的粗大误差；基于9阶Lagrange多项式，插值获得间断的星间距离数据；1.2)采集星载双频GPS接收机得到的卫星轨道数据，包括轨道位置r和轨道速度去除卫星轨道存在的重叠期，进行轨道数据的拼接；截掉由于定轨弱约束造成的卫星轨道数据的开始和结束时段处精度较低的数据；基于3σ准则即莱以特准则，剔除轨道数据中存在的粗大误差；1.3)采集星载加速度计得到的卫星非保守力数据f：基于t检验准则即罗曼诺夫斯基准则，剔除卫星非保守力数据中存在的粗大误差；基于9阶Lagrange多项式，插值获得间断的卫星非保守力数据；步骤二：基于双星能量插值原理反演全球重力场2.1)在地心惯性系中，建立卫星运动方程如下，$\stackrel{··}{r}(r,t^{\prime })=F_e(r,t^{\prime })+F_T(r,t^{\prime })+f(r,t^{\prime }),---\left(1\right)$其中，表示卫星总加速度，r表示卫星的绝对位置矢量，t′表示观测时间；F_e(r，t′)表示地球的引力；F_T(r，t′)表示三体摄动力；f(r，t′)表示作用于卫星的非保守力；在式(1)两边同时乘以轨道速度$\stackrel{·}{r}·\stackrel{··}{r}(r,t^{\prime })=\stackrel{·}{r}·[F_e(r,t^{\prime })+F_T(r,t^{\prime })]+\stackrel{·}{r}·f(r,t^{\prime }),---\left(2\right)$其中，F_e(r，t′)和F_T(r，t′)可表示为$F_{e\left(T\right)}(r,t^{\prime })=\frac{{\partial V}_{e\left(T\right)}}{\partial r},---\left(3\right)$其中，V_e＝V₀+T_e表示地球引力位，表示地球中心引力位，GM表示地球质量M和万有引力常数G之积，表示卫星的地心半径，x，y，z表示位置矢量r的分量，T_e表示扰动位；V_T表示三体摄动能；V_e(T)对时间的一阶导数表示如下：$\frac{{\mathrm{dV}}_{e\left(T\right)}}{{\mathrm{dt}}^{\prime }}=\frac{{\partial V}_{e\left(T\right)}}{\partial r}·\frac{\mathrm{dr}}{{\mathrm{dt}}^{\prime }}+\frac{{\partial V}_{e\left(T\right)}}{{\partial t}^{\prime }}·\frac{{\mathrm{dt}}^{\prime }}{{\mathrm{dt}}^{\prime }},---\left(4\right)$将式(3)代入式(4)中可得$F_{e\left(T\right)}(r,t^{\prime })·\stackrel{·}{r}=\frac{{\mathrm{dV}}_{e\left(T\right)}}{{\mathrm{dt}}^{\prime }}-\frac{{\partial V}_{e\left(T\right)}}{{\partial t}^{\prime }},---\left(5\right)$将式(5)代入式(2)，并两边同时积分；在地心惯性系中，单星扰动位观测方程如下：T_e＝E_k-E_f+V_ω-V_T-V₀-E₀，        (6)其中，表示卫星的动能；表示卫星的耗散能；$V_{\omega }=\int \frac{\partial (V_e+V_T)}{{\partial t}^{\prime }}{\mathrm{dt}}^{\prime }\approx -{\omega }_e(x\stackrel{·}{y}-y\stackrel{·}{x})$表示卫星的位旋转能，ω_e表示地球自转角速度，表示卫星轨道速度的3个分量；E₀表示卫星系统的能量积分常数；据式(6)，在地心惯性系中，双星扰动位差的观测方程建立如下：T_e12＝E_k12-E_f12+V_ω12-V_T12-V₀₁₂-E₀₁₂，        (7)其中，T_e12表示双星扰动位差$T_{e12}(r_1,{\theta }_1,{\lambda }_1,r_2,{\theta }_2,{\lambda }_2)=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\Sigma }}\underset{m=-l}{\overset{l}{\Sigma }}\left\{\right[{\left(\frac{R_e}{r_2}\right)}^{l+1}{\bar{Y}}_{\mathrm{lm}}({\theta }_2,{\lambda }_2)-{\left(\frac{R_e}{r_1}\right)}^{l+1}{\bar{Y}}_{\mathrm{lm}}({\theta }_1,{\lambda }_1)\left]{\bar{C}}_{\mathrm{lm}}\right\}$其中，${\bar{Y}}_{l,m}(\theta ,\lambda )={\bar{P}}_{l,\left|m\right|}\left(\cos \theta \right)Q_m\left(\lambda \right),$$Q_m\left(\lambda \right)=\left\{\begin{array}{cc}\cos \mathrm{m\lambda }& m\ge 0\\ \sin \left|m\right|\lambda & m<0\end{array}\right.;$r₁₍₂₎，θ₁₍₂₎，λ₁₍₂₎表示双星各自的地心半径、地心余纬度和地心经度，R_e表示地球的平均半径；表示规格化的Legendre函数，其中l表示阶数，m表示次数；表示待求的规格化地球引力位系数；表示双星动能差，$E_{f12}=\int ({\stackrel{·}{r}}_2·f_2-{\stackrel{·}{r}}_1·f_1){\mathrm{dt}}^{\prime }$表示双星耗散能差，$V_{\omega 12}=-{\omega }_e(x_{12}{\stackrel{·}{y}}_2-y_2{\stackrel{·}{x}}_{12}-y_{12}{\stackrel{·}{x}}_1+x_1{\stackrel{·}{y}}_{12})$表示双星位旋转能差，V_T12表示双星三体摄动能差，表示双星中心引力位差，E₀₁₂表示双星能量积分常数差，可通过初始位置和速度计算得到；2.2)将GRACE卫星的K波段测量仪的精确星间距离数据ρ₁₂引入到式(7)中的项E_k12，建立双星能量插值卫星观测方程；包括：在地心惯性坐标系中，基于Newton插值模型，单星轨道位置r的泰勒展开表示式，如下$r\left(t\right)=r\left(t_0\right)+\underset{i=1}{\overset{n}{\Sigma }}\left(\begin{array}{c}\lambda \\ i\end{array}\right)\underset{\tau =0}{\overset{i}{\Sigma }}{(-1)}^{i+\tau }\left(\begin{array}{c}i\\ \tau \end{array}\right)r\left(t_{\tau }\right),---\left(8\right)$其中，表示二项式系数，t表示插值点的时间，t₀表示插值点的初始时刻，Δt表示采样间隔，n表示插值点的个数；在式(8)两边同时对t求一阶导数，可得单星轨道速度的展开公式：$\stackrel{·}{r}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\lambda \\ i\end{array}\right)}^{\prime }\underset{\tau =0}{\overset{i}{\Sigma }}{(-1)}^{i+\tau }\left(\begin{array}{c}i\\ \tau \end{array}\right)r\left(t_{\tau }\right).---\left(9\right)$基于式(9)，双星轨道速度差的展开公式表示如下${\stackrel{·}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\lambda \\ i\end{array}\right)}^{\prime }\underset{\tau =0}{\overset{i}{\Sigma }}{(-1)}^{i+\tau }\left(\begin{array}{c}i\\ \tau \end{array}\right)r_{12}\left(t_{\tau }\right),---\left(10\right)$其中，r₁₂＝r₂-r₁和分别表示双星相对轨道位置矢量和相对轨道速度矢量，r₁和r₂分别表示双星绝对轨道位置矢量，和分别表示双星绝对轨道速度矢量；基于式(10)，双星相对动能的展开公式表示如下：$\frac{1}{2}{[{\stackrel{·}{r}}_2\left(t\right)+{\stackrel{·}{r}}_1\left(t\right)]·{\stackrel{·}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}\left(\begin{array}{c}\lambda \\ i\end{array}\right)}^{\prime }\underset{\tau =0}{\overset{i}{\Sigma }}{(-1)}^{i+\tau }\left(\begin{array}{c}i\\ \tau \end{array}\right)\frac{1}{2}[{\stackrel{·}{r}}_2\left(t\right)+{\stackrel{·}{r}}_1\left(t\right)]{·r}_{12}\left(t_{\tau }\right),---\left(11\right)$其中，$\frac{1}{2}[{\stackrel{·}{r}}_2\left(t\right)+{\stackrel{·}{r}}_1\left(t\right)]·r_{12}\left(t_{\tau }\right)$可被改写为：$\frac{1}{2}[{\stackrel{·}{r}}_2\left(t\right)+{\stackrel{·}{r}}_1\left(t\right)]·r_{12}\left(t_{\tau }\right)=\frac{1}{2}[{\stackrel{·}{r}}_2\left(t\right)+{\stackrel{·}{r}}_1\left(t\right)]·[r_{12}^{\left|\right|}\left(t_{\tau }\right)+r_{12}^{\perp }\left(t_{\tau }\right)],---\left(12\right)$其中，表示r₁₂的星星连线方向分量，e₁₂＝r₁₂/|r₁₂|表示由GRACE-A卫星指向GRACE-B卫星的单位矢量；表示r₁₂的垂直于星星连线方向分量；使用GRACE卫星K波段测量仪的高精度星间距离ρ₁₂e₁₂来替换(r₁₂·e₁₂)e₁₂；式(11)可改写为：$\frac{1}{2}{[{\stackrel{·}{r}}_2\left(t\right)+{\stackrel{·}{r}}_1\left(t\right)]·{\stackrel{·}{r}}_{\rho 12}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}\left(\begin{array}{c}\lambda \\ i\end{array}\right)}^{\prime }\underset{\tau =0}{\overset{i}{\Sigma }}{(-1)}^{i+\tau }\left(\begin{array}{c}i\\ \tau \end{array}\right)\frac{1}{2}[{\stackrel{·}{r}}_2\left(t\right)+{\stackrel{·}{r}}_1\left(t\right)]{·r}_{\rho 12}\left(t_{\tau }\right),---\left(13\right)$其中，r_ρ12(t_τ)＝ρ₁₂(t_τ)e₁₂(t_τ)+{r₁₂(t_τ)-[r₁₂(t_τ)e₁₂(t_τ)]e₁₂(t_τ)}；因此，2点、4点、6点和8点星间距离插值公式表示如下${\stackrel{·}{r}}_{\rho 12}\left(t_i\right)=-\frac{1}{2\mathrm{\Delta t}}[r_{\rho 12}\left(t_{i-1}\right)-r_{\rho 12}\left(t_{i+1}\right)],---\left(14\right)$${\stackrel{·}{r}}_{\rho 12}\left(t_i\right)=\frac{1}{12\mathrm{\Delta t}}[r_{\rho 12}\left(t_{i-2}\right)-{8r}_{\rho 12}\left(t_{i-1}\right)+{8r}_{\rho 12}\left(t_{i+1}\right)-r_{\rho 12}\left(t_{i+2}\right)],---\left(15\right)$${\stackrel{·}{r}}_{\rho 12}\left(t_i\right)=-\frac{1}{60\mathrm{\Delta t}}[r_{\rho 12}\left(t_{i-3}\right)-{9r}_{\rho 12}\left(t_{i-2}\right)+{45r}_{\rho 12}\left(t_{i-1}\right),---\left(16\right)$$-45r_{\rho 12}\left(t_{i+1}\right)+{9r}_{\rho 12}\left(t_{i+2}\right)-r_{\rho 12}\left(t_{i+3}\right)]$${\stackrel{·}{r}}_{\rho 12}\left(t_i\right)=\frac{1}{\mathrm{\Delta t}}[\frac{1}{280}{r}_{\rho 12}\left(t_{i-4}\right)-\frac{4}{105}{r}_{\rho 12}\left(t_{i-3}\right)+\frac{1}{5}{r}_{\rho 12}\left(t_{i-2}\right)-\frac{4}{5}{r}_{\rho 12}\left(t_{i-1}\right)---\left(17\right)$$+\frac{4}{5}{r}_{\rho 12}\left(t_{i+1}\right)-\frac{1}{5}{r}_{\rho 12}\left(t_{i+2}\right)+\frac{4}{105}{r}_{\rho 12}\left(t_{i+3}\right)-\frac{1}{280}{r}_{\rho 12}\left(t_{i+4}\right)]$联合式(7)和式(13)，双星能量插值观测方程表示如下T_e12(t)＝E_ρ12(t)-E_f12(t)+V_ω12(t)-V_T12(t)-V₀₁₂(t)-E₀₁₂(t)，    (18)其中，$E_{\rho 12}{\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}\left(\begin{array}{c}\lambda \\ i\end{array}\right)}^{\prime }\underset{\tau =0}{\overset{i}{\Sigma }}{(-1)}^{i+\tau }\left(\begin{array}{c}i\\ \tau \end{array}\right)\frac{1}{2}[{\stackrel{·}{r}}_2\left(t\right)+{\stackrel{·}{r}}_1\left(t\right)]{·r}_{\rho 12}\left(t_{\tau }\right)$表示双星动能差，V_T12(t)＝V_E12(t)+V_S12(t)+V_M12(t)，V_E12(t)表示双星地球固体潮汐能差，V_S12(t)表示双星太阳引力位差，V_M12(t)表示双星月球引力位差，E_f12表示双星耗散能差，V_ω12表示双星位旋转能差，V₀₁₂表示双星中心引力位差，E₀₁₂表示双星能量积分常数差；2.3)基于2点、4点、6点和8点双星能量插值观测方程(18)分别反演120阶GRACE地球重力场。