A method for atmospheric correction of interferograms in continuous deformation monitoring of ground-based synthetic aperture radar

Abstract

The present invention relates to a method for atmospheric correction of continuous deformation monitoring interference pattern of ground-based synthetic aperture radar, comprising the following steps: step 1: input GB-SAR complex image sequence, and form interference image pair; step 2: based on three parameters of thermal signal-to-noise ratio, amplitude discrete index and coherence coefficient, adopt multi-threshold method to obtain high-quality pixels; step 3: sample and thin out the high-quality pixels; step 4: perform spatial phase unwrapping on the interference phase of the thinned high-quality pixels; step 5: establish a polynomial atmospheric phase model of the survey area, and use the polynomial atmospheric phase model to perform atmospheric correction on the unwrapped phase of the high-quality pixels; step 6: adopt a robust fitting method to solve the parameters of the polynomial atmospheric phase model; step 7: optimize the polynomial atmospheric phase model. The present invention can effectively remove the influence of atmospheric disturbance, thereby improving the processing speed and processing accuracy of radar monitoring, and can be applied to ground-based radar monitoring under complex environmental conditions.

Description

Technical Field

[0001] The invention relates to the field of engineering deformation monitoring, and in particular to an atmospheric correction method for interference graphs of continuous deformation monitoring using a ground-based synthetic aperture radar. Background Art

[0002] Deformation monitoring is one of the important technical means of geological disaster prevention and control. By carrying out deformation monitoring work, potential deformation hazards can be discovered in time to reduce economic losses and casualties. Traditional geodetic deformation monitoring technology has the problem of low monitoring point density, which is easy to miss the key parts of the deformed body and break the continuity of deformation information in space.

[0003] Synthetic Aperture Radar (SAR) is a high-resolution, two-dimensional imaging radar that is particularly suitable for large-area surface imaging and is an important development direction of radar. Since the main band used by its sensor is the microwave band, it has good penetration and can penetrate clouds, fog, rain and snow. It has the ability to work all-weather, and is a cutting-edge technology and research hotspot for deformation monitoring.

[0004] The ground-based synthetic aperture radar interferometry (GB-InSAR) technology has greatly facilitated slope deformation monitoring. Compared with traditional deformation monitoring methods, GB-InSAR can work around the clock and in all weather conditions, allowing large-area non-contact, non-invasive and non-destructive measurements without causing damage to the monitored object, and obtaining deformation results with a wide range, long distance and continuous spatial coverage, saving time and effort.

[0005] The heat exchange in the atmosphere causes the airflow to be in constant change, resulting in different phase delays and noises on the time and space scales. As a result, in the continuous monitoring process of the GB-InSAR dam slope deformation, the monitoring results in a short period of time show large phase fluctuations, which are obviously inconsistent with the actual situation. Therefore, it is necessary to apply atmospheric modeling and correction to the continuous deformation monitoring data of ground-based synthetic aperture radar (GB-SAR) to weaken the impact of atmospheric changes on deformation monitoring results and improve the accuracy of deformation monitoring sequences.

[0006] At present, the commonly used GB-SAR atmospheric correction methods include reference point method, meteorological parameter model method, space-time filtering method and function model method. Among them, the reference point method has low accuracy and is suitable for areas with small monitoring range; the meteorological parameter model method is difficult to obtain meteorological parameters with sufficient density; the window size setting of the space-time filtering method has a greater impact on the results; the function model method can obtain better results when selecting a more appropriate function, but it is also susceptible to the deformation area, resulting in deviations in the fitting results. In addition, the above methods have low calculation efficiency when applied to the processing of large-scale continuous GB-SAR image data, and it is difficult to meet the needs of quasi-real-time deformation monitoring. Summary of the invention

[0007] The present invention provides an atmospheric correction method for continuous deformation monitoring interference pattern of ground-based synthetic aperture radar to solve the above technical problems.

[0008] In order to solve the above technical problems, the present invention provides a method for atmospheric correction of continuous deformation monitoring interferogram of ground-based synthetic aperture radar, comprising the following steps:

[0009] Step 1: Input GB-SAR complex image sequence and form interferometric image pairs;

[0010] Step 2: Based on the three parameters of thermal signal-to-noise ratio, amplitude dispersion index and coherence coefficient, a multi-threshold method is used to obtain high-quality pixels;

[0011] Step 3: Performing sampling and thinning processing on the high-quality pixels;

[0012] Step 4: performing spatial phase unwrapping on the thinned high-quality pixel interference phase;

[0013] Step 5: Establish a polynomial atmospheric phase model for the measurement area, and use the polynomial atmospheric phase model to perform atmospheric correction on the unwrapped phase of the high-quality pixel;

[0014] Step 6: Using a robust fitting method to solve the parameters of the polynomial atmospheric phase model;

[0015] Step 7: Optimize the polynomial atmospheric phase model.

[0016] Preferably, in step 1, the first image is selected from the GB-SAR complex image sequence as the main image, and the other images are selected as slave images to form the interference image pair.

[0017] Preferably, step 1 further comprises: filtering and smoothing each interference pattern of the interference image pair.

[0018] Preferably, step 2 further comprises: calculating a mask matrix based on the high-quality pixels, performing mask processing on the filtered interference pattern, and obtaining phase values ​​of discrete high-quality pixels.

[0019] Preferably, in step 3, a uniform grid sampling method is used to perform sampling thinning processing on the high-quality pixels, including: dividing the radar image of the survey area into a square grid, selecting the pixel with the smallest coordinate from the center of each grid as the extracted pixel, and updating the mask matrix.

[0020] Preferably, the side length of the divided grid is less than 1 / 3 of the width of a single complete period interference fringe.

[0021] Preferably, step 4 comprises:

[0022] Step 41: using the updated mask matrix to perform mask calculation on the filtered interference pattern to obtain an interference pattern sequence of discrete high-quality pixels;

[0023] Step 42: Select a high-quality pixel as the starting point of the phase gradient integration, and use the minimum cost flow method to perform spatial phase unwrapping on the thinned interference phase of the high-quality pixel for each interference pattern.

[0024] Preferably, the polynomial atmospheric phase model is a univariate polynomial model of distance, or a bivariate polynomial model of radar plane coordinates.

[0025] Preferably, the robust fitting method includes the least absolute deviation method and the double square method.

[0026] Preferably, step 7 includes: using the unwrapped interferogram after atmospheric correction to calculate the initial linear deformation rate of the high-quality pixels, obtaining the deformation area based on the initial linear deformation rate, eliminating the high-quality pixels in the deformation area, updating the mask matrix, and recalculating the parameters of the polynomial atmospheric phase model of the high-quality pixels in the remaining area.

[0027] Compared with the prior art, the atmospheric correction method for continuous deformation monitoring interferogram of ground-based synthetic aperture radar provided by the present invention has the following advantages:

[0028] 1. The present invention can effectively remove the influence of atmospheric disturbance, thereby improving the processing speed and processing accuracy of radar monitoring, and can be applied to ground-based radar monitoring under complex environmental conditions;

[0029] 2. The present invention greatly improves the efficiency of interferogram phase unwrapping and subsequent atmospheric phase modeling and correction by thinning high-quality pixels, and is particularly suitable for GB-SAR continuous processing and quasi-real-time processing modes that require high timeliness of deformation solution;

[0030] 3. The present invention establishes a polynomial atmospheric phase model of the survey area based on the unwrapped phase of the thinned high-quality pixels. Compared with the reference point method, this method performs polynomial atmospheric correction on the basis of spatial phase unwrapping and does not require prior information of stable areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 The present invention is a flowchart of a method for atmospheric correction of a ground-based synthetic aperture radar continuous deformation monitoring interferogram in a specific implementation manner of the present invention. DETAILED DESCRIPTION

[0032] In order to describe the technical solution of the above invention in more detail, specific embodiments are listed below to demonstrate the technical effects; it should be emphasized that these embodiments are used to illustrate the present invention but not to limit the scope of the present invention.

[0033] The present invention provides a method for atmospheric correction of continuous deformation monitoring interference pattern of ground-based synthetic aperture radar, such as Figure 1 As shown, the following steps are included:

[0034] Step 1: Input GB-SAR complex image sequence and form interferometric image pairs;

[0035] Step 2: Based on the three parameters of thermal signal-to-noise ratio, amplitude dispersion index and coherence coefficient, a multi-threshold method is used to obtain high-quality pixels;

[0036] Step 3: Performing sampling and thinning processing on the high-quality pixels;

[0037] Step 4: performing spatial phase unwrapping on the thinned high-quality pixel interference phase;

[0038] Step 5: Establish a polynomial atmospheric phase model for the measurement area, and use the polynomial atmospheric phase model to perform atmospheric correction on the unwrapped phase of the high-quality pixel;

[0039] Step 6: Using a robust fitting method to solve the parameters of the polynomial atmospheric phase model;

[0040] Step 7: Optimize the polynomial atmospheric phase model.

[0041] Based on the reference point method, the present invention considers the temporal and spatial distribution characteristics of atmospheric disturbances and establishes a polynomial atmospheric phase model based on the thinned high-quality pixels. This method can effectively remove the influence of atmospheric disturbances and improve the processing speed. In addition, compared with the univariate polynomial model, the bivariate polynomial model has higher accuracy. Therefore, this method can be applied to ground-based radar monitoring under complex environmental conditions and can provide technical support for application scenarios such as landslide risk assessment.

[0042] In some embodiments, in step 1, the first image is selected from the GB-SAR complex image sequence as the main image, and the other images are selected as slave images to form the interference image pair, and the interference calculation is performed according to formula (1).

[0043] Intf=M(i,j)S * (i,j) (1)

[0044] Where Intf represents the complex matrix of the interference pattern, M and S represent the master and slave GB-SAR images of the interference pair, respectively, i and j represent the row and column numbers of the pixels in the image, and * represents the conjugate multiplication of the two phase matrices.

[0045] In order to reduce the interference phase noise, improve the quality of the interference pattern, and maintain the interference fringes, step 1 of this embodiment also includes: filtering and smoothing each interference pattern of the interference image pair. Specifically, each interference pattern is filtered using the Goldstein filtering method according to formula (2). The interference phase pattern is transformed into the spectrum domain using Fourier transform, and then the interference filter is smoothed. In the high coherence area, the filter coefficient is small, and the filter strength is small; in the low coherence area, the filter coefficient is large, and the filter strength is large.

[0046]

[0047]

[0048] In the formula, Φ(i,j) represents the phase after filtering, and S(u,v) represents the interference phase diagram The subscript m indicates smoothing; |·| indicates modulo, ⊙ indicates the Hadamard product of matrices, that is, the element-by-element multiplication of two matrices of the same dimension; α∈[0,1] indicates the filter parameter, which is generally 0.5.

[0049] Since there are a large number of false signals and weak signal pixels in GB-SAR image data, the accuracy of GB-SAR image interpretation and interferometry is directly affected. The single threshold method based on thermal signal-to-noise ratio, amplitude dispersion index and coherence coefficient has defects to varying degrees in selecting high-quality pixels. Therefore, the present invention adopts a multi-threshold method to select high-quality pixels (HQP) based on the calculation of three parameters: thermal signal-to-noise ratio, amplitude dispersion index and coherence coefficient.

[0050] (1) Thermal signal-to-noise ratio

[0051] By calculating the thermal signal-to-noise ratio (TSNR) of each pixel, the pixels with a TSNR lower than a given threshold are regarded as false signals and removed.

[0052]

[0053] Where S is the signal power, N t =2σ n 2 is the receiver thermal noise, σ n is the mean error of the receiver thermal noise, S m is the complex number of pixels.

[0054] (2) Coherence coefficient

[0055] High-quality pixels are selected by taking advantage of the premise that the radar scattering characteristics do not change with time. The coherence coefficient of the target pixel is estimated based on the information of neighboring pixels, and points within a certain range that are greater than the set threshold are extracted and retained as high-quality pixels.

[0056]

[0057] Where γ represents the coherence coefficient, M and S represent the master image and slave image of the interference pair, respectively, m and n represent the window size, and * represents the conjugate multiplication of two phase matrices.

[0058] (3) Amplitude dispersion index

[0059] In complex images, the amplitude dispersion index (ADI) is calculated based on the standard deviation of the real and imaginary parts and the time-series amplitude mean.

[0060]

[0061] In the formula, σ A is the time series amplitude standard deviation, m A is the time series amplitude mean, D A is the amplitude dispersion index.

[0062] Appropriate thresholds are set for the thermal signal-to-noise ratio, amplitude dispersion index, and coherence coefficient. When a pixel satisfies the three threshold conditions at the same time, the pixel is considered to be HQP.

[0063] In some embodiments, step 2 further includes: calculating a mask matrix based on the high-quality pixels, performing mask processing on the filtered interference pattern, and obtaining discrete HQP phase values.

[0064] GB-SAR has a high temporal and spatial resolution, resulting in a huge amount of image data. Even after multi-threshold pixel screening, there are still a large number of HQPs, which directly affects the processing efficiency of subsequent spatial phase unwrapping, atmospheric modeling and correction.

[0065] In some embodiments, the present invention uses a uniform grid sampling method to perform sampling thinning processing on the high-quality pixels, including: dividing the radar image of the survey area into a square grid, selecting the pixel with the smallest coordinate from the center of each grid as the extracted pixel, and updating the mask matrix after completing the HQP thinning. In some embodiments, the side length of the divided grid can be determined based on the intensity of the atmospheric change of the interference pattern, that is, the density of the interference pattern fringes, which is usually less than 1 / 3 of the width of a single complete period interference fringe.

[0066] Pixel thinning processing can greatly improve the efficiency of interferogram phase unwrapping and subsequent atmospheric phase modeling and correction while ensuring the correct spatial phase unwrapping of the interferogram. It is particularly suitable for GB-SAR continuous processing and quasi-real-time processing modes that require high timeliness of deformation solution.

[0067] After completing the GB-SAR image pixel thinning, phase unwrapping needs to be performed. In some embodiments, step 4 may include:

[0068] Step 41: Use the updated mask matrix to perform mask calculation on the filtered interference pattern to obtain an interference pattern sequence of discrete high-quality pixels. At this time, the phase value of each interference pattern HQP is still wrapped in the interval (-π, +π), and the wrapped phase of each interference pattern needs to be unwrapped.

[0069] Step 42: To unify the benchmark of the unwrapped phase value, a high-quality pixel in the stable region is selected as the starting point of the phase gradient integration, and the minimum-cost flow (MCF) method is used for each interference pattern to perform spatial phase unwrapping on the thinned high-quality pixel interference phase.

[0070] The present invention establishes a polynomial atmospheric phase model of the measurement area based on the unwrapped phase of high-quality pixels (HQP), and uses the model to perform atmospheric correction on the unwrapped phase of HQP. Compared with the reference point method, this method performs polynomial atmospheric correction on the basis of spatial phase unwrapping, and does not require prior information of stable areas. At the same time, models such as univariate polynomials of distance or bivariate polynomials of radar plane coordinates can be selected according to the complexity of meteorological changes.

[0071] (1) Univariate polynomial model of distance

[0072] When the measurement area is small and the meteorological changes are relatively stable, the unwrapped phase can be established. and the distance to the radar r i A univariate polynomial model of :

[0073]

[0074] In the formula, represents the atmospheric correction model, λ represents the radar wavelength, a is the model coefficient, and k is the model order.

[0075] (2) Bivariate polynomial model of radar plane coordinates

[0076] Unwrapped Phase and HQP pixel coordinates (x i ,y i )The bivariate polynomial model established:

[0077]

[0078] The polynomial model is relatively simple. In actual operation, the calculation results of multiple order models of univariate and bivariate polynomials can be compared and analyzed to preliminarily determine the optimal model and the unwrapped phase. Implement atmospheric corrections.

[0079] In order to weaken the influence of noise and gross phase errors, the present invention adopts a robust fitting method to solve the polynomial atmospheric phase model parameters. In some embodiments, the robust fitting method includes least absolute deviation (LAD) and bisquare method.

[0080] LAD is a nonparametric statistical method. Compared with the traditional least squares method, this method has less impact from outliers and is more robust. It can handle non-normal distribution data and obtain more reliable parameter estimation results without making assumptions about the data. Therefore, LAD is suitable when there are outliers or outliers in the data. The parameters of the polynomial atmospheric phase model are solved according to the principle of formula (8).

[0081]

[0082] The bisquare rule is fitted by minimizing the sum of squares of the residuals and using bisquare weights to reduce the weight of outliers. The commonly used bisquare weight function form is:

[0083]

[0084] Where x is the independent variable of the data point, x0 is the independent variable value of the fitted curve, and R is the parameter that controls the slope of the weight function.

[0085] In the bisquared weighted least squares fit, the calculation of the residuals takes into account the effect of the weight function:

[0086] E=∑w(x i )(y i -f(x i )) 2 (10)

[0087] In the formula, x i and i are the independent variables and dependent variables of the ith data point in the observed data set, f(x i ) is the fitted function for the independent variable x i By minimizing E, the optimal fitting function can be obtained.

[0088] The established polynomial atmospheric phase model is optimized, including: establishing an initial model to implement atmospheric correction, and using the atmospheric corrected unwrapped interferogram to calculate the initial linear deformation rate of the high-quality pixels, obtaining the deformation area based on the initial linear deformation rate, and removing the high-quality pixels in the deformation area. Since the existence of the deformation area will affect the accuracy of atmospheric phase modeling, resulting in deviations in the linear deformation rate calculation of the time series analysis, it is necessary to remove the HQP in the larger deformation area. Update the mask matrix and recalculate the parameters of the polynomial atmospheric phase model of the high-quality pixels in the remaining area to improve the accuracy and reliability of atmospheric modeling.

[0089] In summary, the atmospheric correction method for continuous deformation monitoring interferogram of ground-based synthetic aperture radar provided by the present invention has the following advantages and outstanding effects: the polynomial atmospheric correction method based on high-quality pixels of the present invention can improve the atmospheric correction efficiency of continuous monitoring images, effectively weaken the influence of atmospheric disturbances, and can be used for GB-InSAR deformation monitoring under complex conditions, providing reliable technical support for various safety monitoring.

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

Claims

1. A method for atmospheric correction of continuous deformation monitoring interferogram of ground-based synthetic aperture radar, characterized in that: The steps include: Step 1: Input a GB-SAR complex image sequence and form an interference image pair, and perform filtering and smoothing on each interference pattern of the interference image pair; Step 2: Based on the three parameters of thermal signal-to-noise ratio, amplitude dispersion index and coherence coefficient, a multi-threshold method is used to obtain high-quality pixels, a mask matrix is ​​calculated based on the high-quality pixels, and a mask process is performed on the filtered interference pattern to obtain the phase value of discrete high-quality pixels; Step 3: Sampling and thinning the high-quality pixels, using a uniform grid sampling method to sample and thin the high-quality pixels, including: dividing the radar image of the survey area into square grids, the side length of the divided grids is less than 1 / 3 of the width of a single complete period interference fringe, selecting the pixel with the smallest coordinate from the center of each grid as the extracted pixel, and updating the mask matrix; Step 4: performing spatial phase unwrapping on the thinned high-quality pixel interference phase; Step 5: Establish a polynomial atmospheric phase model for the measurement area, and use the polynomial atmospheric phase model to perform atmospheric correction on the unwrapped phase of the high-quality pixel; Step 6: Using a robust fitting method to solve the parameters of the polynomial atmospheric phase model; Step 7: Optimize the polynomial atmospheric phase model, use the atmospherically corrected unwrapped interferogram to calculate the initial linear deformation rate of the high-quality pixels, obtain the deformation area based on the initial linear deformation rate, remove the high-quality pixels in the deformation area, update the mask matrix, and recalculate the parameters of the polynomial atmospheric phase model of the high-quality pixels in the remaining area.

2. The atmospheric correction method for continuous deformation monitoring interferogram of ground-based synthetic aperture radar according to claim 1, characterized in that: In step 1, the first image is selected from the GB-SAR complex image sequence as the main image, and the other images are selected as slave images to form the interference image pair.

3. The atmospheric correction method for continuous deformation monitoring interferogram of ground-based synthetic aperture radar according to claim 1, characterized in that: Step 4 includes: Step 41: using the updated mask matrix to perform mask calculation on the filtered interference pattern to obtain an interference pattern sequence of discrete high-quality pixels; Step 42: Select a high-quality pixel as the starting point of the phase gradient integration, and use the minimum cost flow method to perform spatial phase unwrapping on the thinned interference phase of the high-quality pixel for each interference pattern.

4. The atmospheric correction method for continuous deformation monitoring interferogram of ground-based synthetic aperture radar according to claim 3, characterized in that: The polynomial atmospheric phase model is a univariate polynomial model of distance, or a bivariate polynomial model of radar plane coordinates.

5. The method for atmospheric correction of continuous deformation monitoring interferogram of ground-based synthetic aperture radar according to claim 4, characterized in that: The robust fitting methods include the least absolute deviation method and the double square method.

Images