JP2008111725A - Leaf area index calculating device, leaf area index calculating method and program thereof - Google Patents

Abstract

There is provided a leaf area index calculation method capable of estimating a leaf area index with high accuracy using aeronautical laser measurement.
A method for storing laser reflection intensity data reflected from the ground by irradiating a laser from the sky and reflecting on the ground, a procedure for storing the reflectance of the ground, and a reflection from the stored reflection intensity data by a tree crown. Of extracting the reflected intensity of the laser (S13), extracting the reflected intensity of the laser reflected by the ground after passing through the canopy from the stored reflected intensity data, and the stored reflected intensity data and reflection Based on the rate, the procedure (S12) for calculating the incident intensity of the laser immediately before reaching the crown, the incident intensity of the laser immediately before reaching the crown, the reflected intensity of the laser passing through the crown, and the laser reflected by the crown And calculating the leaf area index of the tree using the reflection intensity (S15).
[Selection] Figure 8

Description

  The present invention relates to a leaf area index calculating apparatus, a leaf area index calculating method, and a program thereof for calculating an important leaf area index as a stand parameter.

  Understanding forest parameters such as the number of trees, crown area, tree height, breast height diameter, and leaf area index (LAI) is important for quantifying forest conditions and reflecting them in forestry and conservation activities. It is. Here, LAI is an index indicating a certain degree of the leaf of a plant, and is a ratio of the total area of one side of all the leaves existing above the unit area on the ground (for example, patent document) 1).

  As an LAI estimation method, a field survey is generally performed. However, this method can only estimate locally, and it takes time and labor to estimate a wide area LAI.

  Therefore, there is a method of utilizing satellite data when estimating LAI over a wide area. This is based on a regression equation obtained by obtaining a normalized vegetation index (NDVI) from satellite data and performing regression analysis of NDVI and LAI. This method is also based on the premise that an actual measurement value of LAI is acquired and can be expanded over a wide area, but cannot be implemented without actual measurement data.

There is also a method for estimating a wide-area LAI using an aviation laser. In this method, the LAI is estimated using the height (tree height data) from the ground surface obtained by taking the difference between the height of the laser reflection point measured by the aviation laser and the DEM. Specifically, a cumulative frequency distribution is created from a histogram with the horizontal axis as the height from the ground surface and the vertical axis as the frequency, and the LAI is calculated on the assumption that this represents light attenuation.
JP 2005-52045 A

  However, in the conventional LAI estimation method using an aerial laser, since the histogram does not reflect the information inside the tree crown, the calculated LAI remains questionable. That is, there is a problem that LAI cannot be estimated with high accuracy.

  In view of the above problems, the present invention is to provide a leaf area index calculation apparatus, a leaf area index calculation method, and a program thereof that can estimate LAI in a wide area and easily with high accuracy using aviation laser measurement. .

  According to one aspect of the present invention, (b) irradiating a laser from the sky to the ground and storing reflection intensity data of the laser reflected on the ground; (b) storing the reflectance of the ground; C) A step of extracting the reflection intensity of the laser reflected by the tree crown from the stored reflection intensity data, and (d) extracting the reflection intensity of the laser reflected by the ground after passing through the tree crown from the stored reflection intensity data. And (e) calculating the incident intensity of the laser immediately before reaching the crown based on the stored reflection intensity data and reflectance, and (f) a laser reflected on the ground from the stored reflection intensity data. The step of extracting the reflection intensity of the laser, and (g) the incident intensity of the laser just before reaching the crown, the reflection intensity of the laser that passed through the crown, and the reflection of the laser reflected by the crown Degrees, and by using the reflection intensity of the laser reflected by the ground, leaf area index calculating method comprising the steps of calculating the leaf area index tree is provided.

  According to another aspect of the present invention, the computer stores (b) means for irradiating a laser beam from the sky to the ground and storing the reflection intensity data of the laser beam reflected on the ground, and (b) storing the reflectance of the ground. Means, (c) means for extracting the reflection intensity of the laser reflected by the tree crown from the stored reflection intensity data, and (d) the laser reflected from the ground after passing through the tree crown from the stored reflection intensity data. Means for extracting the reflection intensity; (e) means for calculating the incident intensity of the laser immediately before reaching the tree crown based on the stored reflection intensity data and reflectance; and (f) the ground from the stored reflection intensity data. The procedure to extract the reflection intensity of the laser reflected in (1), (g) the incident intensity of the laser just before reaching the crown, the reflection intensity of the laser that passed through the crown, the reflection intensity of the laser reflected by the crown, Using reflection intensity of the laser reflected by the fine ground, a program to function as the means means for calculating the leaf area index tree is provided.

  According to still another aspect of the present invention, (a) a laser irradiation means for irradiating a laser from the sky to the ground, (b) a reflection intensity storage means for storing reflection intensity data of the laser reflected on the ground, and (c) ) Reflectivity storage means for storing the reflectance of the ground, and (d) Extracting the reflection intensity of the laser reflected by the tree crown from the stored reflection intensity data, and passing the tree crown from the stored reflection intensity data Extracts the reflection intensity of the laser reflected on the ground, extracts the reflection intensity of the laser reflected on the ground from the stored reflection intensity data, and based on the stored reflection intensity data and reflectance, the laser immediately before reaching the tree crown The incident intensity of the laser just before reaching the crown, the reflected intensity of the laser that passed through the crown, the reflected intensity of the laser reflected by the crown, and the level reflected by the ground Using reflection intensity of Heather, a calculating means for calculating a leaf area index tree, leaf area calculation unit and an output means for outputting the (e) leaf area index is provided.

  According to the present invention, it is possible to provide a leaf area index calculation device, a leaf area index calculation method, and a program thereof that can estimate LAI with high accuracy using aeronautical laser measurement.

Therefore, it is possible to accurately estimate the CO ₂ reduction amount using this LAI as a stand parameter.

  Next, first and second embodiments of the present invention will be described with reference to the drawings. In the following description of the drawings, the same or similar parts are denoted by the same or similar reference numerals. Further, the embodiments described below exemplify apparatuses and methods for embodying the technical idea of the present invention, and the technical idea of the present invention includes the material, shape, structure, The layout is not specified as follows. The technical idea of the present invention can be variously modified within the scope of the claims.

  As shown in FIG. 1, the leaf area index calculation device according to the first embodiment of the present invention includes a central processing unit (CPU) 1, a display unit 2, an input unit 3, an output unit 4, and a storage device 5. . As the input unit 3, for example, a recognition device such as a keyboard, a mouse, and an OCR, an image input device such as a scanner and a camera, a voice input device such as a microphone, and the like can be used. As the output unit 4, a display device such as a liquid crystal display or a CRT display, a printing device such as an ink jet printer or a laser printer, or the like can be used. The display unit 2 is a liquid crystal display, a CRT display, or the like. The storage device 5 includes a ROM and a RAM. The ROM stores a program executed by the CPU 1. The RAM functions as a temporary data memory or the like that temporarily stores data or the like used during program execution processing in the CPU 1 or is used as a work area. As the storage device 5, for example, a semiconductor memory, a magnetic disk, an optical disk, a magneto-optical disk, a magnetic tape, or the like can be used.

  FIG. 2 is a schematic configuration diagram of the leaf area index calculating apparatus according to the first embodiment of this invention. (First) ground reflection intensity extraction unit 6, incident intensity calculation unit 7, display processing unit 10, (second) ground reflection intensity extraction unit 11, reflectance calculation unit 12, canopy reflection intensity extraction unit 13, canopy passage The intensity extraction unit 14 and the leaf area index calculation unit 15 are realized by the CPU 1 shown in FIG.

  The ground reflection intensity storage unit 8, the reflectance storage unit 9, the orthophoto image storage unit 20, the reflection intensity data storage unit 21, the reflectance storage unit 22, the tree vertex storage unit 23, the DEM storage unit 24, and the threshold storage shown in FIG. The unit 25, the ground reflection intensity storage unit 26, the incident intensity storage unit 27, the canopy reflection intensity storage unit 28, the canopy passage intensity storage unit 29, and the leaf area index storage unit 31 are included in the storage device 5 shown in FIG.

  The orthophoto image storage unit 20 stores an orthophoto image including the forest area A that is an LAI calculation target. An “orthophoto image” is an image obtained by orthogonal transformation by giving three-dimensional measurement data or the like to a photographic image. The orthophoto image includes, for example, coordinate (X, Y) and corresponding color information of red (R), green (G), and blue (B).

  The reflection intensity data storage unit 21 stores reflection intensity data of an aviation laser (hereinafter simply referred to as “laser”). The reflection intensity data can be obtained by irradiating the ground with laser by means of laser irradiation means (not shown) provided in the aircraft and receiving the laser reflected on the ground by laser receiving means (not shown) in the aircraft. . The reflection intensity data is random data obtained at an arbitrary point, and includes three-dimensional coordinates (X, Y, Zd) and the reflection intensity of the laser. The height coordinate Zd indicates the height of a laser reflection point such as the ground or a tree.

The reflection intensity of the laser is indicated by a pulse waveform. The pulse waveform has a plurality of peak values. The first peak value is defined as “first pulse intensity I _first ”, and the last peak value is defined as “last pulse intensity I _last ”.

  The DEM storage unit 24 stores grid data (DEM: Digital Elevation Model) of forest ground. The DEM, for example, is assigned to mesh lattice points obtained by dividing the X-axis and Y-axis of a plane rectangular coordinate system at a constant interval (0.5 m). The DEM includes three-dimensional coordinates (X, Y, Zh). The height coordinate Zh indicates the height of the ground.

The ground reflection intensity extraction unit 6 reads the first pulse intensity I _first included in the plurality of reflection intensity data of the coordinates (X, Y) in the known reflectance region from the reflection intensity data storage unit 21. Here, the “known reflectance region” refers to a region where the reflectance of the ground is known in advance, for example, an asphalt or grassland region. Further, the ground reflection intensity extraction unit 6 calculates the average value of the first pulse intensities I _first , and the average value of the reflected laser reaches the ground in the known reflectance region as shown in FIG. The reflection intensity (hereinafter referred to as “ground reflection intensity”) Ia is extracted. The ground reflection intensity Ia in the known reflectance region is stored in the ground reflection intensity storage unit 8.

  The incident intensity calculation unit 7 reads the ground reflection intensity Ia in the known reflectance region from the ground reflection intensity storage unit 8. Further, the incident intensity calculation unit 7 reads a known reflectance Ra corresponding to the reflection intensity Ia from the reflectance storage unit 9. The known reflectance Ra can be set from a general spectral reflection characteristic graph as shown in FIG. The incident intensity calculation unit 7 uses the ground reflection intensity Ia and the known reflectance Ra in the known reflectance region, and the intensity of the laser incident on the ground as shown in FIG. , Referred to as “incident intensity.”) I0 is calculated.

I0 = Ia / Ra (1)
For example, assuming that the value of the ground reflection intensity Ia is 12.2 and the wavelength of the laser used is about 1.064 [μm], the value of the reflectance Ra is 0.14. The value is 87.12. The calculated incident intensity I0 is stored in the incident intensity storage unit 27.

  The tree vertex storage unit 23 stores the coordinates of tree vertices for each tree in the forest area A (a method for calculating the tree vertices will be described later). In the first embodiment of the present invention, as an example, LAI of trees (crowns) T1 to T6 in the forest area A is calculated in the orthophoto image shown in FIG. For example, the orthophoto image storage unit 20 may store the shapes and coordinates of trees (tree crowns) T1 to T6 in association with the orthophoto image. Note that the calculation target of the LAI is not limited to the trees T1 to T6. For example, the LAI of all the trees in the forest area A may be calculated. The tree vertex storage unit 23 stores the coordinates (Xi, Yi) of the tree vertex Pi (i = 1 to 6) for each of the trees T1 to T6 shown in FIG.

The tree crown reflection intensity extracting unit 13 reads the coordinates (Xi, Yi) of the tree vertex Pi from the tree vertex storage unit 23, and coincides with or is closest to the coordinates (Xi, Yi) of the tree vertex Pi. The intensity I _{first of the} last pulse is extracted as the reflection intensity (hereinafter referred to as “crown reflection intensity”) Ict of the laser reflected by the tree crown as shown in FIG. The tree crown reflection intensity Ict is stored in the tree crown reflection intensity storage unit 28.

  The canopy passage intensity extraction unit 14 refers to the coordinates (Xi, Yi) of the tree vertex Pi stored in the tree vertex storage unit 23, and coordinates (Xi, Yi) of the tree vertex Pi from the reflection intensity data storage unit 21. The height (Z value) Zd of the reflection point of the reflection intensity data at the coordinates (X, Y) closest to or closest to is read. Further, the canopy passage strength extracting unit 14 reads out the ground height Zh at the coordinates (X, Y) closest to the coordinates (Xi, Yi) of the tree vertex Pi of the DEM from the DEM storage unit 24. Further, the canopy passage strength extracting unit 14 calculates a difference value Zk between the height Zd of the reflection point and the height Zh of the ground.

Further, the canopy passage strength extraction unit 14 reads a threshold value (for example, 1 m) from the threshold value storage unit 25 and compares it with the difference value Zk. When the difference value Zk is within the threshold value, it is determined that the laser has passed through the canopy and reflected from the ground, so the last pulse intensity I _last of the reflection intensity data passes through the canopy as shown in FIG. Extracted as _Icb , the reflection intensity of the laser that was attenuated and then reflected off the ground (hereinafter referred to as “canopy passing intensity”). As shown in FIG. 3, it can be said that the canopy passage strength Icb has information inside the canopy because it is attenuated by passing through the canopy. The extracted canopy passage strength Icb is stored in the canopy passage strength storage unit 29.

  The ground reflection intensity extraction unit 11 reads the reflection point height (Z value) Zd of the reflection intensity data in the forest area A from the reflection intensity data storage unit 21. Further, the ground reflection intensity extraction unit 11 reads the DEM ground height Zh in the forest area A from the DEM storage unit 24. Further, the ground reflection intensity extraction unit 11 calculates a difference value Zk between the height Zd of the reflection point and the height Zh of the ground.

Further, the ground reflection intensity extraction unit 11 reads a threshold value (for example, 1 m) from the threshold value storage unit 25 and compares it with the difference value Zk. When the difference value Zk is smaller than the threshold value, it is determined that the laser has directly reached the ground. Therefore, the first pulse intensity I _first of the reflected intensity data is directly reflected on the ground as shown in FIG. The reflection intensity (hereinafter referred to as “ground reflection intensity”) Ib is extracted. The extracted canopy passage strength Icb is stored in the canopy passage strength storage unit 29. Since the canopy passage intensity Icb is obtained from the reflectance Rb in the forest area A, the ground reflection intensity Ib is also extracted from the forest area A.

  The reflectance calculation unit 12 reads the incident intensity I0 from the incident intensity storage unit 27, reads the ground reflection intensity Ib in the forest area A from the ground reflection intensity storage unit 26, and uses the formula (2) to calculate the forest area A The reflectance Rb is calculated.

Rb = Ib / I0 (2)
The reflectance Rb in the forest area A is stored in the reflectance storage unit 22.

  As shown in FIG. 6, the ground reflection intensity Ib, the crown reflection intensity Ict, the canopy passage intensity Icb, and the incident intensity I0 in the forest area A are set for each of the trees T1 to T6.

  The leaf area index calculation unit 15 includes a ground reflection intensity storage unit 26, an incident intensity storage unit 27, a canopy reflection intensity storage unit 28, and a canopy passage intensity storage unit 29. Then, the canopy reflection intensity Ict, the canopy passage intensity Icb and the incident intensity I0 are read out, and the LAI is calculated for each of the trees T1 to T6.

  The ground reflection intensity Ib is expressed by the following formula.

Ib = I0 · Rb (3)
With reference to the Beer-Lambert law, the gap fraction is expressed by the following equation (4).

T (z) = exp (−G · μ · z) (4)
Here, T is the leaf gap, G is the leaf projection coefficient, μ is the leaf group density, and z is the transmission distance from the treetop.

  The canopy passage strength Icb is expressed by the following equation (5).

Icb = (I0-Ict) ・ exp (-2G ・ μ ・ d) ・ Rb (5)
Here, d represents the crown thickness. From equation (5), if the ratio of the canopy passage intensity Icb and the ground reflection intensity Ib in the forest area A is taken, equation (6) is obtained.

Icb / Ib = ((I0−Ict) / I0) • exp (−2G · μ · d) (6)
The leaf area index is expressed by the following formula (7).

LAI = μ · d (7)
Assuming that the leaves are projected on the horizontal plane, G = 1. From Formula (6) and Formula (7), the leaf area index LAI can be obtained by Formula (8).

LAI = -1 / 2 · ln ((Icb · I0) / (Ib · (I0-Ict))) ... (8)
The leaf area index calculation unit 15 performs ground reflection in the forest area A from the ground reflection intensity storage unit 26, the incident intensity storage unit 27, the crown reflection intensity storage unit 28, and the canopy passage intensity storage unit 29 for each of the trees T1 to T6. The intensity Ib, the canopy reflection intensity Ict, the canopy passage intensity Icb, and the incident intensity I0 are read out. Further, the leaf area index calculation unit 15 calculates LAI as shown in Expression (8). The calculated LAI is stored in the leaf area index storage unit 31. As shown in FIG. 7, the leaf area index storage unit 31 stores the LAI value for each of the trees T1 to T6. The output unit 4 records the LAI value on a recording medium, displays it on a screen, and prints it.

  Next, an example of the leaf area index calculation method according to the first embodiment of the present invention will be described with reference to the flowchart of FIG.

  (A) In step S11, the ground reflection intensity extraction unit 11 extracts, from the reflection intensity data storage unit 21, the reflection intensity (ground reflection intensity) Ib of the laser that has directly reached the ground as shown in FIG. The ground reflection intensity Ib in the forest area A is stored in the ground reflection intensity storage unit 26.

  (B) In step S 12, the ground reflection intensity extraction unit 6 extracts the ground reflection intensity Ia in the known reflectance region as shown in FIG. 3 from the reflection intensity data storage unit 21 and stores it in the ground reflection intensity storage unit 8. Let Then, the incident intensity calculation unit 7 reads the ground reflection intensity Ia of the known reflectance region from the ground reflection intensity storage unit 8 and further reads the known reflectance Ra from the reflectance storage unit 9 to obtain the equation (1). Using this, the incident intensity (incident intensity) I0 of the laser immediately before reaching the crowns of the trees T1 to T6 is calculated. The incident intensity I0 is stored in the incident intensity storage unit 27.

  (C) In step S13, the crown reflection intensity extraction unit 13 calculates the reflection intensity (crown reflection intensity) Ict of the laser reflected from the tree of the trees T1 to T6 as shown in FIG. Extracted for each of the trees T1 to T6. The tree crown reflection intensity Ict is stored in the tree crown reflection intensity storage unit 28.

  (D) In step S14, the canopy passage intensity extracting unit 14 reflects the reflection intensity of the laser reflected from the ground after passing through the canopy of the trees T1 to T6 as shown in FIG. (Passing intensity) Icb is extracted for each of the trees T1 to T6. The canopy passage strength Icb is stored in the canopy passage strength storage unit 29.

  (E) In step S15, the leaf area index calculation unit 15 determines the ground reflection intensity in the forest area A from the ground reflection intensity storage unit 26, the incident intensity storage unit 27, the canopy reflection intensity storage unit 28, and the canopy passage intensity storage unit 29. Ib, canopy reflection intensity Ict, canopy passage intensity Icb and incident intensity I0 are read out, and LAI is calculated for each tree T1 to T6 using equation (7). The LAI is stored in the leaf area index storage unit 31.

  Note that the procedures of steps S11 to S14 can be performed independently, and are not particularly limited to the order shown in the flowchart of FIG.

  Next, details of the method for setting the ground reflection intensity Ib in the forest area A in step S11 of FIG. 8 will be described with reference to the flowchart of FIG.

  (A) In step S21, the display unit 2 displays an orthophoto image on the screen. A user inputs an instruction to the forest area A of the orthophoto image as an LAI calculation target to the input unit 3. The display processing unit 10 reads pixel coordinates (Xp, Yp) on the screen indicating the forest area A in response to an instruction from the input unit 3. Then, the display processing unit 10 reads the plane rectangular coordinates (Xo, Yo) on the orthophoto image corresponding to the pixel coordinates (Xp, Yp) and passes them to the ground reflection intensity extracting unit 11.

  (B) In step S22, the ground reflection intensity extraction unit 11 assigns the reflection intensity data in the forest area A. In step S23, the ground reflection intensity extraction unit 11 extracts the height Zd of the reflection point of the laser included in the reflection intensity data.

  (C) In step S <b> 24, the ground reflection intensity extraction unit 11 reads the DEM ground height Zh in the forest area A from the DEM storage unit 24. In step S25, the ground reflection intensity extracting unit 11 takes a difference value Zk between the reflection point height Zd of the reflection intensity data and the DEM ground height Zh.

  (D) In step S26, the ground reflection intensity extraction unit 11 compares the difference value Zk with a threshold value. If the difference value Zk is equal to or greater than the threshold value, the process returns to step S24, and other reflection intensity data in the forest area A is allocated. On the other hand, if the difference value Zk is within the threshold value, the process proceeds to step S27.

(E) In step S27, the ground reflection intensity extraction unit 11 extracts the first pulse intensity I _first of the reflection intensity data as the ground reflection intensity Ib as shown in FIG. In step S28, the ground reflection intensity extraction unit 11 stores the ground reflection intensity Ib in the ground reflection intensity storage unit 26.

  Next, details of the method for setting the crown reflection intensity Ict in step S13 of FIG. 8 will be described with reference to the flowchart of FIG.

  (A) In step S31, the display unit 2 displays an orthophoto image. In response to an instruction from the input unit 3, the display processing unit 10 reads out pixel coordinates (Xp, Yp) indicating the forest area A displayed on the screen of the display unit 2. The display processing unit 10 reads the plane rectangular coordinates (Xo, Yo) of the orthophoto image corresponding to the pixel coordinates (Xp, Yp) and delivers them to the tree crown reflection intensity extracting unit 13.

  (B) In step S32, the tree crown reflection intensity extraction unit 13 determines the tree vertices Pi (P1) of the trees T1 to T6 existing in the forest area A from the tree vertex storage unit 23 based on the plane rectangular coordinates (Xo, Yo). , P2,... P6).

  (C) In step S33, the tree crown reflection intensity extraction unit 13 pays attention to the tree vertex P1 among the tree vertices Pi (P1, P2,... P6).

(D) In step S34, the crown reflection intensity extraction unit 13 determines the intensity I _first of the first pulse of the reflection intensity data having the coordinates (X, Y1) that are coincident with or closest to the coordinates (X1, Y1) of the tree vertex P1. Extracted as canopy reflection intensity Ict.

  (E) In step S35, the canopy reflection intensity extracting unit 13 stores the extracted canopy reflection intensity Ict in the canopy reflection intensity storage unit 28 in association with the tree vertex P1 of the tree T1.

  (F) In step S36, the crown reflection intensity extraction unit 13 determines whether the focused tree vertex Pi is the last tree vertex P6, that is, the tree crown reflection intensity for all tree vertices Pi (P1, P2,... P6). Determine whether Ict is set. If the crown reflection intensity Ict is not extracted for the tree vertices Pi (P1, P2,... P6) of all the trees T1 to T6, i is incremented by 1 in step S37, and the process returns to step S34 to return to the next tree vertex. The crown reflection intensity Ict is extracted for P2. As a result, the crown reflection intensity Ict is extracted for the tree vertices Pi (P1, P2,... P6) of all the trees T1 to T6.

  Next, details of the method for setting the canopy passage strength Icb in step S14 of FIG. 8 will be described with reference to the flowchart of FIG.

  (A) In step S41, the display unit 2 displays an orthophoto image. In response to an instruction from the input unit 3, the display processing unit 10 reads out pixel coordinates (Xp, Yp) indicating the forest area A of the tree group that is the target of the LAI calculation displayed on the screen of the display unit 2. The display processing unit 10 reads the plane rectangular coordinates (Xo, Yo) corresponding to the pixel coordinates (Xp, Yp) and delivers them to the tree crown passing intensity extracting unit 14.

  (B) In step S42, the canopy passage strength extracting unit 14 reads out the tree vertices Pi (P1, P2,... P6) of the trees T1 to T6 existing in the forest area A from the tree vertex storage unit 23.

  (C) In step S43, the canopy passage strength extracting unit 14 pays attention to the tree vertex P1 among the read tree vertices Pi (P1, P2,... P6).

  (D) In step S44, the tree crown passing intensity extracting unit 14 assigns the reflection intensity data having the coordinates (X, Y1) that coincide with or are closest to the coordinates (X1, Y1) of the tree vertex P1.

  (E) In step S45, the canopy passage intensity extracting unit 14 extracts the height Zd of the reflection point of the reflection intensity data.

  (F) In step S46, the canopy passage strength extracting unit 14 reads the ground height Zh of the DEM that is coincident with or closest to the coordinates (X1, Y1) of the tree vertex P1.

  (G) In step S47, the canopy passage intensity extracting unit 14 takes the difference value Zk between the height Zd of the reflection point of the reflection intensity data and the height Zh of the ground of the DEM.

  (H) In step S48, the canopy passage strength extraction unit 14 determines whether the difference value Zk is equal to or greater than a threshold value. If it is larger than the threshold value, the process proceeds to step S49, and reflection intensity data having coordinates next closest to the coordinates (Xi, Yi) of the tree vertex Pi are assigned, and the process returns to step S44. On the other hand, if it is within the threshold value, the process proceeds to step S50, and the canopy passage intensity extracting unit 14 allocates reflection intensity data having Zd.

(R) In step S51, the canopy passage intensity extraction unit 14 extracts the intensity I _last of the last pulse of the reflection intensity data as the canopy passage intensity Icb.

  (Nu) In step S52, the canopy passage strength extraction unit 14 stores (sets) the canopy passage strength Icb in the canopy passage strength storage unit 29 in association with the tree vertex P1.

  (L) In step S53, the canopy passage strength extraction unit 14 pays attention to the tree vertex P6, that is, the crown passage strength Icb for the tree vertices Pi (P1, P2,... P6) of all the trees T1 to T6. Determine whether or not. If the crown passing strength Icb has not been set for all the tree vertices Pi (P1, P2,... P6), the process proceeds to step S54, and the process returns to step S44 for the next tree vertex P2. As a result, the crown passage strength Icb is extracted for the tree vertices Pi (P1, P2,... P6) of all the trees T1 to T6.

  According to the first embodiment of the present invention, conventionally, only local LAI can be obtained by field survey. However, since LAI can be estimated by aviation laser measurement, wide area and quickness can be achieved without entering the site. LAI can be estimated.

  Furthermore, LAI is estimated using the information inside the tree canopy (attenuation of laser reflection intensity) by utilizing the laser intensity such as ground reflection intensity Ib, canopy reflection intensity Ict, canopy passage intensity Icb, and incident intensity I0. Therefore, it is possible to estimate the LAI with high accuracy.

  Further, using the tree vertex Pi stored in the tree vertex storage unit 23, the crown reflection intensity extraction unit 13 extracts the tree crown reflection intensity Ict at the coordinates (X, Y) that coincide with or is closest to the tree vertex Pi. The LAI can be estimated by extracting the canopy passage strength Icb having the coordinates (X, Y) that coincides with or is closest to the tree vertex Pi. A method for determining the tree vertex Pi will be described later.

The LAI obtained in this way can also be used as a parameter when estimating the amount of CO ₂ absorbed in the forest.

(Second Embodiment)
As a second embodiment of the present invention, a method for specifying a tree vertex will be described. The tree vertex Pi included in the tree vertex storage unit 23 illustrated in FIG. 13 is specified as follows.

  In the second embodiment of the present invention, a tree crown is emphasized and reproduced to recognize a vertex of a single tree. In order to recognize the vertices of this single tree, the crown shape is discriminated by newly creating “ridge valley degree” data.

  The ridge valley degree is calculated as follows.

  In the second embodiment of the present invention, using the tree height grid data (0.5 m), which is the difference between DEM and DSM, the search is performed in 8 directions (= all directions) and the ground opening (φ1) And the underground opening (φ2). The search range was 4m, 1m, 2m, 5m, and 10m, and was determined to be 2m that could properly represent the crown shape. Then, the bisector of the angle between the ground opening and the underground opening representing the rough topography in the search range 2 m is defined as the ridge valley degree (φ3).

  Embodiments of tree recognition will be described below with reference to the drawings.

  FIG. 12 is a schematic configuration diagram of a tree vertex recognition apparatus according to the second embodiment of this invention. As shown in FIG. 12, a database 110 storing grid data (DEM: Digital Elevation Model) of the ground in a specific area and a database 111 storing grid data (DSM: Digital Surface Model) of the surface layer (forest) of the specific area. A database 112 storing digital photographs of a specific area, a display unit 13, a tree height calculation unit 114, a file (memory) 115 storing tree height data, and a smoothing filter unit 116 (3 pixels × 3 Pixel), a ridge valley degree calculation unit 117, a file 118 in which the ridge valley degree is stored, a display processing unit 119, a local maximum filter unit 120, a tree vertex candidate extraction unit 121, and tree vertex candidate data are stored. File 122, tree crown upper extraction unit 123, file 124 in which tree crown upper data is stored, and tree vertex determination unit 1 25, a file 126 in which tree vertex data is stored, an output unit 130, and the like.

  The database 110 stores grid data (DEM: Digital Elevation Model: grid-like data at intervals of 0.5 m in the second embodiment of the present invention) of ground in a specific area. This grid data is stored in association with X, Y, and Z (ground height Za) of the ground. That is, the laser data is read, and the ground is restored by generating a TIN for the contour maps connecting the same elevation values.

  The database 111 stores grid data (DSM: Digital Surface Model: 0.5 m interval in the second embodiment of the present invention) of the surface layer (forest) of a specific area. This grid data (in the second embodiment of the present invention, grid-like data with an interval of 0.5 m) stores X, Y, and Z (Zb: tree surface height) of the tree surface in association with each other. . That is, the laser data is read and the surface layer is restored by generating a TIN for the contour map connecting the same elevation values.

  That is, in these grid data, an X value and a Y value corresponding to latitude and longitude and a Z value obtained by laser measurement are assigned to an XY plane rectangular coordinate system having a lattice interval of 0.5 m.

  The tree height calculation unit 114 obtains a difference between the DSM of the database 110 and the DEM of the database 111 and assigns the difference to the X and Y coordinates of the file 115 (memory).

  The smoothing filter unit 116 is also referred to as an average filter, and with respect to the tree height grid data Di of the file 115, an XY coordinate plane 3 × 3 range centering on each data (second embodiment of the present invention). In this embodiment, the average is obtained by 1.5 m × 1.5 m) and stored in the file 126.

  When the tree height grid data Ddi (Dd1, Dd2,...) Is stored in the file 126, the ridge valley degree calculation unit 117 stores any averaged tree to which the tree height grid data Ddi is added. Search range of arbitrary radius of XY coordinate plane centering on each data from high grid data Ddi (hereinafter referred to as tree high grid data Ddi) (search range is 2m which can properly express the characteristics of the solid shape of tree crown) The degree of ridge valley is obtained at every 2 m range below and stored in the file 118.

  The local maximum value filter unit 120 sequentially applies a 3 × 3 range filter to the data Dei of the ridge valley degree of the file 118 and stores the result in the file 127.

  The tree vertex candidate extraction unit 121 compares the crown three-dimensional shape feature data Dfi (XY, degree of ridge valley) of the file 127 with Dei (Xpi, Ypi, degree of ridge valley) of the file 118, and the same is obtained as the vertex candidate Dgi. As a file 122.

  The crown top extraction unit 123 extracts the data Dei at the top of the crown by providing a threshold corresponding to the state of the forest in the ridge valley degree data (crown solid shape feature data Dei) of the file 118, and extracts this data Data Dhi (Dh1, Dh2,...) Is stored in the file 124.

  The tree vertex determination unit 125 determines a superposition of the data Dhi (Dh1, Dh2,...) And the vertex candidate Dgi (Dg1, Dg2,. This is saved in the file 128. When there are a plurality of temporary vertices at the top of the tree crown, the highest tree height at the top of the tree crown is determined as a single tree vertex and stored in the file 129 (Pi).

  The display processing unit 119 displays the determined vertex on the screen with, for example, a circle (identification of a triangle, a cross, a number, etc.). At this time, it may be displayed so as to overlap with the photographic image.

  The output unit 130 attaches a tree number to the coordinates (XY) of the vertex Pi of the file 129 and displays it on a printer (not shown) or a screen. At this time, the tree height may also be displayed.

(Description of operation)
The operation of the tree vertex specifying apparatus configured as described above will be described below. In the second embodiment of the present invention, the display unit 13 displays the DSM of the database 110, the DEM of the database 111, and the digital photograph in an overlapping manner.

  FIG. 13 is a flowchart for explaining the schematic operation of the tree vertex recognition apparatus according to the second embodiment of this invention.

  In the second embodiment of the present invention, the tree height calculation unit 114 reads data (DSM: Digital Surface Model) acquired by aviation laser measurement (S71). Further, the grid data of the ground (hereinafter referred to as DEM: Digital Elevation Model) is read (S72).

  Then, the tree height calculation unit 114 takes the difference Zci (Dci, Zci,...) Between the DSM and DEM data and stores it as the tree height grid data Di in the file 115 (S73).

  A smoothing filter 116 (also referred to as an average filter 16) sequentially averages the tree height grid data Di in a 3 × 3 range and stores it in the file 126 (tree height grid data Ddi (Dd1, Dd2,...). )) (S74).

  Next, the ridge valley degree calculation part 117 calculates the ground opening degree and the underground opening degree for every 2 m range (Ddi) (S75). Next, the ridge valley degree calculation unit 117 uses the underground opening degree and the ground opening degree to create a “ridge valley degree” data, thereby obtaining a ridge valley image that can distinguish the characteristics of the three-dimensional shape of the crown. It is stored in 118 (S76).

  The ridge valley degree is calculated as follows. Using the grid data Ddi of the tree height, search is made in eight directions (= all directions) to obtain the ground opening (φ1) and the underground opening (φ2).

  The search range was 4m, 1m, 2m, 5m, and 10m, and was determined to be 2m that could properly represent the crown shape.

  Then, the bisector of the angle between the ground opening and the underground opening representing the rough topography in the search range 2 m is defined as the ridge valley degree (φ3) (see FIG. 20). The data Dei of the ridge valley degree is displayed by the display processing unit 119 (ridge valley degree image).

  Further, the local maximum value filter unit 120 is sequentially applied to the data Dei of the ridge valley degree (S77), and is stored in the file 127.

  At this time, not only the size of the tree crown but also the 3 × 3 range which is the minimum area size is fixed, and “tree vertex candidate points” are extracted instead of the tree vertices, and tree vertex candidates Dgh (dg1,. (S78). The 3 × 3 range is fixed because, for example, when it is 5 × 5, an experimental result has been obtained that the upper part of the crown where there are no tree vertex candidates is obtained.

  On the other hand, using the property that “ridge valley degree> 0” expresses ridge topography and “ridge valley degree <0” expresses valley topography, the tree vertex candidate extraction unit 121 converts the forest into the ridge valley degree data. The upper part of the tree canopy is extracted by providing a threshold value according to the state (S79), and the upper crown extracting unit 124 is labeled (S80).

  Further, the tree vertex determining unit 125 obtains the highest tree height from the tree vertex candidate points for each of the treetop candidate points and the upper tree crown area, determines that point as the tree vertex (S81), and extracts this as the tree vertex (S82). The determination of the tree vertex will be described later.

  FIG. 14 is a flowchart illustrating processing for obtaining the ridge valley degree of the tree vertex recognition apparatus according to the second embodiment of this invention.

  In the second embodiment of the present invention, an orthophoto image is displayed on the display unit 13 (S90).

  The tree height calculation unit 114 reads DSM (X, Y, Z) and DEM (X, Y, Z) (S91), obtains a difference Zci (S92), and obtains the tree height grid data Di (D1, D2). D2... Are stored in the file 115.

  As shown in FIG. 15, the tree height data Di is data in which coordinates (X, Y: plane rectangular coordinates) and the tree height Zci are associated with each other.

  Next, the smoothing filter unit 116 sequentially averages (3 × 3 range) the grid data Di of a number of tree heights in the file 115 by the filter 16 (S93), and the grid data Ddi of the tree height is stored in the file 126. save.

  Next, the ridge valley degree calculation unit 117 determines the feature of the three-dimensional shape of the tree crown by sequentially obtaining the ridge valley degree in the range of 2 m in each Ddi (S94).

  The ridge valley degree is searched in 8 directions (may be performed at 1m, 5m, 10m depending on the type of tree) to obtain the ground opening and the underground opening, and the difference between the ground opening and the underground opening is 2 It is obtained by dividing by (S94). That is, the feature of the three-dimensional shape of the crown (leaves, branches) of the tree is required.

  Using the tree height grid data Ddi of the file 126, the ground opening (φ1) and the underground opening (φ2) are obtained by searching in eight directions (= all directions) as shown in FIG. The search range is determined to be 2m, which can properly represent the three-dimensional shape of the tree crown.

  FIG. 17 is a plot of the tree height grid data Ddi in an arbitrary search range on an axis composed of distance and height. In addition, black dots in FIG. 17 indicate grid data Ddi at intervals of 0.5 m (2 m range is not shown).

That is, FIG. 17 shows a rough crown shape when one direction is viewed within a range of an arbitrary radius from the sample point of interest. In FIG. 17, φ1 is the ground opening (ground angle), and φ2 is the underground opening (ground angle). 18 (a) and 18 (b), the angle between the bisector of the angle between the ground opening and the underground opening and the distance axis (horizontal axis) is defined as the ridge valley degree (φ3). ing. Further, black dots in FIG. 18 indicate grid data Ddi at intervals of 0.5 m (2 m range is not shown).
In this embodiment, the ridge valley degree (φ3) is obtained with the search range set to 2 m, and is stored in the file 118 as ridge valley degree data (hereinafter, crown three-dimensional shape feature data Dei (Dei: De1, De2...)) (S96). ).

  FIG. 19 is an explanatory diagram representing the crown shape of the locus Ai with the tree height as the Y axis and the distance as the X axis. FIG. 20 is an explanatory diagram expressing the characteristics of the three-dimensional shape of the crown of the trajectory Ai with the ridge valley degree as the Y axis and the distance as the X axis.

  That is, as shown in FIG. 20, when the feature of the three-dimensional shape of the crown of the locus Ai is expressed by the ridge valley degree, the boundary between the crowns appears more conspicuously than in FIG. Numbers (1) to (6) in FIG. 19 and FIG. 20 indicate the crown area read from the photograph with the naked eye.

  In order to improve the ease of recognition as a single tree, the difference between the values of the tree vertices and the edge part is large, the values show the same tendency at the tree vertices, the values show the same tendency at the edges, etc. It is done. As shown in FIG. 20, it can be seen that the “ridge valley degree” data of the present method satisfies the above-mentioned conditions than the conventional “height” data. Thereby, it can contribute to the accuracy improvement at the time of creating a crown crown polygon automatically.

  In FIG. 20, it can be seen that the range from the point where the ridge valley degree has a negative value to the point where it suddenly increases in the + direction and then suddenly becomes a negative value is the crown of a single tree.

  FIG. 22 is a diagram for explaining a method for selecting vertex candidates and tree vertices, and shows that a candidate having the highest tree height is selected at the top of each tree crown. The determination of the tree vertex will be described later in detail.

  Next, the ridge valley degree calculation unit 117 stores the obtained ridge valley degree in the file 118 (XY is associated) as the crown crown shape feature data Dei (see FIG. 21) (S96).

  The above-mentioned ground opening is searched for each element (Ddi) in eight directions from the element, and the maximum inclination angle is obtained within the consideration distance. The average of the maximum inclination angles in each of the 8 directions is expressed as the angle from the zenith. Note that the value of the ground opening depends on the range to be considered and must be set according to the scale of the object to be read.

  The underground opening is searched for each element in eight directions from the element, and the minimum inclination angle is obtained within the considered distance. The average of the minimum inclination angles in each of the 8 directions is expressed as the angle from the zenith. As with the ground opening, the value depends on the range to be considered and must be set according to the scale of the object to be interpreted.

  The ridge valley degree is calculated by subtracting the underground opening from the calculated ground opening and dividing by 2. Ridge terrain has a positive value and valley terrain has a negative value.

  And it progresses to the tree vertex extraction process of FIG. In step S7 of the tree vertex extraction process, the 3 × 3 local maximum value filter unit 120 (using a window of 3 × 3 range) is applied to the feature data Dei (De1, De2,...) Of the crown shape of the file 118. Then, this is saved (Dfi) in the file 127.

  In step S78, the tree vertex candidate extracting unit 121 uses the crown point feature data Dfi (see FIG. 22) of the peak point existing in the crown area ((1), (2),...) As a vertex candidate (hereinafter referred to as Dgi). Extracted and stored in the file 122 (see FIG. 23). Specifically, the value of the data of the file 118 and the value of the data of the same coordinate value of the file 127 are defined in the file 122 as “1” and the others as “0” (see FIG. 23).

  Here, the description of the processes of the local maximum filter unit 120 and the tree vertex candidate extraction unit 121 will be supplemented.

  That is, the local maximum filter unit 120 (also referred to as Local max filter 20) applies a window of 3 × 3 range to the feature data Dei (De1, De2,...) Of the tree crown shape of the file 118, and becomes the maximum in this window. The process of replacing the value with the median value of the window is sequentially performed. These are stored in the file 127 as feature data Dfi of the canopy solid shape.

  Then, the tree vertex candidate calculation unit 21 compares the crown three-dimensional shape feature data Dfi (XY, degree of ridge valley) and Dei (XY, degree of ridge valley) of the file 118, and the same data is stored in the file 122 as the vertex candidate Dgi. Saved.

  On the other hand, in step S79, the upper crown extraction unit 123 uses the property that “ridge valley degree> 0” represents the ridge landform and “ridge valley degree <0” represents the valley terrain. The data Dei at the upper part of the canopy is extracted by setting a threshold value corresponding to the forest condition in the degree data (feature data Dei of the three-dimensional shape of the crown), and this is filed as data Dhi (Dh1, Dh2,...) 124 (see FIG. 24).

  The data Dhi at the top of the tree crown in FIG. 24A is obtained by extracting Dei exceeding the threshold. For example, the data Dei in the area of number 1 is (Xh5, Yh5) is “1”, (Xh6, Yh6) Is represented as “1” (Xh10, Yh10) as “1”.

  FIG. 24B specifically shows FIG. 24A, and the vertex that exceeds the threshold is set to “1”. The range where “1” continues is the upper part of the tree crown.

  In step S80, Xhi and Yhi of the vertex candidates Dhi in the file 124 are polygonized (coordinates Xhi and Yhi and color information), and the display processing unit 119 assigns a label to the grouped area (upper area). (See FIG. 25).

    Next, in step S81, the tree vertex determining unit 125 superimposes and matches the data Dhi (Dh1, Dh2,...) On the top of the tree of the file 124 and the vertex candidate Dgi (Dg1, Dg2,...) Of the file 122. Is determined as a temporary vertex. When there are a plurality of temporary vertices in the upper area of the tree crown, the area with the highest tree height in the upper area of the tree crown (the areas of numbers 1, 2,... In FIGS. 19 and 20). It is determined as a vertex of a single tree, and this is stored in the file 129 (Pi). The display processing unit 119 displays the determined vertex on the screen (see FIG. 26). That is, for example, the vertex candidate Dgi is displayed as shown in FIG. 27, but the vertex-determined image is as shown in FIG.

  The explanation of the above-described vertex determination will be supplemented. For example, as shown in FIG. 29A, in the file 124, the first column is “00110011100”, the second column is “011111011110”, the third column is “001110001000”, and the fourth column is “00000100011”. When the top apex area is obtained, as shown in FIG. 29 (b), a serial number is assigned to each group, with “1” adjacent groups (“0” separated) as one crown top. To go.

  For example, when a plurality of vertices having the same height are present in the crown, the data DQi of the file 128 is displayed once and designated by the operator, and the designated one is determined as the vertex Pi. Save to file 129.

  Next, the explanation about the ridge valley is supplemented. As the characteristics of the ridge valley degree, considering the terrain of the sin curve as shown in FIG. 30, the ridge terrain is represented by “ridge valley degree> 0”, the valley terrain is represented by “ridge valley degree <0”, By the way, it is mentioned that “degree of ridge valley = 0”.

  Here, in order to make it easy to understand, if the ideal tree has a shape as shown in FIG. 30, “ridge valley degree> 0” at the top of the tree, and always “ridge valley degree = 0”. It reaches the edge of the canopy by passing through the point, where it is considered that the “ridge valley degree <0”. Therefore, it can be said that the upper part of the canopy can be defined as “degree of ridge valley> 0”.

  Utilizing this fact, in this work, tree vertices are narrowed down as shown in FIG. In other words, the top of the tree crown is extracted by using the value of the ridge valley degree, and the candidate point having the highest tree height (DSM-DEM) is selected as the tree vertex among the candidate tree top points existing in this extraction. It was.

  The explanation is supplemented about the threshold value of the ridge valley degree.

  When the threshold value is set for the ridge valley degree and the upper part of the tree crown is extracted, it is understood that the tree 5 cannot be extracted in the case of the threshold value α as shown in FIG. On the other hand, in the case of the threshold value β, the trees 5 and 6 are recognized as one. Thus, there is a problem that the tree vertex extracted by the threshold is affected.

Therefore, we decided to examine the optimum threshold for extracting the upper part of the tree crown. Specifically, the threshold value of the ridge valley degree is set to five types “10, 15, 20, 25, 30”, and the one having the best tree vertex extraction result is adopted. Table 1 shows the total results of the tree vertices after narrowing down.

  Looking at Table 1, it can be seen that the result with “degree of ridge valley> 20” was the best. Compared to this, the reason why the number of extracted trees with “ridge ridge degree> 10” is small is that the adjacent tree crown was recognized as one. Also, the reason why the number of extractions with “ridge ridge degree> 30” is small is that there was a crown that could not be extracted because the threshold was set high.

Therefore, it was found that “20” is the optimum threshold for the ridge valley degree for extracting the upper part of the tree crown. The extracted image of the upper part of the tree crown is shown in FIG. 25, and the extracted tree vertex extracted image is shown in FIG.

By calculating the apex of the canopy of a certain forest as described above, a single tree can be specified, so that it can be applied to management of forest assets. For example, when the output unit 130 outputs and prints the vertex data Pi (X, Y, Z) of the file 129, the number and height are known, so that the value of the forest can be automatically grasped.

  In FIG. 33, the crown polygon is created after calculating the tree vertex. In creating the canopy polygon, a canopy polygon (see FIG. 34) is generated by the Wetershed algorithm using the tree vertices extracted in step S82 and the ridge valley degree of the file 118 (S83, S85).

  The threshold value determination method in step S83 is as follows.

  An arbitrary line (trajectory Ai) is drawn in the orthophoto image as in the photographic image of FIG. 14, and the ridge valley degree calculation unit creates a cross-sectional view of the ridge valley degree and displays it (corresponding to FIG. 20).

  Next, an additional distance from the start of the crossing line (corresponding to the X axis in FIG. 20) is followed, and a tree crown and a boundary between the tree crowns in the display unit are searched.

  Then, looking at the additional distance of the boundary position, the operator reads the ridge valley degree (calculated and displayed by the ridge valley degree calculation unit) of the additional distance from FIG.

  This process is repeated to read the degree of ridge valley where the upper crown is formed.

  20 (1) 0 and -50, (2) -50 and -20, (3) -20 and 0, (4) 0 and 0, (5) 0 and 10, (6) Then it became 10 and -10.

  From this result, the ridge valley degree 10 in which (1) to (6) are separated is set as a threshold value.

  If the threshold value is 0, (5) and (6) will be at the top of one tree crown, so 0 is not adopted.

  If the threshold is 30 (2), two upper crowns will be formed, and (5) will not be the upper crown. Therefore, in FIG. 20, a ridge valley degree of about 10 is suitable as a threshold value. The threshold value thus determined is used.

  Although the above embodiment has been described using aerial photographs, gray scale may be assigned to the tree height grid data, DEM, or DSM.

  For example, the difference image between the ground opening and the underground opening is gray, the slope is put in the red channel, and by creating a pseudo color image, the ridge and the top of the mountain are whitish, and the valley and the depression are expressed blackish, The steeper slope is expressed in red. With such a combination of expressions, even one image having a stereoscopic effect (also referred to as a stereoscopic red map) may be used.

  For example, when the value of the ground opening is within a range of about 40 degrees to 120 degrees, 50 degrees to 110 degrees is assigned to 255 gradations corresponding to the gray scale. That is, the ridge portion (convex portion) has a larger value of the ground opening, so the color becomes white.

  The second embodiment of the present invention is preferably realized by a computer system VPS as shown in FIG. The VPS 1 includes a central information processing device (CPU) composed of an appropriate combination of workstation, processor, microcomputer, logic, registers, etc., and a keyboard (KB) for inputting control / operation information necessary for the central information processing device, An information input unit including a mouse, an interactive soft switch, an external communication channel, etc., and an information output unit including a display, a printer, an external communication channel, etc. that display and transmit information output from the central information processing device in a broad sense A ROM that stores information such as an operating system and application programs read by the central information processing apparatus, and a ram that stores information to be processed by the central information processing apparatus and information to be written from the central information processing apparatus at any time. (RAM) and the like. It is possible to integrate and subdivide ROM and RAM as appropriate.

  The application program described above includes the tree height calculation unit 114, the ridge valley degree calculation unit 117, the display processing unit 119, and the local maximum filter unit 120 (also referred to as a Local max filter) according to the second embodiment of the present invention. A tree vertex candidate extraction unit 121, a tree crown upper extraction unit 123, a tree vertex determination unit 125, an output unit 130, and the like.

  In this embodiment, the concept of opening is used. The degree of opening is a quantification of the extent that the point is protruding above the ground and the extent of biting underground. In other words, the ground opening represents the size of the sky that can be seen within the distance L from the sample point of interest, and the underground opening is the basement in the range of distance L when looking up underground with a handstand. Represents the size.

  The opening degree depends on the distance L and the surrounding terrain. In general, the ground opening increases as the point protrudes higher from the surroundings, and takes a larger value at the summit and ridge and smaller at the depression and valley bottom. On the contrary, the underground opening becomes larger as the depth of penetration into the basement is greater, with a large value at the depression and valley bottom, and a small value at the summit and ridge.

  In addition, the opening degree map can extract information suitable for the size of the terrain by specifying the calculation distance, and can display without depending on directionality and local noise.

  In other words, it excels in the extraction of ridge lines and valley lines, and abundant topographical and geological information can be read, and one of the eight directions from the set point A was seen on a certain range of DEM data. An angle vector θi formed by the straight line L1 connecting the point B, which is sometimes the maximum vertex, and the horizontal line is obtained. This angle vector is calculated in eight directions, and the average of these is called the ground opening, and the solid with the air layer pressed against a certain range of DEM data (ground surface: solid) is turned over. Further, the angle formed by the straight line L2 connecting the point C (corresponding to the deepest point) that becomes the maximum vertex when any one of the eight directions from the point A of the inverted DEM data is viewed is obtained. Obtaining and averaging this angle over eight directions is referred to as underground opening ψi.

  That is, the ground opening is a maximum value (vertical) of the slope of the line L1 that generates a topographic cross section for each of the eight directions on the DEM data included in the range from the point of interest to a certain distance and connects each point to the point of interest. (When viewed from the direction). Such processing is performed for eight directions. The angle of inclination is the angle from the zenith (90 degrees for flats, 90 degrees or more for ridges and peaks, and 90 degrees or less for valleys and depressions). Then, a topographic section is generated for each of the eight directions, and the maximum value of the slope of the line connecting each point and the point of interest is obtained (the minimum value when L2 is viewed from the vertical direction in the three-dimensional map of the ground surface). Such processing is performed for eight directions.

  The angle ψi when viewing L2 from the vertical direction in the three-dimensional map of the ground surface is 90 degrees if flat, 90 degrees or less at the ridge or peak, and 90 degrees or more at the valley or depression.

  That is, the ground opening and the underground opening consider two basic points A (iA, jA, HA) and B (iB, jB, HB).

The distance between A and B is
P = {(i _A − i _B ) ² + (j _A − j _B ) ² } ^1/2 … (1)
It becomes.

The relationship between A and B at the sample point is shown with an altitude of 0 m as a reference. The elevation angle θ when viewing the sample point B from the sample point A is θ = tan ⁻¹ {(HB−HA) / P}.
Given in. The sign of θ is positive when (1) HA <HB, and negative when (2) HA> HB.

A set of sample points within the range of the azimuth D distance L from the sample point of interest is described as DSL, and this is referred to as a “DL set of sample points of interest”. here,
DβL: Maximum value of the elevation angle for each element of DSL at the sample point of interest
DδL: The following definition is made as the minimum value of the elevation angle for each element of DSL at the sample point of interest.

Definition 1: The ground angle and underground angle of the D-L set of sample points of interest are each
DφL = 90− DβL
as well as
DψL = 90 + DδL
Means.

  DφL means the maximum value of the zenith angle at which the sky in the direction D can be seen within a distance L from the sample point of interest. Generally speaking, the horizon angle corresponds to the ground angle when L is infinite. Further, DψL means the maximum value of the nadir angle at which the ground in the direction D can be seen within a distance L from the sample point of interest. When L is increased, the number of sample points belonging to DSL increases, so that DβL has a non-decreasing characteristic, while DδL has a non-increasing characteristic. Therefore, both DφL and DψL have non-increasing characteristics with respect to L.

  High angle in surveying is a concept defined with reference to a horizontal plane passing through the sample point of interest, and does not exactly match θ. In addition, if the ground angle and the underground angle are to be strictly discussed, the curvature of the earth must be taken into consideration, and definition 1 is not necessarily an accurate description. Definition 1 is a concept defined on the premise that terrain analysis is performed using DEM.

  The ground angle and the underground angle are concepts about the specified direction D, but the following definition is introduced as an extension of this.

Definition II: The ground opening and the underground opening at the distance L of the sample point of interest are ΦL = (0φL + 45φL + 90φL + 135φL + 180φL + 225φL + 270φL + 315φL) / 8
And ψL = (0ψL + 45ψL + 90ψL + 135ψL + 180ψL + 225ψL + 270ψL + 315ψL) / 8
Means.

  The ground opening represents the extent of the sky that can be seen within the distance L from the sample point of interest, and the underground opening represents the extent of the underground in the range of distance L when looking up underground with a handstand. Represents.

(Other embodiments)
As described above, the present invention has been described according to the first and second embodiments. However, it should not be understood that the description and drawings constituting a part of this disclosure limit the present invention. From this disclosure, various alternative embodiments, examples and operational techniques will be apparent to those skilled in the art.

For example, in the first embodiment of the present invention, the crown reflection intensity extracting unit 13 matches the coordinates (Xi, Yi) of the tree vertex Pi or the last pulse of the reflection intensity data at the closest coordinates (X, Y). In the above description, the intensity I _first is extracted as the canopy reflection intensity Ict, but instead of the tree vertex Pi, the crown upper part extracted by the crown upper data extraction unit 123 described in the second embodiment of the present invention. The data Dhi may be used. That is, the crown reflection intensity extraction unit 13 may extract the intensity I _first of the last pulse at the coordinates (X, Y) near the upper part of the tree crown as the crown reflection intensity Ict using the data Dhi of the upper part of the tree.

Further, in the first embodiment of the present invention, the tree crown passage intensity extracting unit 14 has the last pulse of the reflection intensity data at the coordinates (X, Y) closest to the coordinates (Xi, Yi) of the tree vertex Pi. In the above description, the intensity I _last of the tree is extracted as the canopy passage intensity I _cb , but instead of the tree vertex Pi, the tree crown extracted by the crown upper data extraction unit 123 described in the second embodiment of the present invention. The upper data Dhi may be used. That is, the crown passage intensity extraction unit 14 extracts the last pulse intensity I _last of the reflection intensity data of coordinates (X, Y) near the upper part of the tree crown as the crown passage intensity I _cb using the data Dhi of the upper part of the tree. May be.

  In addition, although the tree vertex determination method has been described in the second embodiment of the present invention, the tree vertex used in the first embodiment of the present invention is determined as described in the second embodiment of the present invention. The determination method may not be limited, and the determination method is not particularly limited. For example, as a method for identifying tree vertices, DSM and DEM are created by using a technique for visually identifying tree vertices such as photographic images or data acquired by aerial laser measurement, and the difference between DSM and DEM is obtained. Thus, a method may be used in which grid data of tree height is calculated, and tree vertices are specified using the tree height or the maximum value of the tree height as a discrimination index.

  As described above, the present invention naturally includes various embodiments not described herein. Therefore, the technical scope of the present invention is defined only by the invention specifying matters according to the scope of claims reasonable from the above description.

It is a specific hardware configuration diagram according to the first embodiment of the present invention. It is a schematic block diagram of the leaf area index | exponent calculation apparatus which concerns on the 1st Embodiment of this invention. It is explanatory drawing explaining a laser reflection intensity. It is a graph of spectral characteristics (reflectance). It is a top view for demonstrating a calculation area. It is a table | surface which shows the extracted canopy inside information (laser reflection intensity). It is a table | surface which shows the calculated leaf area index. It is a flowchart of the leaf area index calculation method according to the first embodiment of the present invention. It is a flowchart of the setting method of ground reflection intensity. It is a flowchart of the setting method of tree crown reflection intensity. It is a flowchart of the setting method of tree crown passage intensity. It is a schematic block diagram of the tree vertex recognition apparatus of the 2nd Embodiment of this invention. It is a flowchart explaining schematic operation | movement of the tree vertex recognition apparatus of the 2nd Embodiment of this invention. It is a flowchart explaining the process which calculates | requires the ridge valley degree of the tree vertex recognition apparatus of the 2nd Embodiment of this invention. It is explanatory drawing of the tree height data Di. It is explanatory drawing explaining calculation of a ridge valley degree. It is explanatory drawing explaining calculation of a ridge valley degree. It is explanatory drawing explaining calculation of a ridge valley degree. It is explanatory drawing expressing the crown shape of locus | trajectory Ai by using a tree height as a Y-axis and distance as an X-axis. It is explanatory drawing expressing the feature of the solid shape of the crown of locus | trajectory Ai by making a ridge valley degree into a Y-axis and making distance into an X-axis. It is explanatory drawing of tree high grid data. It is explanatory drawing explaining calculation of the vertex candidate by ridge valley degree. It is explanatory drawing of the vertex candidate Dgi. It is explanatory drawing of the data Dhi of the upper part of a tree crown. It is explanatory drawing explaining the labeling of a tree crown upper part. It is explanatory drawing of a screen about the determined vertex. It is explanatory drawing until a vertex determination. It is explanatory drawing until a vertex determination. It is explanatory drawing until a vertex determination. It is explanatory drawing of a ridge valley degree. It is explanatory drawing of a ridge valley degree. It is explanatory drawing explaining the relationship between a threshold value and a tree crown part. It is a flowchart explaining other embodiment. It is explanatory drawing explaining the polygon of the tree crown obtained by other embodiment. It is a concrete hardware block diagram of the 2nd Embodiment of this invention.

Explanation of symbols

1 ... CPU
DESCRIPTION OF SYMBOLS 2 ... Display part 3 ... Input part 4 ... Output part 5 ... Memory | storage device 10 ... Display processing part 11 ... Ground reflection intensity extraction part 12 ... I0 calculation part 13 ... Crown reflection intensity extraction part 14 ... Crown passage intensity extraction part 15 ... Leaf Area index calculation unit 16 ... smoothing filter (average filter)
DESCRIPTION OF SYMBOLS 20 ... Orthophoto image memory | storage part 21 ... Reflection intensity data memory | storage part 22 ... Reflectivity memory | storage part 23 ... Tree vertex memory | storage part 24 ... DEM memory | storage part 25 ... Threshold memory | storage part 26 ... Ground reflection intensity memory | storage part 27 ... Incident intensity memory | storage part 28 ... Crown reflection intensity storage unit 29 ... Crown passage intensity storage unit 31 ... Leaf area index storage unit 110 ... Database 111 ... Database 112 ... Database 113 ... Display unit 114 ... Tree height calculation unit 116 ... Smooth filter 117 ... Ridge valley degree calculation unit 118 ... File 119 ... Display processing unit 120 ... Local maximum value filter 121 ... Tree vertex candidate extraction unit 122 ... File 123 ... Upper crown data extraction unit 124 ... File 125 ... Tree vertex determination unit 126 ... File 127 ... File 128 ... File 129 ... File 130 ... Output unit

Claims (9)

Irradiating a laser from the sky to the ground, storing the reflected intensity data of the laser reflected on the ground;
Storing the reflectance of the ground;
Extracting the reflected intensity of the laser reflected by the tree crown from the stored reflected intensity data;
From the stored reflection intensity data, extracting the reflection intensity of the laser reflected through the tree crown and reflected from the ground;
Calculating the incident intensity of the laser immediately before reaching the tree crown based on the stored reflection intensity data and reflectance;
Extracting the reflected intensity of the laser reflected from the ground from the stored reflected intensity data;
Using the incident intensity of the laser immediately before reaching the tree crown, the reflected intensity of the laser that has passed through the crown, the reflected intensity of the laser reflected by the crown, and the reflected intensity of the laser reflected by the ground, the leaf area of the tree A leaf area index calculation method comprising: calculating an index. Calculating the incident intensity of the laser just before reaching the tree crown,
Extracting the reflection intensity of the laser reflected from the ground from the stored reflection intensity data,
2. The leaf area index calculation method according to claim 1, wherein the incident intensity of the laser immediately before reaching the tree crown is calculated by dividing the reflection intensity of the laser reflected by the ground by the reflectance of the ground. The reflection intensity data includes the intensity of the first pulse of the laser, the intensity of the last pulse of the laser, and the height of the reflection point of the laser;
Storing the height of the ground;
Storing a specific threshold;
Comparing the threshold value with the difference value between the height of the reflection point and the height of the ground;
As a result of the comparison, extracting the intensity of the last pulse included in the reflection intensity data whose difference value is within the threshold as the reflection intensity of the laser that has passed through the tree crown;
And extracting the intensity of the first pulse included in the reflection intensity data whose difference value is smaller than the threshold as a result of the comparison as a reflection intensity of the laser reflected by the ground. Item 3. The leaf area index calculation method according to Item 1 or 2.   The reflectance is calculated by dividing the reflection intensity of the laser reflected by the ground by the incident intensity of the laser immediately before reaching the crown of the tree. The leaf area index calculation method described in 1. (Program. Claims 5 to 8 are the same as claims 1 to 4.)
Computer
Means for irradiating a laser from the sky to the ground and storing the reflected intensity data of the laser reflected on the ground;
Means for storing the reflectance of the ground;
Means for extracting the reflection intensity of the laser reflected by the crown of the tree from the stored reflection intensity data;
Means for extracting the reflection intensity of the laser reflected by the ground through the tree crown from the stored reflection intensity data;
Means for calculating the incident intensity of the laser immediately before reaching the tree crown based on the stored reflection intensity data and reflectance;
Means for extracting the reflected intensity of the laser reflected from the ground from the stored reflected intensity data;
Using the incident intensity of the laser immediately before reaching the tree crown, the reflected intensity of the laser that has passed through the crown, the reflected intensity of the laser reflected by the crown, and the reflected intensity of the laser reflected by the ground, the leaf area of the tree A program that functions as a means of calculating an index. The computer,
Means for extracting the reflected intensity of the laser reflected from the ground from the stored reflected intensity data;
The program according to claim 5, which functions as means for calculating the incident intensity of a laser immediately before reaching the tree crown by dividing the reflection intensity of the laser reflected by the ground by the reflectance of the ground. The reflection intensity data includes the intensity of the first pulse of the laser, the intensity of the last pulse of the laser, and the height of the reflection point of the laser;
The computer,
Means for storing the height of the ground;
Means for storing a specific threshold;
Means for comparing the threshold value with a difference value between the height of the reflection point and the height of the ground;
As a result of the comparison, means for extracting the intensity of the last pulse included in the reflection intensity data whose difference value is equal to or greater than the threshold as the reflection intensity of the laser that has passed through the crown,
As a result of the comparison, the intensity of the first pulse included in the reflection intensity data in which the difference value is smaller than the threshold is extracted as a reflection intensity of a laser reflected by the ground. 6. The program according to 6. The computer,
Any one of Claims 5-7 for making it function as a means to calculate the said reflectance by dividing the reflection intensity of the laser reflected in the said ground by the incident intensity of the laser immediately before reaching the crown of the said tree. The program described in the section. Laser irradiation means for irradiating laser from the sky to the ground,
Reflection intensity storage means for storing reflection intensity data of the laser reflected on the ground;
Reflectance storage means for storing the reflectance of the ground;
From the stored reflection intensity data, the reflection intensity of the laser reflected by the crown of the tree is extracted, and from the stored reflection intensity data, the reflection intensity of the laser reflected by the ground passing through the crown is extracted and stored. Based on the reflection intensity data and the reflectance, the incident intensity of the laser immediately before reaching the canopy is calculated, and the reflected intensity of the laser reflected from the ground is extracted from the stored reflection intensity data, and reaches the crown. The leaf area index of the tree is calculated using the incident intensity of the laser just before the reflection, the reflection intensity of the laser that has passed through the crown, the reflection intensity of the laser reflected by the crown, and the reflection intensity of the laser reflected by the ground. Computing means;
A leaf area index calculating device comprising: output means for outputting the leaf area index.