基于残余星间速度原理反演地球重力场的方法

摘要

本发明涉及一种基于残余星间速度原理反演地球重力场的方法；通过将星载激光干涉测距仪的高精度残余星间速度观测量引入GPS接收机的残余轨道速度差分矢量的视线分量建立新型残余星间速度观测方程，进而精确和快速反演地球重力场的方法；该方法地球重力场计算精度高，卫星重力反演速度快，计算机性能要求低，敏感于重力场中高频信号，易于卫星重力反演误差分析；残余星间速度法是建立高精度和高空间分辨率全球重力场模型的关键技术。

主权项

1.一种基于残余星间速度原理反演地球重力场的方法，其特征如下：步骤一：卫星观测数据采集（1）通过星载激光干涉测距仪获取星间速度通过星载GPS接收机获取双星轨道位置（r₁,r₂）和双星轨道速度通过星载加速度计获取作用于双星的非保守力（f₁，f₂）；（2）利用9阶Runge-Kutta线性单步法和12阶Adams-Cowell线性多步法数值模拟公式获取双星参考轨道位置和双星参考轨道速度（3）参考星间速度通过参考轨道速度计算获得其中，表示相对参考轨道速度矢量，表示第一颗参考卫星指向第二颗参考卫星的参考单位矢量，表示相对参考轨道位置矢量，和表示双星的参考轨道位置矢量；（4）参考非保守力通过DTM2000阻力温度模型计算获得；（5）通过国际公布模型DE-405、IERS96和CSR4.0联合计算获取作用于双星的保守力（F₁,F₂）和参考保守力步骤二：残余星间速度观测方程建立在地心惯性系中，基于牛顿插值原理，单星轨道速度的泰勒展开表示如下$\stackrel{·}{r}\left(t\right)=\stackrel{·}{r}\left(t_0\right)+\underset{i=1}{\overset{n}{\Sigma }}\left(\begin{array}{c}\beta \\ i\end{array}\right)\underset{\alpha =0}{\overset{i}{\Sigma }}{(-1)}^{i+\alpha }\left(\begin{array}{c}i\\ \alpha \end{array}\right)\stackrel{·}{r}\left(t_{\alpha }\right)---\left(1\right)$其中，$\left(\begin{array}{c}\beta \\ i\end{array}\right)$表示二项式系数，t表示计算点的时刻，t₀表示插值点的初始时刻，Δt表示采样间隔，n表示插值点的数量；单星参考轨道速度的泰勒展开表示如下基于公式（1）-（2），单星残余轨道速度的泰勒展开表示如下$\delta \stackrel{·}{r}\left(t\right)=\delta \stackrel{·}{r}\left(t_0\right)+\underset{i=1}{\overset{n}{\Sigma }}\left(\begin{array}{c}\beta \\ i\end{array}\right)\underset{\alpha =0}{\overset{i}{\Sigma }}{(-1)}^{i+\alpha }\left(\begin{array}{c}i\\ \alpha \end{array}\right)\delta \stackrel{·}{r}\left(t_{\alpha }\right)---\left(3\right)$其中，基于公式（3）的一阶时间导数，单星残余轨道加速度的泰勒展开表示如下$\delta \stackrel{··}{r}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\beta \\ i\end{array}\right)}^{\prime }\underset{\alpha =0}{\overset{i}{\Sigma }}{(-1)}^{i+\alpha }\left(\begin{array}{c}i\\ \alpha \end{array}\right)\delta \stackrel{·}{r}\left(t_{\alpha }\right)---\left(4\right)$其中，基于公式（4），双星残余轨道加速度差分的泰勒展开表示如下$\delta {\stackrel{··}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\beta \\ i\end{array}\right)}^{\prime }\underset{\alpha =0}{\overset{i}{\Sigma }}{(-1)}^{i+\alpha }\left(\begin{array}{c}i\\ \alpha \end{array}\right)\delta {\stackrel{·}{r}}_{12}\left(t_{\alpha }\right)---\left(5\right)$其中，和分别表示轨道速度差分矢量和轨道加速度差分矢量，和表示双星的轨道速度矢量，和表示双星的参考轨道速度矢量，和表示双星的轨道加速度矢量，和表示双星的参考轨道加速度矢量；双星残余轨道加速度差分的视线分量表示如下$e_{12}\left(t\right)·\delta {\stackrel{··}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\beta \\ i\end{array}\right)}^{\prime }\underset{\alpha =0}{\overset{i}{\Sigma }}{(-1)}^{i+\alpha }\left(\begin{array}{c}i\\ \alpha \end{array}\right)e_{12}\left(t\right)·\delta {\stackrel{·}{r}}_{12}\left(t_{\alpha }\right)---\left(6\right)$其中，e₁₂=r₁₂/|r₁₂|表示第一颗卫星指向第二颗卫星的单位矢量，r₁₂=r₂-r₁表示双星的轨道位置差分矢量，r₁和r₂分别表示双星的轨道位置矢量；引入激光干涉测距系统的高精度残余星间速度进一步提高地球重力场精度，其中和分别表示星间速度和参考星间速度；可被改写为$\delta {\stackrel{·}{r}}_{12}={\delta \stackrel{·}{r}}_{12}^{\left|\right|}+{\delta \stackrel{·}{r}}_{12}^{\perp }---\left(7\right)$其中，${\delta \stackrel{·}{r}}_{12}^{\left|\right|}=(\delta {\stackrel{·}{r}}_{12}·e_{12})e_{12}$和${\delta \stackrel{·}{r}}_{12}^{\perp }=\delta {\stackrel{·}{r}}_{12}-(\delta {\stackrel{·}{r}}_{12}·e_{12})e_{12}$分别表示的视线分量和垂向分量；基于误差传播原理，为了有效降低将公式（7）中的替换为因此，公式（6）可改写为$e_{12}\left(t\right)·\delta {\stackrel{··}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\beta \\ i\end{array}\right)}^{\prime }\underset{\alpha =0}{\overset{i}{\Sigma }}{(-1)}^{i+\alpha }\left(\begin{array}{c}i\\ \alpha \end{array}\right)e_{12}\left(t\right)·\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{\alpha }\right)---\left(8\right)$其中，$\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{\alpha }\right)=\delta {\stackrel{·}{\rho }}_{12}\left(t_{\alpha }\right)e_{12}\left(t_{\alpha }\right)+\{\delta {\stackrel{·}{r}}_{12}\left(t_{\alpha }\right)-[\delta {\stackrel{·}{r}}_{12}\left(t_{\alpha }\right)·e_{12}\left(t_{\alpha }\right)\left]e_{12}\left(t_{\alpha }\right)\right\};$在公式（8）中，取插值点数n=2,4,6,8，得到2点、4点、6点和8点残余星间速度公式$e_{12}\left(t_i\right)·\delta {\stackrel{··}{r}}_{12}\left(t_i\right)=-\frac{e_{12}\left(t_i\right)}{2\mathrm{\Delta t}}·[\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i-1}\right)-\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i+1}\right)]---\left(9\right)$$e_{12}\left(t_i\right)·\delta {\stackrel{··}{r}}_{12}\left(t_i\right)=\frac{e_{12}\left(t_i\right)}{12\mathrm{\Delta t}}·[\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i-2}\right)-8\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i-1}\right)+8\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i+1}\right)-\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i+2}\right)]---\left(10\right)$$e_{12}\left(t_i\right)·\delta {\stackrel{··}{r}}_{12}\left(t_i\right)=-\frac{e_{12}\left(t_i\right)}{60\mathrm{\Delta t}}·[\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i-3}\right)-9\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i-2}\right)+45\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i-1}\right)---\left(11\right)$$-45\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i+1}\right)+9\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i+2}\right)-\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i+3}\right)]$$e_{12}\left(t_i\right)·\delta {\stackrel{··}{r}}_{12}\left(t_i\right)=\frac{e_{12}\left(t_i\right)}{\mathrm{\Delta t}}·[\frac{1}{280}\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i-4}\right)-\frac{4}{105}\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i-3}\right)+\frac{1}{5}\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i-2}\right)-\frac{4}{5}\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i-1}\right)---\left(12\right)$$+\frac{4}{5}\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i+1}\right)-\frac{1}{5}\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i+2}\right)+\frac{4}{105}\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i+3}\right)-\frac{1}{280}\delta {\stackrel{·}{r}}_{\rho 12}\left(t_{i+4}\right)]$在公式（8）中，的具体形式表示如下$\delta {\stackrel{··}{r}}_{12}=\delta g_{12}+\delta T_{12}+\delta F_{12}+\delta f_{12}---\left(13\right)$其中，δT₁₂表示作用于双星的残余地球扰动引力差；表示除地球引力之外的残余保守力差，F₁和F₂表示作用于双星的保守力，和表示参考保守力；表示残余非保守力差，f₁和f₂表示作用于双星的非保守力，和表示参考非保守力；表示残余地心引力差，g₁和g₂表示双星的地心引力，和表示参考地心引力其中，GM表示地球质量M和万有引力常数G的乘积，分别表示双星的地心半径，x₁₍₂₎,y₁₍₂₎,z₁₍₂₎表示轨道位置矢量r₁₍₂₎的3个分量；通过将公式（13）和（14）代入（8），残余星间速度观测方程表示如下$e_{12}\left(t\right)·\delta T_{12}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\beta \\ i\end{array}\right)}^{\prime }\underset{\alpha =0}{\overset{i}{\Sigma }}{(-1)}^{i+\alpha }\left(\begin{array}{c}i\\ \alpha \end{array}\right)e_{12}\left(t\right)·\{\delta {\stackrel{·}{\rho }}_{12}\left(t_{\alpha }\right)e_{12}\left(t_{\alpha }\right)$$+[\delta {\stackrel{·}{r}}_{12}\left(t_{\alpha }\right)-(\delta {\stackrel{·}{r}}_{12}\left(t_{\alpha }\right)·e_{12}\left(t_{\alpha }\right))e_{12}\left(t_{\alpha }\right)]\}---\left(15\right)$$-e_{12}\left(t\right)·[\delta g_{12}+\delta F_{12}\left(t\right)+\delta f_{12}\left(t\right)]$其中，表示残余地球扰动位的一阶梯度，V₁和V₂表示地球扰动位，和表示参考地球扰动位$V(r,\theta ,\lambda )=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\Sigma }}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\Sigma }}({\bar{C}}_{\mathrm{lm}}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{lm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}\left(\cos \theta \right)---\left(16\right)$其中，r,θ，λ分别表示地心半径、地心余纬度和地心经度，R_e表示地球平均半径；表示正规化的缔合Legendre函数，l表示阶数，m表示次数；和表示待估的地球引力位系数；步骤三：地球重力场反演基于残余星间速度观测方程（15），利用激光干涉测距仪的星间速度、GPS接收机的轨道位置和轨道速度、加速度计的非保守力、以及双星的参考非保守力、保守力和参考保守力，求解地球引力位系数和最终通过地球引力位系数的集合建立全球重力场模型。