利用地形和地壳密度数据进行地表重力仿真的方法及装置与流程

1.本发明属于物理大地测量技术领域，特别涉及一种利用地形和地壳密度数据进行地表重力仿真的方法及装置。


背景技术：

2.现今空白区域的重力数据填补一般分为三种方法，一种是模型填补法，另一种方法邻近数据内插法，还有一种是均衡理论填补法。
3.(1)模型填补法
4.模型填补法主要是利用全球重力场模型计算地表重力场数据。超高阶重力场模型的建立和发展，提高了全球重力场的近似精度。然而，该方法存在以下弊端：一、位系数模型是在球近似条件下计算而来的，而地球是一个不规则球体，这个球近似过程本身存在误差，这会导致计算得到的地表重力数据精度较低；二、构造位系数模型的基础数据主要来源于重力卫星测量的数据和地表实测数据，由于重力卫星飞行高度较高，故其只能反映出地球重力数据的低中频信号，其高频和超高频信号无法通过重力卫星的方法得到，再者，构造位系数的地表重力数据分布极不均匀，这就导致利用位系数模型计算的地表重力数据存在很大的系统误差，整体计算精度较差；
5.(2)邻近数据内插法
[0006]“邻近数据内插法”通过一定搜索半径内的重力观测点，根据shepard模型、kriging模型等数学插值模型得到插值半径内的点。这个方法要求在重力数据的空白区内需要事先测量一定分辨率的重力点，这增加了人工作业量，增加了空白区重力数据填补的工作量，而且若数据分布稀疏，插值半径通常远大于网格分辨率的话，该方法无法实现。
[0007]
(3)均衡理论填补法
[0008]“均衡理论填补法”就是基于均衡理论并利用全球高分辨率的地形数据进行重力数据的仿真，但是这种方式将地壳密度看成常密度，这种密度项上的近似会直接带来仿真数据的系统性误差。


技术实现要素：

[0009]
针对现有技术中重力数据空白区填补的缺陷，本发明提出一种利用地形和地壳密度数据进行地表重力仿真的方法及装置，利用均衡理论并结合高分辨率的地形数据和地球物理信息(地壳分层和密度)进行空白区重力数据填补，以提高仿真精度。
[0010]
为了实现上述目的，本发明采用以下的技术方案：
[0011]
本发明提供了一种利用地形和地壳密度数据进行地表重力仿真的方法，包含以下步骤：
[0012]
确定地表重力数据仿真的范围，结合地形数据和地壳密度数据建立该区域地形和地壳密度数据库；
[0013]
根据所述的数据库，计算出该区域某一分辨率下的地形改正；
[0014]
根据所述的数据库，计算出该区域同一分辨率下的层间改正；
[0015]
根据所述的数据库，计算出该区域同一分辨率下的均衡改正；
[0016]
根据均衡理论，结合地形改正、层间改正和均衡改正，构建区域空间重力异常。
[0017]
进一步地，所述地壳密度数据来自crust1.0密度模型，该模型建立地壳各个密度层高程与密度的对应关系。
[0018]
进一步地，在构建区域空间重力异常之后，还包括：
[0019]
以上述建立的该区域某一分辨率下的空间重力异常数据为基础，再利用数学插值方法计算每个离散点的空间重力异常数据。
[0020]
进一步地，在局部平面坐标系下，积分区σ0采用棱柱积分法，地形改正的计算表达式如下：
[0021][0022]
其中，δg
tc
为计算点p的地形改正，(x,y,z)为流动点平面坐标，g为万有引力常量，h和h
p
分别为流动点和计算点高程，hi为地壳各个密度层的高程，i为地壳密度层层数，n 1为计算点所对应的高程起算面所在的密度层，为第i层(x,y)处的密度值，此公式满足n》1；
[0023]
当n＝1时，
[0024][0025]
当n＝0时，
[0026][0027]
进一步地，所述层间改正的计算表达式如下：
[0028][0029]
其中，δgc为计算点p的层间改正，m 1为计算点所对应的高程起算面所在的密度层，为第i层计算点处的密度值。
[0030]
进一步地，所述均衡改正采用airy—heiskanen模型，该模型认为常密度的地壳漂浮在密度为ρm的地幔上并保持平衡，超出海平面的部分或者低于海平面的部分分别在均衡面以下的部分以“山根”或“反山根”的形式进行补偿；在局部平面直角坐标系下，均衡改正的计算表达式如下：
[0031]
δg
ic
＝g
△
ρ
[0032]
[0033][0034]
其中，δg
ic
为计算点p的均衡改正，z2＝-t，z1＝-(t t)，t为均衡面深度，t为补偿深度，
△
ρ为地壳和地幔部分的密度差。
[0035]
进一步地，根据均衡理论，构建区域空间重力异常，均衡异常应满足浮平衡理论，
△gi
≈0，即均衡重力异常
△gi
应满足如下条件：
[0036]
△g空
＝
△gi-δg
c-δg
tc-δg
ic
[0037]
其中，
△g空
为空间重力异常。
[0038]
本发明还提供了一种利用地形和地壳密度数据进行地表重力仿真的装置，包括：
[0039]
数据库建立模块，用于确定地表重力数据仿真的范围，结合地形数据和地壳密度数据建立该区域地形和地壳密度数据库；
[0040]
地形改正计算模块，用于根据所述的数据库，计算出该区域某一分辨率下的地形改正；
[0041]
层间改正计算模块，用于根据所述的数据库，计算出该区域同一分辨率下的层间改正；
[0042]
均衡改正计算模块，用于根据所述的数据库，计算出该区域同一分辨率下的均衡改正；
[0043]
空间重力异常计算模块，用于根据均衡理论，结合地形改正、层间改正和均衡改正，构建区域空间重力异常。
[0044]
与现有技术相比，本发明具有以下优点：
[0045]
传统的均衡理论填补法中，通常将地壳密度看成一个常量2.67g/m3，将地幔处密度视为3.26g/m3，这种密度项的近似会直接带来仿真数据的系统性误差；基于此，本发明中引入地壳密度项，即将地壳密度看作是变量，这种改变更能贴近地壳实际密度分布。理论上讲，根据这种密度分布计算的地形改正、层间改正和均衡改正的精度更高，更加贴近地球的实际，最后的计算结果也更加趋近于“浮平衡”理论，所以，本发明引入变化的地壳密度可以提高均衡理论填补法的仿真精度。
附图说明
[0046]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0047]
图1是本发明实施例的利用地形和地壳密度数据进行地表重力仿真的方法的流程图；
[0048]
图2是本发明实施例的实验区域的地形图；
[0049]
图3是本发明实施例的crust1.0全球地壳模型示意图；
[0050]
图4是本发明实施例的实验区域地形改正分布图；
[0051]
图5是本发明实施例的实验区域层间改正分布图；
[0052]
图6是本发明实施例的实验区域均衡改正分布图；
[0053]
图7是本发明实施例的实验区域空间重力异常分布图；
[0054]
图8是本发明实施例的用于评估本发明方法精度的点位分布图；
[0055]
图9是本发明实施例的crust1.0密度模型分布示意图。
具体实施方式
[0056]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0057]
目前，全球地形模型已达到相当分辨率，例如tandem-x、srtm系列，astergdem和aw3d30，其中srtm的分辨率可以达到12m，研究表明重力异常与地形在中高频部分表现出的强线性关系，这为利用均衡理论和地形数据构建局部区域高分辨率重力异常提供可能。而且，全球地壳模型crust1.0包含丰富的地球内部信息，分辨率已经达到1
°
，这为利用均衡模型和地形数据构造空白区重力数据过程中加入地球物理信息(地壳分层和密度)提供了前提。
[0058]
如图1所示，本实施例的一种利用地形和地壳密度数据进行地表重力仿真的方法，包含以下步骤：
[0059]
步骤s101，确定地表重力数据仿真的范围，结合高精度的地形数据和地壳密度数据建立该区域地形和地壳密度数据库。
[0060]
优选的，地壳密度数据来自crust1.0密度模型(如图9所示)，图中hi(i＝1
…
9)和ρi(i＝1
…
9)分别表示地壳各个密度层的高程和密度，其平面分辨率为1
°
。
[0061]
步骤s102，根据步骤s101中的数据库，计算出该区域某一分辨率下的地形改正，在局部平面坐标系下，局部地形改正计算时，积分区σ0采用棱柱积分法，地形改正的计算表达式如下：
[0062][0063]
其中，δg
tc
为计算点p的地形改正，(x,y,z)为流动点平面坐标，g为万有引力常量，h和h
p
分别为流动点和计算点高程，hi为地壳各个密度层的高程，i为地壳密度层层数，n 1为计算点所对应的高程起算面所在的密度层，为第i层(x,y)处的密度值，此公式满足n》1。
[0064]
当n＝1时，
[0065][0066]
当n＝0时，
[0067]
[0068]
步骤s103，根据步骤s101中的数据库，计算出该区域同一分辨率下的层间改正，层间改正的计算表达式如下：
[0069][0070]
其中，δgc为计算点p的层间改正，m 1为计算点所对应的高程起算面所在的密度层，ρ
ip
为第i层计算点处的密度值。
[0071]
步骤s104，根据步骤s101中的数据库，计算出该区域同一分辨率下的均衡改正，均衡改正采用较为符合地球实际均衡状态的airy—heiskanen模型，该理论认为常密度的地壳漂浮在密度为ρm的地幔上并保持平衡，超出海平面的部分或者低于海平面的部分分别在均衡面以下的部分以“山根”或“反山根”的形式进行补偿；在局部平面直角坐标系下，均衡改正的计算表达式如下：
[0072]
δg
ic
＝g
△
ρ
[0073][0074][0075]
其中，δg
ic
为计算点p的均衡改正，z2＝-t，z1＝-(t t)，t为均衡面深度，t为补偿深度，
△
ρ为地壳和地幔部分的密度差。
[0076]
步骤s105，根据均衡理论，结合步骤s102-s104计算的地形改正、层间改正和均衡改正，构建区域空间重力异常，均衡异常应满足浮平衡理论，
△gi
≈0，即均衡重力异常
△gi
应满足如下条件：
[0077]
△g空
＝
△gi-δg
c-δg
tc-δg
ic
[0078]
其中，
△g空
为空间重力异常。
[0079]
步骤s106，以步骤s105中建立的该区域某一分辨率下的空间重力异常数据为基础，再利用数学插值方法计算每个离散点的空间重力异常数据。
[0080]
下面举一个具体实例，以便更好地理解本发明。
[0081]
步骤s201，确定地表重力数据仿真的范围，结合地形数据和地壳密度数据建立该区域地形和地壳密度数据库。
[0082]
(1)地形数据。
[0083]
实验区分布在非洲东南部津巴布韦境内的马塔贝莱高原。实验区选为6
°×6°
区域，地域分布范围为东经20
°‑
26
°
，南纬28
°‑
34
°
，本实例采用etopo1所发布1
′×1′
的地形数据，etopo1是由ngdc(美国国家地球物理数据中心)开发的地形高程数据，高程信息统计见表1，实验区的海拔高均值在1000m左右，均方差为320m，实验区地形变化剧烈，重力场信号复杂，实验区域的地形图如图2所示。
[0084]
表1高程信息统计表
[0085][0086]
(2)地壳密度数据。
[0087]
crust1.0是分辨率为1
°×1°
的全球三维地壳模型，crust1.0模型将全球划分为1
°×1°
的网络，共计64800个单元。模型将地壳分为9层，第一层到第九层分别是：水层、冰层、上沉积层、中沉积层、下沉积层、上地壳、中地壳、下地壳、地幔。模型中水深数据以及地形数据来自于etopo1，地幔的参数来自llnl-g3dv3模型。实验中需要从全球地壳模型中提取出实验区域的地壳深度以及密度信息。在6
°×6°
的实验区域内，共包含了36
×
9个密度项，其中地壳平均密度最大值为2390kgm-3
，最小值为2080kgm-3
。实际的平均密度值与常密度存在一定差异。crust1.0全球地壳模型如图3所示。
[0088]
步骤s202，根据步骤s201中的数据，利用步骤s102中地形改正公式计算出地形改正值如表2，该区域地形改正分布图如图4所示。
[0089]
表2地形改正值统计表
[0090][0091]
步骤s203，根据步骤s201中的数据，利用步骤s103中层间改正公式计算出层间改正值如表3，该区域层间改正分布图如图5所示。
[0092]
表3层间改正值统计表
[0093][0094]
步骤s204，根据步骤s201中的数据，利用步骤s104中均衡改正公式计算出均衡改正值如表4，该区域均衡改正分布图如图6所示。
[0095]
表4均衡改正值统计表
[0096][0097]
步骤s205，利用步骤s202至步骤s204计算出的数据，结合步骤s105中空间重力异常公式计算出该区域空间重力异常，该区域空间重力异常分布图如图7所示。
[0098]
步骤s206，以步骤s205中建立的该区域某一分辨率下的空间重力异常数据为基础，再利用数学插值方法计算每个离散点的空间重力异常数据。
[0099]
对步骤s205中的数据进行精度评定，用于校核的重力实测数据是从“gravcd-africa”数据集提取，包含1789个该地区实测空间重力异常数据，重力异常数据统计如表5，点位分布图如图8所示。
[0100]
表5重力异常数值统计表
[0101][0102]
利用kriging插值方法，结合本发明计算的该区域内一定分辨率下格网化空间重力异常值，得到图8中每个离散点的重力异常值，并与每个点的真实值进行对比，得到对比结果如表6和表7：
[0103]
表6传统均衡理论填补法空白区域填补差值统计表
[0104][0105]
表7本发明空白区域填补差值统计表
[0106][0107]
对比实验结果表明：变密度空白区域重力数据填补差值的平均值为-35.3848mgal，常密度空白区域重力数据填补差值的平均值为-51.3602mgal，二者均方差相当。这说明加入密度项后，相比于常密度项系统差更小，更加有利于进行空白区的重力数据填补。可见利用crust1.0模型进行空白区域的重力数据填补的系统差更小，且平均值明显优于常密度进行空白区域填补，利用本发明进行重力数据填补的效果更好，更加符合实际要求。
[0108]
与上述利用地形和地壳密度数据进行地表重力仿真的方法相应地，本实施例还提供利用地形和地壳密度数据进行地表重力仿真的装置，包括：
[0109]
数据库建立模块，用于确定地表重力数据仿真的范围，结合地形数据和地壳密度数据建立该区域地形和地壳密度数据库；
[0110]
地形改正计算模块，用于根据所述的数据库，计算出该区域某一分辨率下的地形改正；
[0111]
层间改正计算模块，用于根据所述的数据库，计算出该区域同一分辨率下的层间改正；
[0112]
均衡改正计算模块，用于根据所述的数据库，计算出该区域同一分辨率下的均衡改正；
[0113]
空间重力异常计算模块，用于根据均衡理论，结合地形改正、层间改正和均衡改正，构建区域空间重力异常。
[0114]
本发明根据重力数据高频部分与地形数据的强相关性，利用均衡理论和已发布的全球高分辨率的地形数据、地壳密度数据进行地表重力数据仿真，即仿真出区域内空间重力异常(即重力数据)，以减少地表人工测量的成本。
[0115]
需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。
[0116]
最后需要说明的是：以上所述仅为本发明的较佳实施例，仅用于说明本发明的技术方案，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所做的任何修改、等同替换、改进等，均包含在本发明的保护范围内。