A magnetotelluric three-dimensional hybrid optimization inversion method

Abstract

The application discloses a magnetotelluric three-dimensional hybrid optimization inversion method, and steps include: updating a resistivity model; judging whether a nonlinear conjugate gradient iteration satisfies an exit condition; judging whether the gradient of the objective function of two adjacent nonlinear conjugate gradient iterations is less than a preset threshold; making the objective function jump out of a local optimal solution; judging whether a current particle swarm optimization trigger counter is greater than a maximum number of times; generating an updated resistivity model by using an optimal scaling factor; judging whether the RMS of the particle swarm optimization after optimization reaches an exit condition; and outputting a final resistivity model of a detection area. The application effectively solves the dimension disaster problem of the PSO algorithm applied to three-dimensional large-scale inversion.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration, specifically a three-dimensional hybrid optimization inversion method for magnetotelluric exploration. Background Technology

[0002] Magnetotellurics (MT), a geophysical technique that uses natural alternating electromagnetic fields to detect subsurface electrical structures, has been widely applied in mineral and oil and gas resource exploration and the study of deep crustal structures. The inversion process, which involves constructing a reasonable subsurface resistivity model from observed electromagnetic response data, is the core step in data interpretation.

[0003] However, existing inversion methods face the following serious challenges when dealing with complex geological structures:

[0004] Currently, mainstream inversion algorithms such as the Gauss-Newton method (GN) and the nonlinear conjugate gradient method (NLCG) are essentially local optimization processes. Although they are computationally efficient, in the presence of noise interference in the measured data, sparse sites, and deviations of the initial model from reality, the inversion is prone to getting trapped in local minima, resulting in blurred edges of anomalies and insufficient recovery of physical property amplitudes.

[0005] Due to limitations in observational scope, inversion methods relying solely on electromagnetic data suffer from severe ambiguity. While existing methods can enhance the structural consistency of resistivity models by introducing cross-gradient constraints and utilizing prior tectonic information such as seismic data, they still fail to address the challenge of co-optimizing structural constraints and data fitting in high-dimensional space.

[0006] When intelligent algorithms, such as Particle Swarm Optimization (PSO), are directly applied to large-scale 3D inversion, the computational cost becomes unbearable due to the exponential growth of the mesh parameters to be optimized, and premature convergence is prone to occur in high-dimensional search spaces.

[0007] Existing technologies struggle to effectively combine the efficiency of gradient methods with the robustness of global optimization algorithms. When dealing with complex tectonic regions with hidden fault zones, such as the Luxian earthquake zone, conventional methods often fail to accurately identify the fine boundaries of deep high- and low-resistivity anomalies, thus failing to provide interpretations with high geological reliability. Summary of the Invention

[0008] The purpose of this invention is to provide a three-dimensional hybrid optimization inversion method for magnetotelluric data, comprising the following steps:

[0009] Step 1) Set the initial resistivity model of the detection area Regularization parameters and Cooling factor and Maximum number of triggers for particle swarm optimization and convergence criteria;

[0010] Step 2) Using the objective function, which includes data fitting terms, model smoothing terms, and cross-gradient terms, as constraints, perform nonlinear conjugate gradient iteration to update the resistivity model, denoted as... ;

[0011] Among them, the model smoothing term is used to constrain the spatial smoothness of the resistivity model, and the cross gradient term is used to quantify the structural similarity between the resistivity model and the prior model.

[0012] Step 3) Determine whether the nonlinear conjugate gradient iteration meets the exit condition. If yes, jump to step 9; otherwise, proceed to step 4.

[0013] Step 4) Determine the gradient of the objective function in two adjacent nonlinear conjugate gradient iterations. If the value is less than a preset threshold, proceed to step 5; otherwise, return to step 2.

[0014] Step 5) Adjust the cooling factor attenuation regularization parameters according to the preset parameters. and Reduce the weight constraints of the model smoothing term and cross gradient term in the objective function, so that the objective function can escape the local optimum, and then proceed to step 6).

[0015] Step 6) Determine the current particle swarm optimization trigger counter Is it greater than the maximum number of times? If yes, return to step 2) and continue with the nonlinear conjugate gradient iteration; otherwise, proceed to step 7).

[0016] Step 7) For the current model Preprocessing and distance metric construction are performed, followed by spatiotemporal density clustering analysis to extract scaling factor sequences. ;

[0017] Perform particle swarm optimization in the reduced parameter space, from the scaling factor sequence. Obtain the optimal scaling factor Utilizing the optimal scaling factor Generate updated resistivity model ;

[0018] Step 8) Update the counter And the number of iterations Determine whether the RMS after particle swarm optimization meets the exit condition. If it does, proceed to step 9; otherwise, return to step 2.

[0019] Step 9) Output the final resistivity model of the probe area and terminate the inversion process.

[0020] Furthermore, the cross-gradient term is used to quantify the structural similarity between different physical property models, i.e.:

[0021] (1)

[0022] In the formula, For cross gradient terms; , This represents the gradient between the resistivity model and the prior model.

[0023] Furthermore, in step 2), the objective function of the nonlinear conjugate gradient iteration... As shown below:

[0024] (2)

[0025] In the formula, This represents the affine transformation form of the model parameter m. It is an initial uniform half-space model. This is the model covariance matrix; d represents the magnetotelluric observation data. It is a forward function; These are cross gradient components; The data covariance matrix; , and These correspond to the observed data, model parameters, and the total number of cross gradient components, respectively.

[0026] Furthermore, the gradient of the objective function As shown below:

[0027] (3)

[0028] (4)

[0029] In the formula, It is a forward response For model parameters Jacobian matrix; This is the partial derivative matrix of the cross gradient operator with respect to the model parameters; This refers to the data residuals.

[0030] Furthermore, in step 2), the updated resistivity model As shown below:

[0031] (5)

[0032] In the formula, Step size; Indicating the search direction; This serves as the initial resistivity model and a reference for transforming model parameters. and These are the resistivity model parameters in the transformed domain during the (k-1)th and kth iterations, respectively.

[0033] Furthermore, in step 7), the current resistivity model is... Preprocessing and distance metric construction are performed, followed by spatiotemporal density clustering analysis to extract scaling factors. The steps include:

[0034] Step 7.1) Based on background resistivity and preset resistivity threshold Resistivity model It is divided into low-resistivity subsets and high-resistivity subsets;

[0035] Step 7.2) Use the grid topology index (i, j, k) instead of the geometric center coordinates to calculate the distance;

[0036] The topological spatial distance between two grid cells p and q As shown below:

[0037] (6)

[0038] Attribute distance between two grid cells p and q As shown below:

[0039] (7)

[0040] Step 7.3) Based on the topological spatial distance dist1 and attribute distance dist2, perform clustering independently on the low-resistivity subset and the high-resistivity subset respectively, and divide the preprocessed model unit into n clusters. Each cluster corresponds to a set of resistivity values. and a scaling factor Thus constructing the scaling factor vector This is used for subsequent particle swarm optimization search.

[0041] Furthermore, when performing particle swarm optimization, the objective function As shown below:

[0042] (8)

[0043] In the formula, This is the resistivity model obtained from the current NLCG iteration; The n-dimensional scaling factor vector generated by spatiotemporal density clustering; This represents the scaling factor vector x on the model. The forward response of the model obtained by scaling the resistivity of each cluster.

[0044] Furthermore, when performing particle swarm optimization, each particle in the particle swarm represents a set of scaling factor vectors;

[0045] The particle state is updated as follows:

[0046] (9)

[0047] (10)

[0048] In the formula, , Let be the particle velocities at time t+1 and time t; , The positions of the particles at time t+1 and time t; , Individual learning factors and social learning factors; , The result is a random number (between 0 and 1). , For individual optimal position and global optimal position;

[0049] Updated resistivity model As shown below:

[0050] (11)

[0051] In the formula, The optimal scaling factor; It is a set of resistivity values.

[0052] Furthermore, the exit conditions for nonlinear conjugate gradient iteration are: the root mean square error (RMS) reaches a preset target threshold, the cooling factor of the regularization parameter decays to a preset lower limit, or the total number of inversion iterations reaches a set upper limit.

[0053] Furthermore, the regularization parameter is updated via a cooling factor, i.e.:

[0054] (12)

[0055] (13)

[0056] In the formula, the cooling factor , ; , These are the updated regularization parameters.

[0057] The technical effects of this invention are undeniable. This invention integrates local optimization and global search mechanisms, and its NLCG local inversion stage uses a cost function that includes cross-gradient constraints ( Ensure the model structure is consistent with the prior and converges efficiently; when trapped in a local extremum, trigger the PSO global optimization phase, using a function... The goal is to guide the search to escape local optima. Simultaneously, the ST-DBSCAN density clustering algorithm is introduced to reduce the dimensionality of the high-dimensional model parameters across the entire space to a small number of cluster scaling factors as PSO optimization variables, thereby significantly reducing the computational burden and effectively solving the curse of dimensionality problem when the PSO algorithm is applied to large-scale 3D inversion. Attached Figure Description

[0058] Figure 1 This is a schematic diagram of the rectangular double anomaly body and velocity prior model in Embodiment 12 of the present invention;

[0059] Figure 2 This is a schematic diagram of the cross-sections of conventional inversion and constrained inversion at z = 10 km and x = -1.25 km during the first inversion cooling of Embodiment 12 of the present invention.

[0060] Figure 3 This is a clustering result diagram of Embodiment 12 of the present invention;

[0061] Figure 4 This is a slice diagram of the inversion results of Embodiment 12 of the present invention;

[0062] Figure 5 This is a statistical graph of resistivity parameters from the inversion results of Embodiment 12 of the present invention;

[0063] Figure 6 The RMS and cross gradient value variation curves of the inversion results in Embodiment 12 of the present invention are shown.

[0064] Figure 7 This is a schematic diagram of the stepped double anomaly body and velocity prior model in Embodiment 13 of the present invention;

[0065] Figure 8 This is a slice diagram of the inversion results of Embodiment 13 of the present invention;

[0066] Figure 9 This is a statistical graph of resistivity parameters from the inversion results of Embodiment 13 of the present invention;

[0067] Figure 10 The RMS and cross gradient value variation curves of the inversion results in Embodiment 13 of the present invention are shown.

[0068] Figure 11 This is a flowchart of the method. Detailed Implementation

[0069] The present invention will be further described below with reference to embodiments, but it should not be construed that the scope of the present invention is limited to the following embodiments. Various substitutions and modifications made based on ordinary technical knowledge and common practices in the art without departing from the above-described technical concept of the present invention should be included within the scope of protection of the present invention.

[0070] Example 1:

[0071] See Figure 1 and Figure 11 A three-dimensional hybrid optimization inversion method for magnetotelluric data includes the following steps:

[0072] Step 1) Set the initial resistivity model of the detection area Regularization parameters and Cooling factor and Maximum number of triggers for particle swarm optimization and convergence criteria;

[0073] Step 2) Using the objective function, which includes data fitting terms, model smoothing terms, and cross-gradient terms, as constraints, perform nonlinear conjugate gradient iteration to update the resistivity model, denoted as... ;

[0074] Among them, the model smoothing term is used to constrain the spatial smoothness of the resistivity model, and the cross gradient term is used to quantify the structural similarity between the resistivity model and the prior model.

[0075] Step 3) Determine whether the nonlinear conjugate gradient iteration meets the exit condition. If yes, jump to step 9; otherwise, proceed to step 4.

[0076] Step 4) Determine the gradient of the objective function in two adjacent nonlinear conjugate gradient iterations. If the value is less than a preset threshold, proceed to step 5; otherwise, return to step 2.

[0077] Step 5) Adjust the cooling factor attenuation regularization parameters according to the preset parameters. and Reduce the weight constraints of the model smoothing term and cross gradient term in the objective function, so that the objective function can escape the local optimum, and then proceed to step 6).

[0078] Step 6) Determine the current particle swarm optimization trigger counter Is it greater than the maximum number of times? If yes, return to step 2) and continue with the nonlinear conjugate gradient iteration; otherwise, proceed to step 7).

[0079] Step 7) For the current model Preprocessing and distance metric construction are performed, followed by spatiotemporal density clustering analysis to extract scaling factor sequences. ;

[0080] Perform particle swarm optimization in the reduced parameter space, from the scaling factor sequence. Obtain the optimal scaling factor Utilizing the optimal scaling factor Generate updated resistivity model ;

[0081] Step 8) Update the counter And the number of iterations Determine whether the RMS after particle swarm optimization meets the exit condition. If it does, proceed to step 9; otherwise, return to step 2.

[0082] Step 9) Output the final resistivity model of the probe area and terminate the inversion process.

[0083] Example 2:

[0084] A three-dimensional hybrid optimization inversion method for magnetotelluric data, with the same technical content as in Example 1, further includes a cross-gradient term used to quantify the structural similarity between different physical property models, namely:

[0085] (1)

[0086] In the formula, For cross gradient terms; , This represents the gradient between the resistivity model and the prior model.

[0087] Example 3:

[0088] A three-dimensional hybrid optimization inversion method for magnetotelluric data, with the same technical content as any one of Examples 1-2, further comprising, in step 2), the objective function of the nonlinear conjugate gradient iteration. As shown below:

[0089] (2)

[0090] In the formula, This represents the affine transformation form of the model parameter m. It is an initial uniform half-space model. This is the model covariance matrix; d represents the magnetotelluric observation data. It is a forward function; These are cross gradient components; The data covariance matrix; , and These correspond to the observed data, model parameters, and the total number of cross gradient components, respectively.

[0091] Example 4:

[0092] A three-dimensional hybrid optimization inversion method for magnetotelluric data, with the same technical content as any one of Examples 1-3, further comprising the gradient of the objective function. As shown below:

[0093] (3)

[0094] (4)

[0095] In the formula, It is a forward response For model parameters Jacobian matrix; This is the partial derivative matrix of the cross gradient operator with respect to the model parameters; This refers to the data residuals.

[0096] Example 5:

[0097] A three-dimensional hybrid optimization inversion method for magnetotelluric data, with the same technical content as any one of embodiments 1-4, further comprising, in step 2), the updated resistivity model As shown below:

[0098] (5)

[0099] In the formula, Step size; Indicating the search direction; This serves as the initial resistivity model and a reference for transforming model parameters. and These are the resistivity model parameters in the transformed domain during the (k-1)th and kth iterations, respectively.

[0100] Example 6:

[0101] A three-dimensional hybrid optimization inversion method for magnetotelluric data, with the same technical content as any one of embodiments 1-5, further comprising, in step 7), optimizing the current resistivity model. Preprocessing and distance metric construction are performed, followed by spatiotemporal density clustering analysis to extract scaling factors. The steps include:

[0102] Step 7.1) Based on background resistivity and preset resistivity threshold Resistivity model Divided into low-resistivity subsets ( ) and high-impedance subsets ( );

[0103] Step 7.2) Use the grid topology index (i, j, k) instead of the geometric center coordinates to calculate the distance;

[0104] The topological spatial distance between two grid cells p and q As shown below:

[0105] (6)

[0106] Attribute distance between two grid cells p and q As shown below:

[0107] (7)

[0108] Step 7.3) Based on the topological spatial distance dist1 and attribute distance dist2, perform clustering independently on the low-resistivity subset and the high-resistivity subset respectively, and divide the preprocessed model unit into n clusters. Each cluster corresponds to a set of resistivity values. and a scaling factor Thus constructing the scaling factor vector This is used for subsequent particle swarm optimization search.

[0109] Example 7:

[0110] A three-dimensional hybrid optimization inversion method for magnetotelluric data, with the same technical content as any one of Examples 1-6, further comprising the following steps: when performing particle swarm optimization, the objective function... As shown below:

[0111] (8)

[0112] In the formula, This is the resistivity model obtained from the current NLCG iteration; The n-dimensional scaling factor vector generated by spatiotemporal density clustering; This represents the scaling factor vector x on the model. The forward response of the model obtained by scaling the resistivity of each cluster.

[0113] Example 8:

[0114] A three-dimensional hybrid optimization inversion method for magnetotelluric analysis, with the same technical content as any one of Examples 1-7, further wherein when performing particle swarm optimization, each particle in the particle swarm represents a set of scaling factor vectors.

[0115] The particle state is updated as follows:

[0116] (9)

[0117] (10)

[0118] In the formula, , Let be the particle velocities at time t+1 and time t; , The positions of the particles at time t+1 and time t; , Individual learning factors and social learning factors; , The result is a random number (between 0 and 1). , For individual optimal position and global optimal position;

[0119] Updated resistivity model As shown below:

[0120] (11)

[0121] In the formula, The optimal scaling factor; It is a set of resistivity values.

[0122] Example 9:

[0123] A three-dimensional hybrid optimization inversion method for magnetotelluric data, with the same technical content as any one of Examples 1-8, further wherein the exit condition for nonlinear conjugate gradient iteration is that the root mean square error (RMS) reaches a preset target threshold, the cooling factor of the regularization parameter decays to a preset lower limit, or the total number of inversion iterations reaches a set upper limit.

[0124] Example 10:

[0125] A three-dimensional hybrid optimization inversion method for magnetotelluric data, with the same technical content as any one of Examples 1-9, further wherein the regularization parameter is updated through a cooling factor, i.e.:

[0126] (12)

[0127] (13)

[0128] In the formula, the cooling factor , ; , These are the updated regularization parameters.

[0129] Example 11:

[0130] A three-dimensional hybrid optimization inversion method for magnetotelluric data, with the same technical content as any one of Examples 1-10, further comprising the following core steps and principles in step 7):

[0131] Structural partitioning: First, using ST-DBSCAN clustering, the preprocessed model units are divided into n clusters with spatial and electrical similarity. .

[0132] Parameterization: Then, for each cluster Define a uniform "scaling factor". This scaling factor It is an independent optimization variable, and its physical meaning is that it represents the optimization performance for this cluster. The factor by which the resistivity values ​​of all grid cells within the grid are linearly scaled as a whole.

[0133] Optimization and Mapping: In subsequent particle swarm optimization, the variables to be optimized are derived from the high-dimensional original resistivity model. This is transformed into a low-dimensional scaling factor vector. Obtain the optimal scaling factor. Then, the entire resistivity model is updated using the mapping relationship defined by equation (11) in claim 8.

[0134] Example 12:

[0135] To verify the effectiveness of the method, this invention designed complex theoretical models such as rectangular double anomalies and stepped double anomalies for testing. Experimental results show that the hybrid strategy significantly improves the ability to distinguish and reconstruct subsurface anomalies: cross-gradient constraints make the anomaly boundaries clearer and sharper; the nonlinear global correction in the PSO stage significantly improves the accuracy of resistivity amplitude recovery and effectively weakens the spurious "tailing" phenomenon that easily occurs below low-resistivity anomalies. Compared with traditional inversion methods, this method ultimately achieves lower root mean square error and cross-gradient values, realizing a dual optimization of data fitting accuracy and model structure consistency, and promoting an important shift in magnetotelluric inversion interpretation from "scalar fitting" to "vectorized structural constraints".

[0136] The specific steps of the magnetotelluric three-dimensional hybrid optimization inversion method are as follows:

[0137] Step 1) Initialization: Set the initial resistivity model of the detection area. Regularization parameters and Cooling factor and Maximum number of triggers for Particle Swarm Optimization (PSO) and convergence criteria;

[0138] Step 2) Local Inversion: Perform nonlinear conjugate gradient (NLCG) iteration, introduce cross gradient constraints, and update the current resistivity model. .

[0139] Step 3) Convergence Determination and Cooling: Determine if the NLCG inversion meets the exit conditions (including root mean square error RMS being less than a preset threshold, cooling factor being less than a preset threshold, or reaching the maximum number of iterations); if yes, proceed to Step 7 to output the final model; if no, determine if the objective function descent is less than a certain threshold (e.g., ...). If this condition holds true, it means the inverted solution has fallen into a local optimum. Therefore, reducing the regularization factor lowers the constraints on the model's smoothing and cross-gradient terms, allowing the objective function to escape the local optimum. The regularization parameter is updated by the cooling factor. as well as ,in , .

[0140] Step 4) Global Optimization Trigger: Determine the current PSO trigger counter. Has the maximum number of times been exceeded? ;like > If the global optimization step is successful, return to step 2 and continue with the NLCG local inversion; otherwise, proceed to step 5.

[0141] Step 5) Dimensionality reduction and global optimization: For the current model After dimensionality reduction, spatiotemporal density clustering (ST-DBSCAN) analysis was performed to extract the scaling factor. Perform PSO optimization in the reduced parameter space to obtain the optimal scaling factor. ;use Generate updated resistivity model .

[0142] Step 6) Loop control: Update the counter And the number of iterations Determine whether the RMS optimized by PSO meets the exit conditions. If it does, proceed to step 7; otherwise, return to step 2.

[0143] Step 7) Output: Output the final resistivity model of the probe area and terminate the inversion process.

[0144] The cross gradient component The formula for quantifying the structural similarity between different physical property models is as follows:

[0145] (1)

[0146] It ensures that the resistivity model and the prior model are consistent on the constructed boundary by forcing the spatial gradient fields to be collinear.

[0147] The objective function in the local optimization stage of step 2) is defined as:

[0148] (2)

[0149] In the formula, The model parameter m was subjected to an affine transformation. It is an initial uniform half-space model. This is the model covariance matrix; d represents the magnetotelluric observation data. It is a forward function; These are cross gradient components; The data covariance matrix; , and These correspond to the observed data, model parameters, and the total number of cross-gradient components, respectively.

[0150] The simplified formula for calculating the partial derivatives of the cross gradient term with respect to the model parameters is as follows:

[0151] (3)

[0152] In the formula, Let be the partial derivative matrix of the cross gradient operator with respect to the model parameters; thus, the global gradient expression of the objective function is... Simplified to:

[0153] (4)

[0154] In the formula, It is a forward response For model parameters The Jacobian matrix.

[0155] The update of the model parameters in the NLCG stage follows the formula below.

[0156] (5)

[0157] In the formula, step size Determined through linear search; The initial resistivity model (uniform half-space model) is defined in step 1) of claim 1 and is used only as a reference for the transformation of model parameters; and These are the resistivity model parameters in the transformed domain during the (k-1)th and kth iterations, respectively.

[0158] The specific methods for determining the spatial dimensionality reduction and scaling factor of ST-DBSCAN include:

[0159] 5.1 Electrical binary classification: Based on the background resistivity, the anomalous body is divided into a low-resistivity subset and a high-resistivity subset, and clustering is performed independently for each subset to enhance the homogeneity within the cluster;

[0160] 5.2 Topological Coordinate Representation: The distance is calculated using the grid topological index (i, j, k) instead of the geometric center coordinates, eliminating spatial deviations caused by non-uniform grid size variations. The topological spatial distance between two grid cells p and q is defined as:

[0161] (6)

[0162] Simultaneously, the difference in model parameters (logarithmic resistivity) is introduced as the attribute distance:

[0163] (7)

[0164] The ST-DBSCAN algorithm will take into account dist1 and dist2 to identify dense regions that are similar in both space and electrical properties.

[0165] 5.3 Parameter Correlation: Divide the dimensionality-reduced model parameters into n clusters. Each cluster corresponds to a set of resistivity values. and a scaling factor .

[0166] The objective function (fitness function) in the global optimization phase of PSO Defined as:

[0167] (8)

[0168] In the formula, This is the resistivity model obtained from the current NLCG iteration; n scaling factor vectors generated by ST-DBSCAN clustering.

[0169] The PSO-based model update formula is as follows:

[0170] 6.1) Particle State Update: Each particle in the particle swarm represents a scaling factor vector x, and its velocity v and position are updated according to the following formula. :

[0171] (9)

[0172] (10)

[0173] 6.2) Model parameter reconstruction: using the optimal scaling factor Update the grid resistivity within the i-th cluster:

[0174] (11)

[0175] This exponential mapping relationship enables nonlinear correction of model parameters after global search.

[0176] The hybrid optimization strategy and its exit conditions are as follows:

[0177] Collaborative Mechanism: This strategy achieves initial efficient convergence and structural constraints through the NLCG local optimization phase; when the NLCG iteration meets the preset triggering conditions and the PSO trigger count is reached... At that time, the global optimization phase is triggered;

[0178] Global optimization exit: In the PSO stage, ST-DBSCAN clustering is used to reduce the dimensionality of the high-dimensional model space to the scaling factor space; when the RMS of the model optimized by PSO meets the preset global exit condition (such as the RMS reaching the target threshold or the iteration reaching the upper limit), the global loop is exited.

[0179] Final termination condition: If the model after global optimization still fails to meet the preset convergence criterion, the updated model parameters will be used as the initial values ​​to return to the NLCG local optimization stage for fine iteration; until the system meets any of the following termination conditions: the root mean square error (RMS) reaches the preset target threshold, the cooling factor of the regularization parameter decays to the preset lower limit, or the total number of inversion iterations reaches the set upper limit.

[0180] Example 13:

[0181] Figure 1 This invention is illustrated using a resistivity model comprising high and low resistivity anomalies and a corresponding a priori velocity model as an example for verification. The resistivity model is as follows: An embedded in a uniform half-space of 100 Ω·m A low-resistivity anomalous body of 10 Ω·m and a Two high-resistivity anomalies, each with a diameter of 1000 Ω·m, have dimensions of 15 km × 15 km × 7 km in the x, y, and z directions. The corresponding a priori velocity model is as follows: 5km / s 1km / s and 9km / s.

[0182] The forward modeling used 104 measuring points, evenly distributed across a 35km × 60km area, with both line and point spacing of 5km. The study area was discretized into a 38×38×21 grid, and ModEM was used to simulate the electromagnetic response at 17 frequencies (10⁻²~10²Hz) on the Earth's surface. To simulate actual observation errors, 2% Gaussian random noise was added to the synthesized impedance data. The inversion modeling used the same grid partitioning as the forward modeling, with the initial model set as a uniform half-space of 100Ω·m.

[0183] The key inversion parameters are as follows: regularization factor , PSO is triggered once (60 particles, maximum 100 iterations). Clustering parameters are based on log resistivity. The statistical characteristics are set for the model parameters: In the low-resistivity region, Eps1 = 4.5, Eps2 = 0.103, and MinPts = 8; when In the high-resistivity region, we take Eps1 = 4.5, Eps2 = 0.05, and MinPts = 8. In logarithmic space... Cluster analysis was performed. The low-resistivity anomaly data were relatively concentrated, so a larger Eps2 (0.103) was used to ensure their spatial continuity; the high-resistivity anomaly data were relatively dispersed, so a smaller Eps2 (0.05) was used to ensure uniformity within the cluster. The spatial parameter Eps1 (4.5) was set to match the grid scale.

[0184] Figure 2 This paper presents the preliminary resistivity model slices obtained by the hybrid inversion before PSO optimization and their cross-gradient distribution with the velocity prior model, and compares them with the results of the conventional inversion. (a1)-(a4) correspond to the conventional inversion results, and (b1)-(b4) correspond to the preliminary hybrid inversion results with cross-gradient constraints. The slice locations include horizontal slices (z=7km) and vertical slices (x=-12.5 km). From the resistivity model slices, the boundaries of high-resistivity and low-resistivity anomalies in the hybrid inversion results (b1 and b3) are significantly clearer and sharper than those in the conventional inversion (a1 and a3), and the corresponding cross-gradient distribution further quantifies the structural consistency. The cross-gradient values ​​in (b2 and b4) are generally significantly lower than those in (a2 and a4), especially in the anomaly boundary regions, where the cross-gradients in the hybrid inversion are very small (light gray), indicating a high degree of structural coupling between the resistivity model and the velocity prior model at these locations; while in the conventional inversion (a2 and a4), the cross-gradient values ​​at the boundaries are higher (dark gray), reflecting poor consistency between the electrical structure and the velocity structure. The above comparison clearly shows that the cross gradient constraint term effectively plays a structural guiding role in the early stage of hybrid inversion, prompting the resistivity model to converge in a direction consistent with the structure of the prior velocity model, thereby improving the recovery accuracy of the anomaly boundary and the overall geological rationality of the model, and providing a better starting model for subsequent PSO global optimization.

[0185] Figure 3The results of clustering resistivity model parameters before PSO are presented. Subplot (a) shows a 3D view after clustering, clearly displaying the spatial distribution of the seven different clusters: dark blue and light blue clusters mainly correspond to the background area, while other colored clusters are located inside low- and high-resistivity anomalies. The overall clustering result is highly consistent with the boundary of the inverted model. Subplots (b) and (c) are slices of the xz and yz planes, respectively, further demonstrating the recognition effect of clustering in the vertical and horizontal directions. It can be seen that at the boundary of the anomaly, clustering can clearly distinguish the different electrical parameters on both sides of the boundary, forming sharp cluster boundaries, avoiding parameter mixing, and thus achieving targeted optimization. Data analysis shows that a total of seven clusters were generated, corresponding to seven scaling factor variables, reducing the original mesh size of 38×38×21 (approximately 30,324 model parameters) to only seven, with parameter dimensionality compression exceeding 99.98%. This significantly reduces the search space and computational cost of PSO, improves optimization efficiency and convergence speed, while preserving key structural information of the model, ensuring that subsequent hybrid inversion can effectively recover the details of complex geological boundaries and high and low resistivity anomalies.

[0186] Figure 4 The cross-sectional results in the x, y, and z directions were compared among conventional inversion, constrained inversion, and hybrid inversion. While conventional inversion exhibits good lateral resolution, its longitudinal resolution is significantly insufficient. As a smooth inversion method, it ensures solution stability but also leads to a gradual transition in resistivity at the anomaly boundary, resulting in blurred boundary morphology, especially prone to spurious "tailing" phenomena below low-resistivity anomalies. Constrained inversion achieves high resolution in both the lateral and longitudinal directions, but still struggles to effectively suppress spurious anomalies near the anomaly boundary and suffers from insufficient resistivity recovery. In contrast, hybrid inversion, by fully utilizing the clear structure obtained from cross-gradient constraints, effectively suppresses spurious anomalies near the boundary and significantly improves the accuracy of anomaly resistivity recovery, thereby achieving an overall improvement in the performance of structural and physical property parameter inversion and reducing the ambiguity of the inversion results.

[0187] Figure 5 The resistivity statistical distributions of the three inversion methods are presented, revealing their differences in anomaly recovery and boundary characterization. Figure 5 The conventional inversion of 'a' as a smoothing method not only leads to insufficient recovery of resistivity values ​​for low-resistivity and high-resistivity anomalies (especially high-resistivity anomalies are affected by the insensitivity to airborne magnetotelluric data), but also causes blurring of anomaly boundaries due to the continuous and gradual change in its solution. In contrast, the sub- Figure 5 The constrained inversion of b enhances structural information by introducing cross-gradient constraints, making the resistivity values ​​closer to the true values ​​of 10 Ω·m and 10 3 Ω·m, the boundary is clearer, but incomplete recovery still exists. Figure 5The hybrid inversion of c further optimized the results, with the discontinuity of resistivity distribution at the boundary being the most significant, and the boundary being greatly improved. At the same time, the resistivity of both low and high resistivity anomalies was restored to near the true value, and their statistical distribution was closest to the true model, indicating that the hybrid inversion significantly improved the accuracy of resistivity parameter recovery.

[0188] Figure 6 Three inversion methods, RMS and The curves showing the change with the number of iterations. The root mean square error of the conventional inversion, constrained inversion, and hybrid inversion converged from 4.27 to 1.060, 1.063, and 1.050, respectively. From the convergence process... Figure 6 (a) While conventional inversion converges rapidly, it is prone to getting trapped in local optima, leading to a high final RMS. Constrained inversion, although introducing structural constraints, still struggles to completely suppress spurious boundary anomalies and insufficient resistivity recovery, and is similarly prone to getting trapped in local optima. In contrast, hybrid inversion escapes local optima through PSO global search, demonstrating stronger global optimization capabilities. Furthermore, hybrid inversion consistently maintains lower cross gradient values, confirming the effectiveness of cross gradient constraints in maintaining model structural consistency. Figure 6 b). Figure 6 c further illustrates the decrease in RMS during the PSO optimization phase. With a significantly reduced number of parameters, PSO exhibits efficient search performance, and the RMS stabilizes after 15 iterations. Considering both convergence behavior and computational efficiency, setting the maximum number of PSO iterations to 100 effectively conserves computational resources while ensuring solution quality.

[0189] Example 14:

[0190] like Figure 7 The high- and low-resistivity dual-anomaly model with a stepped boundary, as shown in the embodiment of the present invention, is used for verification. This model embeds two anomaly regions with different resistivity values ​​within a uniform half-space background of 100 Ω·m: a high-resistivity anomaly of 1000 Ω·m and a low-resistivity anomaly of 10 Ω·m, to simulate more realistic irregular geological boundaries and stepped resistivity transition characteristics. The corresponding velocity prior model is set to s0 = 1.65 km / s (background), s1 = 2.12 km / s, and s2 = 1.28 km / s, as shown below. Figure 9 As shown. This velocity prior model aims to provide structured auxiliary constraints to help the inversion process better capture the spatial distribution and boundary features of anomalies. The forward simulation uses 21 frequencies (from 10...). -3 Hz to 10 3The electromagnetic response data (Hz) were obtained, and other forward modeling parameters (such as measurement point layout and grid partitioning) remained consistent with the aforementioned model to ensure experimental comparability. 2% Gaussian random noise was added to the synthetic impedance data to simulate the error effects of actual observations. The inversion process used the same grid partitioning as the forward model, and the initial model remained a uniform half-space of 100 Ω·m.

[0191] The key inversion parameters are set as follows: the regularization factor remains unchanged, and PSO is triggered twice (50 particles, maximum number of iterations 100). The clustering parameters are optimized and adjusted in stages based on the statistical distribution characteristics of log resistivity. In the first clustering, the parameters are set to Eps1=4.5, Eps2=0.08, MinPts=8 and Eps1=4.5, Eps2=0.08, MinPts=8 for different outlier regions, respectively. In the second clustering, the parameters are further refined and set to Eps1=4.5, Eps2=0.05, MinPts=8 and Eps1=4.5, Eps2=0.04, MinPts=8. Cluster analysis was performed in logarithmic space. Given the complex distribution of anomalous data in this stepped model, the first clustering used a relatively large Eps² (0.08) to enhance cluster inclusiveness and spatial continuity, avoiding premature fragmentation. The second clustering gradually reduced Eps² to 0.05 and 0.04 to improve the uniformity of resistivity within clusters and the sharpness of boundaries, ensuring accurate characterization of high-resistivity and low-resistivity anomalous bodies. The spatial neighborhood parameter Eps¹ (4.5) was consistently matched to the grid scale to maintain geometric consistency of the clustering. This phased clustering strategy helps the algorithm gradually converge to sharper stepped boundaries while suppressing noise interference, improving the accuracy and geological plausibility of the model reconstruction.

[0192] Figure 8Three inversion results for complex boundary anomaly models were compared: conventional inversion (Fig. a), constrained inversion (Fig. b), and the hybrid inversion proposed in this paper (Fig. c). These results are presented in slice views at x = -12.5 km, y = -3.75 km, and z = 4 km, where the color bars represent logarithmic resistivity logρ, ranging from 1 to 3, corresponding to resistivity values ​​from 10 Ω·m to 1000 Ω·m. Analysis shows that, compared with conventional and constrained inversion, the proposed hybrid inversion method can still effectively improve the accuracy of anomaly boundary recovery in complex structural backgrounds. For example, in x and y slices, the hybrid inversion clearly depicts the sharp separation between the stepped boundary and the high and low resistivity anomalies, while the conventional inversion shows blurred boundaries and color gradient transitions, and although the constrained inversion shows some improvement, edge diffusion still exists. At the same time, this method significantly improves the restoration effect of resistivity parameters. In z slice, the anomaly shape of the hybrid inversion is closer to the rectangular or stepped features of the real model, avoiding the excessive smoothing of resistivity values ​​common in conventional methods. In addition, it greatly suppresses the false "tailing" phenomenon commonly seen below low resistivity anomalies. For example, obvious downward extension artifacts (the red low resistivity area extends to the background area) can be seen in slices (a1) and (a2). Although the tailing is weakened in constrained inversion (row b), it still remains, while the hybrid inversion (a1) shows a different effect. Figure 8 c1) was almost completely eliminated, ensuring the structural integrity and geological authenticity of the model. This was achieved from the statistical distribution of the resistivity parameter ( Figure 9 As can be seen, the resistivity distribution of the hybrid inversion result (c) is closest to the true model (orange dashed line). Its main anomalies are concentrated in the 10Ω·m, 100Ω·m, and 1000Ω·m ranges. The peak positions correspond more accurately, the distribution width is narrower, and there is no significant shift. This is significantly better than the relatively discrete distribution, peak shift towards lower values, and long-tail bias in the conventional inversion (a), and the results of the constrained inversion (b), which, although the peaks are close, still have a wide and slightly biased distribution. This statistical analysis further quantifies the superiority of the hybrid inversion. The concentration of its parameter distribution and the peak matching degree indicate that the algorithm is more robust to noise interference, thereby improving the overall reliability and accuracy of the inversion.

[0193] Figure 10 The paper demonstrates the RMS and its variation with iteration number for three inversion methods under a complex boundary anomaly model. Figure 10 As shown in Figure a, the RMS values ​​of conventional inversion, constrained inversion, and hybrid inversion converge from the initial value of 9.19 to 1.077, 1.098, and 1.062, respectively, indicating that hybrid inversion can more effectively avoid getting trapped in local extrema, thereby obtaining a better global solution. Figure 10b shows that both constrained inversion and hybrid inversion maintained low cross-gradient values ​​throughout the inversion process, reflecting the effectiveness of both in constraining structural consistency. For this complex model, to balance inversion accuracy and computational efficiency, two PSO optimizations were performed in the hybrid inversion. The decreasing trend of RMS during the optimization process is shown in Figure [Figure number missing]. Figure 10 As shown in c.

Claims

1. A three-dimensional hybrid optimization inversion method for magnetotelluric data, characterized in that, Includes the following steps: Step 1) Set the initial resistivity model of the detection area Regularization parameters and Cooling factor and Maximum number of triggers for particle swarm optimization and convergence criteria; Step 2) Using the objective function, which includes data fitting terms, model smoothing terms, and cross-gradient terms, as constraints, perform nonlinear conjugate gradient iteration to update the resistivity model, denoted as... ; Among them, the model smoothing term is used to constrain the spatial smoothness of the resistivity model, and the cross gradient term is used to quantify the structural similarity between the resistivity model and the prior model. Step 3) Determine whether the nonlinear conjugate gradient iteration meets the exit condition. If yes, jump to step 9; otherwise, proceed to step 4. Step 4) Determine the gradient of the objective function in two adjacent nonlinear conjugate gradient iterations. If the value is less than a preset threshold, proceed to step 5; otherwise, return to step 2. Step 5) Adjust the cooling factor attenuation regularization parameters according to the preset parameters. and Reduce the weight constraints of the model smoothing term and cross gradient term in the objective function, so that the objective function can escape the local optimum, and then proceed to step 6). Step 6) Determine the current particle swarm optimization trigger counter Is it greater than the maximum number of times? If yes, return to step 2) and continue with the nonlinear conjugate gradient iteration; otherwise, proceed to step 7). Step 7) For the current model Preprocessing and distance metric construction are performed, followed by spatiotemporal density clustering analysis to extract scaling factor sequences. ; Perform particle swarm optimization in the reduced parameter space, from the scaling factor sequence. Obtain the optimal scaling factor Utilizing the optimal scaling factor Generate updated resistivity model ; Step 8) Update the counter And the number of iterations Determine whether the RMS after particle swarm optimization meets the exit condition. If it does, proceed to step 9; otherwise, return to step 2. Step 9) Output the final resistivity model of the probe area and terminate the inversion process.

2. The magnetotelluric three-dimensional hybrid optimization inversion method according to claim 1, characterized in that, The cross gradient term is used to quantify the structural similarity between different physical property models, i.e.: (1) In the formula, For cross gradient terms; , This represents the gradient between the resistivity model and the prior model.

3. The magnetotelluric three-dimensional hybrid optimization inversion method according to claim 1, characterized in that, In step 2), the objective function of the nonlinear conjugate gradient iteration is... As shown below: (2) In the formula, This represents the affine transformation form of the model parameter m. It is an initial uniform half-space model. This is the model covariance matrix; d represents the magnetotelluric observation data. It is a forward function; These are cross gradient components; The data covariance matrix; , and These correspond to the observed data, model parameters, and the total number of cross gradient components, respectively.

4. The magnetotelluric three-dimensional hybrid optimization inversion method according to claim 3, characterized in that, gradient of the objective function As shown below: (3) (4) In the formula, It is a forward response For model parameters The Jacobian matrix; This is the partial derivative matrix of the cross gradient operator with respect to the model parameters; This refers to the data residuals.

5. The magnetotelluric three-dimensional hybrid optimization inversion method according to claim 1, characterized in that, In step 2), the updated resistivity model As shown below: (5) In the formula, Step size; Indicating the search direction; This serves as the initial resistivity model and a reference for transforming model parameters. and These are the resistivity model parameters in the transform domain during the (k-1)th and kth iterations, respectively.

6. The magnetotelluric three-dimensional hybrid optimization inversion method according to claim 1, characterized in that, In step 7), the current resistivity model is... Preprocessing and distance metric construction are performed, followed by spatiotemporal density clustering analysis to extract scaling factors. The steps include: Step 7.1) Based on background resistivity and preset resistivity threshold Resistivity model It is divided into low-resistivity subsets and high-resistivity subsets; Step 7.2) Use the grid topology index (i, j, k) instead of the geometric center coordinates to calculate the distance; The topological spatial distance between two grid cells p and q As shown below: (6) Attribute distance between two grid cells p and q As shown below: (7) Step 7.3) Based on the topological spatial distance dist1 and attribute distance dist2, perform clustering independently on the low-resistivity subset and the high-resistivity subset respectively, and divide the preprocessed model unit into n clusters. Each cluster corresponds to a set of resistivity values. and a scaling factor Thus constructing the scaling factor vector This is used for subsequent particle swarm optimization search.

7. The magnetotelluric three-dimensional hybrid optimization inversion method according to claim 1, characterized in that, When performing particle swarm optimization, the objective function As shown below: (8) In the formula, This is the resistivity model obtained from the current NLCG iteration; The n-dimensional scaling factor vector generated by spatiotemporal density clustering; This represents the scaling factor vector x on the model. The forward response of the model obtained by scaling the resistivity of each cluster.

8. The magnetotelluric three-dimensional hybrid optimization inversion method according to claim 1, characterized in that, When performing particle swarm optimization, each particle in the particle swarm represents a set of scaling factor vectors; The particle state is updated as follows: (9) (10) In the formula, , Let be the particle velocities at time t+1 and time t; , The positions of the particles at time t+1 and time t; , Individual learning factors and social learning factors; , It is a random number; , For individual optimal position and global optimal position; Updated resistivity model As shown below: (11) In the formula, The optimal scaling factor; It is a set of resistivity values.

9. The magnetotelluric three-dimensional hybrid optimization inversion method according to claim 1, characterized in that, The exit conditions for nonlinear conjugate gradient iteration are: the root mean square error (RMS) reaches a preset target threshold, the cooling factor of the regularization parameter decays to a preset lower limit, or the total number of inversion iterations reaches a set upper limit.

10. The magnetotelluric three-dimensional hybrid optimization inversion method according to claim 1, characterized in that, The regularization parameter is updated via the cooling factor, i.e.: (12) (13) In the formula, the cooling factor , ; , These are the updated regularization parameters.

Images