A method for calculating bidirectional reflectance factor of forest canopy high-resolution remote sensing image

Abstract

The application discloses a bidirectional reflectance factor calculation method for forest canopy high-resolution remote sensing images, and steps of the method comprise the following steps: 1) determining remote sensing images and corresponding canopy height model information, observation geometric parameters and sky scattering light proportion; collecting forest canopy structure parameters and forest spectrum parameters; 2) calculating a single scattering bidirectional reflectance factor BRF self_sgl only considering that a single scattering occurs inside a target pixel and a single scattering bidirectional reflectance factor BRF adj_sgl affected by adjacent pixels; 3) calculating a multiple scattering bidirectional reflectance factor BRF self_mul only considering that a multiple scattering occurs inside a target pixel and a multiple scattering bidirectional reflectance factor BRF adj_mul affected by adjacent pixels; 4) calculating a total bidirectional reflectance factor BRF self = BRF self_sgl + BRF self_mul , wherein the total bidirectional reflectance factor BRF adj = BRF adj_sgl + BRF adj_mul .

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of remote sensing, and relates to a bidirectional reflectance factor analytical calculation method for high spatial resolution forest canopy. BACKGROUND

[0002] With the increasing richness of remote sensing observation means such as high-resolution satellites, unmanned aerial vehicles, and laser radars, it is possible to obtain high-resolution remote sensing images, which provides an effective approach to accurately describing the spatial variation of forest canopy vegetation parameters. Compared with medium and low spatial resolution remote sensing images, high-resolution remote sensing images provide more abundant texture and structure information, and are widely used in qualitative researches such as tree classification and tree crown extraction. In particular, with the development of unmanned aerial vehicles and laser radars in recent years, some studies are committed to successfully extracting structural parameters such as tree height, crown width, and coverage through stereo image pairs of unmanned aerial vehicle images and laser radar point cloud data. However, while these high spatial resolution forest canopy qualitative researches are rapidly developing, the related quantitative researches still use the vegetation index commonly used in medium and low spatial resolution vegetation remote sensing, or directly migrate the traditional mechanism model. With the increase of the spatial resolution of remote sensing images, the forest canopy presents stronger spatial heterogeneity, and the individual geometric structure characteristics are more obvious. Compared with medium and low spatial resolution remote sensing images, the interaction mechanism with light is quite different. The light shielding and cross-radiation difference between adjacent pixels on high spatial resolution remote sensing images cannot be ignored. The reflected radiation of a remote sensing pixel is not only affected by the incident sunlight, but also affected by the radiation from adjacent pixels. This means that on high-resolution remote sensing images, the radiation brightness information of a pixel not only reflects the radiation information of the pixel itself, but also contains the radiation contribution from adjacent pixels. This adjacent pixel radiation influence leads to the complexity of the radiation brightness information of the remote sensing image, affecting the inversion accuracy of the vegetation parameters.

[0003] Traditional vegetation remote sensing mechanism models, such as the 4-Scale model and the unified model, are mainly used for vegetation parameter inversion of medium and low spatial resolution remote sensing images, mainly focusing on the internal radiation transmission process of the pixel, without considering the radiation influence caused by the mutual shielding and cross-radiation difference of adjacent pixels. With the increase of the spatial resolution of remote sensing images, the forest canopy presents stronger spatial heterogeneity, and the adjacent pixel radiation influence cannot be ignored, which leads to the difficulty of traditional models in accurately describing the radiation transmission process of the forest canopy, reduces the applicability, and further limits the inversion accuracy. In addition, although computer simulation such as the LESS model (LargE-Scale remote sensing data and image Simulation framework) can provide fine radiation transmission analysis, it has high complexity and no analytical solution expression, which is difficult to directly apply to actual vegetation parameter inversion.

[0004] The technical defects of the traditional scheme can be summarized as follows:

[0005] 1) The influence of adjacent pixel radiation is not considered: The traditional radiation transfer model mainly focuses on the inside of the pixel, ignores the radiation influence of adjacent pixels of the forest canopy in high spatial resolution remote sensing images, and leads to the decrease of the vegetation parameter inversion accuracy at high spatial resolution.

[0006] 2) Lack of consideration of the complexity of the forest canopy geometry: The traditional radiation transfer model uses parameters in the average statistical sense, which cannot accurately describe and depict the three-dimensional geometry and spatial heterogeneity of the forest canopy at high spatial resolution, so that the model cannot meet the needs of fine monitoring.

[0007] 3) Limited practicability: Although the computer simulation model can simulate the radiation transfer process of the forest canopy at high spatial resolution, the model has no analytical expression, and the input and simulation process of the model are relatively complex, which leads to the limited practicability of the model and makes it difficult to directly guide the extraction of vegetation parameters. SUMMARY

[0008] In view of the problems in the prior art, the purpose of the present application is to provide a bidirectional reflectance factor calculation method for high spatial resolution forest canopy. The present application aims to fully consider the influence of adjacent pixel radiation, use the canopy height model (CHM) which can accurately describe the geometric structure information of the forest canopy, introduce the path length information, and construct a forest canopy remote sensing mechanism model suitable for high spatial resolution, so as to fully exert the advantages of high-resolution images in depicting the geometric structure and spatial heterogeneity of trees, and improve the vegetation parameter inversion accuracy of high spatial resolution forest canopy. The technical scheme of the present application can effectively improve the vegetation parameter inversion accuracy of high spatial resolution forest canopy by considering the influence of adjacent pixel radiation, combining the geometric optical model and the spectral invariance theory, and giving the analytical expression of the bidirectional reflectance factor of high-resolution remote sensing image of forest canopy. The present application can meet the needs of ecology for fine structure of small-scale stand and community forest canopy, and has wide application prospect.

[0009] The present application establishes a high spatial resolution forest canopy remote sensing model considering the influence of adjacent pixel radiation, and improves the vegetation parameter inversion accuracy at high spatial resolution. The present application introduces the geometric structure information of the forest canopy, more accurately describes the forest canopy structure at high spatial resolution, and meets the needs of ecology. The present application provides an analytical expression for high spatial resolution forest canopy radiation transfer, which is beneficial to the development of corresponding high-precision inversion methods based on the present application.

[0010] The technical scheme of the present application is as follows:

[0011] A bidirectional reflectance factor calculation method for high-resolution remote sensing image of forest canopy, comprising the following steps:

[0012] 1) Determine the basic information of remote sensing image and corresponding canopy height model CHM, including image range, image spatial resolution and canopy height model CHM spatial resolution; collect forest canopy structure parameters and forest spectral parameters, the forest canopy structure parameters include leaf area volume density FAVD, G function, canopy height model CHM, the spectral parameters include leaf reflectivity r L ; Determine the observation geometry parameters and the proportion of sky scattered light, the observation geometry parameters include solar zenith angle θ s , solar azimuth angle observation zenith angle θ v and observation azimuth angle

[0013] 2) Calculate the single scattering bidirectional reflectance factor BRF self_sgl only considering the internal scattering process of the target pixel by using the forest canopy structure parameters, the spectral parameters and the observation geometry parameters;

[0014] 3) Calculate the single scattering bidirectional reflectance factor BRF adj_sgl of each target pixel affected by adjacent pixels by using the forest canopy structure parameters, the spectral parameters and the observation geometry parameters;

[0015] 4) Calculate the multiple scattering bidirectional reflectance factor BRF self_mul only considering the internal scattering process of the target pixel by using the forest canopy structure parameters, the spectral parameters, the observation geometry parameters and the intermediate variable leaf area index LAI obtained from the single scattering process;

[0016] 5) Calculate the multiple scattering bidirectional reflectance factor BRF adj_mul of each target pixel affected by adjacent pixels by using the forest canopy structure parameters, the spectral parameters, the observation geometry parameters and the intermediate variable leaf area index LAI obtained from the single scattering process;

[0017] 6) Add the canopy single scattering bidirectional reflectance factor BRF self_sgl only considering the internal scattering process of the target pixel and the canopy multiple scattering bidirectional reflectance factor BRF self_mul only considering the internal scattering process of the target pixel to obtain the total bidirectional reflectance factor BRF self only considering the internal scattering process of the target pixel;

[0018] 7) Add the canopy single scattering bidirectional reflectance factor BRF adj_sgl considering the influence of adjacent pixel radiation and the canopy multiple scattering bidirectional reflectance factor BRF adj_mul considering the influence of adjacent pixel radiation to obtain the total bidirectional reflectance factor BRF of each target pixel affected by adjacent pixelsadj .

[0019] Furthermore, in step 2), according to ERF self_sgl =K g ·r g +K l ·r l The bidirectional reflectivity factor (BRF) for single-scattering was calculated. self_sgl Among them, K g K represents the proportion of soil area exposed to sunlight. l r represents the proportion of leaf area exposed to light. g r represents the reflectivity of soil under sunlight. l Represents the reflectivity of the leaf when exposed to light;

[0020] 21) Calculate the proportion of sunlit soil area K using a path length distribution model. g The ratio of the area of ​​soil exposed to sunlight, K g K from the large pores between the tree canopy g1 Small pores K in the shaded ground g2 Soil light spots K observed through pores in the tree canopy g3 Composition; where K g1 =P(l(Ω) s =0),l(Ω v =0)) represents the path length l (Ω) in the direction of solar incidence. v The path length l (Ω) in the sensor observation direction is 0. s The joint probability when K is 0; g2 =P(l(Ω) s >0), l(Ω v =0))·P(Ω s |l(Ω s >0), l(Ω v =0), P(l(Ω) s >0), l(Ω v =0)) represents the path length l (Ω) of the sensor observation direction. v When ) is 0 and the path length l (Ω) of the solar incident direction is 0 s The joint probability that Ω is greater than 0, P(Ω) s |l(Ω s >0), l(Ω v =0)) refers to calculating porosity based on the area ratio of shaded soil using path length distribution. Where G(Ω) s () represents the projection ratio of the leaf towards the direction of solar incidence, FAVD is the leaf area volume density, l g2 (Ω s p represents the path length of the shaded soil portion along the direction of solar incidence.l_g2 (Ω s ) is the probability density distribution of path length l g2 (Ω s ), l max , l min represent the maximum and minimum of path length l g2 (Ω s ), respectively; K g3 = P(l(Ω v > 0)) · P(Ω v | l(Ω v > 0)), P(l(Ω v > 0)) · P(Ω v | l(Ω v > 0)) is the proportion of partial porosity of path length l(Ω v ) greater than 0 in the observation direction, l g3 (Ω v ) represents the path length of the tree crown in the vertical observation direction, p l_g3 (Ω v ) is the probability density distribution of path length l g3 (Ω v ), G(Ω v ) represents the projection proportion of leaves in the vertical observation direction;

[0021] 22) The bidirectional reflectance factor BRF l of the illuminated leaf area proportion K l = BRF lc + BRF lt ; by calculating the illuminated tree crown area proportion K c and the shaded tree crown area proportion K t , the illuminated tree crown area proportion K c is converted into the bidirectional reflectance factor BRF lc of the illuminated leaf, and the shaded tree crown area proportion K t is converted into the bidirectional reflectance factor BRF lt of the illuminated leaf; wherein the illuminated tree crown area proportion K c = P(l(Ω s = 0), l(Ω v > 0)) · [1 - P(Ω v | l(Ω s = 0), l(Ω v > 0)], and the shaded tree crown area proportion K t = P(l(Ω s > 0), l(Ω v > 0)) · [1 - P(Ω v | l(Ω s(>0), l(Ω v (>0))]; P(l(Ω s =0), l(Ω v (>0)) refers to the joint probability that the path length in the observation direction is greater than 0 while the path length in the solar incident direction is 0.

[0022] Furthermore, in step 3), first calculate the shielding factor S of the illuminated canopy component c and the shielding factor S of the illuminated soil component g , then derive the new illuminated canopy component K c ′ and the new illuminated soil component K′ g considering the influence of adjacent pixel radiation, and then calculate the single-scattering bidirectional reflectance factor BRF adj_sgl of each target pixel affected by adjacent pixels.

[0023] Furthermore, the specific method for calculating the single-scattering bidirectional reflectance factor BRF adj_sgl of each target pixel affected by adjacent pixels is as follows:

[0024] 31) Calculate the large pore proportion P c (l = 0, Ω s ) = P(F c (l = 0, Ω s )|K c ), the small pore proportion formed after being blocked by adjacent pixels Then calculate the shielding factor S of the illuminated canopy c = [1 - P c (l = 0, Ω s )]·[1 - P c (l > 0, Ω s )]; then calculate the area proportion K c ′ of the illuminated canopy component of the target pixel after being affected by adjacent pixel radiation = K c ·(1 - S c ); l(Ω s |K c ) represents the path length in the solar incident direction of the original illuminated canopy part of the target pixel after being affected by adjacent pixels, p l (Ω s |K c ) represents the probability density distribution of the path length l(Ω s |K c );

[0025] 32) Calculate the large pore proportion P g (l = 0, Ω s)=P(F g (l=0,Ω s )|K g ), the proportion of small apertures formed after being blocked by neighboring pixels and the shading factor S on the soil g =[1-P g (l=0,Ω s )]·[1-P g (l>0,Ω s Then calculate the area ratio K′ of the target pixel's irradiated soil component after being affected by the radiation of neighboring pixels. g =K g ·(1-S g );l(Ω s |K g p represents the path length of the sun's incident direction in the original sunlight-illuminating part of the target pixel after being affected by neighboring pixels. l (Ω s |K g ) represents the path length l (Ω) s |K g The probability density distribution of ).

[0026] 33) The first-scattering bidirectional reflectivity factor BRF is calculated. adj_sgl = (1-β)·(K′) g ·r s +BRF l ′); where β is the proportion of sky-scattered light.

[0027] Furthermore, in step 4), firstly, based on the ratio of the illuminated canopy area K... c , the proportion of shaded canopy area K t and the vertical observation path length distribution F(l>0,Ω) v ), calculated Then, based on the probability p of a photon scattering from a leaf and then colliding again with a leaf inside the canopy, 1 The probability of multiple collisions p m The calculation ignores the soil and only considers the bidirectional reflectivity factor (BRF) of multiple scattering between the canopy layers. self_vmul And considering the bidirectional reflectance factor (BRF) between the soil and the interior of the canopy. self_smul Then, the canopy multiple scattering bidirectional reflectivity factor (BRF) was calculated when only the target pixel was considered. self_mul =BRF self_vmul +BRF self_smul .

[0028] Furthermore, in step 5), the bidirectional reflectance factor (BRF) of multiple scattering caused by neighboring pixels is calculated for each target pixel.adj_mul The method is as follows: First, based on the ratio of the canopy area illuminated by sunlight, K... c , the proportion of shaded canopy area K t and the vertical observation path length distribution F(l>0,Ω) v ), calculated Then, the bidirectional reflectivity factor of multiple scattering from neighboring pixels to the target pixel is calculated. Among them, LAI i LAI represents the leaf area index of the i-th surrounding pixel. self Leaf area index (BRF) representing the current target cell mul_i W represents the multiple scattering BRF calculated for the i-th surrounding pixel when only internal multiple scattering is considered. i The distance-dependent normalized weight matrix is ​​then used; the multiple scattering bidirectional reflectivity factor (BRF) is then calculated. adj_mul =BRF self_mul +BRF others_mul .

[0029] Furthermore, p 1 =0.7exp(k1·LAI)-0.66exp(k2·LAI), Where, k1 = 0.0045exp(1.2555cosθ), k2 = 0.1982lncosθ - 0.7146, i(θ) is the interception probability of the canopy in the direction of θ, where θ is the solar zenith angle.

[0030] Furthermore, Where i0=β·i D +(1-β)·i S ω represents the average probability of the canopy intercepting direct sunlight and sky-scattered light, and r represents the leaf primary scattering albedo. s R represents soil reflectance. dn T represents the albedo at the base of the canopy. dn and T up p(Ω) represents the descending and ascending transmittance of the canopy, respectively. s ) and p(Ω v ) represent the porosity in the direction of solar incidence and the direction of sensor observation, respectively.

[0031] A server is characterized by comprising a memory and a processor, the memory storing a computer program configured to be executed by the processor, the computer program including instructions for performing the steps of the methods described above.

[0032] A computer readable storage medium, having stored thereon a computer program, wherein the computer program is executed by a processor to implement the steps of the method.

[0033] The key points of the present application include:

[0034] ①Combining the geometric optics model and the spectral invariance theory, the forest canopy high-resolution remote sensing pixel BRF is divided into single scattering contribution and multiple scattering contribution, and then the geometric optics model and the spectral invariance theory are used for description respectively, so that the forest canopy radiation transmission process is more comprehensively described.

[0035] ②Considering the radiation influence of adjacent pixels, the shielding factor and the neighborhood convolution algorithm are introduced, the contribution of adjacent pixels to the target pixel BRF is quantitatively expressed, and the BRF of the pixel itself and adjacent pixels is distinguished.

[0036] ③The concept of shielding factor is adopted to quantitatively express the influence of adjacent pixels on the single scattering of the target pixel.

[0037] ④The neighborhood convolution algorithm is used to quantitatively express the influence of adjacent pixels on the multiple scattering of the target pixel.

[0038] ⑤The path length distribution model is introduced, and the high-resolution CHM data is combined, so that the forest canopy spatial distribution and three-dimensional structure under high spatial resolution are more accurately described, and the model precision is improved.

[0039] The advantages of the present application are as follows:

[0040] ①The high spatial resolution vegetation parameter inversion precision is improved: considering the radiation influence of adjacent pixels, the forest canopy radiation transmission process under high spatial resolution conditions can be more accurately simulated, and the forest canopy vegetation parameter inversion precision based on the physical model is improved.

[0041] ②The model effectively distinguishes the reflection radiation contribution of the pixel itself and the radiation influence of adjacent pixels, so that the forest canopy vegetation parameters inverted based on the model can truly represent the vegetation conditions on the high-resolution remote sensing pixel.

[0042] ③The ecological needs are met: the high spatial resolution CHM data is used, the forest canopy geometric structure information is introduced, the model precision for describing the three-dimensional geometric structure of the forest canopy is improved, and the needs of ecology for small-scale clumps and community forest canopy fine structure are met.

[0043] ④The model has strong applicability: the model structure is simple and easy to implement, and compared with the complex model such as LESS, the model is more suitable for practical application. BRIEF DESCRIPTION OF DRAWINGS

[0044] Figure 1 The method flowchart of the present application.

[0045] Figure 2 BRF values calculated by the model versus simulated BRF values by the LESS model;

[0046] Figure 2 (a) Comparison of the calculated values by the APPLE-GO model in the red band versus the simulated values by the LESS model in the red band, considering only the internal scattering of the target pixel;

[0047] Figure 2 (b) Comparison of the calculated values by the APPLE-GO model in the near infrared band versus the simulated values by the LESS model in the near infrared band, considering only the internal scattering of the target pixel;

[0048] Figure 2 (c) Comparison of the calculated values by the APPLE-GO model in the red band versus the simulated values by the LESS model in the red band, considering the influence of the radiation of the neighboring pixels;

[0049] Figure 2 (d) Comparison of the calculated values by the APPLE-GO model in the near infrared band versus the simulated values by the LESS model in the near infrared band, considering the influence of the radiation of the neighboring pixels.

[0050] Figure 3 Comparison of the BRF values calculated by the APPLE-GO model versus the simulated values by the LESS model for the pixel near the edge in Scene 1;

[0051] Figure 3 (a) Schematic diagram of the position of the pixel near the edge in Scene 1,

[0052] Figure 3 (b) Polar coordinate comparison of the BRF values calculated by the APPLE-GO model in the red band versus the simulated values by the LESS model in the red band, considering only the internal scattering of the target pixel (the left side is the calculation result of the APPLE-GO model, and the right side is the simulation result of the LESS model);

[0053] Figure 3 (c) Polar coordinate comparison of the BRF values calculated by the APPLE-GO model in the near infrared band versus the simulated values by the LESS model in the near infrared band, considering only the internal scattering of the target pixel (the left side is the calculation result of the APPLE-GO model, and the right side is the simulation result of the LESS model);

[0054] Figure 3 (d) Polar coordinate comparison of the BRF values calculated by the APPLE-GO model in the red band versus the simulated values by the LESS model in the red band, considering the influence of the radiation of the neighboring pixels (the left side is the calculation result of the APPLE-GO model, and the right side is the simulation result of the LESS model);

[0055] Figure 3(e) The polar coordinate plot of the BRF calculated by the APPLE-GO model considering the influence of the adjacent pixel radiation in the near-infrared band and the simulated value of the LESS model in the near-infrared band (the left side is the calculation result of the APPLE-GO model, and the right side is the simulation result of the LESS model);

[0056] Figure 4 for the reflectance image contrast chart;

[0057] Figure 4 (a) Beijing No. 3 satellite RGB reflectance composite image;

[0058] Figure 4 (b) The reflectance RGB true color composite image calculated by the APPLE-GO model. DETAILED DESCRIPTION

[0059] The application will be further described in detail below with reference to the accompanying drawings, and the examples are only used to explain the application and not to limit the scope of the application.

[0060] The high spatial resolution forest canopy radiation transfer model APPLE-GO proposed in the application is established on the basis of in-depth analysis of the radiation influence of the adjacent pixels of the forest canopy. The application can accurately simulate the radiation influence of the adjacent pixels, effectively distinguish the BRF (bidirectional reflectance factor) of the target pixel itself and the radiation influence of the adjacent pixels, and thus improve the inversion accuracy of vegetation parameters such as leaf area index (LAI) and leaf chlorophyll content (LCC).

[0061] The premise assumptions of the APPLE-GO model are as follows: ① Sensor vertical observation: the model assumes that the sensor vertically observes the forest canopy downward, which is consistent with the observation mode of the unmanned aerial vehicle and most high-resolution satellites; ② Forest canopy CHM is known and has high spatial resolution: the model needs high-precision CHM data to describe the geometric structure of the forest canopy, so as to ensure that the model can accurately calculate the path length distribution and the shielding factor; ③ Leaf area density is known at the modeling scale: the model assumes that the leaf area density is known at the modeling scale, which can be obtained by ground measurement extraction; ④ The leaf at the modeling scale obeys the Poisson distribution; ⑤ The trunk height of the trees in the forest canopy is approximately uniform: the model assumes that the trunk height of the trees in the forest canopy is approximately uniform, which simplifies the model calculation, but can be extended by considering the difference in trunk height; ⑥ The leaf and soil background are regarded as Lambertian bodies; ⑦ The ground is regarded as a flat ground: the model assumes that the leaf and soil background are Lambertian bodies, that is, the reflectivity does not change with the observation angle, which simplifies the model calculation; ⑧ The wood components in the forest canopy scene are not considered.

[0062] The specific construction process and detailed formula of the APPLE-GO model are as follows:

[0063] The APPLE-GO model separates the forest canopy BRF into single scattering and multiple scattering contributions, and models them separately.

[0064] For single scattering, based on the geometric optics model, the path length distribution is introduced to calculate the shadowing factor of the two components of the target pixel, i.e. the illuminated soil and the illuminated canopy, to express the influence of the neighboring pixels on the single scattering of the target pixel.

[0065] For multiple scattering, based on the spectral invariance theory, the neighborhood convolution algorithm is used to consider the influence of the neighboring pixels on the multiple scattering of the target pixel.

[0066] ① Single scattering modeling

[0067] Single scattering refers to the BRF of the target pixel, which is the first collision of the sun photon with the canopy element and is reflected into the sensor. Single scattering is most affected by the canopy structure, which can be reflected in the hotspot effect. Single scattering is the main reason for the directionality of the BRF observed by the sensor. According to the geometric optics theory, single scattering is mainly determined by the BRF contributions of the two parts of the illuminated soil and the illuminated leaf, as shown in equation (1).

[0068] BRF self_sgl =∑K(Ω s )·r=K g ·r g +K l ·r l (1)

[0069] where K(Ω s ) represents the area proportion of the illuminated component in the scene, represents the solar incidence direction (θ s is the solar zenith angle, is the solar azimuth angle), r is the reflectivity of the corresponding component, K g represents the area proportion of the illuminated soil, K l represents the area proportion of the illuminated leaf, r g and r l represent the reflectivity of the two illuminated components, respectively. The illuminated components include the illuminated soil and the illuminated leaf. This study assumes that both parts are Lambertian, so only the area proportion changes with the solar incidence direction, and the reflectivity does not change with the angle.

[0070] When only considering single scattering within the pixel, the illuminated soil area proportion K g is determined by the large pores K g1 between the canopy, the small pores K g2 in the shadow ground, and the soil patches K g3Composition. The APPLE-GO model calculates the area proportion of these three parts by the path length distribution model, which is composed of the following equations (2)-(6).

[0071] K g1 = P(l(Ω s = 0), l(Ω v = 0)) (2)

[0072] In the equation, l represents the path length, l(Ω s = 0) represents the path length of 0 in the direction of the sun incidence, and l(Ω v = 0) represents the path length of 0 in the direction of the sensor observation. The meaning expressed by the equation is the joint probability when the path lengths in both directions are 0. For K g2 , it can be expressed based on the path length distribution as follows:

[0073] K g2 = P(l(Ω s > 0), l(Ω v = 0)) · P(Ω s | l(Ω s > 0), l(Ω v = 0)) (3)

[0074] In the equation, P(l(Ω s > 0), l(Ω v = 0)) represents the joint probability of the path length of 0 in the direction of the sensor observation and the path length greater than 0 in the direction of the sun incidence, which represents the area proportion of the shadow soil in the traditional sense; and P(Ω s | l(Ω s > 0), l(Ω v = 0)) is a conditional probability, which refers to the expression of the porosity obtained by using the path length distribution on the basis of the area proportion of the shadow soil, and can be expressed as:

[0075]

[0076] In the equation, G(Ω s ) represents the projection proportion of the leaf in the direction of the sun incidence, FAVD is the leaf area volume density, l g2 (Ω s ) = l(Ω s | l(Ω s > 0), l(Ω v = 0)) represents the path length of the shadow soil part in the direction of the sun incidence, p l_g2 (Ω s ) = p l (Ω s | l(Ω s > 0), l(Ωv =0)) represents the probability density distribution of this part of the path length, and the upper and lower limits of integration represent the maximum and minimum values ​​of the path length, respectively.

[0077] K g3 =P(l(Ω) v >0))·P(Ω v |l(Ω v >0)) (5)

[0078] In the formula, P(l(Ω) v >0))·P(Ω v |l(Ω v >0)) represents the proportion of porosity in the portion of the tree canopy where the path length is greater than 0 in the observation direction (i.e., the vertical direction). P(Ω) v |l(Ω v >0)) is the conditional probability, which can be expressed as:

[0079]

[0080] Where G(Ω) v () represents the projection ratio of the leaf in the vertical observation direction, FAVD is the leaf area volume density, l g3 (Ω v )=l(Ω v |l(Ω v >0)) represents the path length of the canopy portion in the vertical observation direction, p l_g3 (Ω v ) = p l (Ω v |l(Ω v >0)) represents the probability density distribution of this part of the path length, and the upper and lower limits of integration represent the maximum and minimum values ​​of the path length, respectively.

[0081] The area of ​​illuminated leaves originates from both the illuminated and shaded canopies. Therefore, we first calculate the area ratio of the illuminated and shaded canopies, and then convert it into the BRF contribution of illuminated leaves. The area ratio of the illuminated canopy is K. c The calculation formula is shown in equation (7), where K is the area ratio of the shaded canopy. t As shown in Equation (8), the BRF contribution of the illuminated blades is obtained after conversion, as shown in Equations (9) to (16).

[0082] K c =P(l(Ω) s =0), l(Ω) v >0))·[1-P(Ω v |l(Ω s =0), l(Ω) v >0))] (7)

[0083] K t = P(l(Ω s > 0), l(Ω v > 0)) · [1 - P(l(Ω v | l(Ω s > 0), l(Ω v > 0))] (8)

[0084] where P(l(Ω s = 0), l(Ω v > 0)) means the joint probability of the part with path length greater than 0 (i.e. the canopy part) in the observation direction (i.e. vertical direction) and the part with path length equal to 0 in the sun direction. 1 - P(l(Ω v ) represents the part that is really illuminated after removing the internal porosity.

[0085] The formula of converting the illuminated canopy to the illuminated leaf (refer to Mu X, Hu R, Zeng Y, et al. 2017. Estimating structural parameters of agricultural crops from ground-based multi-angular digital images with a fractional model of sun and shade components. Agricultural and Forest Meteorology, 246:162-177.):

[0086]

[0087] where Γ(Ω s , Ω v ) / π is the scattering phase function, representing the proportion of the light that is scattered to the observation direction after the collision between the light and the leaf, g(θ L ) / 2π is the leaf inclination distribution function (based on the upper surface of the leaf), γ(Ω s → Ω v ) / π is the leaf-scale scattering phase function, multiplied by the leaf albedo r L when the exit and entrance occur in the same side, or multiplied by the leaf transmittance t L when the exit and entrance occur in different sides. P(z, Ω s , Ω v ) is the bidirectional porosity at height z, which is equal to the product of the porosities in the observation and sun directions, multiplied by the hot spot factor, expressed as:

[0088] P c(z, Ω s , Ω v ) = exp(-G(Ω v ) · FAVD · z) · exp(-G(Ω s ) · FAVD · z / cosθ s ) · C HS (12)

[0089]

[0090] where C HS is the hotspot factor. In equation (13), tanθ s is a special form for the normal view, and s L = π 2 · d L / 16 is a form for the spherical distribution of leaves, where d L is the leaf size.

[0091] The formula for converting the shadowed crown to the illuminated leaves:

[0092]

[0093] In equation (15), l represents the average path length in the direction of the sun's incidence. Therefore, the BRF contribution BRF l of the illuminated leaf area fraction K l can be written as the sum of contributions from both the illuminated crown and the shadowed crown:

[0094] BRF l = BRF lc + BRF lt (16)

[0095] Considering the effect of the adjacent pixel's shading on the radiation, the concept of the shading factor is introduced. In the APPLE-GO model, the shading factors of the illuminated crown and the illuminated soil component are defined respectively, and the proportions of the area reduction of these two components caused by the shading effect from the adjacent pixel are calculated. This proportion is called the shading factor. The shading factor is calculated by the path length distribution model, which considers the probabilities of removing the large pores and the small pores in the direction of the sun.

[0096] The proportion of the large pores formed on the illuminated crown after being shaded by the adjacent pixel is calculated as shown in equation (17), and the proportion of the small pores formed on the illuminated crown after being shaded by the adjacent pixel is calculated as shown in equation (18). The shading factor S c of the illuminated crown is expressed in the form of the probabilities of the large pores and the small pores, as shown in equation (19).

[0097] P c (l = 0, Ω s ) = P(Fc (l = 0, Ω s ) | K c ) (17)

[0098]

[0099] S c = [1 - P c (l = 0, Ω s )] · [1 - P c (l > 0, Ω s )] (19)

[0100] In equation (18), l(Ω s | K c ) represents the path length in the direction of the sun's incidence of the portion of the target pixel's illuminated crown that is affected by the neighboring pixels, and p l (Ω s | K c ) represents the probability density distribution of this portion of the path length. The proportion of the target pixel's illuminated crown component that is affected by the radiation of the neighboring pixels, K c ', is:

[0101] K c ' = K c · (1 - S c ) (20)

[0102] Similarly, the proportion of the illuminated soil that is affected by the neighboring pixels is calculated as shown in equation (21), and the proportion of the illuminated soil that is affected by the neighboring pixels is calculated as shown in equation (22). The shadowing factor S g is expressed in the form of the probabilities of large and small pores, as shown in equation (23).

[0103] P g (l = 0, Ω s ) = P(F g (l = 0, Ω s ) | K g ) (21)

[0104]

[0105] S g = [1 - P g (l = 0, Ω s )] · [1 - P g (l > 0, Ω s )] (23)

[0106] In equation (22), l(Ω s | K grepresents the path length in the direction of the incident solar radiation of the target pixel after the illuminated soil fraction is affected by the neighboring pixels, p l (Ω s |K g represents the probability density distribution of this part of the path length. The illuminated soil fraction of the target pixel after being affected by the radiation of the neighboring pixels is K' g That is:

[0107] K′ g = K g ·(1-S g ) (24)

[0108] The new illuminated canopy fraction and the new illuminated soil fraction calculated by formula (20) and formula (24) are respectively replaced by the illuminated canopy and the illuminated soil area fraction originally considered only within the pixel, and the once scattering expression of the target pixel after being affected by the neighboring pixels can be obtained:

[0109] BRF adj_sgl = (1-β) · (K' g · r s + BRF l ') (25)

[0110] ②Multiple scattering modeling

[0111] Multiple scattering refers to the BRF presented by the photons entering the vegetation canopy and finally entering the sensor after multiple collisions between components. With the increase of collision times, multiple scattering tends to be approximately isotropic. Therefore, the heterogeneity of the canopy has less effect on multiple scattering than on once scattering. If only once scattering is considered and multiple scattering is ignored, the numerical value of the model simulated BRF will deviate. In the visible light band where the once scattering albedo of the leaf is low, the influence of multiple scattering is small, and once scattering dominates the reflection radiation of the canopy. However, in the near-infrared band where the once scattering albedo of the leaf is high, the influence of multiple scattering cannot be ignored.

[0112] When only multiple scattering occurring within the pixel is considered, it can be divided into ignoring soil, only considering the contribution of multiple scattering between the canopy (black soil process), and considering the contribution of multiple scattering between the soil and the canopy (white soil process). The model defines a once collision probability p 1 , that is, the probability of the photon being scattered by the leaf and colliding with the leaf in the canopy again, which is related to the incident direction angle, and the calculation formula is shown in formula (26)-(28). The model also defines a multiple collision probability p m , m is the number of collisions, m>1; that is, the probability of the photon being scattered by the leaf and colliding with the leaf in the canopy again, which is independent of the angle, as shown in formula (29)-(30).

[0113] p1 = 0.7exp(k1-LAI) - 0.66exp(k2-LAI) (26)

[0114] where k1and k2are the fitted expressions of the solar zenith angle, which are shown as follows:

[0115] k1= 0.0045exp(l.2555cos0) (27)

[0116] k2= 0.1982ln cos0- 0.7146 (28)

[0117] where 0 is the solar zenith angle.

[0118]

[0119] where i D is the canopy hemispheric interception probability, i.e., the integral of the interception probability of the canopy in all directions in space:

[0120]

[0121] where i(0) is the interception probability of the canopy in the direction of 0.

[0122] In equations (26)-(30), the leaf area index LAI is calculated from the parameters obtained from the first scattering process: the illuminated canopy area fraction K c , the shadowed canopy area fraction K t , and the vertical observation path length distribution (excluding the part with a path length of 0, F(l>0,0 v ).

[0123]

[0124] Based on the single-reverberation probability p 1 and the multiple-reverberation probability p m , the expression for the contribution of multiple scattering between the canopy only (black soil process) is shown in equation (32), and the expression for the contribution of multiple scattering between the soil and the interior of the canopy (white soil process) is shown in equation (33).

[0125]

[0126] In equation (32), i0= β·i D +(1-β)·i S , which represents the average interception probability of the canopy for the direct sunlight and the sky scattered light, and ω is the first scattering albedo of the leaf, which can be approximately equal to the sum of the reflectance and the transmittance of the leaf. In equation (33), r s represents the reflectance of the soil, and R dn represents the bottom albedo of the canopy, and Tdn and T up represent the canopy downward and upward transmittance, p(Ω s ) and p(Ω v ) represent the porosity in the direction of the sun and the sensor, respectively. R dn , T dn and T up are expressed as:

[0127]

[0128]

[0129] BRF of the target pixel only considering the multiple scattering within the canopy self_mul is expressed as the sum of the contribution of the multiple scattering between the canopy and the contribution of the multiple scattering between the soil and the interior of the canopy:

[0130] BRF self_mul = BRF self_vmul + BRF self_smul (37)

[0131] When the influence of the adjacent pixels is considered, the model uses the neighborhood convolution algorithm to obtain the contribution of the multiple scattering of the adjacent pixels to the target pixel by convolving the relative size of the LAI of the neighborhood pixels of the target pixel and the planar Euclidean distance as the weight:

[0132]

[0133] where LAI i represents the leaf area index of the i-th surrounding pixel, LAI self represents the leaf area index of the current target pixel, BRF mul_i represents the multiple scattering BRF calculated by only considering the internal multiple scattering of the i-th surrounding pixel, and W i is a distance-dependent normalized weight matrix.

[0134] The expression of the multiple scattering of the target pixel under the influence of the adjacent pixels is obtained by adding the contribution of the multiple scattering between the pixels and the increment of the multiple scattering under the influence of the adjacent pixels:

[0135] BRF adj_mul = BRF self_mul + BRF others_mul (39)

[0136] ③ Expression of the total BRF of the target pixel:

[0137] The total BRF of the target pixel includes the contribution of the total BRF of the target pixel itself BRF selfand the total BRF contribution of the target pixel considering the influence of adjacent pixels BRF adj The BRF contribution of the target pixel itself BRF self is expressed as the single scattering BRF of the canopy only considering the target pixel BRF self_sgl and the multiple scattering BRF of the canopy only considering the target pixel BRF self_mul The sum of the two parts:

[0138] BRF self = BRF self_sgl + BRF self_mul (40)

[0139] The BRF contribution of the target pixel considering the influence of adjacent pixels BRF adj is expressed as the single scattering BRF of the canopy considering the influence of adjacent pixels BRF adj_sgl and the multiple scattering BRF of the canopy considering the influence of adjacent pixels BRF adj_mul The sum of the two parts:

[0140] BRF adj = BRF adj_sgl + BRF adj_mul (41)

[0141] Summary of input parameters of the APPLE-GO model:

[0142] Forest canopy structure parameters: leaf area volume density (FAVD), G function (G), CHM, etc.

[0143] Spectral parameters: leaf reflectance (r L ), leaf transmittance (t L ), soil reflectance (r s ), etc.

[0144] Observation geometry parameters: solar zenith angle (θ s ), solar azimuth angle observation zenith angle (θ v ), observation azimuth angle , etc.

[0145] Other parameters: proportion of sky scattered light (β), etc.

[0146] Summary of model output:

[0147] The BRF contribution of the target pixel itself: the BRF only considering the scattering process within the pixel.

[0148] The BRF contribution of the target pixel considering the influence of adjacent pixels: the BRF considering the influence of adjacent pixels.

[0149] The flow of the model is shown in Figure 1 , which includes:

[0150] 1. Data preparation:

[0151] Determine the basic information of remote sensing image and corresponding canopy height model (CHM), including image range, image spatial resolution and CHM spatial resolution, etc.

[0152] Collect necessary input data, including forest canopy structure parameters such as leaf area volume density FAVD, leaf reflectance r L and other spectral characteristics data.

[0153] Determine observation geometry parameters such as solar zenith angle θ s , solar azimuth angle , observation zenith angle θ v and observation azimuth angle , and other parameters such as sky diffuse light proportion, etc.

[0154] 2. Calculate the first scattering BRF considering only the target pixel itself: self_sgl Using leaf area volume density FAVD, G function, canopy height model CHM and other structure parameters, leaf reflectance r L and other spectral parameters, and solar zenith angle θ s and other observation geometry parameters, based on equations (1)-(16), the first scattering BRF of the target pixel considering only the scattering within the target pixel itself and without considering the influence of adjacent pixels can be calculated. G function is a well-known parameter in the field of quantitative remote sensing of vegetation, which means the projection coefficient of unit leaf area on the plane perpendicular to the observation direction. Specifically, based on the geometric optical model (as shown in equation 1), the illuminated leaf component (as shown in equations 2-6) and the illuminated soil component (as shown in equations 7-16) are calculated respectively, and the sum of the two is the BRF self_sgl (as shown in equation 1).

[0155] 3. Calculate the first scattering BRF of each target pixel affected by adjacent pixels: self_sgl Using leaf area volume density FAVD, G function, canopy height model CHM and other structure parameters, leaf reflectance r L and other spectral parameters, and solar zenith angle θ sBased on the above-mentioned parameters, the canopy BRF of the first scattering process can be calculated by using the formula (17)~(25). Specifically, the path length of the illumination canopy component and the illumination soil component considering the influence of the adjacent pixel radiation is changed, and then the shielding factor of the illumination canopy component (as shown in formula 17~19) and the shielding factor of the illumination soil component (as shown in formula 21~23) are calculated, and then the new illumination canopy component (as shown in formula 20) and the new illumination soil component (as shown in formula 24) considering the influence of the adjacent pixel radiation are derived, and the sum of the two is the target pixel BRF (as shown in formula 25).

[0156] 4. Calculate the multiple scattering BRF of each target pixel only considering the internal scattering process of the target pixel self_mul : Using the structural parameters such as leaf area volume density FAVD, G function, canopy height model CHM, spectral parameters such as leaf reflectance r L , observation geometry parameters such as solar zenith angle θ s , and the intermediate variable leaf area index LAI (as shown in formula 31) obtained from the first scattering process, the multiple scattering BRF of the target pixel only considering the internal scattering process of the target pixel and not considering the influence of the adjacent pixel radiation can be calculated based on formula (26)~(37). Specifically, the leaf multiple scattering contribution can be calculated based on formula (26)~(30) and formula (32), the soil multiple scattering contribution can be calculated based on formula (33)~(36), and the sum of the two is the target pixel BRF (as shown in formula 37).

[0157] 5. Calculate the multiple scattering BRF of each target pixel considering the influence of the adjacent pixel adj_mul : Using the structural parameters such as leaf area volume density FAVD, G function, canopy height model CHM, spectral parameters such as leaf reflectance r L , observation geometry parameters such as solar zenith angle θ s , and the intermediate variable leaf area index LAI (as shown in formula 31) obtained from the first scattering process, the multiple scattering BRF of each target pixel considering the influence of the adjacent pixel can be calculated based on formula (38)~(39). Specifically, first, the multiple scattering BRF of the target pixel only considering the internal scattering process of the target pixel is calculated, and then the contribution of the adjacent pixel to the multiple scattering of the target pixel is calculated by using the neighborhood convolution algorithm (as shown in formula 38), and then the sum of the two is the target pixel BRF (as shown in formula 39).

[0158] 6. Calculate the total BRF only considering the internal scattering process of the target pixel: the canopy first scattering BRF contribution BRF self_sgl of the target pixel only considering the internal scattering process of the target pixel and the canopy multiple scattering BRF contribution BRF self_mulThe total BRF of each target pixel affected by the neighboring pixels is calculated by adding the canopy single scattering BRF contribution BRF

[0159] 7. The total BRF of each target pixel affected by the neighboring pixels is calculated by adding the canopy single scattering BRF contribution BRF adj_sgl and the canopy multiple scattering BRF contribution BRF adj_mul The total BRF of each target pixel affected by the neighboring pixels is calculated by adding the canopy single scattering BRF contribution BRF

[0160] Modeling the BRF of high spatial resolution forest canopy remote sensing pixels and validation examples:

[0161] To test the performance of the APPLE-GO model under different conditions, the simulated results of the computer simulation model LESS and real remote sensing images were compared.

[0162] 1. Comparison with LESS

[0163] Three different forest canopy scenarios were designed in the computer simulation LESS. The first forest canopy scenario was composed of abstract tree crowns with a height of 5 m in a Poisson distribution. The second forest canopy scenario was composed of abstract tree crowns with a height of 10 m in a distribution closer to the actual distribution. The third forest canopy scenario was a hybrid scenario composed of tree crowns with a height of 5 m and 10 m in a Poisson distribution. In the experiment, for each forest canopy scenario, three spatial resolutions of remote sensing images were set, 2 m, 5 m, and 10 m. The solar zenith angle was set at intervals of 10°, ranging from 0° to 60°, and the solar azimuth angle was set at intervals of 45°, ranging from 0° to 360°. The specific experimental design information is shown in Table 1.

[0164] Table 1: Simulation scenario experimental design table

[0165]

[0166] Figure 2 The comparison between the model estimates and the LESS simulation values is shown in the comparison chart. In the comparison chart, the scatter points represent the results obtained by randomly sampling the pixels in the three virtual forest canopy scenarios under different conditions (sun incident direction, pixel spatial resolution).

[0167] To analyze the BRF variation of high spatial resolution pixel in the vertical observation condition under different solar incident directions, the pixels located at the edge of the scene in Scene 1 were selected, and the BRF calculated by the APPLE-GO model and the BRF simulated by the LESS model were plotted in polar coordinates as shown in Figure 3 .

[0168] 2. Comparison with real remote sensing image

[0169] The APPLE-GO model was applied to a larch forest plot in a certain area, and the spatial resolution of the multispectral remote sensing image was 1.2 m. After inputting the basic parameters of the remote sensing image, the crown structure parameters of the forest, the spectral parameters, and the observation geometry parameters into the APPLE-GO model, the calculated reflectance image was compared with the reflectance image of Beijing No. 3 as shown in Figure 4 (both were resampled to 3 m spatial resolution).

[0170] The comparison results with the simulated data of the LESS model show (as shown in Figure 2 ) that the BRF results calculated by the APPLE-GO model are consistent with the simulated values of the LESS model; the polar coordinate graph shows (as shown in Figure 3 ) that the variation of the calculated values of the APPLE-GO model in each direction is consistent with the simulated values of the LESS model. The BRF of the pixels near the left edge shows obvious differentiation in orientation after being affected by the adjacent pixels, which is well reflected in the simulation results of the APPLE-GO model. Therefore, the APPLE-GO model well describes the influence of adjacent pixels.

[0171] The comparison results with the above satellite remote sensing image show (as shown in Figure 4 ) that the results calculated by the APPLE-GO model are generally close to the reflectance of the Beijing No. 3 satellite.

[0172] In summary, the APPLE-GO model has high simulation accuracy for the BRF of high spatial resolution forest canopy remote sensing pixels.

[0173] Although specific embodiments of the present application are disclosed for illustrative purposes, the purpose is to help understand the content of the present application and to implement it, those skilled in the art can understand that various substitutions, changes and modifications are possible without departing from the spirit and scope of the present application and the appended claims. Therefore, the present application should not be limited to the disclosed content of the best mode, and the scope of the present application claimed is defined by the scope of the claims.

Claims

1. A method for calculating the two-way reflectance factor of high-resolution remote sensing images of forest canopy, comprising the following steps: 1) Determine the basic information of the remote sensing image and the corresponding canopy height model (CHM), including the image range, image spatial resolution, and CHM spatial resolution; collect forest canopy structure parameters and forest spectral parameters. The forest canopy structure parameters include leaf volume density (FAVD), G-function, and CHM. The spectral parameters include leaf reflectance (r). L Determine the observation geometric parameters and the proportion of sky-scattered light, wherein the observation geometric parameters include the solar zenith angle θ. s Sun azimuth Observation of zenith angle θ v and observation azimuth 2) The bidirectional reflectance factor (BRF), which considers only the primary scattering occurring within the target pixel, is calculated using the forest canopy structure parameters, the spectral parameters, and the observation geometric parameters. self_sgl Among them, according to BRF self_sgl =K g ·r g +K l ·r l The bidirectional reflectivity factor (BRF) for single-scattering was calculated. self_sgl ;K g K represents the proportion of soil area exposed to sunlight. l r represents the proportion of leaf area exposed to light. g r represents the reflectivity of soil under sunlight. l 21) Calculate the proportion of soil area illuminated by sunlight using a path length distribution model, K. g The ratio of the area of ​​soil exposed to sunlight, K g K from the large pores between the tree canopy g1 Small pores K in the shaded ground g2 Soil light spots K observed through pores in the tree canopy g3 Composition; where K g1 =P(l(Ω) s =0),l(Ω v =0)) represents the path length l (Ω) in the direction of solar incidence. v The path length l (Ω) in the sensor observation direction is 0. s The joint probability when K is 0; g2 =P(l(Ω) s >0), l(Ω v =0))·P(Ω s |l(Ω s >0), l(Ω v =0), P(l(Ω) s >0), l(Ω v =0)) represents the path length l (Ω) of the sensor observation direction. v When ) is 0 and the path length l (Ω) of the solar incident direction is 0 s The joint probability that Ω is greater than 0, P(Ω) s |l(Ω s >0), l(Ω v =0)) refers to calculating the porosity P(Ω) based on the area ratio of shaded soil and the distribution of path length. s |l(Ω s >0), Where G(Ω) s () represents the projection ratio of the leaf towards the direction of solar incidence, FAVD is the leaf area volume density, l g2 (Ω s p represents the path length of the shaded soil portion along the direction of solar incidence. l_g2 (Ω s ) represents the path length l g2 (Ω s The probability density distribution of ), l max l min They represent path length l g2 (Ω s The maximum and minimum values ​​of K; g3 =P(l(Ω) v >0))·P(Ω v |l(Ω v >0)), P(l(Ω) v >0))·P(Ω v |l(Ω v >0)) is the path length l(Ω) in the observation direction. v The proportion of porosity greater than 0. l g3 (Ω v ) represents the path length of the canopy portion in the vertical observation direction, p l_g3 (Ω v ) represents the path length l g3 (Ω v The probability density distribution of G(Ω) v 2) The proportion of the leaf's projection in the vertical observation direction; 22) The proportion of the illuminated leaf area K l Bidirectional reflectivity factor (BRF) l =BRF lc +BRF lt By calculating the area ratio K of the illuminated tree canopy. c The ratio of the area of ​​the shaded canopy to K t The proportion of the tree canopy area illuminated by sunlight, K c Bidirectional reflectance factor (BRF) of the light-reflecting leaf lc The area ratio K of the shaded tree canopy t Bidirectional reflectance factor (BRF) of the light-reflecting leaf lt Among them, the proportion of tree canopy area illuminated by sunlight K c =P(l(Ω) s =0),l(Ω v >0))· [1-P(Ω v |l(Ω s =0),l(Ω v >0))], the area ratio of the shaded canopy K t =P(l(Ω) s >0), l(Ω v >0))·[1-P(Ω v |l(Ω s >0), l(Ω v >0))];P(l(Ω s =0),l(Ω v >0)) refers to the joint probability that the part with a path length greater than 0 in the observation direction has a path length of 0 in the solar incidence direction; 3) Using the forest canopy structure parameters, the spectral parameters, and the observation geometry parameters, calculate the first-order scattering two-way reflectance factor (BRF) of each target pixel as a result of the influence of neighboring pixels. adj_sgl First, the shading factor S of the light canopy component is calculated. c Shading factor S of soil light content g Then, the new illumination canopy component K′ considering the influence of neighboring pixel radiation is derived. c and the soil component K′ of new sunlight g Then, the bidirectional reflectance factor (BRF) of the primary scattering of each target pixel due to the influence of neighboring pixels is calculated. adj_sgl The specific method is as follows: 31) Calculate the proportion P of large pores formed on the tree canopy after being blocked by neighboring pixels. c (l=0,Ω s )=P(F c (l=0,Ω s )|K c ), the proportion of small apertures formed after being blocked by neighboring pixels Then calculate the shading factor S on the tree canopy. c =[1-P c (l=0,Ω s )]·[1-P c (l>0,Ω s Then calculate the area ratio K of the target pixel's canopy component after being affected by radiation from neighboring pixels. c ′=K c ·(1-S c );l(Ω s |K c p represents the path length of sunlight along the direction of incidence on the original illuminated part of the tree canopy of the target pixel after being affected by neighboring pixels. l (Ω s |K c ) represents the path length l (Ω) s |K c ) probability density distribution; 32) calculate the proportion P of large pores formed on the soil after being blocked by neighboring pixels. g (l=0,Ω s )=P(F g (l=0,Ω s )|K g ), the proportion of small apertures formed after being blocked by neighboring pixels and the shading factor S on the soil g =[1-P g (l=0,Ω s )]·[1-P g (l>0,Ω s Then, calculate the area ratio K of the target pixel's irradiated soil component after being affected by the radiation of neighboring pixels. g ′=K g ·(1-S g );l(Ω s |K g p represents the path length of the sun's incident direction in the original sunlight-illuminating part of the target pixel after being affected by neighboring pixels. l (Ω s |K g ) represents the path length l (Ω) s |K g ) probability density distribution; 33) calculate the first scattering bidirectional reflectivity factor BRF. adj_sgl = (1-β)·(K′) g ·r s +BRF′ l ); β represents the proportion of sky-scattered light; 4) Using the forest canopy structure parameters, the spectral parameters, the observation geometric parameters, and the leaf area index (LAI), an intermediate parameter obtained from the first scattering process, calculate the bidirectional reflectance factor (BRF) considering only the area within the target pixel after multiple scattering. self_mul Among them, the first step is to determine the ratio of the illuminated canopy area K. c , the proportion of shaded canopy area K t and the vertical observation path length distribution F(l>0,Ω) v ), calculated Then, based on the probability p of a photon scattering from a leaf and then colliding again with a leaf inside the canopy, 1 The probability of multiple collisions p m The calculation ignores the soil and only considers the bidirectional reflectivity factor (BRF) of multiple scattering between the canopy layers. self_vmul And considering the bidirectional reflectance factor (BRF) between the soil and the interior of the canopy. self_smul Then, the canopy multiple scattering bidirectional reflectivity factor (BRF) was calculated when only the target pixel was considered. self_mul =BRF self_vmul +BRF self_smul ; 5) Using the forest canopy structure parameters, the spectral parameters, the observation geometric parameters, and the leaf area index (LAI), an intermediate parameter obtained from the first scattering process, calculate the bidirectional reflectance factor (BRF) of each target pixel affected by neighboring pixels. adj_mul ; 6) The canopy primary scattering bidirectional reflectivity factor (BRF) will only consider the scattering process within the target pixel. self_sgl The canopy multiple scattering bidirectional reflectivity factor (BRF) that only considers the scattering process within the target pixel. self_mul Adding them together yields the total bidirectional reflectivity factor (BRF) that only considers the scattering process within the target pixel. self ; 7) The canopy primary scattering bidirectional reflectivity factor (BRF) will take into account the influence of neighboring pixel radiation. adj_sgl and the canopy multiple scattering bidirectional reflectivity factor (BRF) considering the influence of neighboring pixel radiation. adj_mul The summation yields the total bidirectional reflectance factor (BRF) of each target pixel, reflecting the influence of neighboring pixels. adj .

2. The method according to claim 1, characterized in that, In step 5), the bidirectional reflectivity factor (BRF) of multiple scattering caused by neighboring pixels is calculated for each target pixel. adj_mul The method is as follows: First, based on the ratio of the canopy area illuminated by sunlight, K... c , the proportion of shaded canopy area K t and the vertical observation path length distribution F(l>0,Ω) v ), calculated Then, the bidirectional reflectivity factor of multiple scattering from neighboring pixels to the target pixel is calculated. Among them, LAI i LAI represents the leaf area index of the i-th surrounding pixel. self Leaf area index (BRF) representing the current target cell mul_i W represents the multiple scattering BRF calculated for the i-th surrounding pixel when only internal multiple scattering is considered. i The distance-dependent normalized weight matrix is ​​then used; the multiple scattering bidirectional reflectivity factor (BRF) is then calculated. adj_mul =BRF self_mul +BRF others_mul .

3. The method according to claim 1, characterized in that, p 1 =0.7exp(k1·LAI)-0.66exp(k2·LAI), Where, k1 = 0.0045exp(1.2555cosθ), k2 = 0.1982lncosθ - 0.7146, i(θ) is the interception probability of the canopy in the direction of θ, where θ is the solar zenith angle.

4. The method according to claim 3, characterized in that, Where i0=β·i D +(1-β)·i S ω represents the average probability of the canopy intercepting direct sunlight and sky-scattered light, and r represents the leaf primary scattering albedo. s R represents soil reflectance. dn T represents the albedo at the base of the canopy. dn and T up p(Ω) represents the descending and ascending transmittance of the canopy, respectively. s ) and p(Ω v ) represent the porosity in the direction of solar incidence and the direction of sensor observation, respectively.

5. A server, characterized in that, It includes a memory and a processor, the memory storing a computer program configured to be executed by the processor, the computer program including instructions for performing each step of the method of any one of claims 1 to 4.

6. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the method according to any one of claims 1 to 4.

Images