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

CASM.jpg

多源异质数据SRBF重力场全要素建模Fortran代码

多源异质SRBF重力场全要素建模Fortran代码

 二维码 26
发表时间:2025-01-09 11:58作者:章传银网址:http://www.zcyphygeodesy.com
文章附图

[计算目标]

    综合多种异质离散残差观测扰动场元(包括扰动重力、高程异常、空间异常、扰动重力梯度或垂线偏差),采用球面径向基函数逼近方法,估计给定计算面上全要素残差场元格网,从而实现重力场及大地水准面的统一建模。

    本程序是一种适应性超强的重力场全要素建模与质量测评万能工具,观测量无需归算、延拓及格网化。程序既能直接联合多源异质、不同高度、交叉分布、陆海共存的多种重力场观测量,实现区域重力场的全空间全要素高精度建模,又可有效应对复杂情形下各种观测量粗差探测、外部精度测定与计算性能控制系列难题。

幻灯片83.JPG幻灯片84.JPG

[测试入口程序]

    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可执行测试程序和测试输入输出数据。

幻灯片86.JPG幻灯片87.JPG

幻灯片88.JPG幻灯片89.JPG