Space-borne sliding bunching and TOPS mode interference measurement and DSM generation method

Abstract

An interference measurement and DSM generation method of a satellite-borne sliding spotlight and TOPS mode comprises the following steps: in the satellite-borne sliding spotlight and TOPS mode, carrying out imaging processing on original echo data to obtain an SAR (Synthetic Aperture Radar) image; in order to solve the problem of phase jump caused by Doppler center change, phase compensation and resampling are carried out on an SAR image, a non-common frequency spectrum is removed in combination with frequency domain pre-filtering, and the image registration precision and coherence are improved; in interferometric phase processing, filtering is carried out on an interferometric phase of a complex terrain region by combining a prior DEM and a non-neighborhood filtering algorithm, so that the quality of the interferometric phase is improved; for the problem of phase unwrapping in a high-resolution and large-breadth scene, a DCT unwrapping algorithm based on a hard threshold operator and a Kalman filtering unwrapping algorithm based on parameter optimization are adopted, and efficient and stable phase unwrapping is achieved. According to the method, the registration performance, coherence and phase processing precision of the InSAR image in the sliding bunching and TOPS mode can be effectively improved, and generation of the high-precision digital elevation model DSM in a complex scene is realized.

Description

Technical Field

[0001] This invention relates to a method for interferometry and DSM generation of spaceborne sliding beam focusing and TOPS mode, belonging to the field of SAR interferometry technology. Background Technology

[0002] Since its invention in the 1950s, Synthetic Aperture Radar (SAR) has been widely used in land surveying, environmental monitoring, and disaster early warning, thanks to its advantages of all-weather, all-day, and high-resolution Earth observation. With the development of platform and payload technologies, various systems have emerged, including airborne, spaceborne, and multi-satellite networking. Among these, spaceborne SAR has become a research focus due to its wide imaging swath and high coverage efficiency. In the 21st century, the imaging performance and application level of SAR have significantly improved, exhibiting a diversified development trend with higher resolution, wider swath, and multiple polarizations. Typical modes include strip, scanning, TOPS, spotting, and sliding spotting.

[0003] Building upon SAR imaging, Interferometric Synthetic Aperture Radar (InSAR) can acquire high-precision digital surface models (DSMs) through interferometric processing of primary and secondary images, and has been widely used in topographic mapping and deformation monitoring. As SAR technology advances towards higher resolution and wider swaths, InSAR processing is gradually expanding to complex scenarios, such as urban areas and non-uniform terrain, possessing the potential to acquire more detailed topographic information. However, InSAR processing faces new challenges under new imaging modes. Especially in sliding beamforming and TOPS modes, the interferometric phase is easily affected by changes in the Doppler center due to beam pointing rotation and spectral aliasing; simultaneously, its orbit and imaging geometry differ from conventional modes, making it difficult to guarantee registration and positioning accuracy. Furthermore, the large volume of interferometric data under high resolution and wide swath conditions significantly increases the complexity of phase filtering and unwrapping, necessitating optimization of processing efficiency and accuracy.

[0004] Although some research, both domestically and internationally, has focused on TerraSAR-X, COSMO-SkyMed, and Sentinel satellite data, exploring interferometric processing in sliding spotting and TOPS modes, most studies still concentrate on single-module analysis or specific platform applications, with relatively insufficient systematic, end-to-end research. Currently, achieving stable and high-precision InSAR processing for sliding spotting and TOPS modes, and developing a complete, engineering-applicable solution, has become a pressing issue in this field. Summary of the Invention

[0005] The technical problem solved by this invention is to overcome the shortcomings of the prior art and provide an interferometric measurement and DSM generation method for spaceborne sliding spotting and TOPS modes, which realizes stable and high-precision InSAR processing for sliding spotting and TOPS modes. The technical solution of this invention is: A method for interferometry and DSM generation of spaceborne sliding beam focusing and TOPS mode, comprising: Step 1: Acquire raw echo data of spaceborne SAR in sliding spotting and TOPS modes, and determine the corresponding imaging parameters; Step 2: Based on the imaging principle of sliding spotting and TOPS mode, the acquired raw echo data is processed to obtain SAR images. Step 3: Perform phase compensation and resampling for phase jumps caused by changes in the Doppler center of SAR images; Step 4: Perform frequency domain pre-filtering on the resampled signal to remove non-common spectrum regions of the main and auxiliary images; Step 5: For the main and auxiliary images after frequency domain pre-filtering, use a non-neighborhood filtering method based on prior DEM to filter the interferometric phase and reduce noise interference in complex terrain areas. Step 6: For the interferometric phase-filtered signal, phase unwrapping is performed using a Kalman filter phase unwrapping method based on adaptive parameter optimization; Step 7: After completing phase unwrapping, generate a digital surface model (DSM) based on the interference phase.

[0006] Furthermore, step 3 involves phase compensation for the phase jump caused by changes in the Doppler center of the SAR image, specifically as follows: Step 3.1: Calculate the Doppler center variation history of each pixel in the SAR image. Doppler center slope of zero Doppler azimuth line :

[0007]

[0008] in, Indicates the satellite payload speed. For wavelength, This refers to the system beam angular oscillation velocity; For the equivalent velocity of the satellite payload, This is the slant distance to the perigee. Step 3.2: Calculate the Doppler center slope of the zero Doppler azimuth before imaging:

[0009] Inverse phase compensation is applied to the focused SAR image in the azimuth time domain to achieve phase-preserving processing of the image data:

[0010] in, For time-domain phase compensation, This refers to the direction of time.

[0011] Further, resampling after phase compensation is performed, specifically as follows: The sinc function is used as the interpolation kernel to perform pixel-by-pixel interpolation and resampling on the phase-compensated non-baseband signal, so that the auxiliary image corresponds to the main image in space. After resampling, phase inverse compensation is performed on the data based on the offset of each pixel to restore the original phase information of the auxiliary image. The compensated phase is as follows:

[0012] in, for The azimuth offset corresponding to the position, after phase inverse compensation, will yield the resampled auxiliary image.

[0013] Furthermore, step 4 involves frequency domain pre-filtering of the resampled signal, specifically as follows: After the master and slave images are registered, phase compensation is performed on the master and slave images respectively based on the Doppler center change to achieve spectral deskewing; the phase-compensated master and slave images are transformed to the frequency domain, and the frequency domain filtering range is determined according to the geometric offset of the master and slave images and the slope of the Doppler center change to filter out non-common frequency spectrum regions; phase inverse compensation is performed on the filtered image data to output the filtered master and slave images.

[0014] Furthermore, the prior DEM includes SRTM and TanDEM elevation models, which are used to provide terrain reference information during interferometric phase filtering and unwrapping processes.

[0015] Furthermore, step 5, for the pre-filtered primary and secondary images in the frequency domain, uses a non-neighborhood filtering method based on the prior DEM to filter the interferometric phase, specifically as follows: Step 5.1: Invert the terrain based on the prior DEM data and calculate the original terrain phase. Multiply the original interferometric phase by the conjugate of the terrain phase to remove the terrain phase. Step 5.2: Based on noisy observations and , , Let each represent the set of observed data within a pixel block. The initial likelihood values ​​between the pixel blocks are calculated using the following scale-invariant similarity model:

[0016] in: , and This represents the amplitude value of two pixels in the first image. and This represents the amplitude value of the corresponding pixel in the second image; , Represents a set of parameters; , and Represents the phase value of the corresponding pixel; intermediate variable ; And further based on the parameter estimates obtained from the previous iteration and Where i represents the iteration number, the weights are refined using the symmetric Kohlbek-Leibler divergence SD_KL. The formula for calculating SD_KL is:

[0017] in Indicates that under given parameters The probability density function of the observed value O under given conditions; Update the weight accordingly. It is in logarithmic form, and the values ​​of accumulators a, x, and N are refreshed, where a is the weighted sum of squares of magnitude, x is the weighted sum of cross-correlation, and N is the weighted sum; This represents the similarity weight between pixels s and t; Step 5.3: After smoothing and denoising all pixels, calculate the reflectance of pixel s based on the weighted maximum likelihood estimation method. Interference phase With coherence coefficient :

[0018] in, This represents the weighted sum of squares; , These are the main image data and the auxiliary image data, respectively. Indicates the weighted cross-correlation sum. for The complex conjugate; The sum of weights is represented by t; t represents the pixel index within the search window. Step 5.4: Calculate the parameter estimates obtained in this iteration. As the initial value for the next iteration, steps 5.2 to 5.3 are executed repeatedly until the parameters converge or the preset number of iterations is reached, and intermediate filtering results are obtained. Step 5.5: Restore the terrain phase. Restore the terrain phase that was initially removed by the algorithm into the filtering result to obtain the final interferometric phase filtering result.

[0019] Furthermore, step 6, for the interferometric phase-filtered signal, employs a Kalman filter phase unwrapping method based on adaptive parameter optimization for phase unwrapping, specifically as follows: Step 6.1: Select the master and slave image data corresponding to the phase to be untangled, and define the Kalman filter parameters, including the initial error covariance matrix. Process noise covariance matrix Covariance matrix ; Step 6.2: Calculate the estimated phase slope values ​​for the column and row directions respectively using the following formula:

[0020] in, This represents the estimated phase slope in the column direction. This represents the estimated phase slope in the row direction. Represents an image. This is the index in the distance direction. This is an index for the direction of orientation; Calculate the absolute difference between the predicted phase slope and the wrapped phase gradient. If the difference is greater than a set threshold, optimize the slope as follows:

[0021] in, Indicates the column difference. Indicates the row difference; Step 6.3: Combine the error covariance matrix With coherence coefficient matrix Calculate the adaptive weights, including: Weights are set based on covariance:

[0022] in, This represents the weight at position m. This represents the weight at position n; Weights are adjusted based on the coherence coefficient:

[0023] in, This represents the weight of the coherence coefficient at row m-1 and column n; This represents the weight of the coherence coefficient at row m, column n-1; The final prediction weights are:

[0024] , These represent the final predicted weights at indices m and n, respectively. Step 6.4: Use self-final prediction weights to perform state prediction and covariance prediction:

[0025]

[0026] For the phase prior prediction at pixel (m, n), For the phase posterior prediction at pixel (m, n), Predicting the phase prior covariance at pixel (m, n), For the prediction of the phase posterior covariance at pixel (m, n); Step 6.5: Calculate Kalman gain and update state based on the observations, traverse the entire image until all pixels have been processed, and output the final unwrapped phase result.

[0027] Furthermore, a digital surface model (DSM) is generated based on the unwrapped phase results, specifically as follows: Step 8.1: Based on prior DEM data, calculate the absolute slant range difference of the processed area through forward simulation, and then invert the absolute ambiguity number accordingly. ; Step 8.2: Utilize the absolute fuzzy number Untangling phase Convert to absolute phase The transformation relationship is as follows: ; Step 8.3, based on absolute phase The slope distance difference is calculated, and the target elevation value is obtained through geometric transformation. After being projected onto the ground plane, the final DSM is generated.

[0028] In a second aspect, the present invention also proposes a non-volatile storage medium comprising: a computer program product, wherein the method is executed when the computer program product is executed.

[0029] Thirdly, the present invention also proposes a computer program product, which includes a computer program that, when executed by a processor, implements the method described.

[0030] The beneficial effects of this invention compared to the prior art are: This invention establishes an imaging and interferometric processing architecture suitable for sliding spotting and TOPS modes. By combining DEM-based phase compensation, nonlocal iterative filtering, and adaptive Kalman phase unwrapping algorithm, it overcomes the registration error and phase jump problem caused by Doppler time-varying and spectral aliasing. It also significantly suppresses unwrapping errors caused by noise interference and terrain undulation, effectively improving the accuracy and robustness of high-resolution, wide-swath interferometric data processing. Ultimately, it achieves high-precision and high-efficiency InSAR processing suitable for spaceborne sliding spotting and TOPS modes. Attached Figure Description

[0031] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is a schematic diagram of the geometry for spaceborne sliding beam focusing mode imaging provided by the present invention; Figure 3 This is a schematic diagram of the phase compensation and resampling process provided by the present invention; Figure 4 a is a diagram of the coherence coefficients after processing by the phase compensation registration algorithm provided by this invention; Figure 4 b is the phase image after processing by the phase compensation registration algorithm provided by this invention; Figure 5 a is a diagram of the filtered coherence coefficients provided by this invention; Figure 5 b is a comparison diagram of the distribution of the filtered coherence coefficients provided by the present invention; Figure 6 a is used to verify the original SAR image using the filtering algorithm provided by this invention; Figure 6 b is the pre-filtering interference phase diagram provided by this invention; Figure 7 a represents the phase filtering result obtained by the mean filtering algorithm provided by this invention; Figure 7 b represents the phase filtering result obtained by the non-neighborhood filtering algorithm provided by this invention; Figure 7 c is a schematic diagram of the phase filtering result obtained by the mean filtering algorithm based on prior DEM provided by the present invention; Figure 7 d is a schematic diagram of the phase filtering result obtained by the non-neighborhood filtering algorithm based on prior DEM provided by the present invention; Figure 8 a is the unwound phase diagram in the simulated phase diagram provided by the present invention; Figure 8 b is the wrapped phase diagram in the simulated phase diagram provided by the present invention; Figure 9 a is the unwrapping phase diagram of the POEKFPU unwrapping result provided by the present invention; Figure 9 b is a phase error diagram of the POEKFPU unwrapping result provided by the present invention; Figure 9 c is the rewinding phase diagram of the POEKFPU unwinding result provided by the present invention; Figure 10 a is the DSM result diagram generated by interferometry processing of the measured data provided by this invention; Figure 10 b is a DSM elevation accuracy analysis diagram based on the measured data interferometry processing provided by the present invention; Figure 11 a is a reference DEM diagram provided by the present invention; Figure 11 b is the generated DSM result diagram provided by the present invention. Detailed Implementation

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

[0033] Currently, existing spaceborne sliding beamforming and TOPS mode interferometric SAR data processing methods, under the requirements of high-resolution and wide-swath applications, suffer from severe noise in the interferometric phase due to the strong time-varying nature of the Doppler center, significant spectral aliasing, and substantial motion phase errors. This leads to decreased image registration accuracy and increased unwrapping phase jumps, seriously affecting the accuracy of digital elevation model reconstruction. Traditional phase filtering and unwrapping algorithms are sensitive to complex terrain and noise, and have low computational efficiency, making it difficult to meet the requirements of processing accuracy and timeliness in practical engineering.

[0034] This invention provides a method for processing Interferometric Synthetic Aperture Radar (InSAR) data and generating Digital Elevation Models (DSMs) in spaceborne sliding spotting and TOPS modes. The method includes: imaging the raw echo data to obtain SAR images in spaceborne sliding spotting and TOPS modes; addressing phase jump issues caused by Doppler center variations by performing phase compensation and resampling on the SAR images, and combining frequency domain pre-filtering to remove non-common spectra, thereby improving image registration accuracy and coherence; in interferometric phase processing, combining prior DEM and non-neighborhood filtering algorithms to filter the interferometric phase in complex terrain areas, improving interferometric phase quality; and addressing phase unwrapping issues in high-resolution and wide-swath scenarios by employing a DCT unwrapping algorithm based on a hard threshold operator and a parameter-optimized Kalman filter unwrapping algorithm to achieve efficient and stable phase unwrapping. This invention effectively improves the registration performance, coherence, and phase processing accuracy of InSAR images in sliding spotting and TOPS modes, enabling the generation of high-precision Digital Elevation Models (DSMs) in complex scenarios.

[0035] Please see Figure 1 , Figure 1This is a schematic diagram of a phase processing flow provided in an embodiment of the present invention. The method includes: Step 1: Acquire raw echo data of spaceborne SAR in sliding spotting and TOPS modes, and determine the corresponding imaging parameters.

[0036] Specifically, the system acquires raw echo data from a spaceborne SAR system in either sliding spotting mode or TOPS mode, and determines the corresponding system imaging parameters. These imaging parameters include platform parameters, imaging geometric parameters, mode-specific parameters, and signal parameters. The platform parameters include satellite platform position, velocity vector, imaging time window, and pulse repetition frequency. The imaging geometric parameters include slant range, incident angle, azimuth beam pointing angle, and their time-varying sequence. The mode-specific parameters include beam scanning rate and scanning range in sliding spotting mode, and sub-beam switching interval and wavelet arrangement sequence in TOPS mode. The signal parameters include carrier frequency, bandwidth, pulse width, frequency modulation slope, and Doppler parameters.

[0037] Step 2: Based on the imaging principles of sliding spotting and TOPS mode, a unified imaging algorithm is used to process the raw echo data to obtain SAR images.

[0038] Specifically, attached Figure 2 This is a schematic diagram of the imaging geometry for spaceborne sliding beam focusing mode. The acquired echo data undergoes imaging processing, specifically including the following sub-steps: Step 2.1, Azimuth Preprocessing: Based on Beam Scanning Speed With platform speed The relationship between the azimuth and frequency domains is used to perform dealiasing processing on the azimuth signal; range-frequency domain-azimuth time domain echo signal:

[0039] in, For distance frequency, For location and time, Indicates the illumination time at the beam center. As weight, This indicates the horizontal offset of the target in the azimuth direction (along the radar's direction of movement) relative to the radar's crossing point. Indicates the satellite's speed along its trajectory; Indicates signal bandwidth. Represents the distance spectrum signal envelope. This represents the azimuth window function. For the time to synthesize the pore size, For radar carrier frequency, The instantaneous slant range at which the radar reaches the target.

[0040] Directional reference signal: , in, The wavelength of electromagnetic waves, This is the reference slope distance for the sliding contact.

[0041] echo signal With azimuth reference signal Convolution is equivalent to multiplication in the frequency domain, resulting in a preprocessed signal that eliminates azimuth spectral aliasing. Step 2.2, Range Processing: The CS (Chirp Scaling) algorithm is used to perform chirpscaling in the range-Doppler domain to achieve azimuth scaling. Then, range compression is achieved through range-direction matched filtering, and range migration is corrected. The signal after this step is represented as follows:

[0042] in, For azimuth frequency, The center frequency of the Doppler wave. Zero Doppler slant distance, The maximum Doppler frequency.

[0043] Step 2.3, Azimuth Focusing: Using the phase compensation function:

[0044] The hyperbolic form of the azimuth signal is transformed into a quadratic form, and then the azimuth IFFT and frequency domain compensation function are used to finally complete the azimuth matched filtering and achieve signal focusing. Frequency domain compensation function:

[0045] in, This refers to the time of the transformed azimuth.

[0046] Step 2.4, Interference Phase Restoration: Correction is performed on the phase error introduced during the imaging process. The correction function is: To recover the original phase information that can be used for interferometric processing.

[0047] Step 3: Perform phase compensation for phase jumps caused by changes in the Doppler center of the image.

[0048] Specifically, it includes the following sub-steps. Step 3.1: Calculate the Doppler center change history of each pixel. Doppler center slope of zero Doppler azimuth line :

[0049]

[0050] in, Indicates the satellite payload speed. For wavelength, This refers to the system beam angular oscillation speed. For the equivalent velocity of the satellite payload, This is the slant distance to the perigee.

[0051] Step 3.2: Calculate the Doppler center slope of the zero Doppler azimuth before imaging:

[0052] Inverse phase compensation is applied to the focused SAR image in the azimuth time domain to achieve phase-preserving processing of the image data:

[0053] in, For time-domain phase compensation, This refers to the direction of time.

[0054] Step 4, resampling after phase compensation includes the following steps; please refer to the appendix for the resampling procedure. Figure 3 : The sinc function is used as the interpolation kernel to perform pixel-by-pixel interpolation and resampling on the phase-compensated non-baseband signal, so that the auxiliary image corresponds to the main image in space.

[0055] After resampling, phase inverse compensation is performed on the data based on the offset of each pixel to restore the original phase information of the auxiliary image. The compensated phase is as follows:

[0056] in, for The azimuth offset corresponding to the position. After phase inverse compensation, the resampled auxiliary image can be obtained. Please refer to the appendix for the coherence coefficient map and phase map after phase compensation and registration of the simulation data. Figure 4 (a), 4 (b).

[0057] Step 5, frequency domain pre-filtering of the resampled signal includes the following sub-steps: Step 5.1: After the primary and secondary images are registered, phase compensation is performed on the primary and secondary SAR images respectively based on the changes in the Doppler center to achieve spectral deskewing; Step 5.2: Perform Fourier transform on the phase-compensated primary and secondary SAR images to obtain the frequency domain signal; Step 5.3: Determine the frequency domain filtering range based on the geometric offset of the main and auxiliary images and the slope of the Doppler center change, and filter out non-common spectrum regions. Step 5.4: Perform phase inverse compensation on the restored image data to obtain the filtered main and auxiliary images.

[0058] Please refer to the appendix for a comparison of the coherence coefficient plot and coherence coefficient distribution after frequency domain pre-filtering of the simulation data. Figure 5 (a), 5 (b).

[0059] Step 6: For the filtered master and slave images, a non-neighborhood phase filtering method based on the prior DEM is proposed, specifically including the following sub-steps: Step 6.1: Invert the terrain based on the prior DEM data and calculate the original terrain phase. Multiply the original interferometric phase by the conjugate of the terrain phase to remove the terrain phase. Step 6.2: Based on noisy observations and ( , Let each represent the set of observed data within a two-pixel block. The initial likelihood values ​​between the pixel blocks are calculated using the following scale-invariant similarity model:

[0060] in: , and This represents the amplitude value of two pixels in the first image. and This represents the amplitude value of the corresponding pixel in the second image; , and Represents the phase value of the corresponding pixel; intermediate variable ; And further based on the parameter estimates obtained from the previous iteration and (where i represents the iteration number, (representing the parameter set), the weights are refined using the symmetric Kohlbek-Leibler divergence (SD_KL). The formula for calculating SD_KL is:

[0061] in Indicates that under given parameters The probability density function of the observed value O under given conditions; Update the weight accordingly. (Represents the similarity weight between pixels s and t) and its logarithmic form, and refreshes the values ​​of accumulators a, x, and N, where a is the weighted sum of squared magnitudes, x is the weighted sum of cross-correlation, and N is the weighted sum; Step 6.3: After smoothing and denoising all pixels, calculate the reflectance of pixel s based on the weighted maximum likelihood estimation method. Interference phase With coherence coefficient :

[0062] in, , representing the weighted sum of squares; , indicating the weighted cross-correlation sum, for The complex conjugate; , represents the weight sum; t represents the pixel index within the search window; Step 6.4: Estimate the parameters obtained in this iteration. As the initial value for the next iteration, steps 6.2 to 6.3 are executed repeatedly until the parameters converge or the preset number of iterations is reached, and intermediate filtering results are obtained. Step 6.5: Restore the terrain phase. Restore the terrain phase that was initially removed by the algorithm into the filtering result to obtain the final interferometric phase filtering result.

[0063] To analyze the performance and efficiency of filtering algorithms in high-resolution complex scenes, measured data from a certain urban area were selected to generate interferometric phase images. Please refer to the appendix for the original SAR image and the interferometric phase map. Figure 6 .

[0064] The generated interference phase was filtered using different phase filtering algorithms. Please refer to the appendix for the phase filtering results. Figure 7 .

[0065] Table 1 shows the filtering performance of each algorithm, using four evaluation criteria: number of residual points, root mean square error, phase gradient sum, and filtering time. The effectiveness of the algorithms is verified based on these four criteria.

[0066] Table 1

[0067] The noise suppression performance is analyzed below based on the number of residual points, phase gradient, and angle: Due to the complexity of urban features, phase information is also more complex. The non-neighborhood phase filtering algorithm based on prior DEM calculates weights based on image amplitude and phase information during the filtering process and performs mean filtering, avoiding the destruction of phase space continuity. This algorithm has the highest performance compared to other algorithms.

[0068] The root mean square error represents the magnitude of the phase error before and after filtering by the phase filtering algorithm. The non-neighborhood phase filtering algorithm based on prior DEM first calculates the similarity between each pixel during the filtering process, and then calculates weights based on this similarity for filtering. This results in smaller filtering errors for data in complex urban scenes. In contrast, the two mean filtering algorithms show significantly higher filtering errors in complex areas compared to the non-neighborhood algorithm. In terms of filtering efficiency, the mean filtering algorithm is simpler and significantly more efficient than the non-neighborhood filtering algorithm. Non-neighborhood filtering uses an iterative approach to optimize filtering parameters, leading to a significant performance reduction. Introducing prior DEM information reduces scene complexity, and the efficiency of the prior DEM-based non-neighborhood filtering algorithm is significantly better than that of the non-neighborhood filtering algorithm.

[0069] Analysis shows that the proposed non-neighborhood filtering algorithm based on prior DEM has good filtering performance for complex urban scene data, and the efficiency problem of the non-neighborhood filtering algorithm has been improved.

[0070] Step 7: For interferometric phase-filtered signals, an adaptive parameter-optimized Kalman filter phase unwrapping method is proposed, specifically including the following sub-steps: Step 7.1: Select the primary and secondary SAR image data corresponding to the phase to be unwrapped, and define the Kalman filter parameters, including the initial error covariance matrix. Process noise covariance matrix covariance matrix ; Step 7.2: Calculate the estimated phase slope values ​​for the column and row directions respectively using the following formula:

[0071] in, This represents the estimated phase slope in the column direction. This represents the estimated phase slope in the row direction. Represents an image. This is the index in the distance direction. This is the index of the azimuth direction.

[0072] Furthermore, the absolute difference between the predicted phase slope and the winding phase gradient is calculated. If the difference is greater than a set threshold, the slope is optimized as follows:

[0073] in, Indicates the column difference. Indicates the row difference; Step 7.3: Combine the error covariance matrix With coherence coefficient matrix Calculate the adaptive weights, including: Weights are set based on covariance:

[0074] in, This represents the weight at position m. This represents the weight at position n; Weights are adjusted based on the coherence coefficient:

[0075] in, This represents the weight of the coherence coefficient at row m-1 and column n; This represents the weight of the coherence coefficient at row m, column n-1; The final prediction weights are:

[0076] , These represent the final predicted weights at indices m and n, respectively. Step 7.4: Use the final prediction weights to perform state prediction and covariance prediction:

[0077]

[0078] For the phase prior prediction at pixel (m, n), For the phase posterior prediction at pixel (m, n), Predicting the phase prior covariance at pixel (m, n), For the prediction of the phase posterior covariance at pixel (m, n); Step 7.5: Calculate Kalman gain and update state based on the observations, traverse the entire image until all pixels have been processed, and output the final unwrapped phase result.

[0079] The `peaks` function in MATLAB is used to generate phase information of size 512×512, with the phase value range set to [value range missing]. Add 4.0dB of noise and wrap it around. Please see the appendix for the results. Figure 8 .

[0080] For the appendix Figure 8 The wrapped phase data shown is processed using the proposed adaptive parameter-optimized Kalman filter phase unwrapping algorithm. The unwrapped phase is subtracted from the original phase to obtain the unwrapped phase error map. Then, the unwrapped result is rewrapped. The results are shown in the appendix. Figure 9 .

[0081] Step 8: Determine the absolute phase of the unwrapped phase and generate a digital surface model (DSM), specifically including the following sub-steps.

[0082] Step 8.1: Based on prior DEM data, calculate the absolute slant range difference of the processed area through forward simulation, and then invert the absolute ambiguity number accordingly. ; Step 8.2: Utilize the absolute fuzzy number Untangling phase Convert to absolute phase The transformation relationship is as follows: ; Step 8.3: Calculate the slant range difference based on the absolute phase, then convert it through geometric relationships to obtain the target elevation value, and finally generate the DSM after projection onto the ground plane.

[0083] According to step 8, the difference between the generated DSM and the reference DEM is calculated. Analysis of the elevation difference shows that the difference between the measured data and the DEM is generally within 2.5m, verifying the correctness of the proposed interferometric processing workflow under the sliding beam-gathering mode. Please refer to the appendix for the interferometric processing results of the measured data. Figure 10 .

[0084] Please refer to the appendix for a comparison between the reference DEM and the processed DSM. Figure 11 As can be seen, in the sliding beamforming mode, the resolution is improved, and the resulting DSM is also more refined.

[0085] In summary, the sliding beam-gap interferometric SAR data processing method proposed in this embodiment of the invention can effectively improve phase quality and unwrapping accuracy while maintaining high computational efficiency, thus proving that the adaptive interferometric phase processing method based on prior DEM and nonlocal filtering has the advantages of high accuracy and high robustness in complex scenarios.

[0086] This invention establishes an imaging and interferometric processing architecture suitable for sliding spotting and TOPS modes. By combining DEM-based phase compensation, nonlocal iterative filtering, and adaptive Kalman phase unwrapping algorithms, it overcomes registration errors and phase jump problems caused by Doppler time-varying and spectral aliasing. It also significantly suppresses unwrapping errors caused by noise interference and terrain undulations, effectively improving the accuracy and robustness of high-resolution, wide-swath interferometric data processing. Ultimately, it achieves high-precision and high-efficiency InSAR processing suitable for spaceborne sliding spotting and TOPS modes.

[0087] The parts of this invention not described in detail are common knowledge to those skilled in the art.

Claims

1. A method for interferometry and DSM generation of spaceborne sliding beam focusing and TOPS mode, characterized in that... include: Step 1: Acquire raw echo data of spaceborne SAR in sliding spotting and TOPS modes, and determine the corresponding imaging parameters; Step 2: Based on the imaging principle of sliding spotting and TOPS mode, the acquired raw echo data is processed to obtain SAR images. Step 3: Perform phase compensation and resampling for phase jumps caused by changes in the Doppler center of SAR images; Step 4: Perform frequency domain pre-filtering on the resampled signal to remove non-common spectrum regions of the main and auxiliary images; Step 5: For the main and auxiliary images after frequency domain pre-filtering, use a non-neighborhood filtering method based on prior DEM to filter the interferometric phase and reduce noise interference in complex terrain areas. Step 6: For the interferometric phase-filtered signal, phase unwrapping is performed using a Kalman filter phase unwrapping method based on adaptive parameter optimization; Step 7: After completing phase unwrapping, generate a digital surface model (DSM) based on the interference phase.

2. The InSAR processing method based on sliding spotting and TOPS mode according to claim 1, characterized in that: Step 3 addresses the phase jump caused by changes in the Doppler center of the SAR image by performing phase compensation, specifically as follows: Step 3.1: Calculate the Doppler center variation history of each pixel in the SAR image. Doppler center slope of zero Doppler azimuth line : in, Indicates the satellite payload speed. For wavelength, This refers to the system beam angular oscillation velocity; For the equivalent velocity of the satellite payload, This is the slant distance from the perigee. Step 3.2: Calculate the Doppler center slope of the zero Doppler azimuth before imaging: Inverse phase compensation is applied to the focused SAR image in the azimuth time domain to achieve phase-preserving processing of the image data: in, For time-domain phase compensation, This refers to the direction of time.

3. The InSAR processing method based on sliding spotting and TOPS mode according to claim 2, characterized in that: Resampling after phase compensation is as follows: The sinc function is used as the interpolation kernel to perform pixel-by-pixel interpolation and resampling on the phase-compensated non-baseband signal, so that the auxiliary image corresponds to the main image in space. After resampling, phase inverse compensation is performed on the data based on the offset of each pixel to restore the original phase information of the auxiliary image. The compensated phase is as follows: in, for The azimuth offset corresponding to the position, after phase inverse compensation, will yield the resampled auxiliary image.

4. The InSAR processing method based on sliding spotting and TOPS mode according to claim 1, characterized in that: Step 4 involves frequency domain pre-filtering of the resampled signal, specifically as follows: After the master and slave images are registered, phase compensation is performed on the master and slave images respectively based on the Doppler center change to achieve spectral deskewing; the phase-compensated master and slave images are transformed to the frequency domain, and the frequency domain filtering range is determined according to the geometric offset of the master and slave images and the slope of the Doppler center change to filter out non-common frequency spectrum regions; phase inverse compensation is performed on the filtered image data to output the filtered master and slave images.

5. The InSAR processing method based on sliding spotting and TOPS mode according to claim 1, characterized in that: The prior DEM includes SRTM and TanDEM elevation models, which are used to provide terrain reference information during interferometric phase filtering and unwrapping processes.

6. The InSAR processing method based on sliding spotting and TOPS mode according to claim 5, characterized in that: Step 5, for the pre-filtered primary and secondary images in the frequency domain, uses a non-neighborhood filtering method based on the prior DEM to filter the interferometric phase, specifically as follows: Step 5.1: Invert the terrain based on the prior DEM data and calculate the original terrain phase. Multiply the original interferometric phase by the conjugate of the terrain phase to remove the terrain phase. Step 5.2: Based on noisy observations and , , Let each represent the set of observed data within a pixel block. The initial likelihood values ​​between the pixel blocks are calculated using the following scale-invariant similarity model: in: , and This represents the amplitude value of two pixels in the first image. and This represents the amplitude value of the corresponding pixel in the second image; , Represents a set of parameters; , and Represents the phase value of the corresponding pixel; intermediate variable ; And further based on the parameter estimates obtained from the previous iteration and Where i represents the iteration number, the weights are refined using the symmetric Kohlbek-Leibler divergence SD_KL. The formula for calculating SD_KL is: in Indicates that under given parameters The probability density function of the observed value O under given conditions; Update the weight accordingly. It is in logarithmic form, and the values ​​of accumulators a, x, and N are refreshed, where a is the weighted sum of squares of magnitude, x is the weighted sum of cross-correlation, and N is the weighted sum; This represents the similarity weight between pixels s and t; Step 5.3: After smoothing and denoising all pixels, calculate the reflectance of pixel s based on the weighted maximum likelihood estimation method. Interference phase With coherence coefficient : in, This represents the weighted sum of squares; , These are the main image data and the auxiliary image data, respectively. This represents the weighted cross-correlation sum. for The complex conjugate; The sum of weights is represented by t; t represents the pixel index within the search window. Step 5.4: Calculate the parameter estimates obtained in this iteration. As the initial value for the next iteration, steps 5.2 to 5.3 are executed repeatedly until the parameters converge or the preset number of iterations is reached, and intermediate filtering results are obtained. Step 5.5: Restore the terrain phase. Restore the terrain phase that was initially removed by the algorithm into the filtering result to obtain the final interferometric phase filtering result.

7. The InSAR processing method based on sliding spotting and TOPS mode according to claim 1, characterized in that: Step 6, for the interferometric phase-filtered signal, employs a Kalman filter phase unwrapping method based on adaptive parameter optimization for phase unwrapping, specifically as follows: Step 6.1: Select the master and slave image data corresponding to the phase to be untangled, and define the Kalman filter parameters, including the initial error covariance matrix. Process noise covariance matrix Covariance matrix ; Step 6.2: Calculate the estimated phase slope values ​​for the column and row directions respectively using the following formula: in, This represents the estimated phase slope in the column direction. This represents the estimated phase slope in the row direction. Represents an image. This is the index in the distance direction. This is an index for the direction of orientation; Calculate the absolute difference between the predicted phase slope and the wrapped phase gradient. If the difference is greater than a set threshold, optimize the slope as follows: in, Indicates the column difference. Indicates the row difference; Step 6.3: Combine the error covariance matrix With coherence coefficient matrix Calculate the adaptive weights, including: Weights are set based on covariance: in, This represents the weight at position m. This represents the weight at position n; Weights are adjusted based on the coherence coefficient: in, This represents the weight of the coherence coefficient at row m-1 and column n; This represents the weight of the coherence coefficient at row m, column n-1; The final prediction weights are: , These represent the final predicted weights at indices m and n, respectively. Step 6.4: Use self-final prediction weights to perform state prediction and covariance prediction: For the phase prior prediction at pixel (m, n), For the phase posterior prediction at pixel (m, n), Predicting the phase prior covariance at pixel (m, n), For the prediction of the phase posterior covariance at pixel (m, n); Step 6.5: Calculate Kalman gain and update state based on the observations, traverse the entire image until all pixels have been processed, and output the final unwrapped phase result.

8. The InSAR processing method based on sliding spotting and TOPS mode according to claim 7, characterized in that: The Digital Surface Model (DSM) is generated based on the unwrapped phase results, specifically as follows: Step 8.1: Based on prior DEM data, calculate the absolute slant range difference of the processed area through forward simulation, and then invert the absolute ambiguity number accordingly. ; Step 8.2: Utilize the absolute fuzzy number Untangling phase Convert to absolute phase The transformation relationship is as follows: ; Step 8.3, based on absolute phase The slope distance difference is calculated, and the target elevation value is obtained through geometric transformation. After being projected onto the ground plane, the final DSM is generated.

9. A non-volatile storage medium, characterized in that, include: A computer program product that, when executed, performs the method described in any one of claims 1 to 8.

10. A computer program product, characterized in that, The computer program product includes a computer program that, when executed by a processor, implements the steps of the method according to any one of claims 1 to 8.

Images