多源异质数据SRBF重力场全要素建模Fortran代码多源异质SRBF重力场全要素建模Fortran代码 二维码
26
发表时间:2025-01-09 11:58网址:http://www.zcyphygeodesy.com ![]() [计算目标] 综合多种异质离散残差观测扰动场元(包括扰动重力、高程异常、空间异常、扰动重力梯度或垂线偏差),采用球面径向基函数逼近方法,估计给定计算面上全要素残差场元格网,从而实现重力场及大地水准面的统一建模。 本程序是一种适应性超强的重力场全要素建模与质量测评万能工具,观测量无需归算、延拓及格网化。程序既能直接联合多源异质、不同高度、交叉分布、陆海共存的多种重力场观测量,实现区域重力场的全空间全要素高精度建模,又可有效应对复杂情形下各种观测量粗差探测、外部精度测定与计算性能控制系列难题。 [测试入口程序] AllelementSRBFheterogeneous.f90 输入多种异质离散残差观测场元文件observationfl。头文件一行,记录格式:点号/站名,经度(度小数),纬度(度小数),大地高(m),残差观测场元,…,场元类型(0~5),权值,…。记录前5项属性的位置和顺序约定不变。 ▪ 观测场元类型=0扰动重力,=1高程异常,=2空间异常,=3扰动重力梯度,=4垂线偏差。 输入计算面大地高格网文件surfhgtgrdfl。 输入SRBF参数para(1:7):最小、最大阶数,多级次数,SRBF形式(0-径向多级子核函数,1-Poisson小波核函数),Reuter格网等级,SRBF中心作用距离(km),Bjerhammar球面埋藏深度(km)。 输入观测场元类型与权值所在的列序号para(8:9)。 ▪ 当观测权小于1,或超出记录列序号,或文件记录中观测权属性小于零时,程序默认等权。 ▪ 当文件记录中观测权等于零时,该观测量不参与SRBF系数估计,程序可测定该观测量的外部精度指标。 输入可调控场元类型和可调控观测场元贡献率para(10:11)。 输入法方程解算方法para(12):=1-LU分解,=2-Cholesky分解,=3-最小二乘QR分解,=4-最小范数奇异值分解,=5-岭估计法。 [主要调用模块] (1)多源异质SRBF全要素重力场建模 SRBFheterogeneous(observationfl,surfhgtgrdfl,para) 输出计算面上全要素残差重力场元文件SRBFhetero.txt。 ▪ 记录格式:点号,经度(度小数),纬度,计算点大地高(m),残差扰动重力(mGal),残差高程异常(ml),残差空间异常(mGal),残差扰动梯度径向(E),残差垂线偏差南向(″),残差垂线偏差西向(″)。 模块在当前目录下输出计算面上的残差扰动重力SRBFhetero.rga、残差高程异常SRBFhetero.ksi、残差空间异常SRBFhetero.gra、残差扰动重力梯度SRBFhetero.grr和残差垂线偏差向量SRBFhetero.dft格网文件。格网规格同计算面大地高格网。 模块还在当前目录下输出观测量剩余残差点值文件residuals.txt、多种场元球面径向基函数空域曲线文件SRBFspc.txt、谱域曲线文件SRBFdgr.txt和SRBF中心文件SRBFcenter.txt。 ▪ residuals.txt文件中每种类型观测量的统计结果占住一行头文件,格式:场元类型(0~5),原观测量平均值,标准差,最小值,最大值;残差平均值,标准差,最小值,最大值。记录格式:测点号,经度,纬度,测点大地高,残差量,原观测量,场元类型,权值。文件可进一步用于原观测量粗差探测与精度评估。 ▪ SRBFspc.txt头文件格式:球面径向基函数类型(0-径向多级子核函数,1-Poisson小波核函数),次数,最小阶数,最大阶数,补偿深度(km)。记录格式:球面距离(km),扰动重力、高程异常、扰动重力梯度和总垂线偏差归一化径向基函数值。 ▪ SRBFdgr.txt头文件与SRBFspc.txt相同。记录格式:阶数n,扰动重力、高程异常、扰动重力梯度和总垂线偏差n阶归一化径向基函数值。 ▪ SRBFcenter.txt头文件格式:Reuter格网等级,RBF中心点数,子午圈方向单元格网数,平行圈方向最多单元格网数,纬度间隔(′)。记录格式:点号,经度(度小数),大地纬度,单元格网面积差百分比,平行圈方向单元格网经度间隔(′)。 (2)Reuter格网参数计算模块 ReuterGrid(rhd,lvl,Kt,blat,nn,mm,nln,sr,dl,nrd,lon) 输入rhd(4)-用球坐标表示的Reuter格网范围(度小数),即最小、最大经度,最小、最大地心纬度。 输入lvl, nn, mm-Reuter格网等级,子午圈方向最多Reuter中心点数,平行圈方向最多中心点数。 返回参数Blat-第一行平行圈Reuter中心的地心纬度(度小数)。 返回参数Kt-Reuter中心点数量,等于未知数个数。 返回参数nln(nn)-平行圈方向Reuter中心的数量。 返回参数sr(nn) -平行圈方向Reuter单元格网面积与赤道单元格网面积之差的百分比。 返回参数dl(nn) -平行圈方向经度间隔(度小数)。 返回参数nrd(nn,mm) -Reuter中心的排序号。 返回参数lon(nn,mm) -Reuter中心经度。 (3)Reuter中心与观测点空间分布最佳匹配模块 Edgnode(enode,rlatlon,lvl,edgn,lon,blat,nln,gpnt,nn,mm) 模块计算Reuter单元格网中的观测点数,并修正Reuter格网中心点数Kt→gpnt。 返回参数edgn-Reuter格网四周边缘Reuter中心点的数量。 返回参数enode(edgn) -Reuter格网四周边缘Reuter中心点的序号。 返回参数rlatlon(edgn,2) -Reuter格网四周边缘Reuter中心点的球坐标(地心纬度, 经度)。 (4)5种类型场元SRBF曲线计算模块 SRBF5all(RBF,order,krbf,mpn,mdp,minN,maxN,NF,nta) 返回5种类型场元的SRBF曲线RBF(NF+1, 5)。每种SRBF曲线用一维向量表示,向量长度等于NF+1,由SRBF中心作用距离和格网等级计算。 RBF(NF+1,knd):knd=1扰动重力,=2高程异常,=3空间异常,=4重力梯度,=5总垂线偏差 输入SRBF宽度参数nta-nta=(r0-dpth)/r0,dpth为Bjerhammar球面补偿深度,r0为观测点平均地心距。 输入krbf,order-krbf=0径向多极子,krbf=1 Possion小波,order次数。 输入mpn(maxN-minN+1, NF+1), mdp(maxN-minN+1, NF+1)-全部minN~maxN阶勒让德函数及其一阶导数。 (5)计算点在Reuter格网中的位置计算模块 RtGridij(rln,ki,kj,blat,lvl,nn,mm,nln,dl,lon) 输入rln(3)-计算点的球坐标。 返回参数ki,kj-计算点rln在Reuter格网中的位置,用Reuter格网的序号数组元素表示,ki>0, kj>0。 (6)大型正定对称方程组解算模块 Equsolve(BB,xx,nn,BL,knd,bf) 调用mkl_lapack95_ilp64.lib解大型方程组BB.xx=BL,bf(8)为解的性质。 输入参数knd-=1 LU分解,=2 Cholesky分解,=3 最小二乘QR分解,=4 最小范数奇异值分解。 (7)正常重力场元计算模块 GNormalfd(BLH,NFD,GRS) 返回NFD(5)-正常重力位,正常重力,正常重力梯度,正常重力线方向,正常梯度方向。 (8)勒让德函数及其导数计算模块 LegPn_dt1(pn,dp1,n,t) 计算勒让德函数Pn(t)及其对ψ一阶导数t=cos ψ。 (9)大地坐标形式变换包 BLH_RLAT(GRS, BLH, RLAT);RLAT_BLH(GRS, RLAT, BLH) (10)其他辅助模块 LegPn01(mpn,mdp,minN,maxN,NF,dr); PickRecord(line, kln, rec, nn) RBFvalue(RBF(:,1),NF,dr,dln(2),tmp); drln(rln,rlnk,dln); Stat1d(dt,nn,rst) [编译连接] Fortran固定格式代码,Intel Fortran编译器,外部连接库mkl_lapack95_ilp64.lib。 [算法公式] PAGrav4.5说明书 7.10球面径向基函数重力场逼近理论与算法 附件zip压缩包:visual studio_intel fortran 集成环境测试项目、DOS可执行测试程序和测试输入输出数据。 |