基于测量的重力异常数据提取局部重力异常的方法
阅读说明：本技术 基于测量的重力异常数据提取局部重力异常的方法 (Method for extracting local gravity anomaly based on measured gravity anomaly data ) 是由 孟庆奎 舒晴 高维 周坚鑫 徐光晶 张文志 于 2021-09-09 设计创作，主要内容包括：一种基于测量的重力异常数据提取局部重力异常的方法,包括以下步骤：根据重力异常数据Z(x,y)和切割半径r计算与点(x,y)相距r的某4点重力异常的平均值A(x,y)；计算加权参数e；比较区域重力异常与局部重力异常的曲线弯曲方向；对于弯曲方向相反的情形,通过切割区域重力异常修正量进行修正,确定切割半径r区域重力异常；计算并提取局部重力异常L(x,y)＝Z(x,y)-R(x,y)。该方法可以改善和提高局部重力异常的准确度。还提供了一种基于该局部重力异常的油气勘探方法。(A method for extracting local gravity anomaly based on measured gravity anomaly data comprises the following steps: calculating the average value A (x, y) of the gravity anomaly of a certain 4 points which are away from the point (x, y) by r according to the gravity anomaly data Z (x, y) and the cutting radius r; calculating a weighting parameter e; comparing the curve bending directions of the regional gravity anomaly and the local gravity anomaly; for the situation that the bending direction is opposite, correcting through the gravity anomaly correction quantity of the cutting area, and determining that the gravity of the area with the cutting radius r is abnormal; calculating and extracting local gravity anomaly L (x, y) ═ Z (x, y) -R (x, y). The method can improve and improve the accuracy of local gravity anomaly. A method of oil and gas exploration based on the local gravity anomaly is also provided.)
技术领域
本发明属于重力异常数据处理以及油气勘探的技术领域，涉及一种基于测量的重力异常数据提取局部重力异常的方法。
背景技术
一般获得的重力数据为布格重力异常数据，它是当代重力研究中获取最为简单且研究最为广泛的重力异常数据。布格重力异常是由区域重力异常和局部重力异常叠加组合而成。其中，由分布范围较大且影响较深的地质因素(例如：区域地质构造)所引起的重力异常称为区域重力异常，由地质因素(例如：矿产，油气)所引起的小范围的重力异常称为局部重力异常。由此可见，局部重力异常的提取不仅是资源勘探的重要研究内容，同时也有助于地质构造的研究。在此基础下，如何准确地提取出局部重力异常也就成为重要问题。
插值切割法提取局部重力异常是现有技术中常用的方法之一，其以当前计算点场值与周围多点平均值的插值运算为切割算子，通过连续切割，得到切割区域场(区域重力异常)，进而将其从位场异常中减去，便得到切割局部场。根据插值切割算子的不同，依次形成了基于四点圆周插值切割算子、八点窗口插值切割算子和动态改进型插值切割算子的插值切割方法。
然而，现有方法的准确性有待进一步提升。
因此，有必要研究一种基于测量的重力异常数据提取局部重力异常的方法来解决上述的一个或多个技术问题。
发明内容
为解决上述至少一个技术问题，申请人通过研究发现，现有的插值切割法忽视了产生异常的不同地质属性，得到的计算结果经常会产生一定的误差。对此，申请人通过从实际地质属性出发，结合重力位场局部异常(局部重力异常)和区域异常(区域重力异常)的四种叠合关系，对插值切割法加以改进，进而提出一种基于异常约束条件的插值切割法，提高了获取的局部重力异常的准确度。
根据本发明一方面，提供了一种基于测量的重力异常数据提取局部重力异常的方法，其特征在于包括以下步骤：
通过航空、航天、船载或地面平台测量重力异常数据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)确定为：
ΔZ(x，y)为切割区域重力异常修正量，第1次切割的区域重力异常用 R1(x，y)表示，令R1(x，y)作为下次迭代的重力异常数据Z(x，y)，通过该迭代的重力异常数据Z(x，y)确定第2次切割的区域重力异常R2(x，y)；依次迭代下去，最终有
于是切割半径r区域重力异常
计算并提取局部重力异常L(x，y)＝Z(x，y)-R(x，y)。
根据本发明又一方面，计算加权参数e具体为；
e＝c+d(0≤e≤2)
(0≤c≤1；当时，c＝1)
(0≤d≤1；当时，d＝1)
Δ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。
根据本发明又一方面，对于弯曲方向相同的情形，区域重力异常
根据本发明又一方面，切割区域重力异常修正量ΔZ(x，y)通过以下步骤确定：
考虑x方向重力异常剖面，假设局部重力异常附近的区域重力异常为曲率半径为R、曲率中心为o的曲率圆的部分弧段，则x方向切割区域重力异常修正量为：
其中LAC可根据勾股定理求得
曲率半径R可根据曲率K的倒数求得
用一阶差分替代一阶导数，并忽略一阶小量，则
用二阶差分替代二阶导数，并忽略二阶小量，则
同理得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)确定为：
Δ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)；依次迭代下去，最终有
于是切割半径r区域重力异常
计算并提取局部重力异常L(x，y)＝Z(x，y)-R(x，y)。
根据本发明又一优选实施方式，计算加权参数e具体为；
e＝c+d(0≤e≤2)
(0≤c≤1；当时，c＝1)
(0≤d≤1；当时，d＝1)
Δ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。
根据本发明又一优选实施方式，对于弯曲方向相同的情形，区域重力异常
根据本发明又一优选实施方式，切割区域重力异常修正量ΔZ(x，y)通过以下步骤确定：
考虑x方向重力异常剖面，假设局部重力异常附近的区域重力异常为曲率半径为R、曲率中心为o的曲率圆的部分弧段，则x方向切割区域重力异常修正量为：
其中LAC可根据勾股定理求得
曲率半径R可根据曲率K的倒数求得
用一阶差分替代一阶导数，并忽略一阶小量，则
用二阶差分替代二阶导数，并忽略二阶小量，则
同理得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)
(3)计算加权参数e
e＝c+d(0≤e≤2)
(0≤c≤1；当时，c＝1)
(0≤d≤1；当时，d＝1)
Δ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)的加权平均，采用下式计算切割区域场：
(6)对于弯曲方向相反的情形，区域场是Z(x，y)与A(x，y)的加权平均的基础上增加切割区域异常修正量ΔZ(x，y)，此种情况可细分为以下3步计算区域异常：
(6-1)首先考虑x方向异常剖面，假设局部异常附近的区域异常为曲率半径为R、曲率中心为o的曲率圆的部分弧段，则x方向切割区域异常修正量为：
其中LAC可根据勾股定理求得
曲率半径R可根据曲率K的倒数求得
用一阶差分替代一阶导数，并忽略一阶小量，则
用二阶差分替代二阶导数，并忽略二阶小量，则
(6-2)同理得y方向切割区域异常修正量ΔZ(y)
(6-3)取ΔZ(x，y)＝(ΔZ(x)+ΔZ(y))/2
(6-4)计算修正量ΔZ(x，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的重力异常值。
优选地，参见图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)；依次迭代下去，最终有
于是有称这个区域场为切割半径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局部异常与区域异常叠合关系分类表
注：四种分类情形分别对应图1中的(a)、(b)、(c)、(d)。
优选地，对于弯曲方向相反的情形(即分类情形3和4)，通过切割区域异常修正量ΔZ(x，y)来进行补偿，如图2所示。
接下来，进一步介绍理论模型试算的情况。
首先，进行理论模型设计。具体地，设计了图3所示的理论模型，三个深部球体产生区域重力高，浅部球体1产生局部重力低，浅部球体2和3产生局部重力高，理论模型详细参数见表2所示。
表2三维模型异常体的空间位置和物性参数表
接着，进行了理论模型正演。本发明设计的理论模型产生的重力异常为各球体剩余质量的引力位沿重力方向(定义为z方向)的导数之和，即
上式中Vi为球体i产生的引力位，f为引力常数， f＝6.672×10-11m3/(kg·s2),mi为球体i的剩余质量，Ri为球体i 的半径，σi为球体i的剩余密度，ri为测点与球体i的球心的距离，
基于公式(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))更加直观，实现了对插值切割法计算结果的修正，得到了更加准确的重力异常。
本发明可以获得以下一个或多个技术效果：
通过从实际地质属性出发，结合重力位场局部异常(局部重力异常)和区域异常(区域重力异常)的四种叠合关系，对插值切割法加以改进，提高了获取的局部重力异常的准确度。
本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护的范围由所附的权利要求书及其等效物界定。