CN104778693A - Leaf area index calculation method based on projection algorithm and active contour model - Google Patents

Abstract

The invention discloses a leaf area index calculation method based on a projection algorithm and an active contour model. The method comprises steps as follows: firstly, tree point cloud data acquired through laser scanning are zoomed on an upper surface of a sphere in a certain proportion, and an upper spherical image is projected to a plane through stereographic projection and Lambert azimuth equivalent projection; secondly, a level set energy function is constructed, leaf pixels in a stand hemispherical picture are extracted, and according to stand three-dimensional laser point cloud data, topological structure features of point cloud are calculated and leaf point cloud is acquired in a classified manner in combination of a Gaussian mixture model; thirdly, stand canopy structure parameters are estimated from a divided two-dimensional hemispherical picture and classified three-dimensional point cloud data with a Miller formula method and an iterative inversion method; finally, the obtained measurement data are compared with leaf area indexes obtained through manual and actual measurement, and the accuracy and feasibility of the method are proved.

Description

Leaf area index based on projection algorithm and movable contour model calculates

Technical field

The present invention relates to a kind of leaf area index based on projection algorithm and movable contour model to calculate, particularly relate to a kind of broad leaf tree leaf area index method based on LTS (Transport Layer Security, laser scanning) cloud data.

Background technology

Broad leaf tree refers generally to the trees of dicotyledons, and have flat, broader blade, vein reticulates, and leaf is evergreen or fall leaves, and general blade face is broad, the leaf perennial woody plant having various shape with seeds difference.What have is evergreen, and fallen leaves class comes off mostly on autumn and winter season Ye Congzhi.The economic worth of broad leaf tree is large, and be much important commerical tree species, wherein some is famous and precious timber, landscape plant, various fruit etc., also has some broad leaf trees to be used as shade tree or flower garden afforestation tree.The present invention's research is that landscape plant is had a smile on one's face and flowering cherry blade face rebuilds and deformation.Data acquisition is by terrestrial Laser scanner.

Summary of the invention

Goal of the invention: in order to overcome the deficiencies in the prior art, the invention provides a kind of broad leaf tree real blade modeling towards laser point cloud data and deformation method, be intended to by laser scanner technology, build an accurate and feasible three-dimensional live standing tree data collection and analysis platform, merge the fresh approach of graph image, accurate forestry index is obtained by computer automatic analysis, thus the leaf area index that the live standing tree dynamic growth under accurate description Different forest stands condition changes.

Technical scheme: for achieving the above object, the technical solution used in the present invention is:

A kind of computing method of the leaf area index based on projection algorithm and movable contour model, it is characterized in that: the tree point cloud data first obtained laser scanning passes through the upper surface of certain proportional zoom at a ball, by stereographic projection and the equal area projection of Lambert position angle, upper spherical diagram picture is projected in plane again, secondly tectonic level energy collecting flow function, extract standing forest hemisphere figure Leaf pixel, and for standing forest three-dimensional laser point cloud data, the topological features of calculation level cloud also obtains leaf point cloud in conjunction with gauss hybrid models classification; Again, utilize Miller equation and iterative inversion method to neutralize sorted three dimensional point cloud estimation standing forest canopy structural parameter from the rear two-dimentional hemisphere figure of segmentation, finally will obtain measurement data and manually survey leaf area index comparing.

2, the computing method of the leaf area index based on projection algorithm and movable contour model according to claim 1, is characterized in that: comprise following cloud data and obtain and cloud data coordinate conversion step:

(1) scanning obtains whole tree cloud data, first to its conversion carrying out Cartesian coordinates and spherical co-ordinate by setting suitable unified r, can be on the sphere of r to a radius by point cloud compression at random.

$r=\sqrt{x^2+y^2+z^2}$

$\mathrm{\θ}=\arccos [z/\sqrt{(x^2+y^2+z^2)}]$

(2) to step (1) sphere data, through spherical co-ordinate and Cartesian coordinates conversion, the sphere compression process of a cloud is namely completed.

z＝r cosθ

(3) based on step (2) sphere information, Lambert position angle is utilized by spherical projection to plane, the theoretical following formula of projection process sectional drawing.If umbilical point coordinate is (x, y, z), the coordinate of the plane be converted to is (x ', y '), so transformational relation is:

$\left\{\begin{array}{c}{x}^{\′}=x\sqrt{\frac{2}{1-z}}\\ y^{\′}=y\sqrt{\frac{2}{1-z}}\end{array}\right.$

$(x,y.z)=(x^{\′}\sqrt{1-\frac{{\left(x^{\′}\right)}^2+{\left(y^{\′}\right)}^2}{4}},y^{\′}\sqrt{1-\frac{{\left(x^{\′}\right)}^2+{\left(y^{\′}\right)}^2}{4}},-1+\frac{\left(x^{\′}\right)+{\left(y^{\′}\right)}^2}{2})$

3, the computing method of the leaf area index based on projection algorithm and movable contour model according to claim 2, is characterized in that: comprise leaf area index step as follows:

(1) the true leaf area index based on Beer-Lambert law calculates

Digitizing hemisphere camera technique is the method for estimation Efficient leaf area the most conventional in reality.The main points of this method are: set suitable exposure rate, to distinguish blade face element (representing with 0) and sky portion (representing with 1).This differentiation normative reference be one based on fish eyes video camera shooting hemisphere photograph on pixel.In traditional implementation procedure, obtain porosity, zenith angle, the mathematical relation between pitch angle and leaf area index:

-ln P(θ)＝G(θ，α)×L/cos(θ)

What wherein P (θ) represented is porosity, and θ is the zenith angle of incident sunray, and α is Leaf inclination, and L is exactly leaf area index, if further modification, defines below carrying out:

T(θ)＝-ln P(θ)

$K(\mathrm{\θ},\mathrm{\α})=\frac{G(\mathrm{\θ},\mathrm{\α})}{\cos \left(\mathrm{\θ}\right)}$

Wherein G (θ, α) represents the averaging projection of region, blade face in the vertical direction of incident ray, and has to give a definition:

$G(\mathrm{\θ},\mathrm{\α})=\left\{\begin{array}{cc}\cos \mathrm{\α}\cos \mathrm{\θ},& \mathrm{\θ}\≤\mathrm{\π}/2-\mathrm{\α}\\ \cos \mathrm{\α}\cos \mathrm{\θ}(1+\frac{2(\tan x-x)}{\mathrm{\π}}),& \mathrm{\θ}>\mathrm{\π}/2-\mathrm{\α}\end{array}\right.$

Wherein x=cos ^-1(cot α cot θ), so T (θ)=-ln P (θ) formula can be rewritten into:

T(θ)＝K(θ，α)×L

K (θ, α) is extinction coefficient, if the inclined angle alpha of blade face element and porosity P (θ) known, the leaf area index L with some regions of identical zenith angle just can calculate by through type.In view of requiring this theory of identical zenith angle, being incorporated herein this idea of terrestrial latitude, being in the zenith angle approximately equal in Same Latitude region.

T (θ)=K (θ, what α) × L was widely used in Efficient leaf area (LAIe) asks in calculating, consider the covering problem of Main Differences at blade of true leaf area index and effective leaf area index, effectively cut layer operation at this to pre-service trees, each layer carries out Efficient leaf area respectively to be asked for and cumulative can obtain true leaf area index.For upper sphere, add the latitude line of the earth, intervening gaps 15 °, so obtain vertical view.It is the concentric circles of the spacing that nine are not waited, and each circle represents different zenith angles (i.e. trees height).

(2) Leaf angle inclination distribution function

Gap fraction is relevant with leaf area index and Leaf angle inclination distribution.The Leaf angle inclination distribution of frequently seen plants generally can be divided into 5 kinds: the distribution of horizontal leaf angle, the distribution of vertical leaf angle, the distribution of conical surface leaf angle, the distribution of sphere leaf angle, the distribution of ellipsoid leaf angle.Wherein, ellipsoid Leaf angle inclination distribution can regard as general type, can think special shape for other four kinds.

Ellipsoidal harmonics model formulation is:

$\left\{\begin{array}{c}p\left(\mathrm{\&alpha;}\right)=\frac{{2x}^3\sin \mathrm{\&alpha;}}{\mathrm{\&sigma;}{({\cos }^2\mathrm{\&alpha;}+x^2\sin \mathrm{\&alpha;})}^2}\\ \mathrm{\&sigma;}=x+\frac{\arcsin \left(\sqrt{1-x^2}\right)}{\sqrt{1-x^2}},x<1\\ \mathrm{\&sigma;}=x+\frac{\ln (1+\sqrt{1-x^{-2}})/(1-\sqrt{1-x^{-2}})}{2x\sqrt{1-x^{-2}}},x<1\end{array}\right.$

Wherein, p (α) is called Leaf inclination density function, represent that Leaf inclination is the ratio that the blade total area of α accounts for whole canopy area, x is the horizontal semiaxis of ellipsoid and the ratio of vertical semiaxis, x larger expression canopy Leaf inclination more convergence and horizontal distribution, the less expression Leaf inclination of x is more close to vertical distribution.

Have Leaf inclination density function:

In addition, x and average Leaf inclination close and are:

(3) projection function

After Leaf angle inclination distribution presses the simulation of ellipsoid distribution function, projection function formula is expanded into following form by Campbell (1990):

Wherein, ε ₁=(1-x ²) ^1/2, ε ₂=(1-x ^-2) ^1/2, G (θ) is averaging projection's area of θ direction blade, and x is horizontal, the vertical value length ratio of ellipsoid.As x < 1, ellipsoid is upright, and in canopy, most blade tilt is comparatively large, and the average Leaf inclination of canopy is greater than 57 °; As x=1, ellipsoid changes ball into, and this represents that blade occurs that the probability of various Leaf inclination is roughly equal, and the average Leaf inclination of canopy is about 57 °; As x > 1, ellipsoid lies low, and represent that in canopy, most blade tilt is less than normal, average Leaf inclination is less than 57 °.X < 1 is got in this research, obtains projected image

(4) Miller integral formula asks leaf area index

Based on Miller integral formula, Chen and Black supposes that gap fraction is only relevant with incident angle (visual angle), and use gap fraction data integrate of having derived asks the formula of LAI:

$\mathrm{LAI}=2{\&Integral;}_0^{\mathrm{\&pi;}/2}-\ln \left(T\left(\mathrm{\&theta;}\right)\right)\cos \mathrm{\&theta;}\sin \mathrm{\&theta;d\&theta;}$

For the ease of computer disposal, above formula integral formula is converted into difference formula:

$\mathrm{LAI}=2\underset{0}{\overset{n}{\mathrm{\&Sigma;}}}-\ln \left(T\left({\mathrm{\&theta;}}_i\right)\right)\cos {\mathrm{\&theta;}}_i\sin {\mathrm{\&theta;}}_i\mathrm{\&Delta;\&theta;}$

In formula, n represents concentric ring image being divided into n decile according to visual angle, T (θ _i) be the porosity data obtained in i-th ring, Δ θ is the angular field of view of each ring, Δ θ=pi/2 n.

(5) iterative inversion method asks leaf area index

After formula 8 and 9 is substituted into formula 5, arrange and obtain simultaneously containing leaf area index LAI and average Leaf inclination equation:

$\left\{\begin{array}{c}{T}_{\mathrm{sim}}\left(\mathrm{\&theta;}\right)=e^{\mathrm{\&tau;}}\\ \mathrm{\&tau;}=\frac{\mathrm{LAI}\sqrt{{\mathrm{\&tau;}}^2+{\tan }^2\mathrm{\&theta;}}}{\mathrm{\&tau;}+({\sin }^{-1}\sqrt{1-{\mathrm{\&tau;}}^2}/\sqrt{1-{\mathrm{\&tau;}}^2}}\\ \mathrm{\&tau;}=-3+{(\stackrel{\&OverBar;}{\mathrm{\&alpha;}}/9.65)}^{-0.6061}\end{array}\right.$

In formula, T _sim(θ) be the factor of porosity of simulation, LAI is leaf area index, average Leaf inclination, given a pair LAI and value just can calculate a porosity data T (θ) about view angle theta function curve.Iterative inversion method principle uses computing machine to carry out interative computation, by given Different L AI and combined value, the curve that simulation is obtained is close the most with observing the curve that obtains.This algorithm based on least square method, require obtain within the specific limits one group of most suitable LAI with combined value, make the factor of porosity analogue value T of its correspondence _sim(θ) minimum with the deviation extracting the factor of porosity T (θ) obtained from image, as shown in the formula

$L\stackrel{\&OverBar;}{A}I,\widehat{\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}=\arg \underset{A,\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}{\min }\underset{\mathrm{\&theta;}}{\mathrm{\&Sigma;}}{(T\left(\mathrm{\&theta;}\right)-T_{\mathrm{sim}}\left(\mathrm{\&theta;}\right))}^2$

In formula, optimum leaf area index LAI, be optimum average Leaf inclination, T (θ) is actual measurement gap fraction, T _sim(θ) be the factor of porosity of simulation, by formula expansion above.

4, the computing method of the leaf area index based on projection algorithm and movable contour model according to claim 3, is characterized in that: comprise test figure aligning step as follows:

(1) estimation of clumping index

An important hypothesis is had: canopy leaves meets Poisson distribution in space distribution by the method for hemisphere Photographic technique estimation leaf area index (LAI) and average Leaf inclination.In Poisson distribution, assuming that canopy leaves size is unified, spatially free stochastic distribution, but in fact, because vegetation structure is complicated, canopy leaves can not free stochastic distribution, therefore there is deviation by the LAI of the method acquisitions such as hemisphere Photographic technique and actual LAI.In order to address this problem, Nilson proposes the correction model of Poisson distribution, adds a parameter Ω in fact exactly and be used for correcting this deviation on basis, as shown in the formula:

T(θ)＝e ^{-Ω□LAI□G(θ，α)/cosθ}

Researcher further provides new variable " effective leaf area index " LAI _e, it equals parameter Ω and is multiplied by actual leaf area index (LAI), determines that Ω is clumping index (Clumping Index, CI) simultaneously, is actually LAI above with the LAI that hemisphere Photographic technique obtains _e, as shown in the formula:

LAI _e＝LAI□Ω

Clumping index computing method have 3 classes: the distributed intelligence based on canopy pore size, the logarithmic mean based on gap fraction and above-mentioned two class integrated approachs.This research uses the method for Lang (1986), by asking factor of porosity logarithmic mean to obtain clumping index CI, sees following formula:

In formula, Ω (θ) is the clumping index about visual angle, and clumping index changes along with visual angle change.

(2) data statistic

Table 1LAI and average Leaf inclination result data

The average Leaf inclination of table 2

Accompanying drawing explanation

Accompanying drawing 1 three-dimensional coordinate and spherical coordinates corresponding relation

Accompanying drawing 2 projection function image

Accompanying drawing 3 hemisphere Photographic technique lab diagram

Accompanying drawing 4 hemisphere Photographic technique factor of porosity and visual angle functional arrangement (the different n value of a, six groups of data images during b n=9)

Accompanying drawing 5 six groups of laser point cloud data Leaf angle inclination distribution figure (a blade-section b limb part)

Accompanying drawing 6 Leaf inclination density function curve

Accompanying drawing 7 six groups of cloud data LAI and visual angle functional arrangement

Embodiment

2, the computing method of the leaf area index based on projection algorithm and movable contour model according to claim 1, is characterized in that: comprise following cloud data and obtain and cloud data coordinate conversion step:

(1) scanning obtains whole tree cloud data, first to its conversion carrying out Cartesian coordinates and spherical co-ordinate by setting suitable unified r, can be on the sphere of r to a radius by point cloud compression at random.

$r=\sqrt{x^2+y^2+z^2}$

$\mathrm{\&theta;}=\arccos [z/\sqrt{(x^2+y^2+z^2)}]$

(2) to step (1) sphere data, through spherical co-ordinate and Cartesian coordinates conversion, the sphere compression process of a cloud is namely completed.

z＝r cosθ

(3) based on step (2) sphere information, Lambert position angle is utilized by spherical projection to plane, the theoretical following formula of projection process sectional drawing.If umbilical point coordinate is (x, y, z), the coordinate of the plane be converted to is (x ', y '), so transformational relation is:

$\left\{\begin{array}{c}{x}^{\&prime;}=x\sqrt{\frac{2}{1-z}}\\ y^{\&prime;}=y\sqrt{\frac{2}{1-z}}\end{array}\right.$

$(x,y.z)=(x^{\&prime;}\sqrt{1-\frac{{\left(x^{\&prime;}\right)}^2+{\left(y^{\&prime;}\right)}^2}{4}},y^{\&prime;}\sqrt{1-\frac{{\left(x^{\&prime;}\right)}^2+{\left(y^{\&prime;}\right)}^2}{4}},-1+\frac{\left(x^{\&prime;}\right)+{\left(y^{\&prime;}\right)}^2}{2})$

3, the computing method of the leaf area index based on projection algorithm and movable contour model according to claim 2, is characterized in that: comprise leaf area index step as follows:

(1) the true leaf area index based on Beer-Lambert law calculates

Digitizing hemisphere camera technique is the method for estimation Efficient leaf area the most conventional in reality.The main points of this method are: set suitable exposure rate, to distinguish blade face element (representing with 0) and sky portion (representing with 1).This differentiation normative reference be one based on fish eyes video camera shooting hemisphere photograph on pixel.In traditional implementation procedure, obtain porosity, zenith angle, the mathematical relation between pitch angle and leaf area index:

-ln P(θ)＝G(θ，α)×L/cos(θ)

What wherein P (θ) represented is porosity, and θ is the zenith angle of incident sunray, and α is Leaf inclination, and L is exactly leaf area index, if further modification, defines below carrying out:

T(θ)＝-ln P(θ)

$K(\mathrm{\&theta;},\mathrm{\&alpha;})=\frac{G(\mathrm{\&theta;},\mathrm{\&alpha;})}{\cos \left(\mathrm{\&theta;}\right)}$

Wherein G (θ, α) represents the averaging projection of region, blade face in the vertical direction of incident ray, and has to give a definition:

$G(\mathrm{\&theta;},\mathrm{\&alpha;})=\left\{\begin{array}{cc}\cos \mathrm{\&alpha;}\cos \mathrm{\&theta;},& \mathrm{\&theta;}\&le;\mathrm{\&pi;}/2-\mathrm{\&alpha;}\\ \cos \mathrm{\&alpha;}\cos \mathrm{\&theta;}(1+\frac{2(\tan x-x)}{\mathrm{\&pi;}}),& \mathrm{\&theta;}>\mathrm{\&pi;}/2-\mathrm{\&alpha;}\end{array}\right.$

Wherein x=cos ^-1(cot α cot θ), so T (θ)=-ln P (θ) formula can be rewritten into:

T(θ)＝K(θ，α)×L

K (θ, α) is extinction coefficient, if the inclined angle alpha of blade face element and porosity P (θ) known, the leaf area index L with some regions of identical zenith angle just can calculate by through type.In view of requiring this theory of identical zenith angle, being incorporated herein this idea of terrestrial latitude, being in the zenith angle approximately equal in Same Latitude region.

T (θ)=K (θ, what α) × L was widely used in Efficient leaf area (LAIe) asks in calculating, consider the covering problem of Main Differences at blade of true leaf area index and effective leaf area index, effectively cut layer operation at this to pre-service trees, each layer carries out Efficient leaf area respectively to be asked for and cumulative can obtain true leaf area index.For upper sphere, add the latitude line of the earth, intervening gaps 15.So, obtain vertical view.It is the concentric circles of the spacing that nine are not waited, and each circle represents different zenith angles (i.e. trees height).

(2) Leaf angle inclination distribution function

Gap fraction is relevant with leaf area index and Leaf angle inclination distribution.The Leaf angle inclination distribution of frequently seen plants generally can be divided into 5 kinds: the distribution of horizontal leaf angle, the distribution of vertical leaf angle, the distribution of conical surface leaf angle, the distribution of sphere leaf angle, the distribution of ellipsoid leaf angle.Wherein, ellipsoid Leaf angle inclination distribution can regard as general type, can think special shape for other four kinds.

Ellipsoidal harmonics model formulation is:

$\left\{\begin{array}{c}p\left(\mathrm{\&alpha;}\right)=\frac{{2x}^3\sin \mathrm{\&alpha;}}{\mathrm{\&sigma;}{({\cos }^2\mathrm{\&alpha;}+x^2\sin \mathrm{\&alpha;})}^2}\\ \mathrm{\&sigma;}=x+\frac{\arcsin \left(\sqrt{1-x^2}\right)}{\sqrt{1-x^2}},x<1\\ \mathrm{\&sigma;}=x+\frac{\ln (1+\sqrt{1-x^{-2}})/(1-\sqrt{1-x^{-2}})}{2x\sqrt{1-x^{-2}}},x<1\end{array}\right.$

Wherein, p (α) is called Leaf inclination density function, represent that Leaf inclination is the ratio that the blade total area of α accounts for whole canopy area, x is the horizontal semiaxis of ellipsoid and the ratio of vertical semiaxis, x larger expression canopy Leaf inclination more convergence and horizontal distribution, the less expression Leaf inclination of x is more close to vertical distribution.

Have Leaf inclination density function:

In addition, x and average Leaf inclination close and are:

(3) projection function

After Leaf angle inclination distribution presses the simulation of ellipsoid distribution function, projection function formula is expanded into following form by Campbell (1990):

Wherein, ε ₁=(1-x ²) ^1/2, ε ₂=(1-x ^-2) ^1/2, G (θ) is averaging projection's area of θ direction blade, and x is horizontal, the vertical value length ratio of ellipsoid.As x < 1, ellipsoid is upright, and in canopy, most blade tilt is comparatively large, and the average Leaf inclination of canopy is greater than 57 °; As x=1, ellipsoid changes ball into, and this represents that blade occurs that the probability of various Leaf inclination is roughly equal, and the average Leaf inclination of canopy is about 57 °; As x > 1, ellipsoid lies low, and represent that in canopy, most blade tilt is less than normal, average Leaf inclination is less than 57 °.X < 1 is got in this research, obtains projected image

(4) Miller integral formula asks leaf area index

Based on Miller integral formula, Chen and Black supposes that gap fraction is only relevant with incident angle (visual angle), and use gap fraction data integrate of having derived asks the formula of LAI:

$\mathrm{LAI}=2{\&Integral;}_0^{\mathrm{\&pi;}/2}-\ln \left(T\left(\mathrm{\&theta;}\right)\right)\cos \mathrm{\&theta;}\sin \mathrm{\&theta;d\&theta;}$

For the ease of computer disposal, above formula integral formula is converted into difference formula:

$\mathrm{LAI}=2\underset{0}{\overset{n}{\mathrm{\&Sigma;}}}-\ln \left(T\left({\mathrm{\&theta;}}_i\right)\right)\cos {\mathrm{\&theta;}}_i\sin {\mathrm{\&theta;}}_i\mathrm{\&Delta;\&theta;}$

In formula, n represents concentric ring image being divided into n decile according to visual angle, T (θ _i) be the porosity data obtained in i-th ring, Δ θ is the angular field of view of each ring, Δ θ=pi/2 n.

(5) iterative inversion method asks leaf area index

After formula 8 and 9 is substituted into formula 5, arrange and obtain simultaneously containing leaf area index LAI and average Leaf inclination equation:

$\left\{\begin{array}{c}{T}_{\mathrm{sim}}\left(\mathrm{\&theta;}\right)=e^{\mathrm{\&tau;}}\\ \mathrm{\&tau;}=\frac{\mathrm{LAI}\sqrt{{\mathrm{\&tau;}}^2+{\tan }^2\mathrm{\&theta;}}}{\mathrm{\&tau;}+({\sin }^{-1}\sqrt{1-{\mathrm{\&tau;}}^2}/\sqrt{1-{\mathrm{\&tau;}}^2}}\\ \mathrm{\&tau;}=-3+{(\stackrel{\&OverBar;}{\mathrm{\&alpha;}}/9.65)}^{-0.6061}\end{array}\right.$

In formula, T _sim(θ) be the factor of porosity of simulation, LAI is leaf area index, average Leaf inclination, given a pair LAI and value just can calculate a porosity data T (θ) about view angle theta function curve.Iterative inversion method principle uses computing machine to carry out interative computation, by given Different L AI and combined value, the curve that simulation is obtained is close the most with observing the curve that obtains.This algorithm based on least square method, require obtain within the specific limits one group of most suitable LAI with combined value, make the factor of porosity analogue value T of its correspondence _sim(θ) minimum with the deviation extracting the factor of porosity T (θ) obtained from image, as shown in the formula

$L\stackrel{\&OverBar;}{A}I,\widehat{\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}=\arg \underset{A,\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}{\min }\underset{\mathrm{\&theta;}}{\mathrm{\&Sigma;}}{(T\left(\mathrm{\&theta;}\right)-T_{\mathrm{sim}}\left(\mathrm{\&theta;}\right))}^2$

In formula, optimum leaf area index LAI, be optimum average Leaf inclination, T (θ) is actual measurement gap fraction, T _sim(θ) be the factor of porosity of simulation, by formula expansion above.

4, the computing method of the leaf area index based on projection algorithm and movable contour model according to claim 3, is characterized in that: comprise test figure aligning step as follows:

(1) estimation of clumping index

An important hypothesis is had: canopy leaves meets Poisson distribution in space distribution by the method for hemisphere Photographic technique estimation leaf area index (LAI) and average Leaf inclination.In Poisson distribution, assuming that canopy leaves size is unified, spatially free stochastic distribution, but in fact, because vegetation structure is complicated, canopy leaves can not free stochastic distribution, therefore there is deviation by the LAI of the method acquisitions such as hemisphere Photographic technique and actual LAI.In order to address this problem, Nilson proposes the correction model of Poisson distribution, adds a parameter Ω in fact exactly and be used for correcting this deviation on basis, as shown in the formula:

T(θ)＝e ^{-Ω□LAI□G(θ，α)/cosθ}

Researcher further provides new variable " effective leaf area index " LAI _e, it equals parameter Ω and is multiplied by actual leaf area index (LAI), determines that Ω is clumping index (Clumping Index, CI) simultaneously, is actually LAI above with the LAI that hemisphere Photographic technique obtains _e, as shown in the formula:

LAI _e＝LAI□Ω

Clumping index computing method have 3 classes: the distributed intelligence based on canopy pore size, the logarithmic mean based on gap fraction and above-mentioned two class integrated approachs.This research uses the method for Lang (1986), by asking factor of porosity logarithmic mean to obtain clumping index CI, sees following formula:

In formula, Ω (θ) is the clumping index about visual angle, and clumping index changes along with visual angle change.

(2) data statistic

Table 1LAI and average Leaf inclination result data

The average Leaf inclination of table 2

Claims (4)

1. the computing method based on the leaf area index of projection algorithm and movable contour model, it is characterized in that: the tree point cloud data first obtained laser scanning passes through the upper surface of certain proportional zoom at a ball, by stereographic projection and the equal area projection of Lambert position angle, upper spherical diagram picture is projected in plane again, secondly tectonic level energy collecting flow function, extract standing forest hemisphere figure Leaf pixel, and for standing forest three-dimensional laser point cloud data, the topological features of calculation level cloud also obtains leaf point cloud in conjunction with gauss hybrid models classification; Again, utilize Miller equation and iterative inversion method to neutralize sorted three dimensional point cloud estimation standing forest canopy structural parameter from the rear two-dimentional hemisphere figure of segmentation, finally will obtain measurement data and manually survey leaf area index comparing.

2. the computing method of the leaf area index based on projection algorithm and movable contour model according to claim 1, is characterized in that: comprise following cloud data and obtain and cloud data coordinate conversion step:

(1) scanning obtains whole tree cloud data, first to its conversion carrying out Cartesian coordinates and spherical co-ordinate by setting suitable unified r, can be on the sphere of r to a radius by point cloud compression at random.

$r=\sqrt{x^2+y^2+z^2}$

$\mathrm{\&theta;}=\arccos [z/\sqrt{(x^2+y^2+z^2)}]$

(2) to step (1) sphere data, through spherical co-ordinate and Cartesian coordinates conversion, the sphere compression process of a cloud is namely completed.

z＝rcosθ

(3) based on step (2) sphere information, Lambert position angle is utilized by spherical projection to plane, the theoretical following formula of projection process sectional drawing.If umbilical point coordinate is (x, y, z), the coordinate of the plane be converted to is (x ', y '), so transformational relation is:

$\left\{\begin{array}{c}{x}^{\&prime;}=x\sqrt{\frac{2}{1-z}}\\ y^{\&prime;}=y\sqrt{\frac{2}{1-z}}\end{array}\right.$

$(x,y.z)=(x^{\&prime;}\sqrt{1-\frac{{\left(x^{\&prime;}\right)}^2+{\left(y^{\&prime;}\right)}^2}{4}},y^{\&prime;}\sqrt{1-\frac{{\left(x^{\&prime;}\right)}^2+{\left(y^{\&prime;}\right)}^2}{4}},-1+\frac{\left(x^{\&prime;}\right)+{\left(y^{\&prime;}\right)}^2}{2})$

3. the computing method of the leaf area index based on projection algorithm and movable contour model according to claim 2, is characterized in that: comprise leaf area index step as follows:

(1) the true leaf area index based on Beer-Lambert law calculates

Digitizing hemisphere camera technique is the method for estimation Efficient leaf area the most conventional in reality.The main points of this method are: set suitable exposure rate, to distinguish blade face element (representing with 0) and sky portion (representing with 1).This differentiation normative reference be one based on fish eyes video camera shooting hemisphere photograph on pixel.In traditional implementation procedure, obtain porosity, zenith angle, the mathematical relation between pitch angle and leaf area index:

-ln P(θ)＝G(θ，α)×L/cos(θ)

What wherein P (θ) represented is porosity, and θ is the zenith angle of incident sunray, and α is Leaf inclination, and L is exactly leaf area index, if further modification, defines below carrying out:

T(θ)＝-lnP(θ)

$K(\mathrm{\&theta;},\mathrm{\&alpha;})=\frac{G(\mathrm{\&theta;},\mathrm{\&alpha;})}{\cos \left(\mathrm{\&theta;}\right)}$

Wherein G (θ, α) represents the averaging projection of region, blade face in the vertical direction of incident ray, and has to give a definition:

$G(\mathrm{\&theta;},\mathrm{\&alpha;})=\left\{\begin{array}{cc}\cos \mathrm{\&alpha;}\cos \mathrm{\&theta;},& \mathrm{\&theta;}\&le;\mathrm{\&pi;}/2-\mathrm{\&alpha;}\\ \cos \mathrm{\&alpha;}\cos \mathrm{\&theta;}(1+\frac{2(\tan x-x)}{\mathrm{\&pi;}}),& \mathrm{\&theta;}>\mathrm{\&pi;}/2-\mathrm{\&alpha;}\end{array}\right.$

Wherein x=cos ^-1(cot α cot θ), so T (θ)=-lnP (θ) formula can be rewritten into:

T(θ)＝K(θ，α)×L

K (θ, α) is extinction coefficient, if the inclined angle alpha of blade face element and porosity P (θ) known, the leaf area index L with some regions of identical zenith angle just can calculate by through type.In view of requiring this theory of identical zenith angle, being incorporated herein this idea of terrestrial latitude, being in the zenith angle approximately equal in Same Latitude region.

T (θ)=K (θ, what α) × L was widely used in Efficient leaf area (LAIe) asks in calculating, consider the covering problem of Main Differences at blade of true leaf area index and effective leaf area index, effectively cut layer operation at this to pre-service trees, each layer carries out Efficient leaf area respectively to be asked for and cumulative can obtain true leaf area index.For upper sphere, add the latitude line of the earth, intervening gaps 15 °, so obtain vertical view.It is the concentric circles of the spacing that nine are not waited, and each circle represents different zenith angles (i.e. trees height).

(2) Leaf angle inclination distribution function

Gap fraction is relevant with leaf area index and Leaf angle inclination distribution.The Leaf angle inclination distribution of frequently seen plants generally can be divided into 5 kinds: the distribution of horizontal leaf angle, the distribution of vertical leaf angle, the distribution of conical surface leaf angle, the distribution of sphere leaf angle, the distribution of ellipsoid leaf angle.Wherein, ellipsoid Leaf angle inclination distribution can regard as general type, can think special shape for other four kinds.

Ellipsoidal harmonics model formulation is:

$\left\{\begin{array}{c}p\left(\mathrm{\&alpha;}\right)=\frac{2x^3\sin \mathrm{\&alpha;}}{\mathrm{\&sigma;}{({\cos }^2\mathrm{\&alpha;}+x^2\sin \mathrm{\&alpha;})}^2}\\ \mathrm{\&sigma;}=x+\frac{\arcsin \left(\sqrt{1-x^2}\right)}{\sqrt{1-x^2}},x<1\\ \mathrm{\&sigma;}=x+\frac{\ln (1+\sqrt{1-x^{-2}})/(1-\sqrt{1-x^{-2}})}{2x\sqrt{1-x^{-2}}},x<1\end{array}\right.$

Wherein, p (α) is called Leaf inclination density function, represent that Leaf inclination is the ratio that the blade total area of α accounts for whole canopy area, x is the horizontal semiaxis of ellipsoid and the ratio of vertical semiaxis, x larger expression canopy Leaf inclination more convergence and horizontal distribution, the less expression Leaf inclination of x is more close to vertical distribution.

Have Leaf inclination density function:

In addition, x and average Leaf inclination close and are:

(3) projection function

After Leaf angle inclination distribution presses the simulation of ellipsoid distribution function, projection function formula is expanded into following form by Campbell (1990):

$(\left\{\begin{array}{c}G\left(\mathrm{\&theta;}\right)=\frac{{(x^2+{\tan }^2\mathrm{\&theta;})}^{1/2}\cos \mathrm{\&theta;}}{x+\left({\sin }^{-1}{\mathrm{\&epsiv;}}_1\right)/{\mathrm{\&epsiv;}}_1},x\&le;1\\ G\left(\mathrm{\&theta;}\right)=\frac{{(x^2+{\tan }^2\mathrm{\&theta;})}^{1/2}\cos \mathrm{\&theta;}}{x+\frac{1}{2{\mathrm{\&epsiv;}}_{2x}}\ln [(1+{\mathrm{\&epsiv;}}_2)/(1-{\mathrm{\&epsiv;}}_2)]},x>1\end{array}\right.$

Wherein, ε ₁=(1-x ²) ^1/2, ε ₂=(1-x ^-2) ^1/2, G (θ) is averaging projection's area of θ direction blade, and x is horizontal, the vertical value length ratio of ellipsoid.As x < 1, ellipsoid is upright, and in canopy, most blade tilt is comparatively large, and the average Leaf inclination of canopy is greater than 57 °; As x=1, ellipsoid changes ball into, and this represents that blade occurs that the probability of various Leaf inclination is roughly equal, and the average Leaf inclination of canopy is about 57 °; As x > 1, ellipsoid lies low, and represent that in canopy, most blade tilt is less than normal, average Leaf inclination is less than 57 °.X < 1 is got in this research, obtains projected image

(4) Miller integral formula asks leaf area index

Based on Miller integral formula, Chen and Black supposes that gap fraction is only relevant with incident angle (visual angle), and use gap fraction data integrate of having derived asks the formula of LAI:

$\mathrm{LAI}=2{\&Integral;}_0^{\mathrm{\&pi;}/2}-\ln \left(T\left(\mathrm{\&theta;}\right)\right)\cos \mathrm{\&theta;}\sin \mathrm{\&theta;d\&theta;}$

For the ease of computer disposal, above formula integral formula is converted into difference formula:

$\mathrm{LAI}=2\underset{0}{\overset{n}{\mathrm{\&Sigma;}}}-\ln T\left(T\left({\mathrm{\&theta;}}_i\right)\right)\cos {\mathrm{\&theta;}}_i\sin {\mathrm{\&theta;}}_i\mathrm{\&Delta;\&theta;}$

In formula, n represents concentric ring image being divided into n decile according to visual angle, T (θ _i) be the porosity data obtained in i-th ring, Δ θ is the angular field of view of each ring, Δ θ=pi/2 n.

(5) iterative inversion method asks leaf area index

After formula 8 and 9 is substituted into formula 5, arrange and obtain simultaneously containing leaf area index LAI and average Leaf inclination equation:

$\left\{\begin{array}{c}{T}_{\mathrm{sim}}=\left(\mathrm{\&theta;}\right)e^{\mathrm{\&tau;}}\\ \mathrm{\&tau;}=\frac{\mathrm{LAI}\sqrt{{\mathrm{\&iota;}}^2+{\tan }^2\mathrm{\&theta;}}}{\mathrm{\&iota;}+({\sin }^{-1}\sqrt{1-{\mathrm{\&iota;}}^2}/\sqrt{1-{\mathrm{\&iota;}}^2}}\\ \mathrm{\&iota;}=-3+{(\stackrel{\&OverBar;}{\mathrm{\&alpha;}}/9.65)}^{-0.6061}\end{array}\right.$

In formula, T _sim(θ) be the factor of porosity of simulation, LAI is leaf area index, average Leaf inclination, given a pair LAI and value just can calculate a porosity data T (θ) about view angle theta function curve.Iterative inversion method principle uses computing machine to carry out interative computation, by given Different L AI and combined value, the curve that simulation is obtained is close the most with observing the curve that obtains.This algorithm based on least square method, require obtain within the specific limits one group of most suitable LAI with combined value, make the factor of porosity analogue value T of its correspondence _sim(θ) minimum with the deviation extracting the factor of porosity T (θ) obtained from image, as shown in the formula

$L\hat{A}I,\widehat{\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}=\arg \underset{A,\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}{\min }\underset{\mathrm{\&theta;}}{\mathrm{\&Sigma;}}{(T\left(\mathrm{\&theta;}\right)-T_{\mathrm{sim}}\left(\mathrm{\&theta;}\right))}^2$

In formula, optimum leaf area index LAI, be optimum average Leaf inclination, T (θ) is actual measurement gap fraction, T _sim(θ) be the factor of porosity of simulation, by formula expansion above.

4. the computing method of the leaf area index based on projection algorithm and movable contour model according to claim 3, is characterized in that: comprise test figure aligning step as follows:

An important hypothesis is had: canopy leaves meets Poisson distribution in space distribution by the method for hemisphere Photographic technique estimation leaf area index (LAI) and average Leaf inclination.In Poisson distribution, assuming that canopy leaves size is unified, spatially free stochastic distribution, but in fact, because vegetation structure is complicated, canopy leaves can not free stochastic distribution, therefore there is deviation by the LAI of the method acquisitions such as hemisphere Photographic technique and actual LAI.In order to address this problem, Nilson proposes the correction model of Poisson distribution, adds a parameter Ω in fact exactly and be used for correcting this deviation on basis, as shown in the formula:

T(θ)＝e ^{-Ω□LAI□G}( ^{θ，α)/cosθ}

Researcher further provides new variable " effective leaf area index " LAI _e, it equals parameter Ω and is multiplied by actual leaf area index (LAI), determines that Ω is clumping index (Clumping Index, CI) simultaneously, is actually LAI above with the LAI that hemisphere Photographic technique obtains _e, as shown in the formula:

LAI _e＝LAI□Ω

Clumping index computing method have 3 classes: the distributed intelligence based on canopy pore size, the logarithmic mean based on gap fraction and above-mentioned two class integrated approachs.This research uses the method for Lang (1986), by asking factor of porosity logarithmic mean to obtain clumping index CI, sees following formula:

In formula, Ω (θ) is the clumping index about visual angle, and clumping index changes along with visual angle change.