阅读说明：本技术 一种基于小波多尺度相关性分析的机载InSAR轨道误差去除方法 (Airborne InSAR orbit error removing method based on wavelet multi-scale correlation analysis ) 是由 何永红 于 2020-04-09 设计创作，主要内容包括：本发明公开了一种基于小波多尺度相关性分析的机载InSAR轨道误差去除方法,涉及遥感影像的微波遥感技术领域。本发明方案基于轨道误差独立于极化方式,将机载InSAR轨道误差处理理论与方法从单极化模式拓展到双极化模式,对不同极化方式干涉图中的差分干涉相位进行多尺度分析,滤除比轨道误差相位波长短的地形误差相位、噪声相位等成分,并对不同极化方式差分干涉相位轨道误差的高相关性进行降权精化,去除机载InSAR轨道误差相位,该方法理论严密,且不需要其他数据进行支撑,能够有效提高InSAR数据处理精度。(The invention discloses an airborne InSAR orbit error removing method based on wavelet multi-scale correlation analysis, and relates to the technical field of microwave remote sensing of remote sensing images. The scheme of the invention is based on that the orbit error is independent of the polarization mode, the onboard InSAR orbit error processing theory and the onboard InSAR orbit error processing method are expanded from a single polarization mode to a dual polarization mode, the differential interference phases in the interferograms of different polarization modes are subjected to multi-scale analysis, the components such as the terrain error phase and the noise phase which are shorter than the orbit error phase are filtered, the high correlation of the orbit errors of the differential interference phases of different polarization modes is subjected to weight reduction refinement, and the onboard InSAR orbit error phase is removed.)

1. An airborne InSAR orbit error removing method based on wavelet multi-scale correlation analysis is characterized by comprising the following steps:

step S1: acquiring a synthetic aperture radar SAR main image and an auxiliary image of two different polarization modes corresponding to a target area, and generating interferograms of different polarization modes through image pre-filtering, main-auxiliary image registration and other processes; calculating flat ground phase using orbit parametersCalculating terrain phase using orbit parameters and external DEM dataFiltering the interferograms in the two different polarization modes, performing differential interferometry, and removing the flat ground phase and the terrain phase to obtain the differential interferograms in the two different polarization modes

Step S2: differentially interfering the phase patterns of the two different polarization modesRespectively performing wavelet multi-scale analysis, and decomposing and reconstructing signals by using waveletsSeparating high frequency decorrelation noise phase by using a decomposition scale determined by the variation rate of the Root Mean Square Error (RMSE) of the original signal as a high and low frequency optimal decomposition scale MDetermining a maximum decomposition scale L by using the change rate of low-frequency part Average Correlation Coefficient (ACC) in two different polarization modes, and separating out a low-frequency terrain residual error phaseAnd track error phase

Step S3: refining the low-frequency wavelet coefficient of the low-frequency part on the scale L of the differential interferogram in the two different polarization modes according to an average correlation coefficient weighting mode; refining the high-frequency wavelet coefficient in an average correlation coefficient weighting mode for the high-frequency part between the scales L and M; and (5) reconstructing by using the refined wavelet coefficient to obtain the differential interference phase after the track error is corrected.

Wherein the flat ground phaseThe expression is as follows:

wherein λ is the wavelength, R₁And R₂Is a reference slope distance;

the terrain phaseThe expression is as follows:

in the formula, P represents a polarization mode, and may be a linear polarization mode (e.g., HH, HV, VV, HH + VV, and HH-VV or a synthetic polarization mode), B is a base length, θ is an incident angle, α is a base inclination angle, h is a bare terrain phase, and Δ h is a height associated with a scattering target.

2. The method according to claim 1, wherein the step S1 is to obtain two different polarization mode differential interferogramsThe interference phase components are as follows:

wherein P represents a polarization mode such as HH, HV, VV, HH + VV, and HH-VV, or a synthetic polarization mode.

Representing the phase of the terrain residual caused by external DEM errors, in relation to the polarisation mode.

Indicating the noise phase, depending on the polarization mode.

Indicating the track error phase regardless of the polarization mode.

3. The method according to claim 1, wherein the wavelet multi-scale analysis of the differential interference phase in step S2 is calculated by the following formula:

in the formula, HH (x)₁,x₂) And VV (x)₁,x₂) Respectively representing the differential interference phases of HH and VV in two different polarization modes, phi and psi are respectively a scale function and a wavelet function, J is 1-J, J is a wavelet decomposition scale, <, > is an inner product operator,andis a low-frequency wavelet decomposition coefficient of interference phases of HH and VV in two different polarization modes,andthe high-frequency wavelet coefficients with different interference phases and different resolutions in HH and VV polarization modes. E is 1,2,3, each is f (x)₁,x₂) Wavelet decomposition coefficients of high frequency components in vertical, horizontal and diagonal directions.

4. The method according to claim 1, wherein the decomposition scale of step S2 is calculated by the following formula:

the wavelet decomposition reconstruction signal and the original signal Root Mean Square Error (RMSE) change rate determine the optimal decomposition scale of high and low frequencies, and the expression is as follows:

in the formula

The Average Correlation Coefficient (ACC) rate of change determines a maximum decomposition metric, which is expressed as:

in the formula, an Average Correlation Coefficient (ACC) is expressed as:

wherein the content of the first and second substances,<>the average of the windows is opened,is the average value of the high-frequency wavelet coefficient in the window.

5. The method according to claim 1, wherein the mean correlation coefficient weighting manner of step S3 is calculated by the following formula:

low-frequency part refinement mode:

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

high-frequency partial refinement:

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

6. the method according to claim 1, wherein the step S3 of reconstructing the wavelet coefficients of differential interference phase is calculated by the following formula:

in the formulaWhich represents the low-frequency coefficients after the correction,in order to correct the detail coefficients in different directions,the low-frequency coefficient before the correction is carried out,for the different direction detail coefficients after correction, epsilon is 1,2, 3.

Technical Field

The invention relates to the technical field of microwave remote sensing of remote sensing images, in particular to an airborne InSAR orbit error removing method based on wavelet multi-scale correlation analysis.

Background

At present, an airborne InSAR system plays an important role in the fields of forest parameter reproduction, crop monitoring, geological disasters and the like. However, in an Airborne interference (Airborne SAR interference) system, due to the influence of unstable factors such as airflow, position deviation and attitude change occur during the flight of an Airborne vehicle, and the resulting motion error affects the phase center position of a radar antenna and the direction of a radar beam, thereby affecting the imaging quality and the elevation inversion accuracy. In order to establish a strict geometric relationship between an interference phase and a monitored target position, an airborne SAR platform inertial navigation system is required to be capable of accurately recording flight platform motion parameters so as to carry out motion compensation processing in a subsequent data processing process. However, the precision of the navigation positioning system of the carrier can only reach 2-15cm even after DGPS processing, the image compensated by POS data still has residual phase error caused by residual time-varying baseline error, especially for repeated track interferometry, even if a high-precision inertial navigation system is adopted, the track error in each flight process of the carrier is independent, unknown time-varying baseline error is introduced, and in the process of generating an interferogram, the track error caused by baseline variation cannot be offset, and serious track error phases (residual motion error phases) can be generated in the distance direction and the azimuth direction. With the development of high-resolution airborne InSAR, the orbital error phase becomes a determining factor for restricting the improvement of the resolution of the airborne InSAR, and as an important error source, the orbital error phase brings difficulty to the correct interpretation and information extraction of an interferogram and is one of the main limiting factors of the accuracy of interferometry.

In order to eliminate the track error phase, various methods have been developed. The method mainly comprises the following steps: a method based on sub-aperture decomposition (MS), a method based on SAR image self-focusing (WPCA) and a combined calibration method based on multi-baseline InSAR data. However, the method based on the sub-aperture decomposition technique has strong dependence on coherence, and when the monitored ground object geometry has anisotropy, the difference of phase centers can reduce the precision of the method under different incidence angles; a WPCA method is based on, and a high-coherence target point is required in a monitoring scene; the joint calibration method based on multiple baselines requires a large data volume to obtain reliable results, and the universality of the method is greatly limited. Therefore, the airborne InSAR orbit error removal model is constructed, and the problem of airborne SAR data orbit error under the condition of a small amount of data is solved, so that the method has important significance for improving the airborne interference quality.

Disclosure of Invention

In order to solve the problems in the prior art, the invention provides an airborne InSAR orbit error removing method based on wavelet multi-scale correlation analysis. The method expands the InSAR orbit error estimation method from a single polarization data mode (HH or VV) to a dual polarization data mode (HH, HV, VV, VH), solves the problem of the orbit error of the airborne SAR data under the condition of a small amount of data, plays an important role in improving the interpretation precision and application of the InSAR result, and provides a new idea for the compensation of the orbit error of the airborne InSAR.

The invention provides an airborne InSAR orbit error estimation and removal method based on wavelet multi-scale correlation analysis, which comprises the following steps:

step S1: acquiring a synthetic aperture radar SAR main image and an auxiliary image of two different polarization modes corresponding to a target area, and generating interferograms of different polarization modes through image pre-filtering, main-auxiliary image registration and other processes; calculating flat ground phase using orbit parametersCalculating terrain phase using orbit parameters and external DEM dataFiltering the interferograms in the two different polarization modes, performing differential interferometry, and removing the flat ground phase and the terrain phase to obtain the differential interferograms in the two different polarization modes

Step S2: differentially interfering the phase patterns of the two different polarization modesRespectively performing wavelet multi-scale analysis, and separating high-frequency decorrelation noise phase by using a decomposition scale determined by the variation rate of the wavelet decomposition reconstruction signal and the Root Mean Square Error (RMSE) of the original signal as a high-frequency and low-frequency optimal decomposition scale MUsing two different polarizationsDetermining the maximum decomposition scale L by the change rate of the mean correlation coefficient (ACC) of the mode low-frequency part, and separating out the low-frequency terrain residual error phaseAnd track error phase

Step S3: refining the low-frequency wavelet coefficient of the low-frequency part on the scale L of the differential interferogram in the two different polarization modes according to an average correlation coefficient weighting mode; refining the high-frequency wavelet coefficient in an average correlation coefficient weighting mode for the high-frequency part between the scales L and M; and (5) reconstructing by using the refined wavelet coefficient to obtain the differential interference phase after the track error is corrected.

Wherein the flat ground phaseThe expression is as follows:

wherein λ is the wavelength, R₁And R₂Is a reference slope distance;

the terrain phaseThe expression is as follows:

in the formula, P represents a polarization mode, and may be a linear polarization mode (e.g., HH, HV, VV, HH + VV, and HH-VV or a synthetic polarization mode), B is a base length, θ is an incident angle, α is a base inclination angle, h is a bare terrain phase, and Δ h is a height associated with a scattering target.

According to the scheme, under the multi-polarization PolInSAR data mode, the track error signals are irrelevant to the polarization mode, SAR images in different polarization modes are acquired simultaneously, so that interferograms in different polarization modes contain the same track error phase, and after the terrain residual error phase and the noise phase of the differential interference phase are restrained by utilizing wavelet multi-scale analysis, the weight reduction refinement of the track error phase through the high correlation of the track errors in different polarization modes is realized.

In the scheme of the present invention, step S1 specifically includes removing the preprocessing processes such as the flat ground phase and the terrain phase, and finally obtaining the differential interferograms of the two different polarization modes, and analyzing whether the phase component compositions of the differential interferograms of the two different polarization modes and the polarization modes are related to each other, so as to clear up the correlation between the phase components in the interference phase compositions and the polarization modes, and provide theoretical support for the subsequent step S2 of establishing a correction model.

For multi-scale analysis, it is difficult to select a proper wavelet function to completely distinguish the terrain error phase from the track error phase, but the wavelet decomposition technology can primarily separate the terrain error phase and the noise phase with the wavelength smaller than the track error phase. In fact, the noise phase and the terrain error phase are mainly composed of short-wavelength components, and the track error can be better detected through the separation of the components.

The step of performing the multi-scale analysis in step S2 is: the interference phase is first separated into high-frequency phase decorrelation noise phases through two-dimensional wavelet decomposition, and a decomposition scale when the variation rate of a Root Mean Square Error (RMSE) between a literature reconstruction signal and an original signal is about 1 is adopted as a high-low frequency optimal decomposition scale M.

After the interference phase is decomposed and separated into high-frequency phase loss correlation noise phases through two-dimensional wavelet, the difference interference phases with different polarizations are further subjected to multi-scale analysis, so that the terrain error phases can be further removed. To remove the terrain error phase with the wavelength smaller than the orbit error phase, the decomposition scale L is determined by adopting a low-frequency partial average correlation coefficient ACC (averaging correction coefficient) method in a two-polarization mode.

By setting a suitable resolution L, the terrain error phase can be partially removed, but still a residual terrain error phase is mixed in the track error phase.

The step of performing the right reducing correction in step S3 is:

the low-frequency part of the scale L is subjected to weight reduction and refinement,

in the formula (I), the compound is shown in the specification,which represents the low-frequency coefficients after the correction,low frequency coefficient before correction, c_jThe low-frequency correlation coefficient of the two polarization modes is shown,

when the decomposition scale is larger than the optimal decomposition scale, as the terrain residual error and the orbit error have the frequency aliasing condition, part of the orbit error is decomposed into high-frequency components, and the weight reduction and the refinement are carried out on the high-frequency components on the scale when L is larger than or equal to M according to the correlation weighting mode.

In the formulaIn order to correct the detail coefficients in different directions,in order to correct the detail coefficients in different directions,the two-dimensional correlation coefficient in different directions of the high-frequency image of the two polarization modes is shown, and epsilon is 1,2 and 3.

And (5) reconstructing by using the refined wavelet coefficient to obtain the differential interference phase after the track error is corrected and the estimated track error.

Drawings

Fig. 1 is a schematic flow diagram of an airborne InSAR orbit error removal method based on wavelet multi-scale correlation analysis according to an embodiment of the present invention;

FIG. 2 is a schematic diagram of InSAR altimetry principle in the embodiment of the invention;

FIG. 3 is a flow chart of an algorithm of an airborne InSAR orbit error removal method based on wavelet multi-scale correlation analysis in an embodiment of the present invention;

FIG. 4 is a graph showing the experimental results of the variation rate of the measured data HH polarization mode RMSE in the P band in the embodiment of the present invention;

FIG. 5 is a graph showing the result of an ACC rate test of two polarization modes HH-HV of measured data in P band according to an embodiment of the present invention;

FIG. 6 is a graph showing the experimental results of the correlation between the coefficients of the low-frequency and high-frequency wavelets of the P-band HH-HV dual-polarized interferogram in the embodiment of the present invention;

FIG. 7 is a diagram of the lateral track error phase profile results extracted using an HH-HV dual polarization interferogram in an embodiment of the present invention.

FIG. 8 is a diagram of the phase profile of longitudinal track errors extracted using an HH-HV dual polarization interferogram in an embodiment of the present invention.

FIG. 9 is a graph showing the experimental results of track error phase estimation using HH-HV dual polarization interferograms in an embodiment of the present invention.

Detailed Description

In order to facilitate an understanding of the teachings of the present invention, reference will now be made in detail to the following detailed description.

In order to make the purpose, technical scheme and advantages of the invention more clear, the experiment of the invention adopts the full polarization P wave band actual measurement data acquired by the airborne E-SAR of the Germany space station to carry out verification, and the data is acquired by adopting a repeated orbit flight mode under the support of the BioSAR2008 project. The experimental zone being located in the northern part of SwedenProvince, images are acquired in 10 months in 2008, interference base lines are 16 meters, terrain change of the area is about 400 meters, vegetation resources are rich, main tree species are pine trees, spruce trees, birch trees, aspens, twisted-leaf pines, ash trees, yellow willows and Chinese ash trees, and HH and HV dual-polarized differential interferograms and HH and VV dual-polarized differential interferograms are respectively adopted to estimate track error phases.

The method comprises the following specific steps:

step 1: accurately registering main and auxiliary images of SAR phases of different polarization modes of HH, VV and HV, removing the ground phase, and performing interference processing to obtain an original interference pattern; filtering the original interferogram, and differentiating the interferogram and the simulated terrain phase generated by the DEM to obtain a differential interferogram; and then generating differential interference patterns of different polarization modes such as HH, VV, HV and the like by phase unwrapping.

Wherein, the land phase expression is:

the terrain phase expression is:

the differential interference phase composition is obtained as follows:

in the formula, P represents a polarization mode, and is a linear polarization mode of HH, HV, or VV.

Step 2: and performing wavelet multi-scale decomposition and reconstruction on the differential interference image after the unwrapping of a certain polarization mode, such as an HH polarization differential interference image, and taking a decomposition scale determined by the RMSE (RMSE change rate) of a reconstructed signal and an original signal as an optimal scale M for high-low frequency decomposition.

Performing wavelet multi-scale decomposition and reconstruction on the HH polarization differential interference pattern:

HH(x₁,x₂) The differential interference phase representing HH polarization mode, phi and psi are scale function and wavelet function respectively, J is 1-J, J is wavelet decomposition scale, < -, > is inner product operator,andis a low-frequency wavelet decomposition coefficient of interference phases of HH and VV in two different polarization modes,andthe high-frequency wavelet coefficients with different interference phases and different resolutions in HH and VV polarization modes. E is 1,2,3, each is f (x)₁,x₂) Wavelet decomposition coefficients of high frequency components in vertical, horizontal and diagonal directions.

Calculating the RMSE change rate of the reconstructed signal and the original signal in the wavelet multi-scale analysis process:

in the formula

The decomposition scale when the rate of change of the reconstructed signal and the original signal RMSE is close to 1 is taken as the high and low frequency decomposition optimal scale M.

And step 3: and respectively carrying out wavelet multi-scale decomposition on the differential interference patterns such as HH and HV after the two different polarization modes are unwrapped, calculating the ACC change rate of the low-frequency coefficients of the HH and the HV in the two different polarization modes, and determining the maximum scale L of the low-frequency decomposition.

Performing wavelet multi-scale decomposition on HH and VV polarization differential interferograms respectively:

calculating ACC change rates of HH and HV low-frequency coefficients of two different polarization modes:

the Average Correlation Coefficient (ACC) expression is:

and when the ACC change rate of the low-frequency coefficient of the differential interferograms of the two different polarization modes tends to 1, the ACC change rate is taken as a low-frequency decomposition maximum scale L.

And 4, step 4: and performing refined weight reduction processing on a low-frequency part on the HH polarization mode with the scale j of the differential interferogram as L and a high-frequency part on the HH polarization mode with the scale j as M, … and L by adopting a weighting mode.

Low-frequency part refinement mode:

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

high-frequency partial refinement:

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

and 5: and performing wavelet reconstruction by using the wavelet coefficient subjected to the weight reduction processing to reconstruct a differential interference phase subjected to track error correction in the HH polarization mode.

In the formulaWhich represents the low-frequency coefficients after the correction,in order to correct the detail coefficients in different directions,the low-frequency coefficient before the correction is carried out,for the different direction detail coefficients after correction, epsilon is 1,2, 3.

Step 6: in order to verify the effectiveness of the track errors estimated by the HH polarization mode and the HV polarization mode, the track errors are estimated by using the HH polarization mode and the VV polarization mode, and the root mean square error of the track errors estimated by the HH polarization mode and the HV polarization mode and the track errors estimated by the HH polarization mode and the VV polarization mode is calculated.

And 7: and generating the DEM by using the original interferogram and the corrected interferogram, and evaluating the DEM result by using LiDAR as a true value.

From the experimental result graph, the method has better fitting with the original phase at the position with smaller orbit error or close to 0, and has good correction effect at the position with larger orbit error in the azimuth direction, the corrected interference phase is jittered near 0, and most of the orbit error in the azimuth direction existing in the interference fringe pattern is eliminated. HH. The root mean square error between the orbit error estimated by the HV polarization mode and the orbit error estimated by the HH and VV dual-polarization differential interferograms is 0.0065, and the corrected phase has better continuity, consistency and robustness in both the distance direction and the azimuth direction. According to the DEM generated by the corrected interferogram, the InSAR height measurement precision is obviously improved, and the effectiveness of the method is proved.

In summary, the invention provides an airborne InSAR orbit error removal method based on wavelet multi-scale correlation analysis, and the method expands an airborne InSAR orbit error processing theory and method from a single polarization mode (HH or VV) to a dual polarization data mode (HH, HV, VV, VH), thereby making up the problem of airborne InSAR data orbit error under the condition of a small amount of data. The method can solve the restriction of airborne orbit error phase on SAR interferometric measurement precision to a certain extent, and has important significance on InSAR technical theory and measurement method research. Meanwhile, the method has important practical significance for improving the InSAR data processing capacity and enhancing the practicability of the InSAR measurement technology. The method also has important social significance for accurately monitoring the ground settlement characteristics of key areas in China and carrying out potential geological disasters and environmental evaluation, and can be widely applied.

The above description is only exemplary of the present invention and should not be taken as limiting the invention, as any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.