局部重力场SRBF逼近及其性能指标测评Fortran代码

[计算目标]

    由离散残差扰动重力（mGal）、高程异常（m）、空间异常（mGal）、扰动重力梯度（径向，E）或垂线偏差（南向/西向，″）中某种单一类型观测数据，选择球面点质量核函数、Poisson核函数、m次径向多级子核或Poisson小波核函数中的一种球面径向基函数SRBF进行重力场逼近，估计大地水准面外部残差扰动重力、残差高程异常、残差空间异常、残差扰动重力梯度或残差垂线偏差。

    通过输入不同单一类型观测场元、选择不同类型径向基函数及其参数、计算不同类型目标场元，可全面测试径向基函数及其重力场逼近算法的空域、谱域和解析性质。通过将待评估目标场元的观测权置零，或直接将待评估目标场元的测点作为计算点，可有效探测场元观测量粗差，测定其外部精度指标。

    程序本身还可用于各种单一类型观测扰动场元的格网化、解析延拓、类型转换、外部精度测定、重力场全要素建模与大地水准面计算。

    PAGravf4.5对航空、地面和海洋重力采用完全统一的处理方式，航空扰动重力、海洋扰动重力和地面扰动重力属于同一类型的观测场元，无需区分。

[测试入口程序]

    GravityfieldapproachSRBFtest.f90

    输入离散残差观测扰动场元文件observationfl。头文件一行，记录格式：点号/点名 经度（度小数） 纬度（度小数）大地高，…，残差场元，…，权值，…。

    输入计算面大地高格网文件surfhgtgrdfl。

    输入离散检核点位置文件checkpointfl。头文件一行，记录格式：点号/点名 经度（度小数） 纬度（度小数）大地高，…。

    输入SRBF参数para(1:7)：最小、最大阶数，多级次数，SRBF形式（0-径向多级子核函数，1-Poisson小波核函数），Reuter格网等级，SRBF中心作用距离（km），Bjerhammar球面埋藏深度（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-岭估计法。

[主要调用模块]

（1）单一类型残差观测场元SRBF重力场逼近计算

    SRBFestimation(observationfl,surfhgtgrdfl,checkpointfl,para)

    输出与计算面大地高格网规格相同的目标残差扰动场元格网文件SRBFestimate.dat。

    模块还在当前目录下输出观测量剩余残差点值文件residuals.txt、检核点目标残差场元估计结果文件checkreslt.txt、多种场元球面径向基函数空域曲线文件SRBFspc.txt、谱域曲线文件SRBFdgr.txt和SRBF中心文件SRBFcenter.txt。

    residuals.txt头文件格式：原观测量平均值，标准差，最小值，最大值，剩余残差平均值，标准差，最小值，最大值。记录格式：点号，经度，纬度，大地高，权值，原观测量（，垂线偏差西向-仅当原观测量为垂线偏差时存在），残差量（，残差垂线偏差西向）。文件可进一步用于原观测量粗差探测与精度评估。

    checkreslt.txt头文件同checkpointfl，记录格式为在checkpointfl记录的基础上，增加一列目标残差场元估值。

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