CN113761457A - Method for extracting local gravity anomaly based on measured gravity anomaly data - Google Patents

Method for extracting local gravity anomaly based on measured gravity anomaly data 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
Chinese (zh)
Other versions
CN113761457B (en
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/en
Publication of CN113761457A publication Critical patent/CN113761457A/en
Application granted granted Critical
Publication of CN113761457B publication Critical patent/CN113761457B/en
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

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.

Description

Method for extracting local gravity anomaly based on measured 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 Rn(x, y) is determined as:
Figure BDA0003254369000000021
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 cutting1(x, y) represents, let R1(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 iteration2(x, y); successively iterate, finally have
Figure BDA0003254369000000022
So that the gravity in the region of the cutting radius r is abnormal
Figure BDA0003254369000000023
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)
Figure BDA0003254369000000031
(0. ltoreq. c.ltoreq.1; when
Figure BDA0003254369000000032
Then, c is 1)
Figure BDA0003254369000000033
(0. ltoreq. d. ltoreq.1; when
Figure BDA0003254369000000034
Then, d is 1)
Figure BDA0003254369000000035
Figure BDA0003254369000000036
ΔBx＝Z(x+r，y)-Z(x-r，y)
ΔBy＝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
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
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 isy1×Bx1＞0，By1×By2＜0，Bx1×Bx2＜0， By2×Bx2＞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 isy1×Bx1＞0， By1×By2＜0，Bx1×Bx2＜0，By2×Bx2＞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
Figure BDA0003254369000000041
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:
Figure BDA0003254369000000042
wherein L isACCan be obtained according to the Pythagorean theorem
Figure BDA0003254369000000043
The radius of curvature R can be determined from the reciprocal of the curvature K
Figure BDA0003254369000000044
Replacing the first derivative with the first difference and ignoring the first small amount, then
Figure BDA0003254369000000045
Replacing the second derivative with a second order difference and ignoring the second order small quantity, then
Figure BDA0003254369000000046
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 isOCIs the distance from the curvature center 0 to the middle point C between the measuring point A and the measuring point B, LACIs the distance, X, from the measuring point A to the middle point C between the measuring point A and the measuring point BAB、XBEIs 2 times the cutting radius, ZABThe 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 Rn(x, y) is determined as:
Figure BDA0003254369000000071
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 cutting1(x, y) represents, let R1(x, y) as gravity anomaly data Z (x, y) of the next iteration, namely Z (x, y) ═ R1(x, y) determining the region gravity anomaly R of the 2 nd cutting through the iterative gravity anomaly data Z (x, y)2(x, y); successively iterate, finally have
Figure BDA0003254369000000072
So that the gravity in the region of the cutting radius r is abnormal
Figure BDA0003254369000000073
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)
Figure BDA0003254369000000074
(0. ltoreq. c.ltoreq.1; when
Figure BDA0003254369000000075
Then, c is 1)
Figure BDA0003254369000000076
(0. ltoreq. d. ltoreq.1; when
Figure BDA0003254369000000077
Then, d is 1)
Figure BDA0003254369000000078
Figure BDA0003254369000000079
ΔBx＝Z(x，r，y)-Z(x-r，y)ΔBy＝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
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
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 isy1×Bx1＞0，By1×Bx1＜0，By1×Bx1＜0， By2×Bx2＞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 isy1×Bx1＞0， By1×By2＜0，Bx1×Bx2＜0，Ry2×By2＞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
Figure BDA0003254369000000081
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:
Figure BDA0003254369000000091
wherein L isACCan be obtained according to the Pythagorean theorem
Figure BDA0003254369000000092
The radius of curvature R can be determined from the reciprocal of the curvature K
Figure BDA0003254369000000093
Replacing the first derivative with the first difference and ignoring the first small amount, then
Figure BDA0003254369000000094
Replacing the second derivative with a second order difference and ignoring the second order small quantity, then
Figure BDA0003254369000000095
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 isOCIs the distance from the curvature center 0 to the middle point C between the measuring point A and the measuring point B, LACIs the distance, X, from the measuring point A to the middle point C between the measuring point A and the measuring point BAB、XBEIs 2 times the cutting radius, ZABThe 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)
Figure BDA0003254369000000101
(3) Calculating a weighting parameter e
e＝c+d(0≤e≤2)
Figure BDA0003254369000000102
(0. ltoreq. c.ltoreq.1; when
Figure BDA0003254369000000103
Then, c is 1)
Figure BDA0003254369000000104
(0. ltoreq. d. ltoreq.1; when
Figure BDA0003254369000000105
Then, d is 1)
Figure BDA0003254369000000106
Figure BDA0003254369000000111
ΔBx＝Z(x+r，y)-Z(x-r，y)
ΔBy＝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
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) 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 isy1×Bx1＞0，By1×By2＜0， Bx1×Bx2＜0，By2×Bx2＞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 isy1×Bx1＞0，By1×By2＜0， Bx1×Bx2＜0，By2×Bx2＞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:
Figure BDA0003254369000000112
(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:
Figure BDA0003254369000000121
wherein L isACCan be obtained according to the Pythagorean theorem
Figure BDA0003254369000000122
The radius of curvature R can be determined from the reciprocal of the curvature K
Figure BDA0003254369000000123
Replacing the first derivative with the first difference and ignoring the first small amount, then
Figure BDA0003254369000000124
Replacing the second derivative with a second order difference and ignoring the second order small quantity, then
Figure BDA0003254369000000125
(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:
Figure BDA0003254369000000131
see FIG. 2, where LOCIs the distance from the curvature center 0 to the middle point C between the measuring point A and the measuring point B, LACIs the distance, X, from the measuring point A to the middle point C between the measuring point A and the measuring point BAB、XBEIs 2 times the cutting radius, ZABThe 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 used1(x, y) represents, for R1(x, y) repeating the steps (1) to (6) to obtain the region field R of the 2 nd cutting2(x, y); successively iterate, finally have
Figure BDA0003254369000000132
Then there are
Figure BDA0003254369000000133
This 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
Figure BDA0003254369000000141
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
Figure BDA0003254369000000151
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
Figure BDA0003254369000000152
In the above formula ViThe gravity position generated by the sphere i, f is the gravity constant, and f is 6.672 x 10-11m3/(kg·s2),miIs the remaining mass of the sphere i and,
Figure BDA0003254369000000153
Riis the radius of the sphere i, σiIs the residual density of the sphere i, riThe distance between the measuring point and the sphere center of the sphere i,
Figure BDA0003254369000000154
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(Δgcz-Δgll)/var(Δgys-Δgll) (26)
in the formula,. DELTA.gczIndicating local anomalies by interpolation, Δ gysRepresenting local anomalies, Δ g, by constrained interpolationllFor 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.

Claims (7)

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 Rn(x, y) is determined as:
Figure FDA0003254368990000011
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 cutting1(x, y) represents, let R1(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 iteration2(x, y); successively iterate, finally have
Figure FDA0003254368990000012
So that the gravity in the region of the cutting radius r is abnormal
Figure FDA0003254368990000013
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)
Figure FDA0003254368990000021
(0. ltoreq. c.ltoreq.1; when
Figure FDA0003254368990000022
Then, c is 1)
Figure FDA0003254368990000023
(0. ltoreq. d. ltoreq.1; when
Figure FDA0003254368990000024
Then, d is 1)
Figure FDA0003254368990000025
Figure FDA0003254368990000026
ΔBx＝Z(x+r，y)-Z(x-r，y)
ΔBy＝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
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
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 isy1×Bx1＞0，By1×By2＜0，Bx1×Bx2＜0，By2×Bx2＞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 isy1×Bx1＞0，By1×By2＜0，Bx1×Bx2＜0，By2×Bx2＞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
Figure FDA0003254368990000031
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:
Figure FDA0003254368990000032
wherein L isACCan be obtained according to the Pythagorean theorem
Figure FDA0003254368990000041
The radius of curvature R can be determined from the reciprocal of the curvature K
Figure FDA0003254368990000042
Replacing the first derivative with the first difference and ignoring the first small amount, then
Figure FDA0003254368990000043
Replacing the second derivative with a second order difference and ignoring the second order small quantity, then
Figure FDA0003254368990000044
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 isOCIs the distance from the curvature center 0 to the middle point C between the measuring point A and the measuring point B, LACIs the distance, X, from the measuring point A to the middle point C between the measuring point A and the measuring point BAB、XBEIs 2 times the cutting radius, ZABIs 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.
CN202111055317.1A 2021-09-09 2021-09-09 Method for extracting local gravity anomaly based on measured gravity anomaly data Active CN113761457B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111055317.1A CN113761457B (en) 2021-09-09 2021-09-09 Method for extracting local gravity anomaly based on measured gravity anomaly data

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111055317.1A CN113761457B (en) 2021-09-09 2021-09-09 Method for extracting local gravity anomaly based on measured gravity anomaly data

Publications (2)

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

Family

ID=78794292

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111055317.1A Active CN113761457B (en) 2021-09-09 2021-09-09 Method for extracting local gravity anomaly based on measured gravity anomaly data

Country Status (1)

Country Link
CN (1) CN113761457B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115629426A (en) * 2022-10-24 2023-01-20 中国自然资源航空物探遥感中心 Method, system, device and medium for extracting local gravity anomaly

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2007146867A (en) * 2007-12-17 2009-06-27 Горный институт Уральского отделения Российской академии наук (ГИ УрО РАН) (RU) METHOD FOR MULTICOMPONENT GRAVIMETRIC MODELING OF THE GEOLOGICAL ENVIRONMENT
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 (en) * 2018-12-11 2019-04-02 中国人民解放军火箭军工程大学 Gravity datum figure time-varying modification method and system based on gravity anomaly time-varying
CN110031000A (en) * 2019-05-21 2019-07-19 北京理工大学 A kind of evaluation method of Method in Gravity Aided INS region suitability
CN110706275A (en) * 2019-10-16 2020-01-17 中国石油大学(华东) Local gravity anomaly extraction method based on satellite altimetry gravity data
CN111337993A (en) * 2020-03-30 2020-06-26 中国地质科学院 Variable density and variable depth constraint-based gravity density interface inversion method
CN112464520A (en) * 2020-10-28 2021-03-09 中国石油天然气集团有限公司 Local gravity anomaly depth inversion method and device
CN112596106A (en) * 2020-11-09 2021-04-02 中国人民武装警察部队黄金第五支队 Method for gravity-seismic joint inversion of density interface distribution in spherical coordinate system

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2007146867A (en) * 2007-12-17 2009-06-27 Горный институт Уральского отделения Российской академии наук (ГИ УрО РАН) (RU) METHOD FOR MULTICOMPONENT GRAVIMETRIC MODELING OF THE GEOLOGICAL ENVIRONMENT
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 (en) * 2018-12-11 2019-04-02 中国人民解放军火箭军工程大学 Gravity datum figure time-varying modification method and system based on gravity anomaly time-varying
CN110031000A (en) * 2019-05-21 2019-07-19 北京理工大学 A kind of evaluation method of Method in Gravity Aided INS region suitability
CN110706275A (en) * 2019-10-16 2020-01-17 中国石油大学(华东) Local gravity anomaly extraction method based on satellite altimetry gravity data
CN111337993A (en) * 2020-03-30 2020-06-26 中国地质科学院 Variable density and variable depth constraint-based gravity density interface inversion method
CN112464520A (en) * 2020-10-28 2021-03-09 中国石油天然气集团有限公司 Local gravity anomaly depth inversion method and device
CN112596106A (en) * 2020-11-09 2021-04-02 中国人民武装警察部队黄金第五支队 Method for gravity-seismic joint inversion of density interface distribution in spherical coordinate system

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 (en) * 2022-10-24 2023-01-20 中国自然资源航空物探遥感中心 Method, system, device and medium for extracting local gravity anomaly

Also Published As

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

Similar Documents

Publication Publication Date Title
CN108489402B (en) Rapid and fine dereferencing method for surface mine slope rock mass joint scale based on three-dimensional laser scanning
CN110058317B (en) Aviation transient electromagnetic data and aviation magnetotelluric data joint inversion method
CN111323776B (en) Method for monitoring deformation of mining area
CN107227950B (en) Method for evaluating integrity of actual drilling hole track
Piccini Recent developments on morphometric analysis of karst caves
CN106338277B (en) A kind of building change detecting method based on baseline
Jouini et al. Multifractal analysis of reservoir rock samples using 3D X-ray micro computed tomography images
CN113761457A (en) Method for extracting local gravity anomaly based on measured gravity anomaly data
CA2402887A1 (en) Method for characterization of multi-scale geometric attributes
CN110414060A (en) A kind of potential field Boundary Recognition method based on quadravalence spectral moment
Wang et al. Fragmentation calculation method for blast muck piles in open-pit copper mines based on three-dimensional laser point cloud data
CN106873031B (en) A kind of 3 D seismic observation system vertical resolution quantitative analysis evaluation method
CN103149130A (en) Analytical method for particle size in conglomerate core particle structure
CN112785596B (en) Dot cloud picture bolt segmentation and height measurement method based on DBSCAN clustering
CN110954958A (en) Crack and fault prediction method and system
CN112596113A (en) Method for identifying field source position based on intersection points of characteristic values of different gradients of gravity
CN109374436B (en) Hydrate deposit shear band identification method based on CT image
CN1869734A (en) Method for processing varying density terrain correction by heavy prospecting data
CN112036424A (en) Submarine landslide hazard analysis method based on unsupervised machine learning
CN113779888B (en) Ground subsidence risk assessment method, device, equipment and storage medium
CN110703347A (en) Gravity fracture image identification method based on construction background
CN107239629B (en) Fractal dimension analysis method for determining reasonable size of rock structural plane laboratory
CN112859073B (en) Road damage assessment method based on PSInSAR technology
CN117152472B (en) Building deformation measurement method for building design
CN113325474B (en) Method for discriminating biological reef

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