Method and device for large-range underforest terrain mapping based on low-frequency InSAR

Abstract

The present invention discloses a method and device for large-scale underforest terrain mapping based on low-frequency InSAR, including: processing low-frequency InSAR image pairs of a large-scale underforest terrain inversion area to obtain a number of sub-view interferograms; modeling each sub-view interferogram of each scene image, determining the underforest terrain phase according to the complex coherence coefficient of pixels at the same spatial position, and then performing phase-height conversion to obtain the initial underforest terrain; constructing a regional network adjustment model considering orbital system errors and residual forest height errors, using spaceborne lidar ground elevation points as control points, and connecting points of multiple scene images for constraint, and using the least squares collocation method to solve for orbital system error parameters and residual forest height errors, so as to correct the initial underforest terrain to obtain the final underforest terrain of each single scene; finally, fusing multiple scene underforest terrains within the forest coverage area to obtain a large-scale underforest terrain mapping product. The present invention can effectively and robustly map large-scale underforest terrains.

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:

[0006] A method for mapping large-scale understory terrain based on low-frequency InSAR, comprising:

[0007] 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;

[0008] Step 2: Preprocess each sub-view image pair to obtain the corresponding sub-view interferogram;

[0009] Step 3: Select the pixel points representing the same geographical location in all the sub-view interferograms of the same InSAR image pair, and determine the under-forest terrain phase at each same geographical location according to the complex coherence coefficients of these pixel points;

[0010] Step 4: Unwrap, perform absolute phase conversion, and phase-height conversion on each under-forest terrain phase to obtain the initial under-forest terrain;

[0011] Step 5: Construct a block adjustment model considering the InSAR orbit system error and the residual forest height error; use the ground elevation points of spaceborne lidar as control points, and select evenly distributed tie points in the image overlap area for constraint; construct a covariance function and solve the adjustment model by the method of least squares collocation to obtain the orbit system error parameters and the residual forest height error;

[0012] Step 6: Correct the initial under-forest terrain according to the orbit system error parameters and the residual forest height error to obtain the final under-forest terrain of each single scene;

[0013] Step 7: Fuse the under-forest terrains of multiple scenes within the forest-covered area to obtain a large-scale under-forest terrain mapping product.

[0014] Furthermore, perform sub-aperture decomposition on the master-slave image pair, specifically:

[0015] First, perform one-dimensional fast Fourier transform on the master image and the slave image respectively to transform them into the azimuth spectral domain;

[0016] Then, calculate the imaging azimuth angle range corresponding to the required sub-apertures according to the number of required sub-apertures;

[0017] Finally, perform inverse fast Fourier transform on the spectra corresponding to the azimuth angle ranges of each sub-aperture to obtain the sub-view master image and the sub-view slave image corresponding to each sub-aperture.

[0018] Furthermore, the preprocessing of each sub-view image pair includes: interference, removal of flat-earth phase, and phase adaptive filtering.

[0019] Furthermore, the determination of the under-forest terrain phase at each same geographical location according to the complex coherence coefficients of these pixel points includes:

[0020] (1) Select the pixel points representing the same geographical location p in all the sub-view interferograms, and extract the complex coherence coefficients of these pixel points;

[0021] (2) Represent the multiple complex coherence coefficients at the geographical location p as multiple points within the unit circle of the complex plane;

[0022] (3) Perform a linear fit on multiple points within the unit circle of the complex plane:

[0023] ;

[0024] Wherein, represents the real part of the complex coherence coefficient, represents the imaginary part of the complex coherence coefficient, and are the coefficients to be determined for the fitting line;

[0025] (4) Calculate the phases corresponding to the two intersection points where the line intersects the unit circle of the complex plane and ;

[0026] (5) If is satisfied, then determine that the understory terrain phase at geographical location p is ; If is satisfied, then determine that the understory terrain phase at geographical location p is .

[0027] Further, the calculation formula for the phase-height conversion is:

[0028] ;

[0029] Wherein, represents the understory terrain elevation at the target geographical location, represents the radar wavelength, represents the vertical baseline length, represents the distance from the satellite antenna center to the target geographical location, represents the unwrapped absolute phase obtained after unwrapping and absolute phase conversion of the understory terrain phase, represents the incident angle.

[0030] Further, the specific process of step 5 includes:

[0031] Step 5.1, construct a regional network adjustment model considering orbital system errors and residual forest height errors, expressed as:

[0032] ;

[0033] Wherein, is the observation vector of the elevation difference between the initial understory terrain and the spaceborne lidar control points at the same geographical location; is the design matrix composed of range and azimuth coordinates, and X is the matrix of coefficients to be estimated; is a boolean matrix used to control whether to estimate the residual forest height error, 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 spaceborne lidar control points;

[0034] Step 5.2, for the initial understory terrain obtained from each InSAR scene, select spaceborne lidar elevation points evenly distributed within the image coverage as control points to provide absolute elevation constraints for the block adjustment model:

[0035] ;

[0036] In the formula, represents the elevation of the control point ; represents the elevation value of the control point at the corresponding geographical location in the initial understory terrain; and represent the pixel coordinates in the range direction and azimuth direction of the control point at the corresponding geographical location in the SAR image coordinate system. represents the sum of the corresponding orbit error and residual forest height error; in the case of control points, the matrix expression of the adjustment model is:

[0037] ;

[0038] In the formula, is the vector of 6 coefficients in the orbit systematic error parameter vector , is the residual forest height error of the control points, is the observation error corresponding to the m control points;

[0039] Step 5.3, for any initial understory terrain, traverse the results of other initial understory terrains in the geographic coordinate system, find the data with overlapping areas, and then evenly select several tie points within the overlapping areas; based on the equality of the understory terrain elevations obtained from different InSAR data at the same geographical location, establish elevation relative constraint relationships for the block adjustment model; for any tie point establish the following functional relationship:

[0040] ;

[0041] In the formula, J and K respectively represent the results of two initial understory terrain scenes with overlapping areas; and respectively represent the elevation values of the initial understory terrain at the tie point at the corresponding geographical location in J and K, and respectively represent the connecting points after adjustment the sum of the orbital error and the residual forest height error corresponding to the geographical positions in J and K;

[0042] Step 5.4, construct the covariance function;

[0043] First, perform a block adjustment using all spaceborne lidar control points to fit , to estimate the empirical variance and covariance:

[0044] ;

[0045] In the formula, is the empirical variance, is the empirical covariance of the control point pairs with a mutual distance of ; is the total number of control points, and represent different control points in the control point pairs with a mutual distance of ; is the number of control point pairs with a mutual distance of , and are the fitting residuals of different control points , ;

[0046] Then, by selecting the empirical covariance function , and according to the calculated empirical variance and multiple empirical covariances to solve the unknown parameters in the empirical covariance function;

[0047] Step 5.5, according to the solved empirical covariance function , adopt the least squares collocation adjustment method to solve the orbital systematic error parameter vector , the residual forest height signals of the pixels covered by the control points and the pixels not covered by the control points and , respectively:

[0048] ;

[0049] In the formula, , is the covariance matrix of the observation errors, represents the covariance matrix between the observed pixels of the control points, represents the covariance matrix between the observed pixels of the control points and the non-observed pixels, calculated from the empirical covariance function ; , , The orbital system error parameter vectors respectively , the residual forest height error vector of the control point pixels and the residual forest height error vector of the uncovered pixels ;

[0050] Step 5.6. According to the solved orbital system error parameter vectors, the residual forest height signals of the control point pixels, and the residual forest height signals of the uncovered pixels, perform adjustment and correction on the initial understory terrain to obtain the final understory terrain corresponding to the current image, expressed as:

[0051] ;

[0052] In the formula, represents the final understory terrain of the th scene, represents the initial understory terrain of the th scene, represents the solved value of the orbital system error parameter vector of the th scene, is the residual forest height error vector of the uncovered pixels of the th scene.

[0053] Furthermore, select the following Gaussian - Markov model to construct the empirical covariance function:

[0054] ;

[0055] In the formula, represents the distance between two pixels, is the empirical covariance function related to the distance , and are both unknown parameters to be solved in the empirical covariance function.

[0056] Furthermore, the specific process of step 7 includes:

[0057] First, extract all the understory terrain elevation values and coherence amplitude values at the same geographical coordinate position, and calculate the weight factor corresponding to each elevation:

[0058] ;

[0059] In the formula, represents the weight corresponding to the th elevation value , represents the coherence amplitude of the full aperture, is the incident angle, represents the radar wavelength; is the vertical baseline length; represents the standard deviation of the relative elevation error;

[0060] Then, all elevation points are weighted and averaged to obtain the final understory terrain elevation :

[0061] ;

[0062] In the formula, represents the number of all understory terrain elevation values at the same geographical coordinate position.

[0063] A large-scale understory terrain mapping device based on low-frequency InSAR, comprising:

[0064] A data preprocessing module, configured to: perform registration, spectral filtering, sub-aperture decomposition, interference, and flat-earth phase removal processing on the low-frequency InSAR master-slave image pair of the large-scale understory terrain inversion area obtained, to obtain a plurality of sub-view image pairs;

[0065] An initial understory terrain estimation module, configured 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 phase unwrapping, absolute phase conversion, and phase-height conversion on each understory terrain phase to obtain the initial understory terrain;

[0066] A regional network adjustment module, configured to: construct a regional network adjustment model considering the InSAR orbit system error and the residual forest height error; use the spaceborne lidar ground elevation points as control points, and select evenly distributed tie points in the image overlap area for constraint; construct a covariance function and solve the adjustment model by the least squares collocation method to obtain the orbit system error parameters and the residual forest height error, and then correct the initial understory terrain to obtain the final understory terrain of each single scene;

[0067] A fusion and mosaicking module, configured to: fuse the understory terrains of multiple scenes of the large-scale forest to obtain a large-scale understory terrain mapping product.

[0068] Beneficial effects

[0069] The large-scale underforest terrain mapping method based on low-frequency InSAR technology of the present invention realizes obtaining the underforest terrain only by using InSAR data and a small number of control points: Firstly, low-frequency InSAR data of the large-scale underforest terrain inversion area is obtained, and the master image and the slave image of each registered InSAR interferogram are respectively subjected to sub-aperture decomposition; Secondly, operations such as interference, flat-earth phase removal, and filtering are respectively performed on the sub-interferograms after sub-aperture decomposition to obtain multiple sub-view interferograms; Then, the phase of the sub-view interferogram is fitted by a straight line to obtain the underforest terrain phase; The underforest terrain phase is converted to elevation to obtain the initial underforest terrain; Then, considering the change trend of the InSAR orbit error and the statistical law of the residual forest height error, a regional network adjustment model based on least squares registration is constructed, and tie points and spaceborne lidar 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 underforest terrain to obtain the large-scale underforest terrain.

[0070] The beneficial effects of this method are as follows: A technology for obtaining the underforest terrain phase based on sub-aperture decomposition is constructed, overcoming the limitation of the prior art that depends on fully polarized low-frequency InSAR data; A method for regional network adjustment of the underforest terrain considering systematic errors such as orbits and residual forest height errors is proposed, solving problems such as extremely few control points in the bare area and unreliable tie points without considering the penetration depth difference in the prior methods; A robust underforest terrain result fusion strategy under different map sheets and multiple coverages is established, realizing the automatic generation of large-scale underforest terrain products. In terms of accuracy, the underforest terrain products obtained by the method of the present invention are superior to the existing publicly available terrain products, and it is an effective method for autonomous and controllable, large-scale, and high-precision underforest terrain mapping. Brief Description of the Drawings

[0071] Figure 1 It is a flow chart of the method described in the embodiment of the present invention.

[0072] Figure 2 It is the location of the verification area selected by the present invention and the coverage range of the low-frequency InSAR image interferogram.

[0073] Figure 3 It is a graph of variance, covariance distribution, and fitted empirical covariance function.

[0074] Figure 4 It is a difference graph between the initial underforest terrain, the underforest terrain before regional network adjustment, and the underforest terrain after regional network adjustment and the reference DEM respectively. Figure 4 (a) represents the initial underforest terrain in the geographic coordinate system, Figure 4 (b) represents the difference graph between the initial underforest terrain before regional network adjustment in the geographic coordinate system and the reference DEM, Figure 4(c) represents the difference map between the final understory terrain after regional network adjustment in the geographic coordinate system and the reference DEM.

[0075] Figure 5 They are the difference maps between the traditional InSAR DEM products and the understory terrain obtained by the present invention and the airborne DTM. Figure 5 (a) represents the difference map between the DEM product obtained by the existing traditional InSAR technology in the geographic coordinate system and the airborne DTM, Figure 5 (b) is the difference map between the understory terrain obtained by the present invention and the airborne DTM. Detailed implementation manners

[0076] The embodiments of the present invention will be described in detail below. Based on the technical solutions of the present invention, the detailed implementation manners and specific operation processes are given, and the technical solutions of the present invention are further explained and illustrated.

[0077] In order to better illustrate the methods and steps of the present invention, in a test area located in the southern part of Guangdong Province, China, the LT-1 bistatic L-band InSAR data and the spaceborne ICESat-2 data are used to further elaborate on the present invention in detail; note that the specific implementation described here is only used to explain the present invention and is not used to limit the present invention.

[0078] This embodiment provides a method for large-scale understory terrain mapping based on low-frequency InSAR technology. Referring to Figure 1 as shown, it includes the following steps:

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

[0080] In this implementation example, the LT-1 bistatic InSAR data of the test area located in the southern part of Guangdong Province, China is selected. Four adjacent orbit data are selected, and each orbit contains 5 InSAR data, for a total of 20 interferometric pairs. The location of the test area is as Figure 2 shown; the acquired data is subjected to registration and spectral filtering processing, which is a prior art method. For references, see: Jin Guowang, Xu Qing, Zhang Hongmin. Synthetic Aperture Radar Interferometry [M]. National Defense Industry Press, 2014. This embodiment will not elaborate in detail. Then, perform one-dimensional fast Fourier transform on each pair of registered master and slave SAR images respectively to convert them to the azimuth spectral domain; then calculate the imaging azimuth angle range corresponding to the required sub-apertures according to the required number of sub-apertures (greater than 2); finally, perform inverse fast Fourier transform on the spectra corresponding to all sub-aperture azimuth angle ranges to obtain sub-view interferometric pairs.

[0081] Step 2: Perform operations such as interference, flat-earth phase removal, and phase adaptive filtering on each sub-view interference pair to obtain a number of sub-view interferograms.

[0082] Perform interference on each sub-aperture image pair, and then calculate the flat-earth phase using the orbital parameters of the SAR satellite:

[0083] ;

[0084] where represents the radar wavelength, represents the length of the parallel baseline; in order to suppress the influence of noise on the interference phase, the phase after flat-earth removal needs to be filtered. The commonly used filtering method is the adaptive filtering method, which is a mature existing technology and will not be elaborated in detail in this embodiment. Since the bistatic InSAR data is not affected by temporal decorrelation and atmospheric delay, the phase after flat-earth removal and filtering can be considered to contain only the topographic phase, forest signal phase, and system error phases such as orbit.

[0085] Step 3: Expand the complex coherence coefficients at the same positions of the sub-view interferograms within the unit circle of the complex plane, and use the straight-line equation to fit the phases of different sub-view interferograms. Select the phase at the intersection position with the unit circle as the under-forest topographic phase;

[0086] The sub-interferogram consists of two parts: interference amplitude and phase, and is usually expressed in complex format; for the sub-view interferograms of the same interference pair, it is necessary to extract the complex coherence coefficients at the same positions of each sub-view interferogram pixel by pixel and expand them into multiple points within the unit circle of the complex plane; then use a straight line to fit these points:

[0087] ;

[0088] where represents the real part of the complex coherence coefficient, represents the imaginary part of the complex coherence coefficient, and are the coefficients to be determined. Then, calculate the phases and corresponding to the two intersection points of the straight line and the unit circle of the complex plane, and judge the under-forest topographic phase based on the following criterion:

[0089] ;

[0090] Step 4: Unwrap the under-forest topographic phase, perform absolute phase conversion, phase-height conversion, etc. to obtain the initial under-forest topography;

[0091] The wrapped forest terrain phase obtained in step 3 has a value between [-π, π]. In this embodiment, the minimum cost flow method is used to unwrap the phase to obtain the unwrapped relative phase. Since the unwrapped phase is the relative phase with respect to the unwrapping reference point and cannot be directly used for phase-height conversion, control points need to be selected and robust least squares is used to convert it to the absolute phase. The above operations are existing mature technical methods and will not be elaborated in detail in this embodiment. Finally, the absolute phase is converted to elevation using the following formula:

[0092] ;

[0093] In the formula, represents the radar wavelength, represents the vertical baseline length, represents the distance from the satellite antenna center to the target point, represents the unwrapped absolute phase obtained after unwrapping and absolute phase conversion of the forest terrain phase, represents the incident angle. The initial forest terrain is as shown in Figure 4 (a).

[0094] Step 5: Based on the variation trend of the InSAR orbital system error and the statistical law of the residual forest height error, construct a block adjustment model considering the residual forest height error; use the spaceborne lidar ground elevation points as control points and select evenly distributed tie points in the overlapping areas for constraint; construct a covariance function and solve the adjustment model in the least squares collocation method; finally, remove the orbital system error and the residual forest height error from the initial forest terrain to obtain the single-scene forest terrain;

[0095] Step 5.1: Construct a block adjustment model considering the orbital system error and the residual forest height error, expressed as:

[0096] ;

[0097] In the formula, is the observation vector of the elevation difference between the initial forest terrain and the spaceborne lidar control points at the same geographical location; is the design matrix composed of the range and azimuth coordinates, X is the coefficient matrix to be estimated, is a Boolean matrix used to control whether to estimate the residual forest height error (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 spaceborne lidar control points.

[0098] Step 5.2: For the initial forest terrain obtained by each scene of InSAR, select relatively evenly distributed The elevation points of the spaceborne lidar are used as control points to provide absolute elevation constraints for the model. Since it is necessary to model the residual forest height error, control points located in the forest area need to be added. The following functional relationships can be established for each control point:

[0099] ;

[0100] In the formula, represents the elevation of the control point ; represents the elevation value of the corresponding geographical location of the control point on the initial understory terrain; and represent the pixel coordinates in the range and azimuth directions of the corresponding geographical location of the control point in the SAR image coordinate system. represents the sum of the corresponding orbital error and residual forest height error; in the case of control points, the matrix expression corresponding to formula (5) is:

[0101] ;

[0102] In the formula, is the vector of 6 coefficients of the orbital systematic error parameters in, is the residual forest height error of the control points, is the residuals corresponding to the control points;

[0103] Step 5.3. For any initial understory terrain, traverse the results of other initial understory terrains in the geographic coordinate system, find the data with overlapping areas, and then uniformly select positions with relatively flat terrain and high coherence in the overlapping areas as connection points; based on the fact that the understory terrain elevations obtained from different InSAR data at the same geographical location are equal, establish an elevation relative constraint relationship for the model; for any connection point the following functional relationship is established:

[0104] ;

[0105] In the formula, J and K respectively represent the results of two scenes of initial understory terrain with overlapping areas; and respectively represent the initial understory terrain elevation values of the connection point j at the corresponding geographical location in J and K, and respectively represent the sum of the orbital error and residual forest height error of the connection point j at the corresponding geographical location in J and K after adjustment.

[0106] Step 5.4, construct the covariance function. First, perform a regional network adjustment based on the least squares theory using all spaceborne lidar control points to fit , to estimate the empirical variance and covariance:

[0107] ;

[0108] 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 the Euclidean distance and rounded; and are the fitting residuals of different control points, that is, the residual matrix expression of m control points is , where are the unknown parameters solved by the first adjustment fitting. Then, calculate the empirical covariance under different distance intervals , and select the following Gaussian-Markov model as the empirical covariance function according to its distribution characteristics:

[0109] ;

[0110] In the formula, d represents the distance between two pixels, the correlation length , the variance are unknown parameters, and the variance and covariance calculated by formula (9) need to be estimated by the least squares method. The distribution of the estimated variance and covariance is as shown by the dots in Figure 3 , and the fitted empirical covariance function is as shown by the curve in Figure 3 .

[0111] Step 5.5, through the empirical covariance function established by formula (10) and the least squares collocation adjustment method, the orbital system error parameters , the residual forest height signals of the control point pixels and the uncovered pixels and are:

[0112] ;

[0113] In the formula, , is the covariance matrix of the observation error (residual), represents the covariance matrix between the observation pixels of the control points, obtained from the empirical covariance calculated by formula (9); It represents the covariance matrix of the observed pixels and non-observed pixels of the control points, which is calculated by the empirical covariance function, i.e., Equation (10). For the model solution process of the least squares collocation method in the above equation, reference can be made to: Cui Xizhang, Yu Zongchou, Tao Benzao, etc. Generalized Surveying Adjustment (New Edition) [M]. Wuhan Technical University of Surveying and Mapping Press, 2001., which will not be elaborated in detail in this embodiment.

[0114] Step 5.6, according to the solved orbital system error parameter vector, the residual forest height signal of the control point pixels, and the residual forest height signal of the non-covered pixels, perform adjustment and correction on the initial understory terrain to obtain the final understory terrain corresponding to the current image, which is expressed as:

[0115] ;

[0116] In the formula, represents the final understory terrain of the th scene, represents the initial understory terrain of the th scene, represents the solved value of the orbital system error parameter vector of the th scene, is the residual forest height error vector of the non-covered pixels of the th scene.

[0117] Step 6, fuse and mosaic the single-scene understory terrains obtained in Step 5 to obtain a large-scale understory terrain product.

[0118] First, extract all the understory terrain elevation values and coherence amplitude values at the same geographical coordinate position, and calculate the weight factor corresponding to each elevation:

[0119] ;

[0120] In the formula, represents the weight corresponding to the th elevation value, represents the coherence amplitude of the full aperture, is the incident angle, is the vertical baseline length; represents the standard deviation of the relative elevation error.

[0121] Then, perform weighted averaging on all elevation points to obtain the final understory terrain elevation:

[0122] ;

[0123] In the formula, represents the number of all understory terrain elevation values at the same geographical coordinate position.

[0124] 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 terrains of different orbits, which is due to the fact that the data of different orbits are obtained at different times, and the accuracies of satellite orbit determination and baseline estimation are inconsistent. Therefore, block adjustment needs to be used to process it. Figure 4 (b) is the difference map between the initial understory terrain before block adjustment and the reference DEM. Figure 4 (c) is the difference map between the final understory 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 orbits. In order to more prominently highlight the advantages of using low-frequency InSAR for large-scale understory terrain mapping in this embodiment, the airborne high-precision LiDAR DTM is used as a reference and compared with the existing publicly available high-frequency InSAR DEM product (TanDEM-X DEM). Figure 5 (a) is the difference map between the traditional InSAR DEM and the airborne DTM. Figure 5 (b) is the difference map between the understory terrain of this embodiment and the airborne DTM. It is found that Figure 5 (b) the difference with the LiDAR DTM is smaller. The RMSE decreases from 5.87 m to 3.09 m, and the MAE decreases from 4.59 m to 0.34 m, indicating that the accuracy of the understory terrain mapping results using low-frequency InSAR in this embodiment is higher.

[0125] The above embodiments are the preferred embodiments of the present application. Those of ordinary skill in the art can also make various transformations or improvements based on this. Without departing from the general concept of the present application, these transformations or improvements should all fall within the scope protected 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 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.

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