GNSS-R sea surface height measurement method based on variable step size line segment intermediate frequency data processing

Abstract

The present application relates to the technical field of satellite altimetry and ocean surveying and mapping, and in particular to a GNSS-R sea surface height measurement method based on variable step line segment intermediate frequency data processing, which can effectively improve the measurement accuracy and efficiency, comprising obtaining satellite-borne GNSS-R original intermediate frequency data, calculating mirror reflection points through a new variable step line segment (NVSLSM), generating a delay Doppler map (DDM), and performing sea surface height inversion. Compared with the prior art, the NVSLSM is used to process the mirror reflection points, the calculation time is increased by about 2.04 times and 1.47 times compared with the line segment bisection method and the QH method, and on the premise of ensuring the convergence accuracy, the stable solution can be reached through fewer iteration steps, thereby reducing the calculation cost per sample; in the delay Doppler map generation link, the overall calculation efficiency is increased by more than about 90%.

Description

Technical Field

[0001] This invention relates to the fields of satellite altimetry and marine surveying, specifically a GNSS-R sea surface altimetry method based on variable step-size line segment intermediate frequency data processing that can effectively improve the accuracy and efficiency of measurement. Background Technology

[0002] Using ocean gravity field information to assist inertial navigation is one of the key technical approaches to solving the error accumulation problem in underwater inertial navigation systems. Satellite altimetry technology retrieves the ocean gravity field by acquiring sea surface height data. Compared with traditional ship-based measurement methods, this greatly expands the coverage and spatial resolution of gravity field data, providing a richer and more accurate data foundation for gravity-matched navigation. Currently, spaceborne GNSS-R technology has been widely applied in various fields of Earth remote sensing. However, most research still relies on delayed Doppler maps (DDMs), or L1-level data products, generated on-board. It is worth noting that, in addition to such standard data products, some spaceborne GNSS-R missions also synchronously store raw intermediate frequency (IF) data. This data is acquired directly after analog-to-digital conversion without any further on-board digital processing, preserving the original characteristics of the signal to the maximum extent. The acquisition and utilization of this raw IF data provides an important data foundation and verification support for the design optimization and new algorithm development of future spaceborne GNSS-R missions. Higher resolution DDM data products can be generated based on raw intermediate frequency data, which can improve the application effect of GNSS-R in land and ocean remote sensing.

[0003] In spaceborne GNSS-R intermediate frequency (IF) signal processing systems, software receivers have become one of the most prevalent implementation methods. Compared to traditional receivers with fixed logic and dedicated hardware, software receivers can complete acquisition, tracking, correlation, and observation generation processes on general-purpose computing platforms, thus more flexibly adapting to different signal systems, mission modes, and diverse application requirements. Existing research shows that introducing parallel computing architectures can significantly improve the processing throughput of software receivers. However, current GNSS-R signal processing and engineering implementations are mainly geared towards ground or airborne platform scenarios, while the development and optimization of end-to-end processing links for spaceborne platforms are relatively insufficient. Spaceborne missions typically involve higher data generation rates, longer-term continuous observations, and stricter constraints on payload computing power and downlink bandwidth, making the computing speed and data processing capabilities of software receivers one of the core bottlenecks restricting the further development of GNSS-R. Therefore, optimizing efficient IF processing architectures and algorithms for spaceborne conditions has clear engineering necessity and research value. Summary of the Invention

[0004] This invention addresses the shortcomings and deficiencies of existing technologies by proposing a GNSS-R sea surface altimetry method based on variable step-size line segment intermediate frequency data processing. This method accelerates and robustens the solution process for specular points (SPs) in the geometric auxiliary parameter generation stage, thereby improving computational efficiency under large-scale epoch channel conditions and alleviating the data accumulation and real-time processing pressure under continuous on-orbit observation conditions.

[0005] This invention achieves its purpose through the following measures:

[0006] A GNSS-R sea surface altimetry method based on variable step-size line segment intermediate frequency data processing is characterized by the following steps: acquiring raw intermediate frequency data from spaceborne GNSS-R, calculating specular reflection points using a novel variable step-size line segment NVSLSM method, generating a delayed Doppler image (DDM), and performing sea surface altitude inversion. Specifically, the calculation of specular reflection points using the novel variable step-size line segment NVSLSM method includes the following steps:

[0007] Step 1-1: Input the receiver and transmitter positions. The position information is obtained from the navigation information. Let:

[0008] (1),

[0009] in, It is the first The difference between the incident angle and the reflection angle at the next iteration. and The calculation formula is as follows:

[0010] (2),

[0011] (3);

[0012] Step 1-2: Initialize the specular reflection points using formula (2):

[0013] (4),

[0014] (5),

[0015] (6),

[0016] In the WGS84 coordinate system, the major axis The eccentricity is 6,378,137 m. It is 0.08181919084264;

[0017] Steps 1-3: Calculate the initial difference in the number of mirror reflection points:

[0018] (7), (8),

[0019] Wherein, the initial iteration step size and They are 0.4 and 0.6 respectively.

[0020] Steps 1-4: Calculate the step size for subsequent iterations:

[0021] (9);

[0022] Steps 1-5: Recalculate the difference between the angle of incidence and the angle of reflection:

[0023] (10) (11);

[0024] Steps 1-6: If conditions If this condition is met, then the iteration exit condition is satisfied. If it is the final specular reflection point, otherwise return to steps 1-4 and iterate again. for .

[0025] The generation of delayed Doppler maps (DDM) described in this invention is achieved by constructing an FFT-based GPU parallel processing architecture, which includes:

[0026] Step 2-1: Signal generation and open-loop tracking. The initial target motion parameters are provided by the open-loop tracking module. The signal generator generates the transmission signal based on the open-loop tracking information and splits the signal into multiple processing branches: one is sent to the reference signal replication module Replica, one is sent to the pseudo-random code modulation module C / A Code, and one directly participates in the mixing.

[0027] Step 2-2: Receive the echo signal (ReflectSignal) from the target through two or more receiving channels. Each reflected signal is mixed (multiplied) with the corresponding local replica signal (Replica) in a mixer.

[0028] Steps 2-3: C / A code modulation and matched filtering. Two or more parallel C / A Code modules receive the raw signal from the SignalGenerator, apply C / A pseudo-random code modulation, and the modulated signal also enters the mixer. It is multiplied with the reflected signal to realize code correlation processing, and the output signal is processed by FFT.

[0029] Steps 2-4: FFT processing is implemented through at least two parallel FFT modules. The FFT output then enters the IFFT module, which has at least two channels. The outputs of the two or more IFFT channels are combined and then enter the Incoherent & Capture module.

[0030] Steps 2-5: The incoherent accumulation module fuses the multi-channel IFFT output signals, identifies significant target points, and outputs them as DDM distance-Doppler maps.

[0031] In this invention, sea surface height inversion requires re-tracking the position of the specular reflection point on the waveform. This allows for the extraction of time delay values ​​while eliminating path delay errors caused by sea surface roughness. Based on the Der method, the reflected signal corresponding to the specular reflection point is the part that contributes the most to the power waveform, i.e., the time delay corresponding to the largest first derivative of the power waveform, as shown in the following formula:

[0032] (13) (14)

[0033] The path difference between the direct signal and the reflected signal is expressed as:

[0034] (15)

[0035] in, For the observed bistatic delay, and The direct and reflected signal paths are identified. Based on the positions of the receiver and transmitter, the coordinates of the specular reflection point on the WGS84 ellipsoid are calculated. After error correction, the bistatic delay in the model is obtained. The calculation is as follows:

[0036] (16)

[0037] in, The geometric model delay between the direct and reflected signals is used to obtain the height of the specular reflection point relative to the WGS84 ellipsoid:

[0038] (17), among which, It is the angle of incidence of the reflected signal at the point of reflection on the mirror.

[0039] The mean sea level height data in this invention is derived from the DTU21 model, and tidal correction uses the TPXO9 model. The two models are superimposed to form a dynamically changing sea level height reference field. The accuracy and reliability of the inversion results are then verified. The sea level height data is expressed as follows:

[0040] (18), of which, express The actual sea level height of the model after tidal correction. This represents the average sea level height calculated by the model. This indicates the tidal correction calculated by the model.

[0041] To ensure the accuracy and reliability of the sea surface height inversion results and further improve the performance of the inversion model, the present invention also implements the following steps: (1) removing land reflection point data and retaining only the observation values ​​of specular reflection points located in the ocean area; (2) excluding power waveforms with a signal-to-noise ratio of less than 10 dB to reduce the impact of low-quality data on the inversion results;

[0042] (19), among which, It is the maximum power of the power waveform. It is the average power of the power waveform, which is taken from the zero-Doppler slice region in the Delayed Doppler Map (DDM).

[0043] Compared with existing technologies, this invention employs NVSLSM for specular reflection point processing, achieving a computation time improvement of approximately 2.04 times and 1.47 times compared to the line segment bisection method and the QH method, respectively. Its average number of iterations is 9, which is approximately 67.2% and 43.8% less than the line segment bisection method and the QH method, respectively. NVSLSM can achieve a stable solution with fewer iterations while maintaining convergence accuracy, thereby reducing the computational cost per unit sample. In the Delay Doppler Map (DDM) generation stage, an FFT-based GPU parallel processing architecture is constructed. Compared with the CPU serial FFT implementation, the GPU processing time is reduced from 143.00 s to 12.64 s, with an overall computational efficiency improvement of over 90%. When the data block size increases to over 50 ms, the processing time stabilizes in the 8–9 s range, significantly improving the throughput of mid-frequency data processing. Attached Figure Description

[0044] Appendix Figure 1 This is a flowchart of the present invention.

[0045] Appendix Figure 2 This is a schematic diagram of the geometric relationship of the mirror reflection points in this invention.

[0046] Appendix Figure 3 This is a flowchart of GPU data processing in this invention.

[0047] Appendix Figure 4 This is the basic principle of sea surface altitude inversion using spaceborne GNSS-R.

[0048] Appendix Figure 5This is a comparison of the convergence characteristics of three methods for solving the specular reflection point, where the horizontal axis represents the number of iterations and the vertical axis represents the normalization error.

[0049] Appendix Figure 6 This invention presents a comparison of the mirror reflection point positioning errors using three methods based on high-precision numerical solutions in this embodiment. Figure 6 In the diagram, 'a' represents the scatter distribution of the ENU level component error. Figure 6 In this context, b is the Euclidean cumulative distribution function of position error; Figure 6 In this context, c is the cumulative distribution function of the horizontal position error; Figure 6 In this context, d is the cumulative distribution function of the vertical error assignment.

[0050] Appendix Figure 7 These are comparison diagrams of DDMs generated by different strategies in embodiments of the present invention, wherein... Figure 7 In the middle, 'a' corresponds to a convolutional strategy. Figure 7 b corresponds to the CPU-FFT policy. Figure 7 c corresponds to the GPU-FFT strategy.

[0051] Appendix Figure 8 This is a comparison of the time consumption of the GNSS-R key computing module on CPU and GPU platforms in this embodiment of the invention.

[0052] Appendix Figure 9 This invention relates to a comparison and residual statistical analysis of altimetry results in this embodiment, where the reference height is composed of the DTU sea surface model and the tidal correction term; wherein... Figure 9 In the middle, 'a' represents the scatter density distribution of the HALF altimetry results and the reference sea level height; Figure 9 In the middle b, the residual histogram of the HALF relative to the reference height is shown. Figure 9 In the middle, c represents the scatter density distribution of the Der altimetry results and the reference sea level height; Figure 9 In the diagram, d is the residual histogram of Der relative to the reference height. Detailed Implementation

[0053] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0054] The CYGNSS satellite constellation mainly consists of eight small satellites in low Earth orbit. The primary onboard data processor is called the Delay Doppler Map Instrument (DDMI), which uses GPS forward-scattered signals to generate delayed Doppler maps of the reflector. In addition to standard DDM data products, the DDMI can also record raw intermediate frequency (IF) data by digitizing radio frequency signals received by two downward-looking antennas and one upward-looking antenna.

[0055] The processing of raw intermediate frequency (IF) data relies on an extended GNSS software receiver specifically developed for GNSS-R applications. This example addresses the computational bottleneck of spaceborne GNSS-R software receivers by proposing a novel variable step-length segment method (NVSLSM) for solving specular reflection points. This method improves convergence efficiency while maintaining numerical robustness through adaptive step size, adapting to the rapid computation requirements under large-scale epochs and complex geometric conditions. In the data correlation processing stage, this example employs GPU parallelization to accelerate the relevant computational processes, migrating operations such as FFT correlation, originally calculated serially on the CPU, to the GPU, thereby improving the processing efficiency of IF data. The overall technical framework of this invention is as follows: Figure 1 As shown.

[0056] Due to the high dynamism of the spaceborne platform and the susceptibility of reflected signals to sea conditions, direct tracking is difficult. To enable open-loop tracking of reflected signals, signal synchronization under spaceborne conditions needs to be considered. The path delays of direct and reflected signals vary significantly due to the relative motion between the satellite and receiver. The direct and reflected signals received by the receiver at the same time do not originate from the same transmission moment. Furthermore, the path delay of the reflected signal changes continuously over time due to relative motion, resulting in non-strict matching of the Doppler waveforms of reflected signals at adjacent moments. Therefore, synchronization of direct and reflected signals is necessary, and this synchronization relies on high-precision specular reflection points. To accelerate and robusten the specular reflection point calculation process in the spaceborne GNSS-R geometric auxiliary parameter generation stage, this example proposes a novel Noval Variable Step Length Segment Method (NVSLSM). The calculation process is as follows: First, the receiver and transmitter positions are input, obtained from navigation information:

[0057] make: (1), where, It is the first The difference between the incident angle and the reflection angle at the next iteration. and The calculation formula is as follows:

[0058] (2),

[0059] (3),

[0060] Use formula (2) to initialize the specular reflection points:

[0061] (4), (5),

[0062] (6), where, in the WGS84 coordinate system, the major axis The eccentricity is 6,378,137 m. It is 0.08181919084264;

[0063] Calculate the initial difference in the number of mirror reflection points:

[0064] (7), (8), where the initial iteration step size and They are 0.4 and 0.6 respectively.

[0065] Calculate the step size for subsequent iterations:

[0066] (9),

[0067] Recalculate the difference between the angle of incidence and the angle of reflection:

[0068] (10) (11) If the condition If this condition is met, then the iteration exit condition is satisfied. If it is the final specular reflection point, otherwise return to (4) and iterate again, where for .

[0069] The GNSS-R observations after processing the raw intermediate frequency data include one-dimensional and two-dimensional time-delay Doppler correlation power, code delay observations, and carrier phase delay observations. The two-dimensional time-delay correlation power is obtained by cross-correlating the reflected signal with a local replica having a fixed Doppler frequency.

[0070] (12), among which, This is the coherent integration time, which is set to 1ms. It is the reflected signal received by the downward-looking antenna at any given time. It is a local PRN copy code at any given time. It is a fixed Doppler center frequency. It is an estimate of the Doppler frequency at the point of mirror reflection. It is the time delay between the reflected signal and the local pseudocode signal. and different Doppler frequency shifts The relevant values ​​below.

[0071] Under spaceborne conditions, the reflected signal is significantly weakened, resulting in a low signal-to-noise ratio. Therefore, it is necessary to improve the signal-to-noise ratio and reduce thermal noise through coherent and incoherent integration. Considering the squared loss introduced by incoherent integration, the incoherent integration time used in this invention is 1 second.

[0072] The computational resources for spaceborne intermediate frequency data processing are mainly concentrated on the correlation operations of reflected signals and DDM generation. GNSS-R related calculations have prominent data parallelism characteristics. GPUs have large-scale parallel computing units and high-throughput storage hierarchies, which are suitable for parallelizing repetitive complex multiplication and addition, frequency domain point multiplication, and batch FFT operators, thereby significantly improving the effective processing volume per unit time.

[0073] Frameworks that use GPUs for computation, such as Figure 3 As shown. To further improve overall processing efficiency, related calculations, window extraction, and incoherent accumulation are performed on the GPU. Storing intermediate data in GPU memory effectively reduces the additional latency caused by data transfer and improves overall processing throughput. Finally, the GPU directly generates and outputs the DDM result, avoiding the frequent readbacks and repetitive processing of intermediate related quantities in traditional implementations.

[0074] The experimental platform employs workstation-level hardware and software configuration. The GPU is an NVIDIA GeForce RTX 2060 (6 GB GDDR6), used for performing reflection signal correlation calculations and delay-Doppler mapping. The CPU handles data input, task scheduling, and result feedback. The software environment is based on NVIDIA CUDA Toolkit 12.2. In the implementation, correlation calculations, correlation window extraction, and incoherent accumulation are performed on the GPU as much as possible to reduce CPU and GPU data transfer and synchronization overhead and improve end-to-end processing throughput. The experimental data sampling rate is 16MHz, the quantization bit depth is 2 bits, the DDM grid resolution is set to 21×101, the incoherent accumulation length is 1000ms, and the number of parallel processing channels is 4. The above platform and parameter configuration are used to evaluate the processing efficiency and performance of the proposed optimization method in the rapid processing of spaceborne GNSS-R intermediate frequency data.

[0075] The core principle of dual-antenna-based sea surface altimetry lies in the precise extraction of the time delay difference between direct and reflected signals through coordinated processing. This time delay difference is then used to achieve altitude inversion based on the determined geometric relationship between sea surface height and the time delay difference. Accurate extraction of the time delay observation is crucial for achieving high-precision spaceborne GNSS-R altimetry. When using GNSS-R technology to measure sea surface height, sea surface roughness causes the peak power position of the reflected signal's time delay power curve to shift backward. Compared to the autocorrelation function of the direct signal, its waveform trailing edge becomes gentler, introducing uncertainty into GNSS-R altimetry. Therefore, it is necessary to re-track the position of the specular reflection point on the waveform to eliminate path delay errors caused by sea surface roughness while extracting the time delay value.

[0076] (1) Der method. The idea of ​​this algorithm is that the reflected signal corresponding to the mirror reflection point is the part that contributes the most to the power waveform, that is, the time delay that corresponds to the largest first derivative of the power waveform. The formula is as follows:

[0077] (13);

[0078] (2) Half Method: This algorithm is a simplified version of the monostatic radar altitude measurement algorithm. At the leading edge of the waveform, the time delay corresponding to the peak power portion is considered the actual time delay. In GNSS-R, this coefficient is typically chosen to be 0.70, i.e., (14);

[0079] Spaceborne GNSS-R altimetry typically employs the dual-delay difference method, determining the height of the reflecting surface relative to the reference ellipsoid by comparing the difference between the bistatic delay and the geometrically modeled bistatic delay. Since the actual reflection path is affected by atmospheric errors, atmospheric error correction is considered during modeling to reduce its impact on sea level measurements. The ionosphere is corrected using the GIM model, and the troposphere uses the UNB3M model.

[0080] like Figure 4 As shown, the path difference between the direct signal and the reflected signal can be expressed as:

[0081] (15), among which, For the observed bistatic delay, and These represent the direct and reflected signal paths, respectively. Based on the positions of the receiver and transmitter, the coordinates of the specular reflection point on the WGS84 ellipsoid can be calculated. After error correction, the bistatic delay under the model is obtained. The calculation is as follows:

[0082] (16), among which, This is the geometric model delay between the direct and reflected signals; ultimately, the height of the specular reflection point relative to the WGS84 ellipsoid can be obtained as follows:

[0083] (17), among which, It is the angle of incidence of the reflected signal at the point of reflection on the mirror.

[0084] The actual data for the specular reflection points were calculated using Wavpy, an open-source GNSS-R data analysis and simulation software repository. The specular geometry section allows users to set the positions and velocities of the receiver and transmitter, and provides a range of functions for more comprehensive characterization.

[0085] Due to the lack of actual observational data, this example uses a combination of a mean sea level height (SGS) model and a tidal model to construct a validation benchmark. The SGS data is derived from the DTU21 model published by the Technical University of Denmark, while the tidal correction uses the TPXO9 model. The two models are superimposed to form a dynamically changing SGS reference field, which is used to verify the accuracy and reliability of the inversion results. Sea level height data. It can be represented as

[0086] (18), of which, express The actual sea level height of the model after tidal correction. This represents the average sea level height calculated by the model. This indicates the tidal correction calculated by the model.

[0087] To ensure the accuracy and reliability of the sea surface height inversion results and further improve the performance of the inversion model, this invention implements a systematic quality control process. The specific steps are as follows: (1) Remove land reflection point data and retain only the observation values ​​of specular reflection points located within the ocean range; (2) Eliminate power waveforms with a signal-to-noise ratio of less than 10 dB to reduce the impact of low-quality data on the inversion results;

[0088] (19), among which, It is the maximum power of the power waveform. It is the average power of the power waveform, which is taken from the zero-Doppler slice region in the Delayed Doppler Map (DDM).

[0089] This example presents a comparison of the computational efficiency of three methods for solving specular reflection points. Overall, the NVSLSM method performs best in terms of average time consumption, showing a relative improvement of approximately 2.04 times and 1.47 times compared to the line segment bisection method and the QH method, respectively. Simultaneously, its average number of iterations is 9, significantly lower than the line segment bisection method and the QH method, reducing them by approximately 67.2% and 43.8%, respectively. This result demonstrates that NVSLSM can achieve convergence with fewer iterations, thus exhibiting higher computational efficiency and lower computational cost in large-scale processing tasks. Analyzing the dispersion and robustness of the results, the standard deviation (Std) of NVSLSM is 0.076 ms, significantly lower than that of the line segment bisection method and the QH method, indicating smaller fluctuations in computational cost, more concentrated performance distribution, and lower sensitivity to different samples or geometric conditions. In contrast, while the QH method outperforms the line segment bisection method on average, its Std is significantly larger, potentially indicating unstable convergence behavior in some scenarios. Based on the above statistical results, NVSLSM performs better in terms of average efficiency, convergence cost, and stability. It is more suitable as a geometry-aided solution module under the conditions of high data rate and long-term continuous observation in spaceborne GNSS-R, so as to support the subsequent intermediate frequency data processing of spaceborne GNSS-R.

[0090] Figure 5 Convergence curves of the normalized error for three methods during the iteration process are presented. The vertical axis uses a logarithmic scale to highlight the differences in error decay rates among the different algorithms. NVSLSM's error rapidly decreases after the 10th iteration. The magnitude indicates that this method can reach the lower limit of numerical accuracy with fewer iterations; the QH method also shows rapid error decay, but requires a relatively higher number of iterations to achieve the same accuracy; in contrast, the line segment bisection method shows a significantly slower error decay, requiring approximately 30 iterations to approach the lower limit. The magnitude indicates that the method has low convergence efficiency under the same convergence increment and is easily affected by numerical rounding issues. The results show that in solving nonlinear geometric equations such as specular reflection points, using an algorithm with stronger local convergence characteristics can significantly reduce iteration overhead, thus making it more suitable for the engineering needs of spaceborne GNSS-R.

[0091] Based on the spatial distribution of specular reflection points calculated by different algorithms globally, it can be seen that the marine portion of CYGNSS mid-frequency data specular reflection points is mainly concentrated in the western Atlantic-Caribbean-Gulf of Mexico, the northern Indian Ocean-Bay of Bengal, Southeast Asia, and the western Pacific Ocean. Their distribution exhibits a banded characteristic related to orbital coverage, reflecting the sampling characteristics of the observation data under temporal and orbital conditions, as well as the geographical coverage of available reflection events. More importantly, the specular point locations given by the four methods on a global scale show a high degree of consistency, indicating that under this dataset and geometric conditions, the estimation differences of the geometric solutions for specular reflection by each algorithm are generally small. The phenomenon of a single color dominating in the figure mainly stems from overlapping plotting and occlusion of results from multiple methods near the same location, and does not indicate method failure.

[0092] The figure shows the geographic distribution of specular reflection points calculated using the Wavpy method, the line segment bisection method, NVSLSM, and QH methods under the same observation dataset. The specular points are mainly distributed in the ocean area and exhibit a banded distribution along the satellite orbit coverage area. The results from each method are highly consistent spatially, but overlapping points cause local areas to show superposition and occlusion. To evaluate the consistency of positioning among different specular reflection point calculation methods, the estimated three-dimensional Euclidean distance between the specular reflection points and the reference ground truth is used as an error metric.

[0093] (20);

[0094] Statistical results show that the three methods are similar in average error levels. The bias of the line segment bisection method and NVSLSM is almost identical, while QH is slightly higher. The global error magnitudes of the three methods differ only slightly. However, in terms of dispersion and overall error indices, the QH method exhibits superior stability, with a significantly lower STD than the line segment bisection method and NVSLSM, and a smaller RMSE. Furthermore, the maximum error of the QH method is significantly lower than the other two methods, indicating stronger suppression of extreme errors in some samples. In contrast, NVSLSM and the line segment bisection method show high consistency in bias, STD, RMSE, and extreme value statistics. While the QH method shows little difference in average error levels, it performs better in error dispersion and extreme error suppression, with a more concentrated error distribution and less sensitivity to a small number of unfavorable geometries or anomalous samples. In summary, while maintaining consistent convergence criteria, NVSLSM provides support for improving computational efficiency with lower iteration costs while maintaining comparable positioning accuracy to the line segment bisection method. The QH method, on the other hand, demonstrates advantages in error stability and upper bound control.

[0095] To further characterize the directional features of the mirror reflection point positioning error in the local coordinate system, this example projects the three-dimensional position error onto the E, N, and U components in the station-centered coordinate system (East-North-Up, ENU). ENU decomposition reveals the dominant components of the error in the horizontal and vertical directions and their asymmetry, providing more explanatory statistical evidence for subsequent synchronous compensation models and geometric prior error propagation analysis.

[0096] (twenty one), (twenty two),

[0097] From the error amplitude distribution, the errors in the horizontal component of all three methods are significantly greater than those in the vertical component. Among them, the N-axis error is the most prominent: the N-axis bias of the line segment bisection method and NVSLSM is approximately... m, RMSE is approximately m; the N-axis bias of the QH method increases further, while the corresponding RMSE remains at the same order of magnitude. In contrast, the U-axis error is generally smaller: the U-axis RMSE of the segment bisection method and NVSLSM is approximately Furthermore, the bias is close to zero; this indicates that under the current dataset and geometric conditions, the error in solving the mirror point is mainly concentrated in the horizontal plane, while the vertical component is relatively stable.

[0098] In terms of method comparison, the segment bisection method and NVSLSM show a high degree of consistency in the three-directional statistics, indicating that under the same iteration stopping and experimental settings, NVSLSM maintains the same directional accuracy characteristics as the bisection method while providing a feasible basis for efficiency improvement. The QH method shows some differences in the horizontal components: its E-axis bias is smaller and its dispersion is lower, but the U-component shows greater fluctuations and a significantly expanded extreme value range, suggesting that it may be more sensitive to the geometric conditions of some samples or have anomalies in the vertical direction. Overall, the ENU component statistics further verify the differences in the error direction structure of different algorithms: NVSLSM and the bisection method have consistent error structures and stable vertical components, while the QH method introduces more obvious vertical component fluctuations while reducing some horizontal errors.

[0099] like Figure 6 As shown in Figure a, the scatter distribution of the three methods in the horizontal error component plane exhibits significant non-uniformity. The errors are mainly concentrated in several dense areas with a small number of outliers, reflecting the significant differences in error levels under different observation geometry conditions. From an overall perspective, the scatter distribution of horizontal errors in NVSLSM and the line segment bisection method highly overlaps, while the distribution range of the QH method is relatively convergent in some areas, suggesting that it has a certain suppressive effect on the dispersion of horizontal errors.

[0100] Figure 6 b and Figure 6 c in the figure gives the three-dimensional error respectively with horizontal error The CDF curves of the three methods are shown. The overall CDF shapes of the three methods are similar, indicating that the positioning errors are on the same order of magnitude in an average sense. Among them, the QH method's CDF is relatively more to the left in the middle to high percentile range, indicating that it can achieve smaller 3D and horizontal errors on most samples and has a certain suppression effect on samples with larger errors. Meanwhile, the CDFs of NVSLSM and the line segment bisection method almost overlap, further verifying the consistency of their statistical characteristics in positioning accuracy.

[0101] Figure 6 The value d in the middle shows the vertical error. The CDF differences are observed. It can be seen that NVSLSM and the segment bisection method maintain a high degree of consistency in the vertical error distribution, with most samples showing vertical errors concentrated within a relatively small range. In contrast, the vertical error CDF of the QH method shifts significantly to the right and exhibits a longer tail, indicating greater fluctuations and a higher upper bound on the vertical component. This phenomenon suggests that although the QH method has certain advantages in horizontal error statistics, it may be more sensitive to certain geometric conditions in the vertical direction, thus introducing more significant vertical error dispersion. In summary, NVSLSM, while maintaining a consistent accuracy structure with the segment bisection method, can achieve better computational performance through improved iterative efficiency; while the QH method exhibits certain differences in error behavior across different directional components.

[0102] This example analyzes 2 seconds of data from a trajectory, generating a 21*101 DDM using convolution, FFT, and FFT deployed on a GPU. The example provides a comparison of computation time for the same trajectory. It shows that temporal convolution takes 876.70 seconds, FFT on a CPU platform takes 143.00 seconds, while GPU-accelerated FFT takes only 12.64 seconds. While maintaining consistency in the generated results, FFT achieves approximately 6.1 times the speedup compared to direct convolution, while GPU-FFT further improves upon CPU-FFT by approximately 11.3 times, resulting in an overall speedup of nearly 70 times compared to temporal convolution. The fundamental difference lies in the combined impact of algorithm complexity and hardware parallelism. The computational complexity of temporal convolution increases quadratically with data length, while FFT can transform convolution operations into frequency domain multiplication, reducing the overall complexity to... Its advantages are particularly significant in medium- and large-scale data processing. In addition, the DDM generation process involves a large number of batch operations with consistent structure, which is suitable for the parallel computing mode of GPUs. Therefore, GPU-FFT can further improve the algorithmic complexity advantage and achieve an order-of-magnitude performance improvement.

[0103] Figure 7Comparisons of DDM results generated using three implementations—convolution, CPU-FFT, and GPU-FFT—are presented. From the perspective of overall structure and energy distribution characteristics, the three methods maintain a high degree of consistency in the main peak position, delay broadening morphology, and Doppler expansion characteristics, indicating that frequency domain features and temporal domain convolution are theoretically equivalent, and GPU acceleration does not alter the physical structure and amplitude distribution characteristics of the DDM. It should be noted that in the GPU implementation, to meet the data length consistency requirements of the parallel computing framework, the input sequence must be processed to a uniform length. Therefore, there are instances of truncation or padding of individual sampling points in some millisecond data segments; that is, different millisecond segments may have one extra or one missing point. This operation is an engineering-level data alignment process and does not change the main lobe structure or energy concentration region of the DDM, only potentially introducing very weak differences at the edge delay units. Further quantitative analysis shows that, using the DDM generated by convolution as a reference, the correlation coefficients of the DDM generated by CPU-FFT and GPU-FFT are 1 and 0.99, respectively, indicating that the three methods maintain strict consistency at the numerical level, with only minor differences within the machine precision range. This result verifies the equivalence between the frequency domain implementation and the time domain convolution in terms of implementation, and also shows that no additional numerical bias is introduced during GPU parallelization.

[0104] This example illustrates the processing time variation corresponding to the reference signal length for a single migration to the GPU, assuming a fixed total processing time of 2 seconds. It's important to note that the overall data volume processed in each group of experiments remained consistent; the only difference lay in the size of the data block transferred to the GPU each time. It can be observed that when the data block size increased from 1 ms to 20 ms, the processing time decreased from 12.64 s to 10.12 s, indicating that smaller-granularity data migration introduces higher call and transmission overhead. As the data block size further increased to 50 ms and above, the processing time stabilized at approximately 8–9 s, indicating that GPU computing resources were now fully utilized, and the proportion of data transfer and kernel function scheduling overhead in the overall time significantly decreased. In other words, once the block size exceeds a certain threshold, system performance shifts from being computationally driven to being driven by bandwidth or hardware limitations, and the benefits of further increasing the block length gradually diminish. This result demonstrates that, under a fixed total processing volume, the performance of the GPU implementation is not only affected by algorithm complexity but also closely related to the data block size. Increasing the size of data blocks transmitted in a single transaction can effectively reduce the overhead of frequent memory copying and kernel function startup, thereby improving overall throughput efficiency.

[0105] The core computational modules in the GNSS-R DDM generation process mainly include FFT / IFFT, complex conjugate operations, point-to-point complex multiplication, and cumulative averaging. Among them, FFT and IFFT usually account for the main computational overhead and are key links affecting the overall processing efficiency. To systematically evaluate the performance improvement effect of GPU-accelerated FFT in actual processing, this invention constructs two implementation schemes: (1) a CPU software receiver based on Python 3.7; (2) a GPU-accelerated software receiver based on PyCUDA.

[0106] Taking the processing of 1 second of data as an example, this corresponds to a millisecond-by-millisecond computation process of approximately 1000 complex floating-point samples, equivalent to processing 1 second of CYGNSSGPS L1 C / A data. The time consumption of each major computation module is statistically analyzed. From... Figure 8 As can be seen, in the CPU implementation, FFT and IFFT constitute the main sources of time overhead, with the time consumption of a single module being much higher than that of complex multiplication and averaging operations, indicating that frequency domain transformation dominates the overall computation chain. In contrast, in the GPU implementation, the time consumption of FFT and IFFT is significantly reduced, and the time share of each computation module is more balanced, indicating that the parallel architecture effectively alleviates the computational bottleneck caused by frequency domain transformation. Further analysis shows that FFT-like operations have highly regular data access patterns and good parallel decomposability, making them suitable for large-scale parallel computation on GPU architecture; while point-to-point complex multiplication and averaging operations themselves have low computational complexity, and therefore do not form a major bottleneck on the CPU platform. The GPU implementation achieves significant time compression by uniformly mapping frequency domain transformation and batch multiplication operations to a parallel computing framework. The above results indicate that in the GNSS-R DDM generation task, the key to performance optimization lies in the parallel implementation of the frequency domain transformation module. GPU-FFT not only has advantages in theoretical complexity but also shows significant acceleration effects at actual processing scales, providing a computationally feasible foundation for subsequent batch processing and near real-time applications of orbital data.

[0107] To evaluate the effectiveness of the proposed method in sea surface altimetry, this example uses the sea surface height after tidal correction overlaid on the DTU sea surface model as a reference to conduct consistency tests and error statistical analyses on the Half and Der altimetry methods. Figure 9 a and Figure 9 Figure c shows the scatter density distribution between the estimated altitude and the reference altitude. It can be observed that the scatter points of both methods are distributed along the diagonal, with high-density areas concentrated near the diagonal. This indicates that the estimation results can statistically characterize the dominant trend of sea surface height changes. At the same time, the point cloud still has some dispersion on both sides of the diagonal and is accompanied by a small number of outliers, reflecting that the altimeter error may increase significantly under conditions such as complex sea conditions, low glancing angles, or fluctuations in signal-to-noise ratio.

[0108] Figure 9 b and Figure 9 The residual characteristics are further presented in section d. The residual histograms of both methods show a peak-like distribution centered near zero, but also exhibit some heavy-tailed characteristics, indicating that the errors are not strictly Gaussian and that a small number of large error samples exist. Based on full-sample statistics, the mean absolute error (MAE) of HALF is 16.8 m, with a correlation coefficient of 0.703 with the reference height; the MAE of Der is 15.9 m, with a correlation coefficient of 0.709. Compared to HALF, Der reduces the MAE by approximately 5.2% and slightly improves the correlation coefficient, indicating that Der achieves a small improvement in the consistency between the overall error magnitude and the trend of height change. Combining the scatter distribution and residual statistics, it can be seen that both methods can effectively recover the main characteristics of sea surface height change, but significant errors may still occur under some unfavorable observation conditions. Overall, Der achieves a lower MAE while maintaining a similar correlation to HALF, verifying its effectiveness and relative advantage in this dataset.

[0109] This invention addresses two key bottlenecks in spaceborne GNSS-R intermediate frequency data processing: efficient solution for specular reflection points and DDM generation of large-scale intermediate frequency data. Experimental results show that, while maintaining project consistency, a balance between accuracy and efficiency can be achieved through algorithm structure optimization and computational architecture reconstruction.

[0110] In solving for specular reflection points, the differences in spatial error distribution among different algorithms reflect the inherent asymmetry of GNSS-R observation geometry. The horizontal component is mainly constrained by orbital geometry, and its condition number is relatively stable; however, the vertical component is more sensitive to altitude model and initial value perturbations, thus more easily amplifying nonlinear iterative errors. The improved method exhibits more stable convergence characteristics in the three-dimensional direction, indicating that reasonable step size control and convergence criterion design can effectively suppress the spread of vertical errors. This phenomenon shows that under high-dynamic conditions on spaceborne systems, the geometrically constrained enhanced solution framework is superior to strategies that rely solely on local iterative convergence. In the DDM generation stage, the FFT implementation and the convolutional form are theoretically equivalent, and the correlation coefficient reached 1 in the experiment, verifying the strict numerical consistency of the frequency domain implementation. Therefore, the performance difference stems entirely from the computational architecture, not the algorithm itself; GPU parallelization significantly reduces the processing cost per unit time.

Claims

1. A GNSS-R sea surface altimetry method based on variable step-size line segment intermediate frequency data processing, characterized in that, This includes acquiring raw intermediate frequency data from spaceborne GNSS-R, calculating specular reflection points using the novel variable step-size line segment NVSLSM method, generating a delayed Doppler image (DDM), and performing sea surface height inversion. Specifically, the calculation of specular reflection points using the novel variable step-size line segment NVSLSM method includes the following steps: Step 1-1: Input the receiver and transmitter positions. The position information is obtained from the navigation information. Let: （1）， in, It is the first The difference between the incident angle and the reflection angle at the next iteration. and The calculation formula is as follows: （2）， （3）； Step 1-2: Initialize the specular reflection points using formula (2): （4）， （5）， （6）， In the WGS84 coordinate system, the major axis The eccentricity is 6,378,137 m. It is 0.08181919084264; Steps 1-3: Calculate the initial difference in the number of mirror reflection points: （7）， （8）， Wherein, the initial iteration step size and They are 0.4 and 0.6 respectively. Steps 1-4: Calculate the step size for subsequent iterations: （9）； Steps 1-5: Recalculate the difference between the angle of incidence and the angle of reflection: （10）， （11）； Steps 1-6: If conditions If this condition is met, then the iteration exit condition is satisfied. If it is the final specular reflection point, otherwise return to steps 1-4 and iterate again. for .

2. The GNSS-R sea surface altimetry method based on variable step-size line segment intermediate frequency data processing according to claim 1, characterized in that, The generation of delayed Doppler maps (DDMs) is achieved by constructing an FFT-based GPU parallel processing architecture, which includes: Step 2-1: Signal generation and open-loop tracking. The initial target motion parameters are provided by the open-loop tracking module. The signal generator generates a transmission signal based on the open-loop tracking information and splits the signal into multiple processing branches: one branch is sent to the reference signal replication module, one branch is sent to the pseudo-random code modulation module C / A Code, and one branch directly participates in the mixing. Step 2-2: Receive the echo signal, i.e. the reflected signal, from the target through two or more receiving channels. Each reflected signal is mixed with the corresponding local replica signal, i.e., multiplied, in a mixer. Steps 2-3: C / A code modulation and matched filtering. Two or more parallel C / A code modules receive the raw signal from the signal generator, apply C / A pseudo-random code modulation, and the modulated signal also enters the mixer. It is multiplied with the reflected signal to realize code correlation processing, and the output signal is processed by FFT. Steps 2-4: FFT processing is implemented through at least two parallel FFT modules. The FFT output then enters the IFFT module, which has at least two channels. The outputs of the two or more IFFT channels are combined and then enter the incoherent accumulation module. Steps 2-5: The incoherent accumulation module fuses the multi-channel IFFT output signals, identifies significant target points, and outputs them as DDM distance-Doppler maps.

3. The GNSS-R sea surface altimetry method based on variable step-size line segment intermediate frequency data processing according to claim 2, characterized in that, Sea surface height inversion requires re-tracking the position of the specular reflection point on the waveform to extract the time delay value while eliminating the path delay error caused by sea surface roughness. Based on the Der method, the reflected signal corresponding to the specular reflection point is the part that contributes the most to the power waveform, that is, the time delay corresponding to the largest first derivative of the power waveform, as shown in the following formula: （13）， （14）， The path difference between the direct signal and the reflected signal is expressed as: （15）， in, For the observed bistatic delay, and The direct and reflected signal paths are identified. Based on the positions of the receiver and transmitter, the coordinates of the specular reflection point on the WGS84 ellipsoid are calculated. After error correction, the bistatic delay in the model is obtained. The calculation is as follows: （16）， in, The geometric model delay between the direct and reflected signals is used to obtain the height of the specular reflection point relative to the WGS84 ellipsoid: (17), among which, It is the angle of incidence of the reflected signal at the point of reflection on the mirror.

4. The GNSS-R sea surface altimetry method based on variable step-size line segment intermediate frequency data processing according to claim 3, characterized in that, Mean sea level height data is derived from the DTU21 model, and tidal correction uses the TPXO9 model. The two models are superimposed to form a dynamically changing sea level height reference field, which is used to verify the accuracy and reliability of the inversion results. Represented as (18), of which, express The actual sea level height of the model after tidal correction. This represents the average sea level height calculated by the model. This indicates the tidal correction calculated by the model.

5. The GNSS-R sea surface altimetry method based on variable step-size line segment intermediate frequency data processing according to claim 4, characterized in that, To ensure the accuracy and reliability of the sea surface height inversion results and to improve the performance of the inversion model, the following steps were implemented: (1) land reflection point data were removed, and only observations of specular reflection points located in the ocean area were retained; (2) power waveforms with a signal-to-noise ratio of less than 10 dB were excluded to reduce the impact of low-quality data on the inversion results. (19), among which, It is the maximum power of the power waveform. It is the average power of the power waveform, which is taken from the zero-Doppler slice region in the Delayed Doppler Map (DDM).

Images