Ore body rapid positioning method based on time domain spectrum induced polarization middle ladder

Abstract

The invention provides an ore body rapid positioning method based on a time domain spectrum induced polarization middle ladder. The method comprises the following steps: generating down-extended induced polarization middle ladder apparent polarization intensity signal files of a plurality of underground layers based on a spatial domain stable continuation algorithm; three-dimensional normalized total gradient values of different horizons below a measuring point of a measuring line where anomalies are located are calculated based on abnormal distribution in a contour map of lower-extending induced polarization intermediate ladder apparent polarization intensity signal files of multiple underground horizons and de-noised induced polarization intermediate ladder parameter signal files of the lower-extending induced polarization intermediate ladder apparent polarization intensity signal files of the multiple underground horizons; determining the central depth position of the ore body based on the three-dimensional normalized total gradient value; and determining the boundary of the ore body by using a normalized singular spectrum entropy method. According to the method, the center depth position of the ore body is accurately determined based on the three-dimensional normalized total gradient value, the understanding and recognition capability on the internal structure of the ore body is improved, meanwhile, the normalized singular spectrum entropy method is introduced to determine the boundary of the ore body, the form and range of the ore body can be further revealed, and a reference basis is provided for subsequent exploration and mining.

Description

Technical Field

[0001] The invention relates to the technical field of resource exploration, and in particular to a method for quickly locating an ore body based on a time domain spectrum induced polarization ladder. Background Art

[0002] In recent years, with the continuous advancement of geophysical exploration technology, especially in the field of electromagnetic method, induced polarization mid-elevator technology has attracted widespread attention due to its high resolution and good ability to distinguish different geological bodies. By measuring the induced polarization effect and resistivity change generated by the underground medium under the action of the electric field, the induced polarization mid-elevator technology can effectively reveal the existence of underground ore bodies and their spatial distribution characteristics. However, when faced with a multi-layer ore body structure, a single detection method is often difficult to accurately distinguish the ore body information at different levels, especially when the resistivity difference between ore bodies is small or the geological background is complex. How to achieve fast and accurate layered positioning has become a key issue that needs to be solved urgently. Summary of the invention

[0003] To solve the above problems, the purpose of the present invention is to provide a method for rapid ore body positioning based on time domain spectral induced polarization.

[0004] To achieve the above object, the present invention provides the following solutions:

[0005] A method for rapid ore body positioning based on time domain spectrum induced polarization ladder, comprising:

[0006] Step 1: Construct an exploration coordinate system with any point in the study area as the origin, and use the induced polarization ladder method to calculate the apparent resistivity and apparent polarizability values ​​of each measuring point in the study area;

[0007] Step 2: Use the Cole-Cole model to construct a time-domain spectral induced polarization forward model, and input the apparent polarizability value of each measuring point into the time-domain spectral induced polarization forward model to generate the polarization intensity parameters of each measuring point;

[0008] Step 3: Rotate and grid the exploration coordinate system based on the azimuth of the survey line, and draw contour maps of induced polarization gradient parameters and topographic contour maps on the gridded XOY plane; the induced polarization gradient parameters include: apparent resistivity, apparent polarizability value and polarization intensity parameter;

[0009] Step 4: using the energy spectrum density area fractal method to perform signal-noise separation on the contour map of the induced polarization ladder parameter to obtain a denoised induced polarization ladder parameter signal file, and determining the abnormal distribution in the contour map of the denoised induced polarization ladder parameter signal file;

[0010] Step 5: Generate the downward extended induced polarization intensity signal files of multiple layers of underground induced polarization based on the denoised induced polarization parameter signal files based on the spatial domain stable continuation algorithm;

[0011] Step 6: Calculate the three-dimensional normalized total gradient values ​​of different layers below the measuring point of the survey line where the anomaly is located based on the abnormal distribution in the contour map of the downward extended induced polarization gradient polarization intensity signal file of multiple underground layers and the denoised induced polarization gradient parameter signal file, and determine the central depth position of the ore body based on the three-dimensional normalized total gradient value;

[0012] Step 7: Determine the boundary of the ore body using the normalized singular spectral entropy method.

[0013] Preferably, in step 1, the direction perpendicular to the structural strike is used as the direction of the transmitting line and the receiving line, and the survey line is laid, and an exploration coordinate system is constructed with the direction along the survey line as the x direction, the direction perpendicular to the survey line as the y direction, and the vertical downward direction as the z direction.

[0014] Preferably, in step 1, the apparent resistivity calculation formula is:

[0015]

[0016] ρ a is the apparent resistivity of the measuring point, K is the device coefficient, x0 is the abscissa position of the midpoint of MN, y0 is the ordinate position of the midpoint of MN, AB is the length of the transmitting electrode, MN is the length of the receiving electrode, and U0 is the primary field voltage;

[0017] The calculation formula of apparent polarizability is:

[0018]

[0019] Among them, η a is the apparent polarization rate of the measuring point, and ΔU2 is the secondary field voltage after power failure.

[0020] Preferably, in step 2, based on a digital filtering algorithm, 21 filter coefficients are used to obtain a time domain spectral induced polarization forward model in an accumulative manner; wherein the time domain spectral induced polarization forward model is:

[0021]

[0022] Among them, ω r is the frequency variable of the real component of r from 1 to 21 corresponding to time t, a r is the horizontal coordinate corresponding to the 21 filter coefficients, φ r is with a r The corresponding filter coefficient, τ is the time constant, c is the frequency correlation coefficient, ρ0 is the zero-frequency resistivity, R r is the real component, I r Is the imaginary quantity.

[0023] Preferably, in step 3, the rotation of the exploration coordinate system is completed based on the azimuth of the survey line; wherein, the coordinate transformation relation formula is:

[0024] x' = xcosθ - ysinθ

[0025] y' = ycosθ + xsinθ

[0026] x = (x' - y'tanθ)sinθ

[0027]

[0028] wherein, x is the coordinate of the original exploration coordinate system on the x-axis, y is the coordinate of the original exploration coordinate system on the y-axis, x' is the coordinate of the exploration coordinate system on the x-axis, and y' is the coordinate of the rotated exploration coordinate system on the y-axis.

[0029] Preferably, in step 6: calculating the three-dimensional normalized total gradient values of different horizons below the measuring points of the survey line where the anomaly is located based on the anomaly distribution in the contour map of the downward continuation IP mid-gradient apparent polarization intensity signal file and the denoised IP mid-gradient parameter signal file of multiple horizons underground, and determining the central depth position of the ore body based on the three-dimensional normalized total gradient values, including:

[0030] Step 6.1: According to the anomaly distribution range in the contour map of the denoised IP mid-gradient parameter signal file, extract the survey line passing through the center of the anomaly on the Y-axis. On the XOZ plane where the survey line is located, interpolate the three-dimensional normalized total gradient values of different horizons in the Z direction to obtain the three-dimensional normalized total gradient plane file on the XOZ two-dimensional plane where the center of the anomaly is located;

[0031] Step 6.2: Take the position where the maximum value of the three-dimensional normalized total gradient plane file on the Z-axis is located as the central depth position of the ore body.

[0032] Preferably, in step 7: using the normalized singular spectrum entropy method to determine the boundary of the ore body, including:

[0033] Step 7.1: Obtain the number num of actual observed measuring points of the same survey line in the X-axis direction at the central depth position of the ore body;

[0034] Step 7.2: When 3 < num < 5, directly take the boundary of the anomaly field as the boundary of the ore body. When num ≥ 5, use the normalized singular spectrum entropy method to infer the boundary of the ore body in the X-axis direction to obtain the distribution range of the ore body in the X-axis direction;

[0035] Step 7.3: Determine the boundary of the ore body according to the distribution range of the ore body in the X-axis direction.

[0036] Preferably, after step 7, it further includes:

[0037] Step 8: Using a sphere as the forward model, taking the central depth position of the ore body and the distribution range of the ore body in the X-axis direction as spatial position constraints, perform equivalent source fitting on the apparent resistivity and apparent polarizability values ​​of each measuring point to obtain the spatial distribution of the ore body.

[0038] Preferably, in step 8, the numerical calculation formula for the apparent resistivity of the (i, j)th equivalent source at the point (k, l) is:

[0039]

[0040] The numerical calculation formula for the apparent polarizability of the (i, j)th equivalent source at point (k, l) is:

[0041]

[0042] x=x ij -x kl

[0043] y=y ij -y kl

[0044] h0=hz kl

[0045]

[0046] in, represents the apparent resistivity value of the (i,j)th equivalent source at point (k,l), represents the apparent polarizability value of the (i, j)th equivalent source at the point (k, l), h represents the central burial depth of the three-dimensional polarization sphere, r0 represents the radius of the three-dimensional polarization sphere, is the relative resistivity, ρ2 and ρ c are the resistivity of the equivalent sphere and the initial resistivity in the XOZ plane, (x ij ,y ij ) represents the horizontal coordinate of the equivalent source position, (x kl ,y kl , z kl ) represents the coordinates of the observation point, η2 represents the equivalent source sphere polarizability, M ij Represents an intermediate variable.

[0047] The present invention also provides a computer-readable storage medium having a computer program stored thereon, characterized in that when the computer program is executed by a processor, the steps of the above-mentioned method for rapid locating an ore body based on a ladder in time-domain spectral induced polarization are implemented.

[0048] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0049] The present invention provides a method for quickly locating an ore body based on the gradient of induced polarization in time domain spectrum. Compared with the prior art, the present invention accurately determines the central depth position of the ore body based on the abnormal distribution and three-dimensional normalized total gradient value in the gradient parameter signal of induced polarization, improves the understanding and recognition ability of the internal structure of the ore body, and introduces the normalized singular spectrum entropy method to determine the boundary of the ore body, which can further reveal the shape and scope of the ore body and provide a reference basis for subsequent exploration and mining. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0051] Figure 1 A flow chart of a method for rapid ore body positioning based on time domain spectrum induced polarization provided by the present invention;

[0052] Figure 2 It is a statistical diagram of physical property parameters in the embodiment of the present invention.

[0053] Figure 3 Schematic diagram of the spatial distribution of ABMN electrodes in the intermediate gradient device.

[0054] Figure 4 This is the Cole-Cole model fitting the measured induced polarization curve.

[0055] Figure 5 It is the original contour diagram of the ladder parameters of apparent resistivity, apparent polarizability, and apparent polarization intensity in induced polarization.

[0056] Figure 6 It is the double logarithmic relationship diagram of energy spectrum density and area under the SA method.

[0057] Figure 7 It is the contour map of the noise part of the polarization intensity m after SA decomposition.

[0058] Figure 8 It is the contour diagram of the ladder parameters of apparent resistivity, apparent polarizability, and apparent polarization intensity signals in IP.

[0059] Fig. 9 It is a downward extension of the input elevation data, polarization intensity m signal data and the output downward extension contour map of 100, 150, 200 and 300 meters depth.

[0060] Fig.10It is the downward extension contour map of the polarization intensity m signal at altitudes of 1550, 1500, 1450, 1400, 1350, 1300, 1250, 1200, 1100, and 1000 meters.

[0061] Fig.11 It is the three-dimensional normalized total gradient contour map of polarization intensity m at depths of 100, 150, 200, and 300 meters.

[0062] Fig.12 It is the three-dimensional normalized total gradient contour map of the polarization intensity m signal at altitudes of 1550, 1500, 1450, 1400, 1350, 1300, 1250, 1200, 1100, and 1000 meters.

[0063] Fig.13 are contour plots of apparent polarizability and normalized singular spectral entropy for the three cubes. DETAILED DESCRIPTION

[0064] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0065] In order to make the technical problems, technical solutions and beneficial effects to be solved by the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0066] To achieve the above object, the technical solution adopted by the present invention is:

[0067] See also Figure 1 , a method for rapid ore body positioning based on time domain spectrum induced polarization ladder, comprising:

[0068] Step 1: Construct an exploration coordinate system with any point in the study area as the origin, and use the induced polarization ladder method to calculate the apparent resistivity and apparent polarizability values ​​of each measuring point in the study area;

[0069] In step 1, the present invention uses the direction perpendicular to the structural strike as the direction of the transmitting line and the receiving line, and lays the survey line. At the same time, the direction along the survey line is the x direction, the direction perpendicular to the survey line is the y direction, and the vertical downward direction is the z direction to construct an exploration coordinate system, and step 1 includes the following steps:

[0070] S11. Lay out electrodes and transmitting lines according to the specifications, use electrodes as measuring points, and record the elevation data of the measuring points;

[0071] The transmitting line and transmitting electrodes are A and B, and the receiving electrodes M and N are arranged in the middle of the transmitting electrode. Within the range of length, the section of the receiving electrode MN on the straight line where the transmitting electrode AB is located is called the main section. It can also be observed on the side parallel to the AB line. The distance from the side section to the main section is not greater than

[0072] S12. A high-power transmitter transmits signals according to a time window to ensure signal strength, suppress random interference, and record the maximum emission current I that can be kept constant. Because it involves induced polarization information, during observation, the number of time windows should be as large as possible and the width of the time window should be as narrow as possible.

[0073] S13. During the acquisition, the parameters recorded by the receiving end are the spatial position of the receiving MN electrode of the measuring point, the primary field voltage U0, and the secondary field voltage ΔU2 after power failure; the transmitting end records the spatial position of AB and the transmitting current I.

[0074] S14. Calculate the apparent resistivity ρ of the measuring point at the coordinate position (x0, y0) a and apparent polarizability η a , the formula is as follows:

[0075] The calculation formula of apparent resistivity of the measuring point is:

[0076]

[0077] Among them, ρ a is the apparent resistivity; K is the device coefficient, and the calculation formula is as follows;

[0078]

[0079] Where x0 is the horizontal coordinate position of the midpoint of MN, and y0 is the vertical coordinate position of the midpoint of MN. The origin of the coordinate is taken at the midpoint of AB. AB is the length of the transmitting electrode, and MN is the length of the receiving electrode. The specific calculation formula is as follows.

[0080]

[0081]

[0082] The calculation formula of the apparent polarizability of the measuring point is:

[0083]

[0084] Among them, η a is the apparent polarization rate of the measuring point, the primary field voltage U0, and the secondary field voltage ΔU2 after power failure.

[0085] During the actual investigation of the apparent polarizability, the instrument mostly records the apparent charging rate obtained in each time window. The apparent charging rate means the area enclosed by the secondary field voltage ΔU2 decay curve and the time axis in the time window. After the apparent charging rate is normalized by the time interval of the time window, the apparent polarizability can be obtained.

[0086] Step 2: Use the Cole-Cole model to construct a time-domain spectral induced polarization forward model, and input the apparent polarizability value of each measuring point into the time-domain spectral induced polarization forward model to generate the polarization intensity parameters of each measuring point;

[0087] Furthermore, step 2 includes:

[0088] S21: Based on Laplace transform, the frequency domain Cole-Cole model expression is converted into the time domain, and the time domain spectral induced polarization forward model is constructed;

[0089] The frequency domain formula of the Cole-Cole model in the frequency domain is:

[0090]

[0091] Where ρ0 is the zero-frequency resistivity, j is the imaginary unit, ω is the angular frequency, m is the polarizability, τ is the time constant, and c is the frequency-dependent coefficient;

[0092] The discrete form of the digital filtering algorithm proposed by Guptasarma is adopted. Based on the digital filtering algorithm, 21 filter coefficients are used to obtain the time-varying resistivity in a cumulative manner, and the time-domain spectral induced polarization forward model is obtained:

[0093]

[0094] ω r is the frequency variable of the real component of r from 1 to 21 corresponding to time t, a r is the horizontal coordinate corresponding to the 21 filter coefficients, φ r is with a r The corresponding filter coefficients. The real components of the two intermediate variables R r 、Immortal Component I r The calculation formula is as follows:

[0095]

[0096] a r is the horizontal coordinate corresponding to the 21 filter coefficients, φ r is with a r The corresponding filter coefficients are shown in Table 1.

[0097] Table 1

[0098] r Filter abscissa Filter coefficient r Filter abscissa Filter coefficient 1 -3.82704 0.000349998 12 -0.956480 0.042972904 2 -3.56608 -0.000418371 13 -0.695520 0.075423603 3 -3.30512 0.000772828 14 -0.434560 0.139346367 4 -3.04416 -0.000171356 15 -0.173600 0.234486236 5 -2.78320 0.001022172 16 -0.956480 0.366178323 6 -2.52224 0.000897638 17 0.34832 0.284615486 7 -2.26128 0.002208974 18 0.60928 -0.235691746 8 -2.00032 0.003844944 19 0.87024 0.046994188 9 -1.73939 0.006809040 20 1.13120 -0.0005901946 10 -1.47840 0.013029162 21 1.39216 0.000570165 11 -1.21744 0.022661391

[0099] S22: Using the apparent polarizability values ​​of the discharge curves of multiple small time windows observed in the time domain as input, the global optimization algorithm of the 1stOpt software is used to obtain the three polarizability parameters of the frequency domain induced polarization method, and to mine the information reflecting the nature of polarization. In this process, it is ensured that the observation data of each measuring point conforms to the Cole-Cole theoretical model, and the noise data that does not conform to the geophysical laws are removed point by point.

[0100] Notice According to the formula in S21, the polarization intensity m, the time constant and the frequency correlation coefficient c can be calculated using the most optimized algorithm based on the time t of each time window and the apparent polarization rate.

[0101] Step 3: Rotate and grid the exploration coordinate system based on the azimuth of the survey line, and draw contour maps of induced polarization gradient parameters and topographic contour maps on the gridded XOY plane; the induced polarization gradient parameters include: apparent resistivity, apparent polarizability value and polarization intensity parameter;

[0102] The specific process of step 3 is as follows:

[0103] S31: According to the coordinate values ​​of the four measuring points A, B, M, and N measured in practice, combined with the azimuth angle θ of the survey line, based on trigonometric calculations, the exploration coordinate system is rotated to the east-west direction of the survey line x, and the north-south direction of the perpendicular survey line y. The rotated coordinate system is selected for calculation during gridding and data processing, while the original coordinate system is still used for interpretation. The specific conversion relationship between the two coordinate systems is as follows;

[0104] x'=xcosθ-ysinθ

[0105] y'=ycosθ+xsinθ

[0106] x=(x'-y'tanθ)sinθ

[0107]

[0108] Among them, x and y are the exploration geodetic coordinates when the original actual data was collected, x′ and y′ are the coordinates after rotation, and the counterclockwise rotation angle is θ.

[0109] S32: Using the observation data in the rotated coordinate system as the gridding input file, grid the induced polarization gradient resistivity ρ, polarization intensity m, and elevation data of the survey area in the rotated exploration coordinate system, and save the gridding file in the form of a GRD file in the Surfer6Text Grid format. The induced polarization gradient apparent resistivity and polarization intensity grid files generated in this step are collectively referred to as the induced polarization gradient parameter original files.

[0110] The contour diagrams of apparent resistivity and polarization intensity of the IP gradient on the XOY plane are collectively referred to as the original contour diagrams of the IP gradient parameters.

[0111] S33: Draw a topographic contour map of the survey area.

[0112] Step 4: using the energy spectrum density area fractal method to perform signal-noise separation on the contour map of the induced polarization ladder parameter to obtain a denoised induced polarization ladder parameter signal file, and determining the abnormal distribution in the contour map of the denoised induced polarization ladder parameter signal file;

[0113] Furthermore, step 4 includes:

[0114] S41: The energy spectrum density area (SA) fractal filter is used for signal-noise separation. The energy spectrum density area (SA) fractal filter is based on the concept of generalized self-similarity, which converts the data from the distribution in the spatial domain to the frequency domain to analyze its fractal characteristics. The relationship between the spectrum energy density and area of ​​the energy spectrum density area (SA) fractal method is:

[0115] A(≥S)∝S -β

[0116] In the formula, S is a certain threshold of energy spectrum density; A (≥ S) represents the area where the energy spectrum density is greater than the threshold S; β is the anisotropic multidimensional fractal index. The relationship between S and A is that as the S value increases, the corresponding A decreases, and its change law mainly depends on β, which is the anisotropic multidimensional fractal index. A double logarithmic image of the energy spectrum density and the number and area of ​​energy spectrum units with energy spectrum density greater than the threshold S is established, and multi-segment linear fitting is performed based on the least squares method, and the intersection point (threshold) is obtained to form different filters. After the inverse Fourier transform, the field separation is realized, and the decomposition result of noise and signal is obtained.

[0117] The original file of the IP ladder parameters is generated through this step to obtain the IP ladder parameter signal file and the IP ladder parameter noise file.

[0118] S42: for the ladder parameter signal file in the induced polarization, divide and identify the abnormal distribution on the XOY plane. The division method can be selected from statistical methods.

[0119] Step 5: Generate the downward extended induced polarization intensity signal files of multiple layers of underground induced polarization based on the denoised induced polarization parameter signal files based on the spatial domain stable continuation algorithm;

[0120] S51: Analyze the terrain file in the study area, obtain the maximum elevation Hmax and the minimum elevation Hmin in the study area, and calculate the height difference H=Hmax-Hmin.

[0121] S52: To obtain the apparent polarization intensity of the depth layer below Hmin, extend from the terrain surface to the plane at the corresponding altitude, use the retained induced polarization gradient parameter signal file and terrain file as input files, and based on the stable extension algorithm in the spatial domain, obtain the downward extended induced polarization gradient apparent polarization intensity signal file at the specified altitude in the Z-axis direction.

[0122] S53: To obtain the apparent polarization intensity within the range of H depth below the terrain surface, extend from the terrain surface to the depth surface parallel to the terrain surface at the corresponding depth, use the retained induced polarization ladder parameter signal file as the input file, and obtain the downward extended induced polarization ladder apparent polarization intensity signal file from the terrain surface to the underground H depth based on the spatial domain stable extension algorithm. It should be noted that step 5 of the present invention can be completed by patent document 201910708065.4, a gravity magnetic potential field analytical extension method based on the spatial domain.

[0123] Step 6: Calculate the three-dimensional normalized total gradient values ​​of different layers below the measuring point of the survey line where the anomaly is located based on the abnormal distribution in the contour map of the downward extended induced polarization gradient polarization intensity signal file of multiple underground layers and the denoised induced polarization gradient parameter signal file, and determine the central depth position of the ore body based on the three-dimensional normalized total gradient value;

[0124] Furthermore, step 6 includes:

[0125] Step 6.1: According to the anomaly distribution range in the contour map of the gradient parameter signal file of the denoised IP, extract the survey line passing through the center of the anomaly on the Y axis, and on the XOZ plane where the survey line is located, perform numerical interpolation of the three-dimensional normalized total gradients of different layers in the Z direction to obtain the three-dimensional normalized total gradient plane file on the XOZ two-dimensional plane where the center of the anomaly is located;

[0126] Step 6.2: The position of the maximum value of the three-dimensional normalized total gradient plane file on the Z axis is taken as the central depth position of the ore body; it should be noted that the maximum value within the depth interval H below the depth layer Hmin in the Z axis direction and the terrain surface is obtained. The position of the maximum value on the Z axis is the center position of the anomaly that causes the abnormal field, which can be considered as the central depth position of the ore body.

[0127] Step 7: Determine the boundary of the ore body using the normalized singular spectral entropy method.

[0128] Preferably, the step 7: determining the boundary of the ore body by using the normalized singular spectral entropy method comprises:

[0129] Step 7.1: Obtain the number num of actual observation points on the same survey line at the central depth position of the ore body in the X-axis direction;

[0130] Step 7.2: When 3 < num < 5, directly use the boundary of the abnormal field as the boundary of the ore body. When num ≥ 5, use the normalized singular spectrum entropy method to infer the boundary of the ore body in the X-axis direction, and obtain the distribution range of the ore body in the X-axis direction;

[0131] Step 7.3: Determine the boundary of the ore body according to the distribution range of the ore body in the X-axis direction.

[0132] In practical applications, the normalized singular spectrum entropy method proposed by the present invention enhances the numerical characteristics of the abnormal body boundary by utilizing the characteristics that the abnormal field information changes violently at the boundary of the abnormal body, the abnormal field values are distributed chaotically within the specified window, and the entropy value is large. Its mathematical expression formula is:

[0133]

[0134] Among them, S represents the induced polarization intensity parameter signal of a certain horizon, and entropy represents the singular spectrum entropy of relevant quantities within a moving window of size m*n. The normalized singular spectrum entropy method proposed by the present invention calculates the ratio of the singular spectrum entropy of the derivatives of the magnetic anomaly in three directions within the moving window as the NSVDE value of the center point within the window, and takes the maximum value of the NSVDE value as the boundary of the abnormal body.

[0135] Check the number of actual observed measurement points in the X-axis direction of the same survey line within the abnormal distribution range generated based on Step S4, and record the number as num. According to the specification requirements, on the same survey line, three consecutive points must appear for a reliable anomaly. When num ≥ 5, the normalized singular spectrum entropy method can be used to infer the boundary in the X-axis direction at this time. When 3 < num < 5, the boundary of the abnormal field can be directly used as the boundary of the abnormal body at this time; when the field conditions permit, further observations should be encrypted within the abnormal area, preferably so that there are more than 7 measurement points within the abnormal area, and then the normalized singular spectrum entropy method is used to infer the boundary in the X-axis direction.

[0136] Further, after the said Step 7, it further includes:

[0137] Step 8: Use a sphere as the forward model, and use the central depth position of the ore body and the distribution range of the ore body in the X-axis direction as spatial position constraints, and perform equivalent source fitting on the apparent resistivity and apparent polarization rate values of each measurement point to obtain the spatial distribution of the ore body. In Step 8, the diameter of the sphere can be determined according to the thickness information of each layer of the ore body in the preliminary geological data, and the diameter of the sphere must be less than the thickness of the previously known ore body.

[0138] Among them, Step 8 specifically includes:

[0139] S81: Based on the outer rectangle of the abnormal distribution range on the XOY plane in step S4, the horizontal area to be fitted is used for the equivalent source forward modeling; the upper surface of the area to be fitted in the vertical direction in the Z-axis direction is the terrain surface, and the lower bottom surface is the exploration depth interface; the three-dimensional area to be fitted is determined as the three-dimensional space of the equivalent source sphere distribution. The background physical property combination in the three-dimensional space is set according to the lithology type with the largest distribution area in the area, and is recorded as [rho c ,eta c ].

[0140] S82: Use a sphere as the forward model. The diameter of the sphere is determined based on the geological data available in the early stage to ensure that the diameter of the equivalent source sphere is smaller than the thickness of the ore body known in the early stage, so that there are a total of num_source spheres in the three-dimensional space where the equivalent source spheres are distributed.

[0141] S83: Based on the survey line near the abnormal center point of the polarization intensity signal contour map in step S4 as the abnormal center section, in the XOZ plane, the i maximum value points of the normalized total gradient within the abnormal range obtained in step S6, among which the pth extreme value corresponds to the z-axis coordinate z p 0; Assume that the initial thickness of the pth anomaly is Δz p , the distribution range of the pth abnormal body in the Z-axis direction is z p 0-Δz p and z p 0+Δz p In step S7, the distribution range of the abnormal body in the pth layer in the X-axis direction [x p 1,x p 2].

[0142] The Y-axis coordinate of the center point of the anomaly is YY, and the initial distribution range of the anomaly body in the pth layer in the Y-axis direction is set to

[0143] Thus, for the pth anomaly, the initial polarization cube range is determined to be

[0144] The num_mine spheres in the initial cube range of the distribution of the 1st to the pth anomaly are collectively referred to as high polarization equivalent sources. According to the pre-collected statistics of the resistivity and polarizability physical properties of each surrounding rock and ore body in the target area, the physical property combination of the spheres in the high polarization equivalent source area is the physical property combination of high polarization characteristics [rho1,eta1].

[0145] S84: The remaining (num_source-num_mine) spheres represent geological bodies without minerals, collectively referred to as non-polarizable equivalent sources. There may be multiple types of equivalent source property combinations. When setting, you can combine geological information and apparent resistivity morphology inference to reasonably set any of the (num_source-num_mine) spheres as any background equivalent source property combination, that is, [rho ii ,eta ii ], where eta ii Approximately equal to 0, where ii = 2, 3, …, c, so that after the spheres with the same physical property combination are connected, they can have reasonable geological significance.

[0146] S85: Perform forward simulation on the equivalent source distribution of this initial state.

[0147] In the specific forward modeling process, according to the spatial position,

[0148] The apparent resistivity value of the (i, j)th equivalent source at point (k, l) is:

[0149]

[0150] The apparent polarizability of the (i, j)th equivalent source at point (k, l) is:

[0151]

[0152] x=x ij -x kl

[0153] y=y ij -y kl

[0154] h0=hz kl

[0155]

[0156] Where, h represents the center depth of the three-dimensional polarization sphere, r0 represents the radius of the three-dimensional polarization sphere, is the relative resistivity, ρ2 and ρ c , are the resistivity of the equivalent sphere and the initial resistivity in the XOZ plane, respectively.

[0157] After accumulating and summing the apparent resistivity and apparent polarizability generated by each equivalent source within the observation range, the equivalent source forward modeling data, namely, the apparent resistivity file and the apparent polarizability file, are obtained.

[0158] S86: Compare the corresponding points of the simulated data with the observed data. When the mean square error of the simulated data of apparent resistivity and apparent polarizability at the same point and the observed data of apparent resistivity and apparent polarizability are both less than the threshold, the fitting is considered to be completed. At this time, the shape of the connected [rho1,eta1] property combination spheres can be considered to be the shape of the anomaly, which can be considered to be the spatial distribution of the ore body in geological interpretation.

[0159] If it is not less than the threshold, further change the number and distribution form of the modified high polarization equivalent sources and the non-polarization equivalent sources, and repeat the step S85 until the condition of being less than the threshold is met to obtain the spatial distribution of the ore body.

[0160] Based on the abnormal distribution of ladder parameter signals in induced polarization and the three-dimensional normalized total gradient value, the present invention accurately determines the central depth position of the ore body, improves the understanding and recognition ability of the internal structure of the ore body, and introduces the normalized singular spectrum entropy method to determine the boundary of the ore body, which can further reveal the shape and scope of the ore body and provide a reference basis for subsequent exploration and mining.

[0161] As another specific embodiment of the present invention, the method for rapid locating an ore body based on a ladder in time domain spectral induced polarization in the present invention comprises:

[0162] S0: Search for geological and geophysical data within the collection area and define exploration coordinates

[0163] (1) According to the geological data in the area, the mineralized alteration anomaly belt is obviously northwest-oriented, indicating that the silicification in the northwest-oriented dense joint mineralization is widespread and strong. Therefore, the survey line is expected to be vertical, the survey line direction is northeast, and the survey line direction is 128°. The exploration coordinate system is defined as follows: the vertical structural direction, that is, the 128° direction along the survey line, is the x direction; the vertical survey line, that is, the 38° direction, is recorded as the y direction; the vertical downward direction is recorded as the z direction.

[0164] (2) Collect samples of electrical properties measured by previous researchers in the work area, such as Figure 2 The physical properties of each rock and mineral were reasonably generalized in the form of average values, and typical values ​​of the surrounding rock and ore body were obtained.

[0165] Resistivity characteristics: The resistivity of gneiss and monzonite is generally about 1500 ohm-meter, the average resistivity of quartz monzonite is 4500 ohm-meter, and the gold-bearing quartz vein ore in the northwest-trending mineralized alteration body has a higher resistivity due to the strong degree of silicon mineralization, with an average value of 7000 ohm-meter;

[0166] Polarizability characteristics: The polarizability of various surrounding rocks is less than 2%, and the polarizability of quartz vein ore bodies exceeds 5%.

[0167] Based on the geological information and the results of geophysical electrical specimen measurement, the resistivity and polarizability combinations of the three main geological bodies in the area can be summarized as follows: 1 high-resistance and high-polarizability body - quartz vein ore [7000, 5%], 2 medium-resistance and low-polarizability body - quartz monzonite [4500, 1%], 3 low-resistance and low-polarizability body - gneiss, (hornblende) monzonite [1500, 1%]. The quartz vein ore is the ore body corresponding to this embodiment.

[0168] S1: Use a high-power transmitter to transmit positive and negative power supply signals, receive signals to collect induced polarization mid-level scanning data, record the elevation values ​​of the measuring points, and calculate the apparent resistivity and apparent polarizability values ​​of each measuring point.

[0169] The measurement parameters of the IP intermediate ladder are as follows: line spacing of 200 meters, point spacing of 40 meters, and pole spacing of 40 meters. The length of a single line is 1 km, and the length of the intermediate ladder AB emission line is 2 km. In the actual observation process, observations are made in a one-by-five or one-by-four manner. The measurement point is measured in the middle of the AB emission line. During the entire observation process of this area, the power supply current is stably maintained at 6A, which is far beyond the current standard required by the specification and fully supports the normal operation of IP. See the schematic diagram of the spatial distribution of the ABMN electrodes of the intermediate gradient device for details. Figure 3 .

[0170] In terms of technical parameters, the IP ladder is operated with a power supply cycle of 32s, and a total of 18 time windows are collected. The width of the first five time windows is 40ms. Since the signal gradually decreases, the time window width increases to 80ms, 120ms, and the maximum time window width is 160ms. Because the overall value of the polarization rate in the area is small, the increase in delay time will further reduce the value of the polarization rate, which is not conducive to the identification and confirmation of the polarized body. Therefore, combined with the test results, the delay time is selected as 80ms. The specific time windows are shown in Table 1.

[0171] Table 1. Acquisition time window table for ladder measurement in IP

[0172] Serial number 1 2 3 4 5 6 7 8 9 Start time 80 120 160 200 320 400 480 560 640 End time 120 160 200 320 400 480 560 640 720 Serial number 10 11 12 13 14 15 16 17 18 Start time 720 800 960 1120 1280 1440 1600 1760 1920 End time 800 960 1120 1280 1440 1600 1760 1920 2080

[0173] During the acquisition, the parameters recorded by the receiving end are the spatial position and elevation of the receiving MN electrode of the measuring point, the primary field voltage U0, and the secondary field voltage ΔU2 after power failure; the transmitting end records the spatial position of AB and the transmitting current I. After the recording is completed, the recorded data is comprehensively sorted out using an excel sheet. It is worth noting that with the development of current geophysical exploration instruments and equipment, some mathematical operations can be automatically completed during the data acquisition process, and some physical quantities generated by calculations can be directly copied from the instruments and equipment, such as IP and IP1-IP18. Six points on a certain measuring line in actual field observations are used as an example to illustrate, as shown in Table 2.

[0174] Table 2 Original record table of a certain survey line in actual field observation

[0175]

[0176]

[0177]

[0178]

[0179] The rightmost column in the table is the header. The Date and Time rows record the date and time of the observation. TxCurr is the transmitting current, which is constant at 6A in this line. The eight rows AxBxMxNxAyByMyNy represent the actual geodetic coordinate positions of the four electrode points, the two transmitting electrodes AB and the two receiving electrodes MN. The Vp row is the primary field voltage U0. The Cres row is the grounding resistance information obtained during the measurement at the receiving end. These physical quantities are generated by original field observations. The remaining physical quantities are generated by calculation; GeoX and GeoY are the midpoints of the receiving electrode MN, which are the geodetic coordinates of the actual exploration points, generated by simply summing the xy values ​​corresponding to the MxNxMyNy data and multiplying them by 0.5; K is the device coefficient generated by the actual coordinate calculation, and the calculation formula is as follows

[0180]

[0181] Where x0 is the horizontal coordinate position of the midpoint of MN, and y0 is the vertical coordinate position of the midpoint of MN. The origin of the coordinate is taken at the midpoint of AB. AB is the length of the transmitting electrode, and MN is the length of the receiving electrode. The specific calculation formula is as follows:

[0182]

[0183]

[0184] AppRes is the apparent resistivity obtained, and the formula is as follows

[0185]

[0186] Rows IP1-IP18 are the ratios of the secondary field voltage to the primary field voltage recorded in each time window, which is called the charging rate, i.e. ΔU2 / U0; row IP is the ratio of the sum of the results from IP1 to IP18 to the duration of all time windows from IP1 to IP18, i.e. the apparent polarizability.

[0187] S2: The Cole-Cole model is selected as the theoretical model. Based on the optimization algorithm, the apparent parameters of each measuring point, namely the polarization intensity m, time constant tau and frequency correlation coefficient c, are obtained point by point according to the charging rate information corresponding to each time window of the input induced polarization ladder, and the noise data that does not conform to the geophysical laws are removed point by point.

[0188] The Cole-Cole model is the most classic forward model of the induced polarization method. The attenuation curves of each time window that we actually observed can be considered as a comprehensive reflection of underground mine-induced electrical anomalies and human noise. The human noise is random and instantaneous, while the mine-induced electrical anomalies are stable and long-term.

[0189] The Cole-Cole model formula for complex resistivity of the induced polarization model proposed by Pelton is widely used in the industry:

[0190]

[0191] Among them, ρ0 is the zero-frequency resistivity, ω is the angular frequency, m is the polarizability, τ is the time constant, and c is the frequency correlation coefficient.

[0192] The above formula is a frequency domain formula. In order to complete the "frequency-time" conversion of the complex resistivity of the Cole-Cole model, the discrete form of the digital filtering algorithm proposed by Guptasarma is adopted. Based on the digital filtering algorithm, 21 filter coefficients are used to obtain the time-varying resistivity numbered with time in a cumulative manner.

[0193]

[0194] in

[0195]

[0196]

[0197] In addition, note that According to the above formula, the polarization intensity m, the time constant tau and the frequency correlation coefficient c can be calculated using the optimized algorithm based on the time t of each time window and the apparent polarization rate.

[0198] a r is the horizontal coordinate corresponding to the 21 filter coefficients, φ r is with a r The corresponding filter coefficients are shown in Table 3.

[0199] Table 3 Filter coefficients of digital filter algorithm

[0200] r Filter abscissa Filter coefficient r Filter abscissa Filter coefficient 1 -3.82704 0.000349998 12 -0.956480 0.042972904 2 -3.56608 -0.000418371 13 -0.695520 0.075423603 3 -3.30512 0.000772828 14 -0.434560 0.139346367 4 -3.04416 -0.000171356 15 -0.173600 0.234486236 5 -2.78320 0.001022172 16 -0.956480 0.366178323 6 -2.52224 0.000897638 17 0.34832 0.284615486 7 -2.26128 0.002208974 18 0.60928 -0.235691748 8 -2.00032 0.003844944 19 0.87024 0.046994188 9 -1.73939 0.006809040 20 1.13120 -0.0005901946 10 -1.47840 0.013029162 21 1.39216 0.000570165 11 -1.21744 0.022661391

[0201] The global optimization algorithm of 1stOpt software is used to obtain the three polarization rate parameters of the frequency domain induced polarization method, and the essential information reflecting induced polarization is mined, so as to achieve the purpose of identifying underground polarization anomalies. In this process, the observation data of each measuring point is guaranteed to conform to the Cole-Cole theoretical model, and the noise data that does not conform to the geophysical laws is removed point by point.

[0202] The input of the program is the charging rate data of IP1-IP18 and the corresponding time window. The output is the polarization intensity m, time constant tau and frequency correlation coefficient c calculated by the optimized algorithm. The global optimization algorithm avoids the problem of assigning reasonable initial values ​​when calculating the traditional gradient descent algorithm, and can ensure that there is always a solution and the stability of the solution.

[0203] After inversion fitting, the three parameters m, tau, and c of the Cole-Cole model can be obtained. The error between the fitting result and the original observation can be considered to be the influence of noise. Therefore, after fitting, it can be considered that the observed data of the polarizability have been denoised. The parameter m of the Cole-Cole model obtained represents the degree of polarization. In the subsequent mapping process, this parameter is used instead of the observed polarizability (IP) to draw a graph to represent the degree of polarization. The Cole-Cole model is used to fit the measured induced polarization curve, such as Figure 4 As shown in the figure, the blue dots are the original observation data, and the red dots are the fitting results of the Cole-Cole model. The fitting accuracy is over 99%, and the corresponding best estimated parameters m, c, and tau are 0.320, 0.8889, and 0.023, respectively. The difference between the blue and red dots is caused by noise, so the process of back-calculating the three parameters such as polarization intensity from the polarizability is itself a denoising process.

[0204] S3: After gridding the apparent resistivity and polarization intensity m of the IP gradient, the original contour map of the IP gradient parameters on the XOY plane is drawn respectively. According to the GeoX and GeoY calculated above, as well as the survey line orientation, the exploration coordinate system is rotated to the east-west direction of the survey line x and the north-south direction of the perpendicular survey line y based on trigonometric function calculations.

[0205] The rotated coordinate system is used for calculation in gridding and data processing, while the original coordinate system is still used for interpretation. The specific conversion relationship between the two coordinate systems is as follows;

[0206] x'=xcosθ-ysinθ

[0207] y'=ycosθ+xsinθ

[0208] x=(x'-y'tanθ)sinθ

[0209]

[0210] Among them, x and y are the exploration geodetic coordinates when the original actual data was collected, x′ and y′ are the coordinates after rotation, and the counterclockwise rotation angle is θ.

[0211] The apparent resistivity ρ and polarization intensity m of the IP gradient are gridded in the rotated exploration coordinate system, and the original contour map of the IP gradient parameters on the XOY plane is drawn; the elevation data of the gridded survey area is consistent with the grid of the IP gradient parameter contour map, and the topographic contour map of the survey area is generated; the gridded file is saved as a GRD file in the Surfer6 Text Grid format.

[0212] In this example, for the sake of confidentiality and to facilitate the subsequent simplification of the description, the upper left corner of the rotated exploration coordinate system is further set as the 0 point of the horizontal and vertical axes. Figure 5 . Figure 5 From left to right are the original contour maps of apparent resistivity, apparent polarizability, and apparent polarization intensity.

[0213] S4: It is necessary to remove the related noise and separate the signal and noise according to the data relationship between points. Generate a signal contour map.

[0214] SA fractal filtering is used for signal-to-noise separation. SA fractal filtering is based on the concept of generalized self-similarity and converts data from the spatial domain to the frequency domain to analyze its fractal characteristics. The relationship between the spectrum energy density and area of ​​the SA fractal method is:

[0215] A(≥S)∝S -β

[0216] In the formula, S is a certain threshold of energy spectrum density; A (≥ S) represents the area where the energy spectrum density is greater than the threshold S; β is the anisotropic multidimensional fractal index. The relationship between S and A is that as the S value increases, the corresponding A decreases, and its change rule mainly depends on β, which is the anisotropic multidimensional fractal index.

[0217] Taking the polarization intensity m as an example, it is used as input data, based on the double logarithmic relationship between energy spectrum density and area under the SA method, as follows Figure 6 It can be seen that it can be divided into three sections. Linear fitting is performed on the three sections respectively, and the fitting results are shown in Table 4. The first section is considered to be the reflection and embodiment of noise. The contour map corresponding to the noise separated by this method is Figure 7The noise contour results of m separated by this method have a high degree of discreteness, and there is no obvious high and low value area, which reflects the characteristics of random interference. The second and third sections are considered to be the reflection of the anomaly caused by the underground geological body. In the same way, the apparent resistivity is processed to obtain the signal contour map of apparent resistivity, apparent polarizability and apparent polarizability, as shown in Figure 8 . Figure 8 From left to right in the figure are the contour diagrams of the apparent resistivity, apparent polarizability, and apparent polarization intensity signals of the gradient parameters in the induced polarization.

[0218] Table 4 Linear fitting parameters under S-A method

[0219] Index Log(Value) Value Slope Intercept R** 0 -2.692 0.002 -0.26 2.786 0.92 1 -0.492 0.322 -0.463 2.648 0.986 2 2.381 240.563 -0.367 2.053 0.853

[0220] S5 extension: convert the retained IP gradient parameter signals from the terrain surface to each underground elevation along the Z axis in sequence, and obtain the contour map of the retained IP gradient parameter signals at each elevation.

[0221] The downward extension process is specifically divided into two situations. The lowest and highest elevations in the scanned area are between 1616 and 1895 meters, with a height difference of nearly 300 meters. In the case of such a large difference in elevation fluctuations, when extending directly from the terrain surface to the plane, from the highest point to the altitude of 1600 meters, there will be a blank area with an elevation of 300 meters that is difficult to obtain information on. Therefore, for areas below 1600 meters above sea level, the extension is made from the terrain surface to the plane at the corresponding altitude; for areas from the terrain surface to an altitude of 1600 meters, the extension is made downward for 300 meters in a curved surface manner. This ensures that all space from the terrain surface downward has downward extension information. Therefore, the extension method selected is "A method for analytical extension of gravity and magnetic potential fields based on the spatial domain". The output result is a downward extension contour map of the ladder parameter signal at each depth / elevation. Taking the polarization intensity signal data as an example, several depths and elevations are extended downward, respectively, as shown below. Fig. 9 and Fig.10 .

[0222] Fig. 9 a is the measured data of regional elevation, Fig. 9 c is the polarization intensity signal data after denoising by the SA method, Fig. 9 e. Fig. 9 b. Fig. 9 d. Fig. 9 e are the downward extension contour maps obtained by extending the ground elevation to depths of 100, 150, 200, and 300 meters respectively.

[0223] Fig.10a, c, e, g, i, b, d, f, h, and j are the downward extension results of elevations of 1550, 1500, 1450, 1400, 1350, 1300, 1250, 1200, 1100, and 1000 meters respectively.

[0224] S6: Depth of abnormal body: The three-dimensional normalized total gradient method is used to obtain the depth of the abnormal body that causes the abnormality.

[0225] 3D Normalized Total Gradient NTG z The calculation formula is:

[0226]

[0227] Among them, S xz , S yz , S zz The ladder parameter signal S at the height z is marked respectively z are the potential field gradient components along the x, y, and z directions.

[0228] The output result is the three-dimensional normalized total gradient of the gradient parameter signal of the IP at each depth / elevation. Fig.11 and Fig.12 .

[0229] Fig.11 a, b, c, and d represent the three-dimensional normalized total gradient contour maps of the polarization intensity m at depths of 100, 150, 200, and 300 meters, respectively.

[0230] Fig.12 a, c, e, g, i, b, d, f, h, and j are three-dimensional normalized total gradient contour maps at elevations of 1550, 1500, 1450, 1400, 1350, 1300, 1250, 1200, 1100, and 1000 meters above sea level, respectively.

[0231] Combined with the signal contour map generated in step S4, the survey line passing through the center of the anomaly in the signal contour map is extracted on the Y axis, and the three-dimensional normalized total gradient data of the induced polarization gradient parameter signal at each depth / elevation is combined to draw a plane map of the three-dimensional normalized total gradient on the XOZ plane. The maximum value in the Z-axis direction is extracted, and its maximum value represents the depth / elevation value of the center of the anomaly body causing the anomaly source.

[0232] In this example, in the XOY plane, when Y=400 meters, an anomaly can be seen on the signal contour map. On the XOZ plane, the three-dimensional normalized total gradient data of each depth and elevation at X=700 near the center of the anomaly is analyzed. The data tables are shown in Tables 5 and 6.

[0233] Table 5 Calculation results of three-dimensional normalized total gradient at each elevation at X=700, Y=400

[0234] X Y NTG Height(m) 700 400 6.1583 1500 700 400 5.95727 1450 700 400 6.01803 1400 700 400 5.83811 1350 700 400 5.90198 1300 700 400 5.84027 1200 700 400 5.8302 1100 700 400 5.72865 1000

[0235] Table 6 Calculation results of three-dimensional normalized total gradient at each depth at X=700, Y=400

[0236]

[0237] As shown in Tables 5 and 6, the three-dimensional normalized total gradient data has a maximum value at a depth of -150 meters and an altitude of 1400 meters, indicating that there are two abnormal bodies distributed at a depth of 150 meters underground and an altitude of 1400 meters. The center of the abnormal body is at a depth of 150 meters underground and an altitude of 1400 meters. It is also known that the altitude of the anomaly is about 1850 meters, which can be converted to the center of the abnormal body at a depth of 150 meters underground and a depth of about 450 meters underground. There are two ore bodies distributed at a depth of 150 meters underground and a depth of 450 meters underground. The ore bodies have the characteristics of high polarization and are inferred to be ore bodies. Drill holes were set up at abnormal locations for drilling verification. Quartz vein gold ore bodies were found at 144.31-145.53m, with a thickness of 1.22m and a gold grade of 1.17g / t; and silicon-potassium altered rock type mineralization bodies were found at 420m and 450m, with a grade of about 0.5g / t. The drilling verification is highly consistent with the processing and analysis results of the above-mentioned induced polarization ladder data, which proves the practicality of this method. Because the drilling cataloging data is confidential, the corresponding cataloging pictures and test results cannot be provided.

[0238] As another specific implementation of the present invention, after S6, the present invention further includes: performing X-axis direction boundary recognition on the ladder anomaly in the simulated polarization.

[0239] The forward calculation is performed with a cube model, and the intervals of the measuring points of the survey line are 100 meters and 40 meters; the transmitting pole distance is 2000 meters, and the receiving pole distance is 40 meters, which is consistent with Example 1. The apparent resistivity and apparent polarizability anomalies caused by the cube at a plane altitude of z = -120 meters are obtained by numerical simulation. The distribution ranges of cube 1 in the three directions of xyz are [-60, 60], [-140, 140], [-260, -240]; the distribution ranges of cube 2 in the three directions of xyz are [-100, 100], [-140, 140], [-260, -240]; the distribution ranges of cube 3 in the three directions of xyz are [-140, 140], [-140, 140], [-260, -240]. In other words, in the x-axis measuring point direction, for the same measuring line, cube 1 simulates the situation where three measuring points pass through the anomaly, cube 2 simulates the situation where five measuring points pass through the anomaly, and cube 3 simulates the situation where seven measuring points pass through the anomaly.

[0240] For the three cubes, based on the measured point data on each simulated survey line, the point distance is used as the unit size for gridding, that is, 40 meters is used as the grid, the apparent polarizability is gridded, and a grid file is generated.

[0241] For the three cubes, the apparent polarizability grid files are used as input to calculate and generate the normalized singular spectrum entropy grid files. The window size in the x-axis direction is 3, and the window size in the y-axis direction is 5.

[0242] Draw the true boundaries of the three cubes, the contour map of apparent polarizability, and the contour map of normalized singular spectral entropy. Fig.13 . Fig.13 In the figure, the left side shows the contour map of apparent polarizability, and the right side shows the contour map of normalized singular spectral entropy; from top to bottom, they are cubes 1, 2, and 3.

[0243] from Fig.13 It can be seen that for cube 1 with three measuring points passing through the anomaly, at this time, in the X-axis direction, the real boundary and the apparent polarizability anomaly range roughly coincide with each other, and the agreement is relatively good. According to the inferred distribution range of the anomaly in the X-axis direction, the real boundary in the X-axis direction coincides with the maximum value of the normalized singular spectrum entropy. According to the inferred distribution range of the anomaly in the X-axis direction, the real boundary in the X-axis direction corresponds to the real boundary. For cube 2 with five measuring points passing through the anomaly, the real boundary does not agree well with the apparent polarizability anomaly range and the maximum value of the normalized singular spectrum entropy: inferred from the apparent polarizability anomaly range, the boundary is [-60,60]; inferred from the maximum value of the normalized singular spectrum entropy, within the boundary range of [-100,100], the maximum value of the normalized singular spectrum entropy has a certain manifestation, but the boundary characteristics are not obvious enough. The maximum value method of the normalized singular spectrum entropy is available but not easy to use. At this time, when field conditions permit, the number of measuring points should be increased to more than 7, and then the boundary recognition method of normalized singular spectral entropy should be used to infer the boundary in the X-axis direction.

[0244] The present invention also provides a computer-readable storage medium having a computer program stored thereon, characterized in that when the computer program is executed by a processor, the steps of the above-mentioned method for rapid locating an ore body based on ladders in time domain spectral induced polarization are implemented. Compared with the prior art, the beneficial effects of the computer-readable storage medium provided by the present invention are the same as the beneficial effects of the method for rapid locating an ore body based on ladders in time domain spectral induced polarization described in the above-mentioned technical scheme, which will not be repeated here.

[0245] This article uses specific examples to illustrate the principles and implementation methods of the present invention. The above examples are only used to help understand the method and core ideas of the present invention. At the same time, for those skilled in the art, according to the ideas of the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A method for rapid ore body positioning based on time domain spectrum induced polarization ladder, characterized in that: include: Step 1: Construct an exploration coordinate system with any point in the study area as the origin, and use the induced polarization ladder method to calculate the apparent resistivity and apparent polarizability values ​​of each measuring point in the study area; Step 2: Use the Cole-Cole model to construct a time-domain spectral induced polarization forward model, and input the apparent polarizability value of each measuring point into the time-domain spectral induced polarization forward model to generate the polarization intensity parameters of each measuring point; Step 3: Rotate and grid the exploration coordinate system based on the azimuth of the survey line, and draw contour maps of induced polarization gradient parameters and topographic contour maps on the gridded XOY plane; the induced polarization gradient parameters include: apparent resistivity, apparent polarizability value and polarization intensity parameter; Step 4: using the energy spectrum density area fractal method to perform signal-noise separation on the contour map of the induced polarization ladder parameter to obtain a denoised induced polarization ladder parameter signal file, and determining the abnormal distribution in the contour map of the denoised induced polarization ladder parameter signal file; Step 5: Generate the downward extended induced polarization intensity signal files of multiple layers of underground induced polarization based on the denoised induced polarization parameter signal files based on the spatial domain stable continuation algorithm; Step 6: Calculate the three-dimensional normalized total gradient values ​​of different layers below the measuring point of the survey line where the anomaly is located based on the abnormal distribution in the contour map of the downward extended induced polarization gradient polarization intensity signal file of multiple underground layers and the denoised induced polarization gradient parameter signal file, and determine the central depth position of the ore body based on the three-dimensional normalized total gradient value; Step 7: Determine the boundary of the ore body using the normalized singular spectral entropy method.

2. A method for rapid ore body positioning based on time domain spectrum induced polarization according to claim 1, characterized in that: In step 1, the direction perpendicular to the structural strike is used as the direction of the transmitting line and the receiving line, and the survey line is laid. At the same time, an exploration coordinate system is constructed with the direction along the survey line as the x direction, the direction perpendicular to the survey line as the y direction, and the vertical downward direction as the z direction.

3. The method for rapid ore body positioning based on time domain spectrum induced polarization according to claim 2, characterized in that: In step 1, the apparent resistivity calculation formula is: ρ a is the apparent resistivity of the measuring point, K is the device coefficient, x0 is the abscissa position of the midpoint of MN, y0 is the ordinate position of the midpoint of MN, AB is the length of the transmitting electrode, MN is the length of the receiving electrode, and U0 is the primary field voltage; The calculation formula of apparent polarizability is: Among them, η a is the apparent polarization rate of the measuring point, and ΔU2 is the secondary field voltage after power failure.

4. A method for rapid ore body positioning based on time domain spectrum induced polarization according to claim 3, characterized in that: In step 2, based on the digital filtering algorithm, 21 filter coefficients are used to obtain the time domain spectral induced polarization forward model in an accumulative manner; wherein the time domain spectral induced polarization forward model is: Among them, ω r is the frequency variable of the real component of r from 1 to 21 corresponding to time t, a r is the horizontal coordinate corresponding to the 21 filter coefficients, φ r is with a r The corresponding filter coefficient, τ is the time constant, c is the frequency correlation coefficient, ρ0 is the zero-frequency resistivity, R r is the real component, I r Is the imaginary quantity.

5. A method for rapid ore body positioning based on time domain spectrum induced polarization according to claim 4, characterized in that: In step 3, the rotation of the exploration coordinate system is completed based on the azimuth of the survey line; wherein the conversion relationship formula of the coordinate system is: x′=xcosθ-ysinθ y′=ycosθ+xsinθ x′=(x′-y′tanθ)sinθ Wherein, x is the coordinate of the original exploration coordinate system on the x-axis, y is the coordinate of the original exploration coordinate system on the y-axis, x′ is the coordinate of the exploration coordinate system on the x-axis, and y′ is the coordinate of the rotated exploration coordinate system on the y-axis.

6. A method for rapid ore body positioning based on time domain spectrum induced polarization according to claim 5, characterized in that: The step 6: based on the abnormal distribution in the contour map of the downward extended induced polarization gradient polarization intensity signal file of multiple underground layers and the denoised induced polarization gradient parameter signal file, the three-dimensional normalized total gradient values ​​of different layers below the measuring point of the measuring line where the abnormality is located are calculated, and the central depth position of the ore body is determined based on the three-dimensional normalized total gradient value, including: Step 6.1: According to the abnormal distribution range in the isogram of the IP gradient parameter signal file after denoising, extract the survey line passing through the center of the anomaly on the Y-axis. On the XOZ plane where this survey line is located, interpolate the three-dimensional normalized total gradient values at different horizons in the Z direction to obtain the three-dimensional normalized total gradient plane file on the XOZ two-dimensional plane where the anomaly center is located; Step 6.2: Take the position of the maximum value in the three-dimensional normalized total gradient plane file on the Z-axis as the central depth position of the ore body.

7. A method for rapid ore body positioning based on time domain spectrum induced polarization according to claim 6, characterized in that: Step 7 mentioned above: Using the normalized singular spectrum entropy method to determine the boundary of the ore body, including: Step 7.1: Obtain the number num of actual observed measurement points on the same survey line in the X-axis direction at the central depth position of the ore body; Step 7.2: When 3 < num < 5, directly take the boundary of the abnormal field as the boundary of the ore body. When num ≥ 5, use the normalized singular spectrum entropy method to infer the boundary of the ore body in the X-axis direction to obtain the distribution range of the ore body in the X-axis direction; Step 7.3: Determine the boundary of the ore body according to the distribution range of the ore body in the X-axis direction.

8. A method for rapid ore body positioning based on time domain spectrum induced polarization according to claim 7, characterized in that: After Step 7 mentioned above, it further includes: Step 8: Using a sphere as the forward modeling model, taking the central depth position of the ore body and the distribution range of the ore body in the X-axis direction as spatial position constraints, perform equivalent source fitting on the apparent resistivity and apparent polarization rate values of each measurement point to obtain the spatial distribution of the ore body.

9. A method for rapid ore body positioning based on time domain spectrum induced polarization according to claim 8, characterized in that: In Step 8, the calculation formula for the apparent resistivity value of the (i, j)-th equivalent source at the (k, l)-th point is: The calculation formula for the apparent polarization rate value of the (i, j)-th equivalent source at the (k, l)-th point is: x=x ij -x kl y=y ij -and kl h0=hz kl in, represents the apparent resistivity value of the (i,j)th equivalent source at point (k,l), represents the apparent polarizability value of the (i, j)th equivalent source at the point (k, l), h represents the central burial depth of the three-dimensional polarization sphere, r0 represents the radius of the three-dimensional polarization sphere, is the relative resistivity, ρ2 and ρ c are the resistivity of the equivalent sphere and the initial resistivity in the XOZ plane, (x ij ,y ij ) represents the horizontal coordinate of the equivalent source position, (x kl ,y kl , z kl ) represents the coordinates of the observation point, η2 represents the equivalent source sphere polarizability, M ij Represents an intermediate variable.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, it implements the steps in a method for rapid positioning of an ore body based on time-domain spectral IP gradient as described in any one of claims 1-9.

Images