A method for estimating crop water stress

Abstract

The present invention discloses a method for estimating crop water stress, comprising the following steps: obtaining remote sensing data and preprocessing: obtaining multispectral remote sensing images with less coverage of thin clouds and meteorological data of the corresponding study area, and performing preprocessing; obtaining actual evapotranspiration based on the NP method: using a non-parametric evapotranspiration inversion method to obtain actual evapotranspiration, and potential evapotranspiration is obtained through an improved Penman-Montes combination; using the actual evapotranspiration and potential evapotranspiration obtained as input variables, constructing a new crop water index, and monitoring the drought situation in the study area; using the actual evapotranspiration data measured on the surface to directly verify the RS-NP method, using the measured soil moisture data on the surface to calculate the correlation with the new crop water stress index, and performing accuracy assessment and cross-validation on the results of the RS-PM method. The present invention constructs the NM-CWSI method based on the RS-NP method and the RS-PM method, which can quickly and accurately monitor drought changes in areas with complex surfaces.

Description

Technical Field

[0001] The invention relates to a remote sensing estimation method, in particular to a method for estimating crop water stress. Background Art

[0002] Agricultural drought disasters are the most frequent, longest-lasting and most extensive natural disasters in the world. Against the backdrop of global warming, extreme weather occurs frequently, and the frequency and scope of agricultural drought disasters are increasing. This is one of the most serious problems currently facing mankind. For a long time, the traditional way of monitoring agricultural droughts has been to use observation records from meteorological stations for data statistical analysis. However, the number of meteorological stations is limited and cannot cover all areas, so their representativeness is poor. Satellite remote sensing technology has the advantages of wide observation range, large amount of information, fast speed, good real-time performance and strong dynamics, which can well realize the monitoring of drought areas.

[0003] The existing methods for monitoring drought mainly use the temperature vegetation drought index method (TVDI) and the traditional crop stress index estimation method (CWSI). The TVDI method is derived based on the spatial relationship between surface temperature and vegetation index. This index is closely related to soil moisture content, can indicate the degree of vegetation water stress, and better describe agricultural drought conditions. However, the key to this method lies in the determination of the boundary, but the boundary determination is greatly affected by the scatter plot, so the accuracy of the calculation result is unstable. The traditional CWSI method believes that there is a linear relationship between the difference between the canopy temperature of crops under potential evaporation conditions and the atmospheric temperature and the air saturated water vapor pressure difference. The theoretical basis of this method is the canopy energy balance single-layer model. The canopy temperature obtained by remote sensing inversion is combined with meteorological data to comprehensively reflect the deep soil moisture content. It has clear physical meaning and relatively high accuracy. However, the complex calculation process, more input data, poor applicability, and the limitations of meteorological station data limit the real-time application of this model.

[0004] In general, the boundary determination of existing drought monitoring methods is greatly affected by scatter plots, the accuracy of calculation results is unstable, the model is complex, the applicability is poor, and there are data limitations such as meteorological stations. Summary of the invention

[0005] Purpose of the invention: In order to overcome the deficiencies in the prior art, the purpose of the present invention is to provide a method for estimating crop water stress based on a non-parametric evapotranspiration inversion method (NP) that is convenient for high-precision monitoring of drought conditions in a target study area.

[0006] Technical solution: A method for estimating crop water stress according to the present invention comprises the following steps:

[0007] Step 1: Obtain remote sensing data and preprocess: obtain multispectral remote sensing images with less cloud coverage and meteorological data of the corresponding study area, and perform preprocessing;

[0008] Step 2, actual evapotranspiration is obtained based on the NP method: the actual evapotranspiration is obtained by using the nonparametric evapotranspiration inversion method, and the potential evapotranspiration is obtained by the improved Penman-Monteis combination (PM);

[0009] Step 3: Construct a new crop water stress index: Using the actual evapotranspiration and potential evapotranspiration obtained in step 2 as input variables, construct a new crop water stress index to monitor the drought conditions in the study area. The calculation formula of the new crop water stress index is as follows:

[0010]

[0011] Among them, CWSI is a new crop water stress index, ET d is the actual evapotranspiration, i.e. daily evapotranspiration, ET p is the potential evapotranspiration;

[0012] Step 4: Accuracy evaluation and cross-validation: Use actual evapotranspiration data measured on the surface to directly verify the RS-NP (remote sensing non-parametric) method, use surface measured soil moisture data to calculate the correlation with the new crop water stress index, and use MOD16-PET products to evaluate the accuracy and cross-validate the results of the RS-PM method.

[0013] Furthermore, in step 1, the preprocessing used includes image correction, mosaicking, cropping, resampling, and format conversion. The multispectral images of the original remote sensing images are corrected to make the resolution the same. The need for mosaicking is determined according to the size of the study area, and finally they are uniformly cropped to the same study area. Multispectral remote sensing images include Landsat8 data and CLDAS data. The spatial resolution of Landsat8 data is 30 meters, and the spatial resolution of CLDAS data is 8 km.

[0014] Since the NP method involves fewer parameters and avoids the problems of complex impedance parameterization calculation, difficult error source analysis and difficult model correction in the traditional inversion model, the NP method is used to obtain actual evapotranspiration and the PM method is used to obtain potential evapotranspiration. By obtaining actual evapotranspiration and potential evapotranspiration, we prepare for the construction of a new crop water stress index. In step 2, the calculation formula of the non-parametric evapotranspiration inversion method is:

[0015]

[0016] Among them, LE is the latent heat flux, which is the evapotranspiration we are studying. Rn is the net surface radiation, Gs is the soil heat flux, Ts is the surface temperature, ε is the surface emissivity, σ is the Stefan-Boltzmann constant, T0 is the near-surface atmospheric temperature, and Δ is T a The saturated water vapor pressure gradient at , γ is the dry-bulb and wet-bulb constants, and the surface net radiation is quantitatively inverted through remote sensing data.

[0017] Furthermore, the calculation formula of the saturated water vapor pressure gradient Δ is:

[0018]

[0019] Among them, T0 is the near-surface atmospheric temperature.

[0020] Furthermore, the calculation formula of the wet-dry bulb constant γ is:

[0021]

[0022] Among them, C P is the constant pressure specific heat capacity under normal pressure, P is the atmospheric pressure near the surface, ε aw is the ratio of the weights of water vapor and air molecules.

[0023] Furthermore, the calculation formula of the surface emissivity ε is:

[0024]

[0025] Among them, a λ The reflectivity band ρ λ The coefficient converted to emissivity is c, which is the reflectivity band ρ λ The constant that transforms into emissivity, ε v is the emissivity of the vegetation component, ε g is the emissivity of the background component, dε is the reflectivity increment caused by the hole effect caused by multi-scattering in a pixel, f is the vegetation cover, NDVI is the normalized difference vegetation index, NDVIs is the NDVI value of the bare soil pixel, and NDVIv is the NDVI value of the vegetation pixel.

[0026] Furthermore, the calculation formula of vegetation coverage f is:

[0027] f=[(NDVI-NDVI s ) / (NDVI v -NDVI s )] 2 .

[0028] Furthermore, in step 4, the indicators of accuracy assessment and cross-validation include determination coefficient, mean error, root mean square error and relative error percentage.

[0029] Beneficial effects: Compared with the prior art, the present invention has the following significant features:

[0030] 1. It can improve the accuracy of evapotranspiration inversion. The inversion accuracy of RS-NP actual evapotranspiration and RS-PM potential evapotranspiration based on Landsat8 data is high. The correlation coefficient R between the actual evapotranspiration inverted by RS-NP and the measured data at the site is 2 The average error is 0.968 and 8.29W / m 2 , the root mean square error is 21.14W / m 2 The relative error is 5.12%. The correlation coefficient R2 between the potential evapotranspiration inverted by RS-PM and the measured data at the site is 0.872, and the average error is -37.06 W / m 2 , the root mean square error is 83.65W / m 2 , the relative error is 10.61%;

[0031] 2. It can overcome the limitations of traditional drought monitoring methods. The NM-CWSI method is constructed based on the RS-NP method and the RS-PM method, avoiding the problems of complex calculations and poor applicability, so that the method of the present invention can quickly and accurately monitor drought changes in areas with complex surfaces. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 It is the technical roadmap of the present invention;

[0033] Figure 2 is an inversion overall situation diagram of actual evapotranspiration and potential evapotranspiration of the present invention, wherein: (a) actual evapotranspiration, (b) potential evapotranspiration;

[0034] Figure 3 1 is a quantitative comparison diagram of the results of different methods of the present invention, wherein: (a) false color image, (b) NM-CWSI method, (c) theoretical CWSI method, and (d) TVDI method;

[0035] Figure 4 is a scatter plot of the results of different methods of the present invention and soil moisture data, including: (a) NM-CWSI method, (b) theoretical CWSI method, and (c) TVDI method;

[0036] Figure 5 It is the fitting diagram of the NM-CWSI method and the actual soil moisture content at each station of different underlying surfaces of the present invention, among which, (a) forest station, (b) farmland station, (c) grassland station;

[0037] Figure 6 It is the fitting diagram of the NM-CWSI method and the actual soil moisture content of each station at different elevations of the present invention, among which, (a) farmland station, (b) Daman station, and (c) Huangzangsi station. DETAILED DESCRIPTION

[0038] A new crop stress index was calculated for the Heihe River Basin. The method for estimating crop water stress is as follows: Figure 1 As shown, the specific steps include:

[0039] Step 1: Obtain Landsat 8 multispectral remote sensing images and CLADS meteorological data of the Heihe River Basin and perform preprocessing. The preprocessing mainly includes image format conversion, correction, mosaicking and cropping. The required data include Landsat8 data and CLDAS data. Landsat8 data is obtained from the website: https: / / glovis.usgs.gov / , and its spatial resolution is 30 meters. CLDAS data is a land surface data assimilation product, produced and released by the National Meteorological Information Center, with a spatial resolution of about 8km. It is necessary to perform atmospheric correction on Landsat8 data and resample CLDAS meteorological data to the same 30-meter resolution as Landsat8 data. Use the remote sensing software ENVI to perform a series of preprocessing on remote sensing images, correct the multispectral images of the original remote sensing images, determine whether mosaicking is needed according to the size of the study area, and finally crop the multispectral remote sensing images and CLADS images to the same study area.

[0040] Step 2: Obtain actual evapotranspiration and potential evapotranspiration based on the NP method. Since the calculation of the crop stress index requires actual evapotranspiration, in order to build a better model to indicate crop water stress, the RS-NP method is used to obtain actual evapotranspiration based on Landsat 8 data. The formula and principle of the RS-NP method are as follows:

[0041]

[0042] In the formula, LE is the latent heat flux, which is the evapotranspiration we are studying. Rn is the net surface radiation, Gs is the soil heat flux, Ts is the surface temperature, ε is the surface emissivity, σ is the Stefan-Boltzmann constant, T0 is the near-surface atmospheric temperature, and Δ is T a The saturated water vapor pressure gradient at γ is the dry-bulb constant, and the net radiation of the surface is quantitatively inverted through remote sensing data. The saturated water vapor pressure gradient and the dry-bulb constant are usually expressed as:

[0043]

[0044]

[0045] In the formula, C P is the constant pressure specific heat capacity under normal pressure (1.013×10 -3 Mk -1 K-1 ), P is the atmospheric pressure near the surface, ε aw It is the ratio of the weight of water vapor to that of air molecules (approximately 0.622).

[0046] The calculation formula of the surface emissivity ε is as follows:

[0047]

[0048] Among them, a λ and c are the reflectivity band ρ λ The coefficient and constant that transforms to emissivity, ε v and ε g are the emissivity of the vegetation component and the emissivity of the background component, respectively. dε is the reflectivity increment caused by the hole effect caused by multi-scattering in a pixel. f is the vegetation coverage calculation formula as follows:

[0049] f=((NDVI-NDVI s ) / (NDVI v -NDVI s )) 2

[0050] Among them, NDVIs and NDVIv represent the NDVI values ​​of bare soil pixels and vegetation pixels, respectively.

[0051] Based on the surface energy balance, the surface net radiation can be expressed by the surface shortwave radiation and the surface longwave radiation:

[0052] R n =(1-α)R sd +R ld -εσT3 4

[0053] Where α is the surface albedo, R sd is the surface shortwave downward radiation, R ld is the surface longwave downward radiation

[0054]

[0055] I0 is the solar constant (about 1367W / m 2 ), d is the distance from the sun to the earth, and θ is the solar incidence angle.

[0056] τ sw is the atmospheric transmittance, which is calculated as follows:

[0057]

[0058] Kt is the correction coefficient (0<Kt≤1.0), Kt=1 under clear sky conditions, Kt-0.5 under cloudy conditions, W is the atmospheric water vapor content, and the calculation formula is as follows:

[0059] W=0.14e0·P+2.1

[0060] Surface longwave downward radiation R ld The calculation formula is as follows:

[0061] R ld =σε a T0 4

[0062] Where T0 is provided by CLDAS, ε a is the atmospheric emissivity near the surface, which can be calculated by the following formula:

[0063]

[0064] Where e0 is the near-surface water pressure, which can be calculated by the following formula:

[0065]

[0066] Where P is the near-surface atmospheric pressure provided by CLDAS, and q is the specific humidity also from CLDAS data.

[0067] Another parameter required to calculate the available energy is the soil heat flux, which is calculated as follows:

[0068]

[0069] Where NDVI is the normalized difference vegetation index.

[0070] The above are all the formulas of RS-NP model based on Landsat data, in which potential evapotranspiration is obtained by improved Penman-Montes combination.

[0071] The actual evapotranspiration and potential evapotranspiration were obtained based on the RS-NP method and RS-PM method, and the determination coefficient (R 2 ), mean error (ME), root mean square error (RMSE) and relative error percentage (RE) to evaluate its accuracy. Figure 2(a) is an overall picture of actual evapotranspiration inversion. The horizontal axis in the figure is the actual evapotranspiration value obtained by inversion using the RS-NP model, and the vertical axis is the actual evapotranspiration value observed on the surface. The trend line of the scatter plot is close to the 1:1 line, and there is a slight underestimation. The actual evapotranspiration inversion results are generally close to the actual observation values ​​of the experimental area stations, and the correlation coefficient R 2The average error of the inversion is 8.29 W / m 2 , the root mean square error is 21.14W / m 2 The relative error is 5.12%, which is much lower than the 15% to 40% error of mainstream evapotranspiration inversion estimates. Figure 2 (b) is the overall situation of potential evapotranspiration inversion. The horizontal axis of the figure is the potential evapotranspiration value obtained by using the RS-PM model, and the vertical axis is the potential evapotranspiration value of the MOD16 product. The trend line of the scatter plot is close to the 1:1 line, which is slightly underestimated. The overall potential evapotranspiration inversion results are close to the actual observation values ​​of the experimental area stations, and the correlation coefficient R 2 The average error of the inversion is -37.06 W / m 2 , the root mean square error is 83.65W / m 2 The relative error is 10.61%, indicating that the potential evapotranspiration inversion results of the RS-PM model based on Landsat8 have a higher simulation accuracy than the MOD16-PET product and are more consistent with the actual situation. It can be used to build the CWSI model.

[0072] Step 3: Construct a new crop water stress index CWSI, the formula of which is as follows:

[0073]

[0074] In the formula, ET d is the actual evapotranspiration, i.e. daily evapotranspiration (mm / d); RT p is the potential evapotranspiration (mm / d).

[0075] The essence of CWSI is to convert the soil moisture content at the root of crops through the temperature of the vegetation canopy by inverting the ratio of actual daily evapotranspiration to potential evapotranspiration. According to the definition of CWSI, its value range is 0 to 1.

[0076] What is obtained from step 3 is the instantaneous actual evapotranspiration and potential evapotranspiration, but in order to calculate the CWSI index, the instantaneous evapotranspiration must be converted into daily evapotranspiration. The present invention uses the evaporation ratio method to expand the time range of surface evapotranspiration. The method uses the following formula:

[0077]

[0078] Where LE is instantaneous evapotranspiration, R n is the instantaneous shortwave net radiation, R nd is the net radiation during the day, LE dis the daily evapotranspiration during the day. Since the evapotranspiration at night is very small, the daily evapotranspiration during the day is considered to be the daily evapotranspiration for the whole day on a long time scale. nd It can be calculated by the following formula:

[0079]

[0080] In the formula, t 过境 , t 日出 , t 日落 They are the time when the satellite passes through the study area, the sunrise time in the study area, and the sunset time in the study area.

[0081] The LE obtained in the above formula is actually the latent heat flux, in W / m 2 , and ET is the evapotranspiration value, the unit is mm, so it is necessary to convert LE to get ET, thereby obtaining the potential evapotranspiration and actual evapotranspiration data on a daily scale, and putting them into the formula of the CWSI index to obtain the new crop water stress index.

[0082] Figure 3 It is a quantitative comparison chart of other traditional methods and the new crop stress index model proposed in the present invention. As can be seen from Figure 3, the dark blue area is an area with sufficient water and no irrigation is needed temporarily. The light blue area indicates that the crops are suffering from a relatively light degree of water stress and can be irrigated in moderation. The light yellow area indicates that the crops are suffering from a relatively heavy degree of water stress and need timely irrigation. The two CWSI models are consistent and clear in distinguishing between farmland vegetation areas and non-vegetation areas. The NM-CWSI model has clear boundaries for different degrees of stress, while the CWSI theoretical model is less sensitive to small-area differences. The TVDI has blurred boundaries for non-vegetation areas, and there is a reverse trend in some boundary areas. In summary, the NM-CWSI model can better indicate the areas and degrees of crop water stress, and provide corresponding technical support for precision agriculture and smart agriculture.

[0083] Step 4: Accuracy evaluation and cross-validation. The quality of remote sensing inversion results is evaluated qualitatively and quantitatively. Qualitative evaluation refers to whether the crop water stress index inversion result map is consistent with the actual true color image characteristics. Quantitative evaluation refers to the evaluation of indicators calculated by some mathematical statistical formulas. Then select the coefficient of determination (R 2 ), mean error (ME), root mean square error (RMSE) and relative error percentage (RE) were used to evaluate the actual evapotranspiration, potential evapotranspiration and crop water stress index results.

[0084] Coefficient of determination R 2 It is used to evaluate the consistency between the inversion value and the surface observation value. 2The value of is between (0, 1), when R 2 The closer the value of is to 1, the greater the correlation between the inversion value and the measured value, and the better the result. The formula is as follows:

[0085]

[0086] Where Y and Y 真 are the inverted value and the true value respectively, is the mean of the observed values, and N is the number of pixels.

[0087] The mean error ME reflects the average difference between the inversion result and the surface observation. The root mean square error RMSE is very sensitive to the large or small errors of the observation value and is used to measure the degree of dispersion between the calculated value and the observed value. The relative error RE is the ratio of the mean error to the observed value, expressed as a percentage. Generally speaking, the relative error can better reflect the reliability of the inversion result.

[0088]

[0089]

[0090]

[0091] Where Y is the inversion value, Y 真 is the site observation value, and N is the remote sensing evapotranspiration inversion value or the number of site observation values.

[0092] The correlation between the measured soil moisture data on the surface and the NM-CWMI results is calculated. Figure 4 Shown is a scatter plot of the results of the three methods and soil moisture data. Figure 4 (a) is the fitting of the inversion value of the improved NM-CWSI model and soil moisture content. The correlation between the soil moisture content in the Heihe River Basin and the NM-CWSI model is -0.8, and it has passed the significance test with a confidence level of 0.05. This shows that the NM-CWSI model can correctly describe the spatial distribution of the degree of crop water stress with good accuracy. Figure 4 (b) is the fitting of the inversion value of the CWSI theoretical model and the actual soil moisture content. The correlation between the soil moisture in the Heihe River Basin and the CWSI theoretical model is -0.691, and it has passed the significance test at a confidence level of 0.05. This shows that the CWSI theoretical model can also correctly describe the spatial distribution of the degree of crop water stress, but the CWSI inverted by the NM-CWSI model has a higher accuracy in simulating crop water stress in the Heihe River Basin than the CWSI theoretical model. Figure 4(c) The TVDI model inversion value and soil moisture fit, the correlation between soil moisture in the Heihe River Basin and the TVDI model is -0.6459. This shows that the TVDI model can better describe the spatial distribution of the degree of crop water stress, but the CWSI obtained by the NM-CWSI model has a simulation accuracy of -0.834 for crop water stress in the Heihe River Basin, which is higher than the TVDI model. Therefore, overall, the NM-CWSI model can better indicate the degree of crop water stress.

[0093] Step 5: Applicability analysis of the new crop stress index. The applicability of the new crop stress index was explored from two aspects: different underlying surfaces and elevations. Figure 5 is the fitting of the NM-CWSI model and the actual soil moisture content at each station on different underlying surfaces. 2 All of them are above -0.85, indicating that the crop water stress index inverted by the NM-CWSI model has a high correlation with the measured value of soil surface moisture. Figure 6 It is the fitting of NM-CWSI model and actual soil moisture content at farmland sites at different altitudes. 2 From high to low, they are farmland, forest, and grassland. The R 2 Reached -0.91. The correlation between the crop water stress inversion results based on the NM-CWSI model and the measured soil moisture content of the three types of farmland stations at three altitudes is above -0.7. Among them, the inversion result of Daman Station has the highest correlation, while the inversion results of the farmland station and Huangzangsi Station have a lower correlation, and the inversion result of Huangzangsi Station has the lowest accuracy. The altitudes of the three stations are Huangzangsi Station, Daman Station, and Farmland Station from high to low. It can be seen that the altitude has little effect on the inversion accuracy of the NM-CWSI crop water stress estimation model.

[0094] Overall, whether from the perspective of qualitative analysis or quantitative analysis, compared with other methods for monitoring agricultural drought, the new crop stress index proposed in the present invention can more effectively indicate the degree of crop water stress.

Claims

1. A method for estimating crop water stress, characterized in that: The following steps are involved: Step 1: obtain multispectral remote sensing images with less thin cloud coverage and meteorological data of the corresponding study area, and perform preprocessing; Step 2, the actual evapotranspiration is obtained by using the nonparametric evapotranspiration inversion method, and the potential evapotranspiration is obtained by the improved Penman-Montes combination; Step 3: Using the actual evapotranspiration and potential evapotranspiration obtained in step 2 as input variables, a new crop water stress index was constructed to monitor the drought conditions in the study area. The calculation formula of the new crop water stress index is as follows: Among them, CWSI is a new crop water stress index. is the actual evapotranspiration, i.e. the daily evapotranspiration, is the potential evapotranspiration; Step 4: Use the actual evapotranspiration data measured on the surface to directly verify the RS-NP method, use the measured soil moisture data on the surface to calculate the correlation with the new crop water stress index, and conduct accuracy assessment and cross-validation on the results of the RS-PM method; In the step 2, the calculation formula of the non-parametric evapotranspiration inversion method is: in, LE is the latent heat flux, in W / m 2 , and ET is the evapotranspiration value, the unit is mm, so it is necessary to convert LE to get ET. R is the net surface radiation, Gs is the soil heat flux, Ts is the surface temperature, is the surface emissivity, is the Stefan-Boltzmann constant, is the near-surface atmospheric temperature, Δ is The saturated water vapor pressure gradient at , γ is the dry-bulb and wet-bulb constants, and the surface net radiation is quantitatively inverted through remote sensing data.

2. A method for estimating crop water stress according to claim 1, characterized in that: In the step 1, the preprocessing used includes image correction, mosaicking, cropping, resampling, and format conversion. The multispectral image of the original remote sensing image is corrected to make the resolution the same. The need for mosaicking is determined based on the size of the study area, and finally the images are uniformly cropped to the same study area.

3. A method for estimating crop water stress according to claim 1, characterized in that: In the step 1, the multispectral remote sensing image is Landsat8 data, and the meteorological data is CLDAS data.

4. A method for estimating crop water stress according to claim 3, characterized in that: In the step 1, the spatial resolution of the Landsat8 data is 30 meters, and the spatial resolution of the CLDAS data is 8 km.

5. The method for estimating crop water stress according to claim 1, characterized in that: The calculation formula of the saturated water vapor pressure gradient Δ is: in, is the near-surface atmospheric temperature.

6. A method for estimating crop water stress according to claim 1, characterized in that: The calculation formula of the dry-bulb constant γ is: in, is the constant pressure specific heat capacity under normal pressure, P is the atmospheric pressure near the surface, is the ratio of the weights of water vapor and air molecules.

7. A method for estimating crop water stress according to claim 1, characterized in that: The surface emissivity The calculation formula is: Among them, a λ The reflectivity band ρ λ The coefficient of conversion into emissivity is c, which is the reflectivity band ρ λ The constant that transforms into emissivity, ε v is the emissivity of the vegetation component, ε g is the background component emissivity, dε is the reflectivity increment caused by the hole effect caused by multiple scattering in a pixel, f is the vegetation coverage, NDVI is the normalized difference vegetation index, NDVIs is the NDVI value of the bare soil pixel, and NDVIv is the NDVI value of the vegetation pixel.

8. A method for estimating crop water stress according to claim 7, characterized in that: The vegetation coverage f The calculation formula is: 。 9. A method for estimating crop water stress according to claim 1, characterized in that: In step 4, the indicators of accuracy assessment and cross-validation include determination coefficient, mean error, root mean square error and relative error percentage.

Images