Forest stock monitoring method based on SAR remote sensing quantitative model
By using an improved SAR remote sensing quantitative model, combined with topographic slope factors and multidimensional features, the problems of low efficiency, poor accuracy, and high cost in traditional forest stock volume monitoring have been solved, achieving high-precision forest stock volume monitoring and rapid processing.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-12-03
- Publication Date
- 2026-03-31
AI Technical Summary
Traditional methods for monitoring forest stock volume are inefficient, costly, and difficult to cover remote areas. Optical remote sensing is easily affected by cloud and rain weather. Airborne lidar is costly and difficult to apply on a large scale. SAR technology has large errors in tree height inversion under complex terrain. The mechanism for estimating stock volume is missing. Multi-source data fusion is inefficient and the model has poor generalization ability.
A quantitative model based on SAR remote sensing was adopted. The terrain slope factor was introduced to correct the radar incident angle and extinction path through the improved S-RVoG model. Combined with the polarimetric coherence matrix and high-precision DEM data, the forest tree height was inverted using a three-stage inversion method. A volume estimation model was constructed, and multi-dimensional features were automatically exported. The hyperparameters of the random forest were optimized, and unified data resampling and automatic coordinate transformation were achieved.
The accuracy of tree height inversion under complex terrain is improved, the error is reduced by 30%, the accuracy of volume estimation meets the forestry carbon sink measurement standards, and the processing time is greatly shortened to meet the operational timeliness requirements.
Smart Images
Figure CN121763255A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of forest resource monitoring technology, specifically a method for monitoring forest stock volume based on a quantitative SAR remote sensing model. Background Technology
[0002] Traditional forest stock volume monitoring mainly relies on ground surveys, which suffers from low efficiency, high cost, and difficulty in covering remote areas. Optical remote sensing is susceptible to interference from cloud and rain, while airborne lidar (LiDAR) is expensive and difficult to apply on a large scale. Synthetic aperture radar (SAR) has all-weather, day-and-night observation capabilities, but current SAR technology still has the following shortcomings: 1) Tree height inversion is greatly affected by terrain: The traditional polarimetric interferometric synthetic aperture radar (PolInSAR) model does not fully consider the influence of complex terrain slope, resulting in the radar incident angle and extinction path not being corrected, and the high-precision digital elevation model (DEM) data is not fused, resulting in a root mean square error (RMSE) of 3.6m for tree height inversion in steep slope areas (average slope ≥23°).
[0003] 2) Lack of a mechanism for estimating forest volume: Relying on empirical models such as logarithmic fitting of backscattering coefficients in high-density forests (volume > 200m³) is insufficient. 3 Signal saturation is prone to occur in the / ha) model, and the model parameters require extensive ground calibration. It has poor generalization ability and fails to explore the synergistic effect between polarization decomposition features and terrain parameters.
[0004] 3) Inefficient multi-source data fusion: SAR data, ground-based measured data, and UAV LiDAR data have spatiotemporal resolution and coordinate system heterogeneity. Traditional methods require manual intervention for registration and feature extraction, which is time-consuming and difficult to meet the needs of operational applications. Summary of the Invention
[0005] In view of this, the purpose of this invention is to provide a method for monitoring forest stock volume based on a quantitative model of SAR remote sensing, so as to solve the technical problems mentioned in the prior art.
[0006] The method for monitoring forest stock volume based on SAR remote sensing quantitative models includes the following steps: S1. Acquire SAR data and ground-measured data for any monitoring area, wherein the ground-measured data includes sample data, DEM data, and DSM data; S2. Perform image registration, spectral filtering, interferometric coherence, flat phase removal, terrain phase removal, and interferometric fringe pattern filtering on the SAR data to generate a polarization coherence matrix T6. S3. Using an improved S-RVoG model, a terrain slope factor is introduced to correct the radar incident angle and extinction path of the SAR data. Based on the polarization coherence matrix T6 and the DEM data, forest tree height is inverted using a three-stage inversion method. The three-stage inversion method is set as follows: based on the overall least squares method, linear fitting is performed on the complex plane for the complex coherence corresponding to the multi-polarization channels. Each fitted line is extended and intersected with the unit circle of the complex plane to obtain two points. One of these two points is selected as the surface phase point according to the maximum vegetation deviation principle to obtain the surface phase Ф0. Based on the user-defined vertical wavenumber k... z A two-dimensional lookup table of volume coherence, tree height, and extinction coefficient is established using the surface phase Ф0 and topographic information of the DEM data. Then, the forest tree height is obtained by matching the two-dimensional lookup table of volume coherence, tree height, and extinction coefficient from the improved S-RVoG model. S4. Extract the forest tree height obtained from the inversion, as well as the decomposition parameters, terrain slope, local incident angle, and dominant tree species classification results as input features to construct a volume estimation model. Use the constructed volume estimation model to calculate the forest volume per pixel, generate a spatial distribution map of forest volume in a preset coordinate system, and extract key information from the spatial distribution map of forest volume to generate a total volume statistics table for the corresponding monitoring area according to a preset format.
[0007] Optionally, in S2, the preprocessing method includes: performing image registration on the SAR data using any one or more of the correlation coefficient method, maximum spectrum estimation method, and average fluctuation function method; wherein, two images are randomly acquired and used as the main image M1 and the auxiliary image M2; The correlation coefficient method is set as follows: the complex correlation coefficient between the main image M1 and the auxiliary image M2 is calculated through the search window of the step search algorithm, and the complex correlation coefficient is: ; Where I and J are constants, defining the size of the search window respectively; u is the azimuth offset of the auxiliary image M2 in the step search algorithm, v is the range offset of the auxiliary image M2 in the step search algorithm, and i and j are loop variables, representing the coordinates of any pixel within the search window. * The conjugate operation is represented by adjusting the size of the search window and moving it in the auxiliary image M2 to calculate the complex correlation coefficient of the search window at each feature point in the main image M1. When the complex correlation coefficient reaches a maximum value, the current feature point is determined to be the registration point. The maximum spectrum estimation method is set as follows: the interference image between the main image M1 and the auxiliary image M2 is transformed from the spatial domain to the frequency domain by two-dimensional Fourier transform, and the ratio of the maximum value of the frequency domain signal intensity to the noise value is calculated to obtain the signal-to-noise ratio (SNR) of the interference image between the main image M1 and the auxiliary image M2. ; Where max() represents taking the maximum value, |V(f i ,f j The sign represents the amplitude of the interferometric image in the frequency domain after a two-dimensional Fourier transform; when the signal-to-noise ratio (SNR) of the interferometric image reaches its maximum value, it is determined that the main image M1 and the auxiliary image M2 have been successfully registered. The average fluctuation function method is set as follows: the average phase change value f in the region near the registration point is calculated using the average fluctuation function as a registration evaluation index between the main image M1 and the auxiliary image M2. When the average phase change value f reaches a minimum value, the corresponding point is determined to be a registration point on the auxiliary image M2. The average fluctuation function is: ; Where P(i,j) represents the interference phase value at pixel coordinate (i,j) on the phase map, P(i+1,j) represents the interference phase value of the pixel in the same row (j) and next column (i+1), and P(i,j+1) represents the interference phase value of the pixel in the same column (i) and next row (j+1).
[0008] Optionally, before image registration of the SAR data, the acquired SAR data is uniformly resampled to a preset resolution, and bilinear interpolation is used to ensure that the SAR data corresponds one-to-one with the preset resolution.
[0009] Optionally, in S2, the preprocessing method includes: performing spectral filtering on the SAR data, specifically including: In the interferometric image, a filter window of appropriate size is selected based on the requirements, and the spectral data of the interferometric image is obtained through two-dimensional Fourier transform. ; ; Where u is the spatial frequency in the azimuth direction, v is the spatial frequency in the range direction, and FFT is the two-dimensional Fourier transform; The spectrum data With rectangular window Perform convolution smoothing to obtain the smoothed spectral amplitude. : ; Select the smoothed spectral amplitude The median value was used as the standard for the smoothed spectral amplitude. Normalization is performed: ; Select filter parameters between 0 and 1 As the average correlation of effective pixels in the sliding window, and for the spectral function Perform weighted processing: ; A two-dimensional Fourier transform is performed on the spectrum of the weighted interferometric image to obtain the spatially frequency-domain filtered interferometric image. : .
[0010] Optionally, in S2, the preprocessing method includes: performing interferometric coherence processing on the SAR data, specifically including: Obtain the mapping function after registration of the main image M1 and the auxiliary image M2, and establish the correspondence between the main image M1 and the auxiliary image M2 according to the mapping function; resample the auxiliary image M2, and multiply the registered main image M1 and the auxiliary image M2 by conjugate to obtain the interference image.
[0011] Optionally, in S2, the preprocessing method includes: performing de-flattening phase processing on the SAR data, specifically including: Set the orbital position of the main image M1 during imaging as follows: The orbital position during the imaging of the auxiliary image M2 is The position of any pixel in the main image M1 and the auxiliary image M2 on the reference ellipsoid is: Calculate the baseline parallel component B; ; Where D(M, P) is the distance between M and P, and D(S, P) is the distance between S and P, the flat-ground effect reference phase of the main image M1 and the auxiliary image M2 is: ; Where λ is the radar wavelength; For the entire image, first calculate the phase values of the four vertices and the image center, and then use a two-dimensional polynomial to fit the phase values of the remaining points to obtain the flat-ground effect phase. The expression for the two-dimensional polynomial fitting is: ; Where x is the row coordinate of any point in the image, and y is the column coordinate of any point in the image. d is the coefficient. The order of is [1, 2, 3]; The fitted flat-earth effect phase is converted into a reference phase complex value R: ; Where i is the calculated result of the flat-ground effect phase f(x,y); By performing convolutional multiplication on the main image M1, the auxiliary image M2, and the reference phase complex value, an interferometric image with flat-ground phase eliminated is obtained. : .
[0012] Optionally, in S2, the preprocessing method includes: performing terrain-removing phase processing on the SAR data, specifically including: Determining the topographic phase of each pixel in the interferometric image using a DEM : ; in, It is the radar wavelength. It is the ground elevation; Interference image after removing the flat phase The terrain phase is removed to obtain the terrain-phase-removed interferometric image. : .
[0013] Optionally, in S2, the preprocessing method includes: performing interferometric fringe filtering on the SAR data, specifically including: The original image is transformed from the spatial domain to the frequency domain by performing a Fourier transform, and the components in the image that match the frequency of the filter are retained by using the preset frequency components of the filter, while the remaining components in the image are suppressed, to obtain the filtered image F(jw). The filtered image is transformed back from the frequency domain to the spatial domain by performing an inverse Fourier transform, resulting in the processed interference fringe pattern f(t). .
[0014] Optionally, in S2, the method for generating the polarization coherence matrix T6 specifically includes: Target information is captured by transmitting and receiving signals with different polarization methods to the same scatterer, including HH, HV, VH, and VV polarization methods; and a target polarization scattering matrix S is established based on the captured multiple sets of target information. i : ; The target polarization scattering matrix S i The target vector k is obtained after Pauli polarization decomposition. i ; ; The target vector k i The outer product of the product with its conjugate transpose yields the polarization interference coherence matrix T6: ; Where k1 represents the observation data of the main image, k2 represents the observation data of the auxiliary image, and k6 represents the 6×1 target vector representing one polarimetric interferometric pixel synthesized from k1 and k2. H Represented as matrix conjugate calculation, T 11 T is the Hermitian positive semi-definite matrix of the main image. 12 It is the Hermitian positive semi-definite matrix of the auxiliary image. It is a non-Hermitian complex matrix containing polarization and interference information.
[0015] Optionally, in S3, the improved S-RVoG model is set as follows: When the forest canopy height r in the monitoring area h and vertical height h v When the terrain slope is inconsistent, the polarization coherence matrix T of the RVOG model is respectively... 11 and polarization interference matrix 12 An improved S-RVoG model is obtained by introducing a terrain slope factor for correction; The coherence of the RVOG model is as follows: ; ; ; Among them, I V For volume scattering integral, I D For dihedral scattering integral, I G For the surface scattering integral, T v T d T g These are the scattering matrices corresponding to the integral normalization; The polarization vector, σ is the complex coherence coefficient, θ is the extinction coefficient, θ is the radar incident angle, α is the terrain slope angle, kz′ is the effective vertical wavenumber, z′ is the vertical height coordinate within the vegetation layer, and h v h represents the height of the forest canopy. d Let be the height of the dihedral scattering component, and δ(z′) be the Dirac delta function; The coherence of the improved S-RVoG model is as follows: ; in It is scattering from the Earth's surface. For volume scattering, Even-order scattering phase, The ratio of surface scattered power to volume scattered power. The ratio of dihedral (even-order) scattering power to volume scattering power is given by the improved S-RVoG model, which incorporates a trunk layer and a terrain slope factor. and They are respectively: ; ; The volume scattering complex coherence is expressed as: ; Where, k z h is the vertical wavenumber. i r is the vegetation height corresponding to the effective path of the surface scattering component through the canopy to the ground. h This is the normalized attenuation factor, with a value range of [0, 1].
[0016] The beneficial effects that this invention can produce include: 1. The forest stock volume monitoring method based on a SAR remote sensing quantitative model provided in this invention introduces a slope factor into an improved S-RVoG model to correct the radar incident angle and extinction path, establish the correlation between actual vegetation height and vertical height, and add slope correction parameters to the polarimetric coherence matrix and polarimetric interferometry matrix to distinguish the phase contributions of surface scattering, volume scattering, and even-order scattering, making the model more consistent with the physical scattering mechanism under complex terrain. Simultaneously, it integrates 0.5m high-precision DEM data acquired by UAV LiDAR to effectively suppress interferometric phase distortion caused by terrain undulations, and uses it as an external ground truth to calibrate the inversion results. Furthermore, it avoids nonlinear optimization errors through complex coherence linear fitting (overall least squares method), surface phase extraction (maximum vegetation deviation principle), and two-dimensional lookup table (LUT) matching. Finally, in mountainous areas with an average slope ≥23°, the tree height inversion RMSE decreased from 3.6m to 2.5m, and the coefficient of determination (R²) decreased significantly. 2 The error was reduced from 0.65 to 0.82, with an overall reduction of 30%. Specifically, the error in the 20-30° slope section was reduced to 2.3m (a reduction of 15%), and the error in the 10-20° slope section was reduced to 2.1m (a reduction of 18%), significantly enhancing the adaptability to complex terrain.
[0017] 2. This invention abandons the single backscattering coefficient and introduces multi-dimensional core features (including non-intensity features such as PolInSAR tree height, HV / HH polarization ratio, polarization entropy H, and anisotropy A), and mines the synergistic effect of features through random forest nonlinear regression, in areas with a volume > 200m 3In high-density areas of / ha, the RMSE was reduced to 6.4m using traditional methods. 3 / ha,R 2 The value was increased to 0.75, completely resolving the estimation failure problem caused by signal saturation. Model input features (such as slope, local incidence angle, and polarization decomposition parameters) were automatically exported from remote sensing data such as SAR and DEM, eliminating the need to rely on local diameter-volume tables or numerous ground calibration parameters. Simultaneously, the hyperparameters of the random forest (ntree=500, mtry=4, min_node=5) were optimized through 5-fold spatial hierarchical cross-validation (k-means clustering by latitude and longitude to avoid spatial autocorrelation). The final volumetric volume estimation Rt was thus improved. 2 =0.75, RMSE=6.33m 3 / ha, relative RMSE=12.7%, residuals with no significant systematic bias (DW=1.98), accuracy meets the business standards for forestry carbon sink measurement and logging quota setting.
[0018] 3. This invention ensures a one-to-one correspondence between SAR and DEM pixels by uniformly resampling all data to 10m resolution (bilinear interpolation); it automatically completes coordinate system transformation (e.g., WGS84 to CGCS2000) without manual adjustment, reducing human error; it also integrates correlation coefficient method / maximum spectrum estimation method (image registration with an accuracy of 0.1 pixels), spectral adaptive filtering (dynamically adjusting smoothing intensity according to coherence to balance noise reduction and information preservation), orbit parameter removal of flat phase, and DEM removal of terrain phase, forming a standardized processing script for SAOCOM satellite PolInSAR data; furthermore, it extracts multi-dimensional features in batches at once, allowing the random forest model to fix hyperparameters and avoid repeated grid searches, significantly shortening computation time. (Based on 2000km...) 2 Taking the region (9 scenes of SAOCOM data) as an example, the processing time has been shortened from more than 72 hours in the traditional way to 22 hours, and the processing time for a single scene of data (70km×20km) is less than 60 minutes, which meets the operational timeliness requirements of the annual forestry resource change survey "monthly response". Attached Figure Description
[0019] Figure 1 This is a flowchart of the SAOCOM satellite PolInSAR data preprocessing process of the present invention; Figure 2 This is a schematic diagram of the improved S-RVoG model architecture of the present invention; Figure 3 This is a schematic diagram of the complex circle in the improved S-RVoG model of the present invention; Figure 4 This is a sample plot of the sample data of the present invention. Detailed Implementation
[0020] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0021] Please see Figure 1-3 As shown, this invention provides a method for monitoring forest stock volume based on a SAR remote sensing quantitative model, comprising the following steps: Step S1: Obtain SAR data and ground measurement data for any monitoring area. The ground measurement data includes sample data, topographic information, DEM data, and DSM data. Step S2: Preprocess the SAR data. The preprocessing methods include: image registration, spectral filtering, interferometric coherence, flat phase removal, terrain phase removal, and interferometric fringe filtering in sequence to generate the polarization coherence matrix T6. Step S3: Using an improved S-RVoG model, the radar incident angle and extinction path of the SAR data are corrected by incorporating a terrain slope factor. Based on the polarization coherence matrix T6 and DEM data, forest tree height is inverted using a three-stage inversion method. The three-stage inversion method is set as follows: Based on the overall least squares method, the complex coherence corresponding to the multi-polarization channels is fitted with straight lines on the complex plane. Each fitted line is extended and intersected with the unit circle of the complex plane to obtain two points. One of these two points is selected as the surface phase point according to the principle of maximum vegetation deviation, obtaining the surface phase Ф0. Based on the user-defined vertical wavenumber k... z A two-dimensional lookup table of volume coherence, tree height, and extinction coefficient was established using the surface phase Ф0 of the DEM data and topographic information. Then, forest tree height was obtained from the improved S-RVoG model by matching the two-dimensional lookup table of volume coherence, tree height, and extinction coefficient. The principle of maximum vegetation bias is as follows: in the presence of vegetation, the volume scattering of vegetation will cause the total coherence phase to shift from the true phase point at the surface to another point of the volume scattering phase (i.e., produce a "bias"). On the unit circle of the complex plane, the phase change from the volume scattering phase point to the surface phase point is monotonic. The surface phase is the point with the "maximum bias" from the volume scattering cluster; therefore, among the two intersection points, the point with the larger phase value is selected as the surface phase. Step S4: Extract the inverted forest tree height, as well as polarization decomposition parameters, terrain slope, local incident angle, and dominant tree species classification results as input features to construct a volume estimation model. Use the constructed volume estimation model to calculate the forest volume per pixel, generate a spatial distribution map of forest volume in a preset coordinate system, and extract key information from the spatial distribution map of forest volume to generate a total volume statistics table for the corresponding monitoring area according to a preset format.
[0022] Furthermore, before image registration of the SAR data, the acquired SAR data is uniformly resampled to a preset resolution, and bilinear interpolation is used to ensure that the SAR data corresponds one-to-one with the preset resolution.
[0023] Furthermore, in step S2, the image registration method for SAR data adopts any one or more of the correlation coefficient method, maximum spectrum estimation method, and average fluctuation function method; wherein, two images are randomly acquired and used as the main image M1 and the auxiliary image M2; The correlation coefficient method is set as follows: the complex correlation coefficient between the main image M1 and the auxiliary image M2 is calculated through the search window of the step search algorithm. The complex correlation coefficient is: ; Where I and J are constants, defining the size of the search window respectively; u is the azimuth offset of the auxiliary image M2 in the step search algorithm, v is the range offset of the auxiliary image M2 in the step search algorithm, and i and j are loop variables, representing the coordinates of any pixel within the search window. * The conjugate operation is represented by adjusting the size of the search window and moving it in the auxiliary image M2 to calculate the complex correlation coefficient of the corresponding search window at each feature point in the main image M1. When the complex correlation coefficient reaches a maximum value, the current feature point is determined to be the registration point. The maximum spectrum estimation method is set as follows: the interference image between the main image M1 and the auxiliary image M2 is transformed from the spatial domain to the frequency domain through two-dimensional Fourier transform, and the ratio of the maximum value of the frequency domain signal intensity to the noise value is calculated to obtain the signal-to-noise ratio (SNR) of the interference image between the main image M1 and the auxiliary image M2. ; Where max() represents taking the maximum value, |V(f i ,f j The sign represents the amplitude of the interferometric image in the frequency domain after a two-dimensional Fourier transform; when the signal-to-noise ratio (SNR) of the interferometric image reaches its maximum value, the registration of the main image M1 and the auxiliary image M2 is determined to be successful. The average fluctuation function method is set as follows: the average phase change value f in the region near the registration point is calculated using the average fluctuation function as an evaluation index for registration between the main image M1 and the auxiliary image M2. When the average phase change value f reaches its minimum value, the corresponding point is determined to be a registration point on the auxiliary image M2. The average fluctuation function is: ; Where P(i,j) represents the interference phase value at pixel coordinate (i,j) on the phase map, P(i+1,j) represents the interference phase value of the pixel in the same row (j) and next column (i+1), and P(i,j+1) represents the interference phase value of the pixel in the same column (i) and next row (j+1).
[0024] Furthermore, in step S2, the method for performing spectral filtering on the SAR data is set as follows: In the interferometric image, a filter window of appropriate size is selected based on the requirements, and the spectral data of the interferometric image is obtained through two-dimensional Fourier transform. ; ; Where u is the spatial frequency in the azimuth direction, v is the spatial frequency in the range direction, and FFT is the two-dimensional Fourier transform; Spectrum data With rectangular window Perform convolution smoothing to obtain the smoothed spectral amplitude. : ; Select the smoothed spectral amplitude The median value was used as the standard for the smoothed spectral amplitude. Normalization is performed: ; Select filter parameters between 0 and 1 As the average correlation of effective pixels in the sliding window, and for the spectral function Perform weighted processing: ; A two-dimensional Fourier transform is performed on the spectrum of the weighted interferometric image to obtain the spatially frequency-domain filtered interferometric image. : .
[0025] Furthermore, in step S2, the SAR data is subjected to an interferometric coherence processing method, which specifically includes: Obtain the mapping function after registration of the main image M1 and the auxiliary image M2, and establish the correspondence between the main image M1 and the auxiliary image M2 based on the mapping function; resample the auxiliary image M2, and multiply the registered main image M1 and the auxiliary image M2 by their conjugates to obtain the interference image.
[0026] Furthermore, in step S2, the method for removing the flat-land phase from the SAR data is set as follows: Set the orbital position of the main image M1 during imaging as follows: The orbital position during the imaging of the auxiliary image M2 is The position of any pixel in the main image M1 and the auxiliary image M2 on the reference ellipsoid is: Calculate the baseline parallel component B; ; Where D(M, P) is the distance between M and P, and D(S, P) is the distance between S and P, the flat-ground effect reference phase of the main image M1 and the auxiliary image M2 is: ; Where λ is the radar wavelength; For the entire image, first calculate the phase values of the four vertices and the image center, and then use a two-dimensional polynomial to fit the phase values of the remaining points to obtain the flat-ground effect phase. The expression for the two-dimensional polynomial fitting is: ; Where x is the row coordinate of any point in the image, and y is the column coordinate of any point in the image. d is the coefficient. The order of is [1, 2, 3]; The fitted flat-earth effect phase is converted into a reference phase complex value R: ; Where i is the calculated result of the flat-ground effect phase f(x,y); By performing convolutional multiplication on the main image M1, the auxiliary image M2, and the reference phase complex value, an interferometric image with flat-ground phase eliminated is obtained. : .
[0027] Furthermore, in step S2, the SAR data undergoes a terrain-removing phase processing method, specifically including: Determining the topographic phase of each pixel in the interferometric image using a DEM : ; in, It is the radar wavelength. It is the ground elevation; Interference image after removing the flat phase The terrain phase is removed to obtain the terrain-phase-removed interferometric image. : .
[0028] Furthermore, in step S2, the SAR data is subjected to interferometric fringe filtering processing, which specifically includes: The original image is transformed from the spatial domain to the frequency domain by performing a Fourier transform. The components in the image that match the frequency of the filter are retained by the preset frequency components of the filter, while the other components in the image are suppressed, resulting in the filtered image F(jw). The filtered image is transformed back from the frequency domain to the spatial domain by performing an inverse Fourier transform, resulting in the processed interference fringe pattern f(t). .
[0029] Furthermore, in step S2, the method for generating the polarization coherence matrix T6 specifically includes: Target information is captured by transmitting and receiving signals with different polarization methods to the same scatterer. These polarization methods include HH (horizontal transmission, horizontal reception), HV (horizontal transmission, vertical reception), VH (vertical transmission, horizontal reception), and VV (vertical transmission, vertical reception). A target polarization scattering matrix S is then established based on the captured multiple sets of target information. i : ; The target polarization scattering matrix S i The target vector k is obtained after Pauli polarization decomposition. i ; ; The target vector k i The outer product of the product with its conjugate transpose yields the polarization interference coherence matrix T6: ; Where k1 represents the observation data of the main image, k2 represents the observation data of the auxiliary image, and k6 represents the 6×1 target vector representing one polarimetric interferometric pixel synthesized from k1 and k2. H Represented as matrix conjugate calculation, T 11 T is the Hermitian positive semi-definite matrix of the main image. 12 It is the Hermitian positive semi-definite matrix of the auxiliary image. It is a non-Hermitian complex matrix containing polarization and interference information.
[0030] Furthermore, in step S3, the improved S-RVoG model is set as follows: When the forest canopy height r in the monitoring area h and vertical height h v When the terrain slope is inconsistent, the polarization coherence matrix T of the RVOG model is respectively... 11 and polarization interference matrix 12 An improved S-RVoG model is obtained by introducing a terrain slope factor for correction; The coherence of the RVOG model is: ; ; ; Among them, I V For volume scattering integral, I D For dihedral scattering integral, I G For the surface scattering integral, T v T d T g These are the scattering matrices corresponding to the integral normalization; The polarization vector, σ is the complex coherence coefficient, θ is the extinction coefficient, θ is the radar incident angle, α is the terrain slope angle, kz′ is the effective vertical wavenumber, z′ is the vertical height coordinate within the vegetation layer, and h v h represents the height of the forest canopy. d Let be the height of the dihedral scattering component, and δ(z′) be the Dirac delta function; The coherence of the improved S-RVoG model is: ; in It is scattering from the Earth's surface. For volume scattering, Even-order scattering phase, The ratio of surface scattered power to volume scattered power. The ratio of dihedral (even-order) scattering power to volume scattering power is given by the improved S-RVoG model, which incorporates a trunk layer and a terrain slope factor. and They are respectively: ; ; The volume scattering complex coherence is expressed as: ; Where, k z h is the vertical wavenumber. i r is the vegetation height corresponding to the effective path of the surface scattering component through the canopy to the ground. h The normalized attenuation factor has a value range of [0, 1] and is used to compensate for the difference between the calculated result and the actual vertical scattering distribution.
[0031] In a specific embodiment of the present invention, a certain experimental area was selected as the monitoring area, and forest tree height, diameter at breast height (DBH), tree species, and other parameters were uniformly collected from 39 quadrats to generate 39 sets of sample data. It is worth noting that the selected area for the sample data should represent the typical characteristics of the entire study area, such as vegetation type. The quadrat size for each set of sample data is approximately 1 mu (approximately 667 m²). 2 To ensure sample plots do not span multiple plots, the internal environment of each sample plot should be as uniform as possible. If heterogeneity exists within a sample plot (e.g., significant topographical variations), it is necessary to further subdivide it into sub-sample plots and record data for each sub-plot separately. Simultaneously, sample plots should be located in easily identifiable and revisitable locations for subsequent monitoring and long-term research. For example, sample plots should avoid areas with severe disturbances such as roads and rivers, maintaining a certain distance. When acquiring sample data, starting from the center point of the sample plot, use a measuring tape to measure approximately 18.26m in each of the four directions to establish the four corner points of the quadrat plot. Connect the four corner points to form a square quadrat. Figure 4 As shown.
[0032] Specifically, considering the characteristics of L-band SAR and other satellite data, the following selection principles should be considered when choosing sample plots and sampling points: 1) The distance between the sample plots to be measured is greater than 20 times the spatial resolution of the satellite to be tested; 2) The coefficient of variation within a 5×5 window centered on the measurement point is less than 10%; 3) The quadrat interval is greater than the spatial resolution of the satellite to be tested; 4) The layout of sample points in a single quadrat follows a systematic layout method, which involves dividing the quadrat into small grids with equal spacing and selecting sample points that meet the sampling interval at the center of each small grid.
[0033] 5) Avoid placing all sample points at the edge of the sample plot, otherwise they may be affected by the edge effect. You can leave an appropriate buffer zone.
[0034] 6) The sampling points should be representative of the entire sampling area. If there are significant differences within the area, the sampling points should cover these differences.
[0035] In the above, the tree height, diameter at breast height (DBH), volume, and ground sampling method for tree species are set as follows: S101. Selection of measurement points: Ground parameter collection adopts the average tree method, measuring and recording parameters such as average tree height, average diameter at breast height (DBH), tree species, number of trees per hectare, canopy height, leaf area index, leaf size, branch and trunk dimensions, understory surface roughness, slope, and aspect. The survey is conducted in two categories: pure forest sub-compartments and mixed forest sub-compartments.
[0036] For pure stand sub-compartments: Record the average diameter at breast height (DBH), tree height, canopy height, leaf size, and branch dimensions of 3-5 trees. For sub-compartments with clear row and plant spacing, estimate the number of trees per hectare using the row and plant spacing. First, determine the row spacing, then select representative rows based on tree availability, count at least 5 sample trees, and measure their distance (add half the tree spacing to the ends of the first and last sample trees). Multiply the distance by the row spacing to obtain the area occupied by the sample trees, then calculate the number of trees per hectare using the number of sample trees and the area. For sub-compartments with unclear row and plant spacing, use a combination of visual estimation and actual measurement to survey a 100 m wide area in a representative region. 2 A sample circle (radius 5.64 m) was used to record the number of trees with a diameter at breast height (DBH) of 5 cm or more within the specified area, and the number of trees per hectare was estimated. Forest understory surface roughness was measured using a high-resolution camera to photograph the forest surface, obtained by measuring the understory surface of at least one average tree in both the azimuth and distance directions. Understory surface slope and aspect were obtained using a compass, measured by measuring the understory surface of three average trees.
[0037] For mixed forest sub-compartments: Surveys are conducted separately according to tree species composition. When there are more than three constituent tree species, only the top three species are surveyed. For each species, at least three trees are surveyed for their average diameter at breast height (DBH), tree height, canopy height, leaf size, and branch dimensions. For sub-compartments with significant spacing between trees and rows, the number of trees per hectare is estimated based on the spacing, and recorded separately for each constituent tree species. Understory surface roughness is measured using a surface roughness measurement system, obtained by measuring the understory surface of at least one average tree in both the azimuth and distance directions. Understory surface slope and aspect are obtained using a compass, obtained by measuring the understory surface of three average trees.
[0038] S102. Ground Tree Height Data Acquisition: Tree height is measured using UAV lidar or other measuring tools and recorded to one decimal place. When the tree height is less than 10m, the measurement error is less than 3%; when the tree height is greater than or equal to 10m, the measurement error is less than 5%.
[0039] S103. Diameter at Breast Height (DBH) Data Collection: For tree species with a DBH of 5 cm or more within the quadrat, measure the DBH at a point 1.3 meters above the ground using a DBH measuring tape, recording the value to one decimal place. For trees with irregular DBH, measure and record the DBH values from multiple directions, and take the average. Specifically, for DBH measurements: for trees with a DBH less than 20 cm, the measurement error should be less than 0.3 cm; for trees with a DBH greater than or equal to 20 cm, the measurement error should be less than 1.5%.
[0040] S104. Tree Species Data Collection: Record tree species surveys according to the technical guidelines for professional forest resource surveys. If the tree species cannot be determined, collect leaves or take photos for subsequent identification.
[0041] In the above, the method for estimating the accumulation volume is set as follows: S201. Volume Calculation: Volume models were established for 10 tree species: Pinus tabuliformis, Larix gmelinii, Platycladus orientalis, Quercus, Birch, Populus, Robinia pseudoacacia, Aspen, hardwood, and softwood. The theoretical equations for the volume models are as follows: ; in, This refers to the unit accumulation amount for small classes; Mean diameter at breast height; This represents the average tree height. Let be the average density, i.e., the number of plants per hectare; a, b, and c are parameters. Based on the volumetric volume model, the coefficient R is calculated respectively. 2 The estimated values include the standard deviation (SEE), total relative error (TRE), mean systematic error (MSE), mean prediction error (MPE), and mean percentage standard error (MPSE). It is essential to ensure that the MPE of the accumulation model is below 5% and the MPSE is below 15%.
[0042] S202, Small Class Accumulation Estimation: Based on the Average Diameter of Chest (DBC) of the Small Class Survey ), average tree height ( ), average density ( Using three survey factors and the established volumetric model for the corresponding forest type, the unit volume of each arbor forest sub-compartment is calculated; then, based on the sub-compartment area S, the sub-compartment volume M is obtained. When the sub-compartment is a mixed forest type, the unit volume of each tree species can be calculated using the corresponding volumetric model, and then the total volume of the sub-compartment is obtained. If the total volume of the top three tree species is less than 100%, then the unit volume of the sub-compartment is equal to the sum of the unit volume of the three tree species divided by their respective proportions.
[0043] The methods described above for using drones to collect parameters such as tree height, diameter at breast height (DBH), volume, and tree species include both outdoor and indoor data collection. Specifically, the methods include the following steps: S301. Determine the flight area and sample area: Based on the sample area locations provided by other experimental groups and in conjunction with satellite maps, determine the location of the sample plots to be flown, the flight area, etc.
[0044] S302. Flight route planning and parameter settings: Based on project requirements and parameters such as the camera's field of view, determine the flight altitude and ground resolution (GSD), and design the flight route plan, the number of flight strips, and the spacing between flight strips.
[0045] S303. Equipment Inventory Testing and Safety Check: Inventory the instrument payload and check to confirm that all equipment and accessories are complete; check the battery levels of the drone, payload, and gimbal, and charge them. After the experiment, check all equipment for damage or loss, and charge the equipment.
[0046] S304. UAV Data Acquisition and Processing: Process UAV lidar data to obtain information such as tree height, diameter at breast height, volume, and tree species.
[0047] In this embodiment, the method for constructing the accumulation model is set as follows: S401. Extract the original feature layer (22 dimensions in total): For each isocenter (10m×10m grid window), extract the data as shown in Table 1 simultaneously; Table 1: Raw data for the accumulation model
[0048] It should be noted that all variables were uniformly resampled to 10m, and bilinear interpolation was used to ensure a one-to-one correspondence between SAR and DEM pixels.
[0049] S402. Outlier Removal and Missing Value Handling: Based on the 3σ principle, for continuous variables, samples with an absolute deviation from the mean greater than 3 times the standard deviation are removed; variables with a missing rate > 5% (only "tree species" has a 2% missing rate) are filled with the mode of similar quadrats; finally, the sample size n=37 is retained (originally 39 quadrats, 2 outliers were removed).
[0050] S403, Feature Normalization: Except for "Tree Species", perform Z-factor normalization on the other 21 continuous variables. score standardization: x′=(x-μ x ) / σ x ; The μ and σ parameters are retained for equivalent transformation of new images during operational mapping.
[0051] S404, Multicollinearity Diagnosis: Calculate the variance inflation factor (VIF); when VIF > 10, it is considered redundant; when "VV intensity" and "α angle" VIF = 12.4, delete the VV intensity; retain the remaining 20 dimensions.
[0052] S405, Initial screening of feature importance: Gini importance is obtained based on the initial random forest (ntree=500); variables with importance <0.5% are removed, a total of 3 dimensions (aspect, contrast, correlation) are removed; finally, 17 features are used for formal modeling.
[0053] S406. Training-Validation Sample Construction: A 5-fold spatial hierarchical cross-validation method is used: the samples are clustered into 5 clusters based on their latitude and longitude using k-means. One cluster is selected as the validation set in each round, and the rest are used as the training set. This avoids inflated accuracy due to spatial autocorrelation. The training set is used for grid search to optimize hyperparameters. ntree∈{300,500,700},mtry∈{3,4,5},min node ∈{3,5,7}; With the goal of minimizing RMSE, the following was ultimately determined: ntree=500, mtry=4, min node =5; The validation set is used to report unbiased accuracy: R 2 =0.75, RMSE=6.33m 3 / ha, rRMSE=12.7%; The final results showed that there was no significant systematic bias in residuals (DW=1.98).
[0054] S407. Feature Contribution Interpretation: Top 5 in importance: PolInSAR tree height (27.4%), HV / HH (14.8%), slope (11.2%), entropy H (9.6%), α angle (8.1%); When tree height < 15 m, the accumulation increases linearly with tree height; when slope... At that time, the estimated volume of the sediment was too low, and the model had automatically compensated for the terrain effect.
[0055] S408, Accumulation Calculation—Final Prediction Equation: The processed 17-dimensional feature vector is denoted as: X=[H,σ,μ,H ent , A,α,HV,HH,HV / HH,s,θ,z,R,mean GLCM ,var GLCM ,entr GLCM ,T] ᵀ ; The random forest regression model gives the accumulation amount Y of the i-th sample plot. i (Unit: m) 3 Explicit predictions for / ha): ; Among them, T b For the b-th regression tree, Θ b These are the parameters for the corresponding leaf nodes (included in the training model); Input a 17-dimensional feature vector X into the training model i All data is automatically exported from SAR, DEM, and ground survey data, requiring no empirical coefficients; outputs forest stock volume Y per pixel (10 m × 10 m). i The value range is [0, 450] m 3 / ha; it can be used for: Directly generate a county-scale spatial distribution map of timber volume (GeoTIFF, CGCS2000) for forestry departments to use for carbon sequestration measurement and logging quota setting; The total stock volume is summarized and assigned to the corresponding level of the monitoring area (such as county level or forest farm) to obtain a statistical table of total stock volume, which is used for the annual forest resource change survey. As the input for carbon storage conversion, it is calculated at 0.47 tC / m³. 3 The coefficient is converted to carbon tons.
[0056] S409. Connection with Volume Equation: Although random forests are black boxes, their core learning objective is the macroscopic representation of the univariate timber volume equation. For the same batch of sample plots, the measured volume Y is calculated using the traditional univariate volume table method. ref : ; Where N is the number of plants per hectare, and D j For the jth tree, diameter at breast height (DBH) and height (H) j For the height of the j-th tree, f j For tree species shape number.
[0057] In the above, Y i With Y ref A one-to-one mapping is established on 37 sample plots; therefore, random forest essentially approximates Y using 17-dimensional remote sensing feature regression. ref This solves the technical bottleneck of remote sensing's inability to segment data into individual units without explicitly acquiring D, H, and N at the pixel scale.
Claims
1. A method for monitoring forest stock volume based on a SAR remote sensing quantitative model, characterized in that, Includes the following steps: S1. Acquire SAR data and ground-measured data for any monitoring area, wherein the ground-measured data includes sample data, DEM data, and DSM data; S2. Preprocess the SAR data to generate the polarization coherence matrix T6; S3. An improved S-RVoG model is adopted, and a terrain slope factor is introduced to correct the radar incident angle and extinction path of the SAR data. Based on the polarization coherence matrix T6 and the DEM data, the forest tree height is inverted through a three-stage inversion method. The three-stage inversion method is set as follows: based on the overall least squares method, the complex coherence corresponding to the multipolar channel is fitted with a straight line on the complex plane, and each fitted line is extended and intersected with the unit circle of the complex plane to obtain two points; Based on the principle of maximum vegetation deviation, one of the two points is selected as the surface phase point to obtain the surface phase Ф0; according to the user-defined vertical wavenumber k... z A two-dimensional lookup table of volume coherence, tree height, and extinction coefficient is established using the surface phase Ф0 and topographic information of the DEM data. Then, the forest tree height is obtained by matching the two-dimensional lookup table of volume coherence, tree height, and extinction coefficient from the improved S-RVoG model. S4. Extract the forest tree height obtained from the inversion, as well as the decomposition parameters, terrain slope, local incident angle, and dominant tree species classification results as input features to construct a volume estimation model. Use the constructed volume estimation model to calculate the forest volume per pixel, generate a spatial distribution map of forest volume in a preset coordinate system, and extract key information from the spatial distribution map of forest volume to generate a total volume statistics table for the corresponding monitoring area according to a preset format.
2. The forest stock volume monitoring method based on SAR remote sensing quantitative model according to claim 1, characterized in that, In S2, the preprocessing method includes: performing image registration on the SAR data using any one or more of the correlation coefficient method, maximum spectrum estimation method, and average fluctuation function method; wherein, two images are randomly acquired and used as the main image M1 and the auxiliary image M2; The correlation coefficient method is set as follows: the complex correlation coefficient between the main image M1 and the auxiliary image M2 is calculated through the search window of the step search algorithm, and the complex correlation coefficient is: ； Where I and J are constants, defining the size of the search window respectively; u is the azimuth offset of the auxiliary image M2 in the step search algorithm, v is the range offset of the auxiliary image M2 in the step search algorithm, and i and j are loop variables, representing the coordinates of any pixel within the search window. * The conjugate operation is represented by adjusting the size of the search window and moving it in the auxiliary image M2 to calculate the complex correlation coefficient of the search window at each feature point in the main image M1. When the complex correlation coefficient reaches a maximum value, the current feature point is determined to be the registration point. The maximum spectrum estimation method is set as follows: the interference image between the main image M1 and the auxiliary image M2 is transformed from the spatial domain to the frequency domain by two-dimensional Fourier transform, and the ratio of the maximum value of the frequency domain signal intensity to the noise value is calculated to obtain the signal-to-noise ratio (SNR) of the interference image between the main image M1 and the auxiliary image M2. ； Where max() represents taking the maximum value, |V(f i ,f j The sign represents the amplitude of the interferometric image in the frequency domain after a two-dimensional Fourier transform; when the signal-to-noise ratio (SNR) of the interferometric image reaches its maximum value, it is determined that the main image M1 and the auxiliary image M2 have been successfully registered. The average fluctuation function method is set as follows: the average phase change value f in the region near the registration point is calculated using the average fluctuation function as a registration evaluation index between the main image M1 and the auxiliary image M2. When the average phase change value f reaches a minimum value, the corresponding point is determined to be a registration point on the auxiliary image M2. The average fluctuation function is: ； Where P(i,j) represents the interference phase value at pixel coordinate (i,j) on the phase map, P(i+1,j) represents the interference phase value of the pixel in the same row (j) and next column (i+1), and P(i,j+1) represents the interference phase value of the pixel in the same column (i) and next row (j+1).
3. The forest stock volume monitoring method based on SAR remote sensing quantitative model according to claim 1, characterized in that, Before performing image registration on the SAR data, the acquired SAR data is uniformly resampled to a preset resolution, and bilinear interpolation is used to ensure that the SAR data corresponds one-to-one with the preset resolution.
4. The forest stock volume monitoring method based on SAR remote sensing quantitative model according to claim 1, characterized in that, In S2, the preprocessing method includes: performing spectral filtering on the SAR data, specifically including: In the interferometric image, a filter window of appropriate size is selected based on the requirements, and the spectral data of the interferometric image is obtained through two-dimensional Fourier transform. ; ； Where u is the spatial frequency in the azimuth direction, v is the spatial frequency in the range direction, and FFT is the two-dimensional Fourier transform; The spectrum data With rectangular window Perform convolution smoothing to obtain the smoothed spectral amplitude. : ； Select the smoothed spectral amplitude The median value was used as the standard for the smoothed spectral amplitude. Normalization is performed: ； Select filter parameters between 0 and 1 As the average correlation of effective pixels in the sliding window, and for the spectral function Perform weighted processing: ； A two-dimensional Fourier transform is performed on the spectrum of the weighted interferometric image to obtain the spatially frequency-domain filtered interferometric image. : 。 5. The forest stock volume monitoring method based on SAR remote sensing quantitative model according to claim 1, characterized in that, In S2, the preprocessing method includes: performing interferometric coherence processing on the SAR data, specifically including: Obtain the mapping function after registration of the main image M1 and the auxiliary image M2, and establish the correspondence between the main image M1 and the auxiliary image M2 according to the mapping function; resample the auxiliary image M2, and multiply the registered main image M1 and the auxiliary image M2 by conjugate to obtain the interference image.
6. The forest stock volume monitoring method based on SAR remote sensing quantitative model according to claim 1, characterized in that, In S2, the preprocessing method includes: performing flat-land phase removal processing on the SAR data, specifically including: Set the orbital position of the main image M1 during imaging as follows: The orbital position during the imaging of the auxiliary image M2 is The position of any pixel in the main image M1 and the auxiliary image M2 on the reference ellipsoid is: Calculate the baseline parallel component B; ； Where D(M, P) is the distance between M and P, and D(S, P) is the distance between S and P, the flat-ground effect reference phase of the main image M1 and the auxiliary image M2 is: ； Where λ is the radar wavelength; For the entire image, first calculate the phase values of the four vertices and the image center, and then use a two-dimensional polynomial to fit the phase values of the remaining points to obtain the flat-ground effect phase. The expression for the two-dimensional polynomial fitting is: ； Where x is the row coordinate of any point in the image, and y is the column coordinate of any point in the image. d is the coefficient. The order of is [1, 2, 3]; The fitted flat-earth effect phase is converted into a reference phase complex value R: ； Where i is the calculated result of the flat-ground effect phase f(x,y); By performing convolutional multiplication on the main image M1, the auxiliary image M2, and the reference phase complex value, an interferometric image with flat-ground phase eliminated is obtained. : 。 7. The forest stock volume monitoring method based on SAR remote sensing quantitative model according to claim 1, characterized in that, In S2, the preprocessing method includes: performing terrain-removing phase processing on the SAR data, specifically including: Determining the topographic phase of each pixel in the interferometric image using a DEM : ； in, It is the radar wavelength. It is the ground elevation; Interference image after removing the flat phase The terrain phase is removed to obtain the terrain-phase-removed interferometric image. : 。 8. The forest stock volume monitoring method based on SAR remote sensing quantitative model according to claim 1, characterized in that, In S2, the preprocessing method includes: performing interferometric fringe filtering on the SAR data, specifically including: The original image is transformed from the spatial domain to the frequency domain by performing a Fourier transform, and the components in the image that match the frequency of the filter are retained by using the preset frequency components of the filter, while the remaining components in the image are suppressed, to obtain the filtered image F(jw). The filtered image is transformed back from the frequency domain to the spatial domain by performing an inverse Fourier transform, resulting in the processed interference fringe pattern f(t). 。 9. The forest stock volume monitoring method based on SAR remote sensing quantitative model according to claim 1, characterized in that, In S2, the method for generating the polarization coherence matrix T6 specifically includes: Target information is captured by transmitting and receiving signals with different polarization methods to the same scatterer, including HH, HV, VH, and VV polarization methods; and a target polarization scattering matrix S is established based on the captured multiple sets of target information. i : ； The target polarization scattering matrix S i The target vector k is obtained after Pauli polarization decomposition. i ; ； The target vector k i The outer product of the product with its conjugate transpose yields the polarization interference coherence matrix T6: ； Where k1 represents the observation data of the main image, k2 represents the observation data of the auxiliary image, and k6 represents the 6×1 target vector representing one polarimetric interferometric pixel synthesized from k1 and k2. H Represented as matrix conjugate calculation, T 11 T is the Hermitian positive semi-definite matrix of the main image. 12 It is the Hermitian positive semi-definite matrix of the auxiliary image. It is a non-Hermitian complex matrix containing polarization and interference information.
10. The forest stock volume monitoring method based on SAR remote sensing quantitative model according to claim 1, characterized in that, In S3, the improved S-RVoG model is set as follows: When the forest canopy height r in the monitoring area h and vertical height h v When the terrain slope is inconsistent, the polarization coherence matrix T of the RVOG model is respectively... 11 and polarization interference matrix 12 An improved S-RVoG model is obtained by introducing a terrain slope factor for correction; The coherence of the RVOG model is as follows: ； ； ； Among them, I V For volume scattering integral, I D For dihedral scattering integral, I G For the surface scattering integral, T v T d T g These are the scattering matrices corresponding to the integral normalization; The polarization vector, σ is the complex coherence coefficient, θ is the extinction coefficient, θ is the radar incident angle, α is the terrain slope angle, kz′ is the effective vertical wavenumber, z′ is the vertical height coordinate within the vegetation layer, and h v h represents the height of the forest canopy. d Let be the height of the dihedral scattering component, and δ(z′) be the Dirac delta function; The coherence of the improved S-RVoG model is as follows: ； in It is scattering from the Earth's surface. For volume scattering, Even-order scattering phase, The ratio of surface scattered power to volume scattered power. The ratio of dihedral (even-order) scattering power to volume scattering power is given by the improved S-RVoG model, which incorporates a trunk layer and a terrain slope factor. and They are respectively: ； ； The volume scattering complex coherence is expressed as: ； Where, k z h is the vertical wavenumber. i r is the vegetation height corresponding to the effective path of the surface scattering component through the canopy to the ground. h This is the normalized attenuation factor, with a value range of [0, 1].