快速检索        
  武汉大学学报·信息科学版  2018, Vol. 43 Issue (5): 643-650

文章信息

黄谟涛, 刘敏, 马越原, 欧阳永忠, 邓凯亮, 陆秀平, 吴太旗
HUANG Motao, LIU Min, MA Yueyuan, OUYANG Yongzhong, DENG Kailiang, LU Xiuping, WU Taiqi
海域超大范围外部重力场快速赋值模型构建
Model Construction for Computing External Gravity Field in Ultra-Large Sea Area
武汉大学学报·信息科学版, 2018, 43(5): 643-650
Geomatics and Information Science of Wuhan University, 2018, 43(5): 643-650
http://dx.doi.org/10.13203/j.whugis20160351

文章历史

收稿日期: 2016-11-25
海域超大范围外部重力场快速赋值模型构建
黄谟涛1,2,3 , 刘敏2 , 马越原1,2 , 欧阳永忠1 , 邓凯亮1 , 陆秀平1 , 吴太旗1     
1. 海军海洋测绘研究所, 天津, 300061;
2. 信息工程大学地理空间信息学院, 河南 郑州, 450001;
3. 海军工程大学导航工程系, 湖北 武汉, 430033
摘要:针对海域超大范围外部重力场快速赋值的特殊需求,分析了3种传统扰动引力赋值模型的适用性和局限性,分别提出了直接积分改进模型、点质量改进模型和直接积分与点质量混合模型。利用数值计算验证了3组改进模型的合理性和有效性,为实际工程应用奠定了必要的理论基础。
关键词外部重力场     扰动引力     模型构建     超大区域    
Model Construction for Computing External Gravity Field in Ultra-Large Sea Area
HUANG Motao1,2,3, LIU Min2, MA Yueyuan1,2, OUYANG Yongzhong1, DENG Kailiang1, LU Xiuping1, WU Taiqi1     
1. Naval Institute of Hydrographic Surveying and Charting, Tianjin 300061, China;
2. Institute of Geospacial Information, Information Engineering University, Zhengzhou 450001, China;
3. Department of Navigation, Naval University of Engineering, Wuhan 430033, China
First author: HUANG Motao, professor, specializes in marine gravity field. E-mail: Ouyangyz@sohu.com
Corresponding author: DENG Kailiang, PhD. E-mail:dengkailiang036@163.com
Foundation support: The National Natural Science Foundation of China, Nos.41474012, 41374018;the Great Scientific Instrument Development Project of China, No.2011YQ12004503;the National Basic Research Program of China, No.613219; the National Key R & D Program of China, Nos.2016YFC0303007, 2016YFB0501704
Abstract: In accordance with the special requirement of determining the external gravity field quickly in ultra-large sea area, a detail analysis is made on the applicability and limitations of the three traditional models, i.e. spherical harmonic expansion model, direct integral model and point mass model, for computing disturbing gravity. And then three modified computational models, i.e. modified direct integral model, modified point mass model and mixed model of direct integral and point mass, are proposed one by one. A numerical test has been made to prove the reasonableness and validity of the suggested models. It is shown that the singularity in three traditional models have been eliminated in the solutions of the modified models. They can satisfy the practical needs of determining the local gravity field quickly in ultra-large sea areas and full height. The accuracy of the modified point mass model is estimated to be about ±1mGal in smooth sea areas, and to be less than ±3mGal in rugged sea areas. The obtained conclusions are valuable to the practical engineering projects.
Key words: external gravity field     disturbing gravity     model construction     ultra-large area    

地球外部重力场赋值是局部重力场逼近计算的主要研究内容之一,也是大地测量边值问题解算的最终目的[1]。这项研究在空间科学领域具有重要的应用价值[2-3]。在地球外部空间飞行的人造卫星和各类飞行器由于受到地球扰动引力场的作用,其飞行轨道必定会偏离正常的椭圆轨道,即产生所谓的轨道摄动。为了精准确定飞行轨道摄动,必须精确计算飞行器在飞行过程中所受的地球扰动引力矢量,这就是开展外部重力场赋值研究的目的所在[3-4]

当前应用于外部重力场赋值的计算模型主要有3类:(1)基于全球重力位模型的球谐展开式[2-5];(2)基于经典Stokes边值理论的直接积分模型[2, 6];(3)基于现代Bjerhammar换置理论的点质量模型[3, 7-9]。已有研究结果表明,3类模型各具不同的技术特征和适用条件,球谐展开模型主要用于扰动重力的长波部分计算[3-4],但受计算时间的限制,球谐展开式的阶次不宜取得太高;积分模型不需要对地面重力数据作进一步的转换处理,但其在超低空存在奇异性和无法顾及地形效应的缺陷,在一定程度上限制了该模型的应用[3, 9];点质量模型对超低空扰动重力场具有较强的恢复能力,但赋值之前需要完成地面重力数据到点质量的转换,其中涉及较大规模线性方程组的解算问题[3-4, 9]。文献[8]最早将点质量模型引入我国重力场逼近计算研究领域;文献[10]根据内部场源与外部场的频谱响应关系,提出了扰动重力频谱响应质点模型的构建方法。此后有诸多学者陆续开展点质量解算方法和扰动引力快速赋值模式研究[11-14],但已有研究成果主要针对局部区域固定点模式的外部重力场赋值,几乎不涉及海域超大范围流动点模式的外部重力场快速赋值问题。远程飞行器对重力场的保障需求主要体现为扰动引力矢量计算精度和速度两个方面,在确保精度指标要求的前提下,如何突破扰动引力的计算速度指标是该研究领域当前面临的主要挑战。飞行器海域发射阵位具有高度机动性和覆盖范围广的显著特点,因此,与陆地固定发射阵位相比较,远程飞行器对海域外部空间重力场赋值的技术要求更高,特别是数据准备阶段的实时性要求更高。针对这种特定的需求,本文在分析现有赋值模式技术特点的基础上,探讨联合应用多种模式构建一种能够兼备各自模式优良特性的综合模型的可能性,同时从工程化应用角度出发,提出传统赋值模式的改进方案,并对改进效果作出量化评价,为远程飞行器重力场保障提供技术支撑。

1 传统赋值模型及适用性分析 1.1 球谐展开模型

计算地球外部扰动引力三分量的球谐展开式为[2]

$ \begin{array}{*{20}{c}} {{\rm{ \mathsf{ δ} }}{g_r} = - \frac{{GM}}{{r_p^2}}\sum\limits_{n = 2}^N {\left( {n + 1} \right){{\left( {\frac{a}{{{r_p}}}} \right)}^n}} \sum\limits_{m = 0}^n {\left[ {\bar C_{nm}^ * \cos m\lambda + } \right.} }\\ {\left. {{{\bar S}_{nm}}\sin m\lambda } \right]{{\bar P}_{nm}}\left( {\sin \varphi } \right)} \end{array} $ (1)
$ \begin{array}{*{20}{c}} {{\rm{ \mathsf{ δ} }}{g_\varphi } = \frac{{GM}}{{r_p^2}}\sum\limits_{n = 2}^N {{{\left( {\frac{a}{{{r_p}}}} \right)}^n}} \sum\limits_{m = 0}^n {\left[ {\bar C_{nm}^ * \cos m\lambda + } \right.} }\\ {\left. {{{\bar S}_{nm}}\sin m\lambda } \right]\frac{{{\rm{d}}{{\bar P}_{nm}}\left( {\sin \varphi } \right)}}{{{\rm{d}}\varphi }}} \end{array} $ (2)
$ \begin{array}{*{20}{c}} {{\rm{ \mathsf{ δ} }}{g_\lambda } = - \frac{{GM}}{{r_p^2\cos \varphi }}\sum\limits_{n = 2}^N {{{\left( {\frac{a}{{{r_p}}}} \right)}^n}} \sum\limits_{m = 0}^n {m\left[ {\bar C_{nm}^ * \sin m\lambda + } \right.} }\\ {\left. {{{\bar S}_{nm}}\cos m\lambda } \right]{{\bar P}_{nm}}\left( {\sin \varphi } \right)} \end{array} $ (3)

式中,rpφλ为计算点的地心向径、地心纬度和地心经度;GM为地球引力常数;a为地球椭球长半径;Cnm*Snm为完全正常化位系数;Pnm(sinφ)为完全正常化勒让德函数;N为展开式的最高阶次。在实际应用中,球谐展开模型首先用于远程飞行器主动段飞行轨迹扰动引力长波段的计算,然后用于被动段飞行轨迹扰动引力全频段的计算。因此,它是扰动引力计算模型必不可少的组成部分。但如前所述,受计算时间的限制,球谐展开式的阶次不宜取得太高,一般取为N=22或N=36[9]

1.2 直接积分模型

根据文献[2-3],可将扰动引力三分量的数值积分计算式表示为:

$ {\rm{ \mathsf{ δ} }}{g_r} = \frac{R}{{4{\rm{ \mathsf{ π} }}}}\sum\limits_{i = 1}^{{N_1}} {\sum\limits_{j = 1}^{{N_2}} {\Delta {g_{ij}}{{\left[ {\frac{{\partial S\left( {r,\psi } \right)}}{{\partial r}}} \right]}_{{p_{ij}}}}\Delta {\sigma _{ij}}} } $ (4)
$ {\rm{ \mathsf{ δ} }}{g_\varphi } = - \frac{R}{{4{\rm{ \mathsf{ π} }}{r_p}}}\sum\limits_{i = 1}^{{N_1}} {\sum\limits_{j = 1}^{{N_2}} {\Delta {g_{ij}}{{\left[ {\frac{{\partial S\left( {r,\psi } \right)}}{{\partial \psi }}} \right]}_{{p_{ij}}}}\cos {\alpha _{{p_{ij}}}}\Delta {\sigma _{ij}}} } $ (5)
$ {\rm{ \mathsf{ δ} }}{g_\lambda } = - \frac{R}{{4{\rm{ \mathsf{ π} }}{r_p}}}\sum\limits_{i = 1}^{{N_1}} {\sum\limits_{j = 1}^{{N_2}} {\Delta {g_{ij}}{{\left[ {\frac{{\partial S\left( {r,\psi } \right)}}{{\partial \psi }}} \right]}_{{p_{ij}}}}\sin {\alpha _{{p_{ij}}}}\Delta {\sigma _{ij}}} } $ (6)

式中,p为计算点;Δgij代表地面网格平均重力异常;N1N2分别代表纬线和经线方向的数据网格数;ψ为计算点与流动点之间的地心角距;α为计算点与流动点之间的方位角;其他参量的具体含义及其计算式参见文献[2-3]。

在实际应用中,直接积分模型除了需要与球谐展开模型联合使用外,为了提高计算效率,赋值前通常需要对地面重力网格数据进行分级处理,具体做法是:以计算点为中心,将整个数据覆盖区按流动点到计算点由近至远划分成若干个子区,距离计算点越近的子区,数据网格间隔越小,即数据分辨率越高。目前常用的数据网格组合依次为2′×2′、5′×5′、20′×20′和1°×1°。此外,为了减弱边缘效应的影响,一般采用重力异常的逐级余差替代原有的数据网格平均值进行赋值计算,即用δΔg(2′)=Δg(2′)-Δg(5′)替代Δg(2′),用δΔg(5′)=Δg(5′)-Δg(20′)替代Δg(5′),用δΔg(20′)=Δg(20′)-Δg(1°)替代Δg(20′)。此时,数据层之间存在覆盖关系。直接积分赋值模型来源于经典的Stokes边值理论体系,具有较强的理论支撑,数据准备过程比较简便,但由于该模型的核函数存在较强的奇异性,同时无法顾及地形效应的影响,因此不适用于超低空扰动引力场的赋值[3, 9]

1.3 点质量模型

由虚拟点质量计算扰动引力的表达式为[3]

$ {\rm{ \mathsf{ δ} }}{g_r} = - \sum\limits_{i = 1}^{{N_1}} {\sum\limits_{j = 1}^{{N_2}} {\frac{{{r_p} - {R_{ij}}\cos {\psi _{{p_{ij}}}}}}{{l_{{p_{ij}}}^3}}{m_{ij}}} } $ (7)
$ {\rm{ \mathsf{ δ} }}{g_\varphi } = \sum\limits_{i = 1}^{{N_1}} {\sum\limits_{j = 1}^{{N_2}} {\frac{{{R_{ij}}\sin {\psi _{{p_{ij}}}}}}{{l_{{p_{ij}}}^3}}\cos {\alpha _{{p_{ij}}}} \cdot {m_{ij}}} } $ (8)
$ {\rm{ \mathsf{ δ} }}{g_\lambda } = \sum\limits_{i = 1}^{{N_1}} {\sum\limits_{j = 1}^{{N_2}} {\frac{{{R_{ij}}\sin {\psi _{{p_{ij}}}}}}{{l_{{p_{ij}}}^3}}\sin {\alpha _{{p_{ij}}}} \cdot {m_{ij}}} } $ (9)

式中,N1×N2代表点质量的个数;Rij=R-Hij代表相应点质量层所在球面(Bjerhammar球)的半径;R为地球平均半径;Hij为点质量埋藏深度;mij代表埋藏在Bjerhammar球上的网格点质量,其值由相对应的地面网格重力异常数据转换得到,具体解算方程式为:

$ \Delta {g_q} = \sum\limits_{i = 1}^{{N_1}} {\sum\limits_{j = 1}^{{N_2}} {\left( {\frac{{{r_q} - {R_{ij}}\cos {\psi _{{q_{ij}}}}}}{{l_{{q_{ij}}}^3}} - \frac{2}{{{r_q}{l_{{q_{ij}}}}}}} \right){m_{ij}}} } $ (10)

同样,在实际应用中,点质量模型也需要与球谐展开模型联合使用,解算点质量时也需要采用与前面相似的逐级余差方法,具体步骤为:

1) 首先,依据式(10)由Δg(1°)计算点质量M1(1°), 由M1(1°)分别计算Δg(20′)M1、Δg(5′)M1和Δg(2′)M1,求余差δΔg(20′)=Δg(20′)-Δg(20′)M1

2) 依据式(10),由δΔg(20′)计算点质量M2(20′),由M2(20′)分别计算Δg(5′)M2和Δg(2′)M2,求余差δΔg(5′)=Δg(5′)-Δg(5′)M1g(5′)M2

3) 依据式(10),由δΔg(5′)计算点质量M3(5′),由M3(5′)计算Δg(2′)M3,求余差δΔg(2′)=Δg(2′)-Δg(2′)M1g(2′)M2g(2′)M3,由δΔg(2′)计算点质量M4(2′)。

M1(1°)、M2(20′)、M3(5′)、M4(2′)共4层点质量和式(7)~(9)就共同组成完整的点质量赋值模型,此时,点质量层之间存在覆盖关系。该模型除了具有模型结构简单、不存在奇异性和能够精确顾及地形效应等优点外,在计算效果上,这种模式巧妙地将重力异常和点质量的置换过程与函数插值自然结合,从而减弱了数值离散化的影响,因此对超低空外部扰动引力场具有较强的恢复能力[3]。点质量模型的主要缺陷是:在数据准备阶段,需要花费较多的时间来完成多层点质量的解算,不利于实时性保障要求较高场合的应用[11-14]

2 海域发射阵位保障需求与模型改进

由文献[4, 15]可知,考虑到计算精度和速度两个方面的要求,人们一般采用全球模型和局部模型的组合,作为远程飞行器主动段飞行轨迹的重力场赋值模型。主动飞行段虽然时间很短(一般仅为200~300 s),但由于这个阶段的飞行高度较低(一般不超过300 km)[16],因此,对局部重力场赋值模型的构建速度和计算精度要求较高,特别是超低空飞行阶段。如前所述,与陆地固定发射阵位相比较,海域发射阵位的主要特点是高度机动性和覆盖范围广阔(比如需要覆盖整个太平洋和印度洋)。这样的特性一方面要求重力场保障必须有超大范围的基础数据做保证,另一方面要求必须具备实时或准实时化的保障能力,具体体现为数据准备和模型赋值两个阶段的快速保障能力。很显然,要想满足这样的特殊需求,必须对现有的重力场保障模式作必要的改进。除了作为全球模型的球谐展开式是地球外部扰动重力场赋值模型不可或缺的组成部分以外,如何针对直接积分模型和点质量模型的技术特点,通过模型结构优化和多模型联合应用等手段,构造出符合海域发射阵位保障需求的局部重力场快速赋值模型,是本文接下来需要研究解决的核心问题。

2.1 直接积分模型改进

由前面的论述可知,局部重力场赋值模型主要由以发射阵位为中心、一定范围内的基础数据和计算模型自身两个部分组成。直接积分模型使用的基础数据是不同分辨率、对应不同区域大小的重力异常网格平均值。根据当前的技术发展水平和应用需求[3-4],实际应用中一般将基础数据分辨率和相应的覆盖范围确定为:1°×1°→30°×30°;20′×20′→9°×9°;5′×5′→4°×4°;2′×2′→2°×2°。很显然,由于地球表面重力观测成果的基本形式就是不同分辨率的重力异常网格值,因此,无论是陆域还是海域发射阵位,只要一旦确认了发射阵位的具体位置,那么就可以在很短的时间内完成直接积分模型所需基础数据集的构建,不存在数据准备时间上的制约。但是,如前所述,直接积分模型在计算精度方面存在较大的缺陷:一是无法顾及地形效应的影响,二是在超低空高度段存在积分奇异性问题。对于海域发射阵位,因为陆地部分的数据块一般都远离计算点,只会引起较微小的远区地形效应影响,故前者可以忽略不计。关于积分奇异性问题,可以通过下面的途径加以克服。

首先利用已知的地(海)面网格平均重力异常Δg计算相对应的高程异常ζ和垂线偏差两分量(ξ, η)(实际上,远程飞行器重力场保障也同时需要这两组数据),并按式(11)~(13)计算0 m高度上的扰引力三分量[2-3]

$ {\rm{ \mathsf{ δ} }}g_r^0 = - \Delta g - \frac{{2\gamma }}{R}\zeta $ (11)
$ {\rm{ \mathsf{ δ} }}g_\varphi ^0 = - \gamma \cdot \xi $ (12)
$ {\rm{ \mathsf{ δ} }}g_\lambda ^0 = - \gamma \cdot \eta $ (13)

式中,γ为地球平均正常重力。在这个环节,虽然在计算高程异常和垂线偏差时同样会遇到积分奇异性问题,但由于ζ和(ξ, η)的奇异性要比计算扰动引力特别是径向分量的奇异性弱化许多,使用常规的处理方法就能获得比较稳定可靠的计算结果[3, 17],因此这一转换过程是可行有效的。然后,依据式(4)~(6)计算某个固定高度(一般取为h1=3 km或5 km)上的扰动引力三分量(δgr1, δgφ1, δgλ1),由于此时的计算点在地球外部一定高度之上,积分奇异性已经明显减弱,故其计算结果也应当是稳定可靠的。最后,可利用(δgr0, δgφ0, δgλ0)和(δgr1, δgφ1, δgλ1)两组数据,按内插公式(14)~(16)计算得到0 m和h1高度之间任意点hp的扰动引力三分量:

$ {\rm{ \mathsf{ δ} }}g_r^h = {\rm{ \mathsf{ δ} }}g_r^0 + \frac{{{h_p}}}{{{h_1}}}\left( {{\rm{ \mathsf{ δ} }}g_r^1 - {\rm{ \mathsf{ δ} }}g_r^0} \right) $ (14)
$ {\rm{ \mathsf{ δ} }}g_\varphi ^h = {\rm{ \mathsf{ δ} }}g_\varphi ^0 + \frac{{{h_p}}}{{{h_1}}}\left( {{\rm{ \mathsf{ δ} }}g_\varphi ^1 - {\rm{ \mathsf{ δ} }}g_\varphi ^0} \right) $ (15)
$ {\rm{ \mathsf{ δ} }}g_\lambda ^h = {\rm{ \mathsf{ δ} }}g_\lambda ^0 + \frac{{{h_p}}}{{{h_1}}}\left( {{\rm{ \mathsf{ δ} }}g_\lambda ^1 - {\rm{ \mathsf{ δ} }}g_\lambda ^0} \right) $ (16)

以上计算流程既可确保在超低空飞行段也能获得较高精度的扰动引力赋值结果,又可在一定程度上缩减扰动引力的计算时间,因此它是一种比较实用的海域发射阵位重力场保障赋值模式。本文将其称为模型1。

2.2 点质量模型改进

由式(7)~(9)可知,点质量模型使用的基础数据是不同分辨率和不同区域大小的点质量网格值,它们由相对应的重力异常网格值转换而来。这意味着,与前面的直接积分模型相比较,后者在数据准备阶段额外增加了多层点质量数据集的解算过程,这对于机动性要求很强、实时性保障要求很高的海域发射阵位来说,无论如何都是一种不可忽视的时间负担。为了解决此问题,陆域固定发射阵位一般将点质量解算工作提前到数据准备阶段之前完成,并将预先计算得到的点质量和重力异常一同作为基础数据存入诸元准备计算机,这样可大大缩短后端的数据准备时间。但对于海域发射阵位来说,要想将超大范围的重力异常数据一次性严密转换为相对应的点质量数据集,几乎是不可能的,特别是对于分辨率较高的数据层。一方面是因为超广的覆盖范围和超高的分辨率首先意味着要面临超大的计算量,另一方面是因为点质量解算过程同属于数学结构上存在先天性弱点的向下延拓不适定反问题[18],其数值解算过程容易出现失稳是客观存在的,数据覆盖范围越广、分辨率越高,数值解算过程出现失稳的几率也越高。此外,在远程飞行器发射保障中,点质量模型是作为局部重力场模型使用的。这就意味着,当完成了超大范围的多层点质量参数解算后,还需要以发射阵位为中心,按规定的区域大小挑选对应范围内的点质量作为赋值模型的基础数据,这就是所谓的“整体解算与局部应用”矛盾统一体问题。由式(10)知,由重力异常到点质量的整体转换过程已经大大增强了相邻网格点质量之间的关联性,因此,在超大范围的点质量数据集中挑选出一部分点质量,组成局部点质量模型用于外部扰动引力场的赋值,必将引起较大的计算边缘效应[3, 19-20]。针对上述问题,本文依据文献[21]的研究思路,提出了基于“移动窗口”控制的点质量二次迭代解法。

首先将解算点质量方程式(10)右端的系数改写为:

$ {a_{{q_{ij}}}} = \left\{ \begin{array}{l} \frac{{{r_q} - {R_{ij}}\cos {\psi _{{q_{ij}}}}}}{{l_{{q_{ij}}}^3}} - \frac{2}{{{r_q}{l_{{q_{ij}}}}}},\left| {{\varphi _q} - {\varphi _{ij}}} \right| \le \rho \;且\\ \;\;\;\;\;\;\;\;\left| {{\lambda _q} - {\lambda _{ij}}} \right| \le \rho \\ 0,\left| {{\varphi _q} - {\varphi _{ij}}} \right| > \rho \;或\;\left| {{\lambda _q} - {\lambda _{ij}}} \right| > \rho \end{array} \right. $ (17)

式中,ρ表示移动窗口半径。式(17)的物理意义为:人为控制虚拟点质量对地(海)面观测量的作用范围,其效果相当于数值逼近中的“局部基”作用[22]。可利用超松弛迭代法求解由式(17)定义的“移动窗口”控制方法形成的带状稀疏矩阵方程[20],获得点质量的初解。此后,为了达到既强调了局部又不失整体性这一目的,可考虑进一步将窗口以外的点质量对地(海)面重力异常的作用再次转化为窗口以内的点质量修正量。此时,可适当扩大移动窗口半径ρ,按式(17)确定系数矩阵,然后由式(10)计算地面重力异常的预测值,将预测值与实测值的较差作为新的观测量,按照第一次解算过程相同的步骤求解点质量的改正数。实际计算时,可以重复上述迭代过程, 直至重力异常较差小于某个规定的数值(比如0.1 mGal)为止,一般情况下迭代2~3次即可满足要求。第一次点质量解算结果和以后各次改正数的迭加,就构成点质量的最终解算结果。

不难看出,上述对点质量模型的改化过程可望取得3个方面的效果[3, 20-21]:(1)破解海域超大范围点质量参数的解算难题;(2)减弱相距较远网格点质量之间的相关性;(3)有效提高点质量求解过程的数值稳定性。需要补充说明的是,考虑到实际应用保障条件的限制,建议采用如下差异化的点质量解算策略:在基础数据准备的前期,预先完成海域超大范围1°×1°、20′×20′和5′×5′这3层点质量的解算;在基础数据准备阶段,实时完成局部区域2′×2′点质量层的解算,即只实时解算以发射阵位为中心、2°×2°范围内的2′×2′点质量。实验结果表明,上述解算策略完全满足海域实际应用对诸元计算速度的保障需求。本文将由此构成的点质量改进模型称为模型2。关于求解点质量时的窗口半径ρ1和计算地面重力异常预测值时的窗口半径ρ2设计问题,对于不同网格大小的点质量,应有不同的设置方案,本文推荐的参数组合为:(1)1°×1°网格,ρ1=8°,ρ2=20°;(2)20′×20′网格,ρ1=3°,ρ2=6°;(3)5′×5′网格,ρ1=1°,ρ2=3°。

2.3 直接积分和点质量混合模型

针对直接积分和点质量模型各自存在的缺陷,本文在前面已经分别提出了相应的改进策略,并构建了相对应的改进模型。但由前面的分析可知,直接积分模型的不足之处正是点质量模型的优势所在,即点质量模型可有效地顾及地形效应的影响,同时不存在积分奇异性干扰。反过来,后者的不足也正是前者的优势所在,即直接积分模型不需要对基础数据作进一步的转换,可节省可观的数据准备时间。由此可见,直接积分和点质量模型之间是一种优势互补关系。因此,本文提出以分频段赋值的方式构建直接积分和点质量的混合模型:(1)使用直接积分模型完成对应于1°×1°、20′×20′和5′×5′这3种网格数据频段的扰动引力赋值;(2)使用点质量模型完成对应于2′×2′网格数据频段的扰动引力赋值。

上述混合模型虽然也需要实时解算局部区域的2′×2′网格点质量,但不难看出,该模型已经融合了两类传统模型的优点,同时有效弥补了各自的不足。需要补充说明的是,混合模型使用的1°×1°、20′×20′、5′×5′和2′×2′这4层重力异常,都应当只是扣除了22阶或36阶位模型参考场后的基础数据,不再是前面统一使用的重力异常余差,即4层重力异常之间不再存在逐级覆盖关系。因为如果仍然使用重力异常余差作为基础数据,那么直接积分模型赋值部分还会遇到奇异性问题。上述变化对飞行器实际应用保障也会带来一定的影响,因为为了节省扰动引力计算时间,通常采用这样的策略:当计算高度增大到某个规定的量值以后,可考虑从高分辨率到低分辨率,逐步略去各层数据块的影响,最终只在被动飞行段保留低阶位模型参考场的作用。显然,混合模型采用的数据结构无法满足“逐步略去各层数据块影响”的实用性要求。为此,在软件实现时,可采取相应的改进策略:当一旦需要略去上一层高分辨率数据块的影响时,下一层数据块积分区域必须恢复涵盖上一层积分区域。因此计算高度已经增大到一定的量值,故不会出现奇异积分问题。本文将上述混合计算模型称为模型3。

3 数值计算检验与分析

扰动引力赋值模型计算精度主要涉及模型逼近误差和数据传播误差两个方面的问题[3, 23],数值离散化误差、远区截断误差和数据截断误差都属于模型逼近误差的范畴。为了定量估计模型逼近误差的大小,本文采用2 160阶次的全球位模型EGM2008作为标准场开展数值计算检验,由其模拟产生外部扰动引力三分量的真值,同时计算1°×1°、20′×20′、5′×5′和2′×2′这4组地(海)面重力异常的观测量,作为扰动引力赋值模型的基础数据。实验区涵盖了重力场变化比较剧烈的马里亚纳海沟,具体数据覆盖范围为:φ为10°S~60°N;λ为90°E~190°E,参考场的阶次取为N=36,模型1中的固定计算高度取为h1=5 km,4层点质量的埋藏深度依次取为H(1°)=60 km,H(20′)=20 km,H(5′)=5 km,H(2′)=2 km。表 1列出了4组地面重力异常观测量的统计结果,重力异常分布情况如图 1所示。

表 1 地面重力异常统计结果/ (10-5m·s-2) Table 1 Statistic of Anomalies from EGM2008/(10-5m·s-2)
数据块 最小值 最大值 平均值 均方根
2′×2′ -343.4 559.3 -0.05 37.9
5′×5′ -341.6 495.9 -0.05 37.6
20′×20′ -313.7 279.1 -0.04 34.5
1°×1° -234.8 161.7 -0.04 26.4
图 1 重力异常等值线图 Figure 1 Contour Map of Anomalies from EGM2008

考虑到边缘效应的影响,按照实际应用相同的要求将计算区域确定为:φ为5°S~55°N;λ为95°E~185°E,相当于数据覆盖范围向内各缩减5°。在上述范围内分别使用前面提出的模型1至模型3,均匀计算2°×2°网格点不同高度上的扰动引力三分量,同时使用EGM2008位模型计算相应高度上的标准值。将3组模型计算值分别与标准值作比较,结果表明,90%以上的互差值不超过1 mGal(1 mGal=10-5m/s2),互差超过1 mGal的计算点均出现在重力场变化比较剧烈的海沟和海山区域或出现在计算区域的边缘部分。作为两种情形的代表,表 2表 3分别列出了位于海沟区的P1(φ=35°N, λ=143°)和位于计算区边缘的P2(φ=35°N, λ=185°)两个计算点的比对结果。

表 2 海沟区P1计算点比对结果/ (10-5 m·s-2) Table 2 Comparison of the Calculated Results at Point P1 in Oceanic Trench/(10-5 m·s-2)
高度/km 0 0.02 1 3 5 10 30 50 100 300
模型1 Δδgr 0.02 -0.00 -0.38 -1.01 -1.56 -0.22 -0.17 -0.19 -0.17 -0.06
Δδgφ -0.28 -0.36 -1.41 -3.63 -6.03 -5.73 -4.66 -3.77 -2.16 -0.23
Δδgλ 0.45 0.42 0.33 0.30 0.47 0.46 0.40 0.31 0.09 -0.17
模型2 Δδgr -0.01 -0.02 -0.18 -0.27 -0.35 -0.46 -0.36 -0.18 0.08 0.15
Δδgφ -3.31 -3.30 -3.03 -2.68 -2.44 -2.03 -1.21 -0.81 -0.37 -0.02
Δδgλ 0.29 0.29 0.27 0.18 0.07 -0.15 -0.45 -0.46 -0.38 -0.32
模型3 Δδgr -1.37 -1.33 -0.84 -0.89 -0.94 -1.05 -1.36 -1.53 -1.52 -0.71
Δδgφ -3.98 -3.97 -3.86 -3.77 -3.69 -3.49 -2.83 -2.32 -1.40 -0.15
Δδgλ 1.89 1.89 1.83 1.83 1.82 1.80 1.63 1.36 0.69 -0.10
表 3 边缘区P2计算点比对结果/ (10-5m·s-2) Table 3 Comparison of the Calculated Results at Point P2 in Border Area/(10-5m·s-2)
高度/km 0 0.02 1 3 5 10 30 50 100 300
模型1 Δδgr -0.01 0.00 0.09 0.19 0.10 0.09 0.05 0.03 0.02 0.10
Δδgφ -0.20 -0.22 -0.50 -1.30 -2.31 -2.16 -1.70 -1.38 -0.91 -0.43
Δδgλ -0.01 -0.02 -0.05 -0.05 -0.15 -0.14 -0.12 -0.10 -0.04 -0.01
模型2 Δδgr -1.64 -1.61 -0.93 -0.53 -0.42 -0.32 -0.14 -0.09 -0.05 -0.06
Δδgφ -2.53 -2.53 -2.37 -2.15 -2.00 -1.71 -1.07 -0.69 -0.14 0.38
Δδgλ 0.21 0.21 0.18 0.16 0.15 0.13 0.10 0.09 0.06 0.02
模型3 Δδgr -2.47 -2.43 -1.56 -1.24 -1.12 -0.86 -0.30 -0.10 0.01 0.10
Δδgφ -3.67 -3.67 -3.60 -3.41 -3.23 -2.89 -2.13 -1.69 -1.06 -0.44
Δδgλ 0.50 0.50 0.47 0.42 0.38 0.31 0.15 0.08 0.02 -0.00

表 2表 3结果首先可以看出,3组改进模型均有效解决了传统积分模型固有的奇异性问题。相比较而言,在超低空高度段,模型1的计算效果要比模型2和模型3略好一些。这是因为在计算高程异常和垂线偏差时使用了较高阶次的参考场(本文N=360,实际应用时,由于高程异常和垂线偏差是事先完成计算的,故它们使用的参考场阶次可以不同于计算扰动引力时的取值),同时由于ζ和(ξ, η)积分核函数的奇异性程度相对较低,使得采用式(11)~(13)计算0 m高度上的扰动引力精度较高。但从表 2结果也看出,在重力场变化比较剧烈的海沟区,由于受离散化误差的影响,模型1在中低空高度段的计算精度仍远低于3 mGal,这个结果在一定程度上限制了该模型的应用范围。此外,从表 2表 3不难看出,在绝大多数情况下,模型3的计算误差要普遍大于模型2的结果,究其原因,这可能与该模型使用的直接积分和点质量混合模型存在一定的不相容性有关。因为尽管埋藏于地球内部的2′×2′网格点质量能够精确恢复地(海)面的2′×2′网格重力异常,但不能保证2′×2′网格和5′×5′网格重力异常数据之间的平稳过渡和对接,从而可能在一定程度上加大了外部扰动引力计算边缘效应的影响。综合表 2表 3的结果可知,从整体上讲,模型2的计算精度相对稳定可控,具有较好的适用性。计算效率方面,在数据实时准备阶段,模型1需要的时间相对较短,模型2和模型3需要的时间几乎相同;在扰动引力赋值阶段,模型2需要的时间相对较短,模型1和模型3需要的时间几乎相同。需要指出的是,虽然3组模型在两个不同阶段所需时间略有差异,但其量值大小都在工程应用可接受的范围之内,因此在当前技术条件下,该因素基本不影响3组模型在实际应用中的选择。

4 结语

针对海域超大范围外部重力场快速赋值需求,本文提出了3组计算外部扰动引力的改进模型,通过理论分析和数值计算检验,可得出以下结论。

1) 3组改进计算模型均有效克服了传统积分模型固有的奇异性问题,能够满足超大区域和全高度段对局部扰动引力场赋值的要求,具有较好的应用价值。

2) 相比较而言,点质量改进模型具有较为稳定可控的计算精度,在重力场变化比较平缓的海域,该模型的逼近误差(不包括数据传播误差)一般不超过1 mGal;在重力场变化比较剧烈的海域,该模型的计算误差一般可控制在3 mGal以内。

3) 在数据前期准备、实时处理和扰动引力实时赋值3个不同阶段,3组改进计算模型具有不同的技术特点和要求。在实际应用中,应根据用户的具体需求和保障条件,综合考虑数据准备难度、计算精度和计算效率等多种因素,对3组备选模型作出合理和适当的选择。

参考文献
[1] Moritz H. The Computation of the External Gravity Field and the Geodetic Boundary-Value Problem[M]//Orlin H. Extension of Gravity Anomalies to Unsurveyed Areas, Geophysical Monograph Series. Washington D C, USA: American Geophysical Union, 1966
[2] Moritz H. Physical Geodesy[M]. Wien, Austia: Springer, 2005
[3] Huang Motao, Zhai Guojun, Guan Zheng, et al. The Determination and Application of Marine Gravity Field[M]. Beijing: Surveying and Mapping Press, 2005 ( 黄谟涛, 翟国君, 管铮, 等. 海洋重力场测定及其应用[M]. 北京: 测绘出版社, 2005 )
[4] Lu Zhonglian, Wu Xiaoping, Ding Xingbin, et al. Gravimetry in Ballistic Missile[M]. Beijing: Bayi Press, 1993 ( 陆仲连, 吴晓平, 丁行斌, 等. 弹道导弹重力学[M]. 北京: 八一出版社, 1993 )
[5] Wang Jianqiang, Li Jiancheng, Wang Zhengtao, et al. Pole Transform of Spherical Harmonic Function to Quickly Calculate the Disturbing Gravity on Earth-Orbiting Satellites[J]. Geomatics and Information Science of Wuhan University, 2013, 38(9): 1039–1043 ( 王建强, 李建成, 王正涛, 等. 球谐函数变换快速计算扰动引力[J]. 武汉大学学报·信息科学版, 2013, 38(9): 1039–1043. )
[6] Wu Xiaoping. Application of Data in the Determination of the External Disturbing Gravity Field Outside the Earth[J]. Journal of PLA Institute of Surveying and Mapping, 1992(4): 1–10 ( 吴晓平. 在推求地球外部扰动重力场中数据的采用[J]. 解放军测绘学院学报, 1992(4): 1–10. )
[7] Bjerhammar A. Discrete Physical Geodesy[R]. The Ohio State University, Ohio, USA, 1987
[8] Wu Xiaoping. Point-Mass Model of Local Gravity Field[J]. Acta Geodaetica et Cartographica Sinica, 1984, 13(4): 249–258 ( 吴晓平. 局部重力场的点质量模型[J]. 测绘学报, 1984, 13(4): 249–258. )
[9] Huang Motao. Calculation and Study of Gravity Disturbances in Ballistic Computation for Strategic Submarine Missile[D]. Zhengzhou: Institute of Surveying and Mapping, 1991 ( 黄谟涛. 潜地战略导弹弹道扰动引力计算与研究[D]. 郑州: 解放军测绘学院, 1991 )
[10] Huang Jinshui, Zhu Zhuowen. Frequency Response Point Masses Model of Exterior Disturbing Gravity Field[J]. Acta Geophysica Sinica, 1995, 38(2): 182–188 ( 黄金水, 朱灼文. 外部扰动重力场的频谱响应质点模型[J]. 地球物理学报, 1995, 38(2): 182–188. )
[11] Zhao Dongming, Wu Xiaoping. A Substitute Method for Fast Determination of Disturbing Gravity[J]. Journal of PLA Institute of Surveying and Mapping, 2001, 18(B09): 11–13 ( 赵东明, 吴晓平. 扰动引力快速确定的替代算法[J]. 解放军测绘学院学报, 2001, 18(B09): 11–13. )
[12] Zhang Hao. A Study on Fast Computation of Trajectory Disturbing Gravity[D]. Zhengzhou: Information Engineering University, 2007 ( 张嗥. 快速逼近弹道扰动引力的算法研究[D]. 郑州: 信息工程大学, 2007 ) http://cdmd.cnki.com.cn/article/cdmd-90008-2008044481.htm
[13] Zheng Wei, Qian Shan, Tang Guojian. Fast Computation of Disturbing Gravity in Ballistic Missile Guidance[J]. Flying Mechanic, 2007, 25(3): 42–44 ( 郑伟, 钱山, 汤国建. 弹道导弹制导计算中扰动引力的快速赋值[J]. 飞行力学, 2007, 25(3): 42–44. )
[14] Jiang Dong, Wang Qingbin, Zhao Dongming. Analysis on Efficiency of the Methods for Fast Computing Disturbing Gravity[J]. Journal of Geomatics Science and Technology, 2011, 28(6): 411–415 ( 江东, 王庆宾, 赵东明. 空中扰动引力快速赋值算法的效能分析[J]. 测绘科学技术学报, 2011, 28(6): 411–415. )
[15] Chen Guoqiang. Spacecraft Dynamics in Gravity Anomaly Field[M]. Changsha: National University of Defense Technology Press, 1982 ( 陈国强. 异常重力场中飞行器动力学[M]. 长沙: 国防科技大学出版社, 1982 )
[16] Jia Peiran, Chen Kejun, He Li. Long Range Rocket Ballistics[M]. Changsha: National University of Defense Technology Press, 1993 ( 贾沛然, 陈克俊, 何力. 远程火箭弹道学[M]. 长沙: 国防科技大学出版社, 1993 )
[17] Li Jiancheng, Chen Junyong, Ning Jinsheng, et al. Theory of the Earth's Gravity Field Approximation and Determination of China Quasi-Geoid 2000[M]. Wuhan: Wuhan University Press, 2003 ( 李建成, 陈俊勇, 宁津生, 等. 地球重力场逼近理论与中国2000似大地水准面的确定[M]. 武汉: 武汉大学出版社, 2003 )
[18] Huang Motao, Ouyang Yongzhong, Liu Min, et al. Regulation of Point Mass Model for Multi-Source Gravity Data Fusion Processing[J]. Geomatics and Information Science of Wuhan University, 2015, 40(2): 170–175 ( 黄谟涛, 欧阳永忠, 刘敏, 等. 融合海域多源重力数据的正则化点质量方法[J]. 武汉大学学报·信息科学版, 2015, 40(2): 170–175. )
[19] Huang Motao, Guan Zheng. Test and Construction of Disturbing Point Mass Model[J]. Journal of Wuhan Technical University of Surveying and Mapping, 1994, 19(4): 304–309 ( 黄谟涛, 管铮. 扰动点质量模型构制与检验[J]. 武汉测绘科技大学学报, 1994, 19(4): 304–309. )
[20] Huang Motao, Guan Zheng, Ouyang Yongzhong. Accuracy Analysis and Calculation of 1°×1°point Masses in the Area of China[J]. Journal of Wuhan Technical University of Surveying and Mapping, 1995, 20(3): 257–262 ( 黄谟涛, 管铮, 欧阳永忠. 中国地区1°×1°点质量解算与精度分析[J]. 武汉测绘科技大学学报, 1995, 20(3): 257–262. )
[21] Huang Motao. Configuration Optimization of the Disturbing Point Masses Model and Its Sequential Solution[J]. Acta Geodaetica et Cartographica Sinica, 1994, 23(2): 81–89 ( 黄谟涛. 扰动质点赋值模式结构优化及序贯解法[J]. 测绘学报, 1994, 23(2): 81–89. )
[22] Feng Kang. Method of Numerical Computation[M]. Beijing: National Defence Industry Press, 1978 ( 冯康. 数值计算方法[M]. 北京: 国防工业出版社, 1978 )
[23] Zhang Chuanding, Lu Zhonglian, Wu Xiaoping. Truncation Error Formulae for the Disturbing Gravity Vector[J]. Journal of Geodesy, 1998, 72(3): 119–123 DOI:10.1007/s001900050153