Large-range rice planting area extraction method based on GNSS-R remote sensing data

Abstract

The invention discloses a large-range rice planting area extraction method based on GNSS-R remote sensing data, and the method comprises the steps: obtaining the remote sensing data of a target region, and carrying out the preprocessing of the remote sensing data; calculating a first surface reflectance of each reflection point in the target area by using the preprocessed cyclone satellite navigation system receiving data; establishing a relationship between the reflectivity of each grid and the reflectivity of the adjacent grid, and performing interpolation operation on the first surface reflectivity; performing mountain and non-mountain division on the target area according to the slope data, and performing paddy field identification on the divided target area according to the second surface reflectance and the vegetation index; and calculating the rice planting area according to an identification result. According to the method, the water logging characteristics of the paddy field in the special growth period are utilized, the reflectivity calculated by CYGNSS original data is used as an index, and a simpler and more efficient threshold value method is adopted, so that rapid large-scale paddy rice identification is realized.

Description

Technical Field

[0001] The invention relates to the field of agricultural remote sensing technology, and in particular to a method for extracting a large-scale rice planting area based on GNSS-R remote sensing data. Background Art

[0002] Against the backdrop of global climate change and population growth, it is increasingly important to obtain real-time information on large-scale rice field areas, which can provide key data support for formulating scientific and reasonable agricultural policies, adjusting food production layout, and optimizing agricultural structure.

[0003] In terms of rice field area identification, the more commonly used method is to use the backscatter coefficient combined with phenological characteristics for identification, but this method relies on the changes in radar backscatter characteristics during rice growth, which may dilute the unique flooding characteristics of rice fields during the growth cycle. Because some crops have similar backscatter dynamic ranges, it may not be possible to distinguish rice from other crops.

[0004] In recent years, remote sensing technology has made significant breakthroughs in the field of rice identification. Currently, the two main remote sensing data sources widely used for rice identification are optical remote sensing data and synthetic aperture radar (SAR) data, but they have the following problems in application: Optical remote sensing reflects rice growth by capturing vegetation growth status information, such as leaf area index or gross primary productivity, but its application is limited by weather conditions; Although SAR data can penetrate and observe under cloud cover and has a higher tolerance for weather changes, its huge amount of data and complex preprocessing process also place higher demands on computing resources. It has disadvantages in large-scale monitoring and can often only be applied at the county level, making it difficult to meet the problem of rapid identification over a large area.

[0005] In summary, conventional rice identification and planting area extraction have problems in both data source and identification method, and cannot achieve timely update of large-scale rice planting area data and rapid acquisition of information. Summary of the invention

[0006] Purpose of the invention: In view of the above problems, the purpose of the present invention is to provide a method for extracting a large-scale rice planting area based on GNSS-R remote sensing data.

[0007] Technical solution: The method for extracting a large-scale rice planting area based on GNSS-R remote sensing data of the present invention comprises the following steps: Step 1, obtaining remote sensing data of the target area, including the received data of the Cyclone satellite navigation system, the relevant data of the medium resolution imaging spectrometer and the digital elevation model data, and preprocessing the remote sensing data; Step 2, using the pre-processed received data of the Cyclone satellite navigation system to calculate the first surface reflectivity of each reflection point in the target area; Step 3, dividing the target area into grids, establishing a relationship between the reflectivity of each grid and its adjacent grids, and performing an interpolation operation on the first surface reflectivity to obtain a second surface reflectivity; Step 4, dividing the target area into mountainous areas and non-mountainous areas according to the slope data, and identifying rice fields in the divided target area according to the second surface reflectance and vegetation index; Step 5: Calculate the rice planting area based on the recognition results.

[0008] Furthermore, the remote sensing data is preprocessed including: Perform quality control on the data received by the Cyclone satellite navigation system, remove observation data with an incident angle greater than E0, and remove data with a signal-to-noise ratio less than I1 or greater than I2, a receiver gain less than 0 dB, and a reflected signal peak power greater than P0; where I1 is less than I2; The relevant data of the medium-resolution imaging spectrometer are converted in format, projected and resampled, and then the vegetation index is calculated; the calculation formula of the vegetation index is: , In the formula, ρ R is the reflectivity in the red light band, ρ NIR is the reflectivity in the near-infrared band; Remove the noise and outliers in the digital elevation model data, and then calculate the slope data; the slope calculation formula is: , In the formula, z represents the height of the ground surface, w and v Represents geographic horizontal and vertical coordinates.

[0009] Furthermore, the calculation formula of the first surface reflectivity is: , Where λ represents the wavelength of the GPS signal, P rl coh represents the coherent scattering component power, r st represents the distance between the specular reflection point and the GNSS transmitter, r sr Indicates the distance between the mirror reflection point and the receiver, P t Indicates the transmission power of the signal, Gt represents the gain of the transmitting antenna, G r Represents the gain of the receiving antenna.

[0010] Further, step 3 includes: The target area is evenly divided into P×Q grids, and the eight grid cells surrounding each grid point are defined as adjacent grid points; where P and Q are positive real numbers; Remove cells in the grid that contain water surfaces; The reflectivity of each grid point in the target area x With its i The reflectivity of the neighboring grid points y Between y= a i x+b i Linear regression fitting, get the intercept b i With slope a i , any of which x The sampling time interval between the corresponding y does not exceed 24 hours; Calculate the reflectivity correlation coefficient between the target grid point and the adjacent grid points. The calculation formula is: , In the formula, y i is the target grid point i Observed reflectivity; x i is the observed reflectivity at adjacent grid points with a sampling interval not exceeding 24 hours; and Both represent the mean values ​​of all reflectivities of the corresponding grid points in the historical period; The first surface reflectivity is weighted interpolated using the surface reflectivity during the specific rice growth period and winter in the target year to obtain the second surface reflectivity. The calculation formula is: , In the formula, G i For the i The observed reflectivity of adjacent grid points is n represents the total number of nearby grid points with reflectivity data in the target year and target period, w i is the target grid point and i The weighted weight between adjacent grid points is calculated as follows: .

[0011] Further, step 4 includes: Compare target area i grid cells correspond to slope S i With slope threshold S 0 The size of S i > S 0 , then the grid cell is classified as mountainous area, otherwise the grid cell is classified as non-mountainous area; If the grid cell is classified as mountainous, compare the vegetation index NDVI with the vegetation index threshold NDVI0, and the reflectivity difference Δ with the reflectivity difference threshold Δ0. If NDVI>NDVI0, and Δ > Δ0, the grid cell belongs to a rice field, otherwise it does not; = Δ1 - Δ2, Δ1 represents the second surface reflectance of the rice field in a specific growth period, and Δ2 represents the second surface reflectance of the rice field in winter; If the grid cell is classified as non-mountainous area, the vegetation index NDVI is compared with the vegetation index threshold NDVI0, and the second surface reflectance SR1 of the rice field in a specific growth period is compared with the reflectance threshold SR0. If NDVI>NDVI0, and SR1 > If SR0, the grid cell belongs to a rice field, otherwise it does not.

[0012] Further, step 5 includes: The number of grid cells corresponding to the screened rice fields was multiplied by the scale of each grid cell to obtain the rice planting area.

[0013] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages: The present invention utilizes the flooding characteristics of rice fields in special growth periods, takes the reflectivity calculated using CYGNSS original data as an indicator that can well characterize the amount of soil moisture, adopts a simpler and more efficient threshold method, and takes advantage of the relatively small amount of data of CYGNSS to achieve rapid large-scale rice recognition and further extract the rice planting area; simplifies the complex machine learning classification method into an ordinary threshold method, which is simpler and faster; applies CYGNSS technology to rice recognition for the first time, expands the application of CYGNSS in agricultural remote sensing, and provides new technologies and means for large-scale rapid identification of rice fields. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1It is a flowchart of the large-scale rice planting area extraction method based on GNSS-R remote sensing data; Figure 2 It is a flowchart of the large-scale rice planting area extraction method based on GNSS-R remote sensing data; Figure 3 This is the sampling result diagram of CYGNSS data; Figure 4 This is the distribution map of NDVI results; Figure 5 is the slope result diagram; Figure 6 This is the interpolation result map of winter surface reflectivity; Figure 7 This is the interpolation result map of the surface reflectance in a specific growth period; Figure 8 This is the rice planting area identification result map. DETAILED DESCRIPTION

[0015] In order to make the objectives, technical solutions and advantages of the present application more clearly understood, the present application is further described in detail below in conjunction with the accompanying drawings and embodiments.

[0016] The method for extracting large-scale rice planting area based on GNSS-R remote sensing data described in this embodiment is as follows: Figure 1 and Figure 2 As shown, it includes the following steps 1 to 5.

[0017] Step 1: obtain remote sensing data of the target area, including the received data of the Cyclone satellite navigation system, the relevant data of the medium-resolution imaging spectrometer and the digital elevation model data, and pre-process the remote sensing data.

[0018] In this example, the Cyclone Satellite Navigation System (CYGNSS) data is used to calculate the surface reflectivity, such as Figure 3 The figure shows the CYGNSS sampling results obtained on July 1, 2022 using a dot sampling method. The Moderate Resolution Imaging Spectroradiometer (MODIS) data are used to exclude the influence of other water bodies, and the Digital Elevation Model (DEM) data are used to calculate the slope and further classify and identify mountainous and non-mountainous areas. The CYGNSS satellite data comes from the Level 1 V3.0 data provided by NASA, the MODIS data also comes from NASA, and the DEM data comes from the Shuttle Radar Topography Mission SRTM of NASA and the National Geospatial-Intelligence Agency.

[0019] Step 2, using the pre-processed received data of the Cyclone satellite navigation system to calculate the first surface reflectivity of each reflection point in the target area; Step 3, dividing the target area into grids, establishing a relationship between the reflectivity of each grid and its adjacent grids, and performing an interpolation operation on the first surface reflectivity to obtain a second surface reflectivity; Step 4, dividing the target area into mountainous areas and non-mountainous areas according to the slope data, and identifying rice fields in the divided target area according to the second surface reflectance and vegetation index; Step 5: Calculate the rice planting area based on the recognition results.

[0020] Furthermore, the remote sensing data is preprocessed including: Perform quality control on the data received by the Cyclone satellite navigation system, remove observation data with an incident angle greater than E0, and remove data with a signal-to-noise ratio less than I1 or greater than I2, a receiver gain less than 0 dB, and a reflected signal peak power greater than P0; where I1 is less than I2; The relevant data of the medium-resolution imaging spectrometer are format converted, projected, and resampled, and then the vegetation index is calculated. The relevant data include the reflectance of the red light band of 620nm~670nm and the reflectance of the near-infrared band of 841~875nm. Figure 4 This is the distribution map of vegetation index results, where the calculation formula of vegetation index is: , In the formula, ρ R is the reflectivity in the red light band, ρ NIR It is the reflectivity of the near-infrared band; a positive NDVI value indicates the surface vegetation coverage, that is, the greater the coverage, the greater the NDVI value.

[0021] Remove noise and outliers from the digital elevation model data and then calculate the slope data. Figure 5 is the slope result diagram; the slope calculation formula is: , In the formula, z represents the height of the ground surface, w and v Represents geographic horizontal and vertical coordinates.

[0022] In this example, the acquired CYGNSS satellite data is quality controlled. The data group can be quality controlled according to the variable quality flag quality_flags provided in its metadata. Then the observation data with an incident angle greater than 65° is removed to preliminarily ensure that the error caused by the terrain height is within a reasonable range. Then the data with a signal-to-noise ratio less than 2dB or greater than 14dB, a receiver gain Gr less than 0dB, and a peak power of the reflected signal greater than -147dB are removed to complete the preprocessing of the CYGNSS satellite data.

[0023] When calculating the slope, the DEM data can be converted into slope data using geographic information system software to obtain 90m × 90m slope data, which can then be interpolated to a consistent 250m × 250m scale.

[0024] CYGNSS satellite data has a unique sampling method. CYGNSS only collects data information of reflection points. Generally, moist soil has a higher dielectric constant, so the reflectivity will be higher. The surface humidity can be estimated by analyzing the strength of the GPS signal reflected from the surface. That is, the stronger the reflected signal, the higher the surface humidity, and it can even be regarded as a water surface. In the case of land reflection, it can be assumed that the surface is relatively smooth and approximated as a mirror reflection, that is, the surface reflectivity can be calculated based on the power of the coherent scattering component. The reflectivity of each sampling and each reflection point can be calculated by the formula.

[0025] Furthermore, the calculation formula of the first surface reflectivity is: , Where λ represents the wavelength of the GPS signal, P rl coh represents the coherent scattering component power, r st represents the distance between the specular reflection point and the GNSS transmitter, r sr Indicates the distance between the mirror reflection point and the receiver, P t Indicates the transmission power of the signal, G t represents the gain of the transmitting antenna, G r Represents the gain of the receiving antenna.

[0026] Further, step 3 includes: The target area is evenly divided into P×Q grids, and the eight grid cells surrounding each grid point are defined as adjacent grid points; where P and Q are positive real numbers; Remove cells in the grid that contain water surfaces; Select a certain historical period and calculate the reflectivity of each grid point in the target area. x With its i The reflectivity of the neighboring grid points y Between y=a i x+b i Linear regression fitting, get the intercept b i With slopea i , any of which x The sampling time interval between the corresponding y does not exceed 24 hours; Calculate the reflectivity correlation coefficient between the target grid point and the adjacent grid points. The calculation formula is: , In the formula, y i is the target grid point i Observed reflectivity; x i is the observed reflectivity at adjacent grid points with a sampling interval not exceeding 24 hours; and Both represent the mean of all reflectivities of the corresponding grid points in the historical period. In this example, the historical period refers to the period from 2019 to 2021.

[0027] The first surface reflectance of rice in the Jianghuai region of Anhui Province in January, July and August 2022 is weighted interpolated. Since the specific growth period of rice is often greater than 30 days, there are usually multiple adjacent grid points at the target grid point. At this time, a weighted method is used for interpolation to obtain the second surface reflectance. The ground reflectance results for January, July and August 2022 that need to be calculated and interpolated in this example are shown in Figure 6 and Figure 7 The second surface reflectivity calculation formula is: , In the formula, G i For the i The observed reflectivity of adjacent grid points is n represents the total number of nearby grid points with reflectivity data in the target year and target period, w i is the target grid point and i The weighted weight between adjacent grid points is calculated as follows: .

[0028] In the example, CYGNSS satellite data is sampled in the form of reflection points. The linear trajectory of the observation point is the same as the CYGNSS orbit, and the specific observation position is determined by the position of the receiver and the GNSS transmitter. While this sampling method brings many advantages, it also determines that not every grid point on the 250m×250m grid has data, so it is necessary to use the grid points of the existing data to reasonably infer the data on each grid point. In this example, the POBI interpolation method is used to divide the target area data into a 250m×250m grid, and define the 8 grid cells surrounding each grid point as adjacent grid points. The water surface grid is removed, and the global surface water dataset of the European Joint Research Center is used to remove the main water bodies in the target area to reduce the impact of the reflection points on the water surface on the surrounding soil moisture. Then, the surface reflectance data from the historical period of 2019-2021 is used to establish a relationship model between the target grid point and the adjacent grid point, which is used to interpolate the target year data to obtain 250m×250m surface reflectance data consistent with the spatial scale of MODIS data. Specifically, the historical data of the past year or longer sequence are used to perform linear fitting on the reflectance data of each grid point in the target area and its corresponding adjacent grid point with a sampling time interval of no more than 24 hours to obtain the relationship model. y=a i x+b i , and then get the intercept b With slope a Then calculate the reflectivity correlation coefficient between the target grid point and the adjacent grid points, and use the intercept b , slope a And the reflectivity correlation coefficient r The surface reflectance conditions during the specific rice growing period and winter in the target area and target year are weighted interpolated to obtain the surface reflectance distribution results at a scale of 250m×250m in the two periods of the specific rice growing period and winter in the target year.

[0029] Further, step 4 includes: Compare target area i grid cells correspond to slope S i With slope threshold S 0 The size of S i > S 0 , then the grid cell is classified as mountainous area, otherwise the grid cell is classified as non-mountainous area; If the grid cell is classified as mountainous, compare the vegetation index NDVI with the vegetation index threshold NDVI0, and the reflectivity difference Δ with the reflectivity difference threshold Δ0. > Δ0, the grid cell belongs to a rice field, otherwise it does not; = SR1 - SR2, SR1 represents the second surface reflectance of the rice field in a specific growth period, and SR2 represents the second surface reflectance of the rice field in winter; If the grid cell is classified as non-mountainous area, the vegetation index NDVI is compared with the vegetation index threshold NDVI0, and the second surface reflectance SR1 of the rice field in a specific growth period is compared with the reflectance threshold SR0. If NDVI>NDVI0, and SR1 >SR0, the grid cell belongs to a rice field, otherwise it does not.

[0030] In the example, the Jianghuai region of Anhui Province is taken as an example, including Hefei, Wuhu, Huainan, Ma'anshan, Tongling, Anqing, Chuzhou, and Lu'an. The specific growth period can be selected as July and August, and the winter can be selected as January. The slope threshold is 1°, the vegetation index threshold NDVI0 is 0.43, the second surface reflectance threshold is 159.5, and the reflectance difference is 3.9. That is, when the slope of a grid is greater than 1°, it is judged as a mountainous area. If its NDVI is greater than 0.43 and the difference between the reflectance in July and August and the reflectance in January is greater than 3.9, then the grid belongs to a rice field, otherwise it is not; when the slope of a grid is less than 1°, it is judged as a non-mountainous area. If its NDVI is greater than 0.43 and the surface reflectance in July and August is greater than 159.5, then the grid belongs to a rice field, otherwise it is not.

[0031] Since higher surface reflectivity often represents higher soil moisture, the flooding characteristics in the rice growth cycle, that is, the rice fields are soaked for a long time during a specific growth period, are used to set a threshold for the interpolated surface reflectivity through preliminary survey samples, and grid points with values ​​less than the threshold are excluded. Due to the influence of terrain, vegetation coverage, radiation, etc., the seasonal differences in soil moisture in mountainous areas are often large, so the difference in soil moisture between a specific growth period and the dry winter period is used as an indicator for identification. Taking into account the influence of other water bodies, the NDVI data obtained after processing is used to regard grid points with NDVI values ​​less than the threshold as ordinary water bodies, ultimately achieving rice identification.

[0032] Further, step 5 includes: The rice planting area is obtained by multiplying the number of grid cells corresponding to the selected rice fields by the scale of each grid cell. In this example, the scale is 250 meters × 250 meters.

[0033] In this example, the rice field data in the third land survey report of each city in Anhui Province is used as the true value for verification, and the relative error is used as the measurement indicator. After preliminary experiments and field investigations, NDVI0=0.43, SR=159.5, Δ=3.9 are selected as thresholds for rice identification and area extraction. The final average relative error is 4.45%. The relative errors of each city are shown in Table 1. The rice identification results are shown in Table 1. Figure 8 shown.

[0034] Table 1

[0035] From the results, we can see that the accuracy of the results can be basically controlled at about 5%, and the total accuracy of the eight cities in Anhui, Jianghuai and Anhui can reach less than 1%. Among them, the poor accuracy of Ma'anshan may be due to its relatively small planting area and higher sensitivity. In short, this technical solution of the present invention can extract the rice planting area very well.

Claims

1. A large-scale rice planting area extraction method based on GNSS-R remote sensing data is characterized by: The steps include: Step 1, obtaining remote sensing data of the target area, including the received data of the Cyclone satellite navigation system, the relevant data of the medium resolution imaging spectrometer and the digital elevation model data, and preprocessing the remote sensing data; Step 2, using the pre-processed received data of the Cyclone satellite navigation system to calculate the first surface reflectivity of each reflection point in the target area; Step 3, dividing the target area into grids, establishing a relationship between the reflectivity of each grid and its adjacent grids, performing an interpolation operation on the first surface reflectivity to obtain a second surface reflectivity; Step 4, dividing the target area into mountainous areas and non-mountainous areas according to the slope data, and identifying rice fields in the divided target area according to the second surface reflectance and vegetation index; Step 5: Calculate the rice planting area based on the recognition results.

2. The method for extracting large-scale rice planting area based on GNSS-R remote sensing data according to claim 1, characterized in that: Preprocessing of remote sensing data includes: Perform quality control on the data received by the Cyclone satellite navigation system, remove observation data with an incident angle greater than E0, and remove data with a signal-to-noise ratio less than I1 or greater than I2, a receiver gain less than 0 dB, and a reflected signal peak power greater than P0; where I1 is less than I2; The relevant data of the medium-resolution imaging spectrometer are converted in format, projected and resampled, and then the vegetation index is calculated; the calculation formula of the vegetation index is: ， In the formula, ρ R is the reflectivity in the red light band, ρ NIR is the reflectivity in the near-infrared band; Remove the noise and outliers in the digital elevation model data, and then calculate the slope data; the slope calculation formula is: ， In the formula, z represents the height of the ground surface, w and v Represents geographic horizontal and vertical coordinates.

3. The method for extracting large-scale rice planting area based on GNSS-R remote sensing data according to claim 2, characterized in that: The calculation formula of the first surface reflectivity is: ， Where λ represents the wavelength of the GPS signal, P rl coh represents the coherent scattering component power, r st represents the distance between the specular reflection point and the GNSS transmitter, r sr Indicates the distance between the mirror reflection point and the receiver, P t Indicates the transmission power of the signal, G t represents the gain of the transmitting antenna, G r Represents the gain of the receiving antenna.

4. The method for extracting large-scale rice planting area based on GNSS-R remote sensing data according to claim 3, characterized in that: Step 3 includes: The target area is evenly divided into P×Q grids, and the eight grid cells surrounding each grid point are defined as adjacent grid points; where P and Q are positive real numbers; Remove cells in the grid that contain water surfaces; The reflectivity of each grid point in the target area x With its i The reflectivity of the neighboring grid points y Between y=a i x+ b i Linear regression fitting, get the intercept b i With slope a i , any of which x The corresponding y The sampling interval shall not exceed 24 hours; Calculate the reflectivity correlation coefficient between the target grid point and the adjacent grid points. The calculation formula is: ， In the formula, y i is the target grid point i The observed reflectivity, x i is the observed reflectivity at adjacent grid points with a sampling interval not exceeding 24 hours, and represents the mean value of all reflectances at the corresponding grid points; The first surface reflectivity is weighted interpolated using the surface reflectivity during the specific rice growth period and winter in the target year to obtain the second surface reflectivity. The calculation formula is: ， In the formula, G i For the i The observed reflectivity of adjacent grid points is n It represents the total number of nearby grid points with reflectivity data in the target year and target period; w i is the target grid point and i The weighted weight between adjacent grid points is calculated as follows: 。 5. The method for extracting large-scale rice planting area based on GNSS-R remote sensing data according to claim 4, characterized in that: Step 4 includes: Compare target area i grid cells correspond to slope S i With slope threshold S 0 The size of S i > S 0 , then the grid cell is classified as mountainous area, otherwise the grid cell is classified as non-mountainous area; If the grid cell is classified as a mountainous area, the vegetation index NDVI is compared with the vegetation index threshold NDVI0, and the reflectance difference Δ is compared with the reflectance difference threshold Δ0. If NDVI > NDVI0, and Δ> Δ0, the grid cell belongs to a rice field, otherwise it does not. = Δ1 - Δ2, Δ1 represents the second surface reflectance of the rice field in a specific growth period, and Δ2 represents the second surface reflectance of the rice field in winter; If the grid cell is classified as non-mountainous area, the vegetation index NDVI is compared with the vegetation index threshold NDVI0, as well as the second surface reflectance SR1 of the rice field in a specific growth period and the reflectance threshold SR0. If NDVI > NDVI0, and SR1>SR0, the grid cell belongs to a rice field, otherwise it does not.

6. The method for extracting large-scale rice planting area based on GNSS-R remote sensing data according to claim 5, characterized in that: Step 5 includes: The number of grid cells obtained from screening rice fields is multiplied by the scale of each grid cell to obtain the rice planting area.

Images