InSAR nonlinear despeckling method

Abstract

The application discloses an InSAR nonlinear flat-removing effect method and relates to the field of synthetic aperture radar interference processing. The synthetic aperture radar detects a target to be measured to obtain an interference phase diagram and an interference complex image, the interference complex image is pre-filtered, discrete Fourier transform is carried out and summation is performed, and a spectrum amplitude sequence and a phase sequence corresponding to a distance direction are obtained. Then, a normalized histogram is constructed from the spectrum amplitude sequence, a threshold value is adaptively selected according to the histogram distribution, the mean value of a part of the spectrum amplitude sequence which is lower than the threshold value is subtracted from a part of the spectrum amplitude sequence which is higher than the threshold value, a processed spectrum amplitude sequence is obtained, and the spectrum phase sequence is recombined with the processed spectrum amplitude sequence, inverse discrete Fourier transform is performed, a phase angle is taken, a nonlinear flat phase sequence is obtained, and the nonlinear flat phase sequence is expanded to the original size. The nonlinear flat phase sequence expanded to the original size is subtracted from the interference phase diagram, and a final flat-removing result is obtained. The flat-removing result is more accurate and is suitable for different scenes.

Description

Technical Field

[0001] This invention relates to the field of synthetic aperture radar interferometric processing, specifically to an InSAR nonlinear flat-ground effect removal method that compensates for the interferometric phase and removes the nonlinear flat-ground phase in the scene. Background Technology

[0002] Interferometric Synthetic Aperture Radar (InSAR) technology can acquire ground elevation and deformation information over a wide area under all-weather and all-time conditions. It is often used to obtain high-precision digital elevation model (DEM) data and has important application and research value in many fields such as surface deformation detection, moving target detection, marine mapping, forest mapping, flood monitoring, traffic monitoring, and glacier research.

[0003] Based on spatial geometry, even flat terrain with constant elevation will produce dense, periodic fringes in the InSAR interferometric phase map, causing the interferometric phase to fail to reflect actual ground elevation changes. Therefore, it is necessary to eliminate this flat terrain phase, reduce the phase gradient, and facilitate subsequent phase filtering and unwrapping. The process of removing the flat terrain phase is called flat terrain effect removal, or flat terrain phase compensation.

[0004] Currently, commonly used techniques for removing flat-ground effects are mainly divided into three types: those based on orbital parameters, those based on external DEM data, and those based on interference fringe frequencies. In practical scenarios, it is sometimes difficult to obtain the corresponding accurate orbital data or external DEM data. In such cases, the third method can be used, which involves estimating and compensating for the frequency of the brightest fringe in the phase diagram, i.e., frequency shifting, to remove the flat-ground phase. This type of method does not require data other than the interference phase diagram and has a faster computation speed, making it the most commonly used method for removing flat-ground effects. However, since it only compensates for a linear phase that is fixed as a constant, it is only suitable for estimating uniformly distributed flat-ground fringes.

[0005] When the viewing angle is small or the slope angle is large, the interferometric phase fringes will exhibit a clear phenomenon of dense fringes near the edge and sparse fringes far from the edge. In this case, the flat terrain phase is no longer a uniformly distributed linear phase, but rather a scenario of dense fringes near the edge and sparse fringes far from the edge. For such scenarios, if the traditional flat terrain removal method is still used, the estimated flat terrain phase will differ significantly from the true phase, and the interferometric phase map after flattening will still have dense fringes in some areas, affecting subsequent interferometric processing. Although the block-based method based on the frequency shift method can improve the above problems to some extent, the principle of block division for different scenarios is often difficult to determine, and it is easy to introduce local estimation errors caused by abrupt changes in local terrain, which actually increases the workload. Summary of the Invention

[0006] In order to remove the flat-ground phase of scenes with small downward view angles or large slope angles as accurately as possible under unsupervised conditions and reduce the difficulty of subsequent interferometric processing, this invention provides an InSAR nonlinear flat-ground phase removal method, which can remove the flat-ground phase of the scene in an unsupervised, efficient and accurate manner even in the absence of additional orbit and DEM data.

[0007] The InSAR nonlinear flat-land effect removal method includes the following steps:

[0008] Step 1: The synthetic aperture radar (SAR) performs two imaging operations on the same target to obtain two complex images, which are then registered. The two registered SAR complex images are then interferometric to obtain the corresponding interferometric phase diagram.

[0009] In the interferometric phase diagram, the noisy phase Decomposed into flat phase Residual phase with changes in reaction details

[0010] Step 2: Pre-filter the complex interferometric image, and perform discrete Fourier transform on the range data of the pre-filtered complex interferometric image and sum them to obtain the spectral amplitude sequence and phase sequence corresponding to the range direction.

[0011] Based on the Discrete Fourier Transform, the formula for the complex interferometric phase decomposition of the pre-filtered complex interferometric image along the range direction is as follows:

[0012]

[0013] in, It is the complex form of the noisy phase. The imaginary unit is , m is the coordinate position, G represents the number of vectors that make up the complex interference phase, and C... i The weights of the vector. This represents the range frequency of the corresponding vector.

[0014] In each set of distance directions, the complex interference phase sequence is paired. Perform a discrete Fourier transform and sum the results to obtain the spectral amplitude sequence S. i The summation result yields the spectrum S reflecting the overall distance distribution. r And the corresponding spectral amplitude sequence P and phase sequence Φ:

[0015]

[0016] in Let represent the i-th distance-up complex interference phase sequence, u be the coordinate position of the spectral sequence, and ⊙ denote the Shur product operator.

[0017] Step 3: Construct a normalized histogram based on the spectrum amplitude sequence, and adaptively select a threshold based on the histogram distribution.

[0018] The process of adaptively selecting the threshold based on the histogram is as follows:

[0019] (1) Construct a normalized histogram based on the spectrum amplitude sequence, perform a second-order difference operation on the normalized histogram, and find all the maximum points in it;

[0020] Perform second-order difference on the histogram sequence H:

[0021] D(k)=[H(k+1)-H(k)]-[H(k)-H(k-1)],k=2,3,...,N-1 (3)

[0022] Where D is the second-order difference sequence and N is the number of histogram groups, all the maximum points in the histogram can be found by judging that D(k)≥0.

[0023] (2) Starting from the left end, determine the relationship between the frequencies corresponding to the maximum points and the frequencies under a uniform distribution in the histogram. This represents the frequency level of the noise; if This represents the frequency level of the flat-ground phase. Among the maxima that satisfy the flat-ground phase frequency condition, the maxima with the smallest amplitude in the corresponding histogram is selected as the threshold T.

[0024] Step 4: For the portion of the spectral amplitude sequence below the threshold, take its mean as the estimated background noise power value and set this portion to zero; for the portion above the threshold, subtract the estimated background noise power value to further suppress the influence of noise. After the above operations, the processed spectral amplitude sequence is obtained.

[0025] The background noise power estimate is: noise = mean(P(u) < T), where mean(·) represents the mean operation.

[0026] The processed spectral amplitude sequence is as follows:

[0027]

[0028] Step 5: Recombine the processed spectral amplitude sequence and spectral phase sequence, perform an inverse discrete Fourier transform, and obtain the nonlinear flat-ground phase sequence by taking the phase angle, and then expand it to the original size.

[0029] The nonlinear flat-ground phase sequence is as follows:

[0030] Φ p =arg(IFFT(P′⊙Φ)) (5)

[0031] The formula for expanding a nonlinear flat-ground phase sequence to its original size is:

[0032]

[0033] Step 6: Subtract the nonlinear flat phase sequence expanded to the original size from the interferometric phase diagram to obtain the final flattened result.

[0034] The advantages of this invention compared to the prior art are:

[0035] (1) Compared with the traditional frequency shift method, the present invention can estimate the nonlinear flat phase. In small-view systems or complex terrain, the flat phase result is more accurate. The interference fringes after flat phase are sparser, which is beneficial to subsequent phase filtering and unwrapping.

[0036] (2) Compared with the block method, the present invention can adaptively select the threshold and is applicable to different scenarios; and since the present invention takes into account the spectrum of the entire distance upward, it can effectively eliminate the spectrum low peak caused by local terrain changes, and is not easily affected by local areas. The whole process does not require manual intervention. Attached Figure Description

[0037] Figure 1 This is a flowchart of the InSAR nonlinear flat-ground effect removal method of the present invention;

[0038] Figure 2 This is a noisy interferometric phase map generated in an embodiment of the present invention;

[0039] Figure 3 This is a real flat-ground phase map generated in the embodiments of the present invention;

[0040] Figure 4 This is a range-oriented spectrum generated in an embodiment of the present invention;

[0041] Figure 5 This is the normalized histogram generated in the embodiments of this invention;

[0042] Figure 6 This is a nonlinear flat-ground phase map generated in an embodiment of the present invention;

[0043] Figure 7 This is a diagram showing the results of leveling the ground generated in an embodiment of the present invention;

[0044] Figure 8 This is the distance spectrum of the deflating results generated in the embodiments of the present invention;

[0045] Figure 9 This is a flat-ground phase map generated by a traditional frequency-shifting method implementation.

[0046] Figure 10 It is a residual phase map generated by a traditional frequency shift method implementation;

[0047] Figure 11 This is the range spectrum of the deflattening result generated by a traditional frequency shift method implementation;

[0048] Figure 12 This is a flat-ground phase map generated by a block-based implementation method;

[0049] Figure 13 It is the residual phase map generated by the block method implementation;

[0050] Figure 14 This is the distance spectrum of the deflating result generated by the block method implementation example. Detailed Implementation

[0051] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments.

[0052] A nonlinear flat-ground effect removal method for InSAR, flowchart as follows: Figure 1 As shown, it includes the following steps:

[0053] Step 1: The synthetic aperture radar (SAR) performs two imaging operations on the same target to obtain two complex images, which are then registered. The two registered SAR complex images are then interferometric to obtain the corresponding interferometric phase diagram.

[0054] The noisy interferometric phase diagram of the small-angle marine scene obtained by radar parameter simulation is shown in Table 1. Figure 2 As shown, its true flat-ground phase is as follows Figure 3 As shown, the fringes are relatively denser on the left side of the flat-ground phase and relatively sparser on the right side.

[0055] Table 1 Radar Parameters

[0056]

[0057] In the interferometric phase diagram, the density of the fringes is determined by the dominant flat-land phase and the secondary local topographic detail phase and noise. Therefore, the noisy phase can be analyzed based on the additive noise model. Decomposed into flat phase Residual phase with changes in reaction details

[0058]

[0059] Step 2: Pre-filter the complex interferometric image, perform a discrete Fourier transform on the pre-filtered complex interferometric image in the range direction and sum the results to obtain the spectral amplitude sequence and spectral phase sequence corresponding to the range direction.

[0060] According to the Discrete Fourier Transform, the complex interference phase can be decomposed into a combination of multiple linear phases. Considering that the phase on flat terrain is generally composed of parallel fringes along the range direction, the decomposition formula for the complex interference phase distributed along the range direction is given without loss of generality:

[0061]

[0062] in, It is the complex form of the noisy phase. The imaginary unit is , m is the coordinate position, G represents the number of vectors that make up the complex interference phase, and C... i The weights of the vector. This represents the range frequency of the corresponding vector.

[0063] During the estimation process, to prevent abrupt changes in local terrain in the image from interfering with the estimation of the flatland phase, I n Instead of simply selecting a single range direction data point as a sample, the range direction frequency characteristics of the entire image in all directions should be comprehensively considered. Therefore, it is necessary to analyze the spectral sequence S of each range direction. i Summation is performed to obtain the spectrum S that reflects the overall distance distribution. r And the corresponding spectral amplitude sequence P and phase sequence Φ:

[0064]

[0065] in Let represent the i-th distance-up complex interference phase sequence, u be the coordinate position of the spectral sequence, and ⊙ denote the Shur product operator.

[0066] The spectral amplitude sequence P of the simulated image is as follows Figure 4 As shown, it can be seen that its spectral peak is shifted relative to the zero frequency and has a certain width.

[0067] Since the flat-land phase is the main component of the phase, it can be estimated by selecting the top g vector combinations with larger weights based on the magnitude of the vector weights:

[0068]

[0069] Reflected in S r In this context, the magnitude of the vector weight is the same as the magnitude of P, and the range-direction frequency characteristics can be determined by Φ. The key question then becomes how to determine g to estimate the flat-ground phase.

[0070] Step 3: Construct a normalized histogram based on the spectrum amplitude sequence, and adaptively select a threshold based on the histogram distribution.

[0071] For interferometric phase maps in general scenarios, the spectrum often exhibits a single peak, while multiple secondary peaks may appear when the terrain is complex. Traditional frequency-shifting methods rely solely on the location of the peak point, and when the spectral peak is wide or secondary peaks exist, some flat-ground phase information is easily lost, leading to a mismatch between the estimated flat-ground phase and the true flat-ground phase. As the preceding analysis shows, the flat-ground phase, as the main influencing factor of the interferometric phase map, should dominate the spectrum, distributed in the low-frequency region with a high spectral amplitude; while the residual phase plays a secondary role, distributed in the high-frequency region with a generally lower spectral amplitude. Therefore, a threshold can be selected based on the spectral amplitude to distinguish between them.

[0072] The specific process for selecting a threshold based on the spectral amplitude is as follows:

[0073] First, by constructing a normalized histogram, the spectral amplitudes along the entire distance are merged. High-frequency but low-amplitude noise components will cluster at the lower end of the histogram, forming spectral peaks, while low-frequency but high-amplitude flat phase components will be scattered at the higher end. It is important to note that the number of histogram groups should not be too small to ensure effective differentiation between flat phase components and residual phase components. When selecting the threshold, if T is too small, the selected amplitude sequence will contain terrain details and noise, causing interference; if T is too large, only the spectral peaks will be retained, and the entire algorithm will degenerate into a traditional frequency shifting method, resulting in flat phase loss. Considering that flat phases are mostly distributed near spectral peaks, adjacent frequency points have large amplitude differences, and there are few frequency points with the same amplitude level, they will appear as low-frequency convex points in the histogram. The position of this convex point in the histogram is the dividing point between the noise spectrum and the flat phase spectrum. Therefore, this invention selects the maximum point with the smallest corresponding amplitude in the histogram as the selection criterion for T.

[0074] To calculate the location of the maximum point in the histogram, the histogram sequence H needs to be subjected to a second difference:

[0075] D(k)=[H(k+1)-H(k)]-[H(k)-H(k-1)],k=2,3,...,N-1 (5)

[0076] Where D is the second-order difference sequence and N is the number of histogram groups, all the maximum points in the histogram can be found by judging that D(k)≥0.

[0077] It is important to note that due to the randomness of noise, when there are many histogram groups and the noise is severe, the maximum point of minimum amplitude may appear in the noise-dominated region on the left side of the histogram. If this point is chosen as the threshold, the estimated flat phase will contain a large amount of noise. To avoid this phenomenon, the suitability of this point as a threshold can be further determined based on its frequency level.

[0078] For the frequency level of noise, since it represents the main part of the entire histogram, it is usually... That is, frequencies above the uniformly distributed frequency distribution in the histogram correspond to frequencies above the uniformly distributed frequency distribution, while frequencies below the flat-phase frequency distribution correspond to frequencies below the uniformly distributed frequency distribution. In summary, only when the amplitude reaches its minimum maximum value does its frequency fall below or equal to the maximum frequency distribution. Only when the maximum value is satisfied should it be used as the threshold; otherwise, it should be further determined whether the remaining maximum points satisfy the threshold condition, and the point with the smallest amplitude should be used as the threshold. Thus, this invention can adaptively select a suitable threshold based on the number of histogram groups.

[0079] The normalized histogram constructed from the distance-direction spectral magnitude sequence of the simulated image is as follows: Figure 5 As shown, the spectral amplitude is mainly concentrated on the left side of the histogram and has a falling distribution, while the right side has a number of scattered convex points, which are the maximum points that need to be found.

[0080] Step four: For the portion of the spectral amplitude sequence below the threshold, take its mean as the background noise power estimate and set this portion to zero; for the portion above the threshold, subtract the background noise power estimate to further suppress the noise effect, obtaining the processed spectral amplitude sequence.

[0081] The processed spectral amplitude sequence is as follows:

[0082]

[0083] Where noise = mean(P(u) < T) is the estimated background noise value, and mean() represents the mean operation. It should be noted that since the noise is uniformly distributed across the entire spectrum, the portion of the noise with an amplitude greater than the threshold should also be subtracted to reduce the impact of the noise.

[0084] Step 5: Recombine the processed spectral amplitude sequence and spectral phase sequence, perform inverse discrete Fourier transform, and take the phase angle to obtain the nonlinear flat phase sequence. Expand it to the original size and subtract it from the original interferometric phase diagram to obtain the final flattened result.

[0085] The processed spectral amplitude sequence needs to be further restored to its flat-ground phase form. After recombining it with the phase sequence, the spectrum corresponding to the flat-ground phase is obtained. Performing an inverse discrete Fourier transform on it yields the corresponding range-oriented flat-ground complex phase form. Taking the principal phase value gives the flat-ground phase sequence.

[0086] Φ p =arg(IFFT(P′⊙Φ)) (7)

[0087] To obtain the flat-ground phase diagram corresponding to the original interferometric phase diagram, this sequence needs to be further extended. Since the flat-ground phase diagram appears as interference fringes parallel in the range direction, it can be multiplied by a matrix of all 1 columns:

[0088]

[0089] The final de-flattening result is obtained by subtracting the original interference phase diagram from the flat phase diagram and then wrapping them together.

[0090] right Figure 2 The flat-land effect is removed using the method proposed in this invention, and the resulting flat-land phase map is as follows: Figure 6 As shown, the result of leveling the ground is as follows Figure 7 As shown, and the distance spectrum after going to flat ground as shown Figure 8 As shown.

[0091] To verify the effectiveness of the nonlinear deflating method proposed in this invention, the same interferometric phase map was deflated using both the traditional frequency shift method and the block method. The results are as follows: Figures 9-11 and Figures 12-14 As shown.

[0092] The flat-ground phases obtained by the three methods are compared with the true flat-ground phases, and the Phase Standard Deviation (PSD) index is introduced for comparison. The PSD uses the true phase as a reference; the smaller the PSD value, the smaller the deviation between the estimated result and the true value. It can be calculated using the following formula:

[0093]

[0094] Where P and Q are the dimensions of the phase diagram. For the estimated flatland phase, For the true flat ground phase, W[·] is used to wrap the phase, taking the principal phase value of [-π,π].

[0095] The statistical results are shown in Table 2:

[0096] Table 2 Simulation Data Evaluation Results

[0097]

[0098] From Table 2 and Figures 6-14 The results of the three methods for leveling the ground shown can be seen as follows:

[0099] The flat ground fringes estimated using the traditional frequency shift method are uniformly distributed periodic fringes. After removing the flat ground, there are still obvious fringe changes on the left side of the residual phase, and the peak of the range spectrum is shifted to near zero frequency.

[0100] The block-based method estimates the flat-ground phase fringes with significant variations in thickness, resulting in a slower color change in the residual phase map compared to the traditional frequency-shifting method, and narrower spectral peaks.

[0101] Compared to the two existing methods, the method proposed in this invention estimates that the width of the flat-ground phase fringes increases slowly from left to right, the color change in the residual phase is more gradual, and the phase change caused by topographic details is only shown in local areas. The spectral peaks are closer to zero frequency than the previous two methods.

[0102] On the other hand, judging from the PSD value, the nonlinear flattening effect removal method proposed in this invention has the smallest PSD, proving that it best fits the phase of the real flattened ground.

[0103] Based on the above analysis, it can be seen that the InSAR nonlinear deflating effect method proposed in this invention is more adaptable to different scenarios, and the deflating results are significantly improved compared with existing methods.

[0104] The above description is merely a partial embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any equivalent substitutions and modifications made without departing from the spirit and principles of the present invention should be covered within the scope of the present invention.

Claims

1. A nonlinear method for removing flat-land effects in InSAR, characterized in that, It includes the following steps: Step 1: The synthetic aperture radar (SAR) performs two imaging operations on the same target to obtain two complex images, which are then registered. The two registered SAR complex images are then interferometric to obtain the corresponding interferometric phase diagram. Step 2: Pre-filter the complex interferometric image, and perform discrete Fourier transform on the range data of the pre-filtered complex interferometric image and sum them to obtain the spectral amplitude sequence and phase sequence corresponding to the range direction; The formula for complex interferometric phase decomposition of the pre-filtered complex interferometric image is as follows: (1) in, Noisy phase The plural form, The imaginary unit, The coordinates are the position. The number of vectors that make up the phase of the complex interference. For the first The weights of each vector, Represents the corresponding number The range frequency of each vector; Then, the complex interference phase sequence is aligned in the direction of each distance. Perform a discrete Fourier transform and sum the results to obtain the spectral amplitude sequence. The summation result yields the spectral amplitude sequence. and phase sequence : (2) in Representing the A complex interference phase sequence with distance upwards, The coordinates of the spectral sequence. Represents the Shur product operator; Step 3: Construct a normalized histogram based on the spectrum amplitude sequence, and adaptively select a threshold based on the histogram distribution; The process of adaptively selecting the threshold based on the histogram is as follows: Step 301: Construct a normalized histogram based on the spectral amplitude sequence, and process the normalized histogram sequence. Perform second-order difference operations and find all the local maxima. Step 302: Starting from the left end, sequentially determine the relationship between the frequencies corresponding to the maximum points and the frequencies under a uniform histogram distribution. If... This represents the frequency level of the noise; if This represents the frequency level of the flat-ground phase. The number of histogram groups; Step 303: Select the maximum point with the smallest amplitude in the corresponding histogram from the maximum points that satisfy the flat-ground phase frequency condition as the threshold. ; Step 4: For the portion of the spectral amplitude sequence below the threshold, take its mean as the estimated background noise power value and set this portion to zero; for the portion above the threshold, subtract the estimated background noise power value to further suppress the influence of noise. After the above operations, the processed spectral amplitude sequence is obtained. Step 5: Recombine the processed spectral amplitude sequence and spectral phase sequence, perform inverse discrete Fourier transform, and obtain the nonlinear flat-ground phase sequence by taking the phase angle, and then expand it to the original size; Step 6: Subtract the nonlinear flat phase sequence expanded to the original size from the interferometric phase diagram to obtain the final flattening result.

2. The InSAR nonlinear flat-land effect removal method according to claim 1, characterized in that, In step one, the noisy phase in the interference phase diagram is... Decomposed into flat phase Residual phase with changes in reaction details : .

3. The InSAR nonlinear flat-land effect removal method according to claim 1, characterized in that, In step three, the normalized histogram sequence is described. Perform second-order difference operations: (3) in It is a second-order difference sequence, determined by... This allows you to find all the maxima in the histogram.

4. The InSAR nonlinear deflating effect method according to claim 1, characterized in that, In step four, the estimated background noise power value is: , This indicates the operation of taking the average; The processed spectral amplitude sequence is as follows: (4)。 5. The InSAR nonlinear flat-land effect removal method according to claim 1, characterized in that, In step five, the nonlinear flat-ground phase sequence is as follows: (5) This is the processed spectral amplitude sequence; The formula for expanding a nonlinear flat-ground phase sequence to its original size is: (6)。

Images