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

Abstract

本发明涉及空间分层扰动引力场网格模型快速构建方法，可有效解决用于快速、高精度地得到地球外部任意高度层面的扰动引力问题，包括以下步骤：传统积分公式改化为卷积形式，进而利用卷积定理进行快速计算；低阶模型重力异常与剩余重力异常计算，得到剩余重力异常；傅里叶变换及其逆变换，求得剩余扰动引力值；低阶模型扰动引力与最终扰动引力计算，得最终扰动引力；扰动引力精度与计算效率比较，同时与传统方法比较效率差别，本发明方法易操作应用，可有效解决快速、高精度地得到地球外部任意高度层面的扰动引力，有助于空间各种飞行器在运动时的快速实时定位和控制，保证飞行器的安全，具有很强的使用价值。

Description

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

技术领域

本发明涉及航空航天技术，特别是一种空间分层扰动引力场网格模型快速构建方法。

背景技术

随着航空航天技术的发展，各种飞行器的安全已为人们所关注，而影响飞行器安全的一个重要原因在于扰动引力作为地球重力场快速高精度的计算，确保各种飞行器在运动时的快速实时定位，目前计算扰动引力是根据Stokes理论，在视地球为球体且忽略零阶项的情况下，地球外部点的扰动位与地面重力异常有如下的关系式，

   式(1)

式中，Δg为球面上积分流动单元dΣ处的重力异常，S(ρ,ψ)为Stokes函数，其表达式为，

$S(\mathrm{\ρ},\mathrm{\ψ})=\frac{2}{r}+\frac{1}{\mathrm{\ρ}}-\frac{3r}{{\mathrm{\ρ}}^2}-\frac{5R}{{\mathrm{\ρ}}^2}\cos \mathrm{\ψ}-\frac{3R}{{\mathrm{\ρ}}^2}\cos \mathrm{\ψ}\ln \frac{\mathrm{\ρ}-R\cos \mathrm{\ψ}+r}{2\mathrm{\ρ}}$    式(2)

ρ为待求点与地球质心的距离，r为待求点至dΣ的距离，ψ为dΣ的向径与ρ的夹角，R为地球平均半径，Σ为积分区域，这里表示全球范围。

扰动引力是扰动位的偏导数，分别对向径ρ、地心纬度地心经度λ求偏导数，即可得扰动引力在相应方向的分量。经过化简，有

   式(3)

其中和的表达式可以由下式求得

$\left\{\begin{array}{c}\frac{\∂S(\mathrm{\ρ},\mathrm{\ψ})}{\∂\mathrm{\ρ}}=-\frac{2(\mathrm{\ρ}-R\cos \mathrm{\ψ})}{r^3}-\frac{3}{\mathrm{\ρr}}-\frac{1}{{\mathrm{\ρ}}^2}+\frac{6r}{{\mathrm{\ρ}}^3}+\frac{13R}{{\mathrm{\ρ}}^3}\cos \mathrm{\ψ}+\frac{6R}{{\mathrm{\ρ}}^3}\cos \mathrm{\ψ}\ln \frac{\mathrm{\ρ}-R\cos \mathrm{\ψ}+r}{2\mathrm{\ρ}}\\ \frac{\∂S(\mathrm{\ρ},\mathrm{\ψ})}{\∂\mathrm{\ψ}}=\sin \mathrm{\ψ}[-\frac{2\mathrm{\ρ}}{r^3}-\frac{3}{\mathrm{\ρr}}+\frac{5}{{\mathrm{\ρ}}^2}+\frac{5}{{\mathrm{\ρ}}^2}-\frac{3R(r+\mathrm{\ρ})}{{\mathrm{\ρ}}^{2}r(r+\mathrm{\ρ}-R\cos \mathrm{\ψ})}\cos \mathrm{\ψ}+\frac{3}{{\mathrm{\ρ}}^3}\ln \frac{\mathrm{\ρ}-R\cos \mathrm{\ψ}+r}{2\mathrm{\ρ}}]R\end{array}\right.$    式(4)

A为待求点与dΣ连线的方位角。由此，便可由地面重力异常构建空间扰动引力矢量模型。

然而，由于Stokes积分核较复杂且庞大，且每次计算一个点需要对全球的重力异常进行积分，这使该方法计算扰动引力的效率异常低下，严重限制了该方法的应用。因而可利用快速傅里叶变换技术替代积分计算，并利用移去-恢复技术代替积分远区的计算提高计算效率。

前人曾利用Stokes积分法建立过扰动位、大地水准面高等模型，但还未有人利用该方法建立过空间扰动引力模型，更未有人将快速傅里叶变换等技术引入该方法进行快速构建。而扰动引力作为地球重力场领域中的研究重点，深入探讨其快速、高精度的计算方法，有助于空间各种飞行器在运动时的快速实时定位和控制。但至今还未见有效用于快速、准确计算扰动引力的公开报导。

发明内容

针对上述情况，为克服现有技术之缺陷，本发明之目的就是提供一种空间分层扰动引力场网格模型快速构建方法，可有效解决用于快速、高精度地得到地球外部任意高度层面的扰动引力问题。

本发明解决的技术方案是，包括以下步骤：

一、传统积分公式改化为卷积形式：传统的扰动引力计算公式效率较低，为利用FFT技术加快计算效率，需将其变换为卷积形式，进而利用卷积定理进行快速计算；

二、低阶模型重力异常与剩余重力异常计算：利用位系数模型计算低阶模型重力异常值，在原始的重力异常值中扣除该量，得到剩余重力异常，即移去-恢复技术中的“移去”步骤；

三、傅里叶变换及其逆变换：当剩余重力异常计算完毕后，带人已改化为卷积形式的积分公式中，然后利用卷积定理，进行傅里叶的正、逆变换，求得剩余扰动引力值；

四、低阶模型扰动引力与最终扰动引力计算：利用傅里叶变换和逆变换求得剩余扰动引力之后，加上低阶模型扰动引力，即可得最终扰动引力，即移去-恢复技术中的“恢复”步骤；

五、扰动引力精度与计算效率比较：选取不同地形特征的地区，观察快速方法计算的模型的精度，同时与传统方法比较效率差别。

本发明方法易操作应用，可有效解决快速、高精度地得到地球外部任意高度层面的扰动引力，有助于空间各种飞行器在运动时的快速实时定位和控制，保证飞行器的安全，具有很强的使用价值。

附图说明

图1为本发明的空间扰动引力模型构建的技术路线图。

图2为本发明的Stokes积分点位示意图。

图3为本发明的全球地面高程图及地区选取方案。

图4为本发明的低阶模型阶数为360时不同积分半径在不同高度的恢复程度图。

图5为本发明的低阶模型阶数为180时不同积分半径在不同高度的恢复程度图。

图6为本发明的低阶模型阶数为36时不同积分半径在不同高度的恢复程度图。

图7为本发明的低阶模型阶数为18时不同积分半径在不同高度的恢复程度图。

图8为本发明的低阶模型阶数为360时不同积分半径在不同高度的均方根误差图。

图9为本发明的低阶模型阶数为180时不同积分半径在不同高度的均方根误差图。

图10为本发明的低阶模型阶数为36时不同积分半径在不同高度的均方根误差图。

图11为本发明的低阶模型阶数为18时不同积分半径在不同高度的均方根误差图。

具体实施方式

以下结合具体情况对本发明的具体实施方式作详细说明。

本发明在具体实施中，包括以下步骤：

一、传统积分的卷积形式改化：

利用Stokes积分计算扰动引力的传统公式(3)，将式(3)改化(变换)为卷积形式，方法是，先构建一个球面三角形(如图2所示)，ρ、λ分别为待求点球坐标中的球心向径、球心纬度及经度，λ'分别为积分元d∑的地心纬度和经度；P₀为待求点P在球面的投影点，其球面坐标为ψ为P₀与积分元d∑之间的球心角；A为待求点与dΣ连线的方位角；P₀与球面元素d∑，以及北极点N_p构成一个球面三角形，可给出以下几个关系式：

   式(5)

   式(6)

   式(7)

对式(5)两边分别求λ的偏导数，并联立式(6)、式(7)，得：

   式(8)

进一步的，当λ'＝λ-Δλ，则有：

   式(9)

   式(10)

将式(9)代入式(5)，有

   式(11)

将式(8)、(9)、(10)、(11)与式(3)联立，得：

   式(12)

其中，

则式(12)即为卷积形式的Stokes积分公式，当积分半径确定之后，的值即确定，将三角函数值计算并存入动态数组，避免反复计算，提高计算效率；

二、低阶模型重力异常计算：

计算公式为：

  式(13)

式(13)中Δg_M为低阶模型重力异常，γ为地球平均重力，N为模型截断阶数，为扣除正常项的地球重力位模型系数，λ、为对应点地心经纬度值，为正则化的缔合勒让德函数，n、m分别为阶、次数；

三、剩余重力异常计算：

剩余重力异常即为式(12)中的δΔg_o，计算公式为：δΔg_o＝Δg-Δg_M，Δg_M在式(13)中已给出；

四、傅里叶变换及其逆变换：

利用卷积定理进行快速傅里叶变换及其逆变换，卷积定理是：设FFT[f(t)]＝F(s)，FFT[g(t)]＝G(s)，则：

f(t)*g(t)＝IFFT{FFT[f(t)]□FFT[g(t)]}

其中FFT表示傅里叶正变换算子，IFFT表示傅里叶逆变换算子，“*”表示卷积运算符，将此定理用于Stokes积分，即式(12)，得剩余扰动引力，如下，

   式(14)

五、低阶模型扰动引力计算：

利用位系数模型计算低阶模型扰动引力值，计算公式是：

   式(15)

其中R为地球平均半径，f为万有引力常数，M为包括大气在内的地球总质量，其余定义与低阶模型重力异常公式中相同(见上述，不再重复)；

六、最终扰动引力计算：

将步骤四得到的剩余扰动引力与步骤五得到低阶模型扰动引力求和，得最终的扰动引力结果，如下式：

δ_ρ＝ΔT_ρ+δ_Mρ

δ_λ＝ΔT_λ+δ_Mλ

其中δ_ρ、δ_λ分别为最终所求的扰动引力矢量的各个分量值。

本发明方法易操作，实用性强，速度快，效率高，精度高，有效用于空间扰动引力的计算，有助于空间各种飞行器在运动时的快速实时定位和控制，并经实地试验和扰动引力模型的精度与构建效率分析，取得了非常满意的有益技术效果，有关资料如下：

选取10km高度层分别用本发明的快速方法和传统方法计算全球区域扰动引力径向分量，在全球范围选取不同地形特征的三块区域研究快速算法的计算精度。选取区域范围如图3，图中区域①、②、③的经纬度区间分别为：(127°W～176°W，38°S～50°N)、(68°E～164°E，5°S～43°N)、(56°W～41°E，45°N～90°N)，依次表示中低纬度地形变化较平缓地区、中低纬度地形变化较剧烈地区及部分高纬度地区。

将精度统计如下表：

表1快速计算方法与传统Stokes积分计算结果差值比较

从表1可以看出，快速计算方法与传统Stokes积分方法的差值很小，在10^-6mGal以内，因此可认为两者精度相当。

将两种方法的计算耗时统计如下。

表2FFT与传统Stokes积分耗时比较

可以看出快速算法的计算效率为传统方法的数十倍，当积分半径较大时甚至可达数百倍。

表3计算设备硬件性能

另外，考虑到不同参数(包括积分半径和低阶模型的截断阶数)的Stokes方法所计算的模型，其模型精度与计算点的高度有一定关系。以下通过实验，统计这一关系，并制定出不同高度相应的参数选取准则，以进一步提高地球外部扰动引力场模型的构建效率。

将积分半径设置为0’、30’、1°、1°30’、2°、5°，其中0’表示用仅用低阶模型进行计算；将低阶位模型的阶数设置为360、180、36、18。并以截止2160阶的EGM2008模型计算结果为真值，将各种情况计算结果的模型精度随高度变化的情况统计如图3～10。

根据图3～10及表1的情况可得以下结论：

在低空区域随着积分半径的变大，逼近精度随之增加。这是因为低空地区扰动引力场高频信息丰富，积分半径的增大可以有效收集高频信息，使得精度提高；

随着高度的增加，逼近精度随之增加，但是积分半径大的结果的精度逐渐被积分半径小的结果超越。这是因为高空地区，扰动引力场高频信息衰减殆尽，积分半径的增加反而使累积误差增大，最终该误差逐渐掩盖其在低空地区的精度优势。

根据选用不同参数时，模型构建的精度与耗时情况，可制定以下参数选取原则：

(1)、当高度在20km以下时，低阶模型阶数设为360、积分半径为2°，得到的扰动引力模型中误差小于1.2mGal，恢复程度大于94％；

(2)、当高度在20～60km时，低阶模型阶数设为360、积分半径为1°，得到的扰动引力模型中误差小于0.1mGal，恢复程度大于98％；

(3)、当高度在60～100km时，低阶模型阶数设为360、积分半径为30’，得到的扰动引力模型中误差小于0.14mGal，恢复程度大于98％；

(4)、当高度在100～200km时，低阶模型阶数设为180、积分半径为1°，得到的扰动引力模型中误差小于0.16mGal，恢复程度大于97％；

(5)、当高度在200～400km时，低阶模型阶数设为180、积分半径为30’，得到的扰动引力模型中误差小于0.11mGal，恢复程度大于98％；

(6)、当高度在400～800km时，低阶模型阶数设为36、积分半径为30’，得到的扰动引力模型中误差小于0.13mGal，恢复程度大于96％；

(7)、当高度在800km以上时，低阶模型阶数设为18、积分半径为30’，得到的扰动引力模型中误差小于0.04mGal，恢复程度大于97％。

将选取方案制成表格如下。

表4参数选取原则表

由上述资料充分证明，利用本发明提出的方法，可以高效、高精度的计算全球范围不同层面扰动引力模型，为飞行器的飞行安全提供了有力的技术保障，具有很强的实际应用价值。

Claims (1)

1.一种空间分层扰动引力场网格模型快速构建方法，其特征在于，包括以下步骤：

一、传统积分的卷积形式改化：

利用Stokes积分计算扰动引力的传统公式(3)，将式(3)改化为卷积形式，方法是，先构建一个球面三角形，ρ、λ分别为待求点球坐标中的球心向径、球心纬度及经度，λ'分别为积分元d∑的地心纬度和经度；P₀为待求点P在球面的投影点，其球面坐标为ψ为P₀与积分元d∑之间的球心角；A为待求点与dΣ连线的方位角；P₀与球面元素d∑，以及北极点N_p构成一个球面三角形，可给出以下几个关系式：

对式(5)两边分别求λ的偏导数，并联立式(6)、式(7)，得：

进一步的，当λ'＝λ-Δλ，则有：

将式(9)代入式(5)，有

将式(8)、(9)、(10)、(11)与式(3)联立，得：

其中，

则式(12)即为卷积形式的Stokes积分公式，当积分半径确定之后，的值即确定，将三角函数值计算并存入动态数组，避免反复计算，提高计算效率；

二、低阶模型重力异常计算：

计算公式为：

式(13)中Δg_M为低阶模型重力异常，γ为地球平均重力，N为模型截断阶数，为扣除正常项的地球重力位模型系数，λ、为对应点地心经纬度值，为正则化的缔合勒让德函数，n、m分别为阶、次数；

三、剩余重力异常计算：

剩余重力异常即为式(12)中的δΔg_o，计算公式为：δΔg_o＝Δg-Δg_M，Δg_M在式(13)中已给出；

四、傅里叶变换及其逆变换：

利用卷积定理进行快速傅里叶变换及其逆变换，卷积定理是：设FFT[f(t)]＝F(s)，FFT[g(t)]＝G(s)，则：

f(t)*g(t)＝IFFT{FFT[f(t)]□FFT[g(t)]}

其中FFT表示傅里叶正变换算子，IFFT表示傅里叶逆变换算子，“*”表示卷积运算符，将此定理用于Stokes积分，即式(12)，得剩余扰动引力，如下，

五、低阶模型扰动引力计算：

利用位系数模型计算低阶模型扰动引力值，计算公式是：

其中R为地球平均半径，f为万有引力常数，M为包括大气在内的地球总质量，其余定义与低阶模型重力异常公式中相同；

六、最终扰动引力计算：

将步骤四得到的剩余扰动引力与步骤五得到低阶模型扰动引力求和，得最终的扰动引力结果，如下式：

δ_ρ＝ΔT_ρ+δ_Mρ

δ_λ＝ΔT_λ+δ_Mλ

其中δ_ρ、δ_λ分别为最终所求的扰动引力矢量的各个分量值。