A time-dimension based seismic wave velocity-resistivity joint inversion method
By using the seismic wave velocity-resistivity joint inversion method, the problem of low accuracy in the joint inversion of seismic exploration and electromagnetic exploration methods has been solved, realizing high-precision three-dimensional imaging and dynamic monitoring of underground fluid resources, and improving exploration efficiency and reliability.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- NORTHWEST ENGINEERING CORPORATION LIMITED
- Filing Date
- 2025-09-05
- Publication Date
- 2026-07-21
Smart Images
Figure CN121028241B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to a time-dimensional seismic wave velocity-resistivity joint inversion method, belonging to the field of geophysical exploration technology. Background Technology
[0002] Currently, the exploration of underground fluid resources mainly employs two geophysical methods: seismic exploration and electromagnetic exploration. Seismic exploration uses artificial seismic sources to generate seismic waves and monitors their propagation and reflection underground to infer the distribution of seismic wave velocity, thereby determining the distribution area of underground fluid resources. Electromagnetic exploration observes the response of natural or artificial electromagnetic fields underground to infer the conductivity or resistivity distribution of the underground medium, and then determines the distribution area of underground fluid resources based on this conductivity or resistivity distribution.
[0003] However, seismic exploration methods are limited by factors such as strong attenuation of high-frequency components and a rapid decrease in resolution with increasing depth, making it difficult to clearly identify deep geological structures and unable to directly distinguish between seismic wave velocity changes caused by fluids and temperature. While electromagnetic exploration methods are sensitive to deep high-conductivity anomalies (such as thermal fluids), their resolution decreases sharply with decreasing frequency, resulting in coarse imaging of deep targets that cannot meet the requirements for fine imaging in deep exploration. Furthermore, electromagnetic exploration data often exhibits ambiguity, meaning that different geological structures may produce similar electromagnetic responses, making it difficult to accurately determine geological structures based solely on electromagnetic exploration data.
[0004] To compensate for the shortcomings of seismic and electromagnetic exploration methods, existing technologies attempt to provide a simple joint interpretation of the inversion results from the two geophysical methods. For example, comparing low resistivity regions obtained from electromagnetic exploration inversion with low wave velocity regions obtained from seismic exploration inversion, or introducing a resistivity model as an initial or constraint condition in electromagnetic exploration inversion. However, these independent or loosely combined inversion methods have significant drawbacks: First, they lack unified constraints. The models obtained from the inversions of the two geophysical methods are often inconsistent in spatial structure, making it difficult to construct a unified geological model. Furthermore, the ambiguity of results from a single geophysical method often leads to contradictory interpretations. Second, they exhibit strong sequential inversion dependence. For instance, if an inversion is performed first based on electromagnetic exploration and then on seismic exploration, the accuracy of the latter inversion result will overly rely on the accuracy of the former, causing errors to be propagated and amplified. Third, in terms of resolution for deep fractures and faults, the two geophysical methods are difficult to complement each other. For example, seismic exploration may locate a fault but cannot determine whether it is filled with fluid, while electromagnetic exploration, although detecting conductive fluid, may not be able to accurately characterize the fault's geometry due to its low resolution. Finally, in terms of temporal dynamic monitoring, current technologies typically rely on manual comparison of geophysical results from different periods to determine changes in geological structures. This is highly subjective, potentially leading to missed detections or misjudgments, and it cannot clearly define the time-varying anomalies corresponding to different geophysical results. For instance, it cannot determine whether a decrease in seismic velocity corresponds to a decrease in resistivity (increased fluid).
[0005] In summary, existing technologies for the joint inversion of seismic and electromagnetic exploration methods suffer from drawbacks such as low accuracy of inversion results, difficulty in accurately identifying deep geological structures, and inability to clearly define the dynamic changes in geological structures. Summary of the Invention
[0006] This invention provides a time-dimensional seismic wave velocity-resistivity joint inversion method, which can solve the problems of low accuracy of existing inversion results, difficulty in accurately identifying deep geological structures, and inability to clearly identify the dynamic changes of geological structures.
[0007] This invention provides a time-dimensional seismic wave velocity-resistivity joint inversion method, the method comprising:
[0008] S1. Construct a seismic wave velocity distribution model based on the seismic exploration data of the target area, and construct a resistivity distribution model based on the electromagnetic exploration data of the target area;
[0009] S2. Construct constraints based on the seismic wave velocity distribution model and the resistivity distribution model, and determine the objective function for the joint inversion of seismic wave velocity and resistivity based on the constraints; the constraints include prediction bias fitting terms, model structure coupling terms, spatial regularization terms, temporal difference constraints, and model boundary consistency constraints.
[0010] S3. Based on the seismic exploration data and electromagnetic exploration data from multiple periods, update and iterate the seismic wave velocity distribution model and the resistivity distribution model until the objective function is minimized;
[0011] S4. Based on the seismic wave velocity distribution model and resistivity distribution model corresponding to the minimization of the objective function, determine the seismic wave velocity inversion results and resistivity inversion results of the target area in the multiple periods, respectively.
[0012] Optionally, in S2, the objective function for the joint inversion of seismic wave velocity and resistivity is determined based on the aforementioned constraints, specifically including:
[0013] The objective function for joint seismic wave velocity-resistivity inversion is obtained by weighted summation of the prediction deviation fitting term, the model structure coupling term, the spatial regularization term, the temporal difference constraint term, and the model boundary consistency constraint term.
[0014] Optionally, the seismic exploration data includes the observed first arrival and travel times of seismic waves at multiple seismic measuring points, and the electromagnetic exploration data includes the observed apparent resistivity of the subsurface medium at multiple electromagnetic measuring points.
[0015] Optionally, in S2, a prediction deviation fitting term is constructed based on the seismic wave velocity distribution model and the resistivity distribution model, specifically including:
[0016] The calculated first arrival travel time at the multiple seismic measuring points is calculated based on the seismic wave velocity distribution model, and the calculated apparent resistance at the multiple electromagnetic measuring points is calculated based on the resistivity distribution model.
[0017] Based on the calculated first arrival time, the calculated apparent resistivity, the observed first arrival time, and the observed apparent resistivity, a prediction bias fitting term is constructed.
[0018] Optionally, in S2, a model structure coupling term is constructed based on the seismic wave velocity distribution model and the resistivity distribution model, specifically including:
[0019] Based on the gradient vectors at different locations in the seismic wave velocity distribution model and the resistivity distribution model, a model structure coupling term is constructed using the model coupling method.
[0020] Optionally, in S2, a spatial regularization term is constructed based on the seismic wave velocity distribution model and the resistivity distribution model, specifically including:
[0021] Based on the gradient vectors at different locations in the seismic wave velocity distribution model and the resistivity distribution model, a spatial regularization term is constructed using the smoothness constraint method.
[0022] Optionally, in S2, a time difference constraint term is constructed based on the seismic wave velocity distribution model and the resistivity distribution model, specifically including:
[0023] The seismic wave velocity distribution of the target area during the multiple periods is calculated based on the seismic wave velocity distribution model, and the resistivity distribution of the target area during the multiple periods is calculated based on the resistivity distribution model.
[0024] Based on the seismic wave velocity and resistivity distributions corresponding to the multiple periods, a time difference constraint term is constructed.
[0025] Optionally, in S2, a model boundary consistency constraint term is constructed based on the seismic wave velocity distribution model and the resistivity distribution model, specifically including:
[0026] Based on the seismic wave velocity at different locations in the seismic wave velocity distribution model and the resistivity at different locations in the resistivity distribution model, a model boundary consistency constraint term is constructed using a unit step function.
[0027] Optionally, S3 involves updating and iterating the seismic wave velocity distribution model and the resistivity distribution model, specifically including:
[0028] The seismic wave velocity distribution model and the resistivity distribution model are updated synchronously in each iteration;
[0029] Alternatively, the seismic wave velocity distribution model and the resistivity distribution model can be updated alternately in multiple iterations.
[0030] Optionally, the seismic wave velocity distribution model and the resistivity distribution model are updated synchronously in each iteration, specifically including:
[0031] In each iteration, the model parameters of the seismic wave velocity distribution model and the resistivity distribution model for each period are updated synchronously.
[0032] Alternatively, in one iteration, only the model parameters of the seismic wave velocity distribution model and the resistivity distribution model for one period are updated, and in the next iteration, the model parameters of the seismic wave velocity distribution model and the resistivity distribution model for the next period are updated based on the model parameters corresponding to the previous period.
[0033] The beneficial effects that this invention can produce include:
[0034] (1) Simultaneous inversion of two physical property models
[0035] This invention constructs an objective function based on seismic velocity distribution and resistivity distribution models, including prediction bias fitting terms, model structure coupling terms, spatial regularization terms, temporal difference constraint terms, and model boundary consistency constraint terms. The objective function is then used to jointly invert the seismic velocity and resistivity distributions of the target area over multiple periods, enabling simultaneous inversion of both models. By ensuring structural consistency and coupling between the seismic velocity and resistivity distribution models during the inversion process and extending the inversion process in the temporal dimension, high-precision imaging of the three-dimensional structure and dynamic changes of underground fluid resources can be achieved. This improves the detection accuracy and spatiotemporal monitoring capabilities of underground fluid resources, enhancing the safety and efficiency of underground fluid resource exploration and development.
[0036] (2) Improve exploration accuracy and reduce ambiguity
[0037] This invention achieves a high degree of geometric consistency between the inverted seismic wave velocity distribution model and resistivity distribution model by structurally coupling two models representing different physical properties. This effectively reduces the non-uniqueness of single-model inversion, improves the accuracy of characterizing the real underground structure, and is closer to the true values in terms of anomaly location and physical property values than individual inversions. Especially at key structures such as faults and fractures, joint imaging can more clearly delineate their spatial contours and depth of extension. Compared to interpreting the inversion results of each model separately, integrating the inversion results of each model avoids contradictions and can provide a more unified and reliable geological interpretation.
[0038] (3) Improve the quality of three-dimensional structural imaging
[0039] In this invention, the seismic wave velocity distribution is used to constrain the resistivity distribution and refine the boundary of the electrical layer, which makes up for the problem of insufficient resolution of low resistivity layers in electromagnetic exploration methods. At the same time, the resistivity distribution also provides lithology and fluid indicators for seismic imaging, which can distinguish the causes of seismic reflection / velocity anomalies. This makes the three-dimensional model obtained by joint inversion contain richer information, and the imaging clarity and resolution are significantly better than the results of single model inversion.
[0040] (4) Strong dynamic response recognition capability
[0041] This invention, through spatiotemporal dynamic monitoring and combined imaging, can keenly capture the temporal variation characteristics of underground fluid resources. It is more sensitive to minute changes, more accurate in locating them, and can distinguish different types of changes, thus providing strong support for the observation and analysis of underground fluid dynamic processes.
[0042] (5) Wide range of applications and high reliability
[0043] This invention does not rely on specific sites or equipment and is applicable to the exploration of various underground fluid resources. In practical applications, it demonstrates far superior performance compared to traditional single geophysical exploration methods: it can provide a more comprehensive perspective on parameters, more stable inversion results, and reduce the uncertainty of interpretation.
[0044] (6) It has good noise resistance and reliability.
[0045] This invention employs noise adaptation and multi-constraint balancing strategies, which have good noise resistance and reliability for complex field data, and can still produce reasonable models even when there are some missing data or noise. Attached Figure Description
[0046] Figure 1 A flowchart of the time-dimensional seismic wave velocity-resistivity joint inversion method provided in this embodiment of the invention;
[0047] Figure 2 The cross gradient variation diagram provided in the embodiments of the present invention;
[0048] Figure 3 This is a differential cross-sectional view provided for an embodiment of the present invention. Detailed Implementation
[0049] The present invention will now be described in detail with reference to the embodiments, but the present invention is not limited to these embodiments.
[0050] This invention provides a time-dimensional seismic wave velocity-resistivity joint inversion method, such as... Figure 1 As shown, the method includes:
[0051] S1. Construct a seismic wave velocity distribution model based on the seismic exploration data of the target area, and construct a resistivity distribution model based on the electromagnetic exploration data of the target area.
[0052] The seismic exploration data includes the observed first arrival and travel times of seismic waves at multiple seismic monitoring points within the target area, while the electromagnetic exploration data includes the observed apparent resistivity of the subsurface medium at multiple electromagnetic monitoring points within the target area.
[0053] This embodiment first performs preprocessing on seismic exploration data and electromagnetic exploration data, including denoising, distortion correction, and static correction, to ensure the comparability of different types of data. Then, based on the preprocessed seismic exploration data and electromagnetic exploration data, seismic wave velocity distribution models and resistivity distribution models are constructed respectively, ensuring that the seismic wave velocity distribution models and resistivity distribution models can represent the initial distribution of seismic wave velocity and resistivity under a unified three-dimensional grid and coordinate system.
[0054] Seismic wave velocity distribution models and resistivity distribution models can be constructed based on prior geological information, or they can be obtained by inverting the seismic wave velocity distribution and resistivity distribution separately. At the same time, the coordinates of the seismic wave velocity distribution models and resistivity distribution models should be unified and the scale normalized to facilitate mutual comparison and constraints during subsequent joint inversion.
[0055] It is worth noting that in existing technologies, because seismic wave velocity distribution models and resistivity distribution models represent different physical properties, different modeling discrete grids are often required during modeling. For example, seismic wave velocity distribution models require fine grids while resistivity distribution models require large-scale coarse grids, resulting in significant differences in grid scales. This makes subsequent joint inversion inconvenient. However, if the grids are strictly unified to meet joint inversion requirements, it may lead to a mismatch between the model's resolution and that of its exploration data. To address these issues, this embodiment employs a heterogeneous grid unification strategy. Through coordinate transformation or interpolation mapping, the two models representing different physical properties are projected onto a unified computational domain, ensuring grid consistency or correspondence, thus facilitating subsequent joint inversion. For example, a multi-scale model containing nested fine and coarse grids can be constructed, allowing the two models to share nodes in key areas. This preserves the resolution required by each data set while supporting mutual constraints between models, resolving the aforementioned problems in existing technologies.
[0056] S2. Construct constraints based on the seismic wave velocity distribution model and resistivity distribution model, and determine the objective function for the joint inversion of seismic wave velocity and resistivity based on the constraints.
[0057] In this embodiment, the constraints include prediction bias fitting terms, model structure coupling terms, spatial regularization terms, temporal difference constraints, and model boundary consistency constraints.
[0058] This embodiment obtains the objective function for joint inversion of seismic wave velocity and resistivity by weighted summation of the prediction deviation fitting term, model structure coupling term, spatial regularization term, temporal difference constraint term, and model boundary consistency constraint term.
[0059] In this embodiment, the seismic wave velocity distribution model is represented as follows: The resistivity distribution model is expressed as Then the objective function can be expressed as:
[0060] (1)
[0061] In formula (1), Indicated based on seismic wave velocity distribution model and resistivity distribution model Construct the objective function; This represents the prediction bias fitting term; Indicates the coupling terms in the model structure; Represents the space regularization term; Indicates the time difference constraint term; This represents the model boundary consistency constraint term; , , and These represent the weight coefficients of the model structure coupling term, spatial regularization term, temporal difference constraint term, and model boundary consistency constraint term, respectively, and are used to adjust the relative influence of the corresponding terms in the objective function.
[0062] The following details the construction process and expressions of each term in the objective function:
[0063] 1) Prediction bias fitting term
[0064] This embodiment uses a seismic wave velocity distribution model to forward simulate seismic exploration data (i.e., observed first arrival travel times) and a resistivity distribution model to forward simulate electromagnetic exploration data (i.e., observed apparent resistivity). Through the above forward simulation process, the calculated first arrival travel times at multiple seismic measuring points and the calculated apparent resistivity at multiple electromagnetic measuring points can be calculated. The calculated first arrival travel times and calculated apparent resistivity will be used as leading operators in the subsequent model inversion.
[0065] To obtain the first arrival travel time, this embodiment uses the Eikonal equation to approximate the propagation of the first arrival travel time of seismic waves. The Eikonal equation is expressed as:
[0066] (2)
[0067] In formula (2), This indicates the earthquake source in the seismic wave velocity distribution model. The position in the middle The travel time gradient of the seismic wave generated at the point of origin, i.e., the normal vector of the wavefront (seconds / meter). express The Euclidean norm of the gradient during the walk-time; Indicates that seismic waves are in The propagation velocity in the subsurface medium, i.e., the seismic wave velocity distribution model exist The value at that location (meters per second).
[0068] Equation (1) indicates that the travel time field in the subsurface medium satisfies the isotropic Eikonal equation, meaning the magnitude of the travel time gradient is equal to the reciprocal of the wave velocity. By solving equation (1) using a fast travel algorithm, the calculated first arrival travel times at multiple seismic monitoring points can be obtained. .
[0069] To obtain the calculated apparent resistance, this embodiment is based on a resistivity distribution model. The apparent resistivity is determined by inverting the response to electromagnetic wave propagation. For an angular frequency of... The apparent resistivity of a one-dimensional isotropic half-space can be expressed by a magnetotelluric plane wave as:
[0070] (3)
[0071] In formula (3), Represents apparent resistivity (ohm·m); Represents the vacuum permeability (4π×10⁻⁶) -7 Henry / meter); Angular frequency (radians per second); and They represent angular frequencies respectively. The observed electric field amplitude (volts / meter) and magnetic field amplitude (amperes / meter) are shown below.
[0072] Equation (3) shows that the equivalent resistivity response of the underground medium to electromagnetic waves can be calculated by the ratio of the electric field amplitude to the magnetic field amplitude. This is relevant to the resistivity distribution model. In the actual forward modeling, it is necessary to obtain the electric field amplitude and magnetic field amplitude at each electromagnetic measuring point by numerically solving Maxwell's equations, and then obtain the calculated apparent resistance at multiple electromagnetic measuring points according to formula (3). .
[0073] Then, in this embodiment, a prediction deviation fitting term is constructed based on the calculation of the first arrival time, the calculation of the apparent resistance, the observation of the first arrival time, and the observation of the apparent resistivity.
[0074] The prediction bias fitting term is used to quantify the degree of deviation between the calculated data obtained from the forward modeling and the observed data obtained from actual exploration. It includes two parts: first arrival travel time deviation and apparent resistivity deviation, which can be expressed in weighted least squares form as follows:
[0075] (4)
[0076] In formula (4), This represents the prediction bias fitting term; and They represent the first Calculated first arrival travel time (seconds) and observed first arrival travel time (seconds) at each seismic monitoring point; This indicates the total number of seismic monitoring points. and They represent the first Calculated apparent resistivity (ohm·m) and observed apparent resistivity (ohm·m) at each seismic monitoring point; This indicates the total number of electromagnetic measuring points.
[0077] In formula (4), the sum of squares of the differences between the calculated first arrival travel time and the observed first arrival travel time (i.e., the first arrival travel time deviation) represents the seismic wave velocity distribution model. The degree to which the predicted data deviates from the actual data; the sum of squares of the differences between apparent resistivity and observed apparent resistivity (i.e., apparent resistivity deviation) represents the resistivity distribution model. The degree to which the predicted data deviates from the actual data.
[0078] It is worth noting that this embodiment does not limit the weight values corresponding to the first arrival travel time deviation and the apparent resistance deviation. In practical applications, the first arrival travel time deviation and the apparent resistance deviation can be normalized or weighted according to industry standards and exploration requirements to determine their weights (for example, using the reciprocal of the observation error as the weight). Formula (4) only gives an example of assigning the same weight to the first arrival travel time deviation and the apparent resistance deviation.
[0079] 2) Model structure coupling terms
[0080] This embodiment is based on the seismic wave velocity distribution model. and resistivity distribution model The gradient vectors at different positions in the model are obtained, and the model structure coupling term is constructed using the model coupling method.
[0081] Model structure coupling terms are used to characterize the seismic wave velocity distribution model. and resistivity distribution model The structural similarity between them. Model coupling methods include cross-gradient methods, parameter coupling constraint methods, Gramian constraint methods, etc.
[0082] This embodiment uses the cross-gradient method as an example to illustrate the construction of model structure coupling terms. The cross-gradient method calculates the cross product of the gradients of two models, and the norm of the cross product reflects the degree of inconsistency between the gradient fields of the two models. The model structure coupling term based on the cross-gradient method can be expressed as the L2 integral of the cross product vector over the entire model space:
[0083] (5)
[0084] In formula (5), Indicates the coupling terms in the model structure; Seismic wave velocity distribution model exist The gradient vector at that point; Representing the resistivity distribution model exist The gradient vector at the given location; "×" indicates the cross product operation; " represents the squared 2 norm of the vector after the cross product operation; Ω represents the integration domain, i.e., the entire model space.
[0085] According to the cross product property, when and When the directions are the same or opposite (i.e., the two gradient vectors are parallel), the cross product is zero, indicating that the two models are parallel. The structural interfaces at the points are aligned; conversely, if there is an angle between the two gradient vectors, a non-zero cross product is generated, and the larger the value of the cross product, the greater the structural difference between the two models. Therefore, by minimizing the model structural coupling term... This can make the gradient fields of the two models as parallel as possible, thereby enhancing the seismic wave velocity distribution model. and resistivity distribution model The common structural features of the two models. The model structure coupling term ensures consistency in the structural direction of the two models by reducing the angle between the gradient vectors of the two models. When the angle is zero or 180°, the model structure coupling term... It reaches the minimum value.
[0086] This embodiment introduces a model structure coupling term into the objective function. Realize the seismic wave velocity distribution model and resistivity distribution model Consistent coupling for the same geological structure. The cross-gradient method does not rely on prior rock property relationships. By minimizing the difference in gradient directions between the two models, it forces the two models to have the same structural profile (such as fault planes, fracture zones, etc.) on the same geological structure. This ensures that the same geological structure in the subsequent inversion results remains consistent in the two models, improving the reliability of the inversion results.
[0087] Besides the cross-gradient method, this embodiment can also utilize other model coupling methods to construct model structure coupling terms. For example, the parameter coupling constraint method, which uses rock physics relationships to construct the seismic wave velocity distribution model... and resistivity distribution model In related fields, for example, using empirical formulas to initially convert seismic wave velocities into resistivity ranges, and then minimizing the difference between the two in subsequent inversions, this falls under the category of property-based coupling. When reliable lithological or fluid saturation relationships exist, parameter coupling constraint methods can directly constrain parameter amplitudes, but compared to cross-gradient methods, this method has a narrower applicability. Another example is the Gramian constraint method, which uses generalized correlation functions to correlate certain features of different property models (such as gradient fields and anomaly distributions) to achieve structural coupling.
[0088] 3) Spatial regularization term
[0089] This embodiment is based on the seismic wave velocity distribution model. and resistivity distribution model The gradient vectors at different positions in the vector are obtained, and a space regularization term is constructed using the smoothness constraint method.
[0090] Spatial regularization terms are used to constrain the smoothness or prior properties of model parameters in space, thereby addressing the ill-conditioned nature of model parameter inversion and improving the stability of the inversion results. This embodiment uses the Tikhonov regularization form of the second norm as an example, applying it to a seismic wave velocity distribution model... and resistivity distribution model The summation of the squares of the first-order gradients can be used to express the space regularization term as:
[0091] (6)
[0092] In formula (6), Represents the space regularization term; Seismic wave velocity distribution model exist The squared magnitude of the gradient vector at that point; Representing the resistivity distribution model exist The squared magnitude of the gradient vector at a given point; Ω represents the integration domain, i.e., the entire model space.
[0093] Formula (6) penalizes the spatial variation of the two models, tending to produce smooth model solutions and avoiding excessive oscillations or unreasonable spikes. In practical applications, higher-order smoothness constraints (such as the second derivative of the model) or constraints based on geological prior models can also be used to construct spatial regularization terms as needed. Formula (6) only uses a simple first-order gradient squared term as an example to illustrate the smoothing constraint effect of the spatial regularization term.
[0094] 4) Time difference constraint term
[0095] This embodiment is based on a seismic wave velocity distribution model. The seismic wave velocity distribution in the target area over multiple periods was calculated, and based on the resistivity distribution model... Calculate the resistivity distribution of the target area over multiple periods; then, construct time difference constraints based on the seismic wave velocity and resistivity distributions corresponding to the multiple periods.
[0096] The time difference constraint term is used to constrain the model's variation amplitude at different times, ensuring that the model's evolution is stable and reasonable between adjacent time periods during spatiotemporal dynamic monitoring (4D). For those with... The 4D inversion of periodic monitoring data can express the time difference constraint term as the L2 norm of the difference between adjacent periods of the model:
[0097] (7)
[0098] In formula (7), Indicates the time difference constraint term; Seismic wave velocity distribution model In the Period and the Seismic wave velocity distribution during the period (m / s); Representing the resistivity distribution model In the Period and the Resistivity distribution over a period of time (ohm-meter); This represents the sum of squares of the time intervals between adjacent time intervals of all grid cells in the corresponding model; This represents the total for all periods.
[0099] Formula (7) applies a squared penalty to the differences between adjacent periods in the model, encouraging the model parameters to gradually smooth out changes over time. In other words, the time difference constraint term suppresses unrealistic and drastic changes in the results of two adjacent explorations during the inversion process, thus reflecting the true evolution of subsurface properties. In this way, the model can update its parameters accordingly for different monitoring periods, allowing the model characteristics of adjacent periods to be reflected by the time difference constraint term. In connection with this, the time evolution trend of the model is subject to smoothing constraints.
[0100] This embodiment, by adding a time dimension to the spatial joint inversion, enables dynamic monitoring and imaging of subsurface fluid resources changing over time. That is, this embodiment not only acquires the seismic wave velocity and resistivity distributions in three-dimensional space, but also reflects the evolution of these distributions by comparing changes over multiple periods, identifying anomalies in each period, and achieving continuous monitoring of subsurface fluid resources, including resource migration and temperature changes.
[0101] 5) Model boundary consistency constraints
[0102] This embodiment is based on the seismic wave velocity distribution model. Seismic wave velocity and resistivity distribution models at different locations in China The resistivity at different locations is used to construct the model boundary consistency constraint terms using a unit step function.
[0103] Model boundary consistency constraints are used to enhance the consistency between the seismic velocity distribution model and the resistivity distribution model in describing the boundaries of the same geological structure. Geological regions containing underground fluid resources often exhibit common interface characteristics across different physical properties; for example, high-temperature liquid-bearing regions may correspond to lower seismic velocities and lower resistivity. This embodiment utilizes a unit heaviside step function to construct model boundary consistency constraints, ensuring that the two models consistently characterize the boundaries of the same geological structure.
[0104] In this embodiment, threshold values for seismic wave velocity and resistivity are first set. For example, based on typical values in the background or anomalous region, the threshold value for seismic wave velocity is denoted as... And the threshold of resistivity is denoted as Then use the unit step function (when hour ,when hour Convert both models to binary representations, then calculate the sum of squares of the differences between their binary representations to express the boundary consistency constraint term as follows:
[0105] (8)
[0106] In formula (8), This represents the model boundary consistency constraint term; Seismic wave velocity distribution model exist Seismic wave velocity at the location; A threshold representing seismic wave velocity; Representing the resistivity distribution model exist resistivity at that point; The threshold value representing resistivity; Represents the unit step function; Seismic wave velocity distribution model exist Binary representation of the location; Representing the resistivity distribution model exist The binary representation of Ω; Ω represents the integration domain, i.e., the entire model space.
[0107] and The squared difference in the integral measures the difference in the two models' classification of geological bodies at a given location. When both models simultaneously show similar attributes at a location (e.g., both belong to the background region or both belong to the anomalous region), then... and If the output values are the same, the difference between the two is 0; if one model identifies a region as an anomaly while the other identifies it as background, then... and Output different values, with the difference being 1. This generates a penalty. By minimizing... This allows the two models to delineate consistent geological structural boundaries. This helps the two models identify common geological structures, such as faults and reservoir boundaries, thereby improving the reliability of subsequent joint inversion.
[0108] S3. Based on seismic exploration data and electromagnetic exploration data from multiple periods, a seismic wave velocity distribution model was developed. and resistivity distribution model The process is repeated until the objective function is minimized.
[0109] This embodiment establishes the objective function. Then, an optimization algorithm was used to iteratively solve the seismic wave velocity distribution model. and resistivity distribution model The optimal solution for the model parameters, in order to minimize the objective function. This allows for the joint inversion of the two models. In practical applications, the seismic wave velocity distribution model can be updated synchronously in each iteration. and resistivity distribution model The model parameters can also be updated alternately in multiple iterations of the seismic wave velocity distribution model. and resistivity distribution model The model parameters. The specific processes for the two update methods are as follows:
[0110] 1) Synchronously update the model parameters of both models.
[0111] The objective function is computed simultaneously in each iteration. Seismic wave velocity distribution model and resistivity distribution model The gradients (partial derivatives) are used to obtain the seismic wave velocity distribution model. and resistivity distribution model The update direction is determined, and the two sets of model parameters corresponding to the two models are simultaneously corrected according to the selected optimization algorithm (such as the conjugate gradient method, L-BFGS method, or quasi-Newton method). This synchronous update method utilizes the objective function... Various terms (especially model structure coupling terms) and model boundary consistency constraints The coupling information of the two models can be comprehensively considered in a single iteration, thereby accelerating the convergence speed.
[0112] 2) Alternately update the model parameters of the two models.
[0113] In one iteration, the resistivity distribution model is first fixed. With the model parameters unchanged, the seismic wave velocity distribution model is updated. To minimize the objective function using model parameters ; then fix the updated seismic wave velocity distribution model The model parameters remain unchanged, and the resistivity distribution model is updated. To minimize the objective function using model parameters This alternating update process gradually satisfies the joint constraints between the two models. This method optimizes only one set of model parameters at a time, which helps ensure stable convergence. However, each global iteration involves two updates, resulting in slower convergence compared to synchronous updates. In practical applications, the choice between synchronous and alternating updates can be made flexibly based on factors such as problem size and coupling strength.
[0114] For iterative optimization along the time dimension, this embodiment can simultaneously invert the model parameters corresponding to each period, that is, the prediction bias fitting terms for all periods. and time difference constraint Included in the objective function In each iteration, the model parameters for each period are updated synchronously. This means that during iteration, the model is not only adjusted based on the seismic and electromagnetic exploration data of the current period, but also takes into account the constraints of the differences between the model and the model of adjacent periods, thus naturally embedding a smooth prior in the time dimension.
[0115] If computational resources are limited, a sequential inversion method with rolling optimization for each period can also be adopted. First, the baseline model is inverted, and then constraint terms that differ from the previous model are added to the objective function of each subsequent period to guide the inversion. That is, in one iteration, only the model parameters of the two models corresponding to one period are updated, and in the next iteration, the model parameters of the two models corresponding to the next period are updated according to the model parameters corresponding to the previous period.
[0116] It is worth noting that, regardless of whether simultaneous joint inversion or rolling sequential inversion is used in this embodiment, compared with the existing technology which only performs separate inversion for each period based on the monitoring data of each period and then compares the periods, this embodiment can effectively suppress noise artifacts and highlight the true temporal changes. Therefore, this embodiment significantly improves the accuracy of time dimension monitoring compared with the existing technology.
[0117] Time difference constraint The embedding of these technologies ensures the continuous and controllable changes in the time dimension during the inversion process, improving the stability and physical reliability of 4D inversion and imaging.
[0118] Furthermore, this embodiment can adaptively adjust the objective function during the inversion solution process. The weights and parameters corresponding to each term are defined to dynamically balance the constraints of each term for different iteration stages and different regions of the model. For example, increasing the model structure coupling term in the initial stage... The weights are adjusted to quickly unify the structural framework of the two models, and the coupling terms of the model structure are gradually reduced in the mid-to-late stages. The weights are adjusted to finely fit the data; for example, the constraints are reduced in high-noise regions to avoid overfitting the noise. This adaptive mechanism ensures stable convergence of the inversion process without sacrificing detail and realism.
[0119] Furthermore, this embodiment addresses the main noise types in seismic and electromagnetic exploration data (e.g., seismic exploration data mainly contains high-frequency random noise, while electromagnetic exploration data mainly contains instrument errors and ground noise), and in the objective function... This paper introduces a source-specific weighting and robust estimation strategy, which employs a corresponding noise model and weighting matrix for each type of exploration data to improve robustness to outlier data points. For example, Huber loss or L1 norm can be used for electromagnetic exploration data containing outliers to improve noise resistance; while L2 norm can be used to adjust weights to reduce the impact of uncertain picking for random noise in seismic exploration data. The error levels of the two types of exploration data are automatically estimated and used to update the corresponding weights, enabling the joint inversion process to "adapt" to the signal-to-noise characteristics of different data sources, maximizing the extraction of useful information from the exploration data without over-relying on observations with high noise levels.
[0120] Furthermore, this embodiment can also apply constraints to local areas and improve the model mesh during the inversion process. In some cases, the coverage of seismic exploration data and electromagnetic exploration data differs, or the area of interest is limited. A local cross-gradient constraint scheme can be adopted, that is, structural coupling constraints are applied only to overlapping areas or areas of interest, while other areas of the model are allowed to be inverted freely. This ensures the structural uniformity of key areas and reduces the computational overhead and constraint conflicts in irrelevant areas. For the heterogeneous meshes of the two models, multi-grid partitioning inversion can be adopted. For example, a fine mesh can be used in areas with dense seismic exploration data, while electromagnetic exploration data only covers a coarse mesh. Special interpolation processing is performed on the boundary areas when calculating the cross gradient, which can ensure the adaptability of the algorithm to the complex observation meshes in reality.
[0121] S4. Based on the objective function Minimize the corresponding seismic wave velocity distribution model and resistivity distribution model The model parameters were used to determine the seismic wave velocity inversion results and resistivity inversion results for the target area at multiple periods.
[0122] The seismic wave velocity inversion results include the three-dimensional spatial distribution of seismic wave velocities at various periods. Differences in seismic wave velocity distribution between adjacent periods The resistivity inversion results include the three-dimensional spatial distribution of resistivity at various times. Differences in resistivity distribution between adjacent periods .
[0123] To present the inversion results intuitively, the three-dimensional spatial distribution of seismic wave velocity and resistivity can be drawn into a three-dimensional volume, or visualized using two-dimensional profiles (such as planar slices or vertical profiles); the differences in seismic wave velocity distribution and resistivity distribution between adjacent periods can be visualized in the form of isosurfaces or difference cloud maps to show the evolutionary location and magnitude of geological bodies.
[0124] Furthermore, the seismic wave velocity inversion results and resistivity inversion results can be jointly interpreted. For example, by overlaying the three-dimensional spatial distributions of seismic wave velocity and resistivity, anomaly regions where both seismic wave velocity and resistivity decrease simultaneously can be identified. These anomaly regions often correspond to reservoirs containing thermal fluids in underground fluid resource exploration. Through visualization and joint interpretation of the above inversion results, the trend of geological structure changes over time can be intuitively revealed, providing a scientific basis for the assessment, monitoring, exploration, and development of underground fluid resources.
[0125] To verify the effectiveness of the method described in this embodiment, an application test was conducted in a typical geothermal field area. This area contains multiple concealed faults, deep water-conducting fractures, and hydrothermal anomalies, making it difficult for conventional seismic or electromagnetic exploration methods to simultaneously identify changes in their geometric structure and fluid properties. The process of using the method described in this embodiment is as follows:
[0126] 1) Data Acquisition and Preprocessing
[0127] Several two-dimensional seismic survey lines were laid out in the target area, and more than 20 magnetotelluric (MT) measuring points were deployed to acquire first arrival travel time data of seismic waves and MT apparent resistivity data with a frequency range of 0.01Hz to 1000Hz. All data were statically corrected, denoised, and uniformly projected onto the three-dimensional model mesh.
[0128] 2) Initial Model Construction
[0129] Using existing drilling data and geological survey results, initial seismic velocity distribution models and resistivity distribution models were established. The seismic velocity distribution model was constructed independently based on first-arrival travel-time data, while the resistivity distribution model was derived from the inversion projection of the apparent resistivity data from the Median Transit Authority (MT). Both models were interpolated and resampled to a shared 3D computational grid in a unified coordinate system with a spatial resolution of 50m.
[0130] 3) Implement cross-gradient joint inversion
[0131] Construct an objective function that includes a prediction bias fitting term, a model structure coupling term, a spatial regularization term, a temporal difference constraint term, and a model boundary consistency constraint term. The weights of each term are set as follows: the weight of the prediction bias fitting term is 1.0, the weight of the model structure coupling term is 0.5 (high initially, decreasing to 0.2 later), the weight of the spatial regularization term is 0.05, the weight of the temporal difference constraint term is 0.2, and the weight of the model boundary consistency constraint term is 0.1.
[0132] Then, based on data from three periods in the target area—before water injection (T0 period), 6 months after water injection (T1 period), and 12 months after water injection (T2 period)—the conjugate gradient method was used to simultaneously optimize and iterate the objective function, achieving joint inversion of the two models. Iteration continued until the objective function met the convergence threshold (relative error <0.5%), at which point the iteration terminated. The cross-gradient change diagram during the inversion process is shown below. Figure 2 As shown,
[0133] 4) Results Analysis and Verification
[0134] The difference profile of the inversion results is shown below. Figure 3 As shown, by Figure 3 It can be seen that the joint inversion model shows a dipping fault surface at a depth of about 2.1 km underground in period T0. The south side of the fault has a low seismic wave velocity (<3.0 km / s) and low resistivity (<10 Ω·m) anomaly, which is speculated to be a hydrothermal water-conducting channel. This structure is blurred in the seismic wave velocity inversion or resistivity inversion alone, but the joint inversion significantly improves the clarity of boundary identification.
[0135] In the T1-T2 difference model, the resistivity of the anomalous region decreased by approximately 15%, and the seismic wave velocity decreased by approximately 7%, indicating that water injection led to an increase in saturation and temperature, with a significant response from the hydrothermal system. Drilling verification showed that the temperature in this region increased by 9.3°C, and the pore pressure increased by 1.2 MPa, which is highly consistent with the inversion changes, verifying the sensitivity and spatial positioning capability of the method described in this embodiment to changes in the thermal reservoir.
[0136] The above verification demonstrates the following advantages of the method described in this embodiment:
[0137] (1) Simultaneous inversion of two physical property models
[0138] This method constructs an objective function based on seismic velocity distribution and resistivity distribution models, including prediction bias fitting terms, model structure coupling terms, spatial regularization terms, temporal difference constraint terms, and model boundary consistency constraint terms. The objective function is then used to jointly invert the seismic velocity and resistivity distributions of the target area over multiple periods, enabling simultaneous inversion of both models. By ensuring structural consistency and coupling between the seismic velocity and resistivity models during the inversion process and extending the inversion process in the temporal dimension, high-precision imaging of the three-dimensional structure and dynamic changes of underground fluid resources can be achieved. This improves the detection accuracy and spatiotemporal monitoring capabilities of underground fluid resources, enhancing the safety and efficiency of underground fluid resource exploration and development.
[0139] (2) Improve exploration accuracy and reduce ambiguity
[0140] This method structurally couples two models representing different physical properties, ensuring a high degree of geometric consistency between the inverted seismic velocity distribution model and resistivity distribution model. This effectively reduces the non-uniqueness of single-model inversion, improves the accuracy of characterizing the actual subsurface structure, and is closer to the true values in terms of anomaly location and physical property values than individual inversions. Especially at key structures such as faults and fractures, joint imaging can more clearly delineate their spatial contours and depth of extension. Compared to interpreting the inversion results of each model separately, integrating the inversion results of each model avoids contradictions and provides a more unified and reliable geological interpretation.
[0141] (3) Improve the quality of three-dimensional structural imaging
[0142] This method utilizes seismic wave velocity distribution to constrain resistivity distribution and refines the boundaries of electrical layers, compensating for the insufficient resolution of low-resistivity layers in electromagnetic exploration methods. Simultaneously, the resistivity distribution provides lithological and fluid indicators for seismic imaging, distinguishing the causes of seismic reflection / velocity anomalies (e.g., identifying whether high-velocity anomalies are due to bedrock uplift or cooling zones). Therefore, the 3D model obtained through joint inversion contains richer information: both structural details and physical property differences, with more accurate stratigraphic correlation, resulting in significantly better imaging clarity and resolution than single-model inversion results. For example, in a simulation example, joint inversion successfully imaged deep dipping faults that were difficult to identify using single methods and correctly displayed the high-resistivity / high-velocity block structures existing simultaneously on both sides.
[0143] (4) Strong dynamic response recognition capability
[0144] This method, through spatiotemporal dynamic monitoring and combined imaging, can sensitively capture the temporal variation characteristics of underground fluid resources. During water injection and recharge, and production extraction, seismic wave velocity and resistivity distributions often change synchronously. For example, increased temperature and pore fluid volume can decrease both seismic wave velocity and resistivity. This method provides cross-validation of the causes of these changes by fusing these two time-varying information sources. For instance, when a significant decrease in resistivity over time is detected in a certain area, this method simultaneously checks whether there is a corresponding decrease in seismic wave velocity, thus determining that the area is likely due to fluid intrusion. If only a single parameter changes, the joint inversion can also suppress spurious changes and output more reliable results through cross-constraints. Compared with existing technologies, this method is more sensitive to minute changes, more accurate in locating them, and can distinguish between different types of changes (such as seismic wave velocity changes caused by pressure changes and resistivity changes caused by water saturation changes). This is of great value for long-term monitoring and anomaly early warning in geothermal fields, and is also applicable to monitoring shale gas fracturing. Observation and analysis of underground fluid dynamic processes, such as monitoring of geological sealing leaks.
[0145] (5) Wide range of applications and high reliability
[0146] This method is not dependent on specific sites or equipment and is applicable to the exploration of various underground fluid resources. For geothermal exploration, it can be used to identify underground hydrothermal channels, fractured water-conducting zones, and assess the integrity of geothermal reservoir caprocks, helping to optimize drilling site selection and development strategies. For unconventional oil and gas (such as shale gas), it can monitor fracture propagation and proppant injection range during fracturing, and determine the effectiveness and extent of the fracturing by combining changes in seismic wave velocity and resistivity. For carbon dioxide geological storage, it can simultaneously monitor… The changes in seismic wave velocity and resistivity in the injection zone are used to accurately depict the changes in seismic wave velocity and resistivity in the injection zone. Gas diffusion boundaries ensure safe storage. In these applications, this method demonstrates significantly superior performance compared to traditional single geophysical methods: it provides a more comprehensive view of parameters, more stable inversion results, and reduces interpretation uncertainties.
[0147] (6) It has good noise resistance and reliability.
[0148] This method employs noise adaptation and multi-constraint balancing strategies, exhibiting good noise resistance and reliability for complex field data, and is able to produce reasonable models even when data has certain missing parts or noise.
[0149] The above description is merely a few embodiments of this application and is not intended to limit this application in any way. Although this application discloses preferred embodiments as described above, it is not intended to limit this application. Any changes or modifications made by those skilled in the art without departing from the scope of the technical solution of this application using the disclosed technical content are equivalent to equivalent implementation cases and fall within the scope of the technical solution.
Claims
1. A time-dimensional seismic wave velocity-resistivity joint inversion method, characterized in that, The method includes: S1. Construct a seismic wave velocity distribution model based on the seismic exploration data of the target area, and construct a resistivity distribution model based on the electromagnetic exploration data of the target area; S2. Construct constraints based on the seismic wave velocity distribution model and the resistivity distribution model, and determine the objective function for the joint inversion of seismic wave velocity and resistivity based on the constraints; the constraints include prediction bias fitting terms, model structure coupling terms, spatial regularization terms, temporal difference constraints, and model boundary consistency constraints. S3. Based on the seismic exploration data and electromagnetic exploration data from multiple periods, update and iterate the seismic wave velocity distribution model and the resistivity distribution model until the objective function is minimized; S4. Based on the seismic wave velocity distribution model and resistivity distribution model corresponding to the minimization of the objective function, determine the seismic wave velocity inversion results and resistivity inversion results of the target area in the multiple periods, respectively. The objective function is expressed as: (1) in, Represents a model of seismic wave velocity distribution; Represents a resistivity distribution model; Represent the objective function; This represents the prediction bias fitting term; Indicates the coupling terms in the model structure; Represents the space regularization term; Indicates the time difference constraint term; This represents the model boundary consistency constraint term; , , and These represent the weight coefficients of the model structure coupling term, spatial regularization term, temporal difference constraint term, and model boundary consistency constraint term, respectively. The boundary consistency constraint term is expressed as follows: （2） in, This represents the model boundary consistency constraint term; Seismic wave velocity distribution model exist Seismic wave velocity at the location; A threshold representing seismic wave velocity; Representing the resistivity distribution model exist resistivity at that point; The threshold value representing resistivity; Represents the unit step function; Seismic wave velocity distribution model exist Binary representation of the location; Representing the resistivity distribution model exist The binary representation of Ω; Ω represents the integration domain, i.e., the entire model space.
2. The method according to claim 1, characterized in that, The seismic exploration data includes the observed first arrival and travel times of seismic waves at multiple seismic measuring points, and the electromagnetic exploration data includes the observed apparent resistivity of the subsurface medium at multiple electromagnetic measuring points.
3. The method according to claim 2, characterized in that, S2 constructs a prediction deviation fitting term based on the seismic wave velocity distribution model and the resistivity distribution model, specifically including: The calculated first arrival travel time at the multiple seismic measuring points is calculated based on the seismic wave velocity distribution model, and the calculated apparent resistance at the multiple electromagnetic measuring points is calculated based on the resistivity distribution model. Based on the calculated first arrival time, the calculated apparent resistivity, the observed first arrival time, and the observed apparent resistivity, a prediction bias fitting term is constructed.
4. The method according to claim 1, characterized in that, S2 constructs model structure coupling terms based on the seismic wave velocity distribution model and the resistivity distribution model, specifically including: Based on the gradient vectors at different locations in the seismic wave velocity distribution model and the resistivity distribution model, a model structure coupling term is constructed using the model coupling method.
5. The method according to claim 1, characterized in that, S2 constructs a spatial regularization term based on the seismic wave velocity distribution model and the resistivity distribution model, specifically including: Based on the gradient vectors at different locations in the seismic wave velocity distribution model and the resistivity distribution model, a spatial regularization term is constructed using the smoothness constraint method.
6. The method according to claim 1, characterized in that, S2 constructs time difference constraint terms based on the seismic wave velocity distribution model and the resistivity distribution model, specifically including: The seismic wave velocity distribution of the target area during the multiple periods is calculated based on the seismic wave velocity distribution model, and the resistivity distribution of the target area during the multiple periods is calculated based on the resistivity distribution model. Based on the seismic wave velocity and resistivity distributions corresponding to the multiple periods, a time difference constraint term is constructed.
7. The method according to claim 1, characterized in that, S3 involves updating and iterating the seismic wave velocity distribution model and the resistivity distribution model, specifically including: The seismic wave velocity distribution model and the resistivity distribution model are updated synchronously in each iteration; Alternatively, the seismic wave velocity distribution model and the resistivity distribution model can be updated alternately in multiple iterations.
8. The method according to claim 7, characterized in that, In each iteration, the seismic wave velocity distribution model and the resistivity distribution model are updated synchronously, specifically including: In each iteration, the model parameters of the seismic wave velocity distribution model and the resistivity distribution model for each period are updated synchronously. Alternatively, in one iteration, only the model parameters of the seismic wave velocity distribution model and the resistivity distribution model for one period are updated, and in the next iteration, the model parameters of the seismic wave velocity distribution model and the resistivity distribution model for the next period are updated based on the model parameters corresponding to the previous period.