CN104750983A - 一种空间分层网格扰动引力场模型构建与扰动引力快速确定方法 - Google Patents

Abstract

本发明公开了一种空间分层网格扰动引力场模型构建与扰动引力快速确定方法，克服了现有技术中，扰动引力场仍需进一步研究的问题。该发明含有以下步骤：一、高度层面划分，即根据全球外空扰动引力变化特征，将外部空间划分成若干个区间。二、网格单元逼近与网格划分准则的制定，即在不同区间，用基于广域多项式的网格单元分别逼近，并分析精度。三、利用准则进行扰动引力模型构建。四、利用模型进行扰动引力实时解算，结合高速飞行的飞行器的应用需求，使得扰动引力的计算可高效、高精度进行。与现有技术相比本发明计算量不大，不需要特别高级的计算机即可实现，无须较多种类的重力测量数据，逼近精度高。

Description

一种空间分层网格扰动引力场模型构建与扰动引力快速确定方法

技术领域

该发明涉及一种空间分层网格扰动引力模型的构建，特别是涉及一种空间分层网格扰动引力场模型构建与扰动引力快速确定方法。

背景技术

扰动引力作为地球重力场中的研究重点，具有非常重要的理论及应用价值。就目前来看，还未曾有人构建过扰动引力场分层网格模型，也无人对扰动引力场的分层网格划分方法进行过研究。

对于扰动引力场建模，国内外主要集中于以下方法：位系数模型、虚拟点质量模型、Stokes积分、单层密度法、球谐函数换极法、梯度法、向上延拓法及函数逼近等方法。其中函数逼近方法较常用的有，B样条函数、BP-神经网络、有限元分析等。这些方法各有优点，但还是存在一些不足：

一、计算量过大。

如位系数模型，计算某个位置完全至2160阶次的某一个分量的扰动引力值，需要首先读取2161×2162＝4672082个位系数模型，并计算个勒让德函数，另外还有超过9.34×10⁶次的乘法计算，计算量之大可见一斑。

二、需要较多种类的重力测量数据。如梯度法、向上延拓法等。其中梯度法需要在计算区域测量足够多的天文、大地及重力数据，因而该方法也只适用于计算局部地区扰动引力；向上延拓法也需要用到地表延拓位置的垂线偏差值。

三、只能计算局部地区扰动引力。如虚拟点质量法、梯度法等。其中在虚拟点质量法构建全球模型时不得不对超大阶数的矩阵进行运算，尤其是求逆运算，一般情况下计算机无法承担如此重任。

四、逼近精度不高。如B样条函数法。该方法只能利用坐标信息中的高程分量，因而无法对三维的扰动引力场进行精确逼近。

五、灵活性较差或是需要事前已知轨迹函数。如BP-神经网络、有限元分析等。其中BP-神经网络需要利用既定轨道附近点的扰动引力值反复进行学习训练，普遍适用性较差，而且训练耗时较长；有限元分析方法可以通过调整格网单元的分辨率，用于逼近不同地区的扰动引力场，可以利用三个分量的坐标信息，逼近精度可以到达较高水平，但是一般情况下该方法需要事先已知轨迹函数。

发明内容

本发明克服了现有技术中，扰动引力场仍需进一步研究的问题，提供一种的计算效果较好的空间分层网格扰动引力场模型构建与扰动引力快速确定方法。

本发明的技术解决方案是，提供一种具有以下步骤的空间分层网格扰动引力场模型构建与扰动引力快速确定方法：一种空间分层网格扰动引力场模型构建，由以下步骤实现：一、高度层面划分，即根据全球外空扰动引力变化特征，将外部空间划分成若干个区间；二、网格单元逼近与网格划分准则的制定，即在不同区间，用基于广域多项式的网格单元分别逼近，并分析精度，将逼近精度较高且数据生成量较小的方案列出，生成扰动引力分层网格模型构建的准则。

所述高度层面划分的分析过程中，利用全球均方值来估计扰动引力分量的谱能量，公式为

其中，C_n为重力异常阶方差，表示如下，

$C_n={\mathrm{\γ}}^2{(n-1)}^2\underset{m=0}{\overset{n}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{nm}}^2+{\stackrel{\‾}{S}}_{\mathrm{nm}}^2)$

为正则化后的模型位系数，此处使用的是EGM2008位模型。

所述区间为4个，按高度范围划分即20-100KM、100-220KM、220-1000KM和1000KM以上。

所述网格单元逼近与网格层面划分准则的建立过程中，将某块研究区域以一定的经纬度和高度间隔划分成若干个单元，形成一系列的矩形单元，称之为主域；再沿着矩形单元的每条边向外延伸，形成次域单元；其节点分别称之为主节点(8个)、次节点(24个)；采用多项式的形式来逼近待求点的扰动引力，如下：g_i(x_i,y_i,z_i)＝f_ia

其中

f_i＝(1,x_i,y_i,z_i,x_iy_i,x_iz_i,y_iz_i,...)

a＝(a₀,a₁,a₂...)^T

在每个节点处可得到一个多项式的值g_i(x,y,z)，而要达到的结果是主节点上的多项式值与扰动引力值相等，而次节点上的多项式的值与扰动引力的差值的平方和最小，即可建立方程，如下：

$\left\{\begin{array}{c}\min M\left(a\right)={(f^{\′}a-s^{\′})}^T(f^{\′}a-s^{\′})\\ \mathrm{fa}-s=0\end{array}\right.$

其中，fa为主节点上的多项式值，s为主节点上的扰动引力值，f′a为次节点上的多项式值，s′为次节点上的扰动引力值；根据最小二乘平差理论，可求得位置参数向量的值；其中节点(包括主节点和次节点)与待求点的坐标要转换到以单元中心为原点的局部坐标系下，转换公式如下：

其中，是主节点的球坐标，是立方体中心点的球坐标。

一种包含空间分层网格扰动引力场模型构建的扰动引力快速确定方法，包含以下步骤：(1)利用扰动引力分层网格模型构建的准则进行扰动引力模型构建，即设计模型构建流程与数据输出格式；(2)利用构建好的扰动引力模型进行扰动引力实时解算，即设计扰动引力计算的流程。

所述利用扰动引力分层网格模型构建的准则进行扰动引力模型的构建过程如下：首先计算扰动引力场节点值，利用准则自适应地选取所需单元的节点并计算单元内的广域多项式系数；将单元合理编号，并将计算得到的系数以一定的次序存入文件中对应的单元编号后；最终就得到了分层网格模型。

所述利用构建好的扰动引力模型进行扰动引力实时解算包括以下步骤：(1)通过飞行器载体上的GPS接收机、惯性导航及其他的地形或重力匹配导航系统，精确确定飞行器位置坐标信息；(2)对坐标信息进行分析，计算得到飞行器所在位置的层号及单元号；利用载体计算机读取生成的模型数据文件，搜索对应层号和单元号的多项式系数值；(3)将飞行器坐标值代入对应多项式中进行扰动引力解算。

与现有技术相比，本发明空间分层网格扰动引力场模型构建与扰动引力快速确定方法具有以下优点：计算量不大，不需要特别高级的计算机即可实现，无须较多种类的重力测量数据，逼近精度高，设计扰动引力时，结合高速飞行的飞行器的应用需求，使得扰动引力的计算可高效、高精度进行。

附图说明

图1是本发明空间分层网格扰动引力场模型构建与扰动引力快速确定方法中空间分层扰动引力场网格划分准则形成过程图；

图2是本发明空间分层网格扰动引力场模型构建与扰动引力快速确定方法中扰动引力场网格单元示意图；

图3是本发明空间分层网格扰动引力场模型构建与扰动引力快速确定方法中分层网格模型构建流程图；

图4是本发明空间分层网格扰动引力场模型构建与扰动引力快速确定方法中分层网格模型数据输出格式图；

图5是本发明空间分层网格扰动引力场模型构建与扰动引力快速确定方法中扰动引力快速解算流程图。

具体实施方式

下面结合附图和具体实施方式对本发明空间分层网格扰动引力场模型构建与扰动引力快速确定方法作进一步说明：含有以下步骤。

一种空间分层网格扰动引力场模型构建，由以下步骤实现：一、高度层面划分，即根据全球外空扰动引力变化特征，将外部空间划分成若干个区间；二、网格单元逼近与网格划分准则的制定，即在不同区间，用基于广域多项式的网格单元分别逼近，并分析精度，将逼近精度较高且数据生成量较小的方案列出，生成扰动引力分层网格模型构建的准则。

一种包含所述空间分层网格扰动引力场模型构建的扰动引力快速确定方法，包含以下步骤：(1)利用扰动引力分层网格模型构建的准则进行扰动引力模型构建，即设计模型构建流程与数据输出格式；(2)利用构建好的扰动引力模型进行扰动引力实时解算，即设计扰动引力计算流程，结合高速飞行的飞行器的应用需求，使得扰动引力的计算可高效、高精度进行。

下面结合附图对本发明进行详细说明:

一、高度层面划分

根据图1中步骤S1，分析扰动引力传播特性。分析过程中，利用全球均方值来估计扰动引力分量的谱能量，公式见下。

其中，C_n为重力异常阶方差，表示如下，

$C_n={\mathrm{\γ}}^2{(n-1)}^2\underset{m=0}{\overset{n}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{nm}}^2+{\stackrel{\‾}{S}}_{\mathrm{nm}}^2)$

为正则化后的模型位系数，此处使用的是EGM2008位模型。

由于扰动引力径向和水平分量谱形式基本相同，故仅将径向结果(δr部分)列在下表。其中F1表示2～36阶频段，F2表示37～180阶频段，F3表示181～540阶频段，F4表示541～720阶频段，F5表示721～1080阶频段，F6表示1081～2160阶频段，F7表示2～18阶频段，F8表示19～36阶频段；Ei、Di(i＝1,2,…,8)分别表示某一高度层不同频段的能量和所占百分比。

表1   扰动引力径向分量在不同高度不同频段能量及其所占百分比

从上表可以看出随着高度的增加，重力场高频部分能量比迅速减小，如200km以上的层面，仅剩下180阶以下的能量。

根据图1步骤S1，需统计不同层面扰动引力梯度。在这里仅以扰动引力径向分量的径向梯度为例进行统计分析，如下表。

表2   20～2000km扰动引力径向分量径向梯度

高度(km)	最小值(E)	最大值(E)	平均值(E)	标准差(E)
					20	-42.70	96.00	-0.01	3.35
30	-32.00	59.30	-0.01	2.44
					40	-25.20	41.70	-0.01	1.89
50	-20.40	31.40	-0.01	1.52
					60	-16.70	23.90	-0.01	1.26
70	-13.80	18.70	-0.01	1.06
					80	-11.50	14.60	-0.01	0.91
90	-9.70	11.70	-0.01	0.80
					100	-8.18	9.38	-0.01	0.70
110	-7.06	7.62	-0.01	0.63
					120	-6.10	6.24	-0.01	0.56

130	-5.28	5.16	-0.01	0.51
					140	-4.57	4.31	-0.01	0.47
150	-3.99	3.63	-0.01	0.43
					160	-3.52	3.15	-0.01	0.40
170	-3.13	2.84	-0.01	0.37
					180	-2.79	2.56	-0.01	0.35
190	-2.51	2.35	-0.01	0.33
					200	-2.26	2.19	-0.01	0.31
210	-2.05	2.06	-0.01	0.29
					220	-1.87	1.93	-0.01	0.28
250	-1.46	1.12	0.00	0.23
					300	-1.02	0.83	0.00	0.18
400	-0.66	0.55	0.00	0.15
					500	-0.47	0.47	0.00	0.12
600	-0.35	0.40	0.00	0.11
					800	-0.22	0.31	0.00	0.08
1000	-0.17	0.24	0.00	0.07
					1500	-0.10	0.15	0.00	0.04
2000	-0.06	0.09	0.00	0.02

由上表可以看出：

(1)高度越低，扰动引力高频信息越丰富，径向梯度频谱越复杂；

(2)随高度增加，扰动引力径向梯度绝对值最大值减小，全球梯度趋于统一；

(3)20～100km高度内，径向梯度最大绝对值从96E骤减至不足10E，在100km以上时，绝对值均小于10E，220km以上时小于2E，1000km以上时小于0.25E。

因此将外部空间按以下方式划分成4个区间：

表3   高度划分区间

区间编号	高度范围/km
		1	20～100
2	100～220
		3	220～1000
4	>1000

二、网格单元逼近与网格层面划分准则的建立。

利用如图2中所示的网格单元对空间进行剖分：将某块研究区域以一定的经纬度和高度间隔划分成若干个单元，形成一系列的矩形单元，称之为主域；再沿着矩形单元的每条边向外延伸，形成次域单元；其节点分别称之为主节点(8个)、次节点(24个)。

采用多项式的形式来逼近待求点的扰动引力，如下：

g_i(x_i,y_i,z_i)＝f_ia

其中

f_i＝(1,x_i,y_i,z_i,x_iy_i,x_iz_i,y_iz_i,...)

a＝(a₀,a₁,a₂...)^T

在每个节点处可得到一个多项式的值g_i(x,y,z)，而要达到的结果是主节点上的多项式值与扰动引力值相等，而次节点上的多项式的值与扰动引力的差值的平方和最小，即可建立方程，如下：

$\left\{\begin{array}{c}\min M\left(a\right)={(f^{\′}a-s^{\′})}^T(f^{\′}a-s^{\′})\\ \mathrm{fa}-s=0\end{array}\right.$

其中，fa为主节点上的多项式值，s为主节点上的扰动引力值，f′a为次节点上的多项式值，s′为次节点上的扰动引力值。根据最小二乘平差理论，可求得位置参数向量的值。其中节点(包括主节点和次节点)与待求点的坐标要转换到以单元中心为原点的局部坐标系下，转换公式如下：

其中，是主节点的球坐标，是立方体中心点的球坐标。

根据图1中步骤S3、S4，以扰动引力径向分量为例，应用不同大小的网格单元、不同个数系数参数的多项式对不同高度的地区进行逼近，并分析精度。

根据多项式逼近原理，多项式的项数n应该满足P1<n<P1+P2(其中P1为主域节点数，P2为次域节点数)。一般选用的单元为空间六面体(上下面为弧面)结构，所以P1＝8，P2＝24。同时考虑到多项式中三个坐标参数的轮换对称性，次数只宜取到2次或3次，故多项式向量就只有以下4个，相应的多项式系数参数的个数分别为10、14、17、20。

表4

下面，根据上个步骤中对高度层面的划分结果([20,100km]、[100,220km]、[220,1000km]、>1000km)，用不同参数进行逼近，结果如下：

(1)在20～100km高度范围，当梯度绝对值小于20E时，使用10’×10’×20km格网进行逼近，参数个数选为20，精度可优于1mGal；

(2)在20～100km高度范围，当梯度绝对值大于等于20E时，使用5’×5’×10km格网进行逼近，参数个数选为20，精度可优于1mGal；

(3)在100～220km高度范围，当梯度绝对值小于2.5E时，使用60’×60’×120km格网时，20参数的逼近精度可优于1mGal；

(4)在100～220km高度范围，当扰动引力梯度绝对值属于区间[2.5,10E]时，使用20’×20’×40km格网，且参数选为14，其逼近精度即可优于1mGal；

(5)在220km～1000高度范围，当梯度绝对值小于0.5E时，使用2°×2°×240km格网，且参数选为17，逼近精度即可优于1mGal；

(6)在220km～1000高度范围，当梯度绝对值属于区间[0.5,1E]时，使用2°×2°×240km格网，且参数选为20，逼近精度即可优于1mGal；

(7)在220km～1000高度范围，当梯度绝对值属于区间[1,2E]时，使用60’×60’×120km格网，且参数选为14，逼近精度即可优于1mGal；

(8)在1000km以上高度范围，当扰动引力梯度绝对值小于0.1E时，使用5°×5°×600km格网，且参数选为10，逼近精度可优于1mGal；

(9)在1000km以上高度范围，当扰动引力梯度绝对值属于区间[0.1,0.25E]时，使用5°×5°×600km格网，且参数选为14，逼近精度可优于1mGal。

这便使得在精度优于1mGal的前提下，逼近效率达到较高水平，同时数据存储量较少。以下为表格形式的准则。

表5  分层及格网大小选取准则

三、利用准则进行扰动引力模型的构建。

根据上表可以进行扰动引力场分层网格模型构建的流程，如图3。首先计算扰动引力场节点值，利用准则自适应地选取所需单元的节点并计算单元内的广域多项式系数；将单元合理编号，并将计算得到的系数以一定的次序存入文件中对应的单元编号后；最终就得到了分层网格模型。

其中高度区间稍作调整，使区间可以被网格完整剖分。当高度大于2380km时，可直接用18阶位系数模型计算，即可满足要求。模型输出形式为图4。将生成的模型数据提前存入飞行器载体的计算机硬盘中。

四、利用模型快速确定扰动引力。

结合飞行器特点，对扰动引力实时解算过程进行说明。

(1)通过飞行器载体上的GPS接收机、惯性导航及其他的地形或重力匹配导航系统，精确确定飞行器位置坐标信息；

(2)对坐标信息进行分析，计算得到飞行器所在位置的层号及单元号。利用载体计算机读取生成的模型数据文件，搜索对应层号和单元号的多项式系数值。

(3)将飞行器坐标值代入多项式中进行扰动引力解算，其中多项式系数已在上一步中确定。由于解算过程仅为多项式求和运算，所以单点的解算速度可保持在10^-5秒以内，非常适合用于高速飞行的飞行器。

具体流程参见图5。

由此方法构建的模型，可使飞行器的扰动引力计算满足快速实时且高精度的要求，在扰动引力快速计算领域是一大创新，具有较强的实际应用价值。

一、具体算例。

选取3个不同位置处的点进行计算，点的位置分别为：(10°47′30″N，73°27′30″W，高度29km)、(11°42′30″N，59°37′30″W，高度238km)、(0°17′30″S，126°32′30″E，高度119km)，下面将三个点分别记为P1，P2，P3。

首先，通过高度判断层面区间，并结合划分准则(表5)，将对应高度区间的梯度区间列在下表。

表6

点位	高度区间	梯度区间
			P1(高29km)	20-100km	0～20E,>20E
P2(高238km)	220-1180km	0～0.5E,0.5～1E,>1E
			P3(高119km)	100-220km	0～2.5E,>2.5E

其次，通过不同点(P1，P2，P3)所在位置的梯度值判断所属梯度区间，再结合划分准则(表5)，获取对应单元的网格分辨率和参数个数，如下表。

表7

点位	所属梯度区间	网格分辨率	参数个数
				P1	>20E	5′	20
P2	>1E	20′	14
				P3	0～2.5E	20′	14

再次，通过解算相应格网分辨率和参数个数的广域多项式，得到不同点(P1，P2，P3)所在单元的多项式系数如下。

表8

最后，将上表中的系数分别代入表4中的多项式向量4和多项式向量2中，并代入点位(P1，P2，P3)的坐标值，即可得到相应点的扰动引力值。

通过计算，最终解得扰动引力值分别为：-227.493mGal、35.409mGal、-0.853mGal，耗时为4.4×10^-6秒。将EGM2008模型计算的值作为真值，可以发现三个点的扰动引力计算误差分别为：0.346mGal、0.177mGal、0.046mGal，均小于1mGal。

故利用这种模型可快速获取任意点的扰动引力，且精度可优于1mGal。

Claims (7)

1.一种空间分层网格扰动引力场模型构建，其特征在于：由以下步骤实现：

一、高度层面划分，即根据全球外空扰动引力变化特征，将外部空间划分成若干个区间；

二、网格单元逼近与网格划分准则的制定，即在不同区间，用基于广域多项式的网格单元分别逼近，并分析精度，将逼近精度较高且数据生成量较小的方案列出，生成扰动引力分层网格模型构建的准则。

2.根据权利要求1所述的空间分层网格扰动引力场模型构建，其特征在于：所述高度层面划分的分析过程中，利用全球均方值来估计扰动引力分量的谱能量，公式为

其中，C_n为重力异常阶方差，表示如下，

$C_n={\mathrm{\&gamma;}}^2{(n-1)}^2\underset{m=0}{\overset{n}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{\mathrm{nm}}^2+{\stackrel{\&OverBar;}{S}}_{\mathrm{nm}}^2)$

为正则化后的模型位系数，此处使用的是EGM2008位模型。

3.根据权利要求1所述的空间分层网格扰动引力场模型构建，其特征在于：所述区间为4个，按高度范围划分即20-100KM、100-220KM、220-1000KM和1000KM以上。

4.根据权利要求1所述的空间分层网格扰动引力场模型构建，其特征在于：所述网格单元逼近与网格层面划分准则的建立过程中，将某块研究区域以一定的经纬度和高度间隔划分成若干个单元，形成一系列的矩形单元，称之为主域；再沿着矩形单元的每条边向外延伸，形成次域单元；其节点分别称之为主节点(8个)、次节点(24个)；采用多项式的形式来逼近待求点的扰动引力，如下：g_i(x_i,y_i,z_i)＝f_ia

其中

f_i＝(1,x_i,y_i,z_i,x_iy_i,x_iz_i,y_iz_i,...)

a＝(a₀,a₁,a₂...)^T

在每个节点处可得到一个多项式的值g_i(x,y,z)，而要达到的结果是主节点上的多项式值与扰动引力值相等，而次节点上的多项式的值与扰动引力的差值的平方和最小，即可建立方程，如下：

$\left\{\begin{array}{c}\min M\left(a\right)={(f^{\&prime;}a-s^{\&prime;})}^T(f^{\&prime;}a-s^{\&prime;})\\ \mathrm{fa}-s=0\end{array}\right.$

其中，fa为主节点上的多项式值，s为主节点上的扰动引力值，f′a为次节点上的多项式值，s′为次节点上的扰动引力值；根据最小二乘平差理论，可求得位置参数向量的值；其中节点(包括主节点和次节点)与待求点的坐标要转换到以单元中心为原点的局部坐标系下，转换公式如下：

其中，是主节点的球坐标，是立方体中心点的球坐标。

5.一种包含权利要求1所述的空间分层网格扰动引力场模型构建的扰动引力快速确定方法，其特征在于：包含以下步骤：

(1)利用扰动引力分层网格模型构建的准则进行扰动引力模型构建，即设计模型构建流程与数据输出格式；

(2)利用构建好的扰动引力模型进行扰动引力实时解算，即设计扰动引力计算的流程。

6.根据权利要求5所述的扰动引力快速确定方法，其特征在于：所述利用扰动引力分层网格模型构建的准则进行扰动引力模型的构建过程如下：首先计算扰动引力场节点值，利用准则自适应地选取所需单元的节点并计算单元内的广域多项式系数；将单元合理编号，并将计算得到的系数以一定的次序存入文件中对应的单元编号后；最终就得到了分层网格模型。

7.根据权利要求5所述的扰动引力快速确定方法，其特征在于：所述利用构建好的扰动引力模型进行扰动引力实时解算包括以下步骤：(1)通过飞行器载体上的GPS接收机、惯性导航及其他的地形或重力匹配导航系统，精确确定飞行器位置坐标信息；(2)对坐标信息进行分析，计算得到飞行器所在位置的层号及单元号；利用载体计算机读取生成的模型数据文件，搜索对应层号和单元号的多项式系数值；(3)将飞行器坐标值代入对应多项式中进行扰动引力解算。