地球物理大地测量
大型科学
计算平台

CASM.jpg

高精度重力场逼近与大地水准面计算系统PAGravf4.5

高精度重力场逼近与大地水准面计算系统(Precise Approach for Earth Gravity field and Geoid)PAGravf4.5,是一种基于稳态重力场理论的大地测量科学计算Windows程序包。主要由地球重力场数据分析与预处理计算,不同高度各类场元多种地形影响计算,高精度重力场逼近与全要素建模,区域高程基准优化、统一与应用计算,以及大地测量数据文件编辑计算与可视化五大子系统构成。

用户可事先根据计算目的与目标区域实际观测条件,调用PAGravf4.5程序,全面测试算法性能和局部重力场特征,设计高性能计算方案,再按需组织有关功能模块,完成陆海空天重力场数据处理与各种地形影响计算分析工作,进而有效综合多源异质重力场数据,按重力场边值理论确定高精度陆海统一的重力场及大地水准面,精化区域高程基准,开展重力场及高程基准应用服务计算及分析工作。

PAGravf4.5科学目标

        严格遵循地球重力场逼近理论要求,解决地球外部空间各种类型扰动重力场元的多种地形影响问题,实现边值理论框架中区域重力场全要素建模以及外部扰动场元解析关系的循环闭合,梳理完善一些经典概念和传统理念,拓展高程基准的物理大地测量应用能力与计算水平。

    (1)解决地球外部各类场元、多种地形影响的解析相容性与严密统一计算问题,以应对复杂观测情况下地质地球物理重力勘探要求,提升地球重力场数据处理水平。

    (2)在边值理论框架中构造场元关系循环闭合的重力场积分与球面径向基函数逼近算法体系,实现多种异质数据混叠的重力场全要素建模与1cm稳态大地水准面精化。

    (3)发展基于物理大地测量学原理的特色算法,深入挖掘地球重力场数据与方法在区域高程基准优化、统一与应用计算中的潜力,大力拓展其应用服务水平。

PAGravf4.5技术特色

    (1)构造大地水准面外部解析调和的多种地形影响场,研发适合空天地海不同高度、各种类型扰动场元、多种性质地形影响的解析算法体系,以满足复杂观测情形下地球物理重力勘探要求,全面提升地球重力场数据处理与逼近水平。

    (2)研发科学完备的各类扰动场元正反积分运算及径向基函数逼近算法体系,实现由某种边界面一种类型扰动场元,按边值问题解算地球外部各种类型扰动场元的循环闭合,有效提升重力场全要素建模能力及大地水准面精化水平。

    (3)按照物理大地测量理论与技术要求,提出地形影响优选定量准则,梳理重力场与高程基准的几何物理关联;研发基于地球重力场数据和理论方法的特色算法,优化区域高程基准,拓展地球重力场及高程基准成果的应用能力。

    (4)强大的测试分析功能。用户可根据观测条件和计算目标,充分测试验证目标区域重力场与计算技术路线的性能和特点,事先设计高性能计算方案。能解析融合极少天文垂线偏差或GNSS水准数据,具备直接测定离散重力数据外部精度的能力。

课堂教学与技术培训程序
完整的地形影响算法体系

为实现地球外部各种类型重力场元、不同性质的地形影响,PAGravf4.5发展了一套完整的地形影响算法体系。

(1)这套算法公式理论严密,数值积分无计算误差,快速算法的精度可控;

(2)地形影响性质多样,地形影响的重力场元类型可以是任意的;

(3)不同类型场元地形影响之间严格遵循重力场解析关系,算法公式精炼;

(4)充分运用不同性质地形影响之间的解析相容性,算法代码实现短小精悍。

如7.5~7.8节很多算法公式很相似。实际上,一些地形影响算法,只需调整一些参数,就可调用同一套算法代码实现。

重力场积分算法技术特点

(1)固定积分半径。PAGravf4.5通过控制核函数定义域,实现固定半径的重力场积分,包括数值积分与积分快速FFT算法(核函数加窗),以协调统一各种重力场逼近算法。二维FFT采用改化的平面二维核函数,在纬度10°范围内,其计算精度与一维FFT没有明显差异。

(2)计算点和流动点位置。点位大地坐标一律用经纬度和大地高表示,如边界面位置、测点、计算点、积分流动点(面元、体元)的位置一律用点位大地坐标表示,积分格网位置取单元格网中心点大地坐标,积分半径一律采用大地坐标计算。

(3)等重力位边界面。大多数重力场积分公式由Stokes边值问题导出,如Hotine积分、V-M积分,径向梯度积分公式等。Stokes边值问题解要求边界面是等位面,即扰动场元位于重力等位面。

PAGravf4.5程序中,用于边界面的大地高精度只要不低于10m就能满足要求。可用360阶重力位系数模型构造,近地空间可用等正(高)面大地高格网代替。

径向基函数重力场解析逼近

PAGravf4.5提出三个关键性技术措施,使得SRBF重力场逼近算法与观测场元误差无关,避免待定目标场元谱泄漏,实现SRBF重力场解析逼近。

(1)采用边缘效应抑制方法代替法方程正则化。

PAGravf4.5提出了一种通过抑制边缘效应来提高参数估计性能的算法。当球面径向基函数(SRBF)中心点位于计算区域边缘时,将其SRBF系数等于零作为观测方程,以提升SRBF系数参数估计的稳定性和可靠性。

引入边缘效应抑制方法后,法方程不再需要正则化,从而有效避免SRBF逼近算法的解析性质与重力场解析结构受观测量误差影响。

(2)多次累积SRBF逼近法实现重力场最佳逼近。

目标场元是观测场元与滤波器SRBF的卷积。当目标场元与观测场元类型不同时,一个SRBF难以同时与观测场元和目标场元的谱域中心及带宽有效匹配,势必导致目标场元的谱泄漏。而且,除Bjerhammar球埋藏深度(宽度参数)外,SRBF类型、最小最大阶数与SRBF中心分布也都影响重力场逼近性能。因此,仅以埋藏深度为参数的SRBF系数最优估计,不足以保证重力场的最佳逼近。

为解决这一关键问题,PAGravf4.5基于重力场的线性可加性,提出了多次累积SRBF逼近法重力场最佳逼近策略,代替以埋藏深度(或宽度参数)为参数的SRBF系数最优估计方案,不再要求交叉确认最佳埋藏深度(宽度参数)。

每次SRBF逼近采用不同谱域特征的SRBF,多次累积SRBF逼近通过组合多个SRBF谱域中心及带宽,就可充分解析目标场元的谱域信号,避免谱泄漏,从而在空域中实现对目标场元的最佳恢复。每次残差逼近本质是用上次逼近结果作参考重力场,按移去恢复法精化残差目标场元。

(3)多种异质数据融合的协因数阵对角线法。

方差分量估计法SRBF系数参数估计需要迭代计算,受观测误差影响,不但存在没有稳定解的风险,还严重削弱了粗差探测能力。PAGravf4.5提出一种适应各种复杂观测情形的协因数阵对角线标准差法,来实现多种异质观测场元融合的重力场逼近,以代替常见的方差分量估计法。该方法用协因数阵对角线标准差代替迭代过程中残差观测场元的方差,使得参数估计解的性质只与观测量空间分布有关,而不受误差影响,从而有利于融合空间分布存在极端差异的多种类型观测场元(如极少数天文垂线偏差或GNSS水准点数据),有利于精准探测观测量粗差。协因数阵对角线法的法方程无需迭代计算。

技术特色:①能融合极少天文垂线偏差或GNSS水准数据;②可抑制边缘效应;③能选择观测场元并调节其贡献;④重力场逼近算法解析,算法性能与观测场元误差无关;⑤观测场元可位于不同高度,不规则空间分布,无延拓和格网化运算的信息损失;⑥具备直接测定任意高度离散观测场元外部精度的能力(设贡献率κ=0)。

算法与计算技术路线优化

PAGravf4.5的算法体系设计科学严密,原则上可以采用多种方案,计算空天地海任意类型场元的多种性质地形影响;可以由某一边界面上任意一种类型扰动重力场元,按边值理论计算整个地球外部空间各种类型扰动重力场元。场元类型可以是传统类型,也可以是卫星跟踪卫星、卫星轨道摄动等非常见类型重力场元,适用空间可以是大地水准面及其外部整个地球空间。

一种性质的地形影响,也有多种方案可选。如计算陆海完全布格影响,PAGravf4.5就有三种方案和程序可供选用。

针对某特定计算目标,可以选择多个不同PAGravf4.5程序算法、多种不同参数或多条技术路线来实现。实际计算时,应深入考察目标区域重力数据条件和重力场性质,仔细选择、测试和分析合适的PAGravf4.5算法和参数,优化技术路线。

算法性能及参数测试分析

(1)地形影响算法性能测试分析

PAGravf4.5提出的地形影响优选准则,基于物理大地测量基本原理,能大幅降低地形影响分析的复杂性,为有效发挥地形影响在地质地球物理重力勘探和重力场逼近中关键作用提供了具体可行的技术路线。

地形影响的性质随计算区域的地形起伏、重力场结构与重力点分布不同,都会存在明显差异。PAGravf4.5地形影响子系统案例中,选择了某个困难山区,对各种场元各类地形影响的最大最小值与标准差之比进行统计分析,结果显示:不考虑数据分布和重力场结构,局部地形影响对重力数据处理有利,地形Helmert凝聚对重力梯度数据处理有利,而剩余地形影响对精化大地水准面更有利。

若结合具体重力数据资源,就可以分析该区域重力场结构与重力点分布的影响,此时需要进一步按定量准则(1)(2)进行统计分析,结论可能会产生变化。

实际计算前,应针对目标区域的地形、重力场特征与可用重力资源情况,依据地形影响定量准则,全面细致地测试分析地形影响技术路线,保证地形影响算法及参数选择有据可依,才有可能显著提升地形影响数据处理方案的适用性和技术水平。

(2)重力场逼近算法性能测试

大多数重力场逼近算法和参数设置性能,可以用超高阶地球重力场位系数模型自行测试验证。PAGravf4.5许多程序样例,也是以EGM2008重力场位系数模型的540阶为参考重力场,采用541~1800阶作为残差扰动重力场,进行自测试(参考样例)。

PAGravf4.5程序算法测试要点:以部分类型残差扰动场元为观测量,调用待测试的PAGravf4.5程序或功能,得到另一部分类型残差场元的计算值。通过对比分析算法程序计算的残差扰动场元计算值,与残差扰动场元模型值(作为参考值)的差异,评价PAGravf4.5中算法程序和参数设置的技术性能。

PAGravf4.5能由一种位于某个边界面上的场元,计算大地水准面外部整个地球空间的各种类型场元,可以循环计算该边界面上同类型场元。通过比较分析边界面上观测场元与循环计算得到的计算场元之间的差异与相似性,可以分析计算流程中调用的有关程序和功能的算法特点和性能。

(3)径向基函数逼近性能测试

调用PALGrav4.5程序[径向基函数重力场逼近及性能],通过输入不同单一类型观测场元、选择不同类型径向基函数及其参数、计算不同类型目标场元,可全面测试径向基函数及其重力场逼近算法的空域、谱域和解析性质,以揭示目标区域重力场空域谱域精细结构。考察观测场元、径向基函数与目标场元的谱域中心及带宽,以充分解析目标场元谱结构为原则,结合实际观测条件,可优化设计区域重力场最佳逼近方案。

PAGravf4.5程序[多种异构观测场元径向基函数逼近],本身具有很强的测试、验证和分析功能,可以测试验证极少高精度天文垂线偏差或GNSS水准数据的解析融合(令其贡献率κ>1)、多源多代混叠数据的粗差探测与外部精度评估(令κ=0),浅水卫星测高的边界流探测与海面地形分离(令κ<1),以及多种数据空间分布、质量与精度差异大等极端情况下重力场逼近的性能特点和技术水平。

(4)PAGravf4.5自带案例说明

PAGravf4.5自带的地形影响案例,选择平均海拔4000m、地形起伏超过3000m困难山区,以便展示地形影响及其算法细节特征;同样,PAGravf4.5重力场积分案例,选择重力场短波信号非富(移去540阶模型值后的残差扰动重力,空间变化超过300mGal)的复杂特征地区,以便展示局部重力场及其逼近算法的细节特征。

这些案例的统计结果,大致显示了相应算法的基本性能和水平。由于案例本身以介绍计算流程为主要目的,并没有对算法本身及其主要参数进行优化与分析,更大的潜力,还需用户结合具体情况进一步挖掘。

模型重力场元点值计算

 功能:由系统内部地球重力场位系数模型,计算地球空间任意点处的高程异常(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水准高程异常差对实用地面高程异常差的贡献起主导控制作用。