地球物理大地测量
科学
计算软件

CASM.jpg

重力场及大地水准面精化系统PALGrav4.0

重力场及大地水准面精化系统(Precise Approach for Local Gravity field and Geoid) PALGrav4.0是一种用于地球物理大地测量科学计算的Windows程序包,其主要功能包括:地面、海洋、航空及其他近地空间重力场数据处理,地球外部各类扰动重力场元、多种性质地形影响的精密计算,任意高度各种类型扰动重力场积分正反算与解析延拓,稳态局部重力场及大地水准面精化,以及基于重力场数据的区域高程基准优化与应用计算。

PALGrav4.0由40多个win64程序和200余个功能模块组成。为方便自学、教学和技术培训,大多程序配有样例,包括完整计算流程、输入输出数据和计算过程截图。完成全部样例练习约需3个工作日。

PALGrav4.0采用3种自定义格式数据,即离散点值记录文件、格网格值文件和向量格网文件。[离散点值记录标准化]程序,是系统接受外部文本数据的唯一接口。其他程序模块,只接受系统本身产生或转换的自定义格式数据。

PALGrav4.0科学技术特色

PALGrav4.0核心专业算法基于大地测量与稳态地球重力场理论方法,能有效融合大地水准面外部不同高度、交叉分布、多源异质陆海重力场数据,科学目标是高精度局部重力场逼近和1cm水平稳态大地水准面精化。主要技术特点表现为:

(1)构造大地水准面外部多种解析调和的地形影响场,研发地面、海洋、航空及大地水准面外部各类扰动场元、多种性质的地形影响算法体系,进而运用地形影响移去恢复法,提升复杂多源重力场数据融合与局部重力场逼近算法的稳定性和抗差能力。

(2)严格构造重力场边值问题的积分解条件,采用积分运算代替微分运算,研发功能较为完备的重力场积分算法体系,实现大地水准面外部任意高度上各类扰动场元的正反积分运算与向上向下解析延拓,全面提升局部重力场及大地水准面的精化性能。

(3)依据地球重力场与大地测量基本原理,研发了一组基于局部重力场数据和方法的物理大地测量特色算法,优化区域高程基准,拓展局部重力场及高程基准成果的应用服务水平。

(4)基于积分解的唯一性和解析调和性,构造以局部重力场数值模型为参考场的移去恢复法,实现局部重力场及大地水准面的区域尺度逐级控制逼近和小面积局部精化更新,从而消除所谓的大地水准面衔接问题,有效解决大地水准面模型的历史继承问题。

物理大地测量学中地形影响概念

扰动场元的地形影响涉及3个关键要素:①地形或地壳质量调整方式;②扰动场元类型(地形影响对象);③扰动场元所在位置(与地形质量的位置关系)。

按质量调整方式不同,地形影响类型通常包括:局部地形影响,地壳均衡影响,布格质量影响,海水布格补偿影响,地形Helmert凝聚和剩余地形影响等。

当地球外部地形影响调和时,可类似于参考重力场移去恢复方案,采用地形影响移去恢复法,提高局部重力场逼近性能和算法的稳定性。

PALGrav4.0中,各种类型场元的地形影响等于其地形改正的负值,如局部地形影响等于负的局部地形改正。

地形影响优选准则

局部重力场逼近中,场元的地形影响处理有两个基本目的(注意区别于地球物理勘探目的):一是离散扰动场元格网化,另一个是扰动重力场积分时用于分离地形超短波成分。

为提高离散场元格网化精度,要求移去地形影响后,离散场元平滑度有所提高。此时地形影响的优选准则为:移去地形影响后,离散场元的标准差有所下降。

扰动场元积分时要求地形影响仅有超短波成分,因而优选准则为:地形影响移去前后,扰动场元标准差有所减小,且地形影响在数十公里范围内平均值很小。

大地水准面唯一性及PALGrav4.0实现

Stokes边值问题要求大地水准面外部没有质量,需在地球外部扰动位恒定不变条件下将地形质量压缩到大地水准面内部。有一种地形质量压缩方式,压缩后地球外部扰动位不变,而地面到大地水准面间的扰动位等于外部扰动位的解析延拓值,从而获得大地水准面的解析延拓解。

由PALGrav4.0计算的扰动重力地形影响与高程异常地形影响,满足Hotine积分,如扰动重力地形Helmert凝聚(直接影响)与高程异常地形Helmert凝聚(间接影响)满足Hotine积分公式。因此,无论选择的是局部地形影响、地形Helmert凝聚,还是剩余地形影响,由PALGrav4.0程序按地形影响移去恢复法精化的区域大地水准面,都是大地水准面的解析延拓解。

显然,卫星重力场确定或地球位系数模型直接计算的大地水准面高,是大地水准面的解析延拓解。PALGrav4.0通过维持大地水准面解的唯一性,能有效深度融合卫星重力、重力场位系数模型和区域重力场数据,严密实现局部重力场及大地水准面的高精度逼近。

积分解的唯一性与局部参考场数值模型

重力场积分解是唯一的,能保证局部重力场的调和不变性和空间结构的解析不变性。局部重力场的数值积分解在区域尺度上与重力场位系数模型在全球尺度上有着相同的性质和功效。

基于积分解的唯一性和解析调和性,能构造以局部重力场数值模型为参考场的移去恢复法,实现局部重力场及大地水准面的区域尺度逐级控制逼近和小面积局部精化更新,解决大地水准面无缝衔接和历史继承问题。

重力位、大地水准面与高程基准

(1)重力大地水准面是地球重力场边值问题解,重力位为常数,等于正常椭球面(大地水准面高的起算面)的正常重力位。

(2)全球大地位W₀通常用于表示全球高程基准(IERS数值标准),可依据Gauss大地水准面定义,用最新地球重力场模型和海面高观测计算。

(3)零正常高面与零正高面处处重合,即无论是正高还是正常高系统,高程基准面都是同一个,其重力位恒等于高程基准零点的重力位。

(4)全球大地位W₀值与(全球或区域)重力大地水准面的重力位(等于正常椭球面的正常重力位),不存在任何直接的大地测量关系。

有利于高程基准的解析正高系统

(1)令地面到大地水准面间流动点的重力值等于地球外部重力场解析延拓到流动点的重力值,由此得到的正高称为解析正高。全球地面点的解析正高与Helmert正高最大互差不超过1cm。

(2)大地水准面高是大地水准面大地高,是Stokes边值问题在地球坐标系中的解,其度量尺度为地球坐标系的几何尺度;解析正高高差在垂直方向的度量尺度也严格等于地球坐标系的几何尺度。可见解析正高与地球坐标系和大地水准面高的几何尺度一致。

(3)确定解析正高无需地形密度假设,且可用最新地球重力场数据不断精化。从大地测量基准的唯一性、可重复性和可测性角度上考察,解析正高比其他类型正高更适合高程基准目的。

模型重力场元点值计算

 功能:由系统内部地球重力场位系数模型,计算地球空间任意点处的高程异常(m)、空间异常(mGal)、扰动重力(mGal)、垂线偏差向量(sʺ/秒,南向、西向)、(垂向)扰动重力梯度(E)、水平重力梯度向量(E,北向、东向)或扰动位(m²/s²)模型值。

输入:地面或地球外部空间计算点文件。

参数设置:如图1

1)选择计算的模型重力场元类型。

2)输入计算点记录起始行号和大地高属性所在列序号。

3)输入重力场位系数模型最大计算阶数。程序自动选择地球重力场位系数模型最大阶数和输入最大阶数中的最小值作为计算阶数。

(4)利用EGM2008重力场模型和系统默认的正常椭球参数,最大计算阶数取1800阶,计算模型地面高程异常(m)、模型地面扰动重力(mGal)、模型地面垂线偏差(ʺ)和模型地面(垂向)扰动重力梯度(E)格网,如图2~5。

陆海地形质量球谐分析

 功能:由全球陆海数字高程模型,计算陆地和海水完全布格的地形质量(用面密度表示),对其进行球谐分析,将其用相应阶次的规格化球谐系数模型(kg/m²)表示。

输入:(1)全球陆海数字高程模型格网文件。程序要求全球陆海数字高程模型的格网经纬度为经度、地心纬度。

(2)球谐系数阶数n等于陆海数字高程模型在纬度方向上的格网数。如1˚分辨率陆海数字高程模型对应n=180,0.25˚陆海数字高程模型对应n=720。

参数设置:如图1

(1)残差标准差阈值为残差标准差与陆海地形面密度标准差之比的百分数,迭代终止条件为前后两次迭代结果增量的标准差。满足两个条件之一时,迭代计算完成。

(2)当选择输出地形位系数模型时,程序输出的球谐系数为地形质量生成的位系数,称为地形位系数(m)。

输出:全球陆海地形面密度规格化球谐系数文件,残差陆海地形面密度格网数字模型文件,主界面输出的信息(图2)。

说明:(1)由0.25˚×0.25˚全球陆海数字高程模型(图3),计算720阶地形面密度规格化球谐系数。图5是从源0.25˚×0.25˚全球陆海数字高程模型中截取的区域陆海数字高程模型,图4由720阶地形面密度球谐系数计算的区域陆海数字高程模型。

2PALGrav3.0安装目录下\data\地形面密度球谐系数1800.dat,是利用6ʹ×6ʹETOP全陆海数字高程模型,按本程序计算的1800阶地形面密度球谐系数模型文件,地形高度模型残差标准差7.2m

(3)程序可计算直到5400阶地面面密度球谐系数,对于5400阶球谐系数,要达到米级精度地形高度残差标准差,约需30小时迭代计算时间。

广义Hotine数值积分

 功能:由等位面大地高数字模型(m)及其面上的扰动重力数字模型(mGal),按严密广义Hotine积分公式,计算地面及地球外部点的高程异常(m)

输入:(1)地球外部空间计算点文件。

(2)等位面残差扰动重力数字模型(mGal)文件,图2与图3之差。

(3)格网规格完全相同的等位面大地高数字模型(m)文件。等位面大地高数字模型可由“过指定点模型等位面构造”生成。

参数设置:如图1

(1)输入计算点记录起始行号和计算点大地高属性所在列序号。

(2)输入Hotine积分半径dr,10<dr≤300km。

⊙残差扰动重力 = 扰动重力 – 模型扰动重力 (– 扰动重力地形影响)。

输出:残差高程异常点值文件。在源计算点值文件记录的基础上增加一列该点的残差高程异常积分值,保留4位有效数字。

说明:(1)利用EGM2008模型,按“模型重力场元点值计算”,生成大地水准面上2.5'×2.5'残差扰动重力和残差地面高程异常格网模型(1800阶与360阶模型值之差)。

(2)采用200km积分半径,按广义Hotine数值积分,由大地水准面上残差扰动重力,得到残差地面高程异常积分值,如图4。

2020061220200612

(3)以EGM2008模型残差地面高程异常为参考真值,计算Hotine残差地面高程异常积分值与EGM2008模型残差地面高程异常之差,如图6。残差地面高程异常的统计结果如表1

高程基准零点重力位参数计算

功能:利用GNSS水准高程异常、残差高程异常等,估计区域高程基准零点重力位,进而由系统参数确定的椭球面正常重力位、全球大地位,计算区域高程基准零点重力位和高程基准转换参数。程序适用正高和正常高系统。

输入:GNSS水准高程异常(大地水准面高)成果点值文件。

参数设置:如图

(1)GNSS水准成果点值文件记录的起始行号。

(2)输入GNSS水准高程异常(大地水准面高)、残差高程异常(大地水准面高)和权值在记录中的列序号。

输出:区域高程基准参数文件。

20200612

局部地形影响数值积分

 功能:由地面数字高程模型、地面大地高数字模型,按严密积分公式,计算地面及地球外部点高程异常(m)、扰动重力(mGal)、空间异常(mGal)或扰动位(m²/s²)的局部地形影响。

输入:1)格网规格完全相同的地面数字高程模型、地面大地高数字模型,以及空间计算点文件。

2)空间计算点记录应含大地高属性值。程序由地面大地高模型格网构造流动点大地坐标,将空间计算点的经纬度、大地高作为计算点大地坐标。

3)地面数字高程模型属正(常)高系统,地面数字高程模型格网在程序中用于表示地形起伏。

参数设置:如图1

1)选择计算局部地形影响的扰动重力场元类型。

(2)输入计算点记录起始行号和计算点大地高属性所在列序号。

(3)输入局部地形影响积分半径,km。

输出:局部地形影响计算结果点值文件。在源计算点值文件记录的基础上增加一列或若干列局部地形影响积分值,保留4位有效数字。

加权基函数插值格网化

 功能:按输入的格网规格和选定的插值权函数形式及参数值,对离散点值数据进行乘权格网化运算。

输入:待格网化的离散点值文件

参数设置:如图1

(1)输入点值属性记录起始行号,设置待格网化的离散点属性列序号。

(2)权函数形式:余弦函数、高斯函数或指数函数。

(3)权值属性所在的列序号,选择离散点数据非等权插值。

①计算插值点数值时,程序将离散点记录属性中的权值与权函数(插值点与离散点之间距离为自变量的函数)相乘,作为离散点权值。

②不选择“离散点数据非等权插值”时,所有离散点权等于其权函数。

(4)设置格网分辨率和经纬度范围。

(5)输入权函数峰度参数。参数值越大,权函数值随距离衰减越快。

输出:输出格网保留4位有效数字,结果如图2。

说明:(1)加权基函数插值格网化,适合各种空间分布差异大、精度不均匀的地球物理观测数据格网化,自适应采样点密度,无边缘效应,有利于构造高性能地球物理场

(2)权函数峰度越小,插值邻近点数越大,格网化过程的低通滤波能力越强,边缘平滑度越强,对稀疏数据的插值能力也越强。

(3)离散重力场元格网化建议:地形代表性误差是重力数据格网化的主要误差来源。为此,用户可先计算离散重力点的局部地形影响,以局部地形影响为参考属性,调用“指定参考属性观测量权值计算”,计算离散重力点权值,然后调用本模块,构造高性能的格网平均重力场。

区域大地水准面成果精度评估

 功能:基于GNSS水准高程异常与重力场频域误差特性,由GNSS水准残差高程异常,估计并确定用于表达地面高程异常成果精度的两个误差指标和两个误差曲线,即重力地面高程异常差误差、实用地面高程异常内部误差、实用地面高程异常差误差曲线与GNSS水准高程异常差误差曲线。

输入:GNSS水准残差高程异常点值文件

参数设置:如图1

(1)GNSS水准残差点值文件记录的起始行号,残差高程异常列序号。

(2)输入区域GNSS水准点平均间距D。

(3)输入GNSS水准点两两组合后的分组数。

(4)输入GNSS基线大地高差的固定误差和比例系数误差。

输出:用于完整表达地面高程异常成果精度的误差指标及误差曲线文件。

说明:(1)以某地区为例,利用104个GNSS水准点,组合成N= 5356条边,按距离递增排列后分成M= 90组。104个GNSS水准残差高程异常和5356条边GNSS水准残差高程异常差的统计结果如表所示(单位:cm)。

20200612

(2)取GNSS水准点平均间距D= 24km,GNSS大地高差固定误差和比例误差系数分别为a= 15mm、b= 0.2mm/km,程序计算得到每千米正常高差误差估值0.0307cm/km,重力地面高程异常差的误差估值4.159cm,实用地面高程异常内部误差1.6045cm。

(3)由计算结果绘制3条误差曲线,如图2。图中显示:

⊙实用地面高程异常差误差(黑线),既不大于重力地面高程异常差误差(蓝线),也不大于GNSS水准高程异常差误差(棕黄线)。实用地面高程异常差的误差曲线总是在其余两个误差曲线的下方。

20200612

⊙在距离L*=105.8km处,GNSS水准高程异常和重力地面高程异常,对实用地面高程异常的精度贡献相当。小于L*时,GNSS水准高程异常的贡献大,大于L*时,重力地面高程异常的贡献大。

⊙实用高程异常差误差曲线(黑线)的斜率,随距离增大而减小,且不大于GNSS水准高程异常差误差曲线的斜率。当斜率接近零时,实用地面高程异常差的误差逼近重力地面高程异常差的误差。

⊙当距离L接近或小于GNSS水准点平均间距D= 24km时,GNSS水准高程异常差对实用地面高程异常差的贡献起主导控制作用。