Urban green infrastructure carbon storage estimation method and system integrating multi-source data

Abstract

The present invention discloses a method and system for estimating carbon reserves of urban green infrastructure by integrating multi-source data, including performing stereo matching calculation on multispectral images of the study area to generate dense point clouds; classifying the airborne laser radar point clouds and the dense point clouds generated by multispectral images of the study area to obtain ground points and canopy points; obtaining CHM images based on the DEM generated by the ground points and the DSM generated by the canopy points; fusing the CHM images generated by the laser radar in the same area and the CHM images generated by the multispectral images; extracting tree height parameters based on the CHM fused images in combination with edge detection and watershed segmentation algorithms; calculating the carbon reserves of urban green infrastructure according to the allometric growth equations and carbon content of each dominant tree species; and calculating the human activity correction coefficient to correct the carbon reserves. The estimation results obtained by the method proposed in the present invention have high precision and accuracy, and provide scientific data support and decision-making suggestions for urban "energy conservation and emission reduction" and achieving my country's "carbon neutrality" goals.

Description

Technical Field

[0001] The present invention relates to the fields of remote sensing, geographic information and ecological technology, and in particular to a method and system for estimating carbon reserves of urban green infrastructure by integrating multi-source data. Background Art

[0002] Rapid urbanization and socioeconomic development have led to a series of ecological and environmental problems. Maintaining the carbon balance of the urban environment has become a major challenge in responding to human-induced climate change. Urban green infrastructure is a concept related to urban natural ecosystems and green spaces, including various natural elements in urban areas, such as greenways, parks, gardens, green roofs, woodlands, waterways, community farms, forests and wilderness, which play a key role in mitigating urban carbon emissions, regulating urban climate conditions, and improving air purification. Therefore, accurate assessment of urban green infrastructure carbon storage is of great significance for urban management policy formulation, urban green infrastructure construction and protection, and guiding urban green development strategies.

[0003] Carbon storage can usually be calculated by building a carbon storage estimation model based on collected ground measurement data (tree diameter, height, biomass, etc.) or remote sensing data (satellite or aerial images). However, the urban green infrastructure environment is complex, with few homogeneous surface types and pure vegetation cover, making it difficult to accurately estimate using a single data source. At present, commonly used data sources for carbon storage estimation often come from visible light sensors, thermal infrared sensors, and microwave sensors. Based on multispectral image data obtained by visible light sensors, the required characteristic factors (such as vegetation index, terrain characteristics, texture characteristics, etc.) can be extracted and combined with forest structure attributes and spectral information to construct a carbon storage estimation model. Microwave-based sensors (such as synthetic aperture radar (SAR) technology) estimate aboveground carbon storage based on their bands (L band and C band) and polarization (HH, HV, VV and VH). However, satellite remote sensing data are easily affected by atmospheric conditions (clouds, precipitation, fog, etc.) and lack the necessary accuracy in complex urban landscapes; synthetic aperture radar (SAR) based on microwave remote sensing can penetrate clouds for all-weather imaging, but it still faces challenges in terms of signal saturation and forest structure complexity. The emergence of laser radar technology (LiDAR) has greatly improved the ability to generate spatial and temporal data. As an active remote sensing system, it can obtain accurate vegetation canopy structure characteristics and significantly improve the assessment of aboveground carbon storage in urban trees. However, LiDAR has large processing and data volumes, and uneven density may lead to a lack of detailed information in some areas. Therefore, when estimating the carbon storage of urban green infrastructure, it is crucial to select data sources or use multi-source remote sensing data to improve the estimation accuracy of urban green infrastructure carbon storage. Summary of the invention

[0004] In view of the shortcomings of the prior art, the present invention provides a method for estimating carbon storage of urban green infrastructure by integrating multi-source data, comprising the following steps:

[0005] Step 1: measure the tree height and tree species parameters of the sample points in the study area, and obtain the airborne lidar point cloud data and multispectral images of the study area;

[0006] Step 2, preprocessing the airborne laser radar point cloud data and multispectral image obtained in step 1;

[0007] Step 2.1, perform radiation correction, atmospheric correction and geometric correction on the multispectral image, and then perform stereo matching to calculate the three-dimensional coordinates to generate a dense point cloud;

[0008] Step 2.2, using the airborne laser radar point cloud data in step 1 and the dense point cloud in step 2.1 as two point cloud data sources, using the adaptive hybrid model classification method to classify the two point cloud data sources respectively, and obtain the ground points and canopy points of the two point cloud data sources;

[0009] Step 3, based on the ground points obtained by point cloud classification in step 2.2, digital elevation models DEM of two point cloud data sources are generated, based on the canopy points obtained by point cloud classification in step 2.2, digital surface models DSM of two point cloud data sources are generated, and the DSM of the two point cloud data sources are respectively subtracted from the DEM to obtain the canopy height models CHM corresponding to the two point cloud data sources, and then the CHM images corresponding to the two point cloud data sources are obtained;

[0010] Step 4, Fourier transform is used to fuse the CHM images corresponding to the two point cloud data sources in the same area to obtain a CHM fused image;

[0011] Step 5, based on the CHM fusion image, the tree height parameters are extracted by combining edge detection and watershed segmentation algorithm;

[0012] Step 5.1, calculate the gradient magnitude and gradient direction of the CHM fusion image;

[0013] Step 5.2, set the initial water source point and assign a value to the water level depth of the initial water source point;

[0014] Step 5.3, updating the water level depth of the pixel on the CHM fusion image;

[0015] Step 5.4, preliminarily perform watershed segmentation on the CHM fusion image;

[0016] Step 5.5, based on the watershed segmentation result of step 5.4, add a new water source point;

[0017] Step 5.6, perform the operations of steps 5.3-5.4 on the initial water source points set in step 5.2 and the water source points added in step 5.5 to obtain a new watershed segmentation result;

[0018] Step 5.7, find the local maximum gray value in each area obtained by watershed segmentation in step 5.6, and use it as the tree height;

[0019] Step 6: Substitute the tree height parameters of the study area extracted in step 5 into the allometric growth equation of the corresponding tree species to obtain the biomass, and then multiply it by the carbon content of each tree species to obtain the carbon storage of a single tree. Count the carbon storage of all single trees to obtain the preliminary carbon storage of the urban green infrastructure in the study area.

[0020] Step 7: Calculate the human activity correction coefficient of the study area in combination with the quantitative indicators, and correct the preliminary carbon storage obtained in step 6 to obtain the final carbon storage of urban green infrastructure in the study area.

[0021] Furthermore, in step 1, for each sample point, the GPS position, tree height and tree species parameters are recorded.

[0022] Furthermore, the dense point cloud coordinate calculation formula in step 2.1 is as follows:

[0023] (1)

[0024] (2)

[0025] (3)

[0026] In the formula, is the camera focal length, is the stereo baseline, is parallax, To generate the 3D coordinates of points in a dense point cloud, is the coordinate of the principal point of the multispectral image, are the pixel coordinates of the multispectral image.

[0027] Furthermore, the classification process in step 2.2 is specifically performed as follows:

[0028] ① Use the existing high-precision point cloud data to train the regression model so that it can learn the relationship between the surface complexity parameters of the object and the number of triangles in the TIN-GTP-TEN hybrid model;

[0029] ② Input the airborne laser radar point cloud data into the trained regression model. The regression model will automatically select points in the airborne laser radar point cloud data to generate the optimal TIN-GTP-TEN hybrid model, traverse the points to be classified in the airborne laser radar point cloud data that are not involved in the construction of the optimal TIN-GTP-TEN hybrid model, calculate the distance between each point to be classified and its projection point vertically projected onto the triangular surface in the optimal TIN-GTP-TEN hybrid model, calculate the maximum value of the angle between the vector from the point to be classified to its vertex on the triangular surface in the optimal TIN-GTP-TEN hybrid model and the normal vector of the triangular surface, and compare it with the set distance threshold and angle threshold to determine whether the point to be classified is a ground point or a canopy point. The calculation formula of the classification algorithm is as follows:

[0030] (4)

[0031] In the formula, Indicate point P The classification results, Tree represents the canopy point, Ground represents a ground point, Yes P To its projection point vertically projected onto the triangle surface in the optimal TIN-GTP-TEN hybrid model The distance between Yes P To the maximum value of the angle between the vector perpendicularly projected to the vertex of the triangle surface in the optimal TIN-GTP-TEN hybrid model and the normal vector of the triangle surface, Indicate point P Vertically project triangles in the optimal TIN-GTP-TEN hybrid model, is the set distance threshold, is the set angle threshold;

[0032] Perform the above calculation and judgment operations on all points to be classified to complete the classification;

[0033] ③Use the processing method in ② to classify the dense point cloud generated by multispectral images.

[0034] The TIN-GTP-TEN hybrid model is constructed as follows: first, the surface features of the study area are modeled using the existing improved Delaunay triangulation algorithm to generate a TIN surface model, and then the internal space below the TIN surface is divided using GTP. Then, within each GTP unit, the tetrahedron partitioning algorithm is used to decompose each heterogeneous GTP into smaller, homogeneous TENs to ensure the accurate expression of the surface features of the study area.

[0035] Furthermore, the specific process of image fusion using Fourier transform in step 4 is as follows:

[0036] ① Perform Fourier transform on the CHM image generated by airborne lidar point cloud data and the CHM image generated by multispectral image to obtain its frequency domain image. The calculation formula is:

[0037] (6)

[0038] In the formula, represents a CHM image, is the spatial domain coordinate, represents the frequency domain image of the CHM image, is the frequency domain coordinate, represents the imaginary unit, ;

[0039] ② Use filters to separate the low-frequency and high-frequency parts of the frequency domain image;

[0040] Low pass filter The definition is as follows:

[0041] (7)

[0042] In the formula, is the horizontal cutoff frequency in the frequency domain, is the cutoff frequency in the vertical direction in the frequency domain, is the frequency domain coordinate;

[0043] High Pass Filter The definition is as follows:

[0044] (8)

[0045] The low-pass filter and high pass filter With frequency domain image Multiply and separate the desired frequency components, namely:

[0046] (9)

[0047] (10)

[0048] In the formula, represents the frequency domain image of the CHM image, represents the low-frequency image of the CHM image, Represents the high-frequency image of the CHM image;

[0049] ③ The high-frequency and low-frequency parts of two CHM images in the same area of ​​the study area are fused respectively;

[0050] For low-frequency components, the weighted average method is used for fusion. The specific calculation formula is:

[0051] (11)

[0052] In the formula, represents the fused low-frequency image, A low-frequency image representing the CHM image generated from airborne lidar point cloud data, Represents the low-frequency image of the CHM image generated by multispectral imaging, is a weight coefficient between 0 and 1;

[0053] For high-frequency components, the maximum value method is used for fusion. The specific calculation formula is:

[0054] (12)

[0055] In the formula, represents the fused high-frequency image, Indicates taking the maximum value, A high-frequency image representing the CHM image generated from LiDAR data, High-frequency image representing the CHM image generated by multispectral imaging;

[0056] ④ Inverse Fourier transform, merge the fused low-frequency and high-frequency components, and transform them back to the spatial domain through inverse Fourier transform to obtain the CHM fusion image , the specific calculation method is as follows:

[0057] (13)

[0058] In the formula, is the spatial domain coordinate, represents the fused low-frequency image, represents the fused high-frequency image, is the frequency domain coordinate, represents the imaginary unit, .

[0059] Furthermore, the gray value matrix of the CHM fusion image in step 5.1 is , calculate the gradient magnitude and the gradient direction , the specific calculation formula is as follows:

[0060] (14)

[0061] (15)

[0062] In the formula, and They are the gray value matrices exist x Direction and y Partial derivatives in the direction.

[0063] Furthermore, in step 5.2, the Canny edge detection algorithm is used to extract the edges in the CHM fusion image and preliminarily identify the outline of the tree crown. The points in the edge image with pixel grayscale values ​​lower than the surrounding 8 neighboring pixels are taken as the initial water source points, and the water level depth of these initial water source points is calculated. Set to the grayscale value of this point.

[0064] Furthermore, in step 5.3, the gradient direction Iterative simulation updates the convergence of water flow: ① Starting from any water source point, find the gradient direction The next pixel , if the next pixel Water level depth Greater than , then update the next pixel Water level depth Otherwise, do not update ; ② Then from the pixel point Start by looking for the gradient direction The next pixel , if the next pixel Water level depth Greater than , then update the next pixel Water level depth Otherwise, do not update ; ③ Repeat ①-② until the water level depth of all pixels is no longer updated.

[0065] Furthermore, in step 5.4, when the pixel Water level depth Equal to its gray value in DEM When , the pixel point is considered to be the boundary point of different watersheds. Starting from any boundary point, the next boundary point is tracked along the direction of gradient change until the tracked boundary points are connected to form a closed area or the image edge has been tracked. The tracking is ended and all the tracked boundary points are connected in sequence to form the boundary line of the watershed. The above operation is repeated for each boundary point to form multiple segmentation areas, and the watershed segmentation of the CHM fusion image is realized.

[0066] Furthermore, in the step 5.5, the area segmented in step 5.4 is used as the statistical range for determining the average grayscale value or grayscale variance, and the water source point is adjusted based on this range; if the average grayscale value of a certain area after segmentation is significantly lower than the average of other surrounding areas, the pixel points of this area are added as new water source points; if the grayscale variance of a certain area after segmentation is significantly higher than that of other surrounding areas, the pixel points of this area are added as new water source points.

[0067] Furthermore, step 6 includes: finding the allometric growth equation and carbon content applicable to each dominant tree species in the study area in published literature, inputting the tree height into the allometric growth equation to calculate the biomass, and then multiplying it by the corresponding carbon content to obtain the preliminary carbon storage.

[0068] Furthermore, the human activity correction coefficient in step 7 is The calculation formula is as follows:

[0069] (16)

[0070] In the formula, is the human activity factor, obtained by multiplying activity intensity and activity frequency; Indicates different impact factor types, Indicates trimming, Indicates fertilization. It means destruction; S is the urban space factor, obtained from the urbanization rate of the study area; E It is an ecological sensitivity factor obtained through ecological sensitivity analysis and is divided into four sensitivity levels: highly sensitive , generally sensitive , low sensitivity and insensitivity ; A is the social factor, obtained by the ratio of the population of the study area to the area of ​​the study area;

[0071] Multiply the preliminary carbon storage obtained in step 6 by the human activity correction coefficient to obtain the final carbon storage of the study area.

[0072] The present invention also provides a system for estimating carbon reserves of urban green infrastructure by integrating multi-source data, which is used to implement the method for estimating carbon reserves of urban green infrastructure by integrating multi-source data as described above.

[0073] Furthermore, it includes a processor and a memory, the memory is used to store program instructions, and the processor is used to call the stored instructions in the memory to execute the urban green infrastructure carbon storage estimation method that integrates multi-source data as described above.

[0074] Compared with the prior art, the present invention has the following advantages:

[0075] 1) An adaptive hybrid model classification method is proposed, and a TIN-GTP-TEN hybrid model is constructed using machine learning technology to replace the traditional TIN model for point cloud classification, which can more finely represent the changes in complex terrain and effectively improve the accuracy of point cloud classification; 2) The frequency domain fusion technology is used to fuse two types of remote sensing data (multispectral remote sensing images and LiDAR data) from different sensors, effectively making up for the limitation of relying on a single data source; 3) The crown segmentation of the CHM fusion image is performed by combining edge detection with the watershed segmentation algorithm, making the crown edge extraction more accurate; 4) Different The differences in tree species and the impact of human activities were investigated and verified through field surveys to determine the distribution of each dominant tree species in the study area. The carbon storage was estimated based on the allometric growth equation and carbon content of each dominant tree species. Combined with quantitative indicators (human activity factors, urban space factors, ecological sensitivity factors, and social factors), a method for calculating the human activity correction coefficient was proposed to correct the carbon storage. 5) The estimation results obtained by the method proposed in the present invention have high precision and accuracy, which provides a scientific basis for urban greening planning and ecological management, and also provides scientific data support and decision-making recommendations for urban "energy conservation and emission reduction" and the realization of my country's "carbon neutrality" goals. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] In order to more clearly illustrate the technical solutions in the present invention or the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0077] Figure 1 It is a flow chart of urban green infrastructure biomass estimation by integrating multi-source remote sensing data according to an embodiment of the present invention.

[0078] Figure 2 It is a schematic diagram of distance and angle calculation of point cloud classification according to an embodiment of the present invention.

[0079] Figure 3 is a spatial distribution map of carbon storage of urban green infrastructure in an embodiment of the present invention, Figure 3 (a) shows the spatial distribution of carbon storage estimated using only LiDAR data. Figure 3 (b) shows the spatial distribution of carbon storage estimated using only multispectral images. Figure 3 (c) shows the spatial distribution of carbon stocks estimated by fusion of LiDAR data and multispectral images. DETAILED DESCRIPTION

[0080] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention will be clearly and completely described below in combination with the drawings and embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0081] Example 1

[0082] like Figure 1 As shown, an embodiment of the present invention provides a method for estimating carbon storage of urban green infrastructure by integrating multi-source data, including the following steps:

[0083] Step 1: Measure the tree height and tree species parameters of the sample points in the study area, and obtain the airborne lidar point cloud data and multispectral images of the study area.

[0084] Taking environmental factors and research purposes into consideration, this embodiment collected data from 110 representative and evenly distributed tree sample points to verify the accuracy of the urban green infrastructure carbon storage estimation method proposed in the present invention. The sampling points include six main tree species: camphor, cypress, magnolia, bamboo, metasequoia and ginkgo, among which camphor is the main tree species in the study area. For each sample point, the GPS location, tree height and tree species parameters are recorded. A drone equipped with a digital camera and a laser radar system was selected to scan the study area to obtain multispectral data and LiDAR data, and a total of 905 multispectral images of the study area were collected.

[0085] Step 2: preprocess the airborne lidar point cloud data and multispectral images obtained in step 1. The specific execution process includes the following steps:

[0086] Step 2.1, perform radiation correction, atmospheric correction and geometric correction on the multispectral image, and then perform stereo matching to calculate the three-dimensional coordinates to generate a dense point cloud. The specific calculation formula is as follows:

[0087] (1)

[0088] (2)

[0089] (3)

[0090] In the formula, is the camera focal length, is the stereo baseline, is parallax, To generate the 3D coordinates of points in a dense point cloud, is the coordinate of the principal point of the multispectral image, are the pixel coordinates of the multispectral image.

[0091] In step 2.2, the airborne lidar point cloud data in step 1 and the dense point cloud in step 2.1 are used as two point cloud data sources, and the adaptive hybrid model classification method is used to classify the two point cloud data sources respectively to obtain the ground points and canopy points of the two point cloud data sources.

[0092] Combining the advantages of irregular triangulated network (TIN) model in surface modeling and generalized triangular prism (GTP) model and tetrahedron model (TEN) in internal structure modeling, a TIN-GTP-TEN hybrid model is constructed to replace the traditional TIN model. The TIN-GTP-TEN hybrid model uses the tetrahedron model (TEN) as the basic modeling unit and refines the internal structure through the generalized triangular prism model (GTP). The specific construction method is as follows: First, the surface features of the study area are modeled through the existing improved Delaunay triangulation algorithm (Shewchuk, JR (1996). Triangle: Engineering a 2D Quality Mesh Generatorand Delaunay Triangulator. In Applied Computational Geometry TowardsGeometric Engineering, Springer, Berlin, Heidelberg, pp. 203-222.) to generate a TIN surface model, and then the GTP is used to divide the internal space below the TIN surface. Then, inside each GTP unit, the tetrahedron partitioning algorithm is used to decompose each heterogeneous GTP into smaller, homogeneous TENs to ensure the accurate expression of the surface features of the study area. Finally, the TIN-GTP-TEN hybrid model combines the surface TIN with the internal GTP and TEN to form a unified surface model of the features.

[0093] Combined with machine learning technology, a supervised learning-based regression model is used to automatically adjust the parameters of the TIN, GTP, and TEN models. The regression model is trained using existing high-precision point cloud data to learn the relationship between the surface complexity parameters of the object and the number of triangles in the TIN-GTP-TEN hybrid model. The trained model can predict the number of triangles that should be added to the corresponding area based on the input surface complexity parameters of the object, thereby generating more subdivided triangles in areas with higher surface complexity. The trained model does not require manual adjustment of the parameters of TIN, GTP, and TEN. Instead, through automated prediction and adjustment, the TIN-GTP-TEN hybrid model can more accurately represent surface details, thereby improving the accuracy of point cloud classification.

[0094] The airborne laser radar point cloud obtained in step 1 is input into the trained regression model. The regression model will automatically select points in the airborne laser radar point cloud data to generate the optimal TIN-GTP-TEN hybrid model, traverse each unclassified point in the airborne laser radar point cloud that is not involved in the construction of the optimal TIN-GTP-TEN hybrid model, and calculate the distance between each unclassified point and its projection point vertically projected onto the triangular surface in the optimal TIN-GTP-TEN hybrid model (such as Figure 2 midpoint P is the point to be classified, For point P The projection points vertically projected onto the triangular surface in the optimal TIN-GTP-TEN hybrid model are calculated. P arrive The vertical distance of the point to be classified is calculated by calculating the maximum angle between the vector from the point to be classified to its vertical projection to the vertex of the triangle surface in the optimal TIN-GTP-TEN hybrid model and the normal vector of the triangle surface (such as Figure 2 midpoint P is the point to be classified, the triangular surface ABC Yes P Project the surface vertically in the optimal TIN-GTP-TEN hybrid model and calculate the points P and the vector composed of each vertex of the triangle The angle between the normal vector of the triangle surface and the maximum angle is selected), and compared with the set distance threshold and angle threshold to determine whether the point to be classified is a ground point or a canopy point. The calculation formula of the classification algorithm is as follows:

[0095] (4)

[0096] In the formula, Indicate point P The classification results, Tree represents the canopy point, Ground represents a ground point, Yes P To its projection point vertically projected onto the triangle surface in the optimal TIN-GTP-TEN hybrid model The distance between Yes P To the maximum value of the angle between the vector perpendicularly projected to the vertex of the triangle surface in the optimal TIN-GTP-TEN hybrid model and the normal vector of the triangle surface, Indicate point P Vertically project triangles in the optimal TIN-GTP-TEN hybrid model, is the set distance threshold, which is 1.4m in this embodiment. is the set angle threshold, which is 8° in this embodiment;

[0097] The above calculation and judgment operations are performed on all points to be classified to complete the classification.

[0098] The dense point cloud generated by the multispectral image is classified using the above method to obtain canopy points and ground points.

[0099] Step 3: Generate digital elevation models (DEM) of two point cloud data sources based on the ground points obtained by point cloud classification in step 2.2; generate digital surface models (DSM) of two point cloud data sources based on canopy points obtained by point cloud classification in step 2.2; subtract DEM from the DSM of the two point cloud data sources to obtain the canopy height models (CHM) corresponding to the two point cloud data sources, and then obtain the CHM images corresponding to the two point cloud data sources.

[0100] Step 4: Fourier transform is used to fuse the CHM images corresponding to the two point cloud data sources in the same area to obtain a CHM fused image.

[0101] The images acquired by the multispectral sensor can finely reflect the canopy characteristics of ground vegetation, but cannot obtain information on the vertical structure of trees; the point cloud data obtained by the lidar can provide accurate three-dimensional geometric information, but the uneven density distribution will lead to information loss. The fusion of the color and texture data of the multispectral image and the elevation and three-dimensional structure data of the point cloud can effectively make up for the limitations of relying on one data source. This embodiment uses frequency domain fusion technology to fuse the CHM data generated by the airborne lidar point cloud data and the CHM data generated by the multispectral image, and independently processes and optimizes the different features of the two CHM images in the frequency domain. The high-frequency part mainly enhances the edge details of the image generated by the lidar, and the low-frequency part mainly maintains the large-area uniform structure of the multispectral image, thereby improving the quality of the fused image.

[0102] The specific process of image fusion using Fourier transform is as follows:

[0103] ① Perform Fourier transform on the CHM image generated by airborne lidar point cloud data and the CHM image generated by multispectral image to obtain its frequency domain image. The calculation formula is:

[0104] (6)

[0105] In the formula, represents a CHM image, is the spatial domain coordinate, represents the frequency domain image of the CHM image, is the frequency domain coordinate, represents the imaginary unit, ;

[0106] ② Use filters to separate the low-frequency and high-frequency parts of the frequency domain image;

[0107] Low pass filter The definition is as follows:

[0108] (7)

[0109] In the formula, is the horizontal cutoff frequency in the frequency domain, is the cutoff frequency in the vertical direction in the frequency domain, is the frequency domain coordinate;

[0110] High Pass Filter The definition is as follows:

[0111] (8)

[0112] The low-pass filter and high pass filter With frequency domain image Multiply and separate the desired frequency components, namely:

[0113] (9)

[0114] (10)

[0115] In the formula, represents the frequency domain image of the CHM image, represents the low-frequency image of the CHM image, Represents the high-frequency image of the CHM image;

[0116] ③ The high-frequency and low-frequency parts of two CHM images in the same area of ​​the study area are fused respectively;

[0117] For low-frequency components, the weighted average method is used for fusion. The specific calculation formula is:

[0118] (11)

[0119] In the formula, represents the fused low-frequency image, A low-frequency image representing the CHM image generated from airborne lidar point cloud data, Represents the low-frequency image of the CHM image generated by multispectral imaging, is a weight coefficient between 0 and 1;

[0120] For high-frequency components, the maximum value method is used for fusion. The specific calculation formula is:

[0121] (12)

[0122] In the formula, represents the fused high-frequency image, Indicates taking the maximum value, A high-frequency image representing the CHM image generated from LiDAR data, High-frequency image representing the CHM image generated by multispectral imaging;

[0123] ④ Inverse Fourier transform, merge the fused low-frequency and high-frequency components, and transform them back to the spatial domain through inverse Fourier transform to obtain the CHM fusion image , the specific calculation method is as follows:

[0124] (13)

[0125] In the formula, is the spatial domain coordinate, represents the fused low-frequency image, represents the fused high-frequency image, is the frequency domain coordinate, represents the imaginary unit, .

[0126] Step 5: Based on the CHM fusion image, the tree height parameters are extracted by combining edge detection and watershed segmentation algorithm.

[0127] Step 5.1, calculate the gradient magnitude and gradient direction of the CHM fusion image.

[0128] Gray value matrix of CHM fusion image , calculate the gradient magnitude and the gradient direction , the specific calculation formula is as follows:

[0129] (14)

[0130] (15)

[0131] In the formula, and They are the gray value matrices exist x Direction and y Partial derivatives in the direction.

[0132] Step 5.2, set the initial water source point and assign a value to the water level depth of the initial water source point.

[0133] The Canny edge detection algorithm is used to extract the edges in the CHM fusion image and preliminarily identify the outline of the tree crown. The points in the edge image with pixel grayscale values ​​lower than the surrounding 8 neighboring pixels are taken as the initial water source points, and the water level depth of these initial water source points is calculated. Set to the grayscale value of this point.

[0134] Step 5.3, update the water level depth of the pixel on the CHM fusion image.

[0135] By gradient direction Iterative simulation updates the convergence of water flow: ① Starting from any water source point, find the gradient direction The next pixel , if the next pixel Water level depth Greater than , then update the next pixel Water level depth Otherwise, do not update ; ② Then from the pixel point Start by looking for the gradient direction The next pixel , if the next pixel Water level depth Greater than , then update the next pixel Water level depth Otherwise, do not update ; ③ Repeat ①-② until the water level depth of all pixels is no longer updated.

[0136] Step 5.4, preliminarily perform watershed segmentation on the CHM fusion image.

[0137] When the pixel Water level depth Equal to its gray value in DEM When , the pixel point is considered to be the boundary point of different watersheds. Starting from any boundary point, the next boundary point is tracked along the direction of gradient change until the tracked boundary points are connected to form a closed area or the image edge has been tracked. The tracking is ended and all the tracked boundary points are connected in sequence to form the boundary line of the watershed. The above operation is repeated for each boundary point to form multiple segmentation areas, and the watershed segmentation of the CHM fusion image is realized.

[0138] Step 5.5: Adjust the water source point based on the watershed segmentation result in step 5.4.

[0139] The segmented area in step 5.4 is used as the statistical range for judging the average grayscale value or grayscale variance, and the water source point is adjusted based on this range. If the average grayscale value of a certain area after segmentation is significantly lower than the average of other surrounding areas, the pixels in this area are added as new water source points. If the grayscale variance of a certain area after segmentation is significantly higher than other surrounding areas, this indicates that the pixel value distribution in this area is extremely uneven, so the pixels in this area are added as new water source points.

[0140] Step 5.6, perform the operations of step 5.3-step 5.4 on the adjusted water source point to obtain a new watershed segmentation result.

[0141] Perform the operations of steps 5.3-5.4 on the initial water source points set in step 5.2 and the water source points added in step 5.5 to obtain a new watershed segmentation result.

[0142] Step 5.7, find the local maximum grayscale value in each area obtained by the watershed segmentation in step 5.6 and use it as the tree height.

[0143] In this example, a total of 1627 tree height parameters were extracted in the study area.

[0144] Step 6: Substitute the tree height parameters of the study area extracted in step 5 into the allometric growth equations of the corresponding tree species to obtain the biomass, and then multiply it by the corresponding carbon content of each tree species to obtain the carbon storage of a single tree. Count the carbon storage of all single trees to obtain the preliminary carbon storage of the urban green infrastructure in the study area.

[0145] Vegetation carbon storage refers to the total amount of carbon absorbed and stored by vegetation through photosynthesis within a certain time and space range. It is one of the important indicators for measuring the carbon sequestration capacity, health status and biodiversity of ecosystems. In the complex urban environment, the landscape attributes of green infrastructure show obvious heterogeneity. Different tree species play different ecological functions, carbon sequestration capacity, landscape value, economic value and diversity in urban green infrastructure, which are of great importance and cannot be ignored. Urban green infrastructure is closely related to human activities. It is crucial to consider the impact of human activities when estimating carbon storage.

[0146] This embodiment takes into account the different characteristics of each dominant tree species and the impact of human activities in the carbon storage estimation, and evaluates the carbon storage of individual tree species in green infrastructure based on the allometric growth equations and carbon content of each dominant tree species in published literature. The tree height parameter H of the study area extracted in step 5 is input into the allometric growth equation corresponding to each dominant tree species, and the biomass AGB is calculated, and then multiplied by the corresponding carbon content to obtain the carbon storage of a single tree. The carbon storage of all single trees is counted to obtain the preliminary carbon storage of the urban green infrastructure in the study area. The allometric growth equations and carbon content of each dominant tree species are as follows:

[0147] Table 1 Allometric growth equations and carbon content of dominant tree species in the study area

[0148]

[0149] Step 7: Calculate the human activity correction coefficient of the study area in combination with the quantitative indicators, correct the carbon storage obtained in step 6, and obtain the final carbon storage of urban green infrastructure in the study area.

[0150] Taking into account the impact of human activities, a correction coefficient calculation formula that includes human activity factors is designed. The influencing factors include human activity factors (construction, destruction, fertilization), urban space factors (urbanization rate), ecological sensitivity factors, and social factors (population density). According to the intensity and frequency of each influencing factor, a comprehensive impact value is calculated as the human activity correction coefficient of the study area. The preliminary carbon storage estimated in step 6 is corrected to generate a more accurate urban green infrastructure carbon storage estimation result.

[0151] Human activity correction factor The calculation formula is as follows:

[0152] (16)

[0153] In the formula, is the human activity factor, obtained by multiplying activity intensity and activity frequency; Indicates different impact factor types, Indicates trimming, Indicates fertilization. It means destruction; S is the urban space factor, obtained from the urbanization rate of the study area; E It is the ecological sensitivity factor, and its value range is usually between 0 and 1, where 0 means insensitive and 1 means highly sensitive. It is divided into four sensitivity levels: highly sensitive ( )、Generally sensitive( )、Low Sensitivity( ) and insensitive ( ); A It is a social factor, obtained by the ratio of the population of the study area to the area of ​​the study area.

[0154] In this embodiment, through investigation and data analysis, the human activity correction factor is calculated as follows: The calculation of the response activity combines the intensity and frequency. For example, the pruning intensity is set to 0.2 and the frequency is 2 times / year. is 0.2×2=0.4; the fertilization intensity is set to 0.1, the frequency is once a year, and the fertilization is 0.1×1=0.1; the intensity of the damage is set to 0.5, the frequency is 0.1 times / year, and the damage = 0.5×0.1=0.05. According to the statistics bureau data, the urbanization rate of the study area is 86.06%. S Set to 0.8606. Through ecological sensitivity analysis, the ecological sensitivity of the study area is generally sensitive ( ), in this embodiment, E According to the demographic yearbook, the ratio of the population of the study area to the area of ​​the entire study area is 37.6%, so A The human activity factor of the study area in this embodiment is obtained through the above calculations: , urbanization rate ( ), ecological sensitivity ( ) and population density ratio ( ). For the first human activity (pruning) ( ), , for the second human activity (fertilization) ( ), , for the third human activity (destruction) ( ), ,but Therefore, the human activity correction factor of the study area is 0.8932. Multiply the preliminary carbon storage obtained in step 6 by the human activity correction factor 0.8932 to obtain the final carbon storage of the study area.

[0155] In order to verify the effectiveness of the proposed method, based on the tree height and tree species parameters of the measured sample data, the carbon storage data of the measured sample points were calculated using the allometric growth equation and carbon content in Table 1. The tree height and carbon storage data of the measured sample points were fitted and verified with the tree height and carbon storage data estimated and extracted from three data sources: LiDAR, multispectral remote sensing images, and fusion data. The determination coefficient (R 2 The tree height parameter extraction and carbon storage estimation results of the present invention using three data sources were evaluated using three indicators: , root mean square error (RMSE) and bias. The evaluation results are shown in Table 2.

[0156] Table 2 Comparison of indicators of different data sources in tree height extraction and carbon storage estimation

[0157]

[0158] As can be seen from Table 2, there are great differences in the estimation of carbon storage of urban green infrastructure when selecting different data sources. The data fusion method proposed in this paper is effective in tree height extraction (RMSE=1.76 m, R 2=0.75) and carbon storage estimation (RMSE=28.88 kg C, Bias=-2.41, R 2 =0.86) achieved the best results.

[0159] In order to facilitate the visual comparative analysis of the heterogeneity of carbon storage in small-area urban green facility ecosystems due to different data sources and the spatial distribution of dominant tree species in green infrastructure, the carbon storage estimation results of a single data source (LiDAR or multispectral image) and the fusion of LiDAR data and multispectral image based on the method proposed in this invention are interpolated by Kriging to obtain the spatial distribution of carbon storage as shown in the figure. Figure 3 , unit is kg C. Figure 3 (a) shows the spatial distribution of carbon storage estimated using only LiDAR data (LiDAR data are filtered and classified, and the classified canopy points are substituted into the allometric equation to calculate the biomass, which is then multiplied by the corresponding carbon content to obtain the carbon storage, and corrected using the human activity correction factor). LiDAR-based estimates provide the highest spatial resolution of carbon storage, effectively capturing spatial variations, especially those in densely vegetated areas, and highlighting the distribution differences of different vegetation densities. Figure 3 (b) shows the spatial distribution of carbon storage estimated using only multispectral images (stereo matching of multispectral images is performed to generate a dense point cloud, and the classified canopy points are substituted into the allometric equation to calculate the biomass, which is then multiplied by the corresponding carbon content to obtain the carbon storage, and corrected using the human activity correction factor). Compared with lidar, the carbon storage estimates of multi-view spectral images in densely vegetated areas are much lower, which may lead to an underestimation of carbon storage in densely vegetated areas. Figure 3 (c) shows the spatial distribution of carbon storage estimated by fusing LiDAR data and multispectral shadows based on the method proposed in the present invention. The fusion method makes full use of the advantages of the two data sources and can more comprehensively and accurately reflect the spatial distribution of carbon storage, especially in urban areas where vegetation types and densities are unevenly distributed.

[0160] Example 2

[0161] Based on the same inventive concept, the present invention also provides a system for estimating the carbon stock of urban green infrastructure by integrating multi-source data, including a processor and a memory, the memory being used to store program instructions, and the processor being used to call the program instructions in the memory to execute a method for estimating the carbon stock of urban green infrastructure by integrating multi-source data as described above.

[0162] Example 3

[0163] Based on the same inventive concept, the present invention also provides a system for estimating the carbon stock of urban green infrastructure by integrating multi-source data, including a readable storage medium on which a computer program is stored. When the computer program is executed, a method for estimating the carbon stock of urban green infrastructure by integrating multi-source data as described above is implemented.

[0164] In specific implementation, the method proposed in the technical solution of the present invention can be implemented by technical personnel in this field using computer software technology to realize automatic operation process. System devices for implementing the method, such as computer-readable storage media storing the corresponding computer program of the technical solution of the present invention and computer equipment running the corresponding computer program, should also be within the scope of protection of the present invention.

[0165] The specific embodiments described herein are merely examples of the spirit of the present invention. Those skilled in the art may make various modifications or additions to the specific implementation cases described, or replace them with similar methods, without departing from the spirit of the present invention or exceeding the scope defined by the appended claims.

Claims

1. A method for estimating carbon storage of urban green infrastructure by integrating multi-source data, characterized in that: The following steps are involved: Step 1: measure the tree height and tree species parameters of the sample points in the study area, and obtain the airborne lidar point cloud data and multispectral images of the study area; Step 2, preprocessing the airborne laser radar point cloud data and multispectral image obtained in step 1; Step 2.1, perform radiation correction, atmospheric correction and geometric correction on the multispectral image, and then perform stereo matching to calculate the three-dimensional coordinates to generate a dense point cloud; Step 2.2, using the airborne laser radar point cloud data in step 1 and the dense point cloud in step 2.1 as two point cloud data sources, using the adaptive hybrid model classification method to classify the two point cloud data sources respectively, and obtain the ground points and canopy points of the two point cloud data sources; Step 3, based on the ground points obtained by point cloud classification in step 2.2, digital elevation models DEM of two point cloud data sources are generated, based on the canopy points obtained by point cloud classification in step 2.2, digital surface models DSM of two point cloud data sources are generated, and the DSM of the two point cloud data sources are respectively subtracted from the DEM to obtain the canopy height models CHM corresponding to the two point cloud data sources, and then the CHM images corresponding to the two point cloud data sources are obtained; Step 4, Fourier transform is used to fuse the CHM images corresponding to the two point cloud data sources in the same area of ​​the study area to obtain a CHM fused image; Step 5, based on the CHM fusion image, the tree height parameters are extracted by combining edge detection and watershed segmentation algorithm; Step 5.1, calculate the gradient magnitude and gradient direction of the CHM fusion image; Step 5.2, set the initial water source point and assign a value to the water level depth of the initial water source point; The Canny edge detection algorithm is used to extract the edges in the CHM fusion image, and the outline of the tree crown is preliminarily identified. The points in the edge image with pixel grayscale values ​​lower than the surrounding 8 neighboring pixels are taken as the initial water source points, and the water level depth of these initial water source points is calculated. Set to the grayscale value of the point; Step 5.3, updating the water level depth of the pixel on the CHM fusion image; Step 5.4, preliminarily perform watershed segmentation on the CHM fusion image; Step 5.5, based on the watershed segmentation result of step 5.4, add a new water source point; Step 5.6, perform the operations of steps 5.3-5.4 on the initial water source points set in step 5.2 and the water source points added in step 5.5 to obtain a new watershed segmentation result; Step 5.7, find the local maximum gray value in each area obtained by watershed segmentation in step 5.6, and use it as the tree height; Step 6: Substitute the tree height parameters of the study area extracted in step 5 into the allometric growth equation of the corresponding tree species to obtain the biomass, and then multiply it by the carbon content of each tree species to obtain the carbon storage of a single tree. Count the carbon storage of all single trees to obtain the preliminary carbon storage of the urban green infrastructure in the study area. Step 7: Calculate the human activity correction coefficient of the study area in combination with the quantitative indicators, and correct the preliminary carbon storage obtained in step 6 to obtain the final carbon storage of urban green infrastructure in the study area.

2. The method for estimating carbon storage of urban green infrastructure by integrating multi-source data according to claim 1, characterized in that: The specific operations of the classification process in step 2.2 are as follows: ① Use the existing high-precision point cloud data to train the regression model so that it can learn the relationship between the surface complexity parameters of the object and the number of triangles in the TIN-GTP-TEN hybrid model; ② Input the airborne laser radar point cloud data into the trained regression model. The regression model will automatically select points in the airborne laser radar point cloud data to generate the optimal TIN-GTP-TEN hybrid model, traverse the points to be classified in the airborne laser radar point cloud data that are not involved in the construction of the optimal TIN-GTP-TEN hybrid model, calculate the distance between each point to be classified and its projection point vertically projected onto the triangular surface in the optimal TIN-GTP-TEN hybrid model, calculate the maximum value of the angle between the vector from the point to be classified to its vertex on the triangular surface in the optimal TIN-GTP-TEN hybrid model and the normal vector of the triangular surface, and compare it with the set distance threshold and angle threshold to determine whether the point to be classified is a ground point or a canopy point. The calculation formula of the classification algorithm is as follows: （4） In the formula, Indicate point P The classification results, Tree represents the canopy point, Ground represents a ground point, Yes P To its projection point vertically projected onto the triangle surface in the optimal TIN-GTP-TEN hybrid model The distance between Yes P To the maximum value of the angle between the vector perpendicularly projected to the vertex of the triangle surface in the optimal TIN-GTP-TEN hybrid model and the normal vector of the triangle surface, Indicate point P Vertically project triangles in the optimal TIN-GTP-TEN hybrid model, is the set distance threshold, is the set angle threshold; Perform the above calculation and judgment operations on all points to be classified to complete the classification; ③Use the processing method in ② to classify the dense point cloud generated by multispectral images.

3. The urban green infrastructure carbon storage estimation method integrating multi-source data according to claim 2 is characterized by: The TIN-GTP-TEN hybrid model in step 2.2 is constructed as follows: first, the surface features of the study area are modeled using the existing improved Delaunay triangulation algorithm to generate a TIN surface model, and then the internal space below the TIN surface is divided using GTP. Then, within each GTP unit, the tetrahedron partitioning algorithm is used to decompose each heterogeneous GTP into smaller, homogeneous TENs to ensure accurate expression of the surface features of the study area.

4. The urban green infrastructure carbon storage estimation method integrating multi-source data according to claim 1, characterized in that: The specific process of image fusion using Fourier transform in step 4 is as follows: ① Perform Fourier transform on the CHM image generated by airborne lidar point cloud data and the CHM image generated by multispectral image to obtain its frequency domain image. The calculation formula is: （6） In the formula, represents a CHM image, is the spatial domain coordinate, represents the frequency domain image of the CHM image, is the frequency domain coordinate, represents the imaginary unit, ; ② Use filters to separate the low-frequency and high-frequency parts of the frequency domain image; Low pass filter The definition is as follows: （7） In the formula, is the horizontal cutoff frequency in the frequency domain, is the cutoff frequency in the vertical direction in the frequency domain, is the frequency domain coordinate; High Pass Filter The definition is as follows: （8） The low-pass filter and high pass filter With frequency domain image Multiply and separate the desired frequency components, namely: （9） （10） In the formula, represents the frequency domain image of the CHM image, represents the low-frequency image of the CHM image, High frequency image representing the CHM image; ③ The high-frequency and low-frequency parts of two CHM images in the same area of ​​the study area are fused respectively; For low-frequency components, the weighted average method is used for fusion. The specific calculation formula is: （11） In the formula, represents the fused low-frequency image, A low-frequency image representing the CHM image generated from airborne lidar point cloud data, Represents the low-frequency image of the CHM image generated by multispectral imaging, is a weight coefficient between 0 and 1; For high-frequency components, the maximum value method is used for fusion. The specific calculation formula is: （12） In the formula, represents the fused high-frequency image, Indicates taking the maximum value, A high-frequency image representing the CHM image generated from LiDAR data, High-frequency image representing the CHM image generated by multispectral imaging; ④ Inverse Fourier transform, merge the fused low-frequency and high-frequency components, and transform them back to the spatial domain through inverse Fourier transform to obtain the CHM fusion image , the specific calculation method is as follows: （13） In the formula, is the spatial domain coordinate, represents the fused low-frequency image, represents the fused high-frequency image, is the frequency domain coordinate, represents the imaginary unit, .

5. The method for estimating carbon storage of urban green infrastructure by integrating multi-source data according to claim 1, characterized in that: Gray value matrix of CHM-based fusion image in step 5.1 Calculate the gradient magnitude and the gradient direction , the specific calculation formula is as follows: （14） （15） In the formula, and They are the gray value matrices exist x Direction and y Partial derivatives in the direction.

6. The method for estimating carbon storage of urban green infrastructure by integrating multi-source data according to claim 5, characterized in that: In step 5.3, the gradient direction Iterative simulation updates the convergence of water flow: ① Starting from any water source point, find the gradient direction The next pixel , if the next pixel Water level depth Greater than , then update the next pixel Water level depth Otherwise, do not update ; ②Then from the pixel point Start by looking for the gradient direction The next pixel , if the next pixel Water level depth Greater than , then update the next pixel Water level depth Otherwise, do not update ; ③ Repeat ①-② until the water level depth of all pixels is no longer updated.

7. The urban green infrastructure carbon storage estimation method integrating multi-source data according to claim 1, characterized in that: In step 5.4, when the pixel Water level depth Equal to its gray value in DEM When , the pixel point is considered to be the boundary point of different watersheds. Starting from any boundary point, the next boundary point is tracked along the direction of gradient change until the tracked boundary points are connected to form a closed area or the image edge has been tracked. The tracking is ended and all the tracked boundary points are connected in sequence to form the boundary line of the watershed. The above operation is repeated for each boundary point to form multiple segmentation areas, thus realizing the watershed segmentation of the CHM fusion image.

8. The urban green infrastructure carbon storage estimation method integrating multi-source data according to claim 1 is characterized in that: In step 5.5, the area segmented in step 5.4 is used as the statistical range for judging the average gray value or gray variance, and the water source point is adjusted based on this range; if the average gray value of a certain area after segmentation is significantly lower than the average of other surrounding areas, the pixels of this area are added as new water source points; if the gray variance of a certain area after segmentation is significantly higher than that of other surrounding areas, the pixels of this area are added as new water source points.

9. The urban green infrastructure carbon storage estimation method integrating multi-source data according to claim 1, characterized in that: Human activity correction factor in step 7 The calculation formula is as follows: （16） In the formula, is the human activity factor, obtained by multiplying activity intensity and activity frequency; Indicates different impact factor types, Indicates trimming, Indicates fertilization. It means destruction; S is the urban space factor, obtained from the urbanization rate of the study area; E It is an ecological sensitivity factor obtained through ecological sensitivity analysis and is divided into four sensitivity levels: highly sensitive , generally sensitive , low sensitivity and insensitivity ; A It is a social factor, obtained by the ratio of the population of the study area to the area of ​​the study area; Multiply the preliminary carbon storage obtained in step 6 by the human activity correction coefficient to obtain the final carbon storage of the study area.

10. A carbon storage estimation system for urban green infrastructure integrating multi-source data, characterized in that: It includes a processor and a memory, the memory is used to store program instructions, and the processor is used to call the program instructions in the memory to execute a method for estimating carbon reserves of urban green infrastructure by integrating multi-source data as described in any one of claims 1 to 9.

Images