Atmospheric correction methods, devices, equipment and media in high mountain valley areas

Abstract

The present invention provides an atmospheric correction method, device, equipment and medium for high mountain canyon areas, including: obtaining time-series SAR remote sensing images of high mountain canyon areas; obtaining the atmospheric mixed phase in the high mountain canyon areas according to the time-series SAR remote sensing images; and, coupling a weighted global least squares algorithm and a gross error detection algorithm to identify the vertical atmospheric delay phase and atmospheric turbulence phase in the high mountain canyon areas; removing the vertical atmospheric delay phase and atmospheric turbulence phase from the atmospheric mixed phase to obtain the atmospherically corrected deformation phase in the high mountain canyon areas. The present invention effectively solves the problems related to the difficulty in removing the atmospheric phase and extracting the deformation phase in the high mountain canyon areas, and can provide technical support for the fields of geological disaster time-series monitoring, disaster emergency response, etc.

Description

Technical Field

[0001] The present invention relates to the field of atmospheric correction technology, and in particular to an atmospheric correction method, device, equipment and medium for high mountain and canyon areas. Background Art

[0002] The application of synthetic aperture radar interferometry (InSAR) in topographic mapping and monitoring of surface deformation such as earthquake activity, volcanic movement, ground subsidence and landslide has achieved remarkable results, and these applications have also put forward higher requirements for the accuracy of interferometric measurement. Atmospheric effect is one of the main sources of error in InSAR applications and an important factor limiting the acquisition of high-precision terrain and deformation information. For repeated orbit interferometric measurement, since there is a certain time interval between the acquisition of synthetic aperture radar (SAR) images (the revisit period of Sentinel-1A is 12 days), the change of atmospheric environment will cause the atmospheric delay of SAR images to be different each time they are imaged, which will have a great impact on the interferometric phase and seriously interfere with the accuracy of surface deformation calculation. Therefore, the generation of atmospheric effect does not depend on the atmospheric environment when a single image is imaged, but depends on the difference in atmospheric environment when the interferometric image is imaged with the SAR image. The greater the difference, the greater the error caused by atmospheric effect in the interferometric phase.

[0003] For a single SAR image, atmospheric noise is mainly divided into two parts: ionospheric delay and tropospheric delay. Ionospheric delay is the signal advance effect caused by the change of free electrons along the propagation path of the radar signal. This ionospheric influence is more obvious for signals in longer bands (such as L and P bands), and has less impact on C and X bands. Tropospheric delay is the signal delay effect caused by pressure, temperature and relative humidity in the lower troposphere (<5 km), which can cause errors of up to 15 to 20 cm. Since the wavelength of the ionospheric phase is usually long, it can be removed by polynomial fitting and other methods, while the spatial scale of tropospheric delay varies greatly and is relatively difficult to separate. At present, the main research object is landslide deformation with a short wavelength of deformation signal. Therefore, the main research is on the correction of tropospheric delay of InSAR. InSAR tropospheric delay can be divided into the following three components: 1) vertical stratified atmospheric component related to terrain; 2) long-scale component caused by lateral changes in temperature, relative humidity and pressure; 3) atmospheric turbulence component.

[0004] The correlation distance of atmospheric noise can reach tens of kilometers to hundreds of kilometers. In most cases, the spatial correlation noise in an interference pattern is mainly caused by atmospheric noise. Tropospheric delay correction can be roughly divided into three categories. The first category is based on external auxiliary data correction. When the ground base station and satellite observation data are abundant, the correction effect is better, but the effect is limited in valley areas and areas with less data and will introduce additional errors; the second category is based on the characteristics of time series data. The correction is made by spatiotemporal filtering, which takes into account the randomness of tropospheric delay in the time series, but the seasonal variation of terrain-related delay phase under steep terrain cannot be solved well. Different filter window settings will greatly affect the deformation phase extraction effect; the third category is based on the empirical relationship between the phase itself and the elevation to estimate the tropospheric phase. For the study area of ​​high mountain canyons with complex terrain and steep terrain, a single simple empirical model is not suitable for tropospheric delay correction.

[0005] In view of the above problems, there is an urgent need for an optimization algorithm that can carry out time-series InSAR tropospheric delay correction in complex terrain areas such as high mountain canyons. This will provide important support for improving the time-series InSAR monitoring accuracy in complex areas such as high mountain canyons, effectively inverting deformation results such as landslides and subsidence, and enhancing SAR satellite remote sensing disaster monitoring and response in complex terrain areas. Summary of the invention

[0006] In view of this, the purpose of the present invention is to provide an atmospheric correction method, device, equipment and medium for high mountain canyon areas, which effectively solves the related problems such as the difficulty in removing the atmospheric phase and the difficulty in extracting the deformation phase in high mountain canyon areas, and can provide technical support for the fields of geological disaster time series monitoring, disaster emergency response, etc.

[0007] In a first aspect, an embodiment of the present invention provides an atmospheric correction method for a high mountain valley area, comprising:

[0008] Acquire time-series SAR remote sensing images of high mountain and canyon areas;

[0009] Obtain the atmospheric mixing phase in the high mountain valley area based on time-series SAR remote sensing images;

[0010] And, coupling the weighted overall least squares algorithm and the gross error detection algorithm to identify the vertical atmospheric delay phase and atmospheric turbulence phase in the high mountain canyon area;

[0011] The vertical atmospheric delay phase and atmospheric turbulence phase are removed from the atmospheric mixing phase to obtain the atmospherically corrected deformation phase in the high mountain canyon area.

[0012] In one embodiment, a weighted global least squares algorithm and a gross error detection algorithm are coupled to identify the vertical atmospheric delay phase and the atmospheric turbulence phase in the high mountain valley region, including:

[0013] The vertical atmospheric delay phase in the alpine canyon area is identified by using the weighted global least squares algorithm based on the digital elevation model corresponding to the alpine canyon area.

[0014] And, a gross error detection algorithm based on standardized residuals is used to identify the atmospheric turbulence phase in high mountain canyon areas.

[0015] In one embodiment, a weighted global least squares algorithm is used to identify a vertical atmospheric delay phase in a high mountain canyon region based on a digital elevation model corresponding to the high mountain canyon region, including:

[0016] Construct a mathematical empirical model between elevation data and vertical atmospheric delay phase;

[0017] The weighted global least squares algorithm is used to fit the coefficient values ​​of the mathematical empirical model to obtain the vertical atmospheric delay phase inversion model;

[0018] Through the vertical atmospheric delay phase inversion model and based on the digital elevation model corresponding to the high mountain canyon area, the vertical atmospheric delay phase in the high mountain canyon area is identified.

[0019] In one embodiment, a gross error detection algorithm based on standardized residuals is used to identify the atmospheric turbulence phase in a high mountain valley region, including:

[0020] Based on the original time series phase and imaging time of the pixels in the time series SAR remote sensing image, the gross error point set of the pixels in the time series is generated;

[0021] Constructing an indirect adjustment model of the pixel and determining the standardized residual of the indirect adjustment model to determine whether there is a gross error in the pixel based on the standardized residual;

[0022] If yes, the gross error points in the gross error point set corresponding to the pixel are identified and eliminated until there is no gross error in the pixel, and the target point set of the pixel in the time series is obtained;

[0023] Based on the target point set of pixels in time series, the atmospheric turbulence phase in the high mountain canyon area is identified.

[0024] In one implementation, determining whether a pixel has a gross error based on the standardized residual includes:

[0025] If the absolute value of the standardized residual is greater than the given threshold, it is determined that the pixel has a gross error; otherwise, if the absolute value of the standardized residual is less than or equal to the given threshold, it is determined that the pixel does not have a gross error.

[0026] In one implementation, identifying and removing the gross error points in the gross error point set corresponding to the pixel includes:

[0027] Construct the time series normal distribution vector of pixels;

[0028] Determine the normal distribution value under the pre-selected significance level. If any value in the time series normal distribution vector is less than or equal to the normal distribution value, it is determined that there is a gross error point in the gross error point set. The point corresponding to the maximum value in the time series normal distribution vector is recorded as the gross error point, and the gross error point is removed from the gross error point set.

[0029] In one embodiment, the vertical atmospheric delay phase and the atmospheric turbulence phase are removed from the atmospheric mixed phase to obtain the atmospherically corrected deformation phase in the high mountain valley area, including:

[0030] The atmospheric mixing phase is Gaussian filtered in the time domain and space domain respectively, and the vertical atmospheric delay phase and the atmospheric turbulence phase are subtracted from the atmospheric mixing phase after Gaussian filtering to obtain the atmospherically corrected deformation phase in the high mountain canyon area.

[0031] In a second aspect, an embodiment of the present invention further provides an atmospheric correction device for a mountain valley region, comprising:

[0032] Image acquisition module, used to obtain time-series SAR remote sensing images of high mountain and canyon areas;

[0033] The first phase recognition module is used to obtain the atmospheric mixing phase in the high mountain valley area according to the time-series SAR remote sensing image;

[0034] The second phase identification module is used to couple the weighted overall least squares algorithm and the gross error detection algorithm to identify the vertical atmospheric delay phase and the atmospheric turbulence phase in the high mountain canyon area;

[0035] The phase removal module is used to remove the vertical atmospheric delay phase and the atmospheric turbulence phase from the atmospheric mixed phase to obtain the atmospherically corrected deformation phase in the high mountain canyon area.

[0036] In a third aspect, an embodiment of the present invention further provides an electronic device, comprising a processor and a memory, wherein the memory stores computer-executable instructions that can be executed by the processor, and the processor executes the computer-executable instructions to implement any one of the methods provided in the first aspect.

[0037] In a fourth aspect, an embodiment of the present invention further provides a computer-readable storage medium, which stores computer-executable instructions. When the computer-executable instructions are called and executed by a processor, the computer-executable instructions prompt the processor to implement any one of the methods provided in the first aspect.

[0038] The embodiment of the present invention provides an atmospheric correction method, device, equipment and medium for a high mountain canyon area. First, a time-series SAR remote sensing image of the high mountain canyon area is obtained; then, the atmospheric mixed phase in the high mountain canyon area is obtained according to the time-series SAR remote sensing image; and, a weighted overall least squares algorithm and a gross error detection algorithm are coupled to identify the vertical atmospheric delay phase and the atmospheric turbulence phase in the high mountain canyon area; finally, the vertical atmospheric delay phase and the atmospheric turbulence phase are removed from the atmospheric mixed phase to obtain the atmospherically corrected deformation phase in the high mountain canyon area. The above method couples the weighted overall least squares algorithm and the gross error detection algorithm to identify the vertical atmospheric delay phase and the atmospheric turbulence phase in the high mountain canyon area, and then extracts the atmospherically corrected deformation phase in the high mountain canyon area, which effectively solves the problems related to the difficulty in removing the atmospheric phase and extracting the deformation phase in the high mountain canyon area, and can provide technical support for the fields of geological disaster time series monitoring and disaster emergency response.

[0039] Other features and advantages of the present invention will be described in the following description, and partly become apparent from the description, or understood by practicing the present invention. The purpose and other advantages of the present invention are realized and obtained by the structures particularly pointed out in the description, claims and drawings.

[0040] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, preferred embodiments are given below and described in detail with reference to the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] In order to more clearly illustrate the specific implementation methods of the present invention or the technical solutions in the prior art, the drawings required for use in the specific implementation methods or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some implementation methods of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0042] Figure 1 A schematic diagram of a flow chart of an atmospheric correction method for a high mountain valley area provided by an embodiment of the present invention;

[0043] Figure 2 An overall flow chart of an atmospheric correction method for high mountain canyon areas provided by an embodiment of the present invention;

[0044] Figure 3 A schematic diagram of the structure of an atmosphere correction device for high mountain and canyon areas provided by an embodiment of the present invention;

[0045] Figure 4 A schematic diagram of the structure of an electronic device provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0046] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution of the present invention will be clearly and completely described in combination with the embodiments below. Obviously, the described embodiments are part of the embodiments of the present invention, 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.

[0047] At present, the existing atmospheric correction technology cannot be applied to complex terrain areas such as high mountain canyons. Based on this, the present invention implements a method, device, equipment and medium for atmospheric correction in high mountain canyon areas, which effectively solves the problems related to the difficulty in removing the atmospheric phase and extracting the deformation phase in high mountain canyon areas, and can provide technical support for the fields of geological disaster time series monitoring, disaster emergency response, etc.

[0048] To facilitate understanding of this embodiment, firstly, a method for atmospheric correction in a high mountain valley region disclosed in an embodiment of the present invention is described in detail. Figure 1 The flow chart of an atmospheric correction method for high mountain valley areas is shown in FIG. 1 , and the method mainly includes the following steps S102 to S108:

[0049] Step S102, acquiring time-series SAR remote sensing images of high mountain valley areas.

[0050] Step S104, obtaining the atmospheric mixed phase in the high mountain canyon area according to the time-series SAR remote sensing image. Among them, the atmospheric mixed phase is a mixture of deformation phase, atmospheric delay phase and residual noise phase. In one example, the InSAR preprocessing steps of registration and image geocoding can be performed on the main image and the auxiliary image in the time-series SAR remote sensing image, and on the basis of the preprocessing, the steps of baseline combination, differential interference, phase unwrapping, flat ground and terrain phase removal are carried out in sequence to obtain the phase result of the mixed deformation phase and the atmospheric delay phase, which is the atmospheric mixed phase.

[0051] Step S106, coupling the weighted global least squares algorithm and the gross error detection algorithm to identify the vertical atmospheric delay phase and the atmospheric turbulence phase in the high mountain valley area.

[0052] In one example, a weighted global least squares algorithm is used to estimate the vertical atmospheric delay phase, and a Baarda gross error detection solution based on standardized residuals is used to estimate the atmospheric turbulence phase.

[0053] Step S108, removing the vertical atmospheric delay phase and the atmospheric turbulence phase from the atmospheric mixed phase to obtain the atmospherically corrected deformation phase in the high mountain valley area.

[0054] In one example, the atmospheric mixing phase can first be Gaussian filtered in the time domain and the space domain to eliminate the residual noise phase in the atmospheric mixing phase, and the vertical atmospheric delay phase and the atmospheric turbulence phase can be eliminated from the atmospheric mixing phase after Gaussian filtering to obtain the atmospherically corrected deformation phase in the high mountain canyon area.

[0055] The atmospheric correction method for high mountain canyon areas provided in the embodiment of the present invention identifies the vertical atmospheric delay phase and atmospheric turbulence phase in the high mountain canyon areas by coupling the weighted global least squares algorithm and the gross error detection algorithm, and then extracts the deformation phase after atmospheric correction in the high mountain canyon areas, effectively solving the related problems such as the difficulty in removing the atmospheric phase and extracting the deformation phase in the high mountain canyon areas, and can provide technical support for the fields of geological disaster time series monitoring, disaster emergency response, etc.

[0056] For ease of understanding, an embodiment of the present invention provides a specific implementation of an atmosphere correction method for a high mountain valley area.

[0057] For the aforementioned step S102, the acquired time-series SAR remote sensing images may be sentinel or other time-series SAR remote sensing images, and in addition, DEM (Digital Elevation Model) data needs to be acquired. In one example, for a key monitoring area, the time-series SAR remote sensing images and DEM data of the area are acquired to ensure the time baseline and space baseline combination requirements of the time-series SAR remote sensing images, and to obtain DEM data with higher ground resolution and more recent time updates, which, on the one hand, ensures the interferometric coherence of the time-series InSAR, facilitates geographic coding and terrain phase removal, and, on the other hand, facilitates the subsequent elevation data and atmospheric delay phase to construct an inversion model.

[0058] For the aforementioned step S104, the embodiment of the present invention provides a specific implementation method for obtaining the atmospheric mixed phase in the high mountain canyon area based on the time-series SAR remote sensing image. The purpose of this step is to obtain the mixed phase value of the deformation phase and the atmospheric delay phase for atmospheric delay correction. The InSAR preprocessing steps of the main image and the auxiliary image registration and the image geocoding in the time-series SAR remote sensing image are carried out, and on the basis of the preprocessing, the steps of baseline combination, differential interference, phase unwrapping, flat ground and terrain phase removal are carried out in sequence to obtain the phase result of the mixed deformation phase and the atmospheric delay phase, which is convenient for the subsequent correction of the atmospheric delay phase and the extraction of deformation phase information. The mixed phase model is as follows:

[0059] ;

[0060] ;

[0061] in, The atmospheric mixed phase after removing the flat bottom effect, terrain phase and noise phase, which is a mixture of deformation phase, atmospheric delay phase and residual noise phase. is the deformation phase, is the vertical atmospheric delay phase, is the atmospheric delay phase, is the atmospheric turbulence phase, is the residual noise phase, x and r are the azimuth and range directions in the radar coordinate system.

[0062] For the aforementioned step S106, the embodiment of the present invention provides an implementation method of coupling a weighted global least squares algorithm and a gross error detection algorithm to identify the vertical atmospheric delay phase and the atmospheric turbulence phase in the high mountain valley area, see the following steps 1 to 2:

[0063] Step 1: Use the weighted global least squares algorithm to identify the vertical atmospheric delay phase in the high mountain canyon area based on the digital elevation model corresponding to the high mountain canyon area.

[0064] In one example, weighted overall least squares modeling is used to estimate the vertical atmospheric delay phase. The sum of squares of the residuals is minimized, and the weight of the residuals is taken into account. Through iterative optimization, the model parameters are adjusted to minimize the residual between the elevation and the vertical atmospheric delay phase to construct an inversion model, and finally the vertical atmospheric delay phase is estimated.

[0065] In a specific implementation, refer to the following steps 1.1 to 1.3:

[0066] Step 1.1, construct a mathematical empirical model between elevation data and vertical atmospheric delay phase.

[0067] According to the relationship between tropospheric atmospheric delay and elevation, a mathematical empirical model of elevation data and atmospheric delay phase is established. Taking into account the complex terrain in mountainous canyon areas, there will be errors in the coefficient matrix in the construction of the correlation model between elevation and atmospheric phase. The general least squares method is used to solve the estimated final vertical stratified atmospheric phase, which will result in errors. Therefore, the weighted overall least squares method is considered to be used for modeling in order to better solve and estimate the vertical stratified atmospheric phase.

[0068] Among them, the vertical atmospheric delay phase expression generally has the following two modeling methods:

[0069] ;

[0070] ;

[0071] The vertical atmospheric delay phase modeling in the embodiment of the present invention couples two formulas: , is the vertically stratified atmospheric phase, is the elevation, , , , , , is the coefficient to be determined.

[0072] Step 1.2, use the weighted overall least squares algorithm to fit the coefficient values ​​of the mathematical empirical model to obtain the vertical atmospheric delay phase inversion model.

[0073] For the vertical atmosphere, since there are errors in the observed values ​​and coefficient matrix of the model, the least squares solution is not optimal but biased. Therefore, a weighted total least squares model is introduced to estimate the vertical atmosphere.

[0074] The weighted overall least squares solution model is: ;

[0075] Based on the atmospheric delay correction model solution requirements, L is the mixed phase; A is the coefficient matrix for constructing the model, which contains all the known information related to the atmospheric delay and is used to construct the observation model; e and are the error terms of the mixed phase and coefficient matrix, respectively, and X is all the parameters to be determined in the vertical atmospheric phase coupling modeling. , , , The parameter vector to be determined is formed by the combination.

[0076] The corresponding random model: ;

[0077] in is the expectation operator, representing a random vector The expected value of is zero mean, , yes After column vectorization, we get Column vector. is the variance-covariance matrix), which is used to describe the random vector uncertainty. is a known non-negative definite matrix that describes the spatial or structural correlation of the error. The cofactor matrix representing the observed data is:

[0078] ;

[0079] in, may represent the covariance matrix of local or observed data, represents the cross-covariance matrix between the observed and analyzed data, and Represent the covariance matrix and cross-covariance matrix of the analyzed data, respectively. These parameters are used to describe the uncertainty in the data and the correlation between variables.

[0080] Among them, its weight matrix is: .

[0081] Where W represents the weight matrix of local or observed data, may represent the cross-weight matrix between observed and analyzed data, and Represent the weight matrix and cross-weight matrix of the analyzed data respectively. Together, these parameters form the complete weight matrix, which is used to describe the uncertainty in the data and the correlation between variables.

[0082] In summary, by estimating the errors of the observation vector and the coefficient matrix at the same time, the EIV model can be established:

[0083] ;

[0084] in is the random error matrix of the mixed phase, is the random error matrix of the coefficient matrix; vec means that the matrix is ​​straightened by columns, from left to right; ,express The vectorized matrix; is the unit weight mean error before test; They represent the true and column vectorized co-factor matrices of the mixed phase and the coefficient matrix, respectively.

[0085] When the adjustment criteria are met: ;

[0086] The Lagrange multiplication operator-free algorithm is used for iteration. The specific steps are as follows:

[0087] (1) Initialization parameter estimation: The initial parameter estimation is calculated using the least squares method using the cofactor matrix of the observation vector L and the coefficient matrix A. , .

[0088] (2) Constructing Kronecker product: by parameter estimation The Kronecker product of ξ(i) is obtained by calculation, ;

[0089] (3) Update the co-factor matrix: perform observation error co-factor matrix The error cofactor matrix with the coefficient matrix And introduce ξ(i) in step 2 for summation, .

[0090] (4) Calculate the coefficient matrix error: by solving the updated cofactor matrix To calculate the coefficient matrix error The vectorized form of .

[0091] (5) Vectorized inverse operation of matrix error: vectorize the coefficient matrix error Convert back to matrix form by inverse operation .

[0092] (6) Update the coefficient matrix: Update to the coefficient matrix , the original coefficient matrix minus the estimated errors.

[0093] (7) Iteratively update parameter estimation: use the updated coefficient matrix and error cofactor matrix To calculate new parameter estimates , ;

[0094] (8) Check convergence: , if the difference between the new and old parameter estimates is less than the preset threshold ϵ, the iteration ends; otherwise, return to step 2 to continue the iteration.

[0095] According to the above solution model, the optimal parameter vector to be solved is obtained. , , , , and then the vertical atmospheric delay phase inversion model is obtained.

[0096] Step 1.3: Identify the vertical atmospheric delay phase in the alpine valley area based on the digital elevation model corresponding to the alpine valley area through the vertical atmospheric delay phase inversion model. , so that it can be used in the subsequent atmospheric mixing phase Remove it in.

[0097] In one example, the elevation in the digital elevation model is substituted into the vertical atmospheric delay phase inversion model, and the corresponding vertical atmospheric delay phase can be output through the vertical atmospheric delay phase inversion model.

[0098] Step 2: Identify the atmospheric turbulence phase in the high mountain valley area based on the gross error detection algorithm based on the standardized residual.

[0099] This step mainly targets the tropospheric delay effect in the atmospheric delay effect and weakens the turbulence effect. In the high mountain canyon area, the atmospheric delay response is highly correlated with the elevation and is accompanied by a certain turbulence effect. The atmospheric delay phase inversion model can be better constructed through the rough detection method. In the specific implementation, by calculating the standardized value of the residual, detecting and identifying possible outlier atmospheric turbulence phases, and correcting them through the iterative optimization process, the atmospheric turbulence phase is obtained by solving the estimation, and the random characteristics of the turbulent atmospheric delay phase and the drastic impact on the phase are eliminated.

[0100] For details, please refer to steps 2.1 to 2.3 below:

[0101] Step 2.1, based on the original time series phase and imaging time of the pixels in the time series SAR remote sensing image, generate a set of gross error points of the pixels in the time series; construct an indirect adjustment model of the pixels, and determine the standardized residual of the indirect adjustment model to determine whether the pixels have gross errors based on the standardized residuals. In one example, if the absolute value of the standardized residual is greater than a given threshold, it is determined that the pixels have gross errors; or, if the absolute value of the standardized residual is less than or equal to a given threshold, it is determined that the pixels do not have gross errors. In the case of determining that the pixels have gross errors, perform the following step 2.2.

[0102] Taking into account the differences in turbulence effects in the time domain, each pixel is normalized The original phase of the timing and imaging time , get the corresponding normalized timing phase and time . Initialize / update each pixel The set of outlier points in the time series Based on the prior deformation model , construct each pixel The indirect adjustment model is used to obtain the model residuals The prior deformation model uses a linear model to estimate the deformation residual in the simulated time series deformation sequence to construct statistics.

[0103] ;

[0104] ;

[0105] ;

[0106] in, is the observation vector, G is The associated deformation model matrix, and They are and the estimated value vector of x, and e is the observation error vector.

[0107] Standardized residual calculation: , when the standardized residual If the absolute value is greater than a given threshold, there is a gross error in the observation. is the model residual The standard deviation of .

[0108] Step 2.2, identify and remove the gross error points in the gross error point set corresponding to the pixel until there is no gross error in the pixel, and obtain the target point set of the pixel in the time series.

[0109] In one example, a time series normal distribution vector of a pixel can be constructed; the normal distribution value under a pre-selected significance level is determined; if any value in the time series normal distribution vector is less than or equal to the normal distribution value, it is determined that a gross error point exists in the gross error point set, the point corresponding to the maximum value in the time series normal distribution vector is recorded as a gross error point, and the gross error point is removed from the gross error point set.

[0110] In a specific embodiment, see (1) to (4) below:

[0111] (1) Constructing pixels The time series normal distribution vector ;

[0112] The standard deviation of the pixels is the cofactor matrix, is the identity matrix; For pixels The number of images without gross errors in the time series; s is the degree of freedom.

[0113] (2) Perform time series normal distribution vector Gross error judgment. Assuming that there is no gross error, select the significance level Conduct hypothesis testing, is the normal distribution value at this significance level. , assuming that is true, There is no gross error in, otherwise the assumption is false. There are gross errors. is the estimated deformation value.

[0114] (3) If there is no gross error, save the pixel The target point set If there is a gross error, ) are recorded as gross error points, which are removed and then iterated.

[0115] (4) Iteratively obtain the gross error point set of all pixels in the time series and use 99% confidence level for gross error detection.

[0116] Step 2.3, based on the target point set of pixels in the time series, identify the atmospheric turbulence phase in the high mountain canyon area.

[0117] The embodiment of the present invention effectively removes the two major components of the atmospheric delay correction, namely the vertical stratification phase and the turbulence phase.

[0118] For the aforementioned step S108, an embodiment of the present invention also provides an implementation method for removing the vertical atmospheric delay phase and the atmospheric turbulence phase from the atmospheric mixed phase to obtain the atmospherically corrected deformation phase in the high mountain canyon area. The atmospheric mixed phase can be Gaussian filtered in the time domain and the space domain, and the vertical atmospheric delay phase and the atmospheric turbulence phase can be subtracted from the atmospheric mixed phase after Gaussian filtering to obtain the atmospherically corrected deformation phase in the high mountain canyon area.

[0119] For details, please refer to (a) to (c) below:

[0120] (a) Gaussian filtering is performed in the time domain to decompose the phase into low-frequency and high-frequency phases. t is the input time variable of the filter. is the standard deviation of the Gaussian function, which controls the width of the filter. Here, it is the standard deviation of the Gaussian filter in the time domain.

[0121] .

[0122] in ) is the Gaussian function expression in the time domain, is the standard deviation of the Gaussian function in the time domain.

[0123] (b) Gaussian filtering is performed in the spatial domain to decompose the phase into low-frequency and high-frequency in the spatial domain. x, y are spatial coordinates. Here is the standard deviation of the Gaussian filter in the spatial domain.

[0124] .

[0125] in ) is the Gaussian function expression in the spatial domain, is the standard deviation of the Gaussian function in the spatial domain.

[0126] (c) The mixed deformation phase, vertical atmospheric phase, and atmospheric turbulence phase The identification and correction of the local turbulent atmosphere are completed. Specifically, the vertical stratified atmospheric phase and the atmospheric turbulence phase estimated by the above solution are subtracted from the atmospheric mixing phase to obtain the final deformation phase to be corrected.

[0127] In summary, the overall flow chart of the atmospheric correction method for high mountain valley areas provided by the embodiment of the present invention is as follows: Figure 2As shown, it includes: data acquisition, time series processing, atmospheric delay correction and other processes. Among them, the content of data acquisition includes time series SAR data and DEM; the steps of time series processing include primary and secondary registration / geocoding, baseline combination / differential interference / phase unwrapping / deflating, terrain phase, deformation phase mixed atmospheric phase; the step of atmospheric delay correction includes least squares inversion of vertical atmospheric phase, Baarda gross error detection of atmospheric turbulence phase, and obtaining the atmospheric corrected deformation phase. The embodiment of the present invention is applicable to a variety of SAR data at home and abroad. It can be used for time series InSAR tropospheric atmosphere correction in high mountain canyon areas. The method uses a weighted overall least squares algorithm to construct an advanced empirical model for vertical atmospheric phase solution estimation, uses the Baarda gross error detection method to settle and estimate the atmospheric turbulence phase, and couples the atmospheric correction of time series InSAR with the deformation solution of disaster areas, which can effectively solve the problems related to the difficulty of removing atmospheric phase and extracting deformation phase in high mountain canyon areas, has good regional adaptability and phase solution correction characteristics, and can provide technical support for the fields of geological disaster time series monitoring and disaster emergency response.

[0128] Based on the above-mentioned embodiments, the present invention provides an atmospheric correction device for high mountain valley areas. Figure 3 The schematic diagram of the structure of an atmospheric correction device for high mountain canyon areas is shown in FIG. The device mainly includes the following parts:

[0129] Image acquisition module 302, used to acquire time-series SAR remote sensing images of high mountain valley areas;

[0130] The first phase recognition module 304 is used to obtain the atmospheric mixing phase in the high mountain valley area according to the time-series SAR remote sensing image;

[0131] The second phase identification module 306 is used to couple the weighted global least squares algorithm and the gross error detection algorithm to identify the vertical atmospheric delay phase and the atmospheric turbulence phase in the high mountain valley area;

[0132] The phase elimination module 308 is used to eliminate the vertical atmospheric delay phase and the atmospheric turbulence phase from the atmospheric mixed phase to obtain the atmospherically corrected deformation phase in the high mountain valley area.

[0133] The atmospheric correction device for high mountain canyon areas provided by the embodiment of the present invention identifies the vertical atmospheric delay phase and atmospheric turbulence phase in the high mountain canyon areas by coupling the weighted global least squares algorithm and the gross error detection algorithm, and then extracts the deformation phase after atmospheric correction in the high mountain canyon areas, effectively solving the related problems such as the difficulty in removing the atmospheric phase and extracting the deformation phase in the high mountain canyon areas, and can provide technical support for the fields of geological disaster time series monitoring, disaster emergency response, etc.

[0134] In one implementation, the second phase identification module 306 is specifically configured to:

[0135] The vertical atmospheric delay phase in the alpine canyon area is identified by using the weighted global least squares algorithm based on the digital elevation model corresponding to the alpine canyon area.

[0136] And, a gross error detection algorithm based on standardized residuals is used to identify the atmospheric turbulence phase in high mountain canyon areas.

[0137] In one implementation, the second phase identification module 306 is specifically configured to:

[0138] Construct a mathematical empirical model between elevation data and vertical atmospheric delay phase;

[0139] The weighted global least squares algorithm is used to fit the coefficient values ​​of the mathematical empirical model to obtain the vertical atmospheric delay phase inversion model;

[0140] Through the vertical atmospheric delay phase inversion model and based on the digital elevation model corresponding to the high mountain canyon area, the vertical atmospheric delay phase in the high mountain canyon area is identified.

[0141] In one implementation, the second phase identification module 306 is specifically configured to:

[0142] Based on the original time series phase and imaging time of the pixels in the time series SAR remote sensing image, the gross error point set of the pixels in the time series is generated;

[0143] Constructing an indirect adjustment model of the pixel and determining the standardized residual of the indirect adjustment model to determine whether there is a gross error in the pixel based on the standardized residual;

[0144] If yes, the gross error points in the gross error point set corresponding to the pixel are identified and eliminated until there is no gross error in the pixel, and the target point set of the pixel in the time series is obtained;

[0145] Based on the target point set of pixels in time series, the atmospheric turbulence phase in the high mountain canyon area is identified.

[0146] In one implementation, the second phase identification module 306 is specifically configured to:

[0147] If the absolute value of the standardized residual is greater than the given threshold, it is determined that the pixel has a gross error; otherwise, if the absolute value of the standardized residual is less than or equal to the given threshold, it is determined that the pixel does not have a gross error.

[0148] In one implementation, the second phase identification module 306 is specifically configured to:

[0149] Construct the time series normal distribution vector of pixels;

[0150] Determine the normal distribution value under the pre-selected significance level. If any value in the time series normal distribution vector is less than or equal to the normal distribution value, it is determined that there is a gross error point in the gross error point set. The point corresponding to the maximum value in the time series normal distribution vector is recorded as the gross error point, and the gross error point is removed from the gross error point set.

[0151] In one implementation, the phase removal module 308 is specifically configured to:

[0152] The atmospheric mixing phase is Gaussian filtered in the time domain and space domain respectively, and the vertical atmospheric delay phase and the atmospheric turbulence phase are subtracted from the atmospheric mixing phase after Gaussian filtering to obtain the atmospherically corrected deformation phase in the high mountain canyon area.

[0153] The device provided in the embodiment of the present invention has the same implementation principle and technical effects as those of the aforementioned method embodiment. For the sake of brief description, for matters not mentioned in the device embodiment, reference may be made to the corresponding contents in the aforementioned method embodiment.

[0154] An embodiment of the present invention provides an electronic device. Specifically, the electronic device includes a processor and a storage device. The storage device stores a computer program, and when the computer program is executed by the processor, it executes the method described in any one of the above-mentioned implementation methods.

[0155] Figure 4 A structural diagram of an electronic device provided in an embodiment of the present invention, the electronic device 100 includes: a processor 40, a memory 41, a bus 42 and a communication interface 43, wherein the processor 40, the communication interface 43 and the memory 41 are connected via the bus 42; the processor 40 is used to execute an executable module stored in the memory 41, such as a computer program.

[0156] The memory 41 may include a high-speed random access memory (RAM), and may also include a non-volatile memory, such as at least one disk memory. The communication connection between the system network element and at least one other network element is realized through at least one communication interface 43 (which may be wired or wireless), and the Internet, wide area network, local area network, metropolitan area network, etc. may be used.

[0157] The bus 42 may be an ISA bus, a PCI bus, or an EISA bus, etc. The bus may be divided into an address bus, a data bus, a control bus, etc. For ease of representation, Figure 4 Only one bidirectional arrow is used in the diagram, but this does not mean that there is only one bus or only one type of bus.

[0158] Among them, the memory 41 is used to store programs, and the processor 40 executes the program after receiving the execution instruction. The method executed by the device for flow process definition disclosed in any embodiment of the above-mentioned embodiment of the present invention can be applied to the processor 40 or implemented by the processor 40.

[0159] The processor 40 may be an integrated circuit chip with signal processing capabilities. In the implementation process, each step of the above method can be completed by the hardware integrated logic circuit or software instructions in the processor 40. The above processor 40 can be a general-purpose processor, including a central processing unit (CPU), a network processor (NP), etc.; it can also be a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components. The methods, steps and logic block diagrams disclosed in the embodiments of the present invention can be implemented or executed. The general-purpose processor can be a microprocessor or the processor can also be any conventional processor. The steps of the method disclosed in the embodiment of the present invention can be directly embodied as a hardware decoding processor to be executed, or the hardware and software modules in the decoding processor can be executed. The software module can be located in a mature storage medium in the field such as a random access memory, a flash memory, a read-only memory, a programmable read-only memory, or an electrically erasable programmable memory, a register, etc. The storage medium is located in the memory 41, and the processor 40 reads the information in the memory 41 and completes the steps of the above method in combination with its hardware.

[0160] The computer program product of the readable storage medium provided in the embodiment of the present invention includes a computer-readable storage medium storing program code, and the instructions included in the program code can be used to execute the methods described in the previous method embodiments. The specific implementation can be referred to the previous method embodiments, which will not be repeated here.

[0161] If the functions are implemented in the form of software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, or the part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium, including several instructions to enable a computer device (which can be a personal computer, server, or network device, etc.) to perform all or part of the steps of the methods described in various embodiments of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk, etc., which can store program codes.

[0162] Finally, it should be noted that the above-described embodiments are only specific implementations of the present invention, which are used to illustrate the technical solutions of the present invention, rather than to limit them. The protection scope of the present invention is not limited thereto. Although the present invention is described in detail with reference to the above-described embodiments, ordinary technicians in the field should understand that any technician familiar with the technical field can still modify the technical solutions recorded in the above-described embodiments within the technical scope disclosed by the present invention, or can easily think of changes, or make equivalent replacements for some of the technical features therein; and these modifications, changes or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.

Claims

1. An atmospheric correction method for high mountain valley areas, characterized in that: include: Acquire time-series SAR remote sensing images of high mountain and canyon areas; Acquire the atmospheric mixing phase in the high mountain valley area according to the time-series SAR remote sensing image; and, coupling a weighted overall least squares algorithm and a gross error detection algorithm to identify a vertical atmospheric delay phase and an atmospheric turbulence phase in the alpine canyon region; Eliminating the vertical atmospheric delay phase and the atmospheric turbulence phase from the atmospheric mixed phase to obtain the atmospherically corrected deformation phase in the high mountain valley area; Coupling a weighted overall least squares algorithm and a gross error detection algorithm to identify the vertical atmospheric delay phase and the atmospheric turbulence phase in the alpine canyon area, including: using a weighted overall least squares algorithm, based on a digital elevation model corresponding to the alpine canyon area, to identify the vertical atmospheric delay phase in the alpine canyon area; and using a gross error detection algorithm based on standardized residuals to identify the atmospheric turbulence phase in the alpine canyon area; Using a weighted overall least squares algorithm, based on a digital elevation model corresponding to the alpine canyon area, identifying the vertical atmospheric delay phase in the alpine canyon area, including: constructing a mathematical empirical model between elevation data and the vertical atmospheric delay phase; using a weighted overall least squares algorithm to fit the coefficient value of the mathematical empirical model to obtain a vertical atmospheric delay phase inversion model; using the vertical atmospheric delay phase inversion model, based on the digital elevation model corresponding to the alpine canyon area, identifying the vertical atmospheric delay phase in the alpine canyon area; A gross error detection algorithm based on standardized residuals is used to identify the atmospheric turbulence phase in the high mountain canyon area, including: generating a gross error point set of the pixel in the time series based on the time series original phase and imaging time of the pixel in the time series SAR remote sensing image; constructing an indirect adjustment model of the pixel, and determining the standardized residual of the indirect adjustment model to determine whether the pixel has a gross error based on the standardized residual; if so, identifying and eliminating the gross error points in the gross error point set corresponding to the pixel until the pixel has no gross error, thereby obtaining a target point set of the pixel in the time series; based on the target point set of the pixel in the time series, identifying the atmospheric turbulence phase in the high mountain canyon area.

2. The atmospheric correction method for high mountain valley areas according to claim 1, characterized in that: Determining whether the pixel has a gross error based on the standardized residual includes: If the absolute value of the standardized residual is greater than a given threshold, it is determined that there is a gross error in the pixel; or, if the absolute value of the standardized residual is less than or equal to the given threshold, it is determined that there is no gross error in the pixel.

3. The atmospheric correction method for high mountain valley areas according to claim 1, characterized in that: Identifying and removing the gross error points in the gross error point set corresponding to the pixel, including: Constructing a time-series normal distribution vector of the pixel; Determine the normal distribution value under a pre-selected significance level. If any value in the time series normal distribution vector is less than or equal to the normal distribution value, determine that there is a gross error point in the gross error point set, record the point corresponding to the maximum value in the time series normal distribution vector as a gross error point, and remove the gross error point from the gross error point set.

4. The atmospheric correction method for high mountain valley areas according to claim 1, characterized in that: The vertical atmospheric delay phase and the atmospheric turbulence phase are removed from the atmospheric mixed phase to obtain the atmospherically corrected deformation phase in the high mountain valley area, including: The atmospheric mixed phase is Gaussian filtered in the time domain and the space domain respectively, and the vertical atmospheric delay phase and the atmospheric turbulence phase are subtracted from the atmospheric mixed phase after Gaussian filtering to obtain the atmospherically corrected deformation phase in the high mountain canyon area.

5. An atmospheric correction device for high mountain valley areas, characterized in that: include: Image acquisition module, used to obtain time-series SAR remote sensing images of high mountain and canyon areas; A first phase recognition module is used to obtain the atmospheric mixing phase in the high mountain valley area according to the time-series SAR remote sensing image; A second phase identification module is used to couple a weighted overall least squares algorithm and a gross error detection algorithm to identify a vertical atmospheric delay phase and an atmospheric turbulence phase in the alpine canyon area; A phase elimination module, used for eliminating the vertical atmospheric delay phase and the atmospheric turbulence phase from the atmospheric mixed phase to obtain the atmospherically corrected deformation phase in the alpine canyon area; The second phase identification module is specifically used to: identify the vertical atmospheric delay phase in the alpine canyon area by using a weighted global least squares algorithm based on the digital elevation model corresponding to the alpine canyon area; and identify the atmospheric turbulence phase in the alpine canyon area based on a gross error detection algorithm based on standardized residuals; The second phase identification module is specifically used to: construct a mathematical empirical model between elevation data and vertical atmospheric delay phase; use a weighted global least squares algorithm to fit the coefficient value of the mathematical empirical model to obtain a vertical atmospheric delay phase inversion model; through the vertical atmospheric delay phase inversion model, based on the digital elevation model corresponding to the alpine canyon area, identify the vertical atmospheric delay phase in the alpine canyon area; The second phase recognition module is specifically used to: generate a gross error point set of the pixel in the time series based on the time series original phase and imaging time of the pixel in the time series SAR remote sensing image; construct an indirect adjustment model of the pixel, and determine the standardized residual of the indirect adjustment model to determine whether the pixel has a gross error based on the standardized residual; if so, identify and eliminate the gross error points in the gross error point set corresponding to the pixel until the pixel has no gross error, thereby obtaining a target point set of the pixel in the time series; based on the target point set of the pixel in the time series, identify the atmospheric turbulence phase in the high mountain canyon area.

6. An electronic device, characterized in that: The method comprises a processor and a memory, wherein the memory stores computer executable instructions that can be executed by the processor, and the processor executes the computer executable instructions to implement the method according to any one of claims 1 to 4.

7. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores computer-executable instructions. When the computer-executable instructions are called and executed by a processor, the computer-executable instructions prompt the processor to implement the method according to any one of claims 1 to 4.

Images