An inversion method and system for geophysical electromagnetic method data based on entropy constraint

Abstract

This application belongs to the technical field of geophysical methods, and specifically discloses a geophysical electromagnetic method data inversion method and system based on entropy constraint. The method includes: setting conductivity values after dividing the inversion space to construct an initial inversion model; constructing a magnetotelluric regularization inversion objective function based on entropy constraint; using the Gauss-Newton method to minimize the magnetotelluric regularization inversion objective function, updating the initial inversion model, and obtaining a conductivity model to achieve three-dimensional magnetotelluric inversion. This application can achieve the effect of focused inversion to obtain a high-resolution conductivity model with clear boundaries, serving fields such as oil and mineral resource exploration.

Description

Technical Field

[0001] This application belongs to the technical field of geophysical methods, and more specifically, relates to a geophysical electromagnetic method data inversion method and system based on entropy constraint. Background Art

[0002] As an important geophysical method, magnetotelluric method has been widely applied in many fields such as mineral resource exploration, geothermal resource exploration, and deep earth structure detection. With the gradual deepening of electromagnetic exploration into areas with complex geological structures and the increasing demand for fine exploration, high-precision three-dimensional electromagnetic exploration has gradually become a research hotspot. The fine interpretation of large-scale magnetotelluric data requires the development of efficient and stable three-dimensional forward and inversion algorithms.

[0003] High-precision three-dimensional forward algorithms are the basis of inversion algorithms, and the accuracy and speed of forward algorithms will have a greater impact on the efficiency of inversion algorithms. The development of efficient inversion algorithms first requires the formation of high-precision and high-efficiency forward algorithms. The integral equation method and the finite difference method are applied to electromagnetic three-dimensional forward simulation. Different methods have different advantages and disadvantages. The finite difference method has a simple principle and is easy to implement algorithmically, but it is difficult to accurately divide complex geometric structures, which is an irreparable shortcoming. The simulation accuracy of the integral equation method depends severely on the accuracy of the coefficient matrix and the solution of the dyadic Green's function. For large-scale complex problems, solving the linear equations greatly increases the memory requirements, with slow solving speed and reduced solving accuracy. The finite element method can simulate complex terrain, geological bodies, and reflection sources of arbitrary shapes through the use of flexible unstructured grid division technology, with very high calculation accuracy, and is most suitable for fine detection, but the calculation efficiency still needs to be improved. The finite volume method can also discretize the research area using unstructured grids, but its accuracy is still not as good as that of the finite element method.

[0004] The inversion algorithm is an indispensable link in electromagnetic data processing and interpretation. Although the magnetotelluric three-dimensional inversion based on Tikhonov regularization can obtain a stable solution, when there are sharp boundaries in the underground electrical anomaly bodies, the boundaries of the electrical anomaly bodies obtained by the traditional magnetotelluric three-dimensional inversion based on norm regularization are relatively smooth, which brings certain difficulties to subsequent geological and geophysical interpretations. To address this problem, existing non-smooth inversion methods based on non- norm regularization have been proposed. Most non-smooth inversion methods use a certain function as the physical property model and obtain the model solution with sharp boundaries by minimizing the objective function, such as the minimum support functional and the minimum gradient support functional. However, their inversion results depend on the selection of the focusing factor and are greatly affected by human experience, lacking a certain degree of stability. Summary of the Invention

[0005] Aiming at the deficiencies of the existing technology, the purpose of this application is to provide a geophysical electromagnetic method data inversion method and system based on entropy constraint, aiming to solve the problems that the continuous diffusion of the existing magnetotelluric smooth inversion results and the lack of stability in the focusing inversion.

[0006] To achieve the above object, in the first aspect, this application provides a geophysical electromagnetic method data inversion method based on entropy constraint, including the following steps:

[0007] Step 1: After dividing the inversion space, set the conductivity value to construct an initial inversion model;

[0008] Step 2: Based on the initial inversion model, construct a magnetotelluric regularization inversion objective function based on entropy constraint;

[0009] Step 3: Use the Gauss-Newton method to minimize the magnetotelluric regularization inversion objective function, update the initial inversion model, obtain the conductivity model, and realize the three-dimensional magnetotelluric inversion.

[0010] Further preferably, in step 2, the magnetotelluric regularization inversion objective function is:

[0011]

[0012] Among them, is the inversion model parameter, and the inversion model parameter is the logarithm of the conductivity; is the data fitting term; is the model roughness term; is the minimum entropy term; and are regularization parameters; , , is the model parameter of the th inversion unit, , is the number of inversion units, is a very small positive number, with an order of magnitude of .

[0013] Further preferably, in step 2, the approximate spectrum analysis method is used to calculate the regularization parameter:

[0014]

[0015]

[0016] Among them, is a random vector, is the sensitivity of the observed magnetotelluric response data to the inversion model parameter, which is calculated by the adjoint forward method, is the number of iterations, is a positive integer less than 3, is an empirical parameter, ; is the data weighting matrix; is the model roughness matrix; is the diagonal minimum entropy matrix that plays a role in focusing inversion.

[0017] Further preferably, step three specifically includes the following steps:

[0018] Step 3.1: Minimize the magnetotelluric regularization inversion objective function using the Gauss-Newton method with a quasi-quadratic convergence rate, and obtain the normal equation using the Taylor expansion formula;

[0019] Step 3.2: Convert the normal equation into the least squares form, use the least squares algorithm to obtain the inversion model update amount, and use the linear search method to obtain the model update step size, and then obtain the prediction model;

[0020] Step 3.3: When the prediction model meets the preset conditions, use the prediction model as the conductivity model.

[0021] Further preferably, step 3.3 is specifically:

[0022] Step 3.3.1: Conduct three-dimensional magnetotelluric forward modeling of the prediction model based on the vector finite element method to obtain the magnetotelluric response data of the prediction model;

[0023] Step 3.3.2: Use the normalized error to evaluate the fitting of the magnetotelluric response data of the prediction model and the observed magnetotelluric response data. If the normalized fitting difference reaches the set threshold, use the prediction model as the conductivity model; otherwise, go to step 3.3.3:

[0024] Step 3.3.3: Determine whether the current iteration is the maximum number of inversion iterations. If so, use the prediction model as the conductivity model; otherwise, go to step 3.2.

[0025] Further preferably, step one specifically includes the following steps:

[0026] Step 1.1: Discretize the inversion space using unstructured tetrahedral meshes;

[0027] Step 1.2: Set the conductivity values in the unstructured tetrahedral meshes obtained in step 1.1 to construct an initial inversion model.

[0028] Further preferably, the normal equation is:

[0029]

[0030]

[0031] Among them, is the model update amount; is the data weighting matrix; is the model roughness matrix; is the diagonal minimum entropy matrix that plays a role in focusing inversion; and are regularization parameters; is the sensitivity matrix; is the data error between the predicted value and the observed value after the th iteration; is the model parameter after the th iteration;

[0032] In a second aspect, the present application provides a geophysical electromagnetic method data inversion system based on entropy constraint, including:

[0033] A model construction module that sets conductivity values after dividing the inversion space and constructs an initial inversion model;

[0034] A target function construction module for constructing a magnetotelluric regularization inversion target function based on entropy constraint based on the initial inversion model;

[0035] A model update module for minimizing the magnetotelluric regularization inversion target function using the Gauss-Newton method, updating the initial inversion model, and obtaining a conductivity model to achieve three-dimensional magnetotelluric inversion.

[0036] Further preferably, the magnetotelluric regularization inversion target function in the target function construction module is:

[0037]

[0038] Among them, is the inversion model parameter, and the inversion model parameter is the logarithm of conductivity; is the data fitting term; is the model roughness term; is the minimum entropy term; and are regularization parameters; , , is the model parameter of the th inversion unit, , is a very small positive number with a magnitude of .

[0039] Further preferably, the objective function construction module calculates the regularization parameter by using the approximate spectral analysis method:

[0040]

[0041]

[0042] wherein, is a random vector, is the sensitivity of the observed magnetotelluric response data to the inversion model parameters, which is calculated by using the adjoint forward method, is the number of iterations, is a positive integer less than 3, is an empirical parameter, ; is the data weighting matrix; is the model roughness matrix; is the diagonal minimum entropy matrix that plays a role in focusing inversion.

[0043] Further preferably, the model update module includes:

[0044] The normal equation acquisition unit is used to minimize the magnetotelluric regularization inversion objective function by using the Gauss-Newton method with a quasi-quadratic convergence rate, and obtain the normal equation by using the Taylor expansion formula; the prediction model acquisition unit is used to convert the normal equation into the least squares form, adopt the least squares algorithm, obtain the inversion model update amount, and adopt the linear search method to obtain the model update step size, and then obtain the prediction model;

[0045] The condition determination unit is used to use the prediction model as the conductivity model when the prediction model meets the preset conditions.

[0046] Further preferably, in the normal equation acquisition unit, the normal equation is:

[0047]

[0048]

[0049] wherein, is the model update amount; is the data weighting matrix; is the model roughness matrix; is the diagonal minimum entropy matrix that plays a role in focusing inversion; 、 are the regularization parameters; is the sensitivity matrix; is the data error between the predicted value and the observed value after the th iteration; The model parameters after the next iteration; is the prior model.

[0050] Further preferably, the condition determination unit includes:

[0051] A data collector, configured to perform three-dimensional magnetotelluric forward modeling on the prediction model based on the vector finite element method to obtain magnetotelluric response data of the prediction model;

[0052] A prediction model condition determiner, configured to use the normalized error to evaluate the fitting condition between the magnetotelluric response data of the prediction model and the observed magnetotelluric response data. If the normalized fitting difference reaches a set threshold, the prediction model is used as the conductivity model;

[0053] An iteration number determiner, configured to determine whether the current iteration is the maximum number of inversion iterations. If so, the prediction model is used as the conductivity model.

[0054] Further preferably, the model construction module includes:

[0055] A mesh generation unit, configured to mesh the inversion space using unstructured tetrahedral meshes;

[0056] A conductivity filling unit, configured to set conductivity values in the unstructured tetrahedral meshes to construct an initial inversion model.

[0057] In a third aspect, the present application provides an electronic device, including: at least one memory for storing a program; at least one processor for executing the program stored in the memory. When the program stored in the memory is executed, the processor is configured to execute the method described in the first aspect or the further preferred of the first aspect.

[0058] In a fourth aspect, the present application provides a computer-readable storage medium storing a computer program. When the computer program runs on a processor, the processor is caused to execute the method described in the first aspect or any further preferred of the first aspect.

[0059] In a fifth aspect, the present application provides a computer program product. When the computer program product runs on a processor, the processor is caused to execute the method described in the first aspect or any further preferred of the first aspect.

[0060] It can be understood that the beneficial effects of the above second aspect to fifth aspect can refer to the relevant descriptions in the first aspect, and will not be elaborated here.

[0061] Generally speaking, compared with the prior art, the above technical solutions conceived by the present application have the following beneficial effects:

[0062] The present application provides a geophysical electromagnetic method data inversion method based on entropy constraint. During the inversion process, based on the idea of minimum entropy regularization, the minimum entropy stable functional is added to the objective function and expressed in the form of a pseudo - quadratic functional of the model parameters. Combining the characteristics of unstructured grids, a suitable diagonal minimum entropy weighted matrix is constructed by introducing the volume of tetrahedral elements. The high - speed Newton method equation with a quasi - quadratic convergence rate is used to transform the objective function into a least - squares problem, which has a quasi - quadratic convergence rate, to obtain the model update amount, and a line search is used to obtain the optimal model update step size. It can not only improve the resolution ability of the inversion result, achieve a focused inversion effect, and obtain a clear physical property boundary, but also does not require determining reasonable focusing parameters, thus reducing the influence of human factors and having good stability and applicability. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 is a flowchart of the geophysical electromagnetic method data inversion method based on entropy constraint provided by an embodiment of the present application;

[0064] Figure 2(a) is a slice diagram of y = 0m obtained by using a theoretical model provided by an embodiment of the present application;

[0065] Figure 2(b) is a slice diagram of y = 0m obtained by using smooth inversion provided by an embodiment of the present application;

[0066] Figure 2(c) is a slice diagram of y = 0m obtained by using minimum entropy constraint inversion provided by an embodiment of the present application;

[0067] Figure 2(d) is a slice diagram of z = 500m obtained by using a theoretical model provided by an embodiment of the present application;

[0068] Figure 2(e) is a slice diagram of z = 500m obtained by using smooth inversion provided by an embodiment of the present application;

[0069] Figure 2(f) is a slice diagram of z = 500m obtained by using minimum entropy constraint inversion provided by an embodiment of the present application;

[0070] Figure 3 is a schematic diagram of an electronic device provided by an embodiment of the present application. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0071] In order to make the objectives, technical solutions and advantages of the present application clearer, the present application will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0072] As used herein, the term "and / or" describes the associated relationship of associated objects, indicating that three relationships may exist. For example, A and / or B may represent three cases: A exists alone, A and B exist simultaneously, and B exists alone. As used herein, the symbol " / " indicates that the associated objects are in an "or" relationship. For example, A / B means A or B.

[0073] Terms such as "first" and "second" in the description and claims of this application are used to distinguish different objects, rather than to describe a specific order of the objects.

[0074] In the embodiments of this application, words such as "exemplary" or "for example" are used to indicate examples, illustrations, or explanations. Any embodiment or design solution described as "exemplary" or "for example" in the embodiments of this application should not be construed as being more preferred or having more advantages than other embodiments or design solutions. Rather, the use of words such as "exemplary" or "for example" is intended to present the relevant concepts in a specific manner.

[0075] In the description of the embodiments of this application, unless otherwise specified, the meaning of "a plurality" refers to two or more.

[0076] The embodiments of this application will be described below in conjunction with the accompanying drawings in the embodiments of this application.

[0077] As Figure 1 shown, this application provides a geophysical electromagnetic data inversion method based on entropy constraint, including the following steps:

[0078] Step S1: Obtain the observed magnetotelluric response data for inversion, including collecting the actual observed magnetotelluric response data on the ground and obtaining the theoretical observed magnetotelluric response data through three-dimensional magnetotelluric forward modeling using a theoretical model; the theoretical model is a simulation model capable of obtaining the theoretical observed magnetotelluric response data;

[0079] Among them, in this application, three-dimensional magnetotelluric forward modeling uses the vector finite element method with unstructured tetrahedral meshes; specifically:

[0080] Since the frequency band used in the magnetotelluric method is the low-frequency band, the influence of displacement current can be ignored, and the frequency-domain Maxwell equations can be obtained:

[0081] (1)

[0082] (2)

[0083] (3)

[0084] (4)

[0085] Among them, is the magnetic field strength; is the electric field strength; is the magnetic permeability in vacuum; is the angular frequency; is the imaginary unit; is the conductivity; Taking the curl of both sides of formula (2), we can get:

[0086] (5)

[0087] Substituting formula (1) into formula (5), we can get:

[0088] (6)

[0089] Using the secondary field algorithm to solve the 3D magnetotelluric forward problem, let:

[0090] (7)

[0091] Among them, is the background electric field of the layered medium, which is solved by the analytical method; is the value of the secondary electric field; Substituting formula (7) into formula (6), we get the control equation satisfied by the secondary electric field:

[0092] (8)

[0093] Among them, is the conductivity of the background model.

[0094] Adopting the simplest Dirichlet boundary condition: ;

[0095] Among them, is the normal vector at the boundary.

[0096] Using the weighted residual method to solve the above formula, taking the vector basis function as the weight function, we can get the following residual function :

[0097] (9)

[0098] Among them, is the vector basis function of each calculation unit; is the volume of the calculation unit; is the calculation unit domain;

[0099] Then applying the vector Green's first theorem, we can get the discrete form of formula (9):

[0100] (10)

[0101] Among them, and are calculation units for the secondary field and the primary field on the edges; is the total number of calculation units; and are calculation units of the stiffness matrix; It can be expressed as:

[0102]

[0103]

[0104] Among them, and are vector basis functions within element e.

[0105] Finally, formula (10) can be assembled into the following form of a large sparse equation set:

[0106] (11)

[0107] Among them, is the global stiffness matrix; is the right-hand side term of the loading boundary condition; is the secondary field of the edge to be solved. Here, a direct solver based on matrix factorization is used to solve this sparse equation set.

[0108] After obtaining the secondary electric field by solving equation (11), the value of the secondary magnetic field is obtained according to Ampere's law . After obtaining the secondary magnetic field, adding it to the background magnetic field can obtain the total magnetic field value:

[0109] (12)

[0110] (13)

[0111] Among them, is the magnetic field response of the one-dimensional layered medium and can be analytically solved; When the electric field and magnetic field are obtained, the impedance is obtained according to the following calculation formula :

[0112] (14)

[0113] Among them, the subscript 1 and subscript 2 respectively refer to and the electromagnetic field data obtained by the polarization method; and and and are the impedance tensor values; and is the value of the magnetic field component; , is the value of the electric field component; , are two mutually orthogonal directions.

[0114] Step S2: Discretize the inversion space using unstructured tetrahedral meshes; wherein, the inversion space is the underground space corresponding to the measuring points and the preset extension range nearby;

[0115] This application uses the unstructured tetrahedral mesh method to discretize the inversion space. By using unstructured meshes, it is easy to achieve local refinement of key areas (such as measuring points); refinement around the measuring points can obtain very accurate simulation results;

[0116] Step S3: Set the conductivity value in the unstructured tetrahedral meshes obtained in Step S2 to construct an initial inversion model;

[0117] In 3D inversion, the selection of the initial inversion model will affect the inversion effect and convergence speed. Selecting an initial inversion model that is quite different from the true model may lead to non-convergence of magnetotelluric inversion and poor inversion results; therefore, constructing a suitable initial inversion model is of great significance for obtaining high-quality inversion results; generally, a simple initial inversion model can be established with reference to existing drilling, geological, and geophysical prior information, but in areas with less prior information and complex geological conditions, 1D inversion can be carried out first, and the 1D inversion result can be used as the initial model;

[0118] Step S4: Construct a magnetotelluric regularization inversion objective function based on entropy constraint;

[0119] The traditional magnetotelluric 3D inversion based on norm regularization often gives a relatively smooth boundary of the anomaly body. Using the minimum entropy regularization objective function can obtain a solution with a focused inversion effect; the regularization objective function is:

[0120] (15)

[0121] wherein, is the data fitting term; is the model roughness term; is the minimum entropy term; , are the regularization parameters;

[0122] , and are respectively defined as:

[0123] (16)

[0124] (17)

[0125] (18)

[0126] Express as the form of a pseudo - quadratic functional of model parameters:

[0127] (19)

[0128] Among them, is the inversion model parameter. To ensure that is non - negative during the inversion process and endows the model parameters with physical meanings, the conductivity in the logarithmic domain is taken as the model parameter; is the data weighting matrix; represents the forward response of the initial inversion model parameter ; is the observed magnetotelluric response data at the measurement point; is the model roughness matrix; is the diagonal minimum entropy matrix that plays a role in focusing inversion; is the reference model parameter containing prior information; , is the model parameter of the th inversion element, , is the number of inversion elements, is a very small positive number, usually taken as .

[0129] Combined with the unstructured grid characteristics, a suitable diagonal minimum entropy weighting matrix is constructed by introducing the tetrahedron element volume:

[0130] (20)

[0131] Among them, , , is the tetrahedron element volume; is the number of inversion elements, is the model parameter of the th inversion element, is the prior model of the th inversion element. By adding the tetrahedron element volume, the influence of different element volume sizes in the unstructured grid on the inversion result is eliminated.

[0132] According to the volume and distance weights of the inversion elements closest to the inversion element, construct the roughness matrix, which is expressed as follows:

[0133] ​ (21)

[0134] (22)

[0135] Among them, is the element in the th row and th column of the roughness matrix, is the element number of the th unit closest to the th inversion unit in the inversion grid; , are respectively the volumes of the th and th and th units closest to the is the number of adjacent units, is the volume of the th inversion unit; represents the distance from the center of the th inversion unit to the center of the th closest unit; and are respectively the center coordinates of the th inversion unit and the th inversion unit.

[0136] The regularization parameter plays an important role in balancing the data fitting term and the model constraint terms (model roughness term and minimum entropy term). The choice of the regularization parameter affects the stability and convergence rate of the algorithm. The regularization parameter is calculated using the approximate spectral analysis method:

[0137] (23)

[0138] (24)

[0139] Among them, is a random vector, is the sensitivity of the observed magnetotelluric response data to the inversion model parameters, is the number of iterations, is a positive integer less than 3, is an empirical parameter, usually taking ;

[0140] Step S5: Calculate the sensitivity matrix using the adjoint forward method;

[0141] First, obtain the sensitivity matrix of the electromagnetic field components at the measurement points under different polarization modes, and then obtain the sensitivity matrix ; Sensitivity matrix of electromagnetic field components Can be written as:

[0142] (25)

[0143] Taking the partial derivative of both sides with respect to the model parameters gives:

[0144] (26)

[0145] Among them, Is the global stiffness matrix; Is the right - hand side of the loading boundary condition.

[0146] Therefore, the partial derivative of the secondary field with respect to the model parameters can be obtained:

[0147] (27)

[0148] Substituting into the above formula gives:

[0149] (28)

[0150] Among them The order of Is the number of edges, Is the number of inversion units); Since The number of columns of , In order to solve the sensitivity matrix , It is necessary to do Forward modeling calculations, which will consume a large amount of computing time; Since the coefficient matrix Is symmetric, the above formula can be transposed to get:

[0151] (29)

[0152] In this way, when solving Only need to perform (Number of measurement points) quasi - forward modeling calculations. Since the number of measurement points is much smaller than the number of model parameters, it has higher computational efficiency;

[0153] Step S6: Use the LSQR method to solve the least - squares problem of the equivalent Gauss - Newton method equation, obtain the model update amount, linearly search to obtain the optimal model update step size, and update the prediction model;

[0154] Using the Gauss - Newton method with quasi - quadratic convergence speed to minimize the objective function, applying the Taylor expansion formula to the objective function and ignoring the high - order terms can obtain the following normal equation:

[0155]

[0156] (30)

[0157] Among them, is the model update amount; is the data weighting matrix; is the model roughness matrix; is the diagonal minimum entropy matrix that plays a role in focusing inversion; and are regularization parameters; is the sensitivity matrix; is the data error between the predicted value and the observed value after the th iteration; is the model parameter after the th iteration;

[0158] The normal equation is converted into the following least squares form:

[0159] (31)

[0160] The least squares (LSQR) algorithm is used to solve Equation (31) to obtain the model update amount, and the linear search method is used to find the optimized model update step size , and the model update equation can be expressed as:

[0161] (32)

[0162] Step S7: Based on the vector finite element method, perform 3D magnetotelluric forward modeling on the prediction model to obtain the magnetotelluric response data of the prediction model; in the forward modeling, a direct solver, MPI (Message Passing Interface), and OpenMP parallel computing technology are combined to improve the computing efficiency;

[0163] Step S8: Use the normalized error to evaluate the fitting of the model prediction data and the observed data, and judge whether the preset threshold is reached and whether the maximum number of iterations is reached;

[0164] Judge whether the observed magnetotelluric response data and the predicted magnetotelluric response data are fitted. The data normalized fitting difference between the observed magnetotelluric response data and the predicted magnetotelluric response data can be expressed by the following formula:

[0165] (33)

[0166] When the data normalization fitting error reaches the set threshold or the convergence stops, the inversion algorithm terminates and outputs the inversion result; if the set minimum data normalization fitting error is not reached, it is judged whether the maximum number of inversion iterations is reached. If the current iteration number is less than the set number, the next iteration is looped. If the set number is reached, the inversion iteration stops and the inversion result is output;

[0167] Step S9: Iteratively execute steps S5 to S8 until the set termination condition is reached and exit, obtaining the conductivity model to realize the 3D magnetotelluric inversion.

[0168] Embodiment 1

[0169] To verify the effectiveness of the magnetotelluric inversion method based on minimum entropy constraint provided by this application, an embodiment of this application provides a two-low-resistance model with a mountain peak terrain. The height of the mountain peak is 300m relative to the ground, the side length of its upper base is 450m, the side length of its lower base is 1800m, and the resistivity of the background half-space is 100 , the sizes of the two abnormal bodies are both 1000m×1000m×500m, and the resistivities are both 10 , and the center coordinates of the abnormal bodies are (-1000m, 0m, 500m) and (1000m, 0m, 500m) respectively. A total of 121 measuring points are arranged, the number of frequencies is 5 (100Hz~1.668Hz), the number of forward modeling grids is 350418, and 2% Gaussian noise is added to the forward modeling impedance data to simulate the measured data. During the inversion process, the number of grids in the inversion area is 175816, and both the initial model and the prior model are set to 100 of a uniform half-space. The maximum number of Gauss-Newton iterations is set to 40, and the inversion stops when the inversion reaches the maximum number of iterations or the calculated normalized error is less than the specified error; Figures 2(a) and 2(d) are theoretical models, Figures 2(b) and 2(e) are smooth inversion results, Figures 2(c) and 2(f) are minimum entropy inversion results, the white dots in the figures are the positions of the measuring points, and the black rectangular frames are the positions of the abnormal bodies. It can be seen that the minimum entropy constraint inversion can better recover the resistivity value of the abnormal body and the boundary of the abnormal body.

[0170] Embodiment 2

[0171] This application provides a geophysical electromagnetic method data inversion system based on entropy constraint, including:

[0172] A model construction module, used to set the conductivity value after dividing the inversion space and construct an initial inversion model;

[0173] A target function construction module, used to construct a magnetotelluric regularization inversion target function based on entropy constraint based on the initial inversion model;

[0174] A model update module, which is used to minimize the magnetotelluric regularization inversion objective function by using the Gauss-Newton method, update the initial inversion model, obtain the conductivity model, and realize three-dimensional magnetotelluric inversion.

[0175] Further preferably, in the objective function construction module, the magnetotelluric regularization inversion objective function is:

[0176]

[0177] where, is the inversion model parameter, and the inversion model parameter is the logarithm of the conductivity; is the data fitting term; is the model roughness term; is the minimum entropy term; , are the regularization parameters; , , is the model parameter of the th inversion unit, , is the number of inversion units, is a very small positive number, usually taken as .

[0178] Further preferably, the objective function construction module calculates the regularization parameter by using the approximate spectral analysis method:

[0179]

[0180]

[0181] where, is a random vector, is the sensitivity of the observed magnetotelluric response data to the inversion model parameter, which is calculated by using the adjoint forward method, is the number of iterations, is a positive integer less than 3, is an empirical parameter, ; is the data weighting matrix; is the model roughness matrix; is the diagonal minimum entropy matrix that plays a role in focusing inversion.

[0182] Further preferably, the model update module includes:

[0183] The normal equation acquisition unit is used to minimize the magnetotelluric regularization inversion objective function by using the Gauss-Newton method with a quasi-quadratic convergence rate, and obtain the normal equation by using the Taylor expansion formula; the prediction model acquisition unit is used to convert the normal equation into the least squares form, adopt the least squares algorithm to obtain the inversion model update amount, and adopt the linear search method to obtain the model update step size, and then obtain the prediction model;

[0184] The condition determination unit is used to use the prediction model as the conductivity model when the prediction model meets the preset conditions.

[0185] Further preferably, in the normal equation acquisition unit, the normal equation is:

[0186]

[0187]

[0188] Wherein, is the model update amount; is the data weighting matrix; is the model roughness matrix; is the diagonal minimum entropy matrix that plays a role in focusing inversion; 、 are regularization parameters; is the sensitivity matrix; is the data error between the predicted value and the observed value after the th iteration; is the model parameter after the th iteration; is the prior model.

[0189] Further preferably, the condition determination unit includes:

[0190] The data collector is used to perform three-dimensional magnetotelluric forward modeling on the prediction model based on the vector finite element method to obtain the magnetotelluric response data of the prediction model;

[0191] The prediction model condition determiner is used to evaluate the fitting situation between the magnetotelluric response data of the prediction model and the observed magnetotelluric response data by using the normalized error. If the normalized fitting difference reaches the set threshold, the prediction model is used as the conductivity model;

[0192] The iteration number determiner is used to determine whether the current iteration is the maximum number of inversion iterations. If so, the prediction model is used as the conductivity model.

[0193] Further preferably, the model construction module includes:

[0194] The grid meshing unit is used to mesh the inversion space by using unstructured tetrahedral grids;

[0195] A conductivity filling unit is used to construct an initial inversion model in an unstructured tetrahedral mesh.

[0196] In summary, compared with the prior art, the present application has the following advantages:

[0197] The present application provides a magnetotelluric focusing inversion method. When acquiring the observed magnetotelluric response data, based on the secondary field theory, starting from the frequency-domain Maxwell's equations, a higher-order partial differential equation satisfied by the secondary electric field is derived. The weighted residual method is used to transform the differential equation into a finite element equation in integral form; the Dirichlet boundary condition is adopted, assuming that the electric field vanishes on the boundary of the modeling region, setting the condition that the electric field is zero, and expanding the modeling region of the theoretical model, densifying the grid near the measuring points, and gradually increasing the grid size in other regions to complete the setting of the boundary condition and the grid; the magnetotelluric response of complex geological structures can be accurately simulated.

[0198] During the inversion process, based on the idea of minimum entropy regularization, the minimum entropy stability functional is added to the objective function and expressed in the form of a pseudo-quadratic functional of the model parameters. Combining the characteristics of the unstructured grid, a suitable diagonal minimum entropy weighting matrix is constructed by introducing the volume of the tetrahedral element. The high-speed Newton method equation with a quasi-quadratic convergence rate is used to transform the objective function into a least-squares problem, which has a quasi-quadratic convergence rate, obtaining the model update amount, and a line search is used to obtain the optimal model update step size. It can not only improve the resolution ability of the result, achieve the focusing inversion effect, and obtain clear physical property boundaries, but also does not require determining reasonable focusing parameters, thereby reducing the influence of human factors and having good stability and applicability.

[0199] It should be understood that the above system is used to execute the method in the above embodiment. For the corresponding program modules in the system, their implementation principles and technical effects are similar to those described in the above method. The working process of the system can refer to the corresponding process in the above method and will not be elaborated here.

[0200] Based on the method in the above embodiment, as Figure 3 shown, the embodiment of the present application provides an electronic device, which may include: a processor (Processor) 810, a communication interface (Communications Interface) 820, a memory (Memory) 830, and a communication bus 840. Among them, the processor 810, the communication interface 820, and the memory 830 communicate with each other through the communication bus 840. The processor 810 can call the logical instructions in the memory 830 to execute the method in the above embodiment.

[0201] In addition, when the logical instructions in the above-mentioned memory 830 are implemented in the form of software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on such an understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a part of this technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions for causing a computer device (which may be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in various embodiments of this application.

[0202] Based on the method in the above embodiment, an embodiment of this application provides a computer-readable storage medium. The computer-readable storage medium stores a computer program. When the computer program runs on a processor, it causes the processor to execute the method in the above embodiment.

[0203] Based on the method in the above embodiment, an embodiment of this application provides a computer program product. When the computer program product runs on a processor, it causes the processor to execute the method in the above embodiment.

[0204] It can be understood that the processor in the embodiments of this application may be a central processing unit (CPU), or may also be other general-purpose processors, digital signal processors (DSPs), application specific integrated circuits (ASICs), field programmable gate arrays (FPGAs), or other programmable logic devices, transistor logic devices, hardware components, or any combination thereof. The general-purpose processor may be a microprocessor or any conventional processor.

[0205] The method steps in the embodiments of the present application can be implemented in a hardware manner or by a processor executing software instructions. The software instructions can be composed of corresponding software modules, and the software modules can be stored in a random access memory (RAM), flash memory, read-only memory (ROM), programmable ROM (PROM), erasable PROM (EPROM), electrically erasable PROM (EEPROM), register, hard disk, removable hard disk, CD-ROM, or any other form of storage medium well-known in the art. An exemplary storage medium is coupled to the processor so that the processor can read information from the storage medium and write information to the storage medium. Of course, the storage medium can also be a component of the processor. The processor and the storage medium can be located in an ASIC.

[0206] In the above embodiments, it can be implemented in whole or in part by software, hardware, firmware, or any combination thereof. When implemented using software, it can be implemented in whole or in part in the form of a computer program product. The computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, the processes or functions described in the embodiments of the present application are generated in whole or in part. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable systems. The computer instructions can be stored in a computer-readable storage medium or transmitted through the computer-readable storage medium. The computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center in a wired manner (such as coaxial cable, optical fiber, digital subscriber line (DSL)) or a wireless manner (such as infrared, wireless, microwave, etc.). The computer-readable storage medium can be any available medium that the computer can access or a data storage device such as a server or data center that includes one or more integrated available media. The available medium can be a magnetic medium (such as a floppy disk, hard disk, magnetic tape), an optical medium (such as a DVD), or a semiconductor medium (such as a solid state disk (SSD)), etc.

[0207] It can be understood that the various numerical numbers involved in the embodiments of the present application are only for the convenience of description and are not used to limit the scope of the embodiments of the present application.

[0208] Those skilled in the art can easily understand that the above description is only a preferred embodiment of the present application and is not intended to limit the present application. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present application shall be included within the protection scope of the present application.

Claims

1. A geophysical electromagnetic data inversion method based on entropy constraints, characterized in that: The following steps are involved: Step 1: After dividing the inversion space, set the conductivity value and build the initial inversion model; Step 2: Based on the initial inversion model, construct the magnetotelluric regularization inversion objective function based on entropy constraints; Step 3: Use the Gauss-Newton method to minimize the magnetotelluric regularization inversion objective function, update the initial inversion model, obtain the conductivity model, and realize magnetotelluric three-dimensional inversion; The objective function of the magnetotelluric regularization inversion in step 2 for: in, is the inversion model parameter, and the inversion model parameter is the logarithmic conductivity; is the data fitting term; is the model roughness term; is the minimum entropy term; , is the regularization parameter; , , For the The model parameters of the inversion unit, , is the number of inversion units, is a very small positive number with magnitude ; In step 2, the regularization parameter is calculated using the approximate spectral analysis method: in, is a random vector, In order to observe the sensitivity of magnetotelluric response data to the inversion model parameters, the adjoint forward modeling method is used to calculate. is the number of iterations, is a positive integer less than 3, is an empirical parameter, ; is the data weighting matrix; is the model roughness matrix; is the diagonal minimum entropy matrix for focusing inversion.

2. The geophysical electromagnetic method data inversion method according to claim 1, characterized in that: Step three specifically includes the following steps: Step 3.1: Use the Gauss-Newton method with quasi-quadratic convergence rate to minimize the magnetotelluric regularization inversion objective function, and use the Taylor expansion formula to obtain the normal equation; Step 3.2: Convert the normal equation into the least square form, use the least square algorithm to obtain the inversion model update amount, and use the linear search method to obtain the model update step size, and then obtain the prediction model; Step 3.3: When the prediction model meets the preset conditions, the prediction model is used as the conductivity model.

3. The geophysical electromagnetic data inversion method according to claim 2, characterized in that: Step 3.3 specifically includes the following steps: Step 3.3.1: Perform three-dimensional magnetotelluric forward modeling on the prediction model based on the vector finite element method to obtain magnetotelluric response data of the prediction model; Step 3.3.2: Use the normalized error to evaluate the fit between the predicted model magnetotelluric response data and the observed magnetotelluric response data. If the normalized fit error reaches the set threshold, the predicted model is used as the conductivity model; otherwise, go to step 3.3.3; Step 3.3.3: Determine whether the current iteration is the maximum inversion iteration number. If so, use the predicted model as the conductivity model, otherwise go to step 3.

2.

4. The geophysical electromagnetic data inversion method according to claim 1, characterized in that: Step 1 specifically includes the following steps: Step 1.1: Use unstructured tetrahedral grid to divide the inversion space; Step 1.2: In the unstructured tetrahedral mesh obtained in step 1.1, set the conductivity value and build the initial inversion model.

5. The geophysical electromagnetic method data inversion method according to claim 2, characterized in that: The normal equation is: in, is the model update amount; is the data weighting matrix; is the model roughness matrix; is the diagonal minimum entropy matrix for focusing inversion; , is the regularization parameter; is the sensitivity matrix; For the The data error between the predicted value and the observed value after iterations; For the Model parameters after iterations; is the prior model.

6. A geophysical electromagnetic data inversion system based on entropy constraints, characterized in that: include: The model building module sets the conductivity value after dividing the inversion space and constructs the initial inversion model; An objective function construction module is used to construct an entropy-constrained magnetotelluric regularized inversion objective function based on the initial inversion model; The model update module is used to minimize the magnetotelluric regularization inversion objective function using the Gauss-Newton method, update the initial inversion model, obtain the conductivity model, and realize the magnetotelluric three-dimensional inversion; Magnetotelluric regularization inversion objective function in the objective function building module for: in, is the inversion model parameter, and the inversion model parameter is the logarithmic conductivity; is the data fitting term; is the model roughness term; is the minimum entropy term; , is the regularization parameter; , , For the The model parameters of the inversion unit, , is the number of inversion units, is a very small positive number with magnitude ; The objective function building block uses an approximate spectral analysis method to calculate the regularization parameter: in, is a random vector, In order to observe the sensitivity of magnetotelluric response data to the inversion model parameters, the adjoint forward modeling method is used to calculate. is the number of iterations, is a positive integer less than 3, is an empirical parameter, ; is the data weighting matrix; is the model roughness matrix; is the diagonal minimum entropy matrix for focusing inversion.

Images