Transient electromagnetic three-dimensional inversion method based on vector field geometric compression method

Abstract

The invention discloses a transient electromagnetic three-dimensional inversion method based on a vector field geometric compression method, and the method comprises the steps: carrying out the self-adaptive compression of a forward modeling field through a control threshold value according to the geometric features of electromagnetic field data at all moments of a three-dimensional forward modeling calculation result, and only retaining key feature points reflecting the curve form of the change rate of a magnetic field. And in the accompanying forward modeling calculation process, the forward modeling field value at the required moment is recovered through a compression rule. Through the information-recoverable vector field time scale geometric compression method, the data storage requirement can be remarkably reduced, and then the limitation of small models and large storage is broken through. And large-scale three-dimensional inversion calculation is carried out under the condition of limited video memory resources.

Description

Technical Field

[0001] The present invention relates to the field of geophysical exploration, and particularly to a transient electromagnetic three-dimensional inversion method for an information-recoverable vector field time-scale geometric compression method, which is applicable to the field of transient electromagnetic three-dimensional inversion and aims to improve the three-dimensional inversion calculation ability of large-scale complex models. Background Art

[0002] In the electromagnetic field, a vector field that changes with time is an important basis for explaining the evolution law of the electromagnetic field in space. When considering the transient electromagnetic inversion problem, data in which the electromagnetic field changes with time is usually used to invert the spatial distribution of underground conductors, which is of great significance for understanding and interpreting the underground structure.

[0003] However, when performing transient electromagnetic three-dimensional inversion, a major challenge is faced: calculating the sensitivity matrix requires storing the forward field values at all times. The conventional storage method for forward calculation data faces problems such as a large amount of data and high storage requirements. As the model scale increases and the number of time steps increases, these challenges become more severe. In addition, since fast three-dimensional inversion algorithms require the use of a GPU acceleration card for accelerated calculation, an excessive number of time steps and a larger model scale will increase the video memory requirements. When the video memory resources are limited, it will lead to inability to calculate.

[0004] Compressing the original forward electromagnetic field data can reduce the storage space and bandwidth requirements, which is crucial for efficient data transmission and fast three-dimensional inversion in a distributed computing environment. Summary of the Invention

[0005] The purpose of the present invention is to overcome the above-mentioned deficiencies of the prior art and provide a transient electromagnetic three-dimensional inversion method based on a vector field geometric compression method, aiming to improve the three-dimensional inversion calculation ability of large-scale complex models.

[0006] To achieve the above purpose, the present invention adopts the following technical solutions: A transient electromagnetic three-dimensional inversion method based on a vector field geometric compression method includes the following steps: (1). Data acquisition: Using transient electromagnetic method equipment to acquire the time-domain response data of the electromagnetic field, and constructing an electromagnetic field data set containing multiple sampling times; (2). Constructing a three-dimensional inversion objective function; (3). Designing a forward equation according to the three-dimensional forward theory; (4). Information-recoverable vector field compression, compressing the electromagnetic field vector data, and extracting the key point data information on the electric field and magnetic field decay curves; the specific steps are as follows: (4-1). Extraction of vector curve feature points: For all the data in the electromagnetic field attenuation curve, take the first data point of the curve as the starting point, connect the starting point and the third data point to form a line segment, calculate the field value corresponding to the intermediate time through interpolation, and calculate the relative error; (4-2). Setting of compression error threshold: Preset the compression error threshold to control the error range during the compression process; (4-3). Judgment and retention of feature points: Compare the calculated relative error with the compression error threshold. If the error is greater than or equal to the threshold, retain the current data point as a feature point, take this feature point as the new starting point, connect it with the next data point in the vector curve sequence to form a new line segment and perform difference calculation, and compare the error; if it is less than the threshold, consider that this point is not a feature point, retain the original starting point, form a new line segment with the subsequent data points, perform interpolation calculation on multiple field values on the new line segment, calculate their relative errors with the corresponding points on the original curve, compare these relative errors with the compression error threshold, and if any relative error is greater than the compression error threshold, determine the previous data point corresponding to this error as a key feature point and retain it; (4-4). Curve reconstruction: Traverse all data points in sequence and connect these feature points to form a new curve; (4-5). Field value recovery: During the calculation of the product of the Jacobian matrix and the residual vector, for the calculation times between adjacent feature points, calculate the field value data at the required time through linear interpolation; (4-6). Information recovery verification: Restore the compressed data and verify whether the restored data can accurately reflect the characteristics of the original data; evaluate the effectiveness of the compression method by comparing the relative errors between the restored data and the original data.

[0007] In the step (1), the transient electromagnetic method device includes a transmitting coil with a diameter that is powered by a transmitter to excite a primary field, and a receiving coil laid on the ground or mounted on a drone.

[0008] In the step (1), the acquisition of the time-domain response data of the electromagnetic field includes supplying power to the transmitting coil and transmitting a primary field with a duty cycle of 1:1 underground; capturing the change rate of the secondary induced magnetic field generated by the underground eddy current through the receiving coil during the current turn-off period.

[0009] In the step (1), the electromagnetic field data set refers to the change rate of the secondary induced magnetic field generated by the underground eddy current.

[0010] In step (2), the transient electromagnetic three-dimensional inversion is an optimization process of an objective function; first, an initial model for perturbation is set, and then the initial model is updated through iterative calculations until the objective function reaches the convergence condition. At this time, it is considered that the forward modeling result and the actual observed data satisfy the fitting condition; in geophysical inversion, the objective function consists of two parts, namely: φ(m) = φ d + λφ ω （m） , (1) where φ d and φ ω represent the data fitting term and the model constraint term respectively. The subscript ω represents the dimension of the data, m represents the conductivity, λ is the regularization factor for balancing the data term and the model term.

[0011] In step (3), the forward equation is: Kx = b, (2) where K represents the coefficient matrix, x represents the electric and magnetic fields in the control equation of the three-dimensional forward algorithm, and b represents the right-hand side term.

[0012] Calculate J T r: , (3) where r is the residual vector between the observed data and the forward data, J is the Jacobian matrix, v is called the adjoint field, L represents the interpolation operator in time and space, T represents the transpose of a matrix or vector, , , in the time domain method, K represents the coefficient matrix, and the order of K is related to the model scale and the number of calculation times; The residual vector r and the electric and magnetic fields x in the control equation of the three-dimensional forward algorithm in equation (2) are solved through one forward modeling. Then, J T r is used as the source term, and v is calculated through one adjoint forward modeling. Finally, x and v are substituted into equation (3) to calculate J T r.

[0013] In step (4), when calculating J T r, the required electromagnetic field values are restored segment by segment according to the key feature points and data compression rules; the compression algorithm uses the Douglas-Peucker algorithm, the vertical distance limit method, or the angle limit method; the compression uses lossy compression or lossless compression.

[0014] In step (4-2), the compression error threshold is comprehensively considered according to the compression ratio and the scale accuracy. For example, based on the video memory size of the GPU graphics card used by the computing platform, a smaller compression threshold can be set when there is sufficient video memory space, and a larger compression threshold can be used otherwise. Although a smaller compression threshold occupies a larger video memory space, it will obtain a more accurate curve of the electromagnetic field varying with time. A larger compression threshold increases the error between the restored data and the original data while reducing the video memory requirement.

[0015] The beneficial effects of the present invention are as follows: The present invention provides a geometric compression method for the time scale of a vector field with recoverable data information. According to the geometric characteristics of the electromagnetic field data at all moments of the three-dimensional forward calculation results, the forward field is adaptively compressed by controlling the threshold, and only the key feature points reflecting the curve shape of the magnetic field change rate are retained. During the forward calculation process, the forward field values at the required moments are restored through the compression rules. Through this geometric compression method for the time scale of a vector field with recoverable information, the data storage requirement can be significantly reduced, thereby breaking the limitation of "small model, large storage". Large-scale three-dimensional inversion calculations can be carried out under limited video memory resources. The specific description is as follows: 1. Application of the information-recoverable vector field compression technology in three-dimensional inversion: The present invention provides an information-recoverable vector field compression method, which can, while compressing the electromagnetic field vector data, retain the data information of the key feature points of the electromagnetic field decay curve, ensuring that during the inversion process, the forward field values satisfying a certain accuracy at all moments can be restored segment by segment, reducing data redundancy, and improving the three-dimensional inversion calculation ability for large-scale complex models.

[0016] 2. Extraction of key feature points of the vector curve: The key feature points of the vector curve are extracted through a specific algorithm, and these feature points play a key role in the compression and restoration processes, ensuring the physical or geological interpretability of the data.

[0017] 3. Feature point judgment and retention mechanism: The present invention provides a feature point judgment and retention mechanism based on the comparison between the relative error and the compression error threshold to determine which data points should be regarded as key feature points and retained.

[0018] 4. Curve reconstruction and field value restoration technology: The new compressed curve is formed by connecting the extracted feature points, realizing the geometric compression of the vector data. During the forward calculation process, the field value data at all moments is obtained through interpolation calculation, and the maximum relative error between the interpolated data and the original forward field value is controlled. This method can effectively reduce the storage requirement of the forward electromagnetic field data during the calculation process, while maintaining the recoverability of the forward data within a certain error range and the accuracy of the inversion result. Description of the Drawings

[0019] Figure 1 It is a schematic diagram of a method for extracting curve feature points with recoverable information; Figure 2 It is an L-shaped low-resistivity anomaly model; Figure 3 It is the inversion result without data compression; Figures 4 - 7 They respectively correspond to the horizontal slices (z = 110m) of the inversion results in the xoy plane when the compression error thresholds are 2%, 6%, 10% and 20%. Specific implementation manner

[0020] The present invention will be further described below in conjunction with the drawings and embodiments.

[0021] The structures, ratios, sizes, etc. shown in the drawings of this specification are only used to cooperate with the content disclosed in the specification for those skilled in this technology to understand and read, and are not used to limit the limiting conditions under which the present invention can be implemented. Therefore, they do not have technical substance significance. Any modification of the structure, change of the proportional relationship or adjustment of the size, without affecting the effects that the present invention can produce and the purposes that can be achieved, should still fall within the scope covered by the technical content disclosed in the present invention. At the same time, the terms such as "upper", "lower", "left", "right", "middle" and "one" cited in this specification are only for the convenience of clear narration and are not used to limit the scope under which the present invention can be implemented. The change or adjustment of their relative relationships, without substantial change in the technical content, should also be regarded as the scope under which the present invention can be implemented.

[0022] A transient electromagnetic three-dimensional inversion method based on a vector field geometric compression method, comprising the following steps: (1). Data acquisition: Using transient electromagnetic method equipment to collect the time-domain response data of the electromagnetic field, and constructing an electromagnetic field data set containing multiple sampling moments; The transient electromagnetic method equipment includes a transmitting coil with a diameter that supplies power through the transmitter and excites the primary field, and a receiving coil laid on the ground or mounted by a drone.

[0023] Collecting the time-domain response data of the electromagnetic field includes supplying power to the transmitting coil and transmitting a primary field with a duty cycle of 1:1 underground; capturing the change rate of the secondary induced magnetic field generated by the underground eddy current through the receiving coil during the current turn-off period.

[0024] The electromagnetic field data set refers to the change rate of the secondary induced magnetic field generated by the underground eddy current.

[0025] (2). Construct a 3D inversion objective function; 3D transient electromagnetic inversion is a process of optimizing the objective function; First, set an initial model for perturbation, and then update the initial model through iterative calculations until the objective function reaches the convergence condition. At this time, it is considered that the forward modeling result and the actual observed data satisfy the fitting condition; In geophysical inversion, the objective function consists of two parts, namely: φ(m) = φ d + λφ ω （m） , (1) where, φ d and φ ω represent the data fitting term and the model constraint term respectively. The subscript ω represents the dimension of the data, m represents the conductivity, λ is the regularization factor for balancing the data term and the model term.

[0026] (3). According to the 3D forward modeling theory, design the forward modeling equation as: Kx = b, (2) In the formula, K represents the coefficient matrix, x represents the electric and magnetic fields in the control equation of the 3D forward modeling algorithm, and b represents the right - hand side term.

[0027] Calculate J T r: , (3) where, r is the residual vector between the observed data and the forward modeling data, J is the Jacobian matrix, v is called the adjoint field, L represents the interpolation operator in time and space, T represents the transpose of a matrix or vector, , , in the time - domain method, K represents the coefficient matrix, and the order of K is related to the model scale and the number of calculation times; Solve the residual vector r and the electric and magnetic fields x in the control equation of the 3D forward modeling algorithm in equation (2) through one - time forward modeling. Then, take J T r as the source term, calculate v through one - time adjoint forward modeling, and finally substitute x and v into formula (3) to calculate J T r.

[0028] (4). Compress the information - recoverable vector field, compress the electromagnetic field vector data, and extract the key - point data information on the attenuation curves of the electric and magnetic fields; When calculating J T r, recover the required electromagnetic field values section by section according to the key feature points and data compression rules; The compression algorithm adopts the Douglas - Peucker algorithm, the vertical distance limit method, or the angle limit method; The compression adopts lossy compression or lossless compression.

[0029] The specific steps are as follows: (4-1). Extraction of vector curve feature points (as Figure 1 shown): For all data in the electromagnetic field attenuation curve, take the first data point of the curve as the starting point, connect the starting point and the third data point to form a line segment, calculate the field value corresponding to the intermediate time through interpolation, and calculate the relative error; (4-2). Setting the compression error threshold: Preset the compression error threshold to control the error range during the compression process; The compression error threshold is comprehensively considered according to the compression ratio and scale accuracy, such as the video memory size of the GPU graphics card used by the calculation platform. When there is enough video memory space, a smaller compression threshold can be set, and vice versa, a larger compression threshold can be used. Although a smaller compression threshold occupies a larger video memory space, it will obtain a more accurate curve of the electromagnetic field changing with time. A larger compression threshold increases the error between the restored data and the original data while reducing the video memory requirement.

[0030] (4-3). Judgment and retention of feature points: Compare the calculated relative error with the compression error threshold. If the error is greater than or equal to the threshold, retain the current data point as a feature point, take this feature point as the new starting point, connect it with the next data point in the vector curve sequence to form a new line segment and perform difference calculation, and compare the error; If it is less than the threshold, it is considered that this point is not a feature point, retain the original starting point, form a new line segment with the subsequent data points, perform interpolation calculation on multiple field values on the new line segment, and calculate their relative errors with the corresponding points on the original curve. Compare these relative errors with the compression error threshold. If any relative error is greater than the compression error threshold, determine the previous data point corresponding to this error as a key feature point and retain it; (4-4). Curve reconstruction: Traverse all data points in sequence and connect these feature points to form a new curve; (4-5). Field value restoration: During the calculation of the product of the Jacobian matrix and the residual vector, obtain the field value data at the required time by interpolating according to the corresponding time between adjacent feature points; (4-6). Information restoration verification: Restore the compressed data and verify whether the restored data can accurately reflect the characteristics of the original data; Evaluate the effectiveness of the compression method by comparing the relative error between the restored data and the original data. Embodiment

[0031] such as Figure 2As shown in the figure, the abnormal body model parameters are as follows: the size of the transmitting coil loop is 500m×500m, and the transmitting current is 1A. The resistivities of air and the background earth are 106Ω·m and 100Ω·m respectively. The burial depth of the upper surface of the abnormal body is 100m from the ground, and the resistivity of the abnormal body is 10Ω·m. The Dirichlet boundary condition is still applied to the model boundary. The core calculation area is divided into a uniform grid of 20m×20m×20m, and a non-uniform grid is used outside the calculation area. The grid expansion coefficient is 1.3. A total of 11×11 = 121 measuring points are evenly arranged in the loop source. The BEDS-FDTD algorithm is used to perform forward calculation on this model. The entire process requires a total of 12925 time steps, of which the Off-time stage is 6046 time steps. The forward data is taken as the response result within 10 -5 s~10 -2 s after the current is turned off. 28 time channel data are obtained according to the logarithmically equally spaced differences. 3% Gaussian noise is added to the forward data as the measured data for inversion simulation. Both the initial model and the reference model adopt a uniform half-space model with a resistivity of 100Ω·m. The initial regularization factor is 10, and the minimum regularization factor is 10-3. The computing platform uses AutoDL, the CPU is a 10 vCPU Intel Xeon Processor (Skylake, IBRS), the GPU acceleration card is an A100-PCIE-40GB (40GB), the video memory is 40GB, and the host memory is 72GB.

[0032] When the compression error threshold is set to 0, that is, the original forward response data is not compressed. At this time, the video memory occupied is 19.36GB, and this set of results is used as a control; when the compression error threshold is set to 1%, the video memory occupied is 10.08GB, and the reduction rate of video memory occupancy reaches 47.9%; when the compression error threshold is set to 3%, the video memory occupied is 7.47GB, and the reduction rate of video memory occupancy reaches 61.4%; when the compression error threshold is set to 7%, the video memory occupied is 5.79GB, and the reduction rate of video memory occupancy reaches 70.1%; when the compression error threshold is set to 20%, the video memory occupied is 4.4GB, and the reduction rate of video memory occupancy reaches 77.3%.

[0033] Figure 3 The inversion result without data compression is shown. The resistivity value of the abnormal body within the purple frame is 10Ω·m, which is equal to the resistivity of the low-resistance abnormal body set in the model. In addition, the L shape of the abnormal body is basically restored. Due to the influence of the current loop effect, a circular low-resistance area appears on the right side of the abnormal body.

[0034] Figures 4 - 7 Corresponding to the inversion results when the compression error thresholds are 2%, 6%, 10% and 20% respectively in xoyHorizontal slice of the plane (z = 110 m). Through observation, it can be seen that the inversion results under different compression error thresholds can all recover the characteristics of the abnormal body, verifying the effectiveness of the vector field compression method proposed by the present invention.

[0035] Although the specific implementation manners of the present invention have been described above in conjunction with the accompanying drawings, they are not limitations on the protection scope of the present invention. Those skilled in the art should understand that various modifications or deformations that can be made without creative efforts on the basis of the technical solutions of the present invention are still within the protection scope of the present invention.

Claims

1. A transient electromagnetic three-dimensional inversion method based on a vector field geometric compression method, characterized in that: The following steps are involved: (1) Data acquisition: Use transient electromagnetic method equipment to collect time domain response data of the electromagnetic field and construct an electromagnetic field data set containing multiple sampling moments; (2). Construct the three-dimensional inversion objective function; (3) Design the forward modeling equations based on the three-dimensional forward modeling theory; (4) Information recovery vector field compression, compressing electromagnetic field vector data, and extracting key point data information on the electric field and magnetic field attenuation curves; the specific steps are as follows: (4-1). Extraction of feature points of vector curves: For all data in the electromagnetic field attenuation curve, take the first data point of the curve as the starting point, connect the starting point and the third data point to form a line segment, calculate the field value corresponding to the intermediate moment by interpolation, and calculate the relative error; (4-2). Compression error threshold setting: pre-set the compression error threshold to control the error range during the compression process; (4-3). Feature point judgment and retention: Compare the calculated relative error with the compression error threshold. If the error is greater than or equal to the threshold, the current data point is retained as a feature point, and the feature point is used as a new starting point to connect with the next data point in the vector curve sequence to form a new line segment and perform difference calculation to compare the error. If it is less than the threshold, the point is considered not a feature point, the original starting point is retained, and a new line segment is formed with subsequent data points. Multiple field values ​​on the new line segment are interpolated and their relative errors with the corresponding points of the original curve are calculated. These relative errors are compared with the compression error threshold. If any relative error is greater than the compression error threshold, the previous data point corresponding to the error is determined as a key feature point and retained. (4-4). Curve reconstruction: traverse all data points in sequence and connect these feature points to form a new curve; (4-5) Field value recovery: In the process of calculating the product of the Jacobian matrix and the residual vector, the field value data at the required time is obtained by interpolating the adjacent feature points at the corresponding time. (4-6). Information recovery verification: The compressed data is restored to verify whether the restored data can accurately reflect the characteristics of the original data; The effectiveness of the compression method is evaluated by comparing the relative error between the restored data and the original data.

2. The transient electromagnetic three-dimensional inversion method based on the vector field geometric compression method according to claim 1 is characterized in that: In the step (1), the transient electromagnetic method device includes a transmitting coil with a diameter of 0.040 mm / s which is powered by a transmitter and excites a primary field, and a receiving coil which is laid on the ground or mounted on a drone.

3. The transient electromagnetic three-dimensional inversion method based on the vector field geometry compression method as claimed in claim 2 is characterized in that: In the step (1), the acquisition of the time domain response data of the electromagnetic field includes powering the transmitting coil and transmitting it underground to generate a primary field with a duty cycle of 1:1; and capturing the rate of change of the secondary induced magnetic field generated by the underground eddy current through the receiving coil during the current off period.

4. The transient electromagnetic three-dimensional inversion method based on the vector field geometric compression method as claimed in claim 3 is characterized in that: In the step (1), the electromagnetic field data set refers to the rate of change of the secondary induced magnetic field generated by the underground eddy current.

5. The transient electromagnetic three-dimensional inversion method based on the vector field geometric compression method according to claim 1 is characterized in that: In step (2), transient electromagnetic three-dimensional inversion is an objective function optimization process; first, an initial model for disturbance is set, and then the initial model is updated through iterative calculation until the objective function reaches the convergence condition, at which time it is considered that the forward modeling result and the actual observation data meet the fitting condition; in geophysical inversion, the objective function consists of two parts, namely: φ (m) = φ d + λφ ω （m） ， （1） in, φ d and φ ω Represent data fitting terms and model constraint terms, respectively. ω represents the dimension of the data, m represents the conductivity, λ is the regularization factor that balances the data term and the model term.

6. The transient electromagnetic three-dimensional inversion method based on the vector field geometric compression method as claimed in claim 5 is characterized in that: In step (3), the forward equation is: Kx=b, (2) Where K represents the coefficient matrix, x represents the electric and magnetic fields in the three-dimensional forward algorithm control equation, and b represents the right-hand side term.

7. The transient electromagnetic three-dimensional inversion method based on the vector field geometric compression method according to claim 6 is characterized in that: Calculate J T r: , （3） Among them, r is the residual vector of the observed data and the forward data, J is the Jacobian matrix, v is called the adjoint field, L represents the interpolation operator in time and space, T represents the transpose of the matrix or vector, , , in the time domain method, K represents the coefficient matrix; the residual vector r and the electric and magnetic fields x in the three-dimensional forward algorithm control equations in equation (2) are solved by forward modeling once, and then J T r is used as the source term, and v is calculated through a forward modeling. Finally, x and v are substituted into equation (3) to calculate J T r.

8. The transient electromagnetic three-dimensional inversion method based on the vector field geometric compression method according to claim 7 is characterized in that: In step (4), when calculating J T r, the required electromagnetic field values ​​are restored piece by piece according to the key characteristic points and data compression rules; the compression algorithm adopts Douglas-Peucker algorithm, vertical distance limit method or angle limit method; the compression adopts lossy compression or lossless compression.

Images