利用星间距离插值建立全球重力场模型的方法

摘要

本发明公开了一种利用星间距离插值建立全球重力场模型的方法，通过将星载K波段测量仪的精确星间距离引入双星相对轨道位置矢量的星星连线分量，构建星间距离插值卫星观测方程，进而建立高精度和高空间分辨率的全球重力场模型。该方法卫星重力反演精度高，易于感测中高频重力场信号，利于重力卫星系统误差分析，观测方程物理含义明确，计算机性能要求低。

主权项

1.一种利用星间距离插值建立全球重力场模型的方法，包含下列步骤：步骤一：对GRACE卫星观测数据进行预处理，具体包括1.1)采集星载K波段测量仪得到的星间距离ρ₁₂数据：基于t检验准则即罗曼诺夫斯基准则，剔除星间距离数据中存在的粗大误差；基于9阶Lagrange多项式，插值获得间断的星间距离数据；1.2)采集星载双频GPS接收机得到的卫星轨道数据，包括轨道位置r和轨道速度：为了保证卫星轨道数据的精度和连续性，去除卫星轨道存在的重叠期，进行卫星轨道数据的拼接；截掉由于定轨弱约束造成的卫星轨道数据的开始和结束时段处精度较低的数据；基于3σ准则即莱以特准则，剔除卫星轨道数据中存在的粗大误差；1.3)采集星载加速度计得到的卫星非保守力f数据：基于t检验准则即罗曼诺夫斯基准则，剔除卫星非保守力数据中存在的粗大误差；基于9阶Lagrange多项式，插值获得间断的卫星非保守力数据；步骤二：构建星间距离插值观测方程在地心惯性坐标系中，基于Newton插值模型，单星轨道位置r的泰勒展开表示如下$r\left(t\right)=r\left(t_0\right)+\underset{j=1}{\overset{n}{\Sigma }}\left(\begin{array}{c}\alpha \\ j\end{array}\right)\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)r\left(t_{\xi }\right),---\left(1\right)$其中，表示二项式系数，t表示插值点的时间，t₀表示插值点的初始时刻，Δt表示采样间隔，n表示插值点的个数；在(1)式两边同时对时间t求二阶导数，可得单星轨道加速度的展开公式$\stackrel{··}{r}\left(t\right)=\underset{j=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)r\left(t_{\xi }\right),---\left(2\right)$基于(2)式，双星轨道加速度差分的展开公式表示如下${\stackrel{··}{r}}_{12}\left(t\right)=\underset{j=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)r_{12}\left(t_{\xi }\right),---\left(3\right)$其中，r₁₂＝r₂-r₁和分别表示双星相对轨道位置矢量和相对轨道加速度矢量，r₁和r₂分别表示双星绝对轨道位置矢量，和分别表示双星绝对轨道加速度矢量；将(3)式中的投影到星星连线方向LOS可得$e_{12}\left(t\right)·{\stackrel{··}{r}}_{12}\left(t\right)=\underset{j=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)e_{12}\left(t\right)·r_{12}\left(t_{\xi }\right),---\left(4\right)$其中，e₁₂＝r₁₂/|r₁₂|表示由GRACE-A卫星指向GRACE-B卫星的单位矢量；e₁₂(t)·r₁₂(t_ξ)可改写为$e_{12}\left(t\right)·r_{12}\left(t_{\xi }\right)=e_{12}\left(t\right)·[r_{12}^{\left|\right|}\left(t_{\xi }\right)+r_{12}^{\perp }\left(t_{\xi }\right)],---\left(5\right)$其中，表示r₁₂的星星连线方向分量；表示r₁₂的垂直于星星连线方向分量；通过将GRACE卫星K波段测量仪高精度的星间距离ρ₁₂e₁₂替换(r₁₂·e₁₂)e₁₂，(5)式可改写为e₁₂(t)·r₁₂(t_ξ)＝e₁₂(t)·[ρ₁₂(t_ξ)e₁₂(t_ξ)+{r₁₂(t_ξ)-[r₁₂(t_ξ)·e₁₂(t_ξ)]e₁₂(t_ξ)}]，     (6)将(6)式代入(4)式可得$e_{12}\left(t\right)·{\stackrel{··}{r}}_{12}\left(t\right)=\underset{j=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)e_{12}\left(t\right)·r_{\rho 12}\left(t_{\xi }\right),---\left(7\right)$其中，r_ρ12(t_ξ)＝ρ₁₂(t_ξ)e₁₂(t_ξ)+{r₁₂(t_ξ)-[r₁₂(t_ξ)·e₁₂(t_ξ)]e₁₂(t_ξ)}；在(7)式中，的具体形式表示如下${\stackrel{··}{r}}_{12}=g_{12}^0+g_{12}^T+a_{12}+f_{12},---\left(8\right)$其中，表示作用于双星的相对地球扰动引力，表示梯度算子；a₁₂表示除地球引力之外的作用于双星的相对保守力；f₁₂表示作用于双星的相对非保守力；表示作用于双星的相对地球中心引力$g_{12}^0=-\mathrm{GM}(\frac{r_{12}}{{\left|r_2\right|}^3}-\frac{r_1}{{\left|r_1\right|}^3}),---\left(9\right)$其中，GM表示地球质量M和万有引力常数G之积，表示双星各自的地心半径，x₁₍₂₎，y₁₍₂₎，z₁₍₂₎分别表示双星各自位置矢量r₁₍₂₎的三个分量；将(8)式和(9)式代入(7)式，星间距离插值观测方程表示如下$e_{12}\left(t\right)·▿T_{12}\left(t\right)=e_{12}\left(t\right)·\{\underset{j=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)r_{\rho 12}\left(t_{\xi }\right),$(10)$+\mathrm{GM}[\frac{r_2\left(t\right)}{{\left|r_2\left(t\right)\right|}^3}-\frac{r_1\left(t\right)}{{\left|r_1\left(t\right)\right|}^3}]-a_{12}\left(t\right)-f_{12}\left(t\right)\}$其中，表示相对扰动位梯度，T(r，θ，λ)表示地球扰动位$T(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(11\right)$其中，r，θ和λ分别表示卫星的地心半径、余纬度和经度，R_e表示地球的平均半径，L表示球函数展开的最大阶数；表示规格化的Legendre函数，l表示阶数，m表示次数；和表示待求的规格化引力位系数；将(10)式按照泰勒公式展开，3点、5点、7点和9点星间距离插值公式表示如下$e_{12}\left(t_j\right)·▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)·\left\{\frac{1}{{\left(\mathrm{\Delta t}\right)}^2}\right[r_{\rho 12}\left(t_{j-1}\right)-{2r}_{\rho 12}\left(t_j\right)+r_{\rho 12}\left(t_{j+1}\right)]$(12)$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\},$$e_{12}\left(t_j\right)·▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)·\left\{\frac{1}{{\left(\mathrm{\Delta t}\right)}^2}\right[-\frac{1}{12}{r}_{\rho 12}\left(t_{j-2}\right)+\frac{4}{3}{r}_{\rho 12}\left(t_{j-1}\right)$$-\frac{5}{2}{r}_{\rho 12}\left(t_j\right)+\frac{4}{3}{r}_{\rho 12}\left(t_{j+1}\right)-\frac{1}{12}{r}_{\rho 12}\left(t_{j+2}\right)],---\left(13\right)$$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\}$$e_{12}\left(t_j\right)·▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)·\left\{\frac{1}{{\left(\mathrm{\Delta t}\right)}^2}\right[\frac{1}{90}{r}_{\rho 12}\left(t_{j-3}\right)-\frac{3}{20}{r}_{\rho 12}\left(t_{j-2}\right)+\frac{3}{2}{r}_{\rho 12}\left(t_{j-1}\right)$$-\frac{49}{18}{r}_{\rho 12}\left(t_j\right)+\frac{3}{2}{r}_{\rho 12}\left(t_{j+1}\right)-\frac{3}{20}{r}_{\rho 12}\left(t_{j+2}\right)+\frac{1}{90}{r}_{\rho 12}\left(t_{j+3}\right)],---\left(14\right)$$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\}$$e_{12}\left(t_j\right)·▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)·\left\{\frac{1}{{\left(\mathrm{\Delta t}\right)}^2}\right[-\frac{1}{560}{r}_{\rho 12}\left(t_{j-4}\right)+\frac{8}{315}{r}_{\rho 12}\left(t_{j-3}\right)-\frac{1}{5}{r}_{\rho 12}\left(t_{j-2}\right)+\frac{8}{5}{r}_{\rho 12}\left(t_{j-1}\right)$$-\frac{205}{72}{r}_{\rho 12}\left(t_j\right)+\frac{8}{5}{r}_{\rho 12}\left(t_{j+1}\right)-\frac{1}{5}{r}_{\rho 12}\left(t_{j+2}\right)+\frac{8}{315}{r}_{\rho 12}\left(t_{j+3}\right)-\frac{1}{560}{r}_{\rho 12}\left(t_{j+4}\right)],---\left(15\right)$$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\}$步骤三：优选不同点数星间距离插值公式利用“步骤一”中获取的GRACE星载K波段测量仪的星间距离、星载GPS接收机的卫星轨道位置和卫星轨道速度、以及星载加速度计的卫星非保守力观测数据，基于“步骤二”中3点、5点、7点和9点星间距离插值公式(12)～(15)分别反演地球重力场，可得到9点星间距离插值公式(15)有利于120阶GRACE卫星重力反演精度的有效提高；步骤四：建立全球重力场模型利用“步骤一”中获取的GRACE星载K波段测量仪的星间距离、星载GPS接收机的卫星轨道位置和卫星轨道速度、以及星载加速度计的卫星非保守力观测数据，基于“步骤二”中建立的9点星间距离插值公式(15)，计算获得地球引力位系数和最终通过引力位系数的集合建立120阶全球重力场模型WHIGG-GEGM01S；在120阶处，大地水准面累积误差为10.98cm，重力异常累积误差为1.741×10^-6m/s²；步骤五：验证全球重力场模型的正确性和可靠性全球重力场模型的大地水准面标准误差为0.726m，更接近于国际已公布的全球重力场模型EIGEN-GRACE02S的大地水准面标准误差0.735m。