A multi-reference point deformation solution method for time-series InSAR

Abstract

The present application relates to the field of surface deformation monitoring technology and provides a method for multi-reference point deformation solution of time-series InSAR, which includes: performing phase unwrapping on all SAR images to obtain amplitude, coherence, differential phase and unwrapped phase data; using the amplitude, coherence, differential phase and unwrapped phase data to calculate the values ​​of multiple data indicators of pixel points, and calculating the comprehensive score of each pixel point based on the values ​​of all data indicators of each pixel point; sorting all comprehensive scores from large to small, taking the pixel points corresponding to the top comprehensive scores in the sorting results as candidate reference points, and determining stable reference points from all candidate reference points; and obtaining the time-series deformation results of the target area in the target time period under the constraints of multiple reference points based on all reference points. The method of the present application can weaken the influence of atmospheric noise and improve the accuracy of time-series InSAR deformation solution.

Description

Technical Field

[0001] The present application relates to the technical field of surface deformation monitoring, and in particular to a multi-reference point deformation solution method for time-series InSAR. Background Art

[0002] Surface deformation is a slow and irreversible Earth movement phenomenon, usually triggered by natural factors or human activities. When the magnitude of deformation reaches a certain level, it may evolve into a highly destructive geological disaster. Traditional monitoring methods such as total stations, levels, and the Global Navigation Satellite System (GNSS), while excellent in measurement accuracy, have significant limitations in large-scale "dual control" monitoring of potential hazards and risk areas. These limitations include limited monitoring coverage, which makes it difficult to cover large areas. Furthermore, they consume significant manpower and material resources, and their spatial resolution is relatively low. This makes it difficult for ground-based contact point monitoring to fully and promptly reflect potential risks, especially for sudden, less significant geological hazards. Furthermore, traditional contact measurement methods are susceptible to interference from external environmental conditions, making observations significantly more difficult in inaccessible areas or complex terrain.

[0003] Interferometric Synthetic Aperture Radar (InSAR) is an active microwave remote sensing technology that uses the interferometric phase of two Synthetic Aperture Radar (SAR) images to form an interferometer pair, extracting deformation information from the interferometric phase. This technology then simultaneously processes multiple time-interleaved interferometer pairs to obtain temporal deformation information, leading to the development of Multi-Temporal InSAR (MT-InSAR). This technology, which uses SAR satellites to periodically monitor surface deformation, enables all-day, all-weather, large-scale, high-spatial-resolution continuous tracking of the evolution of geological hazard risks. It effectively covers hard-to-reach areas and provides a scientific basis for the delineation of risk zones. Research on the application of MT-InSAR for large-scale geological hazard monitoring not only addresses the shortcomings of traditional monitoring methods but also provides technical support for the establishment of a dual-control and prevention model for geological hazards, further enhancing disaster prevention and mitigation capabilities. This is of great significance for ensuring public safety and preventing and controlling disasters.

[0004] However, during the MT-InSAR time-series deformation solution, the original monitored deformation signal is relative in the spatial dimension. To obtain the absolute value of surface deformation in the spatial dimension, most current methods rely on prior information about the study area. They manually select a stable surface point as the solution reference point and assume that the deformation value at this reference point is zero to obtain the absolute spatial value of the MT-InSAR time-series deformation observation. This method has three major limitations: 1) The reliability of the reference point selection depends largely on the experience of the data processor. Improper reference point selection can easily lead to systematic deviations in the MT-InSAR time-series deformation solution. 2) For the entire study area, the phase continuity assumption cannot be met across the entire image due to decoherence caused by factors such as atmospheric effects or large temporal and spatial baselines. This can lead to unwrapping errors or systematic offsets in areas far from the unwrapping reference point. For large-scale or complex environmental study areas, multiple reference points must be selected for solution correction. 3) Reference points must be manually set based on experience for each study area. In cases where multiple reference points are required for solution, manual selection results in a massive workload and significantly hinders the automation and intelligentization of MT-InSAR time-series deformation processing. To address these limitations, a time-series InSAR deformation monitoring method based on multiple reference point constraints is urgently needed.

[0005] In fact, apart from large-scale deformations induced by events like earthquakes, deformations induced by most other events are localized. This means that within a given monitoring area, the vast majority of surface points are stable. Automatically selecting high-quality interferometric points as reference points for MT-InSAR calculations would overcome the limitations of the single-reference point model and improve the accuracy of MT-InSAR time-series deformation monitoring. However, key technologies such as how to accurately and automatically select high-quality references and construct MT-InSAR calculation models under multi-reference point constraints are currently underdeveloped. This results in low accuracy in time-series InSAR deformation calculations. Summary of the Invention

[0006] The present application provides a multi-reference point deformation solution method for time-series InSAR, which can solve the problem of low accuracy of time-series InSAR deformation solution.

[0007] In a first aspect, an embodiment of the present application provides a multi-reference point deformation solution method for time-series InSAR, the multi-reference point deformation solution method comprising:

[0008] Acquire multiple SAR images of the target area within the target time period, generate differential interferometry pairs based on all the SAR images, perform phase unwrapping on all the differential interferometry pairs, and obtain amplitude, coherence, differential phase, and unwrapped phase data. The earliest SAR image acquired among all the SAR images is the reference SAR image.

[0009] The amplitude, coherence, differential phase, and unwrapped phase data are used to calculate the values ​​of multiple data indicators for each pixel in the reference SAR image. A comprehensive score for each pixel is then calculated based on the values ​​of all the data indicators for each pixel. The data indicators are used to evaluate the quality of the deformation data corresponding to the pixel.

[0010] All comprehensive scores are sorted from large to small, and the pixels corresponding to the top multiple comprehensive scores in the sorting results are selected as candidate reference points. Multiple reference points are determined from all candidate reference points; the stability of the area corresponding to the reference point in the target area is greater than the stability of the areas corresponding to all other candidate reference points;

[0011] The temporal deformation results of the target area in the target time period under the multi-reference point constraint are obtained based on all reference points.

[0012] Optionally, a comprehensive score is calculated for each pixel based on the values ​​of all data indicators for each pixel, including:

[0013] By formula:

[0014] ;

[0015] Calculate the Comprehensive score of pixels ;

[0016] in, Indicates the The scoring coordinate vector corresponding to the pixel point, , Indicates the number of pixels, represents the optimal solution vector, represents the worst solution vector, Indicates the The distance between the pixel point and the worst solution vector, Indicates the The distance between a pixel point and the optimal solution vector:

[0017] ;

[0018] ;

[0019] in, Indicates the The optimal value of a data indicator, Indicates the The worst value of a data indicator, Indicates the weight of each indicator, Indicates the The first pixel The value of the data indicator, Indicates the number of data indicators.

[0020] Optionally, multiple reference points are determined from all candidate reference points, including:

[0021] Divide the target time period into multiple observation groups;

[0022] Calculate statistics for each observation group;

[0023] Determine whether there is an observation group whose statistic is less than or equal to the critical value of the statistic among all observation groups; the observation group includes multiple observation periods;

[0024] If so, the deformation observation period in which the terrain deformation occurs is determined from all observation periods in the observation group that satisfy the statistic less than or equal to the critical value of the statistic, and multiple reference points are determined from all candidate reference points according to the deformation observation period;

[0025] Otherwise, divide each observation group into multiple sub-observation groups, and take all sub-observation groups as multiple observation groups, and return to the step of calculating the statistics for each observation group.

[0026] Optionally, calculate statistics for each observation group, including:

[0027] By formula:

[0028] ;

[0029] ;

[0030] ;

[0031] Calculating statistics ;

[0032] in, represents the sum of the degrees of freedom of all observation periods in the observation group, represents the combined unit weight variance of all observation periods, Indicates the The unit weight variance of the observation data in the observation period is Indicates the degrees of freedom for each observation period, represents the number of observation periods, Expressed as a residual correction vector, represents the observation weight matrix, represents the residual quadratic form.

[0033] Optionally, the deformation observation period in which the terrain deformation occurs is determined from all observation periods in the observation group that satisfy the statistic less than or equal to the statistic critical value, including:

[0034] For every two adjacent observation periods in the observation group, the deformation gap statistics between the two adjacent observation periods are calculated. If the deformation gap statistics is greater than the deformation gap statistics threshold, the two adjacent observation periods are considered to be deformation observation periods.

[0035] Optionally, calculate the deformation gap statistics between two adjacent observation periods, including:

[0036] By formula:

[0037] ;

[0038] ;

[0039] ;

[0040] Calculate the deformation gap statistic between two adjacent observation periods ;

[0041] in, represents the empirical variance between the residuals of two adjacent observation periods, represents the sum of the degrees of freedom of two adjacent observation periods, represents the coordinate difference between two adjacent observation periods, The matrix representing the difference between the observations of two observation periods, The weight matrix representing the difference between the observation results of two observation periods, represents the number of independent observation conditions, Represents the matrix of observation corrections after adjustment, represents the observation weight matrix, represents the residual sum of squares of the first period observations, represents the residual sum of squares of the second period observations, Represents a transpose operation.

[0042] Optionally, multiple reference points are determined from all candidate reference points based on the deformation observation period, including:

[0043] Divide all candidate reference points into stable point sets and unstable point sets;

[0044] Calculate the stability value of the stable point set according to the deformation observation period, the stable point set and the unstable point set;

[0045] Determine whether the stability value is less than the stability threshold;

[0046] If so, all candidate reference points in the stable point set are used as reference points;

[0047] Otherwise, a candidate point in the stable point set is transferred to the unstable point set, and the process returns to the step of calculating the stable value of the stable point set according to the deformation observation period, the stable point set and the unstable point set.

[0048] Optionally, the stability value of the stable point set is calculated based on the deformation observation period, the stable point set, and the unstable point set, including:

[0049] By formula:

[0050] ;

[0051] ;

[0052] ;

[0053] ;

[0054] Calculate the stable value of the stable point set ;

[0055] in, represents the empirical variance between the residuals of two observation periods in the deformation observation period, represents a stable point set statistic, The weight matrix representing the difference between the updated stable point observations, represents the observation degrees of freedom of the stable point set, The weight matrix representing the difference between the observation results of two observation periods, The weight matrix representing the difference between stable point observations, represents the inverse of the cofactor matrix between stable and unstable points, and Transpose each other, The weight matrix representing the difference between observations at non-stationary points, The matrix representing the updated gap value of the unstable point, The matrix representing the difference between the observations of two observation periods, The matrix representing the gap values ​​of the stable points, The matrix representing the gap values ​​of unstable points, Represents a splicing operation.

[0056] Optionally, the temporal deformation results of the target area in the target time period under the multi-reference point constraint are obtained based on all reference points, including:

[0057] All reference points are thinned out to obtain multiple final reference points;

[0058] Construct an InSAR deformation result solution model based on all final reference points;

[0059] The InSAR deformation result calculation model is solved to obtain the temporal deformation results of the target area in the target time period.

[0060] Optionally, the InSAR deformation result solution model is:

[0061] ;

[0062] in, Represents the set of arc observation values ​​between all two adjacent pixel points, represents the set of observations of all final reference points and all GNSS monitoring stations, represents the arc coefficient matrix, represents the coefficient matrix, , represents the coefficient matrix of all final reference points, represents the coefficient matrix of all GNSS monitoring stations, Indicates the result of solving the phase of the reference point. represents the phase results of all final reference points and all GNSS monitoring stations, represents the model observation residual.

[0063] In a second aspect, an embodiment of the present application provides a multi-reference point deformation solution device for time-series InSAR, comprising:

[0064] The phase unwrapping module is used to acquire multiple SAR images of the target area within the target time period, generate differential interferometry pairs based on all the SAR images, and perform phase unwrapping on all the differential interferometry pairs to obtain amplitude, coherence, differential phase, and unwrapped phase data. The earliest SAR image acquired among all the SAR images is the reference SAR image.

[0065] The calculation module is used to calculate the values ​​of multiple data indicators for each pixel in the reference SAR image using amplitude, coherence, differential phase and unwrapped phase data, and calculate the comprehensive score of each pixel based on the values ​​of all data indicators of each pixel; the data indicators are used to evaluate the quality of the deformation data corresponding to the pixel;

[0066] A determination module is used to sort all comprehensive scores from large to small, select the pixels corresponding to the top multiple comprehensive scores in the sorting results as candidate reference points, and determine multiple reference points from all candidate reference points; the stability of the area corresponding to the reference point in the target area is greater than the stability of the areas corresponding to all other candidate reference points;

[0067] The obtaining module is used to obtain the temporal deformation results of the target area in the target time period under the multi-reference point constraints based on all reference points.

[0068] In a third aspect, an embodiment of the present application provides a terminal device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the multi-reference point deformation solution method for time-series InSAR is implemented.

[0069] In a fourth aspect, an embodiment of the present application provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the above-mentioned multi-reference point deformation solution method for time-series InSAR.

[0070] The above solution of the present application has the following beneficial effects:

[0071] In some embodiments of the present application, multiple SAR images of a target area in a target time period are acquired, and differential interferometry pairs are generated based on all the SAR images. Phase unwrapping is performed on all the differential interferometry pairs to obtain amplitude, coherence, differential phase, and unwrapped phase data. The amplitude, coherence, differential phase, and unwrapped phase data are then used to calculate the values ​​of multiple data indicators for each pixel in the reference SAR image. A comprehensive score for each pixel is calculated based on the values ​​of all the data indicators for each pixel. All the comprehensive scores are then sorted from large to small, and the pixels corresponding to the top multiple comprehensive scores in the sorted results are used as candidate reference points. Multiple reference points are determined from all the candidate reference points, and then the temporal deformation results of the target area in the target time period under the constraints of the multiple reference points are obtained based on all the reference points. Among them, candidate reference points are determined based on the comprehensive score, and then the reference points are determined from the candidate reference points. The data quality of the pixel points and the stability of the corresponding area are considered in turn, which effectively improves the effectiveness and accuracy of the reference points. The time series deformation results are obtained based on the precise reference points, which can improve the accuracy of the time series deformation results, thereby effectively improving the accuracy of the time series InSAR deformation solution.

[0072] Other beneficial effects of the present application will be described in detail in the subsequent specific implementation section. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the following briefly introduces the drawings required for use in the embodiments or descriptions of the prior art. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without any creative work.

[0074] Figure 1A flowchart of a multi-reference point deformation solution method for time-series InSAR provided in one embodiment of the present application;

[0075] Figure 2 A schematic diagram of relative error provided in an embodiment of the present application;

[0076] Figure 3 A schematic diagram of the structure of a multi-reference point deformation solution device for time-series InSAR provided in one embodiment of the present application;

[0077] Figure 4 A schematic diagram of the structure of a terminal device provided in one embodiment of the present application. DETAILED DESCRIPTION

[0078] In the following description, specific details such as specific system structures and techniques are provided for purposes of illustration rather than limitation to facilitate a thorough understanding of the embodiments of the present application. However, it will be apparent to those skilled in the art that the present application may be implemented in other embodiments without these specific details. In other cases, detailed descriptions of well-known systems, devices, circuits, and methods are omitted to avoid obscuring the description of the present application with unnecessary detail.

[0079] It should be understood that when used in the present specification and the appended claims, the term "comprising" indicates the presence of described features, integers, steps, operations, elements and / or components, but does not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components and / or collections thereof.

[0080] It will also be understood that the term "and / or" used in this specification and the appended claims refers to and includes any and all possible combinations of one or more of the associated listed items.

[0081] As used in this specification and the appended claims, the term "if" can be interpreted as "when" or "upon" or "in response to determining" or "in response to detecting," depending on the context. Similarly, the phrase "if it is determined" or "if [described condition or event] is detected" can be interpreted as meaning "upon determination" or "in response to determining" or "upon detection of [described condition or event]" or "in response to detecting [described condition or event]," depending on the context.

[0082] In addition, in the description of the present application specification and the appended claims, the terms "first", "second", "third", etc. are only used to distinguish the descriptions and cannot be understood as indicating or implying relative importance.

[0083] References to "one embodiment" or "some embodiments" in this specification mean that a particular feature, structure, or characteristic described in conjunction with that embodiment is included in one or more embodiments of the present application. Thus, phrases such as "in one embodiment," "in some embodiments," "in other embodiments," and "in other embodiments" appearing in various places in this specification do not necessarily refer to the same embodiment, but rather mean "one or more but not all embodiments," unless otherwise specifically emphasized. The terms "including," "comprising," "having," and variations thereof all mean "including but not limited to," unless otherwise specifically emphasized.

[0084] To address the problem of low accuracy of existing time-series InSAR deformation solution, an embodiment of the present application provides a multi-reference point deformation solution method for time-series InSAR. The multi-reference point deformation solution method obtains multiple SAR images of a target area in a target time period, generates differential interferometer pairs based on all SAR images, performs phase unwrapping on all differential interferometer pairs to obtain amplitude, coherence, differential phase and unwrapped phase data, and then uses the amplitude, coherence, differential phase and unwrapped phase data to calculate the values ​​of multiple data indicators for each pixel in the reference SAR image, and calculates a comprehensive score for each pixel based on the values ​​of all data indicators of each pixel. All comprehensive scores are then sorted from large to small, and the pixels corresponding to the top multiple comprehensive scores in the sorted results are used as candidate reference points. Multiple reference points are determined from all candidate reference points, and then the time-series deformation results of the target area in the target time period under the multi-reference point constraint are obtained based on all reference points. Among them, candidate reference points are determined based on the comprehensive score, and then the reference points are determined from the candidate reference points. The data quality of the pixel points and the stability of the corresponding area are considered in turn, which effectively improves the effectiveness and accuracy of the reference points. The time series deformation results are obtained based on the precise reference points, which can improve the accuracy of the time series deformation results, thereby effectively improving the accuracy of the time series InSAR deformation solution.

[0085] Next, the multi-reference point deformation solution method of the time-series InSAR provided in this application is exemplified.

[0086] like Figure 1 As shown, the multi-reference point deformation solution method for time-series InSAR provided in this application includes the following steps:

[0087] Step 11: Acquire multiple SAR images of the target area in the target time period, generate differential interferometry pairs based on all SAR images, perform phase unwrapping on all differential interferometry pairs, and obtain amplitude, coherence, differential phase, and unwrapped phase data.

[0088] The earliest SAR image among all SAR images is the reference SAR image. The target area is the area where surface deformation analysis is required, such as a mountainous area. The target time period is the time period when the terrain deformation of the target area needs to be analyzed, such as before and after a volcanic eruption or during the period of mineral mining. SAR images are obtained by satellites (such as Sentinel-1) and can be obtained by accessing websites or platforms that provide public SAR images. There is a time interval when the satellites take pictures. For multiple SAR images of the target area in the target time period, the acquisition time of different SAR images is different. For example, if the target time period is from March 2 to April 9, the three SAR images acquired in this time period are March 3, March 15, and April 1, respectively, with March 3 being the earliest acquisition time.

[0089] It should be noted that the above-mentioned steps for phase unwrapping all SAR images to obtain unwrapped phase data include preprocessing the SAR images and related data (such as an external reference digital elevation model (DEM)). This includes decompression of the original SAR images, generation of a single-look composite (SLC) file, DEM correction, registration, cropping, and generation of multi-look intensity maps. Next, differential interferometric synthetic aperture radar (D-InSAR) technology is used to obtain a multi-look differential interferogram set. A small baseline set is constructed by setting appropriate spatiotemporal baseline thresholds, and adaptive filtering is applied to the differential interferograms. This results in differential interferograms and phase coherence results corresponding to each interferogram period. Then, a high coherence point is preliminarily selected as the unwrapping reference point for phase unwrapping to obtain the initial unwrapped phase map. At this time, trend term error removal and atmospheric noise reduction are not performed. The singular value decomposition method is used to solve the linear deformation rate of the unwrapped phase map. The cumulative linear deformation at each time point can be obtained by integration in the time domain, and the result is converted from radians to meters (m) to obtain the final unwrapped interferometer data.

[0090] Step 12: Calculate the values ​​of multiple data indicators for each pixel in the reference SAR image using the amplitude, coherence, differential phase, and unwrapped phase data, and calculate a comprehensive score for each pixel based on the values ​​of all data indicators for each pixel.

[0091] The aforementioned metrics are used to evaluate the quality of the deformation data corresponding to a pixel. These metrics include amplitude dispersion, average coherence, phase derivative deviation, residual point score, phase standard deviation, and stacking average deformation rate. The above comprehensive score is used to evaluate the data quality of a pixel by integrating the values ​​of all its metrics.

[0092] It should be noted that the values ​​of various pixel indicators can be calculated using the amplitude, coherence, differential phase, and unwrapped phase data using the aforementioned data indicator calculation formulas. For example, for amplitude dispersion, the amplitude dispersion of each pixel is calculated based on the reflection intensity of all pixels in the amplitude data using the formula for calculating amplitude dispersion. For average coherence, the coherence of all pixels in the coherence data is averaged to obtain the average coherence. For residual point scores, the residual point score of each pixel is calculated using the residual point score calculation formula using the differential phase data. For phase standard deviation, the phase standard deviation of each pixel is calculated using the unwrapped phase data using the phase standard deviation calculation formula. For stacking average deformation rate, the stacking average deformation rate indicator of each pixel is calculated using the unwrapped phase data using the stacking average deformation rate indicator calculation formula. Each type of indicator has different properties and standards, which is not conducive to merging into a comprehensive evaluation indicator. They all need to be converted into positive indicators, that is, larger indicator values ​​indicate better data quality. At the same time, each type of indicator needs to be standardized to remove dimension effects. The forwarding process is:

[0093] For very small indicators, the forward expression is:

[0094] ;

[0095] in, Indicates the value of the data indicator after positive transformation, Indicates the value of the first data indicator, Indicates the The value of the data indicator, Indicates the The value of a data indicator.

[0096] For intermediate indicators, the forward expression is:

[0097] ;

[0098] The best value is , .

[0099] For interval indicators, the forward expression is:

[0100] ;

[0101] The best interval is , .

[0102] Normalize the indicators that have been positive. You can use Min-Max normalization, the formula is as follows:

[0103] ;

[0104] in, For the The minimum value of a data indicator, For the The maximum value of a data indicator.

[0105] The steps for calculating the comprehensive score of each pixel based on the values ​​of all data indicators of each pixel are as follows:

[0106] By formula:

[0107] ;

[0108] Calculate the Comprehensive score of pixels ;

[0109] in, Indicates the The scoring coordinate vector corresponding to the pixel point, , Indicates the number of pixels, represents the optimal solution vector, represents the worst solution vector, Indicates the The distance between the pixel point and the worst solution vector, Indicates the The distance between a pixel point and the optimal solution vector:

[0110] ;

[0111] ;

[0112] in, Indicates the The optimal value of a data indicator, Indicates the The worst value of a data indicator, Represents the weight of each indicator, which can be obtained based on the hierarchical analysis method or entropy weight method. Indicates the The first pixel The value of the data indicator, Indicates the number of data indicators.

[0113] In step 13, all comprehensive scores are sorted from large to small, the pixel points corresponding to the top comprehensive scores in the sorting result are selected as candidate reference points, and multiple reference points are determined from all candidate reference points.

[0114] The stability of the region corresponding to the reference point in the target region is greater than that of the regions corresponding to all other candidate reference points. Stability reflects the degree of regional deformation. The terrain of the region corresponding to the reference point is relatively stable, making it suitable as a reference region for deformation of the entire region. This serves as a benchmark for deformation analysis, ensuring that the deformation results of the entire target region are relative to the reference region, thereby avoiding errors caused by shifts in the entire target region or external factors. The number of candidate reference points is set according to actual needs, typically the first 20% of the ranked results.

[0115] In some embodiments of the present application, the step of determining multiple reference points from all candidate reference points includes:

[0116] The first step is to divide the target time period into multiple observation groups and calculate the statistics of each observation group.

[0117] The above observation groups include multiple observation periods. An observation period is obtained by dividing the target time period into multiple time periods. Usually, there is at least one SAR image in each observation period. Then, at least one observation period is divided into an observation group in chronological order.

[0118] Specifically, through the formula:

[0119] ;

[0120] ;

[0121] ;

[0122] Calculating statistics .

[0123] in, represents the sum of the degrees of freedom of all observation periods in the observation group, represents the combined unit weight variance of all observation periods, Indicates the The unit weight variance of the observation data in the observation period is Indicates the degrees of freedom for each observation period, represents the number of observation periods, Expressed as a residual correction vector, represents the observation weight matrix, Represents the residual quadratic form, which is used to simplify the expression.

[0124] It should be noted that the statistics It conforms to the chi-square distribution. The critical value of the statistic is determined according to the significance level and degrees of freedom of the kaffir distribution. The specific value can be determined by referring to the chi-square distribution table.

[0125] The second step is to determine whether there is any observation group whose statistic is less than or equal to the critical value of the statistic among all observation groups; the observation group includes multiple observation periods.

[0126] If so, it is considered that the accuracy consistency test of all observation groups has passed, and the deformation observation period in which the terrain deformation occurs is determined from all observation periods in the observation group that meet the statistic less than or equal to the statistic critical value, and multiple reference points are determined from all candidate reference points based on the deformation observation period.

[0127] Otherwise, divide each observation group into multiple sub-observation groups, and take all sub-observation groups as multiple observation groups, and return to the step of calculating the statistics for each observation group.

[0128] It should be noted that the above step of determining the deformation observation period in which the terrain deformation occurs from all observation groups that satisfy the statistic less than or equal to the statistic critical value includes: for each two adjacent observation periods in the observation group, calculating the deformation gap statistic between the two adjacent observation periods, if the deformation gap statistic is greater than the deformation gap statistic threshold (usually the 5% significance statistic threshold, i.e. ,in If is the statistical threshold, obtained by looking up the F-distribution statistic table), then the two adjacent observation periods are considered deformation observation periods. If the deformation gap statistic is less than or equal to the deformation gap statistic threshold, then the two adjacent observation periods are not processed and it is assumed that no deformation occurred between the two periods.

[0129] Specifically, through the formula:

[0130] ;

[0131] ;

[0132] ;

[0133] Calculate the deformation gap statistic between two adjacent observation periods ;

[0134] in, represents the empirical variance between the residuals of two adjacent observation periods, represents the sum of the degrees of freedom of two adjacent observation periods, represents the coordinate difference between two adjacent observation periods, The matrix representing the difference between the observation results of two observation periods, where the elements are the deformation of the area corresponding to the pixel point between the two observation periods, called the pixel gap value. The weight matrix representing the difference between the observation results of two observation periods is weighted according to the phase coherence. represents the number of independent observation conditions, Represents the observation value correction matrix after adjustment, where the elements are the observation value correction numbers for each pixel. Represents the observation weight matrix, where the elements are the observation weights for each observation. represents the residual sum of squares of the first period observations, represents the residual sum of squares of the second period observations, Represents a transpose operation.

[0135] The step of determining multiple reference points from all candidate reference points according to the deformation observation period includes:

[0136] First, all candidate reference points are divided into a stable point set and an unstable point set.

[0137] For example, the stable point set and the unstable point set may be randomly divided, for example, all candidate reference points may be divided into the stable point set and the unstable point set in a ratio of 8:2.

[0138] Then, the stability value of the stable point set is calculated based on the deformation observation period, the stable point set and the unstable point set.

[0139] Specifically, through the formula:

[0140] ;

[0141] ;

[0142] ;

[0143] ;

[0144] Calculate the stable value of the stable point set .

[0145] in, represents the empirical variance between the residuals of two observation periods in the deformation observation period, represents a stable point set statistic, The weight matrix representing the difference between the updated stable point observations, represents the observation degrees of freedom of the stable point set, The weight matrix representing the difference between the observation results of two observation periods, The weight matrix representing the difference between stable point observations, represents the inverse of the cofactor matrix between stable and unstable points, and Transpose each other, The weight matrix representing the difference between observations at non-stationary points, The matrix representing the updated gap value of the unstable point, The matrix representing the difference between the observations of two observation periods (i.e. the gap values ​​of all points), The matrix representing the gap values ​​of stable points, including the gap values ​​of all stable points, The matrix representing the gap values ​​of unstable points, including the gap values ​​of all unstable points, Represents a splicing operation.

[0146] Finally, determine whether the stable value is less than the stable threshold (the stable threshold is and Comparison is done by looking up the table. is the F distribution value with significance α, and Table 2 is the F distribution statistics table).

[0147] If so, all candidate reference points in the stable point set are used as reference points.

[0148] Otherwise, a candidate point in the stable point set is transferred to the unstable point set, and the process returns to the step of calculating the stable value of the stable point set according to the deformation observation period, the stable point set and the unstable point set.

[0149] For example, when transferring a candidate point from the stable point set to the unstable point set, starting from the first candidate reference point, the candidate reference points are transferred to the unstable point set in sequence, and it is determined whether the candidate reference points satisfy the following conditions after transfer:

[0150] ;

[0151] in, Indicates the The gap value of the unstable point set after removing it from the stable point set, Indicates the The weight matrix of the unstable point set after removing it from the stable point set, Indicates the The gap value of the unstable point set after removing it from the stable point set.

[0152] If satisfied, it is determined to transfer the candidate reference point to the unstable point set, and return to the step of calculating the stable value of the stable point set according to the deformation observation period, the stable point set and the unstable point set.

[0153] If not, the candidate reference point is returned to the stable point set, and the next candidate reference point is transferred to the unstable point set to determine whether it satisfies the above formula, and so on, until all points are traversed.

[0154] It should be noted that when there is more than one observation group that satisfies the statistic less than or equal to the critical value of the statistic, resulting in more than one pair of deformation observation periods, or when there is more than one pair of observation periods that satisfies the deformation gap statistic greater than the deformation gap statistic threshold, the reference point set of each pair of deformation observation periods is obtained through the above process, and all reference points in the intersection of all reference point sets are used as reference points for participating in subsequent steps.

[0155] Step 14: Obtain the temporal deformation result of the target area in the target time period under the multi-reference point constraint based on all reference points.

[0156] In some embodiments of the present application, the step of obtaining the temporal deformation result of the target area in the target time period under the multi-reference point constraint based on all reference points includes:

[0157] In the first step, all reference points are thinned out to obtain multiple final reference points.

[0158] For example, in order to ensure that there is at most one area corresponding to a reference point within a preset range in the target area, uniform spatial thinning is performed. For the case where there are multiple areas corresponding to reference points within the preset range, the reference points are randomly removed, and only one reference point is retained and the retained reference point is used as the final reference point.

[0159] In the second step, the InSAR deformation result solution model is constructed based on all the final reference points.

[0160] Specifically, the InSAR deformation result solution model is:

[0161] ;

[0162] in, Represents the set of arc observation values ​​between all two adjacent pixel points, represents the set of observations of all final reference points and all GNSS monitoring stations, represents the arc coefficient matrix, represents the coefficient matrix, , represents the coefficient matrix of all final reference points, represents the coefficient matrix of all GNSS monitoring stations, Indicates the phase result of the reference point, which can be obtained by the coherence threshold. If the point satisfies the coherence of 0.3 or above, it can be regarded as the reference point. represents the phase results of all final reference points and all GNSS monitoring stations, represents the model observation residual.

[0163] The third step is to solve the InSAR deformation result calculation model to obtain the temporal deformation results of the target area in the target time period.

[0164] The time series deformation result is obtained by solving The value of . The value of can be directly replaced by the reference point value.

[0165] For example, we can make , , The above formula can be simplified to , the linear equations constructed at this time are large in scale and the elements are relatively sparse. The simplified formula is solved by the Least Squares Minimal Residual Method (LSMR). The problem is projected into a low-dimensional Krylov subspace through Lanczos bidiagonalization, and the residual is gradually minimized. , avoiding direct matrix decomposition, thereby significantly reducing memory consumption and computational complexity, while making the solution more stable, let:

[0166] ;

[0167] In the formula and is a Lanczos bidiagonalized orthogonal matrix, is a bidiagonal matrix, , is the next step orthogonal basis vector The module length, , and are the standard basis vectors, is the optimal solution in the subspace, The optimal solution in the target area is obtained. Finally, the time series deformation results of the target area in the target time period are obtained through block-by-block parallel solving. Specifically, because the program solving blocks are independent, the Python multiprocessing library is used here to decompose the solving process into single module tasks, making full use of the computer CPU resources to synchronously solve multiple single module results.

[0168] It is worth mentioning that candidate reference points are determined based on the comprehensive score, and then reference points are determined from the candidate reference points. The data quality of the pixel points and the stability of the corresponding area are considered in turn, which effectively improves the effectiveness and accuracy of the reference points. Obtaining the time series deformation results based on precise reference points can improve the accuracy of the time series deformation results, thereby effectively improving the accuracy of the time series InSAR deformation solution.

[0169] The method of the present application is illustrated below with reference to a specific example.

[0170] For this study, we selected an orbiting Sentinel-1 SAR dataset covering an island from May 2018 to March 2019. The dataset had an ascending angle of incidence of approximately 40.57 degrees and an azimuth of approximately -12.28 degrees (counterclockwise). TanDEM-X, a digital elevation model acquired by Earth observation satellites with a spatial resolution of 90 meters, was used as the external reference DEM.

[0171] Using the GNSS data of a measurement laboratory, 1 / 3 of the sites were selected to verify the accuracy of the method provided in this application, and the rest were used as multi-reference point constraint benchmarks. The time resolution of the GNSS data was daily, and the monitoring sites covering the study period were selected.

[0172] To demonstrate the advantages of the method provided in this application, the accuracy of the traditional single reference point unwrapping and the Small Baseline Subset Interferometric Synthetic Aperture Radar (SBAS-InSAR) processing results using different sizes of spatiotemporal window filtering to remove atmospheric effects are compared with the method. The relative error distribution of each method station is shown in the figure below. Figure 2As shown in the figure, the horizontal axis represents the method number. Method 1 is the method of this application, and the other methods are SBAS-InSAR methods that use different spatiotemporal filter windows to remove atmospheric effects. The vertical axis represents the relative error of terrain deformation. The relative error of a single GNSS site is the relative error of a single GNSS site, the average of the traditional method is the average relative error of deformation using the traditional method, and the average of the present method is the average relative error of deformation using the present method. The vertical line represents the standard deviation of the solution results. The root mean square error (RMSE) is calculated for each GNSS site data and the results obtained by different methods. To ensure comparability between site data, the calculated RMSE is divided by the maximum deformation value of the GNSS site, eliminating the dimension and obtaining the relative error. Because most of the validation sites are stable sites with observations that are essentially stable near zero, the RMSE divided by the maximum GNSS site observation may result in a relative error greater than 1. While this may result in a large mean relative error, the actual error is likely smaller in magnitude. For this reason, the relative relationship between the methods is sufficient.

[0173] from Figure 2 As can be seen from the figure, the method provided by the present application can greatly improve the accuracy of deformation solution compared with traditional methods. It can be seen from the figure that the relative error mean of the method provided by the present application is the lowest, and its standard deviation is also lower. Compared with the traditional method, the average relative error is reduced by about 35%, indicating that it is more consistent with the data observation of the GNSS verification site. Other traditional SBAS methods have a certain impact on the results when different windows are used, but their errors fluctuate. For different areas, the larger the filter window, the better. In order to obtain the best results, it may be more dependent on the expert experience of the data processor, which is more subjective. Therefore, the method provided by the present application can be used to solve deformation in areas that are greatly affected by the atmosphere or other noise.

[0174] The advantages of the method proposed in this application are: establishing a semi-quantitative selection model for MT-InSAR initial multiple reference points with multi-factor fusion, promoting the initial screening of reference points from the "manual experience" mode to the "semi-quantitative automatic" mode; proposing an automatic selection method for MT-InSAR multiple reference points under the background of non-equal precision observations, realizing the reference point selection from "experience-given" to "automatic identification"; developing MT-InSAR time series deformation monitoring technology under the constraints of multiple reference points, realizing the MT-InSAR solution from "single reference point" to "multiple reference points".

[0175] The following is an exemplary description of the multi-reference point deformation solution device for time-series InSAR provided in this application.

[0176] like Figure 3As shown, an embodiment of the present application provides a multi-reference point deformation solution device for time-series InSAR. The multi-reference point deformation solution device 300 for time-series InSAR includes:

[0177] Phase unwrapping module 301 is used to acquire multiple SAR images of the target area within a target time period, generate differential interferometry pairs based on all the SAR images, and perform phase unwrapping on all the differential interferometry pairs to obtain amplitude, coherence, differential phase, and unwrapped phase data; the earliest acquired SAR image among all the SAR images is used as the reference SAR image;

[0178] Calculation module 302 is used to calculate the values ​​of multiple data indicators for each pixel in the reference SAR image using amplitude, coherence, differential phase, and unwrapped phase data, and to calculate a comprehensive score for each pixel based on the values ​​of all data indicators for each pixel; the data indicators are used to evaluate the quality of the deformation data corresponding to the pixel;

[0179] Determination module 303 is used to sort all comprehensive scores from largest to smallest, select the pixels corresponding to the top multiple comprehensive scores in the sorting result as candidate reference points, and determine multiple reference points from all candidate reference points; the stability of the area corresponding to the reference point in the target area is greater than the stability of the areas corresponding to all other candidate reference points;

[0180] The obtaining module 304 is configured to obtain, based on all reference points, a temporal deformation result of the target area in a target time period under the multi-reference point constraint.

[0181] It should be noted that the information interaction, execution process, etc. between the above-mentioned devices / units are based on the same concept as the method embodiment of this application. Their specific functions and technical effects can be found in the method embodiment section and will not be repeated here.

[0182] Those skilled in the art can clearly understand that, for the convenience and brevity of description, only the division of the above-mentioned functional units and modules is used as an example for illustration. In actual applications, the above-mentioned functions can be distributed and completed by different functional units and modules as needed, that is, the internal structure of the device can be divided into different functional units or modules to complete all or part of the functions described above. The functional units and modules in the embodiment can be integrated into one processing unit, or each unit can exist physically alone, or two or more units can be integrated into one unit. The above-mentioned integrated unit can be implemented in the form of hardware or in the form of software functional units. In addition, the specific names of the functional units and modules are only for the convenience of distinguishing each other, and are not used to limit the scope of protection of this application. The specific working process of the units and modules in the above-mentioned system can refer to the corresponding process in the aforementioned method embodiment, and will not be repeated here.

[0183] like Figure 4 As shown, an embodiment of the present application provides a terminal device. The terminal device D10 of this embodiment includes: at least one processor D100 ( Figure 4 Only one processor is shown in the figure), a memory D101, and a computer program D102 stored in the memory D101 and executable on the at least one processor D100, wherein the processor D100 implements the steps of any of the above-mentioned method embodiments when executing the computer program D102.

[0184] Specifically, when the processor D100 executes the computer program D102, it obtains multiple SAR images of the target area in the target time period, generates differential interference pairs based on all the SAR images, performs phase unwrapping on all the differential interference pairs to obtain amplitude, coherence, differential phase and unwrapped phase data, and then uses the amplitude, coherence, differential phase and unwrapped phase data to calculate the values ​​of multiple data indicators of each pixel point in the reference SAR image, and calculates the comprehensive score of each pixel point based on the values ​​of all data indicators of each pixel point, and then sorts all the comprehensive scores from large to small, and uses the pixel points corresponding to the first multiple comprehensive scores in the sorting results as candidate reference points, and determines multiple reference points from all the candidate reference points, and then obtains the temporal deformation results of the target area in the target time period under the constraints of multiple reference points based on all the reference points. Among them, candidate reference points are determined based on the comprehensive score, and then the reference points are determined from the candidate reference points. The data quality of the pixel points and the stability of the corresponding area are considered in turn, which effectively improves the effectiveness and accuracy of the reference points. The time series deformation results are obtained based on the precise reference points, which can improve the accuracy of the time series deformation results, thereby effectively improving the accuracy of the time series InSAR deformation solution.

[0185] The processor D100 may be a central processing unit (CPU), or may be another general-purpose processor, a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components. A general-purpose processor may be a microprocessor or any conventional processor.

[0186] In some embodiments, the memory D101 may be an internal storage unit of the terminal device D10, such as a hard disk or memory of the terminal device D10. In other embodiments, the memory D101 may also be an external storage device of the terminal device D10, such as a plug-in hard disk, a smart memory card (SMC), a secure digital (SD) card, a flash card, etc. equipped on the terminal device D10. Furthermore, the memory D101 may include both an internal storage unit of the terminal device D10 and an external storage device. The memory D101 is used to store an operating system, application programs, a boot loader, data, and other programs, such as the program code of the computer program. The memory D101 may also be used to temporarily store data that has been output or is about to be output.

[0187] An embodiment of the present application further provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the steps in the above-mentioned various method embodiments can be implemented.

[0188] An embodiment of the present application provides a computer program product. When the computer program product is run on a terminal device, the terminal device can implement the steps in the above-mentioned method embodiments when executing the computer program product.

[0189] If the integrated unit is implemented as a software functional unit and sold or used as a standalone product, it can be stored in a computer-readable storage medium. Based on this understanding, the present application implements all or part of the process steps in the above-mentioned method embodiments by instructing the relevant hardware through a computer program. The computer program can be stored in a computer-readable storage medium. When executed by a processor, the computer program can implement the steps of each of the above-mentioned method embodiments. The computer program includes computer program code, which can be in source code form, object code form, executable file, or some intermediate form. The computer-readable medium can include at least: any entity or device capable of carrying the computer program code to the multi-reference point deformation solution method apparatus / terminal device of the time-series InSAR, recording medium, computer memory, read-only memory (ROM), random access memory (RAM), electrical carrier signal, telecommunication signal, and software distribution medium. Examples include a USB flash drive, mobile hard drive, magnetic disk, or optical disk.

[0190] In the above embodiments, the description of each embodiment has its own focus. For parts that are not described or recorded in detail in a certain embodiment, reference can be made to the relevant description of other embodiments.

[0191] Those skilled in the art will appreciate that the units and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware or a combination of computer software and electronic hardware. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professional and technical personnel can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0192] The above is a preferred embodiment of the present application. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles described in the present application. These improvements and modifications should also be regarded as the scope of protection of the present application.

Claims

1. A multi-reference point deformation solution method for time series InSAR, characterized by: include: Acquire multiple SAR images of the target area within the target time period, generate differential interferometry pairs based on all the SAR images, perform phase unwrapping on all the differential interferometry pairs, and obtain amplitude, coherence, differential phase, and unwrapped phase data. The earliest SAR image acquired among all the SAR images is the reference SAR image. Calculating the values ​​of multiple data indicators for each pixel in the reference SAR image using the amplitude, coherence, differential phase, and unwrapped phase data, and calculating a comprehensive score for each pixel based on the values ​​of all the data indicators for each pixel; the data indicators are used to evaluate the quality of the deformation data corresponding to the pixel; Sorting all comprehensive scores from largest to smallest, selecting pixel points corresponding to the top multiple comprehensive scores in the sorting result as candidate reference points, and determining multiple reference points from all candidate reference points; the stability of the area corresponding to the reference point in the target area is greater than the stability of the areas corresponding to all other candidate reference points; The temporal deformation result of the target area in the target time period under the multi-reference point constraint is obtained based on all reference points.

2. The multi-reference point deformation solution method according to claim 1, characterized in that: Calculating the comprehensive score of each pixel point based on the values ​​of all data indicators of each pixel point includes: By formula: ； Calculate the Comprehensive score of pixels ; in, Indicates the The scoring coordinate vector corresponding to the pixel point, , Indicates the number of pixels, represents the optimal solution vector, represents the worst solution vector, Indicates the The distance between the pixel point and the worst solution vector, Indicates the The distance between a pixel point and the optimal solution vector: ； ； in, Indicates the The optimal value of a data indicator, Indicates the The worst value of a data indicator, Indicates the weight of each indicator, Indicates the The first pixel The value of the data indicator, Indicates the number of data indicators.

3. The multi-reference point deformation solution method according to claim 1, characterized in that: The determining of multiple reference points from all candidate reference points includes: Dividing the target time period into a plurality of observation groups; Calculate statistics for each observation group; Determine whether there is an observation group among all observation groups whose statistic is less than or equal to a critical statistic value; the observation group includes multiple observation periods; If so, determining a deformation observation period in which terrain deformation occurs from all observation periods in the observation group satisfying the statistic being less than or equal to the statistic critical value, and determining multiple reference points from all candidate reference points according to the deformation observation period; Otherwise, each observation group is divided into a plurality of sub-observation groups, and all the sub-observation groups are used as the plurality of observation groups, and the process returns to the step of calculating the statistics of each observation group.

4. The multi-reference point deformation solution method according to claim 3, characterized in that: The calculation of statistics for each observation group includes: By formula: ； ； ； Calculating statistics ; in, represents the sum of the degrees of freedom of all observation periods in the observation group, represents the combined unit weight variance of all observation periods, Indicates the The unit weight variance of the observation data in the observation period is Indicates the degrees of freedom for each observation period, represents the number of observation periods, Expressed as a residual correction vector, represents the observation weight matrix, represents the residual quadratic form.

5. The multi-reference point deformation solution method according to claim 3, characterized in that: Determining the deformation observation period in which the terrain deformation occurs from all observation periods in the observation group satisfying the statistic being less than or equal to the statistic critical value includes: For each two adjacent observation periods in the observation group, a deformation gap statistic between the two adjacent observation periods is calculated. If the deformation gap statistic is greater than a deformation gap statistic threshold, the two adjacent observation periods are considered to be deformation observation periods.

6. The multi-reference point deformation solution method according to claim 5, characterized in that: The calculating of the deformation gap statistic between the two adjacent observation periods includes: By formula: ； ； ； Calculate the deformation gap statistic between two adjacent observation periods ; in, represents the empirical variance between the residuals of two adjacent observation periods, represents the sum of the degrees of freedom of two adjacent observation periods, represents the coordinate difference between two adjacent observation periods, The matrix representing the difference between the observations of two observation periods, The weight matrix representing the difference between the observation results of two observation periods, represents the number of independent observation conditions, Represents the matrix of observation corrections after adjustment, represents the observation weight matrix, represents the residual sum of squares of the first period observations, represents the residual sum of squares of the second period observations, Represents a transpose operation.

7. The multi-reference point deformation solution method according to claim 6, characterized in that: The step of determining a plurality of reference points from all candidate reference points according to the deformation observation period includes: Divide all candidate reference points into stable point sets and unstable point sets; Calculating a stable value of the stable point set according to the deformation observation period, the stable point set and the unstable point set; Determining whether the stability value is less than a stability threshold; If so, all candidate reference points in the stable point set are used as reference points; Otherwise, a candidate point in the stable point set is transferred to the unstable point set, and the process returns to the step of calculating the stable value of the stable point set according to the deformation observation period, the stable point set and the unstable point set.

8. The multi-reference point deformation solution method according to claim 7, characterized in that: The calculating of the stable value of the stable point set according to the deformation observation period, the stable point set and the unstable point set includes: By formula: ； ； ； ； Calculate the stable value of the stable point set ; in, represents the empirical variance between the residuals of two observation periods in the deformation observation period, represents a stable point set statistic, The weight matrix representing the difference between the updated stable point observations, represents the observation degrees of freedom of the stable point set, The weight matrix representing the difference between the observation results of two observation periods, The weight matrix representing the difference between stable point observations, represents the inverse of the cofactor matrix between stable and unstable points, and Transpose each other, The weight matrix representing the difference between observations at non-stationary points, The matrix representing the updated gap value of the unstable point, The matrix representing the difference between the observations of two observation periods, The matrix representing the gap values ​​of the stable points, The matrix representing the gap values ​​of unstable points, Represents a splicing operation.

9. The multi-reference point deformation solution method according to claim 1, characterized in that: The step of obtaining the temporal deformation result of the target area in the target time period under the multi-reference point constraint based on all reference points includes: All reference points are thinned out to obtain multiple final reference points; Construct an InSAR deformation result solution model based on all final reference points; The InSAR deformation result calculation model is solved to obtain the time series deformation result of the target area in the target time period.

10. The multi-reference point deformation solution method according to claim 9, characterized in that: The InSAR deformation result solution model is: ； in, Represents the set of arc observation values ​​between all two adjacent pixel points, represents the set of observations of all final reference points and all GNSS monitoring stations, represents the arc coefficient matrix, represents the coefficient matrix, , represents the coefficient matrix of all final reference points, represents the coefficient matrix of all GNSS monitoring stations, Indicates the result of solving the phase of the reference point. represents the phase results of all final reference points and all GNSS monitoring stations, represents the model observation residual.

Images