A radio astronomy interference suppression method based on image segmentation

Abstract

This invention belongs to the field of radio astronomy technology, specifically relating to a radio astronomy interference suppression method based on image segmentation. The specific process is as follows: the OTSU algorithm is used to binarize the original time-frequency data to obtain initial mask data; adaptive filtering of the initial mask data along the time and frequency axes is performed to obtain density mask data; the number of retained pixels is counted for each frequency channel of the density mask data, and frequency channels with values ​​greater than a set threshold are marked as RFI frequency channels; change point detection is performed on the RFI frequency channels to determine potential RFI time periods; for the mean sequence formed by the mean values ​​within each RFI time period, abnormal mean sequences are marked; for each abnormal mean sequence, its mean is subtracted to achieve RFI removal.

Description

Technical Field

[0001] This invention belongs to the field of radio astronomy technology, specifically relating to a radio astronomy interference suppression method based on image segmentation. Background Technology

[0002] In the field of radio astronomy, RFI (Radio Frequency Interference) refers to any radio frequency signal that affects astronomical observations. The intensity of RFI is usually at least an order of magnitude higher than the weak radiation from noisy celestial objects, thus overwhelming astronomical signals and severely affecting the quality of observational data.

[0003] Currently, astronomers have proposed several effective RFI suppression techniques. SumThreshold, based on a thresholding method, utilizes the intensity difference between RFI and astronomical signals to detect RFI pixels with significantly higher amplitudes than their neighbors. It exhibits excellent accuracy and computational efficiency, but is prone to missing detections in large interference regions. The Hough line-finding algorithm, based on image processing methods, searches based on the straight-line characteristics of interference, achieving ideal detection results for RFIs of specific shapes, but lacks general applicability. Relevance interference suppression algorithms based on deep learning leverage the powerful nonlinear fitting capabilities of network models to accurately and efficiently detect RFIs; however, the training cost is high.

[0004] For interference detection, current algorithms still have limitations in terms of detection accuracy, computational efficiency, and universality. For interference removal, broadband interference is treated with zero-dispersion filtering, which involves subtracting the mean value of each time channel before dispersion removal. This method suffers from RFI leakage and artifacts in non-RFI regions. Narrowband interference is directly removed from pixels marked as RFI during the interference detection stage, which risks losing astronomical signals. Improving the accuracy and computational efficiency of interference detection, increasing robustness to complex interference environments, and reducing the risk of losing astronomical signals during interference removal are the core objectives of designing radio astronomy interference suppression methods. Summary of the Invention

[0005] To address the problems existing in current radio astronomy interference suppression methods, the present invention aims to provide a radio astronomy interference suppression method based on image segmentation. This method can efficiently and accurately detect RFI in various complex environments, remove RFI while reducing the loss of astronomical signals, improve the accuracy and efficiency of interference detection, and reduce the loss of astronomical signals.

[0006] The technical solution for implementing the present invention is as follows:

[0007] A radio astronomy interference suppression method based on image segmentation includes the following steps: S1: Use the OTSU algorithm to binarize the original time-frequency data to obtain the initial mask data. ; S2: Initial mask data Adaptive filtering is performed along the time and frequency axes. Specifically, this adaptive filtering is implemented as follows: sliding windows are constructed along both the time and frequency axes; the density of the center pixel of each window is calculated; a density threshold is calculated based on the statistical distribution of all pixels in the input data; pixels below this threshold are set to zero; and finally, the union of the filtering results along the time and frequency axes is taken to obtain the density mask data. ; S3: Density mask data The number of retained pixels is counted for each frequency channel, and frequency channels with a value greater than a set threshold are marked as RFI frequency channels; S4: Perform change point detection on the RFI frequency channel to determine potential RFI time periods; S5: For the mean sequence formed by the mean values ​​within each RFI time period, mark the abnormal mean sequence; S6: For each abnormal mean sequence, subtract its mean to remove RFI.

[0008] Optionally, the adaptive filtering in step S2 of the present invention is specifically implemented as follows: setting the size to... The window slides along the time axis with the initial mask data. Convolution, set size to The window slides along the frequency axis with the initial mask data. Convolution calculates the horizontal density of the center pixel of the window. and vertical density According to horizontal density and vertical density Calculate the adaptive filtering thresholds for the time axis and the frequency axis; when the density of the pixel at the center of the window is less than the corresponding adaptive filtering threshold, set that pixel to 0; take the union of the two adaptive density comparison results as the final density mask. .

[0009] Optionally, the present invention provides a horizontal adaptive density threshold. and vertical adaptive density threshold for:

[0010]

[0011] in, This represents the total number of pixels in the input data. Indicates the first The horizontal density of the center pixel of each window Indicates the first Vertical density of the center pixel of each window.

[0012] Optionally, in step S3 of the present invention, the noise threshold is defined as follows: the number of pixels retained in each frequency channel after statistical density filtering, the median and the absolute deviation of the median of the number sequence are calculated, and the sum of the median and N times the absolute deviation of the median is defined as the noise threshold.

[0013] Optionally, the noise threshold of the present invention for:

[0014] in, Represents the median of the sequence. Represents the absolute deviation of the median of the sequence. Optionally, in step S4, the present invention performs the FPOP function pruning optimal segmentation algorithm on the RFI frequency channel to detect change points and determine potential RFI time periods.

[0015] Optionally, in step S4 of the present invention, the FPOP algorithm is used to perform change point detection on each RFI channel sequence to obtain the mean within the segment; wherein, when the FPOP algorithm iteratively solves the problem, it uses the squared error loss function as the fitting cost and sets the penalty cost according to the number of time sampling points of the input data and the statistical variance of the difference time series.

[0016] Optionally, the penalty cost described in this invention for:

[0017] in, This represents the total number of time sampling points for the input data. This represents the standard deviation of the background noise.

[0018] Optionally, in step S5, the present invention performs the SumThreshold algorithm on different mean sequences of the RFI frequency channel to mark abnormal mean sequences.

[0019] Optionally, in step S6 of the present invention, only the segmented mean of the time series marked as abnormal is subtracted, and no processing is performed on the time series not marked as abnormal, thereby preserving the background information.

[0020] Beneficial effects: First, it exhibits strong robustness. This invention combines the characteristics of high intensity and high density of RFI pixels, introduces an OTSU adaptive threshold and designs a bidirectional density filter in steps S1 and S2, which can effectively cope with complex and ever-changing RFI environments. It can not only detect RFI of different intensities, but also capture RFI of various shapes.

[0021] Secondly, it offers precise and efficient positioning. This invention employs a refined, layered detection strategy of "pixel coarse screening – channel positioning – time segmentation – intensity determination." In particular, it utilizes the FPOP algorithm for rapid detection of change points, combined with the SumThreshold algorithm for analyzing the mean sequence. This leverages both the sensitivity of change point detection and the ability of the combined threshold algorithm to identify outliers, significantly improving the time-frequency positioning accuracy and efficiency of RFI.

[0022] Third, it preserves astronomical signals to the greatest extent possible. Traditional RFI processing often employs simple pruning strategies, which can easily lead to data loss. In step S6, this invention unifies the removal methods for narrowband and wideband RFI. By subtracting the adaptive segmented mean, it can completely remove RFI while preserving the background noise superimposed on the interference and the effective astronomical signals to the greatest extent possible, thus reducing the loss of astronomical observation data. Attached Figure Description

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

[0024] Figure 1 This is a flowchart illustrating the implementation of the present invention; Figure 2 This is a schematic diagram of the input data in the form of a time-frequency graph received by the present invention; Figure 3 The image shows the result of binarizing the original data using the OTSU algorithm; Figure 4 This is a diagram showing the effect of density filtering output; Figure 5 A diagram illustrating the effect of marking RFI frequency channels according to noise thresholds; Figure 6 The effect of performing FPOP change point detection on the 1644.5MHz frequency channel and using SumThreshold to mark the actual RFI time period; Figure 7 This is a diagram showing the final clean data output after subtracting the mean from the RFI segment. Detailed Implementation

[0025] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0026] It should be noted that, in the absence of conflict, the following embodiments and features can be combined with each other; and, based on the embodiments of this disclosure, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this disclosure.

[0027] It should be noted that various aspects of embodiments within the scope of the appended claims are described below. It will be apparent that the aspects described herein can be embodied in a wide variety of forms, and any particular structure and / or function described herein is merely illustrative. Based on this disclosure, those skilled in the art will understand that one aspect described herein can be implemented independently of any other aspect, and two or more of these aspects can be combined in various ways. For example, any number of aspects set forth herein can be used to implement the device and / or practice the method. Additionally, this device and / or method can be implemented using structures and / or functionalities other than one or more of the aspects set forth herein.

[0028] The data in this embodiment comes from observations of pulsar PSR B0611+22 by the Lovell 76m radio telescope. The data format is high-resolution time-frequency data, such as... Figure 2 As shown, the horizontal axis represents time, the vertical axis represents frequency, and the pixel value represents signal strength. The observation data from the Lovell telescope contains both narrowband continuous RFI of various shapes and transient broadband pulse RFI, which provides a highly challenging verification platform for testing the suppression effect of the method in dynamic and complex RFI scenarios.

[0029] like Figure 1 As shown in the figure, an embodiment of the present invention discloses a radio astronomy interference suppression method based on image segmentation, comprising: S1: The original time-frequency data is binarized using the OTSU threshold segmentation algorithm. The RFI pixels are coarsely extracted by calculating an adaptive segmentation threshold to obtain the initial mask data. ; The OTSU algorithm used in this invention is an adaptive threshold determination method based on maximizing inter-class variance, widely used in image binarization segmentation. The input to OTSU is... Figure 2 The original time-frequency data shown is used to calculate the optimal segmentation threshold, which is denoted as . , all greater than Pixels with certain noise levels are set to 1, otherwise to 0. In noise-dominated electromagnetic interference environments, OTSU's segmentation threshold tends to fall within the high-energy noise range. Therefore, this step can completely preserve the RFI, but it will also leave behind high-intensity noise points, such as... Figure 3 As shown. The initial mask obtained by OTSU segmentation. In the process, the originally dense noise becomes sparse, while the RFI remains dense.

[0030] S2: Initial mask data Density filtering is performed to remove noise points smaller than the density threshold, resulting in density mask data. ; Initial mask data The characteristic is that noisy pixels are scattered throughout the time-frequency map, while RFI pixels are clustered in local areas. Therefore, two density filter sliding windows with a value of 1 are defined: one with a size of... The window slides along the timeline and Convolution yields the number of narrowband RFI pixels within the window. Then the density of the center pixel of the window is Similarly, the size is The window slides along the frequency axis and Convolution yields the number of broadband RFI pixels within the window. The density of the center pixel of the window is .

[0031] To filter out sparse noise, the horizontal and vertical adaptive density thresholds are designed as follows:

[0032]

[0033] in, This represents the total number of pixels in the input data. Each pixel has a horizontal density. and vertical density When the density of the pixel at the center of the window is less than the corresponding adaptive threshold, that pixel is set to 0. The union of the outputs of the two density filters is used as the final density mask. .like Figure 4 As shown, after density filtering, most of the sparse noise is removed, while the dense RFI clusters are retained.

[0034] S3: Density mask data The number of retained pixels is counted for each frequency channel, a noise threshold is defined, and frequency channels with noise values ​​greater than the noise threshold are marked as RFI frequency channels. Density mask data is obtained in step S2. Afterwards, although most of the sparse noise has been filtered out, it is still necessary to further accurately identify which frequency channels are truly contaminated by RFI. This embodiment transforms the RFI frequency channel detection problem into a one-dimensional sequence outlier detection problem by projecting a two-dimensional time-frequency mask onto the frequency axis.

[0035] Specifically, for density mask data For each frequency channel, the number of pixels with a value of 1 on the time axis is accumulated to obtain the pixel count sequence for each frequency channel, denoted as . .exist In the sequence, points with larger values ​​correspond to channels heavily contaminated by RFI, while points with smaller values ​​correspond to relatively clean background channels. Because... The sequence may contain extreme RFI outliers. Traditional statistical methods using the "mean + standard deviation" approach are easily inflated by these large values, leading to an artificially high detection threshold and thus missing some weak interference. Therefore, this embodiment uses a more robust statistical measure to define the noise threshold. :

[0036] in, The median of the sequence reflects the average level of background noise. The median absolute deviation (MAD) of the sequence is used to characterize the range of background noise fluctuations. Choosing 3 times the MAD as a safety margin effectively covers the statistical fluctuations of most background noise. Finally, the sequence is traversed. All values ​​greater than The frequency channel is labeled as the RFI frequency channel. For example... Figure 5 As shown, channels that exceed the noise threshold are marked, and subsequent steps will only process these marked channels, thus significantly reducing the computational load.

[0037] S4: Execute the FPOP algorithm on the RFI frequency channel to detect change points and determine potential RFI time periods; After identifying the contaminated RFI frequency channels in step S3, it is necessary to further pinpoint the start and end times of the RFI in the time dimension. For each RFI frequency channel, extract its intensity data sequence over time. ,in The sequence length is given.

[0038] This embodiment uses the FPOP algorithm to detect change points in the sequence. FPOP is an exact search algorithm based on dynamic programming. Compared with the traditional binary partitioning method or PELT algorithm, FPOP introduces an inequality pruning mechanism, which can significantly reduce computational complexity while ensuring that the global optimum is found. It is very suitable for processing massive amounts of data in radio astronomy observations.

[0039] Specifically, for each RFI channel sequence, the FPOP algorithm aims to find The optimal variable point This minimizes the following objective function:

[0040] in, The fitting cost within the segment, Represents the first point in the sequence. Includes from to the turning point The observed values ​​(i.e., intensity data). The penalty cost required to introduce each variable.

[0041] In this embodiment, considering the statistical characteristics of radio signals, the specific configuration is as follows: Cost function: The squared error loss is used as the fitting cost within the piecewise segments, i.e. ,in, Let be the mean within the segment. This assumes that the mean of the signal remains relatively stable between adjacent change points.

[0042] Cost of punishment: Cost of punishment The settings determine the detection sensitivity. If If it's too small, it will detect a large number of false changes; if If the value is too large, short-term RFIs will be missed. To adapt to different observation durations and noise levels, this embodiment uses an adaptive penalty term setting:

[0043] in, This represents the total number of time sampling points for the input data. Let be the standard deviation of the background noise. To avoid the influence of RFI itself on the noise estimation, a robust estimation method is used here to calculate . That is, it is estimated by calculating the median absolute deviation of the first difference of the original time series.

[0044] By executing the FPOP algorithm described above, the system outputs a series of change point locations. These change points divide the entire time series into several sub-segments, each with a constant mean. These time intervals defined by the change points are identified as potential RFI time periods, awaiting further numerical determination in step S5. For the effectiveness of FPOP change point detection, please refer to [link to relevant documentation]. Figure 6 .

[0045] S5: Perform the SumThreshold algorithm on different mean sequences of the RFI frequency channel to mark outlier mean sequences; After the FPOP change point detection in step S4, each frequency channel marked as RFI is divided into several consecutive time segments on the time axis. Although FPOP can accurately capture the moment when the signal statistical characteristics change abruptly, it does not directly determine whether a certain segment is RFI or normal background fluctuation. Therefore, the SumThreshold algorithm is needed to perform secondary discrimination based on the mean of these segments.

[0046] In specific implementation, the first step is to extract the first... The mean value of each segment obtained in step S4 for each RFI frequency channel We construct the mean sequence for this channel. This sequence reflects the average energy level of this frequency channel over different time periods.

[0047] Subsequently, the SumThreshold algorithm was applied to detect anomalies in the mean sequence. SumThreshold is a combined threshold detection method whose core idea is to integrate the sequence using sliding windows of different scales; the larger the window, the lower the corresponding threshold. This allows the algorithm to detect both short-duration, high-intensity RFIs and long-duration, weak RFIs.

[0048] In this embodiment, a series of incrementally increasing window sizes are set. The corresponding threshold sequence is Threshold The calculation formula is usually related to the benchmark threshold. and window size Related, for example, set as:

[0049] in, Based on the statistical characteristics of background noise, The parameter used to control the threshold decay rate is usually set to 1.5.

[0050] The testing process is as follows: First, scan the mean sequence using the smallest window size. If the mean of a certain segment... Greater than the current benchmark threshold If the value at the point marked as RFI is not found, the segment is temporarily replaced with the current threshold before proceeding to the next iteration. To prevent this strong interference point from affecting the statistics of subsequent larger windows, the values ​​at the RFI locations in the sequence are temporarily replaced with the current threshold value. .

[0051] Secondly, increase the window size accordingly and set the corresponding threshold. Using windows Perform a sliding scan of the sequence and calculate the average value of the data within the window. If the window mean is greater than... If the threshold is not met, all unmarked segments within the window will be marked as RFI. Similarly, the values ​​of newly marked anomalous segments will be replaced with the current threshold before the next iteration. .

[0052] Repeat the "scan-mark-replace" process again until the preset maximum window size is reached.

[0053] Using this multi-scale detection strategy, the algorithm ultimately outputs a Boolean sequence indicating which time segments belong to the anomalous mean sequences. These sequences marked as anomalous correspond to the time periods when the actual RFI occurred; while unmarked segments, even if segmented by FPOP, are considered normal background noise fluctuations and are not processed.

[0054] S6: For each segment of outlier mean sequence, subtract its mean to remove RFI. After the SumThreshold detection in step S5, the system has identified which time segments in each RFI frequency channel are actual radio frequency interference and which are normal background. Step S6 aims to clean and recover the raw data based on these detection results.

[0055] In practice, for each frequency data sequence identified as an RFI channel: First, based on the Boolean sequence of anomaly markers output in step S5, the time segment marked as RFI is locked.

[0056] Secondly, backtracking to step S4, we obtain the FPOP segment mean corresponding to this time segment. The mean Physically, it represents the average intensity of RFI energy superimposed on the background noise.

[0057] Secondly, for each original sampling point within this time segment Perform the subtraction operation: For time segments not marked as RFI by S5, the original data remains unchanged. .

[0058] This embodiment employs a mean subtraction strategy instead of the traditional masking and nulling approach. The core value of this method lies in the fact that target signals in radio astronomy observations are often extremely weak and deeply embedded in noise. Using the traditional nulling method would result in the complete loss of data for that time period, disrupting data continuity and spectral characteristics. This invention, by subtracting the segmented mean of the RFI, merely removes the interference components while fully preserving the high-frequency background noise and potential transient astronomical signals hidden beneath the RFI within that time period. This maximizes the preservation of the statistical distribution characteristics of the observation data, providing a more complete and higher-quality data foundation for subsequent searches and analysis of weak astronomical signals.

[0059] like Figure 7 As shown, through the complete processing of steps S1 to S6 above, both broadband and narrowband RFIs in the original observation data are effectively identified and accurately removed by subtracting the segmented mean. Compared with traditional methods, the interference detection accuracy and efficiency of this embodiment are higher.

[0060] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of the invention. Various variations can be made to the above embodiments of the present invention. All simple and equivalent changes and modifications made based on the claims and description of this application fall within the protection scope of the claims of this patent. All aspects not described in detail in this invention are conventional technical content.

Claims

1. A method for suppressing radio astronomy interference based on image segmentation, characterized in that, Includes the following steps: S1: Use the OTSU algorithm to binarize the original time-frequency data to obtain the initial mask data. ; S2: Initial mask data Adaptive filtering is performed along the time and frequency axes. Specifically, this adaptive filtering is implemented as follows: sliding windows are constructed along both the time and frequency axes; the density of the center pixel of each window is calculated; a density threshold is calculated based on the statistical distribution of all pixels in the input data; pixels below this threshold are set to zero; and finally, the union of the filtering results along the time and frequency axes is taken to obtain the density mask data. ; S3: Density mask data The number of retained pixels is counted for each frequency channel, and frequency channels with a value greater than a set threshold are marked as RFI frequency channels; S4: Perform change point detection on the RFI frequency channel to determine potential RFI time periods; S5: For the mean sequence formed by the mean values ​​within each RFI time period, mark the abnormal mean sequence; S6: For each abnormal mean sequence, subtract its mean to remove RFI.

2. The radio astronomy interference suppression method based on image segmentation according to claim 1, characterized in that, The specific implementation method of the adaptive filtering in step S2 is as follows: set the size to... The window slides along the time axis with the initial mask data. Convolution, set size to The window slides along the frequency axis with the initial mask data. Convolution calculates the horizontal density of the center pixel of the window. and vertical density According to horizontal density and vertical density Calculate the adaptive filtering thresholds for the time axis and the frequency axis; when the density of the pixel at the center of the window is less than the corresponding adaptive filtering threshold, set that pixel to 0; take the union of the two adaptive density comparison results as the final density mask. .

3. The radio astronomy interference suppression method based on image segmentation according to claim 2, characterized in that, Horizontal adaptive density threshold and vertical adaptive density threshold for: in, This represents the total number of pixels in the input data. Indicates the first The horizontal density of the center pixel of each window Indicates the first Vertical density of the center pixel of each window.

4. The radio astronomy interference suppression method based on image segmentation according to claim 1, characterized in that, In step S3, the noise threshold is defined as follows: the number of pixels retained in each frequency channel after statistical density filtering, the median and absolute deviation of the median of the number sequence are calculated, and the sum of the median and N times the absolute deviation of the median is defined as the noise threshold.

5. The radio astronomy interference suppression method based on image segmentation according to claim 4, characterized in that, The noise threshold for: in, Represents the median of the sequence. This represents the absolute deviation of the median of the sequence.

6. The radio astronomy interference suppression method based on image segmentation according to claim 1, characterized in that, In step S4, the FPOP function pruning optimal segmentation algorithm is executed on the RFI frequency channel to detect change points and determine potential RFI time periods.

7. The radio astronomy interference suppression method based on image segmentation according to claim 6, characterized in that, In step S4, the FPOP algorithm is used to detect change points in each RFI channel sequence to obtain the mean within the segment. The FPOP algorithm uses the squared error loss function as the fitting cost during iterative solution and sets the penalty cost based on the number of time sampling points of the input data and the statistical variance of the difference time series.

8. The radio astronomy interference suppression method based on image segmentation according to claim 7, characterized in that, The penalty cost for: in, This represents the total number of time sampling points for the input data. This represents the standard deviation of the background noise.

9. The radio astronomy interference suppression method based on image segmentation according to claim 1, characterized in that, In step S5, the SumThreshold algorithm is applied to different mean sequences of the RFI frequency channel to mark outlier mean sequences.

10. The radio astronomy interference suppression method based on image segmentation according to claim 1, characterized in that, In step S6, only the segmented mean of the time series marked as anomalous is subtracted, while the time series not marked as anomalous are left unprocessed, thus preserving the background information.

Images