A three-dimensional time-frequency domain electromagnetic joint inversion method based on a finite element algorithm

Abstract

The present application relates to the technical field of geophysical electromagnetic data exploration, in particular to a three-dimensional time-frequency domain electromagnetic joint inversion method based on finite element algorithm. The method obtains time domain and frequency domain electromagnetic observation data of the region to be measured respectively; the region to be measured is profiled to obtain a non-structural tetrahedral mesh to construct a time-frequency domain inversion initial model; the time-frequency domain inversion initial model is forward calculated to obtain time domain and frequency domain electromagnetic response data; time domain electromagnetic data residuals and frequency domain electromagnetic data residuals are obtained, and their corresponding adaptive weights are obtained respectively to establish a joint inversion objective function; the joint inversion objective function is iteratively solved to calculate the update vector of the inversion initial model; the inversion initial model is iteratively updated until the set condition is reached, and the time-frequency domain electromagnetic joint inversion is completed. The present application effectively improves the accuracy of three-dimensional inversion through joint inversion, thereby enhancing the practical effect of three-dimensional inversion.

Description

Technical Field

[0001] This invention relates to the field of geophysical electromagnetic data exploration technology, specifically to a three-dimensional time-frequency domain electromagnetic joint inversion method based on the finite element algorithm. Background Technology

[0002] Based on domestic demand for mineral resources and the efforts of Chinese scholars, the exploration of mineral resources is becoming increasingly sophisticated. The main methods employed include time-domain electromagnetic methods and frequency-domain electromagnetic methods. The time-domain electromagnetic method uses coils or grounded electrodes to detect underground anomalies by transmitting a pulsed magnetic field during power outages and observing the secondary eddy current field response at different time intervals during the intermittent period, based on the principle of electromagnetic induction. The frequency-domain electromagnetic method establishes alternating electric and magnetic fields of different frequencies, and based on their distribution patterns and the differences in resistivity and permeability of the anomaly, it uses the skin effect principle to obtain electrical information of geological bodies at different depths through frequency or distance variations.

[0003] However, conventional time-domain and frequency-domain methods each have their own shortcomings. The time-domain transient electromagnetic method is affected by its propagation mechanism and cannot obtain effective information about deep strata in actual exploration, thus having the limitation of relatively small exploration depth. On the other hand, for the frequency-domain controlled-source electromagnetic method, a large detection depth can be achieved by selecting an appropriate offset distance and a lower frequency band. However, the induced current density of low-frequency signals in shallow layers is extremely small, making it insensitive to differences in shallow electrical properties and difficult to maintain the resolution of shallow layer detection. Furthermore, in actual exploration, shallow layer information is often interfered with by high-frequency noise, resulting in the loss of effective high-frequency data. Summary of the Invention

[0004] To address the problems existing in the prior art, this invention provides a three-dimensional time-frequency domain electromagnetic joint inversion method based on the finite element method. This method involves acquiring time-domain and frequency-domain electromagnetic observation data of the region to be measured; constructing an initial time-frequency domain inversion model by subdividing the region into an unstructured tetrahedral mesh; performing forward modeling on the initial time-frequency domain inversion model to obtain time-domain and frequency-domain electromagnetic response data; acquiring the time-domain and frequency-domain electromagnetic data residuals and their corresponding adaptive weights to establish a joint inversion objective function; iteratively solving the joint inversion objective function to calculate the update vector of the initial inversion model; and iteratively updating the initial inversion model until the set conditions are met, thus completing the time-frequency domain electromagnetic joint inversion. This invention effectively improves the accuracy of three-dimensional inversion through joint inversion, thereby enhancing the practical application of three-dimensional inversion.

[0005] This invention adopts the following technical solution: a three-dimensional time-frequency domain electromagnetic joint inversion method based on the finite element algorithm, comprising:

[0006] Time-domain electromagnetic observation data and frequency-domain electromagnetic observation data of the area to be measured are acquired respectively.

[0007] The region to be tested is divided into unstructured tetrahedral meshes, and an initial time-frequency domain inversion model is constructed based on the unstructured tetrahedral meshes.

[0008] The time-frequency domain inversion initial model is subjected to forward modeling calculation using a forward modeling algorithm to obtain time-domain electromagnetic response data and frequency-domain electromagnetic response data.

[0009] The residuals of time-domain electromagnetic data are obtained from time-domain electromagnetic observation data and time-domain electromagnetic response data, and the residuals of frequency-domain electromagnetic data are obtained from frequency-domain electromagnetic observation data and frequency-domain electromagnetic response data.

[0010] Adaptive weights are obtained for the residuals of electromagnetic data in the time domain and the frequency domain, respectively.

[0011] A joint inversion objective function is established based on the time-domain electromagnetic data residuals and the frequency-domain electromagnetic data residuals and their corresponding adaptive weights.

[0012] The joint inversion objective function is solved iteratively based on the quasi-Newton algorithm, and the update vector of the initial inversion model is calculated.

[0013] The initial inversion model is iteratively updated until the set conditions are met, thus completing the joint electromagnetic inversion in the time and frequency domains.

[0014] Furthermore, time-domain electromagnetic observation data and frequency-domain electromagnetic observation data of the area to be measured are obtained separately, specifically as follows:

[0015] In the area to be tested, array transceivers corresponding to the time domain and frequency domain are respectively deployed. The distance between the array transceiver corresponding to the time domain and the area to be tested is a first predetermined distance, and the distance between the array transceiver corresponding to the frequency domain and the area to be tested is a second predetermined distance. The first predetermined distance is less than the second predetermined distance.

[0016] The corresponding time-domain electromagnetic observation data and frequency-domain electromagnetic observation data are acquired in the region under test through array transceivers corresponding to the time domain and frequency domain, respectively.

[0017] Furthermore, the joint inversion objective function is expressed as:

[0018] ;

[0019] In the formula, For the joint inversion objective function, The time-domain data covariance matrix, For the frequency domain data covariance matrix, To constrain the covariance matrix of the model, As a regularization factor, For prediction models, This is the initial model for time-frequency domain inversion. For time-domain electromagnetic observation data, For frequency domain electromagnetic observation data, For time-domain electromagnetic response data, For frequency domain electromagnetic response data, For adaptive weights of time-domain electromagnetic data residuals, For adaptive weights of frequency domain electromagnetic data residuals, express Norm.

[0020] Furthermore, adaptive weights are obtained for the time-domain electromagnetic data residuals and the frequency-domain electromagnetic data residuals, respectively, as follows:

[0021] The expression for the adaptive weights of the time-domain electromagnetic data residuals is:

[0022] ;

[0023] The expression for the adaptive weight of the frequency domain electromagnetic data residual is:

[0024] ;

[0025] in, , representing the residual of electromagnetic data in the time domain; , representing the residual of electromagnetic data in the frequency domain.

[0026] Furthermore, the joint inversion objective function is solved iteratively based on a quasi-Newton algorithm, including:

[0027] The objective function is expanded using Taylor and its minimum value is taken to obtain the three-dimensional joint inversion equation;

[0028] Calculate the gradient of the objective function by taking partial derivatives of the three-dimensional joint inversion equation;

[0029] The quasi-Newton algorithm is used to perform inversion calculations based on the gradient of the objective function to obtain the update vector of the initial inversion model.

[0030] Furthermore, the three-dimensional joint inversion equation is expressed as:

[0031] ;

[0032] in, Denotes the approximate positive definite matrix after the k-th iteration. To update the vector, This is the time-domain sensitivity matrix. This is the frequency domain sensitivity matrix. The time-domain data covariance matrix, This is the transpose of the covariance matrix of the time-domain data. The residual between the observed data and the response data, For the frequency domain data covariance matrix, This is the transpose of the frequency domain data covariance matrix. The model constraint variance matrix, This is the transpose of the model constraint variance matrix. For prediction models, This is the initial model for time-frequency domain inversion.

[0033] Furthermore, the iterative solution process for the joint inversion objective function also includes: calculating the sensitivity matrix in the three-dimensional joint inversion equation using the adjoint forward modeling method, specifically:

[0034] The sensitivity matrix in the three-dimensional joint inversion equation includes: the time-domain sensitivity matrix and the frequency-domain sensitivity matrix;

[0035] Time-domain sensitivity matrix and frequency domain sensitivity matrix They are defined as follows:

[0036] ;

[0037] ;

[0038] in, This represents a time-domain prediction model. This represents a frequency domain prediction model. and The three-dimensional electromagnetic numerical simulation responses in the time domain and frequency domain, respectively, can be defined as follows:

[0039] ;

[0040] ;

[0041] in, Let be the interpolation basis function for the forward solution vector of the i-th iteration and j-th measurement point of the time-domain electromagnetic response. Let be the interpolation basis function for the forward solution vector of the frequency domain electromagnetic response at the j-th measurement point in the i-th iteration. This represents the time-domain forward electromagnetic response solution vector at the j-th measurement point in the i-th iteration. This represents the frequency domain forward electromagnetic response solution vector at the j-th measurement point in the i-th iteration;

[0042] The linear equations for the three-dimensional forward numerical simulation in the time and frequency domains are expressed as follows:

[0043] ;

[0044] ;

[0045] in, This is the time-domain sensitivity matrix. The forward electromagnetic response solution vector in the time domain. The time-domain coefficient matrix, This is the frequency domain sensitivity matrix. The forward electromagnetic response solution vector in the frequency domain. This is the frequency domain coefficient matrix;

[0046] Combining the above equations, the sensitivity matrix can be expressed as:

[0047] ;

[0048] ;

[0049] in, This represents the time-domain sensitivity matrix of the j-th measurement point in the i-th iteration. This represents the frequency domain sensitivity matrix of the j-th measurement point in the i-th iteration. Let represent the prediction model for the i-th iteration;

[0050] By taking the partial derivatives of the model parameters at both ends of the linear equations in the three-dimensional forward numerical simulation in the time and frequency domains, and combining this with the expression for the sensitivity matrix, we obtain:

[0051] ;

[0052] ;

[0053] in, Let be the interpolation basis function for the forward solution vector, defined and This is achieved by solving the following adjoint forward problem:

[0054] ;

[0055] ;

[0056] The sensitivity matrix in the three-dimensional joint inversion equation is obtained as follows:

[0057] ;

[0058] ;

[0059] in, , .

[0060] Furthermore, the setting conditions include:

[0061] The root mean square error (RMS) between the predicted electromagnetic response and the measured electromagnetic response is less than the error threshold.

[0062] The joint inversion objective function converged;

[0063] The joint inversion objective function is solved iteratively to reach the preset maximum number of iterations.

[0064] Furthermore, the initial inversion model is updated in each iteration, including:

[0065] For the joint inversion objective function at the th Prediction model at the next iteration Taking a second approximation and then calculating the partial derivative, we obtain the gradient as follows:

[0066] ;

[0067] in, This indicates the objective function in the prediction model. The approximate value of the gradient at that point. Is the objective function at the th The gradient of the prediction model at the next iteration For Hessian matrix, This is the model update vector at the k-th iteration;

[0068] Using the approximate positive definite matrix after the k-th iteration Replacing the Hessian matrix The update rule for the quasi-Newton algorithm is obtained as follows:

[0069] ;

[0070] in, Given a rank-2 matrix, the update vector for the initial model can be derived by defining its elements:

[0071] ;

[0072] in, Indicates the first The approximate positive definite matrix at the nth iteration This is the model update vector at the k-th iteration. This is the transpose of the model update vector at the k-th iteration. This represents the gradient difference of the objective function at the k-th iteration. This is the transpose of the gradient difference of the objective function at the k-th iteration.

[0073] The beneficial effects of this invention are as follows: In the three-dimensional time-frequency domain joint inversion method based on the finite element algorithm provided by this invention, the electromagnetic forward response is numerically simulated based on the unstructured finite element method. This algorithm can accurately simulate irregular physical property interfaces such as undulating surfaces. For complex geological structure models, this method can still obtain high numerical accuracy. The three-dimensional joint inversion in the time-frequency domain is based on regularization theory to construct the joint inversion objective function, and the inversion equation is derived based on the quasi-Newton numerical optimization algorithm to ensure the stability of the inversion iteration process. This algorithm can fit the data well, making the three-dimensional inversion practical. This invention also provides the idea and feasible algorithm of combining time-domain electromagnetic data and frequency-domain electromagnetic data for inversion. Compared with the separate inversion in the time domain and frequency domain, the joint inversion improves the accuracy of the three-dimensional inversion, obtains the inversion solution with better fitting, and can better enhance the practical effect of the three-dimensional inversion. Attached Figure Description

[0074] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0075] Figure 1 This is a schematic diagram of a three-dimensional time-frequency domain electromagnetic joint inversion method based on the finite element algorithm according to an embodiment of the present invention;

[0076] Figure 2 This is a schematic diagram of an inversion model according to an embodiment of the present invention;

[0077] Figure 3 This is a schematic diagram illustrating the inversion result of a three-dimensional time-frequency domain electromagnetic joint inversion method according to an embodiment of the present invention; wherein, Figure 3 Part a) is a model of a real anomaly under undulating terrain. Figure 3 Part b) shows the time-domain electromagnetic data inversion results. Figure 3 Part c) contains the frequency domain electromagnetic data inversion results. Figure 3 Part d) shows the results of the joint inversion using time-frequency domain electromagnetic data in this invention;

[0078] Figure 4 This is a schematic diagram illustrating the inversion parameters in one of three inversion processes according to an embodiment of the present invention. Detailed Implementation

[0079] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0080] A schematic diagram of a three-dimensional time-frequency domain electromagnetic joint inversion method based on the finite element algorithm according to an embodiment of the present invention is shown below. Figure 1 As shown, it includes:

[0081] Time-domain electromagnetic observation data and frequency-domain electromagnetic observation data of the area to be measured are acquired respectively.

[0082] In this embodiment of the invention, array transceivers corresponding to the time domain and frequency domain are first deployed in the area to be tested. The distance between the array transceiver corresponding to the time domain and the area to be tested is a first predetermined distance, and the distance between the array transceiver corresponding to the frequency domain and the area to be tested is a second predetermined distance. The first predetermined distance is less than the second predetermined distance. Corresponding time domain electromagnetic observation data and frequency domain electromagnetic observation data are acquired in the area to be tested through the array transceivers corresponding to the time domain and frequency domain, respectively.

[0083] Specifically, such as Figure 2 The diagram shows an inversion model in an embodiment of the present invention. The area to be measured is set as an undulating terrain with an irregular low-resistivity vein beneath it. Time-frequency domain array transceivers are deployed in the area to be measured. Long conductor sources in the near-measurement area are used to measure electromagnetic observation data in the time domain, and long conductor sources in the far-measurement area are used to measure electromagnetic observation data in the frequency domain. The offset distance of the time-domain transmitter, i.e., the first set distance between the array transceiver corresponding to the time domain and the area to be measured in this embodiment, can be set to 2 km. The offset distance of the frequency-domain transmitter, i.e., the second set distance between the array transceiver corresponding to the frequency domain and the area to be measured in this embodiment, can be set to 8 km. To make the transmitter conform to the terrain, the long conductor source is divided into multiple short conductors in the actual modeling process. Each segment can be approximated as an electric dipole. A total of 51 measuring lines are deployed in the area to be measured, and 31 measuring points are deployed on each measuring line. The distance between the measuring points and the distance between the measuring lines are both 100 m. The resistivity of the half-space in the calculation area is assumed to be 100 Ω·m, and the resistivity of the plate-like vein is assumed to be 1 Ω·m.

[0084] The region to be tested is divided into unstructured tetrahedral meshes, and an initial time-frequency domain inversion model is constructed based on the unstructured tetrahedral meshes.

[0085] Traditional local meshing techniques mostly use structured cells to discretize space. Although structured meshes can obtain sensitivity information, they have poor fitting accuracy for irregular physical property interfaces. In contrast, unstructured local meshes based on tetrahedral cells can flexibly and accurately discretize arbitrarily complex models. At the same time, time-frequency domain electromagnetic three-dimensional inversion requires the use of unstructured tetrahedral meshes to perform dense meshing of local areas, including the time-frequency domain transmitter and the measurement area, according to the location and influence range of the transceiver device. This is the core part, which improves the accuracy of the sensitivity matrix and electromagnetic response parameters required for computational inversion.

[0086] In this embodiment of the invention, the finite element unstructured tetrahedral partitioning method is used to partition the survey area into free tetrahedral elements, thereby realizing the spatial discretization of the model and the initial model construction. For complex model spatial discretization with arbitrary undulating terrain, the tetrahedral partitioning mesh structure file can be constructed based on the surface elevation information of the survey area and the known underground geological interfaces, thereby generating a model mesh that can be used for finite element algorithms.

[0087] Specifically, the model mesh is first constructed by expanding outward from the outer surface of the core region to the extended region, generating a computational region with sparse mesh between the extended region and the core region. This completes the model mesh construction. Furthermore, to facilitate the addition of the first type of Dirichlet boundary condition, the outer boundary of the extended region is defined at a sufficiently far distance from the transceiver device. This unstructured mesh generation scheme ensures high numerical simulation accuracy in the core region while reducing the required mesh size, effectively improving the speed of the inversion algorithm. In this embodiment, the discrete tetrahedral mesh is established using the mesh generation function of the commercial software Comsol. After generating the mesh using Comsol, the mesh file is exported. Subsequent mesh file conversion processes are all completed by a self-developed program module written in Fortran and Maitlab languages.

[0088] The time-frequency domain inversion initial model is subjected to forward modeling calculation using a forward modeling algorithm to obtain time-domain electromagnetic response data and frequency-domain electromagnetic response data.

[0089] In this embodiment of the invention, the electromagnetic forward response data is obtained by numerical simulation calculation based on the unstructured finite element method. The specific implementation method can refer to any of the existing technologies. This invention does not focus too much on the forward calculation process.

[0090] The residuals of time-domain electromagnetic data are obtained from time-domain electromagnetic observation data and time-domain electromagnetic response data, and the residuals of frequency-domain electromagnetic data are obtained from frequency-domain electromagnetic observation data and frequency-domain electromagnetic response data.

[0091] In this embodiment of the invention, the difference between time-domain electromagnetic observation data and time-domain electromagnetic response data obtained by forward modeling is calculated, and this difference is multiplied by the covariance matrix of the time-domain data. The product is then calculated... The norm value is used as the residual of the time-domain electromagnetic data; similarly, the difference between the frequency-domain electromagnetic observation data and the frequency-domain electromagnetic response data obtained by forward modeling is multiplied by the covariance matrix of the frequency-domain data, and the product is calculated. The value obtained from the norm can be used as the residual of the electromagnetic data in the frequency domain.

[0092] Adaptive weights are obtained for the residuals of electromagnetic data in the time domain and the frequency domain, respectively.

[0093] In this embodiment of the invention, the expression for the adaptive weight of the time-domain electromagnetic data residual is:

[0094] ;

[0095] The expression for the adaptive weighting of the frequency domain electromagnetic data residuals is:

[0096] ;

[0097] in, , representing the residual of electromagnetic data in the time domain; , representing the residual of electromagnetic data in the frequency domain.

[0098] In this embodiment of the invention, when the change of a single data residual term in the time domain or frequency domain is less than 5% of the previous change, or when the search step size is unstable, the weighting factor will be adaptively adjusted based on the current time-frequency domain residual according to the above formula, thereby realizing the adaptive weighting of the time-domain electromagnetic data residual and the adaptive weighting of the frequency-domain electromagnetic data residual.

[0099] A joint inversion objective function is established based on the time-domain electromagnetic data residuals and the frequency-domain electromagnetic data residuals and their corresponding adaptive weights.

[0100] Based on regularization theory, the objective function for the three-dimensional inversion of conventional electromagnetic data can be defined as:

[0101] ;

[0102] in, For data fitting terms, For model constraint terms, and These are observation data and response data, respectively. and These are the prediction model and the initial model for time-frequency domain inversion, respectively. and Let be the constraint covariance matrices of the data and the model, respectively. In this embodiment of the invention, to achieve three-dimensional joint inversion of time-frequency domain electromagnetic data, the observation data is divided into time-domain observation data and frequency-domain observation data, and the response data is divided into time-domain response data and frequency-domain response data. The joint inversion objective function can then be obtained from the L2 expansion, expressed as:

[0103] ;

[0104] In the formula, For the joint inversion objective function, The time-domain data covariance matrix, For the frequency domain data covariance matrix, To constrain the covariance matrix of the model, As a regularization factor, For prediction models, This is the initial model for time-frequency domain inversion. For time-domain electromagnetic observation data, For frequency domain electromagnetic observation data, For time-domain electromagnetic response data, For frequency domain electromagnetic response data, For adaptive weights of time-domain electromagnetic data residuals, For adaptive weights of frequency domain electromagnetic data residuals, express Norm;

[0105] The joint inversion objective function is solved iteratively based on the quasi-Newton algorithm, and the update vector of the initial inversion model is calculated.

[0106] In this embodiment of the invention, the objective function is first expanded using Taylor and its minimum value is taken to obtain the three-dimensional joint inversion equation, which is expressed as:

[0107] ;

[0108] in, Denotes the approximate positive definite matrix after the k-th iteration. To update the vector, This is the time-domain sensitivity matrix. This is the frequency domain sensitivity matrix. The time-domain data covariance matrix, This is the transpose of the covariance matrix of the time-domain data. The residual between the observed data and the response data, For the frequency domain data covariance matrix, This is the transpose of the frequency domain data covariance matrix. The model constraint variance matrix, This is the transpose of the model constraint variance matrix. For prediction models, This is the initial model for time-frequency domain inversion.

[0109] Calculate the gradient of the objective function by taking partial derivatives of the three-dimensional joint inversion equation:

[0110] ;

[0111] in, The residual between the observed data and the response data, and Represents the residual of electromagnetic data in the time domain. Represents the residual of electromagnetic data in the frequency domain. This is the time-domain sensitivity matrix. This is the frequency domain sensitivity matrix. The time-domain data covariance matrix, This is the transpose of the covariance matrix of the time-domain data. For the frequency domain data covariance matrix, This is the transpose of the frequency domain data covariance matrix. The model constraint variance matrix, This is the transpose of the model constraint variance matrix. For prediction models, This is the initial model for time-frequency domain inversion.

[0112] The quasi-Newton algorithm is used to perform inversion calculations based on the gradient of the objective function to obtain the update vector of the initial inversion model.

[0113] In this embodiment of the invention, the obtained three-dimensional joint inversion equation contains a Jacobian matrix term. Therefore, in performing the inversion calculation, this embodiment of the invention needs to use the adjoint forward modeling method to calculate the sensitivity matrix in the three-dimensional joint inversion equation, specifically:

[0114] The sensitivity matrix in the three-dimensional joint inversion equation includes: the time-domain sensitivity matrix and the frequency-domain sensitivity matrix;

[0115] Time-domain sensitivity matrix and frequency domain sensitivity matrix They are defined as follows:

[0116] ;

[0117] ;

[0118] in, This represents a time-domain prediction model. This represents a frequency domain prediction model. and The three-dimensional electromagnetic numerical simulation responses in the time domain and frequency domain, respectively, can be defined as follows:

[0119] ;

[0120] ;

[0121] in, Let be the interpolation basis function for the forward solution vector of the i-th iteration and j-th measurement point of the time-domain electromagnetic response. Let be the interpolation basis function for the forward solution vector of the frequency domain electromagnetic response at the j-th measurement point in the i-th iteration. This represents the time-domain forward electromagnetic response solution vector at the j-th measurement point in the i-th iteration. This represents the frequency domain forward electromagnetic response solution vector at the j-th measurement point in the i-th iteration;

[0122] The linear equations for the three-dimensional forward numerical simulation in the time and frequency domains are expressed as follows:

[0123] ;

[0124] ;

[0125] in, This is the time-domain sensitivity matrix. The forward electromagnetic response solution vector in the time domain. The time-domain coefficient matrix, This is the frequency domain sensitivity matrix. The forward electromagnetic response solution vector in the frequency domain. This is the frequency domain coefficient matrix;

[0126] Combining the above equations, the sensitivity matrix can be expressed as:

[0127] ;

[0128] ;

[0129] in, This represents the time-domain sensitivity matrix of the j-th measurement point in the i-th iteration. This represents the frequency domain sensitivity matrix of the j-th measurement point in the i-th iteration. Let represent the prediction model for the i-th iteration;

[0130] By taking the partial derivatives of the model parameters at both ends of the linear equations in the three-dimensional forward numerical simulation in the time and frequency domains, and combining this with the expression for the sensitivity matrix, we obtain:

[0131] ;

[0132] ;

[0133] in, Let be the interpolation basis function for the forward solution vector, defined and This is achieved by solving the following adjoint forward problem:

[0134] ;

[0135] ;

[0136] The sensitivity matrix in the three-dimensional joint inversion equation is obtained as follows:

[0137] ;

[0138] ;

[0139] in, This represents the time-domain sensitivity matrix for the i-th iteration. This represents the frequency domain sensitivity matrix for the i-th iteration. , .

[0140] By solving the aforementioned adjoint forward modeling problem, the Jacobian matrix corresponding to each grid in the three-dimensional joint inversion equation can be calculated. By combining the Jacobian matrices corresponding to each grid, the overall inversion Jacobian matrix can be obtained.

[0141] The initial inversion model is iteratively updated until the set conditions are met, thus completing the joint electromagnetic inversion in the time and frequency domains.

[0142] In this embodiment of the invention, the joint inversion objective function is first performed on the [number]th [year]. Prediction model at the next iteration Taking a second approximation and then calculating the partial derivative, we obtain the gradient as follows:

[0143] ;

[0144] in, This indicates the objective function in the prediction model. The approximate value of the gradient at that point. Is the objective function at the th The gradient of the prediction model at the next iteration For Hessian matrix, This is the model update vector at the k-th iteration;

[0145] Using the approximate positive definite matrix after the k-th iteration Replacing the Hessian matrix The update rule for the quasi-Newton algorithm is obtained as follows:

[0146] ;

[0147] in, Given a rank-2 matrix, the update vector for the initial model can be derived by defining its elements:

[0148] ;

[0149] in, Indicates the first The approximate positive definite matrix at the nth iteration This is the model update vector at the k-th iteration. This is the transpose of the model update vector at the k-th iteration. This represents the gradient difference of the objective function at the k-th iteration. This is the transpose of the gradient difference of the objective function at the k-th iteration.

[0150] Meanwhile, this invention sets termination conditions for the inversion iterative solution. The inversion ends and the model at this point is output when any of the following conditions are met during the inversion process:

[0151] (1) The root mean square error (RMS) between the electromagnetic response of the predicted model and the electromagnetic response of the measured model is less than the error threshold.

[0152] (2) The joint inversion objective function converges;

[0153] (3) Iterate the joint inversion objective function to reach the maximum number of iterations preset.

[0154] In a specific embodiment of the present invention:

[0155] like Figure 3 As shown, examples of inversion results from three inversion algorithms are given. Figure 3 Part a) is a model of a real anomaly under undulating terrain. Figure 3 Part b) shows the time-domain electromagnetic data inversion results. Figure 3 Part c) contains the frequency domain electromagnetic data inversion results. Figure 3 Part d) shows the results of the joint inversion using time-frequency domain electromagnetic data in this invention. As can be seen, Figure 3 The inversion results in part b) do not accurately characterize the detection range of deep anomalies, while Figure 3 Part c) has insufficient detection resolution for shallow anomalies, while Figure 3 The results in section d) show that by performing joint inversion on the two types of data, the characterization of deep anomalies was significantly improved while maintaining the resolution of shallow anomalies. Figure 4The changes in inversion parameters during the three inversion processes are further shown, among which, Figure 4 The upper left section shows the data fit difference decrease curves corresponding to the three inversion processes. Figure 4 The upper right corner shows the objective function curves corresponding to the three inversion processes. Figure 4 The lower left section shows the regularization term curves corresponding to the three inversion processes. Figure 4 The lower right corner shows the regularization factor curves corresponding to the three inversion processes. As can be seen, all three inversion methods can converge stably.

[0156] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A three-dimensional time-frequency domain electromagnetic joint inversion method based on the finite element algorithm, characterized in that, include: Time-domain electromagnetic observation data and frequency-domain electromagnetic observation data of the area to be measured are acquired respectively. The region to be tested is divided into unstructured tetrahedral meshes, and an initial time-frequency domain inversion model is constructed based on the unstructured tetrahedral meshes. The time-frequency domain inversion initial model is subjected to forward modeling calculation using a forward modeling algorithm to obtain time-domain electromagnetic response data and frequency-domain electromagnetic response data. The residuals of time-domain electromagnetic data are obtained from time-domain electromagnetic observation data and time-domain electromagnetic response data, and the residuals of frequency-domain electromagnetic data are obtained from frequency-domain electromagnetic observation data and frequency-domain electromagnetic response data. Adaptive weights are obtained for the time-domain electromagnetic data residuals and the frequency-domain electromagnetic data residuals, respectively; specifically: The expression for the adaptive weights of the time-domain electromagnetic data residuals is: ； The expression for the adaptive weight of the frequency domain electromagnetic data residual is: ； in, , representing the residual of electromagnetic data in the time domain; , representing the residual of electromagnetic data in the frequency domain; A joint inversion objective function is established based on the time-domain electromagnetic data residuals and the frequency-domain electromagnetic data residuals and their corresponding adaptive weights. The joint inversion objective function is expressed as: ； In the formula, For the joint inversion objective function, The time-domain data covariance matrix, For the frequency domain data covariance matrix, To constrain the covariance matrix of the model, As a regularization factor, For prediction models, This is the initial model for time-frequency domain inversion. For time-domain electromagnetic observation data, For frequency domain electromagnetic observation data, For time-domain electromagnetic response data, For frequency domain electromagnetic response data, For adaptive weights of time-domain electromagnetic data residuals, For adaptive weights of frequency domain electromagnetic data residuals, express Norm; The joint inversion objective function is solved iteratively based on the quasi-Newton algorithm, and the update vector of the initial inversion model is calculated. The initial inversion model is iteratively updated until the set conditions are met, thus completing the joint electromagnetic inversion in the time and frequency domains.

2. The three-dimensional time-frequency domain electromagnetic joint inversion method based on the finite element algorithm according to claim 1, characterized in that: The time-domain electromagnetic observation data and frequency-domain electromagnetic observation data of the area to be measured were acquired separately, as follows: In the area to be tested, array transceivers corresponding to the time domain and frequency domain are respectively deployed. The distance between the array transceiver corresponding to the time domain and the area to be tested is a first predetermined distance, and the distance between the array transceiver corresponding to the frequency domain and the area to be tested is a second predetermined distance. The first predetermined distance is less than the second predetermined distance. The corresponding time-domain electromagnetic observation data and frequency-domain electromagnetic observation data are acquired in the region under test through array transceivers corresponding to the time domain and frequency domain, respectively.

3. The three-dimensional time-frequency domain electromagnetic joint inversion method based on the finite element algorithm according to claim 1, characterized in that, The joint inversion objective function is solved iteratively based on a quasi-Newton algorithm, including: The objective function is expanded using Taylor and its minimum value is taken to obtain the three-dimensional joint inversion equation; Calculate the gradient of the joint inversion objective function by taking partial derivatives of the three-dimensional joint inversion equation; The inversion calculation is performed using a quasi-Newton algorithm based on the gradient of the joint inversion objective function to obtain the update vector of the initial inversion model.

4. The three-dimensional time-frequency domain electromagnetic joint inversion method based on the finite element algorithm according to claim 3, characterized in that: The three-dimensional joint inversion equation is expressed as follows: ； in, Denotes the approximate positive definite matrix after the k-th iteration. To update the vector, This is the time-domain sensitivity matrix. This is the frequency domain sensitivity matrix. The time-domain data covariance matrix, This is the transpose of the covariance matrix of the time-domain data. The residual between the observed data and the response data, For the frequency domain data covariance matrix, This is the transpose of the frequency domain data covariance matrix. The model constraint variance matrix, This is the transpose of the model constraint variance matrix. For prediction models, This is the initial model for inversion.

5. The three-dimensional time-frequency domain electromagnetic joint inversion method based on the finite element algorithm according to claim 4, characterized in that, The iterative solution of the joint inversion objective function also includes: calculating the sensitivity matrix in the three-dimensional joint inversion equation using the adjoint forward modeling method, specifically: The sensitivity matrix in the three-dimensional joint inversion equation includes: the time-domain sensitivity matrix and the frequency-domain sensitivity matrix; Time-domain sensitivity matrix and frequency domain sensitivity matrix They are defined as follows: ； ； in, This represents a time-domain prediction model. This represents a frequency domain prediction model. and The three-dimensional electromagnetic numerical simulation responses in the time domain and frequency domain, respectively, can be defined as follows: ； ； in, Let be the interpolation basis function for the forward solution vector of the i-th iteration and j-th measurement point of the time-domain electromagnetic response. Let be the interpolation basis function for the forward solution vector of the frequency domain electromagnetic response at the j-th measurement point in the i-th iteration. This represents the time-domain forward electromagnetic response solution vector at the j-th measurement point in the i-th iteration. This represents the frequency domain forward electromagnetic response solution vector at the j-th measurement point in the i-th iteration; The linear equations for the three-dimensional forward numerical simulation in the time and frequency domains are expressed as follows: ； ； in, This is the time-domain sensitivity matrix. The forward electromagnetic response solution vector in the time domain. The time-domain coefficient matrix, This is the frequency domain sensitivity matrix. The forward electromagnetic response solution vector in the frequency domain. This is the frequency domain coefficient matrix; combining the above equations, the sensitivity matrix can be expressed as: ； ； in, This represents the time-domain sensitivity matrix of the j-th measurement point in the i-th iteration. This represents the frequency domain sensitivity matrix of the j-th measurement point in the i-th iteration. Let represent the prediction model for the i-th iteration; By taking the partial derivatives of the model parameters at both ends of the linear equations in the three-dimensional forward numerical simulation in the time and frequency domains, and combining this with the expression for the sensitivity matrix, we obtain: ； ； in, Let be the interpolation basis function for the forward solution vector, defined This is achieved by solving the following adjoint forward problem: ； ； The sensitivity matrix in the three-dimensional joint inversion equation is obtained as follows: ； ； in, This represents the time-domain sensitivity matrix for the i-th iteration. This represents the frequency domain sensitivity matrix for the i-th iteration. .

6. The three-dimensional time-frequency domain electromagnetic joint inversion method based on the finite element algorithm according to claim 1, characterized in that, The setting conditions include: The root mean square error (RMS) between the predicted electromagnetic response and the measured electromagnetic response is less than the error threshold. The joint inversion objective function converged; The joint inversion objective function is solved iteratively to reach the preset maximum number of iterations.

7. The three-dimensional time-frequency domain electromagnetic joint inversion method based on the finite element algorithm according to claim 1, characterized in that, The initial inversion model is updated in each iteration, including: For the joint inversion objective function at the th Prediction model at the next iteration Taking a second approximation and then calculating the partial derivative, we obtain the gradient as follows: ； in, This indicates the objective function in the prediction model. The approximate value of the gradient at that point. Is the objective function at the th The gradient of the prediction model at the next iteration For Hessian matrix, , This is the model update vector at the k-th iteration; Using the approximate positive definite matrix after the k-th iteration Replacing Hessian matrix The update rule for the quasi-Newton algorithm is obtained as follows: ； in, Given a rank-2 matrix, the update vector for the initial model can be derived by defining its elements: ； in, Indicates the first The approximate positive definite matrix at the nth iteration This is the model update vector at the k-th iteration. This is the transpose of the model update vector at the k-th iteration. This represents the gradient difference of the objective function at the k-th iteration. This is the transpose of the gradient difference of the objective function at the k-th iteration.

Images