Unsupervised rapid radio storm searching method based on time sequence anomaly detection

Abstract

The invention discloses an unsupervised rapid radio storm searching method based on time sequence anomaly detection. The method comprises the following steps: carrying out achromatic processing on observation data of a radio telescope to obtain a signal intensity time sequence; embedding the signal intensity time sequence into a low-dimensional space and constructing a graph structure, and calculating a sequence abnormal score based on the graph structure to obtain a fast radio burst pulse signal candidate body; and screening the fast radio burst pulse signal candidates by using a barrel counting method, and confirming the arrival time of the fast radio burst pulse signal and the dispersion value corresponding to the fast radio burst pulse signal. According to the method, the fast radio storm searching problem is innovatively solved from the perspective of time sequence anomaly detection, the method has the advantages of being high in robustness and low in false detection rate, the number of candidate bodies is greatly reduced, and the searching efficiency is improved; the consumption of computing resources is greatly reduced by using coarse-grained grid achromatic dispersion; the method does not depend on prior knowledge related to signal shapes, and can be expanded to be used for detecting more unknown universe signals.

Description

Technical Field

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

[0002] Fast Radio Bursts (FRB) are a type of ultra-bright radio pulse signals with a duration of milliseconds, and the instantaneous radiation flux can reach tens of Jy. So far, many properties of FRB remain unknown, so research based on a large amount of astronomical observation data is needed to answer questions such as the origin and radiation mechanism of FRB. Currently, more and more astronomical telescopes have begun to search for and observe FRBs, and have produced a series of scientific results. As many astronomical telescopes observe FRBs, the technical bottleneck of searching for FRBs in massive observation data has gradually emerged: FRBs have one-time and transient characteristics, and real-time searches for FRBs have high requirements for accuracy and speed. Therefore, accurate and efficient FRB search technology is crucial to achieve selective real-time voltage data dumping.

[0003] At present, these single pulse events are mainly searched by professional software based on the theory of achromatic dispersion, such as Presto software and HEIMDALL software. Professional software will output a large number of candidates during the search process. Due to factors such as noise, radio frequency interference (RFI), and system gain changes, the screening of candidates requires professionals with domain knowledge. However, with the continuous increase in the amount of observational data, the traditional manual screening method has become unsustainable and urgently needs to be optimized. In recent years, with the development of computer technology, machine learning algorithms have been widely used in pattern recognition, signal classification and other application directions, and have performed well. Supervised machine learning algorithms have been successfully applied in the field of radio astronomy, such as the identification and classification of pulsar candidates. Supervised learning requires labeled data as training samples. However, the number of FRB sources discovered so far is extremely rare, and the number of FRB pulse signals found is also limited; secondly, labeling data is time-consuming, labor-intensive, and costly.

[0004] The present invention transforms the search for FRB pulse signals into anomaly detection, and establishes an unsupervised fast radio burst search method based on time series anomaly detection, which is used to solve the accuracy and efficiency bottlenecks faced in the current FRB search work. Summary of the invention

[0005] The purpose of the present invention is to provide an unsupervised fast radio burst search method based on time series anomaly detection to solve the problem that the fast radio burst search method in the prior art has low accuracy and efficiency and relies on manual screening of a large number of candidates. The present invention does not rely on prior knowledge of the shape of cosmic signals, can be expanded to detect more unknown cosmic signals, and can be used for the search of fast radio bursts and pulsar signals.

[0006] The objective of the present 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) Preprocess the original observation data of the radio telescope to eliminate the influence of radio frequency interference and noise, and then use an adaptive strategy to determine the N dispersion values ​​of the coarse-grained grid de-dispersion. Normalize the spectrum data after de-dispersion in the frequency dimension to obtain N signal intensity time series. Perform average pooling and maximum pooling on the signal intensity time series to retain different features, and generate 2N signal intensity time series data to be detected;

[0008] (2) For the 2N signal intensity time series data obtained in step (1), use a method based on time series anomaly detection to search for possible anomalies therein and obtain fast radio burst pulse signal candidates;

[0009] (3) Use the bucket counting method to screen all the 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, the step (1) specifically includes:

[0011] Preprocess the original observation data of the radio telescope to eliminate the influence of radio frequency interference and noise; for the value range of the dispersion value [DM left ,DM right ], an adaptive strategy is used to set the dispersion value step size to determine the N dispersion values ​​of the coarse-grained grid dedispersion, and the spectrum data is corrected through the dedispersion operation, and the spectrum data is normalized in the frequency dimension to obtain the signal intensity time series data after dedispersion of different dispersion values; further, each signal intensity time series data is subjected to average pooling and maximum pooling operations to retain different features, and 2N signal intensity time series data to be detected are obtained.

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

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

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

[0015] (2.3) Calculate the sequence anomaly score: q (l q ≥l) is mapped to a path in the graph, the score corresponding to the subsequence is calculated according to the degree of the vertices on the path and the weight of the edge, and then transformed to obtain the anomaly score corresponding to the subsequence;

[0016] (2.4) Sequence anomaly score enhancement: For the same dispersion value and different pooling methods, the anomaly scores calculated for the corresponding time series are fused and enhanced by weighted average; anomaly The anomaly score is subjected to non-maximum suppression, the maximum value within the fixed window is retained, and the non-maximum value is filtered. Each retained anomaly score represents a sequence anomaly, that is, a candidate for a fast radio burst pulse signal, and the corresponding time is the time when the anomaly occurs.

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

[0018] (2.1.1) One-dimensional sliding convolution denoising: For each signal strength time series T to be detected, first use a sliding window of length l to extract subsequences, sliding one sequence point each time to divide T into multiple continuous subsequences of length l; then apply a one-dimensional convolution operation with a convolution kernel size of λ to each subsequence to obtain a vector representation of length l-λ;

[0019] (2.1.2) Use PCA method to reduce dimension: combine the vector representations of all subsequences obtained in step (2.1.1) into a matrix Proj(T,l,λ)∈R (|T|-l)×(l-λ) , and use the PCA method to reduce its dimensionality, retain the first three principal components, and obtain the reduced-dimensional matrix Proj r(T,l,λ)∈R (|T|-l)×3 , this operation is recorded as PCA3;

[0020] (2.1.3) Spatial rotation: Rotate the reduced-dimensional matrix Proj r (T, l, λ), so that the unit vector of the three-dimensional space after dimensionality reduction The vector with the direction of change of the time series value Aligned, then there are two remaining unit vectors The corresponding unit vector after rotation The resulting projected two-dimensional space contains features related to the shape of the sequence.

[0021] Furthermore, the time series value change direction vector is obtained by the following method:

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

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

[0024] (2.2.1) Extract vertices: Take the sampling rate as r, that is, divide the circumference into r equal parts to obtain the sampling angle set For each sampling angle ψ∈Ψ, calculate its radius vector After any two adjacent subsequences are transformed by step (2.1), two consecutive points x in the projected two-dimensional space i ,x i+1 The intersection points between the line segments to obtain the intersection point set For each intersection set I ψ , use nonparametric estimation methods to calculate the probability of each element, and select the intersection point with the largest probability As a graph vertex, to get the graph vertex set

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

[0026] (2.2.3) Constructing a graph structure: Constructing a graph structure based on 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 relationships between patterns, and the edge weights represent the frequency of transitions between patterns.

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

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

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

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

[0031] All the abnormal occurrence times under different dispersion value settings obtained in step (2.4) are divided into bucket widths t bin Divide the data into buckets, count the frequency of each bucket, calculate the corresponding frequency as the confidence level, and select the bucket with a confidence level greater than the threshold value θ. conf The corresponding time is the final anomaly occurrence time, that is, 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.

[0032] The beneficial effects of the present invention are that the present invention innovatively solves the search problem of fast radio bursts from the perspective of time series anomaly detection, has the characteristics of strong robustness and low false detection rate, greatly reduces the number of fast radio burst pulse signal candidates, and improves search efficiency; converts the search for cosmic pulse signals into sequence anomaly detection, does not rely on prior knowledge related to signal shape, and can be expanded to detect more unknown cosmic signals; the search process uses coarse-grained grid dispersion elimination, which greatly reduces the consumption of computing resources compared to traditional fine-grained grid dispersion elimination; at the same time, this method belongs to the category of unsupervised learning, does not require labeled data, and avoids the problem of labeled data in the astronomical field. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0034] Figure 2 It is a schematic diagram of the processing of coarse-grained grid de-dispersion in the present invention;

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

[0036] Figure 4 It is a spatial schematic diagram of the subsequence after embedding and a calculated sequence anomaly score result diagram of the present invention;

[0037] Figure 5 It is a schematic diagram of sequence abnormality enhancement and bucket counting post-processing in the present invention. DETAILED DESCRIPTION

[0038] Here, exemplary embodiments will be described in detail, and examples thereof are shown in the accompanying drawings. When the following description refers to the drawings, unless otherwise indicated, the same numbers in different drawings represent the same or similar elements. Obviously, the drawings used in the following description are only some embodiments of the present invention, and for ordinary technicians in this field, other drawings can also be obtained based on these drawings without paying creative labor. The implementation methods described in the following exemplary embodiments do not represent all implementation methods consistent with the present invention. On the contrary, they are only examples of devices and methods consistent with some aspects of the present invention as detailed in the attached claims.

[0039] The terms used in the present invention are only for the purpose of describing specific embodiments and are not intended to limit the present invention. The singular forms "a", "the" and "the" used in the present invention and the appended claims are also intended to include plural forms unless the context clearly indicates other meanings. It should also be understood that the term "and / or" used herein refers to and includes any or all possible combinations of one or more associated listed items.

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

[0041] The present invention is described in detail below in conjunction with the accompanying drawings. In the absence of conflict, the features of the following embodiments and implementations can be combined with each other.

[0042] The idea of ​​implementing the unsupervised fast radio burst search method based on time series anomaly detection described in the present invention is: de-dispersion processing is performed on the time-frequency data, the spectrum data is corrected, and the spectrum data is normalized in the frequency dimension to obtain signal intensity time series data; the signal intensity time series is embedded in a low-dimensional space; a graph structure (Graph) is constructed based on the data in the low-dimensional space; the sequence to be detected is mapped to a path (Path) in the graph, and the anomaly score is calculated; the calculated anomaly score is processed using non-maximum suppression; the anomaly score obtained under different dispersion value (Dispersion Measure, DM) settings is bucketed and counted to determine the time when the anomaly occurs and the corresponding DM value.

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

[0044] (1) The original observation data of the radio telescope are preprocessed to eliminate the influence of radio frequency interference and noise. Then, an adaptive strategy is used to determine the N dispersion DM values ​​of the coarse-grained grid. 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 respectively subjected to average pooling and maximum pooling processing to retain different features, generating 2N signal intensity time series data to be detected.

[0045] Specifically, Figure 2 As shown in the figure, the original observation data of the radio telescope is preprocessed to eliminate the influence of radio frequency interference and noise; for the DM value interval [DM left ,DM right], an adaptive strategy is used to set the DM step size to determine the N DM values ​​of the coarse-grained grid de-dispersion, and the spectrum data is corrected by the de-dispersion operation, and the spectrum data is normalized in the frequency dimension to obtain the signal intensity time series data after de-dispersion of different DM parameters; further, each signal intensity time series data is subjected to average pooling and maximum pooling operations, such as Figure 3 As shown, 2N signal intensity time series data to be detected are obtained.

[0046] It should be understood that using the coarse-grained grid de-dispersion method to de-disperse the collected time-frequency data can solve the problem of inconsistent arrival times of high-frequency signals and low-frequency signals caused by dispersion. At the same time, compared with the traditional fine-grained grid de-dispersion, it can greatly reduce 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 therein and obtain fast radio burst (FRB) pulse signal candidates.

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

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

[0050] (2.1.1) One-dimensional sliding convolution denoising processing: For each signal strength time series T to be detected, first use a sliding window of length l to extract subsequences, sliding one sequence point each time to divide T into multiple continuous subsequences of length l; then apply a one-dimensional convolution operation with a convolution kernel size of λ to each subsequence to reduce the impact of noise and highlight the shape characteristics of the sequence, and obtain a vector representation of length l-λ.

[0051] (2.1.2) Use PCA method to reduce dimension: combine the vector representations of all subsequences obtained in step (2.1.1) into a matrix Proj(T,l,λ)∈R (|T|-l)×(l-λ) , and use the PCA method to reduce its dimensionality, retaining the first three-dimensional principal components, which capture the most important change trend in the sequence, and obtain the matrix Proj after dimensionality reduction. r (T,l,λ)∈R (|T|-l)×3 , and this operation is denoted as PCA3.

[0052] (2.1.3) Spatial rotation: Rotate the reduced-dimensional matrix Proj r (T, l, λ), so that the unit vector of the three-dimensional space after dimensionality reduction The vector with the direction of change of the time series value Aligned, then there are two remaining unit vectors The corresponding unit vector after rotation The constructed projected two-dimensional space contains features related to the shape of the sequence. After transformation, the point distances of two sequences with similar shapes in the projected two-dimensional space are close.

[0053] Furthermore, the time series value change direction vector It is obtained by the following method: for the signal strength time series T to be detected, the minimum and maximum values ​​are multiplied by the convolution kernel size λ and then multiplied by the unit vector of length l-λ, and then PCA3 dimension reduction is performed to obtain the minimum point and the maximum point, and the vector composed of the minimum point and the maximum point is used as the direction vector of the time series value change It can be expressed 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, max(T) represents the maximum value in the time series T of the signal strength to be detected, and 1 l-λ represents 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. The graph vertices represent the sequence patterns, the edges represent the transition relationship between the patterns, and the edge weights represent the frequency of transitions between the patterns.

[0059] It should be understood that a sequence pattern refers to sequence fragments with similar shapes.

[0060] (2.2.1) Extract vertices: Take the sampling rate as r, that is, divide the circumference into r equal parts to obtain the sampling angle set For each sampling angle ψ∈Ψ, calculate its radius vector After any two adjacent subsequences are transformed by step (2.1), two consecutive points x in the projected two-dimensional space i ,x i+1 The intersection points between the line segments to obtain the intersection point set For each intersection set I ψ , use nonparametric estimation methods to calculate the probability of each element, and take the intersection point with the largest probability As a graph vertex, to get the graph vertex set

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

[0062]

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

[0064]

[0065] Among them, the subscripts 1 and 2 represent the intersection set I and ψ Element x i The first and second dimensions of , σ1 and μ2 represent the mean of the first and second dimensions respectively.

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

[0067] (2.2.3) Constructing the graph structure: Construct the graph structure based on the edge set obtained in step (2.2.2) and the vertex set obtained in step (2.2.1), such as Figure 4 As shown in the figure, 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.

[0068] (2.3) Calculate the sequence anomaly score: q (l q ≥l) is mapped to a path in the graph, and the score corresponding to the subsequence is calculated according to the degree of the vertices on the path and the weight of the edge, and then transformed to obtain the anomaly score corresponding to the subsequence.

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

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

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

[0072] (2.4) Enhancement of sequence anomaly scores: For the same DM value and different pooling methods, the anomaly scores calculated for the corresponding time series are fused and enhanced by weighted average to improve robustness; anomaly The anomaly score is suppressed for non-maximum values, the maximum value in the fixed window is retained, and the non-maximum value is filtered. Each retained anomaly score represents a sequence anomaly, that is, a candidate for the FRB pulse signal, and the corresponding time is the time when the anomaly occurs, such as Figure 5 shown.

[0073] It should be understood that the same anomaly may appear in multiple adjacent columns of length l. q In the subsequence of , the scores of adjacent points in the abnormal subsequence will be 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 occurs.

[0074] Specifically, in order to reduce the impact of noise data on the calculation results of anomaly scores and improve the robustness of the results, the anomaly scores calculated for the same DM value and different pooling methods (i.e., average pooling operation and maximum pooling operation) are fused and enhanced by weighted average, which is conducive to improving the robustness of the results. When the window slides, the same anomaly will appear in multiple windows of length l. q Therefore, in order to determine the time when the anomaly occurs, based on the fusion enhancement, the non-maximum suppression method is used to screen the anomaly points in the fixed window, which are not suppressed and are greater than the preset anomaly threshold θ anomaly The anomaly score represents a sequence anomaly, which is a candidate for the FRB pulse signal, and the corresponding time is the time when the anomaly occurred.

[0075] It should be understood that the processing flow of the non-maximum suppression method is: Input: abnormal score sequence, abnormal score threshold θ anomaly ; Output: The final retained anomaly subset B; Step 1: The median of the anomaly sequence is greater than the threshold θ anomaly 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 positions are within the radius r of this point and the point itself from set A (for one-dimensional sequences, delete the points whose original sequence positions are within the radius r to the left and right of this 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, Figure 5As shown in the figure, all abnormal occurrence times under different DM value settings obtained in step (2.4) are divided into bucket widths t bin Divide the data into buckets, count the frequency of each bucket, calculate the corresponding frequency as the confidence level, and select the bucket with a confidence level greater than the threshold value θ. conf The corresponding time is the final anomaly occurrence time, that is, the arrival time ToA of the FRB pulse signal, and the mean 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 anomaly score in the anomaly score sequence corresponds to a length of l q , associate the starting time point of the sequence, calculate the bucket number of the time point and put it into the corresponding bucket. bin =50ms, t=21.625s is counted into the 432nd bucket.

[0079] In summary, the present invention innovatively solves the problem of fast radio burst search from the perspective of time series anomaly detection. It has the characteristics of strong robustness and low false detection rate, and avoids the problem of traditional methods relying on manual screening of a large number of candidates; this method does not rely on prior knowledge of the shape of cosmic signals and can be expanded to detect more unknown cosmic signals; at the same time, this method uses coarse-grained grid dispersion elimination, which greatly reduces the consumption of computing resources compared to traditional fine-grained grid dispersion elimination; the present invention belongs to the category of unsupervised learning and does not require labeled data, avoiding the problem of labeled data in the astronomical field.

[0080] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some of the technical features may be replaced by equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions 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 following steps are involved: (1) Preprocess the original observation data of the radio telescope to eliminate the influence of radio frequency interference and noise, and then use an adaptive strategy to determine the N dispersion values ​​of the coarse-grained grid de-dispersion. Normalize the spectrum data after de-dispersion in the frequency dimension to obtain N signal intensity time series. Perform average pooling and maximum pooling on the signal intensity time series to retain different features, and generate 2N signal intensity time series data to be detected; (2) For the 2N signal intensity time series data obtained in step (1), use a method based on time series anomaly detection to search for possible anomalies therein and obtain fast radio burst pulse signal candidates; (3) Use the bucket counting method to screen all the 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.

2. The unsupervised fast radio burst search method based on time series anomaly detection according to claim 1 is characterized in that: The step (1) specifically comprises: Preprocess the original observation data of the radio telescope to eliminate the influence of radio frequency interference and noise; for the value range of the dispersion value [DM left ,DM right ], an adaptive strategy is used to set the dispersion value step size to determine the N dispersion values ​​of the coarse-grained grid dedispersion, and the spectrum data is corrected through the dedispersion operation, and the spectrum data is normalized in the frequency dimension to obtain the signal intensity time series data after dedispersion of different dispersion values; further, each signal intensity time series data is subjected to average pooling and maximum pooling operations 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 is characterized in that: The step (2) specifically includes the following sub-steps: (2.1) Construct subsequence embedding representation: First, use a sliding window of length l to divide the signal intensity time series T to be detected into multiple continuous subsequences, and apply local convolution operation to each subsequence to obtain its vector representation; then use the principal component analysis method to reduce the dimension of the matrix composed of the vector representation of the subsequence, retain the first three-dimensional principal components, and obtain the reduced-dimensional matrix; then rotate the reduced-dimensional matrix so that a unit vector in the reduced-dimensional three-dimensional space is aligned with the direction vector of the time series value change, and the two-dimensional space formed by the corresponding unit vectors after the rotation of the remaining two unit vectors contains features related to the sequence shape; (2.2) Constructing a graph structure: Based on the subsequence embedding representation constructed in step (2.1), extract the vertices and edges of the graph to construct a graph structure; where the graph vertices represent sequence patterns, the edges represent the transition relationship between patterns, and the edge weights represent the frequency of transitions between patterns; (2.3) Calculate the sequence anomaly score: q (l q ≥l) is mapped to a path in the graph, the score corresponding to the subsequence is calculated according to the degree of the vertices on the path and the weight of the edge, and then transformed to obtain the anomaly score corresponding to the subsequence; (2.4) Sequence anomaly score enhancement: For the same dispersion value and different pooling methods, the anomaly scores calculated for the corresponding time series are fused and enhanced by weighted average; anomaly The anomaly score is subjected to non-maximum suppression, the maximum value within the fixed window is retained, and the non-maximum value is filtered. Each retained anomaly score represents a sequence anomaly, that is, a candidate for a fast radio burst pulse signal, and the corresponding time is the time when the anomaly occurs.

4. The unsupervised fast radio burst search method based on time series anomaly detection according to claim 3 is characterized in that: The step (2.1) specifically includes the following sub-steps: (2.1.1) One-dimensional sliding convolution denoising: For each signal strength time series T to be detected, first use a sliding window of length l to extract subsequences, sliding one sequence point each time to divide T into multiple continuous subsequences of length l; then apply a one-dimensional convolution operation with a convolution kernel size of λ to each subsequence to obtain a vector representation of length l-λ; (2.1.2) Use PCA method to reduce dimension: combine the vector representations of all subsequences obtained in step (2.1.1) into a matrix Proj(T,l,λ)∈R (|T|-l)×(l-λ) , and use the PCA method to reduce its dimensionality, retain the first three principal components, and obtain the reduced-dimensional matrix Proj r (T,l,λ)∈R (|T|-l)×3 , this operation is recorded as PCA3; (2.1.3) Spatial rotation: Rotate the reduced-dimensional matrix Proj r (T, l, λ), so that the unit vector of the three-dimensional space after dimensionality reduction The vector with the direction of change of the time series value Aligned, then there are two remaining unit vectors The corresponding unit vector after rotation The resulting projected two-dimensional space contains features related to the shape of the sequence.

5. The unsupervised fast radio burst search method based on time series anomaly detection according to claim 3 is characterized in that: The time series value change direction vector is obtained by the following method: For the signal intensity time series T, the minimum and maximum values ​​are multiplied by the convolution kernel size λ and then by the unit vector of length l-λ. The minimum and maximum points are obtained after PCA3 dimensionality reduction. The vector composed of the minimum and maximum points is used as the direction vector of the time series value change.

6. The unsupervised fast radio burst search method based on time series anomaly detection according to claim 3 is characterized in that: The step (2.2) specifically includes the following sub-steps: (2.2.1) Extract vertices: Take the sampling rate as r, that is, divide the circumference into r equal parts to obtain the sampling angle set For each sampling angle ψ∈Ψ, calculate its radius vector After any two adjacent subsequences are transformed by step (2.1), two consecutive points x in the projected two-dimensional space i ,x i+1 The intersection points between the line segments to obtain the intersection point set For each intersection set I ψ , use nonparametric estimation methods to calculate the probability of each element, and select the intersection point with the largest probability As a graph vertex, to get the graph vertex set (2.2.2) Extract edges: For any two adjacent subsequences, the points x in the projected two-dimensional space i ,x i+1 , select the vertex with the closest Euclidean distance (N (i) ,N (i+1) ) is put into the edge set ε, and the weight of the edge is set to the edge corresponding to the two adjacent subsequences during the window sliding process (N (i) ,N (i+1) ) appears, that is, from the graph vertex N during the window sliding process (i) The corresponding sequence pattern is transferred to N (i+1) The number of corresponding sequence patterns; (2.2.3) Constructing a graph structure: Constructing a graph structure based on 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 relationships between patterns, and the edge weights represent the frequency of transitions between patterns.

7. The unsupervised fast radio burst search method based on time series anomaly detection according to claim 3 is characterized in that: The step (2.3) specifically includes the following sub-steps: (2.3.1) Path mapping: For any length l in the signal strength time series T to be detected q (l q ≥l) subsequence After step (2.1), the subsequence embedding representation in the projected two-dimensional space is obtained; then after the edge extraction operation in step (2.2.2), the subsequence Paths in a graph (2.3.2) Anomaly score calculation: According to path P th The degree of the upper vertex and the weight of the edge calculate the score corresponding to the subsequence Where w represents the weight of the edge, deg(N (j) ) represents vertex N (j) Calculate the length of all the l in T q The score of the subsequence is normalized by MinMax, 1.0-MinMax_Norm(Score(P th )) is the abnormal score of the corresponding subsequence, where MinMax_Norm(·) represents MinMax normalization.

8. The unsupervised fast radio burst search method based on time series anomaly detection according to claim 1, characterized in that: The step (3) specifically comprises: All the abnormal occurrence times under different dispersion value settings obtained in step (2.4) are divided into bucket widths t bin Divide the data into buckets, count the frequency of each bucket, calculate the corresponding frequency as the confidence level, and select the bucket with a confidence level greater than the threshold value θ. conf The corresponding time is the final anomaly occurrence time, that is, 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