Joint inversion of magnetotelluric and magnetic data with adaptive dynamic prior information

Abstract

The fusion adaptive dynamic prior information magnetotelluric and magnetic method joint inversion method belongs to the field of geophysical data processing and inversion technology; by constructing a joint inversion objective function including data fitting term, model constraint term and physical coupling term, and using Gauss-Newton method for solving, the core innovation lies in introducing adaptive dynamic prior information extraction technology; the method can automatically generate and optimize prior constraints from data itself and intermediate results, without relying on external preset information, thereby significantly improving the accuracy and reliability of the boundary description and physical property recovery of underground abnormal bodies (such as resistivity, magnetic susceptibility anomaly); solve the problem that the traditional physical property coupling joint inversion technology has strong dependence on prior information and weak application generalization ability in information missing area.

Description

Technical Field

[0001] This invention relates to the field of geophysical data processing and inversion technology, and in particular to a joint inversion method of magnetotelluric and magnetic methods that integrates adaptive dynamic prior information, applicable to joint imaging and anomaly identification of multiple subsurface physical properties (such as resistivity and magnetic susceptibility). Background Technology

[0002] Magnetotelluric sounding and magnetic exploration are two important geophysical exploration methods, sensitive to the resistivity and magnetic characteristics of subsurface media, respectively. However, inversion based on a single method often suffers from problems such as high ambiguity and limited resolution. Joint inversion technology can comprehensively utilize complementary information from two data sources, effectively reducing the ambiguity of the inversion and thus improving the reliability of subsurface structure imaging. Nevertheless, joint inversion still faces two core challenges: first, how to effectively couple the intrinsic relationships between different physical property parameters; and second, how to achieve stability in the inversion process when prior information is lacking. Traditional physical property coupling joint inversion methods usually require pre-setting empirical relationships or cluster centers between physical properties as prior constraints. This strong dependence on prior information is difficult to apply in areas with low exploration levels and scarce known information. Furthermore, accurate identification of anomalies directly determines the reliability of the final inversion results. Therefore, there is an urgent need for an efficient two-dimensional joint inversion method that can adaptively extract prior information and dynamically optimize the inversion process. Summary of the Invention

[0003] The purpose of this invention is to provide a joint inversion method for magnetotellurics and magnetism that integrates adaptive dynamic prior information, so as to solve the problems of existing joint inversion techniques being highly dependent on prior information and having weak generalization ability under conditions of missing information, while achieving efficient and high-precision two-dimensional inversion.

[0004] A joint inversion method combining magnetotelluric and magnetic methods that integrates adaptive dynamic prior information, such as Figure 1 As shown, it includes the following steps:

[0005] Step 1: Input magnetotelluric observation data and magnetic observation data, and set the initial two-dimensional resistivity model and magnetic susceptibility model;

[0006] Step 2: Construct a joint inversion objective function for magnetotelluric and magnetic data based on rock physics constraints. The objective function includes a magnetotelluric data fitting term, a magnetic data fitting term, a model smoothing constraint term, and a rock physics coupling term based on clustering information constraints.

[0007] Step 3: Solve the joint inversion objective function using the Gauss-Newton method to obtain the resistivity and magnetic susceptibility model update values ​​for the current iteration, and then obtain the joint inversion results of the current iteration of magnetotelluric and magnetic methods;

[0008] Step 4: After each joint inversion iteration, perform an adaptive dynamic prior information extraction step. From the current iteration joint inversion results, use an adaptive dynamic prior information extraction technique based on multi-feature fusion and consistency screening to automatically extract the cluster centers and number of clusters in the resistivity and magnetic susceptibility models as the rock physics prior information for the next iteration.

[0009] Step 5: Substitute the extracted cluster centers and the number of clusters into the objective function of the next joint inversion iteration to achieve dynamic closed-loop feedback and update of prior rock physics information;

[0010] Step 6: Determine whether the preset convergence criterion has been met based on the data fitting difference of the current joint inversion resistivity and magnetic susceptibility model. If not, repeat steps 3, 4 and 5 above until the fitting difference requirement is met or the maximum number of iterations is reached.

[0011] In step four, the adaptive dynamic prior information extraction step includes:

[0012] S1: Based on the inversion model obtained in step three, the comprehensive anomaly score of each grid cell is calculated by multi-feature weighted fusion, and multiple sets of weights are combined to enhance the robustness of identification.

[0013] S2: Adaptively determine the optimal number of clusters using the elbow rule, and perform K-means clustering on the abnormal score set;

[0014] S3: Integrate the clustering results under multiple weight combinations, calculate the consistency score of each grid cell, and screen out high-confidence outlier regions accordingly;

[0015] S4: Cluster the highly consistent regions again, extract the centroids of each anomaly, and use them as the prior cluster centers of the attribute parameter space;

[0016] S5: Simultaneously, iteratively update background statistics, including the mean and standard deviation, for anomaly measurement in the next iteration.

[0017] In S1, the calculation of the comprehensive anomaly score integrates spatial anomalies and attribute anomalies:

[0018] S1.1: Spatial anomalies are measured by the normalized Euclidean distance from the grid cells to the center of the model space, and the calculation formula is as follows:

[0019] (1)

[0020] in, For grid cells Spatial anomaly score, For grid cell coordinates, The coordinates of the center of the model space are This represents the total number of rows in the model's grid. This represents the total number of columns in the model grid.

[0021] S1.2: Attribute anomaly is measured by the absolute value of the standard score (Z-score) to determine the degree to which a property parameter value deviates from the global average level. The calculation formula is as follows:

[0022] (2)

[0023] in, For grid cells The attribute abnormality score, These are the attribute parameter values ​​(such as resistivity and magnetic susceptibility) of the grid cell. and These are the mean and standard deviation of the parameters in the entire model, respectively.

[0024] S1.3: The overall abnormal score is adjusted by adjusting the weighting coefficient. and The generation and calculation formula is:

[0025] (3)

[0026] in, For grid cells The overall abnormal score, For spatial anomalies, For attribute anomalies, These are the weighting coefficients for spatial anomalies. The weight coefficients for attribute anomalies satisfy the following conditions: By setting multiple weight combinations Multiple abnormal score sets were obtained. .

[0027] In S5, the iterative update of background statistics includes the following steps:

[0028] S5.1: Initialization: Use the mean and standard deviation of the global model as the initial values ​​for the background statistics;

[0029] S5.2: Based on the current background statistics, calculate the comprehensive anomaly score and use a consensus voting strategy with multiple weights to identify high-confidence anomaly regions;

[0030] S5.3: Use areas not identified as anomalies as the current background and recalculate the background mean and standard deviation;

[0031] S5.4: Determine whether the relative change of the background mean or standard deviation is less than the preset convergence threshold. If so, stop the iteration; otherwise, continue the iteration with the new background statistic.

[0032] S5.5: During the iteration process, the consistency threshold is adaptively adjusted to gradually eliminate the interference of outliers on the background estimation.

[0033] The joint inversion objective function is expressed as follows:

[0034] (4)

[0035] in, For the joint inversion objective function value, This is the resistivity model parameter vector. This is the parameter vector for the magnetic susceptibility model. For including the first The resistivity and magnetic susceptibility vectors of each geological unit. This is the weight matrix for magnetotelluric or magnetic observation data. This is a vector of magnetotelluric or magnetic observation data. For the orthogonal operator, For regularization parameters, The smoothing matrix of the inversion model. As the reference model for inversion, The total regularization parameter for the rock physics coupling term is... This represents the total number of mesh elements in the model. The number of clusters, For elements of the fuzzy membership matrix, Indicates the first A vector composed of the resistivity and magnetic susceptibility model parameters of each grid cell. For the first The prior constraint weights for each cluster are vectors. Including the The resistivity and magnetization vectors of the prior cluster centers.

[0036] In step three, when solving the objective function using the Gauss-Newton method, the model update amount is calculated using the following formula:

[0037] (5)

[0038] in, For the first The amount of resistivity or magnetic susceptibility model update in each iteration. This is the Jacobian matrix of the magnetotelluric or magnetic model. This is the weight matrix for magnetotelluric or magnetic observation data. Let be the regularization parameter for the i-th iteration. The smoothing matrix of the inversion model. For the first The total regularization parameter of the rock physics coupling term at the next iteration. The number of clusters, For the first During the iteration, all grid cells are paired with the... A vector composed of the membership degrees of each cluster. This is a vector of magnetotelluric or magnetic observation data. For the first Forward operators of the next iteration model For the first The inversion model parameter vector obtained from the next iteration. This serves as the reference model for the inversion.

[0039] The beneficial effects of this invention are:

[0040] 1. The adaptive dynamic prior extraction technology proposed in this invention can automatically extract attribute cluster centers from the current model iteration results during the inversion process without relying on externally preset prior information, which significantly improves the generalization ability of this method in areas lacking exploration information.

[0041] 2. This invention improves the robustness and accuracy of anomaly region identification through multi-feature weighted fusion and consistency screening strategies. The iterative optimization mechanism of background statistics effectively eliminates the contamination of background field estimation by anomalies, making the extracted prior information more reliable. Attached Figure Description

[0042] Figure 1 This is a flowchart of the magnetotelluric and magnetic method joint inversion method that integrates adaptive dynamic prior information, as described in this invention.

[0043] Figure 2 This is a schematic diagram of the inversion results of the theoretical model of this invention.

[0044] Figure 3 The diagrams show four sets of K-means clustering and elbow rule analysis of resistivity models with different weight combinations after the inversion model is obtained in this invention. The left side shows the K-means clustering results, and the right side shows the clustering number analysis based on the elbow rule.

[0045] Figure 4 This is a diagram showing the consistency of outlier detection and the final clustering of the resistivity model under a multi-weight combination according to an embodiment of the present invention. The left side is the consistency diagram of outlier detection of the resistivity model under a multi-weight combination, and the right side is the final clustering diagram of the resistivity model.

[0046] Figure 5The diagrams show the K-means clustering and elbow rule analysis of the magnetic susceptibility model obtained by this invention after the inversion model is obtained. The left side shows the K-means clustering results, and the right side shows the clustering number analysis based on the elbow rule.

[0047] Figure 6 This diagram shows the consistency of outlier detection and the final clustering diagram of the magnetic susceptibility model under multi-weighted combinations according to an embodiment of the present invention. The left side shows the consistency diagram of outlier detection in the magnetic susceptibility model under multi-weighted combinations, and the right side shows the final clustering diagram of the magnetic susceptibility model. Detailed Implementation

[0048] Please see Figures 1 to 6 The image shown is an embodiment of the present invention.

[0049] The joint inversion method of magnetotelluric and magnetic methods, which integrates adaptive dynamic prior information, includes the following steps:

[0050] Step 1: Input magnetotelluric observation data and magnetic observation data, and set the initial two-dimensional resistivity model and magnetic susceptibility model.

[0051] Step 2: Construct a joint inversion objective function for magnetotelluric and magnetic data based on rock physics constraints. The objective function includes a magnetotelluric data fitting term, a magnetic data fitting term, a model smoothing constraint term, and a rock physics coupling term based on clustering information constraints.

[0052] Step 3: Solve the joint inversion objective function using the Gauss-Newton method to obtain the resistivity and magnetic susceptibility model update values ​​for the current iteration, and then obtain the joint inversion results of the current iteration of magnetotelluric and magnetic methods.

[0053] In step three, when solving the objective function using the Gauss-Newton method, the model update amount is calculated using the following formula:

[0054] (5)

[0055] in, For the first The amount of resistivity or magnetic susceptibility model update in each iteration. This is the Jacobian matrix of the magnetotelluric or magnetic model. This is the weight matrix for magnetotelluric or magnetic observation data. For the first The regularization parameter at the next iteration The smoothing matrix of the inversion model. For the first The total regularization parameter of the rock physics coupling term at the next iteration. The number of clusters, For the first During the iteration, all grid cells are paired with the... A vector composed of the membership degrees of each cluster. This is a vector of magnetotelluric or magnetic observation data. For the first Forward operators of the next iteration model For the first The inversion model parameter vector obtained from the next iteration. This serves as the reference model for the inversion.

[0056] Step 4: After each joint inversion iteration, perform an adaptive dynamic prior information extraction step. From the current iteration joint inversion results, use an adaptive dynamic prior information extraction technique based on multi-feature fusion and consistency screening to automatically extract the cluster centers and numbers in the resistivity and magnetic susceptibility models as the rock physics prior information for the next iteration.

[0057] In step four, the adaptive dynamic prior information extraction step includes:

[0058] S1: Based on the inversion model obtained in step three, the comprehensive anomaly score of each grid cell is calculated by multi-feature weighted fusion, and multiple sets of weights are combined to enhance the robustness of identification.

[0059] S2: Adaptively determine the optimal number of clusters using the elbow rule, and perform K-means clustering on the abnormal score set;

[0060] S3: Integrate the clustering results under multiple weight combinations, calculate the consistency score of each grid cell, and screen out high-confidence outlier regions accordingly;

[0061] S4: Cluster the highly consistent regions again, extract the centroids of each anomaly, and use them as the prior cluster centers of the attribute parameter space;

[0062] S5: Simultaneously, iteratively update background statistics, including the mean and standard deviation, for anomaly measurement in the next iteration.

[0063] In S1, the calculation of the comprehensive anomaly score integrates spatial anomalies and attribute anomalies:

[0064] S1.1: Spatial anomalies are measured by the normalized Euclidean distance from the grid cells to the center of the model space, and the calculation formula is as follows:

[0065] (1)

[0066] in, For grid cells Spatial anomaly score, ( i, j ) represents the coordinates of the grid cells, ( i c, j c () represents the coordinates of the center of the model space. This represents the total number of rows in the model's grid. This represents the total number of columns in the model grid.

[0067] S1.2: Attribute anomaly is measured by the absolute value of the standard score (Z-score) to determine the degree to which a property parameter value deviates from the global average level. The calculation formula is as follows:

[0068] (2)

[0069] in, For grid cells The attribute abnormality score, These are the attribute parameter values ​​(such as resistivity and magnetic susceptibility) of the grid cell. and These are the mean and standard deviation of the parameters in the entire model, respectively.

[0070] S1.3: The overall abnormal score is adjusted by adjusting the weighting coefficient. and The generation and calculation formula is:

[0071] (3)

[0072] in, For grid cells The overall abnormal score, For spatial anomalies, For attribute anomalies, These are the weighting coefficients for spatial anomalies. The weight coefficients for attribute anomalies satisfy the following conditions: By setting multiple weight combinations Multiple abnormal score sets were obtained. .

[0073] In S5, the iterative update of background statistics includes the following steps:

[0074] S5.1: Initialization: Use the mean and standard deviation of the global model as the initial values ​​for the background statistics;

[0075] S5.2: Based on the current background statistics, calculate the comprehensive anomaly score and use a consensus voting strategy with multiple weights to identify high-confidence anomaly regions;

[0076] S5.3: Use areas not identified as anomalies as the current background and recalculate the background mean and standard deviation;

[0077] S5.4: Determine whether the relative change of the background mean or standard deviation is less than the preset convergence threshold. If so, stop the iteration; otherwise, continue the iteration with the new background statistic.

[0078] S5.5: During the iteration process, the consistency threshold is adaptively adjusted to gradually eliminate the interference of outliers on the background estimation.

[0079] This invention uses guided fuzzy C-means clustering as a coupling constraint to construct the following joint inversion objective function:

[0080] (4)

[0081] in, For the joint inversion objective function value, This is the resistivity model parameter vector. This is the parameter vector for the magnetic susceptibility model. For including the first The resistivity and magnetic susceptibility vectors of each geological unit. This is the weight matrix for magnetotelluric or magnetic observation data. This is a vector of magnetotelluric or magnetic observation data. For the orthogonal operator, For regularization parameters, The smoothing matrix of the inversion model. As the reference model for inversion, The total regularization parameter for the rock physics coupling term is... This represents the total number of mesh elements in the model. The number of clusters, For elements of the fuzzy membership matrix, Indicates the first A vector composed of the resistivity and magnetic susceptibility model parameters of each grid cell. For the first The prior constraint weights for each cluster are vectors. Including the The resistivity and magnetization vectors of the prior cluster centers.

[0082] vector Represents observed geophysical data. This is specified here. For magnetotelluric (MT) data, This is magnetic data. The corresponding model parameters. and These are resistivity and magnetic susceptibility, respectively.

[0083] To minimize the objective function (4), the Gauss-Newton method is used. Assume that the superscript on any variable... This indicates that the variable is in the first... The value at the next iteration, then and They are the first The model parameter estimates and forward modeling response at the nth iteration. Then, at the 1st iteration... In the +1 iteration, the model update amount is obtained by minimizing the following objective function:

[0084] (6)

[0085] in, For the joint inversion objective function value, For the first The amount of resistivity or magnetic susceptibility model update in each iteration. It is the Jacobian matrix of MT and magnetic data. This is the weight matrix for magnetotelluric or magnetic observation data. For regularization parameters, The smoothing matrix of the inversion model. For the overall regularization parameter of the rock physics coupling term, The number of clusters, For the first During the iteration, all grid cells are paired with the... A vector composed of the membership degrees of each cluster. This is a vector of magnetotelluric or magnetic observation data. For the first Forward operators of the next iteration model For the first The inversion model parameter vector obtained from the next iteration. As the reference model for inversion, For the first The iteration contains the first... The resistivity and magnetic susceptibility vectors of each geological unit. Including the The resistivity and magnetization vectors of the prior cluster centers.

[0086] Step 5: Substitute the extracted cluster centers and the number of clusters into the objective function of the next joint inversion iteration to achieve dynamic closed-loop feedback and update of prior rock physics information.

[0087] Step 6: Determine whether the preset convergence criterion has been met based on the data fitting difference of the current joint inversion resistivity and magnetic susceptibility model. If not, repeat steps 3, 4 and 5 above until the fitting difference requirement is met or the maximum number of iterations is reached.

[0088] The final model update amount obtained from the solution is:

[0089] (5)

[0090] in, For the first The amount of resistivity or magnetic susceptibility model update in each iteration. This is the Jacobian matrix of the magnetotelluric or magnetic model. This is the weight matrix for magnetotelluric or magnetic observation data. For the first The regularization parameter at the next iteration The smoothing matrix of the inversion model. For the first The total regularization parameter of the rock physics coupling term at the next iteration. The number of clusters, For the first During the iteration, all grid cells are paired with the... A vector composed of the membership degrees of each cluster. This is a vector of magnetotelluric or magnetic observation data. For the first Forward operators of the next iteration model For the first The inversion model parameter vector obtained from the next iteration. This serves as the reference model for the inversion.

[0091] To reduce reliance on prior information, this invention further proposes an adaptive dynamic prior extraction technique. This technique automatically extracts and updates the feature centers required for attribute clustering from the current inversion result in each inversion iteration, achieving adaptive optimization of prior information. This innovation effectively overcomes the strong dependence on prior information in traditional property coupling joint inversion methods, improving the application generalization ability under conditions of missing information. Through this dynamic feedback mechanism, the system can automatically adjust the prior constraints of the inversion in each iteration, optimizing the recovery process of multiple property parameters.

[0092] The dynamic prior extraction technique proposed in this invention is based on an adaptive anomaly identification algorithm that automatically extracts the prior information required for attribute coupling from the attribute parameter model generated in each inversion iteration. This technique is grounded in statistical pattern recognition and unsupervised machine learning principles, and its theoretical framework and implementation process are as follows:

[0093] Multi-feature weighted fusion and anomaly score calculation:

[0094] Calculate the overall anomaly score for each grid cell, which combines spatial anomalies and attribute anomalies:

[0095] Spatial anomaly metric: Calculate the normalized Euclidean distance from each grid cell to the center of the model space to identify cells located at the model edge or in special structural positions.

[0096] (1)

[0097] in, For grid cells Spatial anomaly score, ( i, j ) represents the coordinates of the grid cells, ( i c , j c () represents the coordinates of the center of the model space. This represents the total number of rows in the model's grid. This represents the total number of columns in the model mesh; the denominator is used for normalization, making... .

[0098] Property anomaly measure: The absolute value of the standard score (Z-score) is used to measure the degree to which the property parameter value deviates from the global average level.

[0099] (2)

[0100] For grid cells The attribute abnormality score, This refers to the attribute parameter value of this mesh cell. and These are the mean and standard deviation of the parameters in the entire model, respectively.

[0101] Overall anomaly score: adjusted by weighting coefficients and Multiple sets of anomaly graphs are generated to enhance robustness.

[0102] (3)

[0103] in, For grid cells The overall abnormal score, For spatial anomalies, For attribute anomalies, These are the weighting coefficients for spatial anomalies. The weight coefficients for attribute anomalies satisfy the following conditions: By setting multiple weight combinations Multiple abnormal score sets were obtained. .

[0104] Adaptive clustering analysis based on the elbow rule:

[0105] For each abnormal score set Perform K-means clustering and automatically determine the optimal number of clusters using the elbow rule. .

[0106] The elbow rule identifies the inflection point by calculating the within-cluster sum of squares (WSS) for different K values:

[0107] (7)

[0108] in, It is the first Clusters, It is the centroid of the cluster. Optimal cluster number. It is the K value corresponding to the maximum slope change in the WSS curve. This method can achieve adaptive classification of anomalies of different shapes and sizes.

[0109] High-consistency region selection and prior information generation:

[0110] By integrating clustering results from multiple weight combinations, a consistency score is calculated for each grid cell to be identified as an anomaly:

[0111] (8)

[0112] in, Each grid cell is identified as having an anomalous consistency score. As an indicator function, when the grid cell In the The value is 1 when a group is identified as an anomaly, and 0 otherwise. For grid cell coordinates, For abnormal score sets, This represents the total number of weighted combinations.

[0113] Subsequently, regions identified as anomalous under most weight combinations (i.e., regions with high consistency scores) are retained, forming high-confidence anomalous regions. These high-consistency anomalous regions are then subjected to final clustering, and the centroid of each cluster is extracted. The centroid encapsulates the comprehensive characteristics of the anomalous body and can be transformed into a priori cluster center in the attribute parameter space. This is then input into the property coupling term of the next inversion iteration.

[0114] Iterative optimization of background statistics:

[0115] To improve the accuracy of anomaly identification, the algorithm introduces an internal iterative optimization mechanism. In each main inversion iteration, the background statistics are updated multiple times based on the current inversion result until convergence. The specific process is as follows:

[0116] Initialization: The first iteration uses the mean and standard deviation of the global model as background statistics.

[0117] Iterative loop:

[0118] 1. Based on the current background statistics, calculate the spatial anomaly and attribute anomaly of each grid cell, and generate a comprehensive anomaly score through multiple sets of weights.

[0119] 2. Employ a consensus voting strategy with multiple weights to identify high-confidence outlier regions (consistency scores exceeding an adaptive threshold).

[0120] 3. Use the regions that were not identified as abnormal (i.e., non-abnormal regions) as the current background and recalculate the background mean and standard deviation.

[0121] 4. If the relative change in the background mean or standard deviation is less than the preset convergence threshold (e.g., 2%), then stop the iteration; otherwise, return to step 1 with the new background statistics.

[0122] Adaptive threshold adjustment: In the early stages of iteration, a more lenient consistency threshold (e.g., 0.4) is used to accommodate potential anomalies; in the later stages, the threshold is increased (e.g., 0.6) to ensure high confidence in the abnormal region, thereby gradually eliminating the interference of anomalies on the background estimation.

[0123] This iterative optimization process effectively overcomes the contamination of background statistics by anomalies, enabling the background estimate to gradually approximate the true geophysical background field and providing a more reliable benchmark for subsequent anomaly feature extraction. The converged background mean and standard deviation will be used for the final connectivity analysis and feature center extraction of the anomaly region.

[0124] The following section verifies the magnetotelluric and magnetic method joint inversion method that integrates adaptive dynamic prior information extraction technology described in this invention.

[0125] To verify the effectiveness and process of this method, a two-dimensional synthetic model containing two typical anomalies was designed. The background region of the model is 6km × 3km in size, with a background dielectric resistivity of 100Ω·m and a magnetic susceptibility of 0.03SI. The model contains two target anomalies ( Figure 2 (a and b in the text): One is a cuboid located in the upper left corner with a resistivity of 1000 Ω·m and a magnetic susceptibility of 0.01 SI; the other is a stepped anomalous body located in the lower right corner with a resistivity of 10 Ω·m and a magnetic susceptibility of 0.05 SI.

[0126] A total of 29 uniformly distributed AMT and magnetic measurement points were established on the surface of the study area, with a spacing of 200 m between points. In the forward modeling simulation, 10 logarithmically spaced frequency points were selected for the AMT method, ranging from 1 Hz to 1000 Hz; the magnetic data consisted of the total magnetic anomaly. The inversion region was divided into 70 × 30 grid cells, with an extended grid added to the periphery to eliminate boundary effects. All inversion calculations used a uniform half-space background model (resistivity 100 Ω·m, magnetic susceptibility 0.03 SI) as the initial model.

[0127] First, the traditional smoothing-constrained joint inversion method was used to invert the magnetotelluric and magnetic data. After eight iterations, the inversion fitting difference tended to stabilize, ultimately settling at 0.446271 and 0.038164, respectively. The final resistivity and magnetic susceptibility models are as follows: Figure 2 As shown in c and d in the figure. The inversion results show that the traditional smooth constraint inversion method can roughly recover the spatial position of the two anomalies, but the characterization of the boundary of the anomalies is relatively blurry, with obvious divergence. The shape and boundary are not clearly reconstructed, and the recovered property amplitude also has a large gap with the real model.

[0128] Next, under identical observation data and initial model conditions, the adaptive dynamic prior constraint joint inversion method proposed in this invention is applied. Based on the current inversion model, this technique automatically identifies high-confidence anomalous regions through multi-feature fusion and cluster analysis, and uses the extracted anomalous body feature centers as prior information for property coupling, feeding them back into the next joint inversion iteration. After 9 iterations, the inversion finally converged, with fitting differences of 0.268075 and 0.130785, respectively. The obtained resistivity and magnetic susceptibility models are as follows: Figure 2 As shown in e and f.

[0129] A comparison of the inversion results from the two methods clearly shows that the model reconstructed by the method of this invention has a significant improvement in spatial resolution. The boundaries of the two anomalies are clearer and sharper, and they better match the contours of the real model. At the same time, the physical property amplitudes of the anomalies are also closer to the theoretical model.

[0130] To illustrate and verify the key steps and effectiveness of adaptive prior information extraction in the inversion process. Figure 3 and Figure 5 The resistivity and magnetic susceptibility inversion models based on the eighth iteration are presented respectively, using multiple sets ( α , β The process of calculating the weighted average of outlier scores and then performing K-means clustering. Figure 3 and Figure 5 The cluster number analysis graph based on the elbow rule on the right shows the elbow rule curve for the corresponding number of clusters K. Its inflection point (red asterisk) clearly indicates that the optimal number of clusters K=2, which is consistent with the number of real outliers in the model, proving the reliability of the adaptive cluster number determination strategy.

[0131] Figure 4 and Figure 6The left-hand side shows the outlier detection consistency map of the resistivity model under multi-weighted combinations, illustrating the distribution of "outlier detection consistency" in the clustering results under these combinations. Colors ranging from black to brown represent consistency scores, with high-consistency regions clearly outlining the spatial locations of the two outliers. The "final connected region analysis" in the right-hand side of the inversion model's final clustering map shows that after final clustering of the high-consistency regions, the two outliers were accurately separated and identified.

[0132] The anomalous mass centers extracted by the above process are dynamically transformed into prior cluster centers in the physical property parameter space. The data is then input into the coupling constraint terms of the next inversion iteration. This dynamic closed-loop mechanism of "inversion-extraction-feedback" enables the method of this invention to extract precise structural constraints step by step from the data itself and intermediate results of the iteration without any external prior information, which significantly improves the recovery ability of complex morphological anomalies and the reliability of the inversion results.

Claims

1. A joint inversion method for magnetotelluric and magnetic methods that integrates adaptive dynamic prior information, characterized in that, Includes the following steps: Step 1: Input magnetotelluric observation data and magnetic observation data, and set the initial two-dimensional resistivity model and magnetic susceptibility model; Step 2: Construct a joint inversion objective function for magnetotelluric and magnetic data based on rock physics constraints. The objective function includes a magnetotelluric data fitting term, a magnetic data fitting term, a model smoothing constraint term, and a rock physics coupling term based on clustering information constraints. Step 3: Solve the joint inversion objective function using the Gauss-Newton method to obtain the resistivity and magnetic susceptibility model update values ​​for the current iteration, and then obtain the joint inversion results of the current iteration of magnetotelluric and magnetic methods; Step 4: After each joint inversion iteration, perform an adaptive dynamic prior information extraction step. From the current iteration joint inversion results, use an adaptive dynamic prior information extraction technique based on multi-feature fusion and consistency screening to automatically extract the cluster centers and number of clusters in the resistivity and magnetic susceptibility models as the rock physics prior information for the next iteration. The adaptive dynamic prior information extraction steps include: S1: Based on the inversion model obtained in step three, the comprehensive anomaly score of each grid cell is calculated by multi-feature weighted fusion, and multiple sets of weights are combined to enhance the robustness of identification. S2: Adaptively determine the optimal number of clusters using the elbow rule, and perform K-means clustering on the abnormal score set; S3: Integrate the clustering results under multiple weight combinations, calculate the consistency score of each grid cell, and screen out high-confidence outlier regions accordingly; S4: Cluster the highly consistent regions again, extract the centroids of each anomaly, and use them as the prior cluster centers of the attribute parameter space; S5: Simultaneously, iteratively update background statistics, including the mean and standard deviation, for anomaly measurement in the next iteration; Step 5: Substitute the extracted cluster centers and the number of clusters into the objective function of the next joint inversion iteration to achieve dynamic closed-loop feedback and update of prior rock physics information; Step 6: Determine whether the preset convergence criterion has been met based on the data fitting difference of the current joint inversion resistivity and magnetic susceptibility model. If not, repeat steps 3, 4 and 5 above until the fitting difference requirement is met or the maximum number of iterations is reached. In S1, the calculation of the comprehensive anomaly score integrates spatial anomalies and attribute anomalies: S1.1: Spatial anomalies are measured by the normalized Euclidean distance from the grid cells to the center of the model space, and the calculation formula is as follows: （1） in, For grid cells Spatial anomaly score, ( i, j ) represents the coordinates of the grid cells, ( i c , j c () represents the coordinates of the center of the model space. This represents the total number of rows in the model's grid. This represents the total number of columns in the model grid. S1.2: Attribute anomaly is measured by the absolute value of the standard score, indicating how much the property parameter value deviates from the global average level. The calculation formula is as follows: （2） in, For grid cells The attribute abnormality score, This refers to the attribute parameter value of this mesh cell. and These are the mean and standard deviation of the parameters in the entire model, respectively. S1.3: The overall abnormal score is adjusted by adjusting the weighting coefficient. and The generation and calculation formula is: （3） in, For grid cells The overall abnormal score, For spatial anomalies, For attribute anomalies, These are the weighting coefficients for spatial anomalies. The weight coefficients for attribute anomalies satisfy the following conditions: By setting multiple weight combinations Multiple abnormal score sets were obtained. .

2. The magnetotelluric and magnetic method joint inversion method integrating adaptive dynamic prior information according to claim 1, characterized in that: In S5, the iterative update of background statistics includes the following steps: S5.1: Initialization: Use the mean and standard deviation of the global model as the initial values ​​for the background statistics; S5.2: Based on the current background statistics, calculate the comprehensive anomaly score and use a consensus voting strategy with multiple weights to identify high-confidence anomaly regions; S5.3: Use areas not identified as anomalies as the current background and recalculate the background mean and standard deviation; S5.4: Determine whether the relative change of the background mean or standard deviation is less than the preset convergence threshold. If so, stop the iteration; otherwise, continue the iteration with the new background statistic. S5.5: During the iteration process, the consistency threshold is adaptively adjusted to gradually eliminate the interference of outliers on the background estimation.

3. The magnetotelluric and magnetic method joint inversion method according to claim 1, characterized in that: The joint inversion objective function is expressed as follows: （4） in, For the joint inversion objective function value, This is the resistivity model parameter vector. This is the parameter vector for the magnetic susceptibility model. For including the first The resistivity and magnetic susceptibility vectors of each geological unit. This is the weight matrix for magnetotelluric or magnetic observation data. This is a vector of magnetotelluric or magnetic observation data. For the orthogonal operator, For regularization parameters, The smoothing matrix of the inversion model. As the reference model for inversion, The total regularization parameter for the rock physics coupling term is... This represents the total number of mesh elements in the model. The number of clusters, For elements of the fuzzy membership matrix, Indicates the first A vector composed of the resistivity and magnetic susceptibility model parameters of each grid cell. For the first The prior constraint weights for each cluster are vectors. Including the The resistivity and magnetization vectors of the prior cluster centers.

4. The magnetotelluric and magnetic method joint inversion method fusing adaptive dynamic prior information according to claim 1, characterized in that: In step three, when solving the objective function using the Gauss-Newton method, the model update amount is calculated using the following formula: （5） in, For the first The amount of resistivity or magnetic susceptibility model update in each iteration. This is the Jacobian matrix of the magnetotelluric or magnetic model. This is the weight matrix for magnetotelluric or magnetic observation data. For the first The regularization parameter at the next iteration The smoothing matrix of the inversion model. For the first The total regularization parameter of the rock physics coupling term at the next iteration. The number of clusters, For the first During the iteration, all grid cells are paired with the... A vector composed of the membership degrees of each cluster. This is a vector of magnetotelluric or magnetic observation data. For the first Forward operators of the next iteration model For the first The inversion model parameter vector obtained from the next iteration. This serves as the reference model for the inversion.

Images