A Dynamic Atomic Gravity Interference Fringe Noise Suppression Method Based on Time-Frequency Peak Filtering

Abstract

This application belongs to the field of precision measurement and signal processing technology, and discloses a dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering. The method includes: phase calculation and unwinding; endpoint expansion of a continuous phase sequence; scaling the extended phase sequence to obtain a normalized signal; frequency modulation of the normalized signal, encoding the frequency modulation into the instantaneous frequency of a unit amplitude analytic signal, to obtain an analytic signal; PWVD transformation of the analytic signal to obtain a time-frequency peak map; finding the maximum frequency axis value at each moment from the time-frequency peak map and converting it into an instantaneous frequency to obtain an instantaneous frequency estimation signal; scaling the instantaneous frequency estimation signal and removing the extended endpoints, then calculating the interference fringe signal using cosine trigonometric functions; and calculating gravity based on the interference fringe signal. This method can recover phase information from a single experiment without requiring extensive superimposed sampling.

Description

Technical Field

[0001] This application belongs to the field of precision measurement and signal processing technology, specifically a dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering. Background Technology

[0002] Currently, atomic gravimeters generally employ a least-squares fitting method based on cosine functions to process data from individual interference fringes. This method can estimate gravity values ​​well under static laboratory conditions, but its applicability under dynamic conditions is poor, and its real-time performance is insufficient.

[0003] Under dynamic conditions, noise suppression is crucial for improving measurement performance, and the processing of dynamic absolute gravity data is a key factor. Besides improving the overall stability and anti-interference capabilities of the atomic gravimeter through hardware improvements, noise suppression can also be addressed from an algorithmic perspective. A common data analysis algorithm for atomic gravimeters involves fitting a cosine function to each interference fringe using the least squares method. This minimizes the variance between the estimated and observed values ​​for each fringe group, thus obtaining an estimate of the gravity value. This method is suitable for laboratory conditions, but its applicability and real-time performance are significantly reduced for mobile measurements. Existing methods suffer from the following problems: Least squares fitting: difficult to process in real-time under dynamic conditions, and lacks sufficient noise resistance. Combined algorithms with EKF: require additional hardware (accelerometers), increasing system complexity and cost. EKF-based methods: rely on a large amount of prior experimental data, resulting in poor portability and versatility. Other existing denoising algorithms: mostly target specific noises (spike noise, power supply harmonics), with limited effectiveness in suppressing random noise. Summary of the Invention

[0004] This application provides a dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering, which solves the problem of difficulty in real-time processing and insufficient noise resistance for moving measurements under dynamic conditions.

[0005] A dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering, according to an embodiment of this application, includes:

[0006] Phase calculation and unwinding, the phase calculation and unwinding includes: performing a preliminary phase estimate on the measurement data of atomic interference fringes, and performing phase unwinding on the preliminary phase estimate to obtain a continuous phase sequence;

[0007] Endpoint extension of a continuous phase sequence;

[0008] The normalized signal is obtained by scaling the continuous phase sequence after endpoint expansion.

[0009] The normalized signal is frequency modulated, and the frequency modulation is encoded into the instantaneous frequency of the unit amplitude analytic signal to obtain the analytic signal.

[0010] Perform PWVD transformation on the analytical signal to obtain the time-frequency peak diagram;

[0011] Find the maximum value of the frequency axis at each moment from the time-frequency peak plot, and convert it into the instantaneous frequency to obtain the instantaneous frequency estimation signal;

[0012] The instantaneous frequency estimation signal is scaled and the extended endpoints are removed. The interference fringe signal is then calculated using the cosine trigonometric function.

[0013] Gravity is calculated based on the interference fringe signal.

[0014] Furthermore, a preliminary phase estimate is performed on the measurement data of the atomic interference fringes, including obtaining a preliminary phase estimate using inverse trigonometric functions, expressed by the formula:

[0015] ,

[0016] in, For preliminary phase estimation, A represents the actual value of the atomic transition probability, A is the center value, and B is the contrast ratio.

[0017] Furthermore, phase dewinding is performed on the preliminary phase estimate to obtain a continuous phase sequence, including: calculating arbitrary phases using a phase dewinding algorithm. Phase value at time :

[0018] ,

[0019] Where Z represents the set of integers, including positive integers, zero, and negative integers. The value of makes the phase difference between adjacent points satisfy:

[0020]

[0021] for Phase of time, for Phase of time, This is a preliminary estimate of the phase.

[0022] Furthermore, endpoint extension of the continuous phase sequence is performed using the following formula:

[0023] ,

[0024] Where the length is M, x (n) is the extended signal with a length of 2p+M, x(1) represents the first number of the unextended signal, x(M) represents the Mth number of the unextended signal, x(np) represents the npth number of the original signal, M is the length of the phase sequence of the unextended signal, p is the number of points to be extended, p≥L+1, and L is the pseudo-Wegener distribution window.

[0025] Furthermore, by scaling the continuous phase sequence after endpoint expansion, a normalized signal is obtained, including:

[0026] ,

[0027] in, The signal after normalization represent The minimum value in, represent The maximum value in, This represents the expanded signal. b are empirical values.

[0028] Furthermore, the analytic signal is:

[0029] ,

[0030] in, It is an integral variable. The representative will arrive Summation, Indicates an analytical signal. The integral variable is The normalized signal.

[0031] Furthermore, the analytic signal is subjected to PWVD transformation using the following formula:

[0032] ,

[0033] in: The result is after PWVD transformation, where w(l) is the window function, k is the frequency sampling point, and the symbol * indicates conjugate. For a sequence length of The analytical signal, Indicates that the sequence length is The conjugate of the analytic signal.

[0034] Furthermore, the maximum frequency value at each moment is found from the time-frequency peak plot and converted into an instantaneous frequency to obtain the instantaneous frequency estimation signal, including:

[0035] Calculate the maximum value of the frequency axis: , Indicates the maximum value of the frequency axis;

[0036] Convert the maximum value of the frequency axis to the instantaneous frequency to obtain the instantaneous frequency estimation signal:

[0037] ,

[0038] in, This represents the number of frequency sampling points. The sampling rate.

[0039] Compared with the prior art, the advantages of this application are as follows: the method of this application can recover phase information from a single experiment without a large number of superimposed samples, which is particularly suitable for dynamic environments (such as aerial, shipborne and vehicle-mounted measurements), and significantly suppresses random noise while ensuring real-time performance.

[0040] Compared with the least squares fitting method, it does not require assuming a fringe function model. Phase information can be recovered through phase unwinding and TFPF processing, resulting in stronger real-time performance and noise resistance.

[0041] Compared with the extended Kalman filter (EKF), it does not rely on prior data and directly recovers the phase from the single interference fringe signal, which is more computationally efficient and more suitable for dynamic scenarios.

[0042] Compared to an auxiliary accelerometer, no additional sensors are added, ensuring the system's lightweight design and stability.

[0043] Compared with wavelet / low-pass filtering, by extracting the peak value in PWVD, an unbiased estimate of the instantaneous frequency is obtained without destroying the original phase information. Attached Figure Description

[0044] Figure 1 A flowchart illustrating the method provided in the embodiments of this application;

[0045] Figure 2The following are the results of the implementation process of the embodiments of this application, where (a) is the power spectral density of the simulated added vibration noise; (b) the vibration signal in (a) is used as noise and superimposed on the simulated interference fringes to construct the simulated fringe signal affected by vibration; (c) the phase information of the atomic interference fringes is obtained after phase extraction and deconvolution of the simulated fringe signal and elimination of π jump; (d) the result of normalizing the demodulated signal; (e) the corresponding time-frequency diagram obtained by converting the deconvolutioned signal into a frequency modulation (FM) analytical signal and analyzing it through pseudo-Wigner-Vell distribution (PWVD); (f) the signal diagram of peak extraction and instantaneous frequency estimation, where the signal energy is concentrated near the instantaneous frequency (IF), and random noise is presented as a dispersed component; (g) the effect diagram of inverse scaling and endpoint truncation, where the unbiased estimate of the signal can be obtained by extracting the peak value in the time-frequency diagram, and the phase information of the corresponding atomic gravity measurement interference fringes is recovered after inverse scaling; (h) the effect diagram of the method of this application before and after recovering the interference fringes. Detailed Implementation

[0046] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0047] This application proposes a dynamic atomic gravity interference fringe noise suppression method based on Time-Frequency Peak Filtering (TFPF). The core idea is to extract and eliminate the phase information in the atomic interference fringes through phase calculation and unwinding. A phase transition is used to maintain phase continuity throughout the measurement process. Frequency modulation coding is employed to convert the unwound signal into a frequency-modulated (FM) analytical signal. Subsequently, pseudo-Wigner-Ville distribution (PWVD) analysis is performed on the analytical signal to obtain its time-frequency distribution. In this distribution, signal energy is concentrated around the instantaneous frequency (IF), while random noise is more discrete. By extracting peak values, an unbiased estimate of the signal can be obtained, effectively suppressing random noise. Finally, the filtered instantaneous frequency signal is inversely scaled and its phase recovered to obtain the interference fringe phase information corresponding to the atomic gravity measurement.

[0048] Compared with the traditional least squares fitting method, the method in this application can recover phase information from a single experiment without a large number of superimposed samples, making it particularly suitable for dynamic environments (such as aerial, shipborne and vehicle-mounted measurements), and significantly suppressing random noise while ensuring real-time performance.

[0049] See Figure 2 As shown in (a) of this application embodiment, a vibration-contaminated fringe signal is constructed for noise suppression using the method of this application embodiment. The ground vibration power spectral density measured by the seismograph serves as the basis for noise generation and is superimposed on the simulated interference fringes to construct the vibration-contaminated fringe signal. The vibration-contaminated fringe signal is as follows: Figure 2 As shown in (b) of the diagram.

[0050] See Figure 1 As shown, this application provides a dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering, comprising:

[0051] S1, Phase Calculation and Unwinding, is used to perform preliminary phase estimation on the measurement data of atomic interference fringes, and to perform phase unwinding on the preliminary phase estimation to obtain a continuous phase sequence.

[0052] Suppose that the measured atomic interference fringes represent the actual values ​​of the atomic transition probabilities. for:

[0053] ;

[0054] The center value A and contrast B can both be obtained through fitting. However, in reality, phase noise is introduced by various components of the gravimeter. Influence, atomic phase value The actual value of the atomic phase should be obtained by adding noise to the phase, and is expressed as:

[0055] ,

[0056] in, Where is the effective wave vector, g is the gravitational acceleration, and T is the Raman pulse time interval. For sweep frequency rate, The actual value of the atomic phase, This represents the total noise of the atomic phase.

[0057] in, The phase noise introduced by the vibration of the reference mirror, hereinafter referred to as vibration phase noise, The sum of phase noise introduced by other noise sources, hereinafter referred to as other phase noise.

[0058] Preliminary phase estimation It can be obtained from inverse trigonometric functions:

[0059] ,

[0060] in, For preliminary phase estimation, The actual values ​​of atomic transition probabilities and atomic phases. There is a π-jump problem, therefore a phase unwinding algorithm is used to calculate... Phase value at time :

[0061] ,

[0062] Where Z represents the set of integers, including positive integers, zero, and negative integers. The value of makes the phase difference between adjacent points satisfy:

[0063] ,

[0064] in, for Phase of time, for The phase at each moment is used to obtain a continuous phase sequence through this step. For example... Figure 2 As shown in (c), phase information of atomic interference fringes is obtained while eliminating π jumps.

[0065] S2, performing endpoint expansion on a continuous phase sequence:

[0066] In time-frequency analysis, the window boundary effect can lead to energy leakage. To avoid signal distortion at both ends, endpoint extension is applied to the continuous phase sequence:

[0067] ,

[0068] Where x(n) is the unextended signal (of length M), x (n) The extended signal (length 2p+M), x(1) represents the first number of the unextended signal, x(M) represents the Mth number of the original signal, and x(np) represents the npth number of the unextended signal. M is the length of the original phase sequence, p is the number of points to be extended (p≥L+1), and L is the pseudo-Wegener distribution window, which is the window function. You can choose from options such as rectangular windows or Hanning windows by adding a sufficiently short window function. This ensures that the instantaneous frequency change of the signal within the short time window is approximately linear, thus satisfying the unbiased estimation condition of TFPF.

[0069] S3, by scaling the continuous phase sequence after endpoint expansion, the normalized signal is obtained:

[0070] To avoid frequency aliasing during frequency modulation, the amplitude of the continuous phase sequence after endpoint expansion is normalized:

[0071] ,

[0072] in, represent The minimum value in, represent The maximum value in the range, where a=0.5, b=0 is commonly used. The result is as follows: Figure 2 As shown in (d), the normalized result of the scaled signal amplitude between 0.5 and 0 is displayed.

[0073] S4, Frequency Modulation Code:

[0074] The normalized signal is frequency modulated, and the frequency modulation is encoded into the instantaneous frequency of the unit amplitude analytical signal.

[0075] Construct an analytic signal This makes it possible to recover the original signal. The normalized signal becomes the instantaneous frequency (IF) of this analytical signal. It is an integral variable. The representative will arrive Taken together.

[0076] ,

[0077] Calculating the instantaneous phase is .

[0078] The instantaneous frequency is the derivative of the phase of the analytic signal:

[0079] ,

[0080] Analyzing signals instantaneous frequency That is, the normalized signal that needs to be encoded. . It is the instantaneous phase.

[0081] S5, PWVD time-frequency analysis, performs PWVD transform on the analytic signal. PWVD transform exhibits time-frequency clustering because signals and noise behave differently in the time-frequency plane. Signal energy is concentrated along a time-varying frequency trajectory (i.e., instantaneous frequency) in the time-frequency plane. Random noise energy is diffusely distributed in the time-frequency plane. The time-frequency peak map can be obtained through PWVD transform, allowing for signal extraction.

[0082] ,

[0083] in: The result is after PWVD transformation, where w(l) is the window function, k is the frequency sampling point, and the symbol * indicates conjugate. For a sequence length of The analytical signal, Indicates that the sequence length is The conjugate of the analytic signal. In the time-frequency plane, the signal energy is concentrated near the instantaneous frequency, while the noise distribution is dispersed. The results are as follows: Figure 2 As shown in (e), random noise is dispersed across the time-frequency plane due to its random nature, resulting in proper separation between signal and noise components.

[0084] S6, Peak Extraction and Instantaneous Frequency Estimation:

[0085] By finding the maximum value of the frequency axis at each time point:

[0086] ,

[0087] And convert it to instantaneous frequency to obtain the IF estimated signal:

[0088] ,

[0089] in, This represents the number of frequency sampling points. The sampling rate is [value]. The results are as follows: Figure 2 As shown in (f), the signal energy is concentrated near the instantaneous frequency (IF), while random noise appears as a dispersed component. By extracting the peak value in the time-frequency domain, an unbiased estimate of the signal can be obtained, thereby effectively suppressing random noise.

[0090] S7, Inverse Scaling and Endpoint Truncation

[0091] Inverse scaling of the instantaneous frequency estimation signal:

[0092] ,

[0093] a and b are values ​​set by the user in the previous text; typically, a=0.5 and b=0 are used. This is the phase signal after inverse scaling.

[0094] And remove extended endpoints:

[0095] ,

[0096] in, For length is The inversely scaled phase signal, For length is The phase signal after removing the extended endpoints. The result is as follows: Figure 2 As shown in (g), the phase information of the interference fringes corresponding to the atomic gravity measurement was obtained.

[0097] S8, recover the interference fringes. The phase signal obtained from S7 is used to calculate the interference fringe signal after time-frequency peak filtering (TFPF) using cosine trigonometric functions. :

[0098] As mentioned above:

[0099] ;

[0100] ;

[0101] in, It can be done Calculations show that The effective wave vector is a fixed quantity. It is the Raman pulse time interval, a parameter used in instrument experiments. The sweep rate is also a known quantity during instrument experiments. (Center value) and contrast Both can be obtained by fitting the interference fringes. Therefore, the gravitational acceleration g can be obtained, as shown in the figure. Figure 2 As shown in (h), the signal-to-noise ratio of the interference fringes is significantly improved.

[0102] The dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering (TFPF) proposed in this application has the following significant advantages compared with existing least squares fitting methods, extended Kalman filtering (EKF), and external accelerometer methods:

[0103] No additional hardware required: Noise suppression is achieved solely through algorithmic processing, without relying on additional accelerometers or integrated navigation systems. This ensures system simplicity and portability, while reducing device size, weight, and power consumption.

[0104] Enhanced adaptability to dynamic environments: Vibration and noise are unavoidable in aircraft, shipborne, and vehicle-mounted platforms. This application achieves phase recovery of single-experiment data by directly acting on the interference fringe signal, avoiding reliance on repeated experiments. It maintains high accuracy even under low signal-to-noise ratio conditions.

[0105] Improved real-time performance: Time-frequency peak filtering outperforms extended Kalman filtering in terms of algorithmic complexity, requires no large amount of prior statistical data, and is suitable for embedded real-time processing. For continuously acquired interference fringe signals, this method can achieve rapid processing and output.

[0106] Significant noise suppression effect: This application utilizes the concentrated energy characteristics of instantaneous frequency to effectively suppress random noise. Compared with traditional filtering methods (low-pass filtering, wavelet transform), which introduce phase distortion, this method maintains the accuracy of gravity measurement.

[0107] High compatibility: This application can be used in parallel with existing data processing methods (least square fitting, Kalman filtering) as an additional noise suppression technique. It is implemented at the software level, making it highly adaptable.

[0108] In practical applications, a research vessel or a vessel with a displacement of 8 tons or more is selected. The core instruments are a cold atom gravimeter and a dual-axis inertial stabilization platform. A high-precision accelerometer is installed beneath the platform for vibration compensation, and a differential GPS antenna is mounted on the top of the hull to acquire position, velocity, and attitude information. In addition, a CG-5 high-precision relative gravimeter is provided as a benchmark calibration device for comparison and verification. The temperature inside the measurement chamber needs to be maintained at 20±2℃, and the humidity at approximately 55%±10%. The sailing speed is controlled at ≤20km / h, the instantaneous maximum vibration acceleration is less than 0.15m / s², and the sampling rate is 1Hz.

[0109] First, approximately 24 hours of static continuous measurements were conducted in the laboratory to evaluate the sensitivity and long-term stability of the atomic gravimeter. Then, the instrument was mounted on the ship, the inertial platform was activated, and attitude calibration was performed to ensure the atomic gravimeter remained vertically aligned. While the ship sailed along the planned route at approximately 10 km / h, atomic interference fringes, accelerometer vibration signals, GPS position, and attitude data were recorded in real time. This application employs a phase calculation and unwinding method. The phase is unwound to eliminate π jumps, and the continuous phase sequence is then extended at its endpoints to suppress PWVD (Shift-Window Short-Time Fourier Transform) boundary effects. Subsequently, the extended phase sequence is scaled to obtain a normalized signal. This normalized signal is then frequency-modulated and encoded into the instantaneous frequency of the unit amplitude analytic signal, yielding the analytic signal. Finally, the analytic signal undergoes PWVD transformation to obtain a time-frequency peak map. The maximum frequency axis value at each moment is found from the time-frequency peak map and converted into an instantaneous frequency, resulting in an estimated instantaneous frequency signal. This estimated instantaneous frequency signal is scaled, and the extended endpoints are removed. Interference fringe signals are calculated using cosine trigonometric functions, and gravity is calculated based on these signals. The absolute gravity value after correction at the CG-5 reference point is compared to calculate the internal and external coincidence accuracy. The entire experiment involved four repeated survey lines, each requiring one round trip to ensure uniform straight-line navigation.

[0110] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A method for suppressing dynamic atomic gravity interference fringe noise based on time-frequency peak filtering, characterized in that, The method includes: Phase calculation and unwinding, the phase calculation and unwinding includes: performing a preliminary phase estimate on the measurement data of atomic interference fringes, and performing phase unwinding on the preliminary phase estimate to obtain a continuous phase sequence; Endpoint extension of a continuous phase sequence; The normalized signal is obtained by scaling the continuous phase sequence after endpoint expansion. The normalized signal is frequency modulated, and the frequency modulation is encoded into the instantaneous frequency of the unit amplitude analytic signal to obtain the analytic signal. Perform PWVD transformation on the analytical signal to obtain the time-frequency peak diagram; Find the maximum value of the frequency axis at each moment from the time-frequency peak plot, and convert it into the instantaneous frequency to obtain the instantaneous frequency estimation signal; The instantaneous frequency estimation signal is scaled and the extended endpoints are removed. The interference fringe signal is then calculated using the cosine trigonometric function. Gravity is calculated based on the interference fringe signal.

2. The dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering according to claim 1, characterized in that, A preliminary phase estimate is made from the measurement data of atomic interference fringes, including obtaining the preliminary phase estimate using inverse trigonometric functions, expressed by the formula: ， in, For preliminary phase estimation, A represents the actual value of the atomic transition probability, A is the center value, and B is the contrast ratio.

3. The dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering according to claim 1, characterized in that, The initial phase estimate is de-wound to obtain a continuous phase sequence, including: calculating arbitrary phases using a phase de-wound algorithm. Phase value at time : ， Where Z represents the set of integers, including positive integers, zero, and negative integers. The value of makes the phase difference between adjacent points satisfy: ， for Phase of time, for Phase of time, This is a preliminary estimate of the phase.

4. The dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering according to claim 1, characterized in that, Endpoint extension of a continuous phase sequence is achieved using the following formula: ， Where the length is M, x (n) is the extended signal with a length of 2p+M, x(1) represents the first number of the unextended signal, x(M) represents the Mth number of the unextended signal, x(np) represents the npth number of the original signal, M is the length of the phase sequence of the unextended signal, p is the number of points to be extended, p≥L+1, and L is the pseudo-Wegener distribution window.

5. The dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering according to claim 1, characterized in that, By scaling the continuous phase sequence after endpoint expansion, a normalized signal is obtained, including: ， in, The signal after normalization represent The minimum value in, represent The maximum value in, This represents the expanded signal. b are empirical values.

6. The dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering according to claim 1, characterized in that, The analytical signal is: ， in, It is an integral variable. The representative will arrive Summation, Indicates an analytical signal. The integral variable is The normalized signal.

7. The dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering according to claim 1, characterized in that, The formula used to perform PWVD transformation on the analytic signal is: ， in: The result is after PWVD transformation, where w(l) is the window function, k is the frequency sampling point, and the symbol * indicates conjugate. For a sequence length of The analytical signal, Indicates that the sequence length is The conjugate of the analytic signal.

8. The dynamic atomic gravity interference fringe noise suppression method based on time-frequency peak filtering according to claim 7, characterized in that, Find the maximum frequency value at each moment from the time-frequency peak plot and convert it to the instantaneous frequency to obtain the instantaneous frequency estimation signal, including: Calculate the maximum value of the frequency axis: , Indicates the maximum value of the frequency axis; Convert the maximum value of the frequency axis to the instantaneous frequency to obtain the instantaneous frequency estimation signal: ， in, This represents the number of frequency sampling points. The sampling rate.

Images