Method and device for forward calculation of Bouguer gravity anomaly using variable density of geological information

Abstract

The present invention belongs to the field of gravity exploration technology, and discloses a method for calculating Bouguer gravity anomaly by forward modeling using geological information variable density. The geological information of the explored area and high-resolution DEM data are used to uniformly forward model the influence value of the intermediate layer material to obtain a fine Bouguer gravity anomaly value, which replaces the traditional Bouguer gravity anomaly value obtained by step-by-step calculation after terrain correction and intermediate layer correction, thereby improving the accuracy of the Bouguer gravity anomaly; the three-dimensional geological body variable density calculation is used, and the deep stratum occurrence is adopted as any angle of field geological measurement or drilling control, which fully utilizes the geological information of the explored area, is more in line with the actual geological facts, and the finely calculated Bouguer gravity anomaly is convenient for effectively distinguishing the anomalies caused by the target body and the interference body, that is, it is convenient to distinguish the regional gravity anomaly from the local gravity anomaly, and can achieve a good gravity exploration effect. The present invention fills the gap in the fine interpretation of large-scale gravity exploration.

Description

Technical Field

[0001] The invention belongs to the technical field of gravity exploration, and in particular relates to a method for forward calculating Bouguer gravity anomaly by using geological information variable density. Background Art

[0002] At present, gravity exploration is very effective in metal ore surveys such as mineralization prediction, delineation of rock mass and structure, and research on the characteristics of metal ore occurrence, and has good application prospects. The object of gravity exploration is the gravity acceleration value at the measuring point. Generally, the relative measurement method is used to measure the step difference between the measuring point and the gravity base point. The gravity acceleration value of each measuring point is obtained by using the gravity acceleration value of the gravity base point and the step difference between the base point and the measuring point, which is referred to as the gravity value.

[0003] The traditional Bouguer gravity anomaly value is obtained by correcting the gravity value through normal field correction (latitude correction), intermediate layer correction, height correction (intermediate layer correction and height correction are collectively called Bouguer correction) and terrain correction. It is the basic data of gravity exploration. The calculation steps and physical meaning of the traditional Bouguer gravity anomaly are as follows: Figure 1 .

[0004] like Figure 1 , the measured gravity value at any point on the ground depends on the latitude, altitude, changes in the surrounding terrain and the uneven distribution of local underground materials, etc. The last factor is an important research object of gravity exploration. The gravity changes caused by other factors are interference. Therefore, the calculation of the Bouguer gravity anomaly at the measuring point is the process of stripping off the interference factors. Assuming that the earth is a smooth ellipsoid with layered density distribution, the density is uniform in the same layer, and the interfaces of each layer are also confocal rotating ellipsoids, then the gravity value of each point on the ellipsoid surface can be calculated based on the earth's gravitational constant, major radius, flattening and rotation angular velocity. The smooth ellipsoid is called a normal ellipsoid, the normal ellipsoid surface is the geoid, and the gravity value on the normal ellipsoid surface is called the normal gravity value.

[0005] Figure 1 a: The measured gravity values ​​of P and Q measuring points are composed of the normal gravity values ​​of the layers below the geoid, the gravity values ​​caused by the materials above the geoid to the terrain undulations, and the gravity values ​​caused by the different heights of the measuring points. The normal ellipsoid is divided into two density interfaces by the normal crust thickness plane. The theoretical density (black words) between the geoid and the normal crust thickness plane is 2.67g / cm 3 The theoretical density below the normal crust thickness plane is 3.27 g / cm 3 There are two types of uneven distribution of materials in the above theoretical model. One is that the actual lower crust interface is the Moho surface. The Moho surface is undulating and generally does not coincide with the normal crust thickness plane. The theoretical density of the crust above the Moho surface (white text) is 2.67g / cm 3Therefore, there is an uneven distribution of material between the normal crust thickness plane and the Moho surface; secondly, there are various geological bodies with local actual density lower or higher than the average density of the crust distributed at different depths, which are the key research objects of our gravity exploration.

[0006] Figure 1 b: Normal gravity correction, that is, eliminating the gravity changes caused by the layers of material below the geoid. After elimination, the gravity values ​​of each ground point are the residual gravity values ​​of the material between the Moho surface and the normal shell thickness plane (the residual density difference is -0.6g / cm 3 ), the gravity value caused by the material between the geoid and the terrain undulation surface, and the gravity value caused by the local uneven density geological body. The calculation formula of the normal gravity value of the CGCS2000 ellipsoid is as follows:

[0007]

[0008] Where: g0 is the normal gravity value, is the latitude of the measuring point.

[0009] Figure 1 c in the formula: Terrain correction means eliminating the change in gravity value caused by terrain fluctuation within a certain range around the measuring point, taking the measuring point elevation surface as the reference surface. For example, areas 1, 2, and 3 around point P and areas 4 and 5 around point Q, where areas 1, 3, and 4 are material deficits and areas 2 and 5 are material surpluses. After terrain correction, the gravity value of each ground measuring point is corrected by the residual gravity value ( Figure 1 b) The gravity value caused by the regular layered material between the geoid and the elevation plane of the measuring point and the gravity value caused by the local uneven density geological body. Figure 2 .

[0010] like Figure 2 , take the rectangular coordinate system XYZ, set the position of the measuring point A as the origin, the Z axis is vertically downward, the X and Y axes are in the horizontal plane where point A is located, dm is the mass element, and its coordinates are (x, y, z). The radius vector from point A to the mass element is represented by r, and the angle between r and the Z axis is θ. The gravitational influence generated by dm at point A, that is, the vertical component of gravity dg is:

[0011]

[0012] If the influence value Δgt of the total mass of area (1) or (2) is calculated, it can be calculated by integration:

[0013]

[0014] Where: G is the gravitational constant, V is the geological volume, (x, y, z) is the coordinates of the center of mass of the geological body, and ρ is the residual density of the geological body.

[0015] In the mass surplus area (such as area (1)), z is negative and ρ is positive. In the mass deficit area (such as area (2)), z is positive and ρ is negative. Therefore, the gravity influence of the terrain on the measuring point is always negative, that is, the terrain undulation reduces the gravity value.

[0016] Figure 1 d in the figure: intermediate layer correction and height correction. Intermediate layer correction is to eliminate the gravity value caused by the regular layered material between the geoid and the elevation plane of the measuring point, and height correction is to eliminate the difference in the normal gravity field with the height change due to the different elevations of the measuring points. The eliminated gravity value is called the Bouguer gravity anomaly value, which includes the regional gravity anomaly caused by the fluctuation of the Moho surface and the local gravity anomaly caused by the local uneven density geological bodies.

[0017] The calculation formula of the intermediate layer correction value Δgm is as follows:

[0018]

[0019] In the formula: ρ is the density, h is the thickness of the intermediate layer, that is, the distance from the elevation plane of the measuring point in Figure d to the geoid, and a is the correction radius of the intermediate layer disk, which is generally equal to the terrain correction radius.

[0020] The calculation formula for the altitude correction Δgh is:

[0021]

[0022] Where: h is the distance from the elevation plane of the measuring point to the geoid, is the latitude of the measuring point.

[0023] Calculation of Bouguer gravity anomaly:

[0024] Δgb=G 观 -g0+Δgh+Δgm+Δgt (6)

[0025] Where: Δgb is the Bouguer gravity anomaly, G 观 is the measured gravity value at the measuring point, g0 is the normal gravity value, Δgh is the height correction value, Δgm is the middle layer correction value, and Δgt is the terrain correction value.

[0026] The above-mentioned traditional Bouguer gravity anomaly calculation method is suitable for the calculation and processing of regional gravity exploration at small and medium scales. If it is used for large-scale gravity exploration or gravity exploration at the mining area scale, there are four obvious defects in the intermediate layer correction and terrain correction:

[0027] 1. The density of traditional intermediate layer correction and terrain correction is a unified value of 2.67g / cm 3, which obviously does not conform to the fact of geological body and density change on the plane. The Bouguer gravity anomaly calculated based on this, the local anomaly information contained in it, is easily submerged by the gravity anomaly caused by density change, and cannot effectively reflect the actual gravity anomaly of the target geological body.

[0028] 2. If there is no vertical stratification calculation for the intermediate layer correction, the geological information of the explored area will not be fully utilized. With the change of the occurrence of strata with different densities in the vertical direction, the density of deep geological bodies has changed. At the same time, a large number of drill cores or mined ores in the mining area can provide a lot of key density information. If the intermediate layer is calculated as a single layer, it will not be conducive to effectively distinguishing between the known geological bodies and the newly discovered geological bodies, and often cannot achieve good gravity exploration results.

[0029] 3. The density of the missing material in the terrain correction is unknown, but is only 2.67 g / cm 3 Calculation using this theoretical density value may lead to errors in terrain correction, which is not conducive to the precise solution of the Bouguer gravity anomaly, and further affects the precise division of regional gravity anomalies and local gravity anomalies.

[0030] 4. If the terrain correction adopts the square domain calculation formula, and the intermediate layer theoretical calculation formula is a disk, there will be inconsistency in the correction range between the terrain correction and the intermediate layer correction, resulting in a loss of accuracy of the Bouguer gravity anomaly.

[0031] Through the above analysis, the problems and defects of the existing technology are as follows: the existing gravity anomaly calculation results have low accuracy and large errors, and sometimes cannot effectively reflect the actual gravity anomaly caused by the target geological body, cannot provide accurate information for large-scale gravity exploration, and the gravity exploration effect is not good.

[0032] At present, some 3D modeling software, such as Geomodel, can perform forward calculation of geological gravity by using various geological information through geological modeling. However, there are still two major problems: First, the forward calculation still uses 2.67g / cm 3 This unified density cannot perform variable density calculations and needs to be improved. Second, the software is complex to operate and expensive, making it difficult to popularize and promote. Summary of the invention

[0033] In view of the problems existing in the prior art, the present invention provides a method for forward calculating the Bouguer gravity anomaly by using geological information variable density.

[0034] The present invention is implemented as follows: a method for forward calculating Bouguer gravity anomaly using geological information variable density, the method for forward calculating Bouguer gravity anomaly using geological information variable density comprises:

[0035] Using the geological information of the explored area, including the boundary, occurrence, depth of the geological body in the geological map or the density change of the geological body vertically in the borehole, according to the three-dimensional distribution state of the geological body, not only block in the plane, but also layer in the vertical direction, and assign the measured density value to different geological bodies to calculate the variable density;

[0036] The high-resolution DEM data is used to uniformly forward calculate the influence value of the intermediate layer material instead of the step-by-step calculation of terrain correction and intermediate layer correction, so as to improve the accuracy of Bouguer gravity anomaly.

[0037] The variable density calculation of three-dimensional geological bodies is used, and the deep strata occurrence is taken as any angle of field geological measurement or drilling control.

[0038] Furthermore, the method specifically comprises:

[0039] The measured gravity values ​​are corrected for the normal field and the height to obtain the gFI-corrected gravity values. According to the density information of various planes and deep parts such as geological maps and boreholes, the geological bodies are layered and divided into blocks in the plane and vertical directions (vertically layered to the geoid, i.e., the horizontal plane with an altitude of 0) using the DEM grid data. The influence values ​​of the layered and divided geological bodies on a certain gravity measuring point are forward calculated using the right cuboid formula. The influence values ​​of all geological bodies are summed to obtain the intermediate material correction value of the gravity value of the measuring point from the material between the geoid and the terrain undulation surface. The above-mentioned gFI-corrected gravity value minus the intermediate material correction value is the precisely calculated Bouguer gravity anomaly value.

[0040] Furthermore, the method for calculating the Bouguer gravity anomaly by forward modeling the geological information variable density of the explored area comprises the following steps:

[0041] Step 1: normal field and height correction are performed on the gravity value of the measuring point; the DEM grid spacing is calculated, and the grid data is vertically layered. After the vertical layering, the grid point is taken as the center, and half the grid spacing is extended in all directions on the plane to form the top surface of a rectangular parallelepiped. The rectangular parallelepiped is constructed by the vertical layering and the top surface, and the elevation of the center point of the rectangular parallelepiped and the half height of the rectangular parallelepiped are calculated;

[0042] Step 2: prepare a three-dimensional geological body (a geological body has a unique density value, so it is also called a density body) interface file using the measured density value; segment the grid layered data with each three-dimensional density body, and obtain the data file for forward calculation belonging to each geological body after segmentation, including the grid spacing, the plane coordinates of the rectangular parallelepiped, the elevation of the center point of the rectangular parallelepiped, the half height of the rectangular parallelepiped, and the measured density value of the rectangular parallelepiped;

[0043] Step 3: Use the variable density forward modeling of the rectangular parallelepiped formula to obtain the impact value of the middle layer material; and perform fine Bouguer gravity anomaly calculation.

[0044] Furthermore, the normal field and height correction of the gravity value of the measuring point in step 1 includes:

[0045] The gravity data of the measuring point includes the three-dimensional coordinates of the measuring point x, y, h and G 观 Four columns of measured gravity values ​​are used to perform normal field correction and height correction to obtain the three-dimensional coordinates of the measuring point x, y, h and the gFI corrected gravity value.

[0046] Furthermore, the stratification of the grid data in step 1 includes:

[0047] The grid files are layered at a certain vertical interval. Based on the forward calculation of the rate of change of gravity value with depth under fixed grid spacing, the DEM elevation of the entire area is divided into two large layers, 400m below the lowest point of the entire area:

[0048] The upper large layer is continuously divided into layers at intervals of 20m to the highest point of the entire area, and the large layer discards the grid points whose actual surface grid elevation is lower than the elevation of the divided layer;

[0049] The next large layer is further divided into two layers at 50m intervals according to the lowest elevation of 400-500m, 500-700m is divided into two layers at 100m intervals, 700-1100m is divided into two layers at 200m intervals, 1100-1500m is divided into one layer at 400m intervals, and 1500-2300m is divided into one layer at 800m intervals. If the lowest point elevation of the DEM grid in the exploration area is greater than 2700m, then 2300m to 0m above sea level is divided into one layer according to the actual remaining interval. If the lowest point elevation is less than 2700m, then 1500m to 0m above sea level is divided into one layer according to the actual remaining interval. All grid points in the large layer are retained;

[0050] Furthermore, because the DEM data is a regular grid, that is, the grid spacing is a unique value, the grid spacing described in step 1 is equal to the distance between two adjacent grid points;

[0051] Furthermore, the elevation of the center point of the grid after stratification and the half-height of the rectangular block described in step 1 are calculated using the stratified elevation values. The elevation of each center point is the sum of the elevation values ​​of the upper and lower layers (if the actual elevation value of the DEM grid point is less than the elevation value of the upper layer, the actual elevation value of the DEM grid point is taken as the elevation value of the upper layer, the same below), and the half-height of each grid cube is the difference between the upper and lower layer elevation values ​​divided by 2.

[0052] Further, the step 2 of preparing the geological interface file using the measured density value includes:

[0053] (1) According to the geological map of the gravity exploration area or other boreholes collected and measured by lithology, the statistical density value of each geological body is obtained; the statistical density value of each geological body includes: the name of the geological body and the statistical density value;

[0054] (2) Prepare geological interface files:

[0055] Simplify the boundary of the geological body, count the interface information of the geological body in the area, enclose the geological body with multiple interfaces, and use the enclosed geological body to split the grid data layered file completed in step 1. The split file matches the measured density value of the enclosed geological body, thus completing the density body preparation.

[0056] The interface information includes: the coordinates of the geological boundary position, the occurrence or the drilling position, and the depth are digitized to form multiple interface information of a geological body. The interface information includes: the x coordinates of the geological boundary point 1 and the x coordinates of the geological boundary point 2 1 ,y 1 、z 1 、x 2 ,y 2 、z 2 , inclination angle A, the relationship between body and surface J, the addition and subtraction relationship between strike and dip.

[0057] After preparation, the data used for forward modeling include: X, Y coordinates of grid points, elevation of the center point of the rectangular parallelepiped, half height of the rectangular parallelepiped, and measured density value.

[0058] Furthermore, the step 3 of using the right cuboid forward modeling to obtain the influence value of the intermediate layer material includes:

[0059] Based on the obtained geological body data file for the forward calculation of the rectangular parallelepiped, combined with the x, y, and h in the gravity value file of the measuring point, the rectangular parallelepiped forward modeling formula is used to perform the forward calculation of a single rectangular parallelepiped in each geological body data file for any measuring point, and the forward modeling values ​​of all rectangular parallelepipeds are added together to obtain the influence value of the geological body on the intermediate layer material of a certain gravity measuring point; the x, y, h, and g of the measuring point are obtained. 中 .

[0060] Further, the step 3 of calculating the fine Bouguer gravity anomaly includes:

[0061] The calculated corrected gravity value gFI is compared with the obtained g 中 Subtract them to get the fine Bouguer gravity anomaly value of the gravity measuring point.

[0062] Another object of the present invention is to provide a device for forward calculating Bouguer gravity anomaly using geological information variable density, the device comprising: an acquisition module, which acquires geological information (geological map and measured density value of related geological body) of the explored area, three-dimensional coordinates of gravity measuring points, measured gravity values ​​of gravity measuring points and DEM grid data; a first determination module, which is used to perform normal field correction and height correction according to the three-dimensional coordinates and measured gravity values ​​of the gravity measuring points to obtain gFI corrected gravity values; a second determination module, which is used to vertically layer the DEM grid data according to the principle of the rate of change of a certain right rectangular parallelepiped with the distance from the measuring point (i.e. depth) and size (grid spacing is constant, i.e. rectangular parallelepiped height, layer spacing) to form a grid layer file; a third determination module, which is used to organize the simplified and organized geological information into a geological interface file, and according to the interface file, the above-mentioned grid layer file is closed and divided, and the measured density of the geological body is matched to finally form a right rectangular parallelepiped forward modeling data file; a fourth determination module, which is used to forward calculate the influence value of the intermediate layer material; and a fifth determination module, which is used to calculate the refined Bouguer gravity anomaly value.

[0063] Another object of the present invention is to provide a computer device, comprising a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of the method for forward calculating the Bouguer gravity anomaly using variable density of geological information.

[0064] In combination with the above technical solutions and the technical problems solved, the advantages and positive effects of the technical solutions to be protected by the present invention are analyzed from the following aspects:

[0065] The present invention can make full use of the geological information of the explored area, including but not limited to the boundary, occurrence, depth of the geological body in the geological map or the density change of the geological body vertically in the borehole, that is, all the geological information that has been mastered; in the terrain correction and intermediate layer correction in the traditional processing, all geological bodies are assumed to be a single density body on the plane, and the occurrence of the geological body is assumed to be two simple and special models of upright or horizontal in the vertical direction. The collected and statistical physical properties only play a qualitative relative role in the interpretation of anomalies. The actual geological situation is that the density of the geological bodies divided on the plane is very different and cannot be calculated and simulated with a unified density. The occurrence of the geological body in the vertical direction is an arbitrary angle, which cannot be replaced by the two special situations of upright or horizontal. Therefore, when the density of the relevant geological bodies has been collected and statistically analyzed, the geological bodies should be divided into blocks on the plane and layers in the vertical direction according to the three-dimensional distribution state of the geological bodies. Different geological bodies are given their measured density values ​​to perform variable density calculation, so as to distinguish and extract effective anomalies and improve the exploration accuracy and effect.

[0066] The present invention adopts high-resolution DEM data to uniformly forward calculate the influence value of the intermediate layer material, replacing the step-by-step calculation of terrain correction and intermediate layer correction, thereby improving the accuracy of Bouguer gravity anomaly;

[0067] The present invention adopts variable density calculation of three-dimensional geological body, which is different from the traditional variable density terrain correction. The traditional variable density terrain correction calculation only considers the change of density of planar geological body, that is, it is considered as vertical strata (steeply inclined 90°) in the vertical direction; while the present invention also integrates the change of deep geological body, that is, the occurrence of deep strata is any angle measured by field geology or controlled by drilling, which is more in line with geological reality.

[0068] The present invention can improve the defects of existing software that cannot fully utilize geological information (occurrence, density changes of different geological bodies) and use unified density calculation, and the Bouguer gravity anomaly calculated in detail is more in line with the actual geological situation.

[0069] The expected benefits and commercial value after the technical solution of the present invention is converted as follows: the present invention is expected to be converted into software for the precise calculation of Bouguer gravity anomalies. The current situation of mineral resources is severe, and large-scale gravity exploration work at the mining area scale is imperative. The method of the present invention is a necessary means to improve exploration accuracy, improve exploration results, and discover deep and marginal mineral resources. After the invention is converted into software, it can be sold or provided with technical services to geological exploration units, and its commercial value is obvious.

[0070] The technical solution of the present invention fills the technical gap in the industry at home and abroad: the present invention fills the method gap of making full use of geological information variable density forward modeling to perform fine Bouguer gravity anomaly. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 Schematic diagram of the calculation process and physical meaning of the traditional Bouguer gravity anomaly provided by an embodiment of the present invention; in the figure: P and Q are surface gravity measurement points;

[0072] Figure 2 is a schematic diagram of terrain influence provided by an embodiment of the present invention;

[0073] Figure 3 It is a flow chart of a method for forward calculation of Bouguer gravity anomaly using geological information variable density provided by an embodiment of the present invention;

[0074] Figure 4 is a schematic diagram of the relationship between gravity measurement points and terrain grid nodes provided by an embodiment of the present invention;

[0075] Figure 5 is a schematic diagram of calculation of a rectangular parallelepiped model provided by an embodiment of the present invention;

[0076] Figure 6It is a flow chart of a method for forward calculation of Bouguer gravity anomaly using right cuboid variable density provided by an embodiment of the present invention;

[0077] Figure 7 It is a gravity profile and a topographic map around a measuring point provided by an embodiment of the present invention;

[0078] Figure 8 It is a gravity profile and a simplified geological map around the measuring point provided by an embodiment of the present invention;

[0079] Fig. 9 It is a schematic diagram of the occurrence change of various geological bodies on the cross section provided by an embodiment of the present invention;

[0080] Fig.10 It is a schematic diagram of the distribution of grid data at an elevation of 1480-1500m provided by an embodiment of the present invention;

[0081] Fig.11 is a schematic diagram of a vertical rectangular parallelepiped of a ΣⅢ geological body provided in an embodiment of the present invention;

[0082] Fig.12 It is a schematic diagram comparing the exploration effects of the traditional method and the improved method provided in the embodiment of the present invention. DETAILED DESCRIPTION

[0083] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0084] The embodiment of the present invention provides a method for calculating Bouguer gravity anomaly by using geological information variable density forward modeling, the method comprising: obtaining geological information of the explored area (geological map and measured density value of related geological body), three-dimensional coordinates of gravity measuring points, measured gravity values ​​of gravity measuring points and DEM grid data (X, Y, Z1); using the three-dimensional coordinates of the gravity measuring points and the measured gravity values, completing normal field and height correction through formulas (1) and (5) to obtain the corrected gravity value gFI; calculating the grid spacing of the DEM grid data, stratifying the DEM grid data according to the above technical solution, and obtaining the vertical longitude after stratification. The cube data file is used to calculate the center point elevation Zm of the rectangular parallelepiped and the half height c of the rectangular parallelepiped at the same time; the geological information is used to establish the geological interface file, and the established multiple geological interfaces are used to segment the rectangular parallelepiped data file to form a rectangular parallelepiped file of a certain geological body, and the measured density of the geological body is matched, and the rectangular parallelepiped file segmentation and measured density matching of all geological bodies are completed in a cycle, that is, the forward modeling data preparation is completed; the three-dimensional coordinates of the gravity measurement points, the grid spacing and the forward modeling data are used to calculate the influence value of the intermediate layer material according to the rectangular parallelepiped forward modeling formula to obtain the intermediate layer material correction value g 中; Compare the calculated corrected gravity value gFI with the obtained g 中 Subtract them to get the fine Bouguer gravity anomaly value of the gravity measuring point.

[0085] The method of obtaining the gFI gravity correction value by normal field correction and altitude correction of the three-dimensional coordinates of the gravity measuring point and the measured gravity value in the first step of this application is a mature solution, that is, there is an existing formula for direct correction, so it will not be described in detail.

[0086] The grid data stratification scheme in the second step of this application is the first innovation of the present invention. First, the density is based on the grid spacing of the DEM data (the grid spacing is the length and width of the cuboid) and the theoretical density is 2.67 g / cm 3 , determine the rate of change of the influence value corresponding to the grid spacing with depth at a certain depth interval. For example, the lowest elevation grid point is subtracted by 400m to divide it into two large layers, and the center point height Z and half height c of the changing cuboid are calculated according to formula (3). The influence value of the cuboid with different heights and half heights on the measuring point directly above is calculated. The influence value is less than or equal to 0.001×10 -5 m / s 2 The center point height and the half height of the cuboid used in the stratification are determined based on the above. For example, the 400m below the lowest elevation grid point is divided into two large layers, and the upper large layer is stratified at intervals of 20m. That is, the influence value of the vertical cuboid with a length, width and height of 10m×10m×20m on the measuring point below a depth of 400m (the height of the center point of the cuboid corresponding to the lowest elevation grid point is the height of the grid point minus 400 / 2) is less than or equal to 0.001×10 -5 m / s 2 ; That is, the principle followed by the stratification is that the stratification depth and stratification interval are determined by the calculation test of the influence value change rate according to the grid spacing. The larger the grid spacing, the smaller the stratification depth and the finer the stratification interval. It is not limited to the stratification data used in the above scheme of this application; when stratifying, the grid points in the previous large layer whose elevation is lower than the lowest stratification elevation will be discarded, and all the grid points in the next large layer will be retained. After completion, a unique undivided upright cuboid data file is formed, including four columns of data of grid point coordinates X and Y, the elevation Zm of the midpoint of the cuboid, and the half-height c of the cuboid, among which the grid point coordinates X and Y are repeated many times because each layer corresponds to these grid points.

[0087] The second major innovation of the present invention is to prepare the geological interface file in the third step, segment the rectangular parallelepiped files belonging to the enclosed geological body according to the geological interface and assign density values ​​to form the forward modeling data file. The geological boundaries in the geological map are simplified into a series of straight lines, and the end point coordinates x1, y1, z1, x2, y2, z2 of the straight lines are extracted, as well as the inclination A of the geological interface, the relationship J between the geological body and the surface (such as 1 on the surface and -1 below the surface) and the addition and subtraction relationship between the strike and the dip (such as 1 counterclockwise and -1 clockwise), and the rectangular parallelepiped files completed in the second step are segmented by multiple geological interface (i.e., multiple lines of geological interface information) data to form all the rectangular parallelepiped files that constitute a three-dimensional geological body, and each rectangular parallelepiped file is assigned the measured statistical density of the geological body; the rectangular parallelepiped files of all geological bodies are completed in a loop to form the forward modeling data file.

[0088] In the fourth step of the present application, the intermediate layer influence value is forward modeled according to the variable density formula of the right cuboid. Its calculation formula is an existing formula, and the idea of ​​using forward modeling to replace the intermediate layer correction and terrain correction in the traditional Bouguer gravity anomaly calculation is the third major innovation of the present invention. Its significance lies in that, first, it solves the defects of the traditional intermediate layer and terrain correction that the variable density calculation cannot be performed and the geological information cannot be fully utilized; second, the density value of the geological body in the "deficient" part of the terrain correction does not actually exist, which is an ideal state, and will inevitably lead to a loss of accuracy, while the forward modeling of the present invention is all carried out using the measured density value, which is more in line with geological facts; third, the traditional intermediate layer calculation formula is a disk model, and the terrain correction is a square domain model. There is a loss of accuracy in the connection between the circular domain and the square domain, while the forward modeling of the present invention unifies the two, and there is no problem of connection between the intermediate layer and the terrain correction, and the accuracy is significantly improved.

[0089] The fifth step of this application is a simple addition and subtraction method, that is, the normal field and height-corrected gFI minus the intermediate layer material impact value gFI to obtain a fine Bouguer gravity anomaly value. Fine calculation of Bouguer gravity anomalies facilitates effective distinction between useful gravity anomalies and interference anomalies, that is, distinguishing anomalies caused by target bodies and interference bodies, improving the effect of gravity exploration, and filling the gap in fine interpretation of large-scale gravity exploration.

[0090] Second, an embodiment of the present invention provides a device for forward calculating Bouguer gravity anomaly using geological information variable density, and the device includes: an acquisition module, which acquires geological information (geological map and measured density value of related geological body) of the explored area, three-dimensional coordinates of gravity measuring points, measured gravity values ​​of gravity measuring points and DEM grid data; a first determination module, which is used to perform normal field correction and height correction according to the three-dimensional coordinates and measured gravity values ​​of the gravity measuring points to obtain gFI corrected gravity values; a second determination module, which is used to vertically layer the DEM grid data according to the principle of the rate of change of a certain right rectangular parallelepiped with the distance from the measuring point (i.e. depth) and size (grid spacing is constant, i.e. rectangular parallelepiped height, layer spacing) to form a grid layer file; a third determination module, which is used to organize the simplified and organized geological information into a geological interface file, and according to the interface file, the above-mentioned grid layer file is closed and divided, and the measured density of the geological body is matched, and finally a right rectangular parallelepiped forward modeling data file is formed; a fourth determination module, which is used to forward calculate the influence value of the intermediate layer material; and a fifth determination module, which is used to calculate the refined Bouguer gravity anomaly value.

[0091] The present invention can improve the defect that the existing software cannot fully utilize the geological information (occurrence, density changes of different geological bodies) to unify the density and calculate the traditional Bouguer gravity anomaly by terrain correction and intermediate layer correction. The finely calculated Bouguer gravity anomaly is more in line with the actual geological conditions.

[0092] like Figure 3 As shown, the method for calculating Bouguer gravity anomaly by forward modeling using geological information variable density provided by an embodiment of the present invention comprises the following steps:

[0093] S101, the gravity value of the measuring point is corrected for normal field and height; the spacing of the DEM grid data is calculated, the layering depth and layering interval are determined experimentally, and the DEM grid data is layered, that is, the rectangular parallelepiped data is obtained, and the elevation of the center point of the rectangular parallelepiped and the half height of the rectangular parallelepiped are calculated;

[0094] S102, using geological information to prepare a geological interface file; using the geological interface to trap and segment the rectangular parallelepiped file after grid stratification, obtaining a rectangular parallelepiped file of a certain geological body, matching the density of the geological body, cyclically completing the trapping and segmentation of the rectangular parallelepiped data of all geological bodies, and obtaining forward modeling data belonging to the geological body after matching the density of the corresponding geological body;

[0095] S103, obtain the influence value of the middle layer material by forward calculation of the upright cuboid; perform fine Bouguer gravity anomaly calculation.

[0096] The method for calculating Bouguer gravity anomaly by forward modeling using geological information variable density provided by an embodiment of the present invention includes:

[0097] The terrain correction and intermediate layer correction in the traditional calculation method are improved to the forward calculation method using geological information variable density, that is, according to the density information of various planes and deep parts such as geological maps and boreholes, the plane, vertical layering and block are used to calculate the influence of the gravity measuring points, and the sum is the comprehensive influence of the material between the geoid and the terrain undulation surface on the gravity value of the measuring point. The improved Bouguer gravity anomaly calculation method is normal field correction, height correction and intermediate material forward correction.

[0098] To facilitate understanding of the right rectangular parallelepiped formula, the right rectangular parallelepiped formula is described below.

[0099] The undulation of the actual terrain surface can be fitted with different surfaces. The commonly used methods are square surface integral method, square midpoint elevation method and square average elevation method. After the terrain surface fitting method is determined, it is extended to a certain depth in the vertical direction to form a three-dimensional model body. Related research shows that when the grid spacing is less than 10m×10m, the three-dimensional bodies formed by the above three fitting methods have basically the same gravity influence value for the same measuring point. Considering the characteristic that the terrain grid data is expressed in the form of node three-dimensional coordinates ( Figure 4 ), and in large-scale gravity exploration, large-scale, high-precision terrain data have generally been surveyed and mapped, and the grid spacing is much less than 10m. The right rectangular parallelepiped method has obvious advantages in terms of computing speed when the amount of data is huge and in terms of ease of computer programming. Therefore, the right rectangular parallelepiped formula is used for the intermediate material forward correction calculation.

[0100] like Figure 5 , O is the origin of coordinates, Z axis is vertical downward (downward positive direction), X and Y are perpendicular to each other to form a horizontal plane, the two pairs of side surfaces of the upright cuboid are parallel to the XOZ and YOZ planes respectively, the top and bottom surfaces are parallel to the XOY plane, the residual density is σ, a certain volume element in the geological body dv = dxdydz, its coordinates in the XYZ coordinate system are (x, y, z), its residual mass is dm, then dm = σdv = σdxdydz, the residual mass element to the calculation point P (x p ,y p , z p ) is r=[(xx p ) 2 +(yy p ) 2 +(zz p ) 2 ] 1 / 2 Therefore, the gravitational potential W generated by the residual mass of the geological body on a certain point is:

[0101]

[0102] Where: G is the gravitational constant, V is the volume of the geological body.

[0103] Since the Z direction is the direction of gravity, the gravity anomaly is the derivative of the gravitational potential of the residual mass along the Z direction:

[0104]

[0105] If the coordinate origin is translated to the calculation point P, the coordinates of the measurement point of the integrand are x p ,y p 、z p are all 0, solve the integral:

[0106]

[0107] Restore the coordinate origin and use point C to represent the center point of the cuboid. The coordinates are (x c ,y c , z c ), the length, width and height are 2a, 2b and 2c respectively, then:

[0108] r=[(x c -x p ) 2 +(y c -y p ) 2 +(z c -z p ) 2 ] 1 / 2

[0109] x 1 =x c -ax p ,y 1 =y c -by p , z 1 =z c -cz p

[0110] x 2 =x c +ax p ,y 2 =y c +ay p , z 2 =z c +cz p

[0111] The above integral form can be expressed as:

[0112]

[0113] Where: x∈[x 1 , x2 ],y∈[y 1 ,y 2 ],z∈[z 1 , z 2 ]

[0114] In order to facilitate programming calculations and reduce loops to increase the calculation speed, formula (10) can be further decomposed into 8 calculation formulas. The calculation results of the 8 formulas correspond to Figure 5 The sum of the terrain correction values ​​of the 8 corner points is the forward modeling value of the right cuboid.

[0115] When using DEM data and the above formula for forward calculation, the following aspects should be noted:

[0116] ① From the formula, we can know that the influence of the rectangular parallelepiped at the same height in the upper and lower halves of the space on the gravity measuring point is equal in magnitude and opposite in sign, that is, the forward value above the plane of the gravity measuring point is negative, and the forward value below the plane of the measuring point is positive, which is different from the terrain correction value, which is only positive;

[0117] ② Figure 5 a and b in the equation correspond to the half-grid spacing of the DEM grid points ( Figure 4 ). Then the half height c of the cuboid can be calculated using the top elevation and bottom elevation;

[0118] ③The three-dimensional coordinates of the DEM grid node are the midpoint coordinates of the top surface of the rectangular parallelepiped ( Figure 4 The z value of the two needs to be converted using the half-height cuboid c.

[0119] ④ In formula (10), ln(y+r), ln(x+r) and When x=y=z=0, x=0&z=0 or y=0&z=0, there is a singularity in the function. In actual programming, a small value can be added to the z coordinate to avoid it.

[0120] In order to facilitate understanding of the embodiments of the present application, the implementation process is described in detail as follows:

[0121] From the measured absolute gravity values, DEM grid data and measured density values ​​(plane or vertical, not limited to geological maps), the following process is required to calculate the fine Bouguer gravity anomaly. The technical route is as follows: Figure 6 .

[0122] 1. Correct the gravity value of the measuring point to normal field and height

[0123] The gravity data of the measuring point includes the three-dimensional coordinates of the measuring point x, y, h and G 观The measured gravity values ​​are four columns, and the normal field correction and height correction are consistent with the traditional calculation method, specifically formula (1)(5). The corrected result is four columns of data, namely the three-dimensional coordinates of the measuring point x, y, h and the gFI corrected gravity value.

[0124] 2. Grid data layer preparation

[0125] The grid file is layered at a certain vertical interval. According to the forward calculation of the rate of change of gravity value with depth under fixed grid spacing, the DEM elevation of the whole area can be divided into two large layers, 400m below the lowest point of the whole area, and the next large layer is divided into two layers at 50m intervals from 400-500m, 100m intervals from 500-700m, 200m intervals from 700-1100m, 400m intervals from 1100-1500m, and one layer from 1500m to absolute 0m elevation according to the actual remaining interval; the previous large layer is continuously layered at 20m intervals to the highest point of the whole area, and the next large layer, that is, all grid points, are adopted. After the previous large layer is out of the ground, the grid points whose actual surface grid elevation is lower than the elevation of the divided layer will be discarded. In the layering process, the center point elevation Zm of the grid point layer and the half height c of the grid cube are calculated.

[0126] 3. Preparation of geological interface files

[0127] The actual density value is obtained according to the geological map of the gravity exploration area or other equally divided lithology acquisition measurement statistics of other boreholes, and the statistical density value of each geological body is obtained, that is, the geological body name and statistical density value. The geological interface file can be prepared in the following two steps.

[0128] Simplify the boundaries of the geological body, and count the interface information of the geological body in the area. The interface information includes x1, y1, z1, x2, y2, z2, dip angle A, relationship J between the body and the surface (the body is located above or below the surface), and strike (strike is calculated by the x and y coordinates of the two points) of boundary points 1 and 2 of the geological body, and the addition and subtraction relationship L between the body and the dip. Prepare the above information for each surface and enclose the geological body with multiple surfaces.

[0129] The above is also the process of establishing a three-dimensional geological upright rectangular model. The geological map needs to be appropriately deleted and selected. The areas close to the gravity measurement points can be layered in detail, while areas far away can be layered roughly. After digitizing the coordinates of the geological boundary position, occurrence or drilling location, depth, etc., write a program to complete the preparation of the above geological interface files.

[0130] 4. Forward modeling data preparation

[0131] The grid layered data is segmented by the closure of each geological interface in step 3. After segmentation, the grid data belonging to the geological body is formed, that is, the data of all upright cuboids belonging to the body, and the density σ of the geological body is given.

[0132] After preparation, the rectangular parallelepiped data of each geological body consists of 5 columns, namely X, Y, Zm, C, and σ.

[0133] 5. Obtain the influence value of the middle layer material by forward calculation of the upright cuboid

[0134] The above three steps form the geological body data file that is finally used for the forward calculation of the rectangular block. Combined with the x, y, and h in the gravity value file of the measuring point, the rectangular block forward modeling formula (10) can be used to perform the forward calculation of a single rectangular block in each geological body data file for a certain measuring point. After the calculation is completed, the forward modeling values ​​of all rectangular blocks are added together, which is the influence value of the geological body on the intermediate layer material of a certain gravity measuring point.

[0135] According to the above steps, the calculation of the influence value of the intermediate layer material of all measuring points and all geological bodies is completed in a loop. After completion, the result file contains four columns, namely x, y, h and g of the measuring point. 中 .

[0136] 6. Calculation of fine Bouguer gravity anomaly

[0137] Subtract the gFI calculated in the first step from the gFI calculated in the fifth step. 中 , that is, to obtain the fine Bouguer gravity anomaly value of the gravity measuring point, and the result is four columns, namely x, y, h and gb of the measuring point.

[0138] The invention has been applied in the large-scale gravity exploration of the III mining area of ​​the project "Research on the Effectiveness of the Technical Methods for Prospecting and Exploration in the Deep Side and Periphery of the Jinchuan Copper-Nickel Deposit". In order to better explain the present invention and facilitate understanding, the present invention is described in detail below in conjunction with the accompanying drawings through specific implementation methods.

[0139] The gravity profile point spacing is 50m, the profile length is 1km, and there are 21 measuring points in total. The terrain conditions within 1km around the measuring points (DEM data with a grid spacing of 10m) are as follows: Figure 7 、Geological map Figure 8 , the geological profile is as follows Fig. 9 , prepare 21 measuring points of x, y, h and G 观 , some data of gravity measurement points are shown in Table 1, DEM grid data are shown in Table 2, and geological information is shown in the geological map ( Figure 8 ) and geological profiles ( Fig. 9 );

[0140] Table 1 Observation data of gravity measurement points and gFI gravity correction values

[0141] Serial number dot number x(degrees) y(degree) h(m) <![CDATA[G 观 (mGal)]]> gFI(mGal) 1 166 38.48626155 102.1309904 1708.006 979967.294 458.623 2 167 38.48669978 102.1312811 1690.140 979970.880 456.661 …… …… …… …… …… …… …… 21 186 38.49335622 102.1381600 1599.376 979986.343 443.548

[0142] Table 2 DEM grid data file format

[0143]

[0144] Step 1: Use formulas (1) and (5) to perform normal field correction and altitude correction. The corrected gFI gravity correction value is shown in the rightmost column of Table 1.

[0145] Step 2: As shown in Table 2, the DEM grid spacing is calculated using the first two columns of coordinates of adjacent grid points. Taking a rectangular parallelepiped with a side length of 10m, the gravity influence value calculated by the experiment of varying depth and layer spacing is less than 0.001×10 -5 m / s 2 (Same as mGal), the maximum elevation of the DEM grid data elevation (the third column) is 1879.658, and the minimum elevation is 1391.712. Round to zero with multiples of 20, then the maximum elevation layer of the grid point surface is 1860, and the minimum elevation layer is 1380. First, use the surface 400m below the minimum elevation, that is, the surface at an altitude of 980 to divide the DEM data into two large layers, because there are grid points with a distance (vertical depth) less than 400m in the upper large layer. Therefore, the altitude 1880-980 is divided into 45 layers of 1860, 1840, 1820...980 with a stratification interval of 20. Select the point whose DEM grid point elevation is greater than the stratification elevation, divide the sum of the DEM actual elevation and the stratification elevation by 2, and calculate the midpoint elevation of the composed rectangular solid. Divide the difference between the DEM actual elevation and the stratification elevation by 2 to calculate the half-height of the composed rectangular solid.

[0146] The layering principle below 980 is the same as above. 980-880 is divided into 930 and 880 layers at 50m intervals, 880-680 is divided into 780 and 680 layers at 100m intervals, 680-280 is divided into 480 and 280 layers at 200m intervals, and 280-0m is the last layer. The calculation method of the midpoint elevation of the upright cuboid and the half-height of the upright cuboid is consistent with the calculation method used in the previous large layer.

[0147] After the stratification is completed, each layer has four columns: X, Y, Zm and c. All stratified data are unified into one file. After completion, the distribution of grids in the 1480-1500m elevation layer is as follows: Fig.10 .

[0148] Step 3: Extract the endpoint coordinates of each geological boundary on the geological map, the occurrence on the geological profile, etc., and establish each geological interface. A total of 9 geological body interface files are established this time, body 1 is Sc-qbs geological body, body 2 is Gn-bm geological body, body 3 is Mi geological body, body 4 is ΣⅢ geological body, body 5 is Q geological body (above 1460 above sea level), body 6 is Gn-ba geological body (1460-0m above sea level), body 7 is ML geological body, body 8 is ΣI geological body, and body 9 is Gn-ba geological body; the geological interface file of ΣⅢ geological body is shown in Table 3. The geological interface files of 9 geological bodies are used to trap and segment the DEM grid layered data, match the measured statistical density values ​​of 9 geological bodies, and form the final forward calculation data files of 9 geological bodies. The forward modeling data file of ΣⅢ geological body (i.e., the upright rectangular parallelepiped file) is shown in Table 3. Fig.11 .

[0149] Table 3 Geological interface files of ΣⅢ geological bodies

[0150]

[0151] Step 4: Using the three-dimensional coordinates x, y, h of the gravity measuring point and the rectangular parallelepiped data of each geological body trap prepared in the third step, the forward modeling values ​​of all rectangular parallelepipeds within 1000m of the plane distance of each gravity measuring point are respectively forward modeled according to formula (10). The sum is the influence value of the intermediate layer material of all geological bodies on the measuring point. The calculation of the influence value of the intermediate layer material of all measuring points is completed in a cycle, which is g 中 .

[0152] Step 5: Subtract the gFI obtained in the first step from the g calculated above. 中 , is the Bouguer gravity anomaly value obtained by careful calculation.

[0153] It should be noted that the embodiments of the present invention can be implemented by hardware, software, or a combination of software and hardware. The hardware portion can be implemented using dedicated logic; the software portion can be stored in a memory and executed by an appropriate instruction execution system, such as a microprocessor or dedicated design hardware. A person of ordinary skill in the art will appreciate that the above-mentioned devices and methods can be implemented using computer executable instructions and / or contained in a processor control code, such as a carrier medium such as a disk, CD or DVD-ROM, a programmable memory such as a read-only memory (firmware), or a data carrier such as an optical or electronic signal carrier. Such code is provided on the carrier medium. The device and its modules of the present invention can be implemented by hardware circuits such as very large-scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, etc., or programmable hardware devices such as field programmable gate arrays, programmable logic devices, etc., can also be implemented by software executed by various types of processors, and can also be implemented by a combination of the above-mentioned hardware circuits and software, such as firmware.

[0154] The present invention has achieved good results in large-scale gravity exploration of the Jinchuan deposit III mining area, as described below:

[0155] In order to demonstrate the exploration effect, the Bouguer gravity anomaly calculated in detail is used for trend analysis. After obtaining the residual gravity anomaly, the exploration effects of the traditional method and the improved method are compared and analyzed. Fig.12 .

[0156] like Figure 8 The ultrabasic rock mass on the surface is roughly located between point numbers 173 and 181, with a dip angle of 60° on the south side and an attitude of 75° on the north side. Fig.12 , the red line is the residual gravity anomaly calculated by the improved method. Compared with the traditional calculation method of the blue line, the anomaly starts at 170 points, but the anomaly calculated by the traditional method ends at 176 points, while the anomaly calculated by the improved method extends to 180 points. The anomaly range and shape of the improved calculation are more consistent with the actual geological conditions; the traditional method calculates a false anomaly display on the north side, which is also inconsistent with the actual geological conditions. As can be seen from the above figure, the exploration effect of the improved method has been significantly improved.

[0157] The above description is only a specific implementation mode of the present invention, but the protection scope of the present invention is not limited thereto. Any modification, equivalent substitution and improvement made by any technician familiar with the technical field within the technical scope disclosed by the present invention and within the spirit and principle of the present invention should be covered by the protection scope of the present invention.

Claims

1. A method for calculating Bouguer gravity anomaly by forward modeling using geological information variable density, characterized in that: include: Using geological information of the explored area, including the boundary, occurrence, depth of geological bodies in geological maps or the density change of geological bodies vertically in the borehole, according to the three-dimensional distribution of geological bodies, not only block in the plane, but also layer in the vertical direction, and assign the measured density values ​​to different geological bodies to calculate the variable density; Using high-resolution DEM data to vertically layer different density bodies, variable density forward modeling is used to calculate the influence value of the intermediate layer material instead of terrain correction and intermediate layer correction step-by-step calculation of the unified density geological body, thus improving the accuracy of Bouguer gravity anomaly; Using three-dimensional geological body variable density calculation, the deep stratum occurrence is adopted as the arbitrary occurrence angle of field geological measurement or drilling control; The method specifically comprises: Step 1, obtain the geological information of the surveyed area, the three-dimensional coordinates of the gravity measuring points, the measured gravity values ​​of the gravity measuring points and the DEM grid data (X, Y, Z1); use the three-dimensional coordinates of the gravity measuring points and the measured gravity values ​​to complete the normal field and height correction formula and obtain the corrected gravity value gFI; Step 2, calculate the grid spacing of the DEM grid data, stratify the DEM grid data according to the three-dimensional distribution state of the geological body, obtain the rectangular parallelepiped data file after stratification, and calculate the center point elevation Zm and half height C of the rectangular parallelepiped at the same time; Step 3, using geological information to establish a geological interface file, using the established multiple geological interfaces to segment the above rectangular parallelepiped data file to form a rectangular parallelepiped file of a geological body, and matching the measured density of the geological body, and completing the segmentation of the rectangular parallelepiped files and the measured density matching of all geological bodies in a cycle, that is, completing the forward data preparation; Step 4: Using the three-dimensional coordinates of the gravity measurement points, grid spacing and forward modeling data, the influence value of the intermediate layer material is calculated according to the right cuboid forward modeling formula to obtain the intermediate layer material correction value g 中 ; Step 5: Compare the calculated corrected gravity value gFI with the obtained g 中 Subtract them to get the fine Bouguer gravity anomaly value of the gravity measuring point.

2. The method for calculating Bouguer gravity anomaly by forward modeling using geological information variable density as claimed in claim 1, characterized in that: The method comprises the following steps: The method of obtaining the gFI gravity correction value by correcting the three-dimensional coordinates of the gravity measuring point and the measured gravity value in step 1 through normal field correction and height correction; the gravity data of the measuring point includes the three-dimensional coordinates of the measuring point x, y, h and G 观 Four columns of measured gravity values ​​are used to perform normal field correction and height correction to obtain the three-dimensional coordinates of the measuring point x, y, h and the gFI corrected gravity value.

3. The method for calculating Bouguer gravity anomaly by forward modeling using geological information variable density as claimed in claim 1, characterized in that: The grid data stratification scheme in step 2 needs to be determined according to the grid spacing of the DEM data. The stratification depth and stratification interval are determined by the influence value change rate calculation test based on the grid spacing. The larger the grid spacing, the smaller the stratification depth and the finer the stratification interval. When stratifying, the grid points in the previous large layer whose elevation is lower than the layer height will be discarded, while all the grid points in the next large layer will be retained. After completion, a unique undivided upright rectangular data file is formed, including four columns of data: grid point coordinates X and Y, the elevation Zm of the midpoint of the rectangular block, and the half-height C of the rectangular block. The grid point coordinates are repeated many times because each layer corresponds to these grid points.

4. The method for calculating Bouguer gravity anomaly by forward modeling using geological information variable density as claimed in claim 1, characterized in that: In step 3, the geological interface file is prepared, the vertical rectangular parallelepiped file belonging to the closed geological body is segmented according to the geological interface and a density value is assigned to form a forward modeling data file.

5. The method for calculating Bouguer gravity anomaly by forward modeling using geological information variable density as claimed in claim 1, characterized in that: In step 4, the intermediate layer influence value is forward modeled according to the variable density of the right cuboid formula, and the forward modeling is used to replace the intermediate layer correction and terrain correction methods in the traditional Bouguer gravity anomaly calculation.

6. The method for calculating Bouguer gravity anomaly by forward modeling using geological information variable density as claimed in claim 1, characterized in that: In step 5, the normal field and height-corrected gFI is subtracted from the intermediate layer material effect value g 中 , and obtain fine Bouguer gravity anomalies.

7. The method for calculating Bouguer gravity anomaly by forward modeling using geological information variable density as claimed in claim 1, characterized in that: The stratification of grid data includes: The grid spacing is equal to the distance between two adjacent grid points; The grid files are layered at a certain vertical interval. Based on the forward calculation of the rate of change of gravity value with depth under a fixed grid spacing, the DEM elevation of the entire area is divided into two large layers, 400m below the lowest point of the entire area: The upper large layer is continuously divided into layers at intervals of 20m to the highest point of the entire area, and the grid points whose actual surface grid elevation is lower than the elevation of the divided layer are discarded; The next large layer is further divided into two layers at 50m intervals according to the lowest elevation of 400-500m, 500-700m is divided into two layers at 100m intervals, 700-1100m is divided into two layers at 200m intervals, 1100-1500m is divided into one layer at 400m intervals, 1500-2300m is divided into one layer at 800m intervals, if the lowest point elevation is greater than 2700m, then 2300m to 0m elevation is divided into one layer according to the actual remaining interval, if the lowest point elevation is less than 2700m, then 1500m to 0m elevation is divided into one layer according to the actual remaining interval, and all grid points in this layer are retained; After stratification, the elevation of the grid center point and the half-height of the upright cuboid are calculated using the stratified elevation values. The elevation of each center point is the sum of the elevation values ​​of the upper and lower layers divided by 2, and the half-height of each grid cube is the difference between the elevation values ​​of the upper and lower layers divided by 2.

8. The method for calculating Bouguer gravity anomaly by forward modeling using geological information variable density as claimed in claim 1, characterized in that: Preparation of geological interface files using measured density values ​​includes: (1) According to the geological map of the gravity exploration area or other boreholes, the actual density values ​​are measured and collected by lithology to obtain the statistical density value of each geological body; the statistical density value of each geological body includes: the name of the geological body and the statistical density value; (2) Prepare geological interface files: Simplify the boundary of the geological body, count the interface information of the geological body in the area, enclose the geological body with multiple faces, and segment the grid data layered file completed in step 1 with the enclosed geological body. The segmented file matches the measured density value of the enclosed geological body, thus completing the density body preparation; The interface information includes: digitization of geological boundary position coordinates, occurrence or drilling position, and depth to form multiple interface information of a geological body; the interface information includes: x1, y1, z1, x2, y2, z2, inclination A, relationship J between body and surface, and addition and subtraction relationship between strike and dip of geological boundary point 1 and point 2.

9. A device for forward calculating Bouguer gravity anomaly using geological information variable density according to any one of claims 1 to 8, characterized in that: The device comprises: The acquisition module acquires the geological information of the explored area, the three-dimensional coordinates of the gravity measurement points, the measured gravity values ​​of the gravity measurement points and the DEM grid data; The first determination module is used to perform normal field correction and height correction according to the three-dimensional coordinates of the gravity measurement point and the measured gravity value to obtain the gFI corrected gravity value; The second determination module is used to vertically layer the DEM grid data according to the principle of the rate of change of a certain right cuboid with the distance and size from the measuring point to form a grid layer file; The third determination module is used to organize the simplified geological information into a geological interface file, and according to the interface file, the above-mentioned grid layer file is segmented and matched with the measured density of the geological body, and finally a right cuboid forward modeling data file is formed; The fourth determination module is used for forward calculation of the impact value of the intermediate layer material; The fifth determination module is used to calculate the refined Bouguer gravity anomaly value.

Images