CN111950108B - Method for calculating gravity gradient tensor of second-degree variable density body - Google Patents
Method for calculating gravity gradient tensor of second-degree variable density body Download PDFInfo
- Publication number
- CN111950108B CN111950108B CN201910398005.7A CN201910398005A CN111950108B CN 111950108 B CN111950108 B CN 111950108B CN 201910398005 A CN201910398005 A CN 201910398005A CN 111950108 B CN111950108 B CN 111950108B
- Authority
- CN
- China
- Prior art keywords
- vertex
- coordinate
- coordinates
- kth
- observation point
- 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.)
- Active
Links
- 230000005484 gravity Effects 0.000 title claims abstract description 50
- 238000000034 method Methods 0.000 title claims abstract description 39
- 238000004364 calculation method Methods 0.000 claims abstract description 40
- 238000003491 array Methods 0.000 claims description 3
- 238000009825 accumulation Methods 0.000 abstract description 2
- 230000008859 change Effects 0.000 abstract description 2
- 230000004044 response Effects 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 2
- 238000007796 conventional method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
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 method for calculating gravity gradient tensor of a two-degree variable density body, and belongs to the technical field of geophysical exploration. The invention comprises the following steps: establishing a plane rectangular coordinate system according to the research area; approximating the cross section of the binary body to be a polygon, and determining the coordinates of the vertex; simulating density change of the binary body by using a binary polynomial function; setting an observation point along an observation line and determining coordinates thereof; designating a current observation point, calculating a gravity gradient tensor contribution value corresponding to each side of the polygon on the current observation point, and adding the gravity gradient tensor contribution values to obtain a gravity gradient tensor value on the current observation point; repeating the above calculation steps until the values of the gravity gradient tensors at all observation points are output. The method can calculate and obtain the accurate analytic solution of the gravity gradient tensor by using an analytic formula deduced in the space domain. Compared with a constant density unit accumulation calculation method, the method has higher efficiency under the same precision requirement.
Description
Technical Field
The invention relates to a calculation method of a weight gradient tensor of a two-degree variable density body, and belongs to the technical field of geophysical exploration.
Background
If a geologic volume has a dimension in the strike direction that is much greater than the dimension perpendicular to its strike direction, then it can be considered a binary volume. When the density distribution of the geologic body is complex, the geologic body is generally divided into a plurality of simple regular-shaped units with uniform density, the gravity response generated by each regular-shaped unit is calculated by using the existing formula, and then the gravity response 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 units. With the increase of the number of the dividing units, the calculation accuracy is higher and higher, but the calculation efficiency is also continuously reduced.
In addition to the above-described method of using the gravity response calculation formula for constant density bodies, gravity response calculation formulas for variable density bodies have also been developed. The gravity anomaly calculation formula which uses a polynomial function of space coordinates to express the change of density and uses a polygon to approximate the shape of a binary body has been deduced, so as to form a corresponding gravity anomaly analysis calculation method.
Gravity gradient data has better resolution than gravity anomaly data. With the advent of gravity gradiometers, one is increasingly studying subsurface geologic structures and finding mineral resources by observing gravity gradient data. The method of 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 polygonal body, but the calculation accuracy is not high due to truncation errors. Therefore, the analysis and calculation method of the gravity gradient tensor of the second-degree variable density polygonal body is developed, and the calculation accuracy of the gravity gradient tensor of the second-degree variable density polygonal body can be obviously improved.
Disclosure of Invention
Aiming at the defects in the prior art, the invention provides a calculation method of a gravity gradient tensor of a two-degree variable density body, 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 space domain. Wherein the density of the polygon model is represented by a binary polynomial function of any order of horizontal and vertical coordinates.
The invention is realized by adopting the following technical scheme:
Step one: establishing a plane rectangular coordinate system in a plane where the section of the binary body is positioned, wherein the horizontal axis x represents the horizontal distance and is taken to be positive to the right; the vertical axis z represents depth, oriented positive;
Step two: approximating the two-degree body section as a polygon with N k sides, and sequentially determining coordinates of each vertex from any vertex in a counterclockwise direction as (x k, zk), wherein k=1, 2, … and N k, Nk +1;
Step three: representing density variation of the dyadic with a binary polynomial function , determining a highest order N x of the x-coordinate and a highest order N z of the z-coordinate, and coefficients D i, j of the polynomial;
step four: setting a plurality of observation points along an observation line and determining coordinates of the observation points, wherein the abscissa of the observation points cannot be equal to the abscissa of the polygon vertexes;
step five: the current observation point is specified, and the coordinates thereof are (x 0, z0),x0 ≠ xk;
step six: calculating a contribution value U k of a gravity gradient tensor U corresponding to the kth (k=1, 2, …, N k) side of the polygon on the current observation point respectively, wherein ,/>;
Step seven: adding the contribution values corresponding to the N k sides to obtain the value of the gravity gradient tensor on the current observation point, wherein,
Uxz(x0,z0)=uxz_1+uxz_2+…+，
Uzx(x0,z0)=uzx_1+uzx_2+…+，
Uxx(x0,z0)=uxx_1+uxx_2+…+，
Uzz(x0,z0)=uzz_1+uzz_2+…+；
Step eight: judging whether observation points which are not calculated exist or not, 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: and outputting the calculated values of the gravity gradient tensor elements U xz、Uzx、Uxx、Uzz on all observation points.
Further, the formula for calculating the gravity gradient tensor contribution value corresponding to the kth edge in the sixth step is as follows:
① When the kth edge is parallel to the z-axis,
（1）
（2）
In the formulas (1) - (2), G is a gravitational constant, that is, g= 6.67408 ×10 -11N×m2/kg2;Nx and N z are the highest orders of the x-coordinate and the z-coordinate in the density function, i and j are orders corresponding to the x-coordinate and the z-coordinate in the non-zero term D i, jxizj, and D i, j is a coefficient of the density polynomial; n and m are also orders of x and z coordinates: n is an order between 0 and i, m is an order between 0 and j; c i n and C j m are combined numbers; (x 0, z0) is the coordinates of the current observation point; functions E xz and E zz represent two piecewise functions related to i, j, n, m:
（3）
（4）
In the above equations (3) - (4), (x k, zk) is the coordinates of the kth vertex, (x k + 1, zk + 1) is the coordinates of the (k+1) th vertex, r k and r k + 1 respectively represent the distances between the kth vertex and the (k+1) th vertex to the current observation point, and rk + 1=[(xk + 1 - x0)2 + (zk + 1 - z0)2]1/2,rk=[(xk - x0)2 + (zk - z0)2]1/2;H represents a recursive sequence, and if the sequence number of the entries in this sequence is denoted by l, the general formula of the recursive sequence H l is:
（5）
K in the above formula (5) represents another recursive sequence, and if s represents the sequence number of the term in this sequence, the general formula of the recursive sequence K s is:
（6）
② When the kth edge is not parallel to the z-axis and the current observation point is on the same line as the kth edge,
（7）
（8）
In the formulas (7) - (8), G is a gravitational constant, that is, g= 6.67408 ×10 -11 N×m2/kg2;Nx and N z are the highest orders of the x-coordinate and the z-coordinate in the density function, i and j are orders corresponding to the x-coordinate and the z-coordinate in the non-zero term D i, jxizj, and D i, j is a coefficient of the density polynomial; (x 0, z0) is the coordinates of the current observation point; and (x k, zk) is the coordinates of the kth vertex, (x k + 1, zk + 1) is the coordinates of the (k+1) th vertex, p and q are the slope and intercept ,p =(zk + 1 - zk)/(xk + 1 - xk),q =(zkxk + 1 - zk + 1xk)/(xk + 1 - xk); of the straight line where the kth edge is located, respectively, and the calculation formula of the number sequence M is:
（9）
In the formula (9), c=1+p 2, d = [p(q - z0) - x0 ]/c;
③ When the kth edge is not parallel to the z-axis and the current observation point is not on the same straight line as the kth edge,
（10）
（11）
In the formulas (10) - (11), G is a gravitational constant, that is, g= 6.67408 ×10 -11 N×m2/kg2;Nx and N z are the highest orders of the x-coordinate and z-coordinate in the density function, i and j are orders corresponding to x and z in the density polynomial, and D i, j is a coefficient of the density polynomial; n and m are also orders of x and z coordinates: n is an order between 0 and i, m is an order between 0 and j; c i n and C j m are combined numbers; (x 0, z0) is the coordinates of the current observation point; the functions F xz and F zz represent two piecewise functions related to i, j, n, m:
（12）
（13）
In the formulas (12) - (13), the (x k, zk) is the coordinate of the kth vertex, (x k + 1, zk + 1) is the coordinate of the kth+1th vertex, p and q are the slope and the intercept ,p =(zk + 1 - zk)/(xk + 1 - xk), q =(zkxk + 1 - zk + 1xk)/(xk + 1 - xk);Q = px0 + q - z0,;l1 and l 2 of the straight line where the kth edge is located respectively, and the order is represented, but the value ranges of the two are different according to the different formulas, and 、/>、/>、/>、/> and/> are both combination numbers; i and J represent two arrays, respectively, wherein the calculation formula of the array I is as follows:
（14）
In the formula (14), q=px 0 + q - z0, b = p(q - z0) - x0, c =1 + p2,rk and r k + 1 represent the distances between the kth vertex and the (k+1) th vertex and the current observation point, respectively; the other column J is a recursive column, and if the sequence number of the term in this column is denoted by l, the general formula of the recursive column J l is:
（15）
Q = px0 + q - z0, a =(q - z0)2 + x0 2, b = p(q - z0) - x0, c =1 + p2,rk and r k + 1 in the above formula (15) represent the distances between the kth vertex and the (k+1) th vertex to the current observation point, respectively; the calculation formula of I 0 is shown as formula (14).
The beneficial effects of the invention are as follows: the gravity gradient tensor calculation method of the two-degree variable density body uses a polygon with a binary polynomial density function, and can flexibly approximate the two-degree geologic body with complex underground shape and variable density in actual exploration. This is because any two-dimensional density distribution can be well approximated by a binary polynomial of sufficiently high order, and any complex boundary can be well approximated by any polygon of sufficiently many sides. In addition, the calculation method can obtain the analytic solution of the gravity gradient tensor. The calculation method is simple, accurate, efficient and high in adaptability, and the inversion algorithm of the gravity gradient data is easy to establish based on the gravity gradient formula.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a diagram showing the geometric relationship between the observation point and the binary according to the present invention;
FIG. 3 is a variable density hexagonal model simulating a lobe;
FIG. 4 is a comparison of the results of the model-generated gravity gradient tensor 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 of the normal density volume superposition method and the CPU time consumed by each of the two methods.
Detailed Description
The foregoing and other objects, features and advantages of the invention will be apparent from the following more particular description of the invention, as illustrated in the accompanying drawings. The flow chart of the invention, as shown in fig. 1, comprises the following steps:
Step one: establishing a plane rectangular coordinate system in a plane where the section of the binary body is positioned, wherein the horizontal axis x represents the horizontal distance and is taken to be positive to the right; the vertical axis z represents depth, oriented positive;
Step two: approximating the two-degree body section as a polygon with N k sides, and sequentially determining coordinates of each vertex from any vertex in a counterclockwise direction as (x k, zk), wherein k=1, 2, … and N k, Nk +1;
Step three: representing density variation of the dyadic with a binary polynomial function , determining a highest order N x of the x-coordinate and a highest order N z of the z-coordinate, and coefficients D i, j of the polynomial;
step four: setting a plurality of observation points along an observation line and determining coordinates of the observation points, wherein the abscissa of the observation points cannot be equal to the abscissa of the polygon vertexes;
step five: the current observation point is specified, and the coordinates thereof are (x 0, z0),x0 ≠ xk;
Step six: calculating a contribution value U k of a gravity gradient tensor U corresponding to the kth (k=1, 2, …, N k) side of the polygon on the current observation point respectively, wherein ,/>;
Step seven: adding the contribution values corresponding to the N k sides to obtain the value of the gravity gradient tensor on the current observation point, wherein,
Uxz(x0,z0)=uxz_1+uxz_2+…+，
Uzx(x0,z0)=uzx_1+uzx_2+…+，
Uxx(x0,z0)=uxx_1+uxx_2+…+，
Uzz(x0,z0)=uzz_1+uzz_2+…+；
Step eight: judging whether observation points which are not calculated exist or not, 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: and outputting the calculated values of the gravity gradient tensor elements U xz、Uzx、Uxx、Uzz on all observation points.
The specific way of calculating the gravity gradient tensor contribution value corresponding to the kth edge in the step six is as follows:
Firstly, judging whether the kth edge is parallel to the z axis, and if x k= xk + 1 is x, describing that the kth edge is parallel; if x k≠ xk + 1, then the description is not parallel;
Then, judging whether the current observation point (x 0, z0) and the kth side are on the same straight line or not, firstly calculating the slope p= (z k + 1 - zk)/(xk + 1 - xk) and the intercept q= (z kxk + 1 - zk + 1xk)/(xk + 1 - xk) of the kth side, and then calculating the value of the discriminant Q=px 0 + q - z0, and if Q=0, indicating that the two sides are on the same straight line; if Q is not equal to 0, the two are not on the same straight line;
Finally, selecting different calculation formulas to calculate the value of u xz_k、uzx_k、uxx_k、uzz_k according to the judgment result, wherein the specific calculation formulas comprise:
① When the kth edge is parallel to the z-axis,
（1）
（2）
In the formulas (1) - (2), G is a gravitational constant, that is, g= 6.67408 ×10 -11N×m2/kg2;Nx and N z are the highest orders of the x-coordinate and the z-coordinate in the density function, i and j are orders corresponding to the x-coordinate and the z-coordinate in the non-zero term D i, jxizj, and D i, j is a coefficient of the density polynomial; n and m are also orders of x and z coordinates: n is an order between 0 and i, m is an order between 0 and j; c i n and C j m are combined numbers; (x 0, z0) is the coordinates of the current observation point; functions E xz and E zz represent two piecewise functions related to i, j, n, m:
（3）
（4）
In the above equations (3) - (4), (x k, zk) is the coordinates of the kth vertex, (x k + 1, zk + 1) is the coordinates of the (k+1) th vertex, r k and r k + 1 respectively represent the distances between the kth vertex and the (k+1) th vertex to the current observation point, and rk + 1=[(xk + 1 - x0)2 + (zk + 1 - z0)2]1/2,rk=[(xk - x0)2 + (zk - z0)2]1/2;H represents a recursive sequence, and if the sequence number of the entries in this sequence is denoted by l, the general formula of the recursive sequence H l is:
（5）
K in the above formula (5) represents another recursive sequence, and if s represents the sequence number of the term in this sequence, the general formula of the recursive sequence K s is:
（6）
② When the kth edge is not parallel to the z-axis and the current observation point is on the same line as the kth edge,
（7）
（8）
In the formulas (7) - (8), G is a gravitational constant, that is, g= 6.67408 ×10 -11 N×m2/kg2;Nx and N z are the highest orders of the x-coordinate and the z-coordinate in the density function, i and j are orders corresponding to the x-coordinate and the z-coordinate in the non-zero term D i, jxizj, and D i, j is a coefficient of the density polynomial; (x 0, z0) is the coordinates of the current observation point; and (x k, zk) is the coordinates of the kth vertex, (x k + 1, zk + 1) is the coordinates of the (k+1) th vertex, p and q are the slope and intercept ,p =(zk + 1 - zk)/(xk + 1 - xk),q =(zkxk + 1 - zk + 1xk)/(xk + 1 - xk); of the straight line where the kth edge is located, respectively, and the calculation formula of the number sequence M is:
（9）
In the formula (9), c=1+p 2, d = [p(q - z0) - x0 ]/c;
③ When the kth edge is not parallel to the z-axis and the current observation point is not on the same straight line as the kth edge,
（10）
（11）
In the formulas (10) - (11), G is a gravitational constant, that is, g= 6.67408 ×10 -11 N×m2/kg2;Nx and N z are the highest orders of the x-coordinate and z-coordinate in the density function, i and j are orders corresponding to x and z in the density polynomial, and D i, j is a coefficient of the density polynomial; n and m are also orders of x and z coordinates: n is an order between 0 and i, m is an order between 0 and j; c i n and C j m are combined numbers; (x 0, z0) is the coordinates of the current observation point; the functions F xz and F zz represent two piecewise functions related to i, j, n, m:
（12）
（13）
In the formulas (12) - (13), the (x k, zk) is the coordinate of the kth vertex, (x k + 1, zk + 1) is the coordinate of the kth+1th vertex, p and q are the slope and the intercept ,p =(zk + 1 - zk)/(xk + 1 - xk), q =(zkxk + 1 - zk + 1xk)/(xk + 1 - xk);Q = px0 + q - z0,;l1 and l 2 of the straight line where the kth edge is located respectively, and the order is represented, but the value ranges of the two are different according to the different formulas, and 、/>、/>、/>、/> and/> are both combination numbers; i and J represent two arrays, respectively, wherein the calculation formula of the array I is as follows:
（14）
In the formula (14), q=px 0 + q - z0, b = p(q - z0) - x0, c =1 + p2,rk and r k + 1 represent the distances between the kth vertex and the (k+1) th vertex and the current observation point, respectively; the other column J is a recursive column, and if the sequence number of the term in this column is denoted by l, the general formula of the recursive column J l is:
（15）
Q = px0 + q - z0, a =(q - z0)2 + x0 2, b = p(q - z0) - x0, c =1 + p2,rk and r k + 1 in the above formula (15) represent the distances between the kth vertex and the (k+1) th vertex to the current observation point, respectively; the calculation formula of I 0 is shown as formula (14).
Embodiment one:
The invention is explained and illustrated below in connection with specific embodiments.
To further illustrate the implementation of the method, the model shown in fig. 3 is used for illustration. The model simulates a valley in the south of california. Meanwhile, a conventional method of dividing a constant density unit is also used to illustrate the effectiveness and accuracy of the present method.
S1: establishing a plane rectangular coordinate system shown in figure 3 in a plane where the section of the binary body is located, wherein the horizontal axis x represents the horizontal distance, the unit is km, and the right direction is positive; the vertical axis z represents depth in km and orientation positive;
S2: as shown in fig. 3, the bi-level cross section is approximately hexagonal, i.e., N k =6, and the polygon vertices are numbered in the counterclockwise direction with vertex coordinates of each of the vertices ：(x1, z1) = (5, 0.1),(x2, z2) = (8, 1.6),(x3, z3) = (10, 1.8),(x4, z4) = (12, 1.6),(x5, z5) = (13.5, 0.6),(x6, z6) = (15, 0.1),(x7, z7) = (5, 0.1);
S3: the density function of the binary is σ(x,z)=-602.4-18.22z+7.16z2+44.55z3+21.25z4+0.027xz+0.0049xz2-0.0091xz3-0.06x2+0.11x2z+0.097x2z2+0.017x3-0.0028x3z,, wherein the unit of density is kg/m 3, and the units of x and z are km; the highest order of the input x-coordinate, N x = 3, the highest order of the z-coordinate, N z = 4, the coefficients of the non-zero terms are D0,0 = -602.4,D0,1 = -18.22,D0,2 = 7.16,D0,3 = 44.55,D0,4 = 21.25,D1,1 = 0.027,D1,2 = 0.0049,D1,3 = -0.0091,D2,0 = -0.06,D2,1 = 0.11,D2,2 = 0.097,D3,0 = 0.017,D3,1 = -0.0028;
S4: the measuring line is positioned on the surface z=0, one observation point is arranged every 0.433km from x=0 km to x=20 km, and the total of 47 observation points are 47 observation point coordinates (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,0),(15.155,0),(15.588,0),(16.021,0),(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: designating the first observation point as the current observation point, and designating the coordinates (0, 0) as (x 0,z0);
S6: for the first side of the polygon, k=1, x 1 ≠ x2, which indicates that the side is not parallel to the z-axis, the slope p=0.5 of the straight line where the side is located, the intercept q= -2.4, the discriminant q= -2.4, and the Q value is not zero, which indicates that the current observation point (0, 0) is not on the same straight line as the side, so that formulas (10) - (13) are used to calculate u xz_1 and u zx_1, and formulas (12) - (15) are used to calculate u zz_1 and u xx_1;
For the second side of the polygon, k=2, x 2 ≠ x3, which indicates that the side is not parallel to the z-axis, the slope p=0.1 of the straight line where the side is located, the intercept q=0.8, the discriminant q=0.8, and the Q value being non-zero indicates that the current observation point (0, 0) is not on the same straight line as the side, so that u xz_2 and u zx_2 are calculated by formulas (10) - (13), and u zz_2 and u xx_2 are calculated by formulas (12) - (15);
For the third side of the polygon, k=3, x 3 ≠ x4, which indicates that the side is not parallel to the z-axis, the slope p= -0.1 of the straight line where the side is located, the intercept q=2.8, the discriminant q=2.8, and the Q value is not zero, which indicates that the current observation point (0, 0) is not on the same straight line as the side, so that u xz_3 and u zx_3 are calculated by using formulas (10) - (13), and u zz_3 and u xx_3 are calculated by using formulas (12) - (15);
For the fourth side of the polygon, k=4, x 4 ≠ x5, which indicates that the side is not parallel to the z-axis, the slope p= -0.6667 of the straight line where the side is located, the intercept q=9.6, the discriminant q=9.6, and the Q value is not zero, which indicates that the current observation point (0, 0) is not on the same straight line as the side, so that formulas (10) - (13) are used to calculate u xz_4 and u zx_4, and formulas (12) - (15) are used to calculate u zz_4 and u xx_4;
For the fifth side of the polygon, k=5, x 5 ≠ x6, which indicates that the side is not parallel to the z-axis, the slope p= -0.3333 of the straight line where the side is located, the intercept q=5.1, the discriminant q=5.1, and the Q value is not zero, which indicates that the current observation point (0, 0) is not on the same straight line as the side, so that formulas (10) - (13) are used to calculate u xz_5 and u zx_5, and formulas (12) - (15) are used to calculate u zz_5 and u xx_5;
for the sixth side of the polygon, k=6, x 6 ≠ x7, which indicates that the side is not parallel to the z-axis, the slope p=0 of the straight line where the side is located, the intercept q=0.1, the discriminant q=0.1, and the Q value being not zero indicates that the current observation point (0, 0) is not on the same straight line as the side, so that u xz_6 and u zx_6 are calculated by formulas (10) - (13), and u zz_6 and u xx_6 are calculated by formulas (12) - (15);
S7: adding the contribution values of the 6 sides to obtain the value of the gravity gradient tensor on the current observation point (0, 0), wherein,
Uxz(0,0)= uxz_1+uxz_2+uxz_3+uxz_4+uxz_5+uxz_6，
Uzx(0,0)= uzx_1+uzx_2+uzx_3+uzx_4+uzx_5+uzx_6，
Uzz(0,0)= uzz_1+uzz_2+uzz_3+uzz_4+uzz_5+uzz_6，
Uxx(0,0)= uxx_1+uxx_2+uxx_3+uxx_4+uxx_5+uxx_6；
S8: judging whether observation points which are not calculated exist, if so, designating the observation points as current observation points, wherein the coordinates of the observation points are (x 0,z0), and returning to the step S6 to perform new calculation; if not, executing the next step;
S9: and outputting calculated values of the gravity gradient tensor elements U xz、Uzx、Uxx、Uzz on all observation points, and drawing a curve as shown in fig. 4.
S10: to verify the validity of the calculation results of the method, the gravity gradient tensor generated by the model of fig. 3 was calculated on the same measurement points by the method of dividing the constant density units. Firstly, dividing the polygonal model in fig. 3 into 10747 constant-density small units by utilizing a grid with the length of 0.1km in the x direction and the length of 0.01km in the z direction, wherein the small units 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 density at the peak and the center point of each small cell is calculated using the density function shown in formula (16), and their average value is taken as the density value of the whole cell; finally, calculating the value of the gradient tensor generated by each polygon small unit on the measuring point by using the existing analytical formula of 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.
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, root mean square errors between the constant density method and the results of the method are recorded under different unit numbers, CPU time consumed by the two methods is recorded, and the results are shown in FIG. 5.
FIG. 2 is a graph showing the geometric relationship between the observation point and the binary. Fig. 3 is a variable density hexagonal model simulating a lobe. FIG. 4 is a comparison of the results of gravity gradient tensors generated 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 agree, which illustrates the correctness of the method. Fig. 5 (a) shows the relationship between the root mean square error and the number of divided units in the constant density method, and (b) shows the CPU time consumed by each of the two methods for different numbers of discrete units. As can be seen from fig. 5 (a), the root mean square error of the constant density method and the method herein results gradually decreases as the number of partitioning units gradually increases. It was found that the precision of the constant density method was improved and the efficiency was gradually reduced in combination with the graph (b): if the root mean square error of the calculation result of the normal density method is less than 0.001E, at least the model in fig. 3 needs to be divided into 7000 small units, and at this time, the calculation of U xz and U zx by the normal density method needs to take 20.4036s, the calculation of U xx and U zz needs 19.3334s, however, the calculation by the normal density method needs 7.5153s and 12.4313s, and the calculation time of the normal density method is 2.71 times and 1.56 times of the calculation time of the normal density method respectively. This demonstrates that the present method is superior to conventional methods in terms of both accuracy and efficiency.
The foregoing is only illustrative of the present invention and is not to be construed as limiting thereof, and all such equivalent modifications, equivalents, and improvements that fall within the spirit and principles of the invention are intended to be included within the scope of the invention.
Claims (1)
1. The method for calculating the gravity gradient tensor of the second-degree variable density body is characterized by comprising the following steps of:
Step one: establishing a plane rectangular coordinate system in a plane where the section of the binary body is positioned, wherein the horizontal axis x represents the horizontal distance and is taken to be positive to the right; the vertical axis z represents depth, oriented positive;
step two: approximating the cross section of the second degree body as a polygon with N k sides, starting from any vertex, sequentially determining the coordinates of each vertex to be (x k,zk) along the anticlockwise direction, wherein the polygon is a closed graph, and the N k +1st vertex is the 1st vertex, wherein k=1, 2, … and N k,Nk +1;
Step three: representing density variation of the dyadic with a binary polynomial function , determining a highest order N x of the x-coordinate and a highest order N z of the z-coordinate, and coefficients D i,j of the polynomial;
Step four: setting a plurality of observation points along the observation line and determining the coordinates of the observation points, wherein the abscissa of the observation points cannot be equal to the abscissa of the polygon vertexes;
Step five: the current observation point is specified, and the coordinates thereof are (x 0,z0),x0≠xk;
Step six: respectively calculating the contribution value U k of the gravity gradient tensor U corresponding to the kth edge of the polygon on the current observation point, wherein k=1, 2, …, N k,
Step seven: adding the contribution values corresponding to the N k sides to obtain the value of the gravity gradient tensor on the current observation point, wherein,
Step eight: judging whether observation points which are not calculated exist or not, 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 calculated values of the gravity gradient tensor elements U xz、Uzx、Uxx、Uzz on all observation points;
in the sixth step, a formula for calculating the gravity gradient tensor contribution value corresponding to the kth edge is as follows:
① When the kth edge is parallel to the z-axis,
In the formulas (1) - (2), G is a gravitational constant, that is, g= 6.67408 ×10 -11N×m2/kg2;Nx and N z are the highest orders of the x-coordinate and the z-coordinate in the density function, i and j are orders corresponding to the x-coordinate and the z-coordinate in the non-zero term D i,jxizj, and D i,j is a coefficient of the density polynomial; n and m are also orders of x and z coordinates: n is an order between 0 and i, m is an order between 0 and j; c i n and C j m are combined numbers; (x 0,z0) is the coordinates of the current observation point; functions E xz and E zz represent two piecewise functions related to i, j, n, m:
In the above equations (3) - (4), (x k,zk) is the coordinates of the kth vertex, (x k+1,zk+1) is the coordinates of the (k+1) th vertex, r k and r k+1 respectively represent the distances between the kth vertex and the (k+1) th vertex to the current observation point, and rk+1＝[(xk+1-x0)2+(zk+1-z0)2]1/2,rk＝[(xk-x0)2+(zk-z0)2]1/2;H represents a recursive sequence, and if the sequence number of the entries in this sequence is denoted by l, the general formula of the recursive sequence H l is:
K in the above formula (5) represents another recursive sequence, and if s represents the sequence number of the term in this sequence, the general formula of the recursive sequence K s is:
② When the kth edge is not parallel to the z-axis and the current observation point is on the same line as the kth edge,
In the formulas (7) - (8), G is a gravitational constant, that is, g= 6.67408 ×10 -11N×m2/kg2;Nx and N z are the highest orders of the x-coordinate and the z-coordinate in the density function, i and j are orders corresponding to the x-coordinate and the z-coordinate in the non-zero term D i,jxizj, and D i,j is a coefficient of the density polynomial; (x 0,z0) is the coordinates of the current observation point; and (x k,zk) is the coordinates of the kth vertex, (x k+1,zk+1) is the coordinates of the (k+1) th vertex, p and q are the slope and intercept ,p＝(zk+1-zk)/(xk+1-xk),q＝(zkxk+1-zk+1xk)/(xk+1-xk); of the straight line where the kth edge is located, respectively, and the calculation formula of the number sequence M is:
In the formula (9), c= 1+p 2,d＝[p(q-z0)-x0 ]/c;
③ When the kth edge is not parallel to the z-axis and the current observation point is not on the same straight line as the kth edge,
In the formulas (10) - (11), G is a gravitational constant, that is, g= 6.67408 ×10 -11N×m2/kg2;Nx and N z are the highest orders of the x-coordinate and z-coordinate in the density function, i and j are orders corresponding to x and z in the density polynomial, and D i,j is a coefficient of the density polynomial; n and m are also orders of x and z coordinates: n is an order between 0 and i, m is an order between 0 and j; c i n and C j m are combined numbers; (x 0,z0) is the coordinates of the current observation point; the functions F xz and F zz represent two piecewise functions related to i, j, n, m:
In the formulas (12) - (13), the (x k,zk) is the coordinate of the kth vertex, (x k+1,zk+1) is the coordinate of the kth+1th vertex, p and q are the slope and the intercept ,p＝(zk+1-zk)/(xk+1-xk),q＝(zkxk+1-zk+1xk)/(xk+1-xk);Q＝px0+q-z0;l1 and l 2 of the straight line where the kth edge is located respectively, and the order is represented, but the value ranges of the two are different according to the different formulas, and and/> are both combination numbers; i and J represent two arrays, respectively, wherein the calculation formula of the array I is as follows:
In the formula (14), q=px 0+q-z0,b＝p(q-z0)-x0,c＝1+p2,rk and r k+1 represent the distances between the kth vertex and the (k+1) th vertex and the current observation point, respectively; the other column J is a recursive column, and if the sequence number of the term in this column is denoted by l, the general formula of the recursive column J l is:
Q＝px0+q-z0,a＝(q-z0)2+x0 2,b＝p(q-z0)-x0,c＝1+p2,rk and r k+1 in the above formula (15) represent the distances between the kth vertex and the (k+1) th vertex to the current observation point, respectively; the calculation formula of I 0 is shown as formula (14).
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 CN111950108A (en) | 2020-11-17 |
| CN111950108B true 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 (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB2468921B (en) * | 2009-03-27 | 2013-06-05 | Qinetiq Ltd | Method and system for the detection of gravitational anomalies |
-
2019
- 2019-05-14 CN CN201910398005.7A patent/CN111950108B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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)
| Title |
|---|
| Analytic Expressions for the Gravity Gradient Tensor of 3D Prisms with Depth-Dependent Density;Li Jiang et al;Surv Geophys;1-27 * |
| Gravity anomalies of 2-D bodies with variable density contrast;Jianzhong Zhang et al;GEOPHYSICS;第66卷(第03期);809-813 * |
| Potential fields and their partial derivatives produced by a 2D homogeneous polygonal source: A summary with some revisions;Zhen Jia et al;GEOPHYSICS;第76卷(第04期);L29-L34 * |
| 密度随深度变化的二度体重力正演公式;张建中 等;石油地球物理勘探;第35卷(第02期);202-207 * |
| 岩石密度空间变化时重力效应的计算方法;郭武林;物探与化探(第04期);226-236 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111950108A (en) | 2020-11-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106646645B (en) | A kind of gravity forward modeling accelerated method | |
| CN110045432B (en) | Gravity field forward modeling method and three-dimensional inversion method under spherical coordinate system based on 3D-GLQ | |
| Li et al. | Three-dimensional gravity modeling in all space | |
| CN109141426B (en) | Method for matching navigation adaptation area by underwater gravity | |
| CN109375280A (en) | Gravitational field quick high accuracy forward modeling method under a kind of spherical coordinate system | |
| CN106777598B (en) | Numerical simulation method for magnetic field gradient tensor of complex magnetic body with arbitrary magnetic susceptibility distribution | |
| CN109254327B (en) | Exploration method and exploration system of three-dimensional ferromagnetic body | |
| CN111158059B (en) | Gravity inversion method based on cubic B spline function | |
| WO2012001388A2 (en) | Gravity survey data processing | |
| CN115292973A (en) | Arbitrarily sampled space wave number domain three-dimensional magnetic field numerical simulation method and system | |
| CN113962077A (en) | Three-dimensional anisotropic strong magnetic field numerical simulation method, device, equipment and medium | |
| CN114943178A (en) | Three-dimensional geological model modeling method and device and computer equipment | |
| Muratov et al. | Grid-characteristic method on unstructured tetrahedral meshes | |
| CN113109883B (en) | Method for forward modeling satellite gravity field based on isoparametric transformation global discrete grid spherical coordinates | |
| CN107942399A (en) | One kind is greatly apart from potential field upward continuation computational methods | |
| CN109444973B (en) | Gravity forward modeling acceleration method under spherical coordinate system | |
| CN106662665B (en) | The interpolation and convolution of rearrangement for the processing of faster staggered-mesh | |
| CN111665556A (en) | Method for constructing stratum acoustic wave propagation velocity model | |
| CN111950108B (en) | Method for calculating gravity gradient tensor of second-degree variable density body | |
| Li et al. | Fast 3D forward modeling of the magnetic field and gradient tensor on an undulated surface | |
| CN107748834A (en) | A kind of quick, high resolution numerical simulation method for calculating fluctuating inspection surface magnetic field | |
| CN115640720B (en) | Self-attraction simulation method based on distance control grid encryption | |
| CN113267287B (en) | Method for reconstructing shock wave overpressure three-dimensional space-time field | |
| CN115270579A (en) | Second-order acoustic wave equation finite difference numerical simulation parameter selection method | |
| Yin et al. | A fast 3D gravity forward algorithm based on circular convolution |
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 |