A soil property high-precision mapping method and system fusing air-ground imaging hyperspectral technology to optimize sampling number

Abstract

The application discloses a kind of fusion space-ground imaging hyperspectral technology optimization sampling number's soil attribute high-precision mapping method and system, belong to soil monitoring technical field.The method includes: collecting regional unmanned aerial vehicle hyperspectral image and pre-processing;Collecting sample point soil sample and obtaining its ground imaging hyperspectral data;Based on ground imaging hyperspectral data, filter optimization sample point subset, determine its effective phosphorus content and train prediction model;Combining vegetation index and maximum interclass variance method segmentation extracts soil area in unmanned aerial vehicle image;Ground imaging hyperspectral and unmanned aerial vehicle soil spectrum are carried out band alignment and reflectivity correction, after the corrected unmanned aerial vehicle hyperspectral data is supplemented, input soil effective phosphorus prediction model and is inverted mapping.The application solves the problem of insufficient coverage of ground imaging data and unmanned aerial vehicle data susceptible to environment, provides key technical support for intelligent agriculture and soil nutrient precision management.

Description

TECHNICAL FIELD

[0001] The application relates to the technical field of soil monitoring, in particular to a soil property high-precision mapping method and system fusing air-ground imaging hyperspectral technology to optimize sampling number. BACKGROUND

[0002] Soil available phosphorus (AP) is a form of phosphorus that can be directly absorbed by plants and is crucial for crop growth and development. However, traditional monitoring methods rely on laboratory chemical analysis (such as the Olsen method or the Bray method), which is time-consuming and labor-intensive, and may cause environmental pollution due to the use of chemical reagents. Moreover, it is difficult to meet the needs of large-scale and high-density sampling. In addition, due to the solubility and easy solubility of soil available phosphorus, it will migrate out of the soil system in a short time, making it challenging to monitor its dynamic changes and make high-precision field-scale predictions.

[0003] In recent years, the development of near-ground sensing technology, especially hyperspectral technology, has made it possible to quickly and non-destructively monitor soil properties. High-spectral sensors mounted on unmanned aerial vehicles can achieve rapid regional coverage, but their data are easily affected by environmental factors such as light and wind speed. Moreover, when there are background elements such as vegetation in the collection area, it will seriously affect the continuity of the image. Ground imaging hyperspectral technology can provide more accurate and controlled spectral data, but it cannot achieve regional continuous monitoring. In addition, the design of soil sampling scheme directly affects the accuracy and cost of DSM. Traditional grid sampling often requires a large number of samples to ensure spatial representativeness at the field scale, resulting in high economic and time costs.

[0004] In summary, developing a soil available phosphorus rapid monitoring scheme that can optimize sample point layout in the presence of vegetation coverage and effectively combine multi-source (ground imaging and unmanned aerial vehicle) hyperspectral data is the key to solving the current demand for dynamic monitoring of soil nutrients in precision agriculture. To address the problems of low monitoring cost efficiency, insufficient accuracy, and imbalance between sampling cost and benefit in existing technologies. SUMMARY

[0005] The purpose of the present application is to provide a soil property high-precision mapping method and system fusing air-ground imaging hyperspectral technology to optimize sampling number, which realizes high-precision and low-cost rapid monitoring of soil available phosphorus by optimizing sample points and data fusion.

[0006] To achieve the above-mentioned purpose, the technical scheme adopted by the present application is as follows:

[0007] A soil property high-precision mapping method fusing air-ground imaging hyperspectral technology to optimize sampling number, comprising the following steps:

[0008] S1: Collecting the unmanned aerial vehicle hyperspectral image with 5cm spatial resolution in the region and pre-processing, the pre-processing including geometric correction, radiation calibration and reflectivity conversion;

[0009] S2: On the day of unmanned aerial vehicle image collection, laying out the training set samples and validation set samples in the region, collecting the soil samples in the sample 1m 2 and mixing uniformly, and placing the soil samples of the sample in the petri dish after pre-processing, obtaining the ground imaging hyperspectral data thereof;

[0010] S3: Based on the ground imaging hyperspectral data, using the conditional Latin hypercube sampling method to select one fourth of the sample optimization sample subset from the training set sample, determining the soil available phosphorus content of the optimization sample subset, and based on the ground imaging hyperspectral data and the corresponding soil available phosphorus content of the optimization sample subset, training to obtain the soil available phosphorus prediction model;

[0011] S4: Image classification is performed on the unmanned aerial vehicle hyperspectral image, and the soil area unmanned aerial vehicle hyperspectral image is segmented by combining the vegetation index and the maximum inter-class variance method, and the spectral reflectivity of the soil area is extracted as the unmanned aerial vehicle soil spectral reflectivity, which is resampled to 1m spatial resolution;

[0012] S5: The ground imaging hyperspectral data and the unmanned aerial vehicle soil spectral reflectivity are aligned in waveband, and the waveband correction coefficient is calculated based on the ground imaging hyperspectral data, and the unmanned aerial vehicle soil spectral reflectivity is corrected to obtain the corrected unmanned aerial vehicle hyperspectral reflectivity data;

[0013] S6: Based on the soil area unmanned aerial vehicle hyperspectral image and the corrected unmanned aerial vehicle hyperspectral reflectivity data, the image is filled to complete, continuous and corrected unmanned aerial vehicle hyperspectral image, which is input into the soil available phosphorus prediction model for inversion, and the soil available phosphorus inversion distribution map of the region is obtained.

[0014] Further, in S2:

[0015] The soil sample is collected, with the sample as the center 1m 2 range, 5 different positions of surface soil samples are collected by using the waffle sampling method and mixed uniformly;

[0016] The training set sample is laid out by using the 50m x 50m grid sampling method;

[0017] The validation set sample is laid out by using the random sampling method;

[0018] The sample is placed in a culture dish, specifically: the sample is homogenized to break up the cohesive soil blocks and aggregates; the treated sample is loaded into a standard culture dish, and a vertical pressure is applied for moderate compaction; a flat scraper is used to level the upper surface of the culture dish.

[0019] Further, the conditional Latin hypercube sampling method in S3 comprises:

[0020] Constructing a standardized covariance matrix of the hyperspectral data of the training set sample wet soil samples;

[0021] Layered cutting and Latin hypercube initialization are performed for each spectral band.

[0022] An optimized sample subset is selected by iteratively optimizing an objective function, which includes a distribution matching term for minimizing the spectral distribution difference between the subset and the full sample set, and a geographical dispersion penalty term for promoting the uniformity of the spatial distribution of the sample points.

[0023] The spectral distribution similarity between the sample subset and the full sample set is measured based on the JS divergence, and when the divergence value is less than a preset threshold, the subset is determined to replace the full sample set for modeling.

[0024] Further, the distribution matching term uses the Kolmogorov-Smirnov statistic to measure the cumulative distribution difference between the subset and the full sample set at each band.

[0025] The geographical dispersion penalty term is calculated based on the Euclidean distance between the sample points.

[0026] Further, in S4:

[0027] The vegetation index includes a normalized vegetation index and a differential vegetation index.

[0028] The maximum inter-class variance method adaptively determines the segmentation threshold according to the gray level characteristics of the image.

[0029] Further, in S5, the band alignment uses cubic spline interpolation and natural boundary method to make the ground imaging hyperspectral data and the UAV soil spectral reflectance have the same band wavelength range and interval.

[0030] Further, the calculation of the band correction coefficient in S5 comprises:

[0031] For each sample and each band, the ratio of its ground imaging hyperspectral reflectance to UAV hyperspectral reflectance is calculated.

[0032] For each band, the ratios of all samples at this band are averaged, and the obtained average value is the correction coefficient of this band.

[0033] Further, the image supplement in S6 is based on the empty area of the unmanned aerial high-spectral image after correction, and the average value of the unmanned aerial high-spectral reflectance data in the range of 9m 2 is replaced.

[0034] Another object of the present application is to provide a soil property high-precision mapping system that fuses air-ground imaging hyperspectral technology to optimize the number of samples, which realizes the soil property high-precision mapping method that fuses air-ground imaging hyperspectral technology to optimize the number of samples when the system is executed, comprising:

[0035] A data acquisition module is configured to acquire unmanned aerial high-spectral images, ground imaging hyperspectral data, and effective phosphorus content of soil samples in a region.

[0036] A sample optimization and model training module is configured to perform conditional Latin hypercube sampling to screen an optimized sample subset and train a soil effective phosphorus prediction model based on the optimized sample data.

[0037] An image processing and segmentation module is configured to process unmanned aerial high-spectral images, segment soil regions, and extract unmanned aerial soil spectral reflectance.

[0038] A data correction module is configured to perform band alignment and reflectance correction on ground imaging and unmanned aerial high-spectral data.

[0039] An inversion mapping module is configured to generate a soil effective phosphorus inversion distribution map using the corrected data and the trained model.

[0040] Further, the system further comprises a data storage module configured to store the soil effective phosphorus prediction model, correction coefficients, and intermediate processing data.

[0041] The present application provides a soil property high-precision mapping method and system that fuses air-ground imaging hyperspectral technology to optimize the number of samples, establishes a soil effective phosphorus high-precision mapping technical framework based on air-ground imaging hyperspectral data fusion from “image preprocessing-spectral correction-sample optimization-cooperative inversion”, effectively solves the industry pain points of high cost, low precision, and low efficiency, provides reliable technical support for rapid, accurate, and low-cost monitoring of soil nutrients in precision agriculture, and has specific beneficial effects, including:

[0042] The unification of high precision and low cost is realized, the prediction precision is guaranteed and improved while the number of field samples is significantly reduced to one fourth by driving conditional Latin hypercube sampling with ground imaging hyperspectral data, the human and material costs of soil investigation are greatly reduced, and the inherent contradiction between sampling cost and model precision of traditional methods such as geostatistical interpolation is solved.

[0043] The application innovates a rapid segmentation method suitable for unmanned aerial vehicle images, realizes rapid, self-adaptive and high-precision segmentation of soil and vegetation by fusing the calculation efficiency of the vegetation index and the objective accuracy of the threshold value determined by the maximum inter-class variance method, and significantly improves the precision of high-precision mapping of soil available phosphorus by using unmanned aerial vehicle hyperspectral under vegetation coverage.

[0044] A technology framework of complementary advantages of multi-source sensing data of air and ground is constructed, the technical advantages of ground imaging hyperspectral (high precision on "point") and unmanned aerial vehicle hyperspectral (wide coverage on "surface") are creatively fused, through band alignment and sub-band correction algorithm, the shortcomings that the unmanned aerial vehicle data is easy to be disturbed by the environment are overcome, and the deficiency that the near-ground sensing data of point measurement cannot be regionally mapped is made up, a collaborative inversion mechanism of "point-surface combination" is formed, and the effect of 1+1>2 is realized.

[0045] A paradigm for rapid monitoring of soil available phosphorus is provided, the advantages of reducing the number of sample points and the rapid data acquisition of the unmanned aerial vehicle hyperspectral platform are combined, the period of high-precision prediction of field scale soil available phosphorus is greatly shortened, and inversion mapping is completed before the change, which provides ideas and methods for low-cost rapid monitoring of easy-to-move soil properties for meter-level precision agriculture. BRIEF DESCRIPTION OF DRAWINGS

[0046] Figure 1 The method flowchart of the application;

[0047] Figure 2 The segmentation of the vegetation and soil parts of the application is shown in the result graph, wherein (a) is the normalized vegetation index NDVI calculation result graph, (c) is the differential vegetation index DVI calculation result graph, (b) is the binary graph after threshold segmentation of (a), and (d) is the binary graph after threshold segmentation of (c);

[0048] Figure 3 The soil available phosphorus inversion distribution graph of the application. DETAILED DESCRIPTION

[0049] In order to make the purpose, technical scheme and advantages of the application more clear and understandable, the application will be further described in detail below in combination with the drawings and examples. It should be understood that the specific examples described herein are only used to explain the application, and are not used to limit the application. In addition, the technical features involved in each embodiment of the application described below can be combined with each other as long as they do not conflict with each other.

[0050] This embodiment provides a high-precision soil property mapping method that optimizes the number of samples by integrating air-to-ground imaging hyperspectral technology. The method includes acquiring UAV hyperspectral images of the area with a spatial resolution of 5 cm and performing preprocessing, including geometric correction, radiometric calibration, and reflectance conversion. On the day of UAV image acquisition, training and validation sampling points are deployed within the area, and sampling points are collected at a depth of 1 m. 2 Soil samples were collected from the sample points, and after preprocessing, the soil samples were placed in petri dishes to obtain ground-based hyperspectral imaging data. Based on the ground-based hyperspectral imaging data, a subset of one-quarter of the sample points in the training set was selected using the conditional Latin hypercube sampling method. The available phosphorus content of the soil in this optimized subset was measured, and a soil available phosphorus prediction model was trained based on the ground-based hyperspectral imaging data and the corresponding available phosphorus content of this optimized subset. Image classification was performed on the UAV hyperspectral images, and the soil area UAV hyperspectral images were segmented using vegetation index and maximum inter-class variance method. A 1m diameter was extracted from each sample point. 2 The spectral data within the specified range is used as the UAV soil spectral reflectance; the ground imaging hyperspectral data and the UAV soil spectral reflectance are band-aligned, and a band correction coefficient is calculated based on the ground imaging hyperspectral data to correct the UAV soil spectral reflectance, resulting in corrected UAV hyperspectral reflectance data; based on the corrected UAV hyperspectral reflectance data, for the missing parts, the surrounding 9m... 2 The average pixel value is filled to obtain a complete and continuous UAV hyperspectral image, which is then input into the soil available phosphorus prediction model for inversion to obtain the soil available phosphorus inversion distribution map of the region. Figure 1 The following is a flowchart of the method of the present invention, which will be described in detail below with reference to the embodiments.

[0051] The study area of the present embodiment occupies a land plot of about 20 mu as a sample area. The Resonon-Pika-L hyperspectral imager is mounted on a UAV platform to carry out aerial survey from 11:30 to 14:00 local time, and obtain the UAV hyperspectral image. After collecting the UAV hyperspectral image in the area, it is pre-processed, which includes geometric correction, radiation calibration and reflectivity conversion. After the flight operation is completed, the POS data (position and orientation system, POS) is imported into the SBGcenter software, and the POS post-difference-free technology is used to geometrically correct the single hyperspectral image obtained for each route. At the same time, according to the calibration file, the radiation calibration is carried out. Then, the data of each route is cut off to remove the turning data between routes. According to the orthophoto of the study area, the hyperspectral images on each route are geographically registered in ArcGIS 10.2 software, and then the multiple hyperspectral images of the study area are spliced together in ENVI 5.3. The reflectivity of the target cloth is used as a reference to convert the pixel gray value (DN) in the image to spectral reflectivity in the software megcube. The merged channel is used, the band interval is set to 2nm to reduce the influence of noise, and the wavelength is about 400-1000nm.

[0052] On the day of UAV image acquisition, the training set samples and the validation set samples are arranged in the area, and 5 surface soil samples at different positions are collected within 1m 2 around each sample point using the quincunx sampling method, and the soil samples of the sample points are pretreated and placed in a culture dish to obtain the ground imaging hyperspectral data. Grid sampling combined with random sampling is used to arrange the samples in the study area, i.e. a regular 50m*50m square is superimposed on the study area, and the center point of the square is taken as a sample, a total of 72, as a training set sample, and 20 samples in the study area are selected by random sampling method, each sample is spaced more than 25m apart, as a validation set sample, a total of 92 samples.

[0053] The soil samples are collected on the bare soil surface, and according to the planned sampling position, the soil surface is arranged flat using a sampling shovel, avoiding stones, crops and other debris, and then 1m 2Within the scope, 5 surface soil samples at different positions were collected by using the quincunx sampling method. All the samples were put into sampling bags, and the actual latitude and longitude coordinates at the time of sampling were recorded on the sampling bags using a carrier phase dynamic real-time differential instrument (Real-time kinematic, RTK), with a calibration error of not more than 5 cm. After the soil samples were taken back to the laboratory, the wet soil samples were put into culture dishes for ground imaging hyperspectral scanning, and the following treatment was carried out: first, the samples were homogenized, the cohesive soil blocks were scattered, and the large particle aggregates were removed to avoid significant particle size differentiation. Then, the treated wet soil was put into a standard culture dish, and vertical pressure was applied for moderate compaction, so that the sample formed a relatively consistent accumulation state and thickness in the culture dish, reducing the large pores and loose areas between particles, thereby reducing the spectral response noise caused by structural differences. Finally, a flat scraper was used to level the surface of the culture dish, trying to eliminate surface cracks, pits and obvious textures, and avoid "smearing effect" or "mirrorization" phenomenon; the leveled soil surface should meet the collection requirements of "small height difference, low roughness and continuous uniformity" to adapt to the pixel-level reflectance acquisition of ground imaging hyperspectral. After completion, it was marked according to the point number, and the remaining soil sample was dried, ground, sieved and then measured by chemical analysis method to determine the available phosphorus content. According to the standard of "NY / T1121.1-2006", the available phosphorus content of each soil sample was determined by ultraviolet / visible spectrophotometer method.

[0054] The ground imaging hyperspectral data was obtained by using a visible-near infrared HPPA (Hyperimager Plant Phenomics Analysis, LICA, Beijing, China) spectrometer, which integrates a variety of advanced technologies such as hyperspectral imaging analysis, RGB true color imaging, wireless automatic control and linear uniform illumination system. Before measurement, a series of pretreatment work was carried out, including zero adjustment, exposure and focusing, and interference from external light sources was avoided to ensure the accuracy and reliability of the data. Then, using the automatic focusing system and the wireless controlled multi-directional moving platform, the imaging hyperspectral of each culture dish with soil sample was collected, and the wavelength of the obtained spectral data was 401.3 to 1000.7 nm, the spectral resolution was 2.7 nm, and there were a total of 223 bands. Finally, the collected hyperspectral images were subjected to black and white correction, absorption test and image cropping, and then the average spectrum of each culture dish was taken as the ground imaging hyperspectral data of the soil sample.

[0055] Based on the ground imaging hyperspectral data, a quarter of the sample points are screened out from the training set sample points by using a conditional Latin hypercube sampling method to obtain an optimized sample point subset, the soil available phosphorus content of the optimized sample point subset is determined, and based on the ground imaging hyperspectral data and the corresponding soil available phosphorus content of the optimized sample point subset, a soil available phosphorus prediction model is trained. Through conditional Latin hypercube sampling, the value range of the reflectivity of each wave band of the sample spectrum is equally divided (divided according to the cumulative distribution function quantile), and the stratification in each spectral dimension is forced to be covered by the sample points. The specific steps are as follows:

[0056] ① Construct a hyperspectral data standardization covariant matrix of the training set sample points.

[0057] Input N candidate points p-dimensional hyperspectral data X ∈ RN×p (for example, N is 72 points, and p is 223 wave bands), and perform z-score standardization on each column, so as to eliminate the dimensional difference, ensure equal weight of each wave band, and avoid that high-variance wave band reflectivity dominates sampling screening. The formula is as follows:

[0058]

[0059] In the formula, z is the standardized value; x is the original wave band reflectivity; μ is the mean of the data set; and σ is the standard deviation of the data set.

[0060] ② Stratification cutting and Latin hypercube initialization are performed on each spectral wave band;

[0061] Firstly, for each wave band j, the cumulative distribution quantile is cut into k layers (k=n, and the target sample amount is, for example, 20 or 36), and the quantile of the jth variable is calculated:

[0062]

[0063] Wherein, Fj is the cumulative distribution function (CDF) of the jth variable, thereby generating an n×p matrix L, each column of which is a random arrangement of {1, 2,..., k}, avoiding sample aggregation caused by wave band correlation, and in addition, each value range layer of each wave band is forced to have a sample point, ensuring that each row and each column value is unique, satisfying the Latin hypercube constraint:

[0064]

[0065] ③ The optimized sample point subset is screened out by iteratively optimizing the objective function, and the objective function includes a distribution matching item for minimizing the spectral distribution difference between the subset and the full sample set, and a geographical dispersion penalty item for promoting the uniformity of the spatial distribution of the sample points.

[0066] Firstly, the objective function is defined as follows:

[0067]

[0068] where,

[0069] Distribution matching term, to minimize the marginal distribution difference between subset S and full sample set X. KS j Kolmogorov-Smirnov statistic, to measure the maximum vertical distance between the cumulative distribution function (CDF) of subset S and full sample set X in the jth band, to ensure the reflectance distribution pattern of subset S in key bands is consistent with that of full sample set X, avoiding missing caused by underestimating high variability areas. The specific formula is as follows:

[0070]

[0071] w j Band weight, the formula is as follows:

[0072]

[0073]

[0074] where VIP j is the variable importance projection of the jth band, to improve the decision weight of the feature band in distribution matching; p is the total number of bands (223 in this example); A is the total number of latent variables; SSY a is the variance explained by the ath latent variable; w ja is the normalized weight of the jth band in the ath latent variable;

[0075] Geographical dispersion penalty term, to suppress the aggregation of sample points in location distribution, to ensure uniform spatial coverage. represents the geographical coordinates of sample i; is the Euclidean distance (unit: meter) between sample i and sample j; y represents the decay radius of spatial correlation, which is 10% of the diagonal length of the field plot in this scheme; λ is a control function, to balance the weight of distribution matching and spatial dispersion; too high will cause the spectral distribution of subset S to deviate from that of full sample set X, and too low will cause uneven spatial distribution of sample points; represents the normalization factor, to eliminate the influence of sample size on the dimension of penalty value; is explained as: when d < y, the function value quickly approaches 1, to punish aggregated samples; when d > 3y, it approaches 0, to exempt from punishment.

[0076] Then replace one sample in the current subset S at random, calculate the degree of change of the objective function , according to the Metropolis criterion to accept the replacement (if <0, accept S^`; if >0, accept S^` with , where​ k is the iteration number, std(X) is the standard deviation of the original hyperspectral data as the initial temperature; 0.95 is the cooling coefficient, is the annealing temperature of the Kth iteration), and is iterated 10000 times. In this process, since the spectral edge points are crucial to representing the boundary of soil variation, an edge point priority replacement strategy is adopted, and the replacement probability is calculated by calculating the Mahalanobis distance MDi of the candidate points After the above operation, the n sample indexes and their geographic coordinates can be output.

[0077] ④Based on the JS divergence to measure the spectral distribution similarity of the sample subset and the full sample set, and when the divergence value is less than the preset threshold, it is determined that the subset can replace the full sample set for modeling.

[0078] The JS divergence is used to measure the distribution similarity of the full sample set and the subset, and the calculation formula is as follows:

[0079]

[0080] When DJS<0.1, the subset can replace the full sample set for modeling (in this embodiment, DJS=0.08<0.1 is calculated).

[0081] According to the screening of 36 and 20 sample points as new training sample sets, respectively, plus the training sample set of the complete 72 sample points without screening, a total of three data sets, respectively, using partial least squares regression PLSR for training, using the three data sets as the training set, and the 20 points as the validation set for accuracy comparison.

[0082] Regarding soil property inversion mapping, in addition to the prediction mapping method based on covariates, spatial interpolation technology can also be used. Ordinary Kriging (OK algorithm) is a commonly used spatial interpolation technology, and its prediction interpolation can be regarded as a regionalized variable, which is composed of a linear weighted sum of observation values. Ordinary Kriging interpolation is performed on all 72 sample points, and the specific form is as follows:

[0083]

[0084] In the formula, is the predicted value at the interpolation point; is the observation value weight assigned to the ith sample point, which is obtained by solving the Kriging equation set based on the semivariogram function; is the observation value of the ith sample point.

[0085] For the above prediction, the determination coefficient R 2 , the root mean square error RMSE and the prediction bias ratio RPD are used as the accuracy comparison standards, and the higher the model accuracy, the higher the R 2The greater the value of RPD is, the smaller the value of RMSE is, and the greater the value of RPD is.When RPD<1.4, it indicates that the model cannot predict the sample; when 1.4<RPD<2, it indicates that the model can roughly estimate the sample, and the prediction ability of the model can be improved by improving the modeling method; and when RPD>2, it indicates that the model has excellent prediction ability.

[0086] The soil property high-precision mapping method for optimizing the number of samples by fusing space-ground imaging hyperspectral technology shows the influence of different covariant groups on the prediction accuracy of soil available phosphorus, uses 72 sample points (scheme two), 36 sample points screened by conditional Latin hypercube sampling (scheme three), and 20 sample points screened by conditional Latin hypercube sampling (scheme four) as the prediction mode of covariant and the prediction mode based on ordinary Kriging interpolation (scheme one), and a total of four prediction schemes, and the results are shown in Table 1.

[0087] The application reveals the optimization mechanism of conditional Latin hypercube sampling and ground imaging hyperspectral for soil available phosphorus prediction by evaluating the prediction performance of partial least squares regression under different data sets. As can be seen from the table results, even if the number of sample points of the training set is reduced, the prediction accuracy based on ground imaging hyperspectral can still reach the level comparable to the full sample. When the sample amount is reduced to 36, the determination coefficient and the root mean square error of the validation set obtained by using the covariant set are 0.61 and 21.9 mg / kg, respectively, which is comparable to the accuracy (R 2 =0.61, RMSE=21.71 mg / kg) under the sampling density of 50 meters grid. 2 And when the sample amount is reduced to 20, the R

[0088] Table 1 Comparison of soil available phosphorus prediction accuracy of different prediction schemes

[0089]

[0090] Compared with the algorithm, the prediction performance of the PLSR modeling based on the covariant is better than that of the scheme of ordinary Kriging interpolation. This is because the OK algorithm only relies on the spatial autocorrelation of the target variable, and does not integrate external information; and the PLSR converts the non-spatial information such as chemical composition and physical structure contained in the spectrum into a prediction basis by introducing the spectral covariant, and expands the information source. At the same time, the PLSR solves the problem of high-dimensional collinearity of the spectrum through principal component dimension reduction, and the OK algorithm has no such mechanism. In addition, the OK algorithm can only capture spatial variation, while the PLSR can indirectly capture non-spatial driving factors through the spectrum, and the dimension reduction makes it less dependent on sample size and more robust.

[0091] The soil available phosphorus prediction model obtained in scheme four is saved for application to subsequent inversion mapping of soil available phosphorus based on unmanned aerial vehicle hyperspectral.

[0092] The unmanned aerial vehicle hyperspectral image is subjected to image classification, the soil region unmanned aerial vehicle hyperspectral image is segmented by combining a vegetation index and a maximum inter-class variance method, and the spectral reflectivity of the soil region is extracted as the unmanned aerial vehicle soil spectral reflectivity. When the unmanned aerial vehicle hyperspectral image is collected, most of the soil surface is exposed, and only a small amount of soybean plants and crop residues at the seedling stage are present. After radiation calibration and geometric correction, a hyperspectral true color image and a spectral diagram of typical ground object points extracted by visual interpretation are obtained. It can be found that the soil spectrum and the straw spectrum are very similar in waveform, and the reflectivity is also small, and it is difficult to find a certain band threshold standard to distinguish between the two. However, the vegetation spectrum and the soil spectrum are obviously distinguished in certain specific band ranges, there are chlorophyll absorption peaks at 450 nm and 650 nm, and there is a reflection peak at 550 nm, and they are obviously distinguished. Therefore, the normalized difference vegetation index (NDVI) and the differential vegetation index (DVI) can be used to distinguish vegetation and soil by using the characteristics of vegetation in the red and near-infrared wavelengths, and the calculation formula is as follows:

[0093]

[0094]

[0095] In the formula, NDVI represents the near-infrared band, and the band with a wavelength of 650 nm of the hyperspectral image is used in the calculation, represents the red light band, and the band with a wavelength of 900 nm of the hyperspectral image is used in the calculation. The calculation of the two indexes is completed in ENVI 5.3 using the band math tool. After the calculation, two gray-scale images of (a) and (c) are obtained, and the pixel values on the images are the single NDVI and DVI values, respectively. Figure 2

[0096] Figure 2 The results of segmenting the vegetation and soil parts are shown in the figures, wherein (a) is the calculation result of the normalized difference vegetation index NDVI, (c) is the calculation result of the differential vegetation index DVI, (b) is the binary image after threshold segmentation of (a), and (d) is the binary image after threshold segmentation of (c).

[0097] According to the gray-scale characteristics of the image, the threshold T is set, when the pixel value in the original image is greater than or equal to the threshold, it is considered to be vegetation, and the pixel is assigned a value of 1, which is displayed as white in the segmentation result image; and when it is less than or equal to the threshold, it is considered to be soil, and the pixel is assigned a value of 0, which is displayed as black in the segmentation result image, as shown in the following formula:

[0098] ​​

[0099] The maximum inter-class variance method traverses the threshold value from the minimum gray value to the maximum gray value according to the gray characteristics of the image, and at a certain gray value, the image is divided into two parts of background and foreground. When the inter-class variance of the foreground and the background is the largest, the probability of misclassification of the two parts of the background and the foreground is the smallest. The threshold value corresponding to the maximum inter-class variance is the adaptive segmentation threshold value. The threshold value left and right is divided into two groups by using the maximum inter-class variance method, and the gray value variance of the two groups is used to measure the difference between the two groups. When the threshold value of the gray value variance of the two groups reaches the maximum, the threshold value is the best threshold value. The (a) and (c) in the (a) and (c) are input into matlab R2016a respectively, and the threshold values obtained are 0.3 and 0.1 respectively. After the segmentation is completed, the visual interpretation method is used to compare the true color image near the sample point with the result after threshold segmentation, and the spectral reflectance curve near the sample point is viewed, the standard curve of soil and vegetation is compared, and whether the classification on the sample point is correct is checked, the segmentation result is evaluated, and the effect diagram after threshold segmentation is seen in the (b) and (d) of the binary image. After comparison, it is found that the segmentation effect of DVI is better because when identifying vegetation, NDVI has better effect when the vegetation coverage is relatively high, and DVI can identify the spectral characteristics of "small" vegetation when the vegetation coverage is relatively low. Figure 2 Figure 2

[0100] According to the binary image of DVI (d), by setting DVI as 0.1 as the segmentation threshold value, when DVI is greater than 0.1, it is vegetation, and when DVI is less than or equal to 0.1, it is soil. Figure 2

[0101] The soil pixel spectrum in a range of 1 m (20 pixels*20 pixels) is selected as the center of each sample point, the average value of the spectral reflectance is taken as the reflectance spectrum of a single sample point, and the whole image is resampled to 1 m spatial resolution.

[0102] ​​​The ground imaging hyperspectral data is band-aligned with the unmanned aerial vehicle soil spectral reflectance, and a band correction coefficient is calculated based on the ground imaging hyperspectral data to correct the unmanned aerial vehicle soil spectral reflectance, so that corrected unmanned aerial vehicle hyperspectral reflectance data are obtained. After the unmanned aerial vehicle soil hyperspectral reflectance is extracted, in order to solve the problem of spectral resolution difference between the ground imaging hyperspectral data and the unmanned aerial vehicle hyperspectral image data covering the entire study area, the extracted ground imaging hyperspectral data and the unmanned aerial vehicle hyperspectral data are aligned by using cubic spline interpolation and natural boundary method, so that the wavelength range and interval of the bands are consistent, the spectral interval is 2.7 nm, the band wavelength is 401.3 nm-1000.7 nm, and there are a total of 223 bands. For all the 92 collected sample points, since both the two kinds of spectra are hyperspectral, the number of bands is large, and therefore the traditional ratio mean method is not suitable for correcting all the bands into one correction coefficient. In the present application, the ratio of single band is used to correct each band of the unmanned aerial vehicle hyperspectral data by using the ground imaging hyperspectral data, and the correction coefficient of each band is calculated as follows: for each sample point and each band, the ratio of the ground imaging hyperspectral reflectance to the unmanned aerial vehicle hyperspectral reflectance is calculated; for each band, the ratio of all the sample points in the band is averaged, and the obtained average value is the correction coefficient of the band.

[0103]

[0104] In the formula, C λ is the correction coefficient of the band with a wavelength of λ; is the unmanned aerial vehicle hyperspectral reflectance of the band λ of the i th sample point; is the ground imaging hyperspectral reflectance of the band λ of the i th sample point; and n is the number of sample points (92 in the present embodiment).

[0105] After the correction coefficient C of each band is obtained, the reflectance of each band of the unmanned aerial vehicle hyperspectral data is adjusted, and the reflectance of each band is multiplied by the corresponding C, so that the corrected unmanned aerial vehicle hyperspectral reflectance data are obtained.

[0106] Based on the soil area unmanned aerial vehicle hyperspectral image and the corrected unmanned aerial vehicle hyperspectral reflectance data, information reconstruction is carried out, the entire image is set to an image with a resolution of 1 m, and then the soil available phosphorus prediction model is inputted for inversion, so that the soil available phosphorus inversion distribution map of the area is obtained. After threshold segmentation and image filling, the average value of the corrected unmanned aerial vehicle hyperspectral reflectance data in the surrounding 9m 2 (3m x 3m) of the missing part is used as the replacement value of the pixel, and finally the corrected soil unmanned aerial vehicle hyperspectral image range is filled to the entire study area, so as to ensure the continuity of the mapping.

[0107] The unmanned aerial vehicle hyperspectral image after the correction and filling is resampled to have a spatial resolution of 1 m, and is input into the soil available phosphorus prediction model to obtain the soil available phosphorus inversion distribution map of the research area. Figure 3 Specifically, first, the unmanned aerial vehicle hyperspectral image of the segmented soil area is subjected to band alignment and reflectance correction, then the missing part of the image is supplemented by the average value of the pixels in the surrounding 3 m x 3 m range, and finally the unmanned aerial vehicle hyperspectral data (each pixel contains spectral information of 223 bands in the wavelength range of 401.3 nm to 1000.7 nm with an interval of 2 nm) after the band alignment and reflectance correction covering the entire research area is input into the pre-trained and saved soil available phosphorus prediction model (in this embodiment, the model is trained based on the partial least squares regression PLSR algorithm). The model maps the spectral characteristics of each pixel to a predicted value of the soil available phosphorus content, and finally generates a digital raster map in which each pixel value represents the predicted concentration of the soil available phosphorus at the position, i.e. the soil available phosphorus inversion distribution map, by traversing all the pixels in the image. The distribution map can intuitively and quantitatively show the spatial distribution of the soil available phosphorus in the research area, and provide direct data support and decision basis for precision fertilization and farmland management.

[0108] The embodiment also provides a soil attribute high-precision mapping system fusing air-ground imaging hyperspectral technology to optimize the number of sampling, which realizes the soil attribute high-precision mapping method fusing air-ground imaging hyperspectral technology to optimize the number of sampling when the system is executed, and includes the following steps.

[0109] A data acquisition module is configured to acquire ground imaging hyperspectral data, unmanned aerial vehicle hyperspectral images, and the available phosphorus content of soil samples in a region.

[0110] A sample optimization and model training module is configured to perform conditional Latin hypercube sampling to screen and optimize a sample point subset, and train a soil available phosphorus prediction model based on the optimized sample point data.

[0111] An image processing and segmentation module is configured to process the unmanned aerial vehicle hyperspectral images, segment the soil area, and extract the unmanned aerial vehicle soil spectral reflectance.

[0112] A data correction module is configured to perform band alignment and reflectance correction on the ground imaging and unmanned aerial vehicle hyperspectral data.

[0113] An inversion mapping module is configured to generate a soil available phosphorus inversion distribution map by using the corrected data and the trained model.

[0114] In addition, the system also includes a data storage module configured to store the soil available phosphorus prediction model, correction coefficients, and intermediate processing data.

[0115] In summary, the present application overcomes the problems of high sampling cost and low precision in the prior art, and realizes fast and accurate monitoring of effective soil phosphorus, and is suitable for farmland scale investigation.

[0116] Embodiments of the application can be provided as methods, systems, or computer program products. Accordingly, the application can take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the application can take the form of a computer program product on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROMs, optical storage devices, etc.) embodying computer-readable program code.

[0117] The application is described in reference to the flow diagrams and / or block diagrams of the methods, apparatus (systems) and computer program products according to embodiments of the application. It should be understood that each flow and / or block in the flow diagrams and / or block diagrams, and combinations of flows and / or blocks in the flow diagrams and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general purpose computer, special purpose computer, embedded processing device, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions specified in the flow diagrams and / or block diagrams block or blocks. Figure 1 one or more flows and / or blocks Figure 1 means for carrying out the function specified by the flow or flows and / or block or blocks.

[0118] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instructions means which implement the function specified in the flow diagrams and / or block diagrams flow or flows and / or block or blocks. Figure 1 one or more flows and / or blocks Figure 1 means for carrying out the function specified by the flow or flows and / or block or blocks.

[0119] These computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flow diagrams and / or block diagrams flow or flows and / or block or blocks. Figure 1 one or more flows and / or blocks Figure 1 means for carrying out the function specified by the flow or flows and / or block or blocks.

[0120] The content not described in detail in the specification of the present application belongs to the prior art known to the person skilled in the art. It is indicated here that the above description helps the person skilled in the art to understand the present application, but does not limit the protection scope of the present application. Any implementation of equivalent replacement, modification, improvement and deletion of the above description without departing from the essential content of the present application falls within the protection scope of the present application.

Claims

1. A high-precision soil property mapping method that integrates air-to-ground imaging hyperspectral technology to optimize the number of samples, characterized in that, Includes the following steps: S1: Acquire UAV hyperspectral images with a spatial resolution of 5cm within the area and perform preprocessing, including geometric correction, radiometric calibration and reflectance conversion; S2: On the day of UAV image acquisition, training and validation sample points are set up in the area, and the sample points are collected at a depth of 1m. 2 Soil samples were collected and mixed thoroughly. After pretreatment, the soil samples from the sampling points were placed in a petri dish to obtain ground imaging hyperspectral data. S3: Based on the ground imaging hyperspectral data, a subset of one-quarter of the sample points in the training set is selected from the sample points using the conditional Latin hypercube sampling method. The available phosphorus content in the soil of the optimized sample point subset is measured. Based on the ground imaging hyperspectral data of the optimized sample point subset and the corresponding available phosphorus content in the soil, a soil available phosphorus prediction model is trained. S4: Perform image classification on the UAV hyperspectral image, segment the UAV hyperspectral image of the soil area by combining vegetation index and Otsu's method, and extract the spectral reflectance of the soil area as the UAV soil spectral reflectance, and resample it to 1m spatial resolution; the vegetation index includes normalized vegetation index and differential vegetation index. The maximum inter-class variance method adaptively determines the segmentation threshold based on the grayscale characteristics of the image. S5: Align the ground imaging hyperspectral data with the UAV soil spectral reflectance by band, and calculate the band correction coefficient based on the ground imaging hyperspectral data to correct the UAV soil spectral reflectance, thereby obtaining the corrected UAV hyperspectral reflectance data; the band alignment is performed using cubic spline interpolation and the natural boundary method to ensure that the ground imaging hyperspectral data and the UAV soil spectral reflectance have the same band wavelength range and spacing; S6: Based on the UAV hyperspectral image of the soil area and the corrected UAV hyperspectral reflectance data, image supplementation is performed to fill the entire image into a complete, continuous, and corrected UAV hyperspectral image, which is then input into the soil available phosphorus prediction model for inversion to obtain the soil available phosphorus inversion distribution map of the area.

2. The method for high-precision mapping of soil properties by optimizing the number of samples using air-to-ground imaging hyperspectral technology according to claim 1, characterized in that, In S2: The soil samples were collected with the sampling point 1m in diameter. 2 Within the specified area, surface soil samples were collected from five different locations using the plum blossom sampling method and mixed evenly. The training set sampling points are arranged using a 50m×50m grid sampling method; The verification sampling points are set up using a random sampling method; The process of placing the sample into the petri dish involves: homogenizing the sample to break up the clumps and agglomerates; placing the treated sample into a standard petri dish and applying vertical pressure to compact it appropriately; and using a straight scraper to level the surface of the petri dish.

3. The method for high-precision mapping of soil properties by optimizing the number of samples using air-to-ground imaging hyperspectral technology according to claim 1, characterized in that, The conditional Latin hypercube sampling method described in S3 includes: Construct a standardized covariate matrix for the hyperspectral data of wet soil samples from the training set; Each spectral band is segmented and initialized with a Latin hypercube; The optimized sample subset is selected by iteratively optimizing the objective function, which includes a distribution matching term to minimize the difference in spectral distribution between the subset and the full sample set, and a geographic dispersion penalty term to promote the uniformity of spatial distribution of the sample points. The similarity of the spectral distribution of a subset of samples to the full sample set is measured by JS divergence. When the divergence value is less than a preset threshold, the subset is determined to replace the full sample set for modeling.

4. The method for high-precision mapping of soil properties by optimizing the number of samples using air-to-ground imaging hyperspectral technology according to claim 3, characterized in that: The distribution matching term uses the Kolmogorov-Smirnov statistic to measure the cumulative distribution difference between the subset and the full sample set in each band. The geographic dispersion penalty term is calculated based on the Euclidean distance between samples.

5. The method for high-precision mapping of soil properties by optimizing the number of samples using integrated air-to-ground imaging hyperspectral technology according to claim 1, characterized in that, The calculation of the band correction coefficients described in S5 includes: For each sample point and each band, calculate the ratio of its ground imaging hyperspectral reflectance to the UAV hyperspectral reflectance; For each band, the ratio of all samples in that band is averaged, and the average value is the correction coefficient for that band.

6. The method for high-precision mapping of soil properties by optimizing the number of samples using air-to-ground imaging hyperspectral technology according to claim 1, characterized in that, The image supplementation described in S6 is based on the missing areas in the corrected UAV hyperspectral image, utilizing the surrounding 9m area. 2 The average value of the corrected UAV hyperspectral reflectance data within the range is used as a substitute.

7. A high-precision soil property mapping system that integrates air-to-ground imaging hyperspectral technology to optimize the number of samples, characterized in that, When executed, this system implements the high-precision soil property mapping method according to any one of claims 1-6, which optimizes the number of samples by fusing air-to-ground imaging hyperspectral technology, including: The data acquisition module is used to acquire UAV hyperspectral images, ground imaging hyperspectral data, and available phosphorus content of soil samples within the region; The sample optimization and model training module is used to perform conditional Latin hypercube sampling to screen and optimize the sample subset, and train a soil available phosphorus prediction model based on the optimized sample data. The image processing and segmentation module is used to process UAV hyperspectral images, segment soil regions, and extract the spectral reflectance of UAV soil. The data correction module is used to perform band alignment and reflectance correction on ground imaging and UAV hyperspectral data; The inversion mapping module is used to generate soil available phosphorus inversion distribution maps using corrected data and trained models.

8. A high-precision soil property mapping system that optimizes the number of samples by integrating air-to-ground imaging hyperspectral technology according to claim 7, characterized in that, The system also includes a data storage module for storing the soil available phosphorus prediction model, correction coefficients, and intermediate processing data.

Images