Tree height estimation method and system based on SAR scattering intensity correction in terrain relief area

Abstract

The application discloses a terrain undulating area tree height estimation method and system based on SAR scattering intensity correction, comprising converting the irradiation area of SAR image pixels in a geographic coordinate system to a SAR pixel coordinate system by using inverse distance weighting; combining the SAR image resolution and the irradiation area of SAR image pixels in the SAR pixel coordinate system to calculate the SAR backscattering intensity after irradiation area correction; taking the measured forest height as prior knowledge, estimating the forest height based on the random volume-surface model, obtaining the best fitting coefficient of the random volume-surface model by using a nonlinear regression analysis method, and then obtaining the forest height of the research area. The application solves the problem of the existing forest height estimation technology based on single SAR backscattering, that is, the estimation deviation caused by insufficient compensation of the terrain effect in hilly, mountainous and other undulating terrains, and can realize large-scale and high-robustness forest height mapping.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of radar remote sensing and forest parameter inversion, and particularly relates to a terrain undulation region tree height estimation method and system based on SAR scattering intensity correction. BACKGROUND

[0002] Forest height is one of the most basic structural parameters of forest ecosystem, and its accurate estimation is of great significance for comprehensively understanding forest structure and function and realizing sustainable management of natural resources. As an active microwave remote sensing system, synthetic aperture radar (SAR) can penetrate the forest canopy and obtain forest vertical structure information. The widely studied polarimetric interferometric SAR (PolInSAR) and tomographic SAR (TomoSAR) technologies based on SAR interference characteristics can extract more detailed forest vertical structure distribution by introducing multi-polarization and multi-orbit observation information, and then realize high-precision forest estimation. However, they rely on repeated orbit interferometric measurement, which is easily affected by time de-coherence in forest-covered areas, and even causes estimation failure due to serious coherence loss.

[0003] Unlike the former two technologies, considering that SAR backscattering intensity has significant correlation with biophysical parameters such as forest height and biomass, the technology based on single SAR backscattering intensity can establish statistical and semi-empirical models for different forest scenes through single transit observation, and then can carry out forest height estimation. Moreover, it has low computational complexity, is suitable for spatial large-scale processing, and has strong multi-source data fusion capability. However, due to the side-looking imaging characteristics of SAR sensors, the measured signal is easily affected by terrain undulation. The same ground object will show different signal characteristics on the SAR image due to the difference in local terrain. Therefore, when using SAR to interpret ground objects or quantitatively invert parameters, the influence of terrain is an unavoidable problem. Especially for mountainous forests with complex terrain, the statistical and semi-empirical models based on single SAR backscattering intensity have a significant decrease in accuracy in complex terrain areas, which is specifically manifested as follows: 1) irradiation area change, and the traditional homomorphic correction method (such as local incidence angle-based radiation correction) is difficult to accurately reflect the heteromorphic relationship between geography and SAR slant range geometry; 2) angle variation effect, and terrain change will change the forest scattering mechanism and microwave penetration path. However, the existing methods cannot compensate for the effect of terrain, and cannot simultaneously consider the irradiation area change and angle variation effect. In summary, the existing technology cannot realize high-precision and robust forest height estimation in hilly areas, and a new method that fuses terrain slope correction and scattering mechanism optimization is urgently needed. SUMMARY

[0004] In order to solve the existing forest height estimation technology based on single SAR backscattering, the estimation deviation problem caused by insufficient compensation of terrain effect in hilly, mountainous and other undulating terrains, and realize large-scale and high-robustness forest height mapping, the present application proposes a terrain undulating area tree height estimation method based on SAR scattering intensity correction, which comprises the following steps:

[0005] Step 1, acquiring HV polarized SAR single view complex image, digital elevation model and forest / non-forest classification data covering the research area;

[0006] Step 2, calculating the irradiation area of SAR image pixels in the geographic coordinate system by using SAR system orbit parameters, SAR image resolution and DEM data;

[0007] Step 3, according to the conversion relationship between the geographic coordinate system and the SAR pixel coordinate system, the irradiation area of SAR image pixels in the SAR pixel coordinate system is obtained by using inverse distance weighting;

[0008] Step 4, combining the SAR image resolution and the irradiation area of SAR image pixels in the SAR pixel coordinate system, the SAR backscattering intensity after irradiation area correction is calculated;

[0009] Step 5, taking the measured forest height as prior knowledge, combining the SAR backscattering intensity after irradiation area correction, the best fitting coefficient of the random volume-ground model is obtained by using nonlinear regression analysis method, and then the forest height of the research area is obtained.

[0010] Further, in the geographic coordinate system, the irradiation area of each SAR image pixel in step 2 is calculated according to the following formula:

[0011] (1)

[0012] In the formula, is the irradiation area of SAR image pixels in the geographic coordinate system, is the resolution of SAR image azimuth direction, is the resolution of SAR image slant range direction, is the projection angle cosine, which is obtained by the dot product of the terrain surface normal vector and the radar slant range plane normal vector:

[0013] (2)

[0014] (3)

[0015] (4)

[0016] In the formula, is the terrain surface normal vector, is the change rate of elevation in the range direction of the SAR image, is the change rate of elevation in the range direction of the SAR image, is the radar slant range plane normal vector, is the velocity vector of the sensor, is the slant range vector.

[0017] Further, the step 3 uses the SAR geo-coding lookup table to find the SAR image pixel After the irradiance area in the geographic coordinate system is converted to the SAR pixel coordinate system, the irradiance area in the radar pixel coordinate is distributed to the 4 neighboring pixels around the SAR image pixel , , , wherein, , , , ; assuming that any pixel on the SAR image is assigned to N irradiance areas in the geographic coordinate system, the corresponding irradiance areas in the geographic coordinate system are weighted and summed to obtain the irradiance area of the pixel in the radar pixel coordinate :

[0018] (5)

[0019] (6)

[0020] (7)

[0021] (8)

[0022] In the formula, is the kth irradiance area in the geographic coordinate system, is the weight of the kth irradiance area, is the SRA image pixel corresponding to the kth irradiance area, is the distance between the pixel and the four surrounding integer index pixels, floor represents the down rounding, and ceil represents the up rounding, represents the azimuth direction integer index obtained by down rounding , respectively represent the azimuth direction integer index obtained by up rounding , represents the azimuth direction integer index obtained by down rounding , Indicates to The directional integer index obtained by rounding up.

[0023] Furthermore, the SAR backscattering intensity after irradiation area correction in step 4 is calculated by the following formula:

[0024] (9)

[0025] In the formula, This represents the SAR backscattering intensity after irradiation area correction. This represents the resolution of the SAR image in the azimuth direction. This represents the slant range resolution of the SAR image. The SLC intensity value is after radiometric calibration. SAR image pixels Irradiation area in SAR pixel coordinates.

[0026] Furthermore, the random volume-surface model in step 5 is represented as follows:

[0027] (10)

[0028] In the formula, Indicates the SAR backscattering intensity. This is the height value along the direction perpendicular to the inclined ground surface. This refers to the surface elevation along a direction perpendicular to the inclined surface. This represents the actual forest height. This indicates the distance from the terrain to the slope. The resulting "pseudo" forest height, The scattering intensity is for the volume scattering mechanism. The scattering intensity is for even-order scattering mechanisms. Extinction coefficient, It is a local incident angle.

[0029] Substituting the SAR backscattering intensity corrected for irradiance area obtained in step 4 into formula (10), we establish its relationship with forest height, namely:

[0030] (11)

[0031] in:

[0032] (12)

[0033] In the formula, For predictor variables, For response variables, represents the fitting coefficient.

[0034] The SAR backscattering intensity after irradiation area correction is converted to a geographic coordinate system by using a geographic coding lookup table, the non-forest area in the backscattering intensity image in the geographic coordinate system is masked by using forest / non-forest classification data, the "pseudo" forest height value caused by the slope of the terrain distance is calculated according to part of the measured forest height value, which is used as a priori value, combined with the backscattering intensity image of the forest area, the best fitting coefficient of formula (11) is obtained by using nonlinear regression analysis technology, and the best fitting coefficient and the SAR backscattering intensity of the study area after irradiation area correction are substituted into formula (11), and then the forest height of the study area is obtained.

[0035] The application also provides a terrain undulating area tree height estimation system based on SAR scattering intensity correction, which is used for realizing the terrain undulating area tree height estimation method based on SAR scattering intensity correction.

[0036] Moreover, the system comprises a processor and a memory, the memory is used for storing program instructions, and the processor is used for calling the stored instructions in the memory to execute the terrain undulating area tree height estimation method based on SAR scattering intensity correction.

[0037] Alternatively, the system comprises a readable storage medium, and the readable storage medium stores a computer program, and the computer program is executed to realize the terrain undulating area tree height estimation method based on SAR scattering intensity correction.

[0038] Compared with the prior art, the application has the following advantages:

[0039] 1) The application considers the difference of radar wave irradiation area caused by the change of terrain slope and the angle dependence of SAR forest backscattering correction in the vertical direction, on the basis of the ground irradiation area related amplitude correction in the first stage, the forest height is further extracted by using the slope correction backscattering semi-empirical method of the ground scattering mechanism in the second stage, the influence of SAR signal distortion of the forest area caused by the terrain undulation is weakened, the interference of the change of terrain slope on the SAR forest height estimation performance is suppressed, and the applicability in the complex terrain area is improved;

[0040] 2) The application takes the spatially distributed sparse LiDAR high-precision tree height data as a key "ground truth" anchor point, realizes the collaborative fusion strategy of "SAR continuous surface + LiDAR discrete point", and can generate an enhanced forest height product with the characteristics of high LiDAR point precision and wide SAR surface coverage;

[0041] 3) The application only depends on single scene SAR backscattering intensity image, avoids the complex interferometric analysis process, and is expected to provide key algorithm support for the production of business, large-scale and high-precision forest height products.

[0042] 4) The present application breaks the dilemma of the traditional method being limited in the area with significant terrain undulations, and can effectively improve the forest height mapping effect and spatial continuity in mountainous, hilly and other areas. BRIEF DESCRIPTION OF DRAWINGS

[0043] In order to more clearly illustrate the technical solutions in the present application or prior art, the following will briefly introduce the drawings needed to be used in the embodiments or prior art description. Obviously, the drawings described below are some embodiments of the present application, and all other embodiments obtained by those of ordinary skill in the art without creative labor based on these drawings also belong to the protection scope of the present application.

[0044] Figure 1 is a flow chart of the terrain undulating area tree height estimation method based on SAR scattering intensity correction of the embodiments of the present application.

[0045] Figure 2 is a contrast chart of the airborne SAR backscattering intensity before and after correction of the embodiments of the present application, wherein (a) is the airborne SAR backscattering intensity chart before correction, and (b) is the airborne SAR backscattering intensity chart after correction.

[0046] Figure 3 is a statistical distribution histogram of the airborne SAR backscattering intensity before and after correction of the embodiments of the present application.

[0047] Figure 4 is a relationship chart of the response variable and the prediction variable in the regression model of the embodiments of the present application, wherein (a) is the relationship chart of the prediction variable and the response variable before correction, and (b) is the relationship chart of the prediction variable and the response variable after correction.

[0048] Figure 5 is a contrast chart of the forest height product before and after correction and the LiDAR forest height reference data of the embodiments of the present application, wherein (a) is the forest height product before correction, (b) is the forest height product after correction, (c) is the LiDAR forest height reference data, and (d) is the RMSE distribution of the forest height before and after correction and the LiDAR forest height reference data. DETAILED DESCRIPTION

[0049] In order to make the purpose, technical solutions and advantages of the present application more clear, the technical solutions of the present application will be further described below in combination with the drawings and embodiments. Obviously, the described embodiments are some embodiments of the present application, but not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative labor also belong to the protection scope of the present application.

[0050] Embodiment 1

[0051] As Figure 1 shown, the embodiment of the present application provides a terrain relief area tree height estimation method based on SAR scattering intensity correction, comprising the following steps:

[0052] Step 1, obtain the HV polarized SAR single view complex (Single-Look Complex, SLC) image covering the study area, digital elevation model (Digital Elevation Model, DEM) and forest / non-forest classification data.

[0053] The SAR data covering the study area is collected by the UAVSAR (Uninhabited Aerial Vehicle Synthetic Aperture Radar) developed by the Jet Propulsion Laboratory (JPL) of the National Aeronautics and Space Administration (NASA), which adopts L-band (wavelength about 24 cm) and has full polarization observation capability, can penetrate the vegetation canopy and is sensitive to the surface and subsurface structure. The digital elevation model uses the FABDEM_V1-2 which removes forests and buildings, and the forest / non-forest classification data uses the national land cover dataset (NLCD).

[0054] Step 2, calculate the irradiation area of SAR image pixels in the geographic coordinate system using SAR system orbit parameters, SAR image resolution and DEM data.

[0055] In the geographic coordinate system, the irradiation area of each SAR image pixel is calculated as follows:

[0056] (1)

[0057] In the formula, is the irradiation area of the SAR image pixel in the geographic coordinate system, is the resolution of the SAR image in the azimuth direction, is the resolution of the SAR image in the range direction, is the cosine of the projection angle, which is obtained by the dot product of the terrain surface normal vector and the radar range plane normal vector:

[0058] (2)

[0059] (3)

[0060] (4)

[0061] where, is the terrain surface normal vector, is the rate of change of elevation in the azimuth direction of the SAR image, is the rate of change of elevation in the range direction of the SAR image, is the radar slant range plane normal vector, is the velocity vector of the sensor, is the slant range vector.

[0062] Step 3, according to the conversion relationship between the geographic coordinate system and the SAR pixel coordinate system, the irradiation area of the SAR image pixel in the SAR pixel coordinate system is obtained by inverse distance weighting.

[0063] The geographic coordinate system regards the Earth as an approximate sphere, takes the center of the Earth as the origin, takes the equatorial plane as the reference surface, and takes the plane of the prime meridian as the initial meridian plane. The coordinates of each point on the Earth are determined by the intersection of the two planes and the axis of rotation of the Earth. Among them, the longitude represents the angle between a point and the plane of the prime meridian, and the latitude represents the angle between a point and the equatorial plane.

[0064] The SAR pixel coordinate system takes the upper left corner of the SAR image as the origin, the horizontal direction as the x-axis, the right as the positive direction, and the vertical direction as the y-axis, the downward as the positive direction. Each pixel point in the coordinate system has a unique coordinate value (slant range coordinate, azimuth coordinate), which is used to represent its position in the radar image.

[0065] The conversion relationship between the SAR image pixel in the geographic coordinate system and the SAR pixel coordinate system can be obtained through the SAR geographic encoding lookup table. The SAR image pixel After the irradiation area in the geographic coordinate system is converted to the SAR pixel coordinate system, it will be deformed, so the irradiation area in the radar pixel coordinate system is distributed to the four neighboring pixels around the SAR image pixel , , , where, , , , . Assuming that any pixel on the SAR image will be distributed to N irradiation areas in the geographic coordinate system, the corresponding irradiation areas in the geographic coordinate system are weighted and summed to obtain the irradiation area of the pixel in the radar pixel coordinate system:

[0066] (5)

[0067] (6)

[0068] (7)

[0069] (8)

[0070] where, is the kth irradiance area in the geographic coordinate system, is the weight of the kth irradiance area, is the SRA image pixel corresponding to the kth irradiance area, is the distance between the pixel and the surrounding integer index pixels, floor represents the down rounding, and ceil represents the up rounding, represents the azimuth integer index obtained by down rounding , respectively represent the azimuth integer index obtained by up rounding , represents the azimuth integer index obtained by down rounding , represents the azimuth integer index obtained by up rounding .

[0071] Step 4, combine the SAR image resolution and the irradiance area of the SAR image pixel in the SAR pixel coordinate system to calculate the SAR backscatter intensity after irradiance area correction.

[0072] The SAR backscatter intensity after irradiance area correction is calculated by the following formula:

[0073] (9)

[0074] where, is the SAR backscatter intensity after irradiance area correction, is the resolution of the SAR image in the azimuth direction, is the resolution of the SAR image in the range direction, is the SLC intensity value after radiation calibration, is the SAR image pixel irradiance area under the SAR pixel coordinate.

[0075] Figure 2 is the comparison chart of the airborne SAR backscatter intensity before and after irradiance area correction, wherein (a) is the airborne SAR backscatter intensity chart before irradiance area correction, and (b) is the airborne SAR backscatter intensity chart after irradiance area correction. From Figure 2It can be seen from Fig. 5 (a) that the areas facing the radar incidence direction appear brighter than the areas facing away from the radar incidence direction. This is because the irradiation area is calculated based on the Earth ellipsoid model, ignoring the non-uniform correspondence between the SAR slant-range image and the geographic coordinate space, resulting in insufficient terrain correction. In addition, the airborne SAR has a wide viewing angle range of about 20°~70°, and the perspective shortening effect of the near-range slope facing the radar is more obvious, appearing brighter, while the far-range slope facing away from the radar is relatively less affected, so there are more shadow areas on the rear slope. From Fig. 5 (b), it can be seen that after the irradiation area correction, the difference in radar scattering intensity between the upslope and the downslope becomes smaller, and the texture characteristics are more uniform, especially in the near range, which has been significantly improved. Figure 2 From Fig. 5 (b), it can be seen that after the irradiation area correction, the difference in radar scattering intensity between the upslope and the downslope becomes smaller, and the texture characteristics are more uniform, especially in the near range, which has been significantly improved. Figure 3 The statistical distribution of the backscattering intensity of the airborne SAR before and after correction is shown, where the variance of the backscattering intensity before and after correction is 8.45 and 4.35 respectively, and the variance is reduced by 48.5%, which shows that after the correction operation, the brightness difference of the image pixel points is smaller.

[0076] Step 5: Using the measured forest height as prior knowledge, combining the SAR backscattering intensity after irradiation area correction, and using nonlinear regression analysis method to obtain the best fitting coefficient of the random volume-ground model, and then obtaining the forest height of the study area.

[0077] In the terrain undulating area, the slope change will change the penetration path of the microwave signal in the forest scene, and thus the scattering mechanism changes. At this time, the scattering process of the forest scene can be described by the slope-based random volume over ground (Slope-RVoG) model, which can be expressed as:

[0078] (10)

[0079] In the formula, represents the SAR backscattering intensity, is the height value along the direction perpendicular to the inclined surface, is the surface elevation along the direction perpendicular to the inclined surface, is the actual forest height, represents the “pseudo” forest height caused by the terrain distance to the slope , is the scattering intensity of the volume scattering mechanism, is the scattering intensity of the double scattering mechanism, is the extinction coefficient, which describes the one-way power loss of the microwave in the forest canopy, is the local incidence angle.

[0080] The SAR backscattering intensity after the irradiation area correction obtained in step 4 is substituted into formula (10) to establish the relationship between the SAR backscattering intensity and the forest height, i.e.

[0081] (11)

[0082] wherein:

[0083] (12)

[0084] in the formula, is a predicted variable, is a response variable, is a fitting coefficient.

[0085] The SAR backscattering intensity after the irradiation area correction is converted to a geographic coordinate system by using a geographic coding lookup table. The GAMMA software is used to calculate the "pseudo" forest height value caused by the terrain distance slope according to part of the measured forest height values. After the non-forest area in the backscattering intensity image is masked by using the forest / non-forest classification data, the "pseudo" forest height value is combined as prior knowledge, and the best fitting coefficient of formula (11) is obtained by using the nonlinear regression analysis technology. The best fitting coefficient and the SAR backscattering intensity after the irradiation area correction of the research area are substituted into formula (11), and then the forest height of the research area is obtained.

[0086] In this embodiment, the regression standard error (SER) describing the dispersion degree of the observation data points around the regression line is used to evaluate the fitting degree of the regression model:

[0087] (13)

[0088] in the formula, is an actual response variable value, is a model predicted value, is a sample number, is a number of prediction variables in the model.

[0089] The smaller the SER is, the higher the fitting degree of the regression equation with the data is, and the stronger the prediction ability is. After the parameters of the nonlinear regression model are determined, the model is transformed from a general framework into a model with specific prediction ability for the current research area and the data characteristics. By applying it to the entire research area, a forest height spatial distribution map covering all forest pixels can be generated. In this embodiment, the optimal fitting coefficient is the fitting coefficient corresponding to the minimum SER.

[0090] Figure 4The relationship between the response variable and the prediction variable in the regression model and the model fitting result are shown, wherein the prediction variable in (a) is the SAR backscattering intensity of the irradiation area without correction, and the prediction variable in (b) is the SAR backscattering intensity of the irradiation area after correction. Figure 4 It can be seen that the dynamic range of the corrected prediction variable and the response variable is smaller, and the sensitivity is higher. The SER index is reduced from 4.31 to 3.40, indicating that the data fitting degree of the corrected regression equation is higher, and the prediction ability is better.

[0091] Figure 5 In (a) and (b), the forest height results estimated by using the SAR backscattering intensity of the irradiation area without correction and the SAR backscattering intensity of the irradiation area after correction are respectively represented. Figure 5 In (c), the LiDAR forest height reference data are represented. Figure 5 As can be seen from (a), (b) and (c), when the terrain undulation is not considered, the forest height estimation result presents a trend of underestimation at a short distance and overestimation at a long distance, and by using the forest height estimation strategy considering the terrain undulation proposed in the application, the obtained result is consistent with the spatial distribution of the LiDAR forest height. In order to further carry out quantitative evaluation, the LiDAR forest height values are binned and counted with 1m as a unit, and the corresponding root mean square error (RMSE) is calculated as the forest height evaluation index. From Figure 5 As can be seen from (d), the error trends of the estimation before and after correction are similar at different forest height levels. Compared with the overall average root mean square error of 6.20m before correction, the overall average root mean square error after correction is 4.49m, and the forest height estimation accuracy is improved by 27.6%.

[0092] Embodiment 2

[0093] Based on the same inventive concept, the application further provides a terrain undulation area tree height estimation system based on SAR scattering intensity correction, comprising a processor and a memory, the memory is used for storing program instructions, and the processor is used for calling the program instructions in the memory to execute the terrain undulation area tree height estimation method based on SAR scattering intensity correction.

[0094] Embodiment 3

[0095] Based on the same inventive concept, the application further provides a terrain undulation area tree height estimation system based on SAR scattering intensity correction, comprising a readable storage medium, and a computer program is stored on the readable storage medium, and the computer program is executed to realize the terrain undulation area tree height estimation method based on SAR scattering intensity correction.

[0096] In specific implementation, the method provided by the technical scheme of the present application can be automatically run by a computer software technology, and the system device of the method, such as a computer readable storage medium storing the corresponding computer program of the technical scheme of the present application and a computer device including the corresponding computer program, should also be within the protection scope of the present application.

[0097] The specific embodiments described herein are merely illustrative of the principles of the present application. Various modifications or changes in light thereof can be made by those skilled in the art to which the present application pertains without departing from the spirit of the present application or exceeding the scope of the appended claims.

Claims

1. A method for tree height estimation in a terrain relief area based on SAR scattering intensity correction, characterized in that, The method comprises the following steps: Step 1, obtaining HV polarized SAR single-view complex images covering the study area, digital elevation model (DEM) data and forest / non-forest classification data; Step 2, calculating the irradiation area of SAR image pixels in the geographic coordinate system by using the orbit parameters of the SAR system, the resolution of the SAR image and the DEM data; Step 3, obtaining the irradiation area of SAR image pixels in the SAR pixel coordinate system by using the inverse distance weighted method according to the conversion relationship between the geographic coordinate system and the SAR pixel coordinate system; Step 4, calculating the SAR backscattering intensity after irradiation area correction by combining the SAR resolution and the irradiation area of SAR image pixels in the SAR pixel coordinate system; Step 5, taking the measured forest height as prior knowledge, combining the SAR backscattering intensity after irradiation area correction, and obtaining the best fitting coefficient of the random volume-surface model by using the nonlinear regression analysis method, and then obtaining the forest height of the study area; The random volume-surface model is expressed as: （10） where denotes the SAR backscatter intensity, is the height value along the direction perpendicular to the tilted terrain, is the terrain elevation along the direction perpendicular to the tilted terrain, is the actual forest height, denotes the "false" forest height caused by the terrain distance to slope gradient, is the scattering intensity of the volume scattering mechanism, is the scattering intensity of the double scattering mechanism, is the extinction coefficient, is the local incidence angle.

2. The method of claim 1, wherein the method is based on SAR intensity correction for a terrain-undulating area tree height estimation. In step 2, the irradiation area of each SAR image pixel in the geographic coordinate system is calculated according to the following formula: （1） wherein is the irradiation area of the SAR image pixel in the geographical coordinate system, is the azimuth resolution of the SAR image, is the slant range resolution of the SAR image, is the cosine of the projection angle, which is obtained by the dot product of the terrain surface normal vector and the radar slant range plane normal vector: （2） （3） （4） wherein is the terrain surface normal vector, is the rate of change of elevation in the azimuth direction of the SAR image, is the rate of change of elevation in the range direction of the SAR image, is the radar slant range plane normal vector, is the velocity vector of the sensor, is the slant range vector.

3. The method of claim 1, wherein the method is based on SAR intensity correction for terrain-roughness regions of tree height estimation. SAR image pixels in step 3 using the SAR geocoding lookup table The irradiance area in the radar pixel coordinate is assigned to the SAR image pixel after the transformation from the geographic coordinate system to the SAR pixel coordinate system The 4 neighboring pixels around , , , where , , , Suppose that any pixel on the SAR image will be assigned to N irradiance areas in the geographic coordinate system, the corresponding irradiance areas in the geographic coordinate system are weighted and summed to obtain the pixel irradiance area in the radar pixel coordinate : （5） wherein is the kth irradiation area in the geographic coordinate system, is the weight of the kth irradiation area.

4. The method of claim 3, wherein the method is based on SAR intensity correction for terrain-roughness regions. 5 The weight in step 3 The formula for calculating the weight is as follows: （6） （7） （8） In the formula, For the SRA image pixel corresponding to the k-th irradiated area, for The distance to the surrounding integer index pixels, where floor represents rounding down and ceil represents rounding up. Indicates to The azimuth integer index obtained by rounding down. They represent respectively to The azimuth integer index obtained by rounding up. Indicates to The azimuth integer index obtained by rounding down. Indicates to The directional integer index obtained by rounding up.

5. The method of claim 1, wherein the method is based on SAR intensity correction for terrain-roughness regions of tree height estimation. In step 4, the SAR backscattering intensity after irradiation area correction is calculated according to the following formula: （9） wherein, is the SAR backscattering intensity after irradiation area correction, is the resolution of the SAR image in azimuth direction, is the resolution of the SAR image in slant range direction, is the SLC intensity value after radiometric calibration, is the SAR image pixel is the irradiation area in the SAR pixel coordinate.

6. The method of claim 1, wherein the method is based on SAR intensity correction for terrain-roughness regions. In step 5, the SAR backscattering intensity after irradiation area correction obtained in step 4 is substituted into formula (10) to establish the relationship between the SAR backscattering intensity and the forest height, that is: （11） In step 5, the SAR backscattering intensity after irradiation area correction is converted to the geographic coordinate system by using the geographic coding lookup table, the non-forest area in the backscattering intensity image in the geographic coordinate system is masked by using the forest / non-forest classification data, the "pseudo” forest height value caused by the terrain distance slope is calculated according to the measured forest height value, which is used as the prior value, the backscattering intensity image in the forest area is combined, and the best fitting coefficient of formula (11) is obtained by using the nonlinear regression analysis technology, the best fitting coefficient and the SAR backscattering intensity after irradiation area correction of the study area are substituted into formula (11), and then the forest height of the study area is obtained. （12） wherein is the predictor variable, is the response variable, is the fitting coefficient.

7. The method of claim 6, wherein the method is based on SAR intensity correction for terrain-roughness regions. The method comprises a processor and a memory, the memory is used for storing program instructions, and the processor is used for calling the program instructions in the memory to execute the method for estimating the tree height in a terrain undulating area based on SAR scattering intensity correction according to any one of claims 1-7. 8.A system for tree height estimation in a terrain-undulated area based on SAR scattering intensity correction, characterized in that, The readable storage medium comprises a computer program stored thereon, and the computer program is executed to realize the method for estimating the tree height in a terrain undulating area based on SAR scattering intensity correction according to any one of claims 1-7. 9.A system for tree height estimation in a terrain-undulated area based on SAR scattering intensity correction, characterized in that, ​

Images