Fourier series-based one-dimensional magnetotelluric anisotropy inversion method and system

Abstract

The invention discloses a Fourier series-based one-dimensional magnetotelluric anisotropy inversion method and system, and the method is used for achieving the inversion calculation of the resistivity of each layer of underground medium in an underground layer model, greatly reducing the number of inversion variables, and remarkably improving the calculation efficiency and reliability. The core of the method is that a Fourier series expansion method is introduced, and the Fourier series expansion method is creatively adopted to approximate resistivity distribution of an underground medium; in addition, the objective function is basically optimized in multiple aspects, for example, a Tikhonov regularization objective function based on an L1 norm is introduced; an alternating direction multiplier method (ADMM) algorithm is introduced and is used for efficiently optimizing the constructed target function; and the upper and lower bounds of the resistivity can be set according to empirical data, so that the search range of a solution space is narrowed, the calculation efficiency is improved, a complicated underground structure can be dealt with better, and a physical unreasonable result caused by an overlarge solution space is avoided.

Description

Technical Field

[0001] The invention belongs to the technical field of one-dimensional magnetotelluric anisotropy inversion, and in particular relates to a one-dimensional magnetotelluric anisotropy inversion method and system based on Fourier series. Background Art

[0002] At present, one-dimensional magnetotelluric (MT) anisotropic inversion mainly relies on the traditional Gauss-Newton algorithm based on the L2 norm. This method estimates the resistivity of each layer by assuming that the underground medium is a layered uniform structure. However, when dealing with models with drastic resistivity changes, the L2 norm-based regularization method may not be able to effectively capture the drastic changes in the underground structure, resulting in insufficient accuracy of the imaging results. In particular, the traditional inversion method needs to process a large number of inversion variables, which easily leads to low computational efficiency, making it difficult to effectively handle complex underground structures, and the solution space is too large, which may produce physically unreasonable results. Summary of the invention

[0003] The purpose of the present invention is to solve the problem that the above-mentioned prior art has low computational efficiency due to a large number of inversion variables, and the accuracy and robustness of the inversion results need to be improved, and further provide a one-dimensional magnetotelluric anisotropy inversion method and system based on Fourier series. Among them, the technical solution of the present invention innovatively proposes to use the Fourier series expansion method to approximate the resistivity distribution of the underground medium, that is, introduces the Fourier series expansion method, and expresses the resistivity of each layer of the underground medium in the underground layer model as a function of the depth of each layer, and then converts the inversion parameters of the traditional resistivity model into coefficients expressed using Fourier series expansion, which greatly reduces the inversion parameter variables, improves the computational efficiency, and can better cope with complex underground structures, avoiding physically unreasonable results due to excessive solution space.

[0004] To this end, the present invention provides the following technical solutions:

[0005] On the one hand, the present invention provides a one-dimensional magnetotelluric anisotropy inversion method based on Fourier series, comprising the following steps:

[0006] S1: Obtain the depth information of the underground layer model in the inversion area;

[0007] S2: A Fourier series expansion method is introduced to express the inversion target parameters of each layer of underground medium in the underground layer model as a function of the depth of each layer, that is, the inversion target parameters are approximated by Fourier series, and the original inversion parameters composed of the inversion target parameters of each layer are converted into inversion parameters constructed by Fourier coefficients;

[0008] S3: Set the inversion constraint function to perform inversion iteration operation to obtain the resistivity of each layer of underground medium in the inversion area;

[0009] Among them, each iterative operation determines the inversion target parameters of each layer of underground media based on the function constructed in step S2, the inversion parameters and the depth information; then, a forward operation is performed based on the inversion target parameters of each layer of underground media to obtain the forward response of the inversion area; and then, based on the inversion constraint function, the forward response and the forward response observation value obtained by the equipment observation in the inversion area, the operation is performed to update the inversion parameters.

[0010] Taking resistivity as an example, in traditional magnetotelluric inversion, it is usually necessary to use a layered model to describe the changes in underground resistivity, which results in a large number of variables that need to be inverted during the inversion process and high computational complexity. For example, the underground layer model is divided into 120 layers, each layer has resistivity in two directions. Then the inversion variables are 240; when the Fourier coefficients proposed in the technical solution of the present invention are used for approximation, the inversion parameters become Fourier coefficients, which greatly reduces the number of parameters of the inversion variables. Therefore, the technical solution of the present invention innovatively proposes to use the Fourier series expansion method to approximate the resistivity distribution of the underground medium, which can effectively solve the problem of a large number of inversion parameter variables.

[0011] Optionally, the inversion target parameter is: resistivity or apparent resistivity or impedance or phase.

[0012] Optionally, a Fourier series expansion method is introduced in step S2 to express the resistivity of each underground medium layer in the underground layer model as a function of the depth of each layer, as follows:

[0013]

[0014] Where, FS x (L i ) represents the resistivity of the underground medium in the x direction The Fourier series function of the relationship with depth, that is, the Fourier series about Approximation of FS y (L i ) represents the resistivity of the underground medium in the y direction The Fourier series function of the relationship with depth, that is, the Fourier series about Approximation of; L i is the depth function of the ith underground medium;

[0015] and is the intermediate variable, and and resistivity and The relationship is:

[0016]

[0017]

[0018] In the formula, ρ l and ρ u The resistivity lower and upper limits of the inversion resistivity model preset in the inversion process are generally empirical values. Optionally, the depth function L i is the logarithmic function of depth, expressed as:

[0019] L i =log 10 (Z i ),i=1,…,nl,

[0020] In the formula, nl represents the total number of underground medium layers, Z i Indicates the depth of the center of the ith underground medium, L i Represents the logarithm of the depth of the center position of the i-th layer of underground medium.

[0021] Optionally, the minimum value of the inversion constraint function is the objective function, and the inversion constraint function is expressed as:

[0022] ψ(X)=ψ d (X)+λψ m (X)

[0023] Where X is the inversion parameter, which is composed of Fourier coefficients; λ represents the regularization factor, ψ d (X) represents the data fitting term, ψ m (X) represents the model smoothing term and satisfies:

[0024]

[0025]

[0026] Where, d f (X) is the forward response of the inversion area obtained by forward modeling based on the resistivity of each layer of underground medium; d obs is the forward response observation value observed by the equipment in the inversion area, that is, the observation data, ρ represents the resistivity, W m is the model roughness matrix, ‖‖1 and ||||2 are the Euclidean 1-norm and 2-norm respectively, W d N d ×N d The diagonal data weight matrix of dimension is represented by:

[0027]

[0028] Where ri is the noise percentage in the ith observation data point, η is the predefined minimum noise level set to avoid possible computational overflow problems, and N d is the number of observation data, is the observation data corresponding to the i-th observation data point.

[0029] As can be seen from the above formula, the technical solution of the present invention proposes a Tikhonov regularization objective function based on the L1 norm. Traditional inversion methods usually rely on the L2 norm to construct a regularized objective function. Although this method performs well when dealing with smooth resistivity distributions, its effect is often limited in underground structures with drastic changes or discontinuities. The present application introduces the L1 norm, which has a natural sparsity characteristic and is more suitable for capturing drastic changes in resistivity in underground media, such as faults, fissures, and highly anisotropic regions. This regularization method based on the L1 norm can more accurately describe the sudden change of resistivity in the objective function, significantly improve the resolution and accuracy of the inversion results, and especially when dealing with complex anisotropic media, can more effectively reveal the true distribution of underground resistivity.

[0030] Optionally, the ADMM algorithm is introduced to minimize the inversion constraint function containing the model smoothness term of L1 regularization, and then the objective function based on the inversion constraint function is expressed as:

[0031]

[0032] wxya m ρ=W,

[0033] In the formula, W is the intermediate variable, N fs The number of harmonics included is given by the objective function. Three variables are iteratively updated during the inversion operation: the inversion parameter X, the intermediate variable W, and the defined variable The update formula is as follows:

[0034]

[0035] U k+1 =U k +τ(W m ρ k+1 -W k+1 ),

[0036] Where, X k+1 , X k are the inversion parameters X, U in the k+1th and kth inversion iterations k+1 , U k is the variable U in the k+1th and kth inversion iterations, α k is the step size for updating the inversion parameter X, gk is the gradient of the inversion constraint function, H k is the Hessian matrix, W k+1 is the intermediate variable W of the k+1th inversion iteration operation, τ is the step coefficient, and its value is In the range, ρ k+1 is the resistivity in the k+1th inversion iteration; u is the Lagrange multiplier used to satisfy the equation constraints in the equation, and β is defined as the coefficient of the quadratic penalty term.

[0037] Among them, the ADMM algorithm is a distributed optimization method, which is particularly good at handling large-scale optimization problems with complex constraints. The ADMM introduced in the present invention decomposes the original objective function into more manageable sub-problems, and alternately updates the inversion variables and auxiliary variables in each iteration, thereby gradually approaching the optimal solution. This alternating update mechanism effectively reduces the solution space in the inversion process, making the solution process more stable and fast.

[0038] Optionally, invert the gradient of the constraint function g k , the Hessian matrix H k It is expressed as:

[0039]

[0040] In the formula, J k is the sensitivity matrix, T is the matrix transpose symbol; T k Represents a matrix at the kth iteration with a dimension of 2nl×N m , N m is the number of inversion parameters X, nl represents the total number of underground media layers, d f is the forward response, d obs is the ground observation value, i.e., the observation data;

[0041] Sensitivity matrix J k The elements of are calculated as:

[0042]

[0043] Where ΔX k (j) represents the change of the independent variable of the underground model, and are the forward responses calculated using the model variables before and after the j-th model variable disturbance, respectively, and the model variables refer to the inversion variables.

[0044] Optionally, the method further comprises: setting boundary constraints of the inversion variables, that is, setting upper and lower limits of the resistivity, thereby constraining the search range of the inversion variables in the inversion operation.

[0045] In the inversion process, the boundary constraints of the restored resistivity model are innovatively introduced to further optimize the search range of the solution space. This boundary constraint can not only greatly reduce the solution space and improve the efficiency of the inversion calculation, but also ensure the physical interpretability of the inversion results and avoid the non-physical deviation of the solution. Especially when facing problems with multiple solutions, this innovative measure can significantly improve the inversion results and make the restored underground resistivity distribution more in line with geological laws. Among them, traditional inversion methods often lead to too broad a solution space due to the lack of an effective constraint mechanism, and then unstable or unreasonable solutions appear, especially in the face of complex geological conditions or large data noise, which is more likely to cause overfitting problems. The technical solution of the present invention effectively limits the resistivity value to a physically acceptable range by imposing reasonable boundary constraints on the resistivity model during the inversion process, thereby avoiding the arbitrariness of the model parameters and the distortion of the inversion results.

[0046] In another aspect, the present invention provides an inversion system based on the above method, comprising:

[0047] A depth information acquisition module is used to acquire the depth information of the underground layer model in the inversion area;

[0048] A processing module is used to introduce a Fourier series expansion method to express the inversion target parameter of each layer of underground medium in the underground layer model as a function of the depth of each layer, that is, to approximate the inversion target parameter by Fourier series, and to convert the original inversion parameter composed of the inversion target parameter of each layer into the inversion parameter constructed by Fourier coefficients;

[0049] The inversion module is used to set the inversion constraint function to perform inversion iteration operations and obtain the resistivity of each layer of underground media in the inversion area;

[0050] Among them, each iterative operation determines the inversion target parameters of each layer of underground media based on the function constructed in step S2, the inversion parameters and the depth information; then, a forward operation is performed based on the inversion target parameters of each layer of underground media to obtain the forward response of the inversion area; and then, based on the inversion constraint function, the forward response and the forward response observation value obtained by the equipment observation in the inversion area, the operation is performed to update the inversion parameters.

[0051] In a third aspect, the present invention provides an electronic terminal, comprising one or more processors and a memory storing one or more computer programs;

[0052] Among them, the processing calls the computer program to implement: a one-dimensional magnetotelluric anisotropy inversion method based on Fourier series.

[0053] In a fourth aspect, the present invention provides a computer-readable storage medium storing a computer program, wherein the computer program is called by a processor to implement: steps of a one-dimensional magnetotelluric anisotropy inversion method based on Fourier series.

[0054] Beneficial Effects

[0055] Compared with the existing method, the advantages of the present invention are:

[0056] 1. The technical solution of the present invention innovatively proposes to use the Fourier series expansion method to approximate the resistivity distribution of the underground medium, that is, to approximate the resistivity through the Fourier series, and convert the original inversion parameters composed of the resistivity of each layer into inversion parameters constructed by Fourier coefficients, so as to achieve the use of a small number of Fourier coefficients to represent the complex resistivity distribution. This not only significantly reduces the number of inversion variables, but also can effectively compress the model parameter space, thereby reducing the computational burden of the inversion process. In addition, since the Fourier series has good approximation properties in mathematical theory, this method can also ensure a high-precision description of the underground resistivity distribution, making the inversion results more stable and reasonable. It should be understood that in addition to resistivity, it is also applicable to other categories of inversion target parameters, and the technical ideas implemented are consistent.

[0057] 2. In a further preferred embodiment of the present invention, a Tikhonov regularization objective function based on the L1 norm is used: that is, a Tikhonov regularization objective function based on the L1 norm is proposed. The introduced L1 norm has a natural sparsity characteristic and is more suitable for capturing the drastic changes in resistivity in underground media, such as faults, cracks and highly anisotropic areas. The constructed objective function can more accurately describe the sudden change of resistivity, significantly improve the resolution and accuracy of the inversion results, and especially when dealing with complex anisotropic media, can more effectively reveal the true distribution of underground resistivity.

[0058] 3. In a further preferred embodiment of the present invention, an alternating direction multiplier method (ADMM) algorithm is innovatively introduced to efficiently optimize the constructed objective function. In the present application, ADMM gradually approaches the optimal solution by decomposing the original objective function into more tractable sub-problems and alternately updating the inversion variables and auxiliary variables in each iteration. The solution space in the inversion process is effectively reduced by the alternating update mechanism, making the solution process more stable and fast. Especially compared to traditional optimization algorithms, ADMM updates the inversion variables and auxiliary variables alternately in sequence when dealing with nonlinear and non-convex problems. The ADMM algorithm can not only handle more complex objective function forms, but also avoid local optimal solutions, ensuring the global optimality of the inversion results. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1It is a schematic diagram of the data processing flow of the inversion method of the present invention;

[0060] Figure 2 It is a schematic diagram of the inversion results of the four models set in the example. DETAILED DESCRIPTION

[0061] The present invention provides a one-dimensional magnetotelluric anisotropy inversion method and system based on Fourier series, which is used to realize the inversion calculation of the target inversion parameters of each underground medium in the underground layer model, greatly reducing the number of inversion variables, and significantly improving the calculation efficiency and reliability. The core is to introduce the Fourier series expansion method, creatively approximate the distribution of the target inversion parameters of the underground medium by the Fourier series expansion method; in addition, the objective function is basically optimized in many aspects, such as: introducing the Tikhonov regularization objective function based on the L1 norm; and introducing the alternating direction multiplier method (ADMM) algorithm, which is used to efficiently optimize the constructed objective function.

[0062] Based on this, the basic technical idea of ​​a one-dimensional magnetotelluric anisotropy inversion method based on Fourier series provided by the present invention is as follows:

[0063] S1: Obtain the depth information of the underground layer model in the inversion area;

[0064] S2: A Fourier series expansion method is introduced to express the inversion target parameters of each layer of underground medium in the underground layer model as a function of the depth of each layer, that is, the inversion target parameters are approximated by Fourier series, and the original inversion parameters composed of the inversion target parameters of each layer are converted into inversion parameters constructed by Fourier coefficients;

[0065] S3: Set the inversion constraint function to perform inversion iteration operation to obtain the target inversion parameters of each layer of underground media in the inversion area.

[0066] The present invention will be further described below in conjunction with the best embodiment.

[0067] Embodiment 1:

[0068] This embodiment is described by taking resistivity as an example, and provides a one-dimensional magnetotelluric anisotropy inversion method based on Fourier series, comprising the following steps:

[0069] Step 1: Obtain the depth information of the underground layer model in the inversion area, that is, calculate the logarithm of the depth of the center position of each layer of underground medium:

[0070] L i =log 10 (Z i ),i=1,…,nl,

[0071] In the formula, nl represents the total number of underground medium layers, Z i Indicates the depth of the center of the ith underground medium, L i Represents the logarithm of the depth of the center position of the i-th layer of underground medium.

[0072] It should be understood that this embodiment preferably calculates the logarithm of the depth L i As a function variable. Among them, since the resistivity of the underground medium usually changes nonlinearly with increasing depth, if the depth is used directly as a function variable, layers with larger or smaller depths may cause excessive numerical differences, especially in deep underground layers, where the magnitude of the change may be much higher than that of the surface layer. By taking the logarithm of the depth, the depth difference will change relatively smoothly, which can avoid numerical instability or convergence difficulties caused by large depth differences during calculations. Selecting the logarithm of the depth as the depth function reduces the scale difference between different depth layers.

[0073] Step 2: The Fourier series expansion method is introduced to express the resistivity of each layer of underground medium in the underground layer model as a function of the depth of each layer. That is, the resistivity is approximated by the Fourier series, and the original inversion parameters composed of the resistivity of each layer are converted into inversion parameters constructed by Fourier coefficients.

[0074] In the one-dimensional underground model, this embodiment only considers the resistivity in the x direction of each layer of underground medium. And the resistivity in the y direction Using FS x (L i ) represents the resistivity of the underground medium in the x direction The Fourier series function of the relationship with depth, that is, the Fourier series about Approximation of FS y (L i ) represents the resistivity of the underground medium in the y direction The Fourier series function of the relationship with depth, that is, the Fourier series about Approximation of , there exists:

[0075]

[0076] In the formula, and is the intermediate variable, and and resistivity and The relationship is:

[0077]

[0078]

[0079] About FS x (L i ) and FS y (L i ) is to introduce the Fourier series expansion method, and express the resistivity of each underground medium in the underground layer model as a function of the depth of each layer. The general form of the Fourier series expansion exemplified in this embodiment is:

[0080]

[0081] Where y(x) is the Fourier series expansion, which represents resistivity in this application; the independent variable x represents the depth / depth function in this application; thus, when the depth is given, the resistivity corresponding to the depth can be obtained, p represents the period, and A n and B n Represent the coefficients of the nth cosine function and sine function, N fs Indicates the number of harmonics contained; A0 is a constant term in the Fourier series expansion coefficient, so N is set fs The value is inverted, then the vector updated in the inversion process is the Fourier coefficient The total composition is 2N fs +1 item.

[0082] It should be understood that in other embodiments, if the above-mentioned operations are performed on non-resistivity inversion target parameters, the target inversion parameters are represented by the Fourier series expansion of y(x), x is the depth / depth function, and the corresponding Fourier series expansion is constructed. This also achieves the conversion of the original inversion parameters composed of the inversion target parameters of each layer into inversion parameters constructed by Fourier coefficients.

[0083] When applied to this embodiment, the logarithm of the depth L i As a function variable, FS x (L i ) and FS y (L i ) is in the Fourier series expansion form:

[0084]

[0085]

[0086] Then these two sets of Fourier coefficients are represented by a vector:

[0087]

[0088] It should be understood that if the underground layer model is divided into 120 layers, the original inversion variable resistivity and 120 each; after the Fourier series expansion is introduced, the inversion variables are all 2Nfs +1 item, N in this example fs It is 12, which greatly reduces the inversion variables. It should also be understood that this embodiment selects the logarithm of the depth as the function variable as the best embodiment. In other feasible embodiments, selecting other depth functions also falls within the protection scope of the present invention and conforms to the technical ideas of the present invention. The overall effect has local differences, which does not affect the implementation of the technical solution of the present invention.

[0089] Step 3: Set the inversion constraint function to perform inversion iteration operation to obtain the resistivity of each layer of underground medium in the inversion area.

[0090] The inversion constraint function set in this embodiment is expressed as:

[0091] ψ(X)=ψ d (X)+λψ m (X)

[0092] Where X is the inversion parameter, which is composed of Fourier coefficients; λ represents the regularization factor, ψ d (X) represents the data fitting term, ψ m (X) represents the model smoothing term and satisfies:

[0093]

[0094]

[0095] Where, d f (X) is the forward response of the inversion area obtained by forward modeling based on the resistivity of each layer of underground medium; d obs is the forward response observation value observed by the equipment in the inversion area, ρ represents resistivity, W m is the model roughness matrix, W d N d ×N d The diagonal data weight matrix of dimension is represented by:

[0096]

[0097] In the formula, r i represents the noise percentage in the i-th observed data point, η is the predefined minimum noise level set, N d is the number of observation data.

[0098] In order to further optimize the inversion results, this embodiment introduces the ADMM algorithm to minimize the inversion constraint function containing the model smoothness term of L1 regularization, and then the objective function based on the inversion constraint function is expressed as:

[0099]

[0100] wxya m ρ=W,

[0101] The augmented Lagrangian function is:

[0102]

[0103] Defining variables

[0104] Based on the objective function, three variables are iteratively updated during the inversion operation: the inversion parameter X, the intermediate variable W, and the variable The update formula is as follows:

[0105]

[0106] U k+1 =U k +τ(W m ρ k+1 -W k+1 ),

[0107] Where, X k+1 , X k is the inversion parameter X in the k+1th and kth inversion iterations, α k is the step size for updating the inversion parameter X, g k is the gradient of the inversion constraint function, H k is the Hessian matrix, W k+1 is the intermediate variable W of the k+1th inversion iteration operation, τ is the step coefficient, and its value is In the range, ρ k+1 is the resistivity in the k+1th inversion iteration operation; u is the Lagrange multiplier used to satisfy the equation constraint in the equation, and β is the quadratic penalty term coefficient.

[0108] The gradient of the inversion constraint function g k , the Hessian matrix H k It is expressed as:

[0109]

[0110] In the formula, J k is the sensitivity matrix, T is the matrix transpose symbol; T k Represents a matrix at the kth iteration with a dimension of 2nl×N m , N m is the number of parameters of the inversion parameter X;

[0111] Since the forward modeling scale of the 1D model is relatively small, and the inversion variables brought by the Fourier series approximation of the underground model are relatively small, and the sensitivity matrix J kThe number of rows is small. Therefore, the perturbation technique can be used to calculate the Jacobian matrix, the sensitivity matrix J k The elements of are calculated as:

[0112]

[0113] Where ΔX k (j) represents the change of the independent variable of the underground model, and are the forward responses calculated using the model variables before and after the jth model variable is disturbed, respectively.

[0114] Based on the above mathematical model, each inversion iteration operation is based on the function constructed in step S2, the inversion parameters and the depth information to determine the resistivity of each layer of underground media; then the forward response of the inversion area is obtained based on the resistivity of each layer of underground media and the forward operation; then the inversion constraint function, the forward response and the forward response observation value observed by the equipment in the inversion area are calculated to update the inversion parameter X, the intermediate variable W and the defined

[0115] This embodiment uses the data fitting difference to determine whether convergence occurs. The data fitting difference is represented by the root mean square error of the data, as follows:

[0116]

[0117] Examples:

[0118] Four structural anisotropic models were established to generate the observed data d obs , these data include two impedance components Z at 18 frequencies ranging from 0.001 to 131 Hz xy and Z yx , where in the one-dimensional underground model, the electrical structure of the earth is assumed to vary only with depth, not with lateral variations. The lateral component of the electric field does not produce a corresponding lateral component of the magnetic field, and vice versa. Therefore, only two elements of the impedance tensor need to be considered, which are:

[0119]

[0120] In the formula, ω is the angular frequency, E represents the electric field intensity, and H represents the magnetic field. That is, after the electric field and magnetic field are obtained by forward modeling, the final data (forward response, i.e., forward response) is obtained by substituting them into the above impedance formula.

[0121] Each model has a total of 36 simulated complex data points, and 1% Gaussian noise is added to these synthetic data. The traditional Gauss-Newton algorithm based on the L2 norm and the Fourier series-based ADMM algorithm (using the L1 norm) are used to invert the four sets of synthetic data. The number of layers for all inversions is fixed at 100, and the thickness of each layer gradually increases from the first to the last layer. When performing inversion experiments using the Fourier series-based ADMM technique, the number of harmonics of the Fourier series expansion is set to 21. The initial model for all experiments is a half-space model with a resistivity of, and the upper and lower resistivity values ​​of these inversions are set to and, respectively. The inversion results of the four models are as follows Figure 2 exhibit.

[0122] In addition, a series of inversion experiments were conducted using synthetic and measured datasets to evaluate the performance of the proposed one-dimensional MT anisotropic inversion method based on the ADMM algorithm and Fourier series expansion. These experiments were performed on a laptop equipped with an 8-core Inteli7-11800H processor (main frequency 4.6GHz) and 32GB of memory, and Python programming was used. During the inversion process, we set three termination conditions:

[0123] 1) reaches the preset RMS value of 1.0,

[0124] 2) The RMS value does not decrease in three consecutive iterations,

[0125] 3) The preset maximum number of iterations is reached. Once any of these conditions is met, the inversion process is completed.

[0126] Through comprehensive testing of synthetic models, the ADMM inversion method based on Fourier series has advantages in robustness, reliability and efficiency in recovering underground resistivity distribution compared with the traditional GN algorithm based on L2 norm objective function. In addition, the influence of the number of harmonics in the Fourier series on the inversion results was studied, verifying the need to use a sufficient number of terms in the Fourier series to accurately interpret the observed data. Finally, the code was successfully applied to interpret actual AMT data from Qinghai Province. Compared with the traditional GN algorithm, the inversion results of the ADMM algorithm proposed in this invention are more consistent with the actual geological and drilling data, further verifying the practicality and reliability of the inversion algorithm proposed in this invention.

[0127] Embodiment 2:

[0128] This embodiment provides an inversion system based on the above inversion method, including: a depth information acquisition module, a processing module, and an inversion module connected in sequence.

[0129] Wherein, the depth information acquisition module is used to obtain the depth information of the underground layer model of the inversion area;

[0130] The processing module is used to introduce the Fourier series expansion method to express the inversion target parameter of each layer of underground medium in the underground layer model as a function of the depth of each layer, that is, to approximate the inversion target parameter by Fourier series, and to convert the original inversion parameter composed of the inversion target parameter of each layer into the inversion parameter constructed by Fourier coefficients;

[0131] The inversion module is used to set the inversion constraint function to perform inversion iteration operations and obtain the resistivity of each layer of underground media in the inversion area;

[0132] Among them, each iterative operation determines the inversion target parameters of each layer of underground media based on the function constructed in step S2, the inversion parameters and the depth information; then, a forward operation is performed based on the inversion target parameters of each layer of underground media to obtain the forward response of the inversion area; and then, based on the inversion constraint function, the forward response and the forward response observation value obtained by the equipment observation in the inversion area, the operation is performed to update the inversion parameters.

[0133] It should be understood that the specific implementation process of each module please refer to the above method content, the present invention will not be repeated here, and the division of the above functional modules is only for example illustration. In some embodiments, some functional modules can be merged, and some functional modules can be split. Each functional module can be implemented in software or hardware or a combination of software and hardware. Among them, the hardware and software equipment includes but is not limited to general-purpose computer equipment, programmable gate arrays, digital signal processors, microprocessors and their corresponding programming or burning software.

[0134] Embodiment 3:

[0135] This embodiment provides an electronic terminal, including one or more processors and a memory storing one or more computer programs;

[0136] Among them, the processing calls the computer program to implement: a one-dimensional magnetotelluric anisotropy inversion method based on Fourier series.

[0137] For the specific implementation process of each step, please refer to the description of the above method.

[0138] It should be understood that in the embodiments of the present invention, the processor may be a central processing unit (CPU), and the processor may also be other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or the processor may also be any conventional processor, etc. The memory may include a read-only memory and a random access memory, and provide instructions and data to the processor. A portion of the memory may also include a non-volatile random access memory. For example, the memory may also store information about the type of device.

[0139] Embodiment 4:

[0140] The present invention provides a computer-readable storage medium storing a computer program, wherein the computer program is called by a processor to implement: steps of a one-dimensional magnetotelluric anisotropy inversion method based on Fourier series.

[0141] For the specific implementation process of each step, please refer to the description of the above method.

[0142] The readable storage medium is a computer-readable storage medium, which may be an internal storage unit of the software and hardware device described in any of the aforementioned embodiments, such as a hard disk or memory of a controller. The readable storage medium may also be an external storage device of the controller, such as a plug-in hard disk, a smart memory card (Smart Media Card, SMC), a secure digital (Secure Digital, SD) card, a flash card (Flash Card), etc. equipped on the controller. Furthermore, the readable storage medium may also include both an internal storage unit of the controller and an external storage device. The readable storage medium is used to store the computer program and other programs and data required by the controller. The readable storage medium may also be used to temporarily store data that has been output or is to be output.

[0143] Based on this understanding, the technical solution of the present invention, in essence or in other words, the part that contributes to the prior art, or all or part of the technical solution can be embodied in the form of a software product, which is stored in a storage medium and includes several instructions for enabling a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention. The aforementioned readable storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk, etc., various media that can store program codes.

[0144] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-readable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program codes. The present application is a device for realizing the function specified in one flow or multiple flows of the flow chart and / or one or multiple blocks of the block diagram with reference to the instructions executed by the processor according to the method, device (system) and computer program product of the embodiment of the present application. These computer program instructions can also be stored in a computer-readable memory that can guide a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer-readable memory produce a product including an instruction device, which realizes the function specified in one flow or multiple flows of the flow chart and / or one or multiple blocks of the block diagram. These computer program instructions may also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, whereby the instructions executed on the computer or other programmable device provide steps for implementing the functions specified in one or more processes in the flowchart and / or one or more boxes in the block diagram.

[0145] It should be emphasized that the examples described in the present invention are illustrative rather than restrictive, and therefore the present invention is not limited to the examples described in the specific embodiments. Any other embodiments derived by those skilled in the art based on the technical solution of the present invention that do not depart from the purpose and scope of the present invention, whether modified or replaced, also fall within the scope of protection of the present invention.

Claims

1. A one-dimensional magnetotelluric anisotropy inversion method based on Fourier series, characterized by: The following steps are involved: S1: Obtain the depth information of the underground layer model in the inversion area; S2: A Fourier series expansion method is introduced to express the inversion target parameters of each layer of underground medium in the underground layer model as a function of the depth of each layer, that is, the inversion target parameters are approximated by Fourier series, and the original inversion parameters composed of the inversion target parameters of each layer are converted into inversion parameters constructed by Fourier coefficients; S3: Set the inversion constraint function to perform inversion iteration operation to obtain the resistivity of each layer of underground medium in the inversion area; Among them, each iterative operation determines the inversion target parameters of each layer of underground media based on the function constructed in step S2, the inversion parameters and the depth information; then, a forward operation is performed based on the inversion target parameters of each layer of underground media to obtain the forward response of the inversion area; and then, based on the inversion constraint function, the forward response and the forward response observation value obtained by the equipment observation in the inversion area, the operation is performed to update the inversion parameters.

2. The method according to claim 1, characterized in that: The inversion target parameter is: resistivity or apparent resistivity or impedance or phase.

3. The method according to claim 2, characterized in that: In step S2, the Fourier series expansion method is introduced to express the resistivity of each underground medium in the underground layer model as a function of the depth of each layer, as follows: Where, FS x (L i ) represents the resistivity of the underground medium in the x direction The Fourier series function of the relationship with depth, that is, the Fourier series about Approximation of FS y (L i ) represents the resistivity of the underground medium in the y direction The Fourier series function of the relationship with depth, that is, the Fourier series about Approximation of; L i is the depth function of the ith underground medium; and is the intermediate variable, and and resistivity and The relationship is: In the formula, ρ l and ρ u The resistivity lower and upper limits of the inversion resistivity model preset during the inversion process.

4. The method according to claim 3, characterized in that: Depth function L i is the logarithmic function of depth, expressed as: L i =log 10 (Z i ),i=1,…,nl, In the formula, nl represents the total number of underground medium layers, Z i Indicates the depth of the center of the ith underground medium, L i Represents the logarithm of the depth of the center position of the i-th layer of underground medium.

5. The method according to claim 3, characterized in that: The inversion constraint function is expressed as: ψ(X)=ψ d (X)+lψ m (X) Where X is the inversion parameter, which is composed of Fourier coefficients; λ represents the regularization factor, ψ d (X) represents the data fitting term, ψ m (X) represents the model smoothing term and satisfies: Where, d f (X) is the forward response of the inversion area obtained by forward modeling based on the resistivity of each layer of underground medium; d obs is the forward response observation value observed by the equipment in the inversion area, that is, the observation data, ρ represents the resistivity, W m is the model roughness matrix, || ||1 and || ||2 are the Euclidean 1-norm and 2-norm respectively, W d N d ×N d The diagonal data weight matrix of dimension is represented by: In the formula, r i represents the noise percentage in the i-th observed data point, η is the predefined minimum noise level set, N d is the number of observation data, is the observation data corresponding to the i-th observation data point.

6. The method according to claim 5, characterized in that The ADMM algorithm is introduced to minimize the inversion constraint function of the model smoothness term containing L1 regularization, and then the objective function based on the inversion constraint function is expressed as: In the formula, W is the defined intermediate variable, N fs The number of harmonics included is given by the objective function. Three variables are iteratively updated during the inversion operation: the inversion parameter X, the intermediate variable W, and the defined variable The corresponding update formula is as follows: U k+1 =U k +τ(W m r k+1 -W k+1 ), Where, X k+1 , X k are the inversion parameters X, U in the k+1th and kth inversion iterations k+1 , U k is the variable U in the k+1th and kth inversion iterations, α k is the step size for updating the inversion parameter X, g k is the gradient of the inversion constraint function, H k is the Hessian matrix, W k+1 is the intermediate variable W of the k+1th inversion iteration operation, τ is the step coefficient, and its value is In the range, ρ k+1 is the resistivity in the k+1th inversion iteration; u is the Lagrange multiplier used to satisfy the equation constraint in the equation, and β is the coefficient of the quadratic penalty term.

7. The method according to claim 6, characterized in that: The gradient of the inversion constraint function g k , the Hessian matrix H k The calculation is expressed as: In the formula, J k is the sensitivity matrix, T is the matrix transpose symbol; T k Represents a matrix at the kth iteration with a dimension of 2nl×N m , N m is the number of parameters of the inversion parameter X, nl represents the total number of underground media layers, d f is the forward response, d obs Observations at the surface, i.e., observation data; Sensitivity matrix J k The elements of are calculated as: Where ΔX k (j) represents the change in the model variables of the underground model, and are the forward responses calculated using the model variables before and after the j-th model variable disturbance, respectively, and the model variables refer to the inversion variables.

8. An inversion system based on the method according to any one of claims 1 to 7, characterized in that: include: A depth information acquisition module is used to acquire the depth information of the underground layer model in the inversion area; A processing module is used to introduce a Fourier series expansion method to express the inversion target parameter of each layer of underground medium in the underground layer model as a function of the depth of each layer, that is, to approximate the inversion target parameter by Fourier series, and to convert the original inversion parameter composed of the inversion target parameter of each layer into the inversion parameter constructed by Fourier coefficients; The inversion module is used to set the inversion constraint function to perform inversion iteration operations and obtain the resistivity of each layer of underground media in the inversion area; Among them, each iterative operation determines the inversion target parameters of each layer of underground media based on the function constructed in step S2, the inversion parameters and the depth information; then, a forward operation is performed based on the inversion target parameters of each layer of underground media to obtain the forward response of the inversion area; and then, based on the inversion constraint function, the forward response and the forward response observation value obtained by the equipment observation in the inversion area, the operation is performed to update the inversion parameters.

9. An electronic terminal, characterized in that: include: one or more processors; a memory storing one or more computer programs; The processor calls the computer program to implement: The steps of the method according to any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that: A computer program is stored, which is called by a processor to implement: The steps of the method according to any one of claims 1 to 7.

Images