2. 武汉大学测绘遥感信息工程国家重点实验室，湖北 武汉430070
2. State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing, Wuhan University, Wuhan 430070, China
1 引 言
地球重力场是地球的基本物理场，是探测地球内部物质密度分布、研究内部结构与地质构造的重要手段。对于局部地区的重力资料而言，可以将重力观测面近似为平面，所有数据处理问题都可以在直角坐标系中进行理论推导和计算。然而，对于区域性的，乃至全球尺度问题，重力观测面已不能近似为平面，而是一个球面问题，需要利用球坐标进行表达与处理。此时，建立在直角坐标系中的数据处理、分析及反演解释方法已不再适用。随着卫星重力观测技术的发展，地球、月球与火星等星体的重力场模型的精度和可靠性都有了大幅度地提高[1, 2, 3, 4, 5, 6, 7]。欲利用重力场研究星体的内部结构与构造，发展球坐标系中的重力异常正演方法显得尤其重要。因为它是在球坐标系中进行重力地形校正、均衡校正、“剥层法”分离异常、重力场合成计算、界面起伏反演、局部和全球密度成像的基础。此外，2002年3月发射的GRACE卫星提供了地球重力场的时变信息，它反映了地球系统的物质质量迁移[8]。在利用GRACE时变重力场研究地球深部动力学过程时，必须扣除大气圈、水圈及地下水等物质的迁移与交换引起的重力变化，因此也亟须发展球坐标系内高精度的重力正演方法。
2 基于球坐标的重力异常正演方法 2.1 方法回顾球坐标系中任意地质体对其外部空间产生的引力位均可由牛顿万有引力积分公式表达
式中，G为万有引力常数；Ω指该地质体；dm与dv为积分质量元与体积元；l为观测点与积分单元的距离；u(r,θ,λ)为观测点处得引力位；ρ(r′,θ′,λ′)为地质体的密度或异常密度分布。
将式(1)对径向r求导，即得到重力或重力异常
在卫星重力的全球地形校正中，可以分别在球谐域与空间域中进行计算。与文献[9]在频率域快速计算位场异常类似，文献[10]将其推广到球谐域进行月球全球的地形重力效应的快速计算中，该方法的计算精度与稳定性受全球地形模型阶次、地形级数展开的次数以及观测高度等的影响，而且能否在局部区域的重力地形校正中成功应用尚值得探讨。在空间域，一般将全球地形剖分成单元体，计算每个单元体产生的引力效应，然后累加求和即得到全球地形的重力效应。其中，采用最普遍的剖分方法是Tesseroid单元体法。其概念最初由文献[11]引入，与传统直角坐标系中的长方体单元类似，Tesseroid表示的是球坐标系下的一个球面棱柱体单元(图 1)，其6个界面由一对固定半径r1、r2且“平行”于参考球面的球面，一对经线平面λ和一对同轴圆锥角φ组成。可知，若单元体内外界面半径r1=R+h1，r2=R+h2，R表示参考球面的半径，当h1=0时，h2可作为经纬度网格下数字地形模型的高程参数。
|
| 图 1 Tesseroid单元体几何示意图 Fig. 1 Geometry of the Tesseroid |
然而，该单元体对应的积分公式(1)不存在严格的解析解，其解决办法主要有4种：① 采用数值积分方法[12]，包括三维高斯-勒让德体积分(3-D Gauss-Legendre cubature)[13]，或者对于公式(1)先对径向r积分转化为二维球面积分问题，再采用数值积分的方法[12]；② 采用直立棱柱体、点元、线元或者面元来近似计算Tesseroid单元体对外部空间任一点的重力效应[12]，该方法计算精度低，采用等质量的直立棱柱体近似计算时，还需要进行直角坐标系的旋转；③ 以单元体几何中心作为泰勒级数展开点近似计算式(1)，具体计算公式参见文献[12, 14]；④ 基于球冠的重力异常正演方法，其基本思路是一个均匀等厚的球冠层(图 2)对北极轴上任一点的重力效应存在解析解，进而利用位场叠加原理与模型剖分的对称性，可以计算任一个Tesseroid单元体对北极轴上任一点的重力效应。但是在计算地形对其他点的重力效应时，需要进行球坐标系旋转与地形模型重构，文献[15]将该方法运用到月球重力地形校正中，并且与文献[16]采用Tesseroid单元体基于泰勒级数展开方法的计算结果进行对比，验证基于球冠的重力异常正演方法的可行性。
|
| 图 2 球冠层几何示意图 Fig. 2 Geometry of the spherical cap |
本文相对于球冠面积分重力正演方法的改进有3点：① 均匀等厚球冠层物质对极轴上某点重力效应计算公式的改进，采用完全的解析解，而不再是球面积分的近似解；② 与文献[15]中分面向与背向计算点的球冠层重力效应计算不同，统一了面向与背向两种情况；③ 球坐标系旋转后的模型重构方法与策略进行改进，提高了计算效率和精度。
2.2.1 均匀等厚球冠层物质对极轴上某点的引力效应如图 2所示为一个等厚而且物质密度分布均匀的球冠层，其对极轴上某点P(r,ψ=0)(如果计算点在球冠层内，则需要将球冠层以该点所在球面为界面分成两层分别计算)的重力效应为[14]
式中，Δρ为均匀等厚球冠层的剩余密度；ψC为球冠层边缘点的极角；l′C为计算点P(r,ψ=0)与球冠层外缘之间的距离
当计算点位于球冠层上、下表面时，式(3)可以分别简化为
经数值计算与对比，式(5)和式(6)与文献[17]中的式(1)和式(3)以及与文献[18]中的式(7)和式(10)计算结果相一致。式(3)可以用于地面与海洋、航空及卫星重力地形校正中，建议在地面与海洋船测重力地形校正中采用式(5)与式(6)计算以提高运算速度。另外需要说明的是，式(3)与式(4)未写成R与h的形式，是因为这样可以使得该公式既适合于球体又适合于椭球体为参考面的Tesseroid单元体对极轴上某点重力效应的计算。
若将球面按经纬度剖分，由北极点至南极点各纬圈的余纬分别为ψC0、ψC1、ψC2、…、ψCN ，其中ψC0=0,ψCN≤π，
，
，N与M分别为纬度与经度方向的网格剖分块数，dlat与dlon分别为纬度与经度方向的网格剖分大小，则由位场叠加原理与模型剖分的对称性，按经纬度剖分的单元体物质对极轴上某点的重力效应为
则全球地形起伏物质对极轴上某点重力效应为
2.2.2 面向与背向计算点的球冠层物质重力效应计算公式的统一
事实上，式(3)适合全球计算，即0<ψC≤π，当ψC=π时，则[19]
令h=r2－r1，若r2≈r,r1≈r
若把均匀等厚球冠层视为一个球面(当球冠层厚度远小于球半径时)，则文献[15]中的式(3)与式(5)可以统一写成
当ψC为π且计算点位于球面时，则得到与式(10)相同的结果。同理，文献[15]中的式(8)与式(9)可以统一成为
此处，l0与ψ0分别对应着r－R与0，l为校正点到参考球面上计算面元的距离。若将l替换为l－(r2－r1)/2，计算精度将得到一定程度的提高。 2.2.3 球坐标系旋转后地形模型重构方法的改进
球坐标系地形校正P点与球心的连线是地轴方向，而在计算其他任意点的地形重力效应时，相当于将地轴沿某一经度线作了旋转。由于地形高程数据是与原坐标系中经纬度一一对应的，那么在旋转后的新坐标系中地形数据不再是规则的按等间隔经纬度剖分的，那么需要重构地形模型。文献[15]将原坐标系的地形数据全部经球坐标转换至新的球坐标系中，在重构新坐标系中经纬度网格点的地形时，通过搜索该点周围的点进而求取平均值。
为提高计算效率，首先将新坐标系中待求点的新球坐标反变换到原坐标系(由于原坐标系中数据分布有规律，可以直接很快地找到该点周围的4个点)，通过周围4个点高程的平均或双线性插值得到新坐标系内待求点的值。坐标反变换时只需将文献[15]中式(15)新坐标系的Z轴与球面的交点(r,θ,λ)变换为(r,θ,λ+π)即可。4点的双线性插值公式[20]为
式中，θi－1≤θ≤θi为余纬度，xi－1≤x≤xi,x=cosθ,λi－1≤λ≤λi为经度，系数为 式中，Δ=(xi－xi－1)(λi－λi－1)。 3 计算稳定性及精度分析
采用的球冠体积分重力正演公式具有严格的解析解，在重力异常正演中，其计算精度取决于原始地形等地质界面模型精度与球坐标系旋转后地形模型重构的精度。若先用球谐函数或者其他函数拟合原始的全球或者局部地形数据，在重构模型时，新坐标系中的值经坐标反变换到原坐标系经过模型直接解算获得。但是如此计算效率太低，往往地球、月球等星体的地形均采用球谐系数表达，其分辨率有限。由于在小于分辨率的尺度上，模型值之间的变化可以看成线性的，因此采用球坐标逆变换与双线性插值方法既可以满足精度又可以提高计算效率。而Tesseroid单元体基于其几何中心点作泰勒级数展开的方法与基于球冠面积分方法的计算精度可能受级数展开次数及Tesseroid单元体表面大小、厚度、深度、计算点高度与表面地理位置的影响。由于在计算极轴上点的重力效应存在解析解，所以探讨各种因素对计算极轴上某点的重力效应的影响。
如图 3与表 1，知Tesseroid单元体基于其几何中心作泰勒级数展开方法的级数展开次数越高，计算精度越高。图 3各分图中1为理论值，2、3与4分别为球冠面积分方法、Tesseroid单元体基于泰勒级数展开零阶近似方法与二阶近似方法计算值与理论值的差异。Tesseroid单元体基于泰勒级数展开零阶近似的方法受单元体大小影响较大；Tesseroid单元体基于泰勒级数展开二阶近似的方法受单元体厚度影响最小，且厚度越厚，误差皆增大；球冠面积分方法受单元度影响最小；Tesseroid单元体基于泰勒级数展开二阶近似的方法受计算点高度影响最小；Tesseroid单元体基于泰勒级数展开方法受单元体球面地理位置影响严重，且越靠近计算点，计算值越偏离理论值。
|
| 图 3 重力异常正演方法计算精度分析 Fig. 3 Accuracy analysis of different methods for gravity anomaly calculation |
| mGal | ||||||||||
| 单元体表面积 | 单元体厚度 | 单元体深度 | 计算点高度 | 表面地理位置 | ||||||
| 平均值 | 均方差 | 平均值 | 均方差 | 平均值 | 均方差 | 平均值 | 均方差 | 平均值 | 均方差 | |
| 球冠面积分 | 1.297 ×10-8 |
1.161 ×10-8 |
2.105 ×10-8 |
2.386 ×10-8 |
3.888 ×10-10 |
9.687 ×10-13 |
3.875 ×10-10 |
3.620 ×10-13 |
6.20 ×10-6 |
1.641 ×10-4 |
| 零阶泰勒 级数展开 |
2.50 ×10-5 |
2.50 ×10-5 |
1.0 ×10-7 |
4.418 ×10-8 |
1.261 ×10-8 |
2.619 ×10-10 |
1.261 ×10-8 |
2.618 ×10-10 |
-4.62 ×10-4 |
1.0 ×10-2 |
| 二阶泰勒 级数展开 |
-2.77 ×10-9 |
4.629 ×10-9 |
-1.76 ×10-8 |
1.013 ×10-8 |
-1.17 ×10-8 |
6.738 ×10-11 |
3.689 ×10-14 |
1.334 ×10-13 |
-1.85 ×10-4 |
4.3 ×10-3 |
为验证各种正演方法在实际应用中的有效性与精度，将其应用于月球重力地形校正中。选取校正参考球面半径为月球平均半径1 737.013 km，地形起伏采用CLTM_s01[22]模型的1~360阶，计算半径为1750 km高度的重力影响值，校正密度选取为2920 kg/m3[23]，并且与球谐域计算结果(地形起伏展开至5次，第5次项解算的重力效应的幅值为±1 mGal)进行比较，结果如图 4。
|
| 图 4 月表地形(单位：m)及各种空间域方法计算的月球地形重力影响值与球谐域计算结果的差值(单位：mGal) Fig. 4 Topography of the Moon (unit: m) and the differences of lunar terrain gravity effection between the forward methods in spatial domain and in spherical harmonic domain (unit: mGal) |
由图 4可知，球冠面积分方法、Tesseroid单元体基于泰勒级数展开零阶近似与二阶近似方法的计算值以及改进的球冠体积分方法，与球谐域计算结果的差异范围(mGal)、均值(mGal)与均方差(mGal)分别为：-6.70~527.62、22.80与36.33，-108.02~13.36、-1.04与6.51，-32.38~32.38、0.84与3.63，-23.08~23.10、0.66与3.37。可见，球冠体积分方法计算精度最高，Tesseroid单元体基于泰勒级数展开二阶近似方法次之，球冠面积分方法精度最低。
4 结 论本文采用球冠层重力正演方法的严格解析解公式，并且统一了按面向与背向计算点的球冠层物质重力效应计算公式，新采用的公式不仅适用于地面与海洋、航空及卫星重力的地形校正、均衡校正、重力场合成等界面起伏的正演计算中，还适用于GRACE时变重力研究中的大气压及水储量的重力效应的正演中。球冠体积分的重力异常正演方法的关键在于球坐标系旋转后新模型的重构。本文将需要重构点的坐标反变换到原坐标系，若原始数据规则，则可以很快找到该点周围的点，进而进行平均或者插值；若数据不规则，则需要先进行网格化成规则模型。相比而言，若先计算好高于分辨率的模型离散数据，进而采用此方案既能满足计算精度又能提高计算效率。
本文通过模型试验与实际应用，表明Tesseroid单元体基于几何中心泰勒级数展开方法的计算精度受级数展开次数、单元体表面大小、厚度、深度、计算点高度与空间地理位置影响，近似的球冠面积分方法的计算精度受模型重构精度、单元体厚度、深度与计算点高度影响，而球冠体积分方法只受原始模型精度的影响。Tesseroid单元体基于几何中心泰勒级数展开方法由于无需进行坐标旋转因而适用于重力三维反演[21]，而球冠体积分方法适用于局部区域或者全球的较高精度的重力或重力异常正演。
| [1] | ZHENG Wei, XU Houze, ZHONG Min, et al. Progress and Present Status of Research on Earth’s Gravitational Field Models [J]. Journal of Geodesy and Geodynamics, 2010, 30(4): 83-91. (郑伟, 许厚泽, 钟敏, 等. 地球重力场模型研究进展和现状[J]. 大地测量与地球动力学, 2010, 30(4): 83-91.) |
| [2] | XIAO Yun, XIA Zheren, WANG Xingtao. Recovering the Earth Gravity Field from Inter-satellite Range-rate of GRACE [J]. Acta Geodaetica et Cartographica Sinica, 2007, 36(1): 19-25. (肖云, 夏哲仁, 王兴涛. 用GRACE星间速度恢复地球重力场[J]. 测绘学报, 2007, 36(1): 19-25.) |
| [3] | YAN Jianguo, LI Fei, PING Jinsong, et al. The Development Situation and Prospect of Lunar Gravity Field Model [J]. Journal of Geomatics, 2010, 35(1): 44-46. (鄢建国, 李斐, 平劲松, 等. 月球重力场模型发展现状及展望[J]. 测绘信息与工程, 2010, 35(1): 44-46.) |
| [4] | YAN Jianguo, LI Fei, PING Jinsong, et al. Lunar Gravity Field Recovery Based on LP Doppler Data [J] . Acta Geodaetica et Cartographica Sinica, 2009, 38(1): 6-11. (鄢建国, 李斐, 平劲松, 等. 利用LP多普勒数据解算月球重力场模型的分析[J]. 测绘学报, 2009, 38(1): 6-11.) |
| [5] | XU Ya, HAO Tianyao. Review on the Progress of Lunar Gravity Field Research and Its Applications [J]. Geochemica, 2010,39(1): 25-31. (徐亚, 郝天姚. 月球重力场研究及其应用进展[J]. 地球化学, 2010, 39(1): 25-31.) |
| [6] | YAN Jianguo, PING Jinsong. A Gravity Field for Mars [J]. Physics, 2009, 38(10): 707-711. (鄢建国, 平劲松. 火星重力场研究现状及发展趋势[J]. 物理, 2009, 38(10): 707-711.) |
| [7] | YAN Jianguo, LI Fei, PING Jinsong. Precision Orbit Determination of MGS Mapping Phase Arcs and Martian Gravity Field Model Solution [J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(5): 484-491. (鄢建国, 李斐, 平劲松. 基于MGS测图段部分弧段的精密定轨及火星重力场模型解算[J]. 测绘学报, 2010, 39(5): 484-491. ) |
| [8] | TAPLEY B D, BETTADPUR S, RIES J C, et al. GRACE Measurements of Mass Variability in the Earth System [J]. Science, 2004, 305(5683): 503-505. |
| [9] | PARKER R L. The Rapid Calculation of Potential Anomalies [J]. Geophysical Journal, 1972, 39: 447-455. |
| [10] | WIECZOREK M A, PHILLIPS R J. Potential Anomalies on a Sphere: Applications to the Thickness of the Lunar Crust [J]. Journal of Geophysical Research ,1998, 103(E1): 1715-1724. |
| [11] | ANDERSON E G. The Effect of Topography on Solutions of Stokes’ Problem [R]. Kensington: University of New South Wales, 1976. |
| [12] | WILD-PEEIFFER F. A Comparison of Different Mass Elements for Use in Gravity Gradiometry [J]. Journal of Geodesy, 2008, 82: 637-653. |
| [13] | ARDALAN A A, SAFARI A. Ellipsoidal Terrain Correction Based on Multi-cylindrical Equal-area Map Projection of the Reference Ellipsoid [J]. Journal of Geodesy, 2004, 78(1-2): 114-123. |
| [14] | HECK B, SEITZ K. A Comparison of Tesseroid, Prism and Point-mass Approches from Mass Reductions in Gravity Field Modelling [J]. Journal of Geodesy, 2007, 81(3): 121-136. |
| [15] | ZHOU Cong, DU Jinsong, LIANG Qing, et al. A Terrain Correction Method for Lunar Gravity in the Domain of the Spherical Cap [J]. Progress in Geophysics, 2010, 25(2): 486-493. (周聪, 杜劲松, 梁青, 等. 基于球冠域的月球重力地形校正方法[J]. 地球物理学进展, 2010, 25(2): 486-493.) |
| [16] | LIANG Q, CHEN C, HUANG Q, et al. Bouguer Gravity Anomaly of the Moon from CE-1 Topography Data: Implications for the Impact Basin Evolution [J]. Science in China Series G: Physics, Mechanics and Astronomy, 2009, 52(12): 1867-1875. |
| [17] | LIU Shiyi. The Utilization of Remnant Spherical Shell Model in the Study of Outer Correction for Gravity Survey[J]. Geophysical and Geochemical Exploration, 1985, 9(1): 25-32. (刘士毅. 利用残球壳模型研究重力测量外部改正方法[J]. 物探与化探, 1985, 9(1): 25-32.) |
| [18] | AN Yulin, ZHANG Minghua, HUANG Jinming, et al. The Computation Scheme and Computation Process for Gravity Correction Values within the Pure Spherical Coordinate System [J]. Geophysical and Geochemical Exploration, 2010, 34(6): 1-9. (安玉林, 张明华, 黄金明, 等. 纯球坐标系内各项重力校正值计算方案和过程[J]. 物探与化探, 2010, 34(6): 1-9.) |
| [19] | VANIEK P, NOVK P, MARTINEC Z. Geoid, Topography, and the Bouguer Plate or Shell [J]. Journal of Geodesy, 2001, 75(4): 210-215. |
| [20] | WANG H S, WU P, WANG Z Y. An Approach for Spherical Harmonic Analysis of Non-smooth Data [J]. Computers and Geosciences, 2006, 32: 1654-1668. |
| [21] | LIANG Qing. Gravity Anomaly Feasures and 3D Density Imaging of the Moon [D]. Wuhan: China University of Geosciences, 2010. (梁青. 月球重力异常特征与三维密度成像研究[D]. 武汉:中国地质大学, 2010.) |
| [22] | PING Jinsong, HUANG Qian, YAN Jianguo, et al. Lunar Topographic Model CLTM_s01 from Chang’E-1 Laser Altimeter [J]. Science in China Series G: Physics, Mechanics and Astronomy, 2008, 52(7): 1105-1114. |
| [23] | DU Jinsong, CHEN Chao, LIANG Qing, et al. The Characteristics of Rock-density Distributions on the Surface and in the Crust of the Moon [J]. Chinese Journal of Geophysics, 2010, 53(9): 2059-2067. (杜劲松, 陈超, 梁青, 等. 月球表层及月壳物质密度分布特征 [J]. 地球物理学报, 2010, 53(9): 2059-2067.) |