CN108038433B - Urban tree carbon content estimation method based on multi-echo airborne laser scanning data - Google Patents

Urban tree carbon content estimation method based on multi-echo airborne laser scanning data Download PDF

Info

Publication number
CN108038433B
CN108038433B CN201711248101.0A CN201711248101A CN108038433B CN 108038433 B CN108038433 B CN 108038433B CN 201711248101 A CN201711248101 A CN 201711248101A CN 108038433 B CN108038433 B CN 108038433B
Authority
CN
China
Prior art keywords
data
tree
crown
als
echo
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.)
Active
Application number
CN201711248101.0A
Other languages
Chinese (zh)
Other versions
CN108038433A (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.)
Xiamen University
Original Assignee
Xiamen University
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 Xiamen University filed Critical Xiamen University
Priority to CN201711248101.0A priority Critical patent/CN108038433B/en
Publication of CN108038433A publication Critical patent/CN108038433A/en
Application granted granted Critical
Publication of CN108038433B publication Critical patent/CN108038433B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/188—Vegetation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/88—Lidar systems specially adapted for specific applications
    • G01S17/89—Lidar systems specially adapted for specific applications for mapping or imaging
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/48—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00
    • G01S7/4802—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2411—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on the proximity to a decision surface, e.g. support vector machines
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
    • 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
    • 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

The invention discloses an urban tree carbon content estimation method based on multi-echo airborne laser scanning data, which comprises the following steps: s1, scanning urban trees by using a multispectral airborne laser scanner to obtain multi-echo ALS point cloud data of a plurality of navigation bands; s2, preprocessing multi-echo ALS point cloud data of a plurality of flight strips to obtain nDSM data and intensity images of all laser channels; s3, classifying the nDSM data, extracting tree type data in the nDSM data, and obtaining a crown height model CHM; s4, estimating shape parameters of the tree by the crown height model CHM, wherein the shape parameters of the tree comprise the crown height and the crown diameter; s5, firstly, constructing an ALS-DBH regression model by using the empirical model, and then selecting a plurality of tree samples on site for verification; s6, estimating carbon content from DBH and crown height of all species. The method can be used for quickly and nondestructively estimating the carbon content of the urban trees based on the ALS tree shape parameters and the abnormal-speed growth model.

Description

Urban tree carbon content estimation method based on multi-echo airborne laser scanning data
Technical Field
The invention relates to an urban tree carbon content estimation method based on multi-echo airborne laser scanning data.
Background
As the population grows, townships become an irresistible trend. As cities expand, ecological and ground attachments change. Ground attachments have changed from the natural environment to various urban structures such as: watertight highways, tall buildings, parks, etc. These transitions present some pressing problems, such as the urban heat island effect. Urban expansion and economic activities also increase energy consumption and lead to greenhouse gas emissions. The vegetation in the city plays an important role in improving the urban climate, for example, tree shading and wind shielding effects can reduce direct solar heat and air infiltration into the house. Meanwhile, the urban vegetation stores a large amount of carbon, so that the urban heat island effect can be reduced to a great extent. Therefore, it is important to estimate and monitor the carbon content of urban trees.
Estimating the carbon content of urban trees requires extensive on-site data, which, if measured manually, would require significant labor and time. In addition, estimating the dry weight of trees is often destructive, and therefore, it is necessary to provide a quick and efficient method for estimating the carbon content of urban trees.
Disclosure of Invention
The invention aims to overcome the defects of the prior art and provides an urban tree carbon content estimation method based on multi-echo airborne laser scanning data.
In order to achieve the purpose, the invention adopts the following technical scheme:
the urban tree carbon content estimation method based on multi-echo airborne laser scanning data comprises the following steps:
s1, scanning urban trees by using an airborne laser scanner comprising green laser, near-infrared laser and short-wave infrared laser to obtain multi-echo ALS point cloud data of a plurality of flight zones;
s2, preprocessing multi-echo ALS point cloud data of a plurality of flight strips to obtain nDSM data and intensity images of all laser channels;
s3, vegetation classification: classifying the nDSM data, and extracting tree category data in the nDSM data to obtain a crown height model CHM;
s4, estimating shape parameters of the tree by the crown height model CHM, wherein the shape parameters of the tree comprise the crown height and the crown diameter;
s5, selection of ALS-DBH regression model: firstly, establishing an ALS-DBH regression model by an empirical model, then obtaining an even number of sampling groups by selecting a plurality of tree samples on site, taking one half of the sampling groups as model groups for model development, taking the other half of the sampling groups as verification groups, carrying out model development on all combination conditions, fitting through a cross verification process, and collecting parameters generated by each combination; by determining the coefficient R2Checking the predictive power and performance of the regression model and the accuracy of the predictions by RMSE of the predicted parameters, comparing R in all combination cases2And RMSE, selected to have a high R in model fitting2And an ALS-DBH regression model with low RMSE in the verification to predict DBH;
s6, prediction of carbon content:
calculating total biomass from DBH and crown height of all species, estimating the amount of carbon storage stored in the tree to be half of total biomass, comparing the predicted carbon storage with the estimated carbon storage by DBH and height measured on site, and by RMSE and R2Evaluation was performed.
Preferably, the wavelengths of the green laser, the near-infrared laser and the short-wave infrared laser in step S1 are 532nm, 1064nm and 1550nm, respectively.
Preferably, step S2 specifically includes:
s21, denoising the multi-echo ALS point cloud data of each flight band respectively;
firstly, removing points with intensity outside 3 sigma by adopting an intensity filter;
secondly, counting the average distance between a certain point and a plurality of points adjacent to the certain point by adopting a statistical filter, and removing the point if the average distance exceeds the sum of the global average distance and a standard deviation;
s22, carrying out amplitude normalization on the intensity of the multi-echo ALS point cloud data of each denoised navigation band;
Figure BDA0001491103560000031
wherein I' is the amplitude normalized intensity, IrawIs the ALS point cloud data raw intensity, a is an index determined by the target geometry, R is the aircraft to target distancerefIs the average reference distance;
s23, generating an intensity image of each channel: fusing the normalized multi-echo ALS point cloud data of a plurality of flight strips together, wherein each flight strip has point cloud data with three wavelengths, and rasterizing the fused point cloud data, so that the pixel value in each grid is the intensity average value of all points in the grid; for the grids without points inside, carrying out linear interpolation on pixel values of adjacent grids to obtain the grids without points inside;
s24, classification of ground points and non-ground points;
fusing the normalized multi-echo ALS point cloud data of a plurality of flight strips together, and classifying ground points and non-ground points of the fused point cloud data by using a slope-based filtering algorithm;
s25, rasterizing the non-ground point data to obtain nDSM data;
dividing the non-ground point cloud data into grids with the side length of L, taking the maximum height value of all points in the grids as the pixel value of the grids, and performing linear interpolation on the pixel value of the grids without points inside from the pixel value of the adjacent grids.
Preferably, in the nsmd data obtained in step S24, if there is a height difference of 2m or more between a pixel value of a certain lattice and its eight neighbors, the height average of its eight neighbors is used as the pixel value of the lattice.
Preferably, step S3 is specifically:
s31, selecting a training sample from input data by taking green laser intensity, near-infrared laser intensity, short-wave infrared laser intensity, nDSM, pNDWI and pNDVI as input, acquiring the class of the training sample through aerial photos of Google Earth, using the training sample for training an SVM classifier, classifying nDSM data by using the trained SVM classifier, and extracting the class of trees in the nDSM data;
wherein the content of the first and second substances,
Figure BDA0001491103560000032
Figure BDA0001491103560000033
CGreen，CNIRand CSWIRGreen laser intensity, near infrared laser intensity, and short wave infrared laser intensity, respectively;
s32, filtering the extracted tree types for a plurality of times by adopting a 3-by-3 mode filter;
s33, empirically testing the filtering times, generating an optimal crown, which refers to a crown shape that can remove isolated pixels and preserve the dominant tree, and generating the crown height model CHM.
Preferably, step S4 is specifically:
s41, obtaining a local maximum pixel set: filtering the CHM by adopting a 3-by-3 local maximum filter, defining pixel points with pixel values higher than eight neighborhoods as local maximum values, and obtaining a local maximum pixel set;
s42, determining treetop: further extracting the local maximum pixel set by using ALS intensity data, and recording pixel points which are also local maximum in an ALS intensity image as treetop points;
s43, taking the detected treetop as a mark, adopting a region growing segmentation algorithm to segment the CHM into independent crowns, wherein each treetop point after segmentation has a closed cut block, and expressing the tree separation effect by absolute precision, the method comprises the following steps:
Figure BDA0001491103560000041
wherein n is1，1Is the detected number of crown blocks, ntotalIs the actual crown number; tree height is defined as the local maximum in each crown blockThe diameter of the crown is defined as the average value of the maximum crown diameter passing through the center of the highest point of the crown and the diameter in the vertical direction.
Preferably, step S5 is specifically:
s51, constructing an ALS-DBH regression model:
by taking the tree height and the crown diameter as independent prediction factors, and constructing an ALS-DBH regression model according to an empirical model to predict DBH, if,
DBHField＝a·CDALS+b·HALS+c；
wherein CD is the diameter of the crown in meters; h is the height of the crown, and the unit is meter; a, b, c are parameters derived from empirical models;
Figure BDA0001491103560000051
Figure BDA0001491103560000052
Figure BDA0001491103560000053
wherein CD is the crown diameter in meters; h is the height of the crown, and the unit is meter; a, b, c refer to parameters derived from empirical models,
Figure BDA0001491103560000054
and
Figure BDA0001491103560000055
the mean values of DBH, CD and H, respectively;
s52, ALS-DBH regression model selection:
selecting a plurality of tree samples on site, dividing the tree samples into a groups a, wherein a is an even number, selecting the a/2 groups for model development, and using the a/2 groups for verification, so that the tree samples exist
Figure BDA0001491103560000056
In a combined situationWill be
Figure BDA0001491103560000057
The combination conditions were fitted by a cross-validation process, and the parameters generated by each combination were collected by determining the coefficient R2Checking the predictive power and performance of the regression model and checking the accuracy of the prediction by RMSE of the prediction parameters;
Figure BDA0001491103560000058
wherein y is a tree shape parameter, yfittedIs a tree shape parameter derived from ALS data,
Figure BDA0001491103560000059
is the shape parameter of trees actually measured in the field and compared
Figure BDA00014911035600000510
R of group ALS-DBH regression model2And RMSE, selected to have a high R in model fitting2And a regression model with low RMSE in the verification to predict DBH.
Preferably, in step S6, the total biomass is the sum of the biomass in leaves, branches, wood and bark, and then,
Figure BDA00014911035600000512
Figure BDA00014911035600000513
Figure BDA00014911035600000514
Figure BDA00014911035600000515
ytotal＝yfoliage+ybranches+ywood+ybark；
then the process of the first step is carried out,
Figure BDA00014911035600000511
where y is the dry weight of the tree component in Kg, D is DBH in cm, β is the model parameter, each tree component has a corresponding β, and e is the error compensation.
Preferably, step S7 is also included;
s7, creating a city carbon storage map:
s71, converting the crown part derived from the ALS data and the carbon content stored in the crown part into vector data in the GIS;
s72, dividing urban land into a plurality of coverage types, and dividing the total carbon storage amount in each coverage type by the total area of the coverage type to obtain the carbon storage amount per unit area of each coverage type;
and S73, performing area division on urban land according to land coverage types, multiplying the total area of each area by the unit area carbon storage amount of the corresponding land coverage type, estimating the carbon content of each area, and finally creating an urban carbon storage map.
After adopting the technical scheme, compared with the background technology, the invention has the following advantages:
1. the data adopted by the method is obtained through multispectral ALS, and the vegetation can be effectively separated according to different attachments with different reflectivities to different wavelengths, wherein the vegetation classification accuracy reaches 90.23 percent, and is 11 percent higher than that of a single spectrum;
2. the CHM is subjected to local maximum filtering to obtain treetops, and the treetops are finally left through maximum intensity image filtering, so that the treetop detection accuracy can be effectively improved;
3. the carbon content model is established by obtaining tree shape parameters (tree height, crown diameter and trunk diameter) through ALS data, so that the carbon content of urban trees can be rapidly calculated without damage, and an urban carbon content map is drawn.
Drawings
FIG. 1 is a block diagram of the overall process of the present invention;
FIG. 2 is a detailed flow diagram of the present invention;
FIG. 3 is a schematic diagram of a floor attachment classification;
FIG. 4 is a diagram illustrating the results of steps S4 according to the present invention, wherein 4(a) is a diagram illustrating the result of filtering the local maximum of the crown height model; FIG. 4(b) is a diagram illustrating the result of filtering the local maximum of the maximum intensity image; FIG. 4(c) shows the result of region growing segmentation;
fig. 5 is a map of urban tree carbon content created by the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
Examples
As shown in fig. 1 and 2, the urban tree carbon content estimation method based on multi-echo airborne laser scanning data includes the following steps:
s1, acquiring multi-echo ALS point cloud data of a plurality of flight belts;
s2, acquiring nDSM data and intensity images of all laser channels;
s3, classifying vegetation to obtain a crown height model CHM;
s4, obtaining the crown height and crown diameter of the tree;
s5, selecting an ALS-DBH regression model;
s6, predicting the carbon content;
s7 creates a city carbon storage map.
The wavelengths of the green laser, the near-infrared laser and the short-wave infrared laser in the step S1 are 532nm, 1064nm and 1550nm, respectively.
Step S2 specifically includes:
s21, denoising multi-echo ALS point cloud data of each flight band respectively:
firstly, removing points with intensity outside 3 sigma by using an intensity filter to remove specular reflection points with high intensity;
secondly, a statistical filter is adopted to count the average distance between a certain point and a plurality of adjacent points, and if the average distance exceeds the sum of the global average distance and one standard deviation, the point is removed.
S22, carrying out amplitude normalization on the intensity of the multi-echo ALS point cloud data of each denoised navigation band, wherein the normalization equation is as follows:
Figure BDA0001491103560000081
wherein I' is the amplitude normalized intensity, IrawIs the ALS point cloud data raw intensity, a is an index determined by the target geometry, R is the aircraft to target distancerefIs the average reference distance. In general, for objects with large surfaces, such as broad-leaved trees, the index a is preferably 2, while for point-like objects, such as coniferous trees, the index a is preferably 3, in the present invention a is 2.
S23, generating an intensity image of each channel: the normalized multi-echo ALS point cloud data of multiple swaths, each with point cloud data of three wavelengths, are fused together, and then each point cloud is rasterized into an intensity image with a 1m ground resolution. Where the size of the pixels is selected based on the dot spacing of the data set, the size of the pixels is selected to be close to the dot distance such that a majority of the pixels have at least one dot and the vertical distribution of the dots is significantly more prominent. Thus, the points are divided into grids with the side length of 1m, and the pixel value in the grids is the average value of the intensities of all the points in the grids; for a lattice without a dot inside, pixel values of adjacent lattices are linearly interpolated. Respective intensity images were generated for channels with wavelengths of 532nm, 1064nm and 1550nm, respectively.
S24, classification of ground points and non-ground points;
fusing the normalized multi-echo ALS point cloud data of a plurality of navigation bands together. And (4) carrying out ground point and non-ground point classification on the fused point cloud data by using a slope-based filtering algorithm.
S25, rasterizing the non-ground point data to obtain nDSM data:
dividing all ground points into grids with the side length of L, wherein L is 1 in the embodiment, and linearly interpolating grids without points in the grids to obtain rasterized data;
non-ground points are also rasterized by a similar method to obtain nDSM, except that the pixel value on the grid is changed into the maximum height of all points in the grid, so that the tree height can be reflected to the maximum;
the above-obtained nmsd contains pits or points with negative pixel values, and these pits are often generated by operations such as rasterization, ground point filtering, and interpolation. If a certain cell has a pixel value with a height difference of more than 2m from its eight neighbors, the cell is filled with its eight neighbors height mean.
Step S3 specifically includes:
s31, taking green laser intensity, near-infrared laser intensity, short-wave infrared laser intensity, nDSM, pNDWI and pNDVI as input, intercepting a training sample from input data, acquiring the class of the sample through aerial photos of Google Earth, using the training sample for training an SVM classifier, classifying nDSM data by using the trained SVM classifier, and extracting the class of trees in the nDSM data;
wherein the content of the first and second substances,
Figure BDA0001491103560000091
Figure BDA0001491103560000092
CGreen，GNIRand CSWIRThe green laser intensity, the near infrared laser intensity and the short wave infrared laser intensity are respectively, and the pNDWI and pNDVI indexes can well distinguish the man-made objects.
Since the study area is a simple residential area, as shown in fig. 3, the target can be divided into six types of ground attachments of a water area, a house, a road, a lawn, trees, and an open area, but the reflection intensity of red light becomes irregular due to water and organic matter in water, so that the water area is removed from the point cloud, and only five types of ground attachments are classified. And intercepting training samples from all input data, and acquiring the category of the samples through aerial photos of the Google Earth.
At least 600 pixels are selected as training samples for each ground attachment, and 1000 pixels occupy the maximum proportion for vegetation. And using the pixel points for training the SVM classifier. The SVM classifier was selected for performing land cover classification with multi-spectral ALS derived data because: due to the popularity of SVMs in a typical ALS-related classification study, the classification results of this study can be compared to previous studies. The kernels often used in the classification of remote sensing data are:
Figure BDA0001491103560000093
and the value of the nuclear width y is 0.167, and the value of the regularization parameter C is 100.
To estimate the quality of the classification result, a total of 522 reference pixels were selected by hierarchical random sampling of the study area and manually labeled as land cover type in order to assess the quality of the classification result. For each category, at least 60 reference pixels are labeled, and for vegetation coverage occupying a large portion of the study area, including over 100 reference pixels, an error matrix is then generated from the labeling errors and the missing errors for estimating the classification result quality. As shown in fig. 2
S32, filtering the extracted tree types for a plurality of times by adopting a 3-by-3 mode filter, so that salt and pepper noise can be effectively reduced;
s33, empirically testing the filtering times, to generate an optimal crown, which refers to a crown shape that can remove isolated pixels and preserve the dominant tree, and generate the crown height model CHM.
In this way, pits in the crown can be filled; while the filtered tree class points are used to mask the nsms to generate nsms of isolated trees. The nsms of an isolated tree are also referred to as Crown Height Models (CHMs), i.e., models that display tree positions in a top-down view by the crown and store height values in pixels.
Step S4 specifically includes:
s41, obtaining a local maximum pixel set: filtering the CHM by adopting a 3-by-3 local maximum filter, defining pixel points with pixel values higher than eight neighborhoods as local maximum values, and obtaining a local maximum pixel set; fig. 4(a) shows the local maximum filtering result of the crown height model.
S42, determining treetop: in order to eliminate the marking error related to local maximum filtering and detect real treetop pixel points, ALS intensity data is used for further extracting a local maximum pixel set, and the pixel points which are also local maximum in an ALS intensity image are marked as treetop points.
ALS intensity values depend not only on the reflectivity of the subject, but also on the distance between the sensor and the subject. The treetop pixel point has a higher intensity value and a higher pixel value. Based on the point, a local maximum pixel set consisting of at least 15 pixel points is further extracted, and only those pixel points which are also local maximum in the maximum intensity image are finally reserved as treetop points. Fig. 4(b) shows the result of the local maximum filtering of the maximum intensity image.
S43, taking the detected treetop as a mark, adopting a region growing segmentation algorithm to segment the CHM into independent crowns, wherein each treetop point after segmentation has a closed cut block, and the absolute precision representation of the tree separation effect has the following steps:
Figure BDA0001491103560000101
wherein n is1，1Is the detected number of crown blocks, ntotalIs the actual crown number; the tree height is defined as the mean of the local maximum values in each crown block, and the crown diameter is defined as the mean of the maximum crown diameter passing through the center of the highest point of the crown and the diameter in the vertical direction. As shown in fig. 4(c), the region growing segmentation result is shown.
Step S5 specifically includes:
s51, constructing an ALS-DBH regression model:
by using ALS to derive tree height and crown diameter as independent prediction factors, and constructing ALS-DBH regression model according to empirical model to predict DBH, if,
DBHField＝a·CDALS+b·HALS+c；
wherein CD is the diameter of the crown in meters; h is the height of the crown, and the unit is meter; a, b, c are parameters derived from empirical models;
Figure BDA0001491103560000111
Figure BDA0001491103560000112
Figure BDA0001491103560000113
wherein CD is the crown diameter in meters; h is the height of the crown, and the unit is meter; a, b, c refer to parameters derived from empirical models,
Figure BDA0001491103560000114
and
Figure BDA0001491103560000115
the mean values of DBH, CD and H, respectively;
s52, ALS-DBH regression model selection:
selecting a plurality of tree samples (40 in the embodiment) on site, dividing the tree samples into a groups, wherein a is an even number, selecting a/2 groups for model development, and using the a/2 groups for verification, wherein the tree samples exist
Figure BDA0001491103560000116
A combination of situations. In this example, a-4, there are 6 combinations, each of the 6 cases was modeled at 0.05 level of significance, and fitted through a cross-validation process, and the parameters generated by each combination were collected byDetermination of the coefficient R2The regression model was checked for prediction ability and performance, and the accuracy of the prediction was checked by RMSE of the prediction parameters.
Wherein the content of the first and second substances,
Figure BDA0001491103560000117
wherein y is a tree shape parameter, yfittedIs a tree shape parameter derived from ALS data,
Figure BDA0001491103560000118
is the shape parameter of trees actually measured in the field and compared
Figure BDA0001491103560000119
R of group ALS-DBH regression model2And RMSE, selected to have a high R in model fitting2And a regression model with low RMSE in the verification to predict DBH.
Step S6 calculates total biomass from DBH and crown height of all species, estimates the amount of carbon storage stored in the tree to be half of total biomass, compares the predicted carbon storage with the estimated carbon storage from DBH and height measured in situ, and passes RMSE and R2Evaluation was performed.
The equation set based on DBH and height for all species was chosen to calculate biomass because there is no attribute or species information in this example. Thus, the total biomass is the sum of the biomass in leaves, branches, wood and bark, and there are,
Figure BDA0001491103560000121
Figure BDA0001491103560000122
Figure BDA0001491103560000123
Figure BDA0001491103560000124
ytotal＝yfoliage+ybranches+ywood+ybark；
then the process of the first step is carried out,
Figure BDA0001491103560000125
where y is the dry weight of the tree component in Kg, D is DBH in cm, β is the model parameter, each tree component has a corresponding β, and e is the error compensation.
S7, creating a city carbon storage map:
s71, converting the crown part derived from the ALS data and the carbon content stored in the crown part into vector data in the GIS;
s72, dividing urban land into a plurality of coverage types, and dividing the total carbon storage amount in each coverage type by the total area of the coverage type to obtain the carbon storage amount per unit area of each coverage type;
s73, performing area division on the urban land according to the land cover type, multiplying the total area of each area by the carbon storage amount per unit area of the corresponding land cover type, estimating the carbon content of each area, and finally creating an urban carbon storage map, such as the map shown in fig. 5, which is the urban tree carbon content map.
By adopting the method, the proportion of correct classification in the classification of the ground point and the non-ground point is 96.4 percent, the correct classification of the trees in the classification of the ground attachments reaches 90.23 percent, and the accuracy of the classical ALS data is only 79 percent. In the tree shape parameter measurement, the root mean square error RMSE of the tree height is 1.21m, and the root mean square error RMSE of the crown diameter is 1.47 m. Correlation coefficient R of DBH in ALS-DBH model prediction20.80, root mean square error RMSE of 4.6 cm. The carbon content RMSE was 142 Kg. This result indicates that ALS-based tree shape parameters and paradoxical growth models can yield consistent performance and accurate estimates, fully revealing that multi-spectral ALS data is carbon storage at the estimated individual tree level and mapping vegetation in urban environmentsThe potential of (2).
The above description is only for the preferred embodiment of the present invention, but the scope of the present invention is not limited thereto, and any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope of the present invention are included in the scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.

Claims (7)

1. The urban tree carbon content estimation method based on multi-echo airborne laser scanning data is characterized by comprising the following steps of:
s1, scanning urban trees by using an airborne laser scanner comprising green laser, near-infrared laser and short-wave infrared laser to obtain multi-echo ALS point cloud data of a plurality of flight zones;
s2, preprocessing multi-echo ALS point cloud data of a plurality of flight strips to obtain nDSM data and intensity images of all laser channels, wherein the step S2 specifically comprises the following steps:
s21, denoising the multi-echo ALS point cloud data of each flight band respectively;
firstly, removing points with intensity outside 3 sigma by adopting an intensity filter;
secondly, counting the average distance between a certain point and a plurality of points adjacent to the certain point by adopting a statistical filter, and removing the point if the average distance exceeds the sum of the global average distance and a standard deviation;
s22, carrying out amplitude normalization on the intensity of the multi-echo ALS point cloud data of each denoised navigation band;
Figure FDA0002483551040000011
wherein I' is the amplitude normalized intensity, IrawIs the ALS point cloud data raw intensity, a is an index determined by the target geometry, R is the aircraft to target distancerefIs the average reference distance;
s23, generating an intensity image of each channel: fusing the normalized multi-echo ALS point cloud data of a plurality of aerial strips, wherein each aerial strip has point cloud data with three wavelengths, and rasterizing the fused point cloud data, so that the pixel value in each grid is the intensity average value of all the points in the grid; for the grids without points inside, carrying out linear interpolation on pixel values of adjacent grids to obtain the grids without points inside;
s24, classification of ground points and non-ground points;
fusing the normalized multi-echo ALS point cloud data of a plurality of flight strips together, and classifying ground points and non-ground points of the fused point cloud data by utilizing a slope-based filtering algorithm;
s25, rasterizing the non-ground point data to obtain nDSM data;
dividing non-ground point cloud data into grids with the side length of L, taking the maximum height value of all points in the grids as the pixel value of the grids, and performing linear interpolation on the pixel value of the grids without points inside from the pixel value of the grids adjacent to the grids;
s3, vegetation classification: classifying the nDSM data, and extracting tree category data in the nDSM data to obtain a crown height model CHM;
s4, estimating shape parameters of the tree by the crown height model CHM, wherein the shape parameters of the tree comprise the crown height and the crown diameter;
s5, selection of ALS-DBH regression model: firstly, establishing an ALS-DBH regression model by an empirical model, then obtaining an even number of sampling groups by selecting a plurality of tree samples on site, taking one half of the sampling groups as model groups for model development, taking the other half of the sampling groups as verification groups, carrying out model development on all combination conditions, fitting through a cross verification process, and collecting parameters generated by each combination; by determining the coefficient R2Checking the predictive power and performance of the regression model and the accuracy of the predictions by RMSE of the predicted parameters, comparing R in all combination cases2And RMSE, selected to have a high R in model fitting2And an ALS-DBH regression model with low RMSE in the verification to predict DBH;
s6, prediction of carbon content:
according to all speciesCalculating total biomass, estimating the amount of carbon stored in the tree to be half of the total biomass, comparing the predicted amount of carbon stored with the amount of carbon stored estimated by the on-site measured DBH and the height, and comparing the estimated amount of carbon stored with RMSE and R2Evaluation was performed.
2. The urban tree carbon content estimation method based on multi-echo airborne laser scanning data according to claim 1, characterized in that: the wavelengths of the green laser, the near-infrared laser and the short-wave infrared laser in the step S1 are 532nm, 1064nm and 1550nm, respectively.
3. The urban tree carbon content estimation method based on multi-echo airborne laser scanning data according to claim 1, characterized in that: in the nsmd data obtained in step S24, if there is a height difference of 2m or more between the pixel value of a certain lattice and its eight neighbors, the height average of its eight neighbors is used as the pixel value of the lattice.
4. The urban tree carbon content estimation method based on multi-echo airborne laser scanning data according to claim 1, wherein step S3 specifically comprises:
s31, taking green laser intensity, near-infrared laser intensity, short-wave infrared laser intensity, nDSM, pNDWI and pNDVI as input, intercepting a training sample from input data, acquiring the class of the sample through aerial photos of Google Earth, using the training sample for training an SVM classifier, classifying nDSM data by using the trained SVM classifier, and extracting the class of trees in the nDSM data;
wherein the content of the first and second substances,
Figure FDA0002483551040000031
Figure FDA0002483551040000032
CGreen，CNIRand CSWIRGreen laser intensity, near infrared laser intensity, and short wave infrared laser intensity, respectively;
s32, filtering the extracted tree types for a plurality of times by adopting a 3-by-3 mode filter;
s33, empirically testing the filtering times to generate an optimal crown, which refers to a crown shape that can remove isolated pixels and retain dominance, and generating the crown height model CHM.
5. The urban tree carbon content estimation method based on multi-echo airborne laser scanning data according to claim 1, wherein step S4 specifically comprises:
s41, obtaining a local maximum pixel set: filtering the CHM by adopting a 3-by-3 local maximum filter, defining pixel points with pixel values higher than eight neighborhoods as local maximum values, and obtaining a local maximum pixel set;
s42, determining treetop: further extracting the local maximum pixel set by using ALS intensity data, and recording pixel points which are also local maximum in an ALS intensity image as treetop points;
s43, taking the detected treetop as a mark, adopting a region growing segmentation algorithm to segment the CHM into independent crowns, wherein each treetop point after segmentation has a closed cut block, and expressing the tree separation effect by absolute precision, the method comprises the following steps:
Figure FDA0002483551040000033
wherein n is1,1Is the detected number of crown blocks, ntotalIs the actual crown number; the crown height is defined as the mean of the local maximum in each crown cut, and the crown diameter is defined as the mean of the maximum crown diameter passing through the center of the highest point of the crown and the diameter in the vertical direction.
6. The urban tree carbon content estimation method based on multi-echo airborne laser scanning data according to claim 1, characterized in that: in step S6, the total biomass is the sum of the biomass in leaves, branches, wood and bark, if any,
Figure FDA0002483551040000034
Figure FDA0002483551040000035
Figure FDA0002483551040000041
Figure FDA0002483551040000042
ytotal＝yfoliage+ybranches+ywood+ybark；
then the process of the first step is carried out,
Figure FDA0002483551040000043
where y is the dry weight of the tree component in Kg, D is DBH in meters, H is the crown height, β is the model parameter, each tree component has a corresponding β, and e is the error compensation.
7. The urban tree carbon content estimation method based on multi-echo airborne laser scanning data according to claim 1, characterized in that: further comprising step S7;
s7, creating a city carbon storage map:
s71, converting the crown part derived from the ALS data and the carbon content stored in the crown part into vector data in the GIS;
s72, dividing urban land into a plurality of coverage types, and dividing the total carbon storage amount in each coverage type by the total area of the coverage type to obtain the carbon storage amount per unit area of each coverage type;
and S73, performing area division on urban land according to land coverage types, multiplying the total area of each area by the unit area carbon storage amount of the corresponding land coverage type, estimating the carbon content of each area, and finally creating an urban carbon storage map.
CN201711248101.0A 2017-12-01 2017-12-01 Urban tree carbon content estimation method based on multi-echo airborne laser scanning data Active CN108038433B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711248101.0A CN108038433B (en) 2017-12-01 2017-12-01 Urban tree carbon content estimation method based on multi-echo airborne laser scanning data

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711248101.0A CN108038433B (en) 2017-12-01 2017-12-01 Urban tree carbon content estimation method based on multi-echo airborne laser scanning data

Publications (2)

Publication Number Publication Date
CN108038433A CN108038433A (en) 2018-05-15
CN108038433B true CN108038433B (en) 2020-06-23

Family

ID=62094967

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711248101.0A Active CN108038433B (en) 2017-12-01 2017-12-01 Urban tree carbon content estimation method based on multi-echo airborne laser scanning data

Country Status (1)

Country Link
CN (1) CN108038433B (en)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109345530A (en) * 2018-10-08 2019-02-15 长安大学 A kind of quantitative evaluation method of all-aluminium piston carbon distribution cleaning effect
CN110147706B (en) * 2018-10-24 2022-04-12 腾讯科技（深圳）有限公司 Obstacle recognition method and device, storage medium, and electronic device
CN111144244B (en) * 2019-12-13 2023-07-11 唐美(北京)文化科技发展有限公司 Furniture wood category identification method and device
CN111337898B (en) * 2020-02-19 2022-10-14 北京百度网讯科技有限公司 Laser point cloud processing method, device, equipment and storage medium
CN113627492A (en) * 2021-07-20 2021-11-09 东软医疗系统股份有限公司 Method for determining size of scanning object, scanning method, scanning device and electronic equipment
CN116596198B (en) * 2023-07-18 2023-09-19 北京林业大学 Urban green land biomass calculation method, device, electronic equipment and storage medium

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103335608A (en) * 2013-07-03 2013-10-02 国家电网公司 Airborne LiDAR three-dimensional data acquisition method for establishing three-dimensional digital power transmission and transformation grid
CN104049245A (en) * 2014-06-13 2014-09-17 中原智慧城市设计研究院有限公司 Urban building change detection method based on LiDAR point cloud spatial difference analysis
CN104155638A (en) * 2014-06-11 2014-11-19 南京林业大学 Tree species classification method based on LiDAR (Light Detection and Ranging) false-vertical waveform model
CN104614729A (en) * 2014-11-20 2015-05-13 中国林业科学研究院资源信息研究所 Method for analyzing elevation matching quality of laser radar flight strip
CN104820830A (en) * 2015-05-08 2015-08-05 南京林业大学 Tree species identification method based on full-waveform LiDAR canopy profile model
CN104867180A (en) * 2015-05-28 2015-08-26 南京林业大学 UAV and LiDAR integrated forest stand characteristic inversion method
CN105913017A (en) * 2016-04-08 2016-08-31 南京林业大学 Corresponding period double high resolution remote sensing image-based forest biomass estimation method
CN106815847A (en) * 2017-01-12 2017-06-09 非凡智慧（宁夏）科技有限公司 Trees dividing method and single tree extracting method based on laser radar point cloud

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015075700A1 (en) * 2013-11-25 2015-05-28 First Resource Management Group Inc. Apparatus for and method of forest-inventory management

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103335608A (en) * 2013-07-03 2013-10-02 国家电网公司 Airborne LiDAR three-dimensional data acquisition method for establishing three-dimensional digital power transmission and transformation grid
CN104155638A (en) * 2014-06-11 2014-11-19 南京林业大学 Tree species classification method based on LiDAR (Light Detection and Ranging) false-vertical waveform model
CN104049245A (en) * 2014-06-13 2014-09-17 中原智慧城市设计研究院有限公司 Urban building change detection method based on LiDAR point cloud spatial difference analysis
CN104614729A (en) * 2014-11-20 2015-05-13 中国林业科学研究院资源信息研究所 Method for analyzing elevation matching quality of laser radar flight strip
CN104820830A (en) * 2015-05-08 2015-08-05 南京林业大学 Tree species identification method based on full-waveform LiDAR canopy profile model
CN104867180A (en) * 2015-05-28 2015-08-26 南京林业大学 UAV and LiDAR integrated forest stand characteristic inversion method
CN105913017A (en) * 2016-04-08 2016-08-31 南京林业大学 Corresponding period double high resolution remote sensing image-based forest biomass estimation method
CN106815847A (en) * 2017-01-12 2017-06-09 非凡智慧（宁夏）科技有限公司 Trees dividing method and single tree extracting method based on laser radar point cloud

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"从高光谱卫星数据中提取植被覆盖区铜污染信息";屈永华等;《光谱学与光谱分析》;20151130;第35卷(第11期);论文第3176-3180页 *
"联合机载LiDA和高分辨率影像估测落叶松单木生物量";赵晨阳;《万方数据》;20170405;论文第1-45页 *

Also Published As

Publication number Publication date
CN108038433A (en) 2018-05-15

Similar Documents

Publication Publication Date Title
CN108038433B (en) Urban tree carbon content estimation method based on multi-echo airborne laser scanning data
CN109190538B (en) Sediment-laden river delta coastal zone evolution analysis method based on remote sensing technology
Chen et al. Assessing tower flux footprint climatology and scaling between remotely sensed and eddy covariance measurements
Chen et al. Isolating individual trees in a savanna woodland using small footprint lidar data
Hou et al. Marine floating raft aquaculture extraction of hyperspectral remote sensing images based decision tree algorithm
Hester et al. Per-pixel classification of high spatial resolution satellite imagery for urban land-cover mapping
Lee et al. Forest vegetation classification and biomass estimation based on Landsat TM data in a mountainous region of west Japan
Ndehedehe et al. Satellite-derived changes in floodplain productivity and freshwater habitats in northern Australia (1991–2019)
Karna Mapping above ground carbon using worldview satellite image and lidar data in relationship with tree diversity of forests
Chen et al. Quantifying the carbon storage in urban trees using multispectral ALS data
Ambinakudige et al. Estimation of area and volume change in the glaciers of the Columbia Icefield, Canada using machine learning algorithms and Landsat images
Kolli et al. Assessment of change in the extent of mangrove ecosystems using different spectral indices in Google Earth Engine based on random forest model
Nelson et al. Spatial statistical techniques for aggregating point objects extracted from high spatial resolution remotely sensed imagery
Westin Quantification of a continuous-cover forest in Sweden using remote sensing techniques
Weng Remote sensing of urban biophysical environments
Hofierka et al. High-resolution urban greenery mapping for micro-climate modelling based on 3D city models
Li Dynamic monitoring algorithm of natural resources in scenic spots based on MODIS Remote Sensing technology
Li et al. Aerodynamic roughness length estimation with lidar and imaging spectroscopy in a shrub-dominated dryland
Huang et al. Comparing the effects of temporal features derived from synthetic time-series NDVI on fine land cover classification
Guyon et al. Applications of multispectral optical satellite imaging in forestry
Neogi et al. Change detection of exposed sandbars around Kaziranga national park
Wei et al. Long-term observation of global nuclear power plants thermal plumes using Landsat images and deep learning
Freitas et al. Evaluation of Shadow Effects on Remotely Sensed Reflectance from Permafrost Thaw Lakes in the Subarctic Forest-Tundra Zone
Maolin Dynamic monitoring algorithm of natural resources in scenic spots based on MODIS Remote Sensing technology
Reis Use of airborne laser scanning to improve selective logging and to assess size-frequency distribution of forest gaps in the Brazilian Amazon

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