CN112577954A - Urban green land biomass estimation method - Google Patents

Urban green land biomass estimation method Download PDF

Info

Publication number
CN112577954A
CN112577954A CN202011295789.XA CN202011295789A CN112577954A CN 112577954 A CN112577954 A CN 112577954A CN 202011295789 A CN202011295789 A CN 202011295789A CN 112577954 A CN112577954 A CN 112577954A
Authority
CN
China
Prior art keywords
biomass
area
image
data
formula
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
CN202011295789.XA
Other languages
Chinese (zh)
Other versions
CN112577954B (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 CN202011295789.XA priority Critical patent/CN112577954B/en
Publication of CN112577954A publication Critical patent/CN112577954A/en
Application granted granted Critical
Publication of CN112577954B publication Critical patent/CN112577954B/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/84—Systems specially adapted for particular applications
    • 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/10048—Infrared 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Theoretical Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Image Processing (AREA)

Abstract

The invention discloses a method for estimating urban greenbelt aboveground biomass, which comprises the steps of firstly taking panchromatic images and hyperspectral image data of a target area, segmenting the data after preprocessing, selecting an area containing vegetation in the data to carry out species identification and aboveground biomass inversion. The aboveground biomass of each plant can be accurately calculated, the calculation result effect is ideal, the accuracy is high, and compared with most of inversion models which do not identify plant species and are based on pixels, the method is quicker, more accurate and more comprehensive.

Description

Urban green land biomass estimation method
Technical Field
The invention relates to the field of measurement and calculation of aboveground biomass, in particular to a method for estimating the biomass on urban green land.
Background
The method has important significance for urban green land resource monitoring and environmental factor investigation. Meanwhile, the information can also be used for mastering the relationship between plants and the environment in the urban green land and the growth, development, updating and succession rules of vegetation in the urban green land, and has important significance for sustainable operation management of urban green land operation, carbon cycle research of an ecological system and understanding of global climate change. The conventional urban green land biomass extraction mainly depends on actual measurement or a statistical analysis method based on the actual measurement method, the precision is often not high, and the method is difficult to be practically popularized on the aspect.
At present, the method for estimating aboveground biomass by remote sensing images mainly comprises the following steps: 1) extracting vegetation indexes, ratios among wave bands and spectral curve characteristics related to the urban greenbelts from the optical remote sensing images to be used as independent variables, using national forest resource secondary survey data or sample area actual measurement data as dependent variables, and establishing a biomass regression model so as to estimate the aboveground biomass of the urban greenbelts, and 2) establishing a model based on the actual measurement data by using radar, laser radar, hyperspectral and hyperspectral images to perform inversion estimation on the aboveground biomass of the urban greenbelts. Both of the above two methods can roughly invert the aboveground biomass of urban greenbelt, but both have problems. In the first method, a plurality of plants are mixed, and only one model is established for inverting the aboveground biomass, so that the accuracy is low. In the second method, although the height information is increased by using radar data, the radar data volume is huge, the operation is complex and difficult to calculate, the estimated data volume calculation is increased, the intensity of the backscattering of the radar is linearly increased along with the increase of biomass, and after a certain biomass level is reached, the backscattering tends to be saturated, which can cause a certain influence on the estimation precision. Meanwhile, the two methods are used for carrying out inversion calculation on biomass on a pixel scale, cannot calculate the biomass for a single plant, and are low in precision.
Chinese patent document CN 107449400A records a forest aboveground biomass measurement system and a forest aboveground biomass measurement method, wherein an unmanned aerial vehicle aerial image is combined with ground actual measurement breast diameter data to establish a model, and an aboveground biomass of a forest is obtained by using a breast diameter-tree crown area regression equation and an aboveground biomass-tree crown area regression equation. Chinese patent document CN 108921885A records a method for jointly inverting forest aboveground biomass by integrating three types of data, preprocessing a high-resolution CCD image and a hyperspectral image, filtering laser radar point cloud data, interpolating to generate a digital terrain model, and normalizing the point cloud data. And then texture features, spectral features and point cloud structure features are extracted based on the three preprocessed data sources respectively, and a model is constructed by combining ground measured data and the extracted vegetation indexes and texture feature variables to predict forest ground biomass. Chinese patent document CN 104361592 a records an aboveground biomass estimation device, which mainly uses synthetic aperture radar SAR images to process, estimates tree height based on scattering data, obtains polygonal forest stands by object segmentation, and estimates aboveground biomass based on pixel reflectivity by establishing a multivariate linear model. However, in the method, trees and shrubs in the target area are modeled according to the same model to calculate the aboveground biomass, and grasslands and forest trees are not modeled to estimate the biomass, so that errors are brought to the measurement and calculation of the aboveground biomass. The requirement of the patent CN 107449400 that the forest species of the forest sample and the model sample are the same is not realistic in practical application and operation. The special CN 108921885A needs to use laser radar to measure the tree height, the data volume is extremely large, all trees and shrubs are mixed together to establish a unified model, and the species of the sample plot to be measured and the species of the forest in the modeling sample plot are also required to be consistent invisibly. None of these methods model the perspective of different plant species and make an estimation of aboveground biomass on the premise of species identification. Meanwhile, the urban green land is different from the urban forest, other construction lands in the green land can interfere with the measurement and calculation of the aboveground biomass, and the method is particularly important for accurately extracting and calculating the vegetation area.
Disclosure of Invention
The invention aims to solve the technical problem of providing an urban green land aboveground biomass estimation method, based on plant outline and species identification premises, the object-oriented urban green land aboveground biomass remote sensing inversion method can rapidly and accurately measure urban green land aboveground biomass and unit area biomass in a large area range, accurately calculate the aboveground biomass of each plant, and has important significance for urban green land resource monitoring and environmental factor investigation.
In order to solve the technical problems, the technical scheme adopted by the invention is as follows:
the method for estimating the biomass on the urban green land comprises the following steps:
acquiring panchromatic images and hyperspectral image data of a target area, and preprocessing the data;
secondly, segmenting the hyperspectral image, screening out vegetation areas in the image data, and forming vegetation image data of a target area, wherein the vegetation image data comprise all trees, shrubs and grasslands in vegetation of the target area;
thirdly, carrying out edge detection on the image of the plant region 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 the data to obtain pure end-member wave spectrum curves in the plant area images, and extracting and storing end-member wave spectrum characteristic curves of different plant species in the green land of a target area;
measuring the aboveground biomass of various plant species in the target area, establishing a relational model of crown widths of different species of trees and shrubs in the green area of the target area and the aboveground biomass of the trees and the shrubs in the green area, and establishing a relational model of the area and the biomass of the grassland;
classifying all data in the 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 of the panchromatic image, and identifying and defining the species category of each contour;
and step six, according to the relation model of the aboveground biomass corresponding to various plants established in the step four, combining the identification result in the step five, and obtaining the aboveground biomass and the aboveground biomass per hectare of the green land image data of the target area through calculation.
The preprocessing of the data in the first step comprises the following steps: and respectively carrying out radiometric calibration, atmospheric correction, geometric correction, mosaic and splicing on the panchromatic image and the hyperspectral image.
The image segmentation in the second step comprises the following specific processes:
and processing the normalized vegetation index NDVI by using a maximum between-class variance algorithm Otsu, and segmenting the hyperspectral image, namely dividing a vegetation area and a non-vegetation area in the image data to obtain target area image data, wherein the image data comprises all trees, shrubs and grasslands in the target area.
The formula for NDVI is:
NDVI＝(ρNIR-ρRED)/(ρNIR+ρRED)；
in the formula, ρNIRReflectivity in the near red band, pREDThe reflectance of the red band.
The formula for Otsu is:
let the gray level of NDVI image data be L (G ═ 1,2,3, …, L), and the number of pixels at gray level i be niThe total pixel number N of the NDVI is shown as follows:
Figure BDA0002785262380000031
the histogram is normalized to:
Figure BDA0002785262380000032
Figure BDA0002785262380000033
where p (i) represents the probability of the occurrence of a picture element with a 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 BDA0002785262380000041
Figure BDA0002785262380000042
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 for differentiating vegetation from non-vegetation in the NDVI image data.
The specific process of the third step is as follows:
carrying out edge detection based on a Sobel operator, and then carrying out vegetation contour identification, wherein the steps are as follows:
firstly, setting a Sobel operator, wherein the formula is as follows:
Figure BDA0002785262380000043
Figure BDA0002785262380000044
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 BDA0002785262380000045
Figure BDA0002785262380000046
finally, the horizontal and vertical gray scale of each pixel of the panchromatic image can be calculated, and the formula is as follows:
Figure BDA0002785262380000051
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;
on the basis of identifying the full-color image contour, taking the hyperspectral image of the plant species in each contour as a data set, and separating the hyperspectral data set G through a high-pass filter template to obtain noise GnSum signal GsThen, there are:
G＝Gn+GS
by maximizing the signal-to-noise ratio of the transformed data, i.e., maximizing the ratio of the signal covariance to the noise covariance, a transformation matrix V is obtained:
Figure BDA0002785262380000052
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 BDA0002785262380000053
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 to obtain a conversion matrix:
V＝[V1,V2,…,Vd]
extracting pure pixels in the hyperspectral image by PPI (maximum likelihood ratio) by using MNF (maximum likelihood ratio) results excluding noise bands to generate a pixel purity image, extracting the pure pixels by setting a threshold value, and calculating the pure pixels in the next step;
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 specific process of the step four is as follows:
survey and record the projection area L of crown of different plant species in urban green landjDiameter at breast height DjAnd above ground biomass YjH, k, b and a are constants, and the trees and shrubs Y of different plant species are establishedj、LjAnd DjRelation model ofType (2):
Lj＝k*Dj+b
Yj＝a*Dj h
the aboveground biomass of trees and shrubs is given by the above formula:
Figure BDA0002785262380000061
let the area of the lawn be RareaC is a constant and the above-ground biomass of the grass is Yg:
Yg＝c*Rarea
And the data classification in the fifth step adopts a spectrum angle classification method, wherein the spectrum angle classification is a classification method for classifying by comparing an unknown spectral line with a sample spectrum curve (an end member spectrum curve in the third step) in an n-dimensional space and selecting the most similar curve as an identification basis, namely a generalized cosine clip angle method. The method takes a spectrum curve as a vector of a spectrum space, calculates an included angle between two vectors, namely calculates an angle between a spectrum to be identified and a reference spectrum, and further determines the similarity of the two vectors, wherein a spectrum angle calculation formula is as follows:
Figure BDA0002785262380000062
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 the end member spectrum in the Pop library, and n is the number of the 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 concrete process of the step six is as follows:
acquiring the crown widths and species types of different arbor and shrub plant species and the area and the type of the grassland according to the identification result in the fifth step, and inverting the aboveground biomass in the area range of the target area by combining the relationship models of the crown widths and the biomass of different species established in the fourth step to acquire the aboveground biomass (U) in the range of the target areabio)：
Figure BDA0002785262380000071
In the formula of UbioIs the biomass in the region of interest, YijIs the biomass of strain i in plant species j, YgBiomass of grassland;
aboveground biomass (U) per unit area within the target areaavr)：
Figure BDA0002785262380000072
In the formula of UavrIs the biomass on the ground per unit area in a target area, UbioIs the biomass in the target region, UareaIs the total area within the target area.
The biomass estimation method further comprises the following steps:
and seventhly, establishing an inversion model to invert the biomass on the green land of the target area, and evaluating the fitting effect and precision of the model.
The inverse model described above takes the coefficient of determination R2The root mean square error RMSE and the relative root mean square error rRMSE evaluation model are calculated according to the following formula:
Figure BDA0002785262380000073
Figure BDA0002785262380000074
Figure BDA0002785262380000075
in the formula, xiThe measured value of the aboveground biomass in the target area is obtained;
Figure BDA0002785262380000076
the measured average value of the upper biomass in the target area is obtained;
Figure BDA0002785262380000077
a model estimate of aboveground biomass within the target region; n is the number of the same plots; i is the target region pattern number.
The invention provides an urban forest aboveground biomass estimation method, which is based on a high-resolution panchromatic image (DN value, textural features and spatial geometric features) and a hyperspectral image (spectral features) and is used for identifying plant outlines and species, and the object-oriented urban green land aboveground biomass remote sensing inversion method is based on hyperspectral data, high-resolution panchromatic remote sensing images and measured data as data sources, and the invention constructs a set of method for urban green land aboveground biomass inversion on the premise of vegetation canopy contour identification, hyperspectral end-member spectral curve feature extraction and species identification based on spectral angle comparison, and combines different species canopy and biomass aboveground models, firstly provides an object-oriented aboveground biomass remote sensing inversion method, can rapidly and accurately measure urban green land aboveground biomass and unit area biomass within a large area range, and each plant is obtained by accurate calculation, the aboveground biomass of each plant, the calculation result is ideal in effect and high in accuracy, and compared with most of pixel-based reverse models, the method is faster, more accurate and more comprehensive. The method provided by the patent has an ideal effect on estimating the biomass on the large-area urban forest land.
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 schematic diagram of a plant contour recognition result of a Sobel operator;
FIG. 4 is a graph of end-member spectral data for different vegetation types after MNF and PPI;
FIG. 5 is a schematic representation of the results of plant species identification based on spectral angle and contour detection;
FIG. 6 shows the identification results and basic information of plant species.
Detailed Description
As shown in fig. 1-6, a method for estimating biomass on urban green land comprises the following steps:
step one, acquiring high-resolution panchromatic images and hyperspectral image data of a target urban area, and preprocessing the data;
secondly, segmenting the hyperspectral image, screening out vegetation areas in the image data, and forming target city area image data which comprise all trees, shrubs and grasslands in the target area;
thirdly, carrying out edge detection on the image of the plant region by utilizing the high-resolution panchromatic 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 the data to obtain pure end-member wave spectrum curves in the images of the plant areas, and extracting and storing end-member wave spectrum characteristic curves of different plant species in the target urban areas;
measuring the aboveground biomass of various plant species in the target urban area, establishing a relation model of crown widths and aboveground biomass of different arbor and shrub species in the urban green area, and establishing a relation model of the area and the aboveground biomass for the grassland;
classifying all data in the 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 of the panchromatic image, and identifying and defining the species category of each contour;
and step six, according to the aboveground biomass relation models corresponding to the various plants established in the step four, combining the identification result in the step five, and obtaining the aboveground biomass in the target area and the aboveground biomass per hectare through calculation.
The preprocessing of the data in the first step comprises the following steps: respectively carrying out radiometric calibration, atmospheric correction, geometric correction, mosaic and splicing on the high-resolution panchromatic image and the hyperspectral image.
The image segmentation in the second step comprises the following specific processes:
and processing the normalized vegetation index NDVI by using a maximum between-class variance algorithm Otsu, and segmenting the hyperspectral image, namely dividing a vegetation area and a non-vegetation area in the image data.
The formula for NDVI is:
NDVI＝(ρNIR-ρRED)/(ρNIR+ρRED)；
in the formula, ρNIRReflectivity in the near red band, pREDThe reflectance of the red band.
The formula for Otsu is:
let the gray level of NDVI image data be L (G ═ 1,2,3, …, L), and the number of pixels at gray level i be niThe total pixel number N of the NDVI is shown as follows:
Figure BDA0002785262380000091
the histogram is normalized to:
Figure BDA0002785262380000092
Figure BDA0002785262380000093
where p (i) represents the probability of the occurrence of a picture element with a 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 BDA0002785262380000094
Figure BDA0002785262380000095
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 for differentiating vegetation from non-vegetation in the NDVI image data.
The specific process of the third step is as follows:
carrying out edge detection based on a Sobel operator, and then carrying out vegetation contour identification, wherein the steps are as follows:
firstly, setting a Sobel operator, wherein the formula is as follows:
Figure BDA0002785262380000101
Figure BDA0002785262380000102
secondly, performing plane convolution on a pixel set Z of the panchromatic high-resolution image and A and B to respectively realize the calculation of image gray values Gx and Gy detected by transverse and longitudinal edges, wherein the formula is as follows:
Figure BDA0002785262380000103
Figure BDA0002785262380000104
finally, the horizontal and vertical gray scale of each pixel of the high-resolution panchromatic image can be calculated, and the formula is as follows:
Figure BDA0002785262380000105
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;
on the basis of high-resolution panchromatic image contour identification, the hyperspectral images of plant species in each contour are used as data sets, and the hyperspectral data sets G are separated through a high-pass filter template to obtain noise GnSum signal GsThen, there are:
G＝Gn+GS
by maximizing the signal-to-noise ratio of the transformed data, i.e., maximizing the ratio of the signal covariance to the noise covariance, a transformation matrix V is obtained:
Figure BDA0002785262380000111
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 BDA0002785262380000112
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 to obtain a conversion matrix:
V＝[V1,V2,…,Vd]
extracting pure pixels in the hyperspectral image by PPI (maximum likelihood ratio) by using MNF (maximum likelihood ratio) results excluding noise bands to generate a pixel purity image, extracting the pure pixels by setting a threshold value, and calculating the pure pixels in the next step;
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 specific process of the step four is as follows:
survey and record the projection area L of the crown of different plant species in the urban green land of the target areajDiameter at breast height DjAnd above ground biomass YjH, k, b and a are constants, and the trees and shrubs Y of different plant species are establishedj、LjAnd DjThe relationship model of (1):
Lj＝k*Dj+b
Yj＝a*Dj h
the aboveground biomass of trees and shrubs is given by the above formula:
Figure BDA0002785262380000121
let the area of the lawn be RareaC is a constant and the above-ground biomass of the grass is Yg:
Yg＝c*Rarea
And the data classification in the fifth step adopts a spectrum angle classification method, wherein the spectrum angle classification is a classification method for classifying by comparing an unknown spectral line with a sample spectrum curve (an end member spectrum curve in the third step) in an n-dimensional space and selecting the most similar curve as an identification basis, namely a generalized cosine clip angle method. The method takes a spectrum curve as a vector of a spectrum space, calculates an included angle between two vectors, namely calculates an angle between a spectrum to be identified and a reference spectrum, and further determines the similarity of the two vectors, wherein a spectrum angle calculation formula is as follows:
Figure BDA0002785262380000122
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 concrete process of the step six is as follows:
acquiring the crown widths and species categories of different arbor and shrub plant species and the areas and types of grasslands according to the identification result in the fifth step in combination with the different established in the fourth stepThe relation model of the species crown width and biomass is used for inverting the aboveground biomass in the area range of the target urban area to obtain the aboveground biomass (U) in the range of the target areabio)：
Figure BDA0002785262380000123
In the formula of UbioIs the biomass of the target region, YijIs the biomass of strain i in species j, YgIs the biomass of the grassland;
aboveground biomass (U) per unit area within a target urban areaaqr)：
Figure BDA0002785262380000124
In the formula of UavrIs the biomass on the ground per unit area in a target area, UbioIs the biomass in the target region, UareaIs the total area within the target area.
The biomass estimation method further comprises the following steps:
and seventhly, establishing an inversion model to invert the aboveground biomass of the target area, and evaluating the fitting effect and precision of the model.
The inverse model described above takes the coefficient of determination R2The root mean square error RMSE and the relative root mean square error rRMSE evaluation model are calculated according to the following formula:
Figure BDA0002785262380000131
Figure BDA0002785262380000132
Figure BDA0002785262380000133
in the formula, xiThe measured value of the aboveground biomass in the target area is obtained;
Figure BDA0002785262380000134
the measured average value of the upper biomass in the target area is obtained;
Figure BDA0002785262380000135
a model estimate of aboveground biomass within the target region; n is the number of the same plots; i is the target region pattern number.
Example 1:
1. the test area of the embodiment is located in Shanghai Fengxian area Wu bridge nursery base (121.41 degrees E, 30.96 degrees N), and the area is about 3.83hm2The vegetation types are rich. Wherein the plants comprise: 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 and multispectral data by means of an unmanned aerial vehicle;
3. 28 square patterns (10X 10 m) were set in the area of the study2) Similarly, the center point coordinates are determined using GPS (Trimble GeoXH6000), which is located by receiving wide area differential signals with an accuracy better than 0.5 m. Recording tree species in the sample plot, counting, measuring crown width, breast diameter and tree height of each tree, and using above-ground biomass as a formulaCalculating by combining the breast height and the tree height, summarizing according to survey data, establishing a relational model of crowns and aboveground biomass of different tree species, and establishing a relation between grassland and aboveground biomass, as shown in table 1:
TABLE 1
Figure BDA0002785262380000141
Note: wherein the projection area of the crown is LjThe grassland area is Rarea。
4. As shown in fig. 1, during data preprocessing, the high-resolution panchromatic images are stitched to generate a complete high-resolution image of the study area. Meanwhile, geometric fine correction is carried out on the image by combining ground actual measurement control point data and adopting a quadratic polynomial model. And carrying out radiometric calibration on the original hyperspectral data by means of sensor radiometric calibration data, and carrying out atmospheric correction by combining an empirical linear model with the spectral data of the ground actual measurement target. Meanwhile, carrying out geometric fine correction on the hyperspectral image by utilizing a ground actual measurement control point;
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. and (5) 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. As shown in fig. 3, Sobel edge detection is performed on the vegetation distribution area of the panchromatic image, and crown contours and grassland contours of different plant species are identified;
8. as shown in fig. 4, on the basis of identifying the outline of the plant species, the plant species which has been marked in the sample plot of the field investigation is marked, the hyperspectral image in the marked outline is processed, including minimum noise separation and pure pixel extraction, the end-member spectral curves of different typical species are extracted and recorded into a warehouse;
9. as shown in fig. 5, 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 is identified by using a method of spectrum angle comparison;
10. as shown in fig. 5 and 6, the outline area of the plant species is calculated, and the identification result is combined to obtain the crown area of trees and shrubs and the distribution area of the grassland, and in this example, there are 23 different plant species, and the total crown area 2028 is obtained, and the grassland 125 blocks are obtained. Meadow area 1019.19m2Area of trees and shrubs 12735.29 m2；
11. As shown in fig. 6, the 23 relation models in step 3 are substituted, calculation is performed to obtain the biomass on the urban green land and the biomass per unit area in the image area, and the biomass on the land of each plant is obtained;
12. using a coefficient of determination (R)2) The effect of model fitting and the accuracy of the estimation were evaluated as Root Mean Square Error (RMSE) and relative Root Mean Square Error (RMSE), as shown in table 2:
TABLE 2 summary table of estimation accuracy of aboveground biomass models
Figure BDA0002785262380000151
In the embodiment, high-resolution panchromatic image data and hyperspectral image data which are acquired simultaneously are combined, a vegetation area is extracted, plant species outlines including tree crown outlines of trees and shrubs and an outline of a grassland are identified, an urban green land aboveground biomass estimation model is established based on extracted characteristic variables, a data source records green land structure information from the aspects of spectrum and space, the data are complementary to each other, the scheme is object-oriented aboveground biomass estimation, the traditional pixel-based inversion method is eliminated, different plant species can be identified, models are established for different plant species, and aboveground biomass is measured. Therefore, the method enhances the capability and the precision of biomass inversion on the urban green land; the data verification results in tables 1 and 2 show that compared with the results of aboveground biomass estimation by other similar remote sensing methods, the above-ground biomass extraction method based on the embodiment has the advantage that the relative root mean square error is reduced by more than 5%.

Claims (9)

1. The method for estimating the biomass on the urban green land is characterized by comprising the following steps of:
acquiring panchromatic images and hyperspectral image data of a target area, and preprocessing the data;
secondly, segmenting the hyperspectral image, screening out vegetation areas in the image data, and forming vegetation image data of a target area, wherein the vegetation image data comprise all trees, shrubs and grasslands in vegetation of the target area;
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-member wave spectrum curves in plant area images, and extracting end-member wave spectrum characteristic curves of different plant species in green lands of a target area and storing the end-member wave spectrum characteristic curves in a warehouse;
measuring the aboveground biomass of various plant species in the target area, establishing a relational model of crown widths of different species of trees and shrubs and the aboveground biomass of the same in the green area of the target area, and establishing a relational model of the elevation area of the grassland and the aboveground biomass;
classifying all data in the 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 of the panchromatic image, and identifying and defining the species category of each contour;
and step six, according to the aboveground biomass relation models corresponding to the various plants established in the step four and the identification result in the step five, the aboveground biomass and the aboveground biomass per hectare of the green land image data of the target area are obtained through calculation.
2. The method for estimating biomass on urban greenfield according to claim 1, wherein the preprocessing of the data in the first step comprises: respectively carrying out radiometric calibration on the panchromatic image and the hyperspectral image, and carrying out atmospheric correction, geometric correction, mosaic and splicing on image data.
3. The method for estimating biomass on urban greenbelt according to claim 2, wherein the image segmentation in the second step comprises the following specific processes:
and processing the normalized vegetation index NDVI by using a maximum between-class variance algorithm Otsu, and segmenting the hyperspectral image, namely dividing a vegetation area and a non-vegetation area in the image data.
The formula for NDVI is:
NDVI＝(ρNIR-ρRED)/(ρNIR+ρRED)；
in the formula, ρNIRReflectivity in the near red band, pREDThe reflectance of the red band.
The formula for Otsu is:
let the gray level of NDVI image data be L (G ═ 1,2,3, …, L), and the number of pixels at gray level i be niThe total pixel number N of the NDVI is shown as follows:
Figure FDA0002785262370000021
the histogram is normalized to:
Figure FDA0002785262370000022
Figure FDA0002785262370000023
where p (i) represents the probability of the occurrence of a picture element with a 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 FDA0002785262370000024
Figure FDA0002785262370000025
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 for differentiating vegetation from non-vegetation in the NDVI image data.
4. The method for estimating biomass on urban green land according to claim 3, wherein the concrete process of the third step is as follows:
carrying out edge detection based on a Sobel operator, and then carrying out vegetation contour identification, wherein the steps are as follows:
firstly, setting a Sobel operator, wherein the formula is as follows:
Figure FDA0002785262370000031
Figure FDA0002785262370000032
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 FDA0002785262370000033
Figure FDA0002785262370000034
finally, the horizontal and vertical gray scale of each pixel of the panchromatic image can be calculated, and the formula is as follows:
Figure FDA0002785262370000035
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;
on the basis of identifying the full-color image contour, taking the hyperspectral image of the plant species in each contour as a data set, and separating the hyperspectral data set G through a high-pass filter template to obtain noise GnSum signal GsThen, there are:
G＝Gn+Gs
by maximizing the signal-to-noise ratio of the transformed data, i.e., maximizing the ratio of the signal covariance to the noise covariance, a transformation matrix V is obtained:
Figure FDA0002785262370000036
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 FDA0002785262370000041
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 PPI (maximum likelihood ratio) by using MNF (maximum likelihood ratio) results excluding noise bands 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.
5. The method for estimating biomass on urban greenbelt according to claim 4, wherein the specific process of step four is as follows:
survey and record the projection area L of crown of different plant species in urban green landjDiameter at breast height DjOn the groundBiomass YjH, k, b and a are constants, and the trees and shrubs Y of different plant species are establishedj、LjAnd DjThe relationship model of (1):
Lj＝k*Dj+b
Yj＝a*Dj h
the aboveground biomass of trees and shrubs is given by the above formula:
Figure FDA0002785262370000042
let the area of the lawn be RareaC is a constant and the above-ground biomass of the grass is Yg：
Yg＝c*Rarea。
6. The method according to claim 5, wherein the data classification in step five is performed by using a spectral angle classification method.
7. The method for estimating biomass on urban green land according to claim 1, wherein the concrete process of the sixth step is as follows:
acquiring the crown widths and species types of different arbor and shrub plant species and the area and the type of the grassland according to the identification result in the fifth step, and inverting the aboveground biomass in the area range of the target area by combining the relationship models of the crown widths and the biomass of different species established in the fourth step to acquire the aboveground biomass (U) in the range of the target areabio)：
Figure FDA0002785262370000051
In the formula of UbioIs the biomass in the region of interest, YijIs the biomass of strain i in plant species j, YgBiomass of grassland;
aboveground biomass (U) per unit area within the target areaavr)：
Figure FDA0002785262370000052
In the formula of UavrIs the biomass on the ground per unit area in a target area, UbioIs the biomass in the target region, UareaIs the total area within the target area.
8. The method for estimating biomass on urban greenfield according to claim 1, further comprising:
and seventhly, establishing an inversion model to invert the biomass on the green land of the target area, and evaluating the fitting effect and precision of the model.
9. The method for estimating biomass on urban greenfield according to claim 8, wherein: the inverse model adopts a decision coefficient R2The root mean square error RMSE and the relative root mean square error rRMSE evaluation model are calculated according to the following formula:
Figure FDA0002785262370000053
Figure FDA0002785262370000054
Figure FDA0002785262370000055
in the formula, xiThe measured value of the aboveground biomass in the target area is obtained;
Figure FDA0002785262370000056
is in the target areaActually measured average value of aboveground biomass;
Figure FDA0002785262370000057
a model estimate of aboveground biomass within the target region; n is the number of the same plots; i is the target area pattern number.
CN202011295789.XA 2020-11-18 2020-11-18 Urban green land biomass estimation method Active CN112577954B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011295789.XA CN112577954B (en) 2020-11-18 2020-11-18 Urban green land biomass estimation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011295789.XA CN112577954B (en) 2020-11-18 2020-11-18 Urban green land biomass estimation method

Publications (2)

Publication Number Publication Date
CN112577954A true CN112577954A (en) 2021-03-30
CN112577954B CN112577954B (en) 2023-08-04

Family

ID=75122949

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011295789.XA Active CN112577954B (en) 2020-11-18 2020-11-18 Urban green land biomass estimation method

Country Status (1)

Country Link
CN (1) CN112577954B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117372503A (en) * 2023-12-08 2024-01-09 水利部交通运输部国家能源局南京水利科学研究院 River and lake shore zone vegetation classification and coverage calculation method and system

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106291582A (en) * 2016-09-28 2017-01-04 中国科学院华南植物园 A kind of divide different forest biomass remote sensing inversion method based on curve of spectrum feature
CN108921885A (en) * 2018-08-03 2018-11-30 南京林业大学 A kind of method of comprehensive three classes data source joint inversion forest ground biomass

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106291582A (en) * 2016-09-28 2017-01-04 中国科学院华南植物园 A kind of divide different forest biomass remote sensing inversion method based on curve of spectrum feature
CN108921885A (en) * 2018-08-03 2018-11-30 南京林业大学 A kind of method of comprehensive three classes data source joint inversion forest ground biomass

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
周利鹏 等: "《基于高光谱影像的舟曲曲瓦沟树种识别》" *
张加龙 等: "《基于遥感的高山松连清固定样地地上生物量估测模型构建》" *
林勇 等: "《高光谱遥感技术在城市绿地调查中的应用及发展趋势》" *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117372503A (en) * 2023-12-08 2024-01-09 水利部交通运输部国家能源局南京水利科学研究院 River and lake shore zone vegetation classification and coverage calculation method and system
CN117372503B (en) * 2023-12-08 2024-03-08 水利部交通运输部国家能源局南京水利科学研究院 River and lake shore zone vegetation classification and coverage calculation method and system

Also Published As

Publication number Publication date
CN112577954B (en) 2023-08-04

Similar Documents

Publication Publication Date Title
CN112541921B (en) Urban green land vegetation information data accurate determination method
CN109146889B (en) Farmland boundary extraction method based on high-resolution remote sensing image
Halme et al. Utility of hyperspectral compared to multispectral remote sensing data in estimating forest biomass and structure variables in Finnish boreal forest
US8111924B2 (en) Remote sensing and probabilistic sampling based method for determining the carbon dioxide volume of a forest
US7639842B2 (en) Remote sensing and probabilistic sampling based forest inventory method
CN110363246B (en) Fusion method of vegetation index NDVI with high space-time resolution
CN112669363B (en) Method for measuring three-dimensional green space of urban green space
CN114821362B (en) Multi-source data-based rice planting area extraction method
CN112577906B (en) Urban green land soil moisture content detection method
Zhang et al. Extraction of tree crowns damaged by Dendrolimus tabulaeformis Tsai et Liu via spectral-spatial classification using UAV-based hyperspectral images
CN111582194A (en) Multi-temporal high-resolution remote sensing image building extraction method based on multi-feature LSTM network
CN112699756B (en) Hyperspectral image-based tea origin identification method and system
CN108898070A (en) A kind of high-spectrum remote-sensing extraction Mikania micrantha device and method based on unmanned aerial vehicle platform
CN110765977A (en) Method for extracting wheat lodging information based on multi-temporal remote sensing data of unmanned aerial vehicle
CN113822141A (en) Automatic glacier and snow extraction method and system based on remote sensing image
CN112577954B (en) Urban green land biomass estimation method
CN116206211A (en) Feature type identification method for integrating time sequence remote sensing information and KL-divergence
Yeom et al. Separability analysis and classification of rice fields using KOMPSAT-2 High Resolution Satellite Imagery
Danoedoro et al. Combining pan-sharpening and forest cover density transformation methods for vegetation mapping using Landsat-8 Satellite Imagery
Kai et al. Effects of resampling image methods in sugarcane classification and the potential use of vegetation indices related to chlorophyll
CN112651295A (en) Urban green land tree identification system and method
CN112577907B (en) Urban green land tree crown loss rate calculation method
Huang et al. Information fusion approach for biomass estimation in a plateau mountainous forest using a synergistic system comprising UAS-based digital camera and LiDAR
Ghiyamat et al. Influence of tree species complexity on discrimination performance of vegetation indices
Iskandar et al. Crown closure segmentation on wetland lowland forest using the mean shift algorithm

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