Vegetation leaf parameter remote sensing inversion method and system based on numerical optimization algorithm

Abstract

This invention discloses a remote sensing inversion method and system for vegetation leaf parameters based on a numerical optimization algorithm, belonging to the field of remote sensing inversion technology. This method calculates vegetation index information from remote sensing observations based on surface reflectance data of various bands, constructs a cost function based on the root mean square error, and then uses a numerical optimization algorithm to solve for the minimum point of the cost function to achieve vegetation leaf parameter inversion. This method can dynamically introduce prior knowledge, adjust the dimensionality of free variables and the initial values ​​of other fixed parameters of the model, and combine the gradient value information of the cost function, which can alleviate ill-conditioned inversion problems to a certain extent. It overcomes the shortcomings of traditional lookup table algorithms where the inversion accuracy is affected by the variable step size setting, and can be quickly applied to different remote sensing satellite data. The method is simple to operate.

Description

Technical Field

[0001] This invention relates to the field of remote sensing inversion technology, and in particular to a method and system for remote sensing inversion of vegetation leaf parameters based on numerical optimization algorithms. Background Technology

[0002] Methods for retrieving vegetation leaf parameters using remote sensing can be broadly categorized into estimation methods based on statistical models and those based on physical models. The former establishes an empirical fitting relationship between measured target parameters and remote sensing observation parameters to achieve regional estimation of the target parameters. This method often relies on measured data and is difficult to apply widely. The latter considers the influence mechanism of each parameter on canopy reflectivity and is more universally applicable. Retrieving vegetation leaf parameters based on physical models typically involves two processes: forward modeling and backward inversion. For continuous grassland canopies, the PROSAIL coupled vegetation radiative transfer model is often used. During backward inversion, a cost function for the target parameters is first constructed using satellite observation data, and then a lookup table algorithm is often used to retrieve the target vegetation leaf parameters.

[0003] Unlike lookup table algorithms, which involve a backward optimization process, optimization algorithms are forward optimization processes that find the optimal value by iteratively solving for the minimum point of the cost function to obtain the target parameter value. Common numerical optimization algorithms include the steepest descent method, the conjugate gradient method, and POWEL algorithm. These algorithms often require solving for the gradient of the cost function. Since the cost function is a composite function of the remote sensing physical model, the gradient calculation process is very complex, but it can be achieved using numerical differentiation or automatic differentiation techniques. When using numerical optimization algorithms to invert vegetation leaf parameters, prior knowledge can be dynamically introduced as constraints, which has great application potential in vegetation leaf parameter inversion compared to traditional optimization algorithms. Summary of the Invention

[0004] In view of this, the purpose of this invention is to provide a method and system for remote sensing inversion of vegetation leaf parameters based on numerical optimization algorithms. This method overcomes the shortcomings of traditional lookup table algorithms, where the inversion accuracy is affected by the variable step size setting, and can be quickly applied to different remote sensing satellite data.

[0005] To achieve the above objectives, the present invention provides the following technical solution:

[0006] The method for remote sensing inversion of vegetation leaf parameters based on numerical optimization algorithms provided by this invention includes the following steps:

[0007] Step 1: Obtain surface reflectance data for each band and calculate vegetation index information from remote sensing observations. The vegetation index includes at least the Normalized Differential Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Normalized Differential Infrared Index (NDII).

[0008] Step 2: Construct a radiative transfer model and set the model parameters, which include model variables, model constants, and display variables;

[0009] Step 3: Construct a cost function RMSE based on the root mean square error for vegetation indices, and set the partial derivative of RMSE with respect to a certain variable and the gradient information ▽T of the cost function;

[0010] Step 4: Use a constrained optimization algorithm to find the minimum point X of the constructed cost function. * The moisture content of the grassland canopy combustibles was obtained.

[0011] Furthermore, the vegetation index in step 1 is calculated according to the following formula:

[0012]

[0013]

[0014]

[0015] Among them, NDVI o EVI o NDII o ρ is the observed vegetation index extracted from satellite data. Blue , ρ Red , ρ NIR , ρ SWIR2 These are the spectral reflectances of blue light, red light, near-infrared light, and short-wave infrared bands extracted from satellite data.

[0016] Furthermore, the model parameters in step 2 are calculated according to the following formula:

[0017]

[0018] ρ(λ)=PROSAIL(C ab (5)

[0019]

[0020]

[0021]

[0022] Where, Cab represents chlorophyll a and b content, EWT is the isohydric thickness, DMC is the dry matter weight, LAI is the leaf area index, psoil is the soil factor, FMC is the combustible water content, ρ(λ) is the spectral reflectance simulated by the PROSAIL model at a wavelength of λ nm, and NDVI is...m EVI m NDII m The vegetation index is simulated by the PROSAIL model.

[0023] Furthermore, the RMSE and gradient information of the cost function in step 3 are calculated according to the following formula:

[0024]

[0025]

[0026] Where T(X) is a cost function with respect to NDVI, EVI, and NDII, and based on the root mean square error. This represents the gradient information of the cost function.

[0027] Furthermore, in step 4, a constrained optimization algorithm is used to find the minimum point X of the constructed cost function. * The specific steps are as follows:

[0028] Step 4.1: Set the initial values ​​for Cab, EWT, FMC, LAI, and psoil;

[0029] Step 4.2: Construct ecological rules and calculate them according to formula (11) to dynamically update the range of values ​​for free variables;

[0030]

[0031] Furthermore, regarding the measured grassland canopy vegetation parameters:

[0032] LAI=0.0218×FMC-1.1024,RMSE=1.05

[0033] EWT=0.0061×ln(DMC)+0.0409,RMSE=0.0024 (12)

[0034] Step 4.3: Calculate X during the k-th iteration using the following formula. k At the initial iteration direction d k :

[0035]

[0036]

[0037] Where k represents the k-th iteration; d k X represents the process in the k-th iteration. k Initial iteration direction; β k-1 This represents the search direction correction factor at the k-th iteration; The cost function T in X represents the value of X. k gradient at;

[0038] Step 4.4: Update d based on the free variable boundary information. k Get X k The feasible directions at the location are shown in the following formula:

[0039]

[0040] in, For component x j The result of the kth iteration; For x j exist Feasible directions at the location; and For component x j The lower and upper boundaries that should be satisfied at the k-th time are dynamically updated according to ecological rules.

[0041] Step 4.5: Determine the numerical iteration search step size λ according to formula (15). k And obtain the output of the (k+1)th iteration according to formula (16), that is, X. k+1 :

[0042]

[0043] X k+1 =X k +λ k d k (16)

[0044] λ represents the search step size; λ k X represents the search step size in the k-th iteration; k u represents the free parameters obtained in the k-th iteration. k X represents the value at the k-th iteration. k The upper boundary that should be satisfied; l k X represents the value at the k-th iteration. k The lower bound that should be satisfied; λ0 represents the initial maximum search step size;

[0045] Step 4.6: When formula (17) is satisfied, stop the iteration and output X. ※ =X k +1;

[0046]

[0047] d k+1 X represents the value at the (k+1)th iteration. k The search direction; Indicates C at the (k+1)th iteration abThe search direction; This indicates the search direction of EWT in the (k+1)th iteration; This indicates the search direction of the FMC in the (k+1)th iteration; This indicates the search direction of LAI in the (k+1)th iteration; This indicates the search direction of psoil in the (k+1)th iteration; ε represents the allowable error limit.

[0048] Furthermore, the value range of ε is 10. -6 ~10 -5 .

[0049] The vegetation leaf parameter remote sensing inversion system based on numerical optimization algorithm provided by the present invention includes a memory, a processor, and a computer program stored in the memory and capable of running on the processor. When the processor executes the program, it implements the above-mentioned method.

[0050] The beneficial effects of this invention are as follows:

[0051] This invention provides a method and system for remote sensing inversion of vegetation leaf parameters based on numerical optimization algorithms, belonging to the field of remote sensing inversion technology. Targeting the MODIS MCD43A4 surface reflectance product, and combining it with the PROSAIL vegetation radiative transfer model, it extracts vegetation indices from satellite observations and model simulations, constructs a cost function based on the root mean square error, and then uses a numerical optimization algorithm to solve for the minimum point of the cost function to achieve vegetation leaf parameter inversion. This method is simple to operate and, compared with traditional inversion methods based on lookup table algorithms, can dynamically introduce prior knowledge, adjust the dimensionality of free variables and the initial values ​​of other fixed parameters of the model, and combine the gradient information of the cost function, which can alleviate ill-conditioned inversion problems to a certain extent. It overcomes the shortcomings of traditional lookup table algorithms, whose inversion accuracy is affected by the variable step size setting, and can be quickly applied to different remote sensing satellite data.

[0052] Other advantages, objectives, and features of the invention will be set forth in part in the description which follows, and in part will be apparent to those skilled in the art from the following examination, or may be learned from practice of the invention. The objectives and other advantages of the invention can be realized and obtained through the following description. Attached Figure Description

[0053] To make the objectives, technical solutions, and beneficial effects of this invention clearer, the following figures are provided for illustration:

[0054] Figure 1 This is a schematic diagram of the overall method flow of the present invention.

[0055] Figure 2This is a schematic diagram of the geographical location of the research plot in a specific embodiment of the present invention.

[0056] Figure 3 This is a scatter plot showing the parameters of vegetation leaves in a specific embodiment of the present invention.

[0057] Figure 4 The result of inverting the moisture content of combustibles based on a numerical optimization algorithm in a specific embodiment of the present invention is shown. Detailed Implementation

[0058] The present invention will be further described below with reference to the accompanying drawings and specific embodiments, so that those skilled in the art can better understand and implement the present invention. However, the embodiments described are not intended to limit the present invention.

[0059] like Figure 1 As shown in the figure, this embodiment of the remote sensing inversion method for vegetation leaf parameters based on numerical optimization algorithms includes the following steps:

[0060] Step 1: Data Preparation

[0061] The measured grassland equal weight thickness (EWT), dry matter weight (DMC), fuel moisture content (FMC), and leaf area index (LAI) data in China were obtained from a field survey conducted in the grasslands of the Qinghai Lake Basin, China, from July 28 to August 2, 2015. The satellite remote sensing data used in this embodiment is the MODIS MCD43A4 surface reflectance product. High-quality surface reflectance information of each band was extracted to produce satellite-observed vegetation indices, namely NDVI, EVI, and NDII, as shown in formulas (1) to (3).

[0062]

[0063]

[0064]

[0065] Among them, NDVI o EVI o NDII o ρ is the observed vegetation index extracted from satellite data. Blue , ρ Red , ρ NIR , ρ SWIR2 These are the spectral reflectances of blue light, red light, near-infrared light, and short-wave infrared bands extracted from satellite data, respectively.

[0066] Step 2: Select PROSAIL as the radiative transfer model for forward modeling. The PROSAIL model couples the vegetation radiative transfer model 4SAIL and the leaf optical model PROSPECT, and is suitable for forward simulation of canopy biaxial reflectance of horizontally uniform vegetation types (such as grasslands and crops). The input parameters of this model include three categories: vegetation structure parameters (leaf structure parameter N, chlorophyll a and b content C). ab Leaf moisture content (EWT) and dry matter weight (DMC), and brown pigment composition (C) bp and carotenoid content C ar The input parameters for the PROSAIL model are shown in Table 1, including leaf tilt angle distribution (LIDFa and LIDFB), geometric observation parameters (solar zenith angle tts, observed zenith angle tto, relative azimuth angle psi), and soil parameter psoil.

[0067] Table 1 Input parameters for the PROSAIL model

[0068]

[0069] Set the chlorophyll a and b content C ab (μg / cm 2 Equal moisture content thickness EWT (g / cm) 2 Dry matter weight (DMC) (g / cm³) 2 Leaf area index (LAI) and soil factor psoil are the free variables in the PROSAIL model, and the other parameters of the model are constants (N=2, C). bp =0, C ar =8μg / cm 2 hspot=0.01, LIDFa=60, LIDFa=0); Using formula (4), the combustible moisture content FMC (%) is used as an explicit variable to further calculate the simulated vegetation index, i.e.:

[0070]

[0071] ρ(λ)=PROSAIL(C ab (5)

[0072]

[0073]

[0074]

[0075] Where ρ(λ) is the spectral reflectance simulated by the PROSAIL model at a wavelength of λnm, and ρ' Blue ,ρ'Red ,ρ' NIR ,ρ' SWIR2 These are the spectral reflectances simulated by the model in the blue (469nm), red (645nm), near-infrared (855nm), and short-wave infrared (1640nm) bands; NDVI m EVI m NDII m The vegetation index is simulated by the PROSAIL model.

[0076] Step 3: Construct a cost function RMSE (Equation 9) based on the root mean square error for NDVI, EVI, and NDII. Use the automatic differentiation tool Tapenade to obtain the partial derivative of RMSE with respect to a certain variable, and calculate the gradient information of the cost function.

[0077]

[0078]

[0079] Step 4: Use a constrained optimization algorithm to find the minimum point X of the constructed cost function. * ;

[0080] Specifically, it includes 6 parts:

[0081] Step 4.1: C ab Initial settings for EWT, FMC, LAI, and psoil;

[0082] The range of initial values ​​for the free parameters is determined by measured data, and the mean of the measured data is selected as the initial value for numerical iteration.

[0083] Furthermore, C ab The initial values ​​are set to 40, EWT to 0.015, FMC to 150, LAI to 4.0, and psoil to 0.5.

[0084] Step 4.2: Use formula (11) to quantitatively describe the relationship characteristics between variables, characterize ecological rules, and thereby constrain the range of values ​​for free variables;

[0085]

[0086] Where x represents the model free variables (i.e., LAI, EWT); This represents the fitting result (fitted LAI or fitted EWT) obtained from the empirical fitting equation of LAI with respect to the measured FMC or EWT with respect to the DMC. Let f(y) represent the RMSE of the empirical fitting equation of LAI with respect to the measured FMC or EWT with respect to the DMC, and let f(y) represent the empirical fitting equation of LAI with respect to the measured FMC or EWT with respect to the DMC.

[0087] Furthermore, regarding the measured grassland canopy vegetation parameters (FMC, LAI, DMC, and EWT) data in China (see...), Figure 3 Empirical fitting equations were constructed for the measured LAI with respect to the measured FMC and for the EWT with respect to the DMC, as shown in the following equations:

[0088]

[0089] Wherein, LAI represents the measured LAI data; f(FMC) represents the empirical fitting equation of LAI with respect to the measured FMC; FMC represents the measured FMC data; EWT represents the measured EWT data; f(DMC) represents the empirical fitting equation of EWT with respect to the measured DMC; DMC represents the measured DMC data; RMSE represents the RMSE of the empirical fitting equation of LAI with respect to the measured FMC or EWT with respect to the DMC;

[0090] Combining formulas (11) to (12), we construct the constraints for numerical iteration, i.e., the ecological rules, as shown in the following equations:

[0091]

[0092] Step 4.3: Calculate X during the k-th iteration. k At the initial iteration direction d k It can be calculated using an unconstrained optimization algorithm; here, the FR conjugate gradient algorithm is used for calculation:

[0093]

[0094] Where, d k X represents the process in the k-th iteration. k Initial iteration direction; β k-1 This represents the search direction correction factor at the k-th iteration; The cost function T in X represents the value of X. k gradient at;

[0095] Step 4.4: Update d based on the free variable boundary information. k Get X k The feasible directions at the location are shown in the following formula:

[0096]

[0097] in, For component x j The result of the kth iteration; For x j exist Feasible directions at the location; and For component x j The lower and upper boundaries that should be satisfied at the k-th time can be dynamically updated according to ecological rules;

[0098] Step 4.5: Determine the numerical iteration search step size λ according to formula (16). k And obtain the output of the (k+1)th iteration according to formula (17), that is, X. k+1 :

[0099]

[0100] X k+1 =X k *λ k d k (17)

[0101] Where λ represents the search step size; λ k X represents the search step size in the k-th iteration; k u represents the free parameters obtained in the k-th iteration. k X represents the value at the k-th iteration. k The upper boundary that should be satisfied; l k X represents the value at the k-th iteration. k The lower boundary that should be satisfied; λ0 represents the initial maximum search step size; in the above, the value of λ0 is 1;

[0102] Step 4.6: When formula (18) is satisfied, stop the iteration and output X. * =X k+1 ;

[0103]

[0104] d k+1 X represents the value at the (k+1)th iteration. k The search direction; Indicates C at the (k+1)th iteration ab The search direction; This indicates the search direction of EWT in the (k+1)th iteration; This indicates the search direction of the FMC in the (k+1)th iteration; This indicates the search direction of LAI in the (k+1)th iteration; This indicates the search direction of psoil in the (k+1)th iteration; ε represents the allowable error limit; the value range of ε is 10. -6 ~10 -5 .

[0105] Based on the above inversion method and system, the minimum point of the constructed cost function is solved, that is, the moisture content of grassland canopy combustibles is inverted. For example... Figure 4 As shown, this is the inversion result of grassland canopy combustible moisture content based on numerical optimization algorithm in this embodiment. The solution for grassland canopy combustible moisture content refers to the minimum point of the cost function obtained by solving.

[0106] The above-described embodiments are merely preferred embodiments provided to fully illustrate the present invention, and the scope of protection of the present invention is not limited thereto. Equivalent substitutions or modifications made by those skilled in the art based on the present invention are all within the scope of protection of the present invention. The scope of protection of the present invention is defined by the claims.

Claims

1. A remote sensing inversion method for vegetation leaf parameters based on numerical optimization algorithms, characterized in that: Includes the following steps: Step 1: Obtain surface reflectance data for each band and calculate vegetation index information from remote sensing observations. The vegetation index includes at least the Normalized Differential Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Normalized Differential Infrared Index (NDII). Step 2: Construct a radiative transfer model and set the model parameters, which include model variables, model constants, and display variables; Step 3: Construct a cost function RMSE based on the root mean square error for vegetation indices, and set the partial derivative of RMSE with respect to a certain variable and the gradient information ▽T of the cost function; Step 4: Use a constrained optimization algorithm to find the minimum point of the constructed cost function. ; Obtain the moisture content of grassland canopy combustibles; In step 4, a constrained optimization algorithm is used to find the minimum point of the constructed cost function. The specific steps are as follows: Step 4.1: Set the initial values ​​for Cab, EWT, FMC, LAI, and psoil; Wherein, Cab is the content of chlorophyll a and b, EWT is the isohydric thickness, LAI is the leaf area index, psoil is the soil factor, and FMC is the combustible moisture content; Step 4.2: Construct ecological rules and calculate them according to formula (11) to dynamically update the range of values ​​for free variables; (11) in, These represent the model's free variables, namely LAI and EWT. This represents the fitting result obtained from the empirical fitting equation of LAI with respect to the measured FMC or EWT with respect to the DMC, i.e., the fitted LAI or the fitted EWT; This represents the RMSE of the empirical fitting equation of LAI with respect to the measured FMC or EWT with respect to the DMC. This represents the empirical fitting equation of LAI with respect to the measured FMC or EWT with respect to the DMC; Furthermore, regarding the measured grassland canopy vegetation parameters: (12) Where DMC is the dry matter weight; Step 4.3: Calculate the results of the k-th iteration using the following formula. Initial iteration direction : (13) in, This indicates the k-th iteration; In the k-th iteration In the initial iteration direction; This represents the search direction correction factor at the k-th iteration; ) indicates that the cost function T is in gradient at; Step 4.4: Update based on free variable boundary information , get The feasible directions at the location are shown in the following formula: (14) in, For components The result of the kth iteration; for exist Feasible directions at the location; and For components The lower and upper boundaries that should be satisfied at the k-th time are dynamically updated according to ecological rules. Step 4.5: Determine the numerical iteration search step size according to formula (15). And obtain the output of the (k+1)th iteration according to formula (16), that is : (15) (16) Indicates the search step size; This represents the search step size in the k-th iteration; This represents the free parameters obtained in the k-th iteration; Indicates the k-th iteration The upper boundary that should be satisfied; Indicates the k-th iteration The lower boundary that should be satisfied; This represents the initial maximum search step size; Step 4.6: When formula (17) is satisfied, stop the iteration and output X. ※ =X k +1; (17) Indicates the (k+1)th iteration The search direction; Indicates C at the (k+1)th iteration ab The search direction; This indicates the search direction of EWT in the (k+1)th iteration; This indicates the search direction of the FMC in the (k+1)th iteration; This indicates the search direction of LAI in the (k+1)th iteration; ε represents the search direction of psoil in the (k+1)th iteration; ε represents the allowable error limit.

2. The method for remote sensing inversion of vegetation leaf parameters based on numerical optimization algorithm as described in claim 1, characterized in that: The vegetation index in step 1 is calculated using the following formula: (1) (2) (3) in, , , For the observed vegetation indices extracted from satellite data, , , , These are the spectral reflectances of blue light, red light, near-infrared light, and short-wave infrared bands extracted from satellite data.

3. The method for remote sensing inversion of vegetation leaf parameters based on numerical optimization algorithm as described in claim 1, characterized in that: The model parameters in step 2 are calculated according to the following formula: (4) (5) (6) (7) (8) in, The spectral reflectance is calculated using the PROSAIL model at a wavelength of λnm. , , The vegetation index is simulated by the PROSAIL model.

4. The method for remote sensing inversion of vegetation leaf parameters based on numerical optimization algorithm as described in claim 1, characterized in that: The cost function RMSE and cost function gradient information in step 3 are calculated according to the following formulas: (9) (10) in, Let be the cost function with respect to NDVI, EVI, and NDII, and based on the root mean square error. This represents the gradient information of the cost function.

5. The method for remote sensing inversion of vegetation leaf parameters based on numerical optimization algorithm as described in claim 1, characterized in that: The range of values ​​for ε is: .

6. A remote sensing inversion system for vegetation leaf parameters based on numerical optimization algorithms, comprising a memory, a processor, and a computer program stored in the memory and capable of running on the processor, characterized in that, When the processor executes the program, it implements the method described in claims 1-5.

Images