CN116503747A - Heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing - Google Patents

Heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing Download PDF

Info

Publication number
CN116503747A
CN116503747A CN202310793690.XA CN202310793690A CN116503747A CN 116503747 A CN116503747 A CN 116503747A CN 202310793690 A CN202310793690 A CN 202310793690A CN 116503747 A CN116503747 A CN 116503747A
Authority
CN
China
Prior art keywords
sentinel
2lai
data
reflectivity
remote sensing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202310793690.XA
Other languages
Chinese (zh)
Other versions
CN116503747B (en
Inventor
徐保东
金文捷
胡琼
王聪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Huazhong Agricultural University
Original Assignee
Huazhong Agricultural University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Huazhong Agricultural University filed Critical Huazhong Agricultural University
Priority to CN202310793690.XA priority Critical patent/CN116503747B/en
Publication of CN116503747A publication Critical patent/CN116503747A/en
Application granted granted Critical
Publication of CN116503747B publication Critical patent/CN116503747B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00Scenes; Scene-specific elements
    • G06V20/10Terrestrial scenes
    • G06V20/188Vegetation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06NCOMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00Machine learning
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06NCOMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00Computing arrangements based on biological models
    • G06N3/02Neural networks
    • G06N3/08Learning methods
    • G06N3/09Supervised learning
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/40Extraction of image or video features
    • G06V10/52Scale-space analysis, e.g. wavelet analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/70Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/77Processing image or video features in feature spaces; using data integration or data reduction, e.g. principal component analysis [PCA] or independent component analysis [ICA] or self-organising maps [SOM]; Blind source separation
    • G06V10/774Generating sets of training patterns; Bootstrap methods, e.g. bagging or boosting
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/70Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/82Arrangements for image or video recognition or understanding using pattern recognition or machine learning using neural networks
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A40/00Adaptation technologies in agriculture, forestry, livestock or agroalimentary production
    • Y02A40/10Adaptation technologies in agriculture, forestry, livestock or agroalimentary production in agriculture
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Software Systems (AREA)
  • General Health & Medical Sciences (AREA)
  • Multimedia (AREA)
  • Artificial Intelligence (AREA)
  • Computing Systems (AREA)
  • Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Databases & Information Systems (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Molecular Biology (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

The invention discloses heterogeneous earth surface leaf area index inversion based on multi-scale remote sensing data, which downloads satellite images of remote sensing satellites with different scales, such as ten-meter-level remote sensing images and hundred-meter-level remote sensing images; obtaining LAI based on the ten-meter remote sensing image; the ten-meter-level LAI result is cooperated with the reflectivity of the hundred-meter-level remote sensing image to construct an LAI-reflection data set; sample selection is carried out based on a sampling strategy of priori knowledge; and (3) carrying out inversion model optimization based on various machine learning models, and obtaining the leaf area index of the heterogeneous earth surface with higher precision. According to the invention, the information of the ten-meter-level remote sensing image and the information of the hundred-meter-level remote sensing image are fully coordinated through a sample selection method based on priori knowledge and a simulated annealing algorithm, so that the inversion accuracy of the heterogeneous earth surface leaf area index is improved.

Description

Heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing
Technical Field
The invention relates to the technical field of inversion of large-area heterogeneous earth surface leaf area indexes, in particular to a heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing, and belongs to the technical field of quantitative remote sensing.
Background
Leaf Area Index (LAI) is an important parameter that characterizes vegetation growth and is defined as one half of the total green leaf area per unit ground area. As is well known, LAI is a basic climate variable (ECV) approved by the Global Climate Observation System (GCOS) and plays an important role in model simulation of the terrestrial ecosystem. Many studies have focused on applications of various LAIs, such as vegetation dynamic monitoring and crop yield estimation. Thus, the generation of a long-term LAI data set worldwide is a prerequisite for capturing the spatiotemporal variation of vegetation growth caused by human activity or climate change. High frequency earth observation is critical to capturing the critical phase of vegetation growth dynamics. Several existing hundred meter resolution (> 100 meters) surface reflectivity products with high temporal resolution (.ltoreq.1 day), including MODIS, VEGETATION, AVHRR, PROBA-V and VIIRS, have been widely used to generate global scale LAI time series. However, the uncertainty of LAI products with hundred meter resolution is generally greater than 0.6, and cannot meet the uncertainty requirements of the GCOS-suggested LAI dataset. Therefore, how to alleviate the surface heterogeneity image to improve the accuracy of the LAI product with hundred meter resolution is a difficulty.
Traditionally, LAI is obtained by field investigation, such as destructive sampling, leaf clamps, and plant canopy analysis. These methods are perhaps the most accurate way to measure LAI in regular field surveys. However, they are time consuming, costly, and collecting a temporally continuous LAI over a large spatial area is not practical. Remote sensing technology provides a fast, reproducible, synchronized and cost-effective method for quantifying large-area biophysical variables, and is widely used to generate global-scale LAI datasets. Currently, some satellites, such as Landsat-5/7/8/9, sentinel-2 and Gaofen-1/6, may provide remote sensing images of non-metered resolution, making it possible to introduce sub-pixel information in a hundred meter resolution grid of pixels. Among them, sentinel-2A emitted 6 months in 2015, provides valuable non-metered resolution information for a number of land applications including crop monitoring and biophysical estimation. The sentinel-2 task can provide multispectral images with high spatial (10-60 meters) and temporal (5 days or less) resolution. In addition, a simplified secondary prototype processor (SL 2P) implemented in the open-source sentinel application platform (SNAP) was developed for the sentinel-2 MSI image for quantifying biophysical variables. Thus, the LAI of Sentinel-2 can be easily generated by the SL2P tool and shown to be highly consistent with ground LAI measurements of different scales in multiple verification activities. However, sentinel-2LAI cannot meet the dynamic monitoring of vegetation growth over a wide range due to the effects of cloud pollution. Therefore, how to integrate sub-pel information provided by the Sentinel-2LAI to improve LAI retrieval with hundred meters resolution remains a challenge to be solved to maximize the role of the high frequency LAI dataset.
Disclosure of Invention
Aiming at the problem of low precision of a long-time-sequence hundred-meter-level LAI product in the prior art, combining the characteristics of the surface heterogeneity image of the hundred-meter-level reflectivity, cooperatively using multi-scale remote sensing data, the invention provides a heterogeneous surface leaf area index inversion method based on multi-scale remote sensing, wherein the ten-meter-level reflectivity image is aggregated to the hundred-meter-level scale according to the upscale rule of 'average before inversion', a LAI-reflection source data set is constructed and obtained, a sample is selected by utilizing a sampling strategy of priori knowledge, and a machine learning inversion model is optimized, so that the heterogeneous surface hundred-meter-level LAI inversion is realized.
The above object of the present invention is achieved by the following technical means:
a heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing comprises the following steps:
step one, acquiring Sentinel-2L 2A data of a research area, MOD09A1 reflectivity data of MODIS, MCD12Q1 land coverage data and high-resolution reference map data of ImagineS;
generating Sentinel-2LAI data based on the Sentinel-2L 2A data;
step three, obtaining an actual Sentinel-2LAI aggregation value based on the Sentinel-2LAI data, wherein the actual Sentinel-2LAI aggregation value, the matched MOD09A1 reflectivity data and the MCD12Q1 land coverage data form a sample point, and each sample point forms an LAI-reflection data set;
step four, selecting representative sample points of the LAI-reflection data set based on priori knowledge;
training the prediction model by using representative sample points to obtain a trained prediction model, inputting MOD09A1 reflectivity data to be predicted to the prediction model to obtain a predicted Sentinel-2LAI aggregation value, and inverting the predicted Sentinel-2LAI aggregation value to obtain a leaf area index of the heterogeneous earth surface.
The actual Sentinel-2LAI aggregate value in step three as described above is obtained based on the following steps:
converting the Sentinel-2LAI data into Sinusidal projection, cutting the Sentinel-2LAI data converted into Sinusidal projection according to the range of a research area, and aggregating the cut Sentinel-2LAI data onto a MODIS pixel grid by using a simple average method to obtain an actual Sentinel-2LAI aggregation value.
The aggregation of the cut Sentinel-2LAI data onto the MODIS pixel grid by the simple averaging method specifically comprises the following steps:
the cut Sentinel-2LAI data is corresponding to the interior of each MODIS pixel grid to average,
if each of the Sentinel-2LAI sub-pel values is between 0 and 10, the Sentinel-2LAI sub-pel value is determined to be a valid value, otherwise the Sentinel-2LAI sub-pel value is determined to be an invalid value,
calculating the effective value proportion of the Sentinel-2LAI sub-pixels in each MODIS pixel grid, if the effective value proportion is more than or equal to 90 percent, calculating the average value of the Sentinel-2LAI sub-pixels in the MODIS pixel grid to obtain the actual Sentinel-2LAI aggregate value,
and (3) obtaining the reflectivity of the MODIS pixels corresponding to the MODIS pixel grid through indexing, and if the effective value proportion is smaller than 90%, masking the pixels corresponding to the MODIS pixel grid with the resolution of 500 meters.
Step four as described above comprises the steps of:
firstly, extracting a set number of sample points to be determined as sample points, and calculating an objective function of the sample points to be determinedOF
Secondly, iteratively replacing a plurality of sample points to be determined, comparing objective functions of the sample points to be determined before and after optimization, accepting the iterative replacement of the sample points to be determined if the objective function after the iterative replacement is smaller than the objective function before the iterative replacement, otherwise refusing to discard the iterative replacement of the sample points to be determined last time;
and finally, outputting each sample point to be determined as a representative sample point if the objective function is smaller than the threshold value or reaches the preset iteration times.
Objective function as described aboveOFBased on the following formula:
wherein ,is the average deviation of reflectivity, +.>Is vegetation type average deviation, ++>Is the nearest neighbor index, i1 is the number of the reflectivity band in the MOD09A1 reflectivity data,ais the number of reflectivity bands in MOD09A1 reflectivity data.
Average deviation of reflectivity as described aboveBased on the following formula:
where n is the number of samples to be validated,index of reflectance-> and />Is the minimum and maximum reflectivity of the ith reflectivity interval, +.>Indicating that the reflectivity of a certain type of reflectivity Band is greater than or equal to +.>And is less than->Is determined by the number of sample points to be confirmed,
average deviation of vegetation typesBased on the following formula: />
wherein ,mrefers to the total number of vegetation types for the sample points to be validated,is the number of sample points to be confirmed belonging to the type t vegetation,k t is the duty cycle of sample points belonging to the type t vegetation in the LAI-reflection dataset,
the nearest neighbor indexNNIBased on the following formula:
wherein min is%d ij ) Representing the geographical location distance between the sample point i to be confirmed and the nearest sample point j to be confirmed,nis the number of sample points to be validated,Ais the area of the investigation region.
The prediction model is one or more of an artificial neural network prediction model, a random forest regression model, a Gaussian process regression model and a support vector machine regression model;
when the prediction models are multiple, constructing a test set with MOD09A1 reflectivity data and actual Sentinel-2LAI aggregation values as sample points, inputting the MOD09A1 reflectivity data in the test set into the prediction model which is completed by training, outputting corresponding predicted Sentinel-2LAI aggregation values, and calculating root mean square error RMSE, average deviation Bias and decision coefficient R between the actual Sentinel-2LAI aggregation values and the predicted Sentinel-2LAI aggregation values of each prediction model according to the actual Sentinel-2LAI aggregation values and the corresponding predicted Sentinel-2LAI aggregation values in the test set 2 Selecting RMSE to be less than or equal to 0.5, and R 2 A prediction model with a Bias of 0.8 or more and 0.01 or less is used as the optimal prediction model.
Compared with the prior art, the invention has the following beneficial effects:
the invention aims to solve the problem of low inversion precision of the long-time sequence hundred-meter-level leaf area index, and cooperates with a ten-meter-level remote sensing image and a hundred-meter-level remote sensing image to generate an LAI-reflection data set, and combines a sampling strategy based on priori knowledge, and a prediction model is optimized to perform leaf area index inversion of the heterogeneous earth surface. Finally, based on the cooperation of the multi-scale remote sensing data, the inversion precision of the hundred-meter leaf area index is improved through a preferable prediction model, and the earth surface heterogeneity effect in leaf area index inversion is relieved.
Drawings
FIG. 1 is a schematic flow chart of the method of the present invention.
Detailed Description
The invention will now be described in further detail with reference to the following examples, which are given by way of illustration and explanation only, and are not intended to be limiting, for the understanding and practice of the invention by those of ordinary skill in the art.
Examples
As shown in FIG. 1, the invention discloses a heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing. The method comprises the following steps:
step one, acquiring Level-1C-Level Sentinel-2L 2A data of a Study area, MOD09A1 reflectivity data of MODIS, MCD12Q1 land coverage data and high-resolution reference map data of ImagineS. Study area was selected as the 200 km by 200 km area on the high score reference plots of kennia, spanish, france and ukraine, respectively. Specifically, the detailed operation steps of the first step are as follows:
image data of Sentinel-2l 1c with low cloud cover was downloaded from the eujec office. Because the TOA reflectivity is provided by the image data of the Sentinel-2L 1C, the Senn 2cor atmospheric correction plug-in provided by the European air office used in the invention performs atmospheric correction on the image data of the Sentinel-2L 1C to generate the Sentinel-2L 2A data. The Sen2Cor atmospheric correction plug-in provides cloud detection for the image and stores the results in the scene classification bands of the L2A data for pixel quality screening to remove the effects of cloud coverage, cloud shadows, and non-vegetation surface coverage.
The Terra MODIS 8 days downloaded from NASA website to synthesize MOD09A1 reflectance data with a resolution of 500 meters. The MOD09A1 reflectance data is additionally quality controlled using a quality control band (QA) built into the MOD09A1 product to remove the effects of snow, partial or complete cloud coverage, cloud shadows, or low pixel quality. This step outputs MOD09A1 reflectance data comprising six reflectance bands (red band: 620-670nm; near infrared band: 841-876nm; green band: 545-565nm; mid infrared band: 1230-1250nm,1628-1652nm,2105-2155 nm), the observation geometry for three solar viewing angles (solar zenith angle, viewing angle zenith angle, and relative azimuth angle).
Downloading the global annual three-level land coverage data of the MCD12Q1, selecting a third scheme of the MCD12Q1 product, namely a scheme based on MODIS leaf area index/photosynthetic effective radiation, and screening vegetation and non-vegetation type pixels in the MOD09A1 product to eliminate unnecessary errors introduced by non-vegetation types. This step outputs MCD12Q1 land cover data including eight vegetation types, respectively grassland, bush, broadleaf farmland, tropical grassland, evergreen broadleaf forest, deciduous broadleaf forest, evergreen conifer forest.
The high-score reference map within the Study area is downloaded from the ImagineS website, aimed at obtaining in-situ LAI measurements to verify LAI values obtained by satellite inversion.
Step two, generating Sentinel-2LAI data based on the Sentinel-2L 2A data
Based on the Sentinel-2L 2A data obtained in the step one, sentinel-2LAI data with 20m resolution is generated by using a SNAP toolbox. Specifically, the detailed operation steps of the second step are as follows:
inputting the Sentinel-2L 2A data obtained in the step one into SNAP software, and inverting by using a Land biological inversion algorithm (Sentinel-2 Land bio-physical processor, SL 2P) to obtain 20m resolution Sentinel-2LAI data. The SL2P algorithm is a biomass inversion algorithm developed by the European space agency for Sentinel-2 image based on a neural network. Specifically, the process of the algorithm includes three main steps: generating a training data set by means of a Radiation Transfer Model (RTM), i.e. PROSAIL, calibrating the neural network and predicting the biophysical variables using the trained neural network. Eight band information (B3, B4, B5, B6, B7, B8A, B11 and B12) and direction information (viewing zenith, solar zenith and relative azimuth) were considered in the PROSAIL model to simulate the corresponding canopy reflectivities. Training was performed for 8 band information and 3 direction information using an ANN neural network. And finally, predicting the Sentinel-2L 2A data by using the calibrated ANN network parameters to obtain the Sentinel-2LAI data.
Step three, obtaining an actual 500m Sentinel-2LAI aggregation value based on Sentinel-2LAI data, wherein the actual 500m Sentinel-2LAI aggregation value, the matched MOD09A1 reflectivity data and the MCD12Q1 land coverage data form a sample point, and each sample point forms an LAI-reflection source data set;
and (3) performing projection conversion on the Sentinel-2LAI data obtained in the step two. The Sentinel-2LAI data was a UTM-WGS84 projection, and to maintain consistency with the MOD09A1 reflectance data projection, the Sentinel-2LAI data was converted to a Sinussoidal projection. The Sentinel-2LAI data converted to Sinussoidal projections was cropped according to the extent of the Study area, study area. And aggregating the cut Sentinel-2LAI data onto a MODIS pixel grid with the resolution of 500 meters by using a simple average method to obtain an actual 500m Sentinel-2LAI aggregation value, and constructing a LAI-reflection data set obtained by matching the actual 500m Sentinel-2LAI aggregation value with MOD09A1 reflectivity data and MCD12Q1 land coverage data. The simple averaging method is based on a MODIS pixel grid with a resolution of 500 meters, and the cut Sentinel-2LAI data with a resolution of 20 meters is corresponding to the interior of each MODIS pixel grid with a resolution of 500 meters to average.
The aggregation rules are as follows: (1) If each 20m resolution Sentinel-2LAI sub-pixel value is between 0 and 10, the 20m resolution Sentinel-2LAI sub-pixel value is judged to be a valid value, otherwise, the 20m resolution Sentinel-2LAI sub-pixel value is judged to be an invalid value. (2) The effective value ratio of 20m Sentinel-2LAI sub-pixels in each 500m resolution MODIS pixel grid is calculated. If the effective value proportion is greater than or equal to 90%, calculating the average value of 20m Sentinel-2LAI sub-pixels in the MODIS pixel grid with the resolution of 500m by using a simple average method, obtaining an actual 500m Sentinel-2LAI aggregate value, and indexing to obtain the MODIS pixel reflectivity corresponding to the MODIS pixel grid with the resolution of 500 m. If the effective value proportion is smaller than 90%, masking the pixel corresponding to the MODIS pixel grid with the resolution of 500 meters.
Step four, selecting representative sample points of the LAI-reflection data set based on priori knowledge;
and step three, obtaining an LAI-reflection data set, wherein the LAI-reflection data set comprises a plurality of sample points, and each sample point consists of an actual 500m Sentinel-2LAI aggregation value, MOD09A1 reflectivity data matched with the actual 500m Sentinel-2LAI aggregation value and MCD12Q1 land coverage data. And carrying out sample selection based on priori knowledge on the LAI-reflection data set.
The sampling is performed by using a priori knowledge-based sampling strategy, and the method iteratively ensures uniform and discrete distribution of the sample in various characteristic attribute variables and geographic spaces through a simulated annealing algorithm. To achieve both sample acquisition efficiency and sample acquisition accuracy, a priori knowledge based on Study is proposedThe area is sampled. Two evaluation indexes are to be constructed to evaluate: (1) Calculating representative of MOD09A1 reflectance data and MCD12Q1 land cover data for selected sample pointsAnd) As an attribute space representative index; (2) For the representative evaluation of the geospatial position, calculating the nearest neighbor index of the sample pointNNI) As a geospatial representative index. Based on->、/>AndNNIcalculating an objective functionOF) And then the sample points are iteratively optimized by using a simulated annealing algorithm, and the specific method is as follows:
and iteratively optimizing the sample points through a simulated annealing algorithm. Firstly, sampling a certain Sample point to be used as a Sample point to be determined, and calculating an objective function of the Sample point to be determinedOF) The method comprises the steps of carrying out a first treatment on the surface of the Secondly, iteratively replacing a plurality of sample points to be determined, comparing objective functions of the sample points to be determined before and after optimization, accepting the iterative replacement of the sample points to be determined if the objective function after the iterative replacement is smaller than the objective function before the iterative replacement, otherwise refusing to discard the iterative replacement of the sample points to be determined last time; and finally, outputting each sample point to be determined as a representative sample point if the objective function is smaller than the threshold value or reaches the preset iteration times. The indexes are calculated as follows:
first, for the reflectance Band of MOD09A1 in the LAI-reflectance dataset, the reflectance corresponding to a certain type of reflectance Band of MOD09A1 of each (n) sample points to be confirmed is arranged from small to large, and is divided into n reflectance sections in sequence. And then calculating the frequency number of the reflectivity wave bands of all the sample points to be confirmed falling into each reflectivity interval. Theoretical assumption that the sample point to be confirmed should fall uniformly in each reflectance interval, i.e. eachEach reflectance interval should have a sample point to be confirmed. Calculating the average deviation of reflectivity based on the following formula
Where n is the number of samples to be validated,index of reflectance-> and />Is the minimum and maximum reflectivity of the ith reflectivity interval, +.>Indicating that the reflectivity of a certain type of reflectivity Band is greater than or equal to +.>And is less than->To be confirmed, and the number of sample points to be confirmed.
And calculating the number ratio of sample points to be confirmed, of which the LAI-reflection data set belongs to the type t vegetation type, aiming at the type t vegetation type of the MCD12Q1 land coverage data in the LAI-reflection data set. And then calculating the number of the duty ratios of the Sample points to be confirmed belonging to the t-th vegetation type of the Sample points to be confirmed in all the Sample points to be confirmed. Finally, calculating the average deviation between the frequency of the sample points to be confirmed falling into each reflectivity interval and the frequency of the sample points to be confirmed falling into each reflectivity interval in theory, namelyThe number and LA of the Sample points Sample to be confirmed belonging to the type t vegetation typesAverage deviation of vegetation types between the numbers belonging to the type t vegetation types in the I-reflection dataset +.>Given by the following relationship: />
wherein ,mrefers to the total number of vegetation types for the sample points to be validated,is the number of sample points to be confirmed belonging to the type t vegetation, +.>Is a sample point to be confirmed belonging to the type t vegetation,k t is the duty cycle of sample points belonging to the type t vegetation in the LAI-reflection dataset.
Nearest neighbor indexNNIThe calculation method of (2) is as follows:
wherein min is%d ij ) Representing the geographical location distance between each Sample point i to be confirmed and its nearest Sample point j to be confirmed in Sample points Sample to be confirmed,d ij is the geographical location distance between the sample point i to be confirmed and the sample point j to be confirmed,nis the number of Sample points to be validated,Ais the area of the Study area.
Objective function [ ]OF) Is characterized by that its attribute space is represented by index #, and its attribute space is represented by index #, its attribute space is represented by index # and />) And a geospatial representative index (>) The method is obtained by common calculation and concretely comprises the following steps: />
wherein ,is the average deviation of reflectivity, +.>Is vegetation type average deviation, ++> and />Are all representative indexes of attribute space, < >>Is the nearest neighbor index, is a geospatial representative index,athe number of the reflectivity bands in the MOD09A1 reflectivity data is 6 by default, i1 is the number of the reflectivity bands in the MOD09A1 reflectivity data, and +.>Is the sum of average deviations of vegetation types in each reflectance band in MOD09A1 reflectance data.
For the red band, the reflectivity of selected sample points within the investigation region is in the interval 0-0.3. For the near infrared band, the selected sample point reflectivities are all between 0-0.6. For the green band, the reflectivity of the selected sample points is between 0 and 0.2. In the sampling results of the three wave bands, the reflectivity difference between the sample point and the LAI-reflection data set is controlled within 5%. The four statistical indicators include MeanDev, stdDev, skDev and KuDev, which represent differences in mean, standard deviation, skewness, and kurtosis of the sample points and LAI-reflectance dataset, respectively. There was little difference from the LAI-reflectance dataset in the sample points selected by representative samples (MeanDeV < 0.005, stdDev < 0.005, skDev < 0.2, kuDev < 0.5).
Step five, inverting leaf area index by optimizing machine learning model
Modeling is carried out by utilizing the representative sample points selected in the step four, and one or more of four machine learning prediction models (artificial neural network prediction model, ANN, random forest regression model, RFR, gaussian process regression model, GPR, support vector machine regression model and SVR) are combined. MOD09A1 reflectance data (six reflectance bands: red band, near infrared band, green band and three mid-infrared bands; three observation geometries of solar viewing angles: solar zenith angle, viewing angle zenith angle and relative azimuth angle) in representative sample points were taken as input data, predicted 500m Sentinel-2LAI aggregate values were taken as output data, and the relationship between the input data and the output data was learned.
The representative sample point randomization obtained in the step four is divided into a training set (80%) of a prediction model and a verification set (20%) of the prediction model. And randomly dividing 20 times to generate 20 groups of training sets and 20 corresponding groups of verification sets, namely, the 20 groups of random training sets and the 20 corresponding groups of verification sets are independently modeled and verified and are used for testing the stability of different machine learning models. In order to verify whether the prediction model has generalization capability, 3 scenes of a research area are selectedAnd 5. Mu.l->And (5) detecting the MOD09A1 reflectivity data corresponding to the high-resolution reference diagram of the measured data. And aggregating the high-resolution reference map onto a 500m MODIS grid by using a simple average method, and performing pixel matching with an MOD09A1 product with cloud and non-vegetation type removed to obtain MOD09A1 reflectivity data and an actual 500m Sentinel-2LAI aggregation value as a test set of a prediction model.
Training of the predictive model is typically based on a training set and a validation set, by which the Root Mean Square Error (RMSE), mean deviation (Bias), and decision coefficient (R) between the actual 500m Sentinel-2LAI aggregate and the predicted 500m Sentinel-2LAI aggregate can be calculated 2 ) RMSE is typically greater than 0, with smaller values representing higher model accuracy. Bias is greater than 0, then the model overestimation is indicated; bias is smallAt 0, then the model underestimation is indicated; bias approaches 0, which means higher model accuracy. R is R 2 The more the value of (2) is generally between 0 and 1, the more it approaches 1, the higher the model accuracy. After performance evaluation by four machine learning, the RMSE is selected to be less than or equal to 0.5, and R 2 And (3) a model with Bias approaching 0 (optionally less than or equal to 0.01) is preferably output, the optimal prediction model obtained by fitting is preferably output, the predicted 500m sentel-2 LAI aggregation value can be obtained by inputting MOD09A1 reflectivity data to be predicted into the optimal prediction model, and further inversion is carried out on the predicted 500m sentel-2 LAI aggregation value to obtain the leaf area index of the heterogeneous earth surface. If not satisfied, the RMSE is less than or equal to 0.5, and R 2 And (3) if the ratio is greater than or equal to 0.8, and if Bias approaches 0, returning to the step one.
In this embodiment, four machine learning model parameters are set as follows:
the input layer of the ANN is set to 9 neurons representing 6 reflectance bands and 3 direction information bands, respectively. The hidden layer has 18 neurons, and the input layer is doubled, so that the performance of the neural network can be maximized. LAI is the only output variable of the output layer. A Levenberg-Marquardt learning algorithm with a squaring loss function is used to learn the nonlinear relationship between the input data and the output data.
In RFR, it is assumed that there are N sample points and M attribute variables. First, each independent decision tree obtains N sample points by a replaced sample as a sample at the root node of the decision tree. Secondly, each root node of the decision tree randomly selects mtry attribute variables (mtry < < M) as splitting attributes of the root node. In the decision tree forming process, each node is split by the last step until the decision tree can not be split any more. The result of random forest regression prediction is the average value of the ntree decision tree results, and mtry=m/3 is set, i.e. the node value is one third of the number of the attribute variables of the input dataset. Let ntree=500, i.e. the number of decision trees is 500.
The GPR, due to its kernel-based nature, makes it possible to learn the training set efficiently by maximizing the marginal likelihood. The GPR super-parameters are adjusted by adopting an automatic optimization super-parameter algorithm, and the algorithm finds the super-parameters which minimize the five-time cross-validation loss value by automatically adjusting the super-parameters.
The SVR principle is to solve for a separation hyperplane that can correctly divide the training dataset and has minimal geometric separation. By using automatic hyper-parameter optimization, hyper-parameters that minimize five-fold cross-validation losses are found.
For the fitting accuracy of the prediction model, RFR average accuracy (bias=0.000, R 2 =0.966 and rmse=0.187), better than ANN average accuracy (bias=0.000, R 2 =0.913 and rmse=0.287), GPR average precision (bias=0.002, R 2 =0.891 and rmse=0.321) and SVR average precision (bias=0.017, R 2 =0.868 and rmse=0.358).
For the validation accuracy of the prediction model, GPR average accuracy (bias= -0.002, r 2 =0.892 and rmse=0.326) and ANN average accuracy (bias=0.002, R 2 =0.903 and rmse=0.305), better than SVR average accuracy (bias= -0.017, R 2 = 0.874 and rmse=0.345) and RFR average accuracy (bias=0.003, R 2 =0.881 and rmse=0.328).
R according to 20 training sets and corresponding 20 verification sets 2 The distribution of RMSE and Bias, the stability of the model is evaluated through numerical fluctuation, wherein the accuracy stability of RFR and ANN is poor, the numerical fluctuation is large, the accuracy stability of GPR and SVR is good, and the numerical fluctuation is small.
The test set has 5836 pixels, and for the test precision of four prediction models, the best precision of the GPR is (bias= -0.308, R 2 =0.826 and rmse=0.593), better than SVR best precision (bias= -0.333, R 2 =0.824 and rmse=0.612), RFR best precision (bias= -0.366, R 2 =0.769 and rmse=0.687) and ANN best precision (bias= -0.316, R 2 = 0.793 and rmse=0.644). R according to test set 2 The distribution of RMSE and Bias, the stability of the model is evaluated through numerical fluctuation, wherein the stability of the accuracy of ANN is poor, the numerical fluctuation is large, the stability of the accuracy of GPR, SVR and RFR is good, and the numerical fluctuation is small. Fitting precision and verification precision of four models are comprehensively evaluatedThe accuracy of the GPR performs best, and is selected as the optimal machine learning model.
It should be noted that the specific embodiments described in this application are merely illustrative of the spirit of the invention. Those skilled in the art may make various modifications or additions to the described embodiments or substitutions thereof without departing from the spirit of the invention or its scope as defined in the accompanying claims.

Claims (7)

1. The heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing is characterized by comprising the following steps of:
step one, acquiring Sentinel-2L 2A data of a research area, MOD09A1 reflectivity data of MODIS, MCD12Q1 land coverage data and high-resolution reference map data of ImagineS;
generating Sentinel-2LAI data based on the Sentinel-2L 2A data;
step three, obtaining an actual Sentinel-2LAI aggregation value based on the Sentinel-2LAI data, wherein the actual Sentinel-2LAI aggregation value, the matched MOD09A1 reflectivity data and the MCD12Q1 land coverage data form a sample point, and each sample point forms an LAI-reflection data set;
step four, selecting representative sample points of the LAI-reflection data set based on priori knowledge;
training the prediction model by using representative sample points to obtain a trained prediction model, inputting MOD09A1 reflectivity data to be predicted to the prediction model to obtain a predicted Sentinel-2LAI aggregation value, and inverting the predicted Sentinel-2LAI aggregation value to obtain a leaf area index of the heterogeneous earth surface.
2. The heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing according to claim 1, wherein the actual Sentinel-2LAI aggregate value in the third step is obtained based on the following steps:
converting the Sentinel-2LAI data into Sinusidal projection, cutting the Sentinel-2LAI data converted into Sinusidal projection according to the range of a research area, and aggregating the cut Sentinel-2LAI data onto a MODIS pixel grid by using a simple average method to obtain an actual Sentinel-2LAI aggregation value.
3. The heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing according to claim 2, wherein the aggregating the clipped Sentinel-2LAI data onto the MODIS pixel grid by using a simple averaging method specifically comprises the following steps:
the cut Sentinel-2LAI data is corresponding to the interior of each MODIS pixel grid to average,
if each of the Sentinel-2LAI sub-pel values is between 0 and 10, the Sentinel-2LAI sub-pel value is determined to be a valid value, otherwise the Sentinel-2LAI sub-pel value is determined to be an invalid value,
calculating the effective value proportion of the Sentinel-2LAI sub-pixels in each MODIS pixel grid, if the effective value proportion is more than or equal to 90 percent, calculating the average value of the Sentinel-2LAI sub-pixels in the MODIS pixel grid to obtain the actual Sentinel-2LAI aggregate value,
and (3) obtaining the reflectivity of the MODIS pixels corresponding to the MODIS pixel grid through indexing, and if the effective value proportion is smaller than 90%, masking the pixels corresponding to the MODIS pixel grid with the resolution of 500 meters.
4. A heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing according to claim 3, wherein said step four comprises the steps of:
firstly, extracting a set number of sample points to be determined as sample points, and calculating an objective function of the sample points to be determinedOF
Secondly, iteratively replacing a plurality of sample points to be determined, comparing objective functions of the sample points to be determined before and after optimization, accepting the iterative replacement of the sample points to be determined if the objective function after the iterative replacement is smaller than the objective function before the iterative replacement, otherwise refusing to discard the iterative replacement of the sample points to be determined last time;
and finally, outputting each sample point to be determined as a representative sample point if the objective function is smaller than the threshold value or reaches the preset iteration times.
5. The heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing of claim 4, wherein said objective functionOFBased on the following formula:
wherein ,is the average deviation of reflectivity, +.>Is vegetation type average deviation, ++>Is the nearest neighbor index, i1 is the number of the reflectivity band in the MOD09A1 reflectivity data,ais the number of reflectivity bands in MOD09A1 reflectivity data.
6. The heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing of claim 5, wherein said average deviation of reflectivityBased on the following formula: />
Where n is the number of samples to be validated,index of reflectance-> and />Is the minimum and maximum reflectivity of the ith reflectivity interval, +.>Indicating that the reflectivity of a certain type of reflectivity Band is greater than or equal to +.>And is less than->Is determined by the number of sample points to be confirmed,
average deviation of vegetation typesBased on the following formula: />
wherein ,mrefers to the total number of vegetation types for the sample points to be validated,is the number of sample points to be confirmed belonging to the type t vegetation,k t is the duty cycle of sample points belonging to the type t vegetation in the LAI-reflection dataset,
the nearest neighbor indexNNIBased on the following formula:
wherein min is%d ij ) Representing the geographical location distance between the sample point i to be confirmed and the nearest sample point j to be confirmed,nis the number of sample points to be validated,Ais the area of the investigation region.
7. The heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing according to claim 6, wherein the prediction model is one or more of an artificial neural network prediction model, a random forest regression model, a gaussian process regression model and a support vector machine regression model;
when the prediction models are multiple, constructing a test set with MOD09A1 reflectivity data and actual Sentinel-2LAI aggregation values as sample points, inputting the MOD09A1 reflectivity data in the test set into the prediction model which is completed by training, outputting corresponding predicted Sentinel-2LAI aggregation values, and calculating root mean square error RMSE, average deviation Bias and decision coefficient R between the actual Sentinel-2LAI aggregation values and the predicted Sentinel-2LAI aggregation values of each prediction model according to the actual Sentinel-2LAI aggregation values and the corresponding predicted Sentinel-2LAI aggregation values in the test set 2 Selecting RMSE to be less than or equal to 0.5, and R 2 A prediction model with a Bias of 0.8 or more and 0.01 or less is used as the optimal prediction model.
CN202310793690.XA 2023-06-30 2023-06-30 Heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing Active CN116503747B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310793690.XA CN116503747B (en) 2023-06-30 2023-06-30 Heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310793690.XA CN116503747B (en) 2023-06-30 2023-06-30 Heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing

Publications (2)

Publication Number Publication Date
CN116503747A true CN116503747A (en) 2023-07-28
CN116503747B CN116503747B (en) 2023-09-15

Family

ID=87330598

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310793690.XA Active CN116503747B (en) 2023-06-30 2023-06-30 Heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing

Country Status (1)

Country Link
CN (1) CN116503747B (en)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2014122256A1 (en) * 2013-02-08 2014-08-14 Universita' Degli Studi Di Milano Method and electronic equipment for determining a leaf area index
US20150339540A1 (en) * 2012-12-26 2015-11-26 Nec Corporation Image measuring method, system, device, and program
CN108662991A (en) * 2018-04-08 2018-10-16 浙江大学 Plot scale leaves of winter wheat area index evaluation method based on remote sensing satellite data
CN109886962A (en) * 2019-03-27 2019-06-14 南京林业大学 The remote sensing inversion method of vegetative coverage and the management measure factor based on LAI and data of multiple angles
CN110579186A (en) * 2019-08-26 2019-12-17 中国农业大学 Crop growth monitoring method based on inversion of leaf area index by inverse Gaussian process
CN111242022A (en) * 2020-01-10 2020-06-05 西安科技大学 High-resolution FAPAR estimation method based on low-resolution remote sensing product downscaling
CN113553697A (en) * 2021-06-24 2021-10-26 中国矿业大学(北京) Long-time-sequence multi-source data-based vegetation disturbance analysis method for coal mining
CN115526098A (en) * 2022-09-14 2022-12-27 国家能源投资集团有限责任公司 Remote sensing calculation method for leaf area index of surface vegetation in mining area and electronic equipment

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150339540A1 (en) * 2012-12-26 2015-11-26 Nec Corporation Image measuring method, system, device, and program
WO2014122256A1 (en) * 2013-02-08 2014-08-14 Universita' Degli Studi Di Milano Method and electronic equipment for determining a leaf area index
CN108662991A (en) * 2018-04-08 2018-10-16 浙江大学 Plot scale leaves of winter wheat area index evaluation method based on remote sensing satellite data
CN109886962A (en) * 2019-03-27 2019-06-14 南京林业大学 The remote sensing inversion method of vegetative coverage and the management measure factor based on LAI and data of multiple angles
CN110579186A (en) * 2019-08-26 2019-12-17 中国农业大学 Crop growth monitoring method based on inversion of leaf area index by inverse Gaussian process
CN111242022A (en) * 2020-01-10 2020-06-05 西安科技大学 High-resolution FAPAR estimation method based on low-resolution remote sensing product downscaling
CN113553697A (en) * 2021-06-24 2021-10-26 中国矿业大学(北京) Long-time-sequence multi-source data-based vegetation disturbance analysis method for coal mining
CN115526098A (en) * 2022-09-14 2022-12-27 国家能源投资集团有限责任公司 Remote sensing calculation method for leaf area index of surface vegetation in mining area and electronic equipment

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
陈丽等: "基于混合像元分解模型的森林叶面积指数反演", 农业工程学报, no. 13 *

Also Published As

Publication number Publication date
CN116503747B (en) 2023-09-15

Similar Documents

Publication Publication Date Title
Huang et al. Assimilation of remote sensing into crop growth models: Current status and perspectives
CN109884664B (en) Optical microwave collaborative inversion method and system for urban overground biomass
CN108921885B (en) Method for jointly inverting forest aboveground biomass by integrating three types of data sources
Halme et al. Utility of hyperspectral compared to multispectral remote sensing data in estimating forest biomass and structure variables in Finnish boreal forest
CN106372730B (en) Utilize the vegetation net primary productivity remote sensing estimation method of machine learning
CN110751094A (en) Crop yield estimation technology based on GEE comprehensive remote sensing image and deep learning method
Babcock et al. Modeling forest biomass and growth: Coupling long-term inventory and LiDAR data
CN114740180A (en) Soil organic carbon estimation method and device based on multi-source remote sensing data
CN114529097B (en) Multi-scale crop phenological period remote sensing dimensionality reduction prediction method
Aboelghar et al. Spectral wheat yield prediction modeling using SPOT satellite imagery and leaf area index
Yuan et al. Research on rice leaf area index estimation based on fusion of texture and spectral information
CN113205014A (en) Time sequence data farmland extraction method based on image sharpening
CN115128013A (en) Soil organic matter content space prediction evaluation method based on partition algorithm
CN117036928A (en) Urban forest carbon sink measuring and calculating method and system based on remote sensing images and machine learning
CN117520783A (en) Crop yield prediction method
CN116503747B (en) Heterogeneous earth surface leaf area index inversion method based on multi-scale remote sensing
US20230108422A1 (en) Methods and systems for use in processing images related to crops
CN114782835B (en) Crop lodging area proportion detection method and device
Li et al. A technique system for the measurement, reconstruction and character extraction of rice plant architecture
Kanna et al. A maize crop yield optimization and healthcare monitoring framework using firefly algorithm through iot
XZ et al. A new soybean NDVI data-based partitioning algorithm for fertilization management zoning.
CN114611699A (en) Soil moisture downscaling method and device, electronic equipment and storage medium
Jihua et al. Design, development and application of a satellite-based field monitoring system to support precision farming
Bhadra et al. PROSAIL-Net: A transfer learning-based dual stream neural network to estimate leaf chlorophyll and leaf angle of crops from UAV hyperspectral images
Jia et al. Quantitative analysis and hyperspectral remote sensing inversion of rice canopy spad in A cold region

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant