CN117688785B - 一种基于种植思想的全张量重力梯度数据反演方法 - Google Patents
一种基于种植思想的全张量重力梯度数据反演方法 Download PDFInfo
- Publication number
- CN117688785B CN117688785B CN202410146681.6A CN202410146681A CN117688785B CN 117688785 B CN117688785 B CN 117688785B CN 202410146681 A CN202410146681 A CN 202410146681A CN 117688785 B CN117688785 B CN 117688785B
- Authority
- CN
- China
- Prior art keywords
- data
- cuboid
- inversion
- adjacent
- elements
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000005484 gravity Effects 0.000 title claims abstract description 98
- 238000000034 method Methods 0.000 title claims abstract description 45
- 230000012010 growth Effects 0.000 claims abstract description 107
- 239000011159 matrix material Substances 0.000 claims abstract description 44
- 230000035945 sensitivity Effects 0.000 claims abstract description 36
- 230000036433 growing body Effects 0.000 claims abstract description 24
- 238000004364 calculation method Methods 0.000 claims description 9
- 238000010348 incorporation Methods 0.000 claims description 5
- 230000009467 reduction Effects 0.000 claims description 4
- 230000008859 change Effects 0.000 claims description 3
- 238000012360 testing method Methods 0.000 description 5
- 238000004590 computer program Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000003384 imaging method Methods 0.000 description 4
- 230000000877 morphologic effect Effects 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 150000003839 salts Chemical class 0.000 description 3
- 230000002159 abnormal effect Effects 0.000 description 2
- 230000005587 bubbling Effects 0.000 description 2
- 239000013078 crystal Substances 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000013499 data model Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明属于地球物理数据解译技术领域，涉及一种基于种植思想的全张量重力梯度数据反演方法，包括：根据长方体元，计算灵敏度矩阵；将剩余密度非零的长方体元作为生长体；搜索生长体的相邻长方体元；结合灵敏度矩阵计算每个相邻长方体元纳入生长体后的密度模型的目标函数值；在相邻长方体元纳入生长体后的密度模型满足预设生长条件的情况下，选择目标函数值最小的密度模型作为新的密度模型继续迭代；否则停止迭代得到最终反演结果；利用单个分量数据分别迭代生长，获得各个分量数据反演结果；根据各个分量数据的反演结果定权，将加权全张量重力梯度数据联合反演，获得最终联合反演结果。其有益效果是，降低反演问题的多解性、提高反演分辨率。
Description
技术领域
本发明涉及地球物理数据解译技术领域，尤其涉及一种基于种植思想的全张量重力梯度数据反演方法。
背景技术
重力梯度数据反演能够获取地下剩余密度结构，剩余密度为地层中异常体的密度与地层密度之间的差异。因而反演能够为深部找矿提供重要依据。重力梯度数据，包含Vxx，Vxy，Vxz，Vyy，Vyz，Vzz数据，表示重力位V在X、Y、Z轴方向上的二阶导数。重力梯度数据具有高分辨率，高频信息多的特点。全张量重力梯度数据表示六种重力分量数据。一些学者利用重力梯度数据进行目标矿体的反演并取得了高分辨率的解译结果。重力梯度数据的三维物性反演在数学上是病态问题，其问题具有固有的多解性，所以要通过引入先验约束条件使该问题摆脱病态。目前，大多数的三维物性反演方法都属于传统的正则化反演方法。通过不同分量重力梯度数据的联合，也可以降低多解性。
除了正则化反演，还存在一类依靠检索思想的反演方法，即通过在可能的解的空间中搜寻最优解。例如，冒泡法与O-R-F（open-reject-filled principal）方法都采用重力异常数据来进行二维反演。但是冒泡法与O-R-F方法都存在着无法同时对不同密度异常体进行反演的问题。而使用源体生长方法从初始的零剩余密度分布估计地下三维的剩余密度分布，能够实现三维的重力异常数据的反演，但反演结果往往呈现出一种边界离散的结果，对于反演结果的内部也无法避免空洞的产生。而基于引入系统搜索最优解的思想提出种植反演的方法，能够实现不同密度下的三维反演，并利用多分量重力梯度数据进行三维反演。但该方法存在着深部反演效果不佳，冗余生长的问题。因此，反演问题的多解性有待进一步降低，反演分辨率也有待进一步提高。
发明内容
要解决的技术问题
鉴于现有技术的上述缺点、不足，本发明提供一种基于种植思想的全张量重力梯度数据反演方法，其解决了如何降低反演问题的多解性、提高反演分辨率的技术问题。
技术方案
为了达到上述目的，本发明采用的主要技术方案包括：
第一方面，本发明提供一种基于种植思想的全张量重力梯度数据反演方法，包括：
设置地下空间均匀剖分网格的数量及尺寸以建立笛卡尔坐标系，得到长方体元；
根据长方体元，计算灵敏度矩阵；
根据先验地质信息设定长方体元的初始剩余密度以得到初始密度模型，将剩余密度非零的长方体元作为生长体；
搜索生长体中长方体元上下、左右和前后的六个相邻长方体元；
根据相邻的长方体元结合灵敏度矩阵计算每个相邻长方体元纳入生长体后的密度模型的目标函数值；
在相邻长方体元纳入生长体后的密度模型满足预设生长条件的情况下，选择目标函数值最小的密度模型作为新的密度模型，继续迭代；
在相邻长方体元纳入生长体后的密度模型不满足预设生长条件的情况下，停止迭代，得到最终反演结果；
利用单个分量数据分别进行迭代生长，获得各个分量数据的反演结果；
根据各个分量数据的反演结果进行定权，将加权的全张量重力梯度数据进行联合反演，获得最终的联合反演结果。
本发明提供的一种基于种植思想的全张量重力梯度数据反演方法，对地下空间剖分，建立笛卡尔坐标系，令原点为O，XOY面为地表面，Z轴垂直向下为正方向。设置地下空间均匀剖分网格的数量及尺寸，每个网格即为长方体元。根据剖分的长方体元，计算灵敏度矩阵，其由单位剩余密度值下的每个长方体元对每个观测点产生的梯度数据组成。根据先验地质信息给定地下空间内剖分长方体元的剩余密度，并在地下空间内设置“种子”，即给定的剩余密度非零的长方体元。种子的数量可以不唯一，剩余密度值也可不相同。地下剖分长方体元的所有剩余密度值构成了密度模型。随着迭代生长的进行，种子会以类似于晶体生长的方式进行增长，称为生长体；将生长体中的长方体元上下左右前后的六个长方体元为相邻状态，迭代伴随着相邻长方体元转化为生长体的过程，其生长体进行了增长，也产生了新的相邻长方体元。根据各个单分量数据反演结果，通过对各梯度分量数据设置权重，使生长体的生长更倾向依赖于具有高权重的梯度分量，进一步提高反演结果的分辨率。改进生长判断准则，实现对不同深度长方体元生长判断的平衡，从而提高了深部目标的成像分辨率。
可选地，所述根据长方体元，计算灵敏度矩阵，包括：
根据如下公式计算每个长方体元产生的梯度分量数据：
；
其中，为引力场在x方向的二阶导数、/>为引力场在y方向的二阶导数、/>为引力场在z方向的二阶导数、/>为引力场在x方向与y方向的二阶导数、/>为引力场在x方向与z方向的二阶导数、/>为引力场在y方向与z方向的二阶导数，上述的六种数据合称为重力全张量梯度数据；/>为计算梯度数据的中间变量、/>为计算梯度数据的中间变量，γ表示万有引力常数；(x，y，z)表示观测点坐标；(ξi，ηj，ζk)表示剖分长方体元的对顶点坐标；ρ表示剩余密度；
在ρ=1g/cm3的情况下，将长方体元对观测点的重力梯度数据构成的矩阵作为灵敏度矩阵。
通过上述公式能够对观测点的重力梯度数据进行计算，进而组成灵敏度矩阵，以供后续预测全张亮重力梯度数据使用。
可选地，所述根据相邻的长方体元结合灵敏度矩阵计算每个相邻长方体元纳入生长体后的密度模型的目标函数值，包括：
根据当前迭代中的密度模型，结合灵敏度矩阵计算每个相邻长方体元分别纳入生长体后，对应的密度模型产生的预测全张量重力梯度数据；
根据当前迭代中的每个密度模型的预测全张量重力梯度数据，计算各个重力梯度数据的数据残差并联合；
根据相邻剖分长方体元，计算其分别纳入生长体后的密度模型的正则化函数值；
根据联合残差与正则化函数，计算不同相邻剖分长方体元纳入生长体后的不同密度模型的目标函数值。
正则化函数本质上体现为对地下剖分模型元素的先验性矩阵。正则化函数的选择能够一定程度上影响反演结果的形态分布。根据当前的种子距离的约束函数作为正则化函数，并通过正则化因子平衡数据不匹配函数与正则化函数。在此约束下，目标会倾向于产生一个内部没有空洞的实体。并能在深度上能够产生适当的约束。
可选地，所述根据当前迭代中的密度模型，结合灵敏度矩阵计算每个相邻长方体元分别纳入生长体后，对应的密度模型产生的预测全张量重力梯度数据，包括：
根据如下公式计算相邻剖分长方体元纳入生长体后的密度模型产生的预测全张量重力梯度数据：
；其中，dpre表示预测全张量重力梯度数据；A表示当前梯度分量数据对应的灵敏度矩阵；mpre表示将本次相邻的一个长方体元纳入原本的生长体后的剩余密度模型分量。
结合灵敏度矩阵和本次相邻的一个长方体元纳入原本的生长体后的剩余密度模型分量对全张量重力梯度数据进行预测，以便于根据预测数据计算数据残差。
可选地，所述根据当前迭代中的每个密度模型的预测全张量重力梯度数据，计算各个重力梯度数据的数据残差并联合，包括：
根据如下公式计算全张量重力梯度数据中各梯度分量数据残差：
；其中，φs(mpre)表示第s个梯度分量数据残差，dobs表示观测数据；
根据如下公式得到联合残差数据：
；其中，φsum(mpre)表示全张量数据的联合残差数据；ps为每种梯度分量残差的权重。
梯度数据的每个分量都包含不同的地质信息，分量太少，则信息量不全面，所以多分量的数据联合反演能够提高反演效果。通过将各个重力梯度分量数据结果进行累加求和，以此平衡各个分量数据对反演结果的影响。
可选地，所述根据相邻剖分长方体元，计算其分别纳入生长体后的密度模型的正则化函数值，包括：
根据如下公式计算每次迭代中根据相邻长方体元预测的各个密度模型的正则化函数值：
；其中，u表示选择的长方体元，U表示从本次地下的剩余密度mpre中选择的生长体的长方体元的集合；θ(mpre)表示与密度模型mpre对应的正则化函数值，mu表示第u个相邻的长方体元的剩余密度值，/>表示避免奇异值的极小值，范围为0~10-15，lu表示第u个相邻的长方体元到生长体中心坐标间的距离。
在迭代过程中，每次迭代要结合正则化函数计算相邻元素的目标函数，选择最小目标函数值进行增长，直到得到最终反演结果。
可选地，所述根据联合残差与正则化函数，计算不同相邻剖分长方体元纳入生长体后的不同密度模型的目标函数值，包括：
根据如下公式计算密度模型长方体元的目标函数值：
；
其中，Г(mpre)表示相邻的剖分长方体元被纳入当前迭代的生长体后，新的密度模型产生的目标函数值，θ(mpre)表示每个相邻长方体的增长产生新的密度模型该密度模型的正则化函数值；μ表示正则化参数。
公式中的正则化函数可以根据不同的选择，使模型在数据拟合的前提下体现出不同的形态特征。
可选地，所述相邻长方体元纳入生长体后的密度模型满足预设生长条件，包括：
相邻长方体元纳入生长体后的密度模型使得联合残差值的降低以及残差改变量满足预设阈值的要求。
可选地，根据下式确定相邻长方体元纳入生长体后的密度模型使得联合残差值的降低：
根据下式确定相邻长方体元纳入生长体后的密度模型使得残差改变量满足预设阈值的要求：
；
其中，表示相邻的长方体元纳入生长体后的密度模型的预测的联合残差值，表示当前的密度模型产生的联合残差值，zt表示第t个相邻的长方体元的Z轴方向上的深度，z0表示观测点到地面的垂直距离，α表示衰减因子，α=3/2，δ表示设置的阈值。
阈值的大小影响反演结果的好坏，较低的阈值往往会造成多余的生长，造成异常体形状上的变形。较高的阈值在数据上的拟合又不甚理想，一个合适的阈值，能够改善反演结果。
可选地，所述根据各个分量数据的反演结果进行定权，将加权的全张量重力梯度数据进行联合反演，包括：
在当前参与反演计算的梯度分量数据残差的权重为1，其余梯度分量数据残差的权重pj为0的情况下进行反演，获得每个梯度分量数据的反演结果：，其中，/>为当前反演结果下的各梯度分量数据残差的和，/>为当前反演结果下第S个梯度分量数据的残差；
根据获得的各个梯度分量的反演结果计算各反演结果产生的方差：；
根据如下公式对每个梯度分量数据单独反演的方差进行定权：
；
将权重应用于全张量重力梯度数据残差求取中，利用全张量重力梯度数据进行联合反演获得最终的联合反演结果；
其中，σ2表示单个梯度分量数据反演结果的方差，表示单个梯度分量数据第r个观测数据，/>表示单个梯度分量数据反演结果产生的第r个重力梯度分量数据，R表示测点总数，R的范围为0~108。
通过对各梯度分量数据设置权重，使生长体的生长更倾向依赖于具有高权重的梯度分量，进一步提高反演结果的分辨率。
有益效果
本发明的有益效果是：本发明的一种基于种植思想的全张量重力梯度数据反演方法，根据各个单分量数据反演结果，对各梯度分量数据设置权重，使生长体的生长更倾向依赖于具有高权重的梯度分量，进一步提高反演结果的分辨率。相较于相关技术，本发明改进生长判断准则，实现对不同深度长方体元生长判断的平衡，从而提高了深部目标的成像分辨率。
附图说明
图1为本发明实施例提供的基于种植思想的全张量重力梯度数据反演方法的流程示意图；
图2为本发明实施例提供的理论模型示意图；
图3为本发明实施例提供的实测数据模型示意图；
图4为本发明实施例提供的不同密度的立方体模型反演结果图；
图5为本发明实施例提供的实测数据的x轴剖面反演结果；
图6为本发明实施例提供的实测数据的y轴剖面反演结果；
图7为本发明实施例提供的实测数据z轴剖面反演结果。
具体实施方式
为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。虽然附图中显示了本发明的示例性实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更清楚、透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。
第一方面，参照图1，本实施例提供了一种基于种植思想的全张量重力梯度数据反演方法，包括：
S101，设置地下空间均匀剖分网格的数量及尺寸以建立笛卡尔坐标系，得到长方体元。
S102，根据长方体元，计算灵敏度矩阵。
S103，根据先验地质信息设定长方体元的初始剩余密度以得到初始密度模型，将剩余密度非零的长方体元作为生长体。
S104，搜索生长体中长方体元上下、左右和前后的六个相邻长方体元。
S105，根据相邻的长方体元结合灵敏度矩阵计算每个相邻长方体元纳入生长体后的密度模型的目标函数值。
S106，在相邻长方体元纳入生长体后的密度模型满足预设生长条件的情况下，选择目标函数值最小的密度模型作为新的密度模型，继续迭代。
S107，在相邻长方体元纳入生长体后的密度模型不满足预设生长条件的情况下，停止迭代，得到最终反演结果。
S108，利用单个分量数据分别进行迭代生长，获得各个分量数据的反演结果。
S109，根据各个分量数据的反演结果进行定权，将加权的全张量重力梯度数据进行联合反演，获得最终的联合反演结果。
本发明实施例提供的一种基于种植思想的全张量重力梯度数据反演方法，对地下空间剖分，建立笛卡尔坐标系，令原点为O，XOY面为地表面，Z轴垂直向下为正方向。设置地下空间均匀剖分网格的数量及尺寸，每个网格即为长方体元。根据剖分的长方体元，计算灵敏度矩阵，其由单位剩余密度值下的每个长方体元对每个观测点产生的梯度数据组成。根据先验地质信息给定地下空间内剖分长方体元的剩余密度，并在地下空间内设置“种子”，即给定的剩余密度非零的长方体元。种子的数量可以不唯一，剩余密度值也可不相同。地下剖分长方体元的所有剩余密度值构成了密度模型。随着迭代生长的进行，种子会以类似于晶体生长的方式进行增长，称为生长体；将生长体中的长方体元上下左右前后的六个长方体元为相邻状态，迭代伴随着相邻长方体元转化为生长体的过程，其生长体进行了增长，也产生了新的相邻长方体元。根据各个单分量数据反演结果，通过对各梯度分量数据设置权重，使生长体的生长更倾向依赖于具有高权重的梯度分量，进一步提高反演结果的分辨率。改进生长判断准则，实现对不同深度长方体元生长判断的平衡，从而提高了深部目标的成像分辨率。
可选地，所述根据长方体元，计算灵敏度矩阵，包括：
根据如下公式计算每个长方体元产生的梯度分量数据：
；
其中，其中，为引力场在x方向的二阶导数、/>为引力场在y方向的二阶导数、为引力场在z方向的二阶导数、/>为引力场在x方向与y方向的二阶导数、/>为引力场在x方向与z方向的二阶导数、/>为引力场在y方向与z方向的二阶导数，上述的六种数据合称为重力全张量梯度数据；/>为计算梯度数据的中间变量、/>为计算梯度数据的中间变量，γ表示万有引力常数；(x，y，z)表示观测点坐标；(ξi，ηj，ζk)表示剖分长方体元的对顶点坐标；ρ表示剩余密度；
在ρ=1g/cm3的情况下，将长方体元对观测点的重力梯度数据构成的矩阵作为灵敏度矩阵。
通过上述公式能够对观测点的重力梯度数据进行计算，进而组成灵敏度矩阵，以供后续预测全张亮重力梯度数据使用。
可选地，所述根据相邻的长方体元结合灵敏度矩阵计算每个相邻长方体元纳入生长体后的密度模型的目标函数值，包括：
根据当前迭代中的密度模型，结合灵敏度矩阵计算每个相邻长方体元分别纳入生长体后，对应的密度模型产生的预测全张量重力梯度数据；
根据当前迭代中的每个密度模型的预测全张量重力梯度数据，计算各个重力梯度数据的数据残差并联合；
根据相邻剖分长方体元，计算其分别纳入生长体后的密度模型的正则化函数值；
根据联合残差与正则化函数，计算不同相邻剖分长方体元纳入生长体后的不同密度模型的目标函数值。
正则化函数本质上体现为对地下剖分模型元素的先验性矩阵。正则化函数的选择能够一定程度上影响反演结果的形态分布。根据当前的种子距离的约束函数作为正则化函数，并通过正则化因子平衡数据不匹配函数与正则化函数。在此约束下，目标会倾向于产生一个内部没有空洞的实体。并能在深度上能够产生适当的约束。
可选地，所述根据当前迭代中的密度模型，结合灵敏度矩阵计算每个相邻长方体元分别纳入生长体后，对应的密度模型产生的预测全张量重力梯度数据，包括：
根据如下公式计算相邻剖分长方体元纳入生长体后的密度模型产生的预测全张量重力梯度数据：
其中，dpre表示预测全张量重力梯度数据；A表示当前梯度分量数据对应的灵敏度矩阵；mpre表示将本次相邻的一个长方体元纳入原本的生长体后的剩余密度模型分量。
结合灵敏度矩阵和本次相邻的一个长方体元纳入原本的生长体后的剩余密度模型分量对全张量重力梯度数据进行预测，以便于根据预测数据计算数据残差。
可选地，所述根据当前迭代中的每个密度模型的预测全张量重力梯度数据，计算各个重力梯度数据的数据残差并联合，包括：
根据如下公式计算全张量重力梯度数据中各梯度分量数据残差：
其中，φs(mpre)表示第s个梯度分量数据残差，dobs表示观测数据；
根据如下公式得到联合残差数据：
其中，φsum(mpre)表示全张量数据的联合残差数据；ps为每种梯度分量残差的权重。
梯度数据的每个分量都包含不同的地质信息，分量太少，则信息量不全面，所以多分量的数据联合反演能够提高反演效果。通过将各个重力梯度分量数据结果进行累加求和，以此平衡各个分量数据对反演结果的影响。
可选地，所述根据相邻剖分长方体元，计算其分别纳入生长体后的密度模型的正则化函数值，包括：
根据如下公式计算每次迭代中根据相邻长方体元预测的各个密度模型的正则化函数值：
其中，u表示选择的长方体元，U表示从本次地下的剩余密度mpre中选择的生长体的长方体元的集合；θ(mpre)表示与密度模型mpre对应的正则化函数值，mu表示第u个相邻的长方体元的剩余密度值，表示避免奇异值的极小值，范围为0~10-15，lu表示第u个相邻的长方体元到生长体中心坐标间的距离。
在迭代过程中，每次迭代要结合正则化函数计算相邻元素的目标函数，选择最小目标函数值进行增长，直到得到最终反演结果。
可选地，所述根据联合残差与正则化函数，计算不同相邻剖分长方体元纳入生长体后的不同密度模型的目标函数值，包括：
根据如下公式计算密度模型长方体元的目标函数值：
其中，Г(mpre)表示相邻的剖分长方体元被纳入当前迭代的生长体后，新的密度模型产生的目标函数值，θ(mpre)表示每个相邻长方体的增长产生新的密度模型该密度模型的正则化函数值；μ表示正则化参数。
公式中的正则化函数可以根据不同的选择，使模型在数据拟合的前提下体现出不同的形态特征。
可选地，所述相邻长方体元纳入生长体后的密度模型满足预设生长条件，包括：
相邻长方体元纳入生长体后的密度模型使得联合残差值的降低以及残差改变量满足预设阈值的要求。
可选地，根据下式确定相邻长方体元纳入生长体后的密度模型使得联合残差值的降低：
根据下式确定相邻长方体元纳入生长体后的密度模型使得残差改变量满足预设阈值的要求：
其中，表示相邻的长方体元纳入生长体后的密度模型的预测的联合残差值，表示当前的密度模型产生的联合残差值，zt表示第t个相邻的长方体元的Z轴方向上的深度，z0表示观测点到地面的垂直距离，α表示衰减因子，α=3/2，δ表示设置的阈值。
阈值的大小影响反演结果的好坏，较低的阈值往往会造成多余的生长，造成异常体形状上的变形。较高的阈值在数据上的拟合又不甚理想，一个合适的阈值，能够改善反演结果。
可选地，所述根据各个分量数据的反演结果进行定权，将加权的全张量重力梯度数据进行联合反演，包括：
在当前参与反演计算的梯度分量数据残差的权重为1，其余梯度分量数据残差的权重pj为0的情况下进行反演，获得每个梯度分量数据的反演结果：，其中，/>为当前反演结果下的各梯度分量数据残差的和，/>为当前反演结果下第S个梯度分量数据的残差；
根据获得的各个梯度分量的反演结果计算各反演结果产生的方差：；
根据如下公式对每个梯度分量数据单独反演的方差进行定权：
将权重应用于全张量重力梯度数据残差求取中，利用全张量重力梯度数据进行联合反演获得最终的联合反演结果；
其中，σ2表示单个梯度分量数据反演结果的方差，表示单个梯度分量数据第r个观测数据，/>表示单个梯度分量数据反演结果产生的第r个重力梯度分量数据，R表示测点总数，R的范围为0~108。
通过对各梯度分量数据设置权重，使生长体的生长更倾向依赖于具有高权重的梯度分量，进一步提高反演结果的分辨率。
下面结合一具体实施例来进一步介绍本申请提出的一种基于种植思想的全张量重力梯度数据反演方法：
本实施例中，根据理论模型数据与实测数据进行数据试验。
理论模型为一对剩余密度、大小不同的立方体模型，模型一的剩余密度为0.5g/cm3，模型二的剩余密度为-0.5g/cm3，模型形状如图2所示。其中，对其产生的重力梯度数据加入5%的高斯随机噪声。
实测数据为美国文顿盐丘地质区域梯度数据，该数据由美国Bell Geospace公司提供。该数据为一个10000×6的矩阵，矩阵的每一行代表对应质量元的重力梯度数据的6个分量。本次实例中对其数据进行抽稀，选择800×6的矩阵进行反演。模型形状如图3所示。
步骤1：对地下空间进行空间剖分，建立笛卡尔坐标系。设置剖分长方体元的尺寸，数量等。
本实施例中，对理论模型试验的地下空间进行网格剖分，将空间划分为20×20×20的长方体元，其单个小长方体元的x，y，z方向间隔均为50m。
对实测数据试验的地下空间进行网格剖分，将空间划分为20×20×20的长方体元；网格元素中z轴方向间隔为50m，x、y轴方向间隔200m。
步骤2：根据地下模型中剖分的长方体元，计算灵敏度矩阵。
本实施例中，根据地下模型中剖分的长方体元，利用梯度分量数据计算公式计算灵敏度矩阵。
步骤3：根据先验的地质信息，给定剖分长方体元的初始剩余密度。包括对种子（生长体）的设定。
本实施例中，如图2、图3所示，根据先验的地质信息分别确定了两个试验的初始模型的位置，密度。
步骤4：根据当前的生长体，找到相邻的长方体元，其中初次迭代生长中的生长体为设置的种子。
本实例中，根据设置的生长体，将生长体中的长方体元上下左右前后的六个长方体元为相邻状态。
步骤5：根据正则化函数计算每个相邻剖分长方体元纳入生长体后的密度模型的目标函数值。
步骤5.1：根据公式计算相邻剖分长方体元纳入生长体后的密度模型的全张量重力梯度数据。
步骤5.2：根据步骤5.1所求的预测数据，计算全张量重力梯度数据中各梯度分量数据残差并联合。
步骤5.3：根据正则化函数计算公式计算相邻剖分长方体元纳入生长体后的密度模型的正则化函数值。
步骤5.4：利用所求的数据残差与正则化函数值，根据目标函数计算公式计算相邻剖分长方体元纳入生长体后的密度模型的目标函数值。
步骤6：判断每个相邻长方体是否满足生长判断条件，如果存在满足生长条件的相邻长方体，选择产生最小目标函数值的长方体元纳入生长体，如果不存在，停止迭代生长，获得反演结果。
步骤7：将每个分量数据进行上述的流程后，获得反演结果。
步骤8：根据反演结果计算残差并定权，决定全张量重力梯度数据中每个分量的权重，并进行上述的迭代生长，以获得最终的反演结果。
本实例中，终止判断中的两个阈值均设置为：0.04。模型试验要经历124次迭代，文顿盐丘实测数据试验要经历90次迭代。
模型结果如图4所示，其中的边框为模型的真实位置。从模型中可以看出本方法能够同时反演不同密度的异常体，且具有一定的抗噪性。
文顿盐丘的实测数据结果如图5、图6、图7所示，选择y =3334.44km，x=442.46km和z=325m三个剖面展示反演结果。得到的成像结果图能够显示地质体的空间位置与形态。从图像可以看出，图像的深部方向上边缘是清晰的，并具有明显的形态特征。这证明了该方法能够高效地反演数据，并且反演结果的边界特点好，形态还原度高的特点。
第二方面，本发明实施例提供一种计算机可读存储介质，其上存储有计算机程序，所述程序执行时实现上述第一方面中任一项所述基于种植思想的全张量重力梯度数据反演方法。
第三方面，本发明实施例提供一种存储设备，包括存储介质和处理器，所述存储介质存储有计算机程序，所述程序被处理器执行时实现上述第一方面中任一项所述基于种植思想的全张量重力梯度数据反演方法。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例，或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质（包括但不限于磁盘存储器、CD-ROM、光学存储器等）上实施的计算机程序产品的形式。
显然，本领域的技术人员可以对本发明进行各种修改和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也应该包含这些修改和变型在内。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行改动、修改、替换和变型。
Claims (10)
1.一种基于种植思想的全张量重力梯度数据反演方法，其特征在于，包括：
设置地下空间均匀剖分网格的数量及尺寸以建立笛卡尔坐标系，得到长方体元；
根据长方体元，计算灵敏度矩阵；
根据先验地质信息设定长方体元的初始剩余密度以得到初始密度模型，将剩余密度非零的长方体元作为生长体；
搜索生长体中长方体元上下、左右和前后的六个相邻长方体元；
根据相邻的长方体元结合灵敏度矩阵计算每个相邻长方体元纳入生长体后的密度模型的目标函数值；
在相邻长方体元纳入生长体后的密度模型满足预设生长条件的情况下，选择目标函数值最小的密度模型作为新的密度模型，继续迭代；
在相邻长方体元纳入生长体后的密度模型不满足预设生长条件的情况下，停止迭代，得到最终反演结果；
利用单个分量数据分别进行迭代生长，获得各个分量数据的反演结果；
根据各个分量数据的反演结果进行定权，将加权的全张量重力梯度数据进行联合反演，获得最终的联合反演结果。
2.根据权利要求1所述的基于种植思想的全张量重力梯度数据反演方法，其特征在于，所述根据长方体元，计算灵敏度矩阵，包括：
根据如下公式计算每个长方体元产生的梯度分量数据：
；
其中，为引力场在x方向的二阶导数、/>为引力场在y方向的二阶导数、/>为引力场在z方向的二阶导数、/>为引力场在x方向与y方向的二阶导数、/>为引力场在x方向与z方向的二阶导数、/>为引力场在y方向与z方向的二阶导数，上述的六种数据合称为重力全张量梯度数据；/>为计算梯度数据的中间变量、/>为计算梯度数据的中间变量，γ表示万有引力常数；(x，y，z)表示观测点坐标；(ξ i ，η j ，ζ k )表示剖分长方体元的对顶点坐标；ρ表示剩余密度；
在ρ=1g/cm3的情况下，将长方体元对观测点的重力梯度数据构成的矩阵作为灵敏度矩阵。
3.根据权利要求1所述的基于种植思想的全张量重力梯度数据反演方法，其特征在于，所述根据相邻的长方体元结合灵敏度矩阵计算每个相邻长方体元纳入生长体后的密度模型的目标函数值，包括：
根据当前迭代中的密度模型，结合灵敏度矩阵计算每个相邻长方体元分别纳入生长体后，对应的密度模型产生的预测全张量重力梯度数据；
根据当前迭代中的每个密度模型的预测全张量重力梯度数据，计算各个重力梯度数据的数据残差并联合；
根据相邻剖分长方体元，计算其分别纳入生长体后的密度模型的正则化函数值；
根据联合残差与正则化函数，计算不同相邻剖分长方体元纳入生长体后的不同密度模型的目标函数值。
4.根据权利要求3所述的基于种植思想的全张量重力梯度数据反演方法，其特征在于，所述根据当前迭代中的密度模型，结合灵敏度矩阵计算每个相邻长方体元分别纳入生长体后，对应的密度模型产生的预测全张量重力梯度数据，包括：
根据如下公式计算相邻剖分长方体元纳入生长体后的密度模型产生的预测全张量重力梯度数据：；
其中，d pre 表示预测全张量重力梯度数据；A表示当前梯度分量数据对应的灵敏度矩阵；m pre 表示将本次相邻的一个长方体元纳入原本的生长体后的剩余密度模型分量。
5.根据权利要求4所述的基于种植思想的全张量重力梯度数据反演方法，其特征在于，所述根据当前迭代中的每个密度模型的预测全张量重力梯度数据，计算各个重力梯度数据的数据残差并联合，包括：
根据如下公式计算全张量重力梯度数据中各梯度分量数据残差：；
其中，φ s (m pre )表示第s个梯度分量数据残差，d obs 表示观测数据；
根据如下公式得到联合残差数据：；
其中，φ sum (m pre )表示全张量数据的联合残差数据；p s 为每种梯度分量残差的权重。
6.根据权利要求4所述的基于种植思想的全张量重力梯度数据反演方法，其特征在于，所述根据相邻剖分长方体元，计算其分别纳入生长体后的密度模型的正则化函数值，包括：
根据如下公式计算每次迭代中根据相邻长方体元预测的各个密度模型的正则化函数值：；
其中，u表示选择的长方体元，U表示从本次地下的剩余密度m pre 中选择的生长体的长方体元的集合；θ(m pre )表示与密度模型m pre 对应的正则化函数值，m u 表示第u个相邻的长方体元的剩余密度值，表示避免奇异值的极小值，范围为0~10-15，l u 表示第u个相邻的长方体元到生长体中心坐标间的距离。
7.根据权利要求5所述的基于种植思想的全张量重力梯度数据反演方法，其特征在于，所述根据联合残差与正则化函数，计算不同相邻剖分长方体元纳入生长体后的不同密度模型的目标函数值，包括：
根据如下公式计算密度模型长方体元的目标函数值：
；
其中，Г(m pre )表示相邻的剖分长方体元被纳入当前迭代的生长体后，新的密度模型产生的目标函数值，θ(m pre )表示每个相邻长方体的增长产生新的密度模型该密度模型的正则化函数值；μ表示正则化参数。
8.根据权利要求1所述的基于种植思想的全张量重力梯度数据反演方法，其特征在于，所述相邻长方体元纳入生长体后的密度模型满足预设生长条件，包括：
相邻长方体元纳入生长体后的密度模型使得联合残差值的降低以及残差改变量满足预设阈值的要求。
9.根据权利要求8所述的基于种植思想的全张量重力梯度数据反演方法，其特征在于，
根据下式确定相邻长方体元纳入生长体后的密度模型使得联合残差值的降低：；
根据下式确定相邻长方体元纳入生长体后的密度模型使得残差改变量满足预设阈值的要求：；
其中，表示相邻的长方体元纳入生长体后的密度模型的预测的联合残差值，/>表示当前的密度模型产生的联合残差值，z t 表示第t个相邻的长方体元的Z轴方向上的深度，z 0 表示观测点到地面的垂直距离，α表示衰减因子，α=3/2，δ表示设置的阈值。
10.根据权利要求1所述的基于种植思想的全张量重力梯度数据反演方法，其特征在于，所述根据各个分量数据的反演结果进行定权，将加权的全张量重力梯度数据进行联合反演，包括：
在当前参与反演计算的梯度分量数据残差的权重为1，其余梯度分量数据残差的权重p j 为0的情况下进行反演，获得每个梯度分量数据的反演结果：，其中，/>为当前反演结果下的各梯度分量数据残差的和，/>为当前反演结果下第S个梯度分量数据的残差；
根据获得的各个梯度分量的反演结果计算各反演结果产生的方差：；
根据如下公式对每个梯度分量数据单独反演的方差进行定权：
；
将权重应用于全张量重力梯度数据残差求取中，利用全张量重力梯度数据进行联合反演获得最终的联合反演结果；
其中，σ2表示单个梯度分量数据反演结果的方差，表示单个梯度分量数据第r个观测数据，/>表示单个梯度分量数据反演结果产生的第r个重力梯度分量数据，R表示测点总数，R的范围为0~108。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410146681.6A CN117688785B (zh) | 2024-02-02 | 2024-02-02 | 一种基于种植思想的全张量重力梯度数据反演方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410146681.6A CN117688785B (zh) | 2024-02-02 | 2024-02-02 | 一种基于种植思想的全张量重力梯度数据反演方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117688785A CN117688785A (zh) | 2024-03-12 |
| CN117688785B true CN117688785B (zh) | 2024-04-16 |
Family
ID=90137503
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202410146681.6A Active CN117688785B (zh) | 2024-02-02 | 2024-02-02 | 一种基于种植思想的全张量重力梯度数据反演方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117688785B (zh) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2020081122A1 (en) * | 2018-10-15 | 2020-04-23 | Illumina, Inc. | Deep learning-based techniques for pre-training deep convolutional neural networks |
| CN115511838A (zh) * | 2022-09-28 | 2022-12-23 | 南京农业大学 | 一种基于群智能优化的植物病害高精度识别方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2008029420A1 (en) * | 2006-09-04 | 2008-03-13 | Geosystem S.R.L. | Method for building velocity models for pre-stack depth migration via the simultaneous joint inversion of seismic, gravity and magnetotelluric data |
| WO2022020102A1 (en) * | 2020-07-22 | 2022-01-27 | The United States Of America, As Represented By The Secretary, Department Of Health And Human Services | Time efficient multi-pulsed field gradient (mpfg) mri without concomitant gradient field artifacts |
-
2024
- 2024-02-02 CN CN202410146681.6A patent/CN117688785B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2020081122A1 (en) * | 2018-10-15 | 2020-04-23 | Illumina, Inc. | Deep learning-based techniques for pre-training deep convolutional neural networks |
| CN115511838A (zh) * | 2022-09-28 | 2022-12-23 | 南京农业大学 | 一种基于群智能优化的植物病害高精度识别方法 |
Non-Patent Citations (4)
| Title |
|---|
| Robust 3D gravity gradient inversion by planting anomalous densities;Leonardo Uieda 等;《GEOPHYSICS》;20120831;第55-66页 * |
| 井地磁异常模量联合反演;李泽林;姚长利;郑元满;;地球物理学报;20181212(第12期);全文 * |
| 某水电站地下厂房初始地应力场反演计算分析;向贵府;刘洋;崔杰;;西南科技大学学报;20171230(第04期);全文 * |
| 重力梯度张量种植反演;曹书锦;朱自强;鲁光银;;中国有色金属学报;20170515(第05期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117688785A (zh) | 2024-03-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111323830B (zh) | 一种基于大地电磁和直流电阻率数据的联合反演方法 | |
| CN112147709B (zh) | 一种基于部分光滑约束的重力梯度数据三维反演方法 | |
| CN109143356A (zh) | 一种自适应混合范数字典学习地震波阻抗反演方法 | |
| CN112363236B (zh) | 一种基于pde的重力场数据等效源延拓与数据类型转换方法 | |
| CN108563837B (zh) | 一种冲积河流水沙模型的模型参数实时校正方法和系统 | |
| CN110133716A (zh) | 基于组合模型加权函数的磁异常数据三维反演方法 | |
| CN108984939A (zh) | 基于3D Gauss-FFT的三维重力场正演方法 | |
| CN115238550A (zh) | 自适应非结构网格的滑坡降雨的地电场数值模拟计算方法 | |
| CN109490978B (zh) | 一种起伏地层的频率域快速高精度正演方法 | |
| CN104316961A (zh) | 获取风化层的地质参数的方法 | |
| CN117688785B (zh) | 一种基于种植思想的全张量重力梯度数据反演方法 | |
| CN105573963A (zh) | 一种电离层水平不均匀结构重构方法 | |
| CN109033588B (zh) | 一种基于空间传播的不确定性量化方法 | |
| CN114200541B (zh) | 一种基于余弦点积梯度约束的三维重磁联合反演方法 | |
| CN104123449A (zh) | 复杂山地区域的分区局部变加密不等距双重网格剖分方法 | |
| CN111859251A (zh) | 一种基于pde的磁测数据等效源上延拓与下延拓方法 | |
| CN111158059A (zh) | 基于三次b样条函数的重力反演方法 | |
| CN112651102B (zh) | 一种起伏地形下的位场全自动极值点深度估算方法 | |
| CN103217715A (zh) | 多尺度规则网格层析反演静校正方法 | |
| CN112748471B (zh) | 一种非结构化等效源的重磁数据延拓与转换方法 | |
| CN112596106B (zh) | 一种球坐标系下重震联合反演密度界面分布的方法 | |
| Vatankhah et al. | A method for 2-dimensional inversion of gravity data | |
| CN111175822A (zh) | 改进直接包络反演与扰动分解的强散射介质反演方法 | |
| CN108132486A (zh) | 一种重磁梯度与地震数据联合界面反演的优化模拟退火法 | |
| CN109188535A (zh) | 地球物理数据处理的方法和装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |