RESEARCH ARTICLE
10.1029/2023RS007697
Key Points: • Meteoroid trajectory and speed
determined with a continuous wave forward scatter radio set-up • Trajectory reconstruction solver is validated against simulated data but the problem is ill-conditioned • Interferometry required for an accuracy of <10% on the velocity and <2° on the inclination for most trajectories
Correspondence to: J. Balis, joachim.balis@aeronomie.be
Citation: Balis, J., Lamy, H., Anciaux, M., & Jehin, E. (2023). Reconstructing meteoroid trajectories using forward scatter radio observations from the BRAMS network. Radio Science, 58, e2023RS007697. https://doi.org/10.1029/2023RS007697
Received 28 FEB 2023 Accepted 25 MAY 2023

Reconstructing Meteoroid Trajectories Using Forward Scatter
Radio Observations From the BRAMS Network
Joachim Balis1,2 , Hervé Lamy1 , Michel Anciaux1 , and Emmanuel Jehin2
1Royal Belgian Institute for Space Aeronomy, Brussels, Belgium, 2STAR Institute, University of Liège, Liège, Belgium
Abstract  In this paper, we aim to reconstruct meteoroid trajectories using a forward scatter radio system
transmitting a continuous wave (CW) with no modulation. To do so, we use the meteor echoes recorded at the receivers of the BRAMS (Belgian RAdio Meteor Stations) network. This system consists, at the time of writing, of a dedicated transmitter and 44 receiving stations located in and nearby Belgium, all synchronized using GPS clocks. Our approach processes the meteor echoes at the BRAMS receivers and uses the time delays as inputs to a nonlinear optimization solver. We compare the quality of our reconstructions with and without interferometric data to the trajectories given by the optical CAMS (Cameras for Allsky Meteor Surveillance) network in Benelux. We show that the general CW forward scatter trajectory reconstruction problem can be solved, but we highlight its strong ill-conditioning. With interferometry, this high sensitivity to the inputs is alleviated and the reconstructed trajectories are in good agreement with optical ones, displaying an uncertainty smaller than 10% on the velocity and 2° on the inclination for most cases. To increase accuracy, the trajectory reconstruction with time delays only should be complemented by information about the signal phase. The use of at least one interferometer makes the problem much easier to solve and greatly improves the accuracy of the retrieved velocities and inclinations. Increasing the number of receiving stations also enhances the quality of the reconstructions.
Plain Language Summary  This paper presents a method for tracking the path and speed of
meteoroids using radio observations. A simple continuous wave signal is reflected by the electrons created when the meteoroid enters the upper atmosphere and creates ionization. The reflected signal, called a meteor echo, is recorded at various locations not co-located with the transmitter. The time delays between several echoes is used to retrieve the meteoroid trajectory and speed. In this work, we apply this method using data from Belgian RAdio Meteor Stations (BRAMS), a Belgian network of 44 receiving stations. The results are compared to accurate observations from Cameras for Allsky Meteor Surveillance, an optical network of cameras surveying the sky for meteors. Our method perfectly works with no measurement errors but small uncertainties on the time delays may significantly impact the accuracy of the reconstruction. Using a large amount of receivers and/or a radio interferometer greatly improves the results. This project is the first to tackle the meteoroid trajectory reconstruction using this type of radio system (with no range information and many receivers at different locations). This novel method is essential to fully exploit the capabilities of BRAMS for future applications such as determination of fluxes or sounding of the upper atmosphere (e.g., wind speed measurements).

© 2023. The Authors. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

1. Introduction
The BRAMS (Belgian RAdio Meteor Stations) network is a Belgian project using forward scatter of radio waves to detect and characterize meteoroids. It comprises a dedicated transmitter located in the South-West of Belgium and 44 receiving stations spread all over the Belgian territory and neighboring countries (see Figure 1 for status in February 2023).
The transmitter emits a circularly polarized continuous radio wave with no modulation at a frequency of 49.97 MHz with a power of 130 W. All the receiving stations are using a 3-element Yagi antenna set up vertically and oriented in azimuth toward the transmitter. At the time of writing, approximately a third of the receiving stations is using analog ICOM-R75 receivers, an external sound card to sample the signal coming from the antenna, and is controlled by the freeware program Spectrum Lab running on a PC (Lamy et al., 2015). The other two thirds use digital SDR-RSP2 receivers controlled by a Linux system running on a Raspberry Pi (Anciaux et al., 2020). All stations are equipped with a Garmin GPS that provides timestamps to the BRAMS data and

BALIS ET AL.

1 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

allows time synchronization between the receiving stations. Additional features of the receiving stations are not described here. Instead, we refer the reader to previous publications (Anciaux et al., 2020; Lamy et al., 2015).

There exist other meteor forward scatter systems but they are either asso-

ciated with meteor radars (sending pulses), or continuous wave (CW)

emitting systems with additional coding (e.g., phase coding). An applica-

tion of such forward scatter systems is wind speed measurements in the

mesosphere/lower-thermosphere (MLT) region, which have seen important

developments in the last few years. Stober and Chau (2015) used data from

the MMARIA (Multistatic/Multifrequency Agile Radar for Investigations of

the Atmosphere) concept, made of two meteor radars and another receiver

providing an additional forward scatter capability. They showed that the

forward scatter observations increased both the number of meteor detections

and the altitudinal coverage. Vierinen et al. (2016) introduced the concept

of coded CW meteor radars and presented their advantages with respect

to pulsed specular meteor radars (SMR). Notably based on this work, the

MMARIA concept was further developed into SIMONe (Spread-spectrum

Figure 1.  Map of the BRAMS network in February 2023. The blue triangle is the transmitter located in Dourbes while the green dots are the 44 active

Interferometric multistatic Meteor radar Observing Network) to extend the capabilities of traditional meteor radar systems by adding either several

receiving stations at the time.

spatially separated receivers to existing transmitting stations, or several

spatially separated transmitters to existing interferometric receiving stations

(Chau et al., 2019). One unique feature of this work is the use of interferomet-

ric transmitting arrays as opposed to interferometric receiving arrays typical in meteor radars. Another example

application is the measurement of ozone concentration in the MLT region, which has been indirectly performed

by the Bologna-Lecce-Modra forward scatter meteor radar, thanks to the detection of meteoroids entering the

Earth's upper atmosphere (Cevolani & Pupillo, 2009).

To the authors' knowledge, no existing CW forward scatter system without modulation has tackled the problem of reconstructing individual meteoroid trajectories and speeds. One of the main difficulties stems from the geometry, which is more complex than in the case of backscatter systems. For instance, as BRAMS does not modulate its CW transmitted signal, the total range traveled by the radio wave between the transmitter (Tx), the reflection point and the receiver (Rx) cannot be estimated. Therefore, the trajectory reconstruction problem is complicated to solve.

In this paper, we present an approach to retrieve meteoroid trajectories with a forward scatter CW set-up such as BRAMS. Two methods are considered: one based only on measurements of time delays between meteor echoes recorded at different receiving stations, and one using the same data but complemented with information from the radio interferometer located in the Humain station which provides the direction of one specular reflection point. In order to assess the quality of the reconstruction, a comparison with data from the optical CAMS-BeNeLux network is provided (Roggemans et al., 2016).

2.  Mathematical Model
The first method (hereafter called Method 1) is based solely on geometrical considerations and relies on the specular condition of the reflection of the radio wave. The specular reflection point is the location along the meteoroid path for which the total distance traveled by the radio wave Si = RTi + RRi is minimum. In this expression, RTi is the distance from the transmitter to the meteor and RRi is the distance from the meteor to the receiver. The condition of specular reflection allows the position of the reflection point to be identified as the point on the ellipsoid, whose foci are the transmitter and the receiver, which is tangent to the meteor path. Because the geometry Tx − Rxi is different for each receiving station Rxi, the corresponding reflection points will be located at various positions along the meteoroid path. This is illustrated in Figure 2 for a reference station Rx0 and another station Rxi. In this example, the specular reflection point P0 for the reference station is created before the corresponding reflection point Pi. The distance between the two points depends on the speed V of the meteoroid. As a consequence, the reference station will detect a meteor echo shortly before receiving station i, the time delay Δti between meteor echoes depending on the meteoroid path and speed.

BALIS ET AL.

2 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

Throughout this work, we assume that the trajectory of the meteoroid body can be described by an uniform rectilinear motion. The rectilinear assumption can be justified by the fact that the influence of the gravitational acceleration and the Coriolis forces is small compared to the magnitude of the meteor velocity: the maximum radius of curvature induced by these accelerations is respectively 2 and 12 times the radius of the Earth (Jeanne, 2020). Atmospheric drag and mass-loss can introduce significant deceleration (Kastinen & Kero, 2022; Schult et al., 2017), even for the small meteoroids detected by BRAMS. However, it is reasonable to assume that the majority of the deceleration occurs at the end of the ablation process. As most of the exploitable echoes come from the beginning and middle of the meteor event, the time delays are thus measured in a region where the constant speed assumption is acceptable.

Figure 2.  Specularity and geometry of a forward scatter set-up. The time delay between the creation of the specular point P0 of the reference station Rx0 and the specular point Pi of the other station Rxi is noted ti. The straight blue arrow indicates the meteor path, traveling at a speed V. The two radio wave paths Tx − P0 − Rx0 and Tx − Pi − Rxi are not coplanar.

2.1.  Forward Problem
Before developing a method for retrieving the trajectory from the time delays observed at the receivers, let us study the forward problem, that is, determine the location of the specular points for a given meteor trajectory. This problem

comes down to finding the intersection between an ellipsoid and a straight

line (Nedeljković, 2006). If the start point is P1 = (x1; y1; z1) and the end point is P2 = (x2; y2; z2), the equations of the meteor trajectory are given by:

 𝑥𝑥𝑥𝑥1−−𝑥𝑥𝑥𝑥12 = 𝑦𝑦𝑦𝑦1−−𝑦𝑦𝑦𝑦12 = 𝑧𝑧𝑧𝑧1−−𝑧𝑧𝑧𝑧12 . (1)

The description of the sought-for ellipsoid is much more complex but, given a coordinate system in which the Cartesian axes are also those of the ellipsoid, it can be described by:
 𝑥𝑎𝑎𝑥22 + 𝑦𝑏𝑏𝑦22 + 𝑧𝑐𝑐𝑧22 = 1, (2)
where a, b, and c are real positive numbers. To facilitate use of this property, a new system of coordinates is introduced whose origin is located at the geometric center point PTR = (xTR; yTR; zTR) between the transmitter and the receiver. We define the y-axis in the same direction as the vector joining these points. The considered ellipsoid is a spheroid (i.e., an ellipsoid of revolution) because the two ranges RTi + RRi are fixed but the location of the reflection point is arbitrary. As a result, we deduce a = c, such that Equation 2 can be expressed as:
 𝑥𝑎𝑎𝑥22 + 𝑦𝑏𝑏𝑦22 + 𝑎𝑧𝑧𝑎22 = 1. (3)

This system can be obtained with the following transformations: a translation to set the location of the ellipsoid center PTR as the origin of the reference frame (conversion between system 0 and ′):
𝑥𝑥′ = 𝑥𝑥0 − 𝑥𝑥𝑇𝑇 𝑇𝑇,  𝑦𝑦′ = 𝑦𝑦0 − 𝑦𝑦𝑇𝑇 𝑇𝑇, (4)
𝑧𝑧′ = 𝑧𝑧0 − 𝑧𝑧𝑇𝑇 𝑇𝑇,
() a rotation around the z-axis by an𝐴a𝐴ngle 𝐴𝐴 = arctan 𝑥𝑥𝑇𝑇𝑇𝑇 (angle between system ′ and ″):
𝑦𝑦𝑇𝑇 𝑇𝑇
𝑥𝑥′′ = cos(𝛼𝛼)𝑥𝑥′ − sin(𝛼𝛼)𝑦𝑦′,  𝑦𝑦′′ = sin(𝛼𝛼)𝑥𝑥′ + cos(𝛼𝛼)𝑦𝑦′, (5)
𝑧𝑧′′ = 𝑧𝑧′,

BALIS ET AL.

3 of 18

BALIS ET AL.

Radio Science

10.1029/2023RS007697

() and a rotation around the x-axis by an𝐴a𝐴ngle 𝐴𝐴 = −arctan 𝑧𝑧𝑇𝑇𝑇𝑇 (angle between system ″ and ‴):
𝑦𝑦𝑇𝑇 𝑇𝑇
𝑥𝑥′′′ = 𝑥𝑥′′,
 𝑦𝑦′′′ = cos(𝛽𝛽)𝑦𝑦′′ − sin(𝛽𝛽)𝑧𝑧′′, (6)
𝑧𝑧′′′ = sin(𝛽𝛽)𝑦𝑦′′ + cos(𝛽𝛽)𝑧𝑧′′.
This last reference frame is the one in which the determination of the specular point location is performed. In the following, we will drop the superscript  ‴ for clarity. Thus, x, y, z should be understood as x′′′, y′′′ and z′′′ for the rest of this subsection.
Using Equation 1, we can deduce x and z values from y and obtain the following system of equations:  𝑧𝑥𝑧𝑥 == 𝑧𝑦𝑥𝑦𝑧𝑦𝑥𝑦1111 −−−− 𝑧𝑦𝑥𝑦𝑦𝑧𝑦𝑥2222((𝑦𝑦𝑦𝑦−−𝑦𝑦𝑦𝑦11))++𝑧𝑥𝑧𝑥11., (7)
Using these values in Equation 3, we obtain a quadratic equation of the form Ay 2 + By + C = 0, whose solutions are given by:
√  𝑦𝑦1,2 = −𝐵𝐵 ± 2𝐵𝐴𝐵𝐴2 − 4𝐴𝐴𝐴𝐴 . (8)
Because we are looking for a tangent point, we already know that there the solution is unique, which means that:
 𝐵𝐵2 − 4𝐴𝐴𝐴𝐴 = 0, (9)
where A, B, and C are functions of a, b, P1, and P2. Replacing a 2 by b 2 − f 2 (where f is half the interfocal distance of the ellipsoid), this gives us a fourth order equation in b. This equation can be solved for b 2 and then for a 2 = b 2 − f 2 with a and b being chosen so that they are positive values. Knowing a and b, we are able to compute A, B, and then y. From Equation 7, we find the full position vector PTan = (x; y; z) in the system  ‴. The final step is to perform an inverse transformation to come back to an easy-to-use and more general coordinate system (i.e., not linked to a specific ellipsoid). It is natural to choose Dourbes, where the BRAMS transmitter Tx is located, as the origin of this global reference system. Like in previous work on the BRAMS network (Lamy & Tétard, 2016), the West-East direction is used as the X-axis and the South-North direction as the Y-axis. The Z-axis is taken out of the plane by the right-hand rule. This reference frame is used for our following computations.

2.2.  Inverse Problem

We have defined a way to compute the location of the specular points, knowing the meteoroid trajectory. In practice however, the trajectory is unknown. To find it, the only inputs at our disposal with Method 1 are the time delays between meteor echoes observed at the receiving stations. The trajectory reconstruction problem is thus inverse.

More formally, a meteoroid trajectory can be defined by the 3D Cartesian coordinates of one specular point (the one corresponding to a reference station) and the three components of the velocity which provide the direction (assuming again a constant speed). This gives a total of six unknowns (respectively called X0, Y0, Z0, Vx, Vy, and Vz) and therefore the need to have at least six equations to avoid solving an underdetermined system. In Method 1, these equations are provided by the specularity condition. Let us write the distance Si, traveled by the radio wave from the transmitter toward a reflection point Pi = (Xi, Yi, Zi) and back to the receiver Rxi = (xi, yi, zi):

√

√

 𝑆𝑆𝑖𝑖 = 𝑋𝑋𝑖2𝑖 + 𝑌𝑌𝑖𝑖2 + 𝑍𝑍𝑖2𝑖 + (𝑋𝑋𝑖𝑖 − 𝑥𝑥𝑖𝑖)2 + (𝑌𝑌𝑖𝑖 − 𝑦𝑦𝑖𝑖)2 + (𝑍𝑍𝑖𝑖 − 𝑧𝑧𝑖𝑖)2, (10)

where, for each receiver Rxi:

⎛ 𝑋𝑋𝑖𝑖 ⎞ ⎛ 𝑋𝑋0 ⎞

⎛ 𝑉𝑉𝑥𝑥 ⎞

⎜ ⎟⎜ ⎟ ⎜ ⎟
 ⎜ 𝑌𝑌𝑖𝑖 ⎟ = ⎜ 𝑌𝑌0 ⎟ + Δ𝑡𝑡𝑖𝑖⎜ 𝑉𝑉𝑦𝑦 ⎟, (11) ⎜ ⎟⎜ ⎟ ⎜ ⎟

⎜ ⎝

𝑍𝑍𝑖𝑖

⎟ ⎠

⎜ ⎝

𝑍𝑍0

⎟ ⎠

⎜ ⎝

𝑉𝑉𝑧𝑧

⎟ ⎠

4 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

BALIS ET AL.

Radio Science

10.1029/2023RS007697

which formalizes the fact that there is a time delay Δti between the creation of the specular point Pi of the receiver Rxi and the specular point P0 of the receiver Rx0 (see Figure 2). The specularity condition implies that the total distance traveled by the radio wave is minimum at the specular reflection point,𝐴h𝐴ence 𝑑𝑑𝑑𝑑𝑖𝑖 = 0 . This time deriv-
𝑑𝑑𝑑𝑑
ative can be expressed as:

 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑖𝑖 = 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 𝑖𝑖 + 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑖𝑖 + 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 𝑖𝑖 . (12)

The partial derivatives with respect to the specular point coordinates can be obtained through (example given for Xi, the same principle can be applied for Yi and Zi):
 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖 = √𝑋𝑋𝑖2𝑖 +𝑋𝑋𝑌𝑌𝑖𝑖𝑖𝑖2 + 𝑍𝑍𝑖2𝑖 + √(𝑋𝑋𝑖𝑖 − 𝑥𝑥𝑖𝑖)2 + 𝑋(𝑋𝑌𝑌𝑖𝑖𝑖𝑖−−𝑥𝑦𝑥𝑦𝑖𝑖𝑖𝑖)2 + (𝑍𝑍𝑖𝑖 − 𝑧𝑧𝑖𝑖)2 . (13)

Using the assumption of constant speed motion, we𝐴𝐴have 𝑑𝑑𝑑𝑑𝑖𝑖 𝐴=𝐴 𝑉𝑉𝑥𝑥 , 𝑑𝑑𝑑𝑑𝑖𝑖 = 𝑉𝐴𝑉𝑦𝐴𝑦 and 𝑑𝑑𝑑𝑑𝑖𝑖 = 𝑉𝑉𝑧𝑧 . Combining Equa-

𝑑𝑑𝑑𝑑

𝑑𝑑𝑑𝑑

𝑑𝑑𝑑𝑑

tions 12 and 13 and using the fact that Si should be minimum at each receiving station, we obtain:

= ⋅ +⋅ + ⋅ √ 2+ 2+ 2

 + ( − ) ⋅ + ( − ) ⋅ + ( − ) ⋅ (14)

√ (

− )2 + ( − )2 + (

− )2

= 0.

As Xi, Yi, and Zi depend on the input time delay Δti through Equation 11, we have an implicit link between the input of our system and the output trajectory, described by our six parameters. This set of equations is highly nonlinear and is therefore particularly difficult to solve.

A nonlinear solver that takes into account additional constraints on the unknowns is then used. On the one hand, the height of all reflection points must lie between for example, 80 and 120 km altitude. Above this interval, the atmospheric density is too low for the radio signals to be reflected by the ionized trail. Below it, the low-mass meteoroids that we are interested in are completely disintegrated into dust. On the other hand, the speed of most meteoroids is larger than ∼11 km/s but smaller than ∼72 km/s. The lower limit corresponds to the Earth's escape velocity, while the upper limit is the maximum geocentric speed at which a Solar System object can enter the Earth's atmosphere. Meteoroids of interstellar origin may pass through our atmosphere with a speed of more than 72 km/s, but these cases are very rare (Froncisz et al., 2020; Hajduková Jr. et al., 2014).

The second method (hereafter called Method 2) is using the same assumptions as Method 1 but includes data from our interferometric radio station located in Humain. Unlike the other receiving stations, it uses 5 antennas in the so-called Jones configuration (Jones et al., 1998; Lamy et al., 2018) and allows to determine the direction of arrival of the meteor echo to within approximately 1°, although errors on azimuth and elevation have different dependencies (Jacobs & Ralston, 1981). The interferometer provides two more equations for the azimuth and elevation of the specular reflection point but does not provide its exact position. With these additional equations, we theoretically need time delays measured between only 3 additional stations and a reference station in order to get at least 6 equations.

3.  Solver Development
The nonlinear optimization solver that we will use is based on an interior-point (IP) method (Karmarkar & Ramakrishnan, 1991). The IP algorithms are used in solving both linear and nonlinear convex optimization problems that contain inequalities as constraints. In our case, these inequality constraints are related to the altitude of the reflection points and the meteoroid speed. There are two important IP algorithms: the barrier approach and the primal-dual method. The primal-dual approach is preferred here due to its efficiency and accuracy. Further information about the IP approach can be found in Boyd and Vandenberghe (2004).

3.1.  Target Definition
The most important quantity to define is the target objective (or cost function) to minimize. This function associates a scalar to every trajectory (defined by our 6 unknowns): the lower its value, the better the reconstruction (given the input time delays that we have at our disposal).

5 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

BALIS ET AL.

Radio Science

10.1029/2023RS007697

For Method 1, the first dimensionless target objective that we defined comes from the specularity condition:

 0 = 1 ∑ ⋅ 1 2, (15)
=1

√

𝐴w𝐴here 𝐴𝐴 =

2
𝑉𝑉𝑥𝑥

+

2
𝑉𝑉𝑦𝑦

+

2
𝑉𝑉𝑧𝑧

is

used

to

make

the

target

objective

dimensionless,

and

the

sum

is

done

on

n

receiv-

ers. This approach is the most intuitive as it uses the minimization of the distance traveled by the radio wave. Nonetheless, it shows two major drawbacks. First, it is very sensitive to the quality of the initial guess. Indeed, as the IP approach is gradient-based, we need to prescribe a first solution to the solver. We observed that this guess needs to be quite close to the correct solution to ensure solver convergence, which is problematic given we generally have no a priori knowledge of the trajectory.

The second drawback is that this target does not explicitly use the measured time delays at our receivers. However, the idea of comparing the measured physical quantities to forward modeled ones through use of their uncertainty distributions is powerful, as it allows the construction of a Bayesian inference approach to model parameter estimation. This approach tends to be robust as it includes the errors at their source, utilizes maximal information (no information is lost through further processing), quickly reveals model limitation to fitting the data, and allows for explicit evaluation of how well the model parameters can be determined for each case (Kastinen & Kero, 2022; Markkanen et al., 2005).

To alleviate these drawbacks, we introduce another target objective for Method 1, defined as:

 𝐽𝐽1 = ∑ 𝑖𝑖=𝑛𝑛1 𝑤𝑤𝑖𝑖( Δm𝑡𝑡𝑖𝑖𝑖a𝑖𝑖𝑖𝑖x𝑖𝑖𝑖𝑖(𝑖𝑖|𝑖𝑖Δ−𝑡𝑡𝑖𝑖Δ𝑖𝑖𝑖𝑖𝑖𝑡𝑖𝑡𝑖𝑖𝑖|𝑖𝑖)𝑖𝑖𝑖𝑖𝑖 )2, (16)

where Δti,obs is the experimental time delay measured at a receiver, while Δti,model is a theoretical time delay coming from the forward problem described in Section 2.1. Indeed, given a trajectory, we can run the forward
problem, compute the time delays that we should observe theoretically at the receivers and compare them with
our measurements. By minimizing the difference between Δti,obs and Δti,model at each receiver, we are actually finding the trajectory which gives the smallest difference between observed time delays and modeled ones, that
is, the “best” trajectory given our set of inputs.

Depending on the geometry as well as on the characteristics of the receiving stations, some echoes will be stronger while others will be fainter. To account for these differences, we introduce a weight wi, which captures the uncertainty that we have on the measured time delays. In our case, we weigh each time delay by the signal-to-noise ratio (SNR) obtained at the receiver, that is, as the echo is stronger, we assume that the time delay measurement will be more reliable.

For Method 2, we use the following target objective:

(

)2

 𝐽𝐽2 = 𝐽𝐽1 + ⎛⎜⎜ 𝜀𝜀 − atan √(𝑋𝑋𝐻𝐻 −𝑥𝑍𝑥𝑍𝐻𝐻𝐻𝐻)−2𝑧+𝑧𝐻(𝐻𝑌𝑌𝐻𝐻 −𝑦𝑦𝐻𝐻 )2 ⎞⎟⎟ + ⎛⎜ 𝜃𝜃 − atan( 𝑋𝑌𝑋𝑌𝐻𝐻𝐻𝐻−−𝑦𝑥𝑦𝑥𝐻𝐻𝐻𝐻 ) ⎞⎟2, (17)

⎜

𝛿𝛿𝛿𝛿

⎟⎜

𝛿𝛿𝛿𝛿

⎟

⎜

⎟⎜

⎟

⎜

⎟⎝

⎠

⎝

⎠

where J1 is computed from Equation 16, ɛ and θ are respectively the elevation and azimuth angle coming from the interferometer, which are computed as described in Jones et al. (1998). PH = (XH; YH; ZH) are the coordinates of Humain specular point and RH = (xH; yH; zH) are the coordinates of the receiver in Humain. δɛ and δθ are the uncertainties associated to the interferometric data.

3.2.  Solver Validation
To validate the trajectory reconstruction solver, 10 synthetic trajectories are simulated. They are drawn from optical data given by the CAMS-BeNeLux network (Jenniskens et al., 2011; Roggemans et al., 2016). Data were provided for 2 clear consecutive nights from 29 to 31 July 2020, in a period without any strong activity from meteor showers. At this date, the BRAMS network had only 28 active receivers. Among the 948 available

6 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

Table 1 CAMS Trajectories for the Solver Validation

No. XH (km) YH (km) ZH (km) Vx (km/s) Vy (km/s) Vz (km/s)

79

44.33 59.11 94.90 −24.59

31.22 −12.70

105 121.59 95.16 99.39 −18.98

34.62 −12.95

188 −52.88 23.12 88.13

0.55

29.19 −5.70

282 −96.09 35.92 88.28 −2.86

36.80 −16.88

477 −77.35 11.34 104.90 −35.37 −39.62 −30.68

532 38.67 51.64 91.67 −26.51

31.40 −11.89

536 21.62 199.76 105.64 −58.06

13.41 −26.14

598

6.18 158.68 103.13 −70.09

−4.60

−5.27

709 −32.21 108.29 97.79 −37.75

23.15 −45.79

773 −92.91 61.60 98.18 −37.61

14.15 −52.09

Note. No. indicates the trajectory identification number from the CAMS data. XH, YH, and ZH are the coordinates of the specular reflection point for the Humain station in a Cartesian referential centered on the Dourbes transmitter. X is directed East-West and counted positive toward East, Y is directed NorthSouth and counted positive toward North. Vx, Vy, and Vz are the velocity components of the meteoroid determined by CAMS and projected in this reference frame.

trajectories, most of them are not located above Belgium and therefore not geometrically suitable to be detected by our BRAMS receiving stations. A selection was made based on the following criteria: (a) each  trajectory must be detected by at least 6 stations, (b) because we want to compare Methods 1 & 2, one of these stations must be the interferometer in Humain, and (c) we restrict ourselves as much as possible to underdense meteor echoes in order to ensure that the specularity condition is valid. The application of these criteria resulted in a selection of 10 suitable CAMS trajectories. The parameters of these trajectories are summarized in Table 1.
The idea behind those simulated trajectories is to verify the proper functioning of the solver when there are no measurement errors, that is, checking that we recover the targeted trajectory assuming we have a perfect knowledge of the time delays at the receivers. To do so, the known CAMS trajectory is used to get the theoretical time delays, through the forward problem described in Section 2.1. The latter are then fed inside the inverse solver, and the output trajectory is compared with the CAMS one.
As shown in Figure 3, the errors are of the same order of magnitude on the 10 reconstructed trajectories: about 0.001° on the inclination angle, a few meters on the reflection point location and about 1 m/s on the velocity. This analysis demonstrates the proper functioning of our solver as it yields the correct solution if the inputs are free of measurement errors.

3.3.  Uncertainty Quantification
As discussed for the Canadian Meteor Orbit Radar (CMOR) network by Mazur et al. (2020), the inverse trajectory reconstruction problem is severely ill-conditioned when no interferometric information is used, all the receiving stations lie close to each other (less than 20 km between two receivers) and are almost on the same plane. The latter is true since the differences in altitude between the receivers are only of a few tens of meters while the trajectories are located at about 100 km in altitude. In this configuration, the output trajectories are highly sensitive to the quality of the input time delays, small errors on the latter leading to large uncertainties on the trajectory parameters.
To verify if this statement holds true with our forward scatter CW problem, we perform a Monte-Carlo uncertainty quantification campaign. For each simulated trajectory, it consists in building a Gaussian distribution of

BALIS ET AL.

Figure 3.  Reconstruction errors in logarithmic scale for 10 synthetic trajectories with Method 1 (in blue) and Method 2 (in red). Top panel: error on the inclination angle. Middle panel: error on the reference specular point position. Bottom panel: error on the velocity.
7 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

BALIS ET AL.

Figure 4.  Uncertainty quantification results for Method 1 with 6 stations (in blue), 9 stations (in red), 12 stations (in yellow). Top panel: uncertainties on the trajectory angles (in the XY plane on the left, XZ plane in the middle, YZ plane on the right). Bottom panel: uncertainties on the inclination angle (on the left), reference specular point position (in the middle) and velocity (on the right).
time delays at a chosen number of BRAMS receivers (6, 9, or 12 stations), centered on the theoretical time delay (which is calculated for the 10 CAMS trajectories with the forward problem). Three values of standard deviation σ of these Gaussian distributions are studied: 1, 5, and 10 ms.
For each of these values, and for each of the 10 simulated trajectories, we draw 1,000 sets of time delays. Then, each of these sets is fed inside our solver, and a corresponding output trajectory is computed, with the target objective given in Equation 16. We can then gather the output solutions into distributions and check their standard deviation. The uncertainty quantification results obtained with Method 1 are shown in Figure 4.
Two main observations can be made: first, it appears that only the angle in the horizontal plane XY is retrieved with an uncertainty of less than 1°. In the vertical directions, uncertainties of about 5 ms on the input time delays already lead to inclination errors of about 5–10°. The same observation holds for the reconstruction of the specular point position and of the velocity. This illustrates the ill-conditioning of the trajectory reconstruction problem without interferometry. Another important observation is that, as the number of stations that are used for the reconstruction increases from 6 to 12, the sensitivity of the output trajectory with respect to the input time delays is decreased. This shows the benefit of having more than 6 receivers detecting the signals, as errors in the determination of the time delays at different receivers can then compensate each other.
The results of the uncertainty quantification campaign run with Method 2 are shown in Figure 5. The target objective used for these trajectory reconstructions is given by Equation 17. In this case, the uncertainties on the elevation δϵ and on the azimuth δθ of the Humain specular point are ∼1°, in agreement with Jones et al. (1998, 2005). The SNRs of the considered meteor events are high since they were also detected by cameras. However, it should be noted that the 1° accuracy deteriorates at low SNR.
As already observed for CMOR by Mazur et al. (2020), the use of the interferometer alleviates the ill-conditioning of the trajectory reconstruction problem. Indeed, it appears that the uncertainties on the inclination, on the reference specular point position and on the velocity are decreased by about an order of magnitude. As an illustration, the inclination angle with about 10 ms of uncertainty on the time delays is known within an uncertainty of less that 1° with the interferometer, while it was of about 10° with Method 1. This statement holds true whatever  the number of receivers taken into account for the reconstructions, further highlighting the benefits brought by interferometric data.
4.  Radio Data Processing
Several methods can be used for the determination of time lags between two receivers. The first possibility is to use the cross-correlation peak between the signals coming at two receivers. This approach has however two significant drawbacks: on the one hand, since the location of the specular point will be different for both receivers,
8 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

BALIS ET AL.

Figure 5.  Uncertainty quantification results for Method 2 with 6 stations (in blue), 9 stations (in red), 12 stations (in yellow). The panels are the same as in Figure 4.

the echoes coming from two receivers can look very different. One can for instance be located at an altitude where the reflection is rather underdense, while the other would be more overdense. On the other hand, the BRAMS receiving stations have different characteristics in terms of noise levels. For echoes with a low SNR, this can lead to a difference in detection points, and thus to a large error on the identification of the cross-correlation peak timing. This problem was already observed for the Advanced Meteor Orbit Radar (AMOR) network by Baggaley et al. (1994). Instead, we prefer to use a physics-based reference point on each signal, which can then be used to compute time delays.

4.1.  Cornu Spiral

To do so, it is important to remember that radio systems work well to detect reflected signals from small meteoroids. As a result of this, most of the echoes which will be of interest for us will be underdense. Thus, the received power during the formation of the trail, considering an infinite trail length and homogeneous ionization, is given by McKinley (1961) and Wislez (2006):

 𝑃𝑃 (𝑡𝑡) = 𝑃𝑃0 12 ([∫−𝑥∞𝑥(𝑡𝑡) cos 𝜋𝜋2𝜋𝜋2 𝑑𝑑𝑑𝑑]2 + [∫−𝑥∞𝑥(𝑡𝑡) sin 𝜋𝜋2𝜋𝜋2 𝑑𝑑𝑑𝑑]2) = (18)

𝑃𝑃0

1 2

(c2os

+

s2in

) ,

𝐴w𝐴here c𝐴o𝐴s and sin are the Fresnel integrals, and x is the Fresnel parameter along the spiral formed by the parametric plot of the Fresnel integrals. The exact expression of P0 can be found in McKinley (1961).
As described by Mazur et al. (2020), a p𝐴𝐴lot of sin a𝐴g𝐴ainst cos produces the Cornu spiral which allows for a graphical representation of the time-varying component of the echo amplitude A from:

√  𝐴𝐴 = (Δcos)2 + (Δsin)2, (19)

𝐴w𝐴here Δc𝐴o𝐴s and Δsin are distances in the parametric plot from the starting point (−0.5; −0.5). The Cornu spiral as well as the variations of echo amplitude A are shown in Figure 6. The specular point is identified as the point (0; 0) on the Cornu spiral (see top panel) and corresponds to x = 0. Looking at the bottom panel, the dimensionless amplitude at the specular point is 0.427.
In practice, we will never observe a smooth amplitude curve as the one shown in Figure 6. Indeed, the raw signal at the time of the meteor echo is always noisy and very often has other superimposed signals on it, such as the direct signal coming from the beacon or interfering aircraft. These affect the meteor signal by adding a spectral component and modifying its amplitude and phase. They can also cause an interference beat if they have slightly different frequencies. Therefore, a pre-processing of the data is necessary to properly exploit the echoes of the

9 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

BALIS ET AL.

Figure 6.  Cornu spiral and corresponding amplitude curve. Top panel: the spiral and its characteristic points. The green circle indicates the start of the spiral (in x = −∞), the purple diamond is the passage at the specular point (x = 0), and the red cross is the end of the spiral (in x = +∞). Bottom panel: Amplitude variation with the Fresnel parameter x. The specular point is highlighted with a purple diamond.
detected meteors. The two following sections describe the methods used to subtract the beacon signal and to filter out the noise.
4.2.  Beacon Signal Subtraction
The basic principle behind the method used to remove the beacon signal is to reconstruct and then subtract it from the full signal. As mentioned above, the transmitter emits a pure cosinusoidal wave, whose amplitude and phase vary at each receiver. The received beacon signal is:
 𝑏𝑏(𝑡𝑡) = 𝐴𝐴cos(2𝜋𝜋𝜋𝜋0𝑡𝑡 + 𝜙𝜙), (20)
where b(t) is the temporal signal, A its amplitude, f0 its frequency and ϕ its phase. However, the lack of stability of the local oscillator of the first BRAMS stations (due to temperature variations in the analog receiver) as well as distortions (due to tropospheric propagation) can affect the initial signal and change the amplitude, frequency and phase parameters.
The subtraction method that we developed is based on the analysis of the Fourier transform of the signal around the beacon frequency. Since the beacon signal is a pure cosinusoidal wave, its Fourier transform is given by:
 𝐵𝐵(𝑓𝑓 ) = 12 [𝛿𝛿(𝑓𝑓 − 𝑓𝑓0)𝑒𝑒𝑖𝑖𝑖𝑖 + 𝛿𝛿(𝑓𝑓 + 𝑓𝑓0)𝑒𝑒−𝑖𝑖𝑖𝑖], (21)
which has Hermitian-symmetry with respect to frequency. The real-valued received signal is measured over a finite number of samples and its discrete Fourier transform (DFT) is computed using the Fast Fourier Transform algorithm (FFT). Because of the inherent symmetry of the spectrum, only positive frequencies need be considered; the negative-frequency terms are redundant. Moreover, the DFT is computed over a finite-time interval, which corresponds to a multiplication of the raw signal by a rectangular time window whose DFT is a sinc function (J. O. Smith, 2011):
sin(𝜋𝜋𝜋𝜋𝜋𝜋 )  𝑊𝑊𝑟𝑟(𝑓𝑓 ) = 𝐴𝐴𝐴𝐴sinc(𝜏𝜏𝜏𝜏 ) = 𝐴𝐴 , (22)
𝜋𝜋𝜋𝜋
where τ is the duration of the finite-time interval.
Thus, the DFT of the meteor signal over a finite-time interval is the DFT of the product between the raw signal and a time window, which means the convolution product between a Dirac delta function δ and a sinc function in the spectral domain. This convolution gives a sinc function centered on the δ position. The sinc is thus centered on the received beacon frequency and has an amplitude equal to that of the received beacon signal. But because the DFT only samples the underlying Fourier transform at multiples of the sampling frequency, its maximum
10 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

BALIS ET AL.

Figure 7.  Example of a DFT on a 1 s interval and corresponding sinc fitting. Top panel: DFT on the whole frequency spectrum. Bottom panel: zoom around the beacon frequency and sinc function fitting. The blue circles are the DFT samples, the red curve is the sinc function fitted on these points and the purple star is the maximum of the fit.
does not necessarily coincide with the maximum of the sinc function. It is thus necessary to fit a sinc to the DFT points in order to retrieve the amplitude A and frequency f of the peak. This procedure is illustrated in Figure 7.
In order to compute the phase ϕ, the following formula is used:
 𝜙𝜙 = 𝜙𝜙0 + 𝛿𝛿𝛿𝛿 𝑁𝑁 − 1 − 2𝜋𝜋𝜋𝜋𝑏𝑏𝑡𝑡𝑠𝑠, (23) 𝑁𝑁
where ϕ0 is the phase at the maximum of the discrete FFT, N is the number of points on which the FFT is computed, and fb is the beacon frequency obtained with the fit of the sinc function. The second term is a correction to get the phase at the maximum of the continuous fitted sinc, where δ is the difference of frequency index between the discrete peak of the FFT and its continuous peak. The last term is a linear time correction where ts represents the timing at which the considered time interval starts.
In order to maximize the frequency resolution, the whole time interval of 300 s contained in each radio file should ideally be used to compute the FFT. However, the longer this time interval, the more likely the signal is subject to small variations of the local oscillator and/or variation of tropospheric conditions. A balance has thus to be found between instabilities and frequency resolution. Experience has shown that, in practice, short time intervals give better results even if the frequency resolution is decreased. Therefore, we decide to split each radio file in small intervals of 1 s in which the reconstruction is done independently.
Two further remarks need to be made. First, to reduce the spectral leakage characteristic of the rectangular window, we decide to use a Hann window (Harris, 1978; Nuttall, 1981) so that the function fitted in the frequency domain is:
 𝑊𝑊ℎ(𝑓𝑓 ) = 𝐴2𝐴 (1si−nc𝜏(𝜏𝜏2𝜏𝜏𝑓𝜏𝑓 )2) . (24)
Second, this whole approach works fine as long as there is no other signal superimposed on the beacon frequency, such as a meteor. To account for this possibility, we use a statistical criterion looking at the Hann window fit quality on each consecutive 1 s interval. If the fit on a given interval is more than 3 median absolute deviation (MAD) worse than the median fit quality over 300 s, the reconstruction in this interval is rejected. In this case, the amplitude, frequency and phase are interpolated from the neighboring intervals. The results obtained with this approach are displayed in Figure 8. While it appears that the subtraction without interpolation (middle panel) introduces some spectral artifact for the meteor highlighted in the red rectangles, the approach with interpolation (bottom panel) preserves its shape.
4.3.  Filtering and Smoothing
After the beacon signal has been subtracted, the noise has to be filtered to isolate the actual echo coming from the meteor. To do this, we use a windowed-sinc band-pass filter defined by the product of a sinc function and a Blackman window. This approach significantly reduces the noise level while keeping a steep transition band and an excellent stopband attenuation, as described in S. Smith (2013).
11 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

Figure 8.  Example of spectrograms before and after beacon signal subtraction. Top panel: Raw spectrogram before processing. Middle panel: Spectrogram after beacon subtraction without interpolation. Bottom panel: Spectrogram after beacon subtraction with interpolation. In the three panels, the power is color-coded in dB FS (Decibels relative to Full Scale). The red rectangles highlight the same meteor in the three spectrograms.

The band-pass filter expression is obtained from the convolution of a low-pass filter and a high-pass filter. The design of such filters requires the selection of two parameters: the cut-off frequency fc and the length of the filter kernel M (number of points of the truncated sinc). M can be computed using:
 𝑀𝑀 ≈ 4 , (25) 𝐵𝐵𝐵𝐵

where BW is the width of the transition band in the filter frequency response (magnitude of its transfer function), measured from where the curve just barely leaves one to where it almost reaches zero. To have a good tradeoff between a sharp transition band and a reasonable length of the filter kernel, BW is always chosen equal to 4 × 10 −3 so that M = 1,000 points. After the selection of fc and M, the filter kernel is calculated with:

( ( ))

 ℎ[𝑖𝑖] = 𝐾𝐾 sin 2𝜋𝜋𝜋𝜋𝑐𝑐 𝑖𝑖 − 𝑀2𝑀 [0.42 − 0.5cos( 2𝜋𝜋𝜋𝜋 ) + 0.08cos( 4𝜋𝜋𝜋𝜋 )].(26)

𝑖𝑖 − 𝑀𝑀

𝑀𝑀

𝑀𝑀

2

Three components can be identified: a constant K chosen to provide unity

gain

at

the

maximum

of

the

window,

the

sinc

function

𝐴w𝐴ith

a

𝑀𝑀 2

shift

to

get

only positive index values and the Blackman window equation.

Figure 9.  Example of amplitude curves of a meteor echo for two windowedsinc filters. The blue and red curves are respectively obtained with a filter passband of 600 and 60 Hz.

Most of the meteor echoes detected are located in a range of 100 Hz around the beacon frequency fb. However, if a too narrow filter is used, the dynamics of the filtered meteor echo will be modified too strongly by the impulse response of the filter. This is illustrated in Figure 9, where the rise of the echo is made unphysically longer by a bandpass of 60 Hz. Heuristically, we find that a bandpass of 600 Hz, that is, fhigh = fb − 300 Hz and flow = fb + 300 Hz is a good trade-off between noise reduction and preservation of the signal dynamics.
The very good frequency performance of the windowed-sinc filter contrast with its poor characteristics in the time domain, notably in terms of ripple and overshoot in the step response. Moreover, as we filter the signal with a bandpass of 600 Hz, it sometimes happen that another parasitic signal in this frequency band remains unfiltered. As a result, oscillations with the meteor

BALIS ET AL.

12 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

echo are created and deteriorate the identification of the timing corresponding to the specular point.

Therefore, a Savitzky-Golay (SG) filter is applied to the data points (Savitzky & Golay, 1964). The idea is the following: for each point in the signal, a polynomial is fitted to that point and the neighboring ones. The value of the polynomial is evaluated at that point and taken as the value of the filtered signal. The filtering operation can be reduced to the calculation of a convolution product of the signal by a kernel adapted to the user's choice: polynomial degree on the one hand, and number of neighboring points used for the filter (i.e., window length) on the other hand.

Figure 10.  Example of amplitude curves of a meteor echo for two Savitzky-Golay windows. The blue curve is obtained without smoothing, while the green and magenta curves are obtained with a smoothing window of respectively 301 and 1,501 points. The purple diamond is the specular point. The dashed black line is the median noise level. The dashed-dotted red line is the median + 3 MAD level.

Heuristically, we choose a polynomial order of 3 and a window length of 301 points. The reason for this choice is that a too small window does not dampen the oscillations enough, while a too large one tends to oversmooth the rise of the meteor echo. The output amplitude curve obtained with two different SG filters for a chosen echo is shown in Figure 10. The specular point (corresponding to a dimensionless amplitude of 0.427) is highlighted on the amplitude curve with a window length of 301 samples. The median and median + 3 MAD noise levels are computed on an interval of 2 s before the start of the meteor echo.

The exact timing (in UTC reference) of the specular point can be identified thanks to the GPS timestamps. A further correction is applied to account for the delays of the different types of receivers. Each has its own band-pass filter with its own response time that needs to be compensated for, in order to compare the timings at several stations. The procedure to obtain these delays will be explained in a future paper describing in detail the BRAMS network.

Finally, the solver described in Section 3 is called with the time delays as inputs, and outputs a reconstructed trajectory. A summary of the complete workflow is given in Figure 11.

5.  Comparison With Optical Trajectories
In this section, we compare reconstructions obtained with the BRAMS data and the optical trajectories given by the CAMS-BeNeLux network, for the 10 cases presented in Table 1. At the time of writing, the CAMS network in the Benelux consists of about 100 cameras that cover the sky above Belgium, The Netherlands, Luxembourg and part of Germany. In the case where several cameras record the same meteor, its trajectory and speed are calculated with a very good accuracy, so that we can use them as a reference. Indeed, CAMS can compute radiants with an uncertainty of less than 2° and speeds with an uncertainty of less than 10% (Jenniskens et al., 2016). Further information about the network can be found in Roggemans et al. (2016).
For each trajectory, the CAMS data file contains the following parameters:
• The time interval during which the meteor was visible by certain cameras. • The geographical position (longitude, latitude and altitude) of the start and end points of the visible trajectory. • The total velocity V∞ corresponding to the appearance of the meteor (associated entry into the atmosphere). • The acceleration parameters a1 and a2 modeled by CAMS (Jacchia et al., 1967). • The uncertainties (standard deviation) on the positions, velocity and acceleration parameters.
CAMS use a deceleration model introduced by Whipple and Jacchia (1957) and Jacchia et al. (1967):
 𝑉𝑉 (𝑡𝑡) = 𝑉𝑉∞ − |𝑎𝑎1𝑎𝑎2|exp(𝑎𝑎2𝑡𝑡), (27)

Figure 11.  Workflow for the trajectory reconstruction.

where t = 0 corresponds to the timing of the first recorded trail. Since we use a constant speed assumption in our reconstruction method, the comparison between our velocity output and the velocity given by CAMS is not

BALIS ET AL.

13 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

Figure 12.  Example of result obtained for CAMS trajectory 79 using Method 1: (a) Projected view in the horizontal plane. (b) 3D view. The green line is the CAMS trajectory. The blue cross is the reference station (Humain). The red plus is the transmitter (Dourbes). The light green diamonds are the stations used for the BRAMS reconstruction. The orange and light blue line represent respectively the radio wave path from the transmitter to Humain specular point, and the radio wave path from this specular point to the receiver in Humain. The purple line is the reconstructed trajectory with the BRAMS signals. The blue line represents the Belgium border.
straightforward. Nonetheless, the beginning points detected by CAMS are around the same altitude as the specular points detected by BRAMS (80−120 km), so that we decide to compare the BRAMS constant speed with the velocity output by CAMS in t = 0.
An example of trajectory reconstruction with Method 1 is shown in Figure 12. This plot presents the projected CAMS trajectory 79 as well as the reconstructed one in the horizontal XY plane and in a 3D frame, where coordinates X, Y, and Z are given in a local Cartesian frame centered on the Dourbes transmitter. Seventeen BRAMS receivers are used to perform the trajectory reconstruction. Despite that the location of Humain specular point is off by about 22 km, the velocity is accurate within 6% of the CAMS data (39.65 km/s instead of 41.72 km/s), and the inclination is off by less than 0.6°.
Figure 13 presents the reconstruction results obtained for trajectory 79, but using Method 2. Since the direction of the reflection point is constrained via equations including the interferometer data, it is now better retrieved and

BALIS ET AL.

14 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

Figure 13.  Example of result obtained for CAMS trajectory 79 using Method 2. The panels, lines and markers are the same as in Figure 12.
this helps greatly the reconstruction of the trajectory. The altitude of the reflection point becomes much more accurate with an error of less than 2 km. The velocity direction is accurate to 0.2° and the speed (41.07 km/s) is very close to the 41.72 km/s measured by CAMS.
Table 2 compares the number of BRAMS receivers detecting exploitable echoes for each of the 10 CAMS trajectories, and gives the reconstructed angles in the horizontal plane and the inclination angles obtained with Method 1 and 2. As it could be expected from the uncertainty quantification campaign of Section 3.3, Method 1 is not able to retrieve accurate vertical inclinations and positions systematically, due to its inherent ill-conditioned nature. It offers nonetheless good information on the orientation in the horizontal plane.
Also in agreement with our uncertainty quantification results, the horizontal angles reconstructions are within the 2° uncertainty of CAMS with both methods 1 and 2, except for meteoroids which are fast and whose trajectories are inclined with respect to the horizontal plane (in particular trajectory 773 which combines a high-speed meteoroid with a large inclination angle). This phenomenon can be explained because, as the meteoroid enters the atmosphere faster, it spends less time in the observable region of BRAMS (i.e., altitudes roughly between 80 and 120 km). Thus, the time delays between the exploitable specular points tend to be smaller. As a result, errors

BALIS ET AL.

15 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

Table 2 Comparison of BRAMS Reconstructions With Methods 1 and 2 and CAMS

on the identification of the specular point location at each receiver lead to larger errors on the measured time delays. As the trajectory is further inclined,

Trajectories

this phenomenon is amplified. In this case, the distance between the specular

No.

Echoes

ΔαXY,1 (°)

ΔαXY,2 (°)

Δα1 (°) Δα2 (°) points (and thus the time delays for a given speed) is indeed smaller than if the

trajectory is more horizontal. Similarly, high speeds and large inclinations also

79

17

0.56

0.09

1.45

0.16

lead to high deceleration rates, making the constant speed model less accurate.

105

13

1.70

1.44

8.69

2.02

188

6

0.13

0.09

7.49

0.32

This observation is further reinforced by looking at the values of Δα1. Indeed, the only cases for which the inclination angle reconstruction is accurate

282

6

0.43

1.47

37.19

1.93

enough with Method 1 are trajectories 79 and 105. These two meteors have

477

17

0.53

0.75

24.95

1.43

either a higher number of exploitable echoes or a lower speed compared to

532

8

0.18

0.11

16.77 0.52 the trajectories with larger reconstruction errors. The trajectory inclination,

536

8

2.25

0.81

22.99

1.46

the meteor velocity and the number of echoes with a sufficient SNR are thus

598

19

2.44

1.29

8.79

3.32

important variables which condition the reachable accuracy of the trajectory

709

11

2.67

2.13

3.38

1.48

reconstruction.

773

12

4.21

2.89

20.67

1.76

Note. ΔαXY,i is the difference between the angles in the horizontal plane reconstructed by CAMS and the method i using BRAMS data. Δαi is the difference between the inclination angles reconstructed by CAMS and the

Table 3 summarizes the reconstruction results obtained with Method 2 on the 10 CAMS trajectories. We can see that the velocity retrieved by our algorithm is always within the 10% error margin of the one given by CAMS. Regarding the inclination, we are always within the 2° accuracy range of

method i using BRAMS data.

CAMS, except for trajectory 598, which is the meteor with the largest veloc-

ity. For the altitude of the specular point, the results are in a lesser agreement

as they are contained within ∼10 km of the value given by CAMS (except for

trajectory 536, which is at 21.91 km). It is interesting to note a similar trend as before: trajectories with a high

speed and inclination (e.g., trajectories 536, 709, and 773) are particularly difficult to reconstruct.

The results given in Table 3 have been obtained with Equation 17. If one has better knowledge about the quality of the data at their disposal (for instance, very good interferometric data compared to less accurate time delays measurements), one could think of giving more weight to the elevation and azimuth terms of this target objective. This might improve the altitude determination of the specular points (since the interferometer gives the direction of one specular point), but this could also deteriorate the velocity reconstruction (as it is mainly driven by the contribution of the time delays).

Table 3 Comparison of BRAMS Reconstructions With Method 2 and CAMS Trajectories

No. ϵ (°)

VCAMS (km/s)

VBRAMS (km/s)

−

[%] 

ΔZH (km)

Δα (°)

79 17.72 105 18.16 188 11.04 282 24.58 477 30.01 532 16.14 536 23.69 598 4.29 709 45.96 773 52.35

41.72 41.55 29.74 40.59 61.33 42.79 65.07 70.44 63.70 65.78

41.07 40.60 28.39 41.32 62.56 41.41 70.13 71.79 65.96 71.49

−1.58 −2.34 −4.75 1.77 1.97 −3.33 7.21 1.88 3.43 7.99

1.98 0.16 7.79 2.02 4.66 0.32 6.75 1.93 10.24 1.43 7.04 0.52 21.91 1.46 8.42 3.32 8.58 1.48 0.39 1.76

Note. The second and third columns give the complement of the zenith angle ϵ and the speed of the meteoroid retrieved by CAMS at their first detection point VCAMS. The fourth column gives the speed of the meteoroid retrieved by BRAMS VBRAMS. The relative errors in terms of velocity, Humain specular point altitude as well as inclination angle (also given in Table 2) are shown in the last three columns.

Let us finally note that this comparison with optical data raises a question regarding the nature of the echo reflections. Indeed, optical systems tend to observe meteors of larger size than radio systems, since they need to be bright enough (the limiting magnitude of CAMS is +5 according to Jenniskens et al. (2021)) in order to be detected. As a result, the selected BRAMS echoes corresponding to the 10 CAMS trajectories are sometimes overdense or transitional. In that case, the specularity often holds as the amplitude rise of the meteor echo is not affected, but creation of secondary reflection points sometimes appear due to the distortion of the meteor trail through wind shear (Wislez, 2006).
6.  Conclusion and Perspectives
In this work, we solved for the first time the general CW forward scatter trajectory reconstruction problem through a nonlinear optimization approach. We studied how the uncertainties on the time delays between meteor echoes at several radio stations impact the reconstruction accuracy of the meteoroid trajectory and speed. Method 1 (using the time delays only) appeared to be strongly ill-conditioned while interferometric data brought by Method 2 alleviated this very high sensitivity to the inputs. We implemented the data processing workflow from the acquisition of the raw signals at the receivers to the reconstructed trajectory. We finally compared our results with data coming from the optical CAMS-BeNeLux network.

BALIS ET AL.

16 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

Since Method 1 is the only approach that can be applied to all archived BRAMS data, its development and improvement are crucial. However, because of its ill-conditioning, it cannot be used systematically as it gives large errors in terms of specular point position and inclination angle on most trajectories. Therefore, we are currently extending the pre-t0 phase approach described for backscatter radars in Mazur et al. (2020). This technique uses the phase of the signal before the specular point in order to retrieve the velocity of the meteoroid. It provides two important benefits. On the one hand, it further constrains the trajectory solver by adding extra inputs, thus reducing the sensitivity of the output trajectory on the input time delays. On the other hand, it allows to get some insight on the deceleration of the meteoroid during its upper atmosphere travel.
More generally, deceleration models will be introduced in the solver described in this paper, even without the extra information coming from the pre-t0 phase. Depending on the distribution of time delays, a constant velocity can provide good fits for a meteoroid with slight deceleration by using a velocity in the middle of the region observed by the BRAMS receivers, rather than the start velocity. Other possibilities such as linear and exponential velocity models will also be tested.
We have shown that Method 2 gives more accurate results. Unfortunately, it is not applicable to the majority of past data since it requires to have at least 3 stations detecting the same meteor as the interferometer, which does not systematically happen because of the specularity condition. We have also demonstrated that increasing the number of receiving stations detecting the same meteoroid decreases the impact of measurement errors on the quality of the reconstruction. Therefore, two priorities are to build a second interferometer in the north of Belgium (where there is already a dense cluster of BRAMS receiving stations, see Figure 1) and to create clusters of radio stations around the existing interferometric station in Humain.
In parallel, the BRAMS network continues to develop. A technical description will be published shortly, also tracing its evolution over time. This future paper will also describe the trade-offs between BRAMS and other existing systems, notably in terms of area coverage, cost and data quality. In addition to the creation of receiving stations clusters and the building of a new interferometer, the use of phase coding on the transmitted CW signal (Vierinen et al., 2016) would in principle allow us to obtain the range or total distance traveled by the wave between the transmitter and the receivers. It has been shown (Lamy et al., 2021) that this extension of Method 1 would greatly facilitate the trajectory reconstruction problem.

Data Availability Statement
Version 1.1 of the BRAMS trajectories package, used for the reconstruction of meteoroid trajectories with the time delays at the BRAMS receiving stations, is preserved at https://doi.org/10.5281/zenodo.7686085, available via Creative Commons Attribution 4.0 International and developed openly at https://gitlab.aeronomie.be/ae/ brams/bramstrajectories (Balis, 2023).

Acknowledgments
The main author is funded by the Solar-Terrestrial Center of Excellence (STCE) (https://stce.be). The BRAMS network is also partly funded by STCE. The BRAMS project is an active Pro-Am collaboration. We would like to thank all operators hosting the receiving stations. We would also like to thank Carl Johannink and Peter Jenniskens for providing us with CAMS data.

References
Anciaux, M., Lamy, H., Martinez Picar, A., Calders, S., Calegaro, A., Ranvier, S., & Verbeeck, C. (2020). The BRAMS receiving station v2.0. In U. Pajer, J. Rendtel, M. Gyssens, & C. Verbeeck (Eds.), Proceedings of the International Meteor Conference, Bollmannsruh, Germany, 3–6 October 2019 (pp. 39–42).
Baggaley, W. J., Bennett, R. G. T., Steel, D. I., & Taylor, A. D. (1994). The Advanced Meteor Orbit Radar facility—AMOR. The Quarterly Journal of the Royal Astronomical Society, 35, 293.
Balis, J. (2023). BRAMS trajectories package (version 1.1) [Software]. Zenodo. https://doi.org/10.5281/zenodo.7686085 Boyd, S., & Vandenberghe, L. (2004). Convex optimization. Cambridge University Press. https://doi.org/10.1017/CBO9780511804441 Cevolani, G., & Pupillo, G. (2009). Ground-based radio observations to probe the ozone content in the meteor region. Annals of Geophysics,
46(2). https://doi.org/10.4401/ag-3399 Chau, J. L., Urco, J. M., Vierinen, J. P., Volz, R. A., Clahsen, M., Pfeffer, N., & Trautner, J. (2019). Novel specular meteor radar systems using
coherent MIMO techniques to study the mesosphere and lower thermosphere. Atmospheric Measurement Techniques, 12(4), 2113–2127. https://doi.org/10.5194/amt-12-2113-2019 Froncisz, M., Brown, P., & Weryk, R. J. (2020). Possible interstellar meteoroids detected by the Canadian Meteor Orbit Radar. Planetary and Space Science, 190, 104980. https://doi.org/10.1016/j.pss.2020.104980 Hajduková, M., Jr., Kornoš, L., & Tóth, J. (2014). Frequency of hyperbolic and interstellar meteoroids. Meteoritics & Planetary Sciences, 49(1), 63–68. https://doi.org/10.1111/maps.12119 Harris, F. (1978). On the use of windows for harmonic analysis with the discrete Fourier transform. Proceedings of the IEEE, 66(1), 51–83. https://doi.org/10.1109/PROC.1978.10837 Jacchia, L., Verniani, F., & Briggs, R. E. (1967). An analysis of the atmospheric trajectories of 413 precisely reduced photographic meteors. Smithsonian Contributions to Astrophysics, 10, 1–139. https://doi.org/10.5479/si.00810231.10-1.1

BALIS ET AL.

17 of 18

1944799x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023RS007697, Wiley Online Library on [20/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Radio Science

10.1029/2023RS007697

Jacobs, E., & Ralston, E. W. (1981). Ambiguity resolution in interferometry. IEEE Transactions on Aerospace and Electronic Systems, AES, 17(6), 766–780. https://doi.org/10.1109/TAES.1981.309127
Jeanne, S. (2020). Méthode d’analyse statistique appliquée au reseau d’observation européen des météores FRIPON (Theses). Université Paris Sciences et Lettres.
Jenniskens, P., Gural, P., Dynneson, L., Grigsby, B., Newman, K., Borden, M., et al. (2011). CAMS: Cameras for allsky meteor surveillance to establish minor meteor showers. Icarus, 216(1), 40–61. https://doi.org/10.1016/j.icarus.2011.08.012
Jenniskens, P., Lauretta, D. S., Towner, M. C., Heathcote, S., Jehin, E., Hanke, T., et al. (2021). Meteor showers from known long-period comets. Icarus, 365, 114469. https://doi.org/10.1016/j.icarus.2021.114469
Jenniskens, P., Nénon, Q., Albers, J., Gural, P., Haberman, B., Holman, D., et al. (2016). The established meteor showers as observed by CAMS. Icarus, 266, 331–354. https://doi.org/10.1016/j.icarus.2015.09.013
Jones, J., Brown, P., Ellis, K., Webster, A., Campbell-Brown, M., Krzemenski, Z., & Weryk, R. (2005). The Canadian Meteor Orbit Radar: System overview and preliminary results. Planetary and Space Science, 53(4), 413–421. https://doi.org/10.1016/j.pss.2004.11.002
Jones, J., Webster, A. R., & Hocking, W. K. (1998). An improved interferometer design for use with meteor radars. Radio Science, 33(1), 55–65. https://doi.org/10.1029/97RS03050
Karmarkar, N. K., & Ramakrishnan, K. G. (1991). Computational results of an interior point algorithm for large scale linear programming. Mathematical Programming, 52(1–3), 555–586. https://doi.org/10.1007/bf01582905
Kastinen, D., & Kero, J. (2022). Radar analysis algorithm for determining meteor head echo parameter probability distributions. Monthly Notices of the Royal Astronomical Society, 517(3), 3974–3992. https://doi.org/10.1093/mnras/stac2727
Lamy, H., Anciaux, M., Ranvier, S., Calders, S., Gamby, E., Martinez Picar, A., & Verbeeck, C. (2015). Recent advances in the BRAMS network. In Proceedings of the International Meteor Conference, Mistelbach, Austria, 27–30 August 2015 (pp. 171–175).
Lamy, H., Balis, J., & Anciaux, M. (2021). Reconstructing meteoroid trajectories using BRAMS data. WGN, Journal of the International Meteor Organization, 49(6), 195–200.
Lamy, H., & Tétard, C. (2016). Retrieving meteoroids trajectories using BRAMS data: Preliminary simulations. In A. Roggemans, & P. Roggemans (Eds.), Proceedings of the International Meteor Conference, Egmond, the Netherlands, 2–5 June 2016 (pp. 143–147).
Lamy, H., Tétard, C., Anciaux, M., Ranvier, S., Picar, A. M., Calders, S., & Verbeeck, C. (2018). First observations with the BRAMS radio interferometer. In M. Gyssens, & J.-L. Rault (Eds.), Proceedings of the International Meteor Conference, Petnica, Serbia, 21–24 September 2017 (pp. 132–137).
Markkanen, J., Lehtinen, M., & Landgraf, M. (2005). Real-time space debris monitoring with EISCAT. Advances in Space Research, 35(7), 1197–1209. (Space Debris). https://doi.org/10.1016/j.asr.2005.03.038
Mazur, M., Pokorný, P., Brown, P., Weryk, R. J., Vida, D., Schult, C., et al. (2020). Precision measurements of radar transverse scattering speeds from meteor phase characteristics. Radio Science, 55(10), e2019RS006987. https://doi.org/10.1029/2019RS006987
McKinley, D. (1961). Meteor science and engineering. McGraw-Hill. Retrieved from https://books.google.be/books?id=BSAIAQAAIAAJ Nedeljković, S. (2006). Meteor forward scattering at multiple frequencies. In L. Bastiaens, J. Verbert, J.-M. Wislez, & C. Verbeeck (Eds.),
Proceedings of the International Meteor Conference, Oostmalle, Belgium, 15–18 September 2005 (pp. 16–24). Nuttall, A. (1981). Some windows with very good sidelobe performance. IEEE Transactions on Acoustics, Speech, and Signal Processing, 29(1),
84–91. https://doi.org/10.1109/TASSP.1981.1163506 Roggemans, P., Johannink, C., & Breukers, M. (2016). Status of the CAMS-BeNeLux network. In Proceedings of the International Meteor
Conference, Egmond, the Netherlands, 2–5 June 2016 (pp. 254–260). Savitzky, A., & Golay, M. J. E. (1964). Smoothing and differentiation of data by simplified least squares procedures. Analytical Chemistry, 36(8),
1627–1639. https://doi.org/10.1021/ac60214a047 Schult, C., Stober, G., Janches, D., & Chau, J. L. (2017). Results of the first continuous meteor head echo survey at polar latitudes. Icarus, 297,
1–13. https://doi.org/10.1016/j.icarus.2017.06.019 Smith, J. O. (2011). Spectral audio signal processing. W3K. Smith, S. (2013). Chapter 16—Windowed-sinc filters. In Digital signal processing: A practical guide for engineers and scientists (pp. 285–296).
Newnes. Stober, G., & Chau, J. L. (2015). A multistatic and multifrequency novel approach for specular meteor radars to improve wind measurements in
the MLT region. Radio Science, 50(5), 431–442. https://doi.org/10.1002/2014RS005591 Vierinen, J., Chau, J. L., Pfeffer, N., Clahsen, M., & Stober, G. (2016). Coded continuous wave meteor radar. Atmospheric Measurement Tech-
niques, 9(2), 829–839. https://doi.org/10.5194/amt-9-829-2016 Whipple, F. L., & Jacchia, L. G. (1957). Reduction methods for photographic meteor trails. Smithsonian Contributions to Astrophysics, 1,
183–206. Wislez, J.-M. (2006). Meteor astronomy using a forward scatter set-up. In C. Verbeeck & J.-M. Wislez (Eds.), Proceedings of the Radio Meteor
School, Oostmalle, Belgium, 10–14 September 2005 (pp. 84–106).

BALIS ET AL.

18 of 18