CN112577907A - Urban green land tree crown loss rate calculation method - Google Patents

Urban green land tree crown loss rate calculation method Download PDF

Info

Publication number
CN112577907A
CN112577907A CN202011295757.XA CN202011295757A CN112577907A CN 112577907 A CN112577907 A CN 112577907A CN 202011295757 A CN202011295757 A CN 202011295757A CN 112577907 A CN112577907 A CN 112577907A
Authority
CN
China
Prior art keywords
image
crown
spectrum
data
green land
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011295757.XA
Other languages
Chinese (zh)
Other versions
CN112577907B (en
Inventor
张浪
易扬
张桂莲
林勇
邢璐琪
仲启铖
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Academy of Landscape Architecture Science and Planning
Original Assignee
Shanghai Academy of Landscape Architecture Science and Planning
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Academy of Landscape Architecture Science and Planning filed Critical Shanghai Academy of Landscape Architecture Science and Planning
Priority to CN202011295757.XA priority Critical patent/CN112577907B/en
Publication of CN112577907A publication Critical patent/CN112577907A/en
Application granted granted Critical
Publication of CN112577907B publication Critical patent/CN112577907B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/84—Systems specially adapted for particular applications
    • G01N21/88—Investigating the presence of flaws or contamination
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/29—Geographical information databases
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/13—Edge detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/136—Segmentation; Edge detection involving thresholding
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/84—Systems specially adapted for particular applications
    • G01N2021/8466—Investigation of vegetal material, e.g. leaves, plants, fruits
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G06T2207/10036—Multispectral image; Hyperspectral image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G06T2207/10041—Panchromatic image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation
    • G06T2207/30188—Vegetation; Agriculture

Abstract

A method for calculating the crown defect rate of trees in urban green land is characterized in that a hyperspectral image, a high-resolution panchromatic remote sensing image and measured data are used as data sources, a set of method for calculating the crown defect rate of the trees in the urban green land is established on the premise of identifying the contour of a vegetation crown, extracting the spectral end member spectral curve characteristics and identifying species based on spectral angle classification, and the method mainly comprises the steps of investigating the type of land features, taking images, preprocessing the images, segmenting a threshold value, identifying the contour, extracting end member spectra, identifying the species of the plants and establishing and calculating a crown defect rate model.

Description

Urban green land tree crown loss rate calculation method
Technical Field
The invention relates to the field of urban green land area vegetation measurement and calculation, in particular to a method for calculating the crown defect rate of urban green land trees.
Background
The measurement and calculation of the crown defect rate of the urban green land vegetation has important significance for monitoring urban green land resources and evaluating quality. The method can be used for mastering the growth conditions and health conditions of plants (arbors and shrubs) in urban greenbelts, and further has important significance for pest control, sustainable management and maintenance of plant diversity of the urban greenbelts. For example, the local standard 'urban forest health economic management code' DB 34/T2199-. And the crown loss rate in the ancient tree health evaluation research progress published in the world forestry research of Liu Yu and Xucheng 2013 is also used as one of the ancient tree health evaluation indexes. The crown defect rate is an index which appears in the tree health degree and structure evaluation index system at high frequency. However, the measurement and calculation of the crown defect rate of urban green land vegetation mainly depends on field actual measurement or a statistical analysis method based on the actual measurement method, and the specific method is to estimate the crown defect condition of each plant (arbor and shrub) through manual investigation and quantify the crown defect condition in a percentage mode. Such observation and calculation methods have strong subjectivity, low precision and huge workload, and are difficult to be practically popularized on the 'surface'.
Disclosure of Invention
The invention aims to solve the technical problem of providing a method for calculating the crown defect rate of trees in urban green lands, which is a crown defect calculation model designed on the premise of plant outline and species identification and used for calculating the crown defect rate of each plant, can quickly and accurately measure the area and the defect rate of the crowns of urban green land vegetation in a large area range, and has important significance for mastering the growth condition and the health condition of plants and further for sustainable operation management and maintenance of urban green lands and plant diversity.
In order to solve the technical problems, the technical scheme adopted by the invention is as follows:
the urban green land vegetation crown defect rate calculation method comprises the following steps:
acquiring panchromatic images and hyperspectral image data of a target area, and preprocessing the data;
secondly, segmenting the hyperspectral image, and distinguishing vegetation and non-vegetation areas in the image data to form a vegetation image data set of the urban green land, wherein the vegetation image data set comprises all trees, shrubs and grasslands in the urban green land;
thirdly, carrying out edge detection on the image of the plant area by utilizing the full-color image, and identifying the outline of the plant species;
combining visible light, red-edge wave bands and near-infrared wave bands of the hyperspectral images, positioning, identifying and gathering the purest pixels in data, acquiring pure end members in the images of the urban green land plant areas, and extracting end member spectrum characteristic curves of different plant species and storing the end members in a warehouse;
surveying and recording typical plant types and geographical positions of urban greenhouses, wherein the plant types comprise different kinds of trees, shrubs and grasslands;
classifying all data in the vegetation image according to the end member wave spectrum characteristic curves of different plant species in the hyperspectral image obtained in the step three, and then overlapping the data with the contour identification result of the panchromatic image to define the species category of each contour;
and step six, establishing a crown damage rate model according to the contour and species identification result, and calculating the crown damage rate of the trees and shrubs.
The concrete process of the step six is as follows:
defining a polygon m at the outermost edge of the projected outline of the crown as
Figure BDA0002785258110000025
The polygon hasPnA vertex Pn(n-1, 2, …, m) is arranged along the positive direction of the boundary, and the coordinates are (x) in sequence1，y1)，(x2，y2)，…，(xn，yn) Establishing
Figure BDA0002785258110000026
The centroid of the polygon m and any two adjacent vertexes form a triangle, the area of the triangle is obtained by the outer product of two plane vectors formed by the three vertexes, and the area of the polygon m at the outermost edge of the crown projection
Figure BDA0002785258110000021
Comprises the following steps:
Figure BDA0002785258110000022
in the polygon at the outermost edge of the crown projection, the exposed part without the blade shading is formed by k (k is 1,2, …, q) polygons ZkComposition of each polygon having OiOne vertex, vertex Oi(i-1, 2, …, j) is arranged along the positive direction of the boundary, and the coordinates are (a) in sequence1，b1)，(a2，b2)，…，(ai，bi) Establishing ZkPolygon area vector diagram, polygon ZkThe centroid of the tree and any two adjacent vertexes form a triangle, the area of the triangle is obtained by the outer product of two plane vectors formed by the three vertexes, and the area S of the exposed part without blade shielding in the polygon of the outermost edge of the crown projectionZComprises the following steps:
Figure BDA0002785258110000023
Figure BDA0002785258110000024
plant and method for producing the sameIdeal crown projection area SQComprises the following steps:
Figure BDA0002785258110000031
the crown defect rate P of the plants is as follows:
Figure BDA0002785258110000032
the preprocessing of the data in the first step includes: and respectively carrying out radiometric calibration, atmospheric correction, geometric correction, mosaic and splicing on the panchromatic image and the hyperspectral image.
The specific process of image segmentation in the second step is to segment the image by using the normalized vegetation index NDVI and the Otsu algorithm.
First, NDVI is calculated, with the formula:
NDVI＝(ρNIR-ρRED)/(ρNIR+ρRED)
in the formula, ρNIRReflectivity in the near red band, pREDReflectance for the red band;
secondly, image segmentation is carried out by utilizing an Otsu algorithm, and the steps are as follows:
let the gray level of the NDVI image be L (G ═ 1,2,3, …, L), and the number of pixels at gray level i be niThe total pixel number of NDVI is shown as follows:
Figure BDA0002785258110000033
the histogram is normalized to:
Figure BDA0002785258110000034
Figure BDA0002785258110000035
wherein p (i) represents the probability of the appearance of the pixel with the gray level i in the image;
and (3) expressing a threshold value by T, and dividing pixels in the normalized image into two types according to the gray level T: c1T and C are not more than2T, the probability of occurrence of two classes is:
Figure BDA0002785258110000036
Figure BDA0002785258110000037
suppose that the mean values of the two types of pixels are m1And m2The average value of global pixels of the NDVI image is mgThen:
w1*m1+w2*m2＝mg
w1+w2＝1
the variance is now:
δ2＝p1(m1-mg)2+p2(m2-mg)2
namely:
δ2＝p1p2(m1-m2)2
in the formula, delta2The inter-class variance is the maximum value, and the corresponding threshold T is the optimal segmentation threshold of the NDVI image.
As shown in fig. 3, the specific process in the third step is as follows:
the method is characterized in that edge detection is carried out based on a Sobel operator to achieve the effect of identifying the canopy contour, and the method comprises the following steps:
firstly, setting a Sobel operator, wherein the formula is as follows:
Figure BDA0002785258110000041
Figure BDA0002785258110000042
secondly, performing plane convolution on a pixel set Z of the panchromatic high-resolution image and A and B to respectively calculate gray values Gx and Gy of the image detected by the transverse edge and the longitudinal edge, wherein the formula is as follows:
Figure BDA0002785258110000043
Figure BDA0002785258110000044
finally, the horizontal and vertical gray scale of each pixel of the panchromatic high-resolution image can be calculated, and the formula is as follows:
Figure BDA0002785258110000045
using an approximation:
|f(Gx，Gy)|＝|Gx|+|Gy|
setting the threshold value as h and the pixel value f (G)x，Gy) If the sum of the edge points is larger than h, the edge points are marked as edge points, and the set of the edge points is the result of contour recognition; then, on the basis of full-color image contour identification, referring to the type and the geographic position of a sample plot plant species, taking a crown high-spectrum image of a certain plant species as a data set, and separating the high-spectrum data set G through a high-pass filter template to obtain noise GNSum signal GSThen, there are:
G＝GN+GS
obtaining a transformation matrix V by maximizing the signal-to-noise ratio of the transformed data, namely maximizing the ratio of the signal covariance to the noise covariance;
Figure BDA0002785258110000051
in the formula: y isNAnd YSNoise and signal respectively of the transformed data, Cov (-) representing covariance, CNAnd CSCovariance of the noise and signal, C, of the data, respectivelyN＝Cov(GN)，CS＝Cov(GS) The above optimization problem can be equivalently:
Figure BDA0002785258110000052
wherein C represents the overall covariance of the data, and C ═ CN+CSAccording to the lagrange multiplier method, the optimal solution of the above equation is:
CV＝λCNV
according to the above formula, the eigenvalues are arranged from large to small, and the eigenvectors corresponding to the first d eigenvalues are taken, so that a conversion matrix can be obtained:
V＝[V1,V2,…,Vd]
extracting pure pixels in the hyperspectral image by using a minimum noise separation result excluding a noise band to generate a pixel purity image, extracting the pure pixels by setting a threshold value, and performing the next calculation;
and extracting end-member spectrum curves of different plant species to generate a complete end-member spectrum curve and recording the complete end-member spectrum curve into an end-member spectrum database.
The data classification in the fifth step adopts a spectrum angle classification method:
identifying the vegetation types of the whole green land area by using a spectral angle mapping method, wherein the specific identification method comprises the steps of extracting an average end member spectral curve in each contour identification result, and comparing data in a wave spectrum database by using the spectral angle mapping method for identification, and the specific method comprises the following steps:
Figure BDA0002785258110000053
in the formula: alpha is the comparison result of the included angle between the space vector t of the unknown wave spectrum and the space vector r of the sample wave spectrum, tiAnd riRespectively representing the values of an unknown spectrum t and a sample spectrum r on the ith wave band, referring to end member spectra in a spectrum library, and n is the number of wave bands. The value of the calculation result is 0-pi/2, and the closer the value is to 0, the closer the test pixel is to the reference spectrum. And calculating a spectrum angle between one pixel spectrum and all reference spectra, wherein the belonging ground object type is the ground object type represented by the alpha minimum reference spectrum in all calculation results.
The method for calculating the crown defect rate of the urban green land vegetation further comprises the following steps:
and seventhly, evaluating the fitting effect and precision of the model by adopting an average error ME, an average relative error MPE and a root mean square error RMSE, wherein the calculation formula is as follows:
Figure BDA0002785258110000061
Figure BDA0002785258110000062
Figure BDA0002785258110000063
in the formula, xiThe method is a measured value of crown defect rate of urban green land plants;
Figure BDA0002785258110000064
estimating a value for a crown defect rate model of urban green land plants; n is the number of the same plots; i is a certain pattern.
The method for calculating the crown defect rate of the urban green land trees, provided by the invention, is based on the fact that high-spectrum data, high-resolution panchromatic remote sensing images and measured data are used as data sources, and a set of method for calculating the crown defect rate of the urban green land trees on the premise of identifying the contour of a vegetation crown, extracting the spectral end member spectral curve characteristics and identifying species based on spectral angle comparison is constructed.
Drawings
The invention is further illustrated by the following examples in conjunction with the accompanying drawings:
FIG. 1 is a full color image and a hyperspectral image after preprocessing in accordance with the present invention;
FIG. 2 is a NDVI calculation and an Otsu threshold segmentation image;
FIG. 3 is a graph of end-member spectral data for different vegetation types after MNF and PPI;
FIG. 4 is a schematic diagram of a crown defect rate model;
FIG. 5 is a schematic diagram of a calculation result of the crown defect rate of urban green land vegetation.
Detailed Description
As shown in fig. 1-5, the method for calculating the crown defect rate of urban green land vegetation comprises the following steps:
acquiring panchromatic images and hyperspectral image data of a target area, and preprocessing the data;
secondly, segmenting the hyperspectral image, and distinguishing vegetation and non-vegetation areas in the image data to form a vegetation image data set of the urban green land, wherein the vegetation image data set comprises all trees, shrubs and grasslands in the urban green land;
thirdly, carrying out edge detection on the image of the plant area by utilizing the full-color image, and identifying the canopy contour of the plant species;
combining visible light, red-edge wave bands and near-infrared wave bands of the hyperspectral images, positioning, identifying and gathering the purest pixels in data, acquiring pure end-member wave spectrum curves in the images of the urban green land plant areas, and extracting and storing end-member wave spectrum characteristic curves of different arbor and shrub plant species;
surveying and recording typical plant types and geographical positions of urban greenhouses, wherein the plant types comprise different kinds of trees, shrubs and grasslands;
classifying all data in the vegetation image according to end-member spectrum characteristic curves of different plant species under the hyperspectral image constructed in the step three, then superposing the data with canopy contour monitoring results of the plant species with the panchromatic image, and identifying and defining the species category of each contour;
step six, establishing a crown damage rate model according to the contour and species identification result, and calculating the crown damage rate of arbors and shrubs;
and seventhly, evaluating the fitting effect and precision of the model by adopting an average error ME, an average relative error MPE and a root mean square error RMSE, wherein the calculation formula is as follows:
Figure BDA0002785258110000071
Figure BDA0002785258110000072
Figure BDA0002785258110000073
in the formula, xiThe method is a measured value of crown defect rate of urban green land plants;
Figure BDA0002785258110000074
estimating a value for a crown defect rate model of urban green land plants; n is the number of the same plots; i is a certain pattern.
As shown in fig. 4, the specific process of the step six is as follows:
defining a polygon m at the outermost edge of the projected outline of the crown as
Figure BDA0002785258110000077
The polygon has PnA vertex Pn(n-1, 2, …, m) is arranged along the positive direction of the boundary, and the coordinates are (x) in sequence1，y1)，(x2，y2)，…，(xm，ym) Establishing
Figure BDA0002785258110000076
The centroid of the polygon m and any two adjacent vertexes form a triangle, the area of the triangle is obtained by the outer product of two plane vectors formed by the three vertexes, and the area of the polygon m at the outermost edge of the crown projection
Figure BDA0002785258110000075
Comprises the following steps:
Figure BDA0002785258110000081
in the polygon at the outermost edge of the crown projection, the exposed part without the blade shading is formed by k (k is 1,2, …, q) polygons ZkComposition of each polygon having OiOne vertex, vertex Oi(i-1, 2, …, j) is arranged along the positive direction of the boundary, and the coordinates are (a) in sequence1，b1)，(a2，b2)，…，(aj，bj) Establishing ZkPolygon area vector diagram, polygon ZkThe centroid of the tree and any two adjacent vertexes form a triangle, the area of the triangle is obtained by the outer product of two plane vectors formed by the three vertexes, and the area S of the exposed part without blade shielding in the polygon of the outermost edge of the crown projectionZComprises the following steps:
Figure BDA0002785258110000082
Figure BDA0002785258110000083
ideal crown projection area S of plantQComprises the following steps:
Figure BDA0002785258110000084
the crown defect rate P of the plants is as follows:
Figure BDA0002785258110000085
as shown in fig. 1, the preprocessing of the data in the first step includes: respectively carrying out radiometric calibration atmospheric correction, geometric correction, mosaic and splicing on the panchromatic image and the hyperspectral image.
As shown in fig. 2, a specific process of image segmentation in the second step is to segment the image by using the normalized vegetation index NDVI and the tsu algorithm.
First, NDVI is calculated, with the formula:
NDVI＝(ρNIR-ρRED)/(ρNIR+ρRED)
in the formula, ρNIRReflectivity in the near red band, pREDReflectance for the red band;
secondly, image segmentation is carried out by utilizing an Otsu algorithm, and the steps are as follows:
let the gray level of the NDVI image be L (G ═ 1,2,3, …, L), and the number of pixels at gray level i be niThe total pixel number of NDVI is shown as follows:
Figure BDA0002785258110000086
the histogram is normalized to:
Figure BDA0002785258110000091
Figure BDA0002785258110000092
wherein p (i) represents the probability of the appearance of the pixel with the gray level i in the image;
and (3) expressing a threshold value by T, and dividing pixels in the normalized image into two types according to the gray level T: c1T and C are not more than2T, the probability of occurrence of two classes is:
Figure BDA0002785258110000093
Figure BDA0002785258110000094
suppose that the mean values of the two types of pixels are m1And m2The average value of global pixels of the NDVI image is mgThen:
w1*m1+w2*m2＝mg
w1+w2＝1
the variance is now:
δ2＝p1(m1-mg)2+p2(m2-mg)2
namely:
δ2＝p1p2(m1-m2)2
in the formula, delta2The inter-class variance is the maximum value, and the corresponding threshold T is the optimal segmentation threshold of the NDVI image.
The specific process in the third step is as follows:
the method is characterized in that edge detection is carried out based on a Sobel operator to achieve the effect of identifying the canopy contour, and the method comprises the following steps:
firstly, setting a Sobel operator, wherein the formula is as follows:
Figure BDA0002785258110000095
Figure BDA0002785258110000096
secondly, performing plane convolution on a pixel set Z of the panchromatic high-resolution image and A and B to respectively calculate gray values Gx and Gy of the image detected by the transverse edge and the longitudinal edge, wherein the formula is as follows:
Figure BDA0002785258110000101
Figure BDA0002785258110000102
finally, the horizontal and vertical gray scale of each pixel of the panchromatic high-resolution image can be calculated:
Figure BDA0002785258110000103
using an approximation:
|f(Gx，Gy)|＝|Gx|+|Gy|
setting the threshold value as h and the pixel value f (G)x，Gy) If the sum of the edge points is larger than h, the edge points are marked as edge points, and the set of the edge points is the result of contour recognition;
then, on the basis of full-color image contour identification, referring to the type and the geographic position of a sample plot plant species, taking a crown high-spectrum image of a certain plant species as a data set, and separating the high-spectrum data set G through a high-pass filter template to obtain noise GNSum signal GSThen, there are:
G＝GN+GS
obtaining a transformation matrix V by maximizing the signal-to-noise ratio of the transformed data, namely maximizing the ratio of the signal covariance to the noise covariance;
Figure BDA0002785258110000104
in the formula: y isNAnd YSNoise and signal respectively of the transformed data, Cov (-) representing covariance, CNAnd CSCovariance of the noise and signal, C, of the data, respectivelyN＝Cov(GN)，CS＝Cov(GS) The above optimization problem can be equivalently:
Figure BDA0002785258110000105
wherein C represents the overall covariance of the data, and C ═ CN+CSAccording to the lagrange multiplier method, the optimal solution of the above equation is:
CV＝λCNV
according to the above formula, the eigenvalues are arranged from large to small, and the eigenvectors corresponding to the first d eigenvalues are taken, so that a conversion matrix can be obtained:
V＝[V1,V2,…,Vd]
extracting pure pixels in the hyperspectral image by using a minimum noise separation result excluding a noise band to generate a pixel purity image, extracting the pure pixels by setting a threshold value, and performing the next calculation;
and extracting end-member spectrum curves of different plant species to generate a complete end-member spectrum curve and recording the complete end-member spectrum curve into an end-member spectrum database.
The data classification in the fifth step adopts a spectrum angle classification method:
identifying the vegetation types of the whole green land area by using a spectral angle mapping method, wherein the specific identification method comprises the steps of extracting an average end member spectral curve in each contour identification result, and comparing data in a wave spectrum database by using the spectral angle mapping method for identification, and the specific method comprises the following steps:
Figure BDA0002785258110000111
in the formula: alpha is the comparison result of the included angle between the space vector t of the unknown wave spectrum and the space vector r of the sample wave spectrum, tiAnd riRespectively representing the values of an unknown spectrum t and a sample spectrum r on the ith wave band, referring to end member spectra in a spectrum library, and n is the number of wave bands. The value of the calculation result is 0-pi/2, and the closer the value is to 0, the closer the test pixel is to the reference spectrum. And calculating a spectrum angle between one pixel spectrum and all reference spectra, wherein the belonging ground object type is the ground object type represented by the alpha minimum reference spectrum in all calculation results.
Example 1:
1. the experimental zone of this example is located in Shanghai Shangxian zone Wu bridge nursery base (121.41 DEG E, 30.96 DEG N) and has an area of about 3.83hm 2. The vegetation types are rich. Wherein the plants comprise: the plant comprises: malus halliana Koehne (Malus halliana Koehne), Loropetalum chinense (R.Br.), Oliver var. rubrum Yueh, Elaeagnus pungens (Thunb.), Acer palmatum (Acer palmatum Thunb.), Viburnum Diatum Thunb (Viburnum Diatum Thunb.), Euonymus japonicus (Euonymus japonica Thunb.), Eurya japonica (Eurya-marrhena Hort.), Liriodendron tulipifera (Liriodendron tulipifera), Gleditsia japonica (Gleditsia trifolia Thailanthus chinensis 'Sunburt'), Hibiscus syrus Linatus (Hisculus), Hibiscus syriacus Linensis (bark), Sabourne mountain ash (Sabourne), Magnolia officinalis (Magnolia officinalis), Magnolia officinalis (Rosa) and so on.
2. Collecting high-resolution images, hyperspectral data and laser radar data by means of an unmanned aerial vehicle;
3. 28 square patterns (10X 10 m) were set in the area of the study2) The coordinate of the center point of the sample is measured by using GPS (Trimble GeoXH6000), the GPS is positioned by receiving wide-area differential signals, the precision is better than 0.5m, and the sample is processed by the GPSRecording vegetation types and conditions thereof in the ground, wherein the vegetation types comprise grassland types and tree and shrub species, and simultaneously measuring the crown width breast height and the tree height of each tree;
4. as shown in fig. 1, during data preprocessing, the high-resolution panchromatic images are spliced to generate a complete high-resolution image of a research area, meanwhile, in combination with ground actual measurement control point data, a quadratic polynomial model is adopted to carry out geometric fine correction on the image, radiation calibration is carried out on original hyperspectral data by means of sensor radiation calibration data, an empirical linear model is utilized to carry out atmospheric correction in combination with ground actual measurement target spectral data, meanwhile, in combination, the ground actual measurement control point is utilized to carry out geometric fine correction on the hyperspectral images, and then embedding and splicing are carried out;
5. as shown in fig. 2, the preprocessed hyperspectral image is used to calculate the NDVI to obtain an NDVI distribution map, and the Otsu method is used to obtain a threshold, which is 0.3 in this example, and the NDVI is segmented to obtain a vegetation area in the distribution map;
6. cutting the preprocessed high-resolution panchromatic image according to the vegetation area in the step 5 to obtain a vegetation distribution area of the panchromatic image;
7. sobel edge detection is carried out on the vegetation distribution area of the panchromatic image, and crown outlines and grassland outlines of different plant species are identified;
8. as shown in fig. 3, on the basis of identifying the outline of the plant species, the marked plant species in the sample plot of the field investigation is calibrated correspondingly, the hyperspectral image in the calibrated outline is processed, including MNF denoising and PPI pure pixel extraction, the end-member spectrum curves of different typical species are extracted and recorded into a warehouse;
9. as shown in fig. 3, other hyperspectral images in the calibration contour are extracted, the spectrum curve of the unknown plant is compared with the end member spectrum curve recorded in the database, and the unknown species are identified by using a method of spectrum angle comparison;
10. as shown in fig. 4 and 5, the contour recognition result is combined with the spectral curve classification result to obtain the contours and species types of grassland, trees and shrubs, the trees and shrubs are screened out, the crown projection boundary is calculated, the polygonal boundary of the exposed part inside the crown projection boundary is obtained, and the crown defect rate is calculated in combination with the crown defect rate calculation model.
Taking tree E as an example, the projection outline of the crown obtained by outline recognition is automatically calculated in a geographic information system GIS to obtain the vertex coordinates of the outline, and the area of the outline obtained by formula calculation is 31.72m2，
In the polygon at the outermost edge of the projection of the crown, the exposed part without the shielding of the blade consists of 1 polygon, the area of the polygon is automatically calculated in a Geographic Information System (GIS) according to the vector principle, and the area of the exposed part is 10.52m2，
The projected area of the ideal crown of the tree E is a circle with the center of the outline as a circular point and the average radius of the outline of 4.3m, and the area is 58.06m2；
The crown defect rate of tree E is:
Figure BDA0002785258110000131
11. carrying out precision evaluation on the calculation result by using the average error ME, the average relative error MPE and the root mean square error RMSE;
evaluation index Accuracy of measurement
ME(％) 76.21％
MPE 0.89
RMSE 0.92
12. In the embodiment, the crown defect rate model is established based on the extracted characteristic variables by combining the high-resolution panchromatic image and the hyperspectral image data acquired at the same time. The data source records vegetation structure information from the angles of texture, spectrum and space, the data of the vegetation structure information are complementary, and the crown defect rate can be rapidly calculated. The method enhances the capacity and the precision of measuring and calculating the crown defect rate of the urban green land vegetation; the data verification results in tables 1 and 2 show that the urban green land vegetation crown defect rate is calculated through the method, and compared with the result of manual investigation, the method is high in precision, convenient and fast.

Claims (7)

1. The urban green land tree crown defect rate calculation method is characterized by comprising the following steps:
acquiring panchromatic images and hyperspectral image data of a target area, and preprocessing the data;
secondly, segmenting the hyperspectral image, and distinguishing vegetation and non-vegetation areas in the image data to form a vegetation image data set of the urban green land, wherein the vegetation image data set comprises all trees, shrubs and grasslands in the urban green land;
thirdly, carrying out edge detection on the image of the plant area by utilizing the full-color image, and identifying the outline of the plant species;
combining visible light, red-edge wave bands and near-infrared wave bands of the hyperspectral images, positioning, identifying and gathering the purest pixels in data, acquiring pure end members in the images of the urban green land plant areas, and extracting end member spectrum characteristic curves of different plant species and storing the end members in a warehouse;
surveying and recording typical plant types and geographical positions of urban greenhouses, wherein the plant types comprise different kinds of trees, shrubs and grasslands;
classifying all data in the hyperspectral vegetation image according to end member spectrum characteristic curves of different plant species in the hyperspectral image obtained in the step three, and then overlapping the data with the outline identification result of the panchromatic image to define the species category of each outline;
and step six, establishing a crown damage rate model according to the contour and species identification result, and calculating the crown damage rate of the trees and shrubs.
2. The method for calculating the crown defect rate of the urban green land trees according to claim 1, wherein the concrete process of the sixth step is as follows:
defining a polygon m at the outermost edge of the projected outline of the crown as
Figure FDA0002785258100000011
The polygon has PnA vertex Pn(n-1, 2, …, m) is arranged along the positive direction of the boundary, and the coordinates are (x) in sequence1，y1)，(x2，y2)，…，(xn，yn) Establishing
Figure FDA0002785258100000012
The centroid of the polygon m and any two adjacent vertexes form a triangle, the area of the triangle is obtained by the outer product of two plane vectors formed by the three vertexes, and the area of the polygon m at the outermost edge of the crown projection
Figure FDA0002785258100000013
Comprises the following steps:
Figure FDA0002785258100000014
in the polygon at the outermost edge of the crown projection, the exposed part without the blade shading is formed by k (k is 1,2, …, q) polygons ZkComposition of each polygon having OiOne vertex, vertex Oi(i-1, 2, …, j) is arranged along the positive direction of the boundary, and the coordinates are (a) in sequence1，b1)，(a2，b2)，…，(ai，bi) Establishing ZkPolygon area vector diagram, polygon ZkThe centroid of the tree and any two adjacent vertexes form a triangle, the area of the triangle is obtained by the outer product of two plane vectors formed by the three vertexes, and the area S of the exposed part without blade shielding in the polygon of the outermost edge of the crown projectionZComprises the following steps:
Figure FDA0002785258100000021
Figure FDA0002785258100000022
ideal crown projection area S of plantQComprises the following steps:
Figure FDA0002785258100000023
the crown defect rate P of the plants is as follows:
Figure FDA0002785258100000024
3. the method for calculating the crown defect rate of urban green land trees according to claim 1, wherein the preprocessing of the data in the first step comprises: and respectively carrying out radiometric calibration, atmospheric correction, geometric correction, mosaic and splicing on the panchromatic image and the hyperspectral image.
4. The method for calculating the crown defect rate of urban green land vegetation according to claim 1, wherein the specific process of image segmentation in the second step is to segment the image by using a normalized vegetation index NDVI and an Otsu algorithm.
First, NDVI is calculated, with the formula:
NDVI＝(ρNIR-ρRED)/(ρNIR+ρRED)
in the formula, ρNIRReflectivity in the near red band, pREDReflectance for the red band;
secondly, image segmentation is carried out by utilizing an Otsu algorithm, and the steps are as follows:
let the gray level of the NDVI image be L (G ═ 1,2,3, …, L), and the number of pixels at gray level i be niThe total pixel number of NDVI is shown as follows:
Figure FDA0002785258100000025
the histogram is normalized to:
Figure FDA0002785258100000026
Figure FDA0002785258100000031
wherein p (i) represents the probability of the appearance of the pixel with the gray level i in the image;
and (3) expressing a threshold value by T, and dividing pixels in the normalized image into two types according to the gray level T: c1T and C are not more than2T, the probability of occurrence of two classes is:
Figure FDA0002785258100000032
Figure FDA0002785258100000033
suppose that the mean values of the two types of pixels are m1And m2The average value of global pixels of the NDVI image is mgThen:
w1*m1+w2*m2＝mg
w1+w2＝1
the variance is now:
δ2＝p1(m1-mg)2+p2(m2-mg)2
namely:
δ2＝p1p2(m1-m2)2
in the formula, delta2The inter-class variance is the maximum value, and the corresponding threshold T is the optimal segmentation threshold of the NDVI image.
5. The method for calculating the crown defect rate of the urban green land trees according to claim 1, wherein the specific process in the third step is as follows:
the method is characterized in that edge detection is carried out based on a Sobel operator to achieve the effect of identifying the canopy contour, and the method comprises the following steps:
firstly, setting a Sobel operator, wherein the formula is as follows:
Figure FDA0002785258100000034
Figure FDA0002785258100000035
secondly, performing plane convolution on a pixel set Z of the panchromatic high-resolution image and A and B to respectively calculate gray values Gx and Gy of the image detected by the transverse edge and the longitudinal edge, wherein the formula is as follows:
Figure FDA0002785258100000041
Figure FDA0002785258100000042
finally, the horizontal and vertical gray scale of each pixel of the panchromatic high-resolution image can be calculated, and the formula is as follows:
Figure FDA0002785258100000043
using an approximation:
|f(Gx，Gy)|＝|Gx|+|Gy|
setting the threshold value as h and the pixel value f (G)x，Gy) If the sum of the edge points is larger than h, the edge points are marked as edge points, and the set of the edge points is the result of contour recognition;
then, on the basis of full-color image contour identification, referring to the type and the geographic position of a sample plot plant species, taking a crown high-spectrum image of a certain plant species as a data set, and separating the high-spectrum data set G through a high-pass filter template to obtain noise GNSum signal GSThen, there are:
G＝GN+GS
obtaining a transformation matrix V by maximizing the signal-to-noise ratio of the transformed data, namely maximizing the ratio of the signal covariance to the noise covariance;
Figure FDA0002785258100000044
in the formula: y isNAnd YSNoise and signal respectively of the transformed data, Cov (-) representing covariance, CNAnd CSCovariance of the noise and signal, C, of the data, respectivelyN＝Cov(GN)，CS＝Cov(GS) To aboveThe formula optimization problem can be equivalently:
Figure FDA0002785258100000045
wherein C represents the overall covariance of the data, and C ═ CN+CSAccording to the lagrange multiplier method, the optimal solution of the above equation is:
CV＝λCNV
according to the above formula, the eigenvalues are arranged from large to small, and the eigenvectors corresponding to the first d eigenvalues are taken, so that a conversion matrix can be obtained:
V＝[V1,V2,…,Vd]
extracting pure pixels in the hyperspectral image by using a minimum noise separation result excluding a noise band to generate a pixel purity image, extracting the pure pixels by setting a threshold value, and performing the next calculation;
and extracting end-member spectrum curves of different plant species to generate a complete end-member spectrum curve and recording the complete end-member spectrum curve into an end-member spectrum database.
6. The method for calculating the crown defect rate of urban green land trees according to claim 1, wherein the data classification in the fifth step adopts a spectral angle classification method:
identifying the vegetation types of the whole green land area by using a spectral angle mapping method, wherein the specific identification method comprises the steps of extracting an average end member spectral curve in each contour identification result, and comparing data in a wave spectrum database by using the spectral angle mapping method for identification, and the specific method comprises the following steps:
Figure FDA0002785258100000051
in the formula: alpha is the comparison result of the included angle between the space vector t of the unknown wave spectrum and the space vector r of the sample wave spectrum, tiAnd riRespectively represent unknown spectra t and samples in the ith wave bandThe value of the spectrum r refers to the end member spectrum in the spectrum library, and n is the number of wave bands. The value of the calculation result is 0-pi/2, and the closer the value is to 0, the closer the test pixel is to the reference spectrum. And calculating a spectrum angle between one pixel spectrum and all reference spectra, wherein the belonging ground object type is the ground object type represented by the alpha minimum reference spectrum in all calculation results.
7. The method for calculating the crown defect rate of the urban green land trees according to claim 1, which is characterized by comprising the following steps: and seventhly, evaluating the fitting effect and precision of the model by adopting an average error ME, an average relative error MPE and a root mean square error RMSE, wherein the calculation formula is as follows:
Figure FDA0002785258100000052
Figure FDA0002785258100000053
Figure FDA0002785258100000054
in the formula, xiThe method is a measured value of crown defect rate of urban green land plants;
Figure FDA0002785258100000055
estimating a value for a crown defect rate model of urban green land plants; n is the number of the same plots; i is a certain pattern.
CN202011295757.XA 2020-11-18 2020-11-18 Urban green land tree crown loss rate calculation method Active CN112577907B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011295757.XA CN112577907B (en) 2020-11-18 2020-11-18 Urban green land tree crown loss rate calculation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011295757.XA CN112577907B (en) 2020-11-18 2020-11-18 Urban green land tree crown loss rate calculation method

Publications (2)

Publication Number Publication Date
CN112577907A true CN112577907A (en) 2021-03-30
CN112577907B CN112577907B (en) 2023-03-31

Family

ID=75122948

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011295757.XA Active CN112577907B (en) 2020-11-18 2020-11-18 Urban green land tree crown loss rate calculation method

Country Status (1)

Country Link
CN (1) CN112577907B (en)

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2009115787A (en) * 2009-04-27 2010-11-10 Государственное учреждение "Научный центр Проблем аэрокосмического мониторинга"-ЦПАМ "АЭРОКОСМОС" (RU) METHOD FOR DETERMINING ECOLOGICAL STATE OF FORESTS
US20140093138A1 (en) * 2011-06-29 2014-04-03 Fujitsu Limited Plant species identification apparatus and method
CN104463164A (en) * 2014-09-03 2015-03-25 中国科学院遥感与数字地球研究所 Tree canopy structure information extraction method based on rib method and crown height ratio
CN105354534A (en) * 2015-09-29 2016-02-24 南京林业大学 Tree species classification method based on multi-source simultaneous high-resolution remote sensing data
CN105825217A (en) * 2016-03-22 2016-08-03 辽宁师范大学 Hyperspectral image interested area automatic extraction method based on active contour model
JP2017163934A (en) * 2016-03-17 2017-09-21 国立大学法人信州大学 Methods of calculating damage classification of pine wilt disease and apparatus for calculating damage classification of pine wilt disease
CN109115951A (en) * 2018-07-31 2019-01-01 东北农业大学 The full nitrogen estimating and measuring method of rice plant based on canopy structure and canopy spectra
CN111091030A (en) * 2018-10-24 2020-05-01 中国测绘科学研究院 Tree species identification method and device, computer equipment and readable storage medium
CN111539580A (en) * 2020-04-30 2020-08-14 上海市园林科学规划研究院 Multi-scheme optimization method for urban greening ecological technology integration application

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2009115787A (en) * 2009-04-27 2010-11-10 Государственное учреждение "Научный центр Проблем аэрокосмического мониторинга"-ЦПАМ "АЭРОКОСМОС" (RU) METHOD FOR DETERMINING ECOLOGICAL STATE OF FORESTS
US20140093138A1 (en) * 2011-06-29 2014-04-03 Fujitsu Limited Plant species identification apparatus and method
CN104463164A (en) * 2014-09-03 2015-03-25 中国科学院遥感与数字地球研究所 Tree canopy structure information extraction method based on rib method and crown height ratio
CN105354534A (en) * 2015-09-29 2016-02-24 南京林业大学 Tree species classification method based on multi-source simultaneous high-resolution remote sensing data
JP2017163934A (en) * 2016-03-17 2017-09-21 国立大学法人信州大学 Methods of calculating damage classification of pine wilt disease and apparatus for calculating damage classification of pine wilt disease
CN105825217A (en) * 2016-03-22 2016-08-03 辽宁师范大学 Hyperspectral image interested area automatic extraction method based on active contour model
CN109115951A (en) * 2018-07-31 2019-01-01 东北农业大学 The full nitrogen estimating and measuring method of rice plant based on canopy structure and canopy spectra
CN111091030A (en) * 2018-10-24 2020-05-01 中国测绘科学研究院 Tree species identification method and device, computer equipment and readable storage medium
CN111539580A (en) * 2020-04-30 2020-08-14 上海市园林科学规划研究院 Multi-scheme optimization method for urban greening ecological technology integration application

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
DUAN DAN-DAN等: "Estimating total leaf nitrogen concentration in winter wheat by canopy hyperspectral data and nitrogen vertical distribution", 《JOURNAL OF INTEGRATIVE AGRICULTURE》 *
李姗姗等: "基于Geoeye-1影像光谱特性的红树林冠层种类识别", 《光谱学与光谱分析》 *
梁志林等: "高光谱遥感城市植被识别方法研究", 《地理空间信息》 *

Also Published As

Publication number Publication date
CN112577907B (en) 2023-03-31

Similar Documents

Publication Publication Date Title
CN112541921B (en) Urban green land vegetation information data accurate determination method
Simonse et al. Automatic determination of forest inventory parameters using terrestrial laser scanning
CN104820830B (en) A kind of wood recognition method based on Full wave shape LiDAR canopy section models
CN112669363B (en) Method for measuring three-dimensional green space of urban green space
CN112577906B (en) Urban green land soil moisture content detection method
CN113033670A (en) Method for extracting rice planting area based on Sentinel-2A/B data
CN114821362B (en) Multi-source data-based rice planting area extraction method
CN113935366B (en) Automatic classification method for single-tree segmentation of point cloud
CN113505635A (en) Method and device for identifying winter wheat and garlic mixed planting area based on optics and radar
Zhang et al. Opportunities of UAVs in orchard management
CN111487643A (en) Building detection method based on laser radar point cloud and near-infrared image
Islam et al. A ground-based platform for reliable estimates of fruit number, size, and color in stone fruit orchards
CN112577954B (en) Urban green land biomass estimation method
CN112577907B (en) Urban green land tree crown loss rate calculation method
CN115223062B (en) Eucalyptus artificial forest area stand accumulation amount time difference correction method based on UAV data
Rinnamang et al. Estimation of aboveground biomass using aerial photogrammetry from unmanned aerial vehicle in teak (Tectona grandis) plantation in Thailand
CN116205678A (en) Agricultural land physical quantity investigation and value quantity estimation method based on remote sensing automatic interpretation
Danoedoro et al. Combining pan-sharpening and forest cover density transformation methods for vegetation mapping using Landsat-8 Satellite Imagery
Hall et al. A method for extracting detailed information from high resolution multispectral images of vineyards
Rudjord et al. Tree species classification with hyperspectral imaging and lidar
Jiménez-Bello et al. Use of remote sensing and geographic information tools for irrigation management of citrus trees
Feng et al. Dection and Health Analysis of Individual Tree in Urban Environment with Multi-Sensor Platform
Iskandar et al. Crown closure segmentation on wetland lowland forest using the mean shift algorithm
Wang et al. Arial images and LiDAR fusion applied in forest boundary detection
Li et al. High-throughput plant height estimation from RGB Images acquired with Aerial platforms: a 3D point cloud based approach

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant