Large-range under-forest topographic mapping method and device based on low-frequency InSAR

Abstract

The invention discloses a low-frequency InSAR (Interferometric Synthetic Aperture Radar)-based large-range under-forest terrain surveying and mapping method and device. The method comprises the following steps: processing a low-frequency InSAR image pair of a large-range under-forest terrain inversion region to obtain a plurality of sub-view interferograms; modeling is carried out on each sub-view interferogram of each image, an under-forest terrain phase is determined according to a complex coherence coefficient of pixels at the same spatial position, and then phase height conversion is carried out to obtain an initial under-forest terrain; a block adjustment model considering orbit system errors and residual forest height errors is constructed, satellite-borne laser radar ground elevation points are used as control points, constraint is carried out on connection points of multi-scene images, and orbit system error parameters and residual forest height errors are obtained through calculation by adopting a least square configuration method. Correcting the initial under-forest topography to obtain the final under-forest topography of each single scene; and finally fusing the multi-scene under-forest terrains in the forest coverage area to obtain a large-range under-forest terrain surveying and mapping product. According to the invention, effective and steady surveying and mapping can be carried out on the large-range under-forest terrain.

Description

Technical Field

[0001] The invention belongs to the field of geodesy, and in particular relates to a large-scale forest understory terrain mapping and a device based on low-frequency InSAR. Background Art

[0002] Interferometric Synthetic Aperture Radar (InSAR) has the advantages of all-day, all-weather and large-scale, and has been widely used in large-scale terrain mapping. However, the existing large-scale or even global DEM products obtained based on high-frequency InSAR technology contain serious tree height signals and cannot be directly used to characterize accurate understory terrain.

[0003] Low-frequency dual-station InSAR technology is not affected by temporal decoherence and atmospheric delay, and can obtain high-quality interferometric data. In addition, low-frequency SAR signals have strong penetration in forest areas, and the elevation obtained is closer to the real forest surface than the existing X / C band InSAR data. However, low-frequency SAR signals will undergo volume scattering during the process of penetrating the canopy, and the measured elevation is between the top of the canopy and the forest surface. The interference of forest scattering on interferometric height measurement must be accurately eliminated to achieve large-scale, high-precision forest terrain mapping. In addition, when mapping large-scale forest terrain, regional block adjustment is required to remove the systematic errors contained in InSAR terrain products. The existing regional block adjustment algorithm is modeled by combining the satellite-borne lidar control points in the bare land area and the connection points between different frames. In large-scale forest areas, there are problems such as very few control points in the bare land area and unreliable connection points that do not consider the difference in penetration depth. Therefore, it is necessary to design a regional block adjustment algorithm that can effectively eliminate the interference of forest signals on interferometric height measurement and adapt to forest areas, so as to achieve large-scale, high-precision forest terrain mapping. Summary of the invention

[0004] The present invention provides a large-scale forest understory terrain mapping and a device based on low-frequency InSAR, which can realize large-scale and high-precision forest understory terrain mapping under the condition of using only low-frequency InSAR data and a small number of control points.

[0005] In order to achieve the above technical objectives, the present invention adopts the following technical solutions: A method for mapping large-scale understory terrain based on low-frequency InSAR, comprising: Step 1, obtaining a low-frequency InSAR master-slave image pair of a large-scale forest terrain inversion area, and performing sub-aperture decomposition on the master-slave image pair to obtain a number of sub-view image pairs; Step 2, preprocessing each sub-view image pair to obtain the corresponding sub-view interference map; Step 3, taking the pixel points representing the same geographical location in all sub-view interferograms of the same InSAR image pair, and determining the understory terrain phase of each same geographical location according to the complex coherence coefficient of these pixel points; Step 4, unwrapping, absolute phase conversion, and phase height conversion are performed on each forest terrain phase to obtain the initial forest terrain; Step 5, construct a regional block adjustment model considering the InSAR orbit system error and residual forest height error; use the satellite-borne laser radar ground elevation points as control points, and select evenly distributed connection points in the image overlap area for constraints; construct a covariance function and solve the adjustment model using the least squares collocation method to obtain the orbit system error parameters and the residual forest height error; Step 6: Correct the initial understory terrain according to the track system error parameters and the residual forest height error to obtain the final understory terrain of each single scene; Step 7: fuse multiple views of forest terrain in the forest coverage area to obtain a large-scale forest terrain mapping product.

[0006] Furthermore, the master-slave image pair is decomposed into sub-apertures, specifically: Firstly, one-dimensional fast Fourier transform is performed on the master image and the slave image respectively to convert them into the azimuth spectrum domain; Then, the imaging azimuth angle range corresponding to the required sub-aperture is calculated according to the required number of sub-apertures; Finally, the frequency spectrum corresponding to the azimuth angle range of each sub-aperture is subjected to inverse fast Fourier transform to obtain the sub-view main image and sub-view slave image corresponding to each sub-aperture.

[0007] Furthermore, the preprocessing of each sub-view image pair includes: interference, deflating phase, and phase adaptive filtering.

[0008] Furthermore, determining the understory terrain phase at the same geographical location according to the complex coherence coefficients of the pixel points includes: (1) Take all pixels representing the same geographic location p in the sub-view interference graph and extract the complex coherence coefficients of these pixels; (2) Representing multiple complex coherence coefficients of the geographic location p as multiple points within the unit circle of the complex plane; (3) Fit multiple points within the unit circle of the complex plane into straight lines: ; In the formula, represents the real part of the complex coherence coefficient, represents the imaginary part of the complex coherence coefficient, and is the coefficient of the fitted line to be determined; (4) Calculate the phases corresponding to the two intersection points where the straight line intersects the unit circle of the complex plane and ; (5) If satisfied , then the forest terrain phase of the geographical location p is determined as If satisfied , then the forest terrain phase of the geographical location p is determined as .

[0009] Furthermore, the calculation formula for the phase height conversion is: ; In the formula, represents the elevation of the forest terrain at the target geographic location, represents the radar wavelength, Indicates the vertical baseline length, Indicates the distance from the center of the satellite antenna to the target geographic location. represents the unwrapped absolute phase obtained by unwrapping the understory terrain phase and converting it to the absolute phase. Represents the angle of incidence.

[0010] Furthermore, the specific process of step 5 includes: Step 5.1, construct a regional block adjustment model considering the track system error and residual forest height error, expressed as: ; In the formula, is the observation vector of the elevation difference between the initial understory terrain at the same geographical location and the satellite-borne lidar control point; is the design matrix composed of range and azimuth coordinates, and X is the coefficient matrix to be estimated; It is a Boolean matrix used to control whether the residual forest height error needs to be estimated, that is, the forest area is 1 and the non-forest area is 0; is the residual forest height error vector; is the observation error vector of the satellite-borne laser radar control point; Step 5.2: For each initial forest terrain acquired by InSAR, select the uniformly distributed Satellite-borne lidar elevation points are used as control points to provide absolute elevation constraints for the regional block adjustment model: ; In the formula, Represents control points 's elevation; Represents control points The elevation value of the corresponding geographical location in the initial understory terrain; and Represents control points The pixel coordinates of the corresponding geographical location in the range and azimuth directions in the SAR image coordinate system, represents the sum of the corresponding orbit error and the residual forest height error; In the case of control points, the matrix expression of the adjustment model is: ; In the formula, is the orbit system error parameter vector The 6 coefficients in for The residual forest height error of each control point is is the observation error corresponding to m control points; Step 5.3: For any initial understory terrain, traverse other initial understory terrain results in the geographic coordinate system, find the data with overlapping areas, and then evenly select several connection points in the overlapping area; based on the fact that the understory terrain obtained by different InSAR data at the same geographical location is equal, establish a relative elevation constraint relationship for the regional network adjustment model; for any connection point The following functional relationship is established: ; Where J and K represent the initial understory terrain results of two scenes with overlapping areas; and Respectively represent the connection points The initial understory terrain elevation values ​​corresponding to the geographical locations in J and K, and Represent the connection points after adjustment The sum of the orbital errors and residual forest height errors corresponding to the geographical locations in J and K; Step 5.4, construct the covariance function; First, a regional block adjustment is performed using all satellite-borne lidar control points to fit , to estimate the empirical variance and covariance: ; In the formula, is the empirical variance, The mutual distance is The empirical covariance of the control point pairs; is the total number of control points, and Represents the mutual distance different control points in the control point pair; The mutual distance is The number of control point pairs, and For different control points , The fitted residuals of Then, by choosing the empirical covariance function , and then according to the calculated empirical variance and multiple empirical covariances To solve the unknown parameters in the empirical covariance function; Step 5.5, based on the empirical covariance function obtained , using the least squares collocation adjustment method to solve the orbit system error parameter vector , residual forest height signals for pixels covered by control points and pixels not covered by control points and , respectively: ; In the formula, , is the covariance matrix of the observation error, represents the covariance matrix between the control point observation pixels, Represents the covariance matrix of the control point observation pixels and non-observation pixels, which is given by the empirical covariance function Calculated; , , are the orbit system error parameter vectors , the residual forest height error vector of the control point pixel and the residual forest height error vector for uncovered pixels The solution value of Step 5.6, based on the calculated orbital system error parameter vector, the residual forest height signal of the control point pixel, and the residual forest height signal of the uncovered pixel, the initial understory terrain is adjusted to obtain the final understory terrain corresponding to the current image, which is expressed as: ; In the formula, Indicates The final understory topography of the scene, Indicates The initial understory topography of the scene, Indicates The orbit system error parameter vector solution value of the scene, For the The residual forest height error vector for pixels not covered by the background.

[0011] Furthermore, the following Gauss-Markov model is selected to construct the empirical covariance function: ; In the formula, represents the distance between two pixels, For the distance The empirical covariance function of and These are all unknown parameters to be solved in the empirical covariance function.

[0012] Furthermore, the specific process of step 7 includes: First, extract all the understory terrain elevation values ​​and coherent amplitude values ​​at the same geographic coordinate location, and calculate the weight factor corresponding to each elevation: ; In the formula, Indicates Elevation value The corresponding weights, represents the full-aperture coherence amplitude, is the incident angle, represents the radar wavelength; is the vertical baseline length; Represents the standard deviation of relative elevation error; Then, all elevation points are weighted averaged to obtain the final understory terrain elevation. : ; In the formula, Represents the number of all understory terrain elevation values ​​at the same geographic coordinate location.

[0013] A large-scale forest understory terrain mapping device based on low-frequency InSAR, comprising: The data preprocessing module is used to perform registration, spectrum filtering, sub-aperture decomposition, interference and flat ground phase removal on the low-frequency InSAR master-slave image pairs acquired in the large-scale understory terrain inversion area to obtain several sub-view image pairs; The initial understory terrain estimation module is used to: take the pixel points representing the same geographical location in all sub-view interferograms of the same InSAR image pair, and determine the understory terrain phase of each same geographical location according to the complex coherence coefficient of these pixel points; then perform unwrapping, absolute phase conversion, and phase height conversion on each understory terrain phase to obtain the initial understory terrain; The regional block adjustment module is used to: construct a regional block adjustment model that takes into account the InSAR orbit system error and residual forest height error; use the satellite-borne lidar ground elevation points as control points, and select evenly distributed connection points in the image overlap area for constraints; construct a covariance function and solve the adjustment model using the least squares collocation method to obtain the orbit system error parameters and residual forest height error, and then correct the initial understory terrain to obtain the final understory terrain for each single scene; The fusion mosaic module is used to fuse multiple views of understory terrain in a large forest area to obtain large-scale understory terrain mapping products.

[0014] Beneficial Effects

[0015] The invention discloses a large-scale forest understory terrain mapping method based on low-frequency InSAR technology, which realizes obtaining forest understory terrain by using only InSAR data and a small number of control points: firstly, low-frequency InSAR data of a large-scale forest understory terrain inversion area is obtained, and sub-aperture decomposition is performed on the main image and the slave image of each registered InSAR interference pair respectively; secondly, interference, flat ground phase removal, filtering and other operations are performed on the sub-interference pairs after sub-aperture decomposition respectively to obtain multiple sub-view interference graphs; then, a straight line is used to fit the phase of the sub-view interference graph to obtain the forest understory terrain phase; the forest understory terrain phase is converted to elevation to obtain the initial forest understory terrain; then, considering the change trend of InSAR orbit error and the statistical law of residual forest height error, a regional block adjustment model based on least squares registration is constructed, and connection points and satellite-borne laser radar ground elevation points are selected for adjustment calculation to solve the unknown parameters of the model; finally, the orbit error and the residual forest height error are calculated according to the error model and removed from the initial forest understory terrain to obtain the large-scale forest understory terrain.

[0016] The beneficial effects of this method are as follows: a phase acquisition technology for forest terrain based on sub-aperture decomposition is constructed, which overcomes the limitation of existing technologies that rely on fully polarized low-frequency InSAR data; a regional network adjustment method for forest terrain that takes into account system errors such as orbits and residual forest height errors is proposed, which solves the problems faced by existing methods such as very few control points in bare land areas and unreliable connection points that do not consider penetration depth differences; a robust forest terrain result fusion strategy under different framing and multiple coverage conditions is established, which realizes the automatic generation of large-scale forest terrain products. In terms of accuracy, the forest terrain products obtained by the method of the present invention are superior to existing public terrain products, and are an effective method for autonomous, controllable, large-scale, and high-precision forest terrain mapping. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 The present invention is a flowchart of the method described in the embodiment of the present invention.

[0018] Figure 2 The location of the verification area and the coverage of the low-frequency InSAR image interference pair selected by the present invention.

[0019] Figure 3 Plots of the variance, covariance distributions and the fitted empirical covariance function.

[0020] Figure 4 These are the difference maps between the initial understory topography, the understory topography before and after block adjustment and the reference DEM. Figure 4(a) represents the initial understory terrain in the geographic coordinate system, Figure 4 (b) shows the difference between the initial understory topography and the reference DEM before block adjustment in the geographic coordinate system. Figure 4 (c) The difference map between the final understory topography and the reference DEM after regional block adjustment in the geographic coordinate system.

[0021] Figure 5 This is a difference map between the traditional InSAR DEM product and the understory terrain obtained by the present invention and the airborne DTM. Figure 5 (a) shows the difference between the DEM product obtained by the existing traditional InSAR technology and the airborne DTM in the geographic coordinate system. Figure 5 (b) is the difference map between the understory terrain obtained by the present invention and the airborne DTM. DETAILED DESCRIPTION

[0022] The following is a detailed description of an embodiment of the present invention. This embodiment is based on the technical solution of the present invention, and provides a detailed implementation method and a specific operation process to further explain the technical solution of the present invention.

[0023] In order to better illustrate the method and steps of the present invention, the present invention is further described in detail using LT-1 dual-station L-band InSAR data and space-borne ICESat-2 data in a test area located in the southern part of Guangdong Province, China; it is noted that the specific implementation described here is only used to explain the present invention and is not intended to limit the present invention.

[0024] This embodiment provides a large-scale forest terrain mapping method based on low-frequency InSAR technology. Figure 1 As shown, the following steps are included:

[0025] Step 1: obtain low-frequency InSAR master-slave image pairs of a large-scale forest terrain inversion area, and perform sub-aperture decomposition on the registered InSAR master image and slave image to obtain several sub-view image pairs.

[0026] This implementation example selects the LT-1 dual-station InSAR data located in the southern experimental area of ​​Guangdong Province, China. The data of 4 adjacent tracks are selected. Each track contains 5 InSAR data, totaling 20 interferometric pairs. The experimental area is located as follows: Figure 2As shown; the acquired data is registered and spectrally filtered, which is a prior art method, and reference can be made to the literature: Jin Guowang, Xu Qing, Zhang Hongmin. Synthetic Aperture Radar Interferometry [M]. National Defense Industry Press, 2014, which will not be elaborated in detail in this embodiment. Afterwards, each pair of registered master and slave SAR images is respectively subjected to one-dimensional fast Fourier transform to convert them into the azimuth spectrum domain; then, the imaging azimuth angle range corresponding to the required sub-aperture is calculated according to the required number of sub-apertures (greater than 2); finally, the spectrum corresponding to the azimuth angle range of all sub-apertures is subjected to inverse fast Fourier transform to obtain the sub-view interference pair.

[0027] Step 2: Perform interference, flat ground phase removal, phase adaptive filtering and other operations on the sub-view interference pairs to obtain several sub-view interference graphs.

[0028] Interferometry is performed on each sub-aperture image pair, and then the flat-earth phase is calculated using the orbital parameters of the SAR satellite: ; In the formula, represents the radar wavelength, Indicates the length of the parallel baseline; in order to suppress the influence of noise on the interference phase, it is necessary to filter the de-flattened phase. The filtering method usually adopted is the adaptive filtering method, which is an existing mature technology and will not be elaborated in detail in this embodiment. Since the dual-station InSAR data is not affected by time decorrelation and atmospheric delay, the de-flattened and filtered phase can be considered to contain only the terrain phase, forest signal phase and orbit and other system error phases.

[0029] Step 3, expand the complex coherence coefficients of the same position of the sub-view interference pattern within the unit circle of the complex plane, and use the linear equation to fit the phases of different sub-view interference patterns, and select the phase at the intersection with the unit circle as the understory terrain phase;

[0030] The sub-interference graph consists of two parts: interference amplitude and phase, which are usually expressed in complex format. For the sub-view interference graphs of the same interference pair, it is necessary to extract the complex coherence coefficient of each sub-view interference graph at the same position pixel by pixel and expand it into multiple points in the unit circle of the complex plane. Then, these points are fitted with a straight line: ; In the formula, represents the real part of the complex coherence coefficient, represents the imaginary part of the complex coherence coefficient, and is the coefficient to be determined. Then, calculate the phases corresponding to the two intersection points where the straight line intersects the unit circle of the complex plane and , the understory terrain phase is determined based on the following criteria: ; Step 4, performing phase unwrapping, absolute phase conversion, phase height conversion, etc. on the understory terrain to obtain the initial understory terrain; The phase of the entangled forest terrain obtained in step 3 is between [-π, π]. This embodiment uses the minimum cost flow method to untangle the phase to obtain the untangled relative phase; because the untangled phase is the relative phase relative to the untangled reference point, it cannot be directly used for phase height conversion. It is necessary to select control points and use robust least squares to convert it to absolute phase. The above operation is an existing mature technical method and will not be elaborated in detail in this embodiment. Finally, the absolute phase is converted to elevation using the following formula: ; In the formula, represents the radar wavelength, Indicates the vertical baseline length, Indicates the distance from the center of the satellite antenna to the target point. represents the unwrapped absolute phase obtained by unwrapping the understory terrain phase and converting it to the absolute phase. represents the angle of incidence. The initial understory terrain is as follows Figure 4 (a) shown.

[0031] Step 5: Based on the changing trend of the InSAR orbit system error and the statistical law of the residual forest height error, a regional block adjustment model considering the residual forest height error is constructed; the satellite-borne lidar ground elevation points are used as control points, and evenly distributed connection points are selected in the overlapping area for constraints; a covariance function is constructed and the adjustment model is solved by the least squares configuration method; finally, the orbit system error and the residual forest height error are removed from the initial understory terrain to obtain the understory terrain of a single scene; Step 5.1, construct a regional block adjustment model considering the track system error and residual forest height error, expressed as: ; In the formula, is the observation vector of the elevation difference between the initial understory terrain at the same geographical location and the satellite-borne lidar control point; is the design matrix composed of range and azimuth coordinates, X is the coefficient matrix to be estimated, is a Boolean matrix that controls whether the residual forest height error needs to be estimated (i.e., 1 for forest areas and 0 for non-forest areas). is the residual forest height error vector, is the observation error vector of the satellite-borne lidar control point.

[0032] Step 5.2: For each initial forest terrain acquired by InSAR, select the terrain with a relatively uniform distribution within the image coverage. Satellite-borne lidar elevation points are used as control points to provide absolute elevation constraints for the model; since the residual forest height error needs to be modeled, control points located in the forest area need to be added. The following functional relationship can be established for each control point: ; In the formula, Represents control points 's elevation; Represents control points The elevation value of the corresponding geographical location in the initial understory terrain; and Represents control points The pixel coordinates of the corresponding geographical location in the range and azimuth directions in the SAR image coordinate system, represents the sum of the corresponding orbit error and the residual forest height error; In the case of control points, the matrix expression corresponding to formula (5) is: ; In the formula, is the orbit system error parameter vector The 6 coefficients in for The residual forest height error of each control point is for The residual corresponding to the control points; Step 5.3: For any initial understory terrain, traverse other initial understory terrain results in the geographic coordinate system, find the data with overlapping areas, and then evenly select locations with relatively flat terrain and high coherence in the overlapping area as connection points; based on the fact that the understory terrain obtained by different InSAR data at the same geographical location is equal, establish a relative elevation constraint relationship for the model; for any connection point The following functional relationship is established: ; Where J and K represent the initial understory terrain results of two scenes with overlapping areas; and They represent the initial understory terrain elevation values ​​corresponding to the geographical location of the connection point j in J and K, respectively. and They represent the sum of the orbit error and residual forest height error of the geographical location corresponding to the connection point j in J and K after adjustment, respectively.

[0033] Step 5.4, construct the covariance function. First, all satellite-borne lidar control points are used to perform a regional network adjustment based on the least squares theory to fit , to estimate the empirical variance and covariance: ; In the formula, is the empirical variance, is the empirical covariance, is the total number of control points, is the number of point pairs within a given distance interval; is the distance between two control points, expressed as Euclidean distance and rounded; and is the fitting residual of different control points, that is, the residual matrix expression of m control points is ,in The unknown parameters solved for the first adjustment fitting. Then, calculate the different distance intervals The empirical covariance under , and according to its distribution characteristics, the following Gauss-Markov model is selected as the empirical covariance function: ; Where d represents the distance between two pixels and the correlation length ,variance As unknown parameters, the variance and covariance calculated by formula (9) need to be estimated by the least squares method. The estimated variance and covariance distribution is as follows: Figure 3 As shown in the middle circle, the fitted empirical covariance function is as follows Figure 3 Shown by the middle curve.

[0034] Step 5.5, using the empirical covariance function established by formula (10) and the least squares collocation adjustment method, the orbit system error parameters , residual forest height signals for control point pixels and uncovered pixels and for: ; In the formula, , is the covariance matrix of the observed errors (residuals), represents the covariance matrix between the control point observation pixels, which is obtained by the empirical covariance calculated by formula (9); The covariance matrix representing the observed pixels and non-observed pixels of the control point is calculated by the empirical covariance function, i.e., formula (10). The model solution process of the least squares collocation method in the above formula can be referred to: Cui Xizhang, Yu Zongchou, Tao Benzao, etc. Generalized Survey Adjustment (New Edition) [M]. Wuhan University of Surveying and Mapping Press, 2001. This embodiment will not be elaborated in detail.

[0035] Step 5.6, based on the calculated orbital system error parameter vector, the residual forest height signal of the control point pixel, and the residual forest height signal of the uncovered pixel, the initial understory terrain is adjusted to obtain the final understory terrain corresponding to the current image, which is expressed as: ; In the formula, Indicates The final understory topography of the scene, Indicates The initial understory topography of the scene, Indicates The orbit system error parameter vector solution value of the scene, For the The residual forest height error vector for pixels not covered by the background.

[0036] Step 6: fuse and mosaic the single-scene forest terrain obtained in step 5 to obtain a large-scale forest terrain product.

[0037] First, extract all the understory terrain elevation values ​​and coherent amplitude values ​​at the same geographic coordinate location, and calculate the weight factor corresponding to each elevation: ; In the formula, Indicates The weight corresponding to the elevation value, represents the full-aperture coherence amplitude, is the incident angle, is the vertical baseline length; Represents the standard deviation of relative elevation error.

[0038] Then, all elevation points are weighted averaged to obtain the final understory terrain elevation: ; In the formula, Represents the number of all understory terrain elevation values ​​at the same geographic coordinate location.

[0039] Figure 4 (a) represents the initial understory terrain in the geographic coordinate system. It is found that there are obvious jumps between the initial understory terrain of different orbits. This is because the data of different orbits are acquired at different times, and the accuracy of satellite orbit determination and baseline estimation is inconsistent. It needs to be processed using regional block adjustment. Figure 4 (b) is the difference map between the initial understory topography and the reference DEM before regional network adjustment. Figure 4(c) is the difference map between the final forest terrain after block adjustment and the reference DEM. It is found that the systematic error is removed after block adjustment, and there is no obvious jump between different tracks. In order to more clearly highlight the advantages of using low-frequency InSAR to perform large-scale forest terrain mapping in this embodiment, airborne high-precision LiDAR DTM is used as a reference and compared with the existing public high-frequency InSAR DEM product (TanDEM-X DEM). Figure 5 (a) is the difference between the traditional InSAR DEM and the airborne DTM. Figure 5 (b) is the difference map between the forest understory terrain of this embodiment and the airborne DTM. Figure 5 (b) The difference of LiDAR DTM is smaller, RMSE decreases from 5.87m to 3.09m, and MAE decreases from 4.59m to 0.34m, indicating that the understory terrain results mapped by low-frequency InSAR in this embodiment are more accurate.

[0040] The above embodiments are preferred embodiments of the present application. Ordinary technicians in this field can also make various changes or improvements on this basis. Without departing from the overall concept of the present application, these changes or improvements should fall within the scope of protection required by the present application.

Claims

1. A method for mapping large-scale forest terrain based on low-frequency InSAR, characterized in that: include: Step 1, obtaining a low-frequency InSAR master-slave image pair of a large-scale forest terrain inversion area, and performing sub-aperture decomposition on the master-slave image pair to obtain a number of sub-view image pairs; Step 2, preprocessing each sub-view image pair to obtain the corresponding sub-view interference map; Step 3, taking the pixel points representing the same geographical location in all sub-view interferograms of the same InSAR image pair, and determining the understory terrain phase of each same geographical location according to the complex coherence coefficient of these pixel points; Step 4, unwrapping, absolute phase conversion, and phase height conversion are performed on each forest terrain phase to obtain the initial forest terrain; Step 5, construct a regional block adjustment model considering the InSAR orbit system error and residual forest height error; use the satellite-borne laser radar ground elevation points as control points, and select evenly distributed connection points in the image overlap area for constraints; construct a covariance function and solve the adjustment model using the least squares collocation method to obtain the orbit system error parameters and the residual forest height error; Step 6: Correct the initial understory terrain according to the track system error parameters and the residual forest height error to obtain the final understory terrain of each single scene; Step 7: fuse multiple views of forest terrain in the forest coverage area to obtain a large-scale forest terrain mapping product.

2. The method for mapping large-scale forest terrain based on low-frequency InSAR according to claim 1, characterized in that: The master-slave image pair is decomposed into sub-apertures as follows: Firstly, one-dimensional fast Fourier transform is performed on the master image and the slave image respectively to convert them into the azimuth spectrum domain; Then, the imaging azimuth angle range corresponding to the required sub-aperture is calculated according to the required number of sub-apertures; Finally, the frequency spectrum corresponding to the azimuth angle range of each sub-aperture is subjected to inverse fast Fourier transform to obtain the sub-view main image and sub-view slave image corresponding to each sub-aperture.

3. The method for mapping large-scale forest terrain based on low-frequency InSAR according to claim 1, characterized in that: The preprocessing of each sub-view image pair includes: interference, removing the flat phase, and phase adaptive filtering.

4. The method for mapping large-scale forest terrain based on low-frequency InSAR according to claim 1, characterized in that: Determining the understory terrain phase at the same geographical location according to the complex coherence coefficients of the pixel points includes: (1) Take all pixels representing the same geographic location p in the sub-view interference graph and extract the complex coherence coefficients of these pixels; (2) Representing multiple complex coherence coefficients of the geographic location p as multiple points within the unit circle of the complex plane; (3) Fit multiple points within the unit circle of the complex plane into straight lines: ； In the formula, represents the real part of the complex coherence coefficient, represents the imaginary part of the complex coherence coefficient, and is the coefficient of the fitted line to be determined; (4) Calculate the phases corresponding to the two intersection points where the straight line intersects the unit circle of the complex plane and ; (5) If satisfied , then the forest terrain phase of the geographical location p is determined as If satisfied , then the forest terrain phase of the geographical location p is determined as .

5. The method for mapping large-scale forest terrain based on low-frequency InSAR according to claim 1, characterized in that: The calculation formula for the phase height conversion is: ； In the formula, represents the elevation of the forest terrain at the target geographic location, represents the radar wavelength, Indicates the vertical baseline length, Indicates the distance from the center of the satellite antenna to the target geographic location. represents the unwrapped absolute phase obtained by unwrapping the understory terrain phase and converting it to the absolute phase. Represents the angle of incidence.

6. The method for mapping large-scale forest terrain based on low-frequency InSAR according to claim 1, characterized in that: The specific process of step 5 includes: Step 5.1, construct a regional block adjustment model considering the track system error and residual forest height error, expressed as: ； In the formula, is the observation vector of the elevation difference between the initial understory terrain at the same geographical location and the satellite-borne lidar control point; is the design matrix composed of range and azimuth coordinates, and X is the coefficient matrix to be estimated; It is a Boolean matrix used to control whether the residual forest height error needs to be estimated, that is, the forest area is 1 and the non-forest area is 0; is the residual forest height error vector; is the observation error vector of the satellite-borne laser radar control point; Step 5.2: For each initial forest terrain acquired by InSAR, select the uniformly distributed Satellite-borne lidar elevation points are used as control points to provide absolute elevation constraints for the regional block adjustment model: ； In the formula, Represents control points 's elevation; Represents control points The elevation value of the corresponding geographical location in the initial understory terrain; and Represents control points The pixel coordinates of the corresponding geographical location in the range and azimuth directions in the SAR image coordinate system, represents the sum of the corresponding orbit error and the residual forest height error; In the case of control points, the matrix expression of the adjustment model is: ； In the formula, is the orbit system error parameter vector The 6 coefficients in for The residual forest height error of each control point is is the observation error corresponding to m control points; Step 5.3: For any initial understory terrain, traverse other initial understory terrain results in the geographic coordinate system, find the data with overlapping areas, and then evenly select several connection points in the overlapping area; based on the fact that the understory terrain obtained by different InSAR data at the same geographical location is equal, establish a relative elevation constraint relationship for the regional network adjustment model; for any connection point The following functional relationship is established: ； Where J and K represent the initial understory terrain results of two scenes with overlapping areas; and Respectively represent the connection points The initial understory terrain elevation values ​​corresponding to the geographical locations in J and K, and Represent the connection points after adjustment The sum of the orbital errors and residual forest height errors corresponding to the geographical locations in J and K; Step 5.4, construct the covariance function; First, a regional block adjustment is performed using all satellite-borne lidar control points to fit , to estimate the empirical variance and covariance: ； In the formula, is the empirical variance, The mutual distance is The empirical covariance of the control point pairs; is the total number of control points, and Represents the mutual distance different control points in the control point pair; The mutual distance is The number of control point pairs, and For different control points , The fitted residuals of Then, by choosing the empirical covariance function , and then according to the calculated empirical variance and multiple empirical covariances To solve the unknown parameters in the empirical covariance function; Step 5.5, based on the empirical covariance function obtained , using the least squares collocation adjustment method to solve the orbit system error parameter vector , residual forest height signals for pixels covered by control points and pixels not covered by control points and , respectively: ； In the formula, , is the covariance matrix of the observation error, represents the covariance matrix between the control point observation pixels, Represents the covariance matrix of the control point observation pixels and non-observation pixels, which is given by the empirical covariance function Calculated; , , are the orbit system error parameter vectors , the residual forest height error vector of the control point pixel and the residual forest height error vector for uncovered pixels The solution value of Step 5.6, based on the calculated orbital system error parameter vector, the residual forest height signal of the control point pixel, and the residual forest height signal of the uncovered pixel, the initial understory terrain is adjusted to obtain the final understory terrain corresponding to the current image, which is expressed as: ； In the formula, Indicates The final understory topography of the scene, Indicates The initial understory topography of the scene, Indicates The orbit system error parameter vector solution value of the scene, For the The residual forest height error vector for pixels not covered by the background.

7. The method for mapping large-scale forest terrain based on low-frequency InSAR according to claim 6, characterized in that: The following Gauss-Markov model is selected to construct the empirical covariance function: ； In the formula, represents the distance between two pixels, For the distance The empirical covariance function of and These are all unknown parameters to be solved in the empirical covariance function.

8. The method for mapping large-scale forest terrain based on low-frequency InSAR according to claim 1, characterized in that: The specific process of step 7 includes: First, extract all the understory terrain elevation values ​​and coherent amplitude values ​​at the same geographic coordinate location, and calculate the weight factor corresponding to each elevation: ； In the formula, Indicates Elevation value The corresponding weight, represents the full-aperture coherence amplitude, is the incident angle, represents the radar wavelength; is the vertical baseline length; Represents the standard deviation of relative elevation error; Then, all elevation points are weighted averaged to obtain the final understory terrain elevation. : ； In the formula, Represents the number of all understory terrain elevation values ​​at the same geographic coordinate location.

9. A large-scale forest terrain mapping device based on low-frequency InSAR, characterized in that: include: The data preprocessing module is used to perform registration, spectrum filtering, sub-aperture decomposition, interference and flat ground phase removal on the low-frequency InSAR master-slave image pairs acquired in the large-scale understory terrain inversion area to obtain several sub-view image pairs; The initial understory terrain estimation module is used to: take the pixel points representing the same geographical location in all sub-view interferograms of the same InSAR image pair, and determine the understory terrain phase of each same geographical location according to the complex coherence coefficient of these pixel points; then perform unwrapping, absolute phase conversion, and phase height conversion on each understory terrain phase to obtain the initial understory terrain; The regional block adjustment module is used to: construct a regional block adjustment model that takes into account the InSAR orbit system error and residual forest height error; use the satellite-borne lidar ground elevation points as control points, and select evenly distributed connection points in the image overlap area for constraints; construct a covariance function and solve the adjustment model using the least squares collocation method to obtain the orbit system error parameters and residual forest height error, and then correct the initial understory terrain to obtain the final understory terrain for each single scene; The fusion mosaic module is used to fuse multiple views of understory terrain in a large forest area to obtain large-scale understory terrain mapping products.

Images