Polar shift period term extraction and reconstruction optimization method and system

Abstract

The invention discloses a polar shift periodic term extraction and reconstruction optimization method, which comprises the following steps: S1, acquiring a polar shift observation sequence, and constructing an original complex polar shift sequence based on the polar shift observation sequence; s2, performing complex singular spectrum analysis on the original complex polar shift sequence to obtain a Chandler term component and an annual term component; s3, respectively carrying out complex Fourier band-pass filtering processing on the original complex polar shift sequence, the Chandler term component and the annual term component; and S4, carrying out complex least square gain matching on the filtering results of the original complex number polar shift sequence, the Chander term component and the annual term component and the Chander term component and the annual term component obtained by complex number singular spectrum analysis to obtain a final polar shift period term. On the basis, robust extraction and continuous reconstruction of a polar shift period term in a polar shift observation sequence are realized through complex singular spectrum analysis and complex Fourier band-pass filtering, and end point distortion and a frequency offset effect are effectively weakened.

Description

Technical Field

[0001] This invention relates to the field of technology, specifically to a method and system for extracting and reconstructing polar shift periodic terms. Background Technology

[0002] Quasi-periodic signals are commonly found in Earth Rotation Parameters (ERP) time series, particularly the Chandler oscillation and annual term of the Polar Motion (PM) component. Accurate extraction and reconstruction of these periodic terms are crucial for building high-precision forecasting models and are also significant for near-real-time applications such as precise orbit determination, precise positioning, and numerical weather prediction.

[0003] Common methods for extracting quasi-periodic terms include least squares fitting, Fourier Transform Band-Pass Filtering (FTBPF), and Singular Spectrum Analysis (SSA), which are widely used in geodesy and geophysical time series processing. SSA, derived from the Karhunen-Loève decomposition concept, is a non-parametric method for principal component analysis of one-dimensional time series. It does not require strict prior frequency assumptions and is not constrained by pure sine assumptions. It can simultaneously extract trend and time-varying periodic terms from noisy sequences, making it particularly suitable for processing periodic information where frequencies drift slowly over time. FTBPF, as a classic frequency domain method, performs global analysis of the entire time series, cleanly separating the signal in the target frequency band. It boasts fast computation speed and high frequency localization accuracy, and is also widely used in extracting quasi-periodic terms from geoscientific signals.

[0004] However, the above methods still have limitations in time series applications: 1) SSA's embedding dimension (lag window) dependence - too large a dimension can easily lead to mode "aliasing", while too small a dimension makes it difficult to achieve gradual separation from weak to strong, resulting in energy leakage and component mixing between periodic terms and noise / trend; 2) FTBPF's boundary fragility - when there is noise or the trend is not fully subtracted, the filter output of the boundary segment is highly sensitive to the endpoint conditions, and is prone to amplitude and phase deviation and endpoint distortion; while boundary information is the key link in constructing continuous reconstruction and carrying out forecasting. Summary of the Invention

[0005] To overcome the limitations of traditional Singular Spectrum Analysis (SSA) and Fourier Transform Band-Pass Filtering (FTBPF) methods in terms of embedding dimension sensitivity, mode mixing, and boundary filtering distortion, this invention provides an optimized method and system for polar motion periodic term extraction and reconstruction. Based on joint analysis of SSA and FTBPF in the complex domain, it achieves robust extraction and continuous reconstruction of polar motion periodic terms in polar motion observation sequences, effectively reducing endpoint distortion and frequency offset effects.

[0006] According to one aspect of the present invention, a method for extracting and reconstructing polar shift periodic terms is provided, comprising: step S1, acquiring a polar shift observation sequence and constructing an original complex polar shift sequence based on the polar shift observation sequence; step S2, performing complex singular spectrum analysis on the original complex polar shift sequence to obtain Chandler term components and annual term components; step S3, performing complex Fourier bandpass filtering on the original complex polar shift sequence, Chandler term components, and annual term components respectively; step S4, performing complex least squares gain matching on the filtering results of the original complex polar shift sequence, Chandler term components, and annual term components, so that the filtering results are consistent with the structure of the Chandler term components and annual term components obtained by complex singular spectrum analysis in step S2 on the complex plane, thereby obtaining the final polar shift periodic term.

[0007] Further, step S1 includes: step S11, obtaining polar motion observation sequences under the same reference frame and unified epoch; step S12, extracting dual-channel polar motion components from the polar motion observation sequences; step S13, constructing an original complex polar motion sequence based on the dual-channel polar motion components.

[0008] Further, step S2 includes: step S21, constructing a polar-shift complex trajectory matrix based on the original complex polar-shift sequence; step S22, performing singular value decomposition on the polar-shift complex trajectory matrix to obtain multiple singular values; step S23, performing a diagonal average operation on the matrix corresponding to each singular value to obtain reconstructed components; step S24, performing a fast Fourier transform on the reconstructed components to obtain the spectrum corresponding to each element in the reconstructed components; step S25, performing spectral peak detection on the spectrum corresponding to each element in the reconstructed components to obtain the peak frequency and calculate the spectral energy; step S26, selecting paired reconstructed components using the peak frequency of each element in the reconstructed components; step S27, reconstructing the annual term component and Chandler term component based on the selected reconstructed components; step S28, performing positive spectrum estimation on the reconstructed Chandler term component and annual term component, searching for a set of spectral peaks within the reference frequency window, selecting the main peak, and constructing a passband for each main peak.

[0009] Further, in step S26, if the relative deviation of the peak frequencies of any two elements in the reconstructed components does not exceed a preset relative tolerance, then the two elements are considered to constitute a pair of reconstructed components. The cumulative spectral energy value of the pair of reconstructed components is calculated, and all the pair of reconstructed components are sorted according to the magnitude of the cumulative spectral energy value. Several pair of reconstructed components with the highest cumulative spectral energy value are selected.

[0010] Further, step S3 includes: step S31, constructing a complex response function for the Chandler term using the frequency band range of the Chandler term, and constructing a complex response function for the annual term using the frequency band range of the annual term; step S32, performing complex Fourier bandpass filtering on the original complex polar shift sequence using the complex response function for the Chandler term and the complex response function for the annual term, respectively, to obtain the Chandler term sequence and the annual term sequence of the filtered original complex polar shift sequence; step S33, performing complex Fourier bandpass filtering on the Chandler term component and the annual term component obtained after complex singular spectrum analysis using the complex response function for the Chandler term and the complex response function for the annual term, respectively, to obtain the Chandler term sequence of the filtered Chandler term component and the annual term sequence of the annual term component.

[0011] Furthermore, in step S31, both the Chandler term complex response function and the anniversary term complex response function are constructed using a complex response function with a soft-edge transition design. The formula for the complex response function is as follows: , in, Represents the complex response function; This indicates the calculation of cosine. Indicates frequency; Indicates the starting point of the left soft edge; Indicates the start of the left passband; Indicates the end of the right passband; This indicates the end point of the right soft edge.

[0012] Further, in step S32, the original complex polar shift sequence is subjected to complex Fourier bandpass filtering, and the corresponding formula is: , In the formula, Represents the primitive complex polar shift sequence. and These represent the annual term sequence and the Chandler term sequence obtained by directly performing complex Fourier bandpass filtering on the original complex polar shift sequence, respectively. Indicates the inverse Fourier transform; This represents the complex response function of the Chandler term constructed using the frequency band range of the Chandler term. This represents the complex response function of the annual term constructed using the frequency band range of the annual term.

[0013] Furthermore, the annual term component and Chandler term component extracted from complex singular spectrum analysis are subjected to complex Fourier bandpass filtering, and the corresponding formulas are as follows: , In the formula, and These represent the annual term sequence and the Chandler term sequence obtained by performing complex singular spectrum analysis followed by complex Fourier bandpass filtering, respectively. and This represents the annual term component and the Chandler term component obtained from the complex singular spectrum analysis in step S2.

[0014] Further, in step S4, complex least squares gain matching is performed, and the corresponding formula is: , in, This represents the final polar shift periodic term obtained by complex least squares gain matching; This indicates the objects to be subjected to gain matching, including the original complex polar shift sequence, the Chandler term component, and the filtering results of the annual term component; Indicates time; Represents the complex gain coefficient. The calculation formula is: , in, This represents the annual and Chandler term components obtained from complex singular spectrum analysis; This indicates the filtering result obtained by directly applying a complex Fourier bandpass filter to the original complex polar shift sequence; * indicates complex conjugate. This indicates a modulo operation on a signal in a complex sequence.

[0015] According to one aspect of the present invention, a polar shift periodic term extraction and reconstruction optimization system is provided, comprising: a polar shift dual-channel construction and complex domain representation module for acquiring a polar shift observation sequence and constructing an original complex polar shift sequence based on the polar shift observation sequence; a complex singular spectrum analysis module for performing complex singular spectrum analysis on the original complex polar shift sequence to obtain Chandler term components and annual term components; a complex Fourier bandpass filtering module for performing complex Fourier bandpass filtering on the original complex polar shift sequence, Chandler term components, and annual term components respectively; and a complex least squares gain matching module for performing complex least squares gain matching on the filtered results of the original complex polar shift sequence, Chandler term components, and annual term components with the Chandler term components and annual term components obtained from the complex singular spectrum analysis to obtain the final polar shift periodic term.

[0016] The above technical solution extends SSA and FTBPF to Complex Singular Spectrum Analysis (CSSA) and Complex Fourier Transform Band-Pass Filtering (CFTBPF), respectively. Based on this, CSSA is performed on the original complex polar shift sequence to obtain the Chandler term component and the annual term component. CFTBPF is then applied to the original complex polar shift sequence, the Chandler term component, and the annual term component to cleanly separate the signal in the target frequency band. Finally, to eliminate the system offset caused by the filter on amplitude and phase, and to ensure that the filtering result is strictly consistent with the signal structure extracted by CSSA in the complex plane, complex least squares gain matching is performed on the filtering result. This achieves robust extraction and continuous reconstruction of the periodic term in the polar shift time series, effectively reducing endpoint distortion and frequency offset effects.

[0017] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0018] (1) Reducing the problem of embedding dimension sensitivity: Embedding dimension sensitivity refers to the requirement for the number of rows (i.e., dimension) of the trajectory matrix when using Singular Spectral Analysis (SSA). If the dimension is too small, it is impossible to distinguish between Chandler terms (~430 days) and annual terms (~365 days); if the dimension is too large, too many components without physical meaning will be generated, resulting in a heavy computational burden. The present invention dynamically determines the dimension based on actual data, thereby reducing the problem of embedding dimension sensitivity.

[0019] (2) Reducing Mode Aliasing: Mode aliasing refers to the mixing of signal components from different physical processes within one or more RC components. For example, in traditional SSA, mode aliasing includes components containing trend + part Chandler, components containing Chandler + part annual, and components containing annual + noise. The specific meaning of aliasing is: components with similar frequencies are difficult to separate; the physical meaning is ambiguous and difficult to interpret; and time-varying signals are incorrectly decomposed. This invention reduces mode aliasing by selecting paired reconstructed components (through energy filtering and pairwise discrimination).

[0020] (3) Reducing Boundary Filtering Distortion: Boundary filtering distortion refers to the spurious oscillations that occur at the boundaries of the processed data in traditional FTBPF. Specifically, this manifests as: spurious fluctuations at the beginning and end of the sequence; amplitude and phase distortion at the boundaries; and unreliable boundary information affecting extrapolation prediction. This invention reduces boundary filtering distortion through soft-edge transition design (which significantly reduces boundary oscillations), complex gain matching (maintaining the proportional relationship between the filtering result and the physical signal; correcting the phase delay introduced by filtering; and improving signal quality at the boundaries). Attached Figure Description

[0021] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the accompanying drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0022] Figure 1 The flowchart illustrates a method for extracting and reconstructing optimization of polar shift periodic terms, as provided in an embodiment of the present invention.

[0023] Figure 2 A comparison chart of Chandler swing time-domain extraction results provided in an embodiment of the present invention.

[0024] Figure 3 A comparison chart of the annual oscillation time-domain extraction results provided in the embodiments of the present invention.

[0025] Figure 4The image shows the time-domain extraction results of long-term trend items provided in an embodiment of the present invention.

[0026] Figure 5 A comparison diagram of the Chandler swing spectrum provided in an embodiment of the present invention.

[0027] Figure 6 A comparison chart of the annual oscillation spectrum provided for embodiments of the present invention.

[0028] Figure 7 A spectrum comparison chart of long-term trend terms provided for embodiments of the present invention. Detailed Implementation

[0029] It should be noted that:

[0030] The terms “comprising” and “having”, and any variations thereof, in the specification, claims, and accompanying drawings of this invention are intended to cover a non-exclusive inclusion, such as a process, method, system, product, or apparatus that includes a series of steps or units, not necessarily limited to those explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0031] The block diagrams shown in the accompanying drawings are merely functional entities and do not necessarily correspond to physically independent entities. That is, these functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor devices and / or microcontroller devices. The flowcharts shown in the accompanying drawings are merely illustrative and do not necessarily include all content and operations / steps, nor do they necessarily have to be performed in the described order. For example, some operations / steps can be decomposed, while others can be combined or partially combined; therefore, the actual execution order may change depending on the specific circumstances.

[0032] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. In addition, the technical features of the various embodiments or individual embodiments provided by the present invention can be arbitrarily combined to form new technical solutions. Such combinations are not bound by the order of steps and / or structural composition patterns, but must be based on the ability of those skilled in the art to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed by the present invention.

[0033] Please refer to the appendix. Figure 1 This invention provides a method for extracting and reconstructing optimization of polar shift periodic terms, specifically including the following steps:

[0034] Step S1: Obtain the polar motion observation sequence and construct the original complex polar motion sequence based on the polar motion observation sequence.

[0035] It should be noted that step S1 implements "polar shift dual-channel construction and complex field representation", and step S1 will be further described below.

[0036] Step S11: Obtain the polar motion observation sequence under the same reference frame and unified epoch.

[0037] It should be noted that polar motion is a type of Earth orientation parameter. Polar motion consists of X-axis and Y-axis components. Polar motion is observed once a day, meaning one X-axis component and one Y-axis component are observed each day, resulting in 30 observations per month. All polar motion observations form a series of XY polar motion components that grows over time; this series of data is called the polar motion observation sequence (i.e., the polar motion time series). This invention processes polar motion data directly in complex space, preserving the phase correlation between the two components.

[0038] The reference frame used in this embodiment is ITRF2020, and the epoch is 2010. The calculation process uses polar motion observation sequences extracted from Earth orientation parameter files from 1985 to May 2025, located within this reference frame and epoch. Understandably, the reference frame, epoch, and the specific time period for the polar motion observation sequences can all be adjusted according to actual needs, and are not limited here.

[0039] Step S12: Extract the polar motion components of the two channels from the polar motion observation sequence.

[0040] It should be noted that the polar shift component of the dual-channel refers to the polar shift... The time series of the polar shift X component (i.e., the data series in which the polar shift X component changes over time) and polar shift The time series of the polar shift Y component (i.e., the data series of the polar shift Y component changing with time), polar shift Time series and polar shift of components The time series of the components together constitute the polar motion observation sequence. This invention uses... Indicated by time polar shift of variables Time series of components; Indicated by time polar shift of variables The time series of the components. Simultaneously record... for ;remember for Then the polar shift will be... Components and polar shift The components are considered as the measured dual channels, that is:

[0041] (1.1)

[0042] Step S13: Construct the original complex polar shift sequence using the polar shift components of the two channels, as follows:

[0043] (1.2)

[0044] in, Represents complex components; For time The original complex polar shift sequence is constructed with the independent variable as the variable.

[0045] Step S2: Perform complex singular spectrum analysis on the original complex polar shift sequence to obtain the Chandler term component and the annual term component.

[0046] It should be noted that step S2 performs Complex Singular Spectrum Analysis (CSSA). Please refer to the appendix for details. Figures 2 to 4 , Figure 2 The diagram illustrates the complex signal of the Chandler term extracted using CSSA, CFTBPF, and CSSA+CFTBPF, decomposed back into the time series plots of the two polar shift components X and Y. Figure 3 The diagram illustrates the time series plots of the complex signals of the annual term extracted using CSSA, CFTBPF, and CSSA+CFTBPF, decomposed back into the two polar shift components X and Y. Figure 4 This illustrates the long-term trend signal, where the long-term trend for all three methods is set to the trend extracted using CSSA. For example... Figure 2As shown, CSSA (black dashed line) represents the signal directly reconstructed through complex singular spectrum analysis. It appears the "cleanest," with smooth fluctuations, because it is essentially composed of a few main reconstructed RC components (such as a pair of conjugate components), filtering out other noise. CFTBPF (red solid line) represents the signal obtained by filtering directly from the raw data using an adaptive soft-sideband pass filter. It contains all the energy within the frequency band, thus the curve is richer in detail and may contain some small signals or noise that were not selected by CSSA. CSSA+CFTBPF (blue solid line): This is a hybrid method. First, CSSA is used to extract the main components, and then CFTBPF filtering is applied to this "cleaned" signal. Its shape is usually between the two mentioned above, containing more information than pure CSSA and smoother than pure CFTBPF. Step S2 will be further described below:

[0047] Step S21: Construct the complex trajectory matrix of polar shifts based on the original complex polar shift sequence, as follows:

[0048] (1.3)

[0049] in, This represents the original complex polar shift sequence at a specific moment in the total time series. Represents the complex trajectory matrix of polar shifts. C represents the polar motion complex trajectory matrix. In the complex field; Indicates the number of sliding windows, and , This represents the total length of the time series (number of sample points). Indicates the window length. The calculation method for is as follows:

[0050] (1.4)

[0051] in, This indicates rounding down to the nearest whole number. This indicates the CSSA window length. Understandably, this invention dynamically determines the dimension L based on actual data, with a minimum of 30 days and a maximum not exceeding the data length - 2 (i.e., the result calculated using Equation 1.4). This yields a more suitable dimension, thereby reducing the problem of embedding dimension sensitivity.

[0052] This embodiment selects (The least common multiple of the corresponding polar shift period) to balance frequency resolution and time smoothness. , The average sampling interval (days) is represented by the adjacent Modified Julian Day. )difference( (This indicates a difference operation) Take the median value ( The result is obtained by taking the median value.

[0053] Step S22, for the polar motion complex trajectory matrix Perform singular value decomposition to obtain multiple singular values, as follows:

[0054] (1.5)

[0055] in, This represents the left singular vector matrix (time-domain eigenvectors that describe the oscillation mode of the signal within a time window, with dimensions L×R). This represents a singular value diagonal matrix, where singular values ​​lie on the diagonal of this matrix (energy intensity matrix, dimension R×R). This represents the right singular vector matrix (phase space directional characteristic matrix, describing the evolution trajectory of the signal in phase space, with dimensions K×R); R represents the matrix transpose operation; R represents the number of singular values.

[0056] It should be noted that, firstly, singular spectrum analysis is performed on the polar motion complex trajectory matrix, and then the first R singular values ​​are selected. The value of R (i.e., the value of R) is actually a value that varies depending on the data used. The method for determining R is as follows: firstly, R is manually defined. max (i.e., the maximum value of R), for example, R max =20, meaning R will not exceed 20. The value of R will be adjusted based on L and K mentioned above, specifically by taking all three values ​​(R...). max The minimum value of (L, K) is the value of R.

[0057] Step S23: Perform a diagonal average operation on the matrix corresponding to the selected singular values ​​to obtain the reconstructed components. .

[0058] Specifically, after singular value decomposition, the first R singular values ​​are selected, and then the matrices corresponding to the selected singular values ​​are processed. Perform a diagonal averaging operation to obtain the corresponding reconstructed component RC. Here, the reconstructed component in complex space... It can obviously also be written in the form of complex polar shift components:

[0059] (1.6)

[0060] in, Represents reconstructed components The complex polar shift components, This represents the polar shift x-direction component in the reconstructed RC time series. This represents the polar shift y-direction component in the reconstructed RC time series, where j is the imaginary unit.

[0061] Step S24: Perform a Fast Fourier Transform on the reconstructed components to obtain the spectrum of each element in the reconstructed components, as follows:

[0062] (1.7)

[0063] in, This represents the signed spectrum obtained after the Fast Fourier Transform. The complex polar shift component represents the reconstructed component. This represents the Fast Fourier Transform.

[0064] Step S25: Perform spectral peak detection on the spectrum of each element in the reconstructed component to obtain the peak frequency, and calculate the corresponding spectral energy.

[0065] In step S25, for each element in the reconstructed component Calculate the signed spectrum, take the positive frequency portion, and calculate the frequency falling within the initial mutually exclusive periodic window. , This indicates the range of periodic values ​​(e.g., the initial mutex window for the Chandler term is...). The initial mutual exclusion window for the anniversary item is Spectral energy within) and the peak frequency within this frequency band ( (Represents the peak element).

[0066] 1) Regarding the spectrum Perform spectrum main peak detection as follows:

[0067] (1.8)

[0068] in, Represents the complex polar shift components Amplitude modulus operation of the signed spectrum below; This represents the peak frequency within the initial frequency band (i.e., the frequency band corresponding to the initial mutually exclusive periodic window), that is, the frequency under the condition of maximum amplitude within the initial frequency band. The corresponding frequency; This represents the minimum frequency value within the initial frequency range. This indicates the maximum frequency value within the initial frequency range.

[0069] 2) In-band spectral energy The calculation formula is as follows:

[0070] (1.9)

[0071] in, express Energy within the target frequency band; Indicates integration operation; Indicates the maximum frequency within the target frequency band; Indicates the minimum frequency within the target frequency band; express Amplitude at a given frequency; This represents the modulo operation of the spectral amplitude of a complex sequence; This represents the time element.

[0072] Step S26: Selecting pairs of reconstructed components using the peak frequency and spectral energy of each element in the reconstructed components.

[0073] In step S26, if the reconstructed components Any two elements and RC j peak frequency satisfy:

[0074] (1.10)

[0075] In the above formula, and They represent and Peak frequency below; express Determine the relative tolerances of the paired reconstructed components (e.g., The value can be 3%, but no limit is specified here. Indicates in and Select the maximum value from the list.

[0076] From the above equation, we know that if the reconstructed components... Any two elements and RC j peak frequency and The relative deviation does not exceed the preset relative tolerance. Then the corresponding and RC j This forms a pair of reconstructed components. All paired reconstructed components are then configured according to... Sort the values ​​of the accumulated spectral energy (i.e., the cumulative spectral energy values) by size, and select the values ​​with the highest accumulated spectral energy values. Each pair of reconstructed components forms a Top-N pair. In this embodiment, two Top-N pairs are selected for the Chandler term, and one Top-N pair is selected for the anniversary term. Understandably, the multi-peak adaptive passband used in this invention selects two Top-N pairs for the Chandler term in order to adapt to the multi-peak characteristics of the Chandler oscillation when constructing the initial target passband.

[0077] Step S27: Reconstruct the anniversary item component and the Chandler item component based on the selected paired reconstruction components, and you will get:

[0078] (1.11)

[0079] in, This indicates the time-based index extracted by CSSA. The Chandler component of the independent variable; This indicates the time-based index extracted by CSSA. For the annual term component of the independent variable.

[0080] Step S28: Perform positive spectrum estimation on the reconstructed Chandler term components and annual term components, find the set of spectral peaks within the reference frequency window, select the main peaks, and construct a passband for each main peak.

[0081] In step S28, within the reference frequency window ( and Calculated from the initial mutually exclusive periodic window given by the Chandler term or the anniversary term. and The 0 in the formula represents the initial frequency band range. The search is conducted within this range to find the set of spectral peaks, and the peak with the highest energy within that range is selected. Each peak (understandably, the peak with the highest energy is the main peak):

[0082] (1.12)

[0083] Furthermore, a relative bandwidth of [value] is constructed for each main peak. The passband is as follows:

[0084] (1.13)

[0085] In the above formula and These represent the lowest and highest frequencies of the newly constructed passband; Indicates the maximum period within the initial target passband; This represents the minimum period within the initial target passband; Indicates relative bandwidth; Indicates selection The largest of them; Indicates selection The smallest of the three. Multiple passbands are allowed to cover possible multi-peak / split spectral structures, such as the common multi-peak Chandler term. It should be noted that the initial mutually exclusive periodic window refers to the periodic interval where the initially given Chandler and annual terms appear, while the initial target passband refers to one or more frequency ranges separated after CSSA processing based on the initial mutually exclusive periodic window, because the Chandler term allows for multiple frequency ranges.

[0086] Step S3: Perform complex Fourier bandpass filtering on the original complex polar shift sequence, the Chandler term component, and the anniversary term component to reconstruct the Chandler term and the anniversary term.

[0087] In step S3, to verify the consistency between the principal periodic components extracted by CSSA (complex singular spectral analysis) and the theoretical bandpass results, and to remove residual noise from CSSA, the original complex polar shift sequence is then processed. and the Chandler item component extracted using CSSA With anniversary items Frequency band ranges using the Chandler and Anniversary terms respectively. Complex Fourier Transform Band-Pass Filtering (CFTBPF) is performed. Simultaneously, the positive frequency component corresponds to forward motion, and the negative frequency component corresponds to backward motion.

[0088] Please refer to the appendix for details. Figures 5 to 7 , Figure 5 This illustrates the comparison of Chandler's oscillation spectrum. Figure 6 This illustrates the comparison of the annual oscillation spectrum. Figure 7 This illustrates the comparison of the spectrum of long-term trend terms. For example... Figure 5 As shown, the RAW curve represents the spectrum of the original polar shift signal. CSSA is the spectrum of the Chandler term after CSSA extraction; CFTBPF is the spectrum obtained by directly bandpass filtering the original data. It exhibits excellent out-of-band suppression, demonstrating the filter's ability to "purify" the signal. CSSA+CFTBPF is the spectrum obtained by first extracting with CSSA and then filtering with CFTBPF. It combines the sharp peaks of CSSA with the out-of-band purity of CFTBPF. Step S3 will be further described below:

[0089] Step S31: Construct the complex response function of the Chandler term using the bandwidth of the Chandler term, and construct the complex response function of the anniversary term using the bandwidth of the anniversary term.

[0090] In step S31, it should be noted that the passband frequency ranges of the anniversary term and the Chandler term are given by equation (1.13). Here, the forward and reverse signals are processed separately. When processing the forward signal (the same steps are performed for the reverse signal), the complex response function of the Chandler term... Complex response function of the anniversary term Both employ the complex response function of a bandpass filter in the frequency domain. Constructed from complex response functions. The soft-edge transition design is adopted, and the formula is:

[0091] (1.14)

[0092] Represents the complex response function; This indicates the calculation of cosine. Indicates frequency; This indicates the starting frequency of the left soft edge (where a soft edge means that the transition from the stopband to the passband and from the passband to the stopband is not directly cut off at the junction of the target passband and the other frequency bands (i.e., the stopband), but rather a more "gentle" transition is made using a linear "slope". This is the starting frequency from the stopband to the passband. This indicates the start of the left passband, i.e., the frequency at which the maximum gain is reached from the start of the passband; This indicates the end of the right passband, i.e., the frequency at which the decay begins in the passband; This indicates the end of the right soft side, i.e., the frequency at which the passband drops to the stopband; for those in the middle... For stopbands outside the frequency range, the complex response function is set to 0, thus suppressing such signals. This soft-edge transition method suppresses edge leakage while allowing multi-passband splicing.

[0093] Step S32: The original complex polar shift sequence is subjected to complex Fourier bandpass filtering using the Chandler term complex response function and the annual term complex response function, respectively, to obtain the Chandler term sequence and the annual term sequence of the filtered original complex polar shift sequence.

[0094] In step S32, the original complex polar shift sequence is processed. Frequency band ranges using the Chandler and Anniversary terms respectively. Perform their respective complex response functions and The filtered signal can be obtained by performing an inverse Fourier transform:

[0095] (1.15)

[0096] In the above formula and These represent directly using the original complex polar shift sequence. The sequence of anniversary items and the sequence of Chandler items obtained by performing CFTBPF; Indicates the inverse Fourier transform; Indicates the frequency band range using the Chandler item. The constructed Chandler term complex response function, Indicates the frequency band range using the anniversary term. Construct a complex response function for the annual term.

[0097] Step S33: The Chandler term complex response function and the anniversary term complex response function are used to perform complex Fourier bandpass filtering on the Chandler term components and anniversary term components obtained after CSSA processing, respectively, to obtain the Chandler term sequence of the filtered Chandler term components and the anniversary term sequence of the anniversary term components.

[0098] In step S33, the frequency band ranges of the anniversary term component and the Chandler term component extracted from CSSA are used respectively. Perform their respective complex response functions and The filtered signal can be obtained by performing an inverse Fourier transform:

[0099] (1.16)

[0100] In the above formula and These represent the sequence of anniversary items and the sequence of Chandler items obtained by first performing CSSA processing and then CFTBPF processing, respectively. and This represents the anniversary item component and the Chandler item component obtained from the CSSA processing in step S2.

[0101] Step S4: Perform complex least squares gain matching between the filtered results of the original complex polar shift sequence, Chandler term component, and annual term component and the Chandler term component and annual term component obtained from complex singular spectrum analysis to obtain the final polar shift periodic term.

[0102] In step S4, after completing the CFTBPF (complex Fourier bandpass filter) process, in order to eliminate the system offset caused by the filter on amplitude and phase, and to make the filtering result strictly consistent with the signal structure extracted by CSSA (complex singular spectrum analysis) in the complex plane, complex least squares gain matching will be performed on the filtering result.

[0103] The objective function for complex least squares gain matching is:

[0104] (1.17)

[0105] in, This represents the minimization operation; This represents the anniversary and Chandler components obtained using CSSA; This indicates the object to be matched for gain (e.g., , , , ); Indicates time; The complex gain coefficient is calculated as follows:

[0106] (1.18)

[0107] in,; This represents the filtering result obtained by directly applying CFTBPF to the original complex polar shift sequence; * indicates complex conjugate. This indicates that the modulus operation is performed on the signal under the original complex polar shift sequence, and the result is calculated from the above formula:

[0108] (1.19)

[0109] in, This represents the signal sequence after complex least-squares gain matching, i.e., the final polar shift periodic term (e.g., , , , (Signal sequence after complex least squares gain matching) Indicates input; This indicates the output.

[0110] In some feasible embodiments, step S5 is also included: calculating metrics and measuring performance. Referring to Tables 1 to 4, this invention evaluates both in-band and out-of-band metrics together, resulting in more diverse and reliable evaluation metrics.

[0111] Table 1 Chandler Full-Domain Aperture Specifications

[0112]

[0113] As shown in Table 1, in the full-frequency domain metrics, after CSSA (Complex Singular Spectrum Analysis), the signal-to-noise ratio (SNR) is only about 2.75 dB, and the RMSE is reduced by 9.5% compared to the original polar shift time series. This indicates that although it can separate the main Chandler components, it still contains a significant amount of out-of-band energy or noise. Meanwhile, CFTBPF (Complex Fourier Bandpass Filter) significantly improves the SNR to 14.20 dB, indicating high purity in frequency band extraction; however, its full-frequency domain RMSE improvement is relatively small, only about 6%. When using CSSA+CFBPF, the SNR further increases to 14.49 dB, exceeding the level when the filter is used alone, but the RMSE improvement is similar to that of CFTBPF, meaning it can maintain the phase consistency of CSSA while improving the purity of the extracted signal.

[0114] Table 2 Annual Full-Frequency Aperture Indicators

[0115]

[0116] As can be seen from Table 2, under the full frequency domain index, when using CSSA, CFTBPF, and CSSA+CFTBPF for the annual term extraction, the signal-to-noise ratio (SNR) is similar under CFTBPF and CSSA+CFTBPF, both showing a significant improvement (approximately 40%) compared to using CSSA alone. However, the reduction in RMSE under CFTBPF and CSSA+CFTBPF is slightly less than that under CSSA alone.

[0117] Table 3. Chandler In-Band Indicators

[0118]

[0119] As shown in Table 3, under the Chandler term's in-band metrics, CSSA maintains an in-band SNR of 2.75 dB, but its in-band RMSE decreases by 30.5%, indicating that it can still effectively capture the main energy within the target frequency band. CFTBPF's in-band RMSE decreases by 100%, meaning the filter extracts the target frequency band signal almost completely, with the smallest signal reconstruction error. Meanwhile, the combined method (CSSA+CFTBPF) has a slightly lower in-band RMSE decrease (approximately 90.9%), but the highest SNR and the lowest leakage ratio (3.43%, slightly better than CFTBPF's 3.66%).

[0120] Table 4 Annual In-Band Indicators

[0121]

[0122] As shown in Table 4, under the in-band performance metrics for the annual term, CSSA reduces the in-band RMSE by approximately 63%, indicating that it can effectively recover the main energy of the annual term within the target frequency band. CFTBPF achieves a 100% reduction in in-band RMSE, meaning the filter almost perfectly preserves the annual term frequency band signal. The in-band RMSE reduction of CSSA+CFTBPF is slightly lower (approximately 87.7%), but the SNR remains at a high level of 15.8 dB, with a band edge leakage of 2.56%, which is basically equivalent to CFTBPF's 2.42%. This indicates that the extraction purity of both methods is almost identical, while the joint method has better time-domain and phase consistency.

[0123] Based on the same technical concept as the aforementioned embodiments, this invention also provides a polar shift periodic term extraction and reconstruction optimization system, comprising: a polar shift dual-channel construction and complex domain representation module for acquiring polar shift observation sequences and constructing an original complex polar shift sequence based on the polar shift observation sequences; a complex singular spectrum analysis module for performing complex singular spectrum analysis on the original complex polar shift sequence to obtain Chandler term components and annual term components; a complex Fourier bandpass filtering module for performing complex Fourier bandpass filtering on the original complex polar shift sequence, Chandler term components, and annual term components respectively; and a complex least squares gain matching module for performing complex least squares gain matching on the filtered results of the original complex polar shift sequence, Chandler term components, and annual term components with the Chandler term components and annual term components obtained from the complex singular spectrum analysis to obtain the final polar shift periodic term.

[0124] In summary, this invention extends SSA and FTBPF to CSSA and CFTBPF, respectively. Based on this, complex singular spectrum analysis is performed on the original complex polar shift sequence to obtain the Chandler term component and the annual term component. Complex Fourier bandpass filtering is then applied to the original complex polar shift sequence, the Chandler term component, and the annual term component to cleanly separate the signal in the target frequency band. Finally, to eliminate the system offset caused by the filter on amplitude and phase, and to ensure that the filtering result is strictly consistent with the signal structure extracted by complex singular spectrum analysis in the complex plane, complex least squares gain matching is performed on the filtering result. This achieves robust extraction and continuous reconstruction of the periodic term in the polar shift time series, effectively reducing endpoint distortion and frequency offset effects.

[0125] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the technical solutions of the embodiments of the present invention.

Claims

1. A method for extracting and reconstructing optimization of polar shift periodic terms, characterized in that, include: Step S1: Obtain the polar motion observation sequence and construct the original complex polar motion sequence based on the polar motion observation sequence; Step S2: Perform complex singular spectrum analysis on the original complex polar shift sequence to obtain the Chandler term component and the annual term component; Step S3: Perform complex Fourier bandpass filtering on the original complex polar shift sequence, the Chandler term component, and the annual term component, respectively. Step S4: Perform complex least squares gain matching between the filtered results of the original complex polar shift sequence, Chandler term component, and annual term component and the Chandler term component and annual term component obtained from complex singular spectrum analysis to obtain the final polar shift periodic term.

2. The method for extracting and reconstructing polar shift periodic terms as described in claim 1, characterized in that, Step S1 includes: Step S11: Obtain the polar motion observation sequence under the same reference frame and unified epoch; Step S12: Extract the polar motion components of the two channels from the polar motion observation sequence; Step S13: Construct the original complex polar shift sequence based on the polar shift components of the dual channels.

3. The method for extracting and reconstructing optimization of polar shift periodic terms as described in claim 1, characterized in that, Step S2 includes: Step S21: Construct a complex polar shift trajectory matrix based on the original complex polar shift sequence; Step S22: Perform singular value decomposition on the polar motion complex trajectory matrix to obtain multiple singular values; Step S23: Perform a diagonal average operation on the matrix corresponding to each singular value to obtain the reconstructed components; Step S24: Perform a fast Fourier transform on the reconstructed components to obtain the spectrum corresponding to each element in the reconstructed components; Step S25: Perform spectrum main peak detection on the spectrum corresponding to each element in the reconstructed component to obtain the peak frequency and calculate the spectrum energy; Step S26: Selecting paired reconstructed components using the peak frequency of each element in the reconstructed components; Step S27: Reconstruct the anniversary item component and the Chandler item component based on the selected reconstruction component; Step S28: Perform positive spectrum estimation on the reconstructed Chandler term components and annual term components, find the set of spectral peaks within the reference frequency window, select the main peaks, and construct a passband for each main peak.

4. The method for extracting and reconstructing polar shift periodic terms as described in claim 3, characterized in that, In step S26, if the relative deviation of the peak frequencies of any two elements in the reconstructed components does not exceed the preset relative tolerance, then the two elements are considered to constitute a pair of reconstructed components. The cumulative spectral energy value of the pair of reconstructed components is calculated, and all the pair of reconstructed components are sorted according to the magnitude of the cumulative spectral energy value. Several pair of reconstructed components with the highest cumulative spectral energy value are selected.

5. The method for extracting and reconstructing polar shift periodic terms as described in claim 1, characterized in that, Step S3 includes: Step S31: Construct the complex response function of the Chandler term using the frequency band range of the Chandler term, and construct the complex response function of the annual term using the frequency band range of the annual term; Step S32: The original complex polar shift sequence is subjected to complex Fourier bandpass filtering using the Chandler term complex response function and the annual term complex response function, respectively, to obtain the Chandler term sequence and the annual term sequence of the filtered original complex polar shift sequence. Step S33: The Chandler term complex response function and the annual term complex response function are used to perform complex Fourier bandpass filtering on the Chandler term components and annual term components obtained after complex singular spectrum analysis, respectively, to obtain the Chandler term sequence of the filtered Chandler term components and the annual term sequence of the annual term components.

6. The method for extracting and reconstructing polar shift periodic terms as described in claim 5, characterized in that, In step S31, both the Chandler term complex response function and the anniversary term complex response function are constructed using a soft-edge transition design. The formula for the complex response function is as follows: ， in, Represents the complex response function; This indicates the calculation of cosine. Indicates frequency; Indicates the starting point of the left soft edge; Indicates the start of the left passband; Indicates the end of the right passband; This indicates the end point of the right soft edge.

7. The method for extracting and reconstructing polar shift periodic terms as described in claim 6, characterized in that, In step S32, the original complex polar shift sequence is subjected to complex Fourier bandpass filtering, and the corresponding formula is: ， In the formula, Represents the primitive complex polar shift sequence. and These represent the annual term sequence and the Chandler term sequence obtained by directly performing complex Fourier bandpass filtering on the original complex polar shift sequence, respectively. Indicates the inverse Fourier transform; This represents the complex response function of the Chandler term constructed using the frequency band range of the Chandler term. This represents the complex response function of the annual term constructed using the frequency band range of the annual term.

8. The method for extracting and reconstructing polar shift periodic terms as described in claim 6, characterized in that, The annual and Chandler term components obtained from complex singular spectrum analysis are then subjected to complex Fourier bandpass filtering, with the corresponding formulas as follows: ， In the formula, and These represent the annual term sequence and the Chandler term sequence obtained by performing complex singular spectrum analysis followed by complex Fourier bandpass filtering, respectively. and This represents the annual term component and the Chandler term component obtained from the complex singular spectrum analysis in step S2.

9. The method for extracting and reconstructing polar shift periodic terms as described in claim 1, characterized in that, In step S4, complex least squares gain matching is performed, and the corresponding formula is: ， in, This represents the final polar shift periodic term obtained by complex least squares gain matching; This indicates the objects to be subjected to gain matching, including the original complex polar shift sequence, the Chandler term component, and the filtering results of the annual term component; Indicates time; Represents the complex gain coefficient. The calculation formula is: ， in, This represents the annual and Chandler term components obtained from complex singular spectrum analysis; This indicates the filtering result obtained by directly applying complex Fourier bandpass filtering to the original complex polar shift sequence; * indicates complex conjugate. This indicates a modulo operation on a signal in a complex sequence.

10. A polar shift periodic term extraction and reconstruction optimization system, characterized in that, include: A polar shift dual-channel construction and complex domain representation module is used to acquire polar shift observation sequences and construct original complex polar shift sequences based on the polar shift observation sequences; The complex singular spectrum analysis module is used to perform complex singular spectrum analysis on the original complex polar shift sequence to obtain the Chandler term component and the annual term component; The complex Fourier bandpass filter module is used to perform complex Fourier bandpass filtering on the original complex polar shift sequence, the Chandler term component, and the annual term component, respectively. The complex least squares gain matching module is used to perform complex least squares gain matching on the filtered results of the original complex polar shift sequence, Chandler term component, and annual term component with the Chandler term component and annual term component obtained from complex singular spectrum analysis, so as to obtain the final polar shift periodic term.

Images