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

Abstract

The present application relates to a kind of non-fisheye digital camera-based vegetation canopy structure measurement method, comprising the following steps: (1) the field of view of non-fisheye digital camera and the actual length and width corresponding to image are calculated;(2) shooting sample is carried out, and image is obtained;(3) the RGB color space of image is converted to LAB color space;(4) the A wave band in LAB color space is preprocessed, and the binary image of background and vegetation is obtained;(5) the FVC, CC, porosity, solar direct radiation and the FAPAR of sky diffuse radiation, the FAPAR of average solar direct radiation in a day, the proportion of diffuse radiation in total radiation, and then the instantaneous FAPAR is calculated;(6) the cumulative pore size distribution of each pixel direction, the vegetation width of ridge row crop and the distance between ridges are calculated;(7) CI, slightly error LAI and LAI_e, and accurate LAI and LAI_e are calculated;(8) average leaf angle is calculated.The present application has low labor cost, no image distortion and wide applicability.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of vegetation canopy structure research and quantitative remote sensing of vegetation, and particularly relates to a method for measuring vegetation canopy structure based on a non-fisheye digital camera. BACKGROUND

[0002] The research of vegetation structure measurement based on photos is a camera optical measurement method, which is derived from the mathematical research of vegetation canopy structure. In the research of vegetation canopy structure, J. Ross, a member of the USSR Academy of Sciences, described the relationship between vegetation canopy structure and radiation in detail 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

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.

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.

[0003] Subsequent scholars A. B. G. Lang and J. M. Chen respectively studied the 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 J M, Cihlar J. Plant canopy gap-size analysis theory for improving optical measurements of leaf-area index [J]. Applied Optics, 1995, 34(27): 6211-6222.

Fournier R A, Hall RJ. Hemispherical Photography in Forest Science: Theory, Methods, Applications [M]. Van Godewijckstraat: Springer, 2017.

【Dale M R T. Lacunarity analysis of spatial pattern: A comparison[J]. Landscape Ecology, 15(5): 467-478.

Plotnick RE, Gardner R H, Hargrove W W, 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.

[0004] Therefore, there is an urgent need for a non-deformation imaging measurement and a measurement method with less human influence to solve the above problems. SUMMARY

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

[0006] To solve the above problems, the vegetation canopy structure measurement method based on a non-fisheye digital camera comprises the following steps:

[0007] (1) determining 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, calculating the field of view of the non-fisheye digital camera according to the corresponding mathematical equation, and further calculating the actual length and width corresponding to the shot vegetation image;

[0008] (2) selecting the shooting angle and shooting time of the image according to the measured vegetation parameters to take samples and obtain the image; the shooting time refers to the time corresponding to a certain solar elevation angle calculated by the solar elevation angle corresponding time calculation program;

[0009] (4) performing Gamma transformation on the RGB color space of the image, and then converting the changed RGB color space image into LAB color space;

[0010] (5) preprocessing the A band in the LAB color space to obtain the binary image of the background and the vegetation;

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

[0012] (7) calculating the leaf chord length and pore length in the vertical direction, the leaf chord length and pore length in the horizontal direction, and the leaf chord length and pore length in a specific direction according to the binary image by using the pore and leaf chord length pixel search method; further calculating the cumulative pore size distribution in each pixel direction; at the same time, calculating the vegetation width and inter-row distance of the ridge row crop by using the pore and leaf chord length pixel search method;

[0013] ⑺According to the cumulative pore size distribution of each pixel direction, the aggregation index of vegetation is calculated by using the M-X aggregation index equation; and then the slightly error real leaf area index and the slightly error effective leaf area index in the image taken in the vertical direction, and the accurate real leaf area index and effective leaf area index in the image taken at the observation zenith angle of 57.5° are calculated according to the aggregation index;

[0014] ⑻According to the vertical leaf chord length, the horizontal leaf chord length and the leaf chord length in a specific direction, the average leaf inclination angle of the vegetation in the image is calculated by using the average leaf inclination angle equation.

[0015] In the step ⑴, the non-fisheye digital camera is loaded on the telescopic rod or the unmanned aerial vehicle.

[0016] In the step ⑵, the sampling conditions are as follows: when the vertical shooting is adopted, the measured vegetation parameters are the cumulative pore size distribution, the vegetation width of ridge-row crops, the distance between ridges, the vegetation coverage, the canopy openness, the slightly error real leaf area index, the slightly error effective leaf area index and the average leaf inclination angle; when the observation zenith angle is 57.5°, the measured vegetation parameters are the accurate real leaf area index and the accurate effective leaf area index; when the absorption photosynthetic radiation ratio of solar direct radiation is measured, the shooting is performed at the zenith angle of the solar incident direction; when the absorption photosynthetic radiation ratio of sky diffuse radiation is measured, the average sampling is required within the zenith angle of 180°; when the average absorption photosynthetic radiation ratio of solar direct radiation in a day is measured, the shooting time is determined by using the solar elevation angle corresponding time calculation program.

[0017] The solar elevation angle corresponding time calculation program in the step ⑵ refers to the equation combination for calculating the solar elevation angle and the time corresponding to the solar elevation angle, and the equation is as follows:

[0018]

[0019] In the formula, θ ⊙_noon is the solar noon elevation angle; is the measured latitude of the place; δ 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 sunrise time of the place; T sunset is the sunset time of the place.

[0020] The pre-processing in the step ⑷ includes the following steps:

[0021] ①Calculate the frequency of A waveband 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 pixel 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 A waveband 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};

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

[0024] Then, take the second local minimum value of the coordinate axis (x) of the pixel value in the positive direction 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 region near the intersection point is the region 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, search for the minimum value of the Gaussian distribution in the frequency of the background value in the green direction using the histogram slope search equation, and take the minimum value as the threshold value;

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

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

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

[0029] In the formula: T is the threshold value; D, E, and F are functions in the threshold value T equation; the points (x1, y1) and (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 composed of the points (x1, y1) and (x2, y2) and the axis (y) composed of the number of pixels;

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

[0031] The equation of the proportion of diffuse radiation in the total radiation in step 8 is the proportion of the sky diffuse radiation in the total radiation, and the equation is calculated as follows:

[0032]

[0033] In the formula, sky is the ratio of the diffuse radiation incident at the top of the canopy to the total incident radiation; R p / R s represents the ratio of the photosynthetically active radiation to the shortwave radiation.

[0034] The pore and leaf chord length pixel search method in step 6 is a kind of algorithm for searching the membership relationship of pixels in the image, that is, whether the pixel belongs to vegetation or background; and the algorithm expression is as follows:

[0035] for i=1,L image

[0036] for j=1,W image

[0037]

[0038] In the formula, = is the mathematical symbol "defined as"; w is the leaf 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; and the search condition is:

[0039]

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

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

[0042] The M-X aggregation index equation in step 8 is an index for describing the aggregation degree of leaves in the image, and the equation form is as follows:

[0043]

[0044] In the formula, Ω E is the M-X aggregation index; Λ(λ) is the heterogeneity index between all pixels in a single row (or single column) of the image; Λ max (λ) is the maximum value of Λ(λ) in the column (or row) statistics, and Λ max (λ)=L image , L imageIt is the number of images of varying lengths; Λ cri (λ) represents the value of the heterogeneity index among all pixels in a single row (or column) of an image when it exhibits a random distribution, and Λ cri (λ)=1;σ 2 (λ) represents the statistical variance among all pixels in a single row of the image; E(λ) represents the statistical expectation among all pixels in a single row of the image; λ i The aperture size calculated for a single row (or column) of pixels in the i-th image; N is the statistical average value of the aperture size calculated based on the pixels in a single row (or column) of the image, within that single row; t The number of apertures calculated for a single row (or column) of an image; F(λ) i ) is λ i The corresponding cumulative distribution function value of pore size.

[0045] The equation for the average leaf tilt angle in step (8) refers to the calculation of the average tilt angle of the leaves in the vegetation canopy, and its equation is:

[0046]

[0047] In the formula: θ l denoted as the average leaf tilt angle; b is the leaf bending factor. Let w be a reference orientation; w is the average chord length of the blade corresponding to the reference orientation; l * The average chord length of the blade is searched in the horizontal direction of the image coordinates (i.e., the direction of length); w * The average chord length of the blade is searched in the vertical direction of the image coordinates (i.e., the direction of width); γ is the average chord length of w and l. * The included angle between them; Ω E MX is the aggregation index.

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

[0049] 1. This 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 membership relationships. This solves the problem of automatically extracting vegetation and background classification thresholds from images and performing preliminary vegetation structure measurements in photogrammetry. Furthermore, the invention proposes methods for calculating the MX vegetation aggregation index, planar leaf tilt angle, and the required location parameters for Fabricated Photosynthetic Radiation Ratio (FAPAR), solving the problem of calculating vegetation structure parameters. Thus, it realizes a method for measuring vegetation canopy structure based on a non-fisheye digital camera.

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

[0051] 3. This invention allows a camera to be mounted on a drone for high-altitude photography. Due to the wide field of view (FOV) of drones, it is easy to acquire overall images of the study area, thus facilitating the analysis of average information on the vegetation canopy structure. Compared to previous single-point measurement algorithms for obtaining average canopy structure information, this reduces manpower and measurement time.

[0052] 4. In image preprocessing, this invention proposes a histogram slope search method. By calculating the slope of the distribution curve in the A band of the LAB color space towards the green direction through the frequency distribution of background values, the minimum value of the background value in the green direction is calculated, and this minimum value is used as a threshold to generate a binary image of the background and vegetation. This solves the problem of automatically or de-classifying image thresholds, and achieves the purpose of automatically distinguishing the background and vegetation in the image and generating a binary image.

[0053] 5. The present invention employs a pixel search method based on pore size and leaf chord length to process binary images, which can conveniently obtain the membership relationship of pixels in the binary image, i.e., whether the pixel belongs to vegetation or background. This allows for the calculation of pore size and leaf chord length in a specific direction, or vegetation width and inter-row distance for ridge crops.

[0054] 6. The MX clustering index equation used in this invention simplifies the problem of choosing the search box size in the Lacunarity-Based clustering index. Mathematically, it makes it easier to calculate the true leaf area and effective leaf area of ​​vegetation. It is not only simple and convenient, but also more practical for image analysis.

[0055] 7. In the calculation of the mean blade tilt angle, compared with the existing polynomial technique that uses empirical algorithms (i.e., statistical method, where the parameters of the equation have no physical meaning; note: the LAI-2000 device uses this method) to calculate the mean blade tilt angle, the physical meaning of each physical quantity in the equation derived in this invention is clearer.

[0056] 8. The equation for the ratio of diffuse radiation to total radiation used in this invention accurately calculates the ratio of diffuse radiation incident on the canopy to the total incident radiation. This ratio is a key influencing factor in calculating the total absorbed photosynthetic radiation ratio (FAPAR). The solution to this problem in this invention helps to calculate FAPAR effectively and correctly. At the same time, the method for calculating the time corresponding to the solar altitude angle facilitates FAPAR photogrammetry.

[0057] 9. By using the method, the cumulative pore size distribution, vegetation width (A1) of row crops, inter-row distance (A2), real leaf area index (LAI), effective leaf area index (LAI_e), average leaf angle (ALA), vegetation canopy aggregation index (CI), vegetation coverage (FVC), crown layer openness (CC) growth season change and absorbed photosynthetic radiation proportion (FAPAR) and other measurement indexes can be conveniently measured.

[0058]

Actual application and verification

[0059] (1) The direct measurement results of the vegetation width (A1) of corn in Zhangye Yinke Oasis, corn in Zhongwei of Ningxia, wheat in Zhongwei of Ningxia, and rice in Zhongwei of Ningxia, and the inter-row distance (A2) of corn in Zhangye Yinke Oasis, corn in Zhongwei of Ningxia, wheat in Zhongwei of Ningxia, and rice in Zhongwei of Ningxia are compared with the measurement method of the present application (see Figure 7 ), and it can be found that the vegetation width (A1) obtained by the direct measurement method (i.e. directly measured on site by using a ruler and other equipment) and the vegetation width (A1) obtained by the measurement method of the present application have high consistency, in which the correlation coefficient (R) is all above 0.86291, and the root mean square error (RMSE) is all below 10.23675 cm. 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 measurement in different regions (i.e. Figure 7 (a)-(d)), there is no obvious rule. Therefore, the measurement method of the present application has high calculation accuracy.

[0060] (2) The real leaf area index (LAI) of corn measured by the direct measurement method is compared with the real leaf area index (LAI) of corn calculated from the image taken under vertical observation by using the present application (see Figure 8 (a)), and it can be found that the root mean square error (RMSE) of the real leaf area index (LAI) measured directly and the real leaf area index (LAI) of the image taken under vertical observation by using the present application is equal to 1.18944, which indicates that the real leaf area index (LAI) measured directly and the real leaf area index (LAI) of the image taken under vertical observation by using the present application have high consistency.

[0061] The effective leaf area index (LAI_e) of corn, wheat and rice measured by the LAI-2000 crown layer analyzer is compared with the effective leaf area index (LAI_e) of corn, wheat and rice calculated from the image taken at an observation zenith angle of 57.5° by using the present application (see Figure 8(b)-(d), it can be found that both have a higher consistency. At the same time, it can be found that the calculation results of wheat are better than that of corn, and the calculation results of corn are better than that of wheat.

[0062] From Figure 8 The results in (a)-(d) can further imply that the key parameters (clumping index) required in calculating the true leaf area index (LAI) and the effective leaf area index (LAI_e) are correct. This also indirectly verifies that the M-X clumping index equation proposed in the present application is correct.

[0063] ⑶ Directly measure the average leaf angle (ALA) of corn in Zhangye Yinke Oasis (see Figure 8 (e) and the average leaf angle (ALA) of corn in Zhongwei, wheat in Zhongwei, and rice in Zhongwei, Ningxia, obtained by LAI-2000 canopy analyzer (see Figure 8 (f)-(h) are compared with the measurement method of the present application, it can be found that there is a high consistency. Among them, the consistency of the average leaf angle calculated by the present application and the average leaf angle obtained by the LAI-2000 canopy analyzer is higher than that of the average leaf angle obtained by the direct measurement method. In the average leaf angle calculation, the measurement of the present application in wheat is better than that in corn, and the measurement in corn is better than that in rice.

[0064] ⑷ Calculate the clumping index (CI) of corn in Zhangye Yinke, corn in Zhongwei, Ningxia, and rice in Zhongwei, Ningxia, by using the method of the present application, as shown in Figure 9 (a), (c) and (e), it can be found that the clumping index (CI) of the vegetation increases with the increase of the Julian day, which is exactly in line with the growth law of the vegetation. However, in Figure 9 (c), there is a decrease after 220 days, which is mainly due to the fact that the corn is in the late stage of maturity, at which time the corn leaves have gradually turned yellow and withered, the vegetation porosity becomes larger, and the clumping degree is weakened, so the clumping index (CI) presents a decreasing trend at this time.

[0065] Calculate the vegetation coverage (FVC) and canopy openness (CC) of corn in Zhangye Yinke, corn in Zhongwei, Ningxia, and rice in Zhongwei, Ningxia, by using the method of the present application, as shown in Figure 9 (b), (d) and (f), it can be found that the vegetation coverage (FVC) and the canopy development degree (CC) present opposite trends. In Figure 9 (d), the vegetation coverage (FVC) presents a decreasing trend after 200 days of Julian day, while the clumping index (CI) presents an increasing trend. This phenomenon and Figure 9(c) The reason described in the above paragraph is the same, and this result implies that the vegetation coverage (FVC) is proportional to the aggregation index (CI), while the canopy development degree (CC) is inversely proportional to the aggregation index (CI). Figure 9 (f), it can be found that the vegetation coverage (FVC) and the canopy development degree (CC) of rice present saturation phenomenon after reaching the Julian day 200.

[0066] ⑸ The direct measurement of the absorbed photosynthetic radiation ratio (FAPAR) of vegetation in Ningxia region by SpectroSense2 multi-channel spectral radiometer and the measurement of the absorbed photosynthetic radiation ratio (FAPAR) by the image of the present application are compared (see Figure 9 (g)), it can be found that the consistency of the two is high, and the distribution points are near the position of the 1:1 line in the figure. Through the research of various vegetation types, it is indirectly proved that the equation of the proportion of diffuse radiation to total radiation proposed in the present application is correct. This also shows that the present application can well measure the absorbed photosynthetic radiation ratio (FAPAR) of vegetation. BRIEF DESCRIPTION OF DRAWINGS

[0067] The specific embodiments of the present application will be further described in detail below in combination with the drawings.

[0068] Figure 1 It is a schematic diagram of the field of view of the camera of the present application and the optical structure of the camera. In the figure: h' is the height of the frame, v' is the width of the frame, d is the diagonal length of the frame, f is the focal length of the camera, 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 position of the camera lens to the top of the canopy, θ o is the zenith angle of the observation direction.

[0069] Figure 2 It is a schematic diagram of the histogram slope search method and the image coordinate system angle calculation of the present application. Among them: (a) is a schematic diagram of the A band pixel distribution frequency in the LAB color space of 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 mathematical abstract graph of the black curve in figure (b), which shows the slope tanα, and the graph 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 the square in the histogram represents the green part value of the A wave band in the LAB color space, i.e. the vegetation; and the black part curve or the square in the histogram represents the green part value of the A wave band in the LAB color space, i.e. the background. In the figure, the symbols (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 included angle between the straight line composed of the points (x1, y1) and (x2, y2) and the axis (y) composed of the pixel groups, and η is the azimuth offset angle in the image.

[0071] Figure 3 The histogram slope search method of the present application is exemplified by taking wheat, corn and rice as examples. In the figures, (a) is the RGB color space image of corn; (b) is the binary image of corn; (c) is the RGB color space image of wheat; (d) is the binary image of wheat; (e) is the RGB color space image of rice; (f) is the binary image of rice; and the symbol FVC is the English abbreviation of vegetation coverage, and T is the threshold value.

[0072] Figure 4 The figure is a schematic diagram of the measurement of the absorbed photosynthetic radiation ratio (FAPAR) of the present application. In the figure, (a) is the angle selected for calculating the absorbed photosynthetic radiation ratio (white sky FAPAR) of the diffuse solar radiation; and (b) is the angle selected for calculating the absorbed photosynthetic radiation ratio (back sky FAPAR) of the average solar direct radiation in a day. The symbols 1-5 represent the selected angles at the time of shooting.

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

[0074] Figure 6 The figure is a schematic diagram of the calculation of the leaf inclination angle based on the leaf chord length of the present application. In the figure, (a) is the three-dimensional situation of the leaf, and (b)-(g) are the geometric relationships of three leaf chord lengths on the plane. In the figure, the symbol is a reference azimuth angle, w is the average leaf chord length in the direction corresponding to the leaf inclination angle, l * is the average leaf chord length in the horizontal direction of the image coordinate, w * is the average leaf chord length in the horizontal direction of the image coordinate, and γ is the included angle between w and l * .

[0075] Figure 7 The verification results of the row crop vegetation width (A1) and the row distance (A2) of the present application. 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 application; (b) is the direct measurement result of the corn row distance (A2) of Zhangye Yingke Oasis and the comparison result of the measurement method of the present application; (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 application; (d) is the direct measurement result of the corn row distance (A2) of Ningxia Zhongwei and the comparison result of the measurement method of the present application; (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 application; (f) is the direct measurement result of the wheat row distance (A2) of Ningxia Zhongwei and the comparison result of the measurement method of the present application; (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 application; (h) is the direct measurement result of the rice row distance (A2) of Ningxia Zhongwei and the comparison result of the measurement method of the present application. In the figure, the D behind the horizontal line represents direct measurement, and the P represents image measurement.

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

[0077] Figure 9Validation of the vegetation canopy index (CI), fractional vegetation cover (FVC), canopy openness (CC) and fraction of absorbed photosynthetically active radiation (FAPAR) for the present invention. (a) is the calculated CI of the corn in Zhangye, Gansu province, (b) is the calculated FVC and CC of the corn in Zhangye, Gansu province, (c) is the calculated CI of the corn in Zhongwei, Ningxia province, (d) is the calculated FVC and CC of the corn in Zhongwei, Ningxia province, (e) is the calculated CI of the rice in Zhongwei, Ningxia province, (f) is the calculated FVC and CC of the rice in Zhongwei, Ningxia province, (g) is the comparison of the FAPAR measured by the SpectroSense2 multi-channel spectroradiometer and the FAPAR calculated by the present invention in the Ningxia province. The symbol D after the horizontal line represents direct measurement, and the symbol P represents image measurement. DETAILED DESCRIPTION

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

[0079] (1) Determine the size of the CCD or COMS frame 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 size of the CCD or COMS frame of the non-fisheye digital camera 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 manufacturer, as shown in Table 1.

[0081] Table 1 Camera reference frame

【Making (some) sense out of sensor sizes [EB / OL].

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

[0083] Frame name Frame width (mm) Frame 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, the field of view FOV h , FOV v , and FOV d of the non-fisheye digital camera are calculated using equations (1-3):

[0085]

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

[0087] The actual length and width corresponding to the captured vegetation image are further calculated.

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

[0089]

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

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

[0092] If the relevant information of the camera is unknown, equations (6-7) are used for calculation:

[0093]

[0094] W s = L s x 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 library of the default digital camera (see Table 2). Among them: the calculation accuracy of equations (6-7) is not as high as that based on equations (4-5).

[0096] Table 2 Relationship between focal length and angle of view 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) Angular eye 180° 14m 114° 20 mm 94° 24 mm 84° 25m 75° 35m 63° 50 mm 46° 70 mm 34° 80 mm 30° 85m 28°30′ 100 mm 24° 135 mm 18° 200 mm 12° 300 mm 8°15′ 400 mm 6°10′ 500 mm 5° 600 mm 4°10′ 800 nm 3°5′

[0098] In the measurement, in order to make the camera shoot the image of the appropriate height, the non-fisheye digital camera can also be loaded on the telescopic rod or the unmanned aerial vehicle.

[0099] (2) According to the selection of vegetation parameters, the shooting angle and shooting time of the image are selected to shoot and sample, and the image is obtained.

[0100] The conditions of shooting sampling are:

[0101] When the vertical observation is made, the vegetation parameters that can be measured are the cumulative pore size distribution, the vegetation width of the ridge row crop, the distance between ridges, the vegetation coverage, the canopy openness, the slightly error real leaf area index, the slightly error effective leaf area index, and the average leaf inclination.

[0102] When the observation zenith angle of 57.5° is used for measurement, the vegetation parameters that can be measured are the accurate real leaf area index and the 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

[0103] When measuring the absorption of photosynthetic radiation of direct solar radiation (i.e., back sky FAPAR), the shooting direction of the sun incident zenith angle needs to be selected.

[0104] When measuring the absorption of photosynthetic radiation of sky diffuse radiation, the average sampling within 180 degrees of the zenith angle is required.

[0105] When measuring the average proportion of absorbed photosynthetic radiation from direct solar radiation throughout the day, it is necessary to use a 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 calculation program.

[0107] When calculating the average absorbed photosynthetic radiation ratio of direct solar radiation throughout the day, time-series FAPAR data is required, and the time calculation program corresponding to the solar altitude angle can calculate the time corresponding to the measurement for acquiring time-series FAPAR data. First, the solar altitude angle throughout the day is calculated using the following equation:

[0108]

[0109]

[0110] In the formula: θ ⊙_noon The solar noon altitude angle; δ is the latitude measured locally; δ is the solar declination; d n For Julian Japan; ω sunrise ω is the local sunrise angle. sunset Let be the local sunset angle. Calculate the solar altitude angle at noon using equation (8), and since the solar altitude angle at sunset and sunrise is 0°, then we can calculate... Figure 4 (b) The solar altitude in steps 1-5. Simultaneously, based on equations (9-11), and then combined with equations (12-13), the sunrise time (T) can be calculated. sunrise ) and sunset time (T sunset ).

[0111]

[0112] Given that the solar altitude angle in the main text is at noon, and the times of sunset and sunrise are known, it is possible to calculate... Figure 4 (b) The time corresponding to the solar altitude of 1-5. Equation (8-13) serves Equation (34), providing it with the average measurement time and measurement altitude angle during the day.

[0113] (3) Perform Gamma transformation on the RGB color space of the image to correct color grayscale differences caused by incorrect exposure during image capture. Then convert the RGB color space image to the LAB color space. The specific process is as follows:

[0114] First, the image is processed using the Gamma function to perform simple grayscale processing to correct the abnormal exposure. The Gamma function is given by equation (14-17).

[0115]

[0116] where r is the red band of the camera, unitless; R is the red band after Gamma transformation, unitless; g is the green band of the camera, unitless; G is the green band after Gamma transformation, unitless; b is the blue band of the camera, unitless; B is the red band after Gamma transformation, unitless; x is a variable which is the value in the bracket of equations (14-16).

[0117] Then, the image in the RGB color space after Gamma transformation is converted into the XYZ color space image by equation (18).

[0118]

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

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

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

[0122] B = 200 x [f(Y / Y n ) - f(Z / Z n )] (21),

[0123]

[0124] where X n is equal to 95.047, Y n is equal to 100.0, and Z n is equal to 108.883 (these are the normalized values under the D 65 light source). L is the luminance from black to white, unitless; A is the luminance from green to red, unitless; B is the luminance from blue to yellow, unitless; t is a variable which is the value in the small bracket of equations (19-21), unitless; f(t) is a variable 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, preprocessing the A band in the LAB color space can produce binary images of the background and vegetation.

[0126] 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 axis (x) and the number of pixels corresponding to that value as the vertical axis (y);

[0128] ② Calculate the derivative f'[N(x)] of the frequency function f[N] composed of the number of pixels (y) of the A-band values ​​in the LAB color space, and calculate the extreme value in the frequency histogram so that the extreme value satisfies the requirement of {f'[N(x)]>0}∧{f'[N(x)]>0}. Then, take the second local minimum value of the coordinate axis (x) of the pixel value in the positive direction as the intersection point of the Gaussian distribution in the frequency of vegetation value and the Gaussian distribution in the frequency of background value. The area near this intersection point is the region where the threshold 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 to search for the minimum value of the Gaussian distribution of the frequency of the background value towards 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, and obtain binary images of the background and vegetation; where 255 in the binary image corresponds to vegetation and 0 corresponds to the background.

[0131] The specific process is as follows:

[0132] Because in the A-band, the frequency of the number of pixels corresponding to soil and vegetation values ​​follows a double Gaussian distribution (i.e., Figure 2 (a)

Yaokai Liu XM, Haoxing Wang & Guangjian Yan. A novel method for extracting green fractional vegetation cover from digital images[J]. Journal of Vegetation Science, 2012, 23(3): 406-418

[0133]

[0134] In the formula: 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 double Gaussian distribution curves of soil and vegetation is... Figure 2 (a) is the point (x6, y6).

[0136] Since the frequency of the primitive values ​​in the histogram is a discrete function rather than a continuous function, local minima cannot be obtained. However, an approximate point (x2, y2) can be obtained, and then... Figure 2 (c) Calculate the slope tanα of the line corresponding to the points (x1, y1) and (x2, y2) in the coordinates of these two points, and then calculate the point (x4, y4). Again, using the points (x2, y2) and (x4, y4), calculate the midpoint (x3, y3) of the line between these two points. Finally, using the coordinate equation of the previously calculated tanα, consider an offset slope tanβ, and together calculate the minimum value of the predicted Gaussian distribution curve of the soil towards the green direction, which is the point (x5, y5).

[0137] The above calculation steps can be converted into the following mathematical equation:

[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 value; D, E, F are functions in the threshold value T equation; (x1, y1) and (x2, y2) are points near the focus in the frequency histogram; β is the offset angle; α is the angle between the straight line composed of (x1, y1) and (x2, y2) and the axis (y) composed of the pixel groups; τ is the control coefficient of tan α, and the default value is 0.7.

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

[0145] (5) The proportions of the background and the vegetation in the binary image are counted, and the vegetation coverage (FVC), the crown cover (CC), the porosity, the absorption photosynthetic radiation ratio of the direct solar radiation (i.e., back sky FAPAR), and the absorption photosynthetic radiation ratio of the sky diffuse radiation (i.e., white sky FAPAR) are calculated. At the same time, according to the absorption photosynthetic radiation ratio corresponding to the shooting time of the direct solar radiation, the average absorption photosynthetic radiation ratio of the direct solar radiation in a day is calculated. For the total absorption photosynthetic radiation ratio calculation, the ratio of the incident diffuse radiation to the total radiation is calculated, and the absorption photosynthetic radiation ratio of the direct solar radiation and the absorption photosynthetic radiation ratio of the sky diffuse radiation are integrated to obtain the instantaneous total absorption photosynthetic radiation ratio (i.e., FAPAR).

[0146] Where: If the photograph is taken at a vertical angle, the vegetation coverage (FVC) and the crown cover (CC) in the image can be calculated. The calculation equations are 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, and the two parameters are obtained by searching the porosity and the leaf chord length pixel search method mentioned below.

[0149] The variation trend of equations (29-30) is shown in Figure 9 (b) and Figure 9 (d).

[0150] If the photograph is taken at a specific angle, the porosity Po (θ o ), i.e. equation (31).

[0151]

[0152] In FAPAR calculation, equation (32) can calculate the instantaneous back sky FAPAR (i.e. the absorbed photosynthetically active radiation ratio of direct solar radiation, its symbol is fAPAR bs ), while equation (32) can calculate the instantaneous white sky FAPAR (i.e. the absorbed photosynthetically active radiation ratio of direct solar radiation of sky diffuse radiation, fAPAR ws )

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

[0153]

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

[0155] Equations (33-34) are integrals, and five photos of different angles are used in the present application to approximate this integral. In equation (33), the angle is selected as Figure 4 (a). And the angle selection of equation (34) is Figure 4 (b).

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

[0157] Consider the total FAPAR of direct solar radiation and sky diffuse radiation as

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

[0159]

[0160] where sky is the proportion of diffuse radiation at the top of the canopy to the total radiation, which is derived from the equation in the Gap Light Analyzer (GLA) manual.

[0161] R p / R s represents the ratio of photosynthetically active radiation and shortwave radiation. According to the research, it can be obtained by table lookup or directly measured, and details can be found in 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

[0162] (6) According to the binary image, the pore and leaf chord length pixel search method is used to calculate the vertical direction of the image leaf chord length and pore size, the horizontal direction of the leaf chord length and pore size, and the leaf chord length and pore size in a certain direction. Further calculation can obtain the cumulative pore size distribution of each pixel direction; at the same time, according to the pore and leaf chord length pixel search method, the vegetation width and inter-row distance of ridge row crops can be calculated.

[0163] Wherein: 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 membership relationship of the image pixel, that is, whether the pixel belongs to vegetation or background. Its 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 size of the aperture length; x and y are the pixel coordinate positions in the image coordinate system (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; and P is the pixel value. The search condition is:

[0168]

[0169] Figure 2 (d) gives a schematic diagram of the angle search. Taking the horizontal direction as the starting point, the image coordinate values in the angle search are calculated using the following equations:

[0170]

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

[0172] In the equations, η is the azimuthal offset angle in the image, i.e., Figure 2 (d). Note that the image starting point (i.e., point (x1, y1) in Figure 5 (a)) and the ending point (i.e., point (x2, y2) in Figure 5 (a)) in the calculation of the blade chord length and the aperture size are also calculated using equations (39-40).

[0173] In the equations, η is the azimuthal offset angle in the image, i.e., Figure 5 (a) and (b) only show a schematic diagram, and the measurer can select the encryption for the truncation according to his / her own requirements. In the diagram, (x1, y1) and (x2, y2) are the starting point and the ending point for the aperture size calculation in the image coordinate system (here, only the aperture size coordinate acquisition is taken as an example), (x s , y s ) is the starting point pixel coordinate in the row, column or specific azimuth in the image coordinate system, and (x e , y e ) is the ending point pixel coordinate in the row, column or specific azimuth in the image coordinate system.

[0174] Equations (37-40) are a calculation example of the aperture size and the blade chord length. Based on the same mathematical principle, the vegetation width Al and the inter-ridge distance A2 of ridge-row crops can be calculated (note that ridge-row crops refer to the planting form of farmland crops, which is divided by the soil during planting. The interval of the division is called the inter-ridge distance, and the width of the divided vegetation part is called the vegetation width, see Figure 5 (b)). The verification results can be seen in Figure 7 .

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

[0176]

[0177] In the formula, 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 is the number of pores calculated by the image element in a certain direction of the image; λ max is the maximum pore size in a certain direction; (x s ,y s ) is the starting point pixel coordinate of the row, column or specific direction in the image coordinate system, (x e ,y e ) is the end point pixel coordinate of the row, column or specific direction in the image coordinate system, and the calculation method can be referred to Figure 5 .

[0178] (2) According to the cumulative pore size distribution of each pixel direction, the aggregation index of vegetation is calculated by using the M-X aggregation index equation; and then the true leaf area index with slight error and the effective leaf area index with slight error in the image vertically shot are calculated according to the aggregation index, and the accurate true leaf area index and effective leaf area index in the image shot at an observation zenith angle of 57.5° are calculated.

[0179] According to the previous research on the heterogeneity of variance and expectation

Dale M R T. Lacunarity analysis of spatial pattern: A comparison [J]. Landscape Ecology, 15 (5): 467-478.

[0180]

[0181] In the formula, Ω E is the M-X aggregation index; Λ (λ) is the heterogeneity index between all image elements in a single row or column; Λ max (λ) is the maximum value of Λ (λ) in the 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 image elements in a single row or column presents a random distribution form, and Λ cri (λ) = 1; σ 2(λ) represents the statistical variance among all pixels in a single row of the image; E(λ) represents the statistical expectation among all pixels in a single row of the image; λ i The aperture size calculated for a single row or column of pixels in the i-th image; F(λ) represents the statistical average value of the aperture size calculated from pixels in a single row or column of the image, within that single row; N represents the number of apertures calculated from pixels in a single row or column of the image; F(λ) i ) is λ i The corresponding cumulative distribution function value of pore size.

[0182] In equation (42), it is analogous to the standardization method of the Lacunarity-Based clustering index, where A max (λ)=L image A cri (λ) = 1, and their determination was based on the literature [Fournier RA, Hall R J. Hemispherical Photography in Forest Science: Theory, Methods, Applications [M]. Van Goghewicz: Springer, 2017.], and their variation trend is as follows. Figure 9 (a) and Figure 9 (c)).

[0183] Based on the aggregation index, the true leaf area index and the effective leaf area index can be calculated. According to the research of F. Baret [Baret F, De Solan B, Lopez-Lozano R, et al. GAI estimates of row crops from downward looking digital photos taken perpendicular to row crops 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.], when the observed zenith angle is 57.5°, the G function is equal to a constant of 0.5, and the osmotic function for calculating the leaf area index can be accurately calculated.

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

[0185]

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

[0187] P = (1 - P0) / (1 - P0) o P0(57.5°) is the porosity when the observation zenith angle is 57.5°, which is the digital camera measurement value when the observation zenith angle is 57.5°, and can be calculated by equation (31).

[0188] Equation LAI e The verification results are shown in Figs. Figure 8 (b) to (d). According to the research of J. Liu, when the observation zenith angle is 0° (i.e., vertical observation), the extinction coefficient at this time is approximately 0.61

Liu J, Pattey E, Admiral S. Assessment of in situ crop LAI measurement using unidirectional view digital photography [J]. Agricultural and Forest Meteorology, 2013, 169: 25-34

[0189]

[0190] P = (1 - P0) / (1 - P0) o P0(0°) is the porosity when the observation zenith angle is 0°, which is the digital camera measurement value when the observation zenith angle is 0° (i.e., vertical observation), and can be calculated by equation (31).

[0191] At this time, equation (46) is substituted into equation (45) to calculate the effective leaf area index of the vertical observation. It is noted that the true leaf area index calculated by equation (46) has a small error, and therefore, in the present application, it is referred to as a slightly erroneous true leaf area index, and the effective leaf area index calculated thereby is referred to as a slightly erroneous effective leaf area index. Although the calculation result has an error, it basically meets the measurement requirements, and the verification results are shown in Fig. Figure 8 (a).

[0192] ⑻According to the leaf chord length in the vertical direction, the leaf chord length in the horizontal direction, and the leaf chord length in a specific direction, the average leaf inclination angle of the vegetation in the image can be calculated by using the average leaf inclination angle equation.

[0193] In the study of vegetation "hot spot" factor, an equation for calculating leaf chord length was proposed

Ma X, Lu L. Mathematical analysis on the component of canopy architecture in row crops. 2019

[0194]

[0195] where w is the leaf chord length; P ol is the leaf distribution form; b is the leaf bending factor; l * is the leaf length; w * is the leaf width; θ l is the average leaf inclination; and φ is the observation azimuth.

[0196] If it is vertical observation, P ol is the leaf distribution form on the stem, which is similar to the aggregation index Ω E Therefore, let P ol = Ω E , the average leaf inclination equation is:

[0197]

[0198] where is the transformation, and all this is done because the accurate leaf length and leaf width are difficult to obtain in the picture. Here is a reference azimuth, and the average chord length of the leaf corresponding to this direction is w, while the average chord length of the leaf searched in the horizontal direction of the image coordinate (i.e., the length direction) is l * , and the average chord length of the leaf searched in the vertical direction of the image coordinate (i.e., the width direction) is w * . γ is the angle between w and l * . Various cases of equation (48) are shown in Figure 6 (b)~(g), and the verification results are shown in Figure 8 (e)~(h).

Claims

1. A method for measuring the structure of vegetation canopy based on a non-fisheye digital camera, comprising the following steps: (1) determining the size of the CCD or COMS frame of the non-fisheye digital camera and the vertical distance from the shooting position of the camera lens to the top of the canopy, calculating the field of view of the non-fisheye digital camera according to the corresponding mathematical equation, and then calculating the actual length and width corresponding to the shot vegetation image; (2) shooting and sampling according to the shooting angle and shooting time of the image selected according to the measured vegetation parameters to obtain the image; the shooting time refers to the time corresponding to a certain solar elevation angle calculated by the solar elevation angle corresponding time calculation program; (3) performing Gamma transformation on the RGB color space of the image, and then converting the changed RGB color space image into LAB color space; (4) preprocessing the A band in the LAB color space to obtain the binary image of the background and the vegetation; (5) calculating the vegetation coverage, canopy openness, porosity, the proportion of absorbed photosynthetic radiation of solar direct radiation, and the proportion of absorbed photosynthetic radiation of sky diffuse radiation by counting the proportion of the two values of background and vegetation in the binary image; At the same time, according to the shooting time corresponding to the proportion of absorbed photosynthetic radiation of solar direct radiation, the average proportion of absorbed photosynthetic radiation of solar direct radiation in a day is calculated, and the proportion of diffuse radiation in total radiation is calculated according to the equation of the proportion of diffuse radiation in total radiation, and then the instantaneous FAPAR is calculated; (6) according to the binary image, the pore and leaf chord length pixel search method is used to calculate the vertical and horizontal leaf chord length and pore length, and the leaf chord length and pore length in a specific direction; and then the cumulative pore size distribution in each pixel direction is calculated; at the same time, the vegetation width and inter-row distance of row 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 equation is used to calculate the aggregation index of the vegetation; and then the true leaf area index and effective leaf area index with slight error in the image shot in the vertical direction, and the accurate true leaf area index and effective leaf area index in the image shot at an observation zenith angle of 57.5° are calculated according to the aggregation index; the M-X aggregation index equation refers to an index for describing the aggregation degree of leaves in the image, and its equation form is: ， ， where Ω E is the M-X aggregation index; Λ (8) according to the vertical leaf chord length, horizontal leaf chord length, and leaf chord length in a specific direction, the average leaf inclination equation can be used to calculate the average leaf inclination of the vegetation in the image. is the heterogeneity index between all pixels in a single row or column of the image; Λ max In the step (1), the non-fisheye digital camera is loaded on the telescopic rod or the unmanned aerial vehicle. is the maximum value of Λ ​ when calculated by column or row, and Λ max ​ = 1 L image , L image is the number of image lengths; Λ cri ​ is the value of the heterogeneity index between all pixels in a single row or column of the image when presented in a random distribution, 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 is the pore size calculated from the pixels in a single row or column of the ith image; is the statistical equilibrium value of the pore size calculated from the pixels in a single row or column of the image in a single row; N t is the number of pores calculated from the pixels in a single row or column of the image; F ​ i is the cumulative distribution function value of the pore size corresponding to ​ i .​​​​​​ ​ 2. The method for measuring vegetation canopy structure based on a non-fisheye digital camera as described in claim 1, characterized in that: ​ 3. The method for measuring vegetation canopy structure based on a non-fisheye digital camera as described in claim 1, characterized in that: The sampling condition in step 2 refers to: when vertical shooting is adopted, the measured vegetation parameters are cumulative pore size distribution, vegetation width of ridge row crops, distance between ridges, vegetation coverage, canopy openness, slightly inaccurate true leaf area index, slightly inaccurate effective leaf area index, and average leaf inclination angle; when the observation zenith angle is 57.5°, the measured vegetation parameters are accurate true leaf area index and accurate effective leaf area index; when the absorption photosynthetic radiation ratio of solar direct radiation is measured, the shooting is performed at the zenith angle of the solar incident direction; when the absorption photosynthetic radiation ratio of sky diffuse radiation is measured, the average sampling within 180 degrees of the zenith angle is required; when the average absorption photosynthetic radiation ratio of solar direct radiation in a day is measured, the shooting time is determined by using the time calculation program corresponding to the solar elevation angle.

4. The method for measuring vegetation canopy structure based on a non-fisheye digital camera as described in claim 1, characterized in that: The solar elevation angle corresponding time calculation program in step 2 refers to a combination of equations for calculating the solar elevation angle and the time corresponding to the solar elevation angle, and the equations are as follows: ， ， ， ， ， ， where: θ ⊙_noon is the solar noon altitude angle; φ l is the locally measured 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 is the local sunset time.

5. The method of claim 1, wherein: The pre-processing in step 7 includes the following steps: ​ x ), the number of pixels corresponding to the value ( y ) as the vertical coordinate of the frequency histogram;​ ii) the frequency function of the number of pixels in the LAB color space having a value of the A band y ) of the A band f [ N ] of the A band f’ [ N (x)]: and the extreme value in the frequency histogram is calculated, such that the extreme value satisfies the requirements of { f’ [ N (x)]>0} and { f’ [ N (x)]>0}. wherein: f’ [ N (x)] is the first discrete value of the derivative of the A-band soil and vegetation dual Gaussian distribution curve to the positive direction of the minimum value; f [ N (x)] is the first discrete value of the derivative of the A-band soil and vegetation dual Gaussian distribution curve to the positive direction of the minimum value; f [ N (x + )] is the second discrete value of the derivative of the A-band soil and vegetation dual Gaussian distribution curve to the positive direction of the minimum value; Then, using the coordinate axis of cell values ​​( x The second local minimum value in the positive direction of the vegetation value is taken 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 this intersection is the region where the threshold to be determined is located. ③ the point of intersection and the point near the point of intersection (x, y) and (x, y) x 1 , y 1 ) as the basis, using the histogram slope search equation to search for the minimum value of the Gaussian distribution in the frequency of the background value towards the green direction, taking the minimum value as the threshold value; x 2 , y 2 ) as the basis, using the histogram slope search equation to search for the minimum value of the Gaussian distribution in the frequency of the background value towards the green direction, taking the minimum value as the threshold value; The histogram slope search equation is as follows: where T is a threshold value; D , E , F are threshold values T are functions in the equation; points ( x 1 , y 1 ) and points ( x 2 , y 2 ) are points near the focus in the frequency histogram; β is an offset angle; τ is a control coefficient of tan α, and α is an included angle between a straight line composed of points ( x 1 , y 1 ) and points ( x 2 , y 2 ) and an axis composed of a plurality of pixel groups ( y ).

4. Classify the A wave band vegetation and background in the LAB color space according to the threshold value to obtain binary images of the background and the vegetation; wherein 255 in the binary image corresponds to the vegetation, and 0 corresponds to the background.

6. The method for measuring vegetation canopy structure based on a non-fisheye digital camera as described in claim 1, characterized in that: The diffuse radiation proportion of total radiation equation in step 8 refers to the proportion of sky diffuse radiation in total radiation, and the equation is calculated as follows: ， wherein: sky is the ratio of diffuse radiation incident at the top of the canopy and total incident radiation; R p / R s represents the ratio of photosynthetically active radiation and shortwave radiation.

7. The method of claim 1, wherein the non-fisheye digital camera is based on a camera with a lens having a focal length of 8 mm and a field of view of 180 degrees. 8 The pore and leaf chord length pixel search method in step 9 refers to an algorithm for searching the membership relationship of pixels in the image; and the algorithm expression is as follows: ， wherein = is the mathematical symbol "defined as"; w is the blade chord length; λ is a pore size length; x and y is a pixel coordinate position in image coordinates; L image is a number of image lengths; W image is a number of image widths; P is a pixel value; The search condition is as follows: ， Taking the horizontal direction as the starting point, the image coordinate value in the angle search is calculated using the following equation: , , wherein: η is the offset angle of the orientation in the image.

8. The method for measuring vegetation canopy structure based on a non-fisheye digital camera as described in claim 1, characterized in that: The average leaf inclination angle equation in step 11 refers to calculating the average value of the leaf inclination angle of the vegetation canopy, and the equation is as follows: ， ， wherein: θ l is the average leaf angle; b is the bending factor of the leaf. φ is a reference orientation; w L is the average chord length of the blade corresponding to the reference orientation; l * the average chord length of the leaflets searched in the horizontal direction of the image coordinates; w * the average chord length of the leaflets searched in the vertical direction of the image coordinates; γ w and l * the angle between the two vectors; Ω E the M-X aggregation index.​

Images