Method and device for reconstructing crust temperature distribution based on magnetic anomaly

Abstract

This application belongs to the field of geophysical magnetic exploration technology, specifically disclosing a method and apparatus for reconstructing crustal temperature distribution based on magnetic anomalies. The method includes: obtaining magnetic susceptibility-temperature curves based on rock physics experiments; constructing a depth-dependent magnetic susceptibility function by combining a preset temperature-depth relationship; using the depth-dependent magnetic susceptibility function as the magnetic susceptibility structure constraint for magnetic anomaly forward and inversion models; inputting magnetic anomaly grid data and the depth-dependent magnetic susceptibility function into a spatial domain forward and inversion framework established based on Cauchy surface integrals; constructing an objective functional; iteratively solving the functional using a reweighted regularized conjugate gradient algorithm; and inverting to obtain the equivalent Curie depth distribution; using the equivalent Curie depth distribution as the first isotherm; combining the surface boundary temperature and the temperature boundary with the Moho discontinuity as the second isotherm; and reconstructing the three-dimensional crustal temperature structure of the region and outputting the temperature distribution results. This application can improve the efficiency of Curie inversion calculations and the geological reliability of the results.

Description

Technical Field

[0001] This application belongs to the field of geophysical magnetic exploration technology, and more specifically, relates to a method and apparatus for reconstructing crustal temperature distribution based on magnetic anomalies. Background Technology

[0002] Crustal temperature structure is a crucial fundamental parameter for studying lithospheric rheology, tectonic evolution, and geothermal resource assessment. Traditional temperature acquisition methods rely on borehole thermometry and surface heat flow observations, which are costly and have sparse spatial coverage, making them unsuitable for continuous modeling at regional scales. Magnetic anomaly data, with its wide coverage and sensitivity to subsurface magnetic structures, is often used to constrain the magnetic crustal floor. Magnetic minerals demagnetize near the Curie temperature, and their corresponding interfaces are often approximated as isothermal surfaces around 580°C. Obtaining the Curie depth through magnetic anomaly inversion can provide indirect constraints on deep temperature structure.

[0003] Among existing Curie estimation methods, spectral analysis is highly sensitive to window division and parameter selection, limiting its applicability; interface inversion methods based on Parker-Oldenburg expansion struggle to stably distinguish the contributions of interface undulations and spatial variations in magnetic susceptibility to anomalies; spatial domain three-dimensional volume inversion or equivalent source methods suffer from large computational scales and ill-conditioned nature, hindering efficient and stable mapping at the regional scale. Furthermore, magnetic susceptibility exhibits significant nonlinearity with temperature. Adopting a uniform magnetic susceptibility assumption easily leads to trade-off errors between interface geometry and magnetic susceptibility structure. Rock physics experiments can provide magnetic susceptibility-temperature curves. By combining the temperature-depth relationship, a structure can be constructed. If explicitly introduced in the inversion... Constraints can improve the physical consistency and stability of results.

[0004] Furthermore, existing technologies mostly remain at the Curie estimation stage, lacking an integrated process that incorporates the Curie line with thermal boundaries such as the Earth's surface and the Moho discontinuity into the temperature modeling framework. Therefore, this application urgently needs to address the problems of low computational efficiency in Curie inversion calculations, unstable trade-offs between interface geometry and magnetic susceptibility structure, and the lack of an integrated reconstruction process from magnetic anomalies to three-dimensional crustal temperature structure in existing technologies. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this application aims to provide a method and apparatus for reconstructing crustal temperature distribution based on magnetic anomalies. This method and apparatus are intended to solve the problems of low Curie inversion calculation efficiency caused by the huge computational load, difficulty in discretization, and long iteration time and slow speed of conventional inversion algorithms due to the calculation of magnetic anomalies.

[0006] To achieve the above objectives, in a first aspect, this application provides a method for reconstructing crustal temperature distribution based on magnetic anomalies, comprising: Acquire magnetic anomaly grid data for the target area; Based on the magnetic susceptibility-temperature curve obtained from rock physics experiments, a depth-dependent magnetic susceptibility function is constructed by combining the preset temperature-depth relationship. The depth-dependent magnetic susceptibility function is used as the magnetic susceptibility structure constraint for the forward and inverse modeling of magnetic anomalies. The magnetic anomaly grid data and the depth-related magnetic susceptibility function are input into the spatial domain forward and inverse modeling framework established based on Cauchy surface integral. An objective functional including data fitting terms and model stabilization terms is constructed. The reweighted regularized conjugate gradient algorithm is used to iteratively solve the problem and obtain the equivalent Curie depth distribution through inversion. Using the equivalent Curie depth distribution as the first isotherm, and combining it with the surface boundary temperature and the temperature boundary with the Moho surface as the second isotherm, the three-dimensional crustal temperature structure of the region is reconstructed and the temperature distribution results are output.

[0007] Secondly, this application also provides a device for reconstructing crustal temperature distribution based on magnetic anomalies, comprising: The data acquisition module is used to acquire magnetic anomaly grid data of the target area; The magnetic susceptibility constraint construction module is used to construct a depth-related magnetic susceptibility function based on the magnetic susceptibility-temperature curve obtained from rock physics experiments and a preset temperature-depth relationship, and to use the depth-related magnetic susceptibility function as the magnetic susceptibility structure constraint for the forward and inverse modeling of magnetic anomalies. The inversion module is used to input the magnetic anomaly grid data and the depth-related magnetic susceptibility function into the spatial domain forward and inversion framework based on Cauchy surface integral, construct an objective functional including data fitting terms and model stabilization terms, and iteratively solve it using the reweighted regularized conjugate gradient algorithm to obtain the equivalent Curie depth distribution. The temperature reconstruction module is used to reconstruct the three-dimensional crustal temperature structure of the region by taking the equivalent Curie depth distribution as the first isotherm, combining it with the surface boundary temperature and the temperature boundary with the Moho surface as the second isotherm, and outputting the temperature distribution results.

[0008] Thirdly, this application provides an electronic device, comprising: at least one memory for storing a program; and at least one processor for executing the program stored in the memory, wherein when the program stored in the memory is executed, the processor is configured to execute the method described in the first aspect or any possible implementation thereof.

[0009] Fourthly, this application provides a computer-readable storage medium storing a computer program that, when run on a processor, causes the processor to perform the method described in the first aspect or any possible implementation thereof.

[0010] Fifthly, this application provides a computer program product that, when run on a processor, causes the processor to execute the method described in the first aspect or any possible implementation thereof.

[0011] It is understood that the beneficial effects of the second to fifth aspects mentioned above can be found in the relevant descriptions in the first aspect mentioned above, and will not be repeated here.

[0012] Overall, the technical solutions conceived in this application have the following beneficial effects compared with the prior art: (1) This application adopts a spatial domain forward modeling operator based on Cauchy surface integral to replace the traditional three-dimensional volume integral form in the form of two-dimensional surface integral to calculate the magnetic anomaly response generated by the interface of sudden change in magnetic susceptibility, which greatly reduces the amount of forward modeling computation. By combining regular grid discretization and plane approximation strategy, the continuous integral is transformed into the summation operation of discrete unit contribution, which further improves the computational efficiency. At the same time, the reweighted regularized conjugate gradient algorithm is used for iterative solution, which significantly improves the efficiency of regional scale in-line mapping while ensuring the inversion accuracy.

[0013] (2) This application constructs a Tikhonov objective functional that includes a data fitting term and a model stabilization term, and balances the relationship between data fitting and model stability through regularization parameters; it adopts a reweighted regularized conjugate gradient algorithm, and ensures the convergence stability of the inversion process through iterative updates of the conjugate search direction and optimization of the line search step size; at the same time, the regularization parameters adopt an adaptive update strategy, and dynamically adjust the regularization strength according to the changes in model bias, further enhancing the robustness of the algorithm to noisy data and initial model bias.

[0014] (3) This application uses a model weighting matrix constructed based on the Fréchet derivative matrix to perform weighted parameterization of interface parameters, balancing the sensitivity differences of different parameters to magnetic anomaly response; at the same time, by introducing a magnetic susceptibility-temperature curve constrained by rock physics experiments that is more consistent with the real geological thermomagnetic characteristics, the trade-off error between interface geometry and magnetic susceptibility structure is significantly reduced, and the physical consistency and stability of the inversion results are improved.

[0015] (4) This application uses the equivalent Curie depth distribution obtained by inversion as the first isothermal surface and assigns it temperature. Combined with the surface boundary temperature and the temperature boundary with the Moho surface as the second isothermal surface, a complete technical chain is established from magnetic anomaly data input to Curie inversion and then to the output of three-dimensional crustal temperature structure. This realizes the integrated reconstruction of magnetic anomaly data to three-dimensional crustal temperature structure and fills the gap in the lack of a complete temperature modeling framework in the existing technology. Attached Figure Description

[0016] Figure 1This is a flowchart illustrating the method for reconstructing crustal temperature distribution based on magnetic anomalies provided in an embodiment of this application. Figure 2 The inner interface Γ and the asymptotic reference plane P and their reference depth in this application are defined as follows: A schematic diagram of the geometric relationship; Figure 3 (a) is a projection view of the equivalent curine geometry of the synthetic model of this application along the X direction; Figure 3 (b) is an equivalent curine depth plane distribution map of the synthetic model of this application; Figure 4 The rock physical magnetic susceptibility (normalized magnetic susceptibility)-temperature curve for this application. and the depth-dependent magnetic susceptibility function constructed from it. Schematic diagram; Figure 5 (a) is a schematic diagram of the synthetic magnetic anomaly obtained by forward modeling based on Cauchy surface integral in this application; Figure 5 (b) is a schematic diagram of the inversion input magnetic anomaly after adding noise to the synthesized magnetic anomaly; Figure 6 (a) is a projection view along the X direction of the equivalent Curie inversion result of this application; Figure 6 (b) is a planar distribution diagram showing the difference in depth between the equivalent Curie inversion result of this application and the real model; Figure 7 (a) is a planar distribution diagram of the equivalent Curie inversion results obtained using the thermomagnetic model assumption; Figure 7 (b) is a planar distribution diagram of the equivalent Curie inversion results obtained using the linear model assumption; Figure 7 (c) is a planar distribution diagram of the equivalent Curie inversion results obtained using the uniform model assumption; Figure 7 (d) is a planar distribution of the depth difference between the equivalent Curie inversion result and the real model under the assumption of the thermomagnetic model; Figure 7 (e) is a planar distribution of the depth difference between the equivalent Curie inversion result and the true model under the assumption of a linear model; Figure 7 (f) is a planar distribution of the depth difference between the equivalent Curie inversion result and the real model under the assumption of a uniform model; Figure 7 (g) represents the real model and different models along the representative profile (x=0). A depth comparison curve of the equivalent Curie inversion results under the assumed conditions; Figure 7 (h) for different Box plot of depth difference distribution in equivalent Curie inversion results under assumed conditions; Figure 8 (a) is a three-dimensional schematic diagram of the regional crustal temperature structure; Figure 8(b) is a temperature plane slice at a depth of 20 km; Figure 8 (c) is a temperature profile along the Y direction at X=0km, where the black curves represent the first isothermal boundary corresponding to the equivalent Curie and the second isothermal boundary corresponding to the Mohorovičić boundary, respectively. Figure 9 This is a schematic diagram of the device for reconstructing crustal temperature distribution based on magnetic anomalies provided in an embodiment of this application; Figure 10 This is a schematic diagram of the structure of the electronic device provided in the embodiments of this application. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0018] In this article, the term "and / or" describes the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, or B existing alone. The symbol " / " in this article indicates that the related objects are in an "or" relationship; for example, A / B means A or B.

[0019] The terms "first" and "second," etc., used in the specification and claims herein are used to distinguish different objects, not to describe a specific order of objects. For example, "first response message" and "second response message," etc., are used to distinguish different response messages, not to describe a specific order of response messages.

[0020] In the embodiments of this application, the terms "exemplary" or "for example" are used to indicate that something is an example, illustration, or description. Any embodiment or design that is described as "exemplary" or "for example" in the embodiments of this application should not be construed as being more preferred or advantageous than other embodiments or design. Specifically, the use of the terms "exemplary" or "for example" is intended to present the relevant concepts in a specific manner.

[0021] In the description of the embodiments of this application, unless otherwise stated, "multiple" means two or more, for example, multiple processing units means two or more processing units, multiple elements means two or more elements, etc.

[0022] The embodiments of this application are described below with reference to the accompanying drawings.

[0023] Reference Figure 1 This application provides a method for reconstructing crustal temperature distribution based on magnetic anomalies, including: S101. Obtain magnetic anomaly grid data for the target area; S102. Based on the magnetic susceptibility-temperature curve obtained from rock physics experiments, a depth-related magnetic susceptibility function is constructed by combining the preset temperature-depth relationship, and the depth-related magnetic susceptibility function is used as the magnetic susceptibility structure constraint for the forward and inverse modeling of magnetic anomalies. S103. Input the magnetic anomaly grid data and depth-related magnetic susceptibility function into the spatial domain forward and inverse modeling framework established based on Cauchy surface integral, construct an objective functional including data fitting terms and model stabilization terms, and use the reweighted regularized conjugate gradient algorithm to iteratively solve the problem, and inversely obtain the equivalent Curie depth distribution. S104. Using the equivalent Curie depth distribution as the first isothermal surface, and combining it with the surface boundary temperature and the temperature boundary with the Moho surface as the second isothermal surface, the three-dimensional crustal temperature structure of the region is reconstructed and the temperature distribution results are output.

[0024] Specifically, firstly, magnetic anomaly data of the target area is acquired via S101. The acquired magnetic anomaly data is then processed into a regular grid, organizing it into a grid data format with uniform spatial resolution. During or after gridding, necessary preprocessing is performed on the data. This preprocessing includes, but is not limited to, one or more operations among outlier removal, missing value imputation, mean removal, linear trend removal, boundary extension, upward extension, pole conversion, frequency domain filtering, and equivalent source reconstruction. After the above processing, regularized magnetic anomaly grid data of the target area is obtained, which serves as the input data for subsequent forward and inverse modeling calculations.

[0025] To verify the feasibility and stability of the method in this application, synthetic model data is used for illustration in the embodiments of this application. (Refer to...) Figure 2 and Figure 3 Construct a 3D synthetic model region, exemplarily setting the horizontal range and mesh interval, such as 400km × 400km, Δx = Δy = 5km. Set the reference depth. (e.g., 20km), and construct the Curie Γ near the reference depth, interface elevation function. It can be composed of several Gaussian undulations superimposed, resulting in spatial undulations in the equivalent Curie depth. An example is setting the observation surface near the Earth's surface (e.g., ), and set up regular observation grid points on the observation surface. .

[0026] It should be noted that in this application direction, direction and The directions can correspond to east-west, north-south, and depth, respectively. Let these be the coordinates of the interface points. These are the coordinates of the observation point.

[0027] Secondly, using S102, based on the magnetic susceptibility-temperature curve from rock physics experiments... Combining temperature-depth relationship Construct a deep correlation magnetic susceptibility function . Reference Figure 4 , These measurements can be derived from κ-T or χ-T measurements obtained from rock physics experiments. It can be a nonlinear curve, a piecewise linear curve, or a fitted function thereof. Temperature-depth relationship. A linear geothermal gradient model can be used. Alternatively, a piecewise gradient model can be used. By... Substitution get This allows the nonlinear characteristics of magnetic susceptibility changing with temperature to be explicitly introduced into subsequent forward and inverse calculations.

[0028] In an optional implementation, The function is depth-only across the entire region; in another alternative implementation, the temperature structure can be updated after the temperature structure is updated. Iterative corrections are made to form a consistent "temperature-magnetic susceptibility" constraint.

[0029] The depth-dependent magnetic susceptibility function describes the nonlinear distribution characteristics of magnetic susceptibility with depth. It is used as a magnetic susceptibility structure constraint in the forward and inverse calculations of magnetic anomalies. In the subsequent forward and inverse processes, it is used as a known fixed physical property parameter to describe the magnetic susceptibility difference on the upper and lower sides of the interface of abrupt change in underground magnetic susceptibility, thereby effectively separating the coupling influence of interface geometric parameters and magnetic susceptibility parameters on the magnetic anomaly response.

[0030] Furthermore, the magnetic anomaly and depth-related magnetic susceptibility model is input into a spatial domain forward and inverse modeling framework based on Cauchy surface integrals to obtain the equivalent Curie depth distribution. (Refer to...) Figure 2 This application uses the subsurface magnetic susceptibility abrupt change interface, i.e., Curie Γ, as the interface model to be inverted, and sets an asymptotic reference plane P or a reference depth. , such that Γ asymptotically approaches P at infinity.

[0031] In forward modeling, the magnetic field components generated by interface Γ are calculated using a Cauchy surface integral expression. Furthermore, a two-dimensional surface integral form replaces the three-dimensional volume integral form, thereby reducing the computational dimensionality and improving the efficiency and numerical stability of forward modeling. To enable numerical calculations, the horizontal projection of interface Γ can be divided into... grid cells A local plane approximation is used within each grid cell, and at the observation point The total field anomaly is obtained by summing the contributions of each unit. . Reference Figure 5 Based on the above forward modeling calculations, synthetic magnetic anomaly data can be obtained. For example, it is possible to The observation data are formed by superimposing ±5% (or ±10%) uniform random noise. and will As inversion input data .

[0032] In the inversion calculation, a Tikhonov objective functional is constructed, which includes data fitting and stabilization terms. For example, the prior model vector... Reference depth The corresponding planar or smooth interface. The analytical Fréchet derivative matrix is ​​obtained by taking the partial derivatives of the model parameters with respect to the Cauchy surface integral kernel function. The weighting matrix of the model is constructed from the comprehensive sensitivity. The interface parameters are then weighted and parameterized. Subsequently, a reweighted regularized conjugate gradient algorithm is used to iteratively solve for the minimum value of the objective functional, and the regularization parameters can be adaptively adjusted according to the model update. The iteration termination condition can be that the data fit reaches a preset tolerance, the model update amount is less than a threshold, or the maximum number of iterations is reached. The final output is the depth distribution of the interface Γ, which is then used as the equivalent Curie depth distribution. . Reference Figure 6 , Figure 6 This is a schematic diagram comparing the equivalent Curie inversion results with the real model; the inverted results can be... The interface is compared with a real model to verify the correctness and stability of the method in this application.

[0033] Furthermore, to illustrate The impact on the inversion results, refer to Figure 7 Inversion comparisons can be performed using different magnetic susceptibility structure assumptions under the same grid, initial model, regularization strategy, and iteration stopping criterion. For example, experiments can be used. Constructed Perform inversion and use linear approximation Perform inversion and use uniform The assumption is to perform an inversion. This involves comparing different assumptions. Compared to the real interface, it can represent The morphology has a first-order control effect on the mean depth level and undulation amplitude of the equivalent Curie; using rock physics Constructed It can reduce systematic bias and improve inversion stability and physical consistency.

[0034] Specifically, through S104, combining the first isotherm surface, surface boundary temperature, and Moho depth and its corresponding temperature, a piecewise temperature gradient model is constructed to reconstruct the regional three-dimensional crustal temperature structure. The equivalent Curie depth distribution obtained in step S103 is then used to... As the first isothermal surface, and given a temperature Further define the depth distribution of the Mohorovičić discontinuity. As a second isothermal surface, and given a temperature The It can be set to 850℃ or determined based on the regional thermal conditions. Surface boundary temperature. It can be set to a constant or a spatial distribution can be given as needed.

[0035] In one exemplary implementation, a piecewise linear temperature gradient model is used to reconstruct the temperature structure: for any horizontal position ,when At that time, the temperature satisfies ;when At that time, the temperature satisfies This yields a three-dimensional temperature volume. . Reference Figure 8 It can output slices at different depths, typical cross-sections, and isothermal surface morphology.

[0036] In another alternative implementation, the steady-state thermal conductivity equation can be solved under given thermal conductivity and internal heat source terms to further improve the physical consistency of the temperature structure.

[0037] Finally, the equivalent Curie depth distribution obtained from the S103 inversion is used as the first isotherm, and a temperature of 580℃ is assigned to this isotherm. The Moho is set as the second isotherm, and a temperature is assigned to the Moho based on the regional thermal state, typically set to 850℃ or determined specifically based on regional geothermal parameters such as geothermal heat flow and geothermal gradient. Simultaneously, the surface boundary temperature is obtained, which can be determined based on measured surface temperature, annual average surface temperature, or remotely sensed surface temperature data.

[0038] Outputs a three-dimensional temperature volume and its cross-sections, slices, and / or isothermal surfaces. The output may include a three-dimensional temperature volume. The results of horizontal slices at specified depths, profiles of specified survey lines, and isothermal surfaces of the 580℃ isotherm and the second isotherm (Moho isothermal boundary) are used to characterize the crustal thermal structure at a regional scale.

[0039] Optionally, the method for obtaining the equivalent Curie depth distribution specifically includes: The interface of abrupt change in underground magnetic susceptibility is used as the interface model to be inverted, and an asymptotic reference plane and its reference depth are set. A forward modeling operator for the magnetic anomaly spatial domain is established based on Cauchy surface integral; the forward modeling operator uses the depth-related magnetic susceptibility function to describe the difference in magnetic susceptibility between the upper and lower sides of the abrupt change in subsurface magnetic susceptibility interface. Determine the depth perturbation vector of the underground magnetic susceptibility abrupt change interface relative to the reference plane, and calculate the magnetic anomaly response generated by the underground magnetic susceptibility abrupt change interface in the form of a two-dimensional surface integral. Construct an objective functional that includes data fitting and stabilization terms, construct a model weighting matrix based on the derivative matrix of the forward operator, and perform weighted parameterization on the interface parameters. Based on the derivative matrix, the minimum value of the objective functional is iteratively solved using the reweighted regularized conjugate gradient algorithm to obtain the depth distribution of the underground magnetic susceptibility abrupt change interface, which serves as the equivalent Curie depth distribution.

[0040] Specifically, in this embodiment, the subsurface magnetic susceptibility abrupt change interface Γ is used as the interface model to be inverted, and an asymptotic reference plane P (reference depth) is set. A forward modeling operator for the magnetic anomaly spatial domain is established based on Cauchy surface integral, and the magnetic anomaly response generated by the interface Γ is calculated by replacing the three-dimensional volume integral with the two-dimensional surface integral form. Construct a Tikhonov objective functional that includes data fitting and stabilization terms, and construct a model weighting matrix based on the Fréchet derivative matrix of the forward operator. The interface parameters are weighted and parameterized. Based on the Fréchet derivative matrix, the minimum value of the Tikhonov objective functional is iteratively solved using the reweighted regularized conjugate gradient algorithm to obtain the depth distribution of the interface Γ, and this is used as the equivalent Curie depth distribution. First, the interface of abrupt change in underground magnetic susceptibility is used as the interface model to be inverted. This interface represents the depth interface at which magnetic minerals (mainly magnetite) in the crust undergo demagnetization. At the same time, an asymptotic reference plane and its reference depth are set. The reference depth is usually preset based on the average Curie depth of the target area or the regional geological background.

[0041] Secondly, a forward modeling operator for the magnetic anomaly spatial domain is established based on the Cauchy surface integral. The core function of this forward modeling operator is to use the aforementioned depth-related magnetic susceptibility function to describe the difference in magnetic susceptibility between the upper and lower sides of the interface of abrupt change in underground magnetic susceptibility, thereby separating the influence of interface undulations and magnetic susceptibility changes on magnetic anomalies.

[0042] Then, the depth perturbation vector of the subsurface magnetic susceptibility abrupt change interface relative to the reference plane is determined, and the magnetic anomaly response generated by the subsurface magnetic susceptibility abrupt change interface is calculated by replacing the traditional three-dimensional volume integral form with a two-dimensional surface integral. Based on this, a Tikhonov objective functional including data fitting terms and model stabilization terms is constructed. The model weighting matrix is ​​constructed based on the Fréchet derivative matrix of the forward operator, and the interface parameters are weighted and parameterized to balance the sensitivity differences of different interface depth parameters to the magnetic anomaly response.

[0043] Finally, based on the Fréchet derivative matrix, the minimum value of the objective functional is iteratively solved using the reweighted regularized conjugate gradient algorithm. The model parameters and search direction are dynamically updated in each iteration until the convergence condition is met, thus obtaining the depth distribution of the subsurface magnetic susceptibility abrupt change interface. This depth distribution is the equivalent Curie depth distribution. Through the above method, this application achieves efficient and stable inversion of the Curie depth.

[0044] Optionally, the forward modeling operator includes: calculating the X-direction component, Y-direction component, and Z-direction component of the magnetic field generated at any observation point through the abrupt change in underground magnetic susceptibility. The calculation methods for the X-direction component, Y-direction component, and Z-direction component of the magnetic field include: Based on the spatial relative positional relationship between any observation point and any point on the interface of abrupt change in underground magnetic susceptibility, and combining the projection component of the uniform magnetization intensity vector in the normal direction outside the interface, the perpendicular component of the uniform magnetization intensity vector, and the non-normalized geometric factor of the normal direction outside the interface, Cauchy-type surface integral integrand terms corresponding to the X-direction component, Y-direction component, and Z-direction component of the magnetic field are constructed respectively. When calculating the magnetic field components in each direction, the coordinate difference between any observation point and any point on the interface in the corresponding direction is used as the position response term, and the interface geometry is coupled with the magnetization information to characterize the influence of the abrupt change in underground magnetic susceptibility on the magnetic field components in each direction. Two-dimensional surface integrals were performed on the underground magnetic susceptibility abrupt interface to obtain the X-direction component, Y-direction component, and Z-direction component of the magnetic field at any observation point. Wherein, the underground magnetic susceptibility abrupt change interface is used to characterize the spatial boundary where the magnetic susceptibility changes abruptly; the projection component of the uniform magnetization intensity vector in the normal direction outside the interface is used to characterize the normal magnetization contribution; the vertical component of the uniform magnetization intensity vector is used to characterize the vertical magnetization contribution; and the non-normalized geometric factor in the normal direction outside the interface is used to characterize the local geometry of the interface. Specifically, the X-direction component, Y-direction component, and Z-direction component of the magnetic field satisfy the following Cauchy surface integral expression: ; ; ; in, For the X-direction component of the magnetic field, This represents the Y-direction component of the magnetic field. This represents the Z-direction component of the magnetic field. This is an interface where the underground magnetic susceptibility changes abruptly. For any point on the interface of abrupt change in underground magnetic susceptibility, For any observation point, Let be the projection vector of the uniform magnetization vector onto the normal direction outside the interface. They are respectively The components in the X, Y, and Z directions are related to the difference in magnetic susceptibility on both sides of the interface and the normal direction of the interface. The vertical component of the uniform magnetization vector; The non-normalized geometric factor for the outward normal of the interface is used to describe the local tilt characteristics of the interface. Through the above three-component expression, this application can fully calculate the total magnetic field vector response generated at any observation point in space by an interface with abrupt changes in magnetic susceptibility, providing an accurate forward modeling basis for subsequent inversion.

[0045] Optionally, the method for determining the non-normalized geometric factor includes: The elevation difference between the subsurface magnetic susceptibility abrupt change interface and the reference depth is defined as an elevation function to characterize the fluctuations of the subsurface magnetic susceptibility abrupt change interface. The negative rate of change of the elevation function along the X direction is determined as the non-normalized geometric factor in the X direction; The negative rate of change of the elevation function along the Y direction is determined as the non-normalized geometric factor in the Y direction; The nonnormalized geometric factor in the Z direction is set to a unit constant (set to 1 in this embodiment); Wherein, the reference depth is the depth of the asymptotic reference plane; the elevation function is determined by the Z coordinate of any point on the abrupt change interface of underground magnetic susceptibility and the reference depth; the non-normalized geometric factors in the X and Y directions respectively characterize the local tilt changes of the abrupt change interface of underground magnetic susceptibility in the X and Y directions.

[0046] Specifically, the non-normalized geometric factor The following conditions must be met: ; ; ; ; in, This is an elevation function of the subsurface magnetic susceptibility abrupt change interface relative to the reference depth. These represent the X, Y, and Z coordinates of any point on the interface of abrupt change in underground magnetic susceptibility. For reference depth.

[0047] This application simplifies the normal calculation process in surface integrals, avoids the additional computational complexity introduced by the normalization factor, and ensures the accurate description of the interface geometry by the forward operator. As a non-normalized geometric factor of the interface outward normal, it directly participates in the calculation of the magnetic field component in the aforementioned Cauchy surface integral expression, and is used to pass the interface normal information to the integral kernel function.

[0048] Optionally, the Cauchy surface integral is achieved through regular grid discretization, and the process of regular grid discretization specifically includes: The horizontal projection of the underground magnetic susceptibility abrupt interface is divided into multiple grid cells, and the local interface is approximated in planar form by the elevation function at the center of each grid cell. The contribution of each grid cell is calculated by forward modeling at the discrete observation points of the target to obtain the total field anomaly.

[0049] Furthermore, the planar approximation of the local interface is achieved using the following formula: ; in, For grid cell indexing, Let the interface depth of the k-th grid cell be . and Let be the non-normalized geometric factors of the k-th mesh element in the X and Y directions, respectively. The coordinates of the cell center For reference depth; Forward modeling is performed using the following formula: ; in, for The overall field is abnormal. Let be the target discrete observation point, and let represent the nth observation point after discretization. The total number of grid cells. This is the forward kernel function.

[0050] Specifically, in this embodiment, the horizontal projection region of the underground magnetic susceptibility abrupt change interface is divided into multiple regular rectangular grid cells, each with uniform dimensions in the X and Y directions. Within each grid cell, a planar approximation strategy is used to simplify the description of the local interface morphology. That is, it is assumed that the interface depth within the cell varies linearly with the horizontal coordinate. This planar approximation uses the elevation function at the center of the cell as a reference and characterizes the tilting features of the interface through a local dip angle coefficient.

[0051] Then, forward modeling is performed at the target discrete observation point, and the magnetic field contribution of each grid cell to the observation point is calculated independently. Finally, the contributions of all grid cells are summed to obtain the total field anomaly at the observation point.

[0052] Optionally, the method for determining the forward kernel function includes: Based on the spatial relative positional relationship between the center of the kth grid cell and the nth target discrete observation point, the non-normalized geometric factor of the kth grid cell, the projection components and vertical components of the uniform magnetization intensity vector in each direction, the cell response of the kth grid cell to the nth target discrete observation point is constructed. The normal response caused by the normal projection component is determined by combining the unit vector of the induced geomagnetic field direction with the coordinate differences in the X, Y, and Z directions between the target discrete observation point and the center of the k-th grid cell; the unit vector of the induced geomagnetic field direction is used to characterize the geomagnetic field direction of the region. The vertical response caused by the vertical magnetization contribution is determined by combining the vertical component of the uniform magnetization intensity vector, the non-normalized geometric factor, and the coordinate differences in the X, Y, and Z directions between the target discrete observation point and the center of the k-th grid cell. The unit response, normal response, and vertical response are superimposed to obtain the forward kernel function.

[0053] The forward kernel function satisfies the following formula: ; in, It is the product of the dimensions of the mesh element in the X and Y directions. The cell center of the k-th point, Let the coordinates be the discrete observation points of any target. The coordinates of the element center include the Z-direction component. This is the unit vector representing the direction of the induced geomagnetic field.

[0054] The forward modeling kernel function expression in this embodiment integrates geometric parameters (element size, element center coordinates, observation point coordinates), physical property parameters (normal magnetization component, volume partial factor, magnetization direction), and geomagnetic field parameters into a single analytical formula, thereby achieving efficient numerical calculation of the forward modeling response.

[0055] Optionally, the process for determining the objective functional is as follows: Based on data weighting matrix and Cauchy integral forward modeling operator Obtain the weighted forward response Based on the data weighting matrix and the observed magnetic anomaly vectors from actual measurements (Inversion input data) yields weighted observation data The data fitting term of the objective functional is determined based on the weighted forward response and weighted observation data; Based on the model weighting matrix Discrete elevation parameter vector of the interface of the inverted target Obtain the weighted model parameters Based on the model weighting matrix and prior model vectors Obtain the weighted prior model parameters The model stabilization term of the objective functional is determined based on the weighted model parameters and the weighted prior model parameters. The objective functional is obtained by adding the data fitting term to the model stabilization term, as shown in the following formula: ; in, Weighted matrix of data, For Cauchy-type integral forward modeling operators, This is the magnetic anomaly observation vector. This is the discrete elevation parameter vector of the interface. This represents the prior model vector.

[0056] The data weighting matrix is ​​determined based on an identity matrix with the same dimension as the magnetic anomaly observation vector; the model weighting matrix... The model weighting matrix is ​​defined in the model space and is used to weight the deviation of the current model parameters from the prior model parameters. The model weighting matrix is ​​determined based on the Fréchet derivative matrix, and the prior model vector is obtained based on prior geological or geophysical information. The Fréchet derivative matrix was obtained analytically for the observation points. The overall site is abnormal regarding the first Unit Elevation partial derivatives satisfy: ; The Fréchet derivative matrix is ​​obtained analytically and is dynamically updated during the inversion iteration process. Using an analytical method avoids computational errors and step size selection issues caused by numerical differencing, and the dynamic updating of F during the inversion iteration process ensures that the weighting matrix always matches the local sensitivity characteristics of the current model. Through the construction of the objective functional and the weighting matrix, this application establishes a balance mechanism between data fitting and model stability.

[0057] Optionally, the weighted parameterization process specifically includes: Through the model weighting matrix Transform the discrete elevation parameter vector m of the interface to obtain the weighted model parameter vector in the weighted space. ; In the weighted space, based on the inverse matrix of the model weighting matrix For Cauchy integral forward modeling operators By performing the transformation, we obtain the weighted forward modeling operator for the weighted space. ; through the data weighting matrix The magnetic anomaly observation vector d is weighted to obtain the weighted observation data vector. ; through the model weighting matrix For prior model vectors We perform weighting to obtain the weighted prior model vector. ; The data fitting term of the weighted spatial objective functional is determined based on the weighted forward operator, the weighted model parameter vector, and the weighted observation data vector; the stabilization term of the model weighted spatial objective functional is determined based on the weighted model parameter vector, the weighted prior model vector, and the regularization coefficient. Based on the data fitting term and stabilization term, the objective functional is rewritten as a weighted spatial objective functional.

[0058] Specifically, the weighted parameterization includes: ; The objective functional is then rewritten in the weighted space as follows: ; in, α is the regularization parameter.

[0059] Optionally, the calculation process of the reweighted regularized conjugate gradient algorithm includes: Based on the weighted model parameter vector of the nth iteration Combined with weighted forward modeling operators With weighted observation data vector The weighted residual vector is calculated. ; Based on the Fréchet derivative matrix of the nth iteration and the weighted residual vector Combined with the regularization parameter of the nth iteration Weighted model parameter vector and weighted prior model parameters The current fastest upward direction is calculated. ; The conjugate coefficients are determined based on the modulus of the current steepest ascent direction in the nth iteration and the modulus of the steepest ascent direction in the previous iteration. The conjugate search direction is updated based on the conjugate coefficient, the current steepest ascent direction, and the conjugate search direction of the previous iteration. ; Based on the conjugate search direction The fastest upward direction The Fréchet derivative matrix of the nth iteration and the regularization parameter of the nth iteration The step size is obtained by line search calculation. ; Based on the The weighted model parameter vector of the next iteration Step length and conjugate search direction The updated weighted model parameter vector for the (n+1)th iteration is calculated. ; Based on the weighted model parameter vector and the inverse of the model weighting matrix obtained in the (n+1)th iteration, the updated interface discrete elevation parameter vector is calculated. Repeat the above iterative steps until the convergence condition is met, and output the final interface elevation parameter vector as the inversion result.

[0060] In the nth iteration, calculate the weighted residual vector: ; Calculate the steepest ascent direction: ; Update conjugate search direction: ; Determine the step size using line search: ; Update the model: ; ; Optionally, the regularization parameters are updated using an adaptive update strategy, which includes: After the nth iteration, based on the weighted model parameter vector of the nth iteration... The weighted model parameter vector for the (n+1)th iteration and weighted prior model vectors Determine the weighted model deviation vector for the nth iteration. Deviation vector from the weighted model in the (n+1)th iteration ; The weighted model deviation vector in the (n+1)th iteration Deviation vector from the weighted model in the nth iteration The norm square ratio is used to determine the deviation ratio coefficient. ; If the deviation ratio coefficient is less than or equal to the preset threshold, then the current regularization parameter remains unchanged for the next iteration; If the deviation ratio coefficient is greater than the preset threshold, the current regularization parameter is divided by the deviation ratio coefficient to obtain the updated regularization parameter for the next iteration; Specifically, the regularization parameter α is adaptively updated as follows: ; ; It should be noted that the regularization parameters in this embodiment are dynamically adjusted using an adaptive update strategy to balance the relative weights of the data fitting term and the model stabilization term at different iteration stages. Specifically, after each iteration, the deviation vector between the weighted model parameters of the current iteration and the prior model vector, as well as the deviation vector between the weighted model parameters of the next iteration and the prior model vector, are first calculated.

[0061] The judgment coefficient for adjusting the regularization parameter is calculated based on the ratio of the square of the magnitude of the deviation vector in the next iteration to the square of the magnitude of the deviation vector in the current iteration. This judgment coefficient reflects the trend of model deviation between two adjacent iterations: if γ is small, it indicates that the model has stabilized; if γ is large, it indicates that the model is still adjusting significantly. Next, the judgment coefficient γ is compared with a preset threshold: if γ is less than or equal to the preset threshold, it indicates that the model deviation has not increased or has converged, and the current regularization parameter is kept unchanged and directly used for the next iteration; if γ is greater than the preset threshold, it indicates that the model deviation has increased significantly, and the current regularization parameter is divided by the judgment coefficient to obtain the updated regularization parameter, which is used for the next iteration. In this embodiment, when the model deviation increases, the regularization parameter is appropriately reduced, thereby reducing the weight of the model stabilization term and allowing the model to update more significantly to fit the observed data; when the model deviation is stable or decreases, the regularization parameter is kept unchanged to maintain the existing balance between data fitting and model stability. Through the aforementioned adaptive update strategy, this application achieves dynamic adjustment of the regularization parameter as the inversion process progresses, avoiding the tedious process of manual parameter tuning and improving the automation and convergence stability of the inversion algorithm.

[0062] Optionally, the reconstructed three-dimensional crustal temperature structure and output of temperature distribution results include: Calculate the temperature gradient between the surface and the equivalent Curie based on the surface boundary temperature and the first isothermal surface. Calculate the temperature gradient between the equivalent Curie and the Mohorovičić discontinuity based on the first and second isotherms; The temperature at any target depth is interpolated according to the piecewise temperature gradient to obtain and output a three-dimensional temperature volume and the cross-section, slice and / or isothermal surface results of the three-dimensional temperature volume.

[0063] Specifically, the reconstruction of the three-dimensional crustal temperature structure includes: Based on surface boundary temperature With the first isothermal surface Calculate the temperature gradient between the Earth's surface and the equivalent Curie; Based on the first isothermal surface With the second isothermal surface Calculate the temperature gradient between the equivalent Curie and the Mohorovičić discontinuity; Based on piecewise temperature gradients for arbitrary depths The temperature at a given location is interpolated to obtain and output a three-dimensional temperature volume. And its cross-sectional, slice and / or isothermal surface results.

[0064] Reference Figure 9 This application also provides a device for reconstructing crustal temperature distribution based on magnetic anomalies, comprising: Data acquisition module 910 is used to acquire magnetic anomaly grid data of the target area; The magnetic susceptibility constraint construction module 920 is used to obtain the magnetic susceptibility-temperature curve based on rock physics experiments, construct a depth-related magnetic susceptibility function in combination with the preset temperature-depth relationship, and use the depth-related magnetic susceptibility function as the magnetic susceptibility structure constraint for the forward and inverse modeling of magnetic anomalies. Inversion module 930 is used to input the magnetic anomaly grid data and depth-related magnetic susceptibility function into a spatial domain forward and inversion framework based on Cauchy surface integral, construct an objective functional including data fitting terms and model stabilization terms, and iteratively solve it using a reweighted regularized conjugate gradient algorithm to invert and obtain the equivalent Curie depth distribution. The temperature reconstruction module 940 is used to reconstruct the three-dimensional crustal temperature structure of the region and output the temperature distribution results by taking the equivalent Curie depth distribution as the first isothermal surface, combining it with the surface boundary temperature and the temperature boundary with the Moho surface as the second isothermal surface.

[0065] In this embodiment, the equivalent Curie depth distribution can be obtained as follows: The subsurface magnetic susceptibility abrupt change interface is used as the interface model to be inverted, and an asymptotic reference plane and its reference depth are set; a forward modeling operator for the magnetic anomaly spatial domain is established based on Cauchy surface integral; the forward modeling operator uses the depth-related magnetic susceptibility function to describe the magnetic susceptibility difference between the upper and lower sides of the subsurface magnetic susceptibility abrupt change interface; the depth perturbation vector of the subsurface magnetic susceptibility abrupt change interface relative to the reference plane is determined, and the magnetic anomaly response generated by the subsurface magnetic susceptibility abrupt change interface is calculated in the form of a two-dimensional surface integral; a target functional including a data fitting term and a stabilization term is constructed, and a model weighting matrix is ​​constructed based on the derivative matrix of the forward modeling operator to perform weighted parameterization processing on the interface parameters; based on the derivative matrix, the minimum value of the target functional is iteratively solved using a reweighted regularized conjugate gradient algorithm to obtain the depth distribution of the subsurface magnetic susceptibility abrupt change interface, which serves as the equivalent Curie depth distribution.

[0066] In this embodiment of the application, the forward modeling operator includes: calculating the X-direction component, Y-direction component, and Z-direction component of the magnetic field generated at any observation point by the interface of abrupt change in underground magnetic susceptibility; the X-direction component, Y-direction component, and Z-direction component of the magnetic field are obtained through the following method: Calculate the X-direction component, Y-direction component, and Z-direction component of the magnetic field generated by the abrupt change in underground magnetic susceptibility at any observation point. Based on the spatial relative positional relationship between arbitrary observation and any point on the interface of abrupt change in underground magnetic susceptibility, and combining the projection component of the uniform magnetization intensity vector in the normal direction outside the interface, the perpendicular component of the uniform magnetization intensity vector, and the non-normalized geometric factor of the normal direction outside the interface, Cauchy-type surface integral integrand terms corresponding to the X-direction component, Y-direction component, and Z-direction component of the magnetic field are constructed respectively. When calculating the magnetic field components in each direction, the coordinate difference between any observation point and any point on the interface in the corresponding direction is used as the position response term, and the interface geometry is coupled with the magnetization information to characterize the influence of the abrupt change in underground magnetic susceptibility on the magnetic field components in each direction. Two-dimensional surface integrals were performed on the underground magnetic susceptibility abrupt interface to obtain the X-direction component, Y-direction component, and Z-direction component of the magnetic field at any observation point. Wherein, the underground magnetic susceptibility abrupt change interface is used to characterize the spatial boundary where the magnetic susceptibility changes abruptly; the projection component of the uniform magnetization intensity vector in the normal direction outside the interface is used to characterize the normal magnetization contribution; the vertical component of the uniform magnetization intensity vector is used to characterize the vertical magnetization contribution; and the non-normalized geometric factor in the normal direction outside the interface is used to characterize the local geometry of the interface.

[0067] In this embodiment of the application, the inversion module is further used to determine the non-normalized geometric factor, and the method for determining the non-normalized geometric factor includes: The elevation difference between the subsurface magnetic susceptibility abrupt change interface and the reference depth is defined as an elevation function to characterize the fluctuations of the subsurface magnetic susceptibility abrupt change interface. The negative rate of change of the elevation function along the X direction is determined as the non-normalized geometric factor in the X direction; The negative rate of change of the elevation function along the Y direction is determined as the non-normalized geometric factor in the Y direction; Set the nonnormalized geometric factor in the Z direction to a unit constant; Wherein, the reference depth is the depth of the asymptotic reference plane; the elevation function is determined by the Z coordinate of any point on the abrupt change interface of underground magnetic susceptibility and the reference depth; the non-normalized geometric factors in the X and Y directions respectively characterize the local tilt changes of the abrupt change interface of underground magnetic susceptibility in the X and Y directions.

[0068] In this embodiment of the application, the Cauchy surface integral is implemented by regular mesh discretization, and the process of regular mesh discretization specifically includes: The horizontal projection of the underground magnetic susceptibility abrupt interface is divided into multiple grid cells, and the local interface is approximated in planar form by the elevation function at the center of each grid cell. The contribution of each grid cell is calculated by forward modeling at the discrete observation points of the target to obtain the total field anomaly.

[0069] In this embodiment, the total field anomaly is calculated based on the discretized nth observation point, the total number of grid cells, and the forward kernel function; the method for determining the forward kernel function includes: Based on the spatial relative positional relationship between the center of the kth grid cell and the nth target discrete observation point, the non-normalized geometric factor of the kth grid cell, the projection components and vertical components of the uniform magnetization intensity vector in each direction, the cell response of the kth grid cell to the nth target discrete observation point is constructed. The normal response caused by the normal projection component is determined by combining the unit vector of the induced geomagnetic field direction with the coordinate differences in the X, Y, and Z directions between the target discrete observation point and the center of the k-th grid cell; the unit vector of the induced geomagnetic field direction is used to characterize the geomagnetic field direction of the region. The vertical response caused by the vertical magnetization contribution is determined by combining the vertical component of the uniform magnetization intensity vector, the non-normalized geometric factor, and the coordinate differences in the X, Y, and Z directions between the target discrete observation point and the center of the k-th grid cell. The unit response, normal response, and vertical response are superimposed to obtain the forward kernel function.

[0070] In this embodiment of the application, the inversion module is further configured to determine the target functional in the following manner: Based on data weighting matrix and Cauchy integral forward modeling operator Obtain the weighted forward response Based on the data weighting matrix and the observed magnetic anomaly vectors from actual measurements Obtain weighted observation data The data fitting term of the objective functional is determined based on the weighted forward response and weighted observation data; Based on the model weighting matrix Discrete elevation parameter vector of the interface of the inverted target Obtain the weighted model parameters Based on the model weighting matrix and prior model vectors Obtain the weighted prior model parameters The model stabilization term of the objective functional is determined based on the weighted model parameters and the weighted prior model parameters. The target functional is obtained by adding the data fitting term to the model stabilization term.

[0071] The data weighting matrix is ​​determined based on an identity matrix with the same dimension as the magnetic anomaly observation vector; the model weighting matrix... The model weighting matrix is ​​defined in the model space and is used to weight the deviation of the current model parameters from the prior model parameters. The model weighting matrix is ​​determined based on the Fréchet derivative matrix, and the prior model vector is obtained based on prior geological or geophysical information. The Fréchet derivative matrix is ​​obtained analytically; the Fréchet derivative matrix is ​​dynamically updated during the inversion iteration process.

[0072] In this embodiment of the application, the inversion module is further used to implement weighted parameterization, the process of which specifically includes: Through the model weighting matrix Transform the discrete elevation parameter vector m of the interface to obtain the weighted model parameter vector in the weighted space. ; In the weighted space, based on the inverse matrix of the model weighting matrix For Cauchy integral forward modeling operators By performing the transformation, we obtain the weighted forward modeling operator for the weighted space. ; through the data weighting matrix The magnetic anomaly observation vector d is weighted to obtain the weighted observation data vector. ; through the model weighting matrix For prior model vectors We perform weighting to obtain the weighted prior model vector. ; The data fitting term of the weighted spatial objective functional is determined based on the weighted forward operator, the weighted model parameter vector, and the weighted observation data vector; the stabilization term of the model weighted spatial objective functional is determined based on the weighted model parameter vector, the weighted prior model vector, and the regularization coefficient. Based on the data fitting term and stabilization term, the objective functional is rewritten as a weighted spatial objective functional.

[0073] In this embodiment of the application, the inversion module is further used to calculate the reweighted regularized conjugate gradient algorithm, the process of which includes: Based on the weighted model parameter vector of the nth iteration Combined with forward modeling operators With weighted observation data vector The weighted residual vector is calculated. ; Based on the Fréchet derivative matrix of the nth iteration and the weighted residual vector Combined with the regularization parameter of the nth iteration Weighted model parameter vector and weighted prior model parameters The current fastest upward direction is calculated. ; The conjugate coefficients are determined based on the modulus of the current steepest ascent direction in the nth iteration and the modulus of the steepest ascent direction in the previous iteration. The conjugate search direction is updated based on the conjugate coefficient, the current steepest ascent direction, and the conjugate search direction of the previous iteration. ; Based on the conjugate search direction The fastest upward direction The Fréchet derivative matrix of the nth iteration and the regularization parameter of the nth iteration The step size is obtained by line search calculation. ; Based on the weighted model parameter vector of the nth iteration Step length and conjugate search direction The updated weighted model parameter vector for the (n+1)th iteration is calculated. ; Based on the weighted model parameter vector and the inverse of the model weighting matrix obtained in the (n+1)th iteration, the updated interface discrete elevation parameter vector is calculated. Repeat the above iterative steps until the convergence condition is met, and output the final interface elevation parameter vector as the inversion result.

[0074] In this embodiment of the application, the regularization parameter is updated using an adaptive update strategy, which is described in the following steps: After the nth iteration, based on the weighted model parameter vector of the nth iteration... The weighted model parameter vector for the (n+1)th iteration and weighted prior model vectors Determine the weighted model deviation vector for the nth iteration. Deviation vector from the weighted model in the (n+1)th iteration ; The weighted model deviation vector in the (n+1)th iteration Deviation vector from the weighted model in the nth iteration The norm square ratio is used to determine the deviation ratio coefficient. ; If the deviation ratio coefficient is less than or equal to the preset threshold, then the current regularization parameter remains unchanged for the next iteration; If the deviation ratio coefficient is greater than the preset threshold, the current regularization parameter is divided by the deviation ratio coefficient to obtain the updated regularization parameter for the next iteration; In this embodiment of the application, the temperature reconstruction module is specifically used for: Calculate the temperature gradient between the surface and the equivalent Curie based on the surface boundary temperature and the first isothermal surface. Calculate the temperature gradient between the equivalent Curie and the Mohorovičić discontinuity based on the first and second isotherms; The temperature at any target depth is interpolated according to the piecewise temperature gradient to obtain and output a three-dimensional temperature volume and the cross-section, slice and / or isothermal surface results of the three-dimensional temperature volume.

[0075] In this embodiment of the application, the data acquisition module is specifically used for: Magnetic anomaly data is acquired and then subjected to regular gridding and preprocessing to obtain magnetic anomaly grid data of the target area. The preprocessing process includes one or more of the following: outlier removal, missing value imputation, mean removal, linear trend removal, boundary extension, upward extension, pole reduction, frequency domain filtering, and / or equivalent source reconstruction.

[0076] Reference Figure 10 Based on the methods in the above embodiments, this application provides an electronic device that may include: a processor 11, a communications interface 12, a memory 13, and a communication bus 14, wherein the processor 11, the communications interface 12, and the memory 13 communicate with each other through the communication bus 14. The processor 11 may call logical instructions in the memory 13 to execute the methods in the above embodiments.

[0077] Furthermore, the logical instructions in the aforementioned memory 13 can be implemented as software functional units and, when sold or used as independent products, can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this application.

[0078] Based on the methods in the above embodiments, this application provides a computer-readable storage medium storing a computer program that, when run on a processor, causes the processor to execute the methods in the above embodiments.

[0079] Based on the methods in the above embodiments, this application provides a computer program product that, when run on a processor, causes the processor to execute the methods in the above embodiments.

[0080] It is understood that the processor in the embodiments of this application can be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, transistor logic devices, hardware components, or any combination thereof. A general-purpose processor can be a microprocessor or any conventional processor.

[0081] The method steps in this application embodiment can be implemented in hardware or by a processor executing software instructions. The software instructions can consist of corresponding software modules, which can be stored in random access memory (RAM), flash memory, read-only memory (ROM), programmable read-only memory (PROM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), registers, hard disks, portable hard disks, CD-ROMs, or any other form of storage medium known in the art. An exemplary storage medium is coupled to the processor, enabling the processor to read information from and write information to the storage medium. Of course, the storage medium can also be a component of the processor. The processor and the storage medium can reside in an ASIC.

[0082] In the above embodiments, implementation can be achieved entirely or partially through software, hardware, firmware, or any combination thereof. When implemented using software, it can be implemented entirely or partially as a computer program product. The computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, all or part of the processes or functions described in the embodiments of this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted through the computer-readable storage medium. The computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wired (e.g., coaxial cable, fiber optic, digital subscriber line (DSL)) or wireless (e.g., infrared, wireless, microwave, etc.) means. The computer-readable storage medium can be any available medium that a computer can access or a data storage device such as a server or data center that integrates one or more available media. The available medium can be a magnetic medium (e.g., floppy disk, hard disk, magnetic tape), an optical medium (e.g., DVD), or a semiconductor medium (e.g., solid-state disk (SSD)).

[0083] It is understood that the various numerical designations used in the embodiments of this application are merely for the convenience of description and are not intended to limit the scope of the embodiments of this application.

[0084] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the scope of protection of this application.

Claims

1. A method for reconstructing the temperature distribution of the earth's crust based on magnetic anomalies, characterized in that, include: Acquire magnetic anomaly grid data for the target area; Based on the magnetic susceptibility-temperature curve obtained from rock physics experiments, a depth-dependent magnetic susceptibility function is constructed by combining the preset temperature-depth relationship. The depth-dependent magnetic susceptibility function is used as the magnetic susceptibility structure constraint for forward and inverse magnetic anomaly modeling. A forward modeling operator for the magnetic anomaly spatial domain is established based on Cauchy surface integral. The forward modeling operator uses the depth-dependent magnetic susceptibility function to describe the magnetic susceptibility difference between the upper and lower sides of the abrupt magnetic susceptibility interface in the subsurface. The magnetic anomaly grid data and the depth-related magnetic susceptibility function are input into the spatial domain forward and inverse modeling framework established based on Cauchy surface integral. An objective functional including data fitting terms and model stabilization terms is constructed. The reweighted regularized conjugate gradient algorithm is used to iteratively solve the problem and obtain the equivalent Curie depth distribution through inversion. Using the equivalent Curie depth distribution as the first isotherm, combined with the surface boundary temperature and the temperature boundary with the Moho surface as the second isotherm, the three-dimensional crustal temperature structure of the region is reconstructed and the temperature distribution results are output. The forward modeling operator includes: calculating the X-direction component, Y-direction component, and Z-direction component of the magnetic field generated at any observation point by the interface of abrupt change in underground magnetic susceptibility; the calculation method for the X-direction component, Y-direction component, and Z-direction component of the magnetic field includes: Based on the spatial relative positional relationship between any observation point and any point on the interface of abrupt change in underground magnetic susceptibility, and combining the projection component of the uniform magnetization intensity vector in the normal direction outside the interface, the perpendicular component of the uniform magnetization intensity vector, and the non-normalized geometric factor of the normal direction outside the interface, Cauchy-type surface integral integrand terms corresponding to the X-direction component, Y-direction component, and Z-direction component of the magnetic field are constructed respectively. When calculating the magnetic field components in each direction, the coordinate difference between any observation point and any point on the interface in the corresponding direction is used as the position response term, and the interface geometry is coupled with the magnetization information to characterize the influence of the abrupt change in underground magnetic susceptibility on the magnetic field components in each direction. Two-dimensional surface integrals were performed on the underground magnetic susceptibility abrupt interface to obtain the X-direction component, Y-direction component, and Z-direction component of the magnetic field at any observation point. The method for determining the non-normalized geometric factor includes: The elevation difference between the subsurface magnetic susceptibility abrupt change interface and the reference depth is defined as an elevation function to characterize the fluctuations of the subsurface magnetic susceptibility abrupt change interface. The negative rate of change of the elevation function along the X direction is determined as the non-normalized geometric factor in the X direction; The negative rate of change of the elevation function along the Y direction is determined as the non-normalized geometric factor in the Y direction; Set the nonnormalized geometric factor in the Z direction to a unit constant.

2. The method of claim 1, wherein, The method for obtaining the equivalent Curie depth distribution specifically includes: The interface of abrupt change in underground magnetic susceptibility is used as the interface model to be inverted, and an asymptotic reference plane and its reference depth are set. A forward modeling operator for the magnetic anomaly spatial domain is established based on Cauchy surface integral; the forward modeling operator uses the depth-related magnetic susceptibility function to describe the difference in magnetic susceptibility between the upper and lower sides of the abrupt change in subsurface magnetic susceptibility interface. Determine the depth perturbation vector of the underground magnetic susceptibility abrupt change interface relative to the reference plane, and calculate the magnetic anomaly response generated by the underground magnetic susceptibility abrupt change interface in the form of a two-dimensional surface integral. Construct an objective functional that includes data fitting and stabilization terms, construct a model weighting matrix based on the derivative matrix of the forward operator, and perform weighted parameterization on the interface parameters. Based on the derivative matrix, the minimum value of the objective functional is iteratively solved using the reweighted regularized conjugate gradient algorithm to obtain the depth distribution of the underground magnetic susceptibility abrupt change interface, which serves as the equivalent Curie depth distribution.

3. The method according to claim 1, characterized in that, The Cauchy surface integral is achieved through regular grid discretization, which specifically includes: The horizontal projection of the underground magnetic susceptibility abrupt interface is divided into multiple grid cells, and the local interface is approximated in planar form within each grid cell by the elevation function at the center of the cell. The contribution of each grid cell is calculated by forward modeling at the discrete observation points of the target to obtain the total field anomaly.

4. The method according to claim 3, characterized in that, The total field anomaly is calculated based on the discretized nth observation point, the total number of grid cells, and the forward kernel function; The method for determining the forward kernel function includes: Based on the spatial relative positional relationship between the center of the kth grid cell and the nth target discrete observation point, the non-normalized geometric factor of the kth grid cell, the projection components and vertical components of the uniform magnetization intensity vector in each direction, the cell response of the kth grid cell to the nth target discrete observation point is constructed. The normal response caused by the normal projection component is determined by combining the unit vector of the induced geomagnetic field direction with the coordinate differences in the X, Y, and Z directions between the target discrete observation point and the center of the k-th grid cell; the unit vector of the induced geomagnetic field direction is used to characterize the geomagnetic field direction of the region. The vertical response caused by the vertical magnetization contribution is determined by combining the vertical component of the uniform magnetization intensity vector, the non-normalized geometric factor, and the coordinate differences in the X, Y, and Z directions between the target discrete observation point and the center of the k-th grid cell. The unit response, normal response, and vertical response are superimposed to obtain the forward kernel function.

5. The method according to claim 2, characterized in that, The construction process of the objective functional includes: The weighted forward response is obtained based on the data weighting matrix and the Cauchy integral forward modeling operator; the weighted observation data is obtained based on the data weighting matrix and the actual measured magnetic anomaly observation vector; and the data fitting term of the objective functional is determined based on the weighted forward response and the weighted observation data. The weighted model parameters are obtained based on the model weighting matrix and the discrete elevation parameter vector of the interface of the inverted target. The weighted prior model parameters are obtained based on the model weighting matrix and the prior model vector. The model stabilization term of the target functional is determined based on the weighted model parameters and the weighted prior model parameters. The target functional is obtained by adding the data fitting term to the model stabilization term. The data weighting matrix is ​​determined based on an identity matrix with the same dimension as the magnetic anomaly observation vector; the model weighting matrix is ​​a weighting matrix defined in the model space, used to weight the deviation of the current model parameters from the prior model parameters; the model weighting matrix is ​​determined based on the Fréchet derivative matrix, and the prior model vector is obtained based on prior geological or geophysical information. The Fréchet derivative matrix is ​​obtained analytically; the Fréchet derivative matrix is ​​dynamically updated during the inversion iteration process.

6. The method according to claim 5, characterized in that, The weighted parameterization process specifically includes: The discrete elevation parameter vector of the interface is transformed by the model weighting matrix to obtain the weighted model parameter vector in the weighted space; In the weighted space, the Cauchy integral forward modeling operator is transformed based on the inverse matrix of the model weighting matrix to obtain the weighted forward modeling operator of the weighted space; the magnetic anomaly observation vector is weighted using the data weighting matrix to obtain the weighted observation data vector; the prior model vector is weighted using the model weighting matrix to obtain the weighted prior model vector. The data fitting term of the weighted spatial objective functional is determined based on the weighted forward operator, the weighted model parameter vector, and the weighted observation data vector; the stabilization term of the model weighted spatial objective functional is determined based on the weighted model parameter vector, the weighted prior model vector, and the regularization coefficient.

7. The method according to claim 6, characterized in that, The calculation process of the reweighted regularized conjugate gradient algorithm includes: The weighted residual vector is calculated based on the weighted model parameter vector of the nth iteration, combined with the weighted forward modeling operator and the weighted observation data vector. Based on the Fréchet derivative matrix of the nth iteration and the weighted residual vector, combined with the regularization parameter, weighted model parameter vector and weighted prior model parameter of the nth iteration, the current steepest ascent direction is calculated. The conjugate coefficient is determined based on the modulus of the current steepest ascent direction in the nth iteration and the modulus of the steepest ascent direction in the previous iteration. The conjugate search direction is then updated based on the conjugate coefficient, the current steepest ascent direction, and the conjugate search direction of the previous iteration. The step size is calculated by line search based on the conjugate search direction, the steepest ascent direction, the Fréchet derivative matrix of the nth iteration, and the regularization parameter of the nth iteration. Based on the weighted model parameter vector, step size, and conjugate search direction of the nth iteration, the updated weighted model parameter vector for the (n+1)th iteration is calculated. Based on the weighted model parameter vector and the inverse of the model weighting matrix obtained in the (n+1)th iteration, the updated interface discrete elevation parameter vector is calculated. The above iterative steps are repeated until the convergence condition is met, and the final interface elevation parameter vector is output as the inversion result.

8. The method according to claim 7, characterized in that, The regularization parameters are updated using an adaptive update strategy, which includes: After the nth iteration, based on the weighted model parameter vector of the nth iteration, the weighted model parameter vector of the (n+1)th iteration, and the weighted prior model vector, the weighted model deviation vector of the nth iteration and the weighted model deviation vector of the (n+1)th iteration are determined respectively. The deviation ratio coefficient is determined by the ratio of the norm squared of the deviation vector of the weighted model in the (n+1)th iteration to that in the nth iteration. If the deviation ratio coefficient is less than or equal to the preset threshold, then the current regularization parameter remains unchanged for the next iteration; If the deviation ratio coefficient is greater than a preset threshold, the current regularization parameter is divided by the deviation ratio coefficient to obtain the updated regularization parameter for the next iteration.

9. The method according to claim 1, characterized in that, The reconstructed three-dimensional crustal temperature structure and output temperature distribution results include: Calculate the temperature gradient between the surface and the equivalent Curie based on the surface boundary temperature and the first isothermal surface. Calculate the temperature gradient between the equivalent Curie and the Mohorovičić discontinuity based on the first and second isotherms; The temperature at any target depth is interpolated according to the piecewise temperature gradient to obtain and output a three-dimensional temperature volume and the cross-section, slice and / or isothermal surface results of the three-dimensional temperature volume.

10. The method according to claim 1, characterized in that, The method for acquiring the magnetic anomaly grid data includes: Magnetic anomaly data is acquired and then subjected to regular gridding and preprocessing to obtain magnetic anomaly grid data of the target area. The preprocessing process includes one or more of the following: outlier removal, missing value imputation, mean removal, linear trend removal, boundary extension, upward extension, pole reduction, frequency domain filtering, and / or equivalent source reconstruction.

11. A device for reconstructing crustal temperature distribution based on magnetic anomalies, characterized in that, include: The data acquisition module is used to acquire magnetic anomaly grid data of the target area; A magnetic susceptibility constraint construction module is used to construct a depth-dependent magnetic susceptibility function based on the magnetic susceptibility-temperature curve obtained from rock physics experiments and a preset temperature-depth relationship. The depth-dependent magnetic susceptibility function is used as the magnetic susceptibility structure constraint for forward and inverse magnetic anomaly modeling. A forward modeling operator for the magnetic anomaly spatial domain is established based on Cauchy surface integral. The forward modeling operator uses the depth-dependent magnetic susceptibility function to describe the magnetic susceptibility difference between the upper and lower sides of the abrupt magnetic susceptibility interface in the subsurface. The inversion module is used to input the magnetic anomaly grid data and the depth-related magnetic susceptibility function into the spatial domain forward and inversion framework based on Cauchy surface integral, construct an objective functional including data fitting terms and model stabilization terms, and iteratively solve it using the reweighted regularized conjugate gradient algorithm to obtain the equivalent Curie depth distribution. The temperature reconstruction module is used to reconstruct the three-dimensional crustal temperature structure of the region and output the temperature distribution results by taking the equivalent Curie depth distribution as the first isothermal surface, combining the surface boundary temperature and the temperature boundary with the Moho surface as the second isothermal surface. The forward modeling operator includes: calculating the X-direction component, Y-direction component, and Z-direction component of the magnetic field generated at any observation point by the interface of abrupt change in underground magnetic susceptibility; the calculation method for the X-direction component, Y-direction component, and Z-direction component of the magnetic field includes: Based on the spatial relative positional relationship between any observation point and any point on the interface of abrupt change in underground magnetic susceptibility, and combining the projection component of the uniform magnetization intensity vector in the normal direction outside the interface, the perpendicular component of the uniform magnetization intensity vector, and the non-normalized geometric factor of the normal direction outside the interface, Cauchy-type surface integral integrand terms corresponding to the X-direction component, Y-direction component, and Z-direction component of the magnetic field are constructed respectively. When calculating the magnetic field components in each direction, the coordinate difference between any observation point and any point on the interface in the corresponding direction is used as the position response term, and the interface geometry is coupled with the magnetization information to characterize the influence of the abrupt change in underground magnetic susceptibility on the magnetic field components in each direction. Two-dimensional surface integrals were performed on the underground magnetic susceptibility abrupt interface to obtain the X-direction component, Y-direction component, and Z-direction component of the magnetic field at any observation point. The method for determining the non-normalized geometric factor includes: The elevation difference between the subsurface magnetic susceptibility abrupt change interface and the reference depth is defined as an elevation function to characterize the fluctuations of the subsurface magnetic susceptibility abrupt change interface. The negative rate of change of the elevation function along the X direction is determined as the non-normalized geometric factor in the X direction; The negative rate of change of the elevation function along the Y direction is determined as the non-normalized geometric factor in the Y direction; Set the nonnormalized geometric factor in the Z direction to a unit constant.

12. An electronic device, characterized in that, include: At least one memory and at least one processor; the memory is used to store a computer program; the processor is used to execute the computer program to implement the method as claimed in any one of claims 1-10.

13. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when run on a processor, causes the processor to perform the method as described in any one of claims 1-10.

Images