A low earth orbit satellite doppler frequency trajectory extraction method and system based on orbit constraint dynamic programming

Abstract

The application discloses a low-orbit satellite Doppler frequency trajectory extraction method and system based on track constraint dynamic programming, which comprises the following steps: firstly, acquiring a time-frequency spectrum of a received signal; then, taking the maximum Doppler acceleration of satellite orbit mechanics as the only scalar prior, determining the maximum allowed jump number of inter-frame frequency index; determining a smoothing penalty coefficient and applying the smoothing penalty coefficient to twice punish the frequency index jump between adjacent time frames; constructing a cost function with a hard constraint of orbit reachable domain on the time-frequency spectrum; finally, under the hard constraint of orbit reachable domain, adopting a Viterbi dynamic programming algorithm to solve a frequency index path which maximizes the cost function, and taking the solved frequency index path as an optimal Doppler frequency trajectory. On the premise of ensuring global optimality, with as little prior information as possible, the method structurally excludes cross-frequency trajectory jump, reduces the calculation complexity, and realizes robust parameter design without experience-based parameter adjustment.

Description

Technical Field

[0001] This invention relates to the field of satellite communication technology, specifically a method and system for extracting the Doppler frequency trajectory of low-Earth orbit satellites based on orbit-constrained dynamic programming. Background Technology

[0002] Low-Earth orbit (LEO) satellite IoT terminals, limited by power and antenna gain, typically operate in environments with extremely low signal-to-noise ratios (SNR), such as -15 dB. To achieve sufficient processing gain, the system generally employs a long-code direct sequence spread spectrum (DSSS) architecture. However, LEO satellites move at high speeds of approximately 7.5 km / s, generating Doppler variation rates on the order of hundreds of Hz / s. Uncompensated high-dynamic Doppler can cause severe energy accumulation losses during the symbol integration period, completely negating the spread spectrum processing gain. Therefore, accurate open-loop Doppler estimation and compensation before despreading is an indispensable and crucial step in achieving reliable communication.

[0003] Existing open-loop Doppler estimation techniques suffer from three main drawbacks. The first is acquisition algorithms based on partially matched filtering or joint time-delay Doppler search. These methods only provide coarse initial frequency values ​​and lack the ability to track continuous frequency change rates. The second is open-loop estimation methods based on FFT spectrum peak finding, including single-segment FFT and piecewise FFT linear fitting. The former treats Doppler as a quasi-static quantity and ignores its rate of change, while the latter is limited by frequency resolution, and independent peak finding in each segment results in large estimation variance at low signal-to-noise ratios. The third is frequency trajectory tracking methods based on time-frequency spectrum (such as Short-Time Fourier Transform (STFT) or Synchronous Compressed Fourier Transform (FSST)). FSST, through instantaneous frequency redistribution, concentrates diffused energy into a high-resolution frequency trajectory, providing a good foundation for subsequent path optimization.

[0004] In the field of time-spectrum frequency trajectory tracking, the Viterbi-type path optimization algorithm based on dynamic programming (DP) is the mainstream framework for this type of method. Kwok and Jones (IEEE TSP, 2000) first applied dynamic programming to frequency trajectory tracking, and many subsequent works have extended it to scenarios such as variable speed machinery fault diagnosis and radar time-frequency interference suppression. However, existing dynamic programming-based frequency trajectory tracking methods all suffer from two common problems that have not been effectively addressed:

[0005] First, there is the reliance on empirically tuned smoothing penalty coefficients. To suppress frequency jumps caused by noise, existing algorithms typically introduce a secondary penalty term for inter-frame frequency jumps into the cost function, and the penalty coefficient needs to be set empirically based on signal characteristics. In low-Earth orbit satellite communication scenarios, channel conditions and dynamic characteristics change drastically, making empirical parameter tuning not only cumbersome but also difficult to guarantee stability and optimality under various operating conditions.

[0006] Secondly, the fully connected state transition graph results in high complexity and a lack of robustness against interference. Existing methods use a fully connected state transition graph, meaning that all frequency states in each time frame must transition to all states in the next frame, leading to a computational complexity of O(n log n) per frame. ,in, This represents the number of frequency states. More importantly, this fully connected structure cannot provide structural immunity to cross-frequency trajectory jumps induced by co-channel interference (CCI). In shared frequency bands such as UHF, co-channel interference signals from adjacent orbital planes form competing frequency trajectories in the time spectrum. When the interference energy exceeds the target signal in some frames, the globally optimal path of existing algorithms is very likely to "jump" to the interference trajectory, especially in the high-leverage region at the edge of the observation window. Even jumps in a few frames can severely bias the slope estimation of the first-order polynomial fitting, leading to subsequent carrier loop lockout.

[0007] To address the aforementioned issues, recent research has attempted to introduce Bayesian filtering (such as particle filtering and extended Kalman filtering) to incorporate physical priors as state transition models for frequency tracking. While these methods inherently possess the ability to inject physical priors, their limitations include: the need for online recursive statistics, reliance on specific assumptions about noise distribution, and the ability to obtain only locally optimal solutions or sampled approximate solutions, failing to guarantee global optimality within a finite time domain.

[0008] In summary, how to structurally eliminate cross-frequency trajectory jumps, reduce computational complexity, and achieve robust parameter design without empirical parameter tuning, while ensuring global optimality and with as little prior information as possible, is a technical challenge that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0009] The purpose of this invention is to overcome the shortcomings of existing technologies and propose a method and system for extracting the Doppler frequency trajectory of low-Earth orbit satellites based on orbital constraint dynamic programming, using the maximum Doppler acceleration of satellite orbital mechanics. As a unique scalar prior, determine the maximum allowed jump of the inter-frame frequency bin index. With smoothing penalty coefficient The lower bound is simultaneously achieved through a hard constraint mechanism within the reachable domain of the same orbit: structurally excluding cross-frequency trajectory jump paths to obtain structural immunity to co-frequency interference; guaranteeing the global optimal solution through Bellman optimality; and reducing the search complexity of each frame from... Down to No experience or parameter tuning is required.

[0010] To achieve the above objectives, the technical solution specifically adopted by the present invention is as follows:

[0011] A method for extracting the Doppler frequency trajectory of low-Earth orbit satellites based on orbit-constrained dynamic programming includes:

[0012] Step 1: Obtain the time spectrum of the received signal. The time spectrum includes multiple time frames and multiple frequency indices. Each frequency index of each time frame corresponds to a time-frequency energy value.

[0013] Step 2: Using the maximum Doppler acceleration of satellite orbital mechanics as the only scalar prior, determine the maximum allowable number of jumps in the inter-frame frequency index based on the maximum Doppler acceleration, the time interval between adjacent time frames of the time spectrum, and the frequency resolution of the time spectrum.

[0014] Step 3: Determine the smoothing penalty coefficient based on the maximum Doppler acceleration and the maximum time-frequency energy value of a single frame in the time spectrum, and apply the smoothing penalty coefficient to perform a secondary penalty on the frequency index jump between adjacent time frames;

[0015] Step 4: Construct a cost function with hard constraints on the track reachability domain in the time spectrum. The cost function includes a path cumulative energy gain term and an inter-frame frequency jump quadratic penalty term weighted by the smoothing penalty coefficient.

[0016] Step 5: Under the hard constraint of the track reachability domain, use the Viterbi dynamic programming algorithm to solve for the frequency index path that maximizes the cost function, and use the solved frequency index path as the optimal Doppler frequency trajectory.

[0017] Preferably, the time spectrum is obtained by performing time-frequency analysis on the received signal.

[0018] Preferably, the maximum allowable number of jumps in the inter-frame frequency index is determined by dividing the product of the maximum Doppler acceleration and the time interval by the frequency resolution, and rounding the quotient up to obtain the integer as the maximum allowable number of jumps in the inter-frame frequency index.

[0019] Preferably, the smoothing penalty coefficient is determined as follows:

[0020] A constraint value is determined based on the maximum allowable number of jumps, and the constraint value is the maximum allowable number of jumps plus one.

[0021] Divide the maximum time-frequency energy value of a single frame by the square of the constraint value to obtain the lower bound of the penalty coefficient;

[0022] The lower bound of the penalty coefficient is multiplied by a preset engineering margin coefficient to obtain the smoothing penalty coefficient. The value of the smoothing penalty coefficient is such that the penalty cost corresponding to any inter-frame frequency index jump that exceeds the maximum allowed number of jumps is greater than the maximum time-frequency energy value of a single frame.

[0023] Preferably, step 4, before constructing the cost function with hard constraints on the orbital reachability domain, further includes a step of normalizing the time spectrum:

[0024] For each time frame of the time spectrum, the time-frequency energy value of each frequency index within the time frame is divided by the maximum time-frequency energy value within the time frame, so that the maximum time-frequency energy value of each time frame after normalization is always equal to 1.

[0025] Preferably, the hard constraint is that the absolute value of the frequency index jump between adjacent time frames does not exceed the maximum allowed number of jumps.

[0026] Preferably, the specific steps of using the Viterbi dynamic programming algorithm to solve for the frequency index path that maximizes the cost function include:

[0027] Initialization steps: For each frequency index of the first time frame, set the state value to the time-frequency energy value corresponding to that frequency index;

[0028] Forward recursion step: For each current frequency index of each current time frame starting from the second time frame, among all predecessor frequency indices that satisfy the hard constraint of the orbit reachability domain, select the predecessor frequency index that has the largest value after subtracting the smoothing penalty coefficient from the state value of the previous time frame and the product of the square of the frequency index jump between the current frames. Record the selected predecessor frequency index as the optimal predecessor index, and update the state value of the current frequency index to the time-frequency energy value of the current frequency index plus the largest value.

[0029] Termination step: Among all the frequency indices in the last time frame, select the frequency index with the largest state value as the optimal path endpoint;

[0030] Reverse backtracking step: Starting from the end point of the optimal path, backtrack frame by frame according to the optimal predecessor index recorded in each time frame until the first time frame is reached to obtain the optimal Doppler frequency trajectory.

[0031] Preferably, step 5 further includes a dynamic maximum allowable jump number adaptive expansion step:

[0032] In the forward recursion step, the maximum allowable number of jumps in subsequent time frames is adjusted in real time based on the estimated instantaneous frequency change rate of the first N0 time frames of the tracked frequency index path.

[0033] When the estimated instantaneous frequency change rate reaches a preset proportion of the nominal value, the maximum allowable number of jumps is increased; when the estimated instantaneous frequency change rate is lower than another preset proportion of the nominal value, the maximum allowable number of jumps is restored.

[0034] The nominal value is the maximum Doppler acceleration multiplied by the time interval and divided by the frequency resolution.

[0035] Preferably, step 5 further includes an edge adaptive trimming step:

[0036] For each time frame of the time spectrum, the time frames corresponding to the first quarter and the last quarter of the samples of each short-time Fourier transform window are removed, and the central region is retained as the effective observation window.

[0037] Preferably, step 5 also includes a full-bandwidth cold-start blind search mode: the region of interest of the time spectrum is extended to the full bandwidth range, the frequency index with the largest time-frequency energy value among all frequency indices is taken as the path starting point in the initialization step, and the maximum allowed number of jumps and the smoothing penalty coefficient are used in the forward recursion step.

[0038] The present invention also provides a system for implementing the above method, comprising:

[0039] The parameter determination module is used to determine the maximum allowable number of jumps in the inter-frame frequency index based on the maximum Doppler acceleration of satellite orbital mechanics as the only scalar prior, the time interval between adjacent time frames of the time spectrum, and the frequency resolution of the time spectrum, and to determine the smoothing penalty coefficient based on the maximum Doppler acceleration and the maximum time-frequency energy value of a single frame of the time spectrum.

[0040] The cost function construction module is used to construct a cost function with hard constraints on the track reachability domain on the time spectrum. The hard constraint is that the absolute value of the frequency index jump between adjacent time frames does not exceed the maximum allowed number of jumps. The cost function includes a path cumulative energy gain term and an inter-frame frequency jump quadratic penalty term weighted by the smoothing penalty coefficient.

[0041] The Viterbi solver module is used to solve for the frequency index path that maximizes the cost function under the hard constraint of the track reachability domain using the Viterbi dynamic programming algorithm, and the solved frequency index path is used as the optimal Doppler frequency trajectory.

[0042] This invention has the following characteristics and beneficial effects:

[0043] This invention introduces a hard constraint on the track reachability region in the state transition diagram (the absolute value of the inter-frame frequency index jump does not exceed...). This hard constraint structurally removes any cross-frequency trajectory jumps that exceed physical possibilities from the search feasible region in a pruning manner. Unlike existing technologies that rely solely on "soft penalties" to suppress jumps at the cost function level, this invention guarantees from the search space level that the optimal frequency trajectory cannot contain non-physical cross-frequency jumps, thus possessing a structural advantage in suppressing co-channel interference. The introduction of this hard constraint mechanism ensures that the global optimal solution is structurally consistent with the physically realizable Doppler frequency trajectory.

[0044] In this invention, the orbit reachability domain jump constraint With smoothing penalty coefficient All are composed of satellite orbital scalars Uniquely determined by a closed analytical expression:

[0045]

[0046] in The lower bound is precisely determined by the criterion that "the penalty for any illegal jump outside the orbital reachability region is strictly greater than the maximum energy gain obtainable in a single frame," ensuring that illegal paths are eliminated in the global optimization. The entire parameter derivation process relies solely on scalars. It achieves robust parameter design without the need for any ephemeris information or experience-based parameter tuning.

[0047] In this invention, the hard constraint of the orbital reachability region will be the first Frame to the The state transition evaluation range for each frequency bin extends from fully connected structures. The number of candidate precursors has been reduced to only those requiring evaluation. neighboring With 10 candidate predecessors, the search complexity per frame is reduced from 1000 to 10000. Down to The overall dynamic programming complexity is from Down to The complexity gain is directly determined by the pruning structure of the state transition diagram, and is independent of the specific operating conditions.

[0048] This invention employs a Viterbi-type forward recursion and backward backtracking to precisely solve the cost function under the hard constraint of the orbital reachability domain. Based on the Bellman optimality principle, the forward recursion retains the optimal sub-path among all feasible paths leading to the current state at each step, while the frequency-indexed path obtained through backward backtracking is the globally maximized path of the cost function under the orbital reachability domain constraint. This global optimality is guaranteed by the algorithm structure itself during the recursion process in each frame, requiring no approximation or iteration.

[0049] This invention extends the region of interest to the full bandwidth in the full-bandwidth blind search mode. ), utilizing hard constraints on the reachability region of the orbit to expand the search space from Compress to This makes the original complexity The computationally infeasible full-bandwidth Doppler blind search becomes engineering feasible. This improvement is directly determined by the reduced algorithm complexity. Attached Figure Description

[0050] Figure 1 This is a schematic diagram illustrating an application scenario of an embodiment of the present invention.

[0051] Figure 2 This is a flowchart illustrating the execution of the Doppler frequency trajectory extraction method based on trajectory-constrained dynamic programming (PC-DP) according to an embodiment of the present invention.

[0052] Figure 3 For embodiments of the present invention, the load-to-dryness ratio is... Signal-to-noise ratio Under the given conditions, the two-dimensional time-frequency energy spectrum obtained after FSST processing and the position marking diagram of the target Doppler frequency trajectory and the co-frequency interference Doppler frequency trajectory are shown.

[0053] Figure 4 The above is a time-series comparison curve of the error between the instantaneous frequency estimation and the true value obtained by the track-constrained dynamic programming method described in this invention and the traditional unconstrained frame-by-frame peak finding (argmax) method during the pilot band duration in an embodiment of this invention.

[0054] Figure 5 To test different load-to-dryness ratios under 500 Monte Carlo simulations. Under the given conditions, the Doppler rate of change (slope) estimation median error comparison curves obtained by the embodiments of the present invention, the unconstrained argmax method, and the weakly constrained dynamic programming method that does not reach the theoretical lower bound are respectively presented.

[0055] Figure 6 In order to decline in Rice Under 114Hz / s complete orbital dynamics Doppler and AWGN channel conditions, the embodiments of the present invention, in conjunction with a second-order Costas phase-locked loop, and four comparative synchronization schemes (traditional second-order phase-locked loop, Kaplan third-order phase-locked loop, single-segment FFT estimation plus phase-locked loop, and segmented FFT linear fitting plus phase-locked loop), yield end-to-end bit error rate performance curves as a function of signal-to-noise ratio.

[0056] Figure 7 To achieve the PC-DP (K=1, ...) of this invention under the denseFSST time grid, ) and unconstrained argmax and pure hard pruning ( A comparison curve of the full-frame extrapolated frequency RMSE as a function of carrier-to-interference ratio under a frequency bin interference spacing of 1.

[0057] Figure 8 To achieve the PC-DP (K=1, ...) of this invention under the denseFSST time grid, ) and unconstrained argmax and pure hard pruning ( Comparison curves of the full-frame extrapolated frequency RMSE as a function of carrier-to-interference ratio under two frequency bin interference intervals.

[0058] Figure 9 To achieve the PC-DP (K=1, ...) of this invention under the denseFSST time grid, ) and unconstrained argmax and pure hard pruning ( 0) Comparison curves of the full-frame extrapolated frequency RMSE as a function of carrier-to-interference ratio under the four frequency bin interference intervals. Detailed Implementation

[0059] The present invention will now be described in detail with reference to specific embodiments. These embodiments will help those skilled in the art to further understand the present invention, but do not limit the invention in any way. It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.

[0060] For ease of description, the mathematical symbols and their physical meanings used in the following embodiments are uniformly explained as follows:

[0061] in, This represents the discrete baseband pilot sequence after downsampling. Indicates the instantaneous Doppler frequency offset of the target satellite; Indicates the Doppler phase of the target satellite; Indicates the phase of co-channel interference signals; Represents the Rician fading coefficient; Indicates the amplitude of co-channel interference; This represents complex Gaussian white noise; Indicates the sampling period after downsampling; The sampling rate is represented by D; the I&D downsampling factor is represented by W; and FSST represents the short-time Fourier transform window length and number of samples. Indicates the shape parameters of the Kaiser window; Indicates the frequency bin width (frequency resolution); Indicates the interval between adjacent time frames in FSST; N represents the total number of time frames; M represents the total number of frequency bins within the region of interest; K represents the maximum allowed number of frequency index jumps between frames. Indicates the smoothing penalty coefficient; Indicates the engineering margin coefficient; This represents the maximum time-frequency energy value of a single frame (which is always 1 after normalization). Indicates the frequency index path of the nth frame; This represents the state value (the maximum cumulative cost to reach the m-th bin in the n-th frame). Represents the backtracking pointer (the optimal predecessor index of the m-th bin in the n-th frame); This represents the maximum Doppler acceleration in the satellite's orbital dynamics. Indicates chip rate; Indicates the length of the spreading code; Indicates the duration of the pilot signal; Represents the Rician fading factor; Indicates the carrier-to-interference ratio; SNR indicates the signal-to-noise ratio. This represents the equivalent noise bandwidth of the Costas phase-locked loop.

[0062] Example 1

[0063] This embodiment uses a low-Earth orbit (LEO) satellite IoT communication scenario as an example, considering an LEO satellite with an orbital altitude of 600 km and a carrier frequency of 400 MHz, whose maximum Doppler variation rate is 114 Hz / s. The receiver adopts a direct-sequence spread spectrum scheme with a spreading code length of... =127 chips, chip rate =25 kcps, after matched filtering and integral clearing (I&D) downsampling, the output sampling rate is... =25 kHz.

[0064] The overall application scenario of this embodiment is as follows: Figure 1 As shown, the low Earth orbit (LEO) satellite operates at an orbital velocity of approximately 7.5 km / s, generating a maximum Doppler variation rate of 114 Hz / s at its overhead transit. The ground-based satellite IoT terminal simultaneously receives downlink direct sequence spread spectrum (DSSS) signals from the target satellite and co-channel interference (C / I as low as 2 dB, frequency offset approximately 50 Hz) from satellites in adjacent orbital planes. After receiving the signal, it undergoes integrator-and-distribution (I&D) downsampling and FSST time-frequency transformation. The proposed orbit-constrained dynamic programming (PC-DP) core algorithm extracts a continuous and smooth Doppler frequency trajectory for open-loop pre-compensation and subsequent second-order Costas phase-locked loop carrier synchronization.

[0065] In this embodiment, the system parameters are as follows:

[0066] The orbital altitude is 600 km; the orbital velocity is 7.5 km / s; the carrier frequency is 400 MHz (UHF); and the maximum Doppler variation rate is... The frequency is 114 Hz / s; the BPSK spreading code is 127 chips; the chip rate is... 25 kcps; I&D downsampling factor D = 4; FSST window function Kaiser = β = 10, W = 2048; frequency bin resolution The Hz frequency is 12.2 Hz; the region of interest (ROI) is... Hz (M=13 bins); pilot duration is 0.5 s; channel model Rician is Co-frequency interference Doppler shift approximately 50 Hz; reachable domain boundary K is 1 bin; smoothing penalty coefficient. It is 5.0; Monte Carlo times It is 500.

[0067] The specific execution flow of the method in this embodiment is as follows.

[0068] Step 1: Obtain the time spectrum of the received signal

[0069] The original baseband signal is downsampled by integral scavenging (I&D) with a downsampling factor of D=4, which means taking the average of the non-overlapping blocks of four consecutive sampling points to obtain the downsampled sampling rate. =25 kHz discrete baseband pilot sequence. The computational complexity of this operation is O(1) per output point, and Parseval's theorem guarantees that it provides AWGN channel performance. The coherent gain. The downsampled pilot signal is modeled as:

[0070]

[0071] The Doppler phase of the target satellite is as follows:

[0072]

[0073] The sampling period is determined by downsampling; it satisfies the following condition at the top crossing time. ,in The value is the Doppler rate of change. for The Rician fading coefficient reflects the Ricean fading characteristics of the satellite-to-ground link. Co-channel interference amplitude is determined by the carrier-to-interference ratio C / I.

[0074]

[0075] The interference phase is independent of the target star, and its frequency is offset from the target star by approximately 50 Hz. The additive noise w[n] follows a complex Gaussian distribution. Its variance It is determined by the system signal-to-noise ratio (SNR).

[0076] The selection of the I&D downsampling factor D=4 is based on the following criteria. In this embodiment, the value of D=4 is uniquely self-consistent under the following four engineering constraints.

[0077] (i) Nyquist Criterion and Chip Synchronization Margin: The RF front-end of the low-orbit IoT terminal performs analog-to-digital conversion at a chip rate of 4 times that of R_c=25kcps, i.e., F_{s,RF}=100 kHz. This satisfies the Nyquist sampling criterion that the maximum spectrum occupancy of the direct-sequence spread spectrum signal is about 2R_c, and also provides a 4× oversampling engineering margin for chip phase estimation.

[0078] (ii) Matched Filter Equivalence: I&D achieves composite matched filtering by using the mean of non-overlapping blocks with D=4 adjacent sampling points, thus reducing the sampling rate to (One sampling point is output per chip), which is the theoretically optimal sampling configuration for a direct-sequence spread spectrum receiver.

[0079] (III) Coherent Integral Gain: According to Parseval's theorem, I&D provides [the following] in an additive white Gaussian noise channel. Coherent integral gain improves the equivalent signal-to-noise ratio from -15 dB to -9 dB, which is crucial for frequency trajectory extraction in deep noise scenarios.

[0080] (iv) Frequency resolution matching: Matching the target Doppler frequency offset to achieve the desired frequency axis resolution (when W=2048) for subsequent FSST. =12.2 Hz) falls within the engineering resolution range. A value greater than 4 will result in an insufficient sampling rate to retain the main lobe information of the spreading code; a value less than 4 will result in insufficient coherence gain of the matched filter and a doubling of the FFT complexity. In summary, D=4 is the only optimal engineering value.

[0081] downsampled pilot sequence Perform a Short-Time Fourier Transform (STFT) with a Kaiser window (shape parameter β=10, window length W=2048). A Kaiser window with β=10 corresponds to... A sidelobe suppression level of dB can effectively suppress spectral leakage between adjacent competing frequency trajectories. After calculating the instantaneous frequency analytical redistribution mapping, the frequency bin width is obtained:

[0082] Two-dimensional FSST energy spectrum matrix The time spectrum comprises multiple time frames and multiple frequency indices, with each frequency index of each time frame corresponding to a time-frequency energy value. Based on the parameters of this embodiment, the region of interest (ROI) is taken as follows: Hz, corresponding to the number of frequency bins .

[0083] The selection of the FSST window length W=2048 and the Kaiser shape parameter β=10 is based on the following criteria. In this embodiment, the values ​​(W=2048, β=10) are determined under the following four considerations. First layer—Frequency-resolvable constraint: [Requirement] ; Project acquisition =12.2 Hz (approximately) (4× safety margin), from The inverse solution yields W=2048 (which is exactly...) This facilitates FFT butterfly optimization.

[0084] The second layer—time grid convention: FSST uses a sample-by-sample sliding frame (hop=1 sample, dense time grid, no inter-frame extraction), with frame intervals... ,exist =0.5 s pilot number of frames The dense mesh provides the highest temporal resolution and smoothest path continuity. The window edge effect region has a fixed duration (approximately W / 2 sampling points) and is adaptively clipped away by subsequent edges, not accumulating as N increases. Since the complexity of PC-DP is... Linear with the number of frames N, a large number of frames N is computationally feasible; and this number of frames makes the smoothness penalty an upper bound. Much larger =5, guaranteeing It falls within the legal range.

[0085] The third layer—Heisenberg time-frequency product compromise: Near the uncertainty limit—doubling W results in a time-frequency product of approximately 0.5, triggering FSST redistribution mapping discontinuities; halving W results in… Increased to 24 Hz, interference spacing They are almost indistinguishable.

[0086] Fourth layer – Sidelobe suppression constraint supply dB sidelobe suppression level, main lobe The dB width is approximately 1.6 frequency bins. Taking a small β (e.g., β=0 for a rectangular window) will result in... dB sidelobes cause the sidelobe energy of the interference frequency trajectory to couple to the target bin; β is taken to be large If the main lobe is smeared with more than 2 bins, the target itself loses its distinguishability. The above four considerations jointly determine... .

[0087] To ensure that the target frequency trajectory does not overflow the ROI within the entire observation window, the selection of the ROI half-width should meet the following criteria: ,in The maximum absolute value of the target satellite's Doppler frequency offset within the entire pilot band is determined by orbital altitude, orbital velocity, carrier frequency, and minimum observation elevation angle. A 1.5 times margin is used to cover orbital element prediction residuals, satellite clock errors, and spectrum occupancy caused by co-channel interference. Under the parameters of this embodiment, the pilot band... ≈38 Hz (overhead time), therefore take It also takes into account the 50 Hz co-frequency interference offset to ensure that both frequency trajectories fall completely within the observation range.

[0088] The above FSST energy spectrum was normalized according to the following criteria: for each frame ,Will Divide by This yields a normalized time-frequency spectrum matrix where the maximum energy of each frame is always 1, which serves as input for subsequent dynamic programming. This normalization optimizes the subsequent smoothing penalty coefficient. The value of this energy level is independent of the absolute energy level and is determined solely by the energy competition ratio; if this step is omitted... The system must be scaled accordingly based on the absolute energy level of each frame, which will introduce implementation complexity.

[0089] Step 2: Determine the maximum allowed number of hops in the inter-frame frequency index.

[0090] Maximum Doppler acceleration in satellite orbital mechanics As a unique scalar prior, based on the maximum Doppler acceleration The time interval between adjacent time frames of the time spectrum and the frequency resolution of the time spectrum Determine the maximum allowed number of jumps K in the inter-frame frequency index.

[0091] Maximum Doppler acceleration in satellite orbital dynamics Calculated from orbital altitude, carrier frequency, and minimum observation elevation angle. For the parameters in this embodiment (h=600 km, =400 MHz), ≈114 Hz / s. The maximum allowed number of jumps in the inter-frame frequency index bin is:

[0092]

[0093] in The interval between adjacent time frames of FSST (in this embodiment) =40 μs), Frequency bin width (in this embodiment) ≈12.2 Hz).

[0094] The minimum physical time required for the Doppler frequency to span a complete frequency bin is:

[0095]

[0096] For any project, a reasonable FSST frame interval (including sample-by-sample sliding) =40 μs, and the commonly used 50%~99% overlap configuration), substituting into the definition, we can obtain the maximum allowable jump K=1 for the inter-frame frequency bin index. That is, between adjacent FSST time frames, the target Doppler frequency changes by a maximum of 1 frequency bin. This value is robust to commonly used engineering frame interval configurations. When The sufficiency of K=1 is demonstrated as follows: when the actual Doppler rate of change satisfies At that time, by The derivation shows that the actual inter-frame frequency bin jump variable is always ≤ K, so K=1 will not exclude any real path from the feasible region.

[0097] Step 3: Determine the smoothing penalty coefficient

[0098] Based on the maximum Doppler acceleration Maximum time-frequency energy value per frame of the time spectrum Determine the smoothing penalty coefficient A smoothing penalty coefficient is applied to apply a secondary penalty to the frequency index jump between adjacent time frames.

[0099] The design principle is: any illegal transitions beyond the orbital reach domain ( The penalty incurred must be strictly greater than the maximum energy gain that frame could potentially obtain, thus rendering it obsolete in the global optimization. This can be formally expressed as... .

[0100] In the normalized FSST energy spectrum ( Substituting K=1, we obtain the theoretical lower bound:

[0101]

[0102] Considering two factors that enhance interference peak energy in engineering implementation: first, interference energy approaches the target signal at low carrier-to-interference ratios; second, noise fluctuations cause the interference bin energy in local frames to momentarily exceed the target bin energy, an engineering margin coefficient is introduced. Covering instantaneous energy fluctuations up to approximately 13 dB, taking .final:

[0103]

[0104] That is, the smoothing penalty coefficient The value of is such that the penalty cost corresponding to any inter-frame frequency index jump exceeding the maximum allowed number of jumps K is borne by . Greater than the maximum time-frequency energy value of a single frame This ensures that illegal transition paths are structurally eliminated during global optimization.

[0105] Introducing engineering margin coefficient Two types of physical factors that instantaneously enhance the energy of interference peaks in coverage engineering implementation: (i) instantaneous load-to-interference ratio fluctuations—at low load-to-interference ratios In the dB range, the difference between the single-frame energy of the interference peak and the target peak is further compressed. It needs to compensate for its energy exceedance during the duration of 1-3 frames; (ii) the instantaneous noise peak—complex Gaussian noise obeys the square envelope after squared envelope. The distribution, with its trailing effect, causes the energy of the interference bin in individual frames to be momentarily higher than that of the target bin by about 8 to 13 dB. The quantitative basis for the value is as follows: In the normalized energy spectrum, the squared envelope of the single-frame noise follows... Distribution, its cumulative distribution function ,Depend on The solution is that the desired level for controlling the false alarm probability in a single frame is... The following are required In engineering =20 (corresponding to approximately 13 dB of transient interference energy exceeding the limit), in a single frame False alarm threshold Sufficient margin is reserved above ≈6.9, and with the K=1 hard constraint forming a double protection, it is sufficient to suppress false alarms caused by transient interference and noise peaks in the cross-frequency trajectory jump within the dense observation window to a negligible level. This derivation relies solely on scalars. No ephemeris, no online interference power estimation, and no parameter search are required.

[0106] Further explanation of the scale and scope of application of λ': (as described above) and All are defined in step 1, per frame energy normalization scale ( Therefore This refers to the engineering margin parameters at the energy scale of the normalized decision frame, rather than fixed constants in terms of absolute energy; if not normalized, It must be scaled accordingly based on the average absolute energy of each frame. The license range is [ ], of which the upper boundary The decision frame number N (total number of FSST time frames) and the budgeted number of allowed spans J are jointly determined, and must satisfy the following conditions: ≤ This ensures that the J legitimate bin transitions of the target frequency trajectory are not frozen due to excessive penalty. In this dense embodiment, when N≈12500, J=5, and K=1... =625, legitimate.

[0107] Step 4: Construct the cost function with hard constraints on the orbit reachability region

[0108] A cost function with a hard constraint on the track reachability domain is constructed on the time spectrum. The hard constraint is that the absolute value of the frequency index jump between adjacent time frames does not exceed the maximum allowed number of jumps K. The cost function includes a path cumulative energy gain term and an inter-frame frequency jump quadratic penalty term weighted by the smoothing penalty coefficient.

[0109] Define frequency index path Inter-frame transfer Under constraints | Within the feasible region where |≤K, construct the joint cost function and take the path that maximizes it as the optimal frequency trajectory:

[0110]

[0111] The first term on the right-hand side of the equation is the path accumulation energy gain term, driving the optimal path towards a high-energy frequency trajectory; the second term is the quadratic penalty term for frequency jumps between adjacent frames, weighted by a smoothing penalty coefficient λ'. and They have the same physical dimensions; K is the maximum allowed number of jumps in the inter-frame frequency bin index, constituting a hard constraint on the track reachability region, i.e., the absolute value of the frequency index jump between adjacent time frames. | No more than K.

[0112] The anti-co-frequency interference mechanism of this invention is achieved through the synergistic implementation of a K=1 track reachability domain hard constraint and a smoothing penalty coefficient λ'. Essentially, it is a deterministic cumulative penalty barrier that grows linearly with the interference frequency spacing. Let the interference frequency trajectory be d frequency bins apart from the target frequency trajectory. The K=1 hard constraint prohibits any frequency jumps with an amplitude greater than 1 bin within a single frame. Therefore, the optimal path cannot instantaneously cross to the interference trajectory within a single frame as in frame-by-frame peak finding (argmax); any crossing to the interference frequency trajectory can only be decomposed into continuous ⌈d⌉-step single-bin jumps. Due to the quadratic penalty term in the cost function, the cumulative smoothing penalty of this multi-step crossing path is at least... This constitutes a deterministic potential barrier that grows linearly with the spacing d; only when the instantaneous cumulative energy gain provided by the disturbance and noise during this crossing process exceeds When the path is crossed, the solution may be able to reach the global optimum; otherwise, it will be structurally suppressed.

[0113] In this embodiment, the Doppler frequency shift between the co-frequency interference on adjacent orbital planes and the target is approximately 50 Hz, corresponding to... accumulating momentum This is significantly higher than the peak noise energy gain per frame under Rician fading (typically less than 3), therefore, in the engineering operating load-to-interference ratio range... The internal structure provides robust structural suppression against cross-frequency trajectory jumps.

[0114] Step 5: Use the Viterbi dynamic programming algorithm to solve for the optimal Doppler frequency trajectory.

[0115] Under the hard constraint of the orbital reachability domain, the Viterbi dynamic programming algorithm is used to solve for the frequency index path that maximizes the cost function, and the solved frequency index path is taken as the optimal Doppler frequency trajectory.

[0116] The specific execution steps are as follows:

[0117] (5-1) Initialization steps: For each frequency index of the first time frame , state value Let this be the time-frequency energy value corresponding to this frequency index, i.e. Backtracking pointer .

[0118] (5-2) Forward recursion steps: For each current time frame starting from the second time frame... Each current frequency index All predecessor frequency indices that satisfy the hard constraints of the orbit reachability domain In the middle, select the state value of the previous time frame. Subtract the smoothing penalty coefficient Inter-frame frequency index jump The index of the largest predecessor frequency after the product of the squares of:

[0119]

[0120] The selected predecessor frequency index j is recorded as the optimal predecessor index, i.e. ; and the status value of the current frequency index Updated to the time-frequency energy value of the current frequency index. Add the aforementioned maximum value:

[0121]

[0122] (5-3) Termination Step: Select a state value from all frequency indices m in the last time frame N. The highest frequency index is used as the endpoint of the optimal path:

[0123]

[0124] (5-4) Backtracking steps: Starting from the endpoint of the optimal path Initially, backtrack frame by frame based on the optimal predecessor index recorded in each time frame, that is, for Take in sequence This continues until the first time frame is traced back to obtain the optimal frequency index path. This is taken as the optimal Doppler frequency trajectory.

[0125] The overall computational complexity of the above steps is: In this embodiment, with K=1 configured, the complexity is simplified to... Compared to unconstrained dynamic programming Approximately 4.3x speedup was achieved with M=13. The computational complexity per frame for each method is compared in Table 1.

[0126] Table 1 Comparison of computational complexity per frame

[0127]

[0128] The overall algorithm flow of the above steps is as follows: Figure 2 As shown, the key decision points are the "hard constraint of the orbital reachability domain" and the "Viterbi forward recursive evaluation range". ".

[0129] Step 6: Proof of Global Optimality

[0130] By the Bellman optimality principle, defined in | The maximum cumulative cost of all feasible paths to the m-th frequency bin in the n-th frame under the constraint |≤K Then any optimal path to (n,m) has its predecessor. The sub-path on must be reachable. The optimal path. Therefore, forward recursion:

[0131]

[0132] By strictly preserving global optimality at each step, the backtracking yields... The path represents the globally maximized cost function within the orbital reachability domain. The orbital reachability domain constraint does not violate optimality; it merely compresses the feasible domain by eliminating non-physical paths.

[0133] Step 7: Frequency Trajectory Post-processing and Carrier Synchronization Interface

[0134] The optimal frequency index sequence output in step 5 of this invention By multiplying by the frequency bin width Obtain the instantaneous Doppler frequency trajectory in time series form. After adaptively cropping the trajectory to remove the boundary effect region of the FSST window function, it can be fed into a first-order least squares polynomial fitting to obtain the fitted Doppler rate of change. With initial frequency offset After reconstructing the instantaneous phase of the entire frame through time integration, inverse open-loop pre-compensation is performed on the original long baseband signal with high dynamic Doppler. The pre-compensated signal is then despread by a matched filter and followed by a second-order Costas phase-locked loop to track the remaining slowly varying frequency offset, achieving continuous carrier locking and data demodulation at an SNR of -15 dB. When interfacing with the Costas loop, the residual frequency offset is limited to within ±20 Hz after PC-DP open-loop pre-compensation, and the equivalent noise bandwidth of the second-order Costas loop is [not specified]. Set to 20~50 Hz to meet the requirements. Traction and tracking requirements with a loss of less than 0.5 dB.

[0135] The specific implementation of edge adaptive pruning is as follows: perform edge adaptive pruning on the FSST output frame, remove the time frames corresponding to the W / 4 samples before and after each STFT window, and use the remaining central region as the effective observation window; the pruning ratio is determined by the main lobe energy concentration of the Kaiser window β=10, which can ensure that the amplitude distortion caused by energy leakage at the window edge is less than 1 dB.

[0136] Implementation effect verification

[0137] like Figure 3 As shown, in In the scenario of coherent co-frequency interference, the FSST energy spectrum is in Two clearly distinguishable energy concentration bands are observed: the target Doppler frequency trajectory (frequency decreases linearly with time, slope approximately -114 Hz / s) and the co-channel interference frequency trajectory (crossing the ROI within 0.2 s of the pilot start). The unconstrained frame-by-frame peak finding (argmax) method exhibits repeated cross-frequency trajectory jumps during interference crossing; The median slope error is within the range of 3 to 30. The maximum fluctuation does not exceed 2 Hz / s, proving that the method has a flat performance curve within this range. The PC-DP method described in this invention maintains continuous and smooth tracking of the target frequency trajectory throughout the entire observation window.

[0138] like Figure 4 As shown, in the same Under these conditions, the unconstrained argmax method exhibits an instantaneous pulse error with an amplitude of -50 Hz in the high lever edge region (first 0.2 s) at the start of the pilot signal, within the effective region. The slope of the first-order fit is relative truth value The deviation reached 86.6 Hz / s. The estimation error of the PC-DP method described in this invention remained consistently within... within, The slope of the first-order fit is The deviation from the true value is less than Fully working Within the phase-locked loop traction zone.

[0139] like Figure 5 As shown, under 500 independent Monte Carlo tests, the PC-DP method described in this invention achieves a median slope error of 3~4 Hz / s at C / I=2 dB, which is about 8 times more accurate than the unconstrained argmax method (about 30 Hz / s), and meets the accuracy threshold of 5 Hz / s (corresponding to the steady-state phase error of the phase-locked loop). The unconstrained argmax method must be up to... Only then can the same threshold be satisfied. Weakly constrained dynamic programming, below the theoretical lower bound, is employed. The improvement achieved relative to argmax is only marginal, which, from the opposite perspective, verifies the hypothesis proposed in this invention. Precision of the lower bound: Below This allows illegal activities even when the instantaneous noise gain exceeds the penalty cost. Jump, this is precisely the failure mode that the lower bound proposed in this invention aims to exclude. All methods in The convergence at that time proves that the gain of the present invention strictly originates from the structural suppression of deterministic co-frequency interference jumps, rather than additional filtering of random background noise.

[0140] like Figure 6 As shown, in Rice decay ( ), 114 Hz / s complete orbital dynamics Doppler, AWGN and Under a combined channel with coherent co-channel interference, the end-to-end bit error rate (BER) performance of the five comparative synchronization schemes varies significantly with the signal-to-noise ratio (SNR). Among them, the traditional second-order phase-locked loop (PLL), the Kaplan third-order PLL, and the single-segment FFT estimation plus PLL scheme cannot accommodate the 114 Hz / s Doppler change rate within the PLL tracking capability, resulting in a BER consistently maintained at approximately 0.5 (complete loss of lock-in). The segmented FFT linear fitting plus PLL scheme... Passing through Threshold, at The error floor is reached; however, in the low signal-to-noise ratio range, the FFT peak finding is contaminated by co-frequency interference, causing estimation bias and shifting the performance waterfall curve to the right. The PC-DP combined with a second-order Costas phase-locked loop described in this invention is unique in... The entire test range meets the requirements. ,exist The relative segmented FFT at the threshold has an equivalent signal-to-noise ratio gain of approximately 3 dB.

[0141] like Figures 7-9 As shown, in the dense time grid, The following comparison of the full-frame extrapolated RMSE for three interference spacings (1 / 2 / 4 frequency bins) shows that the anti-interference barrier... The interference spacing increases linearly with the spacing d: when the interference spacing is near the lower resolvable limit of FSST (1~2 bins, potential barrier) When the interference is low (at the level of noise jitter), the performance of this invention is similar to that of unconstrained argmax with no significant gain; when the spacing increases to the engineering-typical 4 bins (approximately 50 Hz, barrier)... When the gain is much higher than the peak gain of a single frame noise level, this invention significantly outperforms argmax, reducing the full-frame extrapolated RMSE by approximately 40% compared to argmax. This result is consistent with the cumulative penalty barrier analysis, indicating that the structural anti-interference capability of the present invention is fully utilized under typical engineering adjacent track surface interference spacing.

[0142] Example 2: Dynamic K-Adaptive Expansion Mode

[0143] The main difference between this embodiment and embodiment 1 is that the method of determining K in step 2 and the use of dynamic K adaptive expansion mode in the forward recursion process of step 5.

[0144] In engineering scenarios where prior orbital elements are inaccurate or satellite overpass trajectories have significant deviations, the sliding window length is used during the forward recursion process. Frame, based on the previous traced path Estimated instantaneous rate of change of frame The K value is adjusted in real time for subsequent time frames. Specifically:

[0145] When the estimated value reaches a preset percentage of the nominal value, K is increased to handle occasional over-constraint jumps; when the estimated value falls below another preset percentage of the nominal value, K is restored to its original value to maintain complexity gains. The nominal value is... .

[0146] In the specific implementation of this embodiment, when the estimated value When, increase K from 1 to 2; when the estimated value And if this continues for more than 100 frames, K will be restored to 1.

[0147] It should be noted that during the dynamic adjustment process, the smoothing penalty coefficient in step 3... It needs to be based on the adjusted K value. The criteria are redefined to ensure that the penalty for an illegal transition is always strictly greater than the maximum energy gain per frame. That is, when K increases from 1 to 2... It decreased from 0.25 to approximately 0.111. have to .

[0148] This dynamic adjustment mechanism can address engineering scenarios where the prior knowledge of satellite orbital elements is inaccurate or the overpass trajectory has a large deviation, while keeping the overall framework of the hard constraints of the orbital reachability domain unchanged.

[0149] Example 3: Full-bandwidth cold start blind search mode

[0150] The main difference between this embodiment and Embodiment 1 lies in the selection range of the region of interest (ROI) in step 1 and the method of determining the initial path starting point in the initialization step 5.

[0151] In the absence of any frequency prior information during terminal power-on or initial network access, the region of interest (ROI) of the time spectrum is expanded to the full bandwidth, i.e., the ROI is expanded to... Frequency bin number In the initialization step of step 5, the frequency index with the largest time-frequency energy value among all frequency indices is taken as the path starting point. That is, for each frequency index m in the first frame, the state value V(1,m) is set to the time-frequency energy value corresponding to that frequency index. In the forward recursion, the same maximum allowable number of jumps K=1 and smoothing penalty coefficient as in Example 1 are used. .

[0152] Using hard constraints on the reachability region of the orbit to expand the search space from Compression to O(3M) achieves the engineering feasibility of full-bandwidth Doppler blind search. In this embodiment, due to the significant increase in M ​​(from 13 to 2048), the complexity of unconstrained dynamic programming is O(3M). This is computationally infeasible on resource-constrained satellite IoT terminals. The complexity of the PC-DP method in this invention is... The complexity is reduced by about 683 times, making full-bandwidth blind search feasible instead of computationally infeasible.

[0153] Example 4: Parameter Configuration at Different Orbital Altitudes and Carrier Frequencies

[0154] This embodiment illustrates the parameter configuration method of the present invention under different orbital altitudes and carrier frequencies. Those skilled in the art can configure the parameters accordingly based on specific application scenarios and the parameter derivation principles of the present invention. The typical value ranges of the key parameters in this embodiment are as follows: Orbital Altitude carrier frequency (Compatible with UHF / L band satellite IoT links), signal-to-noise ratio (SNR) ≥ -15 dB, spreading code length 63~1023 chips, pilot length 0.2~1 s, downsampling factor FSST window length W = 1024~4096, Kaiser shape parameters Frequency bin width =6~25 Hz, range of interest M=11~33 bins, track reachability jump constraint K=1~3 (K=1 adapts to over-the-top high dynamics, K=2~3 adapts to low elevation angle medium-low dynamics), engineering margin coefficient Penalty coefficient (Calculated from track parameters, no real-time parameter adjustment required), Costas phase-locked loop equivalent noise bandwidth .

[0155] (a) Orbital altitude h = 300 km, carrier frequency (VHF band)

[0156] Maximum Doppler rate of change The orbital altitude, orbital velocity, and carrier frequency are all considered. For a low Earth orbit satellite with h=300km, the orbital velocity is approximately 7.73 km / s, and with a 137 MHz carrier frequency, If the FSST window length W = 1024, ,but Frame interval Substitute:

[0157] . ,Pick have to .

[0158] (ii) Orbital altitude h = 1500 km, carrier frequency =2.4 GHz (L / S band)

[0159] For a low-Earth orbit satellite with h=1500 km and an orbital velocity of approximately 7.12 km / s, at a carrier frequency of 2.4 GHz, If the FSST window length W = 4096, ,but Frame interval Substitute At this point, the number of state transition evaluations per frame is (2K+1)=5, which is 5 / M more complex than unconstrained dynamic programming. ,Pick have to .

[0160] (III) Low elevation angle, low dynamic scene (K=2 or K=3)

[0161] At low observation elevation angles (such as...) In this scenario, the absolute value of the Doppler rate of change is lower than the maximum value at the time of overshoot. In this case, the K value calculated according to the above formula may be greater than 1. This invention also applies. or The complexity-benefit analysis is given in Table 1: the number of state transition evaluations per frame is... M represents relatively unconstrained dynamic programming. In terms of complexity ratio, Under normal ROI conditions The time complexity ratio is (Approximately 4.3 times faster); At that time (Approximately 2.6 times faster); when K=3 is (Approximately 1.86x speedup). Under full-bandwidth cold start blind search mode. , The corresponding complexity ratios are 0.29% / 0.49% / 0.68%, and the speedup compared to unconstrained dynamic programming remains consistently above 146 times. Therefore, this invention... Significant reductions in state transition evaluation can be achieved throughout the entire engineering range, and the complexity gains do not depend on the specific value of K.

[0162] Example 5

[0163] This embodiment also provides a system for implementing the methods described in embodiments 1-4, including:

[0164] The parameter determination module is used to determine the maximum allowable number of jumps in the inter-frame frequency index based on the maximum Doppler acceleration of satellite orbital mechanics as the only scalar prior, the time interval between adjacent time frames of the time spectrum, and the frequency resolution of the time spectrum, and to determine the smoothing penalty coefficient based on the maximum Doppler acceleration and the maximum time-frequency energy value of a single frame of the time spectrum.

[0165] The cost function construction module is used to construct a cost function with hard constraints on the track reachability domain on the time spectrum. The hard constraint is that the absolute value of the frequency index jump between adjacent time frames does not exceed the maximum allowed number of jumps. The cost function includes a path cumulative energy gain term and an inter-frame frequency jump quadratic penalty term weighted by the smoothing penalty coefficient.

[0166] The Viterbi solver module is used to solve for the frequency index path that maximizes the cost function under the hard constraint of the track reachability domain using the Viterbi dynamic programming algorithm, and the solved frequency index path is used as the optimal Doppler frequency trajectory.

[0167] Comparative Example 1: Unconstrained Frame-by-Frame Peak Finding ( )method

[0168] To verify the technical effect of the present invention, Comparative Example 1 is set up. Comparative Example 1 uses the same received signal model, preprocessing method and FSST time spectrum analysis as Example 1, but adopts an unconstrained frame-by-frame peak finding method in the frequency trajectory extraction stage: for each time frame, the frequency bin with the largest FSST energy in the frame is independently selected as the frequency estimate of the frame, without applying any inter-frame continuity constraints or performing dynamic programming path optimization.

[0169] The mathematical model for Comparative Example 1 is:

[0170]

[0171] This method lacks any inter-frame smoothing mechanism and cannot utilize prior information on orbital mechanics to suppress cross-frequency trajectory jumps induced by co-frequency interference.

[0172] like Figure 4 As shown, in Under the given conditions, the unconstrained argmax method in Comparative Example 1 exhibited an instantaneous pulse error with an amplitude of -50 Hz in the high lever edge region (first 0.2 s) at the start of the pilot signal, within the effective region. The first-order fitting slope was -27.4 Hz / s, with a deviation of 86.6 Hz / s relative to the true value of -114 Hz / s. This deviation exceeded the pulling range of the Costas phase-locked loop, leading to the failure of subsequent carrier synchronization.

[0173] like Figure 5 As shown, in 500 independent Monte Carlo trials, Comparative Example 1... The median slope error at that time is approximately This is far higher than the accuracy threshold of 5 Hz / s; it must be reached Only then can the threshold be met.

[0174] This comparative example demonstrates that frame-by-frame independent peak finding lacks inter-frame continuity constraints, and cannot accurately extract continuous Doppler frequency trajectories under conditions of low signal-to-noise ratio and co-frequency interference.

[0175] Comparative Example 2: Using a weakly constrained dynamic programming method with a lower bound than the theoretical limit

[0176] To verify the smoothing penalty coefficient derived in this invention The necessity of a lower bound is established by setting up Comparative Example 2. Comparative Example 2 uses the exact same received signal model, preprocessing method, FSST time-spectrum analysis, and Viterbi-type dynamic programming algorithm framework as Example 1. The only difference is that the value of λ' is lower than the theoretical lower bound of 0.25 proposed in this invention (in Example 1). ).

[0177] The cost function for Comparative Example 2 is:

[0178]

[0179] because The algorithm still allows illegal operations even when the instantaneous noise or interference energy gain exceeds 0.2. Jump, i.e. The value of cannot guarantee the core principle that "the penalty for illegal transitions is strictly greater than the maximum energy gain per frame".

[0180] like Figure 5 As shown, in 500 independent Monte Carlo trials, the comparative example... exist The median slope error within the range only achieves marginal improvement relative to argmax. The slope error at that time is approximately Significantly higher than that of Embodiment 1 of the present invention. This, from the opposite perspective, verifies the implications of the present invention. Precision of the lower bound: Below the theoretical lower bound, illegal transitions are still allowed even when the instantaneous noise gain exceeds the penalty cost, which is precisely the failure mode that the lower bound proposed in this invention aims to exclude.

[0181] Comparative Example 3: Unconstrained Fully Connected Dynamic Programming Method

[0182] To verify the complexity gains and structural anti-interference capabilities of the present invention, Comparative Example 3 was set up. Comparative Example 3 uses the same received signal model, preprocessing method and FSST time-spectrum analysis as Example 1, but uses a typical unconstrained fully connected dynamic programming method in the prior art for frequency trajectory tracking (such as the method described in Kwok & Jones, IEEE TSP, 2000).

[0183] The cost function of this method is:

[0184]

[0185] in The smoothing penalty coefficient is empirically tuned (there is no analytical method to determine this coefficient in the existing technology, and manual parameter tuning is required).

[0186] The difference from Example 1 is that: (1) the state transition diagram is a fully connected structure, that is, each frequency bin in the nth frame performs state transitions to all M frequency bins in the (n+1)th frame, and the complexity of each frame is O(M²). When M=13, 169 transitions need to be evaluated per frame, and when M=2048, approximately 169 transitions need to be evaluated per frame. (2) There is no hard constraint on the reachability domain of the orbit (i.e., no | (with the constraint |≤K), inter-frame jumps of any magnitude are allowed in principle, subject only to a soft penalty term. To suppress; (3) It needs to be adjusted manually based on signal conditions and experience, without any theoretical basis.

[0187] In the decline of Rice ( In a complex channel with a combined interference pattern of 10 dB, 114 Hz / s complete orbital dynamics Doppler, AWGN, and C / I=2 dB coherent co-channel interference, the end-to-end bit error rate performance of Comparative Example 3 with a second-order Costas phase-locked loop is significantly inferior to that of this invention due to the lack of hard constraint protection and the occurrence of cross-frequency trajectory jumps during interference crossing. In terms of complexity, Comparative Example 3 requires evaluation per frame. This invention only requires evaluation, whereas the previous one involved a transfer of knowledge. For this type of transition (when K=1), a reduction of approximately 56 times in single-frame state transition evaluation is achieved (169 / 3≈56.3).

[0188] This comparative example demonstrates that the unconstrained fully connected structure, lacking hard constraints on the track reachability domain, cannot structurally eliminate cross-frequency trajectory jumps when co-frequency interference exists. Furthermore, its fully connected structure causes computational complexity to increase with M², making computation infeasible when M is large.

[0189] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely preferred examples and are not intended to limit the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for extracting the Doppler frequency trajectory of a low-Earth orbit satellite based on orbit-constrained dynamic programming, characterized in that, include: Step 1: Obtain the time spectrum of the received signal. The time spectrum includes multiple time frames and multiple frequency indices. Each frequency index of each time frame corresponds to a time-frequency energy value. Step 2: Using the maximum Doppler acceleration of satellite orbital mechanics as the only scalar prior, determine the maximum allowable number of jumps in the inter-frame frequency index based on the maximum Doppler acceleration, the time interval between adjacent time frames of the time spectrum, and the frequency resolution of the time spectrum. Step 3: Determine the smoothing penalty coefficient based on the maximum Doppler acceleration and the maximum time-frequency energy value of a single frame in the time spectrum, and apply the smoothing penalty coefficient to perform a secondary penalty on the frequency index jump between adjacent time frames; Step 4: Construct a cost function with hard constraints on the track reachability domain in the time spectrum. The cost function includes a path cumulative energy gain term and an inter-frame frequency jump quadratic penalty term weighted by the smoothing penalty coefficient. Step 5: Under the hard constraint of the track reachability domain, use the Viterbi dynamic programming algorithm to solve for the frequency index path that maximizes the cost function, and use the solved frequency index path as the optimal Doppler frequency trajectory.

2. The method according to claim 1, characterized in that, The time spectrum is obtained by performing time-frequency analysis on the received signal.

3. The method according to claim 1, characterized in that, The maximum allowable number of jumps in the inter-frame frequency index is determined as follows: the product of the maximum Doppler acceleration and the time interval is divided by the frequency resolution, and the resulting quotient is rounded up to the nearest integer, which is then used as the maximum allowable number of jumps in the inter-frame frequency index.

4. The method according to claim 1, characterized in that, The smoothing penalty coefficient is determined as follows: A constraint value is determined based on the maximum allowable number of jumps, and the constraint value is the maximum allowable number of jumps plus one. Divide the maximum time-frequency energy value of a single frame by the square of the constraint value to obtain the lower bound of the penalty coefficient; The lower bound of the penalty coefficient is multiplied by a preset engineering margin coefficient to obtain the smoothing penalty coefficient. The value of the smoothing penalty coefficient is such that the penalty cost corresponding to any inter-frame frequency index jump that exceeds the maximum allowed number of jumps is greater than the maximum time-frequency energy value of a single frame.

5. The method according to claim 1, characterized in that, In step 4, before constructing the cost function with hard constraints on the orbital reachability region, the step of normalizing the time spectrum is also included: For each time frame of the time spectrum, the time-frequency energy value of each frequency index within the time frame is divided by the maximum time-frequency energy value within the time frame, so that the maximum time-frequency energy value of each time frame after normalization is always equal to 1.

6. The method according to claim 4, characterized in that, The hard constraint is that the absolute value of the frequency index jump between adjacent time frames does not exceed the maximum allowed number of jumps.

7. The method according to claim 1, characterized in that, The specific steps for using the Viterbi dynamic programming algorithm to find the frequency index path that maximizes the cost function include: Initialization steps: For each frequency index of the first time frame, set the state value to the time-frequency energy value corresponding to that frequency index; Forward recursion step: For each current frequency index of each current time frame starting from the second time frame, among all predecessor frequency indices that satisfy the hard constraint of the orbit reachability domain, select the predecessor frequency index that has the largest value after subtracting the smoothing penalty coefficient from the state value of the previous time frame and the product of the square of the frequency index jump between the current frames. Record the selected predecessor frequency index as the optimal predecessor index, and update the state value of the current frequency index to the time-frequency energy value of the current frequency index plus the largest value. Termination step: Among all the frequency indices in the last time frame, select the frequency index with the largest state value as the optimal path endpoint; Reverse backtracking step: Starting from the end point of the optimal path, backtrack frame by frame according to the optimal predecessor index recorded in each time frame until the first time frame is reached to obtain the optimal Doppler frequency trajectory.

8. The method according to claim 1, characterized in that, Step 5 also includes a dynamic maximum allowed number of jumps adaptive expansion step: In the forward recursion step, the preceding path is indexed based on the already traced frequency. N The estimated instantaneous frequency change rate of 0 time frames is used to adjust the maximum allowable number of jumps in subsequent time frames in real time. When the estimated instantaneous frequency change rate reaches a preset proportion of the nominal value, the maximum allowable number of jumps is increased; when the estimated instantaneous frequency change rate is lower than another preset proportion of the nominal value, the maximum allowable number of jumps is restored. The nominal value is the maximum Doppler acceleration multiplied by the time interval and divided by the frequency resolution.

9. The method according to claim 1, characterized in that, Step 5 also includes an edge adaptive cropping step: For each time frame of the time spectrum, the time frames corresponding to the first quarter and the last quarter of the samples of each short-time Fourier transform window are removed, and the central region is retained as the effective observation window.

10. The method according to claim 1, characterized in that, Step 5 also includes a full-bandwidth cold-start blind search mode: the region of interest of the time spectrum is extended to the full bandwidth range, the frequency index with the largest time-frequency energy value among all frequency indices is taken as the starting point of the path in the initialization step, and the maximum allowed number of jumps and the smoothing penalty coefficient are used in the forward recursion step.

11. A system for implementing the method according to any one of claims 1-10, characterized in that, include: The parameter determination module is used to determine the maximum allowable number of jumps in the inter-frame frequency index based on the maximum Doppler acceleration of satellite orbital mechanics as the only scalar prior, the time interval between adjacent time frames of the time spectrum, and the frequency resolution of the time spectrum, and to determine the smoothing penalty coefficient based on the maximum Doppler acceleration and the maximum time-frequency energy value of a single frame of the time spectrum. The cost function construction module is used to construct a cost function with hard constraints on the track reachability domain on the time spectrum. The hard constraint is that the absolute value of the frequency index jump between adjacent time frames does not exceed the maximum allowed number of jumps. The cost function includes a path cumulative energy gain term and an inter-frame frequency jump quadratic penalty term weighted by the smoothing penalty coefficient. The Viterbi solver module is used to solve for the frequency index path that maximizes the cost function under the hard constraint of the track reachability domain using the Viterbi dynamic programming algorithm, and the solved frequency index path is used as the optimal Doppler frequency trajectory.

Images