Joint inversion method based on high-power artificial source electromagnetic method and micro-motion detection data

Abstract

The application relates to the technical field of geophysical exploration, in particular to a joint inversion method and device based on high-power artificial source electromagnetic method and micro-motion detection data, a storage medium and a program product. The method comprises the following steps: collecting high-power artificial source electromagnetic and micro-motion detection data of a target area; constructing an initial one-dimensional resistivity and shear wave velocity model of each observation point; constructing a horizontal constraint matrix based on all one-dimensional models to limit the physical property mutation of adjacent observation points in the same stratum; interpolating the one-dimensional model to a unified grid to obtain an initial multi-dimensional model; combining electromagnetic data, micro-motion data, horizontal constraints and cross-gradient constraints to establish a unified inversion objective function coupled with multiple constraints; iteratively updating the resistivity and wave velocity model to convergence by synchronous inversion to output an optimal joint inversion model; and delineating a geothermal anomaly area according to the model and a multi-parameter fusion identification criterion.

Description

Technical Field

[0001] This application relates to the field of geophysical exploration technology, specifically to a joint inversion method, device, storage medium, and program product based on high-power artificial source electromagnetic method and micro-motion detection data. Background Technology

[0002] With the increasing demand for urban geothermal resource exploration, geophysical exploration methods have been widely applied. Currently, high-power artificial source electromagnetic methods and micro-motion detection methods are two commonly used geophysical methods in urban geothermal exploration: one is the high-power artificial source electromagnetic method, and the other is the micro-motion detection method. However, the results of the electromagnetic method alone have poor lateral continuity, and the physical property parameters of adjacent measuring points are prone to isolated jumps; the micro-motion method alone has low accuracy in inverting deep strata and is difficult to identify fluid-bearing strata.

[0003] To overcome the limitations of single methods, researchers have proposed joint inversion methods for electromagnetic and seismic data, such as joint inversion methods for ocean-controlled source electromagnetic and magnetotelluric data based on product objective functions, and joint inversion methods for seismic travel time and magnetotelluric sounding. However, most existing joint inversion methods adopt a one-dimensional inversion strategy, that is, inverting each measuring point independently. Although the computational efficiency is high, the inversion results between adjacent measuring points are independent, resulting in discontinuous jumps in physical parameters in the lateral direction, failing to form a continuous and smooth geological profile, and seriously affecting the reliability of geological interpretation. At the same time, resistivity models and velocity models are often built on independent grids, and inconsistent depth sampling and grid mismatch make it difficult to effectively couple the two models in the same spatial coordinate system. In addition, existing methods mainly rely on empirical relationships between rock physical parameters for constraints, and the accuracy is greatly affected by the accuracy of empirical formulas, lacking a constraint mechanism based on structural consistency. To address the issue of lateral discontinuity, some methods have attempted to employ a three-dimensional inversion strategy. However, three-dimensional inversion requires solving a large number of unknown parameters simultaneously, resulting in computational complexity and high memory consumption. This makes it difficult to obtain accurate results within the engineering cycle of urban geothermal exploration, and thus fails to meet actual production needs. Summary of the Invention

[0004] The purpose of this application is to overcome the problems of horizontal discontinuity, grid mismatch, single constraint mechanism, and difficulty in practical application of the existing technology in one-dimensional inversion, and to provide a joint inversion method, device, storage medium and program product based on high-power artificial source electromagnetic method and micro-motion detection data.

[0005] To achieve the above objectives, this application provides a joint inversion method based on high-power artificial source electromagnetic method and micro-motion detection data, including: High-power artificial source electromagnetic data set and micro-motion detection data set of the target exploration area were acquired at multiple observation points, and the data were preprocessed. Each observation point corresponds to a set of artificial source electromagnetic data and a set of micro-motion detection data. Based on the preprocessed data of each observation point, an initial one-dimensional resistivity model and an initial one-dimensional shear wave velocity model corresponding to each observation point are established. A lateral constraint matrix is ​​constructed based on the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points. The lateral constraint matrix is ​​used to constrain the abrupt differences in resistivity and shear wave velocity physical property parameters of the same stratum between adjacent observation points. The initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points are interpolated onto a unified computational grid to obtain the initial multidimensional resistivity model and the initial multidimensional shear wave velocity model. Based on electromagnetic datasets, micro-motion detection datasets, lateral constraint matrices, and cross-gradient constraints, a unified inversion objective function with multiple constraints is constructed. A simultaneous inversion strategy is adopted, which updates the resistivity model and shear wave velocity model simultaneously during the iteration process until the unified inversion objective function converges and the optimal joint inversion model of underground strata is output. Based on the joint inversion model and combined with the multi-parameter fusion identification criteria, the geothermal anomaly zone of the target exploration area is delineated.

[0006] In this embodiment, a lateral constraint matrix is ​​constructed based on the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points. This includes: constructing a layer parameter lateral constraint matrix to constrain the abrupt differences in resistivity and shear wave velocity between adjacent measurement points in the same stratum; constructing a layer thickness lateral constraint matrix to constrain the abrupt differences in the depth of the interface between adjacent measurement points in the same stratum; and unifying the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model onto the same computational grid using grid interpolation technology, and establishing a cross gradient constraint term, which is used to constrain the structural similarity between the resistivity model and the shear wave velocity model.

[0007] In this embodiment of the application, the cross gradient constraint term is defined as: ,in, These are resistivity model parameters. These are the parameters for the shear wave velocity model. For gradient operators, Represents the horizontal coordinate. Represents the vertical coordinates. For the shear wave velocity model of underground media, This is a resistivity model for underground media.

[0008] In this embodiment, a unified inversion objective function with multiple constraints is constructed based on the electromagnetic data set, the micro-motion detection data set, the lateral constraint matrix, and the cross-gradient constraints. This includes: weighting the electromagnetic data fitting term and the micro-motion detection data fitting term, and combining them with the lateral constraint regularization term and the cross-gradient structure constraint term to construct the unified inversion objective function. The unified inversion objective function is expressed as follows: ,in, For electromagnetic data fitting, This is the fitting term for the micro-motion detection data. The layer parameter lateral constraint matrix, The layer thickness lateral constraint matrix is... For the joint model parameter vector, , For regularization parameters, These are the weighting coefficients of the cross-gradient constraint term. For the cross gradient constraint term, This represents the lateral constraint regularization term for layer thickness. Presentation layer parameter lateral constraint regularization term, This represents the cross-gradient structure constraint term.

[0009] In this embodiment, a simultaneous inversion strategy is adopted, updating both the resistivity model and the shear wave velocity model during the iteration process until the unified inversion objective function converges, outputting the optimal subsurface joint inversion model. This includes: using the L-BFGS algorithm for iterative calculation, simultaneously updating the resistivity model parameters and shear wave velocity model parameters in each iteration; calculating the cross gradient between the resistivity model and the shear wave velocity model after each iteration update, and applying lateral constraints to both models respectively; calculating the unified inversion objective function value and the fitting difference to determine if the convergence condition is met; if the convergence condition is met, ending the joint inversion and outputting the optimal subsurface joint inversion model; if the convergence condition is not met, adjusting the regularization parameters and weight coefficients, and iterating again.

[0010] In this embodiment of the application, geothermal anomaly zones in the target exploration area are delineated based on the joint inversion model and multi-parameter fusion identification criteria. This includes: identifying low-resistivity and low-velocity anomaly zones based on the resistivity and shear wave velocity models output by the joint inversion model; analyzing resistivity gradients to locate reservoir boundaries and identify high polarizability characteristics; analyzing velocity anisotropy to identify the main fracture development direction; and performing multi-parameter fusion of the identification results to identify spatial distribution areas with a low-resistivity-low-velocity-high polarization combination pattern, and delineating geothermal anomaly zones by combining regional geological structural characteristics.

[0011] In this embodiment of the application, an initial one-dimensional resistivity model and an initial one-dimensional shear wave velocity model corresponding to each observation point are established based on the preprocessed data of each observation point. This includes: establishing an initial layered model based on regional geological data, borehole data, or single-method inversion results; determining the initial resistivity parameters and shear wave velocity parameters of each layer; and establishing the initial coupling relationship between resistivity and shear wave velocity.

[0012] A second aspect of this application provides a joint inversion device based on high-power artificial source electromagnetic method and micro-motion detection data, the joint inversion device comprising: The data acquisition and preprocessing module is used to acquire high-power artificial source electromagnetic data sets and micro-motion detection data sets of the target exploration area at multiple observation points, and to preprocess the data. Each observation point corresponds to a set of artificial source electromagnetic data and a set of micro-motion detection data. The multidimensional model construction module is used to establish an initial one-dimensional resistivity model and an initial one-dimensional shear wave velocity model corresponding to each observation point based on the preprocessed data of each observation point; it constructs a lateral constraint matrix based on the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points, and the lateral constraint matrix is ​​used to constrain the abrupt differences in resistivity and shear wave velocity physical property parameters of the same stratum between adjacent measurement points; the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points are interpolated onto a unified computational grid to obtain the initial multidimensional resistivity model and the initial multidimensional shear wave velocity model. The objective function construction module constructs a unified inversion objective function with multiple constraints coupled based on electromagnetic datasets, micro-motion detection datasets, lateral constraint matrices, and cross-gradient constraints. The synchronous inversion iteration module is used to simultaneously update the resistivity model and the shear wave velocity model during the iteration process using a synchronous inversion strategy until the unified inversion objective function converges and outputs the optimal joint inversion model of subsurface strata. The geothermal anomaly delineation module delineates geothermal anomaly zones in the target exploration area based on the joint inversion model and multi-parameter fusion identification criteria.

[0013] In this embodiment, the multidimensional model construction module is further used to: construct a layer parameter lateral constraint matrix to constrain the resistivity and shear wave velocity abrupt differences between adjacent measuring points in the same stratum; construct a layer thickness lateral constraint matrix to constrain the abrupt differences in the interface depth of the same stratum between adjacent measuring points; and unify the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model onto the same computational grid through grid interpolation technology, and establish a cross gradient constraint term, which is used to constrain the structural similarity between the resistivity model and the shear wave velocity model.

[0014] In this embodiment of the application, the cross gradient constraint term is defined as: ,in, These are resistivity model parameters. These are the parameters for the shear wave velocity model. For gradient operators, Represents the horizontal coordinate. Represents the vertical coordinates. For the shear wave velocity model of underground media, This is a resistivity model for underground media.

[0015] In this embodiment, the objective function construction module is further configured to: weight and combine the electromagnetic method data fitting term and the micro-motion detection data fitting term, and combine them with the lateral constraint regularization term and the cross-gradient structure constraint term to construct a unified inversion objective function; the unified inversion objective function is expressed as: ,in, For electromagnetic data fitting, This is the fitting term for the micro-motion detection data. The layer parameter lateral constraint matrix, The layer thickness lateral constraint matrix is... For the joint model parameter vector, , For regularization parameters, These are the weighting coefficients of the cross-gradient constraint term. For the cross gradient constraint term, This represents the lateral constraint regularization term for layer thickness. Presentation layer parameter lateral constraint regularization term, This represents the cross-gradient structure constraint term.

[0016] In this embodiment, the synchronous inversion iteration module is further configured to: perform iterative calculations using the L-BFGS algorithm, simultaneously updating the resistivity model parameters and shear wave velocity model parameters in each iteration; after each iteration update, calculate the cross gradient between the resistivity model and the shear wave velocity model, and apply lateral constraints to the two models respectively; calculate the unified inversion objective function value and the fitting difference, and determine whether the convergence condition is met; if the convergence condition is met, terminate the joint inversion and output the optimal subsurface strata joint inversion model; if the convergence condition is not met, adjust the regularization parameters and weight coefficients, and iterate again.

[0017] In this embodiment, the geothermal anomaly delineation module is also used to: identify low-resistivity and low-velocity anomaly regions based on the resistivity model and shear wave velocity model output by the joint inversion model; analyze the resistivity gradient to locate reservoir boundaries and identify high polarizability characteristics; analyze velocity anisotropy to identify the main development direction of fractures; fuse the identification results with multiple parameters to identify spatial distribution areas with a low-resistivity-low-velocity-high polarization combination pattern, and delineate geothermal anomaly zones by combining regional geological structural characteristics.

[0018] In this embodiment of the application, the multidimensional model construction module is also used to: establish an initial layered model based on regional geological data, borehole data or single-method inversion results; determine the initial resistivity parameters and shear wave velocity parameters of each layer, and establish the initial coupling relationship between resistivity and shear wave velocity.

[0019] A third aspect of this application provides a machine-readable storage medium storing instructions that, when executed by a processor, configure the processor to perform any of the above-described joint inversion methods based on high-power artificial source electromagnetic methods and micro-motion detection data.

[0020] The fourth aspect of this application provides a computer program product, including a computer program that, when executed by a processor, implements the joint inversion method based on high-power artificial source electromagnetic method and micro-motion detection data as described above.

[0021] Based on the above technical solution, this application constructs a lateral constraint matrix to constrain the differences in resistivity, shear wave velocity, and layer thickness of the same stratum between adjacent measuring points, effectively suppressing the lateral jumps in physical parameters and interface discontinuities caused by point-by-point independent inversion in traditional one-dimensional inversion. A cross-gradient constraint function is used to force the resistivity model and shear wave velocity model to maintain spatial structural consistency, achieving structural synchronization of the two physical parameters. This constraint does not rely on empirical rock property relationships and has wider applicability, effectively reducing the problem of multiple solutions in geophysical inversion. Electromagnetic data fitting terms, micromotion data fitting terms, lateral constraint terms, cross-gradient constraint terms, and regularization terms are uniformly incorporated into the same objective function, constructing a multi-physics, multi-constraint coupled collaborative inversion framework. A synchronous inversion strategy is adopted, updating the resistivity model and shear wave velocity model simultaneously in one iteration. Compared with traditional three-dimensional joint inversion, this method achieves the lateral continuity effect of approximate two-dimensional inversion while maintaining the computational efficiency of one-dimensional inversion. Attached Figure Description

[0022] Figure 1 This illustration schematically shows a flowchart of a joint inversion method based on high-power artificial source electromagnetic method and micro-motion detection data according to an embodiment of this application; Figure 2 A schematic diagram illustrating lateral constraints according to an embodiment of this application is shown. Figure 3 A schematic diagram illustrating a grid interpolation method for joint inversion according to an embodiment of this application is shown. Figure 4 This illustration schematically shows a combined detection scheme of urban geothermal high-power field source-induced electromagnetic method and micro-motion detection method according to an embodiment of this application; Figure 5A schematic diagram of a joint inversion device based on high-power artificial source electromagnetic method and micro-motion detection data according to an embodiment of this application is shown. Figure 6 The diagram illustrates the internal structure of a computer device according to an embodiment of this application. Detailed Implementation

[0023] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are only for illustration and explanation of the embodiments of this application and are not intended to limit the embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application without creative effort are within the scope of protection of this application.

[0024] It should be noted that if the embodiments of this application involve directional indicators (such as up, down, left, right, front, back, etc.), the directional indicators are only used to explain the relative positional relationship and movement of each component in a certain specific posture (as shown in the figure). If the specific posture changes, the directional indicators will also change accordingly.

[0025] Furthermore, if the embodiments of this application involve descriptions such as "first" or "second," these descriptions are for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated. Therefore, features defined with "first" or "second" may explicitly or implicitly include at least one of those features. Additionally, the technical solutions of various embodiments can be combined with each other, but this must be based on the ability of those skilled in the art to implement them. If the combination of technical solutions is contradictory or impossible to implement, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed in this application.

[0026] In one embodiment, a schematic diagram of the joint inversion method based on high-power artificial source electromagnetic method and micro-motion detection data according to an embodiment of this application is shown. Figure 1 As shown, it includes: Step 101: Acquire high-power artificial source electromagnetic data set and micro-motion detection data set of the target exploration area at multiple observation points, and preprocess the data. Each observation point corresponds to a set of artificial source electromagnetic data and a set of micro-motion detection data.

[0027] First, multiple observation points are deployed within the target exploration area. At each point, high-power artificial source electromagnetic data and micro-motion detection data are collected. The spacing between observation points is determined based on exploration accuracy requirements and geological conditions; smaller spacing results in higher lateral resolution but also increases acquisition costs. High-power artificial source electromagnetic data is obtained by artificially transmitting high-power electromagnetic signals and receiving responses from the subsurface medium. The purpose of high-power transmission is to enhance signal strength, overcome strong electromagnetic interference from power grids and communication equipment in urban environments, and ensure data signal-to-noise ratio and detection depth. Simultaneously, the amplitude and phase of the electromagnetic field response are recorded at multiple transmission frequencies to obtain complex resistivity spectrum data, providing a data foundation for subsequent polarizability parameter extraction. Micro-motion detection data is obtained by deploying seismograph arrays on the surface, either in a ring or along two or more different orientations, continuously collecting weak ground vibration signals caused by natural factors and human activities, and extracting Rayleigh surface wave dispersion curves from these signals. The advantage of this method is that it does not require artificial seismic sources and is adaptable to construction constraints in urban environments. Each observation point corresponds to a set of artificial source electromagnetic data and a set of micro-motion detection data to ensure strict spatial matching between the two types of data, laying the foundation for subsequent joint inversion. After acquisition, both types of data undergo preprocessing: electromagnetic data preprocessing includes power frequency interference suppression, static correction, terrain correction, and denoising; micro-motion data preprocessing includes signal enhancement, dispersion curve extraction, and quality control. The purpose of preprocessing is to remove noise and interference from the raw data, improve the signal-to-noise ratio, and provide high-quality data input for subsequent inversion.

[0028] Step 102: Based on the preprocessed data of each observation point, establish the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model corresponding to each observation point.

[0029] The quality of the initial model directly affects the convergence speed and accuracy of the final result of the inversion. Therefore, establishing an initial one-dimensional layered model that closely approximates the actual underground conditions can reduce the ambiguity of the inversion and improve convergence efficiency.

[0030] The processor establishes an initial one-dimensional resistivity model and an initial one-dimensional shear wave velocity model corresponding to each observation point based on the preprocessed data of each observation point. This includes: establishing an initial layered model based on regional geological data, borehole data, or single-method inversion results; determining the initial resistivity parameters and shear wave velocity parameters of each layer; and establishing the initial coupling relationship between resistivity and shear wave velocity.

[0031] Specifically, based on existing regional geological survey results of the target exploration area, the macroscopic geological characteristics such as the regional stratigraphic sequence, lithological assemblage, and structural framework are understood, thereby determining the stratigraphic sequence and approximate depth range. If there are existing boreholes such as geothermal wells, hydrological wells, or oil exploration wells in the study area or adjacent areas, the stratigraphic layering information (depth, thickness, and lithology of each layer) revealed by the boreholes can be directly used to establish an accurate initial stratigraphic model. In the embodiments of this application, drilling verification results show that the depth error of the stratigraphic interface obtained by inversion is less than 100m, indicating that borehole data has an important constraining effect on the construction of the initial model.

[0032] After establishing the initial stratigraphic model, initial physical properties need to be assigned to each stratum. Based on the lithological characteristics of each stratum in the regional geological data, such as Quaternary sediments, carbonate rocks, and granites, and referring to the known empirical range of rock resistivity values, an initial resistivity is assigned to each stratum. If borehole logging data, such as resistivity logging curves, is available, the logging values ​​can be directly used as the initial resistivity. Similarly, based on the lithology of the strata, referring to the empirical range of rock wave velocities, or using dispersion curves extracted from microseismic data for initial inversion, the initial shear wave velocity of each layer is obtained.

[0033] Since the electromagnetic method is sensitive to resistivity and the micromotion method is sensitive to shear wave velocity, the two parameters belong to different physical fields. In order to achieve joint inversion, it is necessary to establish a coupling relationship between the two in the initial model so that the two parameters can be mutually constrained and updated collaboratively during the inversion process.

[0034] In this application embodiment, the initial coupling relationship is established in the following ways: Empirical constraints are used to establish a preliminary conversion relationship between resistivity and shear wave velocity based on rock physics theories or empirical formulas for the study area. For example, for clastic rock formations, resistivity and wave velocity are often positively correlated; for carbonate rock formations, the relationship is related to porosity and fluid properties.

[0035] The structural consistency assumption assumes that the resistivity interface and the wave velocity interface are basically in the same position in the initial model, that is, the electrical boundary and elastic boundary of the underground geological body coincide, providing initial conditions for the subsequent cross-gradient structural constraints.

[0036] In the absence of boreholes or insufficient geological data, the results of single-method inversion can be used to perform single-method inversion on electromagnetic data and micromotion data separately, such as traditional one-dimensional inversion, and the results can be used as the initial model for joint inversion.

[0037] Step 103: Construct a transverse constraint matrix based on the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points. The transverse constraint matrix is ​​used to constrain the abrupt differences in resistivity and shear wave velocity physical property parameters of the same stratum between adjacent measurement points.

[0038] Step 104: Interpolate the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points onto a unified computational grid to obtain the initial multidimensional resistivity model and the initial multidimensional shear wave velocity model.

[0039] Based on the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points, a lateral constraint matrix is ​​constructed, including: constructing a layer parameter lateral constraint matrix to constrain the abrupt differences in resistivity and shear wave velocity between adjacent measurement points in the same stratum; constructing a layer thickness lateral constraint matrix to constrain the abrupt differences in the depth of the interface between adjacent measurement points in the same stratum; and unifying the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model onto the same computational grid through grid interpolation technology, and establishing a cross gradient constraint term, which is used to constrain the structural similarity between the resistivity model and the shear wave velocity model.

[0040] Specifically, the construction of the lateral constraint matrix includes three aspects. First, a lateral constraint matrix for layer parameters is constructed. This matrix is ​​used to constrain abrupt changes in resistivity and shear wave velocity between adjacent measuring points within the same stratum. For resistivity parameters, the resistivity difference between adjacent measuring points within the same stratum is constrained; for shear wave velocity parameters, the wave velocity difference between adjacent measuring points within the same stratum is constrained. Through this constraint, the physical property parameters change smoothly laterally, avoiding unreasonable abrupt changes.

[0041] Secondly, a lateral constraint matrix for layer thickness is constructed. This matrix is ​​used to constrain abrupt differences in the depth of the same stratum interface between adjacent measuring points, ensuring the lateral continuity of the stratum interface and avoiding unreasonable depth jumps in the same stratum between adjacent measuring points. Figure 2 The mechanism of the above constraints is illustrated. The resistivity ρ and layer thickness h between adjacent measuring points are correlated through a lateral constraint matrix, forcing the parameter differences between adjacent measuring points to be minimized, thereby achieving lateral continuity. It should be noted that, in the embodiments of this application, adjacent measuring points refer to adjacent measuring points on a measuring line.

[0042] Thirdly, to address the issue of poor coupling caused by the mismatch between the resistivity and velocity models' depth grids, this application unifies the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model onto the same computational grid using grid interpolation technology. A unified two-dimensional computational grid is established based on the exploration target and depth range. It simultaneously meets the forward modeling requirements of both electromagnetic and micro-motion detection methods; the one-dimensional resistivity model of each observation point is interpolated to the corresponding nodes of a unified grid to obtain the initial multi-dimensional resistivity model ρ. The one-dimensional shear wave velocity model at each observation point is interpolated to the corresponding nodes of a unified grid to obtain the initial multi-dimensional shear wave velocity model. After interpolation, both models have the same computational grid, meaning that for any grid node... Simultaneously, there exists a resistivity value ρ and wave velocity value ,like Figure 3 As shown, this illustrates how the resistivity ρ, wave velocity v, and layer thickness h at each measurement point (1, 2, ..., b) are interpolated from a one-dimensional model to a two-dimensional unified mesh.

[0043] Furthermore, a cross-gradient constraint term is established to constrain the structural similarity between the resistivity model and the shear wave velocity model. The cross-gradient constraint term is defined as: ,in, These are resistivity model parameters. These are the parameters for the shear wave velocity model. For gradient operators, Represents the horizontal coordinate. Represents the vertical coordinates. For the shear wave velocity model of underground media, This is a resistivity model for the subsurface medium. This constraint ensures that the resistivity model and velocity model have similar structural characteristics, are independent of rock physical properties, and emphasize the similarity of the subsurface structure, thus having wider applicability. Furthermore, the gradient operator is expressed in three-dimensional Cartesian coordinates as follows: + + , The x-direction is the unit vector along the survey line. It is a unit vector in the y-direction, perpendicular to the survey line direction; Let z be the unit vector in the z-direction and depth direction. From this, the gradient of the resistivity model can be derived: + + ,make This represents the rate of change of resistivity in the x-direction. This represents the rate of change of resistivity in the y-direction. This represents the rate of change of resistivity in the z-direction. Therefore: + + Similarly, the gradient of the shear wave velocity model is: + + The cross product of the two vectors is calculated using the determinant, and after expansion, we get... + + , The cross gradient components in the x-direction are... The cross gradient components in the y-direction are... This represents the cross gradient component in the z-direction. In the two-dimensional inversion, it is assumed that the model does not change in the y-direction perpendicular to the survey line, therefore... =0, =0, substitute into the expressions of each component, and finally get = =( - ) .

[0044] Step 105: Based on the electromagnetic dataset, the micro-motion detection dataset, the lateral constraint matrix, and the cross-gradient constraint, construct a unified inversion objective function with multiple constraints coupled.

[0045] Based on electromagnetic datasets, micro-motion detection datasets, lateral constraint matrices, and cross-gradient constraints, a unified inversion objective function with multiple constraints is constructed. This includes: weighting and combining the fitting terms for electromagnetic data and micro-motion detection data, and combining them with lateral constraint regularization and cross-gradient structure constraint terms to construct the unified inversion objective function. The unified inversion objective function is expressed as follows: ,in, For electromagnetic data fitting, This is the fitting term for the micro-motion detection data. The layer parameter lateral constraint matrix, The layer thickness lateral constraint matrix is... For the joint model parameter vector, , For regularization parameters, These are the weighting coefficients of the cross-gradient constraint term. For the cross gradient constraint term, This represents the lateral constraint regularization term for layer thickness. Presentation layer parameter lateral constraint regularization term, This represents the cross-gradient structure constraint term. To standardize shear wave velocity, This is for standardized resistivity.

[0046] In this embodiment, the electromagnetic data fitting term is used to measure the degree of difference between the electromagnetic forward modeling response and the measured electromagnetic data. The input data includes the apparent resistivity, phase, and impedance tensor of the controlled source electromagnetic data, and the model parameters are a resistivity model. A smaller fitting term indicates that the current resistivity model can better explain the measured electromagnetic data. The electromagnetic data fitting term can be further expressed as: ,in, Represents the data weight matrix. The residual vector between the observed data and the forward response data. The layer parameter lateral constraint matrix, The layer thickness lateral constraint matrix is... , For regularization parameters, These are the weighting coefficients of the cross-gradient constraint term in the electromagnetic data fitting term. This is the cross gradient constraint term.

[0047] The fitting term for the micromotion detection data is used to measure the degree of difference between the forward modeling response of the micromotion method and the measured micromotion data. Input data includes the surface wave dispersion curve, shear wave velocity, and HVSR spectral ratio of the micromotion data, while the model parameters are the shear wave velocity model. A smaller fitting term indicates that the current shear wave velocity model can better explain the measured micromotion data. The fitting term for the micromotion detection data can be further expressed as: , These are the weighting coefficients of the cross-gradient constraint terms in the micromotion detection data fitting terms. The lateral constraint regularization terms include layer thickness lateral constraint regularization and layer parameter lateral constraint regularization. The layer thickness lateral constraint regularization term constrains the depth difference of the same stratigraphic interface between adjacent measuring points through the layer thickness lateral constraint matrix, ensuring a smooth lateral change in the stratigraphic interface and avoiding unreasonable depth jumps between adjacent measuring points. The layer parameter lateral constraint regularization term constrains the resistivity and shear wave velocity differences of the same stratigraphic layer between adjacent measuring points through the layer parameter lateral constraint matrix, ensuring continuous lateral change in physical properties and overcoming the problem of discontinuity in parameters between adjacent measuring points in traditional one-dimensional inversion. Regularization parameters are used to control the weights of their respective constraint terms, and their values ​​are dynamically adjusted based on the data signal-to-noise ratio, detection capability, and geological conditions. The cross-gradient structure constraint term constrains the structural similarity between the resistivity model and the shear wave velocity model. By minimizing the cross-gradient structure constraint term, the resistivity model and the velocity model are forced to have consistent structural boundaries, requiring that the locations of resistivity changes and wave velocity changes essentially coincide. This constraint does not rely on empirical relationships of rock physical properties, but emphasizes the similarity of underground structures, thus having wider applicability. Weighting coefficients are used to control the strength of the cross-gradient constraint and can be adjusted based on convergence during the inversion iteration process.

[0048] The four components mentioned above are unified into a single objective function through weighted summation. Specifically, the electromagnetic data fitting term and the micro-motion detection data fitting term ensure the accuracy of the inversion model's fit to the observed data; the lateral constraint regularization term ensures the continuity of the inversion results in the lateral direction; and the cross-gradient structure constraint term ensures the structural consistency between the resistivity model and the velocity model. These four components mutually constrain and synergistically optimize each other, together forming a unified inversion objective function with multi-constraint coupling.

[0049] Step 106: Using a synchronous inversion strategy, the resistivity model and shear wave velocity model are updated simultaneously during the iteration process until the unified inversion objective function converges, and the optimal subsurface strata joint inversion model is output.

[0050] The processor employs a synchronous inversion strategy, simultaneously updating the resistivity model and shear wave velocity model during the iteration process until the unified inversion objective function converges, outputting the optimal joint inversion model of subsurface strata. This includes: using the L-BFGS algorithm for iterative calculation, simultaneously updating the resistivity model parameters and shear wave velocity model parameters in each iteration; calculating the cross gradient between the resistivity model and the shear wave velocity model after each iteration update, and applying lateral constraints to both models respectively; calculating the unified inversion objective function value and the fitting difference to determine if the convergence condition is met; if the convergence condition is met, ending the joint inversion and outputting the optimal joint inversion model of subsurface strata; if the convergence condition is not met, adjusting the regularization parameters and weight coefficients, and iterating again.

[0051] Specifically, the L-BFGS algorithm is used to synchronously update the model parameters. The L-BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shanno) algorithm is a quasi-Newton optimization algorithm suitable for large-scale inversion problems. In each iteration, the processor updates both the resistivity model parameters and the shear wave velocity model parameters simultaneously, rather than alternating them. The resistivity and shear wave velocity parameters change collaboratively in the same iteration step, ensuring that the coupling relationship between them is reflected in the model update in real time, avoiding the parameter asynchrony problems that might occur with alternating updates. After each iteration, the processor executes the following two sub-steps: First, it calculates the cross gradient between the resistivity model and the shear wave velocity model. The magnitude of the cross gradient reflects the directional consistency between the resistivity model gradient and the wave velocity model gradient. When the two gradient directions are parallel or antiparallel, the cross gradient is zero, indicating that the structural boundaries of the two models are completely consistent; when the two gradient directions are inconsistent, the cross gradient is non-zero, requiring penalty and correction in subsequent iterations. Second, lateral constraints are applied to the resistivity model and the shear wave velocity model respectively, including lateral constraints on layer parameters, which constrain the resistivity and wave velocity differences of the same stratum at adjacent measuring points, and lateral constraints on layer thickness, which constrain the differences in the depth of the interface of the same stratum at adjacent measuring points.

[0052] After each iteration, the processor calculates the unified inversion objective function value and the fitting difference. Convergence criteria are determined in two ways: objective function convergence (the change in the objective function value over multiple iterations is less than a preset threshold), and fitting difference (the inversion result is considered satisfactory when the fitting difference between the observed data and the forward response data reaches the expected accuracy). If convergence criteria are not met, the processor adjusts the regularization parameters and weight coefficients based on the current convergence status and then iterates again. Specifically, the regularization parameters are adjusted based on the rate of decrease of the objective function; if convergence is too slow, the regularization parameters are appropriately reduced to enhance the weight of the data fitting terms. The weight coefficients are adjusted based on the changing trend of the cross-gradient terms; if the structural differences between the resistivity model and the wave velocity model are significant, they are increased to strengthen structural constraints. The data weight matrix is ​​adjusted based on the changes in the fitting difference; if the fitting difference for a certain frequency band or measurement point remains large, the weight of that data segment is increased.

[0053] In this embodiment, the L-BFGS algorithm is implemented as follows: the processor stores gradient information and parameter update amounts for a finite number of steps, and uses this historical information to construct an approximate inverse matrix of the Hessian matrix, thereby determining the current update direction. For the joint inversion problem, the gradients of the resistivity model parameters and the shear wave velocity model parameters of the objective function are calculated separately, and then combined into a joint gradient vector. The L-BFGS algorithm determines the update step size and update direction of both models based on this joint gradient vector. This processing method significantly reduces memory usage while ensuring convergence speed, and is suitable for large-scale two-dimensional and three-dimensional inversion calculations. When the convergence condition is met, the processor outputs the optimal subsurface joint inversion model. This joint inversion model includes the optimal resistivity model and the optimal shear wave velocity model. The two models have a unified computational grid and consistent structural boundary positions. This joint inversion model can be directly used for subsequent geothermal anomaly zone delineation, reservoir parameter evaluation, and well location deployment decisions.

[0054] Step 107: Based on the joint inversion model and combined with the multi-parameter fusion identification criteria, delineate the geothermal anomaly zone of the target exploration area.

[0055] Based on the joint inversion model and combined with multi-parameter fusion identification criteria, geothermal anomaly zones in the target exploration area are delineated. This includes: the processor identifying low-resistivity and low-velocity anomaly zones based on the resistivity and shear wave velocity models output by the joint inversion model; analyzing resistivity gradients to locate reservoir boundaries and identify high polarizability characteristics; analyzing velocity anisotropy to identify the main fracture development direction; and fusing the identification results with multiple parameters to identify spatial distribution areas with a low-resistivity-low-velocity-high polarization combination pattern, and combining this with regional geological structural characteristics for constraint to delineate geothermal anomaly zones.

[0056] Based on the resistivity model output by the joint inversion model, low-resistivity anomaly identification is performed: hydrothermal geothermal reservoirs contain high-temperature groundwater, which exhibits low resistivity characteristics due to enhanced ionic conductivity. By analyzing the resistivity model obtained from the inversion, regions with resistivity values ​​lower than the background values ​​of the surrounding rock are identified and marked as potential low-resistivity anomaly zones. Resistivity gradient analysis: There is a significant resistivity difference between the geothermal reservoir and the surrounding rock. By calculating the spatial gradient change of resistivity, the reservoir boundary can be accurately located. Regions with large resistivity gradient values ​​correspond to interfaces with drastic changes in physical properties, usually the top or bottom interface of the reservoir. High polarizability feature identification: High-temperature groundwater in hydrothermal geothermal reservoirs exhibits high polarizability characteristics compared to the surrounding rock due to enhanced ionic conductivity. By analyzing the spatial distribution of polarizability parameters, the location of the geothermal reservoir is further identified.

[0057] Based on the shear wave velocity model output by the joint inversion model, low-velocity anomalies are identified: high-temperature fluids in geothermal reservoirs reduce the elastic modulus of rocks, leading to a decrease in shear wave velocity. By analyzing the shear wave velocity model, regions with wave velocity values ​​lower than the background values ​​of the surrounding rock are identified. Velocity anisotropy analysis: fracture systems in geothermal reservoirs often exhibit directionality, resulting in anisotropic velocity characteristics. By analyzing the variation patterns of shear wave velocity in different directions, the main development directions of fractures can be identified.

[0058] In this embodiment, the polarizability parameter is obtained by inverting the complex resistivity spectrum data collected by the electromagnetic method of a high-power artificial source, and the velocity anisotropy is obtained by azimuth dispersion analysis of arrays in different orientations.

[0059] Specifically, using the amplitude and phase data of electromagnetic field responses at multiple transmission frequencies collected in step 101, a complex resistivity spectrum is constructed. Based on the complex resistivity spectrum model, the measured spectrum is fitted and inverted to extract spectral parameters such as polarizability (charge rate), frequency correlation coefficient, and time constant. The inverted polarizability parameters are then interpolated onto a computational grid unified with the joint inversion model to obtain a polarizability model, which is used for high polarizability feature identification. Using the micromotion signals collected by the ring array or arrays in different azimuths deployed in step 101, Rayleigh surface wave dispersion curves in different azimuths are extracted and inverted to obtain the corresponding shear wave velocity models. By comparing the differences between the shear wave velocity models in different azimuths, the velocity anisotropy parameters are calculated to identify the main development direction of fractures and provide fracture direction constraints for the delineation of geothermal anomaly zones.

[0060] The processor integrates the electrical characteristic identification results with the velocity characteristic identification results to identify spatial distribution areas with a "low resistivity-low velocity-high polarization" combined pattern. The physical meaning of this combined pattern is that the geothermal reservoir simultaneously meets three conditions: low resistivity (containing fluids), low shear wave velocity (softened rock skeleton), and high polarizability (high mineralization in water, ionic conductivity in water). These three conditions corroborate each other, significantly improving the reliability of the identification.

[0061] Based on the identification of the "low resistance-low speed-high polarization" combination mode, the processor further conducts constraint analysis in conjunction with regional geological structural characteristics: it analyzes the relationship between geothermal reservoirs and fault structures, as geothermal resources are usually distributed along fault zones, which provide migration channels and storage space for geothermal water; it analyzes the reservoir characteristics and caprock features of geothermal reservoirs, as in sedimentary basins, the cuttable parts of basement or caprock faults are often favorable development sites for geothermal reservoirs; and it prioritizes the identified geothermal anomaly areas in conjunction with regional geothermal background such as geothermal gradient and heat flow distribution.

[0062] The identification results from the four aspects mentioned above complement and verify each other. Electrical characteristic identification is sensitive to fluids but not to interface location, while velocity characteristic identification is sensitive to interface location but not to fluid properties. Combining the two can simultaneously obtain accurate reservoir location and fluid property information. The triple constraint of the "low resistivity-low velocity-high polarization" combined mode significantly reduces the ambiguity of single parameter identification and improves the accuracy and reliability of geothermal anomaly delineation.

[0063] In one embodiment, the processor first performs threshold segmentation on the resistivity model, extracting regions with resistivity below a certain range from the background value to obtain low-resistivity anomaly zones, such as regions below 30% of the background value. It then performs threshold segmentation on the shear wave velocity model, extracting regions with wave velocities below a certain range from the background value to obtain low-velocity anomaly zones, such as regions below 20% of the background value. Finally, it performs threshold segmentation on the polarizability model, extracting regions with polarizability above a certain range from the background value, such as regions above 50% of the background value, to obtain high-polarizability anomaly zones. The intersection of these three models is then calculated to obtain geothermal anomaly candidate zones that simultaneously meet all three conditions. Finally, combined with regional geological structure maps, the processor determines whether the candidate zones are located near fault zones or favorable structural locations, filters and prioritizes the candidate zones, and ultimately delineates the geothermal anomaly zones with the greatest exploration potential. After delineating the geothermal anomaly zones, the processor can output a spatial distribution map of the geothermal anomaly zones, including anomaly zone boundary coordinates, burial depth range, anomaly intensity, and other information, to guide subsequent risk well placement.

[0064] Through the above technical solution, this application combines high-power artificial source electromagnetic method with micro-motion detection data for joint inversion, fully leveraging the complementary advantages of electromagnetic method's sensitivity to low-resistivity fluids and micro-motion method's high resolution of formation interfaces, effectively overcoming the limitations of single methods in urban environments with strong interference. By constructing a lateral constraint matrix to constrain the differences in physical parameters between adjacent measuring points and employing grid interpolation technology to unify resistivity and velocity models to the same computational grid, the coupling problems of poor lateral continuity and grid mismatch in traditional methods are solved. Based on a unified inversion objective function and synchronous inversion strategy, true collaborative inversion of electromagnetic and seismic data is achieved, significantly improving the identification accuracy of deep geothermal reservoirs. Combining the "low-resistivity-low-velocity-high-polarization" multi-parameter fusion identification criterion to delineate geothermal anomaly zones provides an efficient and reliable solution for the exploration of deep hydrothermal geothermal resources in urban areas.

[0065] Figure 4 This illustration schematically depicts a combined detection scheme of urban geothermal high-power field source-induced electromagnetic method and micro-motion detection method according to an embodiment of this application. Figure 4 The diagram illustrates a comprehensive detection scheme for geothermal resource exploration in an urban environment, as presented in this application. Two acquisition systems are simultaneously deployed: a high-power artificial source electromagnetic method (CSAMT / time-frequency electromagnetic method) and a micro-motion detection method (SPAC array). The electromagnetic method system includes a high-power transmitter Tx positioned 1-3 km outside the city at a distance from the survey area, and electric field measuring electrodes Ex and magnetic probe Hy ​​deployed along the survey line within the survey area. These systems are used to receive artificial source electromagnetic field signals and calculate apparent resistivity and phase. The micro-motion detection system includes a seismograph and a SPAC array, used to continuously acquire natural source surface wave signals and extract Rayleigh surface wave dispersion curves. The diagram also marks the CSAMT and MS labels above the survey area, indicating that the two methods are used for synchronous acquisition within the same survey area and along the same survey line, ensuring spatial consistency of the data. The stratigraphic sequence of the strata below the survey area is shown from the shallow Cenozoic Quaternary to the mid-deep Proterozoic strata. The geothermal reservoir (Q) is marked in an ellipse in the mid-deep carbonate strata, exhibiting low resistivity, low velocity, and high polarization characteristics.

[0066] In one embodiment, a joint inversion device 500 based on high-power artificial source electromagnetic method and micro-motion detection data is provided, such as... Figure 5 As shown, it includes: The data acquisition and preprocessing module 501 is used to acquire high-power artificial source electromagnetic data sets and micro-motion detection data sets of the target exploration area at multiple observation points, and to preprocess the data. Each observation point corresponds to a set of artificial source electromagnetic data and a set of micro-motion detection data.

[0067] The multidimensional model construction module 502 is used to establish an initial one-dimensional resistivity model and an initial one-dimensional shear wave velocity model corresponding to each observation point based on the preprocessed data of each observation point; to construct a lateral constraint matrix based on the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points, the lateral constraint matrix is ​​used to constrain the abrupt differences in resistivity and shear wave velocity physical property parameters of the same stratum between adjacent measurement points; and to interpolate the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points onto a unified computational grid to obtain the initial multidimensional resistivity model and the initial multidimensional shear wave velocity model.

[0068] The objective function construction module 503 constructs a unified inversion objective function with multiple constraints coupled based on the electromagnetic dataset, the micro-motion detection dataset, the lateral constraint matrix, and the cross-gradient constraints.

[0069] The synchronous inversion iteration module 504 is used to simultaneously update the resistivity model and the shear wave velocity model during the iteration process using a synchronous inversion strategy until the unified inversion objective function converges and outputs the optimal joint inversion model of underground strata.

[0070] The geothermal anomaly delineation module 505 delineates geothermal anomaly zones in the target exploration area based on the joint inversion model and multi-parameter fusion identification criteria.

[0071] In one embodiment, the multidimensional model building module 502 is further used to: establish an initial layered model based on regional geological data, borehole data, or single-method inversion results; determine the initial resistivity parameters and shear wave velocity parameters of each layer, and establish the initial coupling relationship between resistivity and shear wave velocity.

[0072] In one embodiment, the multidimensional model construction module 502 is further configured to: construct a layer parameter lateral constraint matrix to constrain the resistivity and shear wave velocity abrupt differences between adjacent measuring points in the same stratum; construct a layer thickness lateral constraint matrix to constrain the abrupt differences in the interface depth of the same stratum between adjacent measuring points; and unify the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model onto the same computational grid using grid interpolation technology, and establish a cross gradient constraint term, which is used to constrain the structural similarity between the resistivity model and the shear wave velocity model.

[0073] The cross gradient constraint term is defined as: ,in, These are resistivity model parameters. These are the parameters for the shear wave velocity model. For gradient operators, Represents the horizontal coordinate. Represents the coordinates in the vertical direction.

[0074] In one embodiment, the objective function construction module 503 is further configured to: weight and combine the electromagnetic method data fitting term and the micro-motion detection data fitting term, and combine the lateral constraint regularization term and the cross-gradient structure constraint term to construct a unified inversion objective function; the unified inversion objective function is expressed as: ,in, For electromagnetic data fitting, This is the fitting term for the micro-motion detection data. The layer parameter lateral constraint matrix, The layer thickness lateral constraint matrix is... For the joint model parameter vector, , For regularization parameters, These are the weighting coefficients of the cross-gradient constraint term. For the cross gradient constraint term, This represents the lateral constraint regularization term for layer thickness. Presentation layer parameter lateral constraint regularization term, This represents the cross-gradient structure constraint term.

[0075] In one embodiment, the synchronous inversion iteration module 504 is further configured to: perform iterative calculations using the L-BFGS algorithm, simultaneously updating the resistivity model parameters and shear wave velocity model parameters in each iteration; after each iteration update, calculate the cross gradient between the resistivity model and the shear wave velocity model, and apply lateral constraints to the two models respectively; calculate the unified inversion objective function value and the fitting difference, and determine whether the convergence condition is met; if the convergence condition is met, terminate the joint inversion and output the optimal subsurface strata joint inversion model; if the convergence condition is not met, adjust the regularization parameters and weight coefficients, and iterate again.

[0076] In one embodiment, the geothermal anomaly delineation module 505 is further configured to: identify low-resistivity and low-velocity anomaly regions based on the resistivity model and shear wave velocity model output by the joint inversion model; analyze the resistivity gradient to locate reservoir boundaries and identify high polarizability characteristics; analyze velocity anisotropy to identify the main development direction of fractures; fuse the identification results with multiple parameters to identify spatial distribution areas with a low-resistivity-low-velocity-high polarization combination pattern, and delineate geothermal anomaly zones by combining regional geological structural characteristics.

[0077] In one embodiment, a machine-readable storage medium is provided, on which instructions are stored, which, when executed by a processor, cause the processor to be configured to perform any of the above-described joint inversion methods based on high-power artificial source electromagnetic methods and micro-motion detection data.

[0078] In one embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the joint inversion method based on high-power artificial source electromagnetic method and micro-motion detection data as described above.

[0079] In one embodiment, a computer device is provided, which may be a terminal, and its internal structure diagram may be as follows: Figure 6 As shown in the figure, the computer device includes a processor A01, a network interface A02, a display screen A04, an input device A05, and a memory (not shown) connected via a system bus. The processor A01 provides computing and control capabilities. The memory includes internal memory A03 and a non-volatile storage medium A06. The non-volatile storage medium A06 stores an operating system B01 and a computer program B02. The internal memory A03 provides an environment for the operation of the operating system B01 and the computer program B02 stored in the non-volatile storage medium A06. The network interface A02 is used for communication with external terminals via a network connection. When the computer program is executed by the processor A01, it implements a joint inversion method based on high-power artificial source electromagnetic methods and micro-motion detection data. The display screen A04 can be a liquid crystal display or an e-ink display. The input device A05 can be a touch layer covering the display screen, buttons, a trackball, or a touchpad mounted on the computer device casing, or an external keyboard, touchpad, or mouse.

[0080] Those skilled in the art will understand that Figure 6 The structure shown is merely a block diagram of a portion of the structure related to the present application and does not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0081] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0082] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0083] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0084] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0085] In a typical configuration, a computing device includes one or more processors (CPU), input / output interfaces, network interfaces, and memory.

[0086] Memory may include non-persistent memory in computer-readable media, such as random access memory (RAM) and / or non-volatile memory, such as read-only memory (ROM) or flash RAM. Memory is an example of computer-readable media.

[0087] Computer-readable media includes both permanent and non-permanent, removable and non-removable media that can store information using any method or technology. Information can be computer-readable instructions, data structures, modules of programs, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other non-transferable medium that can be used to store information accessible by a computing device. As defined herein, computer-readable media does not include transient computer-readable media, such as modulated data signals and carrier waves.

[0088] It should also be noted that the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0089] The above are merely embodiments of this application and are not intended to limit the scope of this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the scope of the claims of this application.

Claims

1. A joint inversion method based on high-power artificial source electromagnetic method and micro-motion detection data, characterized in that, The joint inversion method includes: High-power artificial source electromagnetic data set and micro-motion detection data set of the target exploration area were acquired at multiple observation points, and the data were preprocessed. Each observation point corresponds to a set of artificial source electromagnetic data and a set of micro-motion detection data. Based on the preprocessed data of each observation point, an initial one-dimensional resistivity model and an initial one-dimensional shear wave velocity model corresponding to each observation point are established. A lateral constraint matrix is ​​constructed based on the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points. The lateral constraint matrix is ​​used to constrain the abrupt differences in resistivity and shear wave velocity physical property parameters of the same stratum between adjacent measurement points. The initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points are interpolated onto a unified computational grid to obtain the initial multidimensional resistivity model and the initial multidimensional shear wave velocity model. Based on the electromagnetic dataset, the micro-motion detection dataset, the lateral constraint matrix, and the cross-gradient constraints, a unified inversion objective function with multiple constraints is constructed. A synchronous inversion strategy is adopted, in which the resistivity model and the shear wave velocity model are updated simultaneously during the iteration process until the unified inversion objective function converges, and the optimal subsurface strata joint inversion model is output. Based on the joint inversion model and combined with the multi-parameter fusion identification criteria, the geothermal anomaly zone of the target exploration area is delineated.

2. The joint inversion method according to claim 1, characterized in that, The lateral constraint matrix is ​​constructed based on the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points, including: A lateral constraint matrix for layer parameters is constructed to constrain the abrupt differences in resistivity and shear wave velocity between adjacent measuring points in the same stratum. Construct a lateral constraint matrix for layer thickness to constrain abrupt differences in the depth of the same stratum interface between adjacent measuring points; Furthermore, the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model are unified onto the same computational grid using grid interpolation technology, and a cross gradient constraint term is established. The cross gradient constraint term is used to constrain the structural similarity between the resistivity model and the shear wave velocity model.

3. The joint inversion method according to claim 2, characterized in that, The cross gradient constraint term is defined as follows: in, These are resistivity model parameters. These are the parameters for the shear wave velocity model. For gradient operators, Represents the horizontal coordinate. Represents the vertical coordinates. For the shear wave velocity model of underground media, This is a resistivity model for underground media.

4. The joint inversion method according to claim 1, characterized in that, The construction of a unified inversion objective function based on the electromagnetic dataset, the micro-motion detection dataset, the lateral constraint matrix, and the cross-gradient constraints includes: The electromagnetic method data fitting term and the micro-motion detection data fitting term are weighted and combined, and combined with the lateral constraint regularization term and the cross gradient structure constraint term to construct a unified inversion objective function; The unified inversion objective function is expressed as follows: in, For electromagnetic data fitting, This is the fitting term for the micro-motion detection data. The layer parameter lateral constraint matrix, The layer thickness lateral constraint matrix is... For the joint model parameter vector, , For regularization parameters, These are the weighting coefficients of the cross-gradient constraint term. For the cross gradient constraint term, This represents the lateral constraint regularization term for layer thickness. Presentation layer parameter lateral constraint regularization term, This represents the cross-gradient structure constraint term.

5. The joint inversion method according to claim 1, characterized in that, The simultaneous inversion strategy updates the resistivity model and shear wave velocity model simultaneously during the iteration process until the unified inversion objective function converges, outputting the optimal joint inversion model of subsurface strata, including: The L-BFGS algorithm is used for iterative calculation, and the resistivity model parameters and shear wave velocity model parameters are updated simultaneously in each iteration. After each iteration update, the cross gradient between the resistivity model and the shear wave velocity model is calculated, and lateral constraints are applied to the two models respectively. Calculate the unified inversion objective function value and the fitting difference to determine whether the convergence condition is met; If the convergence condition is met, the joint inversion ends and the optimal subsurface strata joint inversion model is output. If the convergence condition is not met, adjust the regularization parameters and weight coefficients, and iterate again.

6. The joint inversion method according to claim 1, characterized in that, The process of delineating the geothermal anomaly zone of the target exploration area based on the joint inversion model and multi-parameter fusion identification criteria includes: Based on the resistivity model and shear wave velocity model output by the joint inversion model, low-resistivity anomaly regions and low-velocity anomaly regions are identified. Analyze resistivity gradients to locate reservoir boundaries and identify high polarizability features; Analyze velocity anisotropy to identify the main direction of fracture development; The identification results are fused with multiple parameters to identify spatial distribution areas with a combination of low resistance, low speed and high polarization. The geothermal anomaly zone is delineated by combining regional geological structural characteristics with constraints.

7. The joint inversion method according to claim 1, characterized in that, The process of establishing an initial one-dimensional resistivity model and an initial one-dimensional shear wave velocity model corresponding to each observation point based on the preprocessed data of each observation point includes: An initial layered model is established based on regional geological data, borehole data, or inversion results from a single method. Determine the initial resistivity and shear wave velocity parameters for each layer, and establish the initial coupling relationship between resistivity and shear wave velocity.

8. A joint inversion device based on high-power artificial source electromagnetic method and micro-motion detection data, characterized in that, The joint inversion device includes: The data acquisition and preprocessing module is used to acquire high-power artificial source electromagnetic data sets and micro-motion detection data sets of the target exploration area at multiple observation points, and to preprocess the data. Each observation point corresponds to a set of artificial source electromagnetic data and a set of micro-motion detection data. The multidimensional model construction module is used to establish an initial one-dimensional resistivity model and an initial one-dimensional shear wave velocity model corresponding to each observation point based on the preprocessed data of each observation point; it constructs a lateral constraint matrix based on the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points, which is used to constrain the abrupt differences in resistivity and shear wave velocity physical property parameters of the same stratum between adjacent measurement points; and it interpolates the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model of all observation points onto a unified computational grid to obtain the initial multidimensional resistivity model and the initial multidimensional shear wave velocity model. The objective function construction module constructs a unified inversion objective function with multiple constraints coupled based on the electromagnetic dataset, the micro-motion detection dataset, the lateral constraint matrix, and the cross-gradient constraints. The synchronous inversion iteration module is used to simultaneously update the resistivity model and the shear wave velocity model during the iteration process using a synchronous inversion strategy until the unified inversion objective function converges and outputs the optimal subsurface strata joint inversion model. The geothermal anomaly delineation module delineates the geothermal anomaly zone of the target exploration area based on the joint inversion model and multi-parameter fusion identification criteria.

9. The joint inversion apparatus according to claim 8, characterized in that, The multidimensional model building module is also used for: A lateral constraint matrix for layer parameters is constructed to constrain the abrupt differences in resistivity and shear wave velocity between adjacent measuring points in the same stratum. Construct a lateral constraint matrix for layer thickness to constrain abrupt differences in the depth of the same stratum interface between adjacent measuring points; Furthermore, the initial one-dimensional resistivity model and the initial one-dimensional shear wave velocity model are unified onto the same computational grid using grid interpolation technology, and a cross gradient constraint term is established. The cross gradient constraint term is used to constrain the structural similarity between the resistivity model and the shear wave velocity model.

10. The joint inversion apparatus according to claim 9, characterized in that, The cross gradient constraint term is defined as follows: in, These are resistivity model parameters. These are the parameters for the shear wave velocity model. For gradient operators, Represents the horizontal coordinate. Represents the vertical coordinates. For the shear wave velocity model of underground media, This is a resistivity model for underground media.

11. The joint inversion apparatus according to claim 8, characterized in that, The objective function construction module is also used for: The electromagnetic method data fitting term and the micro-motion detection data fitting term are weighted and combined, and combined with the lateral constraint regularization term and the cross gradient structure constraint term to construct a unified inversion objective function; The unified inversion objective function is expressed as follows: in, For electromagnetic data fitting, This is the fitting term for the micro-motion detection data. The layer parameter lateral constraint matrix, The layer thickness lateral constraint matrix is... For the joint model parameter vector, , For regularization parameters, These are the weighting coefficients of the cross-gradient constraint term. For the cross gradient constraint term, This represents the lateral constraint regularization term for layer thickness. Presentation layer parameter lateral constraint regularization term, This represents the cross-gradient structure constraint term.

12. The joint inversion apparatus according to claim 8, characterized in that, The synchronous inversion iteration module is also used for: The L-BFGS algorithm is used for iterative calculation, and the resistivity model parameters and shear wave velocity model parameters are updated simultaneously in each iteration. After each iteration update, the cross gradient between the resistivity model and the shear wave velocity model is calculated, and lateral constraints are applied to the two models respectively. Calculate the unified inversion objective function value and the fitting difference to determine whether the convergence condition is met; If the convergence condition is met, the joint inversion ends and the optimal subsurface strata joint inversion model is output. If the convergence condition is not met, adjust the regularization parameters and weight coefficients, and iterate again.

13. The joint inversion apparatus according to claim 8, characterized in that, The geothermal anomaly delineation module is also used for: Based on the resistivity model and shear wave velocity model output by the joint inversion model, low-resistivity anomaly regions and low-velocity anomaly regions are identified. Analyze resistivity gradients to locate reservoir boundaries and identify high polarizability features; Analyze velocity anisotropy to identify the main direction of fracture development; The identification results are fused with multiple parameters to identify spatial distribution areas with a combination of low resistance, low speed and high polarization. The geothermal anomaly zone is delineated by combining regional geological structural characteristics with constraints.

14. The joint inversion apparatus according to claim 8, characterized in that, The multidimensional model building module is also used for: An initial layered model is established based on regional geological data, borehole data, or inversion results from a single method. Determine the initial resistivity and shear wave velocity parameters for each layer, and establish the initial coupling relationship between resistivity and shear wave velocity.

15. A machine-readable storage medium storing instructions thereon, characterized in that, When executed by the processor, the instruction causes the processor to be configured to perform the joint inversion method based on high-power artificial source electromagnetic method and micro-motion detection data according to any one of claims 1 to 7.

16. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by the processor, it implements the joint inversion method based on high-power artificial source electromagnetic method and micro-motion detection data according to any one of claims 1 to 7.

Images