Dynamic determination method of soil moisture threshold for satellite remote sensing irrigation area interpretation

Abstract

The present invention discloses a soil moisture threshold dynamic determination method for satellite remote sensing interpretation of actual irrigation area, collects time-continuous ground measured soil moisture data, constructs an equation to describe the surface soil moisture attenuation law, and determines the soil moisture jump amplitude under different time intervals after irrigation; uses the measured soil moisture data to evaluate the error of the inversion result, and uses the error transmission principle to optimize the multi-source data combination scheme; combines the soil moisture attenuation equation, the irrigation area identification formula, and the data combination scheme error to construct an irrigation area identification dynamic threshold determination equation considering error transmission, and determines the irrigation area identification threshold under different time intervals and different data source combination schemes; determines the irrigation area identification threshold of soil moisture data at adjacent times within a certain time range, and calculates the irrigation area of ​​multiple rounds within the time range. The present invention can efficiently and accurately obtain the irrigation area identification threshold and support the identification of the irrigation area.

Description

Technical Field

[0001] The invention belongs to the field of remote sensing interpretation of agricultural irrigation area, and in particular relates to a soil moisture threshold dynamic determination method for satellite remote sensing actual irrigation area interpretation. Background Art

[0002] my country's agricultural irrigation water projects are numerous, widespread and large in volume, accounting for more than 60% of the country's total water consumption. Comprehensive coverage of agricultural water metering is difficult and costly, and it is difficult to fully and accurately grasp the situation of agricultural water use, which has become a prominent weak link in water use supervision. In particular, according to the results of the second and third national land and resources surveys, the national land use structure and scale have changed significantly, and the irrigation area of ​​farmland in some areas has changed significantly. Dynamically grasping the changes in the actual irrigation area of ​​agriculture is very necessary to improve the decision-making ability of agricultural water resources allocation and management, strengthen agricultural water management, ensure water for food production, and support high-quality economic and social development.

[0003] Compared with the irrigation area identification technology based on the jump characteristics of soil moisture after irrigation, the irrigation area identification technology based on the difference in crop growth and evapotranspiration is more timely and has the ability to distinguish irrigation times. Soil moisture data with low spatial resolution usually has good temporal continuity and can achieve daily or quasi-daily regional coverage, but it is difficult to reflect the irrigation situation at the scale of plots and cannot be used for accurate monitoring of agricultural irrigation areas. Soil moisture data with high spatial resolution can finely depict the irrigation situation of different plots, but it is difficult to achieve temporal continuous monitoring due to the influence of weather and revisit cycles. At present, the multi-source satellite coordination method has significantly improved the temporal continuity of high spatial resolution data, but due to the unstable time interval and inconsistent inversion errors of multi-source data, it is difficult to determine the threshold when using soil moisture data to identify the irrigation area. Therefore, a threshold dynamic determination method for interpreting the irrigation area of ​​optical remote sensing soil moisture is invented, which is of great significance for the business identification of multi-round irrigation information at the scale of plots and improving the level of agricultural water management. Summary of the invention

[0004] In order to solve the technical problems existing in the background technology, the present invention aims to provide a method for dynamically determining soil moisture thresholds for interpreting actual irrigation areas using satellite remote sensing. First, based on the continuous soil moisture changes monitored by the ground soil moisture station, an equation is constructed to describe the soil moisture attenuation law, and then the soil moisture jump amplitude under different time intervals after irrigation is determined; secondly, the multi-source remote sensing data is inverted, and the inversion results are evaluated for errors using the measured soil moisture data, and the coordinated errors under different data source combination schemes are evaluated based on the error transfer formula to achieve the optimization of the multi-source data combination scheme; finally, the soil moisture attenuation equation, the irrigation area identification formula, and the data combination scheme error are combined to construct an irrigation area dynamic identification threshold determination equation considering error transfer, and determine the irrigation area identification threshold under different time intervals and different data source combination schemes. Through this method, the irrigation area identification threshold can be obtained more efficiently and accurately, supporting high-precision and high-frequency identification of the irrigation area, thereby improving the level of refined management of agricultural water use.

[0005] In order to solve the technical problem, the technical solution of the present invention is:

[0006] A soil moisture threshold dynamic determination method for satellite remote sensing interpretation of actual irrigation area, comprising the following steps:

[0007] S1: Collect the continuous soil moisture measured on the ground, construct an equation to describe the attenuation law of surface soil moisture, and determine the soil moisture jump amplitude under different time intervals after irrigation;

[0008] S2: Invert soil moisture based on multi-source remote sensing data, use measured soil moisture data to evaluate the error of the inversion result, evaluate the coordinated error under different data source combination schemes based on the error transfer formula, and optimize the multi-source data combination scheme;

[0009] S3: Combine the soil moisture attenuation equation, irrigation area identification formula, and data combination scheme error to construct an irrigation area identification dynamic threshold determination equation that takes error transmission into account, and dynamically determine the irrigation area identification threshold under different time intervals and different data source combination schemes;

[0010] S4: Repeat steps S1-S3 to determine the irrigation area recognition threshold within a certain time range, and calculate multiple rounds of irrigation areas within the time range.

[0011] Further, the step S1 specifically includes:

[0012] S11: Collect time-continuous soil moisture and precipitation data, and determine the time range of soil moisture changes caused by irrigation after deducting the effect of precipitation;

[0013] S12: Analyze the variation of surface soil moisture (10 cm) caused by irrigation in multiple time ranges, and obtain the equation for soil moisture decay over time;

[0014] SM t+D =(SM t , D)

[0015] Among them, SM t+D represents the surface soil moisture after the Dth irrigation day, SM t It indicates the surface soil moisture on the day after irrigation;

[0016] S13: Comprehensively consider the soil moisture levels before and after multiple irrigation events and determine the soil moisture reference value SM on the irrigation day. t Compared with the soil moisture reference value SM on the day before irrigation t-1 .

[0017] Further, the step S2 specifically includes:

[0018] S21: Collect multi-source satellite remote sensing images and perform data preprocessing, select appropriate methods, and invert soil moisture

[0019] S22: Use the measured soil moisture data to evaluate the inversion result error. The main indicators include the correlation coefficient R 2 Error standard deviation σ;

[0020]

[0021]

[0022] In the formula represents the inverted soil moisture of the i-th verification point, SM represents the measured soil moisture, δ represents the error between the measured soil moisture and the inverted soil moisture, and n represents the number of verification points;

[0023] S23: The core of irrigation area identification based on soil moisture change is the difference in soil moisture before irrigation (t) and after irrigation (t+k). Combined with remote sensing data collection, the soil moisture difference at adjacent times is calculated by express;

[0024]

[0025] Combination Calculation equation and error propagation equation, The uncertainty of the calculation should be obtained from the uncertainty of soil moisture in the two periods. For the same data source, the uncertainty of the difference in soil moisture is The expression is as follows:

[0026]

[0027] The uncertainty in soil moisture differences for different data sources is:

[0028]

[0029] In the formula, S1 and S2 represent different data types respectively.

[0030] S24: Based on The calculation formula is used to calculate the uncertainty of soil moisture differences under different combinations of collected data sources. The smallest one is the optimal data source combination solution;

[0031] Further, the step S3 specifically includes:

[0032] S31: Identification of irrigation area based on soil moisture is to determine the difference in soil moisture changes in remote sensing inversion Whether the threshold T requirement is met, if so, irrigation is considered to have occurred;

[0033]

[0034] S32: The location and time of irrigation plots in the region are scattered, and the threshold T should also be determined under the most unfavorable conditions. For example, if the time interval between two phases of remote sensing data is 5 days, the most unfavorable condition is that irrigation occurs on the second day, and soil moisture is monitored by remote sensing after decaying for 3 days. The amplitude of soil moisture jump caused by irrigation is reduced, and the difficulty of irrigation identification increases;

[0035] T=SM t+D -SM t-1

[0036] Where SM t-1 represents the soil moisture before irrigation, SM t+D Soil moisture after irrigation D days, where D is determined by the time interval between two periods of soil moisture.

[0037] S33: The soil moisture difference considering the inversion uncertainty is In the comparison with the threshold T, the introduction of uncertainty may result in a reduction in the difference Therefore, irrigation identification needs to consider the negative impact of uncertainty, and at the same time, combined with the attenuation law equation of measured soil moisture, the irrigation area identification equation considering the error characteristics of different time intervals and different data source combinations is obtained;

[0038]

[0039] Further, the step S4 specifically includes the following steps:

[0040] S41: repeating steps S1-S3, calculating soil moisture for optical remote sensing data collected over a period of time, and identifying irrigation areas for soil moisture data from adjacent periods of time;

[0041] S42: Obtaining the number of irrigation times for each pixel based on the accumulated irrigation areas identified in each round;

[0042]

[0043] Among them, IA y is the number of irrigations of the pixel in period y, m is the number of irrigation identification periods; IA x is the irrigation identification result of the pixel in the xth period.

[0044] Compared with the prior art, the advantages of the present invention are:

[0045] 1. The present invention proposes a method for dynamically determining soil moisture thresholds for satellite remote sensing interpretation of actual irrigation areas, which can dynamically determine irrigation identification thresholds under different time intervals and different data source combinations, providing important support for the commercial application of irrigation area identification.

[0046] 2. The present invention uses ground measured data as empirical reference and verification truth value, constructs surface soil moisture attenuation law equation, evaluates the error characteristics of soil moisture inversion, derives the error synthesis equation of multi-source soil moisture combination application, quantitatively evaluates the multi-source data collaboration theory error, and solves the problem of data source optimization under multi-source collaboration conditions.

[0047] 3. The threshold determination method proposed in the present invention can support high-frequency irrigation area identification, realize regional irrigation area identification in different rounds, and support more detailed irrigation process characterization and irrigation water volume estimation. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 The method flow of the present invention;

[0049] Figure 2 An example of a time-continuous soil moisture dataset for an embodiment;

[0050] Figure 3-1 The error evaluation result of the Sentinel-2 data source of the embodiment;

[0051] Figure 3-2 It is the Landsat-8 data source error evaluation result of the embodiment;

[0052] Figure 3-3 The GF-6 data source error evaluation result of the embodiment;

[0053] Figure 3-4 The GF-1 data source error evaluation result of the embodiment;

[0054] Figure 4 An example of interpretation of the irrigation area of ​​optical remote sensing soil moisture in an embodiment;

[0055] Figure 5 This is the irrigation area recognition result of the time period of the embodiment. DETAILED DESCRIPTION

[0056] The specific implementation mode of the present invention is described below in conjunction with embodiments:

[0057] It should be noted that the structures, proportions, sizes, etc. shown in this specification are only used to match the contents disclosed in the specification so that people familiar with this technology can understand and read them, and are not used to limit the conditions under which the present invention can be implemented. Any structural modification, change in proportional relationship or adjustment of size should still fall within the scope of the technical content disclosed in the present invention without affecting the effects and purposes that can be achieved by the present invention.

[0058] At the same time, the terms such as "upper", "lower", "left", "right", "middle" and "one" cited in this specification are only for the convenience of description and are not used to limit the scope of implementation of the present invention. Changes or adjustments to their relative relationships should be regarded as the scope of implementation of the present invention without substantially changing the technical content.

[0059] like Figure 1 As shown, a soil moisture threshold dynamic determination method for satellite remote sensing interpretation of actual irrigation area is specifically applied in the southern part of Hebei Province, including the following steps:

[0060] Step 1: Collect the continuous soil moisture measured on the ground in the embodiment area, construct an equation to describe the attenuation law of surface soil moisture, and determine the soil moisture jump amplitude under different time intervals after irrigation;

[0061] Step 2: Invert soil moisture based on multi-source remote sensing data, use measured soil moisture data to evaluate the error of the inversion result, evaluate the coordinated error under different data source combination schemes based on the error transfer formula, and select the multi-source data combination scheme;

[0062] Step 3: Combine the soil moisture attenuation equation, irrigation area identification formula, and data combination scheme error to construct an irrigation area identification dynamic threshold determination equation that takes error transmission into account, and dynamically determine the irrigation area identification threshold under different time intervals and different data source combination schemes;

[0063] Step 4: Repeat steps S1-S3 to determine the irrigation area recognition threshold within a certain time range, and calculate the irrigation area of ​​multiple rounds within the time range.

[0064] Further, the step S1 specifically includes:

[0065] Step 11: Collect the time-continuous soil moisture content and precipitation data of the embodiment area, and after deducting the influence of precipitation, determine the time range of soil moisture change caused by irrigation, such as Figure 2 As shown;

[0066] Step 12: Analyze the variation of surface soil moisture (10 cm) caused by irrigation in multiple time ranges, and obtain the equation for soil moisture decay over time;

[0067]

[0068] Among them, SM t+D represents the surface soil moisture after the Dth irrigation day, SM t It indicates the surface soil moisture on the day after irrigation;

[0069] Step 13: Comprehensively consider the soil moisture levels before and after multiple irrigation events and determine the soil moisture reference value SM on the irrigation day t Compared with the soil moisture reference value SM on the day before irrigation t-1 In this embodiment, the soil moisture reference value SM caused by irrigation t 32% of the soil moisture reference value SM the day before irrigation t-1 is 12%.

[0070] Further, the step S2 specifically includes:

[0071] Step 21: Collect multi-source remote sensing data such as Gaofen-1 (GF-1), Gaofen-6 (GF-6), Landsat-8, and Sentinel-2, perform radiation calibration, atmospheric correction, geometric correction, and relative radiation correction, select the modified vertical drought index (MPDI) method, combine MPDI with measured soil moisture data, and invert soil moisture

[0072]

[0073] Where R red,v , R nir,v is the reflectivity of vegetation in the red and near-infrared bands, Fv is the vegetation coverage, and M is the slope of the soil line. For the selection of the slope of the soil line, it is necessary to fully consider the triangular feature space of the red and near-infrared bands in different periods, and construct an ideal soil line in the area where both extremely dry and extremely wet points exist. In this embodiment, the ideal soil line is mainly constructed based on the image data under low vegetation coverage conditions from January to March.

[0074] Step 22: Select the measured soil moisture data that are not included in the inversion for the inversion result error evaluation. The main indicators include the correlation coefficient R 2In this embodiment, a total of 80 manual sampling monitoring samples are used for soil moisture accuracy verification, and the error evaluation results of GF-1, GF-6, Landsat-8, and Sentinel-2 data sources are as follows: Figure 3-1 to Figure 3-4 As shown, the R 2 is 0.53, σ is 5.77, and R 2 is 0.65, σ is 4.95, and the R 2 is 0.60, σ is 5.19, and R of Sentinel-2 2 is 0.68, σ is 4.65;

[0075]

[0076] represents the inverted soil moisture of the i-th verification point, SM represents the measured soil moisture, δ represents the error between the measured soil moisture and the inverted soil moisture, and n represents the number of verification points;

[0077] Step 23: The core of the irrigation area identification based on soil moisture change is the difference in soil moisture before irrigation (t) and after irrigation (t+k). Combined with the remote sensing data collection, the soil moisture difference at adjacent times is calculated by express;

[0078]

[0079] Combination Calculation equation and error propagation equation, The uncertainty of the calculation should be obtained from the uncertainty of soil moisture in the two periods. For the same data source, the uncertainty of the difference in soil moisture is The expression is as follows:

[0080]

[0081] The composite uncertainty results for GF-1 when used independently to monitor irrigation area are 8.16, GF-6 is 7.00, Landsat-8 is 7.34, and Sentinel-2 is 6.58;

[0082] The uncertainty in soil moisture differences for different data sources is:

[0083]

[0084] In the formula, S1 and S2 represent different data types. Based on the above formula, the uncertainty synthesis results of GF-1 and GF-6 are 7.60, GF-1 and Landsat-8 are 7.76, GF-1 and Sentinel-2 are 7.41, GF-6 and Landsat-8 are 7.17, GF-6 and Sentinel-2 are 6.79, and Landsat-8 and Sentinel-2 are 6.97;

[0085] Step 24: Based on The calculation formula is used to calculate the uncertainty of soil moisture differences under different combinations of collected data sources. The smallest one is the optimal data source combination solution;

[0086] Further, the step S3 specifically includes:

[0087] Step 31: Identification of irrigation area based on soil moisture is to determine the difference in soil moisture changes in remote sensing inversion Whether the threshold T requirement is met, if so, irrigation is considered to have occurred;

[0088]

[0089] Step 32: The location and time of irrigation plots in the region are scattered, and the threshold T should also be determined under the most unfavorable conditions. For example, if the time interval between two phases of remote sensing data is 5 days, the most unfavorable condition is that irrigation occurs on the second day, and soil moisture is monitored by remote sensing after decaying for 3 days. The amplitude of soil moisture jump caused by irrigation is reduced, and the difficulty of irrigation identification increases;

[0090] T=SM t+D -SM t-1

[0091] Where SM t-1 represents the soil moisture before irrigation, SM t+D Soil moisture after irrigation D days, where D is determined by the time interval between two periods of soil moisture.

[0092] Step 33: The soil moisture difference considering the inversion uncertainty is In the comparison with the threshold T, the introduction of uncertainty may result in a reduction in the difference Therefore, irrigation identification needs to consider the negative impact of uncertainty. At the same time, combined with the attenuation law equation of measured soil moisture, the irrigation area identification equation considering the error characteristics of different time intervals and different data source combinations is obtained. The determination of the irrigation identification threshold within a period of time is shown in Table 1;

[0093] Same data source combination:

[0094]

[0095] Combination of different data sources:

[0096]

[0097] Table 1 Results of dynamic determination of irrigation area threshold for soil moisture change identification

[0098]

[0099]

[0100] Further, the step S4 specifically includes the following steps:

[0101] Step 41: Repeat steps S1-S3 to calculate the soil moisture of optical remote sensing data collected over a period of time, and identify the irrigation area of ​​soil moisture data at adjacent times; Figure 4 and Figure 5 shown.

[0102] Step 42: The irrigation times of each pixel are obtained by summing up the irrigation areas identified in each round;

[0103]

[0104] Among them, IA y is the number of irrigations of the pixel in period y, m is the number of irrigation identification periods, IA x is the irrigation identification result of the pixel in period x. The multi-frequency irrigation area identification results were verified based on field sampling samples, and the average accuracy reached 92%. The specific verification results are shown in Table 2.

[0105] Table 2 Accuracy verification results

[0106] date Verification Points Correct identification Accuracy 2020 / 4 / 12 2020 / 4 / 13 40 37 93％ 2020 / 4 / 16 40 35 88％ 2020 / 4 / 17 40 38 95％ 2020 / 4 / 20 40 36 90％ 2020 / 4 / 22 40 39 98％ 2020 / 4 / 23 40 39 98％ 2020 / 4 / 25 40 37 93％ 2020 / 4 / 27 40 38 95％ 2020 / 4 / 29 40 34 85％ 2020 / 4 / 30 40 36 90％

[0107] The preferred embodiments of the present invention are described in detail above, but the present invention is not limited to the above embodiments, and various changes can be made within the knowledge scope of ordinary technicians in this field without departing from the purpose of the present invention.

[0108] Many other changes and modifications may be made without departing from the concept and scope of the present invention.It should be understood that the present invention is not limited to the specific embodiments, and the scope of the present invention is defined by the appended claims.

Claims

1. A method for dynamically determining soil moisture threshold for satellite remote sensing irrigation area interpretation, characterized in that: The following steps are involved: S1: Collect the continuous soil moisture measured on the ground, construct an equation to describe the attenuation law of surface soil moisture, and determine the soil moisture jump amplitude under different time intervals after irrigation; Step S1 specifically includes the following sub-steps: S11: Collect time-continuous soil moisture and precipitation data, and determine the time range of soil moisture changes caused by irrigation after deducting the effect of precipitation; S12: Analyze the variation of surface soil moisture caused by irrigation in multiple time ranges and obtain the equation for soil moisture decay over time; SM t+D =f(SM t ,D) Among them, SM t+D represents the surface soil moisture after irrigation for D days, SM t It indicates the surface soil moisture on the day after irrigation; S13: Comprehensively consider the soil moisture levels before and after multiple irrigation events and determine the soil moisture reference value SM on the irrigation day. t Compared with the soil moisture reference value SM on the day before irrigation t-1 ; S2: Invert soil moisture based on multi-source remote sensing data, use measured soil moisture data to evaluate the error of the inversion results, evaluate the coordinated errors under different data source combination schemes based on the error transfer formula, and optimize the multi-source data combination scheme; Step S2 specifically includes the following sub-steps: S21: Collect multi-source remote sensing images and perform data preprocessing, and select appropriate methods to invert soil moisture S22: Use the measured soil moisture data to evaluate the inversion results, including the correlation coefficient R 2 , error standard deviation σ; In the formula represents the inverted soil moisture of the i-th verification point, SM represents the measured soil moisture, δ represents the error between the measured soil moisture and the inverted soil moisture, and n represents the number of verification points; S23: The core of irrigation area identification based on soil moisture change is the difference in soil moisture before irrigation (t) and after irrigation (t+k). Combined with remote sensing data collection, the soil moisture difference at adjacent times is calculated by express; Combination Calculation equation and error propagation equation, The uncertainty of the calculation should be synthesized from the uncertainty of soil moisture in the two periods. For the same data source, the uncertainty of soil moisture difference is The expression is as follows: The uncertainty in soil moisture differences for different data sources is: In the formula, S1 and S2 represent different data types respectively; S24: Based on The calculation formula is used to calculate the uncertainty of soil moisture differences under different combinations of collected data sources. The smallest one is the optimal data source combination solution; S3: Combine the soil moisture attenuation equation, irrigation area identification formula, and data combination scheme error to construct an irrigation area identification dynamic threshold determination equation that takes error transmission into account, and dynamically determine the irrigation area identification threshold under different time intervals and different data source combination schemes; Step S3 specifically includes the following sub-steps: S31: Identification of irrigation area based on soil moisture is to determine the difference in soil moisture changes in remote sensing inversion Whether the threshold T requirement is met, if so, irrigation is considered to have occurred; S32: The location and time of the irrigation plots in the region are scattered, and the threshold T is also determined under the most unfavorable conditions; T=SM t+D -SM t-1 Where SM t-1 represents the soil moisture before irrigation, SM t+D Soil moisture after irrigation D days, D is determined by the time interval between two periods of soil moisture; S33: The soil moisture difference considering the inversion uncertainty is In the comparison with the threshold T, the introduction of uncertainty leads to a reduction in the difference Combined with the attenuation law equation of measured soil moisture, the irrigation area identification equation considering the error characteristics of different time intervals and different data source combinations is obtained; S4: Repeat steps S1-S3 to determine the irrigation area recognition threshold of soil moisture data at adjacent times within a certain time range, and calculate multiple rounds of irrigation areas within the time range.

2. The method for dynamically determining soil moisture threshold for satellite remote sensing irrigation area interpretation according to claim 1, characterized in that: Step S4 specifically includes the following sub-steps: S41: repeating steps S1-S3, using multi-source satellite images to invert and calculate soil moisture for a period of time, and identifying the irrigation area of ​​soil moisture data at adjacent times; S42: Obtaining the number of irrigation times for each pixel based on the accumulated irrigation areas identified in each round; Among them, IA y is the number of irrigations of the pixel in period y, m is the number of irrigation identification periods; IA x is the irrigation identification result of the pixel in the xth period.

Images