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

CASM.jpg

多种监测量SRBF法陆地水及负荷形变场时序解算Fortran代码

SRBF法陆地水及负荷形变场时序解算代码

 二维码 39
发表时间:2024-12-21 21:04作者:章传银网址:http://www.zcyphygeodesy.com
文章附图

[计算目标]

    由各种大地测量监测量时间序列,按球面径向基函数SRBF负荷形变场逼近算法,计算地表负荷等效水高(cm),高程异常(大地水准面mm)、地面重力(μGal)、扰动重力(μGal)、地倾斜(mas)、垂线偏差(mas)、地面水平位移(mm)、大地高(mm)、正(常)高(mm)、扰动重力梯度(mE)与水平重力梯度(E)负荷形变效应格网时间序列。

    技术上要求,所有参与计算的监测量,应事先采用相同的地表负荷模型和相容的非潮汐效应(历元归算)算法,统一大地测量监测基准和参考历元。同时,移去负荷形变场的长波部分,以满足负荷形变场局部SRBF逼近条件。

    地面监测量可包括如下6种类型中的一种或多种:1(GNSS水准网点)高程异常变化mm,2(流动重力GNSS或固体潮CORS并置站点)扰动重力变化μGal,3(流动重力或固体潮站点)地面重力变化μGal,4(CORS站点或流动GNSS点)大地高变化mm,5(流动水准网点)正常高变化mm,6(水文监测站)等效水高变化cm。

幻灯片90.JPG

[测试入口程序]

    LoadeffectSRBFheterovariation.f90

    输入参数obstmsqdfl-地面站点监测量记录时序文件名。

    头文件含时序长度,以及按时序长度依次排列的采样历元时刻。记录格式:站点名称,经度,纬度,…,监测量权,监测量类型,…,按时序长度依次排列的时序采样值(缺省值9999.0000)。

    监测量类型 = 1高程异常变化mm,2扰动重力变化μGal,3地面重力变化μGal,4大地高变化mm,5正常高变化mm,6等效水高变化cm。

    输入参数dtmfl-计算面高度格网文件名。

    计算面高度为计算点相对于地面的高度。若计算地面负荷形变场时,输入零值格网。计算面高度格网规格用于指定待估陆地水格网的经纬度范围和空间分辨率。

    输入参数lnth,itern-输入站点平均间距(m)和累计SRBF逼近次数。

    输入参数cntrlknd,kwgh-设置可调控监测量类型(1~6)及其贡献率。

    输入参数knd-法方程解算方法。=1 LU分解,=2 Cholesky分解,=3 最小二乘QR分解,=4 最小范数奇异值分解。

    输入参数hepch,frow-记录时序头文件首个采样时刻列序号,记录时序文件记录中首次采样列序号。

    输入参数kndrow,wghrow-监测量类型和权值列序号。

    输入参数krbf, krbf1,order,order1-主SRBF和累积SRBF1的类型krbf,krbf1及次数。

    输入参数dr, dr1, dpth, dpth1-主SRBF和累积SRBF1中心的作用距离(m)和Bjerhammar球面埋藏深度(m)。

    输入参数minN,maxN, minN,1maxN1-主SRBF和累积SRBF1的勒让德函数最小最大阶数类型。

[主要调用模块]

(1)SRBF法陆地水及负荷形变场解算模块

    LoadestmateSRBF(dtmfl,obstmsqdfl,flv,dtrow,inp)

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

    输入参数dtrow-监测量记录时序中当前历元监测量列序号。

    输出文件-每个历元时间12个文件,包括陆地水负荷等效水高ewh***.dat与监测量残差文件rnt***.txt,以及如下10种负荷形变场格网文件。

      ①高程异常变化SRBFgeoid***.dat

      ②地面重力变化SRBFterrgrav***.dat

      ③扰动重力变化SRBFgravdist***.dat

      ④地倾斜变化向量SRBFgrndtilt***.dat

      ⑤垂线偏差变化向量SRBFvertdefl***.dat

      ⑥水平位移向量SRBFhorzdisp***.dat

      ⑦大地高变化SRBFelliphgt***.dat

      ⑧正常高变化SRBForthohgt***.dat

      ⑨径向重力梯度变化SRBFgradient***.dat

      ⑩水平重力梯度变化向量SRBFhorzgrad***.dat。

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

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

    SRBFdgr.rbf头文件与SRBFspc.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=1 Possion小波,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)正常重力场元计算模块

    GNormalfd(BLH,NFD,GRS)

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

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

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

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

(8)大型法方程解算模块库

    EqHestens(BPB,xx,nn,BPL);EqJordon(BPB,xx,nn,BPL)

    EqCholesky(BPB,xx,nn,BPL);EqueSVD(BPB,xx,nn,BPL)

    RidgeEstimate(BPB,xx,nn,BPL);Equsolve(BPB,xx,nn,BPL,knd,bf)

(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); PickReclong(line, kln, rec, nn)

    RBFvalue(RBF(:,1),NF,dr,dln(2),tmp); drln(rln,rlnk,dln); Stat1d(dt,nn,rst)

[编译连接]

    Fortran固定格式代码,任何fortran编译器,外部连接库mkl_lapack95_ilp64.lib。

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

    7.3大地测量时间序列文件约定格式

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

附件zip压缩包:visual studio_intel fortran 集成环境测试项目、DOS可执行测试程序、全部测试输入输出数据。