Method for automatic identification and removal of daily global drift disturbances in swarm satellite electron density data

Abstract

The application belongs to the field of satellite data processing, and particularly relates to a Swarm satellite electron density data daily global drift disturbance automatic identification and removal method, which comprises normalization processing; short-time Fourier transform is performed on the normalized data to obtain a time-frequency matrix, the time-frequency matrix is classified into strong disturbance orbit data and other orbit data by using a two-dimensional convolutional neural network, and a component containing daily global drift disturbance is identified; for the decomposition result of the strong disturbance orbit, the identified disturbance component is removed; for the disturbance component of the other orbit, a medium-strong disturbance interval is identified and located, and a disturbance latitude range is obtained according to a mapping relationship between an image coordinate and a latitude; for the identified medium-strong disturbance orbit, pre-normalization data thereof are obtained and subjected to non-negative matrix decomposition, a corresponding part in the disturbance component is cut off according to the disturbance latitude range, and the remaining part and the remaining components are added to serve as effective signals and reserved, so that automatic cleaning of massive satellite data is realized.

Description

Technical Field

[0001] This application belongs to the field of satellite data processing, specifically involving a method for automatic identification and removal of daytime global drift disturbances in Swarm satellite electron density data. Background Technology

[0002] Launched by the European Space Agency, the Swarm satellite constellation comprises three satellites with identical payloads. It observes various electromagnetic parameters in the ionosphere to explore the complex structure and dynamics of near-Earth space. The mission is currently operational and publicly provides its observations of core ionospheric parameters such as electron density and magnetic fields. Electron density, a crucial physical quantity of the ionosphere, has wide applications in ionospheric modeling, radio communications, and monitoring of seismic electromagnetic precursors, yielding many significant results. However, current research has not adequately addressed global drift disturbances in Swarm satellite electron density data, and effective removal methods have not yet been proposed. These disturbances, occurring during the day (local time), exhibit strong, varied, and broadband fluctuations that continuously move along orbital latitudes with local time. Their characteristics differ significantly from the geomagnetic storm disturbances and equatorial ionospheric irregularities addressed in mainstream research, severely degrading the quality of this data and thus affecting its scientific use.

[0003] Chinese Patent Publication No. CN 110673206 A discloses a method for detecting seismic anomalies in satellite magnetic field data based on nonnegative matrix factorization. The principle is to utilize the differences in time-frequency distribution characteristics between seismic anomalies (local disturbances) and space environment activities (global disturbances) in satellite magnetic field observation signals. A time-frequency matrix of the signal is constructed using short-time Fourier transform, and further, nonnegative matrix factorization is employed to decompose the mixed signal into components corresponding to different physical sources. This method adaptively separates local anomalous signals related to seismic activity from a complex global background field and identifies earthquake precursor information using statistical methods.

[0004] Chinese Patent Publication No. CN116299562A discloses a filtering method for correcting ionospheric errors in altimeter ranging. The principle is to utilize the inherent adaptability of the dual-frequency ionospheric observation signal, decomposing it into wave components of different scales through ensemble empirical mode decomposition. Based on this, an improved wavelet threshold function is used for overall filtering of the high-frequency components obtained from the decomposition, and stability is enhanced through multiple iterations and ensemble averaging. This method achieves overall suppression of high-frequency noise and systematic errors in the observation data, thereby improving the overall accuracy of the ionospheric correction results.

[0005] The aforementioned technologies have each proven effective in addressing tasks such as static bias calibration, time-frequency separable disturbances, and global noise suppression in satellite observation data. However, none of them can solve the problem of accurately identifying, locating, and separating local, dynamic, and superimposed diurnal drift disturbances in Swarm electron density data. Summary of the Invention

[0006] This application provides an automatic identification and removal method for daytime global drift disturbances in Swarm satellite electron density data, which solves the problem of accurate identification, location and separation of local, dynamic daytime drift disturbances that overlap with valid signals in Swarm electron density data.

[0007] This application provides a method for automatically identifying and removing daytime global drift disturbances in Swarm satellite electron density data, comprising: The selected Swarm satellite daytime observation orbital electron density data were sequentially subjected to wavelet decomposition to remove low-frequency high-amplitude background, as well as normalization processing of direction, length and amplitude. The time-frequency matrix is ​​obtained by performing a short-time Fourier transform on the normalized data. A two-dimensional convolutional neural network is used to classify the time-frequency matrix into strongly disturbed orbit data and other orbit data. Low-frequency target signals are added to the unnormalized data classified as strongly disturbed orbits. For strongly disturbed orbital data with added low-frequency target signals and other orbital data, short-time Fourier transform is performed to obtain their respective time-frequency amplitude matrices. Non-negative matrix decomposition is performed on the respective time-frequency amplitude matrices, and the high-frequency amplitude ratio of each decomposed component is calculated. Based on this, the components containing daytime global drift disturbances are identified. For the decomposition results of strongly disturbed orbits, the identified disturbance components are removed, and the data at the positions corresponding to the added low-frequency target signals in the remaining components are truncated, with the remaining part retained as the effective signal. For the disturbance components of other orbits, the disturbance components of all orbits within a single day are sorted longitudinally by time and plotted. The target detection model of the YOLO framework is used to identify and locate the medium-to-strong disturbance range, and the disturbance latitude range is obtained according to the mapping relationship between image coordinates and latitude. For the identified medium-to-strong perturbation orbits, their unnormalized data is obtained and non-negative matrix decomposition is performed. Based on the perturbation latitude range, the corresponding part of the perturbation component is truncated, and the remaining part is added to the other components to retain as the effective signal.

[0008] Furthermore, nonnegative matrix decomposition is performed on each time-frequency amplitude matrix, including: decomposing the time-frequency amplitude matrix into the sum of the products of R sets of frequency feature vectors and spatial weight vectors through nonnegative matrix decomposition, and reconstructing the corresponding time-domain signal; Based on the frequency eigenvectors obtained from the decomposition Calculate the high-frequency amplitude ratio of each group of decomposed components: , in, For high-frequency amplitude ratio, For the first The frequency characteristic component is at the ... The amplitude at each frequency point and These represent the start and end indices of the selected high-frequency band, respectively. This represents the total number of frequency points. The component group with the largest high-frequency amplitude ratio is selected as the disturbance component group.

[0009] Furthermore, the target detection model using the YOLO framework identifies and locates regions with strong disturbances, and obtains the range of disturbance latitudinals based on the mapping relationship between image coordinates and latitude, specifically including: A basic YOLO object detection framework was built, and a large number of manually annotated samples were used as the dataset for training. The positions of multiple medium-to-strong perturbation orbital perturbations were annotated in each image to obtain the object detection model. The image, which is drawn by sorting all the tracks in a single day longitudinally by time, is input into the trained YOLO object detection model to identify medium and strong perturbations and obtain the probability and bounding box position, i.e. the localization interval. Based on the height and width of the bounding box, obtain the latitude range corresponding to different medium-to-strong disturbance orbital perturbations.

[0010] Furthermore, the YOLO object detection model adopts a lightweight three-part architecture consisting of a backbone network, a feature fusion network, and a decoupled detection head. The backbone network includes an initial convolutional layer, three sets of C3k2 cross-stage local modules, and an SPPF fast spatial pyramid pooling module, used to extract multi-level features from the input image. The feature fusion network is based on a feature pyramid network and a path aggregation network structure, and enhances the model's ability to detect objects at different scales through upsampling and feature concatenation operations. The detection head is an anchorless decoupled detection head used to directly predict the center offset, width, height, and class confidence of the bounding box.

[0011] Furthermore, the predicted bounding box includes: the detector head directly predicts the offset of the bounding box relative to the grid cell, and the prediction of the center coordinates is constrained by the Sigmoid function to ensure that the center of the predicted box is located within the current grid. The decoding formula is: the predicted center coordinates are equal to the product of the original output of the network and the scaling factor, minus the offset, plus the coordinates of the top left corner of the current grid.

[0012] Furthermore, each set of frequency feature vectors and spatial weight vectors forms a frequency feature matrix and a weight feature matrix, respectively. The frequency feature matrix reflects the main frequency features of the original time-frequency amplitude matrix, and the weight feature matrix reflects the weight of the frequency distribution features along the orbital latitude. The number of decomposed features R satisfies the condition that it is less than the product of the number of rows and columns of the time-frequency amplitude matrix divided by the sum of the number of rows and columns.

[0013] Furthermore, the nonnegative matrix decomposition uses KL divergence as the objective function and adds a minimum determinant constraint. The objective function includes the divergence term between each element of the time-frequency amplitude matrix and its reconstructed value, as well as the minimum determinant constraint term. The optimization is performed by the multiplication rule and solved by the gradient descent method to obtain the multiplication update formula of the frequency feature matrix and the weight feature matrix.

[0014] Furthermore, low-frequency targeting signals are added to the unnormalized data classified as strongly disturbed orbital data, specifically including: For strongly disturbed orbital data that overlap with the time-frequency characteristics of useful signals, after classification and identification by a two-dimensional convolutional neural network, a low-frequency target signal is introduced as the anchor point for non-negative matrix decomposition. The low-frequency target signal is extracted from the low-frequency portion of the observation data that does not overlap with the disturbance and spliced ​​to the end of the data before normalization.

[0015] Furthermore, during the nonnegative matrix decomposition process, the low-frequency target signal is used as an anchor point to guide the decomposition algorithm to distinguish between perturbation components and useful signal components. The fixed low-frequency characteristics and known addition positions of the low-frequency target signal are used to calibrate the boundary between the effective signal components and perturbation components in the decomposition results.

[0016] Compared with the prior art, the advantages of this application are as follows: This application adopts a hierarchical processing strategy. For strongly disturbed orbital data that overlaps with useful signals in terms of time-frequency characteristics, it innovatively introduces low-frequency target signals as anchor points for NMF decomposition through CNN classification and identification, effectively separating strong disturbances from useful signals. For medium-strong disturbances that are locally located and constantly moving, based on their visual specificity, it combines the YOLO target detection algorithm to transform one-dimensional signal processing into a two-dimensional image recognition and localization problem, obtaining the latitudinal range of the disturbance and achieving local removal of the disturbance, thus improving the scientific usability of Swarm satellite electron density data. This application constructs a fully automated framework from preprocessing, classification, target addition, and localization removal, combining the advantages of deep learning (CNN, YOLO) and signal processing (wavelet decomposition, NMF), which can realize the automatic cleaning of massive satellite data. Attached Figure Description

[0017] Figure 1A flowchart illustrating an automatic localization and removal method for daytime global drift disturbances in Swarm satellite electron density data, provided in this application embodiment; Figure 2 The original data curve, preprocessed data curve, and normalized data curve of orbit 1 on January 1, 2015, provided for embodiments of this application; Figure 3 The performance evaluation graph of the classification model provided in the embodiments of this application is (a) the loss value, (b) the accuracy, (c) the true class, and (d) the true positive rate. Figure 4 The following is a diagram showing the NMF decomposition effect of strong disturbance in orbit 1 on January 1, 2015, provided for the embodiments of this application: (a) original data, (b) NMF1 disturbance component, and (c) NMF2 remaining component. Figure 5 The NMF decomposition effect diagram of the strong disturbance of orbit 15 on March 4, 2015, provided for the embodiment of this application, (a) original data, (b) NMF1 disturbance component, (c) NMF2 remaining component; Figure 6 The following is a diagram showing the NMF decomposition effect of adding a target signal to the strong disturbance of orbit 15 on March 4, 2015, provided in the embodiments of this application: (a) original data, (b) the added target signal, (c) the NMF1 disturbance component, and (d) the remaining NMF2 component. Figure 7 Visually recognizable effects of moderate to strong disturbances provided for embodiments of this application, (a) January 28, 2015, (b) February 25, 2015, (c) July 7, 2015 and (d) March 8, 2015; Figure 8 The following is a performance evaluation chart of the YOLO target detection model provided in the embodiments of this application: (a) is the loss value, (b) is the accuracy, and (c) is the average accuracy. Figure 9 This is a diagram illustrating the effect of strong perturbation in YOLO target detection and localization provided in an embodiment of this application. Figure 10 The following is a diagram showing the effect of removing medium-strong disturbances based on the positioning interval provided in the embodiments of this application: (a) original data, (b) final retained result, (c) NMF1 disturbance component, and (d) NMF2 remaining component. Detailed Implementation

[0018] 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.

[0019] An automatic identification and removal method for intraday global drift disturbances in Swarm satellite electron density data is proposed. This method combines two-dimensional convolutional neural network classification, a YOLO target detection model, and non-negative matrix factorization, employing a hierarchical processing framework: For situations where strong disturbances and valid signals have wide spatial overlap and similar frequency characteristics, making separation difficult, a two-dimensional convolutional neural network is first used to perform binary classification of the orbits to identify strong disturbances. Then, a target signal is added, using prior low-frequency information as anchor points to guide short-time Fourier transform (NMF) decomposition to accurately separate strong disturbances from valid signals. Finally, the disturbance components are directly removed. For situations where medium-to-strong disturbances and valid signals have local overlap, variable positions, and varying lengths, making overall removal difficult, one-dimensional signal processing is transformed into two-dimensional image recognition. An image is constructed from the disturbance components after the short-time Fourier transform of all orbits for a single day. The YOLO target detection model is used to automatically identify and locate the disturbances within their respective intervals. Finally, only the data corresponding to the disturbance intervals within the disturbance components are removed. This method achieves both effective separation of strong drift disturbances and local removal of moderate to strong disturbances, and the process is fully automated, significantly improving the scientific usability of Swarm satellite electron density data.

[0020] See Figure 1 As shown, a method for automatically identifying and removing daytime global drift disturbances in Swarm satellite electron density data includes: S101, the selected Swarm satellite daytime observation orbital electron density data are sequentially subjected to wavelet decomposition to remove low-frequency high-amplitude background, as well as normalization processing of direction, length and amplitude; S102, the normalized data is subjected to short-time Fourier transform to obtain the time-frequency matrix. The time-frequency matrix is ​​classified into strongly disturbed orbit data and other orbit data using a two-dimensional convolutional neural network. Low-frequency target signals are added to the data before normalization that are classified as strongly disturbed orbits. S103. For the strongly disturbed orbit data with added low-frequency target signal and other orbit data, perform short-time Fourier transform to obtain their respective time-frequency amplitude matrices. Perform non-negative matrix decomposition on their respective time-frequency amplitude matrices, calculate the high-frequency amplitude ratio of each decomposed component, and identify the component containing daytime global drift disturbance based on this. S104, for the decomposition results of strongly disturbed orbits, remove the identified disturbance components, and truncate the data at the positions corresponding to the added low-frequency target signals in the remaining components, and retain the remaining part as the effective signal; S105. For the disturbance components of other orbits, the disturbance components of all orbits within a single day are sorted longitudinally by time and plotted. The target detection model of the YOLO framework is used to identify and locate the medium-strong disturbance range. The disturbance latitude range is obtained according to the mapping relationship between image coordinates and latitude. S106. For the identified medium-to-strong disturbance orbit, obtain its unnormalized data and perform non-negative matrix decomposition. Based on the disturbance latitude range, truncate the corresponding part of the disturbance component, and then add the remaining part to the other components to retain the effective signal.

[0021] In step S101, the local time is selected as Time Using Swarm satellite electron density data as daytime observation data, wavelet decomposition was performed on the single-orbit data to remove low-frequency components as background, thus obtaining preprocessed single-orbit electron density data. Select [ , Data within a latitudinal range, normalized direction (south latitude to north latitude or north latitude to south latitude), length (to... The normalized orbital data is obtained by removing the shorter ones (and padding with zeros) and adjusting the amplitude (up to the range [0, 1]). .

[0022] Perform a short-time Fourier transform on the normalized orbit data obtained in step S101 to obtain a time-frequency amplitude matrix that reflects the time-frequency distribution characteristics of the signal. × This serves as the input for subsequent classification models.

[0023] A two-dimensional convolutional neural network (2D-CNN) is constructed, whose structure includes, in sequence: two identical convolutional blocks, global average pooling, a fully connected layer, and an output layer. Each convolutional block contains a convolutional layer, a batch normalization layer, a ReLU activation function layer, and a max-pooling layer. The specific parameters of the 2D convolutional neural network are shown in the table below: Table 1 shows the parameters of the 2D-CNN network.

[0024] Input size is × ×1, the network performs standardization preprocessing on the input data to accelerate training convergence. Key components of the network, such as convolution operations, activation functions, and loss functions, are shown below.

[0025] Convolution operations are used for feature extraction, employing a 5×5 convolution kernel. Its mathematical expression is: , in, Indicates the first The output feature map of the layer, and For learnable convolutional kernel weights and bias parameters, This is the input feature map for the previous layer. Spatial location coordinates, and These are the output and input channel indices, respectively. This represents the relative coordinate offset within the convolution kernel.

[0026] The activation function introduces nonlinearity using the Modified Linear Unit (ReLU), and its expression is as follows: ; This is the input value for the activation function.

[0027] The model training uses the cross-entropy loss function as the optimization objective, and its formula is as follows: ; in, For batch size, For the number of categories, For the true label of the sample, The predicted class probabilities are obtained from the model. The classification model is then derived. Data identified as strongly perturbed trajectories are taken from the pre-normalization data described in step S101. Low-frequency target signals were then added (after stitching) as data to be decomposed later. .

[0028] A classification model can be obtained by using a large number of manually labeled samples (divided into strongly perturbed orbits and other orbits) as a dataset for supervised training (divided proportionally into training, validation, and test sets). During training, the Adam optimizer can be used, with relevant hyperparameters set, such as 20 training epochs, a batch size of 32, and an initial learning rate of... Since this is a binary classification task, the output layer of the two-dimensional convolutional neural network obtains a normalized probability distribution through the Softmax function, and the category with the higher probability value is used as the final classification result.

[0029] Model performance is primarily measured by the F1 score, as shown in the following formula, supplemented by accuracy. and recall rate .

[0030] ; in, This represents the number of medium-to-strong disturbance regions correctly detected by the model. This represents the number of non-perturbation regions that the model falsely identifies as having moderate to strong perturbations. This represents the number of real-world areas with strong disturbances that were not detected by the model.

[0031] Specifically, low-frequency targeting signals are added to the unnormalized data classified as strongly disturbed orbital data, including: For strongly disturbed orbital data that overlap with the time-frequency characteristics of useful signals, after classification and identification by a two-dimensional convolutional neural network, a low-frequency target signal is introduced as the anchor point for non-negative matrix decomposition. The low-frequency target signal is extracted from the low-frequency portion of the observation data that does not overlap with the disturbance and spliced ​​to the end of the data before normalization.

[0032] Based on the classification results, the data identified as strongly disturbed orbits will be processed as the pre-normalization data in step S101. Low-frequency targeting signal added later (The following data is concatenated and will be used as subsequent data to be decomposed.) .in, Extract the low-frequency portion of the observation data that does not overlap with the disturbance.

[0033] Step S103 specifically includes: performing a short-time Fourier transform on the strongly perturbed orbit data with added target signal and other orbit data obtained in step S102 to generate a time-frequency amplitude matrix. It is decomposed into R sets of frequency eigenvectors using the nonnegative matrix factorization method. and spatial weight vector The sum of the products of each type of vector can form a matrix. The two matrices represent the frequency feature matrix and the position weight feature of the original signal, respectively, as shown in the following formula.

[0034] ; In the formula, for The non-negative matrix is ​​called the initial time-frequency amplitude matrix. for The non-negative matrix is ​​called the frequency characteristic matrix, which reflects the main frequency characteristics of the original time-frequency amplitude matrix. for The non-negative matrix is ​​called the weight feature matrix, which reflects the weights of the frequency distribution characteristics along the orbital latitude. The number of decomposition features should conform to the following: .

[0035] In the nonnegative matrix factorization process, the low-frequency target signal is used as an anchor point to guide the decomposition algorithm to distinguish between perturbation components and useful signal components. The fixed low-frequency characteristics and known addition position of the low-frequency target signal are used to calibrate the boundary between the effective signal components and perturbation components in the decomposition result.

[0036] KL divergence is selected as the objective function for nonnegative tensor decomposition, and a minimum determinant constraint is added to improve the stability of the decomposition results. Nonnegative matrix decomposition uses KL divergence as the objective function and adds a minimum determinant constraint. The objective function includes the divergence term between each element of the time-frequency amplitude matrix and its reconstructed value, as well as the minimum determinant constraint term. Optimization is performed through the multiplication rule, and the gradient descent method is used to solve the problem, resulting in the multiplication update formula for the frequency feature matrix and the weight feature matrix.

[0037] The final objective function is shown below: ; , , in, Representation matrix The One element, Indicates the first one Each element. The coefficients of the minimum determinant constraint.

[0038] By optimizing using the multiplication rule, setting up operators, and solving using gradient descent, the following update formula can be obtained: ; ; in, Represents the element-wise multiplication of a matrix (Hadamard product).

[0039] The overall process of nonnegative matrix factorization is shown below: Input: Time-frequency amplitude matrix Number of features R, maximum number of iterations Convergence threshold .

[0040] Output: Frequency characteristic matrix Spatial weight matrix .

[0041] Step 1: Randomly initialize the matrix , .

[0042] Step 2: Update the feature matrix using iterative formulas and .

[0043] Step 3: If the objective function is less than the preset convergence threshold or the number of iterations reaches the upper limit, stop the iteration and output the matrix. , Otherwise, go back to the previous step.

[0044] Among them, the number of decomposition features Maximum number of iterations and convergence threshold Selection and calculation are performed based on data characteristics and experience. Time-domain reconstruction is then performed on each set of decomposed components to obtain R time-domain signals.

[0045] Each set of frequency feature vectors and spatial weight vectors forms a frequency feature matrix and a weight feature matrix, respectively. The frequency feature matrix reflects the main frequency features of the original time-frequency amplitude matrix, and the weight feature matrix reflects the weight of the frequency distribution features along the orbital latitude. The number of decomposed features R satisfies the condition that it is less than the product of the number of rows and columns of the time-frequency amplitude matrix divided by the sum of the number of rows and columns.

[0046] Based on the decomposition Given a frequency eigenvector, calculate the high-frequency amplitude ratio of each decomposed component group, i.e., the proportion of high-frequency amplitude to the total amplitude of the entire frequency band. Since the frequency characteristics of global drift disturbances are broadband high-frequency signals, select the component group with the largest high-frequency amplitude ratio as the disturbance component group, and the corresponding time-domain signal is the disturbance component. The remaining component groups are referred to as the residual component groups, and the corresponding time-domain signals are the residual components. .

[0047] , In the formula, For high-frequency amplitude ratio, For the first The frequency characteristic component is at the ... The amplitude at each frequency point and These represent the start and end indices of the selected high-frequency band, respectively. The total number of frequency points is denoted as ; the component group with the largest high-frequency amplitude ratio is selected as the disturbance component group.

[0048] Step S104 specifically includes: for the decomposition results of the strongly perturbed orbit obtained from step S103, directly removing the perturbed components. , the remaining components Remove the low-frequency targeting signal added in step S102 Data at the corresponding position is used as a valid signal. It shall be retained.

[0049] Step S105 specifically includes: for the disturbance components of other orbits obtained from step S103, plotting all orbits within a single day along the time axis, with the vertical axis of each sub-plot representing the amplitude and the horizontal axis representing the latitude position along the orbit. , ],get The image of pixels is stored and used as input to the YOLO object detection model.

[0050] A target detection model using the YOLO framework identifies and locates regions with moderate to strong perturbations, and obtains the range of perturbation latitudinal dimensions based on the mapping relationship between image coordinates and latitude. Specifically, this includes: A basic YOLO object detection framework was built, and a large number of manually annotated samples were used as the dataset for training. The positions of multiple medium-to-strong perturbation orbital perturbations were annotated in each image to obtain the object detection model. The image, which is drawn by sorting all the tracks in a single day longitudinally by time, is input into the trained YOLO object detection model to identify medium and strong perturbations and obtain the probability and bounding box position, i.e. the localization interval. Based on the height and width of the bounding box, obtain the latitude range corresponding to different medium-to-strong disturbance orbital perturbations.

[0051] The YOLO object detection model employs a lightweight three-part architecture: a backbone network, a feature fusion network, and a decoupled detection head. The backbone network includes an initial convolutional layer, three sets of C3k2 cross-stage local modules, and an SPPF fast spatial pyramid pooling module, used to extract multi-level features from the input image. The feature fusion network is based on a feature pyramid network and a path aggregation network structure, and enhances the model's ability to detect objects at different scales through upsampling and feature concatenation operations. The detection head is an anchorless decoupled detection head used to directly predict the center offset, width, height, and class confidence of the bounding box.

[0052] Specifically, a YOLO11 object detection framework was built and the model was trained to automatically identify the location range of medium-to-strong perturbations in daily images. A lightweight three-part architecture of "Backbone + Neck + Head" was selected: The backbone network includes an initial convolutional layer, three sets of C3k2 cross-stage local modules, and an SPPF fast spatial pyramid pooling module, which is responsible for extracting multi-level features from the input image; the feature fusion network Neck is based on the FPN (Feature Pyramid Network) and PAN (Path Aggregation Network) structure, which enhances the model's ability to detect targets at different scales through upsampling and feature concatenation operations; the Head is an anchorless decoupled detection head that directly predicts the center offset, width and height of the bounding box, and class confidence.

[0053] The model's training and prediction are based on a bounding box prediction mechanism. It can detect the offset of the bounding box directly predicted by the head relative to the grid cells. Bounding box prediction includes: detecting the offset of the bounding box directly predicted by the head relative to the grid cells; and predicting the center coordinates using the Sigmoid function to ensure the predicted box center is within the current grid cell. The decoding formula is: the predicted center coordinates equal to the product of the Sigmoid function applied to the network's original output and the scaling factor, minus the offset, plus the coordinates of the top-left corner of the current grid cell. (Center coordinates) The prediction uses the Sigmoid function. To constrain the prediction box to ensure its center lies within the current grid, the general decoding formula can be expressed as: , in, It is the original output from the network. It is a scaling factor. It is the offset (default is 0.5). It is the coordinate of the top-left corner of the current grid.

[0054] The model training uses the CIoU loss function, which comprehensively considers the overlap area of ​​the bounding boxes, the distance between the center points, and the consistency of the aspect ratio. Its formula is as follows: , in, It is intersection, union, and comparison. It is Euclidean distance. It is the length of the diagonal of the smallest bounding rectangle. It is a weighting function. Used to measure aspect ratio consistency.

[0055] The dataset used to train the YOLO object detection model was generated manually (divided proportionally into training, test, and validation sets), with bounding boxes labeled for multiple medium-to-strong perturbation regions in each daily image. Relevant hyperparameters for training were set, such as: input image size 640×640, SGD (Stochastic Gradient Descent) optimizer, 200 training epochs, batch size 8, and initial learning rate 0.01.

[0056] The performance evaluation of the YOLO object detection model is based on mean precision (mAP), as shown in the following formula, which is the average precision under different recall rates, supplemented by the F1 score.

[0057] ; AP is the area under the Precision-Recall curve, which is usually achieved by summation.

[0058] Since the core of perturbation localization lies in accurately identifying the latitudinal range, and the vertical axis of the image corresponds to the time series, when evaluating the detection performance on a new dataset, only the degree of overlap of the bounding boxes in the horizontal (latitudinal) direction is considered, ignoring changes in the vertical direction. For predicted boxes... With real frame The formula for calculating its horizontal IoU is: ; ; All coordinates have been normalized. For each predicted bounding box, if its lateral IoU with any unmatched ground truth bounding box of the same category exceeds a preset threshold... If it is true, then it is recorded as a true example.

[0059] The image is input into a trained YOLO object detection model to identify medium to strong perturbations, obtaining the probability and bounding box position, i.e., the localization interval. Since the vertical position of different track patterns in the image is fixed, the horizontal pixel and latitude range […]. , A one-to-one correspondence can be established, based on the height and width of the bounding box, to obtain the latitude range corresponding to different moderate to strong disturbance orbits on that day. .

[0060] In step S106, for the medium-to-strong perturbation orbits identified in step S105, their unnormalized data is obtained. The nonnegative matrix decomposition is performed according to step S103. Based on the perturbation dimensional range obtained in step S105, the data corresponding to the perturbation component's dimensional range is truncated, and this truncated data is added to the remaining components to obtain the effective signal. reserve.

[0061] The model was trained using daytime electron density observation data from three Swarm satellites, Alpha (A), Bravo (B), and Charlie (C), in 2015 and 2016. The daytime electron density observation data from Alpha satellite in 2017 was used as a new independent dataset for validation, demonstrating the effectiveness of the proposed method.

[0062] For electron density data from the Swarm satellite constellation in 2015 (A, B, and C), 2016 (A and B), and 2017 (A), daytime observation data were selected from 6:00 AM to 5:00 PM local time. A five-level wavelet decomposition was performed on the single-orbit data using db4 as the basis function. The low-frequency components of the last level were removed as background to obtain preprocessed single-orbit electron density data. Data within the latitude range of [-60°, 60°] were selected, and the direction (from south latitude to north latitude) and length (to) were normalized. The normalized orbital data is obtained by removing the shorter ones (and padding with zeros) and adjusting the amplitude (up to the range [0, 1]). Taking orbit 1 on January 1, 2015 as an example, its raw data, preprocessed data, and normalized data are as follows: Figure 2 As shown in the figure. The results show that wavelet decomposition can effectively remove low-frequency, high-amplitude background from the data, and the normalized data can obtain data with fixed length, direction, and amplitude range (without changing shape).

[0063] A short-time Fourier transform (window length 64, step size 16, Hanning window function) is performed on the acquired normalized orbit data to obtain the time-frequency amplitude matrix reflecting the time-frequency distribution characteristics of the signal. This serves as the input for subsequent classification models.

[0064] A two-dimensional convolutional neural network (2D-CNN) is constructed, whose structure includes, in sequence: two identical convolutional blocks, global average pooling, a fully connected layer, and an output layer. Each convolutional block contains a convolutional layer, a batch normalization layer, a ReLU activation function layer, and a max-pooling layer. The specific parameters of the network are shown in the table below: Table 2 shows the parameters of the 2D-CNN network.

[0065] Input size is × ×1, the network performs standardization preprocessing on the input data to accelerate training convergence. Key components of the network, such as convolution operations, activation functions, and loss functions, are shown below.

[0066] Convolution operations are used for feature extraction, employing a 5×5 convolution kernel.

[0067] Due to the limited number of strongly perturbed orbits, to maintain dataset balance and prevent spurious results caused by excessively large differences in the number of categories, all strongly perturbed data from three satellites in 2015 and 2016, along with other representative perturbed data, were selected. These were manually labeled as strongly perturbed orbits and other orbits to construct a dataset for supervised training. The dataset was divided into training, validation, and test sets in a 14:3:3 ratio to obtain the classification model. Simultaneously, the orbits of satellite A throughout 2017 were labeled as a new, independent standard dataset for subsequent performance evaluation. The Adam optimizer was used during training, with the following key hyperparameter settings: 20 training epochs, batch size of 32, and initial learning rate. The probability is reduced to half its original value every 20 rounds, and the performance is evaluated on a validation set every 30 iterations. Since this is a binary classification task, the output layer of the convolutional neural network obtains a normalized probability distribution through the Softmax function, and the class with the higher probability value is used as the final classification result.

[0068] The effect of a two-dimensional convolutional neural network is as follows Figure 3 As shown. By Figure 3 (a) and Figure 3 As shown in (b), during model training, the loss and accuracy of the training and validation sets decrease and increase respectively with increasing iteration count, eventually stabilizing without significant overfitting. This demonstrates good convergence, achieving values ​​below 0.08 and above 97%. Figure 3 As shown in (c), the confusion matrix results indicate that the model achieved accuracy, precision, and recall of 98.47%, 98.53%, and 97.63% on the test set, respectively, with a misclassified sample rate of less than 2%, demonstrating stable and reliable classification performance. Figure 3 As shown in (d), the area under the ROC curve (AUC) is 0.996, indicating that the model has a very strong binary classification ability and can accurately identify and distinguish between strongly disturbed orbits and other orbits.

[0069] To further test the model's generalization ability and robustness on real-world unknown data, the trained classification model was directly applied to the new 2017 A* dataset for performance evaluation, as shown in Table 4. The results show that the model's macro-mean precision, recall, and F1 score all exceed 97% on independent datasets spanning multiple years, with an overall error rate of only 1.01%. For non-strongly perturbed orbit categories, the precision, recall, and F1 score all exceed 98.5%, while for strongly perturbed orbit categories, the corresponding metrics are 99.30%, 94.43%, and 96.80%, respectively. Overall, the model maintains high-precision classification performance.

[0070] Table 3 Classification parameters of the 2D-CNN model for new data

[0071] Based on the classification results, the data identified as strongly disturbed orbits will be processed according to the pre-normalization data described in step a. Low-frequency targeting signal added later (The data is then spliced ​​together, with a length of 5% to 10% of the original signal and an amplitude of approximately 1 to 2 times the original), to be used as data for subsequent decomposition. .in, Extracting from the low-frequency portion of the observation data that does not overlap with the disturbance ( The maximum frequency of the data is 1 Hz.

[0072] Short-time Fourier transform (window length 64, step size 16, Hanning window function) is performed on the acquired strongly perturbed orbit data with added target signals and other orbit data to generate a time-frequency amplitude matrix. (size is) The R sets of frequency eigenvectors are decomposed using the nonnegative matrix factorization method. and spatial weight vector The sum of the products of each type of vector can form a matrix, and the two matrices respectively represent the frequency feature matrix and the position weight feature of the original signal.

[0073] Based on experience and practical situation, the number of decomposition features will be... Set to 2, maximum number of iterations Set the convergence threshold to 1000. Let the mean of all elements of the matrix to be decomposed be 0.0001 times. Time-domain reconstruction of each set of decomposed components yields two time-domain signals.

[0074] Based on the decomposition Given a frequency eigenvector, calculate the high-frequency amplitude ratio of each decomposed component group, i.e., the proportion of high-frequency amplitude to the total amplitude of the entire frequency band. Since the frequency characteristics of global drift disturbances are broadband high-frequency signals, select the component group with the largest high-frequency amplitude ratio as the disturbance component group, and the corresponding time-domain signal is the disturbance component. The other group is the residual component group, and the corresponding time-domain signal is the residual component. .

[0075] Taking orbits from January 1, 2015, and March 4, 2015, as examples, Figure 4 , Figure 5 and Figure 6 The separation effects of nonnegative matrix factorization (NMF) on two types of perturbations are demonstrated, as well as the separation effect of NMF on strong perturbations after adding a target signal. Figure 4 It can be seen that, Figure 4 The original data in (a) contains both low-frequency signals and high-frequency medium-strength disturbance signals. After decomposition, it is as follows: Figure 4 In (b) NMF1, the disturbance component contains medium-to-strong disturbance signals from the original data, such as... Figure 4 In (c)NMF2, the remaining component contains the unperturbed signal from the original data. Figure 5 and Figure 6 It can be seen that if non-negative matrix decomposition is performed directly on the strongly disturbed orbit data without the addition of a target signal, the disturbance cannot be separated, and both of the resulting NMF components contain high-frequency disturbance signals. However, if a low-frequency target signal is added and non-negative matrix decomposition is performed on the new data, the target signal can serve as a guide, and the separation of high-frequency strong disturbance signals and low-frequency residual signals can be achieved effectively.

[0076] Taking orbit 15 on March 4, 2015 as an example, such as Figure 6 As shown, Figure 6In (a), after adding a target signal to the strongly disturbed orbital data, non-negative matrix decomposition can effectively separate high-frequency disturbances into the NMF1 disturbance component. Figure 6 In (b) of the above, it is directly removed, retaining the NMF2 component containing the low-frequency non-disturbance signal. Figure 6 (c) can achieve the removal of strong disturbances.

[0077] For the disturbance components of other acquired orbits, all orbits within a single day are plotted longitudinally along time. Each sub-plot has its vertical axis representing the amplitude and its horizontal axis representing the latitude position along the orbit [-60°, 60°]. A pixel-level PNG image is generated and stored, serving as input to the YOLO object detection model.

[0078] On January 28, 2015 ( Figure 7 (a) in the middle, February 25 ( Figure 7 (b) of the middle, July 7 ( Figure 7 (c) in the middle and March 8 ( Figure 7 For example, (d) in the middle, Figure 7 This displays a PNG image constructed from the daytime orbital perturbation components. Figure 7 The first three images (a), (b), and (c) are all perturbation images, demonstrating three typical forms of moderate to strong perturbations. This indicates a significant difference in visual effects between moderately strong perturbation data and undisturbed data on the same day, and a correlation exists between the locations of all orbital perturbations on the same day. This feature is suitable for target detection algorithms.

[0079] A dataset was constructed using all daytime orbital data from 2015 (satellites A, B, and C) and 2016 (satellites A and B) on days with no strong disturbances for supervised training. Multiple moderately to strongly disturbed regions in each daily image were labeled with bounding boxes, and the dataset was divided into training, validation, and test sets in an 8:1:1 ratio to obtain the object detection model. Simultaneously, other orbits obtained from the classification of satellite A throughout 2017 were labeled as a new, independent standard dataset to test the robustness and reliability of the model. The key hyperparameters for training were set as follows: input image size 640×640, SGD (stochastic gradient descent) optimizer, 200 training epochs, batch size 8, initial learning rate 0.01, and data augmentation during training.

[0080] The performance of the YOLO object detection model is as follows: Figure 8 As shown. By Figure 8 As shown in (a), during model training, the losses on both the training and validation sets exhibit a rapid decrease followed by a stabilization. After approximately 100 iterations, the two curves essentially overlap and stabilize at relatively low values, indicating that the model has good convergence in the object detection task and that no significant overfitting is observed. Figure 8 As shown in (b), the PR curve results indicate that the model achieves an AP value of 0.988 for targets in moderately disturbed orbits, maintaining extremely high precision throughout the full recall range, demonstrating excellent recall-precision balance capabilities. Figure 8 As shown in (b) and (c), the mAP performance under different IoU thresholds indicates that the model has an mAP of 0.988 when IoU=0.5 and still reaches 0.910 when IoU=0.8, which verifies its high target localization ability, that is, the overlap between the predicted bounding box and the actual disturbance range is high.

[0081] To visually illustrate the specific localization performance of the YOLO object detection model, four representative days were selected from the test set, such as... Figure 9 As shown, February 18th and September 30th, 2015, respectively, demonstrate the good detection performance when medium-to-strong perturbations occur on one side and both sides. It can be seen that the IoU values ​​for all samples are above 0.8, with most even exceeding 0.9, indicating excellent localization accuracy in scenarios with clear perturbation morphology and regular position. September 9th, 2015, and December 26th, 2016 (Satellite B) show two scenarios: false positives and false negatives. The former is due to false positives caused by the similarity between the background high-frequency perturbation and the target perturbation morphology, while the latter is due to false negatives caused by weak perturbation intensity and insufficient feature recognition. The performance in these scenarios still needs further optimization.

[0082] To further test the model's generalization ability and robustness on real-world unknown data, the trained classification model was directly applied to the 2017 AstroA dataset (2017A) for performance evaluation, as shown in Table 5. The results show that with IoU of 0.5 and 0.8, the model's average recall, accuracy, and F1 score on the 2017A dataset are close to or exceed 95%, with no significant performance degradation compared to the original test set, validating the generalization ability and engineering practicality of the proposed method.

[0083] Table 4 Localization parameters of the YOLO object detection model on the new dataset

[0084] Taking orbit 31 on February 18, 2015 as an example, the prediction results for that day are sorted by vertical position to obtain the corresponding prediction box. Based on the horizontal center position and width, the left and right boundaries are calculated and mapped to latitude. From 60 to 60, we can obtain a latitude range of 28.18 to 47.43, which can be used as the range for removing orbital disturbances.

[0085] For the identified moderately and strongly disturbed orbits, obtain their unnormalized data. , Perform nonnegative matrix decomposition. Based on the perturbation dimensional range obtained in step e, truncate the data corresponding to the perturbation component's dimensional range, and add this truncated data to the remaining components to obtain the effective signal. reserve.

[0086] Taking orbit 31 on February 18, 2015 as an example, Figure 10 This demonstrates the effect of removing strong perturbations from the original data. (By...) Figure 10 (a) and Figure 10 As shown in (c), nonnegative matrix decomposition can separate medium-strong disturbances into NMF1 disturbance components. Based on the obtained latitudinal range, the values ​​within the disturbance range are set to 0 (e.g., ...). Figure 10 (as shown in box (c) in the image), and the remaining NMF2 components (as shown in the image). Figure 10 Adding the results shown in (d) in the diagram yields the final retained result, as shown in the diagram. Figure 10 As shown in (b) of the diagram.

[0087] 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 automatically identifying and removing daytime global drift disturbances in Swarm satellite electron density data, characterized in that, include: The selected Swarm satellite daytime observation orbital electron density data were sequentially subjected to wavelet decomposition to remove low-frequency high-amplitude background, as well as normalization processing of direction, length and amplitude. The time-frequency matrix is ​​obtained by performing a short-time Fourier transform on the normalized data. A two-dimensional convolutional neural network is used to classify the time-frequency matrix into strongly disturbed orbit data and other orbit data. Low-frequency target signals are added to the unnormalized data classified as strongly disturbed orbits. For strongly disturbed orbital data with added low-frequency target signals and other orbital data, short-time Fourier transform is performed to obtain their respective time-frequency amplitude matrices. Non-negative matrix decomposition is performed on the respective time-frequency amplitude matrices, and the high-frequency amplitude ratio of each decomposed component is calculated. Based on this, the components containing daytime global drift disturbances are identified. For the decomposition results of strongly disturbed orbits, the identified disturbance components are removed, and the data at the positions corresponding to the added low-frequency target signals in the remaining components are truncated, with the remaining part retained as the effective signal. For the disturbance components of other orbits, the disturbance components of all orbits within a single day are sorted longitudinally by time and plotted. The target detection model of the YOLO framework is used to identify and locate the medium-to-strong disturbance range, and the disturbance latitude range is obtained according to the mapping relationship between image coordinates and latitude. For the identified medium-to-strong perturbation orbits, their unnormalized data is obtained and non-negative matrix decomposition is performed. Based on the perturbation latitude range, the corresponding part of the perturbation component is truncated, and the remaining part is added to the other components to retain as the effective signal.

2. The method for automatic identification and removal of daytime global drift disturbances in Swarm satellite electron density data according to claim 1, characterized in that, Perform nonnegative matrix decomposition on each time-frequency amplitude matrix, including: decomposing the time-frequency amplitude matrix into the sum of products of R sets of frequency eigenvectors and spatial weight vectors through nonnegative matrix decomposition, and reconstructing the corresponding time-domain signal; Based on the frequency eigenvectors obtained from the decomposition Calculate the high-frequency amplitude ratio of each group of decomposed components: ， in, For high-frequency amplitude ratio, For the first The frequency characteristic component is at the ... The amplitude at each frequency point and These represent the start and end indices of the selected high-frequency band, respectively. This represents the total number of frequency points. The component group with the largest high-frequency amplitude ratio is selected as the disturbance component group.

3. The method for automatic identification and removal of daytime global drift disturbances in Swarm satellite electron density data according to claim 1, characterized in that, The target detection model using the YOLO framework identifies and locates regions with strong disturbances, and obtains the range of disturbance latitudinals based on the mapping relationship between image coordinates and latitude. Specifically, this includes: A basic YOLO object detection framework was built, and a large number of manually annotated samples were used as the dataset for training. The positions of multiple medium-to-strong perturbation orbital perturbations were annotated in each image to obtain the object detection model. The image, which is drawn by sorting all the tracks in a single day longitudinally by time, is input into the trained YOLO object detection model to identify medium and strong perturbations and obtain the probability and bounding box position, i.e. the localization interval. Based on the height and width of the bounding box, obtain the latitude range corresponding to different medium-to-strong disturbance orbital perturbations.

4. The method for automatic identification and removal of daytime global drift disturbances in Swarm satellite electron density data according to claim 3, characterized in that, The YOLO object detection model employs a lightweight three-part architecture: a backbone network, a feature fusion network, and a decoupled detection head. The backbone network includes an initial convolutional layer, three sets of C3k2 cross-stage local modules, and an SPPF fast spatial pyramid pooling module, used to extract multi-level features from the input image. The feature fusion network is based on a feature pyramid network and a path aggregation network structure, and enhances the model's ability to detect objects at different scales through upsampling and feature concatenation operations. The detection head is an anchorless decoupled detection head used to directly predict the center offset, width, height, and class confidence of the bounding box.

5. The method for automatic identification and removal of daytime global drift disturbances in Swarm satellite electron density data according to claim 4, characterized in that, The predicted bounding box includes: the detector head directly predicts the offset of the bounding box relative to the grid cell, and the prediction of the center coordinates is constrained by the Sigmoid function to ensure that the center of the predicted box is located within the current grid. The decoding formula is: the predicted center coordinates are equal to the product of the original output of the network and the scaling factor, minus the offset, plus the coordinates of the top left corner of the current grid.

6. The method for automatic identification and removal of daytime global drift disturbances in Swarm satellite electron density data according to claim 2, characterized in that, Each set of frequency feature vectors and spatial weight vectors forms the frequency feature matrix and weight feature matrix, respectively. The frequency feature matrix reflects the main frequency features of the original time-frequency amplitude matrix, while the weight feature matrix reflects the weights of the frequency distribution features along the orbital latitude. The number of decomposition features R is less than the product of the number of rows and columns of the time-frequency amplitude matrix divided by the sum of the number of rows and columns.

7. The method for automatic identification and removal of daytime global drift disturbances in Swarm satellite electron density data according to claim 6, characterized in that, The nonnegative matrix decomposition uses KL divergence as the objective function and adds a minimum determinant constraint. The objective function includes the divergence term between each element of the time-frequency amplitude matrix and its reconstructed value, as well as the minimum determinant constraint term. The optimization is performed by the multiplication rule and solved by the gradient descent method to obtain the multiplication update formula of the frequency feature matrix and the weight feature matrix.

8. The method for automatic identification and removal of daytime global drift disturbances in Swarm satellite electron density data according to claim 7, characterized in that, Add low-frequency targeting signals to the unnormalized data classified as strongly disturbed orbital data, specifically including: For strongly disturbed orbital data that overlap with the time-frequency characteristics of useful signals, after classification and identification by a two-dimensional convolutional neural network, a low-frequency target signal is introduced as the anchor point for non-negative matrix decomposition. The low-frequency target signal is extracted from the low-frequency portion of the observation data that does not overlap with the disturbance and spliced ​​to the end of the data before normalization.

9. The method for automatic identification and removal of daytime global drift disturbances in Swarm satellite electron density data according to claim 8, characterized in that, In the nonnegative matrix factorization process, the low-frequency target signal is used as an anchor point to guide the decomposition algorithm to distinguish between perturbation components and useful signal components. The fixed low-frequency characteristics and known addition position of the low-frequency target signal are used to calibrate the boundary between the effective signal components and perturbation components in the decomposition result.

Images