A farmland root layer soil salt content inversion method based on remote sensing data fusion and application

Abstract

The present application relates to a kind of farmland root layer soil salt content inversion method and application based on remote sensing data fusion.Different from other remote sensing inversion method can only roughly estimate surface soil salt content, the present application is based on the influence mechanism of crop growth period root layer soil salt on transpiration and evapotranspiration, using the relationship between relative evapotranspiration rate and soil salt stress correction coefficient to establish root layer soil salt content remote sensing inversion model, based on remote sensing image data fusion with different temporal and spatial resolution, inversion obtains the spatial distribution rule of root layer soil salt on regional scale, and combines with the temporal and spatial evolution rule of remote sensing image data analysis for many years.The application of the present application technology can enrich soil water and salt transport theory, perfect remote sensing monitoring method, also can help to deeply understand the temporal and spatial evolution rule of regional scale farmland evapotranspiration rate and root layer soil salt and its driving mechanism, for optimizing planting structure and water and soil resource allocation, promote the sustainable development of saline-alkali agriculture and provide new means and scientific basis for protecting ecological environment.

Description

Technical Field

[0001] This invention relates to farmland root zone soil salinity inversion technology, and more particularly to a method and application for farmland root zone soil salinity inversion based on remote sensing data fusion. Background Technology

[0002] Understanding the spatiotemporal evolution of root zone soil salinity at the regional scale is crucial for promoting the efficient use of water and soil resources, especially for the sustainable development of saline-alkali agriculture. However, due to the spatial variability of soil, techniques for monitoring and assessing root zone soil salinity at the point-scale in farmland are difficult to apply to the regional scale.

[0003] In recent years, remote sensing has gradually become a major means of monitoring and evaluating the spatiotemporal evolution of regional soil water and salt due to its ability to acquire large-area, multi-band, and multi-temporal information. Based on methods such as spectral index, multi-parameter coupling, and machine learning, a large number of remote sensing inversion models have been successfully developed (such as Chinese invention patents: CN11 1783288A, CN115308386A).

[0004] However, existing remote sensing inversion models can only monitor and assess the salinity of surface soil, falling far short of the actual needs of agricultural production to understand the salinity of the entire root zone (root layer) soil. For example, in arid and semi-arid regions, water scarcity and soil salinization are both prominent issues. It is necessary to meet the water requirements of crops to a certain extent while mitigating the harmful effects of salinity and preventing further soil salinization. Therefore, higher demands must be placed on agricultural water resource utilization, namely, achieving efficient water use at both the farmland and regional scales. However, a prerequisite for this is understanding the spatiotemporal evolution of soil salinity in the crop root zone. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention aims to provide a method for inverting farmland root zone soil salinity based on remote sensing data fusion. This method combines a surface energy balance model and an adaptive remote sensing image spatiotemporal fusion model, proposing a fusion method with high spatiotemporal resolution for the relative evapotranspiration rate Λ. It utilizes the approximate relationship between the relative evapotranspiration rate Λ during the peak crop growth period (when sufficient irrigation is generally implemented to achieve substantial yields, i.e., no water stress) and the soil salinity stress correction coefficient to invert regional-scale crop root zone soil salinity. This not only enriches the theory of soil water and salt transport and improves remote sensing monitoring methods, but also helps to deepen the understanding of the spatiotemporal evolution patterns and driving mechanisms of regional-scale farmland evapotranspiration rate and root zone soil salinity. This provides new means and scientific basis for optimizing crop planting structure and water and soil resource allocation, promoting the sustainable development of saline-alkali agriculture, and protecting the ecological environment.

[0006] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0007] A method for retrieving farmland root zone soil salinity based on remote sensing data fusion, characterized by the following steps:

[0008] Step 1: Establish a quantitative relationship between crop relative transpiration rate and root zone soil salinity based on field trial data:

[0009]

[0010] In the formula: T a With T p These represent the actual and potential transpiration rates of the crop, respectively. Average soil salinity in the root zone: This is the soil salinity stress correction factor; This is the critical value for crop salt tolerance. When the soil salt content is higher than this critical value, crops will suffer from salt stress. τ represents the soil salinity corresponding to the permanent wilting of crops under severe stress; τ is the fitting parameter.

[0011] Step 2: Based on the approximate relationship between relative evapotranspiration rate and relative transpiration rate during the peak growth period of crops, the relationship between relative evapotranspiration rate Λ and root zone soil salinity is further established using equation (1). The relationship between the two, namely the remote sensing inversion model of root zone soil salinity:

[0012]

[0013] Step 3: Determine the study period during the peak growth period of crops, acquire remote sensing image data of satellites L and M in the study area during the study period, and use the surface energy balance model to invert the spatial distribution pattern of Λ based on each remote sensing image.

[0014] Step 4: Referring to the adaptive remote sensing image spatiotemporal fusion model, the Λ obtained in Step 3 is fused and calculated to obtain the high-resolution spatial distribution pattern of Λ every day during the study period.

[0015] Step 5: Based on the high-resolution spatial distribution pattern of Λ for each day during the study period obtained in Step 4, obtain the average value of Λ for the study period. The high-resolution spatial distribution pattern is obtained, and the result is obtained by inversion according to equation (2). Spatial distribution patterns.

[0016] Satellite L represents a satellite with high spatial resolution and low temporal resolution, while satellite M represents a satellite with high temporal resolution and low spatial resolution.

[0017] Based on the above scheme, the study period is one or two revisit cycles of a certain satellite L during the peak growth period of crops.

[0018] Based on the above scheme, the specific steps of step 3 are as follows:

[0019] Step 3-1: Set the first and last days of the study period as reference times t. m and t n , get t m and t n The required surface parameters were obtained by inverting each remote sensing image from satellite L at any given time and from satellite M every day during the study period. These parameters mainly include surface albedo α0, surface emissivity ε0, normalized difference vegetation index (NDVI), surface temperature T0, and solar zenith angle (SZA).

[0020] Step 3-2: Estimate the net radiation flux R using the surface parameters obtained in Step 3-1. n :

[0021]

[0022] In the formula: R sd For downlink shortwave radiation; R ld For downward longwave radiation; σ is the Stefan-Boltzmann constant;

[0023] Step 3-3: Calculate vegetation cover f using the NDVI obtained in Step 3-1. c :

[0024]

[0025] In the formula: NDVI soil NDVI value for bare ground; NDVI veg The NDVI value is the value when the vegetation is completely covered.

[0026] Step 3-4, based on T0 and f obtained in steps 3-1 and 3-3 c In addition to meteorological data such as atmospheric pressure, wind speed, and temperature, the frictional wind speed u is iteratively solved. * , heat flux H and Obukhov stability length L;

[0027] Step 3-5: Based on the parameters obtained in steps 3-1, 3-2, 3-3, and 3-4 above, t is obtained by inversion according to the surface energy balance model shown in the following formula. m and t n The reference time is based on the Λ value of remote sensing images from satellites L and M, i.e., Λ L (t m ), Λ L (tn ) and Λ M (t m A) M (t n ), and a certain predicted time t p The Λ value based on satellite M remote sensing imagery, i.e., Λ M (t p ):

[0028]

[0029] In the formula: λ is the latent heat of vaporization of water; ρ w R is the density of water. n G is the net radiation flux; H is the soil heat flux; ... H is the dry =R n -G represents the sensible heat flux under extremely dry conditions; H wet For sensible heat flux under fully humid conditions; ET wet =(H dry -H wet ) / (λρ w ), which represents the evapotranspiration intensity under fully moistened conditions;

[0030] The above t p The time periods other than the reference time are considered within the study period, i.e., times other than t within the study period. m and t n Any day outside of [the specified date].

[0031] Based on the above scheme, the specific steps of step 4 are as follows:

[0032] Step 4-1, using reference time t m and t n Λ obtained from satellite L remote sensing imagery L (t m ) and Λ L (t n The distribution map selects similar pixels for each target pixel within the study area. For any target pixel, the intersection of similar pixels at two reference times is selected as the final similar pixel, specifically:

[0033] For each reference time, Λ has high resolution L , that is, Λ L (t m ) and Λ L (t n The distribution map uses a target pixel as the center and sets a search box of a certain size within its neighborhood. It then compares the difference between each pixel within the box and the target pixel. If the difference is less than a threshold, the pixel is considered a similar pixel to the target pixel. Then, t is selected... m and t nThe intersection of similar pixels of the target pixel at any given time is taken as the final similar pixel of the target pixel;

[0034] The above threshold is determined by the following formula:

[0035] |Λ L (x i y i , t k )-Λ L (x w / 2 y w / 2 , t k )|≤σ×2 / e(k=m or n) (6)

[0036] In the formula: w is the width of the search box; (x w / 2 y w / 2 (x) is the center of the target pixel; i y i ) represents the position of the i-th similar pixel, i = 1, 2, ..., N, where N is the number of similar pixels, including the center pixel; σ is the position of Λ L The standard deviation of all cell values; e is the estimated number of land cover types;

[0037] Step 4-2, calculate the reference time Λ L That is, Λ L (t m ) and Λ L (t n The spatial distance d between the target pixel and the i-th similar pixel in the image i :

[0038]

[0039] Calculate the i-th similar pixel Λ at the reference time L (Λ L (t m ) and Λ L (t n )) and Λ M (Λ M (t m ) and Λ M (t n The correlation coefficient R between the two is... i :

[0040]

[0041] In the formula: Λ L (x i y i , t m ) and Λ L (x i y i , tn ) represent t respectively m and t n Λ of time L That is, Λ L (t m ) and Λ L (t n The value of the i-th similar pixel in the image; Λ M (x i y i , t m ) and Λ M (x i y i , t n ) represent t respectively m and t n Λ of time M That is, Λ M (t m ) and Λ M (t n The value of the i-th similar pixel in the image; μ L and σ L Representing Λ L (x i y i , t m ) and Λ L (x i y i , t n The mean and standard deviation of μ. M and σ M Representing Λ M (x i y i , t m ) and Λ M (x i y i , t n The mean and standard deviation of ().

[0042] Based on the obtained d i and R i Calculate the comprehensive index D of the i-th similar pixel to the target pixel. i :

[0043] D i =(1-R) i )×d i (9)

[0044] D i The normalized reciprocal of W is used as the weight W of the i-th similar pixel. i :

[0045]

[0046] Step 4-3, calculate Λ L (Λ L (t m ) and Λ L (t n )) and Λ M (Λ M (t m ) and Λ M (t n The conversion coefficient V between the i-th similar pixel values ​​in the image i :

[0047]

[0048] Step 4-4, combined with W i and V i According to t m and t n Λ of time L (t m ), Λ L (t n ) and Λ M (t m ), Λ M (t n Distribution plot and t p Λ of time M (t p Distribution plot, predict t p The high spatial resolution Λ at any given time allows us to obtain the distribution pattern of Λ with high spatiotemporal resolution for each day within the study period:

[0049] Λ(x w / 2 y w / 2 , t p ) = T m ×Λ m (x w / 2 y w / 2 , t p )+T n ×Λ n (x w / 2 y w / 2 , t p (12)

[0050] in

[0051] In the formula: k = m or n; Λ(x w / 2 y w / 2 , t p ) represents the final t obtained after fusion. p The Λ value has high spatiotemporal resolution at any given time; when k = m, Λ m (xw / 2 y w / 2 , t p ) is based on t m Time-based remote sensing images from satellites L and M, and t p The time-based prediction of t by satellite M remote sensing images p The value of Λ at time k; when k = n, Λ n (x w / 2 y w / 2 , t p ) is based on t n Time-based remote sensing images from satellites L and M, and t p The time-based prediction of t by satellite M remote sensing images p The value of Λ at time T; m and T n t m and t n Time-based remote sensing imagery relative to t p Time weighting factor at any given moment;

[0052] The above T m With T n Resampled Λ M Image at reference time t m With t n and the predicted time t p The range of change between them is determined by the following formula:

[0053]

[0054] In the formula: k = m or n.

[0055] Another objective of this invention is to provide an application of a method for retrieving farmland root zone soil salinity based on remote sensing data fusion.

[0056] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0057] An application of a method for retrieving farmland root zone soil salinity based on remote sensing data fusion, characterized by:

[0058] The inversion method can be applied to single crop planting structures or multiple crop planting structures within the study area. However, for each major crop, a remote sensing inversion model for root soil salinity should be established separately, i.e., the soil salinity stress correction coefficient of equation (2) should be optimized and determined. The fitting parameter τ in the equation.

[0059] An application of a method for retrieving farmland root zone soil salinity based on remote sensing data fusion, characterized by:

[0060] The inversion method can be used to obtain the spatiotemporal evolution of root soil salinity over many years based on remote sensing image data of crop growth peak over many years.

[0061] The present invention discloses a method for inverting farmland root zone soil salinity based on remote sensing data fusion, and its application has the following beneficial effects:

[0062] Compared to existing remote sensing inversion models that can only roughly estimate the salinity of topsoil, this invention fully considers the impact mechanism of farmland root zone soil salinity on crop transpiration and evapotranspiration, as well as water deficit. This allows for the inversion of salinity across the entire root zone soil, and by combining multi-year remote sensing image data, it further reveals the spatiotemporal distribution patterns of farmland root zone soil salinity at the regional scale. Furthermore, by combining a surface energy balance model with an adaptive remote sensing image spatiotemporal fusion model, this invention obtains a relative evapotranspiration rate Λ with high spatiotemporal resolution, thus meeting the requirements for crop evapotranspiration estimation due to the high spatiotemporal heterogeneity of regional farmland and improving the accuracy of farmland root zone soil salinity inversion. Attached Figure Description

[0063] The present invention includes the following figures:

[0064] Figure 1 Spatial distribution patterns of soil salinity (SSC) in the root zone (0-80cm) of drip-irrigated cotton fields under mulch in the Manas River Basin from 2000 to 2020;

[0065] Figure 2 Location maps of soil sampling points in the Manas River Basin, Xinjiang, from 2018 to 2020: a. Location of Anjihai and Mosuowan irrigation areas and land use distribution in the basin in 2020; b. Distribution of 4km-scale sampling points in Anjihai irrigation area (2018–2020); c. Distribution of 4km-scale sampling points in Mosuowan irrigation area (2020); d. Distribution of 500m-scale sampling points in Anjihai irrigation area (2018–2019); e. Distribution of 100m-scale sampling points in Anjihai irrigation area (2018–2019).

[0066] Figure 3 Comparison between the inverted and measured values ​​of soil salinity in the root zone (0-80cm) of drip-irrigated cotton fields under mulch film in the Anjihai (AJH) and Mosuowan (MSW) irrigation areas from 2018 to early July 2020 (a. Error analysis statistics table; b. 1:1 comparison chart). Detailed Implementation

[0067] The present invention will be further described in detail below with reference to the accompanying drawings.

[0068] In saline-alkali farmland, inhibited crop transpiration and water deficit are mainly caused by low soil moisture content and high soil salinity. During the peak growth period of crops, especially the critical growth period, sufficient irrigation is generally implemented to meet the crop's water requirements in order to obtain a considerable yield. Therefore, water stress can be ignored. However, salinity stress still exists, and the degree of crop water deficit (PWDI) at this time is mainly determined by soil salinity.

[0069]

[0070] In the formula: T a With T p The actual and potential transpiration rates of the crop are respectively (mm / s). -1 ); Average soil salinity in the root zone (g kg) -1 ); Soil salinity stress correction factor:

[0071]

[0072] In the formula: τ is the fitting parameter; Critical salt tolerance value for crops (g kg) -1 When the soil salinity exceeds this critical value, crops will suffer from salt stress. Soil salinity (g / kg) at which crops suffer permanent wilting due to severe stress -1 ).

[0073] During the peak growth period of crops, especially the critical growth period, the vegetation cover is high, soil surface evaporation is weak and accounts for a very small proportion of total evapotranspiration, and the influence of soil water and salt conditions on soil surface evaporation and crop transpiration is very similar. Therefore, it can be assumed that the relative evapotranspiration rate (Λ) is equal to the relative transpiration rate at this time. Then, from equations (1) and (2), we can know that:

[0074]

[0075] Obviously, given the relative evapotranspiration rate Λ, the salinity of the root zone soil can be estimated according to equation (3).

[0076] To date, remote sensing inversion technology for crop evapotranspiration based on the principle of surface energy balance has made significant progress. [1] :

[0077]

[0078] In the formula: λ is the latent heat of vaporization of water, which is generally taken as 2.47 × 10⁻⁶. 6 J kg -1 ;ρ wThe density of water is generally taken as 1.00 × 10⁻⁶. 3 kg m -3 ;R n Net radiative flux (Jm) -2 s -1 The calculation is performed using parameters retrieved from remote sensing images and meteorological data. The calculation method is detailed in the implementation process and steps below; G represents soil heat flux (Jm³). -2 s -1 ), with R n There is a certain correlation, which can be verified through R. n Seek; H dry Sensitive heat flux under extremely dry conditions (J m -2 s -1 ), H dry =R n -G;H wet The sensible heat flux (J / m³) under fully humidified conditions -2 s -1 ), can be calculated using saturated vapor pressure difference and external aerodynamic impedance, etc.; ET wet Evapotranspiration intensity (mm / s) under fully humidified conditions -1 ), ET wet =(H dy -H wet ) / (λρ w ).

[0079] The relative evapotranspiration rate (Λ) exhibits high spatiotemporal variability at the regional scale due to the combined influence of factors such as soil, crops, climate, and human activities. Therefore, its remote sensing inversion requires high spatiotemporal resolution. Typically, high spatial resolution remote sensing images have a large time span (i.e., low temporal resolution), while remote sensing images with a shorter time span often have low spatial resolution. To balance this contradiction, remote sensing fusion technology has emerged, which combines and complements remote sensing data with different spatiotemporal resolutions to simultaneously meet the high spatiotemporal resolution requirement. Currently, adaptive remote sensing image spatiotemporal fusion methods are widely used because they are applicable to surface parameters (such as Λ) that change rapidly over time. This invention uses the relative evapotranspiration rate (Λ) obtained by inverting two types of remote sensing images—satellite L (high spatial resolution, low temporal resolution) and satellite M (high temporal resolution, low spatial resolution)—as an example. The fusion model equation can be summarized as follows: [2] :

[0080] Λ(x w / 2 y w / 2 , t p ) = T m ×Λ m (x w / 2 y w / 2 , t p)+T n ×Λ n (x w / 2 y w / 2 , t p )

[0081]

[0082] In the formula: Λ(x w / 2 y w / 2 , t p ) represents the final t obtained after fusion. p A Λ value with high spatiotemporal resolution at a given time (a day within the study period excluding the first and last two days); Λ m (x w / 2 y w / 2 , t p ) is based on t m Time (Day 1 of the study period) Satellite L and M remote sensing images and t p The time-based prediction of t by satellite M imagery p The value of Λ at time; Λ n (x w / 2 y w / 2 , t p ) is based on t n Satellite L and M imagery at the time of study (last day of the study period) and t p The time-based prediction of t by satellite M imagery p The value of Λ at time T; m and T n t m and t n Time-based remote sensing imagery relative to t p Time weighting factor at any given moment; Λ L and Λ M These represent the Λ values ​​calculated based on satellite L and M remote sensing images, respectively; N is the number of similar pixels; V i The conversion factor is Λ. L and Λ M The slope of the linear regression of the value at the i-th similar pixel in the image; W i V represents the weight of the i-th similar pixel. i W i T m and T n For details on the calculation method, please refer to the implementation process and steps below.

[0083] In summary, by combining the surface energy balance model and the adaptive remote sensing image spatiotemporal fusion model, this invention proposes a fusion method to improve the spatiotemporal resolution of Λ, and then applies Equation (3) to invert and obtain the distribution pattern of crop root soil salinity at the regional scale.

[0084] In this specific implementation, one revisit cycle of satellite L is taken as the study period. The remote sensing images of satellites L and M on the first and last two days of the study period, as well as the remote sensing images of satellite M every day during the period, are used as the basic data. The implementation steps for remote sensing inversion of relative evapotranspiration rate Λ, data fusion, and inversion of the spatial distribution pattern of root zone soil salinity are as follows:

[0085] Step 1: Set the first and last days of the study period as reference times t. m and t n , get t m and t n Remote sensing images from satellite L at a given time and daily remote sensing images from satellite M during the study period were used to retrieve the required surface parameters from each image. [3] The main parameters include surface albedo α0, surface emissivity ε0, normalized difference vegetation index (NDVI), surface temperature T0 (K), and solar zenith angle SZA (°).

[0086] Step 2: Use the relevant parameters obtained in Step 1 to estimate the net radiation flux R in equation (4). n [1] :

[0087]

[0088] In the formula: R sd Downlink shortwave radiation (J m) -2 s -1 ), calculated using the solar zenith angle (SZA) and digital elevation model; R ld For downlink longwave radiation (J m -2 s -1 ), calculated using data such as temperature; σ is the Stefan-Boltzmann constant, typically taken as 5.67 × 10. -8 W m -2 K -4 ;

[0089] Step 3: Calculate vegetation cover f using the Normalized Difference Vegetation Index (NDVI) obtained in Step 1. c [4] :

[0090]

[0091] In the formula: NDVI soil The NDVI value for bare ground is typically taken as 0.05; NDVI veg This is the NDVI value when the vegetation is completely covered, and it is generally taken as 0.7.

[0092] Step 4, based on the surface temperature T0 and vegetation cover f obtained in steps 1 and 3 cIn addition to meteorological data such as atmospheric pressure, wind speed, and temperature, the frictional wind speed u is iteratively solved using atmospheric boundary layer similarity theory. * (ms -1 ), sensible heat flux H (J m -2 s -1 ) and Obukhov stability length (J m -2 s -1 ) [1] ;

[0093] Step 5: Based on the parameters obtained in steps 1-4 above, calculate the relative evapotranspiration rate Λ according to equation (4), and obtain t respectively. m and t n The reference time is based on Λ retrieved from satellite L and M remote sensing images. L and Λ M and the predicted time (t) p (The data obtained from MODIS remote sensing image inversion at times other than the reference time during the study period) M Spatial distribution patterns;

[0094] Step 6, using the above t m and t n Reference time is based on Λ obtained from high spatial resolution satellite L remote sensing imagery. L That is, Λ L (t m ) and Λ L (t n The distribution map selects similar pixels for each target pixel. For any target pixel, the intersection of similar pixels at two reference times is selected as the final similar pixel, and then the weight W of each similar pixel is calculated. i It mainly includes the following three items:

[0095] (1) Select similar pixels

[0096] Λ with high spatial resolution for each reference time L That is, Λ L (t m ) and Λ L (t n The distribution map uses a target pixel as the center and sets a search box of a certain size within its neighborhood. It then compares the difference between each pixel within the box and the target pixel. If the difference is less than a threshold, the pixel is considered a similar pixel. The threshold is determined by Λ. L The standard deviation of all cell values ​​in the distribution map and the estimated number of land cover types in the study area are used to determine [5] :

[0097] |Λ L (x i y i , tk )-Λ L (x w / 2 y w / 2 , t k )|≤σ×2 / e(k=m or n) (8)

[0098] In the formula: W is the width of the search window; (x w / 2 y w / 2 (x) is the center of the target pixel; i y i ) represents the position of the i-th similar pixel (i = 1, 2, ..., N, where N is the number of similar pixels, including the center pixel); σ is Λ L (t m ) or Λ L (t n The standard deviation of all cell values ​​in the matrix is ​​given by , and e is the estimated number of land cover types. Based on this, the intersection of the similar cells selected at two reference times is chosen as the final similar cell for the target cell.

[0099] (2) Calculate the weight W of each similar pixel. i

[0100] Weight W i The contribution of the i-th similar pixel to the target pixel is determined by the position of the similar pixel and the similarity between high-spatial-resolution and low-spatial-resolution pixels. The calculation process is as follows:

[0101] 1) Calculate the reference time Λ L Image is Λ L (t m ) and Λ L (t n The spatial distance d between the target pixel and the i-th similar pixel in the image) i :

[0102]

[0103] 2) Calculate the i-th similar pixel Λ at the reference time. L (Λ L (t m ) and Λ L (t n )) and Λ M (Λ M (t m ) and Λ M (t n The correlation coefficient R between images i (Note: Λ) M The image needs to be resampled to 30m resolution in order to match Λ L (Image matching):

[0104]

[0105] In the formula: Λ L (x i y i , t m ) and Λ L (x i y i , t n ) represent t respectively m and t n Time Λ L Image (i.e., Λ) L (t m ) and Λ L (t n The value of the i-th similar pixel in )); Λ M (x i y i , t m ) and Λ M (x i y i , t n ) represent t respectively m and t n Time Λ M Image (i.e., Λ) M (t m ) and Λ M (t n The value of the i-th similar pixel in the array; μ L and σ L Representing Λ L (x i y i , t m ) and Λ L (x i y i , t n The mean and standard deviation of μ. M and σ M Representing Λ M (x i y i , t m ) and Λ M (x i y i , t n The mean and standard deviation of ().

[0106] 3) Based on spatial distance d i Correlation coefficient R with image i Calculate the composite index D of the i-th similar pixel. i :

[0107] Di =(1-R) i )×d i (11)

[0108] In general, the larger the composite index D value, the smaller the contribution of similar pixels to the target pixel. Therefore, D is used... i The normalized reciprocal of W is used as the weight W of the i-th similar pixel. i :

[0109]

[0110] W i The range is from 0 to 1, and the sum of the weights of all similar pixels is 1;

[0111] (3) Calculate Λ L and Λ M The conversion coefficient V between the i-th similar pixel values ​​in the image i :

[0112]

[0113] Step 7, as shown in equation (5), combines the weights W of each similar pixel. i and conversion factor V i According to t m and t n Λ of time L (i.e. Λ) L (t m ) and Λ L (t n )) and Λ M (i.e. Λ) M (t m ) and Λ M (t n Images and t p Λ of time M (i.e. Λ) M (t p ) and Λ M (t p Imagery, prediction t p The high spatial resolution Λ at each time point yields the distribution pattern of Λ with high spatiotemporal resolution for each day within the study period. Overall, for Λ at two reference times... L (i.e. Λ) L (t m ) and Λ L (t n )), with t p The closer the time, the greater its impact on the target pixel; therefore, time weights need to be set for the two reference times separately. The time weight factor T in equation (5) m With T nResampled Λ M Image at reference time t m With t n and the predicted time t p The range of change between them determines:

[0114]

[0115] Step 8, as shown in equation (3), utilizes the average value of each pixel obtained by fusion within the research period, which has high spatiotemporal resolution Λ. Estimate the average soil salinity in the root zone You can obtain Spatial distribution patterns.

[0116] Example:

[0117] This invention specifically focuses on the Manas River basin in Xinjiang (approximately 11,000 km²). 2 The following is an example of remote sensing inversion of root zone soil salinity in cotton fields irrigated under drip film from 2000 to 2020. Satellites L and M were selected as Landsat (high spatial resolution, low temporal resolution: 16 days, 30 meters) and MODIS (high temporal resolution, low spatial resolution: 1 day, 1 km), respectively. Based on relevant literature and experimental research results, the required cotton salt tolerance threshold for equation (2) was obtained. and The values ​​were 2.99 and 10.48 g / kg, respectively. -1 The parameter τ, optimized from field trials, is 2.58. Based on Landsat and MODIS remote sensing image data of the Manas River Basin from 2000 to early July 2020 (early cotton boll formation stage), the search window width w in equation (8) is set to 5, the estimated number of land cover types e is set to 20, and then the soil salinity of the cotton root layer (0-80cm) during this period is inverted using the technology of this invention. Its spatiotemporal evolution is as follows: Figure 1 As shown (only 5 years of results are displayed). Figure 1 This indicates that over the past 21 years, with the continuous expansion of cotton fields, the increase in the number of years of cultivation, and the widespread application of drip irrigation technology under plastic film, the overall salinity of the root zone soil in cotton fields in the watershed has shown a downward trend, with an average annual decrease of approximately 0.09 g kg. -1 Prior to 2011, due to higher irrigation quotas, the salinity of the cotton field root zone soil decreased rapidly, with an average annual decrease of 0.18 g / kg. -1However, in the past five years (2015-2020), due to the annual reduction in irrigation quotas, the rate of decrease in root zone soil salinity has slowed down and even rebounded. At the beginning of the 21st century, lightly and moderately saline-alkali soils accounted for 59.8% and 39.9% of cotton fields in the watershed, respectively. Moderately saline-alkali soils were mostly distributed in the middle and lower reaches of the watershed. However, by 2020, most of the moderately saline-alkali soils had transformed into lightly saline-alkali soils or even non-saline-alkali soils, resulting in the proportion of lightly saline-alkali soils in the entire watershed rising to 69.7%, and the proportion of non-saline-alkali soils reaching as high as 28.2%. Meanwhile, the average salinity of the cotton field root zone soil decreased to 3.93 g / kg. -1 .

[0118] To verify the inversion results, in early July each year from 2018 to 2020 (early stage of flowering), in the Anjihai Irrigation Area of ​​the middle reaches of the Manas River Basin (… Figure 2 a, with a total area of ​​880 km² 2 Soil samples were collected in layers and their salinity was measured in the Mosuowan Irrigation Area, downstream of the Manas River Basin, in early July 2020, where soil salinization is relatively severe. Figure 2 a, with a total area of ​​1140 km² 2 Sampling and measurements were also conducted in areas with relatively mild soil salinization. In the absence of prior sample information, a regular grid method was used to establish the sampling points, with the number of sampling points estimated based on the proportion of cotton field area to the total irrigated area.

[0119]

[0120] In the formula: n is the number of sampling points; d is the allowable error, taken as 10%; z is the reliability index, taken as 1.64 at a 90% confidence level; p is the proportion of cotton planting area to the total irrigation area, which is 65% and 48% for Anjihai Irrigation District and Mosuowan Irrigation District, respectively. Using the Create Fishnet function of ArcGIS 1v.3, the grid was laid out according to the principle of equal spatial spacing, and the sampling point locations were determined. Calculations showed that the sampling grid size for both Anjihai Irrigation District and Mosuowan Irrigation District was 4km (Level I grid, with 57 and 54 sampling points respectively). Figure 2 (b and 2c). For the Anjihai Irrigation District, where soil salinization is relatively severe and requires special attention, in addition to the Level I grid, Level II (500m × 500m, 81 sampling points) was also conducted in severely salinized areas in 2018 and 2019. Figure 2 d) and Level III (100m × 100m, 36 sampling points, Figure 2e) Nested grid sampling. Based on the pre-determined coordinates of the sampling points, GPS was used for positioning. For each sampling point, soil samples were collected from the middle of the soil layers at 0–5, 5–10, 10–20, 20–40, 40–60, and 60–80 cm using a soil auger next to the drip irrigation tape. The soil salinity was then measured indoors to obtain the average soil salinity of the 0–80 cm root zone.

[0121] The comparison between the inverted values ​​and measured values ​​of root soil salinity at all sampling points in the Anjihai Irrigation District and the Mosuowan Irrigation District is as follows: Figure 3 As shown. Figure 3 This indicates that the technology of this invention can reliably obtain the salinity of the root zone soil in cotton fields. Therefore, the spatiotemporal evolution of salinity in the root zone soil of cotton fields under drip irrigation with plastic film in the Manas River Basin of Xinjiang from 2000 to 2020 is also reasonable and reliable, and can provide a scientific basis for adjusting the planting structure of farmland, allocating water and soil resources, and optimizing irrigation and fertilization systems in this basin.

[0122] The contents not described in detail in this specification are existing technologies known to those skilled in the art.

[0123] References:

[0124] [1] Su Z. The Surface Energy Balance System (SEBS) for estimation of turbulent heat fluxes. Hydrology and Earth System Sciences Discussions, 2002, 6(1): 85-99.

[0125] [2] Zhu X, Chen J, Gao F, Chen

[0126] [3]Huang C,Li Y,Gu J.Improving estimation of evapotranspiration underwater-limited conditions based on SEBS and MODIS data in arid regions.RemoteSensing,2015,7(12):16795-16814.

[0127] [4]Sobrino J,Jiménez J,Paolin L.Land surface temperature retrievalfrom Landsat TM 5.Remote Sensing of Environment,2004,90(4):434-440.

[0128] [5]Gao F,Masek J,Schwaller M.On the blending of the Landsat and MODISsurface reflectance,predicting daily Landsat surface reflectance.IEEETransactions on Geoscience and Remote Sensing,2006,44(8):2207-2218.

Claims

1. A method for farmland root layer soil salt content inversion based on remote sensing data fusion, characterized in that, Includes the following steps: Step 1: Establish a quantitative relationship between crop relative transpiration rate and root zone soil salinity based on field trial data: ； In the formula: and T p These represent the actual and potential transpiration rates of the crop, respectively. This represents the average soil salinity in the root zone. This is the soil salinity stress correction factor; φ L This is the critical value for crop salt tolerance. When the soil salt content is higher than this critical value, crops will suffer from salt stress. φ w This refers to the soil salinity corresponding to the permanent wilting of crops under severe stress. These are the fitting parameters; Step 2: Based on the approximate relationship between relative evapotranspiration rate and relative transpiration rate during the peak growth period of crops, the relative evapotranspiration rate is further established using equation (1). With the salinity of the root zone soil The relationship between the two, namely the remote sensing inversion model of root zone soil salinity: ； Step 3, determine the study period in the peak growth stage of crops, obtain satellite L and M remote sensing image data in the study area during the study period, and use the surface energy balance model to obtain the spatial distribution of based on each remote sensing image. Satellite L represents a satellite with high spatial resolution and low temporal resolution, and satellite M represents a satellite with high temporal resolution and low spatial resolution. Step 4: Referring to the adaptive remote sensing image spatiotemporal fusion model, the data obtained in Step 3 are processed... Perform fusion calculations to obtain the results within the research period. Daily high-resolution spatial distribution patterns; The adaptive remote sensing image spatiotemporal fusion model is as follows: ; In the above two equations: k = m or n ; w is the search box width;​ ( x w / 2 , y w / 2 ) is the center of the target pixel; The final result after fusion t p High spatiotemporal resolution at all times value, t p The time point refers to any day within the research period, excluding the first and last two days; For based on t m Time-based satellite L and M remote sensing images and t p Time-based satellite M-images jointly predicted t p Moment value, t m The timeframe is the first day of the research period; Based on t n Time-lapse satellite L and M imagery and t p Time-based satellite M-images jointly predicted t p Moment value, t n The timeframe is the last day of the research period; T m and T n They are respectively t m and t n Real-time remote sensing images relative to t p Time weighting factor at any given moment; and Calculated based on satellite L and M remote sensing images respectively value; N This represents the number of similar pixels. V i The conversion factor is the factor for... and The first in the image i The slope of linear regression when the values ​​at similar pixels are used; W i For the first i The weights of similar pixels; Step 5, based on the research period obtained in Step 4 High-resolution spatial distribution patterns were obtained daily to determine the study period. average The high-resolution spatial distribution pattern is obtained, and the result is obtained by inversion according to equation (2). Spatial distribution patterns; The specific steps of step 3 are as follows: Step 3-1: Set the first and last days of the study period as reference times. t m and t n , obtain t m and t n Remote sensing images from satellite L at specific times and daily remote sensing images from satellite M during the study period were used to retrieve the required surface parameters, including surface albedo, from each image. α 0. Surface emissivity ε 0. Normalized Difference Vegetation Index NDVI Surface temperature T 0. Solar zenith angle SZA ; Step 3-2 - Estimation of net radiation flux using the surface parameters obtained in Step 3-1 R n ; Step 3-3, using the product from Step 3-1 NDVI Calculating vegetation cover f c ; Step 3-4, based on the results obtained in steps 3-1 and 3-3 T 0 and f c In addition to atmospheric pressure, wind speed, and temperature, the frictional wind speed is iteratively solved. sensible heat flux H and Obukhov stability length L ; Step 3-5: Based on the parameters obtained in steps 3-1, 3-2, 3-3, and 3-4 above, the following parameters are obtained by inversion according to the surface energy balance model: t m and t n Reference time is based on remote sensing imagery from satellites L and M. Value, that is , and and a certain predicted time t p Based on satellite M remote sensing imagery Value, that is ; The above t p The times other than the reference time within the study period, i.e., times other than the reference time within the study period. t m and t n Any day other than; The surface energy balance model is as follows: ； In the above formula: The latent heat of vaporization of water; ρ w The density of water; R n Net radiative flux; G Soil heat flux; H dry = R n – G， The sensible heat flux under extremely dry conditions; H wet For sensible heat flux under fully humid conditions; ET wet For evapotranspiration intensity under fully humid conditions, ; The specific steps of step 4 are as follows: Step 4-1, using the reference time t m and t n Based on satellite L-type remote sensing imagery and The distribution map selects similar pixels for each target pixel within the study area. For any target pixel, the intersection of similar pixels at two reference times is selected as the final similar pixel, specifically: For each reference time with high resolution Each pixel is centered on the target pixel, and a search box of a certain size is set within its neighborhood. The difference between each pixel within the box and the target pixel is compared. If the difference is less than a threshold, the pixel is considered a similar pixel to the target pixel. Then, the search is performed. t m and t n The intersection of similar pixels of the target pixel at any given time is taken as the final similar pixel of the target pixel; The above threshold is determined by the following formula. or n ; In the above formula: ( x i , y i ) is the first i The location of similar pixels i = 1, 2, …, N , N It is the number of similar pixels, including the center pixel; yes The standard deviation of all cell values ​​in the dataset; This is an estimated number of land cover types; Step 4-2, Calculate the reference time Right now and The target pixel in the distribution map and the first i Spatial distance of similar pixels d i Then calculate the reference time. i Similar pixels and Correlation coefficient between R i ; ; and Represent t m and t n time The first in the image i The values ​​of similar pixels, t m and t n time Image and and Represent t m and t n time The first in the image i The values ​​of similar pixels, t m and t n time Image and and Represent and The mean and standard deviation; and Represent and The mean and standard deviation; According to the obtained d i And R i , the comprehensive index of the target pixel and the similar pixel is calculated i D i ;​ ； Will D i The normalized reciprocal as the first i Weights of similar pixels W i ; ； Step 4-3, calculation and the conversion factor between the i similar pixel values V i ; ； Step 4-4, combined W i and V i ,according to t m and t n Moment and Distribution map and t p Moment Distribution map, prediction t p High spatial resolution of time High spatiotemporal resolution was obtained for each day during the study period. The distribution pattern.

2. The farmland root layer soil salt content inversion method based on remote sensing data fusion according to claim 1, characterized in that: The study period is defined as one or two revisit cycles of a satellite L with high spatial resolution and low temporal resolution during the peak growth period of crops.

3. The application of the method for inverting farmland root zone soil salinity based on remote sensing data fusion as described in claim 1, characterized in that: The inversion method is applicable to both single-crop and multi-crop planting structures within the study area, but it is necessary to establish remote sensing inversion models of root soil salinity for each major crop, i.e., optimize and determine the soil salinity stress correction coefficient of equation (2). Fitting parameters in .

4. The application of the method for inverting farmland root zone soil salinity based on remote sensing data fusion as described in claim 1, characterized in that: Based on remote sensing image data of crops during their peak growth period over many years, the spatiotemporal evolution of root zone soil salinity over these years was obtained through inversion.

Images