CN106909750B - A kind of computational methods of broad-leaved Vegetation canopy reflectivity - Google Patents

Abstract

The invention discloses a calculation method and model of broad-leaved vegetation canopy reflectance. The calculation method includes the following steps: S1: input parameters, identify and classify the input parameters, and divide them into leaf parameters, canopy parameters and soil parameters; S2 : calculate the reflectance and transmittance of a single leaf according to the leaf parameters; S3: obtain the extinction coefficient and scattering coefficient of the canopy according to the canopy parameters and the leaf parameters in S2; S4: obtain the canopy extinction and scattering parameters according to the obtained Obtain the relative reflection factor and reflectivity of the canopy; S5 obtain the canopy reflectivity according to the relative reflection factor and reflectivity of the canopy. The present invention couples the PROSPECT model with the SAIL model, cancels the input process of leaf reflectance and transmittance in the vegetation canopy reflectance simulation process under the condition of making full use of available parameters, and the present invention effectively simplifies the vegetation canopy spectral information simulation process The problem of parameter acquisition and algorithm optimization can speed up the calculation process, and the coupling model is beneficial to the inversion of vegetation parameters.

Description

A Calculation Method of Canopy Reflectance of Broad-leaved Vegetation

technical field

The invention relates to the evaluation of the health status of vegetation in road areas of high-grade highways, in particular to a calculation method and model for the reflectivity of broad-leaved vegetation canopies.

Background technique

Vegetation information inversion has always been the most promising research field of quantitative remote sensing, which provides a favorable basis for the development of vegetation ecological environment monitoring and evaluation research. Scholars at home and abroad have proposed many vegetation parameter inversion models, which can be mainly divided into two categories: statistical models and physical models. Among them, the physical model is favored for its excellent wide applicability and stability. The physical inversion model of vegetation information has a fairly strong physical basis and does not depend on the specific types of vegetation and changes in the background environment. However, the accuracy of the physical model in large-scale vegetation information inversion is still poor and the calculation is too large to meet the actual production and application requirements. This is mainly due to two reasons: first, the physical model has high requirements for quantitative remote sensing. It is necessary to invert the surface reflectance through remote sensing images, and at the same time, many difficult-to-obtain accurate parameters are required as input parameters of the physical model; The influence of factors such as soil background reflectance and canopy structure limits the accuracy of biochemical parameter inversion and increases the difficulty of inversion.

The most researched and applied models in recent years are the PROSPECT model for broadleaves, the LIBERTY model for coniferous leaves, and the canopy SAIL model. The PROSPECT model regards broad-leaved leaves as a composite flat plate model composed of N layers of rough flat plates containing leaf biochemical substances and N-1 layers of air. The reflectance and transmittance of this composite slab model. The LIBERTY model regards the cells of coniferous leaves as standard circular cells, and considers that the needles are formed by stacking countless leaf cells in the air, and further obtains the spectral information of infinite cells by obtaining the optical properties of individual cells. The canopy SAIL model is a four-flow approximation to the radiative transfer equation based on the radiative transfer theory. The model inputs the solar zenith angle, solar azimuth angle, observation zenith angle, observation azimuth angle, LAI, soil reflectance, hot spots and Parameters such as the reflectance and transmittance of a single leaf are used to simulate the reflectance of the canopy from 400nm to 2500nm. It is obvious that the output parameters of the broad-leaved and needle-leaved models can be used as the input parameters of the SAIL model to realize the coupling of the models and further realize the simulation of vegetation canopy information. Therefore, scholars at home and abroad have successively carried out research on model coupling and its application. The current application is to use the PROSPECT leaf model to couple with the SAIL canopy reflection model. This coupling model uses the corrected leaf model and canopy model to simulate vegetation canopy. Layer reflectivity, with higher inversion accuracy and faster inversion speed. However, most of the currently used PROSAIL models are couplings of relatively well-developed initial models. The coupling of the PROSPECT model and the SAIL model fails to consider the asymmetry of the upper and lower surfaces of vegetation leaves and the hot spot effect. Although the adaptability However, in the actual application process, problems such as insufficient inversion of vegetation parameters are still exposed, which makes the coupling model subject to certain program restrictions in the actual application process.

Therefore, it is necessary to design a new calculation method and model of broad-leaved vegetation canopy reflectance.

Contents of the invention

The technical problem solved by the present invention is to propose a broad-leaved vegetation canopy reflectance in view of the limitation of the physical model in the actual application process, and in order to solve the parameter acquisition and inversion accuracy problems of the vegetation physical model in the application process. The calculation method and model of the canopy model, the single-leaf reflectance and transmittance obtained by the single-leaf model using the obtained parameters are used as the input parameters of the canopy model, and other input parameters of the canopy model are used to simulate the reflectance of the vegetation canopy, and the vegetation after coupling The canopy model can omit the acquisition process of some complex parameters, and the coupled model is an integrated model, which is more conducive to the inversion of vegetation parameters.

Technical scheme of the present invention is:

A calculation method for broad-leaved vegetation canopy reflectance, comprising the following steps:

S1: parameter identification;

Input model parameters; and the input parameters are preliminarily divided into three categories: leaf parameters, soil parameters and canopy parameters; leaf parameters include pigment content, dry matter content, equivalent water thickness and leaf structure parameters; soil parameters are soil parameters The spectral reflectance of ; the canopy parameter is the remaining other parameters;

S2: Input the PROSPECT model according to the leaf parameters in step S1 to simulate the spectrum of a single leaf, and calculate the spectral reflectance and transmittance of a single leaf;

S3: Calculate the extinction coefficient and the scattering coefficient according to the soil parameters and canopy parameters in step S1;

S4: Input the extinction coefficient and scattering coefficient calculated in step S3 into the SAIL model, and calculate the relevant reflection factor and reflectance of the canopy, including the hemispherical reflectance in the direct sunlight direction, the hemispherical reflectance in the observation direction and the double direct reflectance;

S5: Calculate canopy reflectance.

The step S2 specifically includes the following steps:

S2.1: Calculate the absorption coefficient K; the absorption coefficient is a linear expression of the content of various chemical substances in the leaves, and the calculation formula is:

Among them, c(i) is the concentration of component i in the leaf, and component i is chemical substances such as pigment, water, and dry matter; k(i) is the specific absorption coefficient of component i; N is the structural parameter of the leaf, indicating that the The number of layers, N, can be calculated using the following empirical equation:

N=(0.9*SLA+0.025)/(SLA-0.1)

Among them, SLA is the specific leaf area, which refers to the leaf area per unit dry weight;

S2.2: Calculate the transmittance at the interface:

The transmittance at the interface includes two: one is the transmittance of light entering the blade from the air, and the other is the transmittance of light entering the air from the blade;

Natural light can be regarded as non-polarized light, fully polarized after reflection, and partially polarized after refraction after transmission. According to the Snell-Descarts law, the relationship between the incident angle, the refraction angle and the two refractive indices on the contact surface of the two media can be known, and then the light can be calculated according to the transmittance of the electromagnetic wave passing through the two media given by Schanda. The solid angle α is the average transmittance t _av (α,1,n) of the blade with the incident refractive index n of air with a refractive index of 1, let t ₁₂ =t _av (α,1,n);

According to Stern's research results, there is a relationship between the transmittance t ₂₁ and ^t _{12 of light} entering air with a refractive index of n from a leaf with a refractive index of n, and the interface is calculated accordingly. The transmittance at the interface; the calculation of the transmittance at the interface is an existing technology, and the existing software can be used to complete the calculation;

S2.3: Calculate the total transmittance τ ₁ of the light passing through the flat medium;

For natural light, the total transmittance _τ1 through the flat medium is the integral in the whole hemisphere 2π space, which is related to the absorption coefficient K and the blade thickness D. The calculation formula is:

Among them, t is an intermediate variable;

S2.4: Calculate the reflectance R _α (1) and transmittance T _α (1) of the single-layer flat plate:

Perform spectral simulation of a single-layer plate. The spectrum of a single-layer plate is the core of the spectral simulation of the entire leaf. It mainly includes calculating the reflectance and transmittance of a single-layer plate. The calculation formula is:

S2.5: Calculate the single leaf reflectance R _α (N) and transmittance T _α (N):

Carrying out leaf spectral information simulation, Stokes conducted an in-depth study of the optical phenomenon after light penetrates a finite-layer plate, and obtained a calculation theory for the overall reflectance and transmittance of light passing through a plate with N layers of the same reflectance and transmittance , at this time, the simulation formula of leaf spectral information can be obtained by applying its theory:

in:

The step S3 specifically includes the following steps:

S3.1 Shadow Compensation:

① Calculation of geometric factors related to extinction and scattering: firstly, weighted and discretized according to the distribution probability of general leaf inclination angles to obtain a set of discretized leaf inclination angles; then, calculate the extinction and scattering factors corresponding to each leaf inclination angle; for leaf inclination angle θ _l , the calculation method is as follows:

First, calculate the critical angles β _s and β _o , the calculation formula is:

Among them, is the angle between the normal direction of the horizontal plane and the normal direction of the blade, that is, the leaf inclination angle θ _s is the sun zenith angle, θ _o is the observation zenith angle;

When it is determined that the denominator in the calculation formula is not 0 and the calculation result of cosβ is less than 1, the values of the two critical angles are directly calculated; when the calculation result of cosβ is equal to 1, the two critical angles are equal to π; where cosβ generally refers to β cosine of _s and β _o ;

Then, calculate the extinction coefficient of a single leaf in the direct sunlight direction and the observation direction with the leaf inclination angle _θl , and the calculation formulas are respectively:

Among them, L'=lai/h, lai is the leaf area index, h is the canopy height;

② Calculation of auxiliary azimuth angles β ₁ , β ₂ , β ₃ : the calculation of auxiliary azimuth angles mainly depends on the relative unequal relationship between the two critical angles, and the value selection method is as follows:

If: beta ₁ beta ₂ beta ₃ ψ≤|β _s -β _o | ψ |β _s -β _o | 2π-β _s -β _o |β _s -β _o |<ψ<2π-β _s -β _o |β _s -β _o | ψ 2π-β _s -β _o ψ≥2π-β _s -β _o |β _s -β _o | 2π-β _s -β _o ψ

Among them, ψ is the relative azimuth between the sun direction and the observation direction, that is, the difference between the azimuth angles of the two directions;

After obtaining the three auxiliary azimuth angles, the multiplier related to the reflectance and transmittance of a single blade can be calculated, which is used to calculate the bidirectional scattering coefficient ω;

③ Calculate the two-way scattering coefficient ω(θ _l ):

After the auxiliary azimuth angle is obtained, the two-way scattering coefficient is calculated in conjunction with the single-blade reflectance and transmittance calculated in S2, and the calculation formula is:

Among them, ρ, τ, θ _l represent the reflectivity, transmittance and corresponding leaf inclination angle of the blade respectively; let ρ=R _α (N) in S2; τ=T _α (N).

S3.2 Calculate the scattering coefficient after adding the reflectance and transmittance of the blade:

① The calculation formulas for the backscattering coefficients of diffuse radiation E- and E+ are:

② The formula for calculating the forward scattering coefficient of diffuse radiation E- and E+ is:

③ The formula for calculating the backscattering coefficient of direct solar radiation E _S is:

④ The formula for calculating the forward scattering coefficient of direct solar radiation E _S is:

⑤ The calculation formula of the attenuation coefficient of diffuse radiation E- and E+ is:

⑥ The formula for calculating the backscatter coefficient of radiation E ₀ in the observation direction is:

⑦ The formula for calculating the forward scattering coefficient of radiation E ₀ in the observation direction is:

S3.3 Calculation of general leaf inclination distribution probability (leaf shape distribution elements):

Calculate the general leaf inclination angle distribution under the average leaf inclination angle distribution: the leaf inclination angle refers to the angle between the normal direction of the horizontal plane and the normal direction of the leaves; different types of vegetation canopy correspond to different vegetation leaf inclination angles, corresponding to two leaf distributions at the same time Parameters LIDFa and LIDFb.

In the case of non-ellipsoidal distributions, the general leaf inclination distribution probability can be calculated from the mean leaf inclination, and the calculation process is completed by the following iterations:

x=2θ

y=LIDFa sin(x)+0.5LIDFb sin(2x)

dx＝0.5(y-x+2θ)

x=x+dx

until |dx|<t;

Then F(θ)=2(y+θ)/π.

Among them, θ is the discrete value of the average leaf inclination, and F(θ) is the cumulative leaf inclination, that is, the general leaf inclination distribution probability.

For the ellipsoidal distribution, the parameter LIDFa represents the distribution angle, which is 30 degrees, and the parameter LIDFb is 0. When calculating the distribution probability of the general leaf inclination angle, it can be obtained by using the leaf inclination angle density function made by Campbell (1986);

S3.4 determine the model parameters of the entire canopy;

For any canopy, the leaf inclination will not be fixed, it is a continuous process of increasing from 0 to 90, so when determining the model parameters of the entire canopy, the calculation method is to calculate the general leaf inclination when the leaf inclination is θ _l The product of the inclination distribution probability F(θ _l ) and the model parameters is added, and the calculation formula is:

Z＝∑F(θ _l )Z(θ _l )

Among them, F(θ _l ) is the general leaf inclination angle distribution probability corresponding to the leaf inclination angle θ _l , Z(θ _l ) is the parameter of a single blade, and Z is the modified model of leaf inclination angle Z=∑F(θ _l )Z(θ _l ) to the parameters of the entire canopy after correction; Z(θ _l ) refers to k(θ _l ), K(θ _l ), ω(θ _l ), σ( θ _l ), σ′(θ _l ), s(θ _l ), s′(θ _l ), a(θ _l ), v(θ _l ), u(θ _l ); Z refers to Substitute any one of all coefficients k, K, w, σ, σ', s, s', a, v, u in the four differential equations of the SAIL model for the entire canopy.

Described step S4 specifically comprises the following steps:

S4.1: Calculate the leaf area parameter; the leaf area parameter is an important indicator for evaluating the albedo of the vegetation canopy, and its parameter calculation formula is:

τ _ss =e ^-klai

τ _oo = e ^-Klai

ρ _dd =(e ^mlai -e ^-mlai )/(h ₁ e ^mlai -h ₂ e ^-mlai )

τ _dd =(h ₁ -h ₂ )/(h ₁ e ^mlai -h ₂ e ^-mlai )

Among them, K and k represent the extinction coefficient of the observation direction radiation E ₀ and the direct solar radiation E _s respectively; lai represents the leaf area index, which is the input parameter of the model;

S4.2: Calculate the parameters after adding the hot spot effect, the calculation formula is as follows:

ρ _sd ＝C _S (1-τ _ss τ _dd )-D _S ρ _dd

τ _sd ＝D _S (τ _ss -τ _dd )-C _S τ _ss ρ _dd

ρ _do ＝C _o (1-τ _oo τ _dd )-D _o ρ _dd

τ _do ＝D _o (τ _oo -τ _dd )-C _o τ _oo ρ _dd

ρ _so ＝H _o (1-τ _ss τ _oo )-C _o τ _sd τ _oo -D _o ρ _sd

Among them, the calculation method of each intermediate variable is:

h ₁ =(a+m)/σ

h ₂ =(am)/σ=1/h ₁

C _S =[s'(ka)-sσ]/(k ² -m ² )

C _o =[v(Ka)-uσ]/(K ² -m ² )

D _s =[-s(k+a)-s'σ]/(k ² -m ² )

D _o =[-u(K+a)-vσ]/(K ² -m ² )

H _s ＝(uC _S +vD _s +w)/(K+k)

H _o =(sC _o +s'D _o +w)/(K+k)

Among them, τ _ssoo is the hot spot effect correction parameter; due to the problem of solar radiation, especially the temperature difference between the soil and vegetation in the sparse vegetation area is quite large, which is related to its physical surface and meteorological characteristics, so a temperature correction process is required, and it may also be Provides additional vegetation canopy information. The correction process of the hot spot effect is:

For hotspot correction according to the 2/(k+K) effect, first calculate the parameters according to the given hotspot value:

Given that the initial value of alf is alf=10 ⁶ ,

1) If the hotspot value q is greater than 0, it is considered that there is a hotspot effect. At this time Among them, θ _s is the solar zenith angle, θ _o is the observation zenith angle, is the azimuth of the sun; (if the calculation result of α is greater than 200, its alf value is 200); if the calculation result of alf is 0, the parameters τ _ssoo and sumint under the influence of pure hot spot effect without shadow are calculated according to the following formula:

τ _ssoo = τ _ss

If the calculation result of alf is not 0, it is considered that there is no hot spot effect, and τ _ssoo and sumint are calculated according to the method in 2);

2) If the hot spot value q is 0, it is considered that there is no hot spot effect, and the parameters τ _ssoo and sumint under the no hot spot effect are calculated by the method of circular addition: the specific steps are:

The first step, parameter initialization:

In the second step, τ _ssoo and sumint are calculated by the following iterations:

Step 1. Determine the value of i:

If i is less than 20, then let x2=-log(1-i*fint)/alf;

If i=20, let x2=1;

If i is greater than 20, then end the iteration, let τ _ssoo = f1;;

Step 2. Carry out the following calculations in sequence:

y2＝‐(K+k)*lai*x2+fhot*(‐exp(‐alf*x2))/alf;

f2=exp(y2);

sumint=sumint+(f2‐f1)*(x2‐x1)/(y2‐y1);

x1=x2;

y1=y2;

f1=f2;

Step 3, let i=i+1, and turn to step 1;

Among them, * means multiplication.

S4.3: Calculate bidirectional reflectance:

The two-way reflectance includes two parts, one part is single scattering effect with hot spot effect, and the other part is multiple scattering effect without hot spot effect. The sum of the two parts is the calculation result of the final reflectance effect on the top of the canopy:

ρ _so2 ＝ρ _so +w*lai*sumint

S4.4: Introduce soil reflection:

The soil is located at the bottom of the entire canopy reflection system. After the electromagnetic wave penetrates the canopy, it will hit the soil ground, and then reflect back to the top of the canopy to form a part of the canopy reflectance. The calculation scheme is:

d _n ＝1－r _s * _ρdd

S4.5: Calculate the bi-hemispheric albedo factor:

S4.6: Calculate the hemispheric reflectance in the direction of direct sunlight:

ρ _sd2 = ρ _dd +(τ _sd +τ _ss )*r _s *τ _dd /d _n

S4.7: Calculate the hemispherical reflectance in the observation direction:

ρ _do2 ＝ρ _do +(τ _do +τ _oo )*r _s *τ _dd /d _n

S4.8: Calculate the double direct reflectance:

ρ _sod2 ＝ρ _so +((τ _ss +τ _sd )*τ _do +(τ _sd +τ _ss *r _s *ρ _dd )*τ _oo )*r _s /d _n

ρ _so3 ＝ρ _so2 +ρ _so +τ _ssoo *r _s

where r _s is the soil reflectance.

Described step S5 specifically comprises the following steps:

First of all, according to the data in the database, the basic radiation amount of direct sunlight and atmospheric scattered light can be obtained, which are recorded as E _s and E _d respectively;

Then, calculate the coefficient skyl representing the proportion of scattered radiation in the atmosphere to the total radiation:

skyl＝0.847-1.61*sin(90-θ _s )+1.04*sin(90-θ _s ) ²

Finally, calculate the ratio of the radiation measured in the observation direction to the incident radiation to obtain the canopy reflectance, the formula is:

A calculation model of broad-leaved vegetation canopy reflectance, the formula is as follows:

Among them, R is the reflectivity of the canopy, E _s and E _d are the basic radiation amount of direct sunlight and atmospheric scattered light respectively, skyl is the coefficient indicating the proportion of scattered radiation in the atmosphere to the total radiation; ρ _do2 and ρ _so3 are the observed Directional hemispheric reflectance and direct reflectance; model parameters are solved according to the above calculation method of broad-leaved vegetation canopy reflectance.

Beneficial effect:

Based on the research results of the latest PROSPECT broad-leaf model and SAIL canopy model, the present invention couples the PROSPECT model for simulating broad-leaf leaf spectral information with the SAIL model for simulating canopy spectral information, making full use of available parameters to save leaves In the process of spectral simulation, the process of inputting leaf reflectance and transmittance in the process of vegetation canopy reflectance simulation is canceled, and the leaf model simulation algorithm is added to the canopy SAIL model to complete the coupling process, and the canopy spectrum is directly simulated by vegetation physical and chemical parameters. Coupling with the latest research results of the physical model can greatly improve the accuracy of the physical model simulation spectrum, improve the model's operational capability, and effectively simplify the parameter acquisition problem and optimize the algorithm in the simulation process of vegetation canopy spectral information, speed up the calculation process, and improve the vegetation performance. information retrieval capability.

Description of drawings

Fig. 1 is a schematic diagram of the present invention;

Fig. 2 is the model parameter that the present invention inputs;

Fig. 3 is experimental result of the present invention

Detailed ways

The present invention will be further described below in conjunction with the accompanying drawings.

The present invention adopts the improved SAILH canopy model based on the hot spot effect, firstly classifies the given parameters, and divides them into two parts: leaf parameters and canopy parameters, and then uses the physical and chemical parameters of the leaves to calculate the basic absorption of the leaves, and then The reflectance and transmittance of a single-layer flat plate are derived from the calculation of the transmittance under the action of no interface absorption, and finally the spectral information of the leaves under the given structural parameters is derived by using the layering theory; at the same time, for the canopy parameters, the given leaf inclination angle is first used The distribution model parameters are used to calculate the leaf type distribution elements, and then combined with the calculated geometric elements to calculate the extinction coefficient and scattering coefficient of the canopy model and combined with the leaf type distribution model, a complete set of canopy parameters under the complete leaf distribution is obtained through a cyclic process. Finally, The reflectivity of the canopy is determined by using the obtained two radiation parameters and the sky scattering coefficient.

In this implementation case, a tree in good growth condition and whose canopy spectrum can be easily measured by the highway is taken as an example to further illustrate the technical solution of the present invention. As shown in the accompanying drawings, the calculation method of the improved physical model PROSAILH involved in the present invention is as follows:

After determination, the parameters of the tree are as follows:

Chlorophyll content 45.80 Leaf inclination distribution type Inclined Equivalent water thickness 0.02 LAI 4.30 dry matter content 0.02 solar zenith angle 30.00 Blade structure parameters 1.30 observation zenith angle 10.00 hotspot 8.00 relative azimuth 0.00

Select the experimental wavelength as 400nm-2399nm

S1 parameter classification: After obtaining the relevant parameters of the vegetation, it is divided into two parts. Among them, the leaf parameters include the pigment content of the leaves, the equivalent water thickness, the dry matter content and the maximum incident angle (generally 59 degrees); the canopy parameters include the leaf inclination angle. Distribution parameters, leaf area index, solar zenith angle, observed zenith angle, relative azimuth, and background reflectance of hotspots and soils.

S2: Leaf spectral information calculation: After obtaining the corresponding parameters, it is first necessary to calculate the spectral reflectance and transmittance of a single leaf. The calculation scheme is as follows:

(1) Calculation of the absorption coefficient K: the absorption coefficient is a linear expression of the content of various chemical substances in the leaves, and the calculation formula is:

(2) Calculation of transmittance without interface absorption: natural light can be regarded as non-polarized light, fully polarized after reflection, and partially polarized by refracted light after transmission. According to the Snell-Descarts law, the relationship between the incident angle, the refraction angle and the two refractive indices on the contact surface of the two media can be known, and then the transmission rate of the electromagnetic wave passing through the two media given by Schanda can be calculated. The solid angle is determined by the average transmittance t _av (α,1,n) of the blade with the incident refractive index n of the air with a refractive index of 1. At this time, it is given the code t ₁₂ . According to the research results of Stern, the light from the refractive index n There is a relational formula t ₂₁ =n ⁻² t ₁₂ between the air transmittance t ₂₁ and t ₁₂ of the blade incident refractive index of 1, and the transmittance at the interface is calculated accordingly.

(3) Transmission of light in isotropic media: For natural light, the total transmittance τ through the flat medium is the integral in the 2π space of the entire hemisphere, which is related to the absorption coefficient K and the blade thickness D. The calculation formula is:

(4) Single-layer plate spectrum simulation: The spectrum of a single-layer plate is the core of the spectrum simulation of the entire leaf, mainly including the reflectance and transmittance of a single-layer plate, and its calculation formula is:

(5) Simulation of reflectance and transmittance of a single leaf: Stokes conducted an in-depth study of the optical phenomenon after light penetrates a finite-layer plate, and obtained the overall reflection of light passing through a plate with N layers of the same reflectance and transmittance rate and transmittance, at this time, the simulation formula of leaf spectral information can be obtained by applying its theory:

in:

S3: Calculation of extinction coefficient and scattering coefficient

(6) Calculation of geometric elements:

(7) Calculation of leaf shape distribution elements:

General leaf inclination angle distribution under average leaf inclination angle distribution: Different types of vegetation canopy correspond to different vegetation leaf inclination angles, corresponding to two leaf distribution parameters LIDFa and LIDFb at the same time. In the case of non-ellipsoidal distributions, the general inclination can be calculated from the average leaf inclination, which can be done in a simple iteration with the pseudocode:

x=2θ

y=LIDFa sin(x)+0.5LIDFb sin(2x)

dx＝0.5(y-x+2θ)

x=x+dx

until |dx|<t;

Then F(θ)=2(y+θ)/π

where θ is the discrete value of the mean leaf inclination and F(θ) is the cumulative leaf inclination.

For the ellipsoidal distribution, the parameter LIDFa represents the distribution angle, which is 30 degrees, and the parameter LIDFb is 0. When calculating the general angle distribution, it can be obtained by using the leaf inclination angle density function made by Campbell (1986).

(8) Calculation of extinction and scattering coefficient:

Shadow compensation:

① Calculation of geometric factors related to extinction and scattering: firstly, the leaf inclination distribution is weighted and discretized to obtain a set of discretized center angle value leaf inclination distributions. Afterwards, a set of extinction and scattering factors is calculated corresponding to each leaf inclination angle, and the single calculation scheme is as follows:

Firstly, the critical angles β _s and β _o are calculated, and the calculation formula is:

When it is determined that the denominator in the calculation formula is not 0 and the calculation result of cosβ is less than 1, the values of the two critical angles are directly calculated; when the calculation result of cosβ is greater than 1, the two critical angles are equal to π; where cosβ generally refers to β cosine of _s and β _o ;

After calculating the two temporary critical angles, the extinction coefficients of the direct sunlight direction and the observation direction can be calculated, and the calculation formulas are:

Among them, L'=lai/h, lai is the leaf area index, h is the canopy height;

②Calculation of auxiliary azimuth angles β ₁ , β ₂ , and β ₃ for calculating bidirectional scattering coefficient W: the calculation of auxiliary azimuth angles mainly depends on the relative inequality between the two critical angles, and the value method is as follows:

If: beta ₁ beta ₂ beta ₃ ψ≤|β _s -β _o | ψ |β _s -β _o | 2π-β _s -β _o |β _s -β _o |<ψ<2π-β _s -β _o |β _s -β _o | ψ 2π-β _s -β _o ψ≥2π-β _s -β _o |β _s -β _o | 2π-β _s -β _o ψ

Among them, ψ is the relative azimuth between the sun direction and the observation direction, that is, the difference between the azimuth angles of the two directions;

After obtaining the three auxiliary azimuth angles, the multiplier related to the reflectance and transmittance of a single blade can be calculated, which is used to calculate the bidirectional scattering coefficient W.

③Calculation of two-way scattering coefficient: After obtaining the auxiliary azimuth angle, the two-way scattering coefficient can be calculated in conjunction with the calculated single-blade reflectance and transmittance in S2, and the calculation formula is:

Among them, ρ, τ, θ _l represent the reflectivity, transmittance and corresponding leaf inclination angle of the blade respectively; let ρ=R _α (N) in S2; τ=T _α (N).

Model parameter solution

(9) Calculation of scattering coefficient (addition of blade reflectance and transmittance):

① The calculation formulas for the backscattering coefficients of diffuse radiation E- and E+ are:

② The formula for calculating the forward scattering coefficient of diffuse radiation E- and E+ is:

③ The formula for calculating the backscattering coefficient of direct solar radiation E _S is:

④ The formula for calculating the forward scattering coefficient of direct solar radiation E _S is:

⑤ The calculation formula of the attenuation coefficient of diffuse radiation E- and E+ is:

⑥ The formula for calculating the backscatter coefficient in the observation direction E ₀ is:

⑦ The formula for calculating the forward scattering coefficient in the observation direction E ₀ is:

For any canopy, the leaf inclination will not be fixed, it is a continuous process from 0 to 90, so when determining the model parameters of the entire canopy, the calculation method is the leaf inclination distribution function as the probability function F(θ ) and the model parameters are multiplied and then added, the calculation formula is:

Z＝∑F(θ _l )Z(θ _l )

Among them, F(θ _l ) is the general leaf inclination angle distribution probability corresponding to the leaf inclination angle θ _l , Z(θ _l ) is the parameter of a single blade, and Z is the modified model of leaf inclination angle Z=∑F(θ _l )Z(θ _l ) to the parameters of the entire canopy after correction; Z(θ _l ) refers to k(θ _l ), K(θ _l ), ω(θ _l ), σ(θ _l ), σ′(θ _l ), s(θ _l ), s′(θ _l ), a(θ _l ), v(θ _l ), u(θ _l ); Z refers to Any one of all coefficients k, K, w, σ, σ′, s, s′, a, v, u in the four differential equations of the SAIL model for the entire canopy.

S4: Calculate the relative reflectance factor and reflectance of the canopy, including the hemispherical reflectance in the direct sun direction, the hemispherical reflectance in the observation direction and the double direct reflectance;

(1) Addition of leaf area index: leaf area index is an important indicator for evaluating the albedo of vegetation canopy, and its parameter calculation formula is:

τ _ss =e ^-klai

τ _oo = e ^-Klai

ρ _dd =(e ^mlai -e ^-mlai )/(h ₁ e ^mlai -h ₂ e ^-mlai )

τ _dd =(h ₁ -h ₂ )/(h ₁ e ^mlai -h ₂ e ^-mlai )

(2) Calculation of the parameters of the solution: the parameters after adding the hot spot effect are more complicated, and the calculation formula is as follows:

ρ _sd ＝C _S (1-τ _ss τ _dd )-D _S ρ _dd

τ _sd ＝D _S (τ _ss -τ _dd )-C _S τ _ss ρ _dd

ρ _do ＝C _o (1-τ _oo τ _dd )-D _o ρ _dd

τ _do ＝D _o (τ _oo -τ _dd )-C _o τ _oo ρ _dd

ρ _so ＝H _o (1-τ _ss τ _oo )-C _o τ _sd τ _oo -D _o ρ _sd

The calculation method of each variable is as follows:

h ₁ =(a+m)/σ

h ₂ =(am)/σ=1/h ₁

C _S =[s'(ka)-sσ]/(k ² -m ² )

C _o =[v(Ka)-uσ]/(K ² -m ² )

D _s =[-s(k+a)-s'σ]/(k ² -m ² )

D _o =[-u(K+a)-vσ]/(K ² -m ² )

H _s ＝(uC _S +vD _s +w)/(K+k)

H _o =(sC _o +s'D _o +w)/(K+k)

(3) Correction of hotspot effect: Due to the problem of solar radiation, especially the temperature difference between soil and vegetation in sparsely vegetated areas is quite large, which is related to its physical surface and meteorological characteristics, so a temperature correction process is required, and it may also provide additional vegetation canopy information. The correction process of the hot spot effect is:

For hotspot correction according to the 2/(k+K) effect, first calculate the parameters according to the given hotspot value:

Given that the initial value of alf is alf=10 ⁶ ,

1) If the hotspot value q is greater than 0, it is considered that there is a hotspot effect. At this time Among them, θ _s is the solar zenith angle, θ _o is the observation zenith angle, is the azimuth of the sun; (if the calculation result of α is greater than 200, its alf value is 200); if the calculation result of alf is 0, the parameters τ _ssoo and sumint under the influence of pure hot spot effect without shadow are calculated, and the calculation formula is:

τ _ssoo = τ _ss

If the calculation result of alf is not 0, it is considered that there is no hot spot effect, and τ _ssoo and sumint are calculated according to the method in 2);

2) If the hot spot value q is 0, it is considered that there is no hot spot effect, and the parameters τ _ssoo and sumint under the no hot spot effect are calculated by the method of circular addition: the specific steps are:

The first step gives the initial parameters in the loop:

fhot=lai*sqrt(K*k);

x1=0;

y1=0;

f1=1;

fint=(-exp(-alf))*0.05;

sumint = 0;

The second step calculates the parameters:

Finally let τ _ssoo = f1;

The * in the above code means multiplication;

In the program calculation, in the first 19 cycles, let the intermediate parameter x2 calculate the result according to the given formula and substitute it into the following formula to calculate sumt; in the 20th cycle, set the value of the intermediate parameter x2 to 1 and substitute it into the following formula to calculate sumt sumint, let τ _{ssoo be} equal to f1 calculated in the 20th loop.

(4) Two-way reflectivity calculation: The two-way reflectivity includes two parts, one with hot spot effect and the other without hot spot effect. The sum of the two parts is the final calculation result:

ρ _so2 ＝ρ _so +w*lai*sum

(5) The introduction of soil reflection: the soil is located at the bottom of the entire canopy reflection system. After the electromagnetic wave penetrates the canopy, it will hit the soil ground and reflect back to the top of the canopy after reflection, forming a part of the canopy reflectance. The calculation scheme is: 1－r _s *ρ _dd

(6) Calculation of the reflectivity factor of the two hemispheres:

(7) Calculation of hemispherical reflectance in the direction of direct sunlight:

ρ _sd2 ＝ρ _dd +(τ _sd +τ _ss )*r _s *τ _dd /(1－rs*ρ _dd )

(8) Calculation of hemispherical reflectance in the observation direction:

ρ _do2 ＝ρ _do +(τ _do +τ _oo )*r _s *τ _dd /(1－rs*ρ _dd )

(9) Calculation of double direct reflectance:

ρ _sod2 ＝ρ _so +((τ _ss +τ _sd )*τ _do +(τ _sd +τ _ss *r _s *ρ _dd )*τ _oo )*r _s /(1－rs*ρ _dd )

ρ _so3 ＝ρ _so2 +ρ _so +τ _ssoo *r _s

S5: Calculate canopy reflectance:

According to the data in the database, the basic radiation amounts of direct sunlight and atmospheric scattered light can be obtained, which are E _S and E _d respectively, and then calculate the skyl parameter of the atmospheric scattering coefficient:

sskyl＝0.847-1.61*sin(90-θ _s )+1.04*sin(90-θ _s ) ²

The calculation method of the canopy reflectance is the ratio of the measured radiation amount in the observation direction to the incident radiation amount, and the formula is:

The calculated vegetation reflectance is shown in Figure 3.

The physical and chemical parameters of the vegetation in the highway area selected this time and the spectra used for comparison are all measured values, and the RMSE is used to determine the usability and accuracy of the model simulation results.

After testing, the root mean square error (RMSE) of the spectral information simulated by this model in this example is 0.15706, and the simulation results of the simulation algorithm are available and have high precision.

The described examples are some, not all examples of the present invention, and should not be construed as limiting the present invention. Based on the examples in the present invention, all other implementations obtained by persons of ordinary skill in the art without making innovative efforts belong to the protection scope of the present invention.

Claims (3)

1. a calculation method of broad-leaved vegetation canopy albedo, is characterized in that, comprises the following steps: S1: parameter identification; Input the model parameters; and initially divide the input parameters into three categories: leaf parameters, soil parameters and canopy parameters; S2: Input the PROSPECT model according to the leaf parameters in step S1 to simulate the spectrum of a single leaf, and calculate the spectral reflectance and transmittance of a single leaf; S3: Calculate the extinction coefficient and the scattering coefficient according to the soil parameters and canopy parameters in step S1, the spectral reflectance and transmittance of a single leaf obtained in step S2; S4: Input the extinction coefficient and the scattering coefficient calculated in step S3 into the SAIL model, and calculate the relevant reflection factor and reflectance of the canopy; S5: Calculate canopy reflectance; The step S2 specifically includes the following steps: S2.1: Calculate the absorption coefficient K: <mrow><mi>K</mi><mo>=</mo><mo>&amp;Sigma;</mo><mfrac><mrow><mi>c</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow><mo>*</mo><mi>k</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow><mi>N</mi></mfrac></mrow> Among them, c(i) is the concentration of component i in the leaf; k(i) is the specific absorption coefficient of component i; N is the leaf structure parameter, indicating the number of leaf layers; N is calculated by the following empirical equation: N=(0.9*SLA+0.025)/(SLA-0.1) Among them, SLA is the specific leaf area, which refers to the leaf area per unit dry weight; S2.2: Calculate the transmittance at the interface: Calculate the average transmittance t _av (α, 1, n) of the light from air with a refractive index of 1 incident on a blade with a refractive index n at a solid angle α, let t ₁₂ =t _av (α, 1, n); S2.3: Calculate the total transmittance τ ₁ of the light passing through the flat medium; <mrow><msub><mi>&amp;tau;</mi><mn>1</mn></msub><mo>=</mo><mrow><mo>(</mo><mn>1</mn><mo>-</mo><mi>K</mi><mi>D</mi><mo>)</mo></mrow><mi>exp</mi><mrow><mo>(</mo><mo>-</mo><mi>K</mi><mi>D</mi><mo>)</mo></mrow><mo>+</mo><msup><mrow><mo>(</mo><mi>K</mi><mi>D</mi><mo>)</mo></mrow><mn>2</mn></msup><msubsup><mo>&amp;Integral;</mo><mrow><mi>K</mi><mi>D</mi></mrow><mi>&amp;infin;</mi></msubsup><msup><mi>t</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><msup><mi>exp</mi><mrow><mo>-</mo><mi>t</mi></mrow></msup><mi>d</mi><mi>t</mi></mrow> Among them, t is the intermediate variable, K is the absorption coefficient, and D is the blade thickness; S2.4: Calculate the reflectance R _α (1) and transmittance T _α (1) of the single-layer flat plate: <mrow><msub><mi>R</mi><mi>&amp;alpha;</mi></msub><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow><mo>=</mo><mn>1</mn><mo>-</mo><msub><mi>t</mi><mn>12</mn></msub><mo>+</mo><mfrac><mrow><msup><msub><mi>&amp;tau;</mi><mn>1</mn></msub><mn>2</mn></msup><msubsup><mi>t</mi><mn>12</mn><mn>2</mn></msubsup><mrow><mo>(</mo><msup><mi>n</mi><mn>2</mn></msup><mo>-</mo><msub><mi>t</mi><mn>12</mn></msub><mo>)</mo></mrow></mrow><mrow><msup><mi>n</mi><mn>4</mn></msup><mo>-</mo><msup><msub><mi>&amp;tau;</mi><mn>1</mn></msub><mn>2</mn></msup><msup><mrow><mo>(</mo><msup><mi>n</mi><mn>2</mn></msup><mo>-</mo><msub><mi>t</mi><mn>12</mn></msub><mo>)</mo></mrow><mn>2</mn></msup></mrow></mfrac></mrow> <mrow><msub><mi>T</mi><mi>&amp;alpha;</mi></msub><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><msup><mi>n</mi><mn>2</mn></msup><msup><msub><mi>&amp;tau;</mi><mn>1</mn></msub><mn>2</mn></msup><msubsup><mi>t</mi><mn>12</mn><mn>2</mn></msubsup></mrow><mrow><msup><mi>n</mi><mn>4</mn></msup><mo>-</mo><msup><msub><mi>&amp;tau;</mi><mn>1</mn></msub><mn>2</mn></msup><msup><mrow><mo>(</mo><msup><mi>n</mi><mn>2</mn></msup><mo>-</mo><msub><mi>t</mi><mn>12</mn></msub><mo>)</mo></mrow><mn>2</mn></msup></mrow></mfrac></mrow> Among them, n is the refractive index of the blade; S2.5: Calculate the reflectance R _α (N) and transmittance T _α (N) of a single leaf, that is, the overall reflectance and transmittance after the light passes through N plates with the same reflectance and transmittance: <mfenced open = "" close = ""><mtable><mtr><mtd><mrow><msub><mi>R</mi><mi>&amp;alpha;</mi></msub><mrow><mo>(</mo><mi>N</mi><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>=</mo><mfrac><mrow><msub><mi>R</mi><mi>&amp;alpha;</mi></msub><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow><mrow><mo>(</mo><msup><mi>ab</mi><mrow><mi>N</mi><mo>-</mo><mn>1</mn></mrow></msup><mo>-</mo><msup><mi>a</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><msup><mi>b</mi><mrow><mn>1</mn><mo>-</mo><mi>N</mi></mrow></msup><mo>)</mo></mrow><mo>+</mo><mrow><mo>(</mo><msub><mi>T</mi><mi>&amp;alpha;</mi></msub><mo>(</mo><mn>1</mn><mo>)</mo><msub><mi>T</mi><mn>90</mn></msub><mo>(</mo><mn>1</mn><mo>)</mo><mo>-</mo><msub><mi>R</mi><mi>&amp;alpha;</mi></msub><mo>(</mo><mn>1</mn><mo>)</mo><msub><mi>R</mi><mn>90</mn></msub><mo>(</mo><mn>1</mn><mo>)</mo><mo>)</mo></mrow><mrow><mo>(</mo><msup><mi>b</mi><mrow><mi>N</mi><mo>-</mo><mn>1</mn></mrow></msup><mo>-</mo><msup><mi>b</mi><mrow><mn>1</mn><mo>-</mo><mi>N</mi></mrow></msup><mo>)</mo></mrow></mrow><mrow><msup><mi>ab</mi><mrow><mi>N</mi><mo>-</mo><mn>1</mn></mrow></msup><mo>-</mo><msup><mi>a</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><msup><mi>b</mi><mrow><mn>1</mn><mo>-</mo><mi>N</mi></mrow></msup><mo>-</mo><msub><mi>R</mi><mn>90</mn></msub><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow><mrow><mo>(</mo><msup><mi>b</mi><mrow><mi>N</mi><mo>-</mo><mn>1</mn></mrow></msup><mo>-</mo><msup><mi>b</mi><mrow><mn>1</mn><mo>-</mo><mi>N</mi></mrow></msup><mo>)</mo></mrow></mrow></mfrac></mrow></mtd></mtr></mtable></mfenced> <mrow><msub><mi>T</mi><mi>&amp;alpha;</mi></msub><mrow><mo>(</mo><mi>N</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><msub><mi>T</mi><mi>&amp;alpha;</mi></msub><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow><mrow><mo>(</mo><mi>a</mi><mo>-</mo><msup><mi>a</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><mo>)</mo></mrow></mrow><mrow><msup><mi>ab</mi><mrow><mi>N</mi><mo>-</mo><mn>1</mn></mrow></msup><mo>-</mo><msup><mi>a</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><msup><mi>b</mi><mrow><mn>1</mn><mo>-</mo><mi>N</mi></mrow></msup><mo>-</mo><msub><mi>R</mi><mn>90</mn></msub><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow><mrow><mo>(</mo><msup><mi>b</mi><mrow><mi>N</mi><mo>-</mo><mn>1</mn></mrow></msup><mo>-</mo><msup><mi>b</mi><mrow><mn>1</mn><mo>-</mo><mi>N</mi></mrow></msup><mo>)</mo></mrow></mrow></mfrac></mrow> in: <mrow><mi>&amp;delta;</mi><mo>=</mo><msqrt><mrow><msup><mrow><mo>(</mo><msub><mi>R</mi><mn>90</mn></msub><msup><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow><mn>2</mn></msup><mo>-</mo><msub><mi>T</mi><mn>90</mn></msub><msup><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow><mn>2</mn></msup><mo>-</mo><mn>1</mn><mo>)</mo></mrow><mn>2</mn></msup><mo>-</mo><mn>4</mn><msub><mi>T</mi><mn>90</mn></msub><msup><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow><mn>2</mn></msup></mrow></msqrt><mo>;</mo></mrow> The step S3 specifically includes the following steps: S3.1 Shadow Compensation: ① Calculate the geometric factors related to extinction and scattering: Firstly, it is weighted and discretized according to the distribution probability of general leaf inclination angles to obtain a set of discretized leaf inclination angles (that is, the angle between the normal direction of the horizontal plane and the normal direction of the leaves); after that, the extinction and scattering factors corresponding to each leaf inclination angle are calculated; for Leaf inclination angle θ _l , the calculation method is as follows: First, calculate the critical angles β _s and β _o , the calculation formula is: <mrow><msub><mi>&amp;beta;</mi><mi>s</mi></msub><mo>=</mo><mi>arccos</mi><mrow><mo>(</mo><mo>-</mo><mfrac><mrow><msub><mi>cos&amp;theta;</mi><mi>l</mi></msub><msub><mi>cos&amp;theta;</mi><mi>s</mi></msub></mrow><mrow><msub><mi>sin&amp;theta;</mi><mi>l</mi></msub><msub><mi>sin&amp;theta;</mi><mi>s</mi></msub></mrow></mfrac><mo>)</mo></mrow></mrow> <mrow><msub><mi>&amp;beta;</mi><mi>o</mi></msub><mo>=</mo><mi>arccos</mi><mrow><mo>(</mo><mo>-</mo><mfrac><mrow><msub><mi>cos&amp;theta;</mi><mi>l</mi></msub><msub><mi>cos&amp;theta;</mi><mi>o</mi></msub></mrow><mrow><msub><mi>sin&amp;theta;</mi><mi>l</mi></msub><msub><mi>sin&amp;theta;</mi><mi>o</mi></msub></mrow></mfrac><mo>)</mo></mrow></mrow> Among them, θ _s is the solar zenith angle, θ _o is the observation zenith angle; When it is determined that the denominator in the calculation formula is not 0 and the calculation result of cosβ is less than 1, the values of the two critical angles are directly calculated; when the calculation result of cosβ is equal to 1, the two critical angles are equal to π; where cosβ generally refers to β cosine of _s and β _o ; Then, calculate the extinction coefficient of a single leaf in the direct sunlight direction and the observation direction with the leaf inclination angle _θl , and the calculation formulas are respectively: <mrow><mi>k</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>=</mo><mfrac><mn>2</mn><mi>&amp;pi;</mi></mfrac><msup><mi>L</mi><mo>&amp;prime;</mo></msup><mo>&amp;lsqb;</mo><mrow><mo>(</mo><msub><mi>&amp;beta;</mi><mi>s</mi></msub><mo>-</mo><mfrac><mi>&amp;pi;</mi><mn>2</mn></mfrac><mo>)</mo></mrow><msub><mi>cos&amp;theta;</mi><mi>l</mi></msub><mo>+</mo><msub><mi>sin&amp;beta;</mi><mi>s</mi></msub><msub><mi>tan&amp;theta;</mi><mi>s</mi></msub><msub><mi>sin&amp;theta;</mi><mi>l</mi></msub><mo>&amp;rsqb;</mo></mrow> <mrow><mi>K</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>=</mo><mfrac><mn>2</mn><mi>&amp;pi;</mi></mfrac><msup><mi>L</mi><mo>&amp;prime;</mo></msup><mo>&amp;lsqb;</mo><mrow><mo>(</mo><msub><mi>&amp;beta;</mi><mi>o</mi></msub><mo>-</mo><mfrac><mi>&amp;pi;</mi><mn>2</mn></mfrac><mo>)</mo></mrow><msub><mi>cos&amp;theta;</mi><mi>l</mi></msub><mo>+</mo><msub><mi>sin&amp;beta;</mi><mi>o</mi></msub><msub><mi>tan&amp;theta;</mi><mi>o</mi></msub><msub><mi>sin&amp;theta;</mi><mi>l</mi></msub><mo>&amp;rsqb;</mo></mrow> Among them, L'=lai/h, lai is the leaf area index, h is the canopy height; ②Calculation of auxiliary azimuth angles β ₁ , β ₂ , β ₃ : the value method is as follows: If: β ₁ β β ₂ β β ₃ β ψ≤|β _s -β _o | ψ |β _s -β _o | 2π-β _s -β _o |β _s -β _o |<ψ<2π-β _s -β _o |β _s -β _o | ψ 2π-β _s -β _o ψ≥2π-β _s -β _o |β _s -β _o | 2π-β _s -β _o ψ Among them, ψψ is the relative azimuth between the sun direction and the observation direction, that is, the difference between the azimuth angles of the two directions; ③ Calculate the two-way scattering coefficient ω(θ _l ): After the auxiliary azimuth angle is obtained, the two-way scattering coefficient can be obtained by combining the single-blade reflectance R _α (N) and transmittance T _α (N) and transmittance calculated in S2. The calculation formula is: <mfenced open = "" close = ""><mtable><mtr><mtd><mrow><mi>&amp;omega;</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>=</mo><mfrac><msup><mi>L</mi><mo>&amp;prime;</mo></msup><mrow><mn>2</mn><mi>&amp;pi;</mi></mrow></mfrac><mo>{</mo><mo>&amp;lsqb;</mo><mi>&amp;pi;</mi><mi>&amp;rho;</mi><mo>-</mo><msub><mi>&amp;beta;</mi><mn>2</mn></msub><mrow><mo>(</mo><mi>&amp;tau;</mi><mo>+</mo><mi>p</mi><mo>)</mo></mrow><mo>&amp;rsqb;</mo><mrow><mo>(</mo><mn>2</mn><msup><mi>cos</mi><mn>2</mn></msup><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>+</mo><msup><mi>sin</mi><mn>2</mn></msup><msub><mi>&amp;theta;</mi><mi>l</mi></msub><msub><mi>tan&amp;theta;</mi><mi>s</mi></msub><msub><mi>tan&amp;theta;</mi><mi>o</mi></msub><mi>cos</mi><mi>&amp;Psi;</mi><mo>)</mo></mrow><mo>+</mo><mo>(</mo><mi>&amp;rho;</mi></mrow></mtd></mtr><mtr><mtd><mrow><mo>+</mo><mi>&amp;tau;</mi><mo>)</mo><msub><mi>sin&amp;beta;</mi><mn>2</mn></msub><mo>&amp;lsqb;</mo><mfrac><mrow><mn>2</mn><msup><mi>cos</mi><mn>2</mn></msup><msub><mi>&amp;theta;</mi><mi>l</mi></msub></mrow><mrow><msub><mi>cos&amp;beta;</mi><mi>s</mi></msub><msub><mi>cos&amp;beta;</mi><mi>o</mi></msub></mrow></mfrac><mo>+</mo><msub><mi>cos&amp;beta;</mi><mn>1</mn></msub><msub><mi>cos&amp;beta;</mi><mn>3</mn></msub><msup><mi>sin</mi><mn>2</mn></msup><msub><mi>&amp;theta;</mi><mi>l</mi></msub><msub><mi>tan&amp;theta;</mi><mi>s</mi></msub><msub><mi>tan&amp;theta;</mi><mi>o</mi></msub><mo>&amp;rsqb;</mo><mo>}</mo></mrow></mtd></mtr></mtable></mfenced> Among them, ρ, τ, θ _l represent the reflectivity, transmittance and corresponding leaf inclination angle of the blade respectively; ρ=R _α (N); τ=T _α (N); S3.2 Calculate the scattering coefficient after adding the reflectance and transmittance of the blade: ① The formula for calculating the backscatter coefficient of diffuse radiation E- and E+ is: <mrow><mi>&amp;sigma;</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>=</mo><msup><mi>L</mi><mo>&amp;prime;</mo></msup><mrow><mo>(</mo><mfrac><mrow><mi>&amp;tau;</mi><mo>+</mo><mi>&amp;rho;</mi></mrow><mn>2</mn></mfrac><mo>+</mo><mfrac><mrow><mi>&amp;rho;</mi><mo>-</mo><mi>&amp;tau;</mi></mrow><mn>2</mn></mfrac><msup><mi>cos</mi><mn>2</mn></msup><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow></mrow> ② The formula for calculating the forward scattering coefficient of diffuse radiation E- and E+ is: <mrow><msup><mi>&amp;sigma;</mi><mo>&amp;prime;</mo></msup><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>=</mo><msup><mi>L</mi><mo>&amp;prime;</mo></msup><mrow><mo>(</mo><mfrac><mrow><mi>&amp;tau;</mi><mo>+</mo><mi>&amp;rho;</mi></mrow><mn>2</mn></mfrac><mo>-</mo><mfrac><mrow><mi>&amp;rho;</mi><mo>-</mo><mi>&amp;tau;</mi></mrow><mn>2</mn></mfrac><msup><mi>cos</mi><mn>2</mn></msup><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow></mrow> ③ The formula for calculating the backscattering coefficient of direct solar radiation E _S is: <mrow><mi>s</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mi>&amp;rho;</mi><mo>+</mo><mi>&amp;tau;</mi></mrow><mn>2</mn></mfrac><mi>k</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>-</mo><mfrac><mrow><mi>&amp;rho;</mi><mo>-</mo><mi>&amp;tau;</mi></mrow><mn>2</mn></mfrac><msup><mi>L</mi><mo>&amp;prime;</mo></msup><msup><mi>cos</mi><mn>2</mn></msup><msub><mi>&amp;theta;</mi><mi>l</mi></msub></mrow> ④ The formula for calculating the forward scattering coefficient of direct solar radiation E _S is: <mrow><msup><mi>s</mi><mo>&amp;prime;</mo></msup><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mi>&amp;rho;</mi><mo>+</mo><mi>&amp;tau;</mi></mrow><mn>2</mn></mfrac><mi>k</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>+</mo><mfrac><mrow><mi>&amp;rho;</mi><mo>-</mo><mi>&amp;tau;</mi></mrow><mn>2</mn></mfrac><msup><mi>L</mi><mo>&amp;prime;</mo></msup><msup><mi>cos</mi><mn>2</mn></msup><msub><mi>&amp;theta;</mi><mi>l</mi></msub></mrow> ⑤ The calculation formula of the attenuation coefficient of diffuse radiation E- and E+ is: <mrow><mi>a</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>=</mo><msup><mi>L</mi><mo>&amp;prime;</mo></msup><mrow><mo>(</mo><mn>1</mn><mo>-</mo><mfrac><mrow><mi>&amp;rho;</mi><mo>+</mo><mi>&amp;tau;</mi></mrow><mn>2</mn></mfrac><mo>+</mo><mfrac><mrow><mi>&amp;rho;</mi><mo>-</mo><mi>&amp;tau;</mi></mrow><mn>2</mn></mfrac><msup><mi>cos</mi><mn>2</mn></msup><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow></mrow> ⑥ The formula for calculating the backscatter coefficient of radiation E ₀ in the observation direction is: <mrow><mi>v</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mi>&amp;rho;</mi><mo>+</mo><mi>&amp;tau;</mi></mrow><mn>2</mn></mfrac><mi>K</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>+</mo><mfrac><mrow><mi>&amp;rho;</mi><mo>-</mo><mi>&amp;tau;</mi></mrow><mn>2</mn></mfrac><msup><mi>L</mi><mo>&amp;prime;</mo></msup><msup><mi>cos</mi><mn>2</mn></msup><msub><mi>&amp;theta;</mi><mi>l</mi></msub></mrow> ⑦ The formula for calculating the forward scattering coefficient of radiation E ₀ in the observation direction is: <mrow><mi>u</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mi>&amp;rho;</mi><mo>+</mo><mi>&amp;tau;</mi></mrow><mn>2</mn></mfrac><mi>K</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>l</mi></msub><mo>)</mo></mrow><mo>-</mo><mfrac><mrow><mi>&amp;rho;</mi><mo>-</mo><mi>&amp;tau;</mi></mrow><mn>2</mn></mfrac><msup><mi>L</mi><mo>&amp;prime;</mo></msup><msup><mi>cos</mi><mn>2</mn></msup><msub><mi>&amp;theta;</mi><mi>l</mi></msub></mrow> S3.3 First, for each leaf inclination θ obtained by discretization in ①, calculate its corresponding general leaf inclination distribution probability F(θ): In the case of non-ellipsoidal distributions, the general leaf inclination distribution probability is calculated from the mean leaf inclination, and the calculation process is done by the following iterations: x=2θ y=LIDFa sin(x)+0.5LIDFb sin(2x) dx=0.5(y-x+2θ) x=x+dx until |dx|<t; Then F(θ)=2(y+θ)/π; Among them, LIDFa and LIDFb are the leaf distribution parameters; θ is the discrete value of the average leaf inclination, F(θ) is the cumulative leaf inclination, that is, the general leaf inclination distribution probability; For the ellipsoidal distribution, the general leaf inclination distribution probability F(θ) is obtained by using the Campbell leaf inclination density function; S3.4 determine the model parameters of the entire canopy; Multiply the leaf inclination angle distribution probability corresponding to each leaf inclination angle obtained by discretization in (1) by the model parameters and then add them together to obtain new model parameters. The calculation formula is: Z＝∑F(θ _l )Z(θ _l ) Among them, F(θ _l ) is the general leaf inclination angle distribution probability corresponding to leaf inclination angle θ _l , and Z(θ _l ) refers to k(θ _l ), K(θ _l ), ω(θ _l ), σ(θ _l ), σ′(θ _l ), s(θ _l ), s′(θ _l ), a(θ _l ), v(θ _l ), u(θ _l ) Any one of ; Z correspondingly refers to any one of all coefficients k, K, w, σ, σ′, s, s′, a, v, u in the four differential equations of the SAIL model for the entire canopy. 2. the computing method of broad-leaved vegetation canopy reflectance according to claim 1, is characterized in that, described step S4 specifically comprises the following steps: S4.1: Calculation of leaf area parameters: τ _ss =e ^-klai τ _oo = e ^-Klai ρ _dd =(e ^mlai -e ^-mlai )/(h ₁ e ^mlai -h ₂ e ^-mlai ) τ _dd =(h ₁ -h ₂ )/(h ₁ e ^mlai -h ₂ e ^-mlai ) Among them, lai represents the leaf area index, which is the input parameter of the model; S4.2: Calculate the parameters after adding the hot spot effect, the calculation formula is as follows: ρ _sd ＝C _S (1-τ _ss τ _dd )-D _S ρ _dd τ _sd ＝D _S (τ _ss -τ _dd )-C _S τ _ss ρ _dd ρ _do ＝C _o (1-τ _oo τ _dd )-D _o ρ _dd τ _do ＝D _o (τ _oo -τ _dd )-C _o τ _oo ρ _dd ρ _so ＝H _o (1-τ _ss τ _oo )-C _o τ _sd τ _oo -D _o ρ _sd Among them, the calculation method of each intermediate variable is: <mrow><mi>m</mi><mo>=</mo><msqrt><mrow><msup><mi>a</mi><mn>2</mn></msup><mo>-</mo><msup><mi>&amp;sigma;</mi><mn>2</mn></msup></mrow></msqrt></mrow> h ₁ =(a+m)/σ h ₂ =(am)/σ=1/h ₁ C _S =[s'(ka)-sσ]/(k ² -m ² ) C _o =[v(Ka)-uσ]/(K ² -m ² ) D _s =[-s(k+a)-s'σ]/(k ² -m ² ) D _o =[-u(K+a)-vσ]/(K ² -m ² ) H _s ＝(uC _S +vD _s +w)/(K+k) H _o =(sC _o +s'D _o +w)/(K+k) Among them, τ _ssoo and sumint are hot spot effect correction parameters; S4.3: Calculate bidirectional reflectance: The two-way reflectance includes two parts, one part is single scattering effect with hot spot effect, and the other part is multiple scattering effect without hot spot effect. The sum of the two parts is the calculation result of the final reflectance effect on the top of the canopy: ρ _so2 ＝ρ _so +w*lai*sumint S4.4: Introduce soil reflection: d _n ＝1－r _s *ρ _dd S4.5: Calculate the bi-hemispheric albedo factor: <mrow><msub><mi>&amp;rho;</mi><mrow><mi>dd</mi><mn>2</mn></mrow></msub><mo>=</mo><msub><mi>&amp;rho;</mi><mi>dd</mi></msub><mo>+</mo><mrow><mo>(</mo><msubsup><mi>&amp;tau;</mi><mi>dd</mi><mn>2</mn></msubsup><mo>*</mo><msub><mi>r</mi><mi>s</mi></msub><mo>)</mo></mrow><mo>/</mo><msub><mi>d</mi><mi>n</mi></msub></mrow> S4.6: Calculate the hemispheric reflectance in the direction of direct sunlight: ρ _sd2 = ρ _dd +(τ _sd +τ _ss )*r _s *τ _dd /d _n S4.7: Calculate the hemispherical reflectance in the observation direction: ρ _do2 ＝ρ _do +(τ _do +τ _oo )*r _s *τ _dd /d _n S4.8: Calculate the double direct reflectance: ρ _sod2 ＝ρ _so +((τ _ss +τ _sd )*τ _do +(τ _sd +τ _ss *r _s *ρ _dd )*τ _oo )*r _s /d _n ρ _so3 ＝ρ _so2 +ρ _so +τ _ssoo *r _s where r _s is the soil reflectance. 3. the computing method of broad-leaved vegetation canopy reflectivity according to claim 2, is characterized in that, described step S5 specifically comprises the following steps: First, according to the data in the database, the basic radiation amount of direct sunlight and atmospheric scattered light is obtained, which are recorded as E _s and E _d respectively; Then, calculate the coefficient skyl representing the proportion of scattered radiation in the atmosphere to the total radiation: skyl＝0.847-1.61*sin(90-θ _s )+1.04*sin(90-θ _s ) ² Finally, calculate the ratio of the radiation measured in the observation direction to the incident radiation to obtain the canopy reflectance, the formula is: <mrow><mi>R</mi><mo>=</mo><mfrac><mrow><msub><mi>&amp;rho;</mi><mrow><mi>d</mi><mi>o</mi><mn>2</mn></mrow></msub><mo>*</mo><mi>s</mi><mi>k</mi><mi>y</mi><mi>l</mi><mo>*</mo><msub><mi>E</mi><mi>d</mi></msub><mo>+</mo><msub><mi>&amp;rho;</mi><mrow><mi>s</mi><mi>o</mi><mn>3</mn></mrow></msub><mo>*</mo><mrow><mo>(</mo><mn>1</mn><mo>-</mo><mi>s</mi><mi>k</mi><mi>y</mi><mi>l</mi><mo>)</mo></mrow><mo>*</mo><msub><mi>E</mi><mi>s</mi></msub></mrow><mrow><mi>s</mi><mi>k</mi><mi>y</mi><mi>l</mi><mo>*</mo><msub><mi>E</mi><mi>d</mi></msub><mo>+</mo><mrow><mo>(</mo><mn>1</mn><mo>-</mo><mi>s</mi><mi>k</mi><mi>y</mi><mn>1</mn><mo>)</mo></mrow><mo>*</mo><msub><mi>E</mi><mi>s</mi></msub></mrow></mfrac><mo>.</mo></mrow>