CN115220119B - 一种适用于大规模数据的重力反演方法 - Google Patents

Abstract

本发明公开一种可适用于大规模数据的重力反演方法及处理终端，所述方法包括：步骤1：引入深度加权矩阵构建重力反演的目标函数；步骤2：采用共轭梯度算法求解目标函数的梯度，并且在求解过程中，引入惩罚函数且按相应规则，在迭代过程中求解正则化参数；步骤3：采用快速法计算出的灵敏度矩阵来替代前述步骤中的灵敏度矩阵，以降低对内存的要求，在灵敏度矩阵与其他向量相乘过程中，采用灵敏度矩阵的每一行与向量来相乘。本发明降低了重力反演对内存资源的要求，提高了计算效率，能够很好地适应于大规模数据的重力反演中。

Description

一种适用于大规模数据的重力反演方法

技术领域

本发明涉及重力反演技术领域，具体是一种适用于大规模数据的重力反演方法。

背景技术

重力反演构建三维地质模型，是地质工作者了解地球深部构造、认知地下结构的重要手段。尤其是，面向调查研究程度低地区，比如深海远洋等，采用移动平台重力测量快速获取数据，随后对获取的数据进行处理和解释，能够进一步判明地下地质结构、识别深部构造特征、寻找隐伏矿床等，提高对该地区深部结构的认识，为地质调查工作的深入开展提供服务。同时，随着地质调查研究工作的深入，对三维密度模型精细化程度提出了更高的要求，也就要求更大比例尺的观测数据和更加精细的关于重力密度分布的地下网格模型，基于重力数据构建大范围、高分辨率三维密度模型，涉及的数据量大、计算量大，对计算机内存要求高。在现有的重力反演方法中，往往因计算机硬件资源(例如内存容量)受限制，导致对大规模数据的重力反演无法进行，并且即使可以进行也因数据量庞大而导致反演速率慢。

发明内容

针对现有技术的不足，本发明的目的之一是提供一种可适用于大规模数据的重力反演方法，其能够适用于大规模数据的重力反演的问题；

本发明的目的之二是提供一种处理终端，其能够适用于大规模数据的重力反演的问题。

实现本发明的目的之一的技术方案为：一种可适用于大规模数据的重力反演方法，包括如下步骤：

步骤1：按公式①建立重力反演的目标函数φ_w：

式中，W表示深度加权矩阵，φ_d(m)表示关于m的拟合函数，φ_s(m) 表示关于m的稳定函数，m表示地下密度分布，m＝{m_i}，l＜m_i＜u，m_i表示m中的第i个元素，l和u表示m的取值范围，d表示观测数据，α表示正则化参数，

表示二范数的平方，

对加权矩阵W中第j个元素w_j的表达式如下：

A_ij表示灵敏度矩阵A中的元素，i表示观测点的编号，j表示地下网格的编号，β表示常数；

步骤2：采用共轭梯度算法求出目标函数φ_w的梯度

进一步地，β＝2。

进一步地，共轭梯度算法的具体包括如下步骤：

(1)设定最大迭代次数ite_num，设置j＝0和初始模型m₀，计算出m_w(0)＝W^-1m₀和

(2)第一次迭代设置搜索方向为L₀＝-I₀；

(3)计算步长k，

(4)计算m_w(j+1)＝m_w(j)+kL_j，m_j+1＝Wm_w(j+1)，

(5)若r小于预设的阈值或者当j＞i te_num，则停止计算，否则，设置j＝j+1并继续执行步骤(6)；

(6)计算

新的搜索方向为 L_j＝-I_j+βL_j-1，然后跳转到步骤(3)。

进一步地，采用惩罚函数来对每次迭代得到的模型参数m进行约束，惩罚函数采用如下表达式：

进一步地，若初始模型为0，则设置α₀＝0，否则，按下式计算：

在随后的迭代中，α的取值逐渐减小，其表达式如下：

α_k＝α_k-1*q，0＜q＜1 。

进一步地，执行完步骤2之后，还包括，

步骤3：采用快速法计算出的灵敏度矩阵V_a，并用V_a替代上述灵敏度矩阵A，

快速法计算出的灵敏度矩阵V_a的具体步骤包括：

对于重力数据，引入一个三元组(a_n b_n c_n)，n＝1,2，其表达式为

其中，(x,y,z)表示观测点的坐标，[x₁,x₂]、[y₁,y₂]、[z₁,z₂]表示地下网格取值范围，地下网格在观测点引起的重力异常的表达式为

其中，γ＝6.67×10^-11m³/(kg·s²)，ρ表示地下异常体密度，

观测点坐标已知情况下，g_z是关于变量a_n、b_n、c_n的表达式，

g_z＝f(a_n,b_n,c_n)，n＝1,2

将地下规则划分为大小为d_x×d_y×d_z的长方体，则有

由此，g_z是包含三个变量的函数式，其表达式为

g_z＝f(a₁,b₁,c₁)

在空间坐标系下，竖直向下为正，将地下划分为n_x×n_y×n_z个小长方体，测区的取值范围为[x_min,x_min+n_xd_x]，[y_min,y_min+n_yd_y]，[0,n_zd_z]，观测点规则分布于同一高度平面z＝z_l，其在地面投影位于每个网格的中心，观测点在X方向坐标取值为

观测点在Y方向坐标取值为

观测点在Z方向坐标取值为

①统计a₁,b₁,c₁的值，令a₁＝(x-ξ₁)，b₁＝(y-η₁)，c₁＝(z-ζ₁)，ξ₁、η₁、ζ₁分别表示任一小长方体在X,Y,Z三方向取值的下限，则

由观测点坐标取值和上式，可得到a₁,b₁,c₁的值，用向量A₁、B₁、 C₁表示，将其中的数值从小到大排列，可得

C₁＝(z-(n_z-1)dz,...,z-dz,z)^T

②依据向量A₁、B₁、C₁的值的个数确定V_a矩阵的大小，A₁中包含2n_x-1个值，B₁中包含2n_y-1个值，C₁中包含n_z个值，则V_a大小为 N_k×1，且

N_k＝(2n_x-1)×(2n_y-1)×n_z

③将A₁、B₁、C₁中的值依次取出，计算获得V_a，V_a中包含所有构成灵敏度矩阵A的值，

用V_a替代上述灵敏度矩阵A的具体实现过程包括：

建立灵敏度矩阵A中元素和V_a的关系，对A₁中数值进行编号，记为i_x＝1,…,2n_x-1，对B₁中数值进行编号有i_y＝1,…,2n_y-1，对C₁中数值进行编号有i_z＝1,…,n_z，对V_a中数值进行编号有 i_k＝1,…,(2n_x-1)×(2n_y-1)×n_z，灵敏度矩阵A中的元素为

其中i_obs表示测点编号，j_g表示网格编号，对应的观测点坐标为

对应的小长方体的在三方向取值下限为

由此，计算得到其对应的变量在A₁、B₁、C₁中的位置，进而所得结果在V_a中的位置和

的值，

i_k＝(i_x-1)×(2n_y-1)×n_z+(i_y-1)×n_z+i_z

利用V_a，对反演过程中涉及灵敏度矩阵A的计算进行优化，反演计算中涉及灵敏度矩阵A的计算包括两类：一是A与向量相乘，记为 AV，V表示反演计算中与A相乘的向量，二是A的转置乘A乘向量，记为A^TAV，

对于第一种计算，读取V_a，利用V_a与A中元素的关系，实时从V_a中提取元素构建A中的行向量，与向量V相乘，逐一获得结果向量中的元素，从而避免直接读取和操作A，

对于第二种计算，利用AV计算结果，采用上述相同操作计算 A^TM。

实现本发明的目的之二的技术方案为：一种处理终端，其包括：

存储器，用于存储程序指令；

处理器，用于运行所述程序指令，以执行所述的可适用于大规模数据的重力反演方法的步骤。

本发明的有益效果为：本发明通过将快速计算灵敏度矩阵的方法引入到重力反演中，降低了重力反演对内存资源的要求，且将并行计算方式引入到灵敏度矩阵与向量相乘的计算中，提高了计算效率，能够很好地适应于大规模数据的重力反演中。

附图说明

图1为本发明的流程示意图；

图2为处理终端的结构示意图。

具体实施方案

下面，结合附图以及具体实施方案，对本发明做进一步描述：

如图1所示，一种可适用于大规模数据的重力反演方法，包括如下步骤：

步骤1：按公式①建立重力反演的目标函数φ_w：

式中，W表示深度加权矩阵，φ_d(m)表示关于m的拟合函数，φ_s(m) 表示关于m的稳定函数，m表示地下密度分布，m＝{m_i}，l＜m_i＜u，m_i表示m中的第i个元素，l和u表示m的取值范围，d表示观测数据，α表示正则化参数，

表示二范数的平方。

在本步骤中，通过调整α的值来平衡拟合函数和稳定函数对目标函数φ_w的影响，若α较大，则稳定函数对目标函数φ_w的影响较大，若α较小，则拟合函数对目标函数φ_w的影响较大。

之所以引入深度加权矩阵，原因在于重力异常随着距离的增加衰减，重力反演的结构存在“趋肤效应”，通过引入深度加权矩阵可很好地克服这一影响。

对加权矩阵W中第j个元素w_j的表达式如下：

A_ij表示灵敏度矩阵A中的元素，i表示观测点的编号，j表示地下网格的编号，β表示常数，本实施例中取值为β＝2。

步骤2：采用共轭梯度算法求出目标函数φ_w的梯度

共轭梯度算法的具体包括如下步骤：

(1)设定最大迭代次数ite_num，设置j＝0和初始模型m₀，计算出m_w(0)＝W^-1m₀和

(2)第一次迭代设置搜索方向为L₀＝-I₀；

(3)计算步长k，

(4)计算m_w(j+1)＝m_w(j)+kL_j，m_j+1＝Wm_w(j+1)，

(5)若r足够小(例如小于预设的阈值)或者当j＞i te_num，则停止计算，否则，设置j＝j+1并继续执行步骤(6)；

(6)计算

新的搜索方向为 L_j＝-I_j+βL_j-1，然后跳转到步骤(3)。

在上述共轭梯度算法求解的迭代过程中，为了防止模型参数的值超出范围，采用惩罚函数来对每次迭代得到的模型参数m进行约束，惩罚函数采用如下表达式：

正则化参数α是重力反演计算中的关键参数，若初始模型为0，则设置α₀＝0，否则，按下式计算：

在随后的迭代中，α的取值逐渐减小，其表达式如下：

α_k＝α_k-1*q，0＜q＜1

步骤3：采用快速法计算灵敏度矩阵V_a，采用V_a替代上述灵敏度矩阵A，从而在实际计算过程中，只需要存储和读取V_a，避免直接读取灵敏度矩阵A，从而达到降低了对内存的要求。

同时，对于灵敏度矩阵A与其他向量相乘，实时构建灵敏度矩阵A的每一行，用灵敏度矩阵A的每一行与向量相乘，从而进一步降低对内存的要求，通过引入并行来计算灵敏度矩阵A与其他向量相乘，提高了计算效率。

采用快速法计算灵敏度矩阵V_a的具体步骤如下：

(1)快速计算V_a

对于重力数据，引入一个三元组(a_n b_n c_n)，n＝1,2表达式为

其中，(x,y,z)表示观测点的坐标，[x₁,x₂]、[y₁,y₂]、[z₁,z₂]表示地下网格取值范围，地下网格在观测点引起的重力异常的表达式为

其中，γ＝6.67×10^-11m³/(kg·s²)，ρ表示地下异常体密度，在计算灵敏度矩阵时，一般取单位密度1g/cm³,

观测点坐标已知情况下，g_z可以看做是关于变量a_n、b_n、c_n的表达式，

g_z＝f(a_n,b_n,c_n)，n＝1,2

一般将地下规则划分为大小为d_x×d_y×d_z的长方体，则有

由此，g_z是包含三个变量的函数式，其表达式为

g_z＝f(a₁,b₁,c₁)

在空间坐标系下，竖直向下为正，将地下划分为n_x×n_y×n_z个小长方体，测区的取值范围为[x_min,x_min+n_xd_x]，[y_min,y_min+n_yd_y]，[0,n_zd_z]，观测点规则分布于同一高度平面z＝z_l，其在地面投影位于每个网格的中心，观测点在X方向坐标取值为

观测点在Y方向坐标取值为

观测点在Z方向坐标取值为

①统计a₁,b₁,c₁的值。令a₁＝(x-ξ₁)，b₁＝(y-η₁)，c₁＝(z-ζ₁)，ξ₁、η₁、ζ₁分别表示任一小长方体在X,Y,Z三方向取值的下限，则

由观测点坐标取值和上式，可得到a₁,b₁,c₁的值，用向量A₁、B₁、 C₁表示，将其中的数值从小到大排列，可得

C₁＝(z-(n_z-1)dz,...,z-dz,z)^T

②依据向量A₁、B₁、C₁的值的个数确定V_a矩阵的大小。A₁中包含2n_x-1个值，B₁中包含2n_y-1个值，C₁中包含n_z个值，则V_a大小为 N_k×1，且

N_k＝(2n_x-1)×(2n_y-1)×n_z

③将A₁、B₁、C₁中的值依次取出，计算获得V_a。V_a中包含所有构成灵敏度矩阵A的值，同时，其对内存要求远远低于容纳灵敏度矩阵A的要求，比如，存储大小为10000*100000的灵敏度矩阵A，存储需要约5.36GB空间，而存储组成该矩阵的值的V_a，所需要的空间仅为1.6MB。

(2)将V_a代入到反演计算

①建立灵敏度矩阵A中元素和V_a的关系。对A₁中数值进行编号，记为i_x＝1,…,2n_x-1，对B₁中数值进行编号有i_y＝1,…,2n_y-1，对C₁中数值进行编号有i_z＝1,…,n_z，对V_a中数值进行编号有 i_k＝1,…,(2n_x-1)×(2n_y-1)×n_z。灵敏度矩阵A中的元素为

其中i_obs表示测点编号，j_g表示网格编号，对应的观测点坐标为

对应的小长方体的在三方向取值下限为

由此，计算得到其对应的变量在A₁、B₁、C₁中的位置，进而所得结果在V_a中的位置和

的值。

i_k＝(i_x-1)×(2n_y-1)×n_z+(i_y-1)×n_z+i_z

②利用V_a，对反演过程中涉及灵敏度矩阵A的计算进行优化。反演计算中涉及灵敏度矩阵A的计算包括两类：一是A与向量相乘，记为AV，V表示反演计算中与A相乘的向量，二是A的转置乘A乘向量，记为A^TAV。

对于第一种计算，一般是直接进行矩阵计算，但是针对大数据情况下，很容易出现内存不足，无法计算的情况。为了克服这一情况，直接读取V_a，利用V_a与A中元素的关系，实时从V_a中提取元素构建A 中的行向量，与向量V相乘，逐一获得结果向量中的元素。从而能够避免直接读取和操作A，降低了对内存的要求。同时，采用这种算法，会降低运算效率，增加运算时间，由于每一行向量与列向量相乘是相互独立的，所以采用并行计算实现这一过程，从而提高了计算效率。

对于第二种计算，利用AV计算结果，采用上述同样的思路计算 A^TM。

本发明通过将快速计算灵敏度矩阵的方法引入到重力反演中，降低了重力反演对内存资源的要求，且将并行计算方式引入到灵敏度矩阵与向量相乘的计算中，提高了计算效率，能够很好地适应于大规模数据的重力反演中。

本发明可很好地集成在海洋环境监测和探测设备上，实现对海洋区域的重力密度数据分别的探测，为后续判明地下地质结构、识别深部构造特征、寻找隐伏矿床等提供数据基础。

如图2所示，本发明还提供一种处理终端100，其包括：

存储器101，用于存储程序指令；

处理器102，用于运行所述程序指令，以执行所述可适用于大规模数据的重力反演方法的步骤。

本说明书所公开的实施例只是对本发明单方面特征的一个例证，本发明的保护范围不限于此实施例，其他任何功能等效的实施例均落入本发明的保护范围内。对于本领域的技术人员来说，可根据以上描述的技术方案以及构思，做出其它各种相应的改变以及变形，而所有的这些改变以及变形都应该属于本发明权利要求的保护范围之内。

Claims (7)

1.一种可适用于大规模数据的重力反演方法，其特征在于，包括如下步骤：

步骤1：按公式①建立重力反演的目标函数φ_w：

式中，W表示深度加权矩阵，φ_d(m)表示关于m的拟合函数，φ_s(m)表示关于m的稳定函数，m表示地下密度分布，m＝{m_i}，l＜m_i＜u，m_i表示m中的第i个元素，l和u表示m的取值范围，d表示观测数据，α表示正则化参数，

表示二范数的平方，

对加权矩阵W中第j个元素w_j的表达式如下：

A_ij表示灵敏度矩阵A中的元素，i表示观测点的编号，j表示地下网格的编号，β表示常数；

步骤2：采用共轭梯度算法求出目标函数φ_w的梯度

2.根据权利要求1所述的可适用于大规模数据的重力反演方法，其特征在于，β＝2。

3.根据权利要求1所述的可适用于大规模数据的重力反演方法，其特征在于，共轭梯度算法的具体包括如下步骤：

(1)设定最大迭代次数ite_num，设置j＝0和初始模型m₀，计算出m_w(0)＝W^-1m₀和

(2)第一次迭代设置搜索方向为L₀＝-I₀；

(3)计算步长k，

(4)计算m_w(j+1)＝m_w(j)+kL_j，m_j+1＝Wm_w(j+1)，

(5)若r小于预设的阈值或者当j＞ite_num，则停止计算，否则，设置j＝j+1并继续执行步骤(6)；

(6)计算

新的搜索方向为L_j＝-I_j+βL_j-1，然后跳转到步骤(3)。

4.根据权利要求3所述的可适用于大规模数据的重力反演方法，其特征在于，采用惩罚函数来对每次迭代得到的模型参数m进行约束，惩罚函数采用如下表达式：

5.根据权利要求4所述的可适用于大规模数据的重力反演方法，其特征在于，若初始模型为0，则设置α₀＝0，否则，按下式计算：

在随后的迭代中，α的取值逐渐减小，其表达式如下：

α_k＝α_k-1*q，0＜q＜1。

6.根据权利要求1所述的可适用于大规模数据的重力反演方法，其特征在于，执行完步骤2之后，还包括，

步骤3：采用快速法计算出的灵敏度矩阵V_a，并用V_a替代上述灵敏度矩阵A，

快速法计算出的灵敏度矩阵V_a的具体步骤包括：

对于重力数据，引入一个三元组(a_n b_n c_n)，n＝1,2，其表达式为

其中，(x,y,z)表示观测点的坐标，[x₁,x₂]、[y₁,y₂]、[z₁,z₂]表示地下网格取值范围，地下网格在观测点引起的重力异常的表达式为

其中，γ＝6.67×10^-11m³/(kg·s²)，ρ表示地下异常体密度，

s₁＝-1,s₂＝+1，观测点坐标已知情况下，g_z是关于变量a_n、b_n、c_n的表达式，

g_z＝f(a_n,b_n,c_n)，n＝1,2

将地下规则划分为大小为d_x×d_y×d_z的长方体，则有

由此，g_z是包含三个变量的函数式，其表达式为

g_z＝f(a₁,b₁,c₁)

在空间坐标系下，竖直向下为正，将地下划分为n_x×n_y×n_z个小长方体，测区的取值范围为[x_min,x_min+n_xd_x]，[y_min,y_min+n_yd_y]，[0,n_zd_z]，观测点规则分布于同一高度平面z＝z_l，其在地面投影位于每个网格的中心，观测点在X方向坐标取值为

i＝1,…,n_x，观测点在Y方向坐标取值为

j＝1,…,n_y，观测点在Z方向坐标取值为

k＝1,…,n_z，

①统计a₁,b₁,c₁的值，令a₁＝(x-ξ₁)，b₁＝(y-η₁)，c₁＝(z-ζ₁)，ξ₁、η₁、ζ₁分别表示任一小长方体在X,Y,Z三方向取值的下限，则

由观测点坐标取值和上式，可得到a₁,b₁,c₁的值，用向量A₁、B₁、C₁表示，将其中的数值从小到大排列，可得

C₁＝(z-(n_z-1)dz,...,z-dz,z)^T

②依据向量A₁、B₁、C₁的值的个数确定V_a矩阵的大小，A₁中包含2n_x-1个值，B₁中包含2n_y-1个值，C₁中包含n_z个值，则V_a大小为N_k×1，且

N_k＝(2n_x-1)×(2n_y-1)×n_z

③将A₁、B₁、C₁中的值依次取出，计算获得V_a，V_a中包含所有构成灵敏度矩阵A的值，

用V_a替代上述灵敏度矩阵A的具体实现过程包括：

建立灵敏度矩阵A中元素和V_a的关系，对A₁中数值进行编号，记为i_x＝1,…,2n_x-1，对B₁中数值进行编号有i_y＝1,…,2n_y-1，对C₁中数值进行编号有i_z＝1,…,n_z，对V_a中数值进行编号有i_k＝1,…,(2n_x-1)×(2n_y-1)×n_z，灵敏度矩阵A中的元素为

其中i_obs表示测点编号，j_g表示网格编号，对应的观测点坐标为

对应的小长方体的在三方向取值下限为

由此，计算得到其对应的变量在A₁、B₁、C₁中的位置，进而所得结果在V_a中的位置和

的值，

i_k＝(i_x-1)×(2n_y-1)×n_z+(i_y-1)×n_z+i_z

②利用V_a，对反演过程中涉及灵敏度矩阵A的计算进行优化，反演计算中涉及灵敏度矩阵A的计算包括两类：一是A与向量相乘，记为AV，V表示反演计算中与A相乘的向量，二是A的转置乘A乘向量，记为A^TAV，

对于第一种计算，读取V_a，利用V_a与A中元素的关系，实时从V_a中提取元素构建A中的行向量，与向量V相乘，逐一获得结果向量中的元素，从而避免直接读取和操作A，

对于第二种计算，利用AV计算结果，采用上述相同操作计算A^TM。

7.一种处理终端，其特征在于，其包括：

存储器，用于存储程序指令；

处理器，用于运行所述程序指令，以执行如权利要求1-6任一项所述的可适用于大规模数据的重力反演方法的步骤。

Images