Multi-frame wide-area InSAR deformation result fusion method

Abstract

The invention relates to the technical field of remote sensing surveying and mapping, in particular to a multi-frame wide-area InSAR deformation result fusion method, which aims at the splicing problem of multi-frame InSAR deformation under a large spatial scale, performs optimization on the basis of a double-layer triangulation network and a local extraction method, and performs splicing according to dot-shaped and planar deformation storage characteristics. LOS deformation of each frame of image is converted into a vertical direction; and then on the basis of the double-layer triangulation network model, a genetic algorithm and PSO optimization are used as an improved algorithm for dotted deformation data homonymy point pairing to solve the many-for-one condition during registration. And under the condition of ensuring that the space position offset of the overlapping area is small, the offset of adjacent planar data is calculated by utilizing an optimized local extraction method, so that the deviation of adjacent deformation results in a large range is solved. And finally, generating a wide-area InSAR deformation product with consistent space reference and high availability based on two optimization algorithms.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing mapping, and particularly relates to a method for fusing multi-frame wide-area InSAR deformation results. Background Art

[0002] With the improvement of InSAR technology, it is no longer limited to small-scale and accurate deformation information. Wide-area InSAR monitors large-scale regions, meeting the large-scale deformation requirements of scholars and managers for full coverage of provincial-level regions with the same spatio-temporal reference. In InSAR deformation measurement, the deformation in the corresponding LOS direction is obtained according to the image position, but there are connection deviations in the data from different orbits.

[0003] At the same time, due to the limitation of the SAR image format, wide-area InSAR faces the problem of how to splice multi-frame InSAR deformation results. As the processing area increases, the change rates of the overlapping areas of adjacent images are inconsistent due to the calculation reference points, algorithms, error sizes and distributions for obtaining phase information, and directly splicing will result in large deformation deviations in the overlapping area.

[0004] To ensure the usability of the splicing result and the unity of the adjacent orbit data reference, GNSS data calibration plays a crucial role. However, due to the uneven spatial distribution of GNSS stations, limited openness, and high construction and maintenance costs, the feasibility of GNSS data correction is restricted. A variety of splicing methods have been proposed in the prior art, but there are still technical problems such as low accuracy of the deformation result fusion, high requirements for computer equipment configuration, and low calculation efficiency. Summary of the Invention

[0005] The purpose of the present invention is to provide a method for fusing multi-frame wide-area InSAR deformation results. For point data, a double-layer triangular network method using a genetic algorithm and PSO optimization to extract homologous points is adopted. For planar data, homologous points are judged in space using a distance threshold, so as to solve the technical problems of low fusion accuracy and low efficiency in the existing multi-frame InSAR deformation result splicing methods.

[0006] To achieve the above purpose, the present invention provides a method for fusing multi-frame wide-area InSAR deformation results, including the following steps:

[0007] Step 1: Use SBASInSAR technology to obtain multi-frame deformation results and perform data preprocessing;

[0008] Step 2: Extract homologous points in the overlapping area of multi-frame deformations;

[0009] Step 3: Construct a double-layer triangular network based on wide-area InSAR deformation points;

[0010] Step 4: Use the local extraction method to correct the reference offset of wide-area InSAR planar data;

[0011] Step 5: Wide-area InSAR deformation stitching;

[0012] Step 6: Conduct multi-frame result comparison and gross error detection.

[0013] Optionally, the expression of the differential phase between pixels in Step 1 is as follows:

[0014]

[0015] where is the flat-earth phase component related to distance; is the topographic phase; is the component caused by the topographic displacement in the distance direction between two SAR acquisitions; is the phase caused by atmospheric disturbances; includes degradation factors related to temporal and spatial decorrelation and thermal noise.

[0016] Optionally, before extracting corresponding points in Step 2, the LOS-direction deformation needs to be converted to the vertical direction, and then the deformation data in the overlapping area is extracted according to the topological relationship of the cross-track deformation results. After that, a representative deformation area in the overlapping area is selected as the candidate area for corresponding points;

[0017] In Step 2, the genetic algorithm and PSO optimization are used to improve the extraction of corresponding points for point-like data, and for area-like data, the method of distance threshold is used to judge corresponding points in space.

[0018] Optionally, in Step 2, point pairs with a spatial distance less than 1.2×10 -4 m, the minimum pixel size, are defined as corresponding point pairs, and the calculation formula for the spatial distance is as follows:

[0019]

[0020] where d represents the spatial distance, and (x1, y1) and (x2, y2) represent the spatial coordinates of the deformation points on the master image and the registered image, respectively.

[0021] Optionally, the process of using the genetic algorithm for optimization in Step 2 has six specific implementation steps, namely population initialization, fitness evaluation, roulette wheel selection, crossover, mutation, population update and iteration.

[0022] Optionally, during the process of Step 2, the PSO optimization uses the same fitness function as the genetic algorithm.

[0023] Optionally, during the execution of Step 3, first, the spatial reference correction of the longitude and latitude differences of the corresponding points needs to be performed, and then the deformation rate difference of the corresponding points is calculated and the InSAR result reference offset value is calculated statistically.

[0024] Optionally, the reference deviation in step 4 is determined according to the numerical distribution of the InSAR deformation data. Specifically, first, the proportion of -10 to 10 mm / year in the entire image is statistically calculated. The representative value of the local window is the mode position as the median value, and half of the proportion is taken on both sides, and finally the proportion of the stable interval is formed. The representative value of the local window is the average value of the numerical values within the proportion range.

[0025] Optionally, the formula for calculating the reference offset of the main reference and the registration reference in step 5 is as follows:

[0026]

[0027] where, Δv off represents the reference offset, v mi is the average deformation rate of the coherent target in the main reference Frame; v si is the average deformation rate of the coherent target in the Frame to be reference-corrected; is the reference correction amount, and N represents the number of pairs of homologous points.

[0028] The present invention provides a method for fusing multi-frame wide-area InSAR deformation results. Aiming at the splicing problem of multi-frame InSAR deformation at a large spatial scale, it is optimized on the basis of a double-layer triangular network and a local extraction method, and spliced respectively according to the storage characteristics of point-like and surface-like deformations; first, the LOS deformation of each frame of image is converted into the vertical direction; then, on the basis of the double-layer triangular network model, the genetic algorithm and PSO optimization are used as improved algorithms for pairing homologous points of point-like deformation data to solve the one-to-many situation during registration. While ensuring that the spatial position offset in the overlapping area is small, the offset of adjacent surface-like data is calculated by using the optimized local extraction method to solve the deviation of adjacent deformation results in a large range. Finally, based on the two optimization algorithms, a wide-area InSAR deformation product with consistent spatial reference and strong usability can be generated. Through experimental verification, the present invention reduces the operation cost and can stably obtain the wide-area InSAR deformation fusion result under computers with different configurations. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for use in the description of the embodiments or the prior art. Obviously, the following drawings are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0030] Figure 1 It is a schematic flow chart of the specific steps of a method for fusing multi-frame wide-area InSAR deformation results of the present invention.

[0031] Figure 2 It is a schematic diagram of the research and monitoring area selected in the specific embodiment of the present invention.

[0032] Figure 3 It is a schematic diagram of the spatial correspondence relationship of the image overlapping area in the research and monitoring area of the specific embodiment of the present invention.

[0033] Figure 4 It is a schematic diagram of the corresponding relationship of the same-name points of the double-layer triangular network of a multi-frame wide-area InSAR deformation result fusion method of the present invention.

[0034] Figure 5 It is a schematic diagram of the matching effect of the same-name points in the specific embodiment of the present invention.

[0035] Figure 6 It is a schematic diagram of the comparison of deformation field information in the specific embodiment of the present invention. Detailed implementation manners

[0036] The embodiments of the present invention will be described in detail below. Examples of the embodiments are shown in the accompanying drawings, where the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below by referring to the accompanying drawings are exemplary and are intended to explain the present invention, and should not be construed as a limitation to the present invention.

[0037] The following explains some English term abbreviations used in the present invention:

[0038] SBASInSAR: Small Baseline Subset Interferometric Synthetic Aperture Radar, small baseline subset interferometric synthetic aperture radar;

[0039] LOS: Line-of-sight, along the radar line-of-sight direction;

[0040] The present invention provides a multi-frame wide-area InSAR deformation result fusion method, including the following steps:

[0041] Step 1: Use SBASInSAR technology to obtain multi-frame deformation results and perform data preprocessing;

[0042] Step 2: Extract the same-name points in the multi-frame deformation overlapping area;

[0043] Step 3: Construct a double-layer triangular network based on the wide-area InSAR deformation points;

[0044] Step 4: Use the local extraction method to correct the reference offset of the wide-area InSAR planar data;

[0045] Step 5: Wide-area InSAR deformation splicing;

[0046] Step 6: Perform multi-frame result comparison and gross error detection.

[0047] The specific workflow is as Figure 1 shown, and further explanations will be given below in combination with specific embodiments and implementation steps:

[0048] In this embodiment, it is applied to the deformation monitoring of some areas in the North China Plain from 2019 to 2021. For this purpose, four frames of Sentinel-1A SAR images (2019 - 2021) were collected, with a total of two adjacent orbits. Based on the SARscape platform and SBASInSAR technology, the deformation results of each Frame in the LOS direction were obtained in sequence.

[0049] The monitored research area is as Figure 2 shown.

[0050] Step 1: Extract deformation information of multi-frame images by SBASInSAR technology.

[0051] In order to extract the deformation information of the North China Plain, the present invention obtained Sentinel-1A data (C-band) of ascending orbits from 2019 to 2021 from the European Space Agency respectively. The IW interferometric wide swath was default, with a width of about 250 km for each image, a ground resolution of 5×20 m, and the image range included Frame121, 124, 126, and 129. The SBASInSAR technology was used to obtain the surface deformation information including Beijing and Tianjin. The SBAS-InSAR technology uses a multi-master image interference mode to obtain deformation information, effectively reducing the influence of spatio-temporal decorrelation on the solution, and at the same time can also obtain point-like and surface-like deformation results. This method obtains the surface deformation information through the coherence signal difference of SAR images, uses short spatial baseline differential interferograms to overcome spatial decorrelation, and retains the monitoring points with high correlation. The original differential interferometric phase contains multiple phase components, and the differential phase between adjacent pixels can be expressed by Equation (1). In order to separate the deformation phase, this embodiment collected 30m DEM data and precise orbit data to ensure that the influence of terrain and orbit errors can be largely eliminated during the deformation solution process.

[0052]

[0053] is the flat-earth phase component related to distance; is the terrain phase; is the component caused by terrain displacement in the range direction between two SAR acquisitions; is the phase caused by atmospheric disturbances; includes degradation factors related to spatio-temporal decorrelation and thermal noise. The specific information of the data used in this experiment is shown in Table 1.

[0054] In the data preprocessing stage, it is found that the range of a single frame of image is too large for the extraction model, and the process of the deformation extraction model to identify stable points takes a long time and the recognition accuracy will decrease. Due to the limitations of computer equipment, in order to ensure the acquisition of relatively high-precision deformation information, according to the characteristic that the parameter information of the same orbit is the same, the present invention seamlessly stitches two frames of SAR images in data preprocessing. Then, the stitched data of the same orbit is cropped and batch-processed, improving the calculation efficiency.

[0055] Table 1 Sentinel-1A data information

[0056]

[0057] Step 2: Extracting corresponding points in the multi-frame deformation overlapping area.

[0058] Extracting corresponding points in the overlapping area is an important step in constructing a triangular network. However, the inconsistent satellite incident angles and azimuth angles of adjacent-track images will affect the calculation of the reference deviation. Before extracting corresponding points, it is necessary to unify the projection directions of the deformation rates of each frame of image. The InSAR technology processing projects the three-dimensional deformation information extracted from the phase information onto the one-dimensional deformation along the LOS direction. However, the observation angles and line-of-sight distances of adjacent images result in different deformations and cannot be directly fused. Due to the characteristics that the SAR satellite flies along the polar orbit and performs periodic repeated observations of ground objects and is insensitive to the north-south deformation component, it is found that the main contribution to the deformation along the LOS direction is the vertical direction, followed by the east-west direction. Since the present invention uses only the SAR image data of the ascending orbit and there is no support from ground measured data, the east-west deformation component is selected to be ignored, simplifying the influence of the satellite incident angle on the difference of InSAR deformation results, and converting the deformation along the LOS direction into the vertical direction.

[0059]

[0060] where θ and α respectively represent the satellite incident angle and azimuth angle, d U 、d LOS are respectively the deformation in the vertical direction (U) and the deformation along the LOS direction, θ is the satellite incident angle, d E and d N are respectively the east-west and north-south deformation components.

[0061] To improve the efficiency of homologous point extraction, deformation data in the overlapping area is extracted according to the topological relationship of adjacent track deformation results. Then, representative deformation areas in the overlapping area are selected as candidate areas for homologous points, and the selection conditions are based on data deformation distribution and ground object characteristics. It should be noted that most of the deformations are stable and slight. However, the overlapping area is usually located in the edge area of the image, and is relatively highly affected by incoherent noise and unwrapping errors during the solution process. Therefore, selecting representative overlapping areas with less stable and deformation factors as candidate areas for homologous points according to the deformation distribution information and ground object characteristics of the overlapping area can ensure the amplitude of deformation fluctuations. Study the spatial correspondence relationship of the image overlapping area in the monitoring area as Figure 3 shown

[0062] In the candidate area for homologous points, the homologous points are extracted using distance threshold, genetic algorithm, and PSO optimization respectively. The distance threshold is based on the principle that the spatial positions of homologous points are theoretically the same. However, due to differences in satellite flight angles, observation times, and solution differences, there are distance differences between homologous points. The distance threshold formula judges the spatial distance between point pairs according to the Euclidean distance principle (Formula 4), and defines point pairs with a spatial distance less than the minimum pixel size (1.2×10 -4 m) as homologous point pairs

[0063]

[0064] d represents the spatial distance, and (x1, y1) and (x2, y2) represent the spatial coordinates of the deformation points on the master image and the registered image respectively

[0065] However, only distance information is considered in homologous point pairing. In areas where deformation points are dense, homologous point pairs may form a one-to-many or many-to-one situation, thus affecting the uniqueness and stability of the triangulation network solution. To avoid this situation, the global search advantage of the genetic algorithm (Genetic Algorithm, GA) is utilized to find relatively suitable homologous point pairs through fitness evaluation, increasing the fault tolerance of point pairing. The genetic algorithm optimizes problems by simulating biological evolution principles such as natural selection, crossover, and mutation. There are six specific implementation steps, namely initializing the population, fitness evaluation, selection (roulette wheel selection method), crossover, mutation, population update and iteration. Number the deformation data points in areas A and B of the overlapping area, calculate A(x1, y1) and B(x2, y2) through the distance formula (4). For each individual in the population, calculate the distance d ij of each pair of paired points, where i and j are determined according to the number of deformation points

[0066] f(x) = object(x)(5)

[0067]

[0068] Among them, formula (5) is the fitness value formula, formula (6) is the roulette wheel selection method, and P i represents probability; formula (7) is order crossover, Offspring represents the crossover offspring under a certain probability, Parent represents the population interval, and based on a certain crossover probability, those below the probability are selected for crossover. The value unit of c depends on the number of point pairs; formula (8) is swap mutation, which determines whether the pairing of individuals mutates with a certain probability; finally, a new generation of population is generated, the number of iterations is set, and the steps of fitness evaluation, selection, crossover, mutation, and population update are repeated.

[0069] After extraction using the genetic algorithm, PSO optimization is also added. The core of its principle is to find the optimal solution to the problem through group collaboration and information sharing. The quality of particles is evaluated from a group of particles. The particle position is updated according to its own optimal position (referred to as pbest) and the global optimal position (referred to as gbest). Whenever it moves, the fitness is recalculated, and it is updated when it is in the local optimal or global optimal position. Iterating in this way, it gradually approaches the global optimum. This greatly improves the extraction of one-to-one corresponding points and largely avoids the occurrence of one-to-many point pairs. PSO optimization also uses the sum of the distances between paired points as the fitness function and sets the initial velocity to 0.

[0070] v i (t + 1) = wv i (t) + c1r1(pbest i - x i (t)) + c2r2(gbest - x i (t)) (9)

[0071] x i (t + 1) = x i (t) + v i (t + 1) (10)

[0072]

[0073] In formulas 1 - 9, v i (t) is the velocity of particle i at time t, w is the inertia weight used to balance the global and local search capabilities, c1 and c2 are learning factors, r1 and r2 are random numbers between 0 and 1, pbest i is the individual best position of particle i up to time t, gbest is the global best position of the entire particle swarm up to time t, and x i (t) is the position of particle i at time t. Formula (10) is to update the position of the particle with the updated velocity; formula (11) is the inertia weight adjustment.

[0074] The following is illustrated by specific examples. Example: Assume that in the same geographic coordinate system, the numbers of InSAR deformation points in blocks A and B are 180 and 150 respectively, and the information contained in the deformation points includes deformation rate, longitude and latitude, coherence, solution accuracy, etc. First, number the deformation points in blocks A and B. A = [0, 1, …, 179], B = [0, 1, …, 149]; then calculate the distances between the first point in A and the deformation points in B through the Euclidean distance formula, and each individual is denoted as x i = [x i0 , x i2 , … x i149 . And so on. The pairs of points with the minimum distance between them form pairs of points, and 150 pairs of points will appear. Denote the minimum distance of all pairs of points as d ij , j = [0, 1, …, 149]. However, if a pair forms a pair with two or more points, the number of pairs of points will be greater than 150. The sum of the fitness functions is Then calculate the probability i that the individual x is selected. After that, with a certain crossover probability, randomly select two crossover points c1 and c2, where 0 < c1 < c2 < 149, and retain the middle section from c1 to c2 of x i in block A, denoted as I A . Then, starting from c2 + 1, take numbers in sequence for x i in block B, skipping the existing numbers (minimum distance values) in I A , and fill them into I A to form the offspring C1. According to the above operation, retain the middle section from c1 to c2 of x i in block B, denoted as I B . Then, starting from c2 + 1, take numbers in sequence for x i in block A, skipping the existing numbers (minimum distance values) in I B , and fill them into I B to form the offspring C2. First, judge whether an individual needs to mutate according to a certain probability. If it needs to mutate, search for the position of each individual x i . For the position j, if the generated number is less than the probability, randomly select a position k ≠ j, and swap the positions of x ij and x ik . Put the individuals after crossover and mutation into the population to form a new population. Set the number of iterations and repeat the above steps. However, to prevent the situation of one point corresponding to multiple points from still occurring after the genetic algorithm, use PSO optimization to process the local area. Use the population finally generated by the genetic algorithm as the initial particle swarm for PSO optimization. To ensure the consistency of evaluation, PSO optimization uses the same fitness function as the genetic algorithm. Set the inertia weight w and learning factors c1, c2. For the particle x iRandomly generate r1 and r2 (the particle can be assumed to be the 10th point of block A corresponding to the 12th and 20th points of block B). At the initial stage of the algorithm, the individual best historical position pbest i and the global best historical position gbest i are set based on the initial positions of the particles. Calculate the velocity through formula (9), and combine formula (10) to calculate the position, and substitute it into the fitness formula to obtain the fitness value f(x i ). For the individual extreme value, if f PSO (x i ) < f PSO (pbest i ), then update pbest i+1 = x i ; for the global extreme value, it is necessary to compare the f(x i ) of all particles, find the particle with the minimum fitness. If f PSO (x i ) < f PSO (gbest i ), then update gbest i+1 = x i+1 . The termination condition is the maximum number of iterations. Using genetic algorithm and PSO optimization can pair the same-name points one-to-one to the greatest extent.

[0075] To ensure the one-to-one correspondence of the same-name points and the uniqueness of the triangular network, the present invention sets two barriers. One is based on the calculation accuracy of the deformation points. Assuming that there are 180 points in block A and 150 points in block B, a calculation accuracy matrix of 180×150 is obtained. If the 5th point in block A of the point pair (individual or particle) corresponds to the 6th and 10th points in block B, then obtain P 56 and P 510 's accuracy from the accuracy matrix, and finally retain the point pair with higher accuracy. If the accuracies are the same, then delete the pairing opportunity of the 5th point in block A and block B.

[0076] Step 3: Construct a double-layer triangular network based on the wide-area InSAR deformation points.

[0077] After extracting the same-name points, the spatial coordinate correction and deviation calculation of them are the key to ensuring the adjacent-track image stitching. Calculate the offset between the main Frame and the secondary Frame by constructing a double-layer triangular network. The double-layer triangular network is composed of first establishing triangular networks by the deformation points of each layer respectively, and counting the same-name point pairs in the adjacent-track or same-track overlapping area, which is to ensure the uniqueness of the upper and lower triangular networks (please refer to Figure 4)。First, it is necessary to perform spatial reference correction on the longitude and latitude differences of homologous points, recalculate the longitude and latitude coordinates of homologous points as shown in formula (12), so that the coordinates of homologous points are consistent to achieve the effect of distributed coherent target fusion. Then, by establishing a network of homologous points in the "vertical direction", calculate the deformation rate difference of homologous points and calculate the reference offset value of the InSAR result statistically.

[0078] There are three purposes for the triangular network: 1) Calculate the offset value between the secondary registered image and the primary image; 2) Correct the spatial reference of the primary and secondary images; 3) Detect the accuracy of the triangular network solution of the primary and secondary images.

[0079]

[0080] Formula (12) is used for longitude and latitude correction after confirming homologous points, where (lon A , lat A ) represents the longitude and latitude coordinates of point A, and (lon B , lat B ) represents the longitude and latitude coordinates of point B.

[0081] Step 4: The local extraction method is used to correct the reference offset of the wide-area InSAR areal data.

[0082] In the case of a large amount of point data, the computer processing process for multi-track data splicing is cumbersome and time-consuming. In order to improve the calculation efficiency of the offset, the multi-frame deformation rate splicing method (MFDVS) has been improved. The principle of the MDFVS method is to use the numerical distribution feature - the mode in the dataset to replace the deformation value of the small window. In the comparison between the mode and the average value, since the deformation is mostly in a stable form, the local average value is easily affected by large deformation values and is not representative for a small local area. The difference between the modes of the overlapping local areas is used as the reference deviation, and then the image is corrected. Therefore, the mean of the modes corresponding to the overlapping areas is selected as the observation of the deviation correction equation. However, the local data distribution mostly shows a normal distribution. If only the number with the most occurrences in the local area is used to represent the area, it is found that the representativeness is not strong either. Ideally, the smaller the window obtained, the more representative the mode is. It is proposed to perform value selection according to the numerical distribution of InSAR deformation data. In the study of InSAR deformation analysis, the deformation rate interval of -10 to 10 mm / year is defined as the stable interval. The stable interval of the deformation values of the entire image is statistically analyzed, and the proportion of the stable interval in the entire deformation result is calculated. Then, according to the proportion of the stable interval, the values in the local window are statistically analyzed. The proportion of the stable interval in the local window is based on the mode position as the middle value, and half of the proportion is taken on both sides, and finally the proportion of the stable interval is formed. The representative value of the local window is the average of the values within the proportion range. Taking the vertical deformation as an example, it is assumed that there are M InSAR image datasets in the wide area, and N overlapping areas are generated between them. For the i-th overlapping area, the settlement rate proportions m i and s i composing the overlapping area are statistically analyzed. Theoretically, |m i = s i |. However, there is serious subsidence in some parts of the North China Plain. If there is serious deformation in the overlapping area, due to factors such as the angles and reference points of different images, the deformation difference of the settlement bowl may be too large and is not suitable for estimating the reference offset. Therefore, the seriously deformed areas are excluded, and the overlapping areas are selected according to the stable ground object features.

[0083] Step 5: Wide-area InSAR deformation splicing.

[0084] The genetic algorithm, PSO optimization combined with a double-layer triangular network is used for wide-area InSAR deformation correction. The deformation products of different orbits are spliced through the difference formula (13) between points. The reference offsets of the main reference and the registration reference can be calculated according to the formula.

[0085]

[0086] Among them, Δv off in formula (13) represents the reference offset, and v miis the average deformation rate of coherent targets in the main reference Frame; v si is the average deformation rate of coherent targets in the Frame to be reference-corrected; in formula (14) is the reference correction amount, and N represents the number of pairs of homologous points. Among them, the selection criterion for the main reference Frame is to select the image with more stable ground objects. Stable ground objects such as bare rocks and buildings are not easily deformed.

[0087] Step 6: Comparison of multi-frame results and gross error detection.

[0088] There are differences in shooting time, solution reference points, spatio-temporal changes, etc. among images. However, in terms of the results, the deformation difference in the overlapping area of images in the same orbit is relatively small. The main focus of this embodiment is on the correction between images in adjacent orbits. After the deformation information is projected vertically backward, it is found that there are size differences in the overlapping area. By comparing the deformation results calculated from images in different orbits, it is found that the deformation information in the overlapping area between Frame126 and Frame129 is relatively close, while the gap in the overlapping area between image Frame121 and image Frame124 is relatively large. To detect the existence of gross errors, the same-plane side lengths of the double-layer triangular network and the differences of the double-layer triangular network are detected. The difference in the deformation rate between the upper and lower sides in the triangular network is calculated, and the data characteristics of the overall homologous points are statistically analyzed to judge the reliability of the network solution. The difference between the deformation rate differences of the two triangular network sides is:

[0089] ΔV ms_ij =ΔV m_ij -Δ s_ij (15)

[0090] In formula (15), ΔV ms_ij is the difference in the deformation rate of the triangular network, ΔV m_ij is the difference in the triangular network rate of the deformation result of the main image, and Δ s_ij is the difference in the triangular network rate of the deformation result of the secondary image. The calculated threshold is set to ±10 mm / year. When this value is exceeded, the solution deviation of the double-layer triangular network is relatively large, and the network structure needs to be rechecked and corrected.

[0091] However, the improved local extraction method is based on the information of the deformation data itself for splicing. In order to ensure the standard deviation of the deformation products before and after comparison and splicing and use visual observation of the deformed area after splicing, if the deviation is too large, the data characteristics of the local area are reselected and calculated.

[0092] Furthermore, in this embodiment, a comparative analysis of the results is also carried out:

[0093] 1. Comparison and correction of the overlapping area of the point-like and surface-like deformation results of wide-area InSAR:

[0094] The data of Frame121, 124, 126, and 129 will be used to calculate the deformation of the overlapping area in the Tianjin area, aiming to compare the deformation differences between adjacent-track images. The solution method combining double-layer triangular mesh with genetic algorithm and PSO optimization is more global in extracting corresponding points compared with the distance threshold method, with higher accuracy in finding corresponding points, which is more conducive to establishing a qualified double-layer triangular mesh subsequently. However, due to the large amount of point data of deformation, we need to configure computing devices with higher performance. We also extract corresponding points by only using planar grid data, automatically select the overlapping area through code and crop it. After determining the main image and the secondary image, set the minimum pixel unit value. Here, according to the characteristics of the grid, the minimum pixel unit is set to 1.2×10 -4 . Briefly speaking, the deformation result of SBASInSAR calculates the offset value through the coherence between images. Essentially, it calculates the regional deformation value through points with the minimum pixel unit, and finally forms planar data by interpolation during geocoding.

[0095] By calculating the images of Frame121, 124, 126, and 129, we obtain the deformation data of the North China region from 2019 to 2021. Use the improved double-layer triangular mesh for splicing, and the registration order is clockwise. Taking Frame124 as the main image and Frame121 as the registered image, after obtaining the deviation values of Frame121 and 126 images, then taking Frame124 as the main image, register Frame129 into the Frame124 image, and finally obtain a wide-area InSAR deformation field with better fusion results. The severely deformed areas are mainly distributed in Shengfang Town of Bazhou City, the coastal area of Tianjin City, and the northern part of Cangzhou City, and most areas show a stable state. The settlement is mainly concentrated in the building groups. Due to the existence of a large amount of farmland in the study area, during the planting period, the satellite image signal encounters low reflection signals in water areas, and the deformation information gaps are mainly concentrated on the farmland. The water resources utilization from 2019 to 2021 mainly comes from the South-to-North Water Diversion Project, surface water, and groundwater. The geological reason for the settlement in the study area is that the Quaternary sedimentary layer is widely distributed, and the loose layer is distributed throughout the study area. However, pumping groundwater for domestic and agricultural water use is the main reason for the settlement. In order to prevent the spread of land subsidence, policies such as groundwater extraction bans and water extraction applications have been promulgated in Beijing, Tianjin and other places. From the monthly cumulative deformation amount, the settlement in the North China Plain shows a decreasing trend, especially during the period from 2020 to 2021.

[0096] 2. Accuracy comparison of deformation products before and after splicing

[0097] Figure 5 It is the effect diagram of corresponding point pairing. During the period from 2019 to 2021, for the point data, the deformation rate of the experimental points did not exceed ΔV during corresponding point pairing ms_ij= 10 mm / year. Before calibration, the root mean square error of the corresponding points of the main image and the secondary image in the selected area was ±8.87 mm / year. After calibration, the root mean square error between the corresponding points was only ±6.61 mm / year. This is only the accuracy improvement in the local experimental area. It is believed that more local areas can be selected to obtain higher accuracy results. For the areal data, the overall accuracy (standard deviation) of the Sentinel-1A dataset before and after calibration was 10.8 mm / year and 6.8 mm / year respectively. After calibration, the accuracy was improved by 37%, which was closer to the standard deviation of 6.78 mm / year of the main registration image. After fusion, the standard deviation was 7.0 mm / year.

[0098] 3. Multi-frame deformation rate analysis of Sentinel-1

[0099] Please refer to Figure 6 , Figure 6 for the deformation field information, where a and b are the deformation rates of path40 and path142 respectively, and c is the wide-area InSAR deformation product after offset correction and fusion.

[0100] In summary, the present invention improves the double-layer triangular network and the local extraction method to perform offset correction on the point-like and areal deformation data respectively. The two methods analyze the characteristics and monitoring features of adjacent-track images, construct a double-layer triangular monitoring network for extracting corresponding points and correcting the spatial positions of point-like data. To reduce the computational cost, a 30×30 grid local extraction method is used for areal data, and the characteristic values of the local area are selected to perform offset correction on the deformation of the secondary image. After verification by the embodiments, the improved double-layer triangular network and the improved local extraction method of the present invention respectively correct the InSAR products for the point-like and areal deformation results in a large range, and the corrected results can be seamlessly stitched.

[0101] The above-disclosed are only one or more preferred embodiments of the present invention. Of course, the scope of the rights of the present invention cannot be limited thereby. Those of ordinary skill in the art can understand all or part of the processes of implementing the above embodiments, and the equivalent changes made according to the claims of the present invention still fall within the scope covered by the present invention.

Claims

1. A multi-frame wide-area InSAR deformation result fusion method, characterized in that It includes the following steps: Step 1: Obtain multi-frame deformation results using SBAS InSAR technology and perform data preprocessing; Step 2: Extract homologous points in the overlapping area of multi-frame deformations; Step 3: Construct a double-layer triangular network based on wide-area InSAR deformation points; Step 4: Use the local extraction method to correct the reference offset of wide-area InSAR areal data; Step 5: Wide-area InSAR deformation stitching; Step 6: Conduct multi-frame result comparison and gross error detection.

2. The multi-frame wide-area InSAR deformation result fusion method according to claim 1, wherein The expression of the differential phase between images in Step 1 is as follows: Among them, is the flat ground phase component related to distance; is the terrain phase; is the component caused by terrain displacement in the range direction between two SAR acquisitions; is the phase caused by atmospheric disturbance; includes degradation factors related to time and spatial decorrelation and thermal noise.

3. The multi-frame wide-area InSAR deformation result fusion method according to claim 2, wherein Before extracting homologous points in Step 2, it is necessary to convert the LOS-direction deformation to the vertical direction, then extract the deformation data in the overlapping area according to the topological relationship of the adjacent-track deformation results, and then select a representative deformation area in the overlapping area as the candidate area for homologous points; In Step 2, the genetic algorithm and PSO optimization are used to improve the extraction of homologous points for point-like data, and for areal data, the method of distance threshold is used to judge homologous points in space.

4. The multi-frame wide-area InSAR deformation result fusion method according to claim 3, wherein In step 2, point pairs with a spatial distance less than the minimum pixel size of 1.2×10 -4 m are defined as corresponding point pairs, and the calculation formula for the spatial distance is as follows: Among them, d represents the spatial distance, and (x1, y1) and (x2, y2) respectively represent the spatial coordinates of the deformation points on the master image and the registered image.

5. The multi-frame wide-area InSAR deformation result fusion method according to claim 4, wherein The process of using the genetic algorithm for optimization in Step 2 has six specific implementation steps, namely population initialization, fitness evaluation, roulette wheel selection method, crossover, mutation, population update and iteration.

6. The multi-frame wide-area InSAR deformation result fusion method according to claim 5, wherein During the process of Step 2, the PSO optimization adopts the same fitness function as the genetic algorithm.

7. The multi-frame wide-area InSAR deformation result fusion method according to claim 6, wherein During the execution of Step 3, first, it is necessary to correct the spatial reference of the longitude and latitude difference of homologous points, and then calculate the deformation rate difference of homologous points and calculate the reference offset value of the InSAR result statistically.

8. The multi-frame wide-area InSAR deformation result fusion method according to claim 7, wherein The reference deviation in Step 4 is taken according to the numerical distribution of the InSAR deformation data. Specifically, first, the proportion of -10 to 10 mm / year in the entire image is statistically calculated. The representative value of the local window is the median at the mode position, and half of the proportion is taken on both sides, and finally the proportion of the stable interval is formed. The representative value of the local window is the average value of the numerical values within the proportion range.

9. The multi-frame wide-area InSAR deformation result fusion method according to claim 8, wherein The formula for calculating the reference offset of the master reference and the registered reference in Step 5 is as follows: Among them, Δv off represents the reference offset, and v mi is the average deformation rate of the coherent target in the main reference Frame; v si is the average deformation rate of the coherent target in the Frame to be reference-corrected; is the reference correction amount; N represents the number of pairs of homologous points.

Images