A water conservancy hub insar deformation time series analysis method and system

Abstract

The present application relates to the technical field of deformation time series analysis, in particular to a water conservancy hub InSAR deformation time series analysis method and system, comprising the following steps: screening radar image of water conservancy hub to construct main and auxiliary data sets, generating differential interferogram, extracting permanent scatterers, inverting deformation velocity and mapping to geographic coordinate system, and outputting InSAR deformation time series analysis result. In the present application, by constructing radar image data set covering key structure area, stripping off terrain and atmospheric interference factors, extracting long-term stable scatterer points and calculating their deformation velocity and elevation error, continuous time series analysis of structure deformation is realized, through multi-source data fusion and phase expression optimization, result precision and anti-interference ability are improved, combined with orbit parameter and atmospheric correction data, high-quality deformation information inversion is completed, identification ability of key structure stability change trend is enhanced, and efficient monitoring of water conservancy hub facilities with wide area coverage, high timeliness and low dependence is realized.

Description

Technical Field

[0001] This invention relates to the field of deformation time series analysis technology, and in particular to an InSAR deformation time series analysis method and system for hydraulic engineering projects. Background Technology

[0002] The field of deformation time series analysis technology involves the systematic monitoring and analysis of the deformation process of the earth's surface or engineering structures over time. It covers key aspects such as multi-period acquisition of deformation information, data preprocessing, time series inversion, and deformation trajectory extraction. It primarily uses remote sensing image processing and time series analysis methods to identify dynamic changes and determine trends in monitored objects, and has wide application value in areas such as geological disaster early warning, water conservancy project supervision, and the safe operation of large-scale infrastructure. Among these, the traditional InSAR deformation time series analysis method for water conservancy projects refers to deformation monitoring based on synthetic aperture radar interferometry on key infrastructure such as pumping stations in water conservancy projects. It targets the identification and assessment of structural stability and deformation development trends during the long-term operation of hydraulic structures. Traditional methods primarily rely on high-precision levels, automatic total stations, global navigation satellite systems (GNSS), and optical remote sensing to acquire structural deformation data. While high-precision levels and automatic total stations offer high measurement accuracy, they typically require manual on-site observation, resulting in high workload and susceptibility to weather conditions. Furthermore, their monitoring range is limited and costs are high. Although GNSS technology offers all-weather operation capabilities and a high level of automation, its elevation accuracy is easily affected by interference in complex environments. While optical remote sensing can acquire large-scale non-contact deformation information at a lower cost and with high spatial resolution, its image acquisition is susceptible to weather factors such as clouds, rain, and fog, leading to unstable image quality and limiting its accuracy and applicability. Therefore, it is difficult to achieve continuous, stable, and efficient acquisition and analysis of large-area, multi-time-period deformation information.

[0003] Existing technologies employ high-precision levels and automatic total stations for deformation monitoring. While these methods offer high measurement accuracy, they rely on manual on-site operation, have long work cycles, and are significantly affected by weather conditions, making it difficult to achieve large-scale continuous observation. Although GNSS offers automation and all-weather capabilities, it is affected by building obstruction and electromagnetic interference in complex environments, resulting in significant errors, especially in elevation measurements. While optical remote sensing can acquire images over a wide area, its image quality is unstable in cloudy or rainy conditions, leading to low data availability. All these methods suffer from limited observation range, high costs, and poor timeliness of deformation information, thus restricting the continuity and stability of deformation monitoring of key structures in water conservancy projects during long-term operation. Summary of the Invention

[0004] The purpose of this invention is to address the shortcomings of existing technologies by proposing an InSAR deformation time series analysis method for water conservancy hubs.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a time-series InSAR deformation analysis method for hydraulic engineering projects, comprising the following steps:

[0006] S1: Select multiple C-band synthetic aperture radar images covering the structural areas of sluice gates, pumping stations, and dams, acquire SAR images, and simultaneously collect SRTM digital elevation model data, POD orbit data, and GACOS atmospheric delay data. Organize the images according to their time labels to construct an interferometric radar dataset.

[0007] S2: Based on the SAR images in the interferometric radar dataset, set the main image and auxiliary image, perform registration operation and terrain phase calculation, and after superimposing the data of each interferometric pair, strip the terrain components and construct the differential interferogram structure matrix.

[0008] S3: Based on the differential interferogram structure matrix, calculate the amplitude deviation index of each pixel, filter permanent scatterer points that are less than the preset deviation threshold, record the spatial index and set it as the interferometric modeling control point, analyze the phase change characteristics of adjacent PS points, and establish the interferometric phase expression matrix.

[0009] S4: Based on the interference phase expression matrix, calculate the deformation velocity and elevation error of each permanent scatterer point, establish a coordinate-velocity mapping, and generate directional deformation velocity data;

[0010] S5: Based on the directional deformation velocity data, perform the LOS velocity to vertical velocity conversion, combine the coordinate mapping parameters in the main image, map to the geographic coordinate system, and construct the InSAR deformation time series analysis results of the water conservancy hub.

[0011] As a further embodiment of the present invention, the interferometric radar dataset includes spatially consistent SAR master data, image time-stamped serialized data, SRTM digital elevation model information, POD precision orbit data, and atmospheric delay preprocessing data. The differential interferogram structure matrix includes terrain stripping interferogram sequences, phase difference calculation units, and structure mapping relationship indexes. The interferometric modeling control points include spatial location indexes, permanent scatterer feature identifiers, and phase change measures. The directional deformation velocity data includes deformation velocity values ​​of permanent scatterer points, elevation error information, and LOS directional velocity mapping relationships. The InSAR deformation time series analysis results of the water conservancy hub include vertical deformation velocity distribution, geographic coordinate system structure label list, and deformation monitoring target index.

[0012] As a further aspect of the present invention, the interferometric radar dataset acquisition step specifically comprises:

[0013] S111: Acquire multiple C-band synthetic aperture radar images covering the sluice gate, pumping station, and dam structure area. Set the imaging mode to interferometric wide-swath mode and the polarization mode to VV mode. Arrange the images sequentially according to the imaging time label of each image and match and filter the spatial coverage range. Remove data frames with spatial intersection less than the preset coverage threshold to obtain the regional matching radar image set.

[0014] S112: Based on the region-matching radar image set, synchronously collect the corresponding SRTM digital elevation model data, POD orbit data and GACOS atmospheric delay data according to the time label, construct the binding relationship between the images and the three types of auxiliary data, and generate an image auxiliary data structure set.

[0015] S113: Based on the image auxiliary data structure set, perform pairing operations on image frames of adjacent imaging times, and combine them with the corresponding elevation, orbit and atmospheric delay data structures to obtain the interferometric radar dataset.

[0016] As a further aspect of the present invention, the step of obtaining the differential interferogram structure matrix specifically comprises:

[0017] S211: Based on the SAR images in the interferometric radar dataset, extract the Doppler frequency center value and spatiotemporal baseline length of each image, compare the baseline length of the main image with the preset baseline threshold, retain the image frames that meet the baseline conditions as the main image, and mark the other image frames as auxiliary images, and generate a main and auxiliary image classification identifier set.

[0018] S212: Based on the main and auxiliary image classification identifier set, perform image registration operation on each pair of main and auxiliary images, adjust the image position relationship based on the pixel grid coordinate difference, perform interferometric superposition operation on the registered images, extract the interferometric phase image sequence, and index and encode it according to the main and auxiliary time labels to obtain the interferometric pair image sequence set;

[0019] S213: Based on the interferometric image sequence set and SRTM digital elevation model data, calculate the terrain phase component in each interferogram, peel off the terrain phase from the interferometric phase map pixel by pixel, and perform sequence aggregation on the peeling results to establish a differential interferogram structure matrix.

[0020] As a further aspect of the present invention, the step of obtaining the interference phase representation matrix specifically comprises:

[0021] S311: Obtain the differential interferogram structure matrix, extract the interference amplitude value sequence of each pixel, perform standard deviation and mean statistics on the sequence, and calculate the amplitude deviation index of each pixel to obtain the amplitude deviation index data.

[0022] S312: Based on the amplitude deviation index data, compare the amplitude deviation index of all pixels with the preset deviation threshold pixel by pixel, filter out pixels that are less than the preset deviation threshold, extract the two-dimensional spatial index coordinates of the corresponding pixels in the interferogram, establish a coordinate structure list, and obtain the permanent scatterer spatial index set.

[0023] S313: Based on the permanent scatterer spatial index set, extract the phase value sequence of all adjacent permanent scatterer points in each interferogram, aggregate the phase difference sequence of adjacent point pairs, construct a two-dimensional matrix structure with rows representing point pairs and columns representing time series, and output the interferometric phase expression matrix.

[0024] As a further aspect of the present invention, the step of acquiring directional deformation speed data specifically comprises:

[0025] ;

[0026] in, This represents the amplitude deviation index of the i-th pixel. This represents the amplitude value of the i-th pixel in the j-th interferogram. Let represent the average amplitude of the i-th pixel across all interferograms, and n represent the total number of interferograms.

[0027] As a further aspect of the present invention, the step of acquiring directional deformation speed data specifically comprises:

[0028] S411: Based on the interferometric phase expression matrix, extract the POD orbital data and GACOS atmospheric delay data, perform synchronization operation according to time label and spatial index, perform vector difference operation on the interferometric phase value and orbital error term, strip the orbital term, perform atmospheric delay error field grid resampling and subtract it from the residual phase matrix to obtain the error stripped phase matrix.

[0029] S412: Based on the error stripping phase matrix, perform least squares linear fitting on the stripping phase of each permanent scatterer point in the time series. In the fitting coefficients, the intercept term represents the elevation error and the slope term represents the rate of change in the time dimension. The execution speed corresponds to the elevation, and obtain the deformation and elevation error estimation table.

[0030] S413: Based on the deformation and elevation error estimation table, extract the two-dimensional coordinate values ​​of each permanent scatterer point and the corresponding deformation rate value to form a coordinate-velocity mapping data frame, and uniformly project the velocity components of each point in the structure to the satellite line of sight to generate directional deformation velocity data.

[0031] As a further aspect of the present invention, the steps for obtaining the InSAR deformation time series analysis results of the water conservancy hub are as follows:

[0032] S511: Based on the directional deformation velocity data, extract the satellite incident angle and azimuth angle parameters of each permanent scatterer point, and perform a velocity component conversion operation on all velocity values ​​in the vertical direction to obtain the vertical deformation velocity matrix.

[0033] S512: Based on the vertical deformation velocity matrix, synchronously read the mapping parameter set from image coordinates to geographic coordinates recorded in the main image, perform affine transformation on the image coordinates of each permanent scatterer point, establish the mapping between spatial points and latitude and longitude, and obtain a list of geographic location identifiers;

[0034] S513: Based on the geographic location identifier list and the location index and velocity value in the vertical deformation velocity matrix, construct a structured record unit with geographic coordinates as the primary key and velocity value as the field, aggregate all record units and sort them by time series to establish the InSAR deformation time series analysis results of the water conservancy hub.

[0035] A time-series InSAR deformation analysis system for water conservancy projects includes:

[0036] The radar data construction module is used to perform S1: filter multiple C-band synthetic aperture radar images covering the sluice gate, pumping station and dam structure area, acquire SAR images, simultaneously collect SRTM digital elevation model data, POD orbit data and GACOS atmospheric delay data, organize the sequence according to the image time label, and construct the interferometric radar dataset.

[0037] The interferometric pair generation module is used to perform S2: based on the SAR images in the interferometric radar dataset, set the main image and the auxiliary image, perform registration operation and perform terrain phase calculation, and after superimposing the data of each interferometric pair, strip the terrain components and construct the differential interferogram structure matrix.

[0038] The scatterer modeling module is used to perform S3: calculate the amplitude deviation index of each pixel according to the differential interferogram structure matrix, filter permanent scatterer points that are less than the preset deviation threshold, record the spatial index and set it as the interferometric modeling control point, analyze the phase change characteristics of adjacent PS points, and establish the interferometric phase expression matrix.

[0039] The deformation parameter calculation module is used to execute S4: based on the interference phase expression matrix, calculate the deformation velocity and elevation error of each permanent scatterer point, establish a coordinate-velocity mapping, and generate directional deformation velocity data;

[0040] The time series analysis construction module is used to execute S5: based on the directional deformation velocity data, perform the LOS velocity to vertical velocity conversion, combine the coordinate mapping parameters in the main image, map to the geographic coordinate system, and construct the InSAR deformation time series analysis results of the water conservancy hub.

[0041] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0042] In this invention, a radar image dataset covering key structural areas is constructed, topographic and atmospheric interference factors are removed, long-term stable scatterer points are extracted, and their deformation rate and elevation error are calculated to achieve continuous temporal analysis of structural deformation. The accuracy and anti-interference capability of the results are improved by multi-source data fusion and phase expression optimization. High-quality deformation information inversion is completed by combining orbital parameters and atmospheric correction data, and the deformation results of the structural area are output in the form of geographic coordinate system. This enhances the ability to identify the stability change trend of key structures and achieves efficient monitoring of water conservancy facilities with wide coverage, high timeliness, and low dependence. Attached Figure Description

[0043] Figure 1 This is a flowchart of the main steps of the present invention;

[0044] Figure 2 This is a flowchart of the interferometric radar dataset acquisition process of this invention;

[0045] Figure 3 This is a flowchart of the process for obtaining the differential interferogram structure matrix of the present invention;

[0046] Figure 4 This is a flowchart of the process for obtaining the interference phase representation matrix of the present invention;

[0047] Figure 5 This is a flowchart of the process for acquiring directional deformation velocity data in this invention;

[0048] Figure 6 This is a flowchart illustrating the process of obtaining the InSAR deformation time series analysis results for a water conservancy hub according to the present invention. Detailed Implementation

[0049] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0050] In the description of this invention, it should be understood that the terms "length," "width," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," and "outer," etc., indicating orientation or positional relationships, are based on the orientation or positional relationships shown in the accompanying drawings and are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of the invention. Furthermore, in the description of this invention, "a plurality of" means two or more, unless otherwise explicitly specified.

[0051] Please see Figure 1 A time-series InSAR deformation analysis method for water conservancy projects includes the following steps:

[0052] S1: Select multiple C-band synthetic aperture radar images covering the sluice gate, pumping station, and dam structure areas, set the imaging mode to interferometric wide swath mode and the polarization mode to VV mode, acquire SAR master data with spatial consistency, simultaneously acquire SRTM digital elevation model data, POD orbit data and GACOS atmospheric delay data, and organize them according to the image time label sequence to construct an interferometric radar dataset.

[0053] S2: Based on SAR images in the interferometric radar dataset, analyze the Doppler frequency shift and spatiotemporal baseline length, set the images that meet the minimum baseline requirement as the main images, and the others as auxiliary images, perform registration operations and generate interferometric pair image sequences, combine with the SRTM digital elevation model to calculate the terrain phase, and after superimposing the data of each interferometric pair, strip the terrain components and construct the differential interferogram structure matrix.

[0054] S3: Based on the differential interferogram structure matrix, calculate the amplitude deviation index of each pixel and compare it with the preset deviation threshold. Select permanent scatterer points whose amplitude deviation index is less than the preset deviation threshold, record the corresponding spatial index and set it as the interferometric modeling control point. Analyze the phase change characteristics of adjacent PS points and establish the interferometric phase expression matrix.

[0055] S4: Based on the interferometric phase expression matrix, POD orbital data and GACOS atmospheric delay data are extracted, orbital error terms and atmospheric residual terms are removed, the remaining expression is subjected to least squares fitting operation, the deformation velocity and elevation error values ​​of each permanent scatterer point are calculated, coordinate-velocity mapping is established, and directional deformation velocity data are generated.

[0056] S5: Based on the directional deformation velocity data, perform the LOS-to-velocity conversion to vertical velocity, combine the coordinate mapping parameters in the main image, construct a list of structural labels mapped to the geographic coordinate system, and construct the InSAR deformation time series analysis results of the water conservancy hub.

[0057] The interferometric radar dataset includes spatially consistent SAR master data, image time-stamped serialized data, SRTM digital elevation model information, POD precision orbit data, atmospheric delay preprocessing data, differential interferogram structure matrix including terrain stripping interferogram sequence, phase difference calculation unit, and structure mapping relationship index, interferometric modeling control points including spatial location index, permanent scatterer feature identifier, and phase change measure, directional deformation velocity data including deformation velocity values ​​of permanent scatterer points, elevation error information, and LOS directional velocity mapping relationship, and InSAR deformation time series analysis results of water conservancy projects including vertical deformation velocity distribution, geographic coordinate system structure label list, and deformation monitoring target index.

[0058] Please see Figure 2 Step S1 is as follows:

[0059] S111: Acquire multiple C-band synthetic aperture radar images covering the sluice gate, pumping station, and dam structure area. Set the imaging mode to interferometric wide-swath mode and the polarization mode to VV mode. Arrange the images sequentially according to the imaging time label of each image and match and filter the spatial coverage range. Remove data frames with spatial intersection less than the preset coverage threshold to obtain the regional matching radar image set.

[0060] This embodiment selects a large-scale water conservancy project (including the main dam, ship lock, and a 5-kilometer buffer zone) located in the Yangtze River Basin as the monitoring target, with the geographical range defined as 30.82° to 30.85° north latitude and 111.00° to 111.05° east longitude. The execution process first accesses the European Space Agency (ESA) Copernicus Data Center and retrieves archived data from the Sentinel-1A satellite via the API interface. The specific search criteria are set as follows: imaging mode selected as Interferometric Wide Swath (IW), polarization locked as Vertical-Vertical (VV), orbit direction as Ascending, and time span set from January 1, 2023 to December 31, 2023, returning a total of 32 single-view complex (SLC) image data.

[0061] After obtaining the raw data list, strict spatial coverage filtering is required. First, the Footprint polygon coordinates (WKT format) in the metadata of each SLC image are parsed and projected onto the WGS84 coordinate system along with the preset monitoring target area polygon (ROI). The intersection area of ​​the two is then calculated using polygon Boolean operations. ) and the total area of ​​the monitoring target area ( Preset coverage threshold Only when At that time, the data frame was marked as valid. After screening, two images with incomplete coverage due to orbital offset were removed, retaining 30 fully covered images. Subsequently, the starting imaging time tag (in YYYYMMDDThhmmss format) was extracted from the filename of each image, and an index list was built according to the chronological order. ,in The image is from issue 20230105. The image from issue 20231225 will serve as the base input for subsequent processing.

[0062] S112: Based on the regional matching radar image set, synchronously collect the corresponding SRTM digital elevation model data, POD orbit data and GACOS atmospheric delay data according to the time label, construct the binding relationship between the image and the three types of auxiliary data, and generate an image auxiliary data structure set.

[0063] Based on certainty The system automatically and concurrently accesses three heterogeneous data sources to perform auxiliary data acquisition. First, based on the image center coordinates, it downloads the corresponding SRTMGL1 (1 arcsecond, approximately 30-meter resolution) digital elevation model tiles from the NASA Earthdata server. The data format is HGT, and the coverage area needs to be extended 0.5 degrees beyond the image boundary to prevent resampling edge effects. Second, it determines the imaging time for each Sentinel-1 image. Visit the official ESA POD Hub and search for the publication time in [year]. The precise orbital ephemerides (AUX_POEORB) obtained after [number] days ensure orbital position accuracy better than 5 centimeters; if no precise orbit is published, an error message will be displayed and the process will wait, and the use of predicted orbits is strictly prohibited. Finally, access the Newcastle University GACOS service platform, submit the image imaging time (accurate to the minute) and spatial extent, and download the corresponding tropospheric zenith delay (ZTD) raster data (.ztd format).

[0064] After the data download is complete, a hash structure set named ImgAuxStruct is constructed. This structure uses the image unique identifier (UUID) as the primary key and contains four fields: slc_path (local SLC image path), dem_path (mosaiced and cropped DEM path), orb_path (corresponding orbit file path), and atm_path (time-synchronized GACOS data path). During this process, the SRTM data needs to be preprocessed: HGT files are read, invalid values ​​are filled (NoData is repaired using bilinear interpolation), and the geoid elevation (EGM96) is converted to WGS84 ellipsoidal height using the geoid difference model interpolation formula to ensure that the elevation datum is consistent with the GPS and satellite orbit coordinate systems.

[0065] S113: Based on the image auxiliary data structure set, perform pairing operations on image frames of adjacent imaging times, and combine them with the corresponding elevation, orbit and atmospheric delay data structures to obtain the interferometric radar dataset;

[0066] Based on the ImgAuxStruct set, an interferometric network is constructed using the Short Baseline Set (SBAS) strategy. First, the imaging time of all images is read, and the time interval between adjacent images is calculated. Set the maximum time baseline threshold. The rule is set to allow pairing of only 1 to 3 adjacent images to suppress temporal decorrelation. Based on this rule, for time-ordered image sequences... Generate matching combinations and .For example, (20230105) will be respectively with (20230117) and (20230129) Pairing.

[0067] For each given image combination (main image) With auxiliary images The process involves calling the corresponding auxiliary data structures. SRTM DEM data is projected onto the main image radar geometric coordinate system (Range-Doppler Coordinates) to generate a reference terrain phase; POD orbit data is read to calculate the satellite's position and velocity vectors in the Earth-fixed coordinate system; GACOS data is read and resampled to radar image resolution. Finally, an interferometric radar dataset is established, which consists of... Composed of data packets ( (Total number of pairs) Each data packet encapsulates the primary and secondary image paths, the registered reference DEM array, the satellite state vectors at two time points, and the corresponding atmospheric delay difference map, providing a fully aligned input source for subsequent interferometric processing.

[0068] Please see Figure 3 Step S2 is as follows:

[0069] S211: Based on SAR images in the interferometric radar dataset, extract the Doppler frequency center value and spatiotemporal baseline length of each image, compare the baseline length of the main image with the preset baseline threshold, retain the image frames that meet the baseline conditions as the main image, and mark the other image frames as auxiliary images, and generate a main and auxiliary image classification label set.

[0070] Metadata for paired SAR images is read one by one from the interferometric radar dataset. First, for each image, the Doppler center frequency is calculated using its focusing parameters. The specific calculation employs the maximum energy method: the energy peak location is located in the azimuth spectrum. For Sentinel-1 TOPS mode data, the Doppler frequency variation in the beam overlap region needs to be additionally considered, and this is achieved by extracting the phase difference in the Burst overlap region. The refined estimate ensures that the Doppler center frequency estimation error is less than 10Hz.

[0071] Then, each pair of images (main image) is calculated. auxiliary images The spatiotemporal baseline. Spatial vertical baseline. The calculation is based on the orbital state vector:

[0072]

[0073] in, and These are the position vectors of the primary and secondary satellites at the time of imaging. For side-view average angle, This represents the baseline inclination angle. Vertical baseline. Derived directly from the imaging time difference. A strict baseline screening threshold is set: the absolute value of the vertical baseline. Meters, Doppler center frequency difference Hz. The initially generated pairing list is then filtered a second time to remove interference pairs that exceed a threshold. For example, if a pairing is calculated to... If the image is not valid, it will be marked as "invalid" and will not be included as a primary image in subsequent core processing; only images that meet the criteria will be retained. For image combinations, generate a set of primary and secondary image classification identifiers containing valid pairing indices and baseline parameters.

[0074] S212: Based on the classification and identification set of master and auxiliary images, perform image registration operation on each pair of master and auxiliary images, adjust the image position relationship based on the difference of pixel grid coordinates, perform interferometric superposition operation on the registered images, extract the interferometric phase image sequence, and index and encode it according to the master and auxiliary time labels to obtain the interferometric pair image sequence set;

[0075] Based on the filtered primary and secondary image classification identifier set, high-precision registration is performed on each pair of images. The registration process consists of two stages: First, coarse registration is performed, using orbital parameters and SRTM DEM to calculate the predicted position of the secondary image pixels in the primary image, and calculating pixel-level offsets using cross-correlation functions to correct large-scale translations; then, fine registration is performed using the Enhanced Spectral Diversity (ESD) method. ESD leverages the characteristic that the phase difference in the Burst overlap region of the TOPS mode is extremely sensitive to azimuth shift, minimizing the phase difference in the overlap region... :

[0076]

[0077] Solve for sub-pixel offset in azimuth direction ,in The Doppler frequency difference in the overlapping area is used. This step is performed iteratively until the azimuth registration accuracy reaches 0.001 pixels (approximately 1.4 cm on the ground).

[0078] After registration, the master and slave images are multiplied by conjugate to generate an interferogram. ,in The complex conjugate of the secondary image is used to extract the phase component of the interferogram. The range is obtained as The entangled phase diagram. Finally, based on the time tag of the main image. The time stamps of the auxiliary images are used as filenames (e.g., 20230105_20230117) to serialize and store the generated phase maps, forming an interferometric pair image sequence set.

[0079] S213: Based on the interferometric image sequence set and SRTM digital elevation model data, calculate the terrain phase component in each interferogram, peel off the terrain phase from the interferometric phase map pixel by pixel, and perform sequence aggregation on the peeling results to establish a differential interferogram structure matrix.

[0080] Read the interferometric phase map of each image in the interferometric image sequence set. Using SRTM DEM data already registered to the radar coordinate system, the terrain phase components are calculated. The calculation formula is:

[0081]

[0082] in, The radar wavelength (5.546 cm for Sentinel-1). Slope distance Angle of incidence This represents the DEM elevation value.

[0083] The calculated simulated terrain phase is subtracted pixel by pixel from the original interferometric phase: This operation removed the main topographic stripes, revealing surface deformation and atmospheric delay signals. For all After performing this operation on the amplitude differential interferogram, all differential phase matrices are stacked in time order to construct a dimension-1. The three-dimensional array, namely the differential interferogram structure matrix, provides a clean phase input for time series analysis, in which... This indicates the number of rows in the image, which is the number of pixels in the vertical direction. This indicates the number of columns in the image, which is the number of pixels in the horizontal direction. This represents the number of interferometric image pairs in a time series, i.e., the number of different time points or image pairs involved in the calculation in InSAR (Synthetic Aperture Radar Interferometry) analysis.

[0084] Please see Figure 4 Step S3 is as follows:

[0085] S311: Obtain the differential interferogram structure matrix, extract the interference amplitude value sequence of each pixel, and perform standard deviation and mean statistics on the sequence using the following formula:

[0086] ;

[0087] The amplitude deviation index of each pixel is calculated to obtain the amplitude deviation index data; among which, This represents the amplitude deviation index of the i-th pixel. This represents the amplitude value of the i-th pixel in the j-th interferogram. Let represent the mean amplitude of the i-th pixel across all interferograms, and n represent the total number of interferograms;

[0088] Obtain the differential interferogram structure matrix, and then process each pixel in the matrix. Extract its in Aspect interferogram (here) The amplitude value sequence corresponding to the number of images involved in the calculation (actually the length of the time series). The amplitude value here It is the square root of the absolute radiance value of the SLC image after radiometric calibration.

[0089] Calculate the amplitude deviation index using the formula The specific calculation steps are as follows:

[0090] Calculate the mean Calculate the arithmetic mean of the sequence. ;

[0091] Calculate the standard deviation : ;

[0092] Calculate the deviation index : ;

[0093] Select a strong reflection point on the surface of the sluice gate structure Extract its values ​​across 5 time phases (simplified example). The amplitude sequence of the radar echo. The data is derived from actual radar echo intensity, and the unit is dimensionless digital quantization (DN).

[0094] Table 1: Monitoring Points Amplitude sequence data

[0095]

[0096] As shown in Table 1, substitute the values ​​into the formula to calculate:

[0097] mean .

[0098] Variance calculation:

[0099] ;

[0100] ;

[0101] ;

[0102] ;

[0103] ;

[0104] The summation yields 129.5;

[0105] Divide ;

[0106] Standard deviation ;

[0107] Amplitude Deviation Index .

[0108] The result of 0.0042 is much smaller than the general threshold, indicating that this point has extremely high radiation stability and is suitable as a permanent scatterer. Performing this operation on all pixels in the entire image generates amplitude deviation index data at the same resolution.

[0109] S312: Based on the amplitude deviation index data, compare the amplitude deviation index of all pixels with the preset deviation threshold pixel by pixel, filter out the pixels that are less than the preset deviation threshold, extract the two-dimensional spatial index coordinates of the corresponding pixels in the interferogram, establish a coordinate structure list, and obtain the permanent scatterer spatial index set.

[0110] Set the decision threshold for permanent scatterer candidate points (PSCs). The threshold is selected based on statistical principles; when the signal-to-noise ratio (SNR) of the scatterer is high, (Phase standard deviation). This means that the phase noise standard deviation is less than 0.25 radians, which can guarantee the success rate of subsequent phase unwrapping.

[0111] Traverse the amplitude deviation index data, if a certain pixel's If the condition is met, then mark it as PSC. Extract the coordinates of all pixels that meet the criteria. Stored in a dynamic list. In water conservancy hub areas, concrete dams, metal gates, and rock revetments typically exhibit low... The value is (0.1-0.2), while the water surface area, due to specular reflection and ripples, has a different value. The value is typically greater than 0.6. This step effectively filters out noise points in water surfaces and vegetation-covered areas. The final permanent scatterer spatial index set contains approximately 15,000 discrete point coordinates.

[0112] S313: Based on the spatial index set of permanent scatterers, extract the phase value sequence of all adjacent permanent scatterer points in each interferogram, aggregate the phase difference sequence of adjacent point pairs, and construct a two-dimensional matrix structure with rows representing point pairs and columns representing time series, and output the interferometric phase expression matrix.

[0113] Based on the spatial index set, construct a Delaunay triangulation or a distance-constrained connectivity network, connecting adjacent PSC points into edges (Arc). For each edge (connecting point)... and points ), calculate its in the first Phase difference in amplitude interferogram:

[0114]

[0115] This operation can effectively offset the spatial correlation error shared between two points (such as the long-wavelength component of atmospheric delay).

[0116] Constructing a two-dimensional matrix Its number of rows represents the number of edges generated (e.g., 30,000 edges), and its number of columns represents the number of interferograms (e.g., 29). Matrix elements That is, the first The edge at the 1st The phase difference values ​​on the image. This matrix is ​​the output interference phase representation matrix.

[0117] Please see Figure 5 Step S4 is as follows:

[0118] S411: Based on the interferometric phase expression matrix, POD orbital data and GACOS atmospheric delay data are extracted, and synchronization is performed according to time label and spatial index. Vector difference operation is performed on the interferometric phase value and orbital error term. After stripping the orbital term, the atmospheric delay error field is resampled and subtracted from the residual phase matrix to obtain the error stripped phase matrix.

[0119] External auxiliary data is introduced to refine the interferometric phase representation matrix. First, orbital errors are addressed: although precise orbits are used, baseline-related phase slopes still remain. For each edge, an orbital error phase model is constructed based on the spatial differences between its two endpoints and the satellite baseline parameters.

[0120] ;

[0121] in Let the coordinates be range and azimuth. The linear trend of the residual phase is fitted using the least squares method and then... Subtract from the middle, It is a constant term. These are distance and azimuth coordinates. The coefficient.

[0122] Next, atmospheric delay is addressed. The previously acquired GACOS ZTD data is read, and the tropospheric zenith delay difference between the two endpoints of each edge at each imaging time is calculated. The formula is used to convert it into radar line-of-sight (LOS) phase delay:

[0123] ;

[0124] in Let the angle be the angle of incidence. The calculated... The matrix is ​​directly subtracted from the phase matrix after orbital correction. The residual phase after this step mainly consists of surface deformation and topographic residuals, yielding the error-stripped phase matrix. Experiments show that after introducing GACOS correction, the phase standard deviation of the interferogram is reduced from 1.2 rad to 0.4 rad, significantly improving the detection capability of small deformations.

[0125] S412: Based on the error stripping phase matrix, perform least squares linear fitting on the stripping phase of each permanent scatterer point in the time series. In the fitting coefficients, the intercept term represents the elevation error and the slope term represents the rate of change in the time dimension. The execution speed corresponds to the elevation, and obtain the deformation and elevation error estimation table.

[0126] For each edge in the phase matrix, the error is removed, and a time-series function model is established. Phase difference Expressed as deformation rate and elevation error

[0127] ;

[0128] in, For the first The time baseline of the interferogram For the first Spatial vertical baseline of the interferogram Where u is the radar wavelength and u is the slant range. Angle of incidence Let be the deformation rate to be solved. The elevation error to be solved is denoted as .

[0129] Constructing a system of linear equations ,in To observe the phase vector, , For the residual vector, These are the parameters to be determined.

[0130] Solve using the least squares method: ;

[0131] Select an edge connecting the dam crest and the dam foundation, and set the data for 3 interference pairs ( (Simplified calculation)

[0132] parameter: m, m, ( );

[0133] coefficient rad / m;

[0134] coefficient rad / m / year (assuming T is in years);

[0135] Table 2: Least Squares Inversion Example Data

[0136]

[0137] Construct a system of equations:

[0138] ;

[0139] ;

[0140] ;

[0141] Simplified coefficient matrix :

[0142] Row 1: ;

[0143] Row 2: ;

[0144] Row 3: ;

[0145] Solving this overdetermined system of equations (using numpy.linalg.lstsq in the actual code), we get:

[0146] m / year (45mm / yr). m.

[0147] The results indicate a relative deformation rate difference between the two endpoints of this edge, and that the SRTM DEM exhibits an elevation error of approximately 10 meters at this location. After solving for all edges, the deformation rate at each absolute point is calculated using weighted least squares adjustment or integration at a specific reference point (assuming the bedrock point far from the dam has a velocity of 0). and elevation error Generate a table of deformation and elevation error estimates.

[0148] S413: Based on the deformation and elevation error estimation table, extract the two-dimensional coordinate values ​​of each permanent scatterer point and the corresponding deformation rate value to form a coordinate-velocity mapping data frame, and uniformly project the velocity components of each point in the structure to the satellite line of sight to generate directional deformation velocity data.

[0149] Based on the deformation and elevation error estimation table, extract each permanent scatterer point. Two-dimensional image coordinates And the calculated line-of-sight deformation rate Elevation error Add it back to the original SRTM elevation to obtain the refined elevation. . Build contains The data frame of the quadruplet is output as directional deformation velocity data. At this time, the velocity is still a one-dimensional scalar along the direction of the satellite's line of sight.

[0150] Please see Figure 6 The S5 steps are as follows:

[0151] S511: Based on the directional deformation velocity data, extract the satellite incident angle and azimuth parameters of each permanent scatterer point, perform velocity component conversion calculation on all velocity values ​​in the vertical direction, and obtain the vertical deformation velocity matrix.

[0152] Read the incident angle file from the auxiliary data and extract the incident angle corresponding to each PS point position. Since the deformation of hydraulic facilities mainly manifests as settlement or uplift, assuming the horizontal displacement component is negligible, the line of sight is projected onto the velocity in the vertical direction using trigonometric geometric relationships. The projection formula is:

[0153] ;

[0154] For example, for a certain monitoring point, the measured mm / yr (away from the satellite), angle of incidence at this point .

[0155] calculate ;

[0156] mm / yr;

[0157] This value indicates that the point is subsiding at a rate of approximately 1.9 centimeters per year. Performing this transformation on all points constructs a vertical deformation rate matrix, which accurately reflects the physical deformation pattern of the Earth's surface.

[0158] S512: Based on the vertical deformation velocity matrix, synchronously read the mapping parameter set from the image coordinates to the geographic coordinates recorded in the main image, perform affine transformation on the image coordinates of each permanent scatterer point, establish the mapping between spatial points and latitude and longitude, and obtain a list of geographic location identifiers;

[0159] To map the results in the image coordinate system to the real world, a rigorous geocoding lookup table needs to be established. The geotransformation parameters of the main image are read, including the latitude and longitude of the top-left corner, pixel resolution, and rotation coefficient. The DopplerRange equation is then used, combined with the refined elevation... Perform reverse geocoding.

[0160] For each pixel Solve the following system of equations to find the geodetic coordinates. :

[0161] Distance equation: ;

[0162] Doppler equations: ;

[0163] Earth model equations: ;

[0164] Solve the above system of equations using Newton's iterative method, and obtain the image coordinates. Accurate conversion to WGS84 latitude and longitude coordinates Generate a list of corresponding geographic location identifiers;

[0165] In the above system of equations: Indicates the distance from the radar sensor to the target on the ground; The satellite's three-dimensional coordinates in the WGS84 coordinate system at the time of imaging; The three-dimensional coordinates of the target point on the ground in the WGS84 coordinate system are to be determined. The center frequency of the Doppler signal; The radar wavelength; The velocity vector of the satellite, The velocity vector of the target point on the ground surface (usually assumed to be 0 or to have known motion); The satellite's position vector. The position vector of the target point; With reference to the semi-major axis of the ellipsoid, For the reference ellipsoid's minor semi-axis, The geodetic height of the target point. Let be the radius of curvature of the y-axis.

[0166] S513: Based on the geographical location identifier list and the location index and velocity value in the vertical deformation velocity matrix, construct a structured record unit with geographical coordinates as the primary key and velocity value as the field, aggregate all record units and sort them by time series to establish the InSAR deformation time series analysis results of the water conservancy hub.

[0167] Integrate the vertical velocity data from S511 with the latitude and longitude data from S512. Construct the final GIS-compatible database structure. Each record contains: point ID, longitude, latitude, elevation (refined), annual average settlement rate (mm / yr), and deformation time series (cumulative deformation arranged by date).

[0168] Table 3: Examples of InSAR monitoring results for water conservancy projects

[0169]

[0170] As shown in Table 3, point PS_002 exhibits a significant settlement trend (-12.5 mm / yr), requiring a focused early warning. All records are clustered according to spatial location and sorted by time series to generate standardized Shapefiles or GeoJSON files, which constitute the final InSAR deformation time series analysis results for the water conservancy project. These results directly support the water conservancy department in conducting dam safety assessments, automatically marking areas with settlement rates exceeding -10 mm / yr as red warning zones.

[0171] A time-series InSAR deformation analysis system for water conservancy projects includes:

[0172] The radar data construction module is used to perform S1: filter multiple C-band synthetic aperture radar images covering the sluice gate, pumping station and dam structure area, acquire SAR images, simultaneously collect SRTM digital elevation model data, POD orbit data and GACOS atmospheric delay data, organize the sequence according to the image time label, and construct the interferometric radar dataset.

[0173] The interferometric pair generation module is used to perform S2: based on the SAR images in the interferometric radar dataset, it sets the main image and the auxiliary image, performs registration operation and performs terrain phase calculation, and after superimposing the data of each interferometric pair, it strips the terrain components and constructs the differential interferogram structure matrix.

[0174] The scatterer modeling module is used to execute S3: based on the differential interferogram structure matrix, calculate the amplitude deviation index of each pixel, filter permanent scatterer points that are less than the preset deviation threshold, record the spatial index and set it as the interferometric modeling control point, analyze the phase change characteristics of adjacent PS points, and establish the interferometric phase expression matrix.

[0175] The deformation parameter calculation module is used to execute S4: based on the interference phase expression matrix, calculate the deformation velocity and elevation error of each permanent scatterer point, establish a coordinate-velocity mapping, and generate directional deformation velocity data;

[0176] The time series analysis module is used to execute S5: based on the directional deformation velocity data, it performs the LOS-to-velocity conversion to vertical velocity, combines the coordinate mapping parameters in the main image, maps it to the geographic coordinate system, and constructs the InSAR deformation time series analysis results of the water conservancy hub.

[0177] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments that can be applied to other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for InSAR deformation time series analysis of a water conservancy hub, characterized in that, Includes the following steps: S1: Filter radar images of sluice gates, pumping stations and dam areas, acquire SAR images, and simultaneously collect SRTM digital elevation model data, POD orbit data and GACOS atmospheric delay data. Organize the images into a sequence according to the image time label and construct an interferometric radar dataset. S2: Based on the SAR images in the interferometric radar dataset, set the main image and auxiliary image, perform registration operation and terrain phase calculation, and after superimposing the data of each interferometric pair, strip the terrain components and construct the differential interferogram structure matrix. S3: Based on the differential interferogram structure matrix, calculate the amplitude deviation index of each pixel, filter permanent scatterer points that are less than the preset deviation threshold, record the spatial index and set it as the interferometric modeling control point, analyze the phase change characteristics of adjacent PS points, and establish the interferometric phase expression matrix. S4: Based on the interference phase expression matrix, calculate the deformation velocity and elevation error of each permanent scatterer point, establish a coordinate-velocity mapping, and generate directional deformation velocity data; S5: Based on the directional deformation velocity data, perform the conversion from LOS directional velocity to vertical directional velocity, combine the coordinate mapping parameters in the main image, map to the geographic coordinate system, and construct the InSAR deformation time series analysis results of the water conservancy hub; The specific steps for obtaining the interference phase representation matrix are as follows: S311: Obtain the differential interferogram structure matrix, extract the interference amplitude value sequence of each pixel, perform standard deviation and mean statistics on the sequence, and calculate the amplitude deviation index of each pixel to obtain the amplitude deviation index data. S312: Based on the amplitude deviation index data, compare the amplitude deviation index of all pixels with the preset deviation threshold pixel by pixel, filter out pixels that are less than the preset deviation threshold, extract the two-dimensional spatial index coordinates of the corresponding pixels in the interferogram, establish a coordinate structure list, and obtain the permanent scatterer spatial index set. S313: Based on the permanent scatterer spatial index set, extract the phase value sequence of all adjacent permanent scatterer points in each interferogram, aggregate the phase difference sequence of adjacent point pairs, and construct a two-dimensional matrix structure with rows representing point pairs and columns representing time series, and output the interferometric phase expression matrix. The formula for calculating the amplitude deviation index is as follows: ； Indicates the first The amplitude deviation index of a pixel. Indicates the first Pixel in the Amplitude values ​​in the amplitude interferogram Indicates the first The mean amplitude of a pixel across all interferograms. This indicates the total number of interferograms.

2. The InSAR deformation time series analysis method for water conservancy projects according to claim 1, characterized in that, The interferometric radar dataset includes spatially consistent SAR master data, image time-stamped serialized data, SRTM digital elevation model information, POD precision orbit data, and atmospheric delay preprocessing data. The differential interferogram structure matrix includes terrain stripping interferogram sequences, phase difference calculation units, and structure mapping relationship indexes. The interferometric modeling control points include spatial location indexes, permanent scatterer feature identifiers, and phase change measures. The directional deformation velocity data includes deformation velocity values ​​of permanent scatterer points, elevation error information, and LOS directional velocity mapping relationships. The InSAR deformation time series analysis results of the water conservancy hub include vertical deformation velocity distribution, geographic coordinate system structure label list, and deformation monitoring target index.

3. The InSAR deformation time series analysis method for water conservancy projects according to claim 1, characterized in that, The specific steps for acquiring the interferometric radar dataset are as follows: S111: Acquire multiple C-band synthetic aperture radar images covering the sluice gate, pumping station, and dam structure area. Set the imaging mode to interferometric wide-swath mode and the polarization mode to VV mode. Arrange the images sequentially according to the imaging time label of each image and match and filter the spatial coverage range. Remove data frames with spatial intersection less than the preset coverage threshold to obtain the regional matching radar image set. S112: Based on the region-matching radar image set, synchronously collect the corresponding SRTM digital elevation model data, POD orbit data and GACOS atmospheric delay data according to the time label, construct the binding relationship between the images and the three types of auxiliary data, and generate an image auxiliary data structure set. S113: Based on the image auxiliary data structure set, perform pairing operations on image frames of adjacent imaging times, and combine them with the corresponding elevation, orbit and atmospheric delay data structures to obtain the interferometric radar dataset.

4. The InSAR deformation time series analysis method for water conservancy projects according to claim 1, characterized in that, The specific steps for obtaining the differential interferogram structure matrix are as follows: S211: Based on the SAR images in the interferometric radar dataset, extract the Doppler frequency center value and spatiotemporal baseline length of each image, compare the baseline length of the main image with the preset baseline threshold, retain the image frames that meet the baseline conditions as the main image, and mark the other image frames as auxiliary images, and generate a main and auxiliary image classification identifier set. S212: Based on the main and auxiliary image classification identifier set, perform image registration operation on each pair of main and auxiliary images, adjust the image position relationship based on the pixel grid coordinate difference, perform interferometric superposition operation on the registered images, extract the interferometric phase image sequence, and index and encode it according to the main and auxiliary time labels to obtain the interferometric pair image sequence set; S213: Based on the interferometric image sequence set and SRTM digital elevation model data, calculate the terrain phase component in each interferogram, peel off the terrain phase from the interferometric phase map pixel by pixel, and perform sequence aggregation on the peeling results to establish a differential interferogram structure matrix.

5. The InSAR deformation time series analysis method for water conservancy projects according to claim 1, characterized in that, The specific steps for acquiring the directional deformation velocity data are as follows: S411: Based on the interferometric phase expression matrix, extract the POD orbital data and GACOS atmospheric delay data, perform synchronization operation according to time label and spatial index, perform vector difference operation on the interferometric phase value and orbital error term, strip the orbital term, perform atmospheric delay error field grid resampling and subtract it from the residual phase matrix to obtain the error stripped phase matrix. S412: Based on the error stripping phase matrix, perform least squares linear fitting on the stripping phase of each permanent scatterer point in the time series. In the fitting coefficients, the intercept term represents the elevation error and the slope term represents the rate of change in the time dimension. The execution speed corresponds to the elevation, and obtain the deformation and elevation error estimation table. S413: Based on the deformation and elevation error estimation table, extract the two-dimensional coordinate values ​​of each permanent scatterer point and the corresponding deformation rate value to form a coordinate-velocity mapping data frame, and uniformly project the velocity components of each point in the structure to the satellite line of sight to generate directional deformation velocity data.

6. The InSAR deformation time series analysis method for water conservancy projects according to claim 1, characterized in that, The specific steps for obtaining the InSAR deformation time series analysis results of the water conservancy hub are as follows: S511: Based on the directional deformation velocity data, extract the satellite incident angle and azimuth angle parameters of each permanent scatterer point, and perform a velocity component conversion operation on all velocity values ​​in the vertical direction to obtain the vertical deformation velocity matrix. S512: Based on the vertical deformation velocity matrix, synchronously read the mapping parameter set from image coordinates to geographic coordinates recorded in the main image, perform affine transformation on the image coordinates of each permanent scatterer point, establish the mapping between spatial point location and latitude and longitude, and obtain a list of geographic location identifiers; S513: Based on the geographic location identifier list and the location index and velocity value in the vertical deformation velocity matrix, construct a structured record unit with geographic coordinates as the primary key and velocity value as the field, aggregate all record units and sort them by time series to establish the InSAR deformation time series analysis results of the water conservancy hub.

7. A time-series InSAR deformation analysis system for hydraulic engineering projects, characterized in that, The system is used to implement the InSAR deformation time series analysis method for water conservancy hubs as described in any one of claims 1-6, including: The radar data construction module is used to perform S1: filter multiple C-band synthetic aperture radar images covering the sluice gate, pumping station and dam structure area, acquire SAR images, simultaneously collect SRTM digital elevation model data, POD orbit data and GACOS atmospheric delay data, organize the sequence according to the image time label, and construct the interferometric radar dataset. The interferometric pair generation module is used to perform S2: based on the SAR images in the interferometric radar dataset, set the main image and the auxiliary image, perform registration operation and perform terrain phase calculation, and after superimposing the data of each interferometric pair, strip the terrain components and construct the differential interferogram structure matrix. The scatterer modeling module is used to perform S3: calculate the amplitude deviation index of each pixel according to the differential interferogram structure matrix, filter permanent scatterer points that are less than the preset deviation threshold, record the spatial index and set it as the interferometric modeling control point, analyze the phase change characteristics of adjacent PS points, and establish the interferometric phase expression matrix. The deformation parameter calculation module is used to execute S4: based on the interference phase expression matrix, calculate the deformation velocity and elevation error of each permanent scatterer point, establish a coordinate-velocity mapping, and generate directional deformation velocity data; The time series analysis construction module is used to execute S5: based on the directional deformation velocity data, perform the LOS velocity to vertical velocity conversion, combine the coordinate mapping parameters in the main image, map to the geographic coordinate system, and construct the InSAR deformation time series analysis results of the water conservancy hub.

Images