CN102998713B - 基于功率谱半解析的卫星重力梯度反演方法 - Google Patents

Abstract

本发明涉及一种基于功率谱半解析的卫星重力梯度反演方法；建立了卫星重力梯度仪的重力梯度张量误差和GPS/GLONASS复合接收机的轨道位置误差影响累计大地水准面精度的单独和联合解析误差模型。同时还提出了对卫星重力梯度反演半解析误差模型进行验证的方法；以及基于该卫星重力梯度反演半解析误差模型，通过卫星重力梯度测量精度、卫星轨道位置测量精度和卫星轨道高度对累计大地水准面的影响确定下一代GOCE Follow-On重力梯度卫星参数需求的方法。该卫星重力梯度反演方法精度高，地球重力场解算速度快，卫星观测方程物理含义明确，较大程度简化了计算过程，以及易于开展下一代重力梯度卫星系统需求分析。

Description

基于功率谱半解析的卫星重力梯度反演方法

技术领域

本发明涉及卫星重力梯度学、大地测量学、空间科学、宇航学等交叉技术领域，特别是涉及一种基于功率谱原理精确建立卫星重力梯度反演半解析误差模型和开展GOCE Follow-On需求论证研究的方法。

背景技术

卫星重力梯度测量（SGG）技术的实现是继美国全球定位系统（GPS）星座成功构建之后在大地测量领域的又一项创新和突破，被国际大地测量学界公认为是当前地球重力场探测研究中最高效、最经济和最有发展潜力的方法之一。欧空局（ESA）独立研制的GOCE重力梯度卫星已于2009年3月17日成功发射升空。GOCE采用近圆、极地和太阳同步低轨道，轨道倾角96.5°，轨道离心率0.001，轨道高度250km。为了最大程度减少空间环境扰动导致卫星姿态的变化，GOCE设计为严格对称的八角形棱柱体。GOCE卫星采用卫星跟踪卫星高低（SST-HL）和卫星重力梯度的结合模式，除基于高轨道的GPS/GLONASS卫星对低轨道的GOCE卫星进行精密跟踪定位，利用定位于卫星质心处的重力梯度仪高精度测量卫星轨道高度处引力位的二阶导数，同时利用非保守力补偿技术（Drag-free）精密屏蔽作用于卫星体的大气阻力、太阳光压、地球辐射压、以及轨道高度和姿态控制力等。由于地球重力场信号随卫星轨道高度的增加而急剧衰减(R_e/(R_e+H))^l+1，基于卫星跟踪卫星模式（SST）仅适合于精密确定地球中长波重力场，而卫星重力梯度测量是直接测定地球引力位的二次微分，其结果是将球谐系数放大了l²倍，因此卫星重力梯度测量可抑制地球引力位随高度的衰减效应，进而高精度感测地球中短波重力场信号。

目前国内外科研机构采用的卫星重力梯度反演法主要包括：空域法、时域法、时空域混合法、直接法等。（1）空域法是指不直接处理空间位置相对不规则的卫星轨道采样点的观测值，而将卫星观测值归算到以卫星平均轨道高度为半径的球面上利用快速傅立叶变换（FFT）进行网格化处理，将问题转化为某类型边值问题的解。优点：因网格点数固定从而方程维数一定，且可以利用FFT方法进行快速批量处理，因此极大地降低了计算量；缺点：在进行网格化处理中作了不同程度的近似计算，且不能对色噪声进行处理。（2）时域法是指将卫星观测数据按时间序列处理，卫星星历值直接表示成引力位系数的函数，由最小二乘等方法直接反求引力位系数。优点：直接对卫星观测数据进行处理，不需作任何近似，求解精度较高且能有效处理色噪声；缺点：随着卫星观测数据的增多，观测方程数量剧增，极大地增加了计算量。（3）时空域混合法是指联合空域法的快速性和时域法的精确性反演地球重力场。优点：在保证地球重力场解算精度的前提下，有效改善了计算速度；缺点：相对于单独的空域法和时域法，计算过程较复杂。（4）直接法是指将卫星精密定轨和地球重力场反演合二为一，基于各种卫星观测值同时求解卫星轨道、地面站坐标、地球自转参数、海潮模型和地球重力场模型、以及其它动力学和非动力学参数，通过综合卫星运动学、卫星动力学、大地测量学、地球物理学等多学科的知识建立的一种合乎自然规律的解算方法。优点：不依赖于任何先验的地球重力场模型，理论框架严密，各种地球重力场参数求解精度较高；缺点：整体解算过程较复杂，需要高性能的并行计算机支持。

为了满足本世纪科学和国防对地球重力场精度进一步提高的迫切需求，以及由于GOCE重力梯度卫星的工作寿命预期于2014年前结束，因此目前国内外科研机构正积极开展更高精度的GOCE Follow-On卫星重力梯度测量计划的需求分析和载荷研制。由于现有卫星重力梯度反演法的计算过程较复杂和计算速度较慢，因此，本发明首次建立了卫星重力梯度仪的重力梯度张量误差和GPS/GLONASS复合接收机的轨道位置误差影响累计大地水准面精度的单独和联合半解析误差模型，进而精确和快速地开展了GOCE Follow-On重力梯度卫星的需求论证研究。

发明内容

为了解决上述技术问题，本发明提供了一种基于功率谱半解析的卫星重力梯度反演方法，包含以下步骤：

步骤1，通过星载重力梯度仪采集卫星重力梯度测量数据δV_xyz，基于星载GPS/GLONASS复合接收机采集卫星轨道位置测量数据δr；

步骤2，建立卫星重力梯度的信号功率谱分析模型，并对所采集的卫星重力梯度测量数据进行信号功率谱敏感度分析；

步骤3，通过卫星重力梯度张量误差的功率谱分析和卫星轨道位置误差的功率谱分析建立卫星重力梯度反演半解析误差模型；其中，所述步骤3包括：

步骤3.1，通过卫星重力梯度仪的重力梯度张量误差对累计大地水准面精度的影响建立卫星重力梯度张量的半解析误差模型；

步骤3.2，通过GPS/GLONASS复合接收机的轨道位置误差对累计大地水准面精度的影响建立卫星轨道位置的半解析误差模型；

步骤3.3，通过卫星重力梯度张量的半解析误差模型和卫星轨道位置的半解析误差模型，建立卫星重力梯度和轨道位置的联合半解析误差模型，以此作为卫星重力梯度反演半解析误差模型；

步骤4，使用所述卫星重力梯度反演半解析误差模型，以及采集得到的卫星重力梯度测量数据δV_xyz和卫星轨道位置测量数据δr反演累计大地水准面误差。

本发明还提供了一种对基于功率谱半解析的卫星重力梯度反演方法中所使用的卫星重力梯度反演半解析误差模型进行验证的方法，包括如下步骤：

步骤5，使用所述卫星重力梯度反演半解析误差模型确定300阶GOCE地球重力场精度，并通过与国际公布的GO_CONS_GCF_2_TIM_R2地球重力场模型精度对比，验证所述卫星重力梯度反演半解析误差模型的准确性。

本发明还提供了一种使用基于功率谱半解析的卫星重力梯度反演方法来确定GOCE Follow-On重力梯度卫星参数需求的方法，包括如下步骤：

步骤6，根据地球科学各相关学科需求确定所需的地球重力场反演精度；通过所述卫星重力梯度反演半解析误差模型，依据地球重力场需求精度对GOCEFollow-On重力梯度卫星的关键载荷精度指标和轨道参数进行需求分析。

本发明取得了以下技术效果：

建立卫星重力梯度仪的重力梯度张量误差和GPS/GLONASS复合接收机的轨道位置误差影响累计大地水准面精度的单独和联合半解析误差模型，进而精确和快速反演GOCE Follow-On地球重力场。其优点是：

1）卫星重力梯度反演精度高；

2）地球重力场解算速度快；

3）卫星观测方程物理含义明确；

4）较大程度简化了计算过程；

5）易于开展重力梯度卫星系统需求分析。

附图说明

图1表示卫星重力梯度张量功率谱的敏感度系数表。

图2表示卫星重力梯度张量的敏感度系数|A_ab|（每阶）。

图3表示卫星重力梯度张量的信号幅度谱（每阶）。

图4表示基于GOCE关键载荷匹配精度指标分别估计累计大地水准面精度。

图5表示基于不同卫星重力梯度测量精度估计累计大地水准面精度。

图6表示基于不同卫星轨道位置测量精度估计累计大地水准面精度。

图7表示基于不同卫星轨道高度估计累计大地水准面精度。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及具体实施方式对本发明作进一步的详细描述。

基于功率谱半解析的卫星重力梯度反演方法包含下列步骤：

步骤一：重力梯度卫星数据采集

（1）通过重力梯度卫星的星载重力梯度仪采集卫星重力梯度测量数据δV_xyz；

（2）基于重力梯度卫星的星载GPS/GLONASS复合接收机采集卫星轨道位置测量数据δr。

步骤二：建立卫星重力梯度张量的信号功率谱分析模型

在球坐标系中，地球引力位按球谐函数展开的表达式为

$V(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e}\underset{l=0}{\overset{L}{\mathrm{\Σ}}}{\left(\frac{R_e}{r}\right)}^{1+1}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{lm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right),---\left(1\right)$

其中，GM表示地球质量M和万有引力常数G之积，R_e表示地球的平均半径，L表示地球引力位按球函数展开的最大阶数；表示卫星的地心半径，x,y,z分别表示卫星位置矢量r在直角坐标系中的三个分量，θ和λ表示地心余纬度和经度；表示规格化的Legendre函数，l表示阶数，m表示次数；和表示待求的规格化引力位系数。

球坐标系(r,θ,λ)和直角坐标系(x,y,z)的互换公式表示为，

$\left\{\begin{array}{c}r=\sqrt{x^2+y^2+z^2}\\ \sin \mathrm{\θ}=\frac{\sqrt{x^2+y^2}}{\sqrt{x^2+y^2+z^2}}\\ \cos \mathrm{\θ}=\sqrt{\frac{z}{\sqrt{x^2+y^2+z^2}}}---\left(2\right)\\ \sin \mathrm{\λ}=\sqrt{\frac{y}{x^2+y^2}}\\ \cos \mathrm{\λ}=\frac{x}{\sqrt{x^2+y^2}}\end{array}\right.$

基于公式（1）和（2），在直角坐标系中，地球引力位V分别对x,y,z的二阶导数表示如下

$\mathrm{\Γ}=\left[\begin{array}{ccc}{V}_{\mathrm{xx}}& V_{\mathrm{xy}}& V_{\mathrm{xz}}\\ V_{\mathrm{yx}}& V_{\mathrm{yy}}& V_{\mathrm{yz}}\\ V_{\mathrm{zx}}& V_{\mathrm{zy}}& V_{\mathrm{zz}}\end{array}\right]---\left(3\right)$

其中，地球引力位二阶导数是对称张量，同时在真空情况下满足Laplace方程，表现为无迹性，即V_xx+V_yy+V_zz=0，因此在公式（3）中的9个重力梯度分量中有5个是独立的。

$\left\{\begin{array}{c}{V}_{\mathrm{xx}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}{\left(\frac{R_e^3}{r}\right)}^{l+3}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{lm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ})H^{\mathrm{xx}}\left(\mathrm{\θ}\right)\\ V_{\mathrm{yy}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}{\left({\frac{R_e^3}{r}}_{}^{}\right)}^{l+3}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{lm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ})H^{\mathrm{yy}}\left(\mathrm{\θ}\right)\\ V_{\mathrm{zz}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}{\left(\frac{R_e^3}{r}\right)}^{l+3}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{lm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ})H^{\mathrm{zz}}\left(\mathrm{\θ}\right)\\ V_{\mathrm{xy}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}{\left(\frac{R_e^3}{r}\right)}^{l+3}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}(-{\stackrel{\‾}{C}}_{\mathrm{lm}}\sin \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\cos \mathrm{m\λ})H^{\mathrm{xy}}\left(\mathrm{\θ}\right)\\ V_{\mathrm{xz}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}{\left(\frac{R_e^3}{r}\right)}^{l+3}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{lm}}\mathrm{cism\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ})H^{\mathrm{xz}}\left(\mathrm{\θ}\right)\\ V_{\mathrm{yz}}(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}{\left(\frac{R_e^3}{r}\right)}^{l+3}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}(-{\stackrel{\‾}{G}}_{\mathrm{lm}}\sin \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\mathrm{cism\λ})H^{\mathrm{yz}}\left(\mathrm{\θ}\right)\end{array}---\left(4\right)\right.$

其中，

$\left\{\begin{array}{c}{H}^{\mathrm{xx}}\left(\mathrm{\θ}\right)={\stackrel{\‾}{P}}_{\mathrm{lm}}^{\′\′}\left(\cos \mathrm{\θ}\right)-(n+1){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)\\ H^{\mathrm{yy}}\left(\mathrm{\θ}\right)={\tan }^{-1}\mathrm{\θ}{\stackrel{\‾}{P}}_{\mathrm{lm}}^{\′}\left(\cos \mathrm{\θ}\right)-(n+1+m^2{\sin }^{-2}\mathrm{\θ}){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)\\ H^{\mathrm{zz}}\left(\mathrm{\θ}\right)=(l+1)(l+2){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)\\ H^{\mathrm{xy}}\left(\mathrm{\θ}\right)=n{\sin }^{-1}\mathrm{\θ}[{\stackrel{\‾}{P}}_m^{\′}\left(\cos \mathrm{\θ}\right)-{\tan }^{-1}\mathrm{\θ}{\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)]\\ H^{\mathrm{xz}}\left(\mathrm{\θ}\right)=(l+2){\stackrel{\‾}{P}}_{\mathrm{lm}}^{\′}\left(\cos \mathrm{\θ}\right)\\ H^{\mathrm{yz}}\left(\mathrm{\θ}\right)=m(l+2){\sin }^{-1}\mathrm{\θ}{\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)\end{array}---\left(5\right)\right.$

式中的Legendre函数及一阶导数、二阶导数

$\left\{\begin{array}{c}{\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)={\mathrm{\γ}}_m{2}^{-l}{\sin }^m\mathrm{\θ}\underset{k=0}{\overset{[(l-m)/2]}{\mathrm{\Σ}}}{(-1)}^k\frac{(2l-2k)!}{k!(l-k)!(l-m-2k)!}{\left(\cos \mathrm{\θ}\right)}^{l-m-2k}(m\≤l)\\ {\stackrel{\‾}{P}}_{\mathrm{lm}}^{\′}\left(\cos \mathrm{\θ}\right)=\frac{1}{\sin \mathrm{\θ}}[(l+1)\cos \mathrm{\θ}{\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)-(l-m-1){\stackrel{\‾}{P}}_{l+1,m}\left(\cos \right)\mathrm{\θ}]\\ {\stackrel{\‾}{P}}_{\mathrm{lm}}^{\′\′}\left(\cos \mathrm{\θ}\right)=-l{\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right)+l\cos \mathrm{\θ}{\stackrel{\‾}{P}}_{l-1,m}^{\′}\left(\cos \mathrm{\θ}\right)+\frac{1}{4}{\cos }^2\text{\θ[}{\stackrel{\‾}{P}}_{l-1,m+1}^{\′}\left(\cos \mathrm{\θ}\right)\text{-4}{\stackrel{\‾}{P}}_{l-1,m-1}^{\′}\left(\cos \mathrm{\θ}\right)\text{]}\end{array}\right.$

其中， ${\mathrm{\γ}}_m=\left\{\begin{array}{cc}\sqrt{2(2l+1)\frac{l-\left|m\right|!}{(l+|m\left|\right)!}}& (m\≠0)\\ \sqrt{2l+1}& (m=0)\end{array}.\right.$

卫星重力梯度张量V_ab的功率谱表示为

$P^2\left(V_{\mathrm{ab}}\right)=\underset{l=0}{\overset{L}{\mathrm{\Σ}}}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}{[\frac{1}{4\mathrm{\π}}\∫\∫V_{\mathrm{ab}}(r,\mathrm{\φ},\mathrm{\λ}){\stackrel{\‾}{Y}}_{\mathrm{lm}}(\mathrm{\φ},\mathrm{\λ})\cos \mathrm{\φd\φd\λ}]}^2---\left(6\right)$

其中， ${\stackrel{\&OverBar;}{Y}}_{\mathrm{lm}}(\mathrm{\&phi;},\mathrm{\&lambda;})={\stackrel{\&OverBar;}{P}}_{l\left|m\right|}\left(\sin \mathrm{\&phi;}\right)Q_m\left(\mathrm{\&lambda;}\right),$ $Q_m\left(\mathrm{\&lambda;}\right)=\left\{\begin{array}{cc}\cos \mathrm{m\&lambda;}& m\&GreaterEqual;0\\ \sin \left|m\right|\mathrm{\&lambda;}& m<0\end{array},a,b=x,y,z.\right.$

基于公式（4）和（6）以及球函数的正交归一性，卫星重力梯度张量的信号功率谱表示如下

$P^2\left(V_{\mathrm{ab}}\right)={\left(\frac{\mathrm{GM}}{R_e^3}\right)}^2\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}{A}_{\mathrm{ab}}^2{\left(\frac{R_e}{R_e+H}\right)}^{2l+6}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}^2+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}^2)---\left(7\right)$

其中，A_ab表示敏感度系数，H表示卫星轨道高度，和可由德国慕尼黑工业大学公布的地球重力场模型GO_CONS_GCF_2_TIM_R2获得。

基于公式（7），卫星重力梯度全张量V_xyz的信号功率谱表示如下

P²(V_xyz)=P²(V_xx)+P²(V_yy)+P²(V_zz)+P²(V_xz)    （8）

图1表示基于公式（4）获得的卫星重力梯度张量V_xx,V_yy,V_zz,V_xz,V_xyz的敏感度系数表达式A_xx,A_yy,A_zz,A_xz,A_xyz。图2表示在各阶处卫星重力梯度张量的敏感度系数|A_ab|。研究结果表明：第一，卫星重力梯度垂向张量V_zz是最主要分量，对地球重力场反演精度最敏感；第二，卫星重力梯度水平张量V_xx,V_yy是保证地球重力场反演精度的重要组成部分；第三，卫星重力梯度交叉张量V_xz对地球重力场反演精度的贡献相对较小。

卫星重力梯度张量的信号幅度谱（每阶）如图3所示，细虚线、粗实线、粗虚线、细实线、圆圈线和十字线分别表示Kaula垂向张量V_Kzz和卫星重力梯度张量V_xx,V_yy,V_zz,V_xz,V_xyz的信号幅度谱，其中卫星轨道高度H=250km，引力质量常数GM=3.986004415×10¹⁴N·m²/g²，地球平均半径R_e=6378km。研究结果表明：第一，通过Kaula垂向张量V_Kzz和卫星重力梯度垂向张量V_zz信号幅度谱在各阶处的符合性可有效验证卫星重力梯度张量的信号功率谱公式（7）的正确性；第二，卫星重力梯度垂向分量V_zz信号最强，水平分量V_xx,V_yy次之，交叉分量V_xz最弱；第三，卫星重力梯度对角张量V_xx,V_yy,V_zz是反演高精度和高空间分辨率地球重力场的必备分量。

基于Kaula规则，卫星重力梯度张量的信号功率谱表示如下

$P_K^2\left(V_{\mathrm{ab}}\right)={\left(\frac{\mathrm{GM}}{R_e^3}\right)}^2\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}{A}_{\mathrm{ab}}^2{\left(\frac{R_e}{R_e+H}\right)}^{2l+6}(2l+1)\frac{{10}^{-10}}{l^4}---\left(9\right)$

步骤三：建立卫星重力梯度反演半解析误差模型

3.1）建立卫星重力梯度张量的半解析误差模型

基于公式（4）和（6）以及球函数的正交归一性，卫星重力梯度张量误差δV_ab的功率谱可表示如下

$P^2\left({\mathrm{\&delta;V}}_{\mathrm{ab}}\right)={\left(\frac{\mathrm{GM}}{R_e^3}\right)}^2\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}{A}_{\mathrm{ab}}^2{\left(\frac{R_e}{R_e+H}\right)}^{2l+6}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\right)}^2+-{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\right)}^2---\left(10\right)$

其中，表示地球引力位系数精度。

累积大地水准面误差公式表示如下

${\mathrm{\&sigma;}}_N^L=R_e\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\right)}^2+{\left({\mathrm{\&delta;}\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\right)}^2}---\left(11\right)$

联合公式（10）和（11），卫星重力梯度张量的累积大地水准面误差公式表示如下

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;}{V}_{\mathrm{ab}}\right)=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{A_{\mathrm{ab}}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{ab}}\right)}---\left(12\right)$

其中，σ²(δV_ab)表示卫星重力梯度张量的方差。T表示卫星重力梯度观测总时间，Δt表示卫星重力梯度观测数据的采样间隔；TΔt表示卫星重力梯度观测值的数量，据统计学原理可知，如果卫星重力梯度观测值的数量增加TΔt倍，地球重力场反演精度提高约倍。

基于公式（12）和图1，卫星重力梯度张量V_xx,V_yy,V_zz,V_xz,V_xyz的累积大地水准面误差公式分别表示如下

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;}{V}_{\mathrm{xx}}\right)=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{xx}}\right)}---\left(13\right)$

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;}{V}_{\mathrm{yy}}\right)=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{{[(l+1)(l+2)-\sqrt{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}}]}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{yy}}\right)}---\left(14\right)$

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;}{V}_{\mathrm{zz}}\right)=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{{(l+1)}^2{(l+2)}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{zz}}\right)}---\left(15\right)$

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;}{V}_{\mathrm{xz}}\right)=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{{[(l+1)/5]}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{xz}}\right)}---\left(16\right)$

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;}{V}_{\mathrm{xyz}}\right)=\frac{R_e^4}{\mathrm{GM}\sqrt{T/\mathrm{\&Delta;t}}}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}+{[(l+1)(l+2)-\sqrt{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}}]}^2+{(l+1)}^2{(l+2)}^2+{[(l+1)/5]}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{xyz}}\right)}---\left(17\right)$

3.2）建立卫星轨道位置的半解析误差模型

卫星向心加速度和卫星轨道位置r之间的关系表示如下

$\stackrel{\&CenterDot;\&CenterDot;}{r}=\frac{\mathrm{GM}}{r^2}---\left(18\right)$

在公式（18）两边同除r可得

$\frac{\stackrel{\&CenterDot;\&CenterDot;}{r}}{r}=\frac{\mathrm{GM}}{r^3}---\left(19\right)$

其中，表示卫星重力梯度。

基于功率谱原理，并在公式（19）两边同时微分可得

$P^2\left(\mathrm{\&delta;}{V}_{\mathrm{xyz}}\right)={(-\frac{3\mathrm{GM}}{r^4})}^2{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;r}\right)---\left(20\right)$

其中，σ²(δr)表示卫星轨道位置的方差

$P^2\left(\mathrm{\&delta;}{V}_{\mathrm{xyz}}\right)=\frac{{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{xyz}}\right)}{L_{\max }}---\left(21\right)$

σ²(δV_xyz)表示卫星重力梯度张量的方差，L_max表示地球重力场理论上可反演的最高阶数，由于地球重力场的部分高频信号湮没于观测误差，因此实测最高阶数将低于理论值

$L_{\max }=\frac{\mathrm{\&pi;r}}{D}---\left(22\right)$

其中，表示半波长空间分辨率，表示卫星平均速度。

联合公式（20）～（22），卫星重力梯度张量误差δV_xyz和轨道位置误差δr之间的转换关系表示如下

$\mathrm{\&delta;}{V}_{\mathrm{xyz}}=\sqrt{\frac{9\mathrm{GM}{\mathrm{\&pi;}}^2}{r^5\mathrm{\&Delta;}{t}^2}}\mathrm{\&delta;r}---\left(23\right)$

基于公式（17）和（23），轨道位置误差δr影响累计大地水准面精度的半解析误差模型表示如下：

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;r}\right)=\frac{R_e^4}{\mathrm{GM}\sqrt{T/\mathrm{\&Delta;t}}}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}+{[(l+1)(l+2)-\sqrt{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}}]}^2+{(l+1)}^2{(l+2)}^2+{[(l+1)/5]}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\sqrt{\frac{9\mathrm{GM}{\mathrm{\&pi;}}^2}{r^5\mathrm{\&Delta;}{t}^2}}\mathrm{\&delta;r}\right)}---\left(24\right)$

3.3）建立卫星重力梯度和轨道位置的联合半解析误差模型

基于公式（17）和（24），GOCE卫星重力梯度仪的重力梯度张量误差和GPS/GLONASS接收机的轨道位置误差影响累计大地水准面精度的联合半解析误差模型表示如下

${\mathrm{\&sigma;}}_N(\mathrm{\&delta;}{V}_{\mathrm{xyz}},\mathrm{\&delta;r})=\frac{R_e^4}{\mathrm{GM}\sqrt{T/\mathrm{\&Delta;t}}}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}+{[(l+1)(l+2)-\sqrt{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}}]}^2+{(l+1)}^2{(l+2)}^2+{[(l+1)/5]}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;\&eta;}\right)}---\left(25\right)$

其中，表示GOCE重力梯度卫星关键载荷的总误差，σ²(δV_xyz)表示卫星重力梯度张量方差，表示卫星轨道位置方差。

步骤四：通过卫星重力梯度反演半解析误差模型确定累积大地水准面误差

基于功率谱半解析法，利用2011年的GOCE-Level-1B卫星重力梯度测量数据δV_xyz和卫星轨道位置测量数据δr反演累计大地水准面误差的过程如下

第一，首先以0.1°×0.1°为网格分辨率，在地球表面的经度（0°~360°）和纬度（-90°~90°）范围内绘制网格；其次，按照GOCE卫星轨道在地球表面的轨迹点位置依次加入δη；最后，将分布于地球表面的δη平均归算于划分的网格点δη(φ，λ)处。

第二，将δη(φ，λ)按球谐函数展开为

$\mathrm{\&delta;\&eta;}(\mathrm{\&phi;},\mathrm{\&lambda;})=\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}\left[(C_{{\mathrm{\&delta;\&eta;}}_{\mathrm{lm}}}\cos \mathrm{m\&lambda;}+S_{{\mathrm{\&delta;\&eta;}}_{\mathrm{lm}}}\sin \mathrm{m\&lambda;}){\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\sin \mathrm{\&phi;}\right)\right]---\left(26\right)$

其中，表示δη(φ，λ)按球函数展开的系数

$(C_{{\mathrm{\&delta;\&eta;}}_{\mathrm{lm}}},S_{{\mathrm{\&delta;\&eta;}}_{\mathrm{lm}}})=\frac{1}{4\mathrm{\&pi;}}\&Integral;\&Integral;\left[\mathrm{\&delta;\&eta;}(\mathrm{\&phi;},\mathrm{\&lambda;}){\stackrel{\&OverBar;}{Y}}_{\mathrm{lm}}(\mathrm{\&phi;},\mathrm{\&lambda;})\cos \mathrm{\&phi;d\&phi;d\&lambda;}\right]---\left(27\right)$

δη在各阶处的方差表示为

${\mathrm{\&sigma;}}_l^2\left(\mathrm{\&delta;\&eta;}\right)=\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}(C_{\mathrm{\&delta;}{\mathrm{\&eta;}}_{\mathrm{lm}}}^2+S_{{\mathrm{\&delta;\&eta;}}_{\mathrm{lm}}}^2)---\left(28\right)$

将公式（28）代入（25），可有效和快速地反演全球重力场的精度。

步骤五：卫星重力梯度反演半解析误差模型准确性的验证

基于GOCE关键载荷匹配精度指标分别估计出的累计大地水准面精度如图4所示，实线、虚线和星号线分别表示单独引入GOCE卫星重力梯度仪的重力梯度张量误差V_xyz和GPS/GLONASS复合接收机的轨道位置误差r、以及联合误差V_xyz+r估计累计大地水准面的精度，其中GOCE卫星重力梯度测量精度3×10^-12/s²、卫星轨道位置测量精度1×10^-2m、平均轨道高度H=250km、地球平均半径R_e=6370km、观测时间T=8个月、采样间隔Δt=5s、地球引力常数GM=3.986004415×10¹⁴Nm²/kg。据图中实线和虚线在各阶处的符合性，可验证本发明提出的GOCE各项关键载荷精度指标是匹配的。同时，通过本发明提出的GOCE卫星关键载荷匹配精度指标和欧空局公布的GOCE-Level-1B实测精度指标的符合性，充分证明了本发明建立的卫星重力梯度和卫星轨道位置的单独半解析误差模型是可靠和匹配的。在250阶处，基于卫星重力梯度张量和轨道位置联合半解析误差模型，估计累计大地水准面精度为1.769×10^-1m，其结果与德国慕尼黑工业大学公布的GO_CONS_GCF_2_TIM_R2（采用8个月的GOCE卫星重力梯度观测数据）地球重力场模型精度1.760×10^-1m符合较好，进而证明了本发明建立的联合半解析误差模型是正确的。

在另一优选的实施例中，通过上述步骤五中对卫星重力梯度反演半解析误差模型的准确性进行验证之后，再进行步骤四的通过卫星重力梯度反演半解析误差模型确定累积大地水准面误差。

步骤六：通过卫星重力梯度反演半解析误差模型进行GOCE Follow-On重力梯度卫星需求论证

6.1）卫星重力梯度测量精度影响

基于不同卫星重力梯度测量精度所估计出的累计大地水准面精度如图5所示，基于不同卫星重力梯度测量精度3×10^-12/s²、3×10^-13/s²、3×10^-14/s²和3×10^-15/s²，分别估计了300阶GOCE Follow-On累计大地水准面精度。研究结果表明：在300阶处，基于卫星重力梯度测量精度3×10^-12/s²，估计GOCE Follow-On累计大地水准面精度为6.068×10^-1m；如果采用卫星重力梯度测量精度3×10^-13/s²、3×10^-14/s²和3×10^-15/s²，估计累计大地水准面精度将分别提高10倍、100倍和1000倍。由于星载重力梯度仪的测量精度是决定地球重力场反演精度的最主要因素，因此，如果GOCE Follow-On卫星重力梯度计划采用冷原子干涉重力梯度仪（测量精度10^-13/s²～10^-15/s²），其地球重力场的感测精度较当前GOCE重力梯度卫星至少可提高一个数量级。

6.2）卫星轨道位置测量精度影响

图6表示分别基于卫星轨道位置测量精度10^-2m、10^-3m、10^-4m和10^-5m，估计300阶GOCE Follow-On累计大地水准面精度。研究结果表明：在300阶处，基于卫星轨道位置测量精度10^-2m估计地球重力场精度为6.730×10^-1m；如果卫星轨道位置测量精度分别提高10倍、100倍和1000倍，估计地球重力场精度呈线性升高趋势。由于卫星重力梯度测量对卫星轨道位置精度敏感性较低，而且目前国际全球导航系统（美国GPS、俄罗斯GLONASS、中国北斗、欧洲Galileo等）的最优绝对定轨精度仅为cm级，因此，建议GOCE Follow-On卫星重力梯度计划的定轨精度设计为1～0.1cm。

6.3）卫星轨道高度影响

基于不同卫星轨道高度所估计出的累计大地水准面精度如图7所示，粗实线、粗虚线、细虚线和细实线分别表示基于不同卫星轨道高度200km、250km、300km和350km，估计300阶GOCE Follow-On累计大地水准面精度。研究结果表明：在300阶处，基于卫星轨道高度200km估计累计大地水准面精度为1.049×10^-1m；如果卫星轨道高度升高到250km、300km和350km，估计累计大地水准面精度分别降低了8.639倍、75.491倍和660.819倍。由于地球引力位随卫星轨道高度的升高呈指数衰减，因此，有效降低卫星轨道高度是反演高精度和高空间分辨率地球重力场的关键技术。但随着轨道高度每降低100km，作用于重力卫星的空气阻力约提高一个数量级。因此，虽然GOCE Follow-On重力梯度卫星携带了非保守力补偿系统（Drag-free），但卫星轨道高度的优化选取至关重要。综上所述，建议GOCE Follow-On重力梯度卫星的轨道高度选择在200～300km。

以上具体实施方式仅为本发明的一种实施示例，其描述较为具体和详细，但不能因此而理解为对本发明专利范围的限制。其具体实施步骤顺序和模型参数可根据实际需要进行相应的调整。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。

Claims (8)

1.一种基于功率谱半解析的卫星重力梯度反演方法，其特征在于包含以下步骤：

步骤1，通过星载重力梯度仪采集卫星重力梯度测量数据δV_xyz，基于星载GPS/GLONASS复合接收机采集卫星轨道位置测量数据δr；

步骤2，建立卫星重力梯度的信号功率谱分析模型，并对所采集的卫星重力梯度测量数据进行信号功率谱敏感度分析；

步骤3，通过卫星重力梯度张量误差的功率谱分析和卫星轨道位置误差的功率谱分析建立卫星重力梯度反演半解析误差模型；其中，所述步骤3包括：

步骤3.1，通过卫星重力梯度仪的重力梯度张量误差对累计大地水准面精度的影响建立卫星重力梯度张量的半解析误差模型；

步骤3.2，通过GPS/GLONASS复合接收机的轨道位置误差对累计大地水准面精度的影响建立卫星轨道位置的半解析误差模型；

步骤3.3，通过卫星重力梯度张量的半解析误差模型和卫星轨道位置的半解析误差模型，建立卫星重力梯度和轨道位置的联合半解析误差模型，以此作为卫星重力梯度反演半解析误差模型；

步骤4，使用所述卫星重力梯度反演半解析误差模型，以及采集得到的卫星重力梯度测量数据δV_xyz和卫星轨道位置测量数据δr反演累计大地水准面误差；

其中所述步骤2为：

在球坐标系中，地球引力位按球谐函数展开的表达式为

$V(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{GM}}{R_e}\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\cos \mathrm{m\&lambda;}+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\sin \mathrm{m\&lambda;}){\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\&theta;}\right)---\left(1\right)$

其中，GM表示地球质量M和万有引力常数G之积，R_e表示地球的平均半径，L表示地球引力位按球函数展开的最大阶数；表示卫星轨道位置，x,y,z分别表示卫星位置矢量在直角坐标系中的三个分量，θ和λ表示地心余纬度和经度；表示规格化的Legendre函数，l表示阶数，m表示次数；和表示待求的规格化引力位系数；

球坐标系r,θ,λ和直角坐标系x,y,z的互换公式表示为

$\left\{\begin{array}{c}r=\sqrt{x^2+y^2+z^2}\\ \sin \mathrm{\&theta;}=\frac{\sqrt{x^2+y^2}}{\sqrt{x^2+y^2+z^2}}\\ \cos \mathrm{\&theta;}=\frac{z}{\sqrt{x^2+y^2+z^2}}\\ \sin \mathrm{\&lambda;}=\frac{y}{\sqrt{x^2+y^2}}\\ \cos \mathrm{\&lambda;}=\frac{x}{\sqrt{x^2+y^2}}\end{array}\right.---\left(2\right)$

基于公式(1)和(2)，在直角坐标系中，地球引力位V分别对x,y,z的二阶导数表示如下

$\mathrm{\&Gamma;}=\left[\begin{array}{ccc}{V}_{\mathrm{xx}}& V_{\mathrm{xy}}& V_{\mathrm{xz}}\\ V_{\mathrm{yx}}& V_{\mathrm{yy}}& V_{\mathrm{yz}}\\ V_{\mathrm{zx}}& V_{\mathrm{zy}}& V_{\mathrm{zz}}\end{array}\right]---\left(3\right)$

其中，地球引力位二阶导数是对称张量，同时在真空情况下满足Laplace方程，表现为无迹性，即V_xx+V_yy+V_zz＝0，因此在公式(3)中的9个重力梯度分量中有5个是独立的；

$\left\{\begin{array}{c}{V}_{\mathrm{xx}}(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}{\left(\frac{R_e^3}{r}\right)}^{l+3}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\cos \mathrm{m\&lambda;}+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\sin \mathrm{m\&lambda;})H^{\mathrm{xx}}\left(\mathrm{\&theta;}\right)\\ V_{\mathrm{yy}}(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}{\left(\frac{R_e^3}{r}\right)}^{l+3}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\cos \mathrm{m\&lambda;}+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\sin \mathrm{m\&lambda;})H^{\mathrm{yy}}\left(\mathrm{\&theta;}\right)\\ V_{\mathrm{zz}}(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}{\left(\frac{R_e^3}{r}\right)}^{l+3}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\cos \mathrm{m\&lambda;}+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\sin \mathrm{m\&lambda;})H^{\mathrm{zz}}\left(\mathrm{\&theta;}\right)\\ V_{\mathrm{xy}}(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}{\left(\frac{R_e^3}{r}\right)}^{l+3}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}({-\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\sin \mathrm{m\&lambda;}+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\cos \mathrm{m\&lambda;})H^{\mathrm{xy}}\left(\mathrm{\&theta;}\right)\\ V_{\mathrm{xz}}(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}{\left(\frac{R_e^3}{r}\right)}^{l+3}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\cos \mathrm{m\&lambda;}+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\sin \mathrm{m\&lambda;})H^{\mathrm{xz}}\left(\mathrm{\&theta;}\right)\\ V_{\mathrm{yz}}(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{GM}}{R_e^3}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}{\left(\frac{R_e^3}{r}\right)}^{l+3}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}(-{\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\sin \mathrm{m\&lambda;}+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\cos \mathrm{m\&lambda;})H^{\mathrm{yz}}\left(\mathrm{\&theta;}\right)\end{array}\right.---\left(4\right)$

其中，

$\left\{\begin{array}{c}{H}^{\mathrm{xx}}\left(\mathrm{\&theta;}\right)={\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}^{\&prime;\&prime;}\left(\cos \mathrm{\&theta;}\right)-(l+1){\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\&theta;}\right)\\ H^{\mathrm{yy}}\left(\mathrm{\&theta;}\right)={\tan }^{-1}\mathrm{\&theta;}{\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}^{\&prime;}\left(\cos \mathrm{\&theta;}\right)-(l+1+m^2{\sin }^{-2}\mathrm{\&theta;}){\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\&theta;}\right)\\ H^{\mathrm{zz}}\left(\mathrm{\&theta;}\right)=(l+1)(l+2){\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\&theta;}\right)\\ H^{\mathrm{xy}}\left(\mathrm{\&theta;}\right)=m{\sin }^{-1}\mathrm{\&theta;}[{\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}^{\&prime;}\left(\cos \mathrm{\&theta;}\right)-{\tan }^{-1}\mathrm{\&theta;}{\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\&theta;}\right)]\\ H^{\mathrm{xz}}\left(\mathrm{\&theta;}\right)=(l+2){\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}^{\&prime;}\left(\cos \mathrm{\&theta;}\right)\\ H^{\mathrm{yz}}\left(\mathrm{\&theta;}\right)=m(l+2){\sin }^{-1}\mathrm{\&theta;}{\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\&theta;}\right)\end{array}\right.---\left(5\right)$

式中的Legendre函数及一阶导数、二阶导数

$\left\{\begin{array}{cc}{\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\&theta;}\right)={\mathrm{\&gamma;}}_m{2}^{-l}{\sin }^m\mathrm{\&theta;}\underset{k=0}{\overset{[(l-m)/2]}{\mathrm{\&Sigma;}}}{(-1)}^k\frac{(2l-2k)!}{k!(l-k)!(l-m-2k)!}{\left(\cos \mathrm{\&theta;}\right)}^{l-m-2k}& (m\&le;l)\\ {\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}^{\&prime;}\left(\cos \mathrm{\&theta;}\right)=\frac{1}{\sin \mathrm{\&theta;}}[(l+1)\cos \mathrm{\&theta;}{\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\&theta;}\right)-(l-m-1){\stackrel{\&OverBar;}{P}}_{l+1,m}\left(\cos \mathrm{\&theta;}\right)]& \\ {\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}^{\&prime;\&prime;}\left(\cos \mathrm{\&theta;}\right)=-l{\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\&theta;}\right)+l\cos \mathrm{\&theta;}{\stackrel{\&OverBar;}{P}}_{l-1,m}^{\&prime;}\left(\cos \mathrm{\&theta;}\right)+\frac{l}{4}{\cos }^2\mathrm{\&theta;}[{\stackrel{\&OverBar;}{P}}_{l-1,m+1}^{\&prime;}\left(\cos \mathrm{\&theta;}\right)-4{\stackrel{\&OverBar;}{P}}_{l-1,m-1}^{\&prime;}\left(\cos \mathrm{\&theta;}\right)]& \end{array}\right.$

其中， ${\mathrm{\&gamma;}}_m=\left\{\begin{array}{cc}\sqrt{2(1l+1)\frac{(l-|m\left|\right)!}{(l+|m\left|\right)!}}& (m\&NotEqual;0)\\ \sqrt{2l+1}& (m=0)\end{array}\right.;$

卫星重力梯度张量V_ab的功率谱表示为

$P^2\left(V_{\mathrm{ab}}\right)=\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}{[\frac{1}{4\mathrm{\&pi;}}\&Integral;\&Integral;V_{\mathrm{ab}}(r,\mathrm{\&phi;},\mathrm{\&lambda;}){\stackrel{\&OverBar;}{Y}}_{\mathrm{lm}}(\mathrm{\&phi;},\mathrm{\&lambda;})\cos \mathrm{\&phi;d\&phi;d\&lambda;}]}^2---\left(6\right)$

其中， ${\stackrel{\&OverBar;}{Y}}_{\mathrm{lm}}(\mathrm{\&phi;},\mathrm{\&lambda;})={\stackrel{\&OverBar;}{P}}_{l\left|m\right|}\left(\sin \mathrm{\&phi;}\right)Q_m\left(\mathrm{\&lambda;}\right),$ $Q_m\left(\mathrm{\&lambda;}\right)\left\{\begin{array}{cc}\cos \mathrm{m\&lambda;}& m\&GreaterEqual;0\\ \sin \left|m\right|\mathrm{\&lambda;}& m<0\end{array}\right.,a,b=x,y,z;$

基于公式(4)和(6)以及球函数的正交归一性，卫星重力梯度张量的信号功率谱表示如下

$P^2\left(V_{\mathrm{ab}}\right)={\left(\frac{\mathrm{GM}}{R_e^3}\right)}^2\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}{A}_{\mathrm{ab}}^2{\left(\frac{R_e}{R_e+H}\right)}^{2l+6}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}^2+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}^2)---\left(7\right)$

其中，A_ab表示敏感度系数，H表示卫星轨道高度，和可由德国慕尼黑工业大学公布的地球重力场模型GO_CONS_GCF_2_TIM_R2获得；

基于公式(7)，卫星重力梯度全张量V_xyz的信号功率谱表示如下

P²(V_xyz)＝P²(V_xx)+P²(V_yy)+P²(V_zz)+P²(V_xz)   (8)

基于Kaula规则，卫星重力梯度张量的信号功率谱表示如下

$P_K^2\left(V_{\mathrm{ab}}\right)={\left(\frac{\mathrm{GM}}{R_e^3}\right)}^2\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}{A}_{\mathrm{ab}}^2{\left(\frac{R_e}{R_e+H}\right)}^{2l+6}(2l+1)\frac{{10}^{-10}}{l^4}---\left(9\right).$

2.如权利要求1所述的基于功率谱半解析的卫星重力梯度反演方法，其特征在于：所述步骤3.1为：

基于公式(4)和(6)以及球函数的正交归一性，卫星重力梯度张量误差δV_ab的功率谱表示如下

$P^2\left({\mathrm{\&delta;V}}_{\mathrm{ab}}\right)={\left(\frac{\mathrm{GM}}{R_e^3}\right)}^2\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}{A}_{\mathrm{ab}}^2{\left(\frac{R_e}{R_e+H}\right)}^{2l+6}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\right)}^2+{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\right)}^2---\left(10\right)$

其中，表示地球引力位系数精度；

累积大地水准面误差公式表示如下

${\mathrm{\&sigma;}}_N^L=R_e\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\right)}^2+{\left(\mathrm{\&delta;}{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\right)}^2}---\left(11\right)$

联合公式(10)和(11)，卫星重力梯度张量的累积大地水准面误差公式表示如下

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;}{V}_{\mathrm{ab}}\right)=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{A_{\mathrm{ab}}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{ab}}\right)}---\left(12\right)$

其中，σ²(δV_ab)表示卫星重力梯度张量的方差；T表示卫星重力梯度观测总时间，Δt表示卫星重力梯度观测数据的采样间隔，T/Δt表示卫星重力梯度观测值的数量，据统计学原理可知，如果卫星重力梯度观测值的数量增加T/Δt倍，地球重力场反演精度提高倍；

基于公式(12)，卫星重力梯度张量V_xx,V_yy,V_zz,V_xz,V_xyz的累积大地水准面误差公式表示如下

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;}{V}_{\mathrm{xx}}\right)=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{xx}}\right)}---\left(13\right)$

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;}{V}_{\mathrm{yy}}\right)=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{{[(l+1)(l+2)-\sqrt{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}}]}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{yy}}\right)}---\left(14\right)$

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;}{V}_{\mathrm{zz}}\right)=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{{(l+1)}^2{(l+2)}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{zz}}\right)}---\left(15\right)$

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;}{V}_{\mathrm{xz}}\right)=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{{[l+1/5]}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{xz}}\right)}---\left(16\right)$

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;}{V}_{\mathrm{xyz}}\right)=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}+{[(l+1)(l+2)-\sqrt{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}}]}^2+{(l+1)}^2{(l+2)}^2+{[(l+1)/5]}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{xyz}}\right)}---\left(17\right)$

所述步骤3.2为：

卫星向心加速度和卫星轨道位置r之间的关系表示如下

$\stackrel{\&CenterDot;\&CenterDot;}{r}=\frac{\mathrm{GM}}{r^2}---\left(18\right)$

在公式(18)两边同除r可得

$\frac{\stackrel{\&CenterDot;\&CenterDot;}{r}}{r}\frac{\mathrm{GM}}{r^3}---\left(19\right)$

其中，表示卫星重力梯度；

基于功率谱原理，并在公式(19)两边同时微分可得

$P^2\left(\mathrm{\&delta;}{V}_{\mathrm{xyz}}\right)={(-\frac{3\mathrm{GM}}{r^4})}^2{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;r}\right)---\left(20\right)$

其中，σ²(δr)表示卫星轨道位置的方差

$P^2\left(\mathrm{\&delta;}{V}_{\mathrm{xyz}}\right)=\frac{{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;}{V}_{\mathrm{xyz}}\right)}{L_{\max }}---\left(21\right)$

σ²(δV_xyz)表示卫星重力梯度张量的方差，L_max表示地球重力场理论上可反演的最高阶数，由于地球重力场的部分高频信号湮没于观测误差，因此实测最高阶数将低于理论值

$L_{\max }=\frac{\mathrm{\&pi;r}}{D}---\left(2\right)$

其中，表示半波长空间分辨率，表示卫星平均速度，Δt表示卫星重力梯度观测值的采样间隔；

联合公式(20)～(22)，卫星重力梯度张量误差δV_xyz和轨道位置误差δr之间的转换关系表示如下

$\mathrm{\&delta;}{V}_{\mathrm{xyz}}=\sqrt{\frac{9\mathrm{GM}{\mathrm{\&pi;}}^2}{r^5\mathrm{\&Delta;}{t}^2}}\mathrm{\&delta;r}---\left(23\right)$

基于公式(17)和(23)，轨道位置误差δr影响累计大地水准面精度的半解析误差模型表示如下

${\mathrm{\&sigma;}}_N\left(\mathrm{\&delta;r}\right)=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}+{[(l+1)(l+2)-\sqrt{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}}]}^2+{(l+1)}^2{(l+2)}^2+{[(l+1)/5]}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\sqrt{\frac{9\mathrm{GM}{\mathrm{\&pi;}}^2}{r^5\mathrm{\&Delta;}{t}^2}\mathrm{\&delta;r}}\right)}---\left(24\right);$

所述步骤3.3为：

基于公式(17)和(24)，卫星重力梯度仪的重力梯度张量误差和GPS/GLONASS接收机的轨道位置误差影响累计大地水准面精度的联合半解析误差模型表示如下

${\mathrm{\&sigma;}}_N(\mathrm{\&delta;}{V}_{\mathrm{xyz}},\mathrm{\&delta;r})=\frac{R_e^4}{\mathrm{GM}\left(\sqrt{T/\mathrm{\&Delta;t}}\right)}\sqrt{\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\frac{2l+1}{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}+{[(l+1)(l+2)-\sqrt{\frac{{(l+1)}^3(l+2)(2l+3)}{3(2l+1)}}]}^2+{(l+1)}^2{(l+2)}^2+{[(l+1)/5]}^2}{\left(\frac{R_e+H}{R_e}\right)}^{2l+6}{\mathrm{\&sigma;}}^2\left(\mathrm{\&delta;\&eta;}\right)}---\left(25\right)$

其中，表示重力梯度卫星关键载荷的总误差，σ²(δV_xyz)表示卫星重力梯度张量方差，表示卫星轨道位置方差。

3.如权利要求2所述的基于功率谱半解析的卫星重力梯度反演方法，其特征在于：所述步骤4具体包括以下步骤：

步骤4.1，首先确定网格分辨率为0.1°×0.1°，在地球表面的经度和纬度范围内绘制网格；

步骤4.2，按照重力梯度卫星轨道在地球表面的轨迹点位置依次加入重力梯度卫星关键载荷的总误差δη；

步骤4.3，将分布于地球表面的δη平均归算于划分的网格点δη(φ,λ)处；

步骤4.4，将δη(φ,λ)按球谐函数展开为

$\mathrm{\&delta;\&eta;}(\mathrm{\&phi;},\mathrm{\&lambda;})=\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}\left[(C_{\mathrm{\&delta;}{\mathrm{\&eta;}}_{\mathrm{lm}}}\cos \mathrm{m\&lambda;}+S_{\mathrm{\&delta;}{\mathrm{\&eta;}}_{\mathrm{lm}}}\sin m\mathrm{\&lambda;}){\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\sin \mathrm{\&phi;}\right)\right]---\left(26\right)$

其中，表示δη(φ,λ)按球函数展开的系数

$(C_{\mathrm{\&delta;}{\mathrm{\&eta;}}_{\mathrm{lm}}},S_{\mathrm{\&delta;}{\mathrm{\&eta;}}_{\mathrm{lm}}})=\frac{1}{4\mathrm{\&pi;}}\&Integral;\&Integral;\left[\mathrm{\&delta;\&eta;}(\mathrm{\&phi;},\mathrm{\&lambda;}){\stackrel{\&OverBar;}{Y}}_m(\mathrm{\&phi;},\mathrm{\&lambda;})\cos \mathrm{\&phi;d\&phi;d\&lambda;}\right]---\left(27\right)$

δη在各阶处的方差表示为

${\mathrm{\&sigma;}}_l^2\left(\mathrm{\&delta;\&eta;}\right)=\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}(C_{\mathrm{\&delta;}{\mathrm{\&eta;}}_{\mathrm{lm}}}^2+S_{\mathrm{\&delta;}{\mathrm{\&eta;}}_{\mathrm{lm}}}^2)---\left(28\right)$

将公式(28)代入(25)，即可反演全球重力场的精度。

4.如权利要求3所述的基于功率谱半解析的卫星重力梯度反演方法，其特征在于：所述重力梯度卫星为GOCE重力梯度卫星或者GOCE Follow-On重力梯度卫星。

5.一种对如权利要求1-4中任意一项所述的基于功率谱半解析的卫星重力梯度反演方法中所使用的卫星重力梯度反演半解析误差模型进行验证的方法，其特征在于还包括如下步骤：

步骤5，使用所述卫星重力梯度反演半解析误差模型确定300阶GOCE地球重力场精度，并通过与国际公布的GO_CONS_GCF_2_TIM_R2重力场模型精度对比，验证所述卫星重力梯度反演半解析误差模型的准确性。

6.如权利要求5所述的对卫星重力梯度反演半解析误差模型进行验证的方法，其特征在于：在所述步骤5中还包括，使用卫星重力梯度和卫星轨道位置的单独半解析误差模型计算GOCE卫星关键载荷匹配精度指标，将计算得到的该指标与欧空局GOCE-Level-1B实测精度指标进行对比，验证卫星重力梯度和卫星轨道位置的单独半解析误差模型。

7.一种使用如权利要求1-4中任意一项所述的基于功率谱半解析的卫星重力梯度反演方法来确定GOCE Follow-On重力梯度卫星参数需求的方法，其特征在于还包括如下步骤：

步骤6，根据地球科学各相关学科需求确定所需的地球重力场反演精度；通过所述卫星重力梯度反演半解析误差模型，依据地球重力场需求精度对GOCEFollow-On重力梯度卫星的关键载荷精度指标和轨道参数进行需求分析。

8.如权利要求7所述的确定GOCE Follow-On重力梯度卫星参数需求的方法，其特征在于：所述步骤6包括：

步骤6.1，分析卫星重力梯度测量精度需求：基于不同卫星重力梯度测量精度3×10^-12/s²、3×10^-13/s²、3×10^-14/s²和3×10^-15/s²，使用卫星重力梯度反演半解析误差模型分别估计300阶GOCE Follow-On累计大地水准面精度，通过卫星重力梯度测量精度与累计大地水准面精度的关系确定星载重力梯度仪的测量精度；

步骤6.2，分析卫星轨道位置测量精度需求：分别基于卫星轨道位置测量精度10^-2m、10^-3m、10^-4m和10^-5m，使用卫星重力梯度反演半解析误差模型估计300阶GOCE Follow-On累计大地水准面精度，通过卫星轨道位置测量精度与累计大地水准面精度的关系确定GOCE Follow-On卫星重力梯度计划的定轨精度；

步骤6.3，分析卫星轨道高度需求：基于不同卫星轨道高度200km、250km、300km和350km，使用卫星重力梯度反演半解析误差模型估计300阶GOCEFollow-On累计大地水准面精度，根据卫星轨道高度与累计大地水准面精度之间的关系以及卫星轨道高度与卫星所遇空气阻力之间的关系确定GOCEFollow-On重力梯度卫星的轨道高度。