阅读说明：本技术 基于测量的重力异常数据提取局部重力异常的方法 (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.)

1. A method for extracting local gravity anomaly based on measured gravity anomaly data is characterized by comprising the following steps:

measuring gravity anomaly data Z (x, y) by an aviation, aerospace, shipborne or ground platform, wherein the gravity anomaly data Z (x, y) consists of regional gravity anomaly R (x, y) and local gravity anomaly L (x, y);

determining a cutting radius r;

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 case of opposite bending direction, regional gravity anomaly R_n(x, y) is determined as:

delta Z (x, y) is the correction amount of gravity anomaly in cutting area, and R is used for gravity anomaly in area of 1 st cutting₁(x, y) represents, let R₁(x, y) as gravity anomaly data Z (x, y) of the next iteration, and determining the region gravity anomaly R of the 2 nd cutting according to the gravity anomaly data Z (x, y) of the iteration₂(x, y); successively iterate, finally have

So that the gravity in the region of the cutting radius r is abnormal

Calculating and extracting local gravity anomaly L (x, y) ═ Z (x, y) -R (x, y).

2. The method for extracting local gravity anomaly based on measured gravity anomaly data according to claim 1, wherein the weighting parameter e is calculated;

e＝c+d (0≤e≤2)

(0. ltoreq. c.ltoreq.1; whenThen, c is 1)

(0. ltoreq. d. ltoreq.1; whenThen, d is 1)

ΔB_x＝Z(x+r，y)-Z(x-r，y)

ΔB_y＝Z(x，y+r)-Z(x，y-r)。

3. The method of claim 2, wherein comparing the curve bending direction of the regional gravity anomaly with the curve bending direction of the local gravity anomaly comprises: calculating a comparison parameter

B_x1＝Z(x+r，y+r)-[Z(x+2r，y+r)+Z(x，y+r)]/2

B_y1＝Z(x+r，y+r)-[Z(x+r，y+2r)+Z(x+r，y)]/2

B_x2＝Z(x+2r，y+r)-[Z(x+4r，y+r)+Z(x，y+r)]/2

B_y2＝Z(x+r，y+2r)-[Z(x+r，y+4r)+Z(x+r，y)]/2

When the comparison parameters simultaneously satisfy the following conditions, the curve bending directions indicating the regional gravity anomaly and the local gravity anomaly are opposite: b is_y1×B_x1＞0，B_y1×B_y2＜0，B_x1×B_x2＜0，B_y2×B_x2＞0。

4. The method for extracting local gravity anomaly based on measured gravity anomaly data according to claim 3, wherein the curve bending directions indicating the regional gravity anomaly and the local gravity anomaly are the same when the comparison parameters cannot satisfy the following conditions simultaneously: b is_y1×B_x1＞0，B_y1×B_y2＜0，B_x1×B_x2＜0，B_y2×B_x2＞0。

5. The method for extracting local gravity anomaly based on measured gravity anomaly data according to claim 4, wherein for the same bending direction, the regional gravity anomaly

6. The method for extracting local gravity anomaly based on measured gravity anomaly data according to any one of claims 1-5, characterized in that the cutting zone gravity anomaly correction quantity Δ Z (x, y) is determined by:

considering the x-direction gravity anomaly profile, assuming that the local gravity anomaly nearby region gravity anomaly is a partial arc segment of a curvature circle with a curvature radius of R and a curvature center of o, the x-direction cutting region gravity anomaly correction amount is:

wherein L is_ACCan be obtained according to the Pythagorean theorem

The radius of curvature R can be determined from the reciprocal of the curvature K

Replacing the first derivative with the first difference and ignoring the first small amount, then

Replacing the second derivative with a second order difference and ignoring the second order small quantity, then

Similarly, the gravity abnormal correction quantity delta Z (y) of the y-direction cutting area is obtained

ΔZ(x，y)＝(ΔZ(x)+ΔZ(y))/2；

Wherein L is_OCIs the distance from the curvature center 0 to the middle point C between the measuring point A and the measuring point B, L_ACIs the distance, X, from the measuring point A to the middle point C between the measuring point A and the measuring point B_AB、X_BEIs 2 times the cutting radius, Z_ABIs the difference between the gravity anomaly of the point A and the point B, Z (B) is the gravity anomaly of the point B, Z (A) is the gravity anomaly of the point A, and Z (E) is the pointG, gravity anomaly value of E.

7. A method of oil and gas exploration, comprising the steps of:

extracting local gravity anomaly data according to the method of any one of claims 1-6;

and carrying out oil-gas exploration according to the local gravity anomaly data.

Technical Field

The invention belongs to the technical field of gravity anomaly data processing and oil-gas exploration, and relates to a method for extracting local gravity anomaly based on measured gravity anomaly data.

Background

The commonly obtained gravity data is the Booth gravity anomaly data which is the simplest and most widely researched gravity anomaly data obtained in the current gravity research. The grid gravity anomaly is formed by superposing and combining regional gravity anomaly and local gravity anomaly. Wherein, the gravity anomaly caused by geological factors (such as regional geological structure) with larger distribution range and deeper influence is called regional gravity anomaly, and the gravity anomaly in small range caused by geological factors (such as mineral products and oil and gas) is called local gravity anomaly. Therefore, the extraction of the local gravity anomaly is not only an important research content of resource exploration, but also contributes to the research of geological structures. On the basis, how to accurately extract the local gravity anomaly becomes an important problem.

The interpolation cutting method for extracting local gravity anomaly is one of the commonly used methods in the prior art, and the interpolation operation of the current calculated point field value and the surrounding multipoint average value is used as a cutting operator, a cutting area field (area gravity anomaly) is obtained through continuous cutting, and then the cutting area field is subtracted from the position field anomaly, so that a cutting local field is obtained. According to the difference of the interpolation cutting operators, an interpolation cutting method based on a four-point circumference interpolation cutting operator, an eight-point window interpolation cutting operator and a dynamic improved interpolation cutting operator is formed in sequence.

However, the accuracy of the existing methods is yet to be further improved.

Therefore, there is a need to develop a method for extracting local gravity anomaly based on measured gravity anomaly data to solve one or more of the above-mentioned technical problems.

Disclosure of Invention

In order to solve at least one of the above technical problems, the applicant finds, through research, that the existing interpolation cutting method ignores different geological properties which cause abnormalities, and the obtained calculation result often generates certain errors. In contrast, the applicant improves the interpolation cutting method by combining four superposition relations of local anomaly (local gravity anomaly) and regional anomaly (regional gravity anomaly) of the gravity potential field from actual geological attributes, and further provides an interpolation cutting method based on the anomaly constraint condition, so that the accuracy of the obtained local gravity anomaly is improved.

According to an aspect of the present invention, there is provided a method for extracting local gravity anomaly based on measured gravity anomaly data, characterized by comprising the steps of:

measuring gravity anomaly data Z (x, y) by an aviation, aerospace, shipborne or ground platform, wherein the gravity anomaly data Z (x, y) consists of regional gravity anomaly R (x, y) and local gravity anomaly L (x, y);

determining a cutting radius r;

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 case of opposite bending direction, regional gravity anomaly R_n(x, y) is determined as:

delta Z (x, y) is the correction amount of gravity anomaly in cutting area, and R is used for gravity anomaly in area of 1 st cutting₁(x, y) represents, let R₁(x, y) as gravity anomaly data Z (x, y) of the next iteration, and determining the region gravity anomaly R of the 2 nd cutting according to the gravity anomaly data Z (x, y) of the iteration₂(x, y); successively iterate, finally have

So that the gravity in the region of the cutting radius r is abnormal

Calculating and extracting local gravity anomaly L (x, y) ═ Z (x, y) -R (x, y).

According to another aspect of the present invention, the calculation of the weighting parameter e is specifically;

e＝c+d(0≤e≤2)

(0. ltoreq. c.ltoreq.1; whenThen, c is 1)

(0. ltoreq. d. ltoreq.1; whenThen, d is 1)

ΔB_x＝Z(x+r，y)-Z(x-r，y)

ΔB_y＝Z(x，y+r)-Z(x，y-r)。

According to another aspect of the present invention, the comparing the curve bending direction of the regional gravity anomaly with the local gravity anomaly specifically includes: calculating a comparison parameter

B_x1＝Z(x+r，y+r)-[Z(x+2r，y+r)+Z(x，y+r)]/2

B_y1＝Z(x+r，y+r)-[Z(x+r，y+2r)+Z(x+r，y)]/2

B_x2＝Z(x+2r，y+r)-[Z(x+4r，y+r)+Z(x，y+r)]/2

B_y2＝Z(x+r，y+2r)-[Z(x+r，y+4r)+Z(x+r，y)]/2

When the comparison parameters simultaneously satisfy the following conditions, the curve bending directions indicating the regional gravity anomaly and the local gravity anomaly are opposite: b is_y1×B_x1＞0，B_y1×B_y2＜0，B_x1×B_x2＜0， B_y2×B_x2＞0。

According to still another aspect of the present invention, when the comparison parameter cannot satisfy the following condition at the same time, the curve bending directions indicating the regional gravity anomaly and the local gravity anomaly are the same: b is_y1×B_x1＞0， B_y1×B_y2＜0，B_x1×B_x2＜0，B_y2×B_x2＞0。

According to still another aspect of the present invention, the regional gravity anomaly is for the case where the bending directions are the same

According to still another aspect of the present invention, the cutting zone gravity abnormality correction amount Δ Z (x, y) is determined by:

considering the x-direction gravity anomaly profile, assuming that the local gravity anomaly nearby region gravity anomaly is a partial arc segment of a curvature circle with a curvature radius of R and a curvature center of o, the x-direction cutting region gravity anomaly correction amount is:

wherein L is_ACCan be obtained according to the Pythagorean theorem

The radius of curvature R can be determined from the reciprocal of the curvature K

Replacing the first derivative with the first difference and ignoring the first small amount, then

Replacing the second derivative with a second order difference and ignoring the second order small quantity, then

Similarly, the gravity abnormal correction quantity delta Z (y) of the y-direction cutting area is obtained

ΔZ(x，y)＝(ΔZ(x)+ΔZ(y))/2；

Wherein L is_OCIs the distance from the curvature center 0 to the middle point C between the measuring point A and the measuring point B, L_ACIs the distance, X, from the measuring point A to the middle point C between the measuring point A and the measuring point B_AB、X_BEIs 2 times the cutting radius, Z_ABThe difference value of the gravity anomaly of the measuring point A and the measuring point B is Z (B), the gravity anomaly value of the measuring point B is Z (A), the gravity anomaly value of the measuring point A is Z (E), and the gravity anomaly value of the measuring point E is Z (E).

According to still another aspect of the present invention, there is also provided a method of oil and gas exploration, characterized by comprising the steps of:

extracting local gravity anomaly data according to the method;

and carrying out oil-gas exploration according to the local gravity anomaly data.

The invention can obtain one or more of the following technical effects:

by starting from actual geological attributes and combining four superposition relations of local abnormity (local gravity abnormity) and regional abnormity (regional gravity abnormity) of the gravity potential field, the interpolation cutting method is improved, and the accuracy of the obtained local gravity abnormity is improved.

Drawings

The present invention will be described in further detail with reference to the accompanying drawings and specific embodiments.

FIG. 1 is a diagram illustrating the classification of local anomalies and regional anomalies in a superimposed relationship, in which a thick solid line represents a total field anomaly, a thin solid line between AB represents a regional anomaly, and dotted lines represent two end-point connections of the local anomalies, and (a) - (d) correspond to classification cases 1-4 in Table 1, respectively.

Fig. 2 is a diagram illustrating a cut-region abnormal correction amount (cut-region gravity abnormal correction amount) of a method for extracting a local gravity abnormality based on measured gravity abnormality data according to a preferred embodiment of the present invention.

Fig. 3 is a schematic diagram of a theoretical model of a method for extracting local gravity anomaly based on measured gravity anomaly data according to a preferred embodiment of the present invention.

FIG. 4 is a planar contour plot of forward results of the theoretical model of the present invention.

Fig. 5 is a cross-sectional view of the forward modeling result of the theoretical model of the present invention (y is 0).

Fig. 6 is a comparison plan view of the separation effect of the theoretical model gravity potential field.

Fig. 7 is a comparison cross-sectional view of the separation effect of the gravity potential field of the theoretical model of the present invention (y is 0).

FIG. 8 is a graph of gravity anomaly near the top of the upper mantle calculated and extracted by an interpolation cut method according to the prior art.

FIG. 9 is a graph of the near-top gravity anomaly of the upper mantle calculated and extracted according to the method of the present invention (constrained interpolation cutting).

Detailed Description

The best mode for carrying out the present invention will be described in detail with reference to the accompanying drawings, wherein the detailed description is for the purpose of illustrating the invention in detail, and is not to be construed as limiting the invention, as various changes and modifications can be made therein without departing from the spirit and scope thereof, which are intended to be encompassed within the appended claims.

Example 1

The applicant finds that the existing interpolation cutting method ignores different geological properties generating abnormity, and the obtained calculation result often generates certain errors. In contrast, the applicant improves the interpolation cutting method by combining four superposition relations of local anomaly (local gravity anomaly) and regional anomaly (regional gravity anomaly) of the gravity potential field from actual geological attributes, and further provides an interpolation cutting method based on the anomaly constraint condition, so that the accuracy of the obtained local gravity anomaly is improved.

According to a preferred embodiment of the present invention, there is provided a method for extracting local gravity anomaly based on measured gravity anomaly data, comprising the steps of:

measuring gravity anomaly data Z (x, y) by an aviation, aerospace, shipborne or ground platform, wherein the gravity anomaly data Z (x, y) consists of regional gravity anomaly R (x, y) and local gravity anomaly L (x, y);

determining a cutting radius r;

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 case of opposite bending direction, regional gravity anomaly R_n(x, y) is determined as:

delta Z (x, y) is the correction amount of gravity anomaly in cutting area, and R is used for gravity anomaly in area of 1 st cutting₁(x, y) represents, let R₁(x, y) as gravity anomaly data Z (x, y) of the next iteration, namely Z (x, y) ═ R₁(x, y) determining the region gravity anomaly R of the 2 nd cutting through the iterative gravity anomaly data Z (x, y)₂(x, y); successively iterate, finally have

So that the gravity in the region of the cutting radius r is abnormal

Calculating and extracting local gravity anomaly L (x, y) ═ Z (x, y) -R (x, y).

According to another preferred embodiment of the present invention, the calculation of the weighting parameter e is specifically;

e＝c+d(0≤e≤2)

(0. ltoreq. c.ltoreq.1; whenThen, c is 1)

(0. ltoreq. d. ltoreq.1; whenThen, d is 1)

ΔB_x＝Z(x，r，y)-Z(x-r，y)ΔB_y＝Z(x，y+r)-Z(x，y-r)。

According to another preferred embodiment of the present invention, the comparing the curve bending direction of the regional gravity anomaly with the local gravity anomaly specifically includes: calculating a comparison parameter

B_x1＝Z(x+r，y+r)-[Z(x+2r，y+r)+Z(x，y+r)]/2

B_y1＝Z(x+r，y+r)-[Z(x+r，y+2r)+Z(x+r，y)]/2

B_x2＝Z(x+2r，y+r)-[Z(x+4r，y+r)+Z(x，y+r)]/2

B_y2＝Z(x+r，y+2r)-[Z(x+r，y+4r)+Z(x，y+r)]/2

When the comparison parameters simultaneously satisfy the following conditions, the curve bending directions indicating the regional gravity anomaly and the local gravity anomaly are opposite: b is_y1×B_x1＞0，B_y1×B_x1＜0，B_y1×B_x1＜0， B_y2×B_x2＞0。

According to still another preferred embodiment of the present invention, when the comparison parameter fails to satisfy the following condition at the same time, the curve bending directions indicating the regional gravity anomaly and the local gravity anomaly are the same: b is_y1×B_x1＞0， B_y1×B_y2＜0，B_x1×B_x2＜0，R_y2×B_y2＞0。

According to still another preferred embodiment of the present invention, the regional gravity anomaly is given for the case where the bending directions are the same

According to another preferred embodiment of the present invention, the cutting region gravity abnormality correction amount Δ Z (x, y) is determined by:

considering the x-direction gravity anomaly profile, assuming that the local gravity anomaly nearby region gravity anomaly is a partial arc segment of a curvature circle with a curvature radius of R and a curvature center of o, the x-direction cutting region gravity anomaly correction amount is:

wherein L is_ACCan be obtained according to the Pythagorean theorem

The radius of curvature R can be determined from the reciprocal of the curvature K

Replacing the first derivative with the first difference and ignoring the first small amount, then

Replacing the second derivative with a second order difference and ignoring the second order small quantity, then

Similarly, the gravity abnormal correction quantity delta Z (y) of the y-direction cutting area is obtained

ΔZ(x，y)＝(ΔZ(x)+ΔZ(y))/2；

Wherein L is_OCIs the distance from the curvature center 0 to the middle point C between the measuring point A and the measuring point B, L_ACIs the distance, X, from the measuring point A to the middle point C between the measuring point A and the measuring point B_AB、X_BEIs 2 times the cutting radius, Z_ABThe difference value of the gravity anomaly of the measuring point A and the measuring point B is Z (B), the gravity anomaly value of the measuring point B is Z (A), the gravity anomaly value of the measuring point A is Z (E), and the gravity anomaly value of the measuring point E is Z (E).

There is also provided in accordance with still another preferred embodiment of the present invention a method of oil and gas exploration, characterized by the steps of:

extracting local gravity anomaly data according to the method of any one of claims 1-6;

and carrying out oil-gas exploration according to the local gravity anomaly data.

Example 2

On the basis of embodiment 1, this embodiment further describes the technical principle and practical application effect of the present invention in detail.

There is also provided in accordance with a preferred embodiment of the present invention a method for extracting local gravity anomaly based on measured gravity anomaly data, including the steps of:

(1) the cutting radius R is determined and the gravity anomaly Z (x, y) is measured, which consists of a field of area R (x, y) and a local field L (x, y), i.e. Z (x, y) ═ R (x, y) + L (x, y).

(2) Calculating the average A (x, y) of the gravity anomaly of some 4 points r away from the point (x, y)

(3) Calculating a weighting parameter e

e＝c+d(0≤e≤2)

(0. ltoreq. c.ltoreq.1; whenThen, c is 1)

(0. ltoreq. d. ltoreq.1; whenThen, d is 1)

ΔB_x＝Z(x+r，y)-Z(x-r，y)

ΔB_y＝Z(x，y+r)-Z(x，y-r)

(4) Comparing the bending direction of the curve of the regional abnormality with that of the local abnormality, the method can be specifically divided into 2 steps:

(4-1) calculation of comparative parameters

B_x1＝Z(x+r，y+r)-[Z(x+2r，y+r)+Z(x，y+r)]/2

B_y1＝Z(x+r，y+r)-[Z(x+r，y+2r)+Z(x+r，y)]/2

B_x2＝Z(x+2r，y+r)-[Z(x+4r，y+r)+Z(x，y+r)]/2

B_y2＝Z(x+r，y+2r)-[Z(x+r，y+4r)+Z(x+r，y)]/2

(4-2) when the comparison parameters fail to satisfy the following conditions at the same time, it is indicated that the regional abnormality and the local abnormal bending direction are the same (i.e., classification cases 1 and 2): b is_y1×B_x1＞0，B_y1×B_y2＜0， B_x1×B_x2＜0，B_y2×B_x2＞0。

(4-3) when the comparison parameters satisfy the following conditions at the same time, it is indicated that the regional abnormal and local abnormal curvatures are in opposite directions (i.e., categorizing cases 3 and 4): b is_y1×B_x1＞0，B_y1×B_y2＜0， B_x1×B_x2＜0，B_y2×B_x2＞0。

(5) For the same bend direction, the area field is a weighted average of Z (x, y) and A (x, y), and the cut area field is calculated using the following equation:

(6) for the case of the opposite bending direction, the area field is increased by the cutting area abnormality correction amount Δ Z (x, y) based on the weighted average of Z (x, y) and a (x, y), and this case can be subdivided into the following 3 steps to calculate the area abnormality:

(6-1) first, considering the x-direction abnormal section, assuming that the region abnormality near the local abnormality is a partial arc segment of a curvature circle having a curvature radius R and a curvature center o, the x-direction cutting region abnormality correction amount is:

wherein L is_ACCan be obtained according to the Pythagorean theorem

The radius of curvature R can be determined from the reciprocal of the curvature K

Replacing the first derivative with the first difference and ignoring the first small amount, then

Replacing the second derivative with a second order difference and ignoring the second order small quantity, then

(6-2) similarly obtaining abnormal correction quantity delta Z (y) of the y-direction cutting area

(6-3) taking Δ Z (x, y) ═ Δ Z (x) + Δ Z (y)/2

(6-4) after calculating the correction amount Δ Z (x, y), calculating the cutting region field using the following formula:

see FIG. 2, where L_OCIs the distance from the curvature center 0 to the middle point C between the measuring point A and the measuring point B, L_ACIs the distance, X, from the measuring point A to the middle point C between the measuring point A and the measuring point B_AB、X_BEIs 2 times the cutting radius, Z_ABThe difference value of the gravity anomaly of the measuring point A and the measuring point B is Z (B), the gravity anomaly value of the measuring point B is Z (A), the gravity anomaly value of the measuring point A is Z (E), and the gravity anomaly value of the measuring point E is Z (E).

Preferably, referring to fig. 2, the measuring point D is the current calculation point, the measuring point a is the backward measuring point with the interval cutting radius r, the measuring point B is the forward measuring point with the interval cutting radius r, and the measuring point E is the forward measuring point with the interval 3 times the cutting radius. That is, the measurement points A, D, B are spaced apart by a cutting radius r, B, E by a distance of 2 r. For example, measurement points A and D are separated by a distance r in the x-direction. As shown in fig. 2, the measurement point a can be used as a starting point, which is a starting point (first point) of the cross section.

(7) The area field obtained by the above method is called the area field of the 1 st cut, and R is used₁(x, y) represents, for R₁(x, y) repeating the steps (1) to (6) to obtain the region field R of the 2 nd cutting₂(x, y); successively iterate, finally have

Then there areThis area field is referred to as the cutting radius r area field.

(8) Calculating local anomalies, and subtracting the area field from the gravity anomaly Z (x, y) to obtain a local field:

L(x，y)＝Z(x，y)-R(x，y)。

preferably, the gravity anomaly Z (x, y) is acquired by different platforms such as aeronautical, aerospace, shipboard and terrestrial.

Referring to fig. 1, a schematic diagram of four cases of classifying the overlapping relationship of local anomaly (local gravity anomaly) and regional anomaly (regional gravity anomaly) is shown.

The applicant researches and discovers that the existing interpolation cutting method is only suitable for classification cases 1 and 2, namely, the area anomaly is between the value range of the total field anomaly value and the four-point circumference average value. However, in the case where the bending directions of the local anomaly curves are opposite to the bending directions of the local anomaly curves, namely, in the classification cases 3 and 4, although the calculation result can be infinitely approximated to the average value of the circumference of four points, the calculation result ignores the parts marked by the black vertical lines in (c) and (d) in fig. 1.

TABLE 1 local anomaly and regional anomaly overlap relation classification table

Note: the four classification cases correspond to (a), (b), (c), and (d) in fig. 1, respectively.

Preferably, for the case where the bending directions are opposite (i.e., sort cases 3 and 4), compensation is made by cutting area abnormality correction amount Δ Z (x, y), as shown in fig. 2.

Next, the theoretical model trial calculation will be further described.

First, theoretical model design is performed. Specifically, the theoretical model shown in fig. 3 is designed, the regional gravity generated by the three deep spheres is high, the local gravity generated by the shallow sphere 1 is low, the local gravity generated by the shallow spheres 2 and 3 is high, and the detailed parameters of the theoretical model are shown in table 2.

TABLE 2 spatial location and physical property parameter table of three-dimensional model anomaly

Next, theoretical model forward modeling was performed. The gravity anomaly generated by the theoretical model designed by the invention is the sum of derivatives of the gravitational potential of the residual mass of each sphere along the gravity direction (defined as the z direction), namely

In the above formula V_iThe gravity position generated by the sphere i, f is the gravity constant, and f is 6.672 x 10^-11m³/(kg·s²),m_iIs the remaining mass of the sphere i and,R_iis the radius of the sphere i, σ_iIs the residual density of the sphere i, r_iThe distance between the measuring point and the sphere center of the sphere i,

based on the formula (25) and the designed theoretical model parameters, forward results shown in fig. 4 and 5 can be obtained through forward calculation. As can be seen from FIGS. 4-5, anomaly # 1 belongs to Classification case 1, and anomalies # 2 and 3 belong to Classification case 3.

And further, performing separation and extraction effect comparison analysis on the theoretical model gravity potential field. The designed theoretical model forward modeling result is subjected to gravity potential field separation and extraction by respectively adopting an interpolation cutting method in the prior art and a method (a constrained interpolation cutting method) of the invention, a plurality of cutting radiuses r are adopted in the calculation process, and the r is determined to be 700m through screening and comparison. As shown in fig. 6, the gravity potential field separation plane effect is closer to the theoretical local anomaly (fig. 6d) than the local anomaly (fig. 6f) obtained by the visible constrained interpolation cutting method, and the local anomaly is largely submerged, so that the regional anomaly obtained by the two methods has a high similarity to the theoretical regional anomaly (fig. 6a, 6b and 6 c). In order to more intuitively embody the superiority of the method of the invention, a potential field separation curve (fig. 7) obtained by two methods of a main section (y is 0) is further drawn, and as the No. 1 abnormity is regional gravity high superposition local gravity high (classification situation 1), the two methods obtain consistent results, but aiming at the No. 2 abnormity and the No. 3 abnormity, the separation and extraction results obtained by the method of the invention are obviously improved.

Further, in order to quantitatively analyze the optimization effect of the constrained interpolation cutting method provided by the invention on the interpolation cutting method, an "improvement rate" (abbreviated as IR) is introduced to estimate the optimization quantity of the constrained interpolation cutting method on the local anomaly:

IR＝var(Δg_cz-Δg_ll)/var(Δg_ys-Δg_ll) (26)

in the formula,. DELTA.g_czIndicating local anomalies by interpolation, Δ g_ysRepresenting local anomalies, Δ g, by constrained interpolation_llFor theoretical local anomalies, var represents the mean square error.

For the plane areas with local anomalies of No. 2 and No. 3, the calculated improvement rate is 2.15, so that compared with an interpolation cutting method, the local anomalies extracted by the method are improved by 115%, and the optimization effect is obvious.

Furthermore, the gravity data of China and surrounding areas are used as test objects, and the gravity data range is a sector area with east longitude 60-141 degrees and north latitude 15-59 degrees. The gravity data come from an ultra-high-order Earth gravity field Model EGM2008(Earth gravity Model 2008), the order of the Model is completely 2159 times, the spatial resolution is about 5', the precision can reach 10mGal in the area of more than 80% in China, and the gravity data can be used for small-scale gravity drawing and structural research. The grid spacing for the bump gravity anomaly used in this example was 4 km.

The interpolation cutting method and the constraint interpolation cutting method provided by the invention are respectively utilized to separate the grid gravity anomaly, the cutting radius selects 18 and 19 times of grid spacing, the gravity anomaly value of a 72-76 km reference depth layer is obtained through calculation, and the gravity anomaly value is regarded as the gravity anomaly generated by the medium near the top of the upper mantle (fig. 8 and 9), wherein the AB section position is 43.76 degrees in north latitude.

Comparing the two gravity anomaly plan views (fig. 8(a) and 9(a)) obtained by the two methods, it can be seen that the overall morphology of the whole gravity anomaly is similar, but the amplitude and range of the positive gravity anomaly obtained by the invention (the constrained interpolation cutting method) are relatively reduced and even disappear, and the amplitude and range of the negative gravity anomaly obtained are relatively increased. The comparison results (such as fig. 8(b) and fig. 9(b)) on the AB section are more intuitive, the correction of the interpolation cutting method calculation result is realized, and more accurate gravity anomaly is obtained.

The invention can obtain one or more of the following technical effects:

by starting from actual geological attributes and combining four superposition relations of local abnormity (local gravity abnormity) and regional abnormity (regional gravity abnormity) of the gravity potential field, the interpolation cutting method is improved, and the accuracy of the obtained local gravity abnormity is improved.

It will be understood by those skilled in the art that the present invention is not limited to the embodiments described above, which are described in the specification and illustrated only to illustrate the principle of the present invention, but that various changes and modifications may be made therein without departing from the spirit and scope of the present invention, which fall within the scope of the invention as claimed. The scope of the invention is defined by the appended claims and equivalents thereof.