Document [original]

Astronomy

Astrophysics

A&A 663, A154 (2022)

https://doi.org/10.1051/0004-6361/202142288

A new least squares adjustment approach

for the determination of linear meteor trajectories

A. Margonis1, G. Malissiovas1, F. Neitzel1, and J. Oberst1,2

1Institute of Geodesy and Geoinformation Science, Technische Universität Berlin, 10623 Berlin, Germany

e-mail: [email protected]

2Institute of Planetary Research, German Aerospace Center (DLR), 12489 Berlin, Germany

Received 23 September 2021 / Accepted 7 March 2022

ABSTRACT

Context. Precise astrometric measurements performed on meteor images are required to derive the dynamical parameters of a mete-

oroid. As a consequence, the measurements carried out in this initial step will have a strong impact on the dynamical solution of an

orbiting meteoroid. These measurements relate to the position of the meteor defined by the positions of pixels along its path, as well

as by their uncertainties. Therefore, the use of all available information is of great importance for the subsequent processing steps.

Aims. This paper examines a new geometrical approach for computing the trajectory of a meteor from multi-station observations.

The model considers a more general weighting scheme based on existing stochastic information from the measurements, including the

geometry between each station and the observed meteor.

Methods. We present a novel mathematical model for least squares adjustment of the linear meteor trajectories within the Gauss-

Helmert model, which allows the use of stochastic information from the measured direction vectors from multiple stations. Additionally,

an extended stochastic model is presented that takes into account the geometric relationship between each station and the observed

meteor as a weight component for each group of observations.

Results. The solution of the new approach is demonstrated on a synthetic meteor example, with observations generated from multiple

stations with differing precision. The geometric configuration of the stations has been chosen in such a way that it creates the necessity

to include stochastic information for the observed direction vectors for a realistic solution. The results of the newly developed approach

are compared with those from established methods in the literature. Future investigations and optimisations for developing an even

more improved meteor trajectory model are being addressed.

Key words. meteorites, meteors, meteoroids – methods: analytical – planets and satellites: atmospheres

1. Introduction

Small fragments of comets and asteroids of sizes ranging from

micrometres (µm) to millimetres (mm) moving in various orbits

around the Sun remain unexplored for most of the time. Over

a long period of time, the initial orbits of these dust particles,

known as meteoroids, change due to interactions with other plan-

ets or non-gravitational forces, or both. As a result, a small

fraction of the meteoroid population will inevitably encounter

the orbits of planets with an atmosphere and produce mete-

ors or meteor showers. Detailed measurements of their dynamic

behaviour and properties during their atmospheric entry can be

used for the geometric reconstruction of their trajectories and the

subsequent determination of their heliocentric orbits.

A number of different methods for the trajectory deter-

mination using optical systems have been developed over

recent years. The two most well-established approaches are

the plane-intersection method described by Ceplecha (1987)

and the straight line least squares (SLLS) method published

by Borovicka (1990). The plane intersection method fits a 2D

plane from line-of-sight (LoS) measurements for each station.

The trajectory of a meteor is then calculated by intersecting

the planes generated from at least two stations in space. In the

present paper, we refer to this travelled path by a meteor as

the trajectory line. In the case of multi-station observations, the

best trajectory is calculated as the average line considering the

geometry between each individual trajectory solution. On the

other hand, the SLLS method provides a least squares solution

for the trajectory line by minimising the sum of the squared

orthogonal distances of all measured direction vectors to the

meteor line in space. Even though the latter approach allows each

LoS measurement to influence the solution, measurements from

a single station are given the same weight depending solely on

the station–meteor geometry. Uncertainties of single measure-

ments have an impact on the determined trajectory and must

therefore be taken into account. This problem becomes evident

when measurements provided by several camera systems with

different characteristics (e.g. resolution, optics, etc.) are used for

the estimation of the trajectory.

Recently, several authors noted the need to consider the

standard deviations of measurements in the computation of the

meteor trajectory (Gural 2012;Sansom et al. 2019;Vida et al.

2020). In the study by Vida et al. (2020), the estimated errors of

the measurements contribute to the trajectory solution as follows:

first an estimate of the trajectory line is computed by applying the

plane-intersection method, which is later used as starting values

for the iterative procedure of the SLLS method. By introducing

Gaussian noise to the original measurements with standard devi-

ations based on the errors of the first trajectory estimation, a new

solution is calculated. This iterative process continues until the

solution converges, indicating the best fit, considering dynamic

and geometric parameters of the meteoroid.

Article published by EDP Sciences A154, page 1 of 12

A&A 663, A154 (2022)

Jeanne et al. (2019) presented a data reduction pipeline which

they applied to image data acquired by the Fireball Recovery

and InterPlanetary Observation Network (FRIPON) (Colas et al.

2012). Similarly to previous studies, estimation of the meteor

trajectory is based on the SLLS method. An empirical weight

was introduced for each camera system modifying the objective

function to be minimised under the assumption that systematic

errors, characterised by imperfections present at each station,

dominate the measurement errors on the images. As a result,

the different characteristics of each camera system can be taken

into account when computing the trajectory of a meteor observed

from multiple stations.

Gural (2012) presented a solution for the determination of the

trajectory line of a meteoroid, which incorporates the dynamic

behaviour of the motion of the body. The initial velocity, that

is the velocity before any deceleration occurs, the decelera-

tion parameters, and the position of the meteoroid during its

atmospheric entry are determined by means of a least squares

adjustment. This solution strategy assumes that identical char-

acteristics of the meteoroid will be observed from different

camera systems, allowing the computation of possible time-

offsets between stations due to synchronisation deficiencies.

Empirical errors are estimated for the unknown parameters of

the model in terms of Monte Carlo simulations. This approach

has been used by several authors (Jenniskens et al. 2011;Vida

et al. 2020) and was extended by Egal et al. (2017) by implement-

ing various optimisation techniques for solving the minimisation

problem.

The camera systems of the Desert Fireball Network (DFN)

are able to register a meteor event with submillisecond precision

allowing point-to-point triangulations of multi-station observa-

tions (Howie et al. 2017;Sansom et al. 2019). High temporal

resolution is achieved using a liquid crystal shutter which gen-

erates a unique pattern of short and long pulse widths along

the meteor trail. The pattern is then decoded to binary elements

(0,1) corresponding to a specific time interval within a single

exposure. Each camera station is equipped with a GNSS receiver

which provides absolute timing for every recorded meteor event.

Sansom et al. (2017) presented a sequential Monte Carlo solu-

tion incorporating dynamical as well as physical parameters to

the trajectory solution of the meteor in three dimensions.

In the studies of Ceplecha (1987) and Borovicka (1990), the

standard deviations of the LoS measurements were, in some

sense, disregarded and the stochastic model is based solely

on the geometry between each station and the meteor. In the

most recent studies, a combined solution for the 3D motion of

the meteoroid is provided, including velocity and deceleration

parameters in a single step. In these cases, the standard devia-

tions of the measurements are empirically assigned using Monte

Carlo simulations.

In the present paper, we propose a new geometrical approach

for fitting a straight line to multi-station LoS measurements

inspired by the solution of Borovicka (1990). Our goal is to take

into account the stochastic information of the LoS measurements

and the geometry between each station and the recorded meteor,

resulting to a more realistic stochastic model. The newly pre-

sented solution strategy for the discussed problem is based on the

development of a new functional model, where the LoS measure-

ments are treated as observations. The inclusion of more general

stochastic models allows weights to be computed individually or

as a combination of (a) the precision of the observed LoS mea-

surements and (b) the geometrical configuration between each

station and the recorded meteor.

(

)

VM(aM, bM, cM

)

(

)

(

)

Fig. 1. Schematic illustration of the observation of a linear trajectory

from two locations.

Only overdetermined configurations are considered, where

the solution for the searched parameters is determined by means

of a least squares adjustment. For the considered case of trajec-

tory determination, this leads to a non-linear adjustment problem

with condition equations and general weights for the observa-

tions. We therefore present a least squares solution based on the

Gauss-Newton approach. Such a solution can be determined by

an iteratively linearised Gauss-Helmert (GH) model, as shown in

(Ghilani 2018, p. 497 ff.), (Niemeier 2008, p. 172 ff.), or (Perovi´

2005, p. 203). Furthermore, we avoid the pitfalls mentioned by

Pope (1972) and Lenzmann & Lenzmann (2004) when dealing

with the linearisation of non-linear adjustment problems.

The paper is organised in four sections: Sect. 2describes

the main objective of this work and presents current approaches

for estimating meteor trajectories. Section 3presents a detailed

description of the new approach including the mathematical for-

mulae for the adjustment solution. In Sect. 4we describe a

numerical example and present the results of the trajectory solu-

tion by means of two different methods. Finally, in Sect. 5we

summarise the results presented in the previous section and draw

conclusions.

2. Existing approaches for estimation of a 3D linear

trajectory

In this section, we give a short overview of the two predominant

approaches in the literature for estimating a 3D linear trajectory

of a body, in particular the parameters of a spatial straight line.

2.1. Problem definition

Let us assume a meteoroid observed from a station PS

(XS,YS,ZS) in a 3D Cartesian coordinate system. The rela-

tionship between an image plane and the 3D space can be

mathematically established by performing a geometric calibra-

tion of the camera system that acquired the image. Hence, from

the measured pixel positions of a meteor trail in the image plane,

the direction vectors Vi(ai,bi,ci)and their standard deviations

can be derived, with i=1, ..., nand nbeing the number of mea-

surements from a station. The components ai,bi, and cican be

seen as the coordinate difference of a point on each direction vec-

tor and the point of the station PS. An illustration of the problem

using two stations is presented in Fig. 1.

A154, page 2 of 12

A. Margonis et al.: New meteor trajectory solution

x-direction

y-direction

z-direction

Station 1

Station 2

Linear Trajectory

Fig. 2. Geometric representation of the solution estimated by the plane

intersection method from two stations.

In this study, we assume that the body is moving on a lin-

ear trajectory, that is, a spatial straight line. This assumption was

found to be valid for meteors observed by commercial all-sky

camera systems, because any deceleration cannot be measured

(Jeanne et al. 2019). The goal is to estimate the direction of the

line in which the body moves, a point on that line, and their

uncertainties. The equation of a line passing through a point

on the body’s trajectory PM(XM,YM,ZM) parallel to a direction

vector VM(aM,bM,cM)can be described by

X−XM

=Y−YM

=Z−ZM

,(1)

with X,Y,and Zdenoting any point on this line, as shown in

(Bronshtein et al. 2007, p. 217). The components aM,bM, and

cMcan be seen as the coordinate difference of a point on the

direction vector of the body’s trajectory and the point on the

observed body PM. Therefore, the trajectory line of a meteor

can be obtained by computing, for example, the point coordi-

nates XM,YM, and ZMand the components aM,bM, and cMof the

direction vector.

In the following, we give a short overview of two existing

adjustment approaches for the discussed problem, which take

into account different overdetermined configurations for a least

squares solution.

2.2. Linear trajectory from plane intersection

Ceplecha (1987) was the first to address the discussed adjustment

problem, by fitting planes to the direction vectors of each station

using the method of least squares. The requested parameters of

the trajectory line and point can then be computed by finding the

intersection line of the two planes, as is depicted in Fig. 2.

A plane in 3D can be expressed by the general equation

nxX+nyY+nzZ+d=0,(2)

with X,Y,and Zbeing the 3D coordinates of a point that lies in

the plane; see Bronshtein et al. (2007, p. 214). The plane parame-

ters are denoted with nx,ny,nz, and d.Ceplecha (1987) presented

a direct solution for fitting a plane to the observed direction

vectors, where all vector components are taken as observations

subject to random errors. This leads to the condition equations

nx(ai−eai)+ny(bi−ebi)+nz(ci−eci)+d=0,(3)

with the constraint

x+n2

y+n2

z=1,(4)

allowing all planes in space to be computed. For equally

weighted and uncorrelated observations, a least squares solution

can be derived by minimising the objective function

Ω =

i=1

ai+e2

bi+e2

ci.(5)

A direct least squares solution for this non-linear adjustment

is possible, as presented by Ceplecha (1987). A discussion of

this adjustment problem in the context of total least squares

adjustment, taking into account more general weighting cases,

is provided by Malissiovas et al. (2016) and Malissiovas (2019,

pp. 78–81, 113–118).

2.2.1. Linear solution for fitting a plane in 3D

Here it is important to point out that Vida et al. (2020) pre-

sented an approximate linear solution of the problem, assuming –

contrary to assumptions made by Ceplecha (1987) – that coordi-

nates in xand ydirections are free of errors, leading to the linear

observation equations

ci−eci=nxai+nybi+d,(6)

after taking into account the constraint1

nz=1.(7)

The solution of the approximate substitute problem can be

obtained by minimising the sum of squared errors only in the

z-direction

Ω =

i=1

ci,(8)

which leads to linear normal equations that can be solved

directly, that is without linearisation and iterative computation,

using the rules of linear algebra. This solution clearly refers to

a different adjustment problem, as only the errors of the coordi-

nates in the z-direction are taken into account and will lead to a

different solution for the plane parameters from that presented in

the previous subsection.

2.2.2. Line of intersection of two planes in 3D

If the body has been observed from two stations, then two planes

have to be computed in total. The trajectory line of the body is

simultaneously the line of intersection of the two planes and can

be computed, for example, by taking the cross product of the two

planes

VM(aM,bM,cM)=V1(nx1,ny1,nz1)×V2(nx2,ny2,nz2).(9)

Additionally, a point on the linear trajectory PM(XM,YM,ZM)

can be obtained by computing the intersection point of one direc-

tion vector from one station to the plane of the second. For

the interested reader, the equations for intersection of lines and

planes can be found in Bronshtein et al. (2007, p. 218).

In the case where the body has been observed from multiple

stations, then the weighted mean of all the lines of intersection

can give the solution of the trajectory. This was described in

detail by Ceplecha (1987) and by Vida et al. (2020).

1The restriction of Eq. (7) prevents the computation of planes that are

parallel to the z-axis.

A154, page 3 of 12

A&A 663, A154 (2022)

Fig. 3. Simplified geometric representation of the SLLS method.

2.3. Least squares estimation of a linear trajectory

Borovicka (1990) presented an elegant geometrical way to

express this adjustment problem, by taking into account each

individual observed direction from multiple stations. In that

study, a least squares solution was derived for the requested

trajectory line by minimising the sum of squared orthogonal dis-

tances between the line described by each measured direction

vector and the trajectory line. An illustration of this geometric

arrangement is shown in Fig. 3.

2.3.1. Functional model and resulting observation equations

The equation of a line passing through a station PS(XS,YS,ZS),

with a direction vector Vi(ai,bi,ci), towards the observed body

can be described by

X−XS

=Y−YS

=Z−ZS

,(10)

with i=1, ...nand nbeing the number of direction vectors. The

orthogonal distances between the lines of the measured direction

vectors and the trajectory line Eq. (1) can be expressed by

Di=[(XM−XS)(bMci−cMbi)]

p[(bMci−cMbi)2+(aMci−cMai)2+(aMbi−bMai)2]

−[(YM−YS)(aMci−cMai)]

p[(bMci−cMbi)2+(aMci−cMai)2+(aMbi−bMai)2]

+[(ZM−ZS)(aMbi−bMai)]

p[(bMci−cMbi)2+(aMci−cMai)2+(aMbi−bMai)2]

(11)

Utilizing a concept of ‘zero’ observations subject to ran-

dom errors expressed as the orthogonal distances Di,Borovicka

(1990) presumed the overdetermined system Eq. (11) to be non-

linear observation equations of an adjustment problem. In this

case, the unknown parameters are: the coordinates of a point on

the meteor XM,YM, and ZM; and the components aM,bM, and cM

of the direction vector parallel to the meteor, while the coordi-

nates of each station XS,YS, and ZSand the components ai,bi,

and ciof the direction vectors from each station to the recorded

meteor are taken as fixed values.

Therefore, the orthogonal distances Eq. (11) can be written

equivalently with the general formulation

L−e=Φ(X),(12)

with the observation vector Llisting only zero values and the

accompanied error vector eholding the orthogonal distances Di.

The formal vector Φholds the non-linear functional relationship

between the unknown parameters, which are listed in vector

X=[XM,YM,ZM,aM,bM,cM]T

.(13)

As already explained by Borovicka (1990), it is possible to fix

two unknown parameters, for example by setting XM=XM0and

aM=1, with the superscript ‘0’ denoting the starting value.

It is worth noting that setting parameter aM=1prevents the

calculation of all lines in space. For a more general solution,

the constraint a2

M+b2

M+c2

M=1should be taken into account.

Subsequently, the vector of unknown parameters becomes

X=[YM,ZM,bM,cM]T

.(14)

2.3.2. Stochastic model

For normally distributed errors, e∼(0, σ2

0QLL), where σ2

0is

the theoretical variance of the unit weight, we can construct the

diagonal cofactor matrix

QLL =1

σ2







ΣLL10. . . 0

0ΣLL2. . . 0

.....

0 0 . . . ΣLLm







,(15)

with the diagonal variance–covariance matrices (VCMs) ΣLL1,

ΣLL2, and ΣLLmlisting the variances of the observations for each

station, and mbeing the number of all stations. For example, if

the standard deviations of the observations were known, which

is not the case in this adjustment approach, the VCM for the first

station would be

ΣLL1=







σ2

10. . . 0

0σ2

1. . . 0

.....

0 0 . . . σ2







,(16)

with σ1being the standard deviation for the observations of the

first station. Assuming a non-singular cofactor matrix, the weight

matrix can be obtained from

P=Q−1

LL.(17)

As the direction vectors Vi(ai,bi,ci)are introduced as fixed

(error-free) values in this adjustment approach, their precision

measures cannot be taken into account in the stochastic model.

Instead, Borovicka (1990) considers a different weighting con-

cept that depends on the geometry between each station and

the trajectory line. Following Vida et al. (2020), the weights are

calculated as

Pai=arccos (VM·Vi),

wi=sin2Pai,(18)

where VMis the computed direction vector of the trajectory line

given in Eq. (9) and Pais the perspective angle computed as the

A154, page 4 of 12

A. Margonis et al.: New meteor trajectory solution

dot product between each direction vector Viand VM. Weights wi

are calculated for all available direction vectors of a station. Of

all direction vectors Vi, the one closest to the meteor trajectory is

chosen, yielding the highest weight. This occurs when a meteor

is moving in a line perpendicular to the direction vector, for

instance when the perspective angle is Pa=90◦. Subsequently,

the computed weights from Eq. (18) can be directly introduced

into the weight matrix for the first station as follows

P1=







w10. . . 0

0w1. . . 0

.....

0 0 . . . w1







.(19)

The weight matrices P1,P2,Pmcontaining the weights for each

station can be included in the weight matrix P, which yields







P10. . . 0

0 P2. . . 0

.....

0 0 . . . Pm







.(20)

2.3.3. Least squares solution within the Gauss-Markov model

A least squares solution within the Gauss-Markov (GM) model

can be derived following the Gauss-Newton approach, that is

by linearising the non-linear functional model in Eq. (11) and

obtaining an iterative solution for the unknown parameters by

minimizing the objective function

Ω(e)=eTPe.(21)

Detailed explanations of this iterative procedure are given

by, for example, Niemeier (2008, p. 137 ff.) and Ghilani (2018,

p.207 ff.).

2.4. Discussion and open questions

In this section we review two of the most popular approaches

for the estimation of linear trajectories. These are the line of

intersection of planes from Ceplecha (1987) and the SLLS from

Borovicka (1990). This leads us to two open questions:

– Is there a different geometrical approach that can lead to

an adjustment problem that directly involves the measured

direction vectors from multiple stations?

– Is it possible to incorporate more realistic stochastic mod-

els in the adjustment, where the standard deviations of the

measured direction vectors would be used directly?

3. A new approach for the determination of linear

trajectories

Inspired by the solution of Borovicka (1990), in this section we

present a new geometric point of view for the problem being

discussed. This leads to an adjustment model that allows the

inclusion of the stochastic parameters of the observed direction

vectors directly.

The idea is based on the fact that the direction vectors from

each station should be coplanar and not parallel to the direction

of the meteor. Individual planes would describe the coplanarity

between each individual station with the observed meteor, as

illustrated in Fig. 4.

x-direction

y-direction

z-direction

Station 1

Linear Trajectory

Observed direction

Adjusted direction

Fig. 4. New geometrical perspective based on the coplanarity constraint

between each measured direction vector and the trajectory line.

Two vectors or lines in space would be on the same plane

if they cross each other, that is, they would not be skewed and

not parallel. According to Bronshtein et al. (2007, p. 218), the

determinant of the two vectors should be equal to zero, which

yields



(XM−XS) (YM−YS) (ZM−ZS)

aMbMcM

aibici



=0,(22)

which can be described equivalently by

(XM−XS)

bMcM

bici

−(YM−YS)

aMcM

aici

+(ZM−ZS)

aMbM

aibi

=0,

(23)

(XM−XS)(bMci−cMbi)−(YM−YS)(aMci−cMai)

+(ZM−ZS)(aMbi−bMai)=0.

(24)

In this case, the direction vectors [ai,bi,ci]Tare observations,

while the unknown parameters are the coordinates of a point on

the meteor XM,YM, and ZM, and the components aM,bM, and cM

of the direction vector parallel to the meteor. The coordinates of

each station XS,YS, and ZSare considered as fixed values.

3.1. Functional model and resulting condition equations

As the direction vectors Vi(ai,bi,ci)are introduced as observa-

tions subject to random errors, the observation errors eai,ebi, and

ecihave to be introduced in the functional model, yielding from

Eq. (24) the non-linear condition equations

(XM−XS)bMci−eci−cMbi−ebi

−(YM−YS)aMci−eci−cMai−eai

+(ZM−ZS)aMbi−ebi−bMai−eai =0,

(25)

or written equivalently, following the general form of condition

equations, as

φi(XM,YM,ZM,aM,bM,cM,eai,ebi,eci)=0,(26)

A154, page 5 of 12

A&A 663, A154 (2022)

with i=1, ...nand nbeing the number of direction vectors and φi

representing non-linear differentiable functions of the unknown

parameters and the errors. This can be expressed in matrix

notation as

Φ(X,e)=0,(27)

with the formal vector Φcontaining the implicit non-linear

functional relationship between the error vector

e=ea1,eb1,ec1,· · · ,ean,ebn,ecnT,(28)

which contains the errors of the observations from all stations

and the vector of unknown parameters

X=[XM,YM,ZM,aM,bM,cM]T.(29)

In this study, we apply the same restrictions for the unknowns

as in Borovicka (1990) to obtain a simpler equation system for

the following computations. Therefore, setting XM=XM0and

aM=1, the vector of unknown parameters results in

X=[YM,ZM,bM,cM]T.(30)

3.2. Stochastic model taking into account variances and

covariances of the observations

In this adjustment problem, the cofactor matrix will have the

same structure as in Eq. (15), while for example, the VCM for

the first station contains the variances and covariances of the

observed directions:

ΣLL1=







σ2

a1σa1b1σa1c10. . . 0

σb1a1σ2

b1σb1c10. . . 0

σc1a1σc1b1σ2

c10. . . 0

0 0 0 ...

0 0 0







.(31)

For a non-singular dispersion matrix, the respective weight

matrix P can be obtained from Eq. (17).

3.3. Stochastic model taking into account the geometry

between station and meteor

The presented stochastic model for this adjustment problem can

be extended in such a way that the geometry information between

each station and the observed meteor is also taken into account.

This is allowed, as the weights in the adjustment are relative in

the sense that the relationship between the group of observations

between each station is of importance. Therefore, an individual

variance for the unit weight can be introduced in

QLL =







σ2

ΣLL10. . . 0

σ2

ΣLL2. . . 0

.....

0 0 . . . 1

σ2

ΣLLm







.(32)

The individual variances of the unit weight can be seen as

scaling factors for the observations of each station. Therefore, the

computed weights wifrom Eq. (18), expressed by the geometry

of each station and the recorded meteor, can be introduced in this

stochastic model as a scale factor for each group of observations,

yielding

QLL =







ΣLL10. . . 0

ΣLL2. . . 0

.....

0 0 . . . 1

ΣLLm







.(33)

3.4. Least squares solution within the Gauss-Helmert model

We present a solution for the resulting adjustment problem that

is based on the Gauss-Newton approach. This involves a lin-

ear approximation of the non-linear condition equations and an

approximation of both the unknown parameters and the unknown

errors. The studies of Pope (1972), Lenzmann & Lenzmann

(2004), and Neitzel & Petrovic (2008) have shown theoretically

and numerically that only such a linearisation of the functional

model can lead to a rigorous solution within the GH model.

The first-order Taylor series approximation of the formal

vector Φ(X,e)at the point X0and e0yields

Φ(X,e)≈Φ(X0,e0)+∂Φ(X,e)

∂XX0,e0

(X−X0)

+∂Φ(X,e)

∂eX0,e0

(e−e0)=0.

(34)

Therefore, we compute as a first step the partial derivatives of

the condition equations with respect to the unknown parameters

of the problem and form the Jacobian matrix

A=∂Φ(X,e)

∂XX0,e0

=hjYM|jZM|jbM|jcMi,(35)

with the respective vectors

jYM=





cM0(ai−eai

0)−aM0(ci−eci







,(36)

jZM=





aM0(bi−ebi

0)−bM0(ai−eai







,(37)

jbM=





(XM0−XS)(ci−eci

0)−(ZM0−ZS)(ai−eai







(38)

and

jcM=





(YM0−YS)(ai−eai

0)−(XM0−XS)(bi−ebi







.(39)

A154, page 6 of 12

A. Margonis et al.: New meteor trajectory solution

Subsequently, a second Jacobian matrix is formulated that con-

tains the partial derivatives of the condition equations with

respect to the unknown errors

B=∂Φ(X,e)

∂eX0,e0







B10. . . 0

0 B2. . . 0

.....

0 0 . . . Bm







,(40)

with matrices B1,B2,...,Bmconstructed individually for each

station. For example, the respective matrix for the first station

can be written as

B1=







jea1jeb1jec10 0 0 0 . . .

0 0 0 jea1jeb1jec10. . .

0 0 0 0 0 0 ...







,(41)

with the respective partial derivatives being

jea=bM0(ZM0−ZS)−cM0(YM0−YS),(42)

jeb=cM0(XM0−XS)−aM0(ZM0−ZS),(43)

and

jec=aM0(YM0−YS)−bM0(XM0−XS).(44)

Introducing the Jacobian matrices into Eq. (34) results in the

linearised condition equations as

Φ(X0,e0)+A(X−X0)+B(e−e0)=0.(45)

Following the usual notation, we introduce the vector of

misclosures

w=Φ(X0,e0)−Be0,(46)

and the vector of reduced unknowns x=X−X0into Eq. (45),

resulting in the linearized functional model2

Ax +Be +w=0.(47)

Due to the implicit functional relationship between the unknown

parameters of this type of adjustment problem, a least squares

solution can be derived by searching for stationary points of the

Lagrange function:

K(x,e,k)=eTPe −2kT(Be +Ax +w),(48)

where kdenotes the vector of Lagrange multipliers. The itera-

tive procedure for obtaining a least squares solution within the

GH model can be found in many textbooks and studies in the

geodetic literature, such as those of Niemeier (2008, p. 177) or

Malissiovas (2019, p. 40 ff.). In the following we present this

procedure for solving the discussed problem.

The system of normal equations can be derived by

∂K

∂eT=2P˜

e−2BTˆ

k=0

⇒˜

e=QLLBTˆ

k,(49)

2Similar to the definitions from Wolf (1978)orPerovi´

c(2005, p. 203),

the developed linearised functional model Eq. (47) accompanied by the

stochastic model Eq. (15) results in the well-known GH model.

∂K

∂xT=−2ˆ

kTA=0

⇒ATˆ

k=0,(50)

∂K

∂kT=−2(B˜

e+Aˆ

x+w)=0.(51)

The tilde symbol denotes a prediction, while the circumflex an

estimator. Introducing the error vector from Eqs. (49) into (51)

yields the reduced normal equations

BQLLBTˆ

k+Aˆ

x+w=0.(52)

Furthermore, Eqs. (50) and (52) can be expressed in block matrix

form as







BQLLBTA

AT0



















=





−w







,(53)

that leads to the solution of













=





BQLLBTA

AT0







−1





−w







.(54)

For a non-singular product [BQLLBT]it is possible to explicitly

express an estimate for the vector of reduced unknowns

x=N−1n,(55)

and the vector of Lagrange multipliers

k=−BQLLBT−1(Aˆ

x+w),(56)

after introducing the normal matrix

N=−ATBQLLBT−1A(57)

and the vector

n=ATBQLLBT−1w.(58)

Finally, an estimate for the vector of unknown parameters can be

derived by

Xi=ˆ

xi+X0,(59)

and a prediction for the errors from Eq. (49). A series of lin-

earised adjustments within the GH model have to be computed

because of the iterative nature of the Gauss-Newton approach.

Therefore, the estimated vector ˆ

Xis used as the vector of initial

values for the subsequent iteration step:

i+1=ˆ

Xi,

and the predicted error, excluding for a moment its random

character, as

i+1=˜

ei.(60)

The vector of adjusted observations can be finally computed as

L=L−˜

efinal,(61)

with the subscript ‘final’ denoting the last iteration step of the

procedure.

A154, page 7 of 12

A&A 663, A154 (2022)

Fig. 5. Map illustrating the geometry between the stations and the pro-

jection of the synthetic meteor (ground track) depicted with the red

arrow.

Error estimation

The error estimates of the adjustment results can be easily

computed in terms of the cofactor matrices of the unknown

parameters Qˆ

Xˆ

Xof the errors Q˜

e˜

eand the adjusted observations

Qˆ

Lˆ

L. These cofactor matrices can be obtained using standard

procedures, presented for example in Niemeier (2008, p. 178)

or Malissiovas (2019, p. 44 ff.).

4. Numerical example

A set of synthetic observations were generated using the simula-

tion tool implemented for the study of meteor showers (Margonis

et al. 2018). These observations were used as an input to the

different methods described in this work. The intention of our

synthetic setup is to reveal the limitations of a trajectory solution

when geometric aspects alone are taken into account.

Our hypothetical scenario consists of four observing stations

which simultaneously register a meteor. Two of the stations (A,

D) have a favourable geometry but the measurements are of rel-

atively low precision (e.g. use of low-resolution cameras). The

other two stations (B, C) have a less favourable geometry, as

they are located close to the projection of the trajectory line.

However, these measurements are more precise compared to the

ones made from A and D (e.g. use of high-resolution cameras;

Fig. 5). Three observation vectors from each station were used:

at the beginning, at the middle, and at the end of the meteor trail.

The trajectory solution consists of a point on the meteor trail and

the moving direction.

For the simulated meteor event, the following entry parame-

ters were assigned:

(i) geodetic coordinates of the meteor position at t0=0;

(ii) duration of the meteor ∆tin seconds;

(iii) speed of the meteor in km s−1;

(iv) horizontal coordinates of the trail orientation with respect

to the local plane of a station;

(v) geodetic coordinates of the observing stations.

The geodetic coordinates at the beginning and end point of the

meteor along with those of each station are given in Table 1. The

moving orientation of the meteor was defined with reference to

Table 1. Geodetic coordinates of the synthetic meteor and the stations.

Positions Lat (◦) Lon (◦) Altitude (m)

Station A 52.86894 12.16539 51

Station B 52.01851 12.97838 158

Station C 53.85882 13.26402 251

Station D 52.90209 13.71446 42

Meteor (t= 0) 52.73878 13.33909 53 727

Meteor (t=t∆t) 52.88012 13.29802 46 283

the local coordinate system of station A. The azimuth of the mov-

ing direction was Az =10◦measuring from north to west while

the elevation angle, which corresponds to the entrance angle of

the meteor, was El =25◦. The speed of the meteor was assumed

to be v=19 km s−1. From the known orientation and moving

speed of the meteor, the velocity vector can be calculated.

For our computations, all direction (unit) vectors were trans-

formed to the Earth-centred Earth-fixed (ECEF) coordinate

system, with

X=cos Az sin El,

Y=sin Az sin El,(62)

Z=cos El.

The direction vectors from each station (ai,bi,ci)to the meteor

at time intervals t0,t∆t/2, and t∆twere calculated by simply sub-

tracting the position vectors of the meteor PM(XM,YM,ZM) from

the position vectors of the station PS(XS,YS,ZS).

In our example, we assume that the observations between the

four stations have different standard deviations. The difference in

the quality of the measurements, when derived from optical sys-

tems, may result for example from the use of sensors of different

quality or from following different calibration procedures. Given

in spherical coordinates at first, the direction vectors from sta-

tions A and D were taken with a standard deviation of σA,D=1◦

in both azimuth and elevation, whereas from stations B and

C, σB,C=0.1◦. Assuming uncorrelated spherical coordinates for

simplicity, we can form a diagonal variance–covariance matrix

ΣAzEl.

As all computations in this article were made in a rectangu-

lar coordinate system, the initial stochastic information has to be

converted to the same coordinate system by applying the law of

error propagation. Therefore, the cofactor matrix of the direc-

tion vectors in rectangular coordinates can be computed using

the rules of error propagation, resulting in

Σabc =FΣAzElFT.(63)

Matrix Fcontains the partial derivatives of the non-linear rela-

tionships in Eq. (62) with respect to the spherical coordinates.

The derived VCM Σabc can be used to compute the cofactor

matrix QLL for the adjustment in the GH model from Eq. (15).

In the next step, the calculated variances were used to

generate normally distributed random errors for each station

individually, which were added respectively to the LoS mea-

surements. For the interested reader, the values of the observed

direction vectors, their standard deviations, and the geometrical

weights for each station are listed in Tables A.3 and A.5 of the

Appendix, respectively. Consequently, two least squares solu-

tions were computed for the simulated data. First, by applying

the SLLS method described in Sect. 2.3 and the newly presented

A154, page 8 of 12

A. Margonis et al.: New meteor trajectory solution

0 2 4 6 8 10 12

# measurement

1.5

1.0

0.5

0.0

0.5

1.0

1.5

D (km)

SLLS Stations A & D

Stations B & C

0 2 4 6 8 10 12

# measurement

1.5

1.0

0.5

0.0

0.5

1.0

1.5

D (km)

GH Stations A & D

Stations B & C

0 2 4 6 8 10 12

# measurement

1.5

1.0

0.5

0.0

0.5

1.0

1.5

D (km)

GH ext. Stations A & D

Stations B & C

Fig. 6. Values of the orthogonal distances calculated for each approach following Eq. (11). From left to right: computed values for the SLLS

method, the GH, and the GH with the extended stochastic model are shown, respectively. Highlighted symbols in orange indicate calculated values

derived from low-resolution observations, i.e. measurements obtained from stations A and D. The straight dashed line represents the estimated

trajectory line for each approach.

Table 2. Estimated values for the unknown parameters and their stan-

dard deviations calculated from the SLLS method, the GH, and the GH

with the extended stochastic model, being compared with the simulated

values.

Methods Ym(km)Zm(km)bmcm

Simulated 900.33311 5095.75595 0.41571 –0.22795

SLLS 899.44130 5096.21653 0.47150 –0.28180

Std Dev 0.2592 0.2780 0.0322 0.0364

GH 899.27247 5096.44356 0.41203 –0.24465

Std Dev 0.5734 0.4081 0.0678 0.0514

GH ext. 899.27320 5096.44280 0.41968 –0.22499

Std Dev 0.4287 0.3036 0.0506 0.0381

solution within the GH model from Sect. 3.4. For completeness,

the solution of the GH with the extended stochastic model is

also presented. The values of the estimated parameters and their

standard deviations are shown in Table 2.

The initial values for the unknown parameters Xgiven

in Eq. (29) are calculated using the plane-plane inter-

section method between the first two stations A and B.

For our numerical example, the approximated values X0=

hX0

M,Y0

M,Z0

M,a0

M,b0

M,c0

Miwere equal to X0

M= 3794.49124,

M= 898.85936, Z0

M= 5095.39098 for the point on the estimated

trajectory line, while the vector components were a0

M= 0.87506,

M= 0.29339, c0

M=−0.38496.

An interesting finding is the amount of discrepancy in the

direction of the trajectory line, which can be interpreted as an

angular displacement Drad of the radiant, that is the point in

the sky where the meteor appears to originate. This will be

equal to Drad =3.55◦,Drad =0.88◦, and Drad =0.25◦for the SLLS

method, the GH model and the GH with the extended stochastic

model respectively.

The solutions presented in Table 2and the angular distances

between the simulated and the estimated trajectory lines suggest,

as anticipated, a difference in the results of the two methods.

The solution obtained by the SLLS method is influenced by

the observations acquired by the low-resolution camera systems

in stations A and D to a greater degree. These two stations

have a better viewing geometry compared to stations B and C

which results in a higher weight in the adjustment. The extended

stochastic model offers a hybrid solution, in that the accuracy of

the measurements and the geometrical configuration between the

stations and the meteor influence the estimated trajectory.

The use of different stochastic models for the observations

leads to different standard deviations for the adjusted unknowns.

Beyond the aim for the lowest possible numerical values for the

standard deviations, there is of course the demand for values that

are as realistic as possible. Under this premise, the solution GH

with an extended stochastic model appears to be the most appro-

priate, as both the measurement accuracy and the measurement

configuration are taken into account in the stochastic model.

The effect of each adjustment approach to the estimated

trajectory line is comprehensively shown in Fig. 6. This is a

graphical representation of the orthogonal distances Dbetween

the measured direction vectors and each estimated trajectory

line, computed by Eq. (11) for all adjustment solutions. As

expected, each set of observed direction vectors for a given sta-

tion contributes differently to the estimation of the trajectory in

the GH approach. The juxtaposition of the low resolution mea-

surements of stations A and D, and the high resolution of stations

B and C is clearly projected to the computed orthogonal dis-

tances, depicted in orange and blue, respectively, in Fig. 6. It

can be seen that low-resolution measurements lead to smaller

weights in the adjustment and higher values for the orthogonal

distances, as these measurements would contribute less to the

solution of the trajectory line.

The solution GH with an extended stochastic model can be

regarded as a reference solution, as both the measurement accu-

racy and the measurement configuration are taken into account

in the stochastic model. The extent to which the other solu-

tions, namely GH and SLLS, differ from this reference solution

depends on the number of observations and stations as well as on

the measurement configuration and can therefore not be given in

a generalised form. However, depending on the type and extent of

information on the accuracies and configuration of the measure-

ments, the user is now able to properly take this information into

account in the stochastic model when computing the trajectory

line of a meteor.

5. Conclusion

This study sets out with the aim of assessing the importance

of including more realistic stochastic information for the direc-

tion vectors when estimating the linear flight path of a mete-

oroid. Detailed explanations of the two most popular approaches

in the literature, namely the plane-intersection method of

A154, page 9 of 12

A&A 663, A154 (2022)

Ceplecha (1987) and the SLLS approach of Borovicka (1990),

are presented in Sect. 2accompanied by the corresponding

mathematical formulations. It is noted that especially the latter

approach postulates a geometry-based stochastic model, disre-

garding any existing a priori stochastic information about the

observed direction vectors.

A novel geometric perspective for the discussed adjustment

problem is presented in Sect. 3, based on the coplanarity condi-

tion between each adjusted direction vector and the linear meteor

trajectory. An iterative least squares solution within the GH

model has been described in detail, which is based on the Gauss-

Newton approach and allows the consideration of the a priori

stochastic information of the observed direction vectors directly

in the stochastic model. Additionally, an extended stochastic

model is presented that allows the integration of the geometry

information between the observed group of directions from each

station.

The results of this study show a discrepancy between the

SLLS method and the GH approach described in this work.

Although a single numerical example is presented here, we antic-

ipate that the computation of the trajectory line for various

geometrical configurations between a meteoroid and the observ-

ing stations will yield different results when implementing the

methods presented in this study. It is therefore possible that the

use of the stochastic information leads to a more realistic solution

for the trajectory line of a meteoroid.

A combined solution which takes into account the full 3D

motion of the meteoroid including velocity and deceleration

parameters and incorporates their uncertainties directly in the

adjustment has yet to be investigated. Other effects that may

divert the flight path of a meteoroid from a straight line, such

as the rotation of the Earth or the influence of the Earth’s gravity

field, need to be considered when calculating the true radiant of a

meteor. Furthermore, physical processes during the atmospheric

entry of a meteoroid, such as fragmentation or aerodynamic

effects, may increase the complexity of a trajectory solution.

Acknowledgements. The authors of this study would like to gratefully acknowl-

edge Dr. Jiˇ

ri Borovi˘

cka for his valuable comments and suggestions for the

improvement of this paper. This work was funded in part by the Deutsch

Forschungsgemeinschaft (DFG – German Research Foundation) under the grant

number OB 124/20-1.

References

Borovicka, J. 1990, Bull. Astron. Inst. Czechoslovakia, 41, 391

Bronshtein, I., Semendyayev, K., Musiol, G., & Muehlig, H. 2007, Handbook of

Mathematics, 5th edn. (Berlin Heidelberg New York: Springer)

Ceplecha, Z. 1987, Bull. Astron. Inst. Czechoslovakia, 38, 222

Colas, F., Zanda, B., Vernazza, P., et al. 2012, in Asteroids, Comets, Meteors

(New York: HarperCollins), 1667, 6426

Egal, A., Gural, P. S., Vaubaillon, J., Colas, F., & Thuillot, W. 2017, Icarus, 294,

Ghilani, C. D. 2018, Adjustment Computations: Spatial Data Analysis, 6th edn.

(Hoboken, NJ, USA: John Wiley and Sons)

Gural, P. S. 2012, Meteor. Planet. Sci., 47, 1405

Howie, R. M., Paxman, J., Bland, P. A., et al. 2017, Meteor. Planet. Sci., 52, 1669

Jeanne, S., Colas, F., Zanda, B., et al. 2019, A&A, 627, A78

Jenniskens, P., Gural, P. S., Dynneson, L., et al. 2011, Icarus, 216, 40

Lenzmann, L., & Lenzmann, E. 2004, Allgemeine Vermessungs-Nachr., 111,

Malissiovas, G. 2019, New Nonlinear Adjustment Approaches for Applications in

Geodesy and Related Fields (München: Ausschuss Geodäsie der Bayerischen

Akademie der Wissenschaften (DGK)), C, 841

Malissiovas, G., Neitzel, F., & Petrovic, S. 2016, J. Geodetic Sci., 6, 43

Margonis, A., Christou, A., & Oberst, J. 2018, A&A, 618, A99

Neitzel, F., & Petrovic, S. 2008, Zeitschrift für Geodäsie, Geoinformation und

Landmanagement, 133, 141

Niemeier, W. 2008, Ausgleichungsrechnung, (in German) 2nd edn. [Adjustment

computations] (New York: Walter de Gruyter)

Perovi´

c, G. 2005, Least Squares (Monograph) (Faculty of Civil Engineering,

University of Belgrade)

Pope, A. 1972, In: Proceedings of the 38th Annual Meeting of the American

Society of Photogrammetry, Washington, DC, 449

Sansom, E. K., Rutten, M. G., & Bland, P. A. 2017, AJ, 153, 87

Sansom, E. K., Jansen-Sturgeon, T., Rutten, M. G., et al. 2019, Icarus, 321,

388

Vida, D., Gural, P. S., Brown, P. G., Campbell-Brown, M., & Wiegert, P. 2020,

MNRAS, 492, 5313

Wolf, H. 1978, Zeitschrift für Vermessungswesen, 103, 41

A154, page 10 of 12

A. Margonis et al.: New meteor trajectory solution

Appendix A: LoS measurements and stochastic models

In this Appendix we provide all necessary data for the computations of the numerical example. Interested readers who need further

information are welcome to contact the authors directly. The rectangular coordinates of the meteor and the stations are given in

Table A.1 and Table A.2 respectively. Table A.3 lists the synthetic direction vectors, which are utilized as observations in the

adjustment for the estimation of the meteor line. Table A.4 lists the submatrices for the VCM of each station computed from (63).

For example, for the first station this will be

ΣLL1=





ΣLLA1 0 0

0ΣLLA2 0

0 0 ΣLLA3







Additionally, taking into account the geometrical weights computed from (18) and listed in Table A.5, it is possible to set up the

cofactor matrix of Eq. (33).

Table A.1. Rectangular coordinates of the meteor.

Meteor XM(km) YM(km) ZM(km)

t= 0 3797.10458 900.33312 5095.7559

t=∆t/2 3789.09276 897.00246 5097.5822

t=∆t3781.08094 893.67181 5099.4085

Table A.2. Rectangular coordinates of the four observing stations.

Station XS(km) YS(km) ZS(km)

A 3771.70290 813.08651 5061.79298

B 3832.95433 883.38393 5004.19521

C 3669.49068 864.99810 5127.69427

D 3745.48131 914.05147 5064.01163

Table A.3. Direction (unit) vectors from all stations to the observed synthetic meteor.

Station Time aibici

t00.2401419174778911 0.9038683766787607 0.3540534099682472

At∆t/20.1971236660301171 0.9018866975181360 0.3843730546366086

t∆t0.1253260775060687 0.8973773044059408 0.4230925995973816

t0-0.3576625477767759 0.1701558770149007 0.9182181001456655

Bt∆t/2-0.4212423719599385 0.1308601758910254 0.8973052121409928

t∆t-0.4763628990887521 0.0935718722687177 0.8742555078990843

t00.9358060180567505 0.2633583489248460 -0.2343277120196155

Ct∆t/20.9390976556344460 0.2519004181290621 -0.2337557967800820

t∆t0.9417787313157842 0.2389037634211161 -0.2365963082223801

t00.8322710397522767 -0.2175601835352761 0.5098945802121861

Dt∆t/20.7463624092241152 -0.2999795592073118 0.5941005118285635

t∆t0.6681614200961595 -0.3764661857632456 0.6417425711856406

A154, page 11 of 12

A&A 663, A154 (2022)

Table A.4. VC submatrices of the observed directions for each station.

0.0002513833961047 -0.0000566430726554 -0.0000258995066094

ΣLLA1 -0.0000566430726554 0.0000532340236820 -0.0000974829602497

-0.0000258995066094 -0.0000974829602497 0.0002664324634582

0.0002498276448996 -0.0000447677009112 -0.0000230805650942

ΣLLA2 -0.0000447677009112 0.0000547897748871 -0.0001055989625644

-0.0000230805650942 -0.0001055989625644 0.0002596124364345

0.0002463473555086 -0.0000267890356692 -0.0000161521973202

ΣLLA3 -0.0000267890356692 0.0000582700642781 -0.0001156552218011

-0.0000161521973202 -0.0001156552218011 0.0002500886633666

0.0000021824925413 -0.0000008109641827 -0.0000010004008464

ΣLLB1 -0.0000008109641827 0.0000008636816566 0.0000004759348846

-0.0000010004008464 0.0000004759348846 0.0000004778701630

0.0000022866402467 -0.0000005301593455 -0.0000011514019885

ΣLLB2 -0.0000005301593455 0.0000007595339512 0.0000003605349230

-0.0000011514019885 0.0000003605349230 0.0000005935268045

0.0000022684337791 -0.0000003045682750 -0.0000012686185047

ΣLLB3 -0.0000003045682750 0.0000007777404188 0.0000002491944879

-0.0000012686185047 0.0000002491944879 0.0000007179141313

0.0000003662641861 -0.0000007071191431 0.0000006679811714

ΣLLC1 -0.0000007071191431 0.0000026799100117 0.0000001879859875

0.0000006679811714 0.0000001879859875 0.0000028789103670

0.0000003485674927 -0.0000006789494599 0.0000006686947000

ΣLLC2 -0.0000006789494599 0.0000026976067051 0.0000001793684326

0.0000006686947000 0.0000001793684326 0.0000028797258403

0.0000003340691621 -0.0000006447317689 0.0000006787527110

ΣLLC3 -0.0000006447317689 0.0000027121050358 0.0000001721811841

0.0000006787527110 0.0000001721811841 0.0000028756560281

0.0000885508375082 0.0000357781304326 -0.0001292706444398

ΣLLD1 0.0000357781304326 0.0002160665822785 0.0000337920506503

-0.0001292706444398 0.0000337920506503 0.0002254191804727

0.0001199753599438 0.0000309985209909 -0.0001350717167122

ΣLLD2 0.0000309985209909 0.0001846420598429 0.0000542883102631

-0.0001350717167122 0.0000542883102631 0.0001971010510086

0.0001383948134294 0.0000229718144836 -0.0001306161807863

ΣLLD3 0.0000229718144836 0.0001662226063573 0.0000735938560663

-0.0001306161807863 0.0000735938560663 0.0001791657532056

Table A.5. Geometrical weights calculated for each station.

Station wi

A 0.9558689192

B 0.6198925398

C 0.0187554877

D 0.9149530337

A154, page 12 of 12