Analysis of GNSS-R Code-Level Altimetry Using QZSS C/A, L1C, and BDS B1C Signals and Their Combinations in a Coastal Experiment

The fully completed Quasi-Zenith Satellite System (QZSS) and BeiDou Navigation Satellite System (BDS-3) can provide higher precision code observations owing to their advanced signal modulation. However, the performance of Global Navigation Satellite System-Reflectometry (GNSS-R) code-level altimetry using the new codes is rarely be tested. In this study, interoperable signals with the frequency of 1575.42 MHz, i.e., L1C/A and L1C from QZSS and B1C from BDS, were tested for GNSS-R altimetry simultaneously for the first time. Moreover, a joint weighted GNSS-R altimetry method is proposed to improve retrieval precision based on combinations of the path delays obtained from different signals. To investigate the performance of the signals and their combinations, a ground-based GNSS-R altimetry experiment was conducted and raw intermediate-frequency data were collected over several hours. These data were processed using a software-defined receiver to compute the path delay between the direct and reflected signals. Solutions were derived separately for each code and their combinations. Through comparison with reflector heights obtained from an in situ 26 GHz radar gauge, results showed that the performance of GNSS-R altimetry based on binary offset carrier modulation signals was better than that from binary phase-shift keying. Moreover, results from the joint weighted GNSS-R altimetry method showed further improvement in altimetry precision.


I. INTRODUCTION
C HANGE in global mean sea level is important to the natural environment and to socioeconomic and ecological systems in coastal areas [1]. Traditionally, global sea surface height is measured using buoys, satellite radar altimeters, and tide gauges [2], [3], [4]. However, such observational methods cannot meet the requirements of mesoscale sea surface monitoring with high spatiotemporal spatial resolution [5], [6], especially in coastal areas. Global Navigation Satellite System-Reflectometry (GNSS-R) is a new approach to ocean altimetry with unique advantages, such as all-weather capability, low power consumption, and high-density measurements. At present, there are two usual kinds of GNSS-R altimetry methods: one is GNSS-R altimetry based on SNR observations from geodetic receivers. For most GNSS stations, optimal centimeter-level sea surface height can be achieved, such as in [7], [8], and [9]. However, for this kind of GNSS-R altimetry, one sea level estimation is retrieved from about 20-min to 1-h continuous observations with low elevations, which limits the daily quantity of the available estimations and, thus, the low temporal resolution of the sea level retrievals. Moreover, one sea level estimation is the average over the whole period of the low elevation angles, which results in the low spatial resolution of the retrievals [10]; the other kind is specialized dual-antenna GNSS-R setups based on path delay between the direct and reflected signals. This kind of GNSS-R altimetry uses several millisecond delay waveforms to calculate the path delay and further the sea surface height. Therefore, it can afford sea level estimations with high temporal and spatial resolution.
In GNSS-R altimetry, stable high-precision determination of the path delay between direct and reflected received GNSS signals is critical [11]. The path delay can be obtained directly from the two major types of GNSS observables, i.e., the code and carrier phases. Carrier-phase-delay altimetry is more precise than code-delay altimetry because of the much shorter wavelength used [12]. However, when the reflective surface is rough, carrier-phase-delay altimetry becomes unstable because the carrier phase cannot be tracked continuously [13]. The code delay is used widely in GNSS-R altimetry because of the reasonable autocorrelation and cross-correlation characteristics of the pseudocode itself, which makes it unnecessary to continuously track the code phase when calculating the path delay of the direct reflection signal [14]. In comparison with carrier-phase-delay altimetry, code-delay altimetry has higher stability but relatively poorer precision owing to the long chip length (low chip rate) of the openly available code [15]. However, there is the literature using the method of sliding window of 5 min to fulfill the This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. measurement of sea surface height with RMSE of 5.3 cm by BeiDou Navigation Satellite System (BDS) BIC code [10]. Therefore, in general, GNSS-R code altimetry can be expected to achieve centimeter-level accuracy measurement. Therefore, to investigate the precision of code-delay altimetry, this study focused on GNSS-R altimetry using code observations.
In 1993, the PARIS concept of a passive reflectometry and interferometry system and the use of global positioning system (GPS) reflectometry for sea level monitoring were proposed by Martin-Niera [15]. In 1994, Auber et al. [16] received GPSreflected signals from the ocean during an airborne experiment. To verify the feasibility of GNSS-R altimetry, a shore-based dual-antenna experiment using the GPS L1C/A signal was performed in 1997 in [17] at Rotterdam (The Netherlands). The results showed that accuracy was approximately 3.3 m, which was approximately 1% of the L1C/A code chip length. In 2000, an airborne GNSS-R code-level altimetry experiment was performed using reflected GPS C/A signals and sea surface height measurements with meter-level accuracy were achieved [18]. An airborne GNSS-R code-level altimetry experiment using authorized P(Y) codes was also performed and solutions with centimeter-level accuracy were obtained [19]. In 2014, sea surface topography in the Mediterranean Sea near Italy was retrieved using airborne GNSS-R carrier-phase data, and the experimental results revealed that conditions for such retrievals were optimal at elevations of 10°-30° [20]. To overcome constraints imposed by the signal structure and the nature of a GNSS, PARIS was improved by cross correlating the reflected signals with direct GNSS signals instead of using clean replicas of openly available code [21]. A ground-based experimental demonstration of the new concept was conducted in [22]. The use of BDS GEO satellite in a coastal configuration was first demonstrated in [23] for both wind speed and ocean altimetry. In 2015, Zhang et al. [24] performed a code-level altimetry ground-based experiment using signals from geostationary earth orbit and inclined geosynchronous orbit (IGSO) satellites of the BDS. In 2018, a shipborne GNSS-R code-level experiment based on GPS and BDS signals was conducted and submeterlevel accuracy was realized [10]. In 2020, results of a coastal experiment using BDS-3 B1C and B2a signals showed that GNSS-R code-level altimetry could achieve centimeter-level precision after application of the moving average [25]. The first spaceborne observation of sea surface height through GPS reflectometry with meter-level precision was conducted using the TechDemoSat-1 platform in 2016 [26].
Ground-based GNSS-R altimetry technology can be used as an effective supplement to existing tide gauge stations and other methods used for sea level monitoring in coastal areas. In the construction of traditional tide gauge stations, it is necessary to build tide gauge wells and other infrastructure that require complex construction processes, considerable effort, and substantial material resources. If GNSS-R technology were used for such observations, the only equipment needed could be installed at a suitably high location on the coast, which would be more convenient and incur comparatively low cost. Dense deployment of GNSS-R monitoring equipment in coastal areas could provide important basic data for marine science and geodesy research.
In recent years, although many experiments have been conducted using different platforms to explore the performance of GNSS-R code-level altimetry, most have used GPS L1C/A ranging codes. The modulation of the L1C/A code is binary phase-shift keying and its bandwidth is narrow, which restricts the code precision. As the development of spaceborne GNSS-R missions, such as the CYGNSS, there are also altimetric results based on other GNSS signals. The ocean altimetry performance of spaceborne GNSS-R by processing the reflected waveforms of GPS L1, Galileo E1, and BDS-3 B1 signals from CYGNSS mission was assessed in [27], and meter-level altimetry precision was achieved. First, carrier-phase sea surface altimetry using CYGNSS data transmitted from both GPS and Galileo constellations, and GNSS-R altimetric precision of 3-4 cm at 20 Hz sampling, centimeter-level precision at 1 Hz were realized by comparing with dedicated radar gauge data [28]. The bistatic group delay and carrier-phase delay observations are extracted from 12 tracks of raw CYGNSS data, and the derived surface topography profiles from GPS and Galileo observations show good self-consistence along the coincident tracks [29]. Signals with different modulation and bandwidth will show different performance in GNSS-R code-level altimetry, and their capacity for sea level estimation should further be assessed. Moreover, the lower satellite elevation amplifies the noise in the relative path delay, which ultimately affects the accuracy of ocean altimetry and should be considered to be weaken during the data processing.
In 2020, the full deployment of the BDS-3 constellation was completed [30], which can transmit the B1C signal by medium earth orbit and IGSO satellites, while the L1C signal can be transmitted by GPS and the Quasi-Zenith Satellite System (QZSS) [31], [32], [33]. The B1C and L1C signals based on binary offset carrier modulation can improve the accuracy of GNSS-R code-level altimetry. Additionally, with the rapid development of GNSSs, it is expected that there will be more than 150 satellites in orbit in the future continuously transmitting L-band signals to the ground. Thus, it could be expected that the multiple GNSSs will improve the accuracy and reliability of ocean altimetry results substantially during the same observation period.
Currently, GPS broadcasts only the L1C signal on Block III satellites. The number of satellites in this series is small, and the signal has poor coverage over coastal areas of eastern China. As an enhancement of GPS in the Asia-Pacific region, QZSS modulates the L1C code and the traditional L1C/A code on the L1 frequency band, which means that coastal areas of eastern China can receive the reflected signal from QZSS. It provides the possibility of investigating the use of the L1C code and the L1C/A-code for GNSS-R high-accuracy sea surface retrievals. In this study, the B1C signal from BDS, and the L1C/A and L1C signals from QZSS were selected to test GNSS-R code-level altimetry. Then, a combination method for multiple satellites and multiple signals, based on a weighted least squares iterative algorithm, was proposed to improve retrieval precision.
The rest of this article is organized as follows. In Section II, the principles of the method for GNSS-R code-level altimetry using GNSS signals are introduced. Then, the proposed combination method for multiple satellites and multiple signals is described. In Section III, the details of the experiment are provided. In Section IV, GNSS-R data processing and analysis of the results obtained using the proposed method are presented. Finally, a discussion of the results and a summary of the findings are presented in Sections V and VI, respectively.

A. Ground-Based GNSS-R Code-Level Altimetry
The geometric relationship between a GNSS satellite, receiving antenna, and specular reflection point is shown in Fig. 1. It can be seen that sea surface height can be deduced by subtracting the vertical height (h) between the receiving antenna and the sea level, i.e., the reflector height (RH), from the geodetic height of the antenna. Because the geodetic height of the antenna can be deduced from traditional positioning methods such as PPP, the key to GNSS-R altimetry lies in precise determination of RH. According to the geometry in Fig. 1, RH can be calculated from the path delay between the reflected and direct signals using (1)-(3); therefore, the core of the GNSS-R code-delay sea level measurement technique is precise calculation of the path delay.
In Fig. 1, S is the specular reflection point, h is the vertical distance from the phase center of the reflective antenna to the sea surface, θ is the elevation angle of the satellite, and Δρ is the path delay of the reflected signal relative to the direct signal, which can be expressed as follows [17]: where ρ r and ρ d denote the propagation path lengths of the reflected and direct signals, respectively, and ε ρ is the index of the errors related to the signal propagation path. According to [34], conventional GNSS-R is the most widely used reflectometry technique. Therefore, in this study, the propagation path delays between the reflected and direct signals were calculated based on conventional GNSS-R using code and our software-defined receiver [25].
The direct RHCP and reflected LHCP digital IF signals from QZSS and BDS were collected. We used a receiver with a bandwidth of 20.46 MHz and a central frequency at 0.42 MHz. Although the bandwidth of the GNSS-R instrument (20.46 MHz) is not sufficient to capture the full bandwidth of the B1C signal (32.736 MHz), these observations from BDS satellites can still provide the opportunities to assess the altimetry performance of the navigation signals with different bandwidths and modulations [27]. Then, the received direct and reflected RF signals are separately amplified with two low-noise amplifiers and filtered using a bandpass filter to reduce noise. The signals then entered the downconverter and mixed with a signal provided by one numerically controlled oscillator to be analog IF signals. After filtered by two bandpass filters, the analog signals are quantized with 4 bits using analog-to-digital converters at a sampling rate of 40 MHz for easy storage. The 40 MHz IF signals are used in our software-defined receiver to calculate the path delays. The principles of code-level path-delay calculation in our software-defined receiver were as follows: our receiver generates the local QZSS L1 C/A, L1C, and BDS B1C code sequences and up-converts them to 40 MHz. Then, these clean replicas of the code are cross correlated with direct and reflected IF signals separately to provide delay waveforms. From the waveforms, the time delay between direct and reflected signals can be determined by finding the arrival time of direct and reflected signals. After 10 ms cross-correlated separately both the direct and reflected signals with the pilot and data local replicas, the two obtained waveforms were added incoherently to increase the signal-to-noise ratio. In this article, the direct and reflected raw samples of all obtained IF signals are separately correlated with 10-ms locally generated replica. Therefore, in 1 s, there will be 100 path delays. We adopt the median value of the 100 path delays as the final path delay to avoid any gross error for following RH calculating. At present, two type of feature point in the delay waveform have been proposed as the arrival time of the reflected signal in coastal GNSS-R altimetry. One is the peak point of the waveform; the other is the peak point of the derivative waveform [35]. In [36], it was found that the peak of the waveform can provide better performance than the peak of derivative waveform. Therefore, we used the peak of the waveform as the arrival of the reflected signal. For better precision, cubic spline interpolations were applied to compute the code-level path delays from the peaks positions of the aforementioned waveform. In addition, the bandwidth of the IF recorder was 20.46 MHz, while the full bandwidth of B1C QMBOC is 32.736 MHz. In this study, in order to achieve higher precision, we used signals of QMBOC modulation for GNSS-R altimetry.
According to the geometric relationship between the path delay and the satellite elevation shown in Fig. 1, after accurate determination of the propagation path lengths of the reflected and direct signals, the RH (h) can be obtained as follows: where a is the baseline vector between the direct and reflective antenna phase centers, e is the vector between the direct antenna phase center and the direction of the satellite connection, V 1 V 2 is the projection in the direction of the satellite and direct antenna phase center connection that can be calculated using (3), and ε h is the RH error.

B. Characteristics of QZSS L1C/A, L1C, and BDS B1C Signals
Recently, with the continuous upgrading of GNSS coding structure and modulation mode, a variety of civil ranging codes with better antimultipath performance and higher precision have appeared. For the frequency of 1575.42 (see Table I), QZSS, which is the augment of GPS, can also transmit L1C civil code signals with wider bandwidths of 30.69 MHz, except for traditional L1C/A code. Same with GPS, QZSS L1C/A code, with the bandwidth of 2.046 MHz and the chip length of 1023, adopts BPSK modulation. However, the L1C signal adopts more advanced modulation mode of BOC for data and TMBOC for pilot. Chinese BDS-3, which finished its construction in 2020, also transmitted B1C signals at the central frequency of 1575.42 MHz. QZSS L1C and BDS B1C have the same chip length of 10 230, while the chip rates are 1.023 MHz and 1.023 Mb/s. Different ranging codes are modulated on the data and pilot component. The power ratio of the data and pilot components is 1:3 for both signals. However, the components of B1C adopt BOC and QMBOC modulation. Different with QZSS, the bandwidth of BDS B1C is 32.736 MHz. Compared to QZSS L1C/A, QZSS L1C and BDS B1C signals have wider bandwidth and more advanced modulation mode which is capable of proving more precise measurement. In this article, we choose these three signals with the same central frequency, different bandwidth, and different modulation modes, to assess their capacity on GNSS-R altimetry.

C. Combination Method for Code-Level Altimetry Using BDS and QZSS Signals
The path delay of single signal, which was computed from our software-defined receiver, can be processed further for GNSS-R altimetry using (2). By using path delay of single signal, the decimeter-level precision of sea level can be estimated. However, sea level with decimeter-level precision cannot fulfill the need of marine research, and it is necessary to find method to obtain sea levels with higher precision. Here, in this section, we proposed a combination method using these path-delay values of different satellites and constellations to realize sea level retrieval with higher precision. This method combines the path delays of different signals to realize sea level retrieval using weighted least squares method.
All the path delays obtained from the software-defined receiver using code observation of multiple satellites can be expressed as follows: where t i is the length of time, and prn j (j = 1 . . . n) is the PRN number of GNSS satellites.O stands for all the path delays output by receiver, which is arranged in chronological order, and O i stands for all the path delays from different signals at the ith epoch. All the path delays in each O i are used for combination. In (4), if all used satellites in O i come from the same GNSS system, the proposed method will be the combination of path-delay values of different satellites from one GNSS system. If n satellites come from different GNSS systems, the proposed method will be the combination of path-delay values of different satellites from multi-GNSS systems. In this study, for comparison, the combinations of signals of different satellites from the single QZSS system, single BDS system, and both systems were analyzed, respectively, in Section IV. Then, according to (2), for single epoch t and for all the combined k satellites, the RH derived from the combination of different GNSS signals using each O i can be expressed as follows according (2): Equation (5) is equivalent to the following matrix: In GNSS-R altimetry, an important bias is electromagnetic bias, which increases with increasing incidence angle. Moreover, analysis of the results of a shore-based experiment based on the BDS/QZSS system conducted in a previous study revealed that the code delay in sea surface altimetry is affected considerably by the satellite elevation angle, and that sea surface altimetry is more accurate with high elevation angles than with low elevation angles [25], [37]. In [25], for one boat-borne GNSS-R altimetric experiment, the relationship between the GNSS-R results and sine of the elevation angle can be described as reciprocal. Since influences of elevation angle on GNSS-R altimetry are important, we here make an initial attempt to consider a simple combination model for GNSS-R altimetry using elevation angle weighting.
In application of GNSS-R to sea surface altimetry based on a single GNSS, considering a multisystem GNSS combination, owing to the different orbit hybrid constellations used in BDS and QZSS, the code noise and errors of different satellites will be asymmetric. Therefore, the above standard least squares method in (6) will not be suitable for the combination of different codes. Here, we introduce the weighted least squares iterative algorithm, which can provide a model that is largely unaffected by spotty errors. By introducing the weight coefficient matrix to adjust the weights between different satellites, the optimal combined RH can be estimated more reasonably.
Currently, the most commonly used a priori weighting methods include equal-precision weighting, signal-to-noise ratio weighting, and elevation angle weighting [38], [39]. In the process of delay determination in the former section, the time delay is determined by using the peak value of the delay waveform. However, due to the nonspecular scattering from the sea surface in the actual scenario, the delay obtained from the reflected waveform using the peak has a derivation compared with the actual specular delay [36]. This derivation is sensitive to sea surface roughness. However, the roughness of the reflection surface is related to the incidence angle and the wavelength of the GNSS signal. The higher the satellite elevation angle, the rougher the reflecting surface will be relative to the GNSS signal, which means that the specular point is farther away from the peak point of the delay waveform, that is, the bigger the derivation of the delay. This phenomenon of coastal GNSS-R code-delay altimetry can be called the elevation angle influence. Therefore, considering influences of the elevation angle, the a priori weighting method based on satellite elevation angle was adopted in this study. Existing stochastic models based on satellite altitude angle mainly include the exponential function, trigonometric function, and linear function models. The more commonly used trigonometric function model was adopted here, and the variance of the path-delay observation between satellites can be expressed as follows: where i = 1, 2, . . . , n represents the ith satellite or signal for which the observed variance of the unit weights (σ 0 ) under the same system is the same. Therefore, the corresponding weight coefficient matrix can be written as follows: This matrix is an n × n diagonal matrix, where each weight on the diagonal corresponds to the corresponding observation. According to Li et al. [27], the precision of GNSS-R codedelay altimetry depends on the transmitted GNSS signals, SNR, and the delay determination method. Moreover, the altimetric accuracy depends on the ionospheric and tropospheric delay corrections, the response of the receiving antenna on different signals, and the sea state. All of these impact factors are approximately dependent for different GNSS satellites because of the different reflection geometries, which correspond to the different specular points on the sea surface for different GNSS satellites. Therefore, it is assumed that there is no correlation between the measurements from different GNSS satellites. According to the weighted least squares definition, the RH can be expressed as follows: where A is the abovementioned coefficient matrix. The initial value of sea surface height can be estimated using the standard least squares strategy, i.e., h 0 = (A T A) −1 Al, and the final optimal RH can be obtained after several iterations. It should be noted that although the nominal specular points are different for different PRNs, the distances among them are not far than 10 m for our ground-based experiment. In this study, the heights of these points are considered to be the same. So, in our combinations, the heights' differences of different specular points at a single epoch were ignored.

III. EXPERIMENTAL RESEARCH
To investigate the performance of GNSS-R code delay in sea level monitoring with different signal modulation methods, a dual-antenna ground-based GNSS-R sea surface altimetry experiment was undertaken, in which raw intermediatefrequency data were collected at a sampling rate of 40 MHz. The experimental data were processed and analyzed using the self-developed GNSS-R altimetry software-defined receiver and inversion software. However, we mainly focus on GNSS-R altimetry performance evaluation based on L1C and B1C ranging codes, which are new signals modulated on frequency at 1575.42 MHz in this study. However, only three GPS satellites broadcast L1C signals and no GLONASS satellites broadcast similar signals. Therefore, we focus on the observations from QZSS and BDS system. The GNSS-R sea level results based on the L1C/A and L1C code of QZSS and the B1C code of BDS were obtained simultaneously. The height values of shore-based synchronous observations were obtained using the radar gauge to evaluate the accuracy of the GNSS-R retrievals. Table II presents details of the experimental equipment parameters.
The experiment was conducted on July 10, 2020, on a coastal trestle bridge (37°32 2.60 N, 122°2 44.06 E) in Weihai, Shandong Province, China. During the observation period, the average wind speed was approximately 2 m/s and the sea state was reasonably calm. The experimental location and the installed equipment are shown in Fig. 2(a) and (b). The foremost point of the trestle bridge is approximately 100 m from the coastline. Fig. 2(c) shows that the vertical upward direct antenna and the vertical downward reflected antenna were oriented on the same plumb line, and that they faced the southern sea area to ensure more usable satellites. The height of the antenna above the water is 2-5 m. As shown in Fig. 2(d), the ground-based radar gauge with 26 GHz working frequency was fixed in a position approximately 15 m horizontally from the reflector antenna for external verification. The vertical distance between the reflector antenna and the radar gauge was 1.153 m, which was measured using a Total Station. The output frequency of the radar gauge is 1 Hz, and the ranging accuracy is ±3 mm, which can be used for the validation of the performance of the GNSS-R sea surface height solutions.

IV. RESULTS ANALYSIS
The J01, J02, and J03 satellites of QZSS were captured during the observation period. Meanwhile, five medium earth orbit satellites and one IGSO satellite of BDS were also successfully searched. Therefore, the RH was calculated from the code delay between the direct signal and the reflected signal of the above visible satellites. Taking the visible satellite as J01 as an example, the code delay and related power waveform changes of L1C/A and L1C codes are shown in Fig. 3. The SNRs of the direct and reflected waveforms for L1C/A signals are 18.68 dB and 14.13 dB; for L1C signal, they are 24.11 dB and 20.07 dB, respectively. It can be seen from Fig. 3(a) and (b) that the correlation power waveform of the L1C code is cleaner (less noisy) than that of L1C/A, and that the main peak of the L1C code is sharper. Indeed, the smooth direct and reflected waveforms of L1C are partly due to the 10-ms coherent integration. However, the smoother waveforms and narrower ACF main peaks of L1C comparing with L1 C/A, which indicates better ranging sensitivity, are due to the advanced signal processing technique or wideband GNSS signals [27]. Moreover, the improvement of the ranging sensitivity overcomes the other degradation factors, resulting in an improved ranging performance of L1C signals. Because the satellite signal will suffer from power loss after being reflected by the sea surface, the relative power value of the reflected signal of the L1C/A code and the L1C code is smaller than that of the direct signal. In Fig. 3(b), because the binary offset carrier modulation method is used in the L1C code, there are two secondary peaks near the main peak of the related power waveform, and the bandwidth is greater than that of the L1C/A code.
Further analysis reveals that the received signal power can be affected by the antenna beam angle and gain. Because the beam angle is inversely proportional to the signal gain, when the beam angle increases, the signal gain becomes smaller, and vice versa. In the case of a satellite with a high elevation angle, the signal received by the receiving antenna will be within the range of higher gain, and the received signal will have a larger relative power value; the converse will be true for a satellite with a low elevation angle. Therefore, because the difference between the position of the direct antenna and the reflective antenna will lead to a change in the gain range, the positions of the two antennas should be set appropriately to ensure that the gain range of the signal remains as large as possible.

A. Analysis of QZSS Results
The change sequence of the satellite elevation angle relative to the station during the observation period is shown in Fig. 4.
It can be seen from Fig. 4 that the range of the available satellite altitude angle of J1 is 0°-65°, the range of the available satellite altitude angle of J2 is 65°-0°, and the range of the available satellite altitude angle of J3 is 80°-68°, which is greater than that of J1 and J2. It can also be found that during the observation period, the trajectory of J3 is near apogee, while the trajectory of both J1 and J2 is near perigee; therefore, the change rate of the satellite elevation angle of J3 is less than that of J1 and J2. On the basis of the code path delay output by the GNSS-R altimetry software receiver and the satellite elevation angle, the RH of the L1C/A and L1C code can be calculated using (2), respectively. The resultant height changes are shown in Fig. 5, where blank areas indicate no observed data during the observation period due to the failure of equipment and limit elevation mask. As the amount of the IF data was very large and the recorder cannot work for a long time, so there were several discontinuities. The results illustrated in Fig. 5 reveal that the maximum change of sea level is approximately 2 m. Moreover, it is evident that the height measurement results of the L1C/A code and the L1C code fluctuate above and below the observation value of the radar gauge. The height measurement results reflect the change of tide level, verifying the suitability of height measurements obtained using the L1C and L1C/A signals with our software receiver. The altimetry results shown in Fig. 6 reveal that the dispersion of the height measurement results of the J1, J2, and J3 satellites, obtained using the L1C/A code, is greater than that of the results obtained using the L1C code, and that the system deviation of the height measurement results is also greater. The height measurement results obtained from the visible satellite using the L1C/A code and the L1C code are given in Table III for comparison. According to Li et al. [27], the assessment of the altimetry performance consists of the characterizations of both random error (precision) and the systematic biases (accuracy) of the altimetric results. Main sources of the bias have not been learned and corrected in this study; therefore, we mainly focus on the precision not the accuracy of GNSS-R sea level estimations. In the following statistical tables (see Tables III-VII), although the biases were also listed in the table, the root mean squares (RMSs), standard deviations (STDs), and R 2 values are all calculated using the dataset with the bias removed.  Similarly, the RMS error and the STD of the sea level derived from the J3 L1C/A are 0.634 and 0.597 m, respectively, and those derived from the J3 L1C code are 0.563 and 0.313 m, respectively. The error of the J2 height measurement result is within ±3 m, which is bigger than that of J3. The experimental results also show that the L1C code phase height measurement results of J2 and J3 are more stable than those of the L1C/A code, and that the measurement accuracy of the L1C code can be increased by approximately 24% and 11%, respectively, in comparison with that of the L1C/A code.
It can also be found from Table III that  show that the L1C code height measurement results of J1 and J3 are more stable than those of the L1C/A code, and the accuracy of the L1C code measurement can be increased by approximately 41% and 50%, respectively, in comparison with that of the L1C/A code. Further analysis shows that in the same observation period, the altitude angle of the J3 satellite is notably greater than that of both J1 and J2. In comparison with both J2 and J1, the accuracy of the L1C/A code and the L1C code measurements of J3 can be increased. Chu et al. [37] also conducted GNSS-R altimetry using these QZSS observations, and almost the same RMS values of sea levels were calculated with our results in Table III.

B. Analysis of BDS Results
During the observation period, the change sequences of the satellite elevation angle of the visible five medium earth orbit satellites and one IGSO satellite of BDS relative to the station  were shown in Fig. 7. On the basis of the B1C code path delay output by the altimetry software receiver and the satellite altitude angle, RH can be calculated. The sea surface altimetry result sequence is shown in Fig. 8, and the blank areas indicate no data during the observation period.
The results illustrated in Fig. 8 show that the sea level measurement results based on the B1C code delay fluctuate up and down in comparison with the radar gauge observation value. The height results reflect the sea level changes, verifying the usability of the software receiver based on the B1C code-delay measurement. The height measurement results of the satellite at the same time are shown in Fig. 9 and presented in Table IV for  comparison. It can be found from the results shown in Fig. 7 that there were four visible satellites, i.e., C25, C34, C38, and C43, during the  The measurement precision of the C33 and C41 satellites is substantially lower than that of the C38 and C43 satellites. The statistical results of GNSS-derived sea levels compared with the in situ measurements are given in Table IV. The altimetry precision of the code of C38 and C43 was improved from the first epoch to the second epoch with the increases of satellite elevation angle. Therefore, on the basis of the altimetry result sequence of the different satellites, it can be determined that larger elevation angles produce measurements with higher accuracy in the same observation period. Because the noise in the code path delay is inversely proportional to the satellite altitude angle, the effect of noise in the code path delay is smaller for larger satellite altitude angles. The accuracy of sea surface height measurements will be improved with the increase of the signal-to-noise ratio of the satellite signal. Moreover, by combining the shore-based altimetry results of the L1C/A and L1C codes of QZSS as well as the B1C code of BDS, it can be found that the residual sequence of the shore-based altimetry result has a small range of fluctuation when the satellite elevation angle is large. When the elevation angle is small, the fluctuation phenomenon of the residual sequence is more obvious. Moreover, the direct signal from the zenith direction enters the reflective antenna during the shore-based experiment, and the multipath effect caused by the superposition of signals of different polarization modes affects the result.

C. Analysis of Joint QZSS and BDS Results
Because the BDS B1C signal and the QZSS L1C signal share the same code rate, same modulation method of data plus pilot component combination, and the same signal cycle length of 10 ms, the B1C and L1C signals were selected for the combined BDS and QZSS system and their variance of the unit weights (σ 0 ) was set to be equal.
To verify the precision of the sea surface high survey based on weighted least squares, the GNSS data of the QZSS L1C and the BDS B1C measurements were collected to perform the calculations using (9) separately. The weighted least squares sea surface altimetry results of the QZSS L1C code and the BDS B1C code and their residual sequence are shown in Figs. 10 and 11. By comparing with Figs. 8 and 9, it can be seen that the observation results from a combination of several satellites in the section are more stable and that the measurement precision is improved substantially than that from a single satellite.
The RMS errors and STDs of the combined height measurement results of QZSS are given in Table V. During the period 08:00-10:00 LT, the RMS error of the combined height measurement is 0.568 m. For satellite J3, the height measurement result is broadly similar, but the precision is increased by approximately 18% in comparison with that of J2. During the period 13:30-15:00 LT, the QZSS combination has high precision of 0.46 m, which is slightly lower than that of the height measurement results of J3. In comparison with J1, the  precision is increased by approximately 37%. Additionally, the STD is <0.35 m for the single-system combination and the error range is also within 1.6 m. The results show that the height measurements are more stable with the weighted combination. The RMS errors and STDs of the combined height measurement results of the BDS single system are presented in Table VI. During the period 08:00-10:00 LT, the precision of the combined height measurement is 0.651 m, which is increased by 3%-55% in comparison with the height measurement results of C25, C34, C38, and C43. During the observation period 12:00-15:00 LT, the precision of the combined height measurement is 0.697 m, which is slightly lower than the results of C38 and C43, but increased by approximately 55% and 47% in comparison with that of C33 and C41, respectively. The comparison results show  that the weighted estimation model can effectively improve the precision and stability of code-delay sea surface measurements.
To perform the joint calculations, intermediate-frequency QZSS and BDS measured data were collected under the same shore-based platform conditions. Then, the results were compared with the solution results of the single BDS and QZSS system. The altimetry results and the residual sequence of the dual-system combinations are shown in Fig. 12(a) and (b), respectively. For comparison, the residual sequence of the combined BDS and QZSS is shown in Fig. 12(c). For ease of comparison, fixed values of 2 and 4 m were added to the residual sequences of the QZSS and BDS altimetry results, respectively, shown in Fig. 12(c). Between 12:00 and 13:00, there appears a relatively large negative derivation for altimetry results in Fig. 12. During these epoch, for QZSS, there are retrievals from J03, while for BDS, there are retrievals from C33 and C43. There are negative derivations for both BDS results and QZSS results, for which we speculated that the reason is deterioration in sea state caused by a sudden increase in wind speed. The negative derivations for BDS are larger than QZSS because the elevation angles of C33 (20°-40°) are lower than that of J02. For this  Table VII. The precision of the combined height measurement is approximately 0.49 m, and the corresponding STD is approximately 0.38 m. In comparison, the precision of the QZSS and BDS single-system measurements is approximately 0.53 and 0.67 m, respectively, i.e., lower than that of the combined observation results. The results show that the combination of multiple systems can effectively improve the precision of code-delay sea surface height measurement. It can also be found that the largest error of the combination system is derived from BDS. Previous analysis revealed that measurement precision corresponds to the satellite elevation angle of BDS. The measurement precision of QZSS is slightly higher than that of BDS because the elevation angle of the satellites of QZSS is slightly higher than that of the satellites of BDS.
We also calculated the R 2 parameters with a linear fit between the GNSS-R retrievals and the radar gauge data. For retrievals from single QZSS satellite, except J03 L1 C/A and L1C signals, the R 2 parameters for the other satellites are less than 0.40. However, for J03 L1 C/A and L1C signals, the correlation coefficients are 0.4450 and 0.9297, respectively. This is because of the fewer retrievals of the J01 and J02, as well as the poor meter-level precision. For example, for J03, there are retrievals for all the observation epoch; however, for J01 and J02, there are respective discrete retrievals with no more than 20 min and 1 h. For BDS satellites, C41 shows the worse R 2 value, the reason for which is the same with that of J01 and J02. The correlation coefficients of other BDS results are 0.4459 to 0.8766. However, after combining all the QZSS L1C signals, the R2 parameter reaches up to 0.9252. While for BDS, the R2 parameter reaches up to 0.8891. The R 2 values of the combinations are bigger than that of each satellite, which indicates the superiority of the proposed combination. While for the combination of all the satellites, including all QZSS and BDS satellites, the R 2 parameter is 0.8932, which is bigger than the combination of all BDS satellites but smaller than the combination of all QZSS satellites (see Table VII and Fig. 13). The combination of all the satellites does not show an improvement of R 2 parameter; this is because that there are retrievals of six BDS satellites and only three QZSS satellites are involved in the combination, which makes the R2 parameter of the combination of all the QZSS and BDS satellites to show the similar value to that of BDS satellites. This phenomenon can also be seen in Fig. 13, in which the retrievals show a similar distribution in Fig. 13(b) and (c). However, this value still outperforms 99% values of single satellite, which indicates the correctness of the combination. V. DISCUSSION Because of the limit of the code chip length, the precision of GNSS-R code-level altimetry is generally at the meter or decimeter level at best. Here, we analyzed GNSS-R code-level altimetry using QZSS L1C/A and L1C signals and BDS B1C signals simultaneously, and proposed a simple combination method only considering the influence of elevation angle to improve altimetry precision. From the proposed combination method, the derived sea levels become more stable and the RMS error values showed improvement than that from a single satellite.
For GNSS-R code-delay altimetry, the SNR, the signal bandwidth (the shape of the ACF), can also affect the ranging precision and thus the altimetry performance. However, it is found that the SNRs of reflected signals have an approximate linear relationship with the elevation angle of GNSS satellites. Fig. 14 depicts the SNRs and changes of elevation angles of different satellites and different signals. From this figure, it can be found that with the increase of the elevation angle, the long-term trend of SNR shows an approximate increase. Moreover, in [10], it is found that the precision of derived RHs changes with the reciprocal of the sine of elevation angle. Therefore, here in this article, we think that since the SNRs changed similarly with elevation angle, only using the elevation angle as a weight can also compensate for the changes of SNRs. This is our preliminary results, and the actual influences of SNRs on the retrieval precision will be more complicated for which further work will be needed in our following study.
Because only the influence of elevation angle was considered in the proposed combination method, the improvement was marginal for some combinations due to the influence of many other errors. GNSS signals translating from a GNSS antenna to a receiver, either directly or indirectly following reflection by the sea surface, will be affected by various error factors. Therefore, the actual path delay of a reflected signal relative to a direct signal includes the geometric path delay, atmospheric delay, satellite and receiver clock difference, and several other delays attributable to factors, such as sea surface roughness, antenna errors, hardware delays, electromagnetic bias, and the curvature of the earth [40], [41], [42]. If there is no bias for a sequence, its STD and RMS will be the same. In our study, there are biases for the solutions from different satellites between the in situ measurements. The biases vary from different sea states and elevations. In the next step, a model for computing the bias will be developed. In sea surface altimetry using reflected signals, the clock difference of the satellite and the receiver can be canceled without causing observational errors because the direct and reflected signals received on the ground come from the same satellite and are processed by the same receiver. Antenna error delay refers to the projection of the baseline vector between the direct and reflective antenna phase centers in the direction of the line connecting the satellite and the direct antenna phase center. This projection can be determined from the positional relationship between the reflective antenna, direct antenna, and altitude and azimuth angle of the satellite. When the vertical distance from the receiving antenna to the sea surface is small, it can be assumed that the earth is horizontal; otherwise, the influence of the curvature of the earth must be considered. Analysis on how to determine the true sea surface height when considering the curvature of the earth is available in the literature [43]. Additionally, the effect of atmospheric delay can be ignored for shore-based and airborne platforms because the propagation paths of direct and reflected signals are similar owing to the low observation heights. For spaceborne observation platforms, the propagation paths of direct and reflected signals do not coincide, but the delay caused by the ionosphere and the troposphere must be calculated. When the sea state is poor, there are a large number of diffuse reflection points around the specular reflection point. The cumulative effect on the reflected signal causes the one-dimensional correlation power peak point of the reflected signal shift backward. Additionally, (2) simply indicates that the satellite elevation angle is inversely proportional to the noise in the path delay, while in fact, during data processing, we also found that when the satellite altitude angle is lower than 30°, the path-delay noise will be amplified. Thus, when the satellite altitude angle is low, the calculation result might be affected substantially by noise, and the corresponding precision of the sea surface measurement will be reduced. With the increase of the satellite altitude angle, the path-delay noise is gradually reduced, and the precision of sea surface measurement will be improved correspondingly. Moreover, except elevation angle, the electromagnetic bias in GNSS-R altimetry is also influenced by the signal frequency, wind speed, and wind direction. Therefore, due to the joint influence of above errors, the precision improvement of combining different satellites in a single GNSS was not that obvious. To reduce delay error and to improve retrieval precision, further research should be undertaken to diminish the impact of the errors discussed above.
Bias values were also included in Tables III-VII. From these values (most of them are positive values), we can see that the GNSS-R technology overestimated the RHs comparing with the radar gauge. After combining the satellites within one system, the bias for sea levels estimated from QZSS and BDS signals are 36 and 31 cm, respectively. These values are almost in the same magnitude and may imply the existence of system error. This bias is likely due to the phase center deviation of the antenna and the measure center deviation of radar gauge. However, further confirmation on this bias should be carried out in the future.
In Table VII, the optimal precision of GNSS-R sea level estimation at 1 Hz for our experiment is about 49 cm. This precision is too coarse and it is difficult to meet the actual application needs. Gao et al. [25] verified that the moving average method can be used in GNSS-R altimetry to obtain sea level series with centimeter-level sea levels. Therefore, the tradeoff between altimetry precision and temporal resolution was further evaluated by using the moving average method with different moving window of 1, 5, 10, 30, and 60 min on the final data of QZSS, BDS, and their combination. Table VIII lists the statistical results of code-delay altimetry using a moving average method with different moving windows. From this table, we can see that when the time of moving window is less than 30 min, the precision of all the sea level series can be improved as the increase of the window. However, after the 60-min window was adopted, the precision of retrievals is deteriorated instead. That is, it is not that the longer the moving window, the higher the precision. Moreover, combined with the results of Gao et al. [25], we can also conclude that the turning point of the window time is not fixed on 30 min but changes with the environment of experiment, because there are more BDS results than that of QZSS, and the RMS of BDS-derived sea levels is bigger than that of QZSS (see Table VII). The RMS of the combined sea levels after using a moving average method is between the both, and more inclined to the former. For the experiment in this study, an optimal sea level estimations with RMSE of 11-13 cm, and correlation coefficient bigger than 0.98 was realized. Gao et al. [25] also conducted GNSS-R altimetry using BDS B1C signals using different observation data from another day at the same experimental station. They used the BDS observations from different satellites with this study, and they obtained an optimal sea level series with RMS of 5-9 cm for BDS B1C signals. The precision of sea level derived from the work of Gao et al. [25] is better than the results in this study. However, since the experimental data used are different (obtained from different satellites on different days), it is meaningless to compare the values of the RMS. In this study, since the influences of error of the elevation angle was the only one that was taken care in this study, if most errors would be corrected in the future, centimeter-level sea level estimation from GNSS-R code-delay altimetry is just around the corner.

VI. CONCLUSION
This article tested the GNSS-R altimetry using interoperable signals with frequency of 1575.42 MHz, i.e., L1C/A and L1C from QZSS and B1C from BDS simultaneously for the first time. The results of a dual-antenna shore-based experiment were discussed to explore the precision of sea surface measurement based on BDS and QZSS. Through comparison with results obtained using an in situ radar gauge, sea levels with RMS errors of 0.578-1.563, 0.395-0.726, and 0.634-1.22 m were derived using the BDS B1C signal, the QZSS L1C, and L1 C/A signals, respectively. After using a moving average method, optimal sea levels at 30-min temporal resolution with RMS of 11-13 cm were realized. By comparing the results of QZSS L1C and L1 C/A signals, which are QMBOC and BPSK modulation, respectively, we found that the performance of QMBOC is better than that of BPSK on sea level estimation. In addition, the precision of GNSS-R code-level altimetry is related to elevations. In this study, the signals from BDS MEO were used and their elevations change quickly, so the biases were more obvious. Because of the higher elevation angle of QZSS satellites, the precision of sea level measurements obtained from QZSS signals was slightly higher than that derived from BDS signals.
In order to further improve the precision of GNSS-R codelevel altimetry, we also proposed a preliminary combination method of multiple signals using elevation angle weighting to weaken its influence and improve sea level retrieval precision. Then, the experimental results using the proposed combination method revealed that for most combinations of satellites and signals, the precision of the estimated sea level was improved substantially. After combining all the QZSS and BDS signals, the RMS of estimated sea level is improved from 0.525 and 0.672 m of QZSS signals and BDS signals to 0.489 m. Although the proposed combination model improved the sea level estimation precision, due to other unconsidered influence elements than elevation angle, not all the combinations in this study show obvious precision improvement. This method is only an initial attempt and further research on combining all the influence elements into our model will be on the way. However, this work verified that elevation angle is an influence element for GNSS-R code-level altimetry and will afford a basic idea for future error correction of GNSS-R sea level estimation. Moreover, in this study, the peak of the correlation waveform was identified as the arrival of the reflected signal, and the peak value was obtained using cubic spline interpolations from the peaks positions. However, there are two secondary peaks near the main peak of the related power waveform, such as L1C code, which may lead to misparsing of the peak and the poor retrieval precision of sea level. Therefore, in the future, we will try to fit a theoretical BOC waveform to the full experimental waveform, including secondary peaks, to calculate the arrival of the reflected signals instead of relying only on the data near the main peak.