CN104834320A - 一种空间分层扰动引力场网格模型快速构建方法 - Google Patents
一种空间分层扰动引力场网格模型快速构建方法 Download PDFInfo
- Publication number
- CN104834320A CN104834320A CN201510126443.XA CN201510126443A CN104834320A CN 104834320 A CN104834320 A CN 104834320A CN 201510126443 A CN201510126443 A CN 201510126443A CN 104834320 A CN104834320 A CN 104834320A
- Authority
- CN
- China
- Prior art keywords
- disturbance
- formula
- model
- low
- gravitation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000010276 construction Methods 0.000 title abstract description 6
- 230000005484 gravity Effects 0.000 claims abstract description 42
- 238000004364 calculation method Methods 0.000 claims abstract description 30
- 230000010354 integration Effects 0.000 claims abstract description 23
- 238000000034 method Methods 0.000 claims abstract description 22
- 230000009466 transformation Effects 0.000 abstract description 2
- 238000007796 conventional method Methods 0.000 abstract 1
- 238000011084 recovery Methods 0.000 description 14
- 238000005516 engineering process Methods 0.000 description 8
- 238000010586 diagram Methods 0.000 description 4
- 230000009286 beneficial effect Effects 0.000 description 3
- 230000002159 abnormal effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 238000002806 Stokes method Methods 0.000 description 1
- 230000005856 abnormality Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000004248 saffron Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Landscapes
- Complex Calculations (AREA)
Abstract
本发明涉及空间分层扰动引力场网格模型快速构建方法,可有效解决用于快速、高精度地得到地球外部任意高度层面的扰动引力问题,包括以下步骤:传统积分公式改化为卷积形式,进而利用卷积定理进行快速计算;低阶模型重力异常与剩余重力异常计算,得到剩余重力异常;傅里叶变换及其逆变换,求得剩余扰动引力值;低阶模型扰动引力与最终扰动引力计算,得最终扰动引力;扰动引力精度与计算效率比较,同时与传统方法比较效率差别,本发明方法易操作应用,可有效解决快速、高精度地得到地球外部任意高度层面的扰动引力,有助于空间各种飞行器在运动时的快速实时定位和控制,保证飞行器的安全,具有很强的使用价值。
Description
技术领域
本发明涉及航空航天技术,特别是一种空间分层扰动引力场网格模型快速构建方法。
背景技术
随着航空航天技术的发展,各种飞行器的安全已为人们所关注,而影响飞行器安全的一个重要原因在于扰动引力作为地球重力场快速高精度的计算,确保各种飞行器在运动时的快速实时定位,目前计算扰动引力是根据Stokes理论,在视地球为球体且忽略零阶项的情况下,地球外部点的扰动位与地面重力异常有如下的关系式,
式(1)
式中,Δg为球面上积分流动单元dΣ处的重力异常,S(ρ,ψ)为Stokes函数,其表达式为,
ρ为待求点与地球质心的距离,r为待求点至dΣ的距离,ψ为dΣ的向径与ρ的夹角,R为地球平均半径,Σ为积分区域,这里表示全球范围。
扰动引力是扰动位的偏导数,分别对向径ρ、地心纬度地心经度λ求偏导数,即可得扰动引力在相应方向的分量。经过化简,有
式(3)
其中和的表达式可以由下式求得
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∑的地心纬度和经度;P0为待求点P在球面的投影点,其球面坐标为ψ为P0与积分元d∑之间的球心角;A为待求点与dΣ连线的方位角;P0与球面元素d∑,以及北极点Np构成一个球面三角形,可给出以下几个关系式:
式(5)
式(6)
式(7)
对式(5)两边分别求λ的偏导数,并联立式(6)、式(7),得:
式(8)
进一步的,当λ'=λ-Δλ,则有:
式(9)
式(10)
将式(9)代入式(5),有
式(11)
将式(8)、(9)、(10)、(11)与式(3)联立,得:
式(12)
其中,
则式(12)即为卷积形式的Stokes积分公式,当积分半径确定之后,的值即确定,将三角函数值计算并存入动态数组,避免反复计算,提高计算效率;
二、低阶模型重力异常计算:
计算公式为:
式(13)
式(13)中ΔgM为低阶模型重力异常,γ为地球平均重力,N为模型截断阶数,为扣除正常项的地球重力位模型系数,λ、为对应点地心经纬度值,为正则化的缔合勒让德函数,n、m分别为阶、次数;
三、剩余重力异常计算:
剩余重力异常即为式(12)中的δΔgo,计算公式为:δΔgo=Δg-ΔgM,ΔgM在式(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∑的地心纬度和经度;P0为待求点P在球面的投影点,其球面坐标为ψ为P0与积分元d∑之间的球心角;A为待求点与dΣ连线的方位角;P0与球面元素d∑,以及北极点Np构成一个球面三角形,可给出以下几个关系式:
对式(5)两边分别求λ的偏导数,并联立式(6)、式(7),得:
进一步的,当λ'=λ-Δλ,则有:
将式(9)代入式(5),有
将式(8)、(9)、(10)、(11)与式(3)联立,得:
其中,
则式(12)即为卷积形式的Stokes积分公式,当积分半径确定之后,的值即确定,将三角函数值计算并存入动态数组,避免反复计算,提高计算效率;
二、低阶模型重力异常计算:
计算公式为:
式(13)中ΔgM为低阶模型重力异常,γ为地球平均重力,N为模型截断阶数,为扣除正常项的地球重力位模型系数,λ、为对应点地心经纬度值,为正则化的缔合勒让德函数,n、m分别为阶、次数;
三、剩余重力异常计算:
剩余重力异常即为式(12)中的δΔgo,计算公式为:δΔgo=Δg-ΔgM,ΔgM在式(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λ
其中δρ、δλ分别为最终所求的扰动引力矢量的各个分量值。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510126443.XA CN104834320A (zh) | 2015-03-23 | 2015-03-23 | 一种空间分层扰动引力场网格模型快速构建方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510126443.XA CN104834320A (zh) | 2015-03-23 | 2015-03-23 | 一种空间分层扰动引力场网格模型快速构建方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN104834320A true CN104834320A (zh) | 2015-08-12 |
Family
ID=53812264
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510126443.XA Pending CN104834320A (zh) | 2015-03-23 | 2015-03-23 | 一种空间分层扰动引力场网格模型快速构建方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104834320A (zh) |
Cited By (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108376187A (zh) * | 2018-01-19 | 2018-08-07 | 中国人民解放军92859部队 | 一种海域流动点外部扰动引力垂向分量的无奇异计算方法 |
CN108846195A (zh) * | 2018-06-11 | 2018-11-20 | 中国人民解放军61540部队 | 一种极区扰动重力无奇异性详细计算模型及其建模方法 |
CN112818285A (zh) * | 2021-02-08 | 2021-05-18 | 中国人民解放军92859部队 | 一种计算外部扰动重力北向分量中央区效应的方法 |
CN112949049A (zh) * | 2021-02-08 | 2021-06-11 | 中国人民解放军92859部队 | 一种利用带限思想计算重力异常低阶径向导数的方法 |
CN112965128A (zh) * | 2021-02-08 | 2021-06-15 | 中国人民解放军92859部队 | 一种无奇异性顾及局域保障条件计算外部重力异常的方法 |
CN112965124A (zh) * | 2021-02-08 | 2021-06-15 | 中国人民解放军92859部队 | 一种顾及局域保障条件计算外部重力异常垂直梯度的方法 |
CN112965125A (zh) * | 2021-02-08 | 2021-06-15 | 中国人民解放军92859部队 | 一种基于重力异常计算外部扰动重力东向分量的方法 |
-
2015
- 2015-03-23 CN CN201510126443.XA patent/CN104834320A/zh active Pending
Cited By (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108376187A (zh) * | 2018-01-19 | 2018-08-07 | 中国人民解放军92859部队 | 一种海域流动点外部扰动引力垂向分量的无奇异计算方法 |
CN108846195A (zh) * | 2018-06-11 | 2018-11-20 | 中国人民解放军61540部队 | 一种极区扰动重力无奇异性详细计算模型及其建模方法 |
CN112818285A (zh) * | 2021-02-08 | 2021-05-18 | 中国人民解放军92859部队 | 一种计算外部扰动重力北向分量中央区效应的方法 |
CN112949049A (zh) * | 2021-02-08 | 2021-06-11 | 中国人民解放军92859部队 | 一种利用带限思想计算重力异常低阶径向导数的方法 |
CN112965128A (zh) * | 2021-02-08 | 2021-06-15 | 中国人民解放军92859部队 | 一种无奇异性顾及局域保障条件计算外部重力异常的方法 |
CN112965124A (zh) * | 2021-02-08 | 2021-06-15 | 中国人民解放军92859部队 | 一种顾及局域保障条件计算外部重力异常垂直梯度的方法 |
CN112965125A (zh) * | 2021-02-08 | 2021-06-15 | 中国人民解放军92859部队 | 一种基于重力异常计算外部扰动重力东向分量的方法 |
CN112949049B (zh) * | 2021-02-08 | 2021-11-30 | 中国人民解放军92859部队 | 一种利用带限思想计算重力异常低阶径向导数的方法 |
CN112965125B (zh) * | 2021-02-08 | 2022-08-05 | 中国人民解放军92859部队 | 一种基于重力异常计算外部扰动重力东向分量的方法 |
CN112818285B (zh) * | 2021-02-08 | 2022-09-30 | 中国人民解放军92859部队 | 一种计算外部扰动重力北向分量中央区效应的方法 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104834320A (zh) | 一种空间分层扰动引力场网格模型快速构建方法 | |
CN107219857B (zh) | 一种基于三维全局人工势函数的无人机编队路径规划算法 | |
Li et al. | Rainfall and earthquake-induced landslide susceptibility assessment using GIS and Artificial Neural Network | |
CN103076017B (zh) | 基于可观测度分析的火星进入段自主导航方案设计方法 | |
CN103743402B (zh) | 一种基于地形信息量的水下智能自适应地形匹配方法 | |
CN102707726B (zh) | 一种无人机目标定位方法 | |
CN105509750B (zh) | 一种天文测速与地面无线电组合的火星捕获段导航方法 | |
CN104834017A (zh) | 一种有毒有害气体泄漏扩散事故源定位方法 | |
CN105279581A (zh) | 基于差分进化的geo-uav双基sar路径规划方法 | |
CN107704641A (zh) | 基于实景植被空间分布粗糙度的精细风场模拟方法 | |
CN104751012A (zh) | 沿飞行弹道的扰动引力快速逼近方法 | |
CN104266650B (zh) | 一种基于采样点继承策略的火星着陆器大气进入段导航方法 | |
CN114936471A (zh) | 一种基于并行计算的航天器碰撞预警分层快速筛选方法 | |
CN104596636A (zh) | 声场分离方法 | |
CN104504255A (zh) | 一种螺旋翼升力和阻力力矩的确定方法 | |
CN103616024B (zh) | 一种行星探测进入段自主导航系统可观测度确定方法 | |
CN103344245A (zh) | 火星进入段imu和甚高频无线电组合导航的ud-skf方法 | |
Zhang et al. | Surface area processing in GIS for different mountain regions | |
CN110231619B (zh) | 基于恩克法的雷达交接时刻预报方法及装置 | |
CN112896560A (zh) | 小天体表面安全弹跳移动轨迹规划方法 | |
Wang et al. | Real-Time data driven simulation of air contaminant dispersion using particle filter and UAV sensory system | |
CN104408326B (zh) | 一种对深空探测自主导航滤波算法的评估方法 | |
CN102890743A (zh) | 行星大气进入着陆器落点不确定度分析方法 | |
Leukauf et al. | The impact of a forest parametrization on coupled WRF-CFD simulations during the passage of a cold front over the WINSENT test-site | |
CN105424047A (zh) | 基于路标信息的航天器姿态指向误差辨识方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
EXSB | Decision made by sipo to initiate substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
WD01 | Invention patent application deemed withdrawn after publication | ||
WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20150812 |