Explainable spatio-temporal graph anomaly detection method based on graph filtering

Abstract

The present disclosure provides a kind of based on graph filtering's explainable spatio-temporal graph anomaly detection method, the method comprises: the original telemetry data of received satellite power supply system is preprocessed to obtain normal sample;Graph filtering-based explainable spatio-temporal graph anomaly detection network is constructed;The normal sample is used to train the explainable spatio-temporal graph anomaly detection network, when the number of iterations reaches the maximum iteration number set or the total loss reduction reaches the set threshold value, the explainable spatio-temporal graph anomaly detection network training is completed;Given any received satellite power supply system telemetry data, the explainable spatio-temporal graph anomaly detection network trained is used to carry out anomaly detection and explainable analysis.The present disclosure simultaneously captures spatio-temporal information, provides explainability using graph filtering technology and ingenious network design, is robust to high noise and data loss, can realize high-performance, credible, practical satellite power supply system telemetry signal anomaly detection.

Description

Technical Field

[0001] This disclosure belongs to the field of multidimensional time series anomaly detection, and particularly relates to an interpretable spatiotemporal graph anomaly detection method based on graph filtering. Background Technology

[0002] In recent years, data-driven methods (especially deep learning methods) have been applied to various fields, and using artificial intelligence to solve scientific problems has become a new paradigm. In the field of spacecraft telemetry data processing, numerous data-driven methods have also been used for anomaly detection in satellite power system telemetry signals. Currently, some unresolved issues remain in existing technologies.

[0003] (1) Ignoring spatial information: Most of the current research in this field focuses only on extracting time information and ignores the extraction of spatial information. However, there are complex and obvious mutual coupling relationships between the various sensors in the satellite power system. If the physical, logical and pattern relationships between these sensor nodes can be utilized, the performance of the anomaly detection method can be greatly improved.

[0004] (2) Lack of interpretability: Satellite launch and maintenance costs are enormous, and humans would not entrust such a costly task to an unexplainable black-box model. Meanwhile, interpretability provides traceability, enabling operators to understand their decision-making processes and rationale, accelerating the analysis process of the operations team, and facilitating subsequent anomaly localization and attribution. Furthermore, interpretability can also help people understand its internal mechanisms, providing guidance for model improvement.

[0005] (3) Regarding data characteristics: In order to develop truly practical anomaly detection methods, it is necessary to pay attention to the data characteristics of satellite power system telemetry data, which is also the difficulty of anomaly detection: high noise and missing data. Due to the influence of temporary sensor errors and the complex space transmission environment such as space particle radiation and electromagnetic interference, spacecraft telemetry data always has extremely high noise. At the same time, due to insufficient ground receiving station coverage, space electromagnetic and radiation environment, limited communication channel bandwidth, network channel congestion, limited buffer register storage space, and physical damage to components, satellite telemetry data often has partial missing data.

[0006] Therefore, there is an urgent need in this field for an interpretable anomaly detection method that can simultaneously capture spatiotemporal information and is robust to noise and missing information.

[0007] The information disclosed in the background section is only intended to enhance the understanding of the background of the present invention, and therefore may contain information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0008] The purpose of this disclosure is to develop an interpretable anomaly detection method that can simultaneously capture spatiotemporal information and is robust to noise and missing data. This method utilizes spatiotemporal graph technology to enable the interpretable spatiotemporal graph anomaly detection network to simultaneously capture spatiotemporal information; it provides interpretability to the network by developing a novel adaptive graph filter as the foundation for the network architecture; and it makes the interpretable spatiotemporal graph anomaly detection network robust to high noise and missing data through an autoencoder architecture and adversarial training strategies. To achieve the above objectives, this disclosure adopts the following scheme.

[0009] An interpretable spatiotemporal graph anomaly detection method based on graph filtering for anomaly detection in satellite power system telemetry signals includes the following steps:

[0010] S100: Preprocessing of telemetry signals from the satellite power system, including nearest neighbor interpolation, moving average filtering, maximum-minimum normalization, and graph structure calculation;

[0011] S200: Construct an interpretable spatiotemporal graph anomaly detection network based on graph filtering;

[0012] S300: Train an interpretable spatiotemporal graph anomaly detection network using already received normal samples;

[0013] S400: Given any received satellite power system telemetry data, use the trained interpretable spatiotemporal graph anomaly detection network to perform anomaly detection and interpretable analysis.

[0014] Preferably, in step S100, the preprocessing includes the following sub-steps:

[0015] S101: The received raw telemetry data of the satellite power system is represented as: {l′1, ..., l′ N}, l′ i It is a sequence of length len, representing the observation data received by the i-th sensor, where N is the number of sensors, i ranges from 1 to N, and N is a natural number greater than or equal to 1.

[0016] S102: For {l′1, ..., l′...} N Missing values ​​in} are initially filled using the nearest neighbor imputation method, that is, l′ i The missing value is filled by the preceding or following non-missing value of the missing position, resulting in the padded data representation as: {l″1,...,l″ N}

[0017] S103: For {l″1,…,l″ N Preliminary denoising is performed using a moving average filter, resulting in smoothed data represented as {l″′1,..,l″′ NThe calculation formula is:

[0018]

[0019] Where a is the set sliding window length and b is the number of sliding windows.

[0020] S104: For {l″′1,...,l″′} N Performing a maximum-minimum normalization operation, the normalized data is represented as {l1, ..., l...} N The calculation formula is:

[0021]

[0022] s105: Using {l1, ..., l N Perform graph structure calculations:

[0023]

[0024] Among them o i Representing the i-th sensor node, all sensor nodes form a graph node set o = {o1, ..., o2}. N},(o i o j ) represents node o i and node o j Establish edges between each other, and all edges form the edge set E = {(o1, o2, o3)} of the graph. i ), ..., (o j o N The graph structure is represented as G = (O, E), sim ij =d(l i , l j ) represents the observed sequence l i and l j The Euclidean distance between them, TopM(·) represents the sim ij In the current case i, the set of nodes corresponding to the first m maximum values, A represents the adjacency matrix, used to store the edge information of the graph structure, and can be used to further calculate the Laplacian matrix and the normalized Laplacian matrix:

[0025]

[0026] Where D = diag(d1, ..., d N ), d i =∑ j A ij d i =∑ j A ij , represents the adjacency matrix.

[0027] S106: The processed data l is divided into multiple subsequences using the sliding window method, with each subsequence representing a sample X.

[0028] Preferably, in step S200, the network construction includes the following sub-steps:

[0029] S201: Establish an adaptive quadratic approximation graph filter AQAGF(·):

[0030]

[0031] Where K is the order of the empirically defined polynomial, a k b k c k Both z and X are network learning parameters, X is the sample matrix, I is the identity matrix, and A and U are the eigenvalue matrix and Fourier transform matrix of the Laplacian matrix, respectively. Their calculation methods are as follows:

[0032] S202: Establishing an interpretable spatiotemporal graph anomaly detection network:

[0033]

[0034] Where E(·), G(·), and D(·) represent the encoder, generator, and discriminator, respectively. It is a randomly generated high-dimensional space representation vector, v, and These are the three discrimination scores output by the interpretable spatiotemporal graph anomaly detection network. It is the reconstructed sample output of an interpretable spatiotemporal graph anomaly detection network. It is by The random reconstructed sample matrix generated by the generator, GRL(·) is the gradient reversal function, which means that during network training, the variables are not changed when calculating the loss, but a negative sign is added before the gradient of the variable when the gradient is backpropagated.

[0035] Preferably, in step S300, the network training includes the following sub-steps:

[0036] S301: Input each received normal sample X into the interpretable spatiotemporal graph anomaly detection network to obtain each v.

[0037] S302: Calculate the backpropagation loss value based on the output of the interpretable spatiotemporal graph anomaly detection network. The calculation formula is:

[0038] L=α·L R +L D ,

[0039]

[0040]

[0041] Among them, L, L R and L D These are the total training loss, reconstruction loss, and discrimination loss, respectively, with α and ε being balance coefficients.

[0042] S303: Backpropagate the calculated total loss L to train the interpretable spatiotemporal graph anomaly detection network;

[0043] S304: Repeat steps S301 to S303. When the number of iterations reaches the set maximum number of iterations or the decrease in the total loss L reaches the set threshold, the training of the interpretable spatiotemporal graph anomaly detection network is complete.

[0044] Preferably, in step S400, the practical application includes the following sub-steps:

[0045] S401: Input any received satellite power system telemetry data into the trained interpretable spatiotemporal graph anomaly detection network to obtain the anomaly score matrix, calculated as follows:

[0046]

[0047] Among them, S i，t Here is the anomaly score matrix, and δ is the balance parameter;

[0048] S402: Set the maximum value of the anomaly scores of each sensor at each time point as the final anomaly score determined at that time point. The calculation formula is:

[0049] S403: Based on experience, a threshold is set, and any anomaly scores at each time point exceeding the threshold are marked as anomalies, thus completing anomaly detection. Simultaneously, interpretable analysis is performed on the frequency response function and other data obtained from the network.

[0050] Preferably, in step S202, the encoder, generator, and discriminator each include:

[0051]

[0052] In the above formula, MLP(x) is a multilayer perceptron model, expressed as:

[0053] MLP(x)=σ(b MLP +W MLP x)

[0054] Where σ(·) represents the sigmoid nonlinear activation function, b MLPand W MLP These are the bias terms and weight matrices of the multilayer perceptron model.

[0055] GRU(x) is a gated recurrent neural network unit, expressed as:

[0056] r t =σ(W rg ·[h t-1 x t ]+b rg )

[0057] u t =σ(W ug ·[h t-1 x t ]+b ug )

[0058]

[0059]

[0060] Where r represents resetting the gate's output value, u represents updating the gate's output value, h represents the hidden state, and x represents the input value. The latest hidden state is represented by the subscripts t and t-1, which represent the current signal segment and the previous signal segment, respectively. W represents the weight matrix, b represents the bias term, rg represents the reset gate, ug represents the update gate, tanh(·) represents the hyperbolic tangent function, and ⊙ represents the Hadamard product.

[0061] Compared with the prior art, this disclosure has the following beneficial technical effects:

[0062] 1. It can capture spatiotemporal information simultaneously, not only extracting the time characteristics such as trends and periods of telemetry time series signals, but also utilizing the interrelationship characteristics between various sensors.

[0063] 2. It is capable of performing interpretable anomaly detection, and its anomaly detection results are trustworthy and can support subsequent anomaly localization and attribution.

[0064] 3. It is applicable to telemetry signals from actual satellite power systems and is robust to high noise and missing data. Attached Figure Description

[0065] The accompanying drawings illustrate exemplary embodiments of the present disclosure and, together with the description thereof, serve to explain the principles of the present disclosure. These drawings are included to provide a further understanding of the present disclosure and are incorporated in and constitute a part of this specification.

[0066] Figure 1This is a flowchart illustrating an interpretable spatiotemporal graph network based on graph filtering for anomaly detection in telemetry signals of a satellite power system, provided in one embodiment of this disclosure.

[0067] Figure 2 This is a schematic diagram of the principle of a spatial information capture module, namely an adaptive quadratic approximation graph filter, provided in another embodiment of this disclosure;

[0068] Figure 3 This is a schematic diagram illustrating the training process of an interpretable spatiotemporal graph anomaly detection network based on graph filtering, provided in another embodiment of this disclosure.

[0069] Figure 4 This is a schematic diagram of a process for anomaly detection using a trained interpretable spatiotemporal graph anomaly detection network, provided in another embodiment of this disclosure.

[0070] Figure 5 This is a schematic diagram of dataset partitioning and representative channel visualization provided in another embodiment of this disclosure;

[0071] Figure 6 This is a schematic diagram of anomaly detection results for a representative channel provided in another embodiment of this disclosure;

[0072] Figure 7 This is a schematic diagram illustrating an interpretable result provided by another embodiment of this disclosure. Detailed Implementation

[0073] The following is in conjunction with the appendix Figures 1 to 7 The present disclosure will be further described in detail below with reference to the embodiments. It is to be understood that the specific embodiments described herein are for illustrative purposes only and are not intended to limit the scope of the disclosure. It should also be noted that, for ease of description, only the parts relevant to the present disclosure are shown in the accompanying drawings.

[0074] It should be noted that, where there is no conflict, the embodiments and features described in this disclosure can be combined with each other. The technical solutions of this disclosure will now be described in detail with reference to the accompanying drawings and embodiments.

[0075] Unless otherwise stated, the exemplary implementations / embodiments shown are to be understood as providing exemplary features of various details that provide ways in which the technical concepts of this disclosure can be implemented in practice. Therefore, unless otherwise stated, the features of various implementations / embodiments may be additionally combined, separated, interchanged and / or rearranged without departing from the technical concepts of this disclosure.

[0076] The use of crosshairs and / or shading in the accompanying drawings is generally used to clarify the boundaries between adjacent components. Thus, unless otherwise stated, the presence or absence of crosshairs or shading does not convey or indicate any preference or requirement for the specific material, material properties, dimensions, proportions, commonalities between the illustrated components, or any other characteristics, properties, etc., of the components. Furthermore, in the accompanying drawings, the dimensions and relative dimensions of components may be exaggerated for clarity and / or descriptive purposes. When exemplary embodiments can be implemented differently, a specific process sequence may be performed in a different order than that described. For example, two consecutively described processes may be performed substantially simultaneously or in the reverse order of their description. Furthermore, the same reference numerals denote the same components.

[0077] When a component is referred to as being "on" or "above" another component, "connected to," or "joined to" another component, the component may be directly on, directly connected to, or directly joined to the other component, or there may be intermediate components. However, when a component is referred to as being "directly on" another component, "directly connected to," or "directly joined to" another component, there are no intermediate components. Therefore, the term "connection" can refer to a physical connection, an electrical connection, etc., and may or may not have intermediate components.

[0078] For descriptive purposes, this disclosure may use spatial relative terms such as “below,” “under,” “below,” “down,” “above,” “above,” “higher,” and “side (e.g., in a “sidewall”)” to describe the relationship between one component and another component as shown in the accompanying drawings. In addition to the orientations depicted in the drawings, the spatial relative terms are also intended to encompass different orientations of the device during use, operation, and / or manufacture. For example, if the device in the drawings is flipped, a component described as “below” or “under” another component or feature would subsequently be positioned “above” said other component or feature. Thus, the exemplary term “below” can encompass both “above” and “below” orientations. Furthermore, the device may be otherwise positioned (e.g., rotated 90 degrees or in other orientations), thus interpreting the spatial relative descriptive terms used herein accordingly.

[0079] The terminology used herein is for the purpose of describing particular embodiments and is not intended to be limiting. As used herein, unless the context clearly indicates otherwise, the singular forms “a” and “the” are intended to include the plural forms as well. Furthermore, when the terms “comprising” and / or “including” and variations thereof are used in this specification, it indicates the presence of the stated features, integrals, steps, operations, parts, components, and / or groups thereof, but does not exclude the presence or addition of one or more other features, integrals, steps, operations, parts, components, and / or groups thereof. It should also be noted that, as used herein, the terms “substantially,” “about,” and other similar terms are used as approximate terms rather than as terms of degree, thus explaining the inherent biases in measurements, calculated values, and / or provided values ​​that would be recognized by one of ordinary skill in the art.

[0080] In one embodiment, such as Figure 1 As shown, this disclosure provides an interpretable spatiotemporal graph anomaly detection method based on graph filtering for anomaly detection in satellite power system telemetry signals, comprising the following steps:

[0081] S100: Preprocesses the raw telemetry data of the satellite power system by nearest neighbor interpolation, moving average filtering, maximum and minimum normalization, and graph structure calculation, and obtains normal samples without abnormal data based on the preprocessed signal for training.

[0082] S200: Construct an interpretable spatiotemporal graph anomaly detection network based on graph filtering;

[0083] S300: The interpretable spatiotemporal graph anomaly detection network is trained using the normal samples. When the number of iterations reaches the set maximum number of iterations or the reduction in the total loss reaches the set threshold, the training of the interpretable spatiotemporal graph anomaly detection network is completed.

[0084] S400: Given any received satellite power system telemetry data, use the trained interpretable spatiotemporal graph anomaly detection network to perform anomaly detection and interpretable analysis.

[0085] The above embodiments constitute the complete inventive concept of this disclosure. It captures spatiotemporal information simultaneously through spatiotemporal graph technology, provides interpretability by utilizing graph filtering technology and ingenious network design, and its autoencoder architecture and adversarial training strategy are robust to high noise and missing data. Therefore, it can achieve high-performance, reliable and practical detection of telemetry signal anomalies in satellite power systems.

[0086] In this embodiment, step S100, the preprocessing includes the following sub-steps:

[0087] S101: The received raw telemetry data of the satellite power system is represented as: {l′1, ..., l′N}, l′ i It is a sequence of length len, representing the observation data received by the i-th sensor, where N is the number of sensors, i ranges from 1 to N, and N is a natural number greater than or equal to 1.

[0088] S102: For {l′1, ..., l′...} N For missing values ​​in the data, i.e., data points marked as NaN or blank cells in the original data, the nearest neighbor interpolation method is used for preliminary data completion, i.e., l′ i The missing value is filled by the preceding or following non-missing value of the missing position, resulting in the padded data: {l″1,…,l″ N};

[0089] S103: For {l″1,...,l″ N Preliminary denoising is performed using a moving average filter to obtain smoothed data {l″′1,...,l″′}. N The calculation formula is:

[0090]

[0091] Where the subscripts i and j represent the i-th sensor and the j-th time step, respectively, a is the set sliding window length, and b is the sliding window length. i The number of sliding windows that the space is divided into.

[0092] S104: For {l″′1…,l″′ N Performing a maximum-minimum normalization operation, the normalized data is represented as {l1, ..., l...} N The calculation formula is:

[0093]

[0094] s105: Using {l1, ..., l N Perform graph structure calculations:

[0095]

[0096] Among them o i Represents the i-th sensor node, o i Representing the j-th sensor node, all sensor nodes form a graph node set o = {o1, ..., o2}. N},(o i o j ) represents node o i and node o j Establish edges between each other, and all edges form the edge set E = {(o1, o2, o3)} of the graph. i ), ..., (oj o N The graph structure is represented as G = (O, E), sim is =d(l i , l j ) represents the observed sequence l i and l j The Euclidean distance between them, TopM represents the sim ij In the current case i, the set of nodes corresponding to the first m maximum values, A represents the adjacency matrix, used to store the edge information of the graph structure, and used to further calculate the Laplacian matrix and the normalized Laplacian matrix:

[0097]

[0098] Where D represents the degree matrix, and its calculation method is: D = diag(d1, ..., d2) N ), d i =∑ j A ij d i =∑ j A ij d i This represents the degree of the i-th sensor, i.e., the number of edges established.

[0099] S106: The normalized data {l1, ..., l...} is processed... N The original input multivariate time series data is divided into multiple subsequences using the sliding window method, with each subsequence representing a normal sample X.

[0100] In this embodiment, step S200, constructing an interpretable spatiotemporal graph anomaly detection network based on graph filtering, includes the following sub-steps:

[0101] S201: Establish an adaptive quadratic approximation graph filter AQAGF(·):

[0102]

[0103] Where ReLU(·) is a nonlinear activation function, K is the order of the empirically set polynomial, and a k b k c k Both z and X are network learning parameters, X is the normal sample matrix, I is the identity matrix, and A and U are the eigenvalue matrix and Fourier transform matrix of the Laplacian matrix, respectively. Their calculation methods are as follows:

[0104] In this step, such as Figure 2As shown, the combination of multiple quadratic functions with positive values ​​effectively approximates the shape of any desired filter frequency response function. Each quadratic term can achieve a general single-peak response, and can even realize classic filters such as low-pass, high-pass, band-pass, and band-stop. They extract or suppress corresponding frequency components respectively. By combining these classic filters, the desired frequency response can be obtained. It is clear that the adaptive quadratic approximation graph filter established above has two advantages: First, it is computationally simple. Compared with similar graph filters that are often based on computationally complex polynomials such as Laplace and Bernstein polynomials, the adaptive quadratic approximation graph filter is based only on easily computed quadratic polynomials and has an advantage in fitting the shape of the filter required in real-world applications. Therefore, it can balance computational simplicity and high fitting accuracy. Second, it can provide interpretability. Compared with the spatial information capture modules of similar graph neural networks, the proposed filter can display the extracted and suppressed frequency components through the frequency response function. Combined with the network design, it can show the location of abnormal and normal frequency bands to maintenance personnel, which helps in further anomaly localization and attribution.

[0105] S202: Establishing an interpretable spatiotemporal graph anomaly detection network:

[0106]

[0107] Where E(·), G(·), and D(·) represent the encoder, generator, and discriminator, respectively. It is a randomly generated high-dimensional space representation vector, belonging to intermediate variables, v, and These are the three discrimination scores output by the interpretable spatiotemporal graph anomaly detection network. It is the reconstructed sample output of an interpretable spatiotemporal graph anomaly detection network. It is by The random reconstructed sample matrix generated by the generator, GRL(·) is the gradient reversal function, which means that during network training, the variables are not changed when calculating the loss, but a negative sign is added before the gradient of the variable when the gradient is backpropagated.

[0108] In this step, such as Figure 3 As shown, sample X is first reconstructed using an encoder and a generator. To guide the training of the discriminator and enhance its generalization ability, a randomized high-dimensional space representation is set. and obtain it through the generator Finally, the discriminator is used to analyze X. Perform discrimination and obtain discrimination scores v and v respectively. and In this process, preferably, in order to accelerate adversarial training, iterative optimization is not adopted, but a gradient reversal layer is introduced to optimize the reconstruction part and the discriminator part in a single gradient backpropagation step, thereby achieving simultaneous optimization of the encoder, generator and discriminator.

[0109] In this embodiment, step S300, the network training includes the following sub-steps:

[0110] S301: Input the normal sample X into the interpretable spatiotemporal graph anomaly detection network to obtain each v,

[0111] S302: Calculate the backpropagation loss value based on the output of the interpretable spatiotemporal graph anomaly detection network. The calculation formula is:

[0112] L=α·L R +L D ,

[0113]

[0114]

[0115] Among them, L, L R and L D These represent the total training loss, reconstruction loss, and discrimination loss, respectively. α and ε are balance coefficients, and MAE(·) represents the mean absolute error.

[0116] In this step, such as Figure 3 As shown, the total training loss L is composed of the reconstruction loss L... R and identification loss L D The composition, where α is the balance coefficient between the two. The reconstruction loss is the sum of the input sample X and the reconstructed sample X. The mean absolute error is used to enable the interpretable spatiotemporal graph anomaly detection network to reconstruct data. The first term of the discrimination loss is used to distinguish the discrimination score between normal and abnormal samples, the second term ensures that the discriminator's discrimination score for abnormal samples is greater than that for normal samples, and ε is a balancing parameter between the two terms. The third and fourth terms are set for the same reason.

[0117] S303: Backpropagate the calculated total training loss L to train the interpretable spatiotemporal graph anomaly detection network;

[0118] S304: Repeat steps S301 to S303. Based on experience, the total training loss L obtained after 20 consecutive iterations of training no longer produces a new minimum value, which means that the training of the spatiotemporal graph anomaly detection network is complete.

[0119] In this embodiment, step S400 includes the following sub-steps:

[0120] S401: Input any received satellite power system telemetry data into the trained interpretable spatiotemporal graph anomaly detection network to obtain the anomaly score matrix, calculated as follows:

[0121]

[0122] Where X is the input sample, These are reconstructed samples, where the subscripts i and t represent the i-th sensor and the t-th time step, respectively. Let S represent the reconstructed anomaly score matrix, where D(X) is the anomaly score matrix obtained by the decoder from the input sample. i，t Here is the anomaly score matrix, and δ is the balance parameter;

[0123] In this step, such as Figure 4 As shown, during the testing process, since no network optimization or updates are required, gradient inversion layers and random high-dimensional space representations are not needed. and its generated samples All other inputs and outputs are the same as during network training. Input the test samples to obtain the reconstructed samples. We then use the identification score v to further calculate the anomaly score matrix S. As you can see, we provide anomaly scores for each sensor at each time stamp, which also supports subsequent anomaly localization and attribution.

[0124] S402: Set the maximum value of the anomaly scores of each sensor at each time step as the final anomaly score determined for that time step. The calculation formula is:

[0125] In this step, we use the maximum value of the anomaly score for each sensor as the final anomaly score for that timestamp, which will be used to output the anomaly timestamp of the final detection result.

[0126] S403: Based on experience, a threshold is set, and those with anomaly scores higher than the threshold at each time point are marked as anomalies. Anomaly detection is completed, and at the same time, the frequency response function obtained by the network is interpreted.

[0127] In this embodiment, in step S202, the encoder, generator, and discriminator respectively include:

[0128]

[0129] In the above formula, MLP(x) is a multilayer perceptron model, expressed as:

[0130] MLP(x)=σ(b MLP +W MLP x)

[0131] Where σ(·) represents the sigmoid nonlinear activation function, b MLP and W MLP These are the bias terms and weight matrices of the multilayer perceptron model.

[0132] GRU(x) is a gated recurrent neural network unit, expressed as:

[0133] r t =σ(W rg ·[h t-1 x t ]+b rg )

[0134] u t =σ(W ug ·[h t-1 x t ]+b ug )

[0135]

[0136]

[0137] Where r represents resetting the gate's output value, u represents updating the gate's output value, h represents the hidden state, and x represents the input value. The latest hidden state is represented by the subscripts t and t-1, which represent the current signal segment and the previous signal segment, respectively. W represents the weight matrix, b represents the bias term, rg represents the reset gate, ug represents the update gate, tanh(·) represents the hyperbolic tangent function, and o represents the Hadamard product.

[0138] In this step, such as Figure 3 As shown, the encoder consists of a spatial information aggregation module and a temporal feature extraction module, used to capture spatial and temporal correlations. The spatial information aggregation module employs the adaptive quadratic approximation graph filter, whose frequency response function can intuitively explain which frequency bands the encoder extracts to reconstruct sample data. In the temporal feature extraction module, a gated recurrent neural network unit (GRU) is used for temporal information extraction, and only a multilayer perceptron (MLP) is used as the generator. The discriminator consists of the adaptive quadratic approximation graph filter and the multilayer perceptron (MLP). The former is used to extract discriminative information, and its frequency response function can intuitively explain which frequency bands the discriminator focuses on for discrimination. The latter, the MLP, maps the extracted features to a score.

[0139] In one specific embodiment, to verify the superior anomaly detection performance of the proposed method, two recognized and effective multidimensional time-series anomaly detection evaluation datasets, Mars Science Laboratory (MSL) and Soil Moisture Active & Passive (SMAP), were used. Comparative experiments were conducted with five classically recognized and effective methods in this field: Deep-SVDD, LSTM-VAE, LSTM-NDT, DAGMM, and OmniAnomaly, and five state-of-the-art methods from the past three years: MERLIN, ITAD, MTAD-GAT, USAD, and MAD-GAN. The experimental results are shown in Table 1. The performance of each method was evaluated using three universally accepted metrics in anomaly detection: precision, recall, and F1 score. Precision is the ratio of correctly predicted positive observations to the total number of predicted positive observations, reflecting the percentage of detected anomalies that are true anomalies. Recall is the ratio of correctly predicted positive observations to all observations in the actual category, reflecting the proportion of correctly identified anomalies among all true anomalies. The F1 score is a weighted average of precision and recall, used to characterize the overall performance of anomaly detection models by combining recall and precision. A missing data experiment was conducted, with 30% of the input data manually removed. The results (second to last row) demonstrate that our model maintains a fairly good anomaly detection capability even with missing data. Furthermore, when evaluating the superiority of a method, the F1 score should be the primary consideration, as it is the harmonic mean of precision and recall, reflecting the comprehensive performance of anomaly detection. Clearly, our method (last two rows) exhibits the best overall anomaly detection performance.

[0140] Table 1

[0141]

[0142] In another specific embodiment, to verify the practicality and interpretability of the proposed method, the method was applied to the BIRDS dataset. Figure 5 This diagram shows representative channels and dataset partitioning of the dataset. It reveals that the telemetry data from the measured satellite power system is extremely noisy. Furthermore, due to the significant imbalance in length among the four subsets, and the fact that only Raavana and Uguisu show anomalies, and given that our training data comes entirely from the Tsuru satellite, inconsistent data distribution is unavoidable and a challenge that intelligent methods must face when applied to real-world tasks. The proposed method overcomes these shortcomings by simultaneously capturing spatiotemporal information and demonstrating robustness to noise, achieving high-precision anomaly detection results with an accuracy of 0.8703, a recall of 0.9835, and an F1 score of 0.9235. Figure 6 The results of anomaly detection for representative channels are shown, demonstrating that the proposed method has excellent data reconstruction and anomaly detection capabilities. Figure 7 This is a schematic diagram illustrating the interpretable results of this specific embodiment, showing the frequency response functions of the encoder spatial information aggregation filter and the discriminator spatial information aggregation filter. Interpretive analysis reveals that: reconstructed information is mainly concentrated in the mid-frequency band around 1.0; discriminative information is mainly concentrated in the low-frequency band from 0 to 0.5 and the high-frequency band from 1.5 to 2; and anomaly information may be distributed in the 0 to 0.5 and 1.5 to 2 frequency bands. These interpretable analyses support the further extraction of temporal features from certain frequency bands. These features, along with the reconstructed data, discriminator feature extraction, and high-dimensional spatial vector representation, can be directly used in subsequent analyses such as anomaly analysis, anomaly localization, and attribution.

[0143] In the description of this specification, the references to terms such as "one embodiment / mode," "some embodiments / modes," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment / mode or example is included in at least one embodiment / mode or example of this application. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment / mode or example. Moreover, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments / modes or examples. Furthermore, without contradiction, those skilled in the art can combine and integrate the different embodiments / modes or examples described in this specification, as well as the features of different embodiments / modes or examples.

[0144] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Thus, a feature defined as "first" or "second" may explicitly or implicitly include at least one of that feature. In the description of this application, "multiple" means at least two, such as two, three, etc., unless otherwise explicitly specified.

[0145] Those skilled in the art should understand that the above embodiments are merely for illustrating the present disclosure and are not intended to limit the scope of the disclosure. Those skilled in the art can make other changes or modifications based on the above disclosure, and these changes or modifications still fall within the scope of the present disclosure.

Claims

1. A method for interpretable spatiotemporal graph anomaly detection based on graph filtering, characterized in that, The method includes the following steps: S100: Preprocess the received raw telemetry data from the satellite power system to obtain normal samples; S200: Construct an interpretable spatiotemporal graph anomaly detection network based on graph filtering; S300: The interpretable spatiotemporal graph anomaly detection network is trained using the normal samples. When the number of iterations reaches the set maximum number of iterations or the reduction in the total loss reaches the set threshold, the training of the interpretable spatiotemporal graph anomaly detection network is completed. S400: Given any received satellite power system telemetry data, use the trained interpretable spatiotemporal graph anomaly detection network to perform anomaly detection and interpretable analysis; Specifically, step S200 includes the following steps: S201: Establish an adaptive quadratic approximation graph filter : ， in, It is a non-linear activation function. It is an empirical setting of the polynomial order. , , and These are all network learning parameters. It is a normal sample matrix. It is the identity matrix. and These are the eigenvalue matrix and the Fourier transform matrix of the Laplace matrix, respectively, and their calculation methods are as follows: ; S202: Establish an interpretable spatiotemporal graph anomaly detection network based on graph filtering; wherein, the interpretable spatiotemporal graph anomaly detection network is represented as: ， in, , , These represent the encoder, generator, and discriminator, respectively. It is a randomly generated high-dimensional space representation vector. , and These are the three discrimination scores output by the interpretable spatiotemporal graph anomaly detection network. It is the reconstructed sample output of an interpretable spatiotemporal graph anomaly detection network. It is by Randomly reconstructed sample matrix generated by the generator It is the gradient reversal function, which means that during network training, the variables are not changed when calculating the loss, but a negative sign is added before the gradient of the variable when the gradient is backpropagated.

2. The method according to claim 1, characterized in that, In step S100, the preprocessing includes the following steps: S101: The raw telemetry data of the satellite power system is expressed as follows: , The length is the length of the observed data. The sequence represents the observation data received by the i-th sensor. This represents the number of sensors, where i ranges from 1 to N, and N is a natural number greater than or equal to 1. S102: Yes Missing values ​​in the data are initially filled using the nearest neighbor imputation method to obtain the filled data. ; S103: Completed data Initial denoising is performed using a moving average filter to obtain smoothed data. ; S104: For the smoothed data Perform a maximum-minimum normalization operation to obtain the normalized data. ; S105: Using normalized data Perform graph structure calculations; S106: The normalized data The original input multivariate time series data is divided into multiple subsequences using the sliding window method, with each subsequence representing a normal sample. .

3. The method according to claim 2, characterized in that, In step S105, the graph structure calculation is represented as follows: ， in, Representing the One sensor node, Representing the There are 1 sensor node, and all sensor nodes form a graph node set. , Represents a node and nodes Edging is established between the edges, and all edges form the edge set of the graph. The graph structure is represented as , Represents the observed sequence and The Euclidean distance between them, TopM represents Given the current i-th case, the set of nodes corresponding to the first m maximum values; This represents the adjacency matrix, used to store edge information of the graph structure, and is used for further calculation of the Laplacian matrix. and normalized Laplace matrix ;in, Represents the adjacency matrix. .

4. The method according to claim 1, characterized in that, Step S300 includes the following steps: S301: Normal sample Input the interpretable spatiotemporal graph anomaly detection network to obtain each ; S302: Calculate the backpropagation loss value according to the output of the interpretable spatiotemporal graph anomaly detection network; S303: Through total training loss Perform backpropagation to train the interpretable spatiotemporal graph anomaly detection network; S304: Repeat steps S301 to S303 until the number of iterations reaches the set maximum number of iterations or the total training loss. When the decrease in the value reaches the set threshold, the training of the interpretable spatiotemporal graph anomaly detection network is complete.

5. The method according to claim 4, characterized in that, In step S302, the loss value is calculated as follows: ， in, , and These represent the total training loss, reconstruction loss, and discrimination loss, respectively. and Represents the balance coefficient. It represents the mean absolute error.

6. The method according to claim 5, characterized in that, Step S400 includes the following steps: S401: Input any received satellite power system telemetry data into the trained interpretable spatiotemporal graph anomaly detection network to obtain an anomaly score matrix; S402: Set the maximum value of the abnormal scores of each sensor at each time point as the final abnormal score determined at that time point; S403: Based on experience, a threshold is set, and those with anomaly scores higher than the threshold at each time point are marked as anomalies. Anomaly detection is completed, and at the same time, the frequency response function obtained by the network is interpreted.

7. The method according to claim 6, characterized in that, In step S401, the formula for calculating the anomaly score matrix is: ， in, This is an abnormal score matrix. For balancing parameters.

8. The method according to claim 7, characterized in that, In step S202, the encoder, generator, and discriminator respectively include: ， In the formula For the multilayer perceptron model, the expression is: ， in, This represents the sigmoid non-linear activation function. and These are the bias terms and weight matrices of the multilayer perceptron model.

9. The method according to claim 8, characterized in that, For a gated recurrent neural network unit, the expression is: ， in, This indicates resetting the output value of the door. This indicates updating the output value of the gate. Indicates a hidden state. Indicates the input value. This indicates the latest hidden state, with subscripts t and t-1 representing the current signal segment and the previous signal segment, respectively. Represents the weight matrix. This indicates the bias term; the subscript 'rg' represents the reset gate, and the subscript 'ug' represents the update gate. denoted by , and ⊙ denotes the hyperbolic tangent function.

Images