A multi-base satellite-borne SAR three-dimensional imaging method and device

Abstract

The present invention provides a multi-base satellite-borne SAR three-dimensional imaging method and device, belonging to the technical field of SAR image three-dimensional imaging. The method includes: image preprocessing; geocoding and registration; removing flat ground effect and terrain phase; atmospheric phase error compensation; high-precision tomographic SAR three-dimensional imaging. The device includes an external DEM data acquisition module, a registration module, a simulation module, a compensation module, and an imaging module. The present invention can obtain high-precision three-dimensional imaging results.

Description

Technical Field

[0001] The present invention belongs to the technical field of SAR image three-dimensional imaging, and in particular relates to a multi-base satellite-borne SAR three-dimensional imaging method and device. Background Art

[0002] In recent years, domestic and foreign scholars have devoted themselves to obtaining the three-dimensional electromagnetic backscattering information of the target in the field of SAR three-dimensional imaging, and solving the geometric distortion and information loss problems such as perspective contraction, overlap and shadow in the current two-dimensional SAR images. At present, the existing SAR three-dimensional imaging technologies mainly include curve SAR, array SAR and tomographic SAR technology. Interferometric SAR technology uses two or more SAR images for interference processing to obtain the elevation information of the ground target, but the elevation information obtained by this method does not have the resolution ability. Curve SAR technology obtains the echo data of the ground target in all directions and angles through the curved motion of the platform on the range-azimuth vertical plane to obtain the three-dimensional imaging results. However, this method has strict requirements on the flight trajectory of the SAR platform, the density of sampling points and the location information of the SAR platform. Array SAR three-dimensional imaging technology uses multiple antenna array elements to form multi-angle and multi-baseline coherent observations to achieve three-dimensional imaging. However, this method requires a sufficient number of array antennas and a good array design, has high requirements on hardware, is limited to airborne missions, and is not suitable for satellite missions. Tomographic SAR three-dimensional imaging technology uses multiple SAR images obtained from heavy-orbit or distributed multi-base satellites, and adopts special imaging algorithms to focus imaging in the altitude dimension to obtain high-resolution third-dimensional information of ground targets. It is an important development trend in the current three-dimensional imaging of the ground for multi-base / constellation SAR single-flight small data sets (a small number of images).

[0003] For high-precision tomographic 3D imaging, spectral analysis of angle of arrival estimation and compressed sensing methods of sparse estimation are introduced. The beam forming (BF) spectral estimation algorithm in the spectral analysis method was first used to solve the overlap problem. It is simple in principle and has high efficiency and robustness. It is widely used in tomographic SAR 3D imaging. However, in the case of non-uniform baseline and a small number of image distributions, the side lobes of the scattering power spectrum reconstructed by this method are high and it is difficult to break through the Rayleigh resolution. Therefore, many studies have introduced super-resolution estimation methods in the field of spatial spectrum estimation into tomography, mainly including the Capon algorithm, the multiple signal classification (MUSIC) algorithm, the singular value decomposition (SVD) algorithm, and the relaxation (RELAX) algorithm. These methods can break through the Rayleigh resolution and achieve super-resolution imaging. However, both the Capon and MUSIC methods require the estimation of the covariance matrix of the observed data. In the process of 3D reconstruction of buildings, it is often difficult to achieve a covariance matrix estimation with good robustness. Due to the uneven baseline and less baseline sampling in the observation area, these technologies usually bring some quality problems in the estimated resolution of the scatterer position, making it difficult to distinguish between the ground and the building. Therefore, these methods need to perform multi-view processing in the actual operation process, which reduces the resolution of SAR tomography in the range and azimuth directions and makes the height estimation inaccurate. The RELAX algorithm can achieve better SAR tomography and target detection effects through nonlinear least squares method estimation. In addition, the SVD method is applied to tomographic SAR three-dimensional imaging processing, and the noise suppression ability can be improved through the singular value truncation method (TSVD). With the development of multi-baseline tomographic SAR imaging technology, the resolution requirements for tomography are getting higher and higher. Related studies have shown that the number of scatterers in the same pixel in the urban building environment does not exceed 4, that is, the scatterers are sparse in the tomographic domain. This theory provides the possibility for CS to solve the tomographic SAR imaging problem. Algorithms based on CS theory include SL1MMER algorithm, regularized matching pursuit (ROMP), fast threshold iteration (FIST), etc. These methods use different sparse bases to image the ground and buildings, effectively overcoming the limitations of uneven baselines, improving imaging quality, and achieving high-precision three-dimensional imaging in cases where there is little track data and high noise. However, methods based on compressed sensing solutions rely heavily on hyperparameters, and the computational complexity is greatly increased, which is not conducive to real-time applications. Moreover, these methods have inaccurate hyperparameter estimation under the condition of small data sets, and the tomographic focusing performance is reduced.

[0004] Non-parametric spectrum analysis methods, such as beamforming, capon, and IAA, can obtain the position of building scatterers along the tomographic direction without any prior information. Parametric spectrum estimation methods, such as MUSIC and weighted subspace fitting, can obtain better building heights with prior information. However, in actual experiments, due to uneven baselines and less baseline sampling in the observation area, these technologies usually bring some quality problems in the estimated resolution of scatterer positions, making it difficult to distinguish between the ground and buildings. The compressed sensing method can achieve super-resolution imaging of building areas. Using different sparse bases to image the ground and buildings effectively overcomes the limitations of uneven baselines and improves imaging quality. However, the main disadvantage of these methods is that they require more data sets (more than 15 images) for reliable building height estimation, which is usually not suitable for the small amount of data obtained from a single pass of multi-base / distributed satellite. Due to the uneven distribution of satellites, small data set SAR images lead to reduced focusing performance. Small data sets have the characteristics of short data acquisition time and fast update speed of vertical structures of buildings.

[0005] The multi-baseline data used by tomographic SAR is generally acquired through heavy orbits or a small number of distributed satellites, and the number of images contained is small. This situation will cause under-sampling of the spatial spectrum in the tomographic direction. According to the Nyquist sampling theorem, spectrum aliasing will occur during tomographic imaging, resulting in blurred imaging results. Limited by the baseline span, the resolution of tomographic SAR imaging technology in the altitude direction is poor, usually in the order of tens of meters to tens of meters. The present invention aims to solve the challenges of small data set limitations, small baseline spans, and insufficient Rayleigh resolution, and is aimed at the urgent need for high-precision tomographic synthetic aperture radar (SAR) three-dimensional imaging. Summary of the invention

[0006] In order to solve the above technical problems, the present invention provides a multi-base satellite-borne SAR three-dimensional imaging method and device, which is a small data set tomographic SAR imaging method (DIARA) based on a double iterative adaptive residual method. The double iterative adaptive method does not require any prior knowledge, directly performs height focusing, and screens the height of the strongest scatterer. The outer iterative adaptive method obtains the scattering intensity value of the target scatterer by the least square method using the observation matrix and position information, and the residual function optimizes the signal model. The signal is divided into a main signal composed of several strong scatterers. In each iteration process, the signal of the last strong scatterer needs to be eliminated to obtain the remaining scatterer signal. In this method, the present invention believes that each pixel in the SAR image contains two main scatterers (ground and building). Through this processing, the backscattering characteristics of the main scatterer in the pixel are enhanced, and the noise and sidelobe information are effectively weakened, which makes this method can effectively extract the peaks of the ground and buildings, improve tomographic imaging, and help improve the accuracy of building height estimation.

[0007] In order to achieve the above object, the present invention adopts the following technical solution:

[0008] A multi-base spaceborne SAR three-dimensional imaging method comprises the following steps:

[0009] S1. Preprocess the multi-base satellite-borne SAR image set, obtain external DEM data according to the selected main image area, and generate main and secondary image intensity maps; DEM stands for digital elevation model;

[0010] S2. Generate a query list of coordinate mapping relationship between the main image and the corresponding external DEM data, then use the external DEM data to simulate the SAR intensity map, and realize the geocoding of the SAR main image and the acquisition of radar domain DEM data through the mutual correlation coefficient between the main image intensity map; perform registration between the main and secondary images based on the SAR main image and the radar domain DEM data;

[0011] S3, using satellite orbit data to calculate the parallel baseline and vertical baseline between multi-base spaceborne SAR images; converting the parallel baseline to flat ground phase and removing it in the primary and secondary image interferometry; using the vertical baseline and external DEM data to simulate the terrain phase, and removing it in the primary and secondary image interferometry, finally obtaining the SLC data set; SLC stands for single-look complex image;

[0012] S4, divide the target scene into multiple sub-areas, divide the sub-areas evenly according to a certain number, then calculate the amplitude deviation value of each sub-area, screen the PS points, analyze their phase stability, estimate the initial height, bring it into the phase stability analysis, and use the coherence threshold method to finely screen the PS points, and calculate the phase residual weight of the final finely screened PS points; perform linear phase removal, phase gradient estimation, and weighted phase estimation of the signal based on the obtained height, and finally perform atmospheric phase error compensation; PS represents permanent scatterer;

[0013] S5. Perform residual calculations, obtain the focused signal through double iterative adaptive spectrum estimation, retrieve the height of the strongest scatterer by establishing a search space, obtain the scattering intensity through the least squares method, and calculate the cost function using the scattering intensities obtained twice in a row. If the convergence condition is not met, repeat S5. If the convergence condition is met, a high-precision three-dimensional imaging result is obtained.

[0014] The present invention also provides a multi-base satellite-borne SAR three-dimensional imaging device, comprising the following modules:

[0015] The external DEM data acquisition module pre-processes the multi-base satellite-borne SAR image set, obtains external DEM data according to the selected primary image area, and generates primary and secondary image intensity maps; DEM stands for digital elevation model;

[0016] The registration module generates a query list of coordinate mapping relationships between the main image and the corresponding external DEM data, then uses the external DEM data to simulate the SAR intensity map, and realizes the geocoding of the SAR main image and the acquisition of radar domain DEM data through the mutual correlation coefficient between the main image intensity map; the registration between the main and secondary images is carried out based on the SAR main image and the radar domain DEM data;

[0017] The simulation module uses satellite orbit data to calculate the parallel baseline and vertical baseline between multi-base spaceborne SAR images; converts the parallel baseline into a flat ground phase and removes it in the primary and secondary image interferometer pair; simulates the terrain phase using the vertical baseline and external DEM data, and removes it in the primary and secondary image interferometer pair, and finally obtains the SLC data set; SLC stands for single-look complex image;

[0018] The compensation module divides the target scene into multiple sub-areas, divides the sub-areas evenly according to a certain number, and then calculates the amplitude deviation value of each sub-area, selects the PS points, analyzes their phase stability, estimates the initial height, and brings it into the phase stability analysis. The PS points are finely screened by the coherence threshold method, and the phase residual weights of the final finely screened PS points are calculated; the linear phase removal, phase gradient estimation, and weighted phase estimation of the signal are performed based on the obtained height, and finally the atmospheric phase error compensation is performed; PS represents permanent scatterer;

[0019] The imaging module performs residual calculations and obtains the focused signal through double iterative adaptive spectrum estimation. By establishing a search space, the height of the strongest scatterer is retrieved, and the scattering intensity is obtained through the least squares method. The cost function is calculated using the scattering intensities obtained twice in a row. If the convergence conditions are not met, the function of the imaging module is repeated. If the convergence conditions are met, a high-precision three-dimensional imaging result is obtained.

[0020] The present invention also provides an electronic device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the steps of the above-mentioned multi-base satellite-borne SAR three-dimensional imaging method when executing the program.

[0021] The present invention also provides a non-transitory computer-readable storage medium on which a computer program is stored. When the computer program is executed by a processor, the steps of the above-mentioned multi-base spaceborne SAR three-dimensional imaging method are implemented.

[0022] Beneficial effects:

[0023] 1. The first iteration of the present invention performs high-direction focusing, the second iteration improves the estimation accuracy, and the residual function optimizes the signal model;

[0024] 2. The present invention uses improved tomographic three-dimensional imaging processing to enhance the backscattering characteristics of the main scatterers in the pixels, and effectively weaken the noise and sidelobe information, which makes it possible to effectively extract the peaks of the ground and buildings, improve tomographic imaging, and help improve the accuracy of building height estimation. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 It is a flow chart of the multi-base satellite-borne SAR three-dimensional imaging method of the present invention;

[0026] Figure 2 It is a schematic diagram of the geometric model of synthetic aperture radar interferometry;

[0027] Figure 3 is the phase error distribution map based on sub-region segmentation;

[0028] Figure 4 This is a schematic diagram of the three-dimensional imaging results of a single flight of multi-base satellite-borne SAR. DETAILED DESCRIPTION

[0029] The specific implementation modes of the present invention are described below so that those skilled in the art can understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific implementation modes. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations utilizing the concept of the present invention are protected.

[0030] like Figure 1 As shown, the multi-base satellite-borne SAR three-dimensional imaging method of the present invention comprises the following steps:

[0031] S1. Image preprocessing: Preprocess the images in the multi-base satellite-borne SAR single-look complex image (SLC) data set. Select the main image according to the temporal or spatial characteristics of the image for accurate registration of the main and secondary images and removal of terrain phase. Select appropriate digital elevation model (DEM) data according to the spatial range of the main image. In order to improve the solution efficiency and save computing resources, the main and secondary images and external DEM data are cropped or spliced ​​according to the size of the area of ​​interest, and the intensity map of the cropped main and secondary images is generated to prepare for subsequent registration.

[0032] S2. Geocoding and registration: geocoding the main image and accurately registering the main and secondary images. In geocoding the main image, external DEM data is used to register with the main image, and DEM data is used to simulate the SAR intensity map under the geometric framework of the main image, so as to generate radar domain DEM data, in preparation for the subsequent precise registration of the main and secondary images. The precise registration of the main and secondary images uses radar domain DEM data to reduce the offset caused by terrain factors and improve the accuracy of registration.

[0033] S3. Remove the flat ground effect and terrain phase. Remove the phase components caused by the reference ellipsoid and terrain undulations in the interference phase composed of the primary and secondary images. Calculate the spatial baseline (including vertical baseline and parallel baseline) of the primary and secondary image interference pair using satellite orbit data. Calculate the flat ground phase component through the parallel baseline. Use the vertical baseline and DEM data to simulate the terrain phase component. Finally, remove these two phase components from the interference phase of each interference pair.

[0034] S4, atmospheric phase error compensation, in order to improve the efficiency of atmospheric phase error estimation, the target scene is divided into multiple areas, and the sub-areas are evenly divided according to a certain number. Then the amplitude deviation value is calculated for each sub-area, and the permanent scatterer (PS) points are screened, and the phase stability is analyzed. The initial height is estimated and brought into the phase stability analysis. The coherence threshold method is used to finely screen the PS points, and the phase residual weights of the final screened PS points are calculated. The linear phase removal, phase gradient estimation, and weighted phase estimation of the signal are performed based on the obtained height, and finally the atmospheric phase error compensation is performed.

[0035] S5. High-precision tomographic SAR three-dimensional imaging. For high-precision tomographic SAR three-dimensional imaging, a scatterer model is established. First, residual calculation is performed, and the focused signal is obtained through double iterative adaptive spectrum estimation. By establishing a search space, the height of the strongest scatterer is retrieved, and the scattering intensity is obtained through the least squares method. The two adjacent signals are calculated through the cost function. If the convergence condition is not met, the previous steps are repeated. If the convergence condition is met, a high-precision three-dimensional imaging result is obtained.

[0036] Specifically, the S2 further includes:

[0037] S21, generating a query list of coordinate mapping relationships between the main image and the external DEM data using the main image track data and the external DEM data;

[0038] S22, preliminarily simulate the intensity map under the geometric framework of the SAR main image by querying the table and using the radar scattering model using external DEM data;

[0039] S23, using the spaceborne SAR orbit data, an initial offset estimation is performed between the intensity map of the SAR main image geometric frame simulated by the external DEM data and the main image intensity map;

[0040] S24, using the mutual correlation coefficient method to accurately align the intensity map under the SAR main image geometric framework with the SAR main image intensity map, and resample to generate radar domain DEM data consistent with the resolution and range of the SAR main image;

[0041] S25, using radar domain DEM data and primary and secondary image parameter files, generating an initial query list for primary and secondary image registration taking into account terrain offset;

[0042] S26. Refine the initial query list of S25 using the image domain pixel correlation coefficient method, and resample the secondary image to the main image geometric space through the refined query list to achieve precise registration between the main and secondary images.

[0043] Specifically, the S24 includes:

[0044] The calculation model for primary and secondary image registration using the mutual correlation coefficient is as follows:

[0045] (1)

[0046] in, , Represents the complex value at the row and column number (n, m) of the image coordinates in the primary and secondary image data blocks. M and N represent the size of the data block for calculating the coherence. That is, taking any pixel as the center and calculating according to the window size of N×M, the correlation coefficient value can be obtained. The correlation coefficient of each point is used as an indicator to evaluate the quality of the pixel interference phase.

[0047] Specifically, the S3 further includes:

[0048] S31, calculating the parallel baseline and the vertical baseline according to the orbit data of the primary and secondary images;

[0049] S32, converting the parallel baseline into a flat ground phase, and removing the flat ground phase from the primary and secondary image interference pair phases;

[0050] S33. Use the vertical baseline and external DEM data to simulate the terrain phase, and remove the terrain phase from the primary and secondary image interferometer phase to obtain the SLC data set.

[0051] Specifically, the S3 includes:

[0052] like Figure 2As shown, O is the coordinate of the geocentric origin in the geocentric earth-fixed coordinate system. The geometric model for imaging the same ground target point P, when the satellite is at position S1, the SAR sensor transmits electromagnetic waves to the ground target point P at a certain radar side angle of view or radar wave incident angle θ and receives its backscattered signal. At this time, the one-way distance of the electromagnetic wave along position S1 to the ground target point P is R1. When the satellite arrives at position S2 near position S1 again, the SAR sensor images the ground target point P again. At this time, the one-way distance of the electromagnetic wave along position S2 to the ground target point P is R2. The distance when the SAR sensor images the same ground area at two different times is called the spatial baseline, denoted by B; the angle between the baseline and the horizontal direction is called the baseline inclination, denoted by α; the projection of the spatial baseline on the satellite's line of sight to the distance R1 is called the parallel baseline B ∥ ; The vertical distance from the satellite position S2 point to the line of sight R1 is called the vertical baseline B ⊥ The height of the satellite platform is represented by H, and the geodetic height of the ground target point P is h. The geodetic height is the distance from the ground point along the normal line of the reference ellipsoid (WGS-84 reference ellipsoid) to the reference ellipsoid surface.

[0053] Furthermore, let the slant range difference between the two satellites be ,Depend on Figure 2 The spatial geometric relationship is known, and the present invention can use the triangular cosine law to write the slant range difference of the satellite during two imaging into a mathematical expression related to the spatial baseline as follows:

[0054] (2)

[0055] because and , so formula (2) can be written as:

[0056] (3)

[0057] In order to calculate the reference ellipsoid phase, first project the natural ground target point P onto the reference ellipsoid according to the slant distance R1 of the main image observation, and obtain the reference target point P0. At this time, the distance from the satellite position S1 corresponding to the main image to the reference target point P0 is r1, and r1=R1. The side angle of the SAR sensor to the reference target point P0 is θ0, and the reference target point P0 is connected with the satellite position S2 corresponding to the secondary image, and the distance from S2 to the reference target point P0 is r2. Therefore, the reference ellipsoid phase can be obtained. The mathematical expression of (i.e., the flat-earth phase component) is:

[0058] (4)

[0059] Where λ is the radar wavelength. In addition, the geometric relationship between the spatial baseline and the slant range can be used, and combined with formula (2), it can be obtained that:

[0060] (5)

[0061] Therefore, the flat-earth phase component It can also be written as:

[0062] (6)

[0063] Terrain phase component for:

[0064] (7)

[0065] Among them, δθ represents the change in radar side view angle caused by the change from reference target point P0 to ground target point P. Considering that the earth curvature in the local observation area is very small, the reference ellipsoid can be regarded as a plane. According to the spatial geometric relationship, we can get:

[0066] (8)

[0067] Substituting equation (8) into equation (7) we can continue to convert the terrain phase component Written as:

[0068] (9)

[0069] Specifically, the S4 further includes:

[0070] S41, dividing the SLC data set obtained in S3 into several sub-regions, and then calculating the amplitude deviation values ​​of the sub-regions, and obtaining PS points by screening through threshold setting;

[0071] S42, calculating the initial height of the obtained PS point through phase stability analysis;

[0072] S43, after the initial height is brought into the phase stability model, the PS points are further screened by the coherence threshold method;

[0073] S44, using the residual signal phase to estimate the weight of the PS point after screening to obtain a weight;

[0074] S45. The initial height of the PS point after screening is obtained through S42. At this time, the signal model can refer to formula (11):

[0075] S46, perform linear phase removal to obtain a new signal model, referring to formula (13);

[0076] S46. The new signal model is used for differential calculation to obtain the phase gradient value of each PS point:

[0077] S47, performing weighted phase estimation on the obtained phase gradient value of the PS point and the weight obtained in S44, and then removing the obtained phase from the SAR image.

[0078] Specifically, the S4 includes:

[0079] like Figure 3 As shown, the sub-region is evenly divided into several sub-regions, and the amplitude deviation value is estimated in the sub-regions. and They represent the Pth block in distance and the Qth block in azimuth respectively, and N represents the number of data sets. Represents the phase error value to be estimated for the Pth block in the distance and the Qth block in the azimuth of the Nth image.

[0080] (10)

[0081] in, represents the amplitude deviation value, n represents the image index value, and there are N images in total. (11)

[0082] in, represents the scattering intensity of the scatterer, Represents the pixel complex value of N SLC images, v represents the height, is the spatial frequency, which depends on the vertical baseline of the space and slant distance R, the phase error value is , j represents an imaginary unit, and exp() represents an exponential function. The selected PS point is analyzed using the phase stability model to obtain the height value. The phase stability formula is as follows:

[0083] (12)

[0084] Among them, the phase , arg{} represents the phase value. is the phase stability value, represents the unconstrained optimization solution, Represents the estimated height value.

[0085] Perform linear phase removal on the obtained observation vector;

[0086] (13)

[0087] Phase gradient estimation:

[0088] (14)

[0089] in, represents the phase gradient, represent , represents the conjugation operation, Represents the imaginary part operation.

[0090] The phase gradient is modeled and the gradient is calculated by the weighted least squares method:

[0091] (15)

[0092] in, represents the phase gradient value to be solved, K represents the number of PS points, Represents the phase gradient value of the Kth PS point of the nth SLC. represent The vector of represents transpose, , represents the weight diagonal matrix, all values ​​outside the diagonal are 0, and the weight can be recorded as:

[0093] (16)

[0094] in, represents the weight value of the diagonal elements, j represents the imaginary unit, Represents the phase obtained by the estimated height value, denoted as .

[0095] The phase gradient is integrated to estimate the phase error:

[0096] (17)

[0097] in, represents the estimated value of the phase error value, and p represents the index value when summing the phase gradient.

[0098] The obtained phase error is applied to N SLCs to complete the atmospheric phase compensation.

[0099] Specifically, the S5 further includes:

[0100] S51, estimation of the number of scatterers, usually two scatterers are considered, and if there are less than two scatterers, it will not affect the three-dimensional imaging results;

[0101] S52, establishing a scatterer tomography signal model, first performing a residual calculation to obtain a signal after residual calculation; because there is no signal of a strong scatterer at this time, the signal value subtracted is initialized to 0;

[0102] S53, performing double iterative adaptive spectrum estimation on the signal after residual calculation, and obtaining the height value of the current strongest scatterer by setting the search space;

[0103] S54, obtain the scattering intensity by using the least square method for the height value of the current strongest scatterer, calculate the two adjacent signals by the cost function, if the convergence condition is not met, repeat S53-S54, if the convergence condition is met, obtain a high-precision three-dimensional imaging result.

[0104] Specifically, the S5 includes:

[0105] After the number of scatterers Q is estimated, the signal is modeled first:

[0106] (18)

[0107] Then the residual is calculated, and the scatterer height is obtained by iterative adaptive method. Assume that P is the power of all discrete scattering intensities. A diagonal matrix whose diagonal is the power of the discrete sampled elevation to position, is the tomographic power spectrum, the superscript T represents the transpose of the matrix, ~ is the discrete value of height;

[0108] (19)

[0109] (20)

[0110] Where A represents the steering vector matrix, Represents the pixel value of N images, diag() represents a diagonal matrix, ~ Indicates the spatial frequency of the 1st to Nth images.

[0111] Before iteration, the tomographic power spectrum is transformed using the Fourier beamforming result. Initialize, that is, the data covariance matrix , I is the identity matrix.

[0112] During the iteration process, the data covariance matrix R is calculated as follows:

[0113] (twenty one)

[0114] in, is the conjugate transpose operator, the intermediate matrix .

[0115] The estimated value of the height of the phase center of the scattering target can be obtained by detecting the tomographic power spectrum The peak value is obtained, that is:

[0116] (twenty two)

[0117] Then, the least squares estimate can be used to obtain Estimated value of :

[0118] (twenty three)

[0119] in, .

[0120] The use of iterative methods and residual optimization methods is to improve the estimation accuracy of each scattering target. The convergence condition in the iterative method can be obtained by a cost function Determine, its expression is:

[0121] (twenty four)

[0122] in, is the L2 norm, and the convergence condition can be defined as , i represents the number of iterations. Represents the change in the cost function between two adjacent iterations. represents the scattering intensity of the qth scatterer, represents the height of the qth scatterer, and Q represents the total number of scatterers.

[0123] Figure 4 Schematic diagram of the imaging results.

[0124] In summary, the present invention has good applicability to single-pass multi-base satellite-borne SAR networking observation or few-pass satellite-borne SAR to achieve three-dimensional imaging of ground objects, and can achieve high-efficiency and high-precision three-dimensional imaging of urban areas.

[0125] The present invention also provides a multi-base satellite-borne SAR three-dimensional imaging device, comprising the following modules:

[0126] The external DEM data acquisition module pre-processes the multi-base satellite-borne SAR image set, obtains external DEM data according to the selected primary image area, and generates primary and secondary image intensity maps; DEM stands for digital elevation model;

[0127] The registration module generates a query list of the coordinate mapping relationship between the main image and the corresponding external DEM data, then uses the DEM to simulate the SAR intensity map, and realizes the geocoding of the SAR main image and the acquisition of the radar domain DEM through the mutual correlation coefficient between the SAR main image and the radar domain DEM; the registration between the main and secondary images is carried out based on the SAR main image and the radar domain DEM;

[0128] The simulation module uses satellite orbit data to calculate the parallel baseline and vertical baseline between multi-base spaceborne SAR images; converts the parallel baseline into a flat ground phase and removes it in the primary and secondary image interferometer pair; simulates the terrain phase using the vertical baseline and external DEM data, and removes it in the primary and secondary image interferometer pair, and finally obtains the SLC data set; SLC stands for single-look complex image;

[0129] The compensation module divides the target scene into multiple sub-areas, divides the sub-areas evenly according to a certain number, and then calculates the amplitude deviation value of each sub-area, selects the PS points, analyzes their phase stability, estimates the initial height, and brings it into the phase stability analysis. The PS points are finely screened by the coherence threshold method, and the phase residual weights of the final finely screened PS points are calculated; the linear phase removal, phase gradient estimation, and weighted phase estimation of the signal are performed based on the obtained height, and finally the atmospheric phase error compensation is performed; PS represents permanent scatterer;

[0130] The imaging module performs residual calculations and obtains the focused signal through double iterative adaptive spectrum estimation. By establishing a search space, the height of the strongest scatterer is retrieved, and the scattering intensity is obtained through the least squares method. The cost function is calculated using the scattering intensities obtained twice in a row. If the convergence conditions are not met, the function of the imaging module is repeated. If the convergence conditions are met, a high-precision three-dimensional imaging result is obtained.

[0131] The present invention also provides an electronic device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the steps of the above-mentioned multi-base satellite-borne SAR three-dimensional imaging method when executing the program.

[0132] The present invention also provides a non-transitory computer-readable storage medium on which a computer program is stored. When the computer program is executed by a processor, the steps of the above-mentioned multi-base spaceborne SAR three-dimensional imaging method are implemented.

[0133] It should be understood by those skilled in the art that the embodiments of the present invention may be provided as methods, systems, or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program codes. The solutions in the embodiments of the present invention may be implemented in various computer languages, for example, object-oriented programming language Java and interpreted scripting language JavaScript, etc.

[0134] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0135] These computer program instructions may also be stored in a computer readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture including an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0136] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0137] Although the preferred embodiments of the present invention have been described, those skilled in the art may make other changes and modifications to these embodiments once they have learned the basic creative concept. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications that fall within the scope of the present invention.

[0138] Obviously, those skilled in the art can make various changes and modifications to the present invention without departing from the spirit and scope of the present invention. Thus, if these modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalents, the present invention is also intended to include these modifications and variations.

Claims

1. A multi-base spaceborne SAR three-dimensional imaging method, characterized in that: The steps include: S1. Preprocess the multi-base satellite-borne SAR image set, obtain external DEM data according to the selected main image area, and generate main and secondary image intensity maps; DEM stands for digital elevation model; S2. Generate a query list of coordinate mapping relationship between the main image and the corresponding external DEM data, then use the external DEM data to simulate the SAR intensity map, and realize the geocoding of the SAR main image and the acquisition of radar domain DEM data through the mutual correlation coefficient between the main image intensity map; perform registration between the main and secondary images based on the SAR main image and the radar domain DEM data; S3, using satellite orbit data to calculate the parallel baseline and vertical baseline between multi-base spaceborne SAR images; converting the parallel baseline to flat ground phase and removing it in the primary and secondary image interferometry; using the vertical baseline and external DEM data to simulate the terrain phase, and removing it in the primary and secondary image interferometry, finally obtaining the SLC data set; SLC stands for single-look complex image; S4, divide the target scene into multiple sub-areas, divide the sub-areas evenly according to a certain number, then calculate the amplitude deviation value of each sub-area, screen the PS points, analyze their phase stability, estimate the initial height, bring it into the phase stability analysis, and use the coherence threshold method to finely screen the PS points, and calculate the phase residual weight of the final finely screened PS points; perform linear phase removal, phase gradient estimation, and weighted phase estimation of the signal based on the obtained height, and finally perform atmospheric phase error compensation; PS represents permanent scatterer; S5. Perform residual calculations, obtain the focused signal through double iterative adaptive spectrum estimation, retrieve the height of the strongest scatterer by establishing a search space, obtain the scattering intensity through the least squares method, and calculate the cost function using the scattering intensities obtained twice in a row. If the convergence condition is not met, repeat S5. If the convergence condition is met, a high-precision three-dimensional imaging result is obtained.

2. The multi-base spaceborne SAR three-dimensional imaging method according to claim 1, characterized in that: The S1 includes: S11, selecting a multi-base satellite-borne SAR single-view complex image set according to the target area, and selecting an image as a SAR main image according to the SAR image acquisition time in the multi-base satellite-borne SAR single-view complex image set; S12, selecting external DEM data according to the range of the SAR main image; S13, cropping or splicing the SAR primary and secondary images and the corresponding DEM range according to the area of ​​interest; S14. Obtain a SAR primary image intensity map and a SAR secondary image intensity map by calculating the SAR single-view complex image of the cropped or spliced ​​area.

3. The multi-base spaceborne SAR three-dimensional imaging method according to claim 2, characterized in that: The S2 includes: S21, generating a query list of coordinate mapping relationship between the main image and the DEM using the main image track data and the external DEM data; S22, by querying the table and adopting the radar scattering model, the intensity map under the geometric framework of the SAR main image is preliminarily simulated using external DEM data; S23, using the spaceborne SAR orbit data, an initial offset estimation is performed between the intensity map of the SAR main image in the geometric frame simulated by the external DEM data and the intensity map of the SAR main image; S24, using the mutual correlation coefficient method to accurately align the intensity map of the SAR main image in the geometric frame simulated by the external DEM data with the intensity map of the SAR main image, and resample to generate radar domain DEM data consistent with the resolution and range of the SAR main image; S25, using radar domain DEM data and primary and secondary image parameter files, generating an initial query list for primary and secondary image registration taking into account terrain offset; S26, using the image domain pixel mutual correlation coefficient method to refine the primary and secondary image registration initial query list of S25, and resample the secondary image to the primary image geometric space through the refined primary and secondary image registration initial query list to achieve precise registration between the primary and secondary images.

4. The multi-base spaceborne SAR three-dimensional imaging method according to claim 1, characterized in that: The S3 includes: S31, calculating the parallel baseline and the vertical baseline according to the orbit data of the primary and secondary images; S32, converting the parallel baseline into a flat ground phase, and removing it from the main and secondary image interference pair phases; S33. Use the vertical baseline and external DEM data to simulate the terrain phase, and remove it from the primary and secondary image interferometric phase to obtain the SLC data set.

5. The multi-base spaceborne SAR three-dimensional imaging method according to claim 4, characterized in that: The S4 includes: S41, dividing the SLC data set obtained in S3 into several sub-regions, and then calculating the amplitude deviation values ​​of the sub-regions, and obtaining PS points by screening through threshold setting; S42, calculating the initial height of the PS points obtained by screening through phase stability analysis; S43, after the initial height is brought into the phase stability model, the PS points are further screened by the coherence threshold method; S44, performing weight estimation on the PS points after screening using the residual signal phase to obtain a weight value; S45, using S42 to obtain the initial height of the PS point after screening, performing linear phase removal to obtain a new signal model; S46, performing differential calculation on the new signal model to obtain the phase gradient value of each PS point; S47, performing weighted phase estimation on the obtained phase gradient value of the PS point and the weight value obtained in S44, and removing the obtained phase from the SAR image.

6. The multi-base spaceborne SAR three-dimensional imaging method according to claim 1, characterized in that: The S5 includes: S51. The number of scatterers is set to two.

7. The multi-base spaceborne SAR three-dimensional imaging method according to claim 6, characterized in that: The S5 further includes: S52, establishing a scatterer tomography signal model, and obtaining a signal after residual calculation through a one-step residual calculation.

8. The multi-base spaceborne SAR three-dimensional imaging method according to claim 7, characterized in that: The S5 further includes: S53, performing double iterative adaptive spectrum estimation on the signal after residual calculation, and obtaining the height value of the current strongest scatterer by setting the search space.

9. The multi-base spaceborne SAR three-dimensional imaging method according to claim 8, characterized in that: The S5 further includes: S54, using the height value of the current strongest scatterer to obtain the scattering intensity through the least squares method, and performing cost function calculation, if the convergence condition is not met, repeat S53-S54, if the convergence condition is met, a high-precision three-dimensional imaging result is obtained.

10. A multi-base spaceborne SAR three-dimensional imaging device, characterized in that: Includes the following modules: The external DEM data acquisition module pre-processes the multi-base satellite-borne SAR image set, obtains external DEM data according to the selected primary image area, and generates primary and secondary image intensity maps; DEM stands for digital elevation model; The registration module generates a query list of coordinate mapping relationships between the main image and the corresponding external DEM data, then uses the external DEM data to simulate the SAR intensity map, and realizes the geocoding of the SAR main image and the acquisition of radar domain DEM data through the mutual correlation coefficient between the main image intensity map; the registration between the main and secondary images is carried out based on the SAR main image and the radar domain DEM data; The simulation module uses satellite orbit data to calculate the parallel baseline and vertical baseline between multi-base spaceborne SAR images; converts the parallel baseline into a flat ground phase and removes it in the primary and secondary image interferometer pair; simulates the terrain phase using the vertical baseline and external DEM data, and removes it in the primary and secondary image interferometer pair, and finally obtains the SLC data set; SLC stands for single-look complex image; The compensation module divides the target scene into multiple sub-areas, divides the sub-areas evenly according to a certain number, and then calculates the amplitude deviation value of each sub-area, selects the PS points, analyzes their phase stability, estimates the initial height, and brings it into the phase stability analysis. The PS points are finely screened by the coherence threshold method, and the phase residual weights of the final finely screened PS points are calculated; the linear phase removal, phase gradient estimation, and weighted phase estimation of the signal are performed based on the obtained height, and finally the atmospheric phase error compensation is performed; PS represents permanent scatterer; The imaging module performs residual calculations and obtains the focused signal through double iterative adaptive spectrum estimation. By establishing a search space, the height of the strongest scatterer is retrieved, and the scattering intensity is obtained through the least squares method. The cost function is calculated using the scattering intensities obtained twice in a row. If the convergence conditions are not met, the function of the imaging module is repeated. If the convergence conditions are met, a high-precision three-dimensional imaging result is obtained.

11. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the program, the steps of the multi-base spaceborne SAR three-dimensional imaging method according to any one of claims 1-9 are implemented.

12. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the multi-base spaceborne SAR three-dimensional imaging method described in any one of claims 1-9 are implemented.

Images