A forest canopy density inversion method based on geometric optical model

Abstract

The invention discloses a forest canopy closure inversion method based on a geometric optical model, comprising the steps of: 1) calculating the NDVI value based on the real surface reflectance of the image, extracting the area with the NDVI value between 0-0.05 as the pure bare ground (lighting background) pixel, and the area between 0.95-1 as the pure vegetation cover (lighting canopy) pixel, obtaining the shadow area pixel by visual interpretation, and then obtaining the average spectrum curve of the pixel sample by N-dimensional visualization analysis, and obtaining the end member information corresponding to the light background, light canopy and shadow component of the mixed pixel; 2) obtaining the light background end member abundance of each mixed pixel by using a fully constrained mixed pixel decomposition model; 3) simplifying the light background end member abundance input into a three-component Li-Strahler geometric optical model to invert the forest canopy closure. The invention can improve the accuracy of inverting forest canopy closure from medium-resolution multispectral images.

Description

Technical Field

[0001] The invention relates to a forest canopy closure inversion method, in particular to a forest canopy closure inversion method based on a geometric optical model. Background Art

[0002] There are many methods for extracting forest closure based on different remote sensing data sources. It is generally believed that the methods for inverting forest closure using remote sensing data can be summarized into the following three categories: One is the method based on empirical statistical models (Franklin J. and Strahler AH, 1988; Gao Yunfei, 2012). This method mainly takes forest closure as the dependent variable, selects several factors (such as band DN value, terrain factors, vegetation index, etc.) as independent variables, analyzes the correlation of independent variables, and establishes expressions based on the significance of each factor on the closure, so that the forest closure parameters can be estimated. Commonly used statistical methods include: multivariate stepwise regression (MSR), principal component regression analysis (PCR), partial least squares regression (PLSR), etc. (LEI et al., 2012; Du Xiaoming, 2008; Hu Zhenhua et al., 2017). The other is a semi-empirical model based on data mining, including remote sensing image classification and artificial neural network method, random forest method, support vector machine and other methods (Atzberger.C, 2004; Yang Yanting, 2018). The steps of this type of method generally include determining training samples, training models and estimating canopy density (Siavash Kalbi et al., 2018; Liu Saisai et al., 2020; Wang Rui et al., 2015). The third is to establish a physically meaningful mechanism model at the sub-pixel scale through theoretical analysis and decomposition of mixed pixels, thereby realizing the inversion of canopy density parameters (Chopping M.et al., 2006, 2008; Wang Cong, 2015).

[0003] Among the above-mentioned techniques for estimating forest closure, the empirical statistical model is the simplest and has a high inversion accuracy for a specific study area, but it is less applicable to other areas outside the study area. It is the most commonly used method for estimating forest closure. The semi-empirical model is flexible in operation and has high accuracy. It is one of the effective means of estimating forest closure. The mixed pixel decomposition model has a clear physical meaning and has great potential in the study of forest closure inversion. The use of the mixed pixel decomposition model is theoretically more accurate than the traditional remote sensing classification method (Lv Changchun, 2003), but the decomposition accuracy is still affected by factors such as vegetation and terrain. Although most experimental results have high accuracy, how to extract pure end members, reduce terrain errors, and weaken the influence of vegetation and bare soil needs to be further explored in future studies.

[0004] In remote sensing images, the ground object information contained in each pixel is the mixed spectral information of all ground object types in the pixel. Usually, multiple ground object types can be covered within the pixel scale, which forms the "mixed pixel problem". The mixed pixel problem is a common problem in remote sensing. It is one of the main reasons why ground objects are difficult to identify and have low classification accuracy. It is also the main obstacle to the development of quantitative remote sensing (Yang Wei, 2008). The mixed pixel decomposition process can be understood as the process of decomposing each pixel in the remote sensing image and calculating the percentage (abundance) of each type of ground object (end member) in the pixel coverage area (Lv Changchun, 2003). In order to better decompose the ground object types in mixed pixels, a series of pixel decomposition methods have been proposed under the continuous exploration of discipline pioneers, mainly including linear models, geometric optical models, probability models, random geometric models and fuzzy models (Zhao Yingshi, 2003).

[0005] There are multiple scattering phenomena between the end members of the ground objects, so the mixed pixel decomposition model can be divided into linear model and nonlinear model according to the number of scattering. The linear model does not consider the scattering effect between the ground objects. In the forest canopy estimation model, the pixel binary model is the most common linear model. The pixel binary model simply decomposes the pixel into canopy and non-canopy parts, and extracts the percentage of the canopy coverage of the tree layer in the image to estimate the forest canopy value of the pixel. At present, there are many studies on the inversion of forest canopy using the pixel binary model. Liu Dandan et al. (2018) used the pixel binary model to conduct a remote sensing inversion experiment of forest canopy on Landsat TM images, extracted the canopy coverage and non-canopy coverage information of the study area, established a pixel binary model, and calculated the normalized vegetation index (NDVI) of the study area as a parameter for forest canopy estimation. The experimental results have an accuracy of 89.01%, indicating that the pixel binary model can be effectively applied to forest resource monitoring in large areas. Xu Ding et al. (2013) conducted an experiment to extract canopy density using a binary pixel model and analyzed that the shrub and grass vegetation growing in the gaps between tree canopies in the study area had a significant impact on the model's estimation of forest canopy density. It was necessary to consider methods to weaken the impact of shrub and grass vegetation, or to choose a more complex mixed pixel decomposition model.

[0006] The most representative nonlinear mixed pixel decomposition model is the Li-Strahler geometric optical model. The Li-Strahler geometric optical model was proposed by Academician Li Xiaowen in 1985. The model is based on the assumption of bidirectional reflectance distribution function (Li and Strahler, 1985; 1986; 1992), regards the tree crown as an ellipsoid, and believes that the remote sensing information of mixed pixels can be divided into four components (Zeng Yuan et al., 2008; Fan Weiliang et al., 2010): illuminated tree crown, illuminated background, shadowed tree crown and shadowed background. Then, by analyzing the physical meaning (Gutman G and Ignatov A, 1998), a geometric optical model is established to estimate the forest crown parameters. Due to the clear physical mechanism of the geometric optical model, it has become a hot topic for domestic and foreign scholars (Gu Chengyan, 2018; Wang Cong, 2015). Summary of the invention

[0007] Aiming at the problem that the precision of mixed pixel decomposition using medium and high resolution multispectral images is not high, the present invention provides an effective mixed pixel decomposition method, and on this basis, a physical model is applied to invert forest canopy density parameters, including the following steps:

[0008] 1) Calculate the NDVI value based on the real surface reflectance of the image, extract the area with NDVI value between 0-0.05 as pure bare ground (light background) pixels, and the area with NDVI value between 0.95-1 as pure vegetation cover (light canopy) pixels, obtain the shadow area pixels through visual interpretation, and then obtain the average spectral curve of the pixel samples through N-dimensional visualization analysis to obtain the end member information corresponding to the light background, light canopy and shadow components of the mixed pixel;

[0009] 2) Using the fully constrained mixed pixel decomposition model to obtain the illumination background endmember abundance of each mixed pixel;

[0010] 3) The background endmember abundance input is simplified to a three-component Li-Strahler geometric optical model to invert forest canopy density.

[0011] The present invention adopts the method of directly selecting pixel sample spectral curves based on the NDVI value of the sample area and uses a fully constrained linear spectral decomposition model to perform pixel unmixing, thereby improving the accuracy of mixed pixel decomposition and the accuracy of forest canopy density inversion from medium-resolution multispectral images. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 It is a technical flow chart of the forest canopy density inversion method based on the geometric optical model;

[0013] Figure 2The invention is based on the Li-Strahler geometric optics model and the hybrid pixel decomposition technology to obtain the accuracy verification of the forest canopy density inversion results of Nanning in 2015. DETAILED DESCRIPTION

[0014] Figure 1 This is a technical flow chart of the forest canopy closure inversion method based on the geometric optical model. Based on this process, the inversion of the forest canopy closure of Nanning City using the multispectral data of the Landsat satellite can include the following three steps:

[0015] 1) Calculate the NDVI value based on the real surface reflectance of the image, extract the area with NDVI value between 0-0.05 as pure bare ground (light background) pixels, and the area between 0.95-1 as pure vegetation cover (light crown) pixels, obtain the shadow area pixels through visual interpretation, and then obtain the average spectral curve of the pixel samples through N-dimensional visualization analysis to obtain the end member information corresponding to the light background, light crown and shadow components of the mixed pixel;

[0016] 2) Using the fully constrained mixed pixel decomposition model to obtain the illumination background endmember abundance of each mixed pixel;

[0017] 3) The background endmember abundance input is simplified to a three-component Li-Strahler geometric optical model to invert forest canopy density.

[0018] The forest canopy density inversion process is described as follows:

[0019] 1) Data preparation and endmember extraction: Satellite remote sensing data is used to obtain the three endmember information of mixed pixels, namely the endmember information of illumination background, illumination canopy, and shadow component.

[0020] The medium-resolution image data of the sample area was obtained from the Geospatial Data Cloud website (http: / / www.gscloud.cn / ), and the real surface reflectance image of the sample area was obtained through preprocessing (terrain correction, radiation calibration, atmospheric correction, etc.).

[0021] Calculate the NDVI image of the sample area and obtain the average spectral curve of the illuminated background, illuminated tree crown, and shadow pixel samples. That is, set the NDVI threshold, extract the area with NDVI values ​​between 0 and 0.05 as pure bare ground (illuminated background) pixels, and the area between 0.95 and 1 as pure vegetation cover (illuminated tree crown) pixels, obtain the shadow area pixels through visual interpretation, and extract the average spectral curve of the sample pixels through N-dimensional visualization analysis.

[0022] 2) Extraction of illumination background endmember abundance: The illumination background endmember abundance of each mixed pixel is obtained using the fully constrained mixed pixel decomposition model.

[0023] The fully constrained mixed pixel decomposition model can decompose the area percentage (abundance) occupied by different end members in each pixel according to certain constraints. The specific calculation is as follows:

[0024]

[0025] Among them, ρ j is the pixel value of the jth band; k is the total number of end members; ρ ij is the reflectivity of the ith end member in the jth band; s i is the abundance of the ith end member; e j is the residual corresponding to the j-th band.

[0026] 3) Forest canopy closure inversion: The forest canopy closure is inverted by simplifying the background end-member abundance input into a three-component Li-Strahler geometric optical model.

[0027] The Li-Strahler model assumes that the tree crown is elliptical, analyzes the geometric structure of the ground object and decomposes the pixel into four end members according to the bidirectional reflectance distribution function (BRDF) (Zeng Yuan et al., 2009): illuminated crown (C), illuminated background (G), shadow crown (T) and shadow background (Z). The model describes the surface reflectance on the remote sensing image as the area percentage weighted sum of the four end members, calculated as follows:

[0028] R=K g G+K c C+K t T+K s S

[0029] Among them, R is the image pixel reflectivity, K g , K c , K t , K s are the endmember abundances of the illuminated background, illuminated canopy, shaded canopy, and shaded background, respectively. Further, based on the Li-Strahler geometric optical model, the mixed pixel can be simplified as a three-component model, that is, the mixed pixel is considered to be mainly composed of three components: the illuminated background (C), the illuminated canopy (G), and the shade (S). Then, equation (2) can be simplified to the following form:

[0030] R=K g G+K c C+(1-K c -K g )S

[0031] Based on the above three-component model, the canopy cover can be inverted by extracting the illumination background component, that is, the canopy closure. First, the area percentage of the illumination background component in the mixed pixel is determined according to the Boolean model:

[0032]

[0033] Among them, θ i ,θ v represents the zenith angle of the sun and the observation position, φ is the relative azimuth angle between the sun and the observation position, and O(θ i ,θ v , φ) is the overlapping part of the illumination shadow and the observation shadow, called the overlap function, and its calculation formula is as follows:

[0034] O(θ i ,θ v ,φ)=(t-sin t cos t)(secθ i +secθ v ) / π

[0035] The variable t can be obtained as follows:

[0036]

[0037] Among them, h and r represent the average crown height and average crown radius in the pixel, respectively. These two variables can be obtained from the measured forest structure parameters. The parameter m in the area percentage calculation formula of the light background component is the mean value of the tree cover index, which is an important parameter related to the light background component and forest canopy density. It can be obtained by reverse calculation from the following formula:

[0038]

[0039] Finally, according to the Poisson distribution theory, the forest canopy closure parameter (CC) can be calculated based on the parameter m:

[0040] CC=1-e -πm

[0041] The present invention is described below with a specific example:

[0042] Take the forest canopy density inversion of Nanning, Guangxi Zhuang Autonomous Region as an example. We collected Landsat TM / OLI data covering Nanning in 2015, and completed image preprocessing and end member extraction according to the specific step 1 of the forest canopy density inversion process. Then, based on the Li-Strahler geometric optical model and hybrid pixel decomposition technology of the present invention, the forest canopy density distribution map of Nanning in 2015 was obtained.

[0043] In order to evaluate the inversion effect of forest canopy density parameters, we calculated the R between the inversion results and the measured data of forest canopy density in the plain and mountainous areas of Nanning. 2 (coefficient of determination) and precision, the results are as follows Figure 2 The results for the plain area of ​​Nanning are: R 2 is 0.83, with an accuracy of 93.08%; while the results in the mountainous areas of Nanning are: R 2 The above accuracy verification results show that the forest canopy density inverted by the method of the present invention has a good effect.

Claims

1. A forest canopy density inversion method based on a geometric optical model, characterized in that: Includes steps: 1) Calculate the NDVI value based on the real surface reflectance of the image, extract the area with NDVI value between 0-0.05 as the pure bare ground illumination background pixel, and the area with NDVI value between 0.95-1 as the pure vegetation cover illumination crown pixel, obtain the shadow area pixel by visual interpretation, and then obtain the average spectrum curve of the pixel sample by N-dimensional visualization analysis, and obtain the end member information corresponding to the illumination background, illumination crown and shadow component of the mixed pixel; 2) Using the fully constrained mixed pixel decomposition model to obtain the illumination background endmember abundance of each mixed pixel; 3) Simplify the background endmember abundance input into three components of the Li-Strahler geometric optical model to invert forest canopy density; The step 1) comprises: The medium-resolution image data of the sample area was obtained from the geospatial data cloud website, and the real surface reflectance image of the sample area was obtained through terrain correction, radiation calibration, and atmospheric correction preprocessing; Calculate the NDVI image of the sample area and obtain the average spectral curve of the illumination background, illumination crown, and shadow pixel samples, that is, set the NDVI threshold, extract the area with NDVI values ​​between 0-0.05 as the pure bare ground illumination background pixel, and the area between 0.95-1 as the pure vegetation cover illumination crown pixel, obtain the shadow area pixel through visual interpretation, and then extract the average spectral curve of the sample pixel through N-dimensional visualization analysis; The step 2) comprises: Based on the average spectral curve of the three extracted endmembers, the illumination background endmember abundance K of each mixed pixel in the image is obtained using the fully constrained mixed pixel decomposition model. g 、Light crown K c and the shadow end member abundance K s , The fully constrained mixed pixel decomposition model can decompose the area percentage occupied by different end members in each pixel according to certain constraints, that is, the abundance. The specific calculation is as follows: Among them, ρ j is the pixel value of the jth band; k is the total number of end members; ρ ij is the reflectivity of the ith end member in the jth band; s i is the abundance of the ith end member; e j is the residual corresponding to the j-th band; The step 3) comprises: The obtained illumination background end member abundance K g The average tree height, average crown width, solar zenith angle, solar azimuth angle, observation zenith angle and observation azimuth angle parameters obtained by field measurements were input into the simplified Li-Strahler geometric optical model to obtain the forest canopy density of the sample area.

2. The forest canopy density inversion method based on geometric optical model according to claim 1 is characterized in that: Based on the Li-Strahler geometric optical model, the mixed pixel is simplified as a three-component model, that is, the mixed pixel R is mainly composed of three components: illumination background C, illumination canopy G, and shadow S, as shown in the following formula: R=K g G+K c C+(1-K c -K g )S。 3. The forest canopy density inversion method based on geometric optical model according to claim 2 is characterized in that: The canopy coverage is inverted by extracting the illumination background component, that is, the canopy closure, and the area percentage of the illumination background component in the mixed pixel is determined according to the Boolean model: Among them, m is the mean value of the tree cover index, θ i ,θ v represents the zenith angle of the sun and the observation position, φ is the relative azimuth angle between the sun and the observation position, and O(θ i ,θ v , φ) is the overlapping part of the illumination shadow and the observation shadow, which is called the overlap function.

4. The forest canopy density inversion method based on geometric optical model according to claim 3 is characterized in that: According to the Poisson distribution theory, the forest canopy closure parameter CC is calculated based on the parameter m: CC=1-e -πm 。

Images