Quantitative analysis method of solar shortwave radiation in tree canopy based on computer simulation model

Abstract

The present invention discloses a quantitative analysis method of solar shortwave radiation of a tree canopy based on a computer simulation model, comprising: first, collecting tree laser point cloud data, combining machine vision algorithm to separate branches and leaves and segment single leaves on the tree laser point cloud, and using elliptical and cylindrical geometric units of appropriate size to fit each leaf and skeleton respectively to reconstruct a real model of the tree; secondly, using computer graphics methods, combining the longitude and latitude of the research site and time to simulate the incident light of the sun, and introducing the bidirectional reflection and transmission distribution function of physics and the Monte Carlo ray tracing algorithm to simulate the collision of reflected and transmitted light with leaves in the canopy; finally, according to the instrument measurement, the average roughness and refractive index of leaves of different tree species are obtained, and the shortwave radiation distribution calculation in the tree canopy and the light interception efficiency evaluation of the forest canopy are realized in combination with the ray tracing algorithm. The present invention can accurately obtain the interception amount of direct, reflected and transmitted light energy by the tree canopy.

Description

Technical Field

[0001] The invention relates to the technical field of research on radiation flux in a forest canopy, and in particular to a quantitative analysis method of solar shortwave radiation in a forest canopy based on a computer simulation model. Background Art

[0002] The wavelength of solar shortwave radiation ranges from 300nm to 3000nm, accounting for about 97.5% of the total solar radiation energy, and there is a relatively stable conversion ratio (about 0.45-0.46) between photosynthetically active radiation (PAR) and solar shortwave radiation. Forests are an important ecological environment component and material resource, and solar shortwave radiation is of great significance to scientific issues such as the study of asymmetric light competition of forests, spatial allometric growth phenomena, and light interception capabilities of different canopies.

[0003] When conducting quantitative analysis of forests using optical remote sensing, the methods for studying the interaction mechanism between the canopy and solar shortwave radiation are usually divided into several categories, such as geometric optics models, radiation transfer models, and computer simulation models. The geometric optics model calculates the radiance by calculating the area ratio difference of the four components, namely, visible illumination target, visible illumination background, visible shadow target, and visible shadow background, by emitting light and vegetation structure parameters of the canopy. Due to the complex geometric characteristics of the canopy as a non-rigid object and the limitations of data and computing power on the light transmission model, this type of model usually simplifies and abstracts the description of trees when modeling. For example, the canopy is simplified as an ellipsoid, and the canopy is regarded as a chaotic medium to study the influence of the canopy shape, density, and center thickness on the bidirectional reflectance distribution function (BRDF) of the canopy under illumination. The radiation transfer model takes into account multiple scattering through the radiation transfer equation and is often used in continuous canopy models. For example, the correlation between radiation and leaf area index is established based on the Beer-Lambert-Bouguer law. By stratifying or voxelizing the canopy height, the gradual attenuation of radiation in the canopy and the quantitative inversion of the radiation transfer flux in the canopy are reflected from the calculation of the model equation, or the light interception results produced by different leaf inclination angles and azimuth angles are studied at the scale of a single plant model. At present, there are also improvements to existing models such as the 4SAILT model. In recent years, computer simulation models have also been applied to the analysis of complex forest lighting. They are usually combined with L systems, spatial rendering of plant three-dimensional scenes, and hierarchical heterogeneous canopy porous media models to simulate the real forest canopy structure and carry out light tracing calculations in forest scenes. Ray tracing can usually be divided into two methods: forward ray tracing and reverse ray tracing. Forward ray tracing requires tracing all rays from the light source, as well as the reflected and transmitted rays generated by their intersection with the leaves in the canopy, until the light intensity is small enough or escapes the scene. Reverse ray tracing is based on the principle of ray reversibility. It traces the rays from the viewpoint in reverse, and calculates the reflected and transmitted rays generated after the intersection with the leaves in the canopy according to virtual information such as material, texture, and scene illumination, until the rays reach the light source, so as to calculate the incident radiance at the point of closest collision between the rays and the surface of the object. Typical methods include: the LESS model uses the forward ray tracing method to simulate multi-spectral radiation-related data, and uses the reverse ray tracing method to generate fisheye images or large-scale spectral images. The FluorWPS model uses the forward ray tracing algorithm based on the Monte Carlo law of large numbers to simulate the radiation transfer of sunlight-induced chlorophyll fluorescence (SIF) in the three-dimensional canopy, and can accurately quantify the contribution of each plant part in the radiation process; or it depicts the distribution of leaf inclination angle and leaf density at the voxel scale to accelerate the computational efficiency of the solar radiation transmission model in the canopy; or it uses the box drawing method instead of the traditional Monte Carlo roulette method to calculate the probability of transmission, reflection and scattering of light, thereby reducing the amount of calculation.

[0004] Although there are many mature methods in the study of three-dimensional radiation transfer models of forests, the following problems still exist: 1. Due to the extremely complex three-dimensional topological structure of the tree crown, many voxel-based and simulation data cannot fully describe the details of the tree crown, which makes it impossible to show the fine depiction of light in the tree crown; 2. Laser point clouds can better present the expression morphology and phenotypic characteristics of trees. How to combine forest point cloud data with light transmission methods to collaboratively invert the light radiation value in the tree crown and its distribution in the canopy still needs in-depth research; 3. The leaves in the tree crown are subject to direct illumination, reflection and transmission under the lighting environment. How to effectively present these variables related to the phenotypic characteristics of trees using computer graphics technology and carry out parallel world replication of forest lighting scenes from the perspective of digital twins is a current research hotspot.

[0005] Based on the above problems, it is necessary to study a new method to obtain the distribution of solar shortwave radiation flux in the canopy at different times, and then obtain the amount of direct, reflected and transmitted light energy intercepted by the forest canopy. Summary of the invention

[0006] The technical problem to be solved by the present invention is to provide a quantitative analysis method for solar shortwave radiation of tree canopy based on computer simulation model in view of the shortcomings of the above-mentioned prior art. The quantitative analysis method for solar shortwave radiation of tree canopy based on computer simulation model combines laser point cloud and computer graphics content, adopts pattern recognition algorithm to separate branches and leaves and identify single leaves on tree laser point cloud, and carries out three-dimensional reconstruction of trees based on the acquired tree parameters, tracks light at a fine scale, obtains the distribution of solar shortwave radiation flux in the canopy at different times, and then accurately obtains the interception amount of direct, reflected and transmitted light energy by the tree canopy.

[0007] In order to achieve the above technical objectives, the technical solution adopted by the present invention is:

[0008] The quantitative analysis method of solar shortwave radiation of tree canopy based on computer simulation model includes the following steps:

[0009] Step 1: Collect forest stand point cloud data;

[0010] Step 2: Classify all tree point cloud data in the forest stand into leaf point cloud data and branch point cloud data; perform single leaf segmentation and triangulation on the leaf point cloud data, and then represent each leaf as a set of triangular facets;

[0011] Step 3: Calculate the surface solar irradiance at the forest stand location;

[0012] Step 4: Model the sunlight as parallel rays emitted from above the trees, and obtain the area of ​​a current triangle directly illuminated by the sunlight by intersecting the straight line with the triangle.

[0013] Step 5: Calculate the solar radiation flux of each incident direction vector through the number of incident parallel rays and the surface solar irradiance, calculate the incident radiance of the incident solar ray through the solar radiation flux, the area of ​​the current triangular facet obtained, and the angle between the incident solar ray direction vector and the normal vector of the triangular facet directly illuminated by the incident solar ray, and calculate the direct radiation flux intercepted by multiple triangular facets;

[0014] Step 6: Determine the BRDF value through the direction vector of the incident sunlight, the normal vector of the triangle directly hit by the incident sunlight, the direction vector of the reflected light, and the normal vector of the reflected microplane, and then obtain the reflected radiance of the incident sunlight. After the second intersection, obtain the other triangles that intercept the reflected light, and calculate the reflected radiation flux intercepted by other triangles;

[0015] Step 7: Determine the BTDF value through the direction vector of the incident sunlight, the normal vector of the triangle directly hit by the incident sunlight, the direction vector of the transmitted light, and the normal vector of the transmitted microplane, and then obtain the transmitted radiance of the incident sunlight. After the second intersection, obtain the other triangles that intercept the transmitted light, and calculate the transmitted radiation flux intercepted by other triangles.

[0016] Step 8: Traverse all triangular facets and calculate the solar radiation flux that is ultimately intercepted by all triangular facets within the tree canopy.

[0017] As a further improved technical solution of the present invention, the blade point cloud data is triangulated in step 2, specifically: the point cloud data of each blade is triangulated using Delaunay triangulation, and each blade is represented as a set of triangular facets t l,s ∈T l , where t l,s It is the sth triangle in the lth leaf.

[0018] As a further improved technical solution of the present invention, the step 3 is specifically:

[0019] Surface solar irradiance E at the forest stand location total for:

[0020] E total =E direct +E diffuse (1);

[0021] The direct solar irradiance is E direct =E sc d m τ r τ o τ g τ w τ a sinθs ; E sc is the solar constant; d m is the eccentricity correlation coefficient of the Earth's orbit; τ r , τ o , τ g , τ w and τ a represent the scattering transmittances of Rayleigh, ozone, gas, water, and aerosol, respectively;

[0022] θ s is the local solar altitude angle, expressed as follows:

[0023]

[0024] in, is the local latitude, δ s is the local solar declination angle, τ s is the local solar hour angle, is the incident sun ray vector;

[0025] Horizontal diffuse irradiance E diffuse for:

[0026] E diffuse =E dr +E da +E dm (3);

[0027] Among them, E dr is the Rayleigh scattering of sunlight after entering the atmosphere, E da is the aerosol scattering of sunlight after it enters the atmosphere, E dm It is multiple reflections from the ground and the sky.

[0028] As a further improved technical solution of the present invention, the step 4 is specifically as follows:

[0029] Model the sunlight as parallel rays emitted from above the trees, and determine the sth triangle t in the lth leaf of the canopy point cloud construction. l,s The first intersection, triangle patch t l,s The area is a l,s .

[0030] As a further improved technical solution of the present invention, the step 5 is specifically as follows:

[0031] Calculate the number of solar rays N emitted per square meter based on the fixed interval between adjacent parallel rays inc , then the solar radiation flux of each incident direction vector is:

[0032]

[0033] The incident radiance of the incident sunlight is:

[0034]

[0035] When the incident ray intersects the triangle, it means that there are n inc The incident light beam is reflected by the current triangle patch t of the canopy leaf l,s Intercept; where θ i is the direction vector of the incident sunlight i and the triangle t that intersects the incident sun ray l,s Normal vector ω n The angle between Direct radiation flux intercepted by canopy leaves Where S inc It is the sum of the areas of the triangular faces that intercept the incident sunlight.

[0036] As a further improved technical solution of the present invention, the step 6 is specifically as follows:

[0037] The BRDF values ​​are:

[0038]

[0039] in:

[0040]

[0041]

[0042]

[0043] where θ r is the triangle patch normal vector ω n and the reflected light direction vector ω r The angle between n is the triangle facet that intersects with the incident sun ray l,s Normal vector; ω i is the direction vector of the incident sunlight; θ h is the reflection microplane normal vector ω h and the triangle patch normal vector ω n The angle between h =ω h ·ω n ; α is the roughness of the microsurface; χ + The function is a step function; c = |ω r ·ω h |, η t is the blade refractive index; the reflection microfacet normal vector is defined as

[0044] By determining the direction vector ω of the incident sun ray i The triangular patch t formed by the point cloud of leaves in the crown l,s The first intersection, when the light intersects with the triangular face, indicates that a beam of light is intercepted by the canopy leaves; from the first intersection of the triangular face t l,s The reflected radiance of each incident ray emitted from By performing random importance sampling (Qi et al., 2022) on discrete points on the surface of the reflection lobe pattern, M reflection directions are selected, and the reflection radiance of the mth light is expressed as follows, where m = 1, 2, ..., M:

[0045]

[0046] Among them, I inc (ω i )=I inc ;

[0047] The value of M is 10, and the incident light and the triangle t l,s After the reflected light emitted by the intersection intersects with other triangular facets, the intercepted reflected radiation flux F ref It can be expressed as:

[0048]

[0049] Where S ref Represents the sum of the areas of the triangular patches that intercept the reflected light, θ ref is the direction vector of the reflected light r Along its direction and other normal vectors is ω n The intersecting triangles t l,s The angle between

[0050] As a further improved technical solution of the present invention, the step 7 is specifically as follows:

[0051] The BTDF values ​​are:

[0052]

[0053] Where: η i is the refractive index of the medium on the incident light side, η o is the refractive index of the medium on the side of the transmitted light, ω t is the direction vector of the transmitted light, θ t is the transmitted light direction vector ω t and the triangle normal vector ω of the leaf surface that intersects with it n The angle betweent =ω t ·ω n ; The transmission microfacet normal vector is defined as

[0054]

[0055]

[0056]

[0057] Triangle patch t l,s The transmitted radiance emitted from above is recorded as: By performing random importance sampling on the discrete points on the surface of the left transmission lobe diagram and selecting M directions, the transmission radiance of the mth light is expressed as follows:

[0058]

[0059] Incident light and triangle t l,s If the transmitted light emitted after the intersection is intercepted by other triangular facets, the intercepted transmitted radiation flux F tra It can be expressed as:

[0060]

[0061] Where S tra Represents the sum of the areas of the triangular patches that intercept the transmitted light, θ tra is the transmitted light direction vector ω t With other triangles t l,s After intersection, ω t and the intersecting triangle face normal vector ω n The angle between In formula (17), the area of ​​the triangular facet that intercepts the transmitted light is recorded as a l,s .

[0062] As a further improved technical solution of the present invention, the step 8 is specifically as follows:

[0063] S inc The sum of the areas of the triangles that intercept the incident light. According to the ray tracing algorithm, all triangles t are traversed. l,s ∈T l , the solar radiation flux intercepted by all triangular facets in the canopy is calculated as:

[0064]

[0065] The beneficial effects of the present invention are:

[0066] The present invention reconstructs the real apparent morphological structure of trees from laser point clouds, integrating the multidisciplinary perspectives of surveying and mapping science, computer graphics, and forestry, and integrates the Monte Carlo ray tracing algorithm with the bidirectional reflection and transmission distribution function of physics, to truly simulate the propagation process of light in the tree crown. Field verification has shown that the method of the present invention can accurately obtain the quantitative values ​​of the distribution of solar shortwave radiation in the tree canopy and the interception of the canopy radiation flux. The present invention integrates the digital design concept of forest vegetation communities based on laser point clouds, with the help of computer simulation models of computer science, to achieve a detailed characterization of the distribution of solar shortwave radiation in the tree canopy. The work of the present invention has important implications for studying the absorption of solar radiation by trees under temporal and spatial changes, the coupled correlation between light intensity and canopy morphological structure, and canopy radiation transmission models under different phenotypic parameters. The present invention combines laser point cloud with computer graphics content, adopts pattern recognition algorithm to separate branches and leaves and identify single leaves on tree laser point cloud, and carries out three-dimensional reconstruction of trees based on the acquired tree parameters, and tracks light at a fine scale to obtain the distribution of solar shortwave radiation flux in the canopy at different times, thereby accurately obtaining the interception amount of direct, reflected and transmitted light energy by the tree canopy. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 (a) shows the reflection and transmission diagram of incident sunlight in the entire sphere Ω;

[0068] Figure 1 (b) shows the schematic diagram of the reflection lobe derived from BRDF;

[0069] Figure 1 (c) shows the schematic diagram of the transmission lobe derived from BTDF;

[0070] Figure 2 The solar shortwave radiation irradiance at each moment from 7:00 to 17:00 on the 15th day of each month in April, June, August, and October 2020 under clear sky conditions calculated by the algorithm in this paper is at the experimental plot.

[0071] Figure 3 The polar coordinate diagram shows the calculation results of the bidirectional reflectance distribution function (BRDF) (first row) and the bidirectional transmission distribution function (BTDF) (second row) corresponding to the incident light at different angles.

[0072] Figure 4 Taking the mango tree (left) and the cherry tree (right) as examples, the computer simulation method in this article shows the direct sunlight in the tree canopy at 11:00 on August 15, 2020, and the random reflection and transmission light paths generated based on Monte Carlo.

[0073] Figure 5This is a visualization of the distribution of direct, reflected, and transmitted irradiance in the canopy of four experimental trees at 11:00 on August 15, 2020, as demonstrated by the method used in this article.

[0074] Figure 6 This is a calculated diagram and a local enlarged image of the solar irradiance distribution in the canopy of multiple camphor trees at 11:00 on August 15, 2020.

[0075] Figure 7 (a) is the radiation flux value of the sun intercepted by the mango tree calculated by the method in this paper as it changes over time.

[0076] Figure 7 (b) is the radiation flux value of the sun intercepted by the rubber tree calculated by the method in this paper as it changes over time.

[0077] Figure 7 (c) in the middle is the radiation flux value of the sun intercepted by the crape myrtle tree calculated by the method in this paper as it changes over time.

[0078] Figure 7 (d) in the middle is the radiation flux value of the sun intercepted by the cherry tree calculated by the method in this paper as it changes over time.

[0079] Figure 7 (e) is the radiation flux value of the sun intercepted by a camphor tree forest calculated by the method in this paper as it changes over time. DETAILED DESCRIPTION

[0080] The specific embodiments of the present invention are further described below according to the accompanying drawings:

[0081] Solar shortwave radiation has an important impact on the photosynthesis and growth and development of trees. Generally, the canopy structure of trees and the distribution and optical properties of leaves within the canopy affect the absorption and interception of solar shortwave radiation. Quantitative simulation and characterization of the distribution and interception of solar shortwave radiation within the canopy of trees, and inversion of the changes in the radiation flux of solar shortwave radiation in different canopies of trees under temporal and spatial changes, provide a theoretical basis for forest cultivation, management, quality improvement and efficiency enhancement.

[0082] This paper provides a quantitative analysis method for solar shortwave radiation in the tree canopy based on a computer simulation model. Specifically, it combines laser point cloud with computer graphics content, uses a pattern recognition algorithm to separate branches and leaves and identify single leaves on the tree laser point cloud, and based on the acquired tree parameters, carries out three-dimensional reconstruction of the trees, tracks the light at a fine scale, and obtains the distribution of solar shortwave radiation flux in the canopy at different times, and then calculates the interception of direct, reflected and transmitted light energy by the tree canopy. This work has an enlightening effect on the irradiance inversion and light transmittance analysis in the spatiotemporal changes within the plant canopy from a multidisciplinary perspective.

[0083] 1. Materials and data:

[0084] 1.1. Virtual trees and plot data:

[0085] The dataset of this paper is composed of two parts: public dataset and laser scanner acquisition dataset. The specific data can be found in Table 1. In the first part of the public dataset, the skeleton structure and leaf space information of two trees are obtained by manual measurement, and the virtual tree data are reconstructed. The two trees include a mango tree (Mangifera indica) and a rubber tree (Hevea brasiliensis), and the crowns of the two trees are independent of each other and do not intersect with other crowns. The mango tree comes from an orchard in Ban Bong District, Chonburi Province, Thailand (13.17°N, 101.25°E), and the rubber tree is located in the garden of the Faculty of Agriculture, Kasetsart University, Thailand (13.85°N, 100.75°E). The azimuth and inclination angles of the midrib of the leaf veins on each leaf, and the rotation angle of the leaf relative to the midrib of the leaf veins, etc., were collected using a three-dimensional digitizer.

[0086] In addition, two field tree data were collected at Nanjing Forestry University (32.07°N, 118.81°E), including a crape myrtle tree (Lagerstroemia indica) and a cherry tree (Prunus subg.Cerasus), and a camphor tree (Camphora officinarum) stand was selected in the back mountain. Nanjing Forestry University is located in Nanjing, the capital city of Jiangsu Province. Nanjing has a warm climate, four distinct seasons, and abundant rainfall. It belongs to the warm and humid subtropical monsoon climate. The average annual precipitation here is 1000 to 1600 mm, and the average annual temperature is 15.4℃. The specific tree information is shown in Table 1:

[0087] Table 1, detailed information of experimental trees and plots:

[0088]

[0089] The laser point cloud data of campus trees was collected using the FARO Focus3D X 330 laser scanner, with a maximum scanning distance of 330m, a scanning rate of up to 976,000 points / second, a resolution of more than 70 million color pixels, a scanning wavelength of 1550nm, and a vertical and horizontal field of view of 180° and 360° respectively. The data collected in this experiment had a scanning resolution of 2-3mm dot spacing within a range of 10 meters. During point cloud collection, the scanner was placed about 5 meters away from the tree. With the tree as the center, the scanner collected and spliced ​​data from two stations at two symmetrical positions; for the camphor tree plot data collection, the point cloud data of the entire forest stand was obtained by splicing multi-station scanning data at different locations in the forest. The height of the scanning instrument was about 1m.

[0090] 1.2. Tree 3D reconstruction and triangulation:

[0091] The point cloud data obtained from the laser scanner cannot be used directly. On the one hand, the point cloud is unevenly distributed, that is, there are redundant data, there are holes formed by occlusion, and the existence of noise in the data will also affect the calculation. First, the cloth simulation filtering algorithm is used to filter out the ground points in the lidar point cloud data. At the same time, in order to eliminate outliers, the point cloud data needs to be denoised.

[0092] In order to reconstruct the tree point cloud data in three dimensions, the branch-leaf separation algorithm is used to classify the tree point cloud data into leaf and branch point cloud data. The leaf point cloud data contains a large number of leaves, so a single leaf segmentation algorithm is needed to achieve automatic segmentation of the leaves. Single leaf segmentation uses a three-dimensional point cloud segmentation method, uses a spherical neighborhood model to obtain the center point of the leaf, and selects the DBSCAN algorithm to cluster the center point based on the distribution characteristics of the center point, so as to obtain the center point corresponding to each leaf point cloud surface. Finally, a leaf-leaf separation algorithm based on a three-dimensional watershed is used to achieve individual segmentation of the leaf.

[0093] Because the zenith angle, azimuth, length, width and spatial coordinates of each leaf have been obtained, when a single tree is modeled in 3D, each leaf of the canopy can be approximately expressed as an ellipse in space, with the major axis and minor axis of the ellipse corresponding to the length and width of the leaf, respectively. Initially, the ellipse is in the XOY plane, and the point set of the leaf is obtained by uniformly sampling within the ellipse. In this paper, the sampling interval is set to 0.4 cm, and the set of all points on the surface of the lth leaf is P l ′, the point set after restoring the blade posture can be obtained by the known rotation of the lth blade relative to the X-axis, Y-axis, and Z-axis as follows:

[0094]

[0095] in, is the rotation matrix, C l is the spatial coordinate of the center of the lth leaf. After the spatial rigid transformation of formula (1), the true shape of each leaf and the corresponding point cloud set P are obtained. l , use Delaunay triangulation to analyze the point cloud data P of each leaf l Triangulate, and then represent each leaf as a set of triangular patches t l,s ∈T l , where t l,s is the sth triangular patch in the lth leaf. The branches are reconstructed by a geometric model composed of a series of cylinders, where the bottom diameter of the cylinder is the average diameter of the corresponding branches. The branches above the crown have a smaller diameter, while the branches below the crown have a larger diameter. The sides of these cylinders are uniformly sampled, and each branch cylinder is converted into a discrete point. The point set obtained is Then perform Delaunay triangulation to obtain the branch triangle patch T b Finally, the tree triangle face set is T = {T l ,T b}.

[0096] 2. Ray tracing model:

[0097] 2.1. Calculation of surface solar irradiance in the experimental area:

[0098] The solar irradiance on the ground surface is affected by factors such as the solar altitude angle, ground altitude, and weather. This paper uses the Iqbal model to estimate the surface solar irradiance E of the experimental plot. total , unit is W·m -2 , E total Direct solar irradiance E direct and horizontal diffuse irradiance E diffuse Combined to get:

[0099] E total =E direct +E diffuse (2);

[0100] Direct solar irradiance E direct =E sc d m τ r τ o τ g τ w τ a sinθ s , E sc is the solar constant, which can be taken as 1367 W·m -2 .d m is the eccentricity correlation coefficient of the Earth's orbit, τ r , τ o , τ g , τ w , τ a represents the scattering transmittance of Rayleigh, ozone, gas, water and aerosol respectively, θ s is the solar altitude angle at the observation site. Horizontal diffuse irradiance E diffuse It can be expressed as: E diffuse =E dr +E da +E dm , where E dr and E da is the Rayleigh scattering and aerosol scattering of sunlight after it enters the atmosphere, E dm It is multiple reflections from the ground and the sky.

[0101] 2.2. Modeling of incident sunlight:

[0102] After obtaining the surface solar irradiance of the research plot, the distribution of solar radiation in the tree canopy can be studied based on the interaction between the simulated solar beam and the leaf surface. Since the distance between the sun and the earth reaches 150 million kilometers, the sunlight can be approximately modeled as parallel light, which will be uniformly distributed in a series of planes above the trees. inc The source point on is emitted, and the coordinates of the source point are marked as N inc Represents the total number of source points. Each source point emits a sun ray, and the sun incident ray vector can be obtained from formulas (3) and (4): and the solar altitude angle θ at a certain moment s and azimuth ω s Corresponding.

[0103] The solar altitude angle refers to the angle between the sun's rays at a certain place and the tangent line of the earth's surface perpendicular to the center of the earth. The calculation formula and the corresponding light vector are:

[0104]

[0105] in is the latitude of the experimental site, δ s is the experimental declination, τ s is the solar hour angle of the experimental site. At the same time, the solar altitude angle θ s The solar azimuth angle ω can be calculated s And the corresponding relationship with the light vector is:

[0106]

[0107] At the same time, for a given sun ray vector Can obtain the plane where the source point is located Π inc The plane equation is:

[0108]

[0109] here The uniform setting is 100m, which is higher than the height of all trees. At the same time, the point formula of the straight line simulating the sun's rays emitted from the source point k can be expressed as:

[0110]

[0111] The fixed interval between adjacent parallel rays is set to 0.2 cm. At this time, the number of solar rays emitted is N. inc 2.5*10 5 per square meter, combined with formula (2) to calculate E total Each incident direction vector The solar radiation flux is expressed as The unit is W.

[0112] By determining the incident light and the s-th triangle t in the l-th leaf constructed by the crown point cloud l,s The first intersection, triangle patch t l,s The area is The incident radiance of the sun's rays is The unit is W·m -2 ·sr -1 When the ray intersects a triangle, it indicates that there are n inc The incident ray is reflected by the current triangle patch t of the canopy leaf. l,s Intercept. Among them, θ i is the direction vector of the incident sun light ω i and the intersecting triangle t l,s Normal vector ω n The angle between Direct radiation flux intercepted by the canopy Where S inc It is the sum of the areas of the triangular patches that intercept the incident sunlight.

[0113] 2.3. Calculation of reflected radiation flux from leaves in the tree crown:

[0114] The bidirectional reflectance distribution function (BRDF) defines the incident direction ω i Irradiance E on inc (ω i ) and the reflection direction ω r The radiance I ref (ω r The relationship between the radiance differential dI of light emitted in any direction on the hemispherical surface Ω ref (ω r ) and the incident light radiance differential dE in any direction inc (ω i ), the BRDF function f r (ω r ,ω i ) can be expressed as:

[0115]

[0116] like Figure 1 As shown, Figure 1 It represents the visualization of the bidirectional reflection and transmission distribution function based on the Monte Carlo method in this paper; Figure 1 (a) in the figure shows the reflection and transmission diagram of the incident sunlight in the entire sphere Ω. Figure 1(b) in the figure shows the reflection lobe derived from the BRDF. After sampling the lobe surface at equal intervals, some points (blue) are randomly selected according to the importance sampling strategy and the corresponding reflected light (blue arrow) is generated. Figure 1 (c) in the figure shows the equivalent diagram of the BTDF. The spherical Ω diagram is used to show the normal vector ω of the illuminated triangle. n , the incident light direction vector ω i , the reflected light direction vector ω r , the reflection microplane normal vector ω h and the reflected light direction vector ω r .

[0117] According to the microplane theory, the surface of any object, such as leaves, branches, and tree crowns, is composed of many smaller mirrors after reaching the microscopic scale. The main feature of BRDF is that the microplane normal obeys a statistical distribution. Therefore, the reflected radiance differential dI in formula (7) is ref (ω r ) can be expressed as:

[0118]

[0119] Among them, θ ξ is the incident light direction vector ω i and the reflection microsurface normal vector ω h The angle between r is the triangle patch normal vector ω n and the reflected light direction vector ω r The angle between n is the normal vector calculated for the triangle patch, ω i is the incident light direction vector, ω r is the direction vector of the reflected light, ω h is the normal vector of the reflected microplane, and dA is the area of ​​the current microplane. ξ =ω i ·ω h , cosθ r =ω n ·ω r In formula (8), according to the calculation model of surface roughness, the normal distribution function D(ω h ) is defined as:

[0120]

[0121] Where the reflection microfacet normal vector is defined as θ h is the reflection microplane normal vector ω h and the triangle patch normal vector ωn The angle between h =ω h ·ω n α is the roughness of the micro-surface. The larger α is, the smaller the reflection energy is. Conversely, the smaller α is, the larger the reflection energy is. + The function is a step function, that is, when the input is greater than 0, the output is 1, otherwise the output is 0.

[0122] In formula (8), F r is the Fresnel coefficient, which is expressed according to the Cook-Torrance model:

[0123]

[0124] The variables c and g are used to simplify the formula expression, c = |ω r ·ω h |, η t is the leaf refractive index. The results of the BRDF model usually depend on the setting of leaf surface parameters, the roughness α of the leaf surface and the refractive index η t According to previous literature, the ranges are usually 0.078-0.46 and 1.22-1.68, respectively. Laurel with a thicker cuticle has a smooth surface, and its α ranges from 0.078 to 0.12 and η t The range is 1.22 to 1.51; the hairs on the surface of hazelnut leaves give it a greater roughness (α is 0.46) and a higher refractive index (η t Other literature records include beech leaves (α ranges from 0.29 to 0.32, η t ranged from 1.42 to 1.56), liquidambar and maple leaves (α was about 0.5), wheat (α ranged from 0.18 to 0.46, η t is 1.27), sorghum (α range is 0.27~0.56, η t This study used the OHSP660 of Hangzhou Hongpu Optoelectronics Technology Co., Ltd. to measure the transmittance and reflectance of leaves, and combined the results of multiple measurements with the least squares fitting formula to calculate the average roughness α of each leaf: 0.32 for crape myrtle, 0.22 for cherry blossom, and 0.20 for camphor; and the average refractive index η t : Crape myrtle 1.82, cherry blossom 1.68, camphor 1.57. Mango and rubber trees are virtual trees reconstructed from public datasets, so relevant measurements cannot be made on their leaves.

[0125] Substituting formula (8) into formula (7) yields:

[0126]

[0127] The different solid angles ω at the center of the hemispherical coordinates can be expressed as Where θ is the angle with the zenith direction, is the azimuth. i As the zenith direction, it can be converted into spherical coordinate form ω h (1,θ ξ ,φ) and ω r (1,2θ ξ ,φ), then:

[0128]

[0129] Considering the influence of the geometric attenuation term G caused by the shading or shadow effect of reflected light due to the uneven micro-surface, formula (11) combined with formula (12) is transformed into the final bidirectional reflectance distribution function, which is recorded as:

[0130]

[0131] The expression equation of the geometric attenuation term G is:

[0132]

[0133] By determining the incident ray ω i The triangular patch t formed by the point cloud of leaves in the crown l,s The first intersection, when the light intersects with the triangular face, indicates that a beam of light is intercepted by the canopy leaves. l,s The reflected radiance of each incident ray emitted from Through Figure 1 (b) Random importance sampling is performed on the discrete points on the reflection lobe map surface, and M reflection directions are selected as Figure 1 As shown by the blue arrow in (b), the reflected radiance of the mth light ray is expressed as follows, where m = 1, 2, ..., M:

[0134]

[0135] Represents the direction vector of the mth reflected light; I inc (ω i )=I inc ; M is 10, the incident light and the triangle t l,s After the reflected light emitted by the intersection intersects with other triangular facets, the intercepted reflected radiation flux F ref It can be expressed as:

[0136]

[0137] Where S refRepresents the sum of the areas of the triangular patches that intercept the reflected light, θ ref is the reflected light ω r Along its direction and other normal vectors is ω n The intersecting triangles t l,s The angle between

[0138] 2.4. Calculation of the transmitted radiation flux of leaves in the crown:

[0139] The bidirectional transmission distribution function (BTDF) defines the relationship between the irradiance in the incident direction and the radiance in the transmission direction. The differential radiance of the light transmitted in any direction is proportional to the differential radiance of the incident light in any direction. The BTDF function can be expressed as:

[0140]

[0141] Among them, η i is the refractive index of the medium on the incident light side, η o is the refractive index of the medium on one side of the transmitted light. Since both sides of the blade are air, it is approximately assumed that η i =η o . Fresnel coefficient F r , the normal distribution function D, and the geometric attenuation term G are similar to the expressions in formula (13). t is the direction vector of the transmitted light, θ t Transmitted light from the sun t and the normal vector ω of the leaf triangle patch that intersects it n The angle between t =ω t ·ω n The transmission microfacet normal vector is defined as

[0142]

[0143]

[0144]

[0145] Triangle patch t l,s The transmitted radiance emitted from above is recorded as: Through Figure 1 In (c), random importance sampling is performed on the surface discrete points of the transmission lobe map on the left side, and M directions are selected such as Figure 1 As shown by the green arrow in (c), the transmitted radiance of the mth light is expressed as follows:

[0146]

[0147] Represents the direction vector of the mth transmitted light;

[0148] Incident light and triangle t l,s If the transmitted light emitted after the intersection is intercepted by other triangular facets, the intercepted transmitted radiation flux F tra It can be expressed as:

[0149]

[0150] Where S tra is the sum of the areas of the triangular patches that intercept the light transmitted from other leaves, θ tra is the transmitted light ω t With other triangles t l,s After intersection, ω t and the intersecting triangle face normal vector ω n The angle between In formula (22), the area of ​​the transmitted light intercepted by the triangle is denoted as a l,s .

[0151] According to the ray tracing algorithm, traverse all triangles t l,s ∈T l The solar radiation flux that is ultimately intercepted (absorbed) within the tree canopy can be calculated As shown in the following formula, it includes the radiation flux generated by the direct, reflected and transmitted light intercepted by the canopy leaves, and subtracts the reflected and transmitted radiation flux that escapes after the direct light collides once, S inc The sum of the areas of the triangles that intercept the incident light.

[0152]

[0153] 2.5. Field verification:

[0154] To verify the experimental accuracy of our method, we used a photoelectric global solar radiation transmitter (model RS-RA-N01-AL, Shandong Renke Measurement and Control Technology Co., Ltd., Jinan, Shandong Province) to measure the solar irradiance above and below the tree canopy. The handheld rapid measurement recorder contains a photoelectric global solar radiation transmitter, which is used to measure the solar radiation received from the zenith hemisphere in the spectral range of 300 to 3000nm. Before reading the measurement value, it is necessary to adjust the tripod and check the state of the level bubble on the sensor base to determine whether the sensor is parallel, and the measurement personnel squat to avoid blocking the upper side light to affect the measurement value. The sensor uses a high-precision photosensitive element with high absorption in the full spectrum range and good stability; the sensor uses the standard Modbus-RTU protocol and communicates with the handheld data recorder through the RS-485 serial interface to read the current solar radiation value. The measurement range of the sensor is 0 to 1800W·m-2 .

[0155] An open area (with no shade above) was selected near the experimental trees to replace the solar irradiance E above the canopy. above , and a handheld instrument was used to measure the solar irradiance E below the canopy of the experimental trees. below Since the values ​​at different locations under the canopy are slightly different, it is necessary to uniformly sample and measure at multiple points under the canopy and calculate the average value. At this time, the actual measured solar radiation flux intercepted by the canopy is recorded as

[0156]

[0157] Among them, Area procanopy It is the projection area of ​​the tree crown along the incident direction of sunlight, which can be obtained by projection based on the 3D convex hull algorithm of the tree crown point cloud combined with the solar altitude angle and azimuth at that moment.

[0158] 3. Results and discussion

[0159] 3.1. Calculation results of solar irradiance in the experimental area:

[0160] Using the local solar radiation calculation method described in Section 2.1 of this article, combined with the longitude and latitude of the experimental site, solar altitude angle, solar azimuth angle, and corresponding air quality coefficient parameters, the 15th day of April, June, August, and October 2020 was selected to calculate the hourly solar shortwave radiation radiance of the surface from 7:00 to 17:00 under clear sky conditions. The results are as follows: Figure 2 As shown. At noon, the solar altitude angle is close to the directly overhead position, and the solar shortwave irradiance is large. At other times, the solar altitude angle is low, so the sun needs to travel a longer distance in the atmosphere, and the angle between the sun and the normal line of the horizon increases. At this time, the solar shortwave radiation irradiance reaching the surface is small. Figure 2 The solar shortwave radiation irradiance at each moment from 7:00 to 17:00 on the 15th day of each month in April, June, August, and October 2020 under clear sky conditions calculated by the algorithm in this paper is at the experimental plot.

[0161] 3.2. BRDF model for leaves:

[0162] To show the BRDF and BTDF calculation results of the blade, in two-dimensional polar coordinates, the azimuth angle of 0 to 360° is represented by the angular coordinate, and the zenith angle of 0 to 90° is represented by the radial coordinate (e.g. Figure 3As shown in the figure, the first line shows the calculation results of BRDF under the upper hemisphere distribution, and the second line shows the calculation results of BTDF under the lower hemisphere distribution. Now assume that the azimuth angle of the incident light is 180°, and the zenith angle increases from 0° to 80° in increments of 20°. In the figure, the incident light is represented by a blue ×, and the colors in the two-dimensional polar coordinate system represent the BRDF and BTDF values ​​calculated at all viewpoints in the entire hemisphere under the incident light.

[0163] Figure 3 The polar coordinate diagram shows the calculation results of the bidirectional reflection distribution function (BRDF) (first row) and bidirectional transmission distribution function (BTDF) (second row) corresponding to the incident light at different angles. The parameters are set to roughness α equal to 0.2 and refractive index η t is equal to 1.5, the azimuth angle of the incident light (indicated by the blue ×) is 180°, and the zenith angle starts from 0° from left to right and increases in increments of 20° to 80°. Figure 3 The results presented show that when the zenith angle is 0°, the BRDF value f r The fluctuations in all directions are not much different. As the zenith angle increases, the BRDF peak value increases, and the value on the side of the incident light decreases. The above experiment shows that when we calculate the angle between the light and the leaf normal vector increases, f r The value increases accordingly, and the reflectivity is changed from f in formula (15) r (ω r ,ω i )cosθ i Therefore, the reflectivity still decreases with the increase of the angle between the light and the leaf surface normal vector, that is, when the angle between the light and the leaf surface normal vector is small, more energy is reflected; and the value of BTDF f t When the angle between the sunlight and the leaf normal vector is small, the transmittance is higher. t (ω t ,ω i )cosθ i Decision, such as Figure 4 As shown in the second row, as the zenith angle of the incident light increases, the value f t Smaller.

[0164] 3.3 Visualization results of the interception of solar shortwave radiation by leaves in the canopy:

[0165] This paper reconstructs the trees in three dimensions based on the calculated tree parameters, and uses ray tracing methods and BRDF and BTDF to simulate the calculation of sunlight in the canopy. Figure 4 and Figure 5It shows the computer simulation of the direct sunlight shortwave radiation received by each experimental tree at 11:00 on August 15, 2020 (solar altitude angle 67.23°, solar azimuth angle 139.52°), as well as the path results of the reflected and transmitted light emitted after the direct light collides with the leaves, and the distribution of irradiance in the canopy at that moment.

[0166] Figure 4 Taking the mango tree (left) and the cherry tree (right) as examples, the computer simulation method in this paper shows the direct sunlight in the canopy at 11:00 on August 15, 2020, and the random reflection and transmission light paths generated based on Monte Carlo. Part of the direct light (red) is blocked by the triangular facets of the leaves in the canopy, generating multiple random directions of reflection (blue) and transmission light (green) to continue to intersect with other leaves. Figure 4 It can be seen that because mango trees have a higher canopy leaf area index than cherry trees (3.48 for mango trees and 1.33 for cherry trees), the leaves in the canopy of mango trees intercept more incident direct light (65.1% for mango trees and 47.2% for cherry trees), and according to the direction of the BRDF, 58.1% of the reflected light collides with other leaves in the same canopy, and the rest (41.9%) escape to the sky and the surrounding environment, while only 32.6% of the reflected light in cherry trees collides with other leaves in the same canopy. Similarly, because the mango tree has a denser leaf area density, 66.2% of the transmitted light intersects with the middle and lower leaves of the crown, which enables the shade leaves growing in a weak light environment to produce a certain amount of photosynthesis. Similarly, the cherry tree also has a certain amount of transmitted light (47.9%) intersecting with the middle and lower leaves of the crown, but the proportion is significantly smaller than the mango crown with a higher leaf area density. This is consistent with the actual characteristics of mango trees that like the sun and fear the shade, and require a ventilated and light-transmitting environment to improve fruit quality.

[0167] Figure 5 The method used in this paper shows the visualization distribution of direct, reflected and transmitted irradiance in the canopies of four experimental trees at 11:00 on August 15, 2020; among them, mango trees (first row), rubber trees (second row), crape myrtle trees (third row), and cherry trees (fourth row), the first column is the direct solar irradiance of the four tree crowns, the second column is the reflected solar irradiance of the four tree crowns, and the third column is the transmitted solar irradiance of the four tree crowns. Figure 5It reflects the irradiance of different leaf areas in each tree crown at 11:00 on August 15, 2020. The total solar radiation flux intercepted by the canopy is calculated by the irradiance and the corresponding triangular patch area. At the same time, the proportion of each radiation component in the four experimental trees to the total radiation flux intercepted by the crown can be obtained. The results show that the ratio of direct radiation flux to the total radiation flux intercepted is the largest, among which the crape myrtle tree accounts for 91.87% (323.34W), the cherry tree accounts for 90.20% (1161.66W), the rubber tree accounts for 87.76% (3450.01W), and the mango tree accounts for 83 .24% (834.32W); the second is the transmitted radiation flux, specifically: mango trees account for 10.77% (107.95W), rubber trees account for 8.80% (345.94W), cherry trees account for 7.12% (91.70W), and crape myrtle trees account for 6.20% (21.83W); the lowest proportion is the reflected radiation flux, specifically: mango trees account for 5.99% (60.04W), rubber trees account for 3.44% (135.23W), cherry trees account for 2.68% (34.52W), and crape myrtle trees account for 1.92% (6.77W).

[0168] Corresponding to the crowns of the four trees, since the leaf area index of the cherry tree is smaller than that of the mango tree, only 32.6% of the reflected light emitted collides with other leaves in the same crown, and the solar altitude angle at 11 o'clock is larger, the reflected light is more likely to be reflected above the crown, and the contribution of reflected radiation flux is relatively small (2.65%). The leaves of the mango tree are densely clustered, and a single leaf has a larger area, which can easily intercept the reflected light generated by direct radiation. The mango crown is tower-shaped as a whole, and the lower branches and leaves extend outward, which makes the probability of incident light being intercepted in the lower part of the crown higher than that of other experimental trees, so it can reflect more light to the middle and upper leaves of the crown, generating more reflected radiation flux contribution (5.99%). In addition, a larger total leaf area and a smaller crown projection area mean a higher degree of leaf overlap, resulting in a greater probability of secondary collisions of the sun's transmitted light in subsequent paths, which makes the mango tree with a larger crown leaf area index contribute a relatively large amount of transmitted radiation flux (10.77%), while the crape myrtle tree with the lowest leaf area index has a lower probability of re-interception of transmitted light due to its sparse leaf distribution, resulting in the smallest contribution to transmitted radiation flux (6.46%). For the rubber tree with a larger crown and a moderate leaf area index, its direct, reflected and transmitted proportions are all in the middle of the four trees.

[0169] The ray tracing method in this article is also applicable to data of multiple trees. Figure 6This is a calculation diagram and a local enlarged diagram of the solar irradiance distribution in the canopy of several camphor trees on campus at 11:00 on August 15, 2020. The canopy is divided into three areas: Ⅰ (12.0m-16.1m), Ⅱ (7.9m-12.0m), and Ⅲ (3.8m-7.9m) from the vertical direction according to the height from the crown top to the bottom of the branches. It is common sense that the intercepted solar shortwave radiation is mainly distributed in areas Ⅰ and Ⅱ in the forest canopy. Due to the top-down scanning method of the ground-based lidar, the leaf area of ​​the forest canopy region I calculated by the method in this paper is 99.25m2, but the solar shortwave radiation flux intercepted by the leaves in this layer is 28485.99W, accounting for 36.72% of the total; Region II has a relatively large leaf area of ​​278.47m2, and its solar radiation mainly comes from the oblique sunlight, and the radiation flux generated by the direct light and reflected and transmitted light passing through the pores of the leaf clusters. The solar shortwave radiation flux intercepted by the leaves in this region is 41262.80W, accounting for 53.19% of the total; Region III is blocked by the upper regions I and II, and the solar radiation flux intercepted is significantly smaller than that of regions I and II. The leaf area is 56.22m2, and the solar shortwave radiation flux intercepted is only 7827.44W, accounting for only 10.09% of the total, indicating that the shrub layer under the forest is generally shade-tolerant plants, which mainly rely on transmitted light to absorb energy from nature.

[0170] Figure 7 The radiation flux values ​​of the sun intercepted by four experimental trees (a, mango tree, b, rubber tree, c, crape myrtle tree, d, cherry tree) and e, a camphor tree forest calculated by the method in this paper as time goes by; the first, second and third columns represent the direct, reflected and transmitted radiation flux values ​​intercepted by the tree canopy, respectively. Figure 7The distribution of radiation flux of the experimental trees exposed to sunlight is shown. The solar radiation flux intercepted by the experimental tree canopy from 7:00 to 17:00 is calculated on the 15th day of each month in four months of the year (April, June, August and October) using the algorithm in this paper. The three parts of direct radiation flux, reflected radiation flux and transmitted radiation flux are included. The solar shortwave radiation energy intercepted by the canopy generally presents a single peak change trend over time within a single day. At noon, the solar irradiance reaches its maximum, and the phototropism of leaves makes the ratio of direct solar light intercepted at noon relatively high. Therefore, the solar radiation flux intercepted by the crowns of mango trees, rubber trees and cherry trees reached a peak at 12:00 on June 15, which were 1347.18W, 5437.99W and 1704.34W respectively; the crowns of crape myrtle and several camphor trees reached a peak at 13:00, which were 488.16W and 100079.93W respectively. At the same time, we also found that when the solar altitude angle is small, such as in the morning and afternoon, the reflected light is more easily intercepted by the crown. This is because the light obliquely entering the crown is reflected by the leaves with a smaller zenith angle and is more easily absorbed by the leaves on the other side of the crown. However, since the solar irradiance decreases at a low solar altitude angle, the reflected radiation flux of the crown reaches a minimum in the morning and evening. The transmitted radiation flux intercepted by the crown is related to the propagation distance of the light within the crown and the overlap of the leaves. This allows the crown with a higher leaf area index and a larger crown volume to capture more transmitted light energy. However, due to the changes in the solar altitude and azimuth angles at each moment, as well as the heterogeneity of the spatial structure and vertical distribution of the crown's leaves, the calculation results at each moment also show some irregularity.

[0171] In order to verify the effectiveness of the method in this paper, we used a photoelectric solar radiation transmitter to conduct field verification at 9:00, 12:00, and 15:00 on June 7, 2020, to quantitatively measure the solar irradiance values ​​of the upper and lower canopies, and then combined with the canopy projection area to obtain the radiation flux values ​​at the corresponding time. The specific method is described in Section 2.5. Table 2 shows the canopy radiation flux interception values ​​calculated using the ray tracing method in this paper and the corresponding values ​​measured by field instruments. The data show that the ray tracing algorithm in this paper has a high accuracy rate. The interception efficiency of the same canopy for light is greater at noon than in the morning and afternoon. This is because the phototropism of leaves and the distribution of leaf inclination angles make the canopy photosynthesis most active and the light interception efficiency reaches the maximum value when the solar irradiance is highest at noon.

[0172] Table 2 is a verification comparison of the radiation flux above and below the canopy of each experimental forest at 9:00, 12:00 and 15:00 on June 7, 2020 calculated by the method in this paper and the field measurement values. As shown in Table 2, at 9:00 a.m. and 15:00 p.m., the solar shortwave radiation interception efficiency was 24.53% and 27.06% for the crape myrtle, 30.91% and 31.54% for the cherry tree, and 40.86% and 43.21% for the camphor trees, respectively, and the deviation from the actual measurement value was within 5%. At 12:00 noon, the solar shortwave radiation interception efficiency increased to 32.12% for the crape myrtle, 40.92% for the cherry tree, and 55.36% for the camphor trees. At 12:00 noon, the projection area was close to the vertical projection area, and the interception efficiency showed a great correlation with the leaf area index; due to the deviation in the growth environment of the forest, there will be spatial competition and uneven distribution of forest windows between the canopies, and the solar altitude angle and azimuth will change at each moment, but Table 2 shows that the results of this paper are in good consistency with the field verification.

[0173] Table 2:

[0174]

[0175] The program in this article runs on the Microsoft Windows 10 operating system, and Matlab 2020a is used for code writing and execution. The computer is equipped with an Intel(R) Core(TM) i7-10750H 6-core CPU 2.6GHz, 32GB memory and NVIDIAGeForce GTX1650 graphics processor. Because it involves a large number of light tracing and intersections within the canopy, the running time of the entire algorithm is determined by the number of emitted light rays and the number of triangles contained in the leaves and branches in the canopy. Here we use an axisymmetric bounding box strategy to accelerate the calculation of light intersection, that is, the three-dimensional model of the forest is divided into multiple continuous bounding boxes, and the intersection of the light and each fine triangle in the bounding box is calculated. First, it is necessary to determine the intersection of the light with the larger bounding box to reduce the amount of calculation. The program-related parameters are set to a light sampling spacing of 0.2cm and the number of emitted sunlight rays is 2.5*10 5 per square meter; the point cloud interval before Delaunay triangulation is 0.4cm, so as to ensure that the light can cover each triangular face as much as possible. The program execution time of different trees in the experiment is shown in Table 3, which shows the number of triangular faces contained in the leaves and branches in the experimental crown, the number of sunlight emitted above the scene and the algorithm execution time. At the same time, in order to reflect the time consumption reduced by the axis-aligned bounding box, the acceleration efficiency is now defined as the ratio of the execution time reduced by acceleration to the original execution time. The average acceleration efficiency in each experimental data is 65.6%.

[0176] Table 3:

[0177]

[0178] 4. Conclusion:

[0179] The interception of solar shortwave radiation by trees is closely related to the growth and development of trees. How to accurately depict the propagation process of solar shortwave radiation in the canopy and the distribution of radiation flux in the canopy is of great significance. From the perspective of interdisciplinary studies, this paper applies surveying and mapping remote sensing technology and computer graphics to practical problems in forestry, finely depicts the spatial branch and leaf morphological structure of trees, and accurately reconstructs the real three-dimensional model of trees. Secondly, through the ray tracing algorithm and the bidirectional reflection and transmission distribution function, the propagation process of direct, reflected and transmitted sunlight in the canopy is simulated, and then the radiation flux under different canopy temporal and spatial changes is calculated. Finally, by comparing with the field instrument measurement values, the algorithm in this paper has achieved an accuracy of more than 90%. The experiment found that the larger the leaf area index, the greater the overlap of leaves and the sum of leaf area per unit area, which will increase the direct radiation flux and transmitted radiation flux obtained by the canopy per unit area. When the solar altitude is low in the morning and afternoon, the direct and transmitted radiation flux received by the tree crown will be reduced, but the oblique light is more likely to collide with other leaves in the tree crown after reflection, making the reduction of reflected radiation flux in the tree crown relatively unobvious. At noon, the phototropism of the leaves in the tree crown causes allometric growth and leaf inclination angle distribution to intercept more sunlight, but most of the reflected light escapes into the sky.

[0180] In future research, the algorithm in this paper can be improved in the following aspects: 1. This paper uses a ground laser scanner to collect point cloud data. The bottom-up perspective occlusion of tall trees will result in missing data above the tree canopy, which has a certain impact on the use of pattern recognition algorithms for single leaf segmentation and reconstruction and ray tracing of the upper tree canopy. It is possible to consider combining ground-based and airborne point clouds to obtain a complete canopy point cloud. 2. Since many trees have been finely modeled, all triangular facets in the canopy need to be traversed during ray tracing, and the intersection of a large number of rays is time-consuming. Therefore, how to use appropriate parallel computing, GPU computing, or Monte Carlo random large number laws to accelerate the intersection of rays and triangular facets and reduce the complexity of the algorithm is a further optimization idea.

[0181] The protection scope of the present invention includes but is not limited to the above embodiments. The protection scope of the present invention shall be based on the claims. Any replacement, deformation, and improvement of the technology that can be easily thought of by technicians in this field shall fall within the protection scope of the present invention.

Claims

1. A quantitative analysis method of solar shortwave radiation in tree canopy based on a computer simulation model, characterized in that: The following steps are involved: Step 1: Collect forest stand point cloud data; Step 2: Classify all tree point cloud data in the forest stand into leaf point cloud data and branch point cloud data; perform single leaf segmentation and triangulation on the leaf point cloud data, and then represent each leaf as a set of triangular facets; Step 3: Calculate the surface solar irradiance at the forest stand location; Step 4: Model the sunlight as parallel rays emitted from above the trees, and obtain the area of ​​a current triangle directly illuminated by the sunlight by intersecting the straight line with the triangle. Step 5: Calculate the solar radiation flux of each incident direction vector through the number of incident parallel rays and the surface solar irradiance, calculate the incident radiance of the incident solar ray through the solar radiation flux, the area of ​​the current triangular facet obtained, and the angle between the incident solar ray direction vector and the normal vector of the triangular facet directly illuminated by the incident solar ray, and calculate the direct radiation flux intercepted by multiple triangular facets; Step 6: Determine the BRDF value through the direction vector of the incident sunlight, the normal vector of the triangle directly hit by the incident sunlight, the direction vector of the reflected light, and the normal vector of the reflected microplane, and then obtain the reflected radiance of the incident sunlight. After the second intersection, obtain the other triangles that intercept the reflected light, and calculate the reflected radiation flux intercepted by other triangles; Step 7: Determine the BTDF value through the direction vector of the incident sunlight, the normal vector of the triangle directly hit by the incident sunlight, the direction vector of the transmitted light, and the normal vector of the transmitted microplane, and then obtain the transmitted radiance of the incident sunlight. After the second intersection, obtain the other triangles that intercept the transmitted light, and calculate the transmitted radiation flux intercepted by other triangles. Step 8: Traverse all triangular facets and calculate the solar radiation flux that is ultimately intercepted by all triangular facets within the tree canopy.

2. The method for quantitative analysis of solar shortwave radiation from tree canopy based on computer simulation model according to claim 1, characterized in that: In step 2, the leaf point cloud data is triangulated, specifically: the point cloud data of each leaf is triangulated using Delaunay triangulation, and each leaf is represented as a set of triangular facets t l,s ∈T l , where t l,s It is the sth triangle in the lth leaf.

3. The method for quantitative analysis of solar shortwave radiation from tree canopy based on computer simulation model according to claim 1, characterized in that: The step 3 is specifically as follows: Surface solar irradiance E at the forest stand location total for: AND total =And direct +E diffuse (1); The direct solar irradiance is E direct =E sc d m τ r τ o τ g τ w τ a sinθ s ; E sc is the solar constant; d m is the eccentricity coefficient of the Earth's orbit; τ r , τ o , τ g , τ w and τ a represent the scattering transmittances of Rayleigh, ozone, gas, water, and aerosol, respectively; θ s is the local solar altitude angle, expressed as follows: in, is the local latitude, δ s is the local solar declination angle, τ s is the local solar hour angle, is the incident sun ray vector; Horizontal diffuse irradiance E diffuse for: AND diffuse =And dr +E da +E dm (3); Among them, E dr is the Rayleigh scattering of sunlight after entering the atmosphere, E da is the aerosol scattering of sunlight after it enters the atmosphere, E dm It is multiple reflections from the ground and the sky.

4. The method for quantitative analysis of solar shortwave radiation from tree canopy based on computer simulation model according to claim 2, characterized in that: The step 4 is specifically as follows: Model the sunlight as parallel rays emitted from above the trees, and determine the sth triangle t in the lth leaf of the canopy point cloud construction. l,s The first intersection, triangle patch t l,s The area is a l,s .

5. The method for quantitative analysis of solar shortwave radiation from tree canopy based on computer simulation model according to claim 4, characterized in that: The step 5 is specifically as follows: Calculate the number of solar rays N emitted per square meter based on the fixed interval between adjacent parallel rays inc , then the solar radiation flux of each incident direction vector is: The incident radiance of the incident sunlight is: When the incident ray intersects the triangle, it means that there are n inc The incident light beam is reflected by the current triangle patch t of the canopy leaf l,s Intercept; where θ i is the direction vector of the incident sunlight i and the triangle t that intersects the incident sun ray l,s Normal vector ω n The angle between Direct radiation flux intercepted by canopy leaves Where S inc It is the sum of the areas of the triangular faces that intercept the incident sunlight.

6. The method for quantitative analysis of solar shortwave radiation from tree canopy based on computer simulation model according to claim 5, characterized in that: The step 6 is specifically as follows: The BRDF values ​​are: in: where θ r is the triangle patch normal vector ω n and the reflected light direction vector ω r The angle between n is the triangle facet that intersects with the incident sun ray l,s Normal vector; ω i is the direction vector of the incident sunlight; θ h is the reflection microplane normal vector ω h and the triangle patch normal vector ω n The angle between h =ω h ·ω n ; α is the roughness of the microsurface; χ + The function is a step function; c = |ω r ·ω h |, η t is the blade refractive index; the reflection microfacet normal vector is defined as By determining the incident sun ray direction vector ω i The triangular patch t formed by the point cloud of leaves in the crown l,s The first intersection, when the light intersects with the triangular face, indicates that a beam of light is intercepted by the canopy leaves; from the first intersection of the triangular face t l,s The discrete form of the reflected radiance of each incident ray emitted on is: By performing random importance sampling on the discrete points on the surface of the reflection lobe pattern and selecting M reflection directions, the reflection radiance of the mth light is expressed as follows, where Among them, I inc (ω i )=I inc ; Represents the direction vector of the mth reflected light; Incident light and triangle t l,s After the reflected light emitted by the intersection intersects with other triangular facets, the intercepted reflected radiation flux F ref It can be expressed as: Where S ref is the sum of the areas of the triangular patches that intercept the light reflected from other leaves, θ ref is the direction vector of the reflected light r Along its direction and other normal vectors is ω n The intersecting triangles t l,s The angle between 7. The method for quantitative analysis of solar shortwave radiation from tree canopy based on computer simulation model according to claim 6, characterized in that: The step 7 is specifically as follows: The BTDF values ​​are: Where: η i is the refractive index of the medium on the incident light side, η o is the refractive index of the medium on the side of the transmitted light, ω t is the direction vector of the transmitted light, θ t is the transmitted light direction vector ω t and the triangle normal vector ω of the leaf surface that intersects with it n The angle between t =ω t ·ω n ; The transmission microfacet normal vector is defined as Triangle patch t l,s The transmitted radiance emitted from above is recorded as: By performing random importance sampling on discrete points on the surface of the transmission lobe diagram and selecting M directions, the transmission radiance of the mth light is expressed as follows: in, Represents the direction vector of the mth transmitted light; Incident light and triangle t l,s If the transmitted light emitted after the intersection is intercepted by other triangular facets, the intercepted transmitted radiation flux F tra It can be expressed as: Where S tra Represents the sum of the areas of the triangular patches that intercept the transmitted light, θ tra is the transmitted light direction vector ω t With other triangles t l,s After intersection, ω t and the intersecting triangle face normal vector ω n The angle between In formula (17), the area of ​​the triangular facet that intercepts the transmitted light is recorded as 8. The method for quantitative analysis of solar shortwave radiation from tree canopies based on a computer simulation model according to claim 7, characterized in that: The step 8 is specifically as follows: According to the ray tracing algorithm, traverse all triangles t l,s ∈T l , the solar radiation flux intercepted by all triangular facets in the canopy is calculated as:

Images