A fast three-dimensional inversion method for loop source transient electromagnetic using Gauss-Newton method

Abstract

The present invention discloses a three-dimensional rapid inversion method of loop source transient electromagnetic using Gauss-Newton method, the method comprising steps S1: reading in the initial model file and data file; S2: defining forward parameters and inversion parameters, S3: calculating the initial field distribution of each emission source; S4: successively performing Gauss-Newton inversion iteration process; S5: judging whether the inversion iteration termination condition is reached. The beneficial effect of the present invention is that the method utilizes the definition and setting of the model roughness weight factors in three directions to avoid the deviation caused by the three-dimensional inversion method when processing two-dimensional data. The present invention is suitable for the three-dimensional inversion of geophysical loop source transient electromagnetic detection of underground space material structure. Compared with the existing methods, it can improve the convergence speed, reduce the computer hardware requirements, avoid the deviation of the three-dimensional program inverting two-dimensional data, and can accurately obtain the three-dimensional electrical structure model of the underground to realize the detection of underground space.

Description

Technical Field

[0001] The invention relates to the field of geophysical electromagnetic detection, and in particular to a three-dimensional rapid inversion method of loop source transient electromagnetic using the Gauss-Newton method. Background Art

[0002] Transient Electromagnetic Method (TEM), also known as Time-Domain Electromagnetic Method (TDEM), is an important branch of geophysical electromagnetic detection methods. Generally, a loop source or a grounded conductor source is used to emit a primary electromagnetic field underground. After the primary field is turned off, the conductivity distribution of the underground medium is inferred by observing the pure secondary field induced by the underground medium, so as to achieve the purpose of underground medium detection.

[0003] Over the past few decades, transient electromagnetic methods have been widely used in coal and mineral resource detection, deep earth structure research, groundwater detection, and environmental engineering surveys. With the development of transient electromagnetic technology, various detection modes have evolved to adapt to different application scenarios, such as loop source transient electromagnetic, ground source transient electromagnetic, aerial transient electromagnetic, and more complex semi-aerial transient electromagnetic, well transient electromagnetic, etc. Among all these detection modes, the loop source transient electromagnetic method is currently the most widely used method, which is the focus of this invention.

[0004] At present, the interpretation of loop source transient electromagnetic data generally adopts one-dimensional inversion technology, which is mainly limited by the development of high-dimensional forward and inversion technology of electromagnetic method and the huge computing resource requirements required for high-dimensional inversion. Since the spatial influence range of the loop source is relatively limited, one-dimensional inversion can obtain a more reasonable underground structure model for cases with good stratification of underground structures or shallow detection depth. However, when the detection depth is large, the two-dimensional and three-dimensional structures of the underground structure are more prominent, or the emission source is large, the traditional one-dimensional inversion may cause an erroneous model structure.

[0005] In order to obtain more accurate transient electromagnetic data inversion results, the three-dimensional forward modeling technology of transient electromagnetic has made great progress in recent years. These forward modeling methods and various technical solutions related to them are mainly committed to improving the accuracy and efficiency of forward modeling. Although the three-dimensional forward modeling methods of transient electromagnetic are becoming more and more mature, the three-dimensional inversion technology is still in its infancy. The existing three-dimensional transient electromagnetic inversion algorithm has a slow convergence speed and a low degree of parallelism, resulting in low efficiency. At the same time, the sensitivity matrix calculation has high requirements for computer hardware, so it is difficult to apply it in large-scale underground space detection. Summary of the invention

[0006] In view of this, the present invention aims at solving the existing problems and provides a three-dimensional fast inversion method of loop source transient electromagnetic using Gauss-Newton method, which comprises the following steps:

[0007] S1. Read in the initial model file and the data file, and obtain the parameters of the initial model, the data to be inverted, the parameters of the transceiver device, and the data time channel range parameters;

[0008] The model file and data file include information such as model grid division, model conductivity distribution, transmitter location information, receiving station location information, acquisition time channel and transient electromagnetic data;

[0009] S2. defining forward modeling parameters and inversion parameters, wherein the forward modeling parameters and inversion parameters include time iteration parameters, reference model construction, initialization regularization parameters, model roughness weight factors, inversion iteration termination conditions, and parallel parameters;

[0010] S3. Calculate the initial field distribution of each emission source: Use the pseudo finite volume method based on the structured Yee grid to discretize the Maxwell equations implicitly discretized in the time domain in the space domain to obtain the discrete form of the Helmholtz equation; then process the discretized Maxwell equations at time 0 to obtain the initial field calculation equation:

[0011] (1)

[0012] in is the magnetic field distribution at the initial moment, and Corresponding to the average matrices of dielectric conductivity and magnetic permeability, respectively, is the discretized source term; is the discretized gradient operator matrix, is the discretized curl operator matrix;

[0013] According to the meshing of the initial model parameters read in step S1 and the model conductivity distribution, the matrix in equation (1) is constructed. , , and Matrix; construct the discretized source term in formula (1) according to the transmission source information of the transceiver parameters read in step S1 ; Adjust for different emission sources , by successively solving equation (1), we can obtain the initial field distribution of the emission source;

[0014] S4, according to the time iteration parameter, reference model, initial regularization parameter, model roughness weight factor and parallel parameter settings in step S2, the following Gauss-Newton inversion iteration process is performed successively: the emission source is distributed to different calculation processes, the sensitivity matrix of the implicit calculation model , calculation model update direction , calculate the model update step , build update model , perform forward simulation on the updated model and calculate the data fitting error RMS, store data files and model files;

[0015] S5. Determine whether the inversion iteration termination condition is met. If the termination condition is met, stop the inversion iteration, output the inversion process file, and release the memory usage. If the inversion iteration termination condition is not met, perform step 4 calculation again.

[0016] Furthermore, in the above step S3, the initial magnetic field is solved by directly solving the time domain Maxwell equations at time 0 using a matrix.

[0017] Furthermore, in the above step S4, the sensitivity matrix of the implicit calculation model is The process is:

[0018] The pseudo-finite volume method based on structured Yee grid is used to discretize the Maxwell equations in the space domain after implicit discretization in the time domain, and the discrete form of the Helmholtz equation is obtained:

[0019] (2)

[0020] in is the inverse of the discrete time step, for i The discrete electric field distribution at time , for i The discrete magnetic field distribution at the moment, is the average matrix of magnetic permeability, for i The distribution of emission sources at time +1. ;

[0021] The specific time of the iteration time i in formula (2) is determined by the time range of the data read in step S1;

[0022] The MUMPS direct solver is used to solve the above equation (2) at each moment. i The solution is obtained at each moment, that is, the forward response of the model is obtained; the implicit sensitivity matrix calculation includes calculating the electric field edge sensitivity matrix and vector The product of and the sensitivity matrix at the measuring point With vector The product of

[0023] and vector The product of can be expressed as

[0024] (3)

[0025] By forward modeling the above equation (2), we can get the result of equation (3). With vector After calculating and vector The product of and The sensitivity matrix is ​​completed after multiplying Implicit solution of ;

[0026] The vector involved in the above formula (3) is The model meshing and model conductivity are constructed based on the initial model parameters read in step S1; vector The transient electromagnetic data vector that needs to be inverted is read in step S1;

[0027] The initial model at the starting point of the calculation process in the above step S4 is determined by the initial model parameters read in step S1, and the initial regularization parameter at the starting point of the calculation is determined by the parameter set in step S2, and is gradually reduced during the iteration process.

[0028] Furthermore, in the above step S4, the model update direction is calculated The process is:

[0029] The Gauss-Newton method is used for inversion. k During the Gauss-Newton iteration, the following normal equation is solved to obtain the update direction of the model:

[0030] (4)

[0031] in is the first-order derivative of the inversion objective function. is a Hessian matrix, which is approximated by

[0032] (5)

[0033] in For the k The sensitivity matrix of the secondary model is calculated implicitly using the previous steps; and are the first-order covariance matrices of the data and the model respectively; the model update direction is obtained by solving equation (4) using the conjugate gradient method .

[0034] Furthermore, in the above step S4, the first-order covariance matrix of the model is Defined as

[0035] (6)

[0036] , and They are respectively the roughness weight factors of the model in three different directions. By adjusting these three weight factors, different requirements of continuity in different directions of the model can be met, thus avoiding errors caused by inverting two-dimensional data with a three-dimensional inversion program.

[0037] Furthermore, in the above step S4, the model update step length γ is calculated and the updated model is constructed. The process is:

[0038] Use unit step size, i.e., γ = 1, to test whether the updated model meets the requirements of the objective function reduction. If it does not meet the requirements, gradually reduce the step size until γ meets the requirements of the objective function reduction. The final updated model is .

[0039] Furthermore, in the above step S4, the process of updating the model to perform forward simulation and calculate the data fitting error RMS is:

[0040] Update model The forward calculation is performed using equation (2) to obtain the model response , and then use the following formula to calculate the data fitting error

[0041] (7)

[0042] in i Indicates the data sequence number. N is the total data volume.

[0043] Furthermore, in the above step S4, the implicit calculation of the sensitivity matrix is ​​based on the parallel parameters set in step S2, and all the emission sources are distributed in different processes for parallel calculation. The MUMPS direct solver is used for parallel solution in the sensitivity matrix calculation.

[0044] Furthermore, the termination condition of the inversion iteration in step S5 is (1) the data fitting error RMS reaches the convergence threshold requirement set initially; or (2) the number of inversion iterations in step 4 reaches the maximum number of iterations set initially.

[0045] The invention provides a three-dimensional rapid inversion method of loop source transient electromagnetic using Gauss-Newton method, which has the following beneficial effects:

[0046] This method combines the time-domain finite volume method with a matrix direct solver to achieve transient electromagnetic three-dimensional forward modeling, and uses the Gauss-Newton inversion method to improve the inversion convergence speed.

[0047] This method adopts an implicit sensitivity matrix calculation method and obtains the iterative search direction of the inversion model through the conjugate gradient method, avoiding the large amount of computing time and computer memory requirements caused by the explicit calculation and storage of the sensitivity matrix.

[0048] This method can effectively avoid the deviation caused by the three-dimensional inversion method when processing two-dimensional data by defining and setting the model roughness weight factors in three directions.

[0049] The present invention is particularly suitable for three-dimensional inversion of loop source transient electromagnetic, including devices such as a central loop and a fixed source loop, which can obtain an effective three-dimensional electrical structure model of the underground space and ensure the high efficiency of the calculation, thereby realizing efficient three-dimensional inversion of loop source transient electromagnetic data, meeting the three-dimensional interpretation requirements of transient electromagnetic detection data and improving the prediction accuracy of the underground resistivity model. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] Figure 1 is a method flow chart of an embodiment of the present invention;

[0051] Figure 2 It is a three-dimensional synthetic conductivity test model of the present invention; Figure 2 (a) Xy horizontal slice of the synthetic model, where the black solid line box is 500×500 m 2 Transmitting loop, the black dashed box is the horizontal boundary of the anomaly, and the black dot is the data receiving point; Figure 2 (b) The xz vertical slice of the synthetic model, the solid black frame is the boundary of the abnormal body;

[0052] Figure 3 is the RMS variation process of the three-dimensional inversion iteration process of the three-dimensional synthetic model of the present invention;

[0053] Figure 4 It is a three-dimensional inversion result model diagram of the three-dimensional synthetic model of the present invention; (ac) xy horizontal slice diagram of the inversion result model of the dBz / dt component alone; (d) xz vertical slice diagram of the inversion result model of the dBz / dt component alone; (eg) xy horizontal slice diagram of the inversion result model of the dBx / dt component alone; (h) xz vertical slice diagram of the inversion result model of the dBx / dt component alone; (ik) xy horizontal slice diagram of the inversion result model of the three components jointly; (l) xz vertical slice diagram of the inversion result model of the three components jointly;

[0054] Figure 5 It is the comparison between the observation data of the synthetic model of the present invention and the inversion result data at two measuring points;

[0055] Figure 6 It is the inversion result of the transient electromagnetic measured data of the present invention. Figure 6 (a) is the inversion result of the underground medium conductivity distribution; Figure 6 (b) in the figure is the data fitting error of each measuring point;

[0056] Figure 7 It is the inversion result of transient electromagnetic measured data of the present invention, and the comparison between the measured data and the inversion result data at two measuring points. DETAILED DESCRIPTION

[0057] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the 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 used to limit the present invention.

[0058] Please refer to Figure 1-7 The three-dimensional rapid inversion method of loop source transient electromagnetic using Gauss-Newton method of the present invention comprises the following steps:

[0059] S1. Read in the initial model file and the data file, and obtain the parameters of the initial model, the data to be inverted, the parameters of the transceiver device, and the data time channel range parameters;

[0060] The model files and data files include: model grid division, model conductivity distribution, transmitter location information, receiving station location information, acquisition time channel and transient electromagnetic data.

[0061] S2. Define forward modeling parameters and inversion parameters, wherein the forward modeling parameters and inversion parameters include time iteration parameters, reference model construction, initialization regularization parameters, model roughness weight factors, inversion iteration termination conditions, and parallel parameters.

[0062] S3. Calculate the initial field distribution of each emission source: Use the pseudo finite volume method based on the structured Yee grid to discretize the Maxwell equations implicitly discretized in the time domain in the space domain to obtain the discrete form of the Helmholtz equation; then process the discretized Maxwell equations at time 0 to obtain the initial field calculation equation:

[0063] (1)

[0064] in is the magnetic field distribution at the initial moment, and Corresponding to the average matrices of dielectric conductivity and magnetic permeability, respectively, is the discretized source term; is the discretized gradient operator matrix, is the discretized curl operator matrix.

[0065] Further, the matrix in equation (1) is constructed based on the meshing of the initial model parameters read in step S1 and the model conductivity distribution , , and Matrix; construct the discretized source term in formula (1) according to the transmission source information of the transceiver parameters read in step S1 ; Adjust for different emission sources , solve equation (1) successively to obtain the initial field distribution of the emission source;

[0066] S4, according to the time iteration parameter, reference model, initial regularization parameter, model roughness weight factor and parallel parameter settings in step S2, the following Gauss-Newton inversion iteration process is performed successively: the emission sources are distributed to different calculation processes, and the sensitivity matrix of the implicit calculation model is calculated. , calculation model update direction , calculate the model update step , build update model , perform forward simulation on the updated model and calculate the data fitting error RMS, and store the data file and model file;

[0067] In step 4 above, the sensitivity matrix The calculation is performed in an implicit way. To avoid the huge time consumption and memory requirements of explicit calculation, the implicit solution of the sensitivity matrix involves multiple forward calculations. The pseudo-finite volume method based on structured Yee grid is used to discretize the Maxwell equations in the space domain after implicit discretization in the time domain.

[0068] (2)

[0069] (3)

[0070] in and are the electric field intensity and magnetic induction intensity respectively, and are the magnetic permeability and dielectric constant in vacuum, is the medium conductivity; represents the current density of the external current source.

[0071] The implicit time step iteration method - backward Euler method is used to discretize the Maxwell equations in the time domain, and the time domain discretized forms of equations (2) and (3) are obtained to obtain the discretized form of the Helmholtz equation:

[0072] (4)

[0073] (5)

[0074] The superscript i and i +1 indicates two adjacent moments of electromagnetic field iteration, The present invention uses the finite volume method to discretize the Maxwell equations in the space domain after time domain discretization. The space domain discretization grid adopts the Yee staggered grid method, in which the electric field is defined at the edge and the magnetic field is defined at the center of the surface. The equation system after discretization is expressed by the following formula:

[0075] (6)

[0076] in is the inverse of the discrete time step, for i The discrete electric field distribution at time , for i The discrete magnetic field distribution at the moment, is the average matrix of magnetic permeability, for i The specific time of iteration time i in (2) is determined by the data time range read in step S1. The MUMPS direct solver is used to solve the above equation at the moment i By solving at each moment, the forward response of the model can be obtained.

[0077] Implicit sensitivity matrix calculation mainly involves the electric field edge sensitivity matrix and vector The product of and the sensitivity matrix at the measuring point With vector The product of

[0078] and vector The product of can be expressed as

[0079] (7)

[0080] The result of (7) can be obtained by forward calculation of (6) once. The same method can also be used to obtain With vector After calculating and vector The product of and The sensitivity matrix is ​​completed after multiplying Implicit solution of .

[0081] Furthermore, the vector involved in the above formula (7) The model meshing and model conductivity are constructed based on the initial model parameters read in step S1; vector The transient electromagnetic data vector that needs to be inverted is read in step S1;

[0082] Furthermore, the initial model at the starting point of the calculation process in the above step S4 is determined by the initial model parameters read in step S1, and the initial regularization parameter at the starting point of the calculation is determined by the parameter set in step S2, and is gradually reduced during the iteration process;

[0083] In step 4 above, the model update direction is calculated The present invention adopts the gradient optimization inversion method - Gauss-Newton method for inversion. During the kth Gauss-Newton iteration, the following normal equation is solved to obtain the update direction of the model:

[0084] (8)

[0085] in is the first-order derivative of the inversion objective function. is a Hessian matrix, which can be approximated by

[0086] (9)

[0087] in is the sensitivity matrix of the k-th model, which is calculated implicitly using the above steps. and are the first-order covariance matrices of the data and the model, respectively. The expression form is

[0088] (10)

[0089] in , and are the first-order difference matrices in three directions respectively, , and They are respectively the model roughness weight factors. By adjusting the model roughness weight factors, the smoothness setting of the inversion model in different directions can be achieved, thereby effectively avoiding the deviation caused by the inversion of two-dimensional data by the three-dimensional inversion code.

[0090] The present invention solves equation (8) by the conjugate gradient method (CG) to obtain the model update direction: .

[0091] In step 4 above, the model update step length γ is calculated and the updated model is constructed The present invention adopts the reverse tracking method, also known as the Armijo line search method, to obtain the update step size γ. First, a unit step size, i.e., γ=1, is used to test whether the updated model meets the requirement of reducing the objective function. If it does not meet the requirement, the step size is gradually reduced until γ meets the requirement of reducing the objective function. Finally, the updated model is .

[0092] In step 4 above, the updated model is forward simulated and the RMS of the data fitting error is calculated. Formula (6) is used for forward calculation to obtain the model response , and then use the following formula to calculate the data fitting error

[0093] (11)

[0094] in i Indicates the data sequence number. N is the total data volume.

[0095] In the above step 4, the present invention adopts a multi-level parallel strategy according to the parallel parameters set in step S2 to improve the inversion efficiency and meet the needs of large-scale transient electromagnetic data inversion. Figure 1 As shown, the transmitting source Tx1 to the transmitting source Tx n Assigned to computing process 1, transmitter Tx n+1 To the transmitter Tx 2n Assign it to calculation process 2, and so on, assign all the emission sources that need to be calculated to the corresponding process for calculation, and each process is independent. Then, when MUMPS solves large sparse matrices, it solves them in a multi-threaded parallel manner to maximize the utilization of CPU resources and computer memory, improve inversion efficiency, and shorten calculation time.

[0096] S5, judging whether the inversion iteration termination condition is reached, if the termination condition is reached, stopping the inversion iteration, outputting the inversion process file, and releasing the memory usage, if the inversion iteration termination condition is not reached, performing step 4 calculation again;

[0097] The termination condition of the inversion iteration is determined by the parameters set in step S2, including: (1) the RMS of the data fitting error reaches the convergence threshold requirement set initially; or (2) the number of inversion iterations in step 4 reaches the maximum number of iterations set initially.

[0098] The following is a case study of a synthetic model to demonstrate the present invention and test the application effect of the method of the present invention. The present invention designs a three-dimensional synthetic model to verify the effectiveness of the developed algorithm. Figure 2 The size shown is 400×400×300 m 3The 1 Ωm high-conductivity anomaly is embedded in a 500 Ωm uniform half-space, where the top interface of the anomaly is 200 m above the ground and the size of the emission source is 500×500 m. 2 The horizontal component (dBx / dt and dBy / dt) and vertical component (dBz / dt) of the transient electromagnetic response are used to detect the target. Figure 2 As shown in (a), a total of 169 data acquisition points are arranged on the ground, and each acquisition point collects horizontal and vertical components. The data acquisition time channel range is t=31.6 μs to 10 ms, with a total of 26 time channels with equal logarithmic intervals. 3% Gaussian noise is added to the forward modeling data to simulate the measured data.

[0099] The above models and data are processed step by step according to the method disclosed in the present invention:

[0100] S1. Read the prepared initial model file and data file, and obtain the parameters of the initial model, the data to be inverted, the parameters of the transceiver, and the parameters of the data time channel range. The main information contained in the initial model file is: the grid is divided into 50×50×29, the size of the core area grid in three directions is 25 m, the boundary grid is extended to 13 km, and the initial model is a uniform half space with a conductivity of 500 Ωm. The main information contained in the data file is: the center position of the emission source is (0,0,0,) m, and the size of the emission loop is 500×500 m 2 , the specific locations of 169 measuring points, the time settings of 26 data time channels, and the measurement data of each measuring point in each time channel, totaling 4394 inversion data.

[0101] S2. Define forward modeling parameters and inversion parameters. Time iteration parameters include: (10 -6 , 10 -2 , 10 -6 )s, doubling factor is 3.0, number of equal steps is 30; the reference model is a uniform half-space of 500 Ωm; the initial regularization parameter is 10 -2 ; The model roughness weight factor is 0.5 in all three directions; the inversion iteration termination conditions include: RMS threshold is 1.0, the maximum number of iterations is 25; the parallel parameters include: the number of processes is 1, and the maximum number of threads used by a single process is 8.

[0102] S3. Calculate the initial field distribution of each emission source. For this test case, there is only one emission source, and only 500×500 m 2 The initial field distribution is calculated according to the settings in the data file and the model file by using the transient electromagnetic initial field calculation method developed by the present invention.

[0103] S4, according to the data files, model files and forward and inversion parameter settings read in steps S1 and S2, and using the initial field distribution calculated in step S3, perform Gauss-Newton iterative inversion, that is, gradually calculate the sensitivity matrix of the implicit model , model update direction , model update step length γ, update model Perform calculations. Perform forward simulation on the updated model and calculate the data fitting error RMS, and store the data file and model file.

[0104] S5. This case tests the applicability of the present invention to the vertical component data and horizontal component data of transient electromagnetics. When the vertical component data dBz / dt is used for testing, after 7 Gauss-Newton inversion iterations of the S4 step process, the data fitting error RMS can be reduced to 0.85, reaching the threshold of RMS=1, stopping the inversion iteration and releasing the memory usage; when the horizontal component dBx / dt is used for testing, after 25 Gauss-Newton inversion iterations of the S4 step process, the data fitting error RMS is reduced to 1.74, which does not reach the threshold of RMS=1, reaching the maximum 30 iteration termination conditions, stopping the inversion iteration and releasing the memory usage; when the horizontal component dBx / dt, dBy / dt and vertical component dBz / dt are used for testing at the same time, after 30 Gauss-Newton inversion iterations of the S4 step process, the data fitting error RMS is reduced to 1.45, which does not reach the threshold of RMS=1, reaching the maximum 25 iteration termination conditions, stopping the inversion iteration and releasing the memory usage.

[0105] The changes of RMS during the inversion of three different types of transient electromagnetic data are shown in the figure. Figure 3 The inversion test results are shown in the attached Figure 4 As shown in the figure, from the effects of horizontal and vertical slices of the inversion result model, it can be seen that each model has been effectively reconstructed. Figure 5 It is shown that the observed data at the two measuring points fit well with the inversion result data, which comprehensively illustrates the effectiveness of the present invention for the three-dimensional inversion of the transient electromagnetic synthetic model.

[0106] The following is a practical application case of transient electromagnetic detection to demonstrate the present invention and test the application effect of the method of the present invention. At this time, the transient electromagnetic detection is located on the lake surface, and the detection purpose is to study the underground structure and underground volcanic system. This transient electromagnetic detection has a total of one survey line, with a total of 24 data collection points in the range of x = 100m to 450m, and its location is nearby Figure 6 The black triangle is used in the figure. The size of the transmitting loop is a 18m square loop, and the maximum effective detection time range is 10μs to 6.1ms. The measured component is the vertical component of the transient electromagnetic signal dBz / dt.

[0107] The above measurement data are processed step by step according to the method disclosed in the present invention:

[0108] S1. Read the prepared initial model file and data file, and obtain the parameters of the initial model, the data to be inverted, the parameters of the transceiver, and the parameters of the data time channel range. The main information contained in the initial model file is: the grid is divided into 105×35×39, the size of the core area grid in three directions is 4.5 m, the boundary grid is extended to 6 km, and the initial model is a uniform half-space with a conductivity of 300 Ωm. The main information contained in the data file is: the center position of the emission source is 100 to 450m, and the size of the emission loop is 18×18 m 2 , the specific locations of the 24 emission sources, the time settings of each data time channel, and the measurement data of each measuring point in each time channel, a total of 521 inversion data.

[0109] S2. Define forward modeling parameters and inversion parameters. Time iteration parameters include: (10 -6 , 10 -2 , 10 -6 )s, doubling factor is 3.0, number of equal steps is 30; the reference model is a uniform half-space of 300 Ωm; the initial regularization parameter is 10 -2 ; Model roughness weight factor , as well as , to maintain the continuity in the y direction. 5; The inversion iteration termination conditions include: RMS threshold is 1.0, the maximum number of iterations is 25; The parallel parameters include: the number of processes is 6, and the maximum number of threads used by a single process is 4.

[0110] S3. Calculate the initial field distribution of each emission source. For this test case, the number of emission sources is only 24, and the initial field distribution needs to be calculated 24 times. According to the settings in the data file and the model file, the initial field is calculated using the transient electromagnetic initial field calculation method developed by the present invention.

[0111] S4, according to the data files, model files and forward and inversion parameter settings read in steps S1 and S2, and using the initial field distribution calculated in step S3, perform Gauss-Newton iterative inversion, that is, gradually calculate the sensitivity matrix of the implicit model , model update direction , model update step length γ, update model Perform calculations. Perform forward simulation on the updated model and calculate the data fitting error RMS, and store the data file and model file.

[0112] S5. When the vertical component data dBz / dt is used for testing, after 10 Gauss-Newton inversion iterations in step S4, the data fitting error RMS can be reduced to 0.97, reaching the threshold of RMS=1, stopping the inversion iteration and releasing the memory usage.

[0113] The inversion test results of the transient electromagnetic measured data in this embodiment are shown in the attached Figure 6 As shown, Figure 6 (a) is the inversion result model, which successfully obtains the underground structure and medium distribution of the lake surface in this area. Figure 6 (b) is the data fitting error of each measuring point. Figure 7 The correspondence between the measured data and the inversion result data of the two selected measuring points is shown. From the inversion result model and data fitting, the effectiveness of the present invention for inversion of measured transient electromagnetic data is explained.

[0114] It should be understood that the implementation of the present invention is described above in conjunction with the accompanying drawings, but the present invention is not limited to the above-mentioned specific implementation methods. The above-mentioned specific implementation methods are merely illustrative and not restrictive. For ordinary technicians in this field, improvements or changes can be made based on the above description, and all these improvements and changes should fall within the scope of protection of the claims attached to the present invention.

Claims

1. A three-dimensional fast inversion method of loop source transient electromagnetic using Gauss-Newton method, characterized by: The steps include: S1. Read in the initial model file and data file, and obtain the initial model parameters, data to be inverted, transceiver parameters, and data time channel range parameters; The initial model file and data file include information such as model grid division, model conductivity distribution, transmitter location information, receiving station location information, acquisition time channel and transient electromagnetic data; S2. defining forward modeling parameters and inversion parameters, wherein the forward modeling parameters and inversion parameters include time iteration parameters, reference model, initial regularization parameters, model roughness weight factor, inversion iteration termination condition, and parallel parameters; S3. Calculate the initial field distribution of each emission source: Use the pseudo finite volume method based on the structured Yee grid to discretize the Maxwell equations implicitly discretized in the time domain in the space domain to obtain the discrete form of the Helmholtz equation; then process the discretized Maxwell equations at time 0 to obtain the initial field calculation equation: Among them B 0 is the magnetic field distribution at the initial moment, M eσ and M fμ The average matrices corresponding to the dielectric conductivity and magnetic permeability, S r is the discretized source term; Grad is the discretized gradient operator matrix, Curl is the discretized curl operator matrix, and T represents the transpose of the matrix; The matrix M in equation (1) is constructed based on the meshing of the model and the conductivity distribution of the model read in step S1. eσ 、M fμ , Curl and Grad matrices; construct the discrete source term S in formula (1) according to the transmission source information of the transceiver parameters read in step S1 r ; Adjust S for different emission sources r , by successively solving equation (1), we can obtain the initial field distribution of the emission source; S4, according to the time iteration parameters, reference model, initial regularization parameters, model roughness weight factor and parallel parameters set in step S2, the following Gauss-Newton inversion iteration process is performed successively: the emission sources are distributed to different calculation processes, the sensitivity matrix J of the implicit calculation model, the calculation model update direction δm, the calculation model update step length γ, and the construction of the updated model m k+1 =m k +γδm, where k represents the number of Gauss-Newton iterations, forward simulation is performed on the updated model, the RMS of the data fitting error is calculated, and the data file and model file are stored; S5, judging whether the inversion iteration termination condition is reached, if the termination condition is reached, stopping the inversion iteration, outputting the inversion process file, and releasing the memory usage, if the inversion iteration termination condition is not reached, performing step S4 calculation again; In the above step S4, the process of implicitly calculating the sensitivity matrix J of the model is: The pseudo-finite volume method based on structured Yee grid is used to discretize the Maxwell equations in the space domain after implicit discretization in the time domain, and the discrete form of the Helmholtz equation is obtained: Where β is the inverse of the discrete time step, E i is the discrete electric field distribution at time i, B i is the discrete magnetic field distribution at time i, M eε is the average matrix of dielectric constants, is the emission source distribution at time i+1; The specific time of the iteration time i in formula (2) is determined by the time range of the data read in step S1; The MUMPS direct solver is used to solve the above equation (2) at the i-th moment, that is, to obtain the forward response of the model; the implicit sensitivity matrix calculation includes calculating the electric field edge sensitivity matrix J E and vector v m The product of and the sensitivity matrix at the measuring point With vector v d The product of J E and vector v m The product of can be expressed as Where n is the electric field edge sensitivity matrix J E The number of columns is obtained by forward calculation of the above formula (2) to obtain the result of formula (3). The same method is used to obtain With vector v d After calculating the product of J E and vector v m The product of With v d After multiplying by , the implicit solution of the sensitivity matrix J is completed; The vector v involved in the above formula (3) is m The vector v is constructed based on the meshing of the model read in step S1 and the model conductivity distribution; d The transient electromagnetic data vector of the data to be inverted read in step S1; The initial model at the starting point of the calculation process in the above step S4 is determined by the initial model parameters read in step S1, and the initial regularization parameter at the starting point of the calculation is determined by the parameter set in step S2, and is gradually reduced during the iteration process.

2. The method for three-dimensional rapid inversion of loop source transient electromagnetic using Gauss-Newton method according to claim 1, characterized in that: In the above step S3, the initial magnetic field is solved by directly solving the time domain Maxwell equations at time 0 using a matrix.

3. The method for three-dimensional rapid inversion of loop source transient electromagnetic using Gauss-Newton method according to claim 1, characterized in that: In the above step S4, the process of calculating the model update direction δm is: The Gauss-Newton method is used for inversion. During the kth Gauss-Newton iteration, the following normal equation is solved to obtain the update direction of the model: H k δm k =-g k (4) where g k is the first-order derivative of the inversion objective function; H k is a Hessian matrix, which is approximated by Among them J k is the sensitivity matrix of the k-th model, calculated using the implicit method of the previous steps, is the model parameter; C d and C m are the first-order covariance matrices of the data and the model respectively; the model update direction δm is obtained by solving equation (4) using the conjugate gradient method.

4. The method for three-dimensional rapid inversion of loop source transient electromagnetic using Gauss-Newton method according to claim 3, characterized in that: In the above step S4, the first-order covariance matrix C of the model is m Defined as D x , D y and D z are the first-order difference matrices in the x, y and z directions respectively; α x , α y and α z They are the model roughness weight factors in three different directions, namely x, y and z. By adjusting these three weight factors, different continuity requirements of the model in different directions can be met, avoiding errors caused by inverting two-dimensional data using a three-dimensional inversion program.

5. The method for three-dimensional rapid inversion of loop source transient electromagnetic using Gauss-Newton method according to claim 1, characterized in that: In the above step S4, the model update step length γ is calculated and the updated model m is constructed. k+1 =m k The process of +γδm is: Use unit step size, i.e., γ = 1, to test whether the updated model meets the requirement of reducing the objective function. If it does not meet the requirement, gradually reduce the step size until γ meets the requirement of reducing the objective function. Finally, the updated model is m k+1 =m k +γδm.

6. The method for three-dimensional rapid inversion of loop source transient electromagnetic using Gauss-Newton method according to claim 1, characterized in that: In the above step S4, the implicit calculation of the sensitivity matrix is ​​based on the parallel parameters set in step S2, all the emission sources are distributed in different processes for parallel calculation, and the sensitivity matrix calculation is solved in parallel using the MUMPS direct solver.

7. The method for three-dimensional rapid inversion of loop source transient electromagnetic using Gauss-Newton method according to claim 1, characterized in that: The termination condition of the inversion iteration in step S5 is: (1) the RMS of the data fitting error reaches the convergence threshold requirement set initially; or (2) the number of inversion iterations in step 4 reaches the maximum number of iterations set initially.

Images