Method for extracting large-scale rice planting area based on GNSS-R remote sensing data

Abstract

The present invention discloses a method for extracting the large-scale rice planting area based on GNSS-R remote sensing data, including: obtaining remote sensing data of a target area and preprocessing the remote sensing data; calculating the first surface reflectivity of each reflection point in the target area by using the data received by the Cyclone satellite navigation system after preprocessing; establishing the relationship between the reflectivity of each grid and that of adjacent grids, and performing interpolation operation on the first surface reflectivity; dividing the target area into mountainous and non-mountainous areas according to the slope data, and identifying paddy fields in the divided target area according to the second surface reflectivity and vegetation index; calculating the rice planting area according to the identification result. The present invention utilizes the waterlogging characteristics of paddy fields in a special growth period, uses the reflectivity calculated from the CYGNSS raw data as an index, and adopts a simpler and more efficient threshold method to achieve rapid large-scale rice identification.

Description

Technical Field

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

[0002] In the context of global climate change and population growth, it has become increasingly important to master the real-time information of the large-scale rice paddy area, which can provide key data support for formulating scientific and reasonable agricultural policies, adjusting the grain production layout, and optimizing the agricultural structure.

[0003] In the way of identifying the rice paddy area, it is more commonly used to identify by using the backscattering coefficient in combination with phenological characteristics. However, this method depends on the change of the radar backscattering characteristics during the rice growth period, and these characteristics may dilute the unique waterlogging characteristics of the rice paddy during the growth cycle. Since some crops have a similar backscattering dynamic range, it may not be possible to distinguish rice from other crops.

[0004] In recent years, remote sensing technology has made remarkable breakthroughs in the field of rice identification. At present, 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:

[0005] Optical remote sensing reflects the growth of rice by capturing the growth state information of vegetation, such as leaf area index or gross primary productivity, but its application is limited by weather conditions;

[0006] Although SAR data can penetrate and observe under cloud cover and has a higher tolerance to weather changes, its huge data volume and complex preprocessing process also pose higher requirements for computing resources, and there are disadvantages in large-scale monitoring. It can often only be applied to the county level and is difficult to meet the problem of large-scale and rapid identification.

[0007] In summary, there are problems in both data sources and identification methods for conventional rice identification and planting area extraction, and it is impossible to realize the timely update of large-scale rice planting area data and the rapid acquisition of information. Summary of the Invention

[0008] Object of the Invention: Aiming at the above problems, the object of the present invention is to provide a method for extracting the large-scale rice planting area based on GNSS-R remote sensing data.

[0009] Technical Solution: The method for extracting the large-scale rice planting area based on GNSS-R remote sensing data of the present invention includes the following steps:

[0010] Step 1: Obtain remote sensing data of the target area, including data received by the cyclone satellite navigation system, relevant data of the Moderate Resolution Imaging Spectroradiometer (MODIS), and digital elevation model data, and preprocess the remote sensing data;

[0011] Step 2: Calculate the first surface reflectance of each reflection point in the target area using the data received by the cyclone satellite navigation system after preprocessing;

[0012] Step 3: Divide the target area into grids, establish the relationship between the reflectance of each grid and that of adjacent grids, and perform interpolation on the first surface reflectance to obtain the second surface reflectance;

[0013] Step 4: Divide the target area into mountainous and non-mountainous areas based on the slope data, and identify paddy fields in the divided target area according to the second surface reflectance and vegetation index;

[0014] Step 5: Calculate the paddy planting area according to the identification result.

[0015] Furthermore, the preprocessing of the remote sensing data includes:

[0016] Perform quality control on the data received by the cyclone satellite navigation system, eliminate the observation data with the incident angle greater than E 0 and eliminate the data with the signal-to-noise ratio less than I 1 or greater than I 2 , the receiver gain less than 0 dB, and the peak power of the reflected signal greater than P 0 ; where I 1 is less than I 2 ;

[0017] Perform format conversion, projection conversion, and data resampling on the relevant data of the Moderate Resolution Imaging Spectroradiometer, and then calculate the vegetation index; the calculation formula of the vegetation index is:

[0018] ,

[0019] In the formula, ρ R is the reflectance of the red light band, ρ NIR is the reflectance of the near-infrared band;

[0020] Remove the noise and outliers existing in the digital elevation model data, and then calculate the slope data; the calculation formula of the slope is:

[0021] ,

[0022] In the formula, z represents the surface height, w and vIndicates the geographical horizontal and vertical coordinates.

[0023] Furthermore, the calculation formula for the first surface reflectivity is:

[0024] ,

[0025] In the formula, λ represents the GPS signal wavelength, P rl coh represents the power of the coherent scattering component, r st represents the distance between the specular reflection point and the GNSS transmitter, r sr represents the distance between the specular reflection point and the receiver, P t represents the transmit power of the signal, G t represents the gain of the transmit antenna, G r represents the gain of the receive antenna.

[0026] Furthermore, step 3 includes:

[0027] The target area is evenly divided into a P×Q grid, and the 8 grid cells surrounding each grid point are defined as adjacent grid points; where P and Q are positive real numbers;

[0028] Remove the cells containing the water surface in the grid;

[0029] For the reflectivity of each grid point in the target area x and its i th adjacent grid point's reflectivity y perform y= a i x + b i linear regression fitting to obtain the intercept b i and the slope a i , where the sampling time interval between any x and its corresponding y does not exceed 24 hours;

[0030] Calculate the correlation coefficient of reflectivity between the target grid point and the adjacent grid point. The calculation formula is:

[0031] ,

[0032] In the formula, y i is the i th observed reflectivity of the target grid point;x i is the observed reflectance at adjacent grid points with a sampling time interval not exceeding 24 hours; and both represent the mean value of all reflectances at the corresponding grid points over the historical period;

[0033] The first surface reflectance is weighted and interpolated using the surface reflectance conditions during the specific growth stage of rice and in winter in the target year to obtain the second surface reflectance. The calculation formula is:

[0034] ,

[0035] In the formula, Г i is the observed reflectance of the i th adjacent grid point, n represents the total number of adjacent grid points with reflectance data in the target year and target period, w i is the weighted weight between the target grid point and the i th adjacent grid point. The calculation formula is:

[0036] .

[0037] Furthermore, step 4 includes:

[0038] Compare the slope i corresponding to the S i of the S 0 th grid cell in the target area with the slope threshold S i . S 0 If

[0039] > 0 0 0 0 1 2 1 2 = - 0 0 1 2 1 2 1 represents the second surface reflectance of paddy fields during the specific growth stage, and 2Represents the second surface reflectance of winter paddy fields;

[0040] If the grid cell is divided as non-mountainous, then compare the vegetation index NDVI with the vegetation index threshold NDVI 0 and the second surface reflectance SR of paddy fields at a specific growth stage 1 with the reflectance threshold SR 0 If NDVI > NDVI 0 and SR 1 (>) SR 0 , then the grid cell belongs to paddy fields, otherwise it does not.

[0041] Further, step 5 includes:

[0042] Multiply the number of grid cells corresponding to the screened paddy fields by the scale of each grid cell to obtain the planting area of rice.

[0043] Beneficial effects: Compared with the prior art, the significant advantages of the present invention are:

[0044] The present invention utilizes the waterlogging characteristics of paddy fields at special growth stages, uses the reflectance calculated from CYGNSS raw data, which can well represent the soil moisture content, as an index, adopts a simpler and more efficient threshold method, gives play to the advantage of relatively small data volume of CYGNSS, realizes rapid large-scale rice identification, and further extracts the rice planting area; simplifies the complex machine learning classification method into a common threshold method, which is more simple and fast; applies the CYGNSS technology to rice identification for the first time, expands the application of CYGNSS in agricultural remote sensing, and provides new technologies and new means for rapid large-scale identification of paddy fields. Description of the Drawings

[0045] Figure 1 Is a flowchart of a method for extracting large-scale rice planting areas based on GNSS-R remote sensing data;

[0046] Figure 2 Is a flow block diagram of a method for extracting large-scale rice planting areas based on GNSS-R remote sensing data;

[0047] Figure 3 Is a sampling result diagram of CYGNSS data;

[0048] Figure 4 Is a distribution diagram of NDVI results;

[0049] Figure 5 Is a slope result diagram;

[0050] Figure 6It is the interpolation result map of the surface reflectivity in winter;

[0051] Figure 7 It is the interpolation result map of the surface reflectivity at a specific growth stage;

[0052] Figure 8 It is the recognition result map of the rice planting area. Specific implementation manners

[0053] In order to make the objectives, technical solutions and advantages of the present application clearer, the present application will be further described in detail below with reference to the accompanying drawings and embodiments.

[0054] The method for extracting the large-scale rice planting area based on GNSS-R remote sensing data described in this embodiment has a flowchart as Figure 1 and Figure 2 shown, and includes the following steps 1 to 5.

[0055] Step 1: Obtain the remote sensing data of the target area, including the data received by the Cyclone Global Navigation Satellite System (CYGNSS), the relevant data of the Moderate Resolution Imaging Spectroradiometer (MODIS), and the Digital Elevation Model (DEM) data, and preprocess the remote sensing data.

[0056] In the example, the data received by the Cyclone Global Navigation Satellite System (CYGNSS) is used for the calculation of the surface reflectivity. As Figure 3 shown, it is the CYGNSS sampling result on July 1, 2022 obtained by using the dot sampling method. Obtain the relevant data of the Moderate Resolution Imaging Spectroradiometer (MODIS) to exclude the influence of other water bodies, and obtain the Digital Elevation Model (DEM) data for calculating the slope and further for the classification and recognition of mountains and non-mountains. The CYGNSS satellite data comes from the Level 1 V3.0 data provided by the National Aeronautics and Space Administration (NASA) of the United States. The relevant data of MODIS also comes from NASA, and the DEM data comes from the Shuttle Radar Topography Mission (SRTM) of NASA and the National Geospatial-Intelligence Agency.

[0057] Step 2: Calculate the first surface reflectivity of each reflection point in the target area by using the data received by the preprocessed Cyclone Global Navigation Satellite System.

[0058] Step 3: Divide the target area into grids, establish the relationship between the reflectivity of each grid and its adjacent grids, and perform interpolation on the first surface reflectivity to obtain the second surface reflectivity.

[0059] Step 4: Divide the target area into mountains and non-mountains according to the slope data, and identify the paddy fields in the divided target area according to the second surface reflectivity and the vegetation index.

[0060] Step 5: Calculate the rice planting area according to the recognition result.

[0061] Furthermore, the remote sensing data is preprocessed including:

[0062] The quality control of the received data of the cyclone satellite navigation system is carried out to eliminate the data with an incident angle greater than E 0 Observation data with a signal-to-noise ratio less than I 1 or greater than I 2 , the receiver gain is less than 0dB, and the peak power of the reflected signal is greater than P 0 data; where I 1 Less than I 2 ;

[0063] 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:

[0064] ,

[0065] 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.

[0066] 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:

[0067] ,

[0068] In the formula, z represents the height of the ground surface, w and v Represents geographic horizontal and vertical coordinates.

[0069] 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.

[0070] When calculating the slope, GIS software can be used to convert DEM data into slope data, obtaining slope data of 90 m × 90 m and interpolating it to a consistent scale of 250 m × 250 m.

[0071] CYGNSS satellite data has a unique sampling method. CYGNSS only collects data information of the reflection points. Generally, wet soil has a higher dielectric constant, so the reflectivity is also higher. The surface humidity can be estimated by analyzing the intensity 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 the water surface. In the case of land reflection, it can be assumed that the surface is relatively smooth and approximated as specular reflection, that is, the surface reflectivity can be calculated according to the power of the coherent scattering component. The reflectivity of each sampling and each reflection point can be calculated by the formula.

[0072] Further, the calculation formula of the first surface reflectivity is:

[0073] ,

[0074] In the formula, λ represents the wavelength of the global positioning system signal, P rl coh represents the power of the coherent scattering component, r st represents the distance between the specular reflection point and the GNSS transmitter, r sr represents the distance between the specular reflection point and the receiver, P t represents the transmitting power of the signal, G t represents the gain of the transmitting antenna, G r represents the gain of the receiving antenna.

[0075] Further, step 3 includes:

[0076] The target area is evenly divided into a grid of P × Q, and the 8 grid cells surrounding each grid point are defined as adjacent grid points; where P and Q are positive real numbers;

[0077] Remove the cells containing the water surface in the grid;

[0078] Select a certain historical period, and for each grid point in the target area, the reflectivity x and its i th adjacent grid point's reflectivity y between y = a i x + b iLinear regression fitting to obtain the intercept b i and the slope a i wherein for any one x the sampling time interval between it and the corresponding y does not exceed 24 hours;

[0079] Calculate the correlation coefficient of reflectivity between the target grid point and adjacent grid points. The calculation formula is:

[0080] ,

[0081] In the formula, y i is the i th observed reflectivity of the target grid point; x i is the observed reflectivity at adjacent grid points with a sampling time interval not exceeding 24h; and both represent the mean value of all reflectivities at the corresponding grid points during the historical period. In this example, the historical period refers to the years between 2019 and 2021.

[0082] Perform weighted interpolation on the first surface reflectivity of rice in January, July, and August 2022 in the Jianghuai region of Anhui. Since the specific growth period of rice often exceeds 30 days and there are usually multiple adjacent grid points for the target grid point, weighted interpolation is used at this time to obtain the second surface reflectivity. The results of the ground reflectivity in January, July, and August 2022 that need to be calculated and interpolated in this example are shown in Figure 6 and Figure 7 . The calculation formula for the second surface reflectivity is:

[0083] ,

[0084] In the formula, Г i is the i th observed reflectivity of the adjacent grid point, n represents the total number of adjacent grid points with reflectivity data in the target year and target period, w i is the weighted weight between the target grid point and the i th adjacent grid point. The calculation formula is:

[0085] .

[0086] In the example, the CYGNSS satellite data is sampled in the form of reflection points. The linear trajectory of the observation points is the same as the CYGNSS orbit, and the specific observation position is determined by the positions of the receiver and the GNSS transmitter. While this sampling method brings many advantages, it also determines that not every grid point in the 250 m × 250 m grid has data. Therefore, it is necessary to reasonably infer the data of each grid point using the grid points with existing data. In this example, the POBI interpolation method is adopted. The data of the target area is divided into grids of 250 m × 250 m, and 8 grid cells surrounding each grid point are defined as adjacent grid points. The grid on the water surface is removed, and the global surface water dataset of the European Union Joint Research Centre is used to remove the main water bodies in the target area to reduce the influence of the reflection points on the water surface on the surrounding soil moisture. Then, a relationship model between the target grid points and the adjacent grid points is established using the surface reflectance data from 2019 to 2021 in the historical period for interpolating the data of the target year, and 250 m × 250 m surface reflectance data with the same spatial scale as the MODIS data is obtained. Specifically, for each grid point in the target area and its corresponding adjacent grid points, the reflectance data with a sampling time interval of no more than 24 h is linearly fitted using the data in the historical period of the past year or a longer sequence to obtain the relationship model y = a i x + b i , and then the intercept b and the slope a are obtained. Then, the correlation coefficient of the reflectance between the target grid point and the adjacent grid point is calculated, and the surface reflectance situation during the specific growth period and winter of rice in the target area and target year is weighted and interpolated using the intercept b , the slope a , and the correlation coefficient of the reflectance r . The obtained result is the 250 m × 250 m scale surface reflectance distribution during the two periods of the specific growth period and winter of the target year

[0087] Furthermore, step 4 includes:

[0088] Compare the slope i corresponding to the S i th grid cell in the target area with the slope threshold S 0 . If S i > S 0 , then this grid cell is classified as a mountain area; otherwise, this grid cell is classified as a non-mountain area

[0089] If this grid cell is classified as a mountain area, then compare the vegetation index NDVI with the vegetation index threshold NDVI 0the size, and the difference in reflectance Δ and the reflectance difference threshold Δ 0 the size, if NDVI > NDVI 0 , and Δ > Δ 0 , then this grid cell belongs to a paddy field, otherwise it does not; where Δ = SR 1 - SR 2 , SR 1 represents the second surface reflectance of the paddy field at a specific growth stage, and SR 2 represents the second surface reflectance of the winter paddy field;

[0090] If this grid cell is classified as non-mountainous, then compare the vegetation index NDVI with the vegetation index threshold NDVI 0 the size, and the second surface reflectance SR 1 of the paddy field at a specific growth stage and the reflectance threshold SR 0 the size, if NDVI > NDVI 0 , and SR 1 > SR 0 , then this grid cell belongs to a paddy field, otherwise it does not.

[0091] In the example, taking the Jianghuai region of Anhui as an example, including eight cities: Hefei, Wuhu, Huainan, Ma'anshan, Tongling, Anqing, Chuzhou, and Lu'an. The specific growth stage can be selected as July and August, and winter as January. The slope threshold is 1°, the vegetation index threshold NDVI 0 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 certain grid is greater than 1°, it is determined as a mountain. 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 this grid belongs to a paddy field, otherwise it does not; when the slope of a certain grid is less than 1°, it is determined as non-mountainous. If its NDVI is greater than 0.43 and the surface reflectance in July and August is greater than 159.5, then this grid belongs to a paddy field, otherwise it does not.

[0092] Since a higher surface reflectivity often represents a higher soil moisture content, by taking advantage of the flooding characteristics during the rice growth cycle, i.e., the paddy fields are in a long-term soaking period during specific growth stages, thresholds are set for the interpolated surface reflectivity based on the pre-investigation samples, and grid points with values less than the thresholds are excluded; due to the influence of terrain, vegetation cover, radiation, etc., the seasonal soil moisture differences in mountainous areas are often large, so the difference in soil moisture between specific growth stages and the dry period in winter is used as an index for identification; also considering the influence of other water bodies, using the processed NDVI data, grid points with NDVI values less than the threshold are regarded as ordinary water bodies, and finally the identification of rice is achieved.

[0093] Further, step 5 includes:

[0094] Multiply the number of grid cells corresponding to the selected paddy fields by the scale of each grid cell to obtain the planting area of rice. In the example, the scale is 250 meters × 250 meters.

[0095] In this example, the rice field data from the third national land survey of each city in Anhui Province is used as the true value for verification, and the relative error is used as the measurement index. After preliminary experiments and on-site inspections, as an optimization, NDVI 0 = 0.43, SR = 159.5, Δ = 3.9 are selected as the thresholds for rice identification and area extraction. Finally, the average relative error is 4.45%. The relative errors of each city are shown in Table 1, and the rice identification results are as Figure 8 shown.

[0096] Table 1

[0097]

[0098] It can be seen from the results that the result accuracy can basically be controlled within about 5%. The total accuracy of the eight cities in the Jianghuai region of Anhui can reach within 1%. Among them, the accuracy of Ma'anshan is relatively poor, probably because its planting area is relatively small and thus more sensitive. All in all, 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, divide the target area into grids, establish the relationship between the reflectivity of each grid and the adjacent grids, perform interpolation operation on the first surface reflectivity, and obtain the second surface reflectivity; the details are as follows: 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 x of each grid point in the target area is compared with the reflectivity y of the i-th adjacent grid point. i x+b i Linear regression fitting, get the intercept b i With slope a i , where the sampling time interval between any x and its 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 i-th observed reflectivity of the target grid point; 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, Г i is the observed reflectivity of the ith adjacent grid point, n is the total number of adjacent grid points with reflectivity data in the target year and target period, and w i is the weighted weight between the target grid point and the i-th adjacent grid point, and the calculation formula is: 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: NDVI=(r NIR -r R ) / (ρ NIR +r R ) In the formula, ρ R is the reflectivity of 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, and w and v represent the 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 4 includes: Compare the slope S corresponding to the i-th grid cell in the target area i and the size of the slope threshold S0. If S i >S0, 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; where Δ=Δ1-Δ2, Δ1 represents the second surface reflectance of a rice field in a specific growth period, and Δ2 represents the second surface reflectance of a 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.

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 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