A remote sensing-based method for extracting fruit tree density in forest-grain intercropping areas

Abstract

This invention provides a remote sensing-based method for extracting fruit tree density in forest-grain intercropping areas, belonging to the field of fruit tree density extraction. The method comprises several main steps: First, based on instance segmentation, a deep learning model is used to extract plots from ultra-high spatial resolution images as basic computational units, replacing traditional segmentation methods. Then, based on a time-series model, features from Sentinel-2 time-series images are extracted to distinguish ordinary farmland plots from forest-grain intercropping plots, determining the plot units for subsequent fruit tree density extraction. Next, spatial information features of fruit trees are extracted from GF-2 images based on a semantic segmentation model. Finally, a template filling strategy is designed based on the semantic segmentation results to obtain the planting density of fruit trees in the forest-grain intercropping area under the constraints of the plot results. This extraction method has high accuracy, fast processing speed, and good demonstration application results, providing support for the scientific adjustment of vegetation planting structure in forest-grain intercropping areas.

Description

Technical Field

[0001] This invention relates to the field of agricultural remote sensing monitoring. More specifically, it relates to a method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing. Background Technology

[0002] Agroforestry refers to the creation of multi-layered, multi-functional artificial composite communities through different combinations of various plants. It fully utilizes the spatial and temporal differences created by these plant growth processes to maximize the use of limited agricultural resources. Studies have shown that, under appropriate planting densities, agroforestry can increase yields by 20% to 50%. However, improper intercropping methods can have the opposite effect. As fruit tree canopies grow larger each year, high-density fruit orchards can significantly reduce the yield of underlying food crops. Extracting the planting density of fruit trees in agroforestry areas can provide support for the scientific adjustment of vegetation planting structures in these areas.

[0003] The individual tree is the basic unit of a forest, and the key to extracting forest density lies in the number of trees per stand. Since there is a significant difference in height and remote sensing indices between the center of a tree and its crown boundary, these features can be used to locate the center of an individual tree within a grove. Therefore, two common methods for extracting tree density are currently used: local maxima based on crown height and local maxima based on visible light imagery. The former often uses high-precision UAV imagery or radar data to extract tree height and calculate forest density. However, the high manpower and material costs of acquiring large-area data via UAVs limit the feasibility of this method for large-scale forest density mapping. Using satellite remote sensing data, a functional relationship between remote sensing indices and tree location can be constructed to extract the tree center point and calculate its density. However, the construction process of this functional relationship is complex, and its accuracy is difficult to improve. Furthermore, it is more suitable for extracting individual trees in sparsely planted areas and is not applicable to areas with intercropping of crops and forests where fruit trees grow primarily in dense thickets.

[0004] Therefore, providing a method for extracting fruit tree density suitable for forest-grain intercropping areas is an urgent problem to be solved. Summary of the Invention

[0005] This invention provides a method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing, comprising three aspects:

[0006] Firstly, this invention provides a land parcel extraction method based on instance segmentation, which uses a deep learning model to extract land parcels from ultra-high spatial resolution images instead of traditional segmentation methods. These parcels serve as the basic computational unit for extracting fruit tree density in forest-grain intercropping areas, and include the following steps:

[0007] (1) Preprocessing GF-2 (Gaofen-2) data to obtain ultra-high spatial resolution images includes six steps: radiometric calibration, atmospheric correction, orthorectification, image registration, image fusion and image mosaicking.

[0008] (2) Based on the visual characteristics of farmland plots, a farmland zoning and stratification system is established to divide the farmland with forest-grain intercropping into several relatively consistent subtypes.

[0009] (3) Based on the visual features presented by each type of land parcel, select and improve appropriate instance segmentation models to extract each type of farmland in layers, and merge the subclasses of each farmland to obtain a complete farmland information extraction result map.

[0010] The above-mentioned land parcel extraction based on the instance segmentation idea, in step (1), specifically the preprocessing method for obtaining the ultra-high spatial resolution image is as follows:

[0011] (1-1) Radiometric calibration, which converts the data values ​​(DN) recorded by the sensor into apparent reflectance to eliminate the sensor's own errors, requires the calculation of radiance values ​​based on absolute radiometric calibration coefficients. The formula can be expressed as:

[0012] L λ =Gain·DN+Offset

[0013] Among them, L λ To obtain the converted standard radiance, the payload observation values ​​(DN) and absolute calibration coefficient offset values ​​(Offset) of the GF-2 (Gaofen-2) satellite were obtained from the China Center for Resources Satellite Data and Application. The sensor parameter DN values ​​and Offset values ​​used in this invention are shown in Table 1.

[0014] Table 1 Absolute Radiation Calibration Factors

[0015]

[0016] (1-2) Atmospheric correction converts radiance or surface reflectance into actual surface reflectance, effectively eliminating the influence of atmospheric and illumination factors on ground reflection. This invention utilizes the FLAASH model for atmospheric correction. This model algorithm is based on pixel-level correction of the MODTRAN 4+ radiative transfer model, correcting the cascading effects caused by diffuse reflection. It includes classification maps of cirrus clouds and opaque clouds and can adjust spectral smoothing caused by human intervention.

[0017] (1-3) Orthorectification is the process of generating a planar orthorectified image by eliminating the influence of terrain or the distortion caused by camera orientation. In addition to the conventional geometric correction, the distortion caused by terrain undulation is also corrected according to the DEM. This invention uses the orthorectification tool RPC Orthorectification Workflow and the digital elevation model adopts GMTED 2010 DEM. At the same time, in order to perform image fusion later, the panchromatic image of GF-2 is resampled to 1m resolution and the multispectral image is resampled to 4m resolution.

[0018] (1-4) Image registration: One image (reference image) in the same area is used as a calibration copy for another image. Since the panchromatic image and multispectral image of GF-2 do not correspond completely, the multispectral image is selected as the reference image and the panchromatic image is registered. The output resolution is 1m.

[0019] (1-5) Image fusion: The low-resolution multispectral image is fused with the high-resolution panchromatic image to obtain an image with both high spatial resolution and multispectral characteristics. After orthorectification, the GF-2 image has a panchromatic band resolution of 1m and four multispectral band resolutions of 4m. This invention uses the NNDiffuse Pan Sharpening fusion method to fuse multispectral and panchromatic images.

[0020] (1-6) Image mosaicking: Two or more remote sensing images with geographic coordinates are stitched together to obtain a larger and more complete remote sensing image. Before mosaicking the GF-2 image, the present invention needs to perform image registration. The Sentinel-2 remote sensing image is selected as the reference image, and the GF-2 image is used as the image to be registered. Then, the GF-2 image is mosaicked using the Seamless Mosaic tool.

[0021] The above-mentioned land parcel extraction based on the instance segmentation idea still needs to perform "geographic zoning" in step (2) to improve the interpretation effect of remote sensing images. The specific steps are as follows:

[0022] (2-1) With the help of road network / water system / topographic line, the image is divided into naturally aggregated sub-regions by masking to ensure that the land cover type within each sub-region is relatively uniform. This step cuts the image into several discrete geographic blocks, so that the classification process can be limited to the corresponding sub-regions, thereby avoiding the propagation of classification error. At the same time, dividing the image into blocks can be used for parallel computing to improve extraction efficiency.

[0023] (2-2) Under the constraints of each geographic block, the edge of the target feature is extracted in each region. For farmland with blurred edges, the object boundary is obtained based on texture. Based on the slope calculated based on the terrain data, the differences in geometric and texture features of farmland in different regions are expressed.

[0024] The above-mentioned land parcel extraction based on the instance segmentation idea requires, in step (3), the use of the instance segmentation model to extract different types of farmland in layers. The specific steps are as follows:

[0025] (3-1) The instance segmentation algorithm used in this invention is the Hybrid Task Cascade (HTC) model, with the skeleton being ResNeXt-101-64x4d and Deformable ConvNets v2 (DCN). Sixteen 1024*1024 pixel plot samples drawn by visual interpretation are input into the model and trained.

[0026] (3-2) Using a sliding window, many pixel blocks are cropped from the original image from top to bottom and from left to right for instance segmentation;

[0027] (3-3) Divide each pixel block into a target region and an ignored region. The left and top edges of the ignored region are two pixels away from the grid boundary. When the coordinates (x, y) of the top left corner of the detection box fall within the target region, the detection result is retained.

[0028] (3-4) When the sliding window is on the boundary of the original image, the ignored areas on the edges of the grid in the same direction are merged into the target area.

[0029] Secondly, regarding the extraction of forest-grain intercropping plots, this invention, based on a time series model, distinguishes between ordinary farmland plots and forest-grain intercropping plots based on the features of Sentinel-2 time series images, and determines the plots for subsequent fruit tree density extraction tasks, including the following steps:

[0030] 1. Design time series features suitable for major food crops and major fruit trees, and extract different vegetation phenological change features from long-term Sentinel-2 images;

[0031] 2. Based on the irregular time series model, crop classification is carried out, and based on the classification results, discrimination conditions are designed for the existence of only grain crops, the existence of grain crops and fruit trees, and the existence of only fruit trees, so as to divide farmland into ordinary farmland and forest-grain intercropping plots.

[0032] The extraction of the above-mentioned intercropping plots of forest and grain, in step (1), the time series features designed to be suitable for the main grain crops and the main fruit trees are as follows:

[0033] (1-1) This invention selects three vegetation time-series characteristics to distinguish three plot types: Normalized Difference Vegetation Index (NDVI), Normalized Red Edge Index (NDRE1), and Red Edge Normalized Vegetation Index (NDVIre1). The calculation formula is as follows:

[0034]

[0035] In step (2), the vegetation phenological change characteristics are extracted based on an irregular time series model for the above-mentioned forest-grain intercropping plots. The specific steps are as follows:

[0036] (2-1) This invention extracts irregular VI time series from the original optical image for each pixel location in the image and converts it into a "time-feature" matrix to characterize crop growth patterns. In the "time-feature" matrix, the horizontal axis represents time, and the vertical axis represents the value of the vegetation index;

[0037] (2-2) In this invention, the "time-feature" matrix obtained in the first step is input into the three channels of ResNet-18 to obtain a three-dimensional matrix suitable for ResNet-18, thereby extracting vegetation phenological change features.

[0038] The formula for stretching a VI time series into a "time-feature" matrix is ​​as follows:

[0039]

[0040] Among them, VI ori The original value of VI, VI r VI is the stretched value, and length is the length of the time series.

[0041] (2-3) After obtaining the pixel-level classification results based on temporal features, the present invention designed three discrimination conditions to distinguish between ordinary farmland and forest-grain intercropping plots, namely 1. only grain crops exist (grain crop distribution area is greater than 80%), 2. grain crops and fruit trees exist (grain crop distribution area and fruit tree distribution area are greater than 30%), and only fruit trees exist (fruit tree distribution area is greater than 80%), and then screened out forest-grain intercropping plots.

[0042] Thirdly, regarding the extraction of fruit tree density in forest-grain intercropping areas, this invention extracts spatial information features of fruit trees from GF-2 images based on a semantic segmentation model, and designs a template filling strategy based on the semantic segmentation results to obtain the planting density of fruit trees under the constraints of plot results, including the following steps:

[0043] 1. A semantic segmentation model based on convolutional neural networks was used to extract fruit trees from the preprocessed GF-2 image at 1m resolution, and the semantic segmentation results of the fruit trees were obtained.

[0044] 2. Based on the template filling algorithm, different templates are formed according to the spectral characteristics of fruit tree canopies of different sizes. Then, these generated templates are filled into the fruit tree semantic segmentation results to obtain the center position of the canopy of each fruit tree.

[0045] 3. The planting density of fruit trees is obtained by dividing the number of fruit tree center points within each plot element by the element area (mu).

[0046] The extraction of fruit tree density in the above-mentioned forest-grain intercropping area involves extracting fruit trees from the ultra-high resolution image using a semantic segmentation model in step (1). The specific steps are as follows:

[0047] (1-1) The semantic segmentation model used in this invention is the D-LinkNet model. D-LinkNet uses ResNet trained on the ImageNet dataset as the network encoder and adds a dilated-convolution layer with shortcuts in the center part, which makes the entire network more powerful, has a larger receptive field, and integrates multi-scale information. Twenty 512*512 pixel fruit tree samples drawn by visual interpretation are input into the model and trained.

[0048] In the extraction of fruit tree density in the above-mentioned forest-grain intercropping area, in step (2), the template filling algorithm is used to obtain the center position of the canopy of each fruit tree. The specific steps are as follows:

[0049] (2-1) Randomly select the labeling area and label all individual fruit trees to reduce subjective bias. Create a series of square empty templates with odd lengths to ensure that the grid is centered and assign each fruit tree to the template with the closest area.

[0050] (2-4) After aligning the center of the fruit tree with the center of the template, calculate the spectral average value of the connecting components in each template to generate the original template. The original template series shows the typical characteristics of the crown of a single fruit tree.

[0051] (2-5) A mathematical model is constructed based on the typical features of the original template. Preprocessing is used to obtain better performance. Since the near-infrared response of the fruit tree center is stronger, the spectral characteristics of the template pixels are calculated using the following formula:

[0052] R=C

[0053]

[0054] Where R represents the NIR2 band (908nm), B represents the Red band (658.8nm), G represents the Green band (546.2nm), C is a constant, and BG min and BG maxThese are the minimum and maximum values ​​for bands B and G, respectively. n represents the template length, and distance_to_cent represents the distance to the center of the template, in pixels.

[0055] (2-6) Given a four-connected region C of arbitrary shape on the target raster image, the pixel value in C is denoted as x. i,j Where i,j are the coordinates on the digital grid; and several square sub-connected regions S with different side lengths are given. n Where n is the index of the sub-connected region; let y be the pixel value in each sub-connected region. i,j , where i, j are the coordinates of the bounding box grid of the sub-connected components. The sub-connected components S are filled into C. Each sub-connected component can be reused any number of times. The ultimate goal is to minimize the global value G on the digital grid.

[0056] G = sum(|x i,j -y i,j |)+p+q+r

[0057] Where p is the penalty value for the sub-connected region extending outside C, q is the additional penalty value for unfilled pixels in C, and r is the penalty value for the overlapping part of the sub-connected regions. At most one sub-connected region is placed at each position i,j.

[0058] (2-7) Template matching with the connected parts of the fruit tree forest is achieved through convolution, and the matching results are stored in the fitness map. The fitness map of each template is calculated according to the "big first" principle. The higher the value of the pixel in the fitness map, the more suitable the crown center position is. The local maximum value in the fitness map is considered as the potential crown center. They are confirmed as the crown center after satisfying specific constraints.

[0059] The extraction of fruit tree density in the above-mentioned forest-grain intercropping area involves calculating the fruit tree density in step (3), and the specific steps are as follows:

[0060] (3-1) Reproject the plot results and template filling results to the same projection coordinate system, and then calculate the number of fruit trees inside each plot element in the plot vector file;

[0061] (3-2) Planting a certain number of fruit trees per mu can represent its planting density. Multiplying the area of ​​the plot element by 0.0015 will give the conversion between square meters and areles. Then, by calculating the number of fruit tree center points inside each plot element and dividing by the element area (mu), the fruit tree planting density can be obtained. The result is stored in the attribute table of the plot vector file.

[0062] Compared with existing technologies, the method for extracting fruit tree density in forest-grain intercropping areas provided by this invention achieves at least the following beneficial effects:

[0063] (1) The method for extracting fruit tree density in forest-grain intercropping areas provided by the present invention takes the plot as the smallest basic unit. The plot refers to a complete and closed farmland with the same owner. It is the smallest unit in agricultural production and management. Taking the plot as the smallest basic unit lays a good foundation for subsequent fruit tree density calculation and is more in line with the management of farmland in actual production and planting.

[0064] (2) The method for extracting fruit tree density in forest-grain intercropping areas provided by the present invention is more suitable for fruit tree density surveys in large-scale forest-grain intercropping areas. The remote sensing image on which the present invention is based is a pre-processed ultra-high resolution GF-2 image. Compared with the method of using radar data or ultra-high resolution images obtained by UAVs, the present invention can reduce manpower and material costs to a certain extent by improving image resolution through data fusion.

[0065] (3) The method for extracting fruit tree density in forest-grain intercropping areas provided by the present invention uses a template filling algorithm to obtain the crown center of the fruit tree and completes the task of high-precision forest stand tree count by generating templates. It is applicable to the planting of fruit trees in forest-grain intercropping areas dominated by forests.

[0066] (4) The method for extracting fruit tree density in forest-grain intercropping areas provided by this invention can provide support for scientifically adjusting the planting structure of vegetation in forest-grain intercropping areas, improve the economic benefits of forest-grain intercropping farmland, resolve the conflict between forest and grain land, and guide farmers to scientifically optimize the planting structure. Attached Figure Description

[0067] To more clearly illustrate the technical solutions in this invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0068] Figure 1 This is a flowchart illustrating the method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing technology, as described in this invention.

[0069] Figure 2 This is a model structure diagram of the irregular time series model used in this invention;

[0070] Figure 3 This is a comparison diagram of GF-2 data preprocessing before and after in this invention;

[0071] Figure 4 This is a model structure diagram of the D-LinkNet semantic segmentation model used in this invention;

[0072] Figure 5 This is a diagram showing the internal structure of module B of the D-LinkNet semantic segmentation model used in this invention.

[0073] Figure 6 This is a flowchart of the template filling algorithm used in this invention;

[0074] Figure 7 This is a semantic segmentation result diagram of fruit trees according to an embodiment of the present invention;

[0075] Figure 8 This is a template filling result diagram of an embodiment of the present invention;

[0076] Figure 9 This is the base map for land parcel attribute statistics in this embodiment of the invention. Detailed Implementation

[0077] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0078] The following is combined Figure 1 The present invention describes the process of a remote sensing-based method for extracting fruit tree density in a forest-grain intercropping area, as provided by the present invention, specifically including the following:

[0079] 1. Preprocessing GF-2 data to obtain ultra-high spatial resolution images includes six steps: radiometric calibration, atmospheric correction, orthorectification, image registration, image fusion, and image mosaicking.

[0080] 2. Based on the visual features of farmland plots, a farmland zoning and hierarchical system is established, dividing farmland into several relatively consistent subtypes. Based on the visual features presented by each type of plot, a suitable instance segmentation model is selected and improved to extract each type of farmland in a hierarchical manner. The subtypes of each farmland are then merged to obtain a complete farmland information extraction result map.

[0081] 3. Design time series features suitable for major food crops and major fruit trees, extract different vegetation phenological change features from long-term Sentinel-2 images, classify crops based on irregular time series models, and design discrimination conditions based on the classification results for the presence of only food crops, the presence of both food crops and fruit trees, and the presence of only fruit trees, so as to divide farmland into ordinary farmland and forest-grain intercropping plots.

[0082] 4. Using a semantic segmentation model based on convolutional neural networks, fruit trees in the preprocessed GF-2 image at 1m resolution are extracted to obtain the semantic segmentation results of the fruit trees. Based on the template filling algorithm, different templates are formed according to the spectral characteristics of the canopy of fruit trees of different sizes. Then, these generated templates are filled into the semantic segmentation results of the fruit trees to obtain the center position of the canopy of each fruit tree.

[0083] 5. The planting density of fruit trees is obtained by dividing the number of fruit tree center points within each plot element by the element area (mu).

[0084] The steps described above are explained in detail below.

[0085] 1. Preprocessing of GF-2 imagery: GF-2 data was obtained from the China Resources Satellite Application Center (http: / / 36.112.130.153:7777 / DSSPlatform / productSearch.html). The GF-2 satellite is my country's first domestically developed sub-meter-resolution civilian optical satellite with a tropical week of 69 days. It carries two PMS (panchromatic / multispectral) instruments, and the spectral range it can acquire is shown in Table 2.

[0086] Table 2 GF-2 Satellite Parameters

[0087]

[0088] Preprocessing of GF-2 data includes six steps: radiometric calibration, atmospheric correction, orthorectification, image registration, image fusion, and image stitching.

[0089] (1) Radiometric calibration is the process of converting the data values ​​(DN) recorded by the sensor into apparent reflectance, thereby eliminating the error of the sensor itself. It requires the calculation of the radiance value based on the absolute radiometric calibration coefficient, and the formula can be expressed as:

[0090] L λ =Gain·DN+Offset

[0091] Among them, L λ To obtain the converted standard radiance, the payload observation (DN) and absolute calibration coefficient offset (Offset) values ​​of the GF-2 satellite were obtained from the China Center for Resources Satellite Data and Application. The sensor parameter DN and Offset values ​​used in this invention are shown in Table 1.

[0092] Table 3 Absolute Radiation Calibration Factors

[0093]

[0094]

[0095] (2) Atmospheric correction refers to converting radiance or surface reflectance into actual surface reflectance, which can effectively eliminate the influence of atmospheric and light factors on the reflection of ground objects. This invention uses the FLAASH model for atmospheric correction. The model algorithm is based on the pixel-level correction of the MODTRAN 4+ radiative transfer model, which corrects the cascading effects caused by diffuse reflection. It includes a classification map of cirrus clouds and opaque clouds and can adjust the spectral smoothing caused by human suppression.

[0096] (3) Orthorectification, which is the process of generating a planar orthorectified image to eliminate the influence of terrain or the distortion caused by camera orientation. In addition to the conventional geometric correction, the distortion caused by terrain undulation is also corrected according to the DEM. This invention uses the orthorectification tool RPC Orthorectification Workflow and the digital elevation model adopts GMTED2010DEM. At the same time, in order to perform image fusion later, the panchromatic image of GF-2 is resampled to 1m resolution and the multispectral image is resampled to 4m resolution.

[0097] (4) Image registration: A reference image is used to calibrate another image in the same area. Since the panchromatic image and multispectral image of GF-2 do not correspond completely, the multispectral image is selected as the reference image and the panchromatic image is registered. The output resolution is 1m.

[0098] (5) Image fusion: The low-resolution multispectral image is fused with the high-resolution panchromatic image to obtain an image with both high spatial resolution and multispectral features; After orthorectification, the GF-2 image has a panchromatic band resolution of 1m and four multispectral band resolutions of 4m; This invention uses the NNDiffuse Pan Sharpening fusion method to fuse multispectral and panchromatic images.

[0099] (6) Image mosaicking: stitching together two or more remote sensing images with geographic coordinates to obtain a larger and more complete remote sensing image; before mosaicking the GF-2 image, the present invention needs to perform image registration, select the Sentinel-2 remote sensing image as the reference image, and the GF-2 image as the image to be registered for registration, and then use the Seamless Mosaic tool to mosaic the GF-2 image.

[0100] 2. Extraction of Farmland Plots: Extracting farmland plots requires "geographic zoning" to simulate the human eye's recognition process of farmland targets. This process gradually microscopicates complex macroscopic spatial perception problems, improving the interpretation of remote sensing images. First, with the aid of road networks, water systems, and topographic lines, the image is divided into naturally aggregated sub-regions using a masking method, ensuring relatively uniform land cover types within each sub-region. This step cuts the image into discrete geographic blocks, allowing the classification process to be confined to the corresponding sub-regions, thus avoiding the propagation of classification errors. Furthermore, dividing the image into blocks allows for parallel computation, improving extraction efficiency. Second, under the constraints of each geographic block, the edges of target features within each region need to be extracted. For farmland with blurred edges, texture-based object boundaries are used. Based on slope calculations using topographic data for zoning, the differences in geometric and textural features of farmland in different regions are expressed.

[0101] Farmland plots are extracted using an instance segmentation model, which combines object detection and semantic segmentation. Object detection is performed in the image, and then each pixel is labeled (semantic segmentation), allowing for the extraction of individual instances at a time. The instance segmentation algorithm used in this invention is the Hybrid Task Cascade (HTC) model, with a skeleton of ResNeXt-101-64x4d and Deformable ConvNets v2 (DCN). The algorithm output is post-processed using an overlapping sliding window strategy, and Soft-NMS is used to improve model performance.

[0102] This invention employs an overlapping sliding window strategy, which involves cropping numerous pixel blocks from the original image from top to bottom and left to right for instance segmentation. Each pixel block is divided into a target region and an ignored region. The left and top edges of the ignored region are two pixels from the grid boundary. When the coordinates (x, y) of the top left corner of the detection box fall within the target region, the detection result is retained. When the sliding window is on the boundary of the original image, ignored regions on the edges of the grid in the same direction are merged into the target region. The overlapping sliding window strategy should also satisfy the following rules: 1) The sliding window should be as large as possible to ensure that each farmland object appears completely within at least one sliding window. 2) When the coordinates (x, y) of the top left corner of the sliding window fall within the target region, the detection result should also be retained. 3) When the sliding window is on the boundary of the original image, edge regions on the same direction as the sliding window should be merged into the target region.

[0103] 3. Extraction of Agroforestry Plots: This invention utilizes time-series models based on the features of Sentinel-2 time-series images to distinguish between ordinary farmland plots and agroforestry plots, thus determining the plots for subsequent fruit tree density extraction. First, time-series features suitable for major food crops and major fruit trees are designed to extract different vegetation phenological change features from long-term Sentinel-2 images. This invention selects three vegetation time-series features to distinguish three plot types: Normalized Difference Vegetation Index (NDVI), Normalized Red Edge Index (NDRE1), and Red Edge Normalized Vegetation Index (NDVIre1). The calculation formula is as follows:

[0104]

[0105] Then, crop classification was performed based on an irregular time series model. Based on the classification results, discrimination conditions were designed for fields containing only food crops, fields containing both food crops and fruit trees, and fields containing only fruit trees, thus dividing farmland into ordinary farmland and intercropped forest-food plots. The structure of the irregular time series model is as follows: Figure 2 As shown, two-dimensional convolutional neural network (2DCNN) models are typically used to extract spatial and spectral features. The input is generally a three-dimensional matrix containing R, G, and B channels. Since time-series features are one-dimensional data, the main task of using 2DCNN to process time-series data is to convert the one-dimensional data into a two-dimensional matrix. In this invention, the 2DCNN model automatically learns the features of crop growth patterns by stretching the VI time series into a "time-feature" matrix. The classic residual network ResNet-18 model of 2DCNN is used, and after pre-training on the ImageNet dataset, it is trained. The process mainly includes two steps. First, at each pixel location, the irregular VI time series is extracted from the original optical image and converted into a "time-feature" matrix to represent the crop growth pattern. In the "time-feature" matrix, the horizontal axis represents time, and the vertical axis represents the value of the vegetation index. Second, the "time-feature" matrix obtained in the first step is input into the three channels of ResNet-18 to obtain a three-dimensional matrix suitable for ResNet-18, thereby extracting the features of vegetation phenological changes.

[0106] The formula for stretching a VI time series into a "time-feature" matrix is ​​as follows:

[0107]

[0108] Among them, VI ori The original value of VI, VI r VI is the stretched value, and length is the length of the time series.

[0109] Finally, after obtaining pixel-level classification results based on temporal features, three discrimination conditions were used to distinguish between ordinary farmland and forest-grain intercropping plots: 1. Only grain crops exist (grain crop distribution area is greater than 80%); 2. Grain crops and fruit trees exist (grain crop distribution area and fruit tree distribution area are greater than 30%); 3. Only fruit trees exist (fruit tree distribution area is greater than 80%), and forest-grain intercropping plots were selected.

[0110] 4. Extraction of the crown center of fruit trees in agroforestry plots: First, a semantic segmentation model based on convolutional neural networks is used to extract fruit trees from the preprocessed GF-2 image at 1m resolution, obtaining the semantic segmentation results of the fruit trees. The semantic segmentation model used is the D-LinkNet model.

[0111] The structure of the D-LinkNet network is as follows: Figure 4 As shown, D-LinkNet uses LinkNet with a pre-encoder as its backbone network and adds dilated convolution layers in the central part. Its advantage lies in expanding the receptive field of view through these dilated convolution layers. LinkNet is an effective semantic segmentation neural network, possessing the advantages of skip connections, residual blocks, and an encoder-decoder architecture. D-LinkNet uses a ResNet trained on the ImageNet dataset as its encoder and adds a dilated-convolution layer with shortcuts in the central part, making the entire network more powerful, with a larger receptive field, and capable of fusing multi-scale information. Therefore, the D-LinkNet model has stronger feature extraction capabilities and can effectively classify heavily confused patches.

[0112] The entire network is divided into three modules: A, B, and C. Module A is the encoder part of the network, and module C is the decoder part. The configuration of these two parts is the same as that of the LinkNet network. Module B is unique to D-LinkNet, and its internal structure is as follows: Figure 5 As shown. The parallel structure of module B can effectively fuse multi-scale features. Dilated-convolution can effectively expand the receptive field. The corresponding receptive fields from top to bottom are 31, 15, 7, 3, and 1, respectively. Finally, the results of each branch are added together to obtain the fused features.

[0113] During model training, the D-LinkNet model training parameters were set as follows: the initial learning rate was set to 0.01, and the learning rate was automatically adjusted using the adaptive Adam optimizer; the sample size was 512 pixels × 512 pixels; the training lasted for 100 epochs (training automatically ended on the 72nd epoch); the batch size was set to 16; and the remaining parameters were left as default. Figure 7 The image shown is a semantic segmentation result of fruit trees in an example. The model segmentation results are quite detailed and can effectively extract the fruit trees in the study area.

[0114] Secondly, based on the template filling algorithm, different templates are formed according to the spectral characteristics of fruit tree canopies of different sizes. These generated templates are then filled into the fruit tree semantic segmentation results to obtain the center position of the canopy of each fruit tree.

[0115] The template filling algorithm flowchart is as follows: Figure 6 As shown. Calculating the density of fruit trees within a plot requires determining the accurate central location and number of fruit trees. Due to differences in growth status and location, fruit tree canopies exhibit diversity. Observing individual fruit tree canopies reveals that trees with larger canopies tend to have a higher infrared spectral response at their center and a lower infrared response at the periphery. Conversely, fruit trees with smaller canopies generally have a lower infrared response. Therefore, different templates can be created based on the spectral characteristics of fruit tree canopies of different sizes. These generated templates can then be used to fill fruit tree groves, thus transforming the problem of dividing fruit tree groves into a filling problem.

[0116] The study area contains numerous individual fruit trees of varying sizes, allowing for the statistical extraction of general characteristics and the generation of reasonable templates. First, a labeled area was randomly selected, and all individual fruit trees were labeled to reduce subjective bias. Second, a series of square empty templates of odd lengths were created to ensure grid centering, and each fruit tree was assigned to the template with the closest area. After aligning the center of the fruit tree with the center of the template, the spectral average of the connected components within each template was calculated to generate the original template series. The original template series exhibited typical characteristics of the canopy of individual fruit trees. However, due to the significant shadows in the fruit tree images, the spectral characteristics of the fruit tree forest and individual fruit trees were not entirely identical. Therefore, a mathematical model of the typical characteristics of the original templates needed to be constructed, and preprocessing was performed to obtain better performance. Since the near-infrared response is stronger at the center of the fruit tree, the spectral characteristics of the template pixels were calculated using the following formula:

[0117] R=C

[0118]

[0119] Where R represents the NIR2 band (908nm), B represents the Red band (658.8nm), and G represents the Green band (546.2nm). C is a constant, BGmin and BG max These are the minimum and maximum values ​​for bands B and G, respectively. n represents the template length, and distance_to_cent represents the distance to the template center, in pixels.

[0120] Given a 4-connected region C of arbitrary shape on the target raster image, the pixel value in C is denoted as x. i,j Where i, j are the coordinates on the digital grid. Several square sub-connected regions S with varying side lengths are also given. n , where n is the index of the sub-connected region. Let y denote the pixel value in each sub-connected region. i,j , where i, j are the coordinates of the bounding box grid of the sub-connected components. Sub-connected components S are filled into C, and each sub-connected component can be reused an arbitrary number of times. The ultimate goal is to minimize the global value G on the digital grid.

[0121] G = sum(|x i,j -y i,j |)+p+q+r

[0122] Where p is the penalty value for the sub-connected region extending beyond C, q is the additional penalty value for unfilled pixels in C, and r is the penalty value for the overlapping portion of the sub-connected regions. Also, at most one sub-connected region is placed at each position i,j.

[0123] Template matching with the connected parts of the orchard is achieved through convolution, and the matching results are stored in a fitness map. The fitness map of each template is calculated according to the "big-first" principle; pixels with higher values ​​in the fitness map indicate a more suitable crown center location. Local maxima in the fitness map are considered potential crown centers, and they are confirmed as crown centers after satisfying specific constraints.

[0124] 5. Extraction of planting density of fruit trees in forest-grain intercropping plots: The planting density of fruit trees is obtained by dividing the number of fruit tree center points within each plot element by the element area (mu).

[0125] The results of land parcel, semantic segmentation, and template filling in this embodiment of the invention are shown in the figure below. Figure 8 As shown in the figure. In this invention, the fruit tree canopies have greater overlap and individual tree canopies are smaller. Therefore, when generating the template, adjustments were made to 3*3, 5*5, 7*7, and 9*9 pixels to more accurately identify the center point. Since the semantic segmentation results inevitably contain shadows of the fruit tree canopies, a threshold of 30 was set as the boundary between shadows and canopies based on the shadow reflectance of the fruit trees in the figure. Reflectance below 30 indicates shadows, while reflectance above 30 indicates potential canopy pixels. It can be seen that this invention effectively extracted the center point of the fruit tree canopy in the embodiment.

[0126] Based on this, the plot results and template filling results are reprojected to the same projection coordinate system, and then the number of fruit tree pixels within each feature in the plot vector file is calculated. The number of fruit trees planted per acre can represent the planting density. Multiplying the plot feature area by 0.0015 yields the conversion between square meters and areles. The fruit tree planting density is then obtained by dividing the number of fruit tree center points within each plot feature by the feature area (acres). The results are stored in the plot vector file attribute table, as shown in Table 4.

[0127] Table 4. Results of Land Parcel Attributes

[0128]

[0129] This invention has been implemented on a PC platform. Experiments have demonstrated that this solution can meet the correction requirements for erroneously merged regions in scale-set models. It offers excellent interactivity and fast response during the interactive process, achieving efficient and accurate interactive editing functions. It fully leverages the advantages of scale-set models and effectively solves the problem of erroneously merged regions within them. Furthermore, by combining advanced technologies such as minimum risk segmentation models and topology maps, it enables rapid location and correction of erroneously merged regions. This invention can be widely applied in geographic information systems, remote sensing image processing, and computer vision applications requiring image segmentation for region editing functions in scale-set models, providing strong support for technological advancements in these fields.

[0130] It should be noted that although some specific embodiments of the present invention have been described in detail by way of examples, those skilled in the art should understand that the above examples are for illustrative purposes only and are not intended to limit the scope of the present invention. Those skilled in the art should understand that modifications can be made to the above embodiments without departing from the scope and spirit of the present invention. The scope of the present invention is defined by the appended claims.

Claims

1. A method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing, characterized in that... It includes the following steps: Step 1: Preprocess the GF-2 data to obtain ultra-high spatial resolution images. The preprocessing process includes six steps: radiometric calibration, atmospheric correction, orthorectification, image registration, image fusion, and image mosaicking. Step 2: Based on the visual features of farmland plots, establish a farmland zoning and hierarchical system, divide farmland into several relatively consistent subtypes, and select and improve appropriate instance segmentation models based on the visual features presented by each type of plot to extract each type of farmland in a hierarchical manner, and merge the subtypes of each farmland to obtain a complete farmland information extraction result map. Step 3: Design time series features suitable for major food crops and major fruit trees, extract different vegetation phenological change features from long-term Sentinel-2 images, classify crops based on irregular time series models, and design discrimination conditions based on the classification results: only food crops, food crops and fruit trees, and only fruit trees, to divide farmland into ordinary farmland and forest-grain intercropping plots. Step 4: Using a semantic segmentation model based on convolutional neural networks, fruit trees in the preprocessed GF-2 image at 1m resolution are extracted to obtain the semantic segmentation results of the fruit trees. Based on the template filling algorithm, different templates are formed according to the spectral characteristics of the canopies of fruit trees of different sizes. Then, these generated templates are filled into the semantic segmentation results of the fruit trees to obtain the center position of the canopy of each fruit tree. Step 5: Calculate the planting density of fruit trees by dividing the number of fruit tree center points within each plot element by the element area.

2. The method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing according to claim 1, characterized in that, Methods for preprocessing GF-2 data include: (1) The aforementioned radiometric calibration refers to converting the data values ​​(DN) recorded by the sensor into apparent reflectance, thereby eliminating the sensor's own errors. It requires calculating the radiance value based on the absolute radiometric calibration coefficient, where L... λ For the converted standard radiance, the Offset absolute calibration coefficient offset can be expressed by the formula: L λ =Gain·DN+Offset (2) The atmospheric correction refers to converting radiance or surface reflectance into actual surface reflectance and using the FLAASH model for atmospheric correction. The FLAASH model algorithm is based on the pixel-level correction of the MODTRAN 4+ radiative transfer model. It corrects the cascading effects caused by diffuse reflection, including the classification map of cirrus clouds and opaque clouds, and can adjust the spectral smoothing caused by human suppression. (3) The orthorectification mentioned above refers to the process of eliminating the influence of terrain or the deformation caused by camera orientation to generate a planar orthorectified image. In addition to performing conventional geometric correction, the distortion caused by terrain undulation is also corrected according to the DEM. The orthorectification tool RPC Orthorectification Workflow is used. The digital elevation model adopts GMTED2010DEM. At the same time, in order to perform image fusion later, the panchromatic image of GF-2 is resampled to 1m resolution and the multispectral image is resampled to 4m resolution. (4) The image registration refers to calibrating one image (reference image) in the same area with another image. Since the panchromatic image and multispectral image of GF-2 do not correspond completely, the multispectral image is selected as the reference image and the panchromatic image is registered. The output resolution is 1m. (5) The image fusion refers to fusing a low-resolution multispectral image with a high-resolution panchromatic image to obtain an image that has both high spatial resolution and multispectral features. After orthorectification, the GF-2 image has a panchromatic band resolution of 1m and four multispectral band resolutions of 4m. The NNDiffuse Pan Sharpening fusion method is used to fuse the multispectral and panchromatic images. (6) The image mosaicking refers to stitching together two or more remote sensing images with geographic coordinates to obtain a larger and more complete remote sensing image. Before mosaicking the GF-2 image, image registration is required. The Sentinel-2 remote sensing image is selected as the reference image, and the GF-2 image is used as the image to be registered. Then, the GF-2 image is mosaicked using a seamless mosaicking tool.

3. The method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing according to claim 1, characterized in that, In step 2, the farmland zoning and layering system refers to dividing the image into naturally aggregated sub-regions with the assistance of road networks / water systems / topographic lines through masking, ensuring that the land cover types within each sub-region are relatively uniform. Under the constraints of each geographic block, the edges of target land cover are extracted within each region. For farmland with blurred edges, the object boundary is obtained based on texture. Based on the slope calculated from the topographic data for zoning, the system expresses the differences in geometric and textural features of farmland in different regions.

4. The method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing according to claim 1, characterized in that, In step 2, the instance segmentation model is a combination of object detection and semantic segmentation. The object is detected in the image (object detection), and then each pixel is labeled (semantic segmentation). A single instance can be extracted at a time. The instance segmentation algorithm used is the Hybrid Task Cascade (HTC) model, with the skeleton being ResNeXt-101-64x4d and DeformableConvNets v2 (DCN). The algorithm output is post-processed using an overlapping sliding window strategy, and Soft-NMS is used to improve model performance.

5. The method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing according to claim 4, characterized in that, In the post-processing of the algorithm output using the overlapping sliding window strategy, the overlapping sliding window strategy refers to cropping many pixel blocks from the original image from top to bottom and from left to right for instance segmentation. Each pixel block is divided into a target region and an ignored region. The left and top edges of the ignored region are two pixels away from the grid boundary. When the coordinates (x, y) of the upper left corner of the detection box fall within the target region, the detection result is retained. When the sliding window is on the boundary of the original image, the ignored region on the edge in the same direction as the grid is merged into the target region. At the same time, the overlapping sliding window strategy should meet the following rules: 1) The sliding window should be as large as possible to ensure that each farmland unit appears completely in at least one sliding window; 2) When the coordinates (x, y) of the upper left corner of the sliding window fall within the target region, the detection result should also be retained; 3) When the sliding window is on the boundary of the original image, the edge region in the same direction as the sliding window should be merged into the target region.

6. The method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing according to claim 1, characterized in that, This method utilizes a time-series model to extract vegetation phenological information. The time-series model is a two-dimensional convolutional neural network (2DCNN) model, typically used to extract spatial and spectral features. The input is generally a three-dimensional matrix containing R, G, and B channels. Since time-series features are one-dimensional data, the main task of using 2DCNN to process time-series data is to convert the one-dimensional data into a two-dimensional matrix. The 2DCNN model automatically learns the features of crop growth patterns by stretching the VI time series into a "time-feature" matrix. The method uses the classic ResNet-18 residual network model of 2DCNN, and after pre-training on the ImageNet dataset, it undergoes training. The process includes two steps: First, for each pixel location, the irregular vegetation index (VI) time series is extracted from the original optical image and converted into a "time-feature" matrix to characterize the crop growth pattern. In the "time-feature" matrix, the horizontal axis represents time and the vertical axis represents the value of the vegetation index. Secondly, the "time-feature" matrix obtained in the first step is input into the three channels of ResNet-18 to obtain a three-dimensional matrix suitable for ResNet-18, thereby extracting vegetation phenological change features. The formula for stretching a VI time series into a "time-feature" matrix is ​​as follows: Among them, VI ori The original value of VI, VI r VI is the stretched value, and length is the length of the time series.

7. The method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing as described in claim 6, characterized in that, A method for characterizing crop growth patterns using vegetation index time series, wherein the vegetation indices are the Normalized Difference Vegetation Index (NDVI), the Normalized Red Edge Index (NDRE1), and the Red Edge Normalized Vegetation Index (NDVIre1), and the calculation formulas are as follows:

8. The method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing as described in claim 1, characterized in that, The semantic segmentation model is based on a convolutional neural network. The semantic segmentation model is the D-LinkNet model, which uses Linknet with a pre-encoder as its backbone network and has additional dilated convolutional layers in the central part.

9. The method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing as described in claim 1, characterized in that, A method for obtaining the center of fruit tree canopies using a template filling algorithm is described. This algorithm involves statistically identifying general characteristics of fruit trees and generating appropriate templates. First, a labeled area is randomly selected, and all individual fruit trees are labeled to reduce subjective bias. Second, a series of square empty templates of odd lengths are created to ensure grid centering, and each fruit tree is assigned to the template with the closest area. After aligning the center of the fruit tree with the center of the template, the spectral average of the connected components within each template is calculated to generate the original template. The original template series exhibits the typical characteristics of individual fruit tree canopies. Because fruit tree images contain significant shadows, the spectral characteristics of fruit tree groves and individual fruit trees are not entirely identical. Therefore, a mathematical model needs to be constructed based on the typical characteristics of the original templates, and preprocessing is used to obtain better performance. Since the near-infrared response of the fruit tree center is stronger, the spectral characteristics of the template pixels are calculated using the following formula: Where R represents the NIR2 band (908nm), B represents the Red band (658.8nm), G represents the Green band (546.2nm), C is a constant, and BG min and BG max These are the minimum and maximum values ​​for bands B and G, respectively. n represents the template length, and distance_to_cent represents the distance to the center of the template, in pixels. Given a 4-connected region C of arbitrary shape on the target raster image, the pixel value in C is denoted as x. i,j Where i, j are the coordinates on the digital grid, and several square sub-connected regions Sn with different side lengths are given, where n is the index of the sub-connected region, and the pixel value in each sub-connected region is denoted as y. i,j Fill the sub-connected region S into C, each sub-connected region can be reused an arbitrary number of times, and calculate the global value G on the digital grid, as follows: G=sum(|x i,j -y i,j |)+p+q+r Where p is the penalty value for the sub-connected region extending outside C, q is the additional penalty value for unfilled pixels in C, and r is the penalty value for the overlapping part of the sub-connected regions. At most one sub-connected region is placed at each position i,j. The ultimate goal of the calculation is to minimize G. Template matching with the connected parts of the fruit tree forest is achieved through convolution, and the matching results are stored in the fitness map. The fitness map of each template is calculated according to the "big first" principle. The higher the value of the pixel in the fitness map, the more suitable the crown center position is. The local maximum value in the fitness map is considered as the potential crown center, and they are confirmed as the crown center after satisfying specific constraints.

10. The method for extracting fruit tree density in forest-grain intercropping areas based on remote sensing as described in claim 1, characterized in that, The planting density of fruit trees is obtained by dividing the number of fruit tree center points within each plot element by the element area (mu). The fruit tree coverage and the number of trees planted per mu are calculated, and the results are stored in the plot vector file attribute table. First, the plot results and template filling results are reprojected to the same projection coordinate system. Then, the number of fruit tree pixels within each element in the plot vector file is calculated. Second, the number of trees planted per mu represents the planting density. Multiplying the plot element area by 0.0015 yields the conversion between square meters and areles. Therefore, the fruit tree planting density is obtained by dividing the number of fruit tree center points within each plot element by the element area (mu), and the results are stored in the plot vector file attribute table.

Images