A method for detecting signal loss of a beidou wave buoy and repairing spectral consistency of a pseudo signal

Abstract

The present application relates to the technical field of ocean observation, and discloses a Beidou wave buoy signal loss-of-lock artifact detection and spectrum consistency repair method, which comprises the following steps: extracting an original measurement data sequence, generating a synchronous state data matrix and separating a three-dimensional detrended displacement sequence; constructing a signal quality marker sequence and performing single-wave division, calculating kinematic characteristic parameters based on wave height comparison results; comparing the kinematic characteristic parameters with state markers according to dynamic limit limits, and intercepting a positioning artifact wave band; solving autoregressive model coefficients, implementing bidirectional prediction weighted fusion, replacing abnormal data to generate a time-domain reconstructed vertical sequence; obtaining an energy attenuation slope test theoretical deviation of the time-domain reconstructed vertical sequence, calculating a spectrum check factor mean square error in combination with a horizontal component, and outputting parameters externally when the error is lower than a specified upper limit. The present application combines dynamic physical limits to investigate artifact causes, and ensures the integrity and reliability of wave data through bidirectional deduction and frequency domain double verification.

Description

Technical Field

[0001] This invention relates to the field of marine observation technology, specifically a method for detecting artifacts of signal loss and restoring spectral consistency of BeiDou wave buoys. Background Technology

[0002] Ocean wave observation is fundamental to marine scientific research and marine resource development. Traditional wave observation often employs large buoys based on accelerometers, measuring and integrating motional acceleration to obtain displacement. These devices do not rely on external signals, and data acquisition is relatively stable. With the development of satellite navigation technology, small drifting buoys based on BeiDou positioning have become an important means of acquiring sea state data due to their ease of deployment and wide coverage.

[0003] Unlike accelerometer buoys, BeiDou wave buoys obtain their position coordinates directly by receiving satellite signals. When the buoy encounters wave impacts or momentarily submerges, the antenna entering the water causes a brief loss of BeiDou signal. The instant the signal recovers, a jump in positioning coordinates occurs. This jump data, after bandpass filtering, creates a positioning artifact in the displacement time series that has a large amplitude and shape similar to actual waves. This artifact is deceptive and causes the final output wave statistical parameters to deviate from the true values.

[0004] Currently, the main methods for dealing with positioning artifacts like those from BeiDou buoys rely on direct data removal or statistical threshold determination. Some conventional solutions simply remove artifacts during abnormal periods based on state parameters such as the number of visible satellites or position accuracy factors, and then use linear interpolation to fill in the gaps. This approach leads to a significant loss of effective observational data in sea conditions with frequent signal interruptions, affecting the temporal continuity of wave sequences and reducing data integrity.

[0005] Meanwhile, some solutions directly adopt the quality control standards of accelerometer buoys, identifying abnormal data by setting fixed wave height statistical thresholds. Since the amplitude of positioning artifacts depends on the duration of lock-on loss and the degree of positional jumps, and has no deterministic correlation with background sea conditions, fixed statistical thresholds are insufficient to distinguish positioning artifacts from normal transient waves. This method, relying on a single statistical feature, fails to consider the specific causes of BeiDou artifacts, leading to missed detections and the accidental deletion of normal wave data.

[0006] Furthermore, existing processing methods for anomalous wavebands mostly remain at the level of marking or discarding, lacking verification methods and restoration mechanisms that incorporate wave hydrodynamic boundaries. After filling in anomalous data, current methods do not establish frequency domain verification criteria, making it impossible to verify whether the reconstructed wave surface sequence conforms to the physical laws of deep-water wave evolution in terms of frequency domain characteristics. This restoration approach, lacking physical consistency constraints, leads to unreliable wave spectrum analysis results in the final output. Summary of the Invention

[0007] To address the shortcomings of existing technologies, this invention provides a method for detecting artifacts caused by the loss of lock on BeiDou wave buoy signals and for restoring spectral consistency. This method solves the problem that conventional data processing methods cannot accurately identify artifact intervals when existing wave measurement data is distorted due to short-term loss of lock on BeiDou signals or external interference, and that the restored wave surface sequence deviates from the laws of ocean dynamics, resulting in distortion of wave spectrum characteristics.

[0008] To achieve the above objectives, the present invention provides a method for detecting artifacts in BeiDou wave buoy signal lock-up and repairing spectral consistency, comprising the following steps:

[0009] Read the observation messages and extract the raw measurement data sequence, align the timestamps to generate a synchronization status data matrix;

[0010] Based on the synchronization state data matrix, local tangent plane transformation and high-pass filtering are performed to separate the three-dimensional detrended displacement sequence.

[0011] A signal quality marker sequence is constructed based on a three-dimensional detrended displacement sequence and single-wave segmentation is performed. Then, the single-wave height is compared with the background significant wave height, and the kinematic characteristic parameters are calculated based on the comparison results.

[0012] Based on the kinematic characteristic parameters and the corresponding state markers, the dynamic limit limits are compared and logical judgments are performed. The positioning artifact bands are extracted, the autoregressive model coefficients are calculated, and bidirectional prediction weighted fusion is implemented to replace the original abnormal data in order to repair the signal lock artifacts and generate a time-domain reconstructed vertical sequence.

[0013] The energy decay slope is obtained by performing a fast Fourier transform on the reconstructed vertical sequence in the time domain to verify the theoretical deviation. The mean square error of the spectrum verification factor is calculated by extracting the horizontal component in the three-dimensional detrended displacement sequence and the reconstructed vertical sequence in the time domain. When the error is lower than the specified upper limit, high-quality ocean wave parameters are output.

[0014] Furthermore, in the process of reading observation reports and extracting the original measurement data sequence, aligning the timestamps, and separating the three-dimensional detrended displacement sequence: if and only if the difference between adjacent timestamps is strictly greater than zero, a linear interpolation algorithm is used to map and align the extracted original measurement data sequence to the corresponding sampling time node, and the data is merged to construct a synchronous state data matrix; the geodetic coordinate system parameters in the synchronous state data matrix are retrieved for local tangent plane transformation to separate the northward displacement sequence and the eastward displacement sequence, and the vertical elevation detrending formula is used to subtract the arithmetic mean of the ellipsoidal height within the observation time window from the current ellipsoidal height to extract the original vertical displacement component; high-pass filtering is applied to the northward displacement sequence, the eastward displacement sequence, and the original vertical displacement component to eliminate extremely low frequency drift, and a three-dimensional detrended displacement sequence containing the detrended northward displacement sequence, the detrended eastward displacement sequence, and the detrended vertical displacement sequence is separated.

[0015] Furthermore, the steps of constructing a signal quality indicator sequence based on the three-dimensional detrended displacement sequence and performing single-wave segmentation specifically include: extracting the state parameters under the corresponding timestamp and performing comprehensive evaluation using multivariate logical OR operations to generate a signal quality indicator sequence; obtaining the vertical components in the three-dimensional detrended displacement sequence, and forcibly skipping zero-crossing interpolation when the absolute value of the values ​​of consecutive sampling points is less than the set minimum floating-point number limit, using adjacent zero-crossing moments as time fences for single-wave segmentation to obtain independent wave period segments; and projecting and combining the time periods corresponding to the independent wave period segments onto the signal quality indicator sequence to obtain a set of independent wave units.

[0016] Furthermore, the step of comparing the single wave height with the significant background wave height specifically includes: constructing a statistical time window by extending forward and backward in both directions from the current analyzed wave in the independent wave unit set as the center; calculating the variance of the wave surface displacement within the statistical time window to derive the macroscopic background benchmark and obtain the significant background wave height; determining the extreme span of the wave band for the wave unit in the independent wave unit set to obtain the single wave height; comparing the single wave height with the significant background wave height; and performing extraction when the single wave height is greater than the product of the significant wave height multiple threshold and the significant background wave height to extract a set of suspected abnormal candidate waves; the significant wave height multiple threshold is the dimensionless multiplier ratio limit for judging the deviation of the transient change wave height from the sea state background energy level.

[0017] Furthermore, the steps for calculating kinematic characteristic parameters based on the comparison results specifically include: based on the comparison results, for displacement sampling points within the suspected abnormal candidate wave set, the vertical acceleration is obtained by using the five-point center difference formula for vertical acceleration to obtain the transient derivative; based on the principles of ocean dynamics, the wavelength spatial scale is obtained by using the deep-water dispersion relation wavelength calculation formula, and the wave steepness is calculated by using the wave steepness calculation formula in conjunction with the single-wave height; the vertical acceleration and wave steepness are combined to extract the kinematic characteristic parameters; the kinematic characteristic parameters are checked for limit limits, and a physical anomaly triggering state is generated when the degree of morphological distortion exceeds the dynamic limit limit.

[0018] Furthermore, in the process of extracting the positioning artifact band and calculating the autoregressive model coefficients through logical judgment: multi-dimensional logical comprehensive judgment is performed by combining the physical anomaly trigger state and the state identifiers corresponding to the suspected anomaly candidate wave set to extract the positioning artifact band; the signal missing span of the positioning artifact band is measured, and the outer normal observation sampling points are extracted to form a sample sequence; the variance of the sample sequence is pre-calculated, and when the variance is lower than the set machine floating-point minimum limit, it is determined that the covariance matrix is ​​approaching degradation. The bypass autoregressive calculation process is downgraded to use low-order polynomial smooth interpolation to fill the gaps. When the covariance matrix is ​​not degraded, the autoregressive model coefficients are calculated.

[0019] Furthermore, in the process of implementing bidirectional prediction weighted fusion to replace the original abnormal data to repair signal lock-up artifacts and generate a time-domain reconstructed vertical sequence: the autoregressive model coefficients are called, and the normal data on the left side of the outer side is used as the benchmark to extrapolate to the right. Similarly, the normal data on the right side of the outer side is used as the benchmark to extrapolate to the left, thus obtaining the positive and negative prediction values. The missing span in the positioning artifact band is extracted, and an inverse weight is assigned based on the time interval between each target point to be estimated within the missing span and the boundary of the left and right real data. The positive and negative prediction values ​​are then weighted and fused to replace the original abnormal data, generating a time-domain reconstructed vertical sequence.

[0020] Furthermore, in the process of verifying the theoretical deviation by performing a Fast Fourier Transform (FFT) on the reconstructed vertical sequence in the time domain to obtain the energy attenuation slope: the FFT is used to optimize and locate the spectral peak frequency corresponding to the main energy peak, and the spectral peak frequency is extended and truncated towards the high-frequency direction as a set multiple as the starting point to extract the high-frequency balanced frequency band data; the high-frequency balanced frequency band data is projected onto a double logarithmic coordinate system and a univariate linear regression is performed to obtain the energy attenuation slope, and the spectral slope deviation between the energy attenuation slope and the theoretical value of the attenuation slope is calculated; when the spectral slope deviation is less than the set spectral slope deviation threshold, the reconstructed vertical sequence in the time domain is determined to meet the theoretical deviation test; the theoretical value of the attenuation slope is a pre-set ideal exponential limit characterizing the high-frequency energy attenuation law of deep-water gravity waves with frequency; the spectral slope deviation threshold is the allowable floating-point error limit used to determine whether the high-frequency energy attenuation characteristics of the reconstructed wavefront deviate from the ideal state.

[0021] Furthermore, in the process of calculating the mean square error of the spectral verification factor and outputting parameters: the detrended eastward displacement sequence and the detrended northward displacement sequence in the three-dimensional detrended displacement sequence are extracted, the power spectral density is calculated, and linear equal-weight superposition is performed to construct the horizontal total displacement power spectrum; the vertical hydrodynamic power spectral density is derived from the time-domain reconstructed vertical sequence, and the vertical hydrodynamic power spectral density is combined with the horizontal total displacement power spectrum to calculate the spectral domain physical consistency verification factor; the core frequency band containing the main wave energy is defined to evaluate the global error and obtain the mean square error of the spectral verification factor; when the mean square error of the spectral verification factor is lower than the specified upper limit, high-quality ocean wave parameters are output; the specified upper limit is the maximum tolerance limit used to characterize the degree of agreement between the two-dimensional spectral energy of horizontal dynamics and vertical kinematics.

[0022] This invention also provides a system for detecting artifacts of signal loss and restoring spectral consistency of BeiDou wave buoy signals, comprising:

[0023] The data acquisition and preprocessing module 10 is used to parse the original observation report and perform coordinate system transformation and detrending filtering to generate a three-dimensional detrending displacement sequence;

[0024] The anomaly detection and feature extraction module 20 establishes a data connection with the data acquisition and preprocessing module 10 through a data link. The anomaly detection and feature extraction module 20 is used to receive the three-dimensional detrended displacement sequence and segment it into single waves, then compare the wave height of the single wave with the background significant wave height, and output a set of suspected anomaly candidate waves and corresponding kinematic feature parameters based on the comparison results.

[0025] The decision and reconstruction module 30 establishes a data connection with the anomaly detection and feature extraction module 20 through a data link. The decision and reconstruction module 30 is used to perform joint condition judgment and interpolation substitution based on kinematic feature parameters to generate a temporal reconstruction vertical sequence.

[0026] The verification and output module 40 establishes a data connection with the decision and reconstruction module 30 through a data link. The verification and output module 40 is used to convert the time-domain reconstructed vertical sequence to the frequency domain to verify the physical laws and output high-quality ocean wave parameters.

[0027] This invention provides a method for detecting artifacts in BeiDou wave buoy signal lock-up and repairing spectral consistency. It has the following beneficial effects:

[0028] 1. This invention extracts state parameters to generate a signal quality indicator sequence, and performs multi-dimensional logical comprehensive judgment by combining the kinematic feature parameters composed of vertical acceleration and wave steepness. This step combines the underlying signal quality with the physical limits of wave dynamics to directly locate the cause of artifacts, avoiding the missed detection and misjudgment caused by the traditional single statistical threshold method due to the failure to distinguish the source of data distortion, and improving the accuracy of abnormal artifact detection in Beidou wave measurement data.

[0029] 2. This invention replaces the original abnormal data by performing bidirectional prediction weighted fusion through solving the autoregressive model coefficients, reducing the damage to measurement integrity caused by direct data removal. On this basis, it performs two-dimensional consistency verification by calculating the spectral slope deviation and the mean square error of the spectral verification factor, so that the repaired wave surface sequence remains continuous in the time domain and conforms to the high-frequency attenuation law of waves and the three-dimensional displacement dynamics relationship in the frequency domain. This solves the problem of wave spectrum feature distortion caused by conventional interpolation methods deviating from physical laws.

[0030] 3. This invention can achieve artifact detection and spectral consistency repair by relying only on the original measurement data sequence output by the wave buoy and its built-in state parameters. No additional auxiliary sensors are required. The method requires a single data source and has a moderate computational load. It can be directly deployed on a low-power embedded hardware platform to run in real time, which reduces the manufacturing cost of marine monitoring equipment and improves the engineering applicability of this method. Attached Figure Description

[0031] Figure 1 This is a system architecture diagram of a Beidou wave buoy signal loss artifact detection and spectral consistency repair system according to an embodiment of the present invention;

[0032] Figure 2 This is a flowchart illustrating the overall steps of a method for detecting artifacts in BeiDou wave buoy signal lock-up and repairing spectral consistency, provided in an embodiment of the present invention.

[0033] Figure 3 The above are comparison diagrams of the temporal domain detection and repair effects of typical positioning artifacts in the embodiments of the present invention, wherein (a) is a wavefront displacement time sequence diagram and (b) is a Beidou signal quality marker sequence time sequence diagram.

[0034] Figure 4 This is a comparison chart of the spectral morphology consistency test results before and after repair in an embodiment of the present invention;

[0035] Figure 5 The above is a time-series comparison of significant wave heights between the Beidou buoy and the reference buoy in an embodiment of the present invention, wherein (a) is a panoramic comparison of continuous observations and (b) is a magnified comparison of the local time period in which the positioning artifact is located.

[0036] Figure 6 The above are comparative verification diagrams of the detection method for small-amplitude positioning artifacts in the embodiments of the present invention, wherein (a) is a wavefront displacement time sequence comparison diagram and (b) is a signal quality indicator sequence diagram corresponding to short-time signal loss.

[0037] The module consists of: 10. Data acquisition and preprocessing module; 20. Anomaly detection and feature extraction module; 30. Decision and reconstruction module; and 40. Verification and output module. Detailed Implementation

[0038] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0039] See attached document Figure 1 The system provided in this embodiment is mounted on a BeiDou wave buoy hardware platform. The hardware platform includes a BeiDou satellite navigation receiver, a satellite antenna, and a data processing unit. The BeiDou satellite navigation receiver collects and outputs raw observation messages according to a set frequency. The data processing unit establishes a data connection with the BeiDou satellite navigation receiver, receives the raw observation messages, and allocates computing resources to perform subsequent computational processing. The BeiDou wave buoy signal lock-off artifact detection and spectral consistency repair system provided by this invention includes:

[0040] The data acquisition and preprocessing module 10 is used to parse the original observation report and perform coordinate system transformation and detrending filtering to generate a three-dimensional detrending displacement sequence;

[0041] Anomaly detection and feature extraction module 20, connected to data acquisition and preprocessing module 10, is used to receive three-dimensional detrended displacement sequence and segment single waves, and output a set of suspected anomaly candidate waves and corresponding kinematic feature parameters.

[0042] The decision and reconstruction module 30, connected to the anomaly detection and feature extraction module 20, is used to perform joint conditional judgment and interpolation substitution based on kinematic feature parameters to generate a temporal reconstructed vertical sequence.

[0043] The verification and output module 40, connected to the decision and reconstruction module 30, is used to convert the time-domain reconstructed vertical sequence to the frequency-domain reconstructed physical laws and output high-quality ocean wave parameters.

[0044] See attached document Figure 2 , Figure 2 This is a flowchart of a method for detecting and repairing spectral consistency of BeiDou wave buoy signal lock-up artifacts according to an embodiment of the present invention. The present invention provides a method for detecting and repairing spectral consistency of BeiDou wave buoy signal lock-up artifacts, comprising the following steps:

[0045] Step S1: System initialization and parameter configuration, setting the data sampling frequency of the Beidou satellite navigation receiver and completing the system operation limit parameter configuration, generating a decision threshold parameter set;

[0046] Step S2: Multi-source raw data stream acquisition. Based on the configuration instructions in the decision threshold parameter set, read the Beidou satellite navigation receiver message through the communication interface and extract the raw measurement data sequence containing position and status information.

[0047] Step S3: Time reference alignment and synchronization. Interpolation is performed on the original measurement data sequence based on the hardware second pulse signal to align the timestamps and generate a synchronization status data matrix.

[0048] Step S4: Coordinate transformation and detrending. Extract the position parameters from the synchronous state data matrix, perform local tangent plane transformation and high-pass filtering, and separate the three-dimensional detrending displacement sequence.

[0049] Step S5: Signal quality anomaly pre-labeling, extract the state parameters under the corresponding timestamp of the three-dimensional detrended displacement sequence, compare them with the state boundary values ​​in the decision threshold parameter set, and construct a signal quality label sequence;

[0050] Step S6: Wave unit segmentation, obtain the vertical component in the three-dimensional detrended displacement sequence, use the zero-overpass method to search for the intersection of the wave surface curve and the zero mean line to perform single wave segmentation, and combine the corresponding interval range of the signal quality indicator sequence to obtain an independent wave unit set.

[0051] Step S7: Dynamic estimation of significant wave height. Within the sliding time window, the zero-order spectral moment of the set of independent wave units is calculated to estimate the significant wave height in the background. By comparing the single wave height with the significant wave height in the background, a set of suspected anomalous candidate waves is extracted.

[0052] Step S8: Kinematic feature calculation. The displacement sampling points in the suspected abnormal candidate wave set are numerically differentiated using the central difference scheme. Based on the deep-water dispersion relation, the kinematic feature parameters including vertical acceleration and wave steepness are calculated.

[0053] Step S9: Physical limit violation judgment. Compare the kinematic feature parameters with the dynamic limit limits in the judgment threshold parameter set. When the dynamic limit limits are exceeded, a physical anomaly trigger state is generated.

[0054] Step S10: Multi-dimensional feature fusion judgment, combining the physical anomaly trigger state and the state identifier corresponding to the suspected anomaly candidate wave set to perform logical judgment, and extract the positioning artifact wave band;

[0055] Step S11: Based on the autoregressive model, time series reconstruction is performed by measuring the signal missing span of the localization artifact band, extracting the outer normal observation sampling points, calculating the autoregressive model coefficients, and implementing bidirectional prediction weighted fusion to replace the original abnormal data, thereby generating a time-domain reconstructed vertical sequence.

[0056] Step S12: Spectral morphology consistency check. Perform fast Fourier transform on the reconstructed vertical sequence in the time domain to obtain the power spectral density distribution. Perform univariate linear regression fitting on the high-frequency balanced band in a double logarithmic coordinate system to obtain the energy attenuation slope.

[0057] Step S13: Spectral domain physical consistency closed-loop verification, check the theoretical deviation of energy decay slope, extract the horizontal component in the three-dimensional detrended displacement sequence and calculate the mean square error of the spectral verification factor with the time domain reconstructed vertical sequence. When the error is lower than the specified upper limit, output high-quality ocean wave parameters.

[0058] To further clarify the internal operation principles of the data acquisition and preprocessing module 10, the anomaly detection and feature extraction module 20, the decision and reconstruction module 30, and the verification and output module 40, as well as the deep data transfer relationships between the above steps, the following will provide a detailed explanation of the mathematical model and implementation details of each specific step.

[0059] In this embodiment, to establish a reliable data processing benchmark, the data acquisition and preprocessing module 10 executes step S1: system initialization and parameter configuration. This step aims to establish hardware operating boundaries and algorithm decision benchmarks for the entire anomaly detection and repair process, and is specifically detailed as follows:

[0060] Sub-step S101: Set the data sampling frequency of the BeiDou satellite navigation receiver. Based on the Nyquist sampling theorem and the high-frequency component distribution characteristics of typical ocean gravity waves, in order to fully capture the real physical details of sea surface waves and meet the discrete point requirements of the subsequent time-series reconstruction algorithm, the system sampling frequency is set to be no less than 2Hz.

[0061] Sub-step S102: Initialize the parameter set used for multi-dimensional feature decision-making. Based on preset hardware performance and ocean dynamics experience, construct a decision threshold parameter set expressed by the decision threshold parameter set formula, which is:

[0062] ;

[0063] in, Represents the set of decision threshold parameters; The threshold value for significant wave height is preferably between 1.5 and 2.0. This parameter is determined based on the Rayleigh distribution statistical model and is used to initially screen out large waves or abrupt artifacts from sea surface background fluctuations. This represents the upper limit threshold of the ratio of acceleration to gravitational acceleration, typically around 1.2. Its physical basis is the limit kinematic constraints of actual wave breaking, used to intercept sharp, non-physical jumps. This represents the minimum number of visible satellites required for a valid positioning. Based on the fundamental principle of three-dimensional spatial geometric calculation, this value is fixed to be no less than 4. This indicates the upper limit threshold of the position accuracy factor, which is set to 5.0 based on the general satellite navigation integrity standard; This indicates the maximum allowed continuous interruption duration of the signal. Its value depends on the short-term predictive divergence characteristics of the selected autoregressive model. As a preferred approach, it is set to 1.5 seconds.

[0064] After configuring the system's operational limit parameters, the system enters real-time monitoring and acquisition mode. The data acquisition and preprocessing module 10 executes step S2: multi-source raw data stream acquisition. By capturing spatial location changes and underlying communication status, a dual-channel data source is formed to support subsequent fusion decisions. This is specifically broken down into the following sub-steps:

[0065] Sub-step S201: Extract the spatial sequence containing position information. The observation message output by the BeiDou satellite navigation receiver is read through the buoy's internal communication interface. The three-dimensional spatial coordinates are extracted using the position measurement vector formula, thereby establishing the buoy's trajectory in absolute space. The position measurement vector formula is:

[0066] ;

[0067] in, Indicates time Position measurement vector; Indicates time Latitude; Indicates time longitude; Indicates time The height of the ellipsoid; superscript This represents the matrix transpose operation.

[0068] Sub-step S202: Extract the signal sequence containing state information. Simultaneously parse the communication link state identifier in the same observation message, and extract the underlying integrity parameters using the signal quality measurement vector formula. The technical purpose of extracting this set of parameters is to quantitatively assess the reliability of the instantaneous satellite communication link, providing a fundamental basis for distinguishing between physical waves and electronic signal jump artifacts. The signal quality measurement vector formula is:

[0069] ;

[0070] in, Indicates time The signal quality measurement vector; Indicates time The number of visible satellites; Indicates time The average signal-to-noise ratio of each satellite carrier; Indicates time The position precision factor. The extracted time mentioned above. Position measurement vector and time The signal quality measurement vectors together constitute the original measurement data sequence required by the system.

[0071] Due to the difference in the output update rate between the positioning calculation and the status message at the hardware level, the original measurement data sequence has a time misalignment. To eliminate this deviation, the data acquisition and preprocessing module 10 performs step S3: time base alignment and synchronization. This is specifically broken down into the following sub-steps:

[0072] Sub-step S301: Hardware clock alignment. Receive the hardware second pulse signal output from the BeiDou satellite navigation receiver, using the trigger edge of this second pulse signal as the absolute time reference to eliminate physical bus delay during heterogeneous data transmission within the system.

[0073] Sub-step S302: Heterogeneous data timestamp synchronization. To prevent errors caused by extremely small time differences due to packet loss, the system pre-verifies the time interval between adjacent status messages before performing resampling. A linear interpolation algorithm is used to synchronize the timestamps only if the difference between adjacent timestamps is strictly greater than zero. Signal quality measurement vector mapping aligned to time The corresponding sampling time nodes of the position measurement vectors. The data sequences after merging and aligning are constructed using the synchronization state matrix construction formula, which is:

[0074] ;

[0075] in, Indicates the first Synchronization state matrix at each sampling time; Numeric index subscripts representing discrete sampling points; Indicates the first The absolute timestamp of each sampling moment; Indicates the first Position measurement vector at each sampling time; Indicates the first Signal quality measurement vectors at each sampling time point. By constructing this matrix, a synchronization state data matrix with strictly consistent time nodes is generated.

[0076] Based on the established unified spatiotemporal reference, to eliminate the interference of Earth's curvature and sensor noise on the true undulations of waves, the data acquisition and preprocessing module 10 executes step S4: coordinate transformation and detrending. This is specifically broken down into the following sub-steps:

[0077] Sub-step S401: Local tangent plane coordinate system transformation. Based on the principle that the Earth is approximately an ellipsoid and wave motion manifests as local tangent plane physical undulations, the geodetic coordinate system parameters are retrieved from the synchronous state data matrix. For the transformation from the geodetic coordinate system to the local tangent plane coordinate system, those skilled in the art can use conventional reference ellipsoid projection calculation methods; the transformation process is well-known in the field and will not be elaborated here. After the transformation, the northward displacement sequence and the eastward displacement sequence are separated, and the water surface undulation data with actual physical significance are extracted using the vertical elevation detrending formula. The system incorporates anomaly detection logic when calculating the mean value for a given period to ensure that the total number of valid sampling points participating in the summation is always greater than zero. The vertical elevation detrending formula is:

[0078] ;

[0079] in, Indicates time The original vertical displacement components; Indicates time The height of the ellipsoid; This represents the arithmetic mean of the ellipsoidal height within the observation time window.

[0080] Sub-step S402: High-pass detrending filtering for long-period drift. To eliminate extremely low-frequency water level changes caused by ocean tides and baseline drift caused by receiver thermal noise, a fourth-order Butterworth high-pass filter with a cutoff frequency set to 0.03Hz is simultaneously applied to the converted displacement sequences in the three directions. To avoid edge step oscillations introduced during filter initialization, mirror extension padding is performed at both ends of the input sequence before filtering. After high-pass filtering, a three-dimensional detrended displacement sequence is separated. This three-dimensional detrended displacement sequence includes a detrended northward displacement sequence, a detrended eastward displacement sequence, and a detrended vertical displacement sequence. The detrended vertical displacement sequence is directly input into the subsequent anomaly detection module, while the detrended northward and detrended eastward displacement sequences are cached in storage space for verification of the physical laws of circular motion of water particles in step S13.

[0081] In this embodiment, after removing baseline drift and establishing a unified spatiotemporal reference, the system needs to investigate potential hardware-level failure risks from the communication physical layer. To this end, the anomaly detection and feature extraction module 20, connected to the data acquisition and preprocessing module 10, executes step S5: signal quality anomaly pre-labeling. This involves extracting the state parameters at the timestamps corresponding to the three-dimensional detrended displacement sequence, comparing them with the state boundary values ​​in the decision threshold parameter set, and constructing a signal quality marker sequence. This step is specifically detailed into the following sub-steps:

[0082] Sub-step S501: Extraction of underlying state parameters. The spatial mapping node in the data caching system is retrieved, and for each discrete timestamp, the corresponding number of visible satellites, the average signal-to-noise ratio of each satellite carrier, and the position accuracy factor are extracted. As a preferred approach, the system synchronously starts an independent counter to count the accumulated duration of continuous signal interruption at the current moment. The physical reason for choosing these specific parameters is that they can cross-validate the transient integrity of the link from three independent dimensions: obstruction conditions, electromagnetic interference, and constellation geometry.

[0083] Sub-step S502: Joint decision-making and flag generation of state limits. The extracted state parameters are compared point-by-point with pre-loaded thresholds. To avoid relying on a single dimension for biased judgment, the system uses multi-dimensional logical OR operations for comprehensive evaluation. The signal quality flag sequence calculation formula is used to perform binary quantization of the communication health at each time step, generating a signal quality flag sequence strictly aligned on the time axis. The signal quality flag sequence calculation formula is:

[0084] ;

[0085] in, Indicates time The signal quality flag value has a value of 1, which indicates a risk of hardware lockout or location degradation, and a value of 0, which indicates good communication health. Indicates time The number of visible satellites; This represents the minimum number of visible satellites required to determine a valid location. Indicates time Position accuracy factor; This indicates the upper limit threshold of the position precision factor; Indicates time The duration of continuous no-signal interruption; Indicates the maximum allowed continuous interruption duration of the signal.

[0086] After completing the physical-level link health check, the system needs to discretize the continuous dynamic morphology of sea surface undulations into independently analyzable periodic wave surfaces. Based on this technical objective, the anomaly detection and feature extraction module 20 executes step S6: wave unit segmentation. This involves obtaining the vertical component in the three-dimensional detrended displacement sequence, using the zero-overlap method to search for the intersection of the wave surface curve with the zero mean line for single-wave segmentation, and combining this with the corresponding interval range of the signal quality indicator sequence to obtain a set of independent wave units. This step is further detailed into the following sub-steps:

[0087] Sub-step S601 involves wavefront intersection positioning based on the zero-crossing point method. It searches for adjacent data point pairs in the vertical displacement sequence that change from negative to positive values. To overcome the time resolution bottleneck of discrete sampling, a precise zero-crossing time interpolation formula is used to solve for the theoretical intersection point. Before performing this division operation, the system has built-in still water surface shielding logic. When the absolute value of consecutive sampling points is less than a set minimum floating-point limit (e.g., 10⁻⁶ meters), the current interpolation process is forcibly skipped, eliminating the risk of indeterminate form crashes caused by the denominator approaching zero at the underlying code level. The precise zero-crossing time interpolation formula is:

[0088] ;

[0089] in, Indicates the first A precise zero-crossing moment; Indicates the discrete sampling time before crossing zero; It represents the absolute value of the vertical displacement at the discrete sampling moment before crossing zero; Indicates the discrete sampling time after crossing zero; Represents the vertical displacement at the discrete sampling time after crossing zero; This represents the discrete sampling time interval of the system.

[0090] Sub-step S602: Single wave interception and period calculation. Using two adjacent precise zero-crossing moments as time boundaries, the continuous water level sequence is divided into independent wave period segments. The specific duration of each wave is calculated using the wave unit period calculation formula, which is:

[0091] ;

[0092] in, Indicates the first The period span of a wave unit; Indicates the first A precise zero-time point. After completing the segmentation operation, the system projects the time period corresponding to each independent cycle onto the signal quality indicator sequence, thereby constructing a set of independent wave units containing corresponding state labels.

[0093] Since real ocean background waves are constantly evolving, using a static absolute height threshold can easily lead to misjudgments under extreme sea conditions. To construct a local energy benchmark with environmental adaptability, the anomaly detection and feature extraction module 20 executes step S7: dynamic estimation of salient wave height. Within a sliding time window, the zero-order spectral moment of the set of independent wave units is calculated to infer the salient background wave height. By comparing the individual wave heights with the salient background wave heights, a set of suspected anomalous candidate waves is extracted. This step is specifically broken down into the following sub-steps:

[0094] Sub-step S701: Local sea state benchmark assessment within the sliding window. A statistical time window of 20 to 30 minutes is constructed, extending forward and backward in both directions from the currently analyzed wave as the center. This duration is determined to ensure that the window contains at least 100 complete wave groups to meet the statistical large sample requirement of the Rayleigh distribution, while avoiding the inclusion of non-stationary trends at the tidal scale due to an excessively large window. The variance of wavefront displacement within this window is calculated using the zero-order spectral moment estimation formula:

[0095] ;

[0096] in, Represents the zeroth-order spectral moment within the sliding window; This represents the total number of discrete sampling points within the sliding window; This indicates a summation operation applied to all sampled points within the window. Indicates the first [number]th ... The vertical displacement at each sampling time point. Subsequently, the macroscopic background benchmark is derived using the significant wave height estimation formula, which is:

[0097] ;

[0098] in, This indicates a significant background wave height. Sub-step S702 involves range traversal and initial screening of suspected anomalous waves. For each wave band included in the above set, its extreme span is determined using the single-wave height calculation formula, which is:

[0099] ;

[0100] in, Indicates the first The wave height of a single wave unit; Indicates the first The time window of a wave unit; This indicates an operation to retrieve the maximum value within the corresponding time window; This indicates that the minimum value will be retrieved within the corresponding time window. Indicates time The specific values ​​in the detrended vertical displacement sequence are then used to construct the total horizontal displacement power spectrum. This leads to the construction of a logic for comparing the single-wave height of each wave unit with the significant background wave height using a set multiple. Specifically, when the single-wave height is... Furthermore, a total horizontal displacement power spectrum is constructed that is greater than a significant wave height multiple threshold. Furthermore, the total horizontal displacement power spectrum and the significant background wave height were constructed. When constructing the product of the total horizontal displacement power spectrum, bands exceeding the upper limit of conventional energy statistics are extracted to form a set of suspected anomalous candidate waves.

[0101] After screening candidate wavebands with anomalous height jumps, the positioning artifacts caused by satellite signal jumps are often accompanied by non-physical high-frequency oscillations, exceeding the kinematic limits of water particles constrained by gravity. To provide a basis for judgment from the underlying mechanical mechanism, the anomaly detection and feature extraction module 20 executes step S8: kinematic feature calculation. It uses a central difference scheme to numerically differentiate the displacement sampling points within the suspected anomalous candidate wave set, and calculates kinematic feature parameters including vertical acceleration and wave steepness based on the deep-water dispersion relation. This step is specifically detailed into the following sub-steps:

[0102] Sub-step S801: Higher-order numerical differentiation and extremum search. For each discrete sampling point within the time window, a higher-order central difference algorithm is used to obtain the transient derivative. Compared to conventional unidirectional difference, this method significantly suppresses the amplification effect of the difference process on high-frequency thermal noise through symmetric mask convolution. The solutions are obtained using the five-point central difference formulas for vertical velocity and vertical acceleration, respectively.

[0103] ;

[0104] ;

[0105] in, Indicates time The vertical velocity; Indicates time The vertical acceleration; This represents the discrete sampling time interval. Sub-step S802: Wave steepness dimension feature derivation. Wave steepness is a key dimensionless criterion for assessing whether the wave crest geometry approaches the hydrodynamic breakage critical point. Based on ocean dynamics principles, the wavelength spatial scale is obtained using the deep-water dispersion relation wavelength calculation formula. The deep-water dispersion relation wavelength calculation formula is as follows:

[0106] ;

[0107] in, Indicates the first The deep-water wavelength of each wave unit; Represents the gravitational acceleration constant; This represents the constant pi. Subsequently, the wave steepness calculation formula, combined with single-wave height, is used for calculation. To prevent high-frequency glitches caused by positioning noise from affecting the denominator period... The extremely small value leads to wave steepness overflow. The system adds a physical boundary interception before this value is applied, forcibly classifying wavebands with periods less than 1 second as artifact jumps. The wave steepness calculation formula is:

[0108] ;

[0109] in, Indicates the first The wave steepness of each wave unit. The extreme state variables obtained from the solution will be used as kinematic characteristic parameters and output to downstream modules to support the final joint decision.

[0110] In this embodiment, to accurately isolate non-physical numerical jumps caused by hardware transient loss of lock, the system enters a multi-dimensional in-depth verification and wavefront timing repair stage. As a preferred approach, the decision and reconstruction module 30, connected to the anomaly detection and feature extraction module 20, executes step S9: physical limit violation decision. This step aims to compare the kinematic feature parameters with the dynamic limit limits in the decision threshold parameter set, and generate a physical anomaly trigger state when the dynamic limit limits are exceeded. Since the inherent error jump of the satellite positioning system manifests as a spatial jump of a massless geometric point, its derived transient higher-order derivatives inevitably ignore the mechanical laws of gravity constraining actual water bodies. Based on the above physical mechanism, this step is specifically refined into the following sub-steps:

[0111] Sub-step S901: Acceleration-Gravity Limit Comparison. Retrieve the acceleration extreme value sequence from the previously input kinematic characteristic parameters. Based on actual hydrodynamic constraints, the downward acceleration of free-floating wave particles theoretically cannot exceed Earth's gravity; any data violating this law can be directly identified as a system artifact. The acceleration limit violation judgment formula is used to evaluate the transient acceleration peak value within this band. The acceleration limit violation judgment formula is as follows:

[0112] ;

[0113] in, This indicates that the maximum acceleration value is taken within the corresponding decision time window; This represents the decision time window for the i-th wave unit; Indicates time The absolute value of the vertical acceleration; Indicates time The vertical acceleration; This represents the upper limit threshold of the ratio of acceleration to gravitational acceleration; This represents the gravitational acceleration constant. In this embodiment, The value range is set between 1.15 and 1.25. This range is determined because it strictly tightens the theoretical mechanical upper limit of natural gravity waves while reserving the necessary tolerance for the reasonable basis noise of the sensor hardware itself.

[0114] Sub-step S902, wave steepness breaking limit check. Wave steepness, as a dimensionless quantity reflecting the geometric sharpness of wave crests, also has a theoretical threshold. According to the nonlinear Stokes wave theory, the limiting wave steepness of deep-water gravity waves, i.e., the Michell limit, is approximately 1 / 7 (i.e., 0.142). Wave patterns exceeding this limit will undergo wave breaking, thus failing to maintain continuous wave surface evolution. The degree of morphological distortion is checked using the wave steepness limit exceeding judgment formula. The wave steepness limit exceeding judgment formula is as follows:

[0115] ;

[0116] in, Indicates the first The wave steepness of each wave unit.

[0117] As a preferred approach, when either the aforementioned acceleration limit condition or wave steepness distortion condition is triggered, the decision and reconstruction module 30 determines that the current waveband has exceeded the dynamic limit, and then generates a physical anomaly triggering state for that waveband in the memory bus, denoted as... If none of these are triggered, the physical anomaly state is recorded as follows: .

[0118] Relying solely on physical characteristics can easily lead to the accidental deletion of rare and abnormal waves that actually exist under extreme sea conditions. Cross-domain logical fusion can effectively compensate for this blind spot. Therefore, the decision and reconstruction module 30 executes step S10: multi-dimensional feature fusion decision. This step combines the physical anomaly trigger state with the state identifiers corresponding to the suspected anomaly candidate wave set to perform logical judgment and extract the positioning artifact wavebands. This process is specifically detailed into the following sub-steps:

[0119] Sub-step S1001: Multi-source data time alignment. The previously generated signal quality flag sequence is projected into the decision time window of the corresponding i-th wave unit to ensure strict synchronization between the underlying operating parameters and the surface mechanical characteristics on the absolute time axis. For the specific underlying code implementation of multi-source heterogeneous data underlying clock synchronization and ring buffer alignment, those skilled in the art can use a standard hardware second pulse matching algorithm combined with an interrupt service routine. Its data bus scheduling is a well-known technology in this field and will not be elaborated upon here.

[0120] Sub-step S1002: Multi-dimensional logic synthesis and decision. After obtaining the synchronized aligned feature vectors, a joint conditional judgment is performed using the multi-dimensional feature fusion decision logic formula. The multi-dimensional feature fusion decision logic formula is as follows:

[0121] ;

[0122] in, Indicates the first Multidimensional feature fusion judgment results of wave units; Indicates the first The cumulative value of the signal quality indicator sequence within the decision time window of each wave unit; Represents the logical OR operator; Indicates a physical anomaly trigger state; This represents the logical AND operator. Through this cross-logic, the system accurately locates and extracts the positioning artifact band, providing precise target coordinates for subsequent signal repair.

[0123] The removal of localization artifacts inevitably leaves data gaps in the original sequence, thus disrupting the continuity prerequisite required for frequency domain transformation. Based on this, the decision and reconstruction module 30 executes step S11: time-series reconstruction based on an autoregressive model. The system measures the signal gap span of the localization artifact band, extracts the outer normal observation sampling points, calculates the autoregressive model coefficients, and performs bidirectional prediction weighted fusion to replace the original abnormal data, generating a time-domain reconstructed vertical sequence. This process is specifically detailed into the following sub-steps:

[0124] Sub-step S1101: Model well-being test and sample extraction. The signal gap span of the localization artifact band is precisely measured, and equal time intervals are extended to both sides of this span. Normal observation sampling points on the outer side are extracted as historical regression samples for constructing the prediction equation. To prevent the main program from crashing due to singularities generated by covariance matrix inversion during the underlying solution, the system forcibly pre-calculates the variance of the extracted sample sequence before any matrix division or inversion instructions. When this variance falls below a set machine floating-point minimum limit, the covariance matrix is ​​determined to be approaching degradation. The system immediately bypasses the autoregressive solution process and downgrades to low-order polynomial smooth interpolation to fill gaps, ensuring the algorithm has complete integrity under any still-water conditions.

[0125] Sub-step S1102: Model coefficient calculation and bidirectional weighted fusion. For sample data that passes the well-posedness test, the Yule-Walker equation is called to calculate the prediction parameters, and the autoregressive model time series reconstruction formula is used to extrapolate the data in the missing interval. The autoregressive model time series reconstruction formula is:

[0126] ;

[0127] in, Indicates time The time-domain reconstruction of vertical displacement; This indicates a cumulative summation operation based on the number of backtracking steps from historical sampling points; Indicates the order of the autoregressive model; Indicates the first The autoregressive coefficients corresponding to each discrete sampling point; Indicates time The previous The detrended vertical displacement at each discrete sampling point at the corresponding time. Indicates the discrete sampling time interval; Indicates time The white noise error term.

[0128] As a preferred approach, to eliminate the artificial high-frequency steps caused by forced splicing at the end of the sequence in unidirectional autoregressive predictions, the system strictly implements bidirectional prediction weighted fusion. It not only extrapolates to the right based on the normal data at the outer left end, but also similarly extrapolates to the left based on the normal data at the outer right end. Subsequently, based on the time interval between each estimated target point within the missing span and the boundaries of the left and right true data, inverse weights are assigned, and the predicted values ​​from both directions are weighted and fused. This seamless interpolation replaces the original positioning artifact bands, ultimately outputting a time-domain reconstructed vertical sequence for final frequency domain verification.

[0129] In this embodiment, after wavefront continuity restoration at the preceding time-domain level, simple mathematical interpolation or prediction models inevitably introduce unnatural high-frequency artifacts implicitly. To verify the physical authenticity of the time-series reconstruction results from a fluid dynamics perspective, the verification and output module 40, connected to the decision and reconstruction module 30, executes step S12: spectral morphology consistency check. This step performs a fast Fourier transform on the reconstructed vertical sequence in the time domain to obtain the power spectral density distribution, and performs univariate linear regression fitting on the high-frequency equilibrium band in a double logarithmic coordinate system to obtain the energy attenuation slope. Specifically, this is broken down into the following sub-steps:

[0130] Sub-step S1201: Time-frequency domain transformation and high-frequency band definition. The system receives the input time-domain reconstructed vertical sequence and applies a discrete frequency domain transformation algorithm to perform time-frequency domain mapping. For the specific process of obtaining the Fast Fourier Transform and power spectral density distribution of this sequence, those skilled in the art can use the Welch spectrum estimation algorithm with Hanning window smoothing. Its underlying periodogram averaging processing is a well-known technique in the field and will not be elaborated here. After obtaining the full-band spectrum, accurately stripping the main wave energy is crucial for high-frequency verification. The system optimizes and locates the spectral peak frequency corresponding to the main energy peak, and extends and truncates towards the high-frequency direction starting from a set multiple (e.g., 1.5 to 2.5 times) of this spectral peak frequency, extracting high-frequency balanced band data that has deviated from the main wave energy region for subsequent analysis.

[0131] Sub-step S1202: Fitting the high-frequency attenuation law. Based on Phillips wave spectrum theory, in deep-water gravity wave areas unaffected by shallow-water topographic friction, the high-frequency tail energy of wind waves in an energy equilibrium state exhibits a fixed power-law attenuation law with increasing frequency. If the repaired wave surface data still retains non-smooth microsteps caused by lock-up artifacts, its high-frequency spectrum will exhibit a flat distribution similar to white noise, thus violating the aforementioned established physical law. The frequency domain morphology is verified using the high-frequency equilibrium region spectral density attenuation law formula, which is:

[0132] ;

[0133] in, Represents frequency The power spectral density; Indicates frequency; The sign represents a direct proportion; This represents the frequency raised to the power of negative five. The verification and output module 40 projects the extracted high-frequency equilibrium segment data onto a double logarithmic coordinate system to perform univariate linear regression. It then calculates the spectral slope deviation between the fitted energy decay rate and the theoretical value -5. When the spectral slope deviates Less than the set spectral slope deviation threshold At that time, it was determined that the reconstructed wavefront did not cause high-frequency distortion. As a preferred method, this spectral slope deviation threshold... The value is set to 0.5 to 1.0 (corresponding to an energy decay rate within the allowable tolerance range of -4.5 to -5.5 or relaxed to -4.0 to -6.0), and its definition is based on a comprehensive consideration of the instability of wind energy injection in actual dynamic sea conditions and the inherent high-frequency noise of the buoy sensor.

[0134] Is a one-dimensional vertical spectral morphology test sufficient to fully verify the rationality of wave motion? The answer is no. Based on the mechanical constraint mechanism of the three-dimensional spatial orbit coupling of gravity waves and water particles, the verification and output module 40 further executes step S13: closed-loop verification of spectral domain physical consistency. This step verifies the theoretical deviation of the energy attenuation slope, extracts the horizontal component from the three-dimensional detrended displacement sequence, and calculates the mean square error of the spectral verification factor with the time-domain reconstructed vertical sequence. Specifically, it is broken down into the following sub-steps:

[0135] Sub-step S1301: 3D motion spectrum alignment and verification factor calculation. Based on the fundamental characteristics of the linear theory of micro-amplitude waves, under deep-water conditions, the sum of the variances of the displacement energy in the two orthogonal directions within the horizontal plane should theoretically be strictly equal to the variance of the vertical displacement energy. To verify this characteristic, the system simultaneously extracts the detrended eastward and detrended northward displacement sequences from the 3D detrended displacement sequence, calculates their respective power spectral densities, and performs linear equal-weight superposition to construct the total horizontal displacement power spectrum. Subsequently, combined with the vertical spectrum data obtained in the preceding steps, the multidimensional coupling ratio is calculated using the spectral domain physical consistency verification factor formula. The spectral domain physical consistency verification factor formula is:

[0136] ;

[0137] in, Represents frequency Spectral domain physical consistency check factor; Represents the square root operator; Represents frequency The vertical hydrodynamic power spectral density; Represents frequency The total horizontal displacement power spectral density. Under ideal conditions with no positioning artifacts and the buoy in a deep-water hydrodynamic environment, the theoretical baseline value of this factor constantly approaches 1. To prevent system crashes caused by the denominator approaching zero under extreme still water conditions, the system forcibly checks the frequency before performing division. The system checks whether the total horizontal displacement power spectral density is lower than the preset minimum value of the machine floating point. If it is lower than this minimum value, it is determined that the current ocean lacks real physical wave excitation. The system immediately bypasses the calculation operation at this frequency and forcibly assigns it the theoretical reference value of 1, thereby strictly ensuring the logical integrity of the algorithm instructions under all sea conditions.

[0138] Sub-step S1302: Full-band error assessment and closed-loop decision. To quantify the global physical deviation across the entire major wave frequency band, the system defines the core frequency band containing the main wave energy. The global error is assessed using the mean square error of the spectral check factor, which is calculated as follows:

[0139] ;

[0140] in, This represents the mean square error of the spectral check factor; This indicates the total number of frequency points within the valid test band; This indicates a cumulative summation operation for all frequency points within the valid test band; This indicates the lower limit frequency of the valid test band; Indicates the upper limit frequency of the valid test band; Represents frequency The spectral domain physical consistency verification factor. In this embodiment, the effective verification frequency band is not a fixed constant, but is dynamically calculated based on the peak frequency. The calculation is based on selecting the core frequency band with an energy accumulation ratio of 90% as the judgment interval. This strategy can effectively avoid the dual nonlinear interference of low-frequency instrument drift and high-frequency clutter from a physical perspective.

[0141] Based on the global error results obtained from the multi-dimensional coupling calculations, the system compares the mean square error of the spectral verification factor with a preset upper limit. As a preferred approach, this upper limit is set to 0.05. When the calculated global error is significantly lower than the upper limit, it strongly demonstrates that the three-dimensional trajectory repaired by time-domain interpolation fully conforms to the fundamental laws of ocean hydrodynamics in both the time and frequency domains. Based on this determination, the verification and output module 40 then officially outputs high-quality ocean wave parameters, thus completing a closed-loop data purification and reconstruction verification cycle with rigorous physical self-consistency.

[0142] In this embodiment, based on the hardware platform architecture including a BeiDou satellite navigation receiver, satellite antenna, and data processing unit, system-level testing and verification were conducted using a specific marine surveying engineering task. The observation site was selected in the nearshore waters of Qingdao (coordinates approximately 36.045°N, 120.432°E), with a water depth of approximately 20m. The observation period was from July 13th to July 24th, 2024, generating a total of 576 sets of standard half-hour data. The observation platform used a self-developed BeiDou wave buoy, and to provide an independent hydrodynamic benchmark, a Datawell Waverider accelerometer reference buoy was deployed in the same waters. To ensure strict spatial and temporal alignment of the dual-source heterogeneous data, both systems underwent high-precision time synchronization based on GPS time before entering the water, and the mooring point spacing was controlled within 50 meters, thus eliminating system errors caused by the spatial inhomogeneity of the local wind and wave field. During the test, the buoy experienced various sea states from small to medium waves, with significant background wave height. The distribution range is between 0.4 and 1.4 m.

[0143] With this hardware deployment, the data processing unit executes step S1 (system initialization and parameter configuration) and step S2 (multi-source raw data stream acquisition) according to the configuration instructions in the decision threshold parameter set. Subsequently, the system executes step S3 (time base alignment and synchronization) to generate a synchronization status data matrix. Data from 06:00 to 06:30 on July 22, 2024, is extracted as a typical test example. During this period, the background wave height is significantly high. The amplitude of the background wavefront displacement is approximately 0.78m. The buoy was floating, which is typical of moderate wave conditions. However, a short-term loss of lock-on event lasting approximately 60 seconds was recorded at the bottom of the buoy. During this period, the data acquisition and preprocessing module 10 performed step S4: coordinate transformation and detrending, separating the three-dimensional detrended displacement sequence. Simultaneously, the anomaly detection and feature extraction module 20 performed step S5: signal quality anomaly pre-labeling. When the number of visible satellites at the bottom layer was detected... Less than 5 and position precision factor threshold When the value exceeds 5.0, the signal quality flag sequence at the current moment will be updated. Set to 1.

[0144] See attached document Figure 3 , Figure 3 The layout is arranged vertically, (a) is the wavefront displacement timing sequence, and (b) is the signal quality indicator sequence. The temporal changes. Figure 3 In subgraph (a), at time Near point s, the raw observation data exhibited a sharp-then-gradient displacement jump, with an initial pulse peak amplitude of approximately 5m, followed by a gradually decaying oscillating wake, lasting approximately 60 seconds. This waveform evolution pattern closely matches the typical artifact pattern generated by the internal high-pass filtering mechanism after the BeiDou signal is lost. Considering the abnormal signal quality, the system executes step S6: wave unit segmentation, obtaining an independent set of wave units. Immediately following, within the sliding time window, step S7: dynamic estimation of significant wave height to calculate the background significant wave height is executed. The system extracts a set of candidate waves suspected of being anomalies. Based on this, the system calculates kinematic feature parameters, including vertical acceleration and wave steepness, through step S8: kinematic feature calculation. The parameters of the four suspected anomaly units extracted at this time are shown in Table 1.

[0145] Table 1. List of Detailed Indicators for Abnormal Wave Unit Detection

[0146]

[0147] Since single-dimensional extreme value judgment is prone to failure under complex sea conditions, the decision and reconstruction module 30 executes step S9: physical limit violation judgment and step S10: multi-dimensional feature fusion judgment. By combining the physical anomaly trigger state and the underlying state identifier, the positioning artifact band is accurately extracted. For the removed distorted data, the system executes step S11: time-series reconstruction based on the autoregressive model. As a preferred method, this reconstruction process extracts the coefficients of the autoregressive model from the outer normal observation sampling points and assigns the order of the autoregressive model... The value is set to 8. This value is determined based on the Akaike information content criterion assessment of the historical swell spectrum of the test area, ensuring coverage of the energy of a single main cycle while avoiding the introduction of high-frequency overfitting oscillations. After completing the bidirectional prediction weighted fusion, the system finally generates a time-domain reconstructed vertical sequence, as shown in... Figure 3 As shown by the solid green line in subplot (a), the repaired waveform achieves a smooth and natural transition with the normal sequence before and after the abnormal segment in both amplitude and phase.

[0148] See attached document Figure 4 The smoothing of the time-domain waveform must rely on the physical laws of the frequency domain for closed-loop verification. Therefore, the verification and output module 40 executes step S12: spectral morphology consistency check. The system performs univariate linear regression fitting on the high-frequency balanced band, where the specific range of this high-frequency balanced band is rigorously defined as follows: to Spectral peak frequency This represents the frequency nodes where the main energy of ocean waves is concentrated. Deep-water waves within this range should strictly adhere to the Phillips equilibrium spectrum theory. Energy decay law. Figure 4The blue shaded area represents the high-frequency equilibrium region. The graph visually compares the energy attenuation evolution of the power spectral density in this frequency band before and after the restoration.

[0149] Long-term operational performance directly reflects the system's business availability. (See attached document.) Figure 5 ,exist Figure 5 Subplot (a) shows a panoramic view of 12 days of continuous observation, with fixed system biases between the two independent instruments pre-emptively eliminated; Figure 5 Subplot (b) shows a magnified view of the location artifact on the same day. Verification and output module 40 executes step S13: spectral domain physical consistency closed-loop verification. When the confirmation error is below the specified lower limit, high-quality ocean wave parameters are output. The statistical error matrix of the long-term comparison is shown in Table 2.

[0150] Table 2. Summary of errors in the comparison of wave statistical parameters from long-term continuous observation.

[0151]

[0152] Note: "-" indicates that no specific statistical and comparative calculation of the correlation coefficient of this parameter was performed under this experimental verification dimension. It should be added that, in extreme still water conditions, to prevent the denominator in the relative deviation calculation from approaching zero, the system has a built-in minimum value bypass logic. If the reference wave height is below 0.05m, the error ratio calculation is paused, thus ensuring the crash-resistant integrity of the computational kernel.

[0153] Extreme errors are relatively easy to filter out, but the ability to resist deception in the face of low-amplitude, hidden anomalies is the core criterion for testing a system's robustness. (See attached document) Figure 6 ,exist Figure 6 Subfigure (a) illustrates the missed detection scenarios of wavefront displacement time series and statistical threshold method. Figure 6 Subgraph (b) shows the signal quality indicator sequence corresponding to short-term signal loss. To intercept such anomalies, the system incorporates a detection and early warning threshold coefficient in the dynamic estimation phase. The value was lowered to 1.8. This setting is based on the cutoff property of the Rayleigh distribution transcendental probability function, effectively broadening the initial screening range.

[0154] Based on the above implementation steps and multi-dimensional test data, the following in-depth conclusions and analyses are given regarding the operating mechanism and final effectiveness of the Beidou wave buoy signal lock-up artifact detection and spectral consistency repair system in complex marine environments:

[0155] First, examining the physical boundary characteristics revealed in Table 1, in the typical loss-of-lock event on July 22, the system effectively isolated four anomalous units. Among them, unit 2, as the main anomaly, had a measured single-wave height. Up to 7.28m, equivalent to a significant background wave height It is 9.4 times that of the previous wave. Its peak single-wave acceleration... Soaring to 12.6 m / s 2 It exceeded the set 1.2 times the gravitational acceleration (approximately 11.8 m / s²). 2 Physical threshold; and its single-wave steepness It reaches 0.223, far exceeding Michell's theoretical limit of 0.142 for deep-water wave collapse.

[0156] Based on the aforementioned kinematic hard indicators, Units 1 and 2 were accurately identified as positioning artifacts exhibiting both signal quality degradation and physical limit violations. However, after their corresponding unlock markers returned to normal, Units 3 and 4 did not trigger any over-limit features in terms of wave steepness and acceleration, indicating that the residual effects of the hardware filtering link had dissipated and the wave surface morphology had returned to normal. Since the single wave heights of Units 3 and 4 were still greater than the background baseline and thus included in the set of suspected anomalous candidate waves, this solution accurately identified them as real waves and retained them through the joint logic of the anomaly detection and feature extraction module 20 and the decision and reconstruction module 30, thus avoiding the erroneous deletion phenomenon that is easily caused by the single threshold method that relies solely on wave height statistics.

[0157] The closed-loop verification in the frequency domain further solidifies the reliability of the time-domain restoration. Based on... Figure 4 Quantization calculations show that the unprocessed red anomalous sequence has an attenuation slope of only -2.63 in the high-frequency range, which deviates from the theoretical value of -5. The value reached 2.37, which significantly exceeded the system's set threshold for spectral slope deviation. This non-physical broadband energy expansion directly reflects the destructive impact of the approximately 60-second BeiDou lock-off event on wave morphology. After system reconstruction, the attenuation slope of the green repair curve in this frequency band significantly decreased to -4.04, corresponding to a spectral slope deviation. Reduced to 0.96 (satisfying less than) (Constraints). The high-frequency slope was improved by 1.41 before and after the repair, which confirmed from the perspective of fluid dynamics mechanism that the AR model effectively suppressed pseudo clutter and did not introduce redundant high-frequency artifacts into the original wind and wave spectrum.

[0158] The effectiveness of adjustments at the macro-statistical level determines the final usability of the data product. (Referring to Table 2 and...) Figure 5 In long-term comparative testing, during a continuous 12-day all-weather comparison, the significant wave height produced by this invention was demonstrated. The average deviation from the co-located accelerometer reference was reduced to 0.033m, with a correlation coefficient of 0.99, fully demonstrating the system-level accuracy. (Analysis in...) Figure 5As shown in the anomalous slice at 06:00 in subplot (b), the uninterrupted original distortion caused a significant false surge in wave height from 0.78m to 1.06m, with a false rise rate as high as 47% (as shown by the red trajectory in the figure). However, with the intervention of the system verification and output module 40, the green dashed line representing the reconstruction result smoothly fell back to the background baseline level of 0.75m, highly coinciding with the measurement trajectory of the Waverider reference buoy.

[0159] The deeper methodological superiority is reflected in Figure 6 This is an example of countering small-amplitude concealment artifacts. In this test case, the buoy is... The nearby area experienced a transient loss of lock for only about 3.5 seconds, resulting in a displacement artifact peak of only 1.5m (approximately...). The subsequent decaying oscillation duration is only about 15 seconds. Under such typical transient interference, the traditional statistical threshold method based on 4 times the significant wave height (the interception threshold requires...) The artifact completely failed due to excessively low interference amplitude. Simultaneously, the single-wave steepness of this artifact band... The single-wave acceleration peak value did not reach the 0.142 limit. It also did not exceed the 1.2 times gravitational acceleration constraint, rendering purely physical limitation tests ineffective in this scenario. At this point, the system made a cross-domain call to... Figure 6 The signal quality indicator sequence at the bottom layer of the subgraph (b) The system was marked and activated in conjunction with a low-sensitivity detection warning threshold coefficient. The primary screening network effectively elevates this hidden clutter to the verification pool and removes it.

[0160] These two sets of experimental cases construct a robust complementary closed loop in terms of verification logic: Case 1 (pulse peak value approximately 5m, equivalent to...) This verifies the system's ability to perform end-to-end detection, time-domain restoration, and frequency-domain physical verification of typical large-scale destructive artifacts; while Case 2 (pulse peak approximately 1.5m, equivalent to...) This strongly demonstrates that the cross-domain fusion framework of this scheme possesses extreme sensitivity and detection lower limit that traditional mathematical statistics methods do not have when stripping away hidden micro-artifacts.

Claims

1. A method for detecting artifacts in BeiDou wave buoy signal lock-up and repairing spectral consistency, characterized in that, Includes the following steps: Read the observation messages and extract the raw measurement data sequence, align the timestamps to generate a synchronization status data matrix; Based on the aforementioned synchronization state data matrix, local tangent plane transformation and high-pass filtering are performed to separate the three-dimensional detrended displacement sequence; Based on the three-dimensional detrended displacement sequence, a signal quality indicator sequence is constructed and single-wave segmentation is performed. Then, the single-wave height is compared with the background significant wave height, and the kinematic characteristic parameters are calculated based on the comparison results. Based on the kinematic characteristic parameters and the corresponding state identifier, the dynamic limit is compared and logical judgment is performed. The positioning artifact band is extracted, the autoregressive model coefficients are calculated, and bidirectional prediction weighted fusion is implemented to replace the original abnormal data to achieve the repair of signal lock artifacts and generate a time-domain reconstructed vertical sequence. The energy decay slope is obtained by performing a fast Fourier transform on the reconstructed vertical sequence in the time domain to verify the theoretical deviation. The mean square error of the spectral verification factor is calculated by extracting the horizontal component in the three-dimensional detrended displacement sequence and the reconstructed vertical sequence in the time domain. When the error is lower than the specified upper limit, high-quality ocean wave parameters are output.

2. The method for detecting artifacts and restoring spectral consistency of Beidou wave buoy signal lock-up according to claim 1, characterized in that, In the steps of reading observation reports and extracting raw measurement data sequences, aligning timestamps to generate the synchronization state data matrix, and performing local tangent plane transformation and high-pass filtering based on the synchronization state data matrix to separate the three-dimensional detrended displacement sequence: If and only if the difference between adjacent timestamps is strictly greater than zero, a linear interpolation algorithm is used to map and align the extracted original measurement data sequence to the corresponding sampling time node, and the result is merged to construct the synchronization state data matrix. The geodetic coordinate system parameters in the synchronous state data matrix are retrieved for local tangent plane transformation to separate the northward displacement sequence and the eastward displacement sequence. The original vertical displacement component is extracted by subtracting the arithmetic mean of the ellipsoidal height within the observation time window from the current ellipsoidal height using the vertical elevation detrending formula. High-pass filtering is applied to the northward displacement sequence, the eastward displacement sequence, and the original vertical displacement component to eliminate extremely low-frequency drift, thereby separating the three-dimensional detrended displacement sequence, which includes the detrended northward displacement sequence, the detrended eastward displacement sequence, and the detrended vertical displacement sequence.

3. The method for detecting artifacts of signal loss and restoring spectral consistency of Beidou wave buoys according to claim 1, characterized in that, The steps of constructing a signal quality indicator sequence based on the three-dimensional detrended displacement sequence and performing single-wave segmentation specifically include: The state parameters under the corresponding timestamp are extracted and comprehensively evaluated using multi-way logical OR operations to generate a signal quality flag sequence. The vertical component in the three-dimensional detrended displacement sequence is obtained. When the absolute value of the consecutive sampling points is less than the set minimum floating-point number limit, the zero-crossing interpolation is forcibly skipped. The adjacent zero-crossing moments are used as time fences for single-wave segmentation to obtain independent wave period segments. The time periods corresponding to the independent wave cycle segments are projected and combined onto the signal quality indicator sequence to obtain a set of independent wave units.

4. The method for detecting artifacts and restoring spectral consistency of Beidou wave buoy signal lock-up according to claim 3, characterized in that, The step of further comparing the single wave height with the significant background wave height specifically includes: By taking the current analyzed wave in the set of independent wave units as the center and extending it forward and backward in both directions to construct a statistical time window, the variance of the wave surface displacement within the statistical time window is calculated to derive the macroscopic background benchmark and obtain the significant background wave height. The extreme span of the wave band is determined for the wave unit in the set of independent wave units to obtain the single wave height; The single wave height is compared with the background significant wave height. When the single wave height is greater than the product of the significant wave height multiple threshold and the background significant wave height, extraction is performed to extract a set of suspected abnormal candidate waves. Furthermore, the significant wave height multiple threshold is the dimensionless multiplier ratio limit for judging the deviation of transient abrupt wave height from the sea state background energy level.

5. The method for detecting artifacts and restoring spectral consistency of Beidou wave buoy signal lock-up according to claim 4, characterized in that, The step of calculating the kinematic characteristic parameters based on the comparison results specifically includes: Based on the comparison results, for the displacement sampling points within the suspected abnormal candidate wave set, the vertical acceleration is obtained by solving the transient derivative using the five-point center difference formula for vertical acceleration. Based on the principles of ocean dynamics, the wavelength spatial scale is obtained by using the deep-water dispersion relation wavelength calculation formula, and the wave steepness is calculated by combining the single-wave height with the wave steepness calculation formula. The vertical acceleration and the wave steepness are combined to extract the kinematic feature parameters; The kinematic characteristic parameters are checked for limits. When the degree of morphological distortion exceeds the dynamic limit, a physical anomaly trigger state is generated.

6. The method for detecting artifacts of signal loss and restoring spectral consistency of Beidou wave buoys according to claim 5, characterized in that, In the steps of comparing the kinematic characteristic parameters with the corresponding state markers to perform logical judgments on the dynamic limit limits, extracting the positioning artifact bands, and calculating the autoregressive model coefficients: By combining the physical anomaly triggering state with the state identifiers corresponding to the suspected anomaly candidate wave set, a multi-dimensional logical comprehensive decision is made to extract the positioning artifact wave band; The signal loss span of the positioning artifact band is measured, and the outer normal observation sampling points are extracted to form a sample sequence. The variance of the sample sequence is pre-calculated. When the variance is lower than the set minimum limit of the machine floating point, the covariance matrix is ​​determined to be degenerate. The bypass autoregressive solution process is downgraded to use low-order polynomial smooth interpolation to fill the gaps. When the covariance matrix is ​​not degenerate, the autoregressive model coefficients are obtained.

7. The method for detecting artifacts of signal loss and restoring spectral consistency of Beidou wave buoys according to claim 1, characterized in that, In the step of implementing bidirectional predictive weighted fusion to replace the original abnormal data in order to repair signal lock-up artifacts and generate the time-domain reconstructed vertical sequence: By calling the autoregressive model coefficients, the normal data on the left side of the outer side is used as the benchmark to extrapolate to the right, and similarly, the normal data on the right side of the outer side is used as the benchmark to extrapolate to the left, thus obtaining the predicted values ​​in both directions. Extract the missing span in the positioning artifact band, assign inverse weights based on the time interval between each estimated target point within the missing span and the left and right real data boundaries, and weight and fuse the positive and negative predicted values ​​to replace the original abnormal data to generate the temporal reconstruction vertical sequence.

8. The method for detecting artifacts of signal loss and restoring spectral consistency of Beidou wave buoys according to claim 1, characterized in that, In the step of performing a fast Fourier transform on the reconstructed vertical sequence in the time domain to obtain the energy decay slope and verify the theoretical deviation: The time-domain reconstructed vertical sequence is subjected to fast Fourier transform to optimize and locate the spectral peak frequency corresponding to the main energy peak, and the spectral peak frequency is extended and truncated to the high-frequency direction as a set multiple to extract high-frequency balanced frequency band data. The high-frequency balanced band data is projected onto a double logarithmic coordinate system and a univariate linear regression is performed to obtain the energy attenuation slope. The spectral slope deviation between the energy attenuation slope and the theoretical value of the attenuation slope is calculated. When the spectral slope deviation is less than the set spectral slope deviation threshold, the time-domain reconstructed vertical sequence is determined to meet the theoretical deviation test. Furthermore, the theoretical value of the attenuation slope is a pre-set ideal exponential limit characterizing the attenuation law of high-frequency energy of deep-water gravity waves with frequency. The spectral slope deviation threshold is the allowable floating-point error limit used to determine whether the high-frequency energy attenuation characteristics of the reconstructed wavefront deviate from the ideal state.

9. The method for detecting artifacts of signal loss and restoring spectral consistency of Beidou wave buoys according to claim 1, characterized in that, In the step of extracting the horizontal component from the three-dimensional detrended displacement sequence and calculating the mean square error of the spectral verification factor with the time-domain reconstructed vertical sequence, and outputting high-quality ocean wave parameters when the error is below a specified upper limit: The power spectral density of the detrended eastward displacement sequence and the detrended northward displacement sequence in the three-dimensional detrended displacement sequence is calculated and linearly superimposed with equal weights to construct the total horizontal displacement power spectrum. The vertical hydrodynamic power spectral density is derived from the reconstructed vertical sequence in the time domain. The vertical hydrodynamic power spectral density is combined with the total horizontal displacement power spectral density to calculate the spectral domain physical consistency verification factor. The core frequency band containing the main wave energy is defined to assess the global error and obtain the mean square error of the spectral check factor. When the mean square error of the spectral check factor is lower than the specified upper limit, high-quality ocean wave parameters are output. Furthermore, the specified upper limit is the maximum tolerance limit for characterizing the degree of agreement between the two-dimensional spectrum energy of horizontal dynamics and vertical kinematics.

10. A system for detecting artifacts of signal loss and restoring spectral consistency in Beidou wave buoy signals, characterized in that, The method for detecting artifacts and restoring spectral consistency of BeiDou wave buoy signal lock-up as described in any one of claims 1-9 includes: The data acquisition and preprocessing module (10) is used to parse the original observation report and perform coordinate system transformation and detrending filtering to generate a three-dimensional detrending displacement sequence; The anomaly detection and feature extraction module (20) establishes a data connection with the data acquisition and preprocessing module (10) through a data link. The anomaly detection and feature extraction module (20) is used to receive the three-dimensional detrended displacement sequence and segment it into single waves, then compare the single wave height with the background significant wave height, and output a set of suspected anomaly candidate waves and corresponding kinematic feature parameters based on the comparison results. The decision and reconstruction module (30) establishes a data connection with the anomaly detection and feature extraction module (20) through a data link. The decision and reconstruction module (30) is used to perform joint condition judgment and interpolation substitution based on the kinematic feature parameters to generate a temporal reconstruction vertical sequence. The verification and output module (40) establishes a data connection with the decision and reconstruction module (30) through a data link. The verification and output module (40) is used to convert the time-domain reconstructed vertical sequence to the frequency domain to verify the physical laws and output high-quality ocean wave parameters.

Images