A wheat phenology node information estimation method based on GNSS-IR technology

Abstract

The application discloses a wheat phenology node information estimation method based on GNSS-IR technology, processes satellite observation and orbit data obtained in advance to obtain standard format data, obtains signal-to-noise ratio multipath reflection signals and signal-to-noise ratio reflection components in a low elevation angle range in a system and frequency, carries out time-frequency conversion on the obtained reflection components, extracts frequency characteristics corresponding to frequency peak values, introduces an amplitude attenuation factor to construct a signal-to-noise ratio reflection signal model, fits the signal-to-noise ratio model by using a robust nonlinear least square algorithm to obtain a damping coefficient, extracts wheat phenology information, and uses a turning point of a long time sequence curve of the damping coefficient as a demarcation point of the wheat phenology. The application improves the time and space resolution of the wheat phenology node, and improves the precision and reliability of the estimated wheat phenology node.

Description

TECHNICAL FIELD

[0001] The application belongs to the cross technical field of satellite navigation positioning and remote sensing, and particularly relates to a wheat phenology node information estimation method based on GNSS-IR technology. BACKGROUND

[0002] With the rapid expansion of the application field of Global Navigation Satellite System (GNSS) technology, unlike the traditional GNSS navigation, positioning and timing services, the ground multi-path reflection signals which are often considered as noise interference in signal reception are found to be used for monitoring ground environmental parameters and other remote sensing researches. This technology is called GNSS-IR (Global Navigation Satellite System-Interferometry and Reflectometry). This technology depends on the existing GNSS and has the advantages of all-weather, continuity, low cost, rich signal resources, etc. It can overcome the shortcomings of traditional remote sensing methods such as information assimilation difficulty, low spatiotemporal resolution and high cost. At present, based on the existing observation stations and CORS (Continuously Operating Reference Stations) network, it has provided a wider, cheaper source of ground environmental information at the field scale for the scientific community.

[0003] Wheat is one of the important crops in the world, and timely and accurate measurement of its growth status is of great significance to food security and environmental sustainable development. The advance or delay of the key phenological period of wheat is a direct response to changes in climate conditions and production management measures. Therefore, accurately grasping the wheat phenology information is very necessary for wheat growth monitoring, which is also one of the requirements of precision agriculture.

[0004] However, the traditional wheat growth monitoring at the present stage mainly through field manual observation and remote sensing detection methods all have corresponding use restrictions. The former directly observes crops in the field, which is time-consuming and laborious, and has difficulties and uncertainties in long-term and large-area monitoring of crop phenology changes. Although remote sensing technology can reflect the growth rhythm characteristics of crops at a large scale, it can also provide important supplement for manual phenology observation data of ecological stations. However, this technology is difficult to eliminate the influence of atmospheric interference such as clouds and heavy rain. In addition, due to the limited sampling time resolution of remote sensing images, it is still difficult to obtain the corresponding crop information in a short time interval.

[0005] Therefore, based on the GNSS-IR technology, reliable and high spatio-temporal resolution wheat information monitoring such as wheat phenological nodes can provide an effective supplement for the ability of remote sensing to detect and monitor ground objects, and wheat growth monitoring and yield estimation have important use value and research significance in crop disaster monitoring and food security. SUMMARY

[0006] The present application is aimed at the deficiencies of the above background art, and provides a wheat phenological node information estimation method based on GNSS-IR technology, which realizes rapid monitoring of wheat phenological information, alleviates the technical problems of traditional manual methods being time-consuming and laborious, and insufficient spatio-temporal resolution in remote sensing, and is an effective supplement to existing crop monitoring technologies.

[0007] TECHNICAL SOLUTION The present application provides a wheat phenological node information estimation method based on GNSS-IR technology, which specifically comprises the following steps:

[0008] (1) processing the pre-acquired satellite observation and orbit data to obtain standard format data;

[0009] (2) obtaining signal-to-noise ratio multipath reflection signals and signal-to-noise ratio reflection components in a low elevation angle range for each system and frequency;

[0010] (3) performing time-frequency conversion on the obtained reflection components to extract frequency characteristics corresponding to frequency peaks;

[0011] (4) introducing an amplitude attenuation factor to construct a signal-to-noise ratio reflection signal model;

[0012] (5) fitting the signal-to-noise ratio model using a robust nonlinear least squares algorithm to obtain a damping coefficient;

[0013] (6) extracting wheat phenological period information, and using the turning point of the long time series curve of the damping coefficient as the demarcation point of the wheat phenological period.

[0014] Further, the step (1) is implemented as follows:

[0015] (11) reading the Rinex data, precise orbit data or satellite orbit data to extract the required data, including signal-to-noise ratio, satellite number PRN, azimuth angle, elevation angle, sampling time and other satellite basic information, and converting them into standard format;

[0016] (12) dividing the signal-to-noise ratio complete arc segment into ascending and descending arc segments according to the satellite system, number and frequency;

[0017] (13) monitoring whether the signal-to-noise ratio data exists missing and abnormal according to the satellite system, number and frequency, if the number of consecutive missing epochs is less than the threshold value 10, performing interpolation processing; otherwise, removing;

[0018] (14) According to whether the length of the arc segment where the signal-to-noise ratio data is detected by the satellite system is greater than a threshold value 1 h, if the arc segment is too short and less than the threshold value, it is removed.

[0019] Further, the step (2) is implemented as follows:

[0020] (21) According to the satellite system, the number, and the frequency, the signal-to-noise ratio data in an arc segment is extracted, and the unit is converted from dB-Hz to volts-volts:

[0021] In the formula, SNR v / v is the converted signal-to-noise ratio value, SNR dB-Hz is the original signal-to-noise ratio observation value;

[0022] (22) The direct signal is deducted from the original signal-to-noise ratio arc segment by using a second-order polynomial to fit the signal-to-noise ratio direct signal, to obtain a reflected signal component:

[0023] SNR r = SNR v / v - SNR d (2) In the formula, SNR d is the direct signal, and SNR r is the obtained reflected signal component;

[0024] (23) The data is selected according to the station environment by using SG filtering to reduce noise and selecting a 5°-25° elevation angle and a suitable azimuth angle range:

[0025]

[0026] In the formula, H is the data amount in the window, k is the number of days, w is the window size, i is the i th data in the window, h i is a smoothing coefficient obtained by least square fitting polynomial.

[0027] Further, the step (3) is implemented as follows:

[0028] The signal-to-noise ratio reflected component in the intercepted low elevation angle range is subjected to time-frequency conversion by using the Lomb-Scargle spectrum analysis method, the frequency characteristic f corresponding to the frequency peak value is extracted, and it is converted into a prior reflection height:

[0029]

[0030] In the formula, λ is the signal wavelength, and h is the reflection height.

[0031] Further, the step (4) is implemented as follows:

[0032] On the basis of the traditional signal-to-noise ratio model, an amplitude attenuation factor is introduced to construct a signal-to-noise ratio reflection signal model:

[0033]

[0034] In the formula, A represents the amplitude, is the phase, h0 is the median of the h long time series, D is the amplitude attenuation coefficient, k represents the wave number, and Λ represents the damping coefficient;

[0035] Let t=sinθ, f=2h0 / λ, and formula (5) is converted into a standard cosine function:

[0036]

[0037] Further, the step (6) is implemented as follows:

[0038] (61) The damping coefficient time series is normalized and denoised by using a sliding window filter to reduce its gross error interference and ensure the smoothness of the time series curve:

[0039]

[0040]

[0041] In the formula, Λ represents the original damping coefficient time series; represents the average of the top 20% maximum values in the time series; Λ normal is the normalized result; windowsize is the size of the sliding window, Λ filter (n) is the filtered result;

[0042] (62) Based on the 250 m resolution NDVI data of MODIS, the SG filter is used to smooth the curve, and the wheat phenological nodes including the green-up period, the heading period, and the harvest period are estimated by calculating the extreme points of the curvature of the Logistic simulation curve:

[0043]

[0044] In the formula, t is time, f(t) is the NDVI value in time, a and b are fitting coefficients, d is the initial background value, and c+d is the maximum value of NDVI;

[0045]

[0046] In the formula, K is the curvature of the Logistic simulation curve, z=exp(a+bt), dα is the moving angle of the tangent line when passing through a unit arc length along the time curve, and ds is the arc length.

[0047] (64) The turning point of the damping coefficient long time sequence curve is used as the demarcation point of the wheat phenology:

[0048] Turning point, Inflect s The turning point of the early stage of wheat growth, which represents the jointing stage; Inflect e The turning point of the middle and late stages of wheat growth, which is the mature stage; it satisfies that the second derivative is 0 and the value of the second derivative changes from negative to positive;

[0049] Minimum point, Min s The minimum value of the middle stage of wheat growth, which represents the heading stage; Min e The minimum value of the late stage of wheat growth, which represents the harvest stage; it satisfies that the first derivative is 0 and it is the critical point when the first derivative changes from negative to positive.

[0050] Beneficial effects: Compared with the prior art, the beneficial effects of the present application: Compared with the prior art, the present application adopts a strict data quality screening strategy, based on Lomb-Scargle spectrum analysis and nonlinear least squares algorithm, realizes signal-to-noise ratio model fitting and extracts corresponding characteristic parameters, that is, it improves the time and spatial resolution of wheat phenology node information, and reduces the precision loss and uncertainty caused by manual collection. BRIEF DESCRIPTION OF DRAWINGS

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

[0052] Figure 2 is a flowchart for obtaining standard format data;

[0053] Figure 3 is a flowchart for obtaining signal-to-noise ratio multipath reflection signals by system and frequency. DETAILED DESCRIPTION

[0054] The present application will be further described in detail below in combination with the drawings.

[0055] The application provides a wheat phenology node information estimation method based on GNSS-IR technology, through reading observation value files of a station on the day, satellite orbit products and the like, considering differences of satellite navigation systems such as GPS, BeiDou Navigation Satellite System (BDS) and the like, performing strict data preprocessing work including data screening, arc length selection and the like. And through calculation of a Fresnel zone, selecting a suitable satellite azimuth angle and elevation angle range, introducing Lomb-Scargle spectrum analysis and a nonlinear least square algorithm, constructing a signal-to-noise ratio model, fitting a signal-to-noise ratio curve and extracting a characteristic parameter: a damping coefficient, using a turning point of a long time sequence curve of the damping coefficient as a demarcation point of a wheat phenology period, and taking an average value of arc segment calculation results of multiple satellites of a single system to reduce uncertainty. Figure 1 As shown in the following steps:

[0056] Step 1: processing the satellite observation and orbit data obtained in advance to obtain standard format data, as shown in the following: Figure 2

[0057] Through reading Rinex data, precise orbit data or satellite orbit data, extracting required data including signal-to-noise ratio, satellite number PRN, azimuth angle, elevation angle, sampling time and the like, and converting the same into a standard format, the data is convenient for post-processing. According to the satellite system, number and frequency, the complete arc segment of the signal-to-noise ratio is segmented into ascending and descending arc segments, which is convenient for post-processing. According to the satellite system, number and frequency, the signal-to-noise ratio data is monitored for missing and abnormality, if the number of continuous missing epochs is less than a threshold value 10, interpolation processing is performed; otherwise, the data is removed; according to the system and frequency, the length of the arc segment where the signal-to-noise ratio data is located is detected, if the arc segment is too short and less than a threshold value, the data is removed.

[0058] Step 2: obtaining signal-to-noise ratio multipath reflection signals and signal-to-noise ratio reflection components in a low elevation angle range according to the system, frequency and the like, as shown in the following: Figure 3 The specific implementation process is as follows:

[0059] According to the satellite system, number and frequency, the signal-to-noise ratio data in an arc segment is extracted, and the unit is converted from dB-Hz to volts-volts:

[0060]

[0061] In the formula, SNR v / v is the converted signal-to-noise ratio value, and SNR dB-Hz is the original signal-to-noise ratio observation value.

[0062] ​Further processing the obtained standard format data, using second order polynomial fitting signal-to-noise ratio direct signal, subtracting the direct signal from the original signal-to-noise ratio arc segment to obtain the reflection signal component:

[0063] SNR r =SNR v / v -SNR d (2) In the formula, SNR d is the direct signal, and SNR r is the obtained reflection signal component.

[0064] Subsequently, noise reduction is performed using SG filtering (SG, Savitzky Golay), and data in the height angle range of 5°-25° and the appropriate azimuth angle range according to the station environment are selected:

[0065]

[0066] In the formula, H is the data amount in the window, k is the number of days, w is the window size, i is the i th data in the window, h i is the smoothing coefficient, and the least square fitting polynomial is obtained.

[0067] Step 3: Time-frequency conversion is performed on the obtained reflection component to extract the frequency feature corresponding to the frequency peak value.

[0068] Through Lomb-Scargle spectrum analysis method, time-frequency conversion is performed on the signal-to-noise ratio reflection component in the low height angle range intercepted in the previous step to extract the frequency feature corresponding to the frequency peak value: f, and the prior reflection height is converted:

[0069]

[0070] In the formula, λ is the signal wavelength, and h is the reflection height.

[0071] Step 4: An amplitude attenuation factor is introduced to construct a signal-to-noise ratio reflection signal model.

[0072] According to formula (5), the traditional signal-to-noise ratio model is expanded, an amplitude attenuation factor is introduced, and thus a signal-to-noise ratio reflection signal model is constructed:

[0073]

[0074] In the formula, A represents the amplitude, is the phase, h0 is the median of the long time series of h (in the case of bare soil), D is the amplitude attenuation coefficient, k represents the wave number, and Λ represents the damping coefficient.

[0075] If further t=sinθ, f=2h0 / λ, formula (5) can be converted into a standard cosine function:

[0076]

[0077] Step 5: Obtain the damping coefficient by fitting the signal-to-noise ratio model with the robust nonlinear least squares algorithm.

[0078] The signal-to-noise ratio model is fitted with the robust nonlinear least squares algorithm, i.e. formula (5), where the unknown parameters include A, and Λ. Thus, the signal-to-noise ratio parameter, the damping coefficient, is obtained.

[0079] Step 6: Extract the wheat phenological information, and use the turning point of the long-time series curve of the damping coefficient as the demarcation point of the wheat phenological stage.

[0080] First, the damping coefficient time series is normalized and denoised by sliding window filtering to reduce gross error interference and ensure the smoothness of the time series curve:

[0081]

[0082]

[0083] In the formula, Λ represents the original damping coefficient time series; represents the average of the top 20% of the maximum values in the time series; Λ normal is the normalized result; windowsize is the size of the sliding window, Λ filter (n) is the filtered result.

[0084] Since formula (5) has various unmodeled elevation angle related errors, such as antenna gain, which affect the calculation results of the damping coefficient. However, considering that the antenna gain pattern changes with time is constant and can be considered as a constant, and the surface roughness and dielectric properties change with crop growth, the damping coefficient will also change significantly during the entire wheat growth stage.

[0085] First, based on the 250m resolution NDVI data of MODIS, the SG filter is used to smooth the curve, and the wheat phenological nodes are estimated by calculating the curvature extreme points of the Logistic simulation curve, including the green-up stage, the heading stage, and the harvest stage:

[0086]

[0087] In the formula, t is the time, f(t) is the NDVI value within the time, a and b are the fitting coefficients, d is the initial background value, and c+d is the maximum value of NDVI.

[0088]

[0089] In the formula, K is the curvature of the Logistic simulation curve, z = exp (a + bt), dalpha is the moving angle of the tangent line when the curve passes through a unit arc length along time, and ds is the arc length.

[0090] By observing the long time sequence of the damping coefficient and comparing it with the actual wheat monitoring data, it is found that the turning point of the damping coefficient time sequence is close to the NDVI estimation result, that is:

[0091] (1) Turning point. Inflect s The turning point is the turning point of the early stage of wheat growth, that is, it represents the jointing stage. e The turning point is the turning point of the middle and late stages of wheat growth, that is, it represents the mature stage. At the same time, the second derivative is 0 and the value of the second derivative changes from negative to positive.

[0092] (2) Minimum point. Min s The minimum point is the minimum value of the middle stage of wheat growth, that is, it represents the heading stage. e The minimum point is the minimum value of the late stage of wheat growth, that is, it represents the harvest stage. The first derivative is 0, and it is the critical point where the first derivative changes from negative to positive.

[0093] Based on the above curve change points, that is, the turning point and the minimum point, the jointing stage, the heading stage and the mature stage of the wheat phenology can be obtained, and then the key and reliable phenology information for the near real-time monitoring and yield estimation of wheat can be provided.

[0094] In summary, the present application adopts strict data quality control strategy, based on Lomb-Scargle spectrum analysis, nonlinear least squares algorithm and extended signal-to-noise ratio fitting model, realizes the rapid estimation of the phenology in the growth cycle of wheat within 1-5km 2 around the station, which not only improves the time and spatial resolution of the wheat phenology node information, but also reduces the precision loss and uncertainty caused by manual collection.

Claims

1. A method for estimating information of a wheat phenological node based on GNSS-IR technology, characterized in that, The method comprises the following steps: (1) processing the pre-acquired satellite observation and orbit data to obtain standard format data; (2) obtaining signal-to-noise ratio multipath reflection signals and signal-to-noise ratio reflection components in a low elevation angle range for each system and frequency; (3) performing time-frequency conversion on the obtained reflection components to extract frequency characteristics corresponding to frequency peaks; (4) introducing an amplitude attenuation factor to construct a signal-to-noise ratio reflection signal model; (5) fitting the signal-to-noise ratio model by using a robust nonlinear least squares algorithm to obtain a damping coefficient; (6) extracting wheat phenological period information, and using a turning point of a long-time sequence curve of the damping coefficient as a demarcation point of the wheat phenological period; The step (3) is implemented as follows: By means of a Lomb-Scargle spectrum analysis method, the signal-to-noise ratio reflection components in the low elevation angle range are subjected to time-frequency conversion to extract frequency characteristics f corresponding to frequency peaks, and the frequency characteristics f are converted into a prior reflection height h0: In the formula, λ is a signal wavelength, and h is a reflection height; The step (4) is implemented as follows: The signal-to-noise ratio reflection signal model is constructed by extending a conventional signal-to-noise ratio model and introducing an amplitude attenuation factor: where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude, where A represents the amplitude Let t=sinθ and f=2h0 / λ, and the formula (5) is converted into a standard cosine function: 2.The wheat phenological node information estimation method based on GNSS-IR technology according to claim 1, characterized in that, The step (1) is implemented as follows: (11) the required data are extracted by reading Rinex data, precise orbit data or satellite orbit data, including signal-to-noise ratio, satellite number PRN, azimuth angle, elevation angle, sampling time and other satellite basic information, and the data are converted into a standard format; (12) the signal-to-noise ratio complete arc segments are segmented into ascending and descending arc segments according to the satellite system, number and frequency; (13) whether the signal-to-noise ratio data are missing or abnormal is monitored according to the satellite system, number and frequency, if the number of consecutive missing epochs is less than a threshold value 10, interpolation processing is performed; otherwise, the data are removed; (14) whether the length of an arc segment in which the signal-to-noise ratio data are located is greater than a threshold value 1h is detected, if the arc segment is too short and less than the threshold value, the data are removed. 3.The wheat phenological node information estimation method based on GNSS-IR technology according to claim 1, characterized in that, The step (2) is implemented as follows: (21) the signal-to-noise ratio data in an arc segment are extracted according to the satellite system, number and frequency, and the unit of the data is converted from dB-Hz to volts-volts: where SNR v / v is the transformed signal-to-noise ratio value, SNR dB-Hz is the original signal-to-noise ratio observation; (22) the signal-to-noise ratio direct signal is fitted by using a second-order polynomial, the direct signal is deducted from the original signal-to-noise ratio arc segment to obtain a reflection signal component: SNR r = SNR v / v - SNR d (2) where SNR d is the direct signal, SNR r is the resulting reflected signal component; (23) the data are selected according to the station environment by using SG filtering to reduce noise and selecting a 5°-25° elevation angle and a suitable azimuth angle range: where H is the amount of data within the window, k is the number of days, w is the window size, i is the ith data within the window, h i is the smoothing coefficient, obtained from the least square fitting polynomial.

4. The wheat phenological node information estimation method based on GNSS-IR technology according to claim 1, characterized in that, The step (6) is implemented as follows: (61) the damping coefficient time sequence is subjected to normalization processing and noise reduction processing by using a sliding window filtering to reduce gross error interference and ensure that the time sequence curve is smooth: where Λ represents the original damping coefficient time series; represents the average of the top 20% largest values in the time series; Λ normal is the normalized result; windowsize is the size of the sliding window, Λ filter (n) is the filtered result; (62) based on 250m resolution NDVI data of MODIS, the curve is smoothed by using SG filtering, and wheat phenological nodes including a green-up period, a heading period and a harvest period are estimated by calculating a Logistic simulation curve curvature extreme point: In the formula, t is time, f(t) is NDVI value in time, a and b are fitting coefficients, d is initial background value, and c+d is maximum NDVI value; In the formula, K is curvature of a Logistic simulation curve, z=exp(a+bt), dα is moving angle of a tangent line when passing unit arc length along the time curve, and ds is arc length; (63) Using turning points of long time series curves of damping coefficients as demarcation points of wheat phenological stages: Inflect s is the turning point of the early stage of wheat growth, i.e. represents the jointing stage; Inflect e is the turning point of the middle and late stage of wheat growth, i.e. is the mature stage; simultaneously satisfies that the second derivative is 0 and the second derivative value changes from negative to positive; Min s Min is the minimum value of the wheat growth in the middle stage, which represents the heading stage. e Min is the minimum value of the wheat growth in the late stage, which represents the harvest stage. It satisfies the first derivative being 0 and is the critical point where the first derivative changes from negative to positive.

Images