Method, system, equipment and medium for extracting understory terrain based on interferometric radar images

Abstract

The invention discloses a method, system, device and medium for extracting understory terrain based on interferometric radar images. The method includes: obtaining dual-station InSAR data for interferometric processing and correction; performing azimuth sub-aperture decomposition on the original InSAR image to obtain sub-view SAR images under different spectra; performing coherence amplitude correction, complex coherence coefficient reconstruction and spreading points on the complex plane for the complex coherence coefficients of different sub-view SAR images; performing gross error elimination and least square coherence straight line fitting on the spread points within the complex unit circle to obtain the phase information of the two intersections of the coherence straight line and the unit circle; judging one of the two intersections as a surface phase point based on the effective vertical wave number and the phase of the two intersections; recombining the obtained surface phase and the corrected coherence amplitude into a new complex coherence coefficient, and further obtaining the understory terrain elevation information based on this. The understory terrain of the forest area can be extracted without relying on the assistance of other external data products.

Description

Technical Field

[0001] The present invention relates to the field of forest area terrain mapping remote sensing using synthetic aperture radar interferometry (InSAR), and in particular to a method, system, equipment and medium for extracting understory terrain based on interferometric radar images. Background Art

[0002] Topographic mapping of forest areas has always been the focus and difficulty in the fields of surveying and mapping science and technology and remote sensing science and technology. It is related to the collection of basic geographic information and plays a vital role in many fields such as engineering construction and resource exploration. Optical remote sensing methods are limited in their ability to map the understory of forest areas, while synthetic aperture radar interferometry (InSAR) has become a key technology for topographic mapping of forest areas due to its penetrating properties and its unaffected by clouds and fog. To date, most of the spaceborne InSAR systems launched internationally use the re-entry interferometry mode, and the interferometric baseline is not very suitable for topographic mapping. The development of spaceborne dual-station InSAR systems provides an important opportunity for topographic mapping. The dual-station InSAR system, represented by the TanDEM-X launched in 2010 and the Lutan-1 InSAR system launched in 2022, adopts a "one-transmit, two-receive" observation mode, eliminating the influence of temporal decorrelation, ensuring high interference quality and excellent coherence. This type of dual-station InSAR system is designed to serve terrain mapping and therefore has appropriate height sensitivity.

[0003] However, the existing dual-station InSAR systems generally adopt a single-polarization working mode, and can only provide two observation quantities, coherence and interference phase, after interference processing, which cannot take into account the complex physical model solution requirements. In forest areas, due to the penetrability of electromagnetic waves, the interference phase center height obtained by dual-station InSAR processing is neither at the top of the forest canopy nor at the real ground, but generally located between the two, which limits the accuracy of the real digital terrain product (Digital Terrain Model, DTM). The main reason is that the penetration of electromagnetic waves is affected by multiple factors such as wavelength, forest type, and canopy density. The phase center height is generally located between the surface elevation and the top of the canopy, but the specific location cannot be expressed in a unified measurement method. Therefore, it is necessary to resume accurate physical models and statistical methods to finely deduct the residual vegetation height signal in the forest area during traditional InSAR terrain mapping, so as to invert high-precision understory terrain elevation information.

[0004] In order to remove the residual vegetation height signal in the forest area during dual-station InSAR terrain mapping, the following methods are generally used:

[0005] (1) Assisted by measurement data from multiple platforms. This type of method includes: using space-borne lidar data to establish a physical model to remove the residual forest height in the InSAR signal, and using lidar discrete observation points, optical remote sensing data and other multi-source remote sensing data for intelligent learning to invert the understory terrain. This type of method often requires a lot of data and has high data requirements, and cannot meet the needs of understory terrain extraction in universal scenarios.

[0006] (2) Based on the random volume over ground (RVoG) model used in the traditional inversion of forest terrain, multi-polarization observation data is added to form a full polarization observation space. According to the geometric relationship of the complex coherence coefficient of different polarization modes within the unit circle of the complex plane, the coherence line is fitted to obtain the surface phase point, and further obtain the real surface elevation information to complete the extraction of forest terrain. This method requires the provision of full polarization PolInSAR (Polarization Interferometric Synthetic Aperture Radar) data, while the current space-borne dual-station InSAR data mainly uses a single polarization observation mode to collect terrain mapping data, which cannot meet this data requirement.

[0007] In summary, the existing spaceborne dual-station InSAR for extracting understory terrain in forest areas is heavily dependent on the auxiliary or fully polarized InSAR data from other platforms. The current spaceborne dual-station InSAR data not only has the unique advantages of dual-station observations, but also has the advantages of high resolution. How to mine its own observation space and expand InSAR's own observation information to extract understory terrain needs further research. Therefore, under the premise of using limited observation data, it is necessary to design a method for extracting understory terrain based on interferometric radar images. Summary of the invention

[0008] In view of the above-mentioned deficiencies in the prior art, the purpose of the present invention is to provide a method, system, device and medium for extracting understory terrain based on interferometric radar images, so as to realize a method for extracting understory terrain by fully exploiting the observation space of InSAR itself.

[0009] In a first aspect, a method for extracting understory terrain based on interferometric radar images is provided, comprising the following steps:

[0010] S1: Obtain the dual-station InSAR image data of the forest terrain area to be inverted, perform interferometric processing, obtain the initial InSAR data coherence, perform coherence correction on it to obtain the coherence correction factor, and obtain the corrected coherence amplitude;

[0011] S2: Perform azimuth spectrum filtering on the original InSAR image to achieve sub-aperture decomposition and obtain sub-view SAR images under different spectra;

[0012] S3: respectively obtain the complex coherence coefficients of the sub-view SAR images under different spectra, and perform coherence amplitude correction, complex coherence coefficient reconstruction and point expansion on the complex plane for the complex coherence coefficients under different spectra;

[0013] S4: For any pixel, the points within the complex unit circle are subjected to gross error elimination and least square fitting of the coherence line, the slope and intercept of the coherence line are calculated, and the phase information of the two intersection points of the coherence line and the complex unit circle is further obtained;

[0014] S5: Calculate the effective vertical wave number through the baseline parameter file, and determine one of the two intersection points as the surface phase point according to the positive and negative values ​​of the effective vertical wave number and the phase magnitude relationship between the two intersection points obtained above;

[0015] S6: The surface phase information and the coherence amplitude of the initial InSAR data coherence are reorganized into a new complex coherence coefficient, and then de-flattening, de-terraining, phase filtering, unwrapping, de-orbiting and phase height conversion are performed in sequence to finally obtain the understory terrain elevation information of the understory terrain area to be inverted.

[0016] Furthermore, in step S1, the coherence correction factor includes a signal-to-noise ratio decoherence factor and quantization error decoherence factor Quantization error decoherence factor It represents the measure of the decoherence factor caused in the data processing process, and its value is a constant.

[0017] Furthermore, in step S2, in the process of performing azimuth spectrum filtering on the original InSAR image to realize sub-aperture decomposition, the spectrum center of the sub-view SAR image is expressed as:

[0018] ; ;

[0019] In the formula, represents the spectrum center of the kth sub-view SAR image; N represents the number of sub-view SAR images split along the azimuth direction; Indicates the split bandwidth.

[0020] Furthermore, step S3 specifically includes:

[0021] Obtain the complex coherence coefficient of sub-view SAR images under different spectra , k=1,2,…,N, where N represents the number of sub-view SAR images segmented along the azimuth direction;

[0022] The coherence amplitude correction is performed on the complex coherence coefficients under different spectra. At this time, the coherence correction factor obtained in step S1 is substituted into the following expression for coherence correction:

[0023] ;

[0024] In the formula, the coherence correction factor includes the signal-to-noise ratio decoherence factor and quantization error decoherence factor ; represents the pure body decoherence factor in the interferometric process of SAR images with different sub-views; Indicates the amplitude;

[0025] The reconstruction expression of complex coherence coefficient of SAR images with different sub-views is:

[0026] ;

[0027] In the formula, represents the complex coherence coefficient of the reconstructed k-th sub-view SAR image; It means taking the phase of a complex number; i means the imaginary part;

[0028] The complex coherence coefficient of N reconstructed SAR images with different sub-views Expand the points to the complex plane; within the complex unit circle of the complex plane, any complex coherence coefficient The phase range is , the amplitude range is .

[0029] Furthermore, step S4 specifically includes:

[0030] For a given set of N points within the complex unit circle on the complex plane, the real and imaginary parts are , , centralize the data:

[0031] ;

[0032] ;

[0033] Where N represents the number of sub-view SAR images segmented along the azimuth direction; , denote the mean of the real part and the mean of the imaginary part respectively; , Respectively represent the values ​​of the real and imaginary parts of the k-th point after data centering;

[0034] Perform the least squares straight line fitting on N points, calculate the parameters of the fitting line, and introduce the following auxiliary quantities:

[0035] ;

[0036] In the formula, , and There are three auxiliary quantities;

[0037] Solve to obtain two candidate values ​​for the slope of the fitted line:

[0038] ;

[0039] In the formula, and are two candidate slope values;

[0040] The intercepts corresponding to the above two slopes are:

[0041] ;

[0042] In the formula, and is the intercept corresponding to the above two slopes;

[0043] For outliers, the gross errors are eliminated and they are not used as candidate points for fitting the straight line. Specifically, the residual from the point to the straight line is defined as the vertical distance, and the and Respectively represent the residuals from the point to the two fitted coherent lines:

[0044] Calculate the total residuals from all points to the two fitted coherence lines separately and ;

[0045] Select the fitted coherence line with the smallest total residual as the final fitted coherence line;

[0046] Calculate the mean residual of all points and standard deviation ;

[0047] Suppose the elimination residual exceeds the threshold , the residual of the removed point is greater than All points of

[0048] Refit the coherence line based on the retained points and solve its two intersection points with the complex unit circle and ;

[0049] The phases corresponding to the two intersection points are extracted as follows:

[0050] ;

[0051] In the formula, and are the phases corresponding to the two intersection points respectively; represents the phase of a complex number; i represents the imaginary part.

[0052] Further, in step S5, one of the two intersection points is determined to be a surface phase point based on the positive and negative values ​​of the effective vertical wave number and the phase magnitude relationship between the two intersection points obtained above, and the judgment criterion is:

[0053] ;

[0054] In the formula, represents the effective vertical wave number; It means taking the phase of a complex number; and They represent the phases corresponding to the two intersection points respectively; Indicates the surface phase; express and Conjugate multiplication can be understood as and The phase difference between two phases is obtained; i represents the imaginary part.

[0055] Furthermore, in step S6, the surface phase information and the coherence amplitude of the initial InSAR data coherence after correction are reorganized into a new complex coherence coefficient, which is expressed as follows:

[0056] ;

[0057] In the formula, represents the new complex coherence coefficient; represents the corrected coherence amplitude obtained in step S1; Indicates the surface phase; represents the phase of a complex number; i represents the imaginary part.

[0058] In the second aspect, a forest understory terrain extraction system based on interferometric radar images is provided, comprising:

[0059] The coherence amplitude correction module is used to obtain the dual-station InSAR image data of the forest terrain area to be inverted, and perform interference processing to obtain the initial InSAR data coherence, and perform coherence correction on it to obtain the coherence correction factor, and obtain the corrected coherence amplitude;

[0060] The sub-aperture decomposition module is used to perform azimuth spectrum filtering on the original InSAR image to achieve sub-aperture decomposition and obtain sub-view SAR images under different spectra;

[0061] The sub-view SAR image processing module is used to obtain the complex coherence coefficients of the sub-view SAR images under different spectra, and perform coherence amplitude correction, complex coherence coefficient reconstruction and point expansion on the complex plane for the complex coherence coefficients under different spectra in turn;

[0062] The intersection phase acquisition module is used to remove gross errors and perform least squares fitting of the coherent line on the points within the complex unit circle for any pixel, calculate the slope and intercept of the coherent line, and further obtain the phase information of the two intersection points of the coherent line and the complex unit circle;

[0063] The surface phase point determination module is used to calculate the effective vertical wave number through the baseline parameter file, and determine one of the two intersection points as the surface phase point according to the positive and negative values ​​of the effective vertical wave number and the phase magnitude relationship between the two intersection points obtained above;

[0064] The understory terrain elevation extraction module is used to reorganize the surface phase information and the coherence amplitude of the initial InSAR data coherence after correction into a new complex coherence coefficient, and then perform de-flattening, de-terraining, phase filtering, unwrapping, de-orbiting and phase height conversion in sequence, and finally obtain the understory terrain elevation information of the understory terrain area to be inverted.

[0065] In a third aspect, an electronic device is provided, including:

[0066] Memory on which computer programs or instructions are stored;

[0067] The processor is used to load and execute the computer program or instructions to implement the above-mentioned method for extracting understory terrain based on interferometric radar images.

[0068] In a fourth aspect, a computer-readable storage medium is provided, on which a computer program or instruction is stored. When the computer program or instruction is executed by a processor, the method for extracting understory terrain based on interferometric radar images as described above is implemented.

[0069] The present invention proposes a method, system, device and medium for extracting understory terrain based on interferometric radar images, and constructs a technical solution for extracting understory terrain in forest areas based only on satellite-borne dual-station InSAR data, which gets rid of the dependence of existing methods on other observation data, and achieves high-precision understory terrain extraction in forest areas by expanding InSAR's own observation space. At the data level, the present invention can achieve high-precision understory terrain inversion in forest areas without relying on observation data from other platforms; in terms of accuracy, compared with existing methods, the present invention has higher accuracy and is a robust high-precision understory terrain inversion solution. BRIEF DESCRIPTION OF THE DRAWINGS

[0070] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0071] Figure 1 A flow chart of a method for extracting understory terrain based on interferometric radar images provided in an embodiment of the present invention;

[0072] Figure 2 The location of the verification area and the coverage of the InSAR image provided by the disclosed embodiment of the present invention;

[0073] Figure 3 Schematic diagram of coherent straight line fitting provided in the embodiment disclosed in the present invention, (a) is a schematic diagram of the initial coherent straight line fitting, (b) is a schematic diagram of the coherent straight line fitting after removing the gross error points;

[0074] Figure 4 The result diagram of understory terrain extraction provided by the disclosed embodiment of the present invention; wherein, (a) represents airborne LiDARDTM, (b) represents traditional InSAR DEM, (c) represents understory terrain obtained by the technical method of the present invention, (d) represents airborne LiDARCHM, (e) represents the difference diagram of traditional InSAR DEM and airborne LiDAR DTM, and (f) represents the difference diagram of understory terrain inversion result of the technical method of the present invention and airborne LiDAR DTM;

[0075] Figure 5 The statistical histograms of the difference between the traditional InSAR DEM provided in the disclosed embodiment of the present invention and the understory terrain obtained by the example of the present invention and the LiDAR DTM, wherein (a) represents the probability density statistical histogram of the difference between the traditional InSAR DEM and the airborne LiDAR DTM, and (b) represents the probability density statistical histogram of the difference between the understory terrain inversion result of the technical method of the present invention and the airborne LiDAR DTM;

[0076] Figure 6 The scatter plots of the accuracy verification of the traditional InSAR DEM provided in the disclosed embodiments of the present invention and the understory terrain inversion results obtained by the examples of the present invention and the airborne LiDAR DTM are shown in FIG. The colors of the statistical scatter points represent the corresponding number of pixels, wherein (a) represents the scatter plot of the accuracy verification of the traditional InSAR DEM and the airborne LiDAR DTM, and (b) represents the scatter plot of the accuracy verification of the understory terrain inversion results obtained by the technical method of the present invention and the airborne LiDAR DTM. DETAILED DESCRIPTION

[0077] To make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention will be described in detail below. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other implementation methods obtained by ordinary technicians in this field without creative work belong to the scope of protection of the present invention.

[0078] In order to better illustrate the technical solution of the present invention, taking a test area as an example, a scene of X-band dual-station InSAR data is selected, and airborne LiDAR DTM data is collected as the true value to participate in the verification, and the present invention is further described in detail; it is noted that the specific implementation described here is only used to explain the present invention and is not used to limit the present invention.

[0079] The test area location map and the coverage of the X-band InSAR interferometer and airborne LiDAR DTM are shown in the attached figure. Figure 2 As shown, among the data of the test area selected by the present invention, the airborne LiDAR data is of high precision and can be used as the true value to participate in the result verification of the technical method of the present invention. In the experimental area, the terrain elevation range is 100-400m, the forest coverage rate of the test area exceeds 60%, the average forest height is about 15m, mainly concentrated in 5-25m, and the forest type is mainly northern coniferous forest.

[0080] The basic information of the X-band InSAR interferometer used in the experiment carried out in the example of the present invention is shown in Table 1 below:

[0081] ;

[0082] The embodiment of the present invention provides a method for extracting understory terrain based on interferometric radar images. Figure 1 As shown, the following steps are included:

[0083] S1: Obtain the dual-station InSAR image data of the forest terrain area to be inverted, perform interferometric processing, obtain the initial InSAR data coherence, perform coherence correction on it to obtain the coherence correction factor, and obtain the corrected coherence amplitude.

[0084] The initial InSAR data coherence is expressed as follows:

[0085] (1);

[0086] In the formula, and They represent the main image and auxiliary image involved in the interference, Indicates the conjugate multiplication of the main image and the auxiliary image. Indicates expected value.

[0087] At this point, the coherence of the initial InSAR data obtained needs to be corrected. The coherence correction factor includes the signal-to-noise ratio decoherence factor and quantization error decoherence factor . Signal-to-noise ratio decoherence factor , calculated by the following formula:

[0088] (2);

[0089] (3);

[0090] (4);

[0091] In the formula, and Represent the signal-to-noise ratio of the main image and the auxiliary image, and denote the equivalent zero signal-to-noise ratio of the main image and the auxiliary image, respectively. and Represent the backscattering coefficients of the main image and the auxiliary image respectively.

[0092] Quantization error decoherence factor It represents the decoherence factor caused by data processing, which can be taken as a constant close to 1. Generally speaking, =0.98 can meet the application requirements of most radar images. Of course, other constants close to 1, such as 0.97, 0.99, etc., can also be used in other embodiments.

[0093] S2: Perform azimuth spectrum filtering on the original InSAR image to achieve sub-aperture decomposition, perform azimuth sub-aperture spectrum decomposition on the original InSAR image according to the segmentation bandwidth and the number of segmented sub-view images, and obtain sub-view SAR images under different spectra.

[0094] When observing with synthetic aperture radar, the azimuth incident angle It can be expressed as:

[0095] (5);

[0096] In the formula, represents the Doppler frequency, is the wavelength of the radar electromagnetic wave, is the speed of the SAR sensor. When synthetic aperture radar is used for imaging, the azimuth incident angle The range for full aperture synthesis is: , ], is the range of the azimuth incident angle, and its expression is:

[0097] (6);

[0098] in, is the azimuth resolution, and its expression is:

[0099] (7);

[0100] In the formula, It indicates the Doppler bandwidth. As mentioned above, the Doppler bandwidth of the segmented sub-view SAR image determines the azimuth incident angle and its variation range, thus determining the azimuth resolution.

[0101] According to the above analysis, the number and resolution of the segmented sub-view SAR images are determined by the segmentation bandwidth, so the spectrum center of the segmented sub-view SAR image is They are:

[0102] , (8);

[0103] In the formula, represents the spectrum center of the kth sub-view SAR image; N represents the number of sub-view SAR images split along the azimuth direction; Indicates the split bandwidth.

[0104] In the example provided by the present invention, two spectrum segmentations are performed to obtain a total of 10 sub-view SAR images, and the number of segmented sub-view SAR images N=5 and the segmentation bandwidth are respectively =0.5, the center frequencies of the bandwidth spectrum are: -0.250, -0.125, 0, 0.125, 0.25, and the five sub-views are named sub-view 1, sub-view 2, sub-view 3, sub-view 4, and sub-view 5 respectively; the number of segmented sub-view SAR images N=5, and the segmentation bandwidth is =0.6, the center frequencies of the bandwidth spectrum are: -0.20, -0.10, 0, 0.10, 0.20, and the five sub-views are named sub-view 6, sub-view 7, sub-view 8, sub-view 9, and sub-view 10 respectively.

[0105] S3: Perform multi-view and interference on different sub-view SAR images respectively (the process of multi-view and interference corresponds to the process shown in formula (1)), obtain the complex coherence coefficients of the sub-view SAR images under different spectra, and perform coherence amplitude correction, complex coherence coefficient reconstruction and point expansion on the complex plane for the complex coherence coefficients under different spectra.

[0106] Specifically include: obtaining the complex coherence coefficients of sub-view SAR images under different spectra , k=1,2,…,N, where N represents the number of sub-view SAR images segmented along the azimuth direction. In this embodiment, N=10;

[0107] The coherence amplitude correction is performed on the complex coherence coefficients under different spectra. At this time, the coherence correction factor obtained in step S1 is substituted into the following expression for coherence correction:

[0108] (9);

[0109] In the formula, the coherence correction factor includes the signal-to-noise ratio decoherence factor and quantization error decoherence factor ; represents the pure body decoherence factor in the interferometric process of SAR images with different sub-views; Indicates the amplitude; That is the amplitude of the complex coherence coefficient after correction;

[0110] So far, the complex coherence coefficient reconstruction expression of different sub-view SAR images is:

[0111] (10);

[0112] In the formula, represents the complex coherence coefficient of the reconstructed k-th sub-view SAR image; It means taking the phase of a complex number; i means the imaginary part;

[0113] Then, the complex coherence coefficients of the reconstructed N different sub-view SAR images are Expand the points to the complex plane; within the complex unit circle of the complex plane, any complex coherence coefficient The phase range is , the amplitude range is .

[0114] S4: For any pixel, the points within the complex unit circle are subjected to gross error elimination and least squares fitting of the coherence line, the slope and intercept of the coherence line are calculated, and the phase information of the two intersection points of the coherence line and the complex unit circle is further obtained. In general, the schematic diagram of the coherence line fitting using the sub-view SAR image is as follows: Figure 3 shown.

[0115] For a given set of N points within the complex unit circle on the complex plane, the real and imaginary parts are , , centralize the data:

[0116] (11);

[0117] (12);

[0118] Where N represents the number of sub-view SAR images segmented along the azimuth direction; , denote the mean of the real part and the mean of the imaginary part respectively; , Respectively represent the values ​​of the real and imaginary parts of the k-th point after data centering;

[0119] Perform the least squares straight line fitting on N points, calculate the parameters of the fitting line, and introduce the following auxiliary quantities:

[0120] (13);

[0121] In the formula, , and There are three auxiliary quantities;

[0122] Solve to obtain two candidate values ​​for the slope of the fitted line:

[0123] (14);

[0124] In the formula, and are two candidate slope values;

[0125] The intercepts corresponding to the above two slopes are:

[0126] (15);

[0127] In the formula, and is the intercept corresponding to the above two slopes;

[0128] For outliers with large gross errors, their gross errors are eliminated and they are not used as candidate points for fitting the straight line; specifically, the residual from the defined point to the straight line is the vertical distance:

[0129] (16);

[0130] In the formula, and are the residuals from the point to the two fitted coherence lines;

[0131] Calculate the total residuals from all points to the two fitted coherence lines separately:

[0132] (17);

[0133] In the formula, and is the total residual of all points to two fitted coherence lines;

[0134] The fitted coherence line with the smallest total residual is selected as the final fitted coherence line. The process is as follows:

[0135] (18);

[0136] (19);

[0137] In the formula, b and a are the slope and intercept of the final fitted coherence line respectively;

[0138] Calculate the residual mean and standard deviation of all points;

[0139] (20);

[0140] (twenty one);

[0141] In the formula, and are the residual mean and standard deviation of all points respectively, Represents the residual from the kth point to the final fitted coherence line;

[0142] Suppose the elimination residual exceeds the threshold , the conditions for the points to be removed are:

[0143] (twenty two);

[0144] Keep the relevant points other than the removed points, such as Figure 3 As shown in (a) and (b), when the pixel selected in the example is fitting the coherence line, sub-view 7 (the point corresponding to the 7th sub-view SAR image) and sub-view 8 are outliers in the coherence line obtained by the initial fitting, and there are large gross errors, which affect the quality of the coherence line fitting. After the above gross error judgment and elimination, the coherence line is refitted, and the intersection point of the new coherence line and the unit circle is closer to the real surface phase point, and after being converted into terrain elevation information, it is closer to the real surface elevation;

[0145] Next, find the intersection point between the straight line and the unit circle, refit the coherent straight line based on the retained points, and set the fitted coherent straight line equation to , the unit circle equation is , solve the two intersection points of the fitted coherence line and the complex unit circle as:

[0146] (twenty three);

[0147] Calculate the two intersection points separately Value is , then the two intersection points are and ;

[0148] The phases corresponding to the two intersection points are extracted as follows:

[0149] (twenty four);

[0150] In the formula, and are the phases corresponding to the two intersection points respectively; represents the phase of a complex number; i represents the imaginary part.

[0151] S5: Calculate the effective vertical wave number through the baseline parameter file, and determine that one of the two intersection points is the surface phase point based on the positive and negative values ​​of the effective vertical wave number and the phase size relationship between the two intersection points obtained above.

[0152] Specifically, the effective vertical wave number is calculated first, and the calculation formula is as follows:

[0153] (25);

[0154] In the formula, represents the effective vertical wave number, represents the incident angle of electromagnetic wave, Indicates the radar wave operating wavelength, represents the slope distance, In the example provided by the present invention, the effective vertical wave number The average value is -0.14rad / m. Because the orbit of the satellite-borne InSAR system is high, the incident angle change under the field of view of a single interferometer generally does not exceed 3°, and the corresponding effective vertical wave number There are only minor changes, for example, The range is -0.135 to -0.145.

[0155] At this point, we can judge that one of the two intersection points is the surface phase point by the positive and negative values ​​of the effective vertical wave number and the magnitude relationship between the two intersection points obtained above. The judgment criteria are:

[0156] (26);

[0157] In the formula, represents the effective vertical wave number; It means taking the phase of a complex number; and They represent the phases corresponding to the two intersection points respectively; Indicates the surface phase; express and Conjugate multiplication can be understood as and The phase difference between the two phases is obtained.

[0158] At this point, the complete surface phase is obtained, which can be used to invert the elevation information of the understory terrain in the subsequent steps.

[0159] S6: The surface phase information and the coherence amplitude of the initial InSAR data coherence are reorganized into a new complex coherence coefficient, and then de-flattening, de-terraining, phase filtering, unwrapping, de-orbiting and phase height conversion are performed in sequence to finally obtain the understory terrain elevation information of the understory terrain area to be inverted.

[0160] Among them, the surface phase information and the coherence amplitude after correction of the initial InSAR data coherence are reorganized into a new complex coherence coefficient, which is expressed as follows:

[0161] (27);

[0162] In the formula, represents the new complex coherence coefficient; represents the corrected coherence amplitude obtained in step S1; Indicates the surface phase; represents the phase of a complex number; i represents the imaginary part.

[0163] Then, the reconstructed complex coherence coefficient is successively subjected to deflating, de-terraining, phase filtering, unwrapping, de-orbiting and phase height conversion, and finally the understory terrain elevation information of the forest area is obtained. Figure 4 As shown, Figure 4 (a) is the airborne LiDAR DTM. Figure 4 (b) is the InSAR DEM (Digital Elevation Model) obtained from the original synthetic aperture radar interferometry image. Figure 4 (c) is the forest topography obtained by the technical method of the present invention (also InSAR DEM, for differentiation, the forest topography description is used for differentiation). Because the terrain elevation range is 100-400m, the observation Figure 4 The advantages of the technical method of the present invention cannot be seen from the visual effects in (a)-(c), so the InSAR DEM and the forest terrain obtained by the present invention are differentiated from the LiDAR DTM to obtain Figure 4 (e) and (f) are the difference maps between InSAR DEM and LiDAR DTM, and the difference maps between forest terrain and LiDAR DTM, respectively. Figure 4By comparing the conventional InSAR DEM with the airborne LiDAR canopy height model (ConapyHeight Model, CHM, i.e., forest height) in (d), it can be found that the traditional InSAR DEM contains obvious trend errors in forest height, that is, the higher the forest height, the greater the difference between the InSAR DEM and the LiDAR DTM. In the difference diagram between the understory terrain product obtained by the present invention and the LiDAR DTM, the above trend errors are significantly eliminated, mainly concentrated near 0, which illustrates the effectiveness of the method in this paper. It should be noted that the processes of de-flattening, de-terraining, phase filtering, unwrapping, de-orbiting, and phase height conversion of the reconstructed complex coherence coefficients are all prior art, not the focus of the present invention, and will not be elaborated here.

[0164] In addition, the above two difference graphs are subjected to probability density statistics and presented in the form of probability histograms, as shown in the attached figure. Figure 5 In (a) and (b), it can be seen that the error between the statistical histogram of the difference result of the understory terrain and the LiDAR DTM is significantly eliminated, and the error converges to near 0. Furthermore, the accuracy of the understory terrain obtained by the traditional InSAR DEM and the present invention is statistically evaluated, and the evaluation indicators include the following: (coefficient of determination), RMSE (root mean square error), Std (standard deviation), and Bias (bias), the expressions are:

[0165] (28);

[0166] (29);

[0167] (30);

[0168] (31);

[0169] In the above formula, Indicates the measured value, represents the true value, represents the mean of the true values, represents the average value of the measurement, , n is the number of pixels involved in the statistics. The accuracy verification of the forest terrain obtained by InSAR DEM and the present invention is shown as follows Figure 6 As shown, we can see that the four evaluation indicators have been improved. It increased from 0.9984 to 0.9992, RMSE decreased from 5.07m to 2.09m, STD decreased from 2.95m to 2.09m, and BIAS decreased from 4.12m to -0.02m, which proved that the technical method of the present invention is effective and can improve the inversion accuracy of terrain elevation information in forest areas.

[0170] The embodiment of the present invention further provides a forest terrain extraction system based on interferometric radar images, comprising:

[0171] The coherence amplitude correction module is used to obtain the dual-station InSAR image data of the forest terrain area to be inverted, and perform interference processing to obtain the initial InSAR data coherence, and perform coherence correction on it to obtain the coherence correction factor, and obtain the corrected coherence amplitude;

[0172] The sub-aperture decomposition module is used to perform azimuth spectrum filtering on the original InSAR image to achieve sub-aperture decomposition and obtain sub-view SAR images under different spectra;

[0173] The sub-view SAR image processing module is used to obtain the complex coherence coefficients of the sub-view SAR images under different spectra, and perform coherence amplitude correction, complex coherence coefficient reconstruction and point expansion on the complex plane for the complex coherence coefficients under different spectra in turn;

[0174] The intersection phase acquisition module is used to remove gross errors and perform least squares fitting of the coherent line on the points within the complex unit circle for any pixel, calculate the slope and intercept of the coherent line, and further obtain the phase information of the two intersection points of the coherent line and the complex unit circle;

[0175] The surface phase point determination module is used to calculate the effective vertical wave number through the baseline parameter file, and determine one of the two intersection points as the surface phase point according to the positive and negative values ​​of the effective vertical wave number and the phase magnitude relationship between the two intersection points obtained above;

[0176] The understory terrain elevation extraction module is used to reorganize the surface phase information and the coherence amplitude of the initial InSAR data coherence after correction into a new complex coherence coefficient, and then perform de-flattening, de-terraining, phase filtering, unwrapping, de-orbiting and phase height conversion in sequence, and finally obtain the understory terrain elevation information of the understory terrain area to be inverted.

[0177] It should be understood that the functional unit modules in various embodiments of the present invention may be concentrated in one processing unit, or each unit module may exist physically separately, or two or more unit modules may be integrated in one unit module, and may be implemented in the form of hardware or software.

[0178] An embodiment of the present invention further provides an electronic device, including:

[0179] Memory on which computer programs or instructions are stored;

[0180] The processor is used to load and execute the computer program or instructions to implement the above-mentioned method for extracting understory terrain based on interferometric radar images.

[0181] The embodiment of the present invention also provides a computer-readable storage medium having a computer program or instruction stored thereon. When the computer program or instruction is executed by a processor, the method for extracting understory terrain based on interferometric radar images as described above is implemented.

[0182] It can be understood that the same or similar parts of the above embodiments can be referenced to each other, and the contents not described in detail in some embodiments can refer to the same or similar contents in other embodiments.

[0183] Those skilled in the art will appreciate that the embodiments of the present application may be provided as methods, systems, or computer program products. Therefore, the present application may adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present application may adopt 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.

[0184] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes 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 generate 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.

[0185] 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 comprising 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.

[0186] 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.

[0187] Although the embodiments of the present invention have been shown and described above, it is to be understood that the above embodiments are exemplary and are not to be construed as limitations of the present invention. A person skilled in the art may change, modify, replace and vary the above embodiments within the scope of the present invention.

Claims

1. A method for extracting understory terrain based on interferometric radar images, characterized in that: The steps include: S1: Obtain the dual-station InSAR image data of the forest terrain area to be inverted, perform interferometric processing, obtain the initial InSAR data coherence, perform coherence correction on it to obtain the coherence correction factor, and obtain the corrected coherence amplitude; S2: Perform azimuth spectrum filtering on the original InSAR image to achieve sub-aperture decomposition and obtain sub-view SAR images under different spectra; S3: respectively obtain the complex coherence coefficients of the sub-view SAR images under different spectra, and perform coherence amplitude correction, complex coherence coefficient reconstruction and point expansion on the complex plane for the complex coherence coefficients under different spectra; S4: For any pixel, the points within the complex unit circle are subjected to gross error elimination and least square fitting of the coherence line, the slope and intercept of the coherence line are calculated, and the phase information of the two intersection points of the coherence line and the complex unit circle is further obtained; S5: Calculate the effective vertical wave number through the baseline parameter file, and determine one of the two intersection points as the surface phase point according to the positive and negative values ​​of the effective vertical wave number and the phase magnitude relationship between the two intersection points obtained above; S6: The surface phase information and the coherence amplitude of the initial InSAR data coherence are reorganized into a new complex coherence coefficient, and then de-flattening, de-terraining, phase filtering, unwrapping, de-orbiting and phase height conversion are performed in sequence to finally obtain the understory terrain elevation information of the understory terrain area to be inverted.

2. The method for extracting understory terrain based on interferometric radar images according to claim 1, characterized in that: In step S1, the coherence correction factor includes the signal-to-noise ratio decoherence factor and quantization error decoherence factor Quantization error decoherence factor It represents the measure of the decoherence factor caused in the data processing process, and its value is a constant.

3. The method for extracting understory terrain based on interferometric radar images according to claim 1, characterized in that: In step S2, when the original InSAR image is subjected to azimuth spectrum filtering to realize sub-aperture decomposition, the spectrum center of the sub-view SAR image is expressed as: ， ； In the formula, represents the spectrum center of the kth sub-view SAR image; N represents the number of sub-view SAR images split along the azimuth direction; Indicates the split bandwidth.

4. The method for extracting understory terrain based on interferometric radar images according to claim 1, characterized in that: Step S3 specifically includes: Obtain the complex coherence coefficient of sub-view SAR images under different spectra , k=1,2,…,N, where N represents the number of sub-view SAR images segmented along the azimuth direction; The coherence amplitude correction is performed on the complex coherence coefficients under different spectra. At this time, the coherence correction factor obtained in step S1 is substituted into the following expression for coherence correction: ； In the formula, the coherence correction factor includes the signal-to-noise ratio decoherence factor and quantization error decoherence factor ; represents the pure body decoherence factor in the interferometric process of SAR images with different sub-views; Indicates the amplitude; The reconstruction expression of complex coherence coefficient of SAR images with different sub-views is: ； In the formula, represents the complex coherence coefficient of the reconstructed k-th sub-view SAR image; It means taking the phase of a complex number; i means the imaginary part; The complex coherence coefficient of N reconstructed SAR images with different sub-views Expand the points to the complex plane; within the complex unit circle of the complex plane, any complex coherence coefficient The phase range is , the amplitude range is .

5. The method for extracting understory terrain based on interferometric radar images according to claim 1, characterized in that: Step S4 specifically includes: For a given set of N points within the complex unit circle on the complex plane, the real and imaginary parts are , , centralize the data: ； ； Where N represents the number of sub-view SAR images segmented along the azimuth direction; , denote the mean of the real part and the mean of the imaginary part respectively; , Respectively represent the values ​​of the real and imaginary parts of the k-th point after data centering; Perform the least squares straight line fitting on N points, calculate the parameters of the fitting line, and introduce the following auxiliary quantities: ； In the formula, , and There are three auxiliary quantities; Solve to obtain two candidate values ​​for the slope of the fitted line: ； In the formula, and are two candidate slope values; The intercepts corresponding to the above two slopes are: ； In the formula, and is the intercept corresponding to the above two slopes; For outliers, the gross errors are eliminated and they are not used as candidate points for fitting the straight line. Specifically, the residual from the point to the straight line is defined as the vertical distance, and the and Respectively represent the residuals from the point to the two fitted coherent lines: Calculate the total residuals from all points to the two fitted coherence lines separately and ; Select the fitted coherence line with the smallest total residual as the final fitted coherence line; Calculate the mean residual of all points and standard deviation ; Suppose the elimination residual exceeds the threshold , the residual of the removed point is greater than All points of Refit the coherence line based on the retained points and solve its two intersection points with the complex unit circle and ; The phases corresponding to the two intersection points are extracted as follows: ； In the formula, and are the phases corresponding to the two intersection points respectively; represents the phase of a complex number; i represents the imaginary part.

6. The method for extracting understory terrain based on interferometric radar images according to claim 1, characterized in that: In step S5, one of the two intersection points is determined to be a surface phase point based on the positive and negative values ​​of the effective vertical wave number and the phase magnitude relationship between the two intersection points obtained above. The determination criteria are: ； In the formula, represents the effective vertical wave number; It means taking the phase of a complex number; and They represent the phases corresponding to the two intersection points respectively; Indicates the surface phase; express and Conjugate multiplication; i represents the imaginary part.

7. The method for extracting understory terrain based on interferometric radar images according to claim 1, characterized in that: In step S6, the surface phase information and the coherence amplitude of the initial InSAR data coherence after correction are reorganized into a new complex coherence coefficient, which is expressed as follows: ； In the formula, represents the new complex coherence coefficient; represents the corrected coherence amplitude obtained in step S1; Indicates the surface phase; represents the phase of a complex number; i represents the imaginary part.

8. A forest terrain extraction system based on interferometric radar images, characterized in that: include: The coherence amplitude correction module is used to obtain the dual-station InSAR image data of the forest terrain area to be inverted, and perform interference processing to obtain the initial InSAR data coherence, and perform coherence correction on it to obtain the coherence correction factor, and obtain the corrected coherence amplitude; The sub-aperture decomposition module is used to perform azimuth spectrum filtering on the original InSAR image to achieve sub-aperture decomposition and obtain sub-view SAR images under different spectra; The sub-view SAR image processing module is used to obtain the complex coherence coefficients of the sub-view SAR images under different spectra, and perform coherence amplitude correction, complex coherence coefficient reconstruction and point expansion on the complex plane for the complex coherence coefficients under different spectra in turn; The intersection phase acquisition module is used to remove gross errors and perform least squares fitting of the coherent line on the points within the complex unit circle for any pixel, calculate the slope and intercept of the coherent line, and further obtain the phase information of the two intersection points of the coherent line and the complex unit circle; The surface phase point determination module is used to calculate the effective vertical wave number through the baseline parameter file, and determine one of the two intersection points as the surface phase point according to the positive and negative values ​​of the effective vertical wave number and the phase magnitude relationship between the two intersection points obtained above; The understory terrain elevation extraction module is used to reorganize the surface phase information and the coherence amplitude of the initial InSAR data coherence after correction into a new complex coherence coefficient, and then perform de-flattening, de-terraining, phase filtering, unwrapping, de-orbiting and phase height conversion in sequence, and finally obtain the understory terrain elevation information of the understory terrain area to be inverted.

9. An electronic device, characterized in that: include: Memory on which computer programs or instructions are stored; A processor is used to load and execute the computer program or instruction to implement the method for extracting understory terrain based on interferometric radar images as described in any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program or instruction stored thereon, characterized in that: When the computer program or instruction is executed by a processor, the method for extracting understory terrain based on interferometric radar images as described in any one of claims 1 to 7 is implemented.

Images