A navigation satellite bistatic InSAR high-coherence DEM imaging method based on adaptive region segmentation

Abstract

This invention discloses a high-coherence DEM imaging method for bistatic InSAR based on adaptive region segmentation, belonging to the field of bistatic synthetic aperture radar technology. The specific implementation of this method includes: acquiring a coarse imaging result image and selecting CRC points of the coarse imaging result image; determining the pixel coordinates of each pixel corresponding to the image valley between different resolution units based on the coordinates of each CRC point, obtaining the image valley boundary of each resolution unit; performing boundary gap compensation on the image valley boundaries to obtain the adaptive region segmentation result corresponding to the boundary of each sub-region, the adaptive region segmentation result including multiple sub-regions; performing elevation fitting in each sub-region to generate a SAR image. This implementation can realize a high-coherence imaging algorithm based on adaptive partitioning and elevation fitting suitable for GNSS-based InSAR systems, improving elevation fitting accuracy, ensuring the integrity of resolution units and the continuity of elevation data, and achieving high target deformation inversion accuracy.

Description

Technical Field

[0001] This invention belongs to the field of bistatic synthetic aperture radar technology, specifically relating to a navigation satellite bistatic InSAR high coherence DEM imaging method based on adaptive region segmentation. Background Technology

[0002] A GNSS-based bistatic synthetic aperture radar (BiSAR) utilizes in-orbit navigation satellites as external radiation sources (i.e., transmitters) to transmit signals, with the receiver located near the ground. The receiver can be airborne, vehicle-mounted, or fixed. Due to the abundance of in-orbit navigation satellite resources, the receiver can receive navigation signals from at least eight satellites at any given time, achieving three-dimensional deformation inversion through information fusion from different angles.

[0003] GNSS-based BiSAR uses interferometric techniques to measure target elevation. In GNSS-based InSAR (Interferometric Synthetic Aperture Radar) systems, the accuracy of the Digital Elevation Model (DEM) information limits the accuracy of deformation inversion (i.e., the accuracy of inverting surface deformation based on the phase information of the image). On the one hand, the topological phase error caused by elevation error will greatly reduce the correlation between image pairs, thus affecting the interferometric measurement results. At the same time, elevation error will cause the same target to have different focal positions at different angles, affecting the accuracy of joint processing of multi-angle images. On the other hand, since the two-dimensional resolution of GNSS-based InSAR systems is low, directly using the raw elevation information for imaging will lead to shape distortion of the resolution cells, requiring fitting.

[0004] Traditional fitting algorithms include full-scene fitting algorithms and sub-region fitting algorithms. However, full-scene fitting algorithms have low fitting accuracy and cannot meet the elevation accuracy requirements of PS points (persistent scatterers, i.e., permanent scatterers, large target reflection points); sub-region fitting algorithms require segmentation of resolution cells, which disrupts the integrity of resolution cells and the continuity of elevation data. Summary of the Invention

[0005] In view of this, the present invention provides a navigation satellite bistatic InSAR high coherence DEM imaging method based on adaptive region segmentation. This method proposes the concept of image valleys, utilizes complete resolution cells (CRCs) and image valleys for adaptive sub-region segmentation, and then performs elevation fitting using second-order weighted least squares within each sub-region. Finally, highly correlated SAR images are obtained for interferometry to achieve high-precision deformation inversion results. This solves the problems of low elevation data fitting accuracy and continuity interruption in navigation satellite bistatic imaging algorithms, realizing a high-coherence imaging algorithm based on adaptive partitioning and elevation fitting suitable for GNSS-based InSAR systems. It improves elevation fitting accuracy, ensures the integrity of resolution cells and the continuity of elevation data, compensates for focus position deviations, eliminates resolution cell distortion, exhibits strong correlation between image pairs, and achieves high target deformation inversion accuracy.

[0006] The technical solution for implementing the present invention is as follows:

[0007] A navigation satellite bistatic InSAR high-coherence DEM imaging method based on adaptive region segmentation includes:

[0008] Obtain the coarse imaging result image and select the CRC points of the coarse imaging result image;

[0009] Based on the coordinates of each CRC point, determine the pixel coordinates of each pixel point corresponding to the image valley between different resolution units, and obtain the image valley boundary of each resolution unit;

[0010] Boundary gap compensation is performed on the valley boundaries of the image to obtain adaptive region segmentation results corresponding to the boundaries of each sub-region. The adaptive region segmentation results include multiple sub-regions.

[0011] Elevation fitting is performed in each of the sub-regions to generate SAR images.

[0012] Optionally, selecting the CRC points of the coarse imaging result image includes:

[0013] If we represent the pixel coordinates (x, y) of each pixel in the coarse imaging result image as a target vector, then the blur function between any target vector A and its neighboring target vector B is:

[0014]

[0015] In the above formula, Φ TA and Φ RAλ is the unit vector from the transmitter and receiver to the target vector A, respectively; λ is the wavelength; p is the beam distance after pulse compression; β is the bistatic angle of incidence; Θ is the unit vector along the bisector of angle β; c is the speed of light; m A This is the result of azimuth pulse compression of the beam; w E Ξ and Ξ represent the equivalent angular velocity and the unit vector in the equivalent direction of motion of the satellite relative to the target, respectively; j is an imaginary number representing the phase information of the image;

[0016] Pixels whose amplitude differs from the maximum value of each blur function max(χ(A,B)) by less than 3Db are selected to determine the theoretical 3Db resolution unit:

[0017] The local maximum value within the theoretical 3dB resolution cell range is used as the CRC point:

[0018]

[0019] In the above formula, I c This indicates the amplitude of the coarse imaging result image.

[0020] Optionally, the magnitude of the coordinates (x, y) of each pixel corresponding to the image valley is:

[0021]

[0022] In the above formula, f() represents the amplitude of the coarse imaging result image; (m,n) are the coordinates of each selected CRC point;

[0023] The directions of the coordinates (x, y) of each pixel corresponding to the image valley are:

[0024]

[0025] Optionally, the step of performing boundary gap compensation on the valley boundaries of the image to obtain adaptive region segmentation results corresponding to the boundaries of each sub-region includes:

[0026] Based on the preset image valley threshold value thre M Determine the coarse boundary E1 of each region to be divided:

[0027] E1={(x,y)|M(x,y)>thre M};

[0028] The coarse boundary E1 is tin-processed, and the maximum value of the image valley amplitude M(x,y) is selected along each image valley direction α(x,y) to obtain the tin-processed boundary E2 of each region to be divided:

[0029]

[0030] In the above formula, L represents a straight line passing through the CRC point and along the image valley direction α(x,y);

[0031] The small gaps in the tin-plated boundary E2 of the region to be divided are repaired using the minimum distance criterion, resulting in the repaired boundary E3 of each region to be divided.

[0032] Compensate for the long-distance gaps in the repair boundary E3 of the region to be divided, and determine the compensation boundary E4 of each region to be divided.

[0033] The regions to be divided that do not contain the CRC points are removed. Based on the area of ​​each region to be divided, the regions to be divided are merged to obtain an adaptive region segmentation result.

[0034] Optionally, the step of compensating for long-distance gaps in the repair boundary E3 of the regions to be divided, and determining the compensation boundary E4 of each region to be divided, includes:

[0035] For the breakpoint P of the long-distance gap d0 Record the breakpoint P d0 The previous boundary segment was L. d On line segment L d Find the match with P d0 N adjacent points are denoted as P. d1 ...P dN The reference direction γ of the compensation boundary and the breakpoint P d0 and its N neighboring points P d1 ...P dN The relationship between the valley directions in the image is as follows:

[0036]

[0037] Using θ as the angular tolerance, select point P, which is the point with the maximum value of the image valley amplitude M(x,y) within the range of (γ±θ). l As the next connection point:

[0038] Breakpoint P d0 With connection point P l Connect, and P l As the next breakpoint, repeat the above steps until the region is closed, and obtain the compensation boundary E4.

[0039] Optionally, the step of performing elevation fitting in each of the sub-regions to generate SAR images includes:

[0040] A second-order weighted least squares method is used to perform elevation fitting in each of the sub-regions;

[0041] Based on the fitted elevation of each sub-region, a SAR image is generated using a post-projection algorithm.

[0042] Optionally, the step of employing a second-order weighted least squares method to perform elevation fitting in each of the sub-regions includes:

[0043] The sub-region R Q The elevation values ​​of the K fitted points within the DEM are: DEM(x,y) = a·b;

[0044] In the above formula, a is the weighting matrix to be calculated, and b is the elevation matrix of each pixel, which are expressed as follows:

[0045]

[0046] Using the diagonal weighted matrix W, the second-order weighted least squares method is used to calculate the weighted matrix a to be found, and the estimated value of the weighted matrix a is obtained. for:

[0047]

[0048] In the above formula, the diagonal weighting matrix W is determined according to the amplitude of the coarse imaging result image:

[0049] W=diag(f(x1,y1),f(x2,y2),...,f(x K ,y K ));

[0050] The fitted elevations of each of the sub-regions are obtained as follows:

[0051]

[0052] Optionally, the SAR image I generated using the post-projection algorithm f for:

[0053] In the above formula, I f The generated SAR image; u is the slow time in the beam range direction; s is the range pulse compression result of the echo signal, representing the amplitude information of the image; R(x,y) is the bistatic range synchronized from the transmitter, expressed as: R(x,y)=|P S -[x,y,DEM fit [x,y)]|+|P E -[x,y,DEM fit [(x,y)]|-|P S -P E |;

[0054] In the above formula, PS P represents the satellite's three-dimensional position coordinates. E These are the position coordinates of the satellite's receiving antenna.

[0055] Beneficial effects:

[0056] (1) This invention ensures the continuity of elevation data and the integrity of resolution units, compensates for the deviation of the focusing position, improves the joint processing accuracy of multi-angle images, and eliminates the distortion of resolution units. In the interferometric measurement process, the coherence coefficient is higher than that of traditional algorithms, that is, the correlation between image pairs is high and the coherence phase accuracy is high, which enhances the elevation fitting accuracy and improves the inversion accuracy.

[0057] (2) When generating SAR images, the imaging method of the present invention significantly reduces the elevation fitting error through adaptive region segmentation, and can obtain higher precision fitting elevation. The interference phase error is small and the image pairs have a high coherence coefficient, which lays the foundation for obtaining higher deformation inversion accuracy. Attached Figure Description

[0058] Figure 1 This is a schematic diagram of the main process of the navigation satellite bistatic InSAR high coherence DEM imaging method based on adaptive region segmentation according to an embodiment of the present invention.

[0059] Figure 2 This is a schematic diagram of boundary gap compensation for long-distance gaps according to an embodiment of the present invention.

[0060] Figure 3 This is a schematic diagram illustrating the distribution of CRC points selected from the coarse imaging result image according to an embodiment of the present invention.

[0061] Figure 4 This is a schematic diagram of the valley boundary in an image according to an embodiment of the present invention.

[0062] Figure 5 This is a schematic diagram of the compensation boundary according to an embodiment of the present invention.

[0063] Figure 6 This is a schematic diagram of the adaptive region segmentation result according to an embodiment of the present invention.

[0064] Figure 7 This is a schematic diagram of a SAR image generated using a traditional fitting algorithm.

[0065] Figure 8 This is a schematic diagram of a SAR image generated according to an embodiment of the present invention. Detailed Implementation

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

[0067] Azimuth direction: The direction of the center of the SAR transmit beam.

[0068] Range direction: The direction perpendicular to the azimuth direction of the SAR transmitted beam.

[0069] This invention discloses a navigation satellite bistatic InSAR high-coherence DEM imaging method based on adaptive region segmentation, such as... Figure 1 As shown, it includes the following steps:

[0070] In this embodiment of the invention, CRC points are first selected from the coarse imaging result image; then, based on the calculated image valley results, the image boundary is determined, and the image region is adaptively segmented; finally, elevation fitting is performed in each sub-region to obtain a SAR image with high correlation, elevation continuity, and resolution cell integrity.

[0071] Step 1: Obtain the coarse imaging result image and select the CRC points of the coarse imaging result image.

[0072] In this embodiment of the invention, the coarse imaging result image is generated based on the original elevation data and GNSS echo data of the target scene. The original elevation data can be imported from existing data, and the GNSS echo data comes from the receiver.

[0073] In this embodiment of the invention, the CRC point is selected by the local maximum detection method. The pixel coordinates (x, y) of each pixel in the coarse imaging result image are represented as target vectors. Then, the fuzzy function between the target vector A of any pixel and the target vector B of the neighboring pixel is as shown in the following equation (1):

[0074]

[0075] In the above formula, Φ TA and Φ RA λ is the unit vector from the transmitter and receiver to the target vector A, respectively; λ is the wavelength; p is the beam distance after pulse compression; β is the bistatic angle of incidence; Θ is the unit vector along the bisector of angle β; c is the speed of light; m A This is the result of azimuth pulse compression of the beam; w E Ξ and Ξ represent the equivalent angular velocity and the unit vector in the equivalent direction of motion of the satellite relative to the target, respectively; j is an imaginary number representing the phase information of the image.

[0076] Furthermore, based on the maximum values ​​max(χ(A,B)) of equation (1), pixels with amplitude differences less than 3Db are selected, and the theoretical 3Db resolution unit is determined as shown in equation (2) below:

[0077]

[0078] Furthermore, the amplitude of the coarse imaging result image is expressed using I. c This indicates that the local maximum value within the theoretical 3dB resolution cell range is selected as the CRC point, as shown in equation (3) below:

[0079]

[0080] Step 2: Based on the coordinates of each CRC point, determine the pixel coordinates of each pixel point corresponding to the image valley between different resolution units, and obtain the image valley boundary of each resolution unit.

[0081] In this embodiment of the invention, image valleys are used to identify edge points between resolution units. The magnitude and direction of each pixel (x, y) corresponding to an image valley are defined as shown in equations (4) and (5), respectively:

[0082]

[0083]

[0084] In the above formula, f() represents the amplitude of the coarse imaging result image; (m,n) are the pixel coordinates of each selected CRC point.

[0085] Step 3: Perform boundary gap compensation on the valley boundaries of the image to obtain adaptive region segmentation results corresponding to the boundaries of each sub-region. The adaptive region segmentation results include multiple sub-regions.

[0086] Step 31: Determine the coarse boundaries of each region to be divided based on the preset threshold value of the image valley.

[0087] In this embodiment of the invention, the threshold value of the image valley is thre M You can configure it as needed, for example, thre M The threshold is 10 dB. Based on the image valley threshold value... M The coarse boundary E1 of each region to be divided is determined as shown in equation (6):

[0088] E1={(x,y)|M(x,y)>thre M} (6)

[0089] Step 32: The coarse boundary is tinned to obtain the tinned boundary of each region to be divided.

[0090] In this embodiment of the invention, the coarse boundary E1 of each region to be divided is tinned, that is, along each image valley direction α(x,y), the maximum value point of the image valley amplitude M(x,y) is selected to obtain the tinned boundary E2 of each region to be divided, as shown in the following formula (7):

[0091]

[0092] In the above formula, L represents a straight line passing through the CRC point and along the α(x,y) direction.

[0093] Step 33: Repair the small-distance gaps of the tin-plated boundary of the region to be divided using the minimum distance criterion to obtain the repaired boundary of each region to be divided.

[0094] In this embodiment of the invention, the minimum distance criterion is used to repair the small-distance gaps in the tin-bonded boundary E2 of each region to be divided. For example, the endpoints of the small-distance gaps are connected to form the repair boundary E3 of each region to be divided. The small-distance gaps can be selectively set as needed; for example, a distance of less than 10m can be defined as a small-distance boundary.

[0095] Step 34: Compensate for long-distance gaps in the repair boundaries of the regions to be divided, and determine the compensation boundaries of each region to be divided.

[0096] In this embodiment of the invention, for long-distance gaps other than small-distance gaps, such as... Figure 2 As shown, the breakpoint P for a long gap. d0 Record the breakpoint P d0 The previous boundary segment was L. d On line segment L d Find the match with P d0 N adjacent points are denoted as P. d1 ...P dN The reference direction γ of the compensation boundary and the breakpoint P d0 and its N neighboring points P d1 ...P dN The relationship between the valley directions in the image is shown in equation (8):

[0097]

[0098] Furthermore, using θ as the angular tolerance, the point P where the maximum value of the image valley amplitude M(x,y) within the range of (γ±θ) is selected. l As the next connection point, it is shown in equation (9) below:

[0099]

[0100] Breakpoint P d0 With connection point P l Connection, P l As the next breakpoint, repeat the above operation until the region is closed or M(P) is reached. l The value is very small, resulting in the compensation boundary E4 after long-distance gap compensation. Wherein, in M(P)l In cases where the value is very small, directly set the breakpoint P. d0 Simply connect it to the endpoint of the repair boundary.

[0101] Step 35: Remove the regions to be divided that do not contain CRC points, and merge the regions to be divided according to their area to obtain the adaptive region segmentation result.

[0102] In this embodiment of the invention, based on the compensation boundary E4, regions to be segmented that do not contain CRC points are removed. Then, based on the area of ​​each region to be segmented, regions with areas smaller than a preset area threshold are merged to obtain the final adaptive region segmentation result comprising the various sub-regions. For example, the merging method could be staggered sliding, merging regions smaller than the preset area threshold into adjacent regions larger than the preset area threshold.

[0103] Step 4: Perform elevation fitting in each of the sub-regions to generate SAR images.

[0104] In this embodiment of the invention, a second-order weighted least squares method is used to perform elevation fitting in each sub-region. Wherein, sub-region R... Q The elevation values ​​of the K fitted points within the range are shown in equation (10) below:

[0105] DEM(x,y)=a·b (10)

[0106] In the above formula, a is the weighting matrix to be calculated, and b is the elevation matrix of each pixel, as shown in the following formula (11):

[0107]

[0108] Using the diagonal weighted matrix W, the second-order weighted least squares method is used to calculate the weighted matrix to be found, and the estimated value of the weighted matrix a is obtained. As shown in equation (12):

[0109]

[0110] In the above formula, the diagonal weighting matrix W is determined according to the amplitude of the coarse imaging result image, as shown in the following formula (13):

[0111] W=diag(f(x1,y1),f(x2,y2),...,f(x K ,y K (13)

[0112] Accordingly, the fitted elevations for each sub-region are shown in equation (14) below:

[0113]

[0114] In this embodiment of the invention, a SAR image is generated using a post-projection algorithm based on the fitted elevation of each sub-region.

[0115] In this embodiment of the invention, the post-projection algorithm is a classic time-domain radar imaging algorithm. Based on the fitted elevation of each sub-region, the post-projection algorithm is used to generate a highly coherent SAR image, as shown in the following equation (15):

[0116]

[0117] In the above formula, I f The generated SAR image; u is the slow time in the beam range direction; s is the range pulse compression result of the echo signal, representing the amplitude information of the image; R(x,y) is the bistatic range synchronized from the transmitter, as shown in equation (16) below:

[0118] R(x,y)=|P S -[x,y,DEM fit [x,y)]|+|P E -[x,y,DEM fit [(x,y)]|-|P S -P E | (16)

[0119] In the above formula, P S P represents the satellite's three-dimensional position coordinates. E These are the position coordinates of the satellite's receiving antenna.

[0120] In this embodiment of the invention, HTM mine is used as the target scenario, the actual mine elevation is used as the original elevation data, and the transmitter is the BDMEO19 satellite. The parameters of the BDMEO19 satellite are shown in Table 1 below:

[0121] Table 1

[0122]

[0123] Using echo data from the BDMEO19 satellite on November 22, 29, December 6, and December 13, 2019, experiments were conducted using the adaptive region segmentation-based bistatic InSAR high-coherence DEM imaging method of this invention. The distribution of CRC points selected from the coarse imaging results is as follows: Figure 3 As shown; the image valley boundaries of each resolution unit are as follows Figure 4 As shown; the compensated boundary after boundary gap compensation is as follows. Figure 5 As shown, from Figure 5As can be seen, the boundary gaps of each region to be segmented are compensated and consistent with the CRC point distribution in the SAR image; the adaptive region segmentation result is as follows: Figure 6 As shown, from Figure 6 As can be seen, the distribution of each sub-region is consistent with the distribution of CRC points in the coarse imaging result image, and the edges of the sub-regions avoid the resolution cells, ensuring the integrity of each resolution cell; SAR images using traditional fitting algorithms, such as Figure 7 As shown, the SAR image generated by the fitted elevation determined using the imaging method of this invention is as follows: Figure 8 As shown, from Figure 7 and Figure 8 The comparison shows that the imaging method of the present invention fits the elevation data with high accuracy, and the generated SAR image is continuous and complete.

[0124] Furthermore, using the image from November 22, 2019 as the main image and images from the other three days as secondary images, SAR images were generated by performing elevation fitting using both traditional fitting algorithms and the imaging method of this invention. The average coherence coefficients between each image pair are shown in Table 2 below:

[0125] Table 2

[0126]

[0127] As can be seen from the table above, the coherence coefficients between the image pairs obtained by the imaging method of this invention are much higher than those of traditional fitting algorithms.

[0128] The specific embodiments described above only illustrate the design principles of the present invention. The shapes and names of the components in this description may differ and are not limited. Therefore, those skilled in the art can modify or make equivalent substitutions to the technical solutions described in the foregoing embodiments; and these modifications and substitutions do not depart from the inventive spirit and technical solutions of the present invention, and should all fall within the protection scope of the present invention.

Claims

1. A navigation satellite bistatic InSAR high-coherence DEM imaging method based on adaptive region segmentation, characterized in that, include: Acquire a coarse imaging result image, and select CRC points from the coarse imaging result image, including: If we represent the pixel coordinates (x, y) of each pixel in the coarse imaging result image as a target vector, then the blur function between any target vector A and its neighboring target vector B is: In the above formula, Φ TA and Φ RA λ is the unit vector from the transmitter and receiver to the target vector A, respectively; λ is the wavelength; p is the beam distance after pulse compression; β is the bistatic angle of incidence; Θ is the unit vector along the bisector of angle β; c is the speed of light; m A This is the result of azimuth pulse compression of the beam; w E Ξ and Ξ represent the equivalent angular velocity and the unit vector in the equivalent direction of motion of the satellite relative to the target, respectively; j is an imaginary number representing the phase information of the image; Pixels whose amplitude differs from the maximum value of each blur function max(χ(A,B)) by less than 3Db are selected to determine the theoretical 3Db resolution unit: The local maximum value within the theoretical 3dB resolution cell range is used as the CRC point: In the above formula, I c This indicates the amplitude of the coarse imaging result image; Based on the coordinates of each CRC point, the pixel coordinates of each pixel corresponding to the image valley between different resolution units are determined, thus obtaining the image valley boundary of each resolution unit; the magnitude of the coordinates (x, y) of each pixel corresponding to the image valley is: In the above formula, f() represents the amplitude of the coarse imaging result image; (m,n) are the coordinates of each selected CRC point; The directions of the coordinates (x, y) of each pixel corresponding to the image valley are: Boundary gap compensation is performed on the valley boundaries of the image to obtain adaptive region segmentation results corresponding to the boundaries of each sub-region. The adaptive region segmentation results include multiple sub-regions. Elevation fitting is performed in each of the sub-regions to generate SAR images.

2. The method as described in claim 1, characterized in that, The step of performing boundary gap compensation on the valley boundaries of the image to obtain adaptive region segmentation results corresponding to the boundaries of each sub-region includes: Based on the preset image valley threshold value thre M Determine the coarse boundary E1 of each region to be divided: E1={(x,y)|M(x,y)>thre M }; The coarse boundary E1 is tin-processed, and the maximum value of the image valley amplitude M(x,y) is selected along each image valley direction α(x,y) to obtain the tin-processed boundary E2 of each region to be divided: In the above formula, L represents a straight line passing through the CRC point and along the image valley direction α(x,y); The small gaps in the tin-plated boundary E2 of the region to be divided are repaired using the minimum distance criterion, resulting in the repaired boundary E3 of each region to be divided. Compensation is performed on the long-distance gaps in the repair boundary E3 of the regions to be divided, and the compensation boundary E of each region to be divided is determined. 4； The regions to be divided that do not contain the CRC points are removed. Based on the area of ​​each region to be divided, the regions to be divided are merged to obtain an adaptive region segmentation result.

3. The method as described in claim 2, characterized in that, The step of compensating for long-distance gaps in the repair boundary E3 of the region to be divided, and determining the compensation boundary E4 of each region to be divided, includes: For the breakpoint P of the long-distance gap d0 Record the breakpoint P d0 The previous boundary segment was L. d On line segment L d Find the match with P d0 N adjacent points are denoted as P. d1 ...P dN The reference direction γ of the compensation boundary and the breakpoint P d0 and its N neighboring points P d1 ...P dN The relationship between the valley directions in the image is as follows: Using θ as the angular tolerance, select point P, which is the point with the maximum value of the image valley amplitude M(x,y) within the range of (γ±θ). l As the next connection point: Breakpoint P d0 With connection point P l Connect, and P l As the next breakpoint, repeat the above steps until the region is closed, and obtain the compensation boundary E4.

4. The method as described in claim 3, characterized in that, The step of performing elevation fitting in each of the sub-regions to generate SAR images includes: A second-order weighted least squares method is used to perform elevation fitting in each of the sub-regions; Based on the fitted elevation of each sub-region, a SAR image is generated using a post-projection algorithm.

5. The method as described in claim 4, characterized in that, The method of employing second-order weighted least squares to perform elevation fitting in each of the sub-regions includes: The sub-region R Q The elevation values ​​of the K fitted points within the DEM are: DEM(x,y) = a·b; In the above formula, a is the weighting matrix to be calculated, and b is the elevation matrix of each pixel, which are expressed as follows: Using the diagonal weighted matrix W, the second-order weighted least squares method is used to calculate the weighted matrix a to be found, and the estimated value of the weighted matrix a is obtained. for: In the above formula, the diagonal weighting matrix W is determined according to the amplitude of the coarse imaging result image: W=diag(f(x1,y1),f(x2,y2),...,f(x K ,y K )); The fitted elevations of each of the sub-regions are obtained as follows:

6. The method as described in claim 5, characterized in that, The SAR image I generated using the post-projection algorithm f for: In the above formula, I f The generated SAR image; u is the slow time in the beam range direction; s represents the range pulse compression result of the echo signal, indicating the amplitude information of the image; R(x,y) is the bistatic range synchronized from the transmitter, expressed as: R(x,y)=|P S -[x,y,DEM fit [x,y)]|+|P E -[x,y,DEM fit [(x,y)]|-|P S -P E |; In the above formula, P S P represents the satellite's three-dimensional position coordinates. E These are the position coordinates of the satellite's receiving antenna.

Images