Settlement monitoring method based on small baseline subset deformation model and improved HK formula

Abstract

The application discloses a settlement monitoring method for solving the HK formula by improving a small baseline set deformation model, belongs to the technical field of settlement monitoring, and is used for settlement monitoring. The method comprises the following steps: pre-processing single look complex (SLC) images and digital elevation model (DEM) of a research area, and performing differential interference; constructing a small baseline set deformation model, determining ridge parameters, and performing ridge estimation and solving, so as to realize settlement monitoring. The application can maintain good solving precision when the observation data error is large, effectively monitor the time sequence settlement condition of an experimental area, and is smaller than an existing method whether RMSE or MAE is used for evaluation. Compared with the unimproved HK method, the application improves the RMSE by 12.52% and the MAE by 12.25%, which indicates that the application greatly improves the solving precision of the SBAS-InSAR deformation model ridge estimation and the precision of the SBAS-InSAR technology for monitoring the ground.

Description

Technical Field

[0001] This invention discloses a settlement monitoring method based on an improved HK formula for solving a small baseline set deformation model, belonging to the field of settlement monitoring technology. Background Technology

[0002] Existing technologies typically employ SBAS-InSAR technology to construct deformation models for subsidence prediction. In the normal matrix of the SBAS-InSAR deformation model, eigenvalues ​​gradually decrease; that is, the h-th eigenvalue is always smaller than the previous h-1 eigenvalues. However, the difference between the smallest eigenvalue at the end and the previous eigenvalue is much greater than the decreasing gradient corresponding to the other eigenvalues. Due to the unique characteristics of the deformation model eigenvalues, when using the conventional HK formula to determine the ridge parameter of the SBAS-InSAR deformation model, the result is often influenced by the last minimal eigenvalue, leading to a smaller ridge parameter. This affects the improvement of the model's ill-conditioning and thus also impacts the obtained ground subsidence results. Summary of the Invention

[0003] The purpose of this invention is to provide a settlement monitoring method for solving the HK formula using a small baseline set deformation model, in order to solve the problem in the prior art where the ridge parameters obtained by the conventional HK formula are small, which cannot effectively improve the ill-conditioned nature of the model and thus affect the settlement monitoring results.

[0004] An improved settlement monitoring method based on the small baseline set deformation model and the HK formula is proposed. This method includes preprocessing the single-view complex image (SLC) and digital elevation model (DEM) of the study area, performing differential interferometry, constructing a small baseline set deformation model, determining ridge parameters, and performing ridge estimation to achieve settlement monitoring.

[0005] Determining ridge parameters includes:

[0006] ;

[0007] In the formula, It is a ridge parameter estimate. and It is the intermediate value to be determined. Indicates to Take the average value.

[0008] Determining ridge parameters includes:

[0009] ;

[0010] ;

[0011] ;

[0012] ;

[0013] In the formula, These are settlement observation values. These are estimates from settlement observations. The number of observations. For the number of estimated values, It is a diagonal matrix composed of eigenvalues. yes eigenvalues, It is the eigenvector corresponding to the eigenvalue. It is the eigenvector matrix. yes The amount, It is the coefficient matrix of the small baseline set deformation model.

[0014] Constructing a small baseline set deformation model includes:

[0015] ;

[0016] In the formula, For the parameters to be determined, It's an error.

[0017] The settlement observation values ​​are:

[0018] ;

[0019] In the formula, It is differential interference phase. It is the radar wavelength.

[0020] The differential interference phase is:

[0021] ;

[0022] In the formula, yes Phase of time, yes The phase of a moment.

[0023] The coefficient matrix is ​​as follows:

[0024] ;

[0025] In the formula, It is the first of the coefficient matrix. line, indicating the first One interference pair, It is the main image time. It is the time for the secondary image.

[0026] The ridge estimation solution criterion is:

[0027] ;

[0028] In the formula, Ridge estimation solution criteria The estimated value, For ridge parameters.

[0029] Differentiating the ridge estimation criterion, we obtain the ridge estimation solution formula:

[0030] ;

[0031] ;

[0032] In the formula, For ridge estimation solution formula The estimated value, This is the weight matrix. As a unit array, and Same order.

[0033] Based on the primary and secondary images of the differential interferometric phase, a coefficient matrix is ​​constructed, ridge parameters are determined, and ridge estimation of the deformation model is implemented. The obtained settlement values ​​are projected onto the vertical direction, and the SAR imaging coordinate system is transformed to the WGS84 geographic coordinate system to generate the deformation rate, cumulative settlement, and coordinates of all points in the geographic coordinate system for settlement analysis.

[0034] Compared with existing technologies, the present invention has the following advantages: The present invention can effectively eliminate model ill-conditionedness, avoid the infinite amplification of observation data errors during the solution process, thereby improving the solution accuracy and effectively monitoring the time series subsidence of the region; whether using RMSE or MAE for evaluation, the improved HK method of the present invention is smaller than the existing Median and Max methods. The present invention improves RMSE by 12.52% and MAE by 12.25% compared with the unimproved HK method, indicating that the present invention greatly improves the solution accuracy of ridge estimation of SBAS-InSAR deformation model and the accuracy of SBAS-InSAR technology in monitoring ground subsidence. Attached Figure Description

[0035] Figure 1 This is a flowchart illustrating the technical process of the present invention.

[0036] Figure 2 The cumulative settlement map of the study area obtained using the MAX method;

[0037] Figure 3 Cumulative settlement map of the study area obtained using the MED method;

[0038] Figure 4 Cumulative settlement map of the study area obtained using the MEAN method;

[0039] Figure 5 The root mean square error map of the study area obtained using the MAX method;

[0040] Figure 6 The root mean square error map of the study area obtained using the MED method;

[0041] Figure 7 The root mean square error map of the study area obtained using the MEAN method;

[0042] Figure 8 This is the result of a time series calculation of the settlement value for feature point 1;

[0043] Figure 9 The result is the time series calculation of the settlement value for feature point 2. Detailed Implementation

[0044] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention are described clearly and completely below. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0045] An improved settlement monitoring method based on the small baseline set deformation model and the HK formula is proposed. This method includes preprocessing the single-view complex image (SLC) and digital elevation model (DEM) of the study area, performing differential interferometry, constructing a small baseline set deformation model, determining ridge parameters, and performing ridge estimation to achieve settlement monitoring.

[0046] Determining ridge parameters includes:

[0047] ;

[0048] In the formula, It is a ridge parameter estimate. and It is the intermediate value to be determined. Indicates to Take the average value.

[0049] Determining ridge parameters includes:

[0050] ;

[0051] ;

[0052] ;

[0053] ;

[0054] In the formula, These are settlement observation values. These are estimates from settlement observations. The number of observations. For the number of estimated values, It is a diagonal matrix composed of eigenvalues. yes eigenvalues, It is the eigenvector corresponding to the eigenvalue. It is the eigenvector matrix. yes The amount, It is the coefficient matrix of the small baseline set deformation model.

[0055] Constructing a small baseline set deformation model includes:

[0056] ;

[0057] In the formula, For the parameters to be determined, It's an error.

[0058] The settlement observation values ​​are:

[0059] ;

[0060] In the formula, It is differential interference phase. It is the radar wavelength.

[0061] The differential interference phase is:

[0062] ;

[0063] In the formula, yes Phase of time, yes The phase of a moment.

[0064] The coefficient matrix is ​​as follows:

[0065] ;

[0066] In the formula, It is the first of the coefficient matrix. line, indicating the first One interference pair, It is the main image time. It is the time for the secondary image.

[0067] The ridge estimation solution criterion is:

[0068] ;

[0069] In the formula, Ridge estimation solution criteria The estimated value, For ridge parameters.

[0070] Differentiating the ridge estimation criterion, we obtain the ridge estimation solution formula:

[0071] ;

[0072] ;

[0073] In the formula, For ridge estimation solution formula The estimated value, This is the weight matrix. As a unit array, and Same order.

[0074] Based on the primary and secondary images of the differential interferometric phase, a coefficient matrix is ​​constructed, ridge parameters are determined, and ridge estimation of the deformation model is implemented. The obtained settlement values ​​are projected onto the vertical direction, and the SAR imaging coordinate system is transformed to the WGS84 geographic coordinate system to generate the deformation rate, cumulative settlement, and coordinates of all points in the geographic coordinate system for settlement analysis.

[0075] The technical flowchart of this invention is as follows: Figure 1 As shown, the process includes data processing, differential interferometry, and time series analysis. Data processing includes SLC extraction, DEM preprocessing, followed by geocoding and registration. Differential interferometry involves collecting interferometric image pairs, performing differential interferometry, filtering, and phase unwrapping based on the initial baseline, then baseline optimization, and finally performing differential interferometry, filtering, and phase unwrapping based on the precise baseline. The time series analysis includes atmospheric delay correction, constructing the SABS deformation model, solving using the HK method, and performing time series sedimentation analysis. The HK method solution includes two processing routes: the first route involves constructing the model coefficient matrix A based on the interferometric image pair combination, estimating the parameters of the deformation model using the least squares (LS) method, and calculating the unit weight mean square error; the second route involves constructing the observation matrix L based on the phase unwrapping results, constructing the deformation model normal equations, and calculating eigenvalues. And the eigenvector matrix C, finally the unit weight mean square error and eigenvalues The ridge parameters under the corrected HK formula are calculated by combining the eigenvector C with the eigenvector C. HK parameter estimation is performed.

[0076] The cumulative settlement maps of the study area obtained using the MAX method (using the maximum value of the eigenvalue when determining ridge parameters), the MED method (using the median of the eigenvalue when determining ridge parameters), and the MEAN method (using the method of this invention when determining ridge parameters) are shown below. Figure 2 , Figure 3 , Figure 4As shown, in the selected study area of ​​this example, the distribution of subsidence areas using the three methods is generally consistent. The MED method shows large-scale subsidence or uplift, with significant subsidence fluctuations throughout the study area. Compared to the MEAN method (the method of this invention), the MAX method detects smaller subsidence values ​​and smaller subsidence areas. Further comparison reveals that in areas with poor ground coherence, such as those with large surface deformation gradients at the center of the subsidence funnel in characteristic regions and areas with significant land cover changes, the MAX solution is close to 0. This means the MAX method cannot effectively detect subsidence data, while the MEAN method (the method of this invention) maintains good solution accuracy even with large observation data errors and can effectively monitor ground subsidence in the aforementioned areas.

[0077] The root mean square error maps of the study area obtained using the MAX method, MED method, and MEAN method (the method of this invention) are shown below. Figure 5 , Figure 6 , Figure 7 As shown, the residuals of the MED method are larger and significantly greater than the other two methods within the study area, indicating that the MED method has poor solution accuracy. Comparative analysis of the other two methods reveals that their residual distribution areas are relatively consistent, but the residuals of the MEAN method (the method of this invention) are significantly smaller than those of the MAX method in some regions. Further analysis shows that both the MED and MAX methods have higher RMSE and MAE than the improved HK method of this invention, indicating lower solution accuracy. MEAN improves RMSE by 12.52% and MAE by 12.25% compared to MAX, demonstrating that this invention effectively improves the solution accuracy of the SBAS-InSAR deformation model.

[0078] Figure 2 , Figure 3 , Figure 4 , Figure 5 , Figure 6 , Figure 7 Each section contains a red box, within which are the centers of two settlement funnels. A characteristic point is selected from each funnel, and their time-series settlement curves are plotted as follows: Figure 8 and Figure 9 As shown, the settlement data obtained by MAX are all close to 0 mm, indicating that the MAX method cannot effectively monitor the time-series settlement in this area. The settlement curves of the MED and MEAN methods show a consistent overall trend, but the MED method exhibits greater fluctuations, with RMSE values ​​of 19.1197 mm and 16.3399 mm for the two feature points, respectively. The settlement curve of the MEAN method is relatively smooth and stable, with smaller errors, and RMSE values ​​of 7.2282 mm and 5.6891 mm for the two feature points, respectively. This demonstrates that the improved HK method of this invention can not only accurately monitor ground settlement but also achieves higher accuracy compared to other methods, showing certain application benefits.

[0079] This invention utilizes a novel HK method to determine ridge parameters and obtain ground subsidence information through deformation model solving, among other steps. This avoids the influence of the last specific minimal eigenvalue in the deformation model matrix on the ridge parameters, effectively improving the model's ill-conditioned nature without requiring extensive prior calculations. Compared to existing HK formulas, this invention improves the model's solution accuracy and effectively monitors the time-series subsidence in the region.

[0080] During the preprocessing of single-view complex imagery (SLC) and DEM for the study area, Sentinel-1 SLC data was extracted and subjected to multi-view processing to remove speckle noise; DEM data was acquired, and its missing areas were interpolated. Next, both types of data were initially geocoded to generate a correspondence file between the SAR imaging coordinate system and the WGS-84 coordinate system, while simultaneously simulating a SAR image in the WGS-84 coordinate system. Then, the simulated SAR image was transformed to the SAR coordinate system, and the offset and offset polynomial of corresponding points between the simulated and real SAR images were calculated using an intensity cross-correlation algorithm. Finally, the initial correspondence was optimized based on the offset polynomial to generate a refined query list. A primary image was selected as the reference, and the remaining images were precisely registered with the primary image, with phase gradient removal performed on both the primary and secondary images. Finally, through multiple iterations, the SAR image registration accuracy was improved to better than 0.001 pixels. After registration, image cropping was performed according to the location and research content of the study area.

[0081] In differential interferometry, a spatiotemporal baseline threshold is set to combine all SAR images into interferometric pairs. The spatial baseline is generally set to one-third of the limiting baseline length, and the temporal baseline, for short wavelengths, is generally set to the time interval of the same season. Based on the initial baseline, differential interferometry is performed to obtain the interferometric phase of the interferometric pairs. Adaptive filtering and minimum flow methods are used for filtering and phase unwrapping. Then, based on the initial phase unwrapping results, baseline optimization is performed using ground control point data, followed by another round of differential interferometry, filtering, and phase unwrapping to obtain accurate phase unwrapping information.

[0082] In Matlab software, GACOS data was used to remove the atmospheric delay phase from the unwrapped phase, obtaining the accurate ground subsidence phase. An SBAS-InSAR deformation model was constructed based on the combination of interferometric pairs. Because the Sentinel-1A satellite has high orbit determination accuracy and maintains interferometric baselines well, all interferometric pairs typically belong to the same subset. Based on these factors, a deformation model was constructed.

[0083] In the normal matrix of the SBAS-InSAR deformation model, the eigenvalues ​​decrease progressively: that is, the h-th eigenvalue is always less than the previous h-1 eigenvalues. However, statistical analysis of the eigenvalues ​​of the normal matrices of multiple deformation models reveals that the difference between the smallest eigenvalue at the end and the previous eigenvalue is much greater than the decreasing gradient corresponding to the other eigenvalues. Therefore, when using the conventional HK formula to determine the ridge parameters of the SBAS-InSAR deformation model, it is often affected by the last minimal eigenvalue, resulting in smaller ridge parameters with certain peculiarities, which cannot effectively improve the ill-conditioned nature of the model and thus affect the accuracy of ground settlement calculation. Table 1 shows the analysis of the deformation model in the SBAS-InSAR technology case, displaying all its eigenvalues.

[0084] Table 1. Eigenvalues ​​of the Deformation Model

[0085] ;

[0086] In this embodiment of the invention, A is a coefficient matrix determined by the interferometric pair combination method. Each row corresponds to one interferometric pair, and each column corresponds to one image. In each interferometric pair, the column containing the main image is 1, and the column containing the secondary image is -1.

[0087] ;

[0088] In practical data processing, this invention uses 67 images from the Sentinel-1A satellite from January 3, 2018 to May 31, 2024 as the data source. Based on the latitude and longitude of the study area, 30m resolution SAR MDI / EM data is downloaded, converted to a GAMMA-compatible DEM format, and a DEM parameter file is generated. Interpolation is performed on missing areas. Temporal and spatial baseline thresholds are set to combine all SAR images into interferometric pairs. The spatial baseline is generally set to one-third of the limiting baseline length, and the temporal baseline, for short wavelengths, is generally set to the time interval of the same season. In this application example, a temporal baseline threshold of 84d and a spatial baseline threshold of 200m are set to combine interferometric image pairs, generating a total of 131 differential interferometric pairs. To reduce the influence of decorrelation noise and make the phase information more continuous and consistent, a 5:1 multiview ratio is set for multiview processing. An adaptive filtering method is used, with a filtering window of 32×32 pixels and a coherence threshold of 0.3, for differential interferometry and filtering. Because the phase difference obtained in the differential interferogram is actually the principal value of the phase, in... To restore the true phase difference, which varies periodically, phase unwrapping is required. This involves adding or subtracting an integer multiple of the phase period (2) from the principal value. Only then can the elevation and deformation of ground targets be successfully deduced. This application example uses the minimum cost flow method to perform phase unwrapping on regions that meet the coherence threshold. Based on the initial phase unwrapping results, baseline optimization is performed using ground control point data extracted from the DEM to obtain a precise baseline. Finally, differential interferometry, filtering, and phase unwrapping are performed again based on the precise baseline to obtain accurate phase unwrapping information.

[0089] In this application example, the eigenvalues ​​of the SBAS-InSAR deformation model method matrices constructed in six schemes were statistically analyzed. All eigenvalues ​​gradually decreased, but the smallest eigenvalue decreased by approximately 10 compared to the previous eigenvalue. -15 The gradient is significantly larger than the decreasing gradient between the first 30 eigenvalues. Therefore, when determining the ridge parameters of the SBAS-InSAR deformation model based on the unmodified HK formula, it is often affected by the last minimum eigenvalue, resulting in smaller ridge parameters with certain peculiarities, which cannot effectively improve the ill-conditioned nature of the model. Therefore, this invention improves the conventional HK formula, thereby... Calculate the ridge parameters by taking the average value.

[0090] Based on the ridge estimation solution of the SBAS-InSAR deformation model, i.e., the time series settlement data of each point, the obtained settlement data is projected from the radar line of sight to the vertical direction, and transformed from the SAR imaging coordinate system to the WGS84 geographic coordinate system. Deformation rate, cumulative settlement, and coordinate text data of all points in the geographic coordinate system are generated. The transformed results are output as raster or vector format files and can be imported into other software for further analysis.

[0091] To facilitate comparison of the effectiveness of the new method in determining ridge parameters for ridge estimation, this invention uses 67 images as the data source and processes the data in a specific study area using the above method to obtain time-series subsidence data for that area. To better demonstrate the advantages of this method, a comparative analysis is conducted using the new HK formula method (1), the original HK method (2), and another HK method using the improved MED formula (3). The three HK methods are represented by the following formulas:

[0092] (1);

[0093] (2);

[0094] (3);

[0095] In the formula, Indicates to Take the maximum value. Indicates to Take the median.

[0096] The root mean square error (RMSE) and mean absolute error (MAE) of the three methods are shown in Table 2.

[0097] Table 2 Root Mean Square Error and Mean Absolute Error

[0098] .

[0099] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A settlement monitoring method for solving the HK formula by improving the small baseline set deformation model, characterized in that, The single look complex (SLC) and digital elevation model (DEM) of the study area are preprocessed, and differential interference is performed; A small baseline subset deformation model is constructed, ridge parameters are determined, and ridge estimation calculation is performed to realize settlement monitoring; The determination of the ridge parameters includes: ； wherein is a ridge parameter estimate, and is an intermediate value to be solved, denotes taking the average over the average value; The determination of the ridge parameters includes: ； ； ； ； In the formula, These are settlement observation values. These are estimates of the settlement observations. The number of observations. For the number of estimated values, It is a diagonal matrix composed of eigenvalues. yes eigenvalues, It is the eigenvector corresponding to the eigenvalue. It is the eigenvector matrix. yes The amount, It is the coefficient matrix of the small baseline set deformation model; The construction of the small baseline subset deformation model includes: ； wherein is the parameter to be found, is the error; The settlement observation value is: ； wherein is the differential interference phase, is the radar wavelength; The differential interference phase is: ； wherein is the phase of the time instant, is the phase of the time instant; The coefficient matrix is: ； wherein is the coefficient matrix of the first row, representing the first interfering pair, is the main image time, is the sub image time; The ridge estimation calculation criterion is: ； In the formula, is the estimate value of under the ridge estimation solution criterion is the estimate value of under the ridge estimation solution criterion is the ridge parameter; The ridge estimation calculation formula is obtained by derivation of the ridge estimation calculation criterion: ； ； In the formula, is the ridge estimate solution formula is the estimate value of is the weight matrix, is the unit matrix, and is the same order.

2. The method according to claim 1, wherein, According to the main image and the auxiliary image of the differential interference phase, the coefficient matrix is constructed, the ridge parameters are determined, and the ridge estimation calculation of the deformation model is realized. The obtained settlement value is projected to the vertical direction, converted from the SAR imaging coordinate system to the WGS84 geographic coordinate system, and the deformation rate, cumulative settlement and coordinates of all points in the geographic coordinate system are generated to perform settlement analysis.

Images