A rice field recognition method based on time-series harmonic fitting and terrain temperature mask

Abstract

The present application belongs to the technical field of spaceborne radar data processing and remote sensing mapping, and particularly relates to a paddy field recognition method based on time series harmonic fitting and terrain temperature mask. The present application realizes noise suppression and periodic characteristic extraction of time series scattering through high-value interference removal and time series harmonic fitting, obtains a harmonic mask image, a terrain mask image and a temperature mask image through spatial mask analysis respectively, and completes paddy field recognition in combination with the three masks, and can be suitable for large-scale paddy field recognition. The method provided by the present application effectively removes the space-time noise of radar data, reduces the complexity of time series analysis, realizes large-area high-resolution paddy field recognition without obtaining prior information of rice phenology, can be applied to large-scale and long-time-span rice planting area monitoring, and provides strong information support for policy making of the agricultural department.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of spaceborne radar data processing and remote sensing mapping, and particularly relates to a paddy field recognition method based on time series harmonic fitting and terrain temperature mask. BACKGROUND

[0002] Spaceborne radar can image the earth at all times and in all weather, and plays an irreplaceable role in earth observation under cloudy and foggy weather conditions. The time series scattering data provided by spaceborne radar can be used to identify the growth process differences between rice and non-rice crops, and is therefore widely used in paddy field recognition tasks.

[0003] Existing methods are mainly divided into machine learning methods and phenology recognition methods, both of which need to obtain the rice phenology information of the study area in advance, i.e., the time range of rice transplanting period, tillering period, heading period and harvesting period, and obtain multi-temporal scattering data within these time ranges for feature analysis.

[0004] The machine learning method obtains scattering data at least in two different rice growth periods, constructs a large number of paddy field samples based on ground sampling, and realizes image classification and paddy field recognition by driving machine learning modeling with samples.

[0005] The phenology recognition method obtains scattering data of the rice transplanting period and the heading period according to the rice phenology information, detects pixels with significant scattering enhancement through scattering difference analysis of different periods, and realizes paddy field recognition.

[0006] The premise of achieving high precision by the above two methods is that the rice phenology in the study area has time consistency, so these methods are only suitable for small-scale paddy field recognition (such as city and county level). When the study area is expanded (such as provincial level), the above two methods are difficult to apply, because in different regions far apart in geography, there are differences in light, temperature, precipitation and social and economic conditions, resulting in time inconsistency of rice phenology. SUMMARY

[0007] In view of the above problems or deficiencies, in order to solve the problem that the existing paddy field recognition method is applicable to a small range, the application provides a paddy field recognition method based on time series harmonic fitting and terrain temperature mask. Noise suppression and periodic feature extraction of time series scattering are realized through high-value interference removal and time series harmonic fitting. Harmonic mask map, terrain mask map and temperature mask map are obtained through spatial mask analysis respectively, and paddy field recognition is completed by combining the three masks.

[0008] A paddy field recognition method based on time series harmonic fitting and terrain temperature mask, the specific steps are as follows:

[0009] Step 1. High-value interference removal

[0010] Acquire time-series scattering data of all VH (Vertically transmitting and Horizontally receiving) polarizations provided by the spaceborne radar throughout the year, and convert the data values ​​from logarithmic units (VH)... dB Convert to linear units (VH) linear ).

[0011]

[0012] Data with linear unit values ​​higher than 0.01 were considered high-value interference unrelated to rice growth and removed; the remaining data were considered valid scattering data.

[0013] Step 2. Timing Harmonic Fitting

[0014] For the target location, acquire all the effective scattering data for the whole year to form the temporal scattering S[t] of that point, where t represents time and the value range is normalized to 0 to 1, corresponding to the radar imaging time from the 1st day to the 365th day of the year.

[0015]

[0016] Establish the harmonic formula for the fourth-order component, where a is a first-order constant term, representing the annual time-series scattering mean, A i P represents the magnitude of the i-th cosine term. i This represents the phase of the i-th cosine term. (Through...) Least square regression with fourth-order components yields undetermined coefficients a and A. i P i The value of .

[0017] Step 3. Perform spatial mask analysis on the target location point, including harmonic mask F1, terrain mask F2 and temperature mask F3.

[0018] Harmonic mask diagram F1:

[0019] Calculate the sum of cosine amplitudes A of the harmonic fitting. s .

[0020]

[0021] Where A j A represents the magnitude of the j-th cosine term. s This represents the temporal scattering fluctuation level at that location. More than 500 paddy field sample points were obtained through ground surveys, and the harmonic parameters a and A of these sample points were extracted. s ; and the upper and lower limits U of a were obtained through numerical statistical analysis. a and L a and As The upper and lower limits of the value U A and L A Create a new blank image f1 with all pixel values ​​set to 0. For the target location point, when its harmonic parameters a and A s satisfy:

[0022] L a ≤a≤U a L A ≤A s ≤U A

[0023] The pixel value at position f1 is set to 1 (if the harmonic condition is met), otherwise it remains 0 (if the harmonic condition is not met). Thus, the harmonic mask image F1 is obtained.

[0024] Terrain mask image F2:

[0025] Acquire digital elevation data and calculate the terrain slope using the elevation gradient method. Create a new blank image f2 with all pixel values ​​set to 0. For the target location point, if the terrain slope is less than 5.5°, set the pixel value at the corresponding location in f2 to 1 (slope condition met); otherwise, keep it at 0 (slope condition not met). This yields the terrain mask image F2.

[0026] Temperature mask diagram F3:

[0027] Obtain the monthly average surface temperature data. For the target location point, obtain its 12-month data values ​​to form the time series temperature C[t] of that coordinate point, where t represents time and the value range is normalized to 0 to 1. The corresponding temperature acquisition time is from the 1st to the 12th month of the year.

[0028]

[0029] A harmonic formula for the second-order component is established, where b is a constant term representing the annual temperature average, B represents the amplitude of the cosine term, and Q represents the phase of the cosine term. The values ​​of the undetermined coefficients b, B, and Q are obtained through least squares regression. The time range where the temperature value on the C[t] fitted curve is higher than 18℃ is defined as the high-temperature duration at that location. A new blank image f3 is created, with all pixel values ​​set to 0. For the target location, if its high-temperature duration is greater than 60 days, the pixel value at the corresponding location in f3 is set to 1 (meeting the temperature condition); otherwise, it remains 0 (not meeting the temperature condition). Thus, the temperature mask image F3 is obtained.

[0030] Step 4. Paddy field identification and mapping

[0031] A new blank image f is created, and all pixel values are set to 0. For the target position, when the pixel values of the harmonic mask image F1, the terrain mask image F2 and the temperature mask image F3 are all 1, the pixel value of f at the corresponding position is set to 1 (representing the presence of rice), otherwise it remains 0 (representing the absence of rice). At this point, the rice field recognition result image F is obtained.

[0032] Principles involved in steps 1 to 4:

[0033] High-value interference removal principle:

[0034] The logarithmic unit time series scattering is conducive to visual observation, but the physical meaning of the data value is lost. The data value converted to linear unit represents the surface scattering intensity. VH polarization represents the working mode of the radar antenna, which transmits a vertically polarized signal and receives a horizontally polarized signal. The VH polarization data only appears a data value higher than 0.01 in the scenes of buildings, rocks and radio frequency interference. The data value higher than 0.01 in the rice field is an interference signal unrelated to rice growth. Removing high-value interference will significantly improve the accuracy of time series analysis.

[0035] Principle of time series harmonic fitting:

[0036] Rice can be planted for three rounds in a rice field in a year, and the scattering of each round will show an approximate cosine process of first decreasing and then rising. Therefore, a fourth-order component harmonic formula can be used to fit the annual rice field scattering. The first constant component represents the mean value of the annual rice field scattering, and the three cosine components represent the periodic scattering fluctuations repeated once, twice and three times a year, respectively. Through least squares fitting, the undetermined coefficients in the constant component and the cosine component can be obtained. The harmonic fitting curve of the annual rice field scattering removes random signal disturbances and retains the periodic characteristics of the annual rice field scattering.

[0037] Principle of spatial mask analysis:

[0038] The soil water content of the rice field is relatively high, so its annual time series scattering mean value is lower than that of non-rice objects. The soil moisture and plant coverage change significantly during the growth of rice, so the annual time series scattering fluctuation degree of rice is higher than that of non-rice objects. Obtain the rice field sample points, and the sum of the first constant component a and the cosine amplitude A s Numerical statistical analysis can determine the numerical upper and lower limits of a and A s The a and A s of the rice field should be located within the corresponding numerical upper and lower limit intervals.

[0039] Rice needs to be irrigated with standing water during its growth, so the terrain slope of the rice field should not be greater than 5.5°. The rice field needs to meet the condition that the slope is less than 5.5°.

[0040] The annual air temperature of the rice planting area is a unimodal pattern, and the temperature continuously rises from January to July, and continuously decreases from July to December. The rice generally takes 60 days from the jointing stage to the heading stage, and the air temperature of more than 18 DEG C is required to ensure the photosynthetic accumulation of the plant. Therefore, the rice field needs to meet the condition that the high temperature duration is higher than 60 days.

[0041] Rice field recognition mapping principle:

[0042] When the harmonic parameters, terrain slope and high temperature duration of the target position point meet the conditions of the existence of the rice field at the same time, the position is the rice field.

[0043] In summary, the present application realizes noise suppression and periodic characteristic extraction of time series scattering by high-value interference removal and time series harmonic fitting, obtains harmonic mask map, terrain mask map and temperature mask map respectively through spatial mask analysis, and completes rice field recognition in combination with the three kinds of masks. It can be suitable for large-scale rice field recognition. The existing machine learning method and the phenology recognition method require that the rice phenology of the study area has time consistency, and it is difficult to be popularized in a large area. BRIEF DESCRIPTION OF DRAWINGS

[0044] Figure 1 It is a flowchart of the present application;

[0045] Figure 2 It is a numerical statistical analysis graph of the embodiment. DETAILED DESCRIPTION

[0046] The present application will be further described in detail in combination with the drawings and embodiments.

[0047] A rice field recognition method based on time series harmonic fitting and terrain temperature mask (as shown in Figure 1 The development environment of the present embodiment is GEE (Google earth engine), and the programming language is JavaScript.

[0048] Step 1: Taking the Chinese mainland land south of Qinling and Huaihe River as a research area, searching all ascending track mode VH polarized scattering data of Sentinel-1 in 2020, using image.pow function to realize the conversion of pixel value from logarithmic unit to linear unit, searching the pixel points with pixel value higher than 0.01, and using image.updateMask function to realize the removal of high-value pixels.

[0049] Step 2: Using the imageCollection function to construct time series scattering, using the image.addBands function to establish the fourth-order harmonic component, using the image.linearRegression function to realize the least square fitting of time series scattering and harmonic component, and obtaining the undetermined coefficients a, A i and Pi .

[0050] Step 3, spatial mask analysis of target location points, including harmonic mask map F1, terrain mask map F2, and temperature mask map F3.

[0051] Harmonic mask map F1:

[0052] Calculate the sum of the three cosine amplitudes A s . Obtain 600 rice field sample points through ground investigation, extract the harmonic parameters a and A s . Display a and A in box plot s Numerical statistical distribution, determine the upper and lower limits of the value. Create a blank image f1, set all pixel values to 0. For the target location point, when its harmonic parameters a and A s value meet

[0053] 0.006≤a≤0.025, 0.008≤A s ≤0.029

[0054] Set the pixel value of f1 corresponding position to 1, otherwise keep it as 0; at this point, obtain the harmonic mask map F1.

[0055] Terrain mask map F2:

[0056] Retrieve the NASA DEM digital elevation data, and calculate the terrain slope using the Terrain.slope function. Create a blank image f2, set all pixel values to 0. For the target geographic location, when its terrain slope is less than 5.5°, set the pixel value of f2 corresponding position to 1 (satisfy the slope condition), otherwise keep it as 0 (not satisfy the slope condition). At this point, obtain the terrain mask map F2.

[0057] Temperature mask map F3:

[0058] Retrieve the 2020 ERA5-Land surface air temperature monthly average data, for the target location point, obtain its 12-month pixel value, and construct the time series temperature of the point. Use the image.addBands function to establish the second-order harmonic component, and use the linearRegression function to realize the least squares fitting of the time series scatter and the harmonic component, and obtain the undetermined coefficients b, B and Q. Search for the time range where the temperature curve is higher than 18℃, define it as the high temperature duration. Create a blank image f3, set all pixel values to 0. For the target location point, when its high temperature duration is greater than 60 days, set the pixel value of f3 corresponding position to 1 (satisfy the temperature condition), otherwise keep it as 0 (not satisfy the temperature condition). At this point, obtain the temperature mask map F3.

[0059] Step 4, rice field recognition mapping: create a blank image f, set all pixel values to 0.

[0060] For the target position point, when the pixel values of the harmonic mask figure F1, the terrain mask figure F2 and the temperature mask figure F3 are all 1, the pixel value of f corresponding position is set to 1 (representing the existence of rice), otherwise it remains 0 (representing the non-existence of rice). At this point, the rice field recognition result figure F is obtained.

[0061] In the embodiment, Sentinel-1 data, NASA SEM data and ERA5-Land data in 2020 are processed, Figure 2 The numerical statistical analysis chart of the embodiment is shown in the following table.

[0062] It can be seen from the above embodiment that the rice field extraction of the Chinese land territory south of the Qinling-Huaihe River is realized, and the precision test result reaches 80%. The method provided by the application effectively removes the space-time noise of the radar data, reduces the complexity of the time series analysis, and can realize the large-area high-resolution rice field recognition without obtaining the rice phenology prior information, and can be applied to the rice planting area monitoring in a larger range and a longer time span, thereby providing strong information support for the policy making of the agricultural department.

Claims

1. A method for identifying paddy fields based on time-series harmonic fitting and topographic temperature masking, characterized in that, The specific steps are as follows: Step 1. High-value interference removal; Acquire time-series scattering data of all VH polarizations throughout the year from the spaceborne radar, and convert the data values ​​from logarithmic units VH. dB Convert to linear unit VH linear ; ； Data with a linear unit value higher than 0.01 are removed, and the remaining data are considered valid scattering data; Step 2. Timing harmonic fitting; For the target location, all valid scattering data for the whole year are obtained to form the temporal scattering S[t] of the point, where t represents time and the value range is normalized to 0 to 1, corresponding to the radar imaging time from the 1st day to the 365th day of the year; ； Establish the harmonic formula for the fourth-order component, where a is a first-order constant term, representing the annual time-series scattering mean, A i P represents the magnitude of the i-th cosine term. i Represents the phase of the i-th cosine term; through Least square regression with fourth-order components yields undetermined coefficients a and A. i P i The value; Step 3. Perform spatial mask analysis on the target location points, including harmonic mask F1, terrain mask F2, and temperature mask F3; Harmonic mask diagram F1: Calculate the sum of cosine amplitudes A of the harmonic fitting. s ; ； Where A j A represents the magnitude of the j-th cosine term. s This represents the temporal scattering fluctuation level at that location; over 500 paddy field sample points were obtained through ground surveys, and the harmonic parameters a and A of these paddy field sample points were extracted. s ; and the upper and lower limits U of a were obtained through numerical statistical analysis. a and L a and A s The upper and lower limits of the value U A and L A Create a new blank image f1 with all pixel values ​​set to 0; for the target location point, when its harmonic parameters a and A s satisfy: L a ≤a≤U a ,L A ≤A s ≤U A； The pixel value at position f1 is set to 1, and if the harmonic condition is not met, it remains 0. Thus, the harmonic mask image F1 is obtained. Terrain mask image F2: Acquire digital elevation data and calculate the terrain slope using the elevation gradient method; create a new blank image f2 with all pixel values ​​set to 0; for the target location point, if its terrain slope is less than 5.5°, set the pixel value of the corresponding position in f2 to 1, and keep it at 0 if the slope condition is not met. At this point, the terrain mask image F2 is obtained. Temperature mask diagram F3: Obtain monthly average surface temperature data. For a target location point, obtain its 12-month data values ​​to form the time series temperature C[t] of that coordinate point, where t represents time and the value range is normalized to 0 to 1. The corresponding temperature acquisition time is from the 1st to the 12th month of the year. ； A harmonic formula for the second-order component is established, where b is a constant term representing the annual average temperature, B represents the amplitude of the cosine term, and Q represents the phase of the cosine term. The values ​​of the undetermined coefficients b, B, and Q are obtained through least squares regression. The time range of temperature values ​​above 18℃ on the C[t] fitting curve is defined as the high-temperature duration at that location. A new blank image f3 is created, with all pixel values ​​set to 0. For the target location, when its high-temperature duration is greater than 60 days, the pixel value at the corresponding location of f3 is set to 1; otherwise, it remains 0. Thus, the temperature mask image F3 is obtained. Step 4. Paddy field identification and mapping; Create a new blank image f, and set all pixel values ​​to 0. For the target location point, when the pixel values ​​of the harmonic mask F1, the terrain mask F2, and the temperature mask F3 are all 1, the pixel value at the corresponding position of f represents the presence of rice and is set to 1; otherwise, it represents the absence of rice and is kept at 0. At this point, the rice field identification result image F is obtained.

2. The paddy field identification method based on time-series harmonic fitting and topographic temperature mask as described in claim 1, characterized in that: In step 1, the time-series scattering data of all VH polarizations provided by the spaceborne radar throughout the year are Sentinel-1 satellite data; the image.updateMask function is used to remove data with linear unit values ​​higher than 0.

01.

3. The paddy field identification method based on time-series harmonic fitting and topographic temperature mask as described in claim 1, characterized in that: In step 3, the digital elevation data used is NASADEM data, and the terrain slope is calculated using the Terrain.slope function.

4. The paddy field identification method based on time-series harmonic fitting and topographic temperature mask as described in claim 1, characterized in that: In step 3, the monthly average surface air temperature data is obtained from the ERA5-Land dataset, and the second harmonic components are established using the image.addBands function.

Images