Georeferencing of an Unmanned Aerial System by Means of an Iterated Extended Kalman Filter Using a 3D City Model

In engineering geodesy, the technical progress leads to various kinds of multi-sensor systems (MSS) capturing the environment. Multi-sensor systems, especially those mounted on unmanned aerial vehicles, subsequently called unmanned aerial system (UAS), have emerged in the past decade. Georeferencing for MSS and UAS is an indispensable task to obtain further products of the data captured. Georeferencing comprises at least the determination of three translations and three rotations. The availability and accuracy of Global Navigation Satellite System (GNSS) receivers, inertial measurement units, or other sensors for georeferencing is not or not constantly given in urban scenarios. Therefore, we utilize UAS-based laser scanner measurements on building facades. The building latter are modeled as planes in a three-dimensional city model. We determine the trajectory of the UAS by combining the laser scanner measurements with the plane parameters. The resulting implicit measurement equations and nonlinear equality constraints are covered within an iterated extended Kalman filter (IEKF). We developed a software simulation for testing the IEKF using different scenarios to evaluate the functionality, performance, strengths, and remaining challenges of the IEKF implemented.


Introduction
It is possible nowadays using the increasing capability and ongoing miniaturization of hardware and sensors to construct various kinds of multi-sensor systems (MSS). The MSS are frequently used in engineering geodesy to capture the environment. It is crucial to know the position and orientation of the MSS with respect to a superordinate coordinate system 1 for a further processing of the data captured. The task to determine position and orientation in a superordinate coordinate system is called georeferencing. Georeferencing usually comprises the determination of three translations and three rotations, also known as six Degrees of Freedom (DoF) or pose.
Unmanned aerial vehicles (UAV) have become a widespread, useful and affordable platform for the MSS in the past decade, optimally suited to capture common scenarios in engineering geodesy. The MSS mounted on a UAV are subsequently called unmanned aerial systems (UAS).
Precise georeferencing of a UAS is a challenging task, especially in urban areas. Global Navigation Satellite System (GNSS) data are affected by shadowing and multipath effects and are, therefore, inaccurate or even unavailable. Basically, low-cost GNSS receivers obtain positional accuracies of a few meters. These accuracies may be improved to decimeter or even centimeter level under good GNSS conditions and using satellite-based augmentation systems, differential GNSS, or real-time kinematics. Lightweight and low-cost inertial measurement units (IMU) 2 obtain orientation accuracies of 0.1 • for roll and pitch and slightly worse accuracy, 0.8 • , for heading. Combined GNSS-IMU systems improve the heading accuracy to 0.2°-0.5 • . However, a general problem with IMU data is that it is often seriously affected by drifts.
A further possibility for georeferencing of an MSS is tracking using an external sensor. This is a cumbersome and inefficient task in the case of georeferencing a UAS due to limited measuring ranges and occultations.
In this paper, we present the georeferencing of a UAS, which is equipped with a three-dimensional (3D) laser scanner, a low-cost GNSS receiver, and an IMU, among other sensors. We set ourselves the requirement to determine the pose of the UAS with a precision (as a measure of the accuracy) better than 10 cm in position and 0.1 • for orientation. The accuracy of the pose is decisive for the further results derived, such as digital terrain models or detailed 3D city models. Therefore, we fused the UAS-based laser scanner measurements towards building facades, which are modeled in a 3D city model, and the other sensor data within an iterated extended Kalman filter (IEKF). We developed a software simulation for testing the IEKF developed using different scenarios to evaluate the functionality and performance of the IEKF implemented. All results presented are part of a research project, subsequently called the UAS project.

Georeferencing of MSS/UAS
The strategies to accomplish the georeferencing vary depending on the measurement configuration, the sensors and hardware available, the environment, and the accuracies to be fulfilled. The georeferencing strategies can be generally classified into direct, indirect, and data-driven georeferencing (see, e.g., Schuhmacher and Böhm 2005;Paffenholz 2012;Holst et al. 2015). An alternative classification regarding indoor applications is given in Vogel et al. (2016).
For direct 3 georeferencing, sensors which measure the pose of an MSS directly are integrated into the MSS. These sensors could be, for example, a GNSS receiver (Paffenholz 2012;Talaya et al. 2004), an IMU (Talaya et al. 2004), or an external sensor, such as a total station or a laser tracker (Dennig et al. 2017;Hartmann et al. 2018), which determine the pose by angle and distance measurements to a reflector on the MSS.
Indirect 4 georeferencing comprises methods where the pose of the MSS is determined by measurements towards known targets. These targets may be flat markers with a specific pattern (Abmayr et al. 2008) or simple 3D geometries, such as cylinders or spheres (Elkhrachy and Niemeier 2006). The position of the targets within the superordinate coordinate system is determined using an external sensor, for example, a total station. Indirect georeferencing is commonly used in bundle adjustment or for the georeferencing of static terrestrial laser scanners (TLS).
Data-driven georeferencing conforms basically to the indirect georeferencing. Instead of known targets, the data sets are matched to reference data sets. These reference data sets may be point clouds georeferenced already (Soloviev et al. 2007;Glira et al. 2015), digital surface models or 3D city models (Hebel et al. 2009;Li-Chee-Ming and Armenakis 2013;Unger et al. 2016Unger et al. , 2017. The matching is accomplished, for instance, via point-to-point assignment, for example, an iterative closest point algorithm (Besl and McKay 1992) or point-to-object assignment (see, e.g., Schuhmacher and Böhm 2005). Please note that data-driven georeferencing approaches vary widely. A huge part can be found in the commonly used Simultaneous Localization and Mapping (SLAM) approaches (see, e.g., Nguyen et al. 2006;Lee et al. 2007;Jutzi et al. 2013;Kaul et al. 2015;Nüchter et al. 2015).
In principle, all the georeferencing strategies mentioned are suitable for georeferencing static and kinematic MSS. For a static MSS, the six DoF have to be determined only once, whereas for a kinematic MSS, the six DoF for each measuring epoch should be determined. It may be indispensable to model the pose almost continuously depending on the measuring frequency of the sensors used.

Filtering Techniques for Georeferencing
Georeferencing of a kinematic MSS is commonly covered as part of the state parameter vector within a filtering approach. Filtering is a two-step procedure in which the current state parameter vector is estimated based on the previous state parameter vector and the current observations. In the first step, called prediction, the current state parameter vector is predicted using the previous state parameter vector and a specific system model. The system model takes into account all the controlling variables and environmental noises that affect the current state. In the second step, called the measurement update step, the predictions are modified based on the observation equations which compare the predicted observations with the current original ones.
In recent years, many pose estimation algorithms based on the extended Kalman filter (EKF) have been successfully applied to solve the pose estimation problem of a UAS. Tailanian et al. (2014) focus on the sensor fusion of the local sensors of a UAS in which the GNSS and the IMU have been combined by means of an EKF. In Hol et al. (2007), pose estimation on a six DOF robot using an EKF to fuse this information has been shown. Forster et al. (2013) use the collaboration of multiple UAS for pose estimation to combine multiple SLAM algorithms and create an accurate pose estimation. In their approach, real-time camera pose estimation is accomplished by combining the inertial and vision measurements using nonlinear state estimation approaches.
In the case of using nonlinear observation equations, a linearization (realized by means of Taylor series expansion) has to be applied to overcome the nonlinearity issue (see Denham and Pines 1966). In the case where the effects of the linearization errors tend to affect the efficiency of the filter or its convergence, the re-linearization of the measurement equation around the updated state may reduce these difficulties. Therefore, such a procedure is called IEKF. Researchers commonly use explicit observation equations in the IEKF. This means that the observations are considered as a function of the state parameters. Such an observation model is generally called a Gauss Markov model (GMM).
If the equations relating the observations to the state parameters are condition equations, for example, some 3D points should fulfill the plane equations, then we do not exhibit the typical formulation of the GMM. In other words, the observations and the state parameters are not separable (implicit measurement equations, see Dang 2007). In such a case, we are dealing with a Gauss Helmert model (GHM). In Dang (2007Dang ( , 2008, Steffen and Beder (2007), Steffen (2013) and Ettlinger et al. (2018), an IEKF is used which deals with implicit measurement equations. In Vogel et al. (2018), an IEKF is used for georeferencing and extended with additional nonlinear equality constraints. In Vogel et al. (2019), the approach is further extended by integrating nonlinear inequality constraints. With the ability to handle implicit measurement equations and nonlinear inequality constraints, it is possible to depict almost any mathematical relationship within an IEKF.

Contribution
In this paper, we present an adaption of the IEKF in Vogel et al. (2018) for data-driven georeferencing of a UAS with an accuracy better than 10 cm in position and 0.1 • in orientation. The highlight of the IEKF presented is the fusion of laser scanner measurements towards building facades, modeled as planes in a 3D city model, with the corresponding plane parameters of the 3D city model. This fusion leads to nonlinear implicit measurement equations and additional nonlinear equality constraints, which are covered within the IEKF. To evaluate the functionality and performance of the IEKF, we developed a simulation tool with a simple dynamic system model defining the state of the UAS and other states. On the basis of the states defined, we simulated sensor measurements regarding the sensor specifications and further characteristics. We tested the IEKF for different scenarios and analyzed its sensitivity towards different data characteristics.

Outline
The dedicated sections are as follows: In Sect. 2 we give a detailed description of the UAS, the IEKF implemented and the entire workflow for georeferencing. Section 3 displays the simulation tool developed and gives an overview of the specifically chosen values for the IEKF. We present and discuss the results of the IEKF applied to the simulation in Sect. 4. Finally, conclusions and an outlook are drawn in Sect. 5.

IEKF for Georeferencing of a UAS
The Geodetic Institute (GIH) and the Institute of Photogrammetry and GeoInformation (IPI) of Leibniz University Hannover (LUH) are currently working on the UAS project.
The UAS project deals with the precise determination of the trajectory of a UAS by integrating camera and laser scanner data in combination with generalized object information. Within this paper, we will focus on the usage of laser scanner data in combination with a non-generalized 3D city model. We will exclusively use simulated data to evaluate the functionality and performance of the algorithms developed. The simulation scenarios were chosen in accordance with common measurement scenarios and data characteristics.

General Idea
The UAS moves through an urban area where the building facades are modeled as planes within a 3D city model. Threedimensional city models with a level of detail 1 (block model) and a level of detail 2 (model with differentiated roof structures and boundary surfaces) are freely available for many cities. A detailed classification of the different levels of detail is given in Gröger et al. (2012). The pose of the UAS is roughly known, for example, from measurements of a GNSS receiver and an IMU. A 3D laser scanner captures the environment (red dots in Fig. 1) continuously. These captured 3D points may represent the ground, vegetation, building facades, or other objects. The laser scanner measurements are given in the local sensor coordinate system. The laser scanner measurements are transformed in a superordinate coordinate system according to the roughly known translation t = t x , t y , t z and orientation o = [ , , ] (see Fig. 2). The transformed laser scanner measurements are assigned to planes of the 3D city model based on the distance between the scanner points and the planes of the building model (see Sect. 2.3). Only the points that are close enough to a plane of the 3D city model (green points in Fig. 2) are used as observations afterwards. The final pose parameters are estimated within a GHM, resp. filtering approach, by minimizing the distance between the assigned laser scanner measurements and the planes of the 3D city model to which they are assigned.

UAS
A UAV is equipped with a 3D laser scanner, cameras, a GNSS receiver and an IMU in the UAS project. Figure 3 depicts a simplified sketch of the platform setup (without cameras). The laser scanner 5 scans 16 scan lines which are nearly perpendicular to the sensors' vertical axis (that corresponds to the z-axis in Fig. 3). The divergence between each scan line is 2 • . Thus, the laser scanner has a field of view of 30 • × 360 • . It is possible to set the resolution of the points in the scan lines between 0.1° and 0.4 • . The rotation rate depends on that setting and is between 5 Hz and 20 Hz. We set the resolution to 0.4 • to obtain a higher rotation rate of 20 Hz. Furthermore, it is possible to exclude certain angle areas from the measurement, because they cannot provide any data or the data generated are not needed for further processing. In our case, it would make sense to exclude the angle areas where the laser scanner measures towards the UAV.  ] . The origin of the laser scanner is depicted in the upper left corner. Note the slight difference to the original measurement configuration due to the roughly known pose. The transformed points which are close to one of the planes are assigned to that plane (green dots). The closeness is determined by a distance threshold (dashed blue lines). In this example, 5 points are assigned to plane 1 and 4 points are assigned to plane 2, whereas the other points (red dots) are not assigned. Subsequently, only the green dots are used as observations in the GHM, resp. filtering approach

Assignment Algorithm
A crucial task of our georeferencing process is the assignment of measured 3D points C loc , given in a local coordinate system of the laser scanner, to a plane of the city model, given in the global coordinate system. Therefore, we first have to transform C loc to C glo : The parameters of the translation t describe the position of the origin of the local coordinate system. The rotation matrix R is obtained based on the orientation o = [ , , ] of the local coordinate system according to Luhmann (2013): (1) C glo = t + R ⋅ C loc . ( The assignment is realized by a simple distance criterion, as described in Unger et al. (2016Unger et al. ( , 2017. The Euclidean distance to each plane of the city model is calculated for each 3D point in C glo . Whether the projection of the points into the plane lies within the bounding polygon of the planes is checked for points where the distance to the nearest plane is less than a threshold value d assign . If it is outside, the distance from the point to the boundary polygon of the plane is calculated and replaces the orthogonal distance. Points that project inside the boundary or are closer than d assign are assigned to their closest plane. d assign has to be selected regarding the accuracy of the points as well as the accuracy of the city model. The accuracy of the points depends mainly on the accuracy of the translation and orientation parameters t and o and the accuracy of the measured points C loc . The accuracy of the city model depends on the geometrical accuracy of its vertices and the extent of generalization effects by which the model deviates from the reality captured. The larger the threshold is fixed, the more points are assigned to planes, but the higher the probability is that incorrect assignments will be set that lead to incorrect results.

IEKF (Iterated Extended Kalman Filter)
The IEKF which we adapted for georeferencing of the UAS is given in Vogel et al. (2018). This approach can be applied for many different use cases. All available observations, types of prior information, and requested states need to be linked with the respective uncertainty information. Based on this, all information can be introduced into the IEKF. According to the respective application, several relationships have to be established by means of linear or rather nonlinear functional models. A system model f (⋅) is needed for the requested states x k to describe the physical behavior of the system from epoch k − 1 to the current epoch k . Realization of this  transfer is carried out within the prediction step of the IEKF. Suitable functional relationships for each observation l k have to be formulated within the measurement model h(⋅) . Based on implicit ( h(E(l), x) = ) or explicit ( E(l) − h(x) = ) formulations, relations between the observations available and requested states are considered. E(⋅) is the expected operator of a random vector; here, the expected value of the observation vector l , which can be replaced by E(l) = l + v . This consideration of current observations is carried out within the measurement update step. Available prior information can also be integrated within the measurement model during the measurement update step. In addition, further suitable prior information can also be applied in terms of state constraints by means of additional linear ( Dx k = d ) or nonlinear ( g x k = b ) functions. Here, D is a known matrix and d and b are known vectors with respect to related state constraints. Formulation of such restrictions can be implemented by means of equality constraints within the constraint step.
The basic workflow of this IEKF is depicted in Fig. 4 in a simplified way. Detailed equations for initialization, prediction, measurement update, and constraints are given in the following Sects. 2.4.1-2.4.6. Algorithm 1 given in Sect. 2.5 summarizes the required input, different computations steps and the output.

State Parameter Vector
The state parameter vector x k consists of two parts, as already depicted in Fig. 4. The first part x State,k ∈ ℝ 9 consists of the current position t k , orientation o k and velocity v k of the UAS (see Eq. 6). These state parameters describe the state of the UAS.
The second part consists, on one hand, of the vector x Plane,k ∈ ℝ 4⋅E , which contains the parameters of all the E city model planes in Hesse normal form (see Eq. 7). The Hesse normal form is defined by the 3 × 1 normal vector n e = [n x,e , n y,e , n z,e ] T and the distance to the origin d e with e = 1 … E . On the other hand, it consists of the vector x V,k ∈ ℝ 3⋅M , which contains all the M vertices of the building model (see Eq. 8): The integration of plane parameters and vertices into the state parameter vector is used to identify and correct planes that are not accurately represented in the 3D city model. Although this purpose is not part of this paper, we are already introducing the mathematical relationship.
The complete state parameter vector x k is arranged according to:

Observation Vector
The MSS mentioned in Sect. 2.2 provides discrete 3D laser scanner point clouds (LSC) and 6D pose information by means of a GNSS receiver and an IMU. An LSC consists of a full scan rotation ( 30 • × 360 • ). For the sake of simplification, we will not consider the movement of the UAV during the period of a full scan rotation in the following. This time period is a maximum of 0.05 s for the configuration chosen. It is possible to exclude certain angle areas of the laser scanner, as has already been described in Sect. 2.2. The time period for the non-excluded areas is further reduced by excluding the angular areas in which the laser scanner measures in the direction of the UAV. This is acceptable due to the planned velocity of the UAV of about 1 m s −1 and an angular velocity of 2 • s −1 . In further development steps, the LSC should consist of fewer and temporally closer 3D points, for example, just a half or a quarter scan rotation. In addition, the initial vertices of the 3D city model are introduced as observations into our approach. The observation vector l k consists of three parts for each epoch k = 1 … K . The first part l loc Scan,k consists of the measured LSC in the local sensor coordinate system with a total of N 3D points. l loc Scan,k only contains the 3D points, which were assigned to a plane of the city model 6 (for the assignment algorithm, see Sect. 2.3). These points are assorted on the basis of the planes they are assigned to: l loc Scan,e,k is representing the vector with all N e points assigned to plane e. The total number N of 3D points of the LSC stored in the observation vector is calculated by: The second part l glo Pose,k consists of the direct GNSS and IMU observations for the position and orientation: Finally, the third part l glo V,0 consists of the, in total M, initial vertices of the E model planes 7 : with representing a 3D point which is the vertex of at least one plane.
The observation vector is arranged as follows: We apply Eq. (19) for building the stochastic model of the observation vector. We neglect correlations in all cases for a better discussion of the mainly important issues: The standard deviation of a laser scanner coordinate component is denoted by denote the standard deviation of the GNSS and IMU observations, and V denotes the standard deviation of the initial vertices.

System Equation
Formulating the system equation of the kind: we neglected the control vector u by setting it to zero. The system noise is normally distributed with w k−1 ∼ N(0, ww ) . Equation (20) can be formulated as: where F x,k denotes the transition matrix: Due to the fact that we do not intend to develop an optimal system model within this paper, we applied rather simple linear system equations. Nevertheless, a subsequent adaption of the system model with more complex system equations is easily possible if the future flight characteristics of the UAS require it. For the system equations chosen, the transition matrix F x,k is given by: with Δ is the time period between two consecutive epochs. We apply the following equations for building the variance-covariance matrix (VCM) of the system noise ww : The standard deviations of the system noise for translation, orientation, and velocity are denoted by The system noise for the plane parameters and vertices is given by the standard deviations n,w = n,w , n,w , n,w , d,w and V,w . (20)

Measurement Equation
With the state parameter and observation vector given in Sects.
The rotation matrix R k is calculated based on the current orientation o k according to Eqs. (2)-(5). Subsequently, we obtain N measurement equations of type h I : Figure 5 shows an illustration of the relationship between state parameters, observations, and measurement equations.
In the second type of measurement equation h II , we calculate the difference between the estimated pose t k ;o k and observed pose l glo Pose,k of GNSS and IMU: If the GNSS and IMU data are available, we obtain six measurement equations of type h II . When the GNSS signal is missing, the measurement equation changes to: and we only obtain three measurement equations of this type.
The third type of measurement equation h III is introduced to avoid a datum defect. It computes the difference between the estimated vertices in x V,k and the initial vertices in l glo V,0 . If the vertex V glo m,k is stored in x V,k and V glo m,0 is its initial observation stored in l glo V,0 , the following measurement equations are applied: For a three-dimensional vertex Eq. (33) consists of three measurement equations. Subsequently, we obtain 3 ⋅ M measurement equations of type h III : Altogether, the measurement equations are given by:

Nonlinear Equality Constraint for the State Parameters
We apply two types of nonlinear equality constraints for the state parameters in our IEKF. The prior information which we want to model describes hard constraints which have to be fulfilled. The first type of nonlinear equality constraint g I arises due to the fact that we are using plane parameters within the state parameter vector x k . Here, we must ensure unit normal vectors by means of a length of one. For this, we can make use of the nonlinear equality constraints: with e ∈ {1 … E}.
The nonlinear equality constraint of type g I applies to each of the E planes normal vectors: The right side of the equal sign is stored in vector b I : With the second type of nonlinear equality constraint g II , we ensure that each vertex of a plane is located in the plane. By this, we preserve the topology of the 3D city model. If we assume that V glo m,k is a vertex of plane e, represented by n e and d e , the following equality constraint is applied: This type of nonlinear equality constraint must be fulfilled for each vertex of all planes. In general, 3D city model planes have 4 vertices, but it is also possible for them to have 3 or more than 4 vertices. Thus, the number 4 ⋅ E is just a rough estimate for the total amount of equality constraints of type g II : The right side of the equal sign is stored in vector b II : A comparable procedure can be found in Unger et al. (2016). Altogether, the nonlinear equality constraints for the state parameters are given by:

Initialization
The initial state parameter vector x 0 is created by the GNSS and IMU observations. We assume zero for the velocity in each coordinate component: The initial plane parameters n e,0 and d e,0 were estimated from the planes' vertices stored in the city model by the Drixler algorithm (Drixler 1993). We make use of the vertices of the 3D city model given for the initial vertices x V,0 : We apply the following equations for building the VCM of the initial state parameter vector xx,0 : The standard deviations of the initial state parameters of translation, orientation, and velocity are denoted by The standard deviations of the initial plane parameters and vertices are denoted by n,0 = n,0 , n,0 , n,0 , d,0 and V,0 .

Workflow
The workflow of our algorithm is summarized in Algorithm 1. This is an adaption of the algorithm proposed in Vogel et al. (2018). We highlighted the cross references to equations and sections in blue. The partial derivatives (see lines 9, 22, 23, and 40) were obtained once using the symbolic Math toolbox of Matlab © . Subsequently, we implemented the partial derivatives in a function. Another possibility to obtain the partial derivatives could be the use of INTLAB (Rump 1999). We fixed the stop criterion c stop to 1 × 10 −12 in the measurement update step. This stop criterion was reached after an average of six iterations in the subsequent simulation.

Simulation
Computer simulations are a great tool for analyzing and interpreting engineering systems. Here, we intended to evaluate the functionality and performance of the IEKF implemented. Therefore, we focused on two scenarios which basically conform to the testing environment designated for the UAS Project. Furthermore, we focused on the challenging initialization phase of the IEKF and analyzed a rather short trajectory of 2.5 m length and with 50 epochs. Therefore, we created a model of a fictitious building and a ground plane. Subsequently, we determined a fictitious trajectory beside this building. We assumed a constant velocity for the UAS of 1 m s −1 and simulated laser scanner measurements with a frequency of 20 Hz. Consequently, the epochs simulated are at constant distances of 0.05 m. Each of the epochs consists of a laser scan and the desired values of position and orientation obtained from the determined trajectory at a certain time.
Figures 6, 7 and 8 depict different views of the building and ground plane created. The building created contains roof planes (in dark grey), wall planes (in red), and window planes (in transparent blue). Figure 9 depicts the corresponding 3D city model of the building. As we can see, the ground plane is not part of the 3D city model. The 3D city model contains the vertices of each plane and is given in the superordinate coordinate system. The simulated and correctly georeferenced point cloud (blue dots) of the first epoch and the trajectory of the UAS (magenta line) is depicted in each figure.

Scenarios
We simulated two scenarios in which the laser scanner and IMU measurements were generated differently under certain assumptions regarding measurement accuracy and bias. We repeated the simulation 500 times ( S = 500 ) for both scenarios. In scenario 1, we only added normally distributed noise on the laser scanner measurements, the position (representing the GNSS receiver), and the orientation (representing the IMU). In scenario 2, we systematically perturbed the laser scanner measurements, hitting the windows of the simulated building, and increased the (standard deviation of the) noise for laser scanner measurements, hitting the windows and the ground. The systematic and increased random disturbances of the measurements hitting the windows simulate the infiltration behavior for glass. We included the increased measurement uncertainty and the actual infiltration of the laser into the increased noise of these measurements. The increased noise for measurements hitting the ground simulates possible unevenness and vegetation. Thus, this increased noise is more justified by actual structures in the object space than by an increased uncertainty of the measurement. Furthermore, we added a drift on the -component of the IMU. Therefore, we linked the systematic effect Δ to the epoch number k ∈ {1, … , K} . The -component conforms to the heading of the UAS. All assumed systematic effects Δ and standard deviations of the added noise are depicted in Table 1. Basically, the assumptions of scenario 1 should be consistent with the manufacturer's specifications of the sensors used or planned for our UAS. The assumptions in scenario 2, especially the ones for the laser scanner, are based on experience from test measurements. In our simulation, we assume that the positions and orientations of all sensors in a platform coordinate system or body frame have already been determined in a calibration process. We also assume that the sensors are properly synchronized. To neglect the effect of generalization in the 3D city model, we used the same model for simulation and in the following IEKF algorithm.

Values Chosen for the IEKF
We used a constant distance threshold d assign = 0.3 m for the assignment algorithm (see Sect. 2.3) in our simulation. This chosen threshold is a trade-off between the rather imprecise pose in the first epochs and the more precise pose in the last epochs.
The standard deviations used for building the initial VCM xx,0 (see  are depicted at the top of Table 2. The standard deviations of the initial position t,0 and orientation o,0 conform to the standard deviations of the GNSS and IMU observations. We assume a standard deviation for the initial velocity v,0 which includes the planned maximal velocities of the UAV. For the standard deviations of the initial plane parameters n,0 and d,0 , we met the assumption that these values should be very small. Thus, they correspond roughly to the standard deviation of the initial vertices V,0 , which, in turn, conform to the standard deviation chosen for the observation vector. If the accuracies of the planes' vertices are known, the standard deviations of the plane parameters could also be determined by means of variance propagation. The standard deviations chosen for the system noise (see Eqs. 25-28) can be found in the middle of Table 2. We set the standard deviation for translation t,w , orientation o,w and velocity v,w depending on the time period Δ between two epochs and regarding possible unpredictable movements of the UAS in this time period. We set the system noise to zero for ww,Plane and ww,V , because these parameters are constant over time. Consequently, n,w , d,w , and V,w are zero.
We used the simulated standard deviation LS = 0.02 m of the laser scanner measurements for the stochastic model of the observation vector (see Eq. 19). We also introduced the simulated standard deviations t = 0.5 m and o = 0.2 • for the GNSS and IMU observations. In the case of the initial vertices of the 3D building model, we introduced a small standard deviation of V = 0.1 mm . The 3D city model is, thus, almost introduced as a fixed feature. We are aware that this assumption is too optimistic and that the accuracy of the introduced 3D city model is rather in the range of a few centimeters. However, in the context of this paper, the influence of a less precise 3D city model will not be considered in more detail.
We assume that ll,k and ww do not contain correlations on their own and among themselves.

Results and Discussion
In the following, we depict the results of the IEKF for both scenarios and compare them with the results of a linear Kalman filter (LKF), which only uses the simulated GNSS and IMU observations. The comparison between LKF and IEKF evaluates the benefit of introducing laser scanner observations and a city model. The state parameter vector x k of the LKF solely consists of x State,k (Eq. 6). The observation vector l k consists of l glo Pose,k (Eq. 15). The transition matrix F x,k of the LKF is given by F x,State,k (Eq. 24). Consequently, the VCM of the observation vector and of the system noise have to be adapted. The LKFs' only remaining measurement equation (Eq. 31) has to be transformed in explicit form. For the LKF, no linearization is necessary and no constraints are used. For more information concerning the LKF, see (Simon 2010). Figures 10, 11, 12, 13, 14 and 15 depict the median values (solid lines) and 68%(1 ) confidence intervals (CI) (dashed lines) for the translation parameters t x , t y , t z (Figs. 10, 11, 12) and orientation parameters , , (Figs. 13, 14, 15) in scenario 1. The CI are calculated numerically using the total number of simulated runs S according to Alkhatib et al. (2009). The red lines result from the filtered state parameters x + k of the IEKF after the constraint step. The blue lines result from the filtered state parameters x + k of the LKF, respectively. The true values are plotted in a dashed black line.
The median of the IEKF result for t x (across flight direction) in Fig. 10 is very close to the true translation, whereas we can see slight random deviations between the median of the LKF results and the true values. The latter can be seen for each pose component and is not explicitly mentioned again in the following. The CI of the IEKF and the LKF shrink rapidly in the first 5 epochs. Whereas the CI of the IEKF is very close to the median, the CI of the LKF is significantly broader. The latter can be seen for each pose component and is also not explicitly mentioned again. The IEKF performs well for that pose component because of the measurement constellation. Each measurement distinctively assigned to a plane provides information for that pose component. In other words, each measured plane of the 3D city model is sensitive for that pose component.
The results are similar for t y (in flight direction) in Fig. 11. Again, the median of the IEKF results is very close to the true translation and its CI shrinks within the first epochs. For this pose component, the shrinking of the CI (IEKF) is slightly slower than for t x and lasts until epoch 10. A reason for that slight difference can be found in the weaker measurement constellation for that pose component. Only the triangle-shaped protrusions are sensitive for t y . Consequently, there are significantly fewer observations which provide information for that pose component. Figure 12 shows the results for t z (up-direction). The median of the IEKF results is very constant, but seems to deviate from the true values by a small systematic shift of about 3 cm. This systematic shift is caused by ground points wrongly assigned to a wall plane of the building. This results in a slight rotation around the y-axis (see results for ) in combination with a slight descent of t z . Again, the CI of the IEKF translation shrinks rapidly within the first 6 epochs. The boundaries of the CI in the first epochs are not symmetrical around the median. They are slightly shifted to larger t z values. That means that the IEKF converges slower towards the true translation with initial t z values larger than the true translation. An initial t z value larger than the true translation may cause significantly fewer assignments in the first epochs. Only the roof planes are sensitive for that pose component. The results for , which represents the rotation around the x-axis, are shown in Fig. 13. After 10 epochs, the median of IEKF and LKF obtain comparable results, whereas the median of the IEKF seems to vary slightly more. The CI of the IEKF shrinks within the first 12 epochs. Similar to the t z pose component, there is a shift of the IEKF CI towards larger values. The only sensitive planes for that pose component are the triangle-shaped protrusions and the roof planes. This may explain the larger variations in the median of the IEKF orientation. Figures 14, 15 show the results for (rotation around y-axis) and (rotation around z-axis), respectively. We can see similar characteristics for these pose components. In both cases, the median of the IEKF results is constant but systematically shifted by a small value of about 0.05 • . As has already been mentioned for t z , the systematic shift in is caused by wrongly assigned ground points. Because of the skewed alignment of the laser scanner, the wrongly assigned ground points also cause a slight rotation in . All planes are sensitive for both pose components.
For scenario 2, we depicted the results in a similar way in Figs. 16,17,18,19,20 and 21. In Fig. 16, we can see the results for t x . Similar to scenario 1, the median of the IEKF results comes very close to the true translation. The lower boundary of the CI (IEKF) is farther afield from the median than the upper boundary. It shrinks until epoch 14. An explanation for the one-sided enlarged CI (IEKF) is the systematically extended laser scanner measurements hitting the windows. With an initial t x translation, which is farther afield from the building than the true translation, these systematically extended laser scanner measurements are wrongly assigned to the walls of the building. After some epochs, these false assignments seem to have disappeared. These effects are analyzed later on (see Fig. 29).
The median of the IEKF results for t y in Fig. 17 is close to the true translation. The upper boundary of the CI (IEKF) shrinks until epoch 13, whereas the lower boundary of the CI (IEKF) clearly departs from the median from epoch 14 on. Again, these effects are analyzed later on (see Figs. 28 and 29).
We can see the results for t z in scenario 2 in Fig. 18. The results are comparable to the ones obtained in scenario 1 for that pose component. A significant difference can be found in the upper boundary of the CI (IEKF), which converges to the median much slower than in scenario 1. In scenario 2, it does not converge until epoch 15. Again, initial t z values larger than the true position cause significantly fewer assignments in the first epochs.
We observe comparable results to scenario 1 for the pose parameters and in Figs. 19 and 20. In both figures, we can observe that one boundary of the CI (IEKF) is farther afield from the median than the other and that the CI (IEKF) shrinks slower than in scenario 1.
The pose component is depicted in Fig. 21. The CI (IEKF) shrinks until epoch 12. In addition, we can see that the simulated drift in does not affect the IEKF results, whereas the median of the LKF departs from the true orientation continuously as simulated.
A suitable value to evaluate the performances of the filters is the root-mean-square error (RMSE). In our simulation, the RMSE is the error between filtered state parameters and the true state parameters. The RMSE of the pose parameter t x for run s ∈ {1, … , S} is calculated according to: t x,s,k results from the filtered state parameter vector ̂ + k from the IEKF or the LKF, respectively. t x,k are the true values.  The calculation of the RMSE for the other pose parameters is carried out analogously. Table 3 depicts the characteristic criterions: Minimum, maximum, mean, median, standard deviation (SD), and the lower bound ( ↓ ) and upper bound ( ↑ ) of the 68% and 95% CI of different RMSE. In the last row of Table 3, we display the rate of runs where the IEKF obtained a smaller RMSE than LKF and vice versa. Each criterion is divided into two rows regarding scenario 1 (above) and scenario 2 (below). Table 3 clearly shows that the IEKF obtains better results than the LKF. The median values especially are significantly smaller in each state parameter and each scenario. By contrast, the range of the IEKF results is significantly larger than the range of the LKF results, which can be seen in the SD and 95% CI values. A possible explanation can be found in a larger number of runs, where the IEKF is far away from the true values. We will subsequently call these runs failures. These failures can be caused by a disadvantageous initial pose parameter in the first epoch. By this, a large number of laser scanner measurements are assigned to wrong planes and the IEKF converges to a wrong pose. The interaction between initial pose parameters and performance of the IEKF is analyzed at the end of this section. The failures are also the explanation for the large discrepancy between the mean and median values.     To evaluate the performance of the IEKF over the epochs, we calculated a mean RMSE for each epoch over all S simulation runs. The calculation is as follows for the RMSE of t x in epoch k ∈ {1, … , K}: We obtain the RMSE for the other pose parameters analogously. The results can be found in Figs. 22,23,24,25,26 and 27. Figures 22,23,24,25,26 and 27 generally support the interpretations made before. The RMSE of the translation parameters t x and t z determined by the IEKF decrease with an increasing number of epochs. At the last epoch, the RMSE t x ,Ep. and RMSE t z ,Ep. of the IEKF is significantly smaller than the RMSE obtained by LKF in both scenarios. Again, RMSE t y ,Ep. differs and increases after epoch 10 with an increasing number of epochs. In both scenarios, especially in scenario 2, the RMSE t y ,Ep. obtained by the IEKF is larger than the one obtained by LKF in the last epochs. The RMSE of the orientation parameter show different characteristics, but they are intrinsically quite similar. The RMSE obtained by IEKF starts with large values and significantly decreases after epoch 10. In the last epoch, RMSE ,Ep. and RMSE ,Ep. of IEKF and LKF are quite similar. RMSE ,Ep. of the IEKF undercuts the RMSE ,Ep. of the LKF until epoch 9 (scenario 1), resp. epoch 30 (scenario 2). The RMSE ,Ep. of LKF in scenario 2 is clearly affected by the simulated IMU drift. It is remarkable that (except in RMSE t y ,Ep. ) the RMSE of the IEKF between scenario 1 and scenario 2 is very similar in the last epoch.
For a deeper understanding of the link between initialization and performance of the IEKF, we plotted the initial values of GNSS and IMU of the first epoch in parallel coordinates. The parallel coordinate representation was first used by Inselberg (1985). For this representation, we will show all pose components in six parallel axes. On each axis, all initial values belonging to this component are depicted. This representation is used to identify whether some extreme initial pose values could have had a large impact on the divergence of the IEKF. In addition, we colored the runs where the IEKF failed (failure) in orange and the others (success) in blue. The distinction in failure and success was selected based on the RMSE of the position in the last epoch. The run is classified as a failure when there is a value larger than 10 cm.
In scenario 1 (Fig. 28), the failure rate is 7.6%. For the initial pose parameters t z , , , and , the failures are distributed over the whole value range. For t x , the failures are slightly concentrated in a range between − 9.2 and − 8.55 m and a second range around − 10.6 m. For t y , the failures are even more concentrated in a range between 3.72 and 4.3 m.
Initial values for t y in a range between 3.72 and 4.3 m lead to most of the failures. As previously mentioned, we have a weak measurement constellation for t y . Furthermore, the trajectory is parallel to the y-axis with a true value of t y = 7.5 m in epoch 50, the consequence being that with initial values for t y between the true value (5 m) and 6.19 m, the IEKF may get back to the true values in a later epoch. This possibility is not given for initial values for t y between 3.72 and 4.3 m.
In scenario 2, the failure rate of 20.4% is significantly higher. Again, for the initial pose parameters t z , , , and , the failures are distributed over the whole value range. For t y , the failures are mainly in a range between 3.72 and 5 m. The reason for that has already been explained for scenario 1. For t x , we can see a strong concentration of failures in a range between − 11.37 and − 10.4 m. Thus, especially that range for t x causes the failures in scenario 2. The reason for that can be found in the systematically perturbed laser scanner measurements hitting the windows. These measurements are extended by an additive value Δ = 0.6 m . In combination with a disadvantageous initialization, these perturbed measurements are wrongly assigned to the wall plane and, therefore, cause the failures.
We used GNSS and IMU observations in the measurement update step for all the displayed results in this section. We performed the same experiment without using the GNSS and IMU observations in the measurement update step. The differences in the results are insignificantly small. Therefore, we do not show these results here.

Conclusions and Outlook
In this paper, we presented the mathematical basics and the workflow of an IEKF which makes it possible to determine the trajectory of a UAS better than 5 cm in position (median value) and 0.08 • for orientation (median value). It is not mandatory to have continuous GNSS and IMU observations for the implemented IEKF. The trajectory is mainly obtained using laser scanner measurements of building facades, which are modeled as planes in a 3D city model. The laser scanner measurements and the planes of the 3D city model are combined by implicit measurement equations and nonlinear equality constraints within the IEKF.
To demonstrate the functionality and performance of the IEKF implemented, we developed a simulator, which showed that the IEKF is even suited to handle systematically perturbed observations. Nevertheless, the algorithm demonstrated may be tuned to deal with disadvantageous initial values. A possible starting point can be found in the threshold of the assignment algorithm. Instead of a constant value, this value should be varied regarding the estimated accuracy of the predicted position.
In the future, we plan to evaluate the performance of the IEKF based on real data, especially regarding effects of generalization of the 3D city model and effects of a 3D city model with larger uncertainties. Thus, we can also check whether the assumptions made for the simulation (see Table 1) are also valid in reality. In particular, the performance of the GNSS receiver can deviate strongly from the simulated data with the same precision.
In the IEKF presented, the plane parameters and vertices are introduced as parameters to preserve the topology of the 3D city model. With an increasing number of planes and vertices, the computation becomes inefficient. Here, we plan to reorganize the IEKF presented by means of a dual-state Kalman filter.
As has already been mentioned, we plan to integrate camera measurements into the IEKF. This enables the stabilization of the trajectory in the long term and the improvement of the IEKF's initial behavior. Therefore, we have to integrate the collinearity equations into our approach. Subsequently, we have to integrate the object coordinates of the tie points into the state parameter vector. This emphasizes the need to use a dual-state Kalman filter. By integrating the camera measurements 8 into the IEKF, we must deal with observations captured in different epochs: one or more images observing a tie point were taken in a past epoch and only one image is taken in the current epoch.
Finally, we aim to analyze the benefit of integrating further geometries in addition to the planes of the 3D city model. There are many cylindrical-shaped objects in urban areas, such as lantern, traffic lights or street signs, which may be used in the IEKF.