CN113761457A - 基于测量的重力异常数据提取局部重力异常的方法 - Google Patents

基于测量的重力异常数据提取局部重力异常的方法 Download PDF

Info

Publication number
CN113761457A
CN113761457A CN202111055317.1A CN202111055317A CN113761457A CN 113761457 A CN113761457 A CN 113761457A CN 202111055317 A CN202111055317 A CN 202111055317A CN 113761457 A CN113761457 A CN 113761457A
Authority
CN
China
Prior art keywords
gravity anomaly
gravity
anomaly
local
cutting
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.)
Granted
Application number
CN202111055317.1A
Other languages
English (en)
Other versions
CN113761457B (zh
Inventor
孟庆奎
舒晴
高维
周坚鑫
徐光晶
张文志
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Aero Geophysical Survey and Remote Sensing Center for Natural Resources
Original Assignee
China Aero Geophysical Survey and Remote Sensing Center for Natural Resources
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China Aero Geophysical Survey and Remote Sensing Center for Natural Resources filed Critical China Aero Geophysical Survey and Remote Sensing Center for Natural Resources
Priority to CN202111055317.1A priority Critical patent/CN113761457B/zh
Publication of CN113761457A publication Critical patent/CN113761457A/zh
Application granted granted Critical
Publication of CN113761457B publication Critical patent/CN113761457B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Geophysics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

一种基于测量的重力异常数据提取局部重力异常的方法，包括以下步骤：根据重力异常数据Z(x，y)和切割半径r计算与点(x，y)相距r的某4点重力异常的平均值A(x，y)；计算加权参数e；比较区域重力异常与局部重力异常的曲线弯曲方向；对于弯曲方向相反的情形，通过切割区域重力异常修正量进行修正，确定切割半径r区域重力异常；计算并提取局部重力异常L(x，y)＝Z(x，y)‑R(x，y)。该方法可以改善和提高局部重力异常的准确度。还提供了一种基于该局部重力异常的油气勘探方法。

Description

基于测量的重力异常数据提取局部重力异常的方法
技术领域
本发明属于重力异常数据处理以及油气勘探的技术领域，涉及一种基于测量的重力异常数据提取局部重力异常的方法。
背景技术
一般获得的重力数据为布格重力异常数据，它是当代重力研究中获取最为简单且研究最为广泛的重力异常数据。布格重力异常是由区域重力异常和局部重力异常叠加组合而成。其中，由分布范围较大且影响较深的地质因素(例如：区域地质构造)所引起的重力异常称为区域重力异常，由地质因素(例如：矿产，油气)所引起的小范围的重力异常称为局部重力异常。由此可见，局部重力异常的提取不仅是资源勘探的重要研究内容，同时也有助于地质构造的研究。在此基础下，如何准确地提取出局部重力异常也就成为重要问题。
插值切割法提取局部重力异常是现有技术中常用的方法之一，其以当前计算点场值与周围多点平均值的插值运算为切割算子，通过连续切割，得到切割区域场(区域重力异常)，进而将其从位场异常中减去，便得到切割局部场。根据插值切割算子的不同，依次形成了基于四点圆周插值切割算子、八点窗口插值切割算子和动态改进型插值切割算子的插值切割方法。
然而，现有方法的准确性有待进一步提升。
因此，有必要研究一种基于测量的重力异常数据提取局部重力异常的方法来解决上述的一个或多个技术问题。
发明内容
为解决上述至少一个技术问题，申请人通过研究发现，现有的插值切割法忽视了产生异常的不同地质属性，得到的计算结果经常会产生一定的误差。对此，申请人通过从实际地质属性出发，结合重力位场局部异常(局部重力异常)和区域异常(区域重力异常)的四种叠合关系，对插值切割法加以改进，进而提出一种基于异常约束条件的插值切割法，提高了获取的局部重力异常的准确度。
根据本发明一方面，提供了一种基于测量的重力异常数据提取局部重力异常的方法，其特征在于包括以下步骤：
通过航空、航天、船载或地面平台测量重力异常数据Z(x,y)，该重力异常数据Z(x，y)由区域重力异常R(x，y)和局部重力异常L(x，y)组成；
确定切割半径r；
根据重力异常数据Z(x，y)和切割半径r计算与点(x，y)相距r的某4点重力异常的平均值A(x，y)；
计算加权参数e；
比较区域重力异常与局部重力异常的曲线弯曲方向；
对于弯曲方向相反的情形，区域重力异常Rn(x，y)确定为：
Figure BDA0003254369000000021
ΔZ(x，y)为切割区域重力异常修正量，第1次切割的区域重力异常用 R1(x，y)表示，令R1(x，y)作为下次迭代的重力异常数据Z(x，y)，通过该迭代的重力异常数据Z(x，y)确定第2次切割的区域重力异常R2(x，y)；依次迭代下去，最终有
Figure BDA0003254369000000022
于是切割半径r区域重力异常
Figure BDA0003254369000000023
计算并提取局部重力异常L(x，y)＝Z(x，y)-R(x，y)。
根据本发明又一方面，计算加权参数e具体为；
e＝c+d(0≤e≤2)
Figure BDA0003254369000000031
(0≤c≤1；当
Figure BDA0003254369000000032
时，c＝1)
Figure BDA0003254369000000033
(0≤d≤1；当
Figure BDA0003254369000000034
时，d＝1)
Figure BDA0003254369000000035
Figure BDA0003254369000000036
ΔBx＝Z(x+r，y)-Z(x-r，y)
ΔBy＝Z(x，y+r)-Z(x，y-r)。
根据本发明又一方面，所述比较区域重力异常与局部重力异常的曲线弯曲方向具体包括：计算比较参数
Bx1＝Z(x+r，y+r)-[Z(x+2r，y+r)+Z(x，y+r)]/2
By1＝Z(x+r，y+r)-[Z(x+r，y+2r)+Z(x+r，y)]/2
Bx2＝Z(x+2r，y+r)-[Z(x+4r，y+r)+Z(x，y+r)]/2
By2＝Z(x+r，y+2r)-[Z(x+r，y+4r)+Z(x+r，y)]/2
当比较参数同时满足下列条件时，表明区域重力异常和局部重力异常的曲线弯曲方向相反：By1×Bx1＞0，By1×By2＜0，Bx1×Bx2＜0， By2×Bx2＞0。
根据本发明又一方面，当比较参数不能同时满足下列条件时，表明区域重力异常和局部重力异常的曲线弯曲方向相同：By1×Bx1＞0， By1×By2＜0，Bx1×Bx2＜0，By2×Bx2＞0。
根据本发明又一方面，对于弯曲方向相同的情形，区域重力异常
Figure BDA0003254369000000041
根据本发明又一方面，切割区域重力异常修正量ΔZ(x，y)通过以下步骤确定：
考虑x方向重力异常剖面，假设局部重力异常附近的区域重力异常为曲率半径为R、曲率中心为o的曲率圆的部分弧段，则x方向切割区域重力异常修正量为：
Figure BDA0003254369000000042
其中LAC可根据勾股定理求得
Figure BDA0003254369000000043
曲率半径R可根据曲率K的倒数求得
Figure BDA0003254369000000044
用一阶差分替代一阶导数，并忽略一阶小量，则
Figure BDA0003254369000000045
用二阶差分替代二阶导数，并忽略二阶小量，则
Figure BDA0003254369000000046
同理得y方向切割区域重力异常修正量ΔZ(y)
ΔZ(x，y)＝(ΔZ(x)+ΔZ(y))/2；
其中，LOC为曲率中心0到测点A与测点B中间点C的距离，LAC为测点A到测点A与测点B中间点C的距离，XAB、XBE为2倍的切割半径，ZAB为测点A和测点B重力异常的差值，Z(B)为测点B的重力异常值，Z (A)为测点A的重力异常值，Z(E)为测点E的重力异常值。
根据本发明又一方面，还提供了一种油气勘探方法，其特征在于包括以下步骤：
根据前述的方法提取局部重力异常数据；
根据所述局部重力异常数据进行油气勘探。
本发明可以获得以下一个或多个技术效果：
通过从实际地质属性出发，结合重力位场局部异常(局部重力异常)和区域异常(区域重力异常)的四种叠合关系，对插值切割法加以改进，提高了获取的局部重力异常的准确度。
附图说明
下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1为局部异常与区域异常叠合关系分类示意图，其中，粗实线表示总场异常，AB之间的细实线表示区域异常，虚线为局部异常两端点连线，(a)-(d)分别对应表1中分类情形1-4。
图2为根据本发明的一种优选实施例的基于测量的重力异常数据提取局部重力异常的方法的切割区域异常修正量(切割区域重力异常修正量)示意图。
图3为根据本发明的一种优选实施例的基于测量的重力异常数据提取局部重力异常的方法的理论模型示意图。
图4为本发明的理论模型正演结果平面等值线图。
图5为本发明的理论模型正演结果剖面图(y＝0)。
图6为本发明的理论模型重力位场分离效果对比平面图。
图7为本发明的理论模型重力位场分离效果对比剖面图(y＝0)。
图8为根据现有技术的插值切割法计算和提取得到的上地幔顶部附近重力异常。
图9为根据本发明方法(约束插值切割法)计算和提取得到的上地幔顶部附近重力异常。
具体实施方式
下面结合附图，通过优选实施例来描述本发明的最佳实施方式，这里的具体实施方式在于详细地说明本发明，而不应理解为对本发明的限制，在不脱离本发明的精神和实质范围的情况下，可以做出各种变形和修改，这些都应包含在本发明的保护范围之内。
实施例1
申请人通过研究发现，现有的插值切割法忽视了产生异常的不同地质属性，得到的计算结果经常会产生一定的误差。对此，申请人通过从实际地质属性出发，结合重力位场局部异常(局部重力异常)和区域异常(区域重力异常)的四种叠合关系，对插值切割法加以改进，进而提出一种基于异常约束条件的插值切割法，提高了获取的局部重力异常的准确度。
根据本发明一优选实施方式，提供了一种基于测量的重力异常数据提取局部重力异常的方法，其特征在于包括以下步骤：
通过航空、航天、船载或地面平台测量重力异常数据Z(x,y)，该重力异常数据Z(x，y)由区域重力异常R(x，y)和局部重力异常L(x，y)组成；
确定切割半径r；
根据重力异常数据Z(x，y)和切割半径r计算与点(x，y)相距r的某4点重力异常的平均值A(x，y)；
计算加权参数e；
比较区域重力异常与局部重力异常的曲线弯曲方向；
对于弯曲方向相反的情形，区域重力异常Rn(x，y)确定为：
Figure BDA0003254369000000071
ΔZ(x，y)为切割区域重力异常修正量，第1次切割的区域重力异常用 R1(x，y)表示，令R1(x，y)作为下次迭代的重力异常数据Z(x，y)，即 Z(x，y)＝R1(x，y)，通过该迭代的重力异常数据Z(x，y)确定第2次切割的区域重力异常R2(x，y)；依次迭代下去，最终有
Figure BDA0003254369000000072
于是切割半径r区域重力异常
Figure BDA0003254369000000073
计算并提取局部重力异常L(x，y)＝Z(x，y)-R(x，y)。
根据本发明又一优选实施方式，计算加权参数e具体为；
e＝c+d(0≤e≤2)
Figure BDA0003254369000000074
(0≤c≤1；当
Figure BDA0003254369000000075
时，c＝1)
Figure BDA0003254369000000076
(0≤d≤1；当
Figure BDA0003254369000000077
时，d＝1)
Figure BDA0003254369000000078
Figure BDA0003254369000000079
ΔBx＝Z(x，r，y)-Z(x-r，y)ΔBy＝Z(x，y+r)-Z(x，y-r)。
根据本发明又一优选实施方式，所述比较区域重力异常与局部重力异常的曲线弯曲方向具体包括：计算比较参数
Bx1＝Z(x+r，y+r)-[Z(x+2r，y+r)+Z(x，y+r)]/2
By1＝Z(x+r，y+r)-[Z(x+r，y+2r)+Z(x+r，y)]/2
Bx2＝Z(x+2r，y+r)-[Z(x+4r，y+r)+Z(x，y+r)]/2
By2＝Z(x+r，y+2r)-[Z(x+r，y+4r)+Z(x，y+r)]/2
当比较参数同时满足下列条件时，表明区域重力异常和局部重力异常的曲线弯曲方向相反：By1×Bx1＞0，By1×Bx1＜0，By1×Bx1＜0， By2×Bx2＞0。
根据本发明又一优选实施方式，当比较参数不能同时满足下列条件时，表明区域重力异常和局部重力异常的曲线弯曲方向相同：By1×Bx1＞0， By1×By2＜0，Bx1×Bx2＜0，Ry2×By2＞0。
根据本发明又一优选实施方式，对于弯曲方向相同的情形，区域重力异常
Figure BDA0003254369000000081
根据本发明又一优选实施方式，切割区域重力异常修正量ΔZ(x，y)通过以下步骤确定：
考虑x方向重力异常剖面，假设局部重力异常附近的区域重力异常为曲率半径为R、曲率中心为o的曲率圆的部分弧段，则x方向切割区域重力异常修正量为：
Figure BDA0003254369000000091
其中LAC可根据勾股定理求得
Figure BDA0003254369000000092
曲率半径R可根据曲率K的倒数求得
Figure BDA0003254369000000093
用一阶差分替代一阶导数，并忽略一阶小量，则
Figure BDA0003254369000000094
用二阶差分替代二阶导数，并忽略二阶小量，则
Figure BDA0003254369000000095
同理得y方向切割区域重力异常修正量ΔZ(y)
ΔZ(x，y)＝(ΔZ(x)+ΔZ(y))/2；
其中，LOC为曲率中心0到测点A与测点B中间点C的距离，LAC为测点A到测点A与测点B中间点C的距离，XAB、XBE为2倍的切割半径，ZAB为测点A和测点B重力异常的差值，Z(B)为测点B的重力异常值，Z (A)为测点A的重力异常值，Z(E)为测点E的重力异常值。
根据本发明又一优选实施方式，还提供了一种油气勘探方法，其特征在于包括以下步骤：
根据权利要求1-6任一项所述的方法提取局部重力异常数据；
根据所述局部重力异常数据进行油气勘探。
实施例2
在实施例1的基础上，本实施例进一步通过详细介绍本发明的技术原理与实际应用效果。
根据本发明一种优选实施方式，还提供了一种基于测量的重力异常数据提取局部重力异常的方法，其特征在于包括以下步骤：
(1)确定切割半径r和测量重力异常Z(x，y)，Z(x，y)由区域场R(x，y)和局部场L(x，y)组成，即Z(x，y)＝R(x，y)+L(x，y)。
(2)计算与点(x，y)相距r的某4点重力异常的平均值A(x，y)
Figure BDA0003254369000000101
(3)计算加权参数e
e＝c+d(0≤e≤2)
Figure BDA0003254369000000102
(0≤c≤1；当
Figure BDA0003254369000000103
时，c＝1)
Figure BDA0003254369000000104
(0≤d≤1；当
Figure BDA0003254369000000105
时，d＝1)
Figure BDA0003254369000000106
Figure BDA0003254369000000111
ΔBx＝Z(x+r，y)-Z(x-r，y)
ΔBy＝Z(x，y+r)-Z(x，y-r)
(4)比较区域异常与局部异常曲线弯曲方向，具体可细分为2步：
(4-1)计算比较参数
Bx1＝Z(x+r，y+r)-[Z(x+2r，y+r)+Z(x，y+r)]/2
By1＝Z(x+r，y+r)-[Z(x+r，y+2r)+Z(x+r，y)]/2
Bx2＝Z(x+2r，y+r)-[Z(x+4r，y+r)+Z(x，y+r)]/2
By2＝Z(x+r，y+2r)-[Z(x+r，y+4r)+Z(x+r，y)]/2
(4-2)当比较参数不能同时满足下列条件时，表明区域异常和局部异常弯曲方向相同(即分类情形1和2)：By1×Bx1＞0，By1×By2＜0， Bx1×Bx2＜0，By2×Bx2＞0。
(4-3)当比较参数同时满足下列条件时，表明区域异常和局部异常弯曲方向相反(即分类情形3和4)：By1×Bx1＞0，By1×By2＜0， Bx1×Bx2＜0，By2×Bx2＞0。
(5)对于弯曲方向相同的情形，区域场是Z(x，y)与A(x，y)的加权平均，采用下式计算切割区域场：
Figure BDA0003254369000000112
(6)对于弯曲方向相反的情形，区域场是Z(x，y)与A(x，y)的加权平均的基础上增加切割区域异常修正量ΔZ(x，y)，此种情况可细分为以下3步计算区域异常：
(6-1)首先考虑x方向异常剖面，假设局部异常附近的区域异常为曲率半径为R、曲率中心为o的曲率圆的部分弧段，则x方向切割区域异常修正量为：
Figure BDA0003254369000000121
其中LAC可根据勾股定理求得
Figure BDA0003254369000000122
曲率半径R可根据曲率K的倒数求得
Figure BDA0003254369000000123
用一阶差分替代一阶导数，并忽略一阶小量，则
Figure BDA0003254369000000124
用二阶差分替代二阶导数，并忽略二阶小量，则
Figure BDA0003254369000000125
(6-2)同理得y方向切割区域异常修正量ΔZ(y)
(6-3)取ΔZ(x，y)＝(ΔZ(x)+ΔZ(y))/2
(6-4)计算修正量ΔZ(x，y)后，采用下式计算切割区域场：
Figure BDA0003254369000000131
参见图2，其中，LOC为曲率中心0到测点A与测点B中间点C的距离，LAC为测点A到测点A与测点B中间点C的距离，XAB、XBE为2倍的切割半径，ZAB为测点A和测点B重力异常的差值，Z(B)为测点B的重力异常值，Z(A)为测点A的重力异常值，Z(E)为测点E的重力异常值。
优选地，参见图2，测点D为当前计算点，测点A为间隔切割半径r的后向测点，测点B为间隔切割半径r的前向测点，测点E为间隔3倍切割半径的前向测点。也就是说，测点A、D、B两两之间间隔切割半径r，B、E 之间间隔2r。例如，测点A与D之间的间隔距离为在x方向上的距离r。如图2所示，测点A可作为起算点，其为剖面起始点(第一个点)。
(7)用上述方法得到的区域场称为第1次切割的区域场，用R1(x，y)表示，对R1(x，y)重复使用步骤(1)～(6)，得到第2次切割的区域场 R2(x，y)；依次迭代下去，最终有
Figure BDA0003254369000000132
于是有
Figure BDA0003254369000000133
称这个区域场为切割半径r区域场。
(8)计算局部异常，将区域场与重力异常Z(x，y)相减，得到局部场：
L(x，y)＝Z(x，y)-R(x，y)。
优选地，通过航空、航天、船载和地面等不同平台获取的重力异常 Z(x，y)。
参见图1，其示出了局部异常(局部重力异常)与区域异常(区域重力异常)叠合关系分类四种情形的示意图。
申请人研究发现，现有的插值切割法仅适用于分类情形1和2，即区域异常介于总场异常值和四点圆周平均值的取值区间。然而对于区域异常和局部异常曲线弯曲方向相反的情形，即分类情形3和4，其计算结果虽可无限向四点圆周平均值逼近，但其忽视了图1中的(c)、(d)中黑色竖线所标注的部分。
表1局部异常与区域异常叠合关系分类表
Figure BDA0003254369000000141
注：四种分类情形分别对应图1中的(a)、(b)、(c)、(d)。
优选地，对于弯曲方向相反的情形(即分类情形3和4)，通过切割区域异常修正量ΔZ(x，y)来进行补偿，如图2所示。
接下来，进一步介绍理论模型试算的情况。
首先，进行理论模型设计。具体地，设计了图3所示的理论模型，三个深部球体产生区域重力高，浅部球体1产生局部重力低，浅部球体2和3产生局部重力高，理论模型详细参数见表2所示。
表2三维模型异常体的空间位置和物性参数表
Figure BDA0003254369000000151
接着，进行了理论模型正演。本发明设计的理论模型产生的重力异常为各球体剩余质量的引力位沿重力方向(定义为z方向)的导数之和，即
Figure BDA0003254369000000152
上式中Vi为球体i产生的引力位，f为引力常数， f＝6.672×10-11m3/(kg·s2),mi为球体i的剩余质量，
Figure BDA0003254369000000153
Ri为球体i 的半径，σi为球体i的剩余密度，ri为测点与球体i的球心的距离，
Figure BDA0003254369000000154
基于公式(25)和设计的理论模型参数，经过正演计算可得到图4、图5所示的正演结果。由图4-5可见，1号异常属于分类情形1，2、3号属于分类情形3。
进一步，进行理论模型重力位场分离与提取效果对比分析。分别采用现有技术中的插值切割法和本发明的方法(约束插值切割法)对设计的理论模型正演结果进行了重力位场分离与提取，在计算过程中采用了多个切割半径r，经筛选对比确定r＝700m。重力位场分离的平面效果见图6所示，对比可见约束插值切割法得到的局部异常(图6f)更接近于理论局部异常(图6d)，由于区域异常整体呈现出较大的高值背景，基本淹没了局部异常的存在，故通过两种方法得到的区域异常与理论区域异常表现出很高的相似度(图6a、6b、6c)。为了更加直观体现本发明方法的优越性，进一步绘制了主剖面(y＝0)两种方法得到的位场分离曲线(图7)，由于1号异常为区域重力高叠加局部重力高(分类情形1)，故两种方法得到了一致的结果，但针对2号异常、3号异常，可见本发明方法得到的分离与提取结果得到了明显改善。
进一步，为了定量分析本发明提出的约束插值切割法对插值切割法的优化效果，引入“改善率”(缩写为IR)来估计前者对局部异常的优化量：
IR＝var(Δgcz-Δgll)/var(Δgys-Δgll) (26)
式中Δgcz表示插值切割法得到的局部异常，Δgys表示约束插值切割法得到的局部异常，Δgll为理论局部异常，var表示均方差。
对于2号、3号局部异常的平面区域，经计算改善率为2.15，据此可见相比于插值切割法，本发明方法提取的局部异常改善了115％，优化效果明显。
进一步，以中国及周边区域的重力数据为试验对象，范围为东经60°～141°，北纬15°～59°的扇形区域。重力数据来自超高阶地球重力场模型EGM2008(EarthGravitational Model 2008),该模型的阶次完全至2159次，空间分辨率约为5′，在中国80％以上的面积上的精度可达10mGal之内，可用于小比例尺重力编图和构造研究。本实例使用的布格重力异常的网格间距为4km。
分别利用插值切割法和本发明提出的约束插值切割法对布格重力异常进行了分离，切割半径选取18和19倍网格间距，计算得到72～76km参考深度层的重力异常值，并将其视为上地幔顶部附近介质产生的重力异常(图8，图9)，其中AB剖面位置为北纬43.76°。
对比两种方法得到的重力异常平面图(图8(a)，图9(a))，可见整幅图重力异常的总体形态相似，但本发明(约束插值切割法)得到的正值重力异常幅度和范围相对削减，甚至消失，得到的负值重力异常幅度和范围相对增加。AB剖面上的对比结果(如图8(b)，图9(b))更加直观，实现了对插值切割法计算结果的修正，得到了更加准确的重力异常。
本发明可以获得以下一个或多个技术效果：
通过从实际地质属性出发，结合重力位场局部异常(局部重力异常)和区域异常(区域重力异常)的四种叠合关系，对插值切割法加以改进，提高了获取的局部重力异常的准确度。
本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护的范围由所附的权利要求书及其等效物界定。

Claims (7)

1.一种基于测量的重力异常数据提取局部重力异常的方法，其特征在于包括以下步骤：
通过航空、航天、船载或地面平台测量重力异常数据Z(x,y)，该重力异常数据Z(x，y)由区域重力异常R(x，y)和局部重力异常L(x，y)组成；
确定切割半径r；
根据重力异常数据Z(x，y)和切割半径r计算与点(x，y)相距r的某4点重力异常的平均值A(x，y)；
计算加权参数e；
比较区域重力异常与局部重力异常的曲线弯曲方向；
对于弯曲方向相反的情形，区域重力异常Rn(x，y)确定为：
Figure FDA0003254368990000011
ΔZ(x，y)为切割区域重力异常修正量，第1次切割的区域重力异常用R1(x，y)表示，令R1(x，y)作为下次迭代的重力异常数据Z(x，y)，通过该迭代的重力异常数据Z(x，y)确定第2次切割的区域重力异常R2(x，y)；依次迭代下去，最终有
Figure FDA0003254368990000012
于是切割半径r区域重力异常
Figure FDA0003254368990000013
计算并提取局部重力异常L(x，y)＝Z(x，y)-R(x，y)。
2.根据权利要求1所述的基于测量的重力异常数据提取局部重力异常的方法，其特征在于计算加权参数e具体为；
e＝c+d (0≤e≤2)
Figure FDA0003254368990000021
(0≤c≤1；当
Figure FDA0003254368990000022
时，c＝1)
Figure FDA0003254368990000023
(0≤d≤1；当
Figure FDA0003254368990000024
时，d＝1)
Figure FDA0003254368990000025
Figure FDA0003254368990000026
ΔBx＝Z(x+r，y)-Z(x-r，y)
ΔBy＝Z(x，y+r)-Z(x，y-r)。
3.根据权利要求2所述的基于测量的重力异常数据提取局部重力异常的方法，其特征在于所述比较区域重力异常与局部重力异常的曲线弯曲方向具体包括：计算比较参数
Bx1＝Z(x+r，y+r)-[Z(x+2r，y+r)+Z(x，y+r)]/2
By1＝Z(x+r，y+r)-[Z(x+r，y+2r)+Z(x+r，y)]/2
Bx2＝Z(x+2r，y+r)-[Z(x+4r，y+r)+Z(x，y+r)]/2
By2＝Z(x+r，y+2r)-[Z(x+r，y+4r)+Z(x+r，y)]/2
当比较参数同时满足下列条件时，表明区域重力异常和局部重力异常的曲线弯曲方向相反：By1×Bx1＞0，By1×By2＜0，Bx1×Bx2＜0，By2×Bx2＞0。
4.根据权利要求3所述的基于测量的重力异常数据提取局部重力异常的方法，其特征在于当比较参数不能同时满足下列条件时，表明区域重力异常和局部重力异常的曲线弯曲方向相同：By1×Bx1＞0，By1×By2＜0，Bx1×Bx2＜0，By2×Bx2＞0。
5.根据权利要求4所述的基于测量的重力异常数据提取局部重力异常的方法，其特征在于对于弯曲方向相同的情形，区域重力异常
Figure FDA0003254368990000031
6.根据权利要求1-5任一项所述的基于测量的重力异常数据提取局部重力异常的方法，其特征在于切割区域重力异常修正量ΔZ(x，y)通过以下步骤确定：
考虑x方向重力异常剖面，假设局部重力异常附近的区域重力异常为曲率半径为R、曲率中心为o的曲率圆的部分弧段，则x方向切割区域重力异常修正量为：
Figure FDA0003254368990000032
其中LAC可根据勾股定理求得
Figure FDA0003254368990000041
曲率半径R可根据曲率K的倒数求得
Figure FDA0003254368990000042
用一阶差分替代一阶导数，并忽略一阶小量，则
Figure FDA0003254368990000043
用二阶差分替代二阶导数，并忽略二阶小量，则
Figure FDA0003254368990000044
同理得y方向切割区域重力异常修正量ΔZ(y)
ΔZ(x，y)＝(ΔZ(x)+ΔZ(y))/2；
其中，LOC为曲率中心0到测点A与测点B中间点C的距离，LAC为测点A到测点A与测点B中间点C的距离，XAB、XBE为2倍的切割半径，ZAB为测点A和测点B重力异常的差值，Z(B)为测点B的重力异常值，Z(A)为测点A的重力异常值，Z(E)为测点E的重力异常值。
7.一种油气勘探方法，其特征在于包括以下步骤：
根据权利要求1-6任一项所述的方法提取局部重力异常数据；
根据所述局部重力异常数据进行油气勘探。
CN202111055317.1A 2021-09-09 2021-09-09 基于测量的重力异常数据提取局部重力异常的方法 Active CN113761457B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111055317.1A CN113761457B (zh) 2021-09-09 2021-09-09 基于测量的重力异常数据提取局部重力异常的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111055317.1A CN113761457B (zh) 2021-09-09 2021-09-09 基于测量的重力异常数据提取局部重力异常的方法

Publications (2)

Publication Number Publication Date
CN113761457A true CN113761457A (zh) 2021-12-07
CN113761457B CN113761457B (zh) 2024-05-14

Family

ID=78794292

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111055317.1A Active CN113761457B (zh) 2021-09-09 2021-09-09 基于测量的重力异常数据提取局部重力异常的方法

Country Status (1)

Country Link
CN (1) CN113761457B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115629426A (zh) * 2022-10-24 2023-01-20 中国自然资源航空物探遥感中心 一种提取局部重力异常的方法、***、装置及介质

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2007146867A (ru) * 2007-12-17 2009-06-27 Горный институт Уральского отделения Российской академии наук (ГИ УрО РАН) (RU) Способ многокомпонентного гравиметрического моделирования геологической среды
US20180128940A1 (en) * 2016-11-04 2018-05-10 Robert J. Ferderer Global Inversion of Gravity Data Using the Principle of General Local Isostasy for Lithospheric Modeling
CN109557594A (zh) * 2018-12-11 2019-04-02 中国人民解放***箭军工程大学 基于重力异常时变的重力基准图时变修正方法及***
CN110031000A (zh) * 2019-05-21 2019-07-19 北京理工大学 一种重力辅助惯性导航区域适配性的评价方法
CN110706275A (zh) * 2019-10-16 2020-01-17 中国石油大学(华东) 一种基于卫星测高重力数据的局部重力异常提取方法
CN111337993A (zh) * 2020-03-30 2020-06-26 中国地质科学院 一种基于变密度变深度约束的重力密度界面反演方法
CN112464520A (zh) * 2020-10-28 2021-03-09 中国石油天然气集团有限公司 局部重力异常深度反演方法及装置
CN112596106A (zh) * 2020-11-09 2021-04-02 中国人民武装警察部队黄金第五支队 一种球坐标系下重震联合反演密度界面分布的方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2007146867A (ru) * 2007-12-17 2009-06-27 Горный институт Уральского отделения Российской академии наук (ГИ УрО РАН) (RU) Способ многокомпонентного гравиметрического моделирования геологической среды
US20180128940A1 (en) * 2016-11-04 2018-05-10 Robert J. Ferderer Global Inversion of Gravity Data Using the Principle of General Local Isostasy for Lithospheric Modeling
CN109557594A (zh) * 2018-12-11 2019-04-02 中国人民解放***箭军工程大学 基于重力异常时变的重力基准图时变修正方法及***
CN110031000A (zh) * 2019-05-21 2019-07-19 北京理工大学 一种重力辅助惯性导航区域适配性的评价方法
CN110706275A (zh) * 2019-10-16 2020-01-17 中国石油大学(华东) 一种基于卫星测高重力数据的局部重力异常提取方法
CN111337993A (zh) * 2020-03-30 2020-06-26 中国地质科学院 一种基于变密度变深度约束的重力密度界面反演方法
CN112464520A (zh) * 2020-10-28 2021-03-09 中国石油天然气集团有限公司 局部重力异常深度反演方法及装置
CN112596106A (zh) * 2020-11-09 2021-04-02 中国人民武装警察部队黄金第五支队 一种球坐标系下重震联合反演密度界面分布的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
孙建国: "复杂地表条件下地球物理场数值模拟方法评述", 世界地质, vol. 03, 31 December 2007 (2007-12-31), pages 345 - 362 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115629426A (zh) * 2022-10-24 2023-01-20 中国自然资源航空物探遥感中心 一种提取局部重力异常的方法、***、装置及介质

Also Published As

Publication number Publication date
CN113761457B (zh) 2024-05-14

Similar Documents

Publication Publication Date Title
CN108489402B (zh) 基于三维激光扫描的露天矿山边坡岩体节理规模快速精细取值方法
CN107227950B (zh) 一种实钻井眼轨迹整体性评价方法
CN111323776B (zh) 一种矿区形变的监测方法
Piccini Recent developments on morphometric analysis of karst caves
CN109884710B (zh) 针对激发井深设计的微测井层析成像方法
CN112785596B (zh) 基于dbscan聚类的点云图螺栓分割和高度测量方法
CN113761457A (zh) 基于测量的重力异常数据提取局部重力异常的方法
CA2402887A1 (en) Method for characterization of multi-scale geometric attributes
CN110414060A (zh) 一种基于四阶谱矩的位场边界识别方法
Wang et al. Fragmentation calculation method for blast muck piles in open-pit copper mines based on three-dimensional laser point cloud data
CN106873031B (zh) 一种三维地震观测***垂向分辨率定量分析评价方法
CN103149130A (zh) 砾岩岩心颗粒结构中粒度的分析方法
CN1869734A (zh) 重力勘探数据处理变密度地形校正方法
CN110954958A (zh) 裂缝与断层的预测方法及***
CN112596113A (zh) 一种基于重力不同阶梯度特征值交点的场源位置识别方法
CN109374436B (zh) 一种基于ct图像的水合物沉积物剪切带识别方法
CN113779888B (zh) 地面沉降危险性评估方法、装置、设备及存储介质
CN107179548B (zh) 一种基于真地表的叠前地震成像方法
CN107239629B (zh) 一种岩石结构面实验室合理尺寸确定的分形维数分析方法
CN112859073B (zh) 一种基于PSInSAR技术的道路破损评估方法
CN110989032B (zh) 一种基于倾斜角的重力水平总梯度断裂识别方法
CN104166162B (zh) 基于迭代三参数小波变换的缝洞发育带检测方法
CN113325474B (zh) 生物礁判别方法
CN114152258B (zh) 基于地球物理/几何特征的海洋多场多参数定位融合方法
CN112505752B (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