Measurement of the Sea Surface Height with Airborne GNSS Reﬂectometry and Delay Bias Calibration

: Sea surface height can be measured with the delay between reﬂected and direct global navigation satellite system (GNSS) signals. The arrival time of a feature point, such as the waveform peak, the peak of the derivative waveform, and the fraction of the peak waveform is not the true arrival time of the specular signal; there is a bias between them. This paper aims to analyze and calibrate the bias to improve the accuracy of sea surface height measured by using the reﬂected signals of GPS CA, Galileo E1b and BeiDou B1I. First, the inﬂuencing factors of the delay bias, including the elevation angle, receiver height, wind speed, pseudorandom noise (PRN) code of GPS CA, Galileo E1b and BeiDou B1I, and the down-looking antenna pattern are explored based on the Z-V model. The results show that (1) with increasing elevation angle, receiver height, and wind speed, the delay bias tends to decrease; (2) the impact of the PRN code is uncoupled from the elevation angle, receiver height, and wind speed, so the delay biases of Galileo E1b and BeiDou B1I can be derived from that of GPS CA by multiplication by the constants 0.32 and 0.54, respectively; and (3) the inﬂuence of the down-looking antenna pattern on the delay bias is lower than 1 m, which is less than that of other factors; hence, the effect of the down-looking antenna pattern is ignored in this paper. Second, an analytical model and a neural network are proposed based on the assumption that the inﬂuence of all factors on the delay bias are uncoupled and coupled, respectively, to calibrate the delay bias. The results of the simulation and experiment show that compared to the meter-level bias before the calibration, the calibrated bias decreases the decimeter level. Based on the fact that the specular points of several satellites are visible to the down-looking antenna, the multi-observation method is proposed to calibrate the bias for the case of unknown wind speed, and the same calibration results can be obtained when the proper combination of satellites is selected.


Introduction
Global navigation satellite systems (GNSSs) not only can provide positioning, velocity, and timing for a user but also can be used as sources for remote sensing to explore Earth's physical parameters. This bi-/multistatic observation is known as GNSS reflectometry and was proposed by Marin-Neira in 1993 to provide a high density of sea surface heights through the simultaneous tracking of several observation points [1]. At present, GNSS reflectometry is demonstrated to be capable of measuring sea surface height [2][3][4], wind speed [5][6][7], sea ice [8][9][10], and soil moisture [11][12][13]. Compared to radar altimeters and scatterometers, GNSS reflectometry uses low power and is lower cost because of the lack of a transmitter and can provide higher spatial-temporal resolution with many GNSS satellites. Detailed illustrations of GNSS reflectometry and its application can be found in [14,15]. The success of the UK TDS-1 [16], CYGNSS [17], and Bufeng [18] missions has indicated the capacity to observe global physical parameters using GNSS reflectometry. However, spaceborne GNSS reflectometry is incapable of exploring nearshore areas because of the

Geometry Model
The deterministic specular delay estimated from the waveform can be expressed as (1) where ∆ ion and ∆ trop are biases caused by the ionosphere and troposphere, respectively. For the airborne case, due to the reflected and direct signal, it goes through the same ionosphere path so that the ionosphere biases of the direct and reflected signals are the same, ∆ ion is almost 0 and can be ignored. ∆ trop is considered by using a simple tropospheric delay model as in [30]. (2) where is the elevation angle of the GNSS satellite, h r is the receiver height referring to the reference surface in WGS84, and h trop is the height of the troposphere in the experimental region. ∆ ins is the instrumental delay introduced by the radiofrequency (RF) chains for direct and reflected signals. This system bias can be calibrated by switching the direct and reflected signals to two RF chains in turn in a given period [32] or calibrated with a correlation obtained by precalibration for RF chains.
∆ geo is the propagation delay of the direct and reflected signal, as shown in Figure 1, which is relative to the receiver height and the elevation angle of the GNSS satellite as where is the elevation angle of the GNSS satellite and is accurately computed from the direct signal or IGS ephemeris, H r is the receiver height referring to the sea surface, and d is the distance between the up-and down-looking antennas. Once H r is computed from Equation (3), the sea surface height can be computed as b spec is the bias between the estimated delay by retracking the feature point and the specular delay and is called the specular delay bias. Precisely calibrating the above bias is important for improving the precision of the retrieved sea surface height.

Delay Estimation
The receiver height referring to the sea surface can be obtained using the geometry of GNSS reflectometry and is important for estimating the delay between reflected and direct signals from the measured waveform. At present, two schemes, including retracking given points [27][28][29] and fitting models [30], have been utilized. The method of retracking a given point is used in this paper because of the low complexity of this technique. A retracked point in the reflected waveform that is easy to extract from the waveform and is rarely influenced by the noise and sea state should be defined. As presented in Figure 2, which shows the simulated delay waveforms used the Z-V model [31] for different wind speed, compared to the waveform peak and trailing edge of the waveform, the leading edge is relatively insensitive to the sea state. The peak of the derivative waveform is chosen as the retracked point to estimate the arrival time of the reflected signal and corresponding delay is defined as where W(τ) is the delay waveform of the reflected signal. The peak of the derivative waveform is used in this paper to estimate the delay of reflected signals. The sparse sampling waveform is first interpolated to produce the dense-sampling waveform, and then the interpolated waveform is differentiated and searched for the peak of the derivative waveform. The alternative approach is to first use a polynomial to fit the leading edge of the measured waveform and then analytically compute the zero-crossing point of the second derivative of the fitted polynomial [33]. In this paper, the first method is used to retrack the peak delay of the derivative waveform.

Specular Delay Bias
For an idealized situation in which the impulse response starts to spread from the specular delay to a longer delay, the specular delay is theoretically the same as the delay of the peak of the derivative waveform [27]; however, the actual scenario is complete, so the specular delay has a bias from the delay estimated using the peak of the derivative waveform; this bias should be evaluated and calibrated. In [29], the delay derived using the peak of the derivative waveform is calibrated with a correlation computed from the modeled waveform. To further analyze the specular delay bias, the waveform model of reflected GNSS signals is used [31].
where P t and G t are the transmitted signal power and transmitting antenna gain, respectively; G r is the antenna gain of the receiver; λ is the wavelength of the GNSS signal; R t and R r are the distances from the receiver and the GNSS satellite to the observation area, respectively; Λ(·) is the autocorrelation function of the pseudorandom noise (PRN) code; and σ 2 pq is the scattering coefficient from p polarization to q polarization, which could be computed using Equation (34) in [31].

Elevation Angle
The reflected GNSS signals depend on the surface roughness. The criterion usually used to distinguish smooth and rough surfaces is the Rayleigh criterion, as in [34] For a satellite with a lower elevation angle under the same sea state conditions, the surface is smoother, so the delay waveform is more similar to that of the direct signal and the arrival time of the specular signal is closer to the peak of the delay waveform. Figure 3a shows the changing delay bias versus the elevation angle for GPS CA, Galileo E1b, and BeiDou B1I when the height and wind speed are 3500 m and 5 m/s, respectively, in which the peak delay of the derivative waveform has a negative bias relative to the specular delay and increasing the elevation angle makes the bias lower and ultimately converge to a stable value.

Height
For the diffuse reflection, the reflection zone is jointly determined by the first isorange and glistening zone. The size of the first isorange is computed by 2πab, where a and b are the semimajor and semiminor axes of the first isorange, which are estimated through Equation (3) in [35]. The size of the glistening zone is defined as the equivalent area of the scattering coefficient by where σ pq,max is the peak of the scattering coefficient and computed as max{σ 2 pq } (where max{·} is the maximum operator); r represents the scattering area. As shown in Figure 4, as the receiver height rises, the area of the first iso-range shows the approximately linear growth; however, the equivalent area of the scattering distribution presents the exponential growth. Therefore, the ratio between the equivalent area of the scattering distribution and the area of the first isorange shows an increasing trend. For a lower receiver height, the glistening zone is smaller than the ellipse of the first iso-range and the delay waveform is similar to the ideal triangle function. In this case, the optimal feature point to estimate the specular delay is the waveform peak rather than the peak of the derivative waveform. This illustrates the phenomenon shown in Figure 3b which shows that the lower the receiver height is, the larger the specular delay bias is.   Figure 3c gives the specular delay bias with respect to the wind speed. The absolute bias presents a decreasing trend as the wind speed increases and becomes stable after 12 m/s. For the case of a lower wind speed, the scattered signals come from the zone with the shorter delay, and the delay waveform is more similar to the triangle function than the peak of the derivative waveform; therefore, the delay estimated by the peak of the derivative waveform has more severe bias. Note that for the relatively calm sea surface, specular reflection occurs and coherency properties of specular scattering can be observed [36], and the coherent components coming from the Fresnel zone dominate the reflected signals. In this case, the delay waveform is identical to that of the direct signal, so the unbiased estimation of specular delay can be obtained by retracking the waveform peak.

Pseudorandom Noise Code
The PRN code has an important influence on the estimated delay. A code with a wider bandwidth produces sharper edges, which enables more accurate and precise estimation. The Cramer-Rao bound (CRB) of the ranging precision is inversely proportional to the square of the signal bandwidth [37]. Equation (6) can be written as the convolution between Λ 2 (τ) and the bistatic radar cross-section including the relative position between the transmitter and the receiver, the sea surface parameters, and the antenna footprint. Intuitively, the impact of the PRN code on delay bias is uncoupled from the elevation angle, receiver height, and wind speed. The Z-V model is used to simulate and demonstrate this intuition. Figure 5 shows the simulated biases which are normalized by the bias of GPS CA versus the elevation angle, receiver height, and wind speed. The normalized bias of Galileo E1b presents weak fluctuations for the case of receiver heights lower than 1500 m or wind speeds less than 3 m/s. However, under the other conditions, the biases of Galileo E1b and BeiDou B1I are stable with respect to that of GPS CA. This illustrates that the bias of Galileo E1b and BeiDou B1I can be derived from that of GPS CA by multiplying by the constants 0.32 and 0.54 for Galileo E1b and BeiDou B1I, respectively. The above constants are computed through averaging the normalized biases of the Galileo E1 and BeiDou B1I in Figure 5.

LHCP Antenna
As shown in Equation (6), the antenna gain pattern modulates the signal from different directions. This modulation generates different deformations in the delay waveform from the different satellites. In this paper, only the case of the down-looking and circular-beam antenna with a gain of 12 dB and a beam width of 30 • is analyzed. Figure 6 gives the delay bias considering and ignoring the impact of the antenna gain pattern. The antenna pattern has a more serious influence on the delay bias in the case of high wind speed and low elevation angle. The bias difference between the considered and ignored impact of the antenna gain pattern is under 1 m. Due to the impact of the antenna gain pattern being lower than that of the elevation angel, receiver height, wind speed and PRN code, in the development of the calibration model of the delay bias in this paper, the antenna influence is ignored.

Calibration of the Specular Delay Bias
When the influence of the antenna is ignored, the specular delay bias is jointly determined by the elevation angle, receiver height, wind speed, and PRN code by where S ε is the sine of the elevation angle, v ws is the wind speed, R c is the parameter characterizing the impact of the PRN code, and F(·) represents the mapping relationship between the delay bias and its influencing factors.

Analytical Model
To develop the analytical calibration model, it is assumed that the influences of the elevation angle, receiver height, wind speed, and PRN code on delay bias are independent. It is known that the delay biases of Galileo E1b and BeiDou B1I are stable with respect to that of GPS CA and can be obtained from the delay bias of GPS CA through multiplication with constants. The analytical model is expressed as R c of GPS CA, Galileo E1b, and BeiDou B1I are 1, 0.32 and 0.54, respectively. For the given reference of elevation angle, receiver height, and wind speed, the delay bias of GPS CA is the sines of the referenced elevation angle, receiver height, and wind speed, respectively, which are sin(80 • ), 1000 m, and 2 m/s in this paper. The delay bias of the reference case is derived through simulation using the theoretical model given by Equation (6). Here, f " (S ε ), g h (h r ) and h ws (v ws ) are respectively defined as Equations (12)- (14) are substituted into Equation (10), and the delay bias is expressed as Further work is needed to obtain the expressions of f " (S ε ), g h (h r ) and h ws (v ws ). For the referencing elevation angle, the delay bias is Similarly, for the given receiver height and wind speed, the delay biases are Equation (16) is divided by Equation (10), and f " (S ε ) is expressed as Similarly, g h (h r ) and h ws (v ws ) are respectively expressed as The Monte Carlo method is used to obtain f " (S ε ), g h (h r ) and h ws (v ws ) through the following steps: 1.
the elevation angle, receiver height, and wind speed in the above group are replaced with the given references to produce the other three groups, recorded as the above four groups of parameters are used as the input to simulate the delay waveforms of GPS CA, Galileo E1b, and BeiDou B1I through Equation (6); 4. the delay bias is estimated and recorded as b speci G , b spec_εi G , b spec_hi G and where · is the mean operator and ε s , h s and v s are the interval steps of the elevation angle, receiver height, and wind speed, respectively, which are 2 • , 200 m and 0.2 m/s. Figure 7 gives the results of Monte Carlo for GPS CA, Galileo E1b, and BeiDou B1I, in which f ε (S ε ), g h (h r ) and h ws (v ws ) are inversely proportional to the elevation angle, receiver height, and wind speed and can be best matched by the rational function of In addition, f ε (S ε ), g h (h r ) and h ws (v ws ) of GPS CA, Galileo E1b, and BeiDou B1I have similar changing tendencies and can share the same fitted parameters in Equation (25). The fitted parameters used in this paper are given in Table 1.

Neural Network
Artificial intelligence (AI) methods, such as neural networks [38], can be used to develop predictive models based on empirical data. In this paper, a back propagation (BP) network, as shown in Figure 9, is trained to develop the calibration model of the delay bias. The inputs of neural network are the sine of the elevation angle S ε , height h r , wind speed v ws , and R c . The output is the predicted delay biases. Note that R c of the neural network is quantified as 0, 1, and 2 for GPS CA, Galileo E1b, and BeiDou B1I signal. The data generated in Section 4 are randomly divided into two sets. One set which accounts for 70% of the total data used to train the neural network and the other set accounting for 30% of the total data is used to test the developed neural network. The ability to approximate the analytical function is dependent on the number of hidden neurons and activation functions [39]. The usually used activation functions are purelin, tanh, and sigmoid function. Table 2 gives the RMSEs of the biases in test set for the different combinations of the activation function in hidden and output layer when the hidden node number is 6. From the table, it is seen that when the activation functions in hidden and output layer respectively are tanh and purelin function, the RMSEs are relatively low. Therefore, in this paper the following tanh and purelin functions are used as the activation functions in the hidden and output layers of the neural network.  Figure 9. Structure of the adopted neural network. From Figure 10, it is clear that the RMSEs of the biases in test set for GPS CA, Galileo E1b, and BeiDou B1I are all decreasing as the hidden node number increases and becomes stable after the hidden node number is 6. Therefore, the neuron number of the hidden layer in this paper is 6. Figure 11 gives the simulated delay biases and delay biases computed using the developed neural network model, in which the RMSEs of GPS CA, Galileo E1b, and BeiDou B1 are 0.38, 0.26, and 0.16, respectively, which are slightly better than those of the above analytical model.

Multisatellite Observation
The analytical model and neural network above require information about the elevation angle, receiver height, and wind speed. The elevation angle and receiver height can be obtained from the direct signal. In addition to the existing methods for obtaining wind speed such as scatterometers and nearshore buoys, reflected GNSS signals provide wind speed by best matching the measured delay waveform to the theoretical delay waveform for airborne scenarios [40]. Whether the delay bias can be calibrated for the case of unknown wind speed should be explored. It is known that the specular points of several satellites are visible for a down-looking antenna. Assuming that the wind speed is identical in the coverage of the down-looking antenna, h ws (v ws ) in Equation (16) is the same for different satellites. The Equation (16) is expressed as where b comm = g h (h r )h ws (v ws )b spec_ref . When the down-looking antenna covers N specular points, calibration models of N satellites can be expressed as where A, X and Y are respectively expressed as X in Equation (29) could be solved through least square as Assuming that the measurements of different satellites are independent and that the errors are subject to the normal distribution, the mean and variance of the ith measurement error are e y i = 0 and V e y i = σ 2 y , respectively. The covariance matrix K y of Y is Iσ 2 y . The covariance matrix of X is derived as where e y is the error vector of the measurement Y. From the above equation, it is seen that the height variance σ 2 h measured using multi-observation is not only dependent on the variance of the estimated delay but also influenced by a weight coefficient as (35) where Ξ i = R ci f ε (ε i ). From the above equation, it is seen that the variance of the height is dependent on the elevation angle and PRN code. When the delay waveforms used for the multi-observation have the same PRN type and elevation angle, the variance of the height tends to infinity; in other words, A is under rank, and Equation (32) is an underdetermined equation that has innumerable solutions. To obtain the proper height, two conditions are required:

1.
The elevation angles of the chosen satellites should be as different as possible; 2.
The frequency spectra of the delay waveform should be as different as possible.

Validation
Simulations and experiments are performed to demonstrate the validity of the proposed calibration models.

Simulation
The complex correlation value of the reflected signal in the presence of noise can be modeled as [41] Y si (τ) = P(τ) 2 Z si + √ P n · Z ni (36) where Z si and Z ni are complex Gaussian white noise with real and imaginary parts; P n is the noise power, which is expressed as [42] where T c is the coherent integration time, k is the Boltzmann constant, T R is the receiver noise equivalent temperature, and B w is the receiver bandwidth. To obtain waveforms with different SNRs, the noise power is multiplied by a factor η as P n = ηT 2 c kT R B w . The above complex correlation waveforms are incoherently averaged as (38) where N is the incoherent averaging number. The simulation scenario is shown in Figure 12 and Table 3 Tables 4 and 5 give the bias and RMSE of the retrieved height before and after the calibration. The calibration is unable to improve the RMSEs of the retrieved sea surface height, and the RMSEs after calibration and without calibration are both at the decimeter level. In addition, Galileo E1b and BeiDou B1I yield lower RMSEs than GPS CA because of the wider frequency spectrum. The results of Table 4 show that the uncalibrated biases are −5.41, −1.68, and −2.77 m, and after calibration using the proposed models, the biases decrease to the decimeter level. This illustrates that the proposed models can significantly calibrate the delay bias.  To illustrate the influence of matrix A on the RMSE of the retrieved height, GPS PRN 07, 17, and 18 given in Table 6 are chosen to retrieve the sea surface height using multiobservations. The elevation angles of the three satellites are close, and the R c values are all 1; therefore, the matrix A is perishing in the solution of Equation (33). In this case, the bias and RMSE of the retrieved height are −0.63 and 11.10 m, respectively, which are worse than the RMSE before the calibration. The corresponding Π is 102.74, which is much larger than the 0.96 of the combination of GPS PRN 18, Galileo PRN 06, and BeiDou PRN 10. Comparing the results of two multisatellite observations demonstrates that the geometric configuration is important for achieving proper calibration and preventing the RMSE of the retrieved height from increasing. A suitable configuration should have as different of elevation angles and frequency-spectrum waveforms as possible.

Experiment
The dataset used in this paper was collected during an airborne experiment over the Baltic Sea on 3 December 2015. Two 8-element antennas were installed at the aircraft. One is the up-looking antenna, and the others is the down-looking antenna. RF chains with a sampling rate of 80 MHz, IF of 19.42 MHz, and quantization of 1 bit were used to collect the direct and reflected GNSS signals. The L1 signal of the GPS and E1b signal of Galileo were collected between seconds of the day (SoD) 39,217 and 39,617. During this period, the aircraft flew at an altitude of ∼3 km and a velocity of ∼50 m/s. As shown in Figure 13, PRN01 and 11 of GPS and PRN11 of the Galileo System were visible for both the up-looking and down-looking antennas. The elevation angles of these three satellites are computed using IGS ephemeris and present ranges of 72.33 • ∼72.17 • , 55.71 • ∼53.01 • and 55.90 • ∼55.14 • , respectively. The amplitude of the tides is just a few centimetres in the Baltic Sea, and the salinity, temperature, precipitation, and evaporation also have minor effect on the Baltic sea level [43]; the sea surface height from the DTU10 model [44] is considered the ground truth for comparison with the retrieved sea surface height using the reflected GPS and Galileo signals.
A software receiver, as given in Figure 14, is developed to obtain the delay waveform of the reflected GPS and Galileo signals. The software receiver consists of direct and reflected channels for processing direct and reflected GPS and Galileo signals. The direct and reflected channel share the local replicas from the local signal generator in Figure 14. The direct signals are acquired and tracked to provide a reference for the processing of the reflected signals. The code parallel algorithm in frequency domain is used to accelerate the acquiration procedure of the direct signal. The Phase Locked Loop (PLL) and Delay Locked Loop (DLL) are respectively adopted to dynamically track the carrier and code phase of the direct signal. The detail processing of the direct signal could be seen in [45]. The reflected channel is composed of a set of correlators that perform cross-correlation between the reflected signals and the local replica at different delays. To improve the processing efficiency of the reflected signals, cross-correlation in the frequency domain using the fast Fourier transform (FFT) is implemented. The coherent integration time in the software receiver is the same as the period of the PRN code, i.e., 1 ms and 4 ms for the GPS and Galileo signals, respectively. After cross-correlation, a series of complex correlation waveforms are incoherently averaged to improve the SNR of the waveform. To ensure an output rate of 1 Hz, the incoherent averaging numbers for the waveforms of GPS and Galileo are 1000 and 250, respectively. The ephemeris is used to compute the elevation angles of the satellites. The following steps are conducted to retrieve the sea surface height: 1.
The receiver position and elevation angle of the GNSS satellite are computed from the direct signal and IGS ephemeris; 2.
The peak of the derivative waveform is retracked to estimate the delay from the incoherently averaged waveforms; 3. ∆ trop and b spec are calibrated to obtain ∆ geo ; 4.
Equation (3) is used to compute the receiver height referring to the sea surface; 5.
The sea surface height is retrieved using Equation (4), and a moving average is obtained for the retrieved heights; 6.
A comparison is made with the DTU10 data to obtain the bias and RMSE of the retrieved sea surface height. Figure 15 shows the in situ and retrieved sea surface heights, and both present the same trend. Tables 7 and 8 give the bias and RMSE of the retrieved height, respectively. Similar to the results of the simulation, decimeter-level RMSE can be obtained both before and after the calibration, and the RMSE of Galileo E1b is lower than that of GPS CA. In addition, PRN 01 of GPS yields a better RMSE than PRN 11. This occurs because the larger the sine of the elevation angle is, the higher the retrieved precision is [46]. Meter-level bias is present before the calibration, and decimeter-level bias is achieved after the calibration.
To validate the influence of the selection of satellites on the height using multisatellite observations, the heights of different configurations are obtained. Table 9 shows that the GPS PRNs of 01 and 11 yield different elevation angles; however, Π enlarges, so the bias and RMSE of the retrieved height are unacceptable. When GPS PRN 11 and Galileo PRN 11, which have similar elevation angles, are chosen, proper results are achieved. This is mainly because, as shown in Figure 3, when the elevation angle is over 50 • , the elevation angle has an insignificant effect.

Discussion
This paper proposes three methods for calibrating the delay bias; however, after the calibration, the bias of the retrieved height is still at the decimeter level. Theoretically, fitting a model to the measured waveform may produce unbiased estimation; however, note that the matching procedure is computationally complex and that its retrieved performance depends on the theoretical model used. Therefore, to fit a model, it is important to develop a precise waveform model. The model currently usually used (the Z-V model) was proposed by V. U. Zavorotny et al. [31] based on the Kirchhoff approximation and geometrical optics (KA-GO); this model describes the power distribution of the scattered GNSS signal and has been used to approximate the measured waveform. However, the Z-V model ignores the coherent component in the reflected signal. For some cases, such as calm seas and lakes, strong coherent scattering in the near-forward bistatic geometry takes place. In these scenarios, the Z-V model produces an incorrect result for the fitting of the measured waveform.
The comparison between the retrieved results of the GPS L1 CA and Galileo E1b code clearly shows that the E1b code yields better performance than the CA code because of the wider signal bandwidth. iGNSS-R was proposed because it can use the full-bandwidth signal to obtain higher precision in altimetry. At present, to use wide-bandwidth signals, especially encrypted signals, such as the P(Y) code of GPS, some approaches have been proposed, including partial interferometer [47] and semicodeless technology [48]. In addition, due to the wider frequency spectrum, some other signals from communication satellites have been used to provide more precise altimetry [49]. The proposed methods can also be used to calibrate the bias of the above signals; however, the fitting parameters of the analytical models have to be obtained, and the neural network should be trained again.

Conclusions
This paper analyzed the bias of the delay estimated by retracking the peak of the derivative delay waveform. The results showed that the elevation angle, receiver height, wind speed, PRN code, and down-looking antenna pattern influenced the delay bias. The impact of the PRN code was uncoupled from the elevation angle, receiver height, and wind speed, so the delay biases of Galileo E1b and BeiDou B1I could be obtained from the bias of GPS CA by multiplication with the constants 0.32 and 0.54, respectively. From the assumptions that the effects of the elevation angle, receiver height, and wind speed on the delay bias were uncoupled and coupled, an analytical model and a neural network, respectively, were proposed to calibrate the delay bias. Simulations and experiments were conducted to validate the proposed models. The results showed that through calibration using the proposed models, the biases of the retrieved height were reduced to the decimeter level from the meter level. Because the reflected GNSS signals from several satellites could be received, multisatellite observations were used to calibrate the delay bias under the condition of unknown wind speed. The results of the simulation and experiment illustrated that the calibration performance depended on the chosen satellites, which had different elevations and PRN codes. However, in this paper, after calibration, the decimeter-level bias and RMSE of the retrieved height were achieved, and it is necessary to pursue a measured height with a lower bias and RMSE. Therefore, in the future, work exploring the measurement of the height using wider-spectrum signals, such as L5 of GPS and B3I of BeiDou, will be conducted to obtain better altimetry performance.