Monitoring of Wheat Height Based on Multi-GNSS Reflected Signals

Abstract

Global Navigation Satellite System interferometric reflectometry (GNSS-IR), a new and inexpensive technique, has become available to the broader scientific community for detecting surface environmental information, such as soil moisture, snow depth and vegetation growth. However, there have been limited experiments focusing on the potential of crop height retrieval, especially the performance evaluation of BeiDou Navigation Satellite System (BDS) with other GNSS. Accuracy and reliability are challenging to achieve with traditional methods utilizing a single GNSS, and few measured verification data. In this study, an improved method that includes segmentation processing and multi-GNSS fusion is proposed based on GPS/GLONASS/Galileo/BDS multi-frequency data. Furthermore, experiments were carried out on a farmland in Fengqiu County, Henan Province, China. The results show that the height retrievals from four GNSS were in good agreement with the in situ observations during the whole growth cycle of the wheat after overwintering. Meanwhile, the retrievals based on the proposed method exhibited greater correspondence than the single frequency results, the correlation coefficient was increased and the root-mean-square error (RMSE) was reduced, respectively. Therefore, this study illustrates the feasibility of the proposed method to precisely estimate wheat height and its potential for use in the early warning of wheat lodging based on GNSS-IR.

1. Introduction

Wheat is a type of cereal crop that is widely planted all over the world, and is a staple food for human beings. Yield information comprehensively reflects the growth level of crops in that year, and has important reference value for agricultural food policy and farming planning [1]. The height of wheat is the main index reflecting the growth and development of wheat, and it is a critical factor for predicting the yield during the growth process. At the same time, due to global warming, various extreme weather conditions occur frequently, which may affect the growth of wheat. For example, lodging caused by wind, rain and other meteorological factors threatens the growth of wheat seriously affect the maturity of wheat, reduce production and endanger agricultural production and food security. The most obvious sign of wheat lodging is an unreasonable change in wheat height. Therefore, the measurement or real-time monitoring of wheat height during growth is of great significance for yield prediction and for the early warning of disasters, such as lodging [2].

For a long time, people used the meter ruler to measure the height of wheat. This method is time-consuming, laborious and difficult to use for continuous monitoring. As state-of-the-art technology, the Synthetic Aperture Radar (SAR) and Light Detection and Ranging (LiDAR) have been used more recently to carry out wheat height monitoring. Kumar et al. [3] employed Sentinel-1A Synthetic Aperture Radar analysis of time series data at C-band to estimate winter wheat crop growth parameters, such as leaf area index, vegetation water content, fresh biomass, dry biomass and plant height, using multiple machine learning algorithms. A ground-based multi-sensor phenotyping system equipped with ultrasonic sensors and LiDAR was developed by Yuan et al. [4] to compare it with manual measurements and an unmanned aircraft system. LiDAR provided the best results; the unmanned aircraft system obtained reasonable ones and the ultrasonic sensors did not perform well. These abovementioned methods also have individual restrictions, such as being expensive or having low spatial–temporal resolution [5]. In recent years, GNSS reflected signals have been used opportunistically to retrieve land surface environmental information instead of navigation and positioning. The first use of GNSS-reflected signals on remote sensing was proposed by Hall et al. [6] to determine ocean wind velocity. Subsequently, GNSS reflectometry (GNSS-R) studies primarily pursued the application potential of reflections in the ocean for measurements, such as sea surface wind speed [7,8], ocean altimetry [9,10] and sea ice thickness [11,12]. The above studies focused primarily on using specially designed GNSS receivers to retrieve environmental parameters, which is not conducive to the promotion of the technology. Later, a subset of studies also had success in using geodetic-quality GNSS receivers with a Right-Hand Circularly Polarized (RHCP) antenna to estimate surface environmental monitoring information, such as soil moisture [13,14,15,16,17], changes in vegetation [18,19] and snow depth [20,21]. This technology, known as GNSS interferometric reflectometry (GNSS-IR), was initially proposed by Larson et al. [22] to measure soil moisture fluctuations using signal-to-noise ratio (SNR) interferograms.

The application of GNSS-IR technology in agriculture has always been concerned with measurements, such as soil moisture, vegetation moisture, vegetation index and crop height. Chew et al. [23] studied the relationship between soil moisture and the amplitude, frequency and phase shift of SNR multipath signals. The results showed that phase shift was the most ideal factor to describe fluctuation in soil moisture. Later, Larson and Wan found in the test that the amplitude of SNR had strongly linear relationship with the vegetation water content (VWC) when VWC was less than ~1 kg/m² [18]. Chew et al. [24] verified that the normalized amplitude time series kept below 0.78 for a long period was possibly affected by vegetation. The normalized microwave reflection index (NMRI), calculated from multipath data, was also able to characterize the changes in surface vegetation in [25]. Similarly, the amplitude of the annual growth cycle and the correlation between the NMRI and the normalized difference vegetation index (NDVI) were estimated from 184 GPS sites by Evans et al. [26]. At the same time, considering the impact of vegetation on GNSS reflected signals, it was necessary to mitigate this effect to better estimate soil moisture [27]. A retrieval method of vegetation biomass based on GNSS SNR was proposed by Zheng et al. [28]. Techniques to estimate the vegetation optical depth and infer the VWC using a single GNSS receiver located under the vegetation layer were analyzed [29,30].

In comparison with other aspects, there are few studies on crop height at present. A vegetation-height retrieval algorithm was first proposed using the SNR-dominant period with GPS L1 frequency by Zhang et al. [31], who also found that these vegetation height retrievals were sensitive to changes in crop height of at least one wavelength. Similarly, Zhang et al. [32] comprehensively evaluated the potential of BDS and GPS to retrieve soil moisture and vegetation growth. Empirical mode decomposition (EMD) and wavelet algorithms were compared to estimate the optimal wheat height changes. The conclusion was that the wavelet analysis method showed optimal retrieval ability during the height change process of wheat in the middle and late periods. Li et al. [33] proposed a new retrieval method that used a low-cost GNSS chip and an RHCP antenna to retrieve the peak frequencies of the GNSS SNR series using Lomb–Scargle periodogram (LSP) and simplified the winter wheat as multi-layer equivalent mediums. The results showed that the estimated equivalent medium heights were very close to the in situ observations during the wheat heading and ripening stages.

In previous works, vegetation height retrieval was mostly focused on using GPS- or BDS reflected signals. There are few experiments based on the GLONASS, Galileo or BDS to retrieve the crop growth information such as wheat height in a farmland environment. In addition, the traditional method can only retrieve information in cases where the wheat heights are greater than one wavelength, and also cannot distinguish between bare soil and a height of less than one wavelength. Meanwhile, there are few in situ measured data of the whole growth cycle, and it is difficult to comprehensively compare the retrieval differences between different GNSS and frequencies. Therefore, this study aims to retrieve wheat height measurements based on GPS, GLONASS, Galileo and BDS SNR data obtained with a classical geodetic receiver over an intensively cultivated wheat field in the northeast Henan Province. The differences of individual system and frequency signals on wheat height retrieval are deeply analyzed, and segment processing is performed to obtain optimum retrieval results. An improved method is then proposed to fuse multi-system and multi-frequency data to improve the retrieval accuracy. In addition, the retrieval values are compared with in situ observations for systematic evaluation, resulting in an efficient method for monitoring or providing early warning of wheat height change with high spatial and temporal resolution.

2. Methods

2.1. SNR Characteristics

When the signals generated by GNSS satellites through the atmosphere are reflected by the surface and other ground environment in the process of propagation (as in Figure 1), the signal superposition at the antenna will produce relatively stable interference signals, which is more obvious at the low elevation angle in Figure 2. According to the full model [34], interference SNR data can be given by:

$${\mathrm{SNR}}^2{=\mathrm{A}}_{\mathrm{d}}^2{+\mathrm{A}}_{\mathrm{m}}^2{+2\mathrm{A}}_{\mathrm{d}}{\mathrm{A}}_{\mathrm{m}}\cos \mathsf{\varphi }$$

(1)

where A_d and A_m are the amplitude of the direct and reflected signal, respectively, and $\mathsf{\varphi }$ represents the phase difference between the direct and reflected signals.

The influence of the power of the direct signals is seen in a large arcing trend in the SNR data in Equation (1), which does not contain information about the surface environment. In order to remove this influence and only see the interference pattern, the SNR is firstly converted from dB-Hz to a linear scale (volt/volt) and then detrended with a two-order polynomial fitting, estimated for each satellite track on each day:

$${\mathrm{SNR}}_{\mathrm{V}/\mathrm{V}}{=10}^{{\mathrm{SNR}}_{\mathrm{dB}-\mathrm{Hz}}/20}$$

(2)

where ${\mathrm{SNR}}_{\mathrm{dB}-\mathrm{Hz}}$ represents the raw data and ${\mathrm{SNR}}_{\mathrm{V}/\mathrm{V}}$ denotes the transformed data.

Therefore, the detrended SNR data, that is, the reflected components, have previously been modeled using the following equation [13]:

$$\mathrm{SNR}=\mathrm{A}\cos (\frac{4\mathsf{\pi }\mathrm{h}}{\mathsf{\lambda }}\sin \mathsf{\theta }+\mathsf{\phi })$$

(3)

where θ represents the elevation angle of the satellite, A is the constant amplitude, λ refers to the GNSS signal wavelength (GPS λ_L1 = 0.1902 m, Galileo λ_E1 = 0.1904 m, BDS λ_B1 = 0.1920 m), φ denotes a phase shift and h is the height of the antenna, which is generally replaced by the median of the high reflectance series calculated by the Lomb–Scargle periodogram.

Generally, Equation (2) forms the foundation for remote sensing of soil moisture using GNSS-IR. The parameters in Equation (2), A and φ, are two metrics that have been investigated with regards to their potential for inferring changes in vegetation water content and surface soil moisture.

2.2. Wheat Height Retrieval

Larson et al. proposed an algorithm to retrieve snow depth and sea level height, which demonstrates the potential of the GNSS-IR technique in estimating height changes. Similarly, when wheat grows, the vegetation surface gradually replaces the bare soil surface as the dominant reflecting surface (Figure 3). Consequently, the effective reflected height of the antenna above the reflecting surface decreases. Thus, this method is used to retrieve the wheat height. First, the LSP is applied to estimate the frequency f; then, the f can be converted to the reflected height as follows:

$${\mathrm{h}}_{\mathrm{veg}}=\frac{\mathsf{\lambda }\times \mathrm{f}}{2}$$

(4)

where f denotes the peak frequency of the LSP.

Zhang et al. [31] also proposed an algorithm to retrieve changes in wheat height using the dominant period T of SNR data by wavelet analysis. As the changes in h impact T, this property allows the use of changes in T to infer changes in h, further evaluating relative wheat height. This model is given by:

$${\mathrm{h}}_{\mathrm{veg}}=\frac{\mathsf{\lambda }}{2\times \cos \mathsf{\theta }\times \frac{\mathrm{d}\mathsf{\theta }}{\mathrm{dt}}+\mathrm{T}}$$

(5)

where T is the dominant period of the SNR reflected components and dθ/dt represents the change rate of satellite elevation angle, θ.

However, the selection of wavelet basis will seriously affect the periodic structure of the SNR reflected signals, and the best elevation angle cosθ may change due to different GNSS and frequencies. Therefore, we do not verify the performance of this method, only using Equation (5).

Finally, considering the penetration characteristics of the GNSS signals, the wheat height of multi-satellites (m = 1,2,3…) is retrieved as follows:

$${\mathrm{H}}_{\mathrm{veg}}{=\mathrm{h}}_0-{\mathrm{h}}_{\mathrm{veg}}$$

(6)

$$\mathrm{H}=\frac{{{\displaystyle \sum }}_{\mathrm{m}=1}^{\mathrm{N}}{\mathrm{H}}_{\mathrm{veg}-\mathrm{m}}}{\mathrm{N}}+\mathsf{\lambda }$$

(7)

where h₀ is the median value of the top 15% h data during the wheat growth period, H_veg-m denotes the changes in height and H is refers to final wheat height. Due to the extra wavelength, this method only retrieves wheat height greater than one wavelength.

2.3. Retrieval Steps

The algorithm flow is shown in Figure 4. According to the surrounding environment of the site, the SNR data within the azimuth range of 30–300° and the elevation angle range of 5–25° are selected first. Quality control is then used in the LSP periodogram analysis to determine the proper frequency. If the SNR passes detection, the relevant metrics will be directly calculated; otherwise, the signal will be decomposed through the improved complete ensemble EMD (ICEEMDAN) [35] to obtain the effective frequency information. The SNR metrics A and f are subsequently calculated through a non-linear least square fitting procedure and the LSP.

In addition, the wheat height is retrieved during wheat growth periods by the multi-frequency of the GPS/GLONASS/Galileo/BDS, and the results are systematically analyzed and evaluated at the site. Segmented processing is then carried out to correct the deviation with each GNSS system. Following this, the multi-frequency and multi-system retrieval results are fused to improve the retrieval accuracy and to compensate the shortcomings of traditional methods. Finally, the retrieval values are compared with the in situ observations and a growing degree-days (GDD) model simulation [36], which verify the effectiveness of this method. The GDD model is described in detail in [37].

2.3.1. Quality Control

For the wheat height retrieval, the vegetation canopy changes and the GNSS site environment play an important role in the measurement results. The response of different frequency bands to the field environment is very complex, so quality control is a critical step in a wheat farmland environment. Therefore, two methods are used in periodogram analysis to determine the correct frequency f.

The first type of ratio refers to the ratio of the amplitude corresponding to the highest peak to the average of all other amplitudes:

$${\mathrm{Ratio}}_1{=\mathrm{AMP}}_{\max }{/\mathrm{mean}(\mathrm{AMP}}_{})$$

(8)

where AMP_max is the maximum value of the amplitude sequence in the LSP and AMP denotes the amplitude sequence of the other peaks in the LSP.

The second ratio is given as follows:

$${\mathrm{Ratio}}_2={\mathrm{AMP}}_{\max }/\sec \max ({\mathrm{AMP}}_{})$$

(9)

where secmax(AMP) is the amplitude of the second highest peak in the LSP.

In the data processing, after setting the thresholds of the two parameters, if the values of Ratio₁ and Ratio₂ of the LSP do not exceed the a priori threshold, the frequency f is directly obtained. Otherwise, the ICEEMDAN method [35] will be used to decompose the SNR data, and the fourth layer will be selected for the next step.

2.3.2. Segmented Processing

It can be seen from previous works that it is difficult to distinguish wheat height of less than one wavelength using the traditional method [31,32], which seriously restricts the practical application of GNSS-IR. Therefore, a segmented processing method is proposed to address this issue. Specifically, the wheat height time series is segmented using the normalized amplitude (A_norm) and empirical phenological information [24].

First, whether A_norm is less than 0.78 can be approximately regarded as an indicator of vegetation effect, the height retrieval time series is divided into three sections: before the regreening stage, before the rapid growth stage and after harvest. One wavelength would not be added in the first and last stages to obtain more practical retrieval results for all GNSS frequencies.

Second, the L1 frequency data of each system would also not increase one wavelength from the heading stage to harvest due to the signal attenuation. The winter wheat heading date in the Fengqiu Experimental Station (HNFQ) site is roughly the same, basically around day of year (DOY) 115 according to our investigation.

Therefore, the final H after segment processing can be given as follows:

$$\mathrm{H}=\left\{\begin{array}{c}\frac{{{\displaystyle \sum }}_{\mathrm{m}=1}^{\mathrm{N}}{\mathrm{H}}_{\mathrm{veg}-\mathrm{m}}}{\mathrm{N}}, \mathrm{DOY}\le \mathrm{Day}1\\ \frac{{{\displaystyle \sum }}_{\mathrm{m}=1}^{\mathrm{N}}{\mathrm{H}}_{\mathrm{veg}-\mathrm{m}}}{\mathrm{N}}+\mathsf{\lambda }, \mathrm{Day}1<\mathrm{DOY}<\mathrm{Day}2,\mathrm{except} \mathrm{L}1 \mathrm{frequency}\\ \frac{{{\displaystyle \sum }}_{\mathrm{m}=1}^{\mathrm{N}}{\mathrm{H}}_{\mathrm{veg}-\mathrm{m}}}{\mathrm{N}}, \mathrm{Day}1<\mathrm{DOY}<\mathrm{Day}2,\mathrm{for} \mathrm{L}1 \mathrm{frequency}\\ \frac{{{\displaystyle \sum }}_{\mathrm{m}=1}^{\mathrm{N}}{\mathrm{H}}_{\mathrm{veg}-\mathrm{m}}}{\mathrm{N}}, \mathrm{DOY}\ge \mathrm{Day}3\end{array}\right.$$

(10)

where Day1 is the date when A_norm is first less than 0.78, Day2 denotes the date of heading stage, which is DOY 115 in this study, and Day3 represents the last date when A_norm is less than 0.78.

2.3.3. Multi-Frequency and Multi-System Fusion Retrieval

Considering the sensitivity of different frequency signals to vegetation, it is necessary to fuse multi-frequency and multi-system values to improve the retrieval accuracy. Therefore, an improved method based on the daily mean root-mean-square errors is proposed for fusion.

First, for a single system, the fusion values are the mean of different frequencies (n = 1,2,3):

$${\mathrm{H}}_{\mathrm{sys}}=\frac{{{\displaystyle \sum }}_{\mathrm{n}=1}^{\mathrm{F}}{\mathrm{H}}_{\mathrm{n}}}{\mathrm{F}}$$

(11)

where F denotes the number of frequencies for a single system. Then, the daily mean root-mean-square errors (i = 1, 2, 3, …) with the GDD model simulation values are calculated according to the fused results from each GNSS system (k = 1, 2, 3, 4):

$${\mathsf{\sigma }}_{\mathrm{i}}^{\mathrm{k}}{=\left(\right(\mathrm{H}}_{\mathrm{GDD}-\mathrm{i}}^{}{-\mathrm{H}}_{\mathrm{sys}-\mathrm{i}}^{\mathrm{k}}{)}^2{)}^{-1/2}$$

(12)

where H_GDD indicates the GDD model simulation and H_sys is the retrieval value of a single system.

The weights can be expressed as follows:

$${\mathrm{Q}}_{\mathrm{i}}^{\mathrm{k}}=\frac{{\mathsf{\sigma }}_{\mathrm{i}}^{\mathrm{k}}}{{{\displaystyle \sum }}_{\mathrm{k}=1}^4{\mathsf{\sigma }}_{\mathrm{i}}^{\mathrm{k}}}$$

(13)

Based on these, the final fusion retrieval values are defined as:

$${\mathrm{H}}_{\mathrm{fusion}-\mathrm{i}}={\displaystyle \sum }_{\mathrm{k}=1}^4{\mathrm{H}}_{\mathrm{sys}-\mathrm{i}}^{\mathrm{k}}{\times \mathrm{Q}}_{\mathrm{i}}^{\mathrm{k}}$$

(14)

3. Results

3.1. Experiments

The experiment was carried out at Fengqiu Experimental Station (HNFQ) for the Agroecology Chinese Academy of Sciences, Fengqiu County, Henan Province, China (Figure 5). A Sino M300 Pro GNSS receiver and AT500 GNSS antenna were used to track the GNSS signals from 1 January 2022 to 21 June 2022 (DOY 1–172). The SNR data of DOY 157–167 were lost due to power failure during wheat harvesting. The in situ wheat height measurement had been carried out at five-day interval since 5 March (DOY 64) 2022, measured with a metallic tape ruler from four sampling points in the wheat field. Then, the average of these heights was considered as the in situ wheat height. During the experiment, the measured wheat height at the beginning was 0.2730 m, and the maximum plant height measured on 25 April (DOY 115) was 0.7083 m. After that, the wheat height remained basically unchanged. Then, until DOY 156, when the wheat had been harvested, there was no obvious crop growth at the site.

For this experiment, the sample frequency of the GNSS receiver was set up to 15 Hz, and the types of SNR signals received by the GNSS receiver are shown in

Table 1. Received SNR types for HNFQ.

Satellite System	Frequency Band	SNR Type
GPS	L1	S1
	L2	S2
	L5	S5
GLONASS	G1	S1
	G2	S2
Galileo	E1	S1
	E5b	S7
BDS	B1-2	S2
	B3	S6
	B2b	S7

. Moreover, the low elevation angles (less than 5°) were excluded to avoid obstacle effects surrounding the wheat farmland.

3.2. Wheat Height Retrieval Results

3.2.1. Single Frequency

Figure 6 shows an example of the reflector height from the GPS L1 frequency band during the whole winter wheat growth period. It can be seen from Figure 6 that the reflector height has an obvious downward trend, which is consistent with the growth process of wheat. From January to February (DOY 1–60), the wheat becomes dormant, resulting in a basically consistent reflector height. Then, after March, the wheat regreens and begins to grow gradually. In the stem extension stage (after DOY 75), the wheat grows rapidly, and the reflector height also distinctly increased. In contrast, the wheat gradually heads on DOY 115, which forms a stable wheat reflective canopy and leads to the height remaining relatively unchanged. After the wheat harvest, its height recovers rapidly.

Specifically, the multi-frequency retrieval results from GPS/GLONASS/Galileo/BDS, together with the in situ measurements and the daily wheat height simulations using the GDD model, are illustrated in Figure 7, Figure 8, Figure 9 and Figure 10. The in situ measured heights are recovered 14 times, from 5 March to 13 May (DOY 133), with an interval of five days. In this experiment, the change trends of the wheat height retrievals of each frequency band and each GNSS are consistent with those of the small number of in situ observations and the GDD model simulations. Some characteristics of the wheat growth process are detected, such as rapid growth, remaining relative unchanged, and harvest. As discussed above, the winter wheat enters the rapid growth stage after DOY 60, and reaches a peak of approximately 0.7 m around DOY 120. After that, the wheat height remains relatively unchanged until DOY 156. During the DOY 168–172 period, the wheat height recovers rapidly after harvest. These results show that the in situ GNSS height retrievals are able to detect the change process in wheat height.

It is worth noting that the steady change trend in wheat height did not represent the actual wheat growth during the DOY 1–60 and DOY 166–172 periods, for which the values should be closer to zero. Generally speaking, it is almost impossible to retrieve wheat height of less than one wavelength due to the difficulty in distinguishing the surface of the wheat and the soil [30].

In addition, the retrievals of different frequency bands of different systems are also inconsistent. For the GPS and BDS systems, the retrieval height of L1 and B1-2 are clearly better than these of L2/L5 and B2b/B3. Similarly, for GLONASS and Galileo, the results of the first frequency band (G1/E1) are better than that of the second frequency band (G2/E5b). However, the L1 frequency retrievals, of which the wavelength is nearly 0.19 m, have the problem of being over high after heading stage. The differences between the retrieval and in situ values reach 0.1764, 0.1688, 0.1641 and 0.1599 m, respectively, which are close to one wavelength for these frequencies. One possible reason is the relationship between the signal penetration depth and the GNSS signal wavelength; that is, the higher the frequency, the larger the signal attenuation and the shallower the penetration depth. In early May, the wheat forms a new reflection layer on the top of the wheat crop, resulting in penetration by the GNSS signals. The L1 frequency band may stay at the surface layer of the canopy and may then be reflected, while the L2 frequency penetrates a certain thickness of wheat ears before being received.

3.2.2. Segmented Processing

The segmented processing method in Section 2.3.2 is used to solve the problem above. Figure 11 shows an example of the segmented processing results based on the GPS. Note that, compared with the retrieval results obtained directly, the results of the segmented processing method are more consistent with the in situ measured values and the GDD model simulations, especially the GPS L1 frequency, which is more in line with the changing trend of the actual winter wheat growth process. Compared with Figure 6, the retrieval values are closer to the in situ values during DOY 115–150.

The relative wheat height retrievals of the L1 frequency band using the segmented processing method are compared with the in situ height observations in

Table 2. Relative deviations of L1 frequency band between wheat height retrievals, in situ observations and GDD model.

Band	DOYs	In Situ/m	Original/m	Deviations_ori/m	Segmented/m	Deviations_seg/m
GPS-L1	115	0.7083	0.8557	0.1474	0.6655	−0.0428
	120	0.7000	0.8266	0.1266	0.6364	−0.0636
	124	0.6938	0.8477	0.1509	0.6545	−0.0393
	133	0.7000	0.8655	0.1655	0.6753	−0.0247
GLONASS-G1	115	0.7083	0.8491	0.1483	0.6694	−0.0389
	120	0.7000	0.8604	0.1442	0.6570	−0.0430
	124	0.6938	0.8691	0.1350	0.5416	−0.0522
	133	0.7000	0.8441	0.1378	0.6506	−0.0494
Galileo-E1	115	0.7083	0.8566	0.1408	0.6587	−0.0496
	120	0.7000	0.8442	0.1604	0.6700	−0.0300
	124	0.6938	0.8288	0.1753	0.6787	−0.0151
	133	0.7000	0.8378	0.1441	0.6537	−0.0463
BDS-B1-2	115	0.7083	0.8478	0.1395	0.6575	−0.0508
	120	0.7000	0.8372	0.1372	0.6469	−0.0531
	124	0.6938	0.8523	0.1585	0.6620	−0.0318
	133	0.7000	0.8423	0.1423	0.6520	−0.0480

. There are no in situ measured values at the beginning and after harvest; only a comparison of wheat height after the heading stage is shown here.

The retrieval performance of these two methods is inconsistent. The differences between the four in situ observations and the original retrievals H are obviously larger due to the segmented processing. Specifically, the deviation in the original method ranges from 0.12 m to 0.17 m, and the deviation in this proposed method is from −0.06 m to −0.01 m, exhibiting a negative bias. This is consistent with the previous conclusion: due to the penetration of the L-band signal through the wheat, the retrievals results are significantly smaller than the in situ observations generally. The segmented processing results are closer to the observed growth trend than the GDD simulations. Additionally, the proposed method presents a much better fit for the in situ measurements than the original retrievals.

3.2.3. Multi-Frequency and Multi-System Fusion

To make full use of the abundant GNSS signal multi-frequency and multi-system resources, the method described in Section 2.3.3 is used to fuse the retrievals. The fusion results are depicted in Figure 12.

The retrieval results based on the fusion method in Section 2.3.3 are in good agreement with the in situ observations and the GDD simulations, which fully proves the feasibility and accuracy of wheat height retrieval using multi-frequency and multi-system signals. In contrast, compared with the raw retrievals, the retrievals show great correspondence with the GDD model simulations in DOY 1–60 and 110–157, which proves that the application of GNSS signals for monitoring changes in wheat height has considerable potential. However, it is still necessary to collect more in situ wheat height observations to further analyze and verify the universality and accuracy of this technology.

4. Discussion

In the HNFQ site experiment, the errors between the wheat height retrieval and the GDD model simulations produced by four GNSS system with multi-frequency after segmented processing are illustrated in Figure 13, Figure 14, Figure 15 and Figure 16. Clearly, most of the residuals are concentrated within ±0.15 m, while those from the L5/G2/E5b/B2b frequency band are slightly larger, with a variation of ±0.20 m. The errors based on the L1/G1/E1/B1-2 frequency are slightly better than those of the other frequency, and have the highest accuracy, with most stable residuals smaller than ±0.08 m. In particular, for all frequencies, the residuals increase slowly during the DOY 1–60 period, which may be caused by the rising GDD simulation values and the relatively stable retrieval values. The residuals during the DOY 80–115 period are slightly larger than those at other growth stages. With the rapid wheat growth, the space between the winter wheat becomes larger, and the wheat is well penetrated by the GNSS signals, which explains the above phenomenon. In addition, after the segmented correction, the wheat height residuals of the L1 frequency significantly decrease at the heading stage. In general, this proposed method can effectively overcome the problems of over high and low accuracy.

In terms of correlation, the scatterplot of the wheat height retrieval at HNFQ site using multi-system and multi-frequency fusion is depicted in Figure 17. It shows that the retrievals are consistent with the trend of the in situ values, which are basically distributed on both sides of the regression line. Additionally, we obtained the following score values from the HNFQ site: RMSE = 0.0383 m, mean-absolute error (MAE) = 0.0301 m. The correlation coefficient of the proposed combination method is the largest, reaching 0.9555, which effectively overcomes the uncertainty discussed above and greatly improves the stability and accuracy of wheat height retrieval.

In addition, the performance of the proposed method is compared with that of other single frequency methods of four GNSS.

Table 3. Wheat height retrieval precision statistics.

Methods	Correlation Coefficient	Root-Mean-Square -Error (m)	Mean-Absolute-Error (m)
GPS-L1	0.8868	0.0765	0.0471
GPS-L2	0.6886	0.1053	0.0967
GPS-L5	0.7684	0.0950	0.0780
GLONASS-G1	0.7993	0.0864	0.0712
GLONASS-G2	0.6566	0.1022	0.0940
Galileo-E1	0.7891	0.0853	0.0668
Galileo-E5b	0.6557	0.1143	0.0997
BDS-B1-2	0.8092	0.0899	0.0755
BDS-B3	0.7348	0.1015	0.0907
BDS-B2b	0.7204	0.1068	0.0895
Fusion	0.9555	0.0383	0.0301

shows the accuracy in terms of HNFQ site error statistics based on the multiple retrieval methods.

It can be seen from Table 3 that all the height retrievals are in good agreement with the limited in situ observations, which means that the proposed method performs best. Most of RMSEs and MAEs are nearly 0.095 m and 0.081 m, respectively. Except for the GPS L1 frequency with 0.8868, the retrievals based on a single frequency are almost all below 0.85, proving the difference between the different frequencies based on the GNSS system, while the R2 of the proposed method is above 0.9, reaching 0.955 m. The RMSE and MAE compared with the single frequency results are decreased to a great extent. Specifically, the average correlation coefficient of the proposed method is increased by 27.25%, the average RMSE is decreased by 60.23% and the average MAE is decreased by 62.79%. That is, the height retrieval method can be applied to all GNSS system, especially the proposed multi-system and multi-frequency fusion method, taking full advantage of the difference and complementarity of GNSS signals with different frequencies to improve the retrieval accuracy. It is worth mentioning that due to the single experimental site, the applicability of the proposed model for wheat height retrieval needs to be further verified.

5. Conclusions

Different from traditional navigation, positioning and timing (PNT) services, GNSS has been widely explored as an efficient remote sensing tool, with GNSS-IR technology being used for surface environmental monitoring. Essentially, the contribution of this research is to evaluate the performance of wheat height retrieval using different GNSS. To fully exploit the advantages of multi-mode GNSS combination, a fusion method is proposed to retrieve winter wheat height using multi-frequency and multi-system SNR data and its performance is verified by using in situ measurements from the HNFQ site and simulated values from GDD model.

The long-term experiment in Henan station showed that the retrieval results of a single frequency from four GNSS (GPS/GLONASS/Galileo/BDS) can reflect the height change process of most wheat growth stages after overwintering. In terms of wheat height accuracy, the proposed method performed better than the traditional method using a single frequency. There was good correlation between the in situ values and the retrieval values, with a correlation coefficient of 0.9555, and a corresponding RMSE and MAE of 0.0383 m and 0.0301 m, respectively. These improvements mainly benefited from the data fusion of the multi-frequency and multi-system. Compared with the retrievals from a single frequency, the R2, RMSE and MAE improved on average by 27.25%, 60.23% and 62.79%, respectively. Therefore, retrieving wheat height based on the proposed method may have application in crop monitoring with high time resolution and accuracy.

In the future, we will continue to observe the Henan site for a long period to obtain more information about other relevant crops, such as corn, to confirm the potential and adaptability of this technology. At the same time, future research will focus on using GNSS SNR data to obtain other crop indicators, for example, the division of the growth cycle or crop yield.