利用地质信息变密度正演计算布格重力异常的方法及装置与流程

1.本发明属于重力勘探技术领域，尤其涉及一种利用地质信息变密度正演计算布格重力异常的方法。


背景技术：

2.目前，重力勘探在成矿预测、圈定岩体和构造、研究金属矿赋存特点等金属矿普查方面都是卓有成效的，有着良好的应用前景。重力勘探的观测对象为测点处的重力加速度值，一般采用相对测量的方式，测量测点与重力基点之间的段差，利用重力基点的重力加速度值和基点与测点间的段差，获得各测点的重力加速度值，简称重力值。
3.传统的布格重力异常值由重力值经正常场改正(纬度改正)、中间层改正、高度改正(中间层改正和高度改正合称布格改正)和地形改正后获得，是重力勘探的基础数据。传统布格重力异常的计算步骤及物理意义如图1。
4.如图1，地面上任一点的实测重力值取决于纬度、高度、周围地形变化和地下局部物质的不均匀分布等，而最后一个因素是重力勘探的重要研究对象，其他因素所产生的重力变化都是干扰，因此，测点处布格重力异常的计算就是将干扰因素剥离的过程。假设地球是一个密度成层分布的光滑椭球体，在同一层内密度是均匀的、各层界面也都是共焦旋转椭球面，则椭球面上各点的重力值，可根据地球的万有引力常数、长半径、扁率和自转角速度等计算得出，该光滑椭球体称为正常椭球，正常椭球面为大地水准面，正常椭球面上的重力值称为正常重力值。
5.图1中的a：p、q测点的实测重力值由大地水准面以下分层的正常重力值、大地水准面以上至地形起伏面间物质引起的重力值和测点高度不同而引起的重力值叠加而成。其中正常椭球体以正常地壳厚度平面划分两个密度界面，大地水准面至正常地壳厚度平面间理论密度(黑色字)为2.67g/cm3，正常地壳厚度平面以下理论密度为3.27g/cm3。上述理论模型中存在两大类物质不均匀分布状况，一是实际地壳下界面为莫霍面，莫霍面是起伏不平的，一般不与正常地壳厚度平面重合，莫霍面以上地壳的理论密度(白色字)为2.67g/cm3，因此正常地壳厚度平面与莫霍面间存在物质不均匀分布的情况；二是不同深度分布的各类局部实际密度低于或高于地壳平均密度的地质体，即我们重力勘探重点研究的对象。
6.图1中的b：正常重力值改正即消除大地水准面以下各层物质引起的重力值变化，消除后各地面测点的重力值由莫霍面与正常的壳厚度平面间物质的残余重力值(残余密度差为－0.6g/cm3)、大地水准面以上至地形起伏面间物质引起的重力值及各局部密度不均匀地质体引起的重力值组成。cgcs2000椭球正常重力值计算公式如下：
[0007][0008]
式中：g0为正常重力值，为测点纬度。
[0009]
图1中的c：地形改正即消除以测点高程面为基准面，测点周围一定范围内，地形起伏引起的重力值变化。如p点周围的1、2、3区和q点周围的4、5区，其中1、3、4区为物质亏空，2、5区为物质盈余。地形改正后，各地面测点的重力值由正常场改正残余重力值(图1中的
b)、大地水准面至测点高程平面间规则层状物质引起的重力值及各局部密度不均匀地质体引起的重力值组成。地形影响计算示意图如图2。
[0010]
如图2，取直角坐标系xyz，将测点a所在的位置定为原点，z轴垂直向下，x、y轴在a点所在的水平面内，dm为质量元，其所在位置的坐标为(x,y,z),a点至质量元的矢径用r表示，r与z轴的夹角为θ，dm在a点产生的重力影响，即引力的垂直分量dg为：
[0011][0012]
若计算区域(1)或(2)全部质量的影响值δgt，可用积分计算：
[0013][0014]
式中：g为万有引力常数，v为地质体积，(x，y，z)为地质体质心坐标，ρ为地质体剩余密度。
[0015]
质量盈余区域(如区域(1))中的z为负值，ρ为正值，质量亏损区域(如区域(2))中的z为正值，ρ为负值，因此，地形对测点的重力影响值永远为负值，即地形起伏使重力值减少。
[0016]
图1中的d：中间层改正和高度改正。中间层改正即消除大地水准面至测点高程平面间的规则层状物质引起的重力值，高度改正即消除因测点高程不同而正常重力场随高度变化值的不同。消除后的重力值称为布格重力异常值，包含莫霍面起伏变化引起的区域重力异常和各局部密度不均匀地质体引起的局部重力异常。
[0017]
中间层改正值δgm的计算公式如下：
[0018][0019]
式中：ρ为密度，h为中间层厚度，即图d中测点高程平面至大地水准面的距离，a为中间层圆盘改正半径，一般等于地形改正半径。
[0020]
高度改正δgh的计算公式为：
[0021][0022]
式中：h为测点高程平面至大地水准面的距离，为测点纬度。
[0023]
布格重力异常的计算：
[0024]
δgb＝g
观-g0+δgh+δgm+δgt
ꢀꢀꢀ
(6)
[0025]
式中：δgb为布格重力异常值，g
观
为测点实测重力值，g0为正常重力值，δgh为高度改正值，δgm为中间层改正值，δgt为地形改正值。
[0026]
上述传统的布格重力异常计算方法，适用于中小比例尺区域性重力勘探的计算处理，若用于大比例尺重力勘探或矿区尺度的重力勘探时，则在中间层改正和地形改正中存在四个明显的缺陷：
[0027]
1、传统中间层改正、地形改正的密度均为一统一值2.67g/cm3，平面上明显不符合地质体及密度变化的事实，据此计算的布格重力异常，其中所含的局部异常信息，易被密度变化引起的重力异常所淹没，无法有效反映目标地质体的实际重力异常。
[0028]
2、中间层改正垂向上没有分层计算，则没有充分利用被勘探区域的地质信息，垂
向上随着不同密度地层产状的变化，深部地质体密度已然发生变化，同时，矿区大量的钻孔岩芯或开采矿石等，均能提供大量关键的密度信息，若中间层以一层笼统计算，则不利于有效区分已掌握地质体和新发现地质体，往往不能取得较好的重力勘探效果。
[0029]
3、地形改正的亏空部分，其亏空物质的密度未知，仅以2.67g/cm3这一理论密度值进行计算，可能带来地形改正的误差，不利于布格重力异常的精细求解，进而影响精细划分区域重力异常和局部重力异常。
[0030]
4、地形改正若采用方域计算公式，而中间层理论计算公式为一圆盘状，则地形改正与中间层改正，存在改正范围不一致的情况，带来布格重力异常的精度损失。
[0031]
通过上述分析，现有技术存在的问题及缺陷为：现有的重力异常计算结果精度较低，误差大，有时无法有效反映目标地质体引起的实际重力异常，不能为大比例尺重力勘探提供精确信息，重力勘探效果不佳。
[0032]
现阶段，三维建模的部分软件，如geomodel等，通过地质建模，基本利用了各类地质信息，能够进行地质体重力正演计算，但仍存在两大问题：一是正演计算还是利用2.67g/cm3这一统一密度，未能进行变密度计算，需要改进；二是软件操作复杂，费用贵，难以普及和推广。


技术实现要素：

[0033]
针对现有技术存在的问题，本发明提供了一种利用地质信息变密度正演计算布格重力异常的方法。
[0034]
本发明是这样实现的，一种利用地质信息变密度正演计算布格重力异常的方法，所述利用地质信息变密度正演计算布格重力异常的方法包括：
[0035]
利用被勘探区域的地质信息，包括地质图中地质体的界线、产状、深度或钻孔内垂向上地质体的密度变化，按照地质体的三维分布状态，不仅在平面上分块，还需在垂向上分层，以不同地质体赋予其实测的密度值，来变密度计算；
[0036]
利用高分辨dem数据统一正演计算中间层物质影响值代替地形改正和中间层改正的分步计算，提高布格重力异常的精度；
[0037]
利用三维地质体变密度计算，采用了深部地层产状为野外地质实测或钻孔控制的任意角度。
[0038]
进一步，所述方法具体包括：
[0039]
实测重力值经正常场改正和高度改正后，得到gfi改正重力值；根据地质图、钻孔等各种平面、深部的密度信息，以dem网格数据将地质体按平面、垂向(垂向分层至大地水准面，即海拔0水平面)分层、分块后，利用直立长方体公式正演计算分层、分块的各地质体对某一重力测点的影响值，将所有地质体的影响值求和得到大地水准面至地形起伏面间物质对该测点重力值的中间物质改正值；将上述gfi改正重力值减去该中间物质改正值，即为精细计算的布格重力异常值。
[0040]
进一步，所述利用被勘探区域的地质信息变密度正演计算布格重力异常的方法包括以下步骤：
[0041]
步骤一，测点重力值进行正常场、高度改正；计算dem网格间距，并进行网格数据的垂向分层，垂向分层后以网格点为中心，平面上向四周延伸半网格间距，构成直立长方体顶
面，由垂向分层和顶面构造直立长方体，计算直立长方体中心点高程、所述直立长方体的半高；
[0042]
步骤二，利用实测密度值准备三维地质体(一个地质体具有唯一密度值，因此也称密度体)界面文件；以得到的各三维密度体分割网格分层数据，分割后得到属于各地质体的用于正演计算的数据文件，包括网格间距、直立长方体的平面坐标、直立长方体中心点高程、直立长方体半高、直立长方体的实测密度值；
[0043]
步骤三，直立长方体公式变密度正演计算获得中间层物质影响值；进行精细布格重力异常计算。
[0044]
进一步，步骤一所述测点重力值进行正常场、高度改正包括：
[0045]
测点重力数据包括测点x、y、h三维坐标和g
观
实测重力值四列，进行正常场改正和高度改正得到测点x、y、h三维坐标和gfi改正重力值。
[0046]
进一步，步骤一所述进行网格数据的分层包括：
[0047]
将网格文件按一定的垂向间隔分层，根据正演计算固定网格间距下重力值随深度的变化率，将全区dem高程按照全区高程最低点以下400m分为上下两个大层：
[0048]
所述上一大层以20m间隔连续分层至全区最高点，该大层对实际地表网格高程低于所述划分层高程的网格点进行舍弃；
[0049]
所述下一大层再按最低高程400-500m以50m间隔分两层，500-700m以100m间隔分为两层，700-1100m以200m间隔分为两层，1100-1500m以400m间隔分为一层，1500-2300m以800m间隔分为一层，若被勘探区dem网格最低点高程大于2700m，则2300m至海拔0m标高按实际剩余间隔分为一层，若最低点高程小于2700m，则1500m至海拔0m标高按实际剩余间隔分为一层，该大层中所有的网格点均保留；
[0050]
进一步，因为dem数据为规则网格，即网格间距为唯一值，因此步骤一所述的网格间距等于某两个相邻网格点之间距离；
[0051]
进一步，步骤一所述的分层后网格中心点高程和直立长方体半高，利用分层的高程值进行计算，各中心点高程为上下分层高程值(若dem网格点的实际高程值小于上分层高程值，则取dem网格点的实际高程值为上分层高程值，下同)的和除以2，各网格立方体的半高为上下分层高程值的差除以2。
[0052]
进一步，步骤二所述利用实测密度值准备地质界面文件包括：
[0053]
(1)依照重力勘探区的地质图或其他钻孔分岩性采集测量统计实际密度值，得到每个地质体的统计密度值；所述每个地质体的统计密度值包括：地质体名称和统计密度值；
[0054]
(2)准备地质界面文件：
[0055]
简化地质体的边界，统计区内地质体的各界面信息，以多个界面圈闭地质体，以圈闭的地质体分割步骤一中完成的网格数据分层文件，分割的文件匹配该圈闭地质体的实测密度值，即完成密度体准备。
[0056]
所述各界面信息包括：以地质界线位置坐标、产状或钻孔位置、深度进行数字化，形成某地质体的多个界面信息。所述界面信息包括：地质界线边界点1和点2的x1、y1、z1、x2、y2、z2、倾角a、体与面的关系j、走向与倾向的加减关系。
[0057]
准备完成后用于正演的数据包括：网格点的x、y坐标，直立长方体的中心点高程，直立长方体的半高，实测密度值。
[0058]
进一步，步骤三所述利用直立长方体正演计算获得中间层物质影响值包括：
[0059]
基于得到的用于直立长方体正演计算的地质体数据文件，结合测点重力值文件中的x、y、h，利用直立长方体正演公式进行各地质体数据文件中单个直立长方体对任一测点的正演计算，并将所有直立长方体的正演值相加，得到所述地质体对某一重力测点中间层物质影响值；得到测点的x、y、h和g
中
。
[0060]
进一步，步骤三所述进行精细布格重力异常计算包括：
[0061]
将计算得到的改正重力值gfi与得到的g
中
相减，得到重力测点的精细布格重力异常值。
[0062]
本发明的另一目的在于提供一种利用地质信息变密度正演计算布格重力异常的装置，该装置包括：获取模块，获取被勘探区域的地质信息(地质图及相关地质体的实测密度值)、重力测点的三维坐标、重力测点的实测重力值和dem网格数据；第一确定模块，用于根据重力测点的三维坐标和实测重力值，进行正常场改正和高度改正，获得gfi改正重力值；第二确定模块，用于将dem网格数据，根据一定直立长方体随与测点距离(即深度)和大小(网格间距一定，即长方体高度、分层间隔)的变化率原则垂向分层，形成网格分层文件；第三确定模块，用于将简化整理的地质信息，整理为地质界面文件，根据界面文件，圈闭分割上述网格分层文件，匹配地质体实测密度，最终形成直立长方体正演数据文件；第四确定模块，用于正演计算中间层物质影响值；第五确定模块，用于计算精细的布格重力异常值。
[0063]
本发明的另一目的在于提供一种计算机设备，所述计算机设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行所述利用地质信息变密度正演计算布格重力异常的方法的步骤。
[0064]
结合上述的技术方案和解决的技术问题，从以下几方面分析本发明所要保护的技术方案所具备的优点及积极效果为：
[0065]
本发明能够充分利用被勘探区域的地质信息，包括但不限于地质图中地质体的界线、产状、深度或钻孔内垂向上地质体的密度变化等，即一切已经掌握的地质信息；传统处理时的地形改正和中间层改正中，平面上默认所有地质体为一密度唯一体，垂向上默认地质体的产状为直立或水平状这两类简单特殊的模型，采集统计的物性在异常解释时仅起到定性相对高低的作用，而实际地质情况是，平面上分区的地质体密度差异大，不能以统一的密度来计算和模拟，垂向上地质体的产状为任意的角度，非直立或水平这两类特殊的情况所能代替，因此，在已经采集统计了相关地质体密度的情况下，应按照地质体的三维分布状态，不仅在平面上分块，还需在垂向上分层，以不同地质体赋予其实测的密度值，来变密度计算，才能区分和提取有效异常，提高勘探精度和效果。
[0066]
本发明采用高分辨dem数据统一正演计算中间层物质影响值，代替了地形改正和中间层改正的分步计算，提高了布格重力异常的精度；
[0067]
本发明采用三维地质体变密度计算，即与传统变密度地形改正不同，传统的变密度地形改正计算仅考虑平面地质体密度的变化情况，垂向上即认为是直立地层(陡倾90
°
)；而本发明同时综合了深部地质体的变化情况，即深部地层产状为野外地质实测或钻孔控制的任意角度，这更加符合地质实际。
[0068]
本发明能够改进已有软件不能充分利用地质信息(产状、不同地质体的密度变化)、利用统一密度计算的缺陷，精细计算的布格重力异常更加符合实际的地质情况。
[0069]
本发明的技术方案转化后的预期收益和商业价值为：本发明预期可以转换为精细计算布格重力异常的软件，现阶段矿产资源形势严峻，矿区尺度的大比例尺重力勘探工作势在必行，则本发明的方法是提升勘探精度，改进勘探效果，发现深边部矿产资源的必要手段；将该发明软件化后，能够进行出售或给地质勘探单位进行技术服务，商业价值明显。
[0070]
本发明的技术方案填补了国内外业内技术空白：本发明填补了充分利用地质信息变密度正演进行精细布格重力异常的方法空白。
附图说明
[0071]
图1是本发明实施例提供的传统布格重力异常的计算过程及物理意义示意图；图中：p、q为地表重力测点；
[0072]
图2是本发明实施例提供的地形影响示意图；
[0073]
图3是本发明实施例提供的利用地质信息变密度正演计算布格重力异常的方法流程图；
[0074]
图4是本发明实施例提供的重力测点与地形网格节点之间的相互关系示意图；
[0075]
图5是本发明实施例提供的直立长方体模型计算示意图；
[0076]
图6是本发明实施例提供的利用直立长方体变密度正演计算布格重力异常的方法流程图；
[0077]
图7是本发明实施例提供的重力剖面及测点周围地形图；
[0078]
图8是本发明实施例提供的重力剖面及测点周围地质简化图；
[0079]
图9是本发明实施例提供的剖面上各地质体产状变化示意图；
[0080]
图10是本发明实施例提供的高程1480-1500m网格数据的分布情况示意图；
[0081]
图11是本发明实施例提供的σⅲ地质体的直立长方体示意图；
[0082]
图12是本发明实施例提供的传统方法与改进方法的勘探效果对比示意图。
具体实施方式
[0083]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0084]
本发明实施例提供了一种利用地质信息变密度正演计算布格重力异常的方法，该方法包括：获取被勘探区域的地质信息(地质图及相关地质体的实测密度值)、重力测点的三维坐标、重力测点的实测重力值和dem网格数据(x、y、z1)；利用重力测点的三维坐标和实测重力值，经公式(1)和(5)完成正常场和高度改正，获得改正后的重力值gfi；计算dem网格数据的网格间距，将dem网格数据按照上述技术方案分层，分层后获得直立长方体数据文件，同时计算出直立长方体的中心点高程zm和直立长方体半高c；利用地质信息，建立地质界面文件，利用建立的多个地质界面，分割上述直立长方体数据文件，构成某个地质体直立长方体文件，并匹配该地质体的实测密度，循环完成所有地质体的直立长方体文件分割和实测密度匹配，即完成了正演数据准备；利用重力测点的三维坐标、网格间距和正演数据，根据直立长方体正演公式，进行中间层物质影响值计算，获得中间层物质改正值g
中
；将计算得到的改正重力值gfi与得到的g
中
相减，得到重力测点的精细布格重力异常值。
[0085]
本技术第一步中重力测点的三维坐标和实测重力值经正常场改正和高度改正获得gfi重力改正值的方法，为成熟的方案，即已有公式可供直接改正，因此不再详述。
[0086]
本技术第二步的网格数据分层方案为本发明的第一个创新点，首先，需要根据dem数据的网格间距(网格间距即为长方体的长和宽)，密度采用理论密度2.67g/cm3，以一定深度间隔确定该网格间距对应的影响值随深度的变化率，比如上述的以最低高程网格点下减400m分上下两个大层，变化长方体的中心点高度z和半高c，根据公式(3)，计算不同高度和半高的长方体对正上方测点影响值，以该影响值小于等于0.001
×
10-5
m/s2为准，确定分层中使用的中心点高度和长方体半高，如上述采用的最低高程网格点以下400m分为上下两个大层，上一大层以20m间隔分层，即深度400m以下(最低高程网格点对应的长方体中心点高度为该网格点高程减去400/2)对应长宽高为10m
×
10m
×
20m的直立长方体对测点的影响值小于等于0.001
×
10-5
m/s2；即该分层所遵循的原则即为分层深度和分层间隔根据网格间距进行影响值变化率计算试验确定，网格间距大则对应分层深度更小，分层间隔更细，并不仅限于本技术上述方案中采用的分层数据；分层时，上一大层中网格点高程低于最低分层高程的网格点将被舍弃，而下一大层中所有的网格点均被保留。完成后，形成未分割的唯一的直立长方体数据文件，包括网格点坐标x和y，长方体中点高程zm和长方体的半高c四列数据，其中网格点坐标x和y因每一层均对应这些网格点，则重复出现多次。
[0087]
本技术第三步的地质界面文件准备、根据地质界面分割圈闭地质体所属的直立长方体文件并赋予密度值，形成正演数据文件，为本发明的第二大创新点。将地质图中的地质界线简化为系列直线，提取直线的端点坐标x1，y1，z1，x2，y2，z2，以及该地质界面的倾角a，地质体与面的关系j(如面上为1，面下为－1)及走向与倾向的加减关系(如逆时针为1，顺时针为－1)，由多个地质界面(即为多行地质界面信息)数据圈闭分割上述第二步完成的直立长方体文件，形成组成一个三维地质体的所有直立长方体文件，并给每个直立长方体文件赋予该地质体的实测统计密度；循环完成所有地质体的直立长方体文件，形成正演数据文件。
[0088]
本技术第四步中根据直立长方体公式变密度正演中间层影响值，其计算公式为已有公式，而利用正演计算代替传统布格重力异常计算中的中间层改正和地形改正的思路，为本发明的第三大创新点，其意义在于一是解决了传统中间层和地形改正中不能变密度计算，不能充分利用地质信息的缺陷；二是地形改正的“亏空”部分地质体的密度值实际不存在，为一理想状态，必然带来精度损失，而本发明的正演计算全部利用实测密度值进行，更加符合地质事实；三是传统的中间层计算公式为圆盘模型，而地形改正为方域模型，两者之间有圆域和方域的衔接精度损失，而本发明的正演计算将两者统一进行，不存在中间层和地形改正的衔接问题，精度得到明显改善。
[0089]
本技术第五步为简单的加减法，即将正常场和高度改正后的gfi减去中间层物质影响值gfi，获得精细的布格重力异常值。精细计算布格重力异常，便于有效区分重力有用异常和干扰异常，即区分目标体和干扰体引起的异常，提升重力勘探的效果，填补了大比例尺重力勘探精细解释的空白。
[0090]
第二、本发明实施例提供一种利用地质信息变密度正演计算布格重力异常的装置，该装置包括：获取模块，获取被勘探区域的地质信息(地质图及相关地质体的实测密度值)、重力测点的三维坐标、重力测点的实测重力值和dem网格数据；第一确定模块，用于根
据重力测点的三维坐标和实测重力值，进行正常场改正和高度改正，获得gfi改正重力值；第二确定模块，用于将dem网格数据，根据一定直立长方体随与测点距离(即深度)和大小(网格间距一定，即长方体高度、分层间隔)的变化率原则垂向分层，形成网格分层文件；第三确定模块，用于将简化整理的地质信息，整理为地质界面文件，根据界面文件，圈闭分割上述网格分层文件，匹配地质体实测密度，最终形成直立长方体正演数据文件；第四确定模块，用于正演计算中间层物质影响值；第五确定模块，用于计算精细的布格重力异常值。
[0091]
本发明能够改进已有软件不能充分利用地质信息(产状、不同地质体的密度变化)统一密度分地形改正和中间层改正计算传统布格重力异常的缺陷，精细计算的布格重力异常更加符合实际的地质情况。
[0092]
如图3所示，本发明实施例提供的利用地质信息变密度正演计算布格重力异常的方法包括以下步骤：
[0093]
s101，测点重力值进行正常场、高度改正；计算dem网格数据的间距，试验确定分层深度和分层间隔，进行dem网格数据的分层，即获得直立长方体数据，并计算得到直立长方体的中心点高程，直立长方体的半高；
[0094]
s102，利用地质信息准备地质界面文件；以地质界面圈闭分割网格分层后的直立长方体文件，获得某一地质体的直立长方体文件，匹配该地质体的密度，循环完成所有地质体的直立长方体数据圈闭和分割，匹配对应地质体密度后得到属于所述地质体的正演数据；
[0095]
s103，直立长方体正演计算获得中间层物质影响值；进行精细布格重力异常计算。
[0096]
本发明实施例提供的利用地质信息变密度正演计算布格重力异常的方法包括：
[0097]
将传统计算方法中的地形改正和中间层改正改进为利用地质信息变密度正演计算的方式，即根据地质图、钻孔等各种平面、深部的密度信息，平面、垂向分层、分块进行直立长方体正演计算其对重力测点的影响后，求和即为大地水准面至地形起伏面间物质对测点重力值的综合影响。改进后的布格重力异常计算方法为正常场改正、高度改正和中间物质正演改正。
[0098]
为了便于理解直立长方体公式，下面对直立长方体公式进行描述。
[0099]
实际地形面的起伏变化，可以用不同的面来拟合，通常采用的方法有正方形表面积分法、正方形中点高程法和正方形平均高程法，地形面的拟合方法确定后，向垂向方向延伸一定深度，即组成一个个三维模型体。相关研究表明，网格间距小于10m
×
10m时，上述三种拟合方法组成的三维体，对同一测点的重力影响值基本一致，考虑到地形网格数据的表现形式为节点三维坐标这一特点(图4)，且大比例尺重力勘探中，一般已测绘了大比例尺、高精度地形数据，其网格间距远小于10m。直立长方体法在计算数据量巨大时的运算速度上、便于计算机编程实现上具有明显的优点，因此采用直立长方体公式进行中间物质正演改正计算。
[0100]
如图5，o为坐标原点，z轴铅垂向下(向下正方向)，x、y相互垂直构成水平面，设直立长方体的两对侧面分别平行于xoz和yoz面，顶底面平行于xoy面，剩余密度为σ，地质体内某一体积元dv＝dxdydz，其在xyz坐标系中的坐标为(x，y，z)，它的剩余质量为dm，则dm＝σdv＝σdxdydz，剩余质量元到计算点p(x
p
，y
p
，z
p
)的距离为r＝[(x-x
p
)2+(y-y
p
)2+(z-z
p
)2]
1/2
。所以，地质体的剩余质量对某点产生的引力位w为：
[0101][0102]
式中：g为万有引力常数，v为地质体的体积。
[0103]
因为z方向就是重力的方向，所以重力异常就是剩余质量的引力位沿z方向的导数：
[0104][0105]
将坐标原点平移至计算点p，则被积函数的测点坐标x
p
、y
p
、z
p
均为0，解出该积分：
[0106][0107]
将坐标原点还原，以c点表示长方体的中心点，坐标为(xc，yc，zc)，长、宽、高分别为2a、2b和2c，则式中：
[0108]
r＝[(x
c-x
p
)2+(y
c-y
p
)2+(z
c-z
p
)2]
1/2
[0109]
x1＝x
c-a-x
p
，y1＝y
c-b-y
p
，z1＝z
c-c-z
p
[0110]
x2＝xc+a-x
p
，y2＝yc+a-y
p
，z2＝zc+c-z
p
[0111]
以上积分形式可以表示为：
[0112][0113]
式中：x∈[x1，x2]，y∈[y1，y2]，z∈[z1，z2]
[0114]
为了便于编程计算，减少循环以提高运算速度，可以将(10)式进一步分解为8个计算公式，8个公式的计算结果即对应图5中8个角点的地形改正值，求和后即为直立长方体的正演值。
[0115]
采用dem数据，利用上述公式进行正演计算时，还需注意以下几个方面：
[0116]
①
由公式可知，上下两半空间内相同高度直立长方体对重力测点的影响值大小相等，符号相反，即高于重力测点平面的正演值为负，低于测点平面的正演值为正，其不同于地形改正值，地形改正值仅为正值；
[0117]
②
图5中的a和b对应dem网格点的半网格间距(图4)。则长方体半高c可以利用顶面高程和底面高程计算得出；
[0118]
③
dem网格节点的三维坐标为直立长方体顶面的中点坐标(图4中c
′
点)，而公式中为直立长方体中心的中点坐标，两者的z值需要利用直立长方体的半高c进行转换；
[0119]
④
(10)式中的ln(y+r)、ln(x+r)和当x＝y＝z＝0、x＝0&z＝0或y＝0&z＝0时，函数存在奇点，实际编程中可以将z坐标加一微小值予以避免。
[0120]
为了便于理解本技术实施例，将实现过程详述如下：
[0121]
从实测绝对重力值、dem网格数据和实测密度值(平面或垂向，不限于地质图)，需要通过以下过程实现精细布格重力异常的计算，其技术路线如图6。
[0122]
1、测点重力值进行正常场、高度改正
[0123]
测点重力数据包括测点x、y、h三维坐标和g
观
实测重力值四列，正常场改正和高度改正与传统计算方法相一致，具体为公式(1)(5)。改正后结果为四列数据，即测点x、y、h三
维坐标和gfi改正重力值。
[0124]
2、网格数据分层准备
[0125]
将网格文件按一定的垂向间隔分层，根据正演计算固定网格间距下重力值随深度的变化率，根据变化率可将全区dem高程按照全区高程最低点以下400m分为上下两个大层，下一大层再按最低高程400-500m以50m间隔分两层，500-700m以100m间隔分两层，700-1100m以200m间隔分两层，1100-1500m以400m间隔分一层，1500m至绝对0m标高按实际剩余间隔分为一层；上一大层以20m间隔连续分层至全区最高点，下一大层即所有的网格点均被采用，而上一大层出地表后，实际地表网格高程低于该划分层高程的网格点将被舍弃。分层过程中计算得出网格点该层的中心点高程zm、该网格立方体的半高c。
[0126]
3、地质界面文件准备
[0127]
实际密度值依照重力勘探区的地质图或其他钻孔等分岩性采集测量统计，得出每个地质体的统计密度值，即地质体名称和统计密度值。可以分以下两步准备地质界面文件。
[0128]
简化地质体的边界，统计区内地质体的各界面信息，界面信息有地质体边界点1和点2的x1、y1、z1、x2、y2、z2、倾角a、体与面的关系j(体位于面上、下)和走向(走向以两点的x、y坐标计算)与倾向的加减关系l，准备好各面的上述信息，以多个面圈闭地质体。
[0129]
上述也为三维地质直立长方体模型的建立过程，地质图要做适当的删减和取舍，离重力测点距离近处可详细分层，距离远处则粗略分层，将地质界线位置坐标、产状或钻孔位置、深度等数字化后，编写程序完成上述地质界面文件准备。
[0130]
4、正演数据准备
[0131]
以第3步的各地质界面圈闭分割网格分层数据，分割后形成属于该地质体的网格数据，即属于该体的所有直立长方体的数据，并赋予该地质体的密度σ。
[0132]
准备后每个地质体的直立长方体数据为5列，即x、y、zm、c、σ。
[0133]
5、直立长方体正演计算获得中间层物质影响值
[0134]
上述三步即形成最终用于直立长方体正演计算的地质体数据文件，结合测点重力值文件中的x、y、h，就可以利用直立长方体正演公式(10)进行各地质体数据文件中单个直立长方体对某一测点的正演计算，计算完成后将所有直立长方体的正演值相加，即为该地质体对某一重力测点中间层物质影响值。
[0135]
依照上述步骤循环完成所有测点的、所有地质体的中间层物质影响值计算，完成后结果文件包含四列，即测点的x、y、h和g
中
。
[0136]
6、精细布格重力异常计算
[0137]
将第一步计算获得的gfi减去第五步的g
中
，即获得重力测点的精细布格重力异常值，结果为四列，即测点的x、y、h和gb。
[0138]
该发明已在《金川铜镍矿床深边部及外围找矿勘查技术方法有效性研究》项目的ⅲ矿区大比例尺重力勘查中获得了应用。为了更好地解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。
[0139]
重力剖面点距50m，剖面长度1km，共21个测点。其测点周围1km以内的地形情况(网格间距10m的dem数据)如图7、地质简图如图8，地质剖面情况如图9，准备21个测点的x、y、h和g
观
，重力测点部分数据如表1，dem网格数据如表2，地质信息见地质简图(图8)和地质剖面情况图(图9)；
[0140]
表1重力测点观测数据及gfi重力改正值
[0141]
序号点号x(度)y(度)h(m)g
观
(mgal)gfi(mgal)116638.48626155102.13099041708.006979967.294458.623216738.48669978102.13128111690.140979970.880456.661
……………………………………
2118638.49335622102.13816001599.376979986.343443.548
[0142]
表2 dem网格数据文件格式
[0143][0144]
第一步：利用公式(1)和(5)进行正常场改正和高度改正，改正后的gfi重力改正值如表1中的最右列。
[0145]
第二步：如表2，以相邻网格点前两列坐标计算dem网格间距为以边长为10m的直立长方体，经变化深度和分层间隔试验确定，距离测点400m以外20m的分层间隔所对应的直立长方体(长宽高10m
×
10m
×
20m)，计算得出的重力影响值小于0.001
×
10-5
m/s2(同mgal)，dem网格数据高程(第三列)最大高程为1879.658，最小高程为1391.712，以20的倍数向零取整，则网格点地表最大高程层为1860，最小高程层为1380，先以最小高程下400m面，即海拔980面，将dem数据划分为上下两个大层，因上一大层中均有距离(垂向深度)小于400m的网格点，因此，海拔1880-980以20为分层间隔，分为1860、1840、1820
……
980共45层，选择dem网格点高程大于分层高程的点，以dem实际高程和分层高程之和除以2，计算组成的直立长方体的中点高程，以dem实际高程和分层高程之差除以2，计算组成的直立长方体的半高。
[0146]
980以下的分层原则与以上的相同，980-880以50m间隔分为930和880两层，880-680以100m间隔分为780和680两层，680-280以200m间隔分为480和280两层，280-0m为最后一层。直立长方体中点高程和直立长方体半高的计算方法与上一大层中使用的计算方法相一致。
[0147]
分层完成后，每一层均为x、y、zm和c四列，将所有的分层数据统一至一个文件内，完成后高程1480-1500m层内网格的分布情况如图10。
[0148]
第三步：提取地质图上各地质界线的端点坐标、地质剖面图上的产状等，建立各地质界面，本次共建立9个地质体界面文件，体1为sc-qbs地质体，体2为gn-bm地质体，体3为mi地质体，体4为σⅲ地质体，体5为q地质体(海拔1460以上)，体6为gn-ba地质体(海拔1460-0m)，体7为ml地质体，体8为σi地质体，体9为gn-ba地质体；其中σⅲ地质体的地质界面文件如表3，利用9个地质体的地质界面文件圈闭分割dem网格分层数据，匹配9个地质体的实测统计密度值，形成9个地质体的最终正演计算数据文件，其中σⅲ地质体的正演数据文件(即直立长方体文件)如图11。
[0149]
表3σⅲ地质体的地质界面文件
[0150][0151]
第四步：利用重力测点三维坐标x、y、h和第三步准备的各地质体圈闭的直立长方体数据，根据公式(10)分别正演各重力测点平面距离1000m以内所有直立长方体的正演值，求和后即为所有地质体对该测点的中间层物质影响值，循环完成所有测点的中间层物质影响值计算，即为g
中
。
[0152]
第五步：将第一步获得的gfi减去上述正演计算后的g
中
，则为精细计算获得的布格重力异常值。
[0153]
应当注意，本发明的实施方式可以通过硬件、软件或者软件和硬件的结合来实现。硬件部分可以利用专用逻辑来实现；软件部分可以存储在存储器中，由适当的指令执行系统，例如微处理器或者专用设计硬件来执行。本领域的普通技术人员可以理解上述的设备和方法可以使用计算机可执行指令和/或包含在处理器控制代码中来实现，例如在诸如磁盘、cd或dvd-rom的载体介质、诸如只读存储器(固件)的可编程的存储器或者诸如光学或电子信号载体的数据载体上提供了这样的代码。本发明的设备及其模块可以由诸如超大规模集成电路或门阵列、诸如逻辑芯片、晶体管等的半导体，或者诸如现场可编程门阵列、可编程逻辑设备等的可编程硬件设备的硬件电路实现，也可以用由各种类型的处理器执行的软件实现，也可以由上述硬件电路和软件的结合例如固件来实现。
[0154]
本发明在金川矿床ⅲ矿区的大比例尺重力勘探中，取得了良好的效果，描述如下：
[0155]
为了展示勘探效果，将精细计算的布格重力异常利用一次趋势分析，求取剩余重力异常后，对比分析传统方法与改进方法的勘探效果，如图12。
[0156]
如图8地表的超基性岩体大致位于点号173至181间，南侧倾角60
°
，北侧产状75
°
，本发明计算的如图12，红线为改进的方法计算获得的剩余重力异常，与蓝色线的传统计算方法相比，异常均起始于170点，但传统方法计算的异常于176点结束，而改进方法计算的异常延伸至180点，改进计算的异常范围和形态更加符合实际地质情况；传统方法在北侧计算得出一假异常显示，也与实际地质情况不符。从上图可以看出，改进的方法勘探效果明显提升。
[0157]
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所做的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。