Boundary sensitivity enhanced gradient domain parameterized electromagnetic inversion method

Abstract

The application discloses a boundary sensitivity enhanced gradient domain parameterized electromagnetic inversion method and relates to the technical field of electromagnetic signal inversion, which comprises the following steps: carrying out grid division on a three-dimensional inversion region, establishing inversion model parameters, constructing three-direction difference operators based on the conductivity logarithmic difference value of adjacent grids in the inversion model parameters, and constructing a three-dimensional augmented reversible difference operator based on the three-direction difference operators; converting the inversion model parameters to a difference domain based on the three-dimensional augmented reversible difference operator and a physical scale transformation matrix, constructing a three-dimensional electromagnetic inversion objective function, and deriving a regularization equation by deriving the three-dimensional electromagnetic inversion objective function; the application converts the traditional inversion form taking the grid conductivity logarithmic value as a direct unknown quantity into an inversion form taking the augmented difference domain variable composed of the conductivity logarithmic difference between adjacent grids and a reference parameter as an unknown quantity, and provides innovative theoretical support and an implementable solution for the fine interpretation of geophysical electromagnetic data.

Description

Technical Field

[0001] This invention relates to the field of electromagnetic signal inversion technology, and in particular to a boundary-sensitive enhanced gradient domain parameterized electromagnetic inversion method. Background Technology

[0002] Geophysical inversion technology plays a crucial role in geological exploration and resource detection. By measuring physical parameters and optimizing models, inversion techniques can reveal fine subsurface structures. Electrical conductivity, as a key physical parameter reflecting the composition, porosity, and fluidity of the medium, has its spatial distribution directly related to the detection needs of mineral resources, oil and gas reservoirs, and potential engineering hazards. As exploration scenarios shift from regional surveys to detailed target area investigations, the required accuracy of electrical conductivity inversion has upgraded from macroscopic stratification to quantitative characterization of differences across different grids. This demand for refined grid-based electrical conductivity inversion has directly driven the technological leap from two-dimensional layered inversion to three-dimensional gridded inversion.

[0003] In traditional grid inversion, different electromagnetic methods use the differences in the sensitivity of observed data to grid conductivity to create constraints on grids of different depths and scales. For example, the magnetotelluric method (MT) utilizes natural electromagnetic fields, where low-frequency signals penetrate deep and high-frequency signals penetrate shallow, to create layered constraints on the conductivity of grids at different depths (shallow, medium, and deep). However, its lateral resolution is relatively low. Other methods include the DC resistivity method and the controlled-source audio-frequency magnetotelluric method (CSAMT). Traditional grid inversion suffers from a conflict between the need for finer grids and the issues of ambiguity and computational complexity. While increasing grid density improves accuracy, more grids lead to greater ambiguity and higher computational costs. Furthermore, "smoothing regularization," while suppressing ambiguity, obscures the true conductivity differences between grids. This core contradiction drove the development of this method. By transforming the inversion object from absolute conductivity to the conductivity difference between adjacent grids, it directly focuses on the conductivity variation characteristics between grids without increasing the number of grids, representing a targeted breakthrough in addressing the pain points of traditional inversion. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] Therefore, this invention provides a boundary-sensitive enhanced gradient domain parameterized electromagnetic inversion method to solve the problems of accuracy and computational complexity in underground structure boundary inversion.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution: This invention provides a boundary-sensitive enhanced gradient domain parameterized electromagnetic inversion method, comprising: meshing the three-dimensional inversion region, establishing inversion model parameters, constructing three-directional difference operators based on the logarithmic difference of conductivity between adjacent meshes in the inversion model parameters, and constructing a three-dimensional augmented invertible difference operator based on the three-directional difference operator; transforming the inversion model parameters to the difference domain based on the three-dimensional augmented invertible difference operator and the physical scale transformation matrix, constructing a three-dimensional electromagnetic inversion objective function, differentiating the three-dimensional electromagnetic inversion objective function to obtain the regularization equation; and performing forward modeling calculations based on the forward modeling theory equations. The conventional field and gradient field data are obtained. The Jacobian matrix of the model domain of the conventional field data for the absolute conductivity model is calculated. After transformation by the inverse matrix of the three-dimensional augmented invertible difference operator and the physical scale inverse transformation matrix, the conventional field Jacobian matrix of the physical difference domain is obtained. Then, the conventional field Jacobian matrix of the physical difference domain is linearly transformed to obtain the gradient field Jacobian matrix of the difference domain. Based on the conventional field Jacobian matrix of the physical difference domain and the gradient field Jacobian matrix of the difference domain, the regularization equation is solved by an iterative optimization algorithm to obtain the difference domain update. The current inversion model is updated according to the difference domain update to obtain the three-dimensional inversion result of the conductivity of the underground structure.

[0007] As a preferred embodiment of the boundary sensitivity enhanced gradient domain parameterized electromagnetic inversion method of the present invention, the steps of dividing the three-dimensional inversion region into grids, establishing inversion model parameters, and constructing three-directional difference operators based on the logarithmic difference of conductivity between adjacent grids in the inversion model parameters are as follows: The logarithmic conductivity values ​​of each grid in the three-dimensional inversion region are arranged in the order of z-direction, x-direction, and y-direction and then concatenated into a one-dimensional vector to obtain the inversion model parameters. A z-direction difference operator is constructed based on the logarithmic difference in conductivity of adjacent grids in the z-direction in the inversion model parameters; an x-direction difference operator is constructed based on the logarithmic difference in conductivity of adjacent grids in the x-direction in the inversion model parameters; and a y-direction difference operator is constructed based on the logarithmic difference in conductivity of adjacent grids in the y-direction in the inversion model parameters.

[0008] As a preferred embodiment of the boundary-sensitive enhanced gradient domain parameterized electromagnetic inversion method of the present invention, the specific steps for constructing the three-dimensional augmented invertible difference operator are as follows: Three directional difference operators are applied to the inversion model parameters to obtain three directional difference vectors. The z-direction difference vector, x-direction difference vector, and y-direction difference vector are integrated to obtain a three-dimensional universal difference vector. Based on the actual connectivity between grid nodes in the 3D inversion region grid, a reference node is selected for each connected component. Independent difference edges are selected using the spanning tree algorithm to form a spanning forest covering all grid nodes. The corresponding difference matrix is ​​concatenated with the selection matrix corresponding to the reference node to form an augmented difference operator matrix, thus obtaining the 3D augmented invertible difference operator.

[0009] As a preferred embodiment of the boundary sensitivity-enhanced gradient domain parameterized electromagnetic inversion method of the present invention, the specific steps for transforming the inversion model parameters to the difference domain and constructing the three-dimensional electromagnetic inversion objective function based on the three-dimensional augmented invertible difference operator are as follows: Based on the current value of the difference domain, combined with conventional field data, actual measured conventional electromagnetic data, predicted electromagnetic gradient data, actual measured electromagnetic gradient data, difference domain reference value, inversion model domain reference value, difference space weight matrix and inversion model domain weight matrix, a three-dimensional electromagnetic inversion objective function is constructed. The objective function for the three-dimensional electromagnetic inversion includes a conventional field data fitting term, a gradient field data fitting term, a difference domain constraint term, and an inversion model domain constraint term.

[0010] As a preferred embodiment of the boundary sensitivity enhanced gradient domain parameterized electromagnetic inversion method of the present invention, the difference space weight matrix is ​​obtained by weighting the difference parameters by the inverse of the grid step size, thereby converting the difference in the grid index domain into the rate of change in the physical space.

[0011] As a preferred embodiment of the boundary sensitivity enhanced gradient domain parameterized electromagnetic inversion method of the present invention, the regularization equation is obtained by linearizing the objective function with respect to the difference domain update quantity at the current inversion model, ignoring higher-order terms, taking the derivative and setting the derivative to zero, to obtain a linear regularization equation with respect to the difference domain update quantity.

[0012] As a preferred embodiment of the boundary sensitivity enhanced gradient domain parameterized electromagnetic inversion method of the present invention, the conventional field data and gradient field data are based on the forward modeling theoretical equations. The electromagnetic field solution problem corresponding to the current inversion model is transformed into a linear equation system using the finite element numerical discretization method. The electromagnetic field values ​​of the grid nodes under the current inversion model are obtained by solving the equation system. Then, the simulated observation data corresponding to the observation point position is extracted from the electromagnetic field values ​​of the grid nodes using the interpolation mapping matrix between the actual sampling points and the grid cells to obtain the conventional field data. Finally, the field space gradient at the observation point is calculated based on the difference of the conventional field data of adjacent observation points and the distance between adjacent observation points to obtain the gradient field data.

[0013] As a preferred embodiment of the boundary sensitivity enhanced gradient domain parameterized electromagnetic inversion method of the present invention, the difference domain gradient field Jacobian matrix is ​​obtained by differentiating the forward modeling theoretical equations to obtain the model domain Jacobian matrix of the conventional field data for the absolute conductivity model. Based on the chain rule, the model domain Jacobian matrix is ​​converted into the conventional field Jacobian matrix of the physical difference domain according to the inverse matrix of the three-dimensional augmented invertible difference operator and the physical scale inverse transformation matrix. The spatial gradient operator matrix is ​​constructed based on the inverse matrix of the normalized matrix of the measurement point spacing and the difference matrix of adjacent measurement points. The conventional field Jacobian matrix of the physical difference domain is obtained by linear transformation of the spatial gradient operator matrix.

[0014] As a preferred embodiment of the boundary sensitivity enhanced gradient domain parameterized electromagnetic inversion method of the present invention, the step of solving the regularization equation by using an iterative optimization algorithm refers to solving the regularization equation by using the conjugate gradient method to obtain the difference domain update amount, and using the backtracking method to obtain the current inversion model update step size.

[0015] As a preferred embodiment of the boundary sensitivity enhanced gradient domain parameterized electromagnetic inversion method of the present invention, the specific steps for updating the current inversion model according to the differential domain update amount to obtain the three-dimensional inversion result of the conductivity of the underground structure are as follows: The current inversion model is updated based on the difference domain update amount and the current inversion model update step size; Based on the updated current inversion model, the regularization parameter is dynamically adjusted in a cooling manner according to the ratio of the current data fit difference to the target fit difference. This balances the data fit and model constraints in subsequent iterations, and the iteration is repeated until the data fit difference reaches the expected level, thus obtaining the three-dimensional inversion result of the electrical conductivity of the underground structure.

[0016] The beneficial effects of this invention are as follows: by using the augmented difference domain variables, which are composed of the relative differences between inversion grids and a small number of reference parameters, as the unknowns to be solved in the three-dimensional inversion objective function, and dynamically adjusting the weights of the residual terms of electromagnetic field data and gradient field data, this invention provides theoretical support and an implementable solution for the refined interpretation of geophysical electromagnetic data. Attached Figure Description

[0017] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. 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.

[0018] Figure 1 This is a flowchart of the boundary-sensitive enhanced gradient domain parameterized electromagnetic inversion method.

[0019] Figure 2 A schematic diagram illustrating the construction of an invertible difference transformation matrix under a one-dimensional mesh. Detailed Implementation

[0020] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0021] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0022] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0023] Reference Figures 1-2 As one embodiment of the present invention, this embodiment provides a boundary-sensitive enhanced gradient domain parameterized electromagnetic inversion method, comprising the following steps: S1. Mesh the three-dimensional inversion region, establish inversion model parameters, construct three-directional difference operators based on the logarithmic difference of conductivity between adjacent grids in the inversion model parameters, and construct a three-dimensional augmented invertible difference operator based on the three-directional difference operators.

[0024] S1.1. The logarithmic values ​​of the conductivity of each grid in the three-dimensional inversion region are then processed according to the following order: Direction, again Direction, final The directions are arranged in order and concatenated into a one-dimensional vector to obtain the inversion model parameters.

[0025] Specifically, the parameters of the three-dimensional inversion model are the set of logarithmic values ​​of the grid conductivity corresponding to each measurement point, assuming... Direction One grid, Direction One grid, Direction Each grid.

[0026] Total number of grids in the 3D inversion region for: ; The logarithmic conductivity of each grid cell is denoted as: ,in: ; ; ; One-dimensional global index of the grid for: ; Following the order of z, then x, and finally y, the logarithmic values ​​of the conductivity of all grids are concatenated into a one-dimensional vector. That is, for each fixed... ,according to Arrange, then Loop, then cycle, = , And so on, until all of them are... After the logarithmic conductivity values ​​of each grid are arranged, at this point: ; The resulting one-dimensional vector is: ; ; in, express Direction Each measuring point The j-th measuring point in the direction The inversion conductivity parameter at the k-th grid in the direction, ; .

[0027] S1.2. Construct a z-direction difference operator based on the logarithmic difference of conductivity in the z-direction between adjacent grids in the inversion model parameters, construct an x-direction difference operator based on the logarithmic difference of conductivity in the x-direction between adjacent grids in the inversion model parameters, and construct a y-direction difference operator based on the logarithmic difference of conductivity in the y-direction between adjacent grids in the inversion model parameters.

[0028] Specifically, Directional difference focuses on the same Adjacent in the plane The logarithmic difference of the conductivity of the grid is defined as follows: Directional difference vector: ; in, This represents the parameter difference between adjacent meshes along the z-direction; Indicates the parameters of the inversion model; Indicates the grid index in the x-direction; Indicates the grid index in the y-direction; express Grid index of direction; Indicates the first indivual To position, number indivual To position, number indivual Invert the model parameter values ​​at the grid; Indicates the first indivual To position, number indivual To position, number indivual The inversion model parameter values ​​at the grid.

[0029] ; ; ; Directional difference operator for A dimensional matrix that satisfies .

[0030] for Corresponding index The element is located in the matrix at row number: ; in, express Linear row number of the directional difference term; only in the corresponding The column position is set to -1, in the corresponding The column position is set to 1, and the rest of the columns are 0.

[0031] Build Directional difference operator ; Directional difference focuses on the same The logarithmic difference in conductivity between adjacent x-grids in a plane is defined as follows: Directional difference vector: ; in, express The difference in parameters between adjacent grids in the direction; Indicates the first indivual To the grid, the first indivual To the grid, the first indivual The inversion model parameter values ​​at the grid.

[0032] ; ; ; in, Directional difference operator for A dimensional matrix that satisfies .

[0033] for Corresponding index The element, whose row number in the matrix is: ; It should be noted that only in the corresponding The column position is set to -1, in the corresponding The column position is set to 1, and the rest of the columns are 0.

[0034] Build Directional difference operator ; Directional difference focuses on the same Adjacent in the plane The logarithmic difference of the conductivity of the grid is defined as follows: Directional difference vector: ; in, express The difference in parameters between adjacent grids in the direction.

[0035] ; ; ; in, Directional difference operator for A dimensional matrix that satisfies .

[0036] for Corresponding index The element, whose row number in the matrix is: ; It should be noted that, express The linear row number corresponding to the direction difference term is only present in the corresponding row. The column position is set to -1, in the corresponding The column position is set to 1, and the rest of the columns are 0.

[0037] S1.3. Apply the three directional difference operators to the inversion model parameters to obtain three directional difference vectors. Integrate the z-direction difference vector, x-direction difference vector, and y-direction difference vector to obtain a three-dimensional universal difference vector.

[0038] Specifically, by integrating the difference vectors in the three directions—z, x, and finally y—a three-dimensional universal difference vector is obtained. The three-dimensional universal difference vector can be represented as: ; in, Represents a three-dimensional universal difference vector; This represents the transpose operation of a vector.

[0039] S1.4. Based on the actual connectivity between the grid nodes in the 3D inversion region grid, a reference node is selected for each connected component. The spanning tree algorithm is used to select independent difference edges to form a spanning forest covering all grid nodes. The corresponding difference matrix is ​​concatenated with the selection matrix corresponding to the reference node to form an augmented difference operator matrix, thus obtaining the 3D augmented invertible difference operator.

[0040] To ensure the invertibility and full rank of the difference operator, simply stacking all directional differences is not feasible. In the M nodes of the 3D mesh, B reference values ​​are selected. When the 3D inversion region mesh has multiple connected components not connected by difference edges, a reference node is selected for each connected component, constructing a B×M reference selection matrix. Each row corresponds to a reference, and the corresponding column of the node is 1; then select M-B independent edges from the difference edges to form a spanning tree covering all nodes.

[0041] The following is the simplest case, where the entire model is connected and B=1: Take all z-direction difference edges, and connect each edge to adjacent z-grids at the same (x, y) position: ; Then take the x-direction difference edge on the top layer z=1 and connect adjacent x-grids in the same row: ; Then take the y-direction difference edge at x=1 in the first column of the top layer and connect adjacent y-grids in the same column: ; Total number of difference sides: ; Add a baseline value The total number of rows is M, forming a square matrix.

[0042] Let the selected set of difference edges be . The corresponding (M−1)×M difference matrix is 1×M reference selection matrix The augmented difference operator is then: ; For the general case of multiple connected components, the grid nodes are considered as a graph, and the candidate difference edges of the directed difference between each adjacent grid are considered as edges of the graph; the graph-theoretic spanning tree algorithm is used to select M−B edges, resulting in... Choose a reference node for each connected component and construct a matrix. Each row corresponds to a baseline, and the corresponding column at that node is 1, thus obtaining... It is a square formation and has a full rank.

[0043] In summary, the conventional inversion object Forward mapping to the new inverse object This achieves the transformation of the original model parameters to the difference domain: ; Conversely, it can be seen from Reverse recovery of original model parameters : ; Define the generalized inverse matrix of the augmented difference operator. for: ; but: ; S1.5. Define a scaling transformation matrix from index difference to physical difference, so that the inversion constraints are based on the real physical rate of change, and transform the inversion parameters of the difference domain from the grid index space to the physical space.

[0044] Specifically, let the grid step sizes along the x, y, and z directions in the three-dimensional space be respectively... ,definition The direction is a discrete interval reciprocal diagonal matrix, and the variation in that direction is normalized using the grid step size. for identity matrix of order: ; Directional discrete interval reciprocal matrix: ; Directional discrete interval reciprocal matrix: ; Three-dimensional discrete interval matrix The scaling transformation matrix is ​​constructed based on the grid spacing, converting the index space difference into the rate of change in physical space. Simultaneously, to avoid unnecessary scaling transformations of the reference parameters... Using a block diagonal format: ; in, The identity matrix acting on the reference parameters; block dimension and augmented difference vector The mid-base parameters and the arrangement of the selected independent difference edges are consistent.

[0045] Define the physical scale inverse transformation matrix as: ; in, This is the physical scale transformation matrix, used to convert the difference between adjacent grids in the index space into the rate of change of conductivity in the physical space; It is its inverse matrix, used to convert the rate of change of physical space back to the index space difference component.

[0046] therefore: ; ; make: ; Then we have: ; in, These are augmented difference variables in physical space, including baseline model parameters and independent difference parameters normalized to physical scale.

[0047] S2. Based on the three-dimensional augmented invertible difference operator and the physical scale transformation matrix, the inversion model parameters are transformed into the difference domain to construct the three-dimensional electromagnetic inversion objective function. The regularization equation is obtained by differentiating the three-dimensional electromagnetic inversion objective function.

[0048] S2.1. Based on the current value of the difference domain, combined with conventional field data, actual measured conventional electromagnetic data, predicted electromagnetic gradient data, actual measured electromagnetic gradient data, difference domain reference value, inversion model domain reference value, difference space weight matrix and inversion model domain weight matrix, a three-dimensional electromagnetic inversion objective function is constructed.

[0049] ; Specifically, the objective function is constructed, and in the k-th iteration, the unknown to be solved is the model update amount. The objective function includes data fitting terms and model constraint terms, and is optimized by adjusting the difference domain update amount. To achieve a balance between fitting observational data and ensuring the physical plausibility of the model, the objective function for the three-dimensional electromagnetic inversion is constructed as follows: in, For fitting terms of conventional field data; This is the fitting term for the gradient field data. Because the abrupt change in conductivity at the boundary will cause the gradient value to increase significantly, the gradient data is more sensitive to the boundary of the underground structure. These are the difference domain constraint terms; This is a constraint term in the inversion model domain, based on prior geological information. If the reference electrical conductivity value of a stratum at a certain depth is known... This term forces the current model to... Approximating the reference value reduces ambiguity. If there is no reliable prior information, one can take... Only the difference domain constraints are retained to avoid prior errors misleading the inversion results; by adjusting The weights are used to select the inversion field when When the objective function only contains a gradient field data fitting term, only the gradient field data is inverted; when When the objective function only includes a fitting term for conventional field data, it only inverts conventional field data; when At that time, the joint inversion of the conventional field and the gradient field is realized; These are the internal balance coefficients for the model constraint terms. .

[0050] Among them, the conventional field data fitting term and the gradient field data fitting term constitute the data fitting term, which forces the conventional field data to approximate the observed data; the difference domain constraint term and the inversion model domain constraint term constitute the model constraint term, which, through regularization parameters, [is used to]... Control its overall strength to ensure that the inversion results conform to physical laws.

[0051] It should be noted that the objective function for three-dimensional electromagnetic inversion includes a conventional field data fitting term, a gradient field data fitting term, a difference domain constraint term, and an inversion model domain constraint term.

[0052] S2.2. At the current inversion model, the objective function is linearized first-order with respect to the difference domain update quantity. After ignoring higher-order terms, the derivative is taken and set to zero to obtain the linear regularization equation with respect to the difference domain update quantity.

[0053] Let the current inversion model in the k-th iteration be... The corresponding difference domain parameters are The update amount of the difference domain to be solved is Perform a first-order Taylor expansion of the predicted data terms in the objective function on the current model: ; ; in, For the conventional field obtained based on the forward model, This is the gradient field prediction data obtained based on the forward model of the current model; , This is the Jacobian matrix of the model domain for the absolute conductivity model based on conventional field data.

[0054] Will Substituting the values, we obtain a linearized expression with the difference domain update as the variable: ; ; The difference field Jacobian satisfies: ; ; Define the residual vector of the current regular field data and the residual vector of the gradient field data in the k-th iteration as follows: ; ; in, This is conventional field data, which is the conventional electromagnetic data obtained by the current inversion model through forward modeling; These are standard electromagnetic data measured in practice. The electromagnetic gradient data is obtained by using a spatial difference operator to obtain gradient-type data from the conventional field data corresponding to the current model. The electromagnetic gradient data (the rate of change of field strength along the x / y / z directions) is actually measured and is more sensitive to the boundaries of underground structures, thus improving the inversion resolution. Substituting the first-order linearized relation, the update amount of the difference domain is... The linearized residual function is: ; ; The current residual vector is given by the difference between the current forward response and the observed data. The Jacobian matrix is ​​used to establish the first-order linear relationship between the residual and the update quantity of the difference domain, and also participates in the construction of the coefficient matrix and right-hand side terms of the normal equation.

[0055] Substituting the linearized relationship into the overall objective function of the three-dimensional electromagnetic inversion Update the amount using the difference domain To optimize the variables, the objective function is defined as: ; in, For regular field data weight matrix, This represents the gradient field data weight matrix, with dimensions of [dimensions missing]. , Matching, with the diagonal elements representing the weights of each observation, reduces the interference of noisy data on the inversion results; This is the difference space weight matrix, which can be set based on prior information, depth weights, or structural constraints. If there are no additional priors, it can be set to... ; This is the current value of the difference field at the k-th iteration; This is the reference value for the differential domain, usually set to 0, representing the minimum structural constraint. This means that the underground electrical conductivity changes gradually, which is consistent with the assumption of geological continuity. Large changes are only allowed when the observation data supports it. This is the domain weight matrix for the inversion model. If it is an identity matrix, it constrains the absolute conductivity values; if it is a smoothing matrix, it constrains the continuity of conductivity between adjacent grids to avoid physical meaningless abrupt changes in the model. Take the derivative of the objective function and set the first derivative to zero: ; neglect The term yields the regularized normal equation: ; ; ; Wherein, the coefficient matrix The approximate Hessian matrix of the linearized inversion subproblem in the kth iteration is composed of the conventional field data term, the gradient field data term, the difference domain constraint term, and the inversion model domain constraint term. The Jacobian matrices of the conventional field and gradient field respectively characterize the sensitivity of different observation data to the parameters of the physical difference domain model. The data weight matrix is ​​used to adjust the contribution of different data components in the inversion. The difference domain constraint term and the inversion model domain constraint term are used to constrain the spatial rate of change of conductivity and the absolute conductivity model, respectively, so that the inversion update can maintain physical rationality and numerical stability while satisfying the data fitting. Right end item The update amount of the difference domain model in the k-th iteration is determined by the current data residuals and the model reference term. The driving direction.

[0056] By solving the above equation, we obtain... This leads to the updated model for this iteration.

[0057] S3. Based on the forward modeling equations, perform forward modeling calculations on the current inversion model to obtain the conventional field and gradient field data. Calculate the Jacobian matrix of the model domain for the conventional field data relative to the absolute conductivity model. Transform the model domain using the inverse matrix of the three-dimensional augmented invertible difference operator and the inverse physical scale transformation matrix to obtain the Jacobian matrix of the conventional field in the physical difference domain. Then, perform a linear transformation on the Jacobian matrix of the conventional field in the physical difference domain to obtain the Jacobian matrix of the gradient field in the difference domain.

[0058] S3.1. Based on the forward modeling theory equations, the electromagnetic field problem corresponding to the current inversion model is transformed into a linear equation system using the finite element numerical discretization method. The electromagnetic field values ​​of the grid nodes under the current inversion model are obtained by solving the equations. The simulated observation data corresponding to the observation point position is extracted from the electromagnetic field values ​​of the grid nodes using the interpolation mapping matrix between the actual sampling points and the grid elements to obtain the conventional field data. Then, based on the difference of the conventional field data of adjacent observation points and the distance between adjacent observation points, the field space gradient at the observation point is calculated to obtain the gradient field data.

[0059] Specifically, based on the forward modeling equations, and through numerical discretization methods such as the finite element method, they can be transformed into corresponding linear equations: ; in, This is a coefficient matrix, whose non-zero elements are determined by the grid conductivity. With grid step size The decision is made based on the inversion model parameters of the k-th iteration. Directly related; It is solved by forward modeling. The electromagnetic field value at point b; b is the source term on the right-hand side of the forward modeling equation, which is determined by the source conditions and boundary conditions and is independent of the model parameters.

[0060] The expression for extracting conventional field data is: ; in, The simulated observation data for the k-th iteration; It is the electromagnetic field value vector on the grid node; It is the interpolation mapping matrix between actual sampling points and grid cells, derived from electromagnetic field quantities. The observable data is extracted, the electric field of the grid nodes is converted into observation point data, and after obtaining the simulated observation data of the kth iteration, the simulated observation data at adjacent observation points are spatially differentiated and normalized according to the distance between adjacent observation points to obtain the gradient field data of the kth iteration.

[0061] S3.2. Differentiate the forward modeling equations to obtain the model domain Jacobian matrix of the absolute conductivity model for the conventional field data. Based on the chain rule, and according to the inverse matrix of the three-dimensional augmented invertible difference operator and the inverse transformation matrix of the physical scale, the model domain Jacobian matrix is ​​converted into the conventional field Jacobian matrix of the physical difference domain. Based on the inverse matrix of the normalized matrix of the measurement point spacing and the difference matrix of adjacent measurement points, construct the spatial gradient operator matrix. Perform a linear transformation on the conventional field Jacobian matrix of the physical difference domain through the spatial gradient operator matrix to obtain the gradient field Jacobian matrix of the difference domain.

[0062] Specifically, the process of deriving the Jacobian matrix while maintaining the simplicity of mathematical notation is as follows: Both sides Differentiate: ; The derivatives of the field solution with respect to the model parameters are obtained by rearranging: ; Combination ,right Differentiation yields: ; in, This represents the Jacobian matrix of the model domain for the absolute conductivity model based on conventional field data.

[0063] Substitution The expression, finally the first The Jacobian matrix for the regular field data in the next iteration is: ; in, The Jacobian matrix of the conventional field data for the k-th iteration relative to the absolute conductivity model only reflects the absolute conductivity. Observational data The influence of conductivity difference, however, is not considered in this method. Impact on the data.

[0064] According to the chain rule of derivatives: ; in, Let be the Jacobian matrix of the physical difference field in the k-th iteration. It is a three-dimensional universal difference vector in physical space, obtained by weighted sum of conductivity differences in the z, x, and y directions; The sensitivity of the absolute model to the difference model is given by the inverse matrix of the augmented invertible difference operator. Inverse transformation matrix of physical scale A joint decision.

[0065] After accepting the difference vector in physical space: ; in, This is the total inverse transformation matrix.

[0066] ; The difference domain gradient field Jacobian matrix is ​​obtained through a linear transformation of the field data Jacobian matrix; the electromagnetic gradient data reflects the rate of change of the field strength in space, and the expression for the spatial gradient of the field along the measuring point is: ; in, It is the data difference between adjacent measuring points; It is the distance between adjacent measuring points.

[0067] According to the basic definition, the Jacobian matrix of the gradient field in the difference domain can also be obtained using the chain rule: ; ; Among them, the spatial gradient operator matrix From the spacing normalization matrix sum difference matrix composition: ; Therefore, the first The Jacobian matrix of the gradient field in the difference domain of the next iteration is expressed as: ; in, Let be the Jacobian matrix of the difference field gradient field in the k-th iteration; This is the inverse of the normalized matrix of the distance between measurement points; The difference matrix between adjacent measurement points. The Jacobian matrix of the model domain for the absolute conductivity model is given by the standard field data. This is the total inverse transformation matrix; this matrix correlates the sensitivity of the gradient data with the parameters of the difference domain model, providing constraints on the gradient data for subsequent inversion iterations.

[0068] S4. Based on the Jacobian matrix of the conventional field and the Jacobian matrix of the gradient field in the physical difference domain, the regularization equation is solved by an iterative optimization algorithm to obtain the difference domain update. The current inversion model is updated according to the difference domain update to obtain the three-dimensional inversion result of the conductivity of the underground structure.

[0069] S4.1. Solving the regularization equation using an iterative optimization algorithm refers to solving the regularization equation using the conjugate gradient method to obtain the difference domain update amount, and using the backtracking method to obtain the current inversion model update step size.

[0070] Specifically, PCG is used to solve linear systems: ; Wherein, the coefficient matrix Let be the Hessian matrix of the k-th iteration. For the right-hand driving term of the k-th iteration, That is, the direction of correction of the conductivity difference in this iteration is the solution of the inversion equation.

[0071] When using PCG iterative solutions, there is no need to directly assemble large sparse matrices; iterative calculations are completed only through matrix-vector multiplication, reducing memory usage and computational complexity, and ultimately obtaining the optimal correction amount for the model difference. ; get Then, adopt Backtracking method determines the current inversion model update step size Find a suitable step size to ensure that the objective function truly decreases.

[0072] It should be noted that the objective function is defined with respect to the step size α as follows: ; in, Indicates the update quantity of the difference domain. hour, The overall objective function value obtained from the inversion includes the data fitting term and the model regularization term.

[0073] use The criterion requires the objective function to satisfy: ; in, Before the update The objective function value; For the objective function in gradient at, Ensure the update direction is downward; The Armijo parameter controls the extent of the step size decrease, ensuring that the objective function continues to decrease during iteration and avoiding iteration stagnation caused by an excessively small step size.

[0074] In addition to ensuring the objective function decreases according to the Armijo criterion, there is an additional constraint that the data fit difference does not exceed the target threshold: ; in, , Step size The corresponding data fit between the conventional field and the gradient field is poor; These are the gradient field weighting coefficients, used to balance the fitting contributions of the conventional field and the gradient data; The target fitting difference for the k-th iteration is the maximum fitting difference threshold allowed in this iteration. This avoids underfitting of the data due to excessive regularization or rebound of the fitting difference due to excessive step size, thus ensuring the fitting accuracy of the inversion.

[0075] S4.2. Update the current inversion model based on the difference domain update amount and the current inversion model update step size.

[0076] Specifically, the step size that satisfies the dual constraints is found using the above method. The update is first performed in the difference domain, and then the model is restored to the absolute conductivity model through an invertible transformation to ensure the stability of the difference domain inversion. ; ; Among them, it means The difference domain model parameters for the (k+1)th iteration; The absolute conductivity model parameters for the (k+1)th iteration are represented by the inverse transformation matrix. Update the differential domain The update, converted to an absolute model, is then superimposed onto the current model. superior.

[0077] S4.3. Based on the updated current inversion model, the regularization parameter is dynamically adjusted in a cooling manner according to the ratio of the current data fit difference to the target fit difference. The data fit and model constraints of subsequent iterations are balanced, and the iteration is repeated until the data fit difference reaches the expected level, so as to obtain the three-dimensional inversion result of the electrical conductivity of the underground structure.

[0078] Specifically, based on the ratio of the current data fit difference to the target fit difference, the regularization parameter for the next iteration is adjusted so that the inversion process can dynamically balance data fit and model constraints. ; ; in, The ratio of the fit difference is used to determine the current fit state of the inversion, and it gradually approaches 1; The actual data fit difference for the current iteration; The target fitting difference for the k-th iteration is usually gradually reduced with the number of iterations, guiding the inversion to converge gradually. is the regularization parameter for the k-th iteration, used to balance the weights of the data fitting term and the model regularization term; The adjusted regularization parameters for the next iteration gradually stabilize, resulting in a cooling effect.

[0079] Reference Figure 2 , Figure 2 The difference matrix to be constructed A schematic diagram of generating a simple one-dimensional mesh: This diagram illustrates the difference relationship under a one-dimensional mesh model, showing how, in a single dimension, the difference in conductivity between adjacent meshes can be used to generate a simple one-dimensional mesh. and benchmark value Constructing an invertible transformation .

[0080] In summary, this invention provides theoretical support and a feasible solution for the refined interpretation of geophysical electromagnetic data by using the augmented difference domain variables, which consist of the relative differences between inversion grids and a small number of reference parameters, as the unknowns to be solved in the three-dimensional inversion objective function, and by dynamically adjusting the weights of the residual terms of electromagnetic field data and gradient field data.

[0081] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A boundary-sensitive enhanced gradient domain parameterized electromagnetic inversion method, characterized in that, include: The three-dimensional inversion region is meshed, inversion model parameters are established, and three-directional difference operators are constructed based on the logarithmic difference of conductivity between adjacent meshes in the inversion model parameters. A three-dimensional augmented invertible difference operator is then constructed based on the three-directional difference operators. Based on the three-dimensional augmented invertible difference operator and the physical scale transformation matrix, the inversion model parameters are transformed into the difference domain, a three-dimensional electromagnetic inversion objective function is constructed, and the regularization equation is obtained by differentiating the three-dimensional electromagnetic inversion objective function. Based on the forward modeling equations, forward modeling calculations are performed on the current inversion model to obtain conventional field and gradient field data. The model domain Jacobian matrix of the conventional field data for the absolute conductivity model is calculated. After transformation by the inverse matrix of the three-dimensional augmented invertible difference operator and the physical scale inverse transformation matrix, the conventional field Jacobian matrix of the physical difference domain is obtained. Then, a linear transformation is performed on the conventional field Jacobian matrix of the physical difference domain to obtain the gradient field Jacobian matrix of the difference domain. Based on the Jacobian matrix of the conventional field and the Jacobian matrix of the gradient field in the physical difference domain, an iterative optimization algorithm is used to solve the regularization equation to obtain the difference domain update. The current inversion model is then updated based on the difference domain update to obtain the three-dimensional inversion result of the conductivity of the underground structure.

2. The boundary sensitivity enhanced gradient domain parameterized electromagnetic inversion method as described in claim 1, characterized in that, The process of dividing the three-dimensional inversion region into meshes, establishing inversion model parameters, and constructing three-directional difference operators based on the logarithmic difference in conductivity between adjacent meshes in the inversion model parameters is as follows: The logarithmic conductivity values ​​of each grid in the three-dimensional inversion region are arranged in the order of z-direction, x-direction, and y-direction and then concatenated into a one-dimensional vector to obtain the inversion model parameters. A z-direction difference operator is constructed based on the logarithmic difference in conductivity of adjacent grids in the z-direction in the inversion model parameters; an x-direction difference operator is constructed based on the logarithmic difference in conductivity of adjacent grids in the x-direction in the inversion model parameters; and a y-direction difference operator is constructed based on the logarithmic difference in conductivity of adjacent grids in the y-direction in the inversion model parameters.

3. The boundary sensitivity enhanced gradient domain parameterized electromagnetic inversion method as described in claim 1, characterized in that, The specific steps for constructing the three-dimensional augmented invertible difference operator are as follows: Three directional difference operators are applied to the inversion model parameters to obtain three directional difference vectors. The z-direction difference vector, x-direction difference vector, and y-direction difference vector are integrated to obtain a three-dimensional universal difference vector. Based on the actual connectivity between grid nodes in the 3D inversion region grid, a reference node is selected for each connected component. Independent difference edges are selected using the spanning tree algorithm to form a spanning forest covering all grid nodes. The corresponding difference matrix is ​​concatenated with the selection matrix corresponding to the reference node to form an augmented difference operator matrix, thus obtaining the 3D augmented invertible difference operator.

4. The boundary-sensitivity enhanced gradient domain parameterized electromagnetic inversion method as described in claim 1, characterized in that, The method for transforming the inversion model parameters to the difference domain based on the three-dimensional augmented invertible difference operator and constructing the three-dimensional electromagnetic inversion objective function involves the following steps: Based on the current value of the difference domain, combined with conventional field data, actual measured conventional electromagnetic data, predicted electromagnetic gradient data, actual measured electromagnetic gradient data, difference domain reference value, inversion model domain reference value, difference space weight matrix and inversion model domain weight matrix, a three-dimensional electromagnetic inversion objective function is constructed. The objective function for the three-dimensional electromagnetic inversion includes a conventional field data fitting term, a gradient field data fitting term, a difference domain constraint term, and an inversion model domain constraint term.

5. The boundary-sensitivity enhanced gradient domain parameterized electromagnetic inversion method as described in claim 4, characterized in that, The difference space weight matrix is ​​obtained by weighting the difference parameters with the reciprocal of the grid step size, and converting the difference in the grid index domain into the rate of change in the physical space.

6. The boundary-sensitivity enhanced gradient domain parameterized electromagnetic inversion method as described in claim 1, characterized in that, The regularization equation is obtained by linearizing the objective function with respect to the difference domain update quantity at the current inversion model, ignoring higher-order terms, taking the derivative and setting the derivative to zero, and obtaining the linear regularization equation with respect to the difference domain update quantity.

7. The boundary-sensitivity enhanced gradient domain parameterized electromagnetic inversion method as described in claim 1, characterized in that, The conventional field data and gradient field data are based on the forward modeling theory equations. The electromagnetic field problem corresponding to the current inversion model is transformed into a linear equation system using the finite element numerical discretization method. The electromagnetic field values ​​of the grid nodes under the current inversion model are obtained by solving the equation system. Then, the simulated observation data corresponding to the observation point position is extracted from the electromagnetic field values ​​of the grid nodes using the interpolation mapping matrix between the actual sampling points and the grid cells to obtain the conventional field data. Finally, the field space gradient at the observation point is calculated based on the difference of the conventional field data of adjacent observation points and the distance between adjacent observation points to obtain the gradient field data.

8. The boundary sensitivity enhanced gradient domain parameterized electromagnetic inversion method as described in claim 1, characterized in that, The Jacobian matrix of the gradient field in the difference domain is obtained by differentiating the forward modeling equations to obtain the model domain Jacobian matrix of the absolute conductivity model for the conventional field data. Based on the chain rule, the model domain Jacobian matrix is ​​transformed into the conventional field Jacobian matrix of the physical difference domain according to the inverse matrix of the three-dimensional augmented invertible difference operator and the inverse transformation matrix of the physical scale. The spatial gradient operator matrix is ​​constructed based on the inverse matrix of the normalized matrix of the measurement point spacing and the difference matrix of adjacent measurement points. The Jacobian matrix of the conventional field of the physical difference domain is obtained by linearly transforming the spatial gradient operator matrix.

9. The boundary-sensitivity enhanced gradient domain parameterized electromagnetic inversion method as described in claim 1, characterized in that, The method of using iterative optimization algorithm to solve the regularization equation refers to using the conjugate gradient method to solve the regularization equation, obtaining the difference domain update amount, and using the backtracking method to obtain the current inversion model update step size.

10. The boundary-sensitivity enhanced gradient domain parameterized electromagnetic inversion method as described in claim 1, characterized in that, The steps for updating the current inversion model based on the differential domain update amount to obtain the three-dimensional inversion result of the electrical conductivity of the underground structure are as follows: The current inversion model is updated based on the difference domain update amount and the current inversion model update step size; Based on the updated current inversion model, the regularization parameter is dynamically adjusted in a cooling manner according to the ratio of the current data fit difference to the target fit difference. This balances the data fit and model constraints in subsequent iterations, and the iteration is repeated until the data fit difference reaches the expected level, thus obtaining the three-dimensional inversion result of the electrical conductivity of the underground structure.

Images