西安测绘研究所王云鹏：DQM2022系列超高阶地球重力场模型构建及其精度评估|《测绘学报》2024年53卷第8期

摘要：本文基于椭球谐分析，提出了超高阶地球重力场模型构建的局部积分改进迭代法和全球积分改进迭代法，解决了传统的局部积分改进法存在的局限性，有效提高了改进模型在中国地区的适用精度。选择EGM2008、EIGEN-6C4地球重力场模型作为初始模型，利用中国及周边地区最新的5'×5'实测格网平均重力异常数据，构建了DQM2022系列超高阶地球重力场模型，其完全阶次均为2190。采用地面实测重力异常、GNSS/水准、天文大地垂线偏差等数据对改进模型进行精度评估。结果显示：①相对于初始模型，改进模型表示中国地区重力场的精度显著提高，重力异常精度提高了2.4~2.8mGal，高程异常精度提高了1.0~2.4cm，垂线偏差精度提高了0.07″~0.15″;②在表示中国地区重力场时，基于EIGEN-6C4初始模型构建的改进模型精度最高，高程异常精度约为10.6cm，垂线偏差精度约为2.1″。关键词：局部积分改进迭代法;全球积分改进迭代法;超高阶地球重力场模型;椭球谐分析;重力异常;高程异常;垂线偏差

基金项目

国家自然科学基金(42174001)(41774018);地理信息工程国家重点实验室开放基金(SKLGIE2023-Z-1-1)

王云鹏（1986—），男，硕士，助理研究员，研究方向为物理大地测量。E-mail：wangyp1813@163.com通信作者:刘晓刚E-mail：liuxuyanchu2022@163.com

本文引用格式王云鹏,刘晓刚,李琦,等.DQM2022系列超高阶地球重力场模型构建及其精度评估[J].测绘学报,2024,53(8):1505-1516.doi:10.11947/,LIUXiaogang,LIQi,etal.Constructionofseriesultra-high-degreeEarth'sgravityfieldmodelsDQM2022andtheirprecisionevaluation[J].ActaGeodaeticaetCartographicaSinica,2024,53(8):1505-1516.doi:10.11947/

阅读全文

地球重力场模型是表示地球重力场的一类基本参数的集合[1-4]，通常采用球谐或椭球谐系数来表示，可以方便地计算重力异常、高程异常及垂线偏差等重力场元。随着地球重力场模型阶数和精度的不断提高，其在大地测量学、地球物理学等领域的应用前景越来越广阔。

我国从20世纪70年代初开始开展重力场模型构建的理论研究和模型构制，先后发表过多个系列的地球重力场模型，其中，有代表性的是西安测绘研究所构制的DQM系列模型、武汉测绘科技大学构制的WDM系列模型和中国科学院测量与地球物理研究所构制的IGG系列模型。DQM系列模型先后发表了DQM77、DQM84、DQM94、DQM99和DQM2000等[5-8]，WDM系列模型先后发表了WDM89、WDM92CH和WDM94等[9-10]，IGG系列模型先后发表了IGG71、IGG93和IGG97等[11]。进入21世纪以后，国内学者注重利用卫星重力数据研制低阶重力场模型，如UCAS_Grace01、Tongji-GOGR2019S、WHU-GM-05、WHU-Grace01s、WHU-SWPU-GOGR2022S等[12-16]。目前，国内学者也构建了一些其他的超高阶地球重力场模型，如武汉大学的UGM08、SGG-UGM-1和SGG-UGM-2模型[17-19]，信息工程大学的UGM05、IEU_UGM_2022模型[20-21]，这些模型的主要数据源是现有重力场模型、卫星重力数据和国外实测重力数据，使用的国内重力数据较陈旧，且数据不全面。上述重力场模型的构建方法有最小二乘配置法、联合平差法及局部积分改进法(属于剪切法)等。其中，局部积分改进法是最经济、最便捷的高阶重力场模型构建方法，这种方法已被证明在局部计算中可以获得令人满意的结果。然而，传统的局部积分改进法存在两个问题：①为了减少计算量、降低计算复杂度，传统公式中只考虑了k=0相关项的影响(k是不大于(n-m)/2的自然数，n、m为模型的阶和次)，而忽略了k0相关项的影响[8,19]，导致解算的位系数存在一定误差；②局部积分迭代法对地球重力场只进行了一次逼近，文献[1]研究表明迭代计算能够进一步提高地球重力场模型的精度。

为解决上述问题，本文提出超高阶地球重力场模型构建的局部积分改进迭代法和全球积分改进迭代法(属于剪切法范畴)。推导严密的基于椭球谐分析的剪切法计算公式，避免因计算公式简化引起的误差，并引入迭代计算思想，进一步提高改进模型的精度。采用中国及周边重点区域实测重力数据，构建DQM2022系列超高阶地球重力场模型，并对模型的适用性进行了分析和评价。

本文选择EGM2008、EIGEN-6C4重力场模型作为初始模型，基于椭球谐分析，采用剪切法改进初始模型。EGM2008重力场模型构建时采用地面重力、航空重力、卫星测高及GRACE卫星等数据，在中国区域，使用美国国家地理空间情报局收集的中南部地区重力数据[22]。EIGEN-6C4重力场模型不仅使用了上述数据，还使用了GOCE卫星观测数据[23]。

改进模型的椭球谐系数计算公式为

(1)

式中，n、m为模型的阶和次；和为改进模型的完全正规化椭球谐系数(已减去正常场位系数)；和为初始模型的完全正规化椭球谐系数(已减去正常场位系数)，利用Jekeli系数转换公式计算得到[24-25]；和为椭球谐系数改正数；Pk为权因子，采用文献[5]提出的赋权思想计算权因子，考虑到初始模型的中低阶位系数(主要由卫星重力数据计算得到)精度已足够高，因此，当n≤260时，不改正模型位系数。椭球谐系数改正数的计算公式为

其中

式中，i表示椭球面上第i个纬度带；j表示椭球面上第j个子午带；Nδ和Nλ分别为积分区域内归化纬度方向和经度方向的格网个数；δ为归化纬度的余角；λ为大地经度；a为参考椭球的长半轴；γi为椭球面上的正常重力；为n阶平滑因子[26]；为椭球面格网中心点的地心向径；(cosδ)为完全正规化缔合勒让德函数；是由初始模型计算的椭球面格网平均重力异常；是由实测格网平均重力异常经过各项改正后得到的椭球面格网平均重力异常，计算公式为

式中，为实测格网平均重力异常；为大气改正、、和椭球改正（为实测重力异常向量归算到椭球法线方向的改正，、为椭球法线方向与地心径向之间夹角影响的改正为以格网中心点地心向径代替格网内任意点地心向径的改正);为正常重力二阶梯度改正；为延拓改正。这些改正项的数学模型在文献[27-29]中已经给出，本文重点介绍椭球改正和延拓改正的计算方法，其他改正项不再赘述。、、和的计算公式为

(4)

其中

式中，为重力场球谐模型的截断阶数，取值为2190；为格网中心点的地心向径；θi为格网中心点的地心余纬；、为完全正规化球谐系数(已减去正常场位系数)；(cosθ)为完全正规化缔合勒让德函数；θ为地心余纬；GM为重力场模型的引力常数；e为第一偏心率；J2为2阶0次项正常场位系数；ω为地球自转角速度；、为椭球主曲率半径；hij为格网中心点的大地高。

本文采用泰勒展开的方法计算延拓改正[22]。利用椭球面重力异常，将实测重力异常表示为级数式

(5)

则延拓改正可表示为

(6)

利用球谐展开式表示

(7)

忽略地心径向和椭球法向的差异，对式(7)求k阶导数，代入式(6)，得到

(8)

其中

式中，Kmax为k阶导数的最高阶数，取值为10[22]。

DQM2022系列超高阶地球重力场模型的构建过程如下。

(1)进行坐标和重力系统转换，将实测格网平均重力异常的参考系统归算到初始模型的参考系统。

(2)以初始模型为参考模型，利用式(3)，对实测格网平均重力异常进行各项改正，得到椭球面格网平均重力异常。

(3)利用参考模型，计算椭球面格网平均重力异常，计算公式[22]为

(9)

其中

式中，为重力场椭球谐模型的截断阶数，取值为2160；和为完全正规化椭球谐系数(已减去正常场位系数);、和的定义同式(2)。

(4)基于残余重力异常，利用式(2)计算椭球谐系数改正数，并通过式(1)得到改进模型的椭球谐系数。

(5)以改进模型为参考模型，重复步骤(3)和步骤(4)，直到残余重力异常的均方根趋于稳定为止，一般迭代5~10次即可。

(6)以步骤(5)得到的改进模型作为参考模型，重复步骤(2)—步骤(5)，一般迭代3~5次即可。

(7)迭代结束后，利用Jekeli系数转换公式[24-25]将改进模型的椭球谐系数(2160阶)转换为球谐系数(2190阶)。

本文提出的上述计算步骤包含内外两层循环，内循环主要目的是提高改进模型对重力异常数据的逼近程度，外循环主要目的是提高椭球改正和延拓改正的精度，内外循环相结合，进一步提高了改进模型的精度。延拓改正、椭球改正、位系数解算中涉及超高阶勒让德函数计算、局部或全球积分等密集计算问题，本文采用OpenMP并行计算技术提高计算效率，要求计算设备具备多线程处理能力；模型位系数、勒让德函数、临时数组等大型矩阵的内存消耗较多，而且并行计算使临时数组的存储需求成倍增加，因此需要内存较大的计算设备。本文在模型构建时采用的计算机线程为16个，内存为32GB，局部积分改进迭代法耗时约30min，全球积分改进迭代法耗时约2h。计算过程中，当积分区域为局部地区时，本文称为局部积分改进迭代法；当积分区域为全球时，本文称为全球积分改进迭代法，两种方法均属于剪切法。

本文基于中国及周边地区5'×5'格网平均重力异常数据，采用局部积分改进迭代法和全球积分改进迭代法(全球积分时，因缺乏全球重力实测数据，空白区用初始模型仿真数据填充)，构建了DQM2022系列超高阶地球重力场模型，命名为DQM2022RA、DQM2022RB、DQM2022GA和DQM2022GB(其中，R表示局部积分，G表示全球积分，A表示初始模型为EGM2008，B表示初始模型为EIGEN-6C4)，完全阶次均为2190。改进模型与初始模型的参考椭球参数一致，即DQM2022RA和DQM2022GA的引力常数为3.986004415×1014m3·s-2，长半轴为6378136.3m;DQM2022RB和DQM2022GB的引力常数为3.986004415×1014m3·s-2，长半轴为6378136.46m。

本文分别从内部符合和外部检核两个方面评估DQM2022系列超高阶地球重力场模型的精度。

阶方差体现了重力场模型不同频段的能量分布，反映了重力场模型的能量分布是否与Kaula准则[30]在同一水平。其表达式[31]为

(10)

式中，和为改进模型的球谐系数。

阶误差RMS反映了改进模型与初始模型的符合程度。其表达式[31]为

(11)

其中

式中，和为初始模型的球谐系数。

重力异常和大地水准面的阶次误差和累计误差能够直观地反映模型位系数误差对重力场的影响量级有多大，是评估重力场模型精度最常用的统计量。本文以初始模型EGM2008、EIGEN-6C4为标准模型，利用DQM2022系列超高阶地球重力场模型计算重力异常和大地水准面的阶次误差和累计误差，计算公式[31]为

(12)

式中，R为地球平均半径；和分别为重力异常阶次误差和累计误差；和分别为大地水准面阶次误差和累计误差。

若令为模型构建中使用的5'×5'实测格网平均重力异常，与重力场模型计算的平均重力异常构成差值(本文称为重力异常误差)，其表达式为

(13)

其中

若令ζt为GNSS/水准点的高程异常，与重力场模型计算的高程异常ζmod构成差值(本文称为高程异常误差)，其表达式为

(14)

其中

式中，ζ0为高程异常的零阶项；γ为计算点的正常重力。

若令ξt和ηt为由天文经纬度与大地经纬度导出的天文点实测垂线偏差，与重力场模型计算的垂线偏差ξmod和ηmod构成差值(本文称为垂线偏差误差)，其表达式为

(15)

其中

式中，ξmod为计算点垂线偏差子午分量；ηmod为计算点垂线偏差卯酉分量；ρ为弧度转化为角秒的比例因子(ρ≈206264.806″)；(cosθ)/∂θ为完全正规化缔合勒让德函数的一阶导数[31]。

本文采用的局部重力数据是中国及周边地区的重力异常观测值，包含国内733571个地面重力数据、19668695个船载重力数据、93000条航空重力数据及卫星测高重力数据。基于该数据，在中国及周边地区共构成256479个5'×5'格网平均重力异常。处理过程为：首先，采用均值漂移粗差探测、位系数模型检查及DEM高程检查等方法[32-33]进行粗差探测和剔除，得到纯净的自由空间重力异常；其次，进行中间层改正、空间改正和地形改正得到离散点布格重力异常[27]；然后，采用反距离加权方法内插得到格网中心点的布格重力异常，并以此作为格网平均值；最后，将布格重力异常恢复成自由空间重力异常。格网平均重力异常及高程(正常高)的统计结果见表1。

表1中国及周边地区5'×5'格网平均重力异常及高程(正常高)

5′×5′gridmeangravityanomalyandelevation(normalheight)inChinaanditssurroundingareas

新窗口打开|下载CSV

(1)阶方差和阶误差RMS。利用式(10)计算DQM2022系列超高阶地球重力场模型的阶方差，结果如图1所示。利用式(11)计算DQM2022系列超高阶地球重力场模型的阶误差RMS，结果如图2所示。

图1DQM2022系列超高阶地球重力场模型的阶方差

̓sgravityfieldmodels

图2DQM2022系列超高阶地球重力场模型的阶误差RMS

̓sgravityfieldmodels

由图1可知，DQM2022系列超高阶地球重力场模型的阶方差曲线与Kaula曲线很接近，说明改进模型位系数的量级是正确的。图2阶误差RMS曲线的左端点横坐标位于260阶，这是因为在模型构建过程中，当n≤260时，未改正模型位系数。由图2可知：①DQM2022系列超高阶地球重力场模型的阶误差RMS较小，说明该系列模型与初始模型的符合程度较高；②DQM2022GA模型的阶误差RMS小于DQM2022RA模型，DQM2022GB模型的阶误差RMS小于DQM2022RB模型，说明全球积分模型(DQM2022GA和DQM2022GB)与初始模型的符合程度更好。

(2)阶次误差和累计误差。利用式(12)计算重力异常和大地水准面的阶次误差和累计误差，结果如图3、图4所示。

图3重力异常阶次误差和累计误差

图4大地水准面阶次误差和累计误差

由图3、图4可知：①DQM2022GA和DQM2022GB模型误差曲线位于DQM2022RA和DQM2022RB模型误差曲线下方，说明DQM2022GA和DQM2022GB模型位系数误差对重力异常和大地水准面的影响小于DQM2022RA和DQM2022RB模型；②DQM2022RA和DQM2022RB模型的重力异常累计误差约为1.5mGal，大地水准面累计误差约为2.2cm，DQM2022GA和DQM2022GB模型的重力异常累计误差约为1mGal，大地水准面累计误差约为1.5cm。全球积分模型误差小于局部积分模型的原因是空白区用初始模型仿真数据填充，避免了边界效应。

(3)重力异常的比较。利用式(13)对DQM2022系列超高阶地球重力场模型计算的重力异常数据与模型构建使用的实测数据进行比较，统计结果见表2。

表2中国及周边地区重力异常误差

GravityanomalyerrorinChinaandsurroundingareas

新窗口打开|下载CSV

由表2可知，以标准差为评价指标，在表示中国及周边地区重力异常时，相对于EGM2008模型，DQM2022RA和DQM2022GA模型精度分别提高了2.6和2.4mGal；相对于EIGEN-6C4模型，DQM2022RB和DQM2022GB模型精度分别提高了2.8和2.6mGal。DQM2022系列超高阶地球重力场模型在表示中国及周边地区重力异常时精度相当。

(1)高程异常的比较。本文收集了中国地区实测的一、二、三等GNSS及水准资料，共6227个点，这是目前国内分布最广泛、最权威的GNSS及水准资料，分布情况如图5所示。

注：台湾资料暂缺。

图5中国地区实测GNSS及水准点分布

本文按照文献[34]的统计方式，分5个区域分别进行对比分析。利用式(14)对DQM2022系列超高阶地球重力场模型计算的高程异常数据与实测数据进行比较，统计结果见表3。

表3中国地区高程异常误差

ElevationanomalyerrorinChina

新窗口打开|下载CSV

由表3可知，以标准差为评价指标[35]，在表示中国地区高程异常时：①相对于EGM2008模型，DQM2022RA和DQM2022GA模型在5个区域内的精度均提高了约1cm;②相对于EIGEN-6C4模型，DQM2022RB和DQM2022GB模型在区域2和区域3内的精度提高了约1cm，在其他区域内提高了约3cm，在全国提高了约2cm;③相比其他重力场模型，DQM2022RB和DQM2022GB模型精度最高，在全国的精度为10.6cm。在表示中国地区高程异常时，局部积分模型和全球积分模型精度相当。

(2)垂线偏差的比较。本文收集了中国天文大地锁网的884个天文点，并导出垂线偏差分量ξt和ηt，天文点分布情况如图6所示。

注：台湾资料暂缺。

图6中国地区实测垂线偏差点位分布

利用式(15)对DQM2022系列超高阶地球重力场模型计算的垂线偏差数据与实测数据进行比较，统计结果见表4。

表4中国地区垂线偏差误差

VerticaldeflectionerrorinChina

新窗口打开|下载CSV

由表4可知，以标准差为评价指标，在中国地区：①表示垂线偏差子午分量ξ时，相对于EGM2008模型，DQM2022RA和DQM2022GA模型精度分别提高了0.09″和0.12″，相对于EIGEN-6C4模型，DQM2022RB和DQM2022GB模型精度分别提高了0.07″和0.14″;②表示垂线偏差卯酉分量η时，相对于EGM2008模型，DQM2022RA和DQM2022GA模型精度分别提高了0.14″和0.15″，相对于EIGEN-6C4模型，DQM2022RB和DQM2022GB模型精度分别提高了0.14″和0.15″。在表示中国地区垂线偏差时，全球积分模型精度略高于局部积分模型。

本文基于椭球谐分析，推导了严密的剪切法计算公式，提出了超高阶重力场模型构建的局部积分改进迭代法和全球积分改进迭代法，利用我国最新的5'×5'实测格网平均重力异常数据，构建了DQM2022系列超高阶地球重力场模型，并评估了模型精度，得到以下结论。

(1)相对于EGM2008模型，在中国地区，DQM2022RA模型表示重力异常、高程异常、垂线偏差子午分量和卯酉分量的精度分别提高了约2.6mGal、1.0cm、0.09″、0.14″;DQM2022GA模型表示重力异常、高程异常、垂线偏差子午分量和卯酉分量的精度分别提高了约2.4mGal、1.2cm、0.12″、0.15″。

(2)相对于EIGEN-6C4模型，在中国地区，DQM2022RB模型表示重力异常、高程异常、垂线偏差子午分量和卯酉分量的精度分别提高了约2.8mGal、2.4cm、0.07″、0.14″;DQM2022GB模型表示重力异常、高程异常、垂线偏差子午分量和卯酉分量的精度分别提高了约2.6mGal、2.4cm、0.14″、0.15″。

(3)局部积分改进迭代法存在边界效应的影响，采用局部积分方式得到的地球重力场模型局部适用性较好，而全球积分改进迭代法避免了边界效应，得到的重力场模型在全球的精度更优。

(4)在掌握全球重力数据的条件下，本文提出的全球积分改进迭代法，是一种基于全球已有重力场模型，进而构建全球超高阶重力场模型的最经济、最便捷的方法。下一步，将收集、整理全球1'×1'格网平均重力异常数据，采用全球积分改进迭代法研制我国10800阶DQM系列超高阶地球重力场模型。

初审：侯琳

复审：宋启凡

终审：金君