An unsupervised fast radio burst search method based on time series anomaly detection

Abstract

The application discloses a kind of based on time series anomaly detection's unsupervised fast radio burst searching method, this method includes: to radio telescope observation data is dispersed and is handled to obtain signal intensity time series;Signal intensity time series is embedded into low-dimensional space and constructs graph structure, based on graph structure calculation sequence anomaly divides to obtain fast radio burst pulse signal candidate body;Using bucket counting method to the fast radio burst pulse signal candidate body is screened, confirms the arrival time of fast radio burst pulse signal and its corresponding dispersion value.The application innovatively solves the searching problem of fast radio burst from the angle of time series anomaly detection, with the characteristics of strong robustness, low false alarm rate, greatly reduces the number of candidate body, improves searching efficiency;Using coarse-grained grid dispersion, greatly reduce the consumption of computing resources;Not dependent on signal shape related prior knowledge, can be extended to detect more unknown cosmic signals.

Description

Technical Field

[0001] This invention relates to the field of radio astronomy, and in particular to an unsupervised fast radio burst search method based on time series anomaly detection. Background Technology

[0002] Fast radio bursts (FRBs) are a class of ultra-bright radio pulses with durations on the order of milliseconds and instantaneous radiative fluxes reaching tens of jansky (Jy). To date, many properties of FRBs remain unknown, necessitating research based on vast amounts of astronomical observation data to answer questions about their origin and radiation mechanisms. Currently, an increasing number of astronomical telescopes are searching for and observing FRBs, yielding a series of scientific results. However, as numerous telescopes observe FRBs, the technical bottlenecks in searching for them within massive amounts of observational data are becoming increasingly apparent: FRBs are one-off and transient, and real-time FRB searching demands high precision and speed. Therefore, accurate and efficient FRB search techniques are crucial for achieving selective real-time voltage data dumping.

[0003] Currently, these single-pulse events are primarily searched using specialized software based on achromatic theory, such as Presto and Heimdall. This software generates a large number of candidates during the search process. Due to factors such as noise, radio frequency interference (RFI), and system gain variations, screening these candidates requires professionals with domain knowledge. However, with the continuous increase in observational data, traditional manual screening methods are becoming unsustainable and urgently need optimization. In recent years, with the development of computer technology, machine learning algorithms have been widely applied and performed excellently in pattern recognition, signal classification, and other applications. Supervised machine learning algorithms have been successfully applied in radio astronomy, such as the identification and classification of pulsar candidates. Supervised learning requires labeled data as training samples; however, the number of discovered FRB sources is extremely small, and the number of searched FRB pulse signals is also limited. Furthermore, labeling data is time-consuming, labor-intensive, and costly.

[0004] This invention transforms FRB pulse signal search into anomaly detection, and establishes an unsupervised fast radio burst search method based on time series anomaly detection to solve the accuracy and efficiency bottlenecks currently faced in FRB search. Summary of the Invention

[0005] The purpose of this invention is to address the problems of low accuracy and efficiency in existing fast radio burst (FRB) search methods, which rely on manual screening of a large number of candidate samples. This invention provides an unsupervised FRB search method based on time-series anomaly detection. This invention does not rely on prior knowledge of the shape of cosmic signals and can be extended to detect more unknown cosmic signals, applicable to the search for FRBs and pulsar signals.

[0006] The objective of this invention is achieved through the following technical solution: an unsupervised fast radio burst search method based on time series anomaly detection, comprising the following steps:

[0007] (1) The raw observation data of the radio telescope is preprocessed to eliminate radio frequency interference and noise. Then, an adaptive strategy is used to determine the N dispersion values ​​of the coarse-grained grid dedispersion. The spectrum data after dedispersion processing is normalized in the frequency dimension to obtain N signal intensity time series. The signal intensity time series are subjected to average pooling and maximum pooling respectively to retain different features, generating 2N signal intensity time series data to be detected.

[0008] (2) For the 2N signal intensity time series data obtained in step (1), the time series anomaly detection method is used to search for possible anomalies and obtain fast radio burst pulse signal candidates.

[0009] (3) Use the bucket counting method to screen all fast radio burst pulse signal candidates obtained in step (2) and confirm the arrival time of the fast radio burst pulse signal and its corresponding dispersion value.

[0010] Furthermore, step (1) specifically includes:

[0011] The raw observation data from the radio telescope is preprocessed to eliminate radio frequency interference and noise; for the range of dispersion values ​​[DM] left DM right An adaptive strategy is used to set the dispersion value step size to determine N dispersion values ​​for coarse-grained grid de-dispersion. The spectrum data is then corrected through de-dispersion operation, and the spectrum data is normalized in the frequency dimension to obtain signal intensity time series data after de-dispersion with different dispersion values. Furthermore, average pooling and maximum pooling operations are performed on each signal intensity time series data to retain different features, resulting in 2N signal intensity time series data to be detected.

[0012] Furthermore, step (2) specifically includes the following sub-steps:

[0013] (2.1) Constructing the subsequence embedding representation: First, a sliding window of length l is used to divide the signal intensity time series T to be detected into multiple continuous subsequences, and a local convolution operation is applied to each subsequence to obtain its vector representation; then, the matrix composed of the vector representations of the subsequences is reduced in dimension by using principal component analysis, retaining the first three principal components to obtain the dimension-reduced matrix; then, the dimension-reduced matrix is ​​rotated so that a unit vector in the dimension-reduced three-dimensional space is aligned with the vector of the direction of change of the time series value, and the two-dimensional space formed by the unit vectors corresponding to the remaining two unit vectors after rotation contains features related to the shape of the sequence;

[0014] (2.2) Constructing the graph structure: Based on the subsequence embedding representation constructed in step (2.1), the vertices and edges of the graph are extracted to construct the graph structure; where the graph vertices represent the sequence patterns, the edges represent the transition relationships between patterns, and the edge weights represent the frequency of transitions between patterns.

[0015] (2.3) Calculate the sequence outlier: Take the sequence of length l from T... q (l q Subsequences ≥ l) are mapped to paths in the graph. The scores corresponding to the subsequences are calculated based on the degree of the vertices and the weights of the edges on the paths, and then transformed to obtain the abnormal scores corresponding to the subsequences.

[0016] (2.4) Sequence anomaly enhancement: For time series with the same dispersion value but different pooling methods, the anomalies are enhanced by weighted averaging; for those greater than the threshold θ... anomaly Non-maximum suppression is performed on the outlier segments, retaining the maximum values ​​within a fixed window and filtering out the non-maximum values. Each retained outlier segment represents a sequence outlier, which is a candidate for the fast radio burst pulse signal, and its corresponding time is the time of the outlier occurrence.

[0017] Furthermore, step (2.1) specifically includes the following sub-steps:

[0018] (2.1.1) One-dimensional sliding convolution noise reduction: For each signal intensity time series T to be detected, a sliding window of length l is first used to extract subsequences. Each time, one sequence point is slid to divide T into multiple continuous subsequences of length l. Then, a one-dimensional convolution operation with a kernel size of λ is applied to each subsequence to obtain a vector representation of length l-λ.

[0019] (2.1.2) Dimensionality reduction using PCA: Combine the vector representations of all subsequences obtained in step (2.1.1) into a matrix Proj(T,l,λ)∈R in order. (|T|-l)×(l-λ) The dimensionality of the matrix Proj was reduced using PCA, retaining the first three principal components. r(T,l,λ)∈R (|T|-l)×3 This operation is denoted as PCA3;

[0020] (2.1.3) Spatial rotation: The matrix Proj after rotation and dimensionality reduction r (T,l,λ), so that the unit vector in the reduced three-dimensional space Vector of change direction of time series values Alignment will leave two unit vectors. The corresponding unit vector after rotation The constructed two-dimensional projection space contains features related to the shape of the sequence.

[0021] Furthermore, the vector representing the direction of change of the time series values ​​is obtained using the following method:

[0022] For a signal strength time series T, the minimum and maximum values ​​are multiplied by the convolution kernel size λ, then multiplied by a unit vector of length l-λ. After PCA3 dimensionality reduction, the minimum and maximum points are obtained. The vector formed by the minimum and maximum points is used as the direction vector of the time series value change.

[0023] Furthermore, step (2.2) specifically includes the following sub-steps:

[0024] (2.2.1) Vertex extraction: Take a sampling rate of r, that is, divide the circumference into r equal parts to obtain the set of sampling angles. For each sampling angle ψ∈Ψ, calculate its radius vector. After transformation in step (2.1), any two adjacent subsequences are represented by two consecutive points x in the projected two-dimensional space. i ,x i+1 The intersection points between the line segments are used to obtain the set of intersection points. For each set of intersection points I ψ The probability of each element is calculated using a nonparametric estimation method, and the intersection point with the highest probability is selected. As graph vertices, to obtain the graph vertex set.

[0025] (2.2.2) Edge extraction: For any two adjacent subsequences, the points x in the projected two-dimensional space... i ,x i+1 Choose the vertex with the closest Euclidean distance. (N) (i) N (i+1) Add the edge to the edge set ε, and set the weight of the edge to the weight of the edges corresponding to two adjacent subsequences during the window sliding process (N). (i) N (i+1) The number of times ) appears, that is, from vertex N in the graph during the window sliding process. (i)The corresponding sequence pattern is transferred to N. (i+1) The number of times the corresponding sequence pattern is used;

[0026] (2.2.3) Constructing the graph structure: Construct the graph structure based on the set of edges obtained in step (2.2.2) and the set of graph vertices obtained in step (2.2.1); where the graph vertices represent sequence patterns, the edges represent the transition relationships between patterns, and the edge weights represent the frequency of transitions between patterns.

[0027] Furthermore, step (2.3) specifically includes the following sub-steps:

[0028] (2.3.1) Path mapping: For any length l in the time series T of the signal strength to be detected... q (l q ≥l) subsequence Step (2.1) yields the subsequence embedding representation in the projected two-dimensional space; then, the edge extraction operation in step (2.2.2) yields the subsequence. Paths in a graph structure

[0029] (2.3.2) Anomaly score calculation: based on path P th The score is calculated by using the degree of the upper vertex and the weight of the edge to determine the subsequence. Where w represents the edge weight, deg(N) (j) ) represents vertex N (j) The degree; calculate all lengths l in T. q The scores of the subsequences are calculated and then normalized using MinMax, 1.0 - MinMax_Norm(Score(P)). th )) is the abnormal score of the corresponding subsequence, where MinMax_Norm(·) represents MinMax normalization.

[0030] Furthermore, step (3) specifically includes:

[0031] All anomaly occurrence times under different dispersion values ​​obtained in step (2.4) are calculated according to the bucket width t. bin Divide the data into buckets, count the frequency of each bucket, calculate the corresponding frequency as the confidence level, and select buckets with a confidence level greater than a threshold θ. conf The bucket is defined by the time of the final anomaly, which is the arrival time of the fast radio burst pulse signal. The average of all dispersion values ​​falling into the bucket is calculated as the dispersion value corresponding to the fast radio burst pulse signal.

[0032] The beneficial effects of this invention are as follows: It innovatively solves the search problem for fast radio bursts from the perspective of time-series anomaly detection, exhibiting strong robustness and a low false detection rate, significantly reducing the number of fast radio burst pulse signal candidates and improving search efficiency; it transforms the search for cosmic pulse signals into sequence anomaly detection, independent of prior knowledge related to signal shape, and can be extended to detect more unknown cosmic signals; the search process employs coarse-grained grid dispersion reduction, which significantly reduces computational resource consumption compared to traditional fine-grained grid dispersion reduction; simultaneously, this method falls under the category of unsupervised learning, requiring no labeled data, thus avoiding the data labeling difficulties in the astronomical field. Attached Figure Description

[0033] Figure 1 This is a flowchart of the unsupervised fast radio burst search method based on time series anomaly detection of the present invention;

[0034] Figure 2 This is a schematic diagram of the coarse-grained mesh dedispersion process in this invention;

[0035] Figure 3 This is a schematic diagram of the pooling operation in this invention;

[0036] Figure 4 This is a spatial schematic diagram of the subsequence embedding in this invention and a diagram of the calculated sequence anomaly score.

[0037] Figure 5 This is a schematic diagram of the sequence anomaly enhancement and bin counting post-processing in this invention. Detailed Implementation

[0038] Exemplary embodiments will now be described in detail, examples of which are illustrated in the accompanying drawings. In the following description, when referring to the drawings, unless otherwise indicated, the same numbers in different drawings represent the same or similar elements. It is obvious that the drawings used in the following description are merely some embodiments of the invention, and those skilled in the art can obtain other drawings based on these drawings without any inventive effort. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with the present invention. Rather, they are merely examples of apparatuses and methods consistent with some aspects of the invention as detailed in the appended claims.

[0039] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. The singular forms “a,” “the,” and “the” used in this invention and the appended claims are also intended to include the plural forms unless the context clearly indicates otherwise. It should also be understood that the term “and / or” as used herein refers to and includes any or all possible combinations of one or more of the associated listed items.

[0040] It should be understood that although the terms first, second, third, etc., may be used in this invention to describe various information, this information should not be limited to these terms. These terms are only used to distinguish information of the same type from one another. For example, first information may also be referred to as second information without departing from the scope of this invention, and similarly, second information may also be referred to as first information. Depending on the context, the word "if" as used herein may be interpreted as "when," "when," or "in response to a determination."

[0041] The present invention will now be described in detail with reference to the accompanying drawings. Unless otherwise specified, the features of the following embodiments and implementations can be combined with each other.

[0042] The approach to implementing the unsupervised fast radio burst search method based on time series anomaly detection described in this invention is as follows: Dispersion-reducing processing is performed on time-frequency data to correct the spectral data; the spectral data is normalized in the frequency dimension to obtain signal intensity time-series data; the signal intensity time-series data is embedded into a low-dimensional space; a graph structure is constructed based on the data in the low-dimensional space; the sequence to be detected is mapped to a path in the graph, and anomaly scores are calculated; non-maximum suppression is used to process the calculated anomaly scores; and the anomaly scores obtained under different dispersion measure (DM) settings are binned and counted to determine the time of anomaly occurrence and the corresponding DM value.

[0043] See Figure 1 The unsupervised fast radio burst search method based on time series anomaly detection of the present invention specifically includes the following steps:

[0044] (1) The original observation data of the radio telescope is preprocessed to eliminate radio frequency interference and noise. Then, an adaptive strategy is used to determine the N dispersion DM values ​​of the coarse-grained grid de-dispersion. The spectrum data after de-dispersion processing is normalized in the frequency dimension to obtain N signal intensity time series. The signal intensity time series are subjected to average pooling and maximum pooling respectively to retain different features, generating 2N signal intensity time series data to be detected.

[0045] Specifically, such as Figure 2 As shown, the raw observation data from the radio telescope is preprocessed to eliminate radio frequency interference and noise; for the DM value range [DM... left DM rightAn adaptive strategy is employed to set the DM step size to determine N DM values ​​for coarse-grained grid dedispersion. The spectral data is then corrected through dedispersion operations, and normalization is performed on the spectral data in the frequency dimension to obtain signal intensity time series data after dedispersion with different DM parameters. Furthermore, average pooling and maximum pooling operations are performed on each signal intensity time series data, such as... Figure 3 As shown, 2N time series data of signal intensity to be detected are obtained.

[0046] It should be understood that using coarse-grained grid de-dispersion methods to process the acquired time-frequency data can solve the problem of inconsistent arrival times of high-frequency and low-frequency signals caused by dispersion. At the same time, compared with traditional fine-grained grid de-dispersion, it significantly reduces the consumption of computing resources.

[0047] (2) For the 2N signal intensity time series data to be detected obtained in step (1), an unsupervised time series anomaly detection method is used to search for possible anomalies and obtain Fast Radio Burst (FRB) pulse signal candidates.

[0048] (2.1) Constructing the subsequence embedding representation: First, a sliding window of length l is used to divide the signal intensity time series T to be detected into multiple continuous subsequences, and a local convolution operation is applied to each subsequence to reduce the influence of noise and highlight the shape features of the sequence, so as to obtain its vector representation; then, the matrix composed of the vector representations of the subsequences is reduced in dimension using the principal component analysis (PCA) method, and the first three principal components are retained to obtain the dimension-reduced matrix; then the dimension-reduced matrix is ​​rotated so that a unit vector in the dimension-reduced three-dimensional space is aligned with the vector of the direction of change of the time series value, and the two-dimensional space formed by the unit vectors corresponding to the remaining two unit vectors after rotation contains the features related to the shape of the sequence.

[0049] It should be understood that after the dimension-reduced matrix is ​​rotated, one dimension contains the signal intensity information, and the other two dimensions contain the signal shape information. The two-dimensional space formed by these two dimensions contains features related to the sequence shape.

[0050] (2.1.1) One-dimensional sliding convolution noise reduction: For each signal intensity time series T to be detected, a sliding window of length l is first used to extract subsequences. Each time, one sequence point is slid to divide T into multiple continuous subsequences of length l. Then, a one-dimensional convolution operation with a kernel size of λ is applied to each subsequence to reduce the influence of noise and highlight the shape features of the sequence, resulting in a vector representation of length l-λ.

[0051] (2.1.2) Dimensionality reduction using PCA: Combine the vector representations of all subsequences obtained in step (2.1.1) into a matrix Proj(T,l,λ)∈R in order. (|T|-l)×(l-λ) The dimensionality of the sequence was reduced using PCA, retaining the first three principal components, which capture the most significant trends in the sequence, resulting in the dimensionality-reduced matrix Proj. r (T,l,λ)∈R (|T|-l)×3 This operation is denoted as PCA3.

[0052] (2.1.3) Spatial rotation: The matrix Proj after rotation and dimensionality reduction r (T,l,λ), so that the unit vector in the reduced three-dimensional space Vector of change direction of time series values Alignment will leave two unit vectors. The corresponding unit vector after rotation The constructed two-dimensional projection space contains features related to the shape of the sequence. After transformation, the points of two sequences with similar shapes are close in the two-dimensional projection space.

[0053] Furthermore, the vector of direction of change of time series values The minimum and maximum values ​​are obtained as follows: For the signal intensity time series T to be detected, the minimum and maximum values ​​are multiplied by the convolution kernel size λ, then multiplied by a unit vector of length l-λ. After PCA3 dimensionality reduction, the minimum and maximum points are obtained. The vector formed by the minimum and maximum points is used as the direction vector of the time series value change. It can be represented as:

[0054]

[0055] O min =PCA3(min(T)*λ*1 l-λ )

[0056] O max =PCA3(max(T)*λ*1) l-λ )

[0057] Where, min(T) represents the minimum value in the time series T of the signal strength to be detected, and max(T) represents the maximum value in the time series T of the signal strength to be detected. l-λ Let l represent a unit vector of length l-λ.

[0058] (2.2) Constructing the graph structure: Based on the subsequence embedding representation constructed in step (2.1), the vertices and edges of the graph are extracted to construct the graph structure. Here, the graph vertices represent sequence patterns, the edges represent the transition relationships between patterns, and the edge weights represent the frequency of transitions between patterns.

[0059] It should be understood that sequence patterns refer to sequence segments with similar shapes.

[0060] (2.2.1) Vertex extraction: Take a sampling rate of r, that is, divide the circumference into r equal parts to obtain the set of sampling angles. For each sampling angle ψ∈Ψ, calculate its radius vector. After transformation in step (2.1), any two adjacent subsequences are represented by two consecutive points x in the projected two-dimensional space. i ,x i+1 The intersection points between the line segments are used to obtain the set of intersection points. For each set of intersection points I ψ The probability of each element is calculated using a nonparametric estimation method, and the intersection point with the highest probability is selected. As graph vertices, to obtain the graph vertex set.

[0061] Furthermore, the probability density function f h (x,I ψ The expression for ) is:

[0062]

[0063] Among them, f h (x,I ψ ) represents the probability density function, and n is a set I. ψ The size of h is the bandwidth parameter, and h is selected here. saott In the form of μ(I) ψ ) and σ(I ψ ) represent sets I respectively ψ The mean and standard deviation of all points, where the formula for calculating the standard deviation is:

[0064]

[0065] Wherein, subscripts 1 and 2 represent the set of intersection points I, respectively. ψ element x i The first and second dimensions, σ1 and μ2 represent the mean values ​​of the first and second dimensions, respectively.

[0066] (2.2.2) Edge extraction: For any two adjacent subsequences, the points x in the projected two-dimensional space... i ,x i+1 Choose the vertex with the closest Euclidean distance. (N) (i) N (i+1) Add the edge to the edge set ε, and set the weight of the edge to the weight of the edges corresponding to two adjacent subsequences during the window sliding process (N). (i) N(i+1) The number of times ) appears, that is, from vertex N in the graph during the window sliding process. (i) The corresponding sequence pattern is transferred to N. (i+1) The number of times the corresponding sequence pattern is used.

[0067] (2.2.3) Constructing the graph structure: Based on the edge set obtained in step (2.2.2) and the vertex set obtained in step (2.2.1), construct the graph structure, as follows: Figure 4 As shown in the figure; where the vertices of the graph represent sequence patterns, the edges represent the transition relationships between patterns, and the edge weights represent the frequency of transitions between patterns.

[0068] (2.3) Calculate the sequence outlier: Take the sequence of length l from T... q (l q Subsequences ≥ l) are mapped to paths in the graph. The scores corresponding to the subsequences are calculated based on the degree of the vertices and the weights of the edges on the paths, and then transformed to obtain the abnormal scores corresponding to the subsequences.

[0069] It should be understood that, in T, the length is l q (l q The subsequences of ≥l) are mapped to the paths in the graph, and then their outlier scores are calculated according to the formula in step (2.3.2) below.

[0070] (2.3.1) Path mapping: For any length l in the time series T of the signal strength to be detected... q (l q ≥l) subsequence Step (2.1) yields the subsequence embedding representation in the projected two-dimensional space; then, the edge extraction operation in step (2.2.2) yields the subsequence. Paths in a graph structure

[0071] (2.3.2) Anomaly score calculation: based on path P th The score is calculated by using the degree of the upper vertex and the weight of the edge to determine the subsequence. Where w represents the edge weight, deg(N) (j) ) represents vertex N (j) The degree; calculate all lengths l in T. q The subsequence score is calculated and then normalized using MinMax, 1.0 - MinMax_Norm(Score(P)). th This represents the anomaly score of the corresponding subsequence. The anomaly score result is as follows: Figure 4 As shown, MinMax_Norm(·) represents MinMax normalization.

[0072] (2.4) Sequence outlier enhancement: For time series with the same DM value but different pooling methods, outlier scores are fused and enhanced using a weighted average to improve robustness; for scores greater than the threshold θ... anomaly The outlier components are subjected to non-maximum suppression, retaining the maxima within a fixed window and filtering out the non-maximums. Each retained outlier component represents a sequence anomaly, i.e., an FRB pulse signal candidate, and its corresponding time is the anomaly occurrence time. Figure 5 As shown.

[0073] It should be understood that the same anomaly may appear in multiple adjacent lengths of l. q In the subsequence, the scores of adjacent points in the abnormal subsequence are very high. It is necessary to use the non-maximum suppression method to determine which points in the abnormal subsequence to retain. The time corresponding to the retained points is the time when the anomaly occurred.

[0074] Specifically, to reduce the impact of noisy data on the outlier calculation results and improve the robustness of the results, outlier scores calculated from sequences with the same DM value but different pooling methods (i.e., average pooling and max pooling) are fused and enhanced using a weighted average. This improves the robustness of the results. When the window slides, the same outlier may appear in multiple sequences of length l. q Therefore, in order to determine the time of anomaly occurrence, based on the fusion enhancement, a non-maximum suppression method is used to filter out anomaly segments within a fixed window that are not suppressed and are greater than a preset anomaly segment threshold θ. anomaly The abnormal segment represents a sequence abnormality, namely the FRB pulse signal candidate, and the corresponding time is the time of the abnormality occurrence.

[0075] It should be understood that the processing flow of the nonmaximum suppression method is as follows: Input: outlier sequence, outlier threshold θ anomaly Output: The final set of retained outlier segments B; Step 1: The outlier segment sequence with values ​​greater than the threshold θ anomaly Step 1: Add the point with the largest current value to set A; Step 2: Select the point with the largest current value from set A and add it to set B. At the same time, delete the points in set A whose original sequence position is within the radius r of the point and the point itself (for a one-dimensional sequence, this means deleting the points whose original sequence position is within the radius r to the left and right of the point); Step 3: Repeat Step 2 until set A is empty.

[0076] (3) Use the bucket counting method to screen the FRB pulse signal candidates obtained in step (2) and confirm the arrival time (ToA) of the FRB pulse signal and its corresponding DM value.

[0077] Specifically, such as Figure 5As shown, the occurrence times of all anomalies under different DM value settings obtained in step (2.4) are calculated according to the bucket width t. bin Divide the data into buckets, count the frequency of each bucket, calculate the corresponding frequency as the confidence level, and select buckets with a confidence level greater than a threshold θ. conf The bucket is defined by the time of the final anomaly, which is the arrival time of the FRB pulse signal, ToA. The average value of all dispersion values ​​DM falling into the bucket is calculated as the dispersion value corresponding to the FRB pulse signal.

[0078] It should be noted that each outlier in the outlier sequence corresponds to a sequence of length l. q For each subsequence, associate it with its starting time point, calculate the bucket number for that time point, and place it into the corresponding bucket. For example, t bin =50ms, t=21.625s, then it is counted in the 432nd bucket.

[0079] In summary, this invention innovatively addresses the fast radio burst (FRB) search problem from the perspective of time series anomaly detection, exhibiting strong robustness and a low false positive rate, thus avoiding the problem of traditional methods relying on manual screening of a large number of candidate bursts. This method does not depend on prior knowledge of the shape of cosmic signals, allowing it to be extended to detect more unknown cosmic signals. Furthermore, this method uses coarse-grained grid achromatic dispersion reduction, significantly reducing computational resource consumption compared to traditional fine-grained grid achromatic dispersion reduction. This invention falls under the category of unsupervised learning, requiring no labeled data, thus avoiding the data labeling challenges in the astronomical field.

[0080] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. An unsupervised fast radio burst search method based on time series anomaly detection, characterized in that, The method comprises the following steps: (1) preprocessing the original observation data of the radio telescope, eliminating the influence of radio frequency interference and noise, then determining N dispersion values of coarse-grained grid dispersion by using an adaptive strategy, normalizing the dispersion-processed spectral data in the frequency dimension to obtain N signal intensity time series, and performing average value pooling and maximum value pooling processing on the signal intensity time series respectively to retain different characteristics, and generating 2N signal intensity time series data to be detected; (2) searching for possible abnormalities in the 2N signal intensity time series data obtained in step (1) by using a time series anomaly detection-based method to obtain fast radio burst pulse signal candidates; the step (2) specifically comprises the following sub-steps: (2.1) Constructing the subsequence embedding representation: First, use a length of The sliding window will display the time series of the signal strength to be detected. The sequence is divided into multiple continuous subsequences, and a local convolution operation is applied to each subsequence to obtain its vector representation. Then, principal component analysis is used to reduce the dimension of the matrix composed of the vector representations of the subsequences, retaining the first three principal components to obtain the dimension-reduced matrix. The dimension-reduced matrix is ​​then rotated so that a unit vector in the dimension-reduced three-dimensional space is aligned with the vector of the direction of change of the time series value. The two-dimensional space formed by the unit vectors corresponding to the remaining two unit vectors after rotation contains features related to the shape of the sequence. (2.2) constructing a graph structure: on the basis of the sub-sequence embedding representation constructed in step (2.1), extracting graph vertices and edges to construct a graph structure; wherein the graph vertices represent sequence patterns, the edges represent the transition relationship between patterns, and the edge weights represent the frequency of pattern transition; (2.3) Calculate sequence anomaly score: map the subsequence with length of the graph, calculate the score corresponding to the subsequence according to the degree of the vertex and the weight of the edge on the path, and convert it to obtain the anomaly score corresponding to the subsequence; wherein, ;​ (2.4) Sequence anomaly score enhancement: For the same dispersion value, the anomaly scores calculated by different pooling methods are fused by weighted average to enhance the sequence anomaly score. For the anomaly scores greater than the threshold , non-maximum suppression is performed to retain the maximum value in the fixed window and filter the non-maximum value. Each retained anomaly score represents a sequence anomaly, i.e., a fast radio burst pulse signal candidate, and the corresponding time is the anomaly occurrence time. (3) screening all fast radio burst pulse signal candidates obtained in step (2) by using a bucket counting method, and confirming the arrival time of the fast radio burst pulse signal and the corresponding dispersion value.

2. The unsupervised fast radio burst search method based on time series anomaly detection according to claim 1, wherein, The step (1) specifically comprises: The raw observation data of the radio telescope is preprocessed to eliminate the influence of radio frequency interference and noise; for the value interval [DM left , DM right ] of the dispersion value, an adaptive strategy is used to set the dispersion value step to determine N dispersion values of the coarse-grained grid dispersion, and the spectral data is corrected through the dispersion operation, and the spectral data is normalized in the frequency dimension to obtain the signal intensity time series data after dispersion at different dispersion values; further, the average value pooling and maximum value pooling operations are performed on each signal intensity time series data to retain different features, and 2N signal intensity time series data to be detected are obtained.

3. The unsupervised fast radio burst search method based on time series anomaly detection according to claim 1, wherein, The step (2.1) specifically comprises the following sub-steps: (2.1.1) One-dimensional sliding convolution denoising processing: for each signal intensity time sequence to be detected , first use a sliding window with a length of to extract a subsequence, slide one sequence point at a time, and divide into multiple consecutive subsequences with a length of ; then apply a one-dimensional convolution operation with a convolution kernel size of to each subsequence to obtain a vector representation with a length of ; (2.1.2) Dimension reduction using PCA method: the vector representation of all sub-sequences obtained in step (2.1.1) are combined into a matrix , and dimension reduction is performed using PCA method, retaining the first three principal components, obtaining the reduced matrix , which is denoted as ; (2.1.3) Spatial rotation: rotate the reduced dimension matrix , so that the unit vectors of the three-dimensional space after reduction align with the time series value change direction vectors , then the remaining two unit vectors corresponding to the unit vectors after rotation constitute a projection two-dimensional space containing features related to the sequence shape.

4. The unsupervised fast radio burst search method based on time series anomaly detection according to claim 1, wherein, The time series value change direction vector is obtained by the following method: For signal strength time series The minimum and maximum values are multiplied by the convolution kernel size respectively, and then multiplied by a unit vector of length , and then reduced dimensionally to obtain minimum and maximum points, and a vector composed of the minimum and maximum points is taken as a time series value change direction vector .

5. The unsupervised fast radio burst search method based on time series anomaly detection according to claim 3, wherein, The step (2.2) specifically comprises the following sub-steps: (2.2.1) Extracting vertices: take the sampling rate as , i.e. divide the circumference into equal parts, and get the sampling angle set ; for each sampling angle , calculate its radius vector ; and the intersection point between the line segment formed by the two consecutive points in the projection two-dimensional space after the transformation of any two adjacent subsequences by step (2.1), to obtain the intersection point set ; for each intersection point set , calculate the probability of each element using a non-parametric estimation method, and select the intersection point with the maximum probability as the graph vertex, to obtain the graph vertex set ; wherein represents the probability density function; (2.2.2) Edge extraction: For any two adjacent subsequences, the points in the projected two-dimensional space Choose the vertex with the closest Euclidean distance. ,Will Add to edge set The weight of this edge is set as the weight of the edges corresponding to two adjacent subsequences during the window sliding process. The number of times it appears, that is, from the vertex of the graph during the window sliding process. The corresponding sequence pattern is transferred to The number of times the corresponding sequence pattern is used; (2.2.3) constructing a graph structure: constructing a graph structure according to the edge set obtained in step (2.2.2) and the graph vertex set obtained in step (2.2.1); wherein the graph vertices represent sequence patterns, the edges represent the transition relationship between patterns, and the edge weights represent the frequency of pattern transition.

6. The unsupervised fast radio burst search method based on time series anomaly detection according to claim 1, wherein, The step (2.3) specifically comprises the following sub-steps: (2.3.1) Path mapping: for any sub-sequence of arbitrary length of the time series of signal strength to be detected , its sub-sequence embedding representation in the projected two-dimensional space is obtained by step (2.1), where, ; and then the edge extraction operation in step (2.2.2) is performed to obtain the path of the sub-sequence in the graph structure ; (2.3.2) Anomaly score calculation: based on the path The score is calculated by using the degree of the upper vertex and the weight of the edge to determine the subsequence. Where w represents the weight of the edge, Represents vertices Degree; Calculation All lengths in The scores of the subsequences are calculated and then normalized using MinMax. That is, the abnormal score of the corresponding subsequence, where This indicates MinMax normalization.

7. The unsupervised fast radio burst search method based on time series anomaly detection according to claim 1, wherein, The step (3) specifically comprises: All the abnormal occurrence times of different dispersion values obtained in step (2.4) are binned by the bucket width The frequency of each bucket is counted, and the corresponding frequency is calculated as the confidence. Buckets with a confidence greater than the threshold value are selected The corresponding time of the bucket is the final abnormal occurrence time, i.e. the arrival time of the fast radio burst pulse signal. The mean of all dispersion values falling into the bucket is calculated as the dispersion value corresponding to the fast radio burst pulse signal.

Images