Method for measuring vegetation canopy structure based on non-fisheye digital camera

Abstract

The invention relates to a method for measuring a vegetation canopy structure based on a non-fisheye digital camera. The method comprises the following steps: (1) calculating the actual length and width corresponding to a view field and an image of the non-fisheye digital camera; (2) shooting and sampling to obtain an image; (3) converting an RGB color space of the image into an LAB color space; (4) preprocessing the wave band A in the LAB color space to obtain a binary image of the background and the vegetation; (5) calculating FVC, CC, porosity, FAPAR of direct solar radiation and diffused sky radiation, FAPAR of average direct solar radiation in one day and the proportion of the diffused radiation in the total radiation, and further calculating to obtain instantaneous FAPAR; (6) calculating cumulative pore size distribution, vegetation width of ridge row crops and distance between ridges in each pixel direction; (7) calculating CI, LAI and LAIe with slight errors, and accurate LAI and LAIe; (8) calculating an average leaf inclination angle. The method is low in labor cost, free of image deformation and wide in applicability.

Description

Technical Field

[0001] The invention relates to the field of vegetation canopy structure research and vegetation quantitative remote sensing technology, and in particular to a method for measuring vegetation canopy structure based on a non-fisheye digital camera. Background Art

[0002] Photo-based vegetation structure measurement is a camera optical measurement method that originates from the mathematical study of vegetation canopy structure. In the study of vegetation canopy structure, Soviet academician J. Ross gave a detailed description of the relationship between vegetation canopy structure and radiation in his book "The Radiation Regime and Architecture of Plant Stands" [Ross J. The Radiation Regime and Architecture of Plant Stands [M]. Berlin Heidelberg: Springer, 1981]. In order to describe the relationship between a bunch of leaves and radiation in a canopy in a single direction, a leaf normal direction distribution function based on statistical principles was proposed, which can describe the distribution state of leaves in the canopy with statistical probability. However, this function is difficult to measure in the actual vegetation canopy. For this reason, J. Ross proposed a famous G function (i.e., the projection of the leaf area in the direction of the leaf normal within the unit solid angle). The projection function was proposed to quantify the interception ratio of radiation (or light), and a mathematical and physical relationship between radiation and leaf area was established [Ross J, Nilson T. The extinction of direction radiation in crops. In: Question on Radiation Regime of Plant Stand [R]. Tartu (Russian): Acad. Sci. ESSR, Inst. Phys. Astron, 1965: 25-64.]. In the same period of research on extinction, J. Ross and T. Nilson proposed an exponential optical penetration function based on Beer's law in 1965. This function can well describe the attenuation process of radiation with the optical thickness of vegetation [Ross J, Nilson T. The extinction of direction radiation in crops. In: Question on Radiation Regime of Plant Stand [R]. Tartu (Russian): Acad. Sci. ESSR, Inst. Phys. Astron, 1965: 25-64.]. J. Ross's early research laid the foundation for the study of the interaction between vegetation canopy structure and radiation, and provided a theoretical basis for measuring vegetation canopy structure with optical instruments.

[0003] Subsequent scholars ABG Lang and JM Chen conducted research on common vegetation such as crops and forests [Lang AR G. Simplified estimate of leaf area index from transmittance of the sun's beam. [J]. Agric. For. Meteorol, 1987, 41: 179-186.] [Chen JM, Cihlar J. Plant canopy gap-size analysis theory for improving optical measurements of leaf-area index [J]. Applied Optics, 1995, 34 (27): 6211-6222.], and they respectively provided corresponding vegetation canopy structure theories to calculate leaf area index (LAI). The key issue in LAI calculation is the equation of the aggregation index. In addition to the clustering index proposed by T. Nilson for the exponential penetration function, the existing clustering index methods used in the field of vegetation radiation measurement include Liang and Xiao's clustering index, Chen and Cihlar's clustering index, Pielou Coefficient of Segregation clustering index and Lacunarity-Based clustering index [Fournier RA, Hall RJ. Hemispherical Photography in Forest Science: Theory, Methods, Applications [M]. Van Godewijckstraat: Springer, 2017.]. Since the photos taken by the camera are located in the Cartesian coordinate system, it has no angular coordinate axis (i.e., azimuth and zenith angle). Therefore, among the above four clustering index calculation equations, the Lacunarity-Based clustering index is the most suitable for calculation in photos. However, the Lacunarity-Based clustering index is based on the fractal principle. Although the Lacunarity-Based clustering index has a good advantage in mathematics, it is difficult to control the size of the search box in the fractal, which leads to inaccurate calculation in many cases.This brings great difficulties to the measurement of vegetation structure based on photos [Dale MR T. Lacunarity analysis of spatial pattern: A comparison [J]. Landscape Ecology, 15 (5): 467-478.] [Plotnick RE, Gardner RH, Hargrove WW, et al. Lacunarity analysis: A general technique for the analysis of spatial patterns [J]. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics, 1996, 53 (5): 5461-5468.]. At the same time, based on the research on vegetation canopy structure by J. Ross, ABG Lang, JM Chen, Sylvain G. Leblance, etc., other non-Cartesian coordinate (i.e. polar coordinate) imaging optical shooting technology has developed corresponding equipment. Taking LAI's measurement equipment as an example, this measurement method mainly includes imaging hemispherical photography technology, such as LAI-2000, canopy hemispherical photography camera, etc. In addition to imaging measurement technology, the most commonly used is the non-imaging photon detection pore probability technology, which uses PAR detectors (PAR detectors are mainly composed of photosensitive diodes + filters) to detect the probability of solar patches (pores) in the canopy, such as TRAC, SunScan, etc. Similarly, there are corresponding measurement equipment for vegetation coverage (FVC), canopy openness (CC), average leaf inclination angle (ALA) and absorbed photosynthetic radiation ratio (FAPAR). However, these devices have a common disadvantage: the measurement is a single-point measurement. For example, the fisheye lens has a viewing angle of 180°, a small field of view, and a limited observation range; while the photon detector method depends on the length of the detection rod. This limitation inevitably leads to the inability of a single measurement to reflect the average state of the canopy structure of the sample plot. In order to reflect the average state of the study area, the only way is to increase the number of measurements, which will inevitably lead to a lot of waste of manpower. At the same time, of the two measurement methods mentioned above, the first measurement method, that is, hemispherical photography technology, also has a limitation: the photo deformation is very large, which is particularly obvious when the zenith angle is 0-10° and the zenith angle is greater than 50°. Although the mapping function in the imaging can correct this image deformation, even if the image is processed after shooting, the pixel distortion and mixed pixel problems after interpolation are still difficult to solve, which is much more complicated than non-deformed conventional photos. As for the second measurement method, photon detection, it mainly measures the porosity probability. This method is mainly affected by the actual influence of previous artificial experiments. Therefore, different measurers in the same area may have large differences.

[0004] Therefore, there is an urgent need for a non-degenerative imaging measurement method that is less affected by human factors to solve the above problems. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a method for measuring vegetation canopy structure based on a non-fisheye digital camera with low labor cost, no image deformation and wide applicability.

[0006] In order to solve the above problems, the present invention provides a method for measuring vegetation canopy structure based on a non-fisheye digital camera, comprising the following steps:

[0007] ⑴ Determine the CCD or COMS frame size of the non-fisheye digital camera and the vertical distance from the shooting position of the camera lens to the top of the canopy, calculate the field of view of the non-fisheye digital camera according to the corresponding mathematical equation, and then calculate the actual length and width of the vegetation image taken;

[0008] (2) Selecting the shooting angle and shooting time of the image according to the measured vegetation parameters to perform shooting sampling to obtain the image; the shooting time refers to the time corresponding to a certain solar altitude angle automatically calculated according to the solar altitude angle corresponding time calculation program;

[0009] ⑶ Performing Gamma transformation on the RGB color space of the image, and then converting the transformed RGB color space image into LAB color space;

[0010] (4) Preprocessing the A band in the LAB color space to obtain binary images of the background and vegetation;

[0011] (5) Counting the ratio of the background and vegetation values ​​in the binary image, and calculating the vegetation coverage (FVC), canopy openness (CC), porosity, the ratio of absorbed photosynthetic radiation of direct solar radiation (back sky FAPAR) and the ratio of absorbed photosynthetic radiation of diffuse sky radiation (white sky FAPAR); at the same time, according to the shooting time corresponding to the ratio of absorbed photosynthetic radiation of direct solar radiation, the average ratio of absorbed photosynthetic radiation of direct solar radiation in a day is calculated, and the ratio of diffuse radiation to total radiation is calculated according to the equation of the ratio of diffuse radiation to total radiation, and then the instantaneous FAPAR is calculated;

[0012] (6) According to the binary image, the leaf chord length and pore length in the vertical direction of the image, the leaf chord length and pore length in the horizontal direction, and the leaf chord length and pore length in a specific direction are calculated by using the pore and leaf chord length pixel search method; and then the cumulative pore size distribution in each pixel direction is calculated; at the same time, the vegetation width and inter-ridge distance of the ridge crops are calculated according to the pore and leaf chord length pixel search method;

[0013] (7) According to the cumulative pore size distribution in each pixel direction, the aggregation index of the vegetation is calculated using the MX aggregation index equation; and then the real leaf area index with a slight error and the effective leaf area index with a slight error in the image taken in the vertical direction, as well as the accurate real leaf area index and the effective leaf area index in the image taken at an observation zenith angle of 57.5° are calculated based on the aggregation index;

[0014] ⑻ Based on the blade chord length in the vertical direction, the blade chord length in the horizontal direction, and the blade chord length in a specific orientation, the average leaf inclination angle of the vegetation in the image can be calculated using the average leaf inclination angle equation.

[0015] In the step (1), the non-fisheye digital camera is mounted on a telescopic rod or a drone.

[0016] The conditions for shooting and sampling in the step (2) refer to: when vertical shooting is adopted, the measured vegetation parameters are the cumulative pore size distribution, the vegetation width of the ridge crops, the ridge distance, the vegetation coverage rate, the canopy openness, the true leaf area index with a slight error, the effective leaf area index with a slight error, and the average leaf inclination angle; when the observation zenith angle is 57.5°, the measured vegetation parameters are the accurate true leaf area index and the accurate effective leaf area index; when measuring the ratio of absorbed photosynthetic radiation of direct solar radiation, shooting is performed at the zenith angle of the sun's incidence direction; when measuring the ratio of absorbed photosynthetic radiation of diffuse radiation from the sky, it is necessary to perform average sampling within 180 degrees of the zenith angle; when measuring the ratio of absorbed photosynthetic radiation of the average direct solar radiation in a day, it is necessary to use the time calculation program corresponding to the solar altitude angle to determine the shooting time.

[0017] The solar altitude angle corresponding time calculation procedure in step (2) refers to the combination of equations for calculating the solar altitude angle and the time corresponding to the solar altitude angle, and the equation is as follows:

[0018]

[0019] Where: θ ⊙_noon is the solar noon altitude angle; is the local latitude; δ is the solar declination; d n is the Julian Day; sunrise is the local sunrise hour angle; ω sunset is the local sunset hour angle; T sunrise is the local sunrise time; T sunset The local sunset time.

[0020] The pretreatment in the step (4) comprises the following steps:

[0021] ① Calculate the frequency of the A band value in the LAB color space, and construct a frequency histogram with the pixel value as the horizontal coordinate (x) and the number of pixels corresponding to the value (y) as the vertical coordinate;

[0022] ② Calculate the derivative f'[N(x)] of the frequency function f[N] composed of the number of pixels (y) of the A band value in the LAB color space: And calculate the extreme value in the frequency histogram so that the extreme value satisfies the requirement {f'[N(x)]>0}∧{f'[N(x)]>0};

[0023] Where: f'[N(x)] is the derivative of the double Gaussian distribution curve of soil and vegetation in band A; f[N(x)] is the first discrete value of the minimum value of the double Gaussian distribution curve of soil and vegetation in band A in the positive direction; f[N(x + )] is the second discrete value in the positive direction of the minimum value of the double Gaussian distribution curve of soil and vegetation in band A;

[0024] Then, the second local minimum value of the coordinate axis (x) of the pixel value in the positive direction is taken as the intersection point of the Gaussian distribution in the frequency of the vegetation value and the Gaussian distribution in the frequency of the background value, and the area near the intersection point is the area where the threshold value to be determined is located;

[0025] ③ Based on the intersection point and the points (x1, y1) and (x2, y2) near the intersection point, use the histogram slope search method equation to search for the minimum value of the Gaussian distribution in the frequency of the background value toward the green direction, and use this minimum value as the threshold;

[0026] The histogram slope search method equation is as follows:

[0027] D=y2-y1+tanβ(x1-x2), E=x1y2-x2y2,

[0028] F = x2-x1-tanβ(y1-y2),

[0029] Where: T is the threshold; D, E, and F are the functions in the threshold T equation respectively; point (x1, y1) and point (x2, y2) are points near the focus in the frequency histogram; β is the offset angle; τ is the control coefficient of tanα, and α is the angle between the straight line formed by the point (x1, y1) and the point (x2, y2) and the axis (y) formed by the number of pixels;

[0030] ④ Classifying the vegetation and background of the A band in the LAB color space according to the threshold value to obtain binary images of the background and vegetation; wherein the binary image 255 corresponds to the vegetation, and 0 corresponds to the background;

[0031] The diffuse radiation to total radiation ratio equation in step (5) refers to the ratio of the diffuse radiation of the sky to the total radiation, and the equation is calculated as follows:

[0032]

[0033] Where: sky is the ratio of the diffuse radiation incident on the top of the canopy to the total incident radiation; R p / R s It represents the ratio of photosynthetically active radiation to shortwave radiation.

[0034] The pore and leaf chord length pixel search method in step (6) refers to an algorithm for searching the pixel affiliation in an image, that is, whether the pixel belongs to vegetation or background; the algorithm expression is:

[0035] for i=1,L image

[0036] for j=1,W image

[0037]

[0038] Where: = is the mathematical symbol "defined as"; w is the blade chord length; λ is the pore size length; x and y are the pixel coordinate positions in the image coordinates; L image is the number of image lengths; W image is the number of image widths; P is the pixel value; the search conditions are:

[0039]

[0040] Taking the horizontal direction as the starting point, the image coordinate values ​​in the angle search are calculated using the following equation:

[0041] y i+1 =y i +1, where: η is the offset angle of the orientation in the image.

[0042] The MX aggregation index equation in step (7) is an indicator that describes the aggregation degree of leaves in the image, and its equation form is:

[0043]

[0044] Where: Ω E is the MX aggregation index; Λ(λ) is the heterogeneity index between all pixels in a single row (or column) of the image; Λ max (λ) is the maximum value of Λ(λ) in column (or row) statistics, and Λ max (λ) = L image , L imageis the number of image lengths; Λ cri (λ) is the value when the heterogeneity index between all pixels in a single row (or column) of the image is randomly distributed, and Λ cri (λ)=1;σ 2 (λ) is the statistical variance between all pixels in a single row of the image; E(λ) is the statistical expectation between all pixels in a single row of the image; i The pore size calculated for the pixels in a single row (or column) in the ith image; N is the statistical average value of the pore size in a single row calculated based on the pixels in a single row (or column) in the image; t The number of pores calculated for the pixels in a single row (or column) of the image; F(λ i ) is λ i The corresponding pore size cumulative distribution function value.

[0045] The average leaf inclination angle equation in step (8) refers to calculating the average value of the leaf inclination angle of the vegetation canopy, and the equation is:

[0046]

[0047] Where: θ l is the average leaf inclination angle; b is the bending factor of the blade; is a reference orientation; w is the average chord length of the blade corresponding to the reference orientation; l * w is the average chord length of the blade searched in the horizontal direction of the image coordinate (i.e., the direction of length); * is the average chord length of the blade searched in the vertical direction of the image coordinates (i.e., the direction of the width); γ is the sum of w and l * The angle between E is the MX aggregation index.

[0048] Compared with the prior art, the present invention has the following advantages:

[0049] 1. The present invention uses a histogram slope search method to search for image classification thresholds, and a pore and leaf chord length pixel search method to determine pixel affiliation, thereby solving the problem of automatic extraction of vegetation and background classification thresholds in images and preliminary measurement of vegetation structure in photo measurement technology. At the same time, the invention proposes a MX vegetation aggregation index, a planar leaf inclination angle calculation method, and a calculation of position parameters required for the ratio of absorbed photosynthetic radiation (FAPAR), thereby solving the difficulty in calculating vegetation structure parameters, thereby realizing a vegetation canopy structure measurement method based on a non-fisheye digital camera.

[0050] 2. Since the present invention is directed to a measurement method for measuring vegetation canopy structure using a non-fisheye digital camera, compared to the hemispherical photography technology of a fisheye lens, the present invention avoids the difficult problem of distorted image processing in spherical photos taken with a fisheye lens, making it more widely applicable.

[0051] 3. The present invention can load a camera onto a drone for high-altitude photography. Since the drone has a wide field of view (FOV), it is convenient to obtain the overall image of the study area, and then to analyze the average information of the vegetation canopy structure in the study area. Compared with the previous single-point measurement algorithm to obtain the average information of the canopy structure in the study area, this reduces manpower and measurement time.

[0052] 4. In the image preprocessing, the present invention proposes a histogram slope search method, which calculates the slope of the distribution curve in the green direction in the A band in the LAB color space through the background value frequency distribution, and then calculates the minimum value of the background value in the green direction, and uses this as the threshold to generate a binary image of the background and vegetation, thereby solving the automatic or image classification threshold and achieving the purpose of automatically distinguishing the background and vegetation in the image and generating a binary image.

[0053] 5. The present invention uses a pore and leaf chord length pixel search method to process binary images, and can easily obtain the affiliation of pixels in the binary image, that is, whether the pixel belongs to vegetation or background, so as to calculate the pore size and leaf chord length in a specific direction, or the vegetation width and inter-ridge distance of ridge crops.

[0054] 6. The MX aggregation index equation adopted in the present invention simplifies the difficulty of selecting the search box size in the Lacunarity-Based Aggregation Index mathematically. Therefore, the real leaf area of ​​vegetation and the effective leaf area can be calculated more conveniently, which is not only simple and convenient, but also more practical for image analysis.

[0055] 7. In the calculation of the average blade inclination angle, compared with the existing polynomial technology using empirical algorithms (i.e., statistical methods, the parameters of the equations have no physical meanings, note: the LAI-2000 device uses this method) to calculate the average blade inclination angle, the physical meanings of the physical quantities of the equations derived in the present invention are clearer.

[0056] 8. The equation for the proportion of diffuse radiation to total radiation used in the present invention can effectively calculate the ratio of the incident diffuse radiation of the canopy to the incident total radiation. The ratio of the incident diffuse radiation to the incident total radiation is a key influencing factor in calculating the total absorbed photosynthetic radiation ratio (FAPAR). The solution of the present invention to this problem provides assistance for the effective and correct calculation of FAPAR. At the same time, the calculation method of the time corresponding to the solar altitude angle can facilitate FAPAR photogrammetry.

[0057] 9. The method of the present invention can be used to conveniently measure the cumulative pore size distribution, vegetation width (A1) of ridge crops, distance between ridges (A2), true leaf area index (LAI), effective leaf area index (LAI_e), average leaf inclination angle (ALA), vegetation aggregation index (CI), fractional vegetation cover (FVC), canopy openness (CC) growing season changes and absorbed photosynthetic radiation ratio (FAPAR) and other measurement indicators.

[0058]

Practical application and verification

[0059] ⑴ Compare the vegetation width (A1) of Zhangye Yingke Oasis corn, Ningxia Zhongwei corn, Ningxia Zhongwei wheat, and Ningxia Zhongwei rice, and the ridge spacing (A2) of Zhangye Yingke Oasis corn, Ningxia Zhongwei corn, Ningxia Zhongwei wheat, and Ningxia Zhongwei rice directly with the measurement results of the present invention (see Figure 7 ), it can be found that the vegetation width (A1) obtained by direct measurement method (i.e., direct on-site measurement with a ruler or other equipment) in corn, rice and wheat and the vegetation width (A1) obtained by the measurement method of the present invention have a high consistency, wherein the correlation coefficient (R) is above 0.86291, and the root mean square error (RMSE) is below 10.23675cm. At the same time, it can be found that the measurement of rice is better than that of wheat, and the measurement of wheat is better than that of corn. In the measurements in different regions (i.e. Figure 7 (a)~(d)), there is no obvious regularity. Therefore, the measurement method of the present invention has a higher calculation accuracy.

[0060] (2) Compare the true leaf area index (LAI) of corn measured by the direct measurement method with the true leaf area index (LAI) of corn calculated by the image taken under vertical observation of the present invention (see Figure 8 (a)), it can be found that the root mean square error (RMSE) between the true leaf area index (LAI) measured directly and the true leaf area index (LAI) of the image obtained by the present invention under vertical observation is equal to 1.18944, indicating that the true leaf area index (LAI) measured directly and the true leaf area index (LAI) of the image obtained by the present invention under vertical observation have a high consistency.

[0061] The effective leaf area index (LAI_e) of corn, wheat and rice measured by the LAI-2000 canopy analyzer was compared with the effective leaf area index (LAI_e) of corn, wheat and rice calculated by the image taken at an observation zenith angle of 57.5° in the present invention (see Figure 8(b)~(d)), we can find that: the two have a high consistency. At the same time, we can find that the calculation results of wheat are better than those of corn, and the calculation results of corn are better than those of wheat.

[0062] from Figure 8 The results in (a) to (d) further suggest that the key parameter (aggregation index) required in calculating the true leaf area index (LAI) and the effective leaf area index (LAI_e) is correct. This also indirectly verifies that the MX aggregation index equation proposed in the present invention is correct.

[0063] ⑶ The average leaf inclination angle (ALA) of corn in Zhangye Yingke Oasis was directly measured (see Figure 8 (e) and the average leaf inclination angle (ALA) of Ningxia Zhongwei maize, Ningxia Zhongwei wheat, and Ningxia Zhongwei rice obtained by LAI-2000 canopy analyzer (see Figure 8 (f) to (h)) are compared with the measurement method of the present invention, and it can be found that: they have a high consistency. The consistency between the average leaf inclination angle calculated by the present invention and the average leaf inclination angle obtained by the LAI-2000 canopy analyzer is higher than the average leaf inclination angle obtained by the direct measurement method. In the calculation of the average leaf inclination angle, the measurement of wheat by the present invention is better than that of corn, and the measurement of corn is better than that of rice.

[0064] (4) The aggregation index (CI) of Zhangye Yingke corn, Ningxia Zhongwei corn and Ningxia Zhongwei rice was calculated by the method of the present invention. Fig. 9 As shown in (a), (c) and (e), it can be found that the vegetation aggregation index (CI) increases with the increase of Julian day, which is consistent with the law of vegetation growth. Fig. 9 The Julian days in (c) do decrease after 220 days. The main reason is that the corn is in the late maturity stage. At this time, the corn leaves have gradually turned yellow and shrunk, the pores in the vegetation have become larger, and the degree of aggregation has weakened. Therefore, the aggregation index (CI) shows a decreasing trend at this time.

[0065] The method of the present invention is used to calculate the vegetation coverage (FVC) and canopy openness (CC) of Zhangye Yingke corn, Ningxia Zhongwei corn, and Ningxia Zhongwei rice. Fig. 9 As shown in (b), (d) and (f), it can be found that the vegetation coverage (FVC) and canopy development (CC) show opposite trends. Fig. 9 In (d), the fractional vegetation cover (FVC) shows a decreasing trend after 200 days from the Julian day, while the aggregation index (CI) shows an increasing trend. Fig. 9(c) The reason described in the previous paragraph is the same. This result suggests that the vegetation coverage (FVC) is proportional to the change of the aggregation index (CI), while the canopy development (CC) is inversely proportional to the change of the aggregation index (CI). Fig. 9 (f) It can be found that 200 days after the rice reached the Julian day, the vegetation coverage (FVC) and canopy development (CC) showed saturation.

[0066] (5) Compare the absorbed photosynthetic radiation ratio (FAPAR) of vegetation in Ningxia region directly measured by SpectroSense2 multi-channel spectroradiometer and the absorbed photosynthetic radiation ratio (FAPAR) calculated by the image of the present invention (see Fig. 9 (g)), it can be found that: the two are highly consistent, and the distribution scatter points are all near the 1:1 line position in the figure. Through the study of various vegetation types, it is indirectly proved that the equation of the proportion of diffuse radiation to total radiation proposed in the present invention is correct. This also shows that the present invention can well measure the proportion of absorbed photosynthetic radiation (FAPAR) of vegetation. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] The specific implementation modes of the present invention are further described in detail below in conjunction with the accompanying drawings.

[0068] Figure 1 Schematic diagram of the field of view and optical structure of the camera of the present invention. In the figure: h' is the frame height, v' is the frame width, d is the frame diagonal length, f is the camera focal length, L s is the actual shooting length, W s is the actual shooting width, L v is the vertical shooting length, W v is the vertical shooting width, h is the vertical distance from the shooting position of the camera lens to the top of the canopy, θ o is the zenith angle of the observation direction.

[0069] Figure 2 The schematic diagram of the histogram slope search method and the image coordinate system angle calculation of the present invention. Among them: (a) is a schematic diagram of the A-band pixel distribution frequency in the LAB color space of the soil and background; (b) is a histogram constructed with the pixel value as the horizontal coordinate and the number of pixels corresponding to the value as the vertical coordinate; (c) is a mathematically abstract diagram of the black curve in Figure (b), which shows the slope tanα, and the figure shows the basic idea of ​​the histogram slope search method; (d) is a schematic diagram of the image coordinate system angle calculation.

[0070] In the figure: the gray curve or square in the histogram represents the green value of the A band in the LAB color space, that is, vegetation; and the black curve or square in the histogram represents the green value of the A band in the LAB color space, that is, background. Symbols in the figure: (x1, y1), (x2, y2), (x3, y3), (x4, y4), (x5, y5) and (x6, y6) are points in the coordinate system where the histogram is located, β is the offset angle, α is the angle between the straight line formed by the point (x1, y1) and the point (x2, y2) and the axis (y) formed by the number of pixels, and η is the azimuth offset angle in the image.

[0071] Figure 3 This is a histogram slope search method using wheat, corn and rice as examples. Wherein: (a) is a corn RGB color space image; (b) is a corn binary image; (c) is a wheat RGB color space image; (d) is a wheat binary image; (e) is a rice RGB color space image; (f) is a rice binary image, where the symbol FVC is the English abbreviation for vegetation coverage, and T is a threshold.

[0072] Figure 4 The schematic diagram of the present invention for measuring the absorbed photosynthetic radiation ratio (FAPAR). Among them: (a) is the angle required to calculate the absorbed photosynthetic radiation ratio of the solar diffuse radiation (white sky FAPAR) equation; (b) is the angle selected to calculate the absorbed photosynthetic radiation ratio of the average direct solar radiation of a day (back sky FAPAR). Symbols 1-5 represent the angles selected when shooting.

[0073] Figure 5 The schematic diagram of the leaf chord length pixel search method of the present invention for the truncation selection of continuous vegetation canopy and ridge canopy in the image, where: (a) is a schematic diagram of the truncation of the continuous vegetation canopy in a certain orientation, and (b) is a schematic diagram of the truncation of the ridge canopy in a certain orientation. The gray in the figure represents the truncation of the pore size, and the black line represents the truncation of the leaf chord length.

[0074] Figure 6 The figure is a schematic diagram of blade inclination angle calculation based on blade chord length of the present invention. Wherein: (a) is the three-dimensional situation of the blade, (b) to (g) are the geometric relationship between the three blade chord lengths on the plane. Symbols in the figure: is a reference azimuth, w is The average chord length of the blade corresponding to the direction angle, l * is the average chord length of the blades searched in the horizontal direction of the image coordinates, w * is the average chord length of the blade searched in the horizontal direction of the image coordinate, γ is the sum of w and l * The angle between.

[0075] Figure 7 This is the verification result of the ridge crop vegetation width (A1) and the ridge-to-ridge distance (A2) of the present invention. Among them: (a) is the direct measurement result of the corn vegetation width (A1) of Zhangye Yingke Oasis and the comparison result of the measurement method of the present invention; (b) is the direct measurement result of the corn ridge distance (A2) of Zhangye Yingke Oasis and the comparison result of the measurement method of the present invention; (c) is the direct measurement result of the corn vegetation width (A1) of Ningxia Zhongwei and the comparison result of the measurement method of the present invention; (d) is the direct measurement result of the corn ridge distance (A2) of Ningxia Zhongwei and the comparison result of the measurement method of the present invention; (e) is the direct measurement result of the wheat vegetation width (A1) of Ningxia Zhongwei and the comparison result of the measurement method of the present invention; (f) is the direct measurement result of the wheat ridge distance (A2) of Ningxia Zhongwei and the comparison result of the measurement method of the present invention; (g) is the direct measurement result of the rice vegetation width (A1) of Ningxia Zhongwei and the comparison result of the measurement method of the present invention; (h) is the direct measurement result of the rice ridge distance (A2) of Ningxia Zhongwei and the comparison result of the measurement method of the present invention. Symbols in the figure: D after the horizontal line represents direct measurement, and P represents image measurement.

[0076] Figure 8 The verification results of the real leaf area index (LAI), effective leaf area index (LAI_e) and average leaf inclination angle (ALA) of the vegetation of the present invention are shown in Figure 1. Among them: (a) is the comparison result of corn (LAI) measured by direct measurement method and corn (LAI) calculated by the image taken under vertical observation of the present invention; (b) is the comparison result of corn (LAI_e) measured by LAI-2000 and corn (LAI_e) calculated by the image taken at an observation zenith angle of 57.5° of the present invention; (c) is the comparison result of wheat (LAI_e) measured by LAI-2000 and wheat (LAI_e) calculated by the image taken at an observation zenith angle of 57.5° of the present invention; (d) is the comparison result of rice (LAI_e) measured by LAI-2000. e) is the comparison result of rice (LAI_e) calculated from the image taken at an observation zenith angle of 57.5° in the present invention; (e) is the comparison result of the direct measurement result of the average leaf inclination angle (ALA) of corn in Zhangye Yingke Oasis and the measurement method of the present invention; (f) is the comparison result of the direct measurement result of the average leaf inclination angle (ALA) of corn in Zhongwei, Ningxia and the measurement method of the present invention; (g) is the comparison result of the direct measurement result of the average leaf inclination angle (ALA) of wheat in Zhongwei, Ningxia and the measurement method of the present invention; (h) is the comparison result of the direct measurement result of the average leaf inclination angle (ALA) of rice in Zhongwei, Ningxia and the measurement method of the present invention. Symbols in the figure: D after the horizontal line represents direct measurement, and P represents image measurement.

[0077] Fig. 9The present invention verifies the vegetation aggregation index (CI), vegetation coverage (FVC), canopy openness (CC) growth season changes and the ratio of absorbed photosynthetic radiation (FAPAR) of vegetation. Wherein: (a) is the change of Zhangye Yingke corn aggregation index (CI) calculated by the present invention; (b) is the change of corn vegetation coverage (FVC) and canopy openness (CC) calculated by the present invention in Zhangye area; (c) is the change of Ningxia Zhongwei corn aggregation index (CI) calculated by the present invention; (d) is the change of Ningxia Zhongwei corn vegetation coverage (FVC) and canopy openness (CC) calculated by the present invention; (e) is the change of Ningxia Zhongwei rice aggregation index (CI) calculated by the present invention; (f) is the change of Ningxia Zhongwei rice vegetation coverage (FVC) and canopy openness (CC) calculated by the present invention; (g) is the comparison result of the ratio of absorbed photosynthetic radiation (FAPAR) directly measured by SpectroSense2 multi-channel spectroradiometer and the ratio of absorbed photosynthetic radiation (FAPAR) calculated by the present invention. Symbols in the figure: D after the horizontal line represents direct measurement, and P represents image measurement. DETAILED DESCRIPTION

[0078] A method for measuring vegetation canopy structure based on a non-fisheye digital camera comprises the following steps:

[0079] (1) Determine the CCD or COMS frame size of the non-fisheye digital camera and the vertical distance from the shooting position of the camera lens to the top of the canopy.

[0080] The CCD or COMS frame size of non-fisheye digital cameras has the following parameters: frame height h', frame width v', or frame diagonal length d (see Figure 1 ). These data are provided by the camera when it is manufactured, as shown in Table 1.

[0081] Table 1 Camera reference frames [Making (some) sense out of sensor sizes [EB / OL].

[0082] https: / / www.dpreview.com / articles / 8095816568 / sensorsizes.】

[0083] Frame Name Image width (mm) Image height (mm) 1 / 3 inch sensor 4.8 3.6 One inch sensor 12.8 9.6 4 / 3 inch sensor 18 13.5 APS-C sensor 25.1 16.7 Full-frame sensor 36 24 Medium format sensor 44 33

[0084] Then use equation (1-3) to calculate the field of view FOV of the non-fisheye digital camera h FOV v FOV d :

[0085]

[0086] Where: f is the focal length of the camera; FOV h The field of view corresponding to the frame height; FOV v The field of view corresponding to the frame width; FOV d The field of view corresponding to the diagonal length of the frame.

[0087] Then the actual length and width of the captured vegetation image are calculated.

[0088] The present invention takes the horizontal direction (i.e., the direction of width) as the main shooting direction, and the mathematical principles for other directions (i.e., the direction of height or the direction of diagonal) are the same. Then, the distance from the shooting position of the camera lens to the top of the canopy is determined (which can be determined during measurement), and the actual ground size of the photo can be calculated using equation (4-5).

[0089]

[0090] W s =2tan(0.5×FOV v )×h---------------------(5),

[0091] Where: L s is the actual shooting length, in m; W s is the actual shooting width, in m; h is the distance from the shooting position of the camera lens to the top of the canopy, in m; θ o is the zenith angle of the observation direction, in degrees.

[0092] If the relevant information of the camera is unknown, use equation (6-7) to calculate:

[0093]

[0094] W s =L s ×cosθ o / l w ---------------------(5),

[0095] Where: l w is the aspect ratio of the photo, which can be set in the camera. FOV can be found in the default digital camera library (see Table 2). Where: Equation (6-7) is not as accurate as the calculation based on equation (4-5).

[0096] Table 2 Approximate relationship between focal length and viewing angle of general cameras [Camera lens selection: relationship between camera focal length, field of view angle and depth of field (visible distance) [EB / OL]. https: / / blog.csdn.net / sy95122 / article / details / 80277865.]

[0097] Focal length (f) Field of view (FOV) Focal length (f) Field of view (FOV) Focal length (f) Field of view (FOV) Corner Eyes 180° 14m 114° 20mm 94° 24mm 84° 25m 75° 35m 63° 50mm 46° 70mm 34° 80mm 30° 85m 28°30′ 100mm 24° 135mm 18° 200mm 12° 300mm 8°15′ 400mm 6°10′ 500mm 5° 600mm 4°10′ 800nm 3°5′

[0098] During measurement, in order to enable the camera to capture images at an appropriate height, the non-fisheye digital camera can also be mounted on a telescopic pole or a drone.

[0099] (2) Select the shooting angle and shooting time of the image according to the measured vegetation parameters to perform shooting sampling and obtain the image.

[0100] The conditions for shooting and sampling are:

[0101] When observing vertically, the vegetation parameters that can be measured include cumulative pore size distribution, vegetation width of ridge crops, ridge distance, vegetation coverage, canopy openness, true leaf area index with slight error, effective leaf area index with slight error, and average leaf inclination.

[0102] When the observation zenith angle of 57.5° is used for measurement, the vegetation parameters can be measured as accurate true leaf area index and accurate effective leaf area index [Baret F, De Solan B, Lopez-Lozano R, et al. GAI estimates of row crops from downward looking digital photos taken perpendicular to rows at 57.5° zenith angle: Theoretical considerations based on 3D architecture models and application to wheat crops [J]. Agricultural and Forest Meteorology, 2010, 150 (11): 1393-1401]. This angle is used for measurement because the G function in the infiltration function is equal to 0.5 at this time, and the infiltration function is mathematically solvable.

[0103] When measuring the ratio of absorbed photosynthetic radiation to direct solar radiation (i.e. back sky FAPAR), it is necessary to select the zenith angle of the sun's incidence direction for shooting.

[0104] When measuring the ratio of absorbed photosynthetic radiation to diffuse sky radiation, it is necessary to average the sampling within the 180-degree zenith angle.

[0105] When measuring the ratio of absorbed photosynthetic radiation to the average direct solar radiation in a day, it is necessary to use the time calculation program corresponding to the solar altitude angle to determine the shooting time.

[0106] The shooting time refers to the time corresponding to a certain solar altitude angle which is automatically calculated by the solar altitude angle corresponding time calculation program.

[0107] When calculating the ratio of absorbed photosynthetic radiation to the average direct solar radiation in a day, time series FAPAR data is needed, and the time calculation program corresponding to the solar altitude angle can calculate the time corresponding to the measurement of the time series FAPAR data. First, calculate the altitude angle of the sun in a day, and the equation is as follows:

[0108]

[0109]

[0110] Where: θ ⊙_noon is the solar noon altitude angle; is the local latitude; δ is the solar declination; d n is the Julian Day; sunrise is the local sunrise hour angle; ω sunset is the local sunset hour angle. According to equation (8), the solar altitude angle at noon is calculated, and the solar altitude angle at sunset and sunrise is 0°, then we can calculate Figure 4 (b) The sun altitude of 1-5. At the same time, based on equations (9-11), and then combined with equations (12-13), the sunrise time (T sunrise ) and sunset time (T sunset ).

[0111]

[0112] Considering that the sun's altitude angle in the text is 12 noon, and the time of sunset and sunrise is known, we can calculate Figure 4 (b) The time corresponding to the solar altitude of 1-5. Equation (8-13) serves equation (34) to provide it with the measured time and measured altitude angle in the average process of one day.

[0113] ⑶ Perform Gamma transformation on the RGB color space of the image to correct the color grayscale difference caused by incorrect image exposure. Then convert the RGB color space image to LAB color space. The specific process is as follows:

[0114] First, use the Gamma function to perform simple grayscale processing on the image to deal with abnormal exposure. The Gamma function is equation (14-17).

[0115]

[0116] In the formula: r is the red band captured by the camera, without unit; R is the red band after Gamma transformation, without unit; g is the green band captured by the camera, without unit; G is the green band after Gamma transformation, without unit; b is the blue band captured by the camera, without unit; B is the red band after Gamma transformation, without unit; x is a variable, which is the value in the brackets in equation (14-16).

[0117] Then, equation (18) is used to convert the image in the RGB color space after the Gamma transformation into an image in the XYZ color space.

[0118]

[0119] The XYZ color space is a color transition space, and then the XYZ color space is converted into the LAB color space using equations (19-22).

[0120] L = 116 × f (Y / Y n )-16---------------------(19),

[0121] A=500×[f(X / X n )-f(Y / Y n )]---------------------(20),

[0122] B=200×[f(Y / Y n )-f(Z / Z n )]---------------------(twenty one),

[0123]

[0124] Where: X n Equal to 95.047, Y n =100.0, Z n = 108.883 (these are in D 65 Normalized value under light source). L is the brightness from black to white, unitless; A is the brightness from green to red, unitless; B is the brightness from blue to yellow, unitless; t is a variable, which is the value in parentheses in equation (19-21), unitless; f(t) is the change function.

[0125] (4) The A band in the LAB color space reflects the brightness of green and red. Therefore, this band can be used to distinguish green vegetation from background information. Based on this principle, the A band in the LAB color space is preprocessed to obtain a binary image of the background and vegetation.

[0126] Wherein: the preprocessing includes the following steps:

[0127] ① Calculate the frequency of the A band value in the LAB color space, and construct a frequency histogram with the pixel value as the horizontal coordinate (x) and the number of pixels corresponding to the value (y) as the vertical coordinate;

[0128] ② Calculate the derivative f'[N(x)] of the frequency function f[N] composed of the number of pixels (y) of the A band value in the LAB color space, and calculate the extreme value in the frequency histogram so that the extreme value satisfies the requirements of {f'[N(x)]>0}∧{f'[N(x)]>0}. Then, take the second local minimum in the positive direction of the coordinate axis (x) of the pixel value as the intersection of the Gaussian distribution in the frequency of the vegetation value and the Gaussian distribution in the frequency of the background value. The area near the intersection is the area where the threshold value to be determined is located.

[0129] ③ Based on the intersection point and the points (x1, y1) and (x2, y2) near the intersection point, use the histogram slope search method equation to search for the minimum value of the Gaussian distribution in the frequency of the background value toward the green direction, and use this minimum value as the threshold;

[0130] ④ Classify the vegetation and background in the A band of the LAB color space according to the threshold value to obtain binary images of the background and vegetation; the binary image 255 corresponds to the vegetation, and 0 corresponds to the background.

[0131] The specific process is as follows:

[0132] Since in the A band, the frequency of the number of pixels corresponding to the soil and vegetation values ​​is a double Gaussian distribution (i.e. Figure 2 (a) [Yaokai Liu XM, Haoxing Wang & Guangxian Yan. A novel method for extracting green fractional vegetation cover from digital images [J]. Journal of Vegetation Science, 2012, 23 (3): 406-418]. To this end, a histogram can be constructed with the pixel value as the horizontal axis and the number of pixels corresponding to the value as the vertical axis ( Figure 2 )(b)) to determine the specific location of the classification threshold. In the process of calculating the classification threshold, we must first determine Figure 2The extreme value of the distribution curve in (b) is obtained by taking the derivative of the distribution curve, that is, equation (23):

[0133]

[0134] Where: f'[N(x)] is the derivative of the double Gaussian distribution curve of soil and vegetation in band A; f[N(x)] is the first discrete value of the minimum value of the double Gaussian distribution curve of soil and vegetation in band A in the positive direction, that is: Figure 2 (b) The value of y2 at the midpoint (x2,y2); f[N(x + )] is the second discrete value in the positive direction of the minimum value of the double Gaussian distribution curve of soil and vegetation in band A, that is: Figure 2 (b) the value of y1 at the midpoint (x1, y1);

[0135] If equation (23) satisfies {f'[N(x)]>0}∧{f'[N(x)]>0}, then the local minimum of the soil and vegetation double Gaussian distribution curve is Figure 2 Point (x6, y6) in (a).

[0136] Since the frequency of the histogram image element value is only a discrete function rather than a continuous function, the local minimum cannot be obtained, but the approximate (x2, y2) point can be obtained, and then the local minimum value can be obtained. Figure 2 The coordinates of point (x1, y1) and point (x2, y2) in (c) are used to calculate the slope tana of the straight line corresponding to these two points, and then the point (x4, y4) is calculated; again, using point (x2, y2) and point (x4, y4), the midpoint (x3, y3) of these two points on the straight line is calculated, and finally, the coordinate equation of tanα calculated previously is used to consider an offset slope tanβ, and the minimum value of the estimated soil Gaussian distribution curve in the green direction is calculated together, that is, point (x5, y5).

[0137] The above calculation steps are converted into mathematics, and the following equation is obtained:

[0138]

[0139] D=y2-y1+tanβ(x1-x2)---------------------(25),

[0140] E=x1y2-x2y2---------------------(26),

[0141] F=x2-x1-tanβ(y1-y2)---------------------(27),

[0142]

[0143] Where: T is the threshold; D, E, and F are the functions in the threshold T equation respectively; point (x1, y1) and point (x2, y2) are points near the focus in the frequency histogram; β is the offset angle; α is the angle between the straight line formed by the point (x1, y1) and the point (x2, y2) and the axis (y) formed by the number of pixels; τ is the control coefficient of tanα, and the default value is 0.7.

[0144] The calculated threshold T is approximately -2 to -4 ( Figure 2 (a)). This threshold can well extract vegetation from the image, complete the threshold automatic classification, and obtain the binary image of the background and vegetation. In the binary image value, the vegetation part value is set to 255, and the background part value is set to 0. The calculation result of the histogram slope search method is as follows Figure 3 shown.

[0145] (5) Calculate the ratio of the background and vegetation values ​​in the binary image, and calculate the vegetation coverage (FVC), canopy openness (CC), porosity, the ratio of absorbed photosynthetic radiation of direct solar radiation (i.e., back sky FAPAR) and the ratio of absorbed photosynthetic radiation of diffuse sky radiation (i.e., white sky FAPAR); at the same time, calculate the average absorbed photosynthetic radiation ratio of direct solar radiation in a day according to the shooting time corresponding to the absorbed photosynthetic radiation ratio of direct solar radiation. For the total absorbed photosynthetic radiation ratio, the ratio of incident diffuse radiation to incident total radiation is calculated according to the equation of diffuse radiation to total radiation ratio. The calculated absorbed photosynthetic radiation ratio of direct solar radiation and the absorbed photosynthetic radiation ratio of diffuse sky radiation are combined to obtain the instantaneous total absorbed photosynthetic radiation ratio (i.e., FAPAR).

[0146] If the photo is taken at a vertical angle, the vegetation coverage (FVC) and canopy openness (CC) in the image can be calculated. The calculation equation is as follows:

[0147]

[0148] Where: N v is the number of vegetation pixels in the image; N b is the number of background pixels in the image. These two parameters are obtained by searching the pore and blade chord length pixel search methods mentioned below.

[0149] The changing trend of equation (29-30) is as follows Fig. 9 (b) and Fig. 9 (d).

[0150] If the photo is taken at a specific angle, the porosity P at that specific angle can be calculated.o (θ o ), that is, equation (31).

[0151]

[0152] In the FAPAR calculation, equation (32) can calculate the instantaneous back sky FAPAR (i.e., the ratio of absorbed photosynthetic radiation to direct solar radiation, which is symbolized as fAPAR bs ), and equation (32) can calculate the instantaneous white sky FAPAR (i.e., the ratio of absorbed photosynthetic radiation to direct solar radiation in the diffuse sky radiation, fAPAR ws )

Richter C, Gueymard CA, Lincot D.Solar Energy[M].New York: Springer, 2019

[0153]

[0154] Where: P o (θ s ) is the porosity observed in the sun direction, which is calculated by equation (31).

[0155] Equation (33-34) is an integral. In the present invention, five photos taken at different angles are used to calculate and approximate this integral. In equation (33), the angle is selected by Figure 4 (a). The angle of equation (34) is chosen by Figure 4 (b).

[0156] The time calculation procedure corresponding to the solar altitude angle is the same as step (2).

[0157] Considering direct solar radiation and diffuse sky radiation, the total FAPAR is

[0158] fAPAR=(1-skly)×fAPAR bs +skly×fAPAR ws --------(35)

[0159]

[0160] Where: sky is the proportion of diffuse radiation at the top of the canopy to the total radiation. The equation is derived from the equation in the Gap Light Analyzer (GLA) manual [Frazer GW, Canham CD, Lertzman K P. Gap Light Analyzer (GLA): Imaging software to extract canopy structure and gap light transmission indices from true-colour fisheye photographs, users manual and program documentation.].

[0161] R p / R s It represents the ratio of photosynthetically active radiation to shortwave radiation. According to research, it can be obtained by looking up a table or directly measured. For details, please refer to the literature [Papaioannou G, Nikolidakis G, Asimakopoulos D, et al. Photosynthetically active radiation in Athens [J] Agricultural & Forest Meteorology, 1996, 81 (3-4): 0-298], where the verification results of FAPAR are as follows Fig. 9 (g) as shown.

[0162] (6) Based on the binary image, the pore and leaf chord length pixel search method is used to calculate the leaf chord length and pore size in the vertical direction of the image, the leaf chord length and pore size in the horizontal direction, and the leaf chord length and pore size in a specific direction. Then the cumulative pore size distribution in each pixel direction is calculated; at the same time, the vegetation width and inter-ridge distance of the ridge crop are calculated based on the pore and leaf chord length pixel search method.

[0163] Among them: Taking the leaf chord length and pore size in the horizontal direction of the image (i.e., the length direction of the photo) as an example, the pore and leaf chord length pixel search method refers to an algorithm for searching the pixel affiliation in the image, that is, whether the pixel belongs to vegetation or background. The algorithm expression is:

[0164] for i=1,L image

[0165] for j=1,W image

[0166]

[0167] In equation (37): = is the mathematical symbol “defined as”; w is the blade chord length; λ is the pore size length; x and y are the pixel coordinate positions in the image coordinates (see Figure 2 (d)); L image is the number of pixels in the image length; W image is the number of pixels in the image width; P is the pixel value. The search conditions are:

[0168]

[0169] Figure 2 (d) shows a schematic diagram of angle search. Taking the horizontal direction as the starting point, the image coordinate value in the angle search is calculated using the following equation:

[0170]

[0171] y i+1 =y i +1------------(40),

[0172] Where: η is the azimuth offset angle in the image, that is: Figure 2 (d) Note: The starting point of the image in the calculation of blade chord length and pore size (i.e. point (x1, y1)) is Figure 5 (a)) and the end point (i.e., point (x2, y2) Figure 5 (a)) is also calculated using equations (39-40).

[0173] in Figure 5 (a) and (b) are only for reference. The measurer can choose to encrypt the truncation according to his own requirements. Symbols in the figure: (x1, y1) and (x2, y2) are the starting point and end point of the pore size calculation in the image coordinate system (here only the pore size coordinate acquisition is used as an example), (x s ,y s ) is the starting pixel coordinate of the row, column or specific orientation in the image coordinate system, (x e ,y e ) is the end pixel coordinate of the row, column or feature position in the image coordinate system.

[0174] Equation (37-40) is an example of the calculation of pore size and leaf chord length. Based on the same mathematical principle, the vegetation width A1 and the ridge distance A2 of ridge crops can be calculated (Note: ridge crops refer to the planting form in which farmland crops are separated by soil during planting. The separation interval is called ridge interval, and the interval distance is called ridge distance. Similarly, the width of the separated vegetation part is called vegetation width, see Figure 5 (b)). The verification results can be found in Figure 7 .

[0175] Based on the pore size, the average cumulative pore size distribution function in the image is calculated as:

[0176]

[0177] Where: F(λ) is the average cumulative pore size distribution function in all directions; is the cumulative pore size distribution function in a certain direction; N t The number of pores calculated for the pixels in a certain direction of the image; λ max is the maximum pore size in a certain direction; (x s ,y s ) is the starting pixel coordinate of the row, column or specific position in the image coordinate system, (x e ,y e ) is the pixel coordinate of the end point in the row, column or specific orientation in the image coordinate system. The calculation method can be found in Figure 5 .

[0178] ⑺According to the cumulative pore size distribution in each pixel direction, the MX aggregation index equation is used to calculate the aggregation index of vegetation; then, based on the aggregation index, the true leaf area index with a slight error and the effective leaf area index with a slight error in the image taken in the vertical direction, as well as the accurate true leaf area index and effective leaf area index in the image taken at an observation zenith angle of 57.5° are calculated.

[0179] Based on the previous research on heterogeneity in variance and expectation [Dale MR T. Lacunarity analysis of spatial pattern: A comparison [J]. Landscape Ecology, 15 (5): 467-478.], the present invention proposes the MX aggregation index equation based on the cumulative pore size of vegetation:

[0180]

[0181] Where: Ω E is the MX aggregation index; Λ(λ) is the heterogeneity index between all pixels in a single row or column of the image; Λ max (λ) is the maximum value of Λ(λ) in column or row statistics, and Λ max (λ) = L image , L image is the number of image lengths; Λ cri (λ) is the value when the heterogeneity index between all pixels in a single row or column of the image is randomly distributed, and Λ cri (λ)=1;σ 2(λ) is the statistical variance between all pixels in a single row of the image; E(λ) is the statistical expectation between all pixels in a single row of the image; λ i The pore size calculated for the pixels in a single row or column in the ith image; is the statistical average value of the pore size in a single row calculated based on the pixels in a single row or column in the image; N is the number of pores calculated based on the pixels in a single row or column in the image; F(λ i ) is λ i The corresponding pore size cumulative distribution function value.

[0182] In equation (42), it is analogous to the normalization method of the Lacunarity-Based Aggregation Index, where A max (λ) = L image , A cri (λ)=1. Their determination refers to the literature [Fournier RA, Hall R J. Hemispherical Photography in Forest Science: Theory, Methods, Applications [M]. Van Godewijckstraat: Springer, 2017.], and its changing trend is as follows Fig. 9 (a) and Fig. 9 (c)).

[0183] According to the aggregation index, the real leaf area index and the effective leaf area index can be calculated. According to F. Baret's research [Baret F, De Solan B, Lopez-Lozano R, et al. GAI estimates of row crops from downward looking digital photos taken perpendicular to roWs at 57.5° zenithangle: Theoretical considerations based on 3D architecture models and application to wheat crops [J]. Agricultural and Forest Meteorology, 2010, 150 (11): 1393-1401.], when the observation zenith angle is 57.5°, the G function is equal to a constant of 0.5, and the penetration function for calculating the leaf area index can be accurately calculated.

[0184] The true leaf area index (LAI) and effective leaf area index (LAI e ) is calculated according to equation (44-45):

[0185]

[0186] LAI e =LAI×Ω E -------------------(45),

[0187] Where: P o (57.5°) is the porosity when the observation zenith angle is 57.5°. It is the value measured by the digital camera when the observation zenith angle is 57.5° and can be calculated by equation (31).

[0188] Equation LAI e The verification results are as follows Figure 8 (b)~(d). According to the research of J.Liu, when the observation zenith angle is 0° (i.e. vertical observation), the extinction coefficient is approximately 0.61 [Liu J, Pattey E, Admiral S. Assessment of insitu crop LAI measurement using unidirectional view digital photography [J]. Agricultural and Forest Meteorology, 2013, 169: 25-34], and the penetration function for calculating the leaf area index can still be approximately calculated. At this time, when observing vertically, the true leaf area index is calculated as follows:

[0189]

[0190] Where: P o (0°) is the porosity when the observation zenith angle is 0°. It is the value measured by a digital camera when the observation zenith angle is 0° (ie, vertical observation) and can be calculated by equation (31).

[0191] At this time, substitute equation (46) into equation (45) to calculate the effective leaf area index of vertical observation. Note that the true leaf area index calculated by equation (46) has a small error. Therefore, in the present invention, it is called the true leaf area index with a slight error, and the effective leaf area index calculated is called the effective leaf area index with a slight error. Although the calculation result has an error, it basically meets the measurement requirements. The verification result is as follows Figure 8 (a).

[0192] ⑻Based on the leaf chord length in the vertical direction, the leaf chord length in the horizontal direction, and the leaf chord length in a specific orientation, the average leaf inclination angle equation can be used to calculate the average leaf inclination angle of the vegetation in the image.

[0193] In the study of vegetation "hot spot" factors, an equation for calculating leaf chord length was proposed [Ma X, Lu L. Mathematical analysis on the component of canopy architecture in rowcrops. 2019]:

[0194]

[0195] Where: w is the blade chord length; P ol is the blade distribution form; b is the blade bending factor; l * is leaf length; w * is the leaf width; θ l is the average leaf inclination angle; is the observation azimuth.

[0196] If it is vertical observation, P ol is the distribution pattern of leaves on the stem, which is similar to the aggregation index Ω E , so let P ol =Ω E , then the average blade inclination angle equation is:

[0197]

[0198] In the formula yes The reason for this is that it is difficult to obtain the exact leaf length and width in the image. is a reference direction. The average chord length of the blade corresponding to this direction is w, and the average chord length of the blade searched in the horizontal direction of the image coordinate (i.e., the direction of length) is l * The average chord length of the blade searched in the vertical direction of the image coordinates (i.e., the width direction) is w * γ is the sum of w and l * The angle between them. The various cases of equation (48) are as follows Figure 6 (b) to (g) are shown. At the same time, the verification results are as follows Figure 8 As shown in (e) to (h).

Claims

1. A method for measuring vegetation canopy structure based on a non-fisheye digital camera, comprising the following steps: ⑴ Determine the CCD or COMS frame size of the non-fisheye digital camera and the vertical distance from the shooting position of the camera lens to the top of the canopy, calculate the field of view of the non-fisheye digital camera according to the corresponding mathematical equation, and then calculate the actual length and width of the vegetation image taken; (2) Selecting the shooting angle and shooting time of the image according to the measured vegetation parameters to perform shooting sampling to obtain the image; the shooting time refers to the time corresponding to a certain solar altitude angle automatically calculated according to the solar altitude angle corresponding time calculation program; ⑶ Performing Gamma transformation on the RGB color space of the image, and then converting the transformed RGB color space image into LAB color space; (4) Preprocessing the A band in the LAB color space to obtain binary images of the background and vegetation; (5) Counting the ratio of the background and vegetation values ​​in the binary image, and calculating the vegetation coverage, canopy openness, porosity, the ratio of absorbed photosynthetic radiation of direct solar radiation and the ratio of absorbed photosynthetic radiation of diffuse sky radiation; At the same time, according to the shooting time corresponding to the proportion of absorbed photosynthetic radiation of direct solar radiation, the average proportion of absorbed photosynthetic radiation of direct solar radiation in a day is calculated, and the proportion of diffuse radiation to total radiation is calculated according to the equation of the proportion of diffuse radiation to total radiation, and then the instantaneous FAPAR is calculated; (6) According to the binary image, the leaf chord length and pore length in the vertical direction of the image, the leaf chord length and pore length in the horizontal direction, and the leaf chord length and pore length in a specific direction are calculated by using the pore and leaf chord length pixel search method; and then the cumulative pore size distribution in each pixel direction is calculated; at the same time, the vegetation width and inter-ridge distance of the ridge crops are calculated according to the pore and leaf chord length pixel search method; (7) According to the cumulative pore size distribution in each pixel direction, the aggregation index of the vegetation is calculated using the MX aggregation index equation; and then the real leaf area index with a slight error and the effective leaf area index with a slight error in the image taken in the vertical direction, as well as the accurate real leaf area index and the effective leaf area index in the image taken at an observation zenith angle of 57.5° are calculated based on the aggregation index; ⑻ Based on the blade chord length in the vertical direction, the blade chord length in the horizontal direction, and the blade chord length in a specific orientation, the average leaf inclination angle equation can be used to calculate the average leaf inclination angle of the vegetation in the image.

2. The method of measuring vegetation canopy structure based on a non-fisheye digital camera as claimed in claim 1, characterized in that: In the step (1), the non-fisheye digital camera is mounted on a telescopic rod or a drone.

3. The method of measuring vegetation canopy structure based on a non-fisheye digital camera as claimed in claim 1, characterized in that: The conditions for shooting and sampling in the step (2) refer to: when vertical shooting is adopted, the measured vegetation parameters are the cumulative pore size distribution, the vegetation width of the ridge crops, the ridge distance, the vegetation coverage rate, the canopy openness, the true leaf area index with a slight error, the effective leaf area index with a slight error, and the average leaf inclination angle; when the observation zenith angle is 57.5°, the measured vegetation parameters are the accurate true leaf area index and the accurate effective leaf area index; when measuring the ratio of absorbed photosynthetic radiation of direct solar radiation, shooting is performed at the zenith angle of the sun's incidence direction; when measuring the ratio of absorbed photosynthetic radiation of diffuse radiation from the sky, it is necessary to perform average sampling within 180 degrees of the zenith angle; when measuring the ratio of absorbed photosynthetic radiation of the average direct solar radiation in a day, it is necessary to use the time calculation program corresponding to the solar altitude angle to determine the shooting time.

4. The method of measuring vegetation canopy structure based on a non-fisheye digital camera as claimed in claim 1, characterized in that: The solar altitude angle corresponding time calculation procedure in step (2) refers to the combination of equations for calculating the solar altitude angle and the time corresponding to the solar altitude angle, and the equation is as follows: Where: θ ⊙_noon is the solar noon altitude angle; is the local latitude; δ is the solar declination; d n is the Julian Day; sunrise is the local sunrise hour angle; ω sunset is the local sunset hour angle; T sunrise is the local sunrise time; T sunset The local sunset time.

5. The method of measuring vegetation canopy structure based on a non-fisheye digital camera as claimed in claim 1, characterized in that: The pretreatment in the step (4) comprises the following steps: ① Calculate the frequency of the A band value in the LAB color space, and construct a frequency histogram with the pixel value as the horizontal coordinate (x) and the number of pixels corresponding to the value (y) as the vertical coordinate; ② Calculate the derivative f'[N(x)] of the frequency function f[N] composed of the number of pixels (y) of the A band value in the LAB color space: And calculate the extreme value in the frequency histogram so that the extreme value satisfies the requirement {f'[N(x)]>0}∧{f'[N(x)]>0}; Where: f'[N(x)] is the derivative of the double Gaussian distribution curve of soil and vegetation in band A; f[N(x)] is the first discrete value of the minimum value of the double Gaussian distribution curve of soil and vegetation in band A in the positive direction; f[N(x + )] is the second discrete value in the positive direction of the minimum value of the double Gaussian distribution curve of soil and vegetation in band A; Then, the second local minimum value of the coordinate axis (x) of the pixel value in the positive direction is taken as the intersection point of the Gaussian distribution in the frequency of the vegetation value and the Gaussian distribution in the frequency of the background value, and the area near the intersection point is the area where the threshold value to be determined is located; ③ Based on the intersection point and the points (x1, y1) and (x2, y2) near the intersection point, use the histogram slope search method equation to search for the minimum value of the Gaussian distribution in the frequency of the background value toward the green direction, and use this minimum value as the threshold; The histogram slope search method equation is as follows: D=y2-y1+tanβ(x1-x2),E=x1y2-x2y2,F=x2-x1-tanβ(y1-y2), Where: T is the threshold; D, E, and F are the functions in the threshold T equation respectively; point (x1, y1) and point (x2, y2) are points near the focus in the frequency histogram; β is the offset angle; τ is the control coefficient of tanα, and α is the angle between the straight line formed by the point (x1, y1) and the point (x2, y2) and the axis (y) formed by the number of pixels; ④ Classify the vegetation and background in the A band of the LAB color space according to the threshold value to obtain binary images of the background and vegetation; wherein the binary image 255 corresponds to the vegetation, and 0 corresponds to the background.

6. The method of measuring vegetation canopy structure based on a non-fisheye digital camera as claimed in claim 1, characterized in that: The diffuse radiation to total radiation ratio equation in step (5) refers to the ratio of the diffuse radiation of the sky to the total radiation, and the equation is calculated as follows: Where: sky is the ratio of the diffuse radiation incident on the top of the canopy to the total incident radiation; R p / R s It represents the ratio of photosynthetically active radiation to shortwave radiation.

7. The method of measuring vegetation canopy structure based on a non-fisheye digital camera as claimed in claim 1, characterized in that: The pore and leaf chord length pixel search method in step (6) refers to an algorithm for searching the pixel affiliation in an image, that is, whether the pixel belongs to vegetation or background; the algorithm expression is: for i=1,L image for j=1,W image Where: = is the mathematical symbol "defined as"; w is the blade chord length; λ is the pore size length; x and y are the pixel coordinate positions in image coordinates; L image is the number of image lengths; W image is the number of image widths; P is the pixel value; the search conditions are: Taking the horizontal direction as the starting point, the image coordinate values ​​in the angle search are calculated using the following equation: y i+1 =y i +1, where: η is the offset angle of the orientation in the image.

8. The method of measuring vegetation canopy structure based on a non-fisheye digital camera as claimed in claim 1, characterized in that: The MX aggregation index equation in step (7) is an indicator that describes the aggregation degree of leaves in the image, and its equation form is: Where: Ω E is the MX aggregation index; Λ(λ) is the heterogeneity index between all pixels in a single row or column of the image; Λ max (λ) is the maximum value of Λ(λ) in column or row statistics, and Λ max (λ) = L image , L image is the number of image lengths; Λ cri (λ) is the value when the heterogeneity index between all pixels in a single row or column of the image is randomly distributed, and Λ cri (λ)=1;σ 2 (λ) is the statistical variance between all pixels in a single row of the image; E(λ) is the statistical expectation between all pixels in a single row of the image; i The pore size calculated for the pixels in a single row or column in the ith image; N is the statistical average value of the pore size in a single row calculated based on the pixels in a single row or column in the image; t The number of pores calculated for pixels in a single row or column of the image; F(λ i ) is λ i The corresponding pore size cumulative distribution function value.

9. The method of measuring vegetation canopy structure based on a non-fisheye digital camera as claimed in claim 1, characterized in that: The average leaf inclination angle equation in step (8) refers to calculating the average value of the leaf inclination angle of the vegetation canopy, and the equation is: Where: θ l is the average leaf inclination angle; b is the bending factor of the blade; is a reference position; w is the average chord length of the blade corresponding to the reference orientation; l * w is the average chord length of the blades searched in the horizontal direction of the image coordinates; * is the average chord length of the leaves searched in the vertical direction of the image coordinates; γ is the sum of w and l * The angle between E is the MX aggregation index.

Images