Gravity-magnetoelectricity three-dimensional joint inversion method based on reweighting

Abstract

The invention discloses a gravity-magnetoelectricity three-dimensional joint inversion method based on reweighting, and belongs to the technical field of geophysics. According to the method, the resolution and physical property difference characteristics of data of geophysical gravity, a magnetic method and an electromagnetic method are fully utilized, and a more reliable underground structure model is obtained in a joint inversion mode; according to the method, gravity, magnetic and electric multi-field data are processed at the same time, and the method is suitable for mineral resource exploration, oil and gas reservoir exploration, deep structure research and other scenes under complex geological conditions, and especially has remarkable advantages in hidden ore body boundary recognition, lithologic contact zone positioning and the like; by integrating reweighting, cross gradient and efficient iteration technologies, the manual intervention and parameter debugging cost is greatly reduced, the automatic inversion process from data input to model output is realized, and the engineering application efficiency of geophysical exploration is improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of geophysics, and in particular to a method for gravity-magnetic-electric three-dimensional joint inversion based on reweighting. BACKGROUND

[0002] In geophysical exploration, joint interpretation of different physical methods has been an important research direction to improve the accuracy of underground structure detection. Due to the inherent resolution limitations of traditional single geophysical methods, it is often difficult to accurately depict the complex characteristics of underground media. Gravity, magnetic and magnetotelluric methods each have their own characteristics, for example, gravity and magnetic data provide high lateral resolution and are mainly used for detecting large-scale geological structures; while magnetotelluric data have significant advantages in vertical resolution and can effectively image deep conductive structures. However, due to the essential differences in the response mechanisms of various physical methods to underground structures, when different geophysical methods are independently inverted and interpreted, significant differences in the underground conceptual model and ambiguity and conflicts in the identification of target areas are caused. In order to overcome this limitation, academia and engineering circles have increasingly paid attention to the comprehensive interpretation technology of multiple geophysical data sets, especially joint inversion methods.

[0003] In the prior art, joint inversion methods can be roughly divided into two main strategies: the first strategy enhances the structural similarity between different physical models through structural coupling constraints, such as cross-gradient constraints, cosine constraints, etc.; the second strategy attempts to establish a correlation between different physical parameters, mainly through empirical rock physics relationships or clustering methods. However, both of these two strategies have drawbacks. The coupling effect of the structural coupling strategy is relatively poor, and the direct relationship method may introduce over-coupling. Traditional joint inversion attempts to integrate multi-field data, but lacks effective cross-physical field constraint means, and cannot fully utilize the correlation of different physical parameters in geological structures (such as the consistency of the density, magnetic and electrical boundaries of the same geological body), resulting in low structural matching degree of the inversion model and low computational efficiency.

[0004] Based on the above background, a gravity-magnetic-electric three-dimensional joint inversion method based on reweighting is proposed. SUMMARY

[0005] The purpose of the present application is to provide a gravity-magnetic-electric three-dimensional joint inversion method based on reweighting to solve the problems in the background art.

[0006] To achieve the above purpose, the present application provides a gravity-magnetic-electric three-dimensional joint inversion method based on reweighting, comprising the following steps:

[0007] S1, according to the accuracy requirements and data characteristics of the gravity-magnetic-electric three-dimensional joint inversion model, read the input data file, model file and inversion parameter file to obtain model parameters;

[0008] S2, obtaining the gravity and magnetic inversion model under the weighted domain by using the reweighting function, determining the initial model and the initial model parameter, the gravity and magnetic inversion model including the magnetotelluric resistivity model, the density model and the magnetic susceptibility model;

[0009] S3, updating the magnetotelluric resistivity model by the nonlinear conjugate gradient iterative search, updating the density model and the magnetic susceptibility model by the conjugate gradient iterative search, and then carrying out the joint constraint inversion by the cross gradient coupling;

[0010] S4, sequentially outputting the updated model of S3, carrying out the forward calculation on the output updated model, determining the new model parameter, and evaluating the model updating effect;

[0011] S5, carrying out the cyclic iteration on the updated model, and terminating the iteration when the preset maximum iteration number is reached or the fitting error is less than the specified range.

[0012] Preferably, in S1, the data file and the model file include the data and the initial model of the single physical property inversion; the inversion parameter file includes the maximum iteration number of inversion, the target fitting error threshold, the gravity and magnetic depth weighting coefficient, the initial model coefficient, the cross gradient coefficient, the Lp norm selection parameter and the reweighting weight factor.

[0013] Preferably, in S2, the gravity and magnetic inversion model under the weighted domain is obtained by using the reweighting function, the initial gradient of the model, the search direction and the search step are determined; the reweighting function is represented as:

[0014]

[0015] Wherein, W m is the gravity and magnetic inversion parameter model, is the L p norm weighting function, W mz is the depth weighting function.

[0016] The Lp norm weighting is represented as:

[0017]

[0018] The depth weighting function is represented as:

[0019]

[0020] Wherein, I is the unit matrix, A is the forward operator matrix, diag is the diagonal function, ζ is the normalized weight factor, m is the model matrix, e is a very small constant to prevent matrix singularity, H, z and z0 are the total vertical depth of the underground inversion space, the depth of model discretization and the depth of model top surface respectively.

[0021] Preferably, the initial model is represented as:

[0022] m w0 = W m m0;

[0023] The initial model parameters are:

[0024]

[0025] p0 = g0;

[0026]

[0027] wherein m0 is an initial model, g0 is an initial gradient, p0 is a search direction, and a0 is a search step; A is a forward operator matrix, d is observation data, and λ is a Lagrange multiplier.

[0028] Preferably, the specific steps of S3 are:

[0029] S31, conjugate gradient iteration is used for the density model and the magnetic susceptibility model, which is expressed as:

[0030]

[0031] m wk+1 = m wk + a k+1 p k ;

[0032]

[0033] p k+1 = g k+1 + β k p k ;

[0034] wherein A is a forward operator matrix, k is the number of iterations, and m and m w are model parameters in the actual and weighted domains, respectively, and a and β are step coefficients;

[0035] S32, nonlinear conjugate gradient iteration is used for the magnetotelluric resistivity model, which is expressed as:

[0036] Φ(m k + a k p k ) = min a Φ(m k + a k p k );

[0037] m k+1 = m k + a k p k ;

[0038] pk =-C k g k +β k p k-1 ;

[0039] Where C is the preconditioning factor, and the optimal step sizes α and β are obtained by quadratic and cubic interpolation approximations;

[0040] S33. Joint inversion of different physical property parameters is performed through cross-gradient coupling, and the objective function is expressed as:

[0041]

[0042] in, is the discrete component in the x direction, is the discrete component in the y direction, is the discrete component in the z direction, λ is the Lagrange multiplier, η cross gradient factor, I is the identity matrix;

[0043] At the beginning of the inversion, the initial Lagrange multiplier λ0 and the cross gradient coefficient η0 are given. During the inversion iteration process, λ and η are updated respectively by the ratio of the data fitting term to the model fitting term and the cross gradient term and the number of iterations k; λ k ,η k In the calculation formula, a k and b k It is the attenuation factor of each iteration (its value is less than 1), which is used to realize the attenuation of the proportion of model terms and cross gradient terms during the inversion iteration process, so as to achieve smooth convergence.

[0044] Preferably, the specific steps of S33 are:

[0045] 1) The cross-gradient representation of the physical property model matrix of different geophysical methods is:

[0046]

[0047] Among them, t is the three-dimensional cross gradient vector, m1 and m2 are different physical property model matrices, is the gradient operator;

[0048] For a three-dimensional model, the three components x, y, and z in the three-dimensional cross gradient vector are expressed as:

[0049]

[0050] 2) Based on the differential form, the three components of the three-dimensional cross gradient vector are discretized and expressed as:

[0051]

[0052] where D x , D y , D z are the first-order central difference operators in x, y, z directions respectively; the conventional treatment of the boundary derivatives is to replace them with forward and backward difference, here the derivatives at the boundary are not calculated, i.e. the cross-gradient does not impose structural constraints on the boundary.

[0053] 3) For joint inversion of different physical properties, the cross-gradient constraint is added to the regularization inversion to obtain a new regularization inversion formula based on the cross-gradient constraint, expressed as:

[0054] Φ = Φ d + λΦ m + ηΦ cg ;

[0055] where Φ d is the data fitting term, Φ m is the model fitting term, and Φ cg is the cross-gradient term.

[0056] 4) The three-component discrete in step 2) is substituted into the regularization inversion formula in step 3) to obtain the objective function, expressed as:

[0057]

[0058] 5) The objective function is differentiated with respect to the model parameters, expressed as:

[0059]

[0060] According to the above formula, the gradient can be obtained, and the gradient is updated by iterative solution to approximate the solution of the model m.

[0061] Preferably, in S4, forward calculation is performed on the updated model, and the specific steps are as follows:

[0062] 1) The gravity potential V and the gravity anomaly Δg of the updated density model are calculated, expressed as:

[0063]

[0064] where G is the gravitational constant, and r is the distance between the volume element Q(x, y, z) and the observation point P(X, Y, Z);

[0065] After discretization, we get:

[0066]

[0067] where Δx i = X-x i ; Δyi = Y - y i ; Δz i = Z - z i ;

[0068]

[0069] 2) The magnetic potential and magnetic anomaly calculation is performed on the updated susceptibility model, and when the magnetization is M, the Poisson formula between the gravity potential and the magnetic potential is satisfied and is expressed as:

[0070]

[0071]

[0072] where M x , M y , and M z are the components of the magnetization in the x, y, and z directions respectively, i is the magnetic inclination, and δ is the magnetic declination.

[0073] In addition, the magnetization and the magnetic induction have a relationship as follows:

[0074] B = μ0H;

[0075] where B is the magnetic induction, and H is the magnetization.

[0076] The relationship between the magnetic induction and the magnetization can be derived as follows:

[0077]

[0078] The magnetic anomaly B x , B y , and B z in the direction of the geomagnetic field are expressed as:

[0079] ΔT = B x cos I' cos A' + B y cos I' sin A' + B z sin I';

[0080] where B x is the magnetic induction in the x direction, B y is the magnetic induction in the y direction, B z is the magnetic induction in the z direction, I' is the geomagnetic inclination, and A' is the angle between the magnetic north and the x axis.

[0081] 3) The forward calculation is performed on the updated magnetotelluric resistivity model, and specifically:

[0082] In the magnetotelluric forward process, only the quasi-stable electromagnetic field is studied, and based on this premise, the Maxwell equations can be simplified as:

[0083]

[0084] Where i is an imaginary unit, ω is the angular frequency of the electromagnetic wave, which satisfies the relationship ω=2πf with the frequency f of the electromagnetic wave, σ is the conductivity, and μ is the magnetic permeability of the medium;

[0085] Generally speaking, the relative dielectric constant of underground rocks is usually between 1 and 50, and its conductivity is between 10 -4 To 10s / m, the frequency range used in magnetotelluric method is usually between 10 -4 to 10 3 Hz; therefore, the ratio of displacement current to conduction current is extremely small and can be ignored. Therefore, in the process of magnetotelluric forward modeling, only the quasi-steady electromagnetic field problem is studied. Based on this premise, and assuming that the underground is a uniform half space (the electromagnetic constant μ r , ε r Under the condition of 1), starting from the electric field, taking the curl operation on both sides of the first sub-equation of the equation system, we get:

[0086]

[0087] By expanding the electric field component E and the magnetic field component H along the x, y, and z directions and solving the partial differential equations, the distribution of each component of the electromagnetic field can be obtained;

[0088] Then, by introducing the wave impedance The apparent resistivity and impedance phase are further calculated as follows:

[0089]

[0090] Among them, ρ ij is the apparent resistivity, φ ij is the impedance phase, ω is the angular frequency of the electromagnetic wave, Z ij is the wave impedance.

[0091] Preferably, in S5, during the iterative process of the three-dimensional joint inversion, when the model is updated, the model fitting state is evaluated by calculating the root mean square error based on the pre-input iteration constraint conditions. When the model fitting error enters the preset ideal error range or reaches the maximum allowable number of iterations, the iterative process is stopped; otherwise, the iterative search is continued until the predetermined error range requirements are met.

[0092] Therefore, the present invention provides a method for three-dimensional joint inversion of gravity, magneto-electricity based on reweighting, which has the following beneficial effects:

[0093] (1) The combination of Lp norm and depth weighting function is used in the reweighting strategy, the Lp norm can realize sparse constraint, and the depth weighting can compensate the signal attenuation in the deep part, so that the inversion has advantages in both high resolution in the shallow part and weak signal recognition in the deep part.

[0094] (2) The structural consistency constraint of different physical property models (density, magnetic susceptibility, and resistivity) is realized by cross-gradient coupling, the spatial gradient correlation of multi-physical property parameters is introduced into the joint inversion objective function, the gradient boundary of different physical property parameters is forced to coincide, the inversion result is more consistent with the real distribution characteristics of the geological body, and the structural contradiction problem caused by independent inversion of physical property parameters in the traditional joint inversion is effectively solved; the cross-gradient factor is introduced, the constraint strength is gradually weakened in the iteration process through the adaptive attenuation of the ratio of the data fitting term to the model term, the model meets the structural consistency and avoids excessive smoothing, and the identification ability of the complex geological body is improved.

[0095] (3) The conjugate gradient iteration is used for the density / magnetic susceptibility model, and the nonlinear conjugate gradient iteration is used for the magnetotelluric resistivity model, the pre-condition factor C and the step search of interpolation approximation are combined, the convergence stability of the linear algorithm is utilized, the adaptability to high nonlinear problems (such as resistivity inversion) is enhanced through nonlinear optimization, and the number of iterations is reduced; and the model coefficient and the cross-gradient factor are updated in the iteration process, the weights of the data fitting term and the model constraint term are balanced through the attenuation factor, the iteration divergence caused by unreasonable initial parameter setting is avoided, and the robustness of the algorithm to noise data is improved.

[0096] The technical solutions of the present application will be further described in detail below with the help of the drawings and examples. DESCRIPTION OF DRAWINGS

[0097] Figure 1 The flowchart of the embodiment of the present application is shown in the figure.

[0098] Figure 2 The comparison chart of different depth weighting in the reweighting of the embodiment of the present application is shown in the figure, wherein (a) is L2 regularization without depth weighting, (b) is L2 regularization plus forward operator depth weighting, (c) is L2 regularization plus Dep1 depth weighting function, and (d) is L2 regularization plus Dep2 depth weighting.

[0099] Figure 3 The Lp norm weighting comparison chart in the reweighting of the embodiment of the present application is shown in the figure, wherein (a) is L0 regularization and Dep2 depth weighting, (b) is L1 regularization and Dep2 depth weighting, (c) is L2 regularization and Dep2 depth weighting, and (d) is L02 regularization and Dep2 depth weighting.

[0100] Figure 4Density model obtained by separate inversion and density model obtained by joint inversion of cross gradient, (c) is susceptibility model obtained by separate inversion and susceptibility model obtained by joint inversion of cross gradient, (d) is resistivity model obtained by separate inversion and resistivity model obtained by joint inversion of cross gradient;

[0101] Figure 5 Physical property cross plot of separate inversion and joint inversion of cross gradient of the embodiment of the present application, wherein (a) is the correlation of susceptibility and residual density obtained by separate inversion, (b) is the correlation of apparent resistivity and susceptibility obtained by separate inversion, (c) is the correlation of residual density and apparent resistivity obtained by separate inversion, (d) is the correlation of susceptibility and residual density obtained by joint inversion, (e) is the correlation of apparent resistivity and susceptibility obtained by joint inversion, and (f) is the correlation of residual density and apparent resistivity obtained by joint inversion;

[0102] Figure 6 Cross gradient value plot of separate inversion and joint inversion of cross gradient of the embodiment of the present application, wherein (a) is the cross gradient value of residual density and susceptibility obtained by separate inversion, (b) is the cross gradient value of apparent resistivity and susceptibility obtained by separate inversion, (c) is the cross gradient value of residual density and apparent resistivity obtained by separate inversion, (d) is the cross gradient value of residual density and susceptibility obtained by joint inversion, (e) is the cross gradient value of apparent resistivity and susceptibility obtained by joint inversion, and (f) is the cross gradient value of residual density and apparent resistivity obtained by joint inversion. DETAILED DESCRIPTION

[0103] The technical solutions of the present application will be further described below by means of the accompanying drawings and embodiments.

[0104] In order to make the objects, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the accompanying drawings of the embodiments of the present application. Obviously, the described embodiments are only some of the embodiments of the present application, not all the embodiments of the present application.

[0105] EMBODIMENT

[0106] As shown in the accompanying drawings and embodiments, Figure 1 the present application provides a method for three-dimensional joint inversion of gravity, magnetic and electric based on reweighting, which comprises the following steps:

[0107] S1, input parameters and data preparation: according to the accuracy requirement and data characteristics of the three-dimensional joint inversion model of gravity, magnetic and electric, read the data file, model file and inversion parameter file, obtain the separate physical property inversion response file d MT , Grav , Mag and the model file At the same time, key inversion parameters are extracted, including the maximum number of iterations iter max , target fitting error threshold rms target , gravity and magnetic depth weighting coefficient β, Lagrange multiplier λ, joint inversion weight factor η, Lp norm selection parameter ζ and reweighting weight factor w Lp .

[0108] S2. Construction of reweighted model and calculation of initial variables:

[0109] The reweighted operator is constructed using Lp norm weighting and depth weighting function to obtain the gravity and magnetic inversion model in the weighted domain. Calculate the initial gradient g0, search direction p0, search step α0 and cross gradient value Φ cross .

[0110] S3. Joint model iterative update:

[0111] For the gravity and magnetic part, the conjugate gradient method is used to update the density and magnetic susceptibility model m den ,m mag For the nonlinear part of magnetotelluric, the nonlinear conjugate gradient method is used to iteratively update the resistivity model m res ; Through cross-gradient coupling, joint constraints between multiple physical property parameters are achieved to improve model structure consistency and inversion accuracy.

[0112] S4. Model forward modeling and parameter updating:

[0113] Output updated model file The iteratively updated model is forward calculated (gravity anomaly, magnetic anomaly and magnetotelluric response), and the gradient g is recalculated based on the forward results. k , search step length p k and direction α k , used for the next iteration.

[0114] S5. Iteration termination condition judgment:

[0115] The root mean square error is used to evaluate the goodness of fit of the model; when the fitting error reaches the preset ideal range (rms <rms target ) or the number of iterations reaches the maximum limit (iter>iter max Stop the iteration; otherwise, continue the iterative search until the condition is met.

[0116] like Figure 2 As shown in the figure, a synthetic residual density model consisting of rectangular prisms and stepped anomalies is constructed by evaluating depth weighting. The model parameters are set as follows: the residual density of the stepped anomaly is 1.5 g / cm 3, central prism residual density 3 g / cm 3 , background density 2 g / cm 3 . The model is discretized into 30x20x30 cubic cells with a cell size of 500m x 500m x 500m. The survey grid spacing is 200m with a total of 3721 points. The data is contaminated with 5% noise. The reweighting scheme is constructed using different depth weighting matrices and L2 regularization weighting. The depth weighting functions are selected as the identity matrix, the forward operator matrix, the conventional depth weighting matrix and the improved depth weighting matrix. The results show that for the selection of depth weighting, when it is the identity matrix, the "skin effect" is obvious. The forward operator and the conventional depth weighting can improve the depth positioning but have "tail", while the improved depth weighting method can effectively suppress the "skin effect" and "tail", and the inversion result is more accurate and converges faster.

[0117] As shown in Figure 3 , a susceptibility model is constructed to evaluate different regularization methods. The model contains two step anomalies with susceptibility of 0.5 and 1.0, respectively. The magnetic parameters are set as: field strength 5x10 4 nT, magnetic inclination 45°, magnetic declination 0°. The model is discretized into 30x20x12 cubic cells with a cell size of 500m x 500m x 500m. The survey grid spacing is 200m with a total of 3876 points. The data is contaminated with 5% noise. The reweighting scheme is constructed using the improved depth weighting and different regularization methods. The conjugate gradient algorithm is used for iteration for 10 times. The first 3 iterations are inverted using only depth weighting, and the last 7 iterations use the reweighting scheme. The results show that for the selection of Lp norm weighting, L0 norm weighting inversion abnormal body focusing distortion, unable to accurately fit the target; L1 norm weighting inversion is slightly improved in target fitting, but still not accurate; L2 norm weighting inversion fitting result is more smooth, but there is obvious tailing and diffusion; elastic network regularization combines the advantages of L1 and L2 norms, which can more accurately fit the boundary, reduce divergence, and thus achieve more accurate fitting.

[0118] As shown in Figure 4 , a three-dimensional synthetic model containing five blocks is constructed to evaluate the effectiveness of the cross-gradient joint inversion algorithm. The model features at least one of the adjacent blocks sharing the same physical property (density, susceptibility or resistivity). The calculation domain is discretized into 40x40x12 grids with a horizontal cell size of 200m x 300m. The MT survey grid spacing is 300m with a total of 289 points. The frequency range is selected to be 10 -3 -10 4 Hz. The gravity and magnetic survey point spacing is 300m, with a total of 2500 points. The inversion uses impedance tensor components (Z xx , Z yy , Z xy , Z yx) As the magnetotelluric data, the gravity anomaly and the magnetic anomaly value are taken as the gravity and magnetic data respectively, the maximum iteration number is set to 50, the magnetotelluric inversion target RMS is 0.5, the initial regularization parameter is 1.0, and the gravity-magnetic inversion target is 0.05. The joint inversion strategy based on the cross-gradient constraint effectively overcomes the structure coupling problem of adjacent isotropic bodies in the single inversion, and more fits the real spatial-structure characteristics of the anomaly body.

[0119] As shown in Figure 5 、 Figure 6 , the physical property intersection and the cross-gradient value before and after the single inversion and the joint inversion are analyzed, which shows that the joint inversion significantly reduces the cross-gradient value, and highlights the linear correlation feature trend of the physical property distribution, indicating that the joint inversion enhances the structure correlation.

[0120] Therefore, the method of the three-dimensional joint inversion of gravity-magnetic-electromagnetic based on the reweighting, the innovative reweighting function realizes the inversion optimization of the gravity and magnetic data, the resolution and the physical property difference characteristics of the geophysical gravity, magnetic and electromagnetic data are utilized, the more reliable underground structure model is obtained through the joint inversion mode of optimizing and integrating the multi-source data, improving the model resolution and the reliability, the problems of the poor coupling effect, the weak model homogeneity and correlation and the low calculation efficiency in the traditional joint inversion are effectively solved, the inversion precision and the reliability of the three-dimensional geological model are significantly improved, and the method has important engineering application value.

[0121] Finally, it should be noted that: the above examples are only used to illustrate the technical solutions of the present application but not to limit it, although the present application has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that: it can still modify or equivalently replace the technical solutions of the present application, and these modifications or equivalent replacements also cannot make the modified technical solutions deviate from the spirit and scope of the technical solutions of the present application.

Claims

1. A method for three-dimensional joint inversion of gravity, magneto-electricity and gravity based on reweighting, characterized in that: The following steps are involved: S1. Read the input data file, model file and inversion parameter file according to the accuracy requirements and data characteristics of the gravity, magnetism and electricity three-dimensional joint inversion model to obtain model parameters; S2. Using a reweighting function to obtain a gravity-magnetic inversion model in a weighted domain, and determining an initial model and initial model parameters. The gravity-magnetic inversion model includes a magnetotelluric resistivity model, a density model, and a magnetic susceptibility model; S3, update the magnetotelluric resistivity model through nonlinear conjugate gradient iterative search, update the density model and magnetic susceptibility model through conjugate gradient iterative search, and then perform joint constrained inversion through cross-gradient coupling; S4, outputting the updated model of S3 one by one, performing forward calculation on the output updated model, determining new model parameters, and evaluating the model update effect; S5. Perform cyclic iterations on the updated model. When the preset maximum number of iterations is reached or the fitting error is less than a specified range, the iteration is terminated.

2. The method of three-dimensional joint inversion of gravity, magneto-electricity and gravity based on reweighting according to claim 1, characterized in that: In S1, the data file and model file include the data and initial model of the individual physical property inversion; the inversion parameter file includes the maximum number of inversion iterations, the target fitting error threshold, the gravity and magnetic depth weighting coefficient, the initial model coefficient, the cross gradient coefficient, the Lp norm selection parameter and the reweighting weight factor.

3. The method of three-dimensional joint inversion of gravity, magneto-electricity based on reweighting according to claim 1, characterized in that: In S2, the gravity and magnetic inversion model in the weighted domain is obtained by using the reweighting function to determine the initial gradient, search direction, and search step of the model; the reweighting function is expressed as: Among them, W m is the gravity and magnetic inversion parameter model, For L p Norm weighting function, W mz is the depth weighting function.

4. The method of three-dimensional joint inversion of gravity, magneto-electricity based on reweighting according to claim 3, characterized in that: In S2, the initial model is expressed as: m w0 =W m m0; The initial model parameters are: p0=g0; Among them, m0 is the initial model, g0 is the initial gradient, p0 is the search direction, α0 is the search step size; A is the forward operator matrix, d is the observation data, and λ is the Lagrange multiplier.

5. The method of three-dimensional joint inversion of gravity, magneto-electricity and gravity based on reweighting according to claim 1, characterized in that: The specific steps of S3 are: S31. Conjugate gradient iteration is used for density model and magnetic susceptibility model, which can be expressed as: m wk+1 =m wk +a k+1 p k ; p k+1 g k+1 +β k p k 100. Among them, k is the number of iterations, m and m w are the model parameters in the actual and weighted domains, respectively, α and β are the step size coefficients; S32. The nonlinear conjugate gradient iteration of the magnetotelluric resistivity model is expressed as: Φ(m k +a k p k )=minαΦ(m k +a k p k ); m k+1 =m k +α k p k ; p k D-C k g k +β k p k-1 100. Where C is the preconditioning factor; S33. Joint inversion of different physical property parameters is performed through cross-gradient coupling, and the objective function is expressed as: in, is the discrete component in the x direction, is the discrete component in the y direction, is the discrete component in the z direction, λ is the Lagrange multiplier, η cross gradient factor, I is the identity matrix.

6. The method of three-dimensional joint inversion of gravity, magneto-electricity based on reweighting according to claim 5, characterized in that: The specific steps of S33 are: 1) The cross gradient representation of different physical property model matrices is: Among them, t is the three-dimensional cross gradient vector, m1 and m2 are different physical property model matrices, is the gradient operator; In the three-dimensional cross gradient vector, the x, y, and z components are expressed as: 2) Discretize the three components of the three-dimensional cross gradient vector and express it as: Among them, D x 、D y 、D z are the first-order central difference operators in the x, y, and z directions respectively; 3) Adding the cross-gradient constraint to the regularized inversion, a new regularized inversion formula based on the cross-gradient constraint is obtained, which is expressed as: F=F d +λΦ m +ηΦ cg ; Among them, Φ d is the data fitting term, Φ m is the model fitting term, Φ cg is the cross gradient term; 4) Substitute the three components discretized in step 2) into the regularized inversion formula in step 3) to obtain the objective function, which is expressed as: 5) Derivative the objective function with respect to the model parameters is expressed as:

7. The method of three-dimensional joint inversion of gravity, magneto-electricity based on reweighting according to claim 1, characterized in that: In S4, forward calculation is performed on the updated model, and the specific steps are as follows: 1) Calculate the gravity potential and gravity anomaly Δg for the updated density model, expressed as: Where V is the gravity potential, Δg is the gravity anomaly, G is the gravitational constant, and r is the distance between the volume element Q(x,y,z) and the observation point P(X,Y,Z); After discretization, we get: where, Δx i = X - x i ; Δy i = Y - y i ; Δz i = Z - z i ; 2) Calculate the magnetic potential and magnetic anomaly of the updated magnetic susceptibility model. When the magnetization intensity is M, the Poisson formula between the gravity potential and the magnetic potential is expressed as: Among them, M z 、M y 、M z are the components of the magnetization intensity in the x, y, and z directions, respectively, i is the magnetization inclination, and δ is the magnetization deflection; The component of the magnetic anomaly in the direction of the Earth's magnetic field is expressed as: ΔT=B x so-so-A+B y cosI′sinA′+B z sinI′; Among them, B x is the magnetic induction intensity in the x direction, B y is the magnetic induction intensity in the y direction, B z is the magnetic induction intensity in the z direction, I′ is the geomagnetic inclination, and A′ is the angle between the magnetic north and the x axis; 3) Forward calculation is performed on the updated magnetotelluric resistivity model. The formulas for apparent resistivity and impedance phase are expressed as follows: Among them, ρ ij is the apparent resistivity, φ ij is the impedance phase, ω is the angular frequency of the electromagnetic wave, Z ij is the wave impedance.

8. The method of three-dimensional joint inversion of gravity, magneto-electricity based on reweighting according to claim 1, characterized in that: In S5, during the iterative process of the three-dimensional joint inversion, when the model is updated, the model fitting status is evaluated by calculating the root mean square error based on the pre-input iteration constraint conditions. When the model fitting error enters the preset ideal error range or reaches the maximum allowable number of iterations, the iteration process is stopped; otherwise, the iterative search is continued until the predetermined error range requirements are met.

Images