Articles | Volume 14, issue 5
Research article
31 May 2021
Research article |  | 31 May 2021

Error analyses of a multistatic meteor radar system to obtain a three-dimensional spatial-resolution distribution

Wei Zhong, Xianghui Xue, Wen Yi, Iain M. Reid, Tingdi Chen, and Xiankang Dou

In recent years, the concept of multistatic meteor radar systems has attracted the attention of the atmospheric radar community, focusing on the mesosphere and lower thermosphere (MLT) region. Recently, there have been some notable experiments using such multistatic meteor radar systems. Good spatial resolution is vital for meteor radars because nearly all parameter inversion processes rely on the accurate location of the meteor trail specular point. It is timely then for a careful discussion focused on the error distribution of multistatic meteor radar systems. In this study, we discuss the measurement errors that affect the spatial resolution and obtain the spatial-resolution distribution in three-dimensional space for the first time. The spatial-resolution distribution can both help design a multistatic meteor radar system and improve the performance of existing radar systems. Moreover, the spatial-resolution distribution allows the accuracy of retrieved parameters such as the wind field to be determined.

1 Introduction

The mesosphere and lower thermosphere (MLT) is a transition region from the neutral to the partially ionized atmosphere. It is dominated by the effects of atmospheric waves, including planetary waves, tides and gravity waves. It is also a relatively poorly sampled part of the Earth's atmosphere by ground-based instruments. One widely used approach to sample this region is the meteor radar technique. The ablation of incoming meteors in the MLT region, i.e.,  80–110 km, creates layers of metal atoms, which can be observed from the ground by photometry or lidar (Jia et al., 2016; Xue et al., 2013). During meteor ablation, the trails caused by small meteor particles provide a strong atmospheric tracer within the MLT region that can be continuously detected by meteor radars regardless of weather conditions. Consequently, the meteor radar technique has been a powerful tool for studying the MLT region for decades (Hocking et al., 2001; Holdsworth et al., 2004; Jacobi et al., 2008; Stober et al., 2013; Yi et al., 2018). Most modern meteor radars are monostatic, and this has two main limitations in retrieving the complete wind fields. Firstly, limited meteor rates and relatively low measurement accuracies necessitate that all measurements in the same height range are processed to calculate a “mean” wind. Secondly, classic monostatic radars retrieve winds based on the assumption of a homogenous wind in the horizontal and usually zero wind in the vertical direction.

The latter conditions can be partly relaxed if the count rates are high and the detections are distributed through a representative range of azimuths. If this is the case, a version of a velocity–azimuth display (VAD) analysis can be applied by expanding the zonal and meridional winds using a truncated Taylor expansion (Browning and Wexler, 1968). This is because each valid meteor detection yields a radial velocity in a particular viewing direction of the radar. The radar is effectively a multi-beam Doppler radar where the “beams” are determined by the meteor detections. If there are enough suitably distributed detections in the azimuth in a given observing period, the Taylor expansion approach (using Cartesian coordinates) yields the mean zonal and meridional wind components (u0,v0), the horizontal divergence ux+vy, and the stretching ux-vy and shearing uy+vx deformations of the wind fields from an analysis of the radial velocities. However, because the radar can only retrieve the wind projection in the radial direction as measured from the radar, the vorticity vx-uy of the wind fields is not available. This is common to all monostatic radar systems, and a discussion of measurable parameters in the context of multiple fixed-beam upper-atmosphere Doppler radars is given by (Reid, 1987). Even when relaxing the assumption of a homogeneous wind fields and using the more advanced volume velocity processing (VVP) (Philippe and Corbin, 1979) to retrieve the wind fields, the horizontal gradients of the wind fields cannot be recovered due to the lack of vorticity information. To obtain a better understanding of the spatial variation of the MLT region wind fields, larger area observations (and hence higher meteor count rates) and sampling of the observed area from different viewing angles are needed. An extension of the classic monostatic meteor technique is required to satisfy these needs.

To resolve the limitations outlined above, the concept of multistatic meteor radar systems, such as MMARIA (multistatic and multifrequency agile radar for investigations of the atmosphere) (Stober and Chau, 2015) and SIMO (single-input multiple-output) (Spargo et al., 2019), MIMO (multiple-input multiple-output radar) (Chau et al., 2019; Dorey et al., 1984) have been designed and implemented (Stober et al., 2018). Multistatic systems can utilize the forward scatter of meteor trails, thus providing another perspective for observing the MLT. Multistatic meteor radar systems have several advantages over classic monostatic meteor radars, such as obtaining higher-order wind field information and covering wider observation areas. There have been some particularly innovative studies using multistatic meteor radar systems in recent years. For example, by combining MMARIA and the continuous wave multistatic radar technique (Vierinen et al., 2016), Stober et al. (2018) built a five-station, seven-link multistatic radar network covering an approximately 600 km× 600 km observing region over Germany to retrieve an arbitrary non-homogenous wind field with a 30 km× 30 km horizontal resolution. Chau et al. (2017) used two adjacent classic monostatic specular meteor radars in northern Norway to obtain horizontal divergence and vorticity. Other approaches, such as coded continuous-wave meteor radar (Vierinen et al., 2019) and the compressed sense method in MIMO sparse-signal recovery (Urco et al., 2019) are described in the corresponding references.

Analyzing spatial resolution limits is a fundamental but difficult topic for meteor radar systems. Meteor radar systems transmit and then receive radio waves reflected from meteor trails using a cluster of receiving antennas, commonly five antennas, as in the Jones et al. (1998) configuration. By analyzing the cross-correlations of the signals received on several pairs of antenna, the angle of arrival (AoA) of each return can be determined. The AoA is described by the zenith angle θ and azimuth angle ϕ. By measuring the wave propagation time from the meteor trail, range information can be determined. Most meteor radar systems rely on specular reflections from meteor trails. Thus, by combining the AoA and the range information and then using geometric analysis, the location of a meteor trail can be determined. Accurately locating the meteor trail specular point (MTSP hereinafter) is important since atmospheric parameter retrieval (such as the wind field or the temperature) depends on the location information of meteor trails. The location accuracy, namely the spatial resolution, determines the reliability of the retrieved parameters. For multistatic meteor radar systems that can relax the assumption of a homogenous horizontal wind field, the location accuracy becomes a more important issue because the horizontal spatial resolution affects the accuracy of the retrieved horizontal wind field.

Although meteor radar systems have developed well experimentally in recent years, the reliability of the retrieved atmospheric parameters still requires further investigation for both the monostatic and multistatic meteor radar cases. In an attempt to investigate errors in two radar techniques, Wilhelm et al. (2017) compared 11 years of MLT region wind data from a partial reflection (PR) radar with colocated monostatic meteor radar winds and determined the “correction factors” required to bring the winds into agreement. Spargo et al. (2019) reported a similar study for two locations for data obtained over several years. While the comparisons are interesting, partial reflection radars operating in the medium-frequency (MF) band and lower high-frequency (HF) band produce a height-dependent bias in the measured winds (see, e.g., Reid, 2015), which limits the ability to estimate errors in the meteor winds by comparing them. However, the PR radar technique is one of very few that provides day and night coverage and data rates in the MLT comparable to that of meteor radars.

Meteor radars have largely replaced PR radars for MLT studies and are generally regarded as providing reference-quality winds. It is essential then to know the reliability of atmospheric parameters determined by meteor radars, and to do this some quantitative error analyses are necessary.

Figure 1Schematic diagram of the simplified bistatic configuration used in Hocking's vertical resolution analysis (Hocking, 2018). The two receiving antennas and the transmitting antenna are collinear. The analysis is in a two-dimensional vertical section through the baseline joining the antennas. The radio wave is scattered from a few Fresnel zones that are several kilometers in length around the specular point on the meteor trail and are received by the receiving antennas. The cross-correlation analysis between the receiving antennas can be used to solve for the AoAs. Because the radio wave is reflected from a region a few Fresnel zones in length, the measured phase difference between the receiver antenna pairs to deviates from the ideal phase difference. This deviation from the ideal phase difference is one of the error sources in the PDME. In this work, we solve for the ideal phase difference associated with the AoA directed to the MTSP.


A number of recent studies have discussed AoA measurement errors for meteor radars (Kang, 2008; Vaudrin et al., 2018; Younger and Reid, 2017). These studies focus on the phase errors in receiver antenna pairs; see Younger and Reid (2017) for the monostatic case and Vaudrin et al. (2018) for a more general case that included multistatic meteor radars. Hocking (2018) used another approach and developed a vertical resolution analysis method for the two-dimensional bistatic case. The Hocking method (HM hereinafter) simplifies the error propagation process in the receiving antennas and puts emphasis on how a bistatic meteor radar configuration affects the vertical resolution in a vertical section. It does not consider the radial distance measuring error. In this paper, we consider the more general three-dimensional case and determine the spatial distribution of both the horizontal and vertical resolution uncertainties.

We analyze the multistatic meteor radar resolution distribution in a three-dimensional space for both vertical and horizontal resolution for the first time. This spatial resolution is a prerequisite for evaluating the reliability of retrieved atmospheric parameters, such as the wind field and the temperature.

2 Analytical method

2.1 Brief introduction

The HM will be introduced briefly here to help the reader understand our generalization. In the HM, measurement errors that affect the vertical resolution can be classified into two types: one caused by the zenith angle measuring error δθ and one caused by the pulse length effect on the vertical resolution. The receiving array is a simple antenna pair that is collinear with the baseline (Fig. 1). The HM calculates the vertical resolution in a two-dimensional vertical section that passes through the baseline. The receiver antenna pair is equivalent to one receiver arm in a Jones configuration that is comprised of three collinear antennas usually with a 2λ\2.5λ spacing. The phase difference of the received radio wave between the receiving antenna pair is denoted as ΔΨ. In meteor radar systems, there is generally an “acceptable” phase difference measuring error (PDME hereinafter) δ(ΔΨ). A higher value of δ(ΔΨ) means that more detected signals will be judged as meteor events but with more misidentifications and bigger errors as well. δ(ΔΨ) is set to approximately 30 (Hocking, 2018; Younger and Reid, 2017) in most meteor radar systems. In the HM, the zenith angle measuring error δθ is due to δ(ΔΨ) and δ(ΔΨ) is a constant. Therefore, the error propagation in the receiver is very simple, and δθ is inversely proportional to the cosine of the zenith angle.

Figure 2Schematic diagram of a multistatic meteor radar system using SIMO (single-input and multi-output). There are three transmitters (T1, T2 and T3) and one receiver (RX) in the picture. The transmitter–receiver distance is typically 100–200 km. X0, Y0 and Z0 represent the east, north and up directions of the receiving antenna. Over 90 % of the received energy comes from about 1 km around the specular point of the meteor trail, which is slightly less than the length of the central Fresnel zone (Ceplecha et al., 1998).


We now introduce our analytical method. Our method considers a multistatic system with multiple transmitters and one receiving array in three-dimensional space as shown in Fig. 2. The receiving array is in the Jones configuration, that is, “cross-shaped”, but it may also be “T-shaped” or “L-shaped”. The five receiver antennas are in the same horizontal plane and constitute two orthogonal antenna arms. To avoid a complex error propagation process in the receiving array and to place emphasis on multistatic configurations, the PDMEs in the two orthogonal antenna arms (δ(ΔΨ1) and δ(ΔΨ2)) are constants. Therefore, the AoA measuring errors (including the zenith and azimuth angle measuring errors δθ, δϕ, respectively) can be expressed as simple functions of zenith and azimuth angle. The radial distance is the distance between the MTSP and the receiver, which is denoted as Rs. Rs can be determined by combining the AoA, baseline length di and the radio wave propagation path length R (Stober and Chau, 2015). The geometry is shown in Fig. 4a. α is the angle between the baseline (i.e.,Xi axis) and the line from the receiver to the MTSP (denoted as point A). If α, di and R are known, Rs can be calculated easily using the cosine law as follows:

(1) R s = R 2 - d i 2 2 ( R - d i cos α ) .

A multistatic configuration will influence the accuracy of Rs (denoted as δRs). This is because α, d and R are determined by the multistatic configuration. We consider the error term δRs in our method, which is ignored in the HM. δRs is a function of the AoA measuring errors (δθ and δϕ) and the radio wave propagation path length measuring error (denoted as δR). δR is caused by the measuring error of the wave propagation time δt, which is approximately 21 µs (Kang, 2008). Thus, δR can be set as a constant and the default value in our program is δR=cδt=6.3km. It is worth noting that the maximum unambiguous range for pulse meteor radars is determined by the pulse repetition frequency (PRF) (Hocking et al., 2001; Holdsworth et al., 2004). For multistatic meteor radars utilizing forward scatter, the maximum unambiguous range is c/PRF (where c is the speed of light). For the area where R exceeds the maximum unambiguous range, δR is set to positive infinity.

Figure 3(a) Schematic diagram of the three coordinate systems used in this work. XiYiZi is a class of coordinate systems whose Xi axis points to transmitter i, with i=1,2,3. X0Y0Z0 is the ENU coordinate system to which all errors are compared. (b) Magnified plot of the receiving array. XYZ is fixed on the receiver horizontal plane. The x and y axes are collinear with the two arms of the antenna array.


Figure 4(a) Schematic diagram of the forward-scatter geometry for the radar link between Ti and RX. Point A is the MTSP. (b) Magnified plot of specular point A. The red line represents a radio wave pulse, and S is the half-pulse length. kb is the Bragg vector, which halves the forward-scatter angle β. (c) Schematic diagram of E1 in XYZ, which can be decomposed into three orthogonal vectors.


2.2 Three kinds of coordinate systems and their transformations

To better depict the multistatic system configuration, three kinds of right-hand coordinate systems need to be established, as shown in Fig. 3. These are X0Y0Z0, XiYiZi and XYZ. X0Y0Z0 is the ENU (east–north–up) coordinate system where the X0, Y0 and Z0 axes represent the east, north and up directions, respectively. Another two coordinate systems are established to facilitate different error propagations. All types of errors need to be transformed to the ENU coordinate system X0Y0Z0 in the end. Coordinate system XYZ is established to depict the spatial configuration of the receiving array and has its the origin of XYZ there, as shown in Fig. 3. The z axis is collinear with the antenna boresight and perpendicular to the horizontal plane on which the receiving array lies. The x and y axes are collinear with the arms of the two orthogonal antenna arrays. AoAs will be represented in XYZ for convenience. Inspection of Fig. 4 indicates that it is convenient to analyze the range information in a plane that goes through the baseline and the MTSP. Thus, a coordinate system XiYiZi is established for a transmitter Ti. The coordinate origins of XiYiZi are all on the receiving array. We stipulate that the Xi axis points to transmitter i (Ti). Each pair Ti and receiver RX constitute a radar link, which is referred to as Li. The range-related information for each Li will be calculated in XiYiZi. Different types of errors need to propagate to and be compared in X0Y0Z0, which is convenient for retrieving wind fields.

We stipulate that clockwise rotation satisfies the right-hand corkscrew rule. By rotating ψxX,i, ψyY,i and ψzZ,i about the x, y and z axes, respectively, one can transform XYZ to XiYiZi. It is worth mentioning that XiYiZi is non-unique because any rotation about Xi axis can obtain another satisfactory XiYiZi. Hence, ψxX,i can be set to any value. Similarly, by rotating ψxi,0, ψyi,0 and ψzi,0 about the x, y and z axes, respectively, one can transform XiYiZi to X0Y0Z0. To realize the coordinate transformation between these three coordinate systems, a coordinate rotation matrix AR(ψx,ψy,ψz) is introduced. Using AR, one can transform the coordinate point or vector presentation from one coordinate system to another. The details of the coordinate rotation matrix AR(ψx,ψy,ψz) can be found in Appendix A.

2.3 Two types of measuring errors

The analytical method of the spatial resolution for each radar link is the same. The difference between these radar links is only the value of the six coordinate rotation angles (ψxX,i, ψyY,i and ψzZ,i;ψxi,0, ψyi,0 and ψzi,0) and the baseline distance di. The spatial-resolution-related measurement errors that will cause location errors of the MTSP can be classified into two types: E1 is caused by measurement errors at the receiver, and E2 is due to the pulse length. These two errors are mutually independent. Hence, the total error (Etotal) can be expressed as follows:

(2) E total 2 = E 1 2 + E 2 2 .

E1 is related to three indirect measuring errors. They are zenith, azimuth and radial distance measuring errors, denoted as δθ, δϕ and δRs, respectively. In XYZ, E1 can be decomposed into three orthogonal error vectors using δθ, δϕ and δRs (see Fig. 4c), which we now explain in more detail. PDMEs, i.e., δ(ΔΨ1) and δ(ΔΨ2), are caused by some practical factors, such as phase calibration mismatch and the fact that the specular point is not actually a point but is a few Fresnel zones in length. A meteor radar system calculates phase differences between different pairs of antenna through cross-correlation and then fits them to get the most likely AoAs. Therefore, the system needs to be assigned a tolerance value of δ(ΔΨ1) and δ(ΔΨ2). Different meteor radar systems have different AoA fit algorithms and thus different AoA measuring error distributions. To analyze the spatial resolution for a SIMO meteor radar system as generally as possible and to avoid tedious error propagation at the receiving array, we start the error propagation from δ(ΔΨ1) and δ(ΔΨ2) and set them as constants. AoA measuring errors δθ and δϕ can then be expressed as follows:


where λ is the radio wavelength, D1 and D2 are the length of the two orthogonal antenna arms, and θ and ϕ are the zenith angle and the azimuth angle, respectively. The details can be found in Appendix A. It is worth noting that δθ and δϕ are not mutually independent. The expectation value of their product is not identical to zero unless Eδ2(ΔΨ1)D12 is equal to Eδ2(ΔΨ2)D22. δRs can be expressed as a function of δR, δθ and δϕ as follows:

(5) δ R s = F ( δ R , δ θ , δ ϕ ) = f R ( θ , ϕ ) δ R + f θ ( θ , ϕ ) δ θ + f ϕ ( θ , ϕ ) δ ϕ .

fR(θ,ϕ),fθ(θ,ϕ) and fϕ(θ,ϕ) are the weighting functions of δRs. The details about the weighting function and deduction can be found in Appendix A. Inspection of Fig. 4c indicates that E1 can be decomposed into three orthogonal error vectors in coordinate XYZ, denoted as δRs, Rsδθ and Rssin θδϕ. These three vectors can be expressed in XYZ as follows:


E2 is related to the radio wave propagation path. A pulse might be reflected anywhere within a pulse length (see Fig. 4b). This causes a location error in the MTSP, represented as an error vector DA. D is the median point of the isosceles triangle ΔABC's side BC. The representation of the error vector DA can be solved in XiYiZi by using geometrical relationships as follows:

(9) D A = ( ( 2 - a 1 - a 2 ) x i + d i ( a 2 - 1 ) 2 , ( 2 - a 1 - a 2 ) y i 2 , ( 2 - a 1 - a 2 ) z i 2 ) T ,

where S is the half-pulse length and a1=Rs-SRs. a2=Ri-SRi. di is the baseline length. (xi,yi,zi) is the coordinate value of a MTSP (point A in Fig. 4) in XiYiZi. More details can be found in Appendix A.

2.4 Transformation to ENU coordinates

Thus far, two types of errors in different coordinate systems have been introduced. Now they need to be transformed to ENU coordinates X0Y0Z0 in order to compare different radar links and to analyze the wind fields. E1-related error vectors, which are three orthogonal vectors δRs, Rsδθ, and Rssin θδϕ, are represented in XYZ as Eqs. (6)–(8) and need to be transformed from XYZ to X0Y0Z0 to project δRs, Rsδθ and Rssin θδϕ towards the X0Y0Z0 axis, respectively, and reassemble them to form three new error vectors in the X0Y0Z0 axis. Using the coordinate rotation matrix AR(XYZ,X0Y0Z0)=AR(ψxi,0,ψyi,0,ψzi,0)AR(ψxX,i,ψyY,i,ψzZ,i) and Eqs. (6)–(8), the unit vectors of those three vectors can be represented in X0Y0Z0 as follows:

(10) X 0 ( δ R s ) X 0 ( δ θ ) X 0 ( δ ϕ ) Y 0 ( δ R s ) Y 0 ( δ θ ) Y 0 ( δ ϕ ) Z 0 ( δ R s ) Z 0 ( δ θ ) Z 0 ( δ ϕ ) = A R ( X Y Z , X 0 Y 0 Z 0 ) sin θ cos ϕ cos θ cos ϕ - sin ϕ sin θ sin ϕ cos θ sin ϕ cos ϕ cos θ - sin θ 0 .

X0(δRs),Y0(δRs),Z0(δRs)T, X0(δθ),Y0(δθ),Z0(δθ)T, and X0(δϕ),Y0(δϕ),Z0(δϕ)T are unit vectors of δRs, Rsδθ, and Rssin θδϕ in X0Y0Z0, respectively. The 3×3 matrix on the left-hand side of Eq. (10) is denoted as Pij for i,j=1,2,3.

From Eqs. (6)–(8) and Fig. 4c, we see that the length of those three vectors (the error values) is δRs, Rsδθ, Rssin θδϕ as a function of δR, δθ, δϕ. In order to reassemble them to form new error vectors, transformation of δθ and δϕ into two independent errors δ(ΔΨ1) and δ(ΔΨ2) is needed because δθ and δϕ are not independent. Using Eqs. (3) and (4), one can transform vector δR,δθ,δϕT into three independent measuring errors δR, δ(ΔΨ1) and δ(ΔΨ2). Thus, δRs,Rsδθ,RssinθδϕT can be expressed as follows:

(11) δ R s R s δ θ R s sin θ δ ϕ = f R ( θ , ϕ ) f θ ( θ , ϕ ) f ϕ ( θ , ϕ ) 0 R s 0 0 0 R s sin θ 1 0 0 0 λ 2 π cos θ cos ϕ D 1 - λ 2 π sin ϕ sin θ D 1 0 λ 2 π cos θ sin ϕ D 2 λ 2 π cos ϕ sin θ D 2 δ R δ Δ Ψ 1 δ Δ Ψ 2 .

The product of the first and the second term on the right-hand side of Eq. (11) is a 3×3 matrix, denoted as Wij for i,j=1,2,3. From Eq. (11), we see that the three error values δRs, Rsδθ and Rssin θδϕ are the linear combinations of δR, δ(ΔΨ1),and δ(ΔΨ2) with their corresponding linear coefficients W1j, W2j and W3j. Those three error values can be projected toward new directions (e.g., the X0Y0Z0 axis) by using Pij. It worth noting that in a new direction the same basis's projected linear coefficients from different error values should be used to calculate their sum of squares (SS). Following this, the square root of SS will be used as a new linear coefficient for that basis in the new direction. For example, in X0 directions, basis δ(ΔΨ1)'s projected linear coefficients are X0(δRs)W12, X0(δθ)W22 and X0(δϕ)W32 from δRs, Rsδθ and Rssin θδϕ, respectively. Therefore, the new linear coefficient for δ(ΔΨ1) in the X0 direction is WX0δ(ΔΨ1)=±X0(δRs)W122+X0(δθ)W222+X0(δϕ)W322. Similarly, one can get δR and δ(ΔΨ2)'s new linear coefficients in X0, denoted as WX0δR and WX0δ(ΔΨ2). Thus, the true error value in the X0 direction is WX0δRδR+WX0δ(ΔΨ1)δ(ΔΨ1)+WX0δ(ΔΨ2)δ(ΔΨ2). Because δR, δ(ΔΨ1), and δ(ΔΨ2) are mutually independent, E1 is related to the mean squared error (MSE) values in the X0 direction, denoted as δ(1)X0 and can be expressed as δ(1)X0=± WX0δRδR2+WX0δ(ΔΨ1)δ(ΔΨ1)2+WX0δ(ΔΨ2)δ(ΔΨ2)2.

In short, E1-related errors in ENU coordinate's three axis directions (denoted as δ(1)X0, δ(1)Y0 and δ(1)Z0) can be expressed in the form of a matrix as follows:

(12) δ ( 1 ) 2 X 0 δ ( 1 ) 2 Y 0 δ ( 1 ) 2 Z 0 = P i j 2 W i j 2 δ 2 R δ 2 ( Δ Ψ 1 ) δ 2 Δ Ψ 2 .

The E2-related error vector DA needs transformation from XiYiZi to X0Y0Z0. Therefore, E2-related errors in the ENU coordinate's three axis directions (denoted as δ(2)X0, δ(2)Y0 and δ(2)Z0) can be expressed in the form of a matrix as follows:

(13) δ ( 2 ) X 0 δ ( 2 ) Y 0 δ ( 2 ) Z 0 = ± A R ( ψ x i , 0 , ψ y i , 0 , ψ z i , 0 ) D A .

E1 and E2 are mutually independent. By using Eq. (1), the total MSE values in ENU coordinate's three axis directions (denoted as δtotalX0, δtotalY0 and δtotalZ0) can be expressed in the form of matrix as follows:

(14) δ total 2 X 0 δ total 2 Y 0 δ total 2 Z 0 = δ ( 1 ) 2 X 0 δ ( 1 ) 2 Y 0 δ ( 1 ) 2 Z 0 + δ ( 2 ) 2 X 0 δ ( 2 ) 2 Y 0 δ ( 2 ) 2 Z 0 .

Figure 5(a) The flow chart of the location error calculation process for a point in X0Y0Z0. The notation beside the arrows represents the corresponding equations (black) or coordinate rotation matrix (blue) in the paper.  is the Hadamard product. Thus, E2E2 will yield (δ(2)2X0,δ(2)2Y0,δ(2)2Z0)T. (b) The flow chart of the program to calculate the location errors distributions for a radar link Li. This process includes parameters settings for a radar link, the generation of the sampling grid nodes and the traversing of all the nodes. For each node, the program uses the calculation method described in (a). MC is the multistatic configuration, and IC is the interferometer (receiving antenna) configuration. (c) Schematic diagram of the relationship between the spatial resolution and the total location errors of the MTSP. For a detected point in space, the MSE of MTSP's location errors is ±|δtotalX0|, ±|δtotalY0| and ±|δtotalZ0| in the zonal, meridional and vertical directions, respectively. This means that the actual specular point might occur in a region that forms a 2|δtotalX0|×2|δtotalY0|×2|δtotalZ0| cube and that the detected point is on the centroid of this cube.


In conclusion, for a radar link Li and a MTSP represented as (x0,y0,z0) in the ENU coordinate system X0Y0Z0, as sketched in Fig. 4a, the location errors of this point in east, north and up directions (±δtotalX0, ±δtotalY0 and ±δtotalZ0) can be calculated as follows: firstly, for a point (x0,y0,z0) in X0Y0Z0, use AR to transform it to XiYiZi and denote it as (xi,yi,zi). Following this, in XiYiZi calculate the AoA (θ and ϕ) and the range information (Rs and Ri). Details of AoA and range calculation can be found in Appendix A. It is worth noting that the AoA is given by the angles relative to the axes of XYZ. Secondly, in XYZ use the AoA and Eqs. (3)–(8) to calculate E1's three orthogonal error vectors shown in Fig. 4c; in XiYiZi use the range information and Eq. (9) to calculate E2's error vector DA, as shown in Fig. 4b. Thirdly, project E1's three error vectors to X0Y0Z0 by using Eq. (10) and use Eqs. (11)–(12) to reassemble them to calculate E1-related MSE values in the direction of X0Y0Z0; use Eq. (13) to transform the E2 error vector from XiYiZi to X0Y0Z0. Finally, use Eq. (14) to get the total location errors of a MTSP in (x0,y0,z0). Figure 5a shows the flow chart for the process we have just described.

3 Results and discussion

The program to study the method we have described above is written in the Python language and is presented in the Supplement. To calculate a special configuration of a multistatic radar system, we initially need to set six coordinate transformation angles (ψxX,i, ψyY,i and ψzZ,i;ψxi,0, ψyi,0 and ψzi,0) and the baseline length di for each radar link Li. For example, ψxi,0=ψyi,0=0, ψzi,0=30 and di=250 km means that transmitter Ti is 250 km, 30 east by south of the receiver RX. Further, ψxX,i=5, ψyY,i=0, ψzZ,i=0 means one receiver arm (y axis) points to 60 east by north with a 5 elevation. The detection area of interest for a multistatic meteor radar is usually from 70 to 110 km in height and around 300 km× 300 km in the horizontal. In our program, this area needs to be divided into a spatial grid for sampling. The default value of the sampling grid length is 1 km in height and 5 km in the meridional and zonal directions. After selecting the desired settings, the program steps though the sampling grid nodes and calculates the location errors at each node as described in Fig. 5a. Figure 5b describes the parameter settings and the transversal calculation process. For a given setting of radar link Li, the program will output the squared values of E1-related, E2-related and total MSE values (Etotal2: δtotal2X0,δtotal2Y0,δtotal2Z0; E12: δ(1)2X0, δ(1)2Y0,δ(1)2Z0; E22: δ(2)2X0,δ(2)2Y0,δ(2)2Z0). The location errors can be positive or negative, and thus the spatial resolutions are twice the absolute value of the location errors. For an example, see Fig. 5c. For a detected MTSP, represented as (x0,y0,z0) in X0Y0Z0, with δtotal2X0,δtotal2Y0 and δtotal2Z0 equal to 25, 16 and 9 km2, respectively, the actual position of the MTSP could occur in an area that is ± 5, ± 4 and ± 3 km around (x0,y0,z0) with equal probability. Consequently, the zonal, meridional and vertical resolutions are 10, 8 and 6 km, respectively.

Figure 6The normalized vertical resolution distribution in a vertical section from 50 to 120 km height when the error term “δRs” is ignored. Panels (a–c) are the total, Rsδθ-related and E2-related normalized resolution distributions, respectively. These results are the same as those produced in Hocking's work (Hocking, 2018). The two black arrows represent the positions right above the transmitter (Tx) and the receiver (Rx) and the transmitter–receiver separation is 300 km. The region between the two oblique black lines is the sampling volume for the receiving array because the elevation angle is beyond 30 to reduce influence from potential mutual antenna coupling or from other obstacles in the surrounding area. Except the region at large elevation angles (i.e., 90), the E2-related resolution values are much lower than the Rsδθ-related errors. The Rsδθ-related resolution distribution only depends on the receiving antennas. Thus, the total vertical resolution distribution is nearly unchanged with the variation of the transmitter–receiver distance. Normalized resolution values that exceed 3 km (which correspond 12 km vertical resolution) are not shown.


The HM analyzes the vertical resolution (corresponding to δZ0 in our paper) in a two-dimensional vertical section (corresponding to the X0Z0 plane in our paper). To compare with Hocking's work, ψzi,0 is set to 180, and the other five coordinate transformation angles are all set to zero with d equal to 300 km. The half-pulse length S is set to 2 km and δ(ΔΨ1) to 35. Calculating in the X0Y0 plane only should have degraded our method into Hocking's two-dimensional analysis method but does not because the HM method ignores δRs. In fact, the HM considers only E2 and Rsδθ in the X0Y0 plane. Consequently, we need to further set fR(θ,ϕ), fθ(θ,ϕ) and fϕ(θ,ϕ) to be zero. When this is done, our method degrades into the HM. Hocking's results are shown as the absolute value of vertical location error normalized relative to the half-pulse width |δZ0|/S. Hereinafter, |E|/S is referred to as the normalized spatial resolution, such as δ(1)X0 and δtotalY0, where E represents the location errors in a direction. Thus, the spatial resolutions are 2S times the normalized spatial resolutions.

Our normalized vertical resolution distributions are shown in Fig. 6a and are the same as those presented in Hocking's work (Hocking, 2018). The distribution of Rsδθ-related, E2-related and total normalized vertical resolution distributions are shown in Fig. 6 from left to right, respectively. In most cases, E2 is an order of magnitude smaller than Rsδθ. Only in the region directly above the receiver does E2 have the same magnitude as Rsδθ. In other words, only in the region directly above the receiver can E2 influence the total resolution. E2 is related to the bistatic configuration but Rsδθ is not. Therefore, in the HM, the distribution of the total vertical resolution varies slightly with d. After adding the error term δRt, which is related to the bistatic configuration, the normalized total vertical spatial-resolution distribution changes more obviously with d, as the first two rows in Fig. 7 show. The region between the two black lines represents the sampling volume for the receiver where the elevation angle is beyond 30. As the transmitter–receiver distance becomes longer, resolutions in this sampling volume are not always acceptable. In the first row of Fig. 7, the transmitter–receiver distance is 300 km and about half of the region between the two black lines has normalized vertical resolution values larger than 3 km. Because our analytical method can obtain spatial resolutions in three-dimensional space, the third row of Fig. 7 shows the horizontal section at 90 km altitude for the second row of Fig. 7.

Figure 7The normalized vertical resolution distribution using the analytical method described in this paper. The first and second rows represent a vertical section of height from 50 to 120 km. The third row represents the horizontal section at 90 km, and the receiving array is at the origin, with positive coordinate values representing the eastward or northward directions. The first row has the same parameters settings as Fig. 6 and is used to compare with Fig. 6. The E1-related resolution will change with the transmitter–receiver configuration because it considers the error term “δRs”. Thus, the total vertical resolution will change with the transmitter–receiver configuration. With the transmitter–receiver distance varying from 300 km (the first row) to 150 km (the second row), the total vertical resolution distribution is clearly changed. The third row is the perspective to the horizontal section at 90 km altitude for the second row. Normalized resolution values that exceed 3 km are not shown.


To get a perspective on the spatial-resolution distribution in three-dimensional space, Fig. 8 shows the normalized zonal, meridional and vertical spatial-resolution distributions for a multistatic radar link whose transmitter–receiver separation is 180 km and the transmitter is 30 south by east of the receiver. The classic monostatic meteor radar is a special case of a multistatic meteor radar system whose the baseline length is zero. By setting the transmitter–receiver distance to be zero in our program, a monostatic meteor radar's spatial resolution can also be obtained. In this case, the spatial-resolution distributions are highly symmetrical and correspond to the real characteristics of monostatic meteor radar (this is not shown here but can be found in Fig. S1 in the Supplement). In the discussion above, the receiver and transmitter antennas are all coplanar. By varying ψxX,i, ψyY,i and ψzZ,i in our program, non-coplanar receiver–transmitter antenna situations can also be studied. Slightly tilting the receiver horizontal plane (for example, set ψxX,i=ψyY,i=5) causes the horizontal spatial distributions to change (see Figs. S2 and S3 in the Supplement). In practice, the Earth's curvature and local topography will lead to tilts in the receiver horizontal plane. This kind of tilt should be taken into account for multistatic meteor radar systems, and details relating to the parameter selections for this can be found in the Supplement.

Figure 8The 3D contour plot of the normalized resolution distribution for a multistatic radar link whose baseline length is 180 km and whose transmitter is 30 south by east of the receiver. The black dots represent the position right above the transmitter, and the receiving array is at the origin of the axes. Panels (a–c) are the normalized resolution distributions in the zonal, meridional and vertical directions, respectively. The subplot's four slice circles from bottom to top are the horizontal section at 50, 70, 90 and 110 km height, respectively. The region whose elevation angle of the receiver is less than 30 is not shown, and therefore the slice circles become larger from the bottom to the top. Normalized resolution values that exceed 4 km (which corresponds to 16 km resolution) are not shown.


The AoA error propagation process has been simplified to yield Eqs. (3) and (4) by using constant PDMEs. This is for the sake of providing the most general example of our method. If the analysis of AoA errors were to start from the original received voltage signals (e.g., Vaudrin et al., 2018), the error propagation process would depend on the specific receiver interferometer configuration and the specific signal processing method. The approach used here can be applied to different receiver antenna configurations or new signal processing algorithms. This would involve substitution of δ(ΔΨ1) and δ(ΔΨ2) into other mutually independent measuring errors to suit the experimental arrangement and then establishing a new AoA error propagation to obtain δθ and δϕ. This means rewriting the second and third term in Eq. (11) to the determine a new AoA error propagation matrix and new mutually independent measuring errors, respectively.

It worth noting that except for using the PDMEs as the start of the error propagation, all the analytical processes are built on mathematical error propagations. PDMEs include uncontrolled errors, such as those resulting from the returned wave being scattered from a few Fresnel zones along the meteor trail, phase calibration inaccuracy and noise. However, there are other error sources in practice. For example, aircraft or lightning interference and fading clutter from obstacles can cause further measurement errors in the AoA. These issues are related to actual physical situations and are beyond the scope of this text.

Knowing the valid observational volume for meteor detections and the errors associated with each detection is vital for a meteor radar system as it determines which meteors can be used to calculate wind velocities and also the uncertainties associated with the winds themselves. To reduce the influence of mutual antenna coupling or ground clutter, the elevation angle of a detection should be above a threshold, with 30 being typically used, and this sets the basic valid observational volume. Within this, the normalized vertical resolution varies, and in Figs. 7 and S4 in the Supplement only the areas of normalized vertical resolution with values below 3 km are shown, which we argue represents an acceptable sampling volume. In addition, as the transmitter–receiver distance increases, the sampling volume becomes smaller and the vertical resolution in this volume is reduced. This effect limits the practically usable transmitter–receiver distances for multistatic meteor radars.

The geometry of the multistatic meteor radar case also impacts on the ability of the radar to measure the Doppler shifts associated with drifting meteor trails within the observational volume. This is because the measured Doppler shift is produced by the component of the wind field in the direction of the Bragg vector, which in the multistatic configuration is divergent from the receiver's line of sight (see e.g., Spargo et al., 2019). The smaller the angle between the Bragg vector and the wind fields is, the larger the Doppler shift (and the higher the SNR) will be. This means that within the observational volume, the angular diversity of the Bragg vector should be taken into account in the wind retrieval process. A discussion of wind retrievals is beyond the scope of this text and will be considered in future work.

4 Conclusions

In this study, we have presented preliminary results from an analytical error method analysis of multistatic meteor radar system measurements of angles of arrival. The method can calculate the spatial resolution (the spatial uncertainty) in the zonal, meridional and vertical directions for an arbitrary receiving antenna array configuration in three-dimensional space. A given detected MTSP is located within the spatial-resolution volume with an equal probability. Higher values of spatial-resolution mean that this region needs more meteor counts or longer averaging to obtain a reliable accuracy. Our method shows that the spatial configuration of a multistatic system will greatly influence the spatial-resolution distribution in ENU coordinates and thus will in turn influence the retrieval accuracy of atmospheric parameters such as the wind field. The multistatic meteor radar system's spatial-resolution analysis is a key point in analyzing the accuracy of retrieved wind and other parameters. The influence of the spatial resolution on wind retrieval will be discussed in future work.

Multistatic radar systems come in many types, and the work in this paper considers only single-input (single-antenna transmitter) and multi-output (five-antenna interferometric receiver) pulse radar systems. Although the single-input multi-output (SIMO) pulse meteor radar is a classic meteor radar system, other meteor radar systems, such as continuous-wave radar systems and MISO (multiple-antenna transmitter and single-antenna receiver), also show good experimental results. Using different types of meteor radar systems to constitute a meteor radar network is the future trend, and so we will add the spatial-resolution analysis of other system types using our method in the future. We will also validate and apply the spatial-resolution analysis in the horizontal wind determination to a multistatic meteor radar system that will be soon be installed in China.

Appendix A

A1 Coordinates rotation matrix

For a right-handed rectangular coordinate system XYZ, we rotate clockwise Ψx about the x axis to obtain a new coordinate. We specify that clockwise rotation satisfies the right-hand screw rule. A vector in XYZ, denoted as (x,y,z)T, is represented as (x,y,z)T in the new coordinate. The relationship between (x,y,z)T and (x,y,z)T is as follows:

(A1) x y z = A x ( ψ x ) x y z = 1 0 0 0 cos ψ x sin ψ x 0 - sin ψ x cos ψ x x y z .

Similarly, we rotate Ψy clockwise about the y axis to obtain a new coordinate. The presentation for a vector in new coordinates and the original can be linked by a matrix, Ay(ψy):

(A2) A y ( ψ y ) = cos ψ y 0 - sin ψ y 0 1 0 sin ψ y 0 cos ψ y ,

we rotate Ψz clockwise about the z axis to obtain a new coordinate. The presentation for a vector in new coordinates and original can be linked by a matrix Az(ψz):

(A3) A z ( ψ z ) = cos ψ z sin ψ z 0 - sin ψ z cos ψ z 0 0 0 1 .

For any two-coordinate systems XYZ and XYZ with co-origin, one can always rotate Ψx, Ψy and ψz clockwise in the order of the x, y and z axes, respectively, transforming XYZ to XYZ (Fig. A1). The presentation for a vector in XYZ and XYZ can be linked by a matrix, AR(ψx,ψy,ψz).

(A4) A R ( ψ x , ψ y , ψ z ) = A z ( ψ z ) A y ( ψ y ) A x ( ψ x ) = cos ψ y cos ψ z sin ψ x sin ψ y cos ψ z + cos ψ x sin ψ z - cos ψ y sin ψ z - sin ψ x sin ψ y sin ψ z + cos ψ x cos ψ z sin ψ y - sin ψ x cos ψ y - cos ψ x sin ψ y cos ψ z + sin ψ x sin ψ z cos ψ x sin ψ y sin ψ z + sin ψ x cos ψ z cos ψ x cos ψ y

Figure A1Schematic diagram of Coordinates rotation.


Figure A2The receiving array geometry (only three antennas are shown for clarity).


We call AR(ψx,ψy,ψz) the coordinate rotation matrix.

A2 AoA measuring errors

In coordinate XYZ, AoAs includes zenith angle θ and azimuth angle ϕ. In the plane wave approximation, the radio wave is at angle γ1 and γ2 with an antenna array (Fig. A2). There is a phase difference ΔΨ1 and ΔΨ2 between two antennas (Fig. 1). ΔΨ1 and ΔΨ2 can be expressed as follows:


Using γ1, γ2 the AoA can be expressed as follows:


or in another expression:


substitute cos γ1 and cos γ2 in Eqs. (A7) and (A8) by using Eqs. (A5) and (A6):


Equations (A11) and (A12) link the phase difference with the AoA and expanding θ and ϕ, ΔΨ1 and ΔΨ2 to the first order.


For Eqs. (A13) and (A14), substitute ΔΨ1 and ΔΨ2 using Eqs. (A5), (A6) and (A9), (A10) to the functions of θ, ϕ. Now Eqs. (3) and (4) have been proven. If the zenith angle θ=0, we stipulate that cosϕsinθ and sinϕsinθ are 1.

A3 Radial distance measuring error

Expanding Rs, R and cos α in Eq. (1) to the first order, δRs can be expressed as a function of δR and δ(cos α):

(A15) δ R s = R 2 - 2 R d cos α + d 2 2 ( R - d cos α ) 2 δ R + d ( R 2 - d 2 ) 2 ( R - d cos α ) 2 δ ( cos α ) ,

where α is the angle between Rs and the Xi axis. We denote the zenith and azimuth angles in coordinate XiYiZi as θ and ϕ, respectively. The relationship between α and θ, i.e., ϕ, is

(A16) cos α = sin θ cos ϕ .

Using the coordinate rotation matrix AR(ψxX,i,ψyY,i,ψzZ,i) and sinθcosϕ can be expressed as the function of AoA:

(A17) sin θ cos ϕ = A 11 sin θ cos ϕ + A 12 sin θ sin ϕ + A 13 cos θ .

Aij represent the elements in matrix AR(ψxX,i,ψyY,i,ψzZ,i) for i,j=1,2,3.

Using Eqs. (A16) and (A17), δ(cos α) can be expressed as a function of δθ and δϕ as follows:

(A18) δ cos α = ( A 11 cos θ cos ϕ + A 12 cos θ sin ϕ - A 13 sin θ ) δ θ + ( - A 11 sin θ sin ϕ + A 12 sin θ cos ϕ ) δ ϕ .

Finally, δRs can be expressed as the function of δR,δθ,δϕ as follows.

(A19) δ R s = F ( δ R , δ θ , δ ϕ ) = f R ( θ , ϕ ) δ R + f θ ( θ , ϕ ) δ θ + f ϕ ( θ , ϕ ) δ ϕ
(A20) f R ( θ , ϕ ) = d 2 + R 2 - 2 R d ( A 11 sin θ cos ϕ + A 12 sin θ sin ϕ + A 13 cos θ ) 2 R - d ( A 11 sin θ cos ϕ + A 12 sin θ sin ϕ + A 13 cos θ ) 2
(A21) f θ ( θ , ϕ ) = d ( R 2 - d 2 ) ( A 11 cos θ cos ϕ + A 12 cos θ sin ϕ - A 13 sin θ ) 2 R - d ( A 11 sin θ cos ϕ + A 12 sin θ sin ϕ + A 13 cos θ ) 2
(A22) f ϕ ( θ , ϕ ) = d ( R 2 - d 2 ) ( - A 11 sin θ sin ϕ + A 12 sin θ cos ϕ ) 2 R - d ( A 11 sin θ cos ϕ + A 12 sin θ sin ϕ + A 13 cos θ ) 2

A4 True error of E2

The total length of side AC and side AB represents the pulse width. Side AC equals side CB and they are both equal to half of the pulse width S. In XiYiZi, the presentation of point A is (xi,yi,zi), the receiver is (0,0,0) and Ti is (d,0,0). The distance between Ti and A is Ri=R-Rs. We denote the presentation of point B and C in XiYiZi as (xB,yB,zB) and (xC,yC,zC), respectively. We use vector collinear to establish equations for B and C. Therefore, one can obtain the coordinates of point B and C by the following equations:

(A23) ( x B , y B , z B ) T = R s - S R s ( x i , y i , z i ) T ,

(A24) ( x C - d , y C , z C ) T = R i - S R i ( x i - d , y i , z i ) T .

For isosceles triangle ABC, the perpendicular line AD intersects side CB at the midpoint D. Thus, we obtain the coordinate value of D in XiYiZi as follows:

(A25) ( x D , y D , z D ) = 1 2 ( x B + x c , y B + y c , z b + z c ) = 1 2 ( ( a 1 + a 2 ) x i - a 2 d + d , ( a 1 + a 2 ) , y i ( a 1 + a 2 ) z i ) .

We denote a1=Rs-SRs, a2=Ri-SRi. Finally, one can obtain the error vector of E2 as vector DA in XiYiZi:

(A26) D A = ( ( 2 - a 1 - a 2 ) x i + d ( a 2 - 1 ) 2 , 2 - a 1 - a 2 2 y i , 2 - a 1 - a 2 2 z i ) T .

A5 Calculate AoA and range information in XiYiZi

For a space point (xi,yi,zi) in XiYiZi, which represents a MTSP, Rs can be solved easily as follows:

(A27) R s = ( x i , y i , z i ) R s = x i 2 + y i 2 + z i 2 .

The distance between transmitter Ti and receiver RX is di, as shown in Fig. 4a. Thus, the coordinate value of Ti in XiYiZi is (di,0,0), and Ri can be solved as follows:

(A28) R i = ( x i - d i ) 2 + y i 2 + z i 2 .

Before we calculate the AoAs in XiYiZi, the representation of unit vectors of the x, y and z axes in XiYiZi needs to be known. In XYZ those unit vectors are easily represented as (1,0,0)T, (0,1,0)T, (0,0,1)T. Through the coordinate rotation matrix AR(ψxX,i,ψyY,i,ψzZ,i), one can get those unit vector's representation in XiYiZi as follows:

(A29) n x = ( A 11 , A 21 , A 31 ) T n y = ( A 12 , A 22 , A 32 ) T n z = ( A 13 , A 23 , A 33 ) T .

nx, ny and nz are unit vectors of the x, y and z axes, respectively, and Aij is the element's 3×3 matrix AR(ψxX,i,ψyY,i,ψzZ,i) for i,j=1,2,3. Now the AoA can be obtained as follows:


for 0<θ<180 and 0ϕ<360. When θ=0, we handle it the way as same as in Sect. A2.

Code availability

The program to calculate the 3D spatial resolution distributions is available at (Zhong, 2021).

Data availability

No data sets were used in this article.


The supplement related to this article is available online at:

Author contributions

WZ, XX and WY designed the study. WZ deduced the formulas and wrote the program. WZ wrote the paper for the first version. XX supervised the work and provided valuable comments. IMR revised the paper. All of the authors discussed the results and commented on the paper.

Competing interests

The authors declare that they have no conflict of interest.


Wei Zhong thanks Jia Mingjiao for useful discussions and Zeng Jie for checking the equations in the manuscript.

Financial support

This research has been supported by the B-type Strategic Priority Program of CAS Grant (grant no. XDB41000000), the National Natural Science Foundation of China (grant nos. 41774158, 41974174, 41831071, and 41904135), the CNSA pre-research Project on Civil Aerospace Technologies (grant no. D020105), the Fundamental Research Funds for the Central Universities (grant no. YD3420002004), the Joint Open Fund of Mengcheng National Geophysical Observatory (grant no. MENGO‐202009), and the Open Research Project of Large Research Infrastructures of CAS “Study on the interaction between low/mid-latitude atmosphere and ionosphere based on the Chinese Meridian Project”.

Review statement

This paper was edited by Markus Rapp and reviewed by two anonymous referees.


Browning, K. A. and Wexler, R.: The Determination of Kinematic Properties of a Wind Field Using Doppler Radar, J. Appl. Meteorol., 7, 105–113,<0105:Tdokpo>2.0.Co;2, 1968. 

Ceplecha, Z., Borovička, J., Elford, W. G., ReVelle, D. O., Hawkes, R. L., Porubčan, V., and Šimek, M.: Meteor Phenomena and Bodies, Space Sci. Rev., 84, 327–471,, 1998. 

Chau, J. L., Stober, G., Hall, C. M., Tsutsumi, M., Laskar, F. I., and Hoffmann, P.: Polar mesospheric horizontal divergence and relative vorticity measurements using multiple specular meteor radars, Radio Sci., 52, 811–828,, 2017. 

Chau, J. L., Urco, J. M., Vierinen, J. P., Volz, R. A., Clahsen, M., Pfeffer, N., and Trautner, J.: Novel specular meteor radar systems using coherent MIMO techniques to study the mesosphere and lower thermosphere, Atmos. Meas. Tech., 12, 2113–2127,, 2019. 

Dorey, J., Blanchard, Y., and Christophe, F.: Le projet RIAS, une approche nouvelle du radar de surveillance aérienne, Onde Electr., 64, 15–20, 1984. 

Hocking, W. K.: Spatial distribution of errors associated with multistatic meteor radar, Earth Planets Space, 70, 93,, 2018. 

Hocking, W. K., Fuller, B., and Vandepeer, B.: Real-time determination of meteor-related parameters utilizing modem digital technology, J. Atmos. Sol.-Terr. Phy., 63, 155–169,, 2001. 

Holdsworth, D. A., Reid, I. M., and Cervera, M. A.: Buckland Park all-sky interferometric meteor radar, Radio Sci., 39, RS5009,, 2004. 

Jacobi, C., Hoffmann, P., and Kürschner, D.: Trends in MLT region winds and planetary waves, Collm (52 N, 15 E), Ann. Geophys., 26, 1221–1232,, 2008. 

Jia, M. J., Xue, X. H., Dou, X. K., Tang, Y. H., Yu, C., Wu, J. F., Xu, J. Y., Yang, G. T., Ning, B. Q., and Hoffmann, L.: A case study of A mesoscale gravity wave in the MLT region using simultaneous multi-instruments in Beijing, J. Atmos. Sol.-Terr. Phy., 140, 1–9,, 2016. 

Jones, J., Webster, A. R., and Hocking, W. K.: An improved interferometer design for use with meteor radars, Radio Sci., 33, 55–65,, 1998. 

Kang, C.: Meteor radar signal processing and error analysis, Dissertations & Theses – Gradworks, Department of Aerospace Engineering and Sciences, University of Colorado, US, 2008. 

Philippe, W. and Corbin, H.: On the Analysis of Single-Doppler Radar Data, J. Appl. Meteorol., 18, 532–542,<0532:OTAOSD>2.0.CO;2, 1979. 

Reid, I. M.: Some aspects of Doppler radar measurements of the mean and fluctuating components of the wind-field in the upper middle atmosphere, J. Atmos. Terr. Phys., 49, 467–484,, 1987. 

Reid, I. M.: MF and HF Spaced Antenna radar techniques for investigating the dynamics and structure of the 50 to 110 km height region: A review, Progress in Earth and Planetary Science, 2, 33,, 2015. 

Spargo, A. J., Reid, I. M., and MacKinnon, A. D.: Multistatic meteor radar observations of gravity-wave–tidal interaction over southern Australia, Atmos. Meas. Tech., 12, 4791–4812,, 2019. 

Stober, G., Sommer, S., Rapp, M., and Latteck, R.: Investigation of gravity waves using horizontally resolved radial velocity measurements, Atmos. Meas. Tech., 6, 2893–2905,, 2013. 

Stober, G. and Chau, J. L.: A multistatic and multifrequency novel approach for specular meteor radars to improve wind measurements in the MLT region, Radio Sci., 50, 431–442,, 2015. 

Stober, G., Chau, J. L., Vierinen, J., Jacobi, C., and Wilhelm, S.: Retrieving horizontally resolved wind fields using multi-static meteor radar observations, Atmos. Meas. Tech., 11, 4891–4907,, 2018. 

Urco, J. M., Chau, J. L., Weber, T., Vierinen, J. P., and Volz, R.: Sparse Signal Recovery in MIMO Specular Meteor Radars With Waveform Diversity, IEEE T. Geosci. Remote, 57, 10088–10098,, 2019. 

Vaudrin, C. V., Palo, S. E., and Chau, J. L.: Complex Plane Specular Meteor Radar Interferometry, Radio Sci., 53, 112–128,, 2018. 

Vierinen, J., Chau, J. L., Pfeffer, N., Clahsen, M., and Stober, G.: Coded continuous wave meteor radar, Atmos. Meas. Tech., 9, 829–839,, 2016. 

Vierinen, J., Chau, J. L., Charuvil, H., Urco, J. M., Clahsen, M., Avsarkisov, V., Marino, R., and Volz, R.: Observing Mesospheric Turbulence With Specular Meteor Radars: A Novel Method for Estimating Second-Order Statistics of Wind Velocity, Earth and Space Science, 6, 1171–1195,, 2019.  

Wilhelm, S., Stober, G., and Chau, J. L.: A comparison of 11-year mesospheric and lower thermospheric winds determined by meteor and MF radar at 69  N, Ann. Geophys., 35, 893–906,, 2017. 

Xue, X. H., Dou, X. K., Lei, J., Chen, J. S., Ding, Z. H., Li, T., Gao, Q., Tang, W. W., Cheng, X. W., and Wei, K.: Lower thermospheric-enhanced sodium layers observed at low latitude and possible formation: Case studies, J. Geophys. Res.-Space, 118, 2409–2418,, 2013. 

Yi, W., Xue, X. H., Reid, I. M., Younger, J. P., Chen, J. S., Chen, T. D., and Li, N.: Estimation of Mesospheric Densities at Low Latitudes Using the Kunming Meteor Radar Together With SABER Temperatures, J. Geophys. Res.-Space, 123, 3183–3195,, 2018. 

Younger, J. P. and Reid, I. M.: Interferometer angle-of-arrival determination using precalculated phases, Radio Sci., 52, 1058–1066,, 2017. 

Zhong, W.: 3D Multistatic Meteor Radar Error Anlanyze, OSF [code],, 2021.