CN111950108A - Two-degree variable density body gravity gradient tensor calculation method - Google Patents

Two-degree variable density body gravity gradient tensor calculation method Download PDF

Info

Publication number
CN111950108A
CN111950108A CN201910398005.7A CN201910398005A CN111950108A CN 111950108 A CN111950108 A CN 111950108A CN 201910398005 A CN201910398005 A CN 201910398005A CN 111950108 A CN111950108 A CN 111950108A
Authority
CN
China
Prior art keywords
coordinates
vertex
order
observation point
density
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
CN201910398005.7A
Other languages
Chinese (zh)
Other versions
CN111950108B (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.)
Ocean University of China
Original Assignee
Ocean University of China
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 Ocean University of China filed Critical Ocean University of China
Priority to CN201910398005.7A priority Critical patent/CN111950108B/en
Publication of CN111950108A publication Critical patent/CN111950108A/en
Application granted granted Critical
Publication of CN111950108B publication Critical patent/CN111950108B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Complex Calculations (AREA)

Abstract

The invention relates to a two-degree variable density body weight gradient tensor calculation method, and belongs to the technical field of geophysical exploration. The invention comprises the following steps: establishing a plane rectangular coordinate system according to a research area; approximating the section of the second-degree body to a polygon and determining the coordinates of a vertex; simulating density variation of a binary volume by using a binary polynomial function; setting an observation point along the observation measuring line and determining the coordinate of the observation point; appointing a current observation point, calculating the contribution value of the gravity gradient tensor corresponding to each edge of the polygon on the current observation point, and adding the contribution values to obtain the value of the gravity gradient tensor on the current observation point; and repeating the calculation steps until the values of the gravity gradient tensor at all the observation points are output. The method uses an analytic formula derived in a spatial domain, and can calculate to obtain an accurate analytic solution of the gravity gradient tensor. Compared with the constant density unit accumulation calculation method, the method has higher efficiency under the same precision requirement.

Description

Two-degree variable density body gravity gradient tensor calculation method
Technical Field
The invention relates to a calculation method of a two-degree variable density body weight gradient tensor, and belongs to the technical field of geophysical exploration.
Background
A geobody can be considered as a two-degree volume if its dimension in the direction of travel is much larger than the dimension perpendicular to its direction of travel. When the density distribution of the geologic body is relatively complex, the geologic body is generally divided into a plurality of simple regular body units with uniform density, the gravity response generated by each regular body unit is calculated by using the existing formula, and then the gravity responses generated by the geologic body can be obtained by accumulation. The accuracy of the gravity response data calculated in this way depends on the number of divided cells. With the increase of the number of the division units, the calculation precision is higher and higher, but at the same time, the calculation efficiency is also continuously reduced.
In addition to the above-described method using the gravity response calculation formula for the constant density body, a gravity response calculation formula for a variable density body has been developed. The gravity anomaly calculation formula which expresses the density change by utilizing a polynomial function of a space coordinate and approximates a two-dimensional shape by a polygon is deduced, and a corresponding analytical calculation method of the gravity anomaly is formed.
The gravity gradient data has a better resolution than the gravity anomaly data. With the advent of gravity gradiometers, people are increasingly studying underground geological structures and looking for mineral resources by observing gravity gradient data. The method for calculating the gravity gradient tensor generated by the variable density body is a basic technology for processing and interpreting gravity gradient data. At present, a fourier domain-based method is used for calculating the gravity gradient tensor of the variable density polygon, but the calculation accuracy is not high due to the truncation error. Therefore, the analytic calculation method of the two-degree variable density polygonal body gravity gradient tensor is developed, and the calculation accuracy of the two-degree variable density polygonal body gravity gradient tensor can be obviously improved.
Disclosure of Invention
Aiming at the defects in the prior art, the invention provides a calculation method of a two-degree variable density body weight gradient tensor, which can obtain an accurate analytic solution by deducing a calculation formula of the gravity gradient tensor caused by a variable density polygonal model in a spatial domain. Wherein the density of the polygonal model is represented by a bivariate polynomial function of arbitrary order of horizontal and vertical coordinates.
The invention is realized by adopting the following technical scheme:
the method comprises the following steps: a plane rectangular coordinate system is established in the plane of the section of the two-degree body, and a transverse axisxIndicating a horizontal distance, taking the right as positive; longitudinal axiszIndicating depth, the orientation is positive;
step two: the section of the two-degree body is approximated to be one withN k The polygon of the side is started from any vertex, and the coordinates of each vertex are determined to be (in turn) along the counterclockwise directionx k , z k ) Whereink = 1, 2, …, N k , N k + 1；
Step three: using a function of binary polynomials
Figure DEST_PATH_IMAGE001
Representing the change in density of the second degree volume, determiningxHighest order of coordinatesN x Andzhighest order of coordinatesN z And coefficients of the polynomialD i, j ；
Step four: arranging a plurality of observation points along the observation measuring line and determining coordinates of the observation points, wherein the abscissa of each observation point cannot be equal to the abscissa of the vertex of the polygon;
step five: specifying a current observation point having coordinates ofx 0, z 0)，x 0 ≠ x k ；
Step six: respectively calculating the number of polygons on the current observation pointk(k = 1, 2, …, N k ) Contribution U of gravity gradient tensor U corresponding to edges k Wherein
Figure DEST_PATH_IMAGE002
，
Figure DEST_PATH_IMAGE003
；
Step seven: will be provided withN k Adding the contribution values corresponding to the edges to obtain the value of the gravity gradient tensor at the current observation point, wherein,
U xz (x 0,z 0)=u xz_1+u xz_2+…+
Figure DEST_PATH_IMAGE004
，
U zx (x 0,z 0)=u zx_1+u zx_2+…+
Figure DEST_PATH_IMAGE005
，
U xx (x 0,z 0)=u xx_1+u xx_2+…+
Figure DEST_PATH_IMAGE006
，
U zz (x 0,z 0)=u zz_1+u zz_2+…+
Figure DEST_PATH_IMAGE007
；
step eight: judging whether observation points which are not calculated exist, if so, returning to the step five, designating the observation points as current observation points, and performing new calculation; if not, executing the next step;
step nine: outputting the gravity gradient tensor elements on all observation pointsU xz 、U zx 、U xx 、U zz The calculated value of (a).
Further, the sixth step is to calculatekThe formula of the gravity gradient tensor contribution value corresponding to the edges is as follows:
is as followskEdge andzwhen the axes are parallel to each other, the axes are parallel,
Figure DEST_PATH_IMAGE008
（1）
Figure DEST_PATH_IMAGE009
（2）
in the above-mentioned formulae (1) to (2),Gis a constant of universal gravitation, i.e.G = 6.67408×10-11N×m2/kg2；N x AndN z is a function of densityxCoordinates andzthe highest order of the coordinates is the highest order,iandjis a non-zero termD i, j x i z j InxCoordinates andzthe order to which the coordinates correspond to,D i, j is a coefficient of a density polynomial;nandmis also thatxCoordinates andzorder of coordinates:nis between 0 andiin the first order in between, and in the second order,mis between 0 andjan order in between;C i n andC j m is the number of combinations; (x 0, z 0) Is the coordinate of the current observation point; function(s)E xz AndE zz is represented byi、j、n、mTwo piecewise functions are involved:
Figure DEST_PATH_IMAGE010
（3）
Figure DEST_PATH_IMAGE011
（4）
in the above formulae (3) to (4), (b)x k , z k ) Is the firstkCoordinates of a vertex, (x k + 1, z k + 1) Is the firstk + 1 vertex coordinates, andr k andr k + 1respectively representkA vertex and the secondk+ 1 distance between the vertex and the current observation point, anr k + 1=[(x k + 1 - x 0)2 + (z k + 1 - z 0)2]1/2，r k =[(x k - x 0)2 + (z k - z 0)2]1/2；HRepresents a recursive sequence, if usedlRepresenting items in this columnThe sequence number of (2), then recursion the sequence numberH l The general formula of (A) is as follows:
Figure DEST_PATH_IMAGE012
（5）
in the above formula (5)KRepresents another recursion sequence, if usedsIndicating the sequence number of the item in the sequence, the sequence is recurredK s The general formula of (A) is as follows:
Figure DEST_PATH_IMAGE013
（6）
when it is at firstkEdge andzthe axes are not parallel, and the current point of observation is the secondkWhen the edges are on the same straight line,
Figure DEST_PATH_IMAGE014
（7）
Figure DEST_PATH_IMAGE015
（8）
in the above-mentioned formulae (7) to (8),Gis a constant of universal gravitation, i.e.G = 6.67408×10-11 N×m2/kg2；N x AndN z is a function of densityxCoordinates andzthe highest order of the coordinates is the highest order,iandjis a non-zero termD i, j x i z j InxCoordinates andzthe order to which the coordinates correspond to,D i, j is a coefficient of a density polynomial; (x 0, z 0) Is the coordinate of the current observation point; and (a)x k , z k ) Is the firstkCoordinates of a vertex, (x k + 1, z k + 1) Is the firstk The coordinates of the + 1 vertex points,pandqare respectively the firstkThe slope and intercept of the straight line on which the edges lie,p =(z k + 1 - z k )/(x k + 1 - x k ),q =(z k x k + 1 - z k + 1 x k )/(x k + 1 - x k ) (ii) a And arrays ofMThe calculation formula of (2) is as follows:
Figure DEST_PATH_IMAGE016
（9）
in the above formula (9)c =1 + p 2, d = [p(q - z 0) - x 0]/c；
Third whenkEdge andzthe axes are not parallel, and the current point of observation is the secondkWhen the edges are not in the same straight line,
Figure DEST_PATH_IMAGE017
（10）
Figure DEST_PATH_IMAGE018
（11）
in the above-mentioned formulae (10) to (11),Gis a constant of universal gravitation, i.e.G = 6.67408×10-11 N×m2/kg2；N x AndN z is a function of densityxCoordinates andzthe highest order of the coordinates is the highest order,iandjis in a density polynomialxAndzthe order of the corresponding order is set as,D i, j is a coefficient of a density polynomial;nandmis also thatxCoordinates andzorder of coordinates:nis between 0 andiin the first order in between, and in the second order,mis between 0 andjan order in between;C i n andC j m is the number of combinations; (x 0, z 0) Is the coordinate of the current observation point; function(s)F xz AndF zz is represented byi、j、n、mTwo piecewise functions are involved:
Figure DEST_PATH_IMAGE019
（12）
Figure DEST_PATH_IMAGE020
（13）
in the above formulae (12) to (13), (b)x k , z k ) Is the firstkCoordinates of a vertex, (x k + 1, z k + 1) Is the firstk The coordinates of the + 1 vertex points,pandqare respectively the firstkThe slope and intercept of the straight line on which the edges lie,p =(z k + 1 - z k )/(x k + 1 - x k ), q =(z k x k + 1 - z k + 1 x k )/(x k + 1 - x k )；Q = px 0 + q - z 0,；l 1andl 2also represents the order, but the value ranges of the two can be different due to different formulas,
Figure DEST_PATH_IMAGE021
、
Figure DEST_PATH_IMAGE022
、
Figure DEST_PATH_IMAGE023
、
Figure DEST_PATH_IMAGE024
、
Figure DEST_PATH_IMAGE025
and
Figure DEST_PATH_IMAGE026
are all combination numbers;IandJrespectively represent two number sequences, wherein the number sequencesIThe calculation formula of (2) is as follows:
Figure DEST_PATH_IMAGE027
（14）
in the above formula (14)Q = px 0 + q - z 0, b = p(q - z 0) - x 0, c =1 + p 2，r k Andr k + 1respectively representkA vertex and the secondkThe distance between + 1 vertex and the current observation point; another array of numbersJIs a recursive sequence, if usedlIndicating the sequence number of the item in the sequence, the sequence is recurredJ l The general formula of (A) is as follows:
Figure DEST_PATH_IMAGE028
（15）
in the above formula (15)Q = px 0 + q - z 0, a =(q - z 0)2 + x 0 2, b = p(q - z 0) - x 0, c =1 + p 2，r k Andr k + 1respectively representkA vertex and the secondkThe distance between + 1 vertex and the current observation point;I 0the formula (2) is shown as (14).
The invention has the beneficial effects that: by adopting the two-degree variable density body weight gradient tensor calculation method, a polygon with a binary polynomial density function is used, and the two-degree geologic body with complex underground shape and variable density in actual exploration can be flexibly approximated. This is because any two-dimensional density distribution can be well approximated by a sufficiently high-order bivariate polynomial, and any complex boundary can be well approximated by an arbitrary polygon with a sufficiently large number of edges. In addition, the calculation method can obtain the analytic solution of the gravity gradient tensor. The calculation method of the invention is simple, accurate, efficient and strong in adaptability, and the inversion algorithm of gravity gradient data is easy to establish based on the gravity gradient formula of the invention.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a geometric relationship diagram of observation points and a dyadic volume of the present invention;
FIG. 3 is a variable density hexagonal model simulating a valley;
FIG. 4 is a comparison of the results of the gravity gradient tensor generated by the model of FIG. 3 calculated using the present method and a conventional constant density volume stacking method, respectively;
fig. 5 shows the relation between the root mean square error and the number of divided units in the constant density volume superposition method relative to the present method, and the CPU time consumed by each of the two methods.
Detailed Description
In order to make the aforementioned and other objects, features and advantages of the present invention comprehensible, embodiments accompanied with figures are described in detail below. The flow chart of the present invention, as shown in fig. 1, includes the following steps:
the method comprises the following steps: a plane rectangular coordinate system is established in the plane of the section of the two-degree body, and a transverse axisxIndicating a horizontal distance, taking the right as positive; longitudinal axiszIndicating depth, the orientation is positive;
step two: the section of the two-degree body is approximated to be one withN k The polygon of the side is started from any vertex, and the coordinates of each vertex are determined to be (in turn) along the counterclockwise directionx k , z k ) Whereink = 1, 2, …, N k , N k + 1；
Step three: using a function of binary polynomials
Figure 289120DEST_PATH_IMAGE001
Indicating the change in density of a two-dimensional bodyDeterminingxHighest order of coordinatesN x Andzhighest order of coordinatesN z And coefficients of the polynomialD i, j ；
Step four: arranging a plurality of observation points along the observation measuring line and determining coordinates of the observation points, wherein the abscissa of each observation point cannot be equal to the abscissa of the vertex of the polygon;
step five: specifying a current observation point having coordinates ofx 0, z 0)，x 0 ≠ x k ；
Step six: respectively calculating the number of polygons on the current observation pointk(k = 1, 2, …, N k ) Contribution U of gravity gradient tensor U corresponding to edges k Wherein
Figure DEST_PATH_IMAGE029
，
Figure DEST_PATH_IMAGE030
；
Step seven: will be provided withN k Adding the contribution values corresponding to the edges to obtain the value of the gravity gradient tensor at the current observation point, wherein,
U xz (x 0,z 0)=u xz_1+u xz_2+…+
Figure 906395DEST_PATH_IMAGE004
，
U zx (x 0,z 0)=u zx_1+u zx_2+…+
Figure 623815DEST_PATH_IMAGE005
，
U xx (x 0,z 0)=u xx_1+u xx_2+…+
Figure 19024DEST_PATH_IMAGE006
，
U zz (x 0,z 0)=u zz_1+u zz_2+…+
Figure 833396DEST_PATH_IMAGE007
；
step eight: judging whether observation points which are not calculated exist, if so, returning to the step five, designating the observation points as current observation points, and performing new calculation; if not, executing the next step;
step nine: outputting the gravity gradient tensor elements on all observation pointsU xz 、U zx 、U xx 、U zz The calculated value of (a).
In the sixth step, the first stepkThe specific way of the gravity gradient tensor contribution value corresponding to the edges is as follows:
first, it is judgedkEdge andzwhether the axes are parallel, ifx k = x k + 1If so, the parallel is indicated; if it isx k ≠ x k + 1If so, it is indicated as not parallel;
then, the current observation point is judged (x 0, z 0) And a firstkWhether the edges are on the same straight line or not is calculated firstkSlope of edgep =(z k + 1 - z k )/(x k + 1 - x k ) And interceptq =(z k x k + 1 - z k + 1 x k )/(x k + 1 - x k ) Then calculates the discriminantQ = px 0 + q - z 0A value of, ifQ= 0, explains bothOn the same straight line; if it isQIf not equal to 0, the two are not on the same straight line;
finally, different calculation formulas are selected according to the judgment result to calculateu xz_k 、u zx_k 、u xx_k 、u zz_k The specific calculation formula includes:
is as followskEdge andzwhen the axes are parallel to each other, the axes are parallel,
Figure 655859DEST_PATH_IMAGE008
（1）
Figure 109843DEST_PATH_IMAGE009
（2）
in the above-mentioned formulae (1) to (2),Gis a constant of universal gravitation, i.e.G = 6.67408×10-11N×m2/kg2；N x AndN z is a function of densityxCoordinates andzthe highest order of the coordinates is the highest order,iandjis a non-zero termD i, j x i z j InxCoordinates andzthe order to which the coordinates correspond to,D i, j is a coefficient of a density polynomial;nandmis also thatxCoordinates andzorder of coordinates:nis between 0 andiin the first order in between, and in the second order,mis between 0 andjan order in between;C i n andC j m is the number of combinations; (x 0, z 0) Is the coordinate of the current observation point; function(s)E xz AndE zz is represented byi、j、n、mTwo piecewise functions are involved:
Figure DEST_PATH_IMAGE031
（3）
Figure DEST_PATH_IMAGE032
（4）
in the above formulae (3) to (4), (b)x k , z k ) Is the firstkCoordinates of a vertex, (x k + 1, z k + 1) Is the firstk + 1 vertex coordinates, andr k andr k + 1respectively representkA vertex and the secondk+ 1 distance between the vertex and the current observation point, anr k + 1=[(x k + 1 - x 0)2 + (z k + 1 - z 0)2]1/2，r k =[(x k - x 0)2 + (z k - z 0)2]1/2；HRepresents a recursive sequence, if usedlIndicating the sequence number of the item in the sequence, the sequence is recurredH l The general formula of (A) is as follows:
Figure 308743DEST_PATH_IMAGE012
（5）
in the above formula (5)KRepresents another recursion sequence, if usedsIndicating the sequence number of the item in the sequence, the sequence is recurredK s The general formula of (A) is as follows:
Figure 977622DEST_PATH_IMAGE013
（6）
when it is at firstkEdge andzthe axes are not parallel, and the current point of observation is the secondkWhen the edges are on the same straight line,
Figure 908669DEST_PATH_IMAGE014
（7）
Figure 351414DEST_PATH_IMAGE015
（8）
in the above-mentioned formulae (7) to (8),Gis a constant of universal gravitation, i.e.G = 6.67408×10-11 N×m2/kg2；N x AndN z is a function of densityxCoordinates andzthe highest order of the coordinates is the highest order,iandjis a non-zero termD i, j x i z j InxCoordinates andzthe order to which the coordinates correspond to,D i, j is a coefficient of a density polynomial; (x 0, z 0) Is the coordinate of the current observation point; and (a)x k , z k ) Is the firstkCoordinates of a vertex, (x k + 1, z k + 1) Is the firstk The coordinates of the + 1 vertex points,pandqare respectively the firstkThe slope and intercept of the straight line on which the edges lie,p =(z k + 1 - z k )/(x k + 1 - x k ),q =(z k x k + 1 - z k + 1 x k )/(x k + 1 - x k ) (ii) a And arrays ofMThe calculation formula of (2) is as follows:
Figure 88425DEST_PATH_IMAGE016
（9）
in the above formula (9)c =1 + p 2, d = [p(q - z 0) - x 0]/c；
Third whenkEdge andzthe axes are not parallel, and the current point of observation is the secondkWhen the edges are not in the same straight line,
Figure DEST_PATH_IMAGE033
（10）
Figure DEST_PATH_IMAGE034
（11）
in the above-mentioned formulae (10) to (11),Gis a constant of universal gravitation, i.e.G = 6.67408×10-11 N×m2/kg2；N x AndN z is a function of densityxCoordinates andzthe highest order of the coordinates is the highest order,iandjis in a density polynomialxAndzthe order of the corresponding order is set as,D i, j is a coefficient of a density polynomial;nandmis also thatxCoordinates andzorder of coordinates:nis between 0 andiin the first order in between, and in the second order,mis between 0 andjan order in between;C i n andC j m is the number of combinations; (x 0, z 0) Is the coordinate of the current observation point; function(s)F xz AndF zz is represented byi、j、n、mTwo piecewise functions are involved:
Figure 736444DEST_PATH_IMAGE019
（12）
Figure 166289DEST_PATH_IMAGE020
（13）
in the above formulae (12) to (13), (b)x k , z k ) Is the firstkCoordinates of a vertex, (x k + 1, z k + 1) Is the firstk The coordinates of the + 1 vertex points,pandqare respectively the firstkThe slope and intercept of the straight line on which the edges lie,p =(z k + 1 - z k )/(x k + 1 - x k ), q =(z k x k + 1 - z k + 1 x k )/(x k + 1 - x k )；Q = px 0 + q - z 0,；l 1andl 2also represents the order, but the value ranges of the two can be different due to different formulas,
Figure DEST_PATH_IMAGE035
、
Figure DEST_PATH_IMAGE036
、
Figure 407914DEST_PATH_IMAGE023
、
Figure DEST_PATH_IMAGE037
、
Figure 620721DEST_PATH_IMAGE025
and
Figure 264192DEST_PATH_IMAGE026
are all combination numbers;IandJrespectively represent two number sequences, wherein the number sequencesIThe calculation formula of (2) is as follows:
Figure DEST_PATH_IMAGE038
（14）
in the above formula (14)Q = px 0 + q - z 0, b = p(q - z 0) - x 0, c =1 + p 2，r k Andr k + 1respectively representkA vertex and the secondkThe distance between + 1 vertex and the current observation point; another array of numbersJIs a recursive sequence, if usedlIndicating the sequence number of the item in the sequence, the sequence is recurredJ l The general formula of (A) is as follows:
Figure 599358DEST_PATH_IMAGE028
（15）
in the above formula (15)Q = px 0 + q - z 0, a =(q - z 0)2 + x 0 2, b = p(q - z 0) - x 0, c =1 + p 2，r k Andr k + 1respectively representkA vertex and the secondkThe distance between + 1 vertex and the current observation point;I 0the formula (2) is shown as (14).
The first embodiment is as follows:
the invention will be explained and illustrated with reference to specific embodiments.
In order to further explain the implementation idea and implementation process of the method, the model shown in fig. 3 is used for explanation. The model simulates a valley of the southern state of california. Meanwhile, the conventional method of dividing the constant density cell is also used to illustrate the effectiveness and accuracy of the method.
S1: a rectangular plane coordinate system as shown in FIG. 3 is established in the plane of the section of the two-degree body, and the horizontal axisxRepresenting the horizontal distance, with the unit of km, and taking the right as positive; longitudinal axiszThe depth is expressed in km, and the orientation is positive;
s2: as shown in FIG. 3, the two-degree body is approximately hexagonal in cross-section, i.e.N k = 6, numbering the polygon vertices in the counter-clockwise direction, with vertex coordinates: (x 1, z 1) = (5, 0.1)，(x 2, z 2) = (8, 1.6)，(x 3, z 3) = (10, 1.8)，(x 4, z 4) = (12, 1.6)，(x 5, z 5) = (13.5, 0.6)，(x 6, z 6) = (15, 0.1)，(x 7, z 7) = (5, 0.1)；
S3: the density function of the second degree body isσ(x,z)=-602.4-18.22z+7.16z 2+44.55z 3+21.25z 4+0.027xz+0.0049xz 2-0.0091xz 3-0.06x 2+0.11x 2 z+0.097x 2 z 2+0.017x 3-0.0028x 3 zWherein the unit of density is kg/m3，xAndzthe units of (A) are km; input devicexHighest order N of coordinates x = 3，zHighest order of coordinatesN z = 4, coefficient of non-zero term beingD 0,0 = -602.4，D 0,1 = -18.22，D 0,2 = 7.16，D 0,3 = 44.55，D 0,4 = 21.25，D 1,1 = 0.027，D 1,2 = 0.0049，D 1,3 = -0.0091，D 2,0 = -0.06，D 2,1 = 0.11，D 2,2 = 0.097，D 3,0 = 0.017，D 3,1 = -0.0028；
S4: the measuring line is located on the earth's surfacez Over = 0, fromx = 0km tox = 20km one observation point is provided every 0.433km, for a total of 47 observation points, with observation point coordinates of (0,0), (0.433,0), (0.866,0), (1.299,0), (1.732,0), (2.165,0), (2.598,0), (3.031,0), (3.464,0), (3.897,0), (4.33,0), (4.763,0), (5.196,0), (5.629,0), (6.062,0), (6.495,0), (6.928,0), (7.361,0), (7.794,0), (8.227,0), (8.66,0), (9.093,0), (9.526,0), (9.959,0), (10.392,0), (10.825,0), (11.258,0), (11.691,0), (12.124,0), (12.557,0), (12.99,0), (13.423,0), (13.856,0), (14.289,0), (14.722, 15.155), (16.021, 15.588), (16.454,0), (16.887,0), (17.32,0), (17.753,0), (18.186,0), (18.619,0), (19.052,0), (19.485,0), (19.918, 0);
s5: the first observation point is designated as the current observation point, and the coordinates (0,0) thereof are designated as: (x 0,z 0)；
S6: for the first side of the polygon,k = 1，x 1 ≠ x 2to illustrate the edge andzslope of line in which axes are not parallelp= 0.5, interceptq= -2.4, discriminantQ = -2.4，QThe non-zero value indicates that the current observation point (0,0) is not collinear with the edge, so equations (10) - (13) are used for calculationu xz_1Andu zx_1calculated by equations (12) to (15)u zz_1Andu xx_1；
for the second side of the polygon,k = 2，x 2 ≠ x 3to illustrate the edge andzslope of line in which axes are not parallelp= 0.1, interceptq= 0.8, discriminantQ = 0.8，QThe non-zero value indicates that the current observation point (0,0) is not collinear with the edge, so equations (10) - (13) are used for calculationu xz_2Andu zx_2calculated by equations (12) to (15)u zz_2Andu xx_2；
for the third side of the polygon,k = 3，x 3 ≠ x 4to illustrate the edge andzslope of line in which axes are not parallelp= 0.1, interceptq= 2.8, discriminantQ = 2.8，QThe non-zero value indicates that the current observation point (0,0) is not collinear with the edge, so equations (10) - (13) are used for calculationu xz_3Andu zx_3calculated by equations (12) to (15)u zz_3Andu xx_3；
for the fourth side of the polygon,k = 4，x 4 ≠ x 5to illustrate the edge andzslope of line in which axes are not parallelp= 0.6667, interceptq= 9.6, discriminantQ = 9.6，QThe non-zero value indicates that the current observation point (0,0) is not collinear with the edge, so equations (10) - (13) are used for calculationu xz_4Andu zx_4calculated by equations (12) to (15)u zz_4Andu xx_4；
for the fifth side of the polygon,k = 5，x 5 ≠ x 6to illustrate the edge andzslope of line in which axes are not parallelp= -0.3333, interceptq= 5.1, discriminantQ = 5.1，QThe non-zero value indicates that the current observation point (0,0) is not collinear with the edge, so equations (10) - (13) are used for calculationu xz_5Andu zx_5calculated by equations (12) to (15)u zz_5Andu xx_5；
for the sixth side of the polygon,k = 6，x 6 ≠ x 7to illustrate the edge andzslope of line in which axes are not parallelp= 0, interceptq= 0.1, discriminantQ = 0.1，QThe non-zero value indicates that the current observation point (0,0) is not collinear with the edge, so equations (10) - (13) are used for calculationu xz_6Andu zx_6calculated by equations (12) to (15)u zz_6Andu xx_6；
s7: adding the contribution values of the 6 edges to obtain the value of the gravity gradient tensor at the current observation point (0,0),
U xz (0,0)= u xz_1+u xz_2+u xz_3+u xz_4+u xz_5+u xz_6，
U zx (0,0)= u zx_1+u zx_2+u zx_3+u zx_4+u zx_5+u zx_6，
U zz (0,0)= u zz_1+u zz_2+u zz_3+u zz_4+u zz_5+u zz_6，
U xx (0,0)= u xx_1+u xx_2+u xx_3+u xx_4+u xx_5+u xx_6；
s8: judging whether there is any observation point which is not calculated, if so, designating the observation point as the current observation point, wherein the coordinate of the observation point is (x 0,z 0) Returning to step S6, a new calculation is performed; if not, executing the next step;
s9: outputting the gravity gradient tensor elements on all observation pointsU xz 、U zx 、U xx 、U zz The calculated values of (A) are plotted as shown in FIG. 4.
S10: in order to verify the validity of the calculation result of the method, the gravity gradient tensor generated by the model in the figure 3 is calculated on the same measuring point by utilizing the method of dividing the constant density unit. First, utilizexThe length in the direction is 0.1km,zThe mesh with the length of 0.01km in the direction divides the polygonal model in the figure 3 into 10747 constant-density small cells, wherein the small cells in the middle part are all divided into rectangles, and the boundary part is divided into right-angled triangles or right-angled trapezoids in order to be matched with the shape of the model as much as possible; then, the densities at the apex and the center point of each small cell are calculated using the density function shown in formula (16), and the average of them is taken as the density value of the whole cell; finally, calculating the value of the gradient tensor generated by each small polygon unit on the measuring point by using the existing analytic formula about the constant density polygon gravity gradient tensor, and accumulating the values to obtain the value of the gravity gradient tensor generated by the whole model on the measuring point; the variation curve calculated by this constant density method is shown in fig. 4.
And S11, calculating the root mean square error between the result of the method and the result of the constant density method. Then, the number of the divided units in the constant density method is sequentially adjusted, the root mean square error between the results of the constant density method and the text method under different unit numbers is recorded, and the CPU time consumed by the two methods is recorded, and the result is shown in fig. 5.
FIG. 2 is a geometric relationship diagram of observation points and a dyadic object according to the present invention. Fig. 3 is a variable density hexagonal model simulating a valley. FIG. 4 is a comparison of the results of the gravity gradient tensor produced by the model of FIG. 3 calculated using the present method and a conventional constant density method, respectively. As can be seen from fig. 4, the calculation results of the two methods are consistent, which illustrates the correctness of the method. Fig. 5 (a) is a graph showing the relationship between the root mean square error of the constant density method and the number of divided units, and (b) is a graph showing the CPU time consumed by each of the two methods for different numbers of discrete units. As can be seen from (a) in fig. 5, as the number of the divided units gradually increases, the root mean square error of the result of the constant density method and the method herein gradually decreases. It can be found from the graph (b) that the accuracy of the constant density method is improved and the efficiency is gradually reduced: if the root mean square error of the calculation result of the constant density method is less than 0.001E, the model in FIG. 3 needs to be divided into 7000 small units, and the calculation result of the constant density method is used at this timeU xz AndU zx requires 20.4036s and calculationU xx AndU zz 19.3334s is needed, however, only 7.5153s and 12.4313s are needed for calculation by the method, and the calculation time of the constant density method is 2.71 times and 1.56 times of that of the method respectively. This demonstrates that the method is superior to conventional methods in terms of both accuracy and efficiency.
The above description is only exemplary of the present invention and should not be taken as limiting, and all equivalent modifications, equivalents, and improvements made within the spirit and principle of the present invention should be included in the scope of the present invention.

Claims (2)

1. A two-degree variable density body gravity gradient tensor calculation method is characterized by comprising the following steps:
the method comprises the following steps: a plane rectangular coordinate system is established in the plane of the section of the two-degree body, and a transverse axisxIndicating a horizontal distance, taking the right as positive; longitudinal axiszIndicating depth, the orientation is positive;
step two: the section of the two-degree body is approximated to be one withN k The polygon of the side is started from any vertex, and the coordinates of each vertex are determined to be (in turn) along the counterclockwise directionx k , z k ) Whereink = 1, 2, …, N k , N k + 1；
Step three: using a function of binary polynomials
Figure 156356DEST_PATH_IMAGE001
Representing the change in density of the second degree volume, determiningxHighest order of coordinatesN x Andzhighest order of coordinatesN z And coefficients of the polynomialD i, j ；
Step four: arranging a plurality of observation points along the observation measuring line and determining coordinates of the observation points, wherein the abscissa of each observation point cannot be equal to the abscissa of the vertex of the polygon;
step five: specifying a current observation point having coordinates ofx 0, z 0)，x 0 ≠ x k ；
Step six: respectively calculating the number of polygons on the current observation pointk(k = 1, 2, …, N k ) Contribution U of gravity gradient tensor U corresponding to edges k Wherein
Figure 706024DEST_PATH_IMAGE002
，
Figure 631254DEST_PATH_IMAGE003
；
Step seven: will be provided withN k Adding the contribution values corresponding to the edges to obtain the value of the gravity gradient tensor at the current observation point, wherein,
U xz (x 0,z 0)=u xz_1+u xz_2+…+
Figure 855562DEST_PATH_IMAGE004
，
U zx (x 0,z 0)=u zx_1+u zx_2+…+
Figure 917059DEST_PATH_IMAGE005
，
U xx (x 0,z 0)=u xx_1+u xx_2+…+
Figure 201410DEST_PATH_IMAGE006
，
U zz (x 0,z 0)=u zz_1+u zz_2+…+
Figure 613937DEST_PATH_IMAGE007
；
step eight: judging whether observation points which are not calculated exist, if so, returning to the step five, designating the observation points as current observation points, and performing new calculation; if not, executing the next step;
step nine: outputting the gravity gradient tensor elements on all observation pointsU xz 、U zx 、U xx 、U zz The calculated value of (a).
2. The method of claim 1, wherein the sixth step is performed to calculate the sixth stepkGravity gradient tension with corresponding edgesThe formula for the quantity contribution value is:
is as followskEdge andzwhen the axes are parallel to each other, the axes are parallel,
Figure 579619DEST_PATH_IMAGE008
（1）
Figure 761201DEST_PATH_IMAGE009
（2）
in the above-mentioned formulae (1) to (2),Gis a constant of universal gravitation, i.e.G = 6.67408×10-11N×m2/kg2；N x AndN z is a function of densityxCoordinates andzthe highest order of the coordinates is the highest order,iandjis a non-zero termD i, j x i z j InxCoordinates andzthe order to which the coordinates correspond to,D i, j is a coefficient of a density polynomial;nandmis also thatxCoordinates andzorder of coordinates:nis between 0 andiin the first order in between, and in the second order,mis between 0 andjan order in between;C i n andC j m is the number of combinations; (x 0, z 0) Is the coordinate of the current observation point; function(s)E xz AndE zz is represented byi、j、n、mTwo piecewise functions are involved:
Figure 950874DEST_PATH_IMAGE010
（3）
Figure 850697DEST_PATH_IMAGE011
（4）
in the above formulae (3) to (4), (b)x k , z k ) Is the firstkCoordinates of a vertex, (x k + 1, z k + 1) Is the firstk + 1 vertex coordinates, andr k andr k + 1respectively representkA vertex and the secondk+ 1 distance between the vertex and the current observation point, anr k + 1=[(x k + 1 - x 0)2 + (z k + 1 - z 0)2]1/2，r k =[(x k - x 0)2 + (z k - z 0)2]1/2；HRepresents a recursive sequence, if usedlIndicating the sequence number of the item in the sequence, the sequence is recurredH l The general formula of (A) is as follows:
Figure 416808DEST_PATH_IMAGE012
（5）
in the above formula (5)KRepresents another recursion sequence, if usedsIndicating the sequence number of the item in the sequence, the sequence is recurredK s The general formula of (A) is as follows:
Figure 718476DEST_PATH_IMAGE013
（6）
when it is at firstkEdge andzthe axes are not parallel, and the current point of observation is the secondkWhen the edges are on the same straight line,
Figure 16733DEST_PATH_IMAGE014
（7）
Figure 403852DEST_PATH_IMAGE015
（8）
in the above-mentioned formulae (7) to (8),Gis always ledForce constant, i.e.G = 6.67408×10-11 N×m2/kg2；N x AndN z is a function of densityxCoordinates andzthe highest order of the coordinates is the highest order,iandjis a non-zero termD i, j x i z j InxCoordinates andzthe order to which the coordinates correspond to,D i, j is a coefficient of a density polynomial; (x 0, z 0) Is the coordinate of the current observation point; and (a)x k , z k ) Is the firstkCoordinates of a vertex, (x k + 1, z k + 1) Is the firstk The coordinates of the + 1 vertex points,pandqare respectively the firstkThe slope and intercept of the straight line on which the edges lie,p =(z k + 1 - z k )/(x k + 1 - x k ),q =(z k x k + 1 - z k + 1 x k )/(x k + 1 - x k ) (ii) a And arrays ofMThe calculation formula of (2) is as follows:
Figure 773654DEST_PATH_IMAGE016
（9）
in the above formula (9)c =1 + p 2, d = [p(q - z 0) - x 0]/c；
Third whenkEdge andzthe axes are not parallel, and the current point of observation is the secondkWhen the edges are not in the same straight line,
Figure 664249DEST_PATH_IMAGE017
（10）
Figure 461304DEST_PATH_IMAGE018
（11）
in the above-mentioned formulae (10) to (11),Gis a constant of universal gravitation, i.e.G = 6.67408×10-11 N×m2/kg2；N x AndN z is a function of densityxCoordinates andzthe highest order of the coordinates is the highest order,iandjis in a density polynomialxAndzthe order of the corresponding order is set as,D i, j is a coefficient of a density polynomial;nandmis also thatxCoordinates andzorder of coordinates:nis between 0 andiin the first order in between, and in the second order,mis between 0 andjan order in between;C i n andC j m is the number of combinations; (x 0, z 0) Is the coordinate of the current observation point; function(s)F xz AndF zz is represented byi、j、n、mTwo piecewise functions are involved:
Figure 7823DEST_PATH_IMAGE019
（12）
Figure 181315DEST_PATH_IMAGE020
（13）
in the above formulae (12) to (13), (b)x k , z k ) Is the firstkCoordinates of a vertex, (x k + 1, z k + 1) Is the firstk The coordinates of the + 1 vertex points,pandqare respectively the firstkThe slope and intercept of the straight line on which the edges lie,p =(z k + 1 - z k )/(x k + 1 - x k ), q =(z k x k + 1 - z k + 1 x k )/(x k + 1 - x k )；Q = px 0 + q - z 0,；l 1andl 2also represents the order, but the value ranges of the two can be different due to different formulas,
Figure 191997DEST_PATH_IMAGE021
、
Figure 159953DEST_PATH_IMAGE022
、
Figure 256085DEST_PATH_IMAGE023
、
Figure 967689DEST_PATH_IMAGE024
、
Figure 260306DEST_PATH_IMAGE025
and
Figure 399163DEST_PATH_IMAGE026
are all combination numbers;IandJrespectively represent two number sequences, wherein the number sequencesIThe calculation formula of (2) is as follows:
Figure 982591DEST_PATH_IMAGE027
（14）
in the above formula (14)Q = px 0 + q - z 0, b = p(q - z 0) - x 0, c =1 + p 2，r k Andr k + 1respectively representkA vertex and the secondkThe distance between + 1 vertex and the current observation point; another array of numbersJIs a recursive sequence, if usedlIndicating the sequence number of the item in the sequence, the sequence is recurredJ l The general formula of (A) is as follows:
Figure 232307DEST_PATH_IMAGE028
（15）
in the above formula (15)Q = px 0 + q - z 0, a =(q - z 0)2 + x 0 2, b = p(q - z 0) - x 0, c =1 + p 2，r k Andr k + 1respectively representkA vertex and the secondkThe distance between + 1 vertex and the current observation point;I 0the formula (2) is shown as (14).
CN201910398005.7A 2019-05-14 2019-05-14 Method for calculating gravity gradient tensor of second-degree variable density body Active CN111950108B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910398005.7A CN111950108B (en) 2019-05-14 2019-05-14 Method for calculating gravity gradient tensor of second-degree variable density body

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910398005.7A CN111950108B (en) 2019-05-14 2019-05-14 Method for calculating gravity gradient tensor of second-degree variable density body

Publications (2)

Publication Number Publication Date
CN111950108A true CN111950108A (en) 2020-11-17
CN111950108B CN111950108B (en) 2024-04-16

Family

ID=73336016

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910398005.7A Active CN111950108B (en) 2019-05-14 2019-05-14 Method for calculating gravity gradient tensor of second-degree variable density body

Country Status (1)

Country Link
CN (1) CN111950108B (en)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100250185A1 (en) * 2009-03-27 2010-09-30 Qinetiq Limited Method for detection of gravitational anomalies
CN106855904A (en) * 2017-01-10 2017-06-16 桂林理工大学 A kind of Two bodies gravity anomaly computational methods
CN107024723A (en) * 2017-06-16 2017-08-08 桂林理工大学 A kind of Two bodies the Magnetic Field Numerical Calculation method
CN107402409A (en) * 2017-09-26 2017-11-28 西南石油大学 A kind of three-dimensional irregular stratum fluctuating interface gravity forward modeling method
CN107491411A (en) * 2017-06-23 2017-12-19 中国海洋大学 Gravity anomaly inversion method based on N rank multinomial density functions
CN107561592A (en) * 2017-09-11 2018-01-09 中南大学 A kind of density is the gravitational field forward modeling method of polynomial arbitrary polyhedron

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100250185A1 (en) * 2009-03-27 2010-09-30 Qinetiq Limited Method for detection of gravitational anomalies
CN106855904A (en) * 2017-01-10 2017-06-16 桂林理工大学 A kind of Two bodies gravity anomaly computational methods
CN107024723A (en) * 2017-06-16 2017-08-08 桂林理工大学 A kind of Two bodies the Magnetic Field Numerical Calculation method
CN107491411A (en) * 2017-06-23 2017-12-19 中国海洋大学 Gravity anomaly inversion method based on N rank multinomial density functions
CN107561592A (en) * 2017-09-11 2018-01-09 中南大学 A kind of density is the gravitational field forward modeling method of polynomial arbitrary polyhedron
CN107402409A (en) * 2017-09-26 2017-11-28 西南石油大学 A kind of three-dimensional irregular stratum fluctuating interface gravity forward modeling method

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
JIANZHONG ZHANG ET AL: "Gravity anomalies of 2-D bodies with variable density contrast", GEOPHYSICS, vol. 66, no. 03, pages 809 - 813 *
LI JIANG ET AL: "Analytic Expressions for the Gravity Gradient Tensor of 3D Prisms with Depth-Dependent Density", SURV GEOPHYS, pages 1 - 27 *
ZHEN JIA ET AL: "Potential fields and their partial derivatives produced by a 2D homogeneous polygonal source: A summary with some revisions", GEOPHYSICS, vol. 76, no. 04, pages 29 *
张建中 等: "密度随深度变化的二度体重力正演公式", 石油地球物理勘探, vol. 35, no. 02, pages 202 - 207 *
郭武林: "岩石密度空间变化时重力效应的计算方法", 物探与化探, no. 04, pages 226 - 236 *

Also Published As

Publication number Publication date
CN111950108B (en) 2024-04-16

Similar Documents

Publication Publication Date Title
US7480205B2 (en) 3D fast fault restoration
Li 3-D inversion of gravity gradiometer data
CN109725345B (en) First-arrival forward modeling method and device
CN113553748B (en) Three-dimensional magnetotelluric forward modeling numerical simulation method
Reimond et al. Spheroidal and ellipsoidal harmonic expansions of the gravitational potential of small Solar System bodies. Case study: Comet 67P/Churyumov‐Gerasimenko
CN103149585A (en) Elastic migration seismic wave field construction method and elastic migration seismic wave field construction device
CN108984939A (en) Three-dimensional Gravity field of force forward modeling method based on 3D Gauss-FFT
CN106483559A (en) A kind of construction method of subsurface velocity model
CN111158059B (en) Gravity inversion method based on cubic B spline function
WO2012001388A2 (en) Gravity survey data processing
Golubev et al. Raising convergence order of grid-characteristic schemes for 2D linear elasticity problems using operator splitting
Long et al. Three-dimensional forward modelling of gravity data using mesh-free methods with radial basis functions and unstructured nodes
CN102360084B (en) Universal gravitation interference calculation method based on tetrahedral mass element division for pure gravity rail
CN113109883B (en) Method for forward modeling satellite gravity field based on isoparametric transformation global discrete grid spherical coordinates
CN106662665B (en) The interpolation and convolution of rearrangement for the processing of faster staggered-mesh
CN103675905A (en) Optimal coefficient acquisition method and device, and related wave field simulating method and device
Zhou et al. An iterative factored topography-dependent eikonal solver for anisotropic media
Li et al. Fast 3D forward modeling of the magnetic field and gradient tensor on an undulated surface
WO2018067014A1 (en) Seismic modeling
CN111950108A (en) Two-degree variable density body gravity gradient tensor calculation method
CN112748471B (en) Gravity-magnetic data continuation and conversion method of unstructured equivalent source
CN115270579A (en) Second-order acoustic wave equation finite difference numerical simulation parameter selection method
Wéber Optimizing model parameterization in 2D linearized seismic traveltime tomography
Sava et al. Huygens wavefront tracing: A robust alternative to ray tracing
CN108072899B (en) Self-adaptive implementation method of discontinuous Galerkin finite element seismic numerical simulation algorithm

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