Analysis of the radio environment at prospective radio astronomy sites using Monte Carlo methods

1.   INTRODUCTION

Radio astronomy is the science of studying radio signals from distant cosmic objects, initiated by the pioneering discovery of extraterrestrial radio sources by Jansky in 1932^[1]. However, these astronomical signals are extremely weak, with strengths only a few millionths of that which a communication system can produce^[2], e.g., the system noise temperatures in the L, S, and C bands are typically in the range of 17–50 K. To capture these weak signals, radio telescopes must be extremely sensitive, and world-class radio observatories are constructed in areas with favorable electromagnetic environments; electromagnetic environmental protection zones have been designated around these sites to minimize interference. Nevertheless, with the development of radio technology and the popularity of electronic devices, the complexity of electromagnetic environments is increasing. To improve the efficiency of observations, astronomers have continued to explore and experiment with various methods to overcome these challenges.

Recognizing the importance of monitoring and evaluation of station site interference, the International Engineering Group, in 2003, developed the test protocol^[3] for radio frequency interference (RFI) for radio astronomy sites. This radio wave environment test method has been widely used in the site selection of international radio telescopes, providing a valuable guide for telescope construction. At the same time, to protect the quiet electromagnetic environment around radio observatories, some international radio telescopes have designated electromagnetic wave protection areas. For example, the United States National Astronomical Observatories has designated the United States National Radio Tranquility Zone (National Radio Quiet Zone, NRQZ)^[4], and formulated various radio wave environmental protection policies to regulate other businesses in the quiet zone. The candidate sites of the Square Kilometer Array (SKA)^[5] in Africa and Australia are large and sparsely populated, naturally have a good electromagnetic environment advantage, and the Australian SKA site has set up a 100 km restricted area. In addition, the Atacama Large Millimeter/submillimeter Array (ALMA) in Chile and the Giant Microwave Radio Telescope (GMRT) in India have similarly set up electromagnetic environment protection zones. At present, radio telescopes in China, such as the Inner Mongolia Daily Image instrument, Qinghai Delingha millimeter-wave Telescope, Shanghai Tianma Telescope, and 500 m Aperture Spherical radio Telescope (FAST), have established corresponding electromagnetic environmental protection zones^[6-7].

Testing the radio environment is essential for maintaining observation efficiency and ensuring data quality during the site selection phase for radio telescopes. However, radio astronomy stations, whether under construction or in operation, often have limited access to real-time signal monitoring and assessments of radio environment stability. This limitation hinders a thorough evaluation of the radio environment at these locations.

2.   MONTE CARLO METHODS AND THE NEED FOR RADIO ENVIRONMENT ANALYSIS AT RADIO OBSERVATORY SITES

The origins of Monte Carlo methods can be traced to the Second World War era, with significant advancements driven by figures such as Newmann, Fermi, Ulam, and Metropolis^[8]. During this period, discussions emerged on applying these methods to statistical physics, transport problems, economic modeling, and other areas. However, due to the limitations of computer technology at the time, the practical application of Monte Carlo methods was restricted, mainly to some demonstration projects. As computer technology improved, increasing computational power facilitated the growth of Monte Carlo methods, broadening their application across various modern fields. Initially military-focused, these methods have expanded into fields such as physics, engineering, and economics^[9]. In finance, they are used for risk management and derivative pricing; in physics, they aid in modeling subatomic particle behavior; in biology, they are employed to simulate the dynamics of complex biological systems^[10]. Additionally, Monte Carlo methods have become integral to machine learning, optimization problems, and engineering design^[11]. Their unique stochastic nature, flexibility, and efficiency make them indispensable in contemporary scientific research and engineering applications.

The Monte Carlo method uses random number generation to emulate the statistical characteristics of complex systems. Its core tenet involves estimating a probability distribution or expected value for a given problem, through the statistical analysis of a substantial set of random samples. The principal advantage of this method is its capacity to address analytically intractable problems, particularly those characterized by high dimensionality or intricate boundary conditions. In these contexts, Monte Carlo methods show excellent performance in rapidly yielding approximate solutions. Therefore, these methods offer a flexible, comprehensive, and statistically meaningful solution for predicting long-term interference patterns in radio environments, especially when dealing with complex physical phenomena and multi-variable issues. Its advantages are particularly evident in these contexts.

In light of the fact that episodically interfering signals in the radio spectrum are significant factors that can potentially impact radio astronomy observations, and considering that the pattern of their episodicity and the duration of their occurrence are crucial data for effectively assessing the radio environment at the station and the quality of the expected observational data, the International Telecommunication Union (ITU) has established a third criterion for the protection of Radio Astronomy Services. This criterion specifies the percentage of time during which the interference level exceeds a threshold without causing serious degradation of the quality of radio astronomy observations. The standard ITU-R RA.1513 proposes a percentage-of-time criterion and an analysis of the associated loss of observations, stating that the cumulative loss of radio astronomy observation time due to all disturbances exceeding the threshold level should not exceed 5% in any of the radio astronomy bands designated for primary operations.

However, assessing the stability of the radio environment for astronomical stations, under construction or in operation, is challenging due to the absence of real-time data from radio environment testing methods. Monte Carlo analysis, which relies on extensive historical data, offers distinct advantages for evaluating the long-term stability of the radio environment. By employing Monte Carlo analysis to predict the occurrence frequency of sporadic interference signals from a limited set of spectrum samples, we can determine the proportion of time during which these signals occupy a specific frequency channel. Comparing this proportion to the ITU-R RA.1513 standard, which stipulates that a usable observing channel should have less than 5% interference, allows for a rapid evaluation of the station's environment and the anticipated quality of observations.

Combining the advantages of radio astronomy and Monte Carlo analysis, this paper proposes a Monte Carlo method to assess the radio environment at astronomical sites, aiming at quality assessments of radio astronomy stations based on a large number of data samples. These samples are acquired through continuous monitoring of the radio environment, providing a reference for the selection of radio astronomy sites and the improvement of the radio environment.

3.   PRINCIPLES AND METHODS OF ANALYSIS

The analysis process in this study is organized into three principal stages: Data collection through the radio environment testing system, determination of interference thresholds and filtering of interference signals, and application of a Monte Carlo simulation for predictive analysis.

A radio environmental test system is first used to capture environmental data to provide raw data for further subsequent analysis and assessment. Next, the threshold calculation algorithm is used to calculate the judgment threshold of interfering signals, to provide a criterion for determining the RFI signals, and to collate those signals within a limited number of samples. Finally, the Monte Carlo method is applied to the radio environment data obtained using this observation, to predict the confidence degree of the interference signal over time.

The above steps are used to estimate the confidence level that the proportion of interference time to total time at a given test site does not exceed 5%. This provides theoretical support and a basis for optimizing the radio environment for radio telescope observations.

3.1   Test System Components

The radio environment test system functions to observe the radio environment data, display the spectrum and waterfall spectrum that need to be observed in real time, and store the data. The system comprises a test antenna, low noise amplifier, voltage stabilizing power supply, spectrometer, radio frequency switch, noise source, and control computer. Table 1 shows the models, parameters and technical indicators of some of the main equipment in the system.

Table  1.  Model and parameters of some equipment^[12]

Equipment	Model	Frequency range	Technical specification
R&S logarithmic periodic antenna	HL050	0.85 MHz–26 GHz	Standing wave: ≤2.5 gain/dBi:8–8.5
Low noise amplifier	BFB0003010	10 MHz–7 GHz	Standing wave: ≤2.5 gain/dBi:30
FSP-7 spectrometer	RS-HL050	20 Hz–7 GHz	Phase noise: 127 dBc/Hz
Noise sources	-	9 kHz–3 GHz	CH power: -6.36 dBm
Regulated power supply	UTP3702	D.C.	Voltage stabilization range: 0–25 V

 | Show Table

DownLoad: CSV

The structure of the test system is shown in Fig. 1. The input signal can either be the signal received from the antenna or the signal from the noise source, and this selection is controlled by the RF switch. When the computer controls the radio frequency switch to connect to the noise source, the system automatically switches into calibration mode, in which the system can determine the present working conditions by comparing spectrograms of calibration with normal working conditions. This ensures the accuracy and reliability of the subsequent test data. When the computer accesses the antenna signal, the HL050 antenna is responsible for collecting the electromagnetic signal in the environment, which is amplified and processed by the low-noise amplifier and then input into the spectrometer to be converted into the frequency domain data. The test data are then transmitted to the computer for further storage and analysis.

Figure  1.  The test system structure, incorporating a receiving antenna, noise source, high-frequency switch, low-noise amplifier, spectrum analyzer, and computer.

DownLoad: Full-Size Img PowerPoint

On the operator side, the system uses the LabView graphical programming software and common interface bus control protocol for system control, test data reception, and data storage. Fig. 2 shows the main interface of the upper computer software, which can be divided into three sections from top to bottom.

Figure  2.  Screenshot and explanation of the graphical user interface of the host computer software.

DownLoad: Full-Size Img PowerPoint

The upper right section (Basic Parameter Setting) is responsible for the input of parameters, for reading data, such as frequency band, observation direction, polarized information, integration time, frequency resolution. The middle section is the real-time waterfall chart display area, separated into Average value waterfall chart (left) and Maximum value waterfall chart (right). The third section, at the bottom, is the real-time spectrum display area, in which the red spectrum line is the maximum value spectrum (Maximum spectrum), and the blue spectrum line is the average value spectrum (Average spectrum). After setting the basic parameters, the "START" button will commence data observation and recording.

During the observation process, the software stores the spectrum monitoring data in a plain text file at a set file path for subsequent analysis of the data. The data storage format is shown in Fig. 3, each line in the file is the spectrum signal within the full bandwidth obtained each time. The date at the beginning of the data shows when it was recorded, and is given in Coordinated Universal Time (UTC), recorded in the format of "year-month-day hour: minute: second", and the x at the end of the date gives the power of the frequency band in units of dBm.

Figure  3.  The data storage format, with each row containing spectroscopic data including the time of recording. Each spectrum contains 4000 sampling points, in the frequency range of 1–2 GHz.

DownLoad: Full-Size Img PowerPoint

3.2   Threshold Calculation

Threshold calculation is a key step to extract and determine the RFI in a large number of data samples. Here, in the threshold calculation, the collected power data are first converted into "log-linear" data, which aims to retain the original characteristics and relevance of the data, while also effectively reducing the covariance and anisotropy in the model. This can be calculated as^[13]

where ${x}_{i}$ is the observed power value in units of dBm and X_i is the log-converted data in units of W. The sample mean and the variance of the linear data need to be calculated. The sample mean can be calculated using the numerical mean method with the expression

The estimation of the sample variance is then based on the Simple Thresholding algorithm. Calculating the variance using the Simple Thresholding algorithm involves subtracting each element of the original sample, using the median of the original data, to create a new modified sample of the same size as the original elements. The median of this new sample was then calculated again and multiplied by a constant scale factor ($u=1.4268$) to make this estimate consistent with the expected Gaussian distribution^[14]. Thus, the variance MAD of the sample is expressed as

where med is the median of the original data.

For the calculation of the linear threshold, a sample must be considered to exceed the threshold when the sample value exceeds the sample mean times the variance. Therefore, after obtaining the sample mean and variance, the linear threshold Th (in units of W) is expressed as

In this case, the linear threshold is transformed into a "linear-logarithmic" value to obtain the final threshold result Th* (in units of dBm), given by

If the power level of the selected channel exceeds the threshold Th*, it is considered to be RFI. The analysis flowchart and algorithm pseudo-code are shown in Fig. 4 and Fig. 5.

Figure  4.  Flowchart showing the operation of the thresholding algorithm.

DownLoad: Full-Size Img PowerPoint

Figure  5.  Pseudo-code for the thresholding algorithm.

DownLoad: Full-Size Img PowerPoint

3.3   Introduction to the Monte Carlo Method

The Monte Carlo method, also known as the random sampling method or statistical simulation method, is a method of calculation and numerical simulation using the statistical laws of random numbers. Problem solving using Monte Carlo methods can be roughly divided into two steps^[15].

(1) Build a model of random variables for a physical problem. The solution to the problem is given by the probability of the variable or the mathematical expected value.

(2) Perform numerical simulation experiments with computers, i.e., random sampling of the model to produce values of random variables, so that the probability or average of the amount of change can be derived as an approximate solution to the problem.

Monte Carlo analysis is a statistical method comprising multiple independent experiments, where each experiment's parameters are randomly selected based on a predefined probability distribution. The data from the portion of experimental results that fail to meet the criteria represent the likelihood of such failures occurring in real-world scenarios. In this study, assuming that we take k spectra as the raw data of one experiment, and each spectrum contains j frequency points, the number of frequency points Num present in one experiment is

Comparing the above Num frequency point values with the threshold value, if the proportion of the number of frequency points in the Num frequency point data that do not exceed the interference threshold is more than 95% of the total number of frequency points, the group of data is considered to be in line with the requirements of radio-astronomy observation, and it is recorded as a successful experiment. Otherwise, it is recorded as a failed experiment. By analyzing a large amount of data using the above method, we find that among the N groups of data samples, there are n groups of data samples that are unqualified.

According to statistical theory, the probability value p of any set of data non-compliance tends toward n/N. When N tends to infinity, p tends to the true value. For a finite number of total data samples, N, and a given probability, p, of the occurrence of interference, there exist n sets of data samples that do not meet the requirements. This is to say, the probability of the samples of the interference signals conforms to the Burdock's effort distribution,

However, for practical purposes, it is difficult to obtain the true probability of the occurrence of interference through a finite number of test samples, so the probability distribution, P_n(p), of p is obtained by calculating the number of groups, n, of interference sample data in a finite set of sample data. This is an inverse probability problem, according to Woodward's product law of probability^[16],

where P(n) and $P\left(p\right)$ are the a priori probabilities of n and p, respectively. Since n can be obtained using the method in Fig. 6, P(n) is a constant value. Additionally, since there is a low probability of interference in the observed frequency bands that have passed the preliminary selection and exceed the interference threshold, $P\left(p\right)$ can be considered a very small value in the actual observation, allowing Equation (8) to be rewritten as

Figure  6.  Chart showing data format changes. After dividing the data into N groups, the frequency points of each group are compared with the interference threshold to determine the proportion of frequency points exceeding the threshold. If the proportion is greater than 5%, the data in this group are not considered to meet the observation standards.

DownLoad: Full-Size Img PowerPoint

where k is a constant. With the normalization of the probabilities, the sum of the probabilities is always 1, and therefore,

Integrating both sides of Equation (9) simultaneously and obtaining the integral equality given by

collating these gives

This is further proved using the division integration method, giving

Combining Equation (12) and Equation (13) to find k = N+1, and substituting into Equation (9), the inverse probability ${P}_{n}\left(p\right)$ is

Up to this point, the mathematical expression for the probability distribution of p, given a specific value of n, is shown in Equation (14). Where $p\in \left[\mathrm{0,1}\right]$, in the previous section, we have derived the values of N and n by analyzing the test data from the radio environment, allowing ${P}_{n}\left(p\right)$ to be found.

For analytical purposes, Equation (14) is further treated here. Assuming that the actual value of p is required to be no greater than the set value p_c for C%, then P_n(p) must satisfy the integral inequality,

As a result, ${P}_{n}\left(p\right)$, p_c, and C% provide a basis for determining whether this equation is valid or not. If the inequality is valid, the interference in the experimental frequency band in the radio environment is considered to meet the radio astronomy interference threshold percentage requirement. If not, the interference in the radio environment needs to be managed or the observation data needs to be cleaned.

4.   DATA ANALYSIS

4.1   Calculation of Interference Thresholds

The National Radio Administration Bureau has given a specific division of frequency bands for radio astronomy operations, and according to the Provisions of the People's Republic of China on the Division of Radio Frequencies, the protected frequency bands for radio astronomy in the L-band (1–2 GHz), and the corresponding observation purposes, are shown in Table 2.

Table  2.  L-band radio astronomy protection frequency band

Observation band/MHz	Main operations in the frequency band	Main scientific objectives
1330–1440	Aviation radionavigation, radiolocation, Satellite radionavigation	Neutral chlorine and Pulsar survey, Transient objects, Extragalactic radio sources
1610–1613	Satellite mobile (Earth-space), Aeronautical radio navigation, satellite radio determination (Earth-space)	Pulsars, VLBI extraterrestrial life, Neutral chlorine, Molecular spectral lines
1660–1670	Satellite mobile (Earth-space), Space research (passive), Meteorological aids, Fixed, Mobile (except air mobile)	Pulsars, VLBI extraterrestrial life, Neutral chlorine, Molecular spectral lines
1718–1722	Mobile communications	Pulsars, VLBI extraterrestrial life, Neutral chlorine, Molecular spectral lines

 | Show Table

DownLoad: CSV

This study mainly analyzes the following radio protected frequency bands in the L-band: 1330–1440 MHz, 1610–1613 MHz, 1660–1670 MHz, and 1718–1722 MHz. We observed the L-band radio environment at Yunnan Observatories, Chinese Academy of Sciences, for a long period of time. The integration time of all the data during the observation is uniformly adopted as 3 s. We randomly choose 99600 spectral data points among the 20 days of observation data as samples for the present study. Subsequently, considering the volume of observational data and the speed and memory capacity of computer calculations, we extracted 400 spectra for threshold analysis, giving threshold spectra. The quantity of these data is sufficient for our tests. Using a larger quantity of spectral data increases the data volume for individual experiments, allowing more accurate determination of whether a single experiment fails, consequently enhancing prediction accuracy.

The threshold spectra of the four frequency bands are shown in Fig. 7. Fig. 7A appears rougher because of the relatively large size of the loan. The band thresholds range from −96.5 to −96 dBm for the frequency band 1330–1440 MHz, and the band thresholds range from −96.7 to −96.2 dBm for the frequency bands 1610–1613 MHz, 1660–670 MHz, and 1718–1722 MHz. The band thresholds range from −96.7 to −96.2 dBm. These thresholds provide an important basis for threshold judgment of subsequent data.

Figure  7.  Results of the threshold analysis, with the interference threshold calculated for four radio protection frequency bands by analyzing 400 spectra.

DownLoad: Full-Size Img PowerPoint

Fig. 8 illustrates the distribution of interfering signals in the 400 sample spectra. Fig. 7A, Fig. 7B, Fig. 7C and Fig.7D show the interference in the 1330–1440 MHz, 1610–1613 MHz, 1660–1670 MHz, and 1718–1722 MHz bands in this set of data, respectively, in which the frequency points with interference exceeding the above interference thresholds have been marked using red dots. The percentages of interference frequency points in these bands are 0.84%, 0.67%, 0.87%, and 0.67%, respectively, none of which exceeded 5%. This indicates that this set of data meets the observation requirements.

Figure  8.  Distribution of interference signals extracted from the first set of data, appearing in the four protected radio frequency bands.

DownLoad: Full-Size Img PowerPoint

4.2   Monte Carlo Analysis of Actual Data

The standard ITU-R RA.1513 specifies that, for radio astronomy operations, the maximum allowable value for the percentage of net time lost due to interference is 5%. Therefore, in the results of the Monte Carlo analysis, if the interference in a given data sample exceeds the threshold value by more than 5% of the total data volume of that sample, the data set is considered to have too much interference to meet the observational requirements.

In this study, 20-minute observations are used as a set of data samples, and each set of data samples contains 400 integrated spectra, each with an integration time of 3 s. After exhaustive processing of all observations, we obtain the experimental results shown in Table 3. In a total of 249 sets of data samples, the failure rate varies in different frequency bands. The frequency band 1330–1400 MHz has 4 groups of data samples failing to meet the observation requirements in the test, while three times the number of groups of data samples in the frequency bands 1610–1613 MHz and 1660–1670 MHz fail to meet the observation requirements, and the number of failures in the frequency band 1718–1722 MHz is significantly increased by 22 times.

Table  3.  Experimental results

Frequency/ MHz	Total number of experiments	Number of failed experiments	n/N
1330–1400	249	4	0.016
1610–1613	249	3	0.012
1660–1670	249	3	0.012
1718–1722	249	22	0.088

 | Show Table

DownLoad: CSV

When the number of experimental failures, n, is large, the confidence distribution is concentrated in the interval of higher failure probabilities. This implies a high certainty of interference occurring in this frequency band, causing an observation efficiency insufficient to meet the requirements. Conversely, when the number of experimental failures, n, is small, the confidence distribution is concentrated in the interval of lower failure probabilities, indicating a low probability of RFI occurring in the next experiment, which is more in line with observational requirements.

Combining the probability distribution of interference occurrence mentioned earlier, the deterministic probability distribution represents the relationship between the probability of interference exceeding the threshold and its confidence level. By applying the statistical results from Table 3 to Equation (14), we obtain the probability and confidence level relationship curves for four frequency bands, as shown in Fig. 9. When n is small, the confidence level of interference exceeding the threshold is concentrated in the low-probability interval. For example, the peak of the blue curve in Fig. 9 is at (0.013, 56). This means that the interference occurrence probability for the frequency band 1660–1670 MHz is 1.3%, with a likelihood of 56%.

Figure  9.  Deterministic probability distribution.

The horizontal axis shows the probability of interference occurrence, and the vertical axis shows the confidence level of interference occurrence probability. The consistent results in the 1610–1613MHz and 1660–1670MHz bands cause the curves to overlap.

DownLoad: Full-Size Img PowerPoint

For radio astronomy observations, we require the probability of interference occurrence not to exceed a certain threshold. Therefore, we integrate the values in Equation (14), i.e., Equation (15), to obtain the deterministic probability density curve, as shown in Fig. 10. As mentioned earlier, the ITU-R RA.1513 standard states that for radio astronomy protection bands, the interference should not exceed 5% of the total time with a certainty of not less than 90% to be considered as meeting the requirements for radio astronomy observations. Using the deterministic probability density curve in Fig. 8B, we can draw the following two conclusions:

Figure  10.  The accumulation of deterministic probability distributions. The confidence level represents the certainty when the probability of interference occurrence is lower than a certain value. The consistent results in the 1610–1613MHz and 1660–1670MHz bands cause the curves to overlap.

DownLoad: Full-Size Img PowerPoint

(1) The interference levels in the 1330–1440 MHz, 1610–1613 MHz, and 1660–1670 MHz bands are in accordance with the standards proposed in ITU-R RA.1513, with respective confidence levels of 96.9%, 97.4%, and 97.4% that the time share of their interfering signals does not exceed the stipulated 5% of the total observing time. This indicates that these frequency bands are both stable and reliable for radio astronomy observations.

(2) The 1718–1722 MHz band shows a long period of intermittent interference, with an 88.5% confidence level that the time share of the disturbed signals in this band did not exceed 5% of the total observation time, as specified by ITU. This interference will impact the efficiency and observation quality of radio telescopes. Therefore, it is necessary to further clean up and optimize the radio environment in this frequency band.

5.   CONCLUSIONS

This paper introduces a methodology for evaluating the stability of the radio environment at a radio observatory site. Our approach involves conducting extensive, long-term observations to gather data on the radio environment of the site. Using the interference threshold method, we calculate thresholds for specific frequency bands to filter the temporal distribution of interfering signals. Subsequently, Monte Carlo analysis is employed to ascertain the likelihood that the duration of interference exceeds a certain fraction of the total observation period within the radio environment on-site. This methodology lays a foundation for establishing radio astronomy facilities and assessing their radio environments.

This study also applies Monte Carlo analysis and evaluates the stability of the radio environment for four protected frequency bands within the L-band. The findings indicate that the interference levels in the bands of 1330–1440 MHz, 1610–1613 MHz, and 1660–1670 MHz are in compliance with the standards set forth in the ITU-R RA.1513. However, the band at 1718–1722 MHz exhibits an extended period of intermittent interference, necessitating further mitigation and optimization efforts.