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

CASM.jpg

等位面上逆Stokes_Hotine_Meinesz运算数值积分Fortran代码

逆Stokes_Hotine_Vening Meinesz运算

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

[计算目标]

    由等位面大地高格网(m)及其面上残差高程异常(m)或垂线偏差向量(S,W,″)格网,按逆Stokes、Hotine或Vening-Meinesz运算数值积分方法,计算等位面上其他类型残差扰动重力场元。

    扰动重力场元逆运算积分属Stokes边值问题,要求被积高程异常或垂线偏差向量(S,W)位于等位面,通常用于海洋卫星测高重力场反演计算。

    为实现有限半径积分,通常需采用参考重力场移去恢复法,先移去等位边界面上的模型空间异常/扰动重力,再积分得到计算点处的残差高程异常,最后恢复计算点处的模型高程异常。

    等位面可采用参考重力场模型(不大于360阶)构造,在高度不大于10千米的近地空间,可用等正(常)高面大地高格网表示。

幻灯片56.JPG幻灯片58.JPG

[测试入口程序]

    InvStokesHotineMeinesznumintg.f90

    运算类型参数knd-=0逆Stokes运算,=1逆Hotine运算,=2逆Vening-Meinesz运算。

    输入参数dr-积分半径(m)。

    输入参数calcpntfl-计算点文件名。头文件一行,记录格式:点号/点名 经度(度小数) 纬度(度小数)......。

    输入参数dwmhgrdfl-等位面大地高格网文件名。

    输入参数gravgrdfl-等位面上残差高程异常(knd=0或1)或垂线偏差向量(knd=2)格网文件名。

[主要调用模块]

(1)扰动重力场元逆运算数值积分算法模块

    InversenumIntegral(calcpntfl,dwmhgrdfl,gravgrdfl,knd,dr)

    输出文件reslt.txt。文件记录格式:

      当knd=0时,在输入文件记录的基础上,增加一列由等位面大地高格网内插得到的计算点大地高,和一列该点的残差空间异常积分值。

      当knd=1时,在输入文件记录的基础上,增加一列由等位面大地高格网内插得到的计算点大地高,和一列该点的残差扰动重力积分值。

      当knd=3时,在输入文件记录的基础上,增加一列由等位面大地高格网内插得到的计算点大地高,和3列该点的残差高程异常、残差扰动重力和残差空间异常积分值。

(2)扰动场元径向梯度数值积分算法模块

    real*8 function RaGradientBLH(BLH,gra,dwm,nlat,nlon,hd,dr,GRS)

    输入BLH(3)-计算点的经纬度(度小数)和大地高(m)。

    输入dwm(nlat,nlon)-等位边界面大地高格网,用于精确计算积分面元与计算点的积分距离。

    输入gra(nlat,nlon)-等位面上残差扰动场元格网。

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

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

    返回-残差扰动场元径向梯度积分值。单位等于残差扰动场元gra单位/m,方向为向径方向(由地心指向地球外部)。

(3)逆Vening-Meinesz运算数值积分算法模块

    AnVMdeflBLH(BLH,kai,nta,dwm,nlat,nlon,hd,dr,rst,GRS)

    输入kai(nlat,nlon)-等位面上残差垂线偏差南向分量格网(弧度)。

    输入nta(nlat,nlon)-等位面上残差垂线偏差西向分量格网(弧度)。

    返回rst(3)-计算点的残差高程异常(m)、残差空间异常(m/s2)和残差扰动重力(m/s2)。

(4)场元一阶水平梯度向量估计模块

    HGradientOne(BLH,rga,hgt,nlat,nlon,hd,tx,GRS)

    输入rga(nlat,nlon)-场元数值格网。

    输入hgt(nlat,nlon)-场元所在面的大地高格网。

    返回参数tx(2)-计算点BLH处场元rga的1阶水平梯度向量估值,tx(1)北向,tx(2)东向,单位等于场元gra单位/m。

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

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

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

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

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

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

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

    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);……

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

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

    RLAT_BLH(GRS, RLAT, BLH)

(9)其他辅助模块

    CGrdPntD2(lon,lat,dt,row,col,hd);PickRecord(line, kln, rec, nn)

[编译连接]

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

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

    1.4.1大地测量数据文件格式约定

    7.9.3扰动重力场元逆运算积分公式

    7.1(4)低阶勒让德函数及其一、二阶导数算法

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