A method and device for extracting a single-tree scale young forest by fusing two-dimensional and three-dimensional features, and electronic equipment

Abstract

The application discloses a low-closedness single-tree scale young forest remote sensing extraction scheme fusing two-dimensional and three-dimensional features, generates multi-view core line images based on two-line array stereo images and RPC parameters of GF-7; generates pyramid images for the multi-view core line images and determines a disparity map of the two-line array stereo images based on the pyramid images; calculates point clouds based on the disparity map and performs grid processing on the point clouds to generate a DSM; optimizes the DSM by using laser height measurement data of a satellite, obtains an optimized DSM, and performs image analysis on the optimized DSM by using an isolated forest algorithm to extract pixels containing young forests; fuses a panchromatic backward-looking image and a multi-spectral image of satellite remote sensing images to generate a panchromatic-multi-spectral fusion image; performs shadow extraction on the panchromatic-multi-spectral fusion image by using a shadow detection algorithm to obtain a target shadow area; compares the pixels containing the young forests with the target shadow area to determine a target young forest extraction result, and the accuracy of the young forest extraction result can be improved.

Description

Technical Field

[0001] This invention relates to the field of image processing technology, and in particular to a method, apparatus, and electronic device for remote sensing extraction of single-tree-scale young forests by fusing two- and three-dimensional features. Background Technology

[0002] Afforestation is an important means of ecological restoration. Monitoring afforestation can accurately assess the ecological quality of forest ecosystems and provide a valid basis for taking corresponding measures. Traditional methods often use field sampling surveys of plots to monitor trees and forests. Although plot surveys are highly reliable, the survey process is time-consuming and labor-intensive. In areas with steep terrain and complex geographical environments, it is even more difficult for surveyors to enter.

[0003] With the maturity of satellite technology, using satellite imagery to extract young forests at the individual tree scale has become a feasible technical means. However, compared to the size of young forests, the resolution of available satellite images is relatively low. When the target size is small, a large amount of background information such as bare soil, grass, and shrubs surrounding the young forest can interfere with its identification, easily leading to false positives and false negatives. Furthermore, compared to mature trees, the crown size of young forests is even smaller, making them small-scale and blurry targets even in high spatial resolution imagery. The information provided by the target is very limited, making extraction difficult.

[0004] Therefore, there is an urgent need for those skilled in the art to provide a solution that can effectively process existing high spatial resolution (about 1 meter) satellite images and eliminate interference from various types of background information in order to accurately identify young forests at the single-tree scale. Summary of the Invention

[0005] The purpose of this invention is to provide a remote sensing extraction method, device, and electronic device for low-canopy-density single-tree-scale young forests that integrates two-dimensional and three-dimensional features, which can solve the problem in the prior art that satellite images cannot accurately identify young forests at the single-tree scale.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution:

[0007] This invention provides a remote sensing method for extracting young forests at the single-tree scale with low canopy closure by fusing two-dimensional and three-dimensional features, wherein the method includes:

[0008] Acquire satellite remote sensing images, and generate multi-view collimator images based on the two-line array stereo images and RPC parameters contained in the satellite remote sensing images;

[0009] A pyramid image is generated for the multi-view epipolar image, and a disparity map of the two linear array stereo images is determined based on the pyramid image;

[0010] The point cloud is calculated based on the disparity map, and the calculated point cloud is meshed to generate a digital surface model (DSM).

[0011] The DSM is optimized using satellite laser altimetry data to obtain the optimized DSM;

[0012] The isolated forest algorithm was used to perform image analysis on the optimized DSM and extract pixels containing young forests.

[0013] The panchromatic backview image and the multispectral image of the satellite remote sensing image are fused to generate a panchromatic-multispectral fused image;

[0014] A shadow detection algorithm is used to extract shadows from the panchromatic-multispectral fused image to obtain the target shadow region;

[0015] The pixels containing young forests are compared with the target shadow area to determine the target young forest extraction result; if there is no shadow area near the pixel, the pixel extraction is determined to be incorrect.

[0016] Optionally, the step of acquiring satellite remote sensing imagery and generating multi-view epipolar imagery based on the two-line array stereo imagery and RPC parameters contained in the satellite remote sensing imagery includes:

[0017] Acquire satellite remote sensing images;

[0018] The satellite back-view image contained in the satellite remote sensing image is used as the main image, and the satellite front-view image is used as the auxiliary image, and projection connection points between the main image and the auxiliary image are generated.

[0019] Based on the projection connection points between the main image and the auxiliary image, a two-way coordinate transformation model between the quasi-pixel image and the original image is established, wherein the original image refers to the main image and the auxiliary image, and each original image corresponds to one quasi-pixel image in the multi-view quasi-pixel image;

[0020] Based on the quasi-epidial line image and the bidirectional coordinate transformation model, the image plane coordinates of the quasi-epidial line image are reversed to obtain the image plane coordinates of the original image.

[0021] Based on the image plane coordinates, RPC parameters, and object elevation of the original image, calculate the RPC parameters of the multi-view epipolar image.

[0022] Based on the RPC parameters of the multi-view epipolar image and the coordinate transformation model, the coordinates of the original image points corresponding to the image points of the multi-view epipolar image are calculated point by point.

[0023] The gray levels of the original image points are obtained by bilinear interpolation and used as the gray levels of the multi-view epipolar image to generate the multi-view epipolar image.

[0024] Optionally, the step of using the satellite back-view image included in the satellite remote sensing imagery as the main image and the satellite forward-view imagery as the auxiliary imagery, and generating projection connection points between the main imagery and the auxiliary imagery, includes:

[0025] Use satellite rear-view images as the main images and satellite forward-view images as auxiliary images;

[0026] Using the auxiliary image as the target auxiliary image, the import order of the main image and the auxiliary image is determined, and the average epipolar tilt angle is calculated; wherein, the main image is imported first as a reference, and the auxiliary image is imported later and matched with the main image;

[0027] The range of the multi-view epipolar image is determined based on the average epipolar tilt angle.

[0028] Generate projection connection points between the main image and the auxiliary image that cover the range of the main image quasi-piezoresistive image.

[0029] Optionally, the step of establishing a coordinate transformation model between the quasi-piston image and the original image based on the projection connection points between the main image and the auxiliary image includes:

[0030] Calculate the offset ratio coefficients of the auxiliary image along the direction of the quasi-core image and the auxiliary image perpendicular to the direction of the quasi-core image;

[0031] Calculate the projection trajectory tilt parameter, wherein the projection trajectory tilt parameter represents the angle between the projection trajectory line and the original image line direction;

[0032] Based on the offset ratio coefficient and the projection trajectory tilt parameter, a coordinate transformation model is established between the quasi-core line image and the corresponding original image.

[0033] Optionally, the step of generating a pyramid image for the multi-view epipolar image and determining the disparity map of the two linear array stereo images based on the pyramid image includes:

[0034] Match the top pyramid image of the quasi-nuclear line image of the satellite two-line array stereo image to obtain the disparity map of the top pyramid image;

[0035] The disparity map of the other pyramid images, excluding the top and bottom layers of the multi-view collimator image of the satellite two-line array stereo image, is obtained by matching the pyramid images of the other pyramid images.

[0036] The bottom pyramid image of the quasi-nuclear line image of the satellite two-line array stereo image is matched to obtain the disparity map of the top pyramid image.

[0037] Optionally, the step of calculating the point cloud based on the disparity map and generating a digital surface model (DSM) by meshing the calculated point cloud includes:

[0038] Calculate cloud points based on the disparity map;

[0039] Calculate the planar coordinate range of all point clouds, and based on the planar coordinate range, calculate the length and width of the DSM to be generated; wherein the planar coordinate range includes minimum longitude, maximum longitude, minimum latitude, and maximum latitude;

[0040] Based on the length and width of the DSM to be generated, and the three-dimensional coordinates of each point cloud, calculate the row and column coordinates of the DSM grid points corresponding to each point cloud, and calculate the latitude, longitude coordinates and elevation of the DSM grid points according to the member points of each DSM grid point.

[0041] The digital surface model (DSM) is generated by organizing the latitude, longitude, and elevation coordinates of each DSM grid point into a regular grid.

[0042] Optionally, the step of optimizing the DSM using satellite laser altimetry data to obtain the optimized DSM includes:

[0043] By combining the laser points in the satellite laser altimetry data selected from the field detection points, the elevation values ​​of the DSM grid points corresponding to the geographical locations of the control points and the elevation differences between the control points are calculated.

[0044] The least squares method was used to fit the functional relationship between the elevation difference and the latitude and longitude of the DSM grid points corresponding to the control points:

[0045] Substitute the latitude and longitude coordinates of each grid point on the DSM into the function relationship to obtain the corrected elevation difference, and obtain the corrected digital surface model.

[0046] Arrange the elevation differences in order of magnitude and take the median value as the correction value;

[0047] The elevation value of each grid point in the corrected digital surface model is added to the correction value to obtain the optimized DSM.

[0048] Optionally, the isolated forest algorithm is used to perform image analysis on the optimized DSM, and the step of extracting pixels containing young forests includes:

[0049] The optimized DSM is regarded as a raw dataset, wherein the raw dataset includes multiple cells, and each cell corresponds to a DSM value;

[0050] m data points are randomly selected without replacement from the original dataset as training data, and the training data are placed in the root node of a binary tree;

[0051] If the number of cells in the current node is 1, or the preset maximum tree height has been reached, then the node is marked as a leaf node and the segmentation stops; otherwise, a segmentation value is randomly selected, wherein the segmentation value is generated between the maximum and minimum values ​​of the DSM values ​​of all cells in the current node.

[0052] Based on the segmentation value, the cells of the current node are divided into two subsets: cells with a DSM value less than the segmentation value enter the left subtree, and the rest enter the right subtree.

[0053] Recursively perform the above segmentation process on the left and right subtree nodes to construct an isolated tree;

[0054] Based on the isolated tree, calculate the path length from the root node to the leaf node for each cell in the original dataset, and calculate the anomaly score for each cell based on the path length.

[0055] Based on the abnormal scores, pixels with abnormal scores higher than a preset threshold, or pixels with abnormal scores ranking in the top preset percentage, are identified as pixels containing young forests.

[0056] Optionally, the step of extracting shadows from the panchromatic-multispectral fused image using a shadow detection algorithm to obtain the target shadow region includes:

[0057] The blue-green band difference image is obtained by subtracting the blue band of the panchromatic-multispectral fused image from the green band of the panchromatic-multispectral fused image.

[0058] The shadow region of the near-infrared band image of the panchromatic-multispectral fused image is extracted using the histogram thresholding method to obtain the first shadow region extraction result;

[0059] The shadow region of the blue-green band difference image was extracted using the histogram thresholding method, and the extraction results were then subjected to median filtering and boundary tracing to obtain the boundary line.

[0060] A buffer of a specific width is constructed based on the boundary line to obtain the extraction result of the second shaded area;

[0061] The first shadow region extraction result and the second shadow region extraction result are overlaid to obtain the shadow boundary region;

[0062] The shadow boundary region is filled to obtain the target shadow region.

[0063] This invention provides a remote sensing extraction device for single-tree-scale young forests that integrates two-dimensional and three-dimensional features, wherein the device includes:

[0064] The acquisition module is used to acquire satellite remote sensing images and generate multi-view epipolar images based on the two-line array stereo images and RPC parameters contained in the satellite remote sensing images.

[0065] The disparity map determination module is used to generate a pyramid image for the multi-view epipolar image and determine the disparity map of the two linear array stereo images based on the pyramid image.

[0066] The model generation module is used to calculate point clouds based on the disparity map, and to generate a digital surface model (DSM) by meshing the calculated point clouds.

[0067] An optimization module is used to optimize the DSM using satellite laser altimetry data to obtain an optimized DSM.

[0068] The extraction module is used to perform image analysis on the optimized DSM using the isolated forest algorithm to extract pixels containing young forests;

[0069] The fused image generation module is used to fuse the panchromatic backview image and the multispectral image of the satellite remote sensing image to generate a panchromatic-multispectral fused image;

[0070] The shadow region determination module is used to extract shadows from the panchromatic-multispectral fused image using a shadow detection algorithm to obtain the target shadow region.

[0071] The result determination module is used to compare the pixel containing the young forest with the target shadow area to determine the target young forest extraction result; wherein, if there is no shadow area near the pixel, the pixel extraction is determined to be incorrect.

[0072] This invention provides an electronic device, which includes a processor, a memory, and a program or instructions stored in the memory and executable on the processor. When the program or instructions are executed by the processor, they implement the steps of any of the above-described remote sensing extraction methods for low-canopy-density single-tree-scale young forests that integrate two-dimensional and three-dimensional features.

[0073] This invention provides a readable storage medium storing a program or instructions. When the program or instructions are executed by a processor, they implement the steps of any of the above-described remote sensing extraction methods for low-canopy-density single-tree-scale young forests that integrate two-dimensional and three-dimensional features.

[0074] This application provides a remote sensing extraction scheme for low-canopy-density single-tree-scale young forests that integrates two-dimensional and three-dimensional features. The scheme acquires satellite remote sensing images, generates multi-parallel epipolar images based on the two-line array stereo images and RPC parameters contained in the satellite remote sensing images, generates pyramid images from the multi-parallel epipolar images, and determines the disparity map of the two-line array stereo images based on the pyramid images. Point clouds are calculated based on the disparity map, and the calculated point clouds are meshed to generate a digital surface model (DSM). The DSM is optimized using satellite laser altimetry data to obtain an optimized DSM. An isolated forest algorithm is used to analyze the optimized DSM and extract pixels containing young forests. The panchromatic backview image and multispectral image of the satellite remote sensing images are fused to generate a panchromatic-multispectral fused image. A shadow detection algorithm is used to extract shadows from the panchromatic-multispectral fused image to obtain target shadow regions. Pixels containing young forests are compared with the target shadow regions to determine the target young forest extraction result. The solution provided in this application has two aspects. First, by using DSM (Dense Matching of Stereo Images) optimized from satellite laser altimetry data to extract tree height information, it can eliminate the interference of low-lying herbaceous plants and shrubs on the extraction of young forests. Second, by fusing tree height with the shadow information generated by trees under low canopy closure, the problem of two-dimensional small target detection is transformed into two-dimensional tree shadow detection and three-dimensional tree height feature extraction, which are relatively easier to solve, and can further improve the accuracy of young forest extraction at the single-tree scale. Attached Figure Description

[0075] Figure 1 This is a flowchart illustrating the steps of a remote sensing extraction method for low-canopy-density single-tree-scale young forests that integrates two-dimensional and three-dimensional features, according to an embodiment of this application.

[0076] Figure 2 This is a diagram illustrating the implementation steps of another remote sensing method for extracting young forests at the single-tree scale with low canopy closure by integrating two- and three-dimensional features, according to an embodiment of this application.

[0077] Figure 3 This is a structural block diagram illustrating a remote sensing extraction device for low-canopy-density single-tree-scale young forests that integrates two-dimensional and three-dimensional features, according to an embodiment of this application. Detailed Implementation

[0078] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0079] The following, in conjunction with the accompanying drawings, provides a detailed description of the remote sensing extraction scheme for low-canopy-density single-tree-scale young forests, which integrates two-dimensional and three-dimensional features, through specific embodiments and application scenarios.

[0080] As attached Figure 1 As shown in the figure, the remote sensing extraction method for low canopy closure single-tree scale young forests that integrates two-dimensional and three-dimensional features according to the embodiments of this application includes the following steps:

[0081] Step 101: Acquire satellite remote sensing images, and generate multi-view epipolar images based on the two-line array stereo images and RPC parameters contained in the satellite remote sensing images.

[0082] The remote sensing extraction method for low-canopy-density single-tree-scale young forests with fused two-dimensional and three-dimensional features provided in this application embodiment can be applied to electronic devices, which include a processor and a storage medium. The storage medium stores a computer program related to the remote sensing extraction of low-canopy-density single-tree-scale young forests with fused two-dimensional and three-dimensional features. The processor of the electronic device executes the computer program to implement the remote sensing extraction process for low-canopy-density single-tree-scale young forests with fused two-dimensional and three-dimensional features.

[0083] The satellite used in the embodiments of this application can be any suitable type of satellite, such as the Gaofen-7 (GF-7) satellite.

[0084] In one optional embodiment, the process of acquiring satellite remote sensing imagery and generating multi-view epipolar imagery based on the two-line array stereo imagery contained in the satellite remote sensing imagery and the RPC (Rational Polynomial Coefficient) parameters of the satellite imagery may include the following sub-steps:

[0085] Sub-step 1011: Acquire satellite remote sensing images.

[0086] Sub-step 1012: Use the satellite back-view image contained in the satellite remote sensing image as the main image and the satellite front-view image as the auxiliary image, and generate the projection connection point between the main image and the auxiliary image.

[0087] In one feasible embodiment, using the satellite back-view image contained in the satellite remote sensing image as the main image and the satellite forward-view image as the auxiliary image, and generating the projection tie point between the main image and the auxiliary image may include the following steps:

[0088] Sub-step 1: Use the satellite rearview image as the main image and the satellite forwardview image as the auxiliary image;

[0089] Sub-step 2: Using the auxiliary image as the target auxiliary image, determine the import order of the main image and the auxiliary image, and calculate the average epipolar tilt angle; wherein, the main image is imported first as the reference, and the auxiliary image is imported later and matched with the main image;

[0090] Sub-step 3: Determine the range of the multi-view epipolar image based on the average epipolar tilt angle;

[0091] Sub-step 4: Generate projection connection points between the main image and the auxiliary image that cover the range of the main image quasi-core line image.

[0092] In one optional embodiment, the projection connection point between the main image and the auxiliary image can be calculated according to the following formula:

[0093] P0 = n·epioff + ΔW

[0094] P1=Hn·epioff·tan(θ)+ΔH

[0095] Wherein, ΔH is the outward expansion range of the main image quasi-epidial image relative to the main image in the row direction, ΔW is the outward expansion range of the main image quasi-epidial image relative to the main image in the column direction, θ is the average epipolar tilt angle, P0 is the column coordinate of the starting point of the nth projection trajectory on the main image, P1 is the row coordinate of the starting point of the nth projection trajectory on the main image, H is the length of the main image, and n is the projection trajectory number.

[0096] Sub-step 1013: Based on the projection connection points between the main image and the auxiliary image, establish a two-way coordinate transformation model between the quasi-piston line image and the original image.

[0097] Among them, the original image refers to the main image and the auxiliary image, and each original image corresponds to one of the quasi-pixel images in the multi-view quasi-pixel image.

[0098] In one feasible embodiment, establishing a two-way coordinate transformation model between the quasi-pixel line image and the original image based on the projection connection points between the main image and the auxiliary image may include the following steps:

[0099] S1: Calculate the offset ratio coefficients of the auxiliary image along the direction of the quasi-pixel image and the auxiliary image perpendicular to the direction of the quasi-pixel image;

[0100] S2: Calculate the tilt parameters of the projection trajectory;

[0101] Among them, the projection trajectory tilt parameter represents the angle between the projection trajectory line and the original image line direction;

[0102] S3: Based on the offset scaling factor and the projection trajectory tilt parameter, establish a two-way coordinate transformation model between the quasi-core line image and the corresponding original image.

[0103] The specific steps involved in establishing a two-way coordinate transformation model between the epipolar image and the corresponding original image include the following:

[0104] First, calculate a function with the row coordinates of the projection connection points and the coordinates of the starting point of each projection trajectory, the row coordinates of the image, and the column coordinates of the image as variables:

[0105] y(i,j)=y(i,j0)+a0+a1·x(i,j)+a2·x(i,j)·x(i,j)+a3·y(i,j0)+a3·y(i,j0)·y(i,j0)+a5·x(i,j)·y(i,j0)

[0106] In the formula, i = 1 and 2 are the sequence numbers of the main image and the auxiliary image, y(i,j) is the row coordinate of the j-th projection connection point of image i, x(i,j) is the column coordinate of the j-th projection connection point of image i, j0 represents the starting point number of the projection trajectory where the j-th projection connection point is located, y(i,j0) is the row coordinate of the starting point of the projection trajectory where the j-th projection connection point of image i is located, and a0, a1...a5 are the transformation coefficients to be determined for each original image, which are obtained by the least squares method from the projection connection points on each image.

[0107] The column and row coordinates of the first projection connection point of each original image are represented by CO. i rO i Let i = 1, and 2 be the original image sequence numbers. Let r and c be the row and column coordinates of any point in the epipolar image, then the corresponding row and column coordinates r in the original image are... o c o The calculation method is as follows:

[0108] r o =r·sclc i ·slp+cO i -c·sclr i

[0109] c o =o c +a0+a1·r o +a1·r o ·r o +a3·o c +a4·o c ·o c +a5·r o ·o c

[0110] Among them, o c =c0 i +r·sclc i sclr i and sclr i (i = 1, 2) represent the average scaling factors of the i-th original image along the projection trajectory direction and perpendicular to the projection trajectory direction, respectively; slp is the projection trajectory tilt parameter; r and c are the row and column coordinates of any point in the quasi-epidural image, respectively. o c o Let cO be the row and column coordinates of any point in the quasi-epidural image on the corresponding original image. i Let a0, a1...a5 be the column coordinates of the first projection connection point of the i-th original image, and let a0, a1...a5 be the transformation coefficients to be determined.

[0111] Sub-step 1014: Based on the quasi-epidial image and the two-way coordinate transformation model, the image plane coordinates of the quasi-epidial image are reversed to obtain the image plane coordinates of the original image;

[0112] Sub-step 1015: Calculate the RPC parameters of the multi-view epipolar image based on the image plane coordinates, RPC parameters, and object elevation of the original image.

[0113] Sub-step 1016: Based on the RPC parameters and coordinate transformation model of the multi-view epipolar image, calculate the coordinates of the original image points corresponding to the image points of the multi-view epipolar image point by point;

[0114] Sub-step 1017: Obtain the gray level of the original image points through bilinear interpolation as the gray level of the multi-view epipolar image to generate the multi-view epipolar image.

[0115] Step 102: Generate a pyramid image for the multi-view epipolar image and determine the disparity map of the two linear array stereo images based on the pyramid image.

[0116] In one alternative embodiment, generating a pyramid image for a multi-view epipolar image and determining a disparity map of two linear array stereo images based on the pyramid image may include the following sub-steps:

[0117] Sub-step 1021: Match the top pyramid image of the quasi-nuclear line image of the satellite two-line array stereo image to obtain the disparity map of the top pyramid image;

[0118] Sub-step 1022: Match the pyramid images of the multi-view collimator images of the satellite two-line array stereo image, except for the top and bottom layers, to obtain the disparity map of the other pyramid images;

[0119] Sub-step 1023: Match the bottom pyramid image of the quasi-nuclear line image of the satellite two-line array stereo image to obtain the disparity map of the top pyramid image.

[0120] Step 103: Calculate the point cloud based on the disparity map, and generate a digital surface model (DSM) by meshing the calculated point cloud.

[0121] In an optional embodiment, the step of calculating a point cloud based on a disparity map and generating a digital surface model (DSM) by meshing the calculated point cloud may include the following sub-steps:

[0122] Sub-step 1031: Calculate cloud points based on the disparity map;

[0123] Sub-step 1032: Calculate the planar coordinate range of all point clouds, and based on the planar coordinate range, calculate the length and width of the DSM to be generated;

[0124] The plane coordinate range includes minimum longitude, maximum longitude, minimum latitude, and maximum latitude.

[0125] Sub-step 1033: Based on the length and width of the DSM to be generated and the three-dimensional coordinates of each point cloud, calculate the row and column coordinates of the DSM grid points corresponding to each point cloud, and calculate the latitude, longitude coordinates and elevation of the DSM grid points according to the member points of each DSM grid point.

[0126] Sub-step 1034: Organize the data into a regular grid based on the latitude, longitude coordinates and elevation of each DSM grid point to generate a digital surface model (DSM).

[0127] Step 104: Optimize the DSM using satellite laser altimetry data to obtain the optimized DSM.

[0128] In one optional embodiment, the step of optimizing the DSM using satellite laser altimetry data to obtain the optimized DSM may include the following sub-steps:

[0129] Sub-step 1041: Combine the field inspection points to select the laser points in the satellite laser altimetry data as elevation control points, and solve for the elevation values ​​of the DSM grid points corresponding to the geographical locations of the control points and the elevation difference Δh between the control points;

[0130] A feasible method for determining the elevation difference Δh between the DSM grid points corresponding to the geographical locations of the control points may include the following sub-steps:

[0131] First, based on the laser point and field detection point in the laser altimetry data, calculate the elevation difference dH1 between the laser point and the average elevation of all detection points, calculate the elevation difference dH2 between the laser point and the average elevation of detection points within a 20m radius, and calculate the elevation difference dH3 between the laser point and the nearest detection point.

[0132] Secondly, using publicly available SRTM (Shuttle Radar Topography Mission) DEM data, gross errors were removed from the laser points in the laser altimetry data. Coarse-grained screening of distance attribute parameters, fine screening of echo waveform characteristic parameters, and ground feature screening using GlobeLand30 (30m global land cover data) were performed. The screening results were marked with m_ECP_Flag as the elevation quality control indicator. m_ECP_Flag=1 indicates that the laser point is very suitable for use as an elevation control point, m_ECP_Flag=2 indicates that the laser point can be used as an elevation control point, and m_ECP_Flag=other indicates that it is recommended to reduce the weight of the laser point when using it as an elevation control point.

[0133] Finally, the average value ME, the mean error RMSE, and the number N of dH1, dH2, and dH3 are calculated based on the statistical results of m_ECP_Flag. Based on ME and RMSE, it is determined whether the laser points with m_ECP_Flag=1 and m_ECP_Flag=2 can be used as elevation control points for subsequent experiments.

[0134] Sub-step 1042: Fit the functional relationship between the elevation difference and the latitude and longitude of the DSM grid points corresponding to the control points using the least squares method;

[0135] This functional relationship can be expressed as:

[0136] Δh = f(x, y)

[0137] Sub-step 1043: Substitute the latitude and longitude coordinates of each grid point on the DSM into the function relationship to obtain the corrected elevation difference and obtain the corrected digital surface model;

[0138] Specifically, the corrected elevation difference Δh′ can be obtained by substituting the DSM latitude and longitude x and y into the functional relationship f(x, y), and the corrected digital surface model (H) can be obtained. DSM );

[0139] H DSM =DSM+Δh′

[0140] Sub-step 1044: Arrange the elevation differences Δh in order of magnitude and take the median value as the correction value;

[0141] Sub-step 1045: Add the correction value to the elevation value of each grid point in the corrected digital surface model to obtain the optimized DSM(H′). DSM The optimized model can be represented as:

[0142] H′ DSM =H DSM +Median

[0143] Step 105: Use the isolated forest algorithm to perform image analysis on the optimized DSM and extract pixels containing young forests.

[0144] In one optional embodiment, the method of using the isolated forest algorithm to perform image analysis on the optimized DSM and extracting pixels containing young forests may include the following sub-steps:

[0145] Sub-step 1051: Treat the optimized DSM as a raw dataset.

[0146] Let the original dataset X = {Xn}, where n takes the value from 1 to N. The dataset contains N pixels, and each pixel corresponds to a DSM value.

[0147] Sub-step 1052: Randomly and without replacement draw m data from the original dataset as training data, and put the training data into the root node of a binary tree;

[0148] Sub-step 1053: If the number of pixels in the current node is 1, or the preset maximum tree height is reached, mark the node as a leaf node and stop splitting; otherwise, randomly select a splitting value P.

[0149] Among them, the splitting value is generated between the maximum and minimum DSM values of all pixels in the current node.

[0150] Sub-step 1054: Divide the pixels of the current node into two subsets according to the splitting value. The pixels with DSM values less than the splitting value enter the left subtree, and the rest enter the right subtree.

[0151] Specifically, the data with Q < P can be put into the left subtree F l , and the data with Q ≥ P can be put into the right subtree F r .

[0152] Sub-step 1055: Recursively execute the above splitting process on the left and right subtree nodes to construct an isolation tree, that is, obtain the target tree (iTree tree);

[0153] Sub-step 1056: According to the isolation tree, calculate the path length from the root node to the leaf node for each pixel in the original dataset, and calculate the anomaly score for each pixel based on the path length.

[0154] Specifically, each value in the original dataset X can be traversed through the iTree tree, and the path length L(X i ) during the traversal process can be calculated. Finally, the anomaly score s(X i , Ψ) of X i is calculated according to the following formula:

[0155]

[0156] Among them, c(Ψ) is used to calculate the average path length of the binary search tree, and its function is to normalize the result; the calculation method of H(ψ) is: H(ψ) = ln(i) + ξ (where the Euler constant ζ ≈ 0.58); E(L(X i )) is the average path length of X i in the iTree tree.

[0157] Sub-step 1057: According to the anomaly score, determine the pixels with anomaly scores higher than the preset threshold, or the pixels with anomaly scores ranked in the top preset percentage as the pixels containing young forests.

[0158] In one feasible embodiment, whether a pixel contains young forest can be determined as follows:

[0159] When a certain pixel corresponds to s(x) i When s(x)≈1, the pixel is definitely an outlier, meaning the pixel contains young trees; when s(x)≈1, the pixel contains young trees. i When Ψ)≈0.5, the pixel is almost not an outlier, meaning the pixel contains almost no young trees; when s(x)≈0.5, the pixel is almost not an outlier. i When ψ) << 0.5, the pixel is definitely not an outlier, that is, the pixel does not contain young forest.

[0160] Step 106: Fuse the panchromatic backview image and multispectral image of the satellite remote sensing image to generate a panchromatic-multispectral fused image.

[0161] Step 107: Use a shadow detection algorithm to extract shadows from the panchromatic-multispectral fusion image to obtain the target shadow area.

[0162] In one optional embodiment, the method of extracting shadows from a panchromatic-multispectral fused image using a shadow detection algorithm to obtain the target shadow region may include the following sub-steps:

[0163] Sub-step 1071: Subtract the blue band of the panchromatic-multispectral fusion image from the green band of the panchromatic-multispectral fusion image to obtain a blue-green band difference image;

[0164] Sub-step 1072: Use the histogram thresholding method to extract the shadow region from the near-infrared band image of the panchromatic-multispectral fusion image to obtain the first shadow region extraction result;

[0165] Sub-step 1073: Use the histogram thresholding method to extract the shadow region from the blue-green band difference image, and perform median filtering and boundary tracking on the extraction results to obtain the boundary line;

[0166] Sub-step 1074: Construct a buffer of a specific width based on the boundary line to obtain the extraction result of the second shaded area;

[0167] Sub-step 1075: Overlay the results of the first shadow region extraction and the second shadow region extraction to obtain the shadow boundary region;

[0168] Sub-step 1076: Fill the shadow boundary area to obtain the target shadow area.

[0169] Step 108: Compare the pixels containing the young forest with the target shadow area to determine the target young forest extraction result.

[0170] If there is no shadow area near a pixel, the pixel is considered to have been extracted incorrectly; conversely, if there is a shadow area near a pixel, the pixel is considered to have been extracted correctly.

[0171] In another alternative embodiment, the target young forest can be extracted in the following manner to reduce the false extraction rate of young forests, specifically including the following steps:

[0172] S1: Perform orthorectification and radiometric calibration preprocessing on the satellite panchromatic backsight image; perform RPC orthorectification, radiometric calibration and atmospheric correction preprocessing on the satellite multispectral image; and use the Gram-Schmidt orthogonalization method to fuse the processed satellite panchromatic image and the processed multispectral image to finally generate a panchromatic-multispectral fused image.

[0173] S2, subtract the blue band from the green band of the panchromatic-multispectral fusion image to form a blue-green band difference image;

[0174] S3, extract the shadow region from the near-infrared band image of the panchromatic-multispectral fusion image and the blue-green light band difference image respectively using the histogram thresholding method;

[0175] S4, perform median filtering on the shadow region extraction results of the blue-green light band difference image based on the histogram thresholding method using an n×n window, and trace the boundary of the shadow region to obtain the boundary line;

[0176] Where n is an odd number and n≥3;

[0177] S5. Based on the shadow area boundary tracking results of the blue-green light band difference image, a buffer is made with a radius R (R > n × k, k is the pixel width);

[0178] S6, the shadow region extraction results of the blue-green light band difference image are superimposed with the shadow region extraction results of the near-infrared band image based on the histogram thresholding method to obtain the shadow boundary region, and then the shadow boundary region is filled to obtain the target shadow region;

[0179] S7 combines the pixel extraction results containing young forests with the extraction results of the target shadow area for judgment. If there is no shadow area near the pixel extraction result, the pixel extraction is considered to be incorrect.

[0180] Step S1, which uses the Gram-Schmidt method for image fusion, specifically includes the following sub-steps:

[0181] S1.1 Calculate the overlapping area of ​​the satellite-processed multispectral and panchromatic images, and register and crop the images to make them the same size;

[0182] S1.2, the multispectral image processed by the satellite is weighted according to the spectral response function by a certain weight ω. i A simulated panchromatic band image is generated, which is a simulated low spatial resolution panchromatic band image. The formula for calculating the gray level P of the simulated panchromatic band image is as follows:

[0183]

[0184] Where k is the total number of bands in the satellite-processed multispectral image, and B i is the gray value of the i-th band.

[0185] S1.3, using the simulated panchromatic image as the first component of the Gram-Schmidt orthogonalization, the Gram-Schmidt orthogonalization algorithm is used to perform orthogonal transformation on each band of the satellite-processed multispectral image. The calculation formula is as follows:

[0186]

[0187] Among them, GS T B represents the T-th orthogonal component after Gram-Schmidt orthogonalization transformation. T This indicates that the T-band of the multispectral image processed by GF-7, μ T σ(B) represents the mean gray value of the T-band pixel in the multispectral image after GF-7 processing. T GS l The ) indicates that the T-band of the multispectral image processed by GF-7 is compared with the GS band. l The covariance between them, σ(GS) l GS l ) 2 Indicates GS l The variance is given by i and j, where i and j represent the number of rows and columns of the multispectral image after GF-7 processing, respectively, and M and N represent the number of rows and columns of the entire image.

[0188] S1.4 Replace the first component of the Gram-Schmidt orthogonalization with the satellite-processed panchromatic image, that is, replace the simulated panchromatic band image, to generate a new dataset;

[0189] S1.5, perform Gram-Schmidt inverse transform on the new dataset to generate a panchromatic-multispectral fused image. The Gram-Schmidt inverse transform formula is as follows:

[0190]

[0191] in, This is the T-band of the panchromatic-multispectral fusion image.

[0192] The remote sensing extraction method for low-canopy-density single-tree-scale young forests provided in this application embodiment has two main aspects. First, it utilizes DSM (Digital Sound Mapping) of densely matched stereo images optimized from satellite laser altimetry data to extract tree height information, which can eliminate the interference of low-lying herbaceous plants and shrubs on the extraction of young forests. Second, by fusing tree height with the shadow information generated by trees under low canopy closure, the problem of two-dimensional small target detection is transformed into the relatively easier two-dimensional tree shadow detection and three-dimensional tree height feature extraction, which can further improve the accuracy of single-tree-scale young forest extraction.

[0193] The following is combined with Figure 2 This paper uses a specific example to illustrate the remote sensing extraction method for low-canopy-density single-tree-scale young forests that integrates two-dimensional and three-dimensional features, as provided in this application.

[0194] The method for extracting young forests at the single-tree scale with low canopy closure by integrating two-dimensional and three-dimensional features provided in this embodiment includes the following three steps:

[0195] S1: Acquire remote sensing images from the Gaofen-7 (GF-7) satellite and optimize the densely matched digital surface model (DSM) of the stereo images using laser altimetry data from the GF-7 satellite.

[0196] This step mainly involves optimizing the digital surface model (DSM), which consists of two parts: 1. Image acquisition and processing; 2. DSM optimization. These two parts will be explained in detail below.

[0197] 1. Image Acquisition and Processing

[0198] This study acquires remote sensing imagery from the Gaofen-7 (GF-7) satellite. Based on the GF-7 two-line array stereo imagery and its RPC parameters, a multi-parallel epipolar imagery is generated, along with a corresponding pyramid imagery. The pyramid imagery is then matched to obtain a disparity map of the GF-7 two-line array stereo imagery. Point clouds are then calculated based on the disparity map, and a Data Model (DSM) is generated through point cloud meshing. The purpose of choosing GF-7 satellite remote sensing imagery is as follows: In recent years, some scholars have attempted to use satellite remote sensing data to detect tree mortality at the single-tree scale and estimate survival rates. However, such methods generally require the spatial resolution of the satellite imagery to be at least smaller than the tree crown diameter, placing high demands on spatial resolution. Specifically, a spatial resolution of approximately 0.5 meters is generally required at the single-tree scale. Estimating the survival rate of even smaller saplings requires even higher spatial resolution, making both satellite and UAV data prohibitively expensive. Therefore, how to utilize existing high spatial resolution (approximately 1 meter) satellite remote sensing data to identify dead saplings at the single-tree scale and estimate their survival rate before reaching maturity or the effective afforestation period has significant theoretical and practical implications. Given that GF-7 satellite imagery has high resolution and moderate cost, and that a 1-meter resolution digital surface model (DSM) can be constructed based on its two-line array stereo panchromatic imagery, it was selected as the target image.

[0199] The resolution of the GF-7 satellite remote sensing images used in this specific example is: 0.65m resolution for panchromatic (back-view) two-line array stereo image; and 0.80m resolution for panchromatic (forward-view) two-line array stereo image.

[0200] 2. DSM optimization

[0201] When optimizing the median model of a densely matched digital surface model (DSM) based on GF-7 laser altimetry data, the process begins by selecting laser points from the laser altimetry data as elevation control points based on field detection points, and then calculating the elevation difference Δh between the DSM and the control points. Next, the least squares method is used to fit the coefficient f(x,y) between Δh and the latitude and longitude x and y of the DSM. Substituting x and y into the fitted equation yields the corrected elevation difference Δh′, resulting in the corrected digital surface model (H). DSM Next, the elevation differences Δh are arranged in order of magnitude, and the median is taken as the correction value. Finally, the latitude and longitude x and y are substituted into the coefficient f(x,y), and optimized using the median to obtain the optimized digital surface model (H′). DSM ).

[0202] When screening laser points in laser altimetry data as elevation control points, the height differences dH1 between the laser points in the laser altimetry data and the average elevation of all inspection points, dH2 between the laser points and the average elevation of inspection points within a range of 20 m, and dH3 between the laser points and the elevation value of the nearest inspection point shall be calculated respectively based on the laser points in the laser altimetry data and the field inspection points. Then, the Shuttle Radar Topography Mission (SRTM) DEM data of the global public version is applied to eliminate gross errors of the laser points in the laser altimetry data, conduct coarse-grained screening of ranging attribute parameters, conduct fine screening of echo waveform feature parameters, and conduct feature screening of GlobeLand30 (30 m global land cover data). Finally, the average value, medium error, and quantity of dH1, dH2, and dH3 are calculated according to the screening results, and it is determined whether the screening results can be used as elevation control points for subsequent tests based on the average value and medium error.

[0203] S2. With the help of an outlier detection algorithm, pixels containing young forests are extracted.

[0204] This step mainly uses an outlier detection algorithm to extract pixels of young forests. The specific extraction method is as follows:

[0205] The isolated forest algorithm is used to analyze the optimized DSM image to extract pixels containing young forests. The specific method is as follows: Let the original dataset X = {Xn}, where n ranges from 1 to N, and the dataset contains N pixels, and each pixel has a uniquely corresponding DSM value; randomly select m data from the original dataset as training data and put them into the root node of a binary tree; randomly select a splitting point P (P is generated between the maximum and minimum values of the original dataset); generate a hyperplane according to the splitting point P to divide the data space of the original dataset into two subspaces (put data with Q < P into the left subtree Fl and data with Q ≥ P into the right subtree Fr); recursively repeat the first two steps in the child nodes to continuously construct new child tree nodes until there is only one data in the child tree node or the child tree node has reached the specified height and no longer continues to split, obtaining the iTree tree.

[0206] Then, each value in the original dataset X is traversed through the iTree tree, and the path length L(X i ), is calculated. Finally, the outlier score s(x i , Ψ) of X is calculated according to the path length. The specific formula is as follows: i

[0207]

[0208] Where c(Ψ) is used to calculate the average path length of the binary search tree, and its function is to normalize the result; H(Ψ) is calculated as: H(Ψ)=ln(i)+ζ (where Euler's constant ζ≈0.58); E(L(X i )) is X i Average path length in the iTree tree;

[0209] Finally, based on the anomaly score s(x) i , Ψ) to make a judgment: when s(x i When s(x,Ψ)≈1, the point must be an outlier, meaning the pixel contains young trees; when s(x,Ψ)≈1, the point must be an outlier. i When Ψ)≈0.5, this point is almost not an outlier, meaning that the pixel contains almost no young forests; when s(x)≈0.5, the pixel is almost not an outlier. i When Ψ) << 0.5, the point is definitely not an outlier, meaning that the pixel does not contain young trees.

[0210] S3 combines panchromatic-multispectral fusion imagery from the GF-7 satellite with a shadow detection algorithm to extract shadow areas, thereby reducing the error extraction rate of young forests and outputting young forest extraction results at the single-tree scale.

[0211] This step consists of three parts: 1. Panchromatic-multispectral image fusion; 2. Shadow region extraction; 3. Reducing errors in young forest extraction. The three parts are explained below:

[0212] 1. Panchromatic-Multispectral Image Fusion

[0213] The panchromatic backview image and multispectral image of GF-7 satellite remote sensing imagery are fused through the following steps: the GF-7 panchromatic backview image is preprocessed with orthorectification and radiometric calibration, and the GF-7 multispectral image is preprocessed with RPC (Rational Polynomial Coefficient) orthorectification, radiometric calibration, and atmospheric correction; then, the Gram-Schmidt orthogonalization method is used to fuse the GF-7 processed panchromatic image and multispectral image, finally generating a panchromatic-multispectral fused image.

[0214] 2. Extraction of shaded areas

[0215] The green band of the panchromatic-multispectral fusion image is subtracted from the blue band to form a blue-green band difference image. Then, shadow regions are extracted using the histogram thresholding method based on both the near-infrared image of the panchromatic-multispectral fusion image and the blue-green band difference image. The specific steps for extracting shadow regions using the histogram thresholding method are as follows:

[0216] Next, a median filter is applied to the shadow region extraction results of the blue-green band difference image based on the histogram thresholding method using a window of n×n (n is an odd number and n≥3), and the shadow region boundary is traced to obtain the boundary line; then, based on the shadow region boundary tracing results of the blue-green band difference image, a buffer with a radius R (R>n×k, k is the pixel width) is used.

[0217] Finally, the shadow region extraction results of the blue-green light band difference image and the shadow region extraction results of the near-infrared band image based on the histogram thresholding method are superimposed to obtain the shadow boundary region. Then, the shadow boundary region is filled to obtain the shadow region, which is the target shadow region described in the previous embodiment.

[0218] 3. Reduce erroneous extraction of young forests

[0219] The results of pixel extraction containing young forests are combined with the results of extraction of shaded areas for judgment. If there is no shaded area near the pixel extraction result, the pixel is considered to have been extracted incorrectly.

[0220] This specific example presents a remote sensing extraction method for young forests at the single-tree scale, integrating two-dimensional and three-dimensional features. Firstly, it utilizes DSM (Densely Matched Stereoscopic Imagery) optimized from GF-7 satellite laser altimetry data to extract tree height information, thereby eliminating interference from low-lying herbaceous plants and shrubs in the young forest extraction process. Secondly, by fusing tree height with shadow information generated by trees under low canopy closure, the challenge of two-dimensional small target detection is transformed into the relatively easier task of two-dimensional tree shadow detection and three-dimensional tree height feature extraction, further improving the accuracy of young forest extraction at the single-tree scale.

[0221] Figure 3 The structural block diagram of a remote sensing extraction device for low canopy closure single-tree scale young forest that integrates two-dimensional and three-dimensional features is shown in the embodiment of this application.

[0222] The remote sensing extraction device for low canopy closure single-tree scale young forest that integrates two- and three-dimensional features provided in this application includes the following functional modules;

[0223] The acquisition module 301 is used to acquire satellite remote sensing images and generate multi-view epipolar images based on the two-line array stereo images and RPC parameters contained in the satellite remote sensing images.

[0224] The disparity map determination module 302 is used to generate a pyramid image for the multi-view epipolar image and determine the disparity map of the two linear array stereo images based on the pyramid image.

[0225] The model generation module 303 is used to calculate the point cloud based on the disparity map, and generate a digital surface model (DSM) by meshing the calculated point cloud.

[0226] Optimization module 304 is used to optimize the DSM using satellite laser altimetry data to obtain an optimized DSM;

[0227] The extraction module 305 is used to perform image analysis on the optimized DSM using the isolated forest algorithm to extract pixels containing young forests;

[0228] The fused image generation module 306 is used to fuse the panchromatic backview image and the multispectral image of the satellite remote sensing image to generate a panchromatic-multispectral fused image;

[0229] The shadow region determination module 307 is used to extract shadows from the panchromatic-multispectral fused image using a shadow detection algorithm to obtain the target shadow region;

[0230] The result determination module 308 is used to compare the pixel containing the young forest with the target shadow area to determine the target young forest extraction result; wherein, if there is no shadow area near the pixel, the pixel extraction is determined to be incorrect.

[0231] Optionally, the acquisition module includes:

[0232] The first submodule is used to acquire satellite remote sensing images;

[0233] The second submodule is used to take the satellite rearview image contained in the satellite remote sensing image as the main image and the satellite frontview image as the auxiliary image, and generate the projection connection point between the main image and the auxiliary image.

[0234] The third submodule is used to establish a two-way coordinate transformation model between the quasi-pixel image and the original image based on the projection connection point between the main image and the auxiliary image. The original image refers to the main image and the auxiliary image, and each original image corresponds to one quasi-pixel image in the multi-view quasi-pixel image.

[0235] The fourth submodule is used to perform reverse calculation on the image plane coordinates of the quasi-epidial image based on the quasi-epidial image and the bidirectional coordinate transformation model to obtain the image plane coordinates of the original image.

[0236] The fifth submodule is used to calculate the RPC parameters of the multi-view epipolar image based on the image plane coordinates, RPC parameters, and object elevation of the original image.

[0237] The sixth submodule is used to calculate the coordinates of the original image points corresponding to the image points of the multi-view epipolar image point by point, based on the RPC parameters of the multi-view epipolar image and the coordinate transformation model.

[0238] The seventh submodule is used to obtain the gray level of the original image points through bilinear interpolation as the gray level of the multi-view epipolar image, so as to generate the multi-view epipolar image.

[0239] Optionally, the second submodule is specifically used for:

[0240] Use satellite rear-view images as the main images and satellite forward-view images as auxiliary images;

[0241] Using the auxiliary image as the target auxiliary image, the import order of the main image and the auxiliary image is determined, and the average epipolar tilt angle is calculated; wherein, the main image is imported first as a reference, and the auxiliary image is imported later and matched with the main image;

[0242] The range of the multi-view epipolar image is determined based on the average epipolar tilt angle.

[0243] Generate projection connection points between the main image and the auxiliary image that cover the range of the main image quasi-piezoresistive image.

[0244] Optionally, the third submodule is specifically used for:

[0245] Calculate the offset ratio coefficients of the auxiliary image along the direction of the quasi-core image and the auxiliary image perpendicular to the direction of the quasi-core image;

[0246] Calculate the projection trajectory tilt parameter, wherein the projection trajectory tilt parameter represents the angle between the projection trajectory line and the original image line direction;

[0247] Based on the offset ratio coefficient and the projection trajectory tilt parameter, a two-way coordinate transformation model between the quasi-core line image and the corresponding original image is established.

[0248] Optionally, the disparity map determination module includes:

[0249] The eighth submodule is used to match the top pyramid image of the quasi-nuclear line image of the satellite two-line array stereo image to obtain the disparity map of the top pyramid image.

[0250] The ninth submodule is used to match the pyramid images of other layers besides the top and bottom layers of the multi-view collimator image of the satellite two-line array stereo image to obtain the disparity map of the other pyramid images.

[0251] The tenth submodule is used to match the bottom pyramid image of the quasi-nuclear line image of the satellite two-line array stereo image to obtain the disparity map of the top pyramid image.

[0252] Optionally, the model generation module includes:

[0253] The eleventh submodule is used to calculate cloud points based on the disparity map;

[0254] The twelfth submodule is used to calculate the planar coordinate range of all point clouds, and based on the planar coordinate range, calculate the length and width of the DSM to be generated; wherein, the planar coordinate range includes minimum longitude, maximum longitude, minimum latitude, and maximum latitude;

[0255] The thirteenth submodule is used to calculate the row and column coordinates of the DSM grid points corresponding to each point cloud based on the length and width of the DSM to be generated and the three-dimensional coordinates of each point cloud, and to calculate the latitude and longitude coordinates and elevation of the DSM grid points according to the membership points of each DSM grid point.

[0256] The fourteenth submodule is used to organize regular grid data based on the latitude, longitude coordinates and elevation of each DSM grid point to generate the Digital Surface Model (DSM).

[0257] Optionally, the optimization module is specifically used for:

[0258] By combining the laser points in the satellite laser altimetry data selected from the field detection points, the elevation values ​​of the DSM grid points corresponding to the geographical locations of the control points and the elevation differences between the control points are calculated.

[0259] The least squares method was used to fit the functional relationship between the elevation difference and the latitude and longitude of the DSM grid points corresponding to the control points:

[0260] Substitute the latitude and longitude coordinates of each grid point on the DSM into the function relationship to obtain the corrected elevation difference, and obtain the corrected digital surface model.

[0261] Arrange the elevation differences in order of magnitude and take the median value as the correction value;

[0262] The elevation value of each grid point in the corrected digital surface model is added to the correction value to obtain the optimized DSM.

[0263] Optionally, the extraction module is specifically used for:

[0264] The optimized DSM is regarded as a raw dataset, wherein the raw dataset includes multiple cells, and each cell corresponds to a DSM value;

[0265] m data points are randomly selected without replacement from the original dataset as training data, and the training data are placed in the root node of a binary tree;

[0266] If the number of cells in the current node is 1, or the preset maximum tree height has been reached, then the node is marked as a leaf node and the segmentation stops; otherwise, a segmentation value is randomly selected, wherein the segmentation value is generated between the maximum and minimum values ​​of the DSM values ​​of all cells in the current node.

[0267] Based on the segmentation value, the cells of the current node are divided into two subsets: cells with a DSM value less than the segmentation value enter the left subtree, and the rest enter the right subtree.

[0268] Recursively perform the above segmentation process on the left and right subtree nodes to construct an isolated tree;

[0269] Based on the isolated tree, calculate the path length from the root node to the leaf node for each cell in the original dataset, and calculate the anomaly score for each cell based on the path length.

[0270] Based on the abnormal scores, pixels with abnormal scores higher than a preset threshold, or pixels with abnormal scores ranking in the top preset percentage, are identified as pixels containing young forests.

[0271] Optionally, the shadow area determination module is specifically used for:

[0272] The blue-green band difference image is obtained by subtracting the blue band of the panchromatic-multispectral fused image from the green band of the panchromatic-multispectral fused image.

[0273] The shadow region of the near-infrared band image of the panchromatic-multispectral fused image is extracted using the histogram thresholding method to obtain the first shadow region extraction result;

[0274] The shadow region of the blue-green band difference image was extracted using the histogram thresholding method, and the extraction results were then subjected to median filtering and boundary tracing to obtain the boundary line.

[0275] A buffer of a specific width is constructed based on the boundary line to obtain the extraction result of the second shaded area;

[0276] The first shadow region extraction result and the second shadow region extraction result are overlaid to obtain the shadow boundary region;

[0277] The shadow boundary region is filled to obtain the target shadow region.

[0278] The embodiments provided in this application Figure 3 The remote sensing extraction device for low-canopy-density single-tree-scale young forests, which integrates two-dimensional and three-dimensional features, can achieve... Figure 1 The various processes implemented in the method implementation examples will not be described again here to avoid repetition.

[0279] The remote sensing extraction device for low-canopy-density single-tree-scale young forests provided in this application embodiment has two advantages. First, it uses DSM (Digital Sound Mapping) of densely matched stereo images optimized by satellite laser altimetry data to extract tree height information, which can eliminate the interference of low-lying herbaceous plants and shrubs on the extraction of young forests. Second, by fusing tree height with the shadow information generated by trees under low canopy closure, the problem of two-dimensional small target detection is transformed into two-dimensional tree shadow detection and three-dimensional tree height feature extraction, which are relatively easier to solve, thereby further improving the accuracy of single-tree-scale young forest extraction.

[0280] Optionally, embodiments of this application also provide an electronic device, including a processor, a memory, and a program or instructions stored in the memory and executable on the processor. When the program or instructions are executed by the processor, they implement the processes performed by the aforementioned remote sensing extraction device for low canopy closure single-tree scale young forests that integrates two-dimensional and three-dimensional features, and can achieve the same technical effect. To avoid repetition, they will not be described again here.

[0281] Memory, used to store computer programs;

[0282] When the processor executes the program stored in the memory, it implements the remote sensing extraction method for low-canopy-density single-tree-scale young forests that integrates two-dimensional and three-dimensional features, as shown in the above method embodiments.

[0283] The communication bus mentioned above can be a Peripheral Component Interconnect (PCI) bus or an Extended Industry Standard Architecture (EISA) bus, etc. This communication bus can be divided into address bus, data bus, control bus, etc. For ease of illustration, only one thick line is used to represent it in the diagram, but this does not mean that there is only one bus or one type of bus.

[0284] The communication interface is used for communication between the aforementioned terminal and other devices.

[0285] The memory may include random access memory (RAM) or non-volatile memory, such as at least one disk storage device. Optionally, the memory may also be at least one storage device located remotely from the aforementioned processor.

[0286] The processors mentioned above can be general-purpose processors, including central processing units (CPUs), network processors (NPs), etc.; they can also be digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components.

[0287] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0288] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A remote sensing method for extracting young forests at the single-tree scale with low canopy closure by integrating two-dimensional and three-dimensional features, characterized in that, The method includes: Acquire satellite remote sensing images, and generate multi-view epipolar images based on the two-line array stereo images and RPC parameters contained in the satellite remote sensing images; A pyramid image is generated for the multi-view epipolar image, and a disparity map of the two linear array stereo images is determined based on the pyramid image; The point cloud is calculated based on the disparity map, and the calculated point cloud is meshed to generate a digital surface model (DSM). The DSM is optimized using satellite laser altimetry data to obtain the optimized DSM; The optimized DSM is regarded as a raw dataset, wherein the raw dataset includes multiple cells, and each cell corresponds to a DSM value; m data points are randomly selected without replacement from the original dataset as training data, and the training data are placed in the root node of a binary tree; If the number of cells in the current node is 1, or the preset maximum tree height has been reached, then the node is marked as a leaf node and the segmentation stops; otherwise, a segmentation value is randomly selected, wherein the segmentation value is generated between the maximum and minimum values ​​of the DSM values ​​of all cells in the current node. The current node's cells are divided into two subsets based on the segmentation value. Cells with a DSM value less than the segmentation value are entered into the left subtree, and the rest are entered into the right subtree. Recursively perform the above segmentation process on the left and right subtree nodes to construct an isolated tree; Based on the isolated tree, calculate the path length from the root node to the leaf node for each cell in the original dataset, and calculate the anomaly score for each cell based on the path length. Based on the abnormal scores, pixels with abnormal scores higher than a preset threshold, or pixels with abnormal scores ranking in the top preset percentage, are identified as pixels containing young forests. The panchromatic backview image and the multispectral image of the satellite remote sensing image are fused to generate a panchromatic-multispectral fused image; The blue-green band difference image is obtained by subtracting the blue band of the panchromatic-multispectral fused image from the green band of the panchromatic-multispectral fused image. The shadow region of the near-infrared band image of the panchromatic-multispectral fused image is extracted using the histogram thresholding method to obtain the first shadow region extraction result; The shadow region of the blue-green band difference image was extracted using the histogram thresholding method, and the extraction results were then subjected to median filtering and boundary tracing to obtain the boundary line. A buffer of a specific width is constructed based on the boundary line to obtain the extraction result of the second shaded area; The first shadow region extraction result and the second shadow region extraction result are overlaid to obtain the shadow boundary region. The target shadow area is obtained by filling the shadow boundary region. The pixels containing young forests are compared with the target shadow area to determine the target young forest extraction result; if there is no shadow area near the pixel, the pixel extraction is determined to be incorrect.

2. The method according to claim 1, characterized in that, The steps for acquiring satellite remote sensing imagery and generating multi-view epipolar imagery based on the two-line array stereo imagery and RPC parameters contained within the satellite remote sensing imagery include: Acquire satellite remote sensing images; The satellite back-view image contained in the satellite remote sensing image is used as the main image, and the satellite front-view image is used as the auxiliary image, and projection connection points between the main image and the auxiliary image are generated. Based on the projection connection points between the main image and the auxiliary image, a two-way coordinate transformation model between the quasi-pixel image and the original image is established, wherein the original image refers to the main image and the auxiliary image, and each original image corresponds to one quasi-pixel image in the multi-view quasi-pixel image; Based on the quasi-epidial line image and the bidirectional coordinate transformation model, the image plane coordinates of the quasi-epidial line image are reversed to obtain the image plane coordinates of the original image. Based on the image plane coordinates, RPC parameters, and object elevation of the original image, calculate the RPC parameters of the multi-view epipolar image. Based on the RPC parameters of the multi-view epipolar image and the coordinate transformation model, the coordinates of the original image points corresponding to the image points of the multi-view epipolar image are calculated point by point. The gray levels of the original image points are obtained by bilinear interpolation and used as the gray levels of the multi-view epipolar image to generate the multi-view epipolar image.

3. The method according to claim 2, characterized in that, The step of using the satellite back-view image contained in the satellite remote sensing imagery as the main image and the satellite forward-view imagery as the auxiliary imagery, and generating projection connection points between the main imagery and the auxiliary imagery, includes: Use satellite rear-view images as the main images and satellite forward-view images as auxiliary images; Using the auxiliary image as the target auxiliary image, the import order of the main image and the auxiliary image is determined, and the average epipolar tilt angle is calculated; wherein, the main image is imported first as a reference, and the auxiliary image is imported later and matched with the main image; The range of the multi-view epipolar image is determined based on the average epipolar tilt angle. Generate projection connection points between the main image and the auxiliary image that cover the range of the main image quasi-piezoresistive image.

4. The method according to claim 2, characterized in that, The steps of establishing a two-way coordinate transformation model between the quasi-piston line image and the original image based on the projection connection points between the main image and the auxiliary image include: Calculate the offset ratio coefficients of the auxiliary image along the direction of the quasi-core image and the auxiliary image perpendicular to the direction of the quasi-core image; Calculate the projection trajectory tilt parameter, wherein the projection trajectory tilt parameter represents the angle between the projection trajectory line and the original image line direction; Based on the offset ratio coefficient and the projection trajectory tilt parameter, a two-way coordinate transformation model between the quasi-core line image and the corresponding original image is established.

5. The method according to claim 1, characterized in that, The steps of generating a pyramid image for the multi-view epipolar image and determining the disparity map of the two linear array stereo images based on the pyramid image include: Match the top pyramid image of the quasi-nuclear line image of the satellite two-line array stereo image to obtain the disparity map of the top pyramid image; The disparity map of the other pyramid images, excluding the top and bottom layers of the multi-view collimator image of the satellite two-line array stereo image, is obtained by matching the pyramid images of the other pyramid images. The bottom pyramid image of the quasi-nuclear line image of the satellite two-line array stereo image is matched to obtain the disparity map of the top pyramid image.

6. The method according to claim 1, characterized in that, The step of calculating the point cloud based on the disparity map and generating a digital surface model (DSM) by meshing the calculated point cloud includes: Calculate cloud points based on the disparity map; Calculate the planar coordinate range of all point clouds, and based on the planar coordinate range, calculate the length and width of the DSM to be generated; wherein the planar coordinate range includes minimum longitude, maximum longitude, minimum latitude, and maximum latitude; Based on the length and width of the DSM to be generated, and the three-dimensional coordinates of each point cloud, calculate the row and column coordinates of the DSM grid points corresponding to each point cloud, and calculate the latitude, longitude coordinates and elevation of the DSM grid points according to the member points of each DSM grid point. The digital surface model (DSM) is generated by organizing the latitude, longitude, and elevation coordinates of each DSM grid point into a regular grid.

7. The method according to claim 1, characterized in that, The steps for optimizing the DSM using satellite laser altimetry data to obtain the optimized DSM include: By combining the laser points in the satellite laser altimetry data selected from the field detection points, the elevation values ​​of the DSM grid points corresponding to the geographical locations of the control points and the elevation differences between the control points are calculated. The least squares method was used to fit the functional relationship between the elevation difference and the latitude and longitude of the DSM grid points corresponding to the control points: Substitute the latitude and longitude coordinates of each grid point on the DSM into the function relationship to obtain the corrected elevation difference, and obtain the corrected digital surface model. Arrange the elevation differences in order of magnitude and take the median value as the correction value; The elevation value of each grid point in the corrected digital surface model is added to the correction value to obtain the optimized DSM.

8. A remote sensing extraction device for low-canopy-density single-tree-scale young forests integrating two-dimensional and three-dimensional features, characterized in that, The device includes: The acquisition module is used to acquire satellite remote sensing images and generate multi-view epipolar images based on the two-line array stereo images and RPC parameters contained in the satellite remote sensing images. The disparity map determination module is used to generate a pyramid image for the multi-view epipolar image and determine the disparity map of the two linear array stereo images based on the pyramid image. The model generation module is used to calculate point clouds based on the disparity map, and to generate a digital surface model (DSM) by meshing the calculated point clouds. An optimization module is used to optimize the DSM using satellite laser altimetry data to obtain an optimized DSM. The extraction module treats the optimized DSM as a raw dataset, wherein the raw dataset includes multiple pixels, each pixel corresponding to a DSM value. It randomly extracts m data points without replacement from the raw dataset as training data and places these training data into the root node of a binary tree. If the current node has only one pixel or has reached a preset maximum tree height, the node is marked as a leaf node and the segmentation stops; otherwise, a segmentation value is randomly selected, wherein the segmentation value is generated from the maximum and minimum values ​​of the DSM values ​​of all pixels in the current node. Between these steps, the current node's cells are divided into two subsets based on the segmentation value. Cells with a DSM value less than the segmentation value enter the left subtree, and the rest enter the right subtree. The above segmentation process is recursively performed on the left and right subtree nodes to construct an isolated tree. Based on the isolated tree, the path length from the root node to the leaf node is calculated for each cell in the original dataset, and the anomaly score of each cell is calculated based on the path length. Based on the anomaly score, cells with anomaly scores higher than a preset threshold, or cells with anomaly scores ranking in the top preset percentage, are identified as cells containing young forests. The fused image generation module is used to fuse the panchromatic backview image and the multispectral image of the satellite remote sensing image to generate a panchromatic-multispectral fused image; The shadow region determination module is used to subtract the blue band of the panchromatic-multispectral fusion image from the green band of the panchromatic-multispectral fusion image to obtain a blue-green band difference image; to extract shadow regions from the near-infrared band image of the panchromatic-multispectral fusion image using a histogram thresholding method to obtain a first shadow region extraction result; to extract shadow regions from the blue-green band difference image using a histogram thresholding method; to perform median filtering and boundary tracking on the extraction results to obtain a boundary line; to construct a buffer of a specific width based on the boundary line to obtain a second shadow region extraction result; to overlay the first shadow region extraction result and the second shadow region extraction result to obtain a shadow boundary region; and to fill the shadow boundary region to obtain the target shadow region. The result determination module is used to compare the pixel containing the young forest with the target shadow area to determine the target young forest extraction result; wherein, if there is no shadow area near the pixel, the pixel extraction is determined to be incorrect.

9. An electronic device, characterized in that, The electronic device includes a processor, a memory, and a program or instructions stored in the memory and executable on the processor. The program or instructions are executed by the processor using any of the steps of the remote sensing extraction method for low-canopy-density single-tree-scale young forests that integrates two-dimensional and three-dimensional features according to claims 1-7.

Images