Crop three-dimensional structure parameter extraction method and system fusing LiDAR and optical remote sensing

Abstract

This invention provides a method and system for extracting three-dimensional structural parameters of crops by integrating LiDAR and optical remote sensing. It relates to the field of agricultural remote sensing technology. The method involves spatially overlaying the spatial height distribution information of each LiDAR measurement point in the original LiDAR remote sensing data set with the spectral reflectance intensity information of optical image pixels at corresponding spatial locations in the original optical remote sensing image data set to obtain a joint height-spectral observation vector at each spatial location. Based on the local variation parameter, the set of crop plant center point locations is identified. Using the center point as the growth origin, the boundary of the plant's spatially occupied area is searched. Finally, subsets of the canopy's upper and lower surface point clouds are extracted from the three-dimensional space enclosed by the boundary, and plant height and canopy diameter parameters are calculated. This invention fully utilizes the advantages of LiDAR and optical remote sensing, enabling accurate and efficient extraction of three-dimensional structural parameters of crops.

Description

Technical Field

[0001] This invention relates to the field of agricultural remote sensing technology, and more specifically, to a method and system for extracting three-dimensional structural parameters of crops by integrating LiDAR and optical remote sensing. Background Technology

[0002] Traditional methods for obtaining three-dimensional structural parameters of crops mainly rely on manual field measurements. These methods not only consume a lot of manpower, material resources, and time, but also have a limited measurement range, making it difficult to conduct comprehensive and efficient monitoring of large areas of farmland.

[0003] With the development of remote sensing technology, LiDAR (Light Detection and Ranging) and optical remote sensing technologies are increasingly being applied in agriculture. LiDAR remote sensing can acquire spatial height distribution information of target objects and has high vertical resolution, but its ability to acquire spectral information of target objects is limited. Optical remote sensing can acquire spectral reflectance intensity information of target objects, reflecting their spectral characteristics, but it is insufficient in acquiring spatial height information. Currently, although some studies have attempted to fuse LiDAR and optical remote sensing data, most of these are simply data overlays without fully considering the inherent connections and synergistic effects between the two types of data. This makes it difficult to accurately and efficiently extract the three-dimensional structural parameters of crops, and thus fails to meet the needs of modern agriculture for precise monitoring. Summary of the Invention

[0004] In view of the aforementioned problems, and in conjunction with the first aspect of the present invention, embodiments of the present invention provide a method for extracting three-dimensional structural parameters of crops by integrating LiDAR and optical remote sensing, the method comprising: Acquire the raw LiDAR remote sensing measurement data set and the raw optical remote sensing image data set of the target farmland area; Spatial overlay analysis is performed on the spatial height distribution information of each LiDAR measurement point in the raw LiDAR remote sensing data set and the spectral reflectance intensity information of the optical image pixels at the corresponding spatial location in the raw optical remote sensing data set to obtain the joint observation vector of height spectrum at each spatial location. Based on the local variation parameter between the spatial height distribution information and the spectral reflectance intensity information in the joint height-spectral observation vector, the set of crop plant center point locations is identified in the target farmland area. Taking the center point of each crop plant as the growth origin, the height spectrum joint observation vector of adjacent spatial positions is searched layer by layer in the horizontal direction until the spatial height distribution information in the height spectrum joint observation vector at the search boundary decreases by more than a preset change threshold, thus obtaining the boundary of the plant space occupied by each crop plant. Extract the top and bottom canopy point cloud subsets of each crop plant from the three-dimensional space enclosed by the boundary of the area occupied by each plant. Calculate the plant height parameter of each crop plant based on the difference between the average spatial height value of the top and bottom canopy point cloud subsets. Calculate the canopy diameter parameter of each crop plant based on the maximum horizontal span of the boundary of the area occupied by the plant.

[0005] Furthermore, embodiments of the present invention also provide a system for extracting three-dimensional structural parameters of crops by integrating LiDAR and optical remote sensing, comprising: A processor; a machine-readable storage medium for storing machine-executable instructions of the processor; wherein the processor is configured to perform the above-described method for extracting three-dimensional structural parameters of crops by fusing LiDAR and optical remote sensing by executing the machine-executable instructions.

[0006] In another aspect, embodiments of the present invention also provide a computer program product, the computer program product including machine-executable instructions stored in a computer-readable storage medium, the processor of the crop three-dimensional structure parameter extraction system integrating LiDAR and optical remote sensing reads the machine-executable instructions from the computer-readable storage medium, the processor executes the machine-executable instructions, causing the crop three-dimensional structure parameter extraction system integrating LiDAR and optical remote sensing to perform the above-described crop three-dimensional structure parameter extraction method integrating LiDAR and optical remote sensing.

[0007] Based on the above, firstly, the spatial height distribution information of each LiDAR measurement point in the raw LiDAR remote sensing data set is spatially overlaid with the spectral reflectance intensity information of the corresponding optical image pixels in the raw optical remote sensing data set to obtain a joint height-spectral observation vector at each spatial location. This fully utilizes the advantages of both LiDAR and optical remote sensing, enabling a more comprehensive and accurate description of the target farmland area. Secondly, the set of crop plant center point locations is identified based on the local variation parameter between the spatial height distribution information and the spectral reflectance intensity information in the joint height-spectral observation vector. This improves the accuracy and reliability of plant center point identification and avoids potential misjudgments and omissions in traditional methods. Furthermore, taking the center point of each crop plant as the origin of growth, the joint observation vector of the height spectrum at adjacent spatial locations is searched layer by layer along the horizontal direction until the decrease in the spatial height distribution information in the joint observation vector of the height spectrum at the search boundary exceeds the preset change threshold. This yields the boundary of the spatial area occupied by each crop plant, which can accurately determine the spatial range of the plant. Finally, the point cloud subsets of the upper and lower surfaces of the canopy are extracted from the three-dimensional space enclosed by the boundary of the spatial area occupied by each plant, and the plant height and canopy diameter parameters are calculated accordingly. This allows for the accurate and efficient acquisition of the three-dimensional structural parameters of crops, which helps to improve agricultural production efficiency and quality. Attached Figure Description

[0008] Figure 1 This is a schematic diagram of the execution flow of the method for extracting three-dimensional structural parameters of crops by integrating LiDAR and optical remote sensing, provided in an embodiment of the present invention.

[0009] Figure 2 This is a schematic diagram of exemplary hardware and software components of the crop three-dimensional structure parameter extraction system that integrates LiDAR and optical remote sensing provided in an embodiment of the present invention. Detailed Implementation

[0010] Figure 1 This is a flowchart illustrating a method for extracting three-dimensional structural parameters of crops by integrating LiDAR and optical remote sensing, provided in one embodiment of the present invention. A detailed description follows.

[0011] Step S110: Obtain the set of raw LiDAR remote sensing measurement data and the set of raw optical remote sensing image data for the target farmland area.

[0012] In this embodiment, for a target farmland area requiring extraction of three-dimensional structural parameters of crops, both a raw LiDAR remote sensing measurement data set and a raw optical remote sensing image data set for the target farmland area are acquired simultaneously. The raw LiDAR remote sensing measurement data set is acquired during Earth observation using an airborne LiDAR system. This set contains a large number of discrete LiDAR measurement points, each recording its coordinate information in three-dimensional space. Specifically, it includes the horizontal coordinate value X_i along the horizontal axis, the horizontal coordinate value Y_i along the vertical axis, and the vertical height value Z_i relative to the geoid. The subscript i indicates the sequence number of the i-th LiDAR measurement point in the raw LiDAR remote sensing measurement data set. The spatial distribution density of these LiDAR measurement points is determined by the pulse transmission frequency of the airborne LiDAR system and the flight altitude of the flight platform. The optical remote sensing raw image dataset is acquired by a multispectral sensor mounted on the same or different flight platforms. The aforementioned optical remote sensing raw image dataset contains digital images covering the target farmland area. The digital images are composed of regularly arranged optical image pixels. Each optical image pixel records the reflection intensity value R_j in the red band, the reflection intensity value G_j in the green band, and the reflection intensity value NIR_j in the near-infrared band. The subscript j indicates the sequence number of the j-th optical image pixel in the aforementioned optical remote sensing raw image dataset.

[0013] To ensure the accuracy of subsequent spatial overlay analysis, when acquiring the aforementioned raw LiDAR remote sensing measurement data set and the aforementioned raw optical remote sensing image data set, the platform position and attitude parameters at the acquisition time of each LiDAR measurement point and each optical image pixel are synchronously recorded by the GPS receiver and inertial measurement unit. This allows the aforementioned raw LiDAR remote sensing measurement data set and the aforementioned raw optical remote sensing image data set to be unified into the same planar projection coordinate system through coordinate transformation.

[0014] Step S120: Spatial overlay analysis is performed on the spatial height distribution information of each LiDAR measurement point in the raw LiDAR remote sensing data set and the spectral reflectance intensity information of the optical image pixels at the corresponding spatial location in the raw optical remote sensing image data set to obtain the joint observation vector of height spectrum at each spatial location.

[0015] Step S121: Extract the horizontal x-axis coordinate, horizontal y-axis coordinate, and vertical height value of each LiDAR measurement point from the raw LiDAR remote sensing measurement data set. Combine the horizontal x-axis coordinate and the horizontal y-axis coordinate to generate a planar position identifier for each LiDAR measurement point. Use the vertical height value as the spatial height distribution information corresponding to the planar position identifier.

[0016] In this embodiment, information is first extracted from each LiDAR measurement point in the aforementioned raw LiDAR remote sensing data set. For any LiDAR measurement point in the aforementioned raw LiDAR remote sensing data set, the horizontal x-axis coordinate value X_i, the horizontal y-axis coordinate value Y_i, and the vertical height value Z_i of the LiDAR measurement point are read from the data record of the LiDAR measurement point. The aforementioned horizontal x-axis coordinate value X_i and the aforementioned horizontal y-axis coordinate value Y_i represent the east-west and north-south positions of the aforementioned LiDAR measurement point in the universal transverse Mercator projection coordinate system, respectively. The horizontal axis coordinate value X_i and the horizontal axis coordinate value Y_i are combined into a unique planar location identifier ID_i by string concatenation. The generation rule for the planar location identifier ID_i is as follows: the horizontal axis coordinate value X_i is converted into a string and concatenated with a separator string "", and then the converted string of the horizontal axis coordinate value Y_i is concatenated to the end of the separator string "", forming an identifier string that can uniquely identify a specific planar location in the target farmland area. The vertical height value Z_i is used as the spatial height distribution information corresponding to the planar location identifier ID_i. The spatial height distribution information represents the absolute elevation of the ground surface or ground feature at the location of the LiDAR measurement point from the geoid.

[0017] Step S122: Extract the pixel row index, pixel column index, and red band reflection intensity, green band reflection intensity, and near-infrared band reflection intensity values ​​of each optical image pixel from the optical remote sensing raw image data set, and combine the pixel row index and the pixel column index to generate a pixel plane position identifier for each optical image pixel.

[0018] In this embodiment, information is then extracted from each optical image pixel in the aforementioned optical remote sensing raw image dataset. The digital images in the aforementioned optical remote sensing raw image dataset are arranged according to row and column rules, and each optical image pixel has a unique row and column number in the aforementioned digital image. For any optical image pixel in the aforementioned optical remote sensing raw image dataset, the pixel row index value R_idx_k, pixel column index value C_idx_k, red band reflectance intensity value R_k, green band reflectance intensity value G_k, and near-infrared band reflectance intensity value NIR_k of the aforementioned optical image pixel are read from the data record of the aforementioned optical image pixel, where the subscript k represents the sequence number of the k-th optical image pixel in the aforementioned optical remote sensing raw image dataset. The aforementioned red light band reflection intensity value R_k represents the reflected radiation energy value of the aforementioned optical image pixel in the red light band, the aforementioned green light band reflection intensity value G_k represents the reflected radiation energy value of the aforementioned optical image pixel in the green light band, and the aforementioned near-infrared band reflection intensity value NIR_k represents the reflected radiation energy value of the aforementioned optical image pixel in the near-infrared band. The aforementioned pixel row index value R_idx_k and the aforementioned pixel column index value C_idx_k are combined by string concatenation to form a unique pixel planar location identifier PIX_ID_k. The generation rule for the pixel planar location identifier PIX_ID_k is as follows: the aforementioned pixel row index value R_idx_k is converted into a string and concatenated with a delimiter string "", and then the converted string of the aforementioned pixel column index value C_idx_k is concatenated to the end of the delimiter string "", forming an identifier string that can uniquely identify a specific pixel location in the aforementioned digital image.

[0019] Step S123: Perform a string matching comparison between the planar position identifier of the LiDAR measurement point and the pixel planar position identifier of the optical image pixel. When the planar position identifier and the pixel planar position identifier are the same, it is determined that the LiDAR measurement point and the optical image pixel are in the same spatial position.

[0020] In this embodiment, a matching comparison operation is then performed on the planar location identifier ID_i of each LiDAR measurement point generated in step S121 and the pixel planar location identifier PIX_ID_k of each optical image pixel generated in step S122. Each LiDAR measurement point in the raw LiDAR remote sensing data set is traversed to obtain its planar location identifier ID_i. Each optical image pixel in the raw optical remote sensing data set is traversed to obtain its pixel planar location identifier PIX_ID_k. A character-by-character complete matching comparison is performed between the string of the planar location identifier ID_i and the string of the pixel planar location identifier PIX_ID_k. When every character in the string of the planar location identifier ID_i and the string of the pixel planar location identifier PIX_ID_k are identical and the lengths of the two strings are equal, it is determined that the planar location identifier ID_i and the pixel planar location identifier PIX_ID_k are the same. At this point, it is determined that the LiDAR measurement point and the optical image pixel correspond to the same spatial location in the target farmland area. If the above-mentioned planar position identifier ID_i string and the above-mentioned pixel planar position identifier PIX_ID_k string are different in any one character or the two strings are not of equal length, it is determined that the above-mentioned planar position identifier ID_i and the above-mentioned pixel planar position identifier PIX_ID_k are different. At this time, the above-mentioned LiDAR measurement point and the above-mentioned optical image pixel do not correspond to the same spatial position.

[0021] Step S124: Perform a vector concatenation operation on the vertical height value of the LiDAR measurement point located in the same spatial position and the red light band reflection intensity value, green light band reflection intensity value, and near-infrared band reflection intensity value of the optical image pixel to generate a five-dimensional height-spectral joint observation vector. The first dimension of the height-spectral joint observation vector stores the vertical height value, the second dimension stores the red light band reflection intensity value, the third dimension stores the green light band reflection intensity value, the fourth dimension stores the near-infrared band reflection intensity value, and the fifth dimension stores the ratio vegetation index calculated from the red light band reflection intensity value and the near-infrared band reflection intensity value.

[0022] In this embodiment, the information fusion operation is then performed on the LiDAR measurement point and optical image pixel determined to be in the same spatial position in step S123. The vertical height value Z_i of the LiDAR measurement point is obtained, and the red band reflectance intensity value R_k, green band reflectance intensity value G_k, and near-infrared band reflectance intensity value NIR_k of the optical image pixel in the same spatial position as the LiDAR measurement point are obtained. A five-dimensional vector data structure V is created as the height-spectral joint observation vector. The vertical height value Z_i is stored in the first dimension position V1 of the height-spectral joint observation vector V. The red band reflectance intensity value R_k is stored in the second dimension position V2 of the height-spectral joint observation vector V. The green band reflectance intensity value G_k is stored in the third dimension position V3 of the height-spectral joint observation vector V. The near-infrared band reflectance intensity value NIR_k is stored in the fourth dimension position V4 of the height-spectral joint observation vector V. Calculate the ratio vegetation index value NDVI_value = NIR_k / R_k. Store the calculated ratio vegetation index value NDVI_value in the fifth dimension position V5 of the above-mentioned height-spectral joint observation vector V.

[0023] Step S125: Traverse all spatial locations in the target farmland area, repeatedly perform the string complete matching comparison operation and the vector concatenation operation to generate a spatial distribution matrix of the joint height-spectral observation vector covering all spatial locations in the target farmland area. The row index of the spatial distribution matrix corresponds to the spatial location number in the horizontal axis coordinate direction, and the column index of the spatial distribution matrix corresponds to the spatial location number in the horizontal axis coordinate direction. Each cell of the spatial distribution matrix stores one joint height-spectral observation vector.

[0024] In this embodiment, steps S123 and S124 are then repeated for all spatial locations within the target farmland area. The spatial location range along the horizontal axis is determined based on the minimum and maximum coordinate values ​​X_min and X_max of the target farmland area along the horizontal axis. Similarly, the spatial location range along the vertical axis is determined based on the minimum and maximum coordinate values ​​Y_min and Y_max of the target farmland area along the horizontal axis. This spatial location range along the horizontal axis is then divided into M discrete spatial location numbers at equal intervals according to the spatial resolution delta_X of the LiDAR measurement points. Each spatial location number m corresponds to a horizontal axis coordinate value X_m, where m is a positive integer from 1 to M. The spatial location range along the vertical axis is then divided into N discrete spatial location numbers at equal intervals according to the spatial resolution delta_Y of the LiDAR measurement points. Each spatial location number n corresponds to a horizontal axis coordinate value Y_n, where n is a positive integer from 1 to N. A two-dimensional matrix structure of M rows and N columns is created as the spatial distribution matrix of the joint height-spectral observation vector. The cell in the m-th row and n-th column of the spatial distribution matrix corresponds to the planar position determined by the horizontal coordinate value X_m and the horizontal coordinate value Y_n. For each cell in the spatial distribution matrix, the horizontal coordinate value X_m and the horizontal coordinate value Y_n of the planar position corresponding to the cell are obtained, and a planar position identifier ID_mn is generated for the planar position according to the method in step S121. An optical image pixel matching the planar position identifier ID_mn is searched from the string complete match comparison results in step S123. If a matching optical image pixel is found and a corresponding LiDAR measurement point exists, the joint height-spectral observation vector V_mn for the planar position is generated according to the method in step S124, and the joint height-spectral observation vector V_mn is stored in the cell in the m-th row and n-th column of the spatial distribution matrix. If no matching LiDAR measurement point or optical image pixel exists for the planar position, the cell in the m-th row and n-th column of the spatial distribution matrix is ​​marked as empty. After the above traversal operation, a spatial distribution matrix of the joint height spectral observation vectors covering all spatial locations of the target farmland area is generated. The size of the spatial distribution matrix is ​​M rows multiplied by N columns.

[0025] Step S126: Using each cell in the spatial distribution matrix as the current analysis cell, extract the vertical height value on the first dimension of the joint height-spectral observation vector stored in the current analysis cell, and simultaneously extract the eight adjacent vertical height values ​​on the first dimension of the eight joint height-spectral observation vectors stored in the eight adjacent cells in the horizontal direction of the current analysis cell.

[0026] In this embodiment, the plant center point candidate location identification operation is then performed on each cell of the spatial distribution matrix generated in step S125. The cell in the m-th row and n-th column of the spatial distribution matrix is ​​taken as the current analysis cell CELL_mn. First, it is determined whether the current analysis cell CELL_mn stores a null value. If a null value is stored, the current analysis cell CELL_mn is skipped and the next cell is processed. If the current analysis cell CELL_mn stores a valid height-spectral joint observation vector V_mn, the vertical height value Z_mn is extracted from the first dimension position V1_mn of the height-spectral joint observation vector V_mn. Next, obtain the set of eight adjacent cells in the horizontal direction of the current analysis cell CELL_mn. These eight adjacent cells include the cell CELL_m-1_n-1 in the (m-1)th row and (n-1)th column, the cell CELL_m-1_n in the (m-1)th row and (n-1)th column, the cell CELL_m-1_n in the (m-1)th row and (n+1)th column, the cell CELL_m_n-1 in the (m)th row and (n+1)th column, the cell CELL_m+1_n-1 in the (m+1)th row and (n-1)th column, the cell CELL_m+1_n in the (m+1)th row and (n-1)th column, and the cell CELL_m+1_n in the (m+1)th row and (n+1)th column. For each of the eight adjacent cells, determine whether the cell stores a null value. If it does, skip the cell. If it stores a valid joint height-spectral observation vector, extract the adjacent vertical height values ​​of the cell from the first dimension of the stored joint height-spectral observation vector. After the above extraction operation, obtain the vertical height value Z_mn of the current analysis cell CELL_mn and a set of up to eight adjacent vertical height values ​​Z_adj.

[0027] Step S127: Calculate the absolute value of the height difference between the vertical height value of the current analysis cell and the arithmetic mean of the eight adjacent vertical height values, and use the absolute value of the height difference as the local height protrusion parameter of the current analysis cell.

[0028] In this embodiment, the vertical height value Z_mn extracted in step S126 and the set of adjacent vertical height values ​​Z_adj are then calculated. The count_adj of all valid adjacent vertical height values ​​in the set Z_adj is obtained. If the count_adj is zero, the local height convexity parameter H_rise_mn of the current analysis cell CELL_mn is set to zero and subsequent calculations are skipped. If the count_adj is greater than zero, Z_adj_mean = (Z_adj_1 + Z_adj_2 + ... + Z_adj_count_adj) / count_adj is calculated. Then, delta_Z = |Z_mn - Z_adj_mean| is calculated. Delta_Z is used as the local height convexity parameter H_rise_mn of the current analysis cell CELL_mn.

[0029] Step S128: Extract the calculated value of the ratio vegetation index on the fifth dimension of the joint height-spectral observation vector stored in the current analysis cell, and simultaneously extract the calculated values ​​of the eight adjacent ratio vegetation indices on the fifth dimension of the eight joint height-spectral observation vectors stored in the eight adjacent cells in the horizontal direction of the current analysis cell.

[0030] In this embodiment, the ratio vegetation index calculation value is then extracted from the current analysis cell CELL_mn. The ratio vegetation index calculation value NDVI_mn is extracted from the fifth dimension position V5_mn of the height-spectral joint observation vector V_mn stored in the current analysis cell CELL_mn. Next, a set of eight adjacent cells in the horizontal direction is obtained, which are the same as the eight adjacent cells described in step S126. For each of these eight adjacent cells, it is determined whether the cell stores a null value; if a null value is stored, the cell is skipped; if a valid height-spectral joint observation vector is stored, the adjacent ratio vegetation index calculation value of that cell is extracted from the fifth dimension position of the height-spectral joint observation vector stored in that cell. After the above extraction operation, the ratio vegetation index calculation value NDVI_mn of the current analysis cell CELL_mn and a set of up to eight adjacent ratio vegetation index calculation values ​​NDVI_adj are obtained.

[0031] Step S129: Calculate the vegetation index difference between the calculated value of the ratio vegetation index of the current analysis cell and the arithmetic mean of the calculated values ​​of the eight adjacent ratio vegetation indices, and use the vegetation index difference as the local vegetation vigor parameter of the current analysis cell.

[0032] In this embodiment, the ratio vegetation index calculated value NDVI_mn and the set of adjacent ratio vegetation index calculated values ​​NDVI_adj extracted in step S128 are then calculated. The number of all valid adjacent ratio vegetation index calculated values ​​in the set NDVI_adj is obtained as count_adj_ndvi. If count_adj_ndvi is zero, the local vegetation growth parameter V_growth_mn of the current analysis cell CELL_mn is set to zero and subsequent calculations are skipped. If count_adj_ndvi is greater than zero, NDVI_adj_mean is calculated as NDVI_adj_mean = (NDVI_adj_1 + NDVI_adj_2 + ... + NDVI_adj_count_adj_ndvi) / count_adj_ndvi. Then, delta_NDVI = |NDVI_mn - NDVI_adj_mean| is calculated. Delta_NDVI is used as the local vegetation growth parameter V_growth_mn for the current analysis cell CELL_mn.

[0033] Step S1210: Normalize the local height protrusion parameter and the local vegetation vigor parameter respectively. Multiply the normalized local height protrusion parameter and the normalized local vegetation vigor parameter to obtain the plant center point confidence score of the current analysis cell. When the plant center point confidence score exceeds the preset confidence threshold, mark the spatial position corresponding to the planar position identifier of the current analysis cell as a candidate position of crop plant center point.

[0034] In this embodiment, the H_rise_mn obtained in step S127 and the V_growth_mn obtained in step S129 are then normalized and confidence scores are calculated. First, all cells in the spatial distribution matrix are traversed, collecting all values ​​of H_rise_mn, and the maximum value H_rise_max and minimum value H_rise_min are identified. For the current analysis cell CELL_mn, H_rise_norm_mn = (H_rise_mn - H_rise_min) / (H_rise_max - H_rise_min). Next, all cells in the spatial distribution matrix are traversed, collecting all values ​​of V_growth_mn, and the maximum value V_growth_max and minimum value V_growth_min are identified. For the current analysis cell CELL_mn, V_growth_norm_mn = (V_growth_mn - V_growth_min) / (V_growth_max - V_growth_min). Calculate C_score_mn = H_rise_norm_mn × V_growth_norm_mn. Obtain a pre-set confidence threshold TH_conf. Compare C_score_mn with TH_conf. If C_score_mn > TH_conf, mark the spatial location corresponding to the planar location identifier ID_mn of the current analysis cell CELL_mn as a candidate location for the center point of a crop plant.

[0035] Step S1211: Perform local non-maximum suppression on all cells in the spatial distribution matrix that are marked as candidate locations of the crop plant center point. Compare the confidence scores of the plant center point of adjacent candidate locations. Retain the candidate location with the highest confidence score of the plant center point in the local area as the final location of the crop plant center point. Delete all other candidate locations in the local area that do not have the highest confidence score of the plant center point, and obtain the set of locations of the crop plant center point.

[0036] In this embodiment, local non-maximum suppression is then performed on all candidate crop plant center point locations marked in step S1210 to eliminate multiple candidate locations appearing around the same plant. All cells marked as candidate crop plant center point locations in the spatial distribution matrix are collected, and the planar location identifier and corresponding plant center point confidence score C_score for each candidate cell are obtained. For each candidate location in the candidate location set, a local neighborhood window is defined centered on that candidate location. The size of the local neighborhood window is a preset window radius parameter R_window, and the local neighborhood window covers all cells within a square area centered on the candidate location with a side length of twice R_window plus one. Within the local neighborhood window, all cells marked as candidate crop plant center point locations are extracted, and the plant center point confidence scores C_scores of these candidate locations are compared. From all candidate positions within the aforementioned local neighborhood window, select the candidate position with the highest plant center point confidence score C_score as the retained position within the local neighborhood window. Remove the other candidate positions from the crop plant center point candidate position set. Iterate through all candidate positions, repeatedly performing the operations of defining the local neighborhood window, comparing confidence scores, retaining the highest-scoring position, and deleting other positions, until the non-maximum suppression operation is completed in all local regions. After the aforementioned local non-maximum suppression operation, only the candidate position with the highest plant center point confidence score is retained in each local region as the final crop plant center point position. Collect all retained crop plant center point positions and organize them into a set according to their planar position identifiers, obtaining the crop plant center point position set P_center. Each element in the crop plant center point position set P_center corresponds to the center point position of a detected crop plant in the target farmland area.

[0037] Step S130: Based on the local variation parameter between the spatial height distribution information and the spectral reflectance intensity information in the height-spectral joint observation vector, identify the set of crop plant center point locations in the target farmland area.

[0038] In this embodiment, steps S126 to S1211 together realize the complete process of identifying the set of crop plant center point locations based on the local variation parameter between the spatial height distribution information and the spectral reflectance intensity information in the joint height-spectral observation vector. Steps S126 and S127 calculate the local height protrusion parameter as the local variation parameter of the spatial height distribution information based on the vertical height value. Steps S128 and S129 calculate the local vegetation vigor parameter as the local variation parameter of the spectral reflectance intensity information based on the ratio vegetation index. Step S1210 fuses the two local variation parameters to obtain the plant center point confidence score. Step S1211 obtains the final set of crop plant center point locations through local nonmaximum suppression.

[0039] Step S140: Using the center point of each crop plant as the growth origin, search the joint height spectrum observation vector of adjacent spatial locations layer by layer in the horizontal direction until the spatial height distribution information in the joint height spectrum observation vector at the search boundary decreases by more than a preset change threshold, thereby obtaining the boundary of the plant space occupied by each crop plant.

[0040] In this embodiment, each crop plant center point in the set P_center obtained in step S1211 above is then used as the growth origin, and a layer-by-layer outward expansion search operation is performed to determine the spatial boundary occupied by each crop plant.

[0041] Step S141: Using the center point position of each crop plant in the set of crop plant center point positions as the current growth origin, obtain the initial height spectral joint observation vector stored in the corresponding cell of the spatial distribution matrix of the current growth origin, and extract the initial vertical height value in the initial height spectral joint observation vector as the height reference value of the growth origin.

[0042] In this embodiment, the p-th crop plant center point in the set P_center is first processed and taken as the current growth origin O_p. The cell CELL_p corresponding to the planar position identifier of the current growth origin O_p is then searched in the spatial distribution matrix, and the height-spectral joint observation vector V_p stored in cell CELL_p is obtained. An initial vertical height value Z_p_init is extracted from the first dimension position V1_p of the height-spectral joint observation vector V_p, and this initial vertical height value Z_p_init is taken as the growth origin height reference value H_base_p.

[0043] Step S142: Using the cell of the current growth origin in the spatial distribution matrix as the starting cell of the first layer search, record the planar position identifier of the starting cell of the first layer search as the first layer cell identifier set of the plant occupied area.

[0044] In this embodiment, the cell CELL_p corresponding to the current growth origin O_p in the spatial distribution matrix is ​​then used as the starting cell for the first-layer search. A set of cell identifiers for the plant-occupied region S_region_p is created, and the planar position identifier of the starting cell CELL_p is added to the set of cell identifiers for the plant-occupied region S_region_p. Simultaneously, the starting cell CELL_p is recorded as the first-layer cell identifier set S_layer1_p for the plant-occupied region.

[0045] Step S143: Starting from the first search starting cell, expand outward by one layer to reach the second search cell set, which contains all cells that are horizontally adjacent to the first search starting cell and are not marked as plant-occupied areas.

[0046] In this embodiment, the outward expansion operation is then performed starting from the first-layer search starting cell CELL_p. All cells horizontally adjacent to the first-layer search starting cell CELL_p are obtained. These adjacent cells include the cells above, below, to the left, to the right, and along the four diagonal directions of the first-layer search starting cell CELL_p. Cells that have not yet been added to the plant-occupied area cell identifier set S_region_p are selected from these adjacent cells, and these selected cells form the second-layer search cell set S_candidates_layer2_p.

[0047] Step S144: Extract the vertical height value of the joint observation vector of the height spectrum stored in each cell of the second layer search cell set, and calculate the height reduction ratio parameter of the vertical height value of each cell to the height reference value of the growth origin.

[0048] In this embodiment, each cell in the second-layer search cell set S_candidates_layer2_p is then processed. For the q-th cell CELL_q in S_candidates_layer2_p, the joint height-spectral observation vector V_q stored in CELL_q is obtained, and the vertical height value Z_q is extracted from the first dimension position V1_q of V_q. R_drop_q = (H_base_p - Z_q) / H_base_p is then calculated.

[0049] Step S145: When the height decrease ratio parameter is less than the preset change threshold, add the planar position identifier of the cell to the second layer cell identifier set of the plant-occupied area, and continue to expand outward from the cell to perform the next layer search.

[0050] In this embodiment, the height drop ratio parameter R_drop_q calculated in step S144 is then compared with a preset change threshold TH_height. The preset change threshold TH_height is a pre-set percentage value used to determine whether the drop in vertical height has reached the boundary condition. When the height drop ratio parameter R_drop_q is less than the preset change threshold TH_height, it is determined that the cell CELL_q still belongs to the plant-occupied area of ​​the crop plant corresponding to the current growth origin O_p. The planar position identifier of the cell CELL_q is added to the cell identifier set S_region_p of the plant-occupied area, and the cell CELL_q is marked as a second-layer cell and added to the second-layer cell identifier set S_layer2_p of the plant-occupied area. Then, using the cell CELL_q as a new starting point, the search continues to expand outward by one layer, that is, the operations of steps S143 to S146 are repeated, and adjacent unmarked cells are obtained from the cell CELL_q for further search.

[0051] Step S146: When the height decrease ratio parameter is greater than or equal to the preset change threshold, stop expanding outward from the cell and mark the planar position identifier of the cell as the boundary cell of the plant space occupied area.

[0052] In this embodiment, when the height drop ratio parameter R_drop_q calculated in step S144 is greater than or equal to the preset change threshold TH_height, it is determined that the cell CELL_q is located at the boundary position of the plant-occupied area of ​​the crop plant corresponding to the current growth origin O_p. Further expansion search from cell CELL_q outwards is stopped. The planar position identifier of cell CELL_q is added to a set of boundary cells S_boundary_p for the plant space occupied area, serving as a boundary cell of the plant space occupied area of ​​the crop plant corresponding to the current growth origin O_p.

[0053] Step S147: Repeat the operations of expanding outward by one layer, extracting the vertical height value, calculating the height decrease ratio parameter, comparing the preset change threshold, and marking the boundary cells until the expansion process in all search directions stops due to encountering the boundary cells, and obtain the complete plant space occupied area boundary of the crop plant corresponding to the current growth origin. The complete plant space occupied area boundary includes the set of all planar position identifiers marked as boundary cells.

[0054] In this embodiment, steps S143 to S146 are then repeated. Starting from the first search cell, the process expands outwards layer by layer, with each layer's expansion based on cells marked as part of the plant's occupied area in the previous layer. The expansion process terminates when all cells in all expansion directions of a certain layer stop expanding due to a height decrease ratio parameter being greater than or equal to a preset change threshold TH_height. After the expansion process in all search directions has stopped, the planar position identifiers of all boundary cells in the plant's occupied area boundary cell set S_boundary_p are collected. These planar position identifiers are then connected according to their spatial adjacency to form a closed boundary profile, yielding the complete plant's occupied area boundary corresponding to the current growth origin O_p. This complete plant's occupied area boundary includes the set of all planar position identifiers marked as boundary cells.

[0055] Step S148: Obtain the planar position identifiers of all cells enclosed by the boundary of the complete plant space occupied area, extract the corresponding cell's height-spectral joint observation vector from the spatial distribution matrix according to each planar position identifier, collect the vertical height values ​​in the height-spectral joint observation vectors of all enclosed cells, and form the internal vertical height value sequence of the crop plant.

[0056] In this embodiment, all cells within the internal region enclosed by the boundary of the complete plant space are then acquired. The boundary of the complete plant space forms a closed contour, and the cells contained within this closed contour include the boundary cells themselves and all cells inside the boundary. All cells within the closed contour are traversed, and the planar position identifier of each cell is obtained. Based on each planar position identifier, the corresponding cell is located in the spatial distribution matrix, and the vertical height value is extracted from the first dimension of the height-spectral joint observation vector stored in that cell. All extracted vertical height values ​​are organized into a sequence according to the spatial position order of their corresponding cells, resulting in the internal vertical height value sequence Z_seq_p of the crop plant.

[0057] Step S149: Select the vertical height value with the largest value from the internal vertical height value sequence as the top height value of the canopy of the crop plant, and select the vertical height value with the smallest value from the internal vertical height value sequence as the bottom height value of the canopy of the crop plant.

[0058] In this embodiment, the maximum and minimum value filtering operations are then performed on the aforementioned internal vertical height value sequence Z_seq_p. Each vertical height value in the internal vertical height value sequence Z_seq_p is traversed, and the vertical height value Z_top_p with the largest value in the sequence is found by comparison one by one. This Z_top_p is taken as the top height value of the crop canopy. Similarly, the minimum vertical height value Z_bottom_p in the internal vertical height value sequence Z_seq_p is found by comparison one by one, and this Z_bottom_p is taken as the bottom height value of the crop canopy.

[0059] Step S1410: Based on the canopy top height value and the canopy bottom height value, shift the canopy top height value downward by a preset height step to obtain the canopy upper surface boundary height value, and shift the canopy bottom height value upward by a preset height step to obtain the canopy lower surface boundary height value.

[0060] In this embodiment, a preset height step parameter delta_H is then obtained, which is a pre-set height offset. The canopy upper surface boundary height value H_upper_p is calculated by subtracting the preset height step parameter delta_H from the canopy top height value Z_top_p. The canopy lower surface boundary height value H_lower_p is then calculated by adding the preset height step parameter delta_H to the canopy bottom height value Z_bottom_p.

[0061] Step S1411: Extract all LiDAR measurement points corresponding to vertical height values ​​that are greater than or equal to the boundary height value of the canopy upper surface from the internal vertical height value sequence, and form the above LiDAR measurement points into a subset of the point cloud of the canopy upper surface.

[0062] In this embodiment, the process then iterates through each vertical height value in the internal vertical height value sequence Z_seq_p and its corresponding LiDAR measurement point. For each vertical height value Z_q in the internal vertical height value sequence Z_seq_p, the vertical height value Z_q is compared with the canopy upper surface boundary height value H_upper_p. When the vertical height value Z_q is greater than or equal to the canopy upper surface boundary height value H_upper_p, all information of the LiDAR measurement point corresponding to the vertical height value Z_q is obtained, including the horizontal axis coordinate value, the horizontal axis coordinate value, the vertical height value, and the red light band reflection intensity value, green light band reflection intensity value, and near-infrared band reflection intensity value of the optical image pixels located at the same spatial position as the LiDAR measurement point. The LiDAR measurement point is added to a point cloud data set to obtain the canopy upper surface point cloud subset PointCloud_upper_p.

[0063] Step S1412: Extract all LiDAR measurement points corresponding to vertical height values ​​that are less than or equal to the boundary height value of the lower canopy surface from the internal vertical height value sequence, and form the above LiDAR measurement points into a subset of the lower canopy surface point cloud.

[0064] In this embodiment, the process then iterates through each vertical height value in the internal vertical height value sequence Z_seq_p and its corresponding LiDAR measurement point. For each vertical height value Z_q in the internal vertical height value sequence Z_seq_p, the vertical height value Z_q is compared with the lower canopy surface boundary height value H_lower_p. When the vertical height value Z_q is less than or equal to the lower canopy surface boundary height value H_lower_p, all information of the LiDAR measurement point corresponding to the vertical height value Z_q is obtained, including the horizontal axis coordinate value, the horizontal axis coordinate value, the vertical height value, and the red light band reflection intensity value, green light band reflection intensity value, and near-infrared band reflection intensity value of the optical image pixels located at the same spatial position as the LiDAR measurement point. The LiDAR measurement point is added to a point cloud data set to obtain the lower canopy surface point cloud subset PointCloud_lower_p.

[0065] Step S1413: Extract the red band reflection intensity value, green band reflection intensity value and near-infrared band reflection intensity value corresponding to each LiDAR measurement point in the canopy upper surface point cloud subset, calculate the first vegetation index value for each LiDAR measurement point, and remove optical remote sensing anomalies whose first vegetation index value is lower than a preset vegetation threshold from the canopy upper surface point cloud subset according to the first vegetation index value, so as to obtain the optimized canopy upper surface point cloud subset.

[0066] In this embodiment, the above-mentioned canopy upper surface point cloud subset PointCloud_upper_p is then optimized. Each LiDAR measurement point in the canopy upper surface point cloud subset PointCloud_upper_p is traversed, and the corresponding red band reflectance intensity value R_u, green band reflectance intensity value G_u, and near-infrared band reflectance intensity value NIR_u are obtained. The first vegetation index value VI_u for the above-mentioned LiDAR measurement point is calculated. The calculation process for the first vegetation index value VI_u is as follows: calculate the difference between the near-infrared band reflectance intensity value NIR_u and the red band reflectance intensity value R_u; calculate the sum of the near-infrared band reflectance intensity value NIR_u and the red band reflectance intensity value R_u; divide the difference by the sum to obtain the first vegetation index value VI_u. A preset vegetation threshold TH_veg is obtained, which is a value between 0 and 1. The first vegetation index value VI_u is compared with the preset vegetation threshold TH_veg. When the first vegetation index value VI_u is lower than the preset vegetation threshold TH_veg, the LiDAR measurement point is determined to be an optical remote sensing anomaly. The optical remote sensing anomaly may be a non-vegetated point caused by sensor noise, atmospheric interference, or shadow occlusion. The LiDAR measurement points determined to be optical remote sensing anomalies are removed from the canopy upper surface point cloud subset PointCloud_upper_p. After the above traversal and removal operations, the optimized canopy upper surface point cloud subset PointCloud_upper_opt_p is obtained.

[0067] Step S1414: Extract the red band reflectance intensity value, green band reflectance intensity value, and near-infrared band reflectance intensity value corresponding to each LiDAR measurement point in the canopy lower surface point cloud subset. Calculate the second vegetation index value for each LiDAR measurement point. Based on the second vegetation index value, remove optical remote sensing soil points whose second vegetation index value is higher than a preset soil threshold from the canopy lower surface point cloud subset to obtain the optimized canopy lower surface point cloud subset.

[0068] In this embodiment, the aforementioned sub-set of the canopy lower surface point cloud, PointCloud_lower_p, is then optimized. Each LiDAR measurement point in PointCloud_lower_p is traversed, and the corresponding red band reflectance intensity value R_l, green band reflectance intensity value G_l, and near-infrared band reflectance intensity value NIR_l are obtained. The second vegetation index value VI_l for each LiDAR measurement point is calculated. The calculation process for the second vegetation index value VI_l is as follows: the difference between the near-infrared band reflectance intensity value NIR_l and the red band reflectance intensity value R_l is calculated; the sum of the near-infrared band reflectance intensity value NIR_l and the red band reflectance intensity value R_l is calculated; and the difference is divided by the sum to obtain the second vegetation index value VI_l. A preset soil threshold TH_soil is obtained, which is a value between 0 and 1. The second vegetation index value VI_l is compared with the preset soil threshold TH_soil. When the second vegetation index value VI_l is higher than the preset soil threshold TH_soil, the LiDAR measurement point is determined to be an optical remote sensing soil point, which may be a ground soil exposure point caused by canopy porosity. The LiDAR measurement points determined to be optical remote sensing soil points are removed from the canopy lower surface point cloud subset PointCloud_lower_p. After the above traversal and removal operations, the optimized canopy lower surface point cloud subset PointCloud_lower_opt_p is obtained.

[0069] Step S150: Extract the upper surface point cloud subset and the lower surface point cloud subset of each crop plant from the three-dimensional space enclosed by the boundary of the space occupied by each plant. Calculate the plant height parameter of each crop plant based on the difference between the average spatial height value of the upper surface point cloud subset and the average spatial height value of the lower surface point cloud subset. Calculate the canopy diameter parameter of each crop plant based on the maximum horizontal span of the boundary of the space occupied by the plant.

[0070] In this embodiment, steps S1411 to S1414 have completed the extraction of the point cloud subsets of the upper and lower canopy surfaces from the three-dimensional space enclosed by the boundary of the area occupied by each plant. Then, the plant height and canopy diameter parameters of each crop plant are calculated based on the optimized point cloud subsets.

[0071] Step S151: Obtain the vertical height values ​​of all LiDAR measurement points in the optimized canopy upper surface point cloud subset, sum the vertical height values ​​of all LiDAR measurement points in the optimized canopy upper surface point cloud subset, and divide the sum by the total number of LiDAR measurement points in the optimized canopy upper surface point cloud subset to obtain the average height value of the canopy upper surface; and obtain the vertical height values ​​of all LiDAR measurement points in the optimized canopy lower surface point cloud subset, sum the vertical height values ​​of all LiDAR measurement points in the optimized canopy lower surface point cloud subset, and divide the sum by the total number of LiDAR measurement points in the optimized canopy lower surface point cloud subset to obtain the average height value of the canopy lower surface.

[0072] In this embodiment, firstly, the cumulative sum of the vertical height values ​​of all LiDAR measurement points in the optimized canopy upper surface point cloud subset UpSet_p is calculated. Each LiDAR measurement point in UpSet_p is traversed, the vertical height value is read, and each value is accumulated to obtain Z_upper_sum_p. The total number of LiDAR measurement points in UpSet_p, N_upper_p, is counted. H_canopy_top_p = Z_upper_sum_p / N_upper_p is calculated. Next, the cumulative sum of the vertical height values ​​of all LiDAR measurement points in the optimized canopy lower surface point cloud subset LowSet_p is calculated. Each LiDAR measurement point in LowSet_p is traversed, the vertical height value is read, and each value is accumulated to obtain Z_lower_sum_p. The total number of LiDAR measurement points in LowSet_p, N_lower_p, is counted. H_canopy_bottom_p = Z_lower_sum_p / N_lower_p is calculated.

[0073] Step S152: Calculate the difference between the average height of the upper surface of the canopy and the average height of the lower surface of the canopy, and use the difference as the plant height parameter of the crop.

[0074] In this embodiment, the plant height parameter H_plant_p = H_canopy_top_p - H_canopy_bottom_p of the above-mentioned crop plant is then calculated.

[0075] Step S153: Obtain the set of planar position identifiers of all cells enclosed by the boundary of the complete plant space occupied area, and extract the horizontal axis coordinate value and the horizontal axis coordinate value of each planar position identifier in the set of planar position identifiers.

[0076] In this embodiment, the set of planar location identifiers S_region_p for all cells within the internal region enclosed by the boundary of the complete plant space is then obtained. For each planar location identifier ID_rc in the set S_region_p, the horizontal coordinate value X_rc and the horizontal coordinate value Y_rc are parsed from the planar location identifier ID_rc. The parsing process is as follows: the planar location identifier ID_rc is split according to the delimiter string to obtain two substrings; the first substring is converted into a numerical value to obtain the horizontal coordinate value X_rc, and the second substring is converted into a numerical value to obtain the horizontal coordinate value Y_rc.

[0077] Step S154: Select the maximum and minimum horizontal axis coordinate values ​​from all horizontal axis coordinate values, and calculate the difference between the maximum and minimum horizontal axis coordinate values ​​as the first horizontal span value.

[0078] In this embodiment, all horizontal axis coordinate values ​​X_rc obtained in step S153 are then traversed to find the maximum value X_max_p and the minimum value X_min_p. Span_X_p = X_max_p - X_min_p is then calculated.

[0079] Step S155: Select the maximum and minimum horizontal vertical axis coordinate values ​​from all horizontal vertical axis coordinate values, and calculate the difference between the maximum and minimum horizontal vertical axis coordinate values ​​as the second horizontal span value.

[0080] In this embodiment, all horizontal vertical axis coordinate values ​​Y_rc obtained in step S153 are then traversed to find the maximum value Y_max_p and the minimum value Y_min_p. Span_Y_p = Y_max_p - Y_min_p is then calculated.

[0081] Step S156: Compare the first horizontal span value and the second horizontal span value, and select the largest span value as the canopy diameter parameter of the crop plant.

[0082] In this embodiment, Span_X_p is then compared with Span_Y_p. If Span_X_p ≥ Span_Y_p, then D_canopy_p = Span_X_p; otherwise, D_canopy_p = Span_Y_p.

[0083] Step S157: Extract the red band reflection intensity value and near-infrared band reflection intensity value corresponding to each LiDAR measurement point in the optimized canopy upper surface point cloud subset, calculate the ratio vegetation index value of each LiDAR measurement point, calculate the average value of all ratio vegetation index values, and obtain the average vegetation index value of the canopy upper surface.

[0084] In this embodiment, the ratio vegetation index value is then calculated for each LiDAR measurement point in UpSet_p. Each LiDAR measurement point in UpSet_p is traversed to obtain the red light band reflectance value R_u and the near-infrared band reflectance value NIR_u. RVI_u = NIR_u / R_u is calculated. All calculated RVI_u values ​​are summed to obtain RVI_upper_sum_p, and the total number of LiDAR measurement points in UpSet_p, N_upper_opt_p, is counted. RVI_upper_mean_p = RVI_upper_sum_p / N_upper_opt_p is then calculated.

[0085] Step S158: Extract the corresponding values ​​of the subset of point clouds on the lower surface of the canopy, calculate the average vegetation index value of the lower surface of the canopy, input the average vegetation index value of the upper surface of the canopy and the average vegetation index value of the lower surface of the canopy into a preset leaf area index inversion model, and the leaf area index inversion model outputs the leaf area index parameters of the crop plant.

[0086] In this embodiment, the average vegetation index value RVI_lower_mean_p of the lower canopy surface is then calculated for LowSet_p using the same method as in step S157 above. RVI_upper_mean_p and RVI_lower_mean_p are input into a preset leaf area index inversion model, which outputs LAI_p = exp(ln(RVI_upper_mean_p / RVI_lower_mean_p) × coeff).

[0087] Step S210: Obtain the plant height parameters of all crop plants in the target farmland area, arrange the plant height parameters of all crop plants in ascending order of value to obtain a plant height parameter sorting sequence, divide the plant height parameter sorting sequence into multiple plant height parameter intervals at equal intervals, count the number of crop plants contained in each plant height parameter interval, and obtain the plant frequency distribution value of each plant height parameter interval.

[0088] In this embodiment, all H_plant_p are collected, where p=1, ..., P_total. All H_plant_p are sorted in ascending order to obtain H_sorted, with the minimum value being H_min and the maximum value being H_max. The number of intervals K_height is obtained, and the interval length delta_H_interval=(H_max-H_min) / K_height is calculated. For the k-th interval, its range is [H_min+(k-1)×delta_H_interval, H_min+k×delta_H_interval]. The number of plants in each interval is counted to obtain Freq_H_k.

[0089] Step S220: Using the plant height parameter interval as the horizontal axis coordinate value and the plant frequency distribution value of each plant height parameter interval as the vertical axis coordinate value, generate a histogram of plant frequency distribution of the target farmland area. Extract the plant height parameter interval corresponding to the peak value from the plant frequency distribution histogram as the dominant plant height parameter interval. Calculate the median value of the dominant plant height parameter interval as the representative plant height parameter of the target farmland area.

[0090] In this embodiment, a histogram of the distribution of plant height frequency is generated with the midpoint value of each plant height parameter interval as the horizontal axis and Freq_H_k as the vertical axis. The interval corresponding to the peak value is extracted as the dominant plant height parameter interval, and the median value of the interval is calculated as H_rep = (lower limit value + upper limit value) / 2.

[0091] Step S230: Obtain the canopy diameter parameters of all crop plants in the target farmland area, arrange the canopy diameter parameters of all crop plants in ascending order of value to obtain a canopy diameter parameter sorting sequence, divide the canopy diameter parameter sorting sequence into multiple canopy diameter parameter intervals at equal intervals, count the number of crop plants contained in each canopy diameter parameter interval, and obtain the plant frequency distribution value of each canopy diameter parameter interval.

[0092] In this embodiment, the canopy diameter parameter D_canopy_p of all crop plants in the target farmland area is collected, where the value of p ranges from 1 to P_total. All canopy diameter parameters are sorted in ascending order of value, resulting in a canopy diameter parameter sorting sequence D_sorted. The minimum value in the sorted sequence is denoted as D_min, and the maximum value is denoted as D_max. A preset canopy diameter interval quantity parameter K_diameter is obtained, and the sorted sequence is divided into K_diameter canopy diameter parameter intervals at equal intervals. The interval length delta_D_interval of each canopy diameter parameter interval is calculated as follows: the difference between the maximum value D_max and the minimum value D_min is divided by the canopy diameter interval quantity parameter K_diameter to obtain the interval length delta_D_interval. For the k-th canopy diameter parameter interval, its lower limit is D_min plus k-1 times delta_D_interval, and its upper limit is D_min plus k times delta_D_interval. Iterate through each canopy diameter parameter in the above sorted sequence, determine which canopy diameter parameter interval it falls into, and add the corresponding crop plant count to the corresponding canopy diameter parameter interval. Count the number of crop plants contained in each canopy diameter parameter interval to obtain the plant frequency distribution value Freq_D_k for each canopy diameter parameter interval.

[0093] Step S240: Using the canopy diameter parameter interval as the horizontal axis coordinate value and the plant frequency distribution value of each canopy diameter parameter interval as the vertical axis coordinate value, generate a canopy diameter frequency distribution histogram for the target farmland area. Extract the canopy diameter parameter interval corresponding to the peak value from the canopy diameter frequency distribution histogram as the dominant canopy diameter parameter interval. Calculate the median value of the dominant canopy diameter parameter interval as the representative canopy diameter parameter of the target farmland area.

[0094] In this embodiment, the lower limit or midpoint value of the canopy diameter parameter interval obtained in step S230 is used as the horizontal axis coordinate value, and the plant frequency distribution value Freq_D_k of each canopy diameter parameter interval is used as the vertical axis coordinate value to generate a canopy diameter frequency distribution histogram for the target farmland area. The bar with the largest vertical axis coordinate value is identified from the canopy diameter frequency distribution histogram, and the canopy diameter parameter interval corresponding to this bar is taken as the dominant canopy diameter parameter interval. The median value of the dominant canopy diameter parameter interval is calculated by adding the lower limit and upper limit values ​​of the dominant canopy diameter parameter interval and dividing by 2. The median value is then used as the representative canopy diameter parameter D_rep for the target farmland area.

[0095] Step S250: Obtain the leaf area index parameters of all crop plants in the target farmland area, arrange the leaf area index parameters of all crop plants in ascending order of value to obtain a leaf area index parameter sorting sequence, divide the leaf area index parameter sorting sequence into multiple leaf area index parameter intervals at equal intervals, count the number of crop plants contained in each leaf area index parameter interval, and obtain the plant frequency distribution value of each leaf area index parameter interval.

[0096] In this embodiment, the leaf area index (LAI) parameter LAI_p of all crop plants in the target farmland area is collected, where the value of p ranges from 1 to P_total. All LAI parameters are sorted in ascending order of value to obtain a sorted sequence LAI_sorted. The minimum value in the sorted sequence is denoted as LAI_min, and the maximum value is denoted as LAI_max. A preset LAI interval quantity parameter K_lai is obtained, and the sorted sequence is divided into K_lai intervals at equal intervals. The interval length delta_LAI_interval of each interval is calculated as follows: the difference between the maximum value LAI_max and the minimum value LAI_min is divided by the interval quantity parameter K_lai to obtain the interval length delta_LAI_interval. For the k-th leaf area index (LAI) parameter interval, its lower limit is LAI_min plus k-1 times delta_LAI_interval, and its upper limit is LAI_min plus k times delta_LAI_interval. Iterate through each LAI parameter in the sorted sequence, determine which LAI parameter interval it falls into, and add the corresponding crop plant count to the corresponding LAI parameter interval. Count the number of crop plants contained in each LAI parameter interval to obtain the plant frequency distribution value Freq_LAI_k for each LAI parameter interval.

[0097] Step S260: Using the leaf area index parameter interval as the horizontal axis coordinate value and the plant frequency distribution value of each leaf area index parameter interval as the vertical axis coordinate value, generate a leaf area index frequency distribution histogram for the target farmland area. Extract the leaf area index parameter interval corresponding to the peak value from the leaf area index frequency distribution histogram as the dominant leaf area index parameter interval. Calculate the median value of the dominant leaf area index parameter interval as the representative leaf area index parameter of the target farmland area.

[0098] In this embodiment, the lower limit or midpoint value of the leaf area index (LAI) parameter interval obtained in step S250 is used as the horizontal axis coordinate, and the plant frequency distribution value Freq_LAI_k for each LAI parameter interval is used as the vertical axis coordinate to generate a leaf area index frequency distribution histogram for the target farmland area. The bar with the largest vertical axis coordinate value is identified from the LAI frequency distribution histogram, and the LAI parameter interval corresponding to that bar is taken as the dominant LAI parameter interval. The median value of the dominant LAI parameter interval is calculated by adding the lower limit and upper limit values ​​of the dominant LAI parameter interval and dividing by 2. This median value is then used as the representative LAI_rep parameter for the target farmland area.

[0099] Step S270: Based on the representative plant height parameter, the representative canopy diameter parameter, and the representative leaf area index parameter, generate a comprehensive description file of the three-dimensional structure parameters of crops in the target farmland area. The comprehensive description file of the three-dimensional structure parameters of crops includes plant height distribution feature descriptors, plant canopy width distribution feature descriptors, and plant leaf area distribution feature descriptors of the target farmland area.

[0100] In this embodiment, the representative plant height parameter H_rep obtained in step S220, the representative canopy diameter parameter D_rep obtained in step S240, and the representative leaf area index parameter LAI_rep obtained in step S260 are then integrated. A structured file is created as a comprehensive description file of the crop's three-dimensional structural parameters. This comprehensive description file contains three data fields: the first data field is a plant height distribution feature descriptor, storing the representative plant height parameter H_rep and the plant frequency distribution value Freq_H_k for each plant height parameter interval calculated in step S210; the second data field is a plant canopy width distribution feature descriptor, storing the representative canopy diameter parameter D_rep and the plant frequency distribution value Freq_D_k for each canopy diameter parameter interval calculated in step S230; the third data field is a plant leaf area distribution feature descriptor, storing the representative leaf area index parameter LAI_rep and the plant frequency distribution value Freq_LAI_k for each leaf area index parameter interval calculated in step S250.

[0101] Step S280: Output the comprehensive description file of the three-dimensional structure parameters of the crop to the crop growth monitoring terminal device. The crop growth monitoring terminal device adjusts the irrigation amount, fertilizer amount and planting density parameters in subsequent farmland management operations according to the representative plant height parameter, the representative canopy diameter parameter and the representative leaf area index parameter in the comprehensive description file of the three-dimensional structure parameters of the crop.

[0102] In this embodiment, the comprehensive description file of the crop's three-dimensional structural parameters generated in step S270 is then output to a crop growth monitoring terminal device via a wired or wireless communication network. Upon receiving the comprehensive description file, the crop growth monitoring terminal device parses out the representative plant height parameter H_rep, the representative canopy diameter parameter D_rep, and the representative leaf area index parameter LAI_rep. Based on the representative plant height parameter H_rep, the terminal device determines the average height of the crops in the target farmland area. When H_rep is lower than a preset normal growth height range, it generates adjustment commands to increase irrigation and fertilizer application. Based on the representative canopy diameter parameter D_rep, the terminal device determines the canopy spread of the crops in the target farmland area. When D_rep is lower than a preset normal canopy width range, it generates an adjustment command to increase planting density. The aforementioned crop growth monitoring terminal equipment determines the leaf area development status of crops in the target farmland area based on the aforementioned representative leaf area index parameter LAI_rep. When the aforementioned representative leaf area index parameter LAI_rep is lower than the preset normal leaf area index range, an adjustment command to increase the amount of nitrogen fertilizer applied is generated.

[0103] Step S310: Extract the calculated ratio vegetation index values ​​stored in all cells from the spatial distribution matrix of the height-spectral joint observation vector to generate a spatial distribution map of the ratio vegetation index of the target farmland area.

[0104] In this embodiment, based on the spatial distribution matrix of the joint height-spectral observation vector generated in step S125, further identification and segmentation of plant cluster regions are performed. All cells in the spatial distribution matrix are traversed. For each cell storing a valid joint height-spectral observation vector, the ratio vegetation index (NDVI) calculation value is extracted from the fifth dimension position of that joint height-spectral observation vector. The extracted ratio vegetation index calculation values ​​of each cell are arranged according to the row and column indices of that cell in the spatial distribution matrix, generating a two-dimensional matrix of the same size as the spatial distribution matrix, serving as the ratio vegetation index spatial distribution map NDVI_map. The element in the m-th row and n-th column of the ratio vegetation index spatial distribution map NDVI_map stores the ratio vegetation index calculation value NDVI_mn of the cell in the m-th row and n-th column.

[0105] Step S320: Extract the vertical height values ​​stored in all cells from the spatial distribution matrix of the joint height and spectral observation vector to generate a spatial distribution map of the canopy height of the target farmland area.

[0106] In this embodiment, all cells in the aforementioned spatial distribution matrix are then traversed. For each cell storing an effective height-spectral joint observation vector, the vertical height value Z_value is extracted from the first dimension position of that height-spectral joint observation vector. The extracted vertical height values ​​of each cell are arranged according to the row and column indices of that cell in the spatial distribution matrix to generate a two-dimensional matrix of the same size as the aforementioned spatial distribution matrix, which serves as the canopy height spatial distribution map Height_map. The element in the m-th row and n-th column of the aforementioned canopy height spatial distribution map Height_map stores the vertical height value Z_mn of the cell in the m-th row and n-th column.

[0107] Step S330: Divide the calculated value of the ratio vegetation index of each cell in the ratio vegetation index spatial distribution map by the vertical height value of the corresponding cell in the canopy height spatial distribution map to obtain the vegetation cover density ratio of each cell, generate a vegetation cover density ratio spatial distribution map, identify continuous cell regions in the vegetation cover density ratio spatial distribution map where the vegetation cover density ratio exceeds a preset vegetation density threshold, and mark each continuous cell region as an initial plant cluster region.

[0108] In this embodiment, for the cell in the m-th row and n-th column, Density_mn = NDVI_mn / Z_mn is calculated. A preset vegetation density threshold TH_density is obtained, and cells with Density_mn > TH_density are marked as high-density cells. A connected component labeling algorithm is used to connect adjacent high-density cells into a continuous cell region, and each continuous region is marked as an initial plant cluster region Cluster_t.

[0109] Step S340: Obtain the coordinates of the geometric center point of each initial plant cluster region. Using the coordinates of the geometric center point of each initial plant cluster region as the cluster center point, calculate the Euclidean distance value from each cell in the initial plant cluster region to the cluster center point. Arrange the Euclidean distance values ​​of all cells in each initial plant cluster region in ascending order to obtain the distance value sorting sequence of the initial plant cluster region.

[0110] In this embodiment, for the t-th initial plant cluster region Cluster_t, the coordinates of the geometric center point are calculated as: X_center_t = (X_1 + X_2 + ... + X_Num) / Num, Y_center_t = (Y_1 + Y_2 + ... + Y_Num) / Num. For each cell within this region, the Euclidean distance value Dist_uc = sqrt((X_uc - X_center_t)^2 + (Y_uc - Y_center_t)^2). All Dist_uc values ​​are sorted in ascending order to obtain Dist_sorted_t.

[0111] Step S350: Extract the Euclidean distance value located at a preset percentile position from the distance value sorting sequence as the cluster radius parameter of the initial plant cluster region. Generate a circular boundary for each initial plant cluster region based on the cluster radius parameter and the coordinates of the cluster center point. Calculate the overlap area between the circular boundary of each initial plant cluster region and the boundary of the plant space occupied area. When the overlap area ratio between the circular boundary and the boundary of the plant space occupied area exceeds a preset overlap threshold, mark the initial plant cluster region as a single plant region.

[0112] In this embodiment, the Euclidean distance value located at the preset percentile_R position is extracted from Dist_sorted_t as the cluster radius parameter R_cluster_t. A circular boundary is generated with (X_center_t, Y_center_t) as the center and R_cluster_t as the radius. The overlap ratio Overlap_ratio_t between this circular boundary and the boundary of the plant space-occupied area is calculated. If Overlap_ratio_t > TH_overlap, the initial plant cluster area is marked as a single plant area.

[0113] Step S360: When the overlap area ratio between the circular boundary and the boundary of the plant space-occupying area does not exceed the preset overlap threshold, the initial plant cluster area is marked as a multi-plant merged area. Peak detection segmentation operation based on the local height protrusion parameter is performed on the multi-plant merged area to divide the multi-plant merged area into multiple single plant sub-regions. The boundary of each single plant area and each single plant sub-region is redefined as the corrected plant space-occupying area boundary of the crop plant. The corrected plant space-occupying area boundary is used to replace the plant space-occupying area boundary.

[0114] In this embodiment, when the overlap area ratio Overlap_ratio_t calculated in step S350 is less than or equal to the preset overlap threshold TH_overlap, the initial plant cluster region Cluster_initial_t is marked as a multi-plant merge region Cluster_merge_t. A peak detection segmentation operation based on the local height protrusion parameter is performed on the multi-plant merge region Cluster_merge_t to divide it into multiple single-plant sub-regions. The specific steps of the peak detection segmentation operation based on the local height protrusion parameter are described in detail in subsequent steps S361 to S368. The boundary of each single-plant region Cluster_single_t and the boundary of each single-plant sub-region are redefined as the corrected plant space occupancy region boundary of the crop plant. The original plant space occupancy region boundary obtained in step S147 is replaced with the corrected plant space occupancy region boundary as the basis for subsequent parameter extraction.

[0115] Step S361: Extract the planar position identifiers of all cells from the multi-plant merging region, and read the local height protrusion parameter of the corresponding cell from the spatial distribution matrix according to each planar position identifier to generate the local height protrusion parameter distribution matrix of the multi-plant merging region.

[0116] In this embodiment, the aforementioned multi-plant merging region Cluster_merge_t is first processed. A set of planar location identifiers for all cells in Cluster_merge_t is obtained. For each planar location identifier in the set, the corresponding cell is located in the spatial distribution matrix based on the planar location identifier, and the pre-stored local height convexity parameter H_rise is read from the data structure corresponding to that cell. The local height convexity parameters of each cell are arranged according to the relative position of that cell within the multi-plant merging region, generating a local height convexity parameter distribution matrix H_rise_map_t.

[0117] Step S362: Select cells whose local height convexity parameter is greater than a preset peak threshold from the local height convexity parameter distribution matrix as peak candidate cells, and obtain a set of peak candidate cells.

[0118] In this embodiment, a preset peak threshold TH_peak is then obtained. All cells in the local height convexity parameter distribution matrix H_rise_map_t are traversed, and the local height convexity parameter H_rise of each cell is compared with the preset peak threshold TH_peak. When the local height convexity parameter H_rise is greater than the preset peak threshold TH_peak, the cell is marked as a peak candidate cell. All cells marked as peak candidates are collected to obtain the peak candidate cell set Peak_candidates_t.

[0119] Step S363: Extract the comprehensive confidence score of each peak candidate cell from the peak candidate cell set, take the peak candidate cell with the highest comprehensive confidence score as the first plant peak cell, and record the planar position identifier of the first plant peak cell as the first segmentation origin identifier.

[0120] In this embodiment, for each peak candidate cell in the peak candidate cell set Peak_candidates_t, the comprehensive confidence score C_score corresponding to that cell is obtained. The comprehensive confidence score C_score has already been calculated and stored in step S1210. The peak candidate cell with the largest comprehensive confidence score C_score is found from the peak candidate cell set Peak_candidates_t and designated as the first plant peak cell Peak_first_t. The planar position identifier of the first plant peak cell Peak_first_t is obtained and recorded as the first segmentation origin identifier Origin_first_t.

[0121] Step S364: Using the spatial location corresponding to the first segmentation origin identifier as the first growth origin, perform a spatial region expansion operation based on the height spectrum joint observation vector to obtain the boundary of the first expanded region, and mark all cells within the boundary of the first expanded region as sub-region cells belonging to the first plant.

[0122] In this embodiment, the spatial region expansion operation is then performed using the spatial location corresponding to the first segmentation origin identifier Origin_first_t as the first growth origin. The specific steps of the spatial region expansion operation are described in detail in subsequent steps S3641 to S3646. After the spatial region expansion operation is completed, a first expansion region boundary is obtained. The planar position identifiers of all cells enclosed by the first expansion region boundary are obtained, and these cells are marked as sub-region cells belonging to the first plant, and these cells are added to the first sub-region cell set Subregion_first_t.

[0123] Step S365: Delete all peak candidate cells that have been marked as belonging to the sub-region cells of the first plant from the peak candidate cell set, and select the peak candidate cell with the highest comprehensive confidence score from the remaining peak candidate cells as the peak cell of the second plant.

[0124] In this embodiment, all peak candidate cells that already exist in the first subregion cell set Subregion_first_t are then deleted from the peak candidate cell set Peak_candidates_t. After the deletion operation, an updated peak candidate cell set Peak_candidates_updated_t is obtained. From the updated peak candidate cell set Peak_candidates_updated_t, the peak candidate cell with the highest comprehensive confidence score C_score is selected as the second plant peak cell Peak_second_t.

[0125] Step S366: When the overall confidence score of the second plant peak cell exceeds the preset peak retention threshold, the planar position identifier of the second plant peak cell is recorded as the second segmentation origin identifier. The spatial region expansion operation is performed with the spatial position corresponding to the second segmentation origin identifier as the second growth origin to obtain the boundary of the second expanded region.

[0126] In this embodiment, a preset peak retention threshold TH_keep is then obtained. The comprehensive confidence score C_score_second of the peak cell Peak_second_t of the second plant is obtained, and compared with the preset peak retention threshold TH_keep. When the comprehensive confidence score C_score_second is greater than the preset peak retention threshold TH_keep, the planar position identifier of the peak cell Peak_second_t of the second plant is recorded as the second segmentation origin identifier Origin_second_t. Using the spatial position corresponding to the second segmentation origin identifier Origin_second_t as the second growth origin, the same spatial region expansion operation as in step S364 is performed to obtain the boundary of the second expanded region. All cells within the boundary of the second expanded region are marked as sub-region cells belonging to the second plant, and these cells are added to the second sub-region cell set Subregion_second_t.

[0127] Step S367: Repeat the loop of selecting the peak candidate cell with the highest comprehensive confidence score from the peak candidate cell set, performing spatial region expansion operation, marking sub-region cells, and deleting marked peak candidate cells until the comprehensive confidence score of the remaining peak candidate cells is lower than the peak retention threshold.

[0128] In this embodiment, steps S365 and S366 are then repeated. In each loop, the peak candidate cell with the highest overall confidence score is selected from the current peak candidate cell set, and it is checked whether the overall confidence score of the peak candidate cell exceeds the preset peak retention threshold TH_keep. If it does, a spatial region expansion operation is performed with the peak candidate cell as the growth origin, and the cells within the expanded area are marked as sub-region cells belonging to the current plant, and the marked peak candidate cells are deleted from the peak candidate cell set. If the overall confidence score of the currently selected peak candidate cell is lower than or equal to the preset peak retention threshold TH_keep, the loop stops. After the above loop operation, multiple sub-region cell sets Subregion_s are obtained, where s represents the sub-region number.

[0129] Step S368: Collect all marked sub-region cells, take the boundary of each sub-region cell set as the boundary of each single plant sub-region, take the cell set enclosed by the boundary of each single plant sub-region as the internal cell set of that single plant sub-region, take the maximum value of the vertical height value in the height spectrum joint observation vector stored in all cells of the internal cell set of each single plant sub-region as the plant top height value of that single plant sub-region, and take the minimum value of the vertical height value as the plant bottom height value of that single plant sub-region.

[0130] In this embodiment, the set of all subregion cells Subregion_s obtained in steps S364 to S367 is then collected. For each set of subregion cells Subregion_s, the planar position identifiers of all cells in the set are connected to form a closed polygon boundary, which serves as the boundary of the single plant subregion. The set of all cells enclosed by the boundary of the single plant subregion is taken as the internal cell set of the single plant subregion. Each cell in the internal cell set is traversed, and the vertical height value is extracted from the first dimension of the height-spectral joint observation vector stored in the cell. The vertical height value with the largest value is selected from all extracted vertical height values ​​as the plant top height value Z_top_sub_s of the single plant subregion, and the vertical height value with the smallest value is selected as the plant bottom height value Z_bottom_sub_s of the single plant subregion.

[0131] Step S369: Calculate the plant height estimation parameters for each individual plant sub-region based on the plant top height value and plant bottom height value. Then, perform a weighted average of the plant height estimation parameters for each individual plant sub-region and the original plant height parameters of the corresponding peak candidate cell to obtain the final plant height parameters for that individual plant sub-region.

[0132] In this embodiment, for each sub-region of a single plant, H_est_sub_s = Z_top_sub_s - Z_bottom_sub_s is calculated. The original plant height parameter H_plant_original of the peak candidate cell corresponding to this sub-region is obtained. The preset weighting coefficient alpha is obtained, and H_final_sub_s = alpha × H_est_sub_s + (1-alpha) × H_plant_original is calculated.

[0133] Step S3641: Using the cell corresponding to the first segmentation origin identifier as the starting cell for region expansion, obtain the vertical height value in the height-spectral joint observation vector stored in the starting cell for region expansion as the starting point height reference value.

[0134] In this embodiment, the spatial region expansion operation performed in step S364 above is performed according to the following steps. First, the cell CELL_origin corresponding to the first segmentation origin identifier Origin_first_t is obtained. The vertical height value Z_origin is extracted from the first dimension position of the height-spectral joint observation vector V_origin stored in the cell CELL_origin. The vertical height value Z_origin is used as the starting point height reference value H_start.

[0135] Step S3642: Create a set of extended region cells and add the planar position identifier of the starting cell of the extended region to the set of extended region cells; and create a set of extended boundary cells and add the planar position identifier of the starting cell of the extended region to the extended pending queue.

[0136] In this embodiment, an empty set of region expansion cells, S_expand, is then created. The planar position identifier of the region expansion starting cell, CELL_origin, is added to the region expansion cell set S_expand. An empty set of expansion boundary cells, S_boundary_expand, is then created. An expansion processing queue, Queue_process, is created, and the planar position identifier of the region expansion starting cell, CELL_origin, is added to the tail of the expansion processing queue, Queue_process.

[0137] Step S3643: Take out the planar position identifier of the currently processed cell from the extended processing queue, and obtain the set of planar position identifiers of all adjacent cells in the horizontal direction of the currently processed cell.

[0138] In this embodiment, it is then determined whether the extended queue to be processed, Queue_process, is empty. If the extended queue to be processed, Queue_process, is not empty, a planar position identifier is taken from the head of the queue as the current processing cell identifier, ID_current. The corresponding cell, CELL_current, is located based on the current processing cell identifier, ID_current. All adjacent cells in the horizontal direction of CELL_current are obtained. These adjacent cells include the cell above, below, to the left, to the right, and the cells in the four diagonal directions, totaling eight adjacent cells. The planar position identifiers of these eight adjacent cells are collected to obtain the set of adjacent cell identifiers, Neighbor_IDs.

[0139] Step S3644: Traverse each adjacent cell identifier in the set of planar position identifiers of all adjacent cells, and check whether the adjacent cell identifier already exists in the set of extended area cells. When the adjacent cell identifier does not exist in the set of extended area cells, read the vertical height value in the height-spectral joint observation vector stored in the adjacent cell from the spatial distribution matrix as the adjacent cell height value.

[0140] In this embodiment, each neighboring cell identifier ID_neighbor in the neighboring cell identifier set Neighbor_IDs is then traversed. It is checked whether the neighboring cell identifier ID_neighbor already exists in the region expansion cell set S_expand. If the neighboring cell identifier ID_neighbor does not exist in the region expansion cell set S_expand, the corresponding cell CELL_neighbor is searched in the spatial distribution matrix based on the neighboring cell identifier ID_neighbor. The vertical height value Z_neighbor is extracted from the first dimension of the height-spectral joint observation vector V_neighbor stored in the cell CELL_neighbor, and this vertical height value Z_neighbor is used as the neighboring cell height value.

[0141] Step S3645: Calculate the absolute value of the height difference between the height value of the adjacent cell and the height reference value of the starting point. Divide the absolute value of the height difference by the height reference value of the starting point to obtain the relative height change ratio. Compare the relative height change ratio with the preset regional expansion stop ratio threshold. When the relative height change ratio is less than the regional expansion stop ratio threshold, add the identifier of this adjacent cell to the set of region-expanded cells, and add the identifier of this adjacent cell to the end position of the expansion pending queue.

[0142] In this embodiment, then calculate delta_Z_neighbor = |Z_neighbor - H_start|, and calculate R_change_neighbor = delta_Z_neighbor / H_start. Compare R_change_neighbor with the preset regional expansion stop ratio threshold TH_stop. When R_change_neighbor < TH_stop, add the identifier ID_neighbor of this adjacent cell to S_expand, and add ID_neighbor to the end of Queue_process.

[0143] Step S3646: When the relative height change ratio is greater than or equal to the regional expansion stop ratio threshold, add the identifier of this adjacent cell to the set of expansion boundary cells, and do not add the identifier of this adjacent cell to the expansion pending queue.

[0144] In this embodiment, when the above relative height change ratio R_change_neighbor is greater than or equal to the above preset regional expansion stop ratio threshold TH_stop, add the identifier ID_neighbor of this adjacent cell to the above expansion boundary cell set S_boundary_expand, and do not add the identifier ID_neighbor of this adjacent cell to the above expansion pending queue Queue_process. Repeat the operations of the above Step S3643 to the above Step S3646 until the above expansion pending queue Queue_process becomes an empty queue. At this time, the complete set of region-expanded cells S_expand and the set of expansion boundary cells S_boundary_expand are obtained. Connect the planar position identifiers of all cells in the above expansion boundary cell set S_boundary_expand according to the spatial adjacency relationship to form a closed polygon boundary, and use the above closed polygon boundary as the first expansion region boundary.

[0145] For example, the method may further include: step S410: extracting the horizontal coordinate value, horizontal coordinate value and vertical height value of each LiDAR measurement point from the point cloud subset of the canopy upper surface of each crop plant, and inputting the horizontal coordinate value and horizontal coordinate value of all LiDAR measurement points into the Delaunay triangulation algorithm to generate a triangulation model of the canopy upper surface of the crop plant.

[0146] In this embodiment, after identifying and correcting the spatial occupancy boundaries of all crop plants in the target farmland area, the canopy volume parameters of each crop plant are further calculated. For the p-th crop plant in the target farmland area, the optimized canopy upper surface point cloud subset PointCloud_upper_opt_p is obtained. The horizontal x-axis coordinate, horizontal y-axis coordinate, and vertical height values ​​of each LiDAR measurement point are extracted from the optimized canopy upper surface point cloud subset PointCloud_upper_opt_p. All extracted horizontal x-axis and horizontal y-axis coordinate values ​​are used as input data and fed into the Delaunay triangulation algorithm. The Delaunay triangulation algorithm performs the following operations: constructing a triangular mesh based on the input point set, such that the circumcircle of each triangle in the triangular mesh does not contain any other input points, thereby maximizing the minimum interior angle of all triangles. The Delaunay triangulation algorithm outputs a canopy upper surface triangular mesh model Mesh_upper_p composed of multiple triangular facets.

[0147] Step S420: Extract the vertical height values ​​of the three vertices of each triangle from the triangular mesh model of the canopy upper surface, calculate the arithmetic mean of the vertical height values ​​of the three vertices of each triangle as the representative height value of the triangle, multiply the representative height value of each triangle by the area of ​​the triangle and sum them up to obtain the volume contribution value of the canopy upper surface of the crop plant.

[0148] In this embodiment, each triangular facet Tri_u in Mesh_upper_p is then traversed. For Tri_u, the vertical height values ​​Z1_u, Z2_u, and Z3_u of the three vertices are obtained, and Z_tri_mean_u = (Z1_u + Z2_u + Z3_u) / 3 is calculated. The projected area Area_tri_u of Tri_u on the horizontal plane is calculated. Vol_tri_u = Z_tri_mean_u × Area_tri_u is calculated. All Vol_tri_u are summed to obtain Vol_upper_p.

[0149] Step S430: Extract the horizontal x-axis coordinates, horizontal y-axis coordinates, and vertical height values ​​of each LiDAR measurement point from the subset of the canopy lower surface points of each crop plant. Input the horizontal x-axis coordinates and horizontal y-axis coordinates of all LiDAR measurement points into the Delaunay triangulation algorithm to generate a triangulation model of the canopy lower surface of the crop plant.

[0150] In this embodiment, the same operation is then performed on the optimized canopy lower surface point cloud subset PointCloud_lower_opt_p of the p-th crop plant. The horizontal x-axis coordinate, horizontal y-axis coordinate, and vertical height value of each LiDAR measurement point are extracted from the optimized canopy lower surface point cloud subset PointCloud_lower_opt_p. All extracted horizontal x-axis and horizontal y-axis coordinate values ​​are input into the Delaunay triangulation algorithm to generate a canopy lower surface triangular mesh model Mesh_lower_p composed of multiple triangular patches.

[0151] Step S440: Extract the vertical height values ​​of the three vertices of each triangle from the triangular mesh model of the lower canopy surface, calculate the arithmetic mean of the vertical height values ​​of the three vertices of each triangle as the representative height value of the triangle, multiply the representative height value of each triangle by the area of ​​the triangle and sum them up to obtain the volume contribution value of the lower canopy surface of the crop plant.

[0152] In this embodiment, each triangular facet Tri_v in Mesh_lower_p is then traversed. For Tri_v, the vertical height values ​​Z1_v, Z2_v, and Z3_v of the three vertices are obtained, and Z_tri_mean_v = (Z1_v + Z2_v + Z3_v) / 3 is calculated. The projected area Area_tri_v of Tri_v on the horizontal plane is calculated. Vol_tri_v = Z_tri_mean_v × Area_tri_v is calculated. All Vol_tri_v values ​​are summed to obtain Vol_lower_p.

[0153] Step S450: Calculate the difference between the volume contribution value of the upper surface of the canopy and the volume contribution value of the lower surface of the canopy, and use the difference as the canopy solid volume parameter of the crop plant.

[0154] In this embodiment, V_canopy_p is then calculated as V_upper_p - Vol_lower_p.

[0155] Step S460: Extract the red light band reflectance intensity value, green light band reflectance intensity value, and near-infrared band reflectance intensity value corresponding to each LiDAR measurement point from the point cloud subset of the canopy upper surface of each crop plant. Calculate the first vegetation index value for each LiDAR measurement point by dividing the difference between the near-infrared band reflectance intensity value and the red light band reflectance intensity value by the sum of the near-infrared band reflectance intensity value and the red light band reflectance intensity value.

[0156] In this embodiment, for each LiDAR measurement point in UpSet_p, NDVI_u = (NIR_u - R_u) / (NIR_u + R_u) is then calculated.

[0157] Step S470: Extract the red light band reflectance intensity value, green light band reflectance intensity value, and near-infrared band reflectance intensity value corresponding to each LiDAR measurement point from the subset of points on the lower surface of the canopy of each crop plant. Calculate the second vegetation index value for each LiDAR measurement point by dividing the difference between the near-infrared band reflectance intensity value and the red light band reflectance intensity value by the sum of the near-infrared band reflectance intensity value and the red light band reflectance intensity value.

[0158] In this embodiment, for each LiDAR measurement point in LowSet_p, NDVI_l = (NIR_l - R_l) / (NIR_l + R_l) is then calculated.

[0159] Step S480: Calculate the arithmetic mean of the first vegetation index values ​​of all LiDAR measurement points in the point cloud subset of the upper canopy surface of the crop plant to obtain the average vegetation index value of the upper canopy surface of the crop plant. Calculate the arithmetic mean of the second vegetation index values ​​of all LiDAR measurement points in the point cloud subset of the lower canopy surface of the crop plant to obtain the average vegetation index value of the lower canopy surface of the crop plant.

[0160] In this embodiment, NDVI_upper_mean_p is then calculated as (NDVI_u1 + NDVI_u2 + ... + NDVI_uN) / N_upper_opt_p. NDVI_lower_mean_p is then calculated as (NDVI_l1 + NDVI_l2 + ... + NDVI_lN) / ​​N_lower_opt_p. Step S490: Calculate the difference between the average vegetation index value of the upper surface of the canopy and the average vegetation index value of the lower surface of the canopy. Use this difference as the canopy vertical vegetation index gradient value of the crop plant. Input the canopy vertical vegetation index gradient value, the canopy volume parameter, and the plant height parameter into a pre-trained empirical model. The empirical model outputs the leaf area index parameter of the crop plant.

[0161] In this embodiment, Gradient_NDVI_p is then calculated as NDVI_upper_mean_p - NDVI_lower_mean_p. Gradient_NDVI_p, V_canopy_p, and H_plant_p are input into a pre-trained empirical model, which outputs LAI_p = coeff1 × Gradient_NDVI_p + coeff2 × V_canopy_p + coeff3 × H_plant_p + intercept.

[0162] Step S510: Randomly select a preset number of sample plants from the target farmland area, obtain the measured leaf area index parameter of each sample plant, and combine the canopy body volume parameter, the canopy vertical vegetation index gradient value, the plant height parameter and the measured leaf area index parameter of each sample plant into a sample training data point to obtain a sample training data point set.

[0163] In this embodiment, after calculating the leaf area index (LAI) parameters, a nonlinear regression mapping network is further constructed to improve the estimation accuracy of the LAI parameters. A predetermined number of sample plants are randomly selected from the target farmland area, denoted as N_sample. For each sample plant, the actual LAI parameters are obtained using ground-based measurement methods, resulting in the measured LAI_ground_s parameter, where s represents the sample plant's index, ranging from 1 to N_sample. For each sample plant, the canopy volume parameter V_canopy_s, the canopy vertical vegetation index gradient value Gradient_NDVI_s, and the plant height parameter H_plant_s are obtained. The canopy volume parameter V_canopy_s, the canopy vertical vegetation index gradient value Gradient_NDVI_s, the plant height parameter H_plant_s, and the measured LAI_ground_s parameter are combined into a quintuple data point, which serves as the sample training data point Data_s. Collect all sample training data points to obtain the sample training data point set Dataset_train.

[0164] Step S520: Construct a nonlinear regression mapping network, which includes an input layer, a first hidden layer, a second hidden layer, a third hidden layer, and an output layer. The input layer includes three input nodes corresponding to the canopy volume parameter, the canopy vertical vegetation index gradient value, and the plant height parameter, respectively. The first hidden layer includes twelve computing nodes. Each computing node calculates the activation value of the first hidden layer by weighted summing of the output values ​​of the three input nodes of the input layer and then applying an activation function. The second hidden layer includes eight computing nodes. Each computing node calculates the activation value of the second hidden layer by weighted summing of the twelve activation values ​​of the first hidden layer and then applying an activation function. The third hidden layer includes six computing nodes. Each computing node calculates the activation value of the third hidden layer by weighted summing of the eight activation values ​​of the second hidden layer and then applying an activation function. The output layer includes one output node, which outputs the predicted leaf area index parameter by weighted summing of the six activation values ​​of the third hidden layer.

[0165] In this embodiment, a nonlinear regression mapping network is then constructed. The structure of the nonlinear regression mapping network is as follows: The input layer contains three input nodes. The first input node receives the canopy volume parameter V_canopy, the second input node receives the canopy vertical vegetation index gradient value Gradient_NDVI, and the third input node receives the plant height parameter H_plant. The first hidden layer contains twelve computation nodes. Each computation node in the first hidden layer is connected to one of the three input nodes in the input layer, and each connection corresponds to a weighted summation coefficient. For the i-th computation node in the first hidden layer, the input calculation process is as follows: the three input values ​​of the input layer are multiplied by their corresponding weighted summation coefficients and then summed, plus a bias term, to obtain the linear weighted sum of the node. This linear weighted sum is then input into an activation function for nonlinear transformation. The activation function is a modified linear unit function, whose calculation logic is: when the input value is greater than 0, the input value itself is output; when the input value is less than or equal to 0, 0 is output. After calculation by the activation function, the activation value H1_i of the first hidden layer is obtained. The second hidden layer contains eight computation nodes. Each node in the second hidden layer is connected to one of the twelve computation nodes in the first hidden layer, with each connection corresponding to a weighted summation coefficient. For the j-th computation node in the second hidden layer, the input calculation process is as follows: multiply the twelve activation values ​​of the first hidden layer by their corresponding weighted summation coefficients, sum them up, add a bias term to obtain the linear weighted sum of the node, and then use a modified linear unit activation function to calculate the activation value H2_j of the second hidden layer. The third hidden layer contains six computation nodes. Each node in the third hidden layer is connected to one of the eight computation nodes in the second hidden layer, with each connection corresponding to a weighted summation coefficient. For the k-th computation node in the third hidden layer, the input calculation process is as follows: multiply the eight activation values ​​of the second hidden layer by their corresponding weighted summation coefficients, sum them up, add a bias term to obtain the linear weighted sum of the node, and then use a modified linear unit activation function to calculate the activation value H3_k of the third hidden layer. The output layer contains one output node, which is connected to each of the six computation nodes in the third hidden layer. Each connection corresponds to a weighted summation coefficient. The calculation process for the output node is as follows: the six activation values ​​of the third hidden layer are multiplied by their corresponding weighted summation coefficients and then summed, plus a bias term, to obtain the output value as the parameter LAI_pred for predicting the leaf area index.

[0166] Step S530: Input the canopy volume parameters, canopy vertical vegetation index gradient values, and plant height parameters of each sample training data point in the sample training data point set into the input layer of the nonlinear regression mapping network. After sequential calculation by the first hidden layer, the second hidden layer, the third hidden layer, and the output layer, the predicted leaf area index parameters of the sample training data point are obtained.

[0167] In this embodiment, forward propagation calculation is then performed on each sample training data point Data_s in the aforementioned sample training data point set Dataset_train. The canopy volume parameter V_canopy_s, the canopy vertical vegetation index gradient value Gradient_NDVI_s, and the plant height parameter H_plant_s are extracted from the sample training data point Data_s, and these three values ​​are input into the three input nodes of the input layer of the aforementioned nonlinear regression mapping network. Following the calculation process described in step S520, the aforementioned nonlinear regression mapping network sequentially performs weighted summation and activation function calculations on the twelve computational nodes of the first hidden layer, the eight computational nodes of the second hidden layer, the six computational nodes of the third hidden layer, and the output layer, finally outputting the predicted leaf area index parameter LAI_pred_s for that sample training data point.

[0168] Step S540: Calculate the squared difference between the predicted leaf area index parameter and the measured leaf area index parameter of the training data point. Sum the squared differences of all training data points and divide by the total number of training data points to obtain the mean squared error loss value.

[0169] In this embodiment, for each sample training data point Data_s, Loss_s = (LAI_pred_s - LAI_ground_s)^2 is calculated. MSE_loss = (Loss_1 + Loss_2 + ... + Loss_N_sample) / N_sample is then calculated.

[0170] Step S550: Based on the mean squared error loss value, use the gradient descent optimization algorithm to update all weighted summation coefficients of the first hidden layer, the second hidden layer, the third hidden layer and the output layer in the nonlinear regression mapping network until the mean squared error loss value is less than the preset training stopping threshold, and obtain the trained leaf area index nonlinear regression mapping network.

[0171] In this embodiment, the weighted summation coefficients in the nonlinear regression mapping network are then updated based on the mean squared error loss value MSE_loss. Using a gradient descent optimization algorithm, the partial derivative of the mean squared error loss value MSE_loss with respect to each weighted summation coefficient is calculated to obtain the gradient value of each weighted summation coefficient. The current value of each weighted summation coefficient is subtracted from the current value by a preset learning rate parameter multiplied by the gradient value of that weighted summation coefficient, resulting in the updated value of the weighted summation coefficient. The learning rate parameter is a preset positive number. The updated weighted summation coefficients replace the original weighted summation coefficients, completing one iteration of training. Steps S530 to S550 are repeated, with the updated weighted summation coefficients used for forward propagation and loss calculation in each iteration, and the weighted summation coefficients are updated again. A preset training stopping threshold TH_train is obtained. After each iteration, the currently calculated mean squared error loss value MSE_loss is compared with the preset training stopping threshold TH_train. When the mean squared error loss value MSE_loss is less than the preset training stopping threshold TH_train, the iterative training is stopped, and the current nonlinear regression mapping network is used as the trained leaf area index nonlinear regression mapping network.

[0172] Step S560: Deploy the trained leaf area index nonlinear regression mapping network to the remote sensing data batch processing system. After receiving the new set of raw LiDAR remote sensing measurement data and raw optical remote sensing image data of the target farmland area, the remote sensing data batch processing system extracts the canopy volume parameters, canopy vertical vegetation index gradient values ​​and plant height parameters of each crop plant.

[0173] In this embodiment, the leaf area index nonlinear regression mapping network trained in step S550 is then deployed to a remote sensing data batch processing system. This remote sensing data batch processing system has a data receiving interface for receiving a new set of raw LiDAR remote sensing measurement data and a set of raw optical remote sensing image data for the target farmland area. After receiving the two data sets, the remote sensing data batch processing system extracts the canopy volume parameter V_canopy_new, the canopy vertical vegetation index gradient value Gradient_NDVI_new, and the plant height parameter H_plant_new for each crop plant in the new target farmland area, according to the methods described in steps S110 to S156 and S410 to S490.

[0174] Step S570: Input the canopy volume parameters, canopy vertical vegetation index gradient values ​​and plant height parameters of each crop plant into the leaf area index nonlinear regression mapping network, and the leaf area index nonlinear regression mapping network outputs the final leaf area index batch prediction results for each crop plant.

[0175] In this embodiment, the canopy volume parameter V_canopy_new, the canopy vertical vegetation index gradient value Gradient_NDVI_new, and the plant height parameter H_plant_new of each crop plant extracted in step S560 are then input into the trained leaf area index nonlinear regression mapping network. The leaf area index nonlinear regression mapping network processes each combination of input parameters according to the forward propagation calculation process described in step S530, and outputs the corresponding final leaf area index batch prediction result LAI_final_new.

[0176] Step S580: Generate a spatial distribution raster map of the leaf area index of the target farmland area based on the spatial location coordinates of each crop plant using the batch prediction results of the final leaf area index of all crop plants, and output the spatial distribution raster map of the leaf area index to the crop growth monitoring terminal device for display and storage.

[0177] In this embodiment, the final leaf area index (LAI) batch prediction results for all crop plants in the new target farmland area, along with the spatial coordinates of each crop plant, are collected. A raster matrix matching the spatial extent of the target farmland area is created, with the number of rows and columns determined based on the size of the target farmland area and a preset spatial resolution. For each crop plant, its corresponding raster cell in the raster matrix is ​​determined based on its spatial coordinates, and the final LAI batch prediction result for that plant, LAI_final_new, is assigned to that raster cell. For raster cells not covered by any plants, a preset background value is assigned. After the above assignment operations, a spatial distribution raster map of the leaf area index (LAI) for the target farmland area, LAI_map_new, is generated. The above-mentioned leaf area index spatial distribution raster map LAI_map_new is output to a crop growth monitoring terminal device. The crop growth monitoring terminal device displays the above-mentioned leaf area index spatial distribution raster map LAI_map_new on the screen for users to view, and stores the above-mentioned leaf area index spatial distribution raster map LAI_map_new in the local storage medium for subsequent analysis.

[0178] In one exemplary embodiment, a system for extracting three-dimensional structural parameters of crops by integrating LiDAR and optical remote sensing is provided. This system can be a terminal, server, etc., and its internal structure diagram can be as follows: Figure 2 As shown, this system for extracting three-dimensional structural parameters of crops by integrating LiDAR and optical remote sensing includes a processor, memory, input / output interface, communication interface, display unit, and input device. The processor, memory, and input / output interface are connected via a system bus, and the communication interface, display unit, and input device are also connected to the system bus via the input / output interface. The processor provides computational and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The input / output interface is used for exchanging information between the processor and external devices. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, mobile cellular networks, near-field communication, or other technologies. When the computer program is executed by the processor, it implements a method for extracting three-dimensional structural parameters of crops by integrating LiDAR and optical remote sensing. The display unit is used to form a visually visible image and can be a display screen, projection device, or virtual reality imaging device. The display screen can be an LCD screen or an e-ink screen. The input device can be a touch layer covering the display screen, or a button, trackball, or touchpad set on the shell of a crop three-dimensional structure parameter extraction system that integrates LiDAR and optical remote sensing, or an external keyboard, touchpad, or mouse, etc.

[0179] It should be noted that, in order to simplify the description of the present invention and thus help to understand one or more embodiments of the invention, multiple features may sometimes be grouped into one embodiment, drawing or description thereof in the foregoing description of the embodiments of the present invention.

Claims

1. A method for extracting three-dimensional structure parameters of crops by fusing LiDAR and optical remote sensing, characterized in that, The method includes: Acquire the raw LiDAR remote sensing measurement data set and the raw optical remote sensing image data set of the target farmland area; Spatial overlay analysis is performed on the spatial height distribution information of each LiDAR measurement point in the raw LiDAR remote sensing data set and the spectral reflectance intensity information of the optical image pixels at the corresponding spatial location in the raw optical remote sensing data set to obtain the joint observation vector of height spectrum at each spatial location. Based on the local variation parameter between the spatial height distribution information and the spectral reflectance intensity information in the joint height-spectral observation vector, the set of crop plant center point locations is identified in the target farmland area. Taking the center point of each crop plant as the growth origin, the height spectrum joint observation vector of adjacent spatial positions is searched layer by layer in the horizontal direction until the spatial height distribution information in the height spectrum joint observation vector at the search boundary decreases by more than a preset change threshold, thus obtaining the boundary of the plant space occupied by each crop plant. Extract the top and bottom canopy point cloud subsets of each crop plant from the three-dimensional space enclosed by the boundary of the area occupied by each plant. Calculate the plant height parameter of each crop plant based on the difference between the average spatial height value of the top and bottom canopy point cloud subsets. Calculate the canopy diameter parameter of each crop plant based on the maximum horizontal span of the boundary of the area occupied by the plant. 2.The method of claim 1, wherein, The step involves spatially overlaying the spatial height distribution information of each LiDAR measurement point in the raw LiDAR remote sensing data set with the spectral reflectance intensity information of the corresponding optical image pixels in the raw optical remote sensing data set to obtain a joint height-spectral observation vector at each spatial location. This includes: The horizontal x-axis coordinate, horizontal y-axis coordinate, and vertical height value of each LiDAR measurement point are extracted from the raw LiDAR remote sensing measurement data set. The horizontal x-axis coordinate and the horizontal y-axis coordinate are combined to generate a planar location identifier for each LiDAR measurement point. The vertical height value is used as the spatial height distribution information corresponding to the planar location identifier. Extract the pixel row index, pixel column index, and red band, green band, and near-infrared band reflectance intensity values ​​of each optical image pixel from the optical remote sensing raw image data set. Combine the pixel row index and pixel column index to generate a pixel planar position identifier for each optical image pixel. The planar position identifier of the LiDAR measurement point is compared with the pixel planar position identifier of the optical image pixel by a string matching comparison. When the planar position identifier is the same as the pixel planar position identifier, it is determined that the LiDAR measurement point and the optical image pixel are in the same spatial position. The vertical height value of a LiDAR measurement point located in the same spatial position is vector-stitched with the red light band reflection intensity value, green light band reflection intensity value, and near-infrared band reflection intensity value of the optical image pixel to generate a five-dimensional height-spectral joint observation vector. The first dimension of the height-spectral joint observation vector stores the vertical height value, the second dimension stores the red light band reflection intensity value, the third dimension stores the green light band reflection intensity value, the fourth dimension stores the near-infrared band reflection intensity value, and the fifth dimension stores the vegetation index calculated by the ratio of the red light band reflection intensity value and the near-infrared band reflection intensity value. Traverse all spatial locations in the target farmland area, repeatedly perform the string complete match comparison operation and the vector concatenation operation to generate a spatial distribution matrix of the joint height-spectral observation vector covering all spatial locations in the target farmland area. The row index of the spatial distribution matrix corresponds to the spatial location number in the horizontal axis coordinate direction, and the column index of the spatial distribution matrix corresponds to the spatial location number in the horizontal axis coordinate direction. Each cell of the spatial distribution matrix stores one joint height-spectral observation vector. Using each cell in the spatial distribution matrix as the current analysis cell, extract the vertical height value on the first dimension of the joint height-spectral observation vector stored in the current analysis cell, and simultaneously extract the eight adjacent vertical height values ​​on the first dimension of the eight joint height-spectral observation vectors stored in the eight adjacent cells in the horizontal direction of the current analysis cell. Calculate the absolute value of the height difference between the vertical height value of the current analysis cell and the arithmetic mean of the eight adjacent vertical height values, and use the absolute value of the height difference as a parameter of the local height protrusion of the current analysis cell; Extract the calculated value of the ratio vegetation index on the fifth dimension of the joint height-spectral observation vector stored in the current analysis cell, and simultaneously extract the calculated values ​​of the eight adjacent ratio vegetation indices on the fifth dimension of the eight joint height-spectral observation vectors stored in the eight adjacent cells in the horizontal direction of the current analysis cell; Calculate the vegetation index difference between the calculated vegetation index of the current analysis cell and the arithmetic mean of the calculated vegetation index of the eight adjacent ratios, and use the vegetation index difference as a parameter of the local vegetation vigor of the current analysis cell; The local height protrusion parameter and the local vegetation vigor parameter are normalized respectively. The normalized local height protrusion parameter and the normalized local vegetation vigor parameter are multiplied to obtain the plant center point confidence score of the current analysis cell. When the plant center point confidence score exceeds the preset confidence threshold, the spatial position corresponding to the planar position identifier of the current analysis cell is marked as a candidate position of crop plant center point. Local non-maximum suppression is performed on all cells in the spatial distribution matrix that are marked as candidate locations of the crop plant center point. The confidence scores of the plant center point of adjacent candidate locations are compared. The candidate location with the highest confidence score of the plant center point in the local area is retained as the final location of the crop plant center point. All other candidate locations in the local area that do not have the highest confidence score of the plant center point are deleted, resulting in the set of locations of the crop plant center point. 3.The method of claim 2, wherein, The process involves using the center point of each crop plant as the growth origin and searching horizontally outwards layer by layer for the joint height-spectral observation vector at adjacent spatial locations until the decrease in spatial height distribution information in the joint height-spectral observation vector at the search boundary exceeds a preset change threshold. This process yields the boundary of the spatial area occupied by each crop plant, including: Using the center point of each crop plant in the set of crop plant center point locations as the current growth origin, the initial height-spectral joint observation vector of the current growth origin stored in the corresponding cell of the spatial distribution matrix is ​​obtained, and the initial vertical height value in the initial height-spectral joint observation vector is extracted as the height reference value of the growth origin. The cell of the current growth origin in the spatial distribution matrix is ​​used as the first-level search starting cell, and the planar position identifier of the first-level search starting cell is recorded as the first-level cell identifier set of the plant-occupied area; Starting from the first search cell, expand outwards by one layer to reach the second search cell set, which contains all cells that are horizontally adjacent to the first search cell and are not marked as plant-occupied areas. Extract the vertical height value of the joint observation vector of height spectrum stored in each cell of the second-layer search cell set, and calculate the height reduction ratio parameter of the vertical height value of each cell relative to the height reference value of the growth origin; When the height decrease ratio parameter is less than the preset change threshold, the planar position identifier of the cell is added to the second-level cell identifier set of the plant-occupied area, and the search continues to expand outward from the cell to the next level. When the height decrease ratio parameter is greater than or equal to the preset change threshold, the expansion from that cell outward stops, and the planar position identifier of that cell is marked as the boundary cell of the plant space occupied area; Repeat the operations of expanding outward by one layer, extracting the vertical height value, calculating the height decrease ratio parameter, comparing the preset change threshold, and marking the boundary cells until the expansion process in all search directions stops due to encountering the boundary cells, and obtain the complete plant space occupied area boundary of the crop plant corresponding to the current growth origin. The complete plant space occupied area boundary includes the set of planar position identifiers of all the cells marked as boundary cells. Obtain the planar position identifiers of all cells enclosed by the boundary of the complete plant space. Extract the corresponding cell's height-spectral joint observation vector from the spatial distribution matrix based on each planar position identifier. Collect the vertical height values ​​in the height-spectral joint observation vectors of all enclosed cells to form the internal vertical height value sequence of the crop plant. The vertical height value with the largest value is selected from the internal vertical height value sequence as the top height value of the canopy of the crop plant, and the vertical height value with the smallest value is selected from the internal vertical height value sequence as the bottom height value of the canopy of the crop plant. Based on the canopy top height value and the canopy bottom height value, the canopy top height value is shifted downward by a preset height step to obtain the canopy upper surface boundary height value, and the canopy bottom height value is shifted upward by a preset height step to obtain the canopy lower surface boundary height value; Extract all LiDAR measurement points corresponding to vertical height values ​​that are greater than or equal to the boundary height value of the canopy upper surface from the internal vertical height value sequence, and form the above LiDAR measurement points into a subset of the point cloud of the canopy upper surface; Extract all LiDAR measurement points corresponding to vertical height values ​​that are less than or equal to the boundary height value of the lower canopy surface from the internal vertical height value sequence, and form the above LiDAR measurement points into a subset of the lower canopy surface point cloud; Extract the red band reflectance intensity value, green band reflectance intensity value, and near-infrared band reflectance intensity value corresponding to each LiDAR measurement point in the canopy upper surface point cloud subset. Calculate the first vegetation index value for each LiDAR measurement point. Based on the first vegetation index value, remove optical remote sensing anomalies from the canopy upper surface point cloud subset whose first vegetation index value is lower than a preset vegetation threshold to obtain an optimized canopy upper surface point cloud subset. Extract the red band reflectance intensity value, green band reflectance intensity value, and near-infrared band reflectance intensity value corresponding to each LiDAR measurement point in the canopy lower surface point cloud subset. Calculate the second vegetation index value for each LiDAR measurement point. Based on the second vegetation index value, remove optical remote sensing soil points from the canopy lower surface point cloud subset whose second vegetation index value is higher than a preset soil threshold to obtain an optimized canopy lower surface point cloud subset. 4.The method of claim 3, wherein, The process involves extracting a subset of the canopy upper surface point cloud and a subset of the canopy lower surface point cloud from the three-dimensional space enclosed by the boundary of the area occupied by each plant. The plant height parameter is calculated based on the difference between the average spatial height of the upper surface point cloud subset and the average spatial height of the lower surface point cloud subset. The canopy diameter parameter is calculated based on the maximum horizontal span of the boundary of the area occupied by each plant. The vertical height values ​​of all LiDAR measurement points in the optimized canopy upper surface point cloud subset are obtained. The sum of these vertical height values ​​is then divided by the total number of LiDAR measurement points in the optimized canopy upper surface point cloud subset to obtain the average height value of the canopy upper surface. Similarly, the vertical height values ​​of all LiDAR measurement points in the optimized canopy lower surface point cloud subset are obtained. The sum of these vertical height values ​​is then divided by the total number of LiDAR measurement points in the optimized canopy lower surface point cloud subset to obtain the average height value of the canopy lower surface. Calculate the difference between the average height of the upper surface of the canopy and the average height of the lower surface of the canopy, and use the difference as the plant height parameter of the crop. Obtain the set of planar position identifiers for all cells enclosed by the boundary of the complete plant space occupied area, and extract the horizontal and vertical coordinate values ​​of each planar position identifier in the set of planar position identifiers. The maximum and minimum horizontal axis coordinate values ​​are selected from all horizontal axis coordinate values, and the difference between the maximum and minimum horizontal axis coordinate values ​​is calculated as the first horizontal span value. The maximum and minimum horizontal vertical axis coordinate values ​​are selected from all horizontal vertical axis coordinate values, and the difference between the maximum and minimum horizontal vertical axis coordinate values ​​is calculated as the second horizontal span value. Compare the first horizontal span value and the second horizontal span value, and select the largest span value as the canopy diameter parameter of the crop plant; Extract the red band reflectance intensity value and near-infrared band reflectance intensity value corresponding to each LiDAR measurement point in the optimized canopy upper surface point cloud subset, calculate the ratio vegetation index value of each LiDAR measurement point, calculate the average value of all ratio vegetation index values, and obtain the average vegetation index value of the canopy upper surface. Extract the corresponding values ​​of the subset of point clouds on the lower surface of the canopy, calculate the average vegetation index value of the lower surface of the canopy, input the average vegetation index value of the upper surface of the canopy and the average vegetation index value of the lower surface of the canopy into a preset leaf area index inversion model, and the leaf area index inversion model outputs the leaf area index parameter of the crop plant. 5.The method of claim 1, wherein, After calculating the canopy diameter parameter of each crop plant based on the maximum horizontal span of the area occupied by the plant space, the method further includes: Obtain the plant height parameters of all crop plants in the target farmland area, arrange the plant height parameters of all crop plants in ascending order of value to obtain the plant height parameter sorting sequence, divide the plant height parameter sorting sequence into multiple plant height parameter intervals at equal intervals, count the number of crop plants contained in each plant height parameter interval, and obtain the plant frequency distribution value of each plant height parameter interval. Using the plant height parameter interval as the horizontal axis coordinate value and the plant frequency distribution value of each plant height parameter interval as the vertical axis coordinate value, a histogram of plant frequency distribution of the target farmland area is generated. The plant height parameter interval corresponding to the peak value is extracted from the plant frequency distribution histogram as the dominant plant height parameter interval. The median value of the dominant plant height parameter interval is calculated as the representative plant height parameter of the target farmland area. Obtain the canopy diameter parameters of all crop plants in the target farmland area, arrange the canopy diameter parameters of all crop plants in ascending order of value to obtain the canopy diameter parameter sorting sequence, divide the canopy diameter parameter sorting sequence into multiple canopy diameter parameter intervals at equal intervals, count the number of crop plants contained in each canopy diameter parameter interval, and obtain the plant frequency distribution value of each canopy diameter parameter interval; Using the canopy diameter parameter interval as the horizontal axis coordinate value and the plant frequency distribution value of each canopy diameter parameter interval as the vertical axis coordinate value, a canopy diameter frequency distribution histogram of the target farmland area is generated. The canopy diameter parameter interval corresponding to the peak value is extracted from the canopy diameter frequency distribution histogram as the dominant canopy diameter parameter interval. The median value of the dominant canopy diameter parameter interval is calculated as the representative canopy diameter parameter of the target farmland area. Obtain the leaf area index parameters of all crop plants in the target farmland area, arrange the leaf area index parameters of all crop plants in ascending order of value to obtain the leaf area index parameter sorting sequence, divide the leaf area index parameter sorting sequence into multiple leaf area index parameter intervals at equal intervals, count the number of crop plants contained in each leaf area index parameter interval, and obtain the plant frequency distribution value of each leaf area index parameter interval. Using the leaf area index parameter interval as the horizontal axis coordinate value and the plant frequency distribution value of each leaf area index parameter interval as the vertical axis coordinate value, a leaf area index frequency distribution histogram of the target farmland area is generated. The leaf area index parameter interval corresponding to the peak value is extracted from the leaf area index frequency distribution histogram as the dominant leaf area index parameter interval. The median value of the dominant leaf area index parameter interval is calculated as the representative leaf area index parameter of the target farmland area. Based on the representative plant height parameter, the representative canopy diameter parameter, and the representative leaf area index parameter, a comprehensive description file of the three-dimensional structure parameters of crops in the target farmland area is generated. The comprehensive description file of the three-dimensional structure parameters of crops includes plant height distribution feature descriptors, plant canopy width distribution feature descriptors, and plant leaf area distribution feature descriptors of the target farmland area. The comprehensive description file of the crop's three-dimensional structural parameters is output to the crop growth monitoring terminal device. The crop growth monitoring terminal device adjusts the irrigation amount, fertilizer amount and planting density parameters in subsequent farmland management operations based on the representative plant height parameter, the representative canopy diameter parameter and the representative leaf area index parameter in the comprehensive description file of the crop's three-dimensional structural parameters. 6.The method of claim 2, wherein, The method further includes: Extract the calculated ratio vegetation index values ​​stored in all cells of the spatial distribution matrix of the joint height and spectral observation vector to generate a spatial distribution map of the ratio vegetation index of the target farmland area; Extract the vertical height values ​​stored in all cells from the spatial distribution matrix of the joint height and spectral observation vectors to generate a spatial distribution map of the canopy height of the target farmland area; Divide the calculated value of the vegetation index of each cell in the ratio vegetation index spatial distribution map by the vertical height value of the corresponding cell in the canopy height spatial distribution map to obtain the vegetation cover density ratio of each cell, generate a vegetation cover density ratio spatial distribution map, identify continuous cell regions in the vegetation cover density ratio spatial distribution map where the vegetation cover density ratio exceeds a preset vegetation density threshold, and mark each continuous cell region as an initial plant cluster region; Obtain the coordinates of the geometric center point of each initial plant cluster region. Using the coordinates of the geometric center point of each initial plant cluster region as the cluster center point, calculate the Euclidean distance value from each cell in the initial plant cluster region to the cluster center point. Arrange the Euclidean distance values ​​of all cells in each initial plant cluster region in ascending order to obtain the distance value sorting sequence of the initial plant cluster region. Extract the Euclidean distance value located at a preset percentile position from the distance value sorting sequence as the cluster radius parameter of the initial plant cluster region. Generate a circular boundary for each initial plant cluster region based on the cluster radius parameter and the coordinates of the cluster center point. Calculate the overlap area between the circular boundary of each initial plant cluster region and the boundary of the plant space occupied area. When the overlap area ratio between the circular boundary and the boundary of the plant space occupied area exceeds a preset overlap threshold, mark the initial plant cluster region as a single plant region. When the overlap ratio between the circular boundary and the boundary of the plant space-occupying area does not exceed the preset overlap threshold, the initial plant cluster area is marked as a multi-plant merged area. A peak detection segmentation operation based on the local height protrusion parameter is performed on the multi-plant merged area to divide the multi-plant merged area into multiple single plant sub-regions. The boundary of each single plant area and each single plant sub-region is redefined as the corrected plant space-occupying area boundary of the crop plant. The corrected plant space-occupying area boundary is used to replace the plant space-occupying area boundary.

7. The method of claim 6, wherein the method further comprises: The step of performing peak detection segmentation based on the local height protrusion parameter on the merged multi-plant region, dividing the merged multi-plant region into multiple single-plant sub-regions, includes: Extract the planar position identifiers of all cells in the multi-plant merging region, and read the local height convexity parameter of the corresponding cell from the spatial distribution matrix based on each planar position identifier to generate the local height convexity parameter distribution matrix of the multi-plant merging region. Cells whose local height protrusion parameter is greater than a preset peak threshold are selected from the local height protrusion parameter distribution matrix as peak candidate cells, thus obtaining a peak candidate cell set; Extract the comprehensive confidence score of each peak candidate cell from the peak candidate cell set, take the peak candidate cell with the highest comprehensive confidence score as the first plant peak cell, and record the planar position identifier of the first plant peak cell as the first segmentation origin identifier; Using the spatial location corresponding to the first segmentation origin identifier as the first growth origin, a spatial region expansion operation based on the height spectrum joint observation vector is performed to obtain the boundary of the first expanded region. All cells within the boundary of the first expanded region are marked as sub-region cells belonging to the first plant. Remove all peak candidate cells that have been marked as belonging to the sub-region cells of the first plant from the peak candidate cell set, and select the peak candidate cell with the highest comprehensive confidence score from the remaining peak candidate cells as the peak cell of the second plant. When the overall confidence score of the second plant peak cell exceeds the preset peak retention threshold, the planar position identifier of the second plant peak cell is recorded as the second segmentation origin identifier. The spatial region expansion operation is performed with the spatial position corresponding to the second segmentation origin identifier as the second growth origin to obtain the boundary of the second expansion region. Repeatedly execute the loop of selecting the peak candidate cell with the highest comprehensive confidence score from the peak candidate cell set, performing spatial region expansion operation, marking sub-region cells, and deleting marked peak candidate cells until the comprehensive confidence scores of all remaining peak candidate cells are lower than the peak retention threshold; Collect all marked sub-region cells, take the boundary of each sub-region cell set as the boundary of each single plant sub-region, take the cell set enclosed by the boundary of each single plant sub-region as the internal cell set of that single plant sub-region, take the maximum value of the vertical height value in the height spectrum joint observation vector stored in all cells of the internal cell set of each single plant sub-region as the plant top height value of that single plant sub-region, and take the minimum value of the vertical height value as the plant bottom height value of that single plant sub-region; The plant height estimation parameters for each individual plant subregion are calculated based on the plant top height and plant bottom height values. The plant height estimation parameters for each individual plant subregion are then weighted and averaged with the original plant height parameters of the corresponding peak candidate cells to obtain the final plant height parameters for that individual plant subregion. 8.The method of claim 7, wherein, The step of performing a spatial region expansion operation based on the joint height-spectral observation vector to obtain the boundary of the first expanded region includes: Using the cell corresponding to the first segmentation origin identifier as the starting cell for region expansion, the vertical height value in the height-spectral joint observation vector stored in the starting cell for region expansion is obtained as the starting point height reference value. Create a set of extended region cells and add the planar position identifier of the starting cell of the extended region to the set of extended region cells; and create a set of extended boundary cells and add the planar position identifier of the starting cell of the extended region to the extended pending queue. Take the planar position identifier of the currently processed cell from the extended queue of pending processing, and obtain the set of planar position identifiers of all adjacent cells in the horizontal direction of the currently processed cell; Iterate through each adjacent cell identifier in the set of planar location identifiers of all adjacent cells, and check whether the adjacent cell identifier already exists in the set of extended area cells. When the adjacent cell identifier does not exist in the set of extended area cells, read the vertical height value in the height-spectral joint observation vector stored in the adjacent cell from the spatial distribution matrix as the adjacent cell height value. Calculate the absolute value of the height difference between the adjacent cell height value and the starting point height reference value. Divide the absolute value of the height difference by the starting point height reference value to obtain the relative height change ratio. Compare the relative height change ratio with a preset area expansion stop ratio threshold. When the relative height change ratio is less than the area expansion stop ratio threshold, add the adjacent cell identifier to the area expansion cell set and add the adjacent cell identifier to the end of the expansion pending queue. When the relative height change ratio is greater than or equal to the area expansion stop ratio threshold, the adjacent cell identifier is added to the expansion boundary cell set, but not added to the expansion pending queue. Repeat the operation of retrieving the currently processed cell from the extended processing queue until the extended processing queue becomes an empty queue, thereby obtaining the complete set of extended region cells and the set of extended boundary cells. Connect the planar position identifiers of all cells in the set of extended boundary cells to form a closed polygon boundary, and use the closed polygon boundary as the boundary of the first extended region. 9.A system for extracting three-dimensional structure parameters of crops by fusing LiDAR and optical remote sensing, characterized in that, include: processor; A machine-readable storage medium for storing machine-executable instructions of the processor; The processor is configured to execute the crop three-dimensional structure parameter extraction method that integrates LiDAR and optical remote sensing as described in any one of claims 1 to 8 by executing the machine-executable instructions.

10. A computer program product, characterised in that, The computer program product includes machine-executable instructions stored in a computer-readable storage medium. The processor of the crop three-dimensional structure parameter extraction system integrating LiDAR and optical remote sensing reads the machine-executable instructions from the computer-readable storage medium and executes the machine-executable instructions, causing the crop three-dimensional structure parameter extraction system integrating LiDAR and optical remote sensing to perform the crop three-dimensional structure parameter extraction method integrating LiDAR and optical remote sensing as described in any one of claims 1 to 8.

Images