Multi-angle interference SAR (Synthetic Aperture Radar) three-dimensional reconstruction method based on overlay shadow detection

Abstract

The invention discloses a multi-angle interference SAR (Synthetic Aperture Radar) three-dimensional reconstruction method based on overlay shadow detection, and belongs to the field of interference SAR three-dimensional reconstruction. The method comprises the following steps: firstly, carrying out interference processing and correction on interference SAR data irradiated at each angle to obtain corrected DEM data; meanwhile, calculating a coherence coefficient of the main and auxiliary images at each angle to obtain a coherence coefficient graph; and based on the local frequency and the characteristic value, jointly detecting the overlay shadow area at each angle to obtain an overlay shadow detection value. Then, combining the coherence coefficient graph and the overlay shadow detection value to set the fusion weight of each pixel point in the interference SAR data of each angle, and obtaining a weight graph; and projecting the weight map and the corrected DEM value to a unified coordinate system. And finally, fusing the projected multi-angle interference SAR data based on the fusion weight of the coherence coefficient and the overlay shadow detection mean value, and carrying out interpolation in a null value region to obtain a fused DEM, thereby completing three-dimensional reconstruction of the region to be detected. The elevation precision is improved, and the application range is wide.

Description

Technical Field

[0001] The present invention belongs to the field of interferometric SAR three-dimensional reconstruction, and relates to a multi-angle interferometric SAR three-dimensional reconstruction method based on layover and shadow detection. Background Art

[0002] Interferometric Synthetic Aperture Radar (InSAR) technology can obtain terrain elevation or deformation information through interferometric phase information. With the ability of all-weather and all-day monitoring over a large range, InSAR technology has been widely used in terrain mapping, geological disaster monitoring, urban planning and other fields. However, the side-looking imaging characteristics of the SAR system lead to large areas of layover and shadow phenomena in complex terrain areas (such as steep mountainous areas). Among them, layover means that due to the large terrain relief slope, multiple ground object targets overlap in the SAR image, resulting in phase information confusion and unable to accurately reflect the true terrain; shadow means that due to ground object occlusion, the radar signal cannot reach some areas, which is shown as an information-free area in the SAR image. No effective phase information can be obtained in these two areas, resulting in serious missing or errors in the generated DEM (Digital Elevation Model) data in these areas, which limits the application of InSAR technology in complex terrain areas.

[0003] In order to improve the elevation accuracy and reliability, multi-source InSAR data such as multi-angle data can be introduced and weighted fusion can be carried out at the elevation level. The existing weighted fusion methods are mainly divided into two categories: one is the fusion method based on the coherence coefficient, which uses the coherence coefficient between interferometric SAR data as the weight for fusion [1][3] ; the other is the fusion method based on the elevation error map, which combines the coherence coefficient and the fuzzy height to calculate the elevation error map and performs weighted fusion on multi-source DEM data [4][5] . However, both of these two methods cannot accurately distinguish the layover and shadow areas, resulting in large errors in the fusion results of these areas, limited fusion effects, and unable to effectively solve the problems of missing and error of DEM data in complex terrain areas.

[0004] [1] E.Sansosti et al., “Digital elevation model generation using ascending and descending ERS-1 / ERS-2 tandem data,” International Journal of Remote Sensing, vol.20, no.8, pp.1527–1547, 1999.

[0005] [2]M.Crosetto, “Calibration and validation of SAR interferometry for DEM generation,” ISPRS Journal of Photogrammetry and Remote Sensing, vol.57, no.3, pp.213–227, Dec.2002.

[0006] [3]M.Eineder, “Interferometric DEM reconstruction of Alpine areas–experiences with SRTM data and improved strategies for future missions,” Proceedings of EARSEL 3D workhop, Porto, 2005.

[0007] [4]H.Jiang et al., “Fusion of high resolution DEMs derived from COSMOSkyMed and TerraSAR-X InSAR datasets,” Journal of Geodesy, vol.88, no.6, pp.587–599, Jun.2014.

[0008] [5]Y.Dong et al., “Cascaded multi-baseline interferometry with bistatic TerraSAR-X TanDEM-X observations for DEM generation,” ISPRS Journal of Photogrammetry and Remote Sensing, vol.171, pp.224-237, Jan.2021. Summary of the Invention

[0009] The present invention mainly aims at the problem of DEM data loss and error caused by the layover shadow phenomenon in complex areas, and proposes a multi-angle interferometric SAR three-dimensional reconstruction method based on layover shadow detection. This method is based on multi-angle interferometric SAR data in different directions, introduces layover shadow detection, and fuses the multi-angle interferometric SAR data by combining the mean value of layover shadow detection and the coherence coefficient at the elevation level, and finally obtains high-precision DEM data.

[0010] The multi-angle interferometric SAR three-dimensional reconstruction method based on layover shadow detection includes the following steps:

[0011] Step 1, the synthetic aperture radar collects multi-angle imaging information of the area to be detected, performs interferometric processing on the interferometric SAR data irradiated at each angle, and generates DEM data corresponding to each angle;

[0012] The specific process of performing interferometric processing on the interferometric SAR data at each angle is as follows:

[0013] 1) Register the interferometric SAR single-look complex image pairs irradiated at each angle respectively, and then conjugate multiply to obtain the corresponding interferometric phase diagram;

[0014] 2) Perform flat-earth removal on the interferometric phase diagram based on interferometric spectral shift, and then use a non-local filtering method based on non-linear fringe compensation to filter out the noise phase to obtain the filtered interferometric phase diagram;

[0015] 3) Perform phase unwrapping on the filtered interferometric phase diagram using the statistical cost flow method.

[0016] 4) Elevation inversion: Based on the unwrapped phase, use the Newton iteration method, and combine the interferometric phase equation, slant range equation, and Doppler equation to jointly solve the three-dimensional coordinates of the target point to obtain DEM data.

[0017] The interferometric phase equation is:

[0018]

[0019] The slant range equation of the SAR main image is:

[0020] (|S 1S -P| + |S 1R -P|) = 2r1 (2)

[0021] The Doppler equation of the main antenna is:

[0022]

[0023] Among them, S iS , S iR (i = 1, 2) respectively correspond to the position vectors of the transmitting and receiving antenna phase centers of the SAR image i. V iS , V iR (i = 1, 2) respectively correspond to the velocity vectors of the transmitting and receiving antenna phase centers of the SAR image i. V P and P respectively correspond to the velocity and position vectors of the ground target point. In the geocentric fixed coordinate system, the ground target point is fixed, and V P is a zero vector. r i(i = 1, 2) corresponds to the slant range of SAR image i, λ is the wavelength, f dc1 corresponds to the Doppler center frequency during the imaging of the SAR image of the main antenna, and φ is the absolute interference phase. i = 1 represents the main image, and i = 2 represents the auxiliary image.

[0024] Step two, perform calibration and correction on the generated DEM data to eliminate phase error and slant range error, and obtain the calibrated DEM data;

[0025] The specific process of calibration and correction is as follows:

[0026] First, based on the sensitivity matrix, use the deviation between the elevation value and the control point to perform calibration processing, and correct the phase error and slant range error;

[0027] Then, repeat the elevation inversion step according to the corrected phase and slant range to obtain the calibrated DEM data.

[0028] Step three, calculate the coherence coefficient of the main and auxiliary images at each angle after registration to obtain a coherence coefficient map for subsequent calculation of the fusion weight;

[0029] The calculation method of the coherence coefficient of any two SAR complex images s1 and s2 is as follows:

[0030]

[0031] Among them, M and N respectively represent the length and width of the estimation window.

[0032] Step four, based on the joint detection of local frequency and eigenvalue, detect the layover shadow area at each angle to obtain the layover shadow detection value;

[0033] The local frequency analysis method uses the interference phase feature of the reversed interference phase in the layover shadow area for detection, and the eigenvalue decomposition method is to detect by estimating the number of signal sources. Combine these two methods for layover shadow area detection, specifically:

[0034] Step 401, local frequency estimation: Use the two-dimensional fast Fourier transform and Chirp-z transform to calculate the local frequency, and combine the amplitude threshold segmentation method to distinguish the results of local frequency analysis to obtain the shadow area D1, the layover area L1, and the normal area N;

[0035] Step 402, eigenvalue decomposition: Perform eigenvalue decomposition on the covariance matrix of the signal model, and arrange the decomposed eigenvalues in descending order to obtain the first largest eigenvalue and the second largest eigenvalue of the target; then combine the detection results of local frequency analysis to set a judgment threshold, and perform threshold division on the first largest eigenvalue and the second largest eigenvalue of the target point by point to obtain the detection result of eigenvalue decomposition.

[0036] Define the first and second largest eigenvalues of the pixel point (m, n) as λ1(m, n) and λ2(m, n) respectively, and calculate the areas occupied by these two largest eigenvalues in the three regions D1, L1, and N as thresholds, then the threshold σ D , σ L and σ N1 , σ N2 are expressed as:

[0037]

[0038] Among them, the num(·) operation represents the number of elements that satisfy the conditions in the parentheses.

[0039] Step 403, Joint detection: Use the detection results of local frequency and eigenvalues for joint detection. Traverse each pixel point. When the detection results of the two methods are consistent, directly output the result; if the detection results of the two methods are inconsistent, then make a further judgment:

[0040] If one of the detection results is an aliasing area and the other result is a normal area, then determine whether the pixel point satisfies λ2(m, n) > (σ N2 + σ L ) / 2. If so, mark the pixel point as an aliasing area; otherwise, mark it as a normal area;

[0041] If one of the detection results is a shadow area and the other result is a normal area, then determine whether the pixel point satisfies λ1(m, n) < (σ N1 + σ D ) / 2. If so, mark the pixel point as a shadow area; otherwise, mark it as a normal area.

[0042] Step 404, Obtain the aliasing and shadow detection value according to the joint detection result. The value of the normal area is 1, and the value of the aliasing or shadow area is 0.

[0043] Step Five, Combine the coherence coefficient map and the aliasing and shadow detection value to set the fusion weight of each pixel point in the interferometric SAR data at each angle, and obtain the weight map;

[0044] For the interferometric SAR data at different angles, perform weighted fusion on the pixel points in the same area at the elevation level. The specific process is as follows:

[0045] First, average the aliasing and shadow detection values. For any pixel point, take a window centered on this pixel point with a size of M×N, and average the aliasing and shadow detection values within the window to obtain the average aliasing and shadow detection value of each pixel point. The calculation is as follows:

[0046]

[0047] Among them, dk Represents the joint detection value at the kth angle. The detected overlapping shadow area value is 0, and the pixel values in the remaining area are 1.

[0048] Then, the fusion weight is obtained by combining the coherence coefficient of the corresponding pixel points. The calculation method is as follows:

[0049]

[0050] Among them, γ k represents the coherence coefficient of the kth interferometric SAR data, Md k Represents the mean value of overlapping shadow detection of the kth interferometric SAR data.

[0051] Step 6: Project the weight map and the corrected DEM values into a unified coordinate system, and remove the coordinates of abnormal points;

[0052] Specifically:

[0053] Firstly, according to the longitude and latitude of each pixel point obtained after correction, the common area of the multi-angle interferometric SAR data is selected for projection, and the longitude and latitude intervals of the projection network are calculated to construct the projection network under the unified longitude and latitude coordinate system.

[0054] Then, the elevation value and weight map are projected into a unified longitude and latitude coordinate system, and abnormal pixels with a weight of 0 are removed.

[0055] Step seven, based on the fusion weights of the coherence coefficient and the mean value of the overlap shadow detection, the projected multi-angle interferometric SAR data are fused, and interpolation is performed in the null value area to obtain the fused DEM, thus completing the 3D reconstruction of the area to be detected.

[0056] The fused elevation is expressed as follows:

[0057]

[0058] Among them, h k Indicates the corrected elevation value.

[0059] The elevation of a pixel after fusion is divided into the following cases:

[0060] 1) There are two or more normal pixels in multiple interferometric SAR data, and the fused elevation is equal to the result of weighted fusion based on the coherence coefficient and the mean of overlapping shadow detection.

[0061] 2) There is only one normal pixel in the multiple interferometric SAR data, and the elevation after fusion is equal to the corresponding elevation value of the pixel;

[0062] 3) There are no normal pixels in multiple interferometric SAR data, and the fused elevation is equal to the interpolation of the surrounding fused elevation values.

[0063] The advantages of the present invention are as follows:

[0064] (1) Practicality. The method proposed by the present invention can solve the problem of 3D reconstruction in complex terrain areas. By fusing multi-angle InSAR data in different directions, the elevation accuracy is improved, and the invalid data caused by layover shadow areas is reduced.

[0065] (2) Versatility. The method proposed by the present invention can not only be used for the multi-angle InSAR 3D reconstruction in the embodiment, but also be well migrated to the multi-source InSAR 3D reconstruction in other scenarios. Description of the Drawings

[0066] Figure 1 is the overall flowchart of a multi-angle InSAR 3D reconstruction method based on layover shadow detection according to the present invention;

[0067] Figure 2 is the elevation inversion value after correction of two observation directions in the embodiment of the present invention;

[0068] Figure 3 is the coherence coefficient map of two observation directions in the embodiment of the present invention;

[0069] Figure 4 is the joint detection result of two observation directions in the embodiment of the present invention;

[0070] Figure 5 is the DEM fusion result based on the coherence coefficient in the embodiment of the present invention;

[0071] Figure 6 is the DEM fusion result based on the elevation error map in the embodiment of the present invention;

[0072] Figure 7 is the DEM fusion result of the method of the present invention in the embodiment of the present invention. Detailed Embodiment

[0073] The present invention will be further described in detail below with reference to the drawings and specific examples.

[0074] For the problem of 3D reconstruction of complex terrain, the present invention proposes a multi-angle interferometric SAR 3D reconstruction method based on layover shadow. Based on multi-angle interferometric SAR data in different directions, using the complementarity of their layover shadow areas, the multi-angle interferometric SAR data is fused to improve the elevation accuracy and the quality and reliability of DEM data in complex terrain areas. In addition, weights are set by combining the layover shadow detection mean and the coherence coefficient, fully considering the layover and shadow problems caused by the side-looking imaging characteristics of SAR, so as to reduce the missing and error data brought by the layover shadow area. Using the multi-angle interferometric SAR 3D reconstruction method based on layover shadow detection proposed by the present invention, the 3D coordinates of the multi-angle interferometric SAR are solved, and the implementation process of the proposed multi-angle interferometric SAR 3D reconstruction method in specific applications is described in detail.

[0075] A multi-angle interferometric SAR 3D reconstruction method based on layover shadow detection, the implementation process is as Figure 1 shown, specifically including the following steps:

[0076] Step 1: Perform interferometric processing on the interferometric SAR data irradiated at each angle respectively to generate corresponding DEM data;

[0077] The specific process is:

[0078] Register the single-look complex image pairs of the interferometric SAR irradiated at each angle respectively, and then conjugate multiply to obtain the corresponding interferometric phase diagram; then use the interferometric spectrum shift to remove the flat-earth effect, and use the non-local filtering method based on non-linear fringe compensation to filter out the noise phase; use the statistical cost flow method for phase unwrapping to improve the unwrapping efficiency of the interferometric diagram in complex terrain; finally, perform elevation inversion, and use the Newton iteration method to jointly solve the 3D coordinates of the target point by combining the interferometric phase equation, the slant range equation and the Doppler equation to obtain the DEM data.

[0079] According to the system geometric model, the interferometric phase equation can be listed as:

[0080]

[0081] The slant range equation of the SAR main image is:

[0082] (|S 1S -P| + |S 1R -P|) = 2r1 (11)

[0083] The Doppler equation of the main antenna is:

[0084]

[0085] Among them, S iS , S iR(i = 1, 2) respectively correspond to the position vectors of the transmitting and receiving antenna phase centers of SAR image i. V iS , V iR (i = 1, 2) respectively correspond to the velocity vectors of the transmitting and receiving antenna phase centers of SAR image i. V P And P respectively correspond to the velocity and position vectors of the ground target point. In the geocentric fixed coordinate system, the ground target point is stationary, and V P is a zero vector. r i (i = 1, 2) respectively correspond to the slant ranges of SAR image i, λ is the wavelength, f dc1 corresponds to the Doppler center frequency during the imaging of the SAR image of the main antenna, and φ is the absolute interferometric phase. i = 1 represents the main image, and i = 2 represents the auxiliary image.

[0086] Step 2: Calibrate and correct the generated DEM data to eliminate phase error and slant range error;

[0087] The specific process is as follows:

[0088] First, based on the sensitivity matrix, perform calibration processing using the deviation between the elevation value obtained by inversion and the control points, and correct the phase error and slant range error.

[0089] Then, combine the unwrapped phase and slant range after correction and repeat the elevation inversion operation in Step 1 to obtain the corrected three-dimensional coordinates.

[0090] When L calibrators are deployed in the survey strip and the elevation of each calibrator is known as h i , i = 1, 2, …, L, the difference between the reconstructed elevation and the actual elevation of each calibration point can be obtained, and the matrix equation can be expressed as

[0091] Δ = F·ΔX + M (13)

[0092] Where, is the elevation error at the i-th calibration point; M is the linearized error matrix; is the sensitivity matrix.

[0093] The specific form of the sensitivity matrix is as follows:

[0094]

[0095] Where represents the sensitivity of the ground elevation f to the interferometric parameter X at the i-th calibrator. The interferometric parameter X includes the main antenna height H, the main antenna slant range r1, the baseline length B, the baseline inclination angle θ b , and the interferometric phase φ. Considering the Earth's curvature, it is specifically expressed as follows:

[0096]

[0097] where R e is the local Earth radius, h is the height from the ground target at the calibrator to the Earth ellipsoid surface, θ is the main antenna viewing angle at the calibrator, α is the baseline inclination angle before calibration, r is the slant range from the main antenna to the calibrator, and Δr is the slant range difference between the main and auxiliary antennas to the ground calibrator.

[0098] The interferometric parameter error can be expressed as: ΔX = F + ·Δ. Where F + is the generalized inverse of the sensitivity matrix F. Thus, the interferometric calibration of the spaceborne InSAR system can be achieved.

[0099] Step 3: Calculate the coherence coefficient at each angle;

[0100] For two SAR complex images s1 and s2, the calculation method is as follows:

[0101]

[0102] where M and N represent the length and width of the estimation window respectively.

[0103] Step 4: Jointly detect the layover and shadow regions at each angle based on local frequency and eigenvalue;

[0104] The local frequency analysis method uses the interferometric phase feature of the reversed interferometric phase in the layover and shadow regions for detection; the eigenvalue decomposition method is to detect by estimating the number of signal sources.

[0105] The process of the joint detection method based on local frequency and eigenvalue is as follows:

[0106] First, use the two-dimensional fast Fourier transform and Chirp-z transform to improve the estimation accuracy of local frequency at a relatively low computational cost, and combine the amplitude threshold segmentation method to distinguish the results of local frequency analysis, namely the shadow region D1, the layover region L1, and the normal region N;

[0107] Then, perform eigenvalue decomposition on the covariance matrix of the signal model, and arrange the decomposed eigenvalues in descending order to obtain the first and second largest eigenvalues of the target; then set the joint judgment threshold with the detection result of local frequency analysis as the prior information, and perform threshold division on the first and second largest eigenvalues of the target point by point to obtain the detection result of eigenvalue decomposition;

[0108] Define the first and second largest eigenvalues of the pixel point (m, n) as λ1(m, n) and λ2(m, n) respectively, and the threshold can be calculated by the area occupied by these two largest eigenvalues in the three regions D1, L1, and N. Where the threshold σ D, σ L and σ N1 , σ N2 can be expressed as:

[0109]

[0110] where the num(·) operation represents the number of elements that satisfy the condition in the parentheses.

[0111] Finally, the layover shadow area of the target is jointly detected using the local frequency detection result and the eigenvalue detection result: traverse each pixel point. When the detection results of the two methods are consistent, the result is directly output; if the detection results of the two methods are inconsistent, there are the following four possible situations:

[0112] 1) The pixel is detected as a normal area in the local frequency analysis method, but as a layover area in the eigenvalue decomposition method: if the pixel point λ2(m,n) > (σ N2 + σ L ) / 2, it is marked as a layover area; otherwise, it is marked as a normal area;

[0113] 2) The pixel is detected as a normal area in the local frequency analysis method, but as a shadow area in the eigenvalue decomposition method: if the pixel point λ1(m,n) < (σ N1 + σ D ) / 2, it is marked as a shadow area; otherwise, it is marked as a normal area;

[0114] 3) The pixel is detected as a layover area in the local frequency analysis method, but as a normal area in the eigenvalue decomposition method: if the pixel point λ2(m,n) > (σ N2 + σ L ) / 2, it is marked as a layover area; otherwise, it is marked as a normal area;

[0115] 4) The pixel is detected as a shadow area in the local frequency analysis method, but as a normal area in the eigenvalue decomposition method: if the pixel point λ1(m,n) < (σ N1 + σ D ) / 2, it is marked as a shadow area; otherwise, it is marked as a normal area.

[0116] Step 5: Combine the coherence coefficient map and the joint detection value to set the fusion weight of each pixel point in the interferometric SAR data at each angle, and obtain the weight map;

[0117] For the interferometric SAR data at different angles, due to the different positions of the layover shadow areas, the elevation values obtained by the pixel points in the same area are different, and weighted fusion needs to be performed at the elevation level.

[0118] For any pixel, the mean value of the joint detection in a window of a certain size with the pixel as the center is taken as the judgment weight of the overlapping shadow area of the pixel. Combined with the coherence coefficient of the corresponding pixel, the calculation method of the fusion weight is as follows:

[0119]

[0120] Among them, γ k represents the coherence coefficient at the kth angle, Md k Represents the joint detection mean of the overlapping shadow area at the kth angle:

[0121]

[0122] Among them, d k Represents the joint detection value at the kth angle. The detected overlapping shadow area value is 0, and the pixel value of the remaining area is 1. M×N represents the selected mean window size.

[0123] Averaging the overlapping shadow detection values can reduce the error caused by the overlapping shadow detection method.

[0124] Step 6: Project the weight map of each angle and the corrected DEM value into a unified coordinate system, and remove the coordinates of abnormal points, where the coordinates of abnormal points are pixels with zero fusion weight;

[0125] According to the longitude and latitude of each pixel point obtained after correction, the common area of the multi-angle interferometric SAR data is selected for projection, and the longitude and latitude intervals of the projection network are calculated to construct a projection network under a unified longitude and latitude coordinate system; the elevation value and weight map are projected to a unified longitude and latitude coordinate system, and the coordinates of abnormal points with a weight of 0 are eliminated.

[0126] Step 7: Based on the weight of the coherence coefficient and the joint detection mean, the projected multi-angle interferometric SAR data are fused and interpolated in the null value area to obtain the fused DEM.

[0127] The fused elevation can be expressed as follows:

[0128]

[0129] For the elevation of a pixel after fusion, there may be the following situations:

[0130] 1) There are two or more normal pixels in multiple interferometric SAR data, and the fused elevation is equal to the result of weighted fusion based on the coherence coefficient and the mean of overlapping shadow detection.

[0131] 2) There is only one normal pixel in the multiple interferometric SAR data, and the elevation after fusion is equal to the corresponding elevation value of the pixel;

[0132] 3) There are no normal pixel points in multiple interferometric SAR data, and the elevation after fusion is equal to the interpolation of the elevation values after fusion of the surrounding areas.

[0133] Experimental results

[0134] To illustrate the effectiveness of the present invention, this embodiment conducts experimental verification on multi-angle interferometric SAR data in steep mountainous areas.

[0135] The real data used is multi-angle observation InSAR data in a large-slope complex terrain area in the mountainous area of western Sichuan, with two observation directions: northward and southward. The system parameters of the two observation directions are shown in Table 1.

[0136] Table 1 System parameters of southward and northward data

[0137] Parameter name Northward observation value Southward observation value Center frequency 30 GHz 30 GHz Baseline length 0.5m 0.5m Baseline inclination 45° 45° Carrier altitude 4500m 4510m

[0138] The elevation values obtained by inverting the elevation of the data in these two observation directions using the method of the present invention are as Figure 2 shown. Due to the influence of layover shadows, noise, etc., there are obvious abnormal areas in the images. The coherence coefficient maps corresponding to the interferometric SAR data in these two observation directions are as Figure 3 shown. The layover shadow detection results are as Figure 4 shown, where the white areas represent layover and the black areas represent shadows.

[0139] Corresponding to the set coded row (north-south direction) resolution of 0.2 m and column (east-west direction) resolution of 0.2 m, geocoding is performed on the above three images respectively. To reflect the superiority of the method of the present invention, the fusion results obtained by the method of the present invention are respectively compared with the existing weighted fusion method. As Figure 5 shows the weighted fusion result based on the coherence coefficient, where the pixel points with a coherence coefficient lower than 0.8 are excluded as abnormal data, Figure 6 shows the weighted fusion result based on the elevation error map, Figure 7 shows the weighted fusion result of the method of the present invention. It can be seen from the fusion results that the result of the method of the present invention is smoother than the other two methods.

[0140] To quantitatively evaluate the quality of different DEM products, the root mean square error and maximum elevation error of the elevation before and after fusion under the three methods are calculated respectively, as shown in Table 2. Among them, the abnormal points refer to the data points with a deviation value from the surrounding pixel points greater than the minimum ambiguity height. It can be seen from the table that the method of the present invention has the smallest root mean square error, maximum elevation error and number of abnormal points, proving the superiority of the method of the present invention.

[0141] Table 2 DEM accuracy evaluation

[0142]

[0143]

Claims

1. A multi - angle interferometric SAR three - dimensional reconstruction method based on layover shadow detection, characterized in that, The following steps are involved: Step 1: Synthetic aperture radar collects multi-angle imaging information of the area to be detected, performs interference processing on the interferometric SAR data illuminated at each angle, and generates DEM data corresponding to each angle; Step 2: calibrate and correct the generated DEM data to eliminate phase error and slant range error to obtain corrected DEM data; The specific process of calibration and correction is as follows: Firstly, based on the sensitivity matrix, the deviation between the elevation value and the control point is used for calibration, and the phase error and the slant range error are corrected. Then, the elevation inversion is repeated according to the corrected phase and slant range to obtain the corrected DEM data; Step 3: Calculate the coherence coefficient of the primary and secondary images at each angle after registration to obtain a coherence coefficient map for subsequent calculation of fusion weights; The coherence coefficient of any two SAR complex images s1 and s2 is calculated as follows: Among them, M and N represent the length and width of the estimation window respectively; Step 4: Based on the local frequency and the eigenvalue, the overlapping shadow area at each angle is jointly detected to obtain the overlapping shadow detection value; The local frequency analysis method uses the interference phase characteristics of the interference phase reversal in the overlapping shadow area for detection. The eigenvalue decomposition method uses the estimated number of signal sources for detection. The two methods are combined to detect the overlapping shadow area, specifically: Step 401, local frequency estimation: using two-dimensional fast Fourier transform and Chirp-z transform to calculate the local frequency, and combining the amplitude threshold segmentation method to distinguish the results of local frequency analysis, to obtain the shadow area D1, the overlapped area L1 and the normal area N; Step 402, eigenvalue decomposition: perform eigenvalue decomposition on the covariance matrix of the signal model, and arrange the decomposed eigenvalues in descending order to obtain the first largest eigenvalue and the second largest eigenvalue of the target; then set a judgment threshold in combination with the detection result of the local frequency analysis, perform threshold division on the first largest eigenvalue and the second largest eigenvalue of the target point by point, and obtain the detection result of the eigenvalue decomposition; Define the first and second largest eigenvalues of the pixel point (m, n) as λ1(m, n) and λ2(m, n) respectively, and calculate the areas occupied by these two largest eigenvalues in the three regions D1, L1, and N as thresholds, then the thresholds σ D , σ L and σ N1 , σ N2 are expressed as: The num(·) operation represents the number of elements that satisfy the conditions in the brackets; Step 403, joint detection: performing joint detection using the detection results of the local frequency and the eigenvalue to obtain detection results of the overlapped shadow area and the normal area; Step 404, obtaining an overlap shadow detection value according to the joint detection result, wherein the value of the normal area is 1, and the value of the overlap shadow area is 0; Step 5, combining the coherence coefficient map and the overlap shadow detection value to set the fusion weight of each pixel in the interferometric SAR data at each angle to obtain a weight map; The calculation method of fusion weight is as follows: where γ k represents the coherence coefficient of the k-th interferometric SAR data, and Md k represents the mean value of layover shadow detection of the k-th interferometric SAR data; Step 6: Project the weight map and the corrected DEM values into a unified coordinate system, and remove the coordinates of abnormal points; Specifically: Firstly, according to the longitude and latitude of each pixel point obtained after correction, the common area of the multi-angle interferometric SAR data is selected for projection, and the longitude and latitude intervals of the projection network are calculated to construct the projection network under the unified longitude and latitude coordinate system. Then, the elevation value and weight map are projected into a unified longitude and latitude coordinate system, and abnormal pixels with a weight of 0 are removed; Step 7: Based on the fusion weights of the coherence coefficient and the mean of the layover and shadow detection, fuse the projected multi-angle interferometric SAR data, and perform interpolation in the null value area to obtain the fused DEM, completing the 3D reconstruction of the area to be detected; The elevation after fusion is expressed as follows: Among them, h k represents the corrected elevation value.

2. The multi-angle interferometric SAR three-dimensional reconstruction method based on layover shadow detection according to claim 1, wherein The specific process of performing interferometric processing on the interferometric SAR data of each angle is as follows: 1) Respectively register the single-look complex image pairs of the interferometric SAR irradiated at each angle, and then conjugate multiply to obtain the corresponding interferometric phase map; 2) Perform flat-earth removal on the interferometric phase map based on the interferometric spectrum shift, and then use a non-local filtering method based on non-linear fringe compensation to filter out the noise phase to obtain the filtered interferometric phase map; 3) Perform phase unwrapping on the filtered interferometric phase map using the statistical cost flow method; 4) Elevation inversion: Based on the unwrapped phase, use the Newton iteration method, and combine the interferometric phase equation, the slant range equation, and the Doppler equation to jointly solve the three-dimensional coordinates of the target point to obtain the DEM data.

3. A multi-angle interferometric SAR three-dimensional reconstruction method based on layover shadow detection according to claim 2, characterized in that The interferometric phase equation, the slant range equation, and the Doppler equation are respectively: The interferometric phase equation is: The slant range equation of the SAR main image is: (|S 1S -P|+|S 1R -P|) = 2r1 (6) The Doppler equation of the main antenna is: Among them, S iS , S iR (i = 1, 2) respectively correspond to the position vectors of the transmitting and receiving antenna phase centers of SAR image i; V iS , V iR (i = 1, 2) respectively correspond to the velocity vectors of the transmitting and receiving antenna phase centers of SAR image i; V P and P respectively correspond to the velocity and position vectors of the ground target point. In the geocentric fixed coordinate system, the ground target point is stationary and V P is a zero vector; r i (i = 1, 2) respectively correspond to the slant ranges of SAR image i, λ is the wavelength, f dc1 corresponds to the Doppler center frequency during the imaging of the SAR image of the main antenna, φ is the absolute interference phase; i = 1 represents the main image and i = 2 represents the secondary image.

4. A multi - angle interferometric SAR three - dimensional reconstruction method based on layover shadow detection according to claim 1, characterized in that, The sensitivity matrix is: Among them, represents the sensitivity of the ground elevation f at the i-th calibrator to the interference parameter X; the interference parameter X includes the main antenna height H, the main antenna slant range r1, the baseline length B, the baseline inclination angle θ b and the interference phase φ. Considering the case of the Earth's curved surface, it is specifically expressed as follows: where R e is the local Earth radius, h is the height from the ground target at the calibrator to the Earth ellipsoid surface, θ is the main antenna view angle corresponding to the calibrator, α is the baseline inclination angle before correction, r is the slant range from the main antenna to the calibrator, and Δr is the slant range difference between the main and auxiliary antennas to the ground calibrator.

5. A multi-angle interferometric SAR three-dimensional reconstruction method based on layover shadow detection according to claim 1, wherein The specific process of the joint detection of the local frequency and the eigenvalue is: Traverse each pixel point. When the detection results of the two methods are the same, directly output the result; if the detection results of the two methods are different, then make a further judgment: If one of the detection results is an aliasing region and the other result is a normal region, then determine whether the pixel point satisfies λ2(m,n) > (σ N2 + σ L ) / 2. If so, mark the pixel point as an aliasing region; otherwise, mark it as a normal region; If one of the detection results is a shaded area and the other is a normal area, then determine whether the pixel point satisfies λ1(m,n) < (σ N1 + σ D ) / 2. If so, mark the pixel point as a shaded area; otherwise, mark it as a normal area.

6. A multi - angle interferometric SAR three - dimensional reconstruction method based on layover shadow detection according to claim 1, characterized in that, The calculation process of the fusion weight is: First, average the layover and shadow detection values: For any pixel point, take a window centered on this pixel point with a size of M×N, and average the layover and shadow detection values within the window to obtain the mean of the layover and shadow detection values of each pixel point, and its calculation is as follows: where d k represents the joint detection value at the k-th angle, the detected value of the layover shadow area is 0, and the pixel values of the remaining areas are 1. Then, the mean of the layover and shadow detection values of each pixel point is combined with the coherence coefficient of the corresponding pixel point to calculate the fusion weight.

7. A multi - angle interferometric SAR three - dimensional reconstruction method based on layover shadow detection according to claim 1, characterized in that, For the elevation of a certain pixel point in the fused DEM, it is divided into the following situations: 1) If there are two or more normal pixel points in multiple interferometric SAR data, the elevation after fusion is equal to the result of weighted fusion based on the coherence coefficient and the mean of the layover and shadow detection; 2) If there is only one normal pixel point in multiple interferometric SAR data, the elevation after fusion is equal to the elevation value corresponding to this pixel point; 3) If there are no normal pixel points in multiple interferometric SAR data, the elevation after fusion is equal to the interpolation of the elevation values after fusion of the surrounding areas.

Images