CN117434517B - Forest canopy height estimation method based on extinction coefficient optimization - Google Patents

Forest canopy height estimation method based on extinction coefficient optimization Download PDF

Info

Publication number
CN117434517B
CN117434517B CN202311768830.4A CN202311768830A CN117434517B CN 117434517 B CN117434517 B CN 117434517B CN 202311768830 A CN202311768830 A CN 202311768830A CN 117434517 B CN117434517 B CN 117434517B
Authority
CN
China
Prior art keywords
extinction coefficient
forest
phase
vegetation index
vegetation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202311768830.4A
Other languages
Chinese (zh)
Other versions
CN117434517A (en
Inventor
张国飞
岳彩荣
罗洪斌
王娟
朱俊琳
章皖秋
袁华
姬永杰
张苏云
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southwest Forestry University
Original Assignee
Southwest Forestry University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Southwest Forestry University filed Critical Southwest Forestry University
Priority to CN202311768830.4A priority Critical patent/CN117434517B/en
Publication of CN117434517A publication Critical patent/CN117434517A/en
Application granted granted Critical
Publication of CN117434517B publication Critical patent/CN117434517B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/48—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00
    • G01S7/4802—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/88—Lidar systems specially adapted for specific applications
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/88—Lidar systems specially adapted for specific applications
    • G01S17/89—Lidar systems specially adapted for specific applications for mapping or imaging
    • G01S17/90—Lidar systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/48—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00
    • G01S7/4808—Evaluating distance, position or velocity data
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/70—Determining position or orientation of objects or cameras
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/13—Satellite images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/188—Vegetation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Electromagnetism (AREA)
  • Software Systems (AREA)
  • Multimedia (AREA)
  • General Health & Medical Sciences (AREA)
  • Astronomy & Astrophysics (AREA)
  • Artificial Intelligence (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Medical Informatics (AREA)
  • Health & Medical Sciences (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Image Processing (AREA)

Abstract

The invention relates to the technical field of canopy height estimation, in particular to a forest canopy height estimation method based on extinction coefficient optimization, which comprises the following steps: (1) data preparation; (2) Fitting a complex coherence straight line by adopting a total least square method to obtain a candidate ground phase; (3) Ground phase simulation using LiDAR DTM and satellite orbit data, then referencing ground simulation phaseEstimated earth phaseThe method comprises the steps of carrying out a first treatment on the surface of the (4) obtaining a complex coherence estimate of pure body scattering; (5) Estimating extinction coefficientThe method comprises the steps of carrying out a first treatment on the surface of the (6) Establishing a relation model of leaf area index, vegetation coverage and optical/SAR remote sensing vegetation index and extinction coefficient of the laser radar; (7) constructing an extinction coefficient finite interval: (8) And inverting the limited interval of the forest extinction coefficient through the remote sensing vegetation index, and bringing the limited interval into the third stage of the three-stage inversion algorithm. The method can improve the accuracy of forest canopy height estimation.

Description

Forest canopy height estimation method based on extinction coefficient optimization
Technical Field
The invention relates to the technical field of canopy height estimation, in particular to a forest canopy height estimation method based on extinction coefficient optimization.
Background
The extinction coefficient refers to the attenuation (attenuation coefficient) of electromagnetic waves in a unit distance in a propagation medium, is expressed in dB/m, is an important factor of radiation transmission, and relates to the frequency of microwave signals, forest stand environments and the like. In the RVoG model three-stage solving algorithm, the extinction coefficient is one of important parameters affecting the solving accuracy.
In the RVoG model, the extinction coefficient range is mainlyIs affected by the frequency of the signal. Extinction coefficientIs selected from a wide range of values, i.e. calculated bulk scattering complex coherence +.>Based on pure body scattering complex coherence candidate and body scattering complex coherence observation value in LUT only +.>Bulk scattering complex coherence +.>. In the existing forest canopy height inversion process, many scholars take the extinction coefficient to be a fixed value, the fixed value of the extinction coefficient of needle forests is 0.2 dB/m, and the fixed values of the extinction coefficient of tropical forests are 0.3 dB/m and 0.4 dB/m. Also, the RVoG model has not considered that there is vertical isomerism in forest structures, and it is assumed that an extinction coefficient of a constant value is not suitable for high and dense forest areas such as tropical rain forests.
In practice, the real scattering process of SAR in a forest is very complex, and the extinction coefficient is not only affected by the signal frequency. The extinction coefficient is a function of biophysical properties such as band frequency, polarization mode, forest height, density, etc. It is a function of forest properties, related to penetration depth. Furthermore, the bulk scattering miscibility increases with increasing penetration depth and changes transiently the bulk scattering complex coherence. This means that the higher the extinction coefficient, the less permeable to forest media, and the lower the effective bulk scattering loss, and hence the extinction coefficient is inversely related to penetration depth.
Extinction coefficient is an important parameter in three-stage inversion accuracy, and finding a reliable range of extinction coefficients is a controversial challenge in inversion algorithms. .
Disclosure of Invention
The invention provides a forest canopy height estimation method based on extinction coefficient optimization, which can overcome certain or some defects in the prior art.
The invention relates to a forest canopy height estimation method based on extinction coefficient optimization, which comprises the following steps of:
(1) Preparing data;
calculating 13 polarization scattering mechanisms of the TanDEM-X full polarization SAR interference data, namely: HH. HV, VH, VV, HH +VV, HH-VV, LL, RR, OPT, OPT2, OPT3, PDhigh, PDlow polarization scattering mechanism complex coherence; obtaining gradient angle of distance direction by using LiDAR DTMCorrecting local angle of incidence +.>Vertical wave number->；
(2) Complex coherent straight line CFL fitting to obtain candidate ground phase；
Fitting a complex coherence straight line CFL by adopting a TLS (TLS total least squares) method; two intersection point phases of the complex coherence straight line CFL and the complex coherence unit circle CUC are used as ground phase candidate values;
(3) Ground phase simulation using LiDAR DTM and satellite orbit data, then referencing ground simulation phaseEstimated ground phase->；
(4) Obtaining a complex coherence estimation value of pure body scattering;
(5) Estimating extinction coefficient using a small number of LiDAR CHM samples and F-functions for forest canopy height；
(6) Establishing a relation model of leaf area index, vegetation coverage and optical/SAR remote sensing vegetation index and extinction coefficient of the laser radar through a random forest model;
(7) Inverting the extinction coefficient of the forest according to the relation model, and constructing an extinction coefficient limited interval by the following formula:
；
(8) And inverting the limited interval of the forest extinction coefficient through the remote sensing vegetation index, and bringing the limited interval into a third stage of a three-stage inversion algorithm, so that the accuracy of inverting the forest canopy height by the RVoG model is further improved.
Preferably, in the step (3), the method specifically includes the following steps:
(3.1) fitting a satellite orbit by adopting header file data of SAR satellites;
(3.2) calculating the pixel slant distance by adopting an R-D positioning model and LiDAR DTM;
(3.3) analog ground phase under the effect of the De-level ground phase；
(3.4) reference to the phase of the analog groundEstimated ground phase->。
Preferably, in the step (3.1), specifically:
providing N time position vectors in SAR data header fileAnd velocity vectorS is satellite position, and a polynomial orbit description method and position and speed vectors of discrete points are adopted to fit satellite orbits, so that satellite position vectors and speed vectors at any moment are obtained;
the imaging time of the ith row of the remote sensing image is calculated by adopting the following steps:
；
in the method, in the process of the invention,for reference time->Imaging time for row i,/->Is the azimuth direction; n is the total number of lines>For pulse repetition frequency, +.>And->Obtained from the header file.
Preferably, the position vector uses a third order polynomial to perform satellite orbit simulation, and the velocity vector uses a second order polynomial to perform satellite orbit simulation, as follows:
；
；
in the method, in the process of the invention, is the estimated coefficient.
Preferably, in the step (3.2), specifically:
(a) Calculating radius R of approximate circle of remote sensing image L ；
；
In the method, in the process of the invention,is a first ellipsoid parameter; />Is a second ellipsoid parameter; />The central latitude of the remote sensing image is; cos () is a cosine function, sin () is a sine function;
(b) Calculating a first geocentric angle；
；
In the method, in the process of the invention,is a first geocentric angle; />Is a SAR satellite position vector; r is (r) 0 Obtaining a satellite head file for a first slant distance; cos -1 () Taking an inverse cosine function;
(c) Calculating the slant distance of any pixel;
；
in the method, in the process of the invention,is a picture elementA slant distance; />For line number, ->Is a column number; />The size of the distance-to-ground pixel;
(d) Solving by an indirect positioning method;
and obtaining a geodetic coordinate through LiDAR DEM images, then converting into a space rectangular coordinate system, establishing an R-D positioning model by combining fitted satellite orbits, and carrying out iterative solution to obtain image row-column coordinates (i, j).
Preferably, in the step (3.3), specifically:
in the process of obtaining the pixel slant distanceAfter that, the overall phase deviation +.>Then, the following steps are carried out:
；
in the method, in the process of the invention,is a high blur number +.>Satellite microwave is a satellite moving deviation distance of a time difference from transmitting to receiving, < >>Is the flat ground phase,/->Is an analog ground phase.
Preferably, in the step (3.4), specifically:
first, calculate the analog ground phasePhase with two candidates->Phase difference between->The following formula:
；
second step, comparing: if->When in use, then->The method comprises the steps of carrying out a first treatment on the surface of the If->When in use, then->The following formula:
。
preferably, in step (5), the extinction coefficient is estimated using a F-functionThe formula of (2) is:
；
for bulk scattering complex coherence->Is vertical wave number, < >>For the angle of incidence of the radar wave, the height obtained at this time is +.>The height of the canopy of the forest is required to be converted into +.>I.e. +.>。
Preferably, in the step (6), after denoising, filtering, normalizing and the like are performed on the LiDAR point cloud data, a Leaf Area Index (LAI) and a vegetation coverage (Fraction of Vegetation Cover, fCover) are calculated;
the leaf area index LAI refers to the leaf area of a plant per unit land area, and is calculated by the following formula:
；
in the method, in the process of the invention,is the average value of the scan angle,/">Is->Scan angle values of individual points +.>Is the gap rate>Is an extinction coefficient, ln () is a logarithm, cos () is a cosine;
fCover refers to the vertical projected area of a plant per unit area, calculated from the following formula:
；
in the method, in the process of the invention,is the number of vegetation points->Is the total points.
Preferably, the optical/SAR remote sensing vegetation index comprises:
a. 6 variables generated based on Sentinel-2 multispectral data: the vegetation photosynthetic effective radiation absorption ratio, the normalized vegetation index, the difference vegetation index, the normalized difference index algorithm NDI45, the transformation type normalized vegetation index and the atmospheric resistance vegetation index;
b. 2 radar vegetation index variables generated based on ALOS-2/PALSAR-2: freeman_rvi radar vegetation index freeman_alos, kim_rvi radar vegetation index kim_alos;
c. 2 radar vegetation index variables generated based on TanDEM-X SAR: freeman_rvi radar vegetation index freeman_tanx and kim_rvi radar vegetation index kim_tanx.
The method can accurately estimate the ground phase, thereby improving the accuracy of forest canopy height estimation. According to the invention, the relation model of the remote sensing vegetation indexes such as leaf area index, forest coverage, normalized vegetation index, difference vegetation index, radar vegetation index and the like and the extinction coefficient is researched and explored through a machine learning model, the range of the extinction coefficient value is shortened, and the accuracy of the three-stage algorithm inversion forest canopy height of the RVoG model is further improved.
Drawings
FIG. 1 is a flowchart of a forest canopy height estimation method based on extinction coefficient optimization in an embodiment;
FIG. 2 (a) is a schematic diagram of LiDARDTM after windowing in an embodiment;
FIG. 2 (b) is a schematic diagram of LiDAR DTM after matching with SAR image in the embodiment;
FIG. 3 is a flow chart illustrating the intermediate positioning method according to the embodiment;
FIG. 4 is a schematic diagram of simulated ground phase in an embodiment;
FIG. 5 is a graph showing the distribution of 12 vegetation index in the example;
FIG. 6 is a plot of the peak forest canopy height scatter after the extinction coefficient correction in the example.
Description of the embodiments
For a further understanding of the present invention, the present invention will be described in detail with reference to the drawings and examples. It is to be understood that the examples are illustrative of the present invention and are not intended to be limiting.
Examples
As shown in fig. 1, the embodiment provides a forest canopy height estimation method based on extinction coefficient optimization, which includes the following steps:
(1) Preparing data;
calculating 13 polarization scattering mechanisms of the TanDEM-X full polarization SAR interference data, namely: HH. HV, VH, VV, HH +VV, HH-VV, LL, RR, OPT, OPT2, OPT3, PDhigh, PDlow polarization scattering mechanism complex coherence; obtaining gradient angle of distance direction by using LiDAR DTMCorrecting local angle of incidence (+)>) Vertical wave number->；
(2) Complex coherent straight line CFL fitting to obtain candidate ground phase；
Fitting a complex coherent straight line (Complex Fitting Line, CFL) by using a TLS total least squares method; two intersection phases of the complex coherence straight line CFL and the complex coherence unit circle (Complex Unit Circle, CUC) are used as ground phase candidate values;
(3) Ground phase simulation using LiDAR DTM and satellite orbit data, then referencing ground simulation phaseEstimated ground phase->；
In the step (3), the method specifically comprises the following steps:
(3.1) fitting a satellite orbit by adopting header file data of SAR satellites;
in the step (3.1), specifically:
in the SAR satellite geometric conformation model, S is the satellite position and the position three-dimensional vector thereofAnd velocity vector->The accuracy of the SAR conformational model is affected by the accuracy of the position vector and the speed vector; t is a ground target object, the position vector of which is +.>Projection of +.>，/>The elevation of (2) is。
Providing N time position vectors in SAR data header fileAnd velocity vectorS is the satellite position, and the embodiment adopts a polynomial orbit description method and the position and speed vectors of discrete points to fit the satellite orbit, so that satellite position vectors and speed vectors at any moment are obtained;
the position vector uses a cubic polynomial to simulate the satellite orbit, and the velocity vector uses a quadratic polynomial to simulate the satellite orbit, as follows:
；
；
in the method, in the process of the invention, is the estimated coefficient. The terrSAR X image data head file provides satellite orbit data of 77 moments, and the least square method is adopted to obtain the medium estimation coefficient so as to obtain a satellite position state vector and a satellite speed vector at any moment.
；
；
And inverting the satellite speed vector and the position vector by using a polynomial orbit fitting algorithm. The imaging time of the ith row of the remote sensing image is calculated by adopting the following steps:
；
in the method, in the process of the invention,for reference time->Imaging time for row i (pointing direction,) x +.>Is the azimuth direction; n is the total number of lines (pointing direction),>for pulse repetition frequency, +.>And->Obtained from the header file.
(3.2) calculating the pixel slant distance by adopting an R-D positioning model and LiDAR DTM;
DEM sampling
In general, the spatial resolution (less than 10 m) of the high-resolution remote sensing image is higher than the spatial resolution (30 m) of the original DEM image, and there is a problem of undersampling. In this embodiment, the spatial resolution of the obtained airborne LiDAR DTM is 1m, and there is an oversampling phenomenon, and the LiDAR DTM opens a window of [6*6], as shown in FIG. 2 (a).
The images obtained by projecting the LiDAR DTM windowed in FIG. 2 (a) onto the corresponding SAR image are shown in FIG. 2 (b).
Resolving the pixel slant distance
Inclined distance and azimuth direction of satellite and ground surface target object) Irrespective of column number (+)>) Related to the following.
(a) Calculating radius R of approximate circle of remote sensing image L
；
In the method, in the process of the invention,is a first ellipsoid parameter; />Is a second ellipsoid parameter; for WGS-84 ellipsoids, a first ellipsoial parameter R a 6378137m; second ellipsoid parameter R b 6356752.3141m; />The central latitude of the remote sensing image is; cos () is a cosine function, sin () is a sine function;
(b) Calculating a first geocentric angle；
；
In the method, in the process of the invention,is a first geocentric angle; />Is a SAR satellite position vector; r is (r) 0 Obtaining a satellite head file for a first slant distance; cos -1 () Taking an inverse cosine function;
(c) Calculating the slant distance of any pixel
；
In the method, in the process of the invention,is the pixel slant distance; />For line number, ->Is a column number; />The size of the distance-to-ground pixel;
(d) Solution by indirect localization method
Due to the topographic effect, the SAR image has the phenomena of perspective shrinkage, overlapping masking, shadow and the like, so that the original SAR image coordinates and the DEM geographic coordinates are not in one-to-one correspondence. In this embodiment, an indirect positioning method is used to calculate the mapping relationship between the DEM geographic coordinates and the original SAR image coordinates, and a specific calculation flow is shown in fig. 3. The geodetic coordinates are obtained through LiDAR DEM images, then a space rectangular coordinate system is transformed, an R-D positioning model is established by combining fitted satellite orbits, and iterative solution is carried out, so that image row-column coordinates are obtained。
(3.3) analog ground phase under the effect of the De-level ground phaseAs shown in fig. 4, the left side is the ground phase from which the ground phase is not removed, the middle is the ground phase, and the right side is the ground phase from which the ground phase is removed.
In the step (3.3), specifically:
in the process of obtaining the pixel slant distanceAfterwards, the global phase deviation (++) is removed according to the following formula>) Then, the following steps are carried out:
；
where HOA is a high degree of ambiguity (representing the difference in height resulting in 2 interferometric phase differences)，/>Satellite microwave is a satellite moving deviation distance of a time difference from transmitting to receiving, < >>Is the flat ground phase,/->Is an analog ground phase.
(3.4) reference to the phase of the analog groundEstimated ground phase->。
In the step (3.4), specifically:
first, calculate the analog ground phasePhase with two candidates->Phase difference between->The following formula:
；
second step, comparing: if->When in use, then->The method comprises the steps of carrying out a first treatment on the surface of the If it isWhen in use, then->The following formula:
。
(4) Obtaining a complex coherence estimation value of pure body scattering;
after obtaining accurate phase, the volume scattering complex coherence needs to be solved before the forest canopy height inversion. In the classical three-phase algorithm,as a complex coherent observation of volume scattering +.>Its perpendicular projection on complex coherence fitted line (complex fitting line, CFL) is used as an estimate of the volume scattering complex coherence +.>. The wavelength of the X wave band is short, the penetrability is small, the dominant scattering in the forest is crown body scattering, the body scattering signal mainly comes from the branches and leaves of the crown layer of the tree, the HV polarization channel is influenced by the fluctuation of the topography and the uncertainty of the growth direction of the branches and leaves, and other scattering contributions can be contained. In addition, the vertical projection method changes the mode value and the phase of the effective observation value of the bulk scattering complex coherence, and errors are easy to introduce. The embodiment adopts the optimized estimation method of the pure body scattering of the previous research results (Anhui autumn, 2018; zhang Guofei, etc., 2022), namely, the effective observation value of the pure body scattering is screened out from the polarization complex coherence +.>And obtaining the pure body scattering estimated value by adopting a model value invariable projection method.
(5) Estimation of extinction coefficient using LiDAR data-derived forest canopy height (LiDAR CHM) and F-function with a small number of samples；
According to RVoG model (RVoG model passing vegetation thickness)Extinction coefficient->Effectively bulk scatter->Ground phase->Four parameters simulate pure body scattering complex coherence, which is the prior art, and multiple sets of extinction coefficients are calculated (+.>) Comparing the complex coherence theory value with the complex coherence estimate value, and estimating the extinction coefficient of the effective area by using the function F and the small sample forest canopy height (LiDAR CHM)>；
In step (5), the extinction coefficient is estimated by using F functionThe formula of (2) is:
;
for bulk scattering complex coherence->Is vertical wave number, < >>For radar wave entryFiring angle, height obtained at this time +.>(gradient to volume thickness) to be converted into forest canopy height +.>I.e. +.>。
(6) Establishing a relation model of leaf area index, vegetation coverage and optical/SAR remote sensing vegetation index and extinction coefficient of the laser radar through a random forest model;
in this example, there are 12 vegetation index variables, mainly including:
1. LiDAR generation-based 2 variables
After denoising, filtering, normalizing and the like are performed on LiDAR point cloud data of a research Area by using LiDAR360 software, leaf Area Index (LAI) and vegetation coverage (Fraction of Vegetation Cover, fCover) are calculated.
Leaf Area Index (LAI): refers to the plant leaf area per unit area of land. The leaf area index distribution in the study area is calculated from the following formula and is shown as LAI in fig. 5.
;
In the method, in the process of the invention,bis the average value of the scan angle,b i is the firstiThe scan angle value of the individual points,f v is the gap rate of the material to be processed,is an extinction coefficient, ln () is a logarithm, cos () is a cosine; the leaf area index distribution of the study area is shown as LAI in figure 5.
Vegetation coverage (Fraction of Vegetation Cover, fcovery): the vertical projection area of the plants in the unit area of the research area is calculated by the following formula, and the vegetation coverage distribution of the research area is shown as fCover in figure 5.
In the method, in the process of the invention,is the number of vegetation points->Is the total points.
2. 6 variables generated based on Sentinel-2 multispectral data
Multispectral imaging satellite guard No. 2 (Sentinel-2) has 13 wave bands of Short Wave Infrared (SWIR), near Infrared (NIR), visible light and the like. The Sentinel-2 multispectral image is processed by SNAP software to obtain 6 vegetation index changesMeasuring amount。
Ratio of photosynthetic effective radiation absorption of vegetation (Fraction of absorbed photosynthetically activer radiation, fPAR): defined as the absorption of solar radiation (wavelength=400 nm-700 nm) energy by vegetation per unit area, and the distribution of the photosynthetic active radiation absorption ratio of vegetation in the investigation region is shown as fPAR in fig. 5.
Normalized vegetation index (Normalized Vegetation Index, NDVI): the normalized vegetation index algorithm utilizes the intensity and vitality of the surface vegetation. Spectral features of healthy vegetation show a sudden rise in reflection level at 0.7 μm, while the land without vegetation has a continuous linear course, depending on the surface type. The higher the chlorophyll activity of the plant, the greater the increase in the near infrared (0.78-1 μm) reflection level. In addition to determining the distance between vegetation and other objects, it can also detect the vitality of vegetation. In the red region, the incident solar radiation is not absorbed extensively by pigments in mesophylls, mainly by chlorophyll. In contrast, in the near infrared region, most of the incident radiation is reflected. NDVI constitutes a measure of photosynthetic activity and is closely related to the density and vigor of vegetation. Normalization reduces the effects of topography and atmosphere and enables large-area simultaneous inspection, calculated from the following equation, with the normalized vegetation index distribution of the study area shown in NDVI in fig. 5.
；
Differential vegetation index (Difference Vegetation Index, DVI): in the red region, the incident solar radiation is not absorbed extensively by pigments in mesophylls, mainly by chlorophyll. In contrast, in the near infrared region, most of the incident radiation was reflected, calculated from the following equation, and the study area difference vegetation index distribution is shown in DVI in fig. 5.
；
Normalized differential exponential algorithm NDI45 (Normalized Difference Index 45, NDI 45): the normalized differential index algorithm is more linear than the normalized vegetation index, has lower saturation at higher values, and is calculated by the following formula, and the distribution of the normalized differential index algorithm in the research area is shown as NDI45 in fig. 5.
；
Transformed normalized vegetation index (Transformed Normalised Difference Vegetation index, TNDVI) the transformed normalized vegetation index is the square root of the normalized differential vegetation index. The TNDVI formula is always positive, the variance of the ratio is proportional to the mean value, and the TNDVI formula is calculated by the following formula, and the TNDVI formula is shown in the figure 5.
；
The atmospheric resistance vegetation index (Atmospherically Resistant Vegetation Index, ARVI) is calculated by adopting the following formula by adopting the corrected difference value of the blue light and red light wave band reflectivities, and the distribution of the atmospheric resistance vegetation index in the research area is shown in the ARVI in figure 5.
；
In the method, in the process of the invention,is the reflectivity of near infrared, red and blue bands after correction. />Is a weighting function and is dependent on the aerosol type.
3. ALOS-2/PALSAR-2 based generation of 2 radar vegetation index variables
The L-band full-polarization synthetic aperture radar sensor (PALSAR-2) carried by ALOS-2/PALSAR-2SAR has a satellite revisit period of 14 d.
Freeman_rvi radar vegetation index: f obtained by SAR polarization decomposition d (odd scattering), F s (even scattering) and F v The radar vegetation index constructed by the scattered power component (volume scattering) is calculated by the following formula, and the freeman_RVI radar vegetation index distribution of the research area is shown in the freeman_alos in figure 5.
；
Kim_rvi radar vegetation index: from the following componentsThe backscattering coefficient is calculated and obtained by the following formula, and the Kim_Alos of the research area Kim_RVI radar vegetation index distribution is shown in figure 5.
。
4. 2 radar vegetation index variables generated based on TanDEM-X SAR
Generating 2 radar vegetation indexes based on the TanDEM-X SAR data: the Freeman_RVI radar vegetation index and Kim_RVI radar vegetation index are the same as ALOS-2/PALSAR-2. The Freeman RVI radar vegetation index distribution of the research area is shown in the Freeman tan X in figure 5; the Kim_RVI radar vegetation index distribution of the study area is shown in FIG. 5 as Kim_tan X.
(7) Inverting the extinction coefficient of the forest according to the relation model, and constructing an extinction coefficient limited interval by the following formula:
；
machine learning model
In a real forest scene, the SAR is quite complex in a forest scattering process, and an extinction coefficient is a function of biophysical properties such as band frequency, polarization mode, forest height, density and the like, and is an important factor of radiation transmission. The influence of vegetation indexes calculated based on optical remote sensing, liDAR, SAR and other data on extinction coefficients can have a nonlinear relation; in the embodiment, random forest model estimated extinction coefficient experimental study is adopted. The random forest can better establish the corresponding relation between the characteristic variable and the estimated target from limited training samples.
The random forest is an integrated prediction model in which a plurality of decision trees are arranged in parallel according to random arrangement, each decision tree calculates random vector values of independent sampling, and an average value of predicted result values of a plurality of decision trees is taken as a final result to be output. When the number of decision trees is larger, the generalization error of the random forest gradually converges. Factors such as the number of decision trees, the strength of the individual decision trees, the correlation between decision trees, etc., will affect the generalization error range. The prior researches show that the random forest model has the advantages of high learning speed, wide data processing range and higher prediction accuracy.
The random forest model plays two roles in this embodiment: (1) looking at the correlation of 12 vegetation index variables with extinction coefficient by means of the Gini index (Gini index); (2) and constructing a machine model of 12 vegetation index variables and extinction coefficients by adopting training samples, and estimating a range interval of the extinction coefficients.
(8) And inverting the limited interval of the forest extinction coefficient through the remote sensing vegetation index, and bringing the limited interval into a third stage of a three-stage inversion algorithm, so that the accuracy of inverting the forest canopy height by the RVoG model is further improved. FIG. 6 is a plot of the peak forest canopy height scatter after extinction coefficient correction. The correlation between the estimated value of the forest canopy height corrected based on the extinction coefficient and the measured forest height is improved from 0.75 to 0.82, the inversion accuracy is improved by 0.42m, the Average Absolute Error (AAE) is improved to 1.81m from 2.15m, and the accuracy percentage is improved to 85.68% from 83.05%.
The embodiment can accurately estimate the ground phase, thereby improving the accuracy of forest canopy height estimation. According to the embodiment, the relation model of the remote sensing vegetation indexes such as leaf area indexes, forest coverage, normalized vegetation indexes, difference vegetation indexes and radar vegetation indexes and the extinction coefficient is researched and explored through a machine learning model, the range of the extinction coefficient is shortened, and therefore the accuracy of inversion of forest canopy height by an RVoG model three-stage algorithm is further improved.
The invention and its embodiments have been described above by way of illustration and not limitation, and the invention is illustrated in the accompanying drawings and described in the drawings in which the actual structure is not limited thereto. Therefore, if one of ordinary skill in the art is informed by this disclosure, the structural mode and the embodiments similar to the technical scheme are not creatively designed without departing from the gist of the present invention.

Claims (4)

1. A forest canopy height estimation method based on extinction coefficient optimization is characterized in that: the method comprises the following steps:
(1) Preparing data;
calculating 13 polarization scattering mechanisms of the TanDEM-X full polarization SAR interference data, namely: HH. HV, VH, VV, HH +VV, HH-VV, LL, RR, OPT, OPT2, OPT3, PDhigh, PDlow polarization scattering mechanism complex coherence; obtaining gradient angle alpha of distance direction by LiDARDTM, correcting local incident angle theta-alpha and vertical wave number k z ；
(2) Complex coherent straight line CFL fitting to obtain candidate ground phase
Fitting a complex coherence straight line CFL by adopting a TLS (TLS total least squares) method; two intersection point phases of the complex coherence straight line CFL and the complex coherence unit circle CUC are used as ground phase candidate values;
(3) Ground phase simulation using LiDAR DTM and satellite orbit data, then referencing ground simulation phaseEstimated ground phase->
(4) Obtaining a complex coherence estimation value of pure body scattering;
(5) Estimating an extinction coefficient sigma by adopting a small quantity of LiDAR CHM samples of forest canopy height and an F function;
(6) Establishing a relation model of leaf area index, vegetation coverage and optical/SAR remote sensing vegetation index and extinction coefficient of the laser radar through a random forest model;
(7) Inverting the extinction coefficient of the forest according to the relation model, and constructing an extinction coefficient limited interval by the following formula:
(8) The limited interval of the forest extinction coefficient is inverted through the remote sensing vegetation index and is brought into the third stage of the three-stage inversion algorithm, so that the accuracy of inversion of forest canopy height by the RVoG model is further improved;
in the step (3), the method specifically comprises the following steps:
(3.1) fitting a satellite orbit by adopting header file data of SAR satellites;
(3.2) calculating the pixel slant distance by adopting an R-D positioning model and LiDAR DTM;
(3.3) analog ground phase under the effect of the De-level ground phase
(3.4) reference to the phase of the analog groundEstimated ground phase->
In the step (3.1), specifically:
providing N time position vectors in SAR data header fileAnd velocity vectorS is satellite position, and a polynomial orbit description method and position and speed vectors of discrete points are adopted to fit satellite orbits, so that satellite position vectors and speed vectors at any moment are obtained;
the imaging time of the ith row of the remote sensing image is calculated by adopting the following steps:
wherein t is 0 For reference time, t i Imaging time of the ith row, wherein i is azimuth; n is the total number of rows, PRF is the pulse repetition frequency, t 0 And PRF obtained from the header file;
the position vector uses a cubic polynomial to simulate the satellite orbit, and the velocity vector uses a quadratic polynomial to simulate the satellite orbit, as follows:
X S ＝a 00 +a 01 t+a 02 t 2 +a 03 t 3
Y S ＝a 10 +a 11 t+a 12 t 2 +a 13 t 3
Z S ＝a 20 +a 21 t+a 22 t 2 +a 23 t 3
V X ＝a 01 +2a 02 t+3a 03 t 2
V Y ＝a 11 +2a 12 t+3a 13 t 2
V Z ＝a 21 +2a 22 t+3a 23 t 2
wherein a is 00 、a 01 、a 02 、a 03 、a 10 、a 11 、a 12 、a 13 、a 20 、a 21 、a 22 、a 23 Is an estimated coefficient;
in the step (3.2), specifically:
(a) Calculating radius R of approximate circle of remote sensing image L ；
Wherein R is a Is a first ellipsoid parameter; r is R b Is a second ellipsoid parameter; delta is the latitude of the center of the remote sensing image; cos () is a cosine function, sin () is a sine function;
(b) Calculating a first geocentric angle beta 0 ；
Wherein beta is 0 Is a first geocentric angle;is a SAR satellite position vector; r is (r) 0 Obtaining a satellite head file for a first slant distance; cos -1 () Taking an inverse cosine function;
(c) Calculating the slant distance of any pixel;
wherein R (i, j) is the pixel slant distance;i is a row number, j is a column number; delta r The size of the distance-to-ground pixel;
(d) Solving by an indirect positioning method;
obtaining a geodetic coordinate through a LiDAR DEM image, then converting into a space rectangular coordinate system, establishing an R-D positioning model by combining a fitted satellite orbit, and carrying out iterative solution to obtain an image row-column coordinate (i, j);
in the step (3.3), specifically:
after the pixel skew R (i, j) is found, the overall phase deviation n×2pi is removed according to the following equation:
wherein HOA is a high blur number, S MC Is the satellite movement offset distance of the satellite microwaves from the transmission to the reception time difference,is the flat ground phase,/->Is an analog ground phase;
in the step (3.4), specifically:
first, calculate the analog ground phasePhase with two candidates->Phase difference between->The following formula:
second step, comparingAnd->If->When in use, then->If->When in use, then->The following formula:
2. the method for estimating the height of the canopy of the forest based on the extinction coefficient optimization according to claim 1, wherein the method comprises the following steps: in step (5), the formula for estimating the extinction coefficient σ by using the F function is:
γ v for bulk scattering complex coherence, k z Is the vertical wave number, θ is the incident angle of the radar wave, and the height h is obtained at this time v To be converted into the forest canopy height h in the horizontal state f I.e.
3. The method for estimating the height of the canopy of the forest based on the extinction coefficient optimization according to claim 2, wherein the method comprises the following steps: in the step (6), after denoising, filtering and normalizing the LiDAR point cloud data, calculating a leaf area index LAI and a vegetation coverage fCover;
LAI refers to the plant leaf area per unit land area calculated from the following formula:
where θ is the average value of the scan angle, θ i Is the scan angle value of the i-th point, f v Is the gap rate, σ is the extinction coefficient, ln () is the logarithm, cos () is the cosine;
fCover refers to the vertical projected area of a plant per unit area, calculated from the following formula:
wherein n is veg Is the number of vegetation points, n total Is the total points.
4. A method for estimating a forest canopy height based on extinction coefficient optimization as claimed in claim 3, wherein: the optical/SAR remote sensing vegetation index comprises:
a. 6 variables generated based on Sentinel-2 multispectral data: the vegetation photosynthetic effective radiation absorption ratio, the normalized vegetation index, the difference vegetation index, the normalized difference index algorithm NDI45, the transformation type normalized vegetation index and the atmospheric resistance vegetation index;
b. 2 radar vegetation index variables generated based on ALOS-2/PALSAR-2: freeman_rvi radar vegetation index freeman_alos, kim_rvi radar vegetation index kim_alos;
c. 2 radar vegetation index variables generated based on TanDEM-X SAR: freeman_rvi radar vegetation index freeman_tanx and kim_rvi radar vegetation index kim_tanx.
CN202311768830.4A 2023-12-21 2023-12-21 Forest canopy height estimation method based on extinction coefficient optimization Active CN117434517B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311768830.4A CN117434517B (en) 2023-12-21 2023-12-21 Forest canopy height estimation method based on extinction coefficient optimization

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311768830.4A CN117434517B (en) 2023-12-21 2023-12-21 Forest canopy height estimation method based on extinction coefficient optimization

Publications (2)

Publication Number Publication Date
CN117434517A CN117434517A (en) 2024-01-23
CN117434517B true CN117434517B (en) 2024-02-27

Family

ID=89550229

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311768830.4A Active CN117434517B (en) 2023-12-21 2023-12-21 Forest canopy height estimation method based on extinction coefficient optimization

Country Status (1)

Country Link
CN (1) CN117434517B (en)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102914501A (en) * 2012-07-26 2013-02-06 南京大学 Method for calculating extinction coefficients of three-dimensional forest canopy by using laser-point cloud
CN103376451A (en) * 2012-04-16 2013-10-30 中国科学院电子学研究所 Airborne double-waveband synthetic aperture radar system and method for measuring vegetation thickness utilizing same
CN110569624A (en) * 2019-09-20 2019-12-13 哈尔滨工业大学 Forest three-layer scattering model determining and analyzing method suitable for PolInSAR inversion
CN110703220A (en) * 2019-10-12 2020-01-17 中南大学 Multi-baseline PolInSAR vegetation parameter inversion method considering time decoherence factors
CN110988879A (en) * 2019-12-24 2020-04-10 中南大学 Vegetation parameter inversion method, terminal equipment and storage medium
CN113945926A (en) * 2021-09-17 2022-01-18 西南林业大学 Forest canopy height inversion method improved through underestimation compensation
CN113945927A (en) * 2021-09-17 2022-01-18 西南林业大学 Forest canopy height inversion method through volume scattering optimization
CN115062260A (en) * 2022-06-16 2022-09-16 电子科技大学 Forest biomass PolInSAR estimation method and system suitable for heterogeneous forest and storage medium

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9086488B2 (en) * 2010-04-20 2015-07-21 Michigan Aerospace Corporation Atmospheric measurement system and method
US10261006B2 (en) * 2016-07-21 2019-04-16 Rosemount Aerospace, Inc. Method of estimating cloud particle sizes using LIDAR ratio

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103376451A (en) * 2012-04-16 2013-10-30 中国科学院电子学研究所 Airborne double-waveband synthetic aperture radar system and method for measuring vegetation thickness utilizing same
CN102914501A (en) * 2012-07-26 2013-02-06 南京大学 Method for calculating extinction coefficients of three-dimensional forest canopy by using laser-point cloud
CN110569624A (en) * 2019-09-20 2019-12-13 哈尔滨工业大学 Forest three-layer scattering model determining and analyzing method suitable for PolInSAR inversion
CN110703220A (en) * 2019-10-12 2020-01-17 中南大学 Multi-baseline PolInSAR vegetation parameter inversion method considering time decoherence factors
CN110988879A (en) * 2019-12-24 2020-04-10 中南大学 Vegetation parameter inversion method, terminal equipment and storage medium
CN113945926A (en) * 2021-09-17 2022-01-18 西南林业大学 Forest canopy height inversion method improved through underestimation compensation
CN113945927A (en) * 2021-09-17 2022-01-18 西南林业大学 Forest canopy height inversion method through volume scattering optimization
CN115062260A (en) * 2022-06-16 2022-09-16 电子科技大学 Forest biomass PolInSAR estimation method and system suitable for heterogeneous forest and storage medium

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
Examination of the extinction coefficient in the Beer-Lambert law for an accurate estimation of the forest canopy leaf area index;Taku M;《Forest Science and Technology》;20121231;全文 *
TanDEM-X极化干涉SAR森林冠层高度反演;章皖秋;岳彩荣;颜培东;;东北林业大学学报;20221231(01);全文 *
基于TanDEM-X相干系数的森林高度估测方法;范亚雄;陈尔学;李增元;赵磊;张王菲;金玉栋;蔡丽杰;;林业科学;20200615(06);全文 *
基于X波段极化干涉SAR数据的思茅松林冠层高度反演;陈子怡;章皖秋;岳彩荣;;林业调查规划;20180615(03);全文 *
干涉、极化干涉SAR技术森林高度估测算法研究进展;张王菲;陈尔学;李增元;赵磊;姬永杰;;遥感技术与应用;20171215(06);全文 *
赖格英,曾祥贵,刘影,等.基于ETM和图像融合的优势植被冠层叶面积指数和消光系数的遥感反演.《遥感技术与应用》.2013,全文. *

Also Published As

Publication number Publication date
CN117434517A (en) 2024-01-23

Similar Documents

Publication Publication Date Title
Smith et al. A model relating VHF-band backscatter to stem volume of coniferous boreal forest
CN110865346B (en) Satellite-borne SAR time parameter calibration method based on direct positioning algorithm
Axelsson Beam characteristics of three-dimensional SAR in curved or random paths
CN109932719A (en) RCS high-precision measuring method based on SAR imaging
Ferro-Famil et al. Synthetic aperture radar imaging
Lebsock et al. The feasibility of water vapor sounding of the cloudy boundary layer using a differential absorption radar technique
JP3783058B2 (en) Method and system for inverse estimation of wave direction spectrum from radar image
Laursen et al. Wind direction over the ocean determined by an airborne, imaging, polarimetric radiometer system
CN117434517B (en) Forest canopy height estimation method based on extinction coefficient optimization
Baldini et al. Analysis of dual polarization images of precipitating clouds collected by the COSMO SkyMed constellation
CN113064161B (en) Wave spectrometer cross spectrum calculation method based on double sub-pulse reconstruction
Bringer et al. Studies of terrain surface roughness and its effect on GNSS-R systems using airborne lidar measurements
Lamer et al. Evaluation of gridded scanning ARM cloud radar reflectivity observations and vertical doppler velocity retrievals
Lovell et al. Foliage profiles from ground based waveform and discrete point lidar
Chaubell et al. Backus-Gilbert optimal interpoaltion applied to enhance smap data: Implementation and assessment
Lee Forest parameter estimation using polarimetric SAR interferometry techniques at low frequencies
Bennett et al. The UK NERC fully portable polarimetric ground-based synthetic aperture radar (GB-SAR)
Pan et al. Airborne MMW InSAR interferometry based on time varying baseline and BP algorithm
CN117665809A (en) Method for inverting forest canopy height
Ulander et al. Development and Initial Evaluation of the Tower-based Borealscat-2 P-and L-Band Tomographic SAR
Lamer et al. Evaluation of gridded Scanning ARM Cloud Radar reflectivity observations and vertical Doppler velocity retrievals.
CN117665809B (en) Method for inverting forest canopy height
Liu et al. Simulation of wet atmospheric delay correction for interferometric imaging altimeter based on radiometer
Albinet et al. High-resolution vertical polarimetric imaging of pine forests
Currie et al. High resolution 3-D radar imaging

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant