Interference tomography SAR forest canopy height estimation optimization method and system

Abstract

The invention provides an interference tomography SAR forest canopy height estimation optimization method and system, and belongs to the technical field of forest canopy height estimation.In the first stage, envelope line fitting is conducted on relative reflectivity signals of a forest vertical structure, and initial envelope lines of a forest canopy and the earth surface are obtained; determining a corresponding relative reflectivity loss threshold value K by referring to the LiDAR height; in the second stage, noise interference at the top and the bottom of the canopy of the forest profile is screened and eliminated through a signal continuity threshold value, and a more reasonable envelope line is obtained; and a third stage of grading the forest canopy density based on K value optimization of canopy density segmentation, giving different forest canopy and earth surface K values to different grades to perform relative reflectivity envelope line extraction, and finally obtaining a more accurate forest height according to height index subtraction corresponding to the envelope lines of the canopy and the earth surface. According to the method, the forest canopy height estimation precision is effectively improved by adopting a three-stage error optimization method.

Description

Technical Field

[0001] This invention relates to the field of forest canopy height estimation technology, and specifically to an optimization method and system for estimating forest canopy height using interferometric tomography SAR. Background Technology

[0002] Forest height is a crucial indicator of forest growth and health, and a key parameter for estimating forest biomass and carbon storage. However, ground-based forest height measurements are costly and challenging due to factors such as transportation, topography, and the complexity of forest structure. In recent years, the development of UAVs and airborne LiDAR technology has provided excellent technical means for accurate measurement of forest tree height; however, the high cost of data acquisition makes large-scale forest height measurements difficult. Synthetic Aperture Radar (SAR) is an active remote sensing technology with all-weather, all-day operation, strong penetration, and a large observation range. It can penetrate most forest branches and leaves to reach the ground, thus obtaining vertical structure information of forests through SAR signals, compensating for the shortcomings of other remote sensing methods.

[0003] Forest parameters are extracted using synthetic aperture radar (SAR). Commonly used forest canopy height inversion models include the random volume on ground (RVoG) model, inversion methods utilizing the phase difference between the canopy and the ground, and TomoSAR. Representative models of the RVOG model include the three-stage RVOG inversion method and the maximum likelihood estimation algorithm. However, this model, based on an interferometric complex coherence model, has excessively high requirements for initial value settings, resulting in poor practical performance. Among inversion methods utilizing the phase difference between the canopy and the ground, the classic ESPRIT method obtains the phase center of the canopy and the ground and uses the difference between these phase centers to invert canopy height. However, due to the complexity of forest structure, errors exist in the phase center between the canopy and the ground surface. TomoSAR has been widely used for acquiring forest vertical structure. TomoSAR technology obtains three-dimensional structure by separating the reflected signals of target objects along the elevation direction, showing significant advantages in acquiring forest canopy height. In 2000, Reigber et al. first used L-band airborne datasets to obtain the three-dimensional structure of forests through Fast Fourier Transform (FFT) technology, achieving the first successful experiment in forest applications. This laid the foundation for many researchers to apply tomographic SAR to forests.

[0004] Currently, after acquiring forest vertical structure information using TomoSAR, the envelope method is commonly used to extract the envelope of the relative reflectivity signal of the forest vertical structure, thereby obtaining the forest canopy height. However, to implement TomoSAR technology, SAR data preprocessing, temporal decoherence, image reconstruction, and result analysis are required. Data preprocessing includes flat-ground phase removal, terrain phase correction, and phase calibration of SAR data. Due to differences in terrain, target features, and sensors, even after relevant preprocessing, some residual phase effects will remain. Temporal decoherence needs to consider the baseline time interval. Beamforming methods include non-parametric methods (Beamforming, Capon, etc.) and parametric methods (Music, WSF, etc.). The difference between parametric and non-parametric methods is that non-parametric methods do not require prior knowledge of the number of scatterers or scattering mechanisms, and can directly process the signal. In practical applications, due to various influencing factors such as terrain, baseline time interval, vegetation complexity, and uncertainties in beamforming methods, SAR data suffers from phase focusing errors and residual phase effects, causing the envelope to fall outside the vertical structure information of the forest. Furthermore, the penetration of microwaves through different forest covers amplifies the vertical information range of the forest, leading to an overestimation of forest canopy height.

[0005] In summary, after performing terrain phase correction and phase compensation on the data, some residual phase effects still exist due to various influencing factors such as terrain, baseline time interval, vegetation complexity, and uncertainties in beamforming methods. Therefore, how to further improve the accuracy of TomoSAR forest height estimation is a technical problem that urgently needs to be solved. Summary of the Invention

[0006] The purpose of this invention is to provide an optimization method and system for estimating forest canopy height using interferometric tomography SAR, so as to improve the accuracy of forest height inversion by at least reducing the relative reflection noise interference of TomoSAR.

[0007] To achieve the above objectives, this invention provides an optimization method for estimating forest canopy height using interferometric tomography SAR, the method comprising:

[0008] L-band UAVSAR multi-baseline data were read, and a multi-baseline InSAR covariance matrix was constructed based on the ground phase compensation results. The relative reflectance of the forest vertical profile was reconstructed using Capon spectral estimation.

[0009] Envelope fitting was performed on the relative reflectance signal of the forest vertical profile, and the initial envelope of the forest canopy, the initial envelope of the ground surface, and the signal loss threshold K were determined by combining LiDAR height reference data.

[0010] The relative reflectance signal of the forest vertical profile is subjected to noise filtering based on the continuity of the vertical direction signal to remove noise interference from the top of the canopy and the bottom of the ground, so as to obtain the denoised forest vertical structure reflectance profile and a more reasonable envelope.

[0011] Pauli decomposition was performed on the fully polarimetric SAR data to obtain the canopy density index. Different forest canopy K-values ​​and surface K-values ​​were assigned to different levels of pixels according to the canopy density classification results.

[0012] The relative reflectivity signal is re-extracted from the envelope, and the forest canopy height is obtained by subtracting the height index corresponding to the canopy envelope and the ground surface envelope.

[0013] Optionally, the relative reflectance of the reconstructed forest vertical profile is estimated using Capon spectral estimation, including:

[0014] Preprocessing is performed on L-band UAVSAR multi-baseline data;

[0015] The ground phase is estimated by calling the three-stage multi-baseline RVOG method, and residual phase compensation is performed on the multi-baseline observation data based on the ground phase.

[0016] Calculate the multi-baseline InSAR covariance matrix R based on the compensated multi-baseline observation data;

[0017] The multi-baseline InSAR covariance matrix R is input into the Capon beamforming power estimator to calculate the relative reflectance in the vertical direction of the vegetated area. The specific expression is as follows:

[0018] ;

[0019] In the formula, Let z represent the vertical distribution function of backscattered power estimated using the Capon algorithm; z represents the height variable in the vertical direction; a(z) represents the steering vector at height z. R represents the conjugate transpose of the steering vector a(z); R represents the covariance matrix of the multi-baseline InSAR data. Let represent the inverse of the covariance matrix R.

[0020] Optionally, the initial envelope of the forest canopy, the initial envelope of the land surface, and the signal loss threshold K are determined, including:

[0021] Based on the maximum forest canopy height in the study area Apply upper limit constraints to the inversion results;

[0022] Set threshold K of different time lengths to extract the upper and lower envelopes of the TomoSAR relative reflectivity signal;

[0023] Obtain the estimated values of forest canopy height under different threshold conditions based on the difference between the upper and lower envelopes, and exclude the samples with the retrieved height greater than ;

[0024] Select the optimal threshold K according to the RMSE between the retrieved result and LiDAR - RH100. The specific expression is:

[0025] ;

[0026] where, H represents the estimated value of forest canopy height obtained under the current threshold condition; represents the operation of obtaining the value that makes the target quantity minimum; represents the upper envelope height or canopy side height quantity extracted from the envelope under the current threshold condition; represents the lower envelope height or ground surface side height quantity extracted from the envelope under the current threshold condition; K represents the critical value of signal loss; represents the maximum forest canopy height in the study area.

[0027] Optionally, perform noise screening on the relative reflectance signal of the forest vertical profile based on the signal continuity in the vertical direction, including:

[0028] Set the threshold T and use the threshold T as the maximum critical value for the forest canopy to transition upward to noise and the ground surface to transition downward to noise;

[0029] Compare the relative reflectance function F(H) layer by layer from the top to the bottom of the profile according to the height index H. When F(H) < T, record the corresponding height index as the feedback value H k ; discard the height layers where F(H) ≥ T;

[0030] After completing the screening of all height layers, obtain the feedback sequence {H k} arranged in height order, and calculate the adjacent feedback height difference Δh k = H k+1 −H k ;

[0031] Select the two adjacent feedback heights (H k , H k\* , H k\*+1 ) that produce the maximum Δh

[0032] as the main signal boundaries, and extract the relative reflectance corresponding to the height index between these main signal boundaries as the normal distribution range of the denoised forest.

[0033] Optionally, Pauli decomposition can be performed on the fully polarimetric SAR data to obtain canopy density indices, including:

[0034] Perform Pauli decomposition on fully polarimetric SAR data;

[0035] Pauli decomposition brightness was extracted as an indicator of forest canopy density.

[0036] Based on the canopy density index, the study area was divided into sparse cover, medium cover, and dense cover levels;

[0037] Output the canopy density level label for each pixel.

[0038] Optionally, different forest canopy K-values ​​and surface K-values ​​can be assigned to pixels of different levels according to the canopy density grading results, including:

[0039] Establish threshold sets for canopy envelope extraction and land surface envelope extraction for sparse cover, medium cover, and dense cover levels, respectively;

[0040] The corresponding canopy K value and surface K value are retrieved based on the canopy density level to which each pixel belongs;

[0041] The canopy envelope and the surface envelope were re-extracted using the obtained canopy K-value and surface K-value, respectively.

[0042] Adaptive correction is performed on the envelope position offset, and the corrected canopy envelope height index and surface envelope height index are output.

[0043] Optionally, the corresponding canopy K-value and ground K-value can be retrieved based on the canopy density level of each pixel. The specific expression is as follows:

[0044] ;

[0045] H represents the height of the forest canopy; This represents the operation that minimizes the target quantity; P(·) represents the relative reflectance or the evaluation function corresponding to the relative reflectance. This represents the height difference between the canopy side and the ground side at location (r,x); K(C) represents the height of the surface envelope at location (r,x); C(C) represents the critical value corresponding to the coverage level C; C represents the forest coverage level; C(max) represents the dense coverage level; C(medium) represents the medium coverage level; C(min) represents the sparse coverage level; r,x represent the pixel spatial location index; K=0.1,0.2,0.3,⋯,0.8 indicates that different critical values ​​are searched.

[0046] Optionally, the forest canopy height can be obtained by subtracting the height indices corresponding to the canopy envelope and the ground surface envelope, specifically:

[0047] Read the corrected height index corresponding to the canopy envelope and the height index corresponding to the land surface envelope;

[0048] Perform a difference calculation between the height index corresponding to the canopy envelope and the height index corresponding to the land surface envelope to obtain the pixel-by-pixel forest canopy height. The specific expression is as follows:

[0049] ;

[0050] In the formula, This represents the height difference between the canopy side and the ground side at location (r,x); This represents the height corresponding to the surface envelope at location (r,x);

[0051] Spatial aggregation was performed on all pixel-by-pixel forest canopy heights to output the forest canopy height results for the study area.

[0052] To achieve the above objectives, the present invention also provides an optimization system for estimating forest canopy height using interferometric tomography SAR, the system comprising:

[0053] The relative reflectance calculation module is used to read L-band UAVSAR multi-baseline data, construct the multi-baseline InSAR covariance matrix based on the ground phase compensation results, and reconstruct the relative reflectance of the forest vertical profile using Capon spectrum estimation.

[0054] The relative reflectance loss threshold K determination module is used to perform envelope fitting on the relative reflectance signal of the forest vertical profile, and combine LiDAR height reference data to determine the initial envelope of the forest canopy, the initial envelope of the ground surface, and the signal loss threshold K;

[0055] The noise error elimination module is used to perform noise filtering based on the vertical signal continuity on the relative reflectance signal of the forest vertical profile, remove noise interference from the top of the canopy and the bottom of the ground, and obtain a denoised forest vertical structure reflectance profile and a more reasonable envelope.

[0056] The forest canopy density classification module is used to perform Pauli decomposition on fully polarimetric SAR data to obtain canopy density indices, and assign different forest canopy K values ​​and surface K values ​​to pixels of different levels according to the canopy density classification results.

[0057] The forest height inversion module is used to re-extract the envelope from the relative reflectivity signal and obtain the forest canopy height by subtracting the height index corresponding to the canopy envelope and the ground surface envelope.

[0058] Beneficial effects: Through the above technical solution, the present invention proposes an optimized method and system for estimating the forest canopy height by interferometric层析SAR. First, the multi-baseline RVoG three-stage method is used to estimate the ground phase and remove part of the residual phase, and a multi-baseline InSAR covariance matrix is constructed; on this basis, Capon spectral estimation (beamforming power estimator) is used to reconstruct the relative reflectivity profile in the vertical direction of the vegetation area. Secondly, the envelope fitting is performed on the relative reflectivity signal of the vertical profile. By traversing the signal loss critical value K, the initial upper and lower envelopes of the canopy and the ground surface are extracted, and the canopy height is estimated by the difference between the two envelopes; taking the LiDAR-CHM as the benchmark, the optimal K is determined when the RMSE of the estimated height is the smallest, and the initial envelopes of the canopy and the ground surface are obtained accordingly. Subsequently, aiming at the noise residues above the canopy top and below the ground surface caused by the phase focusing error, a threshold screening denoising based on signal continuity is proposed: set the threshold T as the maximum critical value for the transition from the canopy upward and the ground surface downward to the noise, compare layer by layer along the height index and record the feedback height sequence {H k} that satisfies F(H)<T, and calculate the adjacent interval Δh k =H k+1 −H k . The main signal boundary is determined by the adjacent feedback heights corresponding to the maximum Δh, and the relative reflectivity between them is retained to obtain the true profile and the stable envelope. Finally, in order to suppress the envelope system offset caused by the microwave penetration difference under different canopy density conditions, based on the Pauli decomposition brightness as an index for proxy of the canopy density, the study area is divided into three categories: sparse, medium and dense, and the corresponding segmented K is selected to adaptively correct the envelope, so as to improve the accuracy and spatial consistency of the forest canopy height estimation. The present invention improves the accuracy of forest height estimation by correcting the error of the interferometric层析SAR relative reflectivity envelope, and can effectively alleviate the influence of the phase focusing error on the forest height estimation.

[0059] Other features and advantages of the embodiments of the present invention will be described in detail in the subsequent specific embodiments section. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] The drawings are used to provide a further understanding of the embodiments of the present invention, and constitute a part of the specification, and are used together with the following specific embodiments to explain the embodiments of the present invention, but do not constitute a limitation to the embodiments of the present invention. In the drawings:

[0061] Figure 1 is a schematic flow chart of the optimized method for estimating the forest canopy height by interferometric层析SAR of the present invention.

[0062] Figure 2 is a schematic diagram of the envelope fitting of the present invention and the determination principle of the corresponding signal loss threshold K (a is the envelope fitting, and b is the schematic diagram of the determination of the reflectivity loss threshold).

[0063] Figure 3 It is a schematic diagram of noise rejection based on the vertical signal continuity of the present invention (a is the original envelope; b is the schematic diagram of noise rejection for vertical continuity inspection).

[0064] Figure 4 It is a noise rejection result diagram of coverage classification of the present invention (a is the envelope before optimization; b is the envelope after optimization).

[0065] Figure 5 It is a diagram of the change in canopy height before and after optimization of the present invention (a and b are the inversion result diagram and scatter plot before optimization; c is the inversion result diagram after optimization). Specific implementation manners

[0066] The following will describe in detail the specific implementation manners of the present invention with reference to the accompanying drawings. It should be understood that the specific implementation manners described herein are only for the purpose of illustration and explanation of the present invention, and are not used to limit the present invention. [[ID=?]]

[0067] As Figure 1 shown, the implementation manner of the present invention provides an optimized method for estimating the forest canopy height by interferometric tomography SAR, and the method includes:

[0068] S1: Based on the L-band UAVSAR data, the Capon spectral estimation method is used to reconstruct the relative reflectivity of the forest vertical profile, reconstruct the forest vertical structure characteristics, and obtain the vegetation vertical structure profile. This method is based on the traditional tomography method of minimum variance, and realizes high-resolution imaging of the target ground object by optimizing the weight distribution of the radar echo data.

[0069] S2: Fit the envelope of the relative reflectivity signal in the vertical direction of the forest to obtain the initial envelopes of the forest canopy and the ground surface and the corresponding signal loss critical value K; first, determine the maximum forest canopy height (h max ) of the study area according to the reference LiDAR, estimate the upper and lower envelope heights of TomoSAR by setting different step K values to invert the canopy height under different conditions, and剔除 the samples with the inverted height greater than h max and calculate the RMSE between the inversion result and LiDAR-RH100. When the RMSE is the smallest, determine the K value.

[0070] S3: Use the vertical signal continuity to screen and剔除 the noise interference at the top and bottom of the forest profile canopy to obtain a more accurate envelope; first, set the initial threshold T as the maximum critical value for the forest canopy to transition upward to noise and the ground surface to transition downward to noise; then compare layer by layer from the top to the bottom of the profile according to the height index H. When the relative reflectivity F(H) corresponding to a certain height < T, record the corresponding forest height index as, and the feedback value H k, when F(H) ≥ T, it is not recorded until the screening of all height layers of the pixel is completed, and a feedback sequence {H k} arranged in height order is obtained. Based on this, the adjacent feedback height difference Δh k = H k+1 −H k is calculated. Since the span of the effective signal of the forest vertical structure is significantly greater than the span of the top / bottom noise, the two adjacent feedback heights (H k\* , H k\*+1) that produce the maximum Δh are selected as the main signal boundaries, and the relative reflectance corresponding to the height index between them is extracted as the normal distribution range of the denoised forest, so as to effectively remove the noise above the canopy and below the ground surface, and obtain a stable forest vertical structure reflectance profile for subsequent height inversion.

[0071] S4: Combining the LiDAR canopy height, select the K value based on the canopy density segmentation. Classify the forest canopy density, and different canopy densities are given different forest canopy K values and ground surface K values to extract the envelope of the relative reflectance signal in the vertical direction of the forest. Finally, the more accurate forest height is obtained by subtracting the height indices corresponding to the envelopes of the canopy and the ground surface.

[0072] In a preferred embodiment, in step S1, based on the L-band UAVSAR data, the Capon spectral estimation method is used to reconstruct the relative reflectance in the vertical direction of the forest. This method is based on the traditional tomography method of minimum variance, and by optimizing the weight distribution of the radar echo data, high-resolution imaging of the target ground object is achieved.

[0073] In a preferred embodiment, in step S2, the envelope of the relative reflectance signal of the forest vertical profile is fitted to obtain the initial envelopes of the forest canopy and the ground surface and the corresponding signal loss critical value K; first, the maximum forest canopy height (h max ) of the study area is determined according to the following formula, and different step K values are set to estimate the canopy height under different conditions of the upper and lower envelope height inversion of TomoSAR, and the samples with the inversion height greater than h max are excluded and the RMSE between the inversion result and LiDAR-RH100 is calculated. The K value is determined when the RMSE is the smallest.

[0074] In a preferred embodiment, in step S3, the signal continuity threshold is used to screen and remove the noise interference at the top and bottom of the forest profile canopy to obtain a more accurate envelope; first, a threshold T is set as the maximum critical value for the forest canopy to transition upward to the noise and the ground surface to transition downward to the noise; then, the height index H is compared layer by layer from the top to the bottom of the profile. When F(H) < T, the corresponding height index is recorded as the feedback value H kWhen F(H)≥T, no record is made until all height layers of that pixel have been filtered, resulting in a feedback sequence {H} arranged in height order. k Based on this, the adjacent feedback height difference Δh is calculated. k =H k+1 -H k Since the effective signal span of the forest's vertical structure is significantly larger than the span of the top / bottom noise, the two adjacent feedback heights (H) that generate the maximum Δh are selected. k\* H k\*+1) As the boundary of the main signal, the relative reflectance corresponding to the height index between them is extracted as the normal distribution range of the forest after denoising, thereby effectively removing noise above the canopy and below the ground surface, and obtaining a stable forest vertical structure reflectance profile for subsequent height inversion.

[0075] In a preferred embodiment, in step S4, K-values ​​are selected based on canopy density segments, in conjunction with LiDAR canopy height. Forest canopy density is graded, and different canopy densities are assigned different forest canopy K-values ​​and ground surface K-values. The envelope of the relative reflectance signal in the vertical direction of the forest is extracted, and finally, a more accurate forest height is obtained by subtracting the height indices corresponding to the canopy and ground surface envelopes.

[0076] Therefore, this embodiment proposes an optimization method for estimating forest canopy height using interferometric tomography SAR. By correcting the relative reflectivity envelope error of interferometric tomography SAR, the accuracy of forest height estimation is improved, which can effectively mitigate the impact of phase focusing error on forest height estimation.

[0077] To more clearly explain this application, a specific example of the optimization method for estimating forest canopy height using interferometric tomography SAR is provided below. For example... Figure 2As shown in the figure, this embodiment constructs a forest height inversion method and system of "vertical profile reconstruction - envelope extraction - denoising correction - hierarchical threshold adaptation" based on L-band UAVSAR multi-baseline full-polarization data and LiDAR-CHM constraints. First, the ground phase is estimated by using the multi-baseline RVoG three-stage method and part of the residual phase is removed to construct a multi-baseline InSAR covariance matrix; on this basis, Capon spectral estimation (beamforming power estimator) is used to reconstruct the relative reflectivity profile in the vertical direction of the vegetation area. Secondly, the envelope of the relative reflectivity signal of the vertical profile is fitted, and the initial upper and lower envelopes of the canopy and the ground are extracted by traversing the signal loss critical value K, and the canopy height is estimated by the difference between the two envelopes; taking LiDAR-CHM as the reference, the optimal K is determined when the RMSE of the estimated height is the smallest, and the initial envelopes of the canopy and the ground are obtained accordingly. Subsequently, aiming at the noise residues above the canopy top and below the ground caused by the phase focusing error, threshold screening denoising based on signal continuity is proposed: set the threshold T as the maximum critical value for the canopy to transition upward and the ground to transition downward to noise, compare layer by layer along the height index and record the feedback height sequence {H k} that satisfies F(H)<T, calculate the adjacent interval Δh k =H k+1 −H k , and determine the main signal boundary based on the adjacent feedback height corresponding to the maximum Δh, and retain the relative reflectivity between them to obtain the true profile and stable envelope. Finally, in order to suppress the envelope system offset caused by the microwave penetration difference under different canopy density conditions, based on the Pauli decomposition brightness as a proxy index for canopy density, the study area is divided into three categories: sparse, medium and dense, and the corresponding segmented K is selected to adaptively correct the envelope, so as to improve the estimation accuracy and spatial consistency of forest canopy height. The present invention improves the forest biomass estimation accuracy by correcting the relative reflectivity envelope error of interferometric tomography SAR, and can effectively alleviate the influence of phase focusing error on forest height estimation.

[0078] Regarding step S1: Based on the L-band UAVSAR data, the Capon spectral estimation method is used to reconstruct the relative reflectivity of the forest vertical profile, reconstruct the forest vertical structure characteristics, and obtain the vegetation vertical structure profile. This method is based on the traditional tomography method of minimum variance, and realizes high-resolution imaging of the target ground object by optimizing the weight distribution of radar echo data. The specific implementation steps are as follows:

[0079] (1) Data processing:

[0080] The L-band airborne multi-baseline PolInSAR data and lidar validation data were both derived from publicly available datasets from the AfriSAR project. The PolInSAR datasets underwent polarization calibration, baseline fine registration, and spectral filtering, and were provided in single-look complex format, with each orbit containing SLC data for four polarization channels. The test area had eight orbits; the data underwent multi-look processing, and the RVOG three-stage method was used to estimate the ground phase to remove some phase errors and improve phase focusing. The multi-baseline InSAR covariance matrix was also calculated.

[0081] (2) Capon beamforming power estimator:

[0082] The Capon beamforming power estimator is used to estimate the relative reflectance of the vegetation area in the vertical direction. The Capon spectral estimator is a common nonparametric method in tomographic analysis, capable of obtaining an infinite vertical profile of vegetation without any prior knowledge of the statistical properties of the data. This method is an improved algorithm based on the minimum variance criterion of conventional beamforming. It uses optimal weighting vectors to perform spatial filtering on the signals of each array element to suppress noise interference and enhance the desired signal. The spectral estimation formula is as follows:

[0083] ;

[0084] In the formula, The vertical distribution function of backscattered power is estimated using the Capon algorithm. Let z represent the steering vector at height z, and R represent the covariance matrix of the multi-baseline InSAR data.

[0085] Regarding step S2: The relative reflectance signal of the forest vertical profile is fitted with an envelope to obtain the initial envelopes of the forest canopy and the ground surface, as well as the corresponding signal loss critical value K; for example... Figure 2 As shown, Figure 2 A schematic diagram for determining the envelope fitting and reflectivity loss threshold ( Figure 2 In the diagram (a) shows the envelope fitting and (b) shows the determination of the reflectance loss threshold, the maximum forest canopy height (h) in the study area is first determined according to the following formula. max By setting different time length K values ​​to estimate the upper and lower envelope heights of TomoSAR, canopy heights under different conditions are inverted, and canopies with inverted heights greater than h are removed. max The sample is used to calculate the RMSE between the inversion result and LiDAR-RH100. The K value is determined when the RMSE is minimized. The specific implementation steps are as follows:

[0086] ;

[0087] Regarding step S3: Using the signal continuity threshold to screen and eliminate the noise interference at the top and bottom of the forest profile canopy to obtain a more accurate envelope line, as Figures 3-4 shown, Figure 3 is a schematic diagram of noise elimination based on signal continuity in the vertical direction ( Figure 3 in, (a) is the original envelope line, (b) is the schematic diagram of noise elimination for vertical direction continuity inspection), Figure 4 is a schematic diagram of the change of the envelope line before and after optimization, including the following implementation steps:

[0088] First, set the threshold T (that is, K determined in S2) as the maximum critical value for the forest canopy to transition upward to noise and the ground surface to transition downward to noise; then compare layer by layer from the top to the bottom of the profile according to the height index H. When F(H) < T, record the corresponding height index as the feedback value H k , when F(H) ≥ T, do not record until the screening of all height layers of the pixel is completed, and obtain the feedback sequence {H k} arranged in height order. On this basis, calculate the adjacent feedback height difference Δh k =H k+1 −H k , because the effective signal span of the forest vertical structure is significantly greater than the span of the top / bottom noise, so select the two adjacent feedback heights (H k\* ,H k\*+1) that produce the maximum Δh as the main signal boundary, and extract the relative reflectance corresponding to the height index between them as the denoised normal distribution range of the forest, so as to effectively eliminate the noise above the canopy and below the ground surface, and obtain a stable forest vertical structure reflectance profile for subsequent height inversion.

[0089] Regarding step S4: Optimization selection of the K value based on canopy density segmentation, specifically including the following implementation steps:

[0090] Combined with the LiDAR canopy height, selection of the K value based on canopy density segmentation. First, perform Pauli decomposition on the full-polarization SAR data, use the brightness of the Pauli decomposition as an index of canopy density, divide the study area into three coverage levels of sparse, medium and dense, and use the corresponding critical value K respectively in the extraction of the canopy and ground surface envelope lines under different level conditions to achieve adaptive correction of the envelope line position offset, and calculate the interval between the upper and lower envelope lines to obtain the forest height;

[0091] ;

[0092] In the formula, is the forest coverage level, is the critical value of the forest canopy and the ground surface under different forest coverages. After determining the critical value K, calculate the height H through and The RMSE between different coverage levels is used to determine the forest canopy height for each pixel. When the overall RMSE value for different coverage levels is minimized, the forest canopy height can be obtained. Figure 5 As shown, Figure 5 The results show the forest height predictions before and after optimization.

[0093] In a preferred embodiment, an optimization system for estimating forest canopy height using interferometric tomography SAR is also proposed, the system comprising:

[0094] The relative reflectance calculation module is used to read L-band UAVSAR multi-baseline data, construct the multi-baseline InSAR covariance matrix based on the ground phase compensation results, and reconstruct the relative reflectance of the forest vertical profile using Capon spectrum estimation.

[0095] The relative reflectance loss threshold K determination module is used to perform envelope fitting on the relative reflectance signal of the forest vertical profile, and combine LiDAR height reference data to determine the initial envelope of the forest canopy, the initial envelope of the ground surface, and the signal loss threshold K;

[0096] The noise error elimination module is used to perform noise filtering based on the vertical signal continuity on the relative reflectance signal of the forest vertical profile, remove noise interference from the top of the canopy and the bottom of the ground, and obtain a denoised forest vertical structure reflectance profile and a more reasonable envelope.

[0097] The forest canopy density classification module is used to perform Pauli decomposition on fully polarimetric SAR data to obtain canopy density indices, and assign different forest canopy K values ​​and surface K values ​​to pixels of different levels according to the canopy density classification results.

[0098] The forest height inversion module is used to re-extract the envelope from the relative reflectivity signal and obtain the forest canopy height by subtracting the height index corresponding to the canopy envelope and the ground surface envelope.

[0099] Those skilled in the art will understand that all or part of the steps in the methods of the above embodiments can be implemented by a program instructing related hardware. This program is stored in a storage medium and includes several instructions to cause a microcontroller, chip, or processor to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as a USB flash drive, a portable hard drive, a read-only memory (ROM), a random access memory (RAM), a magnetic disk, or an optical disk.

[0100] The optional embodiments of the present invention have been described in detail above with reference to the accompanying drawings. However, the embodiments of the present invention are not limited to the specific details described above. Within the scope of the technical concept of the embodiments of the present invention, various simple modifications can be made to the technical solutions of the embodiments of the present invention, and these simple modifications all fall within the protection scope of the embodiments of the present invention. It should also be noted that the various specific technical features described in the above specific embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the embodiments of the present invention will not further describe the various possible combinations.

[0101] Furthermore, various different embodiments of the present invention can be combined in any way, as long as they do not violate the spirit of the embodiments of the present invention, they should also be regarded as the content disclosed by the embodiments of the present invention.

Claims

1. An optimization method for estimating forest canopy height using interferometric tomography SAR, characterized in that, The method includes: L-band UAVSAR multi-baseline data were read, and a multi-baseline InSAR covariance matrix was constructed based on the ground phase compensation results. The relative reflectance of the forest vertical profile was reconstructed using Capon spectral estimation. Envelope fitting was performed on the relative reflectance signal of the forest vertical profile, and the initial envelope of the forest canopy, the initial envelope of the ground surface, and the signal loss threshold K were determined by combining LiDAR height reference data. The relative reflectance signal of the forest vertical profile is subjected to noise filtering based on the continuity of the vertical direction signal to remove noise interference from the top of the canopy and the bottom of the ground, so as to obtain the denoised forest vertical structure reflectance profile and a more reasonable envelope. Pauli decomposition was performed on the fully polarimetric SAR data to obtain the canopy density index. Different forest canopy K-values ​​and surface K-values ​​were assigned to different levels of pixels according to the canopy density classification results. The relative reflectivity signal is re-extracted from the envelope, and the forest canopy height is obtained by subtracting the height index corresponding to the canopy envelope and the ground surface envelope.

2. The method for optimizing forest canopy height estimation using interferometric tomography SAR according to claim 1, characterized in that, The relative reflectance of the reconstructed forest vertical profile was estimated using Capon spectral estimation, including: Preprocessing is performed on L-band UAVSAR multi-baseline data; The ground phase is estimated by calling the three-stage multi-baseline RVOG method, and residual phase compensation is performed on the multi-baseline observation data based on the ground phase. Calculate the multi-baseline InSAR covariance matrix R based on the compensated multi-baseline observation data; The multi-baseline InSAR covariance matrix R is input into the Capon beamforming power estimator to calculate the relative reflectance in the vertical direction of the vegetated area. The specific expression is as follows: ； In the formula, Let z represent the vertical distribution function of backscattered power estimated using the Capon algorithm; z represents the height variable in the vertical direction; a(z) represents the steering vector at height z. R represents the conjugate transpose of the steering vector a(z); R represents the covariance matrix of the multi-baseline InSAR data. Let represent the inverse of the covariance matrix R.

3. The optimization method for estimating forest canopy height using interferometric tomography SAR according to claim 1, characterized in that, Determine the initial envelope of the forest canopy, the initial envelope of the land surface, and the signal loss threshold K, including: Based on the maximum forest canopy height in the study area Apply upper limit constraints to the inversion results; Set threshold K of different time lengths to extract the upper and lower envelopes of the TomoSAR relative reflectivity signal; Based on the difference between the upper and lower envelopes, forest canopy height estimates under different threshold conditions were obtained, and inversion heights greater than a certain threshold were removed. The sample; The optimal threshold K is selected based on the RMSE between the inversion results and LiDAR-RH100, and the specific expression is as follows: ； Where H represents the forest canopy height estimate obtained under the current threshold conditions; This indicates the operation that minimizes the target value; This represents the upper envelope height or canopy side height obtained from the envelope under the current threshold conditions; This represents the lower envelope height or ground-side height obtained from the envelope under the current threshold condition; K represents the signal loss threshold. This indicates the maximum forest canopy height in the study area.

4. The optimization method for estimating forest canopy height using interferometric tomography SAR according to claim 3, characterized in that, Perform noise filtering based on vertical signal continuity on the relative reflectance signal of the forest vertical profile, including: Set a threshold T, and use the threshold T as the maximum critical value for the upward transition of the forest canopy to noise and the downward transition of the ground surface to noise; Compare the relative reflectance function F(H) layer by layer from the top to the bottom of the cross-section according to the height index H. When F(H) < T is satisfied, record the corresponding height index as the feedback value H k ; Discard the height layer where F(H) ≥ T; After completing the filtering of all height layers, the feedback sequence {H} is obtained in order of height. k } and calculate the adjacent feedback height difference Δh k =H k+1 -H k ; Choose to generate the maximum Δh k The two adjacent feedback heights (H) k\* H k\*+1 The main signal boundary is used as the boundary, and the relative reflectance corresponding to the height index between the main signal boundaries is extracted as the normal distribution range of the forest after denoising.

5. The method for optimizing forest canopy height estimation using interferometric tomography SAR according to claim 4, characterized in that, Setting a threshold T includes: setting the threshold T and the signal loss critical value K in a unified manner, so that the envelope extraction threshold and the noise continuity screening threshold are used for noise removal and signal boundary discrimination under the same data standard.

6. The optimization method for estimating forest canopy height using interferometric tomography SAR according to claim 1, characterized in that, Pauli decomposition was performed on fully polarimetric SAR data to obtain canopy density indices, including: Perform Pauli decomposition on fully polarimetric SAR data; Pauli decomposition brightness was extracted as an indicator of forest canopy density. Based on the canopy density index, the study area was divided into sparse cover, medium cover, and dense cover levels; Output the canopy density level label for each pixel.

7. The optimization method for estimating forest canopy height using interferometric tomography SAR according to claim 6, characterized in that, Based on the canopy density grading results, different forest canopy K values ​​and surface K values ​​are assigned to pixels of different grades, including: Establish threshold sets for canopy envelope extraction and land surface envelope extraction for sparse cover, medium cover, and dense cover levels, respectively; The corresponding canopy K value and surface K value are retrieved based on the canopy density level to which each pixel belongs; The canopy envelope and the surface envelope were re-extracted using the obtained canopy K-value and surface K-value, respectively. Adaptive correction is performed on the envelope position offset, and the corrected canopy envelope height index and surface envelope height index are output.

8. The method for optimizing forest canopy height estimation using interferometric tomography SAR according to claim 1, characterized in that, The corresponding canopy K-value and surface K-value are retrieved based on the canopy density level of each pixel. The specific expression is as follows: ； H represents the height of the forest canopy; This represents the operation that minimizes the target quantity; P(·) represents the relative reflectance or the evaluation function corresponding to the relative reflectance. This represents the height difference between the canopy side and the ground side at location (r,x); K(C) represents the height of the surface envelope at location (r,x); C(C) represents the critical value corresponding to the coverage level C; C represents the forest coverage level; C(max) represents the dense coverage level; C(medium) represents the medium coverage level; C(min) represents the sparse coverage level; r,x represent the pixel spatial location index; K=0.1,0.2,0.3,⋯,0.8 indicates that different critical values ​​are searched.

9. The optimization method for estimating forest canopy height using interferometric tomography SAR according to claim 8, characterized in that, The forest canopy height is obtained by subtracting the height indices corresponding to the canopy envelope and the ground surface envelope. Specifically: Read the corrected height index corresponding to the canopy envelope and the height index corresponding to the land surface envelope; Perform a difference calculation between the height index corresponding to the canopy envelope and the height index corresponding to the land surface envelope to obtain the pixel-by-pixel forest canopy height. The specific expression is as follows: ； In the formula, This represents the height difference between the canopy side and the ground side at location (r,x); This represents the height corresponding to the surface envelope at location (r,x); Spatial aggregation was performed on all pixel-by-pixel forest canopy heights to output the forest canopy height results for the study area.

10. An optimization system for estimating forest canopy height using interferometric tomography SAR, characterized in that, The system includes: The relative reflectance calculation module is used to read L-band UAVSAR multi-baseline data, construct the multi-baseline InSAR covariance matrix based on the ground phase compensation results, and reconstruct the relative reflectance of the forest vertical profile using Capon spectrum estimation. The relative reflectance loss threshold K determination module is used to perform envelope fitting on the relative reflectance signal of the forest vertical profile, and combine LiDAR height reference data to determine the initial envelope of the forest canopy, the initial envelope of the ground surface, and the signal loss threshold K; The noise error elimination module is used to perform noise filtering based on the vertical signal continuity on the relative reflectance signal of the forest vertical profile, remove noise interference from the top of the canopy and the bottom of the ground, and obtain a denoised forest vertical structure reflectance profile and a more reasonable envelope. The forest canopy density classification module is used to perform Pauli decomposition on fully polarimetric SAR data to obtain canopy density indices, and assign different forest canopy K values ​​and surface K values ​​to pixels of different levels according to the canopy density classification results. The forest height inversion module is used to re-extract the envelope from the relative reflectivity signal and obtain the forest canopy height by subtracting the height index corresponding to the canopy envelope and the ground surface envelope.

Images