Settlement monitoring method based on improved HK formula solution of small baseline set deformation model

Abstract

The invention discloses a settlement monitoring method for improving HK formula solution by a small baseline set deformation model, belongs to the technical field of settlement monitoring, is used for settlement monitoring, and comprises the following steps: preprocessing a single-view complex image SLC and a digital elevation model DEM of a research area, and carrying out differential interference; a small baseline set deformation model is constructed, ridge parameters are determined, ridge estimation calculation is carried out, and settlement monitoring is achieved; according to the invention, good calculation precision can be maintained when the observation data error is large, and the time sequence settlement condition of the experiment area can be effectively monitored; whether the RMSE or the MAE is used for evaluation, the improved HK method is smaller than an existing method, compared with an unimproved HK method, the RMSE of the improved HK method is improved by 12.52%, the MAE of the improved HK method is improved by 12.25%, and therefore the calculation precision of ridge estimation of the SBAS-InSAR deformation model and the ground monitoring precision of the SBAS-InSAR technology are greatly improved.

Description

Technical Field

[0001] The invention discloses a settlement monitoring method for solving HK formula with a small baseline set deformation model, and belongs to the technical field of settlement monitoring. Background Art

[0002] The existing technology usually uses SBAS-InSAR technology to build a deformation model for settlement prediction. In the normal matrix of the SBAS-InSAR deformation model, the eigenvalue gradually decreases, that is, the hth eigenvalue is always smaller than the first h-1 eigenvalues, but the difference between the last smallest eigenvalue and the previous eigenvalue is much larger than the corresponding decrease gradient of other eigenvalues. In view of the particularity of the eigenvalues ​​of the deformation model, when the conventional HK formula is used to determine the ridge parameters of the SBAS-InSAR deformation model, it is often affected by the last minimum eigenvalue, and the obtained ridge parameters are small, which affects the improvement of the model's pathological properties, and thus also has a certain degree of impact on the acquisition of ground settlement results. Summary of the invention

[0003] The purpose of the present invention is to provide a settlement monitoring method for solving the HK formula with a small baseline set deformation model, so as to solve the problem in the prior art that when the conventional HK formula is used to determine the ridge parameters, the ridge parameters obtained are small and the pathological nature of the model cannot be improved well, thus affecting the settlement monitoring results.

[0004] The small baseline set deformation model is used to improve the settlement monitoring method solved by HK formula, including the single-view complex image SLC and digital elevation model DEM preprocessing of the study area, and differential interference; the small baseline set deformation model is constructed, the ridge parameters are determined and the ridge estimation solution is performed to realize settlement monitoring.

[0005] Determining the ridge parameters involves: ; In the formula, is the ridge parameter estimate, and is the intermediate value to be found, Express Take the average.

[0006] Determining the ridge parameters involves: ; ; ; ; In the formula, is the settlement value to be calculated, is an estimate of the settlement to be calculated, is the number of observations, is the number of estimated values, is a diagonal matrix of characteristic roots, yes The characteristic root of is the eigenvector corresponding to the eigenvalue, is the eigenvector matrix, yes The weight, is the coefficient matrix of the small baseline set deformation model.

[0007] Building a small baseline set deformation model includes: ; In the formula, is the parameter to be requested, It's an error.

[0008] The settlement value to be calculated is: ; In the formula, is the differential interference phase, is the radar wavelength.

[0009] The differential interference phase is: ; In the formula, yes The phase of the moment, yes The phase of the moment.

[0010] The coefficient matrix is: ; In the formula, is the coefficient matrix row, indicating the An interference pair, is the main image time, It is the secondary image time.

[0011] The ridge estimate solution is: ; In the formula, It is the ridge estimate solution under the The estimated value of is the ridge parameter.

[0012] Taking the derivative of the ridge estimation solution, we get the ridge estimation solution: ; ; In the formula, is the ridge estimate solution The estimated value of is the weight matrix, is the unit matrix, and Same order.

[0013] According to the main image and the secondary image of the differential interferometry phase, a coefficient matrix is ​​constructed to determine the ridge parameters and realize the ridge estimation solution of the deformation model. The obtained settlement value is projected to the vertical direction and converted from the SAR imaging coordinate system to the WGS84 geographic coordinate system. The deformation rate, cumulative settlement amount and coordinates of all points in the geographic coordinate system are generated for settlement analysis.

[0014] Compared with the prior art, the present invention has the following beneficial effects: the present invention can effectively eliminate the model pathology and avoid the observation data error being infinitely amplified during the solution process, thereby improving the solution accuracy and effectively monitoring the time series settlement of the area; whether using RMSE or MAE evaluation, the improved HK method of the present invention is smaller than the existing Median and Max methods. Compared with the unimproved HK method, the present invention improves RMSE by 12.52% and MAE by 12.25%, indicating that the present invention greatly improves the solution accuracy of the ridge estimation of the SBAS-InSAR deformation model and the accuracy of SBAS-InSAR technology in monitoring ground settlement. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 It is a technical flow chart of the present invention; Figure 2 The cumulative settlement map obtained using the MAX method for the study area; Figure 3 The cumulative settlement map obtained using the MED method for the study area; Figure 4 The cumulative settlement map obtained using the MEAN method for the study area; Figure 5 Root mean square error map obtained using the MAX method for the study area; Figure 6 Root mean square error map obtained using the MED method for the study area; Figure 7 Root mean square error map obtained using the MEAN method for the study area; Figure 8 It is the result of the settlement value time series calculation for feature point 1; Fig. 9 It is the result of the time series calculation of the settlement value for feature point 2. DETAILED DESCRIPTION

[0016] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention is described clearly and completely below. Obviously, the described embodiments are part of the embodiments of the present invention, but not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0017] The small baseline set deformation model is used to improve the settlement monitoring method solved by HK formula, including the single-view complex image SLC and digital elevation model DEM preprocessing of the study area, and differential interference; the small baseline set deformation model is constructed, the ridge parameters are determined and the ridge estimation solution is performed to realize settlement monitoring.

[0018] Determining the ridge parameters involves: ; In the formula, is the ridge parameter estimate, and is the intermediate value to be found, Express Take the average.

[0019] Determining the ridge parameters involves: ; ; ; ; In the formula, is the settlement value to be calculated, is an estimate of the settlement to be calculated, is the number of observations, is the number of estimated values, is a diagonal matrix of characteristic roots, yes The characteristic root of is the eigenvector corresponding to the eigenvalue, is the eigenvector matrix, yes The weight, is the coefficient matrix of the small baseline set deformation model.

[0020] Building a small baseline set deformation model includes: ; In the formula, is the parameter to be requested, It's an error.

[0021] The settlement value to be calculated is: ; In the formula, is the differential interference phase, is the radar wavelength.

[0022] The differential interference phase is: ; In the formula, yes The phase of the moment, yes The phase of the moment.

[0023] The coefficient matrix is: ; In the formula, is the coefficient matrix row, indicating the An interference pair, is the main image time, It is the secondary image time.

[0024] The ridge estimate solution is: ; In the formula, It is the ridge estimate solution under the The estimated value of is the ridge parameter.

[0025] Taking the derivative of the ridge estimation solution, we get the ridge estimation solution: ; ; In the formula, is the ridge estimate solution The estimated value of is the weight matrix, is the unit matrix, and Same order.

[0026] According to the main image and the secondary image of the differential interferometry phase, a coefficient matrix is ​​constructed to determine the ridge parameters and realize the ridge estimation solution of the deformation model. The obtained settlement value is projected to the vertical direction and converted from the SAR imaging coordinate system to the WGS84 geographic coordinate system. The deformation rate, cumulative settlement amount and coordinates of all points in the geographic coordinate system are generated for settlement analysis.

[0027] The technical flow chart of the present invention is as follows: Figure 1As shown, it includes data processing, differential interferometry and time series. Data processing includes SLC extraction, DEM preprocessing, and then geocoding and registration. Differential interferometry includes collecting interference image pair combinations, performing differential interferometry, filtering and phase unwrapping based on the initial baseline, and then optimizing the baseline. Finally, differential interferometry, filtering and phase unwrapping are performed based on the precise baseline. Time series includes atmospheric delay correction, construction of SABS deformation model, HK method solution, and time series sedimentation analysis. The HK method solution includes two processing routes. The first one includes constructing the model coefficient matrix A based on the interference pair combination, using the least squares (LS) method to estimate the parameters of the deformation model and calculate the unit weighted mean error. The second one includes constructing the observation matrix L based on the phase unwrapping results, constructing the deformation model equation and calculating the eigenvalue. And the eigenvector matrix C, finally the unit weight error and eigenvalue Combined with the eigenvector C, the ridge parameter corrected by HK formula is calculated , perform HK parameter estimation.

[0028] The cumulative settlement maps obtained by using the MAX method (using the maximum value of the characteristic part when determining the ridge parameters), the MED method (using the median of the characteristic part when determining the ridge parameters), and the MEAN method (using the method of the present invention when determining the ridge parameters) in the study area are as follows: Figure 2 , Figure 3 , Figure 4 As shown in the figure, in the selected study area of ​​this example, the settlement area distribution of the three methods is generally consistent. Among them, the MED method shows large-scale settlement or uplift, and the settlement fluctuates greatly in the entire study area. Compared with the MEAN method (the method of the present invention), the MAX method monitors smaller settlement values ​​and smaller settlement areas. Further comparison shows that in some areas with poor ground coherence such as the center of the settlement funnel in the characteristic area with large surface deformation gradient and obvious land cover changes, the MAX solution result is close to 0, that is, the MAX method cannot effectively detect settlement data, while the MEAN method (the method of the present invention) can maintain good solution accuracy when the observation data error is large, and can effectively monitor the ground settlement in the above-mentioned areas.

[0029] The root mean square error diagrams obtained by using the MAX method, the MED method, and the MEAN method (method of the present invention) in the study area are shown as follows: Figure 5 , Figure 6 , Figure 7As shown, the residual of the MED method is large and significantly larger than the other two methods in the study area, indicating that the solution accuracy of this method is poor. Comparative analysis of the other two methods found that: the residual distribution areas are relatively consistent, but the residual of the solution result of the MEAN method (the method of the present invention) is significantly smaller than that of the MAX method in some areas. Further analysis found that both the RMSE and MAE of the MED and MAX methods are greater than the improved HK method of the invention, and the solution accuracy is low. Compared with MAX, MEAN has improved by 12.52% in RMSE and 12.25% in MAE, indicating that the invention has effectively improved the solution accuracy of the SBAS-InSAR deformation model.

[0030] Figure 2 , Figure 3 , Figure 4 , Figure 5 , Figure 6 , Figure 7 There are red frame parts in the box. There are two sedimentation funnel centers. Select a characteristic point for each and draw their time series sedimentation curves. Figure 8 and Fig. 9 As shown. It can be found that the settlement data of each period obtained by MAX are close to 0mm, indicating that the MAX method cannot effectively monitor the time series settlement of the area. The overall trend of the settlement curves of the MED method and the MEAN method is consistent, but the MED method fluctuates greatly. The RMSE of the two feature points of the MED method are 19.1197mm and 16.3399mm respectively; the settlement curve of the MEAN method is relatively smooth and stable, and the result error is small. The RMSE of the two feature points is 7.2282mm and 5.6891mm respectively. It shows that the improved HK method of the invention can not only accurately monitor the ground settlement, but also has a higher accuracy than other methods, and has a certain application effect.

[0031] The present invention uses a new HK method to determine the ridge parameters and solve the deformation model to obtain ground subsidence information, avoiding the influence of the last special minimum eigenvalue in the deformation model matrix on the ridge parameters, effectively improving the model's pathological nature, and without the need for a large amount of prior calculations. Compared with the existing HK formula, the invention can improve the model's solution accuracy and effectively monitor the time series subsidence of the area.

[0032] When preprocessing the single-look complex image (SLC) and DEM of the study area, the Sentinel-1 SLC data was extracted and multi-looked to remove its speckle noise; the DEM data was obtained and the missing areas were interpolated. Next, the two data were initially geocoded to generate a correspondence file between the SAR imaging coordinate system and the WGS-84 coordinate system, and a SAR image in the WGS-84 coordinate system was simulated. Then the simulated SAR image was converted to the SAR coordinate system, and the offset and offset polynomial of the same-name points between the simulated SAR image and the real SAR image were calculated according to the intensity cross-correlation algorithm. Finally, the initial correspondence was optimized according to the offset polynomial to generate a refined query list. A main image was selected as the benchmark, and the remaining images were accurately registered with the main image, and the phase gradient of the main image and the secondary image was removed. Finally, the SAR image registration accuracy was made better than 0.001 pixel through multiple iterations. After registration, the image can be cropped according to the location and research content of the study area.

[0033] During differential interferometry, the temporal and spatial baseline threshold is set to combine all SAR images into an interferometric pair. The spatial baseline is generally set to one-third of the limit baseline length, and the temporal baseline is generally set to the time interval of the same season for short wavelengths. Based on the initial baseline, differential interferometry is performed to obtain the interferometric phase of the interferometric pair, and the adaptive filtering method and the minimum cost flow method are used to filter and phase unwrap it. Then, based on the initial phase unwrapping result, the baseline is optimized using the ground control point data, and differential interferometry, filtering and phase unwrapping are performed again to obtain accurate phase unwrapping information.

[0034] In Matlab software, GACOS data is used to remove the atmospheric delay phase in the unwrapped phase to obtain the accurate ground settlement phase. The SBAS-InSAR deformation model is constructed according to the combination of interference pairs. Because the orbit determination accuracy of the Sentinel-1A satellite is high and the interference baseline is well maintained, usually all interference pairs belong to the same subset. Based on the above situation, the deformation model is constructed.

[0035] In the normal matrix of the SBAS-InSAR deformation model, the eigenvalues ​​decrease step by step: that is, the hth eigenvalue is always smaller than the previous h-1 eigenvalues. However, by counting the eigenvalues ​​of the normal matrices of multiple deformation models, it is found that the difference between the smallest eigenvalue at the end and the previous eigenvalue is much larger than the decreasing gradient corresponding to the other eigenvalues. Therefore, when the conventional HK formula is used to determine the ridge parameters of the SBAS-InSAR deformation model, it is often affected by the last minimum eigenvalue. The obtained ridge parameters are small and have certain particularities, which cannot improve the model's pathological nature well, thereby affecting the accuracy of ground subsidence calculations. As shown in Table 1, the deformation model in the SBAS-InSAR technology case is analyzed and all its eigenvalues ​​are displayed.

[0036] Table 1 Characteristic values ​​of deformation model ; In the embodiment of the present invention, A is a coefficient matrix determined by the interference pair combination method, each row corresponds to an interference pair, each column corresponds to an image, and the column where the main image in each interference pair is located is 1, and the column where the secondary image is located is -1: ; In the actual data processing of the present invention, a total of 67 images from the Sentinel-1A satellite from January 3, 2018 to May 31, 2024 are used as the data source. According to the longitude and latitude of the study area, SRTMDEM data with a resolution of 30m is downloaded, converted into the dem format applicable to GAMMA, and a dem parameter file is generated, and the missing areas are interpolated. The time baseline and space baseline thresholds are set to combine all SAR images into interference pairs. The space baseline is generally set to one-third of the limit baseline length, and the time baseline is generally set to the time interval of the same season for short wavelengths. This application example sets the time baseline threshold of 84d and the space baseline threshold of 200m to combine interference image pairs, and a total of 131 differential interference pairs are generated. In order to weaken the influence of decorrelation noise and make the phase information more continuous and consistent, this example sets a multi-view ratio of 5:1 for multi-view processing. An adaptive filtering method is used with a filter window of 32×32 pixels and a coherence threshold of 0.3 for differential interference and filtering. Because the phase difference obtained in the differential interference diagram is actually the main value of the phase, in It changes periodically. To restore the true phase difference, phase unwrapping is required, that is, adding or subtracting integer multiples of the phase period 2 on the basis of the main value. , in order to successfully invert the elevation and deformation of the ground target. This application example uses the minimum cost flow method to perform phase unwrapping on the area that meets the coherence threshold. Based on the initial phase unwrapping results, the baseline is optimized in combination with the 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.

[0037] In this application example, the eigenvalues ​​of the SBAS-InSAR deformation model matrix constructed in the six schemes are counted. All eigenvalues ​​gradually decrease, but the smallest eigenvalue decreases by about 10 compared with the previous eigenvalue. -15 times, which is much larger than the decreasing gradient between the first 30 eigenvalues. Therefore, when the ridge parameters of the SBAS-InSAR deformation model are determined based on the unimproved HK formula, they are often affected by the last minimum eigenvalue. The obtained ridge parameters are small and have certain particularity, which cannot improve the pathological nature of the model. The ridge parameter is calculated by taking the average.

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

[0039] In order to compare the effect of the new method on determining ridge parameters for ridge estimation, the present invention uses 67 images as data sources, uses the above method to process data based on a certain study area, and obtains the time series settlement data of the study area. In order to better demonstrate the advantages of this method, the new HK formula method of formula (1) is used for comparative analysis with the original HK method of formula (2) and another HK method of formula (3) improved by MED. The three HK methods are characterized by the following formulas: (1); (2); (3); In the formula, Express Take the maximum value, Express Take the median.

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

[0041] Table 2 Root mean square error and mean absolute error .

[0042] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the aforementioned embodiments, a person skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some or all of the technical features may be replaced by equivalents, and these modifications or replacements 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 based on the small baseline set deformation model and improved HK formula solution, characterized in that: It includes preprocessing the single-view complex image SLC and digital elevation model DEM of the study area, performing differential interference; building a small baseline set deformation model, determining ridge parameters and performing ridge estimation solution to achieve settlement monitoring.

2. The settlement monitoring method of claim 1 using the small baseline set deformation model and improved HK formula solution is characterized in that: Determining the ridge parameters involves: ； In the formula, is the ridge parameter estimate, and is the intermediate value to be found, Express Take the average.

3. The settlement monitoring method of claim 2 using the small baseline set deformation model and improved HK formula solution is characterized in that: Determining the ridge parameters involves: ； ； ； ； In the formula, is the settlement value to be calculated, is an estimate of the settlement to be calculated, is the number of observations, is the number of estimated values, is a diagonal matrix of characteristic roots, yes The characteristic root of is the eigenvector corresponding to the eigenvalue, is the eigenvector matrix, yes The weight, is the coefficient matrix of the small baseline set deformation model.

4. The settlement monitoring method of claim 3 using the small baseline set deformation model and improved HK formula solution is characterized in that: Building a small baseline set deformation model includes: ； In the formula, is the parameter to be requested, It's an error.

5. The settlement monitoring method of claim 4 using the small baseline set deformation model and improved HK formula solution is characterized in that: The settlement value to be calculated is: ； In the formula, is the differential interference phase, is the radar wavelength.

6. The settlement monitoring method of claim 5 using the small baseline set deformation model and improved HK formula solution is characterized in that: The differential interference phase is: ； In the formula, yes The phase of the moment, yes The phase of the moment.

7. The settlement monitoring method of claim 6 using the small baseline set deformation model and improved HK formula solution is characterized in that: The coefficient matrix is: ； In the formula, is the coefficient matrix row, indicating the An interference pair, is the main image time, It is the secondary image time.

8. The settlement monitoring method of claim 7 using the small baseline set deformation model and improved HK formula solution is characterized in that: The ridge estimate solution is: ； In the formula, It is the ridge estimate solution under the The estimated value of is the ridge parameter.

9. The settlement monitoring method of claim 8 using the small baseline set deformation model and improved HK formula solution is characterized in that: Taking the derivative of the ridge estimation solution, we get the ridge estimation solution: ； ； In the formula, is the ridge estimate solution The estimated value of is the weight matrix, is the unit matrix, and Same order.

10. The settlement monitoring method of claim 9 using the small baseline set deformation model and improved HK formula solution, characterized in that: According to the main image and the secondary image of the differential interferometry phase, a coefficient matrix is ​​constructed to determine the ridge parameters and realize the ridge estimation solution of the deformation model. The obtained settlement value is projected to the vertical direction and converted from the SAR imaging coordinate system to the WGS84 geographic coordinate system. The deformation rate, cumulative settlement amount and coordinates of all points in the geographic coordinate system are generated for settlement analysis.

Images