基于权相交及多次交叉梯度约束的联合反演方法

1.本发明属于地球物理学技术领域，尤其涉及一种基于权相交及多次交叉梯度约束的联合反演方法。


背景技术：

2.随着地球物理勘探逐渐转向地质条件更为复杂的地区，单一数据类型的反演已无法获得可靠的地球物理模型。单个数据反演势必会产生多解性，并且通常难以获得一致的地质
‑
地球物理模型来全面解释反演的单参数模型结果。传统的联合反演在提升反演结果的分辨率和精度上，效果有限，因此也无法获得更加准确的地下信息。
3.重力异常可以突出深部地质体的地质响应，并且随着深度的增加衰减会变慢。重力梯度异常可以突出浅层地质体的地质响应，并且比重力异常具有更高的水平分辨率。重力和重力梯度联合反演可用于更好指导深部矿产勘探和石油勘探，且可以促进对地球内部物理特性的进一步解释和对矿产分布的理解。
4.在新的经济和社会发展形势下，我国已探明的石油和金属矿产资源已经无法满足经济社会发展的需要，加快地质找矿工作的步伐，缓解矿产资源紧缺的矛盾，已经迫在眉睫。然而，浅表矿产越来越少，找矿难度越来越大，找矿向深部进军是必然趋势。通过国内外找矿实践经验表明，现有矿区深部找矿潜力巨大。同时，隐伏矿床的找矿问题也一直备受关注和重视。
5.因此，由过去的500m以浅勘探深度到现在的从深500
‑
1000m可采深度甚至大于1000m采矿深度，高精度、高分辨率的联合反演方法，对于更好指导深部矿产勘探和石油勘探，促进对地球内部物理特性的进一步解释，至关重要。
6.但是现有的反演结果多解性严重，反演精度和分辨率不够高，对于地下地质体位置的圈定不准确，影响其应用。


技术实现要素：

7.本发明通过重力与多个重力梯度分量权相交并进行多次交叉梯度约束，最大程度挖掘重力与重力梯度数据所包含的地下信息，解决了反演的多解性问题，提高了反演的分辨率与精度，更好指导深部找矿和促进对地球内部物理特性的进一步解释。
8.具体的，本发明是通过以下技术方案实现的：
9.提供一种基于权相交及多次交叉梯度约束的联合反演方法，包括以下步骤：
10.步骤1：依据实际地质情况，剖分地下三维空间，设定两个规则立方体的大小和埋深；
11.步骤2：建立相应的核函数矩阵；
12.步骤3：确定重力及重力梯度异常与物性参数的关系，计算地面观测到的重力及重力梯度异常；
13.步骤4：设定观测点坐标及立方体单元位置，确定核函数矩阵的具体计算公式；
14.步骤5：将步骤4中的计算公式中表示的反演问题写成目标函数的形式；
15.步骤6：确定三维交叉梯度函数，并计算其二范数；
16.步骤7：将步骤6中的三维交叉梯度函数分别添加到步骤5中重力目标函数和重力梯度目标函数，得到联合反演目标函数；
17.步骤8：对步骤7中联合反演目标函数在加权参数域利用共轭梯度法并在交叉梯度函数约束下求取最优解。
18.优选的，步骤1中依据实际地质情况，剖分地下三维空间，设定两个规则立方体的大小和埋深，具体为将地下三维空间剖分为m个大小100m
×
100m
×
100m规则立方体单元，背景密度为0，两个密度差均为1g/cm3的规则立方体大小分别为300m
×
300m
×
200m，地面有n(30
×
30)个观测点。
19.优选的，步骤2中建立相应的核函数矩阵a，维数为n
×
m，n＝30
×
30，m＝30
×
30
×
10。
20.优选的，步骤3中地面观测到的重力及重力梯度异常与物性参数的关系为：
21.am＝d，
22.其中列向量d为地面观测数据，包括重力异常δg，重力梯度异常g
xx
、g
xy
、g
xz
、g
yy
、g
yz
和g
zz
；g
xx
表示重力异常在x
‑
方向的二阶偏导数、g
xy
表示重力异常在x
‑
方向和y
‑
方向的偏导数、g
xz
表示重力异常在x
‑
方向和z
‑
方向的偏导数、g
yy
表示重力异常在y
‑
方向的二阶偏导数、g
yz
表示重力异常在y
‑
方向z
‑
方向的偏导数、g
zz
表示重力异常在z
‑
方向的二阶偏导数，其维数为n；列向量m为物性参数，这里代表密度参数ρ，其维数为m；矩阵a表示连接观测数据d和物性参数m的核函数矩阵，维数为n
×
m，计算重力及重力梯度异常。
23.优选的，步骤4中，观测点坐标p坐标为(x1,y1,z1)，立方体单元q的位置为(α,β,γ)，其中α∈(α1,α2)，β∈(β1,β2)，γ∈(γ1,γ2)；
24.步骤2中，核函数矩阵包括重力异常δg，重力梯度异常g
xx
、g
xy
、g
xz
、g
yy
、g
yz
和g
zz
)的具体计算公式如下：
[0025][0026][0027][0028][0029]
[0030][0031][0032]
其中，万有引力常数g＝6.67
×
10
‑
11
m3·
kg
‑1·
s
‑2，ρ为立方体单元的平均剩余密度，也就上式中中的密度差；x
i
＝x
‑
α
i
，y
j
＝y
‑
β
j
，z
k
＝z
‑
γ
k
；μ
ijk
＝(
‑
1)
i
(
‑
1)
j
(
‑
1)
k
。
[0033]
优选的，步骤5中将式中am＝d表示的反演问题写成目标函数的形式：
[0034]
φ＝||am
‑
d||2＝(am
‑
d)
t
(am
‑
d)，
[0035]
引入正则化约束项后，重力和重力梯度目标函数可分别写成：
[0036][0037][0038]
其中gg1代表重力δg，gg2代表重力梯度g
xx
、g
xy
、g
xz
、g
yy
、g
yz
和g
zz
；w
d
为数据加权矩阵，w
m
为模型加权矩阵，也就是所说的深度权；δ为正则化因子。
[0039]
优选的，步骤6：三维交叉梯度函数可写为：
[0040]
其中：
[0041][0042][0043][0044]
将交叉梯度函数的二范数φ
cg
作为约束条件而不是模添加到反演中：
[0045][0046]
优选的，步骤7：将步骤6中的三维交叉梯度函数分别添加到步骤5中重力目标函数和重力梯度目标函数，得到联合反演目标函数，φ'
gg1
和φ'
gg2
，如下所示：
[0047][0048][0049]
其中，λ
gg1
和λ
gg2
分别表示重力和重力梯度交叉梯度函数的系数。
[0050]
优选的，对步骤7中联合反演目标函数φ'
gg1
和φ'
gg2
在加权参数域利用共轭梯度法并在交叉梯度函数约束下求取最优解m
gg1
和m
gg2
，，将第1次交叉梯度约束得到的两个最优解m
1gg1
和m
1gg2
作为下一次交叉梯度的权分别相交传递到第2次交叉梯度，以此类推，完成方
法的整个过程，得到反演结果m'
gg1
和m'
gg2
。
[0051]
与现有技术相比，本发明具有以下有益的技术效果：
[0052]
本发明本发明通过重力与多个重力梯度分量权相交并进行多次交叉梯度约束，最大程度挖掘重力与重力梯度数据所包含的地下信息，一定程度上解决了反演的多解性问题，提高了反演的分辨率与精度，从而更加准确圈定出地下地质体位置，更好指导深部找矿和促进对地球内部物理特性的进一步解释。
附图说明
[0053]
图1是本发明提供的基于权相交及多次交叉梯度约束的联合反演方法的流程图；
[0054]
图2是本发明中地质体三维位置及观测点、立方体单元位置示例；
[0055]
图3是本发明中重力及重力梯度异常图，其中a是重力δg异常；b是重力梯度g
xx
异常；
[0056]
图4是本发明单独反演结果图，其中黑色虚线框为地质体的真实位置，a是重力δg单独反演结果；b是重力梯度g
xx
单独反演结果；
[0057]
图5是本发明中最终反演结果图，a是重力δg最终反演结果；b是重力梯度g
xx
最终反演结果。
具体实施方式
[0058]
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施例对本发明进行详细描述。需要说明的是，在不冲突的情况下，本发明的实施例及实施例中的特征可以相互组合。
[0059]
除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。
[0060]
具体如图1所示，提供一种基于权相交及多次交叉梯度约束的联合反演方法，包括以下步骤：
[0061]
步骤1：依据实际地质情况，剖分地下三维空间，设定两个规则立方体的大小和埋深；
[0062]
设定两个规则立方体的大小和埋深，具体为将地下三维空间剖分为m(30
×
30
×
10)个大小100m
×
100m
×
100m规则立方体单元，背景密度为0，两个密度差均为1g/cm3的规则立方体大小分别为300m
×
300m
×
200m，地面有n(30
×
30)个观测点，左侧立方体顶面埋深300m，右侧立方体顶面埋深400m，相距700m，具体如图2所示；
[0063]
步骤2：建立相应的核函数矩阵a，维数为n
×
m，n＝30
×
30，m＝30
×
30
×
10；
[0064]
步骤3：确定重力及重力梯度异常与物性参数的关系，计算地面观测到的重力及重力梯度异常，地面观测到的重力及重力梯度异常与物性参数的关系为：
[0065]
am＝d，
[0066]
其中列向量d为地面观测数据，包括重力异常δg，重力梯度异常g
xx
、g
xy
、g
xz
、g
yy
、g
yz
和g
zz
；g
xx
表示重力异常在x
‑
方向的二阶偏导数、g
xy
表示重力异常在x
‑
方向和y
‑
方向的偏导数、g
xz
表示重力异常在x
‑
方向和z
‑
方向的偏导数、g
yy
表示重力异常在y
‑
方向的二阶偏导
数、g
yz
表示重力异常在y
‑
方向z
‑
方向的偏导数、g
zz
表示重力异常在z
‑
方向的二阶偏导数，其维数为n；列向量m为物性参数，这里代表密度参数ρ，其维数为m；矩阵a表示连接观测数据d和物性参数m的核函数矩阵，维数为n
×
m，计算重力及重力梯度异常，具体如图3所示；
[0067]
步骤4：设定观测点坐标及立方体单元位置，确定核函数矩阵的具体计算公式；
[0068]
观测点坐标p坐标为(x1,y1,z1)，立方体单元q的位置为(α,β,γ)，其中α∈(α1,α2)，β∈(β1,β2)，γ∈(γ1,γ2)；
[0069]
步骤2中，核函数矩阵包括重力异常δg，重力梯度异常g
xx
、g
xy
、g
xz
、g
yy
、g
yz
和g
zz
)的具体计算公式如下：
[0070][0071][0072][0073][0074][0075][0076][0077]
其中，万有引力常数g＝6.67
×
10
‑
11
m3·
kg
‑1·
s
‑2，ρ为立方体单元的平均剩余密度，也就上式中中的密度差；x
i
＝x
‑
α
i
，y
j
＝y
‑
β
j
，z
k
＝z
‑
γ
k
；μ
ijk
＝(
‑
1)
i
(
‑
1)
j
(
‑
1)
k
。
[0078]
步骤5：将步骤4中的计算公式中表示的反演问题写成目标函数的形式；
[0079]
将式中am＝d表示的反演问题写成目标函数的形式：
[0080]
φ＝||am
‑
d||2＝(am
‑
d)
t
(am
‑
d)，
[0081]
引入正则化约束项后，重力和重力梯度目标函数可分别写成：
[0082][0083][0084]
其中gg1代表重力δg，gg2代表重力梯度g
xx
、g
xy
、g
xz
、g
yy
、g
yz
和g
zz
；w
d
为数据加权矩阵，w
m
为模型加权矩阵，也就是所说的深度权；δ为正则化因子。
[0085]
步骤6：确定三维交叉梯度函数，并计算其二范数；
[0086]
三维交叉梯度函数可写为：
[0087]
其中：
[0088][0089][0090][0091]
将交叉梯度函数的二范数φ
cg
作为约束条件而不是模添加到反演中：
[0092][0093]
步骤7：将步骤6中的三维交叉梯度函数分别添加到步骤5中重力目标函数和重力梯度目标函数，得到联合反演目标函数，φ'
gg1
和φ'
gg2
，如下所示：
[0094][0095][0096]
其中，λ
gg1
和λ
gg2
分别表示重力和重力梯度交叉梯度函数的系数。
[0097]
步骤8：对步骤7中联合反演目标函数在加权参数域利用共轭梯度法并在交叉梯度函数约束下求取最优解，具体如下：
[0098][0099]
将第1次交叉梯度约束得到的两个最优解m
1gg1
和m
1gg2
作为下一次交叉梯度的权分别相交传递到第2次交叉梯度，以此类推，完成方法的整个过程，得到反演结果m'
gg1
和m'
gg2
。
[0100]
结果如图4中单独反演结果图，其中黑色虚线框为地质体的真实位置，a是重力δg单独反演结果，b是重力梯度g
xx
单独反演结果；图5中a是重力δg最终反演结果，b是重力梯度g
xx
最终反演结果。
[0101]
对比例：单独反演及最终反演结果与原始重力、重力梯度异常的均方根误差对比，
[0102] rmse(重力)rmse(重力梯度)单独反演1.42520.0131最终反演0.05220.0079
[0103]
均方根值越大，反演结果越差，精度和分辨率越低，如下为均方根误差(重力：rmse'；重力梯度：rmse”)计算公式：
[0104]
[0105][0106]
a'代表原始重力异常(δg)的核函数矩阵，a”代表原始重力梯度异常(g
xx
)的核函数矩阵；m'代表单独反演及最终反演得到的密度(重力)，m”代表单独反演及最终反演得到的密度(重力梯度)；d'代表原始重力异常(δg)，d”代表原始重力梯度异常(g
xx
)。
[0107]
最后应说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或等同替换，而不脱离本发明技术方案的精神和范围。