A timing insar urban elevation model refinement method based on integral array combination

Abstract

The application relates to a kind of timing InSAR city elevation model refinement methods based on integral combination, comprising: S1, N+1 SAR image is registered to the same image reference coordinate system, and M difference interferograms are generated;S2, generate virtual interferogram and original interferogram merge, and the merged interferogram is divided into multiple subsets;S3, for each interferogram subset, phase is unwrapped using MCF method, and ICA method is used to estimate terrain error;S4, the estimated terrain error is converted into the terrain error phase of next subset, and the phase unwrapping and ICA terrain error estimation operation of next subset are repeated;S5, the sum of all subset estimated terrain error is added to original DEM, which is the refined elevation model.Compared with the prior art, the application is based on integral combination and ICA terrain error estimation, which can improve the interferogram unwrapping accuracy, and can suppress the influence of atmospheric delay and non-modeling error on terrain error estimation, and can efficiently and reliably refine the surface elevation model.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of synthetic aperture radar interferometry, and particularly to a time-series InSAR urban elevation model refinement method based on integer combination. BACKGROUND

[0002] Digital Elevation Model (DEM) is an important basis for smart city construction, and high-precision urban DEM can be used for urban traffic planning, flood disaster warning, and residential environment analysis. At present, urban DEM can be obtained by photogrammetry and LiDAR technology, but these technologies have limitations such as small coverage, high cost, and weather influence. As a space-borne microwave radar observation technology, InSAR has the characteristics of wide coverage, high spatial resolution, low cost, and the ability to penetrate clouds and rain, and can be used for large-scale urban 3D model reconstruction.

[0003] The principle of InSAR technology is to receive radar echo signals of the same area on the ground through satellite antenna repeated imaging, and to obtain the interferometric phase signals from two images to invert the surface deformation and terrain information. In recent years, time-series InSAR uses multiple imaging to obtain phase difference images at multiple time intervals, and estimates terrain errors based on the linear relationship between terrain error phase and vertical baseline of the interferogram to refine DEM. The existing methods for estimating terrain errors in time-series InSAR mainly include two categories: one is to generate single-master interferograms and select permanent scatterer phases, to establish a regression model with linear deformation and terrain error as parameters on the arc segment through the construction of a triangular network, and to search for and solve the terrain error parameters by spectral analysis method; the other is to generate short-baseline interferograms and select high-quality distributed scatterer phases, to establish a parameter regression model of cubic deformation or interval deformation rate and terrain error based on phase unwrapping of the interferogram, and to solve the terrain error parameters by least squares. However, since the InSAR interferometric phase observation contains deformation and terrain information, as well as non-modeled deformation, atmospheric delay, decorrelation noise, and phase ambiguity, the model regression solution accuracy is not high, and the terrain error estimation is deviated due to unwrapping error, which reduces the reliability of urban DEM refinement. In view of the shortcomings of the existing technology, the present application provides a time-series InSAR urban elevation model refinement method based on integer combination. The method first generates short-baseline interferograms by integer combination to improve the phase unwrapping accuracy, then uses independent component analysis (ICA) to separate terrain error phase and suppress the influence of atmospheric delay and non-modeled error, and finally improves the DEM refinement result through iterative estimation. SUMMARY

[0004] The present application aims at providing a time-series InSAR urban elevation model refinement method based on integral combination to overcome the defects of the prior art.

[0005] The object of the present application can be achieved by the following technical solutions.

[0006] A time-series InSAR urban elevation model refinement method based on integral combination, comprising:

[0007] S1, registering N+1 SAR images to the same image reference coordinate system, setting a space-time baseline threshold considering the influence of incoherent noise to generate M differential interferograms;

[0008] S2, generating a virtual interferogram by integral combination and merging with the original interferogram, and dividing the merged interferogram into multiple subsets according to the absolute value of the vertical baseline;

[0009] S3, for each interferogram subset, using the MCF (Minimum Cost Flow) method to unwrap the phase, and using the ICA (Independent Component Analysis) method to estimate the terrain error;

[0010] S4, converting the estimated terrain error into the terrain error phase of the next subset, and repeating steps S2 for phase unwrapping and ICA terrain error estimation operation for the next subset;

[0011] S5, adding the sum of the estimated terrain errors of all subsets to the original DEM (Digital Elevation Model) to obtain the refined DEM.

[0012] Further, the step S1 converts the external medium-low resolution DEM into the SAR image reference coordinate system through geocoding, and obtains the corresponding terrain elevation data h of the SAR image pixel through spatial interpolation, and generates M differential interferograms by differencing the initial terrain elevation data h and the interferogram.

[0013] Further, in step S2, the M differential interferogram phases are combined to generate a virtual interferogram, and the virtual interferogram phase and the corresponding interferogram vertical baseline are:

[0014]

[0015]

[0016] In the formula, is the virtual interferogram generated by combination, is the virtual interferogram vertical baseline, W{·} is a phase wrapping operator, and is the original interferogram, and is the vertical baseline corresponding to the original interferogram, a and b are combination coefficients, on this basis, the original interferogram and the virtual interferogram are combined to generate a new interferogram set, that is,

[0017]

[0018]

[0019] wherein, is the combined interferogram set, is the vertical baseline corresponding to the combined interferogram;

[0020] The combined interferogram is sequentially sorted from small to large according to the absolute value of the vertical baseline , and is divided into multiple subsets according to a preset step, that is,

[0021]

[0022]

[0023] wherein, is the combined interferogram subset, is the corresponding vertical baseline subset, and r is the number of subsets.

[0024] Further, the combination coefficients a and b control the phase noise level amplification caused by the combination of the interferogram, the value range of the combination coefficients a and b is {-2, -1, 1, 2}; at the same time, in order to avoid repeated combination caused by different a and b coefficients, for the same group of interferograms and a virtual interferogram is generated, only the virtual interferogram with the vertical baseline greater than or equal to 0 and the absolute value of the vertical baseline is the smallest is retained; at the same time, a time baseline threshold is set to delete the virtual interferogram with a long initial and tail time interval.

[0025] Further, the step S3 specifically comprises:

[0026] S31, using the MCF method to unwrap the interferogram subset ;

[0027] S32, arranging the , and using the ICA method to decompose to obtain the matrix A and the matrix S;

[0028] S33, calculating the correlation coefficient between each column of the A matrix and the vertical baseline of the interferogram.

[0029] ​​Furthermore, step S31 assumes that the phase set after unwrapping is Will Arranged so that each row represents an interferogram vector, and so on. Each column represents the phase of the same pixel in different interferograms.

[0030] Furthermore, step S32 involves arranging the... ICA decomposition yields:

[0031]

[0032] In the formula, S is the source signal matrix obtained by decomposition. Each row in S represents an independent signal source. The number of rows in the S matrix corresponds to the number of independent signal sources, and the number of columns in the S matrix corresponds to the number of interferogram elements. A is the mixing matrix of signal sources. Each row in matrix A represents the relative contribution coefficient of different independent signals. The number of rows in matrix A corresponds to the number of interferograms, and the number of columns in matrix A corresponds to the number of independent signal sources.

[0033] Further, in step S33, the correlation coefficient between each column of matrix A and the vertical baseline of the interferogram is calculated. Let the maximum absolute value of the correlation coefficient correspond to the k-th column of matrix A, and extract the k-th column vector a from matrix A. t With the k-th row vector s of matrix S t The initial terrain error is then calculated as follows:

[0034]

[0035] In the formula, λ is the radar microwave wavelength, r is the distance from the SAR satellite to the ground, and θ is the satellite incident angle.

[0036] Furthermore, step S4 will be based on The initial terrain error Δh1 obtained from phase set estimation is converted into Terrain error phase of phase set Right now:

[0037]

[0038] From again Subtracting the phase contribution of terrain error Right now:

[0039]

[0040] because Most of the phase contribution from topographic errors has been subtracted and can be unwrapped using the MCF method; meanwhile, since short vertical baselines have lower accuracy in estimating topographic errors, longer vertical baselines... There is a terrain error residual phase, and the terrain error residual amount needs to be calculated by using the ICA terrain error estimation method, so as to obtain Δh2;

[0041] And for Based on to The estimated terrain error is calculated as a terrain error phase That is:

[0042]

[0043] Subtract the terrain error phase contribution from That is:

[0044]

[0045] The MCF method is used to disentangle , and the ICA terrain error estimation method is used to calculate the terrain error residual amount, so as to obtain Δh i .

[0046] The above steps are repeated until all phase subsets in are traversed.

[0047] Further, the final terrain error obtained by summing all subset terrain errors in step S5 is:

[0048]

[0049] The final terrain error obtained is added to the original DEM to obtain a refined terrain height:

[0050] h final = h + Δh all

[0051] Finally, h final is converted into a high-precision DEM in geographic coordinates.

[0052] Compared with the prior art, the present application has the following beneficial effects:

[0053] 1. The present application is based on time series InSAR technology, and a virtual interferogram is generated by introducing an integral combination method, which increases the number of interferogram observations and improves the interferogram phase disentanglement accuracy through short baseline combination and iterative terrain error phase subtraction.

[0054] 2. The present application estimates terrain error by introducing the ICA method, suppresses the influence of non-model deformation and atmospheric delay error, and improves the accuracy of final DEM refinement through multiple iterations. BRIEF DESCRIPTION OF DRAWINGS ​

[0055] Figure 1 is a flowchart of the present application;

[0056] Figure 2 is a schematic diagram of the short baseline interferometric space-time baseline distribution of the present application;

[0057] Figure 3a is a differential interferometric Figure 1 of the present application;

[0058] Figure 3b is a differential interferometric Figure 2 of the present application;

[0059] Figure 3c is a virtual interferogram generated by combining the differential interferometric Figure 1 and the differential interferometric Figure 2 of the present application;

[0060] Figure 4a is a schematic diagram of the original DEM of the present application;

[0061] Figure 4b is a schematic diagram of the refined DEM of the present application;

[0062] Figure 5a is a differential interferometric Figure 1 corrected for topographic residuals of the present application;

[0063] Figure 5b is a differential interferometric Figure 2 corrected for topographic residuals of the present application. DETAILED DESCRIPTION

[0064] The present application will be described in detail below with reference to the accompanying drawings and specific embodiments. The embodiments are implemented on the premise of the technical solutions of the present application, and detailed implementation modes and specific operation processes are given, but the protection scope of the present application is not limited to the following embodiments.

[0065] As shown in Figure 1 is a time-series InSAR urban elevation model refinement method based on integer combination, comprising:

[0066] S1, register N+1 SAR images to the same image reference coordinate system, set the space-time baseline threshold considering the influence of decorrelation noise to generate M differential interferograms;

[0067] S2, generate a virtual interferogram by the integer combination method and merge it with the original interferogram, and divide the merged interferogram into multiple subsets according to the absolute value of the vertical baseline;

[0068] S3, for each interferogram subset, use the MCF method to unwrap the phase, and use the ICA method to estimate the terrain error;

[0069] S4, converting the estimated terrain error to the next subset of terrain error phase, repeating the phase unwrapping and ICA terrain error estimation operation of step S2 for the next subset;

[0070] S5, adding the sum of all subset estimated terrain errors to the original DEM to obtain the refined DEM.

[0071] Step S1 converts the external medium-low resolution DEM to the SAR image reference coordinate system through geocoding, and obtains the corresponding terrain elevation data h of the SAR image pixels through spatial interpolation, and generates M difference interferograms by subtracting the initial terrain elevation data h from the interferograms.

[0072] In step S2, the phases of the M difference interferograms are combined two by two to generate virtual interferograms, and the phase of the virtual interferogram and the corresponding interferogram vertical baseline are:

[0073]

[0074]

[0075] In the formula, is the combined virtual interferogram, is the virtual interferogram vertical baseline, W{·} is the phase wrapping operator, and is the original interferogram, and is the original interferogram corresponding vertical baseline, a and b are combination coefficients, on this basis, the original interferogram and the virtual interferogram are combined to generate a new set of interferograms, that is:

[0076]

[0077]

[0078] In the formula, is the combined interferogram set, is the corresponding combined interferogram vertical baseline;

[0079] According to the absolute value of the vertical baseline from small to large, the and are divided into multiple subsets, that is:

[0080]

[0081]

[0082] In the formula, is the combined interferogram subset, r is the number of subsets corresponding to the vertical baseline.

[0083] The combination coefficients a and b control the phase noise level amplification caused by the combination of the interferograms, and the combination coefficients a and b are in the range of {-2, -1, 1, 2}; at the same time, in order to avoid repeated combination caused by different a and b coefficients, for the same group of interferograms and The virtual interferogram is generated, and only the corresponding vertical baseline is greater than or equal to 0, and the absolute value of the vertical baseline is the smallest virtual interferogram; at the same time, a time baseline threshold is set to delete the virtual interferogram with a long head and tail time interval.

[0084] Step S3 specifically includes:

[0085] S31, using the MCF method to disentangle the subset of interferograms ;

[0086] S32, arranging the after disentanglement, and decomposing the arranged using the ICA method to obtain matrix A and matrix S;

[0087] S33, calculating the correlation coefficient between each column of the A matrix and the vertical baseline of the interferogram.

[0088] Step S31 assumes that the phase set after disentanglement is arranging the according to each row representing an interferogram vector, so that each column in the

[0089] represents the phase of the same pixel in different interferograms. Step S32 arranges the

[0090] after disentanglement, and decomposes the arranged using the ICA method to obtain:

[0091] In the formula, S is the source signal matrix obtained by decomposition, each row in S represents an independent signal source, the number of rows of the S matrix corresponds to the number of independent signal sources, and the number of columns of the S matrix corresponds to the number of interferogram pixels; A is a mixing matrix of the signal source, each row in the matrix A represents the relative contribution coefficient of different independent signals, the number of rows of the matrix A corresponds to the number of interferograms, and the number of columns of the matrix A corresponds to the number of independent signal sources.

[0092] Step S33 calculates the correlation coefficient between each column of the A matrix and the vertical baseline of the interferogram, assumes that the maximum absolute value of the correlation coefficient corresponds to the kth column in the matrix A, and extracts the kth column vector a t of the matrix A and the kth row vector s tThe initial terrain error is calculated as:

[0093]

[0094] where λ is the radar wavelength, r is the slant range, and θ is the satellite incidence angle.

[0095] Step S4 converts the initial terrain error Δh1estimated from the phase set to a terrain error phase The terrain error phase in the phase set That is:

[0096]

[0097] Subtract the terrain error phase contribution from That is:

[0098]

[0099] Since The majority of the terrain error contribution in the phase set has been removed, and the MCF method can be used to resolve the phase; at the same time, since the short vertical baseline estimates the terrain error with low precision, the There is a terrain error residual phase, and the ICA terrain error estimation method is used to calculate the terrain error residual, and Δh2can be obtained.

[0100] For The terrain error is estimated based on to The terrain error phase That is:

[0101]

[0102] Subtract the terrain error phase contribution from That is:

[0103]

[0104] The MCF method is used to resolve the phase , and the ICA terrain error estimation method is used to calculate the terrain error residual, and Δh i can be obtained.

[0105] The above steps are repeated until all phase subsets in are traversed.

[0106] Step S5: The final terrain error is obtained by summing the terrain errors estimated by all subsets:

[0107] ​​​

[0108] The final terrain error obtained is added to the original DEM to obtain the refined terrain elevation:

[0109] h final = h + Δh all

[0110] Finally, h final is converted into a high-precision DEM under geographic coordinates.

[0111] The experiment selects 40 TerraSAR image data in Xiqing District of Tianjin, the image resolution is 3 meters, and the time span is from March 27, 2009 to December 14, 2010. By setting the time baseline and spatial baseline threshold, 114 short baseline interference pairs are selected, and the time and space baseline distribution of the interference pairs is as shown in Figure 2 .

[0112] The 114 short baseline interference pairs are differentiated with external low-resolution DEM data (AW3D30, resolution 30 meters) to generate differential interferograms. In the processing process of the technical solution, step S2 generates a virtual interferogram by combining the interferogram set and merges the original interferogram to generate a total of 2650 interferograms. Figure 3a 、 Figure 3b and Figure 3c show an example of generating a virtual interferogram. After the merged interferogram set is processed by the ICA terrain error estimation in step S3 and the iterative refinement estimation in steps S4 and S5, the refined DEM generated is compared with the original DEM as shown in Figure 4a and Figure 4b . The original differential interferogram terrain residual phase is corrected based on the refined DEM, and the comparison before and after the interferogram correction is as shown in Figure 5a and Figure 5b .

[0113] The above describes the preferred embodiments of the application in detail. It should be understood that those skilled in the art can make many modifications and changes to the embodiments without creative labor based on the concept of the application. Therefore, any technical solution obtained by logical analysis, reasoning or limited experiment based on the existing technology according to the concept of the application should be within the protection scope determined by the claims.

Claims

1. A timing InSAR urban height model refinement method based on integral array combination, characterized in that, The application relates to a method for refining a DEM based on SAR images, and belongs to the technical field of SAR image processing. S1, registering N+1 SAR images to the same image reference coordinate system, setting a space-time baseline threshold considering the influence of incoherent noise to generate M differential interferograms; S2, generating a virtual interferogram by using an integral combination method and merging the virtual interferogram with original interferograms, and dividing the merged interferograms into multiple subsets according to the absolute value of the vertical baseline; S3, for each interferogram subset, performing phase unwrapping by using a minimum cost flow method, and estimating terrain errors by using an independent component analysis method; S4, converting the estimated terrain errors into terrain error phases of the next subset, and repeating steps S3 to perform phase unwrapping and independent component analysis terrain error estimation operations on the next subset; S5, adding the sum of all subset estimated terrain errors to the original DEM to obtain a refined height model; Step S3 specifically comprises the following steps: S31, using a minimum cost flow method on the subset of interferograms to unwrap; S32, regarding the rearranged... The matrix is ​​decomposed using the independent component analysis method. sum matrix ; S33, calculating the correlation coefficient between each column of the matrix and the subset of baselines perpendicular to the interferogram the correlation coefficient between each column of the matrix and the subset of baselines perpendicular to the interferogram Step S31 assumes the set of unwrapped phases as The unwrapped phase set is arranged in a matrix, where each column represents the phase of the same pixel in different interferograms. The unwrapped phase set is arranged in a matrix, where each column represents the phase of the same pixel in different interferograms. The unwrapped phase set is arranged in a matrix, where each column represents the phase of the same pixel in different interferograms. Step S32 arranges the Using independent component analysis for decomposition, we get: ， wherein, is the decomposed source signal matrix, each row of which represents an independent source signal, the number of rows of which corresponds to the number of independent source signals, the number of columns of which corresponds to the number of interference image elements; is the mixing matrix of the source signals, the matrix each row of which represents the relative contribution coefficient of a different independent signal, the matrix the number of rows of which corresponds to the number of interference images, the matrix the number of columns of which corresponds to the number of independent source signals. Step S33: Calculate the matrix The correlation coefficient between each column and the vertical baseline of the interferogram is set, and the matrix corresponding to the maximum absolute value of the correlation coefficient is defined. The Middle Columns, extract matrices respectively The column vector With matrix The row vector The initial terrain error is then calculated as follows: ， wherein is the radar microwave wavelength, is the SAR satellite to ground distance, is the satellite incidence angle.

2. The integer combination-based timing InSAR urban height model refinement method according to claim 1, characterized in that, In step S1, the external medium-low resolution DEM is converted into the SAR image reference coordinate system through geocoding, and the corresponding initial terrain elevation data h of the SAR image pixels is obtained through spatial interpolation; the initial terrain elevation data h is subtracted from the interferogram to generate M differential interferograms.

3. The integer combination based timing InSAR urban height model refinement method according to claim 1, characterized in that, In step S2, the M differential interferogram phases are combined two by two to generate a virtual interferogram, and the virtual interferogram phase and the corresponding interferogram vertical baseline are: ， ， wherein is a virtual interferogram generated from the combination, is a vertical baseline of the virtual interferogram, is a phase wrapping operator, and is an original interferogram, and is a vertical baseline corresponding to the original interferogram, and is a combination coefficient, on the basis of which the original interferogram and the virtual interferogram are merged to generate a new set of interferograms, namely: ， ， wherein is the merged interferogram set, is the corresponding merged interferogram vertical baseline; The merged interferograms are sorted in order of increasing absolute value of the vertical baseline The sorted interferograms are divided into a plurality of subsets, i.e.: and are divided into a plurality of subsets, i.e.: ， ， wherein is the merged interferogram subset ， is the corresponding interferogram vertical baseline subset, is the number of subsets.

4. The integer combination-based timing InSAR urban height model refinement method according to claim 3, characterized in that, The combination coefficient And The control of the phase noise level amplification caused by the combination of the interference diagram is And The value range of ; at the same time, in order to avoid the repeated combination caused by the difference of And The coefficient, for the same group of interference diagrams And Generate a virtual interference diagram, only keep the corresponding vertical baseline Greater than or equal to 0, and the absolute value of the vertical baseline The smallest virtual interference diagram; at the same time, set the time baseline threshold to delete the virtual interference diagram with long head and tail time interval.

5. The integer combination based timing InSAR urban height model refinement method according to claim 3, characterized in that, The step S4 will be based on the initial terrain error resulting from the phase set estimation converted to a terrain error phase in the phase set i.e.: ， Again from Subtracting the terrain error phase contribution i.e.: ， Since The most of phase of topographic error contribution has been removed, the unwrapping can be done by the least cost flow method; meanwhile, since the short vertical baseline estimates the topographic error with low precision, the longer vertical baseline is used to estimate the topographic error There is topographic error residual phase, which needs to use the independent component analysis topographic error estimation method to calculate the topographic error residual amount, and the topographic error residual phase can be obtained ; And for then the terrain error phase is calculated based on to the estimated terrain error i.e.: ， Again from Subtracting the terrain error phase contribution i.e.: ， The terrain error residual is calculated by using the terrain error estimation method of independent component analysis after the minimum cost flow method is used to unwrap ; and ​ The above steps are repeated until all phase subsets in the set are traversed. The above steps are repeated until all phase subsets in the set are traversed.

6. The integer combination-based timing InSAR urban height model refinement method according to claim 5, characterized in that, The final terrain error obtained by summing up the terrain errors of all subsets in step S5 is: , The final terrain error is added to the original DEM to obtain the refined terrain elevation: , h For initial terrain elevation data, finally Geocoding is converted to high-precision DEM in geographic coordinates.

Images