CN111239736A - Method, Apparatus, Equipment and Storage Medium for Surface Elevation Correction Based on Single Baseline - Google Patents

Abstract

The invention discloses a single-baseline-based surface elevation correction method, a single-baseline-based surface elevation correction device, single-baseline-based surface elevation correction equipment and a single-baseline-based storage medium, wherein the method comprises the following steps: acquiring two full-aperture images of a research area, and acquiring each full-aperture image by adopting a sub-aperture decomposition technology to decompose the full-aperture image into two sub-view images; preprocessing the full-aperture image and the sub-view image to obtain a corresponding interference pattern; performing terrain removing processing on the two sub-view interferograms by using DEM data, and removing flat ground and orbit error phases to obtain two sub-view differential interferograms comprising the same troposphere delay error; performing correlation analysis on high-frequency wavelet coefficients of the two sub-view differential interferograms by adopting a wavelet decomposition technology, and extracting troposphere delay error phases; and removing the obtained troposphere delay error phase from the full-aperture interferogram, and calculating according to the satellite orbit parameters to obtain earth surface elevation data. The invention can extract the troposphere delay error phase only by using two full-aperture images and obtain high-precision earth surface elevation data.

Description

Method, Apparatus, Equipment and Storage Medium for Surface Elevation Correction Based on Single Baseline

technical field

The invention relates to the application of synthetic aperture radar interferometry (InSAR) technology in large-scale terrain measurement, belongs to the technical field of surface elevation extraction by spaceborne synthetic aperture radar, and particularly relates to a single baseline-based surface elevation correction method, device and equipment and storage media.

Background technique

Digital elevation model (DEM) is an indispensable basic drawing for national economic construction, social development and national security. Synthetic Aperture Radar (InSAR) technology has the advantages of all-weather, all-day, short observation period and wide range, and has become a powerful tool for obtaining global DEM. However, during the InSAR earth observation, the radar signal will be affected by many factors, which also seriously affects the accuracy of the terrain measurement using InSAR technology. For the currently commonly used heavy-orbit spaceborne SAR systems, the tropospheric delay error is one of the most important error sources affecting the measurement accuracy. Therefore, in order to obtain a large-scale and high-precision InSAR digital elevation model, it is necessary to effectively estimate and remove the tropospheric delay error contained in the interferometric phase.

At present, the tropospheric delay error correction methods commonly used in InSAR measurements mainly include the following three categories: 1) phase-elevation correlation analysis method; 2) external water vapor data-assisted method; 3) time-series InSAR analysis method. Among them, the first type of method usually needs to rely on certain prior information, that is, the tropospheric delay error phase has a certain correlation with the ground height, and this method can only remove the tropospheric delay error phase related to the terrain. When terrain-independent turbulence errors dominate, this method cannot obtain better correction results; compared with the first type of method, the second type of method requires better spatiotemporal consistency between SAR data and external water vapor data, which This limits the popularization of this method. Different from the first two types of methods, in order to use the third type of method to obtain better correction results, a large amount of SAR data covering the same area needs to be covered. However, in practical applications, due to the lack of external data or the small amount of SAR data, it is difficult for the above methods to achieve good tropospheric delay error correction.

Therefore, it is necessary to design a tropospheric delay error correction method without relying on external water vapor assistance data or excess SAR data.

SUMMARY OF THE INVENTION

The technical problem to be solved by the present invention is to provide a single-baseline-based surface elevation correction method, device, equipment and storage medium, which can calibrate the surface phase by extracting the tropospheric delay error phase, thereby obtaining higher-precision surface elevation data .

For realizing the above-mentioned technical purpose, the present invention adopts following technical scheme:

A method for surface elevation correction based on a single baseline, comprising the following steps:

Step 1: Acquire two full-aperture SAR images of the study area, and use sub-aperture decomposition technology to decompose each full-aperture SAR image into two sub-view images corresponding to azimuth incident angles;

Step 2, preprocessing two full-aperture SAR images to obtain a full-aperture interferogram; preprocessing two sub-view images with the same azimuth incident angle to obtain a sub-view interferogram corresponding to the azimuth incident angle;

Step 3: For the two sub-parallax interferograms, both use DEM data and use differential interferometry to perform topographic processing, and then remove the flat-ground phase and orbit error phase processing to obtain two sub-parallax differential interferograms as shown in the following formulas:

where α ₁ and α ₂ represent two different azimuthal incident angles, φ _diff is the differential interference phase in the subparallax interferogram, Δφ _topo is the topographic error phase, φ _atmo is the tropospheric delay error phase, and φ _noise is the random noise error phase;

Step 4: Perform wavelet decomposition on the two sub-parallax differential interferograms, perform correlation analysis on the high-frequency wavelet coefficients of the two sub-parallax differential interferograms at each decomposition scale, and obtain the difference between the two high-frequency wavelet coefficients at the corresponding decomposition scale. The correlation value of , and the wavelet coefficient related to the tropospheric delay error phase is calculated according to the correlation value and any high-frequency wavelet coefficient; then the wavelet coefficient related to the tropospheric delay error phase is subjected to wavelet inverse transformation to obtain the tropospheric delay error phase. ;

In step 5, the full aperture interferogram is unwrapped, and the tropospheric delay error phase obtained in step 4 is removed from the unwrapped phase of the full aperture interferogram; finally, according to the satellite orbit parameters and the full aperture interferogram that removes the tropospheric delay error phase Unwrapped phase, calculated to obtain surface elevation data.

In a better technical solution, the two sub-view interferograms are de-topographically processed using two DEM data obtained by different techniques.

In a better technical solution, the wavelet decomposition of the sub-parallax interferogram is expressed as:

In the formula, m and n represent the number of rows and columns of the interferogram, respectively, x, _y are the coordinates of the pixel in the radar coordinate system, i _x and _{i y} represent the ith row and the ith column of the interferogram, respectively, Φ and Ψ represent the translation function and mother wavelet function of wavelet decomposition, respectively, J represents the optimal decomposition scale, j represents different decomposition scales, v, w represent low-frequency wavelet coefficients and high-frequency wavelet coefficients, ε=1, 2, 3 Corresponding to the horizontal, vertical and diagonal directions respectively;

和

按以下公式进行相关值计算：The correlation analysis method in step 4 is to take the high-frequency wavelet coefficients of the two sub-parallax differential interferograms at each decomposition scale.

and

The correlation value is calculated according to the following formula:

和

时，代入上述相关性计算公式中，即可通过方程组计算得到两个子视差分干涉图在分解尺度j时的相关值f_j和整体偏移系数c_j；High-frequency wavelet coefficients obtained when ε = 1, 2, 3 are obtained from wavelet decomposition

and

When , substituting into the above correlation calculation formula, the correlation value f _j and the overall offset coefficient c _j of the two sub-parallax differential interferograms at the decomposition scale j can be obtained through the equation system calculation;

Then use the correlation f _j of the two sub-parallax differential interferograms when decomposing scale j and the overall offset coefficient c _j to calculate the high-frequency wavelet coefficient corresponding to the tropospheric delay error phase according to the following formula:

代表对流层延迟误差相位所对应的高频小波系数；

表示方位向入射角为α_i子视差分干涉图在分解尺度j时高频小波系数。In the formula,

represents the high-frequency wavelet coefficient corresponding to the tropospheric delay error phase;

Represents the high-frequency wavelet coefficients when the azimuthal incident angle is α _i sub-parallax interferogram at decomposition scale j.

In a better technical solution, when the flat-earth phase is removed, the calculation formula of the flat-earth phase is:

In the formula, φ _flat is the phase of the flat ground, λ is the wavelength of the radar signal, and B _// is the length of the parallel baseline.

In a better technical solution, the method of removing the orbital error phase is as follows: use a polynomial model integrated with the surface elevation h to fit the orbital error phase, solve the unknown parameters in the polynomial model by the least square method, and then convert the orbital error phase from the sub- The unwrapped phase of the corresponding pixel in the parallax interferogram is removed to obtain a new sub-parallax interferogram.

Finally, the orbit error phase φ _orbit is removed from the unwrapped phase of the corresponding pixel in the subparallax interferogram.

In a better technical solution, in step 1, the sub-aperture decomposition technology is used to decompose the full-aperture SAR image as follows:

First, perform one-dimensional Fourier transform on the full-aperture SAR image to obtain the spectrum of the SAR image in the azimuth direction;

Then, the spectrum of the SAR image in the azimuth direction is divided according to different azimuth incident angles, and the spectrum corresponding to the azimuth incident angle is obtained;

Finally, inverse Fourier transform is performed on all the spectrums corresponding to the azimuthal incident angle to obtain an image corresponding to the azimuthal incident angle, that is, the subview image corresponding to the azimuthal incident angle.

In a better technical solution, the preprocessing of the full-aperture SAR image includes: co-registering two full-aperture SAR images, and then pixel-by-pixel conjugate multiplication to generate a full-aperture interferogram;

Preprocessing the two sub-view images corresponding to each azimuth incident angle of the main and auxiliary images includes: multiplying the two sub-view images corresponding to the current azimuth incident angle pixel by pixel to obtain the sub-view corresponding to the current azimuth incident angle. Interferogram.

The present invention also provides a surface elevation correction device based on a single baseline, comprising:

The sub-aperture decomposition module is used to: acquire two full-aperture SAR images, and use the sub-aperture decomposition technology to decompose each full-aperture SAR image into two sub-view images corresponding to the azimuth incident angle;

The data preprocessing module is used to: preprocess the two full-aperture SAR images to obtain the full-aperture interferogram; preprocess the two sub-view images with the same azimuth incident angle to obtain the sub-view interference corresponding to the azimuth incident angle picture;

The differential interferometric phase modeling module is used to: use DEM data for two sub-view interferograms and use differential interferometric technology to perform topographic processing, and then remove the flat phase and orbit error phase processing, and obtain two as shown in the following formulas The subparallax interferogram of :

where α ₁ and α ₂ represent two different azimuthal incident angles, respectively, φ _diff is the differential interference phase in the subparallax interferogram, Δφ _topo is the topographic error phase, φ _atmo is the tropospheric delay error phase, φ _noise is the random noise error phase;

The tropospheric delay error phase extraction module is used to: perform wavelet decomposition on two sub-parallax differential interferograms, perform correlation analysis on the high-frequency wavelet coefficients of the two sub-parallax differential interferograms at each decomposition scale, and obtain two sub-parallax differential interferograms at the corresponding decomposition scale. The correlation value between the high-frequency wavelet coefficients is calculated according to the correlation value and any high-frequency wavelet coefficient to obtain the wavelet coefficient related to the phase of the tropospheric delay error; The tropospheric delay error phase can be obtained;

The surface elevation data correction module is used for: unwrapping the full-aperture interferogram, and removing the obtained tropospheric delay error phase from the unwrapping phase of the full-aperture interferogram; finally, according to the satellite orbit parameters and removing the tropospheric delay error phase The full-aperture interferogram wraps the phase and calculates the surface elevation data.

The present invention also provides a device, comprising a processor and a memory; wherein: the memory is used for storing computer instructions; the processor is used for executing the computer instructions stored in the memory, specifically performing any of the above methods.

The present invention also provides a computer storage medium for storing a program, and when the program is executed, it is used to implement any one of the above-mentioned methods.

beneficial effect

A tropospheric delay error correction method based on sub-aperture correlation analysis of the present invention has the following advantages:

(1) The present invention obtains sub-view interferograms under different azimuth incident angles based on the sub-aperture decomposition technology, extracts the same tropospheric delay error phase through correlation analysis and eliminates it, and weakens the influence of the tropospheric delay error on InSAR height measurement, thereby Improve the measurement accuracy of InSAR technology;

(2) The present invention does not need too much SAR data, only needs two full-aperture SAR images to form a single baseline, and without external water vapor auxiliary data, the tropospheric delay error phase is independent of the azimuth incident angle , the distribution characteristics of the tropospheric delay phase can be obtained;

(3) The present invention uses correlation analysis to extract the same tropospheric delay error phase components from two sub-parallax interferograms, so as to extract the tropospheric delay error phase that eliminates aliasing in the interferograms. Compared with the traditional method, the present invention has The vertical stratified atmosphere related to the terrain and the turbulent atmospheric delay component not related to the ground height can be removed at the same time, thus greatly broadening the application range of InSAR altimetry.

Description of drawings

FIG. 1 is a flowchart of a method according to an embodiment of the present invention.

FIG. 2 is a detailed flowchart of the sub-aperture correlation analysis tropospheric delay error correction algorithm according to the embodiment of the present invention;

FIG. 3 is a schematic diagram of sub-aperture decomposition of a SAR image according to an embodiment of the present invention.

FIG. 4 is a DEM result obtained by the method described in the embodiment of the present invention.

Detailed ways

In order to better illustrate the method and steps of the present invention, the present invention is further described in detail by using the ALOS1PALSAR radar data of two scenes located in the Los Angeles area with an interval of 92 days. It should be understood that the specific implementations described herein are only used to explain the present invention, but not to limit the present invention.

An embodiment of the present invention provides a single-baseline-based surface elevation correction method, as shown in Figures 1 and 2, including the following steps:

Step 1: Acquire two full-aperture SAR images in the study area, and use sub-aperture decomposition technology to decompose each full-aperture SAR image into two sub-view images corresponding to the azimuth incident angle. The basic principle is shown in Figure 3; Two full-aperture SAR images correspond to two sub-view images, which are called forward-view images and rear-view images.

Among them, two full-aperture SAR images: one full-aperture SAR image is the ALOS1PALSAR radar data acquired on July 11, 2010, and this image is the main image in this embodiment; the other full-aperture SAR image is from October 2010 The ALOS1 PALSAR radar data acquired on January 11, this example uses this image as the auxiliary image.

The process of decomposing full-aperture SAR images using sub-aperture decomposition technology is as follows:

First, perform one-dimensional Fourier transform on the full-aperture SAR image to obtain the spectrum of the SAR image in the azimuth direction;

Then, the spectrum of the SAR image in the azimuth direction is divided according to different azimuth incident angles, and the spectrum corresponding to the azimuth incident angle is obtained;

Finally, inverse Fourier transform is performed on all the spectrums corresponding to the azimuthal incident angle to obtain an image corresponding to the azimuthal incident angle, that is, the subview image corresponding to the azimuthal incident angle.

In the specific implementation, the sub-aperture decomposition technology is used to decompose each full-aperture SAR image into two sub-view images with different azimuth incident angles. Specifically, by using two band-pass filters with different azimuth incident angles, different A subview image of the azimuthal viewing angle.

Step 2, preprocessing two full-aperture SAR images to obtain a full-aperture interferogram; preprocessing two sub-view images with the same azimuth incident angle to obtain a sub-view interferogram corresponding to the azimuth incident angle;

The preprocessing of full-aperture SAR images includes: co-registration of two full-aperture SAR images, and then pixel-by-pixel conjugate multiplication to generate a full-aperture interferogram;

Preprocessing the two sub-view images of each azimuth incident angle includes: performing interference processing on the two sub-view images corresponding to the current azimuth incident angle to obtain a sub-view interferogram corresponding to the current azimuth incident angle.

Step 3: For the two sub-parallax interferograms, both use DEM data and use differential interferometry to perform topographic processing, and then remove the flat-ground phase and orbit error phase processing to obtain two sub-parallax differential interferograms as shown in the following formulas:

where α ₁ and α ₂ represent two different azimuthal incident angles, φ _diff is the differential interference phase in the subparallax interferogram, Δφ _topo is the topographic error phase, φ _atmo is the tropospheric delay error phase, and φ _noise is the random noise error phase.

Specifically, the space-borne repetitive orbit InSAR will be affected by various factors such as atmosphere and surface deformation during earth observation. Therefore, the InSAR interference phase appears as the superposition of various phase components.

In order to ensure the phase quality of the interferogram, this embodiment selects two SAR images with short time baselines, and the surface deformation is almost zero, so the influence of the surface deformation can be ignored; in addition, compared with the tropospheric delay, the ionospheric delay The magnitude of is generally relatively small, so the effect of ionospheric delay is ignored. Therefore, for the InSAR terrain measurement in the embodiment of the present invention, without considering the influence of surface deformation and ionospheric delay, the phase φ _int of the interferogram under different observation angles can be expressed as:

where:

α ₁ and α ₂ represent two different azimuthal incident angles;

φ _topo is the topographic phase related to topographic measurement, which is related to the bare surface elevation h and the height Δh of the surface scattering target;

可以根据卫星的轨道参数将其进行有效的去除；λ是雷达波长，B_//是平行基线；φ _flat is the flat-earth phase, which can be expressed as:

It can be effectively removed according to the orbital parameters of the satellite; λ is the radar wavelength, B _// is the parallel baseline;

φ _orbit is the orbit error phase, which is caused by inaccurate orbit parameters. In conventional InSAR processing, polynomial model fitting method is usually used to remove it;

φ _noise is the phase of random noise error, which can be weakened by interferogram filtering;

主要产生原因为两次成像过程中微波信号受到不同水汽环境的干扰，干涉后以相位的形式体现在干涉图中。φ _atmo is the tropospheric delay error phase, which can be expressed as:

The main reason is that the microwave signal is interfered by different water vapor environments during the two imaging processes, and is reflected in the interferogram in the form of phase after the interference.

According to the above analysis of the phase components in the interferogram, only the tropospheric delay error in the interferogram is difficult to be effectively removed by a fixed model, which has become the main source of error in the topographic measurement of spaceborne heavy-orbit SAR.

Studies have shown that the change of the azimuthal incident angle will lead to the change of the backscattering characteristics, so different sub-view interferograms have different interference phase centers. In addition, since the sub-view images under different azimuth incident angles share approximately the same orbital parameters, from the analysis of the above-mentioned flat-ground phase and orbit error phase, it can be obtained that each sub-view interference phase contains the same flat-ground phase and orbit error phase. . Since the synthetic aperture time of SAR sensors is generally very short (usually several seconds), and considering the small beam width angle of spaceborne SAR sensors, it can be considered that the atmospheric components that different sub-aperture radar waves pass through during transmission are approximately the same, As a result, different sub-aperture images contain the same tropospheric delay error, and these same tropospheric delay errors are converted into the same tropospheric delay error phase after interference processing. Although each aperture image contains only a portion of the full-aperture SAR image spectrum, it contains the same tropospheric delay error phase. This feature makes it possible to extract and remove the same tropospheric delay errors.

In this step 3, in order to extract the same tropospheric delay error phase, it is necessary to first remove the terrain phase, flat phase and orbit error phase components corresponding to the bare ground elevation in each sub-view interferogram.

First, the differential interference technology is used to remove the terrain, that is, the terrain phase is simulated by using external DEM data, and the terrain phase is removed from the interference phase:

为采用外部DEM数据h模拟得到的地形相位，△φ_topo为地形误差相位，由外部DEM数据的不准确所导致；B_⊥代表垂直基线的长度，λ为雷达信号波长，R₁和R₂分别表示主、辅影像的天线中心到地面目标点的距离。In the formula,

is the terrain phase simulated by the external DEM data h, Δφ _topo is the terrain error phase, which is caused by the inaccuracy of the external DEM data; B _⊥ represents the length of the vertical baseline, λ is the radar signal wavelength, R ₁ and R ₂ respectively Indicates the distance from the antenna center of the main and auxiliary images to the ground target point.

Although the interferograms at different azimuth observation angles correspond to different phase centers, in the process of removing the topography by differential interferometry, if the same external DEM data is used to remove the topographic phases in the interferograms of different polarizations, it will lead to various sub-interferences. The subparallax interferograms corresponding to the visual images contain similar terrain error phases, which will seriously interfere with the subsequent extraction of the same tropospheric delay error phase. Therefore, in this embodiment, the external DEM data obtained by different technical means are used to perform terrain removal processing, and the terrain phases in the corresponding sub-view interferograms are respectively removed, which can enhance the discrimination of the two finally obtained sub-parallax differential interferograms. Improve estimation accuracy of tropospheric delay error.

根据卫星的轨道参数λ、R₂、R₁计算平地相位，进而将平地相位从子视干涉相位中去除：Then, use the flat-earth phase calculation formula

Calculate the flat-earth phase according to the satellite's orbital parameters λ, R ₂ , and R ₁ , and then remove the flat-earth phase from the sub-view interference phase:

Then, use the polynomial model integrated into the surface elevation h to fit the orbit error phase, solve the unknown parameters in the polynomial by the least square method, and then remove the orbit error phase from the unwrapped phase of the corresponding pixel in the subparallax interferogram. , to obtain a new subparallax interferogram.

其中x,y是像素点在雷达坐标系的坐标，h_i为第i个像素对应的地表高程，a₁,a₂,a₃为待求参数；Specifically in this embodiment, the orbit error phase φ _orbit is fitted by using the polynomial a ₀ +a ₁ x+a ₂ y+a ₃ h, and the orbit error phase model is obtained as:

where x, y are the coordinates of the pixel in the radar coordinate system, hi is the surface elevation corresponding to the _ith pixel, and a ₁ , a ₂ , and a ₃ are the parameters to be determined;

In order to accurately estimate the values of the unknown parameters a ₁ , a ₂ , a ₃ , the above orbital error phase model can be used to fit the unwrapped subparallax interferometric phase, and then the least squares method can be used to solve the unknown parameters a ₁ , a ₂ , a ₃ ; then, according to the coordinates x, y and the surface elevation hi corresponding to each _pixel , the orbit error phase φ _orbit can be calculated through the above orbit error phase model.

It should also be noted that: for the nonlinear track error phase, in practical applications, a high-order polynomial can be selected to replace the polynomial a ₀ +a ₁ x+a ₂ y+a ₃ h in this embodiment for calculation.

The orbit error phase φ _orbit is removed from the unwrapping phase of the corresponding pixel in the sub-parallax interferogram, and a new sub-parallax interferogram can be obtained:

Among them, φ _diff is the differential interference phase in the sub-parallax differential interferogram.

Step 4: Perform wavelet decomposition on the two sub-parallax differential interferograms, perform correlation analysis on the high-frequency wavelet coefficients of the two sub-parallax differential interferograms at each decomposition scale, and obtain the difference between the two high-frequency wavelet coefficients at the corresponding decomposition scale. The correlation value of , and the wavelet coefficient related to the tropospheric delay error phase is calculated according to the correlation value and any high-frequency wavelet coefficient; then the wavelet coefficient related to the tropospheric delay error phase is subjected to wavelet inverse transformation to obtain the tropospheric delay error phase. .

The complexity of the spatial distribution of tropospheric delay errors makes it difficult to separate them effectively by simple spatial filtering. However, different phase components in the differential interferogram appear as different wavelength components in the frequency domain, so wavelet transform can decompose the differential interferogram into different wavelength components, allowing us to study them at different scales. Therefore, the embodiments of the present invention use wavelet transform to decompose the differential interferogram into different scales, and estimate the magnitude of the tropospheric delay error phase in each decomposition scale.

The wavelet decomposition of the subparallax interferogram is expressed as:

In the formula, m and n represent the number of rows and columns of the interferogram, respectively, x, _y are the coordinates of the pixel in the radar coordinate system, i _x and _{i y} represent the ith row and the ith column of the interferogram, respectively, Φ and Ψ represent the translation function and mother wavelet function of wavelet decomposition, respectively, J represents the optimal decomposition scale, j represents different decomposition scales, v, w represent low-frequency wavelet coefficients and high-frequency wavelet coefficients, ε=1, 2, 3 Corresponding to the horizontal, vertical and diagonal directions respectively;

)分别进行如下的整体相关性分析计算：In the two subparallax interferograms, the phase of the tropospheric delay error is the same as analyzed above, and the tropospheric delay error is usually expressed as long-wavelength high-frequency information in the frequency domain, so their high-frequency wavelet coefficients have a certain similarity, then For the high-frequency coefficients of different sub-view interferograms at each decomposition scale (represented as

) respectively perform the following overall correlation analysis and calculation:

和

时，代入上述相关性计算公式中，即可通过方程组计算得到两个子视差分干涉图在分解尺度j时的相关值f_j和整体偏移系数c_j；High-frequency wavelet coefficients obtained when ε = 1, 2, 3 are obtained from wavelet decomposition

and

When , substituting into the above correlation calculation formula, the correlation value f _j and the overall offset coefficient c _j of the two sub-parallax differential interferograms at the decomposition scale j can be obtained through the equation system calculation;

Then use the correlation f _j of the two sub-parallax differential interferograms when decomposing scale j and the overall offset coefficient c _j to calculate the high-frequency wavelet coefficient corresponding to the tropospheric delay error phase according to the following formula:

代表对流层延迟误差相位所对应的高频小波系数；

表示方位向入射角为α_i子视差分干涉图的在分解尺度j时高频小波系数。对流层延迟误差相位对应的高频小波系数时，具体可取两个方位向入射角中任意一个对应的子视差分干涉图得到的高频小波系数

进行计算。In the formula,

represents the high-frequency wavelet coefficient corresponding to the tropospheric delay error phase;

Represents the high-frequency wavelet coefficients at the decomposition scale j of the subparallax interferogram whose azimuthal incident angle is α _i . When the high-frequency wavelet coefficient corresponding to the tropospheric delay error phase, the high-frequency wavelet coefficient obtained from the sub-parallax interferogram corresponding to any one of the two azimuth incident angles can be used.

Calculation.

Finally, in order to extract the tropospheric delay error phase contained in the subparallax interferogram, the wavelet coefficients related to the tropospheric delay error phase obtained above are subjected to wavelet inverse transformation, and the tropospheric delay error phase contained in the phase of the full aperture interferogram can be obtained. .

In step 5, the full aperture interferogram is unwrapped, and the tropospheric delay error phase obtained in step 4 is removed from the unwrapped phase of the full aperture interferogram; finally, according to the satellite orbit parameters and the full aperture interferogram that removes the tropospheric delay error phase To plot the entanglement phase, the surface elevation data h can be calculated by the following formula:

In the formula, θ is the antenna incident angle of each pixel point of the main image, α is the baseline inclination angle, and φ' _topo represents the full-aperture interferogram wrapping phase with the tropospheric delay error phase removed.

In a more preferred embodiment, it also includes performing noise filtering processing on the full-aperture interferogram wrapping phase φ' _topo , that is, removing the random noise error phase, and then extracting the random noise error phase and the tropospheric delay error phase from the full-aperture interferogram. Remove it from the unwrapping phase, record the terrain phase at this time as φ” _topo , and then use the formula to calculate the surface elevation data with higher precision:

As shown in Figure 4, Figure 4(a)-(c) respectively represent the surface elevation data InSAR DEM1 obtained by directly using the full-aperture SAR image processing (without tropospheric delay error correction), and using the embodiment of the present invention to remove the tropospheric delay error The surface elevation data InSAR DEM2 and the surface elevation data SRTM DEM provided by USGS after calculating the phase; then, the SRTM DEM is used as the reference data to carry out qualitative analysis on InSAR DEM1 and InSAR DEM2 respectively, and the elevation errors of InSAR DEM1 and InSAR DEM2 are obtained respectively. As shown in Figures 4(d) and 4(e) respectively, it can be seen that the implementation method of the present invention can effectively remove the tropospheric delay error in the InSAR topographic measurement.

The present invention also provides a single-baseline-based surface elevation correction device, corresponding to the above method embodiments, including:

The sub-aperture decomposition module is used to: acquire two full-aperture SAR images, and use the sub-aperture decomposition technology to decompose each full-aperture SAR image into two sub-view images corresponding to the azimuth incident angle;

The data preprocessing module is used to: preprocess the two full-aperture SAR images to obtain the full-aperture interferogram; preprocess the two sub-view images with the same azimuth incident angle to obtain the sub-view interference corresponding to the azimuth incident angle picture;

The differential interferometric phase modeling module is used to: use DEM data for two sub-view interferograms and use differential interferometric technology to perform topographic processing, and then remove the flat phase and orbit error phase processing, and obtain two as shown in the following formulas The subparallax interferogram of :

where α ₁ and α ₂ represent two different azimuthal incident angles, respectively, φ _diff is the differential interference phase in the subparallax interferogram, Δφ _topo is the topographic error phase, φ _atmo is the tropospheric delay error phase, φ _noise is the random noise error phase;

The tropospheric delay error phase extraction module is used to: perform wavelet decomposition on the two sub-parallax differential interferograms, perform correlation analysis on the high-frequency wavelet coefficients of the two sub-parallax differential interferograms at each decomposition scale, and take the correlation greater than the preset value. The high-frequency wavelet coefficients of the threshold are all used as the wavelet coefficients related to the phase of the tropospheric delay error; then the wavelet coefficients related to the phase of the tropospheric delay error are subjected to wavelet inverse transformation to obtain the tropospheric delay error phase;

The surface elevation data correction module is used to: remove the obtained tropospheric delay error phase from the full aperture interferogram; finally, calculate the surface elevation data according to the satellite orbit parameters and the full aperture interferogram from which the tropospheric delay error phase is removed.

The present invention further provides a device including a processor and a memory; wherein: the memory is used to store computer instructions; the processor is used to execute the computer instructions stored in the memory, and specifically perform the steps described in the above method embodiments.

The present invention also provides a computer storage medium for storing a program, and when the program is executed, it is used to implement the method described in the above method embodiments.

The above embodiments are the preferred embodiments of the application, and those of ordinary skill in the art can also carry out various transformations or improvements on this basis. Without departing from the general concept of the application, these transformations or improvements should belong to the present application. within the scope of the application for protection.

Claims (10)

1. A single-baseline-based surface elevation correction method is characterized by comprising the following steps:

step 1, acquiring two full-aperture SAR images of a research area, and decomposing each full-aperture SAR image into two sub-view images corresponding to azimuth incident angles by adopting a sub-aperture decomposition technology;

step 2, preprocessing the two full-aperture SAR images to obtain a full-aperture interferogram; preprocessing two sub-view images with the same azimuth incident angle to obtain a sub-view interferogram corresponding to the azimuth incident angle;

step 3, for the two sub-view interferograms, DEM data are used, a differential interference technology is adopted for terrain removing processing, then flat ground phase and orbit error phase removing processing are carried out, and two sub-view interferograms shown in the following formulas are obtained:

in the formula, α₁And α₂Representing two different azimuthal angles of incidence, phi_diffIs the differential interference phase in the sub-parallax partial interference pattern, △ phi_topoFor the phase of the terrain error, phi_atmoFor tropospheric delay error phase, phi_noiseIs a random noise error phase;

step 4, performing wavelet decomposition on the two sub-view differential interferograms, performing correlation analysis on the high-frequency wavelet coefficients of the two sub-view differential interferograms under each decomposition scale to obtain correlation values between the two high-frequency wavelet coefficients under the corresponding decomposition scale, and calculating to obtain wavelet coefficients related to troposphere delay error phases according to the correlation values and any one high-frequency wavelet coefficient; then, performing inverse wavelet transform on the wavelet coefficient related to the troposphere delay error phase to obtain a troposphere delay error phase;

step 5, carrying out unwrapping processing on the full-aperture interferogram, and removing the troposphere delay error phase obtained in the step 4 from the unwrapping phase of the full-aperture interferogram; and finally, calculating to obtain the earth surface elevation data according to the satellite orbit parameters and the wrapped phase of the full-aperture interference diagram with the troposphere delay error phase removed.

2. The method of claim 1, wherein the two sub-view interferograms are topographically processed using two DEM data obtained from different techniques.

3. The method of claim 1, wherein the wavelet decomposition of the sub-view differential interferogram is represented as:

in the formula, m and n respectively represent the number of rows and columns of the interference pattern, x and y are coordinates of pixel points in a radar coordinate system, i_xAnd i_yI th representing interference patterns, respectively_xRow and ith_yColumns, Φ and Ψ respectively represent the translation function and mother wavelet function of the wavelet decomposition, J represents the optimal decomposition scale, J represents different decomposition scales, v, w respectively represent the low-frequency wavelet coefficient and high-frequency wavelet coefficient, and ∈ 1,2,3 respectively corresponds to the horizontal, vertical and diagonal directions;

the correlation analysis method in the step 4 is that the high-frequency wavelet coefficient of the two sub-view difference interferograms under each decomposition scale is taken

And

the correlation value calculation is performed according to the following formula:

high-frequency wavelet coefficient when epsilon is 1,2,3 obtained by wavelet decomposition

And

then, the correlation value f of the two sub-view differential interferograms in the decomposition scale j can be obtained by calculation of an equation set by substituting the correlation value f into the correlation calculation formula_jAnd global offset coefficient c_j；

Reuse of the correlation f of two sub-view differential interferograms at the decomposition scale j_jAnd global offset coefficient c_jCalculating a high-frequency wavelet coefficient corresponding to the troposphere delay error phase according to the following formula:

in the formula (I), the compound is shown in the specification,

representing a high-frequency wavelet coefficient corresponding to the troposphere delay error phase;

indicating an azimuthal angle of incidence of α_iThe sub-parallax partial interferograms are high frequency wavelet coefficients at decomposition scale j.

4. The method of claim 1, wherein when removing the flat phase, the calculation formula of the flat phase is:

in the formula, phi_flatIs the phase of the flat ground, λ is the wavelength of the radar signal, B_//Is the parallel baseline length.

5. The method of claim 1, wherein the track error phase is removed by: and fitting the orbit error phase by using a polynomial model fused into the earth surface elevation h, resolving unknown parameters in the polynomial model by using a least square method, and then removing the orbit error phase from the unwrapping phase of the corresponding pixel point in the sub-parallax partial interferogram to obtain a new sub-parallax partial interferogram.

6. The method of claim 1, wherein the step 1 of decomposing the full aperture SAR image by using the sub-aperture decomposition technique comprises:

firstly, carrying out one-dimensional Fourier transform on a full-aperture SAR image to obtain a frequency spectrum of the SAR image in the azimuth direction;

then, dividing the frequency spectrum of the SAR image in the azimuth direction according to different azimuth direction incidence angles to obtain a frequency spectrum corresponding to the azimuth direction incidence angle;

and finally, performing inverse Fourier transform on all frequency spectrums corresponding to the azimuth incident angles to obtain images corresponding to the azimuth incident angles, namely the sub-view images corresponding to the azimuth incident angles.

7. The method of claim 1, wherein preprocessing the full aperture SAR image comprises: co-registering the 2 full-aperture SAR images, and then carrying out pixel-by-pixel conjugate multiplication to generate a full-aperture interferogram;

the preprocessing of the 2 sub-view images corresponding to each azimuth incident angle of the main image and the auxiliary image comprises the following steps: and carrying out pixel-by-pixel conjugate multiplication on the 2 sub-view images corresponding to the current azimuth incident angle to obtain a sub-view interference image corresponding to the current azimuth incident angle.

8. A single baseline-based surface elevation correction apparatus, comprising:

a sub-aperture decomposition module to: acquiring two full-aperture SAR images, and decomposing each full-aperture SAR image into two sub-view images corresponding to azimuth incident angles by adopting a sub-aperture decomposition technology;

a data pre-processing module to: preprocessing the two full-aperture SAR images to obtain a full-aperture interferogram; preprocessing two sub-view images with the same azimuth incident angle to obtain a sub-view interferogram corresponding to the azimuth incident angle;

a differential interferometric phase modeling module to: for the two sub-view interferograms, DEM data are used, a differential interference technology is adopted for terrain removing processing, then flat ground phase and orbit error phase removing processing are carried out, and two sub-parallax partial interferograms shown in the following formulas are obtained:

in the formula, α₁And α₂Respectively representing two different azimuthal angles of incidence, phi_diffIs the differential interference phase in the sub-parallax partial interference pattern, △ phi_topoFor the phase of the terrain error, phi_atmoFor tropospheric delay error phase, phi_noiseIs a random noise error phase;

a tropospheric delay error phase extraction module for: performing wavelet decomposition on the two sub-view differential interferograms, performing correlation analysis on high-frequency wavelet coefficients of the two sub-view differential interferograms under each decomposition scale to obtain correlation values between the two high-frequency wavelet coefficients under the corresponding decomposition scale, and calculating to obtain wavelet coefficients related to troposphere delay error phases according to the correlation values and any high-frequency wavelet coefficient; then, performing inverse wavelet transform on the wavelet coefficient related to the troposphere delay error phase to obtain a troposphere delay error phase;

a surface elevation data correction module to: carrying out unwrapping processing on the full-aperture interferogram, and removing the obtained troposphere delay error phase from the unwrapping phase of the full-aperture interferogram; and finally, calculating to obtain the earth surface elevation data according to the satellite orbit parameters and the wrapped phase of the full-aperture interference diagram with the troposphere delay error phase removed.

9. An apparatus comprising a processor and a memory; wherein: the memory is to store computer instructions; the processor is configured to execute the computer instructions stored by the memory, in particular to perform the method according to any one of claims 1 to 7.

10. A computer storage medium storing a program which, when executed, performs the method of any one of claims 1 to 7.

Images