Monitoring of Wheat Height Based on Multi-GNSS Reﬂected Signals

: Global Navigation Satellite System interferometric reﬂectometry (GNSS-IR), a new and inexpensive technique, has become available to the broader scientiﬁc 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 veriﬁcation 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 coefﬁcient 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.


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 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 multifrequency 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.

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: where A d and A m are the amplitude of the direct and reflected signal, respectively, and φ represents the phase difference between the direct and reflected signals.
Remote Sens. 2022, 14, x FOR PEER REVIEW 3 of 19 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 multisystem 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.

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: where Ad and Am are the amplitude of the direct and reflected signal, respectively, and ϕ 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: where SNR dB-Hz represents the raw data and SNR V/V denotes the transformed data. Therefore, the detrended SNR data, that is, the reflected components, have previously been modeled using the following equation [13]: 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.

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: where f denotes the peak frequency of the LSP. 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: where SNR dB-Hz represents the raw data and SNR V/V denotes the transformed data. Therefore, the detrended SNR data, that is, the reflected components, have previously been modeled using the following Equation [13]: 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.

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: 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: 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: where h0 is the median value of the top 15% h data during the wheat growth period, Hvegm 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.

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. 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: 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: where h 0 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.

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. Remote Sens. 2022, 14, x FOR PEER REVIEW 6 of 19 In addition, the wheat height is retrieved during wheat growth periods by the multifrequency 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].

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 In addition, the wheat height is retrieved during wheat growth periods by the multifrequency 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].

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 Remote Sens. 2022, 14, 4955 7 of 18 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: 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: 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 1 and Ratio 2 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.

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: 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.

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): where F denotes the number of frequencies for a single system. Then, the daily mean rootmean-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): 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: Based on these, the final fusion retrieval values are defined as:

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. where F denotes the number of frequencies for a single system. Then, the daily mean rootmean-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): where HGDD indicates the GDD model simulation and Hsys is the retrieval value of a single system. The weights can be expressed as follows: Based on these, the final fusion retrieval values are defined as:

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. Moreover, the low elevation angles (less than 5 • ) were excluded to avoid obstacle effects surrounding the wheat farmland.  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. 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. Moreover, the low elevation angles (less than 5°) were excluded to avoid obstacle effects surrounding the wheat farmland.

Satellite System
Frequency Band SNR Type 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 Figures 7-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.

Single Frequency
Remote Sens. 2022, 14, x FOR PEER REVIEW 10 of 19 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.   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.  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.  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

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 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. Figure 11. Wheat height retrieval with the GPS system using segmented processing.
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. 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 Figure 11. Wheat height retrieval with the GPS system using segmented processing.
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. 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.

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. 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.

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.

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 Figures 13-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 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.

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 Figures 13-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. 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.   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 shows the accuracy in terms of HNFQ site error statistics based on the multiple retrieval methods.   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 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. 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. Remote Sens. 2022, 14, x FOR PEER REVIEW 16 of 19 Figure 17. Correlation between the proposed method and wheat height.
In addition, the performance of the proposed method is compared with that of other single frequency methods of four GNSS. Table 3 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.

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.
Author Contributions: M.S. proposed the method and designed the experiments; K.C. collected and carried out initial processing of the experiment data; M.S. and K.C. wrote the paper; F.S. modified the paper and provided supervision; all of the authors contributed the discussion. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: All data included in this study are available upon request by contact with the corresponding author.