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

CASM.jpg

区域地表负荷SRBF逼近法全要素负荷形变效应计算Fortran代码

SRBF逼近法负荷形变效应计算源码

 二维码 20
发表时间:2024-12-15 10:46作者:章传银网址:http://www.zcyphygeodesy.com
文章附图

[计算目标]

    由移去全球负荷球谐系数模型值的区域陆地水、海平面变化、江河湖库水、冰川雪山等地表环境负荷等效水高残差格网,按谱域球面径向基函数(SRBF)逼近与负荷效应综合算法,计算负荷残差等效水高逼近值(cm,用于质量评估),以及高程异常(mm)、地面重力(μGal)、扰动重力(μGal)、地倾斜(SW南向/西向mas)、垂线偏差(SW南向/西向mas)、水平位移(EN东向/北向mm)、地面径向(大地高mm)、地面正(常)高(mm)、扰动重力梯度(径向mE)与水平重力梯度(NE北向/西向mE)的残差负荷形变效应。

    程序要求地表负荷残差格网范围必须大于计算点分布范围,以吸收边缘效应。计算海平面变化负荷效应时,输入计算点的高度为正(常)高;计算陆地水负荷效应时,输入计算点相对于地面的高度。

幻灯片82.JPG幻灯片85.JPG

[地球物理模型]

    地球负荷勒夫数文件Love_load_cm.dat。采用球对称无旋转弹性地球模型REF6371计算的负荷勒夫数(来源于区域地面回弹计算器REAR1.0,2015.11)。

[测试入口程序]

    SurfaceEwhSRBFloadeffect.f90

    输入计算点空间位置文件calcpntfl,记录格式:点号/点名 经度(度小数) 纬度(度小数) 高度(相对于陆海面,m)......

    输入负荷等效水高残差(cm)格网文件ewhgridfl。

    输入高阶负荷勒夫数文件loadlovfl。

    输入累积逼近次数para(1)和法方程解算方法para(2)(1-LU分解,2-Cholesky分解,3-最小二乘QR分解,4,5-最小范数奇异值分解)。

    输入主SRBF参数para(3:9):SRBF形式(0-径向多级子核函数,1-Poisson小波核函数),多级次数,最小、最大阶数,Bjerhammar球面埋藏深度(km),Reuter等级,SRBF中心作用距离(km),Reuter格网等级。

    输入累积逼近SRBF参数para(10:16)。

[主要调用模块]

(1)残差负荷等效水高SRBF逼近法全要素负荷效应计算

    EwhSRBFLoadeffectpnt(calcpntfl,ewhgridfl,loadlovfl,para)

    离散计算点全要素残差负荷形变效应文件reslt.txt。头文件同输入计算点文件,文件记录在输入文件记录的基础上,增加残差负荷等效水高逼近值(cm,用于质量评估)和14列地表环境负荷效应的残差值,包括高程异常(mm)、地面重力(μGal)、扰动重力(μGal)、地倾斜(SW南向/西向mas)、垂线偏差(SW南向/西向mas)、水平位移(EN东向/北向mm)、地面径向(大地高mm)、地面正(常)高(mm)、扰动重力梯度(径向mE)与水平重力梯度(NE北向/西向mE)负荷效应的残差值。

    模块还在当前目录下输出全要素场元球面径向基函数空域曲线文件SRBFspc.rbf与全要素场元球面径向基函数谱域曲线文件SRBFdgr.rbf。

    SRBFspc.rbf头文件格式:球面径向基函数类型(0-径向多级子核函数,1-Poisson小波核函数),次数,最小阶数,最大阶数,补偿深度(km)。记录格式:球面距离(km),高程异常、地面重力、扰动重力、地倾斜、垂线偏差、水平位移、径向位移、正(常)高、扰动梯度和水平梯度归一化径向基函数值。

    SRBFdgr.rbf头文件与*spc.rbf相同。记录格式:阶数n,高程异常、地面重力、扰动重力、地倾斜、垂线偏差、水平位移、径向位移、正(常)高、扰动梯度和水平梯度归一化径向基函数值。

(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)全要素负荷效应SRBF曲线计算模块

    SRBF11all(RBF,flv,order,krbf,mpn,mdp,mp2,minN,maxN,NF,nta)

    返回全部11种要素的SRBF曲线RBF(NF+1, 11)。每种SRBF曲线用一维向量表示,向量长度等于NF+1,由SRBF中心作用距离和格网等级计算。

    RBF(NF+1,knd):knd=1等效水高,2高程异常,3地面重力,4扰动重力,5地倾斜,6垂线偏差,7水平位移,8大地高,9正常高,10扰动梯度,11水平梯度

    输入flv(:,3)-负荷勒夫数。

   输入SRBF宽度参数nta-nta=(r0-dpth)/r0,dpth为Bjerhammar球面补偿深度,r0为观测点平均地心距。

    输入krbf,order-krbf=0径向多极子,krbf-1Possion小波,order次数。

   输入mpn(maxN-minN+1, NF+1), mdp(maxN-minN+1, NF+1), mp2(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)大型正定对称方程组解算模块

    subroutine 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)正常重力场元计算模块

    normdjn(GRS,djn);GNormalfd(BLH,NFD,GRS)

    返回NFD(5)-正常重力位,正常重力,正常重力梯度,正常重力线方向,正常梯度方向。

(8)勒让德函数及其导数计算模块

    LegPn_dt2(pn,dp1,dp2,n,t)

    计算勒让德函数Pn(t)及其对ψ一、二阶导数t=cos ψ。

(9)大地坐标形式变换包

    BLH_RLAT(GRS, BLH, RLAT);BLH_XYZ(GRS, BLH, XYZ)

    RLAT_BLH(GRS, RLAT, BLH)

(10)其他辅助模块

    LegPn02(mpn,mdp,mp2,minN,maxN,NF,dr); PickRecord(str0, 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。

[算法公式]ETideLoad4.5参考说明书

    7.2.3大地测量数值格网文件

    8.7负荷SRBF逼近与负荷效应SRBF综合

附件rar压缩包:visual studio_intel fortran 集成环境测试项目、DOS可执行测试程序、地球物理模型和测试输入输出数据。