Instantaneous frequency rapid extraction method suitable for high-dynamic star chain signal

Abstract

The invention discloses an instantaneous frequency rapid extraction method suitable for a high-dynamic star chain signal. The method comprises the following steps: realizing accurate tracking of a linear time-varying frequency based on an improved Rife interpolation algorithm; frequency hopping caused by multi-spectral line aliasing is eliminated through a time-frequency fuzzy region identification and compensation mechanism; a united Kalman filtering correction model is adopted, a signal domain observation value and an orbital domain theoretical Doppler value are fused, and the overall precision and robustness of frequency estimation are improved; and finally, according to a star chain satellite rapid matching and recognition algorithm, based on Doppler change rate consistency judgment and different-satellite Doppler differential matching, the satellite identity can still be accurately confirmed under a non-cooperative condition even if an error exists in the initial position of a receiver.

Description

Technical Field

[0001] This invention relates to the field of Starlink signal matching technology, and specifically to a method for rapid extraction of instantaneous frequency of highly dynamic Starlink signals. Background Technology

[0002] Starlink, as a representative of large-scale Low Earth Orbit (LEO) satellite communication systems, has shown application potential in non-cooperative navigation and positioning. However, its downlink signal is characterized by high dynamics, low signal-to-noise ratio, and multi-spectral line aliasing. Starlink signals operate in the Ku / Ka band, employing phased array antennas and frequency division multiple access (FDMA) technology. Its frequency domain structure exhibits a dual-layer characteristic, with narrowband beacon lines and wideband orthogonal frequency division multiplexed signals coexisting. The primary operating channel frequency band for Starlink signals is 10.7–12.7 GHz, with nine narrowband beacon lines distributed near each operating channel, and adjacent spectral lines spaced approximately 44 kHz apart. When multiple Starlink satellites pass overhead simultaneously, the beacon signals from different satellites intersect in the time and frequency domain, resulting in multi-spectral line aliasing. Due to the low orbital altitude and high speed of LEO satellites, their Doppler frequency offset exhibits significant linear time-varying characteristics, making traditional signal processing methods for medium and high Earth orbit satellites difficult to apply directly. Furthermore, under non-cooperative reception conditions, the signal power is low and the signal-to-noise ratio is often below 10dB, which further increases the difficulty of frequency estimation and satellite identification.

[0003] The existing technology is as follows: Traditional spectrum monitoring methods typically capture signals and perform spectrum analysis using broadband receivers. These methods are suitable for preliminary signal detection and macroscopic analysis, but they lack specific processing procedures for the specific structure of Starlink signals and cannot effectively solve the core problems of multi-spectral interference and satellite identification.

[0004] CN120750413A proposes a Ku-band Starlink downlink signal interference suppression method. This method involves segmenting the OFDM signal, performing a fast Fourier transform, and using the FCME algorithm to locate the interfering frequency for spectrum shifting. Although this method targets Starlink signals, its core purpose is interference suppression, and it does not address the issues of high-precision instantaneous frequency extraction and satellite identification.

[0005] CN120742375A presents a satellite navigation and positioning method that fuses observation parameters from a Global Navigation Satellite System (GNSS) and an Inertial Navigation System, and uses a positioning accuracy evaluation model to assess the reliability of the results. While this method involves multi-source data fusion, it primarily targets traditional GNSS signals and does not consider the unique high dynamic range and multi-spectral aliasing characteristics of Starlink signals. The disadvantages of existing technologies are as follows: First, high-dynamic signals are difficult to track stably. The high speed of Starlink satellites causes the Doppler frequency offset of the signal to exhibit strong linear time-varying characteristics. Traditional static frequency estimation methods cannot adapt to this dynamic change, and time-frequency analysis with a fixed window length also struggles to balance time and frequency resolution, resulting in severe time-spectrum spread, unstable frequency tracking, and large estimation errors.

[0006] Second, spectral line aliasing leads to large frequency estimation errors. Starlink signals use a frequency division multiple access (FDMA) structure, with its nine spectral lines distributed at 44kHz intervals. When multiple satellites pass overhead simultaneously, the beacon signals of different satellites intersect in the time and frequency domains, forming a complex spectral line aliasing phenomenon. Existing methods lack effective spectral line identification and alignment mechanisms, and frequency jumps easily occur in the aliasing region, significantly reducing estimation accuracy.

[0007] Third, there is a lack of effective non-cooperative satellite identification mechanisms. Under non-cooperative conditions, the receiver cannot know the key parameters of the signal in advance, making it difficult to effectively correlate the observed Doppler shift with satellite orbit data. There is a lack of feature extraction and matching algorithms to identify specific satellite sources from the received signal. Moreover, existing satellite identification methods mostly rely on known signal structures or cooperative modes, which cannot achieve accurate identification of Starlink satellites. Existing filtering methods are mostly based on single-source observations and lack a dual-source fusion mechanism in the signal and orbit domains, resulting in weak noise resistance and poor convergence. Summary of the Invention

[0008] The purpose of this invention is to provide a method for rapid extraction of instantaneous frequencies suitable for high-dynamic Starlink signals. It achieves accurate tracking of linear time-varying frequencies based on an improved Rife interpolation algorithm; eliminates frequency jumps caused by multi-spectral aliasing through a time-frequency ambiguity region identification and compensation mechanism; employs a joint Kalman filter correction model, fusing signal domain observations and orbital domain theoretical Doppler values ​​to improve the overall accuracy and robustness of frequency estimation; and finally, based on a Starlink satellite rapid matching and identification algorithm, using Doppler rate of change consistency judgment and alien Doppler differential matching, it can accurately confirm satellite identity even under non-cooperative conditions and when there are errors in the receiver's initial position.

[0009] The technical solution adopted in this invention is as follows: A method for rapid extraction of instantaneous frequency suitable for high dynamic starlink signals, the method comprising: Acquire Starlink signals; An adaptive sliding time-frequency extraction window is constructed, and the instantaneous frequency in the Starlink signal is extracted according to the window; Frequency compensation involves identifying time-frequency ambiguity regions by constructing a joint decision model of local time-frequency density gradient and Renyi entropy, determining the error of the time-frequency ambiguity regions using linear interpolation, and simultaneously compensating for the instantaneous frequency based on the error. Error suppression involves fusing the compensated instantaneous frequency with signal domain observations and orbital domain theoretical Doppler using a Kalman filter to output a high-precision Doppler frequency. Satellite matching involves constructing a satellite matching feature library using TLE files, inputting high-precision Doppler frequencies into the satellite matching feature library, determining the satellite identity through matching degree checks, and outputting the satellite ID.

[0010] Furthermore, the acquisition of the Starlink signal specifically includes: establishing an LFM equivalent model of the Starlink signal, the complex baseband expression of the Starlink signal being as follows: (1) In the formula, For signal amplitude, The initial frequency, The coefficient of frequency change. t For time, j It is the imaginary unit.

[0011] Furthermore, the instantaneous frequency of the Starlink signal includes the reference carrier frequency, spectral line shift, Doppler time-varying frequency, and noise fluctuation, as shown in the following equation: (2) In the formula, This is the baseband reference frequency after down-conversion. Multiples of the spectral line The spectral spacing of the beacon signal. For the initial Doppler frequency shift, For Doppler change rate, For noise fluctuations, the instantaneous frequency satisfies the following formula: (3) In the formula, For the estimated instantaneous frequency, This refers to the actual frequency.

[0012] Furthermore, the construction of the adaptive sliding time-frequency extraction window specifically includes: The adaptive sliding window mechanism is as follows: (4) In the formula, For window indexing, The sliding step size, For window length, s The representative signal indicates that the window length is dynamically adjusted as follows: (5) In the formula, As the baseline window length, As an adaptive factor, The maximum expected rate of change of the Doppler flow of the system. For the first Frequency of time, For the first Frequency of time.

[0013] Furthermore, the extraction of the instantaneous frequency specifically includes: The signal and noise components in the Starlink signal are separated as follows: Based on the ranking statistics estimation of the noise basis, let The accumulated power spectra are sorted in ascending order as follows: (6) Select the first sorted The noise characteristics are estimated at a single point, as shown in the following formula: (7) In the formula, The proportion of low-power point noise. Given the total number of data points, an adaptive detection threshold is constructed based on the statistical characteristics of noise to distinguish between signal and noise components, as shown in the following formula: (8) In the formula, Let k be the noise mean, and k be the time interval. For standard deviation estimation, The adaptive detection factor is designed as a function of the signal-to-noise ratio, as shown in the following equation: (9) In the formula, As the benchmark detection factor, To adjust the range, For reference signal-to-noise ratio, Control the slope of the transition zone. SNR Signal-to-noise ratio; The Rife generalized multispectral interpolation method is used to optimize the divided signal components, as shown in the following equation: (10) In the formula, P represents the number of adjacent spectral lines used, and P represents the spectral line amplitude. ，i For the first i Spectral lines, k 0 is the initial value for the iteration, and the weight coefficient is... , The spectral characteristic optimization design based on the window function is as follows: (11) In the formula, The total number of data points. For the Fourier transform of the window function, iterate over equation (10) to get the... The next iteration is as follows: (12) Iteration initial value The iteration termination condition is calculated from equation (10): (13) In the formula, The preset accuracy threshold, For signal components; The iterative process is convergent, and the convergence rate satisfies: (14) In the formula, This represents the true deviation. It is the Lipschitz constant; The instantaneous frequency is estimated as follows: (15) In the formula, r is the number of iterations. For window length, The total number of data points. The sampling rate.

[0014] Furthermore, the frequency compensation specifically includes: The construction of the joint decision model of local time-frequency density gradient and Renyi entropy specifically includes: Local time-frequency density calculation, within the time window With frequency window Below, for the time-frequency matrix Integrating the square of the amplitude yields the local time-frequency density as follows: (16) In the formula, The initial time, The initial frequency, The width of the time window. The width of the frequency window. For frequency; Renyi entropy is calculated as follows: (17) In the formula, To normalize the time-frequency probability, Let be the order of entropy; Joint decision-making is achieved by setting an adaptive threshold based on the characteristics of local time-frequency density gradient and Renyi entropy to realize joint decision-making in fuzzy regions. If the density gradient at a certain time-frequency point is greater than the density threshold... Furthermore, Renyi's entropy is greater than the entropy threshold. If the value is 0, then the point is determined to be a fuzzy region.

[0015] Furthermore, the setting of the adaptive threshold specifically includes: threshold and An adaptive method using statistical analysis of noise floor regions is employed. First, noise floor regions are identified by selecting areas in the time-frequency matrix whose amplitude is less than 0.1 times the peak amplitude. Then, noise floor characteristics are statistically analyzed, and the mean density gradient of these regions is calculated. with standard deviation and Renyi entropy mean with standard deviation ,according to and Obtain the density threshold and entropy threshold .

[0016] Furthermore, the compensation process specifically includes determining the fuzzy region and performing compensation using linear interpolation of effective points before and after the fuzzy region, assuming the time range of the fuzzy region is... The most recent valid moment before searching the fuzzy region The most recent valid time after the region Then the compensation frequency at any time t within the fuzzy region is: (18) In the formula, for The compensation frequency at any given time, for Frequency of time, for Frequency of time.

[0017] Furthermore, the error suppression specifically includes: Choosing frequency and rate of change of frequency as the state vectors, as shown in the following equation: (19) In the formula, for Frequency of time, for Rate of change of frequency at any given moment; Assume the state equations adopt a constant velocity model, as follows: (20) In the formula, the state transition matrix F is the transition matrix, and x(k-1) is the value at the previous time step. For sliding window intervals, The noise is process noise and follows a zero-mean Gaussian white noise distribution. Q is the process noise covariance matrix. The noise fluctuation at the frequency is greater than the rate of frequency change. The observation equations adopt a dual-source observation model, including signal domain observations and orbital domain observations, as follows: (twenty one) In the formula, The frequency observed in the signal domain, i.e., the frequency after Reif frequency compensation. For orbital domain observations, i.e., theoretical Doppler frequencies, For the observation matrix, To observe noise.

[0018] Federal Kalman filtering is used to fuse signal domain observations and orbital domain observations, including local filtering and global fusion: Local filtering applies Kalman filtering to each observation separately to obtain a local state estimate. , With covariance , Global fusion is achieved by weighted fusion of local estimates, with the weights determined by the inverse of the covariance matrix, as shown below: (twenty two) (twenty three) In the formula, For the estimation of the fused global state, This is the global covariance matrix.

[0019] Furthermore, the satellite matching specifically includes: Construct a satellite feature library to be matched: Set up the receiver The set of satellites that receive signals at all times is If the set contains M2 elements, then the receiver actually receives... pseudorange rate of the signal from the first satellite As shown in the following formula: (twenty four) In the formula, and For the first The true position and velocity of the satellite, and For the receiver's actual position and velocity, and The effects on the ionosphere and troposphere, At the speed of light, and For receiver and satellite clock drift; The first feature in the feature library to be matched The pseudorange rate corresponding to each satellite is as follows: (25) In the formula, , , For receiver location, For receiver speed, For satellite position, This represents the estimation error of the satellite's velocity; Identity matching: Coarse matching is performed using the Doppler rate of change consistency criterion. For the preprocessed signal, sliding window incoherent accumulation is employed to estimate the instantaneous frequency. , For a sliding window index, the rate of change for each second of data is calculated using a linear fit as follows: (26) In the formula, Multiply the number of sampling points by the total number of frames. , The sampling interval after downsampling; Theoretical rate of change calculation, theoretical Doppler sequence for candidate satellites The theoretical rate of change was obtained using the same block division and fitting method. According to the threshold Perform a matching decision, if If so, the satellite will be included in the fine-match candidate set; Construct an alien Doppler difference library for a fine-grained matching candidate set, and process the candidate set. Two satellites in China The theoretical Doppler difference is calculated as follows: (27) In the formula, These are the theoretical difference parameters. and These are the theoretical parameters for the i-th and j-th satellites, respectively. Calculate the corrected frequency and temporary satellite The measured difference is as follows: (28) In the formula, These are the measured differential parameters. These are the measured correction parameters. These are theoretical reference parameters; Perform a matching decision, if If the satellite is a match, its ID is output. These are the theoretical difference parameters. To determine the decision threshold, Doppler differential matching is performed on satellites that intersect in all time intervals until all satellites are traversed to obtain the final result.

[0020] In summary, due to the adoption of the above technical solution, the beneficial effects of this application are: 1. Compared with traditional static frequency estimation methods or STFT with a fixed window length, the improved Rife interpolation algorithm used in this invention can adaptively track the linear time-varying frequency of the starlink signal caused by high dynamic motion, fundamentally overcoming the shortcomings of spectrum spread in traditional methods.

[0021] 2. This invention stabilizes the instantaneous frequency estimation error within 10Hz through time-frequency ambiguity compensation and joint Kalman filtering correction, which is far lower than the estimation error of existing methods under the same high dynamic and low signal-to-noise ratio conditions, laying a solid foundation for subsequent high-precision navigation and positioning.

[0022] 3. This invention addresses the aliasing problem caused by the time-frequency crossover of signals from multiple Starlink satellites by proposing an innovative method that uses a joint decision based on local time-frequency density gradient and Renyi entropy. This method can accurately and adaptively identify time-frequency ambiguity regions and effectively correct frequency jumps caused by energy mixing through linear interpolation compensation. This mechanism solves the specific interference problem that existing spectrum monitoring methods cannot handle, ensuring the continuity and reliability of the frequency estimation curve throughout the entire observation period.

[0023] 4. This invention designs a two-stage satellite matching and identification mechanism. First, it uses the consistency of Doppler rate of change for rapid screening, which greatly narrows down the candidate range. Then, it uses the differential Doppler characteristics of alien stars to cleverly eliminate the influence of common errors such as ionospheric, tropospheric and clock errors. The method of this invention does not rely on the precise position of the receiver. Even with only an initial position estimate with errors, it can still achieve accurate confirmation of satellite identity, solving the core pain point of target identification under non-cooperative conditions. Attached Figure Description

[0024] Figure 1 This is a flowchart illustrating a method for rapid extraction of instantaneous frequencies for high-dynamic starlink signals according to the present invention. Figure 2 This is a schematic diagram of the instantaneous frequency extraction process in the method of the present invention; Figure 3 This is a schematic diagram of the satellite matching process in the method of the present invention. Detailed Implementation

[0025] The present invention will now be described in detail with reference to the accompanying drawings.

[0026] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0027] The method of this invention revolves around the entire process of Starlink signal reception, instantaneous frequency extraction, and satellite identification, achieving high-precision signal processing and satellite matching in non-cooperative environments through a modular architecture. For example... Figure 1 As shown, the method of the present invention includes six functional modules. The overall process follows the logic of signal preprocessing, accurate frequency extraction, error correction, satellite matching and result output. Each module works in concert to finally output high-precision Doppler frequency shift and satellite identification.

[0028] Example Starlink signal instantaneous frequency extraction and filtering correction: The instantaneous frequency extraction and filtering correction algorithm proposed in this invention adopts a modular pipeline design, such as... Figure 2 As shown, high-precision frequency estimation of Starlink signals is achieved through multi-stage processing. The entire process uses the baseband signal as input, which is the signal captured by the Starlink RF signal receiving and preprocessing module, down-converted, and bandpass filtered before being output. After key steps such as adaptive threshold detection, time-frequency analysis, ambiguity region compensation, and Kalman filter correction, the final output is high-precision Doppler frequency and satellite orbit information that meet navigation and positioning requirements, as detailed below: The acquisition of Starlink signals specifically includes: Since the Doppler frequency offset of Starlink signals changes linearly with time, its instantaneous frequency conforms to the characteristics of a linear frequency modulated (LFM) signal. First, an LFM equivalent model of the Starlink signal is established. Ignoring noise and spectral shift, the complex baseband expression of the Starlink signal is: (1) In the formula, For signal amplitude, The initial frequency, This is the frequency change rate coefficient.

[0029] The instantaneous frequency of the Starlink signal can be decomposed into the reference carrier frequency, spectral line shift, Doppler time-varying frequency, and noise fluctuation, as shown in the following equation: (2) In the formula, This is the baseband reference frequency after down-conversion. Multiples of the spectral line The spectral spacing of the beacon signal. For the initial Doppler frequency shift, For Doppler change rate, For noise fluctuations, the core objective of instantaneous frequency estimation is to extract the noise-fluctuating signal from the noisy signal. The estimation error should satisfy the following formula: (3) In the formula, For the estimated instantaneous frequency, This refers to the actual frequency.

[0030] This embodiment addresses the high dynamic characteristics of Starlink signals, where traditional fixed-length windows cannot simultaneously meet the requirements of time and frequency resolution. Therefore, this embodiment proposes an adaptive sliding window mechanism, specifically constructed as follows: (4) In the formula, For window indexing, The sliding step size, The window length is set as follows, and the window length is dynamically adjusted: (5) In the formula, As the baseline window length, As an adaptive factor, The maximum expected rate of change of the Doppler flow of the system. For the first Frequency of time, For the first The frequency of time; this design allows the window length to be automatically reduced to maintain time resolution when the frequency changes drastically, and the window length to be increased to improve the frequency estimation accuracy when the frequency changes slowly.

[0031] The extraction of instantaneous frequency specifically includes: Before frequency estimation, it is first necessary to separate the signal components and noise components in the starlink signal. Based on the ranking statistics, the noise basis is estimated, and then... The accumulated power spectra are sorted in ascending order as follows: (6) Select the first sorted The noise characteristics are estimated at a single point, as shown in the following formula: (7) In the formula, The proportion of low-power point noise. Given the total number of data points, an adaptive detection threshold is constructed based on the statistical characteristics of noise to distinguish between signal and noise components, as shown in the following formula: (8) In the formula, The noise mean. For standard deviation estimation, The adaptive detection factor is designed as a function of the signal-to-noise ratio: (9) In the formula, As the benchmark detection factor, To adjust the range, For reference signal-to-noise ratio, The slope of the transition region is controlled; this allows for the use of a conservative threshold to control false alarms at low signal-to-noise ratios and a lenient threshold to ensure detection probability at high signal-to-noise ratios.

[0032] The Rife generalized multi-spectral line interpolation method is used to optimize the divided signal components, making full use of the information from multiple spectral lines near the peak of the Starlink signal, as shown in the following equation: (10) In the formula, The weighting coefficient is the number of adjacent spectral lines used. , The spectral characteristic optimization design based on the window function is as follows: (11) In the formula, The Fourier transform of the window function ensures that the interpolation process fully considers the characteristics of the window function used in practice, improving the accuracy of the estimation. To further improve the estimation accuracy, an iterative refinement mechanism is introduced for the . Next iteration: (12) Iteration initial value The iteration termination condition is calculated from equation (10): (13) In the formula, Given a preset accuracy threshold, under appropriate conditions, this iterative process converges, and the convergence rate satisfies: (14) In the formula, This represents the true deviation. It is the Lipschitz constant; The final instantaneous frequency estimate is: (15) In the formula, This represents the number of iterations.

[0033] The instantaneous frequency extraction process provided in this embodiment has significant advantages in engineering implementation. Each processing module adopts a pipelined architecture, facilitating parallel processing on FPGA or DSP platforms. The adaptive mechanism reduces reliance on manual parameter adjustments and improves the system's automation level. Experimental results show that this process can achieve real-time processing on typical embedded platforms, meeting the practical application requirements of Starlink signal processing.

[0034] Frequency compensation specifically includes: When multiple Starlink satellites pass overhead simultaneously, their beacon signals intersect in the time and frequency domains, creating time-frequency ambiguity regions. This results in energy mixing in the time spectrum, causing abrupt changes in Reif frequency estimation. To address this, this embodiment designs a joint decision method using local time-frequency density gradient and Renyi entropy to accurately identify ambiguity regions. Furthermore, linear interpolation compensation is used to eliminate estimation errors caused by aliasing.

[0035] The core characteristics of time-frequency ambiguity regions are high energy concentration and large distribution uncertainty. Based on these two characteristics, a joint decision model of local time-frequency density gradient and Renyi entropy is constructed to identify ambiguity regions in three steps.

[0036] The first step is the calculation of local time-frequency density. Local time-frequency density is used to quantify the degree of energy concentration around a point in the time-frequency domain. Higher energy concentration corresponds to a higher density value. The calculation of local time-frequency density takes place within a time window. With frequency window Below, for the time-frequency matrix Integrating the square of the amplitude yields the local time-frequency density as follows: (16) In the formula, The initial time, The initial frequency, The width of the time window. This represents the width of the frequency window. The second step is to calculate the Renyi entropy. Ambiguous regions, due to the mixture of multiple signals and noise, have greater distribution uncertainty and thus higher Renyi entropy values. Unambiguous regions, with a single signal and less distribution uncertainty, correspond to lower Renyi entropy values. The Renyi entropy is expressed as follows: (17) In the formula, To normalize the time-frequency probability, Let be the order of entropy; The third step is joint decision-making. Based on the characteristics of local time-frequency density gradient and Renyi entropy, an adaptive threshold is set to achieve joint decision-making for fuzzy regions. If the density gradient at a certain time-frequency point is greater than the density threshold... Furthermore, Renyi's entropy is greater than the entropy threshold. If the value is 0, then the point is determined to be a fuzzy region.

[0037] Specifically, the setting of the adaptive threshold includes: threshold and An adaptive method using statistical analysis of the noise floor region is employed. First, the noise floor region is identified by selecting areas in the time-frequency matrix whose amplitude is less than 0.1 times the peak amplitude. At this point, there is no signal, only noise. Then, the noise floor characteristics are statistically analyzed, and the mean density gradient of the noise floor region is calculated. with standard deviation and Renyi entropy mean with standard deviation Finally, a threshold was set to ensure coverage of 99.7% of the noise floor area, based on... and Obtain the density threshold and entropy threshold .

[0038] After determining the fuzzy region, compensation is performed, specifically including: determining the fuzzy region, and using linear interpolation of effective points before and after it for compensation. Linear interpolation can avoid overfitting while ensuring accuracy. Let the time range of the fuzzy region be... The most recent valid moment before searching the fuzzy region The most recent valid time after the region Then the compensation frequency at any time t within the fuzzy region is: (18) In the formula, for The compensation frequency at any given time, for Frequency of time, for Frequency of time.

[0039] After Reif frequency estimation and time-frequency compensation, the frequency estimation error has been reduced to within 10Hz. However, in low signal-to-noise ratio scenarios, noise can still cause error fluctuations. Therefore, this embodiment integrates signal domain frequency estimation and orbital domain theoretical Doppler, designs a joint Kalman filter correction method, and constructs a two-dimensional state model of frequency and frequency change rate based on the high dynamic characteristics of Starlink signals. By integrating dual-source observations, error suppression is achieved. The specific error suppression includes: The state vector definition takes into account the linear change in Doppler frequency offset of the Starlink signal over a short period of time. Frequency and the rate of frequency change are chosen as the state vectors to track the time-varying frequency while avoiding excessive model complexity caused by higher-order derivatives. The state vector is defined as follows: (19) In the formula, for Frequency of time, for Rate of change of frequency at any given moment; Assuming the rate of change of frequency is constant over a very short time, the state equation adopts a constant velocity model, as follows: (20) In the formula, the state transition matrix ,in For sliding window intervals, The noise is process noise and follows a zero-mean Gaussian white noise distribution. Q is the process noise covariance matrix. The noise fluctuation at the frequency is greater than the rate of frequency change. The observation equations adopt a dual-source observation model, including signal domain observations and orbital domain observations, as follows: (twenty one) In the formula, The frequency observed in the signal domain, i.e., the frequency after Reif frequency compensation. For orbital domain observations, i.e., theoretical Doppler frequencies, For the observation matrix, To observe noise.

[0040] Due to the different reliability of dual-source observations, a federated Kalman filter is used to fuse signal domain and orbital domain observations, including local filtering and global fusion: Local filtering applies Kalman filtering to each observation separately to obtain a local state estimate. , With covariance , Global fusion is achieved by weighted fusion of local estimates, with the weights determined by the inverse of the covariance matrix, as shown below: (twenty two) (twenty three) In the formula, For the estimation of the fused global state, This is the global covariance matrix.

[0041] Local filtering can be updated independently, avoiding the collapse of the overall filter due to the failure of a single observation. For example, when orbital domain observations fail, only signal domain observations are used. Through inverse covariance weighting, observations with higher confidence are automatically assigned greater weights, resulting in a global covariance matrix. Always positive definite to ensure filter convergence.

[0042] Satellite matching specifically includes: like Figure 3As shown, satellite identification using motion information in signals first requires the construction of a feature library for matching. This feature library is constructed using publicly available TLE files. By comparing the measured Doppler rate of change with the theoretical Doppler rate of change, preliminary coarse matching results are obtained. Then, a database of alien Doppler differences among candidate satellites is constructed. The measured Doppler differences are compared with the theoretical Doppler differences. The satellite identity is finally confirmed through a matching degree check, and the satellite ID is output.

[0043] In navigation and positioning scenarios for high-speed moving platforms, the platform's precise location is unknown, but its initial position has a certain margin of error. Given, using The set of low-Earth orbit satellites that can be observed at the corresponding time can be obtained. The set has N elements. Based on the TLE file and the SGP4 model, the value of the set within a certain time period can be obtained. The trajectory of the elements in the middle, using these trajectories and Then a matching feature library can be constructed.

[0044] Construct a satellite feature library to be matched: Set the receiver time The set of satellites that received the signal is If the set contains M elements, then the receiver actually receives... pseudorange rate of the signal from the first satellite As shown in the following formula: (twenty four) In the formula, and For the first The true position and velocity of the satellite, and For the receiver's actual position and velocity, and The effects on the ionosphere and troposphere, At the speed of light, and For receiver and satellite clock drift; The first feature in the feature library to be matched The pseudorange rate corresponding to each satellite is as follows: (25) In the formula, , , For receiver location, For receiver speed, For satellite position, This represents the estimation error of the satellite's velocity; Identity matching: Coarse matching is performed using the Doppler rate of change consistency criterion. This coarse matching aims to narrow down the satellite candidate set, based on the consistency between the measured frequency rate of change and the theoretical Doppler rate of change, thus avoiding the high complexity of traditional trajectory matching. For the preprocessed signal, a sliding window incoherent accumulation method is used to estimate the instantaneous frequency. , For a sliding window index, the rate of change for each second of data is calculated using a linear fit as follows: (26) In the formula, Multiply the number of sampling points by the total number of frames. , The sampling interval after downsampling; Theoretical rate of change calculation, theoretical Doppler sequence for candidate satellites Using the same block division and fitting method, we obtained According to the threshold Perform a matching decision, if If the satellite is found to be a candidate for fine matching, then the coarse matching success rate is passed.

[0045] Due to the presence of the confusion curve, the results obtained are not accurate for all elements. Therefore, in the results obtained, the coarse matching results only retain the matching results of one satellite, that is, the matching result corresponding to the largest difference between the minimum and second smallest values. The matching results of other satellites are put into the next step, as follows: A differential library of alien Doppler data for a finely matched candidate set is constructed. Utilizing the approximate influence of the ionosphere / troposphere within the same field of view and the synchronization of clocks between satellites of the same system, non-orbital errors are eliminated, and the candidate set is then refined. Two satellites in China The theoretical Doppler difference is calculated as follows: (27) In the formula, These are the theoretical difference parameters. and The first i and j Theoretical parameters of the satellite; Calculate the corrected frequency and temporary satellite The measured difference is as follows: (28) In the formula, These are the measured differential parameters. These are the measured correction parameters. These are theoretical reference parameters; Perform a matching decision, if If the satellite is a match, its ID is output. These are the theoretical difference parameters. To determine the decision threshold, Doppler differential matching is performed on satellites that intersect in all time intervals until all satellites are traversed to obtain the final result.

[0046] This article uses specific embodiments to illustrate the principles and implementation methods of the present invention. The descriptions of the embodiments above are only for the purpose of helping to understand the method and core ideas of the present invention. It should be noted that those skilled in the art can make several improvements and modifications to the present invention without departing from the principles of the present invention, and these improvements and modifications also fall within the protection scope of the claims of the present invention.

Claims

1. A method for rapid extraction of instantaneous frequency suitable for high dynamic range starlink signals, characterized in that, The method includes: Acquire Starlink signals; An adaptive sliding time-frequency extraction window is constructed, and the instantaneous frequency in the Starlink signal is extracted according to the window; Frequency compensation involves identifying time-frequency ambiguity regions by constructing a joint decision model of local time-frequency density gradient and Renyi entropy, determining the error of the time-frequency ambiguity regions using linear interpolation, and simultaneously compensating for the instantaneous frequency based on the error. Error suppression involves fusing the compensated instantaneous frequency with signal domain observations and orbital domain theoretical Doppler using a Kalman filter to output a high-precision Doppler frequency. Satellite matching involves constructing a satellite matching feature library using TLE files, inputting high-precision Doppler frequencies into the satellite matching feature library, determining the satellite identity through matching degree checks, and outputting the satellite ID.

2. The method for rapid extraction of instantaneous frequency for high dynamic starlink signals according to claim 1, characterized in that, The acquisition of the Starlink signal specifically includes: establishing an LFM equivalent model of the Starlink signal, the complex baseband expression of the Starlink signal is as follows: （1） In the formula, For signal amplitude, The initial frequency, The coefficient of frequency change. t For time, j It is the imaginary unit.

3. The method for rapid extraction of instantaneous frequency for high dynamic starlink signals according to claim 1, characterized in that, The instantaneous frequency of the Starlink signal includes the reference carrier frequency, spectral line shift, Doppler time-varying frequency, and noise fluctuation, as shown in the following formula: （2） In the formula, This is the baseband reference frequency after down-conversion. Multiples of the spectral line The spectral spacing of the beacon signal. For the initial Doppler frequency shift, For Doppler change rate, For noise fluctuations, the instantaneous frequency satisfies the following formula: （3） In the formula, For the estimated instantaneous frequency, This refers to the actual frequency.

4. The method for rapid extraction of instantaneous frequency for high dynamic starlink signals according to claim 1, characterized in that, The construction of the adaptive sliding time-frequency extraction window specifically includes: The adaptive sliding window mechanism is as follows: （4） In the formula, For window indexing, The sliding step size, For window length, s The representative signal indicates that the window length is dynamically adjusted as follows: （5） In the formula, As the baseline window length, As an adaptive factor, The maximum expected rate of change of the Doppler flow of the system. For the first Frequency of time, For the first Frequency of time.

5. The method for rapid extraction of instantaneous frequency for high dynamic starlink signals according to claim 1, characterized in that, The extraction of the instantaneous frequency specifically includes: The signal and noise components in the Starlink signal are separated as follows: Based on the ranking statistics estimation of the noise basis, let The accumulated power spectra are sorted in ascending order as follows: （6） Select the first sorted The noise characteristics are estimated at a single point, as shown in the following formula: （7） In the formula, The proportion of low-power point noise. Given the total number of data points, an adaptive detection threshold is constructed based on the statistical characteristics of noise to distinguish between signal and noise components, as shown in the following formula: （8） In the formula, Let k be the noise mean, and k be the time interval. For standard deviation estimation, The adaptive detection factor is designed as a function of the signal-to-noise ratio, as shown in the following equation: （9） In the formula, As the benchmark detection factor, To adjust the range, For reference signal-to-noise ratio, Control the slope of the transition zone. SNR Signal-to-noise ratio; The Rife generalized multispectral interpolation method is used to optimize the divided signal components, as shown in the following equation: （10） In the formula, P represents the number of adjacent spectral lines used, and P represents the spectral line amplitude. ，i For the first i Spectral lines, k 0 is the initial value for the iteration, and the weight coefficient is... , The spectral characteristic optimization design based on the window function is as follows: （11） In the formula, The total number of data points. For the Fourier transform of the window function, iterate over equation (10) to get the... The next iteration is as follows: （12） Iteration initial value The iteration termination condition is calculated from equation (10): （13） In the formula, The preset accuracy threshold, For signal components; The iterative process is convergent, and the convergence rate satisfies: （14） In the formula, This represents the true deviation. It is the Lipschitz constant; The instantaneous frequency is estimated as follows: （15） In the formula, r is the number of iterations. For window length, The total number of data points. The sampling rate.

6. The method for rapid extraction of instantaneous frequency for high dynamic starlink signals according to claim 1, characterized in that, The frequency compensation specifically includes: The construction of the joint decision model of local time-frequency density gradient and Renyi entropy specifically includes: Local time-frequency density calculation, within the time window With frequency window Below, for the time-frequency matrix Integrating the square of the amplitude yields the local time-frequency density as follows: （16） In the formula, The initial time, The initial frequency, The width of the time window. The width of the frequency window. For frequency; Renyi entropy is calculated as follows: （17） In the formula, To normalize the time-frequency probability, Let be the order of entropy; Joint decision-making is achieved by setting an adaptive threshold based on the characteristics of local time-frequency density gradient and Renyi entropy to realize joint decision-making in fuzzy regions. If the density gradient at a certain time-frequency point is greater than the density threshold... Furthermore, Renyi's entropy is greater than the entropy threshold. If the value is 0, then the point is determined to be a fuzzy region.

7. The method for rapid extraction of instantaneous frequency for high dynamic starlink signals according to claim 6, characterized in that, The setting of the adaptive threshold specifically includes: threshold and An adaptive method using statistical analysis of noise floor regions is employed. First, noise floor regions are identified by selecting areas in the time-frequency matrix whose amplitude is less than 0.1 times the peak amplitude. Then, noise floor characteristics are statistically analyzed, and the mean density gradient of these regions is calculated. with standard deviation and Renyi entropy mean with standard deviation ,according to and Obtain the density threshold and entropy threshold .

8. The method for rapid extraction of instantaneous frequency for high dynamic starlink signals according to claim 6, characterized in that, The compensation process specifically includes determining the fuzzy region and performing compensation using linear interpolation of effective points before and after the fuzzy region, where the time range of the fuzzy region is set as follows: The most recent valid moment before searching the fuzzy region The most recent valid time after the region Then the compensation frequency at any time t within the fuzzy region is: ）（18） In the formula, for The compensation frequency at any given time, for Frequency of time, for Frequency of time.

9. A method for rapid extraction of instantaneous frequency for high dynamic starlink signals according to claim 7, characterized in that, The error suppression specifically includes: Choosing frequency and rate of change of frequency as the state vectors, as shown in the following equation: （19） In the formula, for Frequency of time, for Rate of change of frequency at any given moment; Assume the state equations adopt a constant velocity model, as follows: （20） In the formula, the state transition matrix F is the transition matrix, and x(k-1) is the value at the previous time step. For sliding window intervals, The noise is process noise and follows a zero-mean Gaussian white noise distribution. Q is the process noise covariance matrix. The noise fluctuation at the frequency is greater than the rate of frequency change. The observation equations adopt a dual-source observation model, including signal domain observations and orbital domain observations, as follows: （21） In the formula, The frequency observed in the signal domain, i.e., the frequency after Reif frequency compensation. For orbital domain observations, i.e., theoretical Doppler frequencies, For the observation matrix, To observe noise; Federal Kalman filtering is used to fuse signal domain observations and orbital domain observations, including local filtering and global fusion: Local filtering applies Kalman filtering to each observation separately to obtain a local state estimate. , With covariance , Global fusion is achieved by weighted fusion of local estimates, with the weights determined by the inverse of the covariance matrix, as shown below: （22） （23） In the formula, For the estimation of the fused global state, This is the global covariance matrix.

10. The method for rapid extraction of instantaneous frequency for high dynamic starlink signals according to claim 1, characterized in that, The satellite matching specifically includes: Construct a satellite feature library to be matched: Set up the receiver The set of satellites that receive signals at all times is If the set contains M2 elements, then the receiver actually receives... pseudorange rate of the signal from the first satellite As shown in the following formula: （24） In the formula, and For the first The true position and velocity of the satellite, and For the receiver's actual position and velocity, and The effects on the ionosphere and troposphere, At the speed of light, and For receiver and satellite clock drift; The first feature in the feature library to be matched The pseudorange rate corresponding to each satellite is as follows: （25） In the formula, , , For receiver location, For receiver speed, For satellite position, This represents the estimation error of the satellite's velocity; Identity matching: Coarse matching is performed using the Doppler rate of change consistency criterion. For the preprocessed signal, sliding window incoherent accumulation is employed to estimate the instantaneous frequency. , For a sliding window index, the rate of change for each second of data is calculated using a linear fit as follows: （26） In the formula, Multiply the number of sampling points by the total number of frames. , The sampling interval after downsampling; Theoretical rate of change calculation, theoretical Doppler sequence for candidate satellites The theoretical rate of change was obtained using the same block division and fitting method. According to the threshold Perform a matching decision, if If so, the satellite will be included in the fine-match candidate set; Construct an alien Doppler difference library for a fine-grained matching candidate set, and process the candidate set. Two satellites in China The theoretical Doppler difference is calculated as follows: （27） In the formula, These are the theoretical difference parameters. and These are the theoretical parameters for the i-th and j-th satellites, respectively. Calculate the corrected frequency and temporary satellite The measured difference is as follows: （28） In the formula, These are the measured differential parameters. These are the measured correction parameters. These are theoretical reference parameters; Perform a matching decision, if If the satellite is a match, its ID is output. These are the theoretical difference parameters. To determine the decision threshold, Doppler differential matching is performed on satellites that intersect in all time intervals until all satellites are traversed to obtain the final result.

Images