1、重力积场积分核函数的解析性质重力场积分公式中的积分核函数变量是被积场量所在边界面的积分流动面元到计算点的直线距离可用流动点与计算点之间的球面角距中表示。由于重力场积分核函数存在不收敛问题导致重力场积分算法普遍存在频谱泄漏问题。当计算点也位于边界面上时重力场积分还存在奇异性问题而奇异积分值的数值计算结果一般随边界面上被积场量的空间分辨率不同出现明显跳变导致局部重力场数值积分计算结果容易产生不确定性。Poisson积分式是推导Stokes边值积分解和Hotine边值积分解的理论基础也是推导Vening-Meinesz积分公式的依据。然而由于Possion核函数存在严重的高阶振荡与不收敛问题特别是球边界面上计算点的Possion核无定义从而导致空域重力场积分的数值计算结果存在难以控制的不确定性。为应对此类问题通常假设(残差)扰动场量随空间变化具有各向同性的随机统计性质或者借助边值问题谱域逼近解的近似改化重力场积分核函数。这类改化算法在以往数十年间一度成为重力场数值积分的重要技术手段然而由于普遍引入非重力场解析关系假设如观测场量空间统计假设或数据驱动参数估计法等这些改进方法的普适性缺乏确切的重力场理论依据PAGravf4.5因此不再推荐重力场积分核函数改化算法。重力场积分数值计算虽存在某些缺陷然而重力场积分公式显式直观地表达了重力场量之间的解析函数关系因而空域重力场边值理论仍然也一直是认识和理解重力场逼近理论和物理大地测量学原理的重要理论基础。2、SRBFSpherical Radial Basis Function球面径向基函数重力场逼近算法基于单一边界面上单一类型观测扰动场量边值条件的空域重力场边值理论无法适应空天地海多种异构重力场探测技术的需求。谱域重力场逼近理论不涉及边界面和边值条件问题可直接由多种类型观测场量按最小二乘法逼近重力场因而成为现代重力场逼近理论的主流方法。SRBFSpherical Radial Basis Function球面径向基函数重力场逼近算法的核心原理是将地球局部重力场表达为一系列定义在球面上的基函数的线性组合通过求解这些基函数的系数来最优地拟合观测数据。你可以把它想象成用许多不同大小、位置和形状的“拼图块”即基函数去拼凑出一个复杂的重力场图像。这种方法尤其适用于构建高精度的局部重力场模型因为它弥补了传统球谐函数Spherical Harmonics 在表达局部细节方面的不足。下面介绍球面径向基函数谱域重力场逼近理论和方法。与地球外部扰动位球谐级数形式对比这样好理解一点。3、适合重力场逼近的球面径向基函数用于重力场逼近的SRBF应满足Laplace方程。常见的点质量核函数、Poisson核函数、径向多极子核函数和Poisson小波核函数都是具有调和性质的球面径向基核函数。依据扰动重力场量定义扰动位的球面径向基函数展开式(最右边表达式)导出其他各种类型扰动场量的球面径向基函数参数化形式。对于局部重力场逼近目的通常借助参考重力场模型移去中低阶扰动重力场的模型值由区域观测残差场量逼近残差场元。此时式(10.16)~式(10.21)中阶数n的最小最大值范围(重力场的谱域带)和自变量中k定义域与选择的参考重力场、目标区域重力场结构及其空间分辨率(SRBF 节点数K或球面格网等级Q)都密切相关需采用实际数据测试分析后确定。可以采用全球和区域一致的球面等面积Reuter格网根据给定的格网等级Q构造球面径向基函数SRBF节点(中心)再通过自适应算法对比观测扰动场量的空间分布情况使得观测点空间分布与SRBF节点空间分布处处一致。Reuter格网等级Q用于表示球面等面积格网的空间分辨率其作用相当于全球重力位系数模型的最大阶数Nmaxn。Reuter格网又称为Reuter网格是一种用于在球面上特别是地球表面均匀分布点的特殊网格系统。它由Reuter在1982年提出在重力场建模尤其是球面径向基函数SRBF方法中常被用来预先确定基函数的数量和中心位置。1它是一个覆盖全球球面的点集。2点并非按等经纬度间隔排列那样会两极密集、赤道稀疏而是通过特殊算法调整每个纬度圈上的点数使得各点之间的球面距离尽可能均匀。3从极地看点的分布类似向日葵花盘或斐波那契螺旋的变体但从各个方向观测都呈现均匀性。4实际效果赤道附近点稍稀疏两极附近点稍密集因为球面面积在赤道更大但每个点所代表的“面积单元”大致相等。当观测点为规则格网时SRBF节点也呈现规则分布;当观测点分布不规则时SRBF节点分布也不规则测点分布密度大的地方SRBF节点分布密度也大。可见SRBF中心点空间分布与观测点空间分布处处一致。4、谱域SRBF法局部重力场逼近计算地球重力场逼近问题本质上是由一种类型观测重力场量求解地球外部空间其他各种类型重力场元的线性空间变换问题。可见移去常数因子GM/4π和相同变量因子1/r后不改变扰动场量之间的解析关系式(10.16)~式(10.21)因而可改写为由离散残差扰动重力mGal、高程异常m、空间异常mGal、扰动重力梯度径向E或垂线偏差南向/西向″中某种单一类型观测数据选择球面点质量核函数、Poisson核函数、m次径向多级子核或Poisson小波核函数中的一种球面径向基函数SRBF进行重力场逼近估计大地水准面外部残差扰动重力、残差高程异常、残差空间异常、残差扰动重力梯度或残差垂线偏差。通过输入不同单一类型观测场元、选择不同类型径向基函数及其参数、计算不同类型目标场元可全面测试径向基函数及其重力场逼近算法的空域、谱域和解析性质。通过将待评估目标场元的观测权置零或直接将待评估目标场元的测点作为计算点可有效探测场元观测量粗差测定其外部精度指标。程序本身还可用于各种单一类型观测扰动场元的格网化、解析延拓、类型转换、外部精度测定、重力场全要素建模与大地水准面计算。PAGravf4.5对航空、地面和海洋重力采用完全统一的处理方式航空扰动重力、海洋扰动重力和地面扰动重力属于同一类型的观测场元无需区分。5、入口函数GravityfieldapproachSRBFtest.f90输入离散残差观测扰动场元文件observationfl。头文件一行记录格式点号/点名 经度度小数 纬度度小数大地高…残差场元…权值…。输入计算面大地高格网文件surfhgtgrdfl。输入离散检核点位置文件checkpointfl。头文件一行记录格式点号/点名 经度度小数 纬度度小数大地高…。输入SRBF参数para(1:7)最小、最大阶数多级次数SRBF形式0-径向多级子核函数1-Poisson小波核函数Reuter格网等级SRBF中心作用距离kmBjerhammar球面埋藏深度km。输入观测场元和待估目标场元类型para(8:9)。输入观测文件记录中观测场元列序号para(10:11)当观测场元为垂线偏差时Para(10:11)为垂线偏差南、西方向的列序号否则para(11)缺省。输入观测文件记录中观测权属性列序号para(12)。▪ 当观测权属性列序号小于1或超出记录列序号或文件记录中观测权属性小于零时程序默认等权。▪ 当文件记录中观测权等于零时该观测量不参与SRBF系数估计程序可测定该观测量的外部精度指标。输入法方程解算方法para(13)1-LU分解2-Cholesky分解3-最小二乘QR分解4-最小范数奇异值分解5-岭估计法。6、程序编制程序编制完成以后一运行程序就闪退通过调试锁定了闪退出在函数Equsolve上。deepseek搜索了一下问题可能出在缺DLL运行库或CAD对内存大小有限制。于是编了一个纯VB Form程序调用函数结果VS提示缺少mkl_avx2.dll在执行程序路径补充此DLL程序运行成功。程序主界面如下图。各种结果文件与章老师程序一致。检核点结果文件高程异常相差0.0005m。程序在当前目录下输出观测量剩余残差点值文件。头文件格式:原观测量平均值标准差最小值最大值剩余残差平均值标准差最小值最大值。记录格式:点号经度纬度大地高权值原观测量(垂线偏差西向一仅当原观测量为垂线偏差时存在)残差量(残差垂线偏差西向)。可进一步用于原观测量粗差探测与精度评估。程序还在当前目录下输出多种场元球面径向基函数空域曲线文件*spc。头文件格式:球面径向基函数类型(0-径向多级子核函数1-Poisson小波核函数)次数最小阶数最大阶数补偿深度(km)。多种场元球面径向基函数谱域曲线文件*dgr。头文件与*spc相同。记录格式:阶数n扰动重力、高程异常、扰动重力梯度和总垂线偏差n阶归一化径向基函数值。多种场元球面径向基函数中心文件文件*center。头文件格式:Reuter格网等级RBF中心点数子午圈方向单元格网数平行圈方向最多单元格网数纬度间隔(’)。记录格式:点号经度(度小数)大地纬度单元格网面积差百分比平行圈方向单元格网经度间隔(’)参考文献局部重力场SRBF逼近及其性能指标测评Fortran代码 - 章传银PAGravf4.5参考说明书.pdf