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

CASM.jpg

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

格林积分法陆地水及负荷形变时序解算代码

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

[计算目标]

    由各种大地测量监测量时间序列,以负荷格林函数积分为动力学约束,计算地表负荷等效水高(cm),高程异常(大地水准面mm)、地面重力(μGal)、扰动重力(μGal)、地倾斜(mas)、垂线偏差(mas)、地面水平位移(mm)、大地高(mm)、正(常)高(mm)、扰动重力梯度(mE)与水平重力梯度(E)负荷形变效应格网时间序列。

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

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

幻灯片87.JPG

[测试入口程序]

    GeodeicGreenestmLoadeffect.f90

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

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

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

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

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

    输入参数lnth,dr-输入站点平均间距(m)和格林函数积分半径(m)。

    输入参数lvb,itern-设置Laplace平滑算子权值与累计逼近次数。

    输入参数nfar-设置边缘效应抑制参数。

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

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

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

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

    输入参数kndrow-监测量类型列序号。

    输入参数wghrow-监测量权值列序号。

[主要调用模块]

(1)格林积分法陆地水及负荷形变场解算模块

    Loadeffectestmgreen(dtmfl,obstmsqdfl,GF,dtrow,inp)

    输入参数GF(8000,9)-高程异常(10-13)、地面重力10-17、扰动重力10-18、地倾斜10-14、垂线偏差10-19、水平位移10-12、径向位移10-11、径向重力梯度10-15与水平重力梯度10-15的负荷间接效应格林函数值,GF(i,1:9)的作用距离为100i(m)。

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

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

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

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

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

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

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

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

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

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

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

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

    其中,***为历元时间。

(2)监测量格林积分约束观测方程构造模块

    CoefLoadobs(BLH,hd,BB,nx,GF,dr,GRS,kobs,sk,kk)

    输入BLH(3)-监测站点的经纬度(度小数)和高度(m)。

    输入dr, hd(6)-格林积分半径(m)和格网规格参数(最小最大经度,最小最大纬度,经度间隔,纬度间隔)。

    输入nx=nlat×nlon-待估未知参数个数。

    输入GRS(6)-gm, ae, j2, omega, 1/f, 缺省值。

    输入kobs-监测量类型(1~5)。

    返回BB(nx)-观测方程系数向量。

    返回sk(kk)-kk为BB(nx)<nx中不为零的个数,sk(kk)为不为零的系数序号向量。

(3)残差负荷效应格林积分法计算模块

    rntGreenintegral(BLH,ewh,hd,nlat,nlon,GF,direct,indrct,GRS,dr)

    输入ewh (nlat,nlon)-地表环境负荷等效水高格网(cm)。

    返回参数direct(10) -计算点处高程异常(mm)、地面重力(μGal)、扰动重力(μGal)、地倾斜(SW南向/西向mas)、垂线偏差(SW南向/西向mas)、扰动重力梯度(径向mE)与水平重力梯度(NE北向/西向mE)的残差负荷直接效应。

    返回参数indrct(14) -计算点处高程异常(mm)、地面重力(μGal)、扰动重力(μGal)、地倾斜(SW南向/西向mas)、垂线偏差(SW南向/西向mas)、水平位移(EN东向/北向mm)、地面径向(大地高mm)、地面正(常)高(mm)、扰动重力梯度(径向mE)与水平重力梯度(NE北向/西向mE)的残差负荷间接效应。其中,地面正(常)高的残差负荷间接效应indrct(11)=0。

(4)Laplace空间平滑约束参数方程构造模块

    Laplacesmth(BPB,s,nlat,nlon,lvb)

    输入参数lvb-Laplace平滑算子权值。

    输入参数s(nlat*nlon) -提供nlat*nlon实数向量空间。

    输入输出参数BPB(nlat*nlon,nlat*nlon) -输入法方程系数阵,输出增加空间平滑约束后的法方程系数阵。

(5)格网边缘远区零约束参数方程构造模块

    Farzonexx(BPB,nlat,nlon,nfar,wgh)

    输入参数nfar-设置边缘效应抑制参数。

    输入参数wgh-将格网边缘待估参数对应的法方程对角线+wgh。

    输入输出参数BPB(nlat*nlon,nlat*nlon)-输入法方程系数阵,输出格网边缘远区零约束后的法方程系数阵。

(6)读取负荷间接影响格林函数

    LGrnFunc(loadgrfl,GF)

    输入loadgrfl-负荷间接效应格林函数文件名LoadGreen.txt。

    返回GF(8000,9)-高程异常(10-13)、地面重力10-17、扰动重力10-18、地倾斜10-14、垂线偏差10-19、水平位移10-12、径向位移10-11、径向重力梯度10-15与水平重力梯度10-15的负荷间接效应格林函数值,GF(i,1:9)的作用距离为100i(m)。

(7)正常重力场元计算模块

    GNormalfd(BLH,NFD,GRS)

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

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

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

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

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

    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)

(10)其他辅助模块

    PickReclong(line, kln, rec, nn);Stat1d(dt,nn,rst);Gauss2D(lon,lat,dt,row,col,hd)

[编译连接]

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

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