本发明属于重力勘探技术领域，涉及一种通过航空重力数据物性层析反演计算方法，特别是涉及一种可以应用于大数据计算框架spark的航空重力数据物性层析反演计算方法。
背景技术：
重力勘探是一种传统的地球物理勘探方法，是矿产资源和地质构造研究的主要技术手段。例如，《重力勘探及其应用》(科技创新导报，2013.08，第240页)公开了一种重力勘探方法，通过对地表采集的重力异常数据进行数学计算，可以反演出地下不同深度的岩石密度分布。传统的重力反演方法是通过将采集的数据和模型计算的数据的差值最小二乘意义下最小化建立目标函数，进而反演求解。该方法通常面临欠定方程求解问题，解不稳定、计算量巨大、计算结果精细度低，难以将大面积区域采集的大规模航空重力数据进行反演计算。
技术实现要素：
本发明的目的在于克服现有技术的缺陷，提供一种计算步骤简单、计算过程稳定、计算结果准确度高，特别是可以应用于大数据计算的航空重力数据物性层析处理方法
为实现上述目的，本发明采用了如下技术方案：
一种应用于大数据计算的航空重力数据物性层析计算方法，包括以下步骤：
a1：通过按照地表网格采集到的航空重力数据d1，计算概率函数ηq：
其中，n为航空重力数据d1中数据采集点的数量，δg(x,y,z)为航空重力数据d1中位于地表(x,y,z)采集点的重力异常测量值，δgu(x,y,z)为地表(x,y,z)采集点重力异常的格林函数；
a2：计算地下网格的密度mq：
mq＝δρ·ηq(2)
其中，δρ为迭代计算的步进量；
a3：按照地表网格计算重力正演数据d2；
其中，δg(x,y,z)为重力正演数据d2中位于地表(x,y,z)采集点的重力异常正演值，g为万有引力常数，q为地下网格剖分形成的小立方体的数量，(xq,yq,zq)为地下网格剖分形成的某个小立方体的中心点坐标，δρq是所述小立方体的密度，δvq是所述小立方体的体积；
a4：对重力正演数据与航空重力数据之间的差值d2-d1，应用公式(1)计算得到概率函数偏差ηq′；
a5：更新地下网格密度mq：
mq＝mq+δρ·ηq′(4)
a6：重复步骤a3至a5，直到重力正演数据d2与航空重力数据d1均方误差恒定或最小，输出地下网格密度mq。
本发明提出的一种应用于大数据计算的航空重力数据物性层析计算方法，一是将重力层析反演问题由传统的欠定方程求解方法，转变成为正演迭代计算方法，通过对概率函数的迭代式密度转化实现航空重力数据的物性分布计算，该方法计算步骤简单、计算过程稳定、计算结果准确度高；二是上述迭代计算和复杂多路处理过程非常适用于大数据分布式并行计算框架spark，从而提升运算速度，可达到了接近实时的处理效果。
附图说明
图1是采用传统的单cpu计算方法与采用实施例2所述方法的计算时间对比图；
图2是采用实施例2所述方法针对数据模型的三维反演效果叠合示意图。
具体实施方式
以下结合附图1，进一步说明本发明一种应用于大数据计算的航空重力数据物性层析计算方法的具体实施方式。本发明一种应用于大数据计算的航空重力数据物性层析计算方法不限于以下实施例的描述。
实施例1：
本实施例给出一种应用于大数据计算的航空重力数据物性层析计算方法的计算原理及步骤。
重力层析反演的基本原理是将地下反演空间剖分成大小相等的小立方体，在每个立方体内，认为密度恒定。重力层析反演公式的推导要从重力正演计算公式出发，重力正演计算公式为：
其中，δg(x,y,z)为地表(x,y,z)处的重力异常值，g为万有引力常数，(xq,yq,zq)为地下网格剖分形成的某个小立方体的中心点坐标，δρq为所述小立方体的密度，δvq为所述小立方体的体积，q为地下网格剖分形成的小立方体的数量，该公式表示地下q个小立方体对地表(x,y,z)处产生的重力异常贡献的和(即测量或计算获得的重力数据)。
对于单位质量的小立方体，其δρqδvq为单位值，其产生的重力数据可表示为：
公式(5)是重力异常的格林函数。
结合公式(3)和(5)，可建立如下的方程：
其中，n为重力数据中数据采集点的数量，ηq为概率函数，表示地下某个小立方体在对地表重力异常(即重力数据)的贡献值的概率。
对于按照地表网格采集的航空重力数据，通过公式(1)进行迭代计算，可以将概率函数ηq转化成为地下网格的密度ρq，具体步骤如下：
a1：通过按照地表网格采集到的航空重力数据d1，计算概率函数ηq：
其中，n为航空重力数据d1中数据采集点的数量，δg(x,y,z)为航空重力数据d1中位于地表(x,y,z)采集点的重力异常测量值，δgu(x,y,z)为地表(x,y,z)采集点重力异常的格林函数；
a2：计算地下网格的密度mq：
mq＝δρ·ηq(2)
其中，δρ为迭代计算的步进量，为一个预设的密度值；当δρ取值较小时，计算结果将更为精准，但是需要进行更多的迭代计算次数；
a3：按照地表网格计算重力正演数据d2；
其中，δg(x,y,z)为重力正演数据d2中位于地表(x,y,z)采集点的重力异常正演值，g为万有引力常数，q为地下网格剖分形成的小立方体的数量，(xq,yq,zq)为地下网格剖分形成的某个小立方体的中心点坐标，δρq是所述小立方体的密度，δvq是所述小立方体的体积；
a4：对重力正演数据与航空重力数据之间的差值d2-d1，应用公式(1)计算得到概率函数偏差ηq′；
a5：更新地下网格密度mq：
mq＝mq+δρ·ηq′(4)
a6：重复步骤a3至a5，直到重力正演数据d2与航空重力数据d1均方误差恒定或最小，输出地下网格密度mq。
本实施例给出的一种应用于大数据计算的航空重力数据物性层析计算方法，一是计算步骤简单、计算过程稳定、计算结果准确度高；二是将重力层析反演问题由传统的欠定方程求解方法，转变成为正演迭代计算方法，非常适用于大数据计算框架spark进行运算加速。
实施例2：
本实施例给出一种采用大数据计算框架spark进行实施例1所述航空重力数据物性层析计算方法的具体实施方式。
apachespark是专为大规模数据处理而设计的快速通用的计算引擎。spark是ucberkeleyamplab(加州大学伯克利分校的amp实验室)所开源的类似hadoopmapreduce的通用并行框架。spark可以利用dag进行调度执行规划，可以处理多步复杂的作业，并且对于每一个任务启动一个线程，而不是启动进程，大大缩短了任务启动时间。实施例1中所述航空重力数据物性层析计算方法可通过大数据框架spark进行加速计算，其中的迭代计算非常适合spark处理，因此可以获得非常好的加速效果。具体步骤如下：
b1：对地下网格剖分形成的小立方体总数进行切分，构建rdd，并持久化到内存；
b2：通过按照地表网格采集到的航空重力数据d1，计算概率函数ηq：
其中，n为航空重力数据d1中数据采集点的数量，δg(x,y,z)为航空重力数据d1中位于地表(x,y,z)采集点的重力异常测量值，δgu(x,y,z)为地表(x,y,z)采集点重力异常的格林函数；
b3：计算地下网格的密度mq：
mq＝δρ·ηq(2)
b4：迭代计算地表网格重力正演数据d2和地下网格密度mq，具体步骤如下：
b4.1：对rdd进行map算子，按照地表网格计算各自的重力正演数据d2；
其中，δg(x,y,z)为重力正演数据d2中位于地表(x,y,z)采集点的重力异常正演值，g为万有引力常数，q为地下网格剖分形成的小立方体的数量，(xq,yq,zq)为地下网格剖分形成的某个小立方体的中心点坐标，δρq是所述小立方体的密度，δvq是所述小立方体的体积；
b4.2：对rdd进行reduce算子，合并重力正演数据d2；
b4.3对重力正演数据与航空重力数据之间的差值d2-d1，应用公式(1)计算得到概率函数偏差ηq′；
b4.4更新地下网格密度mq：
mq＝mq+δρ·ηq′(4)
b4.5：重复步骤b4.1至b4.4，直到重力正演数据d2与航空重力数据d1均方误差恒定或最小；
b5：输出地下网格密度mq。
本实施例给出的一种基于大数据计算的航空重力数据物性层析处理方法，采用大数据框架spark，通过分布式内存动态数据集rdd实现内存数据共享，过程中生成的正演计算结果保存在内存不是将其写入磁盘，同时使用成本更低的shuffle操作对重力数据计算结果进行合并，有效减少磁盘io的时间、提升了运算效率。
实施例3：
本实施例给出实施例2所述实施方式的计算效率的仿真结果。
采用具有18个节点的cpu计算集群对实施例2所述实施方式进行测试，并通过相同的单cpu计算进行对比。每个节点搭配intelxeone5-2680v2cpu，浮点计算能力为211.2gflops，ddr3内存120g，节点间通过infiniband连接。
如图1所示，是单cpu计算和spark框架下并行计算的计算时间对比。从图中可以看出，随着地下网格规模的增加，spark框架下并行计算具有明显的优势。
如图2所示，是针对数据模型的三维反演效果叠合示意图。其中，a、b为纵向剖面，c、d为水平切片；图片中黑色线框是实际的数据模型，灰度色块是反演的结果。从不同剖面的叠合示意图中可以看出，本方法具有非常好的反演效果。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。