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

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


背景技术：

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


技术实现要素：

3.针对现有技术的不足，本发明的目的之一是提供一种可适用于大规模数据的重力反演方法，其能够适用于大规模数据的重力反演的问题；
4.本发明的目的之二是提供一种处理终端，其能够适用于大规模数据的重力反演的问题。
5.实现本发明的目的之一的技术方案为：一种可适用于大规模数据的重力反演方法，包括如下步骤：
6.步骤1：按公式
①
建立重力反演的目标函数φw：
[0007][0008]
式中，w表示深度加权矩阵，φd(m)表示关于m的拟合函数，φs(m) 表示关于m的稳定函数，m表示地下密度分布，m＝{mi}，l＜mi＜u，mi表示m中的第i个元素，l和u表示m的取值范围，d表示观测数据，α表示正则化参数，表示二范数的平方，
[0009]
对加权矩阵w中第j个元素wj的表达式如下：
[0010][0011]aij
表示灵敏度矩阵a中的元素，i表示观测点的编号，j表示地下网格的编号，β表示常数；
[0012]
步骤2：采用共轭梯度算法求出目标函数φw的梯度
[0013]
进一步地，β＝2。
[0014]
进一步地，共轭梯度算法的具体包括如下步骤：
[0015]
(1)设定最大迭代次数ite_num，设置j＝0和初始模型m0，计算出m
w(0)
＝w-1
m0和
[0016]
(2)第一次迭代设置搜索方向为l0＝-i0；
[0017]
(3)计算步长k，
[0018]
(4)计算m
w(j+1)
＝m
w(j)
+klj，m
j+1
＝wm
w(j+1)
，
[0019]
(5)若r小于预设的阈值或者当j＞i te_num，则停止计算，否则，设置j＝j+1并继续执行步骤(6)；
[0020]
(6)计算新的搜索方向为 lj＝-ij+βl
j-1
，然后跳转到步骤(3)。
[0021]
进一步地，采用惩罚函数来对每次迭代得到的模型参数m进行约束，惩罚函数采用如下表达式：
[0022][0023]
进一步地，若初始模型为0，则设置α0＝0，否则，按下式计算：
[0024][0025]
在随后的迭代中，α的取值逐渐减小，其表达式如下：
[0026]
αk＝α
k-1
*q，0＜q＜1 。
[0027]
进一步地，执行完步骤2之后，还包括，
[0028]
步骤3：采用快速法计算出的灵敏度矩阵va，并用va替代上述灵敏度矩阵a，
[0029]
快速法计算出的灵敏度矩阵va的具体步骤包括：
[0030]
对于重力数据，引入一个三元组(a
n b
n cn)，n＝1,2，其表达式为
[0031][0032]
其中，(x,y,z)表示观测点的坐标，[x1,x2]、[y1,y2]、[z1,z2]表示地下网格取值范围，地下网格在观测点引起的重力异常的表达式为
[0033][0034]
其中，γ＝6.67
×
10-11
m3/(kg
·
s2)，ρ表示地下异常体密度，
观测点坐标已知情况下，gz是关于变量an、bn、cn的表达式，
[0035]gz
＝f(an,bn,cn)，n＝1,2
[0036]
将地下规则划分为大小为d
x
×dy
×dz
的长方体，则有
[0037][0038]
由此，gz是包含三个变量的函数式，其表达式为
[0039]gz
＝f(a1,b1,c1)
[0040]
在空间坐标系下，竖直向下为正，将地下划分为n
x
×
ny×
nz个小长方体，测区的取值范围为[x
min
,x
min
+n
xdx
]，[y
min
,y
min
+n
ydy
]，[0,n
zdz
]，观测点规则分布于同一高度平面z＝z
l
，其在地面投影位于每个网格的中心，观测点在x方向坐标取值为观测点在y方向坐标取值为观测点在z方向坐标取值为
[0041]
①
统计a1,b1,c1的值，令a1＝(x-ξ1)，b1＝(y-η1)，c1＝(z-ζ1)，ξ1、η1、ζ1分别表示任一小长方体在x,y,z三方向取值的下限，则
[0042][0043]
由观测点坐标取值和上式，可得到a1,b1,c1的值，用向量a1、b1、 c1表示，将其中的数值从小到大排列，可得
[0044][0045][0046]
c1＝(z-(n
z-1)dz,...,z-dz,z)
t
[0047]
②
依据向量a1、b1、c1的值的个数确定va矩阵的大小，a1中包含2n
x-1个值，b1中包含2n
y-1个值，c1中包含nz个值，则va大小为 nk×
1，且
[0048]
nk＝(2n
x-1)
×
(2n
y-1)
×
nz[0049]
③
将a1、b1、c1中的值依次取出，计算获得va，va中包含所有构成灵敏度矩阵a的值，
[0050]
用va替代上述灵敏度矩阵a的具体实现过程包括：
[0051]
建立灵敏度矩阵a中元素和va的关系，对a1中数值进行编号，记为i
x
＝1,
…
,2n
x-1，对b1中数值进行编号有iy＝1,
…
,2n
y-1，对c1中数值进行编号有iz＝1,
…
,nz，对va中数值进行编号有 ik＝1,
…
,(2n
x-1)
×
(2n
y-1)
×
nz，灵敏度矩阵a中的元素为其中i
obs
表示测
点编号，jg表示网格编号，对应的观测点坐标为对应的小长方体的在三方向取值下限为由此，计算得到其对应的变量在a1、b1、c1中的位置，进而所得结果在va中的位置和的值，
[0052][0053]
ik＝(i
x-1)
×
(2n
y-1)
×
nz+(i
y-1)
×
nz+iz[0054][0055]
利用va，对反演过程中涉及灵敏度矩阵a的计算进行优化，反演计算中涉及灵敏度矩阵a的计算包括两类：一是a与向量相乘，记为 av，v表示反演计算中与a相乘的向量，二是a的转置乘a乘向量，记为a
t
av，
[0056]
对于第一种计算，读取va，利用va与a中元素的关系，实时从va中提取元素构建a中的行向量，与向量v相乘，逐一获得结果向量中的元素，从而避免直接读取和操作a，
[0057][0058]
对于第二种计算，利用av计算结果，采用上述相同操作计算 a
t
m。
[0059]
实现本发明的目的之二的技术方案为：一种处理终端，其包括：
[0060]
存储器，用于存储程序指令；
[0061]
处理器，用于运行所述程序指令，以执行所述的可适用于大规模数据的重力反演方法的步骤。
[0062]
本发明的有益效果为：本发明通过将快速计算灵敏度矩阵的方法引入到重力反演中，降低了重力反演对内存资源的要求，且将并行计算方式引入到灵敏度矩阵与向量相乘的计算中，提高了计算效率，能够很好地适应于大规模数据的重力反演中。
附图说明
[0063]
图1为本发明的流程示意图；
[0064]
图2为处理终端的结构示意图。
具体实施方案
[0065]
下面，结合附图以及具体实施方案，对本发明做进一步描述：
[0066]
如图1所示，一种可适用于大规模数据的重力反演方法，包括如下步骤：
[0067]
步骤1：按公式
①
建立重力反演的目标函数φw：
[0068][0069]
式中，w表示深度加权矩阵，φd(m)表示关于m的拟合函数，φs(m) 表示关于m的稳定函数，m表示地下密度分布，m＝{mi}，l＜mi＜u，mi表示m中的第i个元素，l和u表示m的取值范围，d表示观测数据，α表示正则化参数，表示二范数的平方。
[0070]
在本步骤中，通过调整α的值来平衡拟合函数和稳定函数对目标函数φw的影响，若α较大，则稳定函数对目标函数φw的影响较大，若α较小，则拟合函数对目标函数φw的影响较大。
[0071]
之所以引入深度加权矩阵，原因在于重力异常随着距离的增加衰减，重力反演的结构存在“趋肤效应”，通过引入深度加权矩阵可很好地克服这一影响。
[0072]
对加权矩阵w中第j个元素wj的表达式如下：
[0073][0074]aij
表示灵敏度矩阵a中的元素，i表示观测点的编号，j表示地下网格的编号，β表示常数，本实施例中取值为β＝2。
[0075]
步骤2：采用共轭梯度算法求出目标函数φw的梯度
[0076][0077]
共轭梯度算法的具体包括如下步骤：
[0078]
(1)设定最大迭代次数ite_num，设置j＝0和初始模型m0，计算出m
w(0)
＝w-1
m0和
[0079]
(2)第一次迭代设置搜索方向为l0＝-i0；
[0080]
(3)计算步长k，
[0081]
(4)计算m
w(j+1)
＝m
w(j)
+klj，m
j+1
＝wm
w(j+1)
，
[0082]
(5)若r足够小(例如小于预设的阈值)或者当j＞i te_num，则停止计算，否则，设置j＝j+1并继续执行步骤(6)；
[0083]
(6)计算新的搜索方向为 lj＝-ij+βl
j-1
，然后跳转到步骤(3)。
[0084]
在上述共轭梯度算法求解的迭代过程中，为了防止模型参数的值超出范围，采用惩罚函数来对每次迭代得到的模型参数m进行约束，惩罚函数采用如下表达式：
[0085][0086]
正则化参数α是重力反演计算中的关键参数，若初始模型为0，则设置α0＝0，否则，按下式计算：
[0087][0088]
在随后的迭代中，α的取值逐渐减小，其表达式如下：
[0089]
αk＝α
k-1
*q，0＜q＜1
[0090]
步骤3：采用快速法计算灵敏度矩阵va，采用va替代上述灵敏度矩阵a，从而在实际计算过程中，只需要存储和读取va，避免直接读取灵敏度矩阵a，从而达到降低了对内存的要求。
[0091]
同时，对于灵敏度矩阵a与其他向量相乘，实时构建灵敏度矩阵a的每一行，用灵敏度矩阵a的每一行与向量相乘，从而进一步降低对内存的要求，通过引入并行来计算灵敏度矩阵a与其他向量相乘，提高了计算效率。
[0092]
采用快速法计算灵敏度矩阵va的具体步骤如下：
[0093]
(1)快速计算va[0094]
对于重力数据，引入一个三元组(a
n b
n cn)，n＝1,2表达式为
[0095][0096]
其中，(x,y,z)表示观测点的坐标，[x1,x2]、[y1,y2]、[z1,z2]表示地下网格取值范围，地下网格在观测点引起的重力异常的表达式为
[0097][0098]
其中，γ＝6.67
×
10-11
m3/(kg
·
s2)，ρ表示地下异常体密度，在计算灵敏度矩阵时，一般取单位密度1g/cm3,观测点坐标已知情况下，gz可以看做是关于变量an、bn、cn的表达式，
[0099]gz
＝f(an,bn,cn)，n＝1,2
[0100]
一般将地下规则划分为大小为d
x
×dy
×dz
的长方体，则有
[0101][0102]
由此，gz是包含三个变量的函数式，其表达式为
[0103]gz
＝f(a1,b1,c1)
[0104]
在空间坐标系下，竖直向下为正，将地下划分为n
x
×
ny×
nz个小长方体，测区的取值范围为[x
min
,x
min
+n
xdx
]，[y
min
,y
min
+n
ydy
]，[0,n
zdz
]，观测点规则分布于同一高度平面z＝z
l
，其在地面投影位于每个网格的中心，观测点在x方向坐标取值为
观测点在y方向坐标取值为观测点在z方向坐标取值为
[0105]
①
统计a1,b1,c1的值。令a1＝(x-ξ1)，b1＝(y-η1)，c1＝(z-ζ1)，ξ1、η1、ζ1分别表示任一小长方体在x,y,z三方向取值的下限，则
[0106][0107]
由观测点坐标取值和上式，可得到a1,b1,c1的值，用向量a1、b1、 c1表示，将其中的数值从小到大排列，可得
[0108][0109][0110]
c1＝(z-(n
z-1)dz,...,z-dz,z)
t
[0111]
②
依据向量a1、b1、c1的值的个数确定va矩阵的大小。a1中包含2n
x-1个值，b1中包含2n
y-1个值，c1中包含nz个值，则va大小为 nk×
1，且
[0112]
nk＝(2n
x-1)
×
(2n
y-1)
×
nz[0113]
③
将a1、b1、c1中的值依次取出，计算获得va。va中包含所有构成灵敏度矩阵a的值，同时，其对内存要求远远低于容纳灵敏度矩阵a的要求，比如，存储大小为10000*100000的灵敏度矩阵a，存储需要约5.36gb空间，而存储组成该矩阵的值的va，所需要的空间仅为1.6mb。
[0114]
(2)将va代入到反演计算
[0115]
①
建立灵敏度矩阵a中元素和va的关系。对a1中数值进行编号，记为i
x
＝1,
…
,2n
x-1，对b1中数值进行编号有iy＝1,
…
,2n
y-1，对c1中数值进行编号有iz＝1,
…
,nz，对va中数值进行编号有 ik＝1,
…
,(2n
x-1)
×
(2n
y-1)
×
nz。灵敏度矩阵a中的元素为其中i
obs
表示测点编号，jg表示网格编号，对应的观测点坐标为对应的小长方体的在三方向取值下限为由此，计算得到其对应的变量在a1、b1、c1中的位置，进而所得结果在va中的位置和的值。
[0116][0117]
ik＝(i
x-1)
×
(2n
y-1)
×
nz+(i
y-1)
×
nz+iz[0118][0119]
②
利用va，对反演过程中涉及灵敏度矩阵a的计算进行优化。反演计算中涉及灵敏度矩阵a的计算包括两类：一是a与向量相乘，记为av，v表示反演计算中与a相乘的向量，二是a的转置乘a乘向量，记为a
t
av。
[0120]
对于第一种计算，一般是直接进行矩阵计算，但是针对大数据情况下，很容易出现内存不足，无法计算的情况。为了克服这一情况，直接读取va，利用va与a中元素的关系，实时从va中提取元素构建a 中的行向量，与向量v相乘，逐一获得结果向量中的元素。从而能够避免直接读取和操作a，降低了对内存的要求。同时，采用这种算法，会降低运算效率，增加运算时间，由于每一行向量与列向量相乘是相互独立的，所以采用并行计算实现这一过程，从而提高了计算效率。
[0121][0122]
对于第二种计算，利用av计算结果，采用上述同样的思路计算 a
t
m。
[0123]
本发明通过将快速计算灵敏度矩阵的方法引入到重力反演中，降低了重力反演对内存资源的要求，且将并行计算方式引入到灵敏度矩阵与向量相乘的计算中，提高了计算效率，能够很好地适应于大规模数据的重力反演中。
[0124]
本发明可很好地集成在海洋环境监测和探测设备上，实现对海洋区域的重力密度数据分别的探测，为后续判明地下地质结构、识别深部构造特征、寻找隐伏矿床等提供数据基础。
[0125]
如图2所示，本发明还提供一种处理终端100，其包括：
[0126]
存储器101，用于存储程序指令；
[0127]
处理器102，用于运行所述程序指令，以执行所述可适用于大规模数据的重力反演方法的步骤。
[0128]
本说明书所公开的实施例只是对本发明单方面特征的一个例证，本发明的保护范围不限于此实施例，其他任何功能等效的实施例均落入本发明的保护范围内。对于本领域的技术人员来说，可根据以上描述的技术方案以及构思，做出其它各种相应的改变以及变形，而所有的这些改变以及变形都应该属于本发明权利要求的保护范围之内。