2. 大地测量与地球动力学国家重点实验室, 武汉 430077;
3. 中国科学院大学地球科学学院, 北京 100049;
4. 北京环球信息应用开发中心, 北京 100094;
5. 61287部队, 成都 610036
2. Stake Key Laboratory of Geodesy and Earth's Dynamics, Wuhan 430077, China;
3. College of Earth Science, University of Chinese Academy of Sciences, Beijing 100049, China;
4. Global Information Application and Development Center of Beijing, Beijing 100094, China;
5. 61287 Troops, Chengdu 610036, China
地壳均衡理论自诞生以来就一直应用于重力场数据处理中,均衡异常的优点在于变化平缓,更有利于内插和推估.早在1937年,以普拉特均衡重力异常等于零为条件,Dubovskii解出了一个完全至6阶的重力场模型;Heiskanen和Uotila分别于1938和1962年应用均衡重力异常以及空间重力异常对重力场模型参数进行 了推估;SjÖberg推导出了Airy-Heiskanen(A-H)地形/均衡重力位和相应的重力异常(1998).
Rapp(1977, 1982)采用球面模型计算均衡补偿时,当补偿深度为50 km,Airy模型与观测值的符合程度最好,至少达到150阶.OSU89B模型在处理没有重力观测值区域的重力异常估值时采用了由高程数据基于地形/均衡假说的球面模型计算的异常值(Pavlis et al.,1990).
地形/均衡重力场模型是研究地壳均衡理论的应用方向之一,其不仅能够反映地球内部构造与地形质量分布的对应关系,而且为估计地球的扰动质量提供了一个可行的思路.Tsoulis(2004)应用均衡模型基于凝聚到球面的质量分布假设得到了完全至180阶的重力场模型,该模型反映了相应分辨率的大地水准面高,并对地球内部结构有一定的解释.Makhloof(2007)给出了三种地形/均衡质量模型(Pratt-Hayford(P-H)、A-H、混合模型)引力位的1、2阶导数的影响公式,并成功将基于均衡理论的球近似积分公式应用于引力张量分量向下延拓的数值模拟计算中.Bagherbandi(2012)的研究表明地形/均衡补偿40~60阶次时的地球重力场模型较仅仅采用地形数据时构建的模型有所改善.在EGM 系列重力场模型的构建中,重力空白区也是利用全 球均衡模型推估的空间重力异常填补的.
均衡理论的应用是基于地形数据和地壳模型展开的,原有的研究受限于数据分辨率、算法等原因,主要集中于区域和低阶重力场模型构建,均衡理论应用于超高阶地球重力场模型构建的研究较少.
与传统重力场模型相比,地形/均衡重力场模型未顾及实测重力数据,因此,应用前需要进行实测修正.以重力空白区填补为例,由于均衡假设推估的异常与实测值有一定的差异,直接填补还包含有一个区域性系统误差.顾及该系统差有两种办法,一是利用实测一定分辨率格网数据求残差来确定系统差;二是利用现有重力位模型,截断至一定分辨率进行拟稳分析,滤除大尺度系统差,提高填补精度.
本文初步尝试应用三种均衡补偿机制构建超高阶次地球重力场模型,对比分析补偿深度的调整对不同均衡补偿机制的影响;并以EGM2008为参考模型,截断至不同阶次为拟稳条件,分别对地形/均衡重力位模型进行系统差纠正,滤除大尺度系统差,对比分析拟稳模型的精度.
2 地形/均衡重力场位模型 2.1 A-H均衡模型A-H均衡模型原理示意图(Hofmann-Wellenhof et al.,2005, 莫里茨 1992)见图 1. 图中h为地形高,h'为海底深度,T0为补偿深度,t为山根,t'为反山根;ρc为地壳的平均密度,ρm为地幔的平均密度,ρw为海水的平均密度.为了简化模型,将地表仅分为陆地和海洋两部分.在地形模型SRTM(Shuttle Radar Topography Mission)中,认为正值为陆地,负值为海洋,为了统一计算公式,定义等效地形高h″(Rummel et al.,1988):
式中,Re为地球平均半径.
考虑到地球质体对外部点的引力位,令Q(rQ,θQ,λQ)为地形或补偿质量的一个微分单元(从极限的角度看可以认为是一个点),P(rP,θP,λP)即为球外或球面上的计算点,见图 2.
式中(nm,nm)为重力位系数,(n,m)为阶和次,nm是正常化的伴随Legendre多项式, 为地球平均密度,dsQ=sinθQdθQdλQ为积分面元,(ATop(Q),ACom(Q))为地形质量和补偿质量对外部点的影响,将二者积分后代入上下限,并依泰勒级数展开至三阶项(hQ/Re)3,则有A-H位模型(Anm,Anm)计算公式(Rummel et al.,1988; Tsoulis,2001):
2.2 面凝聚位模型
在A-H均衡模型的基础上,做如下假设:可见的地形质量凝聚在海平面上,即地形质量中流动点Q球坐标的径向分量rQ为地球平均半径Re.同时,在地形质量以下,补偿深度为T0处的补偿质量与其径向方向上的地形质量在数值上相等,符号相反,故面凝聚位模型(Lnm,Lnm)计算式为(Tsoulis et al.,2004)
式中, χL(n) = 3R<sup>n+2e-3(Re-T0)n+2 4(2n+1)π eR<sup>n+3e . 2.3 Airy-凝聚组合均衡模型
A-H均衡模型的补偿计算中,需要依据地形数据计算出相应的山根或者反山根,山根即是描述陆地地区补偿面向下延伸到地幔部分的质量,而反山根则是用于描述海洋地区的盈余部分质量.由于山根t的值与h″成正比,故在海洋地区,当h″足够大时,反山根将向上隆起接近海底甚至高出海底,这显然与实际情况不符.因此,本节组合A-H均衡模型和面凝聚均衡模型,以期得到更加符合实际情况的位模型,即陆地部分采用A-H均衡模型,而海洋部分则应用面凝聚均衡模型,综合二者得到组合均衡位模型(Hirt et al.,2012).
故均衡重力位模型的统一表达式为
3 轮胎调和分析及拟稳分析 3.1 轮胎调和分析
分析式(6)不难发现,被积函数f(θ,λ)是流动点Q的函数,Q以离散的等间隔球面格网形式分布在球面上,通常该格网是以东经0°和北纬90°为起始点,自东向西,由北到南的顺序给出,且取格网中点值为格网均值,因此,解算时通常将积分形式转换成离散求和形式,即基于球面边值的调和分析过程.
传统的调和分析中,由于平滑因子无法准确求得,其解算结果的逼近效果有限.另外,由于数据分辨率与模型阶次满足尼奎斯特法则,高分辨率数据对应高阶次重力场模型,因此,高阶次模型的解算存在递推稳定性和效率问题,应用轮胎调和分析可以有效地解决这一问题.
自洽的轮胎调和分析方法利用球面到“轮胎面”的映射关系,能够克服平滑因子不精确引起的误差,并且采用跨阶次递推方法,实现稳定、快速、高速地恢复重力场模型.其将高精度的调和分析分解为样条分析、傅里叶分析和调和分析三步,能够实现以99.9%的精度还原地球重力场,实际计算中,式(6) 即采用轮胎调和分析方法进行解算(张传定等,2004).
3.2 拟稳分析重力场模型逼近和恢复是物理大地测量学的重要部分,而现有的技术手段只能通过有限阶次的级数形式表示,随着观测手段、技术理论的发展,恢复模型的阶次在不断提高,尽管如此,其也仅仅是实际地球重力场模型的一个近似表示.以精度较高的重力场模型为基准,对计算所得均衡位模型进行系统差的滤除工作,以期提高均衡位模型的精度,参考拟稳平差原理,本文称其为拟稳分析.
EGM2008模型是国际上发布的基于地面重力异常数据、卫星跟踪卫星和卫星测高数据等构制的 高精度、高分辨率地球重力场模型(章传银等,2009). 以截断至一定阶次的EGM2008模型重力异常为拟稳基准进行拟稳分析,能够充分顾及现有的高精度全球重力异常模型和基于均衡理论推估的重力异常模型,从而消除大尺度系统差,提高推估重力异常的精度.
以下为拟稳分析具体过程.重力异常计算公式为(陆仲连,1996)
式中,G为万有引力常数.计算位模型对应的5'、30'、1°、2°分辨率的重力异常,以EGM2008重力异常为基准,计算得30'、1°、2°分辨率的均衡位模型重力异常系统残差,将获得的三组残差都采用最小曲率曲面拟合法(Smith et al.,1990)加密为5',分别用于校正均衡位模型5'分辨率的重力异常,将校正后的重力异常作为已知边界条件反演完全至2160阶次的重力场位模型,即拟稳模型.最后,对所得的三个拟稳模型进行对比分析.
位模型阶方差计算公式为(Tsoulis,2001)
大地水准面和重力异常相对阶误差计算式为(刘晓刚,2011;吴星,2009)
式中,( <sup>EGM<sub>nm, <sup>EGM<sub>nm)为参考模型位系数,δN为大地水准面高差,δΔg为重力异常差.
模型相对累计大地水准面高误差和累计重力异常误差计算公式为(VÖlgyesi et al 1992)
4 数值实验与精度分析
SRTM提供了陆地和海深数据,采用平均为5'分辨率的SRTM模型,依据式(6)进行轮胎调和分析(张传定等,2004),得完全至2160阶次的重力场位 模型系数.其中地壳的平均密度为ρc=2670 kg·m-3, 地幔的平均密度为ρm=3270 kg·m-3,海水的平均密度为ρw=1025 kg·m-3,地球平均半径为Re=6371 km,地球的平均密度为=5500 kg·m-3.
4.1 均衡位模型对比分析从式(6)可以看出,补偿深度的选择对位系数模型的影响较大.因此,在考虑常用的补偿深度30 km的基础上,分析了补偿深度的变化对两种补偿模型的影响,比较了相同的补偿深度在两种补偿模型中的贡献.
为了全面对比分析均衡模型的贡献特征,图 1—9分别给出了相同补偿深度不同均衡模型、同一均衡模型不同补偿深度的位模型与EGM2008(Pavlis et al.,2008)参考模型的阶方差、相对大地水准面阶误差、相对重力异常阶误差的对比分析情况.图中,Layer和Airy分别指代面凝聚位模型和A-H位模型.
由图 3—5可得到如下结论:
(1) 两种均衡模型对补偿深度都较为敏感,补偿深度相同时,在540阶次的范围内,A-H位模型的阶方差曲线在面凝聚位模型的下方,其更接近参考模型的阶方差曲线,这是由于A-H位模型的补偿机制较面凝聚位模型更接近实际地球质量分布,更能反映较为真实的扰动质量分布所致;
(2) 从重力场扰动场元的角度分析,无论补偿深度为80 km还是40 km,A-H位模型的有效阶次在高阶部分都要比面凝聚位模型高些,补偿深度为40 km时两者在低阶部分较为接近.
4.2 补偿深度的对比分析补偿深度是分析均衡模型对重力场贡献的重要参数,其对两种均衡补偿机制的影响见图 6—9.
对比分析图 6—9,可知:
(1) 对任一补偿模型,补偿深度越浅,模型的阶方差曲线对应的能量谱在540阶次以内越小,补偿效果越明显,而补偿深度的变化对两个模型超高阶部分的补偿效果没有影响;
(2) 面凝聚位模型中,补偿深度为40 km时的相对大地水准面阶误差精度在380阶次以内要高于补偿深度为80 km时的模型,随着阶次的增大,模型的精度逐渐降低,到一定阶次时,补偿深度贡献的差异几乎为零;
(3) A-H位模型中,补偿深度为50 km时的相对大地水准面阶误差的精度在440阶次以内要高于补偿深度为80 km时的模型.
综合补偿深度和模型补偿机制,本研究中,对参数补偿深度的设置是以10 km为间隔进行的,通过对两种补偿机制下不同补偿深度构建的地形/均衡位模型与参考EGM2008对比,选定相对最优的补 偿深度值作为后续研究的基础,即补偿深度为50 km 的A-H位模型和补偿深度为40km的面凝聚位模型.图 10、表1和表2给出了相应重力异常差的空间分布以及统计结果,模型计算所得重力异常均为地球表面的值,且不同分辨率的重力异常是通过模型截断至相应阶次分别计算所得,从而得到模型重力异常差,空间分布图以30'分辨率的数据为例图示.
对比分析图 10、表1和表2,可知:
从模型重力异常估值与EGM2008模型差值分布的角度看,A-H位模型的推估效果优于面凝聚位模型,A-H位模型不同分辨率空间异常差值平均值的绝对值要小于面凝聚位模型.两种推估模型的整体趋势较为接近,且不同分辨率空间异常差值平均 值的绝对值在1 mGal附近,均方差大约为20 mGal.
前文选取的EGM2008参考模型在构制过程中,采用地形模型对缺乏重力数据的区域进行了推估,为了避免地形数据的相关性对本文研究的影响,本小节选取GOCO02s纯卫星重力模型作为参考,对A-H位模型、面凝聚位模型和Airy-凝聚组合模型进行精度对比分析.
GOCO02s模型是完全至250阶次的全球重力场位模型,其融合了GOCE、GRACE、CHAMP以及SLR等卫星数据.引入以下统计量对上述模型进行分析,即阶次n1和n2间的模型相对阶误差、大地水准面高和重力异常的相对均方根:
分析不同补偿深度对均衡位模型不同频段的影响,设定组合模型中陆地部分补偿深度为50 km,海洋部分的补偿深度分别为20 km、30 km和40 km(即Combined 20、Combined 30和Combined 40), 然后分别计算上述六个模型在2~36、37~50、51~120 和121~250阶等四个区间的各统计量分布值,以GOCO02s 计算值为参考,分析各均衡位模型的精度,具体计算结果如表3—5所示.
分析表3—5,可以得到如下结论:
(1) 补偿深度为40 km的面凝聚位模型在前120阶的精度要高于补偿深度为50 km的A-H位模型,而在121~250阶的精度较低;
(2) 补偿深度为40 km的组合模型在2~36阶和51~120阶的精度高于其他各均衡位模型,而补 偿深度为30 km的组合模型在37~50阶的精度最好,从重力异常均方根的角度看,补偿深度为20 km的组合模型的精度最好;
(3) 组合模型的精度与补偿深度密切相关,调整补偿深度,可以得到相应频段的最佳均衡位模型,其精度要高于其他组合模型和单一均衡位模型.
4.4 EGM2008模型拟稳后误差分析对补偿深度为50 km的A-H位模型进行EGM2008 拟稳分析,以30'、1°、2°三个分辨率的EGM2008模型重力异常作为拟稳基准,分别将相 应分辨率的A-H位模型的重力异常拟稳到EGM2008 上,并分析这几种分辨率数据的拟稳结果.
为了更清晰地说明问题,下面仅图示出30'和2°分辨率得到的拟稳模型阶方差,并将二者在不同的频段分别与EGM2008和A-H位模型进行对比.
图 11中两种拟稳模型阶方差在低阶部分与EGM2008一致性较好;且30'分辨率对应的阶方差与参考阶方差曲线直到360阶都有较好的一致性,较校正前的一致性范围有所提高;而2°分辨率对应的阶方差与参考阶方差曲线在90阶内的一致性较好,但是在90阶之后,该模型对应的阶方差曲线与参考曲线的偏差较大.
拟稳模型的相对累计大地水准面高误差和累计重力异常误差可以从另外的角度反映不同基准对均衡位模型的影响,具体计算结果见图 13—14以及表6和表7.
从图 13—14可以看出,拟稳基准的分辨率越高,拟稳的效果越好.不同分辨率能够改善的最大阶次也有所不同,从相对累计大地水准面的角度看,30'、1°和2°分辨率的基准能够改善的模型阶次分别为140、88和11;而从相对累计重力异常的角度看,相应的阶次分别为220、111和22.大地水准面主要覆盖低模型的低阶部分,而重力异常主要覆盖至模型的中高阶部分,可见,拟稳基准能够提高原模型在中高阶部分的精度,累计大地水准面高误差分别趋向于各自的最大值,而累计重力异常误差在超高阶部分的误差仍存在于三种拟稳模型之中.
从表6和表7可以看出,随着模型阶次的增加,累计误差都在增大;在各分辨率拟稳基准能够改善的模型阶次范围内,分辨率越高,两种统计累计误差越小,30'分辨率对应的模型在140阶时累计大地水准面误差达到0.068 m,而在220阶时的累计重力异常误差为1.101 mGal.
依据尼奎斯特法则,从图 14可以看出,三种拟稳基准能够改善分辨率对应阶次的A-H位模型,且分辨率越高,在低阶部分改善的效果越明显.
5 结论与展望本文研究了以5'分辨率的SRTM模型为基础,依据不同的均衡补偿机制,得到的完全至2160阶次的地形/均衡位系数模型,并以不同分辨率的EGM2008模型重力异常为基准,进行了拟稳分析研究.
数值实验表明,A-H位模型和面凝聚位模型对应的最佳补偿深度分别为50 km和40 km,此时二者在不同频段的优劣性不一,前250阶二者精度低于Airy-面凝聚组合模型.此外,补偿深度的选取对均衡补偿在超高阶部分几乎没有影响.
地形/均衡位模型的有效阶次与扰动场元的能量分布有关,不同的扰动场元对均衡补偿机制的响应亦有所区别,响应的频段有一定差异.相对累计重力异常差能够覆盖到重力场中高频信号,其显示拟稳模型在高阶部分具有一定的改善效果.30'分辨率模型重力异常拟稳后的位模型在140阶时相对累计大地水准面误差为0.068 m,在220阶时相对累计重力异常误差为1.101 mGal.
从各模型与参考模型的对比分析中不难发现,均衡位模型在低阶和超高阶部分的逼近效果较差,调节补偿深度,位模型的响应阶次约在37~440频段内,因此,在实际应用过程中,应充分考虑均衡理论的这一特性,综合均衡补偿机制和其他各类实测数据,才能更好地为地球重力场模型的恢复服务.
本文中所采用的数据仅仅为地形数据,没有顾及地壳密度和厚度的实际变化情况,故研究存在一定的假设基础,为进一步分析拟稳模型在重力数据空白区的填补应用提供基础,若能在均衡位模型的构建中考虑更加接近地球实际密度分布的情况,将不仅能够检验不同均衡模型的适用性,而且能够为均衡理论应用于实际重力测量数据处理提供理论基础.
[1] | Bagherbandi M, Sjöberg L E. 2012. A synthetic Earth gravity model based on a topographic-isostatic model. Stud. Geophys. Geod., 56(4): 935-955. |
[2] | Hirt C, Kuhn M, Featherstone W E, et al. 2012. Topographic/isostatic evaluation of new-generation GOCE gravity field models. J. Geophys. Res., 117, B05407, doi:10.1029/2011JB008878. |
[3] | Hofmann-Wellenhof B, Moritz H. 2005. Physical Geodesy. Graz: Springer. |
[4] | Liu X G. 2011. Theory and methods of the earth's gravity field recovery from GOCE data [Ph. D. thesis](in Chinese). Zhengzhou: PLA Information Engineering University. |
[5] | Lu Z L. 1996. Theory and Methods of the Earth's Gravity Field (in Chinese). Beijing: PLA Publishing House. |
[6] | Makhloof A A, Ilk K. 2007. Effects of topographic-isostatic masses on gravitational functionals at the Earth's surface and at airborne and satellite altitudes. J. Geod., 82(2): 93-111. |
[7] | Moritz H. 1992. The Figure of the Earth: Theoretical Geodesy and the Earth's Interior (in Chinese). Chen J Y, Zuo C H, Trans. Beijing: Surveying and Mapping Press. |
[8] | Pavlis N K, Holmes S A, Kenyon S C, et al. 2008. An Earth gravitational model to degree 2160: EGM 2008. Vienna, Austria. |
[9] | Pavlis N K, Rapp R H. 1990. The development of an isostatic gravitational model to degree 360 and its use in global gravity modelling. Geophys. J. Int., 100(3): 369-378. |
[10] | Rapp R H. 1977. Determination of the potential coefficients to degree 52 from 5° mean gravity anomalies. Bull. Geod., 51(4): 301-323. |
[11] | Rapp R H. 1982. Degree variances of the Earth's potential, topography and its isostatic compensation. Bull. Geod., 56(2): 84-94. |
[12] | Rummel R, Rapp R, Sunkel H, et al. 1988. Comparisons of global topographic isostatic models to the Earth's observed gravity field. Columbus: Department of Geodetic Science and Surveying, The Ohio State University. |
[13] | Sjöberg L E. 1998. The exterior Airy/Heiskanen topographic-isostatic gravity potential, anomaly and the effect of analytical continuation in Stokes' formula. J. Geod., 72(11): 654-662. |
[14] | Smith W H F, Wessel P. 1990. Gridding with continuous curvature splines in tension. Geophysics, 55(3): 293-305. |
[15] | Tsoulis D. 2001. A comparison between the Airy-Heiskanen and the Pratt-Hayford isostatic models for the computation of potential harmonic coefficients. J. Geod., 74(9): 637-643. |
[16] | Tsoulis D, Stary B. 2004. An isostatically compensated gravity model using spherical layer distributions. J. Geod., 78(7-8): 418-424. |
[17] | Völgyesi L, Tóth G Y. 1992. Optimal topographic-isostatic crust models for global geopotential interpretation. Periodica Polytechnica Civ. Eng., 36(2): 207-241. |
[18] | Wu X. 2009. Theory and methods of satellite gradiometry data processing [Ph. D. thesis] (in Chinese). Zhengzhou: PLA Information Engineering University. |
[19] | Zhang C D, Xu H Z, Wu X. 2004. Torus problem of the harmonic analysis in the Earth's gravity field. // Proceedings of Geodesy and Geodynamics Development (in Chinese). Wuhan: Science and Technology Press. |
[20] | Zhang C Y, Guo C X, Chen J Y, et al. 2009. EGM 2008 and its application analysis in Chinese mainland. Acta Geodaetica et Cartographica Sinica (in Chinese), 38(4): 283-289. |
[21] | 刘晓刚. 2011. GOCE卫星测量恢复地球重力场模型的理论与方法[博士论文]. 郑州: 解放军信息工程大学. |
[22] | 陆仲连. 1996. 地球重力场理论与方法. 北京: 解放军出版社. |
[23] | 莫里茨 H. 1992. 地球形状——理论大地测量学和地球内部物理学. 陈俊勇,左传惠译. 北京: 测绘出版社. |
[24] | 吴星. 2009. 卫星重力梯度数据处理理论与方法[博士论文]. 郑州: 解放军信息工程大学. |
[25] | 张传定,许厚泽,吴星. 2004. 地球重力场调和分析中的“轮胎”问题. 大地测量与地球动力学进展. 武汉: 科学技术出版社. |
[26] | 章传银,郭春喜,陈俊勇等. 2009. EGM 2008地球重力场模型在中国大陆适用性分析. 测绘学报, 38(4): 283-289. |