CN112987118B - 一种利用带限思想计算重力异常高阶径向导数的方法 - Google Patents

一种利用带限思想计算重力异常高阶径向导数的方法 Download PDF

Info

Publication number
CN112987118B
CN112987118B CN202110180193.3A CN202110180193A CN112987118B CN 112987118 B CN112987118 B CN 112987118B CN 202110180193 A CN202110180193 A CN 202110180193A CN 112987118 B CN112987118 B CN 112987118B
Authority
CN
China
Prior art keywords
gravity
order radial
order
radial derivative
gravity anomaly
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.)
Active
Application number
CN202110180193.3A
Other languages
English (en)
Other versions
CN112987118A (zh
Inventor
邓凯亮
黄谟涛
陈欣
黄辰虎
吴太旗
陆秀平
熊雄
黄贤源
刘玉红
范瑾煜
董斌
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
92859 TROOPS PLA
Original Assignee
92859 TROOPS PLA
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by 92859 TROOPS PLA filed Critical 92859 TROOPS PLA
Priority to CN202110180193.3A priority Critical patent/CN112987118B/zh
Publication of CN112987118A publication Critical patent/CN112987118A/zh
Application granted granted Critical
Publication of CN112987118B publication Critical patent/CN112987118B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00Measuring gravitational fields or waves; Gravimetric prospecting or detecting
    • G01V7/02Details
    • G01V7/06Analysis or interpretation of gravimetric records

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及一种利用带限思想计算重力异常高阶径向导数的方法,其主要技术特点是:将地球外部重力异常Poisson积分式的解析核函数表示为球谐级数展开式,通过求导方法得到重力异常高阶径向导数球谐级数表达式;利用移去‑恢复技术将重力异常高阶径向导数球谐级数表达式截断为与重力观测值频谱范围相一致的带限求和式,并利用全球重力位模型补偿重力异常高阶径向导数远区截断误差,计算得到重力异常高阶径向导数。本发明设计合理,能够精确计算重力异常高阶径向导数,提高了重力异常高阶径向导数的计算精度,解决了重力异常高阶径向导数计算结果的不稳定性问题,可广泛用于物理大地测量领域。

Description

一种利用带限思想计算重力异常高阶径向导数的方法
技术领域
本发明属于物理大地测量技术领域,尤其是一种利用带限思想计算重力异常高阶径向导数的方法。
背景技术
高精度的重力异常高阶径向导数(即四阶径向导数、五阶径向导数和六阶径向导数)是推求地球内部重力和潜器水下重力辅助导航的基础信息,同时可反演近地表异常物体和确定矿产资源分布规律。
计算重力异常高阶径向导数的解析核函数在球面边界存在奇异性,导致计算结果不稳定;位场理论的球谐表达式具有较好的稳定性,为克服计算重力异常高阶径向导数的解析核函数的奇异性问题提供了新途径。在实际应用中受观测数据覆盖范围限制,无法做到全球覆盖,实际计算过程中重力异常高阶径向导数的全球积分式需要进行适用观测数据保障条件的改化,以保证计算结果的可靠性。
综上所述,如何克服重力异常高阶径向导数计算结果奇异性问题以提高重力异常高阶径向导数的计算精度是目前迫切需要解决的问题。
发明内容
本发明的目的在于克服现有技术的不足,提出一种设计合理、精度高且稳定性强的利用带限思想计算重力异常高阶径向导数的方法。
本发明解决其技术问题是采取以下技术方案实现的:
一种利用带限思想计算重力异常高阶径向导数的方法,包括以下步骤:
步骤1、将地球外部重力异常Poisson积分式的解析核函数表示为球谐级数展开式,通过求导方法得到重力异常高阶径向导数球谐级数表达式;
步骤2、利用移去-恢复技术将重力异常高阶径向导数球谐级数表达式截断为与重力观测值频谱范围相一致的带限求和式,并利用全球重力位模型补偿重力异常高阶径向导数远区截断误差,计算得到重力异常高阶径向导数。
而且,所述步骤1得到的重力异常高阶径向导数球谐级数表达式为:
Figure BDA0002941955180000011
式中,Δg为外部空间计算点
Figure BDA0002941955180000012
重力异常;ΔgR为球面上流动点
Figure BDA0002941955180000013
处的已知观测重力异常;R为地球椭球平均半径;r为计算点地心向径;
Figure BDA0002941955180000014
为计算点的纬度和经度;
Figure BDA0002941955180000015
为流动点的纬度和经度;σ为单位球面;dσ为单位球面的面积元;ψ为计算点至流动点之间的球面角距;
Figure BDA0002941955180000021
是计算点至流动点之间的空间距离;K(r,ψ)为积分核函数;Pn(cosψ)为n阶勒让德多项式级数。
而且,所述步骤2的具体实现方法为:利用移去-恢复技术,从重力异常观测值中移去参考重力异常,得到残差重力异常;从积分式核函数中移去与参考场对应阶次的核函数球谐表达式,得到截断核函数,使之与残差重力异常的频谱匹配;基于截断核函数和残差重力异常的局域积分得到残差重力异常高阶径向导数;利用全球重力场位模型高阶信息进行远区效应补偿;恢复参考重力异常高阶径向导数对重力异常高阶径向导数球谐级数表达式改化处理,得到重力异常高阶径向导数的计算公式并进行计算。
而且,所述重力异常高阶径向导数的计算公式为:
Figure BDA0002941955180000022
式中δΔgR为ΔgR的残差重力异常;g4pref、g5pref和g6pref分别是参考重力异常四阶、五阶和六阶径向导数;
Figure BDA0002941955180000023
Figure BDA0002941955180000024
分别是重力异常四阶、五阶和六阶径向导数的改化核函数;
Figure BDA0002941955180000025
Figure BDA0002941955180000026
分别是四阶、五阶和六阶径向导数的远区补偿。
本发明的优点和积极效果是:
本发明设计合理,其将地球外部重力异常Poisson积分式的解析核函数表示为球谐级数展开式并通过求导方法得到重力异常高阶径向导数球谐级数表达式;利用移去-恢复技术将重力异常高阶径向导数球谐级数表达式截断为与重力观测值频谱范围相一致的带限求和式,同时利用全球重力位模型补偿重力异常高阶径向导数远区截断误差,精确计算重力异常高阶径向导数,提高了重力异常高阶径向导数的计算精度,解决了重力异常高阶径向导数计算结果的不稳定性问题,可广泛用于物理大地测量领域。
具体实施方式
一种利用带限思想计算重力异常高阶径向导数的方法,包括以下步骤:
步骤1、将地球外部重力异常Poisson积分式的解析核函数表示为球谐级数展开式,通过求导方法得到重力异常高阶径向导数球谐级数表达式。
在本步骤中,地球外部重力异常Poisson积分式的解析计算式为:
Figure BDA0002941955180000027
Figure BDA0002941955180000028
式中,Δg为外部空间计算点
Figure BDA0002941955180000031
重力异常;ΔgR为球面上流动点
Figure BDA0002941955180000032
处的已知观测重力异常;R为地球椭球平均半径;r为计算点地心向径;
Figure BDA0002941955180000033
为计算点的纬度和经度;
Figure BDA0002941955180000034
为流动点的纬度和经度;σ为单位球面;dσ为单位球面的面积元;ψ为计算点至流动点之间的球面角距;
Figure BDA0002941955180000035
是计算点至流动点之间的空间距离;K(r,ψ)为积分核函数。
重力异常高阶径向导数可表示为:
Figure BDA0002941955180000036
将解析核函数K(r,ψ)用球谐级数展开式表示:
Figure BDA0002941955180000037
式中Pn(cosψ)为n阶勒让德(Legendre)多项式级数。
对式(4)求高阶径向导数,本发明取四阶、五阶和六阶为高阶;并令r=R,可得球面上的核函数径向导数为:
Figure BDA0002941955180000038
将式(5)代入式(3),可得到球谐级数展开表达的重力异常高阶径向导数积分计算式:
Figure BDA0002941955180000039
步骤2、依据各类重力异常观测经滤波处理后均表现为一类有限频谱带宽信号的特性,顾及远区效应影响,利用移去-恢复技术将重力异常高阶径向导数球谐级数表达式截断为与重力观测值频谱范围相一致的带限求和式,同时利用全球重力位模型补偿重力异常高阶径向导数远区截断误差,计算得到重力异常高阶径向导数。
本步骤的具体实现过程如下:
由式(6)表达的重力异常高阶径向导数计算模型是理论计算式,计算重力异常高阶径向导数要求全球积分,但在实际应用中受观测数据覆盖范围限制,无法做到全球覆盖,实际计算过程中重力异常高阶径向导数的全球积分式需要进行适用观测数据保障条件的改化,以保证计算结果的可靠性。本发明顾及实测数据局域保障条件,引入全球重力场位模型,利用移去-恢复技术,从重力异常观测值中移去参考重力异常,以得到残差重力异常;从积分式核函数中移去与参考场对应阶次的核函数球谐表达式,以得到截断核函数,使之与残差重力异常的频谱匹配;基于截断核函数和残差重力异常的局域积分得到残差重力异常高阶径向导数;利用全球重力场位模型高阶信息进行远区效应补偿,以削弱远区截断误差的影响;恢复参考重力异常高阶径向导数对式(6)改化处理,得到计算重力异常高阶径向导数的计算公式,具体计算公式为:
Figure BDA0002941955180000041
式中δΔgR为ΔgR的残差重力异常;g4pref、g5pref和g6pref分别是参考重力异常四阶、五阶和六阶径向导数;
Figure BDA0002941955180000042
Figure BDA0002941955180000043
分别是重力异常四阶、五阶和六阶径向导数的改化核函数;
Figure BDA0002941955180000044
Figure BDA0002941955180000045
分别是四阶、五阶和六阶径向导数的远区补偿。下面分别给出具体表达式。
式(7)中,残差重力异常δΔgR的计算式为:
δΔgR=ΔgR-ΔgRref (8)
式中,ΔgRref是球面上的参考重力异常,具体计算式为:
Figure BDA0002941955180000046
式中L是参考场位模型阶次;Δgn(θ,λ)是重力异常n阶拉普拉斯面球谐函数,具体表达式为:
Figure BDA0002941955180000047
式中(θ,λ)是计算点的余纬和经度;GM为地球引力常数;
Figure BDA0002941955180000048
为完全规格化缔合勒让德函数;
Figure BDA0002941955180000049
Figure BDA00029419551800000410
为完全规格化地球位系数。
式(7)中参考重力异常四阶、五阶和六阶径向导数g4pref、g5pref和g6pref的具体表达式为:
Figure BDA00029419551800000411
式(7)中重力异常四阶、五阶和六阶径向导数的改化核函数
Figure BDA0002941955180000051
Figure BDA0002941955180000052
的具体表达式为:
Figure BDA0002941955180000053
式中N为与重力异常观测值滤波尺度对应的最高频谱阶次。
式(7)中四阶、五阶和六阶径向导数的远区补偿
Figure BDA0002941955180000054
Figure BDA0002941955180000055
的具体表达式为:
Figure BDA0002941955180000056
其中:
Figure BDA0002941955180000057
Figure BDA0002941955180000058
下面采用全球重力场位模型EGM2008作为标准场开展数值计算检验及分析比较。
选取重力场变化比较剧烈的马里亚纳海沟作为试验区,具体覆盖范围为:(
Figure BDA0002941955180000059
11°N~14°N;λ:143°E~146°E)。选取r=R+h,R=6371km,由EGM2008模型分别计算3组分别对应于h0=0km、h6=6km、h10=10km高度面上的1′×1′网格剩余重力异常“真值”
Figure BDA00029419551800000510
(i对应0km,6km,10km)和剩余重力异常四阶、五阶和六阶径向导数“真值”
Figure BDA00029419551800000511
(j=4,5,6)。表1列出了3组不同高度面1′×1′网格剩余重力异常“真值”的统计结果,表2列出了相对应3组高度面四阶、五阶和六阶径向导数“真值”的统计结果。
表1不同高度面EGM2008模型重力异常统计结果(mGal)
高度面(km) 最小值 最大值 平均值 均方根值
0 -78.48 132.75 -0.05 26.36
6 -41.18 74.21 -0.04 16.22
10 -30.45 52.29 -0.04 12.00
表2不同高度面EGM2008模型重力异常径向偏导数统计结果
Figure BDA0002941955180000061
为了检验本发明的计算效果,首先以前面选定的3个高度面(h0=0km、h6=6km、h10=10km)上的位模型剩余重力异常
Figure BDA0002941955180000062
作为观测量,依据本发明公式(7)计算相对应高度面上的1′×1′网格四阶、五阶和六阶径向偏导数
Figure BDA0002941955180000063
将计算值
Figure BDA0002941955180000064
与相对应的“真值”
Figure BDA0002941955180000065
作比较,可获得不同高度面不同阶次偏导数计算模型的精度评估信息,具体比对统计结果列于表3。这里积分半径统一取为ψ0=0.5°,故计算区域外围0.5°范围内的比对结果不参加精度评估统计计算(下同)。为了对比分析评价径向偏导数积分模型改化前后的计算效果,本试验同时给出了采用传统算法(公式3)完成相同参量计算获得的精度评估结果,具体见表4。
表3利用本发明计算不同高度面重力异常高阶径向导数精度检核
Figure BDA0002941955180000066
表4利用原始模型计算不同高度面重力异常高阶径向导数精度检核
Figure BDA0002941955180000071
从表3检核结果可以看出,依据本发明计算重力异常四阶至六阶径向导数,均可获得比较满意的符合精度。从表3并结合表2统计结果可以看出,本发明的绝对精度(互差均方根值)都随计算高度面的增高和偏导数阶数的增大而提升,其相对精度(互差均方根值/径向导数均方根值)的变化趋势则正好相反,均随计算高度面的增高和偏导数阶数的增大而降低。计算高度面越高,其相对精度的降低幅度越明显。这个结果显然跟高度面越高,高阶偏导数的绝对量值越小有关,同时跟偏导数阶数越高,积分模型离散化误差的影响越大有关,这是符合理论分析预期的结果。进一步对比表3和表4统计结果可知,本发明结果精度明显优于传统算法,表明本发明实用易行,具有较高的应用价值。
需要强调的是,本发明所述的实施例是说明性的,而不是限定性的,因此本发明包括并不限于具体实施方式中所述的实施例,凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式,同样属于本发明保护的范围。

Claims (2)

1.一种利用带限思想计算重力异常高阶径向导数的方法,其特征在于:包括以下步骤:
步骤1、将地球外部重力异常Poisson积分式的解析核函数表示为球谐级数展开式,通过求导方法得到重力异常高阶径向导数球谐级数表达式;
步骤2、利用移去-恢复技术将重力异常高阶径向导数球谐级数表达式截断为与重力观测值频谱范围相一致的带限求和式,并利用全球重力位模型补偿重力异常高阶径向导数远区截断误差,计算得到重力异常高阶径向导数;
所述步骤2的具体实现方法为:利用移去-恢复技术,从重力异常观测值中移去参考重力异常,得到残差重力异常;从积分式核函数中移去与参考场对应阶次的核函数球谐表达式,得到截断核函数,使之与残差重力异常的频谱匹配;基于截断核函数和残差重力异常的局域积分得到残差重力异常高阶径向导数;利用全球重力场位模型高阶信息进行远区效应补偿;恢复参考重力异常高阶径向导数对重力异常高阶径向导数球谐级数表达式改化处理,得到重力异常高阶径向导数的计算公式并进行计算;
所述重力异常高阶径向导数的计算公式为:
Figure FDA0003589384650000011
式中δΔgR为ΔgR的残差重力异常;g4pref、g5pref和g6pref分别是参考重力异常四阶、五阶和六阶径向导数;
Figure FDA0003589384650000012
Figure FDA0003589384650000013
分别是重力异常四阶、五阶和六阶径向导数的改化核函数;
Figure FDA0003589384650000014
Figure FDA0003589384650000015
分别是四阶、五阶和六阶径向导数的远区补偿。
2.根据权利要求1所述的一种利用带限思想计算重力异常高阶径向导数的方法,其特征在于:所述步骤1得到的重力异常高阶径向导数球谐级数表达式为:
Figure FDA0003589384650000016
式中,Δg为外部空间计算点
Figure FDA0003589384650000017
重力异常;ΔgR为球面上流动点
Figure FDA0003589384650000018
处的已知观测重力异常;R为地球椭球平均半径;r为计算点地心向径;
Figure FDA0003589384650000019
为计算点的纬度和经度;
Figure FDA00035893846500000110
为流动点的纬度和经度;σ为单位球面;dσ为单位球面的面积元;ψ为计算点至流动点之间的球面角距;
Figure FDA00035893846500000111
是计算点至流动点之间的空间距离;K(r,ψ)为积分核函数;Pn(cosψ)为n阶勒让德多项式级数。
CN202110180193.3A 2021-02-08 2021-02-08 一种利用带限思想计算重力异常高阶径向导数的方法 Active CN112987118B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110180193.3A CN112987118B (zh) 2021-02-08 2021-02-08 一种利用带限思想计算重力异常高阶径向导数的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110180193.3A CN112987118B (zh) 2021-02-08 2021-02-08 一种利用带限思想计算重力异常高阶径向导数的方法

Publications (2)

Publication Number Publication Date
CN112987118A CN112987118A (zh) 2021-06-18
CN112987118B true CN112987118B (zh) 2022-05-27

Family

ID=76392837

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110180193.3A Active CN112987118B (zh) 2021-02-08 2021-02-08 一种利用带限思想计算重力异常高阶径向导数的方法

Country Status (1)

Country Link
CN (1) CN112987118B (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102661716A (zh) * 2012-04-20 2012-09-12 武汉理工大学 基于光纤陀螺技术的桥梁和隧道线形及刚度检测方法与系统

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102236109B (zh) * 2010-05-06 2013-04-24 中国石油天然气集团公司 一种重磁干扰区的方差系数干扰校正方法
CN102590856A (zh) * 2011-01-11 2012-07-18 中国科学院地质与地球物理研究所 基于小波频谱分析的位场异常分离方法
FI126877B (en) * 2014-07-08 2017-07-14 Helsingin Yliopisto Formation of Gravity Estimate
CN108387951B (zh) * 2018-01-19 2019-11-12 中国人民解放军92859部队 一种基于重复线校正海空重力仪格值的新算法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102661716A (zh) * 2012-04-20 2012-09-12 武汉理工大学 基于光纤陀螺技术的桥梁和隧道线形及刚度检测方法与系统

Also Published As

Publication number Publication date
CN112987118A (zh) 2021-06-18

Similar Documents

Publication Publication Date Title
CN108415879B (zh) 基于向上延拓的航空重力最小二乘向下延拓解析方法
CN110632674B (zh) 一种航空重力测量数据的弱信息提取方法
CN108761510B (zh) 利用基于地形改正的重力场模型进行水准高差测量的方法
CN107504974B (zh) 地形分块与地形测点加权的地形匹配定位方法
CN113189559B (zh) 一种星载成像高度计遥感数据海底地形反演方法
CN108319566B (zh) 基于向上延拓的航空重力点对点向下延拓解析方法
CN111257956A (zh) 一种基于Matlab的区域似大地水准面精化方法
CN112949049B (zh) 一种利用带限思想计算重力异常低阶径向导数的方法
CN112965124B (zh) 一种顾及局域保障条件计算外部重力异常垂直梯度的方法
CN112987118B (zh) 一种利用带限思想计算重力异常高阶径向导数的方法
CN113419280B (zh) 基于改进椭圆拟合的叠前裂缝密度估算方法
CN112965127B (zh) 一种基于重力异常计算外部扰动重力径向分量的方法
CN113900069A (zh) 一种基于干涉成像高度计的垂线偏差计算方法及其系统
CN112818285B (zh) 一种计算外部扰动重力北向分量中央区效应的方法
CN114089432B (zh) 一种利用卫星测高数据反演海洋重力梯度的频率域方法
CN115236759B (zh) 一种确定地球重力场的六边形网格剖分方法
CN115406401A (zh) 减小矿区测量高程异常差值的方法
CN112965125B (zh) 一种基于重力异常计算外部扰动重力东向分量的方法
CN112965128B (zh) 一种无奇异性顾及局域保障条件计算外部重力异常的方法
CN112965123B (zh) 一种基于重力异常计算外部扰动重力北向分量的方法
CN112965126B (zh) 一种计算外部扰动重力东向分量中央区效应的方法
CN111721272A (zh) 一种基于椭球面计算的工程表面测量方法
CN112836378B (zh) 基于Poisson理论计算外部重力异常垂直梯度中央区效应的方法
CN113985491B (zh) 一种基于多源数据的重力场模型精化方法及系统
CN117128921A (zh) 一种陆地区域高程异常的确定方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant