Interference image pair screening method and application thereof

Abstract

The invention relates to the technical field of image processing, in particular to an interference image pair screening method and application thereof. The method comprises the following steps of: constructing an interference baseline network undirected graph by taking an acquisition moment of an SAR (Synthetic Aperture Radar) image in an interference image pair as a node and a connecting line between the interference image pairs as an edge; all target edges, smaller than a preset threshold value, on each graph node in the interference baseline network undirected graph are removed in sequence according to the sequence of the collection moments, and the preset threshold value is the mean value of coherence coefficients of all the edges connected with the graph nodes; and taking the residual interference image pair in the interference baseline network undirected graph after the target edge is removed as a screened target interference image pair. The objective of the invention is to solve the problem of how to obtain a high-precision InSAR.

Description

Technical Field

[0001] This application relates to the field of image processing technology, and particularly relates to a method for screening interferometric image pairs and its application. Background Art

[0002] InSAR (Interferometric Synthetic Aperture Radar) technology is an active microwave remote sensing system. Its basic working principle is to use a radar sensor to emit electromagnetic waves to the ground, and after receiving the scattered echo signals from ground objects, it forms an image. According to the phase information between the SAR (Synthetic Aperture Radar) sensor and the ground target recorded in the SAR image, the SAR images of the same area on repeated orbits are subjected to interferometric superposition processing, and the phase differences are compared to infer surface elevation or surface deformation information. InSAR can stack SAR images at multiple times to form multiple interferograms, obtain time-series deformation information, and capture the temporal evolution of the surface. It has extensive applications in fields such as landslide identification, glacier displacement monitoring, subsidence monitoring in mining areas, and urban deformation monitoring.

[0003] Changes in environmental factors such as surface vegetation growth, human activities, water vapor changes, and geomorphic changes will cause decorrelation noise between SAR images, affecting the quality of interferometric image pairs. The quality of interferometric image pairs is crucial for the high-precision inversion of surface deformation. Selecting interferometric image pairs with poor quality or a small number will reduce the accuracy of InSAR deformation inversion. Selecting a large number of interferometric image pairs will increase the computational cost and may not necessarily obtain high-precision surface deformation.

[0004] In view of this, this application aims to propose a selection method for screening high-quality interferometric image pairs from a large amount of interferometric image pair data, so as to obtain high-precision InSAR results. Summary of the Invention

[0005] The main purpose of this application is to provide a method for screening interferometric image pairs, aiming to solve the problem of how to obtain high-precision InSAR.

[0006] To achieve the above object, a method for screening interferometric image pairs provided by this application includes:

[0007] S10, taking the acquisition time of the SAR images in the interferometric image pair as nodes and the connection lines between the interferometric image pairs as edges, constructing an undirected graph of the interferometric baseline network;

[0008] S20, in the order of the sequence of the acquisition times, successively remove all target edges with values less than a preset threshold on each graph node in the undirected graph of the interference baseline network, where the preset threshold is the mean of the coherence coefficients of all the edges connected to the graph node;

[0009] S30, use the remaining interferogram pairs in the undirected graph of the interference baseline network after removing the target edges as the filtered target interferogram pairs.

[0010] Optionally, before the S10, the following steps are further included:

[0011] S40, calculate the atmospheric delay phase value in the initial interferogram pair;

[0012] S50, subtract the atmospheric delay phase value from the initial interferogram to obtain the interferogram pair for constructing the undirected graph of the interference baseline network.

[0013] Optionally, in the S40, the calculation expression of the atmospheric delay phase value is:

[0014]

[0015] In the formula, ZTDk is the atmospheric delay phase value, representing the total zenith delay of the vertical stratification component and the horizontal turbulence component set, k is the coordinate position; T represents the turbulence signal, x k is the site coordinate vector in the local geocentric coordinate system; L0 represents the stratification component delay at sea level; represents the remaining unmodeled residuals, including unmodeled stratification and turbulence signals; is the altitude scale;

[0016] Among them, the altitude scale The calculation expression is:

[0017]

[0018] In the formula, represents the altitude scale, represents the lowest altitude, represents the highest altitude.

[0019] Optionally, the calculation expression of the coherence coefficient is:

[0020]

[0021] In the formula, represents the amplitude and phase of the master image, represents the amplitude and phase of the slave image, represents the complex conjugate of, used to calculate the coherence.

[0022] In addition, to achieve the above object, the present application also provides a method for surface deformation inversion, and the method for surface deformation inversion includes the following steps:

[0023] Obtain a target interferometric pair obtained by screening based on the screening method of the interferometric pair described above;

[0024] Select a target interferogram generated by superimposing SAR images corresponding to two different acquisition times in the target interferometric pair, and calculate the target interferometric pair corresponding to the target interferogram. The target interferometric pair is composed of the cumulative deformation amount in the radar line-of-sight direction, the remaining topographic phase in the differential interferogram, the atmospheric delay phase, and the sum of the decorrelation noise;

[0025] Calculate the velocity vector according to the target interferometric pair and the acquisition time;

[0026] Calculate the minimum norm solution of the velocity vector by using the singular value decomposition method, and integrate the minimum norm solution of the velocity vector to obtain the deformation amount corresponding to the target interferogram;

[0027] According to the numerical interval where the deformation amount is located, determine the layover area, shadow area, foreshortening area, and non-deformed area in the SAR image, and mask the layover area and the shadow area to complete the surface deformation inversion.

[0028] Optionally, the calculation expression of the velocity vector is:

[0029]

[0030] In the formula, represents the i-th target interferogram, represents the acquisition time corresponding to the i-th target interferogram, represents the (i - 1)-th target interferogram, the acquisition time corresponding to the (i - 1)-th target interferogram.

[0031] Optionally, the step of calculating the minimum norm solution of the velocity vector by using the singular value decomposition method and integrating the minimum norm solution of the velocity vector to obtain the deformation amount corresponding to the target interferogram includes:

[0032] Let the interferometric pair value of the i-th interferogram be expressed as:

[0033]

[0034] In the formula, T A represents the acquisition time A, T B represents the acquisition time B, v pIndicates the surface deformation rate.

[0035] Integrals over each time period within the time interval between the master and slave images are represented by an M×N matrix:

[0036]

[0037] The generalized inverse matrix of matrix B is calculated using the singular value decomposition method, and the minimum norm solution of the velocity vector is calculated based on the generalized inverse matrix;

[0038] The minimum norm solution of the velocity vector is integrated to obtain the deformation amount corresponding to the target interferogram.

[0039] In addition, to achieve the above object, the present application also provides an application of the screening method for the interferogram pair as described above in surface deformation inversion.

[0040] In addition, to achieve the above object, the present application also provides a computer system, the computer system includes: a memory, a processor, and a computer program stored on the memory and executable on the processor, and when the computer program is executed by the processor, it implements the steps of the screening method for the interferogram pair as described in any one of the above or the surface deformation inversion method as described in any one of the above.

[0041] In addition, to achieve the above object, the present application also provides a computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by the processor, it implements the steps of the screening method for the interferogram pair as described in any one of the above or the surface deformation inversion method as described in any one of the above.

[0042] The present application has at least the following beneficial effects:

[0043] 1. A method based on graph theory is used to construct an undirected graph of the interferometric baseline network, and all target edges less than a preset threshold on each graph node in the undirected graph of the interferometric baseline network are removed, where the preset threshold is the mean of the coherence coefficients of all edges connected to the graph node, which not only removes low-quality interferogram pairs but also ensures the integrity of the interferometric baseline network;

[0044] 2. The atmospheric delay phase in the interferogram pair data is corrected to improve the accuracy of the InSAR data structure;

[0045] 3. For the screened interferogram pairs, surface deformation information is inverted based on the SBAS-InSAR technology, and geometric distortions in the SAR images are considered during the inversion process. The areas with geometric distortions are visually marked and then masked to ensure the accuracy of the InSAR results. Description of the Drawings

[0046] Figure 1 Schematic flowchart of the first embodiment of the interference image pair screening method of the present application;

[0047] Figure 2 Schematic diagram of the undirected graph involved in the embodiment of the present application;

[0048] Figure 3 Schematic diagram of the directed simple graph involved in the embodiment of the present application;

[0049] Figure 4 Schematic diagram of the undirected graph of the interference baseline network before removing the target edge involved in the embodiment of the present application;

[0050] Figure 5 Schematic diagram of the undirected graph of the interference baseline network after removing the target edge involved in the embodiment of the present application;

[0051] Figure 6 Schematic diagram of the phase standard deviation before and after GACOS correction involved in the embodiment of the present application;

[0052] Figure 7 Radar visibility geometric distortion recognition map involved in the embodiment of the present application;

[0053] Figure 8 Interference image pair connection graph based on the fully connected method involved in the embodiment of the present application;

[0054] Figure 9 Interference image pair connection graph based on the average coherence coefficient threshold method involved in the embodiment of the present application;

[0055] Figure 10 Interference image pair connection graph based on the small baseline set method involved in the embodiment of the present application;

[0056] Figure 11 Interference image pair connection graph based on the interference image pair screening method proposed in the present application involved in the embodiment of the present application.

[0057] Figure 12 Schematic diagram of the architecture of the hardware operating environment of the computer system involved in the embodiment of the present application.

[0058] The realization, functional characteristics and advantages of the purpose of the present application will be further described in conjunction with the embodiments and with reference to the accompanying drawings. Specific embodiments

[0059] To better understand the above technical solutions, the exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided so that the present disclosure can be more thoroughly understood and the scope of the present disclosure can be fully communicated to those skilled in the art.

[0060] First Embodiment

[0061] In order to reduce the impact of decoherence noise on InSAR deformation inversion, a large number of studies have been carried out on the optimal selection of InSAR interferogram pairs. At present, there are mainly three common methods for the optimal selection of InSAR interferogram pairs: The first method is based on the past experience of InSAR data processors, and manually eliminates the interferogram pairs with poor quality by means of manual visual interpretation, only retaining the relatively ideal interferogram pairs. This method is time-consuming and laborious, with low efficiency, and it is difficult to meet the needs of processing long-time series InSAR data with a large amount of data. Moreover, this method relies on the prior knowledge of data processors and has a certain degree of subjectivity. The second method is the method of screening interferogram pairs using the average coherence coefficient of all interferogram pairs as a threshold, that is, the average coherence coefficient threshold method. Existing studies have shown that this method can achieve relatively ideal results in urban areas with little change in coherence. However, in mountainous areas with complex terrain, coherence often changes with seasons. In summer, there is abundant precipitation and high vegetation coverage, resulting in low coherence; in winter, it is relatively dry and the vegetation coverage is low, resulting in high coherence. Due to the large difference in coherence, the subsequent disconnection of the interferometric baseline network may occur after overall averaging using this method, making it impossible to perform subsequent surface deformation information inversion. The third method is to perform threshold constraints on the temporal and spatial baselines, that is, the short spatio-temporal baseline method. The quality of this method depends on the setting of the temporal and spatial baseline thresholds. When the threshold is set too large, it will increase temporal and spatial decoherence and generate redundant low-coherence interferogram pairs; when the threshold is set too small, it will lead to a reduction in the number of connected interferogram pairs and reduce the inversion accuracy of surface deformation information. Therefore, this method has great instability, and different regions may have different optimal thresholds. The setting of the temporal and spatial baseline thresholds still depends more on past experience.

[0062] In this embodiment, aiming at the deficiencies of the above InSAR interferogram pair optimization selection method, a screening method for interferogram pairs based on graph theory is proposed in this embodiment:

[0063] Referring to Figure 1 , in this embodiment, the screening method for the interferogram pairs includes the following steps:

[0064] S10. Taking the acquisition time of SAR images in the interferometric image pair as nodes and the connection lines between interferometric image pairs as edges, construct an undirected graph of the interferometric baseline network;

[0065] In this embodiment, an interferometric image pair consists of at least two SAR images of the same area at different times and / or different viewpoints. Each SAR image records a time stamp during acquisition. Taking this time stamp as the acquisition time of the SAR image, using it as a graph node, and taking the connection lines between interferometric image pairs as edges, construct an undirected graph of the interferometric baseline network.

[0066] It should be noted that in the undirected graph of the interferometric baseline network, a target node that has a connection line with another node, that is, an interferometric image pair with interference, and the connection line between the target node and other nodes can be one or more. For example, if the nodes connected to node A are node B and node C, then node A and node B form an interferometric image pair, and node A and node C also form an interferometric image pair.

[0067] It should be noted that in this embodiment, graph theory is used to screen interferometric image pairs. Graph theory is a branch of mathematics that studies the objective world in the form of graphs. The graph in graph theory consists of some points and the connection lines connecting these points. The points represent things or objects, and the connection lines between points are called edges, representing the relationships between things or objects (used to represent interferometric image pairs in this embodiment). Each edge has a corresponding value, called the weight of the edge, representing the importance between two points (i.e., the importance of this interferometric image pair).

[0068] In graph theory, the graph G=(V,E) refers to two sets (V,E). The set V is the point set of the graph, and the set E is the edge set of the graph. The graph is a data structure composed of the relationship between the point set V and the edge set E. Refer to Figure 2 the schematic diagram of the undirected graph shown. Suppose there is a connection line from x to y in graph G, denoted as x→y, where point x is called the starting point of the edge and point y is called the ending point of the edge. If both x→y and y→x exist at the same time, then these two edges are combined into one. At this time, graph G is called an undirected graph.

[0069] In addition, refer to Figure 3 the schematic diagram of the directed simple graph shown. Suppose in graph G, there is only one edge from any starting point x to y, and there is no edge from itself to itself, then graph G is also called a simple graph or a directed simple graph.

[0070] S20. In the order of the sequence of the acquisition times, sequentially remove all target edges with values less than a preset threshold on each graph node in the undirected graph of the interferometric baseline network, where the preset threshold is the average value of the coherence coefficients of all edges connected to the graph node;

[0071] In this embodiment, in the undirected graph of the interference baseline network constructed, each graph node corresponds to an interferometric image pair, and each interferometric image pair has its corresponding acquisition time. They are sorted in the order of the acquisition time before and after. The graph node associated with the interferometric image pair with the earliest acquisition time is marked as the first node. Starting from the first node, the average weight of all the edges connected to this node is used as the threshold, and all the edges below the threshold on this node are removed for optimization and screening.

[0072] In this embodiment, each edge of the undirected graph of the interference baseline network corresponds to a preset weight value, and this preset weight value is the coherence coefficient of the graph node.

[0073] The coherence coefficient is a standardized covariance function, and its value range is [0, 1]. 0 indicates complete incoherence, and 1 indicates complete coherence. The larger the coherence coefficient, the higher the coherence. The echo signals between two SAR images are more similar, and the interferometric image pairs in the formed interferogram can accurately reflect the distance difference between the two echoes, and more accurate surface deformation information can be retrieved.

[0074] Optionally, coherence reflects the similarity degree between two radar echo signals. The coherence r of two zero-mean circular Gaussian complex random signals y1 and y2 can be defined as:

[0075]

[0076] Under ideal conditions, the mathematical expectation E{} in the above formula can be obtained by calculating the ensemble average for each pixel of a large number of interferometric image pairs acquired simultaneously under the same conditions. However, in actual situations, this is usually difficult to achieve. Therefore, it is assumed that in a window of N pixels, the random process is stationary and ergodic. At this time, the spatial average of these N pixels can replace the ensemble average to obtain the coherence coefficient :

[0077]

[0078] The average of the correlation coefficients corresponding to each edge is calculated to obtain the threshold 。

[0079] S30. The remaining interferometric image pairs in the undirected graph of the interference baseline network after removing the target edges are used as the screened target interferometric image pairs.

[0080] As an example, referring respectively to Figure 4 and Figure 5 which show the schematic diagram of the undirected graph of the interference baseline network before removing the target edges and the schematic diagram of the undirected graph of the interference baseline network after removing the target edges. Starting from node A, there are 6 edges connected to node A, and the weights of each edge are C AB =0.45, CAC = 0.22, C AD = 0.27, C AE = 0.39, C AF = 0.42, C AG = 0.28, calculated C avg = 0.29. Therefore, sides AC, AD, and AG are excluded, that is, three interfering image pairs AC, AD, and AG are excluded. The remaining AB, AE, and AF are all used as the target interfering image pairs after screening, and then enter the next node B until the last node G ends.

[0081] It should be noted that in the process of screening interfering image pairs involved in this embodiment, there is no exclusion of nodes (that is, the SAR images themselves are not excluded), and what is excluded are the low-quality interfering image pairs in each interfering image pair that interferes with the node, ensuring the integrity of the interferometric baseline network.

[0082] In the technical solution provided in this embodiment, the connection line of the interfering image pair is used as the edge of the graph, and the coherence coefficient of the interfering image pair is used as the weight of the edge. Starting from the first node, the average weight of all the edges connected to the node is used as the threshold, and all the edges below the threshold on the node are excluded until the last node ends. In the method of this application, each node has an edge connected, which not only excludes the low-quality interfering image pairs but also ensures the integrity of the interferometric baseline network.

[0083] Second Embodiment

[0084] Based on the first embodiment, in this embodiment, the InSAR technology is similar to geodetic technologies such as the Global Navigation Satellite System (GNSS) and very long baseline interferometry (VLBI). They all use radar signals to obtain information. When the radar signal passes through the atmosphere, a delay phenomenon will occur, resulting in a shift in the electromagnetic wave carrier phase, which is called the atmospheric delay effect.

[0085] According to the atmospheric stratification theory, the InSAR atmospheric delay can be divided into two categories, namely ionospheric atmospheric delay and tropospheric atmospheric delay. Existing research has proved that the phase delay effect of the ionosphere on electromagnetic waves is inversely proportional to the wavelength. For example, the ALOS / PALSAR data in the L band (wavelength 23.6 cm) is affected by a greater phase delay in the ionosphere than the Sentinel-1 data in the C band (wavelength 5.6 cm). In addition, using the multi-temporal InSAR data processing method can also weaken the ionospheric delay effect to a certain extent. Therefore, in this embodiment, the atmospheric delay effect of the troposphere on InSAR is considered, and the ionospheric atmospheric delay effect is not considered.

[0086] Furthermore, the troposphere mainly contains two components: dry air and water vapor. Dry air will contract under atmospheric pressure and temperature changes, thus causing a delay effect on radar signals. However, since the delay effect caused by dry air is relatively stable and changes relatively little, it can generally be ignored. Therefore, water vapor, as the main factor causing atmospheric delay in InSAR, varies greatly and is difficult to predict. Water vapor changes over time and space. When radar signals pass through water vapor, refraction occurs, resulting in changes in the signal propagation path and direction, thus generating phase delay errors and affecting the accuracy of InSAR results.

[0087] Therefore, in this embodiment, the GACOS atmospheric model is introduced to correct the atmospheric delay phase in the data. Specifically as follows:

[0088] S40, calculate the atmospheric delay phase value in the initial interferogram pair;

[0089] S50, subtract the atmospheric delay phase value from the initial interferogram to obtain the interferogram pair for constructing the undirected graph of the interference baseline network.

[0090] Furthermore, in S40, the calculation expression of the atmospheric delay phase value is:

[0091]

[0092] In the formula, ZTDk is the atmospheric delay phase value, representing the total zenith delay of the vertical stratification component and the horizontal turbulence component set, k is the coordinate position; T represents the turbulence signal, x k is the site coordinate vector in the local geocentric coordinate system; L0 represents the stratification component delay at sea level; represents the remaining unmodeled residuals, including unmodeled stratification and turbulence signals; is the height scale;

[0093] Among them, the height scale The calculation expression of is:

[0094]

[0095] In the formula, represents the height scale, represents the lowest altitude, represents the highest altitude.

[0096] It should be noted that not all GACOS products can effectively correct the atmospheric delay error of each interferometric image pair. On the contrary, additional noise may be introduced during GACOS correction. This is because GACOS relies on meteorological data from weather forecasting models. If there is intense atmospheric activity during a certain period, the forecasting model may have difficulty capturing all atmospheric changes, resulting in the introduction of new noise in InSAR atmospheric delay correction. Therefore, the standard deviation (STD) of the phase of the interferometric image pair should also be used as an indicator to evaluate the effect of atmospheric correction.

[0097] Exemplarily, in a specific embodiment, referring to Figure 6 the schematic diagram of the standard deviation of the phase before and after GACOS correction shown in the figure, after 435 interferometric image pairs are selected and corrected by GACOS, the STD of 330 interferometric image pairs decreases, accounting for 75.8% of all interferometric image pairs, and the STD of 105 interferometric image pairs increases, accounting for 24.2% of all interferometric image pairs. The maximum increase is 0.08 rad. Generally speaking, GACOS can effectively improve the error caused by atmospheric delay.

[0098] The third embodiment

[0099] In this embodiment, a method for further inverting the surface deformation information of the optimized interferometric image pair is provided, and the steps are as follows:

[0100] Step S100: Obtain the target interferometric image pair obtained after screening based on the screening method of the interferometric image pair.

[0101] Step S200: Select the target interferogram generated by superimposing the SAR images corresponding to two different acquisition times in the target interferometric image pair, and calculate the target interferometric image pair corresponding to the target interferogram. The target interferometric image pair consists of the cumulative deformation amount in the radar line-of-sight direction, the residual topographic phase in the differential interferogram, and the sum of the decorrelation noise.

[0102] Specifically, for the i-th (i = 1, 2,..., N + 1) target interferogram generated by the interference superposition of two SAR images at times TA and TB, the interferometric image pair of the pixel with azimuth coordinate x and range coordinate y can be expressed as:

[0103]

[0104] In the formula, is the cumulative deformation amount in the radar line-of-sight direction, is the residual topographic phase in the differential interferogram, is the atmospheric delay phase, is the decorrelation noise.

[0105] It should be noted that the atmospheric delay phase value removed in the second embodiment is the global atmospheric influence, but there will still be local atmospheric residuals in the interferometric image pair. Therefore, the atmospheric delay phase is still included in the formula of this interferometric image pair. .

[0106] Step S300: Calculate the velocity vector according to the target interferometric image pair and the acquisition time.

[0107] Optionally, in order to obtain time-series deformation information with physical significance, the phase in the above formula is expressed as the product of the velocity vector and time between two acquisition times, that is: and time, namely:

[0108]

[0109] After arrangement, the calculation expression of the velocity vector is:

[0110]

[0111] In the formula, represents the i-th target interferogram, represents the acquisition time corresponding to the i-th target interferogram, represents the (i - 1)-th target interferogram, the acquisition time corresponding to the (i - 1)-th target interferogram.

[0112] Step S400: Calculate the minimum norm solution of the velocity vector by using the singular value decomposition method, and integrate the minimum norm solution of the velocity vector to obtain the deformation amount corresponding to the target interferogram.

[0113] Furthermore, this step specifically includes:

[0114] Step S401: Assume that the interferometric image pair value of the i-th interferogram is expressed as:

[0115]

[0116] In the formula, T A represents the acquisition time A, T B represents the acquisition time B, and v p represents the surface deformation rate.

[0117] Step S402: Represent the integral of each time period over the time interval between the master and slave images by an M×N matrix:

[0118]

[0119] Step S403: Calculate the generalized inverse matrix of matrix B using the singular value decomposition method, and calculate the minimum norm solution of the velocity vector according to the generalized inverse matrix;

[0120] Step S404: Integrate the minimum norm solution of the velocity vector to obtain the deformation amount corresponding to the target interferogram.

[0121] In this embodiment, optionally, since the InSAR technology uses the method of side-looking active emission and reception of microwave signals for imaging, the signals emitted by the SAR sensor during side-looking imaging will inevitably be affected by terrain undulations, resulting in differences between the recorded situation in the image and the actual surface situation. This difference is caused by the imaging geometry and is called the geometric distortion of the SAR image, including layover, shadow, and foreshortening. During the imaging process of the SAR sensor, even subtle changes in elevation can cause large-scale distortions in the image. Among them, the effects of the local incidence angle and the incidence angle in the radar line of sight direction are the most significant. Therefore, the deformation amount includes the local incidence angle and the incidence angle in the radar line of sight direction.

[0122] Step S500: Determine the layover area, shadow area, foreshortening area, and non-deformed area in the SAR image according to the numerical interval where the deformation amount is located, and mask the layover area and the shadow area to complete the inversion of surface deformation.

[0123] In this embodiment, the local incidence angle and the incidence angle in the radar line of sight direction are used as the deformation amount, and the SAR area is judged according to the numerical intervals of the two.

[0124] Optionally, let the local incidence angle be , and the incidence angle in the radar line of sight direction be , and the numerical intervals they are in can be as follows:

[0125]

[0126] Exemplarily, in a specific embodiment, the interferometric baseline network composed of 30 Sentinel-1A images from January 10, 2022 to December 24, 2022 in the study area is regarded as an undirected graph. Combining with external DEM data, based on the LSM (Layover and Shadow Map, LSM) algorithm and the R index, the geometric distortion of radar visibility in the study area is identified, as shown in Figure 7 shown.

[0127] From Figure 7It can be seen that the shadow area (Shadow) and layover area (Layover) of the Sentinel-1A ascending orbit satellite in the study area are mainly concentrated in the eastern and northern regions, showing a strip-shaped distribution. The shadow accounts for 0.2% of the entire study area, and the layover accounts for 4.4% of the entire study area. In actual InSAR data processing, the deformation information obtained in the shadow and layover areas is unavailable, while the foreshortening area and the non-deformed area (Good visibility) are available.

[0128] The fourth embodiment

[0129] Based on any of the foregoing embodiments, the effectiveness and accuracy of the screening of the foregoing interferometric image pairs are verified in this embodiment.

[0130] In this embodiment, the traditional fully connected method (without optimizing the selection of interferometric image pairs), the average coherence coefficient threshold method, the small baseline set method (the time baseline is set to 120 days, the spatial baseline is set to 50%, and there are 214 image pairs in total) are compared with the method of the present application. For the connection of interferometric image pairs of the four methods, please refer to Figures 8 - 11 the interferometric image pair connection diagram based on the fully connected method, the interferometric image pair connection diagram based on the average coherence coefficient threshold method, the interferometric image pair connection diagram based on the small baseline set method, and the interferometric image pair connection diagram based on the screening method of the interferometric image pairs proposed in the present application as shown.

[0131] From Figures 8 - 11 It can be seen that the average coherence coefficient threshold method overall improves the coherence of the interferometric baseline network. However, since the study area is a complex mountainous area, the coherence changes significantly with seasons. The coherence of the interferometric image pairs in January, February, March, April, November, and December is relatively high, while the coherence of the interferometric image pairs in other months is relatively low. After overall averaging, the interferometric image pairs in May and June are basically all eliminated, and the interferometric baseline network is disconnected, making it impossible to perform subsequent surface deformation inversion. The method proposed in the present application regards the entire interferometric baseline network as a graph composed of nodes and edges from the perspective of graph theory. Among them, the SAR image acquisition time is used as the nodes of the graph, the connection lines of the interferometric image pairs are used as the edges of the graph, and the coherence coefficient of the interferometric image pairs is used as the weight of the edges. Starting from the first node, the average weight of all the edges connected to this node is used as the threshold, and all the edges below the threshold on this node are eliminated until the last node is reached. In the method proposed in the present application, each node has edges connected, which not only eliminates the low-quality interferometric image pairs but also ensures the integrity of the interferometric baseline network.

[0132] Since the average coherence coefficient threshold method cannot be used for subsequent surface deformation inversion, we continue to conduct a quantitative comparison using the fully connected method, the small baseline set method, and the interferogram optimization selection method based on graph theory. The coherence of surface deformation inversion, the proportion of effective interferograms, and the inversion error are used as indicators for comparison. The InSAR technique captures and inverses surface deformation information by comparing the phase differences of radar images. Coherence quantifies the degree of consistency of phase information between radar images. High coherence of surface deformation inversion means good consistency of phase information and can obtain more accurate surface deformation information. The proportion of effective interferograms refers to the proportion of interferograms participating in the inversion calculation in all interferograms during surface deformation inversion. A higher proportion of effective interferograms means that more interferograms in the monitoring area participate in the inversion calculation, which helps to improve the reliability and accuracy of monitoring. RMSE refers to the error during surface deformation inversion. The lower this value, the higher the model fitting degree and the more accurate the surface deformation information obtained by inversion. The average values of all pixels of the three indicators, namely the coherence of surface deformation inversion, the proportion of effective interferograms, and the inversion error RMSE, for the three methods are shown in Table 1:

[0133] Table 1. Overall average values of the coherence of surface deformation inversion, the proportion of effective interferograms, and RMSE for the three methods

[0134] It is easy to know that the full-connection method has the worst accuracy. The coherence of surface deformation inversion is 0.18, the proportion of effective interferograms is 70.41%, and the RMSE is 3.52 rad. The full-connection method contains a large number of low-quality interferogram pairs, with the lowest coherence and proportion of effective interferograms and the largest inversion error during surface deformation inversion, indicating that more interferogram pairs do not necessarily yield reliable deformation inversion results and increase the computational cost. Thus, it can be seen that the optimized selection of interferogram pairs is crucial for InSAR data processing. Compared with the full-connection method, the coherence of surface deformation inversion using the small baseline subset method is increased by 0.13, the proportion of effective interferograms is increased by 11.01%, and the RMSE is reduced by 0.89 rad, indicating that the small baseline subset method can effectively suppress temporal and spatial decorrelation and improve the accuracy of surface deformation inversion. The optimized selection method of interferogram pairs based on graph theory proposed in this application has the best accuracy. The coherence of surface deformation inversion is 0.40, the proportion of effective interferograms is 88.13%, and the RMSE is only 2.31 rad. Compared with the method without optimized selection of interferogram pairs, the coherence of surface deformation inversion is increased by 0.22, the proportion of effective interferograms is increased by 17.72%, and the RMSE is reduced by 1.21 rad. Compared with the small baseline subset method, the coherence of surface deformation inversion is increased by 0.09, the proportion of effective interferograms is increased by 6.71%, and the RMSE is reduced by 0.32 rad, and the number of interferogram pairs participating in the inversion calculation is reduced by 40 pairs. In summary, the optimized selection method of InSAR interferogram pairs based on graph theory proposed in this application can effectively suppress the influence of decoherence noise on InSAR deformation inversion, select as few high-quality interferogram pairs as possible from a large amount of data, and improve the reliability and accuracy of InSAR surface deformation inversion.

[0135] As an implementation solution, Figure 12 It is a schematic architecture diagram of the hardware operating environment of the computer system involved in the solution of the embodiment of this application.

[0136] As Figure 12As shown in the figure, the computer system may include: a processor 1001, such as a CPU, a memory 1005, a user interface 1003, a network interface 1004, and a communication bus 1002. Among them, the communication bus 1002 is used to realize the connection and communication between these components. The user interface 1003 may include a display screen (Display), an input unit such as a keyboard (Keyboard). Optionally, the user interface 1003 may also include a standard wired interface and a wireless interface. The network interface 1004 may optionally include a standard wired interface and a wireless interface (such as a WI-FI interface). The memory 1005 may be a high-speed RAM memory or a stable memory (non-volatile memory), such as a disk memory. Optionally, the memory 1005 may also be a storage device independent of the aforementioned processor 1001.

[0137] Those skilled in the art can understand that Figure 12 the computer system architecture shown in the figure does not constitute a limitation on the computer system, and may include more or fewer components than shown in the figure, or combine certain components, or have different component arrangements.

[0138] As Figure 12 shown, in the memory 1005 as a storage medium, there may be included an operating system, a network communication module, a user interface module, and a computer program. Among them, the operating system is a program for managing and controlling the hardware and software resources of the computer system, and for running the computer program and other software or programs.

[0139] In Figure 12 the computer system shown in the figure, the user interface 1003 is mainly used to connect to a terminal and perform data communication with the terminal; the network interface 1004 is mainly used to connect to a background server and perform data communication with the background server; the processor 1001 may be used to call the computer program stored in the memory 1005.

[0140] In this embodiment, the computer system includes: a memory 1005, a processor 1001, and a computer program stored on the memory and executable on the processor, where:

[0141] When the processor 1001 calls the computer program stored in the memory 1005, it performs the following operations:

[0142] S10, taking the acquisition time of the SAR images in the interferometric image pair as nodes and the connection lines between the interferometric image pairs as edges, constructing an undirected graph of the interferometric baseline network;

[0143] S20, in the order of the acquisition times, sequentially remove all target edges with values less than a preset threshold on each graph node in the undirected graph of the interference baseline network, where the preset threshold is the mean of the coherence coefficients of all the edges connected to the graph node;

[0144] S30, use the remaining interferometric image pairs in the undirected graph of the interference baseline network after removing the target edges as the screened target interferometric image pairs.

[0145] When the processor 1001 calls the computer program stored in the memory 1005, the following operations are performed:

[0146] S40, calculate the atmospheric delay phase value in the initial interferometric image pair;

[0147] S50, subtract the atmospheric delay phase value from the initial interferometric image to obtain the interferometric image pair for constructing the undirected graph of the interference baseline network.

[0148] When the processor 1001 calls the computer program stored in the memory 1005, the following operations are performed:

[0149] Obtain the target interferometric image pairs obtained after screening based on the screening method of the interferometric image pairs described above;

[0150] Select the target interferogram generated by superimposing the SAR images corresponding to two different acquisition times in the target interferometric image pairs, and calculate the corresponding target interferometric image pairs of the target interferogram, where the target interferometric image pairs are composed of the cumulative deformation amount in the radar line-of-sight direction, the residual topographic phase in the differential interferogram, the atmospheric delay phase, and the decorrelation noise;

[0151] Calculate the velocity vector according to the target interferometric image pairs and the acquisition times;

[0152] Calculate the minimum norm solution of the velocity vector using the singular value decomposition method, and integrate the minimum norm solution of the velocity vector to obtain the deformation amount corresponding to the target interferogram;

[0153] According to the numerical interval where the deformation amount is located, determine the layover area, shadow area, foreshortening area, and non-deformed area in the SAR image, and mask the layover area and the shadow area to complete the surface deformation inversion. When the processor 1001 calls the computer program stored in the memory 1005, the following operations are performed:

[0154] Let the interferometric image pair value of the i-th interferogram be expressed as:

[0155]

[0156] In the formula, T ADenote the acquisition time A, T B Denote the acquisition time B, v p Denote the surface deformation velocity.

[0157] Integrate each time period over the time interval between the master and slave images and represent it as an M×N matrix:

[0158]

[0159] Use the singular value decomposition method to calculate the generalized inverse matrix of matrix B, and calculate the minimum norm solution of the velocity vector according to the generalized inverse matrix;

[0160] Integrate the minimum norm solution of the velocity vector to obtain the deformation amount corresponding to the target interferogram.

[0161] In addition, those of ordinary skill in the art can understand that all or part of the processes in the methods of the above embodiments can be completed by instructing relevant hardware through a computer program. The computer program includes program instructions, and the computer program can be stored in a storage medium, and the storage medium is a computer-readable storage medium. The program instructions are executed by at least one processor in the computer system to implement the process steps of the embodiments of the above methods.

[0162] Therefore, the present application also provides a computer-readable storage medium, the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, it implements each step of the method for screening interference image pairs or the method for inverting surface deformation as described in the above embodiments.

[0163] Among them, the computer-readable storage medium can be various computer-readable storage media such as a USB flash drive, a mobile hard disk, a read-only memory (ROM), a magnetic disk, or an optical disc that can store program codes.

[0164] It should be noted that since the storage medium provided in the embodiments of the present application is the storage medium used for implementing the methods of the embodiments of the present application, based on the methods introduced in the embodiments of the present application, those skilled in the art can understand the specific structure and deformation of the storage medium, so it will not be elaborated here. Any storage medium used in the methods of the embodiments of the present application belongs to the scope to be protected by the present application.

[0165] Those skilled in the art should understand that the embodiments of the present application can be provided as a method, a system or a computer program product. Therefore, the present application can take the form of a complete hardware embodiment, a complete software embodiment or an embodiment combining software and hardware aspects. Moreover, the present 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.

[0166] The present application is described with reference to the flowcharts and / or block diagrams of methods, devices (systems) and computer program products according to the embodiments of the present application. It should be understood that each flow and / or block in the flowchart and / or block diagram can be implemented by computer program instructions, and the combination of flows and / or blocks in the flowchart and / or block diagram can also be implemented. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor or other programmable data processing devices to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing devices generate means for implementing the specified functions in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0167] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer-readable memory generate a manufactured article including instruction means, and the instruction means implements the specified functions in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0168] These computer program instructions can also be loaded onto a computer or other programmable data processing device, so that a series of operation steps are executed on the computer or other programmable device to generate a computer-implemented process, and thus the instructions executed on the computer or other programmable device provide steps for implementing the specified functions in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0169] It should be noted that in the claims, any reference signs placed between parentheses shall not be construed as limiting the claim. The word "comprising" does not exclude the presence of elements or steps not listed in the claim. The word "a" or "an" preceding an element does not exclude the presence of a plurality of such elements. The present application may be implemented by means of hardware comprising several distinct elements, and by means of a suitably programmed computer. In a unit claim listing several means, several of these means may be embodied by the same item of hardware. The use of the words first, second, and third, etc. does not denote any order. These words may be interpreted as names.

[0170] Although the preferred embodiments of the present application have been described, additional changes and modifications can be made by those skilled in the art once they learn the basic inventive concept. Therefore, the appended claims are intended to be construed to include the preferred embodiments as well as all changes and modifications falling within the scope of the present application.

[0171] Obviously, those skilled in the art can make various changes and modifications to the present application without departing from the spirit and scope of the present application. Thus, if these modifications and variations of the present application fall within the scope of the claims of the present application and their equivalent technologies, the present application is also intended to include these changes and modifications.

Claims

1. A method for screening interference image pairs, characterized in that: The method comprises the following steps: S10, using the acquisition time of SAR images as nodes and the lines between interferometric image pairs as edges, an undirected graph of the interferometric baseline network is constructed; S20, in the order of the acquisition moments, sequentially removing all target edges on each graph node in the interference baseline network undirected graph that are smaller than a preset threshold, wherein the preset threshold is the mean value of the coherence coefficients of all edges connected to the graph nodes; S30, taking the remaining interference image pairs in the interference baseline network undirected graph after removing the target edge as the screened target interference image pairs.

2. The method for screening interference image pairs according to claim 1, characterized in that: Before S10, the following steps are also included: S40, calculating the atmospheric delay phase value in the initial interferometric image pair; S50, subtracting the atmospheric delay phase value from the initial interference image to obtain the interference image pair used as an undirected graph for constructing an interference baseline network.

3. The method for screening interference image pairs according to claim 2, characterized in that: In S40, the calculation expression of the atmospheric delay phase value is: ； Where ZTDk is the atmospheric delay phase value, which represents the total zenith delay of the vertical stratified component and the horizontal turbulent component, k is the coordinate position; T represents the turbulence signal, x k is the site coordinate vector in the local geocentric coordinate system; L0 represents the layered component delay at sea level; represents the remaining unmodeled residual, including the unmodeled stratification and turbulence signals; is a height scale; Among them, the height scale The calculation expression is: ； In the formula, Indicates the height scale, Indicates the minimum altitude, Indicates the maximum altitude.

4. The method for screening interference image pairs according to claim 1, characterized in that: The calculation expression of the coherence coefficient is: ； In the formula, represents the amplitude and phase of the main image, represents the amplitude and phase of the auxiliary image, express The complex conjugate of is used to calculate the coherence, and N represents the total number of samples.

5. A surface deformation inversion method, characterized in that: The surface deformation inversion method comprises the following steps: Acquire a target interference image pair obtained after screening based on the interference image pair screening method according to claim 1; Selecting a target interference image generated by superimposing two SAR images corresponding to different acquisition times in the target interference image pair, and calculating a target interference image pair corresponding to the target interference image pair, wherein the target interference image pair is composed of the cumulative deformation amount in the radar line of sight direction, the residual terrain phase in the differential interference image, and the sum of the decorrelation noise; According to the target interference image pair and acquisition time, the velocity vector is calculated. The minimum norm solution of the velocity vector is calculated by using a singular value decomposition method, and the minimum norm solution of the velocity vector is integrated to obtain a deformation variable corresponding to the target interference pattern; According to the numerical range of the deformation variable, the overlapping area, shadow area, perspective contraction area and undeformed area in the SAR image are determined, and the overlapping area and the shadow area are masked to complete the surface deformation inversion.

6. The surface deformation inversion method according to claim 5, characterized in that: The calculation expression of the velocity vector is: ； In the formula, represents the i-th target interference graph, represents the acquisition time corresponding to the i-th target interference pattern, represents the i-1th target interference graph, The acquisition time corresponding to the i-1th target interference pattern.

7. The surface deformation inversion method according to claim 5, characterized in that: The step of using the singular value decomposition method to calculate the minimum norm solution of the velocity vector, integrating the minimum norm solution of the velocity vector, and obtaining the deformation variable corresponding to the target interference pattern comprises: Let the interference image pair value of the i-th interference pattern be The expression is: ； Where, T A Indicates obtaining time A, T B Indicates the acquisition time B, v p Indicates the rate of surface deformation; The integral of each time period over the time interval of the master and slave images is represented by an M×N matrix: ； A singular value decomposition method is used to calculate a generalized inverse matrix of the matrix B, and a minimum norm solution of the velocity vector is calculated according to the generalized inverse matrix; The minimum norm of the velocity vector is solved as an integral to obtain the deformation variable corresponding to the target interference pattern.

8. An application of the interference image pair screening method as claimed in claim 1 in surface deformation inversion.

9. A computer system, characterized in that: The computer system includes: a memory, a processor, and a computer program stored in the memory and executable on the processor. When the computer program is executed by the processor, the steps of the interference image pair screening method described in any one of claims 1 to 4 or the surface deformation inversion method described in any one of claims 5 to 7 are implemented.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, which, when executed by the processor, implements the steps of the method for screening interference image pairs as described in any one of claims 1 to 4 or the method for inverting surface deformation as described in any one of claims 5 to 7.

Images