169 km以远地壳质量的重力校正值高精度计算及其数值特征

引用本文

安玉林, 郭良辉, 张明华. 169 km以远地壳质量的重力校正值高精度计算及其数值特征[J]. 物探与化探, 2015,39(1): 1-11.
AN Yu-Lin, GUO Liang-Hui, ZHANG Ming-Hua. High precision computation and numerical value characteristics of gravity emendation values arising from mass of the Earth's crust at the distance over 169 km from the observation point[J]. Geophysical and Geochemical Exploration, 2015,39(1): 1-11. 复制到剪切板

Doi:10.11720/wtyht.2015.1.01
Permissions

《物探与化探》编辑部所有

169 km以远地壳质量的重力校正值高精度计算及其数值特征

安玉林¹, 郭良辉¹, 张明华²

1.中国地质大学地球物理与信息技术学院,北京 100083

2.中国地质调查局发展研究中心,北京 100037

作者简介: 安玉林(1941-),男,教授,1982年毕业于武汉地院北京研究生部,主要从事重力和磁力勘探的教学与科研工作. anyulin@cugb.edu.cn

收稿日期: 2013-12-16

基金: SinoProbe专项课题( SinoProbe-02-01,SinoProbe-01-05 ); 国家自然科学基金面上项目( 41374093 ); 北京高等学校青年英才

摘要

以前国内重力勘探教科书中,关于2.0 km以远地壳质量重力校正值的计算仅限于2.0~166.7 km圆形环带以内,并且采用的是直角坐标系内成立的计算公式。近年,中国地质调查局推出直角坐标系公式和球坐标系公式一起应用的重力校正值计算程序,但校正值计算涉及范围仍然局限于2.0~166.7 km圆形环带内。笔者曾推导出球壳型六面体重力场△ g( z_i)公式和其他与重力校正计算相关的公式,现用这些公式开展纯球坐标系内地壳质量的重力校正值高精度计算及其数值特征研究。取得的成果是:①全球陆地和海洋表面、尺度约40 km正方形网格上,169 km以远地壳全部质量重力校正值计算;②中国陆地2'×2'网格上,169 km以远地壳全部质量重力校正值计算;③西藏雅江大转弯3°×2°小区地表、尺度约0.556 km正方形网格上,169 km以远地壳全部质量重力校正值计算。通过对上述全球和局部地区169 km以远地壳质量的重力校正值分布特征分析,得到如下结论:①全球重力校正值的最大值、最小值和平均值分别为106.990×10^-5 m/s²(87.877°E,32.271°N),-41.146×10^-5 m/s²(166.122°E,28.327°N)和-16.439×10^-5 m/s²,其数值分布特征与全球高程/海深分布特征基本一致。②在局部地区,169~1 272 km大环带的地壳质量的陆地地形校正值分布特征与该区高程分布特征基本一致。这说明,在地形高程差异大的地区,重力校正值中存在与地形高程正相关的高频成分,与以前众多专家的认识大不相同。实际上,该高频成分是由计算区本身相邻计算点之间存在较大的高程差值引起的。③无论局部地区及其周围陆高或海深变化多么大,1 272 km以远地壳质量的重力校正值均近似为数值很小的常数,可以不计算。④当局部地区及其周围高程或海深变化均很平缓时,169 km以远地壳全部质量的重力校正值也近似为常数,也可以不计算。此成果对于完善地壳质量重力校正值高精度计算有重要意义。

关键词: 169 km以远; 地壳质量; 重力校正值; 高精度计算; 数值特征

中图分类号:P631.2 文献标志码:A 文章编号:1000-8918(2015)01-0001-011

High precision computation and numerical value characteristics of gravity emendation values arising from mass of the Earth's crust at the distance over 169 km from the observation point

AN Yu-Lin¹, GUO Liang-Hui¹, ZHANG Ming-Hua²

1. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China

2. Development and Research Center of China Geological Survey, Beijing 100037, China

Abstract

In the previous textbooks of gravitational prospecting, the computation of gravity emendation values arising from mass of the Earth's crust at the distance over 2 km from the observation point is confined only to an area of circular ring 2.0~166.7 km in breadth and the computation expressions coming into existence in the orthogonal coordinate system are adopted. China Geological Survey has generalized in recent years the computation program of gravity emendation value that adopts the computation expressions both in the orthogonal coordinate system and in the spherical coordinate system; nevertheless, the computation of gravity emendation values remains confined to an area of circular ring 2.0~166.7 km in breadth and fails to take into account the gravitation arising from mass of the Earth's crust at the distance over 166.7 km from the observation point. 10 years ago, the authors deduced alone gravity expression of the hexahedron with a spherical crust form and other expressions related to the computation of gravity emendation values. After this, adopting these expressions, the authors conducted researches on high precision computation and numerical value characteristics of gravity emendation values arising from mass of the Earth's crust in pure spherical coordinate system, with the following computation achievements obtained: ① high precision computation result of "Gravity Emendation Value arising from Mass of the Earth's Crust at Distance Over 169 km From Observation Point" ( GEVMECDO 169 km FOP ) at the square grid 40 km in breadth on whole Earth's land and ocean surface; ② high precision computation achievements of GEVMECDO 169 km FOP at longitude and latitude grid of China's continent; ③ high precision computation achievements of GEVMECDO 169km FOP at the square grid 0.556km in breadth in local 3°×2° area along the large curved portion of the Yarlung Zang Zangbo River. Based on analyzing numerical values distribution characteristics of GEVMECDO 169 km FOP on the whole Earth's and local surfaces, the authors have arrived at the following conclusions d: ① The maximum, minimum and mean values of gravity emendation values on the whole Earth's surface are 106.990 mgal (87.877°E, 32.271°N), -41.146 mgal (166.122°E,28.327°N) and -16.439 mGal respectively; the distribution characteristics of these values are in the main consistent with the distribution characteristics of the Earth's land altitude/sea depth values. ② The distribution characteristics of gravity terrain emendation values in a local area arising from mass of the Earth's crust in large circular rings 169-1 272 km in breadth are in the main consistent with the distribution characteristics of land altitude values in this local area. These data suggest that there are high frequency components in gravity emendation values of the area with larger altitude differences,and these high frequency components with terrain altitude values have positive correlation. These facts differ remarkably from the opinions put forward by many researchers. In fact, these high frequency components result from larger altitude differences of neighboring points in a local area. ③ No mater how large the differences of altitudes /sea depths in a local area and its peripheries are, these GEVMCDO 1 272 km FOP are approximately close to a very small constant and may not be computed; ④ Where the altitude or sea depth of local area and its peripheries changes smoothly, the gravity emendation values of the whole mass of the crust over the distance of 169km are approximately close to a constant and hence may not be calculated. The results achieved by the authors have important significance for perfecting high precision computation of gravity emendation values arising from mass of the Earth's crust. The authors suggest introducing these achievements into new textbooks of gravitational prospecting and also propose to study and extend this research topic with high precision and high speed with the purpose of convenient application of gravity emendation values arising from mass of the Earth's crust on the basis of these achievements.

Keyword: distance over 169 km; mass of the Earth's crust; gravity emendation value; high precision computation; numerical characteristics

Show Figures

笔者阐述的是“ 计算点处2.0 km以远全球地壳质量引起的重力校正值的高精度计算方法” 和“ 计算点外2.0 km以远地壳质量的4大环带区的各自重力校正值的数值特征” , 采用的计算重力校正值的“ 地壳质量单元模型” 是由“ 地球的余纬θ ₁与θ ₂、经度λ ₁与λ ₂和径向半径R₁与R₂” 等6个参数确定的。这6个参数确定了6个扇形面, 共同界定了一个“ 地壳质量单元模型” 。这种单元模型可以形成无缝对接组合, 高精度地模拟全球地壳质量分布。

我国陆地上山地多, 平原少, 在重力找矿过程中, 专家们非常重视对实测重力值进行各项校正计算。除了计算2.0~166.7 km范围内地壳质量重力校正值外, 早在20世纪80年代, 就关注到“ 广义重力校正” , 接受“ 格莱尼异常” 概念, 认为有必要计算出全球地壳质量的重力校正值^{[1, 2, 3, 4]}。为了进行“ 广义重力校正” , 文献[1-4]采用的“ 地壳质量单元模型” 尽管名称不同, 实际上就是“ 常密度球壳型六面体” 。文献[1]最先刊出该模型在测点i的重力场△ g(z_i)解析表达式; 雷受旻^[2]、陈凤英^[3]、杨学样^[4]等分别独立推导了该表达式; 雷受旻^[2]、刘士毅^[3]、杨学样^{[4, 5]}等还阐述了运用该表达式进行“ 广义重力校正” 的多种技术问题。为了减小全球地壳质量重力地形校正值计算量, 刘立言等^[6]设计计算了重力地形改正(166.7 km至全球)改算图, 供技术人员应用; 陈邦彦^[7]亲自主持计算了南海中央海盆格莱尼重力异常, 探讨了该异常的应用前景。

但由于多种原因, 全球重力校正值计算并没有被广泛采用, 甚至连“ 常密度球壳型六面体” 的Δ g(z_i)解析表达式及其较为概略的推导过程, 也没能编入国内“ 重力勘探” 教科书中, 致使其鲜为人知。

21世纪初, 重力资料的定量处理和解释技术研究得到加强。2002年, 柴玉璞承担了原石油天然气集团总公司物探局下达的“ 复杂地形条件下重力场外部校正技术研究” 项目。安玉林作为协作者, 与柴玉璞一起, 详细推导了与“ 常密度球壳型六面体” Δ g(z_i)表达式相关的所有公式。2006年, 在国家油气专项“ 青藏高原油气资源战略选区调查与评价” 中, 张明华承担了“ 青藏高原比如县、可可西里盆地重力资料处理” 子专题, 作为协作者, 安玉林采用亲自推导的公式, 研制了“ 2.0~169.0 km 地壳质量高精度地形校正和均衡校正值计算程序” , 并用于比如、可可西里, 乃至青藏高原N36° 以南重力资料校正中^[8]。

近年来, 以陈超为首的科研集体参与国家嫦娥绕月探测工程, 开展全月月壳质量重力校正值计算、重力异常圆滑处理和密度横向变化反演等研究, 取得了一系列成果^{[9, 10, 11, 12, 13, 14]}。他们计算全月重力异常校正值采用的Tesseroid单元模型^[15]及其Δ g(z_i)解析表达式, 就是“ 常密度球壳型六面体” 模型及其△ g(z_i)解析表达式, 并认为依据该公式研制的重力校正值计算方法具有较高精度。为了提高计算效率, 他们还给出该模型体的简化的△ g(z_i)表达式^[12]和双线性内插出需要重构点的坐标的方法^[14]。

2012年以来, 郭良辉承担了 SinoProbe专项课题( SinoProbe-02-01, SinoProbe-01-05)、国家自然科学基金面上项目(41374093)、北京高等学校青年英才计划和中国地质大学(北京)“ 重力全球地形改正实验” 等课题, 本文阐述的就是上述课题研究成果。

目前普通计算机硬件发展达到相当高的水平, 可提供巨大存储空间和高速计算能力。故此, 实现“ 高精度计算” 应该成为全球重力校正值计算程序研制的首要课题。

笔者研制的“ 计算点处169.0 km以远全球地壳质量引起的重力校正值的计算方法” 之所以具有高精度, 是因为采取了5种技术方案。

(1)采用“ 常密度球壳型六面体” 作为单体模型, 用其△ g(z_i)解析表达式作为地壳质量引起的重力校正值的正演计算公式^[8]。该模型优于其他单体模型^[15], 可以形成无缝对接的组合模型, 精确模拟地壳质量分布, 而且其正演表达式严格精确。

(2)以重力测点为北极点, 建立以其为中心的新的球面经纬度坐标系, 实现了纯球坐标系中的重力改正值计算。而且, 球面上点的原经纬度坐标向新经纬度坐标转换的公式^[8]严格采用球面三角公式。

(3)采用文献[8]提出的以重力测点为中心的球面环半径设计公式, 把2.0 km以远球面, 划分为宽度有规律增加的107个小环带。进而把107个小环带组合为16个二级环带, 4个二级环带为一组再合成4个一级环带。一级环带分别是测点外2.0→ 169.0 km, 169.0→ 1 271.5 km, 1 271.5→ 5 938.9 km和5 938.9→ 2 0354.7 km。2.0→ 169.0 km内45个小环均分为36个扇形, 总共1 620个扇形单元。而169 km以远的12个二级环带, 从内到外, 分别分割为36、40、45、50、55、60、60、60、52、40、24、4个扇形(见4.2节表1)。地壳质量如此剖分, 精细地表示出了其空间分布。

(4)对应于169 km以远的12个二级环带, 提取“ 全球12组正方形网格结点上高程/海深值” (见3.2节表1), 参与每个测点重力校正值计算, 保证全球(无论是靠近两极, 还是在赤道附近的)测点的重力校正值具有几乎相同的计算精度。

(5)采用测点的实际高程值和落在每个“ 球壳形六面体” 顶面扇形内的3~8个高程/海深值的平均值, 计算测点重力校正值; 而不是预先计算出球面上选定结点处的重力校正值, 然后采用双线性内插方法, 内插出测点处的重力校正值^[16]。

采用上述方案实现的全球地壳质量重力校正值“ 高精度” 计算过程, 其运算速度也相当之快, 原因如下:①第2技术方案中, 采用的新、旧球坐标间转换公式简洁, 计算用时极少; ②第4技术方案可以大大减少求取“ 球壳型六面体” 径向厚度平均值所需的高程/海深值的搜索量; ③针对全球面上、不同位置和范围大小的重力测区的一般情况, 研制了普遍适用的二种“ 校正值计算前结点数缩减程序” , 用于剔除“ 12个正方形网格结点上计算中用不到的高程、海深值” , 大大减少了计算所需的搜索量; ④抛弃仅在陆地上方便采用、但不适合于全球的中间层校正值计算过程; ⑤将大计算区划分为n个小计算区, 在多线程计算机上进行n个单线程程序并行计算, 最高可提速近n倍。

在实现全球地壳质量重力校正值高精度计算基础上, 笔者研究了计算点外169 km以远地壳质量的3个一级环带区各自重力校正值的数值特征。

以前, 即使采用球坐标系地壳质量模型计算重力校正值, 也仅仅算到166.7 km^[16], 因为多数专家认为, 166.7 km以远地壳质量在重力测区内产生的重力值呈现光滑低缓背景异常特征, 可借助高通滤波加以消除。然而, 在高程/海深值变化较大的重力测区, 这种认识是错误的。在柴玉璞与安玉林协作开展的研究中, 张家界重力测区重力校正实践^[17]表明, 在高程/海深值变化较大的重力测区, 166.7 km以远地壳质量的重力校正值具有显著的高频成分(由测区内高程/海深值急剧变化引起的), 其幅度大到严重影响异常形态的程度。此次的研究成果再次证实了这一结论。

全球地壳质量高精度重力校正值计算, 可以有效消除实测重力值内包含的低缓背景和高频成分, 提取出高精度重力异常。如果计算出整个中国0.5 km× 0.5 km网格上、2 km以远地壳质量的高精度重力校正值, 供重力勘探专业科研和生产单位调用, 必将产生可观的经济效益。

1 与地壳质量分布有关的重力校正值计算

按照爱黎地壳均衡理论, 绘制如下与地壳质量分布有关的重力校正值计算示意图(图1)。

	Figure Option View Download New Window
	图1 观测点P处地壳质量重力校正值计算示意

首先说明, 图1中, 代表地壳平均厚度的30 km等厚层H₀应该是等厚球壳形状, 为了方便, 这里画成了水平层形状, 但并不影响对问题的阐述和理解。相对应的, 陆地高程线、海底深度线、与陆地高程和海底深度关联的均衡补偿线等, 也均画成依附于水平线的形状。

与地壳质量分布有关的校正值, 是指陆地多余质量校正值、海洋亏损质量校正值、陆地下补偿质量校正值和海洋下补偿质量校正值等四种, 前二种可称为地形校正, 后两种可称为均衡校正。研制实际重力校正值计算程序时, 地形校正分为传统意义的“ 地形校正” 和“ 水体校正” 。

图1中, R₀线近似代表地球水准面, R₀是该面与地心的距离; R₀-H₀线代表平均厚度地壳的球形下界面。等厚层H₀的密度ρ =2.67 g/cm³。相对应的, 陆地多余质量密度Δ ρ =2.67 g/cm³; 海洋亏损质量密度Δ ρ =1.64 g/cm³; 陆地下补偿质量Δ ρ =-0.60 g/cm³, 深度补偿系数t=4.45; 海洋下补偿质量密度Δ ρ =+0.60 g/cm³, 深度补偿系数t=2.73。

假设计算点P_i分布位置是: 在陆地, P_i位于地形线上; 在海洋, P_i位于海平面(R₀)上。

图中4个矩形柱体代表4个球壳型六面体。h_AV代表单个球壳型六面体顶面包含的多个地形高程的平均值或多个海底深度的平均值, R₁和R₂代表单个球壳型六面体的两个小球面。位于上面的两个“ 六面体” , R₁=R₀, R₂=R₀± h_AV; 位于下面的两个“ 六面体” , R₁=R₀-H₀, R₂=R₁± t_AV。对应陆地下低密度补偿层, t_AV=4.45h_AV, 对应海洋下高密度补偿层, t_AV=2.73h_AV。

经过严格分析, 得到:总校正值=(水体校正值绝对值+低密度补偿层校正值绝对值 )-(陆地地层校正值绝对值+高密度补偿层校正值绝对值)。特别指出:①在图1和总校正值表达式中, 抛弃了以前被广泛采用的“ 中间层校正” 环节, 而是直接计算每个常密度球壳型六面体的重力校正值; ②文内的“ 重力地形校正值” 仅指陆地地壳质量产生的校正值。

2 纯球坐标系内地壳质量重力校正值计算方案和计算公式

与地壳质量分布有关的纯球坐标系内重力校正值计算方案中, 每个重力校正值计算点P_i都要被当成北极点, 而每个场源点Q_i原来的经纬度坐标, 都要根据球面三角公式, 变换成以P_i为北极的新的经纬度坐标。

以计算2.0~169.0 km 的大环带区重力校正值^[8]时所采用的计算方案为例(图2)。

	Figure Option View Download New Window
	图2 计算2~168.971 km 的重力校正值的新方案

该示意图中心点为P_i点, 围绕该点有许多同心圆, 又有16条从该点出发的、22.5° 等夹角射线。计算同心圆半径(即环半径)a_m值的公式为

$\begin{matrix} a_{m} = a_{m - 1} + Δa \cdot b^{m - 1}, m = 2, 3, \dots, M, \end{matrix}$ （1）

式中:a₁=2 000 m, Δ a=400 m, b=1.08, 环数M=46。这套参数在西藏雅江大转弯计算区计算的2.0~168.971 km环区地壳质量产生的重力校正值, 精度相当高。

该计算方案由M个环和N个(图2是16个)扇组成。在两个球面同心圆和两条射线之间的“ 场源体” 称为球壳型六面体。该六面体由6个扇形截面限定, 分别对应参数θ ₁、θ ₂、λ ₁、λ ₂、R₁、R₂, 其中, θ ₁、θ ₂、λ ₁、λ ₂是以P_i点为北极点的新的余纬坐标和经度坐标, 而夹角Δ λ =λ ₂-λ ₁=360/N。

为了适应地形起伏巨大地区重力校正值的高精度计算, 式(1)的基本参数定为:a₁=2 000 m, Δ a=400 m, b=1.08, M=108。

	Figure Option View Download New Window
	图3 球面三角公式推导方法示意

为了使每一个重力测点P_i均成为“ 北极点” , 需要计算出每一个地形结点Q_i相对于该“ 北极点” 的新的经度和纬度。为此, 按照图3, 采用文献[18]中的球面三角形余弦定理、正弦定理和五元素公式, 可推导出“ 球面三角公式” 的边(α )、角(β )为:

$\begin{matrix} \begin{matrix} α = \arccos [sinφ \cdot \sin φ_{0} + cosφ \cdot \cos φ_{0} \cdot \cos (λ - λ_{0})], (2) \\ β = \arctan [\frac{cosφ \cdot \sin (λ - λ_{0})}{sinφ \cdot \cos φ_{0} - cosφ \cdot sinφ \cdot \cos (λ - λ_{0})}], (3) \end{matrix} \end{matrix}$

其中, 式(2)直接由余弦定理得到。而由正弦定理和五元素公式, 立即得到sinβ = $\begin{matrix} \frac{sinb \cdot sinα}{sina} \end{matrix}$ 和cosβ = $\begin{matrix} \frac{cosb \cdot sinc - sinb \cdot cosc \cdot cosα}{sina} \end{matrix}$ , 继而求出tanβ = $\begin{matrix} \frac{sinβ}{cosβ} \end{matrix}$ , 取反正切即得式(3)。当然, 上面二式也可仅由余弦定理和五元素公式导出。

公式中, (λ ₀, φ ₀)和(λ , φ )分别是重力测点P_i和场源点Q_i的地理经度和纬度; λ -λ ₀是Q_i点与P_i的经度差。a单位若取“ ° ” , 它相当于后面式(6)中的余纬θ , 而(90° -a)是Q_i点相对P_i点的“ 纬度” φ ; a单位若取“ 弧度” , 则a× R₀ 就是P_i与Q_i的球面距离。β 单位若取“ ° ” , 它是以地理北极点N、P_i点和Q_i点为顶点的球面三角形的一个角, 该角顶点在P_i, 边为P_iN和P_iQ_i; 计算中, 以P_iN为“ 0° 线” , 则β 就是Q_i相对P_i点的“ 经度” ; 如果λ -λ ₀> 0, 则β 取值为0° ~180° , 如果λ -λ ₀< 0, 则β 取值为0° ~-180° (即180° ~360° ); 可见, 这里的“ 经度” 是顺时针增加, 与地理经度反时针增加相反。

图4中, z_i=R₀+h_i, R₀为地球平均半径, h_i为i点高程/海深值。常密度球壳型六面体的Δ g(z_i)解析表达式推导步骤简述如下。

(1)见图4, 首先推导出“ 球壳型等厚层(即球冠)” Δ g(z_i)解析表达式, 其初始公式是

$\begin{matrix} Δg (z_{i}) = Gρ \overset{2 π}{\int_{0}} \overset{θ}{\int_{0}} \overset{R_{2}}{\int_{R_{1}}} \frac{(z_{i} - rcosθ)}{{r^{2} - 2 z_{i} rcosθ + (z_{i})^{2}}^{\frac{3}{2}}} \cdot r^{2} sinθdrdθdλ \end{matrix}$ (4)

	Figure Option View Download New Window
	图4 “ 球壳型等厚层” 重力场Δ g(z_i)表达式推导示意

采用文献[18]中的积分公式, 经过较长推导过程, 可求得“ 球壳型等厚层” Δ g(z_i)解析表达式。该表达式就包含在式(6)中, 即:在该式中, 令θ ₂=θ 、θ ₁=0, 再去掉系数(1/N)所得到的公式。

(2)由余纬角分别为θ ₂和θ ₁(θ ₂> θ ₁) 的两个“ 球壳型等厚层” 解析表达式相减, 可求得“ 球壳型等厚同心园环” 的Δ g(z_i)解析表达式

$\begin{matrix} △ g (z_{i})_{球壳型等厚同心圆环} = Δg (z_{i}, θ_{2}) - Δg (z_{i}, θ_{1}) \end{matrix}$ （5）

式(5)就是式(6)去掉系数(1/N)所得到的表达式。

(3)把“ 球壳型等厚同心圆环” 分为N等份, 得到N个球壳型六面体, 它们在计算点P(0, 0, z_i)的重力值相等, 故该“ 六面体” 在P(0, 0, z_i)点的重力Δ g(z_i)解析表达式为^[8]:

$\begin{matrix} \begin{matrix} Δg (z_{i}) = - V (z_{i}) = \frac{1}{N} \cdot \frac{2}{3} πGρ \cdot z_{i} \{[{(\frac{r}{z_{i}})}^{2} + \frac{r}{z_{i}} \cdot cosθ + 3 co s^{2} θ - 2] \sqrt[]{{(\frac{r}{z_{i}})}^{2} - 2 \cdot \frac{r}{z_{i}} \cdot cosθ + 1} - \\ {3 \cdot cosθ \cdot si n^{2} θ \cdot \ln [\frac{r}{z_{i}} - cosθ + \sqrt[]{{(\frac{r}{z_{i}})}^{2} - 2 \frac{r}{z_{i}} cosθ + 1}]|}_{r = R_{1}, θ = θ_{1}}^{r = R_{2}, θ = θ_{2}} = \\ Δg (z_{i}, r = R_{2}, θ_{2}, θ_{1}) - Δg (z_{i}, r = R_{1}, θ_{2}, θ_{1}) \end{matrix} \end{matrix}$ （6）

本论文阐述的重点是, 169 km 以远全球地壳质量重力校正值计算程序研究、全球和局部地区校正值计算及不同二级环带校正值的数值特征分析。上面式(1)、(2)、(3)、(6)是本研究所采用的最为重要的计算公式。

为了保证计算结果的高精度和所得结论的可靠性, 本次研究方案中采用的环半径参数值仍然是a₁=2 000 m, Δ a=400 m, b=1.08, M=108; 所涉及环号为46~108, 共62个小环带, 12个二级环带。

3 实施高精度和高速度计算的技术方案

3.1 12个正方形网格结点数组提取方法

为了适应“ 引言” 中所述高精度计算方案(3)和方案(4)的要求, 提高计算速度, 对应169 km以远的12个二级环带, 地形高程/海深结点数组的正方形网格密度并不相同。由内到外, 二级环带所需网格密度应该递减。经试算确定, 沿经度线和纬度线的网格间距(球面距离)分别是10、20、30、40、60、90、130、180、250、350、480、660 km。确定上述正方形网格点经、纬度坐标的方法如下。

(1)设沿经度线的球面距离Δ l(如10 km), 则有纬度间隔Δ φ =180.0× (Δ l/R₀)/π (单位为“ ° ” )。故纬度点

$\begin{matrix} \begin{matrix} φ_{m - 1} = 90 - (m - 1) \cdot Δφ, \\ m = 2, 3, \dots, inT [k (90 / Δφ) + 1] \end{matrix} \end{matrix}$ （7）

式中, in代表商数取整(即舍去商数的小数部分取整数值, 下同)。

(2)设沿纬度线φ _m的球面距离也是Δ l(如10 km), 其切割园半径r=R₀· cosφ _m, 则切割园上经度间隔Δ λ _m=180.0× (Δ l/Δ r_m)/π (单位也是“ ° ” )。故经度点

$\begin{matrix} \begin{matrix} Δ λ_{m (n - 1)} = (- 180.0) + (n - 1) \cdot Δ λ_{m}, \\ m = 2, 3, \dots, inT [k (360 / Δ λ_{m}) + 1] \end{matrix} \end{matrix}$ （8）

按照这种算法得到的正方形网格, 其网格面积非常接近(Δ l)²。

采用式(7)、(8)编制了“ 全球正方形网格结点取数程序” , 从全球高程/海深值高密度数据中, 经内插取出了12个不同尺度正方形网格场源结点数组, 各数组均包括所有点的地理经度、纬度、陆高/海深值等3列数据。这些数组既可以供全球重力校正值计算使用, 也可以供地球任何一个局部地区重力校正值计算使用。

3.2 高精度计算方案综合数据表

如前文所述, 把计算每个测点重力校正值涉及的169 km以远地壳质量的46~108环, 分为3个一级环带, 每个一级环带又分为4个二级环带, 每个二级环带均包括多个小环。表1给出一级环号、二级环号和所有小环的环号、环半径、扇数、扇宽、陆高/海深值数组、网格尺度和数据个数等。由该表可以算出, 计算1个点的重力校正值所涉及的“ 球壳型六面体” 的个数为5 610。

表1 高精度计算方案综合数据

3.3 方形网格结点数的算前缩减方法

表1中第一大环第一个二级环带的方形网格结点数约510万(5 103 980)之多, 致使求取该二级环带内“ 球壳型六面体” 平均厚度的搜索量非常大。当计算区范围有限或分割成较小计算区时, 远离算区的绝大多数结点值并不参与求平均计算, 故应该结合具体算区大小, 在重力校正值正式计算之前, 剔除多余结点值, 以便进一步提高计算速度。为此, 针对该二级环带, 采取如下两种技术方案。

假设有限算区内所有计算点P的最小纬度、最大纬度、最小经度、最大经度分别为PW_min 、PW_max、PJ_min(算区最左端)、PJ_max(算区最右端)。

第一方案是“ 球面圆形区缩减法” :①求出中值PW00=(PW_min+PW_max)/2 和PJ00=(PJ_min+PJ_max)/2; ②采用式(2)分别计算出点(PJ00, PW00)与点(PJ_min, PW_min)和点(PJ_max, PW_max)之间的球面距离RP_min和RP_max, 并选取二者最大值 R_max=MAX(RP_min, RP_max); ③以公式R=R_max+300 km(式中300 km是一略大于该二级环带最大环半径292.013 km, 而又方便使用的整数)确定搜索半径, 并以(PJ00, PW00)点为中心确定搜索范围; ④再次采用式(2), 分别求取5 103 980个结点Q与(PJ00, PW00)点的距离QP; ⑤如果某一结点的QP> R, 则“ 剔除该结点” , 否则, 保留它, 并在已经被保留结点的数量KK上加1, 如此重复, 直至搜索完第5 103 980个结点Q为止, 记下KK值, 它是缩减了的结点数组的数据个数。

第二方案是“ 球面扇形区缩减法” :①求Δ θ =180.0× (300 km/R₀)/π , 单位为“ ° ” ; ②如果所有结点中某结点Q的纬度Q_ω 满足条件{Qω < PW_min-Δ θ -0.1⁰}或者{Qω > PW_max+Δ θ +0.1⁰}, 则“ 剔除该结点” , 否则, 保留它。此二步完成了沿南北方向“ 剔除结点” 计算。此后, 再进一步完成沿东西方向“ 剔除结点” 计算。此算法包含了第一方案中的③、④、⑤、⑥4个步骤, 它们涉及计算区东西方向的最左端和最右端的计算点经度以及二者中心经度等与地球东、西半球分界线(即东经180° 与西经180° 重合线)的位置关系。

上面二方案中, 第一方案保留的结点数多于第二方案。我们研制了包含12个结点数组和两个算前缩减方案的统一程序, 使用该程序大大提高了重力校正值计算速度。

4 全球和中国陆地重力校正值计算成果图及其数值特征

4.1 全球重力校正值计算成果图及其数值特征

计算点个数为320 358, 则需要计算的“ 球壳型六面体” 的总个数有17亿之多(5 610× 320 358)。由于一次计算320 358个点所需时间较长, 故分为10个分段计算(表2)。一台笔记本电脑可以同时计算3~4个分段, 以便提高速度。

表2 全球320358个重力校正值分段计算

根据式(1)、(2)、(3)、(6)和表2, 研制了用于全球重力校正值计算的Fortran程序。全球重力校正值计算成果见图5。

图5中显示的重力校正值的分布特征是:绝大部分陆地为正值, 高程越大, 正值越大, 中国青藏高原处最大; 绝大部分海洋为负值, 深度越大, 负值绝对值越大。这一特征, 经与全球高程和海深等值线图对比, 一目了然。全球重力校正值的最大值、最小值和平均值分别为106.990× 10^-5 m/s²(87.877° E, 32.271° N)、-41.146× 10^-5 m/s²(166.122° E, 28.327° N)和-16.439× 10^-5 m/s²。

	Figure Option View Download New Window
	图5 全球约40 km方形网格点上169 km以远地壳质量重力校正值等值线

4.2 中国陆地重力校正值计算成果及其数值特征

图6~图9为全国陆地经纬度2'× 2'网格上地壳质量高精度重力校正值计算成果中的几个重要图件。

图6可见, 169~1 272 km一级环带内, 浅于大地水准面的陆地地壳质量的高精度重力地形校正值与陆地地形高程呈现明显的正相关关系, 该特征在青藏高原表现得最为显著。

图7可见, 169 km 以远地壳质量高精度重力总校正值与地形高程宏观分布呈正相关关系。其中由地形高程差异引起的重力校正值的高频成分被掩盖于大的总校正值中, 不如图8明显, 但是, 这一特征与图5的特征很相似。

图8可见, 1 272 km 以远地壳质量的“ 重力校正值” 呈现出“ 西北高东南低的十分低缓的倾斜面特征” , 其最大值仅为2.854× 10^-5 m/s²(是图7中最大值107.153× 10^-5 m/s²的1/37), 最大差值是3.328× 10^-5 m/s²。但是, 由地形高程差异引起的重力校正值的高频成分依然清晰可见。

	Figure Option View Download New Window
	图6 169~1 272 km环带地壳质量在全国陆地2'× 2'经纬网格上的重力地形校正值等值线

	Figure Option View Download New Window
	图7 169 km以远地壳质量在全国陆地范围2'× 2'网格上的重力总校正值等值线

	Figure Option View Download New Window
	图8 1 272 km以远地壳质量在全国陆地范围2'× 2'网格上的重力校正值等值线

	Figure Option View Download New Window
	图9 5 939 km以远地壳质量在全国陆地范围 2'× 2'网格上的重力校正值等值线图

图9可见, 5 939 km以远地壳质量的“ 重力校正值” 几乎是一个平均值为0.273× 10^-5 m/s²的平面(其最大差值是0.056× 10^-5 m/s²), 其值可以忽略。但是, 正是因为数值微小, 其中由地形高程差异引起的重力校正值的高频成分更加清晰可见。

由图6~图9给出的高精度重力校正值数值特征, 能够得到下述结论。

(1)完全可以不计算5 939 km 以远地壳质量引起的几乎是微小平面状的重力校正值。

(2)中国陆地范围如此之大, 1 272 km以远地壳质量重力校正值等值线图依然如此低缓。那么, 更小的区域中, 重力校正值之间的最大差值必定更小, 等值线图必定更平缓(见图12), 故也可以不加计算。

(3)地形高程差异大的地区, 重力校正值中总是存在高频成分, 并与地形高程值呈正相关关系。这一数值特征, 与以前众多专家的认识, 即认为“ 166.7 km以远地壳质量在重力测区内产生的重力值呈现光滑低缓背景异常特征, 可借助高通滤波加以消除” , 大不相同。

实际上, 重力校正值中的高频成分, 正是由计算区内相邻计算点之间存在较大的高程差值所引起。在测区重力观测过程中, 当地形起伏较大时, 相邻测点的高程值差异也必然较大, 而观测值中必然既含169 km以内, 也含169 km以外地壳质量引起的重力高频观测数值, 因此, 需要借助计算全球地壳质量的重力校正值, 加以消除。否则, 将在重力异常中存有超出观测精度的较大误差, 对异常形态产生严重干扰, 在测点距离较大时, 有可能形成较低频的虚假异常。

5 西藏雅江大转弯地区重力校正值计算成果图及其数值特征

西藏雅江(93° E~96° E, 29° N~31° N]大转弯地区的地形变化非常剧烈。计算了该区0.005 0° × 0.005 0° 精细经纬度网格上240 001点的重力校正值。该区重力校正值计算成果如图10~图13所示。

图10中, 169~1 272 km 一级环区地壳质量的地形校正值的变化达16.08× 10^-5 m/s², 其等值线特征完全反映了西藏雅江大转弯地区地形高程变化和雅江及其支流水系的分布特征。

图11重力校正值的局部变化特征是图10地形校正值数值变化的反应。由图10和图11可以得出4.2节中第(3)条结论。

图12中大、小极值之差很小(0.151 mgal), 故可不计算1 272 km以远的地壳质量的重力校正值。这进一步支持了4.2节中第(2)条结论。

图13重力校正值等于图11与图12之和, 它几乎是图11的平行曲面。

	Figure Option View Download New Window
	图10 西藏雅江大转弯169~1 272 km大环地壳质量地形校正值等值线

	Figure Option View Download New Window
	图11 西藏雅江大转弯169~1 272 km大环地壳质量重力校正值等值线

	Figure Option View Download New Window
	图12 西藏雅江大转弯1 272 km以远地壳质量重力校正值等值线

	Figure Option View Download New Window
	图13 西藏雅江大转弯169 km以远地壳质量重力校正值等值线

6 结论和建议

(1)全球重力校正值的最大值、最小值和平均值分别为106.990× 10^-5 m/s²(87.877° E, 32.271° N), -41.146× 10^-5 m/s²(166.122° E, 28. 327° N) 和-16.439× 10^-5 m/s²; 其数值分布特征与全球陆高/海深值分布特征基本一致。

(2)在局部地区内, 169~1 272 km一级环带的地壳质量“ 重力地形校正值” 的分布特征, 与该区高程值分布特征基本一致。这说明, 在地形高程差异大的地区, 重力校正值中存在与地形高程正相关的高频成分。实际上, 该高频成分是由计算区本身相邻计算点之间存在较大的高程差值所引起。

在测区重力观测过程中, 当地形起伏较大时, 相邻测点的高程值差异也必然较大; 而观测值中必然既含169 km以内, 也含169 km以外地壳质量引起的重力高频观测数值; 因此, 需要借助计算全球地壳质量的重力校正值, 加以消除。否则, 将在重力异常中存有超出观测精度的较大误差, 对异常形态产生严重干扰, 甚至出现虚假异常。

(3)无论局部地区及其周围陆高/海深变化多么大, 1 272 km以远地壳质量的重力校正值均近似为数值很小的常数, 可以不用计算。

(4)当局部地区及其周围陆高/海深变化均很平缓时, 169 km以远地壳质量的重力校正值也近似为常数, 也可以不计算。

我们一贯认为, 计算精度是科学计算程序的“ 生命” , 完善纯球坐标系内重力校正值计算方法和程序时, 应该把保证“ 高精度” 放在第一位。当前, 一般计算机已经具有多核, 其硬盘空间达到1T(1 024 G), 内存达4G或更多, 单个CPU速度超过2 G/s, 可见, 计算机硬件已经基本解决了存储空间和计算速度问题。另外, 计算机语言和数学计算方法水平很高, 又可采用双精度数据运算, 这又为实现高精度计算提供了保证。因此, 以前不能实现的科学计算和推演的问题, 现在已经能够解决。理论、方法和实用技术研究不应该再满足于传统做法, 不能再采用明显不符合理论要求的技术软件。

基于上述认识, 强烈建议:在重磁资料处理、反演和解释领域, 大力支持和发展能在纯球坐标系内应用的理论方法和实用化程序。如果能计算出整个中国陆地和领海的0.5 km× 0.5 km网格上、2 km以远地壳质量的高精度重力校正值, 供重力勘探专业科研和生产单位调用, 必将产生可观的经济效益。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	DZ/T0082 -1993. 中国地质矿产部. 区域重力调查技术规定[S]. [本文引用:1]
[2]	雷受旻. 重力广义地形改正值和均衡改正值的一种计算方法[J]. 海洋地质与第四纪地质, 1884, 4(4): 101-111. [本文引用:3]
[3]	刘士毅. 利用残球壳模型研究重力测量外部改正方法[J]. 物探与化探, 1985, 9(1) : 25-32. [本文引用:3]
[4]	杨学样. 关于确定山区区域重力测量中间层校正系数的几点意见[J]. 物探与化探, 1990, 14(5): 396-400. [本文引用:3]
[5]	杨学样. 布格改正和地形改正的误差—关于区域重力测量中地形改正最大半径的讨论[J]. 地壳形变与地震, 1992, 12(2): 1-6. [本文引用:1]
[6]	刘立言, 刘定和. 重力地形改正(166. 7 km至全球)改算图[J]. 物探与化探, 1986, 10(5): 381-383. [本文引用:1]
[7]	陈邦彦. 南海中央海盆格莱尼重力异常[C]//中国地球物理学会第六届学术年会论文集, 1990: 307. [本文引用:1]
[8]	安玉林, 张明华, 黄金明, 等. 纯球坐标系内各项重力校正值计算方案和计算过程[J]. 物探与化探, 2010, 34(6): 697-705. [本文引用:5]
[9]	梁青, 陈超. 基于嫦娥一号地形数据的月球布格重力异常与撞击盆地演化[J]. 中国科学 G 辑: 物理学, 力学, 天文学, 2009, 39(10): 1379-1386. [本文引用:1]
[10]	周聪, 杜劲松, 梁青, 等. 基于球冠域的月球重力地形校正方法[J]. 地球物理学进展, 2010, 25(2): 486-493. [本文引用:1]
[11]	杜劲松, 陈超, 梁青, 等. 月球表层及月壳物质密度分布特征[J]. 地球物理学报, 2012, 53(9): 2059-2067. [本文引用:1]
[12]	杜劲松, 陈超, 梁青, 等. 球冠体积分的重力异常正演方法及其与 Tesseroid 单元体泰勒级数展开方法的比较[J]. 测绘学报, 2012, 41(3): 339-346. [本文引用:2]
[13]	杜劲松, 陈超, 梁青, 等. 基于球冠滑动平均的球面重力异常分离方法[J]. 武汉大学学报: 信息科学版, 2012, 37(7): 864-868. [本文引用:1]
[14]	杜劲松, 陈超, 梁青, 等. 月球重力异常及其计算方法[J]. 武汉大学学报: 信息科学版, 2012, 37(11): 1369-1373, 1385 [本文引用:2]
[15]	Heck B, Seitz K. A Comparison of tesseroid, prismand point-mass approches from mass reductions in gravity field modelling[J]. Journal of Geodesy, 2007, 81(3): 121-136. [本文引用:2]
[16]	DZ-T0082 -2006. 区域重力调查规范: 附录G[S]. [本文引用:2]
[17]	Yang Z J, Chai Y P. Spherical external gravity correction[J]. Applied Geophysics, 2005, 2(3): 131-134. [本文引用:1]
[18]	《数学手册》编写组. 数学手册[M]. 北京: 人民教育出版社, 1979: 48-49, 254, 263. [本文引用:1]

0082

0.0

... 概念,认为有必要计算出全球地壳质量的重力校正值^[1,2,3,4] ...

1884

0.0

... 概念,认为有必要计算出全球地壳质量的重力校正值^[1,2,3,4] ...

... 雷受旻^[2]、陈凤英^[3]、杨学样^[4]等分别独立推导了该表达式 ...

... 雷受旻^[2]、刘士毅^[3]、杨学样^[4,5]等还阐述了运用该表达式进行#cod#x0201c ...

1985

0.0

刘士毅. 利用残球壳模型研究重力测量外部改正方法[J]. 物探与化探, 1985, 9(1) : 25-32.

Gravitational formulae of two remnant spherical shells are given in this paper and are used in stone slab correction and topographic correction for gravity survey. The results of the correction have been compared with the corresponding values calculated by finite and infinite plane formulae commonly used in the past, which shows that, in using the plane formula to conduct correction, the error resulting from the imperfection of the model can not be neglected in many cases.The variability in sign of the topographic correction values is also dealt with in this paper with the aid of the remnant spherical shell model.In view of the wide utilization of the microcomputer, the author suggests that the remnant spherical model should be used in stone slab correction and topographic correction.

No.1 Geophysical and Geochemical Prospecting Party, Bureau of Geology and Mineral Resources of Inner Mongol Autonomous Region

本文给出了二个残球壳的引力公式。将这两个公式用于重力测量中中间层改正和地形改正,并将改正结果与以往常用的无限和有限平板公式的相应计算结果进行了对比。对比结果表明,过去所用的平板改正公式因其模型不完善所引起的误差在许多情况下是不可忽略的。作者还利用残球壳模型论述了地形改正值可以变符号的原因。鉴于当令微型计算机已相当普及,作者建议应该采用残球壳模型进行中间层改正和全球地形改正。

... 概念,认为有必要计算出全球地壳质量的重力校正值^[1,2,3,4] ...

... 雷受旻^[2]、陈凤英^[3]、杨学样^[4]等分别独立推导了该表达式 ...

... 雷受旻^[2]、刘士毅^[3]、杨学样^[4,5]等还阐述了运用该表达式进行#cod#x0201c ...

1990

0.0

杨学样. 关于确定山区区域重力测量中间层校正系数的几点意见[J]. 物探与化探, 1990, 14(5): 396-400.

No. 5 Metallargical Exploration Party of Jiangxi Province

... 概念,认为有必要计算出全球地壳质量的重力校正值^[1,2,3,4] ...

... 雷受旻^[2]、陈凤英^[3]、杨学样^[4]等分别独立推导了该表达式 ...

... 雷受旻^[2]、刘士毅^[3]、杨学样^[4,5]等还阐述了运用该表达式进行#cod#x0201c ...

1992

0.0

... 雷受旻^[2]、刘士毅^[3]、杨学样^[4,5]等还阐述了运用该表达式进行#cod#x0201c ...

1986

0.0

刘立言, 刘定和. 重力地形改正(166. 7 km至全球)改算图[J]. 物探与化探, 1986, 10(5): 381-383.

No. 2 Exploration Geophysics Research Party, MGMR

... 为了减小全球地壳质量重力地形校正值计算量,刘立言等^[6]设计计算了重力地形改正(166 ...

1990

0.0

... 陈邦彦^[7]亲自主持计算了南海中央海盆格莱尼重力异常,探讨了该异常的应用前景 ...

2010

0.0

安玉林, 张明华, 黄金明, 等. 纯球坐标系内各项重力校正值计算方案和计算过程[J]. 物探与化探, 2010, 34(6): 697-705.

The computation project and computation process for gravity correction values within the pure spherical coordinate system are expounded for the first time in this paper. All important formulae in the scheme are educed separately by the authors. The excellences of this computing scheme are as follows: ① By using observational point P as the center of the computation region and by using spherical surface distances am=am-1+Δa·bm-1(m=2,3,…,M) as the radius, the computation region is plotted out as M girdles. Then, By using Δλ0 as the angle interval of radials from P point, all girdles are plotted out as integral N parts. This scheme for plotting out the computation region is extremely precise and its computation precision is very high. ② Under the condition of complex terrain and remarkable different altitudes, the computation precision can be improved by minimizing values of Δa,b and Δλ0. ③ The gravity correction values of very large surveying area will not be confined by the 6°cingulum and the standard map sheet and can directly be computed by using altitude values with longitude and latitude.

1.School of Geophysics and Information Technology，China University of Geosciences, Beijing 100083, China；2. Development and Research Center of China Geological Survey, Beijing 100037, China

在国内，首次阐述了“纯球坐标系内各项重力校正值的计算方案和计算过程”。该计算方案有如下优点：①以观测点P为中心，以球面距离am=am-1+Δa·bm-1(m=2,3,…,M)为半径将计算区划分为M个环带，再以Δλ0为由P点出发射线的角度间隔将所有环带等分为N （整数）份的计算区划分方案是最严密拼接，因而校正值的计算精度相当高。②在地形复杂和相对高差大的情况下，通过减小Δa、b和Δλ0 数值，可进一步提高计算精度。③可以直接采用带有经、纬度坐标的高程数据，进行大面积测区的各项重力校正值计算，而不受6°带和标准图幅限制。

... 以南重力资料校正中^[8] ...

... g(z_i)解析表达式作为地壳质量引起的重力校正值的正演计算公式^[8] ...

... 而且,球面上点的原经纬度坐标向新经纬度坐标转换的公式^[8]严格采用球面三角公式 ...

... 0 km 的大环带区重力校正值^[8]时所采用的计算方案为例(图2) ...

... g(z_i)解析表达式为^[8]: ...

2009

0.0

... 近年来,以陈超为首的科研集体参与国家嫦娥绕月探测工程,开展全月月壳质量重力校正值计算、重力异常圆滑处理和密度横向变化反演等研究,取得了一系列成果^{[9,10,11,12,13,14]} ...

2010

0.0

周聪, 杜劲松, 梁青, 等. 基于球冠域的月球重力地形校正方法[J]. 地球物理学进展, 2010, 25(2): 486-493.

Terrain correction is needed when we do research on the internal structure of the Earth, the Lunar, and Mars by use of satellite gravity measurements, then calculate global Bouguer gravity anomalies. It is effective to carry out terrain correction in spherical coordinates. In this paper, we propose a method of terrain correction in the domain of the spherical cap, and analyze the transformation of different spherical coordinates and the reconstruction of the terrain model that the method involved. Then we conduct a theoretical verification . Based on the CE-1 topographic data and the lunar gravity, the lunar Bouguer gravity anomalies are obtained via our terrain correction method, which are compared with the Tesseroid approach. Furthermore, we analyze the spatial resolution and advantages of our method.

(Institute of Geophysics & Geomatics, China University of Geosciences, Wuhan 430074, China)

利用卫星重力测量资料研究地球、月球与火星等星体的内部构造时，需要进行重力地形校正，计算全球布格重力异常，而在球坐标中实现地形校正计算是一种有效的途径.本文提出球坐标系中的球冠域地形校正计算方法，给出了该方法涉及的球坐标系之间坐标转换方法和球冠域内地形模型重构方法，并进行理论验证．作者利用嫦娥一号激光测高数据对月球重力进行地形校正，获得了月球全球布格重力异常，并与球坐标系中Tesseroid 单元体地形校正方法对比，分析了球冠域地形校正方法的计算精度、空间分辨率及其优缺点．

2012

0.0

杜劲松, 陈超, 梁青, 等. 月球表层及月壳物质密度分布特征[J]. 地球物理学报, 2012, 53(9): 2059-2067.

The variations of rock-density on the lunar surface and inside the crust indicate the component of lunar crust and even helpful in comprehending lunar tectonics and evolution. Based on the distributions of metal elements Fe, Th and Mg on the lunar surface derived from the measurement of Lunar Prospector mission, the global petrologic distributions of plagioclase, Mg-rich rock, basalt and KREEP rock can be interpreted mineralogically and geochemically. In this study, the rock-density distribution on the lunar surface was calculated according to petrologic distributions and densities of the rock samples, meteorites and estimations from recent achievements as well. Values of density in the near surface were adjusted with the relationship between the contents of Fe and rock-density. The density model of lunar crust was determined using separated lunar Bouguer gravity anomaly data, which were calculated from topography model CLTM-s01 and gravity model SGM100h. This model represents the lateral variations of density within the depth of 40 km. The distributions of density in the lunar crust indicate that the lunar crust, at least the upper crust, in most area, is dominated by light silicate rock with rich aluminium, calcium and magnesium elements. It is speculated that the primary lunar crust is composed globally of light and calcium magnesium silicate rocks. The basalt and KREEP rock, as a thin layer, covers the original lunar crust.

月球表层与月壳岩石密度的横向与径向的变化,反映了月表及内部成分以及月球演化等特征.本文利用月球勘探者号伽马射线谱仪探测的月表Fe, Th与Mg元素分布数据,依据前人给出的元素含量与岩石类型的关系,对月球表层进行了岩性填图,并结合岩石样品与陨石的密度测试数据建立初始密度模型,采用铁元素与岩石密度的关系对其进行修正,从而建立了月表物质密度分布模型.基于嫦娥一号激光测高数据和日本SELENE计划发布的月球重力模型,计算出月球布格重力异常,进而反演得到月壳0~40 km深度范围内岩石平均密度分布模型.分析表明,大部分区域上,月壳至少月壳上部岩石成分主要以轻质的富含铝、钙、镁质的硅酸盐类岩石为主.由此推测,原始月壳极有可能是由轻质的、富含钙、镁质硅酸盐类岩石构成的全球性月壳.现今的玄武岩与克里普岩只是覆盖于原始的月壳之上的岩层,且厚度不大.

2012

0.0

... g(z_i)表达式^[12]和双线性内插出需要重构点的坐标的方法^[14] ...

2012

0.0

2012

0.0

... g(z_i)表达式^[12]和双线性内插出需要重构点的坐标的方法^[14] ...

2007

0.0

... 他们计算全月重力异常校正值采用的Tesseroid单元模型^[15]及其#cod#x00394 ...

... 该模型优于其他单体模型^[15],可以形成无缝对接的组合模型,精确模拟地壳质量分布,而且其正演表达式严格精确 ...

0082

0.0

... 而不是预先计算出球面上选定结点处的重力校正值,然后采用双线性内插方法,内插出测点处的重力校正值^[16] ...

... 7 km^[16],因为多数专家认为,166 ...

2005

0.0

... 在柴玉璞与安玉林协作开展的研究中,张家界重力测区重力校正实践^[17]表明,在高程/海深值变化较大的重力测区,166 ...

1979

0.0