A method and device for monitoring surface deformation at a CO2 injection site

Abstract

The application provides a CO2 injection site surface deformation monitoring method and device. The method obtains deformation time series monitored by each GNSS station in the target site within a preset time period and a plurality of InSAR deformation interferograms, and maps to a unified space-time grid to obtain GNSS deformation values and InSAR deformation values of each space-time grid point; for the deformation time series corresponding to the GNSS deformation value of each space-time grid point and the space field corresponding to the InSAR deformation value of each time point, the high-frequency and low-frequency time components and the high-frequency and low-frequency space components are obtained, so as to obtain a space-time high-frequency fusion field and a space-time low-frequency fusion field; based on the space-time high-frequency fusion field and the space-time low-frequency fusion field of each space-time grid point, the fusion deformation value of each space-time grid point is determined to obtain the continuous slope distribution of the surface at each time point. The method can obtain a high-precision and high-continuity fusion deformation field and surface slope distribution.

Description

Technical Field

[0001] This application relates to the field of carbon capture, utilization and storage (CCUS) monitoring technology, and more specifically, to a method and apparatus for monitoring surface deformation at a CO2 injection site. Background Technology

[0002] Currently, monitoring surface deformation at CO2 storage sites primarily relies on single remote sensing or ground-based measurement methods. Among these, spaceborne synthetic aperture radar interferometry (InSAR) can acquire large-area, high spatial resolution relative deformation fields. Through differential processing of radar images at different time phases, it can identify uplift or subsidence areas around the injection well at the centimeter level. Global Navigation Satellite Systems (GNSS), on the other hand, employ continuously operating fixed receiving stations to acquire high temporal resolution three-dimensional absolute displacement sequences for single points at millisecond-level sampling rates. GNSS is often used to calibrate the absolute benchmark of InSAR results and capture transient deformation events.

[0003] However, during CO2 injection through buried pipes, existing monitoring methods cannot obtain surface deformation fields with both high temporal and spatial resolutions due to the low temporal resolution, omission of plume transport events during revisit periods, and interference from atmospheric noise with millimeter-level uplift signals by the single InSAR technology. At the same time, the single GNSS technology suffers from sparse spatial deployment and the inability to reconstruct the non-uniform uplift field morphology formed by plume transport. Summary of the Invention

[0004] The purpose of this application is to provide a method and apparatus for monitoring surface deformation at CO2 injection sites. By mapping GNSS and InSAR data to a unified spatiotemporal grid and performing multi-scale decomposition and frequency fusion, it can simultaneously utilize the high temporal resolution of GNSS and the high spatial resolution of InSAR to obtain a high-precision and highly continuous fused deformation field and surface slope distribution, thereby significantly improving the spatiotemporal resolution capability of CO2 injection site caprock deformation monitoring and the reliability of risk assessment.

[0005] Firstly, a method for monitoring surface deformation at a CO2 injection site is provided, which may include: The deformation time series and multiple InSAR deformation interferograms monitored by each GNSS station at the target site within a preset time period are obtained. The data of the two are mapped to a unified spatiotemporal grid to obtain the GNSS deformation value and InSAR deformation value of each spatiotemporal grid point. The deformation time series corresponding to the GNSS deformation values ​​of each spatial grid point is decomposed in time dimension to obtain the high-frequency and low-frequency time components; the spatial field corresponding to the InSAR deformation values ​​of each time point is decomposed in space dimension to obtain the high-frequency and low-frequency spatial components. Based on high-frequency time components and high-frequency spatial components, a spatiotemporal high-frequency fusion field is obtained; based on low-frequency time components and low-frequency spatial components, a spatiotemporal low-frequency fusion field is obtained. Based on the spatiotemporal high-frequency fusion field and spatiotemporal low-frequency fusion field of each spatiotemporal grid point, the fusion deformation value of each spatiotemporal grid point is determined to obtain the continuous slope distribution of the land surface at each time point.

[0006] In one possible implementation, the acquisition of high-frequency and low-frequency time components includes: For each spatial grid point, its GNSS deformation values ​​at all time points are arranged in chronological order to form a one-dimensional deformation time series; A multi-resolution analysis algorithm is used to decompose the deformation time series, and the approximation coefficients and detail coefficients are separated layer by layer according to the preset number of decomposition layers. The component corresponding to the last layer of approximation coefficients is taken as the low-frequency time component, and the component after superimposing all layers of detail coefficients is taken as the high-frequency time component.

[0007] In one possible implementation, the acquisition of high-frequency spatial components and low-frequency spatial components includes: For each time point, the InSAR deformation values ​​of all spatial grid points at that time point are arranged in the order of the rows of the spatial grid points to form a spatial field, which is a two-dimensional discrete signal matrix. A two-dimensional discrete wavelet transform is performed on the two-dimensional discrete signal matrix, and the low-frequency subband coefficients and the high-frequency subband coefficients in three directions are obtained by decomposing the matrix layer by layer. The spatial field obtained by reconstructing the low-frequency subband coefficients of the final layer is taken as the low-frequency spatial component, and the spatial field superimposed after reconstructing all the high-frequency subband coefficients is taken as the high-frequency spatial component.

[0008] In one possible implementation, obtaining the spatiotemporal high-frequency fusion field includes: For each spatiotemporal grid point, extract the high-frequency temporal component value and the high-frequency spatial component value of that spatiotemporal grid point; calculate the local variance of the high-frequency temporal component and the local variance of the high-frequency spatial component within the neighborhood window centered on that spatiotemporal grid point; and calculate the fusion weight of that spatiotemporal grid point based on the above two local variances. The high-frequency time component value and the high-frequency spatial component value of the spatiotemporal grid point are weighted and summed according to the fusion weight to obtain the spatiotemporal high-frequency fusion value of the spatiotemporal grid point; the spatiotemporal high-frequency fusion values ​​of all spatiotemporal grid points constitute the spatiotemporal high-frequency fusion field.

[0009] In one possible implementation, obtaining the spatiotemporal low-frequency fused field includes: For each spatiotemporal grid point, the low-frequency time component value and low-frequency spatial component value of that spatiotemporal grid point are extracted; the covariance cross-fusion method is used to calculate the local estimated covariance matrix of the low-frequency time component and the low-frequency spatial component at that spatiotemporal grid point respectively; Based on the local estimated covariance matrix, the optimal fusion coefficient for the spatiotemporal grid point is obtained by solving a system of linear equations that minimize the variance of the fused estimation error. The low-frequency time component value and low-frequency spatial component value of the spatiotemporal grid point are linearly combined according to the optimal fusion coefficient to obtain the spatiotemporal low-frequency fusion value of the spatiotemporal grid point; the spatiotemporal low-frequency fusion values ​​of all spatiotemporal grid points constitute the spatiotemporal low-frequency fusion field.

[0010] In one possible implementation, the fused deformation value for each spatiotemporal grid point is determined, including: For each spatial grid point, the spatiotemporal low-frequency fusion values ​​of the spatial grid point at all time points are arranged in chronological order to serve as the approximation coefficient sequence of the spatial grid point; the spatiotemporal high-frequency fusion values ​​of the spatial grid point at the same time series are used as the detail coefficient sequence of the spatial grid point. Perform a one-dimensional discrete wavelet inverse transform on the approximation coefficient sequence and the detail coefficient sequence to reconstruct the complete deformation value of the spatial grid point at each time point; traverse all spatial grid points to obtain the fused deformation value of each spatiotemporal grid point.

[0011] In one possible implementation, the distribution of continuous surface slope at each time point is obtained, including: For each time point, extract the two-dimensional spatial deformation matrix of that time point from the fused deformation field; Based on this two-dimensional spatial deformation matrix, the deformation gradient of each spatial grid point in the first and second directions is calculated; Based on the deformation gradient in two directions, the composite slope value of each spatial grid point is calculated. This composite slope value represents the degree of inclination of the ground surface at that spatial grid point. After spatial smoothing of the composite slope values ​​of all spatial grid points, the spatial grid points whose smoothed composite slope values ​​exceed the preset safe slope are clustered into connected domains. Areas whose connected domain area exceeds the preset risk area threshold are marked as cap uplift risk areas, and the boundary coordinates and average slope values ​​of the risk areas are output.

[0012] Secondly, a CO2 injection site surface deformation monitoring device is provided, which may include: The acquisition unit is used to acquire the deformation time series and multiple InSAR deformation interferograms monitored by each GNSS station at the target site within a preset time period, and to map the data of the two to a unified spatiotemporal grid to obtain the GNSS deformation value and InSAR deformation value of each spatiotemporal grid point. The decomposition unit is used to perform time-dimensional decomposition on the deformation time series corresponding to the GNSS deformation values ​​of each spatial grid point to obtain the high-frequency and low-frequency time components; and to perform spatial-dimensional decomposition on the spatial field corresponding to the InSAR deformation values ​​of each time point to obtain the high-frequency and low-frequency spatial components. The acquisition unit is also used to acquire a spatiotemporal high-frequency fusion field based on high-frequency time components and high-frequency spatial components; and to acquire a spatiotemporal low-frequency fusion field based on low-frequency time components and low-frequency spatial components. The determination unit is used to determine the fusion deformation value of each spatiotemporal grid point based on the spatiotemporal high-frequency fusion field and the spatiotemporal low-frequency fusion field of each spatiotemporal grid point, so as to obtain the continuous slope distribution of the land surface at each time point.

[0013] Thirdly, an electronic device is provided, which includes a processor, a communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other through the communication bus; Memory, used to store computer programs; When a processor executes a program stored in memory, it implements any of the steps described in the first aspect above.

[0014] Fourthly, a computer-readable storage medium is provided, wherein a computer program is stored therein, and when executed by a processor, the computer program implements the steps of any of the methods described in the first aspect above.

[0015] This application provides a method and apparatus for monitoring surface deformation at a CO2 injection site. The method acquires deformation time series and multiple InSAR deformation interferograms monitored by GNSS stations at the target site within a preset time period. The data is mapped to a unified spatiotemporal grid to obtain GNSS deformation values ​​and InSAR deformation values ​​at each spatiotemporal grid point. The deformation time series corresponding to the GNSS deformation values ​​at each spatial grid point is decomposed into high-frequency and low-frequency time components. The spatial field corresponding to the InSAR deformation values ​​at each time point is decomposed into high-frequency and low-frequency spatial components. Based on the high-frequency time and high-frequency spatial components, a spatiotemporal high-frequency fused field is obtained. Based on the low-frequency time and low-frequency spatial components, a spatiotemporal low-frequency fused field is obtained. Based on the spatiotemporal high-frequency and spatiotemporal low-frequency fused fields at each spatiotemporal grid point, the fused deformation value at each spatiotemporal grid point is determined to obtain the continuous slope distribution of the surface at each time point. This method maps GNSS and InSAR data to a unified spatiotemporal grid and performs multi-scale decomposition and frequency fusion. It can simultaneously utilize the high temporal resolution of GNSS and the high spatial resolution of InSAR to obtain a high-precision and highly continuous fused deformation field and surface slope distribution, thereby significantly improving the spatiotemporal resolution of CO2 injection site caprock deformation monitoring and the reliability of risk assessment. Attached Figure Description

[0016] To more clearly illustrate the technical solutions of the embodiments of this application, the accompanying drawings used in the embodiments of this application will be briefly introduced below. It should be understood that the following drawings only show some embodiments of this application and should not be regarded as a limitation of the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0017] Figure 1 A flowchart illustrating a method for monitoring surface deformation at a CO2 injection site, provided in an embodiment of this application; Figure 2 A schematic diagram of a surface deformation monitoring device for a CO2 injection site provided in this application embodiment; Figure 3 This is a schematic diagram of the structure of an electronic device provided in an embodiment of this application. Detailed Implementation

[0018] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of the embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application. Unless otherwise defined, the technical or scientific terms used in this application should have the ordinary meaning understood by those skilled in the art. The words "first," "second," and similar terms used in this application do not indicate any order, quantity, or importance, but are only used to distinguish different components. The words "comprising" or "including," etc., mean that the element or object preceding the word covers the element or object listed after the word and its equivalents, but do not exclude other elements or objects. The words "connected," "coupled," or "connected," etc., are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. "Up," "down," "left," "right," etc., are only used to indicate relative positional relationships. When the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0019] The surface deformation monitoring method for CO2 injection sites provided in this application aims to address the shortcomings of single monitoring methods (such as using only GNSS or only InSAR) in terms of spatiotemporal resolution, continuity, and accuracy, and overcome the limitations of traditional fusion methods in dealing with non-stationary, multi-scale deformation fields. This method, through multi-scale decomposition, adaptive weighted fusion, and dynamic decomposition strategies, can reconstruct the spatiotemporal evolution of surface deformation caused by CO2 injection with high accuracy and reliability, and intelligently identify caprock uplift risk areas based on the reconstructed fused deformation field.

[0020] To achieve the above objectives, the technical solution adopted in this application includes a series of carefully designed data processing and fusion steps. In a typical implementation environment, this method relies on a network of continuously operating Global Navigation Satellite System (GNSS) reference stations distributed around the CO2 injection site and its surroundings, as well as ascending and / or descending orbit imagery data from spaceborne Synthetic Aperture Radar (InSAR). GNSS stations can provide high temporal resolution (e.g., hourly or daily) but relatively sparse spatially distributed deformation time series; InSAR can provide high spatial resolution (e.g., tens of meters per pixel) but non-uniform temporal sampling (depending on the satellite revisit period) deformation interferograms. The two types of data complement each other through the technical solution of this application.

[0021] The preferred embodiments of this application are described below with reference to the accompanying drawings. It should be understood that the preferred embodiments described herein are for illustration and explanation only and are not intended to limit this application. Furthermore, the embodiments and features in the embodiments of this application can be combined with each other without conflict.

[0022] Figure 1 This is a flowchart illustrating a method for monitoring surface deformation at a CO2 injection site, provided as an embodiment of this application. Figure 1 As shown, the method may include: Step S110: Obtain the deformation time series and multiple InSAR deformation interferograms monitored by each GNSS station at the target site within a preset time period, map the data of the two to a unified spatiotemporal grid, and obtain the GNSS deformation value and InSAR deformation value of each spatiotemporal grid point.

[0023] Obtain available data for the target CO2 injection site area within a preset monitoring period (e.g., from six months before injection begins to one year after injection ceases). Available data includes: Time series data of surface deformation monitored by each GNSS station. The data for each station may include a timestamp t and the three-dimensional (east E, north N, vertical U) or at least vertical deformation value at that time.

[0024] Multiple InSAR deformation interferograms. Each interferogram corresponds to the spatial deformation field at a specific time baseline (e.g., from the main image time t_m to the secondary image time t_s), and is usually in the form of geocoded raster images.

[0025] Subsequently, the two types of data are resampled and mapped to a unified spatiotemporal grid (spatiotemporal reference). Specifically, a unified regular spatial grid (e.g., with a resolution consistent with or slightly lower than that of the InSAR data, such as 20m × 20m) and a unified time axis (e.g., in days, covering the entire monitoring period) are determined. For each spatiotemporal grid point (x, y, t): GNSS deformation value: estimated from sparse GNSS station observations through spatiotemporal interpolation (such as ordinary kriging or spatiotemporal inverse distance weighting), denoted as G(x,y,t).

[0026] InSAR Deformation Values: For each interferogram, the deformation values ​​within its covered time period are distributed to various time points on a unified time axis through either time-linear interpolation or time-weighted interpolation based on a deformation model. Spatially, the interferogram pixel values ​​are directly sampled onto a unified spatial grid using bilinear interpolation, denoted as I(x,y,t).

[0027] Thus, the registered GNSS deformation value G(x,y,t) and InSAR deformation value I(x,y,t) are obtained at each spatiotemporal grid point.

[0028] In a preferred embodiment, prior to resampling, to improve data quality, this application may further include a step of anomaly detection and repair of the collected raw data, including: For GNSS time series: For each GNSS station, a sliding window (e.g., a window width of 30 days) is used. At each window location, the local mean μ_local and local standard deviation σ_local of the deformation values ​​within the window are calculated. Data points deviating from the local mean by more than three times the local standard deviation (i.e., |value-μ_local|>3σ_local) are marked as outliers. For the marked outlier moments, spatiotemporal kriging interpolation is used for repair. Specifically, using the deformation values ​​of the station at adjacent moments (e.g., three moments before and after), and the deformation values ​​of neighboring stations at the same moment (e.g., within a 5-kilometer radius), a spatiotemporal variogram model is constructed to perform optimal unbiased estimation of the outlier points, replacing the original outliers.

[0029] For InSAR interferograms: Isolated erroneous pixel blocks often appear in InSAR data due to phase unwrapping errors. This method identifies these errors through connected component analysis: pixel blocks with abrupt changes in pixel value (difference from surrounding pixels exceeding an integer multiple of 2π) and an area smaller than a preset spatial connected component threshold (e.g., less than 50 consecutive pixels) are deleted entirely. The resulting holes are filled using Gaussian process regression with surrounding valid pixels (all valid points within 5 pixels of the hole boundary). Gaussian process regression can provide an estimate of the uncertainty of the filled value, thus ensuring the statistical reasonableness of the filling result.

[0030] Step S120: Perform time dimension decomposition on the deformation time series corresponding to the GNSS deformation values ​​of each spatial grid point to obtain the high-frequency and low-frequency time components; perform spatial dimension decomposition on the spatial field corresponding to the InSAR deformation values ​​of each time point to obtain the high-frequency and low-frequency spatial components.

[0031] This step requires spatiotemporal multi-scale decomposition of the resampled data, with the aim of decomposing the complex deformation signal into different components with more physical meaning.

[0032] (1) Time dimension decomposition: For each fixed spatial grid point (i.e., fixed planar position coordinates (x, y)) in the unified spatiotemporal grid, perform the following operations: The GNSS deformation values ​​(i.e., the obtained G(x,y,t)) of the spatial grid point at all time points t_1, t_2, ..., t_T within a preset time period are arranged in chronological order to form a one-dimensional array. This array is the deformation time series of the spatial grid point, denoted as S_G(t), where t is the discrete time index.

[0033] A one-dimensional discrete wavelet basis (such as the Daubechies wavelet db4) and a preset number of decomposition levels L (e.g., L=3) are selected. The Mallat pyramid algorithm is executed on the deformed time series S_G(t) to decompose it level by level. In each level of decomposition, the approximation coefficients of the previous level are further decomposed into new approximation coefficients (low-frequency part) and detail coefficients (high-frequency part) for the current level. Specifically: Level 1 decomposition: S_G(t) → approximation coefficient cA1 + detail coefficient cD1.

[0034] Level 2 decomposition: cA1 → approximation coefficient cA2 + detail coefficient cD2.

[0035] ... Lth level decomposition: cA(L-1) → approximation coefficient cAL + detail coefficient cDL.

[0036] Next, the approximation coefficients cAL of the last layer (the Lth layer) are taken and reconstructed using one-dimensional discrete wavelet reconstruction to obtain a reconstructed sequence with the same length as the original time series. This sequence is the low-frequency time component of the spatial grid points, denoted as the low-frequency time component G_low(x,y,t). This component reflects the long-term trend and regional deformation caused by CO2 injection, such as the slow surface uplift caused by the overall increase in pore pressure.

[0037] Simultaneously, the detail coefficients (cD1, cD2, ..., cDL) of all layers (layer 1 to layer L) are reconstructed using one-dimensional discrete wavelet reconstruction. Then, the reconstructed sequences from each layer are superimposed point-by-point to obtain a new sequence of the same length as the original time series. This sequence represents the high-frequency time component of the spatial grid points, denoted as the high-frequency time component G_high(x,y,t). This component captures useful details in short-term fluctuations of deformation, transient disturbances, or monitoring noise, such as the rapid elastic response caused by changes in injection and sampling pulse pressure.

[0038] By traversing all spatial grid points, the GNSS low-frequency time component G_low(x,y,t) and high-frequency time component G_high(x,y,t) corresponding to each spatiotemporal grid point can be obtained.

[0039] (2) Spatial dimension decomposition: For each fixed time point t (i.e., a certain time layer in the unified spatiotemporal grid), perform the following operations: Constructing the spatial field: Extract the InSAR deformation values ​​(i.e., I(x,y,t)) of all spatial grid points (planar positions (x,y)) at time point t. Arrange these values ​​in spatial row order (e.g., the x-direction corresponds to the column index) to form a two-dimensional discrete signal matrix, denoted as M_I(t). This matrix is ​​the spatial field at that time point.

[0040] A two-dimensional discrete wavelet basis (such as the two-dimensional tensor product wavelet db4 corresponding to the one-dimensional wavelet) and a preset number of decomposition levels M (e.g., M=3) are selected. The two-dimensional Mallat decomposition algorithm is performed on the two-dimensional discrete signal matrix M_I(t). In each decomposition level, the low-frequency subband coefficients of the previous level are further decomposed into four subbands: A low-frequency subband coefficient (approximate component, denoted as cA_j).

[0041] The high-frequency subband coefficients (detail components) are in three directions: horizontal (cH_j), vertical (cV_j), and diagonal (cD_j). Here, j represents the decomposition layer number.

[0042] Next, the low-frequency subband coefficients cA_M of the final layer (layer M) are taken. Two-dimensional discrete wavelet reconstruction (inverse transform) is performed on these coefficients to obtain a reconstruction matrix with the same size as the original spatial field matrix. This matrix represents the low-frequency spatial component at time point t, denoted as the low-frequency spatial component I_low(x,y,t). This component characterizes the regional background trend and large-scale smooth changes of the deformation field.

[0043] Simultaneously, high-frequency subband coefficients (i.e., cH_1, cV_1, cD_1, cH_2, ..., cD_M) of all decomposition layers (layers 1 to M) and all directions (horizontal, vertical, and diagonal) are collected. Two-dimensional discrete wavelet reconstruction is performed on each high-frequency subband coefficient, and then all reconstructed matrices are superimposed pixel-by-pixel to form a new matrix of the same size as the original spatial field. This matrix represents the high-frequency spatial component at time point t, denoted as the high-frequency spatial component I_high(x,y,t). This component captures the edges of the deformation field, local anomalies, the fine structure of the plume front, and spatial discontinuities.

[0044] By iterating through all time points, we can obtain the InSAR low-frequency spatial component I_low(x,y,t) and high-frequency spatial component I_high(x,y,t) corresponding to each spatiotemporal grid point.

[0045] In a preferred embodiment that embodies dynamic adaptability, the number of decomposition layers L and M are not fixed, but are adaptively determined based on the dynamic characteristics of CO2 injection.

[0046] Adaptive Time-Dimension Decomposition Layer L: Obtain the cumulative injection volume curve of the CO2 injection well over time. Calculate the first derivative of this curve to obtain the instantaneous injection rate r(t). Analyze the fluctuation characteristics of r(t). Specifically, monitor the fluctuation period of r(t) (e.g., find the period corresponding to the dominant frequency by calculating its power spectral density). If the monitored fluctuation period is less than a preset time threshold (e.g., less than 7 days indicates the presence of rapid injection pulses or pressure oscillations), increase the time-dimensional decomposition layer L (e.g., from 3 layers to 5 layers). Increasing the number of layers allows for the decomposition of higher-order high-frequency components, thereby more finely characterizing the instantaneous elastic deformation response caused by rapid pressure changes.

[0047] Adaptation of the number of spatial dimension decomposition layers M: Calculate the average spatial gradient of the InSAR deformed spatial field at each time point. I(x,y,t). This gradient value reflects the degree of drastic local changes in the deformation field. If the average spatial gradient is greater than a preset gradient threshold (e.g., greater than 5 mm / 100 m), it indicates that complex fine structures (such as fingering) may have formed at the edge of the CO2 plume or that local folding has occurred in the caprock. In this case, the number of spatial decomposition layers M is increased (e.g., from 3 to 4 or 5 layers) to capture more high-frequency spatial information at more scales, thereby more accurately characterizing the fine morphology of the plume leading edge.

[0048] Step S130: Obtain the spatiotemporal high-frequency fusion field based on the high-frequency time component and the high-frequency spatial component; obtain the spatiotemporal low-frequency fusion field based on the low-frequency time component and the low-frequency spatial component.

[0049] A. Constructing a spatiotemporal high-frequency fusion field (high-frequency component fusion): 1. Component extraction: For each spatiotemporal grid point (x,y,t), extract its high-frequency time component value G_high(x,y,t) and high-frequency spatial component value I_high(x,y,t).

[0050] 2. Calculate the fusion weights: First, define a spatial neighborhood window W centered on the current grid point (x,y,t) (e.g., 11×11 pixels; in the time dimension, it can also be formed by taking one time point before and after the current grid point to create a spatiotemporal window). Within this window: Calculate the local variance of the GNSS high-frequency time component σ_G_high²=Var(G_high(p)), where p∈W.

[0051] Calculate the local variance of the InSAR high-frequency spatial components σ_I_high²=Var(I_high(p)), where p∈W.

[0052] Components with high local variance indicate drastic changes within that local region, containing richer and more reliable high-frequency details, and therefore should be assigned higher weights. Based on this principle, the base weights are calculated as follows: w_G_raw=σ_G_high² / (σ_G_high²+σ_I_high²+ε); w_I_raw=σ_I_high² / (σ_G_high²+σ_I_high²+ε); Here, ε is a very small positive number (such as 1e-8) used to prevent the denominator from being zero.

[0053] Furthermore, the calculation of the fusion weights can also incorporate real-time flow data from the CO2 injection well as an external modulation factor to enhance the physical response characteristics of the method, including: Obtain the instantaneous CO2 injection rate v(t) at time t corresponding to the current spatiotemporal grid point. Set an injection rate threshold T_v (e.g., 30% of the maximum design rate of the injection well).

[0054] When v(t) > T_v (high injection rate stage, deformation response is mainly characterized by rapid time changes): increase the global baseline value of the GNSS high-frequency time component. Specifically, set a time advantage factor α_t > 1 (e.g., α_t = 1.5). The updated weights are: w_G_final=(α_t*w_G_raw) / (α_t*w_G_raw+w_I_raw); w_I_final=w_I_raw / (α_t*w_G_raw+w_I_raw); This method greatly enhances the time sensitivity of the fusion field to the response of the fast feather front.

[0055] When v(t) ≤ T_v (low injection or cessation of injection phase, deformation is mainly characterized by slow diffusion and spatial reshaping): increase the global baseline value of the InSAR high-frequency spatial components and set the spatial dominance factor β_s > 1 (e.g., β_s = 1.5). The updated weights are: w_G_final=w_G_raw / (w_G_raw+β_s*w_I_raw); w_I_final=(β_s*w_I_raw) / (w_G_raw+β_s*w_I_raw); This approach enhances the ability to depict the details of spatial deformation during slow diffusion.

[0056] 3. Weighted Summation: The high-frequency component values ​​are weighted and summed according to the final weights to obtain the spatiotemporal high-frequency fusion value F_high(x,y,t): F_high(x,y,t)=w_G_final*G_high(x,y,t)+w_I_final*I_high(x,y,t); By traversing all spatiotemporal grid points, the set of all F_high(x,y,t) constitutes the spatiotemporal high-frequency fusion field.

[0057] B. Constructing a spatiotemporal low-frequency fusion field (low-frequency component fusion): The low-frequency components represent the main trend of deformation and the regional background field, and their fusion focuses more on achieving unbiased optimal estimation statistically.

[0058] 1. Component extraction: For each spatiotemporal grid point (x,y,t), extract its low-frequency time component value G_low(x,y,t) and low-frequency spatial component value I_low(x,y,t).

[0059] 2. Calculate the optimal fusion coefficients: Employ the Covariance Intersection (CI) fusion method. This method does not require knowledge of the cross covariance between the two components and is suitable for estimations from different sources (GNSS and InSAR) with unknown correlations.

[0060] First, estimate the local covariance of each component at the spatiotemporal grid point. This can be approximated by analyzing the fluctuations of G_low and I_low within a neighborhood window W (similar to high-frequency fusion but with a larger window, such as 21×21 pixels): Σ_G=Var(G_low(p)) For p∈W, it is denoted as a scalar value σ_G_low² (in practice, isotropic properties can be assumed).

[0061] Σ_I=Var(I_low(p)) For p∈W, it is denoted as σ_I_low².

[0062] According to CI fusion theory, the estimated covariance after fusion is Σ_F=[ωΣ_G] - ¹+(1-ω)Σ_I - ¹] - ¹. The optimal fusion coefficient ω is obtained by minimizing the trace of Σ_F (i.e., minimizing the total variance of the fusion estimate). That is, solving for: ω*=argmin([ωσ_G_low - ²+(1-ω)σ_I_low - ²] - ¹); where 0≤ω*≤1.

[0063] 3. Linear combination: The low-frequency component values ​​are linearly combined according to the optimal fusion coefficient to obtain the spatiotemporal low-frequency fusion value F_low(x,y,t): F_low(x,y,t)=ω*×G_low(x,y,t)+(1-ω*)×I_low(x,y,t); By traversing all spatiotemporal grid points, the set of all F_low(x,y,t) constitutes the spatiotemporal low-frequency fusion field.

[0064] Step S140: Determine the fusion deformation value of each spatiotemporal grid point based on the spatiotemporal high-frequency fusion field and the spatiotemporal low-frequency fusion field of each spatiotemporal grid point.

[0065] For each fixed spatial grid point (x, y), perform the following steps, including: The spatiotemporal low-frequency fusion values ​​F_low(x,y,t) of the spatial grid point at all time points are arranged in time order t=1,2,...,T, and used as the approximate coefficient sequence A(x,y,t) for one-dimensional discrete wavelet inverse reconstruction of the spatial grid point.

[0066] The spatiotemporal high-frequency fusion value F_high(x,y,t) of the same spatial grid point in the same time series is used as the detail coefficient sequence D(x,y,t) of that spatial grid point.

[0067] For the approximation coefficient sequence A(x,y,t) and the detail coefficient sequence D(x,y,t), a one-dimensional discrete wavelet inverse transform (IDWT) is performed using the wavelet basis and level L used in the time-dimensional decomposition described above. This inverse transform reconstructs the low-frequency approximation and high-frequency details layer by layer, ultimately restoring the complete deformation signal at each time point.

[0068] Output fused deformation value: The result of the inverse transformation is the complete fused deformation value F_final(x,y,t) of the spatial grid point at each time point. By traversing all spatial grid points, the fused deformation value of each spatiotemporal grid point (each spatial point, each time point) can be obtained. The set of all F_final(x,y,t) constitutes the final fused deformation field. This deformation field combines the high temporal resolution of GNSS and the high spatial resolution of InSAR, while reducing noise and system errors through multi-scale decomposition and adaptive fusion.

[0069] Step S150: Based on the fused deformation field formed by the fused deformation values ​​of each spatiotemporal grid point, obtain the continuous slope distribution of the land surface at each time point.

[0070] Based on the reconstructed high-quality fused deformation field, the impact of surface deformation on the integrity of the caprock is further assessed, and potential risk areas are identified.

[0071] For each time point t, perform the following steps, including: Extract the two-dimensional spatial deformation matrix F_final(x,y,t) at this time point from the fused deformation field.

[0072] Based on this matrix, calculate the deformation gradient of each spatial grid point (x, y) in the first direction (e.g., east-west) and the second direction (e.g., north-south). The gradient can be calculated using the central difference method. The east-west gradient Grad_E(x,y)=[F_final(x+1,y)-F_final(x-1,y)] / (2*Δx); The north-south gradient Grad_N(x,y)=[F_final(x,y+1)-F_final(x,y-1)] / (2*Δy); Where Δx and Δy are the physical spacing of the spatial grid in the east-west and north-south directions, respectively.

[0073] The composite slope value Slope(x,y) represents the degree of inclination of the ground surface at that grid point, and is synthesized from the gradients in two directions: Slope(x,y)=arctan(sqrt(Grad_E(x,y)²+Grad_N(x,y)²)) This value is usually expressed as a decimal for the slope of the ground, or converted to a slope angle (degrees).

[0074] The slope map obtained by direct calculation may contain isolated slope spikes caused by residual noise from fusion. Spatial smoothing can be performed using median filtering or Gaussian filtering to effectively preserve the edges of the main deformation gradients while eliminating isolated abnormal spikes.

[0075] Next, risk areas can be identified, including: Spatial grid points whose smoothed synthetic slope values ​​exceed the preset safe slope threshold (e.g., arctan(0.01), corresponding to a slope angle of approximately 0.57°; the specific value needs to be set according to the lithology of the strata, mechanical parameters of the caprock, etc.) are marked as potential high-risk points.

[0076] Perform connectivity analysis on all high-risk points and cluster spatially adjacent points (e.g., using the 8-neighbor rule) into a risk region.

[0077] For each connected component (risk area), calculate its area (number of pixels × area of ​​a single pixel). Areas with a connected component area smaller than a preset risk area threshold (e.g., 10,000 square meters) are removed, as they are considered not to constitute a risk of overall caprock damage.

[0078] The reserved areas (connected regions with an area greater than a threshold) are marked as caprock uplift risk zones. The boundary coordinates (such as the set of polygonal points constituting its outline) and the average slope value within each risk zone are output. This information can provide direct basis for dynamic adjustments to CO2 injection schemes, triggering of early warning systems, and site safety assessments.

[0079] Corresponding to the above method, this application also provides a CO2 injection site surface deformation monitoring device, such as... Figure 2 As shown, the device includes: The acquisition unit 210 is used to acquire the deformation time series and multiple InSAR deformation interferograms monitored by each GNSS station at the target site within a preset time period, and to map the data of the two to a unified spatiotemporal grid to obtain the GNSS deformation value and InSAR deformation value of each spatiotemporal grid point. Decomposition unit 220 is used to perform time dimension decomposition on the deformation time series corresponding to the GNSS deformation values ​​of each spatial grid point to obtain the high-frequency and low-frequency time components; and to perform spatial dimension decomposition on the spatial field corresponding to the InSAR deformation values ​​of each time point to obtain the high-frequency and low-frequency spatial components. The acquisition unit 210 is also used to acquire a spatiotemporal high-frequency fusion field based on high-frequency time components and high-frequency spatial components; and to acquire a spatiotemporal low-frequency fusion field based on low-frequency time components and low-frequency spatial components. The determination unit 230 is used to determine the fusion deformation value of each spatiotemporal grid point based on the spatiotemporal high-frequency fusion field and the spatiotemporal low-frequency fusion field of each spatiotemporal grid point, so as to obtain the continuous slope distribution of the land surface at each time point.

[0080] The functions of each functional unit of the CO2 injection site surface deformation monitoring device provided in the above embodiments of this application can be realized through the above methods and steps. Therefore, the specific working process and beneficial effects of each unit in the CO2 injection site surface deformation monitoring device provided in the embodiments of this application will not be repeated here.

[0081] This application also provides an electronic device, such as... Figure 3 As shown, it includes a processor 310, a communication interface 320, a memory 330, and a communication bus 340, wherein the processor 310, the communication interface 320, and the memory 330 communicate with each other through the communication bus 340.

[0082] Memory 330 is used to store computer programs; When the processor 310 executes the program stored in the memory 330, it performs the following steps: The deformation time series and multiple InSAR deformation interferograms monitored by GNSS stations at the target site within a preset time period are acquired. These data are mapped to a unified spatiotemporal grid to obtain the GNSS deformation value and InSAR deformation value at each spatiotemporal grid point. The deformation time series corresponding to the GNSS deformation value at each spatial grid point are decomposed into high-frequency and low-frequency time components. The spatial field corresponding to the InSAR deformation value at each time point is decomposed into high-frequency and low-frequency spatial components. Based on the high-frequency time and high-frequency spatial components, a spatiotemporal high-frequency fused field is obtained. Based on the low-frequency time and low-frequency spatial components, a spatiotemporal low-frequency fused field is obtained. Based on the spatiotemporal high-frequency and spatiotemporal low-frequency fused fields at each spatiotemporal grid point, the fused deformation value at each spatiotemporal grid point is determined to obtain the continuous slope distribution of the land surface at each time point.

[0083] The communication bus mentioned above can be a Peripheral Component Interconnect (PCI) bus or an Extended Industry Standard Architecture (EISA) bus, etc. This communication bus can be divided into address bus, data bus, control bus, etc. For ease of illustration, only one thick line is used to represent it in the diagram, but this does not mean that there is only one bus or one type of bus.

[0084] The communication interface is used for communication between the aforementioned electronic devices and other devices.

[0085] The memory may include random access memory (RAM) or non-volatile memory (NVM), such as at least one disk storage device. Optionally, the memory may also be at least one storage device located remotely from the aforementioned processor.

[0086] The processors mentioned above can be general-purpose processors, including central processing units (CPUs), network processors (NPs), etc.; they can also be digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components.

[0087] The implementation methods and beneficial effects of the various components of the electronic device in the above embodiments for solving the problem can be found in [reference needed]. Figure 1 The steps in the illustrated embodiments are used to implement the electronic device. Therefore, the specific working process and beneficial effects of the electronic device provided in this application will not be repeated here.

[0088] In another embodiment provided in this application, a computer-readable storage medium is also provided, which stores instructions that, when executed on a computer, cause the computer to perform any of the CO2 injection site surface deformation monitoring methods described in the above embodiments.

[0089] In another embodiment provided in this application, a computer program product containing instructions is also provided, which, when run on a computer, causes the computer to execute any of the CO2 injection site surface deformation monitoring methods described in the above embodiments.

[0090] Those skilled in the art will understand that the embodiments in this application can be provided as methods, systems, or computer program products. Therefore, the embodiments in this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the embodiments in this application can take the form of a computer program product implemented 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.

[0091] This application describes embodiments of methods, apparatus (systems), and computer program products according to embodiments of this application with reference to flowchart illustrations and / or block diagrams. 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 illustrations. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0092] 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.

[0093] 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.

[0094] Although preferred embodiments have been described in this application, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments as well as all changes and modifications falling within the scope of the embodiments of this application.

[0095] Obviously, those skilled in the art can make various modifications and variations to the embodiments of this application without departing from the spirit and scope of the embodiments of this application. Therefore, if these modifications and variations to the embodiments of this application fall within the scope of the claims in this application and their equivalents, then this application also intends to include these modifications and variations.

Claims

1. A method for monitoring surface deformation at a CO2 injection site, characterized in that, The method includes: The deformation time series and multiple InSAR deformation interferograms monitored by each GNSS station at the target site within a preset time period are obtained. The data of the two are mapped to a unified spatiotemporal grid to obtain the GNSS deformation value and InSAR deformation value of each spatiotemporal grid point. The deformation time series corresponding to the GNSS deformation values ​​of each spatial grid point is decomposed in time dimension to obtain the high-frequency and low-frequency time components; the spatial field corresponding to the InSAR deformation values ​​of each time point is decomposed in space dimension to obtain the high-frequency and low-frequency spatial components. Based on high-frequency time components and high-frequency spatial components, a spatiotemporal high-frequency fusion field is obtained; based on low-frequency time components and low-frequency spatial components, a spatiotemporal low-frequency fusion field is obtained. Based on the spatiotemporal high-frequency fusion field and spatiotemporal low-frequency fusion field of each spatiotemporal grid point, the fusion deformation value of each spatiotemporal grid point is determined to obtain the continuous slope distribution of the land surface at each time point.

2. The method as described in claim 1, characterized in that, The acquisition of high-frequency and low-frequency time components includes: For each spatial grid point, its GNSS deformation values ​​at all time points are arranged in chronological order to form a one-dimensional deformation time series; A multi-resolution analysis algorithm is used to decompose the deformation time series, and the approximation coefficients and detail coefficients are separated layer by layer according to the preset number of decomposition layers. The component corresponding to the last layer of approximation coefficients is taken as the low-frequency time component, and the component after superimposing all layers of detail coefficients is taken as the high-frequency time component.

3. The method as described in claim 1, characterized in that, The acquisition of high-frequency and low-frequency spatial components includes: For each time point, the InSAR deformation values ​​of all spatial grid points at that time point are arranged in the order of the rows of the spatial grid points to form a spatial field, which is a two-dimensional discrete signal matrix. A two-dimensional discrete wavelet transform is performed on the two-dimensional discrete signal matrix, and the low-frequency subband coefficients and the high-frequency subband coefficients in three directions are obtained by decomposing the matrix layer by layer. The spatial field obtained by reconstructing the low-frequency subband coefficients of the final layer is taken as the low-frequency spatial component, and the spatial field superimposed after reconstructing all the high-frequency subband coefficients is taken as the high-frequency spatial component.

4. The method as described in claim 1, characterized in that, Acquiring the spatiotemporal high-frequency fusion field includes: For each spatiotemporal grid point, extract the high-frequency temporal component value and the high-frequency spatial component value of that spatiotemporal grid point; calculate the local variance of the high-frequency temporal component and the local variance of the high-frequency spatial component within the neighborhood window centered on that spatiotemporal grid point; and calculate the fusion weight of that spatiotemporal grid point based on the above two local variances. The high-frequency time component value and the high-frequency spatial component value of the spatiotemporal grid point are weighted and summed according to the fusion weight to obtain the spatiotemporal high-frequency fusion value of the spatiotemporal grid point; the spatiotemporal high-frequency fusion values ​​of all spatiotemporal grid points constitute the spatiotemporal high-frequency fusion field.

5. The method as described in claim 4, characterized in that, Acquiring the spatiotemporal low-frequency fusion field includes: For each spatiotemporal grid point, the low-frequency time component value and low-frequency spatial component value of that spatiotemporal grid point are extracted; the covariance cross-fusion method is used to calculate the local estimated covariance matrix of the low-frequency time component and the low-frequency spatial component at that spatiotemporal grid point respectively; Based on the local estimated covariance matrix, the optimal fusion coefficient for the spatiotemporal grid point is obtained by solving a system of linear equations that minimize the variance of the fused estimation error. The low-frequency time component value and low-frequency spatial component value of the spatiotemporal grid point are linearly combined according to the optimal fusion coefficient to obtain the spatiotemporal low-frequency fusion value of the spatiotemporal grid point; the spatiotemporal low-frequency fusion values ​​of all spatiotemporal grid points constitute the spatiotemporal low-frequency fusion field.

6. The method as described in claim 1, characterized in that, Determine the fused deformation value for each spatiotemporal grid point, including: For each spatial grid point, the spatiotemporal low-frequency fusion values ​​of the spatial grid point at all time points are arranged in chronological order to serve as the approximation coefficient sequence of the spatial grid point; the spatiotemporal high-frequency fusion values ​​of the spatial grid point at the same time series are used as the detail coefficient sequence of the spatial grid point. Perform a one-dimensional discrete wavelet inverse transform on the approximation coefficient sequence and the detail coefficient sequence to reconstruct the complete deformation value of the spatial grid point at each time point; traverse all spatial grid points to obtain the fused deformation value of each spatiotemporal grid point.

7. The method as described in claim 1 or 6, characterized in that, Obtain the continuous surface slope distribution at each time point, including: For each time point, extract the two-dimensional spatial deformation matrix of that time point from the fused deformation field; Based on this two-dimensional spatial deformation matrix, the deformation gradient of each spatial grid point in the first and second directions is calculated; Based on the deformation gradient in two directions, the composite slope value of each spatial grid point is calculated. This composite slope value represents the degree of inclination of the ground surface at that spatial grid point. After spatial smoothing of the composite slope values ​​of all spatial grid points, the spatial grid points whose smoothed composite slope values ​​exceed the preset safe slope are clustered into connected domains. Areas whose connected domain area exceeds the preset risk area threshold are marked as cap uplift risk areas, and the boundary coordinates and average slope values ​​of the risk areas are output.

8. A surface deformation monitoring device for a CO2 injection site, characterized in that, The device includes: The acquisition unit is used to acquire the deformation time series and multiple InSAR deformation interferograms monitored by each GNSS station at the target site within a preset time period, and to map the data of the two to a unified spatiotemporal grid to obtain the GNSS deformation value and InSAR deformation value of each spatiotemporal grid point. The decomposition unit is used to perform time-dimensional decomposition on the deformation time series corresponding to the GNSS deformation values ​​of each spatial grid point to obtain the high-frequency and low-frequency time components; and to perform spatial-dimensional decomposition on the spatial field corresponding to the InSAR deformation values ​​of each time point to obtain the high-frequency and low-frequency spatial components. The acquisition unit is also used to acquire a spatiotemporal high-frequency fusion field based on high-frequency time components and high-frequency spatial components; and to acquire a spatiotemporal low-frequency fusion field based on low-frequency time components and low-frequency spatial components. The determination unit is used to determine the fusion deformation value of each spatiotemporal grid point based on the spatiotemporal high-frequency fusion field and the spatiotemporal low-frequency fusion field of each spatiotemporal grid point, so as to obtain the continuous slope distribution of the land surface at each time point.

9. An electronic device, characterized in that, The electronic device includes a processor, a communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other through the communication bus; Memory, used to store computer programs; A processor, when executing a program stored in memory, implements the method of any one of claims 1-7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, implements the method described in any one of claims 1-7.

Images